MOEADGGR.m

classdef MOEADGGR < ALGORITHM
% <multi/many> <real/binary/permutation>
% MOEA/D-GGR
% delta --- 1 --- The probability of selecting candidates from neighborhood
% p --- 1 --- The parameter p of GNS
% Tm --- 0.1 --- The mating neighborhood size
% Tr --- 0.1 --- The replacement neighborhood size
% nr --- 0.1 --- The maximum number of replacement
% 
%   Author: Ruihao Zheng
%   Last modified: 17/08/2022

    methods
        function main(Algorithm,Problem)
            %% Parameter setting
            [delta, p, Tm, Tr, nr] = Algorithm.ParameterSet(1, 1, 0.1, 0.1, 0.1);

            %% initialization
            % Generate the weight vectors
            [W, Problem.N] = UniformPoint(Problem.N, Problem.M);
            Tm = ceil(Problem.N * Tm);
            Tr = ceil(Problem.N * Tr);
            nr = ceil(Problem.N * nr);
            % Detect the mating and replacement neighbors of each solution
            B = pdist2(W, W);
            [~,B] = sort(B, 2);
            Bm = B(:, 1:Tm);
            Br = B(:, 1:Tr);
            % Generate random population
            Population = Problem.Initialization();
            % Initialize the reference point
            z = min(Population.objs, [], 1);
            % Dertermine the scalar function
            switch p
                case 1
                    type = 1.2;
                case inf
                    type = 2;
                otherwise
                    type = 7 + min(999, p)*0.001;
            end
            h = (prod(W,2).^(-1/Problem.M));
            % Repair the boundary weight generated by "UniformPoint"
            zero_counter = sum(W<=1e-6, 2);
            [~, I] = max(W, [], 2);
            for i = 1 : Problem.N
                W(i, I(i)) = W(i, I(i)) - zero_counter(i)*1e-6;
            end
            % Add boundary subproblems to the closest non-boundary
            % subproblem's neighborhood
            Br = mat2cell(Br, ones(1,size(Br,1)), size(Br,2));
            subp_boundary = find(zero_counter~=0)';
            for i = 1 : length(subp_boundary)
                tmp = setdiff(B(subp_boundary(i),:), subp_boundary, 'stable');
                tmp = tmp(1);
                Br{tmp} = union(Br{tmp}, subp_boundary(i));
            end
            
            %% Optimization
            while Algorithm.NotTerminated(Population)
                % For each solution
                for i = 1 : Problem.N
                    % Choose the parents
                    if rand < delta
                        P = Bm(i, randperm(Tm));
                    else
                        P = randperm(Problem.N);
                    end

                    % Generate an offspring
                    Offspring = OperatorGAhalf_2(Population(P(1:2)));
                    
                    % Update the reference point
                    z = min(z, Offspring.obj);
                    
                    % find the most suitable problem for the offspring
                    g_O = calSubpFitness(type, Offspring.obj, z, W) .* h;
                    [~, Rj] = min(g_O);
                    
                    % Update the neighbors
                    R = Br{Rj}(randperm(length(Br{Rj})));
                    g_old = calSubpFitness(type, Population(R).objs, z, W(R, :)) .* h(R);
                    g_new = g_O(R);
                    Population(R(find(g_old>=g_new, nr))) = Offspring;
                end
            end
        end
    end
end