% the main.m
clc;
clear all
%%
% Global parameters setting
format long;
addpath('Benchmark_Solar_Cell');
fun = @evaluate_normal_fitness;
algorithmDir = 'GWO';
popuSize = 30;
FES = 50000;


problemIndex = [1:6];
maxIter = FES / popuSize;
maxIter = ceil(maxIter);

totaltime = 30;
analysisChart = [];
addpath('../');
for f = problemIndex
    resultChart = [];
    time = 1;
    %     nFES = 0;
    func_flag = f;
    PV_Xrange;
    problem = func_flag;

    while time <= totaltime
        tic
        optimalChart = [];
        %         IAEI2 = [];
        IAEI3 = [];
        rand('seed', sum(100 * clock));
        D = length(Xmax);
        ProblemSize = D;
        lu = [Xmin;Xmax];


        % Population initialization
        popu = repmat(lu(1, :), popuSize, 1) + rand(popuSize, D) .* (repmat(lu(2, :) - lu(1, :), popuSize, 1));

        % Population evaluation
        nFES=0;
        IAEI1 = [];
        for i=1:popuSize
            [popuFitness(i,:),IAEI1(i,:)] =  fun(popu(i,:),func_flag);
            nFES = nFES + 1;
            FE_best{func_flag}(time,nFES) =  popuFitness(i,:);
        end
        %         IAEI2 = [IAEI2;IAEI1];
        optimal = min(popuFitness);
        optimalChart = [optimalChart;repmat(optimal,popuSize,1)];
        record(1) = min(popuFitness);
        %         [popuFitness] = Evaluation(popu, benchmark, f);
        Fitness = popuFitness;
        %         Fitness = Fitness';
        %         popuFitness = popuFitness';

        Alpha_score = 1e+30;
        Alpha_pos = zeros(1, ProblemSize);
        Beta_pos = zeros(1,ProblemSize);
        Beta_score = 1e+30;
        Delta_pos = zeros(1,ProblemSize);
        Delta_score = 1e+30;
        for i = 1 : popuSize
            if Fitness(i) < Alpha_score
                Alpha_score = Fitness(i);
                Alpha_pos = popu(i, :);
            end
            if Fitness(i)>Alpha_score && Fitness(i)<Beta_score
                Beta_score = Fitness(i); % Update beta
                Beta_pos = popu(i,:);
            end
            if Fitness(i)>Alpha_score && Fitness(i)>Beta_score && Fitness(i)<Delta_score
                Delta_score = Fitness(i); % Update delta
                Delta_pos = popu(i,:);
            end
            %             if nFES == FES; break; end
        end

        %         [~, SortedIndex] = sort(Fitness, 'ascend');
        %         popu = popu(SortedIndex,:);


        %         mkdir traj
        G=1;
        while nFES <= FES
            %% Search machine
            algorithmDir = 'GWO';
            searchMachine = @GWO;
            %             [Pop, PopOld2, popu] = searchMachine(Fitness, PopOld, PopOld2, lu, PopSize,ProblemSize);
            [popu]=GWO_operator(popuSize,popu,Xmin,Xmax,maxIter,Alpha_pos,Beta_pos,Delta_pos,nFES);
            [popu] = BoundaryDetection(popu,lu);
            G =G +1;
            %% Boundary detection
            %                   [popu] = BoundaryDetection(popu,lu);


            %% Population evaluation
            popuFitness = [];
            for i=1:popuSize
                [popuFitness(i,:),IAEI1(i,:)] =  fun(popu(i,:),func_flag);
                nFES = nFES + 1;
                FE_best{func_flag}(time,nFES) =  popuFitness(i,:);
            end
            %                   [popuFitness] = Evaluation(popu, benchmark, f);
            %                   popuFitness = popuFitness';


            %% Update gBest
            for i = 1 : popuSize
                if popuFitness(i) < Alpha_score
                    Alpha_score = popuFitness(i);
                    Alpha_pos = popu(i, :);
                    IAEI3 = IAEI1(i,:);
                end
                if Fitness(i)>Alpha_score && Fitness(i)<Beta_score
                    Beta_score = Fitness(i); % Update beta
                    Beta_pos = popu(i,:);
                end
                if Fitness(i)>Alpha_score && Fitness(i)>Beta_score && Fitness(i)<Delta_score
                    Delta_score = Fitness(i); % Update delta
                    Delta_pos = popu(i,:);
                end
                optimal = Alpha_score;
                optimalChart = [optimalChart;optimal];
            end
            %             IAEI2 = [IAEI2;IAEI1];
            [Fitness, I] = min([Fitness, popuFitness], [], 2);


            %%
            % For minimization problems
            optimal = min(Fitness);
            record(G) = min(Fitness);

        end
        cputime(problem,time) = toc;
        fprintf('%s problem %5.0f time %5.0f |%5.0f -----> %9.16f\n', algorithmDir, problem,time,nFES,optimal);
        optimalChart = optimalChart(1:50:FES);
        resultChart = [resultChart,optimalChart];
        best_index = find(min(Fitness)==Fitness);
        best_individual{problem,time} = popu(best_index(1),:);
        outcome(problem, time) = min(Fitness);
        %         IAEI{problem,time} = IAEI2(end,:);
        IAEI{problem,time} = IAEI3;
        Record{problem,time}=record;
        time = time + 1;
    end
    analysisChart = [analysisChart;resultChart(end,:)];
    path = ['./',algorithmDir];
    sheetName = ['F',num2str(problem)];
    xlswrite(path,resultChart,[sheetName,'_D',num2str(D)]);
    mean_value=mean(outcome(problem,:),2);
    std_value = std(outcome(problem,:),0,2);
end
mean_ = mean(analysisChart,2);
std_ = std(analysisChart,0,2);
analysisChart = [analysisChart,mean_];
analysisChart = [analysisChart,std_];
xlswrite(path,analysisChart,['Statistic']);
for j =1:6
    final_data (j,1) = max(outcome(j,:));
    final_data (j,2) = min(outcome(j,:));
    final_data (j,3) = mean(outcome(j,:));
    final_data (j,4) = std(outcome(j,:));
end
save('GWO','outcome','Record','best_individual','FE_best','cputime','final_data','IAEI')
