function [theResult, theFinalCost, theInitialCost] = ...
			genetic(theFcn, theParameters, theOptions, empty, varargin)

% genetic -- Genetic optimization driver.
%  genetic('demo') demonstrates itself with 10 parameters.
%  genetic(N) demonstrates itself with N parameters (default = 10).
%   The procedure combines random vectors genetically, thereby
%   redistributing their elements, in hopes of reducing the
%   standard-deviation of each.
%  genetic('theFcn', theParameters, theOptions, [], varargin) minimizes
%   theFcn cost-function, using theParameters, a two-column array
%   of [xmin xmax], the permitted ranges.  TheOptions vector is:
%   [#-of-parents; #-of-children; maximum-iterations; maximum-cost].
%   The cost-function will be called repeatedly, using all of the
%   evolving parameters at once, followed by the "varargin" arguments,
%   if any, until halted by reaching the maximum number of iterations
%   or by falling below maximum-cost.  (The probabilistic genetic
%   actions take place in the "breed" function.)  The Fcn input
%   is one array, whose columns represent individual solutions,
%   plus any additional arguments specified by "varargin".
 
% Copyright (C) 2001 Dr. Charles R. Denham, ZYDECO.
%  All Rights Reserved.
%   Disclosure without explicit written consent from the
%    copyright owner does not constitute publication.
 
% Version of 29-Mar-2001 17:08:53.
% Updated    14-Jun-2001 07:26:58.

if nargin < 1
	help(mfilename)
	theFcn = 'demo';
end

if isequal(theFcn, 'demo')
	theFcn = 10;
end

if nargin < 2 & ischar(theFcn)
	theFcn = eval(theFcn);
end

if ~ischar(theFcn)
	nParameters = theFcn(1);
	nParents = 3*nParameters;
	nChildren = 10*nParents;
	if length(theFcn) > 1
		nParents = theFcn(2);
	end
	if length(theFcn) > 2
		nChildren = theFcn(3);
	end
	theFcn = 'std';
	theParameters = zeros(nParameters, 2);
	theParameters(:, 2) = 1;
	theMaximumIterations = 100;
	theMaximumCost = 1e-4;
	theOptions = [nParents; nChildren; theMaximumIterations; theMaximumCost];
	tic
	[result, theFinalCost, theInitialCost] = ...
			feval(mfilename, theFcn, theParameters, theOptions);
	t = toc;
	disp([' ## Elapsed time: ' num2str(ceil(10*t)/10) ' s.'])
	hold off
	plot(1:length(theFinalCost), theFinalCost, '-bo', ...
				1:length(theInitialCost), theInitialCost, '-ro')
	xlabel('Index')
	ylabel('Standard-deviation')
	legend final initial
	set(gca, 'XLim', [0 length(theFinalCost)+1])
	set(gcf, 'Name', ['Genetic ' int2str(nParameters)])
	return
end

nParents = theOptions(1);
nChildren = theOptions(2);
theMaximumIterations = theOptions(3);
theMaximumCost = theOptions(4);

[m, n] = size(theParameters);
nParameters = m;

sz = [1 nParents+nChildren];
xmin = repmat(theParameters(:, 1), sz);
xmax = repmat(theParameters(:, 2), sz);
xdif = xmax - xmin;

xmin_p = xmin(:, 1:nParents);
xdif_p = xmax(:, 1:nParents) - xmin(:, 1:nParents);
xmin_c = xmin(:, 1:nChildren);
xdif_c = xmax(:, 1:nChildren) - xmin(:, 1:nChildren);

xval = rand(nParameters, nParents) .* xdif_p + ...
			xmin_p;

% Evaluate the parents.

if isempty(varargin)
	theInitialCost = feval(theFcn, xval);
else
	theInitialCost = feval(theFcn, xval, varargin{:});
end

% Breed many generations by genetic-crossover
%  of the parameter-vectors.  Once debugged,
%  we will add a small amount of genetic-drift,
%  to the mix, plus an even smaller amount of
%  gross genetic-mutation.

% Need to trap invalid costs, such as NaN or negative.

iter = 0;
cost = theInitialCost;

while iter < theMaximumIterations & all(cost > theMaximumCost)
	if iter > 0 & rem(iter, 100) == 0
		disp([' ## iter: ' int2str(iter)])
		cost
	end
	iter = iter+1;
	xval = (xval(:, 1:nParents) - xmin(:, 1:nParents)) ./ ...
				xdif(:, 1:nParents);
	xnew = breed(xval, nChildren);
	if (0)
		xnew = drift(xnew);
		xnew = mutate(xnew);
	end
	xnew = xnew .* xdif(:, 1:nChildren) + ...
				xmin(:, 1:nChildren);
	if isempty(varargin)
		newCost = feval(theFcn, xnew);
	else
		newCost = feval(theFcn, xnew, varargin{:});
	end

% We need to compute a fitness, then use
%  it as the probability of survival.

	cost = [cost newCost];
	xval = [xval xnew];
	
	fitness = exp(-cost);
	cumulative = cumsum(fitness);
	cumulative = cumulative ./ cumulative(end);   % Normalize.
	cumulative(end) = 1;
	
	index = zeros(size(cumulative));
	
% The following is complicated by a "while" loop,
%  protecting the search-scheme, until we discover
%  how to make it bullet-proof.  We want "f" never
%  to be empty.
	
	okay = ~~0;
	while ~okay
		okay = ~~1;
		r = sort(rand(size(cumulative)));
		for i = 1:length(r)
			f = find(r(i) <= cumulative);   % Very slow: 40%
			if isempty(f)
				okay = ~~0;
				break
			else
				index(i) = f(1);
			end
		end
	end
	f = find(diff(index) == 0);
	if any(f)
		index(f) = [];
	end
	
	cost = cost(index);
	xval = xval(:, index);
	
	[cost, index] = sort(cost);
	
	xval = xval(:, index);
	xval = xval(:, 1:nParents);
	cost = cost(1:nParents);
end

if nargout > 0
	theResult = xval;
	theFinalCost = cost;
else
	result = xval
	cost
end


% ---------- breed ---------- %


function progeny = breed(parents, nChildren)

% breed -- Parents breed children.
%  breed(parents, nChildren) breds nChildren from
%   the given parents.

DRIFT_PROB = 1e-4;     % See the "drift" function below.
DRIFT_AMOUNT = 1e-4;   % See the "drift" function below.
MUTATE_PROB = 1e-4;    % See the "mutate" function below.

[m, n] = size(parents);
progeny = zeros(m, nChildren);
for i = 1:2:nChildren
	k = floor(rand(1, 2) .* n + 1);
	progeny(:, [i i+1]) = crossover(parents(:, k));   % Slow: 10%
end
progeny = drift(progeny, DRIFT_PROB, DRIFT_AMOUNT);
progeny = mutate(progeny, MUTATE_PROB);


% ---------- crossover ---------- %


function y = crossover(x)

% crossover -- Genetic crossover.
%  crossover(x) combines the two columns of x
%   by genetic crossover.

[m, n] = size(x);

i = (rand(m, 1) < 0.5);
j = ~i;

y = zeros(size(x));
y(:, 1) = x(:, 1) .* i + x(:, 2) .* j;   % Slow: 8%
y(:, 2) = x(:, 1) .* j + x(:, 2) .* i;   % Slow: 8%


% ---------- drift ---------- %


function y = drift(x, prob, amount)

% drift -- Genetic drift by small amounts.
%  drift(x, prob, amount) adds genetic-drift
%   to the columns of x, by the given small
%   probability, in the given amount.

if nargin < 2, prob = 1e-4; end
if nargin < 3, amount = 1e-4; end

r = rand(size(x));
s = 2 * (rand(size(x)) - 0.5);
y = x;
y(r < prob) = y(r < prob) + s(r < prob) .* amount;
y = max(min(y, 1), 0);


% ---------- mutate ---------- %


function y = mutate(x, prob)

% mutate -- Genetic mutation by small probabilities.
%  mutate(x, prob) mutates the columns of x by the
%   given small probability.

if nargin < 2, prob = 1e-4; end

r = rand(size(x));
s = rand(size(x));
y = x;
y(r < prob) = s(r < prob);