Merge commit '726f9914af03df2b743f4541a4d0dae243ef267b' into HEAD

Author: yashikno

Jenkinsfile

@@ -592,15 +592,14 @@ pipeline {
                     sh """#!/bin/bash
                         set -x
                         make doxygen
+                        cd doc/api/html
+                        git init
+                        git checkout -b gh-pages
                         git config user.email "mc.stanislaw@gmail.com"
                         git config user.name "Stan Jenkins"
-                        git checkout --detach
-                        git branch -D gh-pages
-                        git push https://${GIT_USERNAME}:${GIT_PASSWORD}@github.com/stan-dev/math.git :gh-pages
-                        git checkout --orphan gh-pages
-                        git add -f doc
+                        git add -f .
                         git commit --author='Stan BuildBot <mc.stanislaw@gmail.com>' -m "auto generated docs from Jenkins"
-                        git subtree push --prefix doc/api/html https://${GIT_USERNAME}:${GIT_PASSWORD}@github.com/stan-dev/math.git gh-pages
+                        git push https://${GIT_USERNAME}:${GIT_PASSWORD}@github.com/stan-dev/math.git gh-pages --force
                         """
                 }
             }

stan/math/fwd/fun.hpp

@@ -116,6 +116,7 @@
 #include <stan/math/fwd/fun/tcrossprod.hpp>
 #include <stan/math/fwd/fun/tgamma.hpp>
 #include <stan/math/fwd/fun/to_fvar.hpp>
+#include <stan/math/fwd/fun/trace_dot.hpp>
 #include <stan/math/fwd/fun/trace_quad_form.hpp>
 #include <stan/math/fwd/fun/trigamma.hpp>
 #include <stan/math/fwd/fun/trunc.hpp>

stan/math/fwd/fun/trace_dot.hpp

@@ -0,0 +1,37 @@
+#ifndef STAN_MATH_FWD_FUN_TRACE_DOT_HPP
+#define STAN_MATH_FWD_FUN_TRACE_DOT_HPP
+
+#include <stan/math/fwd/core.hpp>
+#include <stan/math/fwd/fun/multiply.hpp>
+#include <stan/math/prim/err.hpp>
+#include <stan/math/prim/fun/trace.hpp>
+
+namespace stan {
+namespace math {
+
+/**
+ * Compute the trace of the product of two matrices with
+ * forward-mode autodiff support.
+ *
+ * @tparam EigMat1 A type either inheriting from `Eigen::DenseBase` or a
+ * `var_value` with an inner type inheriting from `Eigen::DenseBase`
+ * @tparam EigMat2 A type either inheriting from `Eigen::DenseBase` or a
+ * `var_value` with an inner type inheriting from `Eigen::DenseBase`
+ *
+ * @param A first matrix (m x n)
+ * @param B second matrix (n x m)
+ * @return trace of A * B
+ * @throw std::invalid_argument if A and B have incompatible dimensions
+ */
+template <typename EigMat1, typename EigMat2,
+          require_all_eigen_t<EigMat1, EigMat2>* = nullptr,
+          require_any_vt_fvar<EigMat1, EigMat2>* = nullptr>
+inline return_type_t<EigMat1, EigMat2> trace_dot(EigMat1&& A, EigMat2&& B) {
+  check_size_match("trace_dot", "A.cols()", A.cols(), "B.rows()", B.rows());
+  check_size_match("trace_dot", "A.rows()", A.rows(), "B.cols()", B.cols());
+  return trace(multiply(std::forward<EigMat1>(A), std::forward<EigMat2>(B)));
+}
+
+}  // namespace math
+}  // namespace stan
+#endif

stan/math/fwd/prob.hpp

@@ -5,5 +5,6 @@
 #include <stan/math/fwd/fun/Eigen_NumTraits.hpp>
 
 #include <stan/math/fwd/prob/std_normal_log_qf.hpp>
+#include <stan/math/fwd/prob/student_t_qf.hpp>
 
 #endif

stan/math/fwd/prob/student_t_qf.hpp

@@ -0,0 +1,94 @@
+#ifndef STAN_MATH_FWD_PROB_STUDENT_T_QF_HPP
+#define STAN_MATH_FWD_PROB_STUDENT_T_QF_HPP
+
+#include <stan/math/fwd/meta.hpp>
+#include <stan/math/fwd/fun/digamma.hpp>
+#include <stan/math/fwd/fun/exp.hpp>
+#include <stan/math/fwd/fun/hypergeometric_2F1.hpp>
+#include <stan/math/fwd/fun/hypergeometric_pFq.hpp>
+#include <stan/math/fwd/fun/inv_inc_beta.hpp>
+#include <stan/math/fwd/fun/log.hpp>
+#include <stan/math/fwd/fun/sqrt.hpp>
+#include <stan/math/fwd/fun/value_of.hpp>
+#include <stan/math/fwd/fun/value_of_rec.hpp>
+#include <stan/math/prim/meta.hpp>
+#include <stan/math/prim/prob/student_t_lpdf.hpp>
+
+namespace stan {
+namespace math {
+
+template <typename T_p, typename T_nu, typename T_mu, typename T_sigma,
+          require_all_stan_scalar_t<T_p, T_mu, T_sigma, T_nu>* = nullptr,
+          require_any_fvar_t<T_p, T_nu, T_mu, T_sigma>* = nullptr>
+inline auto student_t_qf(const T_p& p, const T_nu& nu, const T_mu& mu,
+                         const T_sigma& sigma) {
+  static constexpr const char* function = "student_t_qf";
+  using FvarT = return_type_t<T_p, T_mu, T_sigma, T_nu>;
+  using T_partials = partials_type_t<FvarT>;
+
+  auto p_val = value_of(p);
+  auto nu_val = value_of(nu);
+  auto mu_val = value_of(mu);
+  auto sigma_val = value_of(sigma);
+
+  check_nonnegative(function, "Degrees of freedom parameter", nu_val);
+  check_positive(function, "Scale parameter", sigma_val);
+  check_bounded(function, "Probability parameter", p_val, 0.0, 1.0);
+
+  if (unlikely(p_val == 0.0)) {
+    return FvarT{NEGATIVE_INFTY, 0.0};
+  } else if (unlikely(p_val == 1.0)) {
+    return FvarT{INFTY, 0.0};
+  } else if (unlikely(p_val == 0.5)) {
+    return FvarT{mu_val, 0.0};
+  }
+
+  const auto p_val_flip = p_val < 0.5 ? p_val : 1.0 - p_val;
+  const double p_sign = value_of_rec(p_val) < 0.5 ? -1.0 : 1.0;
+  auto sqrt_nu_val = sqrt(nu_val);
+  auto ibeta_arg = inv_inc_beta(0.5 * nu_val, 0.5, 2.0 * p_val_flip);
+  auto rtn_val
+      = mu_val
+        + p_sign * sigma_val * sqrt_nu_val * sqrt(-1.0 + 1.0 / ibeta_arg);
+
+  FvarT rtn(rtn_val, 0.0);
+
+  if constexpr (is_autodiff_v<T_p>) {
+    rtn.d_ += p.d_ * exp(-student_t_lpdf(rtn_val, nu_val, mu_val, sigma_val));
+  }
+
+  if constexpr (is_autodiff_v<T_nu>) {
+    const auto half_nu = nu_val / 2.0;
+    Eigen::Matrix<T_partials, -1, 1> hyper_arg_a{{0.5, half_nu, half_nu}};
+    Eigen::Matrix<T_partials, -1, 1> hyper_arg_b{
+        {1.0 + half_nu, 1.0 + half_nu}};
+    const auto hyper_arg
+        = hypergeometric_pFq(hyper_arg_a, hyper_arg_b, ibeta_arg);
+    const auto hyper2f1 = hypergeometric_2F1(1.0, (1.0 + nu_val) / 2.0,
+                                             (2.0 + nu_val) / 2.0, ibeta_arg);
+    const auto digamma_a1 = digamma(half_nu);
+    const auto digamma_a2 = digamma((1.0 + nu_val) / 2.0);
+
+    const auto arg_1 = (4.0 * hyper_arg * sqrt(1.0 - ibeta_arg)) / nu_val;
+    const auto arg_2 = -2.0 * hyper2f1 * (-1.0 + ibeta_arg)
+                       * (log(ibeta_arg) - digamma_a1 + digamma_a2);
+
+    const auto num1 = sigma_val * (-2.0 + (2.0 - arg_1 + arg_2) / ibeta_arg);
+    const auto den1 = 4.0 * sqrt_nu_val * sqrt(-1.0 + 1.0 / ibeta_arg);
+    rtn.d_ += nu.d_ * p_sign * num1 / den1;
+  }
+
+  if constexpr (is_autodiff_v<T_mu>) {
+    rtn.d_ += mu.d_;
+  }
+
+  if constexpr (is_autodiff_v<T_sigma>) {
+    rtn.d_ += sigma.d_ * p_sign * sqrt_nu_val * sqrt(-1.0 + 1.0 / ibeta_arg);
+  }
+
+  return rtn;
+}
+}  // namespace math
+}  // namespace stan
+
+#endif

stan/math/prim/fun.hpp

@@ -328,6 +328,7 @@
 #include <stan/math/prim/fun/to_row_vector.hpp>
 #include <stan/math/prim/fun/to_vector.hpp>
 #include <stan/math/prim/fun/trace.hpp>
+#include <stan/math/prim/fun/trace_dot.hpp>
 #include <stan/math/prim/fun/trace_gen_inv_quad_form_ldlt.hpp>
 #include <stan/math/prim/fun/trace_gen_quad_form.hpp>
 #include <stan/math/prim/fun/trace_inv_quad_form_ldlt.hpp>

stan/math/prim/fun/hypergeometric_pFq.hpp

@@ -4,6 +4,7 @@
 #include <stan/math/prim/meta.hpp>
 #include <stan/math/prim/err/check_not_nan.hpp>
 #include <stan/math/prim/err/check_finite.hpp>
+#include <stan/math/prim/fun/to_ref.hpp>
 #include <stan/math/prim/fun/to_row_vector.hpp>
 #include <boost/math/special_functions/hypergeometric_pFq.hpp>
 
@@ -25,34 +26,53 @@ namespace math {
 template <typename Ta, typename Tb, typename Tz,
           require_all_vector_st<std::is_arithmetic, Ta, Tb>* = nullptr,
           require_arithmetic_t<Tz>* = nullptr>
-inline return_type_t<Ta, Tb, Tz> hypergeometric_pFq(const Ta& a, const Tb& b,
-                                                    const Tz& z) {
-  plain_type_t<Ta> a_ref = a;
-  plain_type_t<Tb> b_ref = b;
+inline return_type_t<Ta, Tb, Tz> hypergeometric_pFq(Ta&& a, Tb&& b, Tz&& z) {
+  decltype(auto) a_ref = to_ref(std::forward<Ta>(a));
+  decltype(auto) b_ref = to_ref(std::forward<Tb>(b));
   check_finite("hypergeometric_pFq", "a", a_ref);
   check_finite("hypergeometric_pFq", "b", b_ref);
   check_finite("hypergeometric_pFq", "z", z);
-
   check_not_nan("hypergeometric_pFq", "a", a_ref);
   check_not_nan("hypergeometric_pFq", "b", b_ref);
   check_not_nan("hypergeometric_pFq", "z", z);
 
-  bool condition_1 = (a_ref.size() > (b_ref.size() + 1)) && (z != 0);
-  bool condition_2 = (a_ref.size() == (b_ref.size() + 1)) && (std::fabs(z) > 1);
+  const bool condition_1 = (a_ref.size() > (b_ref.size() + 1)) && (z != 0);
+  const bool condition_2
+      = (a_ref.size() == (b_ref.size() + 1)) && (std::fabs(z) > 1);
 
   if (condition_1 || condition_2) {
-    std::stringstream msg;
-    msg << "hypergeometric function pFq does not meet convergence "
-        << "conditions with given arguments. "
-        << "a: " << to_row_vector(a_ref) << ", "
-        << "b: " << to_row_vector(b_ref) << ", "
-        << "z: " << z;
-    throw std::domain_error(msg.str());
+    [&]() STAN_COLD_PATH {
+      std::stringstream msg;
+      msg << "hypergeometric function pFq does not meet convergence "
+             "conditions with given arguments. "
+             "a: "
+          << to_row_vector(a_ref) << ", "
+          << "b: " << to_row_vector(b_ref) << ", "
+          << "z: " << z;
+      throw std::domain_error(msg.str());
+    }();
+  }
+  // For plain vectors, we can use Eigen's Map to avoid unnecessary copies
+  using a_ref_t = decltype(a_ref);
+  using b_ref_t = decltype(b_ref);
+  constexpr bool is_a_plain_vec
+      = std::is_same_v<std::decay_t<a_ref_t>, plain_type_t<a_ref_t>>;
+  constexpr bool is_b_plain_vec
+      = std::is_same_v<std::decay_t<b_ref_t>, plain_type_t<b_ref_t>>;
+  if constexpr (is_a_plain_vec && is_b_plain_vec) {
+    // We use type erasure not do a hard copy here
+    using map_t = Eigen::Map<Eigen::VectorXd>;
+    auto map_a = map_t(const_cast<double*>(a_ref.data()), a_ref.size());
+    auto map_b = map_t(const_cast<double*>(b_ref.data()), b_ref.size());
+    return boost::math::hypergeometric_pFq(map_a, map_b, z);
+  } else {
+    // We need pointers to `a` and `b`'s data here so we hard evaluate.
+    decltype(auto) a_eval = eval(a_ref);
+    decltype(auto) b_eval = eval(b_ref);
+    return boost::math::hypergeometric_pFq(
+        std::vector<double>(a_eval.data(), a_eval.data() + a_eval.size()),
+        std::vector<double>(b_eval.data(), b_eval.data() + b_eval.size()), z);
   }
-
-  return boost::math::hypergeometric_pFq(
-      std::vector<double>(a_ref.data(), a_ref.data() + a_ref.size()),
-      std::vector<double>(b_ref.data(), b_ref.data() + b_ref.size()), z);
 }
 }  // namespace math
 }  // namespace stan

stan/math/prim/fun/trace.hpp

@@ -19,8 +19,8 @@ namespace math {
  */
 template <typename T, require_eigen_t<T>* = nullptr,
           require_not_st_var<T>* = nullptr>
-inline value_type_t<T> trace(const T& m) {
-  return m.trace();
+inline auto trace(T&& m) {
+  return make_holder([](auto&& m_) { return m_.trace(); }, std::forward<T>(m));
 }
 
 }  // namespace math

stan/math/prim/fun/trace_dot.hpp

@@ -0,0 +1,43 @@
+#ifndef STAN_MATH_PRIM_FUN_TRACE_DOT_HPP
+#define STAN_MATH_PRIM_FUN_TRACE_DOT_HPP
+
+#include <stan/math/prim/meta.hpp>
+#include <stan/math/prim/err.hpp>
+#include <stan/math/prim/fun/Eigen.hpp>
+
+namespace stan {
+namespace math {
+
+/**
+ * Compute the trace of the product of two matrices,
+ * \f$ \text{tr}(A \cdot B) = \sum_{i,j} A_{ij} B_{ji} \f$.
+ *
+ * This is more efficient than computing the full product and
+ * taking the trace, as it avoids forming the intermediate matrix.
+ *
+ * @tparam EigMat1 A type either inheriting from `Eigen::DenseBase` or a
+ * `var_value` with an inner type inheriting from `Eigen::DenseBase`
+ * @tparam EigMat2 A type either inheriting from `Eigen::DenseBase` or a
+ * `var_value` with an inner type inheriting from `Eigen::DenseBase`
+ *
+ * @param A first matrix (m x n)
+ * @param B second matrix (n x m)
+ * @return trace of A * B
+ * @throw std::invalid_argument if A and B have incompatible dimensions
+ */
+template <typename EigMat1, typename EigMat2,
+          require_all_eigen_vt<std::is_arithmetic, EigMat1, EigMat2>* = nullptr>
+inline auto trace_dot(EigMat1&& A, EigMat2&& B) {
+  check_size_match("trace_dot", "A.cols()", A.cols(), "B.rows()", B.rows());
+  check_size_match("trace_dot", "A.rows()", A.rows(), "B.cols()", B.cols());
+  return make_holder(
+      [](auto&& A_, auto&& B_) {
+        return A_.cwiseProduct(B_.transpose()).sum();
+      },
+      std::forward<EigMat1>(A), std::forward<EigMat2>(B));
+}
+
+}  // namespace math
+}  // namespace stan
+
+#endif

stan/math/prim/meta.hpp

@@ -71,6 +71,7 @@
 #include <stan/math/prim/meta/ad_promotable.hpp>
 #include <stan/math/prim/meta/append_return_type.hpp>
 #include <stan/math/prim/meta/base_type.hpp>
+#include <stan/math/prim/meta/common_container_type.hpp>
 #include <stan/math/prim/meta/contains_std_vector.hpp>
 #include <stan/math/prim/meta/contains_tuple.hpp>
 #include <stan/math/prim/meta/error_index.hpp>