%% Simple Ammonia Loop - Evaluation of the combined weights
% This program can be used to calculate the norm of the difference in the
% weights of the significant PLS regression components. It requires
% precalculated data which it loads for the recalculation of the norms. It
% allows to use different parameters in the calculation of the norms, and
% hence, sthe evaluation of the convergence of the norm. It applies PLS
% regression to all dependent variable combined.
%
% This program is complementary to the program "ProgNormCal_ind.m".
%
% 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 = 5;
% 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';     % 'Rea', 'Sep'
Ver.Ind    = 'MolFrac'; % 'MolFlow', 'MolFrac', or 'PartPress'
Ver.dep    = true;    	% true or false (Incorporation of dependency)
Ver.diff   = true;    	% true or false (Difference of absolute values)
Ver.plot   = true;

% Loading of the needed parameters
ParSamp

% Definition of the sampling grid options
opt.val = 'y';
opt.type = 'LHS';
opt.normtyp = 'fro';
opt.comptresh = 0.02;
opt.normtresh = 0.05;

% 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);
optpls = statset('UseParallel',true);

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

switch Ver.model
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.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,U.ini.out] = fsamp(SampAn,U.ini.tot,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.')

% Load of the data
switch Ver.model
    case 'Rea'
        load(['Localsave\Paper\Rea_20180902_AddSamp_dep_diff_all_LHSaug.mat'],'Utrain','xtrain','Np','n')
    case 'Sep'
        load(['Localsave\Paper\Sep_20180902_AddSamp_dep_diff_all_LHSaug.mat'],'Utrain','xtrain','Np','n')
end

% Reassignment of variables
U.tra.tot = Utrain;
x.tra.tot = xtrain;

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

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

%% a) Calculation of the Norm

% Calculation of several values
% Np.a = 25;
% Np.c = 5;
opt.Np = size(x.tra.red,1);
Np.tot = Np.a:Np.c:opt.Np;
n.PLS  = (opt.Np-Np.a)/Np.c+1;
if exist('xload','var')
    n.norm = 5;
else
    n.norm = n.norm/Np.c;
end
opt.comptresh = 0.02;
% When the separation section is used, then we have to use the following
% term.
% [1:8 10:end]

% Predefinition of the matrices
if true
    xload      = cell(n.PLS,1);
    yload      = xload;
    xloads     = xload;
    load_used  = xload;
    xscores    = xload;
    yscores    = xload;
    stats      = xload;
    beta       = xload;
    pctvar     = xload;
    sumpct     = xload;
    diff_loads = xload;
    diff_mat   = cell(n.in,1);
    
    % 1-norm
    norm_comb_1           = xload;
    norm_comb_mean_1      = xloads;
    
    % 2-norm
    norm_comb_2           = xload;
    norm_comb_mean_2      = xloads;
    
    % inf-norm
    norm_comb_i           = xload;
    norm_comb_mean_i      = xloads;
    
    % frob-norm
    norm_comb_f           = xload;
    norm_comb_mean_f      = xloads;
end
for k = 1:n.PLS
        % Scaling of the independent variables
        [U.tra.scaled,Var.U.mu,Var.U.std] = zscore(U.tra.sel(1:Np.tot(k),:));

        % Scaling of the dependent variables
        [x.tra.scaled,Var.x.mu,Var.x.std]  = zscore(x.tra.red(1:Np.tot(k),:));
        
        % Calculation of the Partial least square regression
        [xload{k},yload{k},xscores{k},yscores{k},...
                 beta{k},pctvar{k},~,stats{k}] = plsregress(U.tra.scaled(:,:),x.tra.scaled,n.in,'options',optpls);
         
        load_used{k} = stats{k}.W;
        % Rearrangement and normalization of the loads
        xloads{k} = load_used{k}./(ones(n.in,1)*sqrt(diag(stats{k}.W'*stats{k}.W)'));
        sumpct{k} = sum(pctvar{k}(2,:) > opt.comptresh);
        
        % Definition of the index of the maximum absolute and
        % adjustment of the loads
        [~,in.m] = max(abs(xloads{k}));
        in.m = in.m+(0:n.in:n.in*(n.in-1));
        in.c = sign(xloads{k}(in.m)) ~= sign(xloads{max(k-1,1)}(in.m));
        xloads{k}(:,in.c) = -xloads{k}(:,in.c);

        % Calculation of the norm and the difference
        if k > 1
            diff_loads{k} = (xloads{k}-xloads{k-1});
            
            % 1-norm
            norm_comb_1{k}   = norm(diff_loads{k}(:,1:min(sumpct{k},sumpct{k-1})),1);
            norm_comb_mean_1{k} = mean(cell2mat(norm_comb_1(max(k-n.norm,1)+1:k)));
            
            % 2-norm
            norm_comb_2{k}   = norm(diff_loads{k}(:,1:min(sumpct{k},sumpct{k-1})),2);
            norm_comb_mean_2{k} = mean(cell2mat(norm_comb_2(max(k-n.norm,1)+1:k)));
            
            % inf-norm
            norm_comb_i{k}   = norm(diff_loads{k}(:,1:min(sumpct{k},sumpct{k-1})),inf);
            norm_comb_mean_i{k} = mean(cell2mat(norm_comb_i(max(k-n.norm,1)+1:k)));
            
            % Frobenius norm
            norm_comb_f{k}   = norm(diff_loads{k}(:,1:min(sumpct{k},sumpct{k-1})),'fro');
            norm_comb_mean_f{k} = mean(cell2mat(norm_comb_f(max(k-n.norm,1)+1:k)));
        end
    
    if mod(Np.tot(k),100) == 0
        fprintf('%i points are sampled.\n',Np.tot(k))
    end
    
end

%% b) Plot of the results

if Ver.plot
    
    close all
    % Definition of integer values
    itg.fig = 1;
    itg.lab = 1;
    itg.fda = 1;
    itg.axe = 1;
    itg.ann = 1;

    % Definition of figure and values
    fig.dimw12 = 12.3;  % 12.5 for ESCAPE
    fig.dimh12 = 5.5;
    fig.dimw21 = 8.4;
    fig.dimh21 = 12;
    fig.dimw = 8.4;
    fig.dimh = 6;
    fig.Font = [10 8];
    fig.label = [];
    fig.ann   = [];
    fig.dat   = [];

    % Definition of a colour string
    rk = {'k','k','r','r','k','k','r','r','k','k','b','y','g'};%[0.93 0.69 0.13]}
    ls = {'-','--','-','--','-','--','-','--','-','--','-','--','-','--','-','--'};

    fig.in(itg.fig) = figure(itg.fig);
    set(fig.in(itg.fig),'units','centimeters','position', [5+itg.fig,5+itg.fig, fig.dimw, fig.dimh]);
    fig.axes(itg.axe) = gca;
    hold on
    fig.dat(itg.fda) = plot(Np.tot(2:n.PLS),cell2mat(norm_comb_mean_f(2:n.PLS)),'r'); itg.fda = itg.fda+1;
    fig.dat(itg.fda) = plot([0 Np.tot(n.PLS)],opt.normtresh*[1 1],'k-.'); itg.fda = itg.fda+1;
%     legend('$\Delta\mathbf{W}_{k,\Delta p}$','$\Delta\mathbf{W}_{k,\Delta T}$','$\Delta\mathbf{W}_{k,\xi}$','$\Delta\mathbf{W}_{k,Q_{HEx5}}$')
    set(fig.axes(itg.axe).Legend,'Location','Best','FontSize',fig.Font(1),'Interpreter','latex','Box','off')
    fig.label(itg.lab) = xlabel('Number of sampling points $n_p$'); itg.lab = itg.lab+1;
    fig.label(itg.lab) = ylabel('Norm $\left|\left|\Delta\mathbf{W}_s^k\right|\right|_{F}^{av}$'); itg.lab = itg.lab+1;
%     fig.axes(itg.axe).YScale = 'log';
%     fig.axes(itg.axe).XScale = 'log';
    if strcmp(Ver.model,'Rea')
%         fig.axes(itg.axe).YTick = 0:0.1:0.7;
%         fig.axes(itg.axe).YTickLabel = {'0.0';'0.1';'0.2';'0.3';'0.4';'0.5';'0.6';'0.7';'0.8';'0.9';'1.0'};
        fig.axes(itg.axe).YTick = 0:0.2:1.2;
        fig.axes(itg.axe).YTickLabel = {'0.0';'0.2';'0.4';'0.6';'0.8';'1.0';'1.2'};
        fig.ann(itg.ann) = annotation('textbox',[.20 .24 .05 .05],'String','$\gamma$'); itg.ann = itg.ann+1;
    elseif strcmp(Ver.model,'Sep')
%         fig.axes(itg.axe).YTick = 0:0.5:3.0;
%         fig.axes(itg.axe).YTickLabel = {'0.0';'0.5';'1.0';'1.5';'2.0';'2.5';'3.0'};
        fig.axes(itg.axe).YTickLabel = {'0.0';'0.2';'0.4';'0.6';'0.8';'1.0';'1.2';'1.4'};
        fig.ann(itg.ann) = annotation('textbox',[.20 .20 .05 .05],'String','$\gamma$'); itg.ann = itg.ann+1;
    end
    xlim([0 opt.Np])

    itg.axe = itg.axe+1;
    itg.fig = itg.fig+1;
 
    fig.in(itg.fig) = figure(itg.fig);
    set(fig.in(itg.fig),'units','centimeters','position', [5+itg.fig,5+itg.fig, fig.dimw, fig.dimh]);
    fig.axes(itg.axe) = gca;
    hold on
    fig.dat(itg.fda) = plot(Np.tot(2:n.PLS),cell2mat(sumpct(2:n.PLS)),'r'); itg.fda = itg.fda+1;
    legend('$\Delta\mathbf{W}_{k,\Delta p}$','$\Delta\mathbf{W}_{k,\Delta T}$','$\Delta\mathbf{W}_{k,\xi}$','$\Delta\mathbf{W}_{k,Q_{HEx5}}$')
    set(fig.axes(itg.axe).Legend,'Location','Best','FontSize',fig.Font(1),'Interpreter','latex','Box','off')
    fig.label(itg.lab) = xlabel('Number of sampling points $n_p$'); itg.lab = itg.lab+1;
    fig.label(itg.lab) = ylabel('Significant weights $n_s$'); itg.lab = itg.lab+1;
    ylim([2 7])
    xlim([0 opt.Np])
    itg.axe = itg.axe+1;
    itg.fig = itg.fig+1;   

    set(fig.ann, 'Interpreter', 'Latex');
    set(fig.ann, 'Fontsize', fig.Font(1));
    set(fig.ann, 'LineStyle', 'none');
    set(fig.label, 'Fontsize', fig.Font(1));
    set(fig.label, 'Interpreter', 'Latex');
    set(fig.axes,'FontSize', fig.Font(2));
    set(fig.axes,'Box','on');
    set(fig.axes,'TickLabelInterpreter', 'Latex');
    set(fig.axes,'LineWidth',.75);
    set(fig.axes,'XMinorTick', 'on');
    set(fig.axes,'YMinorTick', 'on');
    set(fig.dat,'LineWidth',.75)
end