function [children] = crossover(pop_keep, ma,pa, N_total_children, method)
% BINARY_GA - Implements a couple of crossover algorithms
% 
% Inputs: 
%   pop_keep the population of most fit individuals 
%   ma = indices into pop_keep for the "mother" chromosome
%   pa = indices into pop_keep for the "father" chromosome
%        each pairing produced two children, we repeatidly cycle through 
%        the parents performing crossover until we have enough children 
%   N_total_children = number of total children that we must generate 
%   method = method used to perform crossover
%   
% Outputs:
% 
% 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.
% 
%-----

[N_keep,N_bits] = size(pop_keep);
nPairs          = length(ma); 

children = zeros(N_total_children,N_bits); % initialize storage
switch method 
 case 1, % single point crossover
  crossover_points = ceil( (N_bits-1) * rand(1,N_total_children) );

  child_number = 1;
  while( child_number<N_total_children )
    
    if( child_number <= nPairs ) % we cycle through parents breading 2 children repeatedly until we get enough children 
      ind = child_number;
    else
      ind = mod( child_number, nPairs ); 
    end    
    mother = pop_keep( ma( ind ), : ); 
    father = pop_keep( pa( ind ), : ); 
    
    indsL = 1:crossover_points(child_number);
    indsR = (crossover_points(child_number)+1):N_bits;

    children(child_number,:) = [ mother(indsL), father(indsR) ];
    child_number = child_number + 1; 
    
    if( child_number<=N_total_children )
      children(child_number,:) = [ father(indsL), mother(indsR) ]; 
      child_number = child_number + 1; 
    end 
  end
 case 2, % double point crossover
  crossover_points_1 = ceil( (N_bits-1) * rand(1,N_total_children) );
  crossover_points_2 = ceil( (N_bits-1) * rand(1,N_total_children) );

  child_number = 1;
  while( child_number<N_total_children )
    
    if( child_number <= nPairs ) % we cycle through parents breading 2 children repeatedly until we get enough children 
      ind = child_number;
    else
      ind = mod( child_number, nPairs ); 
    end    
    mother = pop_keep( ma( ind ), : ); 
    father = pop_keep( pa( ind ), : ); 
    
    indsL = min( crossover_points_1(child_number), crossover_points_2(child_number) );
    indsR = max( crossover_points_1(child_number), crossover_points_2(child_number) );

    children(child_number,1:indsL)          = mother(1:indsL);
    children(child_number,(indsL+1):indsR)  = father((indsL+1):indsR);
    children(child_number,(indsR+1):N_bits) = mother((indsR+1):N_bits);
    child_number = child_number + 1; 
    
    if( child_number<=N_total_children )
      children(child_number,1:indsL)          = father(1:indsL);
      children(child_number,(indsL+1):indsR)  = mother((indsL+1):indsR);
      children(child_number,(indsR+1):N_bits) = father((indsR+1):N_bits);
      child_number = child_number + 1; 
    end 
  end
 case 3, % uniform crossover 
  mask = round( rand(N_total_children,N_bits) ); % 0's get bits from the mother and 1's get bits from the father

  child_number = 1;
  while( child_number<N_total_children )
    
    if( child_number <= nPairs ) % we cycle through parents breading 2 children repeatedly until we get enough children 
      ind = child_number;
    else
      ind = mod( child_number, nPairs ); 
    end    
    mother = pop_keep( ma( ind ), : ); 
    father = pop_keep( pa( ind ), : ); 
    
    indsM = find( mask(child_number,:)==0 );
    indsF = find( mask(child_number,:)==1 );

    children(child_number,indsM) = mother(indsM);
    children(child_number,indsF) = father(indsF);
    child_number = child_number + 1; 
    
    if( child_number<=N_total_children )
      children(child_number,indsM) = father(indsM);
      children(child_number,indsF) = mother(indsF);
      child_number = child_number + 1; 
    end 
  end
 otherwise
  error('unknown value for method');
end