% solve whole domain IHCP with Tikhonov regularization
% assume dimensionless variables
A = 15;
order = 'zero';
x = [ 1. ];   % measurement locations
t = [ -0.24: 0.06: 1.8 ]';
T = [ 0.000000 0.000000 0.000000 0.000000 0.000000 0.000007 ...
      0.000374 0.002171 0.006323 0.013381 0.023656 0.037319 ...
      0.054465 0.075145 0.099389 0.127201 0.157880 0.189293 ...
      0.219594 0.247679 0.272930 0.295005 0.313715 0.328954 ...
      0.340665 0.348823 0.353763 0.356545 0.358089 0.358943 ...
      0.359415 0.359677 0.359821 0.359901 0.359945 ]';

dt = t(2) - t(1);
rng(129);
noise = randn( size(t) );
sigma_pct = 0.005;

sigma = sigma_pct*max(T);
Y = T + sigma*noise;
% can't have negative times
tzero = t - t(1);
phi = fdX22B10T0( x, tzero, A );
dphi = diff( phi );
n = length( dphi );  % exclude first data point (I.C.)
X = zeros( n, n );
for ic = 1:n
    X( ic:n, ic ) = dphi(1:n-ic+1);
end
cond(X); 
if ( strcmp(order,'first') )
    H = -eye(n,n);
    for in = 1:n-1
        H(in,in+1)=1;
    end
    H(1,1:3)= [-3 4 -1]/2;
    H(n,n-2:n)= [1 -4 3]/2;
    alpha_vec = [ 5e-3 4e-3 3e-3 2e-3 1e-3 5e-4  ];
elseif (strcmp(order, 'second') )
    H = -2*eye(n,n);
    for in =2:n-1
        H(in,in+1)=1;
        H(in,in-1)=1;
    end
    H(1,1:4)= [2 -5 4 -1];
    H(n,n-3:n)= [-1 4 -5 2];
    alpha_vec = [ 5e-3 1e-3 5e-4  ];
else
    H = eye(n,n);
%    alpha_vec = [ 5 4 3 2.75 2.5 2 1.75 1.5 1.25 1 0.9 0.75 .65 0.5  ]*1e-3;
    alpha_vec = [1.7 1.6 1.5 1.4 1.300 1.2 1.1  1 0.9 0.75 .65 0.6  0.46 ...
                 0.2  0.1 0.05 0.045 0.035 0.03 0.025  0.02 0.015 ... 
                 0.0125 0.01 0.009 0.007 0.005 0.004 0.003 0.0025 0.002 ...
                 0.00175 0.0015 0.00125 0.001 0.0009]*1e-3;
end

q_exact = 0.6*q_1pulse( t(2:n+1)-dt/2, 0, 1.2);
figure(99); clf;
hold on
legstr = [ "ro", "b+", "md", "ks", "gv" ];
for ia = 1:length(alpha_vec)
    alpha = alpha_vec(ia);
    [ xLcurve(ia,1), yLcurve(ia,1) ] = Lcurve_func( alpha,  X, H, Y(2:end) );
end
loglog( xLcurve, yLcurve, 'k.' ,...
           'linewidth', 1.0 )
ylabel('log(||Y-Xq||)'); xlabel('log(||q||)')
iq = find( alpha_vec == 0.0017, 1);
plot(xLcurve(iq), yLcurve(iq),'ro','linewidth',1.0)
text(1.07*xLcurve(iq), yLcurve(iq), sprintf('\\alpha=%8.2e min(R_q)',alpha_vec(iq)))
iT = find( alpha_vec == 0.0012, 1);
plot(xLcurve(iT), yLcurve(iT),'ro','linewidth',1.0)
text(1.07*xLcurve(iT), 1.01*yLcurve(iT), sprintf('\\alpha=%8.2e MDP',alpha_vec(iT)))
iTE = find( alpha_vec <= 1.31e-3, 1);
plot(xLcurve(iTE), yLcurve(iTE),'ro','linewidth',1.0)
text(1.07*xLcurve(iTE), 1.01*yLcurve(iTE), sprintf('\\alpha=%8.2e MDP E.V.',alpha_vec(iTE)))
iG=find( alpha_vec == 0.00046, 1);
plot(xLcurve(iG), yLcurve(iG),'ro','linewidth',1.0)
text(1.07*xLcurve(iG), 1.01*yLcurve(iG), sprintf('\\alpha=%8.2e GCV',alpha_vec(iG)))
iL=find( alpha_vec == 0.0000125, 1);
plot(xLcurve(iL), yLcurve(iL),'ro','linewidth',1.0)
text(1.07*xLcurve(iL), 1.07*yLcurve(iL), sprintf('\\alpha=%8.2e',alpha_vec(iL)))

box on

%      fig=figure(99);
%      print(fig,fullfile(pwd, 'Figure4_16.tif'),'-dtiff','-r600')
     
     
function [x y] = Lcurve_func( alpha,  X, H, Y)
    n = size(X,1);
    F = (X'*X + alpha*H'*H) \ X';
    qhat = F*Y;
    That = X*qhat;
    x = qhat'*qhat;
    y = (Y-That)'*(Y-That);
end