diff --git a/Documentation_texfol/.gitignore b/Documentation_texfol/.gitignore
new file mode 100644
index 0000000000000000000000000000000000000000..7604ca46acf0379031a19fcf0ec0f4bc4348ab7c
--- /dev/null
+++ b/Documentation_texfol/.gitignore
@@ -0,0 +1,7 @@
+/*.aux
+/*.log
+/*.out
+/*.pdf
+/*.synctex*
+/*.toc
+/auto/
\ No newline at end of file
diff --git a/Documentation_texfol/Graphiques/.gitignore b/Documentation_texfol/Graphiques/.gitignore
new file mode 100644
index 0000000000000000000000000000000000000000..d0f356253e7b577ba9b540ca650a621460cca87d
--- /dev/null
+++ b/Documentation_texfol/Graphiques/.gitignore
@@ -0,0 +1,7 @@
+call_graph.pdf
+degeneracy.pdf
+input_output.pdf
+m3.pdf
+processes.pdf
+regions.pdf
+window.pdf
diff --git a/Documentation_texfol/Graphiques/GNUmakefile b/Documentation_texfol/Graphiques/GNUmakefile
new file mode 100644
index 0000000000000000000000000000000000000000..b6cade49662a5ca170a2d19560289b2cc158793c
--- /dev/null
+++ b/Documentation_texfol/Graphiques/GNUmakefile
@@ -0,0 +1,14 @@
+sources = window processes m3 call_graph input_output degeneracy
+
+objects := $(addsuffix .pdf, ${sources})
+
+%.pdf: %.gv
+ dot -Tpdf -o $@ $<
+
+%.pdf: %.odg
+ unoconv --doctype=graphics $<
+
+all: ${objects}
+
+clean:
+ rm -f ${objects}
diff --git a/Documentation_texfol/Graphiques/call_graph.gv b/Documentation_texfol/Graphiques/call_graph.gv
new file mode 100644
index 0000000000000000000000000000000000000000..0689b3c6bd6e621e6ba313a3aabb2d1e7b500f71
--- /dev/null
+++ b/Documentation_texfol/Graphiques/call_graph.gv
@@ -0,0 +1,9 @@
+digraph call_graph
+{
+ main -> {get_snapshot successive_overlap non_successive_overlap};
+ main -> dispatch_snapshot;
+ non_successive_overlap -> {interpolate_eddy weight};
+ dispatch_snapshot -> write_eddy;
+ interpolate_eddy -> write_eddy;
+ successive_overlap -> weight;
+}
\ No newline at end of file
diff --git a/Documentation_texfol/Graphiques/comparaison_201.pdf b/Documentation_texfol/Graphiques/comparaison_201.pdf
new file mode 100644
index 0000000000000000000000000000000000000000..af99153ff1bfeaed1fd2d53dc91da512f6ccea14
Binary files /dev/null and b/Documentation_texfol/Graphiques/comparaison_201.pdf differ
diff --git a/Documentation_texfol/Graphiques/comparaison_202.pdf b/Documentation_texfol/Graphiques/comparaison_202.pdf
new file mode 100644
index 0000000000000000000000000000000000000000..117a887a8c2b4b5c556b06155d785d85b2ab1d13
Binary files /dev/null and b/Documentation_texfol/Graphiques/comparaison_202.pdf differ
diff --git a/Documentation_texfol/Graphiques/degeneracy.odg b/Documentation_texfol/Graphiques/degeneracy.odg
new file mode 100644
index 0000000000000000000000000000000000000000..67bfc175c74d929776111358f4934e6c7152708d
Binary files /dev/null and b/Documentation_texfol/Graphiques/degeneracy.odg differ
diff --git a/Documentation_texfol/Graphiques/extrema.pdf b/Documentation_texfol/Graphiques/extrema.pdf
new file mode 100644
index 0000000000000000000000000000000000000000..11c52d384a8d315484765e67eb7790c97868a922
Binary files /dev/null and b/Documentation_texfol/Graphiques/extrema.pdf differ
diff --git a/Documentation_texfol/Graphiques/input_output.odg b/Documentation_texfol/Graphiques/input_output.odg
new file mode 100644
index 0000000000000000000000000000000000000000..001bc6e7995dea4372b2484bcac344ed828e94ed
Binary files /dev/null and b/Documentation_texfol/Graphiques/input_output.odg differ
diff --git a/Documentation_texfol/Graphiques/m3.odg b/Documentation_texfol/Graphiques/m3.odg
new file mode 100644
index 0000000000000000000000000000000000000000..3d45c1ea9f203661100066bbb9bfc517068fc70c
Binary files /dev/null and b/Documentation_texfol/Graphiques/m3.odg differ
diff --git a/Documentation_texfol/Graphiques/processes.odg b/Documentation_texfol/Graphiques/processes.odg
new file mode 100644
index 0000000000000000000000000000000000000000..29bf08c3978e29f1e86ccbd41daebf7f8ee05f52
Binary files /dev/null and b/Documentation_texfol/Graphiques/processes.odg differ
diff --git a/Documentation_texfol/Graphiques/regions.py b/Documentation_texfol/Graphiques/regions.py
new file mode 100755
index 0000000000000000000000000000000000000000..223d196f70e6b008063abaff6de19553b5375e7b
--- /dev/null
+++ b/Documentation_texfol/Graphiques/regions.py
@@ -0,0 +1,27 @@
+#!/usr/bin/env python3
+
+import netCDF4
+from mpl_toolkits import basemap
+import jumble
+from matplotlib import pyplot as plt
+import glob
+import itertools
+
+m = basemap.Basemap(llcrnrlon=-30, llcrnrlat=-90, urcrnrlon=60, urcrnrlat = 0)
+plt.figure()
+m.drawcoastlines()
+m.drawmeridians(range(0,361, 15), labels = [0, 0, 0, 1])
+m.drawparallels(range(-90,91, 15), labels = [1, 0, 0, 0])
+colors = itertools.cycle(['blue', 'orange', 'green', 'red', 'purple', 'brown',
+ 'pink', 'gray', 'olive', 'cyan'])
+
+for filename in glob.glob("h_region_*.nc") + ["h_outermost.nc"]:
+ with netCDF4.Dataset(filename) as f:
+ jumble.draw_bbox(f.variables["lon"][0], f.variables["lon"][-1],
+ f.variables["lat"][0], f.variables["lat"][-1],
+ colors = [next(colors)], label = filename)
+
+plt.legend(loc = "lower left", bbox_to_anchor = (1, 0.5))
+plt.subplots_adjust(right = 0.7)
+plt.savefig("regions.pdf")
+##plt.show()
diff --git a/Documentation_texfol/Graphiques/test_get_eddy.pdf b/Documentation_texfol/Graphiques/test_get_eddy.pdf
new file mode 100644
index 0000000000000000000000000000000000000000..be29ff9ad0c4aab965b3790c5a86acd26d5a4129
Binary files /dev/null and b/Documentation_texfol/Graphiques/test_get_eddy.pdf differ
diff --git a/Documentation_texfol/Graphiques/test_get_snapshot.pdf b/Documentation_texfol/Graphiques/test_get_snapshot.pdf
new file mode 100644
index 0000000000000000000000000000000000000000..73830a06361538153a4be4f6f5b7d026152f143e
Binary files /dev/null and b/Documentation_texfol/Graphiques/test_get_snapshot.pdf differ
diff --git a/Documentation_texfol/Graphiques/window.odg b/Documentation_texfol/Graphiques/window.odg
new file mode 100644
index 0000000000000000000000000000000000000000..61ff61792cbf4325f84c1494feb0b7c838e2dea7
Binary files /dev/null and b/Documentation_texfol/Graphiques/window.odg differ
diff --git a/Documentation_texfol/documentation.tex b/Documentation_texfol/documentation.tex
new file mode 100644
index 0000000000000000000000000000000000000000..884bb9c7857d4d972bedc2b559643ae8cb574c84
--- /dev/null
+++ b/Documentation_texfol/documentation.tex
@@ -0,0 +1,1332 @@
+\documentclass[a4paper,french]{article}
+
+\usepackage[utf8x]{inputenc}
+
+\usepackage{amsmath}
+
+\usepackage[T1]{fontenc}
+\usepackage{lmodern}
+
+\usepackage{algorithmic}
+\usepackage{babel}
+\usepackage{graphicx}
+\usepackage{xcolor}
+
+\usepackage[np]{numprint}
+
+\usepackage{hyperref} \hypersetup{pdftitle={Documentation du programme
+ en Fortran de détection de tourbillons}, pdfauthor={Lionel Guez},
+ hypertexnames=false}
+
+\usepackage{algorithm}
+
+\renewcommand{\algorithmicfor}{\textbf{pour}}
+\renewcommand{\algorithmicend}{\textbf{fin}}
+\renewcommand{\algorithmicdo}{\textbf{faire}}
+\renewcommand{\algorithmicto}{\textbf{à}}
+\renewcommand{\algorithmicif}{\textbf{si}}
+\renewcommand{\algorithmicthen}{\textbf{alors}}
+\renewcommand{\algorithmicelse}{\textbf{sinon}}
+\renewcommand{\algorithmiccomment}[1]{\textcolor{blue}{(*~#1~*)}}
+\renewcommand{\algorithmicwhile}{\textbf{tant que}}
+\renewcommand{\algorithmicrepeat}{\textbf{répéter}}
+\renewcommand{\algorithmicuntil}{\textbf{jusqu'à}}
+
+ %---------------
+
+\graphicspath{{Graphiques/}}
+
+\author{Lionel GUEZ}
+
+\title{Documentation du programme en Fortran de détection de
+ tourbillons}
+
+ %---------------
+
+\begin{document}
+
+\maketitle
+\tableofcontents
+
+\section{Divers}
+
+Pour le critère de fusion ou de division des tourbillons, utiliser la
+vitesse géostrophique. Chaque ligne de niveau de la hauteur correspond
+à une ligne de courant de la vitesse. La norme de la vitesse et la
+composante orthoradiale doivent être à peu près constantes sur une
+ligne, prendre la moyenne sur la ligne de l'une ou l'autre. Chercher
+la ligne de vitesse maximale en partant de l'extremum et en
+s'éloignant. La zone intérieure à cette ligne est le c\oe{}ur du
+tourbillon. Les extremums à deux dates correspondent au même
+tourbillon si les c\oe{}urs se recouvrent.
+
+But du programme : produire une base de données de tourbillons et un
+graphe. Les sommets du graphe sont les tourbillons. Les arcs relient
+les tourbillons dont les c\oe{}urs se recouvrent à deux dates
+successives.
+
+Le graphe est orienté acyclique. L'ordre topologique est tout trouvé :
+numérotation par dates successives. Le graphe non orienté associé est
+une forêt.
+
+Pour comparer les tourbillons à deux dates, il semble coûteux de
+n'utiliser que les listes de tourbillons aux deux dates. Cela implique
+une double boucle sur les listes de tourbillons, à l'intérieur d'une
+boucle sur les dates, soit potentiellement
+$8000 \times 8000 \times 5000$ comparaisons. Pour éviter cette
+comparaison aveugle, on peut conserver l'information sur la position
+des tourbillons dans des cartes aux différentes dates. La conservation
+en mémoire vive à 5000 dates semble aussi trop lourde. Il faut se
+restreindre aux cinq (max\_delta + 1) dates nécessaires à la
+profondeur de vue.
+
+Lorsqu'on compare des tourbillons à distance temporelle delta, on veut
+tester s'ils n'ont pas de prédécesseur ou successeur visible à une
+date intermédiaire. Il faut donc, si on crée un arc entre deux dates
+successives, ou delta arcs à distance delta (par interpolation de
+tourbillons invisibles), enregistrer ce fait en prévision du passage
+au delta supérieur. Mais d'un autre côté, on ne veut pas que le fait
+de mettre des arcs à l'étape delta ait une importance pour la suite de
+l'étape delta. (Et accessoirement, on ne veut pas que l'ordre dans
+lequel on passe les tourbillons à l'étape delta ait une importance.)
+On est donc obligé d'enregistrer, pour chaque tourbillon, non
+seulement s'il a des prédécesseurs visibles et s'il a des successeurs
+visibles, mais aussi à quelle distance temporelle.
+
+On ne cherche des superpositions qu'entres tourbillons visibles, en
+excluant les tourbillons interpolés.
+
+Possibilité de parallélisme dans l'extraction des tourbillons mais au
+prix d'une conservation de ssh, uv, v en mémoire vive à plusieurs
+dates.
+
+On pourrait écrire l'algorithme principal en explicitant les
+get\_snpashot, en explicitant le cas m < n\_proc et en supprimant la
+variable k\_end\_main\_loop. Ce serait mieux du point de vue de la
+performance, mais je pense que ce serait finalement moins clair.
+
+Cf. figure (\ref{fig:call_graph}).
+\begin{figure}[htbp]
+ \centering
+ \includegraphics[width=\textwidth]{call_graph}
+ \caption{Arbre des appels du programme de détection de tourbillons.}
+ \label{fig:call_graph}
+\end{figure}
+
+Cas où les composantes out\_cont ou speed\_cont d'un tourbillon sont
+vides, c'est-à-dire égales à null\_ssh\_contour(). On peut distinguer
+les cas du tableau (\ref{tab:null_ssh_contour}).
+\begin{table}[htbp]
+ \centering
+ \begin{tabular}{lllllll|l}
+ out\_cont & suff\_amp & extr\_map & radius4 & speed\_cont
+ & max\_speed & fraction & note \\
+ \hline
+ null & F & < 0 & 0 & null & \np{e4} & 5 \% & 1 \\
+ $\ne$ null & F & < 0 & 0 & null & \np{e4} & 6 \% \\
+ $\ne$ null & T & > 0 & 1 & null & NaN & 3 \% & 2 \\
+ $\ne$ null & T & > 0 & 1 & null & $\ne \np{e4}$, $\ne$ NaN & 47 \% \\
+ $\ne$ null & T & > 0 & $\ge 2$ & null & $\ne \np{e4}$, $\ne$ NaN & 21 \%
+ & 3
+ \end{tabular}
+ \caption{Cas où les composantes out\_cont ou speed\_cont d'un tourbillon sont
+ vides. Colonne fraction : pour la région 5, nombre
+ de cas sur 260 extremums au total. Note 1 : pas de bon contour à
+ innermost\_level. Note 2 : out\_cont est près d'une côte, champ de
+ vitesse non défini, max\_speed vaut $10^4$ dans le fichier
+ DBF. Note 3 : pas mieux que out\_cont pour la vitesse maximale.}
+ \label{tab:null_ssh_contour}
+\end{table}
+Sur la région 5, sur 260 extremums, je trouve les nombres de cas
+suivants pour les 5 lignes du tableau (\ref{tab:null_ssh_contour}) :
+12, 15, 7, 122, 54. Au total : 210 cas où les composantes out\_cont ou
+speed\_cont d'un tourbillon sont vides, sur 260 tourbillons, soit 81
+\% des tourbillons.
+
+\section{Calcul de la vitesse azimutale par rapport à l'extremum}
+
+Notons $C$ la position de l'extremum et $M$ un point où l'on veut
+calculer la vitesse azimutale, $\vec k_C$ et $\vec k$ les vecteurs
+verticaux en $C$ et $M$. Le vecteur unitaire azimutal en $M$ par
+rapport à $C$ est :
+\begin{equation*}
+ \vec u = \frac{1}{\| \vec k_C \wedge \vec k \|} \vec k_C \wedge \vec k
+\end{equation*}
+C'est un vecteur orthogonal au grand cercle passant par $C$ et
+$M$. La composante azimutale de la vitesse $\vec V$ est donc
+$\vec V \cdot \vec u$. En pratique, on peut la calculer en projetant
+les trois vecteurs sur la base $(\mathbf{x_2}, \mathbf{y_2},
+\mathbf{z_p})$ (cf. résumé du chapitre 1 de Holton (2004 k0727)) :
+\begin{equation*}
+ \vec V =
+ \begin{bmatrix}
+ - u \sin \lambda - v \sin \phi \cos \lambda \\
+ u \cos \lambda - v \sin \phi \sin \lambda \\
+ v \cos \phi
+ \end{bmatrix},
+ \vec k_C =
+ \begin{bmatrix}
+ \cos \phi_C \cos \lambda_C \\
+ \cos \phi_C \sin \lambda_C \\
+ \sin \phi_C
+ \end{bmatrix},
+ \vec k =
+ \begin{bmatrix}
+ \cos \phi \cos \lambda \\
+ \cos \phi \sin \lambda \\
+ \sin \phi
+ \end{bmatrix}
+\end{equation*}
+En résumé, il faut calculer un produit vectoriel, normaliser le
+résultat et calculer un produit scalaire. (La permutation circulaire
+du produit mixte ne simplifie rien : il faut toujours calculer le
+produit vectoriel pour la normalisation.) Un inconvénient de cette
+méthode est que, si $C$ et $M$ sont proches,
+$\| \vec k_C \wedge \vec k \| \ll 1$. C'est-à-dire que le produit
+vectoriel doit être calculé par différence petite par rapport aux
+termes de la soustraction, donc avec une mauvaise précision.
+
+Si $C$ et $M$ sont proches, une méthode peut-être plus précise
+consisterait à considérer la projection stéréographique oblique de centre
+$C$. On peut projeter $M$ et projeter $\vec V$, calculer la composante
+azimutale du vecteur projeté dans le plan de projection, et inverser
+la projection à partir de cette composante azimutale. Les calculs
+semblent assez laborieux. La projection de la vitesse fait intervenir
+le jacobien de la projection stéréographique oblique. Je n'ai
+d'ailleurs pas trouvé les formules toutes prêtes.
+
+Dans l'approximation plane, on confond la sphère au voisinage de $C$
+avec le plan tangent en $C$ (le plan de la la projection
+stéréographique oblique de centre $C$). Les formules de la projection
+stéréographique oblique se réduisent, à l'ordre 1 en
+$\delta \phi = \phi - \phi_C$ et
+$\delta \lambda = \lambda - \lambda_C$, à :
+\begin{align*}
+ & x \approx a \cos \phi_C\ \delta \lambda \\
+ & y \approx a\ \delta \phi
+\end{align*}
+Admettons que la vitesse soit identique à sa projection, à ce niveau
+d'approximation. Alors, en notant $\rho$ la
+distance entre l'extremum et $M$, la vitesse azimutale est :
+\begin{equation*}
+ V_\theta = \frac{x v - y u}{\rho}
+\end{equation*}
+
+
+\section{Bons contours}
+
+Soit un extremum de SSH. Un \og bon contour\fg{} de cet extremum, est
+une ligne de niveau de SSH fermée contenant cet extremum et ne
+contenant aucun autre extremum de SSH d'amplitude significative.
+
+Supposons que l'amplitude minimale soit nulle.
+
+Je crois qu'il ne peut exister qu'un bon contour au plus à un niveau
+donné. Je crois aussi, sans l'avoir démontré, que tous les points du
+maillage intérieurs au bon contour ont une SSH comprise entre la SSH
+de l'extremum et la SSH du contour.
+
+Appelons \og chemin de sortie\fg{} un des quatre chemins partant de
+l'extremum, sur le maillage, en ligne droite à longitude constante ou
+à latitude constante (sans rebroussement, sans sauter de point ni
+rester sur le même point), allant jusqu'au premier point du maillage
+strictement extérieur au contour. Je crois, sans l'avoir démontré,
+que, pour chacun des quatre chemins de sortie, la suite des valeurs de
+SSH aux points rencontrés à l'intérieur du contour (tous les points
+sauf le dernier) varie de façon monotone. J'utilise cette hypothèse
+dans outermost\_possible\_level.
+
+Admettons les propriétés ci-dessus. Considérons l'un des quatre
+chemins de sortie. Indiçons, mettons avec $l$ variant de 1 (pour
+l'extremum) à $n$ (pour le dernier point du chemin, le seul
+strictement à l'extérieur), les points rencontrés sur ce
+chemin. Notons $h_c$ la valeur de SSH sur le contour et $h_l$ la
+valeur de SSH en un point $l$. Supposons :
+\begin{equation}
+ \label{eq:non_degeneracy}
+ \forall l < n, h_l \ne h_{l + 1}
+\end{equation}
+Soit $l < n$ (un point quelconque du chemin sauf le dernier),
+c'est-à-dire que nous supposons le point $l$ à l'intérieur du
+contour. Si par exemple l'extremum est un maximum alors :
+\begin{equation*}
+ h_l \ge h_c
+\end{equation*}
+Si, de plus :
+\begin{equation*}
+ l + 1 = n
+\end{equation*}
+alors le contour passe entre $l$ et $l + 1$ :
+\begin{equation*}
+ h_{l + 1} < h_c \le h_l \lor h_l \le h_c < h_{l + 1}
+\end{equation*}
+On peut ainsi écrire :
+\begin{equation*}
+ l + 1 = n \Rightarrow h_{l + 1} < h_c \lor h_l = h_c < h_{l + 1}
+\end{equation*}
+Par ailleurs :
+\begin{align*}
+ & h_{l + 1} < h_c \Rightarrow l + 1 = n \\
+ & h_l < h_{l + 1} \Rightarrow l + 1 = n
+\end{align*}
+On a donc le critère :
+\begin{equation*}
+ l + 1 = n \Leftrightarrow h_{l + 1} < h_c \lor h_l < h_{l + 1}
+\end{equation*}
+De même, on montre que si l'extremum est un minimum alors :
+\begin{equation*}
+ l + 1 = n \Leftrightarrow h_{l + 1} > h_c \lor h_l > h_{l + 1}
+\end{equation*}
+Utilisation pour la procédure inside.
+
+Sans l'hypothèse (\ref{eq:non_degeneracy}), on peut avoir le cas
+pathologique où $h_c = h_l = h_{l + 1}$.
+\includegraphics{degeneracy}
+Le contour peut passer n'importe où entre $l$ et $l + 1$, cela dépend
+du fonctionnement de \verb+Contour_531+.
+
+Le plus éloigné des bons contours, si on le cherchait avec une
+précision infinie, contiendrait (dans la surface intérieure à ce
+contour) forcément les points voisins de l'extremum. Donc le chemin de
+sortie de ce bon contour le plus éloigné serait de longueur supérieure
+à 3.
+
+\section{Problèmes de dégénérescence}
+
+Possibilité que deux points adajacents aient exactement la même valeur
+de SSH. J. Tierny mentionne ce type de problème dans son
+\href{http://www-pequan.lip6.fr/~tierny/stuff/teaching/tierny_topologicalDataAnalysis.pdf}{cours}
+(§ 1.1.2 Range representation).
+
+J'ai pensé ajouter une valeur aléatoire :
+\begin{verbatim}
+call random_number(harvest)
+ssh(i, j) = ssh(i, j) + (1. + harvest) * 10. * epsilon(0.)
+\end{verbatim}
+Cela n'enlève pas forcément la dégénérescence. En simple précision par
+exemple, la valeur ajoutée est de l'ordre de \np{e-6}. Si par exemple
+pour deux positions successives de même SSH, les valeurs de harvest
+sont \np{0.53} et \np{0.50}, la différence entre ces deux valeurs de
+harvest est perdue par la troncature lors de l'ajout à ssh, si bien
+que la dégénérescence n'est pas levée.
+
+Autre idée : ne modifier la ssh que pour des extremums dégénérés. Par
+exemple, pour un maximum dégénéré, ajouter epsilon. Ou bien, pour un
+maximum dégénéré, le diminuer à la valeur autour la plus proche ?
+
+Il y aurait aussi l'idée (citée par J. Tierny) de ne pas toucher au
+champ de SSH mais de changer la définition de l'ordre. Mais il
+faudrait faire ces changements jusque dans les procédures de recherche
+de contour.
+
+Si j'admets la possibilité de dégénérescence alors les procédures
+impactées sont get\_1\_outerm, local\_extrema, inside et
+outermost\_possible\_level.
+
+Je cherche un bon contour en ne spécifiant qu'un seul point, qui doit
+être intérieur au bon contour, et un ensemble de points qui doivent
+être extérieurs. Si deux points voisins ont la même valeur de SSH et
+sont tous les deux extremums (non stricts donc), et si je les retiens
+tous les deux comme extremums valides, il ne peut exister de bon
+contour autour de l'un des deux. Si ce sont des extremums de signes
+différents, je peux me satisfaire de ne pas trouver de bon contour. Si
+ce sont des extremums de même signe, je voudrais trouver un bon
+contour autour des deux. Il faut donc n'en retenir qu'un des deux
+comme valide. Je peux prendre arbitrairement le premier rencontré dans
+l'ordre de parcours du tableau.
+
+Il faut prévoir dans get\_1\_outerm la possibilité de ne pas
+trouver de bon contour. La logique est impactée. \`A strictement
+parler, la valeur de innermost\_level ne permet plus de prévoir si
+l'amplitude va être suffisante, et on ne sait plus quelle valeur
+intiale de level\_good prendre dans get\_1\_outerm. Je garde
+néanmoins innermost\_level comme valeur intiale de level\_good. La
+conséquence est que je risque de ne pas détecter de contour autour
+d'un extremum dégénéré alors qu'on pourrait peut-être en trouver un en
+prenant une valeur initiale de level\_good plus proche de
+l'extremum. On peut aussi imaginer le cas bizarre de trois maximums
+(par exemple) égaux adjacents. Dans ce cas, l'algorithme choisi
+conserve deux de ces maximaux et il n'y a pas de bon contour possible,
+alors qu'on pourrait en principe en trouver. On se désintéresse de ces
+cas pathologiques.
+
+\section{Entrées et sorties}
+
+Format shapefile. Un type de géométrie par fichier donc il faut
+séparer les contours et les positions des extremums. Par ailleurs, on
+n'enregistre pas de contours pour les tourbillons interpolés. On
+pourrait créer un shapefile séparé pour les extremums interpolés
+(points), avec des méta-données supplémentaires sur les contours
+fictifs associés. En notant $m$ le numéro de processus MPI :
+\begin{itemize}
+\item \verb+extremum_interpolated_$m.shp+ : points
+\item \verb+extremum_interpolated_$m.dbf+ : valeur de SSH à l'extremum,
+ indice de date, indice de tourbillon à cette date, aire du contour
+ de SSH le plus éloigné, valeur de SSH sur le contour de SSH le plus
+ éloigné, aire du contour de vitesse maximale, valeur de vitesse sur
+ le contour de vitesse maximale
+\end{itemize}
+Mais il me semble plus simple et clair d'écrire les tourbillons
+interpolés avec les autres, en utilisant une forme NULL pour le
+polygone correspondant à un contour interpolé.
+\begin{itemize}
+\item \verb+extremum_$m.shp+ : points
+\item \verb+extremum_$m.dbf+
+ : valeur de SSH, indice de date, indice de tourbillon à cette date,
+ interpolé (logique), cyclone (logique), suff\_amp (logique), valeur
+ de vitesse sur le contour de vitesse maximale
+\item \verb+outermost_contour_$m.shp+ : polygones
+\item \verb+outermost_contour_$m.dbf+ : aire, valeur de SSH, indice de
+ date, indice de tourbillon à cette date, twice, radius4
+\item \verb+max_speed_contour_$m.shp+ : polygones
+\item \verb+max_speed_contour_$m.dbf+ : aire, valeur de SSH, indice de date, indice de tourbillon à cette date
+\end{itemize}
+Soit un total de 9 fichiers par processus en comptant les fichiers
+shx.
+
+L'indice de date et l'indice de tourbillon dans
+\verb+outermost_contour_$m.dbf+ et \verb+max_speed_contour_$m.dbf+
+sont en principe inutiles puisqu'ils sont aussi dans
+\verb+extremum_$m.dbf+ et que les tourbillons sont stockés dans le
+même ordre dans les trois shapefiles. La redondance me semble
+préférable pour la solidité et le confort d'utilisation. La
+méta-donnée logique \og interpolé \fg{} dans \verb+extremum_$m.dbf+
+est nécessaire parce qu'il est possible de ne pas trouver de contour
+extérieur autour d'un extremum détecté.
+
+J'avais d'abord mis la valeur de vitesse dans le fichier
+\verb+max_speed_contour_$m.dbf+ mais cette vitesse peut être associée
+au contour extérieur, avec un contour vide dans \verb+max_speed_contour_$m+.
+
+Les réels dans un fichier shp sont stockés en double précision donc
+les polygones occupent chacun : 16 B $\times$ nombre moyen de côtés
+d'un polygone. Soit, en comptant 30 côtés en moyenne par polygone,
+environ \np{0.5} KiB. (Ce nombre de côtés moyen m'a été communiqué par
+Rémi. On peut aussi l'obtenir grossièrement en supposant un rayon
+d'environ 100 km, soit un périmètre d'environ 600 km et une résolution
+de 25 km.) Comme chaque fichier shp peut occuper au maximum 2 GiB, il
+peut contenir au maximum \np{4e6} polygones. On ne peut donc mettre à
+coup sûr les polygones pour 16 ans dans un seul fichier, même en
+séparant cyclones et anticyclones. Faire un fichier par an ? Séparer
+en outre cyclones et anticyclones ?
+
+Une possibilité simple de stockage de graphe est la liste d'adjacences
+(en texte), c'est-à-dire, pour chaque sommet du graphe, une ligne :
+\begin{verbatim}
+indice date, indice tourbillon, degré sortant d
+\end{verbatim}
+suivie de d lignes :
+\begin{verbatim}
+indice date, indice tourbillon, poids
+\end{verbatim}
+Estimation de la taille du fichier, en comptant 10 caractères par
+entier, 13 caractères par réel, \np{5e7} sommets et deux arcs par
+sommet en moyenne : 5 GiB (au lieu de \np{0.5} GiB en binaire). Mais
+cela implique de connaître tous les arcs sortants d'un sommet au
+moment de l'écriture pour ce sommet. Dans notre algorithme, les arcs
+sortants sont découverts en plusieurs étapes (delta = 1 à 5) et nous
+ne voulons pas stocker le graphe en mémoire vive. Nous choisissons
+donc de stocker le graphe au format liste d'arcs (en texte),
+c'est-à-dire, pour chaque sommet du graphe, une ligne :
+\begin{verbatim}
+indice date, indice tourbillon, indice date, indice tourbillon, poids
+\end{verbatim}
+Fichier \verb+edgelist_$m.txt+. L'estimation de la taille du fichier, avec
+les mêmes hypothèses (en particulier deux arcs par sommet) est de 6
+GiB.
+
+Problème des sommets isolés. L'information est déjà présente, de façon
+implicite, dans les fichiers évoqués ci-dessus. En effet, on a le
+nombre de tourbillons à chaque date d'après la liste des identifiants
+de tourbillons dans les shapefiles, et on a la liste des tourbillons
+non isolés d'après la liste des arcs. Mais la liste des sommets isolés
+n'est donc pas très facilement trouvée. On ajoute par conséquent les
+fichiers de sortie suivants :
+\begin{description}
+\item[number\_eddies\_\$m.csv] Une colonne indice de date, une colonne
+ nombre de tourbillons visibles, une colonne nombre de tourbillons
+ interpolés.
+\item[isolated\_nodes\_\$m.txt] Sur chaque ligne : indice de date, indice
+ de tourbillon. Peut être lu comme une simple liste
+ d'adjacence. Remarque : par construction, les tourbillons interpolés
+ ne sont jamais isolés.
+\end{description}
+
+Les entrées-sorties sont dans : l'algorithme principal directement
+(lecture de longitude et latitude), get\_snapshot (lecture de SSH, u,
+v), successive\_overlap (écriture dans
+\verb+edgelist_$m.txt+, arcs entre k\_begin(m) et k\_end\_main\_loop(m)),
+non\_successive\_overlap directement (écriture dans
+\verb+edgelist_$m.txt+),
+write\_eddy (écriture dans extremum\_\$m, outermost\_contour\_\$m,
+max\_speed\_contour\_\$m), dispatch\_snapshot directement (écriture
+dans isolated\_nodes\_\$m, number\_eddies\_\$m). Les processus
+écrivent dans des fichiers indicés par le numéro de processus donc
+deux processus ne peuvent écrire dans le même fichier. Cf. figure
+(\ref{fig:input_output}).
+\begin{figure}[htbp]
+ \centering
+ \includegraphics[width=\textwidth]{input_output}
+ \caption{Entrées et sorties du programme en Fortran.}
+ \label{fig:input_output}
+\end{figure}
+
+\section{Déclarations pour l'algorithme principal}
+
+On suppose que le maillage est le même à toutes les dates et qu'il est
+régulier.
+
+Type dérivé polyline : défini dans \verb+Contour_531+. Mémoire
+utilisée par un scalaire de ce type :
+\begin{equation*}
+\mathrm{taille}(\mathtt{polyline}) = 2 (\mathtt{n\_points} + 1) \mathrm{mots}
+\end{equation*}
+soit 8 (n\_points + 1) B.
+
+Je définis un type dérivé ssh\_contour, associé à une ligne de niveau
+de SSH. Mémoire utilisée par un scalaire de ce type :
+\begin{align*}
+ \mathrm{taille}(\mathtt{ssh\_contour})
+ & = \mathrm{taille}(\mathtt{polyline}) + 2\ \mathrm{mots} \\
+ & = 2 (\mathtt{n\_points} + 2) \mathrm{mots}
+\end{align*}
+soit 8 (n\_points + 2) B. J'ai hésité à ajouter un champ speed,
+scalaire réel, pour une valeur de vitesse associée au contour. La
+vitesse utilisée peut être la composante azimutale par rapport à un
+centre défini indépendamment du contour. De ce point de vue, la
+définition de ce type dérivé serait un peu gênante. Par ailleurs, il
+ne serait pas très satisfaisant de laisser ce champ non défini pour le
+contour le plus extérieur, pour lequel nous n'avons pas besoin de
+calculer une vitesse.
+
+Les tailles ci-dessous sont calculées pour max\_delta = 4, 16 ans,
+soit environ 5000 dates, 8000 tourbillons par date, 30 côtés par
+polygone (sans différence entre contour extérieur et contour de
+vitesse maximale), 720 latitudes par 1440 longitudes.
+\begin{description}
+\item[max\_delta] Scalaire entier. Intervalle maximal d'indices de
+ dates auxquelles on cherche des recouvrements. Normalement 4.
+\item[longitude, latitude] Vecteurs de réels
+\item[ssh, u, v] Tableaux de réels, deux dimensions (nlon,
+ nlat). Environ 4 MiB chacun.
+\item[eddy] La comparaison entre out\_cont\%ssh et ssh\_extr donne la
+ distinction anticyclone -- cyclone mais il est quand même plus
+ confortable d'avoir cette information directement dans les
+ sorties. Pour cela, j'ajoute un champ cyclone. De même, le fait que
+ l'amplitude de ssh soit ou non en dessous du seuil est codé dans
+ extr\_map mais il me semble plus confortable de l'avoir en plus ici,
+ dans le champ suff\_amp. Mémoire utilisée par un scalaire de ce type
+ :
+ \begin{equation*}
+ \mathrm{taille}(\mathtt{eddy})
+ = 11\ \mathrm{mots} + \mathrm{taille}(\mathtt{out\_cont})
+ + \mathrm{taille}(\mathtt{speed\_cont})
+ \end{equation*}
+ Notons $n_o$ le nombre de points dans le contour extérieur
+ (out\_cont\%n\_points) et $n_s$ le nombre de points dans le contour
+ de vitesse maximale (speed\_cont\%n\_points). On obtient :
+ \begin{align*}
+ \mathrm{taille}(\mathtt{eddy})
+ & = [2 (n_o + n_s) + 19]\ \mathrm{mots} \\
+ & = 4 [2 (n_o + n_s) + 19]\ \mathrm{B}
+ \end{align*}
+soit environ \np{0.5} KiB.
+\item[snapshot\%list\_vis] On a besoin d'encapsuler le vecteur de type
+ eddy dans un scalaire de type snapshot pour pouvoir considérer un
+ vecteur de snapshots. Notons $n_{ve}$ le nombre de tourbillons à une
+ date donnée (number\_vis\_eddies). Mémoire utilisée par ce champ :
+ \begin{equation*}
+ n_{ve} \langle \mathrm{taille}(\mathtt{eddy}) \rangle
+ = n_{ve}
+ [2 (\langle n_o \rangle + \langle n_s \rangle) + 19]\ \mathrm{mots}
+ \end{equation*}
+ soit environ 4 MiB.
+\item[snapshot\%extr\_map] Notons $n_\lambda$ le nombre de longitudes,
+ $n_\phi$ le nombre de latitudes. Taille : $n_\lambda n_\phi$ mots,
+ soit environ 4 MiB.
+\item[snapshot\%ind\_extr] Taille : $2 n_{ve}$ mots, soit environ
+ \np{0.06} MiB.
+\item[flow] : Vecteur de type snapshot, de taille max\_delta + 1,
+ chaque élément correspondant à une date. Environ 40 MiB.
+\item[dist\_lim] Scalaire entier constant, = 3° / résolution.
+\item[weight\_delta] Scalaire réel.
+\end{description}
+On a :
+\begin{align*}
+ \mathrm{taille}(\mathtt{snapshot})
+ & = n_{ve} \langle \mathrm{taille}(\mathtt{eddy}) \rangle
+ + (n_\lambda n_\phi + 2 n_{ve} + 2)\ \mathrm{mots} \\
+ & = [n_{ve} (2 (\langle n_o \rangle + \langle n_s \rangle) + 21)
+ + n_\lambda n_\phi + 2]\ \mathrm{mots}
+\end{align*}
+D'où la taille d'un scalaire de type snapshot : environ 8 MiB.
+
+J'avais d'abord pensé mettre dans le type eddy les champs ind\_lon et
+ind\_lat. Ce n'est pas très cohérent parce que les indices se réfèrent
+au maillage global, alors que les autres champs sont des propriétés
+intrinsèques. Ainsi, j'étais embêté pour définir un tourbillon à
+partir d'une procédure travaillant sur une section du maillage
+global. J'ai donc choisi donc de remplacer ces champs par un champ
+ind\_extr dans snapshot. Et j'ajoute à la place des champs longitude,
+latitude dans le type eddy.
+
+i désigne un indice de tourbillon, j un indice de position dans la
+fenêtre temporelle (entre 1 et max\_delta + 1) et k un indice de
+date. Cf. figure (\ref{fig:window}) et algorithme
+(\ref{alg:main_sequential}).
+\begin{figure}[htbp]
+ \centering
+ \includegraphics{window}
+ \caption{Indice de position dans la fenêtre glissante et indice de
+ date.}
+ \label{fig:window}
+\end{figure}
+\begin{algorithm}[htbp]
+ \begin{algorithmic}
+ \STATE \COMMENT{\{max\_delta $\ge 1$; n\_dates $\ge$ max\_delta +
+ 1\}}
+
+ \STATE \COMMENT{prologue}
+
+ \STATE entrer(longitude, latitude)
+ \STATE \COMMENT{longitude et latitude dans l'ordre croissant}
+ \FOR{k = 1 \TO max\_delta + 1}
+
+ \STATE flow(k) = get\_snapshot(k)
+
+ \ENDFOR
+ \FOR{k = 2 \TO max\_delta + 1}
+
+ \STATE appel de successive\_overlap(k, k, flow)
+ \ENDFOR
+ \FOR{delta = 2 \TO max\_delta}
+
+ \FOR{k = delta + 1 \TO max\_delta + 1}
+
+ \STATE appel de non\_successive\_overlap(k, k, delta, flow)
+
+ \ENDFOR
+ \ENDFOR
+
+ \STATE \COMMENT{boucle principale}
+
+ \FOR{k = max\_delta + 2 \TO n\_dates}
+
+ \STATE appel de dispatch\_snapshot(flow(1), k - max\_delta - 1)
+ \STATE flow(:max\_delta) = flow(2:)
+
+ \STATE flow(max\_delta + 1) = get\_snapshot(k)
+
+ \STATE appel de successive\_overlap(max\_delta + 1, k, flow)
+
+ \FOR{delta = 2 \TO max\_delta}
+
+ \STATE appel de non\_successive\_overlap(max\_delta + 1, k, delta,
+ flow)
+
+ \ENDFOR
+ \ENDFOR
+ \STATE \COMMENT{épilogue}
+
+ \FOR{j = 1 \TO max\_delta + 1}
+
+ \STATE appel de dispatch\_snapshot(flow(j), n\_dates - max\_delta - 1
+ + j)
+ \ENDFOR
+ \end{algorithmic}
+ \caption{Algorithme principal, séquentiel}
+ \label{alg:main_sequential}
+\end{algorithm}
+
+\section{Décomposition de domaine}
+
+Cf. figure (\ref{fig:processes}).
+\begin{figure}[htbp]
+ \centering
+ \includegraphics[width=\textwidth]{processes}
+ \caption{Raccordement des processus. Les comparaisons représentées
+ ici en couleurs sont faites par le processus m. En rouge dans la
+ boucle principale, en bleu, vert et marron dans l'épilogue.}
+ \label{fig:processes}
+\end{figure}
+Il faut, pour m < n\_proc :
+\begin{equation}
+ \label{eq:stitching}
+ \mathtt{k\_end}(m)
+ = \mathtt{k\_begin}(m + 1) + \mathtt{max\_delta} - 1
+\end{equation}
+
+k\_end - k\_begin + 1 appels à get\_snapshot, pour les dates k\_begin
+à k\_end. k\_end - k\_begin + 1 appels à dispatch\_snapshot, pour les
+dates k\_begin à k\_end. Si m > 1 alors max\_delta envois. Si m <
+n\_proc alors max\_delta réceptions.
+
+Nombre d'appels à get\_snapshot avec lecture. Si m < n\_proc : k\_end -
+k\_begin - max\_delta + 1. Si m = n\_proc : k\_end - k\_begin + 1.
+
+successive\_overlap est appelé pour tous les indices k compris entre
+k\_begin + 1 et k\_end\_main\_loop, seulement. Nombre d'appels à
+successive\_overlap : k\_end\_main\_loop - k\_begin. Donc si m <
+n\_proc : k\_end - k\_begin - max\_delta + 1 appels. Si m = n\_proc : k\_end
+- k\_begin appels.
+
+Nombre d'appels à non\_successive\_overlap dans le prologue :
+\begin{equation*}
+ \sum_{\delta = 2} ^{\mathtt{max\_delta}}
+ (\mathtt{max\_delta} - \delta +1)
+ = \frac{\mathtt{max\_delta} (\mathtt{max\_delta} - 1)}{2}
+\end{equation*}
+Dans la boucle principale : (max\_delta - 1)(k\_end\_main\_loop -
+k\_begin - max\_delta). Dans l'épilogue, si m < n\_proc :
+$\frac{\mathtt{max\_delta} (\mathtt{max\_delta} - 1)}{2}$. Total si m
+< n\_proc : (max\_delta - 1)(k\_end - k\_begin - max\_delta +
+1). Total si m = n\_proc : (max\_delta - 1)(k\_end - k\_begin -
+$\frac{\mathtt{max\_delta}}{2}$).
+
+Si k\_end - k\_begin était le même pour m = n\_proc et m <
+n\_proc, le processus m=n\_proc aurait des plus grands nombres
+d'appels aux trois procédures. Notons $\Delta k$ = k\_end - k\_begin,
+considéré comme une fonction de m. Si, pour m < n\_proc :
+\begin{equation}
+ \label{eq:charge}
+ \Delta k(m) \approx \Delta k(\mathtt{n\_proc}) + \mathtt{max\_delta}
+\end{equation}
+alors la charge devrait être à peu près équilibrée.
+
+On peut prendre, pour tout m :
+\begin{equation*}
+ \mathtt{k\_begin}(m)
+ = 1 + E[(m - 1) \mathtt{n\_dates} / \mathtt{n\_proc}]
+\end{equation*}
+(On prend la partie entière au lieu de nint pour pouvoir démontrer des
+encadrements.) D'où, avec la contrainte (\ref{eq:stitching}), pour m
+< n\_proc :
+\begin{equation*}
+ \mathtt{k\_end}(m)
+ = E(m\ \mathtt{n\_dates} / \mathtt{n\_proc}) + \mathtt{max\_delta}
+\end{equation*}
+Alors on peut montrer que $\Delta k(m) - \Delta k(\mathtt{n\_proc})$
+vaut max\_delta ou max\_delta - 1. On a donc l'approximation
+(\ref{eq:charge}). Chaque processus fait aussi le même nombre de
+lectures de ssh, u, v.
+
+Pour n\_proc $\ne 1$ et m < n\_proc, le processus m + 1 ne doit pas
+relire une date déjà lue par le processus m. On doit donc avoir :
+\begin{equation}
+ \label{eq:min_delta_k}
+ \mathtt{k\_begin}(m + 1) > \mathtt{k\_begin}(m) + \mathtt{max\_delta}
+\end{equation}
+Ce qui est équivalent à :
+\begin{equation*}
+ E(m \mathtt{n\_dates} / \mathtt{n\_proc})
+ - E[(m - 1) \mathtt{n\_dates} / \mathtt{n\_proc}]
+ > \mathtt{max\_delta}
+\end{equation*}
+Or on peut montrer que :
+\begin{equation*}
+ E(m \mathtt{n\_dates} / \mathtt{n\_proc})
+ - E[(m - 1) \mathtt{n\_dates} / \mathtt{n\_proc}]
+ > E(\mathtt{n\_dates} / \mathtt{n\_proc}) - 1
+\end{equation*}
+Il suffit donc que :
+\begin{equation*}
+ E(\mathtt{n\_dates} / \mathtt{n\_proc}) \ge \mathtt{max\_delta} + 1
+\end{equation*}
+
+L'équation (\ref{eq:stitching}) et l'inégalité (\ref{eq:min_delta_k})
+impliquent que, pour m < n\_proc :
+\begin{equation}
+ \label{eq:delta_k}
+ \mathtt{k\_begin} + \mathtt{max\_delta}
+ \le \mathtt{k\_end} - \mathtt{max\_delta}
+\end{equation}
+donc les appels de get\_snapshot dans le prologue font des
+lectures. Les appels de get\_snapshot dans l'épilogue font des
+réceptions. L'inégalité (\ref{eq:delta_k}) implique en outre que, pour
+m < n\_proc, la boucle principale fait au moins une itération. Si m <
+n\_proc : tous les appels à get\_snapshot de la boucle principale sauf
+le dernier font des lectures ; le dernier appel à get\_snapshot de la
+boucle principale fait une réception.
+
+Si $m \ge 2$ alors :
+\begin{equation*}
+ \mathtt{k\_end\_main\_loop}(m - 1) = \mathtt{k\_begin}(m)
+\end{equation*}
+
+Exemples : cf. tableaux (\ref{tab:m2}) et (\ref{tab:m3}) et figure
+(\ref{fig:m3}).
+\begin{table}[htbp]
+ \centering
+ \begin{tabular}{llll}
+ m & k\_begin & k\_end & k\_end\_main\_loop \\
+ \hline
+ 1 & 1 & 9 & 6 \\
+ 2 & 6 & 10 & 10
+ \end{tabular}
+ \caption{max\_delta = 4, n\_dates = 10, n\_proc = 2. m = 1 : lecture
+ dates 1 à 5, écriture date 1, réception date 6, écriture dates 2 à
+ 4, réception dates 7 à 9, écriture dates 5 à 9. m = 2 : lecture 6
+ à 10, envoi 6 à 9, écriture 10.}
+ \label{tab:m2}
+\end{table}
+\begin{table}[htbp]
+ \centering
+ \begin{tabular}{llll}
+ m & k\_begin & k\_end & k\_end\_main\_loop \\
+ \hline
+ 1 & 1 & 9 & 6 \\
+ 2 & 6 & 14 & 11 \\
+ 3 & 11 & 15 & 15
+ \end{tabular}
+ \caption{max\_delta = 4, n\_dates = 15, n\_proc = 3}
+ \label{tab:m3}
+\end{table}
+\begin{figure}[htbp]
+ \centering
+ \includegraphics[width=\textwidth]{m3}
+ \caption{max\_delta = 4, n\_dates = 15, n\_proc = 3. En bleu le
+ prologue, en noir la boucle principale, en rouge l'épilogue. Pour
+ un processus donné, les actions sur une même colonne sont dans une
+ même boucle de l'algorithme principal.}
+ \label{fig:m3}
+\end{figure}
+
+Le dispatch\_snapshot dans la dernière boucle sur j n'est pas
+forcément une écriture. Exemple : n\_proc = 2, max\_delta = 4,
+n\_dates = 10, m = 2.
+
+Si m = 1 alors write\_eddy est appelé pour les dates k\_begin à
+k\_end. Si m > 1 alors write\_eddy est appelé pour les dates
+k\_begin + max\_delta à k\_end. Avec l'équation (\ref{eq:stitching}),
+une date donnée est bien écrite par un seul processus.
+
+\section{Algorithme principal, parallèle}
+
+\begin{algorithmic}[1]
+ \STATE \COMMENT{\{max\_delta $\ge 1$; n\_dates $\ge$ max\_delta + 1;
+ n\_proc
+ $\le E\left(\frac{\mathtt{n\_dates}}{\mathtt{max\_delta} + 1}
+ \right)$\}}
+
+ \STATE \COMMENT{prologue}
+
+ \STATE m = numéro de processus
+
+ \IF{m = 1}
+
+ \STATE entrer(longitude, latitude)
+ \STATE \COMMENT{longitude et latitude dans l'ordre croissant}
+
+ \STATE envoyer longitude, latitude à tous les processus
+
+ \ELSE
+
+ \STATE recevoir longitude, latitude du processus 1
+ \ENDIF
+
+ \STATE k\_begin = 1 + E[(m - 1) n\_dates / n\_proc]
+
+ \IF{m < n\_proc}
+
+ \STATE k\_end = E(m n\_dates / n\_proc) + max\_delta
+
+ \STATE k\_end\_main\_loop = k\_end - max\_delta + 1
+ \ELSE
+
+ \STATE k\_end = n\_dates
+ \STATE k\_end\_main\_loop = k\_end
+ \ENDIF
+
+ \FOR{k = k\_begin \TO k\_begin + max\_delta}
+
+ \STATE flow(k - k\_begin + 1) = get\_snapshot(k)
+ \COMMENT{lecture}
+
+ \ENDFOR
+ \FOR{k = k\_begin + 1 \TO k\_begin + max\_delta}
+
+ \STATE appel de successive\_overlap(k - k\_begin + 1, k, flow)
+ \ENDFOR
+ \FOR{delta = 2 \TO max\_delta}
+
+ \FOR{k = k\_begin + delta \TO k\_begin + max\_delta}
+
+ \STATE appel de non\_successive\_overlap(k - k\_begin + 1, k, delta,
+ flow)
+
+ \ENDFOR
+ \ENDFOR
+
+ \STATE \COMMENT{boucle principale}
+
+ \FOR{k = k\_begin + max\_delta + 1 \TO k\_end\_main\_loop}
+
+ \STATE appel de dispatch\_snapshot(flow(1), k - max\_delta - 1)
+
+ \STATE flow(:max\_delta) = flow(2:)
+
+ \STATE flow(max\_delta + 1) = get\_snapshot(k)
+
+ \STATE appel de successive\_overlap(max\_delta + 1, k, flow)
+
+ \FOR{delta = 2 \TO max\_delta}
+
+ \STATE appel de non\_successive\_overlap(max\_delta + 1, k, delta,
+ flow)
+
+ \ENDFOR
+ \ENDFOR
+ \STATE \COMMENT{épilogue}
+
+ \FOR{k = k\_end\_main\_loop + 1 \TO k\_end}
+
+ \STATE \COMMENT{m < n\_proc et k $\ge$ k\_end - max\_delta + 2}
+ \STATE appel de dispatch\_snapshot(flow(1), k - max\_delta - 1)
+
+ \STATE flow(:max\_delta) = flow(2:)
+
+ \STATE flow(max\_delta + 1) = get\_snapshot(k)
+ \COMMENT{réception} \STATE \COMMENT{stitching}
+
+ \FOR{delta = k - k\_end + max\_delta \TO max\_delta}
+
+ \STATE appel de non\_successive\_overlap(max\_delta + 1, k, delta,
+ flow)
+
+ \ENDFOR
+ \ENDFOR
+
+ \FOR{j = 1 \TO max\_delta + 1}
+
+ \STATE appel de dispatch\_snapshot(flow(j), k\_end - max\_delta - 1
+ + j)
+
+ \ENDFOR
+\end{algorithmic}
+
+\section{Résumé du problème de géométrie algorithmique}
+
+\begin{algorithmic}
+ \STATE recherche de tous les extremums locaux
+
+ \FOR{chaque extremum local}
+
+ \STATE recherche du \og bon contour le plus extérieur\fg{}
+ \ENDFOR
+\end{algorithmic}
+
+\og Bon contour\fg{} : contient l'extremum considéré et aucun
+autre. \og Le plus extérieur\fg{} : plus grande différence de niveau
+avec l'extremum considéré, ou encore amplitude maximale.
+
+Recherche du bon contour le plus extérieur : par dichotomie.
+\begin{algorithmic}
+ \STATE $a_1 = 0$
+
+ \STATE $a_2$ = recherche d'une borne supérieure pour
+ l'amplitude en parcourant les points sur un rayon à partir de
+ l'extremum jusqu'à ce que le sens de variation de SSH change
+
+ \WHILE{$a_2 - a_1 >$ précision}
+
+ \STATE recherche du bon contour à l'amplitude
+ $\frac{a_1 + a_2}{2}$
+
+ \IF{bon contour trouvé}
+
+ \STATE $a_1 = \frac{a_1 + a_2}{2}$
+
+ \ELSE
+
+ \STATE $a_2 = \frac{a_1 + a_2}{2}$
+
+ \ENDIF
+
+ \ENDWHILE
+\end{algorithmic}
+
+Recherche d'un bon contour à une amplitude donnée. On cherche tous les
+contours à cette amplitude. Pour chaque contour, pour chaque extremum,
+on regarde si l'extremum est à l'intérieur du contour (test de point
+intérieur à polygone).
+
+\section{Sous-algorithmes}
+
+\subsection{get\_snapshot}
+
+Pour chaque tourbillon du résultat, le contour le plus extérieur doit
+être fermé, doit contenir l'extremum de SSH du tourbillon considéré et
+ne doit contenir aucun autre extremum de SSH avec une amplitude
+significative. Amplitude significative : supérieure à un certain
+minimum, entre l'extremum de SSH et le contour le plus
+lointain.
+
+Dans une première version qui n'envisageait pas d'amplitude minimum,
+j'avais une procédure get\_eddy qui définissait un tourbillon, et qui
+trouvait donc le contour le plus extérieur et le contour de maximum de
+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.
+
+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
+peut être artificiellement diminuée. (Ce n'est pas le même cas qu'un
+tourbillon près d'une côte.) C'est-à-dire qu'on pourrait trouver un
+contour plus extérieur avec un domaine élargi. Il pourrait être
+intéressant d'éliminer un tel tourbillon. Comment le détecter ? Le
+fait d'aller jusqu'au bord du domaine dans outermost\_possible\_level
+n'est pas un critère suffisant, parce que outermost\_possible\_level
+ne teste que la direction est. Il faudrait tester après l'appel de
+get\_1\_outerm si le cadre le plus extérieur touche le bord du
+domaine.
+
+\subsection{max\_speed\_contour\_ssh}
+
+Il manque des tests pour lesquels direction vaut 1 et 4.
+
+\subsection{set\_max\_speed}
+
+Notons $d$ la distance dans l'espace des indices entre l'extremum
+cible et le countour extérieur. Par définition de la partie entière,
+les points à distance $E(d)$ de l'extremum cible sont à l'intérieur du
+contour. Certains points à distance $E(d) + 1$ peuvent être à
+l'extérieur. D'où : radius = $E(d) + 1$.
+
+Dans le cas où radius vaut 1, je calcule max\_speed. Ce calcul ne sera
+utile que dans weight si le tourbillon a une intersection avec avec un
+tourbillon à une autre date. Donc, pour certains tourbillons, ce
+calcul est inutile. Difficile d'éviter ces calculs inutiles. Il
+faudrait les reporter à weight mais alors il faudrait avoir accès dans
+weight à u et v, c'est-à-dire qu'il faudrait conserver u, v à
+plusieurs dates en dehors de get\_snpashot. Cela remettrait en
+question tout l'algorithme.
+
+\subsection{successive\_overlap}
+
+j, k, scalaires entiers, données. flow : donnée-résultat. Cherche les
+arcs entre flow(j - 1) et flow(j), correspondant aux dates k - 1 et
+k. Écrit ces arcs dans \verb+graph.txt+. Met à jour flow(j -
+1)\%list\_vis\%delta\_out et
+flow(j)\%list\_vis\%delta\_in. Accélération possible en ne comparant à
+i1 que les trois tourbillons les plus proches (distance entre les
+extremums), comme dans la fonction \verb+Eddy_Intersection+ du
+programme pour Matlab.
+
+On ne peut éviter la copie de polylines due à l'alternative contour
+extérieur -- contour de vitesse maximale. Pour ne pas faire la copie
+dans polygon\_1 et polygon\_2, il faudrait faire des copies de
+out\_cont dans speed\_cont dans get\_snapshot, on ne gagnerait rien en
+quantité de copie, et on perdrait en mémoire vive occupée.
+
+\subsection{non\_successive\_overlap}
+
+j, k, delta $\in \{2, \dots, \mathtt{max\_delta}\}$ : scalaires
+entiers, données. flow : donnée-résultat. Cherche les arcs entre
+flow(j - delta) et flow(j), correspondant aux dates k - delta et
+k. Écrit ces arcs dans \verb+graph.txt+, écrit les tourbillons
+interpolés dans les shapefiles via interpolate\_eddy et
+write\_eddy. Met à jour flow(j - delta)\%list\_vis\%delta\_out,
+flow(j)\%list\_vis\%delta\_in, flow\%number\_eddies. Cf. algorithme
+(\ref{alg:non_successive_overlapping}).
+\begin{algorithm}[htbp]
+ \begin{algorithmic}
+ \STATE \COMMENT{recouvrements à dates non successives}
+
+ \FOR{i1 = 1 \TO flow(j - delta)\%number\_vis\_eddies}
+
+ \IF{flow(j - delta)\%list\_vis(i1)\%suff\_amp}
+
+ \STATE polygon\_1 = merge(flow(j -
+ delta)\%list\_vis(i1)\%speed\_cont\%polyline, flow(j -
+ delta)\%list\_vis(i1)\%out\_cont\%polyline, flow(j -
+ delta)\%list\_vis(i1)\%max\_speed $\ne$ 0)
+
+ \STATE i1\_lon = flow(j - delta)\%ind\_extr(1, i1)
+
+ \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))}
+
+ \IF{flow(j)\%list\_vis(i2)\%suff\_amp et (flow(j -
+ delta)\%list\_vis(i1)\%delta\_out $\ge$ delta ou
+ flow(j)\%list\_vis(i2)\%delta\_in $\ge$ delta) et flow(j -
+ delta)\%list\_vis(i1)\%cyclone $\Leftrightarrow$
+ flow(j)\%list\_vis(i2)\%cyclone}
+
+ \STATE polygon\_2 =
+ merge(flow(j)\%list\_vis(i2)\%speed\_cont\%polyline,
+ flow(j)\%list\_vis(i2)\%out\_cont\%polyline,
+ flow(j)\%list\_vis(i2)\%max\_speed $\ne$ 0)
+
+ \IF{polygon\_1 et polygon\_2 se recouvrent}
+
+ \STATE weight\_delta = weight(flow(j - delta)\%list\_vis(i1),
+ flow(j)\%list\_vis(i2))
+
+ \STATE flow(j - delta + 1)\%number\_eddies += 1
+
+ \STATE appel de interpolate\_eddy(flow(j -
+ delta)\%list\_vis(i1), flow(j)\%list\_vis(i2), j - delta, j, j - delta + 1,
+ k - delta + 1, flow(j - delta + 1)\%number\_eddies)
+
+ \STATE écrire dans graph.txt : k - delta, i1, k - delta + 1,
+ flow(j - delta + 1)\%number\_eddies, weight\_delta
+
+ \FOR{j\_interp = j - delta + 2 \TO j - 1}
+
+ \STATE flow(j\_interp)\%number\_eddies += 1
+
+ \STATE appel de interpolate\_eddy(flow(j -
+ delta)\%list\_vis(i1), flow(j)\%list\_vis(i2), j - delta, j, j\_interp, k -
+ j + j\_interp, flow(j\_interp)\%number\_eddies)
+
+ \STATE écrire dans graph.txt : k - j + j\_interp - 1,
+ flow(j\_interp - 1)\%number\_eddies, k - j + j\_interp,
+ flow(j\_interp)\%number\_eddies, weight\_delta
+ \ENDFOR
+
+ \STATE écrire dans graph.txt : k - 1, flow(j - 1)\%number\_eddies, k,
+ i2, weight\_delta
+
+ \STATE flow(j - delta)\%list\_vis(i1)\%delta\_out = delta
+
+ \STATE flow(j)\%list\_vis(i2)\%delta\_in = delta
+ \ENDIF
+ \ENDIF
+ \ENDFOR
+ \ENDIF
+ \ENDFOR
+ \end{algorithmic}
+ \caption{Sous-algorithme non\_successive\_overlap(j, k, delta,
+ flow)}
+ \label{alg:non_successive_overlapping}
+\end{algorithm}
+
+\subsection{local\_extrema}
+
+a(ind\_extr(1, k), ind\_extr(2, k)) est un extremum de a par rapport
+aux huit points voisins. extr\_map vaut 0 en un point où il n'y a pas
+d'extremum, le numéro de l'extremum dans ind\_extr en un point
+extremum. Accélération possible, au prix de tableaux supplémentaires,
+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.
+
+Cf. aussi fontion \verb+skimage.feature.peak_local_max+ dans le paquet
+Python scikit-image. La fonction applique
+\verb+scipy.ndimage.maximum_filter+ et compare le résultat au tableau
+de départ. \verb+scipy.ndimage.maximum_filter+ applique
+\verb+scipy.ndimage.maximum_filter1d+ dans chaque dimension. Cf. aussi
+Douglas (1996 k0968).
+
+\subsection{get\_1\_outerm}
+
+Au lieu d'initialiser l'amplitude maximale à 4 m plus du bruit, comme
+dans le programme pour Matlab :
+\begin{algorithmic}
+ \STATE appel de random\_number(level\_try)
+
+ \STATE level\_try = e\%ssh\_extr $\pm$ [4 + (2 * level\_try -
+ 1) * precision] \COMMENT{un peu de bruit sur l'amplitude, autour de
+ 4}
+\end{algorithmic}
+et de supposer qu'il n'y a pas de bon contour avec cette amplitude, on
+cherche une amplitude maximale possible. Dans le cas le plus simple où
+noise\_around est faux, on purrait parcourir les points dans une
+direction jusqu'à ce que le sens de variation de SSH change. Dans mes
+tests, le nombre d'itérations dans la recherche par dichotomie est
+typiquement une dizaine, pour une précision de \np{0.5} mm. En
+considérant simplement l'extremum de ssh sur la fenêtre, la taille de
+l'intervalle initial [innermost\_level, level\_try] est plus
+grande. Sur mes tests, elle est plus grande typiquement d'un facteur 2
+environ seulement, donc le nombre d'itérations change peu.
+
+Pour l'amplitude minimale, si je prends simplement l'amplitude
+correspondant à innermost\_level alors l'algorithme échoue lorsque
+l'extremum est à un point du bord et que le point à innermost\_level
+est au bord. Dans ce cas, good\_contour échoue parce que Contour\_531
+ne ferme pas les contours passant par le bord. J'avais essayé, pour
+éviter le problème, de prendre une amplitude minimale légèrement
+inférieure mais cela pose des problèmes de cohérence dans
+l'algorithme.
+
+\subsection{interpolate\_eddy}
+
+Tous les arguments sont des données. e1, e2 : scalaires de type
+eddy. j1, j2, indices de position dans la fenêtre correspondant à e1,
+e2. j : indice de position dans la fenêtre auquel on veut
+interpoler. k, i : indices de date et de tourbillon. Cf. algorithme
+(\ref{alg:interpolate_eddy}).
+\begin{algorithm}[htbp]
+ \begin{algorithmic}
+ \STATE e = interpolation entre e1 et e2
+
+ \STATE appel de write\_eddy(e, k, i)
+ \end{algorithmic}
+ \caption{Sous-algorithme interpolate\_eddy(e1, e2, j1, j2, j, k, i)}
+ \label{alg:interpolate_eddy}
+\end{algorithm}
+On peut mettre le champ suff\_amp à faux pour un tourbillon interpolé,
+en cohérence avec un champ out\_cont vide. Il faut alors un
+champ interpolated pour distinguer les tourbillons interpolés de ceux
+pour lesquels un extremum a été détecté mais sans obtenir de contour
+extérieur.
+
+\subsection{good\_contour}
+
+Trouve, si elle existe, une ligne de niveau à un niveau donné,
+contenant un point donné et ne contenant pas un ensemble de points
+donnés. Le résultat est un polygone à zéro point s'il n'existe pas de
+ligne de niveau satisfaisante. Correspond à une partie de la fonction
+extraction\_eddies dans le programme pour Matlab.
+
+\subsection{write\_eddy}
+
+On pourrait passer en arguments d'entrée les numéros des
+champs dans les fichiers DBF mais cela ferait 13
+arguments.
+
+\subsection{weight}
+
+D'autant plus proche de 0 que les tourbillons sont
+ressemblants. Comment calculer ce poids ? Comment calculer le poids
+faisant intervenir un ou deux tourbillons interpolés ? Pour faire au
+plus simple, j'ai supposé que l'on utilisait pour les arcs faisant
+intervenir les tourbillons interpolés le poids associé aux tourbillons
+visibles entre lesquels on interpole. Normalement, ce choix doit
+donner un poids plus faible à ces arcs.
+
+\section{Tests}
+
+Pour faire des tests : données au 1\ier{} janvier 2006. Tests sur
+différents domaines. Cf. figure (\ref{fig:regions}).
+\begin{figure}[htbp]
+ \centering
+ \includegraphics[width=\textwidth]{regions}
+ \caption{Régions pour les tests.}
+ \label{fig:regions}
+\end{figure}
+Régions 1, 2 et 3 : 29 $\times$ 17, 57 $\times$ 33; 120 $\times$ 120,
+contenant respectivement 8, 27 et 216 extremums.
+
+Pour la région 1, coin en :
+\begin{equation*}
+ \lambda = \np{5.125}^\circ, \phi = -\np{36.375}^\circ
+\end{equation*}
+soit $(\np[rad]{0.0894}, - \np[rad]{0.6349})$. Indices base 1 du coin
+: 21, 215. Pas de grille : \np{0.25}°, soit \np[rad]{0.0044}. Nombre
+de points de grille : $29 \times 17$. 5 minimums et 3 maximums dans ce
+cadre. Cf. figure (\ref{fig:extrema}).
+\begin{figure}[htbp]
+ \centering
+ \includegraphics[width=\textwidth]{extrema}
+ \caption{Région 1, lignes de niveau de SSH (sans signification
+ particulière). Maximums : points rouges, minimums : points bleus.}
+ \label{fig:extrema}
+\end{figure}
+Minimum cible en :
+\begin{equation*}
+ \lambda = \np{9.625}^\circ, \phi = - \np{33.875}^\circ
+\end{equation*}
+Indices base 1 de ce minimum : 39, 225. SSH en ce mimimum :
+\np{0.2759}. Extremum numéro 6 dans ce cadre. Contour le plus
+intérieur : SSH = \np{0.293300003} m. Contour le plus extérieur : SSH
+= \np{0.460744} m. Contour de vitesse maximale : SSH = \np{0.3815}
+m. Cf. figures (\ref{fig:test_get_eddy}) et (\ref{fig:test_get_snapshot}).
+\begin{figure}[htbp]
+ \centering
+ \includegraphics[width=\textwidth]{test_get_eddy}
+ \caption{Projection stéréographique centrée sur un minimum de
+ SSH. Champ de vitesse, contour de SSH le plus extérieur en bleu,
+ contour de SSH de vitesse maximale en vert. Le minimum de SSH est
+ marqué par un point.}
+ \label{fig:test_get_eddy}
+\end{figure}
+\begin{figure}[htbp]
+ \centering
+ \includegraphics[width=\textwidth]{test_get_snapshot}
+ \caption{Région 1, sans minimum d'amplitude. Contours extérieurs et
+ contours de vitesse maximale. En rouge les anticyclones, en bleu
+ les cyclones. Croix : contour extérieur sans contour de maximum de
+ vitesse distinct. Cercles : contour extérieur et contour de
+ maximum de vitesse distinct. Les contours extérieurs 1, 2 et 8
+ passent près de l'extremum et il n'y a donc pas de contour de
+ vitesse maximale correspondant. Les contours extérieurs 1 et 7
+ touchent le bord et sont donc peut-être artificiellement
+ sous-estimés.}
+ \label{fig:test_get_snapshot}
+\end{figure}
+
+Comparaison sur la région 1 avec le programme Matlab, pour les
+tourbillons 4 et 5. Les coordonnées des extremums sont bien les
+mêmes. Les SSH sur les contours extérieurs sont pratiquement les mêmes
+: différence de l'ordre de \np{0.1} mm. Du coup les contours
+extérieurs sont quasi-identiques. Les rayons des contours extérieurs
+sont presque les mêmes aussi : j'ai \np{87.7} km au lieu de \np{87.5}
+km et \np{90.80} km au lieu de \np{90.81} km. Les SSH sur les contours
+de vitesse max diffèrent de 2 cm environ. Cf. figure
+(\ref{fig:comparaison_Matlab}).
+\begin{figure}[htbp]
+ \centering
+ \includegraphics[width=\textwidth]{comparaison_201}
+ \includegraphics[width=\textwidth]{comparaison_202}
+ \caption{Comparaison avec le programme Matlab. Pour chacun des deux
+ tourbillons, le contour sans marqueur est à la valeur de SSH
+ donnée par le programme Matlab.}
+ \label{fig:comparaison_Matlab}
+\end{figure}
+J'obtiens comme vitesse maximale : $\np[m s^{-1}]{0.29}$ au lieu de
+$\np[m s^{-1}]{0.30}$ et $\np[m s^{-1}]{0.23}$ au lieu de
+$\np[m s^{-1}]{0.26}$.
+
+Test sur la région 3, sans minimum d'amplitude. 14 extremums sans
+contour extérieur (suff\_amp faux), dont 5 sur les bords et 11
+dégénérés. Avec minimum d'amplitude : 34 extremums pour lesquels
+suff\_amp est faux. Je vérifie que, pour chaque extremum, si suff\_amp
+est faux dans le cas sans minimum, suff\_amp reste faux dans le cas
+avec minimum. Il y a donc 20 extremums qui passent à suff\_amp faux à
+cause du minimum d'amplitude. Je vérifie que, pour ces 20 extremums,
+le contour extérieur dans le cas sans minimum est le même que dans le
+cas avec minimum. Entre le cas sans minimum et le cas avec minimum,
+seulement trois contours extérieurs grossissent (je le vois en
+comparant les fichiers outermost\_contour\_1.dbf) : dans le cas avec
+minimum, ils englobent un extremum d'amplitude insuffisante. Dans le
+cas avec minimum, 48 contours extérieurs sont calculés deux fois.
+
+\end{document}