disp('Running MATLAB script file UNC2000LinkScaSta1.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 ;
addpath \matlabR12\work\SteveFun -end ;
addpath \matlabR12\work\SteveSmo -end ;


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

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

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


vitest = [1 2 3 4 5 6 7 8] ;
               %  1 - means and t-test
               %  2 - Kolmogorov-Smirnov
               %  3 - variances and F-test
               %  4 - mean removed Kolmogorov-Smirnov
               %  5 - standardized Kolmogorov-Smirnov
               %  6 - median removed Kolmogorov-Smirnov
               %  7 - median and MAD st'dized Kolmogorov-Smirnov
               %  8 - Quartile st'dized Kolmogorov-Smirnov



%  Set common parameters
%  
nloc = 256 ;
%nloc = 19;
          %  number of locations to consider, for each scale
alpha = 0.05 ;
nminew = 8 ;
          %  minimum number in each window
nrepeat = 3 ;
          %  number of times to repeat movie (while showing)
fps = 2 ;
          %  frames per second, for movie (both show and save)

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

figh = figure(1) ;
clf ;



%  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' ;  
  end;
  
  

	%  loop through tests
	%
	for itest = vitest ;
	
	
	
		if itest == 1 ;    %  1 - means and t-test
          testtit = 'Mean t-test' ;      
		elseif itest == 2 ;    %  2 - Kolmogorov-Smirnov
          testtit = 'Kolm-Smir' ;      
		elseif itest == 3 ;    %  3 - variances and F-test
            testtit = 'Var. F-test' ;      
		elseif itest == 4 ;    %  4 - mean removed Kolmogorov-Smirnov
            testtit = 'Mean R. K-S' ;      
		elseif itest == 5 ;    %  5 - standardized Kolmogorov-Smirnov
            testtit = 'Std''zed K-S' ;  
    elseif itest == 6;    %  6 - median removed Kolmogorov-Smirnov
            testtit = 'Median R. K-S' ;    
		elseif itest == 7 ;    %  7 - median and MAD st'dized Kolmogorov-Smirnov
            testtit = 'MAD S''d K-S' ;  
		elseif itest == 8 ;    %  8 - Quartile st'dized Kolmogorov-Smirnov
            testtit = 'Q S''d K-S' ;  
		end ;    %  end of itest if-block
		
	
		
		
		%  loop through data sets
		%
		for idat = vidat ;
		
		
		
			%  Read in Data
			%
      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]) ;


      vx = bincts;
      %??? fix this, in terms of data, and adjust for time
      
      %  set parameters
      %
			n = nbin ;
			%  left comes from load (secs)
			%  right comes from load (secs)
			binwid = (right - left) / nbin ;
			%vgap = [binwid 0.02 50] ;
			vgap = logspace(-3,2,6) ;
			    %  corresponds to time ~ 0.02 secs (hopefully more than round trip time)
		  ngap = length(vgap) ;
			mingap = min(vgap) ;
			nwid = log2(nbin) - 1 ;
                          %  number of widths to show at each scale
      tminew = nminew * binwid ;
                          %  minimum time span for each window
			tgrid = linspace(left,right,n)' ;



			
			

				
				
      %  loop through gaps
      %
      clear moviestruct ;
      for igap = 1:ngap ;

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


        %  loop through wid's
        %
        mout = zeros(nwid,nloc) ;

        if ipart == 1 ;    %  just making static plot, so use full size
        
          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) ;

        elseif ipart == 2 ;    %  just making static plot, so use full size

          minwid = log10(tminew) ;
          maxwid = log10((right - left - mingap) / 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 + mingap / 2) ;
          maxloc = (right - tminew - mingap / 2) ;
        
        end ;
          
          
        vloc = linspace(minloc,maxloc,nloc) ;
        for iwid = 1:nwid ;

          wid = vwid(iwid) ;

          %if iwid == 1 ;
            disp(['  Working on wid ' num2str(iwid) ' = ' ...
                           num2str(wid) ...
                           '    for ' simultit ', ' testtit ', ' ...
                           dattit ', gap ' num2str(gap)]) ;
          %end ;
                           
          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 ;
						%mglobal = 2 * nloc - 1 ;
						mglobal = 2 * min((right - left)/tminew,nloc) - 1 ;
						salpha = 1 - (1 - alpha)^(1/mglobal) ; 
          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 = vx(flagleft) ;
            dataright = vx(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 ;    %  Kolmogorov-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 Kolmogorov-Smirnov

              dataleftmr = dataleft - mean(dataleft) ;
              datarightmr = dataright - mean(dataright) ;
              
              h = kstest2(datarightmr,dataleftmr,salpha) ;
                  %  1 for significant difference found
                  %  0 otherwise
                  %  Note:  have not adjusted significance for removal of mean

                
            elseif itest == 5 ;    %  Standardized Kolmogorov-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
                  %  Note:  have not adjusted significance for standardization

                  
            elseif itest == 6 ;    %  median removed Kolmogorov-Smirnov

              dataleftmr = dataleft - median(dataleft) ;
              datarightmr = dataright - median(dataright) ;
              countleft    =  length(dataleft);
              countright     =  length(dataright);
              ncoefficient      =  2*...
                    countleft * countright /(countleft + countright);
              
              % The sqrt(2) comes from the fact that I have 2 Brownian
              % bridges on the interval of length .5. The same applies
              % for the sqrt in the calculation of p-value.

              %              
              % Calculate F1(x) and F2(x), the empirical (i.e., sample) CDFs.
							%
							
							binEdges    =  [-inf ; sort([dataleftmr;datarightmr]) ; inf];
							
							binCounts1  =  histc (dataleftmr , binEdges);
							binCounts2  =  histc (datarightmr , binEdges);
							
							sumCounts1  =  cumsum(binCounts1)./sum(binCounts1);
							sumCounts2  =  cumsum(binCounts2)./sum(binCounts2);
							
							sampleCDF1  =  sumCounts1(1:end-1);
							sampleCDF2  =  sumCounts2(1:end-1);
							
							
							%
							% Compute the test statistic of interest.
							%
							
              deltaCDF  =  abs(sampleCDF1 - sampleCDF2);
              KSstatistic   =  max(deltaCDF);
              lambda =  (sqrt(ncoefficient) + 0.12 + 0.11/sqrt(ncoefficient))...
                                            * KSstatistic ;
              
              
              %
							%  Use the asymptotic Q-function to approximate the 2-sided P-value.
							%
              j       =  [1:101]';
              p  =  sqrt(2 * sum((-1).^(j-1).*exp(-2*lambda*lambda*j.^2)));

              h = 0;
							if p <= salpha, 
                h = 1; 
              end ;
                          %  1 for significant difference found
                          %  0 otherwise
                          %  Note:  HAVE adjusted significance for removal of median

                
            elseif itest == 7 ;    %  median and MAD st'dized Kolmogorov-Smirnov

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

              gpmadleft = gpmad(dataleftst) ;
              gpmadright = gpmad(datarightst) ;

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

              if gpmadright > 0 ;
                datarightst = datarightst / gpmadright ;
              end ;
              
              h = kstest2(datarightst,dataleftst,salpha) ;
                  %  1 for significant difference found
                  %  0 otherwise
                  %  Note:  have not adjusted significance for standardization

                  
            elseif itest == 8 ;    %  Quartile st'dized Kolmogorov-Smirnov

              medianleft = median(dataleft) ;
              medianright = median(dataright) ;
              
              dataleftst = dataleft - medianleft ;
              datarightst = dataright - medianright ;

              Qoneleft =  prctile(dataleftst,25) ;
              Qoneright =  prctile(datarightst,25) ;

              Qthreeleft =  prctile(dataleftst,75) ;
              Qthreeright =  prctile(datarightst,75) ;

              flagoneleft = dataleftst < 0 ;
              flagoneright = datarightst < 0 ;

              if Qoneleft < 0 ;    %  then rescale guys on Q1 side
                dataleftst(flagoneleft) = phiinv(0.25) * ...
                                    dataleftst(flagoneleft) / Qoneleft ;
              end ;

              if Qoneright < 0 ;    %  then rescale guys on Q1 side
                datarightst(flagoneright) = phiinv(0.25) * ...
                                    datarightst(flagoneright) / Qoneright ;
              end ;

              if Qthreeleft > 0 ;    %  then rescale guys on Q3 side
                dataleftst(~flagoneleft) = phiinv(0.75) * ...
                                    dataleftst(~flagoneleft) / Qthreeleft ;
              end ;

              if Qthreeright > 0 ;    %  then rescale guys on Q3 side
                datarightst(~flagoneright) = phiinv(0.75) * ...
                                    datarightst(~flagoneright) / Qthreeright ;
              end ;
              
              countleft    =  length(dataleftst);
              countright     =  length(datarightst);
              ncoefficient      =  4*...
                     countleft * countright /(countleft + countright);

              % The sqrt(4) comes from the fact that I have 4 Brownian
              % bridges on the subintervals of length .25. The same applies
              % for the sqrt in the calculation of p-value.
    
                     
              %
							% Calculate F1(x) and F2(x), the empirical (i.e., sample) CDFs.
							%
							
							binEdges    =  [-inf ; sort([dataleftst;datarightst]) ; inf];
							
							binCounts1  =  histc (dataleftst , binEdges);
							binCounts2  =  histc (datarightst , binEdges);
							
							sumCounts1  =  cumsum(binCounts1)./sum(binCounts1);
							sumCounts2  =  cumsum(binCounts2)./sum(binCounts2);
							
							sampleCDF1  =  sumCounts1(1:end-1);
							sampleCDF2  =  sumCounts2(1:end-1);
							
							
							%
							% Compute the test statistic of interest.
							%
							
              deltaCDF  =  abs(sampleCDF1 - sampleCDF2);
              KSstatistic   =  max(deltaCDF);
              lambda =  (sqrt(ncoefficient) + 0.12 + 0.11/sqrt(ncoefficient))...
                                            * KSstatistic ;
              
        
              %
							%  Use the asymptotic Q-function to approximate the 2-sided P-value.
							%
              j       =  [1:101]';
              p  =  (2 * sum((-1).^(j-1).*exp(-2*lambda*lambda*j.^2)))^(0.25) ;

              h = 0;
							if p <= salpha, 
                h = 1; 
              end ;
                          %  1 for significant difference found
                          %  0 otherwise
                          %  Note:  HAVE adjusted significance for standardization

                          
            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 = ['ScaSta, ' dattit ', ' simultit ...
                              ', gap = ' num2str(gap) ...
                              ', ' testtit] ;
          title(titstr) ;
          xlabel('loc') ;
          ylabel('log10(wid)') ;
            
			
				if ipart == 1 ;    %  then save as a postscript file
            
            orient landscape ;
            outstr = ['\matlabR12\work\NetWorkData\UNCLink2000\UNC2000LinkScaSta1' ...
                             'p' num2str(ipart) ...
                             's' num2str(isimul) ...
                             'd' num2str(idat) ...
                             't' num2str(itest) ...
                             'g' num2str(igap) ...
                                                      ] ;

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


  			elseif ipart == 2 ;    %  then save movie frame
    			
    			moviestruct(igap) = getframe(figh) ;

			
				end ;    %  of ipart if-block

				
				
      end ;    %  of igap loop

      

			if ipart == 2 ;    %  then save and run movie

			
			  %  first reorder movie to make play back and forth
			  %
        vmorder = [(1:ngap),((ngap-1):-1:2)] ;
        moviestruct = moviestruct(vmorder) ;
        %  reorder frames, to loop back to beginning
			  
        
        %  play movie on screen
        %
			  movie(figh,moviestruct,nrepeat,fps) ;
			  
			  
			  %  save movie as avi file
			  %
        outstr = ['\matlabR12\work\NetWorkData\UNCLink2000\UNC2000LinkScaSta1' ...
                         'p' num2str(ipart) ...
                         's' num2str(isimul) ...
                         'd' num2str(idat) ...
                         't' num2str(itest) ...
                                                      ] ;
			  movie2avi(moviestruct,outstr,'compression','MSVC', ...
			                               'keyframe',fps,'fps',fps) ;
			

				
			end ;


	

    end ;    %  of idat loop

	
	
	end ;    %  of itest loop

	
	
end ;    %  of isimul if-block