(* http://www.hrz.uni-dortmund.de/A1/software/mathematica.html *)
(* Mathematica available on HRZ server zx3 *)

<<Statistics`Master`
<<Graphics`Master`

(* function to be optimized *)
f[x_] := ( ndim = Dimensions[x][[1]]; 
           (1 + Sum[(x[[i]]-i)^2, {i, 1, ndim}])^-1 )

(* produces a vector of N(0,1) random variates *)
norm = NormalDistribution[0, 1]
GAUSS[x_, s_] := ( ndim = Dimensions[x][[1]];
                   s*Table[Random[norm], {i, 1, ndim}] )

(* ************************************************************************** *)

OnePlusOneES[individual_, genStop_, sigmaStop_, a_:1.2, lEvo_:50, pOpt_:1/5] := 
   ( objParent = individual[[2]];
     parentFitness = individual[[1]];
     sigma = individual[[3]];
 (*  gen = 0;  *)                                                 (* del *)
     While[ True, 
       success = 0;
       Do[ objOffspring = objParent + GAUSS[objParent, sigma];
           offspringFitness = f[objOffspring];
           fitnlist = Append[fitnlist, parentFitness];            (* del *)
           sigmalist = Append[sigmalist, sigma];                  (* del *)
 (*        gen = gen + 1;                               *)        (* del *)
           If[ offspringFitness >= parentFitness,       
               objParent = objOffspring;               
               parentFitness = offspringFitness;      
 (*            Print["gen = ", gen,                     *)        (* del *)
 (*                  "   Fitness = ", parentFitness];   *)        (* del *)
               success = success + 1]  
           , {lEvo} ];
       pSuccess = success/lEvo;
       If[pSuccess > pOpt, sigma = sigma*a];
       If[pSuccess < pOpt, sigma = sigma/a];
 (*    Print["   "];                                           *) (* del *)
 (*    Print["sigma = ", sigma, "   P-success = ", pSuccess];  *) (* del *)
       If[ gen > genStop || sigma < sigmaStop, Break[] ] ];
       {parentFitness, objParent, sigma} )

(* ************************************************************************** *)
(* ************************************************************************** *)

(* Experiment 1 *)
(* Initialization *)

sigmaInitial = 0.1
objInitial = {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}
parentFitness = f[objInitial]
individual = {parentFitness, objInitial, sigmaInitial}

fitnlist = {parentFitness}
sigmalist = {sigmaInitial}

xx = OnePlusOneES[individual, 4000, 0.0001]

pf1=LogListPlot[1-fitnlist, AxesLabel -> {"g", "1-f"}]
ps1=LogListPlot[sigmalist, AxesLabel -> {"g", "sigma"}]

Show[pf1, ps1]

(* ************************************************************************** *)

(* Experiment 2: changing a to a=1.0 *)
(* Initialization *)

sigmaInitial = 0.1
objInitial = {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}
parentFitness = f[objInitial]
individual = {parentFitness, objInitial, sigmaInitial}
fitnlist = {parentFitness}
sigmalist = {sigmaInitial}

xx = OnePlusOneES[individual, 4000, 0.0001, 1.0]

pf2=LogListPlot[1-fitnlist, AxesLabel -> {"g", "1-f"}]
ps2=LogListPlot[sigmalist, AxesLabel -> {"g", "sigma"}]

Show[pf1,ps1, pf2, ps2]

(* ************************************************************************** *)