pop_single.m

classdef pop_single < handle
% =insert documentation here=

%      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 ^_^
    
    % all properties are public
    properties 
        algorithm          % type of optimization algorithm used 
        funfcn             % objective function(s)
        individuals        % members of the population
        fitnesses          % corresponding fitnesses        
        size               % population size
        lb                 % lower bounds
        ub                 % upper bounds
        orig_size          % original size of the input
        dimensions         % dimensions
        funevals   = 0;    % number of function evaluations made        
        iterations = 0;    % iterations so far performed
        options            % options structure (see function [set_options] for info)            
        pop_data           % structure to store intermediate data                           
        % contents for single-objective optimization:
        %      pop_data.parent_population
        %      pop_data.offspring_population
        %      pop_data.function_values_parent
        %      pop_data.function_values_offspring
    end
        
    % public methods
    methods (Access = public)
        
        % constructor
        function pop = pop_single(varargin)

            % TODO: too much stuff is done in constructor. This messes stuff
            % up for subclasses (e.g., pop_multi). Divide constructor?
          
            % default check
            error(nargchk(2, 7, nargin));
            
            % input is ( new [pop_data] structure, previous [population] object, options )
            % (subsequent call from GODLIKE)
            % = = = = = = = = = = = = = = = = = = = = = = = = = = 
            
            if (nargin == 3) 
                
                % assign new pop_data structure
                pop.pop_data = varargin{1};
                
                % simply copy previous object
                pop.funfcn     = varargin{2}.funfcn;       pop.iterations = varargin{2}.iterations;
                pop.algorithm  = varargin{2}.algorithm;    pop.lb         = varargin{2}.lb;
                pop.funevals   = varargin{2}.funevals;     pop.ub         = varargin{2}.ub;
                pop.dimensions = varargin{2}.dimensions;   pop.orig_size  = varargin{2}.orig_size;  
                
                % copy individuals and fitnesses
                pop.individuals = pop.pop_data.parent_population;
                pop.fitnesses   = pop.pop_data.function_values_parent;
                
                % size and options might have changed
                pop.size = size(pop.individuals, 1);%#ok
                pop.options = varargin{3}; 
                % replicate [ub] and [lb]
                pop.lb = repmat(pop.lb(1, :), pop.size, 1);    
                pop.ub = repmat(pop.ub(1, :), pop.size, 1);
                
                % Some algorithms need some lengthier initializing
                pop.initialize_algorithms;
                
                % return
                return
            end
            
            % input is ( funfcn, popsize, lb, ub, dimensions, options )
            % (initialization call from GODLIKE)
            % = = = = = = = = = = = = = = = = = = = = = = = = = = 
             
            % parse input
            % · · · · · · · · · · · · · · · · · · · · · · 
                        
            % assign input
            pop.funfcn  = varargin{1};   pop.ub         = varargin{4};
            pop.size    = varargin{2};   pop.orig_size  = varargin{5};
            pop.lb      = varargin{3};   pop.dimensions = varargin{6};
            pop.options = varargin{7};
            
            % cast funfcn to cell if necessary
            if ~iscell(pop.funfcn), pop.funfcn = {pop.funfcn}; end
                        
            % replicate [lb] and [ub] to facilitate things a bit
            % (and speed it up some more)
            pop.lb = repmat(pop.lb, pop.size, 1);   pop.ub = repmat(pop.ub, pop.size, 1);
            
            % set optimization algorithm
            pop.algorithm = pop.options.algorithm;
                             
            % Initialize population             
            % · · · · · · · · · · · · · · · · · · · · · · 
            
            % initialize population 
            pop.individuals = pop.lb + rand(pop.size, pop.dimensions) .* (pop.ub-pop.lb);
                                      
            % insert copy into info structure
            pop.pop_data.parent_population = pop.individuals;
                        
            % temporarily copy parents to offspring positions 
            pop.pop_data.function_values_offspring = [];
            pop.pop_data.offspring_population      = pop.individuals;
            
            % evaluate function for initial population (parents only)             
            pop.evaluate_function;   
                        
            % copy function values into fitnesses properties
            pop.fitnesses = pop.pop_data.function_values_offspring;
            pop.pop_data.function_values_parent = pop.fitnesses;
                        
            % delete entry again
            pop.pop_data.function_values_offspring = [];
            
            % some algorithms need some lengthier initializing
            pop.initialize_algorithms;
                                  
        end % function (constructor)  
        
        % single iteration
        function iterate(pop, times, FE)
            % [times] and [FE] are only used for the MultiStart algorithm
                        
            % select proper candiadates 
            if strcmpi(pop.algorithm, 'GA')
                pool = ...               % binary tournament selection for GA
                pop.tournament_selection(pop.size, 2); 
            else
                pool = 1:pop.size;       % whole population otherwise
            end 
            
            % create offspring
            if nargin == 1
                pop.create_offspring(pool);
            else
                pop.create_offspring(pool, times, FE);  
            end
            
            % if the algorithm is MS, this is the only step
            if strcmpi(pop.algorithm, 'MS')
                % adjust iterations
                pop.iterations = pop.iterations + times;
                % then return                
                return            
            end
            
            % carefully evaluate objective function(s)
            try
                pop.evaluate_function;   
            catch userFcn_ME                
                pop_ME = MException('pop_single:function_doesnt_evaluate',...
                    'GODLIKE cannot continue: failure during function evaluation.');
                userFcn_ME = addCause(userFcn_ME, pop_ME);
                rethrow(userFcn_ME);
            end
            
            % replace the parents
            pop.replace_parents;         
                        
            % increase number of iterations made
            pop.iterations = pop.iterations + 1;
                        
        end % function (single iteration)              
                
    end % methods
    
    % % protected/hidden methods
    methods (Access = protected, Hidden)
        
        % tournament selection (only for GA)
        function pool = tournament_selection(pop, pool_size, tournament_size)
            
            % initialize mating pool
            pool = zeros(pool_size, 1);
            
            % total number of competitors
            rnd_inds = zeros(pool_size*tournament_size,1);
            
            % create random indices outside the loop (faster)
            for i = 1:floor(pool_size*tournament_size/pop.size)
                offset = pop.size*(i-1);
                [dummy, rnds] = sort(rand(pop.size,1));
                rnd_inds(offset+1:min(end,offset+pop.size), :) = rnds(1:min(end,nnz(~rnd_inds)));
            end
            
            % fill the mating pool
            for i = 1:pool_size    
                
                % select [tournament_size] individuals at random
                inds = rnd_inds(1:tournament_size);
                rnd_inds = rnd_inds(tournament_size+1:end);
                
                % let them compete according to
                % (xj < yj) if fit(xj) < fit(yj)
                [best, ind] = min(pop.fitnesses(inds));                
                
                % insert the index of the best one in the pool
                pool(i) = inds(ind);                
            
            end % for
            
        end % function (tournament selection)        
        
        % generate new generation
        function create_offspring(pop, pool, times, FE)
            
            % get the size of the pool
            pool_size = length(pool);  
            
            % rename some stuff 
            parent_pop = pop.individuals(pool, :);
            parent_fit = pop.fitnesses(pool, :);
            
            % initialize
            newpop = zeros(pop.size, pop.dimensions);           % empty new population            
            newfit = NaN(pop.size, pop.options.num_objectives); % placeholder for the sites to 
                                                                % evaluate the function                        
            % determine which algorithm to use
            type = upper(pop.algorithm);
            
            % generate offspring with selected algorithm
            switch type
                
                % Multistart
                case 'MS'  % Multistart
                                                                
                    % number of function evaluations still allowed
                    allowed_FE = pop.options.MaxFunEvals - FE;
                    
                    % set relaxed options
                    options = optimset(...
                        'MaxIter', times, ...
                        'display', 'off', ...
                        'TolFun' , 10*pop.options.TolFun,...
                        'TolX'   , 10*pop.options.TolX);
                    
                    % reinitialize newpop
                    newpop = parent_pop;
                                        
                    % loop through the population
                    for i = 1:pop.size                        
                        % reset options 
                        options = optimset(options, 'MaxFunEvals', allowed_FE);                        
                        % optimize this individual
                        [newpop(i, :), newfit(i, :), ef, output] =...
                            fminsearch(pop.funfcn{1}, parent_pop(i, :), options);
                        % update function evaluations
                        pop.funevals = pop.funevals + output.funcCount;
                        % number of function evaluations still allowed
                        allowed_FE = allowed_FE - pop.funevals;                        
                        % if the maximum has been exceeded, exit
                        if allowed_FE < 0
                            % first insert result into pop
                            pop.pop_data.offspring_population      = newpop;
                            pop.pop_data.function_values_offspring = newfit;
                            % replace the parents 
                            pop.replace_parents;
                            % then exit                            
                            break                        
                        end                        
                    end
                                    
                % Differential Evolution 
                case 'DE'  % Differential Evolution 
                                        
                    % I love DE for its elegance and simplicity, and 
                    % yet powerful optimization qualities
                    
                    % rename some stuff   
                    Flb = pop.options.DE.Flb;
                    Fub = pop.options.DE.Fub;
                    crossconst = pop.options.DE.CrossConst; 
                    
                    % Neoteric Differential Evolution                    
                    for i = 1:pop.size   
                        % random indices
                        base = round(rand*(pool_size-1))+1; % RANDI is slower
                        d1   = round(rand*(pool_size-1))+1;
                        d2   = round(rand*(pool_size-1))+1;
                        % d2 may not be equal to d1
                        while (d1 == d2), d2 = round(rand*(pool_size-1))+1; end
                        % DE operator
                        if rand < crossconst || round(rand*(pool_size-1))+1 == i;
                            % DE operator when rnd < Cr
                            F = rand*(Fub-Flb) + Flb;
                            newpop(i, :) = parent_pop(base,:) +  ...
                                F*(parent_pop(d1,:) - parent_pop(d2,:));
                        else
                            % insert random parent otherwise
                            rnd_ind      = round(rand*(pool_size-1))+1;
                            newpop(i, :) = parent_pop(rnd_ind, :);
                            newfit(i, :) = parent_fit(rnd_ind, :);
                        end                        
                    end % for
                                    
                % Particle Swarm Optimization 
                case 'PSO' % Particle Swarm Optimization 
                    
                     % first, make sure the pool is large enough
                    temp_pop = zeros(pop.size, pop.dimensions);
                    if numel(parent_pop) ~= pop.size*pop.dimensions                        
                        % insert all values from the pool
                        temp_pop(1:size(parent_pop,1), :) = parent_pop; %#ok                        
                        % insert random members from [parent_pop]
                        for i = size(parent_pop,1)+1:pop.size %#ok                            
                            randinds = round(rand*(size(parent_pop,1)-1)+1); %#ok                            
                            temp_pop(i, :) = parent_pop(randinds, :);
                        end                        
                        % equate the two
                        parent_pop = temp_pop;
                    end     

                    % Creating offspring with PSO is pretty simple:
                    newpop = parent_pop + pop.pop_data.velocities;
                                        
                    % Since the velocity can not be zero, each individual
                    % changes during the creation of offspring, so the
                    % function values of ALL new individuals have to be 
                    % re-calculated.
                                        
                % Genetic Algorithm
                case 'GA'  % Genetic Algorithm
                                        
                    % Generating offspring in GA requires quite many 
                    % operations. Compared to DE or PSO, it's really 
                    % quite messy:
                    
                    % rename some stuff
                    Coding     = pop.options.GA.Coding;
                    MutProb    = pop.options.GA.MutationProb;
                    CrossProb  = pop.options.GA.CrossProb;
                    NumBits    = pop.options.GA.NumBits;                    
                    if strcmpi(Coding, 'Binary')
                        Binary = true;  Real = false;
                    else
                        Binary = false; Real = true;
                    end
                                        
                    % save signs
                    signs = sign(parent_pop);
                    
                    % initialize some arrays that keep track of the signs
                    child_signs = zeros(2, pop.dimensions);
                    new_signs   = zeros(pop.size, pop.dimensions);
                    
                    % convert to binary 
                    if Binary                               
                        % find correct multiplier
                        multiplier = 1; 
                        temp_pop   = round(parent_pop);           % initialize
                        parent_pop = abs(parent_pop);             % take absolute value
                        temp_pop   = abs(temp_pop);               % take absolute value
                        while (max(temp_pop(:)) <= 2^(NumBits)) && ~all(temp_pop(:) == 0)
                            multiplier = multiplier*10;           % adjust multiplier
                            temp_pop   = round(parent_pop*multiplier);% convert to integers
                        end
                        multiplier = multiplier/10;               % correct multiplier
                        temp_pop   = round(parent_pop*multiplier);% convert to integers
                        
                        % see if selected number of bits cannot represent population
                        if (multiplier == 0.1) && ~all(temp_pop(:) == 0)
                            error('pop_single:numbits_insufficient', ...
                                  ['Maximum value in population can not be represented by the\n',...
                                  'selected number of bits. Increase ''NumBits'' option, or\n',...
                                  'rescale your problem.']);
                        end                                                                       
                        % convert each column separately
                        bit_representation = false(pool_size, pop.dimensions*NumBits);
                        % NOTE: MATLAB's DEC2BIN() should be avoided, as it would 
                        % be called in a loop and is not a builtin. Moreover, the
                        % output of DEC2BIN() is a string, which is pretty 
                        % inconvenient for the mutation operator. Therefore, do 
                        % the conversion manually
                        for i = 1:pop.dimensions 
                            % convert this column to bits
                            bits = temp_pop(:, i);
                            bits = bits*pow2(1-NumBits:0);
                            bits = floor(bits);            % oddly enough, this is much faster 
                            bits = bits - fix(bits/2)*2 ;  % than using REM(bits,2)...
                            bits = logical(bits);
                            % append in output matrix           
                            bit_representation(:, NumBits*(i-1)+1:NumBits*(i-1)+NumBits) = bits;                            
                        end                        
                        % equate the two
                        parent_pop = bit_representation;
                        % redefine newpop
                        newpop = false(pop.size, pop.dimensions*NumBits);                        
                        % initialize children
                        children = false(2, size(parent_pop,2)); %#ok
                        % and define convenient conversion array
                        % (starts mattering for large population sizes)
                        convert_to_dec = repmat(2.^(NumBits-1:-1:0), pop.size, 1);
                                             
                    % convert to array of strings in case of real-representation
                    elseif Real
                        % do everything in one go with INT2STR()    
                        % (avoid NUM2STR, as its horrifically slowin a loop)
                        % (note that the signs are still included in the array)
                        real_representation = int2str(abs(parent_pop)*1e18);
                        % convert to array of doubles
                        real_representation = real_representation - '0';                       
                        % equate the two
                        parent_pop = real_representation;
                        % initialize children
                        children = zeros(2, size(parent_pop, 2)); %#ok
                        % redefine newpop
                        newpop = zeros(pop.size, size(parent_pop, 2)); %#ok
                    end
                                                            
                    % perform crossover
                    for i = 1:2:pop.size-1
                        
                        % select two parents at random
                        parent_inds  = round(rand(2,1)*(pool_size-1)+1);
                        parents      = parent_pop(parent_inds, :);
                        if Binary, parent_signs = signs(parent_inds, :); end
                                                
                        % crossover if a random number is less than [CrossProb]
                        % otherwise, just insert the two parents into the new 
                        % population 
                        if (rand < CrossProb)
                                                        
                            % select random crossover point
                            crosspoint = round(rand*(size(parents,2)-1))+1; %#ok
                                                                                                                
                            % perform crossover
                            children(1, :) = [parents(1,1:crosspoint),parents(2,crosspoint+1:end)];
                            children(2, :) = [parents(2,1:crosspoint),parents(1,crosspoint+1:end)];
                                                         
                            % also keep track of the signs 
                            if Binary
                                index = ceil(crosspoint/NumBits);    
                                child_signs(1, :) = [parent_signs(1,1:index),...
                                    parent_signs(2,index+1:end)];
                                child_signs(2, :) = [parent_signs(2,1:index),...
                                    parent_signs(1,index+1:end)];
                            end
                            
                            % insert children 
                            newpop(i:i+1, :)    = children;
                            if Binary, new_signs(i:i+1, :) = child_signs; end
                                                        
                        else
                            newpop(i:i+1, :)    = parents;   
                            newfit(i:i+1, :)    = parent_fit(parent_inds, :);
                            if Binary, new_signs(i:i+1, :) = parent_signs; end
                        end % if   
                    end % for
                    
                    % if the population size is an uneven number, the last entry 
                    % is still open. Just stick a random parent there
                    if mod(pop.size,2)
                        index = round(rand*(pool_size-1)+1);
                        newpop(end, :)   = parent_pop(index, :);
                        newfit(end, :)   = parent_fit(index, :);
                        if Binary, new_signs(end,:) = signs(index,:); end
                    end                    
                    
                    % mutation operator
                    mutate = rand(pop.size, size(parent_pop,2)) < MutProb; %#ok
                    % If any individual mutates, the function has to be re-evaluated
                    newfit(sum(mutate,2)>0,:) = NaN;                     
                    % Binary coding - simply flip bits
                    if Binary, newpop(mutate) = ~newpop(mutate); end
                    % Real coding - select a new number from [0,9]
                    if Real
                        % don't mutate spaces
                        space_inds = newpop(mutate) == (' '-'0');
                        mutate(mutate) = ~space_inds;
                        % don't mutate signs
                        sign_inds = newpop(mutate) == ('-'-'0');
                        mutate(mutate) = ~sign_inds;
                        % random new digits                        
                        rnd_inds = round(rand(nnz(mutate),1)*9);
                        % convert to strings
                        newpop(mutate) = rnd_inds;                     
                    end
                                        
                    % convert back to real numbers                  
                    if Binary                        
                        % initialize
                        temp_pop = zeros(pop.size, pop.dimensions);                        
                        % convert columnwise again
                        for i = 1:pop.dimensions                            
                            % convert column to decimal representation   
                            temp_pop(:, i) = sum(convert_to_dec.*newpop(:, 1:NumBits), 2);                           
                            % delete entries
                            newpop(:, 1:NumBits) = [];                             
                        end                        
                        % divide by multiplier, and re-assign signs
                        temp_pop = temp_pop/multiplier.*new_signs;                        
                        % assign newpop
                        newpop = temp_pop;                        
                    elseif Real
                        % initialize
                        temp_pop = zeros(pop.size, pop.dimensions);
                        % assign space character
                        space = ' '-'0';   
                        % append one "space" to the end 
                        newpop(:, end+1) = space;
                        % then convert back to double, column per column
                        for i = 1:pop.dimensions
                            % trim leading "spaces"
                            while all(newpop(:,1)==space), newpop = newpop(:, 2:end); end
                            % first find one that does not begin with a "space"                           
                            non_space = find(newpop(:,1) ~= space, 1);
                            % find indices for the next "space"
                            space_ind = find(newpop(non_space,:) == space, 1);
                            space_ind = space_ind-1;
                            % use power trick forthe conversion
                            powers = 10.^(space_ind-1:-1:0);
                            powers = powers(ones(pop.size,1),:);
                            % remove residual spaces
                            ttemp_pop = newpop(:,1:space_ind);
                            ttemp_pop(ttemp_pop == space) = 0;
                            % insert in final array
                            temp_pop(:, i) = sum(ttemp_pop.*powers,2)/1e18;
                            % adjust newpop
                            newpop = newpop(:, space_ind+1:end);
                        end
                        % assign newpop
                        % CG: removed 'signs.*' because it's not needed:
                        % Real representation preserves signs? 
                        newpop = temp_pop; 
                    end
                                     
                % Adpative Simulated Annealing
                case 'ASA' % Adpative Simulated Annealing
                    
                    % Generating the new population is very straightforward 
                    % in ASA. It only requires changing ONE of the decision 
                    % variables per individual, which can easily be vectorized.
                    
                    % first, make sure the pool is large enough
                    temp_pop = zeros(pop.size, pop.dimensions);
                    if numel(parent_pop) ~= pop.size*pop.dimensions
                        
                        % insert all values from the pool
                        temp_pop(1:size(parent_pop,1), :) = parent_pop; %#ok
                        
                        % insert random members from [parent_pop]
                        for i = size(parent_pop,1)+1:pop.size %#ok                            
                            randinds = round(rand*(size(parent_pop,1)-1)+1); %#ok                            
                            temp_pop(i, :) = parent_pop(randinds, :);
                        end
                        
                        % equate the two
                        parent_pop = temp_pop;
                    end              
                    
                    % Use Bolzmann distribution to create new points              
                    rands  = sqrt(pop.pop_data.temperature)*randn(pop.size, pop.dimensions);                                      
                    newpop = parent_pop + rands;
                    % at least one dimension always changes, so ALL function 
                    % values have to be recomputed
                    
            end % switch
                        
            % check constraints and boundaries after 
            % offspring generation
            [newpop, newfit] = pop.honor_bounds(newpop, newfit);
                         
            % insert result into pop
            pop.pop_data.offspring_population      = newpop;
            pop.pop_data.function_values_offspring = newfit;
            
        end % function (create offspring)
        
        % selectively replace the parent population with members 
        % from the offspring population (single-objective optimization)
        function replace_parents(pop)
            
            % rename for clarity
            new_fits = pop.pop_data.function_values_offspring;
            new_inds = pop.pop_data.offspring_population;
            
            % operation depends on the algorithm again
            switch upper(pop.algorithm)
                
                % Multistart algorithm
                case 'MS'
                    % MS just replaces everything
                    pop.pop_data.parent_population = new_inds;
                    pop.pop_data.function_values_parent = new_fits;
                
                % Differential Evolution 
                case 'DE' % Differential Evolution                     
                    % DE and GA both use simple greedy replacement 
                    better_inds = new_fits < pop.fitnesses;
                    pop.pop_data.parent_population(better_inds, :) = new_inds(better_inds, :);
                    pop.pop_data.function_values_parent(better_inds, :) = new_fits(better_inds, :);                
                % Genetic Algorithm
                case 'GA'  % Genetic Algorithm
                    % DE and GA both use simple greedy replacement 
                    better_inds = new_fits < pop.fitnesses;
                    pop.pop_data.parent_population(better_inds, :) = new_inds(better_inds, :);
                    pop.pop_data.function_values_parent(better_inds, :) = new_fits(better_inds, :);
                    
                % Particle Swarm Optimization
                case 'PSO' % Particle Swarm Optimization
                                        
                    % PSO simply replaces all parents
                    pop.pop_data.parent_population = new_inds;
                    pop.pop_data.function_values_parent = new_fits;
                    
                    % update the neighbor bests
                    % (this implementation is fast, not intuitive)
                    % add one NaN to the new_fits array
                    new_fits  = [new_fits; NaN];
                    % copy neighbors
                    neighbors = pop.pop_data.neighbors;
                    % let those that are zero refer to the NaN entry
                    neighbors(neighbors == 0) = size(new_fits,1); %#ok      
                    % find the best ones
                    [neighbor_best, ind] = min(new_fits(neighbors),[],2);                    
                    % find those that are better
                    better_neighbors = neighbor_best <  pop.pop_data.neighbor_best_fits;
                    % no better ones might be found
                    if any(better_neighbors)
                        % insert function values
                        pop.pop_data.neighbor_best_fits(better_neighbors) = ...
                            neighbor_best(better_neighbors);
                        % insert individuals
                        for i = (find(better_neighbors)).'
                            pop.pop_data.neighbor_best_inds(i, :) = ...
                                new_inds(neighbors(i, ind(i)), :);
                        end
                    end
                    % chop off additional NaN-entry again
                    new_fits = new_fits(1:end-1);
                    
                    % update the local bests
                    new_locals = new_fits < pop.pop_data.local_best_fits;
                    pop.pop_data.local_best_fits(new_locals, 1) = new_fits(new_locals, 1);
                    pop.pop_data.local_best_inds(new_locals, :) = new_inds(new_locals, :);                    
                    
                    % update the global best
                    if (min(new_fits) < pop.pop_data.global_best_fit)                        
                        [pop.pop_data.global_best_fit, ind] = min(new_fits);                        
                        pop.pop_data.global_best_ind = new_inds(ind, :);
                    end
                    
                    % create random matrices
                    r1  = rand(pop.size, 1);  r1 = r1(:, ones(1,pop.dimensions)); 
                    r2  = rand(pop.size, 1);  r2 = r2(:, ones(1,pop.dimensions)); 
                    r3  = rand(pop.size, 1);  r3 = r3(:, ones(1,pop.dimensions)); 

                    % update velocities                     
                    pop.pop_data.velocities = ...
                        pop.options.PSO.omega  *pop.pop_data.velocities + ...
                        pop.options.PSO.eta1   *r1.*(pop.pop_data.neighbor_best_inds - new_inds)+...
                        pop.options.PSO.eta2   *r2.*...
                                          bsxfun(@minus, pop.pop_data.global_best_ind, new_inds)+...
                        pop.options.PSO.eta3   *r3.*(pop.pop_data.local_best_inds - new_inds);
                    
                   % check the bounds
                   pop.honor_bounds([]);
                                                                                  
                % Adaptive Simulated Annealing
                case 'ASA' % Adaptive Simulated Annealing
                    
                    % rename some stuff                    
                    T       = pop.pop_data.temperature;
                    T0      = pop.options.ASA.T0;
                    nrg     = pop.pop_data.function_values_offspring;
                    prevnrg = pop.pop_data.function_values_parent;    
                    cool    = pop.options.ASA.CoolingSchedule;                      
                    iters   = pop.iterations - pop.pop_data.iters;
                                        
                    % reject or accept the new population, according to 
                    % the probabilistic rule
                    nrgdiff  = (prevnrg - nrg);                % energy difference
                    nrgdiff  = nrgdiff / max(abs(nrgdiff(:))); % rescale the differences
                    ind      = nrgdiff > 0;                    % always accept better ones
                    probind  = ~ind & rand(pop.size, 1) < exp( nrgdiff/T );
                                                               % accept worse ones based on 
                                                               % probabalistic rule     
                    swapinds = ind|probind;                    % indices to be swapped
                                                              
                    % apply cooling schedule
                    pop.pop_data.temperature = max(eps,cool(T, T0, iters));
                                                              
                    % replace the individuals
                    pop.pop_data.parent_population(swapinds, :) = new_inds(swapinds, :);                    
                    
                    % also replace the function values
                    pop.pop_data.function_values_parent(swapinds, :) = new_fits(swapinds, :);                    
                    
            end % switch     
            
            % copy individuals and fitnesses to respective properties            
            pop.fitnesses   = pop.pop_data.function_values_parent;       
            pop.individuals = pop.pop_data.parent_population;   
            
        end % function
        
        % evaluate the objective function(s) correctly
        function evaluate_function(pop)
                   
            % multi-objective optimization overloads this function for
            % some initialization, but returns here later.
            
            % find evaluation sites
            if isempty(pop.pop_data.function_values_offspring)
                sites = 1:pop.size; % only in pop-initialization
                fvs = zeros(length(sites), 1);  % Single-objective only
            else
                % TODO: check if this works correctly; set as logicals
                % here, but used as index below
                sites = ~isfinite(pop.pop_data.function_values_offspring(:, 1));
                fvs = zeros(nnz(sites), ...
                            size(pop.pop_data.function_values_offspring, 2));
            end
            num_pop = nnz(sites);

            pop_inputs = ...
              pop.pop_data.offspring_population(sites, :);

            % evaluate all functions for each population member in
            % parallel or in serial
            if pop.options.UseParallel
              parfor pop_num=1:num_pop
                fvs(pop_num, :) = ...
                    evaluate_one_function(pop, pop_inputs(pop_num, :));
              end
            else
              for pop_num=1:num_pop
                fvs(pop_num, :) = ...
                    evaluate_one_function(pop, pop_inputs(pop_num, :));
              end
            end
            pop.pop_data.function_values_offspring(sites, :) = fvs;
                          
            % update number of function evaluations
            % (count each function call as one even though it returns
            % multiple objectives)
            pop.funevals = pop.funevals + num_pop*numel(pop.funfcn); 

        end % function 

        % Evaluate a single iterations of function(s), so that this can be run in
        % parallel.  
        function pop_output = evaluate_one_function(pop, pop_input)
        
          % TODO: move this to pop_multi and keep a simpler one here. Test
          % with demo.
            
          if pop.options.obj_columns
            pop_output = ...
                feval(pop.funfcn{1}, pop_input);
          else
              for ii = 1:numel(pop.funfcn)
                pop_output(1, ii) = ...
                    feval(pop.funfcn{ii}, pop_input);
              end
          end
          
        end % function
          
        % check boundaries
        function [newpop, newfit] = honor_bounds(pop, newpop, newfit)
                                    
            % find violation sites
            outsiders1 = false; outsiders2 = false;
            if ~isempty(newpop)
                outsiders1 = newpop < pop.lb;
                outsiders2 = newpop > pop.ub;
            end
            
            % PSO requires more elaborate check
            if strcmpi(pop.algorithm, 'PSO')
                
                % rename for clarity
                velocity = pop.pop_data.velocities; % extract velocities
                velUb = (pop.ub - pop.lb)/5;        % upper bounds on velocity
                velLb = (pop.ub - pop.lb)/1e50;     % lower bounds on velocity
                          
                % bounce against bounds
                if any(outsiders1(:) | outsiders2(:))
                    newpop(outsiders1)   = newpop(outsiders1) - velocity(outsiders1);
                    newpop(outsiders2)   = newpop(outsiders2) - velocity(outsiders2);
                    velocity(outsiders1) = -velocity(outsiders1);
                    velocity(outsiders2) = -velocity(outsiders2);                    
                end
                
                % limit velocity
                Velsign    = sign(velocity);
                outsiders1 = abs(velocity) > abs(velUb);
                outsiders2 = abs(velocity) < abs(velLb);
                if any(outsiders1(:)) || any(outsiders2(:))
                    velocity(outsiders1) = Velsign(outsiders1).*velUb(outsiders1);
                    velocity(outsiders2) = Velsign(outsiders2).*velLb(outsiders2);
                end
                
                % re-insert velocity
                pop.pop_data.velocities = velocity;                
                
            % boundary violations in all other algorithms 
            % are simply reinitialized     
            else  
                reinit = pop.lb + rand(pop.size, pop.dimensions).*(pop.ub-pop.lb);
                if any(outsiders1(:) | outsiders2(:))
                    newpop(outsiders1) = reinit(outsiders1);
                    newpop(outsiders2) = reinit(outsiders2);
                    % also remove any function values
                    newfit(any(outsiders1,2), :) = NaN;
                    newfit(any(outsiders2,2), :) = NaN;
                end
            end % if            
        end % function
        
        % initialize algorithms
        function initialize_algorithms(pop)
            
            % PSO
            if strcmpi(pop.algorithm, 'PSO')
                                
                % initialize velocities
                % (average velocities about 20% of [[lb]-[ub]] interval)
                pop.pop_data.velocities = randn(pop.size, pop.dimensions) .* (pop.ub-pop.lb)/5;
                                
                % rename for clarity
                NumNeighbors = pop.options.PSO.NumNeighbors;
                
                % initialize neighbors                 
                switch lower(pop.options.PSO.NetworkTopology)
                    
                    case 'star'
                        % star topology - in each star, there is one 
                        % focal particle, to which the other members 
                        % of the star are connected. 
                        
                        % initialize
                        all_particles = (1:pop.size).';
                        num_stars     = floor(pop.size/NumNeighbors);
                        pop.pop_data.neighbors = zeros(pop.size, NumNeighbors-1);
                        
                        % initialize stars
                        if num_stars ~= 0                            
                            [dummy, focals] = sort(rand(pop.size,1));                            
                            focals = all_particles(focals(1:num_stars));
                            all_particles(focals) = [];
                            % select [NumNeighbors] random & unique neighbors
                            % for each focal particle
                            for i = 1:num_stars
                                % select new neighbors
                                [dummy, inds] = sort(rand(size(all_particles,1),1)); %#ok
                                new_neighs = all_particles(inds(1:NumNeighbors-1));
                                % adjust array
                                for j = 1:NumNeighbors-1
                                    all_particles(all_particles == new_neighs(j)) = [];
                                end
                                % assign new neighbors to focal particle
                                pop.pop_data.neighbors(focals(i), :) = new_neighs;
                                % assign focal particle to new neighbors
                                pop.pop_data.neighbors(new_neighs, 1) = focals(i);
                            end
                        else
                        end
                        
                        % population might be badly scaled for selected
                        % number of stars. Correct for this
                        % TODO - it works; those particles simply have no neighbors
                        
                        
                    case 'ring'
                        
                        % initialize
                        all_particles = (1:pop.size).';
                        num_rings     = floor(pop.size/NumNeighbors);
                        pop.pop_data.neighbors = zeros(pop.size, 2);
                        
                        % form the ring
                        for i = 1:num_rings
                            % randomly select [NumNeighbors] particles
                            [dummy, inds] = sort(rand(size(all_particles,1),1)); %#ok
                            new_neighs = all_particles(inds(1:NumNeighbors));
                            % insert circularly shifted arrays
                            pop.pop_data.neighbors(new_neighs, :) = ...
                                [circshift(new_neighs,1), circshift(new_neighs,-1)];
                             % adjust array                             
                             for j = 1:NumNeighbors
                                 all_particles(all_particles == new_neighs(j)) = [];
                             end                             
                        end                            
                            
                        % population might be badly scaled for selected
                        % number of rings. Correct for this
                        % TODO - it works; those particles simply have no neighbors
                
                            
                    case 'fully_connected'
                        % fully connected swarm - all particles have 
                        % ALL other particles as neighbor
                        
                        % initialize
                        pop.pop_data.neighbors = zeros(pop.size, pop.size-1);
                        
                        % fill the neighbors                        
                        for i = 1:pop.size
                            pop.pop_data.neighbors(i, :) = [1:i-1, i+1:pop.size];                            
                        end                        
                        
                end % switch
                                
                % find global best solution
                [global_best, index] = min(pop.fitnesses);
                pop.pop_data.global_best_ind = pop.individuals(index, :);
                pop.pop_data.global_best_fit = global_best;
                
                % initially, local best solutions are the function values themselves
                pop.pop_data.local_best_inds = pop.individuals;
                pop.pop_data.local_best_fits = pop.fitnesses;
                
                % find the neighbor best 
                pop.pop_data.neighbor_best_fits = zeros(pop.size,1);
                pop.pop_data.neighbor_best_inds = zeros(pop.size, pop.dimensions);
                for i = 1:pop.size
                    neighbors = pop.pop_data.neighbors(i, :);
                    neighbors = neighbors(neighbors ~= 0);
                    if isempty(neighbors), continue, end
                    [neighbor_best, ind] = min(pop.fitnesses(neighbors));
                    pop.pop_data.neighbor_best_fits(i, 1) = neighbor_best;
                    pop.pop_data.neighbor_best_inds(i, :) = pop.individuals(ind, :);
                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 
                
    end % methods (private)
end % classdef