disp('Running MATLAB script file UNC2000LinkScaSpa1.m') ;
%
%    Scale based Stationarity analysis of UNC 2000 Link data
%    This uses the versions binned to 2^20 bins


%addpath \matlab\stevesmo ;
%addpath \matlab\stevefun ;


ipart = 1 ;    %  1 - Make several images (for different gaps)

%visimul = [0 1 2 3] ;
visimul = [4 1] ;
                %  0 - no multiple comparison adjustment
                %  1 - adjust for multiple comparisons
                %  2 - Bonferroni adjustment
                %  3 - mixture
                %  4 - global simultaneous

                
vidat = [1] ;  
              %  1 - Peak
              %  2 - Off Peak


vitest = [3 4 5 2 1] ;
               %  1 - means and t-test
               %  2 - Kolomogorov-Smirnov
               %  3 - variances and F-test
               %  4 - mean removed Kolomogorov-Smirnov
               %  5 - standardized Kolomogorov-Smirnov



%  Set common parameters
%  
alpha = 0.05 ;
nminew = 8 ;
          %  minimum number in each window
          

%  Set Colormap
%
comap = [.35, .35, .35; ...
         0,    0,   1; ...
         1,    0,   0] ; 
colormap(comap) ;

figure(1) ;
clf ;



%  loop through data sets
%
for idat = vidat ;


  mfilepre = '\matlabR12\work\NetWorkData\UNCLink2000\MatlabData\' ;

  if idat == 1 ;    %  Peak Data
    mfile = ['UNC2000LinkPeakIP' 'bin'] ;
    dattit = 'Peak IP' ;
  elseif idat == 2 ;    %  Off Peak data
    mfile = ['UNC2000LinkOffIP' 'bin'] ;
    dattit = 'Off Peak IP' ;
	end ;    %  of idat if-block
	
  mfilestr = [mfilepre mfile] ;
	
  eval(['load ' mfilestr]) ;


  
  %  set parameters
  %
	n = nbin ;
	%  left comes from load (secs)
	%  right comes from load (secs)
	binwid = (right - left) / nbin ;
	vgap = [binwid 0.1 50] ;
	    %  corresponds to time ~ 0.02 secs (hopefully more than round trip time)
	%vgap = [1] ;
	%nwid = 10 
	nwid = log2(nbin) - 1 ;
                  %  number of widths to show at each scale
	%nloc = 400 ;
	nloc = 256 ;
                  %  number of locations to consider, for each scale
  tminew = nminew * binwid ;
                  %  minimum time span for each window
	tgrid = linspace(left,right,n)' ;




	%  loop through simul levels
	%
	for isimul = visimul ;
	
	if isimul == 0 ;   
            simultit = 'not simultaneous' ;      
	elseif isimul == 1 ; 
            simultit = 'simultaneous' ;      
	elseif isimul == 2 ; 
            simultit = 'Bonferroni' ;    
	elseif isimul == 3 ; 
            simultit = 'mixture' ;  
	elseif isimul == 4 ; 
            simultit = 'global' ;  
            mglobal = 2 * nloc - 1 ;
            salphaglobal = 1 - (1 - alpha)^(1/mglobal) ; ;
	end;
	
	
	
	%  loop through tests
	%
	for itest = vitest ;
	
	
	
		if itest == 1 ;    %  1 - means and t-test
            testtit = 'Mean t-test' ;      
		elseif itest == 2 ;    %  2 - Kolomogorov-Smirnov
            testtit = 'Kolm-Smir' ;      
		elseif itest == 3 ;    %  3 - variances and F-test
            testtit = 'Var. F-test' ;      
		elseif itest == 4 ;    %  4 - mean removed Kolomogorov-Smirnov
            testtit = 'Mean R. K-S' ;      
		elseif itest == 5 ;    %  5 - standardized Kolomogorov-Smirnov
            testtit = 'Std''zed K-S' ;      
		end ;    %  end of itest if-block
		
	
		
		
			
			if ipart == 1 ;
			
        %  loop through gaps
        %
        for igap = 1:length(vgap) ;

          gap = vgap(igap) ;
          go2 = gap / 2 ;
                %  gap over 2


          %  loop through wid's
          %
          mout = zeros(nwid,nloc) ;
          
          minwid = log10(tminew) ;
          maxwid = log10((right - left - gap) / 2) ;
          vwid = logspace(minwid,maxwid,nwid+1)' ;
          vwid = vwid(1:nwid) ;
                %  cut off last one, since doesn't add to plot
          maxwid = log10(vwid(nwid)) ;
                %  readjust maxwid
          
          minloc = (left + tminew + go2) ;
          maxloc = (right - tminew - go2) ;
          vloc = linspace(minloc,maxloc,nloc) ;
          for iwid = 1:nwid ;

            wid = vwid(iwid) ;

            disp(['  Working on wid ' num2str(iwid) ' = ' ...
                             num2str(wid) ...
                             '    for ' simultit ', ' testtit ', ' ...
                             dattit ', gap ' num2str(gap)]) ;

            locflag = left + wid + go2 <= vloc ;
            locflag = locflag & (vloc <= right - wid - go2) ;
            
            vindloc = (1:nloc) ;
            vindloc = vindloc(locflag) ;
            

            if isimul == 0 ;    %  don't bother about simultaneous inference
              salpha = alpha ;
            elseif isimul == 1 ;
              %m = (right - left) / (2 * wid) ;
              m = (right - left) / wid ;
              m = min(length(vindloc),m+1) ;
              salpha = 1 - (1 - alpha)^(1/m) ;

            elseif isimul == 2 ;
              salpha = alpha / length(vindloc) ;
     
            elseif isimul == 3 ;
              m = (right - left) / (2*wid) ;
              m = min(length(vindloc),m) ;
              salpha = 1 - (1 - (alpha/nwid))^(1/m) ;

            elseif isimul == 4 ;
              salpha = salphaglobal ;

            end ;

            
            
            for iloc = vindloc ;

              loc = vloc(iloc) ;
                
              flagleft = (loc - go2 - wid <= tgrid) ;
              flagleft = flagleft & (tgrid <= loc - go2) ;
                  %  ones where data inside left group

              flagright = (loc + go2 <= tgrid) ;
              flagright = flagright & (tgrid <= loc + go2 + wid) ;
                  %  ones where data inside right group

              dataleft = bincts(flagleft) ;
              dataright = bincts(flagright) ;


              
              if itest == 1 ;    %  means and t-test

                meanleft = mean(dataleft) ;
                meanright = mean(dataright) ;
                varleft = var(dataleft) ;
                varright = var(dataright) ;
                
                if meanleft == meanright ;                      

                  h = 0 ;
                      %  no significant difference

                else ;
                  
                  if varright == 0 ;
                    if varleft == 0 ;
                      h = 1 ;
                         %  significant difference

                    else ;    %  do 1 sample t-test
                      h = ttest(dataleft,meanright,salpha) ;
                         %  1 for significant difference found
                         %  0 otherwise
                    end ;

                  else ;

                    if varleft == 0 ;    %  do 1 sample t-test
                      h = ttest(dataright,meanleft,salpha) ;
                         %  1 for significant difference found
                         %  0 otherwise

                    else ;    %  do full two sample t-test

                    dfleft = length(dataleft) - 1;
										dfright = length(dataright) - 1;
										dfe  = dfleft + dfright;
										msleft = dfleft * varleft ;
										msright = dfright * varright ;
										
										difference = meanleft - meanright;
										pooleds    = sqrt((msleft + msright) * (1/(dfleft + 1) + 1/(dfright + 1)) / dfe);
										
										ratio = difference / pooleds ;
										p = 1 - tcdf(ratio,dfe) ;
										    %  one tail p-value
										p = 2 * min(p, 1-p) ;
										    %  two tail p-value

 								    h = 0;
	  								if p <= salpha, 
                        h = 1; 

                        end ;
                          %  1 for significant difference found
                          %  0 otherwise
                    
                    end ;
                
                  end ;
                    
                end ;
                  
              
                    
              elseif itest == 2 ;    %  Kolomogorov-Smirnov

                h = kstest2(dataright,dataleft,salpha) ;
                    %  1 for significant difference found
                    %  0 otherwise


              elseif itest == 3 ;    %  variances and F-test
              
                varleft = var(dataleft) ;
                varright = var(dataright) ;
                
                if varright == 0 ;
                  if varleft == 0 ;
                    h = 0 ;
                       %  no significant difference

                  else ;    %  do 1 sample t-test
                    h = 1 ;
                       %  significant difference

                  end ;

                else ;

                  if varleft == 0 ;    %  do 1 sample t-test
                    h = 1 ;
                       %  significant difference

                  else ;    %  do variance F-test

                    dfleft = length(dataleft) - 1;
										dfright = length(dataright) - 1;
	
										Fstat = varleft / varright ;
	
										p = 1 - fcdf(Fstat,dfleft,dfright) ;
										    %  one tail p-value
										p = 2 * min(p, 1-p) ;
										    %  two tail p-value

								    h = 0;
    								if p <= salpha, 
                      h = 1; 
                    end ;
                            %  1 for significant difference found
                            %  0 otherwise
                  
                  end ;
              
                end ;
                  
             
              elseif itest == 4 ;    %  mean removed Kolomogorov-Smirnov

                dataleftmr = dataleft - mean(dataleft) ;
                datarightmr = dataright - mean(dataright) ;
                
                h = kstest2(datarightmr,dataleftmr,salpha) ;
                    %  1 for significant difference found
                    %  0 otherwise

              
              elseif itest == 5 ;    %  Standardized Kolomogorov-Smirnov

                dataleftst = dataleft - mean(dataleft) ;
                datarightst = dataright - mean(dataright) ;

                varleft = var(dataleftst) ;
                varright = var(datarightst) ;

                if varleft > 0 ;
                  dataleftst = dataleftst / varleft ;
                end ;

                if varright > 0 ;
                  datarightst = datarightst / varright ;
                end ;
                
                h = kstest2(datarightst,dataleftst,salpha) ;
                    %  1 for significant difference found
                    %  0 otherwise


              end ;    %  end of itest if-block

              mout(iwid,iloc) = h + 1 ;
                    %  2 for significant difference found
                    %  1 otherwise


            end ;    %  of iloc loop

            
            
          end ;    %  of iwd loop

          
          

          %  Make a plot
          %
          mout = mout + 1 ;
                %  add one, to give entry numbers in colormap
         
          image([minloc,maxloc],[minwid,maxwid],mout) ;
            set(gca,'YDir','normal') ;
              titstr = ['UNC 2000 ScaSta, ' dattit ', ' simultit ...
                                ', gap = ' num2str(gap) ...
                                ', ' testtit] ;
            title(titstr) ;
            xlabel('loc') ;
            ylabel('log10(wid)') ;
       

            orient landscape ;
            outstr = ['\matlabR12\work\NetWorkData\UNCLink2000\UNC2000LinkScaSta1' ...
                             'd' num2str(idat) ...
                             'p' num2str(ipart) ...
                             's' num2str(isimul) ...
                             't' num2str(itest) ...
                             'g' num2str(igap) ...
                                                      ] ;

          eval(['print -dpsc ' outstr]) ;  
          disp('    Finished Plot') ;
    			%gppause ;


        end ;    %  of igap loop

			
			
			end ;    %  of ipart if-block


	
		end ;    %  of itest loop
	
		
		
	end ;    %  of isimul if-block
	

	
end ;    %  of idat loop
	