function [children] = crossover(pop_keep, ma,pa, N_total_children, method, wfn, bounds)
% CROSSOVER - 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
%           (the first three methods are the same as thoes described in the chapter on binary genetic algorithms)
%   
% 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_vars] = size(pop_keep);
nPairs          = length(ma); 

children = zeros(N_total_children,N_vars); % initialize storage
switch method 
 case 1, % single point crossover
  crossover_points = ceil( (N_vars-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 ) + 1; 
    end    
    mother = pop_keep( ma( ind ), : ); 
    father = pop_keep( pa( ind ), : ); 
    
    indsL = 1:crossover_points(child_number);
    indsR = (crossover_points(child_number)+1):N_vars;

    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_vars-1) * rand(1,N_total_children) );
  crossover_points_2 = ceil( (N_vars-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 ) + 1; 
    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_vars) = mother((indsR+1):N_vars);
    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_vars) = father((indsR+1):N_vars);
      child_number = child_number + 1; 
    end 
  end
 case 3, % uniform crossover 
  mask = round( rand(N_total_children,N_vars) ); % 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 ) + 1; 
    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
 case 4, % blending method (3.9)
  crossover_points = ceil( (N_vars-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 ) + 1; 
    end    
    mother = pop_keep( ma( ind ), : ); 
    father = pop_keep( pa( ind ), : ); 
    
    indsL = 1:crossover_points(child_number);
    indsR = (crossover_points(child_number)+1):N_vars;

    beta = rand();
    
    children(child_number,:) = [ beta*mother(indsL), (1.-beta)*father(indsR) ];
    inds = children(child_number,:) < 0.; children(child_number,inds) = 0.; % respect variable boundaries 
    inds = children(child_number,:) > 1.; children(child_number,inds) = 1.;
    child_number = child_number + 1; 
    
    if( child_number<=N_total_children )
      children(child_number,:) = [ (1.-beta)*father(indsL), beta*mother(indsR) ]; 
      inds = children(child_number,:) < 0.; children(child_number,inds) = 0.; % respect variable boundaries 
      inds = children(child_number,:) > 1.; children(child_number,inds) = 1.;
      child_number = child_number + 1; 
    end 
  end
 case 5, % blending method (3.13) or the books suggested method: we blend two variables at a specific point and cross all other variables
  alpha_locations = ceil( N_vars * 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 ) + 1; 
    end    
    mother = pop_keep( ma( ind ), : ); 
    father = pop_keep( pa( ind ), : ); 
    
    indsL = 1:(alpha_locations(child_number)-1);
    indsC = alpha_locations(child_number);
    indsR = (alpha_locations(child_number)+1):N_vars;

    beta = rand();
    nv = mother(indsC) - beta*(mother(indsC)-father(indsC));
    if( nv < 0.0 ) nv = 0.0; end % restrict boundaries 
    if( nv > 1.0 ) nv = 1.0; end 
    
    children(child_number,:) = [ mother(indsL), nv, father(indsR) ];
    child_number = child_number + 1; 
    
    if( child_number<=N_total_children )
      nv = mother(indsC) + beta*(mother(indsC)-father(indsC));
      if( nv < 0.0 ) nv = 0.0; end % restrict boundaries 
      if( nv > 1.0 ) nv = 1.0; end 
      children(child_number,:) = [ father(indsL), nv, mother(indsR) ]; 
      child_number = child_number + 1; 
    end 
  end
 otherwise
  error('unknown value for method');
end