%%####################################################################
%%#### Particle swarm optimization
%%#### With linkage operator
%%#### Deepak devicharan july 2003
%%####################################################################

%%## to apply this to different equations do the following
%%## generate initial particles in a search space close to actual soln
%%## fool around with no of iterations, no of particles, learning rates

%%## for a truly generic PSO do the following
%%## increase the number of particles , increase the variance
%%## i.e let the particles cover a larger area of the search space
%%## then fool around as always with the above thins


%declare the parameters of the optimization


cancer_input = [-0.927 -0.919 -0.919 -0.968 -1.013 
    -0.919 -0.919 -0.968 -1.013 -1.035 
    -0.919 -0.968 -1.013 -1.035 -0.985 
    -0.968 -1.013 -1.035 -0.985 -1.003 
    -1.013 -1.035 -0.985 -1.003 -0.998 
    -1.035 -0.985 -1.003 -0.998 -1.007 
    -0.985 -1.003 -0.998 -1.007 -1.000 
    -1.003 -0.998 -1.007 -1.000 -0.992 
    -0.998 -1.007 -1.000 -0.992 -0.988 
    -1.007 -1.000 -0.992 -0.988 0.870 
    -1.000 -0.992 -0.988 -0.870 -0.874 
    -0.992 -0.988 -0.870 -0.874 -0.879 
    -0.988 -0.870 -0.874 -0.879 -0.873
    -0.870 -0.874 -0.879 -0.873 -0.929
    -0.874 -0.879 -0.873 -0.929 -0.922 
    -0.879 -0.873 -0.929 -0.922 -0.869 
    -0.873 -0.929 -0.922 -0.869 -0.958
    -0.929 -0.922 -0.869 -0.958 -0.899 
    -0.922 -0.869 -0.958 -0.899 -0.893 
    -0.869 -0.958 -0.899 -0.893 -0.944 
    -0.958 -0.899 -0.893 -0.944 -0.936
    -0.899 -0.893 -0.944 -0.936 -0.934 
    -0.893 -0.944 -0.936 -0.934 -0.887 
    -0.944 -0.936 -0.934 -0.887 -0.943 
    -0.936 -0.934 -0.887 -0.943 -0.940 
    ];
load cancer_output;

max_iterations  = 100;
no_of_particles = 20;
dimensions      = 2;

%delta_min       = -0.003;
%delta_max       = 0.003;

c1 = 1.3;
c2 = 1.3;

%initialise the particles and teir velocity components

for count_x = 1:no_of_particles
    for count_y = 1:dimensions
        particle_position(count_x,count_y) = rand*10;
        particle_velocity(count_x,count_y) = rand*1000;
        p_best(count_x,count_y) = particle_position(count_x,count_y);
    end
end

%initialize the p_best_fitness array
for count = 1:no_of_particles
    p_best_fitness(count) = -1000;
end


%particle_position
%particle_velocity

%main particle swrm routine
for count = 1:max_iterations
    
    %find the fitness of each particle
    %change fitness function as per equation requiresd and dimensions
    for count_x = 1:no_of_particles
        x = particle_position(count_x,1);
        y = particle_position(count_x,2);
%---------------------------------------------------
% 產生訓練樣本與測試樣本
% 特別注意：此工具箱用于分類時，只能處理2類分類，且目標值必須為 1 或 -1。

        %n1 = [rand(3,5),rand(3,5)+1];
        %x1 = [1*ones(1,5),-1*ones(1,5)];
        for i = 1:300
            for j = 1:9
                n1(j,i) = cancer_input(i,j);        
            end
        end
        for i = 1:300
            x1(i) = cancer_output(i)-3;
        end
    
        %n2 = [rand(3,5),rand(3,5)+1];
        %x2 = [1*ones(1,5),-1*ones(1,5)];
        for i = 1:67
            for j = 1:9
                n2(j,i) = cancer_input(300+i,j);        
            end
        end
        for i = 1:67
            x2(i) = cancer_output(300+i)-3;
        end

        xn_train = n1;          % 訓練樣本，每一列為一個樣本
        dn_train = x1;          % 訓練目標，行向量

        xn_test = n2;           % 測試樣本，每一列為一個樣本
        dn_test = x2;           % 測試目標，行向量

%---------------------------------------------------
% 參數設置

        trnX = xn_train';
        trnY = dn_train';
        tstX = xn_test';
        tstY = dn_test';

        ker = 'rbf';        % 核函數 k = exp(-(u-v)*(u-v)'/(2*p1^2))
        global p1 ;
        p1 = x;             % p1 is width of rbfs (sigma)
        C = y;             % 折衷系數

%---------------------------------------------------
% 訓練與測試

        [nsv,alpha,bias] = svc(trnX,trnY,ker,C);                        % 訓練

        actfunc = 0;                                                    % 1 為實際輸出，0 為取sign輸出 
        predictedY = svcoutput(trnX,trnY,tstX,ker,alpha,bias,actfunc);  % 測試

%---------------------------------------------------
% 結果統計

        Result = ~abs(predictedY-tstY)               % 正確分類顯示為1
        Percent = sum(Result)/length(Result)   % 正確分類率
        soln = 1-Percent
        
        %x = particle_position(count_x,1);
        %y = particle_position(count_x,2);
        %z = particle_position(count_x,3);
        %soln = x^2 - 3*y*x + z;
        
        %x = particle_position(count_x); 
        %soln = x^2-2*x+1;
        
       % x = particle_position(count_x);
   % soln = x-7;
        
        if soln~=0 
            current_fitness(count_x) = 1/abs(soln)+0.0001;
        else
            current_fitness(count_x) =1000;
        end
    end
    
    %decide on p_best etc for each particle
    for count_x = 1:no_of_particles
        if current_fitness(count_x) > p_best_fitness(count_x)
            p_best_fitness(count_x) = current_fitness(count_x);
            for count_y = 1:dimensions
                p_best(count_x,count_y) = particle_position(count_x,count_y);
            end
        end
    end
        
    %decide on the global best among all the particles
    [g_best_val,g_best_index] = max(current_fitness);
    
    %g_best contains the position of teh global best
    for count_y = 1:dimensions
        g_best(count_y) = particle_position(g_best_index,count_y);      
    end

    %update the position and velocity compponents
    for count_x = 1:no_of_particles
        for count_y = 1:dimensions
            p_current(count_y) = particle_position(count_x,count_y);
        end

        for count_y = 1:dimensions
            particle_velocity(count_y) = particle_velocity(count_y) +  c1*rand*(p_best(count_y)-p_current(count_y)) + c2*rand*(g_best(count_y)-p_current(count_y));
            particle_positon(count_x,count_y) = p_current(count_y) +particle_velocity(count_y);
        end
    end
                                 
              
end

g_best        
current_fitness(g_best_index)