%%% Version: May 29th, 2017
%%%
%%% call by: Wajima_bloodCoagulation_snakeVenom_MOR_40h
%%%
%%% Model order reduction for brown snake venom- fibrinogen setting
%%%
%%% Setting: Reduced model for initial decline and recovery dynamics of fibrinogen
%%%
%%% by Jane Knoechel 
%%%

clear all; 
addpath('BloodCoagulationModel')

%%% load the input-response indices and variables for model specification
%%%
file    = dir('Wajima_bloodCoagulation_snakeVenom_40h.mat');
load(file.name)

%%% set variables for figures
%%%

LabelX = 'time [h]';

%%% compute maximal value of indices and sort
%%%

ir_max                          = max(ir,[],1);
[ir_sortedMax,Index_sortedMax]  = sort(ir_max);

%%% visualise obtained ordered indices
%%%
FigNr =10;

figure(FigNr)
semilogy(model.I.nrOfStates:-1:1,ir_sortedMax,'*')

disp([num2cell((1:model.I.nrOfStates)') model.I.stateName(Index_sortedMax(model.I.nrOfStates:-1:1))' num2cell(ir_sortedMax(model.I.nrOfStates:-1:1))'])


%%% set reference solution
%%%
t_ref       = model.t_ref;
x_ref       = model.x_ref;
X_init      = model.X_init;
options     = model.solver.options;
y_ref       = model.h(x_ref); 
I           = model.I;

%%% Inicialise variables necessary for model reduction
%%%

TOL                  = 2e-1; % model reduction relative error tolerance

options.Jacobian    = [];
envStates           = [];
negStates           = [];
selectedStates      = I.stateName;

odeFcn               = @(t,X,model)  Wajima_bloodCoagulation_ode(t,X,model.par,model);


fprintf('\n Starting model order reduction - state Nr. ');
    
%%% start of model reduction based on ordering of max of indices
%%%
start = tic;

model.I.selectedStates =selectedStates;

for i = 1:I.nrOfStates % iteratively considering each state for elimination
    fprintf('%s-',I.stateName{Index_sortedMax(i)});%Index_sortedMax(i)); 
    
    currentState = I.stateName(Index_sortedMax(i));
    
    %%% reset options for ode solver
    options.NonNegative = 1:length(model.I.selectedStates);
    
    %%% first consider the state to belong to the environment
    %%%
    model.I.environStates  = [envStates,currentState];
    
    %%% create transformation matrix
    %%%
    model.transf = create_transf_matrices(model,'Galerkin'); 
    
    %%% set initial value for reduced model
    %%%
    X_init_red = model.transf.M*X_init;
        
    %%% compute solution of red model with environmental states
    %%%
    
    [~,x_red_tmp] = ode15s(@(t,X) odeFcn(t,X,model),t_ref,X_init_red,options);
    
    %%%  map reduced solution back to full original system to compute the output 
    %%%
    x_red_env = x_red_tmp*model.transf.invM';
    y_red_env = model.h(x_red_env);
    
    %%% reset environment states for calculation of system where this state
    %%% is neglected
    %%%
    model.I.environStates = envStates;
         
    %%% second consider the state to be negligible 
    %%%
    
    model.I.neglectedStates = [negStates,currentState];
    model.I.selectedStates  = setdiff(selectedStates,currentState);
    
    %%% reset options for ode solver
    options.NonNegative = 1:length(model.I.selectedStates);
    
    %%% create transformation matrix
    %%%
    model.transf = create_transf_matrices(model,'Galerkin'); 
    
    %%% set initial value for reduced model
    %%%
    X_init_red = model.transf.M*X_init;
    
    %%% compute solution of reduced model with neglected states
    %%%

    [~,x_red_tmp] = ode15s(@(t,X) odeFcn(t,X,model),t_ref,X_init_red,options);
    %%%  map reduced solution back to full original system to compute the output 
    %%%
    x_red_neg = x_red_tmp*model.transf.invM';
    y_red_neg = model.h(x_red_neg);
    
    %%% reset neglected states for next calculation of system where this state
    %%% is set to be environment
    %%%
    model.I.neglectedStates = negStates;
    model.I.selectedStates  = selectedStates;
    
    %%% compute error between the reference and reduced solutions
    %%%
    
    err_rel_env = sqrt(trapz(t_ref,(y_ref-y_red_env).^2)/trapz(t_ref,y_ref.^2));
    err_rel_neg = sqrt(trapz(t_ref,(y_ref-y_red_neg).^2)/trapz(t_ref,y_ref.^2));
    
    rel_sol_error = sqrt(trapz(t_ref,(y_red_env-y_red_neg).^2)/trapz(t_ref,y_red_env.^2));
    
    %%% consider smallest error to chose elimination process
    %%%
    if err_rel_env < TOL && err_rel_neg < TOL % both errors are below threshold
        if err_rel_env > err_rel_neg % state is to be neglected
            negStates       = [negStates,currentState];
            selectedStates  = setdiff(selectedStates,currentState);
            % reset resulting transfer matrix
            model.I.selectedStates  = selectedStates;
        elseif err_rel_env < err_rel_neg && rel_sol_error > options.RelTol
            envStates   = [envStates,currentState];
        else
            %%% in the case that the error in the neglected case is only
            %%% larger due to tolerance of ode solver the state is still neglected 
            negStates       = [negStates,currentState];
            selectedStates  = setdiff(selectedStates,currentState);
             % reset resulting transfer matrix
            model.I.selectedStates  = selectedStates;
        end
    elseif err_rel_env < TOL || err_rel_neg < TOL % only one error is below threshold
        if err_rel_env < TOL
            envStates = [envStates,currentState];
        else
            negStates = [negStates,currentState];
            selectedStates = setdiff(selectedStates,currentState);
             % reset resulting transfer matrix
            model.I.selectedStates  = selectedStates;
        end
    end
end
tictoc = toc(start); hours=floor(tictoc/(60*60)); minutes = floor(tictoc/60)-60*hours; seconds = round(tictoc-60*minutes-60*60*hours);

dynStates = setdiff(selectedStates,envStates);
%%% produce format strings to print selected and environmental states
%%%
dynString = repmat('%s, ',1,length(dynStates));
envString = repmat('%s, ',1,length(envStates));

fprintf('\n\n Model order reduction based on indices - done : Total elapsed time: %d h %d min %2d sec',hours, minutes, seconds);
fprintf('\n\n Elimination-reduced model consists of %d dynamical states ',length(dynStates));
fprintf(['\n these are : ' dynString(1:end-2) ' '],dynStates{:});
fprintf('\n and %d environmental states ',length(envStates));
fprintf(['\n these are : ' envString(1:end-2) ' '],envStates{:});
fprintf('\n')

%%% reduce the obtained elimination-reduced model further by using lumping
%%% by Dokoumetzidis and Aarons
%%% 
redmodel = model;

redmodel.I.environStates   = envStates;
redmodel.I.neglectedStates = negStates;
redmodel.I.selectedStates  = [dynStates,envStates];

lumpedredmodel = redmodel;
%%% get transformation matrix for first reduction step such that not
%%% constructed in each step of lumping algorithm
%%%
redmodel.transf = create_transf_matrices(redmodel,'Galerkin');

%%% visualise elimination-reduced model output
%%%

%%% reset options for ode solver
redmodel.solver.options.NonNegative = 1:length(redmodel.I.selectedStates);
redmodel.solver.options.Jacobian    = [];
     
%%% set initial value for reduced model
%%%
X_init_red = redmodel.transf.M*X_init;

%%% compute solution of reduced model with neglected states
%%%
[~,x_tmp] = ode15s(@(t,X) odeFcn(t,X,redmodel),t_ref,X_init_red,redmodel.solver.options);

%%%  map reduced solution back to full original system to compute the output
%%%
x_red      = x_tmp*redmodel.transf.invM';
y_red      = redmodel.h(x_red);

%%% Compute relative error between original model and red model for two
%%% time windows : pre and post nadir
%%%

[y_refmin,Index_min]=min(y_ref);

t_refNadir = t_ref(Index_min);

err_rel_preNadir=sqrt(trapz(t_ref(t_ref<t_refNadir),(y_ref(t_ref<t_refNadir)-y_red(t_ref<t_refNadir)).^2)/trapz(t_ref(t_ref<t_refNadir),y_ref(t_ref<t_refNadir).^2));
err_rel_postNadir=sqrt(trapz(t_ref(t_ref>t_refNadir),(y_ref(t_ref>t_refNadir)-y_red(t_ref>t_refNadir)).^2)/trapz(t_ref(t_ref>t_refNadir),y_ref(t_ref>t_refNadir).^2));

err_rel=sqrt(trapz(t_ref,(y_ref-y_red).^2)/trapz(t_ref,y_ref.^2));

%%% apply lumping
%%%
lumpedredmodel.I.GalerkinMatrix = redmodel.transf.M;
lumpedredmodel.odeFcn            = odeFcn;

%%% reset options for ode solver
lumpedredmodel.solver.options.NonNegative = [];
lumpedredmodel.solver.options.Jacobian    = [];

%%% apply lumping algorithm of Dokoumetzidis and Aarons
%%% 
start = tic;

[L,invL] = AutomatedLumping_DokoumetzidisAarons(lumpedredmodel,TOL);

tictoc = toc(start); hours=floor(tictoc/(60*60)); minutes = floor(tictoc/60)-60*hours; seconds = round(tictoc-60*minutes-60*60*hours);
fprintf('\n\n Further model reduction by lumping - done : Total elapsed time: %d h %d min %2d sec',hours, minutes, seconds);
fprintf('\n\n Reduced model consists of %d states ',size(L,1)-length(envStates));
%%% print lumped states
%%%
GraficLumpedModelIndices = [3,1,5,2,4];
for i=1:size(L,1)-length(envStates)
    LumpedStates(i,1:sum(L(i,:)))=redmodel.I.selectedStates(logical(L(i,:)));
    LumpedString = repmat('%s, ',1,sum(L(i,:)));
    fprintf(['\nThe %d lump consists of the following states in the elim.-red. model : ' LumpedString(1:end-2) ' '],GraficLumpedModelIndices(i),LumpedStates{i,1:sum(L(i,:))});
end
fprintf('\n')


lumpedredmodel.I.LumpingMatrix = L;

lumpedredmodel.transf = create_transf_matrices(lumpedredmodel,'CombinedLumpingGalerkin');
%%% visualisation of reduced model dynamics
%%%


%%% set initial value for reduced model
%%%
X_init_red = lumpedredmodel.transf.M*X_init;

%%% compute solution of reduced model with neglected states
%%%
[~,x_tmp] = ode15s(@(t,X) odeFcn(t,X,lumpedredmodel),t_ref,X_init_red,lumpedredmodel.solver.options);

%%%  map reduced solution back to full original system to compute the output
%%%
x_lump_red      = x_tmp*lumpedredmodel.transf.invM';
y_lump_red      = lumpedredmodel.h(x_lump_red);


SF_nm_to_um = 1/1000;

FigNr = 5;

figure(FigNr)
%%% plot reference output
plot([-1;t_ref],[y_ref(1); y_ref]*SF_nm_to_um,'Color',lumpedredmodel.output.color.Fg)
hold on;
%%% plot the elimination-reduced systems output
plot([-1;t_ref],[y_red(1); y_red]*SF_nm_to_um,'--','Color',lumpedredmodel.output.color.Fg)
hold on;
%%% plot the lumped reduced systems output
plot([-1;t_ref],[y_lump_red(1); y_lump_red]*SF_nm_to_um,'-.','Color',lumpedredmodel.output.color.Fg)
xlabel(LabelX); ylabel('concentration [\muM]');
xlim([-0.1 1]);ylim([0 9])
axis square; legend('Fg','Fg_{red}','Fg_{lumped-red}','Location','eastoutside');
fett(gcf);

%%% inset
%%%
axes('Position',[0.311360123647604 0.423314794215795 0.35 0.35])
box on
plot([-1;t_ref],[y_ref(1); y_ref]*SF_nm_to_um,'Color',lumpedredmodel.output.color.Fg)
hold on;
%%% plot the elimination-reduced systems output
plot([-1;t_ref],[y_red(1); y_red]*SF_nm_to_um,'--','Color',lumpedredmodel.output.color.Fg)
hold on;
%%% plot the lumped reduced systems output
plot([-1;t_ref],[y_lump_red(1); y_lump_red]*SF_nm_to_um,'-.','Color',lumpedredmodel.output.color.Fg)
set(gca,'XLim',[-1 model.timespan(end)])
set(gca,'YLim',[0 9])
axis square
fett(gcf)