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

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

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

%%% 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
%%%

Tolerance                  = [1e-3,1e-2,1e-1]; % model reduction relative error tolerance
for j=1:3
    TOL = Tolerance(j);
    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('%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')
    
    Order_redModel(j) =length(dynStates);
    Order_env_redModel(j) =length(envStates);
    
end

FigNr = 2;

figure(FigNr)
%%% barplot order against tolerance
bar(Tolerance,Order_redModel,'basevalue',1e-10)
set(gca,'XTick',Tolerance)
set(gca,'FontSize',24);
ylabel('order of reduced model')