Повторная выборка 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(x) через интерполяцию от ближайшего f(o) точки. Обратите внимание, что наш набор данных таков, что все значения x лежат между минимумом (1) и максимумом (7) o ценности. Это условие / требование, которое я изложил ранее.

Наш метод тяги будет использовать векторизованный двоичный поиск, используя thrust::upper_bound, чтобы найти «точку вставки», где каждый желаемый x значение вписывается в o последовательность. Это дает нам нашего правого соседа и левого соседа (справа-1) для интерполяции. Как только мы узнаем точку вставки, было бы тривиальным расширением этого алгоритма выбрать, например, два левый и два правильные соседи (или больше), если мы хотим использовать что-то кроме линейной интерполяции.

Точка вставки затем дает нам наших левых и правых соседей, и мы используем эту информацию для передачи соответствующим образом созданным thrust::transform операция x вектор (желаемые точки интерполяции) вместе с thrust::tuple (с помощью thrust::zip_iterator) который обеспечивает:

• индекс правого соседа
• функциональное значение правого соседа
• указатель левого соседа
• функциональное значение левого соседа

С этими величинами плюс желаемый индекс (x), интерполяция проста.

РЕДАКТИРОВАТЬ: Вдохновленный другим ответом, я решил включить метод, который избегает параллельного двоичного поиска, но вместо этого использует метод суммы префикса для определения индексов вставки для x данные в o данные. Этот метод предполагает и то и другое x а также o последовательности отсортированы.

Мы начнем с операции merge_by_key. Мы будем сливаться x с o, чтобы установить порядок (это кажется более эффективным, чем бинарный поиск). x а также o количества будут «ключами», а значения будут все 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

Когда мы делаем префиксную сумму (включающее сканирование) для объединенных значений, мы получаем:

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

Затем мы можем выполнить операцию copy_if, чтобы извлечь только индексы вставки, связанные с x vals (чьи объединенные значения равны нулю) для создания той же последовательности индекса вставки, которая была получена на шаге 1:

 d_i:  1, 1, 2, 2, 2, 2

Оставшаяся часть метода 2 может затем использовать точно такой же оставшийся интерполяционный код (thrust::transform), как это было использовано в методе 1.

Вот полностью проработанный пример, показывающий оба метода:

$ 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-итератор, передаваемый функтору transform (interpolation), и изменим сам функтор для реализации желаемой арифметики.

Также обратите внимание, что этот метод предполагает ввод o последовательность уже отсортирована. Если это не так, то необходимо добавить сортировку по ключу o (ключи) с f (ценности). x последовательность не должна быть отсортирована для метода 1, но должна быть отсортирована для метода 2 (объединение требует, чтобы обе последовательности были отсортированы).