用CUDA压缩“稀疏数据”(CCL:连通元件标记约简)

Question

我有一个5百万个32位整数的列表(实际上是一个2048x2560图像)，这是90%的零。非零单元是完全不连续的标签(例如，2049,8195,1334300,34320923,4320932) (这是我们自定义连接组件标记CCL算法的输出)。我是一个NVIDA特斯拉K40工作，所以我希望它，如果这需要任何前缀扫描工作，它使用洗牌，选票或任何更高的CC功能。

我不需要一个完整的例子，只是一些建议。

为了说明，这里有一个博客被我们的CCL算法标记。

​

​

其他气泡将有一个不同的独特标签(例如13282)。但它们都将被零包围，并且在形状上都是椭球。(我们为椭球优化了CCL，这就是我们不使用库的原因)。但是一个副作用是blob标签不是连续的数字。我们不关心它们编号的顺序，但是我们希望一个标记为#1，另一个标记为#2，最后一个被标记为#n，其中n是图像中的小块数。

我说的标签#1是什么意思？我的意思是，所有的2242个细胞应该被一个1取代，所有的13282个细胞都有一个#2，等等。

我们的CCL的最大blob数等于2048x2560。所以我们知道数组的大小。

事实上，罗伯特·克罗夫拉一天前已经对此给出了一个很好的答案。这并不准确，但我现在明白了如何应用这个答案。所以我不需要任何帮助了。但是他在时间和精力上都很慷慨，让我用例子重写问题，所以我就这么做了。

Answer 1

一种可能的办法是使用以下顺序：

1. thrust::transform -将输入数据转换为所有1或0： 0 27 42 0 18 99 94 91 0-输入数据0 1 1 1 0-这将是我们的“掩码矢量”
2. thrust::inclusive_scan -将掩码矢量转换为一个渐进序列： 0 1 1 0 1 1 1 0-掩码向量0 1 2 2 3 4 5 6 6-序列向量
3. 另一个用来掩盖不增加值的thrust::transform： 0 1 1 1 0-掩码向量0 1 2 2 3 4 5 6 6-序列向量

注意，我们可以将前两个步骤与thrust::transform_inclusive_scan结合起来，然后将第三个步骤作为一个带有稍微不同的转换函子的thrust::transform来执行。这种修改使我们不必创建临时的“掩码”矢量。

下面是一个完整的示例，展示了使用thrust::transform_inclusive_scan进行“修改”的方法

$ cat t635.cu
#include <iostream>
#include <stdlib.h>

#include <thrust/device_vector.h>
#include <thrust/host_vector.h>
#include <thrust/transform.h>
#include <thrust/transform_scan.h>
#include <thrust/generate.h>
#include <thrust/copy.h>


#define DSIZE 20
#define PCT_ZERO 40

struct my_unary_op
{
  __host__ __device__
  int operator()(const int data) const
  {
    return (!data) ?  0:1;}
};

struct my_binary_op
{
  __host__ __device__
  int operator()(const int d1, const int d2) const
  {
    return (!d1) ? 0:d2;}
};

int main(){

// generate DSIZE random 32-bit integers, PCT_ZERO% are zero
  thrust::host_vector<int> h_data(DSIZE);
  thrust::generate(h_data.begin(), h_data.end(), rand);
  for (int i = 0; i < DSIZE; i++)
    if ((rand()%100)< PCT_ZERO) h_data[i] = 0;
    else h_data[i] %= 1000;
  thrust::device_vector<int> d_data = h_data;
  thrust::device_vector<int> d_result(DSIZE);
  thrust::transform_inclusive_scan(d_data.begin(), d_data.end(), d_result.begin(), my_unary_op(), thrust::plus<int>());
  thrust::transform(d_data.begin(), d_data.end(), d_result.begin(), d_result.begin(), my_binary_op());
  thrust::copy(d_data.begin(), d_data.end(), std::ostream_iterator<int>(std::cout, ","));
  std::cout << std::endl;
  thrust::copy(d_result.begin(), d_result.end(), std::ostream_iterator<int>(std::cout, ","));
  std::cout << std::endl;
  return 0;
}

$ nvcc -o t635 t635.cu
$ ./t635
0,886,777,0,793,0,386,0,649,0,0,0,0,59,763,926,540,426,0,736,
0,1,2,0,3,0,4,0,5,0,0,0,0,6,7,8,9,10,0,11,
$

对更新的响应，这些新的信息使问题更难解决，在我看来。组织语法技术出现在脑海中，但没有对32位整数(标签)所占范围的任何限制或对特定标签在数据集中重复的次数的任何限制，直方图技术似乎是不切实际的。这促使我考虑对数据进行排序。

这样的方法应该能奏效：

1. 使用thrust::sort对数据进行排序。
2. 使用thrust::unique删除重复项。
3. 删除重复项的排序数据现在为输出集0、1、2、.提供了排序。让我们把这叫做“地图”。我们可以使用parallel binary-search technique将原始数据集中的每个标签转换为其映射的输出值。

在我看来，这个过程相当“昂贵”。我建议重新考虑上游标签的运作，看看是否可以重新设计，以产生一个更适合有效下游处理的数据集。

总之，这里有一个充分发挥作用的例子：

$ cat t635.cu
#include <iostream>
#include <stdlib.h>

#include <thrust/device_vector.h>
#include <thrust/host_vector.h>
#include <thrust/transform.h>
#include <thrust/generate.h>
#include <thrust/sort.h>
#include <thrust/unique.h>
#include <thrust/copy.h>


#define DSIZE 20
#define PCT_ZERO 40
#define RNG 10

#define nTPB 256

// sets idx to the index of the first element in a that is
// equal to or larger than key

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

__global__ void find_my_idx(const int *a, const unsigned len_a,  int *my_data, int *my_idx, const unsigned len_data){
  unsigned idx = (blockDim.x * blockIdx.x) + threadIdx.x;
  if (idx < len_data){
    unsigned sp_a;
    int val = my_data[idx];
    bsearch_range(a, val, len_a, &sp_a);
    my_idx[idx] = sp_a;
    }
}


int main(){

// generate DSIZE random 32-bit integers, PCT_ZERO% are zero
  thrust::host_vector<int> h_data(DSIZE);
  thrust::generate(h_data.begin(), h_data.end(), rand);
  for (int i = 0; i < DSIZE; i++)
    if ((rand()%100)< PCT_ZERO) h_data[i] = 0;
    else h_data[i] %= RNG;
  thrust::device_vector<int> d_data = h_data;
  thrust::device_vector<int> d_result = d_data;
  thrust::sort(d_result.begin(), d_result.end());
  thrust::device_vector<int> d_unique = d_result;
  int unique_size = thrust::unique(d_unique.begin(), d_unique.end()) - d_unique.begin();
  find_my_idx<<< (DSIZE+nTPB-1)/nTPB , nTPB >>>(thrust::raw_pointer_cast(d_unique.data()), unique_size, thrust::raw_pointer_cast(d_data.data()), thrust::raw_pointer_cast(d_result.data()), DSIZE);

  thrust::copy(d_data.begin(), d_data.end(), std::ostream_iterator<int>(std::cout, ","));
  std::cout << std::endl;
  thrust::copy(d_result.begin(), d_result.end(), std::ostream_iterator<int>(std::cout, ","));
  std::cout << std::endl;
  return 0;
}
$ nvcc t635.cu -o t635
$ ./t635
0,6,7,0,3,0,6,0,9,0,0,0,0,9,3,6,0,6,0,6,
0,2,3,0,1,0,2,0,4,0,0,0,0,4,1,2,0,2,0,2,
$

Answer 2

我的回答类似@RobertCrovella给出的答案，但我认为使用thrust::lower_bound代替自定义二进制搜索更简单。(现在它是纯推力，后端可以互换)

1. 复制输入数据
2. 对复制的数据进行排序
3. 根据已排序的数据创建唯一的列表
4. 在唯一列表中找到每个输入的下界

我在下面列举了一个完整的例子。有趣的是，通过对thrust::unique的另一次调用，预挂排序步骤可以使进程变得更快。根据输入数据的不同，这可以极大地减少排序中的元素数量，这是这里的瓶颈。

#include <iostream>
#include <stdlib.h>
#include <thrust/device_vector.h>
#include <thrust/host_vector.h>
#include <thrust/transform.h>
#include <thrust/generate.h>
#include <thrust/sort.h>
#include <thrust/unique.h>
#include <thrust/binary_search.h>
#include <thrust/copy.h>

int main()
{
  const int ndata = 20;
  // Generate host input data
  thrust::host_vector<int> h_data(ndata);
  thrust::generate(h_data.begin(), h_data.end(), rand);
  for (int i = 0; i < ndata; i++)
  {
    if ((rand() % 100) < 40)
      h_data[i] = 0;
    else
      h_data[i] %= 10;
  }

  // Copy data to the device
  thrust::device_vector<int> d_data = h_data;
  // Make a second copy of the data
  thrust::device_vector<int> d_result = d_data;
  // Sort the data copy
  thrust::sort(d_result.begin(), d_result.end());
  // Allocate an array to store unique values
  thrust::device_vector<int> d_unique = d_result;
  {
    // Compress all duplicates
    const auto end = thrust::unique(d_unique.begin(), d_unique.end());
    // Search for all original labels, in this compressed range, and write their
    // indices back as the result
    thrust::lower_bound(
      d_unique.begin(), end, d_data.begin(), d_data.end(), d_result.begin());
  }

  thrust::copy(
    d_data.begin(), d_data.end(), std::ostream_iterator<int>(std::cout, ","));
  std::cout << std::endl;
  thrust::copy(d_result.begin(),
               d_result.end(),
               std::ostream_iterator<int>(std::cout, ","));
  std::cout << std::endl;
  return 0;
}