%%% Version: June 5th, 2017
%%%
%%% [L_out,invL_out] = AutomatedLumping_DokoumetzidisAarons(model,TOL)
%%%
%%% This function uses an automated lumping algorithm to produce for the
%%% system specified in the structure model an reduced system.
%%%
%%% Input: model - structure with the following necessary fields:
%%%                    -
%%%
%%% Output: L_out - lumping matrix
%%%
%%% Sources: A. Dokoumetzidis and L.Aarons "Proper lumping in systems 
%%%          biology models" (2009)
%%%
%%% Authors: Jane Knoechel 
%%%

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

function [L_out,invL_out] = AutomatedLumping_DokoumetzidisAarons(model,TOL)
%%% set required variables for automated lumping
%%%
X_init  = model.X_init;
odeFcn  = model.odeFcn;
t_ref   = model.t_ref;
options = model.solver.options;

%%% get state space dimension of model system
%%%
if isfield(model.I,'selectedStates')
    if isfield(model.I,'environStates')
        n = length(setdiff(model.I.selectedStates,model.I.environStates));
    else
        n = length(model.I.selectedStates);
    end
else
    n = model.I.nrOfStates;
end
%%% initialise variables for lumping algorithm
%%%
L_old                = eye(n);
noOfPossCombinations = factorial(n)/(factorial(2)*factorial(n-2));
err_rel              = zeros(1,noOfPossCombinations);

L = L_old;

if isfield(model.I,'environStates')
    %%% extend lumping matrix by the identity for the environmental
    %%% states, assumption that environmental states are the final
    %%% entries of the state vector
    %%%
    
    L = [L_old zeros(size(L_old,1),length(model.I.environStates)); zeros(length(model.I.environStates),n) eye(length(model.I.environStates))];
end

%%% compute reference system
%%% 
model.I.LumpingMatrix  = L;
model.transf           = create_transf_matrices(model,'CombinedLumpingGalerkin');
X_init_red             = model.transf.M*X_init;

[~,x_temp]           = ode15s(@(t,X) odeFcn(t,X,model),t_ref,X_init_red,options);
y_ref                = model.h(x_temp*model.transf.invM');

%%% start of lumping algorithm: iteratively lumping in pairs and using the assumption
%%% that optimal lumped model is found by not unlumping states once lumped
%%%
for i=1:n-1
    Combinations            = nchoosek(1:(n+1-i),2);
    L_pre                   = eye(n-i,n-i);
    L_poss                  = zeros(n-i,n-i+1,size(Combinations,1));
    L                       = zeros(n-i,n,size(Combinations,1));
    
    for k=1:size(Combinations,1)
        %%% compute possible lumping matrices for ith step, the complete
        %%% lumping matrix is then the combination of the new lumping
        %%% matrix and the old: L = L_ith*L_(i-1)th
        %%%
        L_poss(:,:,k)           = [L_pre(:,1:Combinations(k,2)-1),...
                                    L_pre(:,Combinations(k,1)),...
                                    L_pre(:,Combinations(k,2):end)];
        L(:,:,k)                = L_poss(:,:,k)*L_old;
        
        %%% if model is specified by dynamical and environmental states
        %%% (constant), then environmental states will not be considered for 
        %%% lumping and thus need to be accounted for seperately 
        %%%
        L_tmp = L(:,:,k);
        if isfield(model.I,'environStates') 
            %%% extend lumping matrix by the identity for the environmental
            %%% states, assumption that environmental states are the final
            %%% entries of the state vector
            %%%
            
            L_tmp = [L(:,:,k) zeros(size(L(:,:,k),1),length(model.I.environStates)); zeros(length(model.I.environStates),n) eye(length(model.I.environStates))];
        end
        model.I.LumpingMatrix = L_tmp ;
        
        model.transf = create_transf_matrices(model,'CombinedLumpingGalerkin');
        X_init_red   = model.transf.M*X_init;
        %%% solve reduced system
        %%%
        [~,x_tmp]    = ode15s(@(t,X) odeFcn(t,X,model),t_ref,X_init_red,options);
        y_red        = model.h(x_tmp*model.transf.invM');
        %%% compute relative L2 error
        %%% Note: since for specific lumping matrix ode system might run
        %%%       into numerical issues use try and catch to set error for
        %%%       these cases
        %%%
        try
            err_rel(k)        = relativeErrorL2(t_ref,y_ref,y_red);
        catch
            err_rel(k)        = Inf;
        end
    end
    %%% compute minimal error and corresponding index in ith combination
    %%% possibilities
    %%%
    [~,Ind] = min(err_rel);
    
    %%% only if error is below the user defined threshold, ith step is
    %%% accepted otherwise the algorithm is stoped
    %%%
    if min(err_rel)<TOL
        L_old = L(:,:,Ind);
    else
        break;
    end
    clear err_rel
end
%%% assign variables to output 
%%%
if isfield(model.I,'environStates')
    %%% extend lumping matrix by the identity for the environmental
    %%% states, assumption that environmental states are the final
    %%% entries of the state vector
    %%%
    
    L_out = [L_old zeros(size(L_old,1),length(model.I.environStates)); zeros(length(model.I.environStates),n) eye(length(model.I.environStates))];
    invL_out = pinv(L_out);
else
    L_out       = L_old;
    invL_out    = pinv(L_out);
end
end
%%% +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
%%% END: MAIN FUNCTION
%%% +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++


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

function Error = relativeErrorL2(t,X,Y)

Error = sqrt(trapz(t,(X-Y).^2)/trapz(t,X.^2));

end