 %% VS_P
 % Theta_Emax:optimal result
% Emax:optimal predication
% MSEo:MSE of optimal
% MSEMAX:Maximum MSE
% S_msemax :sample with MSEMAX
% Theta_MSEmax:sample with MSEMAX in the set window
% st_model:st at optimal result

close all; clear all; clc

load('FW_3_26')% input samples
load('FW_3_P_30');
%% 计算average lift
S=A(:,1:3);
Y=A(:,4);
theta=[2 2 2]; lob = [1e-10 1e-10 1e-10]; upb = [10 10 10];%parameters in kriging
[dmodel, perf] = dacefit(S, Y, @regpoly0, @corrgauss, theta, lob, upb);
%% 计算power consumption
SP=B(:,1:3);
YP=B(:,4)+0.3;
[dmodelP, perfP] = dacefit(SP, YP, @regpoly0, @corrgauss, theta, lob, upb);%kriging matemodel
%% 计算lift to power consumption
[EP MSEEP] = predictor(S, dmodelP);
SLP=A(:,1:3);
YLP=A(:,4)./EP;
[dmodelLP, perf] = dacefit(SLP, YLP, @regpoly0, @corrgauss, theta, lob, upb);%kriging matemodel
%% 计算std
% the=[10	2 4	
% 14	1.7 6	
% 17	3.1 17	
% 11	4.2 18	
% 12	2.2 27];
% the=[81.75607254	54.94864698	60.53461817
% 82.5769001	56.26391853	64.00934361
% 78.98268277	54.77895887	62.15143848
% 77.59011286	56.75605019	62.17085242
% 71.46860829	58.39279628	56.69781872
% 71.46168926	58.94953005	53.23838895
% 72.9343253	57.05069709	62.91978814
% 74.20237589	52.21711027	68.74948382
% 71.59559549	57.09138258	59.80612652
% 64.06502182	58.37746173	57.80382705];
% [ st_model ] = stmodel_optimal( the(5,:), dmodel );



N=10000;
D=3;
t=80;
LB = [4,-15,0];
UB = [20,60,6];
% [ ss ] = sample_uniform( LB, UB, N, D );
[ ss ] = gridsamp([4 -15 0;20 60 6], t);
[YX MSE] = predictor(ss, dmodel);

[a]=find(YX==max(YX(:)));
Emax=max(YX(:))
Theta_Emax=[ss(a,:)]
% [YXo MSEo] = predictor(Theta_Emax, dmodel);
[z1,z2,MSEo,z3] = predictor(Theta_Emax, dmodel);
MSEMAX=max(MSE);
[b]=find(MSE==max(MSE(:)));
Theta_MSEmax=[ss(b,:)];
% t=800;
% X = gridsamp([50 50;68 68], t);
% %X = gridsamp([50 50;68 68], t);%缩小搜索范围
% [YX MSE] = predictor(X, dmodel);
% X1 = reshape(X(:,1),t,t); X2 = reshape(X(:,2),t,t);
% YX = reshape(YX, size(X1));
% MYX=max(YX(:));
% MMSE=max(MSE);
% MM= reshape(MSE, size(X1));
% [a,b]=find(YX==max(YX(:)));
% [c,d]=find(MM==max(MM(:)));
% tho2=X1(1,a);
% tho1=X2(b,1);
% Opt=[tho2,tho1];
% Mpt=[X2(d,1),X1(1,c)];
% MMM=MM(b,a);
%% 计算std  estimation of lift to Power consumption ratio st at given Point
[YXLP MSELP] = predictor(ss, dmodelLP);

[aLP]=find(YXLP==max(YXLP(:)));
EmaxLP=max(YXLP(:))
Theta_EmaxLP=[ss(aLP,:)]
% [YXo MSEo] = predictor(Theta_Emax, dmodel);
[z1LP,z2LP,MSEoLP,z3LP] = predictor(Theta_EmaxLP, dmodelLP);
MSEMAXLP=max(MSELP);
[bLP]=find(MSELP==max(MSELP(:)));
Theta_MSEmaxLP=[ss(bLP,:)];

%% PSO search for MAX eigenvalue and Max mse in both total region and setted windows
% LB = [4 -15 0];
% UB = [30,60,6];
w=0.15;
m=3; % number of variables，搜索空间维数
[ best_variables,S_msemax ] = fun_VS_4( LB, UB, w, m, dmodel );
[best3 MSEbest3] = predictor(best_variables, dmodel);
[best33 MSEbest33] = predictor(Theta_Emax, dmodel);
%% PSO search for L/P in both total region and setted windows
% LB = [4 -15 0];
% UB = [30,60,6];
w=0.15;
m=3; % number of variables，搜索空间维数
[ best_variablesLP,S_msemaxLP ] = fun_VS_4( LB, UB, w, m, dmodelLP );
[best3LP MSEbest3LP] = predictor(best_variablesLP, dmodelLP);
[best33LP MSEbest33LP] = predictor(Theta_EmaxLP, dmodelLP);
%% 计算可靠度
mean3 = [Theta_EmaxLP(1,1) Theta_EmaxLP(1,2) Theta_EmaxLP(1,3)];   %%均值
SIGMA3 = [0.15 0 0; 0 0.4 0; 0 0 0.05];  %%协方差矩阵,假设各维度独立，故可设为对角阵
r3 = mvnrnd(mean3,SIGMA3,1000000);  %%最后一个参数表示生成该类100个样本点
[YXRB3 MSERB] = predictor(r3, dmodelLP);
RB4=find(YXRB3>0.12);
RB3LP=length(RB4);

disp(Theta_MSEmax); 
disp(MSEMAX); 
disp(S_msemax); 
disp(Theta_Emax); 
disp(MSEo); 
disp(Emax); 
disp(best_variables);


disp(['Theta_MSEmaxLP = ',num2str(Theta_MSEmaxLP)]); 
disp(['MSEMAXLP = ',num2str(MSEMAXLP)]); ; 
disp(['S_msemaxLP = ',num2str(S_msemaxLP)]); 
disp(['Theta_EmaxLP = ',num2str(Theta_EmaxLP)]); ; 
disp(['MSEoLP = ',num2str(MSEoLP)]);  
disp(['EmaxLP = ',num2str(EmaxLP)]);  
disp(['best_variablesLP = ',num2str(best_variablesLP)]); ;
disp(['RB3LP = ',num2str(RB3LP)]); ;
 % Theta_Emax:optimal result
% Emax:optimal predication
% MSEo:MSE of optimal
% MSEMAX:Maximum MSE
% S_msemax :sample with MSEMAX
% Theta_MSEmax:sample with MSEMAX in the set window
% st_model:st at optimal result
