%% Simple Ammonia Loop - Analysis of Norm convergence.
% This program can be used to sample points of the section of the simple
% ammonia loop and other models, if desired. It uses a maximum number of
% points and samples this number using incremental sampling. PLS regression 
% is subseqeuently applied to all dependent variables simultaenously and
% the data is saved.
%
% 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';     % 'Ros', '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;      % true or false (Saving of the data)

% Loading of the needed parameters
ParSamp

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

Np.a = 25;
Np.b = 5;
Np.c = 5;
Np.max = 5e3;
Np.tot = Np.a:Np.c:opt.Np;
n.PLS  = (opt.Np-Np.a)/Np.c+1;
n.norm = 5*Np.c;

% 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 'Ros'
    % Definition of the anonymous function
    SampAn  = @(initRos,U)  SampRose(initRos,U);
    
    % Definition of the varied vector
    U.ini.tot = [0 0];
    
    % Definition of the calculated values and indices used
    ncal = [];
    ind.in  = 1:2;
    ind.out = 1;
          
    % Definition of the range for plotting
    lbmub = zeros(3,numel(U.ini.tot));
    lbmub(1,:) = [-2 -1];
    lbmub(2,:) = U.ini.tot;
    lbmub(3,:) = [2   3];
                      
    % Definition of the initial values
    init = [];
    
    % Definition of the calculation function
    CalAn = @(U,x) CalRose(U,x);
    
    % Rearrangement of 
    Ver.Ind = 'MolFlow';
    
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(1,:) = U.ini.tot.*[1   1  0.5 0.5*ones(1,itg.j-1)] ...
                          - [4e5 10 zeros(1,itg.j)];
    lbmub(2,:) = U.ini.tot;
    lbmub(3,:) = U.ini.tot.*[1   1  1.5 1.5*ones(1,itg.j-1)] ...
                          + [4e5 10 zeros(1,itg.j)];
    
%     lbmub = zeros(3,numel(U.ini.tot));
%     lbmub(1,:) = U.ini.tot.*[1   1  0.8*ones(1,itg.j)] ...
%                           - [4e5 10 zeros(1,itg.j)];
%     lbmub(2,:) = U.ini.tot;
%     lbmub(3,:) = U.ini.tot.*[1   1  1.2*ones(1,itg.j)] ...
%                           + [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

% Predefinition of the matrix
U.tra.tot = [];
U.tra.un  = [];
tic
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
opt.Np = Np.max;
toc
% 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);

% 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) Evaluation of the partial least square performance

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

% Predefinition of the matrices
if true
    xload      = cell(n.PLS,n.red);
    yload      = xload;
    xloads     = xload;
    load_used  = xload;
    xscores    = xload;
    yscores    = xload;
    stats      = xload;
    beta       = xload;
    pctvar     = xload;
    sumpct     = xload;
    diff_loads = xload;
    statsWl    = xload;
    statsWh    = xload;
    diff_mat   = cell(n.in,n.red);
    
    % 1-norm
    norm_comb_1           = xload;
    norm_comb_mean_1      = xloads;
    norm_comb_mean_1(1,:) = num2cell(ones(1,n.red));
    
    % 2-norm
    norm_comb_2           = xload;
    norm_comb_mean_2      = xloads;
    norm_comb_mean_2(1,:) = num2cell(ones(1,n.red));
    
    % inf-norm
    norm_comb_i           = xload;
    norm_comb_mean_i      = xloads;
    norm_comb_mean_i(1,:) = num2cell(ones(1,n.red));
    
    % frob-norm
    norm_comb_f           = xload;
    norm_comb_mean_f      = xloads;
    norm_comb_mean_f(1,:) = num2cell(ones(1,n.red));
    
    norm_loads = cell(n.PLS,n.red,n.in);
    norm_loads_mean = norm_loads;
    optpls = statset('UseParallel',true);
end

tic
for k = 1:n.PLS
    
    for j = 1:n.red     
        
        % 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),j));
        
        % Calculation of the Partial least square regression
        [xload{k,j},yload{k,j},xscores{k,j},yscores{k,j},...
                 beta{k,j},pctvar{k,j},~,stats{k,j}] = plsregress(U.tra.scaled,x.tra.scaled,n.in,'options',optpls);
         
        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);
        if k > 1
            diff_loads{k,j} = (xloads{k,j}-xloads{k-1,j});
            
            % 1-norm
            norm_comb_1{k,j}   = norm(diff_loads{k,j}(:,1:min(sumpct{k,j},sumpct{k-1,j})),1);
            norm_comb_mean_1{k,j} = mean(cell2mat(norm_comb_1(max(k-n.norm,1)+1:k,j)));
            
            % 2-norm
            norm_comb_2{k,j}   = norm(diff_loads{k,j}(:,1:min(sumpct{k,j},sumpct{k-1,j})),2);
            norm_comb_mean_2{k,j} = mean(cell2mat(norm_comb_2(max(k-n.norm,1)+1:k,j)));
            
            % inf-norm
            norm_comb_i{k,j}   = norm(diff_loads{k,j}(:,1:min(sumpct{k,j},sumpct{k-1,j})),inf);
            norm_comb_mean_i{k,j} = mean(cell2mat(norm_comb_i(max(k-n.norm,1)+1:k,j)));
            
            % inf-norm
            norm_comb_f{k,j}   = norm(diff_loads{k,j}(:,1:min(sumpct{k,j},sumpct{k-1,j})),'fro');
            norm_comb_mean_f{k,j} = mean(cell2mat(norm_comb_f(max(k-n.norm,1)+1:k,j)));
            for l = 1:n.in
                norm_loads{k,j,l} = norm(diff_loads{k,j}(:,l),'fro');
                norm_loads_mean{k,j,l} = mean(cell2mat(norm_loads(max(k-n.norm,1)+1:k,j,l)));
                diff_mat{l,j}(:,k) = diff_loads{k,j}(:,l);
            end
        end
        test2{j}(:,k) = xloads{k,j}(:,1);
        test3{j}(:,k) = mean(test2{j}(:,max(k-n.norm,1)+1:k),2);
    end
    
    if mod(Np.tot(k),100) == 0
        fprintf('%i points are sampled.\n',Np.tot(k))
    end
    
end
toc

l = 1;
m = 1;
for k = 1:n.PLS
    for j = 1:n.red    
        statsWl{k,j} = xloads{k,j}(m,1);
        statsWh{k,j} = xloads{k,j}(m,sumpct{end}+1);
    end
end


itg.fig = 1;
for j = 1:n.red
    figure(itg.fig);
    subplot(2,n.red,j)
    hold on
    for l = 1:n.in%sumpct{end,j}
        loglog(Np.tot(2:n.PLS),cell2mat(norm_loads_mean(1:n.PLS,j,l)))
        hold on
        if j == 1
            ylabel('Norm of the Difference for each component')
        end
    end
%     loglog([0 opt.Np],opt.normtresh*[1 1],'k--')
    legend('1','2','3','4','5','6','7','8','9')
    
    subplot(2,n.red,n.red+j)    
    plot(Np.tot(2:n.PLS),cell2mat(norm_comb_mean_1(2:n.PLS,j)))   
    hold on
    plot(Np.tot(2:n.PLS),cell2mat(norm_comb_mean_2(2:n.PLS,j))) 
    plot(Np.tot(2:n.PLS),cell2mat(norm_comb_mean_i(2:n.PLS,j)))  
    plot([0 opt.Np],opt.normtresh*[1 1],'k--')
    if j == 1
        ylabel('Norm of the Difference for significant loads')
    end
end

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

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