%% INVERSION: BAYESIAN LINEAR REGRESSION
%

%% 1 - INITIALIZE UQLAB
%
clearvars
rng(100,'twister')
uqlab

%% 2 - FORWARD MODEL
% Define the forward model as a UQLab MODEL object:
% 1st forward model
ModelOpts1.Name = 'Armour failure';
ModelOpts1.mFile = 'uq_ArmourFailure';

forwardModels(1).Model = uq_createModel(ModelOpts1);
forwardModels(1).PMap = [1 2 3 4 5];

% 2nd forward model
ModelOpts2.Name = 'Toe failure';
ModelOpts2.mFile = 'uq_ToeFailure';

forwardModels(2).Model = uq_createModel(ModelOpts2);
forwardModels(2).PMap = [1 2 5 7 8];

% 3rd forward model
ModelOpts3.Name = 'Overtopping failure';
ModelOpts3.mFile = 'uq_OvertoppingFailure';

forwardModels(3).Model = uq_createModel(ModelOpts3);
forwardModels(3).PMap = [4 5 6 9];

%% 3 - PRIOR DISTRIBUTION OF THE MODEL PARAMETERS
% Specify these distributions as a UQLab INPUT object:
  
  PriorOpts.Marginals(1).Name = 'Delta';
  PriorOpts.Marginals(1).Type = 'Gaussian';
  PriorOpts.Marginals(1).Parameters = [1.439 0.055];
  
  PriorOpts.Marginals(2).Name = 'Dn';
  PriorOpts.Marginals(2).Type = 'Gaussian';
  PriorOpts.Marginals(2).Parameters = [1.687 0.084];
  
  PriorOpts.Marginals(3).Name = 'KD';
  PriorOpts.Marginals(3).Type = 'Gaussian';
  PriorOpts.Marginals(3).Parameters = [-0.7985 0.352];
  
  PriorOpts.Marginals(4).Name = 'cot';
  PriorOpts.Marginals(4).Type = 'Gaussian';
  PriorOpts.Marginals(4).Parameters = [2.5 0.125];
  
  PriorOpts.Marginals(5).Name = 'Hs';
  PriorOpts.Marginals(5).Type = 'Gumbel';
  PriorOpts.Marginals(5).Parameters = [-0.7985 0.352];
  
  PriorOpts.Marginals(6).Name = 'Tp';
  PriorOpts.Marginals(6).Type = 'Gumbel';
  PriorOpts.Marginals(6).Parameters = [-9.0662 1.2568];
  
  PriorOpts.Marginals(7).Name = 'Nod';
  PriorOpts.Marginals(7).Type = 'Gaussian';
  PriorOpts.Marginals(7).Parameters = [4 1.2];
  
  PriorOpts.Marginals(8).Name = 'hd';
  PriorOpts.Marginals(8).Type = 'Gaussian';
  PriorOpts.Marginals(8).Parameters = [0.68 0.05];
  
  PriorOpts.Marginals(9).Name = 'WL';
  PriorOpts.Marginals(9).Type = 'Gaussian';
  PriorOpts.Marginals(9).Parameters = [1.6127 0.6474];

myPriorDist = uq_createInput(PriorOpts);

%% 4 - MEASUREMENT DATA
%
% Provide the set of measurements that are stored in the 'Data' vector to
% the measurement vector |y|

AllData = csvread('Book1.csv');
y1 = AllData(:,1);
y2 = AllData(:,3);
y3 = AllData(:,2);

% Data group 1
myData(1).y = y1;
myData(1).Name = 'Armour failure';
myData(1).MOMap = [ 1;... % Model ID
                    1];   % Output ID
% Data group 2
myData(2).y = y2;
myData(2).Name = 'Toe Failure';
myData(2).MOMap = [ 2;... % Model ID
                    1];   % Output ID
                
% Data group 3
myData(3).y = y3;
myData(3).Name = 'Overtopping Failure';
myData(3).MOMap = [ 3;... % Model ID
                    1];   % Output ID
                
%% 5 - DISCREPANCY MODEL
%
% To infer the discrepancy variance, uniform priors are put on the
% discrepancy parameters:
%
% * $\sigma^2_{\mathrm{1}} \sim \mathcal{U}(0, 0.01)$
% * $\sigma^2_{\mathrm{2}} \sim \mathcal{U}(0, 0.01)$
%
% Create two UQLab INPUT objects based on these distributions:
SigmaOpts.Marginals(1).Name = 'Sigma2a';
SigmaOpts.Marginals(1).Type = 'Uniform';
SigmaOpts.Marginals(1).Parameters = [0 0.01];

SigmaDist1 = uq_createInput(SigmaOpts);

SigmaOpts.Marginals(1).Name = 'Sigma2t';
SigmaOpts.Marginals(1).Type = 'Uniform';
SigmaOpts.Marginals(1).Parameters = [0 0.01];

SigmaDist2 = uq_createInput(SigmaOpts);

SigmaOpts.Marginals(1).Name = 'Sigma2o';
SigmaOpts.Marginals(1).Type = 'Uniform';
SigmaOpts.Marginals(1).Parameters = [0 0.01];

SigmaDist3 = uq_createInput(SigmaOpts);

%%
% Assign these distributions to the discrepancy model options:
DiscrepancyOpts(1).Type = 'Gaussian';
DiscrepancyOpts(1).Prior = SigmaDist1;
DiscrepancyOpts(2).Type = 'Gaussian';
DiscrepancyOpts(2).Prior = SigmaDist2;
DiscrepancyOpts(3).Type = 'Gaussian';
DiscrepancyOpts(3).Prior = SigmaDist3;

%% 6 - BAYESIAN ANALYSIS WITH KNOWN MEASUREMENT ERROR
%
%% 6.1 Solver options
% In this example, the Bayesian calibration is solved with 
% an adaptive Metropolis algorithm using $10^2$ chains of length $10^3$
% iterations each
Solver.Type = 'MCMC';
Solver.MCMC.Sampler = 'AIES';
Solver.MCMC.Steps = 200;
Solver.MCMC.NChains = 100;
%% 6.2 Bayesian analysis
%
% The options of the Bayesian analysis are specified with the following
% structure:
BayesOpts.Type = 'Inversion';
BayesOpts.Name = 'Bayesian model1';
BayesOpts.Prior = myPriorDist;
BayesOpts.ForwardModel = forwardModels;
BayesOpts.Data = myData;
BayesOpts.Discrepancy = DiscrepancyOpts;
BayesOpts.Solver = Solver;

%%
% Run the Bayesian inversion analysis:
myBayesianAnalysisKnownDisc = uq_createAnalysis(BayesOpts);

%%
% Print out a report of the posterior sample :
uq_print(myBayesianAnalysisKnownDisc)
