pop_multi.m

classdef pop_multi < pop_single
% POP_MULTI         Class definition for a population to be
%                   used for multi-objective optimization 
%
% POP_MULTI is a SubClass of POP_SINGLE. The class 
% constructor works in the same way as that of POP_SINGLE,
% with the exception that an additional property is set:
%
%   pop.num_objectives      (number of objectives)
%
% All inputs and other properties are the same as for 
% POP_SINGLE -- type 'help pop_single' for more information. 
%
% The method ITERATE is now suited to optimize multi-
% objective problems. To that end, several other (hidden) 
% methods have been implemented: 
%
%   NON_DOMINATED_SORT()
%
%       a general implementation of NSGA-II's non-dominated 
%       sorting procedure. Sorts the current population
%       according to the domination level of the individuals.
%
%
%   pool = TOURNAMENT_SELECTION(pool_size, tournament_size) 
%
%       a general tournament selection procedure, that takes
%       [tournament_size] individuals randomly selected from
%       the offspring population and lets them compete with
%       the rankings and crowding distances as determining 
%       factors. The winning individual of each tournament
%       is inserted into [pool] until that [pool] contains
%       [pool_size] individuals. 
%
%   UPDATE_ALGORITHMS()
%    
%       Called from NON_DOMINATED_SORT(), updates some 
%       globaly changing values for the different algorithms. 
%       In pop_single, this is done in REPLACE_PARENTS(), but
%       as that step is not executed here, an extra method is
%       required. This updates for instance the [lbest],
%       [nbest] and [gbest] for PSO, and the temperature for
%       ASA. 
%
%
% See also pop_single, GODLIKE. 

%      Author : Rody P.S. Oldenhuis
% Affiliation : Delft University of Technology
%               Faculty of Aerospace Engineering
%               Dep. of Astrodynamics & Satellite Systems 
%     Contact : oldnhuis@dds.nl
%   Licensing/
%    (C) info : Frankly I don't care what you do with it, 
%               as long as I get some credit when you copy 
%               large portions of the code ^_^

% last edited 28/07/2009
    
    % additional properties are also public
    properties
        num_objectives     % number of objectives
        % contents of pop_data for multi-objective
        % optimization:
        %      pop_data.parent_population
        %      pop_data.offspring_population
        %      pop_data.function_values_parent
        %      pop_data.function_values_offspring
        %      pop_data.front_number
        %      pop_data.crowding_distance        
    end
    
    % public methods
    methods (Access = public)
                
        % simple constructor: create pop_single object, and 
        % just add the number of objectives
        function pop = pop_multi(varargin)
            % construct pop_single object
            pop = pop@pop_single(varargin{:});             
            % just set the amount of objectives
            pop.num_objectives = pop.options.num_objectives;              
        end % function (constructor)
        
        % perform one multi-objective iteration
        function iterate(pop)
            
            % one iteration is:
            pop.non_dominated_sort;      % non-dominated sort
            if ~strcmpi(pop.algorithm, 'PSO')
                pool = ...               % binary tournament selection (half pop.size, not for PSO)
                    pop.tournament_selection(round(pop.size/2), 2);
            else
                pool = 1:pop.size;       % the whole population for PSO
            end
            pop.create_offspring(pool);  % create new offspring
            pop.evaluate_function;       % evaluate objective function(s)            
            
            % increase number of iterations made
            pop.iterations = pop.iterations + 1; 
            
        end % function (iterate)
        
    end % public methods
    
    % protected/hidden methods
    methods (Access = protected, Hidden)
        
        % non-dominated sort, and crowding distance assignment
        function non_dominated_sort(pop)         
            
            % this function sorts the population according to non-domination.
            % At the same time, it will compute the crowding distance for every
            % individual. It will store these values for all individuals in the 
            % data structure [pop.pop_data]. 
            
            % NOTE: NON_DOMINATED_SORT is by far the most computationally costly
            % function of both POP_SINGLE and POP_MULTI. Any effort regarding
            % increasing efficiency should be focussed on this method. 
            
            % determine if this is the first call or a subsequent call, 
            % and create fitnesses and number of variables accordingly
            if ~isempty(pop.pop_data.function_values_offspring)  
                if strcmpi(pop.algorithm, 'PSO')
                    inds = [pop.pop_data.local_best_inds;
                        pop.pop_data.offspring_population; ];
                    fits = [pop.pop_data.local_best_fits;
                        pop.pop_data.function_values_offspring];
                else
                    inds = [pop.pop_data.parent_population;
                        pop.pop_data.offspring_population; ];
                    fits = [pop.pop_data.function_values_parent;
                        pop.pop_data.function_values_offspring];
                end
                N    = 2*pop.size;
            else  
                inds = pop.individuals; 
                fits = pop.fitnesses;
                N    = pop.size;
            end
                             
            % pre-calculate some stuff for crowding distances
            crowding_dists = zeros(N, 1);              % initialize            
            [sorted_fitnesses, indices] = sort(fits);  % sort fitnesses            
            crowding_dists(indices(1, :))   = inf;     % always include boundaries
            crowding_dists(indices(end, :)) = inf;     % always include boundaries
                                         
            % count the number of solutions that dominate every other solution
            guy_is_dominated_by = inf(N, 1);            
            for guy = 1:N
                
                % extract its fitnesses
                this_guy = fits(guy, :);
                
                % This is where the largest portion of the computational cost is,
                %
                %                      -> !!BY FAR!! <-
                %
                % BSXFUN() has been used to compare everything in one go, which
                % is also much faster than using a nested for. Note that strictly 
                % speaking, TWO comparisons are required, since dominance is 
                % defined as
                %
                % (xi <= yi) for all i, AND (xi < yi) for at least one value of i
                %
                % I coded the two comparisons for completeness, but in practice 
                % two equal function values are almost never encountered, so I 
                % left it commented. To use the stricter version, just uncomment 
                % the second comparison (and the addition to the [dominated] array.
                less_eq   = bsxfun(@le, fits([1:guy-1, guy+1:N], :), this_guy);
                %less      = bsxfun(@lt, fits([1:guy-1, guy+1:N], :), this_guy);
                dominated = all(less_eq, 2);% & any(less, 2);
                
                % insert the domination count in dominatrix
                guy_is_dominated_by(guy) = nnz(dominated);
                
                % compute this guy's crowding distance
                for m = 1:pop.num_objectives                    
                    % current sorting indices
                    sort_inds = indices(:, m);                    
                    % find this guy's index
                    guys_ind = find(sort_inds == guy);                    
                    % compute crowding distance
                    if isfinite(crowding_dists(guy))
                        crowding_dists(guy) = crowding_dists(guy) + ...
                            (sorted_fitnesses(guys_ind+1, m)-sorted_fitnesses(guys_ind-1, m))/...
                            (sorted_fitnesses(end,m)-sorted_fitnesses(1,m));
                    else break
                    end % if
                end % for                
            end % for
            
            % create new population
            new_pop = zeros(pop.size, pop.dimensions);
            new_fit = zeros(pop.size, pop.num_objectives);
            fronts  = zeros(pop.size, 1);
            not_filled = pop.size;
            for i = 0:max(guy_is_dominated_by)                
                % extract current front
                front = guy_is_dominated_by == i;                
                % number of entries
                entries = nnz(front);                
                % see if it still fits in the population 
                if entries <= not_filled                     
                    % if so, insert all individuals from this front in the population
                    % and keep track of their fitnesses                    
                    new_pop(pop.size-not_filled+1:pop.size-not_filled+entries, :) = inds(front, :);
                    new_fit(pop.size-not_filled+1:pop.size-not_filled+entries, :) = fits(front, :);
                    fronts (pop.size-not_filled+1:pop.size-not_filled+entries, 1) = i;
                    % adjust number of entries that have not yet been filled
                    not_filled = not_filled - entries;                  
                % if it does not fit, insert the remaining individuals based on 
                % their crowding distance
                else break 
                end % if                
            end % for
              
            % (for the first iteration, the WHOLE current population will fit 
            % in the new population. Skip that case)
            if (N ~= pop.size)                
                % sort last front encountered by crowding-distance
                front_inds = inds(front, :);       front_fits = fits(front, :);
                [sorted_front, indices] = sort(crowding_dists(front), 'descend');                
                % extract individuals & fitnesses in proper order
                sorted_inds = front_inds(indices, :); sorted_fits = front_fits(indices, :);                    
                % insert the remaining individuals in the new population 
                new_pop(pop.size-not_filled+1:end, :) = sorted_inds(1:not_filled, :);
                new_fit(pop.size-not_filled+1:end, :) = sorted_fits(1:not_filled, :);
                fronts (pop.size-not_filled+1:end, 1) = i; 
            end  

            % insert results in data structure                
            pop.pop_data.parent_population      = new_pop;            
            pop.pop_data.function_values_parent = new_fit;
            pop.pop_data.front_number           = fronts;
            pop.pop_data.crowding_distance      = crowding_dists;  
            
            % copy individuals and fitnesses to respective class properties
            pop.fitnesses   = pop.pop_data.function_values_parent;       
            pop.individuals = pop.pop_data.parent_population;  
            
            %  some algorithms need additional operations             
            pop.update_algorithms;                       
            
        end % function (non-dominated sort)
        
        % tournament selection with crowding distances and rankings 
        % as competitive factors
        function pool = tournament_selection(pop, pool_size, tournament_size)
            
            % initialize mating pool
            pool = zeros(pool_size, 1);
            
            % fill the mating pool
            for i = 1:pool_size    
                
                % select [tournament_size] individuals at random
                equal = true;
                while equal % make sure they're not equal
                    inds  = round( rand(tournament_size, 1)*(pop.size-1) + 1); % random indices
                    equal = numel(unique(inds)) ~= numel(inds);        % check if any are equal
                end
                
                % let them compete according to
                % (xj < yj) if (rank(xj) < rank(yj))
                %           or (rank(xj) == rank(yj) and distance(xj) > distance(yj)
                ranks     = pop.pop_data.front_number(inds);
                distances = pop.pop_data.crowding_distance(inds);
                best = inds(1);         
                % As default
                % TODO: this should be calculated as a ranking, so that there is
                % always a solution.
                for j = 1:tournament_size                    
                    % compare ranks
                    less_rank = ranks(j) < [ranks(1:j-1); ranks(j+1:end)];                                       
                    % rank is less than all others
                    if all(less_rank), best = inds(j); break; end                    
                    % compare distances and rank equality
                    more_dist = ranks(j) == [ranks(1:j-1); ranks(j+1:end)] &...
                                distances(j) >= [distances(1:j-1); distances(j+1:end)];                    
                    % rank is equal, but distance is less than all others
                    if any(~less_rank & more_dist), best = inds(j); break; end                        
                end % for
                
                % insert the index of the best one in the pool
                pool(i) = best;                
            
            end % for
        end % function (tournament selection)        
         
        % initialize algorithms
        function initialize_algorithms(pop)
            
            % PSO
            if strcmpi(pop.algorithm, 'PSO')
                                
                % initialize velocities
                % (velocities are 20% of [[lb]-[ub]] interval)
                pop.pop_data.velocities = (-(pop.ub-pop.lb)/2 + ...
                    rand(pop.size, pop.dimensions) .* (pop.ub-pop.lb)) / 5;
                
                % rename for clarity
                num_objectives = pop.options.num_objectives;                
                NumNeighbors   = pop.options.PSO.NumNeighbors;
                
                % initialize neighbors     
                % TODO - variable number of neighbors
                pop.pop_data.neighbors = round(rand(pop.size, NumNeighbors)*(pop.size-1)+1);
                
                % global best solution is just taken to be the first one                
                pop.pop_data.global_best_ind = pop.individuals(1, :);
                pop.pop_data.global_best_fit = inf(1, num_objectives);
                
                % initially, local best solutions are the individuals themselves
                pop.pop_data.local_best_inds = pop.individuals;
                pop.pop_data.local_best_fits = inf(pop.size, num_objectives);
                
                % set the neighbor bests 
                for i = 1:pop.size                    
                    pop.pop_data.neighbor_best_fits(i, 1:num_objectives) = inf;
                    pop.pop_data.neighbor_best_inds(i, :) = ...
                        pop.individuals(pop.pop_data.neighbors(i, 1), :);
                end
                
            % ASA
            elseif strcmpi(pop.algorithm, 'ASA')
                
                % if the initial temperature is left empty, estimate 
                % an optimal one. This is simply the largest quantity
                % in (ub - lb), divided by 4, and squared. This 
                % ensures that during the first few iterations, 
                % particles are able to spread over the entire search 
                % space; 4 = 2*2*std(randn(inf,1)).
                if isempty(pop.options.ASA.T0) && (pop.iterations == 0)
                    % only do it upon initialization
                    
                    % find the maximum value
                    sqrtT0dv4 = max(pop.ub(1,:) - pop.lb(1, :))/4;
                    
                    % set T0
                    pop.options.ASA.T0 = sqrtT0dv4^2;
                    
                end
                
                % initialize temperature
                pop.pop_data.temperature = pop.options.ASA.T0;
                
                % initialize iterations
                pop.pop_data.iters = pop.iterations;
                
            end% if
        end % function 
        
        % Update globally changing variables associated with each algorithm
        function update_algorithms(pop)
            
            
            
            
            % select which algorithm we have
            switch upper(pop.algorithm)
            
                % Genetic Algorithm
                case 'GA' % Genetic Algorithm
                    % do nothing
                    
                % Differential Evolution
                case 'DE' % Differential Evolution
                    % do nothing
                
                % Partice Swarm Optimization
                case 'PSO' % Partice Swarm Optimization
                    
                    % set of current Pareto solutions
                    Paretos = find(pop.pop_data.front_number == 0);
                    
                    % compute new local and global bests, compute new
                    % neighbors and the new best neighbors
                    for i = 1:pop.size
                        
                        % new local bests
                        if all(pop.fitnesses(i, :) <= pop.pop_data.local_best_fits(i, :) ) || ...
                           any(Paretos == i)
                            % Replace local best when it's part of the new Pareto front,
                            % OR it dominates the previous local best
                            pop.pop_data.local_best_fits(i, :) = pop.fitnesses(i, :);
                            pop.pop_data.local_best_inds(i, :) = pop.individuals(i, :);
                        end
                        
                        % New neighbors are [NumNeighbors] particles, 
                        % randomly drawn from the Pareto front                        
                        new_neighs = Paretos(round(rand(...
                            pop.options.PSO.NumNeighbors,1)*(size(Paretos,1)-1))+1);
                        pop.pop_data.neighbors(i, :) = new_neighs;
                        
                        % update the neighbor bests
                        [dummy, max_dist] = max(pop.pop_data.crowding_distance(new_neighs));
                        best_neighbor     = new_neighs(max_dist);                        
                        pop.pop_data.neighbor_best_fits(i, :) = ...
                            pop.fitnesses(best_neighbor, :);
                        pop.pop_data.neighbor_best_inds(i, :) = ...
                            pop.individuals(best_neighbor, :);
                                                
                    end % for 
                    
                    % update the global best                    
                    [maxdist, index] = max(pop.pop_data.crowding_distance(Paretos));
                    if all(pop.fitnesses(Paretos(index), :) <= pop.pop_data.global_best_fit)
                        % replace global best when it dominates the previous one
                        pop.pop_data.global_best_fit = pop.fitnesses(Paretos(index), :);
                        pop.pop_data.global_best_ind = pop.individuals(Paretos(index), :);
                    end % if
                    
                % Adaptive Simulated Annealing
                case 'ASA' % Adaptive Simulated Annealing
                    
                    % rename some stuff 
                    T       = pop.pop_data.temperature;
                    T0      = pop.options.ASA.T0;
                    cool    = pop.options.ASA.CoolingSchedule;                      
                    iters   = pop.iterations - pop.pop_data.iters;
                    
                    % essentially, only the temperature is lowered
                    pop.pop_data.temperature = max(eps, cool(T, T0, iters));
                                                            
            end % switch            
        end % method (update_algorithms)

        % overload from pop_single
        function evaluate_function(pop)

            % ugly hack because of pop_single constructor running
            % evaluations before pop_multi construction is complete
            if isempty(pop.num_objectives)
              num_objectives = pop.options.num_objectives;
            else
              num_objectives = pop.num_objectives;
            end
          
            if isempty(pop.pop_data.function_values_offspring)
              pop.pop_data.function_values_offspring = ...
                  NaN(pop.size, num_objectives);
            end

            % call super
            evaluate_function@pop_single(pop);
        end

    end % protected/hidden methods 
    
end % classdef