Discussion :

[DrWaste]

Bonjour,

J'ai un code utilisant xarray pour extraire des données depuis une source netCDF (ECMWF ERA5).
De manière simplifiée, il y a quatre dimensions (x, y, z et t) et trois variables (r, h et g).
Toujours en simplifiant un peu, je dois extraire toutes les valeurs de r, h et g pour toutes les dimensions z et t pour chaque couple x/y.
J'ai un code série qui marche, dont une représentation simplifiée est :

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
 
import xarray as xr
 
def extract_data(data,z,x,y,t):
  out_str = ''
  for t_val in t:
    for z_val in z:
      out_str += data.sel(x,y,t_val,z_val).data
  return(out_str)
 
def data_get_write(dir,r,h,g,x,y,z,t):
  r_str = extract_data(r,x,y,z,t)
  with open(dir+'r_file', 'w') as f_r:
    f_r.write(r_str)
  h_str = extract_data(h,x,y,z,t)
  with open(dir+'h_file', 'w') as f_h:
    f_h.write(h_str)
  g_str = extract_data(g,x,y,z,t)
  with open(dir+'g_file', 'w') as f_g:
    f_g.write(g_str)
  return(<operation>)
 
if __name__ == '__main__':
  ds=xr.open_dataset('data.nc')
  r=ds['r']
  h=ds['h']
  g=ds['g']
  z=ds['z']
  t=ds['t']
  dir_list=[]
  x_list=[]
  y_list=[]
  for val_x in x:
    for val_y in y:
      dir=(<operation>)
      dir_list.append(dir)
      x_list.append(x)
      y_list.append(y)
  list_len=len(dir_list)
  r_list=[r]*list_len
  h_list=[h]*list_len
  g_list=[g]*list_len
  z_list=[z]*list_len
  t_list=[t]*list_len
 
  results = map(data_get_write, r_list, h_list, g_list, x_list, y_list, z_list, t_list)

Je l'ai écrit de cette manière (créer des listes pour chaque paramètre puis les passer à la fonction via map) dans l'objectif de paralléliser l'exécution en utilisant concurrent.futures.ThreadPoolExecutor ou concurrent.futures.ProcessPoolExecutor en remplaçant l'appel map de cette manière :

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
4
 
    max_workers = min(32, int(arguments["-n"]) + 4)
    with concurrent.futures.ProcessPoolExecutor(max_workers=max_workers) as executor:
      results = executor.map(data_get_write, r_list, h_list, g_list, x_list, y_list, z_list, t_list)

Avec ThreadPoolExecutor le problème c'est que la version parallélisée est moins rapide que la version série...
Sur mon jeu de données de test (netcdf de 600Mo, dimensions x et y avec une vingtaine de valeurs et t à 6 valeurs) et mon serveur de test, la version série prend environ 2'40", la version parallèle 6'30". J'ai fait varier le nombre de cœurs alloués de 1 à 20, il n'y a aucun impact sur le temps d'exécution (ou du moins juste à la marge, plus ou moins quelques secondes, moins de 10%). Il n'y a pas plus de variabilité entre les runs en faisant changer le nombre de cœurs alloués qu'en relançant plusieurs fois de suite le même run.
Ça s'explique, si je ne m'abuse, par le verrou global qui empêche les threads concurrents d'accéder aux données en lecture de manière réellement concurrente puisqu'elles sont passées par référence et pas par valeur. Et ça semble être confirmé en basculant sur ProcessPoolExecutor même si la baisse de vitesse d'exécution est surprenante.

Avec ProcessPoolExecutor on constate un gain, ça passe à 1'25" mais à nouveau, quel que soit le nombre de cœurs, il n'y pas de modification significative du temps d'exécution. On pourrait penser qu'il s'agit d'un bottleneck sur le filesystem mais il n'y a pas non plus de modification en fonction du système de stockage utilisé, que ce soit le disque local au serveur ou un montage NFS.

Il doit donc y avoir quelque chose de fondamental, soit dans mon choix de bibliothèque de parallélisation ou dans mon code même, qui l'empêche de scaler avec la surface d'exécution. C'est peut-être très simple mais je suis complètement novice, autant en lecture / extraction de données depuis du netCDF qu'en parallélisation de code, donc il ne serait pas surprenant que ce soit un cas de pebkac. Toute aide / suggestion de piste pour améliorer mon code est la bienvenue, merci !

[fred1599]

Bonsoir,

Il ne serait pas déjà possible d'optimiser votre algorithme pour éviter les boucles (en vectorisant) ?

Avant de faire une version parallélisée, il faudra déjà optimiser la version non paralléliser, vérifier si les temps d'exécution conviennent, déterminer les goulots d'étranglement puis seulement à partir de là, paralléliser les tâches qui prennent du temps.

[wiztricks]

Salut,

La loi de Amdhal s'applique (à vous de faire les calculs pour savoir si vos attentes sont réalistes).
Passer de 2mn40 à 1mn25, n'est pas si mal.

- W

[DrWaste]

Bonjour,

Merci pour la suggestion de vectorisation. Effectivement, pour la création des listes initiales qui seront passées à map ça aurait du sens et ce serait plus élégant, à l'occasion je le ferai.

Et aussi pour la suggestion d'améliorer le code hors de l'aspect parallélisation. Couplé à d'autres discussions que j'ai eues, ça m'a motivé pour ajouter des perf_counter() partout dans mon code et regarder le temps d'exécution de chaque étape. Si il y a un conseil que je peux donner, c'est celui là : si on veut identifier ce qui prend du temps dans un script de ce type, des perf_counter() partout avec des print de différence à la fin c'est vraiment assez parlant.

J'ai commencé par le main et on voit assez bien le coupable (en abscisses les étapes et en ordonnées le temps d'exécution cumulatif) :



Toutes les boucles au début sont epsilonesques (d'où le fait que je ne vais pas me précipiter pour les vectoriser autrement que pour l'élégance).
L'étape qui contribue la quasi totalité du temps c'est qui correspond à

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

results = executor.map()

.

Je suis donc allé mettre des perf_counter() de plus en plus profond dans les fonctions appelées et ce qui prend vraiment du temps c'est l'extraction des données depuis le tableau xarray qui est de cette forme :

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

dataset.sel(x=var1,y=var2)

.

Je me suis demandé si on ne pourrait pas accélérer cette étape et je suis tombé sur une page pas très rassurante. L'issue est marquée résolue mais en lisant les commentaires c'est parce qu'il n'y a pas de solution, pas parce que le code a été amélioré...

J'avais choisi d'utiliser xarray plutôt que netCDF4 pour avoir une dépendance plus générique mais je me suis demandé si ce ne serait pas plus rapide avec netCDF4. Il est donc temps de faire un petit test, j'ai importé timeit puis :

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
 
>>> code_seg='''\
... import netCDF4 as nc
... dn = nc.Dataset('ERA/data.nc')
... zn=dn['z'][:]
... for i in range(10000):
...     res=zn[0,:,0,0]
... '''
>>> timeit.timeit(code_seg, number=10**3)
53.229370541870594
 
>>> code_seg='''\
... import xarray as xr
... dn = xr.open_dataset('ERA/data.nc')
... zn=dn['z']
... for i in range(10000):
...     res=zn.sel(valid_time='2019-08-03T06:00:00.000000000',latitude='36.25',longitude='351')
... '''
>>> timeit.timeit(code_seg, number=10**3)

À l'heure où j'écris ces lignes, j'attends encore le résultat avec xarray, ça fait plus d'une heure que ça tourne...

Je crois que j'ai une piste pour améliorer mon code et je la dois au moins en partie à vos suggestions, merci beaucoup

[wiztricks]

Je crois que j'ai une piste pour améliorer mon code et je la dois au moins en partie à vos suggestions, merci beaucoup

A tout hasard, si on cherche avec les mots clefs "parallel netcdf python" on tombe sur des bibliothèques comme celle ci.

- W

[DrWaste]

Hello,

N'ayant aucune donnée, difficile d'affirmer ni vérifier ce que vous dîtes, moins encore vous proposer des solutions qui fonctionnent du premier coup, étant donné qu'on ne peut pas tester.

On ne connaît même pas les longueurs de t et z. Sûre que si elles sont de petites tailles, l'impact côté vectorisation est négligeable.

En plus vous faîtes de la concaténation non optimisée (utiliser les listes, append et join pour améliorer).

Mais bon encore une fois, sans donnée de test, difficile de vous aider plus...

Les longueurs de t e z sont effectivement assez réduites. Environ max 800 pour t actuellement, habituellement plutôt aux alentours de 10. Et z est à environ 30.

Un jeu de donnée ressemble à cette pièce jointe : data.zip. Le fichier netcdf ici est tout petit par rapport aux jeux tests qu'on utilise habituellement, t = 1, x et y (lat et lon) à moins de 10, un exemple typique est un à deux ordres de grandeur supérieurs sur x et y.

Ça permet de voir que ce qui prend vraiment du temps, de très très loin, c'est le select des données, comme sur mon graphique précédent.

A tout hasard, si on cherche avec les mots clefs "parallel netcdf python" on tombe sur des bibliothèques comme celle ci.

- W

Oui, merci en effet, et je l'envisagerai sur un autre projet qui utilise ces même données mais à un ordre de grandeur supérieur où la parallélisation de la lecture en input commencera à avoir un intérêt. Sur le projet en cours c'est négligeable, ce qui impacte vraiment c'est le select des données dans la structure en mémoire.

D'ailleurs, à cet effet, j'ai écrit un petit script tout bête qui utilise timeit pour comparer le temps de différentes méthodes. Actuellement je charge les données dans un xarray et je fais un select en utilisant les étiquettes des dimensions. Ce script compare cette méthode à l'extraction du tableau de la variable seule, en utilisant netCDF4 et xarray, sur la structure netcdf jointe plus haut demo-timeit.zip, il lui faut les bibliothèques timeit, xarray et netCDF4. On voit que l'extraction des données en utilisant les étiquettes est beaucoup beaucoup plus lente que la lecture directe dans le tableau de données, l'écart est assez saisissant. Avec xarray on a un facteur d'environ 300x ! Le problème c'est qu'extraire les données en utilisant le tableau présume qu'on connaît l'ordre des dimensions dans l'imbrication des tableaux, dans l'exemple donnée c'est 1) t, 2) p, 3) x et 4) y. Si cet ordre change je dois changer l'ordre des arguments passés dans le tableau res=zn[0,:,0,0]. Et comme je boucle sur t et z, je ne peux pas tester et construire cet ordre avant d'être sur cette boucle assez interne de mon code, ou alors prévoir les 16 combinaisons possibles et selon l'ordre des dimensions appeler la fonction correspondante mais ça veut dire autant de fonctions qui font quasiment exactement la même chose. J'ai l'impression qu'une méthode plus élégante / efficace doit exister mais je bute encore sur ce point. Pour l'instant je teste si l'ordre est celui que j'ai défini par défaut et je génère une erreur critique mais ça crée une fragilité dans mon code.

[fred1599]

N'ayant aucune donnée, difficile d'affirmer ni vérifier ce que vous dîtes, moins encore vous proposer des solutions qui fonctionnent du premier coup, étant donné qu'on ne peut pas tester.

On ne connaît même pas les longueurs de t et z. Sûre que si elles sont de petites tailles, l'impact côté vectorisation est négligeable.

En plus vous faîtes de la concaténation non optimisée (utiliser les listes, append et join pour améliorer).

Mais bon encore une fois, sans donnée de test, difficile de vous aider plus...