Flirting With Models

#include <vector_functions.h>
 
#include <thrust/host_vector.h>
#include <thrust/device_vector.h>
 
#include <thrust/generate.h>
#include <thrust/transform.h>
#include <thrust/count.h>
 
#include <cstdlib>
#include <math.h>
 
float2 bottomLeft;
float2 topRight;
 
__device__
float f(float x) {
	return (x*x*__cosf(x));
}
 
struct above_below_function {
	__device__ 
	float operator()(float2 p) const {
		return (p.y < f(p.x)) ? 1 : 0;
	}
};
 
float2 make_random_float2(){
	return make_float2( (rand() / (float)RAND_MAX) * (float)(topRight.x - bottomLeft.x) + (float)bottomLeft.x,
						(rand() / (float)RAND_MAX) * (float)(topRight.y - bottomLeft.y) + (float)bottomLeft.y);
}
 
int main() {
	srand(time(0));
 
	bottomLeft.x = 0;
	bottomLeft.y = 0;
 
	topRight.x = 1;
	topRight.y = 1;
 
	int numPoints = 10000000;
 
	// find the area of our bounding square
	float area = (topRight.x - bottomLeft.x) * (topRight.y - bottomLeft.y);
 
	// generate our points
	thrust::host_vector<float2> h_locations(numPoints);		
	thrust::generate(h_locations.begin(), h_locations.end(), make_random_float2);
 
	// create two device vectors -- one for our locations, and another to create a count
	thrust::device_vector<float2> d_locations = h_locations;
	thrust::device_vector<unsigned int> d_boolean(numPoints);
 
	// for each element, if the y is below f(x), then put a 1 in the array.  otherwise, put a 0
	thrust::transform(d_locations.begin(), d_locations.end(), d_boolean.begin(), above_below_function());
 
	int total = thrust::count(d_boolean.begin(), d_boolean.end(), 1);
 
	std::cout << "Area under curve: " << area * (total / (float)numPoints) << std::endl;
 
	return 0;
}

5 Responses to “Monte-Carlo Integration using Thrust”

• Nathan Says:
  September 21st, 2009 at 12:26 am

  Neat example!

  Have you given any thought to how you might move the random number generation over to the GPU? One approach is to use something like the parallel Mersenne Twister (cf. code in the CUDA SDK) to generate the samples.

  However, I’m more interested whether hash functions are sufficient for most applications of RNGs. This Microsoft Research paper [1] argues that good hash functions make adequate random number generators (using trivial input to the hash). The idea seems sound, and the authors provide some empirical evidence to support the claim.

  Any thoughts?

  [1] http://research.microsoft.com/apps/pubs/default.aspx?id=70502
• Corey Says:
  September 21st, 2009 at 7:26 am

  Truthfully, I haven’t gotten nearly as much time to play with CUDA and Thrust as I would like. From the quick searching I did, it seemed like there were plenty of RNG on GPU examples out there, so I would imagine that it wouldn’t be difficult at all to get the random number generation to take place on the GPU.

  As for the hash question? Well, a good hash-function should return a hash value that is uniformly distributed over the range. In this sense, one could say that a hash function might make a good uniform RNG. Are there any particular implementations of this idea that you are thinking about?
• Nathan Says:
  September 21st, 2009 at 12:53 pm

  I was interested in whether there were cheap hash functions that were “good enough” for a broad class of randomized algorithms. The MD5 method used in the paper is good, but probably overkill for most applications.

  I’m thinking of relatively cheap hash functions like these:
  http://burtleburtle.net/bob/hash/integer.html
  http://www.concentric.net/~Ttwang/tech/inthash.htm

  Another technique would be to seed a small, serial per-thread RNG with hash(thread_id) and then generate a sequence within each thread.

  I guess I should just implement this scheme, run it through DIEHARD, and then see what pops out
• Corey Says:
  September 21st, 2009 at 2:00 pm

  Well, I won’t pretend that this is my area of expertise.

  Certainly, having the ability to create a uniform between [0,1] would be extremely powerful (with inverse CDFs giving you access to any other distribution you want), especially if it is fast/cheap. However, I guess what it comes down to is ‘how good is good enough’? If your problem domain does not require extreme accuracy in your PRNGs, then you could probably get away with cheap integer hashing. For something with higher requirements, I would need to see a proof that these cheap hash-functions truly create enough randomness (which the articles never really claim).

  One big issue I see right off the bat is that these integer hashes will spit out an integer as a result — not quite as useful if you want a value between [0,1]. Though, [http://www.concentric.net/~Ttwang/tech/inthash.htm] provides an interesting method of binning based on power of 2 table sizes. I suppose you could change your table size based on the accuracy you wanted your random values to be…

  Again, I think it may come down to the proof of randomness. Though, given how fast these functions are compared to more complex, more random alternatives, it may be worth it for applications that really only need PRNGs, and not RNGs…

  Keep me updated!
• Nathan Says:
  September 21st, 2009 at 5:25 pm

  Sure thing!

Leave a Reply