% pj, 09/03/2004: NOTE: end year is hardcoded to 2010 below


close all; clear; 

site_list = input('site list (in single quotes):   ');
comp = input('neu direction component (in single quotes):   ');
factor = input('input edit criteria factor:   ');
[sites] = textread(site_list,'%4c','commentstyle','shell');
if comp == 'n'
    neu_txt = 'north';
    column = 1;
elseif comp == 'e'
    neu_txt='east';
    column = 2;
elseif comp == 'u'
    neu_txt='up';
    column = 3;
end 

array_name = site_list(1:4);

% get all times for all sites to be included in stack
times=[];
goodsites=[];
count=0;
nsites = length(sites(:,1));
for s = 1:nsites    
  site = sites(s,:);
  file=[site '_resids_unf.dat'];


  % pj, 09/02/2004: if site is skipped in cleaning, skip it here
  % code by pfang
  [fexist,tmp]=size(dir(file));
  if (fexist)
     eval(['load ' file]);
     array=[site '_resids_unf'];
     data=[eval(array)];
     times = [times; data(:,1)];
     count = count+1;
     goodsites(count,:) = sites(s,:);
  % 
  end

end

time=unique(times);
ndays = length(time);

% fill the residual matrices, change units to mm
%posmatrix = zeros(ndays,nsites);
%sigmatrix = ones(ndays,nsites).*9999999999;


% pj, 09/02/2004
ngoodsites=length(goodsites(:,1));


posmatrix = zeros(ndays,ngoodsites);
sigmatrix = ones(ndays,ngoodsites).*9999999999;

for s=1:ngoodsites
%for s = 1:nsites 
  clear data row

    %site = sites(s,:);
    site = goodsites(s,:);

    array=[site '_resids_unf'];
    data=[eval(array)];
    data(:,2:7) = data(:,2:7).*1000;

    rms_n = std(data(:,2));
    istart = find(data(:,1) > 1990.0) ;
    istart = istart(1);
    istop = find(data(:,1) < 2011.0);
    istop = istop(length(istop));

    [c,idata,itime] = intersect(data(:,1),time);
    posmatrix(itime,s) = data(idata,1+column);
    sigmatrix(itime,s) = data(idata,4+column);
    
end

% ITERATION 1

% compute number of sites per day; apply stack only when at least 3 sites
%npos1 = round(nsites-sum(sigmatrix./9999999999,2));
npos1 = round(ngoodsites-sum(sigmatrix./9999999999,2));

inostack = find(npos1 < 3);
posmatrix(inostack,:) =[];
sigmatrix(inostack,:) =[]; 
npos1(inostack,:) =[]; 
time_stack1 = time;
time_stack1(inostack,:) =[]; 

% compute the stack as the weighted mean
weight = 1./sigmatrix.^2;
mn = sum(posmatrix.*weight,2)./sum(weight,2);

%for i = 1:nsites
for i = 1:ngoodsites

  mn_matrix(:,i) = mn;
end

fltmatrix = posmatrix - mn_matrix;
stat1 = sum(fltmatrix.*weight,2);
stat2 = sum(weight,2);
stat3 = sum(fltmatrix.^2.*weight,2);
nrms = sqrt( stat3./(npos1-1) );
wrms = sqrt( npos1./stat2 ).*nrms;
sig_mn = wrms./sqrt(npos1);

% clean posmatrix of outliers
edit_criterion = factor .* sig_mn;
for r = 1:size(time_stack1)
   iedits = find(abs(fltmatrix(r,:)) >= edit_criterion(r));
   pos = posmatrix(r,:);
   sig = sigmatrix(r,:);
   pos(iedits) = 0;
   sig(iedits) = 9999999999;
   posmatrix(r,:) = pos;
   sigmatrix(r,:) = sig;
end

% ITERATION 2

%% compute number of sites per day; apply stack only when at least 3 sites
clear inostack weight mn mn_matrix fltmatrix stat1 stat2 stat3 nrms wrms sig_mn

%npos = round(nsites-sum(sigmatrix./9999999999,2));
npos = round(ngoodsites-sum(sigmatrix./9999999999,2));

ncleaned = npos1-npos;
inostack = find(npos < 3);
posmatrix(inostack,:) =[];
sigmatrix(inostack,:) =[]; 
npos(inostack,:) =[]; 
time_stack = time_stack1; 
time_stack(inostack,:) =[]; 

% compute the stack as the weighted mean
weight = 1./sigmatrix.^2;
mn = sum(posmatrix.*weight,2)./sum(weight,2);

%for i = 1:nsites
for i = 1:ngoodsites
  mn_matrix(:,i) = mn;
end


fltmatrix = posmatrix - mn_matrix;
stat1 = sum(fltmatrix.*weight,2);
stat2 = sum(weight,2);
stat3 = sum(fltmatrix.^2.*weight,2);
nrms = sqrt( stat3./(npos-1) );
wrms = sqrt( npos./stat2 ).*nrms;
sig_mn = wrms./sqrt(npos);

% save stack in m, but then go back to mm
stack(:,1) = time_stack;
stack(:,2) = mn./1000;
stack(:,3) = sig_mn./1000;
stack(:,4) = npos;
eval(['save stack_' array_name '_' comp '.dat stack -ascii'])
stack(:,2) = mn;
stack(:,3) = sig_mn;

%compute slope of the stack
p  = polyfit(stack(:,1),stack(:,2),1);
slope = p(1);

figure
orient tall
subplot(4,1,1)
plot(stack(:,1), stack(:,2),'.')
xlabel 'years'
ylabel 'stack (mm)'
axis([1991 max(time) min(stack(:,2))*1.5 max(stack(:,2))*1.5])
text(max(time)-.1*(max(time)-min(time)),max(stack(:,2)).*1.25,neu_txt)
txt=sprintf('slope = %6.2f mm/yr',slope);
text(min(time)+.1*(max(time)-min(time)),max(stack(:,2)).*1.25,txt)
eval(['print -dpsc stack_' array_name '_' comp '.ps'])
subplot(4,1,2)
plot(stack(:,1), wrms,'.')
axis([1991 max(time) 0 max(wrms)*1.2])
text(max(time)-.1*(max(time)-min(time)),max(wrms).*1.1,neu_txt)
xlabel 'years'
ylabel 'postfit scatter (mm)'
subplot(4,1,3)
plot(stack(:,1), stack(:,4),'.')
axis([1991 max(time) 2.5 max(npos1)+.5])
xlabel 'years'
ylabel '# sites in stack'
subplot(4,1,4)
plot(time_stack1, ncleaned,'.')
axis([1991 max(time) -0.5 max(ncleaned)+1])
xlabel 'years'
ylabel '# sites cleaned'
txt=sprintf('edit criteria = %6.0f x sig',factor);
text(1992, max(ncleaned)+.5,txt)
eval(['print -dpsc stack_' array_name '_' comp '.ps'])

% subtract the stack from each timeseries
%for s = 1:nsites    
for s = 1:ngoodsites    
  clear data c idata istack data_filt

  %site = sites(s,:);
  site = goodsites(s,:);

  array=[site '_resids_unf'];
  data=[eval(array)];
  data(:,2:7) = data(:,2:7).*1000;
  [c,idata,istack] = intersect(data(:,1),time_stack);
  data_filt = [data(:,1) data(:,1+column) data(:,4+column)];
  data_filt(idata,2) = data(idata,1+column) - stack(istack,2);
  data_filt(idata,3) = sqrt(data(idata,4+column).^2 + stack(istack,3).^2);
  [wrms_data,chi2]=wstat2(data(:,1+column),data(:,4+column),1);
  [wrms_filt,chi2_filt]=wstat2(data_filt(:,2),data_filt(:,3),1);

  figure
  subplot(2,1,1)
  plot(data(:,1), data(:,1+column), '.')
  eval(['title ' site ':' neu_txt])
  ylabel 'un-filtered (mm)'
  axis([min(data(:,1)) max(data(:,1)) min(data(:,1+column))*1.3 max(data(:,1+column))*1.3])
  txt=sprintf('wrms = %6.2f mm',wrms_data);
  text(min(data(:,1))+.1*(max(data(:,1))-min(data(:,1))),max(data(:,1+column)).*1.1,txt)
  xlabel 'years'
  subplot(2,1,2)
  plot(data(:,1), data_filt(:,2), '.')
  axis([min(data(:,1)) max(data(:,1)) min(data(:,1+column))*1.3 max(data(:,1+column))*1.3])
  txt=sprintf('wrms = %6.2f mm',wrms_filt);
  text(min(data(:,1))+.1*(max(data(:,1))-min(data(:,1))),max(data(:,1+column)).*1.1,txt)
  xlabel 'years'
  ylabel 'filtered (mm)'
  eval(['print -dpsc ' site '_' comp '_unf_flt.ps'])
end
