%**************************************************************************************************
%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;
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 = 1
    
    % Define the dimension of the problem
    n = 30;
    % Define the size of population
    popsize = 100;
    
    switch problem
        
        case 1
            
            % lu: define the upper and lower bounds of the variables
            lu = [-100 * ones(1, n); 100 * ones(1, n)];
            % Load the data for this test function
            load sphere_func_data
            A = []; M = []; a = []; alpha = []; b = [];
            
        case 2
            
            lu = [-100 * ones(1, n); 100 * ones(1, n)];
            load schwefel_102_data
            A = []; M = []; a = []; alpha = []; b = [];
            
        case 3
            
            lu = [-100 * ones(1, n); 100 * ones(1, n)];
            load high_cond_elliptic_rot_data
            A = []; a = []; alpha = []; b = [];
            
            if n == 2, load elliptic_M_D2,
            elseif n == 10, load elliptic_M_D10,
            elseif n == 30, load elliptic_M_D30,
            elseif n == 50, load elliptic_M_D50,
            end
            
        case 4
            
            lu = [-100 * ones(1, n); 100 * ones(1, n)];
            load schwefel_102_data
            A = []; M = []; a = []; alpha = []; b = [];
            
        case 5
            
            lu = [-100 * ones(1, n); 100 * ones(1, n)];
            load schwefel_206_data
            M = []; a = []; alpha = []; b = [];
            
        case 6
            
            lu = [-100 * ones(1, n); 100 * ones(1, n)];
            load rosenbrock_func_data
            A = []; M = []; a = []; alpha = []; b = [];
            
        case 7
            
            lu = [0 * ones(1, n); 600 * ones(1, n)];
            load griewank_func_data
            A = []; a = []; alpha = []; b = [];
            
            if n == 2, load griewank_M_D2,
            elseif n == 10, load griewank_M_D10,
            elseif n == 30, load griewank_M_D30,
            elseif n == 50, load griewank_M_D50,
            end
            
        case 8
            
            lu = [-32 * ones(1, n); 32 * ones(1, n)];
            load ackley_func_data
            A = []; a = []; alpha = []; b = [];
            
            if n == 2, load ackley_M_D2,
            elseif n == 10, load ackley_M_D10,
            elseif n == 30, load ackley_M_D30,
            elseif n == 50, load ackley_M_D50,
            end
            
        case 9
            
            lu = [-5 * ones(1, n); 5 * ones(1, n)];
            load rastrigin_func_data
            A = []; M = []; a = []; alpha = []; b = [];
            
        case 10
            
            lu = [-5 * ones(1, n); 5 * ones(1, n)];
            load rastrigin_func_data
            A = []; a = []; alpha = []; b = [];
            if n == 2, load rastrigin_M_D2,
            elseif n == 10, load rastrigin_M_D10,
            elseif n == 30, load rastrigin_M_D30,
            elseif n == 50, load rastrigin_M_D50,
            end
            
        case 11
            
            lu = [-0.5 * ones(1, n); 0.5 * ones(1, n)];
            load weierstrass_data
            A = []; a = []; alpha = []; b = [];
            if n == 2, load weierstrass_M_D2,
            elseif n == 10, load weierstrass_M_D10,
            elseif n == 30, load weierstrass_M_D30,
            elseif n == 50, load weierstrass_M_D50,
            end
            
        case 12
            
            lu = [-pi * ones(1, n); pi * ones(1, n)];
            load schwefel_213_data
            A = []; M = []; o = [];
            
        case 13
            
            lu = [-3 * ones(1, n); 1 * ones(1, n)];
            load EF8F2_func_data
            A = []; M = []; a = []; alpha = []; b = [];
            
        case 14
            
            lu = [-100 * ones(1, n); 100 * ones(1, n)];
            load E_ScafferF6_func_data
            if n == 2, load E_ScafferF6_M_D2,
            elseif n == 10, load E_ScafferF6_M_D10,
            elseif n == 30, load E_ScafferF6_M_D30,
            elseif n == 50, load E_ScafferF6_M_D50,
            end
            A = []; a = []; alpha = []; b = [];
            
        case 15
            
            lu = [-5 * ones(1, n); 5 * ones(1, n)];
            load hybrid_func1_data
            A = []; M = []; a = []; alpha = []; b = [];
            
        case 16
            
            lu = [-5 * ones(1, n); 5 * ones(1, n)];
            load hybrid_func1_data
            if n == 2, load hybrid_func1_M_D2,
            elseif n == 10, load hybrid_func1_M_D10,
            elseif n == 30, load hybrid_func1_M_D30,
            elseif n == 50, load hybrid_func1_M_D50,
            end
            A = []; a = []; alpha = []; b = [];
            
        case 17
            
            lu = [-5 * ones(1, n); 5 * ones(1, n)];
            load hybrid_func1_data
            if n == 2, load hybrid_func1_M_D2,
            elseif n == 10, load hybrid_func1_M_D10,
            elseif n == 30, load hybrid_func1_M_D30,
            elseif n == 50, load hybrid_func1_M_D50,
            end
            A = []; a = []; alpha = []; b = [];
            
        case 18
            
            lu = [-5 * ones(1, n); 5 * ones(1, n)];
            load hybrid_func2_data
            if n == 2, load hybrid_func2_M_D2,
            elseif n == 10, load hybrid_func2_M_D10,
            elseif n == 30, load hybrid_func2_M_D30,
            elseif n == 50, load hybrid_func2_M_D50,
            end
            A = []; a = []; alpha = []; b = [];
            
        case 19
            
            lu = [-5 * ones(1, n); 5 * ones(1, n)];
            load hybrid_func2_data
            if n == 2, load hybrid_func2_M_D2,
            elseif n == 10, load hybrid_func2_M_D10,
            elseif n == 30, load hybrid_func2_M_D30,
            elseif n == 50, load hybrid_func2_M_D50,
            end
            A = []; a = []; alpha = []; b = [];
            
        case 20
            
            lu = [-5 * ones(1, n); 5 * ones(1, n)];
            load hybrid_func2_data
            if n == 2, load hybrid_func2_M_D2,
            elseif n == 10, load hybrid_func2_M_D10,
            elseif n == 30, load hybrid_func2_M_D30,
            elseif n == 50, load hybrid_func2_M_D50,
            end
            A = []; a = []; alpha = []; b = [];
            
        case 21
            
            lu = [-5 * ones(1, n); 5 * ones(1, n)];
            load hybrid_func3_data
            if n == 2, load hybrid_func3_M_D2,
            elseif n == 10, load hybrid_func3_M_D10,
            elseif n == 30, load hybrid_func3_M_D30,
            elseif n == 50, load hybrid_func3_M_D50,
            end
            A = []; a = []; alpha = []; b = [];
            
        case 22
            
            lu = [-5 * ones(1, n); 5 * ones(1, n)];
            load hybrid_func3_data
            if n == 2, load hybrid_func3_HM_D2,
            elseif n == 10, load hybrid_func3_HM_D10,
            elseif n == 30, load hybrid_func3_HM_D30,
            elseif n == 50, load hybrid_func3_HM_D50,
            end
            A = []; a = []; alpha = []; b = [];
            
        case 23
            
            lu = [-5 * ones(1, n); 5 * ones(1, n)];
            load hybrid_func3_data
            if n == 2, load hybrid_func3_M_D2,
            elseif n == 10, load hybrid_func3_M_D10,
            elseif n == 30, load hybrid_func3_M_D30,
            elseif n == 50, load hybrid_func3_M_D50,
            end
            A = []; a = []; alpha = []; b = [];
            
        case 24
            
            lu = [-5 * ones(1, n); 5 * ones(1, n)];
            load hybrid_func4_data
            if n == 2, load hybrid_func4_M_D2,
            elseif n == 10, load hybrid_func4_M_D10,
            elseif n == 30, load hybrid_func4_M_D30,
            elseif n == 50, load hybrid_func4_M_D50,
            end
            A = []; a = []; alpha = []; b = [];
            
        case 25
            
            lu = [2 * ones(1, n); 5 * ones(1, n)];
            load hybrid_func4_data
            if n == 2, load hybrid_func4_M_D2,
            elseif n == 10, load hybrid_func4_M_D10,
            elseif n == 30, load hybrid_func4_M_D30,
            elseif n == 50, load hybrid_func4_M_D50,
            end
            A = []; a = []; alpha = []; b = [];
            
    end
    
    % Choose the search strategy
    Strategy=0;
    switch Strategy
        case 0
            circle=3000; % 0: traditional JADE
        case 1
            circle=2971; % 1: randomly 
        case 2
            circle=2679; % 2: parallelly 
        case 3
            circle=2971; % 3: memory-based
    end
    
    % Record the best results
    FbestChart= [];
    %Generate chaotic sequences
    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.01;
        
        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 = benchmark_func(popold, problem, o, A, M, a, alpha, b);
        
        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 
            %% Chaotic local search
            if Strategy>0;
                
                lb=lu(1,:);ub=lu(2,:);
                
                z=0.5;
                g=k_circle;
                
                % CLS strategies 
                switch Strategy
                    % randomly
                    case 1;
                        j=randi([1,12],1,1);
                        %                         j=14;% Single chaos
                        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 = benchmark_func(temp_X, problem, o, A, M, a, alpha, b);
                        
                        if fit_temp<valBest(1);
                            popold(indBest(1),:)=temp_X;
                            valBest(1)=fit_temp;
                        end
                        
                        % parallelly
                    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 = benchmark_func(temp_X, problem, o, A, M, a, alpha, b);
                        
                        
                        [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 = benchmark_func(temp_X, problem, o, A, M, a, alpha, b);
                        
                        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.979;

            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 = benchmark_func(ui, problem, o, A, M, a, alpha, b);
            
            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\1.xlsx'],FbestChart,['Sheet',num2str(problem)]);
%     xlswrite(['C:\Users\周天楽\Desktop\1.xlsx'],mean(FbestChart(circle,:)),['Sheet',num2str(problem)],'A4953');%3002 2973 2681
%     xlswrite(['C:\Users\周天楽\Desktop\1.xlsx'],std(FbestChart(circle,:)),['Sheet',num2str(problem)],'B4953');
    
end
toc;
