function [pop_decoded,pop_costs,best_fn_values,avg_fn_values] = binary_GA(wfn, bounds, ...
         N_pop, N_gene, X_rate, mu, max_number_of_iterations, parent_selection_method, crossover_method )
% BINARY_GA - Implements a binary 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
%   bounds = [ number_of_variables x 2 ] input matrix where the ith row is the valid bounds for the ith variable
%   N_gene = number of bits used to encode a scalar representation 
%   
% Outputs:
%   x_bin = [ number_of_samples x N_var * N_gene ] binary encoding of each variable in x
% 
% 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.
% 
%-----

% 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 )
  N_gene = 100; % the discreteness of the floating point representations
end
if( nargin < 5 )
  X_rate = 0.5; % selection rate 
end
if( nargin < 6 ) 
  mu = 0.02; % the mutation rate
end
if( nargin < 7 )
  max_number_of_iterations = 20; 
end
if( nargin < 8 ) 
  parent_selection_method = 3; % method to select mates 
end
if( nargin < 9 )
  crossover_method = 1; % method used to do cross over
end

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

% Generate the initial population: 
N_bits = N_gene*N_var; 
pop    = round( rand(N_pop,N_bits) ); 

N_keep = floor( X_rate * N_pop );
  
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 chromosomes
  pop_decoded = var_decode(pop,bounds);
  
  % Find cost for each chromosome
  pop_costs = wfn( pop_decoded );
  
  % compute statistics on this population:
  best_fn_values(number_of_iterations+1) = min( pop_costs );
  avg_fn_values(number_of_iterations+1) = mean( pop_costs );
  
  % Sort the costs and order the chromosomes so that the most fit chromosomes (most negative function evaluation) are at the top:
  [pop_costs,inds] = sort(pop_costs(:),1,'ascend');
  pop              = pop(inds,:);
  pop_decoded      = pop_decoded(inds,:);
  
  % 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
  
  % 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
  children = crossover( pop_keep, ma, pa, N_pop-N_keep, crossover_method );
  
  pop = [ pop_keep; children ]; % the new populaton before mutation
  
  % Mutation
  pop = mutate( pop, mu );
end