%%% Version: May 24th, 2017
%%%
%%% transf = create_transf_matrices(model,type)
%%%
%%% Input : model       - with the following necessary fields:
%%%                     I  - Indexing for the state vector and fields:
%%%                         selectedStates - Names of selected States from
%%%                                          original model states
%%%         type        - specifies the type of transformation
%%%
%%% Output : transf     - structure containing the following fields:
%%%                     M - transformation matrix
%%%                     invM - inverse of the transformation matrix M
%%%
%%% This function creates a transformation matrix for the original system
%%% into a reduced system.
%%%
%%% Authors: Jane Knoechel 
%%%
%%% 
%%% NOTE: work in progress

%%% +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
%%% BEGIN: MAIN FUNCTION
%%% +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
function transf = create_transf_matrices(model,type)
I = model.I;

switch type
    case 'Galerkin'
        if isfield(I,'selectedStates')
            check_for_consistency_states_for_reduction(I)
             
            selectedStates = I.selectedStates;
            
            %%% initializing the matrix
            M       = zeros(length(selectedStates),I.nrOfStates);
            %%% compute the corresponding indices for selected states
            I_selectedStates = zeros(length(selectedStates),1);
            for i = 1:length(selectedStates)
                I_selectedStates(i) = I.(selectedStates{i});
            end
            %%% assign ones the selected states such that these are mapped
            %%% to the identity
            
            M(1:length(selectedStates),I_selectedStates) = eye(length(selectedStates));
            invM    = M';
        else
            %%% 
            ErrorMsg = sprintf('To transform the system by Galerkin projection, you need to specify model.Ix.selectedStates');
            report_error_and_terminate(ErrorMsg,1)
        end
    case 'Lumping'
        %%% create lumping matrix
        if isfield(I,'toBeLumpedStates') || isfield(I,'LumpingMatrix')
            %%% two possiblilities for lumping matrix
            if isfield(I,'toBeLumpedStates')
                %%% -----------------------------------------------------------
                %%% 1) construct matrix based on specified strings in
                %%% Ix.toBeLumpedStates (cell matrix with coloumns corresponding to lump)
                check_for_consistency_states_for_reduction(I)
                 
                toBeLumpedStates = I.toBeLumpedStates; 
                
                %%% compute the order of the reduced model by checking how
                %%% many states are specified per lump
                nrOfLumps=size(toBeLumpedStates,2);
                
                tmp_dim = 0;
                for i= 1:nrOfLumps
                    tmp_dim=tmp_dim+length(find(~cellfun('isempty',I.toBeLumpedStates(:,i))));
                end
                
                dim_red_model = I.nrOfStates - tmp_dim + nrOfLumps;
                %%% initializing the matrix
                M       = zeros(dim_red_model,I.nrOfStates);
                
                statesNotLumped = setdiff(I.stateName,toBeLumpedStates(~cellfun('isempty',I.toBeLumpedStates)));
                %%% first step compute identity matrix for all states
                %%% which are not lumped to that end compute indices of
                %%% these states
                I_statesNotLumped = zeros(length(statesNotLumped),1);
                for i = 1:length(statesNotLumped)
                    I_statesNotLumped(i) = I.(statesNotLumped{i});
                end
                M(1:length(statesNotLumped),I_statesNotLumped) = eye(length(statesNotLumped));
                %%% second step compute lumping 
                for i= 1:nrOfLumps
                    % compute index of lumped state in current lump
                    I_currentLump = zeros(length(find(~cellfun('isempty',I.toBeLumpedStates(:,i)))),1);
                    currentLumpStates = toBeLumpedStates(~cellfun('isempty',I.toBeLumpedStates(:,i)));
                    for j = 1:length(currentLumpStates)
                        I_currentLump(j) = I.(currentLumpStates{j});
                    end
                    M(length(statesNotLumped)+i,I_currentLump)=ones(length(currentLumpStates),1);
                end
                
                invM = pinv(M); % pseudo inverse is used to construct inverse of lumping matrix, could be extended in the future
            else
                %%% -----------------------------------------------------------
                %%% 2) for automated lumping algorithm faster to provide
                %%% lumping matrix
                M = I.LumpingMatrix;
                invM = pinv(M); % pseudo inverse is used to construct inverse of lumping matrix, could be extended in the future
            end
        else
            %%% 
            ErrorMsg = sprintf('To transform the system by lumping of states, you need to specify model.Ix.tobeLumpedStates');
            report_error_and_terminate(ErrorMsg,1)
        end
    case 'CombinedLumpingGalerkin'
        %%% note we assume that first the galerkin projection is applied
        %%% and then lumping
        
        if isfield(I,'selectedStates') || isfield(I,'GalerkinMatrix')
            if isfield(I,'selectedStates') && ~isfield(I,'GalerkinMatrix')
                check_for_consistency_states_for_reduction(I)
                
                %%% construct galerking matrix
                selectedStates = I.selectedStates;
                
                %%% initializing the matrix
                G       = zeros(length(selectedStates),I.nrOfStates);
                %%% compute the corresponding indices for selected states
                I_selectedStates = zeros(length(selectedStates),1);
                for i = 1:length(selectedStates)
                    I_selectedStates(i) = I.(selectedStates{i});
                end
                %%% assign ones the selected states such that these are mapped
                %%% to the identity
                
                G(1:length(selectedStates),I_selectedStates) = eye(length(selectedStates));
                invG    = G';
            else
                G       = I.GalerkinMatrix;
                invG    = G';
            end
            %%% construct lumping matrix
            if isfield(I,'toBeLumpedStates') || isfield(I,'LumpingMatrix')
                %%% two possiblilities for lumping matrix
                if isfield(I,'toBeLumpedStates')
                    %%% -----------------------------------------------------------
                    %%% 1) construct matrix based on specified strings in
                    %%% Ix.toBeLumpedStates (cell matrix with coloumns corresponding to lump)
                    toBeLumpedStates = I.toBeLumpedStates;
                    
                    %%% compute the order of the reduced model by checking how
                    %%% many states are specified per lump
                    nrOfLumps=size(toBeLumpedStates,2);
                    
                    tmp_dim = 0;
                    for i= 1:nrOfLumps
                        tmp_dim=tmp_dim+length(find(~cellfun('isempty',I.toBeLumpedStates(:,i))));
                    end
                    
                    dim_red_model = length(selectedStates) - tmp_dim + nrOfLumps;
                    %%% initializing the matrix
                    L       = zeros(dim_red_model,I.nrOfStates);
                    
                    statesNotLumped = setdiff(selectedStates,toBeLumpedStates(~cellfun('isempty',I.toBeLumpedStates)));
                    %%% first step compute identity matrix for all states
                    %%% which are not lumped to that end compute indices of
                    %%% these states
                    I_statesNotLumped = zeros(length(statesNotLumped),1);
                    for i = 1:length(statesNotLumped)
                        I_statesNotLumped(i) = strcmp(selectedStates,statesNotLumped(i));
                    end
                    L(1:length(statesNotLumped),I_statesNotLumped) = eye(length(statesNotLumped));
                    %%% second step compute lumping
                    for i= 1:nrOfLumps
                        % compute index of lumped state in current lump
                        I_currentLump       = zeros(length(find(~cellfun('isempty',I.toBeLumpedStates(:,i)))),1);
                        currentLumpStates   = toBeLumpedStates(~cellfun('isempty',I.toBeLumpedStates(:,i)));
                        for j = 1:length(currentLumpStates)
                            I_currentLump(j) = I.(currentLumpStates{j});
                        end
                        L(length(statesNotLumped)+i,I_currentLump)=ones(length(currentLumpStates),1);
                    end
                    
                    invL = pinv(L); % pseudo inverse is used to construct inverse of lumping matrix, could be extended in the future
                else
                    %%% -----------------------------------------------------------
                    %%% 2) for automated lumping algorithm faster to provide
                    %%% lumping matrix
                    L = I.LumpingMatrix;
                    invL = pinv(L); % pseudo inverse is used to construct inverse of lumping matrix, could be extended in the future
                end
                %%% assign output matrices
                M = L*G; % x_red = M*x = L*G*x;
                invM = invG*invL;
            else
                %%%
                ErrorMsg = sprintf('To transform the system by lumping of states, you need to specify model.I.tobeLumpedStates or I.LumpingMatrix');
                report_error_and_terminate(ErrorMsg,1)
            end
        else
            %%% 
            ErrorMsg = sprintf('To transform the system by Galerkin projection, you need to specify model.I.selectedStates or I.GalerkinMatrix');
            report_error_and_terminate(ErrorMsg,1)
        end
    otherwise
        %%% 
        ErrorMsg = sprintf('The Type of transformation is not yet included or was misspelled, check spelling or chose different type.');
        report_error_and_terminate(ErrorMsg,1)
end

%%% assign variables to output structure (transf)
%%%
transf.M    = M ;
transf.invM = invM;

end
%%% +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
%%% END: MAIN FUNCTION
%%% +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++


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

function [] = report_error_and_terminate(message,terminate)
%%% +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

%%% reports an error message and terminates the program

if terminate
    error(message);
end

end

function [] = check_for_consistency_states_for_reduction(I)
%%% +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

%%% check for consistency such that no mispecification are included in the
%%% reduced model

selectedStates      = I.selectedStates;

if isfield(I,'environStates')
    envStates           = I.environStates;
    %%% envStates need to be a subset of the selectedStates
    %%%
    
    consistencyEnvStates = setdiff(envStates,selectedStates);
    
    if ~isempty(consistencyEnvStates)
        ErrorMsg = sprintf('Environmental states need to be a subset of the selected states otherwise will be neglected.');
        report_error_and_terminate(ErrorMsg,1)
    end
end
if isfield(I,'toBeLumpedStates')
    toBeLumpedStates    = I.toBeLumpedStates;
    %%% toBeLumpedStates need to be a subset of the selected States
    %%%
    
    consistencyLumpStates = setdiff(toBeLumpedStates,selectedStates);
    
    if ~isempty(consistencyLumpStates)
        ErrorMsg = sprintf('The to be lumped states need to be a subset of the selected states.');
        report_error_and_terminate(ErrorMsg,1)
    end
    if isfield(Ix,'environStates')
        %%% the toBeLumpedStates should not be included in the envStates
        %%%
        consistencyLumpStatesEnv = setdiff(toBeLumpedStates,envStates);
        if length(consistencyLumpStatesEnv)<length(toBeLumpedStates)
            ErrorMsg = sprintf('The to be lumped states cannot be a subset of the environmental states.');
            report_error_and_terminate(ErrorMsg,1)
        end
    end
end

end