スムーズ機能とアップサンプリング

最も簡単なタスクから始めましょう。スムーズ テスト関数で、シェイプのメッシュからシェイプ[6, 7]の 1 つへのアップサンプリングがどのように機能するかを次に示します。[20, 21]

これは単純なタスクですが、出力にはすでに微妙な違いがあります。一見すると、3 つの出力はすべて妥当です。基礎関数に関する事前の知識に基づくと、注目すべき特徴が 2 つあります。 の中間のケースでは、griddataデータが最も歪んでいます。y == -1プロットの境界 (xラベルに最も近い) に注意してください。関数は厳密にゼロになるはずです ( はy == -1滑らかな関数のノード線であるため) が、 の場合はそうではありませんgriddata。また、x == -1プロットの境界 (後ろ、左側) にも注意してください。基礎関数は で局所的最大値を持ちます (境界付近で勾配がゼロであることを意味します) が[-1, -0.5]、griddata出力はこの領域で明らかにゼロ以外の勾配を示しています。影響は微妙ですが、それでもバイアスであることに変わりはありません。

悪の機能とアップサンプリング

少し難しいタスクは、evil 関数でアップサンプリングを実行することです。

3 つの方法の間には明らかな違いが現れ始めています。表面プロットを見ると、出力に明らかな偽の極値が現れていますinterp2d(プロットされた表面の右側にある 2 つの山に注目してください)。一見するとgriddata、 と は同様の結果を生成するように見えますが、その付近に局所的最小値を生成する機能は、基礎となる関数には存在しません。RBFInterpolator[0.4, -0.4]

しかし、 にははるかに優れている重要な側面が 1 つありますRBFInterpolator。それは、基礎となる関数の対称性を尊重していることです (もちろん、サンプル メッシュの対称性によっても可能になります)。 からの出力はgriddata、滑らかなケースですでに弱く見えるサンプル ポイントの対称性を破壊します。

滑らかな関数と散在したデータ

ほとんどの場合、散在したデータに対して補間を実行する必要があります。このため、これらのテストがより重要になると考えています。上記のように、サンプル ポイントは対象ドメイン内で擬似的に均一に選択されました。現実的なシナリオでは、測定ごとに追加のノイズが発生する可能性があるため、生のデータを補間することが理にかなっているかどうかを検討する必要があります。

スムーズ関数の出力:

すでにちょっとしたホラーショーが起こっています。私は、最低限の情報量を残すために、プロット専用interp2dに から までの出力を切り取りました。明らかに、基礎となる形状の一部は存在しますが、この方法が完全に機能しない大きなノイズ領域があります。 の2番目のケースは形状をかなりうまく再現しますが、等高線プロットの境界に白い領域があることに注意してください。これは、 が入力データポイントの凸包内でのみ機能するという事実によるものです（言い換えれば、[-1, 1]griddatagriddata外挿）。凸包の外側にある出力ポイントについては、デフォルトの NaN 値を維持しました。2^これらの機能を考慮すると、RBFInterpolator最もパフォーマンスが良いようです。

悪意のある機能と散在するデータ

そして私たち全員が待ち望んでいた瞬間がやってきました。

が諦めるのは、それほど驚くことではありませんinterp2d。実際、 の呼び出し中に、スプラインを構築できないと文句を言うinterp2dフレンドリーな がいくつかあることを覚悟しておく必要があります。RuntimeWarning他の 2 つの方法に関しては、RBFInterpolator結果が外挿されるドメインの境界付近でも、 が最良の出力を生成するようです。

それでは、優先順位の高い順に 3 つの方法について少し説明しましょう (最悪の方法は誰にも読まれにくいものになります)。

scipy.interpolate.RBFInterpolator

クラス名の RBF はRBFInterpolator「ラジアル基底関数」の略です。正直に言うと、この投稿の調査を始めるまでこのアプローチを検討したことはありませんでしたが、将来的にはこれを使用することになるでしょう。

スプラインベースの方法 (後述) と同様に、使用方法は 2 つのステップに分かれています。まず、RBFInterpolator入力データに基づいて呼び出し可能なクラス インスタンスを作成し、次に指定された出力メッシュに対してこのオブジェクトを呼び出して補間結果を取得します。スムーズ アップサンプリング テストの例:

import scipy.interpolate as interp

sparse_points = np.stack([x_sparse.ravel(), y_sparse.ravel()], -1)  # shape (N, 2) in 2d
dense_points = np.stack([x_dense.ravel(), y_dense.ravel()], -1)  # shape (N, 2) in 2d

zfun_smooth_rbf = interp.RBFInterpolator(sparse_points, z_sparse_smooth.ravel(),
                                         smoothing=0, kernel='cubic')  # explicit default smoothing=0 for interpolation
z_dense_smooth_rbf = zfun_smooth_rbf(dense_points).reshape(x_dense.shape)  # not really a function, but a callable class instance

zfun_evil_rbf = interp.RBFInterpolator(sparse_points, z_sparse_evil.ravel(),
                                       smoothing=0, kernel='cubic')  # explicit default smoothing=0 for interpolation
z_dense_evil_rbf = zfun_evil_rbf(dense_points).reshape(x_dense.shape)  # not really a function, but a callable class instance

sparse_points = np.stack([x_sparse.ravel(), y_sparse.ravel()], -1)  # shape (N, 2) in 2d
z_dense_smooth_griddata = interp.griddata(sparse_points, z_sparse_smooth.ravel(),
                                          (x_dense, y_dense), method='cubic')  # default method is linear

z_dense_smooth_griddata = interp.griddata((x_sparse.ravel(), y_sparse.ravel()),
                                          z_sparse_smooth.ravel(), (x_dense, y_dense), method='cubic')

# reminder: x_sparse and y_sparse are of shape [6, 7] from numpy.meshgrid
zfun_smooth_interp2d = interp.interp2d(x_sparse, y_sparse, z_sparse_smooth, kind='cubic')   # default kind is 'linear'
# reminder: x_dense and y_dense are of shape (20, 21) from numpy.meshgrid
xvec = x_dense[0,:] # 1d array of unique x values, 20 elements
yvec = y_dense[:,0] # 1d array of unique y values, 21 elements
z_dense_smooth_interp2d = zfun_smooth_interp2d(xvec, yvec)   # output is (20, 21)-shaped array

kind = 'cubic'
if kind == 'linear':
    kx = ky = 1
elif kind == 'cubic':
    kx = ky = 3
elif kind == 'quintic':
    kx = ky = 5
# bisplrep constructs a spline representation, bisplev evaluates the spline at given points
bisp_smooth = interp.bisplrep(x_sparse.ravel(), y_sparse.ravel(),
                              z_sparse_smooth.ravel(), kx=kx, ky=ky, s=0)
z_dense_smooth_bisplrep = interp.bisplev(xvec, yvec, bisp_smooth).T  # note the transpose

if not rectangular_grid:
    # TODO: surfit is really not meant for interpolation!
    self.tck = fitpack.bisplrep(x, y, z, kx=kx, ky=ky, s=0.0)

tx, ty, c, o = _fitpack._surfit(x, y, z, w, xb, xe, yb, ye, kx, ky,
                                task, s, eps, tx, ty, nxest, nyest,
                                wrk, lwrk1, lwrk2)