在不使用 insert 的情况下插入 Numpy 数组

Interpolate a Numpy array without using insert

所以我正在尝试做一些插值,我只是想知道我正在使用的代码是否有可能更高效。

所以问题是这样的。我有一组值,其中包含来自许多不同站点的数据。它的尺寸为 N x 85,N 可能因世界地区而异。

目前他们之间的时间间隔并不均匀。列 0:51 之间有 3 小时的时间差,然后列 52:84 之间有 6 小时的时间差。我只想在它们之间做一个简单的线性插值,这样所有的行之间都有 3 小时的间隔,所以基本上只是计算行之间的均值 52:84 并将它们插入正确的地方。

这是一个小规模示例的快速代码示例,这是我目前一直在做的但效率不高(插入功能是我想避免的)。

import numpy as np
np.set_printoptions(linewidth=np.nan)

array = np.random.rand(5, 10) * 10
print(array)

interpolation_array = np.empty((5, 4))

for i, j in enumerate(list(range(5, 9))):
    interpolation_array[:, i] = np.mean(array[:, j:(j+2)], axis=1)
print(interpolation_array)

# This is the line that is not memory efficient
final_array = np.insert(array, list(range(6, 10)), interpolation_array, axis=1)

print(final_array)

因此,正如@hpaulj 所建议的,我使用 Numba(LLVM 编译器比 Numpy 更快)实现了一个实现,它创建了一个大矩阵来写入所有值,然后使用并行循环来填充它并执行简单的线性插值。这是我使用的代码。

import numpy as np
from numba import jit, prange

@jit(nopython=True)
def numba_approach(array, c_start):  # c_start is column to start interpolating in (zero indexed)
    num_rows = array.shape[0]
    num_cols = array.shape[1]
    num_interp_columns = num_cols - 1 - c_start

    final_array = np.empty((num_rows, num_cols + num_interp_columns))

    # Populate the Portion That is not interpolated
    for i in prange(num_rows):
        for j in range(c_start + 1):
            final_array[i, j] = array[i, j]

    z = 1
    for j in prange(c_start + 2, num_cols + num_interp_columns, 2):
        for i in range(num_rows):
            final_array[i, j] = array[i, j - z]
        z += 1

    # Interpolate
    for j in prange(c_start + 1, num_cols + num_interp_columns - 1, 2):
        for i in range(num_rows):
            final_array[i, j] = (final_array[i, j - 1] + final_array[i, j + 1]) / 2

    return final_array

这给了我大约 4 倍的加速,这是相对重要的,因为这段代码经常 运行。

基准:

%timeit numpy_approach(test_array)

100 loops, best of 3: 2.16 ms per loop

%timeit numba_approach(test_array, 47)

1000 loops, best of 3: 446 µs per loop

我认为正确的方法是使用适当的线性插值,例如 numpy.interp。必须明确定义每个值的对应时间。这允许泛化到值之间的任何时间间隔。

import numpy as np
import matplotlib.pylab as plt

times = [0, 3, 6, 9, 12, 18, 24, 30, 36]
values = np.random.rand(len(times))

times_regular = np.arange(0, times[-1]+1, 3)

values_regular = np.interp(times_regular, times, values)

plt.plot(times_regular, values_regular, 's', label='evently spaced');
plt.plot(times, values, 'd-', label='measure');
plt.xlabel('time'); plt.ylabel('value'); plt.legend();

基于数组的解决方案,因为只需要中点,可以是:

data_every6hr = np.array([[7, 5, 6, 9, 8, 5],
                          [7, 9, 6, 5, 6, 9],
                          [5, 6, 7, 9, 8, 8],
                          [5, 9, 8, 5, 7, 6]], dtype=float)

# Perform the interpolation for every line
intermediate_values = (data_every6hr[:, 1:] + data_every6hr[:, :-1])/2

# Insert the interpolated values before each column:
data_every3hr = np.insert(data_every6hr,
                          range(1, data_every6hr.shape[1]),
                          intermediate_values,
                          axis=1)

print(data_every3hr)
#array([[7. , 6. , 5. , 5.5, 6. , 7.5, 9. , 8.5, 8. , 6.5, 5. ],
#       [7. , 8. , 9. , 7.5, 6. , 5.5, 5. , 5.5, 6. , 7.5, 9. ],
#       [5. , 5.5, 6. , 6.5, 7. , 8. , 9. , 8.5, 8. , 8. , 8. ],
#       [5. , 7. , 9. , 8.5, 8. , 6.5, 5. , 6. , 7. , 6.5, 6. ]])

data_every6hr 只是输入数组中数据间隔 6 小时的部分。