%**************************************************************************************************
%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.
%**************************************************************************************************

%**************************************************************************************************
warning off;
clc;
clear all;
tic;



format long;
format compact;

'ACoS-JADE'

% 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)
       
    %Set the golobal optimum
    optimum = problem*100;
    
    % Define the dimension of the problem
    n = 30;
    
    % Set the population size
    popsize = 100;
    
    % Define the boundary
    lu = [-100 * ones(1, n); 100 * ones(1, n)];

    % Record the best results
    outcome = []; 

    %Main body which was developed by the authors
    time = 1;
    
    % The total number of runs
    totalTime = 1;
    
    while time <= totalTime
   
        rand('seed', sum(100 * clock));

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

        
        valParents =cec14_func(popold', problem);
        valParents = valParents';
        
        c = 1/10;
        p = 0.05;

        CRm = 0.5;
        Fm = 0.5;
        
        % Initialize the selection ration in the Eigen coordiante system as sita = (0.5, 0.5, ...0.5).
        sita = 0.5.*ones(1,popsize);

        % Record the values and indices of the best solutions
        [~, indBest] = sort(valParents, 'ascend');
         
        % Compute the weights and mueff
        cov_length = popsize./2;
        weights = log(cov_length + 0.5) - log(1 : cov_length)'; 
        mueff = sum(weights)^2/sum(weights.^2);
        weights = weights/sum(weights);
        
        % Initialize the weighted mean value of the distribution: p_center
        p_sort = popold(indBest(1:popsize./2,:),:);
        p_center =(p_sort' * weights)';
        
        %Initialzie the archive population and fitness value in SVM-EA
        pop_pool = popold;
        fit_pool = valParents;
        center_pool = p_center;
     
        % Initialize covariance matrix C and B as unity matrix.
        B = eye(n);
        D = eye(n);
        BD = B * D;
        C = BD * transpose(BD);        
            
        FES = popsize;
        
        % The arvhive population in JADE itself
        popachieve = [];
        fitachieve = [];
        G = 0;
        
        while FES <  10000*n %& min(fit)>error_value(problem)
             
             % The old population becomes the current population
             pop = popold; 
             valTemp  = valParents ;

            % Generate CR according to a normal distribution with mean CRm, and std 0.1
            % Generate F according to a cauchy distribution with location parameter Fm, and scale parameter 0.1
            [F_original, CR_original] = randFCR(popsize, CRm, 0.1, Fm, 0.1);
            [F_Eigen, CR_Eigen] = randFCR(popsize, CRm, 0.1, Fm, 0.1);
            
            % The general population
             popAll = [pop; popachieve];
             
            %==============================================================
            % ============In the original coordinate system================
            %============================================================== 
            
            r0 = [1 : popsize];
            [r1, r2] = gnR1R2(popsize, size(popAll, 1), r0);
            [~, indBest] = sort(valParents, 'ascend');
         
            % Find the p-best solutions
            pNP = max(round(p * popsize), 2); % choose at least two best solutions
            randindex = ceil(rand(1, popsize) * pNP); % select from [1, 2, 3, ..., pNP]
            randindex = max(1, randindex); % to avoid the problem that rand = 0 and thus ceil(rand) = 0
            pbest = pop(indBest(randindex), :); % randomly choose one of the top 100p% solutions

            %  = == == == == == == == ==Mutation == == == == == == == == == 
            vi1 = pop + F_original(:, ones(1, n)) .* (pbest - pop + pop(r1, :) - popAll(r2, :));
            
            % == == == == == == == == ==Crossover == == == == == == == == =
            mask = rand(popsize, n) > CR_original(:, ones(1, n));     % mask is used to indicate which elements of ui comes from the parent
            rows = (1:popsize)'; cols = floor(rand(popsize, 1) * n) + 1;  % choose one position where the element of ui doesn't come from the parent
            jrand = sub2ind([popsize n], rows, cols); mask(jrand) = false;  
            ui1 = vi1;  ui1(mask) = pop(mask);
            
            % == == == == == == == Boundary repair == == == == == == == ==
            ui1 = boundConstraint(ui1, pop, lu);     
            
       
            
            %==============================================================
            % ============In the Eigen coordinate system===================
            %==============================================================

            r0 = [1 : popsize];
            [r1, r2] = gnR1R2(popsize, size(popAll, 1), r0);
            [~, indBest] = sort(valParents, 'ascend');
            
           % Find the p-best solutions
            pNP = max(round(p * popsize), 2); % choose at least two best solutions
            randindex = ceil(rand(1, popsize) * pNP); % select from [1, 2, 3, ..., pNP]
            randindex = max(1, randindex); % to avoid the problem that rand = 0 and thus ceil(rand) = 0
            pbest = pop(indBest(randindex), :); % randomly choose one of the top 100p% solutions

           % == == == == == == == == == Mutation == == == == == == == == == 
            vi2 = pop + F_Eigen(:, ones(1, n)) .* (pbest - pop + pop(r1, :) - popAll(r2, :));
            
           % == == == === == == == == = Crossover == == == == == == == == = 
            J_= mod(floor(rand(popsize, 1)*n), n) + 1;
            J = (J_-1)*popsize + (1:popsize)';
            crs = rand(popsize, n) < CR_Eigen(:, ones(1, n));            
                      
            vi2 = vi2 * B;   ui2 = pop * B; % turn into the Eigen coordiante system
            ui2(J) = vi2(J); ui2(crs) = vi2(crs);
            ui2 = ui2 * inv(B);  % return back
            
            % == == == == == == == Boundary repair == == == == == == == ==
            ui2 = boundConstraint(ui2, pop, lu);
            
           % Select the coordiante system based on the selection vecotr sita to implement the offspring 
           system_selected = zeros(1,popsize);  
            for i = 1 : popsize 
                if rand<sita(i)
                   ui(i,:) = ui1(i,:);
                   F_new(i,:) = F_original(i,:);
                   CR_new(i,:) = CR_original(i,:);
                   system_selected(i) = 1;
                else
                   ui(i,:) = ui2(i,:);
                    F_new(i , :) = F_Eigen(i,:);
                  CR_new( i,:) = CR_Eigen(i,:);   
                 end
            end
            
           % Evaluation the trial population  
            valOffspring = cec14_func(ui', problem);
            valOffspring = valOffspring';

            FES = FES + popsize;

            % Set the archive popsize as T*popsize   
            T  = 3; 
            
            % Update the archieve population based on the rule of "first-in and first-out"
            pop_pool = [pop_pool;ui];
            fit_pool = [fit_pool;valOffspring];
          
            if size(pop_pool,1) > T*popsize
                pop_pool(1:popsize,:) = [];
                fit_pool(1:popsize,:) = [];
            end
            
            G = G + 1;
         
            % Sort the individul in the archive population in ACoS
            [~, 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,:);   
             
            % 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);
        
            % Compute the learning ratio cmu
             cmu =min(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)';
            
            % 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)';
            
           % == == ==   == == ==  == ==Selection  == == ==  == == == == ==  
           F_success = [];
           CR_success = [];
  
           pcount_achi = [];
           fitcount_achi = [];
            
           sita_success= [];
           value_differ = [];
           count_achi = 0;
           number_success = [];
            
           success_record = zeros(1,popsize);
           
           for i = 1:popsize
                if valTemp(i) >= valOffspring(i)
                    pop(i,:) = ui(i,:);
                    valParents(i,:) = valOffspring(i,:);  
                    
                    count_achi = count_achi + 1;
                    pcount_achi(count_achi,:) = popold(i,:);
                    fitcount_achi(count_achi,:) = valTemp(i,:);
                    
                    F_success(count_achi) = F_new(i);
                    CR_success(count_achi) = CR_new(i);
                    
                    number_success(count_achi) = system_selected(i);
                    
                    success_record(i) = 1;
                end 
           end
                   

           popold = pop;
           valTemp = valParents;
            
            
          % Update Fm and CRm
           if count_achi~= 0;
              goodF = F_success;
              goodCR = CR_success;
            
              CRm = (1 - c) * CRm + c * mean(goodCR);
              Fm = (1 - c) * Fm + c * sum(goodF .^ 2) / sum(goodF);
           end
            
           % Upadate  the 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);

            
            %  Archive mechanism in JADE itself
            [popachieve, fitachieve] = update_achive(popachieve, fitachieve, pcount_achi, fitcount_achi, popsize);
            
        end

          outcome = [outcome min(valParents-optimum)];

           time = time + 1;

     end

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

end
toc;
