From acc96ac159c9d1943c864491732dc7ca79e03eba Mon Sep 17 00:00:00 2001
From: Lionel GUEZ <guez@lmd.ens.fr>
Date: Thu, 2 Aug 2018 15:05:42 +0200
Subject: [PATCH] Performance improvement: start looking for outermost contour
 at amplitude min_amp if the difference between innermost level from
 local_extrema and ssh extremum is smaller. This is not only a performance
 improvement but also an improvement of results: the algorithm sometimes finds
 a contour at just min_amp when before it did not find any: the bisection
 algorithm converged just below min_amp. Also, as a result, if an outermost
 contour is found, it is now necessarily of sufficient amplitude.

---
 Documentation_texfol/documentation.tex | 43 +++++---------------------
 get_1_outerm.f                         |  9 +++---
 set_all_outerm.f                       | 14 ++++-----
 3 files changed, 20 insertions(+), 46 deletions(-)

diff --git a/Documentation_texfol/documentation.tex b/Documentation_texfol/documentation.tex
index 3c1df6de..d8bde8d8 100644
--- a/Documentation_texfol/documentation.tex
+++ b/Documentation_texfol/documentation.tex
@@ -964,34 +964,7 @@ vitesse. Mais avec la possibilité d'une amplitude minimale, il faut
 d'abord se préoccuper de tous les contours extérieurs, qui ne sont
 plus indépendants les uns des autres, pour ne pas chercher inutilement
 des contours de maximum de vitesse. Dans le nouvel algorithme, il n'y
-a donc plus de procédure get\_eddy. Cf. tableau (\ref{tab:outermost_contour}).
-\begin{table}[htbp]
-  \centering
-  \begin{tabular}{l|l|l|l|l|l}
-    flat\_extr & \multicolumn{2}{c|}{F} & \multicolumn{3}{c}{V} \\
-    suff\_amp & \multicolumn{2}{c|}{?} & F & \multicolumn{2}{c}{V} \\
-    extr\_map > 0 & \multicolumn{2}{c|}{V} & F & \multicolumn{2}{c}{V} \\
-    noise\_around & \multicolumn{3}{c|}{?} & F & V \\
-    2\ieme & \multicolumn{2}{c|}{V} & \multicolumn{2}{c|}{F} & V \\
-    suff\_amp & F & V & F & \multicolumn{2}{c}{V} \\
-    extr\_map > 0 & F & V & F & \multicolumn{2}{c}{V}
-  \end{tabular}
-  \caption{Logique pour calculer les contours les plus extérieurs dans
-    get\_snapshot. L'ordre des lignes du tableau suit l'ordre des
-    instructions. suff\_amp n'est d'abord défini que sur une partie
-    des extremums et est donc d'abord différent de extr\_map, mais il
-    contient la même information que extr\_map à la fin de la
-    séquence. noise\_around n'a besoin d'être défini que sur une
-    partie des extremums. Le point d'interrogation 
-    signifie que la valeur n'est pas calculée, elle ne nous intéresse
-    pas. La ligne \og 2\ieme\fg{} donne les cas où le programme passe
-    par la deuxième ligne d'appel à get\_1\_outerm.}
-  \label{tab:outermost_contour}
-\end{table}
-Si un extremum A n'est pas plat, on considère que ce n'est pas du
-bruit, il ne peut pas être englobé dans le contour extérieur d'un
-autre extremum B, même s'il s'avère finalement qu'on ne trouve pas de
-contour extérieur autour de A.
+a donc plus de procédure get\_eddy.
 
 Si le domaine considéré n'est pas global, comment traiter un
 tourbillon au bord du domaine ? La taille du contour le plus extérieur
@@ -1070,8 +1043,9 @@ flow(j)\%list\_vis\%delta\_in, flow\%number\_eddies. Cf. algorithme
 
     \STATE i1\_lat = flow(j - delta)\%ind\_extr(2, i1)
 
-    \FOR{i2 > 0 dans abs(flow(j)\%extr\_map(i1\_lon - dist\_lim:i1\_lon
-      + dist\_lim, i1\_lat - dist\_lim:i1\_lat + dist\_lim))}
+    \FOR{i2 $\ne$ 0 dans abs(flow(j)\%extr\_map(i1\_lon -
+      dist\_lim:i1\_lon + dist\_lim, i1\_lat - dist\_lim:i1\_lat +
+      dist\_lim))}
 
     \IF{flow(j)\%list\_vis(i2)\%suff\_amp et (flow(j -
       delta)\%list\_vis(i1)\%delta\_out $\ge$ delta ou
@@ -1138,11 +1112,10 @@ en utilisant le fait que deux maximums ne peuvent être côte à côte et
 deux minimums ne peuvent être côte à côte.
 
 J'avais d'abord codé le sens de l'extremum (maximum ou minimum) par le
-signe dans extr\_map. Mais j'ai besoin de coder dans une carte le fait
-que l'amplitude soit significative ou non. J'ajoute donc un argument
-vecteur cyclone à local\_extrema, qui prend moins de place qu'une
-carte supplémentaire. Pour le sens de l'extremum, je n'ai pas besoin
-que l'information soit sur une carte.
+signe dans extr\_map. Mais c'est plus clair et plus simple sans ce
+codage. J'ajoute donc un argument vecteur cyclone à local\_extrema,
+qui prend moins de place qu'une carte supplémentaire. Pour le sens de
+l'extremum, je n'ai pas besoin que l'information soit sur une carte.
 
 Cf. aussi fontion \verb+skimage.feature.peak_local_max+ dans le paquet
 Python scikit-image. La fonction applique
diff --git a/get_1_outerm.f b/get_1_outerm.f
index 3082da67..811dc666 100644
--- a/get_1_outerm.f
+++ b/get_1_outerm.f
@@ -11,7 +11,8 @@ contains
     ! dichotomy on level of ssh. If the procedure cannot find an
     ! outermost contour, it return a null ssh_contour. An outermost
     ! contour is not found if, and only if, there is no good contour
-    ! at innermost level.
+    ! at innermost level. If a contour is found than its level is
+    ! farther from the extremum than innermost level.
 
     use contour_531, only: polyline
     use derived_types, only: ssh_contour, missing_ssh
@@ -23,9 +24,9 @@ contains
     real, intent(in):: coord_extr(:) ! (2)
 
     real, intent(in):: innermost_level
-    ! ssh level of the innermost contour around the target extremum,
-    ! in m. By construction, innermost_level < extremum for a maximum,
-    ! > extremum for a minimum.
+    ! ssh level of the innermost contour to consider around the target
+    ! extremum, in m. Assume: innermost_level < extremum for a
+    ! maximum, > extremum for a minimum.
 
     real, intent(in):: outside_points(:, :) ! (2, :) longitude and
     ! latitude, in rad, of all the significant extrema, except the
diff --git a/set_all_outerm.f b/set_all_outerm.f
index 8b633762..0c20225e 100644
--- a/set_all_outerm.f
+++ b/set_all_outerm.f
@@ -82,7 +82,7 @@ contains
     selection = argwhere(cyclone)
     n_cycl = size(selection)
     selection = selection(indexx(s%list_vis(selection)%ssh_extr))
-    sorted_extr(:n_cycl) = selection(n_cycl:1:- 1) ! (descending order of ssh)
+    sorted_extr(:n_cycl) = selection(n_cycl:1:- 1) ! descending order of ssh
 
     selection = argwhere(.not. cyclone)
     selection = selection(indexx(s%list_vis(selection)%ssh_extr))
@@ -95,17 +95,17 @@ contains
        llc = max(s%ind_extr(:, i) - max_radius, 1)
        urc = min(s%ind_extr(:, i) + max_radius, [nlon, nlat])
 
+       ! No need to consider contours with amplitudes < min_amp:
+       if (abs(s%list_vis(i)%ssh_extr - innermost_level(i)) < min_amp) &
+            innermost_level(i) = s%list_vis(i)%ssh_extr &
+            + merge(min_amp, - min_amp, s%list_vis(i)%cyclone)
+
        s%list_vis(i)%out_cont = get_1_outerm(s%list_vis(i)%cyclone, &
             s%list_vis(i)%coord_extr, innermost_level(i), &
             nearby_extr(s%extr_map(llc(1):urc(1), llc(2):urc(2)), s%list_vis, &
             i), ssh(llc(1):urc(1), llc(2):urc(2)), &
             corner = corner + (llc - 1) * step, step = step)
-       if (s%list_vis(i)%out_cont%n_points == 0) then
-          s%list_vis(i)%suff_amp = .false.
-       else
-          s%list_vis(i)%suff_amp = abs(s%list_vis(i)%out_cont%ssh &
-               - s%list_vis(i)%ssh_extr) >= min_amp
-       end if
+       s%list_vis(i)%suff_amp = s%list_vis(i)%out_cont%n_points /= 0
     end do
 
   end subroutine set_all_outerm
-- 
GitLab