%% PCE METAMODELING: ESTIMATION OF STATISTICAL MOMENTS
%
% This example showcases an application of sparse
% polynomial chaos expansion (PCE) to the estimation of
% the statistical moments of a computational model with random inputs.
% The model is a simply supported beam model that computes
% deflections at mid-span of a beam subjected to a uniform random load.
% The convergence behavior of the PCE-based and pure Monte-Carlo-based 
% estimation of the statistical moments is compared.

%% 1 - INITIALIZE UQLAB
%
% Clear all variables from the workspace, set the random number generator
% for reproducible results, and initialize the UQLab framework:
clearvars
rng(100,'twister')
uqlab

%% 2 - COMPUTATIONAL MODEL
%
% The simply supported beam model is shown in the following figure:
uq_figure
[I,~] = imread('SimplySupportedBeam.png');
image(I)
axis equal
set(gca,'visible','off')

%% 
% The forward model computes the deflection of the beam $V$
% at mid-span location, which reads:
%
% $$ V = \frac{ 5 p L^4 }{32 E b h^3}$$
%
% This computation is carried out by the function
% |uq_SimplySupportedBeam(X)| supplied with UQLab.
% The input variables of this function are gathered into the $N \times M$
% matrix |X|, where $N$ and $M$ are the numbers of realizations
% and input variables, respectively.
% The input variables are given in the following order:
%
% # $b$: beam width $(m)$
% # $h$: beam height $(m)$
% # $L$: beam length $(m)$
% # $E$: Young's modulus $(MPa)$
% # $p$: uniform load $(kN/m)$

%%
% Create a MODEL object from the function:
ModelOpts.mFile = 'uq_SimplySupportedBeam';
ModelOpts.isVectorized = true;

myModel = uq_createModel(ModelOpts);

%% 3 - PROBABILISTIC INPUT MODEL
%
% The simply supported beam model has five inputs,
% modeled by independent lognormal random variables.
% The detailed model is given in the following table:

%%
% Define an INPUT object with the following marginals:
Input.Marginals(1).Name = 'b'; % beam width
Input.Marginals(1).Type = 'Lognormal';
Input.Marginals(1).Moments = [0.15 0.0075]; % (m)

Input.Marginals(2).Name = 'h'; % beam height
Input.Marginals(2).Type = 'Lognormal';
Input.Marginals(2).Moments = [0.3 0.015]; % (m)

Input.Marginals(3).Name = 'L'; % beam length
Input.Marginals(3).Type = 'Lognormal';
Input.Marginals(3).Moments = [5 0.05]; % (m)

Input.Marginals(4).Name = 'E'; % Young's modulus
Input.Marginals(4).Type = 'Lognormal';
Input.Marginals(4).Moments = [3e10 4.5e9]; % (Pa)

Input.Marginals(5).Name = 'p'; % uniform load
Input.Marginals(5).Type = 'Lognormal';
Input.Marginals(5).Moments = [1e4 2e3]; % (N/m)

%%
% Create the INPUT object:
myInput = uq_createInput(Input);

%% 4 - POLYNOMIAL CHAOS EXPANSION (PCE) METAMODEL
%
% Select PCE as the metamodeling tool:
MetaOpts.Type = 'Metamodel';
MetaOpts.MetaType = 'PCE';

%%
% Select the sparse-favouring least-square minimization LARS
% for the coefficient calculation strategy:
MetaOpts.Method = 'OLS';

%% 
% Specify a sparse truncation scheme (hyperbolic norm with $q = 1$):
MetaOpts.TruncOptions.qNorm = 1;

%%
% Specify the range of the degrees to be compared
% by the adaptive algorithm:
MetaOpts.Degree = 2;

MetaOpts.Display = 'verbose';
%%
% The degree with the lowest Leave-One-Out cross-validation (LOO)
% error estimator is chosen as the final metamodel.

%% 5 - MONTE-CARLO- VS PCE-BASED ESTIMATION OF MOMENTS
%
% The moments of a PCE can be analytically calculated
% from its coefficients without additional sampling.
% In this section, the mean and standard deviation estimates obtained
% from Monte-Carlo (MC) samples of increasing size are compared to the ones
% based on the PCE of the same samples.
%
% The simply supported beam model with lognormal input distributions allows
% for the explicit calculation of mean and standard deviation, which are
% provided for reference:
mean_exact = 0.008368320689566;
std_exact = 0.002538676671701;

%%
% Specify a set of experimental design sizes to test:
NED = 200;

%%
% Initialize the arrays in which the results will be stored:
mean_MC = zeros(size(NED));
std_MC = zeros(size(NED));
mean_PCE = zeros(size(NED));
std_PCE = zeros(size(NED));

%% 
% Loop over the experimental design sizes. In each step, generate an 
% experimental design, create a PCE metamodel,
% and compute the corresponding mean and standard deviation:
for ii = 1: length(NED)
    
    % Get a sample from the probabilistic input model with LHS sampling
    X_ED = uq_getSample(NED(ii),'LHS');
    
    % Evaluate the full model on the current experimental design
    Y_ED = uq_evalModel(myModel,X_ED);
    
    % Calculate the moments of the experimental design
    mean_MC(ii) = mean(Y_ED);
    std_MC(ii) = std(Y_ED, 0, 1);
    
    % Use the sample as the experimental design for the PCE
    MetaOpts.ExpDesign.X = X_ED;
    MetaOpts.ExpDesign.Y = Y_ED;
    
    % Create and add the PCE metamodel to UQLab
    myPCE = uq_createModel(MetaOpts);
    
    % Calculate the mean and standard deviation from the PCE coefficients
    mean_PCE(ii) = myPCE.PCE.Moments.Mean;
    std_PCE(ii) = sqrt(myPCE.PCE.Moments.Var);

end

if ~exist('TOL', 'var')
    TOL = 1e-2;
end

CC = myPCE.PCE.Coefficients;
Basis = full(myPCE.PCE.Basis.Indices);  %full sparse into full matrix

[P , M] = size(Basis);

% Build format of printing
MyFormat = '     [';
for i=1:M
    MyFormat = strcat(MyFormat, '%3i');
end
MyFormat= strcat(MyFormat, ']\t\t\t%12.7f\n');
fprintf('\n%%--- List of coefficients (sorted by amplitude) ---%%\n');
fprintf('     [ alpha_1   ...   alpha_M]\t\t  Coefficient\n');

% Print the coefficients according to the amplitude 
% and if they are greater than a threshold
[SortCoefs, ind] = sort(abs(CC),'descend');
TheStd = sqrt(sum(CC(2:end).^2));

for i=1:P
    cc = CC(ind(i));    
    if (abs(cc)> TOL*TheStd )
        tmp = Basis(ind(i),:);
        fprintf(MyFormat,tmp,cc);
    end
end

