stan changes

stan/math/fwd/fun/accumulator.hpp

@@ -9,4 +9,92 @@
 #include <vector>
 #include <type_traits>
 
+namespace stan {
+namespace math {
+
+/**
+ * Class to accumulate values and eventually return their sum.  If
+ * no values are ever added, the return value is 0.
+ *
+ * This class is useful for speeding up autodiff of long sums
+ * because it uses the <code>sum()</code> operation (either from
+ * <code>stan::math</code> or one defined by argument-dependent lookup.
+ *
+ * @tparam T Type of scalar added
+ */
+template <typename T>
+class accumulator<T, require_fvar_t<T>> {
+ private:
+  std::vector<T> buf_;
+
+ public:
+  /**
+   * Add the specified arithmetic type value to the buffer after
+   * static casting it to the class type <code>T</code>.
+   *
+   * <p>See the std library doc for <code>std::is_arithmetic</code>
+   * for information on what counts as an arithmetic type.
+   *
+   * @tparam S Type of argument
+   * @param x Value to add
+   */
+  template <typename S, typename = require_stan_scalar_t<S>>
+  inline void add(S x) {
+    buf_.push_back(x);
+  }
+
+  /**
+   * Add each entry in the specified matrix, vector, or row vector
+   * of values to the buffer.
+   *
+   * @tparam S type of the matrix
+   * @param m Matrix of values to add
+   */
+  template <typename S, require_matrix_t<S>* = nullptr>
+  inline void add(const S& m) {
+    buf_.push_back(stan::math::sum(m));
+  }
+
+  /**
+   * Recursively add each entry in the specified standard vector
+   * to the buffer.  This will allow vectors of primitives,
+   * autodiff variables to be added; if the vector entries
+   * are collections, their elements are recursively added.
+   *
+   * @tparam S Type of value to recursively add.
+   * @param xs Vector of entries to add
+   */
+  template <typename S>
+  inline void add(const std::vector<S>& xs) {
+    for (size_t i = 0; i < xs.size(); ++i) {
+      this->add(xs[i]);
+    }
+  }
+
+#ifdef STAN_OPENCL
+
+  /**
+   * Sum each entry and then push to the buffer.
+   * @tparam S A Type inheriting from `matrix_cl_base`
+   * @param x An OpenCL matrix
+   */
+  template <typename S,
+            require_all_kernel_expressions_and_none_scalar_t<S>* = nullptr>
+  inline void add(const S& xs) {
+    buf_.push_back(stan::math::sum(xs));
+  }
+
+#endif
+
+  /**
+   * Return the sum of the accumulated values.
+   *
+   * @return Sum of accumulated values.
+   */
+  inline T sum() const { return stan::math::sum(buf_); }
+};
+
+}  // namespace math
+}  // namespace stan
+
 #endif

stan/math/fwd/fun/mdivide_right.hpp

@@ -12,117 +12,89 @@
 
 namespace stan {
 namespace math {
-
+/*
 template <typename EigMat1, typename EigMat2,
-          require_all_eigen_vt<is_fvar, EigMat1, EigMat2>* = nullptr,
-          require_vt_same<EigMat1, EigMat2>* = nullptr>
-inline Eigen::Matrix<value_type_t<EigMat1>, EigMat1::RowsAtCompileTime,
-                     EigMat2::ColsAtCompileTime>
-mdivide_right(const EigMat1& A, const EigMat2& b) {
-  using T = typename value_type_t<EigMat1>::Scalar;
-  constexpr int R1 = EigMat1::RowsAtCompileTime;
-  constexpr int C1 = EigMat1::ColsAtCompileTime;
-  constexpr int R2 = EigMat2::RowsAtCompileTime;
-  constexpr int C2 = EigMat2::ColsAtCompileTime;
-
-  check_square("mdivide_right", "b", b);
-  check_multiplicable("mdivide_right", "A", A, "b", b);
-  if (b.size() == 0) {
-    return {A.rows(), 0};
+          require_all_eigen_vt<is_fvar, EigMat1, EigMat2>* = nullptr>
+inline auto
+mdivide_right(const EigMat1& b, const EigMat2& A) {
+  std::cout << "\nUsing 1: " << "\n";
+  using A_fvar_inner_type = typename value_type_t<EigMat2>::Scalar;
+  using b_fvar_inner_type = typename value_type_t<EigMat1>::Scalar;
+  using inner_ret_t = return_type_t<A_fvar_inner_type, b_fvar_inner_type>;
+  constexpr auto R1 = EigMat1::RowsAtCompileTime;
+  constexpr auto C1 = EigMat1::ColsAtCompileTime;
+  constexpr auto R2 = EigMat2::RowsAtCompileTime;
+  constexpr auto C2 = EigMat2::ColsAtCompileTime;
+
+  check_square("mdivide_right", "A", A);
+  check_multiplicable("mdivide_right", "b", b, "A", A);
+  if (A.size() == 0) {
+    using ret_t = decltype(mdivide_right(b.val(), A.val()).eval());
+    return promote_scalar_t<fvar<inner_ret_t>, ret_t>{b.rows(), 0};
   }
 
-  Eigen::Matrix<T, R1, C1> val_A(A.rows(), A.cols());
-  Eigen::Matrix<T, R1, C1> deriv_A(A.rows(), A.cols());
-  Eigen::Matrix<T, R2, C2> val_b(b.rows(), b.cols());
-  Eigen::Matrix<T, R2, C2> deriv_b(b.rows(), b.cols());
+  Eigen::Matrix<A_fvar_inner_type, R2, C2> val_A(A.rows(), A.cols());
+  Eigen::Matrix<A_fvar_inner_type, R2, C2> deriv_A(A.rows(), A.cols());
 
-  const Eigen::Ref<const plain_type_t<EigMat1>>& A_ref = A;
+  const auto& A_ref = to_ref(A);
   for (int j = 0; j < A.cols(); j++) {
     for (int i = 0; i < A.rows(); i++) {
       val_A.coeffRef(i, j) = A_ref.coeff(i, j).val_;
       deriv_A.coeffRef(i, j) = A_ref.coeff(i, j).d_;
     }
   }
 
-  const Eigen::Ref<const plain_type_t<EigMat2>>& b_ref = b;
-  for (int j = 0; j < b.cols(); j++) {
-    for (int i = 0; i < b.rows(); i++) {
+  Eigen::Matrix<b_fvar_inner_type, R1, C1> val_b(b.rows(), b.cols());
+  Eigen::Matrix<b_fvar_inner_type, R1, C1> deriv_b(b.rows(), b.cols());
+  const auto& b_ref = to_ref(b);
+  for (Eigen::Index j = 0; j < b.cols(); j++) {
+    for (Eigen::Index i = 0; i < b.rows(); i++) {
       val_b.coeffRef(i, j) = b_ref.coeff(i, j).val_;
       deriv_b.coeffRef(i, j) = b_ref.coeff(i, j).d_;
     }
   }
-
-  Eigen::Matrix<T, R1, C2> A_mult_inv_b = mdivide_right(val_A, val_b);
-
-  return to_fvar(A_mult_inv_b,
-                 mdivide_right(deriv_A, val_b)
-                     - A_mult_inv_b * mdivide_right(deriv_b, val_b));
+  auto A_mult_inv_b = mdivide_right(val_b, val_A).eval();
+  promote_scalar_t<fvar<inner_ret_t>, decltype(A_mult_inv_b)>
+ret(A_mult_inv_b.rows(), A_mult_inv_b.cols()); ret.val() = A_mult_inv_b; ret.d()
+= mdivide_right(deriv_b, val_A)
+      - multiply(A_mult_inv_b, mdivide_right(deriv_A, val_A));
+  return ret;
 }
 
 template <typename EigMat1, typename EigMat2,
           require_eigen_vt<is_fvar, EigMat1>* = nullptr,
-          require_eigen_vt<std::is_arithmetic, EigMat2>* = nullptr>
-inline Eigen::Matrix<value_type_t<EigMat1>, EigMat1::RowsAtCompileTime,
-                     EigMat2::ColsAtCompileTime>
-mdivide_right(const EigMat1& A, const EigMat2& b) {
-  using T = typename value_type_t<EigMat1>::Scalar;
-  constexpr int R1 = EigMat1::RowsAtCompileTime;
-  constexpr int C1 = EigMat1::ColsAtCompileTime;
-
-  check_square("mdivide_right", "b", b);
-  check_multiplicable("mdivide_right", "A", A, "b", b);
-  if (b.size() == 0) {
-    return {A.rows(), 0};
-  }
-
-  Eigen::Matrix<T, R1, C1> val_A(A.rows(), A.cols());
-  Eigen::Matrix<T, R1, C1> deriv_A(A.rows(), A.cols());
-
-  const Eigen::Ref<const plain_type_t<EigMat1>>& A_ref = A;
-  for (int j = 0; j < A.cols(); j++) {
-    for (int i = 0; i < A.rows(); i++) {
-      val_A.coeffRef(i, j) = A_ref.coeff(i, j).val_;
-      deriv_A.coeffRef(i, j) = A_ref.coeff(i, j).d_;
-    }
-  }
-
-  return to_fvar(mdivide_right(val_A, b), mdivide_right(deriv_A, b));
-}
+          require_eigen_vt<is_var_or_arithmetic, EigMat2>* = nullptr>
+          inline auto
+          mdivide_right(const EigMat1& b, const EigMat2& A) {
+            using T_return = return_type_t<EigMat1, EigMat2>;
+            check_square("mdivide_right", "A", A);
+            check_multiplicable("mdivide_right", "b", b, "A", A);
+            if (A.size() == 0) {
+              using ret_type = decltype(A.transpose().template
+cast<T_return>().lu().solve(b.template
+cast<T_return>().transpose()).transpose().eval()); return ret_type{b.rows(), 0};
+            }
+            return A.transpose().template cast<T_return>().lu().solve(b.template
+cast<T_return>().transpose()).transpose().eval();
+          }
 
 template <typename EigMat1, typename EigMat2,
-          require_eigen_vt<std::is_arithmetic, EigMat1>* = nullptr,
+          require_eigen_vt<is_var_or_arithmetic, EigMat1>* = nullptr,
           require_eigen_vt<is_fvar, EigMat2>* = nullptr>
-inline Eigen::Matrix<value_type_t<EigMat2>, EigMat1::RowsAtCompileTime,
-                     EigMat2::ColsAtCompileTime>
-mdivide_right(const EigMat1& A, const EigMat2& b) {
-  using T = typename value_type_t<EigMat2>::Scalar;
-  constexpr int R1 = EigMat1::RowsAtCompileTime;
-  constexpr int C1 = EigMat1::ColsAtCompileTime;
-  constexpr int R2 = EigMat2::RowsAtCompileTime;
-  constexpr int C2 = EigMat2::ColsAtCompileTime;
-
-  check_square("mdivide_right", "b", b);
-  check_multiplicable("mdivide_right", "A", A, "b", b);
-  if (b.size() == 0) {
-    return {A.rows(), 0};
-  }
-
-  Eigen::Matrix<T, R2, C2> val_b(b.rows(), b.cols());
-  Eigen::Matrix<T, R2, C2> deriv_b(b.rows(), b.cols());
-
-  const Eigen::Ref<const plain_type_t<EigMat2>>& b_ref = b;
-  for (int j = 0; j < b.cols(); j++) {
-    for (int i = 0; i < b.rows(); i++) {
-      val_b.coeffRef(i, j) = b_ref.coeff(i, j).val_;
-      deriv_b.coeffRef(i, j) = b_ref.coeff(i, j).d_;
-    }
-  }
-
-  Eigen::Matrix<T, R1, C2> A_mult_inv_b = mdivide_right(A, val_b);
-
-  return to_fvar(A_mult_inv_b, -A_mult_inv_b * mdivide_right(deriv_b, val_b));
-}
-
+          inline auto
+          mdivide_right(const EigMat1& b, const EigMat2& A) {
+            using T_return = return_type_t<EigMat1, EigMat2>;
+            check_square("mdivide_right", "A", A);
+            check_multiplicable("mdivide_right", "b", b, "A", A);
+            if (A.size() == 0) {
+              using ret_type = decltype(A.transpose().template
+cast<T_return>().lu().solve(b.template
+cast<T_return>().transpose()).transpose().eval()); return ret_type{b.rows(), 0};
+            }
+            return A.transpose().template cast<T_return>().lu().solve(b.template
+cast<T_return>().transpose()).transpose().eval();
+          }
+*/
 }  // namespace math
 }  // namespace stan
 #endif

stan/math/fwd/fun/sum.hpp

@@ -45,7 +45,7 @@ inline value_type_t<T> sum(const T& m) {
   if (m.size() == 0) {
     return 0.0;
   }
-  const Eigen::Ref<const plain_type_t<T>>& m_ref = m;
+  const auto& m_ref = to_ref(m);
   return {sum(m_ref.val()), sum(m_ref.d())};
 }
 

stan/math/mix.hpp

@@ -5,6 +5,11 @@
 #include <stan/math/mix/fun.hpp>
 #include <stan/math/mix/functor.hpp>
 
+#include <stan/math/rev/core.hpp>
+#include <stan/math/rev/meta.hpp>
+#include <stan/math/rev/fun.hpp>
+#include <stan/math/rev/functor.hpp>
+
 #ifdef STAN_OPENCL
 #include <stan/math/opencl/rev.hpp>
 #endif
@@ -14,11 +19,6 @@
 #include <stan/math/fwd/fun.hpp>
 #include <stan/math/fwd/functor.hpp>
 
-#include <stan/math/rev/core.hpp>
-#include <stan/math/rev/meta.hpp>
-#include <stan/math/rev/fun.hpp>
-#include <stan/math/rev/functor.hpp>
-
 #include <stan/math/prim.hpp>
 
 #endif

stan/math/mix/functor/derivative.hpp

@@ -1,8 +1,8 @@
 #ifndef STAN_MATH_MIX_FUNCTOR_DERIVATIVE_HPP
 #define STAN_MATH_MIX_FUNCTOR_DERIVATIVE_HPP
 
-#include <stan/math/fwd/core.hpp>
 #include <stan/math/prim/fun/Eigen.hpp>
+#include <stan/math/fwd/core.hpp>
 #include <stan/math/rev/core.hpp>
 #include <vector>
 

stan/math/mix/meta.hpp

@@ -6,12 +6,12 @@
 #include <stan/math/fwd/core.hpp>
 #include <stan/math/fwd/meta/is_fvar.hpp>
 #include <stan/math/fwd/meta/partials_type.hpp>
+#include <stan/math/fwd/meta.hpp>
 
 #include <stan/math/rev/core.hpp>
 #include <stan/math/rev/meta/is_var.hpp>
 #include <stan/math/rev/meta/partials_type.hpp>
 
-#include <stan/math/fwd/meta.hpp>
 #include <stan/math/rev/meta.hpp>
 #include <stan/math/prim/meta.hpp>
 

stan/math/prim/fun/eigenvalues.hpp

@@ -7,14 +7,11 @@
 namespace stan {
 namespace math {
 
-template <typename T>
-Eigen::Matrix<std::complex<T>, -1, 1> eigenvalues(
-    const Eigen::Matrix<T, -1, -1>& m) {
+template <typename Mat, require_eigen_t<Mat>* = nullptr>
+inline auto eigenvalues(const Mat& m) {
   check_nonzero_size("eigenvalues", "m", m);
   check_square("eigenvalues", "m", m);
-
-  Eigen::EigenSolver<Eigen::Matrix<T, -1, -1>> solver(m);
-  return solver.eigenvalues();
+  return m.eigenvalues().eval();
 }
 
 }  // namespace math

stan/math/prim/fun/mdivide_left_tri.hpp

@@ -26,20 +26,24 @@ namespace math {
 template <Eigen::UpLoType TriView, typename T1, typename T2,
           require_all_eigen_t<T1, T2> * = nullptr,
           require_all_not_eigen_vt<is_var, T1, T2> * = nullptr>
-inline Eigen::Matrix<return_type_t<T1, T2>, T1::RowsAtCompileTime,
-                     T2::ColsAtCompileTime>
-mdivide_left_tri(const T1 &A, const T2 &b) {
+inline auto mdivide_left_tri(const T1 &A, const T2 &b) {
   using T_return = return_type_t<T1, T2>;
   check_square("mdivide_left_tri", "A", A);
   check_multiplicable("mdivide_left_tri", "A", A, "b", b);
+  using ret_type = decltype(A.template cast<T_return>()
+                                .eval()
+                                .template triangularView<TriView>()
+                                .solve(b.template cast<T_return>().eval())
+                                .eval());
   if (A.rows() == 0) {
-    return {0, b.cols()};
+    return ret_type(0, b.cols());
   }
 
   return A.template cast<T_return>()
       .eval()
       .template triangularView<TriView>()
-      .solve(b.template cast<T_return>().eval());
+      .solve(b.template cast<T_return>().eval())
+      .eval();
 }
 
 /**