%**************************************************************************************************
%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;
clear all
clc;
tic;
format long;
format compact;
'SVM_jDE'

problemSet = [1:30];
 
for problemIndex= 1:30   
    
    problem = problemSet(problemIndex)
   
    %Set the golobal optimum
    optimum = problem * 100.0;
    
    %Define the dimension of the problem
    n = 30;
    
    % Set the population size
    popsize = 100;

    tau1 = 0.1;
    tau2 = 0.1;

    % Define the boundary
    lu = [-100 * ones(1, n); 100 * ones(1, n)];
    
    % Record the best results
    outcome = []; 

    %Main body which was provided by Dr. J. Zhang
    time = 1; 
    
    % The total number of runs
    totalTime = 1;
    
    LB = repmat(lu(1,:), popsize, 1);
    UB = repmat(lu(2,:), popsize, 1);
    
    while time <= totalTime
        
        %Initialize the setting of F and CR.
         F = 0.5 * ones(popsize, 1);
        CR = 0.9 * ones(popsize, 1);
         
        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));

        fit = cec14_func(popold',problem);
        fit = fit';
        
        
        % Initialzie the archive population and fitness value 
        pop_pool = popold;
        fit_pool = fit;
        
        
        % Compute the weights and mueff
        cov_length = popsize;  
        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
        [val, index] = sort(fit);
        p_sort = popold(index,:);
        p_center =(p_sort' * weights)';
     
        % Initialize the covariance matrix C and B as the unity matrix
        B = eye(n);
        D = eye(n);
        BD = B * D;
        C = BD * transpose(BD);        
        
        FES = popsize;
       
        % Set the maxium fitness evalution numbers
        max_FES = n*10000;
        
        % Initialize the probability vector sita = (0.5, 0.5, ...0.5).
        sita = 0.5*ones(popsize, 1);
        
        % Initialize the generation number.
        G = 0;
        
        while FES < n * 10000
            
            % The old population becomes the current population
            pop = popold;      
            fitold = fit;
           
            Fold = F;
            CRold = CR;
            
            F_original = F;
            CR_original = CR;
            F_Eigen = F;
            CR_Eigen = CR;
            
            IF  = rand(popsize, 1) < tau1;
            ICR = rand(popsize, 1) < tau2;

            F_original(IF)  = 0.1 + 0.9 * rand(sum(IF), 1);
            CR_original(ICR) = 0.0 + 1.0 * rand(sum(ICR), 1);

            
            %==============================================================
            % ============In the original coordinate system================
            %============================================================== 
            
            r0 = [1:popsize];
            [r1, r2, r3] = gnR1R2R3(popsize, r0);

            %  == == == == == == == == = Mutation == == == == == == == == =
            vi1  = pop(r1, :) + F_original(:, ones(1, n)).* (pop(r2, :) - pop(r3, :));

          

            % == == == == == == == == == 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 == == == == == == == ==
            BL  = ui1<LB; ui1(BL) = 2*LB(BL) - ui1(BL);
            BLU = ui1(BL)>UB(BL); BL(BL) = BLU; ui1(BL) = UB(BL);
            BU  = ui1>UB; ui1(BU) = 2*UB(BU) - ui1(BU);
            BUL = ui1(BU)<LB(BU); BU(BU) = BUL; ui1(BU) = LB(BU);

    
            
            %==============================================================
            % ============In the Eigen coordinate system===================
            %==============================================================
            IF  = rand(popsize, 1) < tau1;
            ICR = rand(popsize, 1) < tau2;

            F_Eigen(IF)  = 0.1 + 0.9 * rand(sum(IF), 1);
            CR_Eigen(ICR) = 0.0 + 1.0 * rand(sum(ICR), 1);

            r0 = [1:popsize];
            [r1, r2, r3] = gnR1R2R3(popsize, r0);

            %  == == == == == == == == Mutation == == == == == == == == ==
            vi2  = pop(r1, :) + F_Eigen(:, ones(1, n)) .* (pop(r2, :) - pop(r3, :));

            %== == == === == == == == = 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 coordinate system
            ui2(J) = vi2(J); ui2(crs) = vi2(crs);
            ui2 = ui2 * B'; % return back
            
            % == == == == == == == Boundary repair == == == == == == == ==
            BL  = ui2<LB; ui2(BL) = 2*LB(BL) - ui2(BL);
            BLU = ui2(BL)>UB(BL); BL(BL) = BLU; ui2(BL) = UB(BL);
            BU  = ui2>UB; ui2(BU) = 2*UB(BU) - ui2(BU);
            BUL = ui2(BU)<LB(BU); BU(BU) = BUL; ui2(BU) = LB(BU);
      
           % Select the coordiante system to yeild the offspring based
           % on the probability vector sita
           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

           fitSet = cec14_func(ui',problem);
           fitSet = fitSet';
           FES = FES + popsize ;
         
           % Set the archive popsize as T*popsize   
            T = 3;
            pop_pool = [pop_pool;ui];
            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,:);   
            
            % 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 = 1/3*mueff./n^2;
 
            % 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 == == == == == == == == =
             success_record = zeros(1,popsize);
             for i = 1 : popsize
                if fitold(i) >= fitSet(i)
                    
                    pop(i, :) =ui(i,:);
                    fit(i, :) = fitSet(i,:);
                    Fold(i,:) = F_new( i,:); 
                    CRold(i, :) = CR_new(i,:);
                    success_record(i) = 1;
                end
             end
             F = Fold;
             CR = CRold;
             popold = pop;
             fitold = fit;
            
             % Updating of 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);
       
        end

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

    end

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

end
toc;





