% -------------------------------------------------------------
% GA FOR MIXTURE MODELS, The function implementing the BLX crossover
% Copyright (c) 2005 - 2006 Jussi Tohka (jussi.tohka@tut.fi)
% Institute of Signal Processing, Tampere University of Technology, Finland
% and
% Laboratory of Neuro Imaging, David Geffen School of Medicine,
% University of California, Los Angeles, CA, USA 

% The algorithm is described in 
% J. Tohka, E. Krestyannikov, I. Dinov, D. Shattuck, U. Ruotsalainen,
% A.W. Toga. Genetic algorithms for finite mixture model based tissue
% classification in brain MRI. In  Proc. of European Medical and
% Biological Engineering Conference, IFMBE Proceedings vol. 11, pp. 4077
% - 4082, Prague, Czech Republic, 2005.
% ----------------------------------------------------------------
% Permission to use, copy, modify, and distribute this software 
% for any purpose and without fee is hereby
% granted, provided that the above copyright notice appear in all
% copies.  The author, Tampere University of Technology and University of
% California, Los Angeles make no representations
% about the suitability of this software for any purpose.  It is
% provided "as is" without express or implied warranty.
% ---------------------------- ---------------------------------
% 
% --------------------------------------------------------
% Input: old_pop     : old population by selection routine
%              n     : how many individuals to generate. must be an
%                      even number
%          alpha     : exploration
%     upperlimit     : upperbound for each parameter
%     lowerlimit     : lowerbound for each parameter 
%     sum_constraint : indexces of variables that should sum to one                     
%     pph            : type of handling equality constraint. if 1, it is assumed that the last variable is subject to the constraint

function new_pop = mixturega_BLXxover(old_pop,n,alpha,upperlimit,lowerlimit,sum_constraint,pph)
 
  sz = size(old_pop);
  % sz(2) = popsize
  % sz(1) = number of variables  
  if pph == 0
    
  %  mates = unidrnd(sz(2),n,2);
    mates =  ceil(sz(2) .* rand([n 2]));
    new_pop = zeros(sz(1),n);  
    b = rand(sz(1),n);
    b = b*(1 + 2*alpha) - alpha;
    new_pop = old_pop(:,mates(:,1)) + b.*(old_pop(:,mates(:,2)) - old_pop(:,mates(:,1)));
    new_pop = max(min(new_pop,repmat(upperlimit,n,1)'),repmat(lowerlimit,n,1)');
    new_pop(sum_constraint,:) = new_pop(sum_constraint,:)./(repmat(sum(new_pop(sum_constraint,:)),length(sum_constraint),1));
    % make sure there is no NaNs
    new_pop = isnan(new_pop)*(1/length(sum_constraint)) + (~isnan(new_pop)).*new_pop;
  else
    m = sz(1) -1;
    l = length(sum_constraint) - 1;
  %  mates = unidrnd(sz(2),n,2);
    mates =  ceil(sz(2) .* rand([n 2]));
    new_pop = zeros(sz(1),n);  
    b = rand(sz(1) - 1,n);
    b = b*(1 + 2*alpha) - alpha;
    new_pop(1:m,:) = old_pop(1:m,mates(:,1)) + b.*(old_pop(1:m,mates(:,2)) - old_pop(1:m,mates(:,1)));
    new_pop(1:m,:) = max(min(new_pop(1:m,:),repmat(upperlimit(1:m),n,1)'),repmat(lowerlimit(1:m),n,1)');
    
    new_pop(sum_constraint(1:l),:) =  new_pop(sum_constraint(1:l),:) - repmat((sum(new_pop(sum_constraint(1:l),:)) > 1).*(sum(new_pop(sum_constraint(1:l),:)) - 1)/l,l,1);
    mi = min(new_pop(sum_constraint(1:l),:));  
    new_pop(sum_constraint(1:l),:) = new_pop(sum_constraint(1:l),:) + sign(new_pop(sum_constraint(1:l),:)).*repmat((mi < 0).*mi,l,1);
    new_pop(m + 1,:) = 1 - sum(new_pop(sum_constraint(1:l),:));
  end