Reduce Boost dependency footprint

stan/math/prim/core/init_threadpool_tbb.hpp

@@ -3,7 +3,7 @@
 
 #include <stan/math/prim/err/invalid_argument.hpp>
 
-#include <boost/lexical_cast.hpp>
+#include <charconv>
 
 #ifndef TBB_INTERFACE_NEW
 #include <tbb/tbb_stddef.h>
@@ -21,6 +21,7 @@
 #endif
 
 #include <cstdlib>
+#include <cstring>
 #include <thread>
 
 namespace stan {
@@ -48,9 +49,10 @@ inline int get_num_threads() {
 #ifdef STAN_THREADS
   const char* env_stan_num_threads = std::getenv("STAN_NUM_THREADS");
   if (env_stan_num_threads != nullptr) {
-    try {
-      const int env_num_threads
-          = boost::lexical_cast<int>(env_stan_num_threads);
+    int env_num_threads = 0;
+    const char* end = env_stan_num_threads + std::strlen(env_stan_num_threads);
+    auto result = std::from_chars(env_stan_num_threads, end, env_num_threads);
+    if (result.ec == std::errc() && result.ptr == end) {
       if (env_num_threads > 0) {
         num_threads = env_num_threads;
       } else if (env_num_threads == -1) {
@@ -61,7 +63,7 @@ inline int get_num_threads() {
                          "The STAN_NUM_THREADS environment variable is '",
                          "' but it must be positive or -1");
       }
-    } catch (const boost::bad_lexical_cast&) {
+    } else {
       invalid_argument("get_num_threads(int)", "STAN_NUM_THREADS",
                        env_stan_num_threads,
                        "The STAN_NUM_THREADS environment variable is '",

stan/math/prim/fun/grad_2F1.hpp

@@ -13,7 +13,6 @@
 #include <stan/math/prim/fun/sign.hpp>
 #include <stan/math/prim/fun/hypergeometric_2F1.hpp>
 #include <cmath>
-#include <boost/optional.hpp>
 
 namespace stan {
 namespace math {

stan/math/prim/fun/hypergeometric_2F1.hpp

@@ -17,7 +17,7 @@
 #include <stan/math/prim/fun/sqrt.hpp>
 #include <stan/math/prim/fun/square.hpp>
 #include <stan/math/prim/fun/hypergeometric_pFq.hpp>
-#include <boost/optional.hpp>
+#include <optional>
 
 namespace stan {
 namespace math {
@@ -43,7 +43,7 @@ namespace internal {
  * @return Gauss hypergeometric function
  */
 template <typename Ta1, typename Ta2, typename Tb, typename Tz,
-          typename RtnT = boost::optional<return_type_t<Ta1, Ta2, Tb, Tz>>,
+          typename RtnT = std::optional<return_type_t<Ta1, Ta2, Tb, Tz>>,
           require_all_arithmetic_t<Ta1, Ta2, Tb, Tz>* = nullptr>
 inline RtnT hyper_2F1_special_cases(const Ta1& a1, const Ta2& a2, const Tb& b,
                                     const Tz& z) {
@@ -149,7 +149,7 @@ inline RtnT hyper_2F1_special_cases(const Ta1& a1, const Ta2& a2, const Tb& b,
  */
 template <typename Ta1, typename Ta2, typename Tb, typename Tz,
           typename ScalarT = return_type_t<Ta1, Ta2, Tb, Tz>,
-          typename OptT = boost::optional<ScalarT>,
+          typename OptT = std::optional<ScalarT>,
           require_all_arithmetic_t<Ta1, Ta2, Tb, Tz>* = nullptr>
 inline return_type_t<Ta1, Ta2, Tb, Tz> hypergeometric_2F1(const Ta1& a1,
                                                           const Ta2& a2,
@@ -168,15 +168,15 @@ inline return_type_t<Ta1, Ta2, Tb, Tz> hypergeometric_2F1(const Ta1& a1,
   // Check whether value can be calculated by any special-case rules
   // before estimating infinite sum
   OptT special_case_a1a2 = internal::hyper_2F1_special_cases(a1, a2, b, z);
-  if (special_case_a1a2.is_initialized()) {
-    return special_case_a1a2.get();
+  if (special_case_a1a2.has_value()) {
+    return special_case_a1a2.value();
   }
 
   // Check whether any special case rules apply with 'a' arguments reversed
   // as 2F1(a1, a2, b, z) = 2F1(a2, a1, b, z)
   OptT special_case_a2a1 = internal::hyper_2F1_special_cases(a2, a1, b, z);
-  if (special_case_a2a1.is_initialized()) {
-    return special_case_a2a1.get();
+  if (special_case_a2a1.has_value()) {
+    return special_case_a2a1.value();
   }
 
   Eigen::Matrix<double, 2, 1> a_args(2);

stan/math/prim/prob/bernoulli_logit_glm_rng.hpp

@@ -8,8 +8,7 @@
 #include <stan/math/prim/fun/as_array_or_scalar.hpp>
 #include <stan/math/prim/fun/inv_logit.hpp>
 #include <stan/math/prim/fun/scalar_seq_view.hpp>
-#include <boost/random/bernoulli_distribution.hpp>
-#include <boost/random/variate_generator.hpp>
+#include <random>
 #include <vector>
 
 namespace stan {
@@ -42,8 +41,6 @@ namespace math {
 template <typename T_x, typename T_alpha, typename T_beta, class RNG>
 inline typename VectorBuilder<true, int, T_alpha>::type bernoulli_logit_glm_rng(
     const T_x& x, const T_alpha& alpha, const T_beta& beta, RNG& rng) {
-  using boost::bernoulli_distribution;
-  using boost::variate_generator;
   using T_x_ref = ref_type_t<T_x>;
   using T_alpha_ref = ref_type_t<T_alpha>;
   using T_beta_ref = ref_type_t<T_beta>;
@@ -76,9 +73,8 @@ inline typename VectorBuilder<true, int, T_alpha>::type bernoulli_logit_glm_rng(
 
   for (size_t m = 0; m < M; ++m) {
     double theta_m = alpha_vec[m] + x_beta(m);
-    variate_generator<RNG&, bernoulli_distribution<>> bernoulli_rng(
-        rng, bernoulli_distribution<>(inv_logit(theta_m)));
-    output[m] = bernoulli_rng();
+    std::bernoulli_distribution bernoulli_rng(inv_logit(theta_m));
+    output[m] = bernoulli_rng(rng);
   }
 
   return output.data();

stan/math/prim/prob/bernoulli_logit_rng.hpp

@@ -6,8 +6,7 @@
 #include <stan/math/prim/fun/inv_logit.hpp>
 #include <stan/math/prim/fun/scalar_seq_view.hpp>
 #include <stan/math/prim/fun/size.hpp>
-#include <boost/random/bernoulli_distribution.hpp>
-#include <boost/random/variate_generator.hpp>
+#include <random>
 
 namespace stan {
 namespace math {
@@ -29,8 +28,6 @@ namespace math {
 template <typename T_t, class RNG>
 inline typename VectorBuilder<true, int, T_t>::type bernoulli_logit_rng(
     const T_t& t, RNG& rng) {
-  using boost::bernoulli_distribution;
-  using boost::variate_generator;
   ref_type_t<T_t> t_ref = t;
   check_finite("bernoulli_logit_rng", "Logit transformed probability parameter",
                t_ref);
@@ -40,9 +37,8 @@ inline typename VectorBuilder<true, int, T_t>::type bernoulli_logit_rng(
   VectorBuilder<true, int, T_t> output(N);
 
   for (size_t n = 0; n < N; ++n) {
-    variate_generator<RNG&, bernoulli_distribution<> > bernoulli_rng(
-        rng, bernoulli_distribution<>(inv_logit(t_vec[n])));
-    output[n] = bernoulli_rng();
+    std::bernoulli_distribution bernoulli_rng(inv_logit(t_vec[n]));
+    output[n] = bernoulli_rng(rng);
   }
 
   return output.data();

stan/math/prim/prob/bernoulli_rng.hpp

@@ -5,8 +5,7 @@
 #include <stan/math/prim/err.hpp>
 #include <stan/math/prim/fun/scalar_seq_view.hpp>
 #include <stan/math/prim/fun/size.hpp>
-#include <boost/random/bernoulli_distribution.hpp>
-#include <boost/random/variate_generator.hpp>
+#include <random>
 
 namespace stan {
 namespace math {
@@ -28,8 +27,6 @@ namespace math {
 template <typename T_theta, class RNG>
 inline typename VectorBuilder<true, int, T_theta>::type bernoulli_rng(
     const T_theta& theta, RNG& rng) {
-  using boost::bernoulli_distribution;
-  using boost::variate_generator;
   static constexpr const char* function = "bernoulli_rng";
   ref_type_t<T_theta> theta_ref = theta;
   check_bounded(function, "Probability parameter", value_of(theta_ref), 0.0,
@@ -40,9 +37,8 @@ inline typename VectorBuilder<true, int, T_theta>::type bernoulli_rng(
   VectorBuilder<true, int, T_theta> output(N);
 
   for (size_t n = 0; n < N; ++n) {
-    variate_generator<RNG&, bernoulli_distribution<> > bernoulli_rng(
-        rng, bernoulli_distribution<>(theta_vec[n]));
-    output[n] = bernoulli_rng();
+    std::bernoulli_distribution bernoulli_rng(theta_vec[n]);
+    output[n] = bernoulli_rng(rng);
   }
 
   return output.data();

stan/math/prim/prob/beta_proportion_rng.hpp

@@ -6,7 +6,6 @@
 #include <stan/math/prim/fun/max_size.hpp>
 #include <stan/math/prim/fun/scalar_seq_view.hpp>
 #include <stan/math/prim/fun/to_ref.hpp>
-#include <boost/random/variate_generator.hpp>
 
 namespace stan {
 namespace math {

stan/math/prim/prob/beta_rng.hpp

@@ -6,10 +6,8 @@
 #include <stan/math/prim/fun/log_sum_exp.hpp>
 #include <stan/math/prim/fun/max_size.hpp>
 #include <stan/math/prim/fun/scalar_seq_view.hpp>
-#include <boost/random/gamma_distribution.hpp>
-#include <boost/random/uniform_real_distribution.hpp>
-#include <boost/random/variate_generator.hpp>
 #include <cmath>
+#include <random>
 
 namespace stan {
 namespace math {
@@ -35,9 +33,6 @@ namespace math {
 template <typename T_shape1, typename T_shape2, class RNG>
 inline typename VectorBuilder<true, double, T_shape1, T_shape2>::type beta_rng(
     const T_shape1 &alpha, const T_shape2 &beta, RNG &rng) {
-  using boost::variate_generator;
-  using boost::random::gamma_distribution;
-  using boost::random::uniform_real_distribution;
   using T_alpha_ref = ref_type_t<T_shape1>;
   using T_beta_ref = ref_type_t<T_shape2>;
   static constexpr const char *function = "beta_rng";
@@ -53,29 +48,24 @@ inline typename VectorBuilder<true, double, T_shape1, T_shape2>::type beta_rng(
   size_t N = max_size(alpha, beta);
   VectorBuilder<true, double, T_shape1, T_shape2> output(N);
 
-  variate_generator<RNG &, uniform_real_distribution<>> uniform_rng(
-      rng, uniform_real_distribution<>(0.0, 1.0));
+  std::uniform_real_distribution<> uniform_rng(0.0, 1.0);
   for (size_t n = 0; n < N; ++n) {
     // If alpha and beta are large, trust the usual ratio of gammas
     // method for generating beta random variables. If any parameter
     // is small, work in log space and use Marsaglia and Tsang's trick
     if (alpha_vec[n] > 1.0 && beta_vec[n] > 1.0) {
-      variate_generator<RNG &, gamma_distribution<>> rng_gamma_alpha(
-          rng, gamma_distribution<>(alpha_vec[n], 1.0));
-      variate_generator<RNG &, gamma_distribution<>> rng_gamma_beta(
-          rng, gamma_distribution<>(beta_vec[n], 1.0));
-      double a = rng_gamma_alpha();
-      double b = rng_gamma_beta();
+      std::gamma_distribution<> rng_gamma_alpha(alpha_vec[n], 1.0);
+      std::gamma_distribution<> rng_gamma_beta(beta_vec[n], 1.0);
+      double a = rng_gamma_alpha(rng);
+      double b = rng_gamma_beta(rng);
       output[n] = a / (a + b);
     } else {
-      variate_generator<RNG &, gamma_distribution<>> rng_gamma_alpha(
-          rng, gamma_distribution<>(alpha_vec[n] + 1, 1.0));
-      variate_generator<RNG &, gamma_distribution<>> rng_gamma_beta(
-          rng, gamma_distribution<>(beta_vec[n] + 1, 1.0));
-      double log_a = std::log(uniform_rng()) / alpha_vec[n]
-                     + std::log(rng_gamma_alpha());
-      double log_b
-          = std::log(uniform_rng()) / beta_vec[n] + std::log(rng_gamma_beta());
+      std::gamma_distribution<> rng_gamma_alpha(alpha_vec[n] + 1, 1.0);
+      std::gamma_distribution<> rng_gamma_beta(beta_vec[n] + 1, 1.0);
+      double log_a = std::log(uniform_rng(rng)) / alpha_vec[n]
+                     + std::log(rng_gamma_alpha(rng));
+      double log_b = std::log(uniform_rng(rng)) / beta_vec[n]
+                     + std::log(rng_gamma_beta(rng));
       double log_sum = log_sum_exp(log_a, log_b);
       output[n] = std::exp(log_a - log_sum);
     }

stan/math/prim/prob/binomial_rng.hpp

@@ -5,8 +5,7 @@
 #include <stan/math/prim/err.hpp>
 #include <stan/math/prim/fun/max_size.hpp>
 #include <stan/math/prim/fun/scalar_seq_view.hpp>
-#include <boost/random/binomial_distribution.hpp>
-#include <boost/random/variate_generator.hpp>
+#include <random>
 
 namespace stan {
 namespace math {
@@ -31,8 +30,6 @@ namespace math {
 template <typename T_N, typename T_theta, class RNG>
 inline typename VectorBuilder<true, int, T_N, T_theta>::type binomial_rng(
     const T_N& N, const T_theta& theta, RNG& rng) {
-  using boost::binomial_distribution;
-  using boost::variate_generator;
   using T_N_ref = ref_type_t<T_N>;
   using T_theta_ref = ref_type_t<T_theta>;
   static constexpr const char* function = "binomial_rng";
@@ -50,10 +47,9 @@ inline typename VectorBuilder<true, int, T_N, T_theta>::type binomial_rng(
   VectorBuilder<true, int, T_N, T_theta> output(M);
 
   for (size_t m = 0; m < M; ++m) {
-    variate_generator<RNG&, binomial_distribution<> > binomial_rng(
-        rng, binomial_distribution<>(N_vec[m], theta_vec[m]));
+    std::binomial_distribution<> binomial_rng(N_vec[m], theta_vec[m]);
 
-    output[m] = binomial_rng();
+    output[m] = binomial_rng(rng);
   }
 
   return output.data();

stan/math/prim/prob/categorical_logit_rng.hpp

@@ -6,8 +6,7 @@
 #include <stan/math/prim/fun/cumulative_sum.hpp>
 #include <stan/math/prim/fun/softmax.hpp>
 #include <stan/math/prim/fun/Eigen.hpp>
-#include <boost/random/uniform_01.hpp>
-#include <boost/random/variate_generator.hpp>
+#include <random>
 
 namespace stan {
 namespace math {
@@ -27,16 +26,14 @@ namespace math {
  */
 template <class RNG>
 inline int categorical_logit_rng(const Eigen::VectorXd& beta, RNG& rng) {
-  using boost::uniform_01;
-  using boost::variate_generator;
   static constexpr const char* function = "categorical_logit_rng";
   check_finite(function, "Log odds parameter", beta);
 
-  variate_generator<RNG&, uniform_01<> > uniform01_rng(rng, uniform_01<>());
+  std::uniform_real_distribution<> uniform01_rng(0, 1);
   Eigen::VectorXd theta = softmax(beta);
   Eigen::VectorXd index = cumulative_sum(theta);
 
-  double c = uniform01_rng();
+  double c = uniform01_rng(rng);
   int b = 0;
   while (c > index(b)) {
     b++;

stan/math/prim/prob/categorical_rng.hpp

@@ -4,28 +4,25 @@
 #include <stan/math/prim/meta.hpp>
 #include <stan/math/prim/err.hpp>
 #include <stan/math/prim/fun/cumulative_sum.hpp>
-#include <boost/random/uniform_01.hpp>
-#include <boost/random/variate_generator.hpp>
+#include <random>
 
 namespace stan {
 namespace math {
 
 template <class RNG>
 inline int categorical_rng(
     const Eigen::Matrix<double, Eigen::Dynamic, 1>& theta, RNG& rng) {
-  using boost::uniform_01;
-  using boost::variate_generator;
   static constexpr const char* function = "categorical_rng";
   check_simplex(function, "Probabilities parameter", theta);
 
-  variate_generator<RNG&, uniform_01<> > uniform01_rng(rng, uniform_01<>());
+  std::uniform_real_distribution<> uniform01_rng(0, 1);
 
   Eigen::VectorXd index(theta.rows());
   index.setZero();
 
   index = cumulative_sum(theta);
 
-  double c = uniform01_rng();
+  double c = uniform01_rng(rng);
   int b = 0;
   while (c > index(b, 0)) {
     b++;

stan/math/prim/prob/cauchy_rng.hpp

@@ -5,8 +5,7 @@
 #include <stan/math/prim/err.hpp>
 #include <stan/math/prim/fun/max_size.hpp>
 #include <stan/math/prim/fun/scalar_seq_view.hpp>
-#include <boost/random/cauchy_distribution.hpp>
-#include <boost/random/variate_generator.hpp>
+#include <random>
 
 namespace stan {
 namespace math {
@@ -32,8 +31,6 @@ namespace math {
 template <typename T_loc, typename T_scale, class RNG>
 inline typename VectorBuilder<true, double, T_loc, T_scale>::type cauchy_rng(
     const T_loc& mu, const T_scale& sigma, RNG& rng) {
-  using boost::variate_generator;
-  using boost::random::cauchy_distribution;
   static constexpr const char* function = "cauchy_rng";
   using T_mu_ref = ref_type_t<T_loc>;
   using T_sigma_ref = ref_type_t<T_scale>;
@@ -50,9 +47,8 @@ inline typename VectorBuilder<true, double, T_loc, T_scale>::type cauchy_rng(
   VectorBuilder<true, double, T_loc, T_scale> output(N);
 
   for (size_t n = 0; n < N; ++n) {
-    variate_generator<RNG&, cauchy_distribution<> > cauchy_rng(
-        rng, cauchy_distribution<>(mu_vec[n], sigma_vec[n]));
-    output[n] = cauchy_rng();
+    std::cauchy_distribution<> cauchy_rng(mu_vec[n], sigma_vec[n]);
+    output[n] = cauchy_rng(rng);
   }
 
   return output.data();