%% 加载神经网络的数据集 每行一个样本 每列一个特征 输入P 输出T
%样本数据就是前面问题描述中列出的数据
load input.mat % 【需要修改】导入输入数据集
load output.mat %【需要修改】导入输出数据集
P_train=input(1:75,1:5); %【需要修改】第1-75行作为训练数据，输入数据集第1-5列是输入的特征
T_train=output(1:75,1); %【需要修改】第1-75行作为训练数据，输出数据集第1列作为最终预测的特征
P_test=input(76:100,1:5); %【需要修改】第76-100行作为测试数据，输入数据集第1-5列是输入的特征
T_test=output(76:100,1); %【需要修改】第76-100行作为测试数据，输出数据集第1列作为最终预测的特征
%% 归一化
% [p_train,ps_train]=mapminmax(P_train',0,1);
% [t_train,ps_output]=mapminmax(T_train',0,1);
% [t_test,ts_output]=mapminmax(T_test',0,1);
% p_test=mapminmax('apply',P_test',ps_train);
%% 获取输入层节点、输出层节点个数
% 初始隐层神经元个数
inputnum=size(P_train,2);
hiddennum=2*inputnum+1;
outputnum=size(T_train,2);
w1num=inputnum*hiddennum; % 输入层到隐层的权值个数
w2num=outputnum*hiddennum;% 隐层到输出层的权值个数
N=w1num+hiddennum+w2num+outputnum; %待优化的变量的个数
% warning('off')
%% 定义遗传算法参数
NIND=50; %个体数目
MAXGEN=50; %最大遗传代数
PRECI=10; %变量的二进制位数
GGAP=0.93; %代沟
px=0.84; %交叉概率
pm=0.02; %变异概率
trace=zeros(N+1,MAXGEN);%寻优结果的初始值（初始化操作）
FieldD=[repmat(PRECI,1,N);repmat([-0.5;0.5],1,N);repmat([1;0;1;1],1,N)];%区域描述器，这一项主要是用来控制个体由2进制变成10进制的规则（个人理解就是数据由2进制转化为10进制并对样本进行归一化的过程）
Chrom=crtbp(NIND,PRECI*N);%初始种群 这里会生成一个40 x（10*N）的矩阵，即每一个种群中都包含了N个2进制的数用来表示N个个体
%% 优化
gen=0; %代计数器
X=bs2rv(Chrom,FieldD); %计算初始种群的十进制转换（根据上述规则将X的2进制转化为10进制）
ObjV=Objfun(X,P_train,T_train,hiddennum,P_test,T_test); %计算目标函数值（将X带入适应度函数进行判断）
while gen<MAXGEN
fprintf('%d\n',gen)
FitnV=ranking(ObjV); %分配适应度值
SelCh=select('sus',Chrom,FitnV,GGAP); %选择
SelCh=recombin('xovsp',SelCh,px); %重组
SelCh=mut(SelCh,pm); %变异
X=bs2rv(SelCh,FieldD); %子代个体的十进制转换
ObjVSel=Objfun(X,P_train,T_train,hiddennum,P_test,T_test); %计算子代的目标函数值（经过交叉变异重组后再次计算适应度）
[Chrom,ObjV]=reins(Chrom,SelCh,1,1,ObjV,ObjVSel); %重插入子代到父代，得到新种群
X=bs2rv(Chrom,FieldD);%上述while中的计算都是用的2进制编码形式，接着将新的X转换为10进制
gen=gen+1; %代计数器增加
%获取每代的最优解及其序号，Y为最优解个体的值,I为个体的序号
[Y,I]=min(ObjV);
trace(1:N,gen)=X(I,:); %记下每代的最优值
trace(end,gen)=Y; %记下每代的最优值
end
%% 画进化图
figure(1);
plot(1:MAXGEN,trace(end,:));
grid on
xlabel('遗传代数')
ylabel('误差的变化')
title('进化过程')
bestX=trace(1:end-1,end);
bestErr=trace(end,end);
[~,bestT_sim]=BpFunction(bestX,P_train,T_train,hiddennum,P_test,T_test);
fprintf(['最优初始权值和阈值:\nX=',num2str(bestX'),'\n最小误差err=',num2str(bestErr),'\n'])
figure
plot(T_test(:,1),'--','linewidth',1.2)
hold on
plot(bestT_sim(:,1),'linewidth',1.2)
legend('期望值','预测值')
xlabel('测试样本编号'),ylabel('指标值')
title('GA-BP测试集预测值和期望值的对比')
set(gca,'fontsize',12)
figure
plot(T_test(:,1)-bestT_sim(:,1),'r-','linewidth',1.2)
xlabel('测试样本编号'),ylabel('预测偏差')
title('GA-BP神经网络测试集的预测误差')
set(gca,'fontsize',12)
YReal=T_test(:,1);
YPred=bestT_sim(:,1);
n=size(T_test,1);
mae = mean(abs(YReal - YPred));
mape = mean(abs((YReal - YPred)./YReal));
rmse = sqrt(sum((YReal - YPred).^2)./n);
r2 = 1 - (sum((YPred - YReal).^2) / sum((YReal - mean(YReal)).^2));
display(['mae=',num2str(mae),',mape=',num2str(mape),',rmse=',num2str(rmse),',r2=',num2str(r2)]);
