close all
clear
clc
fuhe = xlsread('XSR2.xlsx');
%归一化
fuhe = fuhe';
Y = mapminmax(fuhe,0,1);
fuhe = Y';
for u = 5:100
train_x = [fuhe(1:88,1:4)];
train_y = [fuhe(1:88,u)];
test_x = [fuhe(89,1:4)];
test_y = [fuhe(89,u)];
%参数初始化
NP=20; % 蜂群规模
FoodNumber=NP/2; % 蜜源（解）数量
limit=100; % 当有蜜源连续没被更新的次数超过limit时，该蜜源将被重新初始化
maxCycle=200; % 最大迭代次数
% 待优化参数信息
D=2; % 待优化参数个数，次数为c和g两个
ub=ones(1,D)*50; % 参数取值上界
lb=ones(1,D)*(0.01); % 参数取值下界
runtime=2; % 可用于设置多次运行（让ABC算法运行runtime次）以考察程序的稳健性

BestGlobalMins=ones(1,runtime); % 全局最小值初始化，这里的优化目标为SVR预测结果中的平均平方误差（MSE），初始化为最差值1
BestGlobalParams=zeros(runtime,D); % 用于存放ABC算法优化得到的最优参数
for r=1:runtime % 运行ABC算法runtime次
    % 初始化蜜源
    Range = repmat((ub-lb),[FoodNumber 1]);
    Lower = repmat(lb, [FoodNumber 1]);
    Foods = rand(FoodNumber,D) .* Range + Lower;
 % 计算每个蜜源（解）得目标函数值，fobj为计算SVR预测的平均平方误差（MSE）的函数
    ObjVal=ones(1,FoodNumber); 
 for k = 1:FoodNumber
        ObjVal(k) = svfobj(Foods(k,:),train_x,train_y,test_x,test_y);
 end
 Fitness=calculateFitness(ObjVal); % 计算适应度函数值
    
    trial=zeros(1,FoodNumber); % 用于记录第i个蜜源有连续trail(i)次没被更新过

    % 标记最优蜜源（解）
    BestInd=find(ObjVal==min(ObjVal)); 
    BestInd=BestInd(end);
    GlobalMin=ObjVal(BestInd); % 更新全局最优目标函数值
    GlobalParams=Foods(BestInd,:); % 更新全局最优参数为最优蜜源
     iter=1; % 迭代开始
    while ((iter <= maxCycle)) % 循环条件

%%%%%%%%%%%%%%%%%%%%%引领蜂搜索解的过程%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        for i=1:(FoodNumber) % 遍历每个蜜源（解）
            Param2Change=fix(rand*D)+1; 
            neighbour=fix(rand*(FoodNumber))+1; 
            % 需要保证选择的相邻蜜源不是当前蜜源（i）
            while(neighbour==i)
                neighbour=fix(rand*(FoodNumber))+1;
            end
            sol=Foods(i,:); % 提取当前蜜源（解）对应的的参数
            % 新的蜜源：v_{ij}=x_{ij}+\phi_{ij}*(x_{kj}-x_{ij})
            sol(Param2Change)=Foods(i,Param2Change)+(Foods(i,Param2Change)-Foods(neighbour,Param2Change))*(rand-0.5)*2;
            % 确保参数取值范围不越界
            ind=find(sol<lb);
            sol(ind)=lb(ind);
            ind=find(sol>ub);
            sol(ind)=ub(ind);
            % 计算新蜜源的目标函数值和适应度函数值
            ObjValSol=svfobj(sol,train_x,train_y,test_x,test_y);
            FitnessSol=calculateFitness(ObjValSol);
            % 更新当前蜜源的相关信息
            if (FitnessSol>Fitness(i))  %适应度值越大越好，适应度值越大 mse越小
                Foods(i,:)=sol;
                Fitness(i)=FitnessSol;
                ObjVal(i)=ObjValSol;
                trial(i)=0; % 如果当前蜜源被更新了，则对应的trial归零
            else
                trial(i)=trial(i)+1; % 如果当前蜜源没有被更新，则trial(i)加1
            end
        end
        % 跟随蜂搜索解的过程 %
        % 计算解（蜜源）的概率
        prob=(0.9.*Fitness./max(Fitness))+0.1;
        % 循环初始化
        i=1;
        t=0;
        while(t<FoodNumber) % 循环条件
            if(rand<prob(i)) % 若随机概率小于当前解（蜜源）的概率
                t=t+1; % 循环计数器加1
                
                Param2Change=fix(rand*D)+1; 
                neighbour=fix(rand*(FoodNumber))+1; 
                % 需要保证选择的相邻蜜源不是当前蜜源（i）
                while(neighbour==i)
                    neighbour=fix(rand*(FoodNumber))+1;
                end
                sol=Foods(i,:); % 提取当前蜜源i（解）对应的的参数
                % 新的蜜源：v_{ij}=x_{ij}+\phi_{ij}*(x_{kj}-x_{ij})
                sol(Param2Change)=Foods(i,Param2Change)+(Foods(i,Param2Change)-Foods(neighbour,Param2Change))*(rand-0.5)*2;
                % 防止参数越界
                ind=find(sol<lb);
                sol(ind)=lb(ind);
                ind=find(sol>ub);
                sol(ind)=ub(ind);
                % 计算新蜜源的目标函数值和适应度函数值
                ObjValSol=svfobj(sol,train_x,train_y,test_x,test_y);
                FitnessSol=calculateFitness(ObjValSol);
                % 更新当前蜜源的相关信息
                if (FitnessSol>Fitness(i))
                    Foods(i,:)=sol;
                    Fitness(i)=FitnessSol;
                    ObjVal(i)=ObjValSol;
                    trial(i)=0; % 如果当前蜜源被更新了，则对应的trial归零。
                else
                    trial(i)=trial(i)+1; % 如果当前蜜源没有被更新，则trial(i)加1。在统计每个蜜源未被更新的次数，为下个步骤的侦查蜂做准备。
                end
            end
    
            i=i+1; % 更新i
            if (i==(FoodNumber)+1) % 若值超过蜜源数量，则i重新初始化
                i=1;
            end   
        end 
        % 记住最优蜜源
         ind=find(ObjVal==min(ObjVal));
         ind=ind(end);
         if (ObjVal(ind)<GlobalMin)
            GlobalMin=ObjVal(ind);
            GlobalParams=Foods(ind,:);
         end
        
% 侦查蜂搜索解的过程 %
        % 找出连续最多次都没有被更新的蜜源        
        ind=find(trial==max(trial)); 
        ind=ind(end);
        % 如果连续没有更新的次数大于限定次数，则由侦查蜂重新初始化该蜜源
        if (trial(ind)>limit) 
            Bas(ind)=0;
            sol=(ub-lb).*rand(1,D)+lb;
            ObjValSol=svfobj(sol,train_x,train_y,test_x,test_y);
            FitnessSol=calculateFitness(ObjValSol);
            Foods(ind,:)=sol;
            Fitness(ind)=FitnessSol;
            ObjVal(ind)=ObjValSol;
        end
        
        iter=iter+1;

    end % 一次ABC算法完结
     BestGlobalMins(r)=GlobalMin; % 记录本次ABC算法的最优目标函数值
    BestGlobalParams(r,:)=GlobalParams; % 记录本次ABC算法的最优参数
end % end of runs

bestc=GlobalParams(1);
Data1(u,:) = bestc
bestg=GlobalParams(2);
Data2(u,:) = bestg
disp('打印选择结果');
str=sprintf('Best c = %g，Best g = %g',bestc,bestg);
disp(str)
%% 利用回归预测分析最佳的参数进行SVM网络训练
cmd_cs_svr=['-s 3 -t 2',' -c ',num2str(bestc),' -g ',num2str(bestg)];
model_cs_svr=libsvmtrain(train_y,train_x,cmd_cs_svr); % SVM模型训练
%% SVM网络回归预测
[output_test_pre,acc,~]=libsvmpredict(test_y,test_x,model_cs_svr); % SVM模型预测及其精度
test_pre = output_test_pre;
Data(u,:) = test_pre
end