%% INVERSION: SURROGATE MODEL ACCELERATED CALIBRATION
%% 1 - INITIALIZE UQLAB
clearvars
rng(100,'twister')
uqlab

%% 2 - FORWARD MODEL
FILELOCATION = fullfile(uq_rootPath, 'Examples', 'SimpleDataSets',...
    'WDD_Experiment');

load(fullfile(FILELOCATION,'wdd_experiment_doe40.mat'), 'X', 'Y');
% load('wdd_experiment_doe40.mat', 'X', 'Y');


load(fullfile(FILELOCATION,'wdd_experiment_val30'), 'Xval', 'Yval');
% load('wdd_experiment_val30', 'Xval', 'Yval');

%% 3 - INPUT MODEL
% Scaled prototype experimental dam
InputOpts.Marginals(1).Name = 'E1'; %Elastic modulus of the dam
InputOpts.Marginals(1).Type = 'Gaussian';
InputOpts.Marginals(1).Moments = [1.5e9 0.75e8]; % (Pa)

InputOpts.Marginals(2).Name = 'D1'; %The density of the dam
InputOpts.Marginals(2).Type = 'Gaussian';
InputOpts.Marginals(2).Moments = [2200 110]; % (kg/m3)

InputOpts.Marginals(3).Name = 'E2'; %Elastic modulus of the foundation
InputOpts.Marginals(3).Type = 'Gaussian';
InputOpts.Marginals(3).Moments = [3e10 1.5e9]; % (Pa)

InputOpts.Marginals(4).Name = 'D2'; %The density of the foundation
InputOpts.Marginals(4).Type = 'Gaussian';
InputOpts.Marginals(4).Moments = [2400 120]; % (kg/m3)

myInput = uq_createInput(InputOpts);
%% 4 - SURROGATE MODEL
MetaOpts.Type = 'Metamodel';
MetaOpts.MetaType = 'PCE';

%% 
MetaOpts.ExpDesign.X = X;
MetaOpts.ExpDesign.Y = Y;
%%
MetaOpts.Degree = 1:10;
%%
MetaOpts.ValidationSet.X = Xval;
MetaOpts.ValidationSet.Y = Yval;

%%
myPCE = uq_createModel(MetaOpts);
uq_print(myPCE)
uq_display(myPCE)

%%
%uq_evalModel(myPCE,[1.5E9 2200 3e10 2400]);
%% 5 - MEASUREMENT DATA
myData.y = [137,157,304,380,488,620];  % (hz)
myData.Name = 'Realfrequency';
%%
%Known residual variance
% DiscrepancyOpts.Type = 'Gaussian';
% DiscrepancyOpts.Parameters = 0.5 ;
% DiscrepancyOpts.Parameters = [0.5, 0.5, 0.5, 0.5, 0.5, 0.5] ;
% Bayes.Discrepancy = DiscrepancyOpts ;
%% Unknown residual variance_type1
% 0.01 is a first guest, values will have to be adapted 
DiscrepancyPriorOpts.Name = 'Prior of discrepancy parameter';
DiscrepancyPriorOpts.Marginals(1).Name = 'Sigma2, com1';
DiscrepancyPriorOpts.Marginals(1).Type = 'Uniform';
DiscrepancyPriorOpts.Marginals(1).Parameters = [0, 0.01 * 137^2];
DiscrepancyPriorOpts.Marginals(2).Name = 'Sigma2, com2';
DiscrepancyPriorOpts.Marginals(2).Type = 'Uniform';
DiscrepancyPriorOpts.Marginals(2).Parameters = [0, 0.01 * 157^2];
DiscrepancyPriorOpts.Marginals(3).Name = 'Sigma2, com3';
DiscrepancyPriorOpts.Marginals(3).Type = 'Uniform';
DiscrepancyPriorOpts.Marginals(3).Parameters = [0,0.01 *  304^2];
DiscrepancyPriorOpts.Marginals(4).Name = 'Sigma2, com4';
DiscrepancyPriorOpts.Marginals(4).Type = 'Uniform';
DiscrepancyPriorOpts.Marginals(4).Parameters = [0,0.01 *  380^2];
DiscrepancyPriorOpts.Marginals(5).Name = 'Sigma2, com5';
DiscrepancyPriorOpts.Marginals(5).Type = 'Uniform';
DiscrepancyPriorOpts.Marginals(5).Parameters = [0, 0.01 * 488^2];
DiscrepancyPriorOpts.Marginals(6).Name = 'Sigma2, com6';
DiscrepancyPriorOpts.Marginals(6).Type = 'Uniform';
DiscrepancyPriorOpts.Marginals(6).Parameters = [0, 0.01 * 620^2];
myDiscrepancyPrior = uq_createInput(DiscrepancyPriorOpts);
DiscrepancyOpts.Type = 'Gaussian';
DiscrepancyOpts.Prior = myDiscrepancyPrior;
BayesOpts.Discrepancy = DiscrepancyOpts

%% Unknown residual variance_type2
% assuming that you have a relative error of 10 
% could be represented by multiplying the upper boundaries by 0.01.
% 
% DiscrepancyPriorOpts.Marginals.Name = 'Sigma2';
% DiscrepancyPriorOpts.Marginals.Type = 'Uniform';
% DiscrepancyPriorOpts.Marginals.Parameters = [0, 620^2 * 0.01];
% myDiscrepancyPrior = uq_createInput(DiscrepancyPriorOpts);
% DiscrepancyOpts.Type = 'Gaussian';
% DiscrepancyOpts.Prior = myDiscrepancyPrior;

%% Known residual variance
% DiscrepancyOpts.Type = 'Gaussian';
% DiscrepancyOpts.Parameters = [0.5 0.5 0.5 0.5] ;
% Bayes.Discrepancy = DiscrepancyOpts ;

%%
Solver.MCMC.Steps = 500;
Solver.MCMC.NChains = 20;

%% 6 - BAYESIAN ANALYSIS
BayesOpts.Type = 'Inversion';
BayesOpts.Data = myData;
BayesOpts.Prior = myInput;
%%
BayesOpts.ForwardModel.Model = myPCE;

%%
% Solver.Type = 'MCMC';
% Solver.MCMC.Sampler = 'AIES'; % AM, HMC, AIES,MH
%%
%Solver.MCMC.Visualize.Parameters = Parameters;
%Solver.MCMC.Visualize.Interval = Interval;
%%
%BayesOpts.ForwardModel.Model = myForwardModel;
%%
% myBayesianAnalysis_fullModel = uq_createAnalysis(BayesOpts);

%%
% uq_print(myBayesianAnalysis_fullModel)

%%
% Run the Bayesian inversion analysis:
myBayesianAnalysis_surrogateModel = uq_createAnalysis(BayesOpts);
%%
uq_postProcessInversion(myBayesianAnalysis_surrogateModel,'priorPredictive',2000);
uq_display(myBayesianAnalysis_surrogateModel, 'trace', [1 2 3 4])
%%
% Print out a report of the results:
uq_print(myBayesianAnalysis_surrogateModel)
uq_display(myBayesianAnalysis_surrogateModel)

%%
% Comparing |myBayesianAnalysis_fullModel| and
% |myBayesianAnalysis_surrogateModel| it can be seen that the results are
% practically identical. The small differences come from the randomness of
% the MCMC algorithms.

%%
% The number of original forward model calls in MCMC with the full model
% was $N = 30{,}000$ compared to the $N = 50$ model evaluations necessary
% to compute the PCE surrogate. In cases where the original forward model
% is computationally expensive, accelerating MCMC with surrogate models
% result in significant reduction of the total computational costs.
%
% In this example, a PCE surrogate was employed, but generally any 
% surrogate model available in UQLab (e.g., Kriging, LRA, SVR) can be used.
