disp('Running MATLAB script file gpchetat.m') ;
%
%    FOR DEVELOPMENT AND TESTING OF MATLAB FUNCTION gpchetat,
%    General Purpose Curve Heights at Event Times Analysis

ipart = 2 ;    %  1 - example data from egcheta1.m
               %  2 - example from code decay project
               %  3 - example based on volcano data
               %  4 - example based on volcano data, modified 
               %               "detrending", i.e. subtracting mean


if ipart == 1 ;    %  then work with simple egcheta1.m example

  %  Create curve and artificial data
  %
%  tgrid = linspace(0,.5,401)' ;
%  tgrid = linspace(.4,1,401)' ;
          %  used these to check truncation action
  tgrid = linspace(0,1,401)' ;
  cur = phi((tgrid - 0.5) / 0.1) + .2 * tgrid + .2 ;
          %  Gaussian cdf, "slanted and shifted upwards"
  x = [0.13; 0.20; 0.31; 0.41; 0.72] ;
  o = [0.28; 0.38; 0.64; 0.76; 0.86] ;

%  gpcheta([tgrid,cur],[x,o],1,2) ;
%  gpcheta([tgrid,cur],[x,o],1,0) ;
%  gpcheta([tgrid,cur],[x,o],2) ;
%  gpcheta([tgrid,cur],[x,o],3) ;
          %  used these to test variations
  gpcheta([tgrid,cur],[x,o]) ;
    subplot(2,1,1) ;
    title('Simple Example, from chetaeg1.m') ;


    orient tall ;
  print -dps \matlab\steve\ps\gpc1.ps ;


elseif ipart == 2 ;   %  then work with code decay data
                      %  this part copied from perfass1.m

  %  8 - Module oa/OAxtavt
  modstr = 'oa/OAxtavt' ;

  %  Read in the data from file mr_(8)
  %
  disp(['    Loading MR Data Set for Module number ' ... 
                            num2str(8) ',  ' modstr]) ;
    frmtstr = '%g %g %g %g %g %g %g %g %g %g' ;
    insize = [10, inf] ;
    instr =  ['\matlab\steve\decay\data\mr_' num2str(8)] ;
  fid = fopen(instr,'r') ;
    indata = fscanf(fid,frmtstr,insize) ;
    indata = indata' ;
          %  since data rows are read in as columns
  fclose(fid) ;
  disp(['    Read in ' num2str(size(indata,1)) ' data points']) ;


  %  create vectors for use below
  %
  vdatefd = indata(:,2) ;     %  date of first delta
  vdateld = indata(:,3) ;     %  date of last delta
  vacp = indata(:,9) ;        %  code in A/C/P field

  indata = 0 ;
          %  to save space


  %  work with times
  %
  xstr = 'year (avg 1st, last del)' ;
  vdate = (vdatefd + vdateld) / 2 ;
          %  use date of first delta as default "date"

  mint = min(vdate) ;
  maxt = max(vdate) ;
  delt = (maxt - mint) / 401 ;
          %  spacing on time axis


  %  work with mr's, and bin
  %
  nall = length(vdate) ;
  cflag = (vacp == 2) ;
          %  one where the MR is corrective
  vdatec = vdate(cflag) ;
          %  dates of corrective MRs
  ncorr = length(vdatec) ;
  pflag = (vacp == 4) ;
          %  one where the MR is perfective
  vdatep = vdate(pflag) ;
          %  dates of perfective MRs
  nperf = length(vdatep) ;
    seed = 76528735 ;
    rand('seed',seed) ;
  randind = ceil(nall * rand(ncorr,1)) ;
          %  vector of random indices, among 1,...,nall
  vdater = vdate(randind) ;
          %  vector of randomly chosen MRs, from set of all MRs

  [vcts,tgrid] = gplbinr(vdate,[mint; maxt]) ;
  vctsp = gplbinr(vdatep,[mint; maxt]) ;
  vctsc = gplbinr(vdatec,[mint; maxt]) ;



  %  Now do smoothing to get time damped change rates
  %
  h = - (11/12) / log(.5) ;
          %  based on idea of "Fault potential haf life should be 11 months"

  kge = -h * log(.001 * h) ;
          %  endpoint of grid for kernel evaluation
  kergrid = 0:delt:kge ;
          %  grid for evaluation of kernel
  nkg = length(kergrid) ;
  vker = exp(-kergrid / h) ;
  vker = vker / (sum(vker) * delt) ;
          %  this makes vker work like a density, ie. sum(. * delt) = 1

  tdcr = conv(vcts,vker) ;
  tdcr = tdcr(1:401) ;
  tdcrp = conv(vctsp,vker) ;
  tdcrp = tdcrp(1:401) ;



  %  Next get Proportion of Perfective Maintenance
  %  and plug into gpcheta
  %
  ppm = tdcrp ./ tdcr ;
  gpcheta([tgrid,ppm],[vdatec, vdater]) ;
    subplot(2,1,1) ;
    title(['Code Decay Example, Module ' modstr]) ;


    orient tall ;
  print -dps \matlab\steve\ps\gpc2.ps ;


elseif ipart == 3 ;   %  then work with Allen Glazner's volcano data
                      %  copied from volcano3.m

  ivoldat = 2 ;    %  1 - work with "many.txt"
                   %  2 - work with "one.txt"

  h1 = 0.018 ;    %  bandwidth for volcano data
  h2 = 0.02 ;    %  bandwidth for glaciation data

  %  read in volcano data
  %
  if ivoldat == 1 ;
    infstr = '\matlab\steve\data\many.txt' ;
    outstr = '\matlab\steve\others\manglac' ;
  elseif ivoldat == 2 ;
    infstr = '\matlab\steve\data\one.txt' ;
    outstr = '\matlab\steve\others\oneglac' ;
  end ;

  frmtstr = '%g' ;
  insize = [1, inf] ;
    fid = fopen(infstr,'r') ;
      vdata = fscanf(fid,frmtstr,insize) ;
      vdata = vdata' ;
          %  since data rows are read in as columns
    fclose(fid) ;
  disp(['    Read in ' num2str(size(vdata,1)) ' Volcano data points']) ;


  %  read in glaciation data
  %
  infstr = '\matlab\steve\data\specmap.txt' ;

  frmtstr = '%g %g' ;
  insize = [2, inf] ;
    fid = fopen(infstr,'r') ;
      gdata = fscanf(fid,frmtstr,insize) ;
      gdata = gdata' ;
          %  since data rows are read in as columns
    fclose(fid) ;
  disp(['    Read in ' num2str(size(gdata,1)) ' Glaciation data points']) ;



  left = 0 ;
  right = .782 ;


  %  get uniformly distributed data 
  %
  nv = length(vdata) ;
  seed = 93893749 ;
  rand('seed',seed) ;
  rdata = left + rand(nv,1) * (right - left) ;


  %  smooth glaciation data
  %
  [gsmooth,tgrid] = gpnpr(gdata,h2,[left,right],0,1) ;
            %  0 - binned implementation
            %  1 - truncate points outside range


%  gpcheta([tgrid,gsmooth],[vdata,rdata],3) ;
          %  use the above to fiddle test based on sd's
  gpcheta([tgrid,gsmooth],[vdata,rdata]) ;
    subplot(2,1,1) ;
    title('Volcano Example, from volcano3.m') ;


    orient tall ;
  print -dps \matlab\steve\ps\gpc3.ps ;



elseif ipart == 4 ;   %  then work with Allen Glazner's volcano data
                      %  copied from volcano3.m
                      %  but subtract large bandwidth smooth to detrend

  ivoldat = 2 ;    %  1 - work with "many.txt"
                   %  2 - work with "one.txt"

  h1 = 0.018 ;    %  bandwidth for volcano data
  h2 = 0.02 ;    %  bandwidth for glaciation data

  %  read in volcano data
  %
  if ivoldat == 1 ;
    infstr = '\matlab\steve\data\many.txt' ;
    outstr = '\matlab\steve\others\manglac' ;
  elseif ivoldat == 2 ;
    infstr = '\matlab\steve\data\one.txt' ;
    outstr = '\matlab\steve\others\oneglac' ;
  end ;

  frmtstr = '%g' ;
  insize = [1, inf] ;
    fid = fopen(infstr,'r') ;
      vdata = fscanf(fid,frmtstr,insize) ;
      vdata = vdata' ;
          %  since data rows are read in as columns
    fclose(fid) ;
  disp(['    Read in ' num2str(size(vdata,1)) ' Volcano data points']) ;


  %  read in glaciation data
  %
  infstr = '\matlab\steve\data\specmap.txt' ;

  frmtstr = '%g %g' ;
  insize = [2, inf] ;
    fid = fopen(infstr,'r') ;
      gdata = fscanf(fid,frmtstr,insize) ;
      gdata = gdata' ;
          %  since data rows are read in as columns
    fclose(fid) ;
  disp(['    Read in ' num2str(size(gdata,1)) ' Glaciation data points']) ;



  left = 0 ;
  right = .782 ;


  %  get uniformly distributed data 
  %
  nv = length(vdata) ;
  seed = 93893749 ;
  rand('seed',seed) ;
  rdata = left + rand(nv,1) * (right - left) ;


  %  smooth glaciation data
  %
  [gsmooth,tgrid] = gpnpr(gdata,h2,[left,right],0,1) ;
            %  0 - binned implementation
            %  1 - truncate points outside range

  hfact = 3 ;
  gtrend = gpnpr(gdata,hfact*h2,[left,right],0,1) ;
  gsmooth = gsmooth - gtrend ;



%  gpcheta([tgrid,gsmooth],[vdata,rdata],3) ;
          %  use the above to fiddle test based on sd's
  gpcheta([tgrid,gsmooth],[vdata,rdata]) ;
    subplot(2,1,1) ;
    title('Volcano Example, from volcano3.m') ;


    orient tall ;
  print -dps \matlab\steve\ps\gpc4.ps ;

end ;




