%% RELIABILITY ANALYSIS

% Initialization
clc;
uqlab;
rng(1,'twister');

% Define the model type (UQLink) and name
ZS.Type = 'UQLink';

% Define the executable ZSOIL command and path
ZS.Command = '"C:\Program Files\ZSoil\ZSoil 2020 v20.06 x64\Z_Soil.exe" "C:\Users\gjdes.KH2-AS-MAR19\Documents\MATLAB\UQLab\UQDay\COL_NEUMANN_K0\RF\IS_RF.inp" /E';
ZS.ExecutionPath = 'C:\Users\gjdes.KH2-AS-MAR19\Documents\MATLAB\UQLab\UQDay\COL_NEUMANN_K0\RF';

% Define the template input file name
ZS.Template = 'IS_RF.inp.tpl.inp';

% Parser user to retrieve results
ZS.Output.Parser = 'readSF';
ZS.Output.FileName = {"IS_RF.his"};

% Create the UQLink wrapper
myUQLinkModel = uq_createModel(ZS);

% Forward model
mOpts.mFile = 'RFmodel';
mOpts.isVectorized = false;
mOpts.Parameters.z = [0.125/2:0.125:2.5]; % Discretization of the RF.
mOpts.Parameters.mu = 25; % Mean.
mOpts.Parameters.sigma = 2.5; % STD.
mOpts.Parameters.theta = 0.5; % Correlation length.
mOpts.Parameters.myUQLinkModel = myUQLinkModel;
myModel = uq_createModel(mOpts);

% Probabilistic model
for el = 1:20
    IOpts.Marginals(el).Name = 'cu';              
    IOpts.Marginals(el).Type = 'Gaussian';
    IOpts.Marginals(el).Moments = [0 1];
end
IOpts.Marginals(21).Name = '$\alpha$';              
IOpts.Marginals(21).Type = 'Uniform';
IOpts.Marginals(21).Parameters = [0 5];
myInput = uq_createInput(IOpts);

%% RELIABILITY

AKOptions.Type = 'Reliability';
AKOptions.Method = 'AKMCS';
AKOptions.Input = myInput;
AKOptions.Model = myModel;
AKOptions.AKMCS.MaxAddedED = 1;
AKOptions.LimitState.Threshold = 1;
AKOptions.LimitState.CompOp = '<';
AKAnalysis = uq_createAnalysis(AKOptions);

% Results AK
uq_print(AKAnalysis);
uq_display(AKAnalysis);
