(* :Title: Differential Evolution *)
(* :Context: DifferentialEvolution` *)

(* :Author: Mark Fisher *)

(* :Package version: 1.1, August 2000 *)

(* :Mathematica Version 4.0 *)

(* :Discussion:
This package implements the differential evolution algorithm for
continuous function optimization. The algorithm is due to Kenneth
Price and Rainer Storn.

Motto: "Only the low-cost survive."

The algorithm starts with a randomly chosen population of n real
vectors, each with d >= 2 components. A cost function is given. (Cost
is "negative fitness.") The cost of each member of the current
population is determined. Each member then becomes a "target." The
target is "mated" with a noisy vector randomly constructed from other
members of the current generation to produce an offspring. The
offspring replaces the target in the next generation if its cost is
lower.

The mate is the sum of a "base" vector (whose construction is
described below) and the scaled difference between two randomly chosen
vectors. (The scale is given by the "weighting factor.") The offspring
is constructed from the target and the mate by randomly choosing
"genes" (i.e. components) from either the target or the mate. The
probability that a gene comes from the mate is determined by the
"cross-over rate." (At least one gene always comes from the mate, even
if the cross-over rate is zero.) There are two ways to choose the base
vector: (1) the base vector is chosen randomly from the current
generation; (2) the base vector is a linear combination of the target
and the "best" vector of the current generation (in terms of cost).
The linear combination is determined by the "use-best weighting
factor."

There is no convergence criterion built in. One simply runs the
algorithm for a given number of generations. Convergence could be
based on the dispersion of the costs or the population.

*)

(* :References:
(1) Kenneth Price and Rainer Storn, "Differential Evolution,"
Dr. Dobb's Journal, April 1997.
(2) Kenneth Price and Rainer Storn, "Differential Evolution---A simple
and efficient adaptive scheme for global optimization over continuous
spaces," Working Paper TR-95-012, International Computer Science
Institute, Berkeley, CA.
*)

BeginPackage["DifferentialEvolution`",
	{"Statistics`DescriptiveStatistics`"}]

DifferentialEvolution::usage = "DifferentialEvolution[costfun, pop0, opts]
takes a function (costfun) of a d-dimensional point and an inital
population (pop0) of n points. (Thus pop0 is an n by d matrix.) The initial
population should have about 10 * d points chosen randomly.
DifferentialEvolution returns the minimum value of the function after
NumberOfGenerations. (The default setting is NumberOfGenerations -> 200.)
The option UseBest controls whether to use the best member of the
population in constructing the vector to mate with the target. Other
options are WeightingFactor, CrossoverRate, BestWeightingFactor,
PrintSummary, WriteHistory, and EvolutionHistory."

NumberOfGenerations::usage = "NumberOfGenerations is an option for
DifferentialEvolution. It specifies the number of generations to run."

WeightingFactor::usage = "WeightingFactor is an option for
DifferentialEvolution. It scales the difference between the two
randomly chosen vectors to be added to the base vector in constructing
the mate."

CrossoverRate::usage = "CrossoverRate is an option for
DifferentialEvolution. It determines the probability that an
offspring's gene comes from the mate rather than the target. At least
on gene always comes from the mate, even when the cross-over rate is
zero."

$EvolutionHistory::usage = "$EvolutionHistory is the history of values
and populations of points from the most recent use of the function
DifferentialEvolution."

EvolutionHistory::usage = "EvolutionHistory is an option for
DifferentialEvolution. It determines whether the entire history of values
and populations of points from the most recent use of the function
DifferentialEvolution is kept. The default is EvolutionHistory -> False.
For EvolutionHistory -> True, the history is assigned to the variable
$EvolutionHistory."

$GenerationNumber::usage = "$GenerationNumber is the internal
generation counter. Its value can be examined in a Dialog session."

UseBest::usage = "UseBest is an option for DifferentialEvolution. The
default setting is UseBest -> False, which combines three randomly
chosen vectors to mate with the target. UseBest -> True combines (1)
the best of the population, (2) the target, and (3) two randomly
chosen vectors to mate with the target."

BestWeightingFactor::usage = "BestWeightingFactor is an option for
DifferentialEvolution. It controls the extend to which the best of the
population is introduced into the vector that is mated with the
target. It is only used with the option UseBest -> True."

PrintSummary::usage = "PrintSummary is an option for
DifferentialEvolution. It specifies whether to print summary
statistics of each generation. The default is PrintSummary -> False."

WriteHistory::usage = "WriteHistory is an option for
DifferentialEvolution. It specifies whether to write the population
and cost history to files. The default is WriteHistory -> False."

HistoryFileName::usage = "HistoryFileName is an option for
DifferentialEvolution. It is used in conjunction with WriteHistory ->
True to specify the name of the file where the population and cost
history will be stored. The default is HistoryFileName ->
\"history.m\"."

Begin["`Private`"]

Options[DifferentialEvolution] = Sort @ {NumberOfGenerations -> 200,
WeightingFactor -> .8, CrossoverRate -> .5, UseBest -> False,
BestWeightingFactor -> .9, PrintSummary -> False, WriteHistory -> False,
HistoryFileName -> "history.m", EvolutionHistory -> False}

DifferentialEvolution::filename = "The file name `1` is not a string."

DifferentialEvolution::toosmall = "The initial population must have at least
4 members."

DifferentialEvolution[costfun_, pop0_?(MatrixQ[#, NumberQ]&),
		opts___?OptionQ] :=
	Module[{maxgen, wf, cr, usebest, lambda, hist, fun,
		print, write, file, n, d, cost0, cost, pop, pos, min},
	{maxgen, wf, cr, usebest, lambda, print, write, file, hist} =
		{NumberOfGenerations, WeightingFactor, CrossoverRate, UseBest,
			BestWeightingFactor, PrintSummary, WriteHistory, 
			HistoryFileName, EvolutionHistory} /. {opts} /. 
			Options[DifferentialEvolution];
	{usebest, print, write, hist} = TrueQ /@ {usebest, print, write, hist};
	If[!StringQ[file], Message[DifferentialEvolution::filename, file];
		Return[$Failed]];
	{n, d} = Dimensions[pop0];
	If[n < 4, Message[DifferentialEvolution::toosmall]; Return[$Failed]];
	cost0 =  costfun /@ pop0;
	$GenerationNumber = 0;
	If[print, PrintSummary[{cost0, pop0}]];
	If[write, Export[file, Flatten /@ Transpose[{cost0, pop0}], "Table"]];
	fun = If[hist, NestList, Nest];
	return = fun[NextGen[#, costfun, n, d, wf, cr, lambda, usebest, 
		print, write, file]&, {cost0, pop0}, maxgen];
	If[hist, 
		$EvolutionHistory = return; {cost, pop} = Last[return],
		{cost, pop} = return];
	min = Min[cost];
	pos = Position[cost, min][[1, 1]];
	{min, pop[[pos]]}
	]

NextGen[{cost_, pop_}, costfun_, n_, d_, wf_, cr_, lambda_,
		usebest_, print_, write_, file_] :=
	Module[{a, b, c, best, base, mate, genes, offspring, score, 
		cost1, pop1, strm},
	If[usebest,
	(* then *)
		{a, b} = ChooseAB[n];
		best = Position[cost, Min[cost]][[1,1]];
		base = Table[lambda pop[[best]], {n}] + (1 - lambda) pop,
	(* else *)
		{a, b, c} = ChooseABC[n];
		base = pop[[c]]
		];
	mate = base + wf (pop[[a]] - pop[[b]]);
	genes = Genes[n, d, cr];
	offspring = genes pop + (1 - genes) mate;
	score = costfun /@ offspring;
	{cost1, pop1} = Transpose @ MapThread[If[#1 <= #3, {#1, #2}, {#3, #4}]&,
			{score, offspring, cost, pop}];
	$GenerationNumber++;
	If[print, PrintSummary[{cost1, pop1}]];
	If[write, 
		ExportAppend[file, Flatten /@ Transpose[{cost1, pop1}], "Table"]];
	{cost1, pop1}
	]

ExportAppend[file_String, expr_, format_String] :=
	Module[{stream},
	stream = OpenAppend[file];
	WriteString[stream, ExportString[expr, format]];
	Close[stream]
	]

PrintSummary[{cost_, pop_}] := 
	Module[{labels, vals, pos, vec},
	labels = {"Generation", "Minimum", "Mean", "Standard deviation",
		"Best vector"};
	vals = Prepend[Through[{Min, Mean, StandardDeviation}[cost]],
		$GenerationNumber];
	pos = Position[cost, vals[[2]]][[1, 1]];
	vec = pop[[pos]];
	Print @ TableForm[Transpose[{labels, Append[vals, vec]}],
		TableDepth -> 2]
	]

ChooseAB = Compile[{{n, _Integer}},
	Transpose @ Table[
		Module[{a=0, b=0},
		While[a = Random[Integer, {1, n}]; a == pos, Null];
		While[b = Random[Integer, {1, n}]; b == pos || b == a, Null];
		{a, b}],
	{pos, n}]
	]

ChooseABC = Compile[{{n, _Integer}},
	Transpose @ Table[
	Module[{a=0, b=0, c=0},
		While[a = Random[Integer, {1, n}]; a == pos, Null];
		While[b = Random[Integer, {1, n}]; b == pos || b == a, Null];
		While[c = Random[Integer, {1, n}]; c == pos || c == b || c == a, Null];
	{a, b, c}],
	{pos, n}]
	]

(* 1 => from target, 0 => from mate *)
(* "last" gene always comes from mate, even when cr = 0 *)
Genes[n_, d_, cr_] :=
	MapThread[
		RotateRight, (* randomly pick first gene *)
		{
		Transpose @ Append[Table[Ceiling[Random[] - cr], {d-1}, {n}],
			Table[0, {n}]], (* last gene from mate *)
		Table[Random[Integer, {1, d}], {n}] (* rotate by this amount *)
		}]

End[]
EndPackage[]

(* :Examples: *)

(* First example:
costfun[{x_, y_, z_}] := x^2 + y^2 + z^2
xb = 50
pop0 = Table[xb (2 Table[Random[], {3}] - 1), {30}];

DifferentialEvolution[costfun, pop0]

(* to see how this works: *)

{vals, points} = Transpose[$EvolutionHistory];

(* rate of improvement *)
logdiffs = Apply[Subtract, Partition[Log[Min /@ vals], 2, 1], {1}];
ListPlot[logdiffs, PlotJoined -> True, PlotStyle -> Red, PlotRange -> All];

(* population evoloution *)
graphs = Graphics3D[{PointSize[.02],Point /@ #}]& /@ Take[points, 40];
Show[#, PlotRange -> Table[{-xb, xb}, {3}]]& /@ graphs;

*)

(* Second example:

f[{x0_, x1_}] := 100*(x0^2 - x1)^2 + (1 - x0)^2;
xb = 2.048;
pop0 = Table[xb (2 Table[Random[], {2}] - 1), {30}];

DifferentialEvolution[f, pop0, NumberOfGenerations -> 200,
	EvolutionHistory -> True]

(* to see how this works: *)

{vals, points} = Transpose[$EvolutionHistory];
mins = Min /@ vals;

(* rate of improvement *)
logdiffs = Apply[Subtract, Partition[Log[mins], 2, 1], {1}];
ListPlot[logdiffs, PlotJoined -> True, PlotStyle -> Red, PlotRange -> All];

(* population evolution *)

pos = Flatten @ MapThread[First @ Position @ ## &, {vals, mins}];
minpoints = MapThread[#1[[#2]]&, {points, pos}];
graphs = Graphics[{PointSize[.02],Point /@ #}]& /@ points;
minpts = Graphics[{PointSize[.02],Red, Point[#]}]& /@ minpoints;
cp = Graphics[
    ContourPlot[f[{x0, x1}], {x0, -xb, xb}, {x1, -xb, xb},
      ContourShading -> True, Contours ->Prepend[40Range[10], 10],
      PlotPoints -> 45, DisplayFunction -> Identity,
	ColorFunction -> (GrayLevel[.5 + #^(1/4)/2]&)]];
Show[{cp, Graphics[{Line[{{-xb,1},{xb,1}}], Line[{{1,-xb},{1,xb}}]}],#},
	PlotRange -> {{-xb, xb}, {-xb, xb}}, Frame -> True, AspectRatio -> 1,
	DisplayFunction -> $DisplayFunction]& /@ 
		Take[Transpose[{graphs, minpts}], 101];

*)