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

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);
J(1,:) = (-t_obs(:).*exp(-est_para(1)*t_obs(:)))';
J(2,:) = (-t_obs(:).*exp(-est_para(2)*t_obs(:)))';
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
     
    % step 1: compute the profile likelihood
    theta_range = 0:0.01:20;
    likelihood_profile_of_theta_1 = zeros(1,length(theta_range));
    est_para_for_profile_of_theta_1 = zeros(2,length(theta_range));
    optimization_options = optimoptions(@lsqnonlin,'Algorithm','levenberg-marquardt','MaxFunEvals',1500,'Display','off');
    for i=1:length(theta_range)
        theta1 = theta_range(i);
        est_theta2 = lsqnonlin(@(x) my_obj_fun([theta1;x], y_obs, t_obs), [0], [], [], optimization_options); % x is estimated_parameter
        likelihood_profile_of_theta_1(i) = 0.5*sum((y_obs - toy_system_evaluate([theta1;est_theta2],t_obs)).^2);
        est_para_for_profile_of_theta_1(:,i) = [theta1;est_theta2];
    end

    likelihood_profile_of_theta_2 = zeros(1,length(theta_range));
    est_para_for_profile_of_theta_2 = zeros(2,length(theta_range));
    optimization_options = optimoptions(@lsqnonlin,'Algorithm','levenberg-marquardt','MaxFunEvals',1500,'Display','off');
    for i=1:length(theta_range)
        theta2 = theta_range(i);
        est_theta1 = lsqnonlin(@(x) my_obj_fun([x;theta2], y_obs, t_obs), [0], [], [], optimization_options); % x is estimated_parameter
        likelihood_profile_of_theta_2(i) = 0.5*sum((y_obs - toy_system_evaluate([est_theta1;theta2],t_obs)).^2);
        est_para_for_profile_of_theta_2(:,i) = [est_theta1;theta2];
    end

    % step 2: define thresholds for acceptable parameters (very generous definition, different from Steiert et al, 2012)
    acceptable_est_para = [];
    threshold_for_theta_1 = max(exp(-likelihood_profile_of_theta_1))*0.9;
    acceptable_indicator_1 = exp(-likelihood_profile_of_theta_1)>=threshold_for_theta_1;
    acceptable_est_para = [acceptable_est_para, est_para_for_profile_of_theta_1(:,acceptable_indicator_1)];
    threshold_for_theta_2 = max(exp(-likelihood_profile_of_theta_2))*0.9;
    acceptable_indicator_2 = exp(-likelihood_profile_of_theta_2)>=threshold_for_theta_2;
    acceptable_est_para = [acceptable_est_para, est_para_for_profile_of_theta_2(:,acceptable_indicator_2)];

    h = figure(iter); set(h,'position',[100,0,850,800])
    figure(h); subplot(3,2,1); semilogy(theta_range, likelihood_profile_of_theta_1,'k','linewidth',3); line(xlim,[1,1]*(-log(threshold_for_theta_1)),'linewidth',3,'color',[0.6,0.6,0.6],'linestyle',':')
    xlabel('\theta_1'); ylabel('optimized cost function value'); title('(a) loglikelihood profile of \theta_1')
    figure(h); subplot(3,2,2); semilogy(theta_range, likelihood_profile_of_theta_2,'k','linewidth',3); line(xlim,[1,1]*(-log(threshold_for_theta_2)),'linewidth',3,'color',[0.6,0.6,0.6],'linestyle',':')
    xlabel('\theta_2'); ylabel('optimized cost function value'); title('(b) loglikelihood profile of \theta_2')
    figure(h); subplot(3,2,3); % plot(est_para_for_profile_of_theta_1(1,:), est_para_for_profile_of_theta_1(2,:),'g-',est_para_for_profile_of_theta_1(1,find(acceptable_indicator_1==1)), est_para_for_profile_of_theta_1(2,find(acceptable_indicator_1==1)),'b.');
    hold off; plot(est_para_for_profile_of_theta_1(1,:), est_para_for_profile_of_theta_1(2,:),'.', 'color', [0.7,0.7,0.7], 'markersize',15,'linewidth',3);
    hold on;  plot(est_para_for_profile_of_theta_1(1,find(acceptable_indicator_1==1)), est_para_for_profile_of_theta_1(2,find(acceptable_indicator_1==1)),'.', 'color', [0,0,0], 'markersize',15,'linewidth',3);
    xlabel('\theta_1'); ylabel('\theta_2'); title('(c) optimized parameters of loglikelihood profile of \theta_1')
    figure(h); subplot(3,2,4); % plot(est_para_for_profile_of_theta_2(2,:), est_para_for_profile_of_theta_2(1,:),'g-',est_para_for_profile_of_theta_2(2,acceptable_indicator_2==1), est_para_for_profile_of_theta_2(1,acceptable_indicator_2==1),'b.');
    hold off; plot(est_para_for_profile_of_theta_2(2,:), est_para_for_profile_of_theta_2(1,:),'.', 'color', [0.7,0.7,0.7], 'markersize',15,'linewidth',3);
    hold on;  plot(est_para_for_profile_of_theta_2(2,find(acceptable_indicator_2==1)), est_para_for_profile_of_theta_2(1,find(acceptable_indicator_2==1)),'.', 'color', [0,0,0], 'markersize',15,'linewidth',3);
    xlabel('\theta_2'); ylabel('\theta_1'); title('(d) optimized parameters of loglikelihood profile of \theta_2')
    drawnow

    % step 3: For each possible new experiment, simulate the system using all the acceptable parameter settings,
    % compute the range of the simulated predictions, and suggest the experiment with largest range.
    possible_t = 10.^(-5:(log10(5)+5)/1000:log10(5));
    all_simulated_obs=[];
    for i=1:size(acceptable_est_para,2)
        simulated_obs = toy_system_evaluate(acceptable_est_para(:,i),possible_t);
        all_simulated_obs = [all_simulated_obs;simulated_obs'];
    end
    score = max(all_simulated_obs)-min(all_simulated_obs);
    [~,ind] = max(score);
    t_new = possible_t(ind)

    % h = figure(2); set(h,'position',[100,100,850,300])
    figure(h); subplot(3,2,5); semilogx(possible_t, all_simulated_obs,'color',[1,1,1]*0.5,'linewidth',3);
    xlim([min(possible_t), max(possible_t)])
    xlabel('Candidate time point'); ylabel('x_3'); title('(e) Simulated dyanmics using acceptabel parameters')
    figure(h); subplot(3,2,6); semilogx(possible_t, score,'k','linewidth',3); line(t_new*[1,1],ylim,'linewidth',3,'color',[0.6,0.6,0.6],'linestyle',':')
    xlim([min(possible_t), max(possible_t)])
    xlabel('Candidate time point'); ylabel('Range score'); title('(f) Suggest new time point using range score')    
    drawnow
    
    t_obs = [t_obs;t_new];
    y_obs = toy_system_evaluate(true_para,t_obs);   % obs data if noise free
    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);
    J=[];
    J(1,:) = (-t_obs(:).*exp(-est_para(1)*t_obs(:)))';
    J(2,:) = (-t_obs(:).*exp(-est_para(2)*t_obs(:)))';
    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(iter+1);
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 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

