%%% This is a code of DE-NPS-ML algorithm (Diversity enhancement-based Differential Evolution algorithm with linearly decreasing p-Lehmer's mean and Cauchy distribution's scale parameters, Expert Systems with Applicatiuons, https://doi.org/10.1016/j.eswa.2026.132947) 
%%% The first version of DE-NPS algorithm (Zhenghao Song, Liangliang Sun, Natalja Matsveichuk, Yuri Sotskov, ShenglongJiang, Yang Yu, Differential evolution based on individual information parametersetting and diversity measurement of aggregated distribution, Swarm Evol.Comput. 92 (2025) 101793), on which the DE-NPS-ML variant is based, has been obtained from (https://github.com/shaoguaiguai).
%%% The subroutines gnR1R2.m, updateArchive.m, and boundConstraint.m were obtained from the same link as the original DE-NPS algorithm (https://github.com/shaoguaiguai).
%%% DE-NPS has been based on L-SHADE (Ryoji Tanabe and Alex Fukunaga: Improving the Search Performance of SHADE Using Linear Population Size Reduction,  Proc. IEEE Congress on Evolutionary Computation (CEC-2014), Beijing, July, 2014.)
%%% DE-NPS has also been based on JADE's source code from Dr. Q. Zhang's website (http://dces.essex.ac.uk/staff/qzhang).
%%% cec2017 problems are currently available from Prof. P.N. Suganthan's website https://github.com/P-N-Suganthan/CEC2017-BoundContrained.


clc;
clear all;
format long;
format compact;
for problem_size = 30:30 %%%AP [10 30 50]
    
    max_nfes = 10000 * problem_size;
    
    rand('seed', sum(100 * clock));
    
    %%%AP val_2_reach = 10^(-8);
    max_region = 100.0;
    min_region = -100.0;
    lu = [-100 * ones(1, problem_size); 100 * ones(1, problem_size)];
    fhd=@cec17_func;
    optimum = 100:100:3000;
    
myCluster = parcluster('local');
myCluster.NumWorkers = 13;  
    for func = 1:30
        
        totalTime = 51;
        parfor run_id = 1:totalTime
            std0 = 0.03;%%%APstd
            std1 = 0.01;%%%APstd

            p_best_rate = 0.11;
            arc_rate = 1.4;
            memory_size = 4;%%%AP ;
            pop_size = round(25 * (log(problem_size)) * sqrt(problem_size)); %%%AP
            
            max_pop_size = pop_size;
            min_pop_size = 4.0;
            

            popold = repmat(lu(1, :), pop_size, 1) + rand(pop_size, problem_size) .* (repmat(lu(2, :) - lu(1, :), pop_size, 1));
            pop = popold; 
            
            fitness = feval(fhd,pop',func);
            fitness = fitness';
            
            nfes = 0;
            [bsf_fit_var,gbestid] = min(fitness);
            bsf_solution = pop(gbestid,:);
            
            
            memory_sf = 0.5 .* ones(memory_size, 1);
            memory_cr = 0.8 .* ones(memory_size, 1);
            memory_pos = 1;
            
            archive.NP = arc_rate * pop_size; % the maximum size of the archive
            archive.pop = zeros(0, problem_size); % the solutions stored in te archive
            archive.funvalues = zeros(0, 1); % the function value of the archived solutions
            
            itter_record_name = ['L-SHADE_fid_',num2str(func),'_',num2str(problem_size),'D_',num2str(run_id),'.dat'];

            Co = zeros(problem_size,problem_size);
            counter = zeros(pop_size,1);
            

            while nfes < max_nfes
                pop = popold; 
                [temp_fit, sorted_index] = sort(fitness, 'ascend');
                
                mem_rand_index = ceil(memory_size * rand(pop_size, 1));
                mu_sf = memory_sf(mem_rand_index);
                mu_cr = memory_cr(mem_rand_index);
                

                cr = normrnd(mu_cr, 0.1);%%%AP cr = normrnd(mu_cr, 0.1);
                term_pos = find(mu_cr == -1);
                cr(term_pos) = 0;
                if nfes < 0.4 * max_nfes
                    cr = max(cr,0.6);
                    %cr = min(cr, 1);
                else
                    cr = min(cr, 1);
                    cr = max(cr, 0);
                end
                

                std = std0 + (std1-std0)*nfes/max_nfes;%%%APstd
                if nfes < 0.2*max_nfes
    
                    sf = sqrt(2).*pi^(-1/3).*(1-mu_sf.^2).*exp(-mu_sf.^2/2) + std * sin(pi * (rand(pop_size, 1) - 0.8));%%%AP sf = sqrt(2).*pi^(-1/3).*(1-mu_sf.^2).*exp(-mu_sf.^2/2) + 0.1 * sin(pi * (rand(pop_size, 1) - 0.8));
                    
                    sf = min(sf, 0.6);
                else
                    sf = mu_sf + std * tan(pi * (rand(pop_size, 1) - 0.5));%%%AP sf = mu_sf + 0.1 * tan(pi * (rand(pop_size, 1) - 0.5));
                    
                    pos = find(sf <= 0);
                    
                    while ~ isempty(pos)
                        sf(pos) = mu_sf(pos) + std * tan(pi * (rand(length(pos), 1) - 0.5));%%%AP sf(pos) = mu_sf(pos) + 0.1 * tan(pi * (rand(length(pos), 1) - 0.5));
                        pos = find(sf <= 0);
                    end
                end
                sf = min(sf, 1);
                
                r0 = [1 : pop_size];
                popAll = [pop; archive.pop];
                [r1, r2] = gnR1R2(pop_size, size(popAll, 1), r0);
                
                pNP = max(round(p_best_rate * pop_size), 2); 
                randindex = ceil(rand(1, pop_size) .* pNP); 
                randindex = max(1, randindex); 
                pbest = pop(sorted_index(randindex), :); 
                
                vi = pop + sf(:, ones(1, problem_size)) .* (pbest - pop + pop(r1, :) - popAll(r2, :));
                vi = boundConstraint(vi, pop, lu);
                
                mask = rand(pop_size, problem_size) > cr(:, ones(1, problem_size)); 
                rows = (1 : pop_size)'; cols = floor(rand(pop_size, 1) * problem_size)+1; 
                jrand = sub2ind([pop_size problem_size], rows, cols); mask(jrand) = false;
                ui = vi;
                t1 = rand;
                bsf_solution_mask = bsf_solution(ones(pop_size, 1), :);
                a2=-1+nfes*((-1)/max_nfes);
                l = (a2-1)*rand + 1;
                
                if t1 > 0.1
                    ui(mask) = pop(mask);
                else
                    dis = abs(pop(mask)-ui(mask));
                    ui(mask)  = dis*exp(1.*l).*cos(l.*2*pi) + pop(mask);
                end
                
                ui = boundConstraint(ui, pop, lu);%%%AP 

                children_fitness = feval(fhd, ui', func);
                children_fitness = children_fitness';
                

                for i = 1 : pop_size
                    nfes = nfes + 1;
                    
                    if children_fitness(i) < bsf_fit_var
                        bsf_fit_var = children_fitness(i);
                        bsf_solution = ui(i, :);
                    end
                    
                    if nfes > max_nfes; break; end
                end

                
                dif = abs(fitness - children_fitness);
                
                
                %% I == 1: the parent is better; I == 2: the offspring is better
                I = (fitness > children_fitness);
                goodCR = cr(I == 1);
                goodF = sf(I == 1);
                dif_val = dif(I == 1);
                
                archive = updateArchive(archive, popold(I == 1, :), fitness(I == 1));
                
                [fitness, I] = min([fitness, children_fitness], [], 2);
                
                popold = pop;
                popold(I == 2, :) = ui(I == 2, :);
                
                num_success_params = numel(goodCR);
                
                if num_success_params > 0
                    sum_dif = sum(dif_val);
                    dif_val = dif_val / sum_dif;
                    
                    ahr = 4-2*nfes/max_nfes;
                    memory_sf(memory_pos) = (dif_val' * (goodF .^ ahr)) / (dif_val' * (goodF.^(ahr-1)));
                    %%%AP memory_sf(memory_pos) = (dif_val' * (goodF)) / sum(dif_val');
                    
                    if max(goodCR) == 0 || memory_cr(memory_pos)  == -1
                        memory_cr(memory_pos)  = -1;
                    else
                        memory_cr(memory_pos) = (dif_val' * (goodCR .^ ahr)) / (dif_val' * (goodCR.^(ahr-1)));
                        %%%APmemory_cr(memory_pos) = (dif_val' * (goodCR)) / sum(dif_val');
                    end
                    
                    memory_pos = memory_pos + 1;
                    if memory_pos > memory_size;  memory_pos = 1; end
                end
                for i = 1:pop_size
                    if I(i) == 1
                        counter(i) = counter(i) + 1;
                    else
                        counter(i) = 0;
                    end
                end
                pop_1 = popold -repmat(mean(popold),pop_size,1);
                covM = (pop_1'*pop_1)./(size(pop_1,1)-1);
                [V,D] = eig(covM);
                for i=1:problem_size
                    if D(i,i) < 0.0001
                        Co(i,i) = Co(i,i) + 1;
                    end
                end
                for i=1:problem_size
                    if Co(i,i) > 20 * problem_size
                        for z=1:pop_size
                            if counter(z) > 5*problem_size && z~=gbestid
                                [temp_fit_chil, sorted_index_chil] = sort(fitness, 'ascend');
                                z1 = randperm(numel(sorted_index_chil),1);
                                z2 = randperm(numel(sorted_index_chil),1);
                                if sorted_index_chil(z)/pop_size > 0.5
                                    
                                    popold(z,i) = popold(z1,i) - rand * (popold(z2,i) - bsf_solution_mask(z,i));
                                else
                                    popold(z,i) = bsf_solution_mask(z,i) + rand * (popold(z1,i) - popold(z2,i));
                                end
                                popold(z,:) = boundConstraint(popold(z,:), pop, lu);%%%AP 
                                fitness(z) = feval(fhd,popold(z,:)',func);
                                counter(z) = 0;
                            end
                            
                        end
                        Co(i,i) = 0;
                    end
                end

                plan_pop_size = round((((min_pop_size - max_pop_size) / max_nfes) * nfes) + max_pop_size);
                
                if pop_size > plan_pop_size
                    reduction_ind_num = pop_size - plan_pop_size;
                    if pop_size - reduction_ind_num <  min_pop_size; reduction_ind_num = pop_size - min_pop_size;end
                    
                    pop_size = pop_size - reduction_ind_num;
                    for r = 1 : reduction_ind_num
                        [valBest indBest] = sort(fitness, 'ascend');
                        worst_ind = indBest(end);
                        popold(worst_ind,:) = [];
                        pop(worst_ind,:) = [];
                        fitness(worst_ind,:) = [];
                    end
                    
                    archive.NP = round(arc_rate * pop_size);
                    
                    if size(archive.pop, 1) > archive.NP
                        rndpos = randperm(size(archive.pop, 1));
                        rndpos = rndpos(1 : archive.NP);
                        archive.pop = archive.pop(rndpos, :);
                    end
                end
                
            end
            
            bsf_error_val = bsf_fit_var - optimum(func);
                    
    FB(func,run_id) = bsf_error_val; % AP
    LWF(func,run_id) = nfes; % AP
    XBEST(func,run_id,:) = bsf_solution; % AP
    %save wyniki FB LWF XBEST 
    func
    run_id
    bsf_error_val
    nfes

        end 
        save results FB LWF XBEST 

    end 
    save results FB LWF XBEST 
end