Merge pull request #3412 from stan-dev/boost-random

Replace boost/random with stdlib random

Author: andrjohns

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_;
     }