Replace boost/random with stdlib random

src/stan/io/random_var_context.hpp

@@ -3,9 +3,9 @@
 
 #include <stan/io/var_context.hpp>
 #include <stan/io/validate_dims.hpp>
-#include <boost/random/uniform_real_distribution.hpp>
 #include <algorithm>
 #include <limits>
+#include <random>
 #include <string>
 #include <vector>
 
@@ -50,8 +50,7 @@ class random_var_context : public var_context {
       for (size_t n = 0; n < num_unconstrained_; ++n)
         unconstrained_params_[n] = 0.0;
     } else {
-      boost::random::uniform_real_distribution<double> unif(-init_radius,
-                                                            init_radius);
+      std::uniform_real_distribution<double> unif(-init_radius, init_radius);
       for (size_t n = 0; n < num_unconstrained_; ++n)
         unconstrained_params_[n] = unif(rng);
     }

src/stan/mcmc/hmc/base_hmc.hpp

@@ -6,9 +6,9 @@
 #include <stan/callbacks/structured_writer.hpp>
 #include <stan/mcmc/base_mcmc.hpp>
 #include <stan/mcmc/hmc/hamiltonians/ps_point.hpp>
-#include <boost/random/uniform_01.hpp>
 #include <cmath>
 #include <limits>
+#include <random>
 #include <stdexcept>
 #include <string>
 #include <vector>
@@ -35,7 +35,7 @@ class base_hmc : public base_mcmc {
         integrator_(),
         hamiltonian_(model),
         rand_int_(rng),
-        rand_uniform_(rand_int_),
+        rand_uniform_(),
         nom_epsilon_(0.1),
         epsilon_(nom_epsilon_),
         epsilon_jitter_(0.0) {}
@@ -196,7 +196,9 @@ class base_hmc : public base_mcmc {
     this->epsilon_ = this->nom_epsilon_;
     if (this->epsilon_jitter_)
       this->epsilon_
-          *= 1.0 + this->epsilon_jitter_ * (2.0 * this->rand_uniform_() - 1.0);
+          *= 1.0
+             + this->epsilon_jitter_
+                   * (2.0 * this->rand_uniform_(this->rand_int_) - 1.0);
   }
 
  protected:
@@ -207,7 +209,7 @@ class base_hmc : public base_mcmc {
   BaseRNG& rand_int_;
 
   // Uniform(0, 1) RNG
-  boost::uniform_01<BaseRNG&> rand_uniform_;
+  std::uniform_real_distribution<> rand_uniform_;
 
   double nom_epsilon_;
   double epsilon_;

src/stan/mcmc/hmc/hamiltonians/dense_e_metric.hpp

@@ -5,8 +5,7 @@
 #include <stan/math/prim.hpp>
 #include <stan/mcmc/hmc/hamiltonians/base_hamiltonian.hpp>
 #include <stan/mcmc/hmc/hamiltonians/dense_e_point.hpp>
-#include <boost/random/variate_generator.hpp>
-#include <boost/random/normal_distribution.hpp>
+#include <random>
 
 namespace stan {
 namespace mcmc {
@@ -42,13 +41,12 @@ class dense_e_metric : public base_hamiltonian<Model, dense_e_point, BaseRNG> {
 
   void sample_p(dense_e_point& z, BaseRNG& rng) {
     typedef typename stan::math::index_type<Eigen::VectorXd>::type idx_t;
-    boost::variate_generator<BaseRNG&, boost::normal_distribution<> >
-        rand_dense_gaus(rng, boost::normal_distribution<>());
+    std::normal_distribution<> rand_dense_gaus;
 
     Eigen::VectorXd u(z.p.size());
 
     for (idx_t i = 0; i < u.size(); ++i)
-      u(i) = rand_dense_gaus();
+      u(i) = rand_dense_gaus(rng);
 
     z.p = z.inv_e_metric_.llt().matrixU().solve(u);
   }

src/stan/mcmc/hmc/hamiltonians/diag_e_metric.hpp

@@ -4,8 +4,7 @@
 #include <stan/callbacks/logger.hpp>
 #include <stan/mcmc/hmc/hamiltonians/base_hamiltonian.hpp>
 #include <stan/mcmc/hmc/hamiltonians/diag_e_point.hpp>
-#include <boost/random/variate_generator.hpp>
-#include <boost/random/normal_distribution.hpp>
+#include <random>
 
 namespace stan {
 namespace mcmc {
@@ -42,11 +41,10 @@ class diag_e_metric : public base_hamiltonian<Model, diag_e_point, BaseRNG> {
   }
 
   void sample_p(diag_e_point& z, BaseRNG& rng) {
-    boost::variate_generator<BaseRNG&, boost::normal_distribution<> >
-        rand_diag_gaus(rng, boost::normal_distribution<>());
+    std::normal_distribution<> rand_diag_gaus;
 
     for (int i = 0; i < z.p.size(); ++i)
-      z.p(i) = rand_diag_gaus() / sqrt(z.inv_e_metric_(i));
+      z.p(i) = rand_diag_gaus(rng) / sqrt(z.inv_e_metric_(i));
   }
 };
 

src/stan/mcmc/hmc/hamiltonians/softabs_metric.hpp

@@ -4,8 +4,7 @@
 #include <stan/math/mix.hpp>
 #include <stan/mcmc/hmc/hamiltonians/base_hamiltonian.hpp>
 #include <stan/mcmc/hmc/hamiltonians/softabs_point.hpp>
-#include <boost/random/variate_generator.hpp>
-#include <boost/random/normal_distribution.hpp>
+#include <random>
 
 namespace stan {
 namespace mcmc {
@@ -83,13 +82,12 @@ class softabs_metric : public base_hamiltonian<Model, softabs_point, BaseRNG> {
   }
 
   void sample_p(softabs_point& z, BaseRNG& rng) {
-    boost::variate_generator<BaseRNG&, boost::normal_distribution<> >
-        rand_unit_gaus(rng, boost::normal_distribution<>());
+    std::normal_distribution<> rand_unit_gaus;
 
     Eigen::VectorXd a(z.p.size());
 
     for (idx_t n = 0; n < z.p.size(); ++n)
-      a(n) = sqrt(z.softabs_lambda(n)) * rand_unit_gaus();
+      a(n) = sqrt(z.softabs_lambda(n)) * rand_unit_gaus(rng);
 
     z.p = z.eigen_deco.eigenvectors() * a;
   }

src/stan/mcmc/hmc/hamiltonians/unit_e_metric.hpp

@@ -3,8 +3,7 @@
 
 #include <stan/mcmc/hmc/hamiltonians/base_hamiltonian.hpp>
 #include <stan/mcmc/hmc/hamiltonians/unit_e_point.hpp>
-#include <boost/random/variate_generator.hpp>
-#include <boost/random/normal_distribution.hpp>
+#include <random>
 
 namespace stan {
 namespace mcmc {
@@ -37,11 +36,10 @@ class unit_e_metric : public base_hamiltonian<Model, unit_e_point, BaseRNG> {
   }
 
   void sample_p(unit_e_point& z, BaseRNG& rng) {
-    boost::variate_generator<BaseRNG&, boost::normal_distribution<> >
-        rand_unit_gaus(rng, boost::normal_distribution<>());
+    std::normal_distribution<> rand_unit_gaus;
 
     for (int i = 0; i < z.p.size(); ++i)
-      z.p(i) = rand_unit_gaus();
+      z.p(i) = rand_unit_gaus(rng);
   }
 };
 

src/stan/mcmc/hmc/nuts/base_nuts.hpp

@@ -128,7 +128,7 @@ class base_nuts : public base_hmc<Model, Hamiltonian, Integrator, BaseRNG> {
       bool valid_subtree = false;
       double log_sum_weight_subtree = -std::numeric_limits<double>::infinity();
 
-      if (this->rand_uniform_() > 0.5) {
+      if (this->rand_uniform_(this->rand_int_) > 0.5) {
         // Extend the current trajectory forward
         this->z_.ps_point::operator=(z_fwd);
         rho_bck = rho;
@@ -164,7 +164,7 @@ class base_nuts : public base_hmc<Model, Hamiltonian, Integrator, BaseRNG> {
         z_sample = z_propose;
       } else {
         double accept_prob = std::exp(log_sum_weight_subtree - log_sum_weight);
-        if (this->rand_uniform_() < accept_prob)
+        if (this->rand_uniform_(this->rand_int_) < accept_prob)
           z_sample = z_propose;
       }
 
@@ -328,7 +328,7 @@ class base_nuts : public base_hmc<Model, Hamiltonian, Integrator, BaseRNG> {
     } else {
       double accept_prob
           = std::exp(log_sum_weight_final - log_sum_weight_subtree);
-      if (this->rand_uniform_() < accept_prob)
+      if (this->rand_uniform_(this->rand_int_) < accept_prob)
         z_propose = z_propose_final;
     }
 

src/stan/mcmc/hmc/nuts_classic/base_nuts_classic.hpp

@@ -84,7 +84,7 @@ class base_nuts_classic
     util.H0 = this->hamiltonian_.H(this->z_);
 
     // Sample the slice variable
-    util.log_u = std::log(this->rand_uniform_());
+    util.log_u = std::log(this->rand_uniform_(this->rand_int_));
 
     // Build a balanced binary tree until the NUTS criterion fails
     util.criterion = true;
@@ -101,7 +101,7 @@ class base_nuts_classic
       ps_point* z = 0;
       Eigen::VectorXd* rho = 0;
 
-      if (this->rand_uniform_() > 0.5) {
+      if (this->rand_uniform_(this->rand_int_) > 0.5) {
         z = &z_plus;
         rho = &rho_plus;
         util.sign = 1;
@@ -133,7 +133,7 @@ class base_nuts_classic
         subtree_prob = n_valid_subtree ? 1 : 0;
       }
 
-      if (this->rand_uniform_() < subtree_prob)
+      if (this->rand_uniform_(this->rand_int_) < subtree_prob)
         z_sample = z_propose;
 
       n_valid += n_valid_subtree;
@@ -227,7 +227,8 @@ class base_nuts_classic
       double accept_prob
           = static_cast<double>(n2) / static_cast<double>(n1 + n2);
 
-      if (util.criterion && (this->rand_uniform_() < accept_prob))
+      if (util.criterion
+          && (this->rand_uniform_(this->rand_int_) < accept_prob))
         z_propose = z_propose_right;
 
       Eigen::VectorXd& subtree_rho = left_subtree_rho;

src/stan/mcmc/hmc/static/base_static_hmc.hpp

@@ -59,7 +59,7 @@ class base_static_hmc
 
     double acceptProb = std::exp(H0 - h);
 
-    if (acceptProb < 1 && this->rand_uniform_() > acceptProb)
+    if (acceptProb < 1 && this->rand_uniform_(this->rand_int_) > acceptProb)
       this->z_.ps_point::operator=(z_init);
 
     acceptProb = acceptProb > 1 ? 1 : acceptProb;

src/stan/mcmc/hmc/static_uniform/base_static_uniform.hpp

@@ -4,9 +4,9 @@
 #include <stan/callbacks/logger.hpp>
 #include <stan/mcmc/hmc/base_hmc.hpp>
 #include <stan/mcmc/hmc/hamiltonians/ps_point.hpp>
-#include <boost/random/uniform_int_distribution.hpp>
 #include <cmath>
 #include <limits>
+#include <random>
 #include <string>
 #include <vector>
 
@@ -46,7 +46,7 @@ class base_static_uniform
     double sum_prob = 1;
     double sum_metro_prob = 1;
 
-    boost::random::uniform_int_distribution<> uniform(0, L_ - 1);
+    std::uniform_int_distribution<> uniform(0, L_ - 1);
     int Lp = uniform(this->rand_int_);
 
     for (int l = 0; l < Lp; ++l) {
@@ -61,7 +61,7 @@ class base_static_uniform
       sum_prob += prob;
       sum_metro_prob += prob > 1 ? 1 : prob;
 
-      if (this->rand_uniform_() < prob / sum_prob)
+      if (this->rand_uniform_(this->rand_int_) < prob / sum_prob)
         z_sample = this->z_;
     }
 
@@ -79,7 +79,7 @@ class base_static_uniform
       sum_prob += prob;
       sum_metro_prob += prob > 1 ? 1 : prob;
 
-      if (this->rand_uniform_() < prob / sum_prob)
+      if (this->rand_uniform_(this->rand_int_) < prob / sum_prob)
         z_sample = this->z_;
     }
 

src/stan/mcmc/hmc/xhmc/base_xhmc.hpp

@@ -82,7 +82,7 @@ class base_xhmc : public base_hmc<Model, Hamiltonian, Integrator, BaseRNG> {
       double ave_subtree = 0;
       double log_sum_weight_subtree = -std::numeric_limits<double>::infinity();
 
-      if (this->rand_uniform_() > 0.5) {
+      if (this->rand_uniform_(this->rand_int_) > 0.5) {
         this->z_.ps_point::operator=(z_plus);
         valid_subtree = build_tree(this->depth_, z_propose, ave_subtree,
                                    log_sum_weight_subtree, H0, 1, n_leapfrog,
@@ -105,7 +105,7 @@ class base_xhmc : public base_hmc<Model, Hamiltonian, Integrator, BaseRNG> {
       ++(this->depth_);
 
       double accept_prob = std::exp(log_sum_weight_subtree - log_sum_weight);
-      if (this->rand_uniform_() < accept_prob)
+      if (this->rand_uniform_(this->rand_int_) < accept_prob)
         z_sample = z_propose;
 
       // Break if exhaustion criterion is satisfied
@@ -224,7 +224,7 @@ class base_xhmc : public base_hmc<Model, Hamiltonian, Integrator, BaseRNG> {
 
     double accept_prob
         = std::exp(log_sum_weight_right - log_sum_weight_subtree);
-    if (this->rand_uniform_() < accept_prob)
+    if (this->rand_uniform_(this->rand_int_) < accept_prob)
       z_propose = z_propose_right;
 
     return std::abs(ave_subtree) >= x_delta_;

src/stan/services/pathfinder/multi.hpp

@@ -17,7 +17,7 @@
 #include <tbb/parallel_for.h>
 #include <tbb/concurrent_vector.h>
 #include <tbb/concurrent_queue.h>
-#include <boost/random/discrete_distribution.hpp>
+#include <random>
 #include <string>
 #include <vector>
 
@@ -230,14 +230,11 @@ inline int pathfinder_lbfgs_multi(
   Eigen::Array<double, Eigen::Dynamic, 1> weight_vals
       = stan::services::psis::psis_weights(lp_ratios, tail_len, logger);
   stan::rng_t rng = util::create_rng(random_seed, stride_id);
-  using discrete_dist_t
-      = boost::random::discrete_distribution<Eigen::Index, double>;
-  boost::variate_generator<stan::rng_t&, discrete_dist_t> rand_psis_idx(
-      rng, discrete_dist_t(boost::iterator_range<double*>(
-               weight_vals.data(), weight_vals.data() + weight_vals.size())));
+  std::discrete_distribution<Eigen::Index> rand_psis_idx(weight_vals.cbegin(),
+                                                         weight_vals.cend());
   Eigen::Matrix<Eigen::Index, -1, 1> psis_draw_idxs(num_multi_draws);
   for (size_t i = 0; i <= num_multi_draws - 1; ++i) {
-    psis_draw_idxs.coeffRef(i) = rand_psis_idx();
+    psis_draw_idxs.coeffRef(i) = rand_psis_idx(rng);
   }
   /**
    * The sort helps two main things

src/stan/services/pathfinder/single.hpp

@@ -16,6 +16,7 @@
 #include <tbb/parallel_for.h>
 #include <tbb/concurrent_queue.h>
 #include <tbb/task_group.h>
+#include <random>
 #include <string>
 #include <vector>
 #include <atomic>
@@ -213,12 +214,12 @@ inline elbo_est_t est_approx_draws(LPF&& lp_fun, RNG&& rng,
                                    size_t num_samples, const EigVec& alpha,
                                    const std::string& iter_msg, Logger&& logger,
                                    bool calculate_lp = true) {
-  boost::variate_generator<stan::rng_t&, boost::normal_distribution<>>
-      rand_unit_gaus(rng, boost::normal_distribution<>());
+  std::normal_distribution<> rand_unit_gaus;
   const auto num_params = taylor_approx.x_center.size();
   size_t lp_fun_calls = 0;
-  Eigen::MatrixXd unit_samps
-      = generate_matrix(rand_unit_gaus, num_params, num_samples);
+  Eigen::MatrixXd unit_samps = generate_matrix(
+      [&rng, &rand_unit_gaus]() { return rand_unit_gaus(rng); }, num_params,
+      num_samples);
   Eigen::Array<double, Eigen::Dynamic, 2> lp_mat(num_samples, 2);
   lp_mat.col(0) = (-taylor_approx.logdetcholHk)
                   + -0.5

src/test/unit/io/random_var_context_test.cpp

@@ -2,7 +2,6 @@
 #include <stan/io/empty_var_context.hpp>
 #include <stan/services/util/create_rng.hpp>
 #include <gtest/gtest.h>
-#include <boost/random/uniform_real_distribution.hpp>
 #include <test/test-models/good/services/test_lp.hpp>
 #include <test/unit/util.hpp>
 

src/test/unit/mcmc/hmc/nuts/softabs_nuts_test.cpp

@@ -338,15 +338,15 @@ TEST(McmcSoftAbsNuts, transition_test) {
 
   stan::mcmc::sample s = sampler.transition(init_sample, logger);
 
-  EXPECT_EQ(3, sampler.depth_);
-  EXPECT_EQ((2 << 3) - 1, sampler.n_leapfrog_);
+  EXPECT_EQ(5, sampler.depth_);
+  EXPECT_EQ((2 << 4) - 1, sampler.n_leapfrog_);
   EXPECT_FALSE(sampler.divergent_);
 
-  EXPECT_FLOAT_EQ(0.74693149, s.cont_params()(0));
-  EXPECT_FLOAT_EQ(-0.74414188, s.cont_params()(1));
-  EXPECT_FLOAT_EQ(0.60859376, s.cont_params()(2));
-  EXPECT_FLOAT_EQ(-0.74102008, s.log_prob());
-  EXPECT_FLOAT_EQ(0.99934167, s.accept_stat());
+  EXPECT_FLOAT_EQ(0.2057635, s.cont_params()(0));
+  EXPECT_FLOAT_EQ(0.87303215, s.cont_params()(1));
+  EXPECT_FLOAT_EQ(0.21624902, s.cont_params()(2));
+  EXPECT_FLOAT_EQ(-0.42564371, s.log_prob());
+  EXPECT_FLOAT_EQ(0.99925238, s.accept_stat());
   EXPECT_EQ("", debug.str());
   EXPECT_EQ("", info.str());
   EXPECT_EQ("", warn.str());