7

2 つの 2D 配列 x(ni, nj) と y(ni,nj) があり、1 つの軸を補間する必要があります。niごとに最後の軸に沿って補間したい。

私が書いた

import numpy as np
from scipy.interpolate import interp1d

z = np.asarray([200,300,400,500,600])
out = []
for i in range(ni):
    f = interp1d(x[i,:], y[i,:], kind='linear')
    out.append(f(z))
out = np.asarray(out)

ただ、この方法は配列サイズが大きすぎるとループするので効率が悪く遅いと思います。このような多次元配列を補間する最速の方法は何ですか? ループなしで線形およびキュービック補間を実行する方法はありますか? ありがとう。

4

2 に答える 2

14

あなたが提案する方法にはpythonループがあるため、値が大きいniと遅くなります。とは言っても、大きくなる予定がない限り、niあまり心配する必要はありません。

次のコードでサンプル入力データを作成しました。

def sample_data(n_i, n_j, z_shape) :
    x = np.random.rand(n_i, n_j) * 1000
    x.sort()
    x[:,0] = 0
    x[:, -1] = 1000
    y = np.random.rand(n_i, n_j)
    z = np.random.rand(*z_shape) * 1000
    return x, y, z

そして、この2つのバージョンの線形補間でそれらをテストしました:

def interp_1(x, y, z) :
    rows, cols = x.shape
    out = np.empty((rows,) + z.shape, dtype=y.dtype)
    for j in xrange(rows) :
        out[j] =interp1d(x[j], y[j], kind='linear', copy=False)(z)
    return out

def interp_2(x, y, z) :
    rows, cols = x.shape
    row_idx = np.arange(rows).reshape((rows,) + (1,) * z.ndim)
    col_idx = np.argmax(x.reshape(x.shape + (1,) * z.ndim) > z, axis=1) - 1
    ret = y[row_idx, col_idx + 1] - y[row_idx, col_idx]
    ret /= x[row_idx, col_idx + 1] - x[row_idx, col_idx]
    ret *= z - x[row_idx, col_idx]
    ret += y[row_idx, col_idx]
    return ret

interp_1Dave's answer に従って、コードの最適化されたバージョンです。interp_2Python ループを一切回避する線形補間のベクトル化された実装です。このようなコーディングには、numpy でのブロードキャストとインデックス作成について十分に理解する必要があり、最適化されていないものinterp1dもあります。代表的な例は、値を補間するビンを見つけることです。ビンが見つかると、interp1d確実に早期にループから抜け出します。上記の関数は、値をすべてのビンと比較しています。

n_iしたがって、結果は、何が何であるかに大きく依存し、補間する値n_jの配列の長さにも大きく依存します。が小さくてが大きいz場合は からの利点が期待でき、逆の場合はからの利点が期待できます。小さいものは に、長いものは に有利です。n_jn_iinterp_2interp_1zinterp_2interp_1

私は実際にさまざまな と を使用して両方のアプローチの時間をn_i計測n_jしzまし(5,)た(50,)。

ここに画像の説明を入力

ここに画像の説明を入力

したがって、 forzの形状は、いつでもを使用し、その他の(5,)場合は使用する必要があるようです。当然のことながら、しきい値はshapeの場合とは異なり、現在は約です。の場合はコードに固執し、そうでない場合は上記のように変更する必要があると結論付けたくなるかもしれませんが、そのステートメントには多くの推定があります。自分でさらに実験したい場合は、グラフを作成するために使用したコードを次に示します。interp_2n_j < 1000interp_1z(50,)n_j < 100n_j * len(z) > 5000interp_2

n_s = np.logspace(1, 3.3, 25)
int_1 = np.empty((len(n_s),) * 2)
int_2 = np.empty((len(n_s),) * 2)
z_shape = (5,)

for i, n_i in enumerate(n_s) :
    print int(n_i)
    for j, n_j in enumerate(n_s) :
        x, y, z = sample_data(int(n_i), int(n_j), z_shape)
        int_1[i, j] = min(timeit.repeat('interp_1(x, y, z)',
                                        'from __main__ import interp_1, x, y, z',
                                        repeat=10, number=1))
        int_2[i, j] = min(timeit.repeat('interp_2(x, y, z)',
                                        'from __main__ import interp_2, x, y, z',
                                        repeat=10, number=1))

cs = plt.contour(n_s, n_s, np.transpose(int_1-int_2))
plt.clabel(cs, inline=1, fontsize=10)
plt.xlabel('n_i')
plt.ylabel('n_j')
plt.title('timeit(interp_2) - timeit(interp_1), z.shape=' + str(z_shape))
plt.show()
于 2013-01-28T23:05:22.667 に答える
2

最適化の 1 つは、次のように結果配列を 1 回割り当てることです。

import numpy as np
from scipy.interpolate import interp1d

z = np.asarray([200,300,400,500,600])
out = np.zeros( [ni, len(z)], dtype=np.float32 ) 
for i in range(ni):
    f = interp1d(x[i,:], y[i,:], kind='linear')
    out[i,:]=f(z)

これにより、 の呼び出しで発生する、実装で発生するメモリのコピーを節約できますout.append(...)。

于 2013-01-28T18:38:14.393 に答える