Minimizing with differential evolution

Question

A differential evolution algorithm is given here. I would like to get this kind of animation. I thought I could use NMinimize, given DifferentialEvolution as an option, but it turns out that does not work as I espected.

Is it possible to extract intermediate step in DifferentialEvolution, or do I have to implement algorithm myself?

f[x_, y_] := 
  -20 E^(-0.2 Sqrt[0.5 (x^2 + y^2)]) - E^(0.5 (Cos[2 π x] + Cos[2 π y])) + E + 20

p1 = 
  Plot3D[f[x, y], {x, -5, 5}, {y, -5, 5}, 
  PerformanceGoal -> "Quality",
  ColorFunction -> "WatermelonColors", 
  Mesh -> None, 
  BoxRatios -> {1, 1, 1}];   

p2 = 
   DensityPlot[f[x, y], {x, -5, 5}, {y, -5, 5}, 
     ColorFunction -> "WatermelonColors", 
     PlotPoints -> 200, 
     PerformanceGoal -> "Quality", 
     Frame -> False, 
     PlotRangePadding -> None];

p3 = Plot3D[0, {x, -5, 5}, {y, -5, 5}, PlotStyle -> Texture[p2], Mesh -> None];

Show[p1, p3, PlotRange -> {0, 15}]

When I use StepMonitor to track iterations as follows, it does not work.

{fit, intermediates} = 
  Reap[NMinimize[{f[x, y], -5 <= x <= 5, -5 <= y <= 5}, {x, y}, 
    MaxIterations -> 1000, 
    Method -> {"DifferentialEvolution", "InitialPoints" -> Tuples[Range[-5, 5], 2]}, 
    StepMonitor :> Sow[{x, y}]]];

Table[
  ListPlot[Take[intermediates[[1, i ;; i + 10]]], 
    Frame -> True, ImageSize -> 350, AspectRatio -> 1], 
  {i, 10, 1000, 100}]

EDIT Here is the result when we used @Michael E2 solution. Cool!!

f[x_, y_] := -20 E^(-0.2 Sqrt[0.5 (x^2 + y^2)]) - 
  E^(0.5 (Cos[2 \[Pi] x] + Cos[2 \[Pi] y])) + E + 20

p1 = Plot3D[f[x, y], {x, -5, 5}, {y, -5, 5}, 
   PerformanceGoal -> "Quality", ColorFunction -> "WatermelonColors", 
   Mesh -> None, BoxRatios -> {1, 1, 1}];

p2 = DensityPlot[f[x, y], {x, -5, 5}, {y, -5, 5}, 
   ColorFunction -> "WatermelonColors", PerformanceGoal -> "Quality", 
   Frame -> False, PlotRangePadding -> None];

p3 = Plot3D[-0.5, {x, -5, 5}, {y, -5, 5}, PlotStyle -> Texture[p2], 
   Mesh -> None];

p4 = Show[p1, p3, PlotRange -> {-0.5, 15}]


  {fit, intermediates} = 
   Reap[NMinimize[{f[x, y], -5 <= x <= 5, -5 <= y <= 5}, {x, y}, 
     MaxIterations -> 40, 
     Method -> {"DifferentialEvolution", 
       "InitialPoints" -> Tuples[Range[-5, 5], 2]}, 
     StepMonitor :> 
      Sow[{Optimization`NMinimizeDump`vecs, 
        Optimization`NMinimizeDump`vals}]]]; // Quiet
sim = Table[
   Show[p4, 
    ListPointPlot3D[{Append[#, 0] & /@ intermediates[[1, i, 1]]}, 
     PlotRange -> {{-5, 5}, {-5, 5}, {-5, 5}}, Boxed -> False, 
     PlotStyle -> Directive[AbsolutePointSize[3], Black]]], {i, 40}];

  Export[NotebookDirectory[] <> "DE2DSim.gif", sim, 
  "AnimationRepetitions" -> 50];

p1 = Plot3D[f[x, y], {x, -5, 5}, {y, -5, 5}, 
   PerformanceGoal -> "Quality", ColorFunction -> "WatermelonColors", 
   Mesh -> None, PlotStyle -> Opacity[0.4], BoxRatios -> {1, 1, 1}];

p2 = DensityPlot[f[x, y], {x, -5, 5}, {y, -5, 5}, 
   ColorFunction -> "WatermelonColors", PerformanceGoal -> "Quality", 
   Frame -> False, PlotRangePadding -> None];

p3 = Plot3D[-0.5, {x, -5, 5}, {y, -5, 5}, PlotStyle -> Texture[p2], 
   Mesh -> None];

p4 = Show[p1, p3, PlotRange -> {-0.5, 15}]


{fit, intermediates} = 
   Reap[NMinimize[{f[x, y], -5 <= x <= 5, -5 <= y <= 5}, {x, y}, 
     MaxIterations -> 40, 
     Method -> {"DifferentialEvolution", 
       "InitialPoints" -> Tuples[Range[-5, 5], 2]}, 
     StepMonitor :> 
      Sow[{Optimization`NMinimizeDump`vecs, 
        Optimization`NMinimizeDump`vals}]]]; // Quiet
sim = Table[
   Show[p4, 
    ListPointPlot3D[{Append[#, 0] & /@ intermediates[[1, i, 1]], 
      Partition[
       Flatten@Riffle[intermediates[[1, i, 1]], 
         f @@@ intermediates[[1, i, 1]]], 3]}, Boxed -> False, 
     PlotStyle -> {Directive[AbsolutePointSize[3], Black], 
       Directive[AbsolutePointSize[3], Red]}]], {i, 40}];

Export[NotebookDirectory[] <> "DE2DSim.gif", sim, 
  "AnimationRepetitions" -> 50];

Multicolumn[
 Table[Show[p4, 
   ListPointPlot3D[{Append[#, 0] & /@ intermediates[[1, i, 1]]}, 
    PlotRange -> {{-5, 5}, {-5, 5}, {-5, 5}}, Boxed -> False, 
    PlotStyle -> Directive[AbsolutePointSize[3], Black]]], {i, 1, 30, 
   2}], 5, Appearance -> "Horizontal"]

Answer 1

Here's a way:

Block[{f},
 f[x_, y_] := -20 E^(-0.2 Sqrt[0.5 (x^2 + y^2)]) - 
   E^(0.5 (Cos[2 \[Pi] x] + Cos[2 \[Pi] y])) + E + 20;
 {fit, intermediates} = 
  Reap[NMinimize[{f[x, y], -5 <= x <= 5, -5 <= y <= 5}, {x, y}, 
    MaxIterations -> 30, 
    Method -> {"DifferentialEvolution", 
      "InitialPoints" -> Tuples[Range[-5, 5], 2]}, 
    StepMonitor :> 
     Sow[{Optimization`NMinimizeDump`vecs, 
       Optimization`NMinimizeDump`vals}]]];
 ]

Manipulate[
 Graphics[{
   PointSize[Medium],
   Point[intermediates[[1, n, 1]], 
    VertexColors -> 
     ColorData["Rainbow"] /@ 
      Rescale[intermediates[[1, n, 2]], 
       MinMax[intermediates[[1, All, 2]]]]]
   },
  PlotRange -> 5, Frame -> True],
 {n, 1, Length@intermediates[[1]], 1}
 ]

You can find out about things like Optimization`NMinimizeDump`vecs by inspecting the code for Optimization`NMinimizeDump`CoreDE.