function [Fbest,BestChart,BestFitness,farmlayout,farmlayout_NA]=CJADE(popsize,wf,iterations,t,NA_type,tn,Strategy,circle,Chaos_p,ws_folder,algorithmDir)

radius=0.01;
LEP=24;
goodCR = zeros(0);
goodF = zeros(0);

D = wf.turbine_num; down = 1;up = wf.rows *wf.cols;

BestChart=zeros(iterations,1);
BestFitness=zeros(iterations,1);
farmlayout = zeros(iterations,wf.rows*wf.cols);
farmlayout_NA = zeros(iterations,wf.rows*wf.cols);

rand('seed', sum(100 * clock));

[popold,lu] = windfarm_init(popsize, wf.turbine_num,wf);
[valParents,~, ~]= wf_fitness(wf,popold);

%% parameters
n=D;
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

ns = [];
nf = [];
pfit = ones(1, 12);
normfit=zeros(1,12);

%% the values and indices of the best solutions
[valBest, indBest] = sort(valParents, 'descend');

FES = 0; k_circle=1;

CR_circle=zeros(1,circle);
F_circle=zeros(1,circle);
success_num =zeros(LEP,12);
fail_num = zeros(LEP,12);
r_count = 0;
invere_flag = 0;
temp_lep = 1;
an = n;
for iter = 1:iterations
    %% Chaotic local search
    if Strategy>0
        lb=lu(1,:);
        ub=lu(2,:);
        randPopuList = randperm(popsize);
        randPopuList = setdiff(randPopuList,1,'stable');
        indiR1 = popold(randPopuList(1),:);
        indiR2 = popold(randPopuList(2),:);
        z=0.5;
        g=k_circle;
        radius = 0.5 * (1-iter/iterations);
        switch Strategy
            case 1
                %% CJADE-R
                j=randi([1,12],1,1);
                temp_X=popold(indBest(1),:)+radius*(ub-lb).*(Chaos_p(j,g)-z);%r=0.1
                temp_X = windfarm_constraint(temp_X, wf.NA_loc, D,down,up);
                [fit_temp,power_order]= wf_fitness(wf,temp_X);

                if fit_temp>valBest(1)
                    r_count= r_count +1;
                    popold(indBest(1),:)=temp_X;
                    valBest(1)=fit_temp;
                end

            case 2
                %% CJADE-P
                temp_X=rand(12,n);
                for j=1:12
                    temp_X(j,:)=popold(indBest(1),:)+radius*(ub-lb)*(Chaos_p(j,g)-z);
                    temp_X = windfarm_constraint(temp_X, wf.NA_loc, D,down,up);
                    [fit_temp,power_order]= wf_fitness(wf,temp_X);

                end

                [num_Fbest, num_X]=max(fit_temp);

                if num_Fbest>valBest(1)
                    r_count= r_count +1;
                    popold(indBest(1),:)=temp_X(num_X,:);
                    valBest(1)=num_Fbest;
                end

            case 4
                % Stochastic universal sampling
                rr = rand;
                
                %  cumulative probability
                normfit1 = cumsum(pfit)/sum(pfit);
                index = find(rr<normfit1);              
                select = index(1);
                
%                 x_index = floor(popold(indBest(1),:) -1 / wf.rows);
%                 y_index = floor(popold(indBest(1),:) -1  - x_index* wf.rows);
%                 x_index = x_index + radius*(wf.rows)*(Chaos_p(select,g)-z);
%                 y_index = y_index + radius*(wf.cols)*(Chaos_p(select,g)-z);
%                 temp_X =  x_index *12+y_index +1;
                temp_X=popold(indBest(1),:)+radius*(ub-lb)*(Chaos_p(select,g)-z);
                
                temp_X = windfarm_constraint(temp_X, wf.NA_loc, D,down,up);
                [fit_temp,power_order]= wf_fitness(wf,temp_X);
                
                
                if fit_temp>valBest(1)
                    r_count= r_count +1;
                    
                    popold(indBest(1),:)=temp_X;
                    valBest(1)=fit_temp;
                    success_num(temp_lep,select) = 1;
                    fail_num(temp_lep,~select)=1;



                else
                    success_num(temp_lep,:) = 0;
                    fail_num(temp_lep,:)=1;
                end
                ns = success_num;
                nf = fail_num;

                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

                temp_lep = temp_lep +1;
                if temp_lep > LEP
                    temp_lep = 1;
                end

            case 3
                % Stochastic universal sampling
                rr = rand;
                normfit1 = cumsum(pfit)/sum(pfit);
                index = find(rr<normfit1);              
                select = index(1);

                lpcount=[];
                temp_X=popold(indBest(1),:)+(indiR1-indiR2)*Chaos_p(select,g)*0.01;%r=0.1

                temp_X = windfarm_constraint(temp_X, wf.NA_loc, D,down,up);
                [fit_temp,~]= wf_fitness(wf,temp_X);

                
                if fit_temp>valBest(1)
                    popold(indBest(1),:)=temp_X;
                    r_count= r_count +1;
                    tlpcount = zeros(1, 12);
                    tlpcount(select)=valBest(1)-fit_temp;
                    success_num(temp_lep,select)=1;
                    lpcount=[lpcount;tlpcount];

                    valBest(1)=fit_temp;
                end
                temp_lep = temp_lep+1;
                ns =[ns; sum(lpcount, 1)];
            
                %success and failure memory
                if k_circle+1 >= LEP
                    for i = 1 : 12
                        if sum(sum(ns, 1))== 0
                            pfit(i) = 1/12;%to avoid the possible null success rates
                        else
                            pfit(i) = sum(ns(:, i)) / sum(sum(ns, 1))  + 1/12;
                        end
                    end
                    if temp_lep >LEP
                        temp_lep = 1;
                    end
                    if size(ns,1)>LEP
                         ns(1, :) = [];
                    end


                end
        end



    end
    %%
    CR_circle(k_circle)=CRm;
    F_circle(k_circle)=Fm;
%     [popold(indBest(1),:),valBest(1),pfit,temp_lep] = Chaotic_selection(popold(indBest(1),:),valBest(1),down,up,Chaos_p,wf,D,pfit,temp_lep,Strategy,iter,iterations,success_num,fail_num,LEP);

    pop = popold; % the old population becomes the current population

    if  ~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 = windfarm_constraint(vi, wf.NA_loc, D,down,up);
    %
    %             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);


    ui = windfarm_constraint(ui, wf.NA_loc, D,down,up);
    [valOffspring,~]= wf_fitness(wf,ui);


    % == == == == == == == == == == == == == == == Selection == == == == == == == == == == == == ==
    % I == 1: the parent is better; I == 2: the offspring is better
    [valParents, I] = max([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, 'descend');

    BestChart(iter) = valBest(1) / wf.power_total;
    BestFitness(iter) = valBest(1);
    Fbest  = valBest(1);
    [best_farmlayout,best_farmlayout_NA]  = gene_layout_by_indices_one(wf,popold(indBest(1),:));
    farmlayout(iter,:) = best_farmlayout;
    farmlayout_NA(iter,:) = best_farmlayout_NA;

    %             display(['The best optimal value of DE is : ', num2str(min(valParents))]);
    
%     fprintf('%s NA %d Turbine Num:%d Wind %s run: %d iteration: %d  eta %f fitness %f\n',algorithmDir,NA_type,tn,ws_folder,t,iter,(Fbest / wf.power_total),Fbest)
end


end