Interpolation linéaire rapide en Num Py /Scipy " le long d'un chemin"

Une interpolation linéaire entre deux valeurs y1,y2 à des endroits x1 et x2, en ce qui concerne le point xi est tout simplement:

yi = y1 + (y2-y1) * (xi-x1) / (x2-x1)

avec quelques expressions num Py vectorisées nous pouvons sélectionner les points pertinents de l'ensemble de données et appliquer la fonction ci-dessus:

I = np.searchsorted(altitudes, location)

x1 = altitudes[I-1]
x2 = altitudes[I]

time = np.arange(len(alltemps))
y1 = alltemps[time,I-1]
y2 = alltemps[time,I]

xI = location

yI = y1 + (y2-y1) * (xI-x1) / (x2-x1)

le problème est que certains points se situent sur les limites de (ou même en dehors de) la portée connue, qui devrait être prise en compte:

I = np.searchsorted(altitudes, location)
same = (location == altitudes.take(I, mode='clip'))
out_of_range = ~same & ((I == 0) | (I == altitudes.size))
I[out_of_range] = 1  # Prevent index-errors

x1 = altitudes[I-1]
x2 = altitudes[I]

time = np.arange(len(alltemps))
y1 = alltemps[time,I-1]
y2 = alltemps[time,I]

xI = location

yI = y1 + (y2-y1) * (xI-x1) / (x2-x1)
yI[out_of_range] = np.nan

Heureusement, Scipy fournit déjà l'interpolation ND, qui s'occupe aussi facilement des temps de décalage, par exemple:

from scipy.interpolate import interpn

time = np.arange(len(alltemps))

M = 150
hiketime = np.linspace(time[0], time[-1], M)
location = np.linspace(altitudes[0], altitudes[-1], M)
xI = np.column_stack((hiketime, location))

yI = interpn((time, altitudes), alltemps, xI)

Voici un code de référence (sans pandas en fait, peu, j'ai fait notamment de la solution de l'autre):

import numpy as np
from scipy.interpolate import interp1d, interpn

def original():
    return np.array([interp1d(altitudes, alltemps[i, :])(loc)
                                for i, loc in enumerate(location)])

def OP_self_answer():
    return np.diagonal(interp1d(altitudes, alltemps)(location))

def interp_checked():
    I = np.searchsorted(altitudes, location)
    same = (location == altitudes.take(I, mode='clip'))
    out_of_range = ~same & ((I == 0) | (I == altitudes.size))
    I[out_of_range] = 1  # Prevent index-errors

    x1 = altitudes[I-1]
    x2 = altitudes[I]

    time = np.arange(len(alltemps))
    y1 = alltemps[time,I-1]
    y2 = alltemps[time,I]

    xI = location

    yI = y1 + (y2-y1) * (xI-x1) / (x2-x1)
    yI[out_of_range] = np.nan

    return yI

def scipy_interpn():
    time = np.arange(len(alltemps))
    xI = np.column_stack((time, location))
    yI = interpn((time, altitudes), alltemps, xI)
    return yI

N, sigma = 1000., 5

basetemps = 70 + (np.random.randn(N) * sigma)
midtemps = 50 + (np.random.randn(N) * sigma)
toptemps = 40 + (np.random.randn(N) * sigma)
trend = np.sin(4 / N * np.arange(N)) * 30
trend = trend[:, np.newaxis]
alltemps = np.array([basetemps, midtemps, toptemps]).T + trend
altitudes = np.array([500, 1500, 4000], dtype=float)
location = np.linspace(altitudes[0], altitudes[-1], N)

funcs = [original, interp_checked, scipy_interpn]
for func in funcs:
    print(func.func_name)
    %timeit func()

from itertools import combinations
outs = [func() for func in funcs]
print('Output allclose:')
print([np.allclose(out1, out2) for out1, out2 in combinations(outs, 2)])

Avec le résultat suivant sur mon système:

original
10 loops, best of 3: 184 ms per loop
OP_self_answer
10 loops, best of 3: 89.3 ms per loop
interp_checked
1000 loops, best of 3: 224 µs per loop
scipy_interpn
1000 loops, best of 3: 1.36 ms per loop
Output allclose:
[True, True, True, True, True, True]

Scipy est interpn souffre un peu en termes de vitesse par rapport à la méthode la plus rapide, mais pour sa généralité et sa facilité de il est certainement le chemin à parcourir.

je vais proposer un peu de progrès. Dans le second cas (interpolation "le long d'un chemin"), nous faisons beaucoup de différentes fonctions d'interpolation. Une chose que nous pourrions essayer est de faire juste une fonction d'interpolation (une qui fait interpolation dans la dimension d'altitude sur tous les temps comme dans le premier cas ci-dessus) et d'évaluer cette fonction encore et encore (d'une manière vectorisée). Cela nous donnerait beaucoup plus de données que nous ne le voulons (cela nous donnerait une matrice de 1 000 x 1 000 au lieu d'une matrice de 1 000 éléments). vecteur.) Mais alors notre résultat cible serait juste le long de la diagonale. Donc, la question Est, est-ce que l'appel d'une seule fonction sur un chemin d'arguments plus complexes s'exécute plus rapidement que de faire de nombreuses fonctions et de les appeler avec des arguments simples?

La réponse est oui!

la clé est que la fonction d'interpolation retourne par scipy.interpolate.interp1d est capable d'accepter un numpy.ndarray comme entrée. Ainsi, vous pouvez effectivement appeler la fonction d'interpolation plusieurs fois à la vitesse C en vous alimentant dans une entrée vectorielle. C'est-à-dire: c'est beaucoup plus rapide que d'écrire une boucle for qui appelle la fonction interpolation encore et encore sur une entrée scalaire. Alors que nous calculons beaucoup de points de données que nous finissons par jeter, nous gagnons encore plus de temps en ne construisant pas beaucoup de différentes fonctions d'interpolation que nous utilisons à peine.

old_way = interped_along_path = np.array([interp1d(altitudes, finaltemps.values[i, :])(loc) 
                                                      for i, loc in enumerate(location)])
# look ma, no for loops!
new_way = np.diagonal(interp1d(altitudes, finaltemps.values)(location)) 
# note, `location` is a vector!
abs(old_way - new_way).max()
#-> 0.0

et encore:

%%timeit
res = np.diagonal(interp1d(altitudes, finaltemps.values)(location))
#-> 100 loops, best of 3: 16.7 ms per loop

donc cette approche nous donne un facteur de 10 de mieux! Quelqu'un peut-il faire mieux? Ou suggérer une approche totalement différente?