% window_toad
% Dirk Husmeier
% August 2003
% This program moves a window of specified length
% and with a specified step size  along the 
% sequence alignment. For each window position the 
% program computes and plots the
% - average proportion of topology-defining sites,
% - average proportion of sites supporting topology 1,
% - average proportion of sites supporting topology 2,
% - average proportion of sites supporting topology 3.

clf
clear

load topos.in  % Vector of topology labels
load topoDefSites.in  % Vector of topo def sites

Nmax= input('Length of the alignment?');
WIN= input('Window size?')
STEP= input('Step size?')
Nsteps= floor((Nmax-WIN)/STEP)-1; % Number of different window positions

for i=1:Nsteps
    t_left=(i-1)*STEP+1;
    t_right= t_left-1+WIN;
    t_position(i)=t_left-1+WIN/2;
    % Binary indicator vector of topo def sites in window
    TopoDefSitesInWin=(topoDefSites>=t_left).*(topoDefSites<t_right);
    % Number and proportion of topo def sits in window
    N_topoDefSites=sum(TopoDefSitesInWin);
    % N_topoDefSites= ...
    % (topoDefSites>=t_left)*(topoDefSites<t_right)';
    P_topoDefSites(i)= 100*N_topoDefSites/WIN;
    % Number and proportion of individual topos
    N1= (topos==1)*TopoDefSitesInWin';
    N2= (topos==2)*TopoDefSitesInWin';
    N3= (topos==3)*TopoDefSitesInWin';
    P1(i)= 100*N1/N_topoDefSites;
    P2(i)= 100*N2/N_topoDefSites;
    P3(i)= 100*N3/N_topoDefSites;
end

xmin=0;
xmax=t_position(Nsteps);
ymin=0.9*min(P_topoDefSites);
ymax=1.1*max(P_topoDefSites);
subplot(2,2,1)
   plot(t_position,P_topoDefSites)
   axis([xmin xmax ymin ymax]);
   set(gca,'Fontsize',11)
   title('Topology-defining sites')
ymin=0;
ymax=70;
subplot(2,2,2)
   plot(t_position,P1)
   hold on 
   plot([1 t_position(Nsteps)],[100/3 100/3],'r--')
   axis([xmin xmax ymin ymax]);
   set(gca,'Fontsize',11)
   title('Topology 1')
subplot(2,2,3)
   plot(t_position,P2)
   hold on 
   plot([1 t_position(Nsteps)],[100/3 100/3],'r--')
   axis([xmin xmax ymin ymax]);
   set(gca,'Fontsize',11)
   title('Topology 2')
subplot(2,2,4)
   plot(t_position,P3)
   hold on 
   plot([1 t_position(Nsteps)],[100/3 100/3],'r--')
   axis([xmin xmax ymin ymax]);
   set(gca,'Fontsize',11)
   title('Topology 3')