带权重的推力直方图

Question

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

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

用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 对同一单元格内的粒子质量求和：

密度.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