Thrust Histogram with weights

Question

I want to compute the density of particles over a grid. Therefore, I have a vector that contains the cellID of each particle, as well as a vector with the given mass which does not have to be uniform.

I have taken the non-sparse example from Thrust to compute a histogram of my particles. However, to compute the density, I need to include the weight of each particle, instead of simply summing the number of particles per cell, i.e. I'm interested in rho[i] = sum W[j] for all j that satify cellID[j]=i (probably unnecessary to explain, since everybody knows that).

Implementing this with Thrust has not worked for me. I also tried to use a CUDA kernel and thrust_raw_pointer_cast, but I did not succeed with that either.

EDIT:

Here is a minimal working example which should compile via nvcc file.cu under CUDA 6.5 and with Thrust installed.

#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
     *
     *
     */
}

I hope the comment at the end of the file also makes clear what I mean by computing the density. If there is any question open, please feel free to ask. Thanks!

There still seems to be a problem in understanding my problem, which I am sorry for! Therefore I added some pictures. Consider the first picture. For my understanding, a histogram would simply be the count of particles per grid cell. In this case a histogram would be an array of size 36, since there are 36 cells. Also, there would be a lot of zero entries in the vector, since for example in the upper left corner almost no cell contains a particle. This is what I already have in my code.

Now consider the slightly more complicated case. Here each particle has a different mass, indicated by the different size in the plot. To compute the density I can't just add the number of particles per cell, but I have to add the mass of all particles per cell. This is what I'm unable to implement.

Answer 1

What you described in your example does not look like a histogram but rather like a segmented reduction.

The following example code uses thrust::reduce_by_key to sum up the masses of particles within the same cell:

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);
}

compile using

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

output

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