可能有多种方法。例如,您可以以线程并行的方式使用基本的cuda binary search。我将演示 thrust 的实现。
出于本次讨论的目的,我假设两个数据集(已知样本点和所需样本点)是任意放置的(即,我不假设任何一个样本集是均匀分布的)。但是,我将规定或要求所需的样本点完全包含在已知的样本点内。我相信这是明智的,因为通常的线性插值需要在所需采样点的任一侧都有一个已知的采样点。
因此我们将使用这样的数据集:
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 序列的“插入点”。这为我们提供了我们的右邻居和左邻居 (right-1) 用于插值。一旦我们知道了插入点,选择例如这个算法将是一个简单的扩展。 two 左边和 two 右边的邻居(或更多),如果我们想使用线性插值以外的东西。
然后插入点为我们提供左右邻居,我们使用此信息将x 向量(所需的插值点)和thrust::tuple(通过@987654338 @) 提供:
有了这些数量,加上所需的索引 (x),插值就很简单了。
编辑: 受另一个答案的启发,我决定包含一种避免并行二进制搜索的方法,而是使用前缀和方法来识别 x 数据的插入索引在o 数据中。此方法假定 x 和 o 序列都已排序。
我们将从 merge_by_key 操作开始。我们将x 与o 合并,以建立排序(这似乎比二分查找更有效)。 x 和 o 数量将是“键”,o 的值全为 1,x 的值全为 0。然后使用我们的示例数据,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 vals(其合并的 vals 为零)相关的插入索引,以生成与步骤 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 序列已排序。如果不是,则有必要添加o(键)和f(值)的键排序。 x 序列不需要为方法 1 排序,但必须为方法 2 排序(合并要求两个序列都排序)。