%%% Version: June 27th, 2017
%%%
%%% call by: Wajima_bloodCoagulation_snakeVenom_MOR_randomOrdering
%%%
%%% Model order reduction for brown snake venom- fibrinogen setting
%%%
%%% Setting: Reduced model for initial decline and recovery dynamics of
%%% fibrinogen based on random ordering of state variables
%%%
%%% by Jane Knoechel 
%%%

clear all; 
addpath('BloodCoagulationModel')

%%% set random seed for repeatability
%%%
rng(20032017)

%%% set variables for figures
%%%

LabelX = 'time [h]';

%%% define brown snake venom-fibrinogen system
%%%
I   = Wajima_bloodCoagulation_indexing;
I.neglectedStates = {}; 

X0  = Wajima_bloodCoagulation_initialValues(I);

%%% assign variables to model structure
%%%

model.I        = I;
model.X0       = X0;

%%% the model can simulate different setting, where however the parameters
%%% are different (Wajima et al). In vivo all reactions are taking place,
%%% in vitro however the production of the factors is not active and thus
%%% assumed to be zero - this is realised by the string set in
%%% model.scenario
model.scenario = 'in_vivo'; % possible scenarios are in_vitro and in_vivo

%%% parameteres of the model were adapted by Gulati to simulate the influence of the
%%% snake venoms on the blood coagulation model. The original parameters
%%% are thus accessible by setting model.author to Wajima and the modified
%%% parameters by setting it to Gulati
model.author   = 'Gulati'; % possible authors are Wajima or Gulati


par = Wajima_bloodCoagulation_parameters(model);

%%% set parameter of intrinsic pathway to zero since otherwise during
%%% numerical solution of perturbed system with ode solver unphysiological
%%% activation of intrinsic pathway (set either v43 or v42 to zero)
%%%

par(I.v43) = 0;

%%% define input state
%%%
dose_snake_venom = 0.0015; % amount in mg
SF_mg_to_nmol = 1e-3/2e5*1e9; 

I.input  = I.AVenom;
u_ref = zeros(size(X0));
u_ref(I.input) = SF_mg_to_nmol*dose_snake_venom; 

%%% define output state
%%%
I.output = I.Fg;
model.h = @(x) x(:,I.output);

%%% assign variables to model structure
%%%
model.par   = par;
model.u_ref = u_ref;

%%% define system of ordinary differential equations
%%%
model.ode   = @(t,X) Wajima_bloodCoagulation_ode(t,X,par,model);


%%% -----------------------------------------------------------------------
%%% solve the model for reference input
%%%
timespan = [0 40];
X_init   = model.X0 + model.u_ref;

model.jac = @(t,X) Wajima_bloodCoagulation_jacobian(t,X,par,model);

options        = odeset('RelTol',1e-9,'AbsTol',1e-12,'NonNegative',1:I.nrOfStates,'Jacobian',model.jac);

[t_ref,x_ref]  = ode15s(model.ode,timespan,X_init,options);

y_ref       = model.h(x_ref); 


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


%%% start of model reduction based on ordering of max of indices
%%%
start = tic;

%%% start MOR for 10 random orderings
%%%
for j=1:10
    
envStates           = [];
negStates           = [];
selectedStates      = I.stateName;
model.I.environStates       = negStates;
model.I.neglectedStates     = envStates;
model.I.selectedStates      = selectedStates;
    
%%% compute random ordering for MOR
%%%
Index(j,:) = randperm(model.I.nrOfStates);


fprintf('\n Starting model order reduction - state Nr. ');
    
for i = 1:I.nrOfStates % iteratively considering each state for elimination
    fprintf('%d-',Index(j,i)); 
    
    currentState = I.stateName(Index(j,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-reducted 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);

end
