带权重的推力直方图

Question

我想计算网格上粒子的密度。因此，我有一个包含每个粒子的 cellID 的矢量，以及一个具有给定 mass 的矢量，它不必是均匀的。

我使用了 Thrust 中的非稀疏示例来计算我的粒子的直方图。但是，要计算密度，我需要包括每个粒子的权重，而不是简单地将每个单元格的粒子数相加，即我对满足 j 的所有 rho[i] = sum W[j] 感兴趣16=]（可能不需要解释，因为每个人都知道）。

用 Thrust 实现这个对我没有用。我还尝试使用 CUDA 内核和 thrust_raw_pointer_cast，但我也没有成功。

编辑：

这是一个最小的工作示例，应该在 CUDA 6.5 下通过 nvcc file.cu 编译并安装 Thrust。

#include <thrust/device_vector.h>
#include <thrust/sort.h>
#include <thrust/copy.h>
#include <thrust/binary_search.h>
#include <thrust/adjacent_difference.h>


// Predicate
struct is_out_of_bounds {
    __host__ __device__ bool operator()(int i) {

        return (i < 0); // out of bounds elements have negative id;
    }
};



// cf.: https://code.google.com/p/thrust/source/browse/examples/histogram.cu, but modified
template<typename T1, typename T2>
void computeHistogram(const T1& input, T2& histogram) {
    typedef typename T1::value_type ValueType; // input value type
    typedef typename T2::value_type IndexType; // histogram index type

    // copy input data (could be skipped if input is allowed to be modified)
    thrust::device_vector<ValueType> data(input);

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


    // there are elements that we don't want to count, those have ID -1;
    data.erase(thrust::remove_if(data.begin(), data.end(), is_out_of_bounds()),data.end());



    // number of histogram bins is equal to the maximum value plus one
    IndexType num_bins = histogram.size();

    // find the end of each bin of values
    thrust::counting_iterator<IndexType> search_begin(0);
    thrust::upper_bound(data.begin(), data.end(), search_begin,
            search_begin + num_bins, histogram.begin());

    // compute the histogram by taking differences of the cumulative histogram
    thrust::adjacent_difference(histogram.begin(), histogram.end(),
            histogram.begin());

}




int main(void) {



    thrust::device_vector<int> cellID(5);
    cellID[0] = -1; cellID[1] = 1; cellID[2] = 0; cellID[3] = 2; cellID[4]=1;
    thrust::device_vector<float> mass(5);
    mass[0] = .5; mass[1] = 1.0; mass[2] = 2.0; mass[3] = 3.0; mass[4] = 4.0;

    thrust::device_vector<int> histogram(3);
    thrust::device_vector<float> density(3);


    computeHistogram(cellID,histogram);

    std::cout<<"\nHistogram:\n";

    thrust::copy(histogram.begin(), histogram.end(),
                    std::ostream_iterator<int>(std::cout, " "));
    std::cout << std::endl;
    // this will print: " Histogram 1 2 1 " 
    // meaning one element with ID 0, two elements with ID 1 
    // and one element with ID 2

    /* here is what I am unable to implement:
     *
     *
     * computeDensity(cellID,mass,density);
     *
     * print(density):        2.0 5.0 3.0
     *
     *
     */
}

我希望文件末尾的评论也能清楚地说明我计算密度的意思。如果有任何问题未解决，请随时提问。谢谢！

我的问题理解起来好像还是有问题，很抱歉！因此我添加了一些图片。考虑第一张图片。据我了解，直方图只是每个网格单元的粒子数。在这种情况下，直方图将是一个大小为 36 的数组，因为有 36 个单元格。此外，向量中会有很多零条目，因为例如在左上角几乎没有单元格包含粒子。这是我的代码中已有的。

现在考虑稍微复杂一点的情况。这里每个粒子都有不同的质量，由图中不同的大小表示。要计算密度，我不能只添加每个单元格的粒子数，但我必须添加每个单元格所有粒子的质量。这是我无法实现的。

Answer 1

您在示例中描述的内容看起来不像直方图，而更像是分段缩减。

以下示例代码使用 thrust::reduce_by_key 对同一单元格内的粒子质量求和：

density.cu

#include <thrust/device_vector.h>
#include <thrust/sort.h>
#include <thrust/reduce.h>
#include <thrust/copy.h>
#include <thrust/scatter.h>
#include <iostream>

#define PRINTER(name) print(#name, (name))
template <template <typename...> class V, typename T, typename ...Args>
void print(const char* name, const V<T,Args...> & v)
{
    std::cout << name << ":\t\t";
    thrust::copy(v.begin(), v.end(), std::ostream_iterator<T>(std::cout, "\t"));
    std::cout << std::endl << std::endl;
}

int main()
{

    const int particle_count = 5;
    const int cell_count  = 10;

    thrust::device_vector<int> cellID(particle_count);
    cellID[0] = -1; cellID[1] = 1; cellID[2] = 0; cellID[3] = 2; cellID[4]=1;
    thrust::device_vector<float> mass(particle_count);
    mass[0] = .5; mass[1] = 1.0; mass[2] = 2.0; mass[3] = 3.0; mass[4] = 4.0;

    std::cout << "input data" << std::endl;
    PRINTER(cellID);
    PRINTER(mass);

    thrust::sort_by_key(cellID. begin(), cellID.end(), mass.begin());

    std::cout << "after sort_by_key" << std::endl;
    PRINTER(cellID);
    PRINTER(mass);

    thrust::device_vector<int> reduced_cellID(particle_count);
    thrust::device_vector<float> density(particle_count);

    int new_size = thrust::reduce_by_key(cellID. begin(), cellID.end(),
                                         mass.begin(),
                                         reduced_cellID.begin(),
                                         density.begin()
                                        ).second - density.begin();                                        
    if (reduced_cellID[0] == -1)
    {
        density.erase(density.begin());
        reduced_cellID.erase(reduced_cellID.begin());
        new_size--;
    }
    density.resize(new_size);
    reduced_cellID.resize(new_size);

    std::cout << "after reduce_by_key" << std::endl;

    PRINTER(density);
    PRINTER(reduced_cellID);

    thrust::device_vector<float> final_density(cell_count);
    thrust::scatter(density.begin(), density.end(), reduced_cellID.begin(), final_density.begin());
    PRINTER(final_density);
}

编译使用

nvcc -std=c++11 density.cu -o density

输出

input data
cellID:     -1   1  0   2   1   

mass:       0.5  1  2   3   4   

after sort_by_key
cellID:     -1   0  1   1   2   

mass:       0.5  2  1   4   3   

after reduce_by_key
density:        2   5   3   

reduced_cellID: 0   1   2   

final_density:  2   5   3   0   0   0   0   0   0   0

带权重的推力直方图

Thrust Histogram with weights

cuda

gpu

histogram

thrust