使用 Thrust 从单个数组构建多个直方图

Using Thrust to construct multiple histograms from a single array

假设我有一个数组,用于存储多个事件的时间戳。例如,T1_e1, T2_e1,.....,T1_e2, T2_e2, T3_e2,.....T1_eN, T2 ,eN,..

我知道Thrust 提供了计算相邻差异的函数,但这里我需要为多个事件执行此操作。基本上,从单个输入数组构建多个直方图。

所以输出会有 N 个不同的直方图(每个事件一个),如下所示: e1 的直方图 bins,e2 的直方图 bins,e3 的直方图 bins,....eN 的直方图 bins。

Input1 (timestamps): 100, 101, 104, 105, 101,104, 106, 111, 90, 91, 93,  94,95
Input2 (events):    4123,4123,4123,4123,2129,2129,2129,2129,300,300,300,300,300

output: 4123:(1,2),(2,0),(3,1),(4,0),(5,0)
        2129:(1,0),(2,1),(3,1),(4,0),(5,1)
        300: (1,2),(2,1),(3,0),(4,),(5,0)  

bin 的数量将是固定的,即每个直方图 5 个 bin。 关于元组:(x,y) -> x 是属于同一事件的两个连续时间戳之间的差异。 y 是计数。

如果我们考虑事件4123,第一个元组是(1,2),因为101和100的差是1,105和104是1。所以有两个时间戳差属于这个bin , 因此 (1,2).

有人可以建议最有效的方法吗?到目前为止,似乎我将不得不编写自己的代码。但如果有现成的解决方案,我想先尝试一下。

这是计算稀疏直方图的一种可能方法(即仅非零分箱)。使用推力散射将稀疏直方图转换为密集直方图(包括零和非零区间)应该不难。

  1. 使用 thrust::adjacent_difference 和一个特殊函子 (my_adj_diff) 计算时间戳差异,该函子仅计算类似事件的相邻差异。

  2. 使用 thrust::remove_if 删除每个事件序列(由第 1 步中的仿函数创建)开始处的零值(每个事件一个)。

  3. 将事件和时间戳差异合并为一个整数,以便使用 thrust::transform 进行直方图绘制。在我的示例中,这只是将事件乘以 100,并添加时间戳差异(假设是一组小于 100 最大 bin 的 bin)。

  4. 从推力histogram example code.

  5. 使用稀疏直方图的方法

这是一个完整的示例。数据与您的预期输出不完全匹配(仅限非零值),因为您的预期输出有误。

$ cat t678.cu
#include <thrust/adjacent_difference.h>
#include <thrust/iterator/zip_iterator.h>
#include <thrust/iterator/constant_iterator.h>
#include <thrust/device_vector.h>
#include <thrust/host_vector.h>
#include <thrust/copy.h>
#include <thrust/remove.h>
#include <thrust/sort.h>
#include <thrust/inner_product.h>
#include <thrust/reduce.h>
#include <thrust/functional.h>


#include <iostream>

#define ZIP(X,Y) thrust::make_zip_iterator(thrust::make_tuple(X,Y))
#define SCALE 100

struct my_adj_diff
{
  template <typename T>
  __host__ __device__
  T operator()(T &d2, T &d1) const
    {
      if (thrust::get<1>(d1) == thrust::get<1>(d2)) {
        thrust::get<0>(d2) -= thrust::get<0>(d1);}
      else {
        thrust::get<0>(d2) = 0;}
      return d2;
    }
};

struct my_is_zero
{
  template <typename T>
  __host__ __device__
  bool operator()(const T &d1) const
    {
      return (thrust::get<0>(d1) == 0);
    }
};

struct my_combine
{
  template <typename T>
  __host__ __device__
  T operator()(const T &d1, const T &d2) const
    {
      return (d1*SCALE)+d2;
    }
};

// sparse histogram using reduce_by_key
// modified from: https://github.com/thrust/thrust/blob/master/examples/histogram.cu
template <typename Vector1,
          typename Vector2,
          typename Vector3>
void sparse_histogram(const Vector1& input,
                            Vector2& histogram_values,
                            Vector3& histogram_counts)
{
  typedef typename Vector1::value_type ValueType; // input value type
  typedef typename Vector3::value_type IndexType; // histogram index type
  thrust::device_vector<ValueType> data(input);

  // sort data to bring equal elements together
  thrust::sort(data.begin(), data.end());

  // number of histogram bins is equal to number of unique values (assumes data.size() > 0)
  IndexType num_bins = thrust::inner_product(data.begin(), data.end() - 1,
                                             data.begin() + 1,
                                             IndexType(1),
                                             thrust::plus<IndexType>(),
                                             thrust::not_equal_to<ValueType>());

  // resize histogram storage
  histogram_values.resize(num_bins);
  histogram_counts.resize(num_bins);

  // compact find the end of each bin of values
  thrust::reduce_by_key(data.begin(), data.end(),
                        thrust::constant_iterator<IndexType>(1),
                        histogram_values.begin(),
                        histogram_counts.begin());

}


int main(){

  int tstamps[] = { 100, 101, 104, 105, 101,104, 106, 111, 90, 91, 93,  94,95 };
  int mevents[] = {4123,4123,4123,4123,2129,2129,2129,2129,300,300,300,300,300};
  int dsize = sizeof(tstamps)/sizeof(int);

  thrust::host_vector<int>   h_stamps(tstamps, tstamps+dsize);
  thrust::host_vector<int>   h_events(mevents, mevents+dsize);
  thrust::device_vector<int> d_stamps = h_stamps;
  thrust::device_vector<int> d_events = h_events;
  thrust::device_vector<int> diffs(dsize);
  // compute timestamp differences by event
  thrust::adjacent_difference(ZIP(d_stamps.begin(), d_events.begin()), ZIP(d_stamps.end(), d_events.end()), ZIP(d_stamps.begin(), d_events.begin()), my_adj_diff());
  d_stamps[0] = 0;  // fix up first event for adjacent_difference
  int sz1 = thrust::remove_if(ZIP(d_stamps.begin(), d_events.begin()), ZIP(d_stamps.end(), d_events.end()), my_is_zero()) - ZIP(d_stamps.begin(), d_events.begin());
  d_stamps.resize(sz1);
  d_events.resize(sz1);
  // pack events and timestamps into a single vector - assumes max bin (time difference) is less than SCALE
  thrust::device_vector<int> d_data(sz1);
  thrust::transform(d_events.begin(), d_events.end(), d_stamps.begin(), d_data.begin(), my_combine());
  // compute histogram
  thrust::device_vector<int> histogram_values;
  thrust::device_vector<int> histogram_counts;
  sparse_histogram(d_data, histogram_values, histogram_counts);
  thrust::copy(histogram_values.begin(), histogram_values.end(), std::ostream_iterator<int>(std::cout, ","));
  std::cout << std::endl;
  thrust::copy(histogram_counts.begin(), histogram_counts.end(), std::ostream_iterator<int>(std::cout, ","));
  std::cout << std::endl;
}
$ nvcc t678.cu -o t678
$ ./t678
30001,30002,212902,212903,212905,412301,412303,
3,1,1,1,1,2,1,
$

请注意,您需要使用 CUDA 7.0(不是 CUDA 7.0 RC 或任何更早版本的 CUDA),或者在尝试使用带 [=11= 的 zip 迭代器时下载最新的 thrust master 分支 from github, because older versions of thrust have an issue ].