function [Fbest,BestChart,BestFitness,farmlayout,farmlayout_NA]=AGPSO(popsize,wf,iterations,t,NA_type,tn,wt,algorithmDir,omega)


D = wf.turbine_num; down = 1;up = wf.rows *wf.cols;
%%  parameter settings for L-SHADE
p_best_rate = 0.11;
arc_rate = 1.4;
memory_size = 5;
problem_size = D;
pop_size = 18 * D;
max_pop_size = pop_size;
min_pop_size = 4.0;

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


disp(pop_size)
% Population initialization
[popold,lu] = windfarm_init(pop_size, wf.turbine_num,wf);

popsize = pop_size;

% Population evaluation
[popuFitness,power_order,lp_power_accum]= wf_fitness(wf,popold);
optimal = max(popuFitness);
for i = 1:popsize
    pop_power(i,power_order(i,:))  = pop_power(i,power_order(i,:)) + (lp_power_accum(i,:) / sum(lp_power_accum(i,:)));
end

% popuFitness = feval(fhd,popu',func);
% popuFitness = popuFitness';

vel = zeros(popsize, D);
pBest = popold;
pBestFit = popuFitness;
pBest_power = power_order;
[~, gBestId] = max(pBestFit);
gBest = pBest(gBestId,:);


flag = zeros(popsize,1);
sg = 7;
pm = 0.01;
sum_count = 3;
count = 0;
r = 4;
his_popu = zeros(r,popsize,D);


nfes = 0;
run_funcvals = [];


%%  parameter settings for Hybridization
First_calss_percentage=0.5;

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


[bsf_fit_var, bsf_index] = max(popuFitness);
bsf_solution = popu(bsf_index, :);

run_funcvals = [run_funcvals;ones(pop_size,1)*bsf_fit_var];


memory_sf = 0.5 .* ones(memory_size, 1);
memory_cr = 0.5 .* ones(memory_size, 1);
memory_pos = 1;

archive.NP = arc_rate * pop_size; % the maximum size of the archive
archive.pop = zeros(0, problem_size); % the solutions stored in te archive
archive.funvalues = zeros(0, 1); % the function value of the archived solutions

memory_1st_class_percentage = First_calss_percentage.* ones(memory_size, 1); % Class#1 probability for Hybridization


%% Initialize CMAES parameters
sigma = 0.5;          % coordinate wise standard deviation (step size)
xmean = rand(problem_size,1);    % objective variables initial point
mu = pop_size/2;               % number of parents/points for recombination
weights = log(mu+1/2)-log(1:mu)'; % muXone array for weighted recombination
mu = floor(mu);
weights = weights/sum(weights);     % normalize recombination weights array
mueff=sum(weights)^2/sum(weights.^2); % variance-effectiveness of sum w_i x_i

% Strategy parameter setting: Adaptation
cc = (4 + mueff/problem_size) / (problem_size+4 + 2*mueff/problem_size); % time constant for cumulation for C
cs = (mueff+2) / (problem_size+mueff+5);  % t-const for cumulation for sigma control
c1 = 2 / ((problem_size+1.3)^2+mueff);    % learning rate for rank-one update of C
cmu = min(1-c1, 2 * (mueff-2+1/mueff) / ((problem_size+2)^2+mueff));  % and for rank-mu update
damps = 1 + 2*max(0, sqrt((mueff-1)/(problem_size+1))-1) + cs; % damping for sigma usually close to 1

% Initialize dynamic (internal) strategy parameters and constants
pc = zeros(problem_size,1);
ps = zeros(problem_size,1);   % evolution paths for C and sigma
B = eye(problem_size,problem_size);                       % B defines the coordinate system
D = ones(problem_size,1);                      % diagonal D defines the scaling
C = B * diag(D.^2) * B';            % covariance matrix C
invsqrtC = B * diag(D.^-1) * B';    % C^-1/2
eigeneval = 0;                      % track update of B and D
chiN=problem_size^0.5*(1-1/(4*problem_size)+1/(21*problem_size^2));  % expectation of


%% main loop
Hybridization_flag=1; % Indicator flag if we need to Activate CMAES Hybridization
max_it = 200
for iter = 1:max_it

    popu = popold; % the old population becomes the current population
    [~, sorted_index] = sort(popuFitness, 'ascend');

    mem_rand_index = ceil(memory_size * rand(pop_size, 1));
    mu_sf = memory_sf(mem_rand_index);
    mu_cr = memory_cr(mem_rand_index);
    mem_rand_ratio = rand(pop_size, 1);

    %% for generating crossover rate
    cr = normrnd(mu_cr, 0.1);
    term_pos = find(mu_cr == -1);
    cr(term_pos) = 0;
    cr = min(cr, 1);
    cr = max(cr, 0);

    %% for generating scaling factor
    if(iter <= max_it/2)
        sf=0.45+.1*rand(pop_size, 1);
        pos = find(sf <= 0);

        while ~ isempty(pos)
            sf(pos)=0.45+0.1*rand(length(pos), 1);
            pos = find(sf <= 0);
        end
    else
        sf = mu_sf + 0.1 * tan(pi * (rand(pop_size, 1) - 0.5));

        pos = find(sf <= 0);

        while ~ isempty(pos)
            sf(pos) = mu_sf(pos) + 0.1 * tan(pi * (rand(length(pos), 1) - 0.5));
            pos = find(sf <= 0);
        end
    end
    sf = min(sf, 1);

    %% for generating Hybridization Class probability
    Class_Select_Index=(memory_1st_class_percentage(mem_rand_index)>=mem_rand_ratio);
    if(Hybridization_flag==0)
        Class_Select_Index=or(Class_Select_Index,~Class_Select_Index);%All will be in class#1
    end

    %%
    r0 = [1 : pop_size];
    popAll = [popu; archive.pop];
    [r1, r2] = gnR1R2(pop_size, size(popAll, 1), r0);

    pNP = max(round(p_best_rate * pop_size), 2); %% choose at least two best solutions
    randindex = ceil(rand(1, pop_size) .* pNP); %% select from [1, 2, 3, ..., pNP]
    randindex = max(1, randindex); %% to avoid the problem that rand = 0 and thus ceil(rand) = 0
    pbest = popu(sorted_index(randindex), :); %% randomly choose one of the top 100p% solutions

    vi=[];
    temp=[];
    disp(sum(Class_Select_Index))
    if(sum(Class_Select_Index)~=0)
        vi(Class_Select_Index,:) = popu(Class_Select_Index,:) + sf(Class_Select_Index, ones(1, problem_size)) .* (pbest(Class_Select_Index,:) - popu(Class_Select_Index,:) + popu(r1(Class_Select_Index), :) - popAll(r2(Class_Select_Index), :));
    end

    if(sum(~Class_Select_Index)~=0)
        for k=1:sum(~Class_Select_Index)
            temp(:,k) = xmean + sigma * B * (D .* randn(problem_size,1)); % m + sig * Normal(0,C)
        end
        vi(~Class_Select_Index,:) = temp';
    end

    if(~isreal(vi))
        Hybridization_flag=0;
        continue;
    end


%     vi = boundConstraint(vi, popu, lu);


    mask = rand(pop_size, problem_size) > cr(:, ones(1, problem_size)); % mask is used to indicate which elements of ui comes from the parent
    rows = (1 : pop_size)'; cols = floor(rand(pop_size, 1) * problem_size)+1; % choose one position where the element of ui doesn't come from the parent
    jrand = sub2ind([pop_size problem_size], rows, cols); mask(jrand) = false;
    ui = vi; ui(mask) = popu(mask);

%     children_fitness = feval(fhd, ui', func);
%     children_fitness = children_fitness';

     ui = windfarm_constraint(ui, wf.NA_loc, D,down,up);
    [offsPbestFitness,offs_power_order,offs_lp_power_accum]= wf_fitness(wf,ui);
    optimal = max(popuFitness);

    

    %%%%%%%%%%%%%%%%%%%%%%%% for out
    for i = 1 : pop_size
        nfes = nfes + 1;
        if (children_fitness(i) < bsf_fit_var && isreal(ui(i, :)) && sum(isnan(ui(i, :)))==0 && min(ui(i, :))>=-100 && max(ui(i, :))<=100)
            bsf_fit_var = children_fitness(i);
            bsf_solution = ui(i, :);
            bsf_index = i;
        end

        if nfes > max_nfes;
            break;
        end
    end
    %%%%%%%%%%%%%%%%%%%%%%%% for out

    [bsf_fit_var, bsf_index] = max(offsPbestFitness);
 
    bsf_solution = ui(bsf_index, :);

    run_funcvals = [run_funcvals;ones(pop_size,1)*bsf_fit_var];

    dif = abs(popuFitness - children_fitness);


    %% I == 1: the parent is better; I == 2: the offspring is better
    Child_is_better_index = (popuFitness > children_fitness);
    goodCR = cr(Child_is_better_index == 1);
    goodF = sf(Child_is_better_index == 1);
    dif_val = dif(Child_is_better_index == 1);
    dif_val_Class_1 = dif(and(Child_is_better_index,Class_Select_Index) == 1);
    dif_val_Class_2 = dif(and(Child_is_better_index,~Class_Select_Index) == 1);

    archive = updateArchive(archive, popold(Child_is_better_index == 1, :), popuFitness(Child_is_better_index == 1));

    [popuFitness, Child_is_better_index] = min([popuFitness, children_fitness], [], 2);

    popold = popu;
    popold(Child_is_better_index == 2, :) = ui(Child_is_better_index == 2, :);

    num_success_params = numel(goodCR);

    if num_success_params > 0
        sum_dif = sum(dif_val);
        dif_val = dif_val / sum_dif;

        %% for updating the memory of scaling factor
        memory_sf(memory_pos) = (dif_val' * (goodF .^ 2)) / (dif_val' * goodF);

        %% for updating the memory of crossover rate
        if max(goodCR) == 0 || memory_cr(memory_pos)  == -1
            memory_cr(memory_pos)  = -1;
        else
            memory_cr(memory_pos) = (dif_val' * (goodCR .^ 2)) / (dif_val' * goodCR);
        end

        if (Hybridization_flag==1)% if the Hybridization is activated
            memory_1st_class_percentage(memory_pos) = memory_1st_class_percentage(memory_pos)*L_Rate+ (1-L_Rate)*(sum(dif_val_Class_1) / (sum(dif_val_Class_1) + sum(dif_val_Class_2)));
            memory_1st_class_percentage(memory_pos)=min(memory_1st_class_percentage(memory_pos),0.8);
            memory_1st_class_percentage(memory_pos)=max(memory_1st_class_percentage(memory_pos),0.2);
        end

        memory_pos = memory_pos + 1;
        if memory_pos > memory_size;  memory_pos = 1; end
    end

    %% for resizing the population size
    plan_pop_size = round((((min_pop_size - max_pop_size) / max_nfes) * nfes) + max_pop_size);

    if pop_size > plan_pop_size
        reduction_ind_num = pop_size - plan_pop_size;
        if pop_size - reduction_ind_num <  min_pop_size; reduction_ind_num = pop_size - min_pop_size;end

        pop_size = pop_size - reduction_ind_num;
        for r = 1 : reduction_ind_num
            [valBest, indBest] = sort(popuFitness, 'ascend');
            worst_ind = indBest(end);
            popold(worst_ind,:) = [];
            popu(worst_ind,:) = [];
            popuFitness(worst_ind,:) = [];
            Child_is_better_index(worst_ind,:) = [];
        end

        archive.NP = round(arc_rate * pop_size);

        if size(archive.pop, 1) > archive.NP
            rndpos = randperm(size(archive.pop, 1));
            rndpos = rndpos(1 : archive.NP);
            archive.pop = archive.pop(rndpos, :);
        end

        %% update CMA parameters
        mu = pop_size/2;               % number of parents/points for recombination
        weights = log(mu+1/2)-log(1:mu)'; % muXone array for weighted recombination
        mu = floor(mu);
        weights = weights/sum(weights);     % normalize recombination weights array
        mueff=sum(weights)^2/sum(weights.^2); % variance-effectiveness of sum w_i x_i
    end

    %% CMAES Adaptation
    if(Hybridization_flag==1)
        % Sort by fitness and compute weighted mean into xmean
        [~, popindex] = sort(popuFitness);  % minimization
        xold = xmean;
        xmean = popold(popindex(1:mu),:)' * weights;  % recombination, new mean value

        % Cumulation: Update evolution paths
        ps = (1-cs) * ps ...
            + sqrt(cs*(2-cs)*mueff) * invsqrtC * (xmean-xold) / sigma;
        hsig = sum(ps.^2)/(1-(1-cs)^(2*nfes/pop_size))/problem_size < 2 + 4/(problem_size+1);
        pc = (1-cc) * pc ...
            + hsig * sqrt(cc*(2-cc)*mueff) * (xmean-xold) / sigma;

        % Adapt covariance matrix C
        artmp = (1/sigma) * (popold(popindex(1:mu),:)' - repmat(xold,1,mu));  % mu difference vectors
        C = (1-c1-cmu) * C ...                   % regard old matrix
            + c1 * (pc * pc' ...                % plus rank one update
            + (1-hsig) * cc*(2-cc) * C) ... % minor correction if hsig==0
            + cmu * artmp * diag(weights) * artmp'; % plus rank mu update

        % Adapt step size sigma
        sigma = sigma * exp((cs/damps)*(norm(ps)/chiN - 1));

        % Update B and D from C
        if nfes - eigeneval > pop_size/(c1+cmu)/problem_size/10  % to achieve O(problem_size^2)
            eigeneval = nfes;
            C = triu(C) + triu(C,1)'; % enforce symmetry
            if(sum(sum(isnan(C)))>0 || sum(sum(~isfinite(C)))>0 || ~isreal(C))
                Hybridization_flag=0;
                continue;
            end
            [B,D] = eig(C);           % eigen decomposition, B==normalized eigenvectors
            D = sqrt(diag(D));        % D contains standard deviations now
            invsqrtC = B * diag(D.^-1) * B';
        end

    end

end %%%%%%%%nfes

%% Violation Checking
if(max(bsf_solution)>100)
    fprintf('%d th run, Above Max\n', run_id)
end

if(min(bsf_solution)<-100)
    fprintf('%d th run, Below Min\n', run_id)
end

if(~isreal(bsf_solution))
    fprintf('%d th run, Complix\n', run_id)
end

bsf_error_val = abs(bsf_fit_var - optimum);
if bsf_error_val < val_2_reach
    bsf_error_val = 0;
end

if(sum(isnan(bsf_solution))>0)
    fprintf('%d th run, NaN\n', run_id)
end

fprintf('%d th run, best-so-far error value = %1.8e\n', run_id , bsf_error_val)
outcome = [outcome bsf_error_val];

%% From Noor Code ( print files )
errorVals= [];
for w = 1 : progress
    bestold = run_funcvals(RecordFEsFactor(w) * max_nfes) - optimum;
    if abs(bestold)>1e-8
        errorVals(w)= abs(bestold);
    else
        bestold=0;
        errorVals(w)= bestold;
    end
end
allerrorvals(:, run_id) = errorVals;
end

