disp('Running MATLAB script file VarHeavyTail2figB.m') ;

%    Makes Figure B for the Variable Heavy Tails paper

%    Copied from UNC2001RespSiz1.m
%    with parameters set to:
%        ipart = 43  -  Full Data Pareto Q-Q
%        vitimeb = 10 ;    %  vector of indices of time blocks

%
%    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.
%
%    Because of the size of the data, they are summarized as:
%         Summary statistics:
%              rsp01power.xls
%              rsp01morestats.xls  (with trimmed down version in 
%                                           rsp01morestatsraw.xls
%                                   and just data saved in
%                                           rsp01morestatsraw.txt)
%         Quantiles of body of distribution:
%              responsesizequantiles.xls
%         Individual responses:
%              20010420_1300.raw
%                     to
%              20010429_800.raw
%                (21 time blocks)
%
%    For further details, see the files:
%         responsesizetailsnotes.txt
%         daydatekey.txt


icolor = 1 ;    %  0 - black and white only
                %  1 - original full color


%vitimeb = 1:21 ;    %  vector of indices of time blocks
vitimeb = 10 ;    %  vector of indices of time blocks
                    %  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




figure(1) ;
clf ;


%  Load summary data (gives filenames, etc.)
%
%load('\Research\NetworkData\ResponseSize\DataUNC2001\UNC2001RespSiz1DatSumm.mat') ;
load('..\ResponseSize\DataUNC2001\UNC2001RespSiz1DatSumm.mat') ;


%  Load Quantiles
%
%fid = fopen('\Research\NetworkData\ResponseSize\DataUNC2001\responsesizequantilesRaw.txt','rt') ;
fid = fopen('..\ResponseSize\DataUNC2001\responsesizequantilesRaw.txt','rt') ;
        %  'rt' is for "read only" and "text"
fstr = '%f' ;
[mq,cnt] = fscanf(fid,fstr) ;
fclose('all') ;
mq = reshape(mq,21,1000)' ;
        %  use this for 1000x21 matrix, since reshape works column wise


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


  %  first get part that appeared earlier
  %  (i.e. for given quantiles)
  %
  vqdata = mq(:,virowraw(itimeb)) ;
      %  unpack vector of small fixed set of quantiles

  n = vnumresp(itimeb) ;
      %  unpack full sample size

  nq = length(vqdata) ;
  vccdf = flipud((1/n) + (0:(nq-1))' / nq) ;


  %  load large data values, and integrate
  %
  filename = filenamecell{itimeb} ;
%  mtailresp = load(['\Research\NetworkData\ResponseSize\DataUNC2001\' filename]) ;
  mtailresp = load(['..\ResponseSize\DataUNC2001\' filename]) ;
  vtailsize = mtailresp(:,1) ;
  ntail = length(vtailsize) ;

  qdata = [vqdata; sort(vtailsize)] ;
  vccdf = [vccdf; flipud((1:ntail)' / n)] ;
  
  [qdata, vindsort] = sort(qdata) ;
  vccdf = vccdf(vindsort) ;
  nq = length(qdata) ;
  

  %  Use lines copied and modified from qqSM.m
  %
%    vqalign = [0.99; 0.99999] ;
%    vqalign = [0.9; 0.9999] ;
%    vqalign = [0.9; 0.999] ;
%    vqalign = [0.9; 0.99] ;
%    vqalign = [0.5; 0.99] ;
  vqalign = [0.8; 0.99] ;
  maxstep = 1000 ;
  relaccthreshold = 10^(-6) ;
  if icolor == 0 ;
    colorcell = {'k' 'k' [0.2 0.2 0.2]} ;
            %  1st - data
            %  2nd - line
            %  3rd - simulated versions
    ltypestr = '--' ;
    simltypestr = ':' ;
  else ;
    colorcell = {'r' 'g' 'b'} ;
            %  1st - data
            %  2nd - line
            %  3rd - simulated versions
    ltypestr = '-' ;
    simltypestr = '-' ;
  end ;
  titlestr = ['HTTP R. S., Pareto Q-Q, ' num2str(vqalign(1)) ...
                     ' & ' num2str(vqalign(2)) ... 
                     ' quant. matched, n = ' num2str(n)] ;
%  titlefontsize = 18 ;
  titlefontsize = 12 ;
  diststr = 'Pareto-log' ;
%  xlabelstr = ', Q' ;
%  xlabelstr = [diststr, xlabelstr] ;
  xlabelstr = 'log_{10}(Pareto Quantiles)' ;
%  ylabelstr = 'Data Q' ;
%  ylabelstr = ['log ' ylabelstr] ;
  ylabelstr = 'log_{10}(Observed Quantiles)' ;
%  labelfontsize = 18 ;
  labelfontsize = 12 ;
%  parfontsize = 18 ;
  parfontsize = 12 ;
  vshowq = [0.1, 0.5, 0.9, 0.99, 0.999, 0.9999, 0.99999] ;
  if icolor == 0 ;
    nsim = 40 ;
  else ;
    nsim = 100 ;
  end ;
  nsimplotval = 900 ;


  %  estimate Pareto parameters
  %
  p1 = vqalign(1) ;
  p2 = vqalign(2) ;
  vq = cquantSM(vqdata,vqalign,0) ;
       %  0 to indicate presorted data
       %  (only work with small fixed quantiles)
  q1hat = vq(1) ;
  q2hat = vq(2) ;

  numera = log(1 - p1) - log(1 - p2);
  oldsigma = 0 ;
        %  common and starting values

  for istep = 1:maxstep ;

    q1z = q1hat + oldsigma ;
    q2z = q2hat + oldsigma ;

    alpha = numera / (log(q2z) - log(q1z)) ;
    sigma = q1z * (1 - p1)^(1/alpha) ;

    relacc = abs(sigma - oldsigma) / sigma ;

    disp(['     For iteration step ' num2str(istep) ...
                ', relacc = ' num2str(relacc) ...
                ', alpha = ' num2str(alpha) ...
                ', sigma = ' num2str(sigma)]) ;

    
    if relacc < relaccthreshold ;
      errflag = 0 ;
      disp(['     For iteration step ' num2str(istep) ...
                  ', relacc = ' num2str(relacc) ...
                  ', alpha = ' num2str(alpha) ...
                  ', sigma = ' num2str(sigma)]) ;
      break ;
    else ;
      errflag = 1 ;
      oldsigma = sigma ;
    end ;

  end ;

  if errflag ~= 0 ;
    disp('!!!   Warning:  Pareto fit may be unstable   !!!') ;
  end ;

  ididqalign = 1 ;


  %  for these estimated parameters, get theoretical quantiles
  %
  qtheory = (vccdf.^(-1/alpha) - 1) * sigma ;
  pgrid = (1:n)' / (n + 1) ;
  qtheoryall = ((1 - pgrid).^(-1/alpha) - 1) * sigma ;


  %  put on log scale
  %
  qdata = log(qdata) ;
  qtheory = log(qtheory) ;
    q1hat = log(q1hat) ;
    q2hat = log(q2hat) ;
  qtheoryall = log(qtheoryall) ;

  left = qtheory(1) ; 
  right = qtheory(nq) ; 
  bottom = qdata(1) ;
  top = qdata(nq) ;

  par1str = '\alpha' ;
  par1val = alpha ;
  par2str = '\sigma' ;
  par2val = sigma ;


  %  Make main graphic in current axes
  %
  plot(qtheory,qdata, ...
      [colorcell{1} '-'], ...
      'LineWidth',3) ;
    axis equal ;
    axis([left,right,bottom,top]) ;

%    th = title(titlestr) ;
%    set(th,'FontSize',titlefontsize) ;

    xlh = xlabel(xlabelstr) ;    
    ylh = ylabel(ylabelstr) ;    
    set(xlh,'FontSize',labelfontsize) ;
    set(ylh,'FontSize',labelfontsize) ;


  %    overlay 45 degree line
  %
  minmin = min(left,bottom) ;
  maxmax = max(right,top) ;
  hold on ;
    plot([minmin,maxmax],[minmin,maxmax], ...
         [colorcell{2} ltypestr], ...
         'LineWidth',2) ;
  hold off ;


  %  Add simulated realizations (if needed)
  %
  simseed = 34875094 ;
  rand('seed',simseed) ;
      %  set seed

  for isim = 1:nsim ;    %  loop through simulation steps

    disp(['      working on simulated data set ' num2str(isim) ...
                                          ' of ' num2str(nsim)]) ;

    vsimdat = sigma * ((1 - rand(n,1)).^(-1/alpha) - 1) ;
    vsimdat = log(vsimdat) ;

    vsimdat = sort(vsimdat) ;
        %  sort this column

    nspvo3 = floor(nsimplotval / 3) ;
          %  one third of total, to put at each end
    vindleft = (1:nspvo3)' ;
          %  left end  (include all)
    vindright = ((n-nspvo3+1):n)' ;
          %  right end  (include all)

    nspvlo = nsimplotval - length(vindleft) ...
                         - length(vindright) ;
          %  number of grid points left over (for use in middle)
    vindmid = round(linspace(nspvo3+1,n-nspvo3,nspvlo)') ;

    vind = [vindleft; vindmid; vindright] ;
          %  vector of indices, full

    qtheoryp = qtheoryall(vind,:) ;
    vsimdatp = vsimdat(vind,:) ;

    hold on  ;
      plot(qtheoryp,vsimdatp,simltypestr,'Color',colorcell{3}) ;
    hold off ;

  end ;    %  of loop through simulated data sets


  %  replot important stuff
  %
  hold on ;
    plot(qtheory,qdata, ...
         [colorcell{1} '-'], ...
         'LineWidth',3) ;
    plot([minmin,maxmax],[minmin,maxmax], ...
         [colorcell{2} ltypestr], ...
         'LineWidth',2) ;
  hold off ;


  %  Add text for parameters and time block
  %
    tx = left + 0.1 * (right - left) ;
    ty = bottom + 0.9 * (top - bottom) ;
  th = text(tx,ty,[par1str ' = ' num2str(par1val)]) ;
  set(th,'FontSize',parfontsize) ;

    ty = bottom + 0.8 * (top - bottom) ;
  th = text(tx,ty,[par2str ' = ' num2str(par2val)]) ;
  set(th,'FontSize',parfontsize) ;

%    tx = left + 0.1 * (right - left) ;
%    ty = bottom + 0.9 * (top - bottom) ;
%  th = text(tx,ty,[daystrcell{itimeb} ', ' timestrcell{itimeb}]) ;
%  set(th,'FontSize',parfontsize) ;


  %  show fit quantiles
  %
  hold on ;
    sch = plot([q1hat; q2hat], [q1hat; q2hat], 'ko') ;
  hold off ;


  %  show additional quantiles
  %
  for isq = 1:length(vshowq) ;

    sq = vshowq(isq) ;
  
    if  0 < sq   &  sq < 1  ;    %  then have valid probability

      [temp,qi] = min(abs((1 - vccdf) - sq)) ;
          %  gets index of this quantile
          %  i.e. where prob is closest to sq prob

      qehat = qdata(qi) ;
      qthat = qtheory(qi) ;

      qhx = qthat + 0.02 * (right - left) ;

      hold on ;
        plot([qthat], [qehat], 'k+') ;
        th = text(qhx,qehat,[num2str(sq) ' quantile']) ;
          set(th,'FontSize',parfontsize) ;
      hold off ;

    else ;    %  then don't have valid probability, give error message

      disp('!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!') ;
      disp('!!!   Warning from qqSM:                       !!!') ;
      disp('!!!   Invalid entry (not in (0,1)) in vshowq   !!!') ;
      disp('!!!   So won''t show this quantile              !!!') ;
      disp('!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!') ;

    end ;    %  of sq if-block  


  end ;    %    of isq loop through quantiles to show




%      orient landscape ;
%    print('-dpsc2',['UNC2001RS1allQQpar' num2str(itimeb) '.ps']) ;


end ;    %  of itimeb loop




if icolor == 0 ;
  pstr = 'VarHeavyTail2figBbw' ;
else ;
  pstr = 'VarHeavyTail2figB' ;
end ;
  orient portrait ;
  set(gcf,'PaperSize',[6.0, 4.0]) ; 
  set(gcf,'PaperPosition',[0.25, 0.25, 5.5, 3.5]) ; 
print('-dpsc', [pstr '.ps']) ;




