% GSA code v1.1.
% Generated by Esmat Rashedi, 2010.
% "	E. Rashedi, H. Nezamabadi-pour and S. Saryazdi,
% “GSA: A Gravitational Search Algorithm? Information sciences, vol. 179,
% no. 13, pp. 2232-2248, 2009."

% Gravitational Search Algorithm.
function [Fbest,BestChart]=GSA(fhd,N,max_it,ElitistCheck,min_flag,Rpower,P_index,limit,p)

%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;
iteration=1;
count_list = zeros(1,N);
success_list = zeros(1,N);
history_fitness = zeros(1,N);
% G0 = 100;
% G = repmat(G0,N,1);
% engineering funchtion
% [bl,bu,D,fobj]=Get_Functions_details(P_index);
% X= repmat(bl, N, 1) + rand(N, D).*(repmat(bu - bl, N, 1));
%get allowable range and dimension of the test function.
% [down,up,dim]=test_functions_range(F_index);

D = 30;down = -100;up = 100;

%random initialization for agents.
X=initialization(D,N,up,down);

%create the best so far chart and average fitnesses chart.
BestChart=zeros(3000,1);
V=zeros(N,D);
a = zeros(N,D);


while FES <  D*10000 
    % while iteration <= max_it
    %Checking allowable range. The time complexity is O(n)
    X=space_bound(X,up,down);

    %Evaluation of agents. The time complexity is O(n)
    fitness=feval(fhd,X',P_index);
    
    % 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, ~]=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(['iteration: ' num2str(iteration) '    fitness: ' num2str(Fbest) ' problem: ' num2str(P_index)]);
    
    BestChart(iteration) = Fbest;
    %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);   

    %Calculation of accelaration in gravitational field. eq.7-10,21. 
    % The time complexity is O(n^3)
    a = 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);
    
    FES=FES+N;
    iteration=iteration+1;
    history_fitness = fitness;
    
end

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

