Merge branch 'stan-dev:develop' into feature/issue-3121-beta-neg-binomial-rng

Author: lingium

stan/math/prim/meta/is_eigen.hpp

@@ -20,6 +20,28 @@ template <typename T>
 struct is_eigen
     : bool_constant<is_base_pointer_convertible<Eigen::EigenBase, T>::value> {};
 
+namespace internal {
+// primary template handles types that have no nested ::type member:
+template <class, class = void>
+struct has_internal_trait : std::false_type {};
+
+// specialization recognizes types that do have a nested ::type member:
+template <class T>
+struct has_internal_trait<T,
+                          std::void_t<Eigen::internal::traits<std::decay_t<T>>>>
+    : std::true_type {};
+
+// primary template handles types that have no nested ::type member:
+template <class, class = void>
+struct has_scalar_trait : std::false_type {};
+
+// specialization recognizes types that do have a nested ::type member:
+template <class T>
+struct has_scalar_trait<T, std::void_t<typename std::decay_t<T>::Scalar>>
+    : std::true_type {};
+
+}  // namespace internal
+
 /**
  * Template metaprogram defining the base scalar type of
  * values stored in an Eigen matrix.
@@ -28,7 +50,9 @@ struct is_eigen
  * @ingroup type_trait
  */
 template <typename T>
-struct scalar_type<T, std::enable_if_t<is_eigen<T>::value>> {
+struct scalar_type<T,
+                   std::enable_if_t<is_eigen<T>::value
+                                    && internal::has_scalar_trait<T>::value>> {
   using type = scalar_type_t<typename std::decay_t<T>::Scalar>;
 };
 
@@ -40,10 +64,41 @@ struct scalar_type<T, std::enable_if_t<is_eigen<T>::value>> {
  * @ingroup type_trait
  */
 template <typename T>
-struct value_type<T, std::enable_if_t<is_eigen<T>::value>> {
+struct value_type<T,
+                  std::enable_if_t<is_eigen<T>::value
+                                   && internal::has_scalar_trait<T>::value>> {
   using type = typename std::decay_t<T>::Scalar;
 };
 
+/**
+ * Template metaprogram defining the base scalar type of
+ * values stored in an Eigen matrix.
+ *
+ * @tparam T type to check.
+ * @ingroup type_trait
+ */
+template <typename T>
+struct scalar_type<T,
+                   std::enable_if_t<is_eigen<T>::value
+                                    && !internal::has_scalar_trait<T>::value>> {
+  using type = scalar_type_t<
+      typename Eigen::internal::traits<std::decay_t<T>>::Scalar>;
+};
+
+/**
+ * Template metaprogram defining the type of values stored in an
+ * Eigen matrix, vector, or row vector.
+ *
+ * @tparam T type to check
+ * @ingroup type_trait
+ */
+template <typename T>
+struct value_type<T,
+                  std::enable_if_t<is_eigen<T>::value
+                                   && !internal::has_scalar_trait<T>::value>> {
+  using type = typename Eigen::internal::traits<std::decay_t<T>>::Scalar;
+};
+
 /*! \ingroup require_eigens_types */
 /*! \defgroup eigen_types eigen  */
 /*! \addtogroup eigen_types */

stan/math/prim/prob/beta_neg_binomial_cdf.hpp

@@ -101,7 +101,7 @@ inline return_type_t<T_r, T_alpha, T_beta> beta_neg_binomial_cdf(
         1.0);
     auto C = lgamma(r_plus_n + 1.0) + lbeta(a_plus_r, b_plus_n + 1.0)
              - lgamma(r_dbl) - lbeta(alpha_dbl, beta_dbl) - lgamma(n_dbl + 2.0);
-    auto ccdf = stan::math::exp(C) * F;
+    auto ccdf = stan::math::exp(C + stan::math::log(F));
     cdf *= 1.0 - ccdf;
 
     if constexpr (!is_constant_all<T_r, T_alpha, T_beta>::value) {
@@ -116,15 +116,17 @@ inline return_type_t<T_r, T_alpha, T_beta> beta_neg_binomial_cdf(
       if constexpr (!is_constant<T_r>::value || !is_constant<T_alpha>::value) {
         auto digamma_r_alpha = digamma(a_plus_r);
         if constexpr (!is_constant<T_r>::value) {
-          auto partial_lccdf = digamma(r_plus_n + 1.0)
-                               + (digamma_r_alpha - digamma_n_r_alpha_beta)
-                               + (dF[2] + dF[4]) / F - digamma(r_dbl);
-          partials<0>(ops_partials)[i] += partial_lccdf * chain_rule_term;
+          partials<0>(ops_partials)[i]
+              += (digamma(r_plus_n + 1)
+                  + (digamma_r_alpha - digamma_n_r_alpha_beta)
+                  + (dF[2] + dF[4]) / F - digamma(r_dbl))
+                 * chain_rule_term;
         }
         if constexpr (!is_constant<T_alpha>::value) {
-          auto partial_lccdf = digamma_r_alpha - digamma_n_r_alpha_beta
-                               + dF[4] / F - digamma(alpha_dbl);
-          partials<1>(ops_partials)[i] += partial_lccdf * chain_rule_term;
+          partials<1>(ops_partials)[i]
+              += (digamma_r_alpha - digamma_n_r_alpha_beta + dF[4] / F
+                  - digamma(alpha_dbl))
+                 * chain_rule_term;
         }
       }
 
@@ -135,10 +137,11 @@ inline return_type_t<T_r, T_alpha, T_beta> beta_neg_binomial_cdf(
           partials<1>(ops_partials)[i] += digamma_alpha_beta * chain_rule_term;
         }
         if constexpr (!is_constant<T_beta>::value) {
-          auto partial_lccdf = digamma(b_plus_n + 1.0) - digamma_n_r_alpha_beta
-                               + (dF[1] + dF[4]) / F
-                               - (digamma(beta_dbl) - digamma_alpha_beta);
-          partials<2>(ops_partials)[i] += partial_lccdf * chain_rule_term;
+          partials<2>(ops_partials)[i]
+              += (digamma(b_plus_n + 1) - digamma_n_r_alpha_beta
+                  + (dF[1] + dF[4]) / F
+                  - (digamma(beta_dbl) - digamma_alpha_beta))
+                 * chain_rule_term;
         }
       }
     }

stan/math/rev/core/arena_matrix.hpp

@@ -482,6 +482,7 @@ namespace internal {
 template <typename T>
 struct traits<stan::math::arena_matrix<T>> {
   using base = traits<Eigen::Map<T>>;
+  using Scalar = typename base::Scalar;
   using XprKind = typename Eigen::internal::traits<std::decay_t<T>>::XprKind;
   enum {
     PlainObjectTypeInnerSize = base::PlainObjectTypeInnerSize,

test/unit/math/prim/meta/value_type_test.cpp

@@ -1,5 +1,6 @@
-#include <stan/math/prim/meta.hpp>
 #include <test/unit/util.hpp>
+#include <stan/math/prim/meta.hpp>
+#include <unsupported/Eigen/KroneckerProduct>
 #include <gtest/gtest.h>
 #include <vector>
 
@@ -8,16 +9,16 @@ TEST(MathMetaPrim, value_type_vector) {
   using std::vector;
 
   EXPECT_SAME_TYPE(vector<double>::value_type,
-                   value_type<vector<double> >::type);
+                   value_type<vector<double>>::type);
 
   EXPECT_SAME_TYPE(vector<double>::value_type,
-                   value_type<const vector<double> >::type);
+                   value_type<const vector<double>>::type);
 
-  EXPECT_SAME_TYPE(vector<vector<int> >::value_type,
-                   value_type<vector<vector<int> > >::type);
+  EXPECT_SAME_TYPE(vector<vector<int>>::value_type,
+                   value_type<vector<vector<int>>>::type);
 
-  EXPECT_SAME_TYPE(vector<vector<int> >::value_type,
-                   value_type<const vector<vector<int> > >::type);
+  EXPECT_SAME_TYPE(vector<vector<int>>::value_type,
+                   value_type<const vector<vector<int>>>::type);
 }
 
 TEST(MathMetaPrim, value_type_matrix) {
@@ -33,5 +34,14 @@ TEST(MathMetaPrim, value_type_matrix) {
                    value_type<Eigen::RowVectorXd>::type);
 
   EXPECT_SAME_TYPE(Eigen::RowVectorXd,
-                   value_type<std::vector<Eigen::RowVectorXd> >::type);
+                   value_type<std::vector<Eigen::RowVectorXd>>::type);
+}
+
+TEST(MathMetaPrim, value_type_kronecker) {
+  Eigen::Matrix<double, 2, 2> A;
+  const auto B
+      = Eigen::kroneckerProduct(A, Eigen::Matrix<double, 2, 2>::Identity());
+  Eigen::Matrix<double, 4, 1> C = Eigen::Matrix<double, 4, 1>::Random(4, 1);
+  EXPECT_TRUE((std::is_same<double, stan::value_type_t<decltype(B)>>::value));
+  Eigen::MatrixXd D = B * C;
 }