function [pop_natural,pop_costs,best_fn_values,avg_fn_values] = continuous_GA(wfn, bounds, ...
          N_pop, X_rate, mu, max_number_of_iterations, parent_selection_method, crossover_method, use_hybrid_GA )
% CONTINUOUS_GA - Implements a continuous genetic algorithm for function minimization
% 
% Inputs: 
%   wfn = a function handle accepting a double vector input and returning a scalar output representing a function to minimize
%         inputs to wfn are in the "natural domain of the problem" and are not normalized variables
%         (which are used in the genetic operators below)
%   bounds = [ number_of_variables x 2 ] input matrix where the ith row is the valid bounds for the ith variable
%   
% Outputs:
%   x_bin = [ number_of_samples x N_var ] continuous encoding of each variable in x
%   pop_natural = final iteration input variables for each member of the population 
%   pop_costs = final iteration cost of each member of the population 
%   best_fn_values = best function value (smallest) at each iteration
%   avg_fn_value = average function value at each iteration
% 
% Written by:
% -- 
% John L. Weatherwax                2006-08-28
% 
% email: wax@alum.mit.edu
% 
% Please send comments and especially bug reports to the
% above email address.
% 
%-----

plot_iterations = true;

% Parse the parameters of the genetic algorithm (if not provided defaults are used)
% 
if( nargin < 3 ) 
  N_pop = 100; % the population size
end
if( nargin < 4 )
  X_rate = 0.5; % selection rate 
end
if( nargin < 5 ) 
  mu = 0.02; % the mutation rate
end
if( nargin < 6 )
  max_number_of_iterations = 20; 
end
if( nargin < 7 ) 
  parent_selection_method = 3; % method to select mates 
end
if( nargin < 8 )
  crossover_method = 1; % method used to do cross over
end
% Parameters used with a hybrid GA:
if( nargin < 9 )
  use_hybrid_GA = 0; % by default we don't use this feature (turn this on if you want to use it)
end
hybrid_GA_n_iters     = 5; % every this many outer iterations we will stop and then do local solves
hybrid_GA_n_solutions = 3; % the number of the best solutions we will try to refine local

% Check inputs:
N_var = size(bounds,1); 
assert( size(bounds,2)==2, 'bounds input must have two columns' );

% Generate the initial population of normalized variables:
pop = rand(N_pop,N_var); 

N_keep = floor( X_rate * N_pop );

number_of_print_statements = 10; % print this many "iteration ..." lines
iter_mod = floor( max_number_of_iterations / number_of_print_statements );

number_of_iterations = 0;
best_fn_values = zeros(max_number_of_iterations,1);
avg_fn_values  = zeros(max_number_of_iterations,1);
while( true )
  % Decode normalized variables into natural variables (the variables from the problem definition): 
  pop_natural = var_decode(pop,bounds);
  
  % Find cost for each variable: 
  pop_costs = wfn( pop_natural );
  
  % Sort the costs and order the population so that the most fit member (most negative function evaluation) are at the top:
  [pop_costs,inds] = sort(pop_costs(:),1,'ascend');
  pop              = pop(inds,:);

  % compute statistics on this population:
  best_fn_values(number_of_iterations+1) = pop_costs(1);
  avg_fn_values(number_of_iterations+1) = mean( pop_costs(1:N_keep) ); % for some objective function, mutation can produce very poor solutions that skew the value of the mean thus we evaluate this only over the samples that will be kept on this itertion

  if plot_iterations & mod(number_of_iterations,iter_mod) == 0 
    fprintf('iteration %8d; best_fn= %16.6f; avg_fn= %16.6f\n',number_of_iterations,best_fn_values(number_of_iterations+1),avg_fn_values(number_of_iterations+1));
  end
  
  % check for convergence for quick exit (since population is sorted by the objective function at this point): 
  number_of_iterations = number_of_iterations + 1; 
  if( number_of_iterations >= max_number_of_iterations )
    break;
  end

  % Use a hybrid GA if requested: 
  if( use_hybrid_GA & mod(number_of_iterations,hybrid_GA_n_iters) == 0 ) % time to perform local solves
    addpath('../SHIP'); 
    nelder_eps = 1.e-6; 
    for ii=1:hybrid_GA_n_solutions
      %
      % Call the UNCONSTRAINED optimization routine "nelder".  
      % 
      % This is an unconstrained routine and as such can output optimal solutions that don't
      % necessarily satisfy the bounds imposed by the problem domain.  This requires us to determine
      % ourselves which simplex vertex has the smallest CONSTRAINED functional values.
      %
      simplex0                 = pop( ii + [0:N_var], : ) + nelder_eps * rand( N_var+1, N_var ); % the initial simplex guess: [N_var+1 x N_var]
      simplex0_nat             = var_decode(simplex0,bounds); % the initial simplex guess in natural coordinates; and in input required by SHIP optimization functions: [N_var+1,N_var] 
      wfn_convert_to_row_input = @(x) ( wfn( x(:)' ) ); % this function wrapper converts inputs vectors to row vectors 
      simplex1_nat             = nelder( simplex0_nat', wfn_convert_to_row_input, nelder_eps, 100 )'; % output is [N_var+1 x N_var].  UNCONSTRAINED variable output
                                                                                                      
      simplex1                 = var_encode(simplex1_nat,bounds); % [N_var+1,N_var].  Solutions now in range [0,1] and truncated to within the problem domain
      simplex1_nat             = var_decode(simplex1,bounds); % truncated variables in the natural domain 
      fs                       = wfn( simplex1_nat );
      [fx_min,fx_min_ind]      = min(fs); % find the vertex with the smallest constrained function value
      x_min                    = simplex1(fx_min_ind,:);
      if( fx_min < pop_costs(ii) )
        pop(ii,:)     = x_min;
        pop_costs(ii) = fx_min;
      end
    end
    [pop_costs,inds] = sort(pop_costs(:),1,'ascend'); % resort (the ordering of the best chromosome may have changed because of this)
    pop              = pop(inds,:);
  end
  
  % Natural Selection: (could also perform thresholding): 
  pop_keep       = pop(1:N_keep,:);
  pop_keep_costs = pop_costs(1:N_keep);
  C_n_keep_p_1   = pop_costs(N_keep+1);

  % Compute indices for the mother (ma) and father (pa): 
  [ma,pa] = select_mates( pop_keep_costs, parent_selection_method, C_n_keep_p_1 );
  
  % Mate the chromosomes to produce children: 
  children = crossover( pop_keep, ma, pa, N_pop-N_keep, crossover_method, wfn, bounds );
  
  pop = [ pop_keep; children ]; % the new populaton before mutation
  
  % Mutation (keep 2 chromosomes untouched ... this is elitism): 
  pop = mutate( pop, mu, 2 );  
end