%%% Version: May 29th, 2017
%%%
%%% call by: Wajima_bloodCoagulation_PTtest_SnakeVenom_MOR_1h
%%%
%%% Model order reduction for brown snake venom- fibrinogen setting
%%%
%%% Setting: Reduced model for the first hour after envenomation 
%%%         (initial decline of fibrinogen)
%%%
%%% by Jane Knoechel 
%%%

clear all; 
addpath('BloodCoagulationModel')

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

%%% set variables for figures
%%%

LabelX = 't [h]';

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

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

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

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

stateName = model.I.stateName;

sortedStatesBasedIndices = stateName(Index_sortedMax(model.I.nrOfStates:-1:1));
LegendStates             = Wajima_bloodCoagulation_legendLabel(sortedStatesBasedIndices);

figure(FigNr)
plot(1:19,ir_sortedMax(model.I.nrOfStates:-1:model.I.nrOfStates-18),'*')
set(gca,'XTick',1:19)
set(gca,'XTickLabel',LegendStates(1:19))
set(gca,'yscale','log')
xticklabel_rotate([],45,[],'Fontsize',24)
set(gca,'FontSize',24);
ylim([1e-5 1e4])
ylabel('||ir||_\infty')

FigNr = 2;
figure(FigNr)
plot(1:19,ir_sortedMax(36:-1:18),'*')
set(gca,'XTick',1:19)
set(gca,'XTickLabel',LegendStates(20:38))
set(gca,'yscale','log')
xticklabel_rotate([],45,[],'Fontsize',24)
set(gca,'FontSize',24);
ylim([1e-5 1e4]);
ylabel('||ir||_\infty')

%%% 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                  = 1e-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-1 % iteratively considering each state for elimination
    fprintf('%d-',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.Ix.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.Ix.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 - done : Total elapsed time: %d h %d min %2d sec',hours, minutes, seconds);
fprintf('\n\n 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')

%%% visualisation of reduced model dynamics
%%%

redmodel = model;

redmodel.I.environStates   = envStates;
redmodel.I.neglectedStates = negStates;
redmodel.I.selectedStates  = selectedStates;

%%% create transformation matrix
%%%
redmodel.transf = create_transf_matrices(model,'Galerkin');

%%% 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,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);

FigNr = 2;

figure(FigNr)
%%% plot reference output
plot(t_ref,y_ref,'Color',redmodel.output.color.Fg)
hold on;
%%% plot the reduced systems output
plot(t_ref,y_red,'--','Color',redmodel.output.color.Fg)
xlabel(LabelX); ylabel('concentration [nM]');
xlim([-0 model.timespan(end)]); 
axis square; legend('Fg','Fg_{red}','Location','eastoutside');
fett(gcf);