function [Fbest,BestChart,BestFitness,farmlayout,farmlayout_NA] = CLPSO_wf(popuSize,wf,max_it,t,NA_type,tn,wt,algorithmDir)
optimalChart = zeros(max_it,1);

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


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

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

% Population initialization
[popu,lu] = windfarm_init(popuSize, wf.turbine_num,wf);
[popuFitness,power_order]= wf_fitness(wf,popu);
popuFitness = popuFitness';
% Population evaluation
nFES=0;

optimal = max(popuFitness);
BestChart(1,1) =optimal;


%%
% Local parameters setting
% Velocity  initialization
mV = 0.25 * (lu(2,:) - lu(1,:));
maxVel = repmat(mV, popuSize,1);
minVel = -maxVel;

vel = minVel + 2 .* maxVel .* rand(popuSize, D);
% Personal best initialization
pBest = popu;
pBestFit = popuFitness;
% Global best initialization
[gBestFit,gBestId]=max(pBestFit);
gBest = pBest(gBestId,:);
% Inertia weight initialization
omega0 = 0.9;
omega1 = 0.2;
% Number of generations of not improving
flag = zeros(popuSize,1);
% m dimensions are randomly chosen to learn from the gbest
m = 5 .* ones(popuSize,1);
% Assign pBest individual for each dimension
ai=zeros(popuSize,D);
pBestSelected=1:popuSize;
pBestSelected=repmat(pBestSelected',1,D);
[pBestSelected, ai] = AssignPBest(popuSize, D, ai, pBestFit, pBestSelected, m);


for iter = 1:max_it
    %% Search machine
    
    searchMachine = @CLPSO;
    inerOmega = omega0 - (omega0 - omega1) * iter / max_it;
    [popu, vel, flag] = searchMachine(popu, vel, mV, pBest, pBestFit, pBestSelected, ai, gBest, inerOmega, flag, m);

    popu = windfarm_constraint(popu, wf.NA_loc, D,down,up);
    %% Boundary detection + Population evaluation + pBest update
    indiFitness = zeros(popuSize,D);
    for i = 1:popuSize
        if (sum(popu(i,:) > lu(2,:) ) + sum(popu(i,:) < lu(1,:))) == 0
            
            [indiFitness(i,:),~,~]= wf_fitness(wf,popu(i,:));

            
            nFES = nFES + 1;
            
            if indiFitness(i) > optimal
                optimal = indiFitness(i);
     
            else
                optimalChart = [optimalChart;optimal];
            end

            if pBestFit(i) >= indiFitness(i)
                flag(i) = flag(i)+1;
            else
                pBest(i,:) = popu(i,:);
                pBestFit(i) = indiFitness(i);
            end

            if pBestFit(i) > gBestFit
                gBestFit = pBestFit(i);
                gBest  = pBest(i,:);
            end
        end
    end

    optimal = gBestFit;
    BestChart(iter) = optimal / wf.power_total;
    BestFitness(iter) = optimal;
    Fbest  = optimal;
    [best_farmlayout,best_farmlayout_NA]  = gene_layout_by_indices_one(wf,gBest);
    farmlayout(iter,:) = best_farmlayout;
    farmlayout_NA(iter,:) = best_farmlayout_NA;


%     fprintf('%s NA %d Turbine Num:%d Wind %s run: %d iteration: %d  eta %f fitness %f\n',algorithmDir,NA_type,tn,wt,t,iter,(Fbest / wf.power_total),Fbest)
    
end