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

'SaDE'

for problemIndex = 1:30
    
        problem = problemIndex
        optimum = problem*100;

        % Define the dimension of the problem
         n = 30;

         popsize = 50;

        % Record the best results
         outcome = []; 

        %Main body which was developed by the authors

          time = 1;
    
        % The total number of runs
          totalTime = 1;
    
          lu = [-100 * ones(1, n); 100 * ones(1, n)];
        % Initialize 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 population
          fit = cec14_func(p', problem);
          fit = fit';
         
          % The learning generation
          learn_l = 50;
         
          FES = popsize;
         
          G = 0;
         
          success_number = zeros(learn_l,4);
          all_number = zeros(learn_l,4);
                   
          % 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 = p(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);        
        
         % Initialize F and CR
          CR = zeros(popsize,1);
           F = zeros(popsize,1);
          CR_CML = zeros(popsize,1);
          F_CML = zeros(popsize,1);  
          
           % Initialzie the archive population and fitness value 
           popachive = p;
           fitachive = fit;
           
          % Initialize the selection probability vector sita
           sita = 0.5*ones(popsize,1);
                  
      while FES < n * 10000 %& min(fit)>error_value(problem)
             
             pTemp = p;
             fitTemp = fit;
             G = G + 1;
             
            %==============================================================
            % ============In the original coordinate system================
            %============================================================== 
             
            % Generate F
             F =  0.65 + 0.3 * randn(popsize,1);    
             F = min(1, F);                          % truncation

            % we don't want F = 0. So, if F<=0, we regenerate F (instead of trucating it to 0)
             pos = find(F <= 0 );
             while ~ isempty(pos)
                   F(pos) = 0.65+0.3*randn(length(pos), 1);
                   F = min(1, F);                      % truncation
                   pos = find(F <= 0 );
              end

            if G <= learn_l;
                trial_ratio =  0.25.*ones(1,4); 
                CRm = 0.5*ones(1,3);
            end
              
            [~ , indexBest] = min(fit);
            pbest = p(indexBest,:);
            
            % Record the trial vector generaiton strategy
            trial_record = zeros(popsize,4);
            trial_record_CML = zeros(popsize,4);
                 
         for i = 1:popsize
               % Select the trial vector strategy
               ratio_m = rand;
               index = floor(rand(1,5) * popsize) +1;
                 
             if ratio_m < trial_ratio(1)
                 trial_record(i,1) = 1;    % rand 1 bin
                 CR(i,:) = CRm(1) + 0.1*randn ;
                 while CR(i,:)<= 0 || CR(i)>1
                     CR(i,:) = CRm(1)+0.1*randn; 
                 end
                 v1 = p(index(1), :)  + F(i,:).*(p(index(2), :) - p(index(3), :));          
            elseif  trial_ratio(1)< ratio_m  &&  ratio_m <=sum(trial_ratio(1:2))
                 trial_record(i,2) = 1;  % rand 2 bin
                 CR(i,:) = CRm(2)+0.1*randn;
                 while CR(i,:)<= 0 || CR(i)>1
                     CR(i,:) = CRm(2) + 0.1*randn;
                 end
                v1 = p(index(1), :) +  F(i,:).* (p(index(2), :) - p(index(3), :)) + F(i,:).*(p(index(4), :) - p(index(5), :));    
            elseif  sum( trial_ratio(1:2)) < ratio_m  &&  ratio_m <= sum(trial_ratio(1:3)) 
                 trial_record(i,3) = 1;     % rand to best 2
                  CR(i,:) = CRm(3)+ 0.1*randn;
                  while CR(i,:)<= 0 || CR(i)>1
                     CR(i,:) =  CRm(3)+ 0.1*randn;
                 end
                v1 = p(i, :) + F(i,:)*(pbest-p(i,:))+ F(i,:).* (p(index(1), :) - p(index(2), :)) + F(i,:).*(p(index(3), :) - p(index(4), :)) ;   
             else   
                trial_record(i,4) = 1;  %current to rand
                CR(i,:) = 1;
                 v1 = p(i,:) +  rand.* (p(index(1), :) - p(i, :)) + F(i,:).*(p(index(2), :) - p(index(3), :));   
             end
            
              vi1(i,:) = v1;
         end   
          
          % == == == == =  Crossover == == == == =
            mask = rand(popsize, n) > CR(:, 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) = p(mask);
            
            % Boundary respair        
             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===================
            %==============================================================
             
            F_CML =  0.65 + 0.3 * randn(popsize,1);
            F_CML = min(1, F_CML);                          % truncation

            % we don't want F = 0. So, if F<=0, we regenerate F (instead of trucating it to 0)
            pos = find(F_CML <= 0 );
            while ~ isempty(pos)
                F_CML(pos) = 0.65 +0.3*randn(length(pos), 1);
                F_CML = min(1, F_CML);                      % truncation
                pos = find(F_CML <= 0 );
            end    
                
        for i = 1:popsize
            % To select the trial vector generation strategy
             ratio_m = rand; 
             index = floor(rand(1,5) * popsize) +1;
               
            if ratio_m < trial_ratio(1)
                 trial_record_CML(i,1) = 1;    % rand 1 bin
                 CR_CML(i,:) =  CRm(1)+ 0.1*randn;
                 while CR_CML(i,:)<=0 || CR_CML(i)>1
                     CR_CML(i,:) = CRm(1)+ 0.1*randn;
                 end
                 v1 = p(index(1), :)  + F_CML(i,:).*(p(index(2), :) - p(index(3), :));          
            elseif  trial_ratio(1)< ratio_m  &&  ratio_m <=sum(trial_ratio(1:2))
                 trial_record_CML(i,2) = 1;  % rand 2 bin
                 CR_CML(i,:) =  CRm(2)+ 0.1*randn;
                 while CR_CML(i,:)<=0 || CR_CML(i)>1
                     CR_CML(i,:) =  CRm(2)+ 0.1*randn;
                 end
                v1 = p(index(1), :) +  F_CML(i,:).* (p(index(2), :) - p(index(3), :)) + F_CML(i,:).*(p(index(4), :) - p(index(5), :));    
            elseif   sum(trial_ratio(1:2)) < ratio_m  &&  ratio_m <= sum(trial_ratio(1:3)) 
                 trial_record_CML(i,3) = 1;     % rand to best 2
                  CR_CML(i,:) =  CRm(3)+ 0.1*randn;
                  while CR_CML(i,:)<=0 || CR_CML(i)>1
                     CR_CML(i,:) = CRm(3)+ 0.1*randn;
                 end
                v1 = p(i, :) + F_CML(i,:)*(pbest - p(i,:))+ F_CML(i,:).* (p(index(1), :) - p(index(2), :)) + F_CML(i,:).*(p(index(3), :) - p(index(4), :)) ;   
            else   
                trial_record_CML(i,4) = 1;  %current to rand
                CR_CML(i,:) = 1;
                 v1 = p(i,:) +  rand.* (p(index(1), :) - p(i, :)) + F_CML(i,:).*(p(index(2), :) - p(index(3), :));   
            end
            
            vi2(i,:) = v1;
        end   
        
          % == == == == =  Crossover == == == == =
            J_= mod(floor(rand(popsize, 1)*n), n) + 1;
            J = (J_-1)*popsize + (1:popsize)';
            crs = rand(popsize, n) < CR_CML(:, ones(1, n));      
       
            vi2 = vi2 * B;
            ui2 = p * B;
            ui2(J) = vi2(J);
            ui2(crs) = vi2(crs);
            % returen
            ui2 = ui2 * B';
            
            
            % Boundary respair
             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
           number_count = zeros(1,popsize);
           for i = 1 : popsize
             if rand < sita(i)  
                 ui(i,:) = ui1(i,:);
                 CR_new(i,:) = CR(i,:);
                 trial_new(i,:)= trial_record(i,:);
                 number_count(i) = 1;
             else
                 ui(i,:) = ui2(i,:);
                 CR_new(i,:) = CR_CML(i,:);
                 trial_new(i,:)= trial_record_CML(i,:);
             end
           end
               
         % Evaluate the fitness
             fitSet = cec14_func(ui', problem);
             fitSet = fitSet';
                
             FES = FES + popsize;
          
         % Update the archive population
            popachive = [popachive;ui];
            fitachive = [fitachive;fitSet];
            
            T = 3;
            if size(popachive,1) >= T*popsize
                popachive(1:popsize,:) = [];
                fitachive(1:popsize,:) = [];
            end
            
            % Sort the individul in the archive population
            [~, index2] = sort(fitachive);
            mutant_sort = popachive(index2,:);
            p_select =  mutant_sort(1 : size(mutant_sort,1)./2,:);   
            
            % Compute the weights and mueff
            cov_length =size(mutant_sort,1)./2;
            weights = log(cov_length + 0.5) - log(1 : cov_length)';
            mueff = sum(weights)^2/sum(weights.^2);
            weights = weights/sum(weights);
            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;       
            C = (1 - cmu) * C + cmu .* C_u;
            C = triu(C) + triu(C, 1)';
            
           if mod(G,3) ==0;
            % Eigen decomposition
             [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
            p_center = (p_select' * weights)';
             
           
            count_achi = 0;
            trial_success = [];
            
            success_record = zeros(1,popsize);
              
            for i = 1 : popsize

                if   fitSet(i) <= fit(i)
                    
                    pTemp(i, :) = ui(i,:);
                    fitTemp(i, :) = fitSet(i,:);
                     count_achi = count_achi + 1;  % 加一则成功的进行了一次子代向父代的超越
                    trial_success(count_achi,:) = trial_new(i,:);
                    trial_success(count_achi,trial_new(i,:)~=0) = CR_new(i,:);   
                    
                    success_record(i) = 1;
                end
                      
             end   
             
            CRtemp = zeros(1,3);
            if count_achi ~= 0
              for i = 1: 3
                  CRtemp(1,i) =  mean( trial_success(:,i));
                if CRtemp(1,i) ~= 0
                  CRm(1,i) = CRtemp(1,i);
                end
              end
            end
            
            
            if count_achi ~= 0 
              success_number(1,:)= [];
              success_number(learn_l,:) = [  length(find(trial_success(:,1)~=0))  length(find(trial_success(:,2)~=0)) length(find(trial_success(:,3)~=0)) length(find(trial_success(:,4)~=0))];
            end
           
            if count_achi ~= 0 
            all_number(1,:) = [];  
            all_number(learn_l,:) =  [ length(find(trial_new(: , 1) ~= 0))  length(find(trial_new(: , 2) ~= 0)) length(find(trial_new(: , 3) ~= 0)) length(find(trial_new(: ,4) ~= 0)) ];
            end
             
            SKG = zeros(1,4);
            if G >= learn_l
                for i = 1 : 4
                    if sum(all_number(:,i)) == 0;
                     SKG(1,i) = 0.01;
                    else
                    SKG(1,i) = sum(success_number(:,i))/sum(all_number(:,i)) + 0.01;
                    end
                end
                for i = 1:4
                    trial_ratio(1,i) = SKG(1,i)/sum(SKG) ;
                end
            end
                        
           fit = fitTemp;
           p = pTemp;
           
           % Update the probability vector sita.
            para = 0.05;
             alfa = 2;
          for i = 1:popsize
                if number_count(i)==1 && success_record(i) == 1
                    sita(i) = sita(i) +para*(1-sita(i))*exp(-alfa*sita(i));
                elseif number_count(i)==0 && success_record(i) == 1 
                    sita(i) = sita(i) - para*(sita(i))*exp(alfa*sita(i)-alfa*1);
                elseif number_count(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;
