%% Simple Ammonia Loop - Evaluation of the neural network fit as function 
%% of the norm of the weights (individually calculated).
% This program evaluates the correlation between the norm of the
% significant weights and the surrogate model fit. To this end, it samples
% data (or loads data calculated in "ProgSampAna_woFit.m"). Subsequently,
% it performs PLS regression at the different sampling steps and fits 10
% neural networks to the data. PLS regression is applied to all dependent 
% variables individually. The program "ProgSampAna_wFit_all.m" is complementary 
% and applieds PLS regression simultaenously.
%
% Written by Julian Straus, Trondheim 09.12.2018
%
%  The components in this system are defined as 
%   x(1) - Hydrogen
%   x(2) - Nitrogen
%   x(3) - Ammonia
%   x(4) - Argon
%   x(5) - Methane

% Clearing of the workspace and the memory, closing of all open
% windows/figures
clc, clear variables, close all

% Start of the timer for the overall function file
time.in.overall = tic;

% Definition of the number of reactors in the CSTR cascade. Currently, the
% initialization requires 10 CSTRs for each reactor. If this should be
% changed, then the initialization in InitSamp has to be adjusted.
Reac1_n  = 10;
Reac2_n  = 10;

chem = 3;
% Definition of the inlet values for the pipe
def.P.FPipe = [2000; 27e5; 10;...                           Feed Stream Q [mol/s], p [Pa], and T [C]
%                [65; 24; 10]/100];...                        Feed Stream mole fractions 3 Chemicals
               [70; 23; 5; 1; 1]/100];...                   Feed Stream mole fractions 5 Chemicals
%                [65; 21; 2; 1; 1; 2; 3; 5]/100];...          Feed Stream mole fractions 8 Chemicals
%                ones(chem,1)/chem];...                       Feed Stream mole fractions 3/5/8 Chemicals

% Definition of the inlet values for the makeup section
def.M.FMake = [1809.17; 27.32e5; 7.495;...                  Feed Stream Q [mol/s], p [Pa], and T [C]
              [73.02; 24.28; 1.62; 0.49; 0.59]/100];...     Feed Stream mole fractions [-]
def.M.HEx1  = [9274.61; 10e5; -10;...                       Water Stream HEx1 Q [mol/s], p [Pa], and T [C]
               ones(5,1)];...                               Water Stream HEx1 mole fractions [-]
def.M.FMix  = [6561.38; 127.50e5; 25.82;...              	Mixer Stream Q [mol/s], p [Pa], and T [C]
              [64.77; 20.25; 2.66; 6.98; 5.34]/100];...     Mixer Stream Q mole fractions [-]

% Definition of the inlet values for the reaction section
def.R.FReact = [8274.29; 140e5; 51.62;...                   Feed Stream Q [mol/s], p [Pa], and T [C]
                0.6658; 0.2113; 0.0243; 0.0557; 0.0429];... Feed Stream mole fractions [-]
def.R.Val1   = 832.4;...                                    Split Valve Q [mol/s]
def.R.Temp   = [379.1; 410.9];...                           Temperature Reactor Inlets [C]
    
% Calculation of the inlet values for the separation section
def.S.FSep = [7539.40; 131.2e5; 53.73;...                   Feed Stream Q [mol/s], p [Pa], and T [C]
              0.5728; 0.1791; 0.1385; 0.0619; 0.0477];...   Feed Stream Mole fractions [-]
def.S.HEx6 = [11184.22; 3.828e5; 10;...                     Water Stream HEx6 Q [mol/s], p [Pa], and T [C]
              ones(5,1)];...                                Water Stream HEx6 mole fractions

% Calculation of version histories
Ver.Lid   = 'n';
Ver.model = 'Rea';      % 'Pip', 'Rea', 'Sep'
Ver.Ind   = 'MolFrac';  % 'MolFlow', 'MolFrac', or 'PartPress'
Ver.add   = true;       % true or false (Adaptive sampling)
Ver.dep   = true;       % true or false (Incorporation of dependency)
Ver.diff  = true;       % true or false (Difference of absolute values)
Ver.save  = true;       % This works only for the pipe model atm for saving
Ver.load  = true;       % Load of the training points

% Loading of the needed parameters
ParSamp

% Calculation of the sampling points
opt.val  = 'y';
opt.type = 'LHS';
opt.Np   = 5e3;

% Calculation of the initial conditions and the anonymous functions
[initPipe,     ConAnPipe,...
 initMakeUp,   ConAnMakeUp,...
 initSep,      ConAnSep,...
 initReact,    ConAnReact]...
= InitSamp(par,scl,con,ther,def,itg,dirmain,Ver);

% Calculation of the options used in each solver 'iter-detailed'
options = optimset('Algorithm','trust-region-dogleg','Display','off',...
                   'Jacobian','on','MaxIter',1e3,'MaxFunEvals',1e5,...
                   'TolCon',1e-8,'TolFun',1e-8);

% Addition of the  Path direction of the unit operations
UO = [dirmain,'\Util4. Unit Operations\Simulation'];
addpath(UO)

% Definition of a local safe directory
dirsave = 'LocalSave\';

% Addition of an additional extra line break
fprintf('\n\n')


%% Prerequesite a) Definition of the chosen Ver.model

switch Ver.model
case 'Pip'
    % Definition of the anonymous function
    SampAn  = @(initSep,U)  SampPipe(initSep,U,options,par.P,scl.P,con,ther,def.P,itg);
    
    % Definition of the varied vector
    U.ini.tot = [def.P.FPipe(2:3)' def.P.FPipe(4:end)'*def.P.FPipe(1)];
    
    % Definition of the calculated values and indices used
    ncal = [];
    ind.in = 1:itg.k-1;
    ind.out  = itg.k+1:2*itg.k;
          
    % Definition of the range for plotting
    lbmub = zeros(3,numel(U.ini.tot));
    lbmub(2,:) = U.ini.tot;
%     lbmub(1,:) = U.ini.tot.*[1   1  0.7 0.6*ones(1,itg.j-1)] ...
%                           - [3e5 5 zeros(1,itg.j)];
%     lbmub(3,:) = U.ini.tot.*[1   1  1.3 1.4*ones(1,itg.j-1)] ...
%                           + [3e5 5 zeros(1,itg.j)];
    
%     lbmub(1,:) = U.ini.tot.*[1   1  0.8*ones(1,itg.j)] ...
%                           - [4e5 10 zeros(1,itg.j)];
%     lbmub(3,:) = U.ini.tot.*[1   1  1.2*ones(1,itg.j)] ...
%                           + [4e5 10 zeros(1,itg.j)];
    
    lbmub(1,:) = U.ini.tot.*[1   1  0.5 0.5*ones(1,itg.j-1)] ...
                          - [4e5 10 zeros(1,itg.j)];
    lbmub(3,:) = U.ini.tot.*[1   1  1.5 1.5*ones(1,itg.j-1)] ...
                          + [4e5 10 zeros(1,itg.j)];
    
    % Definition of the initial values
    init = initPipe;
    
    % Definition of the calculation function
    CalAn = @(U,x) CalPipe(U,x,ind,par.P,scl.P,con,opt,Ver,On1j);
    
case 'Rea'
    % Definition of the anonymous function the the varied vector
    SampAn  = @(initReact,U)  SampReact(initReact,U,options,par.R,scl.R,con,ther,def.R,itg);
    U.ini.tot = [def.R.FReact(2:3)' def.R.FReact(4:8)'*def.R.FReact(1) par.R{3}.ratio(1:2)' par.R{17}.Tref];

    % Definition of the calculated values and indices used
    ncal = [];
    ind.in  = 1:10;
    ind.out = 9:16;
    ind.mH  = [2:3 5:8];
    ind.HEx5in  = 74:81;
    ind.HEx5out = 1:8;
    
    % Definition of the range for plotting
    lbmub = zeros(3,numel(U.ini.tot));
    lbmub(1,:) = U.ini.tot.*[1 1 0.8 0.8 0.8 0.8 0.8 1 1 1] - [10e5 20 zeros(1,5) 5 20 15];
    lbmub(2,:) = U.ini.tot;
    lbmub(3,:) = U.ini.tot.*[1 1 1.2 1.2 1.2 1.2 1.2 1 1 1] + [10e5 20 zeros(1,5) 5 20 15];
    
    % Definition of the initial values
    init = initReact;
    
    if Ver.dep
        % Definition of the dependency
        opt.dep = [3 4 def.R.FReact(5)/def.R.FReact(4) .1];...
    end

    % Definition of the calculation function
    CalAn = @(U,x) CalReact(U,x,ind,par.R,scl.R,con,opt,Ver,On1j);
    
case 'Sep'
    % Definition of the anonymous function
    SampAn  = @(initSep,U)  SampSep(initSep,U,options,par.S,scl.S,con,ther,def.S,itg);
    
    % Definition of the varied vector
    U.ini.tot = [def.S.FSep(2:3)' def.S.FSep(4:8)'*def.S.FSep(1) ...
                  def.S.HEx6(1:3)' par.S{3}.ratio(2) ...
                  par.S{9}.Tref par.S{11}.Tref par.S{15}.Tref];
    
    % Definition of the calculated values and indices used
    ncal = [];
    ind.in = 1:14;
    ind.out  = [65:72 170:177 27];
    ind.mH  = [2:3 5:8];
    ind.HEx = [90 114 146];
    ind.HEx8in  = 33:40;
    ind.HEx8out = 90:97;
    ind.HEx9in  = 90:97;
    ind.HEx9out = 114:121;
    ind.HEx10in  = 138:145;
    ind.HEx10out = 146:153;
          
    % Definition of the range for plotting
%     lbmub = zeros(3,numel(U.ini.tot));
%     lbmub(1,:) = U.ini.tot.*[1   1  0.7 0.7 0.7 0.7 0.7 0.7 1   1 1   1 1 1] ...
%                           - [10e5 15 zeros(1,6)               2e5 10 10 5 5 5];
%     lbmub(2,:) = U.ini.tot;
%     lbmub(3,:) = U.ini.tot.*[1   1  1.3 1.3 1.3 1.3 1.3 1.3 1   1 1   1 1 1] ...
%                           + [10e5 15 zeros(1,6)               2e5 10 10 5 5 5];
    
    lbmub(1,:) = U.ini.tot.*[1    1  0.85 0.8 0.8 0.8 0.8 0.8 1 1 1  1 1 1] ...
                          - [10e5 25 zeros(1,6)               0 5 40 4 4 8];
    lbmub(2,:) = U.ini.tot;
    lbmub(3,:) = U.ini.tot.*[1    1  1.15 1.2 1.2 1.2 1.2 1.2 1 1 1  1 1 1] ...
                          + [10e5 25 zeros(1,6)               0 5 40 4 4 8];
                      
    % Definition of the initial values
    init = initSep;
    
    if Ver.dep
        % Definition of the dependency
        opt.dep = [3 4 def.S.FSep(5)/def.S.FSep(4) .1];...
    end
    
    % Definition of the calculation function
    CalAn = @(U,x) CalSep(U,x,ind,par.S,scl.S,con,opt,Ver,On1j);
    
end

%% Prerequesite b) Initialization of the system

% Calculation of the initial function
[x.ini.tot,~,~,exitflag.ini,time.ini] = fsamp(SampAn,lbmub(2,:),init,ncal);
        
% Calculation of the necessary variables
[U.ini,x.ini] = CalAn(U.ini,x.ini);

% Calculation of the number of varied variables
n.in  = numel(ind.in);
n.out = numel(ind.out);
n.red = numel(x.ini.red);

% Printing of the results
disp('The initialization has been performed.')

%% Part a) Calculation of the fitting points

if Ver.load
switch Ver.model
    case 'Pip'
        load('Localsave\Paper\Pip5_20180826_AddSamp_diff_wo_5e3.mat','U','x')
    case 'Rea'
        load('Localsave\Paper\Rea_20180902_AddSamp_diff_wo_5e3.mat','U','x')
    case 'Sep'
        load('Localsave\Paper\Sep_20180227_AddSamp_diff_wo_5e3.mat','U','x')
end
else
    for k = 1:n.PLS
        % Definition of thenumber of sampling points
        if k == 1
            opt.Np = Np.a;

            % Call of the grid definition function
            Ufull = GridDef(lbmub,ind.in,opt);
        else
            opt.Np = Np.b;
            % Call of the grid definition function
            Ufull = GridDef_add(lbmub,ind.in,opt,U.tra.un);
        end
        % Readjustment of the output values
        U.inc.un  = Ufull.grid;
        U.inc.tot = Ufull.out;

        % Combination of the matrices
        U.tra.un        = [U.tra.un; U.inc.un];
        U.tra.tot       = [U.tra.tot; U.inc.tot];

        if rem(Np.a+Np.b*(k-1),100) == 0
            fprintf('%4i points were defined.\n',Np.a+Np.b*(k-1))
        end
    end
    
    % Sampling of the data
    time.in.train = tic;
    [x.tra.tot,cal.tra,b.tra,exitflag.tra,time.tra,U.tra.out] = fsamp(SampAn,U.tra.tot,x.ini.tot',ncal);
    time.train = toc(time.in.train);
end
% Calculation of the necessary variables
[U.tra,x.tra,OnN1.tra,nsamp.tra] = CalAn(U.tra,x.tra);

% Printing of the results
disp('The training space has been sampled.')


%% Part b) Load of the validation points

% Load of the sampled validation points
switch Ver.model
    case 'Pip'
        load([dirsave Ver.model num2str(numel(U.ini.tot)-2) '_20180224_1e4RND_dep.mat'],'Uval','xval')
    case 'Rea'
        load([dirsave 'Rea_20180227_1e4RND_dep.mat'],'Uval','xval')
    case 'Sep'
        load([dirsave 'Sep_20180227_1e4RND_dep.mat'],'Uval','xval')
end
% Printing of the results
disp('The validation space has been loaded.')

% Readjustment of the output values
U.val.tot = Uval;
x.val.tot = xval;

% Calculation of the necessary variables
[U.val,x.val,OnN1.val,nsamp.val] = CalAn(U.val,x.val);

% Definition of the labels for the table and the scaling factor for the
% error
switch Ver.model
    case 'Sep'
        label.RM = {'Pressure G'; 'Temperature G'; 'Pressure L'; 'Temperature L'; 'Temperature W'; ...
                    'Gas H2'; 'Gas N2'; 'Gas NH3'; 'Gas Ar'; 'Gas CH4'; 'Duty 8'; 'Duty 9'; 'Duty 10'};
        scl.err = [ones(nsamp.val,1)*[10^-5 1 10^-5 1 1] U.val.tot(:,3:7) ones(nsamp.val,3)*10^-3];
%         scl.err = 1./x.val.red*100;

    case 'Rea'
        label.RM = {'Pressure'; 'Temperature'; 'Extent of reaction'; 'Duty'};
        label.Leg = {'p_{out}','T_{out}','EoR','W'};
        scl.err = [10^-5 1 1 10^-3];
%         scl.err = 1./x.val.red*100;

    case 'Pip'
        label.RM = {'Pressure'};
        label.Leg = {'p_{out}'};
        scl.err = 10^-5;
%         scl.err = 1./U.val.samp(:,1)*100;

end


%% Part c) Evaluation of the partial least square performance

% Definition of some options
opt.comptresh = 0.02;
opt.Np = 2000;

% Calculation of several values
Np.a = 25;
Np.c = 5;
Np.d = [5 25 100];
Np.NNt = [100 1000];
Np.tot = Np.a:Np.c:opt.Np;
Np.totNN = [Np.a:Np.d(1):Np.NNt(1) ...
            Np.NNt(1)+Np.d(2):Np.d(2):Np.NNt(2) ...
            Np.NNt(2)+Np.d(3):Np.d(3):opt.Np];

n.PLS  = (opt.Np-Np.a)/Np.c+1;
n.NN   = numel(Np.totNN);
n.norm = 5*Np.c;

ind.NN = [1:Np.d(1)/Np.c:(Np.NNt(1)-Np.a)/Np.c+1 ...
          (Np.NNt(1)-Np.a+Np.d(2))/Np.c+1:Np.d(2)/Np.c:(Np.NNt(2)-Np.a)/Np.c+1 ...
          (Np.NNt(2)-Np.a+Np.d(3))/Np.c+1:Np.d(3)/Np.c:n.PLS];

normtyp = 'fro';

% Selection of the chosen data points
switch Ver.Ind
case 'MolFlow'
    U.ini.sel = U.ini.samp;
    U.tra.sel = U.tra.samp;
    U.val.sel = U.val.samp;
case 'MolFrac'
    U.ini.sel = U.ini.Molfrac;
    U.tra.sel = U.tra.Molfrac;
    U.val.sel = U.val.Molfrac;
case 'PartPress'
    U.ini.sel = U.ini.PartPress;
    U.tra.sel = U.tra.PartPress;
    U.val.sel = U.val.PartPress;
end

% Calculation of the indices for training, validation and testing set and
% rearranging of the training set in the case of a complete grid
trainR = 75/100;
valR   = 15/100;
testR  = 15/100;
if opt.val == 'n'
    CG = 2^n.in;
    [trainInd,valInd,testInd] = dividerand(nsamp.tra-CG,trainR,valR,testR);
else
    [trainInd,valInd,testInd] = dividerand(nsamp.tra,trainR,valR,testR);
    CG = 0;
end  

% Predefinition of the matrices
if true
    xload    = cell(n.PLS,n.red);
    xloads    = xload;
    beta      = xload;
    pctvar    = xload;
    load_used = xload;
    sumpct    = xload;
    diffloads = xload;
    stats     = xload;
    norm_comb = xload;
    net       = xload;
    norm_comb_mean      = xload;
    norm_comb_mean(1,:) = num2cell(ones(1,n.red));
    norm_loads          = cell(n.PLS,n.red,n.in);
    norm_loads_mean     = norm_loads;
end

m = 1;

for k = 1:n.PLS
    
    % Adjustment of the sampling intervall
    if Np.tot(k) == Np.NNt(1)
        m = 2;
    elseif Np.tot(k) == Np.NNt(2)
        m = 3;
    end
    
    % Scaling of the independent variables
    [U.tra.scaled,Var.U.mu,Var.U.std] = zscore(U.tra.sel(1:Np.tot(k),:));
    U.val.scaled = (U.val.sel-OnN1.val*Var.U.mu)./(OnN1.val*Var.U.std);

    % Scaling of the dependent variables
    [x.tra.scaled,Var.x.mu,Var.x.std]  = zscore(x.tra.red(1:Np.tot(k),:));
    scl.red.var = OnN1.val*Var.x.std;
    scl.red.add = OnN1.val*Var.x.mu;
    
    for j = 1:n.red     
        % Calculation of the Partial least square regression
        [xload{k,j},~,~,~,beta{k,j},pctvar{k,j},~,stats{k,j}] = ...
            plsregress(U.tra.scaled(1:Np.tot(k),:),x.tra.scaled(1:Np.tot(k),j),n.in);
             
        load_used{k,j} = stats{k,j}.W;
             
        % Rearrangement and normalization of the loads
        xloads{k,j} = load_used{k,j}./(ones(n.in,1)*sqrt(diag(stats{k,j}.W'*stats{k,j}.W)'));
        sumpct{k,j} = sum(pctvar{k,j}(2,:) > opt.comptresh);

        % Definition of the index of the maximum absolute and
        % adjustment of the loads
        [~,in.m] = max(abs(xloads{k,j}));
        in.m = in.m+(0:n.in:n.in*(n.in-1));
        in.c = sign(xloads{k,j}(in.m)) ~= sign(xloads{max(k-1,1),j}(in.m));
        xloads{k,j}(:,in.c) = -xloads{k,j}(:,in.c);
        
        if k > 1
            diffloads{k,j} = xloads{k,j}-xloads{k-1,j};
            norm_comb{k,j}   = norm(diffloads{k,j}(:,1:min(sumpct{k,j},sumpct{k-1,j})),normtyp);
            norm_comb_mean{k,j} = mean(cell2mat(norm_comb(max(k-n.norm,1)+1:k,j)));
            for l = 1:n.in
                norm_loads{k,j,l} = norm(diffloads{k,j}(:,l),normtyp);
                norm_loads_mean{k,j,l} = mean(cell2mat(norm_loads(max(k-n.norm,1)+1:k,j,l)));
            end
        end
        
        if mod(Np.tot(k),Np.d(m)) == 0
            % Fitting of 10 neural networks to the data
            for l = 1:10
                % Calculation of the new independent variables
                C.tra = U.tra.scaled(1:Np.tot(k),:);
                C.val = U.val.scaled;

                % Calculation of network parameter
                hiddenLayerSize = [5 5];
                neti = cascadeforwardnet(hiddenLayerSize,'trainlm');
                neti.trainParam.max_fail = 100;
                neti.trainParam.epochs = 5e3;
                neti.trainParam.min_grad = 1e-8;

                % Set up Division of Data for Training, Validation, Testing
                neti.divideFcn = 'divideind';
                neti.divideParam.trainInd = [1:CG trainInd+CG];
                neti.divideParam.valInd = valInd+CG;
                neti.divideParam.testInd = testInd+CG;

                % Train the Network
                [net{k,l,j},tr] = train(neti,C.tra',x.tra.scaled(1:Np.tot(k),j)','UseParallel','yes'); % );% ,'UseGPU','yes');%

                % Calculation of the results using PLS
                x.ann{k,j}(:,l) = net{k,l,j}(C.val')'.* ...
                                  scl.red.var(:,j)+scl.red.add(:,j);
            end
            % Calculation of the absolute error
            err.val{k,j} = (x.val.red(:,j)-mean(x.ann{k,j},2)).*scl.err(:,j);


            % Calculation of the maximum, mean, and median error
            err.comb{k,j} = [max(abs(err.val{k,j}))
                                 mean(abs(err.val{k,j}))
                                 median(abs(err.val{k,j}))];
            err.MSE{k,j} = sum(err.val{k,j}.^2)/nsamp.val;
            err.RMSE{k,j} = sqrt(err.MSE{k,j});
        end
        
    end
    
    % Print of the number of fitted points
    if mod(Np.tot(k),100) == 0
        fprintf('%i points are fitted.\n',Np.tot(k))
    end
end

if Ver.save
    if strcmp(Ver.model,'Pip')
        save(['Localsave\Paper\Pip' num2str(n.in - 2) '_'...
              datestr(date,'yyyymmdd') '_AddSamp_diff_w_5e3.mat'],...
              'U','x','n','lbmub','opt','err','Np','net')
    elseif strcmp(Ver.model,'Rea')
        save(['Localsave\Paper\Rea_'...
              datestr(date,'yyyymmdd') '_AddSamp_diff_w_5e3.mat'],...
              'U','x','n','lbmub','opt','err','Np','net')
    end
end

%% PLot of the results of this analysis

for j = 1:n.red
    figure(1)
    subplot(3,n.red,j)
    hold on
    for l = 1:sumpct{end,j}
        plot(Np.tot(2:n.PLS),cell2mat(norm_loads_mean(1:n.PLS,j,l)))
        hold on
        plot([0 opt.Np],[0.2 0.2],'k--')
        ylim([0 1.5])
        ylabel('Norm of the Difference for each component')
    end
    subplot(3,n.red,n.red+j)    
    plot(Np.tot(2:n.PLS),cell2mat(norm_comb_mean(2:n.PLS,j)))    
    hold on
    plot([0 opt.Np],[0.2 0.2],'k--')
%     ylim([0 1.5])
    ylabel('Norm of the Difference for significant loads')
    
    err.test = cell2mat(err.comb(j,:));
    subplot(3,n.red,2*n.red+j)     
    plot(Np.totNN,err.test(2,:))  
%     ylim([0 .1])
    ylabel('Absolute Error')
    xlabel('Sample Number')
    
    figure(2)
    err.test = cell2mat(err.comb(j,:));
    subplot(1,n.red,j)
    loglog(cell2mat(norm_comb_mean(ind.NN(2:end),j,1)),err.test(2,2:end),'k')
    ylabel('Error')
    xlabel('Norm')
%     ylim([0 0.4])
end

% Calculation of the overall time
time.overall = toc(time.in.overall);
