% Example in Doc sec 2.2.1 to 2.2.5 modified by Sec. 2.5 
fprintf(' Example in Doc, sec 2.2.1 to 2.2.5; modified by Sec. 2.5\n ')
clear all 
clearvars
uqlab
rng(100)
PriorOpts.Name = 'Model parameters prior';
PriorOpts.Marginals(1).Name = 'b';
PriorOpts.Marginals(1).Type = 'Constant';
PriorOpts.Marginals(1).Parameters = 0.15; % (m)
PriorOpts.Marginals(2).Name = 'h';
PriorOpts.Marginals(2).Type = 'Constant';
PriorOpts.Marginals(2).Parameters = 0.3; % (m)
PriorOpts.Marginals(3).Name = 'L';
PriorOpts.Marginals(3).Type = 'Constant';
PriorOpts.Marginals(3).Parameters = 5; % (m)
PriorOpts.Marginals(4).Name = 'E';
PriorOpts.Marginals(4).Type = 'Lognormal';
PriorOpts.Marginals(4).Moments = [30000 4500] ; % (MPa)
PriorOpts.Marginals(5).Name = 'p';
PriorOpts.Marginals(5).Type = 'Constant';
PriorOpts.Marginals(5).Parameters = 0.012; % (kN/m)

myPriorDist = uq_createInput(PriorOpts);

ModelOpts1.Name = 'Beam bending deflection';
ModelOpts1.mFile = 'uq_SimplySupportedBeam';
ForwardModels(1).Model = uq_createModel(ModelOpts1);
% ForwardModels(1).PMap = [1 2 3 4 5];

ModelOpts2.Name = 'Beam bending deflection 2';
ModelOpts2.mFile = 'uq_SimplySupportedBeam';
ForwardModels(2).Model = uq_createModel(ModelOpts2);
% ForwardModels(2).PMap = [1 2 3 4 6];


V_mid = [12.84, 13.12, 12.13, 12.19, 12.67]/1000; % (m)

% Data group 1
myData.y = V_mid;
myData.Name = 'Beam mid-span deflection';
myData.MOMap = [ 1 1 1 2 2;... % Model ID
                 1 1 1 1 1]; % Output ID


%

     % Discrepancy group 1
DiscrepancyPriorOpts1.Name = 'Prior of sigma_1ˆ2';
DiscrepancyPriorOpts1.Marginals.Name = 'Sigma2';
DiscrepancyPriorOpts1.Marginals.Type = 'Uniform';
DiscrepancyPriorOpts1.Marginals.Parameters = [0, 1e-4];
myDiscrepancyPrior1 = uq_createInput(DiscrepancyPriorOpts1);
DiscrepancyOpts.Type = 'Gaussian';
DiscrepancyOpts.Prior = myDiscrepancyPrior1;
% D i s c r e p a n c y group 2
%DiscrepancyOpts(2).Type = 'Gaussian';
%DiscrepancyOpts(2).Prior = myDiscrepancyPrior1;
         
              
BayesOpts.Type = 'Inversion';
BayesOpts.Name = 'Bayesian multiple models';
BayesOpts.Prior = myPriorDist;
BayesOpts.ForwardModel = ForwardModels;
BayesOpts.Data = myData;
BayesOpts.Discrepancy = DiscrepancyOpts;

myBayesianAnalysis = uq_createAnalysis(BayesOpts);



uq_postProcessInversion(myBayesianAnalysis, 'burnIn', 0.3);

uq_print(myBayesianAnalysis)

uq_display(myBayesianAnalysis)


