CUDA 从非均匀采样中重新采样



我想从非制服采样中重新采样(插值(序列。我不认为tex有效,因为它基本上是假设你的样本是均匀的插值?做搜索会太慢吗?

我应该做推力吗?任何指针都值得赞赏。任何例子都会有很大帮助。

更新:

假设带有圆圈标记的线是我的样本。我知道每个圆点的值。显然,样品均匀地分布在水平轴上。

现在,我想知道采样线下方行上每个 x 标记的值。x 标记沿线均匀分布。

---o--------o----o------o------o------o------(采样(

-

-X-----X-----X-----X-----X-----X-----X---(已知插值(

所以我想知道如何使用 CUDA 获取每个 x 标记位置的值?显然,使用 C/C++ 的最基本算法是针对每个 x 标记位置,搜索两个最近的圆位置,然后进行线性插值。但在这种情况下,您需要首先对两个序列进行排序,然后循环访问 x 标记,然后对于每个 x 标记,您进行搜索。这听起来很夸张。

我想知道我们应该如何在 CUDA 中做到这一点?谢谢。

可能有很多方法。 例如,您可以以线程并行方式使用基本的 cuda 二进制搜索。 我将演示一个推力实现。

出于此讨论的目的,我将假设两个数据集(已知采样点和所需采样点(都是任意放置的(即,我不假设任何一个样本集都是均匀分布的(。 但是,我将规定或要求所需的采样点完全包含在已知采样点内。 我相信这是明智的,因为通常的线性插值需要在所需采样点的两侧都有一个已知的采样点。

因此,我们将使用如下数据集:

   o:  1,3,7
f(o):  3,7,15
   x:  1.5, 2.5, 4.5, 5.0, 6.0, 6.5
f(x):    ?,   ?,   ?,   ?,   ?,   ?

我们看到f是已知的函数值,对应于 f(o) = 2o+1 ,在这种情况下是一条直线(尽管此方法不需要已知的样本点来拟合任何特定函数(。 x表示我们希望根据已知值(f(o)(插值函数值的索引。 我们的愿望是通过从最近的f(o)点插值来计算f(x)。 请注意,我们的数据集使得x的所有值都介于最小值 (1( 和最大值 (7( o值之间。 这是我之前说过的规定/要求。

我们的推力方法是使用矢量化的二叉搜索,使用 thrust::upper_bound 来定位每个所需x值适合o序列的"插入点"。 这为我们提供了右邻居和左邻居(右-1(进行插值。 一旦我们知道了插入点,如果我们想使用线性插值以外的其他东西,那么选择例如两个左邻和两个右邻(或更多(将是该算法的一个微不足道的扩展。

然后,插入点为我们提供左邻和右邻,我们使用此信息将x向量(所需的插值点(以及提供以下thrust::tuple(通过thrust::zip_iterator(传递给精心制作的thrust::transform

操作:
  • 右邻索引
  • 右邻功能值
  • 左邻索引
  • 左邻函数值

有了这些量,再加上所需的指数(x(,插值就很简单了。

编辑:受到另一个答案的启发,我决定包括一种避免并行二进制搜索的方法,而是使用前缀和方法来识别o数据中x数据的插入索引。 此方法假定xo序列已排序。

我们将从merge_by_key操作开始。 我们将xo合并,以建立排序(这似乎比二叉搜索更有效(。 xo数量将是"键",值将为 1 表示o,所有 0 表示x。 然后使用我们的示例数据,merge_by_key将生成以下内容:

o keys:  1,3,7
o vals:  1,1,1
x keys:  1.5,2.5,4.5,5.0,6.0,6.5
x vals:    0,  0,  0,  0,  0,  0
merged keys:  1, 1.5, 2.5,   3, 4.5, 5.0, 6.0, 6.5,   7
merged vals:  1,   0,   0,   1,   0,   0,   0,   0,   1

当我们对合并的 val 进行前缀和(包含扫描(时,我们得到:

ins. ind.:    1,   1,   1,   2,   2,   2,   2,   2,   3

然后,我们可以执行copy_if操作,仅提取与x val(其合并的 val 为零(关联的插入索引,以生成步骤 1 中生成的相同插入索引序列:

 d_i:  1, 1, 2, 2, 2, 2
然后,方法

2的其余部分可以使用与方法1完全相同的剩余插值码(thrust::transform(。

下面是一个完整的示例,显示了这两种方法:

$ cat t1224.cu
#include <thrust/device_vector.h>
#include <thrust/binary_search.h>
#include <thrust/transform.h>
#include <thrust/copy.h>
#include <thrust/iterator/zip_iterator.h>
#include <thrust/iterator/permutation_iterator.h>
#include <thrust/iterator/transform_iterator.h>
#include <iostream>
#include <thrust/merge.h>
#include <thrust/iterator/constant_iterator.h>
#include <thrust/scan.h>
struct interp_func
{
  template <typename T>
  __host__ __device__
  float operator()(float t1, T t2){  // m = (y1-y0)/(x1-x0)  y = m(x-x0) + y0
    return ((thrust::get<1>(t2) - thrust::get<3>(t2))/(thrust::get<0>(t2) - thrust::get<2>(t2)))*(t1 - thrust::get<2>(t2)) + thrust::get<3>(t2);
    }
};
using namespace thrust::placeholders;
int main(){
  // sample data
  float o[] = {1.0f, 3.0f, 7.0f}; // unevenly spaced sample points for function f
  float f[] = {3.0f, 7.0f, 15.0f}; // f(o) = 2o+1
  float x[] = {1.5f, 2.5f, 4.5f, 5.0f, 6.0f, 6.5f}; // additional desired sample points for f
  int so = sizeof(o)/sizeof(o[0]);
  int sx = sizeof(x)/sizeof(x[0]);
  // setup data on device
  thrust::device_vector<float> d_o(o, o+so);
  thrust::device_vector<float> d_f(f, f+so);
  thrust::device_vector<float> d_x(x, x+sx);
  thrust::device_vector<int>   d_i(sx); // insertion indices
  thrust::device_vector<float> d_r(sx); // results
  // method 1: binary search
  // perform search for insertion indices
  thrust::upper_bound(d_o.begin(), d_o.end(), d_x.begin(), d_x.end(), d_i.begin());
  // then perform linear interpolation based on left and right neighbors
  std::cout << "Method 1 insertion indices:" << std::endl;
  thrust::copy(d_i.begin(), d_i.end(), std::ostream_iterator<int>(std::cout, ","));
  std::cout << std::endl;
  thrust::transform(d_x.begin(), d_x.end(), thrust::make_zip_iterator(thrust::make_tuple(thrust::make_permutation_iterator(d_o.begin(), d_i.begin()), thrust::make_permutation_iterator(d_f.begin(), d_i.begin()), thrust::make_permutation_iterator(d_o.begin(), thrust::make_transform_iterator(d_i.begin(), _1-1)), thrust::make_permutation_iterator(d_f.begin(), thrust::make_transform_iterator(d_i.begin(), _1-1)))), d_r.begin(), interp_func());
  // output results
  std::cout << "Interpolation points:" << std::endl;
  thrust::copy(d_x.begin(), d_x.end(), std::ostream_iterator<float>(std::cout, ","));
  std::cout << std::endl << "Interpolated values:" << std::endl;
  thrust::copy(d_r.begin(), d_r.end(), std::ostream_iterator<float>(std::cout, ","));
  std::cout << std::endl << "Expected values:" << std::endl;
  for (int i = 0; i < sx; i++) std::cout << 2*x[i]+1 <<  ",";
  std::cout << std::endl;
  //method 2: merge + prefix sum
  thrust::device_vector<float> d_kr(sx+so);
  thrust::device_vector<int> d_vr(sx+so);
  thrust::device_vector<int> d_s(sx+so);
  thrust::merge_by_key(d_o.begin(), d_o.end(), d_x.begin(), d_x.end(), thrust::constant_iterator<int>(1), thrust::constant_iterator<int>(0), d_kr.begin(), d_vr.begin());
  thrust::inclusive_scan(d_vr.begin(), d_vr.end(), d_s.begin());
  thrust::copy_if(d_s.begin(), d_s.end(), d_vr.begin(), d_i.begin(), _1 == 0);
  std::cout << "Method 2 insertion indices:" << std::endl;
  thrust::copy(d_i.begin(), d_i.end(), std::ostream_iterator<int>(std::cout, ","));
  std::cout << std::endl;
  // remainder of solution method would be identical to end of method 1 starting with the thrust::transform
  return 0;
}
$ nvcc -o t1224 t1224.cu
$ ./t1224
Method 1 insertion indices:
1,1,2,2,2,2,
Interpolation points:
1.5,2.5,4.5,5,6,6.5,
Interpolated values:
4,6,10,11,13,14,
Expected values:
4,6,10,11,13,14,
Method 2 insertion indices:
1,1,2,2,2,2,
$

同样,一旦我们知道插入点,选择 2 个右邻和 2 个左邻域进行更多涉及的插值将是一个微不足道的扩展。 我们只需修改传递给转换(插值(函子的 zip 迭代器,并修改函子本身以实现所需的算术。

另请注意,此方法假定输入o序列已排序。 如果不是,则有必要添加具有f(值(的o(键(的按键排序。 对于方法 1,不需要对x序列进行排序,但必须对方法 2 进行排序(合并要求对两个序列进行排序(。

最佳方法的细节取决于所涉及的大小(即它是大量短序列还是单个巨大序列等(,但在高级别上,您可以仅使用输入序列的(并行潜在 O(N((排序和并行前缀总和来做到这一点。特别是,您可以避免任何二分搜索。查看现代GPU的"intervalExpand"背后的想法:https://nvlabs.github.io/moderngpu/intervalmove.html

简而言之,在伪代码中:

1:  sort the input sequence
2:  for each input point seq[i]: 
      let count[i] = number of output points in the interval [seq[i], seq[i+1])
3:  let indices = exclusive prefix-sum of count
4:  use intervalExpand() to go from seq, count, indices to the desired output.  

您可以在步骤4中坚持使用所需的任何插值公式,包括线性,三次等。 重要的是,intervalExpand会告诉你每个输出索引,哪些是正确的输入索引夹在输出中间。

同样,如果您正在处理大量小序列,则二进制搜索实际上可能运行更快且更容易编写。 否则,您应该能够使用 modernGPU 库中的模板化代码相对轻松地执行此操作。

希望有帮助。

最新更新