%**************************************************************************************************
%Reference:  J. Zhang and A. C. Sanderson, "JADE: adaptive differential evolution
%                     with optional external archive," IEEE Trans. Evolut. Comput., vol. 13,
%                     no. 5, pp. 945-958, 2009.
%
% Note: We obtained the MATLAB source code from the authors, and did some
%           minor revisions in order to solve the 25 benchmark test functions,
%           however, the main body was not changed.
%*************************************************************************************************
clc;
clear all;
% mex cec13_func.cpp -DWINDOWS
tic;

format long;
format compact;

'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.
for   problem = 14
    
    varargin{:}=problem;
    
    % Define the dimension of the problem
    n = 30;
    
    fhd=str2func('cec13_func');
    
    popsize = 50;
    
    lu = [-100 * ones(1, n); 100 * ones(1, n)];
     
    Strategy=3;
    switch Strategy
        case 0
            circle=5000;
        case 1
            circle=2971; % Maximum numbef of iterations 4839 & 3572
        case 2
            circle=2679;
        case 3
            circle=5883;
    end
    % Record the best results
    FbestChart= [];
    Chaos_p=GenerateChaos(circle);
    %Record the number of success or failure
    ns = [];
    nf = [];
    pfit = ones(1, 12);
    normfit=zeros(1,12);
    LEP=50;
    
    
    %Main body which was developed by the authors
    time = 1;
    
    % The total number of runs
    totalTime = 30;
    
    while time <= totalTime
        
        radius=0.0001;
        
        outcome = [];
        
        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=feval(fhd,popold',varargin{:});
        valParents=valParents';
        
        c = 1/10;
        p = 0.05;
        
        CRm = 0.5;
        Fm = 0.5;
        
        Afactor = 1;
        
        archive.NP = Afactor * popsize;     % the maximum size of the archive
        archive.pop = zeros(0, n);          % the solutions stored in te archive
        archive.funvalues = zeros(0, 1);    % the function value of the archived solutions
        
        %% the values and indices of the best solutions
        [valBest, indBest] = sort(valParents, 'ascend');
        
        FES = 0; k_circle=1;
        
        CR_circle=zeros(1,circle);
        F_circle=zeros(1,circle);
        
        while FES < n*10000 %& min(fit)>error_value(problem)
            %% Chaotic local search
            if Strategy>0;
                
                lb=lu(1,:);ub=lu(2,:);
                
                z=0.5;
                g=k_circle;
                
                switch Strategy
                    % 随机混沌
                    case 1;
                        j=randi([1,12],1,1);
%                         j=14;
                        temp_X=popold(indBest(1),:)+radius*(ub-lb).*(Chaos_p(j,g)-z);%r=0.1
                        
                        tp = temp_X>ub;
                        tm = temp_X<lb;
                        temp_X=(temp_X.*(~(tp+tm)))+ub.*tp+lb.*tm;
                        
                        fit_temp=feval(fhd,temp_X',varargin{:});
                        fit_temp=fit_temp';
                        
                        if fit_temp<valBest(1);
                            popold(indBest(1),:)=temp_X;
                            valBest(1)=fit_temp;
                        end
                        
                        %并联混沌
                    case 2;
                        temp_X=rand(12,30);
                        for j=1:12;
                            temp_X(j,:)=popold(indBest(1),:)+radius*(ub-lb)*(Chaos_p(j,g)-z);
                            
                            tp = temp_X(j,:)>ub;
                            tm = temp_X(j,:)<lb;
                            temp_X(j,:)=(temp_X(j,:).*(~(tp+tm)))+ub.*tp+lb.*tm;
                        end
                        fit_temp=feval(fhd,temp_X',varargin{:});    
                        fit_temp=fit_temp';
                        
                        [num_Fbest, num_X]=min(fit_temp);
                        
                        if num_Fbest<valBest(1);
                            popold(indBest(1),:)=temp_X(num_X,:);
                            valBest(1)=num_Fbest;
                        end
                        
                        % stochasitc universal sampling
                    case 3;
                        rr = rand;
                        j=1;
                        partsum = 0;
                        for i=1:12;
                            normfit(i) = pfit(i) / sum(pfit);
                        end
                        while partsum<rr;
                            partsum = partsum + normfit(j);
                            j=j+1;
                        end
                        select=j-1;
                        
                        lpcount=[];
                        npcount=[];
                        temp_X=popold(indBest(1),:)+radius*(ub-lb)*(Chaos_p(select,g)-z);%r=0.1
                        
                        tp = temp_X>ub;
                        tm = temp_X<lb;
                        temp_X=(temp_X.*(~(tp+tm)))+ub.*tp+lb.*tm;
                        
                        fit_temp=feval(fhd,temp_X',varargin{:});   
                        fit_temp=fit_temp';
                        
                        if fit_temp<valBest(1);
                            popold(indBest(1),:)=temp_X;
                            valBest(1)=fit_temp;
                            
                            tlpcount = zeros(1, 12);
                            tlpcount(select)=1;
                            lpcount=[lpcount;tlpcount];
                            
                            tnpcount = ones(1, 12);
                            tnpcount(select)=0;
                            npcount=[npcount;tnpcount];
                        else
                            tlpcount = zeros(1,12);
                            lpcount=[lpcount;tlpcount];
                            tnpcount = ones(1,12);
                            npcount=[npcount;tnpcount];
                        end
                        
                        ns = [ns; sum(lpcount, 1)];
                        nf = [nf; sum(npcount, 1)];
                        
                        %success and failure memory
                        if k_circle+1 >= LEP;
                            for i = 1 : 12
                                if (sum(ns(:, i)) + sum(nf(:, i))) == 0
                                    pfit(i) = 0.01;%to avoid the possible null success rates
                                else
                                    pfit(i) = sum(ns(:, i)) / (sum(ns(:, i)) + sum(nf(:, i))) + 0.01;
                                end
                            end
                            if size(ns,1)>LEP;
                                ns(1, :) = [];
                            end
                            if size(nf,1)>LEP;
                                nf(1, :) = [];
                            end
                        end    
                end
%                 radius=radius*0.986;
            end
            %%
            CR_circle(k_circle)=CRm;
            F_circle(k_circle)=Fm;
            
            
            pop = popold; % the old population becomes the current population
            
            if FES > 1 && ~isempty(goodCR) && sum(goodF) > 0 % If goodF and goodCR are empty, pause the update
                CRm = (1 - c) * CRm + c * mean(goodCR);
                Fm = (1 - c) * Fm + c * sum(goodF .^ 2) / sum(goodF); % Lehmer mean
            end
            
            % 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, CR] = randFCR(popsize, CRm, 0.1, Fm, 0.1);
            popAll = [pop; archive.pop];
            [r1, r2] = gnR1R2(popsize, size(popAll, 1));
            
            % 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 == == == == == == ==
            vi = pop + F(:, ones(1, n)) .* (pbest - pop + pop(r1, :) - popAll(r2, :));
            vi = boundConstraint(vi, pop, lu);
            
            % == == == == == == == 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;
            ui = vi; ui(mask) = pop(mask);
            
            valOffspring=feval(fhd,ui',varargin{:});   
            valOffspring=valOffspring';
            
            switch Strategy
                case 0
                    FES=FES+popsize;
                case 1
                    FES=FES+popsize+1;
                case 2
                    FES=FES+popsize+12;
                case 3
                    FES=FES+popsize+1;
            end
            
            % == == == == == == == == == == == == == == == Selection == == == == == == == == == == == == ==
            % I == 1: the parent is better; I == 2: the offspring is better
            [valParents, I] = min([valParents, valOffspring], [], 2);
            popold = pop;
            
            archive = updateArchive(archive, popold(I == 1, :), valParents(I == 1));
            
            popold(I == 2, :) = ui(I == 2, :);
            
            goodCR = CR(I == 2);
            goodF = F(I == 2);
            
            k_circle=k_circle+1;
            
            [valBest indBest] = sort(valParents, 'ascend');
            %             display(['The best optimal value of DE is : ', num2str(min(valParents))]);
            outcome = [outcome; min(valParents)];
        end
        
        
        FbestChart = [FbestChart outcome];
        time
        time = time + 1;
        
    end
%         xlswrite(['C:\Users\周天楽\Desktop\JADE 13 D50.xlsx'],FbestChart,['Sheet',num2str(problem)]);
%         xlswrite(['C:\Users\周天楽\Desktop\JADE 13 D50.xlsx'],mean(FbestChart(circle,:)),['Sheet',num2str(problem)],'A5002');%3002 2973 2681
%         xlswrite(['C:\Users\周天楽\Desktop\JADE 13 D50.xlsx'],std(FbestChart(circle,:)),['Sheet',num2str(problem)],'B5002');
    
end
toc;
