%% Bayesian Study on my Kriging Model
% Sophy Oliver Sep-2019

clear all

%% Step 1: Start the UQLab framework

rng(100,'twister')
uqlab;

%% Step 2: Create the Model

% This model is a Kriging meta model of a more expensive model. 
% The more expensive model consists of an ocean biogeochemical model which
% attempts to model and replicate the chemistry of the real ocean using 
% PDEs and 6 input parameters, then returns a single value indicating the
% difference (or discrepency) between the modelled ocean and the real
% ocean. So we want to find the 6 parameters which result in the lowest
% output from the kriging meta model.

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

%% Step 3: Create prior parameter info

% The model takes in 6 parameter values, each with a uniform distribution
% between a lower and upper bound

a = [150, 4, 0.0001, 0.1, 0, 0.4];
b = [200, 48, 0.5, 4, 10, 1.8];
name = {'ro2ut', 'Ic', 'KPHY', 'muZOO', 'kZOO', 'b'};
for i=1:6
    PriorOpts.Marginals(i).Name = name{i};
    PriorOpts.Marginals(i).Type = 'Uniform';
    PriorOpts.Marginals(i).Parameters = [a(i) b(i)];
end
myPriorDist = uq_createInput(PriorOpts);

%% Step 4: Create prior discrepancies info

% Not really sure if my discrepancies are known or not...
% Should this be the discrepency between the real model and the meta model?

DiscrepancyOptsKnown.Type = 'Gaussian';
DiscrepancyOptsKnown.Parameters = 1;

%% Step 5: Data

% For the original expensive model the output is never less than zero, and
% we want to find the 6 parameters which result in zero (as the the output
% value is the difference bwteen the modelled ocean and the real ocean).

% For the Kriging meta model, as it's not an exact replica of the true model
% it sometimes gives a result below zero. So here below I'm not sure if the
% target should be zero, or the minimum value returned by the meta model
% (~-50).

myData.Name = 'Target';
myData.y = -50;

%% Step 6: Solver

Solver.Type = 'MCMC';
Solver.MCMC.Sampler = 'AIES';
Solver.MCMC.NChains = 50;
Solver.MCMC.Steps = 200;

Solver.MCMC.Visualize.Parameters = 1:6;
Solver.MCMC.Visualize.Interval = 20;

%% Step 7: Run Bayes Analysis

BayesOpts.Type = 'Inversion';
BayesOpts.Data = myData;
BayesOpts.Discrepancy = DiscrepancyOptsKnown;
BayesOpts.Solver = Solver;

myBayesianAnalysis = uq_createAnalysis(BayesOpts);
uq_print(myBayesianAnalysis)
uq_display(myBayesianAnalysis)
