uqlab

% LR
InputOpts.Marginals(1).Name = 'LR';
InputOpts.Marginals(1).Type = 'Uniform';
InputOpts.Marginals(1).Moments  = [3 20];

% E2
InputOpts.Marginals(2).Name = 'RT';
InputOpts.Marginals(2).Type = 'Uniform';
InputOpts.Marginals(2).Moments  = [20.1 30];

% E3
InputOpts.Marginals(3).Name = 'A1';
InputOpts.Marginals(3).Type = 'Uniform';
InputOpts.Marginals(3).Moments  = [0 90];

InputOpts.Marginals(4).Name = 'A2';
InputOpts.Marginals(4).Type = 'Uniform';
InputOpts.Marginals(4).Moments  = [0 90];

InputOpts.Marginals(5).Name = 'E1';
InputOpts.Marginals(5).Type = 'Gaussian';
InputOpts.Marginals(5).Moments = [102000, 10200];
% InputOpts.Marginals(5).Cov = 0.1;
myInput = uq_createInput(InputOpts);

MetaOpts.Type = 'Metamodel';
MetaOpts.MetaType = 'PCE';

MetaOpts.ExpDesign.X = X;
MetaOpts.ExpDesign.Y = Y;

MetaOpts.Degree = 2:10;

MetaOpts.ValidationSet.X = Xval;
MetaOpts.ValidationSet.Y = Yval;

myPCE = uq_createModel(MetaOpts);

uq_print(myPCE)

YPCE = uq_evalModel(myPCE,Xval);

uq_figure

cmap = uq_colorOrder(2);
uq_histogram(Yval, 'FaceColor', cmap(1,:))
hold on
uq_histogram(YPCE, 'FaceColor', cmap(2,:))
hold off

xlabel('Critical buckling pressure (MPa)')
ylabel('Counts')
uq_legend(...
    {'True model response', 'PCE prediction'},...
    'Location', 'northwest')

uq_figure

uq_plot(Yval, YPCE, '+')
hold on
uq_plot([min(Yval) max(Yval)], [min(Yval) max(Yval)], 'k')
hold off

axis equal
axis([min(Yval) max(Yval) min(Yval) max(Yval)])

xlabel('$\mathrm{Y_{true}}$')
ylabel('$\mathrm{Y_{PCE}}$')

fprintf('PCE metamodel validation error: %5.4e\n', myPCE.Error.Val)
fprintf('PCE metamodel LOO error:        %5.4e\n', myPCE.Error.LOO)

% 
function Y = LSF(X) 
    Y = uq_evalModel(myPCE,X)-5;
end
% -------------------------------------------

% 
% rng(100,'twister')
% uqlab  
% % 
% % myModel = uq_createModel(ModelOpts);
% myModel = uq_createModel(ModelOpts);
% Iopts.Marginals(1).Name = 'E1' ;
% Iopts.Marginals(1).Type = 'Gaussian' ;
% Iopts.Marginals(1).Moments = [102000 10200] ;
% myInput = uq_createInput(Iopts) ;
% 
% % Constraint function
% RBDOOpts.LimitState.Model = myModel ;
% % RBDOOpts.LimitState.Model = myPCEModel;
% RBDOOpts.Cost.mString = '3.1415.*100.*X(:,1).*(200 + 100./X(:,2)).*(100./X(:,2))';
% % Design variables
% RBDOOpts.Input.DesVar(1).Name = 'LR';
% RBDOOpts.Input.DesVar(1).Type = 'constant';
% RBDOOpts.Input.DesVar(2).Name = 'RT';
% RBDOOpts.Input.DesVar(2).Type = 'constant';
% RBDOOpts.Input.DesVar(3).Name = 'A1';
% RBDOOpts.Input.DesVar(3).Type = 'constant';
% RBDOOpts.Input.DesVar(4).Name = 'A2';
% RBDOOpts.Input.DesVar(4).Type = 'constant';
% RBDOOpts.Input.EnvVar = myInput;
% 
% RBDOOpts.Optim.Bounds = [3 20.1 0 0; 20 30 90 90 ];
% RBDOOpts.Optim.StartingPoint = [4 25 45 45];
% RBDOOpts.TargetBeta = 3;
% RBDOOpts.Optim.MaxIter = 500;
% 
% RBDOOpts.Type = 'RBDO';
% RBDOOpts.Method = 'QMC';
% 
% RBDO_PCEOpts = RBDOOpts;
% 
% RBDO_PCEOpts.Metamodel.Type = 'PCE';
% 
% PCEOpts.Degree = 2:10 ;
% PCEOpts.ExpDesign.X = X;
% PCEOpts.ExpDesign.Y = Y;
% RBDO_PCEOpts.Metamodel.PCE = PCEOpts;
% 
% RBDO_PCEOpts.Metamodel.Enrichment.Convergence = {'stopSign','stopDf'};
% RBDO_PCEOpts.Metamodel.Enrichment.ConvThreshold = 0.001;
% RBDO_PCEOpts.Metamodel.Enrichment.MaxAdded = 50;
% RBDO_PCEOpts.Metamodel.Enrichment.MOStrategy = 'mean';
% 
% myRBDO_PCE = uq_createAnalysis(RBDO_PCEOpts);
% 
% uq_print(myRBDO_PCE)
% uq_display(myRBDO_PCE)