disp('Running MATLAB script file UNC2001RespSiz2.m') ;
%
%    For making analysis of Response Data, based on UNC traces
%    gathered in April 2001, by Don Smith and Co.
%    and filtered by Felix Hernandez Campos.
%
%    This works with the full data sets
%
%    For further details, see the files:
%         FullRSDataAddress.txt
%         daydatekey.txt


ipart = 24 ;    %  1  -  full data read, and compute simple summaries
               %  2  -  full data read, and Mice and Elephants
               %  3  -  full data read, Elephants
               %  4  -  simulated exponential (same mean)
               %            version of Mice and Elephants (ipart = 2)
               %  5  -  simulated exponential (same mean)
               %            version of Elephants (ipart = 3)
               %  11  -  full data read, log-log scatterplot
               %            for non-0 dur, time vs. size
               %  21  -  full data read, extremal index analysis, response size
               %  22  -  full data read, extremal index analysis, response times
               %  23  -  full data read, extremal index analysis, response rates
               %  24  -  full data read, extremal index analysis, response inverse rates
               %  25  -  full data read, extremal index analysis > 10k, response size
               %  26  -  full data read, extremal index analysis > 10k, response times
               %  27  -  full data read, extremal index analysis > 10k, response rates
               %  28  -  full data read, extremal index analysis > 10k, response inverse rates
               %  31  -  full data read, time series view, response size
               %  32  -  full data read, time series view, response times
               %  33  -  full data read, time series view, response rates
               %  34  -  full data read, time series view, response inverse rates
               %  41  -  full data read, start time. lanbda * (EK+1) estimation
               %  42  -  full data read, start time. log cdf & derivative
               %  43  -  full data read, start time. cluster dist'n cdf
               %  100  -  full data read, and save >10k, non-0 dur.

%vitimeb = 1:21 ;    %  vector of indices of time blocks
%vitimeb = 11 ;    %  vector of indices of time blocks
vitimeb = [11 19] ;    %  vector of indices of time blocks
%  Maybe only Thursday afternoon actually saved as file
                    %  1 - Monday, April 23, 8:00-12:00
                    %  2 - Monday, April 23, 13:00-17:00
                    %  3 - Monday, April 23, 19:30-23:30
                    %  4 - Tuesday, April 24, 8:00-12:00
                    %  5 - Tuesday, April 24, 13:00-17:00
                    %  6 - Tuesday, April 24, 19:30-23:30
                    %  7 - Wednesday, April 25, 8:00-12:00
                    %  8 - Wednesday, April 25, 13:00-17:00
                    %  9 - Wednesday, April 25, 19:30-23:30
                    %  10 - Thursday, April 26, 8:00-12:00
                    %  11 - Thursday, April 26, 13:00-17:00
                    %  12 - Thursday, April 26, 19:30-23:30
                    %  13 - Friday, April 20, 8:00-12:00
                    %  14 - Friday, April 20, 13:00-17:00
                    %  15 - Friday, April 20, 19:30-23:30
                    %  16 - Saturday, April 21, 8:00-12:00
                    %  17 - Saturday, April 21, 13:00-17:00
                    %  18 - Saturday, April 21, 19:30-23:30
                    %  19 - Sunday, April 29, 8:00-12:00
                    %  20 - Sunday, April 29, 13:00-17:00
                    %  21 - Sunday, April 29, 19:30-23:30




%  Load earlier saved indexing, and summary info
%
load('DataUNC2001\UNC2001RespSiz1DatSumm.mat') ;


%  loop through time blocks
%
for itimeb = vitimeb ;    %  loop through time blocks


  if  ipart == 1  | ...
      ipart == 2  | ...
      ipart == 3  | ...
      ipart == 4  | ...
      ipart == 5  | ...
      ipart == 11  | ...
      ipart == 21  | ...
      ipart == 22  | ...
      ipart == 23  | ...
      ipart == 24  | ...
      ipart == 25  | ...
      ipart == 26  | ...
      ipart == 27  | ...
      ipart == 28  | ...
      ipart == 31  | ...
      ipart == 32  | ...
      ipart == 33  | ...
      ipart == 34  | ...
      ipart == 41  | ...
      ipart == 42  | ...
      ipart == 43  | ...
      ipart == 100  ;    %  full data read

    %  Read in main data
    %
    filename = filenamecell{itimeb} ;
    ncfilename = length(filename) ;
        %  number of characters in name
    bigfilename = filename(1:(ncfilename-4)) ;
        %  cutoff last 4 characters
    bigfilename = [bigfilename '.rsps.txt'] ;

    tic ;
    mresp = load(['DataUNC2001\' bigfilename]) ;
    loadtime = toc ;

    n = size(mresp,1) ;



    if ipart == 1 ;    %  full data read, and compute simple summaries


      disp(' ') ;
      disp(['  Loaded ' num2str(n) ' by ' ...
                        num2str(size(mresp,2)) ' matrix in ' ...
                        num2str(loadtime) ' secs']) ;

      flagnon0dur = mresp(:,3) > 0 ;
      nnon0dur = sum(flagnon0dur) ;

      flaggt1500 = (mresp(:,1) > 1500) ;
      ngt1500 = sum(flaggt1500) ;

      flaggt3000 = (mresp(:,1) > 3000) ;
      ngt3000 = sum(flaggt3000) ;

      flaggt10k = (mresp(:,1) > 10000) ;
      ngt10k = sum(flaggt10k) ;

      thruput = mresp(flagnon0dur,1) ./ mresp(flagnon0dur,3) ;
      thruputgt1500 = mresp(flaggt1500,1) ./ mresp(flaggt1500,3) ;
      thruputgt1500n0 = mresp(flaggt1500&flagnon0dur,1) ./ ...
                          mresp(flaggt1500&flagnon0dur,3) ;
      thruputgt3000 = mresp(flaggt3000,1) ./ mresp(flaggt3000,3) ;
      thruputgt3000n0 = mresp(flaggt3000&flagnon0dur,1) ./ ...
                          mresp(flaggt3000&flagnon0dur,3) ;
      thruputgt10k = mresp(flaggt10k,1) ./ mresp(flaggt10k,3) ;
      thruputgt10kn0 = mresp(flaggt10k&flagnon0dur,1) ./ ...
                         mresp(flaggt10k&flagnon0dur,3) ;

      disp(' ') ;
      disp(['  Proportion of zero durations is: ' ...
                  num2str((n - nnon0dur) / n)]) ;


      %  WRITE ANSWERS INTO A FILE
      fid = fopen(['UNC2001RespSiz2SummT' num2str(itimeb) '.txt'],'wt') ;
                %  'wt' is for "delete contents of this file and open 
                %               for writing" (with 't' for "text").
                %  Use 'ay' to append to what is already there, 'rt' for reading.

        titlstr1 = 'Output from the MATLAB Script UNC2001RespSiz2.m,   ' ;
        titlstr1 = [titlstr1,date] ;
      cntbytes = fprintf(fid,'%1s\n\n\n',titlstr1) ;


        titlstr2 = ['     Responses from  ' daystrcell{itimeb} ', ' ...
                                            timestrcell{itimeb}] ;
      cntbytes = fprintf(fid,'%1s\n\n\n',titlstr2) ;


        forstr = '     Total Number of responses = %10.0f\n\n\n\n' ;
      cntbytes = fprintf(fid,forstr,n) ;



        forstr = '     Number of non-0 duration (>1 packet) responses = %10.0f\n\n' ;
      cntbytes = fprintf(fid,forstr,nnon0dur) ;

        forstr = '     Proportion of non-0 duration responses = %6.4f\n\n' ;
      cntbytes = fprintf(fid,forstr,(nnon0dur)/n) ;

        forstr = '     Max Throughput = %10.4g\n\n' ;
      cntbytes = fprintf(fid,forstr,max(thruput)) ;

        forstr = '     Min Throughput = %10.4g\n\n\n\n' ;
      cntbytes = fprintf(fid,forstr,min(thruput)) ;



        forstr = '     Number of responses > 1500 bytes = %10.0f\n\n' ;
      cntbytes = fprintf(fid,forstr,ngt1500) ;

        forstr = '     Proportion of responses > 1500 bytes = %6.4f\n\n' ;
      cntbytes = fprintf(fid,forstr,ngt1500/n) ;

        forstr = '     Number of 0 Duration responses > 1500 bytes = %10.0f\n\n\n' ;
      cntbytes = fprintf(fid,forstr,sum(flaggt1500 & (~flagnon0dur))) ;

        forstr = '     Max Throughput > 1500 bytes = %10.4g\n\n' ;
      cntbytes = fprintf(fid,forstr,max(thruputgt1500)) ;

        forstr = '     Min Throughput > 1500 bytes = %10.4g\n\n\n\n' ;
      cntbytes = fprintf(fid,forstr,min(thruputgt1500)) ;

        forstr = '     Max Throughput > 1500 bytes & non-0 = %10.4g\n\n' ;
      cntbytes = fprintf(fid,forstr,max(thruputgt1500n0)) ;

        forstr = '     Min Throughput > 1500 bytes & non-0 = %10.4g\n\n\n\n' ;
      cntbytes = fprintf(fid,forstr,min(thruputgt1500n0)) ;



        forstr = '     Number of responses > 3000 bytes = %10.0f\n\n' ;
      cntbytes = fprintf(fid,forstr,ngt3000) ;

        forstr = '     Proportion of responses > 3000 bytes = %6.4f\n\n\n' ;
      cntbytes = fprintf(fid,forstr,ngt3000/n) ;

        forstr = '     Number of 0 Duration responses > 3000 bytes = %10.0f\n\n\n' ;
      cntbytes = fprintf(fid,forstr,sum(flaggt1500 & (~flagnon0dur))) ;

        forstr = '     Max Throughput > 3000 bytes = %10.4g\n\n' ;
      cntbytes = fprintf(fid,forstr,max(thruputgt3000)) ;

        forstr = '     Min Throughput > 3000 bytes = %10.4g\n\n\n\n' ;
      cntbytes = fprintf(fid,forstr,min(thruputgt3000)) ;

        forstr = '     Max Throughput > 3000 bytes & non-0 = %10.4g\n\n' ;
      cntbytes = fprintf(fid,forstr,max(thruputgt3000n0)) ;

        forstr = '     Min Throughput > 3000 bytes & non-0 = %10.4g\n\n\n\n' ;
      cntbytes = fprintf(fid,forstr,min(thruputgt3000n0)) ;



        forstr = '     Number of responses > 10 kbytes = %10.0f\n\n' ;
      cntbytes = fprintf(fid,forstr,ngt10k) ;

        forstr = '     Propotion of responses > 10 kbytes = %6.4f\n\n\n' ;
      cntbytes = fprintf(fid,forstr,ngt10k/n) ;

        forstr = '     Number of 0 Duration responses > 10 kbytes = %10.0f\n\n\n' ;
      cntbytes = fprintf(fid,forstr,sum(flaggt10k & (~flagnon0dur))) ;

        forstr = '     Max Throughput > 10k bytes = %10.4g\n\n' ;
      cntbytes = fprintf(fid,forstr,max(thruputgt10k)) ;

        forstr = '     Min Throughput > 10k bytes = %10.4g\n\n\n\n' ;
      cntbytes = fprintf(fid,forstr,min(thruputgt10k)) ;

        forstr = '     Max Throughput > 10k bytes & non-0 = %10.4g\n\n' ;
      cntbytes = fprintf(fid,forstr,max(thruputgt10kn0)) ;

        forstr = '     Min Throughput > 10k bytes & non-0 = %10.4g\n\n\n\n' ;
      cntbytes = fprintf(fid,forstr,min(thruputgt10kn0)) ;



      fclose('all') ;



    elseif  ipart == 2  | ...
            ipart == 4  ;     %  full data read, and Mice and Elephants
                              %  (or exponential modification)


      vstart = mresp(:,2) ;
          %  vector of start times
      vstart = vstart - min(vstart) ;
          %  set so that smallest is 0

      if ipart == 2 ;
        vend = vstart + mresp(:,3) ;
            %  vector of end times
      elseif ipart == 4 ;
        meanlength = mean(mresp(:,3)) ;
        rand('state',14987540986) ;
        randlength = -meanlength * log(rand(n,1)) ;
            %  exponential random variables with same mean
        vend = vstart + randlength ;
            %  vector of end times
      end ;

      mresp = [] ;
          %  to save space


      %  set up graphic preliminaries
      %
      figure(1) ;
      clf ;
      newfrac = 1/4 ;
      nplot = 5000 ;
      ntotal = length(vstart) ;
      rand('state',1923475987435) ;
      

      %  make first page, over full range of data
      %
      left = 0 ;
      right = 4 * 60 * 60 ;
      bottom = 0 ;
      top = 1 ;
      vind = randperm(n)' ;
          %  random re-ordering of indices
      vind = vind(1:nplot) ;
          %  random indices for plotting
      randht = rand(nplot,1) ;
          %  random height for y-coordinate

      plot([vstart(vind)'; vend(vind)'], ...
                      [randht'; randht'],'k-') ;
        axis([left,right,bottom,top]) ;
        if ipart == 2 ;
          title(['Mice & Elephants Plot, ' ...
                      daystrcell{itimeb} ', ' ...
                      timestrcell{itimeb} ', over 4 hours = ' ...
                      num2str(right - left) ' sec'],'FontSize',18) ;
        elseif ipart == 4 ;
          title(['Sim''d Exponential M & E Plot, ' ...
                      daystrcell{itimeb} ', ' ...
                      timestrcell{itimeb} ', over 4 hours = ' ...
                      num2str(right - left) ' sec'],'FontSize',18) ;
        end ;
        xlabel('time (sec)','FontSize',18) ;
        ylabel('Random Height','FontSize',18) ;
          tx = left + 0.2 * (right - left) ;
          ty = bottom + 0.9 * (top - bottom) ;
        if ipart == 2 ;
          text(tx,ty,['Random sample of ' num2str(nplot)],'FontSize',18) ;
        elseif ipart == 4 ;
          text(tx,ty,['Same mean, Random sample of ' num2str(nplot)],'FontSize',18) ;
        end ;
          ty = bottom + 0.8 * (top - bottom) ;
        text(tx,ty,['out of ' num2str(n) ' original responses'],'FontSize',18) ;

      hold on ;
        newleft = (right + left) / 2 - newfrac * (right - left) / 2 ;
        newright = (right + left) / 2 + newfrac * (right - left) / 2 ;
            %  new boundaries at center +- 1/2 of new fraction of all
        plot([[newleft; newleft] , [newright; newright]], ...
             [[bottom; top] , [bottom; top]],'g-') ;
      hold off ;

      if ipart == 2 ;
        pstr = ['UNC2001RS2MiceEleT' num2str(itimeb) 'S'] ;
      elseif ipart == 4 ;
        pstr = ['UNC2001RS2MiceEleT' num2str(itimeb) 'ExpS'] ;
      end ;
      orient landscape ;
      print('-dpsc2',[pstr '1.ps']) ;


      %  iterate through additional pages
      %
      nstep = 6 ;
          %  This is tuned for Thursday Afternoon
          %  If care about other time blocks, may want to tune this
      for istep = 2:nstep ;

        %  reduce data
        %
        keepflag = vstart <= newright ;
            %  ones where vstart <= right, i.e. not "off on right"
        keepflag = keepflag & (newleft <= vend) ;
            %  ones where also not "off on left"
        vstart = vstart(keepflag) ;
        vend = vend(keepflag) ;
        nkeep = length(vstart) ;
        numbracket = sum((vstart <= newleft) & (newright <= vend)) ;
            %  number that bracket full interval
        disp(['   kept ' num2str(nkeep) ' central responses']) ;

        if nkeep < nplot ;
          disp('!!!   Not enough responses left for full plot,') ;
          disp('!!!       terminating execution') ;
          break ;
        end ;


        %  make next plot
        %
        left = newleft ;
        right = newright ;
        bottom = 0 ;
        top = 1 ;
        vind = randperm(nkeep)' ;
            %  random re-ordering of indices
        vind = vind(1:nplot) ;
            %  random indices for plotting
        randht = rand(nplot,1) ;
            %  random height for y-coordinate

        plot([vstart(vind)'; vend(vind)'], ...
                        [randht'; randht'],'k-') ;
          axis([left,right,bottom,top]) ;
          if ipart == 2 ;
            title(['Zoomed Mice & Elephants Plot, ' ...
                        daystrcell{itimeb} ', ' ...
                        timestrcell{itimeb}, ', over ' ...
                        num2str(right - left) ' sec'],'FontSize',18) ;
          elseif ipart == 4 ;
            title(['Sim''d Exp. Zoomed M & E Plot, ' ...
                        daystrcell{itimeb} ', ' ...
                        timestrcell{itimeb}, ', over ' ...
                        num2str(right - left) ' sec'],'FontSize',18) ;
          end ;
            xlabel('time (sec)','FontSize',18) ;
          ylabel('Random Height','FontSize',18) ;
          tx = left + 0.2 * (right - left) ;
          ty = bottom + 0.9 * (top - bottom) ;
        if ipart == 2 ;
          text(tx,ty,['Random sample of ' num2str(nplot)],'FontSize',18) ;
        elseif ipart == 4 ;
          text(tx,ty,['Same mean, Random sample of ' num2str(nplot)],'FontSize',18) ;
        end ;
          ty = bottom + 0.8 * (top - bottom) ;
        text(tx,ty,['from ' num2str(nkeep) ' crossing this range'],'FontSize',18) ;
          ty = bottom + 0.7 * (top - bottom) ;
        text(tx,ty,['out of ' num2str(n) ' original responses'],'FontSize',18) ;
          ty = bottom + 0.6 * (top - bottom) ;
        text(tx,ty,['(' num2str(numbracket) ' of which cover this range)'],'FontSize',18) ;

        hold on ;
          newleft = (right + left) / 2 - newfrac * (right - left) / 2 ;
          newright = (right + left) / 2 + newfrac * (right - left) / 2 ;
              %  new boundaries at center +- 1/2 of new fraction of all
          if istep < nstep ;
            plot([[newleft; newleft] , [newright; newright]], ...
                 [[bottom; top] , [bottom; top]], 'g-') ;
          end ;
        hold off ;


        if ipart == 2 ;
          pstr = ['UNC2001RS2MiceEleT' num2str(itimeb) 'S'] ;
        elseif ipart == 4 ;
          pstr = ['UNC2001RS2MiceEleT' num2str(itimeb) 'ExpS'] ;
        end ;
        orient landscape ;
        print('-dpsc2',[pstr num2str(istep) '.ps']) ;


      end ;



    elseif  ipart == 3  | ...
            ipart == 5  ;     %  full data read, largest Mice and Elephants
                              %  (or exponential modification)


      vstart = mresp(:,2) ;
          %  vector of start times
      vstart = vstart - min(vstart) ;
          %  set so that smallest is 0

      if ipart == 3 ;
        vend = vstart + mresp(:,3) ;
            %  vector of end times
      elseif ipart == 5 ;
        meanlength = mean(mresp(:,3)) ;
        rand('state',14987540986) ;
        randlength = -meanlength * log(rand(n,1)) ;
            %  exponential random variables with same mean
        vend = vstart + randlength ;
            %  vector of end times
      end ;



      [temp,sortind] = sort(mresp(:,3)) ;
      vind = (1:n)' ;
      vind = vind(sortind) ;
          %  indices of data, sorted in increasing duration order
      
      sortind = [] ;
      mresp = [] ;
          %  to save space


      %  set up graphic preliminaries
      %
      figure(1) ;
      clf ;
      vnplot = [5000, 2000, 1000, 500, 200, 100] ;
      maxnplot = max(vnplot) ;
      ntotal = length(vstart) ;
      rand('state',750937457093) ;


      %  make main plots
      %
      left = 0 ;
      right = 4 * 60 * 60 ;
      bottom = 0 ;
      top = 1 ;
      randht = rand(maxnplot,1) ;
          %  random height for y-coordinate
      randht = 1 - randht ;
          %  this simple transformation makes letters 
          %  appear more clearly for this seed

      for nplot = vnplot ;

        vindplot = vind(n-nplot+1:n) ;
        randhtplot = randht(maxnplot-nplot+1:maxnplot) ;

        plot([vstart(vindplot)'; vend(vindplot)'], ...
                        [randhtplot'; randhtplot'],'k-') ;
          axis([left,right,bottom,top]) ;
          if ipart == 3 ;
            title(['Elephants Plot, ' ...
                        daystrcell{itimeb} ', ' ...
                        timestrcell{itimeb} ', over 4 hours = ' ...
                        num2str(right - left) ' sec'],'FontSize',18) ;
          elseif ipart == 5 ;
            title(['Sim''d Exp. Elephants Plot, ' ...
                        daystrcell{itimeb} ', ' ...
                        timestrcell{itimeb} ', over 4 hours = ' ...
                        num2str(right - left) ' sec'],'FontSize',18) ;
          end ;
          xlabel('time (sec)','FontSize',18) ;
          ylabel('Random Height','FontSize',18) ;
            tx = left + 0.2 * (right - left) ;
            ty = bottom + 0.8 * (top - bottom) ;
          text(tx,ty,['Largest ' num2str(nplot) ' of ' ...
                         num2str(n) ' responses'],'FontSize',18) ;
          if ipart == 5 ;
            ty = bottom + 0.6 * (top - bottom) ;
            text(tx,ty,['Simulated from Exponential with same mean'],'FontSize',18) ;
          end ;

        hold on ;
          newleft = (right + left) / 2 - newfrac * (right - left) / 2 ;
          newright = (right + left) / 2 + newfrac * (right - left) / 2 ;
              %  new boundaries at center +- 1/2 of new fraction of all
          plot([[newleft; newleft] , [newright; newright]], ...
               [[bottom; top] , [bottom; top]],'g-') ;
        hold off ;

        if ipart == 3 ;
          pstr = ['UNC2001RS2EleT' num2str(itimeb) ...
                       'p' num2str(nplot)] ;
        elseif ipart == 5 ;
          pstr = ['UNC2001RS2EleT' num2str(itimeb) ...
                       'ExpP' num2str(nplot)] ;
        end ;
        orient landscape ;
        print('-dpsc2',[pstr '.ps']) ;

      end ;



    elseif ipart == 11 ;     %  full data read, log-log scatterplot for non-0 dur

      figure(1) ;
      clf ;

      disp(' ') ;
      disp(['  Loaded ' num2str(n) ' by ' ...
                        num2str(size(mresp,2)) ' matrix in ' ...
                        num2str(loadtime) ' secs']) ;

      flagnon0dur = mresp(:,3) > 0 ;

      vx = mresp(flagnon0dur,1) ;
      vy = mresp(flagnon0dur,3) ;
      nplot = length(vx) ;
      mresp = [] ;
          %  to save space


      %  Set up preliminaries
      %
      timestr = [daystrcell{itimeb} ', ' timestrcell{itimeb}] ;
      titlefontsize = 18 ;
      textfontsize = 18 ;
      labelfontsize = 18 ;
      mindat = min([vx; vy]) ;
      maxdat = max([vx; vy]) ;
          %  overall extrema
      rangedat = maxdat - mindat ;
      rangelogdat = log10(maxdat) - log10(mindat) ;

      xstr = 'Response Size (bytes)' ;
      ystr = 'Response Duration (sec)' ;

      plot(log10(vx),log10(vy),'k.') ;
        axis([log10(mindat),log10(maxdat),log10(mindat),log10(maxdat)]) ;
        axis square ;
        title('log - log Scatterplot, for all non-0 Durations','FontSize',titlefontsize) ;
        xlabel('log_{10}(x)','FontSize',labelfontsize) ;
        ylabel('log_{10}(y)','FontSize',labelfontsize) ;
        hold on ;
            tx = log10(mindat) + 0.05 * rangelogdat ;
            ty = log10(mindat) + 0.85 * rangelogdat ;
          text(tx,ty,timestr,'FontSize',textfontsize) ;
            tx = log10(mindat) + 0.05 * rangelogdat ;
            ty = log10(mindat) + 0.7 * rangelogdat ;
          text(tx,ty,['n = ' num2str(nplot) ' out of ' num2str(n)],'FontSize',textfontsize) ;
            tx = log10(mindat) + 0.05 * rangelogdat ;
            ty = log10(mindat) + 0.55 * rangelogdat ;
          text(tx,ty,['y ~ ' ystr],'FontSize',textfontsize) ;
            tx = log10(mindat) + 0.05 * rangelogdat ;
            ty = log10(mindat) + 0.4 * rangelogdat ;
          text(tx,ty,['x ~ ' xstr],'FontSize',textfontsize) ;
        hold off ;


        pstr = ['UNC2001RS2logSPallT' num2str(itimeb)] ;
        orient landscape ;
        print('-dpsc2',[pstr '.ps']) ;



    elseif  ipart == 21  | ...
            ipart == 22  | ...
            ipart == 23  | ...
            ipart == 24  | ...
            ipart == 25  | ...
            ipart == 26  | ...
            ipart == 27  | ...
            ipart == 28  ;     %  full data read, extremal index analysis,
                               %  based on the estimates (8.11) 
                               %  from Embrecths, Kueppelberg, Mikosch (1997)
                               %  page 429, section 8.1

      figure(1) ;
      clf ;

      disp(' ') ;
      disp(['  Loaded ' num2str(n) ' by ' ...
                        num2str(size(mresp,2)) ' matrix in ' ...
                        num2str(loadtime) ' secs']) ;


      %  set up preliminaries
      %
      vstart = mresp(:,2) ;
          %  vector of start times
      vstart = vstart - min(vstart) ;
          %  set so that smallest is 0

      vsizes = mresp(:,1) ;
      vtimes = mresp(:,3) ;

      if  ipart == 21  | ...
          ipart == 22  | ...
          ipart == 23  | ...
          ipart == 24  ;    %  then keep only non-zero entries
          
        flagn0 = (vsizes > 0) & (vtimes > 0) ;
            %  ones where want to keep data

      elseif  ipart == 25  | ...
              ipart == 26  | ...
              ipart == 27  | ...
              ipart == 28  ;    %  then keep only non-zero entries & > 10k

        flagn0 = (vsizes > 10^4) & (vtimes > 0) ;
            %  ones where want to keep data

      end ;

      vstart = vstart(flagn0) ;
      vsizes = vsizes(flagn0) ;
      vtimes = vtimes(flagn0) ;
      n = sum(flagn0) ;

      [vstart, vind] = sort(vstart) ;      
          %  sort according to start times
      vstart = linspace(1 / (2 * n), 1 - (1 / (2 * n)), n) ;
          %  normalize to be equally spaced on the interval [0,1] 


      if ipart == 21 ;    %  analyze response sizes

        vdata = vsizes(vind) ;
            %  sizes in first column
            %  sorted by start times
        datastr = 'HTTP Resp. Sizes' ;
        datoutstr = 'siz' ;

      elseif ipart == 22 ;    %  analyze response times

        vdata = vtimes(vind) ;
            %  durations in third column
            %  sorted by start times
        datastr = 'HTTP Resp. Times' ;
        datoutstr = 'tim' ;

      elseif ipart == 23 ;    %  analyze response rates

        vdata = vsizes(vind) ./ vtimes(vind) ;
            %  durations in third column
            %  sorted by start times
        datastr = 'HTTP Resp. Rates' ;
        datoutstr = 'rat' ;

      elseif ipart == 24 ;    %  analyze response inverse rates

        vdata = vtimes(vind) ./ vsizes(vind) ;
            %  durations in third column
            %  sorted by start times
        datastr = 'HTTP Resp. Inverse Rates' ;
        datoutstr = 'irat' ;

      elseif ipart == 25 ;    %  analyze response sizes, > 10k

        vdata = vsizes(vind) ;
            %  sizes in first column
            %  sorted by start times
        datastr = 'HTTP Resp. Sizes, > 10k' ;
        datoutstr = 'siz10k' ;

      elseif ipart == 26 ;    %  analyze response times, > 10k

        vdata = vtimes(vind) ;
            %  durations in third column
            %  sorted by start times
        datastr = 'HTTP Resp. Times, > 10k' ;
        datoutstr = 'tim10k' ;

      elseif ipart == 27 ;    %  analyze response rates, > 10k

        vdata = vsizes(vind) ./ vtimes(vind) ;
            %  durations in third column
            %  sorted by start times
        datastr = 'HTTP Resp. Rates, > 10k' ;
        datoutstr = 'rat10k' ;

      elseif ipart == 28 ;    %  analyze response inverse rates, > 10k

        vdata = vtimes(vind) ./ vsizes(vind) ;
            %  durations in third column
            %  sorted by start times
        datastr = 'HTTP Resp. Inverse Rates, > 10k' ;
        datoutstr = 'irat10k' ;

      end ;
      [vsortdata,vindsort] = sort(vdata) ;
          %  sorted in increasing order of data values
      [tempinvsort,vrank] = sort(vindsort) ;
          %  get ranks of data, by inverting sort indices


      vthreshprop = [0.2, 0.1, 0.05, 0.02, 0.01, 0.005, 0.002, 0.001, ...
                       0.0005, 0.0002, 0.0001, 0.00005, 0.00002, 0.00001, ...
                       0.000005, 0.000002, 0.000001] ;
          %  vector of thresholds, as tail proportions
      flagthresh = round(n * vthreshprop) >= 5 ;
          %  one where can have at least 5 data points larger
      vthreshprop = vthreshprop(flagthresh) ;
          %  keep only these thresholds

      vnblocks = logspace(2,log10(n/100),30) ;
          %  vector of numbers of blocks


      %  Main data loop
      %
      vN = [] ;
      mtheta2 = [] ;
      for threshprop = vthreshprop ;

        N = round(n * threshprop) ;
            %  total number of points above threshold
        vN = [vN; N] ;
            %  
        flag =  (vrank > n - N)  ;
            %  ones at data values above threshold
        un = vsortdata(n - N) ;
            %  threshold value (value of last point before threshold)
        vstartthresh = vstart(flag) ;
            %  keep only values where data is above threshold


        vtheta2 = [] ;
        for nblocks = vnblocks ;

          inblocks = round(nblocks) ;
              %  integer version of number of blocks

          disp(['    Working on nblocks = ' num2str(inblocks) ...
                        ' & threshprop = ' num2str(threshprop)]) ;

          binedges = linspace(0,1,inblocks+1) ;
              %  edges of blocks
          vcounts = histc(vstartthresh,binedges) ;
              %  count number of thresholded starting times within blocks
          vcounts = vcounts(1:inblocks) ;
              %  cut off last block (only counts # = last edge)
          K = sum(vcounts > 0) ;
              %  number of vcounts > 0
          theta2 = K / N ;
              %  from formula (8.11)

          vtheta2 = [vtheta2; theta2] ;


          disp(['        Results:  N = ' num2str(N) ...
                                 ', K = ' num2str(K) ...
                                 ', theta2 = ' num2str(theta2)]) ;

        end ;    %  of loop through nblocks


        mtheta2 = [mtheta2, vtheta2] ;


      end ;    %  of loop through threshprop


      %  Make main graphics result
      %
      subplot(2,2,1) ;    %  functions of log10(threshprop)
        plot(log10(vthreshprop'),mtheta2','-') ;
          title(['Extremal Index, ' datastr]) ;
          xlabel('log_{10}(Tail Threshold Proportion)') ;
          ylabel('{\theta}2') ;

      subplot(2,2,3) ;    %  functions of log10(number of blocks)
        plot(log10(vnblocks'),mtheta2,'-') ;
          title(['Extremal Index, ' daystrcell{itimeb} ', ' timestrcell{itimeb}]) ;
          xlabel('log_{10}(Number of Blocks)') ;
          ylabel('{\theta}2') ;

      subplot(2,2,2) ;
        surf(mtheta2)
          title(['Extremal Index, time block ' num2str(itimeb)]) ;

      subplot(2,2,4) ;
        plot(log10(vthreshprop'),log10(vN),'r-',...
                     'LineWidth',2) ;
          title(['# above threshold, N, n = ' num2str(n)]) ;
          xlabel('log_{10}(Tail Threshold Proportion)') ;
          ylabel('log_{10}(N)') ;

      pstr = ['UNC2001RS2ExtInd1T' num2str(itimeb) datoutstr] ;
      orient landscape ;
      print('-dpsc2',[pstr '.ps']) ;



    elseif  ipart == 31  | ...
            ipart == 32  | ...
            ipart == 33  | ...
            ipart == 34  ;     %  full data read, zooming time series view


      figure(1) ;
      clf ;

      disp(' ') ;
      disp(['  Loaded ' num2str(n) ' by ' ...
                        num2str(size(mresp,2)) ' matrix in ' ...
                        num2str(loadtime) ' secs']) ;


      %  set up preliminaries
      %
      vstart = mresp(:,2) ;
          %  vector of start times
      vstart = vstart - min(vstart) ;
          %  set so that smallest is 0

      vsizes = mresp(:,1) ;
      vtimes = mresp(:,3) ;

      flagn0 = (vsizes > 0) & (vtimes > 0) ;
          %  ones where want to keep data

      vstart = vstart(flagn0) ;
      vsizes = vsizes(flagn0) ;
      vtimes = vtimes(flagn0) ;
      n = sum(flagn0) ;

      [vstart, vind] = sort(vstart) ;      
          %  sort according to start times
      vstart = linspace(1 / (2 * n), 1 - (1 / (2 * n)), n) ;
          %  normalize to be equally spaced on the interval [0,1] 


      if ipart == 31 ;    %  analyze response sizes

        vdata = vsizes(vind) ;
            %  sizes in first column
            %  sorted by start times
        datastr = 'HTTP Resp. Sizes' ;
        datoutstr = 'siz' ;

      elseif ipart == 32 ;    %  analyze response times

        vdata = vtimes(vind) ;
            %  durations in third column
            %  sorted by start times
        datastr = 'HTTP Resp. Times' ;
        datoutstr = 'tim' ;

      elseif ipart == 33 ;    %  analyze response rates

        vdata = vsizes(vind) ./ vtimes(vind) ;
            %  durations in third column
            %  sorted by start times
        datastr = 'HTTP Resp. Rates' ;
        datoutstr = 'rat' ;

      elseif ipart == 34 ;    %  analyze response inverse rates

        vdata = vtimes(vind) ./ vsizes(vind) ;
            %  durations in third column
            %  sorted by start times
        datastr = 'HTTP Resp. Inverse Rates' ;
        datoutstr = 'irat' ;

      end ;


      vnzoom = logspace(log10(n),3,4) ;
          %  vector of subsample sizes to zoom through
      top = max(log10(vdata)) ;
      bottom = min(log10(vdata)) ;
      range = top - bottom ;
      top = top + 0.05 * range ;
      bottom = bottom - 0.05 * range ;
      istart = 1 ;
      iend = n ;
      for izoom = 1:4 ;
      
        nzoom = vnzoom(izoom) ;

        subplot(4,1,izoom) ;
          plot(vstart(istart:iend),log10(vdata(istart:iend))) ;
            title([daystrcell{itimeb} ', ' timestrcell{itimeb} ...
                     ', Time Series, nzoom = ' num2str(nzoom)]) ;
            xlabel('time index') ;
            ylabel(datastr) ;
            axis([vstart(istart),vstart(iend),bottom,top]) ;

        if izoom < 4 ;

          nzoom = vnzoom(izoom+1) ;

          istart = ceil((n - nzoom) / 2) ;
          iend = ceil((n + nzoom) / 2) ;

          hold on ;
            plot([vstart(istart); vstart(istart)], ...
                     [bottom; top],'-r') ;
            plot([vstart(iend); vstart(iend)], ...
                     [bottom; top],'-r') ;
          hold off ;
        end ;

      end ;    %  of izoom loop



      pstr = ['UNC2001RS2TS1T' num2str(itimeb) datoutstr] ;
      orient tall ;
      print('-dpsc2',[pstr '.ps']) ;



    elseif  ipart == 41  | ...
            ipart == 42  | ...
            ipart == 43  ;     %  full data read, start time analysis


      figure(1) ;
      clf ;

      disp(' ') ;
      disp(['  Loaded ' num2str(n) ' by ' ...
                        num2str(size(mresp,2)) ' matrix in ' ...
                        num2str(loadtime) ' secs']) ;


      %  set up preliminaries
      %
      vstart = mresp(:,2) ;
          %  vector of start times
      vstart = vstart - min(vstart) ;
          %  set so that smallest is 0

      mresp = [] ;
          %  to save space

      [vstart, vind] = sort(vstart) ;      
          %  sort according to start times

      maxtime = vstart(n) ;


      if ipart == 41 ;    %  full data read, start time, 
                          %        lanbda * (EK+1) estimation


        tgrid = linspace(maxtime/2,maxtime,200) ;

%        Noft = n * cprobSM(vstart,tgrid,0) ;
%            %  scale up "probabilities", 0 for presorted
%        this was too slow with these big data sets, hence replace by:
%
        vprob = linspace(1/(n+1),1 - (1/(n+1)),n) ;
        vstartdiff = vstart(2:n) - vstart(1:(n-1)) ;
        flaglastsame = [0; (vstartdiff == 0)] ;
            %  ones where the last element is the same as this one
        vstartlastdiff = vstart(~flaglastsame) ;
        vproblastdiff = vprob(~flaglastsame) ;
        NofT = interp1(vstartlastdiff,vproblastdiff,tgrid) ;
            %  essentially main line from cprobSM

        plot(tgrid,Noft ./ t,'r-','LineWidth',2) ;
          title(['Est. of \lambda E(K+1), HTTP Starts, ' ...
                   daystrcell{itimeb} ', ' timestrcell{itimeb} ...
                   ', n = ' num2str(n)],'FontSize',18) ;
          xlabel('time','FontSize',18) ;
          ylabel('N(t) / t','FontSize',18) ;


        pstr = ['UNC2001RS2startavg1T' num2str(itimeb) datoutstr] ;
        orient landscape ;
        print('-dpsc2',[pstr '.ps']) ;



      elseif ipart == 42  ;     %  full data read, start time, log cdf
                                %        and derivative


        vintarr = vstart(2:n) - vstart(1:(n-1)) ;

        maxia = max(vintarr) ;
        tgrid = linspace(1/201,maxia,201) ;

%        vlccdf = cprobSM(vstart,tgrid,0) ;
%        this was too slow with these big data sets, hence replace by:
%
        vprob = linspace(1/(n+1),1 - (1/(n+1)),n) ;
        vstartdiff = vstart(2:n) - vstart(1:(n-1)) ;
        flaglastsame = [0; (vstartdiff == 0)] ;
            %  ones where the last element is the same as this one
        vstartlastdiff = vstart(~flaglastsame) ;
        vproblastdiff = vprob(~flaglastsame) ;
        vlccdf = interp1(vstartlastdiff,vproblastdiff,tgrid) ;
            %  essentially main line from cprobSM

        vlccdf = log(1 - vlccdf) ;


        tcent = mean([tgrid(1:200), tgrid(2:201)],2) ;
        tdif = tgrid(2:201) - tgrid(1:200) ;
        vlambda = (vslope(2:201) - vslope(1:200)) ./ tdif ;
            %  vector of estimated lambdas

        vek = ((n / maxtime) ./ vlambda) - 1 ;
            %  vector of estimated EKs



        subplot(3,1,1) ;
          plot(tgrid,vlccdf,'b-','LineWidth',2) ;
            title(['1 - CDF of HTTP Start Interarrivals, ' ...
                     daystrcell{itimeb} ', ' timestrcell{itimeb}], ...
                     'FontSize',18) ;
            xlabel('Interarrival Time','FontSize',18) ;
            ylabel('log(CCDF)','FontSize',18) ;


        subplot(3,1,2) ;
          plot(tcent,vlambda,'b-','LineWidth',2) ;
            title(['Slope of 1 - CDF, gives estimate of \lambda' ...
                     ', n = ' num2str(n)],'FontSize',18) ;
            xlabel('Interarrival Time','FontSize',18) ;
            ylabel('Estimated \lambda','FontSize',18) ;


        subplot(3,1,3)
          plot(tcent,vek,'b-','LineWidth',2) ;
            title(['Estimated Expected Number of cluster points' ...
                     ', n = ' num2str(n)],'FontSize',18) ;
            xlabel('Interarrival Time','FontSize',18) ;
            ylabel('Estimated EK','FontSize',18) ;



        pstr = ['UNC2001RS2startCCDF1T' num2str(itimeb) datoutstr] ;
        orient tall ;
        print('-dpsc2',[pstr '.ps']) ;



      elseif ipart == 43  ;     %  full data read, start time, cdf
                                %        of cluster dist'n


        vintarr = vstart(2:n) - vstart(1:(n-1)) ;

        maxia = max(vintarr) ;
        tgrid = linspace(1/401,maxia,401) ;

%        vlccdf = cprobSM(vstart,tgrid,0) ;
%        this was too slow with these big data sets, hence replace by:
%
        vprob = linspace(1/(n+1),1 - (1/(n+1)),n) ;
        vstartdiff = vstart(2:n) - vstart(1:(n-1)) ;
        flaglastsame = [0; (vstartdiff == 0)] ;
            %  ones where the last element is the same as this one
        vstartlastdiff = vstart(~flaglastsame) ;
        vproblastdiff = vprob(~flaglastsame) ;
        vlccdf = interp1(vstartlastdiff,vproblastdiff,tgrid) ;
            %  essentially main line from cprobSM

        vlccdf = log(1 - vccdf) ;


        tcent = mean([tgrid(1:400), tgrid(2:401)],2) ;
        tdif = tgrid(2:401) - tgrid(1:400) ;
        vlambda = (vintarr(2:401) - vintarr(1:400)) ./ tdif ;
            %  vector of estimated lambdas

        vek = ((n / maxtime) ./ vlambda) - 1 ;
            %  vector of estimated EKs


        viacent = mean([vintarr(1:400), vintarr(2:401)],2) ;
        vintccdf = viacent ./ tdif ;
            %  Areas of trapezoids
        vintccdf = cumsum(flipud(vintccdf)) ;
        vintccdf = flipud(vintccdf) ;
            %  numerical integral of tail of ccdf

        vintclustcdf = ((vlambda ./ tcent) + ...
                           log(vlambda .* (vek + 1) .* vintccdf)) ...
                                          ./ vek ;
            %  vector of integral of tail of cluster cdf
        tcentd = mean([tcent(1:399), tcent(2:400)],2) ;
        vclustcdf = (vintclustcdf(2:400) - vintclustcdf(1:399)) ./ ...
                              (tcent(2:400) - tcent(1:399)) ;
            %  numerically differentiate to get tail cdf of cluster dist'n
        vclustcdf = 1 - vclustcdf ;


        plot(tcentd,vclustcdf,'b-','LineWidth',2) ;
          title(['CDF of Cluster Distribution, ' ...
                   daystrcell{itimeb} ', ' timestrcell{itimeb}], ...
                   'FontSize',18) ;
          xlabel('Time','FontSize',18) ;
          ylabel('CDF','FontSize',18) ;


        pstr = ['UNC2001RS2startClustCDF1T' num2str(itimeb) datoutstr] ;
        orient landscape ;
        print('-dpsc2',[pstr '.ps']) ;



      end ;



    elseif ipart == 100 ;     %  full data read, and save >10k, non-0 dur.

      disp(' ') ;
      disp(['  Loaded ' num2str(n) ' by ' ...
                        num2str(size(mresp,2)) ' matrix in ' ...
                        num2str(loadtime) ' secs']) ;

      flagnon0dur = mresp(:,3) > 0 ;
      flaggt10k = (mresp(:,1) > 10000) ;

      saveflag = flagnon0dur & flaggt10k ;

      mresp = mresp(saveflag,:) ;


      %  save pared responses in a .mat file
      %
      savestr =['DataUNC2001\UNC2001RSallBigT' num2str(itimeb) '.mat'] ;
      save(savestr,'mresp') ;


      disp(['  saved ' num2str(sum(saveflag)) ' flows > 10k with non-0 duration,']) ;
      disp(['    for time block ' num2str(itimeb), ...
                  daystrcell{itimeb} ', ' timestrcell{itimeb}]) ;



    end ;    %  of ipart if-block


  end ;    %  of ipart if-block


end ;    %  of itimeb loop through time blocks


beep ;
