Гистограмма тяги с весами

Question

Я хочу вычислить плотность частиц по сетке. Поэтому у меня есть вектор, содержащий cellID каждой частицы, а также вектор с заданным mass, который не обязательно должен быть однородным.

Я вычислил гистограмму моих частиц с помощью нерезкого примера от Thrust. Однако, чтобы вычислить плотность, мне нужно включить вес каждой частицы вместо простого суммирования числа частиц на ячейку, то есть меня интересует rho[i] = sum W[j] для всех j, которые удовлетворяют cellID[j]=i (вероятно, не нужно объяснять, поскольку все знает это).

Реализация этого с Thrust не сработало для меня. Я также пытался использовать ядро ​​CUDA и thrust_raw_pointer_cast, но мне тоже этого не удалось.

EDIT:

Вот минимальный рабочий пример, который должен составить через nvcc file.cu под CUDA 6.5 и с 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