Optimisation/Calculs parallèles d'une simple boucle basée sur Pandas

Question

J'ai cette boucle simple qui traite un grand ensemble de données.

for i in range (len(nbi.LONG_017)):
    StoredOCC = []
    for j in range (len(Final.X)):
        r = haversine(nbi.LONG_017[i], nbi.LAT_016[i], Final.X[j], Final.Y[j])
        if (r < 0.03048):
            SWw = Final.CUM_OCC[j]
            StoredOCC.append(SWw) 
    if len(StoredOCC) != 0:
        nbi.loc[i, 'ADTT_02_2019'] = np.max(StoredOCC)

len(nbi.LONG_017) est de 3000 et len(Final.X) est de 6 millions de points de données.

Je me demandais s'il existe un moyen efficace de mettre en œuvre ce code ou d'utiliser le calcul parallèle pour le rendre plus rapide?

J'ai utilisé le code fourni ici : Formule de Haversine en Python (Roulement et Distance entre deux points GPS) pour ma fonction haversine:

def haversine(lon1, lat1, lon2, lat2):
   """
   Calcule la distance de grand cercle entre deux points 
   sur la terre (spécifié en degrés décimaux)
   """
   # convertir les degrés décimaux en radian 
   lon1, lat1, lon2, lat2 = map(radians, [lon1, lat1, lon2, lat2])

   # formule de haversine 
   dlon = lon2 - lon1 
   dlat = lat2 - lat1 
   a = sin(dlat/2)**2 + cos(lat1) * cos(lat2) * sin(dlon/2)**2
   c = 2 * asin(sqrt(a)) 
   r = 6372.8 # Rayon de la terre en kilomètres. Utilisez 3956 pour les miles
   return r * c

Answer 1

Avant de parler de parallélisation, vous pouvez travailler sur l'optimisation de votre boucle. La première façon serait d'itérer sur les données au lieu des valeurs incrémentales sur la longueur et ensuite accéder aux données à chaque fois :

#exemple jouet
np.random.seed(1)
size_nbi = 20
size_Final = 100

nbi = pd.DataFrame({'LONG_017':np.random.random(size=size_nbi)/100+73, 
                    'LAT_016':np.random.random(size=size_nbi)/100+73,})

Final = pd.DataFrame({'X':np.random.random(size=size_Final)/100+73,
                      'Y':np.random.random(size=size_Final)/100+73, 
                      'CUM_OCC':np.random.randint(size_Final,size=size_Final)})

En utilisant votre méthode, vous obtenez environ 75 ms pour ces tailles de dataframes :

%%timeit
for i in range (len(nbi.LONG_017)):
    StoredOCC = []
    for j in range (len(Final.X)):
        r = haversine(nbi.LONG_017[i], nbi.LAT_016[i], Final.X[j], Final.Y[j])
        if (r < 0.03048):
            SWw = Final.CUM_OCC[j]
            StoredOCC.append(SWw) 
    if len(StoredOCC) != 0:
        nbi.loc[i, 'ADTT_02_2019'] = np.max(StoredOCC)
# 75.6 ms ± 4.05 ms per loop (mean ± std. dev. of 7 runs, 10 loops each)

Maintenant, si vous changez légèrement les boucles, itérer sur les données elles-mêmes et utiliser une liste pour le résultat puis mettre le résultat comme une colonne en dehors des boucles, vous pouvez descendre à 5ms donc 15 fois plus rapide, juste en l'optimisant (et vous pouvez trouver une meilleure optimisation) :

%%timeit
res_ADIT = []
for lon1, lat1 in zip (nbi.LONG_017.to_numpy(),
                       nbi.LAT_016.to_numpy()):
    StoredOCC = []
    for lon2, lat2, SWw in zip(Final.X.to_numpy(), 
                               Final.Y.to_numpy(), 
                               Final.CUM_OCC.to_numpy()):
        r = haversine(lon1, lat1, lon2, lat2)
        if (r < 0.03048):
            StoredOCC.append(SWw) 
    if len(StoredOCC) != 0:
        res_ADIT.append(np.max(StoredOCC))
    else:
        res_ADIT.append(np.nan)
nbi['ADIT_v2'] = res_ADIT
# 5.23 ms ± 305 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)

Maintenant, vous pouvez aller encore plus loin et utiliser numpy et la vectorisation de votre deuxième boucle, et également passer directement les valeurs en radians au lieu de devoir les map à chaque boucle :

# liste vide pour le résultat
res_ADIT = []

# travailler avec des tableaux en radians
arr_lon2 = np.radians(Final.X.to_numpy())
arr_lat2 = np.radians(Final.Y.to_numpy())
arr_OCC = Final.CUM_OCC.to_numpy()

for lon1, lat1 in zip (np.radians(nbi.LONG_017), #pass directly the radians
                       np.radians(nbi.LAT_016):

    # faire toutes les soustractions de manière vectorielle
    arr_dlon = arr_lon2 - lon1
    arr_dlat = arr_lat2 - lat1 

    # pareil ici en utilisant des fonctions numpy
    arr_dist = np.sin(arr_dlat/2)**2 + np.cos(lat1) * np.cos(arr_lat2) * np.sin(arr_dlon/2)**2
    arr_dist = 2 * np.arcsin(np.sqrt(arr_dist)) 
    arr_dist *= 6372.8

    # extraire les valeurs de CUM_OCC qui répondent au critère
    r = arr_OCC[arr_dist<0.03048]

    # vérifier qu'il y a au moins un élément
    if r.size>0:
        res_ADIT.append(max(r))
    else :
        res_ADIT.append(np.nan)

nbi['AUDIT_np'] = res_ADIT

Si vous faites un timeit de ceci, vous descendez à 819 µs ± 15.6 µs per loop (mean ± std. dev. of 7 runs, 1000 loops each) donc environ 90 fois plus rapide que votre solution initiale pour la petite taille des dataframes et tous les résultats sont les mêmes

print(nbi)
     LONG_017    LAT_016  ADTT_02_2019  ADIT_v2  AUDIT_np
0   73.004170  73.008007           NaN      NaN       NaN
1   73.007203  73.009683          30.0     30.0      30.0
2   73.000001  73.003134          14.0     14.0      14.0
3   73.003023  73.006923          82.0     82.0      82.0
4   73.001468  73.008764           NaN      NaN       NaN
5   73.000923  73.008946           NaN      NaN       NaN
6   73.001863  73.000850           NaN      NaN       NaN
7   73.003456  73.000391           NaN      NaN       NaN
8   73.003968  73.001698          21.0     21.0      21.0
9   73.005388  73.008781           NaN      NaN       NaN
10  73.004192  73.000983          93.0     93.0      93.0
11  73.006852  73.004211           NaN      NaN       NaN

Vous pouvez jouer un peu avec le code et augmenter la taille de chaque dataframe jouet, et vous verrez comment en augmentant surtout la taille des données Final, le gain devient intéressant par exemple avec size_nbi = 20 et size_Final = 1000, vous obtenez un gain de 400 fois plus rapide avec la solution vectorisée. et ainsi de suite. Pour votre taille complète de données (3K * 6M), vous avez encore besoin d'un peu de temps et sur mon ordinateur cela prendrait environ 25 minutes selon mon estimation alors que votre solution est en centaines d'heures. Maintenant, si cela n'est pas suffisant, vous pouvez penser à la parallélisation ou également utiliser numba

Answer 2

Cette méthode utilisant Balltree prendra environ une minute sur ma machine (en fonction du rayon, plus il est petit, plus c'est rapide), avec les tailles que vous mentionnez (7000 et 6M)

import numpy as np
import sklearn
import pandas as pd

Je génère des données, merci Ben j'ai utilisé votre code

#échantillon fictif
np.random.seed(1)
size_nbi = 7000
size_Final = 6000000

nbi = pd.DataFrame({'LONG_017':np.random.random(size=size_nbi)/10+73, 
                    'LAT_016':np.random.random(size=size_nbi)/10+73,})

Final = pd.DataFrame({'X':np.random.random(size=size_Final)/10+73,
                      'Y':np.random.random(size=size_Final)/10+73, 
                      'CUM_OCC':np.random.randint(size_Final,size=size_Final)})

nbi_gps = nbi[["LAT_016", "LONG_017"]].values
final_gps = Final[["Y", "X"]].values

Créer un Balltree

%%time

from sklearn.neighbors import BallTree
import numpy as np

nbi =  np.radians(nbi_gps)
final = np.radians(final_gps)

tree = BallTree(final, leaf_size=12, metric='haversine')

Qui a pris Temps écoulé: 23.8 s

Interroger avec

%%time

rayon = 0.0000003

StoredOCC_indices = tree.query_radius(nbi, r=rayon, return_distance=False, count_only=False)

(moins d'une seconde)

Et obtenir le maximum des éléments qui vous intéressent

StoredOCC = [ np.max( Final.CUM_OCC[i] ) for i in StoredOCC_indices]

Cela crée automatiquement des Nan pour les listes vides, ce qui est pratique. Cela a pris Temps écoulé: 3.64 s

Pour ce rayon le calcul complet sur 6 millions a pris moins d'une minute sur ma machine.