function experiment_design_uncertainty_iterative_process

all_uncertainty = NaN;
all_error = NaN;

t_obs = [1;3];                                  % design of initial experiment
true_para = [14;11.5];                          % true parameter
y_obs = toy_system_evaluate(true_para,t_obs);   % obs data if noise free

% step 1: parameter optimization based on current data
optimization_options = optimoptions(@lsqnonlin,'Algorithm','levenberg-marquardt','MaxFunEvals',1500);
est_para = lsqnonlin(@(x) myfun(x, y_obs, t_obs), [0;0], [], [], optimization_options) % x is estimated_parameter
pred = toy_system_evaluate(est_para,t_obs);

% step 2: compute Jacobian for the current experimental design at estimated parameter, and compute uncertainty
J=[];
J(1,:) = (-t_obs(:).*exp(-est_para(1)*t_obs(:)))';
J(2,:) = (-t_obs(:).*exp(-est_para(2)*t_obs(:)))';
uncertainty = sum(1./(svd([J]).^2))  % evalue of J'J is the same as square of singular value of J

all_uncertainty = [all_uncertainty; sum(1./(svd([J]).^2))];
all_error = [all_error; min(sum((true_para - est_para).^2), sum((true_para - flipud(est_para)).^2))];

for iter=1:5 % let's do 5 rounds of exp design and see what the algorithm will suggest
    
    % step 3: for each possible new experiment compute the Jacobian for current+possible experimental design at estimated parameter and compute uncertainty
    possible_t = 10.^(-5:(log10(5)+5)/1000:log10(5));
    score = zeros(size(possible_t));
    for i=1:length(possible_t)
        new_t = unique([t_obs(:);possible_t(i)]);
        new_J = [];
        new_J(1,:) = (-new_t(:).*exp(-est_para(1)*new_t(:)))';
        new_J(2,:) = (-new_t(:).*exp(-est_para(2)*new_t(:)))';
        score(i) = sum(1./(svd([new_J]).^2));
    end
    [~,ind] = min(score);
    t_new = possible_t(ind);
    best_new_uncertainty = score(ind);

    % visualize the uncertainty across all possible new designs, where we can immediate see a few things that make perfect sense:
    if iter==1
        figure(1); 
        loglog(possible_t, score,'k','linewidth',3)
        xlim([0,5])
        xlabel('Candidate time points')
        ylabel('parameter uncertainty')
        line(t_new*[1,1], ylim,'linewidth',3,'color',[0.6,0.6,0.6],'linestyle',':')
        title(['Iteration ',num2str(iter)])
        legend({'Uncertainty','Best new time point'})
    else
        h = figure(2);  
        set(h,'Position',[64 108 930 640]);
        subplot(2,2,iter-1)
        loglog(possible_t, score,'k','linewidth',3)
        xlim([0,5])
        xlabel('Candidate time points')
        ylabel('parameter uncertainty')
        line(t_new*[1,1], ylim,'linewidth',3,'color',[0.6,0.6,0.6],'linestyle',':')
        title(['(',char(96-1+iter), ') Iteration ',num2str(iter)])
        if iter==2
            legend({'Uncertainty','Best new time point'})
        end
    end
    
    
    % step 4: now that the best new experiment is selected (t=0.0008), let's actually do it, obtain more data
    t_obs = [t_obs;t_new];
    y_obs = toy_system_evaluate(true_para,t_obs);

    % back to step 1: parameter optimization based on current data
    optimization_options = optimoptions(@lsqnonlin,'Algorithm','levenberg-marquardt','MaxFunEvals',1500);
    est_para = lsqnonlin(@(x) myfun(x, y_obs, t_obs), [0;0], [], [], optimization_options) % x is estimated_parameter
    pred = toy_system_evaluate(est_para,t_obs);

    % step 2: compute Jacobian for the current experimental design at estimated parameter, and compute uncertainty
    J=[];
    J(1,:) = (-t_obs(:).*exp(-est_para(1)*t_obs(:)))';
    J(2,:) = (-t_obs(:).*exp(-est_para(2)*t_obs(:)))';
    uncertainty = sum(1./(svd([J]).^2))  % evalue of J'J is the same as square of singular value of J

    all_uncertainty = [all_uncertainty; sum(1./(svd([J]).^2))];
    all_error = [all_error; min(sum((true_para - est_para).^2), sum((true_para - flipud(est_para)).^2))];

end


h = figure(3);
set(h, 'position', [100 100 1100 400])
subplot(1,3,1);
semilogy([0,0,1:length(t_obs)-2],t_obs, '-ok', 'linewidth',3)
title('(a) Time point selected each iteration')
xlim([0,5])
xlabel('Iterations')
grid on

subplot(1,3,2);
plot([0,0,1:length(t_obs)-2],all_error, '-ok', 'linewidth',3)
title('(b) Error in estimated parameter')
xlim([0,5])
xlabel('Iterations')
grid on

subplot(1,3,3);
semilogy([0,0,1:length(t_obs)-2],all_uncertainty, '-ok', 'linewidth',3)
title('(c) Uncertainty after each iteration')
xlim([0,5])
xlabel('Iterations')
grid on
return


function [t_new, possible_t, score] = design_new_exp(est_para, t)

possible_t = 10.^(-5:(log10(5)+5)/1000:log10(5));
score = zeros(size(possible_t));
for i=1:length(possible_t)
    new_t = unique([t(:);possible_t(i)]);
    new_J = [];
    new_J(1,:) = (-new_t(:).*exp(-est_para(1)*new_t(:)))';
    new_J(2,:) = (-new_t(:).*exp(-est_para(2)*new_t(:)))';
    score(i) = sum(1./(svd([new_J]).^2));
end
[~,ind] = min(score);
t_new = possible_t(ind);


% possible_t = 0.00001:0.00001:5;
% score = zeros(size(possible_t));
% J(1,:) = (-t(:).*exp(-est_para(1)*t(:)))';
% J(2,:) = (-t(:).*exp(-est_para(2)*t(:)))';
% for i=1:length(possible_t)
%     J_new = [-possible_t(i)*exp(-est_para(1)*possible_t(i)); ...
%              -possible_t(i)*exp(-est_para(2)*possible_t(i))]; 
%     score(i) = sum(1./(svd([J,J_new]).^2));
% end
% [~,ind] = min(score);
% t_new = possible_t(ind);
% 



function f = myfun(est_para, observed_data, t)
predicted_data = toy_system_evaluate(est_para,t);
f = observed_data - predicted_data;




function obs = toy_system_evaluate(para,t)

if length(para)==2
    para = para(:);
end

N = size(para,2);

for i=1:N
    theta1 = para(1,i);
    theta2 = para(2,i);
    for j=1:length(t)
        obs(j,i) = exp(-theta1*t(j)) + exp(-theta2*t(j));
    end
end

