Comment faire en sorte que scipy.interpolate donne un résultat extrapolé au-delà de la plage d'entrée ?

Question

Je tente de porter un programme qui utilise un interpolateur fait maison (développé par un collègue mathématicien) pour utiliser les interpolateurs fournis par scipy. Je voudrais utiliser ou envelopper l'interpolateur scipy pour qu'il ait un comportement aussi proche que possible de l'ancien interpolateur.

Une différence clé entre les deux fonctions est que dans notre interpolateur original - si la valeur d'entrée est au-dessus ou en dessous de la plage d'entrée, notre interpolateur original extrapolera le résultat. Si vous essayez cela avec l'interpolateur scipy, une ValueError est lancée. Considérez ce programme comme exemple:

import numpy as np
from scipy import interpolate

x = np.arange(0,10)
y = np.exp(-x/3.0)
f = interpolate.interp1d(x, y)

print f(9)
print f(11) # Cause une ValueError, car c'est plus grand que max(x)

Y a-t-il un moyen sensé de faire en sorte que, au lieu de planter, la dernière ligne fasse simplement une extrapolation linéaire, prolongeant les gradients définis par les deux premiers et les deux derniers points jusqu'à l'infini.

Notez que dans le logiciel réel, je n'utilise pas réellement la fonction exp - elle est là pour illustration seulement!

Answer 1

Il peut être plus rapide d'utiliser l'indexation booléenne avec de gros ensembles de données, car l'algorithme vérifie si chaque point est en dehors de l'intervalle, alors que l'indexation booléenne permet une comparaison plus facile et plus rapide.

Par exemple:

# Modules nécessaires
import numpy as np
from scipy.interpolate import interp1d

# Données d'origine
x = np.arange(0,10)
y = np.exp(-x/3.0)

# Classe d'interpolation
f = interp1d(x, y)

# Plage de sortie (assez grande)
xo = np.arange(0, 10, 0.001)

# Approche avec l'indexation booléenne

# Générer un tableau de sortie vide pour les valeurs "y"
yo = np.empty_like(xo)

# Les valeurs inférieures au minimum de "x" sont extrapolées en même temps
low = xo < f.x[0]
yo[low] =  f.y[0] + (xo[low]-f.x[0])*(f.y[1]-f.y[0])/(f.x[1]-f.x[0])

# Les valeurs supérieures au maximum de "x" sont extrapolées en même temps
high = xo > f.x[-1]
yo[high] = f.y[-1] + (xo[high]-f.x[-1])*(f.y[-1]-f.y[-2])/(f.x[-1]-f.x[-2])

# Les valeurs à l'intérieur de la plage d'interpolation sont interpolées directement
inside = np.logical_and(xo >= f.x[0], xo <= f.x[-1])
yo[inside] = f(xo[inside])

Dans mon cas, avec un ensemble de données de 300000 points, cela signifie une accélération de 25,8 à 0,094 secondes, ce qui est plus de 250 fois plus rapide.

Answer 2

J'ai fait cela en ajoutant un point à mes tableaux initiaux. De cette façon, j'évite de définir des fonctions personnalisées, et l'extrapolation linéaire (dans l'exemple ci-dessous : extrapolation vers la droite) semble correcte.

import numpy as np
from scipy import interp as itp

xnew = np.linspace(0,1,51)
x1=xold[-2]
x2=xold[-1]
y1=yold[-2]
y2=yold[-1]
right_val=y1+(xnew[-1]-x1)*(y2-y1)/(x2-x1)
x=np.append(xold,xnew[-1])
y=np.append(yold,right_val)
f = itp(xnew,x,y)

Answer 3

Je n'ai pas assez de réputation pour commenter, mais au cas où quelqu'un rechercherait un wrapper d'extrapolation pour une interpolation linéaire 2D avec scipy, j'ai adapté la réponse qui a été donnée ici pour l'interpolation 1D.

def extrap2d(interpolator):
xs = interpolator.x
ys = interpolator.y
zs = interpolator.z
zs = np.reshape(zs, (-1, len(xs)))
def pointwise(x, y):
    if x < xs[0] or y < ys[0]:
        x1_index = np.argmin(np.abs(xs - x))
        x2_index = x1_index + 1
        y1_index = np.argmin(np.abs(ys - y))
        y2_index = y1_index + 1
        x1 = xs[x1_index]
        x2 = xs[x2_index]
        y1 = ys[y1_index]
        y2 = ys[y2_index]
        z11 = zs[x1_index, y1_index]
        z12 = zs[x1_index, y2_index]
        z21 = zs[x2_index, y1_index]
        z22 = zs[x2_index, y2_index]

        return (z11 * (x2 - x) * (y2 - y) +
        z21 * (x - x1) * (y2 - y) +
        z12 * (x2 - x) * (y - y1) +
        z22 * (x - x1) * (y - y1)
       ) / ((x2 - x1) * (y2 - y1) + 0.0)

    elif x > xs[-1] or y > ys[-1]:
        x1_index = np.argmin(np.abs(xs - x))
        x2_index = x1_index - 1
        y1_index = np.argmin(np.abs(ys - y))
        y2_index = y1_index - 1
        x1 = xs[x1_index]
        x2 = xs[x2_index]
        y1 = ys[y1_index]
        y2 = ys[y2_index]
        z11 = zs[x1_index, y1_index]
        z12 = zs[x1_index, y2_index]
        z21 = zs[x2_index, y1_index]
        z22 = zs[x2_index, y2_index]#

        return (z11 * (x2 - x) * (y2 - y) +
        z21 * (x - x1) * (y2 - y) +
        z12 * (x2 - x) * (y - y1) +
        z22 * (x - x1) * (y - y1)
        ) / ((x2 - x1) * (y2 - y1) + 0.0)
    else:
        return interpolator(x, y)
def ufunclike(xs, ys):
    if  isinstance(xs, int) or isinstance(ys, int) or isinstance(xs, np.int32) or isinstance(ys, np.int32):
        res_array = pointwise(xs, ys)
    else:
        res_array = np.zeros((len(xs), len(ys)))
        for x_c in range(len(xs)):
            res_array[x_c, :] = np.array([pointwise(xs[x_c], ys[y_c]) for y_c in range(len(ys))]).T

    return res_array
return ufunclike

Je n'ai pas beaucoup commenté et je suis conscient que le code n'est pas super propre. Si quelqu'un voit des erreurs, veuillez me le faire savoir. Dans mon cas d'utilisation actuel, cela fonctionne sans problème :)

Answer 4

Je crains qu'il n'y ait pas de moyen facile de le faire dans Scipy à ma connaissance. Vous pouvez, comme je suis assez sûr que vous le savez, désactiver les erreurs de limites et remplir toutes les valeurs de la fonction au-delà de la plage avec une constante, mais cela ne aide vraiment pas. Voir cette question sur la liste de diffusion pour quelques idées supplémentaires. Peut-être pourriez-vous utiliser une sorte de fonction morceau par morceau, mais cela semble être une douleur majeure.

Answer 5

Le code ci-dessous vous donne le module d'extrapolation simple. k est la valeur à laquelle l'ensemble de données y doit être extrapolé en fonction de l'ensemble de données x. Le module numpy est requis.

 def extrapol(k,x,y):
        xm=np.mean(x);
        ym=np.mean(y);
        sumnr=0;
        sumdr=0;
        length=len(x);
        for i in range(0,length):
            sumnr=sumnr+((x[i]-xm)*(y[i]-ym));
            sumdr=sumdr+((x[i]-xm)*(x[i]-xm));

        m=sumnr/sumdr;
        c=ym-(m*xm);
        return((m*k)+c)