function [ma,pa] = select_mates(pop_costs, method, C_n_keep_p_1)
% SELECT_MATES - Implements a couple of mate selection algorithms
% 
% Inputs: 
%   
% 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 = length(pop_costs);

switch method 
 case 1, % pairing from top to bottom
  ma = 1:2:N_keep;
  pa = 2:2:N_keep;
 case 2, % random pairing
  ma = ceil( N_keep*rand(1,N_keep) );
  pa = ceil( N_keep*rand(1,N_keep) );
 case 3, % weighted random pairing (rank weighting)
  probs = ( N_keep - (1:N_keep) + 1 ) ./ sum( 1:N_keep );
  ma = zeros(1,N_keep); % I'm sure there are faster ways to do this but ... 
  pa = zeros(1,N_keep);
  for ii=1:N_keep, 
    ma(ii) = find( mnrnd( 1, probs ) );
    pa(ii) = find( mnrnd( 1, probs ) );
  end
 case 4, % weighted random pairing (cost weighting) ... I'm sure there are faster ways to do this but ... 
  probs = abs( ( pop_costs - C_n_keep_p_1 ) / sum( pop_costs - C_n_keep_p_1 ) );
  ma = zeros(1,N_keep);
  pa = zeros(1,N_keep);
  for ii=1:N_keep, 
    ma(ii) = find( mnrnd( 1, probs ) );
    pa(ii) = find( mnrnd( 1, probs ) );
  end
 case 5, % tournament selection ... I'm sure there are faster ways to code this but ... 
  ma = zeros(1,N_keep);
  pa = zeros(1,N_keep);
  for ii=1:N_keep, 
    inds = ceil( N_keep*rand(1,2) ); % draw two individuals one will become the mother 
    [mv,mindex] = min( pop_costs(inds) );
    ma(ii) = inds(mindex); % we always keep the fitest of the group
    
    inds = ceil( N_keep*rand(1,2) ); % draw two individuals one will become the father 
    [mv,mindex] = min( pop_costs(inds) );
    pa(ii) = inds(mindex); % we always keep the fitest of the group
  end
 otherwise
  error('unknown value for method');
end

% send out the same number of parents 
minLen = min( length(ma), length(pa) );
ma = ma(1:minLen);
pa = pa(1:minLen); 

% make sure we don't have the same father and mother: 
inds = find( ma == pa );
nSame = length(inds);
if( nSame > 0 )
  ma(inds) = ceil( nSame*rand(1,nSame) ); % pick a mate at random 
end