%**************************************************************************************************
%Authors: Zhi-Zhong Liu and Yong Wang
%Last Edited: Janaury 17, 2017
%Email: zhizhongliu@csu.edu.cn and ywang@csu.edu.cn
%Reference:   Z.-Z. Liu, Y. Wang, S. Yang, and K. Tang. An adaptive framework to tune the coordinate systems in nature-inspired optimization algorithms. 
%              IEEE Transactions on Cybernetics, in press.
%**************************************************************************************************

clc;
clear all;
tic;

format long;
format compact;

'ACoS-PSO'

% Choose the problems to be tested. Please note that for test functions F7
% and F25, the global optima are out of the initialization range. For these
% two test functions, we do not need to judge whether the variable violates
% the boundaries during the evolution after the initialization.
problemSet = [1 : 30];
for problemIndex = 1:30

    problem = problemSet(problemIndex)
    
    % record the optimum of the problem
    optimum = problem.*100;

    % Define the dimension of the problem
    n = 30;
    
    %Define the boudary
    lu = [-100 * ones(1, n); 100 * ones(1, n)];

    % Record the best results
    outcome = [];
    outcome_fitness=[];
    
    % the setting of the population size.
    if n==30
        popsize = 40;  % n = 30;
    else
        popsize = 60;  % n = 50;
    end
 
   
    % Main body
    time = 1;

    % The total number of runs
    totalTime = 1;
    
    % The boundary
    LB = repmat(lu(1,:), popsize, 1);
    UB = repmat(lu(2,:), popsize, 1);

    while time <= totalTime

        rand('seed', sum(100*clock));

        % Initialize the main population
        p = repmat(lu(1, :), popsize, 1) + rand(popsize, n) .* (repmat(lu(2, :) - lu(1, :), popsize, 1));

        % Evaluate the objective function values
        fit =  cec14_func(p', problem);
        fit = fit';
        
        % Initialize the archive population and fitness value.
        pop_pool = p;
        fit_pool = fit;
       
        FES_max = n * 10000;
   
        % Record history best 
        p_h_best = p;
        fit_h_best = fit;
        
        % Record the number of function evaluations (FES)
        FES = popsize;
        
        % Initialize the veloticy of the particle
        v = repmat(lu(1, :),  popsize, 1) + rand(popsize, n) .* (repmat(lu(2, :) - lu(1, :), popsize, 1));
        
        % One half best indididuals is selected to compute p_center
        cov_length = popsize/2;
 
        % Compute the weights and mueff.
        weights = log(cov_length + 0.5) - log(1 : cov_length)';
        mueff = sum(weights)^2/sum(weights.^2);
        weights = weights/sum(weights);
        
        % The initialization of p_center
        [val, index] = sort(fit);
        p_sort = p(index(1:cov_length),:);
        p_center =(p_sort'* weights)';
     
        % Initial Covarriance C and B as a unity matrix.
        B = eye(n);
        D_c = eye(n);
        BD = B * D_c;
        C = BD * transpose(BD);        
        
        % Initialize the selection vector sita = (0.5, 0.5, ...0.5).
        sita = 0.5*ones(popsize,1);
        
        % Initialize the geneartion number.
        G = 0;
        
        while FES < FES_max 
                
             % *******************      PSO-cf    **************************      
              c1=1.4962;             
              c2=1.4962;  
              w=0.7298;            

%           % *******************      PSO-w     **************************
%                w = 0.9  - 0.5 * FES./FES_max;  
%                c1 =  2.0;
%                c2 =   2.0 ;         

            % Find the global history best individual
              [ ~, index] = min(fit_h_best);
              global_best = p_h_best(index,:);
                
            %velocity renew
            % in the original  coordiante system
             v1 = w.*v + c1.*(p_h_best - p).*rand(popsize,n) + c2.*(global_best(ones(1,popsize),:) - p).*rand(popsize,n);
            % in the Eigen coordinate system
             v2 = w.*v + c1.*(p_h_best - p)*B.*rand(popsize,n)*B' + c2.*(global_best(ones(1,popsize),:) - p)*B.*rand(popsize,n)*B';
            
            % Selecte the coordinate system based on the probability vector sita and record the information
             system_selected = zeros(popsize,1);
             for i = 1:popsize
                if rand< sita(i)
                  v(i,:) = v1(i,:);
                  system_selected(i) = 1;
                else 
                  v(i,:) = v2(i,:);
                end
             end
             
             % respair the velocity
             BL  = v <LB; v(BL) =LB(BL);   
             BU  = v>UB; v(BU) = UB(BU);
           
             % Obtain the new particle swarms
             p_new = p+ v;
                 
             % Respair the particles  
             BL  = p_new <LB;  p_new(BL) = LB(BL);           
             BU  = p_new>UB; p_new(BU) = UB(BU);
              
             FES = FES + popsize;
            
             % Evaluate the new paricles
             fitSet = cec14_func(p_new', problem);
             fitSet = fitSet';
            
            % Set the archive popsize as T*popsize
             T = 3;
            
            % Update the archive based on the  "first-in and first-out" rule
            pop_pool = [pop_pool;p_new];
            fit_pool = [fit_pool;fitSet];
             
            if size(pop_pool,1) > T*popsize
                pop_pool(1:popsize,:) = [];
                fit_pool(1:popsize,:) = [];
            end
            
            %Sort the individul in the Archive population 
            [~, index2] = sort(fit_pool);
            mutant_sort = pop_pool(index2,:);   
            
            % one half best individual is selected 
            cov_length= size(mutant_sort,1)./2;
            p_select =  mutant_sort(1 : cov_length,:);   
            
            % Comput the weight and mueff
            weights = log(cov_length + 0.5) - log(1 : cov_length)';
            mueff = sum(weights)^2/sum(weights.^2);
            weights = weights/sum(weights);
           
            % Compute the learning ratio cmu
             cmu = min(1/3.*mueff./n^2,1);
  
            %  Employ the Rank-mu-update method to estimate C_u
            p_to_mean =  p_select - p_center(ones(1,cov_length),:);       
            C_u =p_to_mean'*diag(weights)*p_to_mean;
            
           % Utilize the cumulative information to estimate C
           C = (1 - cmu ) * C + cmu .* C_u ;    
           C = triu(C) + triu(C, 1)';
            
           G = G+1;
           
           % Eigen decomposition of C
           if mod(G,3)==0
               [B, D] = eig(C);
           % limit condition of C to 1e20 + 1
              if max(diag(D)) > 1e20*min(diag(D))
                tmp = max(diag(D))/1e20 - min(diag(D));
                C = C + tmp*eye(n);
                [B, D] = eig(C);
             end
           end
            
           % Compute the weighted center of population distribution
            p_center = (p_select' * weights)';
          
            % Record the performance of the offspring in the coordinate system
           success_record = zeros(popsize,1);
           for i = 1 : popsize
                if fitSet(i,:) < fit_h_best(i,:);
                    p_h_best(i,:) = p_new(i,:);
                    fit_h_best(i,:) = fitSet(i,:);
                    success_record(i) = 1;
                end
              
           end
           
           p = p_new;
           fit = fitSet;  
           
           % Upadate the selection probability vector: sita
           para = 0.05;
           alfa = 2;
           
           for i = 1:popsize
                if system_selected(i)==1 && success_record(i) == 1
                    sita(i) = sita(i) +para*(1-sita(i))*exp(-alfa*sita(i));
                elseif system_selected(i)==0 && success_record(i) == 1 
                    sita(i) = sita(i) - para*(sita(i))*exp(alfa*sita(i)-alfa*1);
                elseif system_selected(i)==1 &&success_record(i) ==0
                    sita(i) = sita(i) -0.1* para*(1-sita(i))*exp(-alfa*sita(i));
                else 
                    sita(i) = sita(i) + 0.1*para*(sita(i))*exp(alfa*sita(i)-alfa*1);
                end
           end
           
            sita = max(min(sita,1),0);
         
        end

        outcome = [outcome min(fit_h_best-optimum)];
        time = time + 1;

    end

    sort(outcome)
    mean(outcome)
    std(outcome)

end
toc;
    
