%****************************************************************************************************
% Author: Biao Xu and Yong Wang
% Last Edited: 23/04/2015
% Email: xubiao1189@163.com; ywang@csu.edu.cn

% Reference: Y. Wang, B. Xu, G. Sun, and S. Yang. A two-phase differential evolution for
% uniform designs % in constrained experimental domains. IEEE Transactions on Evolutionary Computation,
% in press, DOI: 10.1109/TEVC.2017.2669098.
%****************************************************************************************************

clc;
clear;
tic;

format long;
format compact;

NP = 200;

NS = 20;

sample_numbers = 100;

numb = floor(sample_numbers/(NP/NS));

for problem = 1 : 6
    
    problem
    
    rand('seed', sum(100 * clock));
    
    store = [];
    
    for time = 1 : 1
        
        time
        
        switch problem
            
            case 1
                
                lu=[ones(1, 9), 0.8 * ones(1, 5); 175 * ones(1, 9), 2.5 * ones(1, 5)];
                
            case 2
                
                lu = [78 33 27 27 27; 102 45 45 45 45];
                
            case 3
                
                lu = [-10 -10 -10 -10 -10 -10 -10; 10 10 10 10 10 10 10];
                
            case 4
                
                lu = [0 0 -0.55 -0.55; 1200 1200 0.55 0.55];
                
            case 5
                
                lu = [-10 -10 -10 -10 -10 -10 -10 -10 0;  10 10 10 10 10 10 10 10 20];
                
            case 6
                
                lu = [0 0 0 100 6.3 5.9 4.5; 1000 40 40 300 6.7 6.4 6.25];
                
            case 7
                
                lu = [-8 -8; 8 8];
                
        end
        
        switch problem
            
            case 1
                
                boxset=xlsread('Crash-Box_50k.xlsx');
                
            case 2
                
                boxset=xlsread('G04_50k.xlsx');
                
            case 3
                
                boxset=xlsread('G09_50k.xlsx');
                
            case 4
                
                boxset=xlsread('G05_50k.xlsx');
                
            case 5
                
                boxset=xlsread('G18_50k.xlsx');
                
            case 6
                
                boxset=xlsread('G21_50k.xlsx');
                
        end
        
        Normtest = normalization(boxset, lu);
        
        Results = [];
        
        n = size(lu, 2);
        
        count = 0;
        
        % Initialize the population
        p = ones(NP, 1) * lu(1, :) + rand(NP, n) .* (ones(NP, 1) * (lu(2, :) - lu(1, :)));
        
        % Evaluate the population based on the first fitness function
        fit = fitness(p, problem);
        
       %% Phase I
        while (1)
            
            k = 0;
            
            % Determine the reference point
            referPoint = (lu(2, :) - lu(1, :)) .* rand(1, n) + lu(1, :);
            
            referPoint_Norm = normalization(referPoint, lu);
            
            p_  = [];
            fit_ = [];
            
            while k <floor(NP/NS)
                
                p_Norm = normalization(p, lu);
                
                % Find the individual in the population which is nearest to
                % the reference point
                [nearestVal, nearestIndex] = min(sum((ones(size(p_Norm, 1), 1) * referPoint_Norm -  p_Norm).^2, 2));
                
                % Compute the distance from the individuals in the population
                % to the individual whose subscript is equal to "nearestIndex"
                [neighborVal, neighborIndex] = sort(sum((ones(size(p_Norm, 1), 1) *  p_Norm(nearestIndex, :) -  p_Norm).^2, 2));
                
                % Determine the subpopulation in the kth iteration
                subpop = p(neighborIndex(1 : NS), :);
                fitSubpop = fit(neighborIndex(1 : NS), :);
                
                p(neighborIndex(1 : NS), :) = [];
                fit(neighborIndex(1 : NS), :) = [];
                
                F = 0.9;
                CR = 0.9;
                
                % Implement DE for each subpopulation
                offSubpop = DE(subpop, NS, lu, n, F, CR);
                
                % Evaluate the offspring subpopulation based on the first
                % fitness function
                fitOffSubpop = fitness(offSubpop, problem);
                
                % Implement selection between the parent subpopulation and
                % the offspring subpopulation
                new_offspring_index = find(fitSubpop >= fitOffSubpop);
                subpop(new_offspring_index, :) = offSubpop(new_offspring_index, :);
                fitSubpop(new_offspring_index, :) = fitOffSubpop(new_offspring_index, :);
                
                % Select numb feasible solutions from each subpopulation
                feasiIndex = find(fitSubpop == 0);
                if  length(feasiIndex) >= numb
                    Results = [Results; subpop(feasiIndex(1 : numb), :)];
                else
                    Results = [];
                end
                
                k = k + 1;
                p_(((k - 1) * NS + 1) : k * NS, :) = subpop;
                fit_(((k - 1) * NS + 1) : k * NS, :) = fitSubpop;
                
            end
            
            p = [p; p_];
            fit = [fit; fit_];
            
            if length(Results) >= sample_numbers
                
                break;
                
            else
                
                Results=[];
                
            end
            
        end
        
       %% Phase II
        sample_set = Results(1 : sample_numbers, :);
        
        Norm_sample_set = normalization(sample_set, lu);
        
        [~, Space_D] = SpaceCompute(Norm_sample_set);
        
        % Evaluate the population based on the second fitness function 
        FF2_old = min(Space_D);
        
        FES = 1;
        
        while (1)
            
            new_pop = DE(sample_set, sample_numbers, lu, n, F, CR);
            
            fit_new_pop = fitness(new_pop, problem);
            
            feasiIndex = find(fit_new_pop == 0);
            
            if ~isempty(feasiIndex)
                
                L_temp = length(feasiIndex);
                
                temp_solution = new_pop(feasiIndex, :);
                
                for i = 1 : L_temp
                    
                    FES = FES + 1;
                    
                    old_sample_set = sample_set;
                    
                    sample_set = [sample_set; temp_solution(i, :)];
                    
                    Norm_sample_set = normalization(sample_set, lu);
                    
                    [Space_Sec, Space_D] = SpaceCompute(Norm_sample_set);
                    
                    [Order_Space_D, OSD_index] = sort(Space_D);
                    
                    [Min_Spa, MSV_index] = min(Space_Sec(OSD_index(1 : 2)));
                    
                    sample_set(OSD_index(MSV_index), :) = [];
                    
                    Norm_sample_set = normalization(sample_set, lu);
                    
                    [~, Space_D] = SpaceCompute(Norm_sample_set);
                    
                    % Evaluate the offspring population based on the second fitness function 
                    FF2_new = min(Space_D);
                    
                    if FF2_new > FF2_old
                        
                        FF2_old = FF2_new;
                        
                        count = 0;
                        
                    else
                        
                        sample_set = old_sample_set;
                        
                        count = count + 1;
                        
                    end
                    
                end
                
                if count >= 100
                    
                    break;
                    
                end
                
            end
            
        end
        
       %% Performance Criterion
        
        NormD = normalization(sample_set, lu);
        
        [M, ~] = size(Normtest);
        
        dd = zeros(M, 1);
        
        for j = 1 : M
            
            Euler_D = Euler_Distance(Normtest(j, :), NormD);
            
            dd(j) = min(Euler_D);
            
        end
        
        md = max(dd);
        
        store = [store; md]
        
        
    end
    
end

toc;
