 %% FW_2
 % Opt: optimal result
% MYX: the optimal predication
% MMSE: the maximum MSE
% Mpt: the point with MMSE
clc;
clear all;
close all;
load('FW_2_41')% input samples
load('FW_2_P_20');
%% 计算average lift
S=A(:,1:2);
Y=A(:,3);
theta=[2 2]; lob = [1e-1 1e-1]; upb = [10 10];%parameters in kriging
[dmodel, perf] = dacefit(S, Y, @regpoly0, @corrgauss, theta, lob, upb);%kriging matemodel
%% 计算power consumption
SP=B(:,1:2);
YP=B(:,3);
theta=[2 2]; lob = [1e-1 1e-1]; upb = [10 10];%parameters in kriging
[dmodelP, perfP] = dacefit(SP, YP, @regpoly0, @corrgauss, theta, lob, upb);%kriging matemodel
%% 计算lift to power consumption
[EP MSEEP] = predictor(S, dmodelP);
SLP=A(:,1:2);
YLP=A(:,3)./EP;
theta=[2 2]; lob = [1e-1 1e-1]; upb = [10 10];%parameters in kriging
[dmodelLP, perf] = dacefit(SLP, YLP, @regpoly0, @corrgauss, theta, lob, upb);%kriging matemodel
%% 计算std  estimation for st at given point
%  the=[14 30];
%  [ st_model ] = stmodel_optimal( the, dmodel );


t=1000;
X = gridsamp([2 0;25 60], 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);
Aa=[Opt, MYX, MMSE, Mpt]
% Opt: optimal result
% MYX: the optimal predication
% MMSE: the maximum MSE
% Mpt: the point with MMSE

%% 计算std  estimation of power consumption st at given point
%  the=[14 30];
%  [ st_model ] = stmodel_optimal( the, dmodel );


tP=1000;
XP = gridsamp([2 0;25 60], t);
%X = gridsamp([50 50;68 68], t);%缩小搜索范围
[YXP MSEP] = predictor(XP, dmodelP);
X1P = reshape(XP(:,1),t,t); X2P = reshape(XP(:,2),t,t);
YXP = reshape(YXP, size(X1P));
MYXP=max(YXP(:));
MMSEP=max(MSEP);
MMP= reshape(MSEP, size(X1));
[aP,bP]=find(YXP==max(YXP(:)));
[cP,dP]=find(MMP==max(MMP(:)));
tho2P=X1P(1,aP);
tho1P=X2P(bP,1);
OptP=[tho2P,tho1P];
MptP=[X2P(dP,1),X1P(1,cP)];
MMM=MM(bP,aP);
AaP=[OptP, MYXP, MMSEP, MptP]
%% 计算std  estimation of lift to LPower consumLPtion ratio st at given LPoint
%  the=[14 30];
%  [ st_model ] = stmodel_oLPtimal( the, dmodel );


tLP=1000;
XLP = gridsamp([2 0;25 60], t);
%X = gridsamLP([50 50;68 68], t);%缩小搜索范围
[YXLP MSELP] = predictor(XLP, dmodelLP);
X1LP = reshape(XLP(:,1),t,t); X2LP = reshape(XLP(:,2),t,t);
YXLP = reshape(YXLP, size(X1LP));
MYXLP=max(YXLP(:));
MMSELP=max(MSELP);
MMLP= reshape(MSELP, size(X1));
[aLP,bLP]=find(YXLP==max(YXLP(:)));
[cLP,dLP]=find(MMLP==max(MMLP(:)));
tho2LP=X1LP(1,aLP);
tho1LP=X2LP(bLP,1);
OptLP=[tho2LP,tho1LP];
MptLP=[X2LP(dLP,1),X1LP(1,cLP)];
MMM=MM(bLP,aLP);
AaLP=[OptLP, MYXLP, MMSELP, MptLP]

figure(1), mesh(X1, X2, YX)
hold on,
plot3(S((1:length(Y)),1),S((1:length(Y)),2),Y(1:length(Y)),'.k', 'MarkerSize',10)
xlabel('Flapping frequency(Hz)','fontsize',12,'fontname','Times New Roman');
ylabel('Angle of attack (degree)','fontsize',12,'fontname','Times New Roman');
zlabel('Lift','fontsize',12,'fontname','Times New Roman');
hold on


% plot3(tho1,tho2,MYX,'.r', 'MarkerSize',14)
% hold on,
% plot3(S(22,1),S(22,2),Y(22),'.r', 'MarkerSize',10)
% hold off
figure(2), mesh(X1, X2, reshape(MSE, size(X1)) )
xlabel('Flapping frequency(Hz)','fontsize',12,'fontname','Times New Roman');
ylabel('Angle of attack (degree)','fontsize',12,'fontname','Times New Roman');
zlabel('MSE','fontsize',12,'fontname','Times New Roman');
hold on
% plot3(tho1,tho2,MMM,'.r', 'MarkerSize',14)
% [emodel perf] = dacefit(S, Y, @regpoly0, @correxp, 2)%其他类型
% [y, dy] = predictor(S(1,:), dmodel)
% [y, dy, mse, dmse] = predictor([51 51], dmodel)

figure(3), mesh(X1P, X2P, YXP)
hold on,
plot3(SP((1:length(YP)),1),SP((1:length(YP)),2),YP(1:length(YP)),'.k', 'MarkerSize',10)
xlabel('Flapping frequency(Hz)','fontsize',12,'fontname','Times New Roman');
ylabel('Angle of attack (degree)','fontsize',12,'fontname','Times New Roman');
zlabel('Power consumption','fontsize',12,'fontname','Times New Roman');
hold on

figure(4), mesh(X1P, X2P, reshape(MSEP, size(X1P)) )
xlabel('Flapping frequency(Hz)','fontsize',12,'fontname','Times New Roman');
ylabel('Angle of attack (degree)','fontsize',12,'fontname','Times New Roman');
zlabel('MSE of power consumption','fontsize',12,'fontname','Times New Roman');
hold on

figure(5), mesh(X1LP, X2LP, YXLP)
hold on,
plot3(SLP((1:length(YLP)),1),SP((1:length(YLP)),2),YLP(1:length(YLP)),'.k', 'MarkerSize',10);
set(gca,'fontsize',14,'fontname','Times New Roman');
xlabel('Flapping frequency(Hz)','fontsize',14,'fontname','Times New Roman');
ylabel('Angle of attack ({\circ})','fontsize',14,'fontname','Times New Roman');
zlabel('Lift to power consumption ratio','fontsize',14,'fontname','Times New Roman');
hold on

figure(6), mesh(X1LP, X2LP, reshape(MSELP, size(X1LP)) );
set(gca,'fontsize',14,'fontname','Times New Roman');
xlabel('Flapping frequency(Hz)','fontsize',14,'fontname','Times New Roman');
ylabel('Angle of attack ({\circ})','fontsize',14,'fontname','Times New Roman');
zlabel('MSE','fontsize',14,'fontname','Times New Roman');
hold on
%% 计算可靠度
% mean2 = [15 30];   %%均值
% SIGMA2 = [0.25 0;0 1];  %%协方差矩阵,假设各维度独立，故可设为对角阵
% r2 = mvnrnd(mean2,SIGMA2,10000000);  %%最后一个参数表示生成该类100个样本点
% [YXRB2 MSERB] = predictor(r2, dmodel);
% RB1=find(YXRB2>0.24);
% RB2=length(RB1);

%% PSO search for MAX lift and Max mse in both total region and setted windows
LB=[4 -15]; %lower bounds of variables
UB=[24 60]; %upper bounds of variables
w=0.15;
m=2; % number of variables，搜索空间维数
[ best_variables,S_msemax ] = fun_VS_2( LB, UB, w, m, dmodel )
[best2 MSEbest2] = predictor(best_variables, dmodel);
[best22 MSEbest22] = predictor(Opt, dmodel);
%% PSO search for L/P and Max mse in both total region and setted windows

[ best_variablesLP,S_msemaxLP ] = fun_VS_2( LB, UB, w, m, dmodelLP )
[best2LP MSEbest2LP] = predictor(best_variablesLP, dmodelLP);
[best22LP MSEbest22LP] = predictor(OptLP, dmodelLP);
%% 计算可靠度
mean2 = [best_variablesLP(1,1) best_variablesLP(1,2)];   %%均值
SIGMA2 = [0.15 0;0 0.4];  %%协方差矩阵,假设各维度独立，故可设为对角阵
r2 = mvnrnd(mean2,SIGMA2,10000000);  %%最后一个参数表示生成该类100个样本点
[YXRB2 MSERB] = predictor(r2, dmodelLP);
RB1=find(YXRB2>0.125);%%找到他的定义
RB2LP=length(RB1);
%% 输出升力
disp(Mpt);
disp(MMSE);
disp(S_msemax);
disp(best_variables);
disp(Opt);
disp(MYX);

%% 输出L/P
disp(['MptLP = ',num2str(MptLP)]);
disp(['MMSELP = ',num2str(MMSELP)]);
disp(['S_msemaxLP = ',num2str(S_msemaxLP)]);
disp(['best_variablesLP = ',num2str(best_variablesLP)]);
disp(['OptLP = ',num2str(OptLP)]);
disp(['MYXLP = ',num2str(MYXLP)]);
disp(['RB2LP= ',num2str(RB2LP)]);
 % Opt: optimal result
% MYX: the optimal predication
% MMSE: the maximum MSE
% Mpt: the point with MaxMSE
%能耗加尾缀 P
%比值加尾缀 LP
