
function outcome=ADCCDE()
    
% Reference: Ying Xia Wu, Yu Hong Wang, Sheng Xin Zhang, Li Ming Zheng, Shao Yong Zheng, 
% Adaptive dual-module cooperation and competition differential evolution, 
% Engineering Applications of Artificial Intelligence. 181,115602, (2026).

    clc;
    clear;
    rand('seed', sum(100 * clock));
  
    totaltime = 1;
    problem_size = 30;

for problem_size = problem_size
    pop_size = round(sqrt(problem_size) * log(problem_size) * 25);
    max_pop_size = pop_size;
    min_pop_size = 4;
    Max_FES = problem_size * 10000;
    val_2_reach = 10^(-8);
    lu = [-100 * ones(1, problem_size); 100 * ones(1, problem_size)];
for func = 5 %[1,3:30]
  optimum = func * 100.0;

  %% Record the best results
  outcome = []; 

  fprintf('\n-------------------------------------------------------\n')
  fprintf('Function = %d, Dimension size = %d\n', func, problem_size)
  G_Max = 0;
  pb = 0.4;
  ps = .5;
  freq_inti = 0.5;
  if problem_size == 10
      G_Max = 2163;
  end
  if problem_size == 30
      G_Max = 2745;
  end
  if problem_size == 50
      G_Max = 3022;
  end
  if problem_size == 100
      G_Max = 3401;
  end
  

  for run_id = 1 : totaltime
%   parfor run_id = 1 : totaltime    
    warning off;
    pop_size = round(sqrt(problem_size) * log(problem_size) * 25);
    %% Initialize the main population
    popold = repmat(lu(1, :), pop_size, 1) + rand(pop_size, problem_size) .* (repmat(lu(2, :) - lu(1, :), pop_size, 1));
    pop = popold; % the old population becomes the current population

    fitness = cec17_func(pop',func);
    objfitness = fitness';

    popsize = pop_size;
    arc_rate = 1;
    arc_size = round(popsize * arc_rate);
    memory_size = 5;
    memory_pos = 1;
    memory_sf = 0.5 * ones(1, memory_size);
    memory_cr = 0.8 * ones(1, memory_size);
    SEL = round(ps*pop_size);
    pop_sf = zeros(1, popsize);
    pop_cr = zeros(1, popsize);
    g_p_best_rate = 0.25;
    p_best_rate = g_p_best_rate;
    p_num = round(popsize * p_best_rate);

    archive = [];
    arc_ind_count = 0;
    archive2.NP = arc_rate * pop_size; % the maximum size of the archive
    archive2.pop = zeros(0, problem_size); % the solutions stored in te archive
    archive2.funvalues = zeros(0, 1); % the function value of the archived solutions
    [bsf_fitness, minindex] = min(objfitness);
    bsf_solution = pop(minindex, :);
    FES = popsize;
    gen = 0;
    gg=0;  %%% generation counter used For Sin
    memory_sf2 = 0.5 .* ones(memory_size, 1);
    memory_cr2 = 0.5 .* ones(memory_size, 1);

    memory_freq = freq_inti*ones(memory_size, 1);
    memory_pos2 = 1;      
    flag1 = false;
    flag2 = false;

    goodF1all = [];
    goodF2all =[];
    badF1all = [];
    badF2all = [];
    goodF1 = [];
    goodF2 = [];
    badF1 = [];
    badF2 = [];
    Dis_M = 0.5;  % initial psi for ESA
    %% main loop
    while FES < Max_FES
    gen = gen + 1;
    child = zeros(popsize, problem_size);
    [temp_fit, sorted_array] = sort(objfitness, 'ascend');
    
    std_p=std(pop);
    [~,idx_d]=sort(std_p);

    rank_d=zeros(1,problem_size);
    rank_d(idx_d)=1:problem_size;

    rank_i=zeros(1,pop_size);
    rank_i(sorted_array)=1:pop_size;
      
    Distance1 = zeros(1,pop_size); Distance2 = zeros(1,pop_size);
    
 for VN = 1:2   % VN == 1 is jSO, VN == 2 is L-SHADE_cnEpSin
   if VN == 1         
    for i = 1 : popsize
        rmz = ceil(rand * memory_size);
        if rmz == memory_size
            mu_sf = 0.9; mu_cr = 0.9;
        else
            mu_sf = memory_sf(rmz);
            mu_cr = memory_cr(rmz);
        end
        if mu_cr < 0
            pop_cr(i) = 0;
        else
            pop_cr(i) = mu_cr + 0.1 * sqrt(-2 * log(rand)) * sin(2 * pi * rand); 
            if pop_cr(i) > 1
                pop_cr(i) = 1;
            elseif pop_cr(i) < 0
                pop_cr(i) = 0;
            end
        end
        if FES < 0.25 * Max_FES && pop_cr(i) < 0.7
            pop_cr(i) = 0.7;
        end
        if FES < 0.5 * Max_FES && pop_cr(i) < 0.6
            pop_cr(i) = 0.6;
        end
        
        pop_sf(i) = mu_sf + 0.1 * tan(pi * (rand - 0.5)); 
        while pop_sf(i) <= 0
            pop_sf(i) = mu_sf + 0.1 * tan(pi * (rand - 0.5));  
        end
        
        if pop_sf(i) > 1
            pop_sf(i) = 1;
        end
        if FES < 0.6 * Max_FES && pop_sf(i) > 0.7
            pop_sf(i) = 0.7;
        end
        
        p_best_ind = sorted_array(ceil(rand * p_num));
        while FES < 0.5 * Max_FES && p_best_ind == i
            p_best_ind = sorted_array(ceil(rand * p_num));   
        end

        jF = pop_sf(i);
        if FES < 0.2 * Max_FES
            jF = jF * 0.7;
        elseif FES < 0.4 * Max_FES
            jF = jF * 0.8;
        else
            jF = jF * 1.2;
        end
        r1 = ceil(rand * popsize);
        while r1 == i
            r1 = ceil(rand * popsize);
        end
        r2 = ceil(rand * (popsize + arc_ind_count));
        while r2 == i || r2 == r1
            r2 = ceil(rand * (popsize + arc_ind_count));
        end
        random_variable = ceil(rand * problem_size);
        if r2 > popsize
            r2 = r2 - popsize;
            rrr = rand(1, problem_size);
            from_el = rrr < pop_cr(i);
            from_el(random_variable) = 1;
            child(i, :) = from_el .* (pop(i, :) + jF * (pop(p_best_ind, :) - pop(i, :)) + pop_sf(i) * (pop(r1, :) - archive(r2, :))) + (1 - from_el) .* (pop(i, :)); 
        else
            rrr = rand(1, problem_size);
            from_el = rrr < pop_cr(i);
            from_el(random_variable) = 1;
            child(i, :) = from_el .* (pop(i, :) + jF * (pop(p_best_ind, :) - pop(i, :)) + pop_sf(i) * (pop(r1, :) - pop(r2, :))) + (1 - from_el) .* (pop(i, :));
        end
        %---repair boundary
        less_min = child(i, :) < lu(1, :);  
        child(i, :) = (1 - less_min) .* child(i, :) + less_min .* ((pop(i, :) + lu(1, :)) / 2);
        more_max = child(i, :) > lu(2, :);
        child(i, :) = (1 - more_max) .* child(i, :) + more_max .* ((pop(i, :) + lu(2, :)) / 2);
        %-------------------  
    end  %end popsize ,end for  
       ui = child;
   else
            sorted_index=sorted_array;
            mem_rand_index = ceil(memory_size * rand(pop_size, 1));
            mu_sf = memory_sf2(mem_rand_index);
            mu_cr = memory_cr2(mem_rand_index);
            mu_freq = memory_freq(mem_rand_index);
            
            %% for generating crossover rate
            cr = normrnd(mu_cr, 0.1);
            term_pos = find(mu_cr == -1);
            cr(term_pos) = 0;
            cr = min(cr, 1);
            cr = max(cr, 0);
            
            %% for generating scaling factor
            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) + 0.1 * tan(pi * (rand(length(pos), 1) - 0.5));
                pos = find(sf <= 0);
            end
         
            freq = mu_freq + 0.1 * tan(pi*(rand(pop_size, 1) - 0.5));
            pos_f = find(freq <=0);
            while ~ isempty(pos_f)
                freq(pos_f) = mu_freq(pos_f) + 0.1 * tan(pi * (rand(length(pos_f), 1) - 0.5));
                pos_f = find(freq <= 0);
            end
            
            sf = min(sf, 1);
            freq = min(freq, 1);
            
            LP = 20;
            flag1 = false;
            flag2 = false;
            nfes=FES; max_nfes=Max_FES;
            if(nfes <= max_nfes/2)
                flag1 = false;
                flag2 = false;
                if (gg <= LP)
                    p1 = 0.5;
                    p2 = 0.5;
                    c=rand;
                    if(c < p1)
                        sf = 0.5.*( sin(2.*pi.*freq_inti.*gg+pi) .* ((G_Max-gg)/G_Max) + 1 ) .* ones(pop_size,problem_size);
                        flag1 = true;
                    else
                        sf = 0.5 *( sin(2*pi .* freq(:, ones(1, problem_size)) .* gg) .* (gg/G_Max) + 1 ) .* ones(pop_size,problem_size);
                        flag2 = true;
                    end
                    
                else
                    %% compute the probability as used in SaDE 
                    ns1 = size(goodF1,1);
                    ns1_sum = 0;
                    nf1_sum = 0;
                    %               for hh = 1 : size(goodF1all,2)
                    for hh = gg-LP : gg-1
                        ns1_sum = ns1_sum + goodF1all(1,hh);
                        nf1_sum = nf1_sum + badF1all(1,hh);
                    end
                    sumS1 = (ns1_sum/(ns1_sum + nf1_sum)) + 0.01;
                    
                    
                    ns2 = size(goodF2,1);
                    ns2_sum = 0;
                    nf2_sum = 0;
                    for hh = gg-LP : gg-1
                        ns2_sum = ns2_sum + goodF2all(1,hh);
                        nf2_sum = nf2_sum + badF2all(1,hh);
                    end
                    sumS2 = (ns2_sum/(ns2_sum + nf2_sum)) + 0.01;
                    
                    p1 = sumS1/(sumS1 + sumS2);
                    p2 = sumS2/(sumS2 + sumS1);
                    
                    if(p1 > p2)
                        sf = 0.5.*( sin(2.*pi.*freq_inti.*gg+pi) .* ((G_Max-gg)/G_Max) + 1 ) .* ones(pop_size,problem_size);
                        flag1 = true;
                    else
                        sf = 0.5 *( sin(2*pi .* freq(:, ones(1, problem_size)) .* gg) .* (gg/G_Max) + 1 ) .* ones(pop_size,problem_size);
                        flag2 = true;
                    end
                end
            end
            
            r0 = [1 : pop_size];
            popAll = [pop; archive2.pop];
            [r1, r2] = gnR1R2(pop_size, size(popAll, 1), r0);
            
            pNP = max(round(p_best_rate * pop_size), 2); %% choose at least two best solutions
            randindex = ceil(rand(1, pop_size) .* pNP); %% select from [1, 2, 3, ..., pNP]
            randindex = max(1, randindex); %% to avoid the problem that rand = 0 and thus ceil(rand) = 0
            pbest = pop(sorted_index(randindex), :); %% randomly choose one of the top 100p% solutions
            
            vi = pop + sf(:, ones(1, problem_size)) .* (pbest - pop + pop(r1, :) - popAll(r2, :));
            vi = boundConstraint(vi, pop, lu);
            
           %% Bin crossover according to the Eigen coordinate system
            J_= mod(floor(rand(pop_size, 1)*problem_size), problem_size) + 1;
            J = (J_-1)*pop_size + (1:pop_size)';
            crs = rand(pop_size, problem_size) < cr(:, ones(1, problem_size));
            if rand<pb
                %% coordinate ratation             
                %%%%% Choose neighbourhood region to the best individual
                best = pop(sorted_index(1), :);
                Dis = pdist2(pop,best,'euclidean'); % euclidean distance                
                %%%% Sort
                [Dis_ordered idx_ordered] = sort(Dis, 'ascend');
                Neighbour_best_pool = pop(idx_ordered(1:SEL), :); %%% including best also so start from 1
                Xsel = Neighbour_best_pool;              
                xmean = mean(Xsel);
                % covariance matrix calculation
                C =  1/(SEL-1)*(Xsel - xmean(ones(SEL,1), :))'*(Xsel - xmean(ones(SEL,1), :));
                C = triu(C) + transpose(triu(C,1)); % enforce symmetry
                [R,D] = eig(C);
                % limit condition of C to 1e20 + 1
                if max(diag(D)) > 1e20*min(diag(D))
                    tmp = max(diag(D))/1e20 - min(diag(D));
                    C = C + tmp*eye(problem_size);
                    [R, D] = eig(C);
                end
                TM = R;
                TM_=R';
                Xr = pop*TM;
                vi = vi*TM;
                %% crossover according to the Eigen coordinate system
                Ur = Xr;
                Ur(J) = vi(J);
                Ur(crs) = vi(crs);
                %%
                ui = Ur*TM_;      
            else            
                ui = pop;
                ui(J) = vi(J);
                ui(crs) = vi(crs);        
            end
            %%%%%%%%     
   end

     if VN == 1
          u1 = ui;
          for m = 1 : popsize
              Distance1(m) =  norm(u1(m,:) - pop(m,:)); 
          end
          sf1 = pop_sf;
          cr1 = pop_cr;
      elseif VN == 2
          u2 = ui;
          for m = 1 : popsize
              Distance2(m) =  norm(u2(m,:) - pop(m,:));  
          end 
          sf2 = sf;
          cr2 = cr;
     end      
 end %%%%%%%%% end VN=1：2
 
%%%%%%%%%%%%%%%%%% 
ui_shou=[];ui_tan=[];pop_sf_shou=[];pop_cr_shou=[];pop_sf_tan=[];pop_cr_tan=[];
   for i = 1 : popsize
          
         [~,Idx1] = min([Distance1(i),Distance2(i)]);
     
      if Idx1 == 1
         ui_shou(i,:) = u1(i,:); pop_sf_shou(i) = sf1(i); pop_cr_shou(i) = cr1(i);
         ui_tan(i,:) = u2(i,:); pop_sf_tan(i) = sf2(i); pop_cr_tan(i) = cr2(i);
      else
         ui_shou(i,:) = u2(i,:); pop_sf_shou(i) = sf2(i); pop_cr_shou(i) = cr2(i);
         ui_tan(i,:) = u1(i,:); pop_sf_tan(i) = sf1(i); pop_cr_tan(i) = cr1(i);
      end
   end

    afa = 1 + (nfes/max_nfes);
    ui_coll = ui_tan;
    for i = 1 : pop_size
         ddd = floor((pop_size - rank_i(i) +1) * problem_size^afa/pop_size);
         udim = rank_d <= ddd;
         ui_coll(i,udim' == 1) = ui_shou(i,udim' == 1);
    end

    if Dis_M > 0.55 + 0.1 * (nfes/max_nfes) 
        child=ui_tan; 
        pop_sf = pop_sf_tan;pop_cr = pop_cr_tan;

    elseif Dis_M > 0.45 + 0.1 * (nfes/max_nfes) && Dis_M <= 0.55+ 0.1 * (nfes/max_nfes)
        child = ui_coll;
        pop_sf = sf1;pop_cr = cr1;
    else
        child=ui_shou; 
        pop_sf = pop_sf_shou;pop_cr = pop_cr_shou;
    end   

%%%%%%%%%%%%%%%%%%%%%%%%%%  end
     child_fitness = zeros(1,popsize);
    for i = 1 : popsize
      child_fitness(i) = cec17_func(child(i,:)',func);  
    end
    
     child_fitness = child_fitness';
     FES = FES + popsize;
    
      dist = inf * ones(pop_size,1);
      for m = 1 : pop_size
          dist(m) = norm(child(m,:) - pop(m,:));
      end
      
      [~,idxx] = sort(dist);
      rank_d = zeros(pop_size,1);
      rank_d(idxx) = [1:pop_size]';
      distR = rank_d./pop_size;   
     
      dif = abs(objfitness - child_fitness);
          
        %% I == 1: the parent is better; I == 2: the offspring is better
        I = (objfitness > child_fitness);
        goodCR = pop_cr(I == 1)';
        goodF = pop_sf(I == 1)';
        goodFreq = freq(I == 1);
        dif_val = dif(I == 1);
        good_distR = distR(I == 1);

        %% chnage here also
        %% recored bad too
        badF = pop_sf(I == 0);

        %% Change Noor
        if flag1 == true
            goodF1 = goodF;
            goodF1all = [goodF1all size(goodF1,1)];

            badF1 = badF;
            badF1all = [badF1all size(badF1,1)];

            %% Add zero for other one  or add 1 to prevent the case of having NaN
            goodF2all = [goodF2all 1];
            badF2all = [badF2all 1];

        end
        if flag2 == true
            goodF2 = goodF;
            goodF2all = [goodF2all size(goodF2,1)];

            badF2 = badF;
            badF2all = [badF2all size(badF2,1)];

            %% Add zero for other one
            goodF1all = [goodF1all 1];
            badF1all = [badF1all 1];
        end
        archive2 = updateArchive(archive2, pop(I == 1, :), objfitness(I == 1));
      
    if FES < Max_FES
        [min_child_fitness, mindex] = min(child_fitness);
        if min_child_fitness < bsf_fitness
            bsf_fitness = min_child_fitness;
            bsf_solution = child(mindex(1), :);
        end
    else
        elp = Max_FES - (FES - popsize);
        elp_child_fitness = child_fitness(1:elp);
        elp_child = child(1:elp, :);
        [min_elp_child_fitness, mindex] = min(elp_child_fitness);
        if min_elp_child_fitness < bsf_fitness
            bsf_fitness = min_elp_child_fitness;
            bsf_solution = elp_child(mindex, :);
        end
        break;
    end
    %------------Add individuals to archive----------------
    equal_fit = (objfitness == child_fitness);
    pop(equal_fit, :) = child(equal_fit, :);
    
    less_fit = (child_fitness < objfitness);
    dif_fitness = abs(objfitness(less_fit) - child_fitness(less_fit));
    objfitness(less_fit) = child_fitness(less_fit);
    success_sf = pop_sf(less_fit);
    success_cr = pop_cr(less_fit);
   if arc_size > 1
    addinto_archive = pop(less_fit, :);  
    addinto_archive_num = size(addinto_archive, 1); 
    if (arc_ind_count + addinto_archive_num) <= arc_size
        archive = [archive; addinto_archive];
        arc_ind_count = arc_ind_count + addinto_archive_num;
    else
        releft_num = addinto_archive_num + arc_ind_count - arc_size;
        into_num = arc_size - arc_ind_count;
        archive = [archive; addinto_archive(1:into_num, :)];
        arc_ind_count = arc_size;
        addinto_archive(1:into_num, :) = [];
        for jj = 1:releft_num
            r11 = ceil(rand * arc_size);
            archive(r11, :) = addinto_archive(jj, :);
        end
    end
   end %end if arc_size > 1
    pop(less_fit, :) = child(less_fit, :);
    %-----------------------------------------------
    num_success_params = length(success_sf);
    if num_success_params > 0
        old_sf = memory_sf(memory_pos);
        old_cr = memory_cr(memory_pos);
        memory_sf(memory_pos) = 0;
        memory_cr(memory_pos) = 0;

        sumA = sum(dif_fitness);
        weightA = dif_fitness / sumA;  %col vector
        all_sf_product = weightA' .* success_sf .* success_sf; %succexx_sf is row vector
        tem_sf_product = weightA' .* success_sf;
        memory_sf(memory_pos) = sum(all_sf_product);
        temp_sum_sf = sum(tem_sf_product);
        all_cr_product = weightA' .* success_cr .* success_cr;
        tem_cr_product = weightA' .* success_cr;
        memory_cr(memory_pos) = sum(all_cr_product);
        temp_sum_cr = sum(tem_cr_product);
        
        memory_sf(memory_pos) = memory_sf(memory_pos) / temp_sum_sf;
        if temp_sum_cr == 0 || memory_cr(memory_pos) == -1
            memory_cr(memory_pos) = -1;
        else
            memory_cr(memory_pos) = memory_cr(memory_pos) / temp_sum_cr;
        end
        memory_sf(memory_pos) = (memory_sf(memory_pos) + old_sf) / 2;
        memory_cr(memory_pos) = (memory_cr(memory_pos) + old_cr) / 2;
        
        memory_pos = memory_pos + 1;
        if memory_pos > memory_size
            memory_pos = 1;
        end
    end
    %%
       if num_success_params > 0
                sum_dif = sum(dif_val);
                dif_val = dif_val / sum_dif;
                
                %% for updating the memory of scaling factor
                memory_sf2(memory_pos2) = (dif_val' * (goodF .^ 2)) / (dif_val' * goodF);
                
                %% for updating the memory of crossover rate
                if max(goodCR) == 0 || memory_cr2(memory_pos2)  == -1
                    memory_cr2(memory_pos2)  = -1;
                else
                    memory_cr2(memory_pos2) = (dif_val' * (goodCR .^ 2)) / (dif_val' * goodCR);
                end
                
                %% for updating the memory of freq
                if max(goodFreq) == 0 || memory_freq(memory_pos2)  == -1
                    memory_freq(memory_pos2)  = -1;
                else
                    memory_freq(memory_pos2) = (dif_val' * (goodFreq .^ 2)) / (dif_val' * goodFreq);
                end
                
               %% for updating psi
                a = 0.9;
                Dis_M = (1 - a) * Dis_M + a * (dif_val' * (good_distR .^ 2)) / (dif_val' * good_distR); % Lehmer mean
    
                memory_pos2 = memory_pos2 + 1;
                if memory_pos2 > memory_size;  memory_pos2 = 1; end
       end
    
    %-----------Calculate the next population size-----------------------
    plan_pop_size = round((min_pop_size - max_pop_size) / (Max_FES) * FES + max_pop_size);
    if popsize > plan_pop_size
        reduction_ind_num = popsize - plan_pop_size;
        if (popsize - reduction_ind_num) < min_pop_size
            reduction_ind_num = popsize - min_pop_size;
        end
        [objfitness, sindexx] = sort(objfitness, 'ascend');
        objfitness((popsize - reduction_ind_num + 1) : popsize) = [];
        pop = pop(sindexx, :);
        pop((popsize - reduction_ind_num + 1) : popsize, :) = [];
        popsize = popsize - reduction_ind_num;
        pop_size=popsize;
        SEL = round(ps*pop_size);
        arc_size = popsize * arc_rate;
        if arc_ind_count > arc_size
            arc_ind_count = arc_size;
        end
        p_best_rate = g_p_best_rate * (1 - 0.5 * FES / Max_FES);
        p_num = round(popsize * p_best_rate);
        if p_num <=1
            p_num = 2;
        end 
    end  
    end

    bsf_error_val = min(objfitness) - optimum;
    if bsf_error_val < val_2_reach
        bsf_error_val = 0;
     end

    fprintf('ADCCDE %d th run, best-so-far error value = %1.8e\n', run_id , bsf_error_val)
    outcome = [outcome bsf_error_val];   
     
  end %% end 1 run
  
  fprintf('\n')
  fprintf('ADCCDE_mean error value = %1.8e, std = %1.8e\n', mean(outcome), std(outcome))
  
end %% end 1 function run
end
end



