%%% Version: September 3rd, 2017
%%%
%%% This function computes the input-response indices for a given model
%%% based on simulations of the model's equations for the times specified
%%% in tstarValues AND t_ref.
%%%
%%% Sources: Knoechel, Kloft and Huisinga, "Reducing complex 
%%% systems pharmacology models based on a novel input-response index"
%%% J Pharmacokinet Pharmacodyn, 2017 (submitted)
%%%
%%% Authors: Wilhelm Huisinga & Jane Knoechel
%%%

%%% +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
%%% BEGIN: MAIN FUNCTION
%%% +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

function [ir] = compute_input_response_indices_new(model)

I = model.I;

%%% generate figure for interims plotting
figure(3); clf; pause(0.1);


%%% -----------------------------------------------------------------------
%%% determine time span for the solution of the reference trajectory,
%%% the input-perturbed trajectories as well as for the
%%% state-perturbed trajectories starting at the various tstar values.
%%%
tstarValues = model.tstarValues;
t_ref       = model.t_ref';

if (min(tstarValues)<min(t_ref)) || (max(tstarValues)>max(t_ref))
    message(' tstarValues not allowed to be outside the timespan!'); 
    terminate = 1;
    report_error_and_terminate(message,terminate);
end;

%%% all times, including tstar values and reference simulation times
allTimes = unique(sort([t_ref,tstarValues]));


%%% -----------------------------------------------------------------------
%%% compute reference trajectory for all times and the corresponding
%%% reference output
%%%
X0    = model.X0;
u_ref = model.u_ref;

X.init   = X0 + u_ref;

[~,x_ref]  = ode15s(model.ode,allTimes,X.init,model.solver.options);
check_whether_ODE_solver_reported_warnings;

y_ref = model.h(x_ref);


%%% -----------------------------------------------------------------------
%%% compute input-perturbed trajectories starting at 0
%%%
model.Delta_u_per = u_ref*(model.perturbFac-1);

ir_summand = zeros(length(tstarValues),I.nrOfStates);
ir_sum     = zeros(length(tstarValues),I.nrOfStates);

for k = 1:model.nrOfPerturb
    
    tic
    fprintf('\n Perturbation %d/%d - tstar value Nr. ',k,model.nrOfPerturb);
        
    %%% solve ODE with perturbed input at 0. Delta_u_per denotes the
    %%% difference of the perturbed input to the reference input
    %%%
    Delta_u_per = model.Delta_u_per(:,k); 
    
    X.init      = X0 + u_ref + Delta_u_per;
    
    [~,x_per] = ode15s(model.ode,allTimes,X.init,model.solver.options);
    check_whether_ODE_solver_reported_warnings
            
    %%% -----------------------------------------------------------------------
    %%% compute state-perturbed trajectories starting at tstar
    %%%
    for tsNr = length(tstarValues):-1:1
        
        tstar = tstarValues(tsNr);
        
        fprintf('%d-',tsNr);        
        %%% do not simulate trajectories for tstar identical to the end time.
        %%% In this case, the corresponding ir index is zero by definition.
        if tstar == t_ref(end)
            continue;
        end;
               
        %%% determine all times larger or equal to the current tstar value
        tspan = allTimes(allTimes>=tstar)';
        
        for i = 1:I.nrOfStates 
                                
            %%% apply perturbation u_xi of the ith state to the state of the
            %%% reference trajectory at tstar
            u_xi_tstar    = zeros(I.nrOfStates,1);
            u_xi_tstar(i) = x_per(tstar==allTimes,i)-x_ref(tstar==allTimes,i);
            
            %%% if state is approximately unperturbed then skip solving the ODE
            if abs(u_xi_tstar(i)) < model.solver.options.AbsTol
                continue;
            end;
            
            %%% determine state-perturbed initial value and check for 
            %%% non-negativity (relevant in particular for very small
            %%% values)
            X0_refper =  x_ref(tstar==allTimes,:)' + u_xi_tstar;

                        
            %%% solve ODE with perturbed initial conditions at tstar
            %%%
            [~,x_refper] = ode15s(model.ode,tspan,X0_refper,model.solver.options);
            check_whether_ODE_solver_reported_warnings            
                        
            %%% if tspan contains only two elements, i.e., a start and end
            %%% time, the solver will return also all intermediate
            %%% integration steps, rather than just the solution at the two
            %%% time points. Below, we correct for this. 
            if length(tspan) == 2
                x_refper = x_refper([1 end],:);
            end;

            Z = norm(Delta_u_per,2);
                    
            y_refper = model.h(x_refper);
            ir_summand(tsNr,i) = trapz( tspan, ( y_refper - y_ref(allTimes>=tstar,:) ).^2 ) /( Z^2);
            
        end; % for i = 1:I.nrOfStates 
        
    end; % for t = 1:length(tstarRange)-1
        
    %%% update input response index
    ir_sum = ir_sum + ir_summand;
    
    tictoc = toc; minutes = floor(tictoc/60); seconds = round(tictoc-60*minutes);
    fprintf(' (elapsed time: %d min %2d sec)', minutes, seconds);
    
    %%% do some interims plotting
    interims_plot(k,ir_summand,model);
    
end; % for k = 1:model.nrOfPerturb

%%% -----------------------------------------------------------------------
%%% compute input response indices
%%%
ir = sqrt( 1/(allTimes(end)^2*model.nrOfPerturb) * ir_sum );

end
%%% +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
%%% END: MAIN FUNCTION
%%% +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++


%%% +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
%%% BEGIN: LOCAL SUB-ROUTINES
%%% +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

function [] = report_error_and_terminate(message,terminate)
%%% +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

%%% reports an error message and terminates the program

fprintf('\n\n\n');
fprintf(2,'>> %s  << ',message);
fprintf('\n\n\n');
fprintf('   :-)');
fprintf('\n\n\n');
if terminate
    error(' ');
end;

end


%%% +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
function [] = check_whether_ODE_solver_reported_warnings()
%%% +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

%%% check, whether ODE solver reported any warnings, e.g., unmet local
%%% tolerance etc.

[~,warning_ODE_solver] = lastwarn;
if ~isempty(warning_ODE_solver)
    if ~strcmp(warning_ODE_solver,'MATLAB:lang:cannotClearExecutingFunction')
        message = sprintf('ODE solver had problems: %s!',warning_ODE_solver);
        report_error_and_terminate(message,0);
    end;
end;

end


%%% +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
function [] = interims_plot(k,ir_summand,model)
%%% +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
I           = model.I;
output      = model.output;
tstarValues = model.tstarValues;

%%% compute the indices in the right time scale of interest
if isfield(model,'SF_h_to_s')
    SF_h_to_s   = model.SF_h_to_s;
    time        = tstarValues*SF_h_to_s;
    xlabeltext  = 't [s]';
else %%% default time scale is in hours for blood coagulation model
    time        = tstarValues;
    xlabeltext  = 't [h]' ;
end

ir = sqrt(ir_summand); 

%%% graphical output of the input-response indices for selected states
%%%
FigNr = 3;
linestyle = '-';      % positive perturbation
if (model.perturbFac(k)-1)<0
    linestyle = '--'; % negative perturbation
end;
if length(tstarValues)<5 % only stars, if less than 5 tstar values
    linestyle = '*';
end;

figure(FigNr);
for ii=1:length(output.states)
    state = output.states{ii};
    plot(time,ir(:,I.(state)),linestyle,'Color',output.color.(state));
    hold on;
end
title(sprintf('%d/%d perturbations',k,model.nrOfPerturb));
xlabel(xlabeltext); ylabel('ir index');
xlim([0 time(end)]); 
axis square; set(gca,'yscale','lin'); legend(output.legend,'Location','eastoutside');
fett(gcf); pause(0.1);

end

