%PCE METAMODELING: 
%This example showcases how to perform polynomial chaos expansion (PCE) metamodeling using existing data sets. The data sets come from a finite element model of a truss structure and are retrieved from different MAT-files. The files consist of an experimental design of size $200$ and a validation basis of size $10^4$.
%More information about the truss structure model can be found in the README.txt file located in the same folder as the truss data set.

%Contents
%1 - INITIALIZE UQLAB
%2 - RETRIEVE DATA SETS
%3 - INPUT MODEL
%4 - POLYNOMIAL CHAOS EXPANSION (PCE) METAMODEL
%5 - VALIDATION
%1 - INITIALIZE UQLAB
%Clear all variables from the workspace and initialize the UQLab framework:
clear all
close all
clearvars
uqlab
%2 - RETRIEVE DATA SETS
%The experimental design and the validation basis are stored in two separate files in the following location:
%FILELOCATION = fullfile(uq_rootPath, 'Examples', 'SimpleDataSets',...
    %'Truss_Matlab_FEM');
%Read the experimental design data set file and store the contents in matrices:
%load(fullfile(FILELOCATION,'Truss_Experimental_Design.mat'), 'X', 'Y');
%Read the validation basis data set file and store the contents in matrices:
%load(fullfile(FILELOCATION,'Truss_Validation_Basis.mat'), 'Xval', 'Yval');

opts = spreadsheetImportOptions("NumVariables", 11);

% Specify sheet and range
opts.Sheet = "DOE for metamodel";
opts.DataRange = "A2:K81";

% Specify column names and types
opts.VariableNames = ["SMYS", "UTS", "P", "dmin", "dc", "D", "t", "Wc", "Lc", "SL", "Sc"];
opts.VariableTypes = ["double", "double", "double", "double", "double", "double", "double", "double", "double", "double", "double"];

% Import the data
X = readtable("C:\Users\r01am21\OneDrive - University of Aberdeen\Documents\RESEARCH BOX\SFERM\UQLab\DOE STATIC PFEM EVD geo CSK100DP 141022rev SFERM2023.xlsx", opts, "UseExcel", false);

X = table2array(X)
clear opts
opts = spreadsheetImportOptions("NumVariables", 11);

% Specify sheet and range
opts.Sheet = "DOE for metamodel";
opts.DataRange = "A82:K98";

% Specify column names and types
opts.VariableNames = ["SMYS", "UTS", "P", "dmin", "dc", "D", "t", "Wc", "Lc", "SL", "Sc"];
opts.VariableTypes = ["double", "double", "double", "double", "double", "double", "double", "double", "double", "double", "double"];

% Import the data
Xval = readtable("C:\Users\r01am21\OneDrive - University of Aberdeen\Documents\RESEARCH BOX\SFERM\UQLab\DOE STATIC PFEM EVD geo CSK100DP 141022rev SFERM2023.xlsx", opts, "UseExcel", false);
Xval = table2array(Xval)
clear opts

opts = spreadsheetImportOptions("NumVariables", 1);

% Specify sheet and range
opts.Sheet = "DOE for metamodel";
opts.DataRange = "L2:L81";

% Specify column names and types
opts.VariableNames = "Pb";
opts.VariableTypes = "double";

% Import the data
Y = readtable("C:\Users\r01am21\OneDrive - University of Aberdeen\Documents\RESEARCH BOX\SFERM\UQLab\DOE STATIC PFEM EVD geo CSK100DP 141022rev SFERM2023.xlsx", opts, "UseExcel", false);
Y = table2array(Y)
clear opts

opts = spreadsheetImportOptions("NumVariables", 1);

% Specify sheet and range
opts.Sheet = "DOE for metamodel";
opts.DataRange = "L82:L98";

% Specify column names and types
opts.VariableNames = "Pb";
opts.VariableTypes = "double";

% Import the data
Yval = readtable("C:\Users\r01am21\OneDrive - University of Aberdeen\Documents\RESEARCH BOX\SFERM\UQLab\DOE STATIC PFEM EVD geo CSK100DP 141022rev SFERM2023.xlsx", opts, "UseExcel", false);
Yval = table2array(Yval)
clear opts
%3 - INPUT MODEL
%Because PCE requires a choice of polynomial basis, a probabilistic input model needs to be defined. Specify the marginals of the probabilistic input model:
% Specify the marginals as follows:The probabilistic input model consists of eleven (11) independent random variables.

InputOpts.Marginals(1).Name = 'SMYS';  % Tensile Yield Strength
InputOpts.Marginals(1).Type = 'Gaussian';
InputOpts.Marginals(1).Parameters = [624.48 32.97];  % (MPa)

InputOpts.Marginals(2).Name = 'UTS';  % Ultimate Tensile Strength
InputOpts.Marginals(2).Type = 'Gaussian';
InputOpts.Marginals(2).Parameters = [732 28.52];  % (MPa)

InputOpts.Marginals(3).Name = 'P';  % Internal Pressure
InputOpts.Marginals(3).Type = 'Gaussian';
InputOpts.Marginals(3).Parameters = [25, 0.559];  % (MPa)

InputOpts.Marginals(4).Name = 'dmin';  % dmin
InputOpts.Marginals(4).Type = 'Gaussian';
InputOpts.Marginals(4).Parameters = [3.77 2.060];  % (mm)

InputOpts.Marginals(5).Name = 'dc';  %  dmax or dcluster
InputOpts.Marginals(5).Type = 'Gaussian';
InputOpts.Marginals(5).Parameters = [4.790 1.544];  % (mm)

InputOpts.Marginals(6).Name = 'D';  %  Outer Diameter
InputOpts.Marginals(6).Type = 'Gaussian';
InputOpts.Marginals(6).Parameters = [458.6 0.691];  % (mm)

InputOpts.Marginals(7).Name = 't';  % Pipe wall thickness
InputOpts.Marginals(7).Type = 'Gaussian';
InputOpts.Marginals(7).Parameters = [7.98 0.071];  % (mm)

InputOpts.Marginals(8).Name = 'Wc';  % Parameter 1 Effective corrosion cluster width
InputOpts.Marginals(8).Type = 'Gaussian';
InputOpts.Marginals(8).Parameters = [129.42 6.471];  % (mm)

InputOpts.Marginals(9).Name = 'Lc';  % Parameter 2 Effective corrosion cluster length
InputOpts.Marginals(9).Type = 'Gaussian';
InputOpts.Marginals(9).Parameters = [210 10.5];  % (mm)

InputOpts.Marginals(10).Name = 'SL';  % Parameter 3 Interacting longitudinal distance
InputOpts.Marginals(10).Type = 'Gaussian';
InputOpts.Marginals(10).Parameters = [30 1.5];  % (mm)

InputOpts.Marginals(11).Name = 'Sc';  % Parameter 4 Interacting circumferential distance
InputOpts.Marginals(11).Type = 'Gaussian';
InputOpts.Marginals(11).Parameters = [21.08 1.054];  % (mm)


%Create an INPUT object based on the specified marginals:
myInput = uq_createInput(InputOpts);

%Specify an INPUT object: an INPUT object myInput can be specified with:
MetaOpts.Input = myInput;
% 4.4A Select PCE as the metamodeling tool: Sparse regression methods (LARS, OMP, SP, BCS)
MetaOpts.Type = 'Metamodel';
MetaOpts.MetaType = 'PCE';
% Specification of 
MetaOpts.Degree = 14;
MetaOpts.Method = 'LARS';
MetaOpts.LARS.KeepIterations = 1;
%Use experimental design loaded from the data files:
MetaOpts.ExpDesign.X = X;
MetaOpts.ExpDesign.Y = Y;
%Set the maximum polynomial degree to 20:
%MetaOpts.Degree = 1:20;
%Provide the validation data set to get the validation error:
MetaOpts.ValidationSet.X = Xval;
MetaOpts.ValidationSet.Y = Yval;
%Create the metamodel object and add it to UQLab:
myPCE_LARS = uq_createModel(MetaOpts)
%Specify a MODEL object: a MODEL object myModel can be specified with:
MetaOpts.FullModel = myPCE_LARS;
%Print a summary of the resulting PCE metamodel:
%uq_print(myPCE_sparseLARS)
%5 - VALIDATION
%Evaluate the PCE metamodel at the validation set points:
YPCE_LARS = uq_evalModel(myPCE_LARS,Xval);
%Print the validation and leave-one-out (LOO) cross-validation errors:
%fprintf('PCE metamodel validation error: %5.4e\n', myPCE_sparseLARS.Error.Val)
%fprintf('PCE metamodel LOO error:        %5.4e\n', myPCE_sparseLARS.Error.LOO)