tic % 計時器
%% 清空環境變數
close all
clear
clc
format compact
%% 資料提取
% 載入測試資料wine,其中包含的資料為classnumber = 3,wine:178*13的矩陣,wine_labes:178*1的列向量
load wine.mat
% 選定訓練集和測試集
% 將第一類的1-30,第二類的60-95,第三類的131-153做為訓練集
train_wine = [wine(1:30,:);wine(60:95,:);wine(131:153,:)];
% 相應的訓練集的標籤也要分離出來
train_wine_labels = [wine_labels(1
 
:30);wine_labels(60:95);wine_labels(131:153)];
% 將第一類的31-59,第二類的96-130,第三類的154-178做為測試集
test_wine = [wine(31:59,:);wine(96:130,:);wine(154:178,:)];
% 相應的測試集的標籤也要分離出來
test_wine_labels = [wine_labels(31:59);wine_labels(96:130);wine_labels(154:178)];
%% 資料預處理
% 資料預處理,將訓練集和測試集歸一化到[0,1]區間
[mtrain,ntrain] = size(train_wine);
[mtest,ntest] = size(test_wine);

dataset = [train_wine;test_wine];
% mapminmax為MATLAB
 
自帶的歸一化函式
[dataset_scale,ps] = mapminmax(dataset',0,1);
dataset_scale = dataset_scale';

train_wine = dataset_scale(1:mtrain,:);
test_wine = dataset_scale( (mtrain+1):(mtrain+mtest),: );
%% FA優化引數
% 引數向量 parameters [n N_iteration alpha betamin gamma]
% n為種群規模，N_iteration為迭代次數
para=[10,50,0.5,0.2,1];

% 待優化引數上下界 Simple
 
 bounds/limits for d-dimensional problems
d=2; % 待優化引數個數
Lb=[0.01,0.01]; % 下界
Ub=[100,100]; % 上界

% 引數初始化 Initial random guess
u0=Lb+(Ub-Lb).*rand(1,d);

% 迭代尋優 Display results
[bestsolutio,bestojb]=ffa_mincon_svm(@objfun_svm,u0,Lb,Ub,para,train_wine_labels,train_wine,test_wine_labels,test_wine);
%% 列印引數選擇結果
bestc=bestsolutio(1);
bestg=bestsolutio(2);

disp('列印選擇結果');
str=sprintf('Best c = %g，Best g = %g',bestc,bestg);
disp(str)
%% 利用最佳的引數進行SVM網路訓練
cmd_gwosvm = ['-c ',num2str(bestc),' -g ',num2str(bestg)];
model_gwosvm = svmtrain(train_wine_labels,train_wine,cmd_gwosvm);
%% SVM網路預測
[predict_label,accuracy] = svmpredict(test_wine_labels,test_wine,model_gwosvm);
% 列印測試集分類準確率
total = length(test_wine_labels);
right = sum(predict_label == test_wine_labels);
disp('列印測試集分類準確率');
str = sprintf( 'Accuracy = %g%% (%d/%d)',accuracy(1),right,total);
disp(str);
%% 結果分析
% 測試集的實際分類和預測分類圖
figure;
hold on;
plot(test_wine_labels,'o');
plot(predict_label,'r*');
xlabel('測試集樣本','FontSize',12);
ylabel('類別標籤','FontSize',12);
legend('實際測試集分類','預測測試集分類');
title('測試集的實際分類和預測分類圖','FontSize',12);
grid on
%% 顯示程式執行時間
toc

% 螢火蟲演算法主程式開始 Start FA
function [nbest,fbest]=ffa_mincon_svm(costfhandle,u0, Lb, Ub, para,train_wine_labels,train_wine,test_wine_labels,test_wine)
% 檢查輸入引數 Check input parameters (otherwise set as default values)
if nargin<5
    para=[20 100 0.25 0.20 1];
end
if nargin<4
    Ub=[];
end
if nargin<3
    Lb=[];
end
if nargin<2
    disp('Usuage: FA_mincon(@cost,u0,Lb,Ub,para)');
end

% n=number of fireflies
% MaxGeneration=number of pseudo time steps
% ------------------------------------------------
% alpha=0.25;      % Randomness 0--1 (highly random)
% betamn=0.20;     % minimum value of beta
% gamma=1;         % Absorption coefficient
% ------------------------------------------------
n=para(1);
MaxGeneration=para(2);
alpha=para(3);
betamin=para(4);
gamma=para(5);

% 檢查上界向量與下界向量長度是否相同 Check if the upper bound & lower bound are the same size
if length(Lb) ~=length(Ub)
    disp('Simple bounds/limits are improper!')
    return
end

% 計算待優化引數維度 Calcualte dimension
d=length(u0);

% 初始化目標函式值 Initial values of an array
zn=ones(n,1)*10^100;
% ------------------------------------------------
% 初始化螢火蟲位置 generating the initial locations of n fireflies
[ns,Lightn]=init_ffa(n,d,Lb,Ub,u0);

for k=1:MaxGeneration % 迭代開始
% 更新alpha（可選）This line of reducing alpha is optional
 alpha=alpha_new(alpha,MaxGeneration);

% 對每個螢火蟲計算目標函式值 Evaluate new solutions (for all n fireflies)
for i=1:n
    zn(i)=costfhandle(ns(i,:),train_wine_labels,train_wine,test_wine_labels,test_wine);
    Lightn(i)=zn(i);
end

% 根據亮度排序 Ranking fireflies by their light intensity/objectives
[Lightn,Index]=sort(zn);
ns_tmp=ns;
for i=1:n
    ns(i,:)=ns_tmp(Index(i),:);
end

%% 找出當前最優 Find the current best
nso=ns;
Lighto=Lightn;
nbest=ns(1,:);
Lightbest=Lightn(1);

% 另存最優值 For output only
fbest=Lightbest;

% 向較優方向移動 Move all fireflies to the better locations
[ns]=ffa_move(n,d,ns,Lightn,nso,Lighto,alpha,betamin,gamma,Lb,Ub);
end

% ----- All the subfunctions are listed here ------------
% 初始化螢火蟲位置 The initial locations of n fireflies
function [ns,Lightn]=init_ffa(n,d,Lb,Ub,u0)
ns=zeros(n,d);
if ~isempty(Lb) % 如果引數界限不為空 if there are bounds/limits
   for i=1:n
       ns(i,:)=Lb+(Ub-Lb).*rand(1,d); % 則在取值範圍內隨機取值
   end
else % 如果沒有設定引數界限
    for i=1:n
        ns(i,:)=u0+randn(1,d); % 在原有引數上加白噪聲
    end
end
% 初始化目標函式 initial value before function evaluations
Lightn=ones(n,1)*10^100;

% Move all fireflies toward brighter ones
function [ns]=ffa_move(n,d,ns,Lightn,nso,Lighto,alpha,betamin,gamma,Lb,Ub)
% 引數取值範圍絕對值 Scaling of the system
scale=abs(Ub-Lb);

% 更新螢火蟲 Updating fireflies
for i=1:n
    % The attractiveness parameter beta=exp(-gamma*r)
    for j=1:n
        r=sqrt(sum((ns(i,:)-ns(j,:)).^2));
        % Update moves
        if Lightn(i)>Lighto(j) % 如果i比j亮度更強 Brighter and more attractive
            beta0=1;
            beta=(beta0-betamin)*exp(-gamma*r.^2)+betamin;
            tmpf=alpha.*(rand(1,d)-0.5).*scale;
            ns(i,:)=ns(i,:).*(1-beta)+nso(j,:).*beta+tmpf;
        end
    end % end for j
end % end for i

% 防止越界 Check if the updated solutions/locations are within limits
[ns]=findlimits(n,ns,Lb,Ub);

% This function is optional, as it is not in the original FA
% The idea to reduce randomness is to increase the convergence,
% however, if you reduce randomness too quickly, then premature
% convergence can occur. So use with care.
% alpha引數更新函式 
function alpha=alpha_new(alpha,NGen)
% alpha_n=alpha_0(1-delta)^NGen=10^(-4);
% alpha_0=0.9
delta=1-(10^(-4)/0.9)^(1/NGen);
alpha=(1-delta)*alpha;

% 防止越界 Make sure the fireflies are within the bounds/limits
function [ns]=findlimits(n,ns,Lb,Ub)
for i=1:n
    % Apply the lower bound
    ns_tmp=ns(i,:);
    I=ns_tmp<Lb;
    ns_tmp(I)=Lb(I);
    % Apply the upper bounds
    J=ns_tmp>Ub;
    ns_tmp(J)=Ub(J);
    % Update this new move
    ns(i,:)=ns_tmp;
end

%% SVM_Objective Function
function f=objfun_svm(cv,train_wine_labels,train_wine,test_wine_labels,test_wine)
% cv為長度為2的橫向量，即SVM中引數c和v的值

cmd = [' -c ',num2str(cv(1)),' -g ',num2str(cv(2))];
model=svmtrain(train_wine_labels,train_wine,cmd); % SVM模型訓練
[~,fitness]=svmpredict(test_wine_labels,test_wine,model); % SVM模型預測及其精度
f=1-fitness(1)/100; % 以分類預測錯誤率作為優化的目標函式值

（廣告）歡迎掃描關注微信公眾號：Genlovhyy的資料小站(Gnelovy212)