function [sol, fval, evals] = genetic(varargin)
%GENETIC                         Optimization using genetic algorithm
%
%   GENETIC(func, popsize, lb, ub) tries to find the global optimum of 
%   the fitness-function [func], using a basic genetic algorithm. The 
%   population size set by [popsize], and the boundaries for each 
%   dimension are set by the vectors [lb] and [ub], respectively. 
%
%   [sol, fval, evals] = GENETIC(...) returns the trial vector found to 
%   yield the global minimum in [sol], and the corresponding function value 
%   by [fval]. The total amount of function evaluations that the algorithm
%   performed is returned in [evals].
%
%   The function [func] must accept a number of input arguments [N], equal to 
%   the number of elements in the vectors [lb] or [ub]. The function must be 
%   vectorized so that inserting a matrix of [popsize]x[N] will return a 
%   vector of size [popsize] containing the corresponding function values
%   for the [N] trial vectors.
%
%   The default control parameters GENETIC uses, are
%
%       crossover probability = 0.25    Probability that an individual will 
%                                       perform a crossover
%
%       mutation probability  = 0.01    Probability that an individual gene
%                                       will mutate
%
%       convergence = 100   Maximum amount of iterations the best
%                           individual ever found should remain the best
%                           individual before the algorithm converges. This
%                           is the default convergence criterion -- see
%                           HEURSET for other options.
%
%   These values, as well as additional options, may be given manually in 
%   any extra input variables. Additionally, GENETIC may be called as 
%   GENETIC(PROBLEM), where [PROBLEM] is a complete problem-structure 
%   constructed by HEURSET.
%
%   Some optimizations require repeated evaluation of a function that takes 
%   in the order of hours per evaluation. In those cases, it might be 
%   convenient to interrupt the optimization after a few days and use the 
%   results thus far obtained. Unfortunately, usually after you interrupt a 
%   function, its results are lost. For that reason, GENETIC creates two 
%   global variables, [GENETIC_bestind] and [GENETIC_bestfval], which 
%   store the best individual and function value thus far found. When the
%   optimization is interrupted, the intermediate results from the
%   optmization are still available. Once an optimization has succesfully
%   converged, these variables are deleted from the workspace again.
%   
%   See also heurset, diffevolve, swarm, multiswarm, simanneal, godlike.

% Author: Rody P.S. Oldenhuis
% Delft University of Technology
% E-mail: oldenhuis@dds.nl
% Last edited: 02/Feb/2008


%%  initialise & parse input   
          
    % initial values
    skippop  = false;
    problem  = heurset;  
    [pop, func, popsize, lb, ub, grace, display, maxfevals, convmethod,...
     convvalue, mutationprob, crossprob] = parseprob(problem);
    
    % define temporary globals
    global GENETIC_bestfval GENETIC_bestind
 
    % common inputs
    if (nargin >= 4)
        func     =  varargin{1};
        popsize  =  varargin{2}; 
        lb       =  varargin{3};
        ub       =  varargin{4};        
    end
    
    % additional inputs
    if (nargin >= 5) 
        if ~isstruct(varargin{5})
            options = heurset(varargin{5:end});
        else
            options = varargin{5};                
        end
        [dum1, dum2, dum3, dum4, dum5, grace, display, maxfevals, convmethod,...
         convvalue, mutationprob, crossprob] = parseprob(options);
    end
    
    % if called from GODLIKE
    if (nargin == 2)  
        problem = varargin{2};  
        % errortrap
        if ~isstruct(problem)
            error('PROBLEM should be a structure. Type `help heurset` for details.')
        end 
        [pop, func, popsize, lb, ub, grace, display, maxfevals, convmethod,...
         convvalue, mutationprob, crossprob] = parseprob(problem);
        
        % make sure the options are correct        
        convmethod = 'maxiters'; 
        grace   = 0;
        display = false;
        skippop = true;
    end
    
    % if given a full problem structure
    if (nargin == 1)
        problem  = varargin{1};        
        % errortrap
        if ~isstruct(problem)
            error('PROBLEM should be a structure. Type `help heurset` for details.')
        end 
        [dum1, func, dims, popsize, lb, ub, grace, display, maxfevals, convmethod,...
                  convvalue, mutationprob, crossprob] = parseprob(problem);
    end  
    
    % initialize convergence method, setting the default values if omitted
    if strcmpi(convmethod, 'exhaustive')
        convergence  = 1;
        maxiterinpop = convvalue;
        maxiters     = inf; 
    elseif strcmpi(convmethod, 'maxiters')
        convergence  = 2;
        maxiterinpop = inf;
        maxiters     = convvalue; 
    elseif strcmpi(convmethod, 'achieveFval')
        convergence  = 3;          
        %errortrap
        if isempty(convvalue)
            error('Please define function value to achieve.')
        end        
        maxiterinpop = inf;
        maxiters     = inf;        
    else
        convergence  = 1;
        maxiterinpop = convvalue;
    end  
    
    % problem dimensions
    dims = size(lb, 2);
    
    % errortraps
    if ( (size(lb, 2) ~= 1 || size(ub, 2) ~= 1) &&...
         (size(lb, 2) ~= dims || size(ub, 2) ~= dims) )
        error('Upper- and lower boundaries must be the same size.')
    end
    if (popsize < 2)
        error('GENETIC needs a population size of at least 2.')
    end  
    
%%  initial values

    % iteration-zero values
    oldbestfit   = inf;
    oldbestind   = zeros(1, dims);
    iters        = 0;
    evals        = 0;
    improvement  = maxiterinpop;
    firsttime    = true;
    converging1  = false;
    converging2  = false;
    
    % boundaries
    if (numel(lb) == 1)
        mins   = repmat(lb, popsize, dims);
        maxs   = repmat(ub, popsize, dims); 
    end
    if (numel(lb) == numel(ub) && size(lb, 2) == dims)
        mins   = repmat(lb, popsize, 1);
        maxs   = repmat(ub, popsize, 1);                        
    end
    range  = (maxs - mins); 
    
    % Display strings
    if display        
        fprintf(1, ['\nNote that when GENETIC is interrupted, the best solution and function\n',...
                    'value are saved in global variables GENETIC_bestind and GENETIC_bestfval.\n\n']);
        fprintf(1, 'Optimizing with ');
        switch convergence
            case 1
            	fprintf(1, 'exhaustive convergence criterion...\n');
            case 2
                fprintf(1, 'maximum iterations convergence criterion...\n');
            case 3
                fprintf(1, 'goal-attain convergence criterion...\n');
        end
        fprintf(1, 'Current best solution has cost of   ------------- ');
        overfevals1  = '\n\nCould not find better solution within given maximum';
        overfevals2  = '\nof function evaluations. Increase MaxFevals option.\n\n';
        fitbackspace = '\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b';
        converge1bck = [fitbackspace, '\b\b\b\b\b\b\b\b\b\b\b'];
        converge2bck = [fitbackspace, '\b\b\b\b\b\b\b'];
        convergestr  = ', possibly converging: ';
        itersstr     = ', iterations left: ';
    end
    
%%  initialize pop
    if (~skippop)
        pop = repmat(ub, popsize, 1) - repmat(ub-lb, popsize, 1).* rand(popsize, dims);
    end
    
%%  Genetic Algorithm loop

    % replication matrices for vectorization
    poprep = ones(1, popsize);
    popvec = (1:popsize)';
    
    % loop while one of the convergence criteria is not met
    while (improvement > 0 && iters <= maxiters)
    
        % evaluate and scale fitnesses
        try
            fitnesses  = feval(func, pop);
        catch
            error('genetic:fevalerror', ...
                 ['Genetic cannot continue because the supplied cost function ',...
                  'gave the following error:\n %s'], lasterr);
        end
        if (numel(fitnesses) ~= popsize)
            error('genetic:function_not_vectorized_or_incorrect_dimensions', ...
                 ['The user-supplied cost function does not appear to get enough arguments,\n',...
                  'or is not vectorized properly. Make sure the dimensions of the limits [Lb]\n',...
                  'and [Ub] are the same as the required input vector for the function, or that\n',...
                  'the function is vectorized properly.'])            
        end
        fitproper  = fitnesses;
        fitnesses  = fitnesses ./ sum(fitnesses);       
        fitcum     = cumsum(fitnesses);
        [bestindfit, inds] = min(fitproper);
        bestind    = pop(inds, :);

        % increase function evaluations counter
        evals = evals + numel(fitnesses);
                    
        % improving solutions
        if (bestindfit <= oldbestfit)
            
            % new best individual
            oldbestind  = bestind;
            oldbestfit  = bestindfit; 
            improvement = maxiterinpop;
            
            % assign also the globals
            GENETIC_bestfval = bestindfit;
            GENETIC_bestind  = bestind;
            
            % display progress