classdef PPS2 < ALGORITHM

    methods
        function main(Algorithm,Problem)
            %% Parameter setting
            [delta, nr] = Algorithm.ParameterSet(0.9, 2);

            [W, Problem.N] = UniformPoint(Problem.N, Problem.M);
            T = ceil(Problem.N / 10);
            Tmax = ceil(Problem.maxFE / Problem.N);

            B = pdist2(W, W, 'cosine');
            [~, B] = sort(B, 2);
            B = B(:, 1:T);
            
            Population = Problem.Initialization();
            Z = min(Population.objs, [], 1);

            max_change = 1.0;
            L = 20;
            push_stage = true;

            ideal_list = zeros(Tmax, Problem.M);
            nadir_list = zeros(Tmax, Problem.M);

            alpha = 0.95;
            % Tc = 800;
            Tc = 0.9 * ceil(Problem.maxFE/Problem.N);
            tao = 0.1;
            cp = 2;


            epsilon_0 = 0;
            epsilon_k = 0;
            NS = [];
            
            
            while Algorithm.NotTerminated(Population)
                k = ceil(Problem.FE / Problem.N); % 当前代数
                ideal_list(k, :) = min(Population.objs, [], 1);
                nadir_list(k, :) = max(Population.objs, [], 1);

                % 1. 更新epsilon
                CV = sum(max(Population.cons, 0), 2);
                if k > L
                    max_change = cal_max_change(k, L, ideal_list, nadir_list);
                    if max_change <= alpha && push_stage == true
                        push_stage = false;
                        
                        epsilon_0 = max(CV); % 最大约束违反程度
                        epsilon_k = epsilon_0;
                    else
                        rfk = sum(CV <= 0) / Problem.N; % 可行解比例
                        epsilon_k = update_epsilon(rfk, alpha, tao, k, Tc, cp, epsilon_0, epsilon_k);
                    end
                else
                    epsilon_k = 0;
                end


                % 2. 遍历每个子问题，产生子代，选择更新
                for i = 1 : Problem.N
                    % 2.1 选择父代
                    r = rand();
                    if r < delta
                        S = B(i, :);
                    else
                        S = 1:Problem.N;
                    end
                    P = S(randperm(length(S), 2));

                    % 2.2 产生子代
                    OffSpring = OperatorDE(Problem, Population(i), Population(P(1)), Population(P(2)));

                    % 2.3 选择更新
                    Z = min(Z, OffSpring.objs);

                    %  侧重调整不同目标的权重，使得往较好的方向优化 (W不变，.*和./无区别)
                    g_old = max(abs(Population(P).objs - Z) .* W(P, :), [], 2);
                    g_new = max(abs(OffSpring.obj - Z) .* W(P, :), [], 2);

                    cv_old = sum(max(Population(P).cons, 0), 2);
                    cv_new = sum(max(OffSpring.con, 0), 2);

                    if push_stage == true 
                        Population(P(find(g_new <= g_old, nr))) = OffSpring;
                    else
                        condition = (cv_old <= epsilon_k & cv_new <= epsilon_k & g_new <= g_old) | ...
                        (cv_new == cv_old & g_new <= g_old) | ...
                        (cv_new < cv_old);
                        Population(P(find(condition, nr))) = OffSpring;
                    end
                end

                % 3. 环境选择
                NS = NDSelect([NS, Population], Problem.N);
                
                if Problem.FE >= Problem.maxFE
                    Population = NS;
                end

            end
        end
    end
end

% 计算理想点/对立点最大变化率
function result = cal_max_change(k, l, ideal_list, nadir_list)
    rzk = max(abs(ideal_list(k, :) - ideal_list(k-l, :) / max(abs(ideal_list(k-l, :)), 1e-6)));
    rnk = max(abs(nadir_list(k, :) - nadir_list(k-l, :) / max(abs(nadir_list(k-l, :)), 1e-6)));
    result = max(rzk, rnk);
end

% Pull阶段，更新epsilon
function result = update_epsilon(rfk, alpha, tao, k, Tc, cp, epsilon_0, epsilon_k)
    if rfk < alpha
        result = (1 - tao) * epsilon_k;
    else
        result = epsilon_0 * (1 - k / Tc) ^ cp;
    end
end

% 环境选择
function Population = NDSelect(Population, N) 
    % 只要可行解
    isFeasible = all(Population.cons <= 0, 2);
    Population = Population(isFeasible);
    if isempty(Population)
        return;
    end
    
    [FrontNo, ~] = NDSort(Population.objs, inf);
    Next = FrontNo == 1;
    if sum(Next) > N
        % CrowDis = CrowdingDistance(Population(Next).objs);
        % [~, index] = sort(CrowDis, 'descend');
        % Next = Next(index(1:N));

        Temp = find(Next);
        Del = Truncation(Population(Next).objs, sum(Next) - N);
        Next(Temp(Del)) = false;
    end

    Population = Population(Next);
end

% 推荐，质量很高(比CrowdingDistance效果好)
function Del = Truncation(PopObj, K)
	N = size(PopObj, 1);
	Del = false(1, N);
	Distance = pdist2(PopObj, PopObj);

	% 删除拥挤的
	while sum(Del) < K
		Remain = find(~Del);
		Temp = Distance(Remain, Remain);
		Temp = sort(Temp, 2);
		[~, Rank] = sortrows(Temp);
		Del(Remain(Rank(1))) = true; 
	end
end