GPU-based sampling from an arbitrary distribition.

The sample_cdf() device function will draw random numbers from an arbitrary
disribution given a cumulative distribution function in the form of
a list of x,y points, beginning with y=0 and ending with y=1.  For an
example of how to convert a ROOT histogram to this form, see the unit
test in test_sample_cdf.py

1 files changed, 24 insertions, 0 deletions

diff --git a/src/random.h b/src/random.h
index b084c19..6abc415 100644
--- a/src/random.h
+++ b/src/random.h
@@ -19,6 +19,30 @@ __device__ float3 uniform_sphere(curandState *s)
 	return make_float3(c*cosf(theta), c*sinf(theta), u);
 }
 
+// Draw a random sample given a cumulative distribution function
+// Assumptions: ncdf >= 2, cdf_y[0] is 0.0, and cdf_y[ncdf-1] is 1.0
+__device__ float sample_cdf(curandState *rng, int ncdf, 
+			    float *cdf_x, float *cdf_y)
+{
+  float u = curand_uniform(rng);
+
+  // Find u in cdf_y with binary search
+  // list must contain at least 2 elements: 0.0 and 1.0
+  int lower=0;
+  int upper=ncdf-1;
+  while(lower < upper-1) {
+    int half = (lower+upper) / 2;
+    if (u < cdf_y[half])
+      upper = half;
+    else
+      lower = half;
+  }
+  
+  float frac = (u - cdf_y[lower]) / (cdf_y[upper] - cdf_y[lower]);
+  return cdf_x[lower] * frac + cdf_x[upper] * (1.0f - frac);
+}
+
+
 extern "C"
 {