How to do weighted average of two device_vectors with points selected using map in thrust?

Question

I have two device_vector P & Q (say of size 100). I have two device_vector maps(MapP & MapQ of size say 10) for P & Q which has the indices of points to be selected from P & Q. I have a device_vector D for weight.

I need to compute (P*D+Q)/(D+1) for all points from P & Q which have been selected using the respective maps.

My method is as below. It works, but is too cumbersome. Can anyone suggest a better way to do it?

#include <thrust/device_vector.h>
#include <thrust/random.h>
#include <thrust/sequence.h>
#include <thrust/execution_policy.h>
#include <thrust/transform.h>
#include <thrust/functional.h>
#include <thrust/iterator/transform_iterator.h>
#include <thrust/iterator/permutation_iterator.h>

thrust::device_vector<float> random_vector(const size_t N, 
                                         unsigned int seed = thrust::default_random_engine::default_seed)
{
    thrust::default_random_engine rng(seed);
    thrust::uniform_real_distribution<float> u01(0.0f, 10.0f);
    thrust::device_vector<float> temp(N);
    for(size_t i = 0; i < N; i++) {
        temp[i] = u01(rng);
    }
    return temp;
}

// note: functor inherits from unary_function
struct increment : public thrust::unary_function<int,int>
{
  __host__ __device__
  int operator()(int x) const
  {
    return x + 1;
  }
};

int main(int argc, char * argv[])
{

int N=atoi(argv[1]);

thrust::device_vector<float> P = random_vector(N,1);
thrust::device_vector<float> Q = random_vector(N,9);

thrust::device_vector<int> D(N);
thrust::sequence(thrust::device, D.begin(), D.begin() + N, 1);

thrust::device_vector<float> temp(10);

thrust::device_vector<int> MapP(10); // map
thrust::device_vector<int> MapQ(10); // map

MapP[0]=0;MapP[1]=5;MapP[2]=4;MapP[3]=2;MapP[4]=7;MapP[5]=1;MapP[6]=9;MapP[7]=3;MapP[8]=6;MapP[9]=8;
MapQ[0]=10;MapQ[1]=15;MapQ[2]=12;MapQ[3]=14;MapQ[4]=11;MapQ[5]=17;MapQ[6]=13;MapQ[7]=19;MapQ[8]=18;MapQ[9]=16;


// The weighted average is (D*P+Q)/(D+1)
// We compute D*P first

//thrust::transform(thrust::device, P.begin(), P.end(), D.begin(), temp.begin(), thrust::multiplies<float>()); // use permutation iterator

thrust::transform(thrust::device, thrust::make_permutation_iterator(P.begin(),MapP.begin()),
                                  thrust::make_permutation_iterator(P.end(),MapP.end()),
                  thrust::make_permutation_iterator(D.begin(),MapP.begin()), 
                  temp.begin(), thrust::multiplies<float>());


// Then we add D*p to Q

//thrust::transform(thrust::device, temp.begin(), temp.end(), Q.begin(), temp.begin(), thrust::plus<float>()); // use permutation iterator

thrust::transform(thrust::device, temp.begin(), temp.end(),
                  thrust::make_permutation_iterator(Q.begin(),MapQ.begin()), 
                  temp.begin(), thrust::plus<float>());


// Then we divide by D+1

//thrust::transform(thrust::device, temp.begin(), temp.end(), thrust::make_transform_iterator(D.begin(), increment()), temp.begin(),  thrust::divides<float>());

thrust::transform(thrust::device, temp.begin(), temp.end(),
                  thrust::make_permutation_iterator(D.begin(),MapP.begin()), 
                  temp.begin(), thrust::divides<float>());


// replace contents of P with the weighted sum using pts in map M

thrust::copy(thrust::device, temp.begin(), temp.end(), thrust::make_permutation_iterator(P.begin(),MapP.begin())); // use permutation iterator

return 0;
}

What does `(P*D+Q)/(D+1)` mean mathematically? Is this an element-wise operation? If P and D are vectors then P*D is either a scalar or a matrix depending on whether `*` denotes the inner or outer product. And division isn't defined.... — talonmies, Nov 23 '18 at 14:07

Robert Crovella · Accepted Answer · 2018-11-23T14:20:36.117

I'll assume you want element-wise operation on the vectors, since that is the behavior of your supplied demonstrator code.

Note that when passing the end of a permutation iterator, we do not use the end of the source vector:

thrust::make_permutation_iterator(P.end(),MapP.end()),
                                  ^^^^^

but instead the beginning:

thrust::make_permutation_iterator(P.begin(),MapP.end()),

Refer to the thrust quick start guide for an example of this.

Also note that both in your question and your code you refer to dividing by D+1, but your code is actually dividing by D, not D+1.

Regarding your question, everything can be done with a single call to thrust::transform with an appropriately defined functor. Since we need to pass multiple vectors to thrust::transform in this realization, we will introduce thrust::zip_iterator

$ cat t332.cu
#include <thrust/device_vector.h>
#include <thrust/random.h>
#include <thrust/sequence.h>
#include <thrust/execution_policy.h>
#include <thrust/transform.h>
#include <thrust/functional.h>
#include <thrust/iterator/transform_iterator.h>
#include <thrust/iterator/permutation_iterator.h>
#include <thrust/iterator/zip_iterator.h>

thrust::device_vector<float> random_vector(const size_t N,
                                         unsigned int seed = thrust::default_random_engine::default_seed)
{
    thrust::default_random_engine rng(seed);
    thrust::uniform_real_distribution<float> u01(0.0f, 10.0f);
    thrust::device_vector<float> temp(N);
    for(size_t i = 0; i < N; i++) {
        temp[i] = u01(rng);
    }
    return temp;
}

// The weighted average is (D*P+Q)/(D+1)
struct w_avg
{
template <typename T>
  __host__ __device__
  float operator()(T x) const
  {
    return (thrust::get<0>(x)*thrust::get<1>(x)+thrust::get<2>(x))/(thrust::get<1>(x)+1.0f);
  }
};

int main(int argc, char * argv[])
{

int N=atoi(argv[1]);

thrust::device_vector<float> P = random_vector(N,1);
thrust::device_vector<float> Q = random_vector(N,9);

thrust::device_vector<int> D(N);
thrust::sequence(thrust::device, D.begin(), D.begin() + N, 1);


thrust::device_vector<int> MapP(10); // map
thrust::device_vector<int> MapQ(10); // map

MapP[0]=0;MapP[1]=5;MapP[2]=4;MapP[3]=2;MapP[4]=7;MapP[5]=1;MapP[6]=9;MapP[7]=3;MapP[8]=6;MapP[9]=8;
MapQ[0]=10;MapQ[1]=15;MapQ[2]=12;MapQ[3]=14;MapQ[4]=11;MapQ[5]=17;MapQ[6]=13;MapQ[7]=19;MapQ[8]=18;MapQ[9]=16;


// The weighted average is (D*P+Q)/(D+1)

thrust::transform(thrust::device, thrust::make_zip_iterator(thrust::make_tuple(
                                                            thrust::make_permutation_iterator(P.begin(),MapP.begin()),
                                                            thrust::make_permutation_iterator(D.begin(),MapP.begin()),
                                                            thrust::make_permutation_iterator(Q.begin(),MapQ.begin()))),
                                  thrust::make_zip_iterator(thrust::make_tuple(
                                                            thrust::make_permutation_iterator(P.begin(),MapP.end()),
                                                            thrust::make_permutation_iterator(D.begin(),MapP.end()),
                                                            thrust::make_permutation_iterator(Q.begin(),MapQ.end()))),
                                                            thrust::make_permutation_iterator(P.begin(),MapP.begin()),
                                  w_avg());


for (int i = 0; i < 5; i++) {
  std::cout << P[i] << std::endl;}
return 0;
}
$ nvcc -o t332 t332.cu
$ ./t332 100
4.02976
3.75275
5.32832
8.53189
8.46641
$

Note that the functor in the above code divides by D+1. It would be trivial to change it to divide by D, instead, to match your code (but not your stated intent).

How to do weighted average of two device_vectors with points selected using map in thrust?

1 Answers1