%VSP velocity inversion example for Chapter 5
%(originally suggested by Gary Pavlis, Indiana University).
clear
rand('state',0);
randn('state',0);

%load particular realization of the data vector used in book
load vsp

%Number of model parameters
k=1;
M=k*N;

%System matrix
G=zeros(N,M);
for i=1:N
  G(i,1:k*i) = ones(1,k*i)*dd;
end

%tril(ones(N,N))*dd;

mls=(G\t);

[U,s,V]=svd(G);

%higher order Tikh regularization
[L1,W] = get_l(N,1);
[U1,sm1,X1,V1]=cgsvd(G,L1);
disp('GCV for first order:')
pause 
figure(2)
bookfonts;
[reg_corner1,gdvfun,reg_param]=gcv(U1,sm1,t,'Tikh');
[m1,rho,eta]=tikhonov(U1,sm1,X1,t,reg_corner1);
xlim([1 1e4]);
%print -deps2 c5fgcv1curve.eps

[L2,W] = get_l(N,2);
[U2,sm2,X2,V2]=cgsvd(G,L2);
disp('GVC for second order:')
pause
figure(3)
bookfonts;
[reg_corner2,gdvfun,reg_param]=gcv(U2,sm2,t,'Tikh');
[m2,rho,eta]=tikhonov(U2,sm2,X2,t,reg_corner2);
xlim([1 1e4]);
%print -deps2 c5fgcv2curve.eps

for i=1:N
  j=(i-1)*20+10;
  struesamp(i)=mean(strue(j-9:j+10));
end


%plot/compare model results
disp('GCV first order and true model:')
pause
figure(4)
bookfonts;
plotconst(m1*1000,0,maxdepth)
hold on
plot(depth,strue*1000,'k--')
hold off
ylabel('Slowness (s/km)')
xlabel('Depth (m)')
disp(['1st order reg corner is:  ',num2str(reg_corner1)]);
%print -deps2 c5fgcv1model.eps
disp('GCV second order and true model:')
pause

figure(5)
bookfonts;
plotconst(m2*1000,0,maxdepth)
hold on
plot(depth,strue*1000,'k--')
hold off
ylabel('Slowness (s/km)')
xlabel('Depth (m)')
disp(['2nd order reg corner is:  ',num2str(reg_corner2)]);
%print -deps2 c5fgcv2model.eps

dinterp=10:20:990;
strueresamp=spline(depth,strue,dinterp)';

disp(['1st order model 2-norm misfit is ',num2str(norm(m1-strueresamp))]);
disp(['2nd order model 2-norm misfit is ',num2str(norm(m2-strueresamp))]);




