% Example in Doc sec 2.2.1 to 2.2.5 modified folliowing 2.3.1 to 2.3.4.2, oneDiscrepany"
% 
fprintf(' Example in Doc, sec 2.2.1 to 2.2.5. modified folliowing\n')
fprintf('2.3.1 to 2.3.4.2 one Discrepandy ,\n')
clear all 
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);

ModelOpts.Name = 'Forward model';
ModelOpts.mFile = 'uq_SimplySupportedBeamTwo';
myForwardModel = uq_createModel(ModelOpts);

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.MOMap = [ 1 1 1 1 1  1 1 1 1 1 ;... % Model ID
                 1 1 1 1 1  2 2 2 2 2 ]; % Output ID
myData.MOMapNotUnique = [1];
VmomapPre =  cat (1,V(:,1),V(:,2))
Vmomap = transpose(VmomapPre)

myData.Name = 'Beam quarter and midspan deflection'
myData.y = Vmomap

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;

BayesOpts.Type = 'Inversion';
BayesOpts.Data = myData;
BayesOpts.Name = 'Testing using one Distribution';
BayesOpts.Prior = myPriorDist;
%BayesOpts.ForwardModel = myForwardModel;
BayesOpts.Discrepancy = DiscrepancyOpts;

% 
uq_listModels

myBayesianAnalysis = uq_createAnalysis(BayesOpts);
%uq_postProcessInversion(myBayesianAnalysis,'priorPredictive',1000);

%

uq_print(myBayesianAnalysis)
uq_display(myBayesianAnalysis)



%
%uq_postProcessInversion(myBayesianAnalysis,'priorPredictive',1000);



%uq_display(myBayesianAnalysis,'predDist',1)
