Use CCCL random for device PRNG draws

cpp/src/decisiontree/batched-levelalgo/kernels/builder_kernels.cuh

@@ -13,9 +13,12 @@
 #include <cuml/common/utils.hpp>
 
 #include <raft/core/handle.hpp>
-#include <raft/random/rng.cuh>
+#include <raft/core/math.hpp>
 
 #include <cub/cub.cuh>
+#include <cuda/std/array>
+#include <cuda/std/limits>
+#include <cuda/std/random>
 
 namespace ML {
 namespace DT {
@@ -170,13 +173,11 @@ CUML_KERNEL void excess_sample_with_replacement_kernel(
   subsequence = fnv1a32(subsequence, uint32_t(treeid));
   subsequence = fnv1a32(subsequence, uint32_t(nodeid));
 
-  raft::random::PCGenerator gen(seed, subsequence, uint64_t(0));
-  raft::random::UniformIntDistParams<IdxT, uint64_t> uniform_int_dist_params;
-
-  uniform_int_dist_params.start = 0;
-  uniform_int_dist_params.end   = n;
-  uniform_int_dist_params.diff =
-    uint64_t(uniform_int_dist_params.end - uniform_int_dist_params.start);
+  cuda::std::philox4x64 rng(static_cast<cuda::std::philox4x64::result_type>(seed));
+  rng.set_counter(
+    cuda::std::array<cuda::std::philox4x64::result_type, cuda::std::philox4x64::word_count>{
+      0, 0, static_cast<cuda::std::philox4x64::result_type>(subsequence), 0});
+  cuda::std::uniform_int_distribution<IdxT> col_dist(IdxT{0}, static_cast<IdxT>(n - 1));
 
   IdxT n_uniques = 0;
   IdxT items[MAX_SAMPLES_PER_THREAD];
@@ -199,8 +200,7 @@ CUML_KERNEL void excess_sample_with_replacement_kernel(
       // but this ensures some forward progress in order to generate at least 'k' unique random
       // samples.
       if (mask[thread_local_sample_idx] == 0 and cta_sample_idx < n_parallel_samples)
-        raft::random::custom_next(
-          gen, &items[thread_local_sample_idx], uniform_int_dist_params, IdxT(0), IdxT(0));
+        items[thread_local_sample_idx] = col_dist(rng);
       else if (mask[thread_local_sample_idx] ==
                0)  // indices that exceed `n_parallel_samples` will not generate
         items[thread_local_sample_idx] = n;  // n is outside valid range [0, n-1]
@@ -297,17 +297,20 @@ CUML_KERNEL void algo_L_sample_kernel(int* colids,
   if (tid >= work_items_size) return;
   const uint32_t nodeid = work_items[tid].idx;
   uint64_t subsequence  = (uint64_t(treeid) << 32) | uint64_t(nodeid);
-  raft::random::PCGenerator gen(seed, subsequence, uint64_t(0));
-  raft::random::UniformIntDistParams<IdxT, uint64_t> uniform_int_dist_params;
-  uniform_int_dist_params.start = 0;
-  uniform_int_dist_params.end   = k;
-  uniform_int_dist_params.diff =
-    uint64_t(uniform_int_dist_params.end - uniform_int_dist_params.start);
+  cuda::std::philox4x64 rng(static_cast<cuda::std::philox4x64::result_type>(seed));
+  rng.set_counter(
+    cuda::std::array<cuda::std::philox4x64::result_type, cuda::std::philox4x64::word_count>{
+      0, 0, static_cast<cuda::std::philox4x64::result_type>(subsequence), 0});
+  cuda::std::uniform_real_distribution<float> uniform01(cuda::std::numeric_limits<float>::min(),
+                                                        1.0f);
+  cuda::std::uniform_int_distribution<IdxT> sample_dist(IdxT{0}, static_cast<IdxT>(k - 1));
+  constexpr double max_W = 1.0 - static_cast<double>(cuda::std::numeric_limits<float>::epsilon());
   float fp_uniform_val;
   IdxT int_uniform_val;
   // fp_uniform_val will have a random value between 0 and 1
-  gen.next(fp_uniform_val);
-  double W = raft::exp(raft::log(fp_uniform_val) / k);
+  fp_uniform_val = uniform01(rng);
+  double W       = raft::exp(raft::log(fp_uniform_val) / k);
+  W              = W < max_W ? W : max_W;
 
   size_t col(0);
   // initially fill the reservoir array in increasing order of cols till k
@@ -321,15 +324,16 @@ CUML_KERNEL void algo_L_sample_kernel(int* colids,
   // randomly sample from a geometric distribution
   while (col < n) {
     // fp_uniform_val will have a random value between 0 and 1
-    gen.next(fp_uniform_val);
+    fp_uniform_val = uniform01(rng);
     col += static_cast<int>(raft::log(fp_uniform_val) / raft::log(1 - W)) + 1;
     if (col < n) {
       // int_uniform_val will now have a random value between 0...k
-      raft::random::custom_next(gen, &int_uniform_val, uniform_int_dist_params, IdxT(0), IdxT(0));
+      int_uniform_val                   = sample_dist(rng);
       colids[tid * k + int_uniform_val] = col;  // the bad memory coalescing here is hidden
       // fp_uniform_val will have a random value between 0 and 1
-      gen.next(fp_uniform_val);
+      fp_uniform_val = uniform01(rng);
       W *= raft::exp(raft::log(fp_uniform_val) / k);
+      W = W < max_W ? W : max_W;
     }
   }
 }
@@ -348,13 +352,16 @@ CUML_KERNEL void adaptive_sample_kernel(int* colids,
   const uint32_t nodeid = work_items[tid].idx;
 
   uint64_t subsequence = (uint64_t(treeid) << 32) | uint64_t(nodeid);
-  raft::random::PCGenerator gen(seed, subsequence, uint64_t(0));
+  cuda::std::philox4x64 rng(static_cast<cuda::std::philox4x64::result_type>(seed));
+  rng.set_counter(
+    cuda::std::array<cuda::std::philox4x64::result_type, cuda::std::philox4x64::word_count>{
+      0, 0, static_cast<cuda::std::philox4x64::result_type>(subsequence), 0});
+  cuda::std::uniform_int_distribution<uint32_t> toss_dist;
 
   int selected_count = 0;
   for (int i = 0; i < N; i++) {
-    uint32_t toss = 0;
-    gen.next(toss);
-    uint64_t lhs = uint64_t(M - selected_count);
+    uint32_t toss = toss_dist(rng);
+    uint64_t lhs  = uint64_t(M - selected_count);
     lhs <<= 32;
     uint64_t rhs = uint64_t(toss) * (N - i);
     if (lhs > rhs) {

cpp/src/decisiontree/batched-levelalgo/quantiles.cuh

@@ -13,14 +13,15 @@
 #include <raft/core/handle.hpp>
 #include <raft/core/nvtx.hpp>
 #include <raft/core/resource/comms.hpp>
-#include <raft/random/rng_device.cuh>
 #include <raft/util/cuda_utils.cuh>
 
 #include <rmm/device_uvector.hpp>
 #include <rmm/exec_policy.hpp>
 
 #include <cub/cub.cuh>
 #include <cuda/std/algorithm>
+#include <cuda/std/array>
+#include <cuda/std/random>
 #include <thrust/iterator/counting_iterator.h>
 #include <thrust/iterator/transform_iterator.h>
 #include <thrust/scan.h>
@@ -53,13 +54,12 @@ static __global__ void sampleOwnedColumnsKernel(T* out,
 
   std::uint64_t global_row = sample_idx;
   if (static_cast<std::uint64_t>(sample_count) != global_rows) {
-    raft::random::UniformIntDistParams<std::uint64_t, std::uint64_t> uniform_int_dist_params;
-    uniform_int_dist_params.start = 0;
-    uniform_int_dist_params.end   = global_rows;
-    uniform_int_dist_params.diff  = global_rows;
-    raft::random::PCGenerator gen(seed, static_cast<uint64_t>(sample_idx), uint64_t(0));
-    raft::random::custom_next(
-      gen, &global_row, uniform_int_dist_params, std::uint64_t(0), std::uint64_t(0));
+    cuda::std::philox4x64 rng(static_cast<cuda::std::philox4x64::result_type>(seed));
+    rng.set_counter(
+      cuda::std::array<cuda::std::philox4x64::result_type, cuda::std::philox4x64::word_count>{
+        0, 0, static_cast<cuda::std::philox4x64::result_type>(sample_idx), 0});
+    cuda::std::uniform_int_distribution<std::uint64_t> row_dist(std::uint64_t{0}, global_rows - 1);
+    global_row = row_dist(rng);
   }
 
   auto sample_end = ::cuda::std::lower_bound(

cpp/src/explainer/kernel_shap.cu

@@ -9,8 +9,8 @@
 #include <raft/core/handle.hpp>
 #include <raft/util/cudart_utils.hpp>
 
-#include <curand.h>
-#include <curand_kernel.h>
+#include <cuda/std/array>
+#include <cuda/std/random>
 
 namespace ML {
 namespace Explainer {
@@ -110,14 +110,17 @@ CUML_KERNEL void sampled_rows_kernel(IdxT* nsamples,
 
   // First k threads of block generate samples
   if (threadIdx.x < k_blk) {
-    curandStatePhilox4_32_10_t state;
-    curand_init((unsigned long long)seed, (unsigned long long)tid, 0, &state);
-    int rand_idx = (int)(curand_uniform(&state) * ncols);
+    cuda::std::philox4x64 rng(static_cast<cuda::std::philox4x64::result_type>(seed));
+    rng.set_counter(
+      cuda::std::array<cuda::std::philox4x64::result_type, cuda::std::philox4x64::word_count>{
+        0, 0, static_cast<cuda::std::philox4x64::result_type>(tid), 0});
+    cuda::std::uniform_int_distribution<int> col_dist(0, static_cast<int>(ncols) - 1);
+    int rand_idx = col_dist(rng);
 
     // Since X is initialized to 0, we quickly check for collisions (if k_blk << ncols the
     // likelihood of collisions is low)
     while (atomicExch(&(X[2 * blockIdx.x * ncols + rand_idx]), 1) == 1) {
-      rand_idx = (int)(curand_uniform(&state) * ncols);
+      rand_idx = col_dist(rng);
     }
   }
   __syncthreads();

cpp/src/genetic/genetic.cu

@@ -22,6 +22,9 @@
 #include <rmm/mr/per_device_resource.hpp>
 #include <rmm/resource_ref.hpp>
 
+#include <cuda/std/array>
+#include <cuda/std/random>
+
 #include <device_launch_parameters.h>
 
 #include <algorithm>
@@ -58,19 +61,19 @@ CUML_KERNEL void batched_tournament_kernel(const program_t progs,
   int idx = blockIdx.x * blockDim.x + threadIdx.x;
   if (idx >= n_tours) return;
 
-  raft::random::detail::PhiloxGenerator rng(seeds[idx], idx, 0);
-
-  int r;
-  rng.next(r);
+  cuda::std::philox4x64 rng(static_cast<cuda::std::philox4x64::result_type>(seeds[idx]));
+  rng.set_counter(
+    cuda::std::array<cuda::std::philox4x64::result_type, cuda::std::philox4x64::word_count>{
+      0, 0, static_cast<cuda::std::philox4x64::result_type>(idx), 0});
+  cuda::std::uniform_int_distribution<int> program_dist(0, n_progs - 1);
 
   // Define optima values
-  int opt           = r % n_progs;
+  int opt           = program_dist(rng);
   float opt_penalty = parsimony * progs[opt].len * (2 * criterion - 1);
   float opt_score   = progs[opt].raw_fitness_ - opt_penalty;
 
   for (int s = 1; s < tour_size; ++s) {
-    rng.next(r);
-    int curr           = r % n_progs;
+    int curr           = program_dist(rng);
     float curr_penalty = parsimony * progs[curr].len * (2 * criterion - 1);
     float curr_score   = progs[curr].raw_fitness_ - curr_penalty;
 

cpp/src/umap/simpl_set_embed/algo.cuh

@@ -33,7 +33,6 @@
 #include <thrust/sort.h>
 #include <thrust/system/cuda/execution_policy.h>
 
-#include <curand.h>
 #include <math.h>
 #include <stdint.h>
 

cpp/src/umap/simpl_set_embed/optimize_batch_kernel.cuh

@@ -9,10 +9,11 @@
 
 #include <cuml/manifold/umapparams.h>
 
-#include <raft/random/rng.cuh>
 #include <raft/util/cuda_utils.cuh>
 #include <raft/util/cudart_utils.hpp>
 
+#include <cuda/std/array>
+#include <cuda/std/random>
 #include <cuda_runtime_api.h>
 
 #include <stdint.h>
@@ -263,11 +264,13 @@ CUML_KERNEL void optimize_batch_kernel_reg(T const* head_embedding,
     /**
      * Negative sampling stage
      */
-    raft::random::detail::PhiloxGenerator gen((uint64_t)seed, (nnz_t)row, 0);
+    cuda::std::philox4x64 rng(static_cast<cuda::std::philox4x64::result_type>(seed));
+    rng.set_counter(
+      cuda::std::array<cuda::std::philox4x64::result_type, cuda::std::philox4x64::word_count>{
+        0, 0, static_cast<cuda::std::philox4x64::result_type>(row), 0});
+    cuda::std::uniform_int_distribution<int> tail_dist(0, tail_n.d - 1);
     for (int p = 0; p < n_neg_samples; p++) {
-      int r;
-      gen.next(r);
-      nnz_t t                  = r % tail_n;
+      nnz_t t                  = tail_dist(rng);
       T const* negative_sample = tail_embedding + (t * n_components);
       T negative_sample_reg[n_components];
       for (int i = 0; i < n_components; ++i) {
@@ -428,11 +431,13 @@ CUML_KERNEL void optimize_batch_kernel(T const* head_embedding,
     /**
      * Negative sampling stage
      */
-    raft::random::detail::PhiloxGenerator gen((uint64_t)seed, (nnz_t)row, 0);
+    cuda::std::philox4x64 rng(static_cast<cuda::std::philox4x64::result_type>(seed));
+    rng.set_counter(
+      cuda::std::array<cuda::std::philox4x64::result_type, cuda::std::philox4x64::word_count>{
+        0, 0, static_cast<cuda::std::philox4x64::result_type>(row), 0});
+    cuda::std::uniform_int_distribution<int> tail_dist(0, tail_n.d - 1);
     for (int p = 0; p < n_neg_samples; p++) {
-      int r;
-      gen.next(r);
-      nnz_t t                  = r % tail_n;
+      nnz_t t                  = tail_dist(rng);
       T const* negative_sample = tail_embedding + (t * n_components);
       if constexpr (use_shared_mem) {
         dist_squared = rdist<T>(current_buffer, negative_sample, n_components);
@@ -582,11 +587,13 @@ CUML_KERNEL void optimize_sequential_kernel_vertex_per_thread(T const* head_embe
       int n_neg_samples =
         static_cast<int>(T(epoch - _epoch_of_next_negative_sample) / epochs_per_negative_sample);
 
-      raft::random::detail::PhiloxGenerator gen(seed, static_cast<nnz_t>(e), 0);
+      cuda::std::philox4x64 rng(static_cast<cuda::std::philox4x64::result_type>(seed));
+      rng.set_counter(
+        cuda::std::array<cuda::std::philox4x64::result_type, cuda::std::philox4x64::word_count>{
+          0, 0, static_cast<cuda::std::philox4x64::result_type>(e), 0});
+      cuda::std::uniform_int_distribution<int> tail_dist(0, tail_n.d - 1);
       for (int p = 0; p < n_neg_samples; p++) {
-        int r;
-        gen.next(r);
-        nnz_t t                  = r % tail_n;
+        nnz_t t                  = tail_dist(rng);
         T const* negative_sample = tail_embedding + (t * N_COMPONENTS);
 
 #pragma unroll
@@ -746,11 +753,13 @@ CUML_KERNEL void optimize_sequential_kernel_vertex_per_warp(T const* head_embedd
       int n_neg_samples =
         static_cast<int>(T(epoch - _epoch_of_next_negative_sample) / epochs_per_negative_sample);
 
-      raft::random::detail::PhiloxGenerator gen(seed, static_cast<nnz_t>(e), 0);
+      cuda::std::philox4x64 rng(static_cast<cuda::std::philox4x64::result_type>(seed));
+      rng.set_counter(
+        cuda::std::array<cuda::std::philox4x64::result_type, cuda::std::philox4x64::word_count>{
+          0, 0, static_cast<cuda::std::philox4x64::result_type>(e), 0});
+      cuda::std::uniform_int_distribution<int> tail_dist(0, tail_n.d - 1);
       for (int p = 0; p < n_neg_samples; p++) {
-        int r;
-        gen.next(r);
-        nnz_t t = r % tail_n;
+        nnz_t t = tail_dist(rng);
 
         partial_dist = T(0);
 #pragma unroll

cpp/src/umap/simpl_set_embed/optimize_inverse_kernel.cuh

@@ -10,12 +10,14 @@
 #include <cuml/manifold/umapparams.h>
 
 #include <raft/linalg/unary_op.cuh>
-#include <raft/random/rng.cuh>
 #include <raft/util/cuda_utils.cuh>
 #include <raft/util/cudart_utils.hpp>
 
 #include <rmm/device_uvector.hpp>
 
+#include <cuda/std/array>
+#include <cuda/std/random>
+
 namespace UMAPAlgo {
 namespace SimplSetEmbed {
 namespace Algo {
@@ -128,11 +130,13 @@ CUML_KERNEL void optimize_inverse_kernel(T* head_embedding,
     T next_neg = epoch_of_next_negative_sample[row];
     int n_neg  = int((epoch - next_neg) / epochs_per_negative_sample);
 
-    raft::random::detail::PhiloxGenerator gen(seed, row, 0);
+    cuda::std::philox4x64 rng(static_cast<cuda::std::philox4x64::result_type>(seed));
+    rng.set_counter(
+      cuda::std::array<cuda::std::philox4x64::result_type, cuda::std::philox4x64::word_count>{
+        0, 0, static_cast<cuda::std::philox4x64::result_type>(row), 0});
+    cuda::std::uniform_int_distribution<nnz_t> vertex_dist(0, n_vertices - 1);
     for (int p = 0; p < n_neg; ++p) {
-      int r;
-      gen.next(r);
-      nnz_t t           = abs(r) % n_vertices;
+      nnz_t t           = vertex_dist(rng);
       const T* negative = tail_embedding + t * n_components;
 
       // Compute distance to negative sample

python/cuml/tests/test_random_forest.py

@@ -285,7 +285,7 @@ def test_tweedie_convergence(max_depth, split_criterion):
 
     # model trained on tweedie data with tweedie criterion must perform
     # no worse on tweedie loss, allowing tiny numeric variation.
-    deviance_tolerance = 1e-5
+    deviance_tolerance = 1e-4
     assert (
         mse_tweedie_deviance + deviance_tolerance >= tweedie_tweedie_deviance
     )