等效于C ++ Cuda中的一维numpy.interpCUDA中为Lerp

问题描述

我有两个向量“ 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，以便可以使用权重向量来解决此问题。

解决方法

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

从我所见，numpy函数还允许x数组大致上是任意的。

为了进行适当的插值，我们必须找到每个xp值所属的x的“段”。我唯一想到的方法是执行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,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
$

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