%% Simple Ammonia Loop - Sampling using PLS regression individually
% This program can be used to sample points of the section of the simple
% ammonia loop and other models, if desired. It allows predetermined
% sampling and incremental sampling using the PLS regression based sampling
% algorithm of "J. Straus and S. Skogestad. A new termination criterion for 
% sampling for surrogate model generation using partial least squares 
% regression, Comp. Chem. Eng. 121, 2019, 75-85".
%
% PLS regression is applied to all dependent variables individually. The
% program "ProgSammp_all.m" is complementary and applieds PLS regression
% simultaenously.
%
% Surrogate models based on aritifical neural networks are fitted after the
% sampling.
%
% 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.TraLoad = false;    % Loading saved data for the training space or not
Ver.ValLoad = true;     % Loading saved data for the validation space or not
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)

% 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;

Np.a = 25;
Np.b = 5;
Np.c = 5;
Np.max = 5e3;
if ~Ver.add
    Np.a = Np.max;
end
Np.tot = Np.a:Np.b:Np.max;

n.PLS  = (Np.max-Np.a)/Np.b+1;
n.norm = 5*Np.c;
n.norm = 5*Np.c;      % has to be 2<= n.norm < 5 for without filtering

% 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);

% 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
    % 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
    % Definition of the range for plotting
    lbmub = zeros(3,numel(U.ini.tot));
    lbmub(1,:) = U.ini.tot.*[1   1  0.6 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.4 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)] ...
%                           - [3e5 5 zeros(1,itg.j)];
%     lbmub(2,:) = U.ini.tot;
%     lbmub(3,:) = U.ini.tot.*[1   1  1.2*ones(1,itg.j)] ...
%                           + [3e5 5 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.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.')
    
%% Part a) Calculation of the fitting points
% If multiple runs should be conducted, uncomment the following line the
% lines at the end of this section.
% for om = 1:20
    
% Predefinition of the matrices
if true
    U.tra.un        = [];
    U.tra.tot       = [];
    U.tra.scl       = [];
    U.tra.samp      = [];
    U.tra.Molfrac   = [];
    U.tra.PartPress = [];
    x.tra.tot       = [];
    x.tra.scl       = [];
    x.tra.red       = [];
    OnN1.tra        = [];
    exitflag.tra    = [];
    time.tra        = [];
    nsamp.tra       = 0;
    time.inctot     = 0;

    xloads          = cell(Np.max,n.red);
    pctvar          = xloads;
    sumpct          = xloads;
    diff_loads      = xloads;
    stats           = xloads;
    norm_comb       = xloads;
    norm_comb_mean  = xloads;
    xloadsum        = cell(1,n.red);
    normat = ones(n.norm,n.red);
end

time.in.train = tic;
time.PLS = 0;
if Ver.TraLoad
    switch Ver.model
        case 'Pip'
            load([dirsave Ver.model num2str(numel(U.ini.tot)-2) '_20170509_1e4LHS_dep_add.mat'],'Uval','xval')
        case 'Rea'
            load([dirsave 'Rea_20180227_1e4RND_dep_diff.mat'],'Uval','xval')
        case 'Sep'
            load([dirsave 'Rea_20180227_1e4RND_dep_diff.mat'],'Uval','xval')
    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);

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.tot = Ufull.out;
    U.inc.un  = Ufull.grid;
    
    % Sampling of the data
    time.in.inctot = tic;
    [x.inc.tot,cal.inc,b.inc,exitflag.inc,time.inc,U.inc.out] = fsamp(SampAn,U.inc.tot,x.ini.tot',ncal);
    time.inctot = toc(time.in.inctot) + time.inctot;
    
    % Calculation of the necessary variables
    [U.inc,x.inc,OnN1.inc,nsamp.inc] = CalAn(U.inc,x.inc);
    
    % Combination of the variables
    if true
        U.tra.un        = [U.tra.un; U.inc.un];
        U.tra.tot       = [U.tra.tot; U.inc.out];
        U.tra.samp      = [U.tra.samp; U.inc.samp];
        U.tra.Molfrac   = [U.tra.Molfrac; U.inc.Molfrac];
        U.tra.PartPress = [U.tra.PartPress; U.inc.PartPress];
        x.tra.tot       = [x.tra.tot; x.inc.tot];
        x.tra.red       = [x.tra.red; x.inc.red];
        exitflag.tra    = [exitflag.tra; exitflag.inc];
        time.tra        = [time.tra; time.inc];
        OnN1.tra        = [OnN1.tra; OnN1.inc];
        nsamp.tra       = nsamp.tra+nsamp.inc;
        
        switch Ver.model
        case {'ReacHysys','HexHYSYS'}
            U.tra.scl       = [U.tra.scl; U.inc.scl];
            x.tra.scl       = [x.tra.scl; x.inc.scl];
        end
    end

    % 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
    
    % Scaling of the independent and dependent variables
    U.tra.scaled = zscore(U.tra.sel);
    x.tra.scaled = zscore(x.tra.red);
    time.in.PLS = tic;
    if k == 1
        for l = 1:n.red     
            % Calculation of the Partial least square regression
            [~,~,~,~,...
                     ~,pctvar{k,l},~,stats{k,l}] = plsregress(U.tra.scaled,x.tra.scaled(:,l),n.in);
            
            
            xload = stats{k,l}.W;
            % Rearrangement and normalization of the loads
            xloads{k,l} = xload./(ones(n.in,1)*sqrt(diag(stats{k,l}.W'*stats{k,l}.W)'));
            sumpct{k,l} = sum(pctvar{k,l}(2,:) > opt.comptresh);
        end
    else
        for j = Np.b*(k-2)+Np.c:Np.c:Np.b*(k-1)+Np.c-1
            for l = 1:n.red     
                % Calculation of the Partial least square regression
                [~,~,~,~,~,pctvar{j,l},~,stats{j,l}] = ...
                    plsregress(U.tra.scaled(1:Np.a+j,:),x.tra.scaled(1:Np.a+j,l),...
                               n.in,'options',optpls);

                xload = stats{j,l}.W;
                
                % Rearrangement and normalization of the loads
                xloads{j,l} = xload./(ones(n.in,1)*sqrt(diag(stats{j,l}.W'*stats{j,l}.W)'));
                sumpct{j,l} = sum(pctvar{j,l}(2,:) > opt.comptresh);
                
                % Definition of the index of the maximum absolute and
                % adjustment of the loads
                [~,in.m] = max(abs(xloads{j,l}));
                in.m = in.m+(0:n.in:n.in*(n.in-1));
                in.c = sign(xloads{j,l}(in.m)) ~= sign(xloads{max(j-Np.c,1)}(in.m));
                xloads{j,l}(:,in.c) = -xloads{j,l}(:,in.c);
                
                % Calculation of the norm and the difference
                diff_loads{j,l} = xloads{j,l}-xloads{max(j-Np.c,1),l};
                norm_comb{j,l} = norm(diff_loads{j,l}(:,1:min(sumpct{j,l},sumpct{max(j-Np.c,1),l})),opt.normtyp);
                norm_comb_mean{j,l} = mean(cell2mat(norm_comb(max(j+1-n.norm,0)+1:j,l)));
            end
        end
        if all(cell2mat(norm_comb_mean(j,:)) < opt.normtresh)
            break
        end
    end
    time.PLS = time.PLS+toc(time.in.PLS);
    if rem(Np.a+Np.b*(k-1),100) == 0
        fprintf('%4i points were sampled, but do not suffice.\n',Np.a+Np.b*(k-1))
    end
end
end
% Removal of the cells not necessary
itg.t      = Np.b*(k-1)+Np.c;
xloads     = xloads(1:itg.t,:);
pctvar     = pctvar(1:itg.t,:);
sumpct     = sumpct(1:itg.t,:);
stats      = stats(1:itg.t,:);
diff_loads = diff_loads(1:itg.t,:);
norm_comb  = norm_comb(1:itg.t,:);
norm_comb_mean  = norm_comb_mean(1:itg.t,:);
    
time.train = toc(time.in.train);
% Printing of the results
if Np.a+itg.t-1 > Np.max
    fprintf('The training space required more than %4i points and the loop broke.\n',Np.max)
else
    fprintf('The training space has been sampled with %4i points.\n',Np.a+itg.t-Np.c)
end

% If multiple runs should be conducted, uncomment the following lines
% Np.tra(om) = Np.a+itg.t-Np.c;
% end

%% Part b) Calculation of the validation points
% Definition of the grid options

if Ver.ValLoad
    % 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')
            dum = Uval;
            clear Uval            
            Uval.tot = dum;
            Uval.out = Uval.tot;
        case 'Rea'
            load([dirsave 'Rea_20180227_1e4RND_dep.mat'],'Uval','xval')
            dum = Uval;
            clear Uval            
            Uval.tot = dum;
            Uval.out = Uval.tot;
        case 'Sep'
            load([dirsave 'Sep_20180227_1e4RND_dep.mat'],'Uval','xval')
            dum = Uval;
            clear Uval            
            Uval.tot = dum;
            Uval.out = Uval.tot;
    end

else
    optval = opt;
    optval.type = 'RND';
    optval.val = 'y';
    optval.Np = 1e4;

    % Call of the grid definition function
    Uval = GridDef(lbmub,ind.in,optval);
    dum = Uval.out;    
    clear Uval
    Uval.tot = dum;

    % Sampling of the data
    time.in.Val = tic;
    [xval,cal.val,b.val,exitflag.val,time.val,Uval.out] = fsamp(SampAn,Uval.tot,x.ini.tot',ncal);
    time.Val = toc(time.in.Val);
end

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

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

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

%% Part c) Scaling and Calcuation of values
% This part of the code scales and calculates the values for the PLS
% regression.

% 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 = 70/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  

% Scaling of the independent variables
[U.tra.scaled,Var.U.mu,Var.U.std] = zscore(U.tra.sel);
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);
scl.red.var = OnN1.val*Var.x.std;
scl.red.add = OnN1.val*Var.x.mu;

%% Part d) PLS on the values of the reaction section and fitting of an ANN
% This part of this script calculates the extent of reaction for the
% different inlet conditions and then fits with PLS.

% n_comp = 9;
% n_fit = 3;

itg.fig = l;%*2+k*2*length(var_comp); % +it*100;
    
% Define the value, PLS should be fitted to (as row vector)
%   1 = Outlet p [Pa]
%   2 = Outlet T [C]
%   3 = Extent of reaction [kmol/s]
% var_fit = 1:n.red;
var_fit = [13];
% close all
% clear err, label
clear xload

label.row = {'Max' 'Mean' 'Median'};
switch Ver.model
    case 'Ros'
        label.RM = {'Y'};
        label.col2 = {'Y'};
        label.Leg = {'y','y_{lin}'};
        scl.err.val = 1;
        scl.err.tra = 1;
        hiddenLayerSize = [5];
    case 'Pip'
        label.RM = {'Pressure'};
        label.col2 = {'p'};
        label.Leg = {'p_{out}','p_{outin}'};
        scl.err.val = 10^-5;
        scl.err.tra = 10^-5;
%         scl.err.val = 1./U.val.samp(:,1)*100;
%         scl.err.val = 1./x.val.red*100;
%         scl.err.tra = 1./x.tra.red*100;
%         scl.err.val = 1./(max(x.val.red)-min(x.val.red))*100;
        hiddenLayerSize = [5 5];
    case 'Rea'
        label.RM = {'Pressure'; 'Temperature'; 'Extent of reaction'; 'Duty'};
        label.col2 = {'p'; 'T'; 'EoR'; 'P_HEx'};
        label.Leg = {'p_{out}','p_{outin}','T_{out}','T_{outin}','EoR','EoR_{lin}','W','W_{lin}'};
        scl.err.val = [10^-5 1 1 10^-3];
        scl.err.tra = [10^-5 1 1 10^-3];
%         scl.err.val = 1./x.val.red*100;
%         scl.err.val = 1./(max(x.val.red)-min(x.val.red))*100;
        hiddenLayerSize = [5 5];
    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'};
        label.col2 = {'p_Gas'; 'T_Gas'; 'p_Liq'; 'T_Liq'; 'T_H2O'; ...
                      'H2'; 'N2'; 'NH3'; 'Ar'; 'CH4'; 'P_HEx8'; 'P_HEx9'; 'P_HEx10'};
        scl.err.val = [ones(nsamp.val,1)*[10^-5 1 10^-5 1 1] U.val.tot(:,3:7) ones(nsamp.val,3)*10^-3];
        scl.err.tra = [ones(nsamp.tra,1)*[10^-5 1 10^-5 1 1] U.tra.tot(:,3:7) ones(nsamp.tra,3)*10^-3];
%         scl.err.val = 1./x.val.red*100;
%         scl.err.val = 1./(max(x.val.red)-min(x.val.red))*100;
        hiddenLayerSize = [5 5];
end

% Predefinition of vectors
expvar = cell(1,n.red);

for l = var_fit
% var_comp = [sumpct{itg.t,l} n.in+1];
var_comp = n.in+1; [1:n.in-1 n.in+1];

for j = var_comp
    if j < n.in+1
        label.Title{j} = [num2str(j),' Components'];
        label.col{j} = ['Comp',num2str(j)];
    else
        label.Title{j} = 'Full Components';
        label.col{j} = 'FullComp';
    end
end
for k = var_comp
    
    % Assignment of the varied variables
    if k < n.in
        n_comp = k;
    else
        n_comp = n.in;
    end
    n_fit = l;

    % Calculation of the Partial least square regression
    [xload{itg.t,l},yload,xscores,yscores,...
                par.red,pctvar{itg.t,l},~,stats{itg.t,l}] = plsregress(U.tra.scaled,x.tra.scaled(:,n_fit),n_comp,'options',optpls);

    % Plotting of the cumulated sums
    sumpct{itg.t,l} = sum(pctvar{itg.t,l}(2,:) > opt.comptresh);
    if k == max(var_comp)
        expvar{l} = cumsum(pctvar{end,l},2);
    end
    load_used = stats{itg.t,l}.W;
    xloads{itg.t,l} = load_used./(ones(n.in,1)*max(abs(load_used)));

    % Calculation of the new independent variables
    switch k
        case num2cell(1:n.in)
            mult{k,l} = xloads{itg.t,l};pinv(xloads{itg.t,l}'); 
        case n.in+1
            mult{k,l} = 1; 
    end
    C.tra = U.tra.scaled*mult{k,l};
    C.val = U.val.scaled*mult{k,l};

    % Calculation of network parameter
    net = cascadeforwardnet(hiddenLayerSize,'trainlm');
    net.trainParam.max_fail = 100;
    net.trainParam.epochs = 10e3;
    net.trainParam.min_grad = 1e-9;


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

    % Train the Network
    [net,tr] = train(net,C.tra',x.tra.scaled(:,n_fit)','UseParallel','yes'); % );% ,'UseGPU','yes');%

    % Calculation of the results using PLS
    x.val.ann(:,n_fit) = net(C.val')'.* ...
                         scl.red.var(:,n_fit)+scl.red.add(:,n_fit);
    x.tra.ann(:,n_fit) = net(C.tra')'.* ...
                         scl.red.var(1:nsamp.tra,n_fit)+scl.red.add(1:nsamp.tra,n_fit);
    x.val.lin(:,n_fit) = ([OnN1.val U.val.scaled]*par.red).* ...
                         scl.red.var(:,n_fit)+scl.red.add(:,n_fit);


    rnge.PLS = 1:nsamp.val;

    % Calculation of the absolute error
    err.lin{l,k} = (x.val.red(:,n_fit)-x.val.lin(:,n_fit)).*scl.err.val(:,n_fit);
    err.val{l,k} = (x.val.red(:,n_fit)-x.val.ann(:,n_fit)).*scl.err.val(:,n_fit);
    err.tra{l,k} = (x.tra.red(:,n_fit)-x.tra.ann(:,n_fit)).*scl.err.tra(:,n_fit);
    err.norm{l,k} = abs((x.val.red(:,n_fit)-x.val.ann(:,n_fit)))/...
                    (max(abs(x.val.red(:,n_fit)))-min(abs(x.val.red(:,n_fit))));

    % Calculation of the maximum, mean, and median error
    err.comb.val{l,k} = [max(abs(err.val{l,k}));
                         mean(abs(err.val{l,k}));
                         median(abs(err.val{l,k}))];
    err.comb.lin{l,k} = [max(abs(err.lin{l,k}));
                         mean(abs(err.lin{l,k}));
                         median(abs(err.lin{l,k}))];
    MSE.val{l,k} = sum(err.val{l,k}.^2)/nsamp.val;
    MSE.tra{l,k} = sum(err.tra{l,k}.^2)/nsamp.tra;
    RMSE.val{l,k} = sqrt(MSE.val{l,k});
    RMSE.tra{l,k} = sqrt(MSE.tra{l,k});
    
    % Plotting of the variables and comparison with the PLS fitted values
    figure(itg.fig); itg.fig = itg.fig+1;
    plot(x.val.red(:,n_fit),x.val.ann(:,n_fit),'ro')
    hold on
    plot(x.tra.red(:,n_fit),x.tra.red(:,n_fit),'ko')
    plot([min(x.val.red(:,n_fit)) max(x.val.red(:,n_fit))],...
         [min(x.val.red(:,n_fit)) max(x.val.red(:,n_fit))],'k')
    ylabel([label.RM{n_fit},' Surrogate']);
    xlabel([label.RM{n_fit},' Simulated']);
    
    % Definition of the net for saving
    combnet{k,l} = net;

end
end

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