Merge commit '6025c339b0fad05ac568a678d38cc6e74cddf57b' into HEAD

Author: yashikno

make/compiler_flags

@@ -162,6 +162,10 @@ endif
 ## makes reentrant version lgamma_r available from cmath
 CXXFLAGS_OS += -D_REENTRANT
 
+ifeq (gcc,$(CXX_TYPE))
+CXXFLAGS_WARNINGS += -Wno-int-in-bool-context -Wno-attributes
+endif
+
 ################################################################################
 # Setup OpenCL
 #

stan/math/fwd/fun.hpp

@@ -93,6 +93,7 @@
 #include <stan/math/fwd/fun/proj.hpp>
 #include <stan/math/fwd/fun/quad_form.hpp>
 #include <stan/math/fwd/fun/quad_form_sym.hpp>
+#include <stan/math/fwd/fun/read_fvar.hpp>
 #include <stan/math/fwd/fun/rising_factorial.hpp>
 #include <stan/math/fwd/fun/round.hpp>
 #include <stan/math/fwd/fun/sin.hpp>

stan/math/fwd/fun/Eigen_NumTraits.hpp

@@ -3,6 +3,7 @@
 
 #include <stan/math/prim/fun/Eigen.hpp>
 #include <stan/math/prim/core.hpp>
+#include <stan/math/fwd/fun/read_fvar.hpp>
 #include <stan/math/fwd/core.hpp>
 #include <stan/math/fwd/core/std_numeric_limits.hpp>
 #include <limits>
@@ -163,5 +164,17 @@ struct ScalarBinaryOpTraits<std::complex<stan::math::fvar<T>>,
   using ReturnType = std::complex<stan::math::fvar<T>>;
 };
 
+namespace internal {
+
+/**
+ * Enable linear access of inputs when using read_fvar.
+ */
+template <typename EigFvar, typename EigOut>
+struct functor_has_linear_access<
+    stan::math::read_fvar_functor<EigFvar, EigOut>> {
+  enum { ret = 1 };
+};
+
+}  // namespace internal
 }  // namespace Eigen
 #endif

stan/math/fwd/fun/quad_form.hpp

@@ -14,47 +14,43 @@ namespace math {
  * Symmetry of the resulting matrix is not guaranteed due to numerical
  * precision.
  *
- * @tparam Ta type of elements in the square matrix
- * @tparam Ra number of rows in the square matrix, can be Eigen::Dynamic
- * @tparam Ca number of columns in the square matrix, can be Eigen::Dynamic
- * @tparam Tb type of elements in the second matrix
- * @tparam Rb number of rows in the second matrix, can be Eigen::Dynamic
- * @tparam Cb number of columns in the second matrix, can be Eigen::Dynamic
+ * @tparam EigMat1 type of the first (square) matrix
+ * @tparam EigMat2 type of the second matrix
  *
  * @param A square matrix
  * @param B second matrix
  * @return The quadratic form, which is a symmetric matrix of size Cb.
  * @throws std::invalid_argument if A is not square, or if A cannot be
  * multiplied by B
  */
-template <typename Ta, int Ra, int Ca, typename Tb, int Rb, int Cb,
-          require_any_fvar_t<Ta, Tb>...>
-inline Eigen::Matrix<return_type_t<Ta, Tb>, Cb, Cb> quad_form(
-    const Eigen::Matrix<Ta, Ra, Ca>& A, const Eigen::Matrix<Tb, Rb, Cb>& B) {
+template <typename EigMat1, typename EigMat2,
+          require_all_eigen_t<EigMat1, EigMat2>* = nullptr,
+          require_not_eigen_col_vector_t<EigMat2>* = nullptr,
+          require_any_vt_fvar<EigMat1, EigMat2>* = nullptr>
+inline promote_scalar_t<return_type_t<EigMat1, EigMat2>, EigMat2> quad_form(
+    const EigMat1& A, const EigMat2& B) {
   check_square("quad_form", "A", A);
   check_multiplicable("quad_form", "A", A, "B", B);
-  return multiply(transpose(B), multiply(A, B));
+  return multiply(B.transpose(), multiply(A, B));
 }
 
 /**
  * Return the quadratic form \f$ B^T A B \f$.
  *
- * @tparam Ta type of elements in the square matrix
- * @tparam Ra number of rows in the square matrix, can be Eigen::Dynamic
- * @tparam Ca number of columns in the square matrix, can be Eigen::Dynamic
- * @tparam Tb type of elements in the vector
- * @tparam Rb number of rows in the vector, can be Eigen::Dynamic
+ * @tparam EigMat type of the matrix
+ * @tparam ColVec type of the vector
  *
  * @param A square matrix
  * @param B vector
  * @return The quadratic form (a scalar).
  * @throws std::invalid_argument if A is not square, or if A cannot be
  * multiplied by B
  */
-template <typename Ta, int Ra, int Ca, typename Tb, int Rb,
-          require_any_fvar_t<Ta, Tb>...>
-inline return_type_t<Ta, Tb> quad_form(const Eigen::Matrix<Ta, Ra, Ca>& A,
-                                       const Eigen::Matrix<Tb, Rb, 1>& B) {
+template <typename EigMat, typename ColVec, require_eigen_t<EigMat>* = nullptr,
+          require_eigen_col_vector_t<ColVec>* = nullptr,
+          require_any_vt_fvar<EigMat, ColVec>* = nullptr>
+inline return_type_t<EigMat, ColVec> quad_form(const EigMat& A,
+                                               const ColVec& B) {
   check_square("quad_form", "A", A);
   check_multiplicable("quad_form", "A", A, "B", B);
   return dot_product(B, multiply(A, B));

stan/math/fwd/fun/quad_form_sym.hpp

@@ -13,49 +13,45 @@ namespace math {
  *
  * Symmetry of the resulting matrix is guaranteed.
  *
- * @tparam TA type of elements in the symmetric matrix
- * @tparam RA number of rows in the symmetric matrix, can be Eigen::Dynamic
- * @tparam CA number of columns in the symmetric matrix, can be Eigen::Dynamic
- * @tparam TB type of elements in the second matrix
- * @tparam RB number of rows in the second matrix, can be Eigen::Dynamic
- * @tparam CB number of columns in the second matrix, can be Eigen::Dynamic
+ * @tparam EigMat1 type of the first (symmetric) matrix
+ * @tparam EigMat2 type of the second matrix
  *
  * @param A symmetric matrix
  * @param B second matrix
  * @return The quadratic form, which is a symmetric matrix of size CB.
  * @throws std::invalid_argument if A is not symmetric, or if A cannot be
  * multiplied by B
  */
-template <typename TA, int RA, int CA, typename TB, int RB, int CB,
-          require_any_fvar_t<TA, TB>...>
-inline Eigen::Matrix<return_type_t<TA, TB>, CB, CB> quad_form_sym(
-    const Eigen::Matrix<TA, RA, CA>& A, const Eigen::Matrix<TB, RB, CB>& B) {
-  using T = return_type_t<TA, TB>;
+template <typename EigMat1, typename EigMat2,
+          require_all_eigen_t<EigMat1, EigMat2>* = nullptr,
+          require_not_eigen_col_vector_t<EigMat2>* = nullptr,
+          require_any_vt_fvar<EigMat1, EigMat2>* = nullptr>
+inline promote_scalar_t<return_type_t<EigMat1, EigMat2>, EigMat2> quad_form_sym(
+    const EigMat1& A, const EigMat2& B) {
+  using T_ret = return_type_t<EigMat1, EigMat2>;
   check_multiplicable("quad_form_sym", "A", A, "B", B);
   check_symmetric("quad_form_sym", "A", A);
-  Eigen::Matrix<T, CB, CB> ret(multiply(transpose(B), multiply(A, B)));
-  return T(0.5) * (ret + transpose(ret));
+  promote_scalar_t<T_ret, EigMat2> ret(multiply(B.transpose(), multiply(A, B)));
+  return T_ret(0.5) * (ret + ret.transpose());
 }
 
 /**
  * Return the quadratic form \f$ B^T A B \f$ of a symmetric matrix.
  *
- * @tparam TA type of elements in the symmetric matrix
- * @tparam RA number of rows in the symmetric matrix, can be Eigen::Dynamic
- * @tparam CA number of columns in the symmetric matrix, can be Eigen::Dynamic
- * @tparam TB type of elements in the vector
- * @tparam RB number of rows in the vector, can be Eigen::Dynamic
+ * @tparam EigMat type of the (symmetric) matrix
+ * @tparam ColVec type of the vector
  *
  * @param A symmetric matrix
  * @param B vector
  * @return The quadratic form (a scalar).
  * @throws std::invalid_argument if A is not symmetric, or if A cannot be
  * multiplied by B
  */
-template <typename TA, int RA, int CA, typename TB, int RB,
-          require_any_fvar_t<TA, TB>...>
-inline return_type_t<TA, TB> quad_form_sym(const Eigen::Matrix<TA, RA, CA>& A,
-                                           const Eigen::Matrix<TB, RB, 1>& B) {
+template <typename EigMat, typename ColVec, require_eigen_t<EigMat>* = nullptr,
+          require_eigen_col_vector_t<ColVec>* = nullptr,
+          require_any_vt_fvar<EigMat, ColVec>* = nullptr>
+inline return_type_t<EigMat, ColVec> quad_form_sym(const EigMat& A,
+                                                   const ColVec& B) {
   check_multiplicable("quad_form_sym", "A", A, "B", B);
   check_symmetric("quad_form_sym", "A", A);
   return dot_product(B, multiply(A, B));

stan/math/fwd/fun/read_fvar.hpp

@@ -0,0 +1,52 @@
+#ifndef STAN_MATH_FWD_FUN_READ_FVAR_HPP
+#define STAN_MATH_FWD_FUN_READ_FVAR_HPP
+
+#include <stan/math/fwd/meta.hpp>
+
+namespace stan {
+namespace math {
+
+/**
+ * Functor for extracting the values and tangents from a matrix of fvar.
+ * This functor is called using Eigen's NullaryExpr framework.
+ */
+template <typename EigFvar, typename EigOut>
+class read_fvar_functor {
+  const EigFvar& var_mat;
+  EigOut& val_mat;
+
+ public:
+  read_fvar_functor(const EigFvar& arg1, EigOut& arg2)
+      : var_mat(arg1), val_mat(arg2) {}
+
+  inline decltype(auto) operator()(Eigen::Index row, Eigen::Index col) const {
+    val_mat.coeffRef(row, col) = var_mat.coeffRef(row, col).val_;
+    return var_mat.coeffRef(row, col).d_;
+  }
+
+  inline decltype(auto) operator()(Eigen::Index index) const {
+    val_mat.coeffRef(index) = var_mat.coeffRef(index).val_;
+    return var_mat.coeffRef(index).d_;
+  }
+};
+
+/**
+ * Function applying the read_fvar_functor to extract the values
+ * and tangets of a given fvar matrix into separate matrices.
+ *
+ * @tparam EigFvar type of the Eigen container of fvar.
+ * @tparam EigOut type of the Eigen containers to copy to
+ * @param[in] FvarMat Input Eigen container of fvar.
+ * @param[in] ValMat Output Eigen container of values.
+ * @param[in] DMat Output Eigen container of tangents.
+ */
+template <typename EigFvar, typename EigOut>
+inline void read_fvar(const EigFvar& FvarMat, EigOut& ValMat, EigOut& DMat) {
+  DMat = EigOut::NullaryExpr(
+      FvarMat.rows(), FvarMat.cols(),
+      read_fvar_functor<const EigFvar, EigOut>(FvarMat, ValMat));
+}
+
+}  // namespace math
+}  // namespace stan
+#endif

stan/math/fwd/fun/softmax.hpp

@@ -8,36 +8,39 @@
 namespace stan {
 namespace math {
 
-template <typename T>
-inline Eigen::Matrix<fvar<T>, Eigen::Dynamic, 1> softmax(
-    const Eigen::Matrix<fvar<T>, Eigen::Dynamic, 1>& alpha) {
+template <typename ColVec,
+          require_eigen_col_vector_vt<is_fvar, ColVec>* = nullptr>
+inline auto softmax(const ColVec& alpha) {
   using Eigen::Dynamic;
   using Eigen::Matrix;
+  using T = typename value_type_t<ColVec>::Scalar;
 
   Matrix<T, Dynamic, 1> alpha_t(alpha.size());
   for (int k = 0; k < alpha.size(); ++k) {
-    alpha_t(k) = alpha(k).val_;
+    alpha_t.coeffRef(k) = alpha.coeff(k).val_;
   }
 
   Matrix<T, Dynamic, 1> softmax_alpha_t = softmax(alpha_t);
 
   Matrix<fvar<T>, Dynamic, 1> softmax_alpha(alpha.size());
   for (int k = 0; k < alpha.size(); ++k) {
-    softmax_alpha(k).val_ = softmax_alpha_t(k);
-    softmax_alpha(k).d_ = 0;
+    softmax_alpha.coeffRef(k).val_ = softmax_alpha_t.coeff(k);
+    softmax_alpha.coeffRef(k).d_ = 0;
   }
 
   for (int m = 0; m < alpha.size(); ++m) {
     T negative_alpha_m_d_times_softmax_alpha_t_m
-        = -alpha(m).d_ * softmax_alpha_t(m);
+        = -alpha.coeff(m).d_ * softmax_alpha_t.coeff(m);
     for (int k = 0; k < alpha.size(); ++k) {
       if (m == k) {
-        softmax_alpha(k).d_
-            += softmax_alpha_t(k)
-               * (alpha(m).d_ + negative_alpha_m_d_times_softmax_alpha_t_m);
+        softmax_alpha.coeffRef(k).d_
+            += softmax_alpha_t.coeff(k)
+               * (alpha.coeff(m).d_
+                  + negative_alpha_m_d_times_softmax_alpha_t_m);
       } else {
-        softmax_alpha(k).d_
-            += negative_alpha_m_d_times_softmax_alpha_t_m * softmax_alpha_t(k);
+        softmax_alpha.coeffRef(k).d_
+            += softmax_alpha_t.coeff(k)
+               * negative_alpha_m_d_times_softmax_alpha_t_m;
       }
     }
   }

stan/math/fwd/fun/trace_quad_form.hpp

@@ -11,29 +11,16 @@
 namespace stan {
 namespace math {
 
-template <int RA, int CA, int RB, int CB, typename T>
-inline fvar<T> trace_quad_form(const Eigen::Matrix<fvar<T>, RA, CA> &A,
-                               const Eigen::Matrix<fvar<T>, RB, CB> &B) {
+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_quad_form(const EigMat1& A,
+                                                       const EigMat2& B) {
   check_square("trace_quad_form", "A", A);
   check_multiplicable("trace_quad_form", "A", A, "B", B);
-  return trace(multiply(transpose(B), multiply(A, B)));
+  return B.cwiseProduct(multiply(A, B)).sum();
 }
 
-template <int RA, int CA, int RB, int CB, typename T>
-inline fvar<T> trace_quad_form(const Eigen::Matrix<fvar<T>, RA, CA> &A,
-                               const Eigen::Matrix<double, RB, CB> &B) {
-  check_square("trace_quad_form", "A", A);
-  check_multiplicable("trace_quad_form", "A", A, "B", B);
-  return trace(multiply(transpose(B), multiply(A, B)));
-}
-
-template <int RA, int CA, int RB, int CB, typename T>
-inline fvar<T> trace_quad_form(const Eigen::Matrix<double, RA, CA> &A,
-                               const Eigen::Matrix<fvar<T>, RB, CB> &B) {
-  check_square("trace_quad_form", "A", A);
-  check_multiplicable("trace_quad_form", "A", A, "B", B);
-  return trace(multiply(transpose(B), multiply(A, B)));
-}
 }  // namespace math
 }  // namespace stan
 

stan/math/opencl/kernel_generator/binary_operation.hpp

@@ -187,8 +187,7 @@ ADD_BINARY_OPERATION(addition_, operator+, common_scalar_t<T_a COMMA T_b>, "+");
 ADD_BINARY_OPERATION(subtraction_, operator-, common_scalar_t<T_a COMMA T_b>,
                      "-");
 ADD_BINARY_OPERATION_WITH_CUSTOM_CODE(
-    elewise_multiplication_, elewise_multiplication,
-    common_scalar_t<T_a COMMA T_b>, "*",
+    elt_multiply_, elt_multiply, common_scalar_t<T_a COMMA T_b>, "*",
     using view_transitivity = std::tuple<std::true_type, std::true_type>;
     inline std::pair<int, int> extreme_diagonals() const {
       std::pair<int, int> diags0
@@ -200,7 +199,7 @@ ADD_BINARY_OPERATION_WITH_CUSTOM_CODE(
     });
 
 ADD_BINARY_OPERATION_WITH_CUSTOM_CODE(
-    elewise_division_, elewise_division, common_scalar_t<T_a COMMA T_b>, "/",
+    elt_divide_, elt_divide, common_scalar_t<T_a COMMA T_b>, "/",
     inline std::pair<int, int> extreme_diagonals() const {
       return {-rows() + 1, cols() - 1};
     });
@@ -246,7 +245,7 @@ ADD_BINARY_OPERATION_WITH_CUSTOM_CODE(
  */
 template <typename T_a, typename T_b, typename = require_arithmetic_t<T_a>,
           typename = require_all_kernel_expressions_t<T_b>>
-inline elewise_multiplication_<scalar_<T_a>, as_operation_cl_t<T_b>> operator*(
+inline elt_multiply_<scalar_<T_a>, as_operation_cl_t<T_b>> operator*(
     T_a&& a, T_b&& b) {  // NOLINT
   return {as_operation_cl(std::forward<T_a>(a)),
           as_operation_cl(std::forward<T_b>(b))};
@@ -263,7 +262,7 @@ inline elewise_multiplication_<scalar_<T_a>, as_operation_cl_t<T_b>> operator*(
 template <typename T_a, typename T_b,
           typename = require_all_kernel_expressions_t<T_a>,
           typename = require_arithmetic_t<T_b>>
-inline elewise_multiplication_<as_operation_cl_t<T_a>, scalar_<T_b>> operator*(
+inline elt_multiply_<as_operation_cl_t<T_a>, scalar_<T_b>> operator*(
     T_a&& a, const T_b b) {  // NOLINT
   return {as_operation_cl(std::forward<T_a>(a)), as_operation_cl(b)};
 }

stan/math/opencl/kernel_generator/matrix_vector_multiply.hpp

@@ -15,7 +15,7 @@ namespace math {
 template <typename T_matrix, typename T_vector,
           typename = require_all_kernel_expressions_t<T_matrix, T_vector>>
 inline auto matrix_vector_multiply(T_matrix&& matrix, T_vector&& vector) {
-  return rowwise_sum(elewise_multiplication(
+  return rowwise_sum(elt_multiply(
       std::forward<T_matrix>(matrix),
       colwise_broadcast(transpose(std::forward<T_vector>(vector)))));
 }