等效于 C++ Cuda 中的一维 numpy.interp（CUDA 中的 Lerp）

Question

我有两个向量 'xp' 和 'fp' 分别对应于数据的 x 和 y 值。第三个向量“x”是我想要评估插值的 x 坐标。我使用 NumPy 的 interp 函数在 python 中的结果符合预期。

import numpy as np 

xp = np.array([1.0, 1.5, 2.0, 2.5, 3.5, 4.0, 4.5, 7.0, 8.0, 9.0, 10.0, 14.0, 17.0, 18.0, 20.0])
yp = xp**2
x = np.array([3,5])
np.interp(x,xp,yp) #array([  9.25,  26.  ])

我的问题是如何在 cuda 内核中复制这个算法？

这是我的尝试：

大小 --> len(xp) == len(fp), x_size --> len(x)。

__global__ void lerp_kernel(double* xp, double* yp, double* output_result, 
                          double* x, unsigned int size, unsigned int x_size)
{
  for( unsigned int idx = blockIdx.x*blockDim.x + threadIdx.x ; idx < size ; idx += 
  blockDim.x*gridDim.x )

{     
 
 if(idx > size - 2){    
  idx = size - 2; 
  }
  double dx = xp[idx + 1] - xp[idx];
  double dy = yp[idx + 1] - yp[idx];
  double dydx = dy/dx ;
  double lerp_result = yp[idx] + (dydx * (x[idx] - xp[idx]));
  output_result[idx] = lerp_result;
  
   }
  }

我认为我犯的一个错误是我没有在 xp 中搜索包含 x 的索引范围（在 python 中使用类似 numpy.searchsorted 的东西）。我不确定如何在 CUDA 中实现这部分。如果有人知道在 cuda 中做 lerp 的更好方法，请告诉我。

我在 Cg (https://developer.download.nvidia.com/cg/lerp.html) 的文档中看到了 lerp 函数，但是这些函数需要 x 向量的权重向量（0-1 之间的分数）。我不知道如何重新调整 x 以便我可以使用权重向量来解决这个问题。

Answer 1

要模仿numpy.interp 的行为需要几个步骤。我们将至少做出一个简化假设：numpy.interp 函数expects 你的xp 数组正在增加（我们也可以说“排序”）。否则，它特别提到需要进行（内部）排序。我们将跳过这种情况并假设您的 xp 数组正在增加，如您在此处所示。

据我所知，numpy 函数还允许x 数组或多或少是任意的。

为了进行正确的插值，我们必须找到每个x 值所属的xp 的“段”。我能想到的唯一方法是执行binary search。 （还要注意thrust 有方便的二进制搜索）

那么这个过程将是：

1. 对x 中的每个元素使用二分搜索，在xp 中找到相应的“段”（即端点）
2. 使用基于已识别线段的直线方程 (y=mx+b) 来计算端点之间的插值

这是一个例子：

$ cat t40.cu
#include <iostream>

typedef int lt;

template <typename T>
__device__ void bsearch_range(const T *a, const T key, const lt len_a, lt *idx){
  lt lower = 0;
  lt upper = len_a;
  lt midpt;
  while (lower < upper){
    midpt = (lower + upper)>>1;
    if (a[midpt] < key) lower = midpt +1;
    else upper = midpt;
    }
  *idx = lower;
  return;
  }


template <typename T>
__global__ void my_interp(const T *xp, const T *yp, const lt xp_len, const lt x_len, const T *x, T *y){

  for (lt i = threadIdx.x+blockDim.x*blockIdx.x; i < x_len; i+=gridDim.x*blockDim.x){
    T val = x[i];
    if ((val >= xp[0]) && (val < xp[xp_len - 1])){
      lt idx;
      bsearch_range(xp, val, xp_len, &idx);
      T xlv = xp[idx-1];
      T xrv = xp[idx];
      T ylv = yp[idx-1];
      T yrv = yp[idx];
     // y  =      m                *   x       + b
      y[i] = ((yrv-ylv)/(xrv-xlv)) * (val-xlv) + ylv;
    }
  }
}

typedef float mt;
const int nTPB = 256;
int main(){

  mt xp[] = {1.0, 1.5, 2.0, 2.5, 3.5, 4.0, 4.5, 7.0, 8.0, 9.0, 10.0, 14.0, 17.0, 18.0, 20.0};
  lt xp_len = sizeof(xp)/sizeof(xp[0]);
  mt *yp = new mt[xp_len];
  for (lt i = 0; i < xp_len; i++) yp[i] = xp[i]*xp[i];
  mt x[] = {3,5};
  lt x_len = sizeof(x)/sizeof(x[0]);
  mt *y = new mt[x_len];
  mt *d_xp, *d_x, *d_yp, *d_y;
  cudaMalloc(&d_xp, xp_len*sizeof(xp[0]));
  cudaMalloc(&d_yp, xp_len*sizeof(yp[0]));
  cudaMalloc(&d_x,   x_len*sizeof( x[0]));
  cudaMalloc(&d_y,   x_len*sizeof( y[0]));
  cudaMemcpy(d_xp, xp, xp_len*sizeof(xp[0]), cudaMemcpyHostToDevice);
  cudaMemcpy(d_yp, yp, xp_len*sizeof(yp[0]), cudaMemcpyHostToDevice);
  cudaMemcpy(d_x, x, x_len*sizeof(x[0]), cudaMemcpyHostToDevice);
  my_interp<<<(x_len+nTPB-1)/nTPB, nTPB>>>(d_xp, d_yp, xp_len, x_len, d_x, d_y);
  cudaMemcpy(y, d_y, x_len*sizeof(y[0]), cudaMemcpyDeviceToHost);
  for (lt i = 0; i < x_len; i++) std::cout << y[i] << ",";
  std::cout << std::endl;
}
$ nvcc -o t40 t40.cu
$ cuda-memcheck ./t40
========= CUDA-MEMCHECK
9.25,26,
========= ERROR SUMMARY: 0 errors
$

我并不是说上述代码没有缺陷或适用于任何特定目的。我的目标是演示一种方法，而不是提供经过全面测试的代码。特别是，可能需要仔细测试边缘情况（例如，在整个范围的边缘或范围之外的值）。