addpath '/home/jscales/matlab/regularization'  % path to P-CH's reg package
% This package is freey available from netlib.org.


% parmeters of synthetic noise
exactvar = 1;  % a scaling parameter for the noise
corlen = 0;    % 0 corresponds to white noise
dist = 'U';    % U is uniform, N is normal

ncols = 200;   % size of model used in reconstruction.  can be changed
               % arbitrarily.  The "true" model is picewise continuous
               % with a low-velocity zone.  The true data were simulated 
	       % using 1000 layers, as if piecewise continuous.
	          

nrows = 100;   % number of data
ztop = 0;      % physical dimensions of the model
zbottom = 40;
deltaz = (zbottom - ztop)/ncols;
firstrec = ztop + deltaz/2;
lastrec = zbottom ;
deltarec = (lastrec - firstrec)/nrows;
rec = firstrec:deltarec:lastrec-deltarec;

% create the velocity model.  The exact model has a linear increase
% in velocity and a low velocity zone in the middle (irrespective
% of parameterization).

vmax = 4;
vmin = .5;
dv = (vmax - vmin)/ncols
% exactvelmodel(1:ncols) = 2;
exactvelmodel=vmin:dv:vmax-dv;
mid = floor(ncols/2);
deltam = floor(.1 * ncols);
exactvelmodel(mid-deltam:mid+deltam) = 1;
exactslowmodel = 1 ./ exactvelmodel;

% this is just for zero-offset, straight ray tomography
jacobian =  makejacobian(nrows,ncols,deltaz,rec);  

% starting model
backgroundvelmodel(1:ncols) = 2;
backgroundslowmodel = 1./backgroundvelmodel;

backgroundtt = jacobian * backgroundslowmodel';

% creat syntheic travel times.  samplevar is the sample variance of the 
% synthetic noise, as opposed to the variance of the probability model.
[times,samplevar] = maketimes(jacobian,exactslowmodel,1000,100,dist,0, ...
			      exactvar,corlen);     
deltatt = times - backgroundtt;

%%%%%%%%%%%%%end of setup%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%% first object lesson is how awful the pseudo-inverse model is using
%% all the singular values

updatedslowmodel = backgroundslowmodel + ( pinv(jacobian) * deltatt)' ;
updatedvelmodel = 1./updatedslowmodel;
predictedtimes = jacobian * updatedslowmodel';
set(gca,'FontSize',16)
plot(updatedvelmodel,'k-','LineWidth',2)
title('pseudo-inverse model')
xlabel('depth layer')
ylabel('velocity (km/s)')
hold on
axis([0 200 0 5])
plot(exactvelmodel,'r:','LineWidth',2)

figure
set(gca,'FontSize',16)
plot(times,'r-','LineWidth',2)
hold on
xlabel('receiver number')
ylabel('time (ms)')
plot(predictedtimes, 'k-','LineWidth',2)
legend('data','predicted for pseudoinv',4)


%% so even though the model fits the data well, it's lousy.
%% we can do better by regularizing the problem a bit.  The problem
%% is always to know how much regularization to use.  Here
%% we use the L-curve approach.  This is a popular method that doesn't
%% require knowledge of the errors in the data.  It jointly optimizes the
%% data variance reduction and the model (semi)norm.

% a compact singular value decomposition
[U,s,V] = csvd(jacobian);

figure
set(gca,'FontSize',16);
[reg_corner,rho,eta,reg_param] = l_curve(U,s,deltatt);
loglog(rho,eta,'LineWidth',2)
set(gca,'FontSize',16)
xlabel('residual norm')
ylabel('model norm')
text(1,1.5,strcat('L-curve:  corner = ',num2str(reg_corner)), 'FontSize',16)

% Now we compute the Tikhonov model associated with the regularization
% parameter given by the L-curve

lambda =  reg_corner;  
tikupdate = tikhonov(U,s,V,deltatt,lambda) ;
tikhonovmodel = 1./(tikupdate + backgroundslowmodel');
figure
set(gca,'FontSize',16);plot(tikhonovmodel,'LineWidth',2)
hold on
plot(exactvelmodel,':r','LineWidth',2)
axis([0 ncols 0 6])
%title('Tikhonov + L-curve')
xlabel('Depth layer')
ylabel('Velocity (km/s)')
legend('Tikhonov','true model',2)

figure
% now plot the preicted times
plot(times,'r-','LineWidth',2)
hold on
set(gca,'FontSize',16);
plot(jacobian * (backgroundslowmodel + tikupdate')','k-','LineWidth',1)
legend('data','predicted for Tikhonov',4)
xlabel('Receiver number')
ylabel('Time (ms)')
hold off
% Here is the first estimate of sigma
residual = times - jacobian * (backgroundslowmodel + tikupdate')';
sigmaest = sqrt(var(residual))


%% here we look explicitly at the trade off of model norm and data fit.
%% finding the simplest model that fits the data is called
%% occam's inversion.  if you can't achieve a chi-square of 1,
%% then loop over svalues used and find an optimal value based on where
%% the curve asymptotes.

figure
for i=1:20
  svalsused = i;
  % truncated svd
  tsvdslowupdate = backgroundslowmodel' + tsvd(U,s,V,deltatt,svalsused);
  tsvdvelupdate = 1./tsvdslowupdate;
  tsvdresidual = times - jacobian * tsvdslowupdate;
  tsvdsigmas(i) = sqrt(var(tsvdresidual));
end
set(gca,'FontSize',16);
semilogy(tsvdsigmas,'r*','LineWidth',2)
hold on
semilogy(tsvdsigmas,'r-','LineWidth',2)
ylabel('residual \sigma')
xlabel('number of singular values')
linex = [0 20];
liney = [samplevar samplevar];
line(linex,liney)
text(5,1,strcat('true noise variance = ',num2str(samplevar)),'FontSize',16)

figure
svalused = 5;
% truncated svd
tsvdslowupdate = backgroundslowmodel' + tsvd(U,s,V,deltatt,svalused);
tsvdvelupdate = 1./tsvdslowupdate;
set(gca,'FontSize',16);
plot(tsvdvelupdate,'bo','LineWidth',2)
%title(strcat('truncated svd.  number of svals used = ', num2str(svalused)))
tsvdresidual = times - jacobian * tsvdslowupdate;
sqrt(var(tsvdresidual))
hold on
plot(exactvelmodel,'k-','LineWidth',2)
legend('TSVD model','true model')
xlabel('Depth layer')
ylabel('velocity (km/s')
hold off
updatedsigma = sqrt(var(tsvdresidual))

figure
set(gca,'FontSize',16)
% plot(times,'r.','LineWidth',2)
hold on
title('predicted versus measured times')
xlabel('receiver number')
ylabel('time (ms)')
predictedtimes = jacobian*tsvdslowupdate;
plot(predictedtimes, 'r.','LineWidth',2)
errorbars = samplevar * ones(1,length(times));
errorbar(1:100,times,errorbars)
legend('predicted for truncated SVD','data',4)
axis([0 100 0 30])