% Example in Doc sec 2.2.1 to 2.2.5 modified by Sec. 2.5  and Sec 2.3
% afterwards
fprintf(' Example in Doc, sec 2.2.1 to 2.2.5; modified by Sec. 2.5 and also Sec. 2.3\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)
PriorOpts.Marginals(6).Name = 'P';
PriorOpts.Marginals(6).Type = 'Constant';
PriorOpts.Marginals(6).Parameters = 0.05; % ( kN )
myPriorDist = uq_createInput(PriorOpts);

ModelOpts1.Name = 'Model: Beam quarter and midspan deflection';
ModelOpts1.mFile = 'uq_SimplySupportedBeamTwo';
ForwardModels(1).Model = uq_createModel(ModelOpts1);
ForwardModels(1).PMap = [1 2 3 4 5];
ModelOpts2.Name = 'Beam elongation';
ModelOpts2.mString = 'X(:,5).*X(:,3)./(X(:,1).*X(:,2).*X(:,4))';
ForwardModels(2).Model = uq_createModel(ModelOpts2);
ForwardModels(2).PMap = [1 2 3 4 6];


U = [0.485; 0.466; 0.486]/1000; % (m)
% Data group 1

V = [[ 8.98; 8.66; 8.85; 9.19; 8.64],... % L /4 (m)
[12.84; 13.12; 12.13; 12.19; 12.67]]/1000; % L /2 (m)

myData(1).y = V
myData(1).Name = 'Data: Beam quarter and midspan deflection'
myData(1).MOMap = [ 1 1;... % Model ID
                   1 2] ; % Output ID


%  Data group 2
myData(2).y = U;
myData(2).Name = 'Beam elongation';
myData(2).MOMap = [ 2;... % Model ID 
     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(1).Type = 'Gaussian';
DiscrepancyOpts(1).Prior = myDiscrepancyPrior1;
% D i s c r e p a n c y group 2
DiscrepancyOpts(2).Type = 'Gaussian';
DiscrepancyOpts(2).Parameters = 1e-6; % knownd i s c r e p a n c yv a r i a n c e
         
              
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)

uq_postProcessInversion(myBayesianAnalysis,'priorPredictive',1000);



uq_print(myBayesianAnalysis)
uq_display(myBayesianAnalysis,'predDist',1)

