clearvars
rng(100,'twister')
uqlab
ModelOpts.mFile = 'RBDOLSF';
myModel = uq_createModel(ModelOpts);
InputOpts.Marginals(1).Name = 'H' ;  % significant wave height operating
InputOpts.Marginals(1).Type = 'Weibull';
InputOpts.Marginals(1).Moments = [3.19 1.691];
InputOpts.Marginals(2).Name = 'T';  % wave period operating
InputOpts.Marginals(2).Type = 'Lognormal';
InputOpts.Marginals(2).Moments = [10.21 1.838];
InputOpts.Marginals(3).Name = 'V';  % surface current operating
InputOpts.Marginals(3).Type = 'Weibull';
InputOpts.Marginals(3).Moments = [0.284 0.159];
InputOpts.Marginals(4).Name = 'He' ;  % significant wave height extreme
InputOpts.Marginals(4).Type = 'Weibull';
InputOpts.Marginals(4).Moments = [17.20 3.199];
InputOpts.Marginals(5).Name = 'Te';  % wave period extreme
InputOpts.Marginals(5).Type = 'Lognormal';
InputOpts.Marginals(5).Moments = [18 5.4];
InputOpts.Marginals(6).Name = 'Ve';  % surface current extreme
InputOpts.Marginals(6).Type = 'Weibull';
InputOpts.Marginals(6).Moments = [2.48 0.2604];
InputOpts.Copula(1).Type = 'DVine'; % dependece for operating condition
InputOpts.Copula(1).Families = {'Clayton', 'Clayton', 'Independent'};
InputOpts.Copula(1).Parameters = {0.04,0.05,[]};
InputOpts.Copula(1).Rotations = [0 270 0];
InputOpts.Copula(1).Variables = [1 2 3];
InputOpts.Copula(1).Structure = 1:3;
InputOpts.Copula(2).Type = 'DVine'; %dependence for extreme condition
InputOpts.Copula(2).Families = {'Gumbel', 'Gumbel', 'Gumbel'};
InputOpts.Copula(2).Parameters = {3.31, 2.29,2.31};
InputOpts.Copula(2).Rotations = [0 270 180];
InputOpts.Copula(2).Variables = [4 5 6];
InputOpts.Copula(2).Structure = 1:3;
myInput = uq_createInput(InputOpts);
RBDOOpts.Input.EnvVar = myInput;
RBDOOpts.Cost.mString = '(3.142.*((((0.22225.*X(:,1))+(X(:,1).^2)).*650)+(((0.22225.*X(:,2))+(X(:,2).^2)).*500)+(((0.22225.*X(:,3))+(X(:,3).^2)).*805)))';
RBDOOpts.LimitState.Model = myModel;
RBDOOpts.Input.DesVar(1).Name = 't1';
RBDOOpts.Input.DesVar(1).Type = 'Gaussian';
RBDOOpts.Input.DesVar(1).Std = 0.1;
RBDOOpts.Input.DesVar(2).Name = 't2';
RBDOOpts.Input.DesVar(2).Type = 'Gaussian';
RBDOOpts.Input.DesVar(2).Std = 0.1;
RBDOOpts.Input.DesVar(3).Name = 't3';
RBDOOpts.Input.DesVar(3).Type = 'Gaussian';
RBDOOpts.Input.DesVar(3).Std = 0.1;
RBDOOpts.Optim.Bounds = [0.045 0.025 0.017; 0.065 0.045 0.025];
RBDOOpts.Optim.StartingPoint = [0.05 0.03 0.02];
RBDOOpts.TargetBeta = 3;
RBDOOpts.Optim.MaxIter = 500;
RBDOOpts.Type = 'RBDO';
RBDOOpts.Method = 'Two-Level';
RBDOOpts.Reliability.Method = 'Subset' ;
RBDOOpts.Reliability.SubsetSim.p0 = 0.2 ;
RBDOOpts.Reliability.Simulation.BatchSize = 1000 ;
RBDO_PCKOpts = RBDOOpts;
RBDO_PCKOpts.Metamodel.Type = 'PCK';
PCKOpts.ExpDesign.NSamples =25;
PCKOpts.ExpDesign.Sampling = 'LHS';
RBDO_PCKOpts.Metamodel.PCK = PCKOpts;
%RBDO_PCKOpts.Metamodel.Enrichment.MaxAddedED = 1000;
RBDO_PCKOpts.Constraints.Enrichment.MaxAdded = 1000;
RBDO_PCKOpts.Metamodel.Enrichment.LearningFunction='U';
RBDO_PCKOpts.Metamodel.Enrichment.SampleSize=100000;
RBDO_PCKOpts.Optim.Method = 'HGA';
myRBDO_PCK= uq_createAnalysis(RBDO_PCKOpts);
uq_print(myRBDO_PCK)
uq_display(myRBDO_PCK)
uq_print(myInput);
uq_display(myInput)