Skip to content
Snippets Groups Projects
Commit 8aaeb033 authored by Lionel GUEZ's avatar Lionel GUEZ
Browse files

Starting to break everything, targetting two separate programs for 

extraction of eddies and overlapping.

Procedure dispatch_snapshot will be called by the overlap program and
will not write eddies. So remove the calls to write_eddy from
dispatch_snapshot and move them to test_get_snapshot and
test_read_snapshot.

Do not create files isolated_nodes_1.csv and number_eddies_1.csv in
test_get_snapshot and test_read_snapshot, that will be done by the
overlap program.
parent cd7598f2
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
Documentation_texfol/documentation.tex
+ 787
− 754
View file @ 8aaeb033
......@@ -50,88 +50,39 @@
\maketitle
\tableofcontents
\section{Divers}
\section{Extraction des tourbillons}
\subsection{Divers}
Chaque ligne de niveau de la hauteur correspond à une ligne de courant
de la vitesse.
La zone intérieure à la ligne de vitesse maximale est le c\oe{}ur du
tourbillon. On cherche les recouvrements entre ces c\oe{}urs, l'idée
étant que 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 pathologiques. out\_cont est vide si et seulement s'il n'y a pas
de bon contour à innermost\_level. Le champ speed de extremum.dbf vaut
missing\_speed si et seulement si radius4 == 0 ou speed\_outerm dans
set\_max\_speed est NaN. Si radius4 $\ge$ 2 et speed dans extremum.dbf
vaut missing\_speed alors speed\_outerm dans set\_max\_speed est NaN,
e\%speed\_cont est non vide dans max\_speed\_contour.shp et
mean\_speed sur e\%speed\_cont a donné un NaN. En d'autres termes, les
deux appels à mean\_speed ont tous les deux produit un
NaN. (e\%speed\_cont est non vide dans max\_speed\_contour.shp et
speed dans extremum.dbf vaut missing\_speed) si et seulement si
(radius4 $\ge$ 2 et speed dans extremum.dbf vaut missing\_speed).
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}).
But du programme : produire une base de données de tourbillons. Le
programme ne traite qu'une date.
Cas pathologiques. \verb+out_cont+ est vide si et seulement s'il n'y a
pas de bon contour à \verb+innermost_level+. Le champ speed de
extremum.dbf vaut \verb+missing_speed+ si et seulement si radius4 == 0
ou \verb+speed_outerm+ dans \verb+set_max_speed+ est NaN. Si radius4
$\ge$ 2 et speed dans extremum.dbf vaut \verb+missing_speed+ alors
\verb+speed_outerm+ dans \verb+set_max_speed+ est NaN,
\verb+e%speed_cont+ est
non vide dans \verb+max_speed_contour.shp+ et \verb+mean_speed+ sur
\verb+e%speed_cont+ a donné un NaN. En d'autres termes, les deux
appels à \verb+mean_speed+ ont tous les deux produit un
NaN. (\verb+e%speed_cont+ est
non vide dans \verb+max_speed_contour.shp+ et speed dans extremum.dbf
vaut \verb+missing_speed+) si et seulement si (radius4 $\ge$ 2 et
speed dans extremum.dbf vaut \verb+missing_speed+).
Cas où les composantes \verb+out_cont+ ou \verb+speed_cont+ d'un
tourbillon sont vides, c'est-à-dire égales à
\verb+null_ssh_contour()+. On peut distinguer les cas du tableau
(\ref{tab:null_ssh_contour}).
\begin{table}[htbp]
\centering
\begin{tabular}{lllll|l}
out\_cont & valid & radius4 & max\_speed
\verb+out_cont+ & valid & radius4 & \verb+max_speed+
& fraction & note \\
\hline
null & F & 0 & \np{e4} & 5 \% & 1 \\
......@@ -152,11 +103,11 @@ les cas du tableau (\ref{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
12, 15, 7, 122, 54. Au total : 210 cas où les composantes \verb+out_cont+ ou
\verb+speed_cont+ d'un tourbillon sont vides, sur 260 tourbillons, soit 81
\% des tourbillons.
\section{Calcul de la vitesse azimutale par rapport à l'extremum}
\subsection{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
......@@ -225,8 +176,7 @@ distance entre l'extremum et $M$, la vitesse azimutale est :
V_\theta = \frac{x v - y u}{\rho}
\end{equation*}
\section{Bons contours}
\subsection{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
......@@ -247,7 +197,7 @@ 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.
dans \verb+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
......@@ -304,7 +254,7 @@ 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}
\subsection{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
......@@ -333,8 +283,8 @@ 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.
impactées sont \verb+get_1_outerm+, \verb+local_extrema+, inside et
\verb+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
......@@ -348,28 +298,25 @@ 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
Il faut prévoir dans \verb+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
parler, la valeur de \verb+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
intiale de \verb+level_good+ prendre dans \verb+get_1_outerm+. Je garde
néanmoins \verb+innermost_level+ comme valeur intiale de \verb+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
prenant une valeur initiale de \verb+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}
\subsection{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. Les
tourbillons interpolés sont écrits avec les autres, en utilisant une
forme NULL pour le polygone correspondant à un contour interpolé.
séparer les contours et les positions des extremums.
\begin{itemize}
\item \verb+extremum_$m.shp+ : points
\item \verb+extremum_$m.dbf+
......@@ -378,7 +325,7 @@ forme NULL pour le polygone correspondant à un contour interpolé.
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
date, indice de tourbillon à cette date, 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
......@@ -400,10 +347,7 @@ L'indice de date et l'indice de tourbillon dans
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é.
préférable pour la solidité et le confort d'utilisation.
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
......@@ -417,56 +361,9 @@ 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 ?
Nous stockons 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, en
comptant 10 caractères par entier, 13 caractères par réel, \np{5e7}
sommets et deux arcs par sommet en moyenne, 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}
\subsection{Structures de données}
On suppose que le maillage est le même à toutes les dates et qu'il est
régulier.
On suppose que le maillage est régulier.
Type dérivé polyline : défini dans \verb+Contour_531+. Mémoire
utilisée par un scalaire de ce type :
......@@ -493,19 +390,16 @@ calculer une vitesse.
Nous avons besoin que la composante extr\_map de snapshot soit étendue
en longitude dans le cas d'un domaine périodique pour son utilisation
par set\_max\_speed et successive\_overlap.
par \verb+set_max_speed+.
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
Les tailles ci-dessous sont calculées pour 8000 tourbillons, 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
\item[eddy] La comparaison entre \verb+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. Mémoire utilisée par
......@@ -539,10 +433,6 @@ soit environ \np{0.5} KiB.
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*}
......@@ -563,147 +453,724 @@ 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\}}
\subsection{Résumé du problème de géométrie algorithmique}
\STATE \COMMENT{prologue}
\begin{algorithmic}
\STATE recherche de tous les extremums locaux
\STATE entrer(longitude, latitude)
\STATE \COMMENT{longitude et latitude dans l'ordre croissant}
\FOR{k = 1 \TO max\_delta + 1}
\FOR{chaque extremum local}
\STATE appel de get\_snapshot(flow(k), k)
\STATE recherche du \og bon contour le plus extérieur\fg{}
\ENDFOR
\end{algorithmic}
\ENDFOR
\FOR{k = 2 \TO max\_delta + 1}
\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.
\STATE appel de successive\_overlap(k, k, flow)
\ENDFOR
\FOR{delta = 2 \TO max\_delta}
Recherche du bon contour le plus extérieur : par dichotomie.
\begin{algorithmic}
\STATE $a_1 = 0$
\FOR{k = delta + 1 \TO max\_delta + 1}
\STATE $a_2$ = max ou min de SSH sur le domaine - SSH de l'extremum cible
\STATE appel de non\_successive\_overlap(k, k, delta, flow)
\WHILE{$a_2 - a_1 >$ précision}
\ENDFOR
\ENDFOR
\STATE recherche du bon contour à l'amplitude
$\frac{a_1 + a_2}{2}$
\STATE \COMMENT{boucle principale}
\IF{bon contour trouvé}
\FOR{k = max\_delta + 2 \TO n\_dates}
\STATE $a_1 = \frac{a_1 + a_2}{2}$
\STATE appel de dispatch\_snapshot(flow(1), k - max\_delta - 1)
\STATE flow(:max\_delta) = flow(2:)
\ELSE
\STATE appel de get\_snapshot(flow(max\_delta + 1), k)
\STATE $a_2 = \frac{a_1 + a_2}{2}$
\STATE appel de successive\_overlap(max\_delta + 1, k, flow)
\ENDIF
\FOR{delta = 2 \TO max\_delta}
\ENDWHILE
\end{algorithmic}
\STATE appel de non\_successive\_overlap(max\_delta + 1, k, delta,
flow)
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).
\ENDFOR
\ENDFOR
\STATE \COMMENT{épilogue}
Les procédures géométriques utilisées :
\begin{itemize}
\item polygon\_area\_2d (bibliothèque Geometry)
\item polygon\_contains\_point (bibliothèque Geometry)
\item find\_contours\_reg\_grid (bibliothèque Contour\_531)
\end{itemize}
travaillent dans le plan
et non sur la sphère. Ces procédures sont appelées par :
\begin{itemize}
\item spherical\_polyline\_area
\item inside\_4
\item good\_contour
\end{itemize}
Dans le plan, par deux points passe un unique segment, et la donnée
d'un triangle est équivalente à la donnée de l'ensemble de ses trois
sommets. Sur la sphère, par deux points passent deux arcs de grand
cercle, et plusieurs triangles peuvent être associés à un ensemble de
trois sommets. Par ailleurs, les procédures géométriques dans le plan
ne peuvent pas admettre en arguments d'entrée des longitudes définies
modulo $2 \pi$. Pour que les valeurs de longitudes en entrée des
procédures géométriques soient contraintes sans ambiguïté, nous
imposons que tout domaine de longitudes en entrée des procédures
géométriques soit de longueur strictement inférieure à $\pi$.
\begin{equation}
\label{eq:length_pi}
\forall (\lambda_1, \lambda_2) \in \mathbb{R}^2, \forall k \in \mathbb{Z},
|\lambda_2 - \lambda_1| < \pi
\Rightarrow |\lambda_2 + 2 k \pi - \lambda_1| \ge \pi
\end{equation}
(Toute modulation de $2 \pi$ d'une longitude viole la contrainte.) Les
procédures qui appellent les procédures géométriques doivent se
charger d'ajouter le multiple de $2 \pi$ pour ramener les longitudes
dans un intervalle de longueur strictement inférieure à $\pi$. La
contrainte de domaine de longitudes s'applique aux procédures
\verb+get_1_outerm+ et \verb+set_max_speed+.
\FOR{j = 1 \TO max\_delta + 1}
\subsection{Sous-algorithmes}
\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}
\subsubsection{get\_1\_outerm}
\section{Décomposition de domaine}
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)
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}
\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 [\verb+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.
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.
Pour l'amplitude minimale, si je prends simplement l'amplitude
correspondant à \verb+innermost_level+ alors l'algorithme échoue lorsque
l'extremum est à un point du bord et que le point à \verb+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.
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.
\subsubsection{get\_snapshot}
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.
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.
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}$).
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.
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.
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 \verb+outermost_possible_level+
n'est pas un critère suffisant, parce que \verb+outermost_possible_level+
ne teste que la direction est. Il faudrait tester après l'appel de
\verb+get_1_outerm+ si le cadre le plus extérieur touche le bord du
domaine.
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
\subsubsection{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.
\subsubsection{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 c'est plus clair et plus simple sans ce
codage. J'ajoute donc un argument vecteur cyclone à \verb+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).
Cf. figure (\ref{fig:periodicity}).
\begin{figure}[htbp]
\centering
\includegraphics{periodicity}
\caption{Cas de périodicité dans local\_extrema. Le tableau
représenté est extr\_map. Pour définir la case (2, j), il faut
avoir copié la case (m - 1, j - 1) (qui peut être non nulle). Pour
définir la case (m - 1, j), il faut avoir copié la case (2, j)
(qui peut être non nulle). (La case (2, j - 1) doit aussi avoir
été copiée mais cela a été fait au passage sur la colonne
précédente.)}
\label{fig:periodicity}
\end{figure}
Cas de périodicité : field(i + m - 2, :) = field(i, :).
\subsubsection{max\_speed\_contour\_ssh}
Il manque des tests pour lesquels direction vaut 1 et 4.
\subsubsection{set\_all\_outerm}
Notons ici :
\begin{align*}
& u = \mathtt{urc(1)} \\
& l = \mathtt{llc(1)} \\
& s = \mathtt{step(1)} \\
& e = \mathtt{s\%list\_vis(i)\%coord\_extr(1)} \\
& w = \mathtt{corner\_window(1)} \\
& m = \mathtt{max\_radius(1)}
\end{align*}
Pour que la contrainte de domaine de longitudes dans \verb+get_1_outerm+
soit satisfaite, il faut que, dans set\_all\_outerm, au moment de
l'appel à \verb+get_1_outerm+ :
\begin{equation}
\label{eq:max_radius}
(u - l) s < \pi
\end{equation}
et que $e$ et outside\_points(1, :) soient dans l'intervalle
$[w, w + (u - l) s]$. Pour que l'inégalité (\ref{eq:max_radius}) soit
satisfaite, il suffit que :
\begin{equation*}
2 m s < \pi
\end{equation*}
ce que nous imposons en amont de l'appel de set\_all\_outerm. Si
periodic est vrai alors :
\begin{equation*}
[w, w + (u - l) s] = [e - m s, e + m s]
\end{equation*}
donc $e$ appartient bien à cet intervalle. Si periodic est faux,
$[w, w + (u - l) s]$ peut être plus réduit que $[e - m s, e + m s]$
mais nous savons aussi que $e$ est dans l'intervalle [corner(1),
corner(1) + (nlon - 1) s], donc $e$ est bien toujours dans
$[w, w + (u - l) s]$.
Il reste le problème de outside\_points(1, :). outside\_points(1, :)
contient des coordonnées prises dans s\%list\_vis\%coord\_extr(1). Si
le domaine est périodique et si l'extremum cible est proche en
longitude de corner(1) par exemple alors outside\_points(1, :) peut
contenir des longitudes proches de corner(1) + (nlon - 1) s. Il faut
leur retrancher $2 \pi$. Nous posons le problème plus généralement de
la façon suivante. Soit $\lambda \in \mathbb{R}$. On cherche $k \in
\mathbb{Z}$ tel que :
\begin{equation*}
\lambda + 2 k \pi \in [w, w + (u - l) s]
\end{equation*}
Supposons que $k$ existe (c'est le cas si $\lambda$ est une longitude
retournée par nearby\_extr). Alors $k$ est unique (cf. propriété
(\ref{eq:length_pi})). On a :
\begin{equation*}
\lambda + 2 k \pi \in [w, w + 2 \pi[
\end{equation*}
C'est-à-dire :
\begin{equation*}
k - 1 < \frac{w - \lambda}{2 \pi} \le k
\end{equation*}
Ou encore :
\begin{equation*}
k = \left\lceil \frac{w - \lambda}{2 \pi} \right\rceil
\end{equation*}
\subsubsection{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 \verb+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 \verb+get_snpashot+. Cela remettrait en
question tout l'algorithme.
\subsubsection{spherical\_polyline\_area}
La procédure de calcul des contours donne une suite de sommets. Les
courbes reliant ces sommets ne sont en fait pas imposés
intrinsèquement par cette procédure. Par exemple, il n'est ni plus ni
moins juste de supposer que ces courbes sont des portions de grands
cercles ou de supposer qu'elles donnent des segments en projection
plate carrée. Néanmoins, dans le reste du programme, pour décider si
un contour entoure un extremum, ou pour calculer l'intersection de
deux tourbillons, nous supposons que les contours sont des polygones
en $(\lambda, \phi)$ (c'est-à-dire que les courbes reliant les sommets
donnent des segments en projection plate carrée).
Un domaine $S$ quelconque à la surface de la Terre a pour aire :
\begin{equation*}
a^2 \iint_S \cos \phi\ \ud \phi\ \ud \lambda
= a^2 \iint_S \ud \lambda\ \ud(\sin \phi)
\end{equation*}
Nous avons une procédure qui calcule la surface $\iint \ud x\ \ud y$
d'un polygone en $(x, y)$. Nous calculons la surface du polygone en
$(\lambda, \sin \phi)$ passant par les sommets du
contour. C'est-à-dire que nous supposons que les courbes reliant les
somments donnent des segments en projection conforme cylindrique.
Rémi fait un autre calcul : la surface du polygone en
$(\lambda \cos \phi, \phi)$ passant par les sommets du contour. Le
changement de variable :
\begin{align*}
& x = \lambda \cos \phi \\
& y = \phi
\end{align*}
a pour jacobien $\cos \phi$ donc on a bien :
\begin{equation*}
\iint \cos \phi\ \ud \phi\ \ud \lambda = \iint \ud x\ \ud y
\end{equation*}
Les deux surfaces obtenues, celle du polygone en
$(\lambda, \sin \phi)$ et celle du polygone en
$(\lambda \cos \phi, \phi)$ sont très proches si les sommets sont très
proches. C'est le cas puisque les sommets sont séparés au maximum de
\np{0.25}° en longitude et latitude. On peut aussi le voir en
considérant la formule de calcul de surface du polygone :
\begin{align*}
S & = \frac{1}{2} \sum_i x_i (y_{i + 1} - y_{i - 1}) \\
& =
\frac{1}{2} \sum_i \lambda_i \cos \phi_i (\phi_{i + 1} - \phi_{i - 1}) \\
& \approx \frac{1}{2} \sum_i \lambda_i (\sin \phi_{i + 1} - \sin \phi_{i - 1})
\end{align*}
si les valeurs successives des $\phi_i$ sont proches.
\subsubsection{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{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}. 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 pour le test \verb+Get_1_outerm+ 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{Test Set\_max\_speed. 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{Test Get\_snapshot\_region\_1. 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 (valid faux), dont 5 sur les bords et 11
dégénérés. Avec minimum d'amplitude : 34 extremums pour lesquels
valid est faux. Je vérifie que, pour chaque extremum, si valid
est faux dans le cas sans minimum, valid reste faux dans le cas
avec minimum. Il y a donc 20 extremums qui passent à valid 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.
Sur la région 3, avec un mm de minimum d'amplitude, le contour
extérieur 168 a un diamètre de 27 points, environ 7 degrés en
longitude. \`A 40° S, donc un diamètre de 600 km environ. On ne peut
donc pas prendre max\_radius = 12. Ayant fait le test avec max\_radius
= 12, je constate aussi que le temps d'exécution est très sensible à
max\_radius : un facteur \np{2.7} environ entre max\_radius = 12 et
max\_radius = 20. Le temps d'exécution semble donc proportionnel à
max\_radius$^2$, ce qui peut se comprendre.
\section{Recouvrements}
\subsection{Divers}
La zone intérieure à la ligne de vitesse maximale est le c\oe{}ur du
tourbillon. On cherche les recouvrements entre ces c\oe{}urs, l'idée
étant que les extremums à deux dates correspondent au même tourbillon
si les c\oe{}urs se recouvrent.
But du programme : produire 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 (\verb+max_delta+ + 1) dates nécessaires à la
profondeur de vue.
Lorsqu'on compare des tourbillons à distance temporelle delta (entre 1
et \verb+max_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.
On pourrait écrire l'algorithme principal en explicitant les
\verb+get_snpashot+, en explicitant le cas m < \verb+n_proc+ et en supprimant la
variable \verb+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 recherche de recouvrements.}
\label{fig:call_graph}
\end{figure}
Les procédures géométriques utilisées :
\begin{itemize}
\item polygon\_area\_2d (bibliothèque Geometry)
\item gpc\_polygon\_clip\_f (bibliothèque GPC\_F)
\end{itemize}
Ces procédures sont appelées par :
\begin{itemize}
\item spherical\_polyline\_area
\item \verb+successive_overlap+
\end{itemize}
Cf. discussion autour de l'équation~(\ref{eq:length_pi}).
\subsection{Entrées et sorties}
On n'enregistre pas de contours pour les tourbillons interpolés. Les
tourbillons interpolés sont écrits avec les autres, en utilisant une
forme NULL pour le polygone correspondant à un contour interpolé.
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é.
Nous stockons 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, en
comptant 10 caractères par entier, 13 caractères par réel, \np{5e7}
sommets et deux arcs par sommet en moyenne, 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), \verb+get_snapshot+ (lecture de
SSH, u, v), \verb+successive_overlap+ (écriture dans
\verb+edgelist_$m.txt+, arcs entre \verb+k_begin+(m) et
\verb+k_end_main_loop+(m)), \verb+non_successive_overlap+ directement
(écriture dans \verb+edgelist_$m.txt+), \verb+write_eddy+ (écriture
dans \verb+extremum_$m+, \verb+outermost_contour_$m+,
\verb+max_speed_contour_$m+), \verb+dispatch_snapshot+ directement
(écriture dans \verb+isolated_nodes_$m+,
\verb+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 de recouvrement.}
\label{fig:input_output}
\end{figure}
\subsection{Déclarations pour l'algorithme principal}
On suppose que le maillage est le même à toutes les dates et qu'il est
régulier.
Nous avons besoin que la composante extr\_map de snapshot soit étendue
en longitude dans le cas d'un domaine périodique pour son utilisation
par \verb+successive_overlap+.
Les tailles ci-dessous sont calculées pour \verb+max_delta+ = 4 et une
durée de 16 ans, soit environ 5000 dates.
\begin{description}
\item[max\_delta] Scalaire entier. Intervalle maximal d'indices de
dates auxquelles on cherche des recouvrements. Normalement 4.
\item[flow] : Vecteur de type snapshot, de taille \verb+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}
i désigne un indice de tourbillon, j un indice de position dans la
fenêtre temporelle (entre 1 et \verb+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 appel de get\_snapshot(flow(k), 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 appel de get\_snapshot(flow(max\_delta + 1), 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}
\subsection{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 < \verb+n_proc+ :
\begin{equation}
\label{eq:stitching}
\mathtt{k\_end}(m)
= \mathtt{k\_begin}(m + 1) + \mathtt{max\_delta} - 1
\end{equation}
k\_end - \verb+k_begin+ + 1 appels à \verb+get_snapshot+, pour les dates \verb+k_begin+
à k\_end. k\_end - \verb+k_begin+ + 1 appels à \verb+dispatch_snapshot+, pour les
dates \verb+k_begin+ à k\_end. Si m > 1 alors \verb+max_delta+ envois. Si m <
\verb+n_proc+ alors \verb+max_delta+ réceptions.
Nombre d'appels à \verb+get_snapshot+ avec lecture. Si m < \verb+n_proc+ : k\_end -
\verb+k_begin+ - \verb+max_delta+ + 1. Si m = \verb+n_proc+ : k\_end - \verb+k_begin+ + 1.
\verb+successive_overlap+ est appelé pour tous les indices k compris entre
\verb+k_begin+ + 1 et \verb+k_end_main_loop+, seulement. Nombre d'appels à
\verb+successive_overlap+ : \verb+k_end_main_loop+ - \verb+k_begin+. Donc si m <
\verb+n_proc+ : k\_end - \verb+k_begin+ - \verb+max_delta+ + 1 appels. Si m = \verb+n_proc+ : k\_end
- \verb+k_begin+ appels.
Nombre d'appels à \verb+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 : (\verb+max_delta+ - 1)(\verb+k_end_main_loop+ -
\verb+k_begin+ - \verb+max_delta+). Dans l'épilogue, si m < \verb+n_proc+ :
$\frac{\mathtt{max\_delta} (\mathtt{max\_delta} - 1)}{2}$. Total si m
< \verb+n_proc+ : (\verb+max_delta+ - 1)(k\_end - \verb+k_begin+ - \verb+max_delta+ +
1). Total si m = \verb+n_proc+ : (\verb+max_delta+ - 1)(k\_end - \verb+k_begin+ -
$\frac{\mathtt{max\_delta}}{2}$).
Si k\_end - \verb+k_begin+ était le même pour m = \verb+n_proc+ et m <
\verb+n_proc+, le processus m=\verb+n_proc+ aurait des plus grands nombres
d'appels aux trois procédures. Notons $\Delta k$ = k\_end - \verb+k_begin+,
considéré comme une fonction de m. Si, pour m < \verb+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 :
< \verb+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
vaut \verb+max_delta+ ou \verb+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
Pour \verb+n_proc+ $\ne 1$ et m < \verb+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}
......@@ -727,19 +1194,20 @@ Il suffit donc que :
\end{equation*}
L'équation (\ref{eq:stitching}) et l'inégalité (\ref{eq:min_delta_k})
impliquent que, pour m < n\_proc :
impliquent que, pour m < \verb+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
donc les appels de \verb+get_snapshot+ dans le prologue font des
lectures. Les appels de \verb+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.
m < \verb+n_proc+, la boucle principale fait au moins une
itération. Si m < \verb+n_proc+ : tous les appels à
\verb+get_snapshot+ de la boucle principale sauf le dernier font des
lectures ; le dernier appel à \verb+get_snapshot+ de la boucle
principale fait une réception.
Si $m \ge 2$ alors :
\begin{equation*}
......@@ -751,7 +1219,7 @@ Exemples : cf. tableaux (\ref{tab:m2}) et (\ref{tab:m3}) et figure
\begin{table}[htbp]
\centering
\begin{tabular}{llll}
m & k\_begin & k\_end & k\_end\_main\_loop \\
m & \verb+k_begin+ & k\_end & \verb+k_end_main_loop+ \\
\hline
1 & 1 & 9 & 6 \\
2 & 6 & 10 & 10
......@@ -765,7 +1233,7 @@ Exemples : cf. tableaux (\ref{tab:m2}) et (\ref{tab:m3}) et figure
\begin{table}[htbp]
\centering
\begin{tabular}{llll}
m & k\_begin & k\_end & k\_end\_main\_loop \\
m & \verb+k_begin+ & k\_end & \verb+k_end_main_loop+ \\
\hline
1 & 1 & 9 & 6 \\
2 & 6 & 14 & 11 \\
......@@ -784,16 +1252,17 @@ Exemples : cf. tableaux (\ref{tab:m2}) et (\ref{tab:m3}) et figure
\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,
Le \verb+dispatch_snapshot+ dans la dernière boucle sur j n'est pas
forcément une écriture. Exemple : \verb+n_proc+ = 2, \verb+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.
Si m = 1 alors \verb+write_eddy+ est appelé pour les dates
\verb+k_begin+ à k\_end. Si m > 1 alors \verb+write_eddy+ est appelé
pour les dates \verb+k_begin+ + \verb+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}
\subsection{Algorithme principal, parallèle}
\begin{algorithmic}[1]
\STATE \COMMENT{\{max\_delta $\ge 1$; n\_dates $\ge$ max\_delta + 1;
......@@ -891,237 +1360,44 @@ une date donnée est bien écrite par un seul processus.
\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$ = max ou min de SSH sur le domaine - SSH de l'extremum cible
\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).
Les procédures géométriques utilisées :
\begin{itemize}
\item polygon\_area\_2d (bibliothèque Geometry)
\item polygon\_contains\_point (bibliothèque Geometry)
\item find\_contours\_reg\_grid (bibliothèque Contour\_531)
\item gpc\_polygon\_clip\_f (bibliothèque GPC\_F)
\end{itemize}
travaillent dans le plan
et non sur la sphère. Ces procédures sont appelées par :
\begin{itemize}
\item spherical\_polyline\_area
\item inside\_4
\item good\_contour
\item successive\_overlap
\end{itemize}
Dans le plan, par deux points passe un unique
segment, et la donnée d'un triangle est équivalente à la donnée de
l'ensemble de ses trois sommets. Sur la sphère, par deux points
passent deux arcs de grand cercle, et plusieurs triangles peuvent être
associés à un ensemble de trois sommets. Par ailleurs, les procédures
géométriques dans le plan ne peuvent pas admettre en arguments
d'entrée des longitudes définies modulo $2 \pi$. Pour que les valeurs
de longitudes en entrée des procédures géométriques soient contraintes
sans ambiguïté, nous imposons que tout domaine de longitudes en entrée
des procédures géométriques soit de longueur strictement inférieure à
$\pi$.
\begin{equation}
\label{eq:length_pi}
\forall (\lambda_1, \lambda_2) \in \mathbb{R}^2, \forall k \in \mathbb{Z},
|\lambda_2 - \lambda_1| < \pi
\Rightarrow |\lambda_2 + 2 k \pi - \lambda_1| \ge \pi
\end{equation}
(Toute modulation de $2 \pi$ d'une longitude viole la contrainte.) Les
procédures qui appellent les procédures géométriques doivent se
charger d'ajouter le multiple de $2 \pi$ pour ramener les longitudes
dans un intervalle de longueur strictement inférieure à $\pi$. La
contrainte de domaine de longitudes s'applique aux procédures
get\_1\_outerm, set\_max\_speed et successive\_overlap.
\section{Sous-algorithmes}
\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{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.
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{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{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 valid à 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{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 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
\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).
Cf. figure (\ref{fig:periodicity}).
\begin{figure}[htbp]
\centering
\includegraphics{periodicity}
\caption{Cas de périodicité dans local\_extrema. Le tableau
représenté est extr\_map. Pour définir la case (2, j), il faut
avoir copié la case (m - 1, j - 1) (qui peut être non nulle). Pour
définir la case (m - 1, j), il faut avoir copié la case (2, j)
(qui peut être non nulle). (La case (2, j - 1) doit aussi avoir
été copiée mais cela a été fait au passage sur la colonne
précédente.)}
\label{fig:periodicity}
\end{figure}
Cas de périodicité : field(i + m - 2, :) = field(i, :).
\subsection{max\_speed\_contour\_ssh}
\STATE appel de dispatch\_snapshot(flow(j), k\_end - max\_delta - 1
+ j)
Il manque des tests pour lesquels direction vaut 1 et 4.
\ENDFOR
\end{algorithmic}
\subsection{Sous-algorithmes}
\subsubsection{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
\subsection{non\_successive\_overlap}
\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 valid à faux pour un tourbillon interpolé,
en cohérence avec un champ \verb+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.
\subsubsection{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,
\verb+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]
......@@ -1200,152 +1476,20 @@ flow(j)\%list\_vis\%delta\_in, flow\%number\_eddies. Cf. algorithme
\label{alg:non_successive_overlapping}
\end{algorithm}
\subsection{set\_all\_outerm}
Notons ici :
\begin{align*}
& u = \mathtt{urc(1)} \\
& l = \mathtt{llc(1)} \\
& s = \mathtt{step(1)} \\
& e = \mathtt{s\%list\_vis(i)\%coord\_extr(1)} \\
& w = \mathtt{corner\_window(1)} \\
& m = \mathtt{max\_radius(1)}
\end{align*}
Pour que la contrainte de domaine de longitudes dans get\_1\_outerm
soit satisfaite, il faut que, dans set\_all\_outerm, au moment de
l'appel à get\_1\_outerm :
\begin{equation}
\label{eq:max_radius}
(u - l) s < \pi
\end{equation}
et que $e$ et outside\_points(1, :) soient dans l'intervalle
$[w, w + (u - l) s]$. Pour que l'inégalité (\ref{eq:max_radius}) soit
satisfaite, il suffit que :
\begin{equation*}
2 m s < \pi
\end{equation*}
ce que nous imposons en amont de l'appel de set\_all\_outerm. Si
periodic est vrai alors :
\begin{equation*}
[w, w + (u - l) s] = [e - m s, e + m s]
\end{equation*}
donc $e$ appartient bien à cet intervalle. Si periodic est faux,
$[w, w + (u - l) s]$ peut être plus réduit que $[e - m s, e + m s]$
mais nous savons aussi que $e$ est dans l'intervalle [corner(1),
corner(1) + (nlon - 1) s], donc $e$ est bien toujours dans
$[w, w + (u - l) s]$.
Il reste le problème de outside\_points(1, :). outside\_points(1, :)
contient des coordonnées prises dans s\%list\_vis\%coord\_extr(1). Si
le domaine est périodique et si l'extremum cible est proche en
longitude de corner(1) par exemple alors outside\_points(1, :) peut
contenir des longitudes proches de corner(1) + (nlon - 1) s. Il faut
leur retrancher $2 \pi$. Nous posons le problème plus généralement de
la façon suivante. Soit $\lambda \in \mathbb{R}$. On cherche $k \in
\mathbb{Z}$ tel que :
\begin{equation*}
\lambda + 2 k \pi \in [w, w + (u - l) s]
\end{equation*}
Supposons que $k$ existe (c'est le cas si $\lambda$ est une longitude
retournée par nearby\_extr). Alors $k$ est unique (cf. propriété
(\ref{eq:length_pi})). On a :
\begin{equation*}
\lambda + 2 k \pi \in [w, w + 2 \pi[
\end{equation*}
C'est-à-dire :
\begin{equation*}
k - 1 < \frac{w - \lambda}{2 \pi} \le k
\end{equation*}
Ou encore :
\begin{equation*}
k = \left\lceil \frac{w - \lambda}{2 \pi} \right\rceil
\end{equation*}
\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{spherical\_polyline\_area}
La procédure de calcul des contours donne une suite de sommets. Les
courbes reliant ces sommets ne sont en fait pas imposés
intrinsèquement par cette procédure. Par exemple, il n'est ni plus ni
moins juste de supposer que ces courbes sont des portions de grands
cercles ou de supposer qu'elles donnent des segments en projection
plate carrée. Néanmoins, dans le reste du programme, pour décider si
un contour entoure un extremum, ou pour calculer l'intersection de
deux tourbillons, nous supposons que les contours sont des polygones
en $(\lambda, \phi)$ (c'est-à-dire que les courbes reliant les sommets
donnent des segments en projection plate carrée).
Un domaine $S$ quelconque à la surface de la Terre a pour aire :
\begin{equation*}
a^2 \iint_S \cos \phi\ \ud \phi\ \ud \lambda
= a^2 \iint_S \ud \lambda\ \ud(\sin \phi)
\end{equation*}
Nous avons une procédure qui calcule la surface $\iint \ud x\ \ud y$
d'un polygone en $(x, y)$. Nous calculons la surface du polygone en
$(\lambda, \sin \phi)$ passant par les sommets du
contour. C'est-à-dire que nous supposons que les courbes reliant les
somments donnent des segments en projection conforme cylindrique.
Rémi fait un autre calcul : la surface du polygone en
$(\lambda \cos \phi, \phi)$ passant par les sommets du contour. Le
changement de variable :
\begin{align*}
& x = \lambda \cos \phi \\
& y = \phi
\end{align*}
a pour jacobien $\cos \phi$ donc on a bien :
\begin{equation*}
\iint \cos \phi\ \ud \phi\ \ud \lambda = \iint \ud x\ \ud y
\end{equation*}
Les deux surfaces obtenues, celle du polygone en
$(\lambda, \sin \phi)$ et celle du polygone en
$(\lambda \cos \phi, \phi)$ sont très proches si les sommets sont très
proches. C'est le cas puisque les sommets sont séparés au maximum de
\np{0.25}° en longitude et latitude. On peut aussi le voir en
considérant la formule de calcul de surface du polygone :
\begin{align*}
S & = \frac{1}{2} \sum_i x_i (y_{i + 1} - y_{i - 1}) \\
& =
\frac{1}{2} \sum_i \lambda_i \cos \phi_i (\phi_{i + 1} - \phi_{i - 1}) \\
& \approx \frac{1}{2} \sum_i \lambda_i (\sin \phi_{i + 1} - \sin \phi_{i - 1})
\end{align*}
si les valeurs successives des $\phi_i$ sont proches.
\subsection{successive\_overlap}
\subsubsection{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.
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.
\verb+out_cont+ dans \verb+speed_cont+ dans \verb+get_snapshot+, on ne
gagnerait rien en quantité de copie, et on perdrait en mémoire vive
occupée.
\subsection{weight}
\subsubsection{weight}
D'autant plus proche de 0 que les tourbillons sont
ressemblants. Comment calculer ce poids ? Comment calculer le poids
......@@ -1355,115 +1499,4 @@ 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.
\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.
\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}. 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 pour le test \verb+Get_1_outerm+ 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{Test Set\_max\_speed. 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{Test Get\_snapshot\_region\_1. 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 (valid faux), dont 5 sur les bords et 11
dégénérés. Avec minimum d'amplitude : 34 extremums pour lesquels
valid est faux. Je vérifie que, pour chaque extremum, si valid
est faux dans le cas sans minimum, valid reste faux dans le cas
avec minimum. Il y a donc 20 extremums qui passent à valid 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.
Sur la région 3, avec un mm de minimum d'amplitude, le contour
extérieur 168 a un diamètre de 27 points, environ 7 degrés en
longitude. \`A 40° S, donc un diamètre de 600 km environ. On ne peut
donc pas prendre max\_radius = 12. Ayant fait le test avec max\_radius
= 12, je constate aussi que le temps d'exécution est très sensible à
max\_radius : un facteur \np{2.7} environ entre max\_radius = 12 et
max\_radius = 20. Le temps d'exécution semble donc proportionnel à
max\_radius$^2$, ce qui peut se comprendre.
\end{document}
Tests/test_get_snapshot.f
+ 12
− 27
View file @ 8aaeb033
......@@ -5,17 +5,16 @@ program test_get_snapshot
use, intrinsic:: ISO_FORTRAN_ENV
! Libraries:
use jumble, only: new_unit
use netcdf, only: nf90_nowrite
use netcdf95, only: nf95_open, find_coord, nf95_inquire_dimension, &
nf95_get_var, nf95_close
use nr_util, only: assert, deg_to_rad, twopi, pi
use shapelib, only: shpfileobject, shpclose
use derived_types, only: snapshot
use dispatch_snapshot_m, only: dispatch_snapshot
use get_snapshot_m, only: get_snapshot
use init_shapefiles_m, only: init_shapefiles
use write_eddy_m, only: write_eddy
implicit none
......@@ -23,8 +22,7 @@ program test_get_snapshot
TYPE(shpfileobject) hshp_extremum ! shapefile extremum_$m
TYPE(shpfileobject) hshp_outermost ! shapefile outermost_contour_$m
TYPE(shpfileobject) hshp_max_speed ! shapefile max_speed_contour_$m
integer unit_isolated, unit_number_eddies
integer i
real:: min_amp = 0.
! minimum amplitude of ssh, between extremum and outermost contour,
......@@ -32,7 +30,7 @@ program test_get_snapshot
real:: min_radius = 20.
! minimum radius of the equal-area disk for an outermost contour, in km
integer ncid, dimid, varid, nlon, nlat
real lon_min, lon_max, lat_min, lat_max ! longitudes and latitudes, in degrees
real step(2) ! longitude and latitude steps, in rad
......@@ -67,7 +65,7 @@ program test_get_snapshot
call nf95_get_var(ncid, varid, lat_max, start = [nlat])
call nf95_close(ncid)
call assert(lon_max > lon_min .and. lat_max > lat_min, &
"test_get_snapshot coordinate order")
step = [(lon_max - lon_min) / (nlon - 1), (lat_max - lat_min) / (nlat - 1)] &
......@@ -85,32 +83,19 @@ program test_get_snapshot
step = step, periodic = periodic)
call init_shapefiles(hshp_extremum, hshp_outermost, hshp_max_speed)
call new_unit(unit_isolated)
open(unit_isolated, file = "isolated_nodes_1.csv", status = "replace", &
action = "write")
write(unit_isolated, fmt = *) '"date index" "eddy index"' ! title line
call new_unit(unit_number_eddies)
open(unit_number_eddies, file = "number_eddies_1.csv", status = "replace", &
action = "write")
write(unit_number_eddies, fmt = *) '"date index" ' &
// '"number of visible extrema" "number of interpolated eddies"'
! (title line)
call dispatch_snapshot(s, hshp_extremum, hshp_outermost, hshp_max_speed, &
unit_isolated, unit_number_eddies, k = 1, m = 1, k_begin = 1, &
max_delta = 4)
! Write snapshot:
do i = 1, s%number_vis_extr
call write_eddy(s%list_vis(i), 1, i, hshp_extremum, hshp_outermost, &
hshp_max_speed)
end do
print *, "Number of extrema:", s%number_vis_extr
CALL shpclose(hshp_extremum)
print *, 'Created shapefile "extremum_1".'
CALL shpclose(hshp_outermost)
print *, 'Created shapefile "outermost_contour_1".'
CALL shpclose(hshp_max_speed)
print *, 'Created shapefile "max_speed_contour_1".'
close(unit_isolated)
print *, 'Created "isolated_nodes_1.csv".'
close(unit_number_eddies)
print *, 'Created "number_eddies_1.csv".'
end program test_get_snapshot
Tests/test_read_snapshot.f
+ 13
− 25
View file @ 8aaeb033
......@@ -2,22 +2,23 @@ program test_read_snapshot
use, intrinsic:: ISO_FORTRAN_ENV
use derived_types, only: snapshot
use dispatch_snapshot_m, only: dispatch_snapshot
use init_shapefiles_m, only: init_shapefiles
use jumble, only: new_unit
! Libraries:
use nr_util, only: deg_to_rad
use read_snapshot_m, only: read_snapshot
use shapelib, only: shpfileobject, shpclose
use shapelib_03, only: shp_open_03
use derived_types, only: snapshot
use init_shapefiles_m, only: init_shapefiles
use read_snapshot_m, only: read_snapshot
use write_eddy_m, only: write_eddy
implicit none
type(snapshot) s
TYPE(shpfileobject) hshp_extremum ! shapefile extremum_1
TYPE(shpfileobject) hshp_outermost ! shapefile outermost_contour_1
TYPE(shpfileobject) hshp_max_speed ! shapefile max_speed_contour_1
integer unit_isolated, unit_number_eddies, i
integer i
real:: corner(2) = [0.125, - 59.875]
! longitude and latitude of the corner of the whole grid, in degrees
......@@ -45,32 +46,19 @@ program test_read_snapshot
call init_shapefiles(hshp_extremum, hshp_outermost, hshp_max_speed)
call new_unit(unit_isolated)
open(unit_isolated, file = "isolated_nodes_1.csv", status = "replace", &
action = "write")
write(unit_isolated, fmt = *) '"date index" "eddy index"' ! title line
call new_unit(unit_number_eddies)
open(unit_number_eddies, file = "number_eddies_1.csv", status = "replace", &
action = "write")
write(unit_number_eddies, fmt = *) '"date index" ' &
// '"number of visible eddies" "number of interpolated eddies"'
! (title line)
call dispatch_snapshot(s, hshp_extremum, hshp_outermost, hshp_max_speed, &
unit_isolated, unit_number_eddies, k = 1, m = 1, k_begin = 1, &
max_delta = 4)
! Write snapshot:
do i = 1, s%number_vis_extr
call write_eddy(s%list_vis(i), 1, i, hshp_extremum, hshp_outermost, &
hshp_max_speed)
end do
print *, "Number of extrema:", s%number_vis_extr
CALL shpclose(hshp_extremum)
print *, 'Created shapefile "extremum_1".'
CALL shpclose(hshp_outermost)
print *, 'Created shapefile "outermost_contour_1".'
CALL shpclose(hshp_max_speed)
print *, 'Created shapefile "max_speed_contour_1".'
close(unit_isolated)
print *, 'Created "isolated_nodes_1.csv".'
close(unit_number_eddies)
print *, 'Created "number_eddies_1.csv".'
print *, "s%ind_extr:"
......
depend.mk
+ 3
− 3
View file @ 8aaeb033
dispatch_snapshot.o : write_eddy.o send_snapshot.o derived_types.o
dispatch_snapshot.o : send_snapshot.o derived_types.o
get_1_outerm.o : spherical_polyline_area.o good_contour.o derived_types.o
get_snapshot.o : set_all_outerm.o set_max_speed.o receive_snapshot.o nearby_extr.o get_var.o derived_types.o
nearby_extr.o : derived_types.o
......@@ -11,11 +11,11 @@ set_max_speed.o : spherical_polyline_area.o mean_speed.o max_speed_contour_ssh.o
spherical_polygon_area.o : spherical_polyline_area.o
successive_overlap.o : weight.o spherical_polyline_area.o spherical_polygon_area.o derived_types.o
test_get_1_outerm.o : get_1_outerm.o derived_types.o
test_get_snapshot.o : init_shapefiles.o get_snapshot.o dispatch_snapshot.o derived_types.o
test_get_snapshot.o : write_eddy.o init_shapefiles.o get_snapshot.o derived_types.o
test_local_extrema.o : local_extrema.o
test_nearby_extr.o : nearby_extr.o derived_types.o
test_read_eddy.o : write_eddy.o read_field_indices.o read_eddy.o init_shapefiles.o derived_types.o
test_read_snapshot.o : read_snapshot.o init_shapefiles.o dispatch_snapshot.o derived_types.o
test_read_snapshot.o : write_eddy.o read_snapshot.o init_shapefiles.o derived_types.o
test_set_all_outerm.o : set_all_outerm.o get_var.o derived_types.o
test_set_max_speed.o : set_max_speed.o derived_types.o
test_spherical_polygon_area.o : spherical_polygon_area.o
......
dispatch_snapshot.f
+ 2
− 19
View file @ 8aaeb033
......@@ -4,27 +4,14 @@ module dispatch_snapshot_m
contains
subroutine dispatch_snapshot(s, hshp_extremum, hshp_outermost, &
hshp_max_speed, unit_isolated, unit_number_eddies, k, m, k_begin, &
max_delta)
! Library:
use shapelib, only: shpfileobject
subroutine dispatch_snapshot(s, unit_isolated, unit_number_eddies, k, m, &
k_begin, max_delta)
use derived_types, only: snapshot
use send_snapshot_m, only: send_snapshot
use write_eddy_m, only: write_eddy
type(snapshot), intent(in):: s
TYPE(shpfileobject), intent(inout):: hshp_extremum ! shapefile extremum_$m
TYPE(shpfileobject), intent(inout):: hshp_outermost
! shapefile outermost_contour_$m
TYPE(shpfileobject), intent(inout):: hshp_max_speed
! shapefile x_speed_contour_$m
integer, intent(in):: unit_isolated
! logical unit for file isolated_nodes_$m.csv
......@@ -40,11 +27,7 @@ contains
!------------------------------------------------------------------
if (m == 1 .or. k >= k_begin + max_delta) then
! Write snapshot:
do i = 1, s%number_vis_extr
call write_eddy(s%list_vis(i), k, i, hshp_extremum, hshp_outermost, &
hshp_max_speed)
if (s%list_vis(i)%delta_in == huge(0) .and. s%list_vis(i)%delta_out &
== huge(0)) write(unit_isolated, fmt = *) k, i
end do
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment