Schnelle lineare Interpolation in Numpy / Scipy "entlang eines Weges"

Angenommen, ich habe Daten von Wetterstationen bei 3 (bekannten) Höhen auf einem Berg. Speziell erfasst jede Station jede Temperatur eine Temperaturmessung an ihrer Stelle. Ich habe zwei Arten von Interpolation, die ich gerne durchführen würde. Und ich möchte in der Lage sein, jedes schnell durchzuführen.

Also lasst uns einige Daten einrichten:

import numpy as np from scipy.interpolate import interp1d import pandas as pd import seaborn as sns np.random.seed(0) 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) alltemps = np.array([basetemps, midtemps, toptemps]).T # note transpose! trend = np.sin(4 / N * np.arange(N)) * 30 trend = trend[:, np.newaxis] altitudes = np.array([500, 1500, 4000]).astype(float) finaltemps = pd.DataFrame(alltemps + trend, columns=altitudes) finaltemps.index.names, finaltemps.columns.names = ['Time'], ['Altitude'] finaltemps.plot() 

Großartig, so sehen unsere Temperaturen so aus: Rohtemperatur-Daten

Interpoliere alle Zeiten für die gleiche Höhe:

Ich denke, das ist ziemlich einfach. Sagen Sie, ich möchte die Temperatur auf einer Höhe von 1.000 für jedes Mal bekommen. Ich kann einfach gebaut in scipy Interpolation Methoden:

 interping_function = interp1d(altitudes, finaltemps.values) interped_to_1000 = interping_function(1000) fig, ax = plt.subplots(1, 1, figsize=(8, 5)) finaltemps.plot(ax=ax, alpha=0.15) ax.plot(interped_to_1000, label='Interped') ax.legend(loc='best', title=finaltemps.columns.name) 

Temperatur mit statischen interp

Das funktioniert gut. Und lasst uns über die Geschwindigkeit sehen:

 %%timeit res = interp1d(altitudes, finaltemps.values)(1000) #-> 1000 loops, best of 3: 207 µs per loop 

Interpolieren Sie "entlang eines Weges":

So jetzt habe ich ein zweites, verwandtes Problem. Sagen Sie, ich kenne die Höhe einer Wanderpartei als Funktion der Zeit, und ich möchte die Temperatur an ihrem (bewegten) Ort durch lineare Interpolation meiner Daten durch die Zeit berechnen. Besonders die Zeiten, in denen ich die Lage der Wanderpartei kenne, sind die gleichen Zeiten, an denen ich die Temperaturen an meinen Wetterstationen kenne. Ich kann das ohne zu viel Aufwand machen:

 location = np.linspace(altitudes[0], altitudes[-1], N) interped_along_path = np.array([interp1d(altitudes, finaltemps.values[i, :])(loc) for i, loc in enumerate(location)]) fig, ax = plt.subplots(1, 1, figsize=(8, 5)) finaltemps.plot(ax=ax, alpha=0.15) ax.plot(interped_along_path, label='Interped') ax.legend(loc='best', title=finaltemps.columns.name) 

Temperatur mit beweglichem Interp

So funktioniert das wirklich nett, aber es ist wichtig zu beachten, dass die obige Schlüsselliste das Listenverständnis benutzt, um eine enorme Menge an Arbeit zu verbergen. Im vorigen Fall erstellt scipy eine einzige Interpolationsfunktion für uns und bewertet sie einmal auf einer großen Menge an Daten. In diesem Fall konstruiert scipy tatsächlich N individuelle Interpolationsfunktionen und bewertet jedes einmal auf einer kleinen Datenmenge. Das fühlt sich inhärent ineffizient an. Es gibt eine für Loop lauern hier (in der Liste Verständnis) und darüber hinaus fühlt sich das nur schlaff.

Es ist nicht überraschend, das ist viel langsamer als der vorherige Fall:

 %%timeit res = np.array([interp1d(altitudes, finaltemps.values[i, :])(loc) for i, loc in enumerate(location)]) #-> 10 loops, best of 3: 145 ms per loop 

Also das zweite Beispiel läuft 1.000 langsamer als das erste. Im Einklang mit der Vorstellung, dass das schwere Heben ist die "machen eine lineare Interpolation Funktion" Schritt … was geschieht 1.000 Mal im zweiten Beispiel, sondern nur einmal in der ersten.

Also, die Frage: Gibt es einen besseren Weg, um das zweite Problem zu lösen? Zum Beispiel gibt es einen guten Weg, um es mit einer 2-dimensinoalen Interpolation einzurichten (was vielleicht den Fall behandeln könnte, wo die Zeiten, an denen die Wanderparties bekannt sind, nicht die Zeiten sind, an denen die Temperaturen abgetastet wurden)? Oder gibt es einen besonders glatten Weg, um Dinge hier zu behandeln, wo die Zeiten sich ausrichten? Oder andere?

3 Solutions collect form web for “Schnelle lineare Interpolation in Numpy / Scipy "entlang eines Weges"”

Für einen festen Zeitpunkt können Sie die folgende Interpolationsfunktion nutzen:

 g(a) = cc[0]*abs(a-aa[0]) + cc[1]*abs(a-aa[1]) + cc[2]*abs(a-aa[2]) 

Wo a die Höhen des Wanderers ist, aa der Vektor mit den 3 Messhöhen und cc ist ein Vektor mit den Koeffizienten. Es gibt drei Dinge zu beachten:

  1. Für gegebene Temperaturen ( alltemps ), die aa entsprechen, kann die Bestimmung von cc durch Lösen einer linearen Matrixgleichung unter Verwendung von np.linalg.solve() .
  2. g(a) ist einfach zu vektorisieren für eine (N,) dimensionale a und (N, 3) dimensionale cc (einschließlich np.linalg.solve() ).
  3. g(a) heißt ein erster Ordnung univariate Spline-Kernel (für drei Punkte). Mit abs(a-aa[i])**(2*d-1) würde die Spline-Reihenfolge zu d ändern. Dieser Ansatz könnte eine vereinfachte Version eines Gaußschen Prozesses im Maschinellen Lernen interpretieren.

Also der Code wäre:

 import matplotlib.pyplot as plt import numpy as np import seaborn as sns # generate temperatures np.random.seed(0) N, sigma = 1000, 5 trend = np.sin(4 / N * np.arange(N)) * 30 alltemps = np.array([tmp0 + trend + sigma*np.random.randn(N) for tmp0 in [70, 50, 40]]) # generate attitudes: altitudes = np.array([500, 1500, 4000]).astype(float) location = np.linspace(altitudes[0], altitudes[-1], N) def doit(): """ do the interpolation, improved version for speed """ AA = np.vstack([np.abs(altitudes-a_i) for a_i in altitudes]) # This is slighty faster than np.linalg.solve(), because AA is small: cc = np.dot(np.linalg.inv(AA), alltemps) return (cc[0]*np.abs(location-altitudes[0]) + cc[1]*np.abs(location-altitudes[1]) + cc[2]*np.abs(location-altitudes[2])) t_loc = doit() # call interpolator # do the plotting: fg, ax = plt.subplots(num=1) for alt, t in zip(altitudes, alltemps): ax.plot(t, label="%d feet" % alt, alpha=.5) ax.plot(t_loc, label="Interpolation") ax.legend(loc="best", title="Altitude:") ax.set_xlabel("Time") ax.set_ylabel("Temperature") fg.canvas.draw() 

Messung der Zeit gibt:

 In [2]: %timeit doit() 10000 loops, best of 3: 107 µs per loop 

Update: Ich habe das ursprüngliche Listenverständnis in doit() , um die Geschwindigkeit um 30% zu importieren (Für N=1000 ).

Weiterhin, wie zum Vergleich angefordert, @ moarningsuns Benchmark-Code-Block auf meiner Maschine:

 10 loops, best of 3: 110 ms per loop interp_checked 10000 loops, best of 3: 83.9 µs per loop scipy_interpn 1000 loops, best of 3: 678 µs per loop Output allclose: [True, True, True] 

Beachten Sie, dass N=1000 eine relativ kleine Zahl ist. Mit N=100000 werden die Ergebnisse erzeugt:

 interp_checked 100 loops, best of 3: 8.37 ms per loop %timeit doit() 100 loops, best of 3: 5.31 ms per loop 

Dies zeigt, dass dieser Ansatz besser für große N als der interp_checked Ansatz interp_checked .

Eine lineare Interpolation zwischen zwei Werten y1 , y2 an den Orten x1 und x2 , bezogen auf Punkt xi ist einfach:

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

Mit einigen vektorisierten Numpy-Ausdrücken können wir die relevanten Punkte aus dem Datensatz auswählen und die obige Funktion anwenden:

 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) 

Das Problem ist, dass einige Punkte auf den Grenzen von (oder sogar außerhalb) des bekannten Bereichs liegen, die berücksichtigt werden sollten:

 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 

Zum Glück bietet Scipy bereits eine ND-Interpolation, die auch so einfach die Mismatching-Zeiten kümmert, zum Beispiel:

 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) 

Hier ist ein Benchmark-Code (ohne pandas eigentlich, Bit ich habe die Lösung von der anderen Antwort):

 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)]) 

Mit folgendem Ergebnis auf meinem System:

 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's interpn leidet etwas im Hinblick auf die Geschwindigkeit im Vergleich zu der sehr schnellsten Methode, aber für seine Allgemeingültigkeit und Benutzerfreundlichkeit ist es definitiv der Weg zu gehen.

Ich werde ein bisschen Fortschritt anbieten. Im zweiten Fall (Interpolation "entlang eines Pfades") machen wir viele verschiedene Interpolationsfunktionen. Eine Sache, die wir ausprobieren können, ist, nur eine Interpolationsfunktion zu machen (eine, die Interpolation in der Höhendimension über alle Zeiten wie im ersten Fall oben) und diese Funktion über und über (auf vektorisierte Weise) auswertet. Das würde uns mehr Daten geben, als wir wollen (es würde uns eine 1.000 x 1.000 Matrix anstelle eines 1.000-Element-Vektors geben). Aber dann wäre unser Ziel Ergebnis nur entlang der Diagonale. Also die Frage ist, ruft eine einzelne Funktion auf dem Weg komplexere Argumente laufen schneller als viele Funktionen und rufen sie mit einfachen Argumenten?

Die Antwort ist ja!

Der Schlüssel ist, dass die interpolierende Funktion, die von scipy.interpolate.interp1d wird, in der Lage ist, ein numpy.ndarray als seine Eingabe zu akzeptieren. So können Sie die Interpolationsfunktion oftmals mit C-Geschwindigkeit effektiv anrufen, indem Sie einen Vektor-Eingang eingeben. Dh das ist so, viel schneller als das Schreiben einer for-Schleife, die die Interpolationsfunktion immer wieder auf einem Skalar-Eingang aufruft. Wenn wir also viele viele Datenpunkte berechnen, die wir abwischen, sparen wir noch mehr Zeit, indem wir nicht viele verschiedene interpolierende Funktionen konstruieren, die wir kaum nutzen.

 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 

und doch:

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

Also dieser Ansatz bekommt uns einen Faktor von 10 besser! Kann jemand besser machen? Oder schlagen Sie einen ganz anderen Ansatz vor?

  • Probleme mit scipy.optimize mit Matrix als Eingabe, Grenzen, Einschränkungen
  • Python theano Optimierungsfehler durch: local_dot_to_dot22
  • Schneiden von spärlichen Matrizen in Scipy - welche Typen funktionieren am besten?
  • Was ist ein "Skalar" in numpy?
  • Griddata Laufzeitfehler - Python / SciPy (Interpolation)
  • Matplotlib - Schneller Weg, um viele Subplots zu erstellen?
  • Ermittlung eines spärlichen Matrixquotienten
  • Speichern / laden scipy sparse csr_matrix im tragbaren Datenformat
  • Wie bekomme ich Wert aus einer theano tensor Variable, die von einer freigegebenen Variablen unterstützt wird?
  • Polyfit () bekam ein unerwartetes Keyword-Argument 'w'
  • Fitting (ein Gaußer) mit Scipy vs. ROOT et al
  • Python ist die beste Programmiersprache der Welt.