
% Gravitational Search Algorithm.
function [Fbest,BestChart,BestFitness,farmlayout,farmlayout_NA]=GSA(N, ElitistCheck, min_flag, Rpower,limit,p,wf,max_it,t,NA_type,tn,ws_folder,algorithmDir)

%V:   Velocity.
%a:   Acceleration.1
%M:   Mass.  Ma=Mp=Mi=M;
%dim: Dimension of the test function.
%N:   Number of agents.
%X:   Position of agents. dim-by-N matrix.
%R:   Distance between agents in search space.
%[low-up]: Allowable range for search space.
%Rnorm:  Norm in eq.8.
%Rpower: Power of R in eq.7.



Rnorm=2;
FES=0;

count_list = zeros(1,N);
success_list = zeros(1,N);
history_fitness = zeros(1,N);

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);


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

%random initialization for agents.
% X=initialization(D,N,up,down);
[X,lu] = windfarm_init(N, wf.turbine_num,wf);
X = windfarm_constraint(X, wf.NA_loc, D,down,up);
%create the best so far chart and average fitnesses chart.

V=zeros(N,D);
a = zeros(N,D);

pop_power = zeros(N,wf.cols * wf.rows);
% while FES <  D*10000
for iteration = 1:max_it
    % while iteration <= max_it
    %Checking allowable range. The time complexity is O(n)


    %Evaluation of agents. The time complexity is O(n)
    %     fitness=feval(fhd,X',P_index);
    
%     X1 = zeros(size(X));
%     for llll = 1:N
%         X1(llll,:)  = 2:26;
%     end
    
    [fitness,power_order,lp_power_accum]= wf_fitness(wf,X);
%     for i = 1:N
%         pop_power(i,power_order(i,:))  = pop_power(i,power_order(i,:)) + (lp_power_accum / sum(lp_power_accum))';
%     end
%     X = repalce_worst(X, N, power_order,wf.NA_loc,pop_power);

    % The difference of fitness. The time complexity is O(n)
    if iteration == 1
        diff_fitness = zeros(1,N);
    else
        diff_fitness = history_fitness - fitness;
    end
    % update the failed individual counter. The time complexity is O(n)

    int = diff_fitness > 0;
    for l = 1 : N
        if int(l)
            count_list(l) = count_list(l) + 1;
            success_list(l) = 0;
        else
            count_list(l) = 0;
            success_list(l) = success_list(l) + 1;
        end
    end

    % update the Bestfitness  best_X is index of best.
    % The time complexity is 2*O(n);
    if min_flag==1
        [best, ~]=min(fitness); %minimization.
        %         best_list = [best_list;best_X];
    else
        [best, best_index]=max(fitness); %maximization.
    end

    if iteration==1
        Fbest=best;
    end
    if min_flag==1
        if best<Fbest  %minimization.
            Fbest=best;
        end
    else
        if best>Fbest  %maximization
            Fbest=best;
        end
    end
%     disp([ 'NA ' num2str(NA_type) ' Turbine Num ' num2str(tn) ' Wind Profile ' num2str(wt) ' run: ' num2str(t) ' iteration: ' num2str(iteration) '    eta: ' num2str(Fbest / wf.power_total) ' Fitness: ' num2str(Fbest) ]);
%     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,iteration,(Fbest / wf.power_total),Fbest)
    %X = repalce_worst(X, N,power_order,wf.NA_loc);
    
    BestChart(iteration) = Fbest / wf.power_total;
    BestFitness=Fbest;
    [best_farmlayout,best_farmlayout_NA]  = gene_layout_by_indices_one(wf,X(best_index,:));
    farmlayout(iteration,:) = best_farmlayout;
    farmlayout_NA(iteration,:) = best_farmlayout_NA;

    %Calculation of M. eq.14-20 The time complexity is *O(4n) max min sum individual
    [M]=massCalculation(fitness,min_flag);

    % Self-adaptive Gravitational constant. The time complexity is O(2n)
    
    [G,count_list,success_list] = Gconstant(iteration,max_it,N,count_list,success_list,limit,p,a);
%     G=GSA_Gconstant(iteration,max_it);
    %Calculation of accelaration in gravitational field. eq.7-10,21.
    % The time complexity is O(n^3)
%     a =GSA_Gfield(M,X,G,Rnorm,Rpower,fitness,ElitistCheck,iteration,max_it,lu,wf);
    a = GSA_Gfield(M,X,G,Rnorm,Rpower,ElitistCheck,iteration,max_it);
    %Agent movement. eq.11-12 The time complexity O(n);

    [X,V]=move(X,a,V);
    X = windfarm_constraint(X, wf.NA_loc, D,down,up);
    
    FES=FES+N;

    history_fitness = fitness;

end

end
%iteration
% x = 1:3000;
% plot(x,count);


