% solve whole domain IHCP with Tikhonov regularization
% Use the EXPECTED VALUE of delta-T_RMS in the Morozov Discrepancy
%    Principle to determine the optimal regularization coefficient
% 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 = [  2 1.75 1.5 1.25 1 0.9 0.75 .65 0.5  ]*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);
    deltaT_RMS_vec(ia,1) = sqrt( ET_func( alpha,  X, H, sigma, q_exact )/n ) ;
end
semilogx( alpha_vec, deltaT_RMS_vec(:,1), 'k-' ,...
        [ alpha_vec(1) alpha_vec(end) ], [sigma sigma], 'k--', 'linewidth', 1.0 )
xlabel('\alpha'); ylabel('\itE\rm(\it\DeltaT_{RMS}\rm)')
ymin = 1.5e-3;
ylim( [ ymin 2.3e-3] )
box on

% find the optimal alpha
falpha = @(alpha) ( sqrt(ET_func( alpha, X, H, sigma, q_exact )/n) - ...
                        sigma)^2;
options = optimset('TolX', 1e-8, 'TolFun', 1e-10);
alpha_opt = fminbnd( falpha, 1e-4, 1e-2, options );

DT_min = sqrt(ET_func( alpha_opt, X, H, sigma, q_exact )/n); % should = sigma

% annotations
ax = gca;
x = [alpha_opt alpha_opt];  % data x-coordinates
y = [DT_min ymin];  % data y-coordinates

% Convert to normalized figure coordinates
pt1 = ax.Position(1:2) + [(x(1)-ax.XLim(1))/diff(ax.XLim)*ax.Position(3), ...
                          (y(1)-ax.YLim(1))/diff(ax.YLim)*ax.Position(4)];
pt2 = ax.Position(1:2) + [(x(2)-ax.XLim(1))/diff(ax.XLim)*ax.Position(3), ...
                          (y(2)-ax.YLim(1))/diff(ax.YLim)*ax.Position(4)];

annotation('arrow', [pt1(1) pt2(1)], [pt1(2) pt2(2)]);
text( 1.02*alpha_opt, 0.00155, sprintf('\\alpha_{opt}=%8.4g', alpha_opt) )
text( 0.6*alpha_opt, sigma*1.015, sprintf('\\sigma_{Y} = %8.4f', sigma) )

% generate second plot of qhat and q_exact
F_opt = (X'*X + alpha_opt*H'*H) \ X';
qhat = F_opt*Y(2:n+1);
fprintf('optimal alpha = %10.4g\n',alpha_opt)
figure(98); clf
plot( t(2:n+1)-dt/2, qhat, 'ro', t, 0.6*q_1pulse(t,0,1.2),'k--', 'linewidth', 1.0)
xlabel('time'); ylabel('heat flux');
text(1.1,0.5,sprintf('Tikhonov %s order',order) )
legend('optimal estimate', 'Exact', 'location', 'northwest' )
%      fig=figure(99);
%      print(fig,fullfile(pwd, 'Figure4_14a.tif'),'-dtiff','-r600')
%      
%      fig=figure(98);
%      print(fig,fullfile(pwd, 'Figure4_14b.tif'),'-dtiff','-r600')
     
function [ET ET_bias ET_rand] = ET_func( alpha,  X, H, sigma, q)
    n = size(X,1);
    F = (X'*X + alpha*H'*H) \ X';
    R = (eye(n,n) - X*F);
    ET_bias = (X*q)'*R'*R*X*q;
    ET_rand = sigma^2 * trace( R'*R );
    ET = ET_bias + ET_rand;
end