%% 1 - INITIALIZE UQLAB
%
% Clear all variables from the workspace 
% and initialize the UQLab framework:
clearvars
uqlab

%% 2 - RETRIEVE DATA SETS
%
% The experimental design and the validation basis are stored
% in two separate files in the following location:
FILELOCATION = fullfile(uq_rootPath, 'Examples', 'SimpleDataSets',...
    'Cang');

%%
% Read the experimental design data set file and store the contents 
% in matrices:
load(fullfile(FILELOCATION,'Cang_DOE6.mat'), 'X', 'Y');

%%
% Read the validation basis data set file and store the contents
% in matrices:
load(fullfile(FILELOCATION,'Cang_VAL2.mat'), 'Xval', 'Yval');

%% 3 - INPUT MODEL
%% KS
% InputOpts.Marginals(1).Name = 'a';
% InputOpts.Marginals(1).Type = 'Gumbel'; %Uniform\Gaussian\Lognormal\Gumbel,etc
% InputOpts.Marginals(1).Moments = [6.808e+01 2.909e+01];
% 
% InputOpts.Marginals(2).Name = 'b';
% InputOpts.Marginals(2).Type = 'Gumbel' ;
% InputOpts.Marginals(2).Moments = [8.120e+00 1.391e+00];
% 
% InputOpts.Marginals(3).Name = 'c';
% InputOpts.Marginals(3).Type = 'Laplace';
% InputOpts.Marginals(3).Moments = [4.965e-02 9.829e-04];
%% BIC/auto
%%%%Error:Calculation of parameters from moments is not defined for marginal type: Rayleigh !
InputOpts.Marginals(1).Name = 'a';
InputOpts.Marginals(1).Type = 'Rayleigh '; %Uniform\Gaussian\Lognormal\Gumbel,etc
InputOpts.Marginals(1).Moments = [ 6.392e+01, 3.341e+01];

InputOpts.Marginals(2).Name = 'b';
InputOpts.Marginals(2).Type = 'Uniform' ;
InputOpts.Marginals(2).Moments = [8.264e+00, 1.490e+00];

InputOpts.Marginals(3).Name = 'c';
InputOpts.Marginals(3).Type = 'Laplace';
InputOpts.Marginals(3).Moments = [4.965e-02, 9.829e-04];

%% BIC/ks（kernel smoothing）
%%%%Error: the specified input does not seem to be either a string nor a recognized object!
% InputOpts.Marginals(1).Name = 'a';
% InputOpts.Marginals(1).Type = 'ks ';    %Uniform\Gaussian\Lognormal\Gumbel,etc
% InputOpts.Marginals(1).Moments = [ 6.750e+01 2.786e+01];
% 
% InputOpts.Marginals(2).Name = 'b';
% InputOpts.Marginals(2).Type = 'Uniform' ;
% InputOpts.Marginals(2).Moments = [8.264e+00, 1.490e+00];
% 
% InputOpts.Marginals(3).Name = 'c';
% InputOpts.Marginals(3).Type = 'Laplace';
% InputOpts.Marginals(3).Moments = [4.965e-02, 9.829e-04];
%%
% % InputOpts.Marginals(4).Name = 'F4';
% % InputOpts.Marginals(4).Type = 'Weibull' ;
% % InputOpts.Marginals(4).Moments = [2.106e-01 7.730e-02];
% % 
% % InputOpts.Marginals(5).Name = 'F5';
% % InputOpts.Marginals(5).Type = 'Logistic';
% % InputOpts.Marginals(5).Moments = [3.730e+01 5.456e+01];
% % 
% % InputOpts.Marginals(6).Name = 'F6';
% % InputOpts.Marginals(6).Type = 'Logistic';
% % InputOpts.Marginals(6).Moments = [5.412e-01 1.036e-02];
%% 
myInput = uq_createInput(InputOpts);
uq_print(myInput)
uq_display(myInput)      
% %uq_display(myInput, [1,2], 'marginals');


%%
X1=X(:,1);
X2=X(:,2);
X3=X(:,3);
% % X4=X(:,4);
% % X5=X(:,5);
% % X6=X(:,6);
%% Specification of inference options for each marginal
% InputOpts.Marginals.Type = 'ks';  %kernel smoothing (ks)【不适用于任何已知分布类型】\auto
% InputOpts.Inference.Criterion = 'BIC';  %KS\BIC
% 
% InputOpts.Marginals(1).Inference.Data = X1;
% InputOpts.Marginals(2).Inference.Data = X2;
% InputOpts.Marginals(3).Inference.Data = X3;
% % InputOpts.Marginals(4).Inference.Data = X4;
% % InputOpts.Marginals(5).Inference.Data = X5;
% % InputOpts.Marginals(6).Inference.Data = X6;

% InputHat = uq_createInput(InputOpts);
% uq_print(InputHat);

%%
% X1=X(:,1);
% X2=X(:,2);
% X3=X(:,3);
% X4=X(:,4);
% X5=X(:,5);
% subplot(3,3,1);
% plotmatrix(X1,X2);
% subplot(3,3,2);
% plotmatrix(X1,X3);
% subplot(3,3,3);
% plotmatrix(X1,X4);
% subplot(3,3,4);
% plotmatrix(X1,X5);
% subplot(3,3,5);
% plotmatrix(X1,Y);
% subplot(3,3,6);
% plotmatrix(X2,Y);
% subplot(3,3,7);
% plotmatrix(X3,Y);
% subplot(3,3,8);
% plotmatrix(X4,Y);
% subplot(3,3,9);
% plotmatrix(X5,Y);

%%
%X_MC = uq_getSample(35,'LHS');

%% 4 - POLYNOMIAL CHAOS EXPANSION (PCE) METAMODEL
MetaOpts.Type = 'Metamodel';
MetaOpts.MetaType = 'PCE';
MetaOpts.ExpDesign.X = X;
MetaOpts.ExpDesign.Y = Y;
MetaOpts.Degree = 1:10;
MetaOpts.Method = 'LARS' ;
%MetaOpts.Method = 'OMP';
%MetaOpts.Method = 'Quadrature';%算不动,不适合，误差极大
%MetaOpts.Method = 'OLS';%算不动
%MetaOpts.Method = 'BCS';
%MetaOpts.Method = 'SP';
%MetaOpts.TruncOptions.qNorm = 0.5;
PCEOpts.TruncOptions.qNorm = [0.3 0.5 0.75 1.0 ];
PCEOpts.TruncOptions.MaxInteraction =2;
MetaOpts.ValidationSet.X = Xval;
MetaOpts.ValidationSet.Y = Yval;
myPCE = uq_createModel(MetaOpts);
uq_print(myPCE)
%%
%%
% SSEOpts.Type = 'Metamodel';
% SSEOpts.MetaType = 'SSE';
% SSEOpts.ExpDesign.X = X;
% SSEOpts.ExpDesign.Y = Y;
% SSEOpts.ValidationSet.X = Xval;
% SSEOpts.ValidationSet.Y = Yval;
% SSEOpts.ExpOptions.Type = 'Metamodel';
% SSEOpts.ExpOptions.MetaType = 'PCE';
% SSEOpts.ExpOptions.Degree = 0:5;
% 
% 
% MetaOpts.Method = 'LARS' ;
% MetaOpts.TruncOptions.qNorm = 0.75;
% SSEOpts.Refine.NExp = 3;
% SSEOpts.PostProcessing.outputMoments = true;
% mySSE = uq_createModel(SSEOpts);
% uq_print(mySSE)
% %myFigures = uq_display(mySSE,1);

%% Kriging
% MetaOpts.Type = 'Metamodel';
% MetaOpts.MetaType = 'Kriging';
% MetaOpts.ExpDesign.X = X;
% MetaOpts.ExpDesign.Y = Y;
% MetaOpts.EstimMethod = 'ML';
% MetaOpts.Optim.Method = 'CMAES';
% MetaOpts.ValidationSet.X = Xval;
% MetaOpts.ValidationSet.Y = Yval;
% mySurrogate = uq_createModel(MetaOpts);
% uq_print(mySurrogate)

%% PCK
% MetaOpts.Type = 'Metamodel';
% MetaOpts.MetaType = 'PCK';
% MetaOpts.Mode = 'sequential';
% MetaOpts.ExpDesign.X = X;
% MetaOpts.ExpDesign.Y = Y;
% MetaOpts.Kriging.Optim.Method = 'CMAES';
% MetaOpts.PCE.Degree = 1:5;
% MetaOpts.ValidationSet.X = Xval;
% MetaOpts.ValidationSet.Y = Yval;
% mySurrogate = uq_createModel(MetaOpts);
% uq_print(mySurrogate)

%% SVM
% MetaOpts.Type = 'Metamodel';
% MetaOpts.MetaType = 'SVR';
% 
% MetaOpts.ExpDesign.X = X;
% MetaOpts.ExpDesign.Y = Y;
% MetaOpts.Loss = 'l1-eps';
% MetaOpts.Kernel.Isotropic = 0;
% MetaOpts.Optim.Method = 'BFGS'; %HGA
% %MetaOpts.QPSolver = QPSolver;
% MetaOpts.Kernel.Family = 'Matern-5_2'; %'Gaussian';
% MetaOpts.ValidationSet.X = Xval;
% MetaOpts.ValidationSet.Y = Yval;
% mySurrogate = uq_createModel(MetaOpts);
% uq_print(mySurrogate)
    
%% 5 - VALIDATION

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('$\mathrm{Y}$')
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)