[Wolfram Mathematica] Карта минимумов ландшафта

[Mysterious Light]

Есть такое построение:

пусть дана красивая (гладкая, без выраждений корней любых производных) действительнозначная функция $a$ действительных аргументов, которая мыслится как поверхность (график) $a(\mathbf{x})$ в пространстве $(\mathbf{x},a)$
пусть $M = \{ \mathbf{x} | \operatorname{d}a(\mathbf{x}) = 0, \operatorname{d}^2a(\mathbf{x}) > 0 \}$ — множество всех локальных минимумов функции,
пусть для каждого $\mathbf{x}_0\in M$ определён колодец $W(\mathbf{x}_0) = \{ \mathbf{x} | C(\mathbf{x},\mathbf{x}_0) \}$,
где $C(x,y)$ означает существование непрерывной монотонноубывающей кривой на поверхности $a$ между $x$ и $y$, то есть найдётся такая непрерывная $c: [0,1] \to \operatorname{dom} a$, что $c(0)=x$, $c(1)=y$ и $a\circ c$ монотонно убывает.
Тогда все колодцы $\{ W(\mathbf{x}) | \mathbf{x}\in M \}$ образуют покрытие $\operatorname{dom}a$ с точностью до их границ. В каждом колодце находится ровно один локальный минимум.
Между тем, их границы предсталяются кривыми на одну размерность меньше $a$, на которых можно поискать их локальные минимумы $J=\{\mathbf{x} | \min_{\partial W(\matbf{y})} a(\mathbf{x}), \exists \mathbf{y}\in M \}$. Очевидно, это седловые точки $\operatorname{d} a = 0$, и в касательном(-ых) направлении(-ях) $\operatorname{d}^2a>0$. Однако, будучи пограничной точкой между двумя колодцами, в соответствующем направлении $\operatorname{d}^2a<0$.
Кстати, пару $\langle M, J \rangle$ можно рассматривать как мультиграф с вершинами $M$ и рёбрами $J$, где каждое ребро — пограничная точка, соединяющая соответствующие два колодцаю

Решил взяться за задачу построить такое разбиение.
В приложенном файле найдёте постраение графа. Однако, иногда оно даёт неправильный результат. Буду очень благодарен, если кто-то найдёт причину (а мне кажется, что причина где-то в методе работы FindMinimum) и ещё скажет, как можно по-красивее или по-проще переписать программу.
Ещё хотелось бы построить $\bigcup \{ \partial W(\mathbf{x}) \| \mathbf{x}\in M \}$ в виде контура. Как это сделать — не представляю.

(Оффтоп)

Что-то я не нашел, как прикреплять файл к сообщению.
Выкладываю на файлообменник: http://rghost.ru/60196002

(сам код)

Generate a new landscape function and limit it in the square [0..1]*[0..1] with (k Exp[\[Alpha]/x])-like functions

\[Alpha] = 1/10;
k = 50;
altitude[x_, y_] = 
  k (Exp[\[Alpha]/(x + \[Alpha])] + Exp[\[Alpha]/(y + \[Alpha])] + 
      Exp[\[Alpha]/(1 - x + \[Alpha])] + Exp[\[Alpha]/(1 - y + \[Alpha])] - 
      4) + Sum[RandomInteger[{-100, 100}]/(i + j)
      Sin[\[Pi] i x] Sin[\[Pi] j y], {i, 1, 5}, {j, 1, 5}];

dplot = DensityPlot[altitude[x, y], {x, 0, 1}, {y, 0, 1}, PlotRange -> All, 
  ColorFunction -> (Blend[{Darker[Green], Green, Lighter[Green], Yellow, 
       Orange, Darker[Red]}, #] &), MaxRecursion -> 3, PlotPoints -> 20]
Plot3D[altitude[x, y], {x, 0, 1}, {y, 0, 1}, PlotRange -> All]

Find numerically roots of f and filter distinct solutions.

FindDistinctRoots[f_, {x_, x1_, x2_, dx_}, {y_, y1_, y2_, dy_}] := 
  DeleteDuplicates[
   Select[
    Join @@ Table[
      {x, y} /. FindRoot @@ {f, {{x, x0}, {y, y0}}},
      {x0, x1, x2, dx}, {y0, y1, y2, dy}],
    x1 <= #[[1]] <= x2 && y1 <= #[[2]] <= y2 &],
   4 Abs[#1[[1]] - #2[[1]]] < dx && 4 Abs[#1[[2]] - #2[[2]]] < dy &];

Solve \[Del]altitude=0

grad[x_, y_] = 
  If[-\[Alpha]/2 < x < 1 + \[Alpha]/2 && -\[Alpha]/2 < y < 1 + \[Alpha]/2, 
   Evaluate[D[altitude[x, y], {{x, y}}]], {x - 0.5, y - 0.5}];
suspect = FindDistinctRoots[grad[x, y], {x, 0, 1, 0.1}, {y, 0, 1, 0.1}] // 
  Quiet

For each point, find a matrix G of all second derivatives, diagonalized them and find all junctions which are points with \[Del]=0 and G contains one negative eigenvalue and the rest positive.

G[f_, {x0_, y0_}] := 
  Apply[Derivative[##][f][x0, 
     y0] &, {{{2, 0}, {1, 1}}, {{1, 1}, {0, 2}}}, {2}];
Diagonalization[g_] := {Eigenvalues[g], 
   Inverse[Transpose[Normalize /@ Eigenvectors[g]]]};
JunctionDirection[g_] := 
  With[{js = Last /@ Select[Transpose[Diagonalization[g]], First[#] < 0 &]},
   If[Length[js] == 1, js[[1]], Message[JunctionDirection::notAJunction]; None]
   ];

Off[JunctionDirection::notAJunction];
juncDirPos = 
  DeleteCases[
   Transpose[{JunctionDirection /@ (G[altitude, #] & /@ suspect), 
     suspect}], {None, _}];
Prepend[juncDirPos, {"direction", "position"}] // Grid

minima = Select[suspect, 
  And @@ (# > 0 & /@ First[Diagonalization[G[altitude, #]]]) &]

Show[
 dplot,
 Graphics[{PointSize[Large], Blue, Point[minima], Black, Thick,(*Point[Last/@
   juncDirPos],*)Line[{#2 - #1/20, #2 + #1/20}] & @@@ juncDirPos}]
 ]

NearestMinimum[{x0_, y0_}] := 
  Nearest[minima, 
    FindArgMin[altitude[x, y], {{x, x0}, {y, y0}}, Method -> "Newton"]][[1]];

NearestMinimum[{x0_, y0_}, {dx_, dy_}] := 
  Nearest[minima, 
    FindArgMin[altitude[x, y], {{x, x0, x0 + dx}, {y, y0, y0 + dy}}]][[1]];

Hash-function of point is a sum of its coordinates. There might not be collitions.
Vertices are minima, edges are junctions.

vertices = Total /@ minima
edges = Cases[
  juncDirPos, {dir_, pos_} :> 
   Tooltip[Sort[
     Total[NearestMinimum[pos - dir 10^-3, dir 10^-4]] \[UndirectedEdge] 
      Total[NearestMinimum[pos + dir 10^-3, dir 10^-4]]], pos]]
Graph[vertices, 
 DeleteDuplicates[DeleteCases[edges, Tooltip[x_ \[UndirectedEdge] x_, _]], 
  First[#1] == First[#2] &], VertexCoordinates -> minima]
Show[dplot, %, PlotRange -> {{0, 1}, {0, 1}}]
Show[
 dplot,
 Graphics[{
   PointSize[Large], Blue, Point[minima], Black, Thick,(*Point[Last/@
   juncDirPos],*)Line[{#2 - #1/20, #2 + #1/20}] & @@@ juncDirPos,
   Thin, Gray,
   Cases[juncDirPos, {dir_, pos_} :> 
     Line[{NearestMinimum[pos - dir 10^-3, dir 10^-4], pos, 
       NearestMinimum[pos + dir 10^-3, dir 10^-4]}]]
   (*Cases[%%%,Tooltip[m1_\[UndirectedEdge]m2_,pos_]:>Line[{m1,pos,m2}]]*)
   }
  ]]

[Mysterious Light]

Поправка ко вступительному посту: колодец локального минимума состоит не из точек, до которых можно прийти из локального минимума монотонной непрерывной функцией, а из точек, от которых можно спуститься наискорейшим (градиентным) спуском именно к данному минимуму.

Созрел конкретный вопрос:
Как в Wolfram Mathematica сделать градиентный спуск? FindMinimum[... Method -> "ConjugateGradient"] не работает правильно. Экспериментально было подобрано, что "QuasiNewton" в большем числе случаев находит правильный минимум.