I want to compute mean and std using thrust and I found this code. I am trying to use complex values and I am having some problems.
Here is the code:
#include <thrust/device_vector.h>
#include <thrust/host_vector.h>
#include <thrust/transform_reduce.h>
#include <thrust/functional.h>
#include <thrust/extrema.h>
#include <cmath>
#include <float.h>
typedef struct
{
float re,im;
} mycomplex;
// structure used to accumulate the moments and other
// statistical properties encountered so far.
template <typename T>
struct summary_stats_data
{
T n;
T min;
T max;
T mean;
T M2;
// initialize to the identity element
void initialize()
{
n.re = mean.re = M2.re = 0;
n.im = mean.im = M2.im = 0;
min = std::numeric_limits<T>::max();
max = std::numeric_limits<T>::min();
}
float varianceRe() { return M2.re / ( n.re - 1 ); }
float varianceIm() { return M2.im / ( n.im - 1 ); }
float variance_nRe() { return M2.re / n.re; }
float variance_nIm() { return M2.im / n.im; }
};
// stats_unary_op is a functor that takes in a value x and
// returns a variace_data whose mean value is initialized to x.
template <typename T>
struct summary_stats_unary_op
{
__host__ __device__
summary_stats_data<T> operator()(const T& x) const
{
summary_stats_data<T> result;
result.n.re = 1;
result.n.im = 1;
result.min = x;
result.max = x;
result.mean = x;
result.M2.re = 0;
result.M2.im = 0;
return result;
}
};
// summary_stats_binary_op is a functor that accepts two summary_stats_data
// structs and returns a new summary_stats_data which are an
// approximation to the summary_stats for
// all values that have been agregated so far
template <typename T>
struct summary_stats_binary_op
: public thrust::binary_function<const summary_stats_data<T>&,
const summary_stats_data<T>&,
summary_stats_data<T> >
{
__host__ __device__
summary_stats_data<T> operator()(const summary_stats_data<T>& x, const summary_stats_data <T>& y) const
{
summary_stats_data<T> result;
// precompute some common subexpressions
T n;
n.re = x.n.re + y.n.re;
n.im = x.n.im + y.n.im;
T delta;
delta.re = y.mean.re - x.mean.re;
delta.im = y.mean.im - x.mean.im;
T delta2;
delta2.re = delta.re * delta.re;
delta2.im = delta.im * delta.im;
//Basic number of samples (n), min, and max
result.n = n;
result.min.re = thrust::min( x.min.re, y.min.re );
result.min.im = thrust::min( x.min.im, y.min.im );
result.max.re = thrust::max( x.max.re, y.max.re );
result.max.im = thrust::max( x.max.im, y.max.im );
result.mean.re = x.mean.re + delta.re * y.n.re / n.re;
result.mean.im = x.mean.im + delta.im * y.n.im / n.im;
result.M2.re = x.M2.re + y.M2.re;
result.M2.im = x.M2.im + y.M2.im;
result.M2.re += delta2.re * x.n.re * y.n.re / n.re;
result.M2.im += delta2.im * x.n.im * y.n.im / n.im;
return result;
}
};
template <typename Iterator>
void print_range(const std::string& name, Iterator first, Iterator last)
{
typedef typename std::iterator_traits<Iterator>::value_type T;
std::cout << name << ": ";
thrust::copy(first, last, std::ostream_iterator<T>(std::cout, " "));
std::cout << "\n";
}
int main(void)
{
typedef mycomplex T;
const int N = 4;
// initialize host array
thrust::host_vector<T> h_x(N);
h_x[ 0 ].re = h_x[ 0 ].im = 4.0f;
h_x[ 1 ].re = h_x[ 1 ].im = 7.0f;
h_x[ 2 ].re = h_x[ 2 ].im = 13.0f;
h_x[ 3 ].re = h_x[ 3 ].im = 16.0f;
// Copy host_vector H to device_vector D
thrust::device_vector<T> d_x = h_x;
// setup arguments
summary_stats_unary_op<T> unary_op;
summary_stats_binary_op<T> binary_op;
summary_stats_data<T> init;
init.initialize();
// compute summary statistics
summary_stats_data<T> result = thrust::transform_reduce( d_x.begin(), d_x.end(), unary_op, init, binary_op );
std::cout <<"******Summary Statistics Example*****"<<std::endl;
print_range("The data", d_x.begin(), d_x.end());
std::cout <<"Count : "<< result.n.re << std::endl;
std::cout <<"Minimum : "<< result.min.re <<std::endl;
std::cout <<"Maximum : "<< result.max.re <<std::endl;
std::cout <<"Mean : "<< result.mean.re << std::endl;
std::cout <<"Variance : "<< result.varianceRe() << std::endl;
std::cout <<"Standard Deviation : "<< std::sqrt(result.variance_nRe()) << std::endl;
return 0;
}
And I am receiving errors :
....include/c++/4.4.7/limits(284): error: no suitable constructor exists to convert from "int" to "mycomplex"
....include/c++/4.4.7/limits(282): error: no suitable constructor exists to convert from "int" to "mycomplex"
...include/c++/4.4.7/bits/stream_iterator.h(191): error: no operator "<<" matches these operands
The code was not originally designed to work with a complex value (i.e. an arbitrary struct). It was designed to work correctly for POD data types for which the following types of assignments make sense:
min = std::numeric_limits<T>::max();
max = std::numeric_limits<T>::min();
Since you have defined your template type to be mycomplex, which is apparently a struct definition you created, you cannot make a direct assignment from a scalar value to a struct. Perhaps a more correct description would be to say that std::numeric_limits<T>::max() doesn't know what entity to return when T = mycomplex.
Rather than define your own complex structure, from the standpoint of code re-use, I would suggest picking up the CUDA header cuComplex.h and using that. Using that framework with float-based complex data, you could do:
min = make_cuFloatComplex(std::numeric_limits<float>::max(), std::numeric_limits<float>::max());
Since the underlying code is thrust code, it might be more sensible to use thrust::complex<float> instead (whose functions are outlined in thrust/complex.h). This could potentially impact the rest of the code you have written/modified as well.
I see that you have converted most of the rest of the arithmetic to work on your particular struct definition, so possibly the only remaining thing would be to handle the << operator for cout. This could be done by overloading << for your data type, or breaking down the line that has the error into data types that cout understands (i.e. output the real and imaginary parts separately).