From d00537bc3d488fecc4f065bbd240bc9f6ed16334 Mon Sep 17 00:00:00 2001 From: Ben Date: Mon, 26 Oct 2020 13:20:14 -0400 Subject: [PATCH 01/12] `var` versions of quad_form and trace_quad_form (Issue #2101) --- stan/math/prim/err/check_square.hpp | 2 +- stan/math/prim/meta/is_vector.hpp | 10 ++ stan/math/rev/fun/quad_form.hpp | 102 ++++++++++++++++++ stan/math/rev/fun/trace_quad_form.hpp | 30 +++++- test/unit/math/mix/fun/quad_form_test.cpp | 14 +++ .../math/mix/fun/trace_quad_form_test.cpp | 15 +++ test/unit/math/test_ad.hpp | 8 ++ 7 files changed, 179 insertions(+), 2 deletions(-) diff --git a/stan/math/prim/err/check_square.hpp b/stan/math/prim/err/check_square.hpp index ecd509ef91b..2aabff38eee 100644 --- a/stan/math/prim/err/check_square.hpp +++ b/stan/math/prim/err/check_square.hpp @@ -17,7 +17,7 @@ namespace math { * @param y Matrix to test * @throw std::invalid_argument if the matrix is not square */ -template * = nullptr> +template * = nullptr> inline void check_square(const char* function, const char* name, const T_y& y) { check_size_match(function, "Expecting a square matrix; rows of ", name, y.rows(), "columns of ", name, y.cols()); diff --git a/stan/math/prim/meta/is_vector.hpp b/stan/math/prim/meta/is_vector.hpp index 21940ce5c39..4a4c38b5dfe 100644 --- a/stan/math/prim/meta/is_vector.hpp +++ b/stan/math/prim/meta/is_vector.hpp @@ -100,6 +100,11 @@ STAN_ADD_REQUIRE_CONTAINER(eigen_col_vector, is_eigen_col_vector, template struct is_col_vector : internal::is_col_vector_impl {}; +STAN_ADD_REQUIRE_UNARY(col_vector, is_col_vector, + require_eigens_types); +STAN_ADD_REQUIRE_CONTAINER(col_vector, is_col_vector, + require_eigens_types); + /** \ingroup type_trait * If the input type T is an eigen matrix with 1 column at compile time this * has a static member with a value of true. Else this has a static @@ -122,6 +127,11 @@ STAN_ADD_REQUIRE_CONTAINER(eigen_row_vector, is_eigen_row_vector, template struct is_row_vector : internal::is_row_vector_impl {}; +STAN_ADD_REQUIRE_UNARY(row_vector, is_row_vector, + require_eigens_types); +STAN_ADD_REQUIRE_CONTAINER(row_vector, is_row_vector, + require_eigens_types); + /** \ingroup type_trait * If the input type T is an eigen matrix with 1 column or 1 row at compile time * this has a static member with a value of true. Else this has a static diff --git a/stan/math/rev/fun/quad_form.hpp b/stan/math/rev/fun/quad_form.hpp index d930ac33f50..fd4f72a2bf8 100644 --- a/stan/math/rev/fun/quad_form.hpp +++ b/stan/math/rev/fun/quad_form.hpp @@ -93,6 +93,65 @@ class quad_form_vari : public vari { quad_form_vari_alloc* impl_; }; + +/** + * Return the quadratic form \f$ B^T A B \f$. + * + * Symmetry of the resulting matrix is not guaranteed due to numerical + * precision. + * + * @tparam T1 type of the first (square) matrix + * @tparam T2 type of the second matrix + * + * @param A square matrix + * @param B second matrix + * @return The quadratic form, which is a symmetric matrix. + * @throws std::invalid_argument if A is not square, or if A cannot be + * multiplied by B + */ +template * = nullptr, + require_any_var_matrix_t* = nullptr> +inline auto quad_form_impl(const Mat1& A, const Mat2& B, bool symmetric) { + check_square("quad_form", "A", A); + check_multiplicable("quad_form", "A", A, "B", B); + + using return_t + = promote_var_matrix_t; + + auto arena_A = to_arena(A); + auto arena_B = to_arena(B); + + check_not_nan("multiply", "A", value_of(arena_A)); + check_not_nan("multiply", "B", value_of(arena_B)); + + auto arena_res = to_arena(value_of(arena_B).transpose() * + value_of(arena_A) * + value_of(arena_B)); + + if(symmetric) { + arena_res = (arena_res + arena_res.transpose()).eval(); + } + + return_t res = arena_res; + + reverse_pass_callback([arena_A, arena_B, res]() mutable { + auto C_adj_B_t = (res.adj() * value_of(arena_B).transpose()).eval(); + + if (!is_constant::value) + forward_as>(arena_A).adj() += + value_of(arena_B) * C_adj_B_t; + + if (!is_constant::value) + forward_as>(arena_B).adj() += + value_of(arena_A) * C_adj_B_t.transpose() + + value_of(arena_A).transpose() * value_of(arena_B) * res.adj(); + }); + + return res; +} } // namespace internal /** @@ -157,6 +216,49 @@ inline var quad_form(const EigMat& A, const ColVec& B, bool symmetric = false) { return baseVari->impl_->C_(0, 0); } +/** + * Return the quadratic form \f$ B^T A B \f$. + * + * Symmetry of the resulting matrix is not guaranteed due to numerical + * precision. + * + * @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. + * @throws std::invalid_argument if A is not square, or if A cannot be + * multiplied by B + */ +template * = nullptr, + require_not_col_vector_t* = nullptr, + require_any_var_matrix_t* = nullptr> +inline auto quad_form(const Mat1& A, const Mat2& B, bool symmetric = false) { + return internal::quad_form_impl(A, B, symmetric); +} + +/** + * Return the quadratic form \f$ B^T A B \f$. + * + * @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 * = nullptr, + require_col_vector_t* = nullptr, + require_any_var_matrix_t* = nullptr> +inline var quad_form(const Mat& A, const Vec& B, bool symmetric = false) { + return internal::quad_form_impl(A, B, symmetric)(0, 0); +} + } // namespace math } // namespace stan #endif diff --git a/stan/math/rev/fun/trace_quad_form.hpp b/stan/math/rev/fun/trace_quad_form.hpp index 21b81356e5f..89226b7f4e4 100644 --- a/stan/math/rev/fun/trace_quad_form.hpp +++ b/stan/math/rev/fun/trace_quad_form.hpp @@ -74,7 +74,8 @@ class trace_quad_form_vari : public vari { } // namespace internal template > + require_all_eigen_t* = nullptr, + require_any_st_var* = nullptr> inline return_type_t trace_quad_form(const EigMat1& A, const EigMat2& B) { using Ta = value_type_t; @@ -93,6 +94,33 @@ inline return_type_t trace_quad_form(const EigMat1& A, new internal::trace_quad_form_vari(baseVari)); } +template * = nullptr, + require_any_var_matrix_t* = nullptr> +inline var trace_quad_form(const Mat1& A, const Mat2& B) { + check_square("trace_quad_form", "A", A); + check_multiplicable("trace_quad_form", "A", A, "B", B); + + arena_t arena_A = A; + arena_t arena_B = B; + + var res = (value_of(arena_B).transpose() * + value_of(arena_A) * + value_of(arena_B)).trace(); + + reverse_pass_callback([arena_A, arena_B, res]() mutable { + if (!is_constant::value) + forward_as>(arena_A).adj() += + res.adj() * value_of(arena_B) * value_of(arena_B).transpose(); + + if (!is_constant::value) + forward_as>(arena_B).adj() + += res.adj() * (value_of(arena_A) + value_of(arena_A).transpose()) * value_of(arena_B); + }); + + return res; +} + } // namespace math } // namespace stan #endif diff --git a/test/unit/math/mix/fun/quad_form_test.cpp b/test/unit/math/mix/fun/quad_form_test.cpp index 3d9a9166337..f51d6309604 100644 --- a/test/unit/math/mix/fun/quad_form_test.cpp +++ b/test/unit/math/mix/fun/quad_form_test.cpp @@ -40,13 +40,27 @@ TEST(MathMixMatFun, quadForm) { stan::test::expect_ad(f, a00, a02); stan::test::expect_ad(f, a02, a22); + stan::test::expect_ad_matvar(f, a00, a00); + stan::test::expect_ad_matvar(f, a00, a02); + stan::test::expect_ad_matvar(f, a02, a22); + stan::test::expect_ad(f, a11, b11); stan::test::expect_ad(tols, f, a22, b22); stan::test::expect_ad(f, a22, b23); stan::test::expect_ad(tols, f, a44, b42); + stan::test::expect_ad_matvar(f, a11, b11); + stan::test::expect_ad_matvar(tols, f, a22, b22); + stan::test::expect_ad_matvar(f, a22, b23); + stan::test::expect_ad_matvar(tols, f, a44, b42); + stan::test::expect_ad(f, a00, v0); stan::test::expect_ad(f, a11, v1); stan::test::expect_ad(f, a22, v2); stan::test::expect_ad(tols, f, a44, v4); + + stan::test::expect_ad_matvar(f, a00, v0); + stan::test::expect_ad_matvar(f, a11, v1); + stan::test::expect_ad_matvar(f, a22, v2); + stan::test::expect_ad_matvar(tols, f, a44, v4); } diff --git a/test/unit/math/mix/fun/trace_quad_form_test.cpp b/test/unit/math/mix/fun/trace_quad_form_test.cpp index f070584b3cc..8661911f5bb 100644 --- a/test/unit/math/mix/fun/trace_quad_form_test.cpp +++ b/test/unit/math/mix/fun/trace_quad_form_test.cpp @@ -12,11 +12,17 @@ TEST(MathMixMatFun, traceQuadForm) { stan::test::expect_ad(f, a11, a11); stan::test::expect_ad(f, a11, a1); + stan::test::expect_ad_matvar(f, a11, a11); + stan::test::expect_ad_matvar(f, a11, a1); + Eigen::MatrixXd m00(0, 0); Eigen::VectorXd v0(0); stan::test::expect_ad(f, m00, m00); stan::test::expect_ad(f, m00, v0); + stan::test::expect_ad_matvar(f, m00, m00); + stan::test::expect_ad_matvar(f, m00, v0); + stan::test::ad_tolerances tols; tols.hessian_fvar_hessian_ = 1e0; tols.hessian_hessian_ = 1e0; @@ -31,13 +37,22 @@ TEST(MathMixMatFun, traceQuadForm) { stan::test::expect_ad(tols, f, a22, b22); stan::test::expect_ad(tols, f, a22, a2); + stan::test::expect_ad_matvar(tols, f, a22, a22); + stan::test::expect_ad_matvar(tols, f, a22, b22); + stan::test::expect_ad_matvar(tols, f, a22, a2); + Eigen::MatrixXd a44(4, 4); a44 << 2, 3, 4, 5, 6, 10, 2, 2, 7, 2, 7, 1, 8, 2, 1, 112; Eigen::MatrixXd b42(4, 2); b42 << 100, 10, 0, 1, -3, -3, 5, 2; stan::test::expect_ad(tols, f, a44, b42); + stan::test::expect_ad_matvar(tols, f, a44, b42); + // size mismatch exceptions stan::test::expect_ad(f, a44, b22); stan::test::expect_ad(f, a44, a2); + + stan::test::expect_ad_matvar(f, a44, b22); + stan::test::expect_ad_matvar(f, a44, a2); } diff --git a/test/unit/math/test_ad.hpp b/test/unit/math/test_ad.hpp index 494093ef66a..627844d0a96 100644 --- a/test/unit/math/test_ad.hpp +++ b/test/unit/math/test_ad.hpp @@ -1881,11 +1881,15 @@ template void expect_ad_matvar(const F& f, const EigMat& x) { ad_tolerances tols; expect_ad_matvar_v(tols, f, x); + + stan::math::recover_memory(); } template void expect_ad_matvar(const ad_tolerances& tols, const F& f, const EigMat& x) { expect_ad_matvar_v(tols, f, x); + + stan::math::recover_memory(); } /** @@ -1908,6 +1912,8 @@ void expect_ad_matvar(const ad_tolerances& tols, const F& f, const EigMat1& x, expect_ad_matvar_vv(tols, f, x, y); expect_ad_matvar_vv(tols, f, x, y); expect_ad_matvar_vv(tols, f, x, y); + + stan::math::recover_memory(); } /** @@ -1927,6 +1933,8 @@ template void expect_ad_matvar(const F& f, const EigMat1& x, const EigMat2& y) { ad_tolerances tols; expect_ad_matvar(tols, f, x, y); + + stan::math::recover_memory(); } /** From ae11340c670dddabd8aa11b381fb5883cd891a54 Mon Sep 17 00:00:00 2001 From: Stan Jenkins Date: Mon, 26 Oct 2020 17:38:56 +0000 Subject: [PATCH 02/12] [Jenkins] auto-formatting by clang-format version 6.0.0-1ubuntu2~16.04.1 (tags/RELEASE_600/final) --- stan/math/prim/meta/is_vector.hpp | 12 ++++------- stan/math/rev/fun/quad_form.hpp | 30 +++++++++++++-------------- stan/math/rev/fun/trace_quad_form.hpp | 15 +++++++------- 3 files changed, 26 insertions(+), 31 deletions(-) diff --git a/stan/math/prim/meta/is_vector.hpp b/stan/math/prim/meta/is_vector.hpp index 4a4c38b5dfe..e6e907ef6a9 100644 --- a/stan/math/prim/meta/is_vector.hpp +++ b/stan/math/prim/meta/is_vector.hpp @@ -100,10 +100,8 @@ STAN_ADD_REQUIRE_CONTAINER(eigen_col_vector, is_eigen_col_vector, template struct is_col_vector : internal::is_col_vector_impl {}; -STAN_ADD_REQUIRE_UNARY(col_vector, is_col_vector, - require_eigens_types); -STAN_ADD_REQUIRE_CONTAINER(col_vector, is_col_vector, - require_eigens_types); +STAN_ADD_REQUIRE_UNARY(col_vector, is_col_vector, require_eigens_types); +STAN_ADD_REQUIRE_CONTAINER(col_vector, is_col_vector, require_eigens_types); /** \ingroup type_trait * If the input type T is an eigen matrix with 1 column at compile time this @@ -127,10 +125,8 @@ STAN_ADD_REQUIRE_CONTAINER(eigen_row_vector, is_eigen_row_vector, template struct is_row_vector : internal::is_row_vector_impl {}; -STAN_ADD_REQUIRE_UNARY(row_vector, is_row_vector, - require_eigens_types); -STAN_ADD_REQUIRE_CONTAINER(row_vector, is_row_vector, - require_eigens_types); +STAN_ADD_REQUIRE_UNARY(row_vector, is_row_vector, require_eigens_types); +STAN_ADD_REQUIRE_CONTAINER(row_vector, is_row_vector, require_eigens_types); /** \ingroup type_trait * If the input type T is an eigen matrix with 1 column or 1 row at compile time diff --git a/stan/math/rev/fun/quad_form.hpp b/stan/math/rev/fun/quad_form.hpp index fd4f72a2bf8..9e8093b0e60 100644 --- a/stan/math/rev/fun/quad_form.hpp +++ b/stan/math/rev/fun/quad_form.hpp @@ -117,9 +117,9 @@ inline auto quad_form_impl(const Mat1& A, const Mat2& B, bool symmetric) { check_multiplicable("quad_form", "A", A, "B", B); using return_t - = promote_var_matrix_t; + = promote_var_matrix_t; auto arena_A = to_arena(A); auto arena_B = to_arena(B); @@ -127,27 +127,26 @@ inline auto quad_form_impl(const Mat1& A, const Mat2& B, bool symmetric) { check_not_nan("multiply", "A", value_of(arena_A)); check_not_nan("multiply", "B", value_of(arena_B)); - auto arena_res = to_arena(value_of(arena_B).transpose() * - value_of(arena_A) * - value_of(arena_B)); + auto arena_res = to_arena(value_of(arena_B).transpose() * value_of(arena_A) + * value_of(arena_B)); - if(symmetric) { + if (symmetric) { arena_res = (arena_res + arena_res.transpose()).eval(); } - + return_t res = arena_res; reverse_pass_callback([arena_A, arena_B, res]() mutable { auto C_adj_B_t = (res.adj() * value_of(arena_B).transpose()).eval(); if (!is_constant::value) - forward_as>(arena_A).adj() += - value_of(arena_B) * C_adj_B_t; + forward_as>(arena_A).adj() + += value_of(arena_B) * C_adj_B_t; if (!is_constant::value) - forward_as>(arena_B).adj() += - value_of(arena_A) * C_adj_B_t.transpose() - + value_of(arena_A).transpose() * value_of(arena_B) * res.adj(); + forward_as>(arena_B).adj() + += value_of(arena_A) * C_adj_B_t.transpose() + + value_of(arena_A).transpose() * value_of(arena_B) * res.adj(); }); return res; @@ -251,14 +250,13 @@ inline auto quad_form(const Mat1& A, const Mat2& B, bool symmetric = false) { * @throws std::invalid_argument if A is not square, or if A cannot be * multiplied by B */ -template * = nullptr, +template * = nullptr, require_col_vector_t* = nullptr, require_any_var_matrix_t* = nullptr> inline var quad_form(const Mat& A, const Vec& B, bool symmetric = false) { return internal::quad_form_impl(A, B, symmetric)(0, 0); } - + } // namespace math } // namespace stan #endif diff --git a/stan/math/rev/fun/trace_quad_form.hpp b/stan/math/rev/fun/trace_quad_form.hpp index 89226b7f4e4..46a6636c3c1 100644 --- a/stan/math/rev/fun/trace_quad_form.hpp +++ b/stan/math/rev/fun/trace_quad_form.hpp @@ -74,7 +74,7 @@ class trace_quad_form_vari : public vari { } // namespace internal template * = nullptr, + require_all_eigen_t* = nullptr, require_any_st_var* = nullptr> inline return_type_t trace_quad_form(const EigMat1& A, const EigMat2& B) { @@ -104,18 +104,19 @@ inline var trace_quad_form(const Mat1& A, const Mat2& B) { arena_t arena_A = A; arena_t arena_B = B; - var res = (value_of(arena_B).transpose() * - value_of(arena_A) * - value_of(arena_B)).trace(); + var res + = (value_of(arena_B).transpose() * value_of(arena_A) * value_of(arena_B)) + .trace(); reverse_pass_callback([arena_A, arena_B, res]() mutable { if (!is_constant::value) - forward_as>(arena_A).adj() += - res.adj() * value_of(arena_B) * value_of(arena_B).transpose(); + forward_as>(arena_A).adj() + += res.adj() * value_of(arena_B) * value_of(arena_B).transpose(); if (!is_constant::value) forward_as>(arena_B).adj() - += res.adj() * (value_of(arena_A) + value_of(arena_A).transpose()) * value_of(arena_B); + += res.adj() * (value_of(arena_A) + value_of(arena_A).transpose()) + * value_of(arena_B); }); return res; From e6fb3f233e24a18d7b5a94724a7f151a57f13e7a Mon Sep 17 00:00:00 2001 From: Ben Date: Mon, 26 Oct 2020 13:42:29 -0400 Subject: [PATCH 03/12] Updated quad_form docs (Issue #2101) --- stan/math/rev/fun/quad_form.hpp | 17 ++++++++++------- 1 file changed, 10 insertions(+), 7 deletions(-) diff --git a/stan/math/rev/fun/quad_form.hpp b/stan/math/rev/fun/quad_form.hpp index fd4f72a2bf8..33b39d46473 100644 --- a/stan/math/rev/fun/quad_form.hpp +++ b/stan/math/rev/fun/quad_form.hpp @@ -100,12 +100,13 @@ class quad_form_vari : public vari { * Symmetry of the resulting matrix is not guaranteed due to numerical * precision. * - * @tparam T1 type of the first (square) matrix - * @tparam T2 type of the second matrix + * @tparam Mat1 type of the first (square) matrix + * @tparam Mat2 type of the second matrix * * @param A square matrix * @param B second matrix - * @return The quadratic form, which is a symmetric matrix. + * @param symmetric indicates whether the output should be made symmetric + * @return The quadratic form * @throws std::invalid_argument if A is not square, or if A cannot be * multiplied by B */ @@ -222,11 +223,12 @@ inline var quad_form(const EigMat& A, const ColVec& B, bool symmetric = false) { * Symmetry of the resulting matrix is not guaranteed due to numerical * precision. * - * @tparam EigMat1 type of the first (square) matrix - * @tparam EigMat2 type of the second matrix + * @tparam Mat1 type of the first (square) matrix + * @tparam Mat2 type of the second matrix * * @param A square matrix * @param B second matrix + * @param symmetric indicates whether the output should be made symmetric * @return The quadratic form, which is a symmetric matrix. * @throws std::invalid_argument if A is not square, or if A cannot be * multiplied by B @@ -242,11 +244,12 @@ inline auto quad_form(const Mat1& A, const Mat2& B, bool symmetric = false) { /** * Return the quadratic form \f$ B^T A B \f$. * - * @tparam EigMat type of the matrix - * @tparam ColVec type of the vector + * @tparam Mat type of the matrix + * @tparam Vec type of the vector * * @param A square matrix * @param B vector + * @param symmetric indicates whether the output should be made symmetric * @return The quadratic form (a scalar). * @throws std::invalid_argument if A is not square, or if A cannot be * multiplied by B From cb7b47826d2f0aea358d10fecd13c6b8922eb653 Mon Sep 17 00:00:00 2001 From: Ben Date: Fri, 30 Oct 2020 11:15:08 -0400 Subject: [PATCH 04/12] Convert `mat` to `var` types (Issue #2101) --- stan/math/rev/fun/quad_form.hpp | 11 +++-- stan/math/rev/fun/to_var_value.hpp | 29 ++++++++++++ stan/math/rev/fun/trace_quad_form.hpp | 68 ++++++++++++++++++++++++++- 3 files changed, 102 insertions(+), 6 deletions(-) diff --git a/stan/math/rev/fun/quad_form.hpp b/stan/math/rev/fun/quad_form.hpp index 92a3a5ef3d6..46d066ae787 100644 --- a/stan/math/rev/fun/quad_form.hpp +++ b/stan/math/rev/fun/quad_form.hpp @@ -3,6 +3,7 @@ #include #include +#include #include #include #include @@ -122,9 +123,9 @@ inline auto quad_form_impl(const Mat1& A, const Mat2& B, bool symmetric) { * value_of(A) * value_of(B).eval()), Mat1, Mat2>; - auto arena_A = to_arena(A); - auto arena_B = to_arena(B); - + auto arena_A = to_arena(to_var_value_if::value>(A)); + auto arena_B = to_arena(to_var_value_if::value>(B)); + check_not_nan("multiply", "A", value_of(arena_A)); check_not_nan("multiply", "B", value_of(arena_B)); @@ -141,11 +142,11 @@ inline auto quad_form_impl(const Mat1& A, const Mat2& B, bool symmetric) { auto C_adj_B_t = (res.adj() * value_of(arena_B).transpose()).eval(); if (!is_constant::value) - forward_as>(arena_A).adj() + forward_as>(arena_A).adj() += value_of(arena_B) * C_adj_B_t; if (!is_constant::value) - forward_as>(arena_B).adj() + forward_as>(arena_B).adj() += value_of(arena_A) * C_adj_B_t.transpose() + value_of(arena_A).transpose() * value_of(arena_B) * res.adj(); }); diff --git a/stan/math/rev/fun/to_var_value.hpp b/stan/math/rev/fun/to_var_value.hpp index f1021fa3fc3..0764738308d 100644 --- a/stan/math/rev/fun/to_var_value.hpp +++ b/stan/math/rev/fun/to_var_value.hpp @@ -27,6 +27,35 @@ to_var_value(const T& a) { return res; } +/** + * Forward var_values on. + * + * @tparam T type of the input + * @param a matrix to convert + */ +template * = nullptr> +auto to_var_value(T&& a) { + return std::forward(a); +} + +/** + * If the condition is true, calls `to_var_value` on the argument, otherwise + * forward the argument through. + * + * @tparam T type of argument + * @param a argument + * @return argument copied/evaluated on AD stack + */ +template * = nullptr> +inline auto to_var_value_if(T&& a) { + return std::forward(a); +} + +template * = nullptr> +inline auto to_var_value_if(const T& a) { + return to_var_value(a); +} + } // namespace math } // namespace stan diff --git a/stan/math/rev/fun/trace_quad_form.hpp b/stan/math/rev/fun/trace_quad_form.hpp index 46a6636c3c1..1a116a6f8a4 100644 --- a/stan/math/rev/fun/trace_quad_form.hpp +++ b/stan/math/rev/fun/trace_quad_form.hpp @@ -94,9 +94,75 @@ inline return_type_t trace_quad_form(const EigMat1& A, new internal::trace_quad_form_vari(baseVari)); } +/** + * Compute trace(B^T A B). + * + * This overload handles combinations of + * `Eigen::Matrix, R, C>` and + * `var_value` and + * copies the `Eigen::Matrix, R, C>` argument + * to a `var_value`. + * + * @tparam Mat1 type of the first matrix + * @tparam Mat2 type of the second matrix + * + * @param A matrix + * @param B matrix + * @return The trace of B^T A B + * @throw std::domain_error if A is not square + * @throw std::domain_error if A cannot be multiplied by B + */ +template * = nullptr, + require_eigen_st* = nullptr> +inline var trace_quad_form(const Mat1& A, const Mat2& B) { + return trace_quad_form(A, to_var_value(B)); +} + +/** + * Compute trace(B^T A B). + * + * This overload handles combinations of + * `Eigen::Matrix, R, C>` and + * `var_value` and + * copies the `Eigen::Matrix, R, C>` argument + * to a `var_value`. + * + * @tparam Mat1 type of the first matrix + * @tparam Mat2 type of the second matrix + * + * @param A matrix + * @param B matrix + * @return The trace of B^T A B + * @throw std::domain_error if A is not square + * @throw std::domain_error if A cannot be multiplied by B + */ +template * = nullptr, + require_var_matrix_t* = nullptr> +inline var trace_quad_form(const Mat1& A, const Mat2& B) { + return trace_quad_form(to_var_value(A), B); +} + +/** + * Compute trace(B^T A B). + * + * This overload works with `var_value>` types and + * `Eigen::Matrix` types + * + * @tparam Mat1 type of the first matrix + * @tparam Mat2 type of the second matrix + * + * @param A matrix + * @param B matrix + * @return The trace of B^T A B + * @throw std::domain_error if A is not square + * @throw std::domain_error if A cannot be multiplied by B + */ template * = nullptr, - require_any_var_matrix_t* = nullptr> + require_any_var_matrix_t* = nullptr, + require_all_not_eigen_st* = nullptr> inline var trace_quad_form(const Mat1& A, const Mat2& B) { check_square("trace_quad_form", "A", A); check_multiplicable("trace_quad_form", "A", A, "B", B); From 202b4f4c8bcb75128037c6142423e5cd984dfa65 Mon Sep 17 00:00:00 2001 From: Ben Date: Fri, 30 Oct 2020 11:55:37 -0400 Subject: [PATCH 05/12] Switched trace_quad_form to to_var_value_if and added to_var_value_if tests (Issue #2101) --- stan/math/rev/fun/trace_quad_form.hpp | 64 +--------- test/unit/math/rev/fun/to_var_value_test.cpp | 121 +++++++++++++++++++ 2 files changed, 126 insertions(+), 59 deletions(-) diff --git a/stan/math/rev/fun/trace_quad_form.hpp b/stan/math/rev/fun/trace_quad_form.hpp index 1a116a6f8a4..d953d29e9d4 100644 --- a/stan/math/rev/fun/trace_quad_form.hpp +++ b/stan/math/rev/fun/trace_quad_form.hpp @@ -97,59 +97,6 @@ inline return_type_t trace_quad_form(const EigMat1& A, /** * Compute trace(B^T A B). * - * This overload handles combinations of - * `Eigen::Matrix, R, C>` and - * `var_value` and - * copies the `Eigen::Matrix, R, C>` argument - * to a `var_value`. - * - * @tparam Mat1 type of the first matrix - * @tparam Mat2 type of the second matrix - * - * @param A matrix - * @param B matrix - * @return The trace of B^T A B - * @throw std::domain_error if A is not square - * @throw std::domain_error if A cannot be multiplied by B - */ -template * = nullptr, - require_eigen_st* = nullptr> -inline var trace_quad_form(const Mat1& A, const Mat2& B) { - return trace_quad_form(A, to_var_value(B)); -} - -/** - * Compute trace(B^T A B). - * - * This overload handles combinations of - * `Eigen::Matrix, R, C>` and - * `var_value` and - * copies the `Eigen::Matrix, R, C>` argument - * to a `var_value`. - * - * @tparam Mat1 type of the first matrix - * @tparam Mat2 type of the second matrix - * - * @param A matrix - * @param B matrix - * @return The trace of B^T A B - * @throw std::domain_error if A is not square - * @throw std::domain_error if A cannot be multiplied by B - */ -template * = nullptr, - require_var_matrix_t* = nullptr> -inline var trace_quad_form(const Mat1& A, const Mat2& B) { - return trace_quad_form(to_var_value(A), B); -} - -/** - * Compute trace(B^T A B). - * - * This overload works with `var_value>` types and - * `Eigen::Matrix` types - * * @tparam Mat1 type of the first matrix * @tparam Mat2 type of the second matrix * @@ -161,14 +108,13 @@ inline var trace_quad_form(const Mat1& A, const Mat2& B) { */ template * = nullptr, - require_any_var_matrix_t* = nullptr, - require_all_not_eigen_st* = nullptr> + require_any_var_matrix_t* = nullptr> inline var trace_quad_form(const Mat1& A, const Mat2& B) { check_square("trace_quad_form", "A", A); check_multiplicable("trace_quad_form", "A", A, "B", B); - arena_t arena_A = A; - arena_t arena_B = B; + auto arena_A = to_arena(to_var_value_if::value>(A)); + auto arena_B = to_arena(to_var_value_if::value>(B)); var res = (value_of(arena_B).transpose() * value_of(arena_A) * value_of(arena_B)) @@ -176,11 +122,11 @@ inline var trace_quad_form(const Mat1& A, const Mat2& B) { reverse_pass_callback([arena_A, arena_B, res]() mutable { if (!is_constant::value) - forward_as>(arena_A).adj() + forward_as>(arena_A).adj() += res.adj() * value_of(arena_B) * value_of(arena_B).transpose(); if (!is_constant::value) - forward_as>(arena_B).adj() + forward_as>(arena_B).adj() += res.adj() * (value_of(arena_A) + value_of(arena_A).transpose()) * value_of(arena_B); }); diff --git a/test/unit/math/rev/fun/to_var_value_test.cpp b/test/unit/math/rev/fun/to_var_value_test.cpp index 3b722fc95c2..7afc4e94206 100644 --- a/test/unit/math/rev/fun/to_var_value_test.cpp +++ b/test/unit/math/rev/fun/to_var_value_test.cpp @@ -2,6 +2,127 @@ #include #include +TEST(AgradRevMatrix, to_var_value_types) { + using stan::math::to_var_value; + using stan::math::var; + using stan::math::var_value; + + using mat_var = Eigen::Matrix; + using vec_var = Eigen::Matrix; + using row_vec_var = Eigen::Matrix; + + using var_mat = var_value; + using var_vec = var_value; + using var_row_vec = var_value; + + var a = 2.0; + mat_var b = Eigen::MatrixXd(2, 2); + vec_var c = Eigen::VectorXd(2); + row_vec_var d = Eigen::RowVectorXd(2); + + var_mat e = Eigen::MatrixXd(2, 2); + var_vec f = Eigen::VectorXd(2); + var_row_vec g = Eigen::RowVectorXd(2); + + auto av = to_var_value(a); + auto bv = to_var_value(b); + auto cv = to_var_value(c); + auto dv = to_var_value(d); + + auto ev = to_var_value(e); + auto fv = to_var_value(f); + auto gv = to_var_value(g); + + test::expect_same_type(); + test::expect_same_type(); + test::expect_same_type(); + test::expect_same_type(); + + test::expect_same_type(); + test::expect_same_type(); + test::expect_same_type(); + + stan::math::recover_memory(); +} + +TEST(AgradRevMatrix, to_var_value_if_types) { + using stan::math::to_var_value_if; + using stan::math::var; + using stan::math::var_value; + + using mat_var = Eigen::Matrix; + using vec_var = Eigen::Matrix; + using row_vec_var = Eigen::Matrix; + + using var_mat = var_value; + using var_vec = var_value; + using var_row_vec = var_value; + + var a = 2.0; + mat_var b = Eigen::MatrixXd(2, 2); + vec_var c = Eigen::VectorXd(2); + row_vec_var d = Eigen::RowVectorXd(2); + + var_mat e = Eigen::MatrixXd(2, 2); + var_vec f = Eigen::VectorXd(2); + var_row_vec g = Eigen::RowVectorXd(2); + + double h = 4.0; + Eigen::MatrixXd i(2, 2); + Eigen::VectorXd j(2); + Eigen::RowVectorXd k(2); + + { + auto av = to_var_value_if(a); + auto bv = to_var_value_if(b); + auto cv = to_var_value_if(c); + auto dv = to_var_value_if(d); + + auto ev = to_var_value_if(e); + auto fv = to_var_value_if(f); + auto gv = to_var_value_if(g); + + test::expect_same_type(); + test::expect_same_type(); + test::expect_same_type(); + test::expect_same_type(); + + test::expect_same_type(); + test::expect_same_type(); + test::expect_same_type(); + } + + { + auto av = to_var_value_if(a); + auto bv = to_var_value_if(b); + auto cv = to_var_value_if(c); + auto dv = to_var_value_if(d); + + auto ev = to_var_value_if(e); + auto fv = to_var_value_if(f); + auto gv = to_var_value_if(g); + + auto hv = to_var_value_if(h); + auto iv = to_var_value_if(i); + auto jv = to_var_value_if(j); + + test::expect_same_type(); + test::expect_same_type(); + test::expect_same_type(); + test::expect_same_type(); + + test::expect_same_type(); + test::expect_same_type(); + test::expect_same_type(); + + test::expect_same_type(); + test::expect_same_type(); + test::expect_same_type(); + } + + stan::math::recover_memory(); +} + TEST(AgradRevMatrix, to_var_value_matrix_test) { Eigen::MatrixXd val(2, 3); val << 1, 2, 3, 4, 5, 6; From d5942cf57626f9540e7f1f77a7589241bd8643f7 Mon Sep 17 00:00:00 2001 From: Stan Jenkins Date: Fri, 30 Oct 2020 16:01:15 +0000 Subject: [PATCH 06/12] [Jenkins] auto-formatting by clang-format version 6.0.0-1ubuntu2~16.04.1 (tags/RELEASE_600/final) --- stan/math/rev/fun/quad_form.hpp | 2 +- test/unit/math/rev/fun/to_var_value_test.cpp | 12 ++++++------ 2 files changed, 7 insertions(+), 7 deletions(-) diff --git a/stan/math/rev/fun/quad_form.hpp b/stan/math/rev/fun/quad_form.hpp index 46d066ae787..38361f781c5 100644 --- a/stan/math/rev/fun/quad_form.hpp +++ b/stan/math/rev/fun/quad_form.hpp @@ -125,7 +125,7 @@ inline auto quad_form_impl(const Mat1& A, const Mat2& B, bool symmetric) { auto arena_A = to_arena(to_var_value_if::value>(A)); auto arena_B = to_arena(to_var_value_if::value>(B)); - + check_not_nan("multiply", "A", value_of(arena_A)); check_not_nan("multiply", "B", value_of(arena_B)); diff --git a/test/unit/math/rev/fun/to_var_value_test.cpp b/test/unit/math/rev/fun/to_var_value_test.cpp index 7afc4e94206..9fa2c766944 100644 --- a/test/unit/math/rev/fun/to_var_value_test.cpp +++ b/test/unit/math/rev/fun/to_var_value_test.cpp @@ -6,15 +6,15 @@ TEST(AgradRevMatrix, to_var_value_types) { using stan::math::to_var_value; using stan::math::var; using stan::math::var_value; - + using mat_var = Eigen::Matrix; using vec_var = Eigen::Matrix; using row_vec_var = Eigen::Matrix; - + using var_mat = var_value; using var_vec = var_value; using var_row_vec = var_value; - + var a = 2.0; mat_var b = Eigen::MatrixXd(2, 2); vec_var c = Eigen::VectorXd(2); @@ -32,7 +32,7 @@ TEST(AgradRevMatrix, to_var_value_types) { auto ev = to_var_value(e); auto fv = to_var_value(f); auto gv = to_var_value(g); - + test::expect_same_type(); test::expect_same_type(); test::expect_same_type(); @@ -49,11 +49,11 @@ TEST(AgradRevMatrix, to_var_value_if_types) { using stan::math::to_var_value_if; using stan::math::var; using stan::math::var_value; - + using mat_var = Eigen::Matrix; using vec_var = Eigen::Matrix; using row_vec_var = Eigen::Matrix; - + using var_mat = var_value; using var_vec = var_value; using var_row_vec = var_value; From 1a45c65eaa7c5fb8dd24a75b1ff895638c712a24 Mon Sep 17 00:00:00 2001 From: Ben Date: Mon, 2 Nov 2020 14:09:52 -0500 Subject: [PATCH 07/12] Updated includes (Issue #2101) --- stan/math/rev/fun/trace_quad_form.hpp | 1 + 1 file changed, 1 insertion(+) diff --git a/stan/math/rev/fun/trace_quad_form.hpp b/stan/math/rev/fun/trace_quad_form.hpp index d953d29e9d4..8976d4871b3 100644 --- a/stan/math/rev/fun/trace_quad_form.hpp +++ b/stan/math/rev/fun/trace_quad_form.hpp @@ -4,6 +4,7 @@ #include #include #include +#include #include #include #include From a1eaa0cef5b00a89f9c5104b81fe75bdae43488c Mon Sep 17 00:00:00 2001 From: Ben Date: Thu, 12 Nov 2020 09:12:36 -0500 Subject: [PATCH 08/12] Broke out control flow (Issue #2101) --- stan/math/rev/fun/quad_form.hpp | 89 +++++++++++++++++++++++++-------- 1 file changed, 68 insertions(+), 21 deletions(-) diff --git a/stan/math/rev/fun/quad_form.hpp b/stan/math/rev/fun/quad_form.hpp index 38361f781c5..4019cd4969a 100644 --- a/stan/math/rev/fun/quad_form.hpp +++ b/stan/math/rev/fun/quad_form.hpp @@ -123,35 +123,82 @@ inline auto quad_form_impl(const Mat1& A, const Mat2& B, bool symmetric) { * value_of(A) * value_of(B).eval()), Mat1, Mat2>; - auto arena_A = to_arena(to_var_value_if::value>(A)); - auto arena_B = to_arena(to_var_value_if::value>(B)); + if(!is_constant::value && !is_constant::value) { + auto arena_A = to_arena(to_var_value_if::value>(A)); + auto arena_B = to_arena(to_var_value_if::value>(B)); - check_not_nan("multiply", "A", value_of(arena_A)); - check_not_nan("multiply", "B", value_of(arena_B)); + check_not_nan("multiply", "A", value_of(arena_A)); + check_not_nan("multiply", "B", value_of(arena_B)); - auto arena_res = to_arena(value_of(arena_B).transpose() * value_of(arena_A) - * value_of(arena_B)); + auto arena_res = to_arena(value_of(arena_B).transpose() * value_of(arena_A) + * value_of(arena_B)); - if (symmetric) { - arena_res = (arena_res + arena_res.transpose()).eval(); - } + if (symmetric) { + arena_res += arena_res.transpose().eval(); + } + + return_t res = arena_res; + + reverse_pass_callback([arena_A, arena_B, res]() mutable { + auto C_adj_B_t = (res.adj() * value_of(arena_B).transpose()).eval(); + + forward_as>(arena_A).adj().noalias() += value_of(arena_B) * C_adj_B_t; - return_t res = arena_res; + forward_as>(arena_B).adj().noalias() + += value_of(arena_A) * C_adj_B_t.transpose() + + value_of(arena_A).transpose() * value_of(arena_B) * res.adj(); + }); + + return res; + } else if(!is_constant::value) { + auto arena_A = to_arena(A); + auto arena_B = to_arena(to_var_value_if::value>(B)); + + check_not_nan("multiply", "A", value_of(arena_A)); + check_not_nan("multiply", "B", value_of(arena_B)); + + auto arena_res = to_arena(value_of(arena_B).transpose() * value_of(arena_A) + * value_of(arena_B)); + + if (symmetric) { + arena_res += arena_res.transpose().eval(); + } - reverse_pass_callback([arena_A, arena_B, res]() mutable { - auto C_adj_B_t = (res.adj() * value_of(arena_B).transpose()).eval(); + return_t res = arena_res; - if (!is_constant::value) - forward_as>(arena_A).adj() - += value_of(arena_B) * C_adj_B_t; + reverse_pass_callback([arena_A, arena_B, res]() mutable { + auto C_adj_B_t = (res.adj() * value_of(arena_B).transpose()); - if (!is_constant::value) - forward_as>(arena_B).adj() - += value_of(arena_A) * C_adj_B_t.transpose() - + value_of(arena_A).transpose() * value_of(arena_B) * res.adj(); - }); + forward_as>(arena_B).adj().noalias() + += value_of(arena_A) * C_adj_B_t.transpose() + + value_of(arena_A).transpose() * value_of(arena_B) * res.adj(); + }); - return res; + return res; + } else { + auto arena_A = to_var_value_if::value>(A); + auto arena_B = to_arena(B); + + check_not_nan("multiply", "A", value_of(arena_A)); + check_not_nan("multiply", "B", value_of(arena_B)); + + auto arena_res = to_arena(value_of(arena_B).transpose() * value_of(arena_A) + * value_of(arena_B)); + + if (symmetric) { + arena_res += arena_res.transpose().eval(); + } + + return_t res = arena_res; + + reverse_pass_callback([arena_A, arena_B, res]() mutable { + auto C_adj_B_t = (res.adj() * value_of(arena_B).transpose()); + + forward_as>(arena_A).adj().noalias() += value_of(arena_B) * C_adj_B_t; + }); + + return res; + } } } // namespace internal From 39c5824b92050093fdf1292224c55fb546f47e7d Mon Sep 17 00:00:00 2001 From: Ben Date: Thu, 12 Nov 2020 10:59:41 -0500 Subject: [PATCH 09/12] Broke out more flow control, removed to_var_value_if (Issue #2101) --- stan/math/rev/fun/quad_form.hpp | 67 +++++++++++------ stan/math/rev/fun/to_var_value.hpp | 18 ----- stan/math/rev/fun/trace_quad_form.hpp | 73 ++++++++++++++---- test/unit/math/rev/fun/to_var_value_test.cpp | 78 -------------------- 4 files changed, 101 insertions(+), 135 deletions(-) diff --git a/stan/math/rev/fun/quad_form.hpp b/stan/math/rev/fun/quad_form.hpp index 4019cd4969a..16fbc1c7b5e 100644 --- a/stan/math/rev/fun/quad_form.hpp +++ b/stan/math/rev/fun/quad_form.hpp @@ -119,13 +119,13 @@ inline auto quad_form_impl(const Mat1& A, const Mat2& B, bool symmetric) { check_multiplicable("quad_form", "A", A, "B", B); using return_t - = promote_var_matrix_t; + = promote_var_matrix_t; if(!is_constant::value && !is_constant::value) { - auto arena_A = to_arena(to_var_value_if::value>(A)); - auto arena_B = to_arena(to_var_value_if::value>(B)); + arena_t> arena_A = A; + arena_t> arena_B = B; check_not_nan("multiply", "A", value_of(arena_A)); check_not_nan("multiply", "B", value_of(arena_B)); @@ -142,22 +142,32 @@ inline auto quad_form_impl(const Mat1& A, const Mat2& B, bool symmetric) { reverse_pass_callback([arena_A, arena_B, res]() mutable { auto C_adj_B_t = (res.adj() * value_of(arena_B).transpose()).eval(); - forward_as>(arena_A).adj().noalias() += value_of(arena_B) * C_adj_B_t; + if(is_var_matrix::value) { + arena_A.adj().noalias() += value_of(arena_B) * C_adj_B_t; + } else { + arena_A.adj() += value_of(arena_B) * C_adj_B_t; + } - forward_as>(arena_B).adj().noalias() - += value_of(arena_A) * C_adj_B_t.transpose() - + value_of(arena_A).transpose() * value_of(arena_B) * res.adj(); + if(is_var_matrix::value) { + arena_B.adj().noalias() + += value_of(arena_A) * C_adj_B_t.transpose() + + value_of(arena_A).transpose() * value_of(arena_B) * res.adj(); + } else { + arena_B.adj() + += value_of(arena_A) * C_adj_B_t.transpose() + + value_of(arena_A).transpose() * value_of(arena_B) * res.adj(); + } }); return res; } else if(!is_constant::value) { - auto arena_A = to_arena(A); - auto arena_B = to_arena(to_var_value_if::value>(B)); + arena_t> arena_A = value_of(A); + arena_t> arena_B = B; - check_not_nan("multiply", "A", value_of(arena_A)); - check_not_nan("multiply", "B", value_of(arena_B)); + check_not_nan("multiply", "A", arena_A); + check_not_nan("multiply", "B", arena_B.val()); - auto arena_res = to_arena(value_of(arena_B).transpose() * value_of(arena_A) + auto arena_res = to_arena(value_of(arena_B).transpose() * arena_A * value_of(arena_B)); if (symmetric) { @@ -169,21 +179,26 @@ inline auto quad_form_impl(const Mat1& A, const Mat2& B, bool symmetric) { reverse_pass_callback([arena_A, arena_B, res]() mutable { auto C_adj_B_t = (res.adj() * value_of(arena_B).transpose()); - forward_as>(arena_B).adj().noalias() - += value_of(arena_A) * C_adj_B_t.transpose() - + value_of(arena_A).transpose() * value_of(arena_B) * res.adj(); + if(is_var_matrix::value) { + arena_B.adj().noalias() + += arena_A * C_adj_B_t.transpose() + + arena_A.transpose() * value_of(arena_B) * res.adj(); + } else { + arena_B.adj() + += arena_A * C_adj_B_t.transpose() + + arena_A.transpose() * value_of(arena_B) * res.adj(); + } }); return res; } else { - auto arena_A = to_var_value_if::value>(A); - auto arena_B = to_arena(B); + arena_t> arena_A = A; + arena_t> arena_B = value_of(B); check_not_nan("multiply", "A", value_of(arena_A)); - check_not_nan("multiply", "B", value_of(arena_B)); + check_not_nan("multiply", "B", arena_B); - auto arena_res = to_arena(value_of(arena_B).transpose() * value_of(arena_A) - * value_of(arena_B)); + auto arena_res = to_arena(arena_B.transpose() * value_of(arena_A) * arena_B); if (symmetric) { arena_res += arena_res.transpose().eval(); @@ -192,9 +207,13 @@ inline auto quad_form_impl(const Mat1& A, const Mat2& B, bool symmetric) { return_t res = arena_res; reverse_pass_callback([arena_A, arena_B, res]() mutable { - auto C_adj_B_t = (res.adj() * value_of(arena_B).transpose()); + auto C_adj_B_t = (res.adj() * arena_B.transpose()); - forward_as>(arena_A).adj().noalias() += value_of(arena_B) * C_adj_B_t; + if(is_var_matrix::value) { + arena_A.adj().noalias() += arena_B * C_adj_B_t; + } else { + arena_A.adj() += arena_B * C_adj_B_t; + } }); return res; diff --git a/stan/math/rev/fun/to_var_value.hpp b/stan/math/rev/fun/to_var_value.hpp index 0764738308d..529ada2aa93 100644 --- a/stan/math/rev/fun/to_var_value.hpp +++ b/stan/math/rev/fun/to_var_value.hpp @@ -38,24 +38,6 @@ auto to_var_value(T&& a) { return std::forward(a); } -/** - * If the condition is true, calls `to_var_value` on the argument, otherwise - * forward the argument through. - * - * @tparam T type of argument - * @param a argument - * @return argument copied/evaluated on AD stack - */ -template * = nullptr> -inline auto to_var_value_if(T&& a) { - return std::forward(a); -} - -template * = nullptr> -inline auto to_var_value_if(const T& a) { - return to_var_value(a); -} - } // namespace math } // namespace stan diff --git a/stan/math/rev/fun/trace_quad_form.hpp b/stan/math/rev/fun/trace_quad_form.hpp index 8976d4871b3..4af4d35b249 100644 --- a/stan/math/rev/fun/trace_quad_form.hpp +++ b/stan/math/rev/fun/trace_quad_form.hpp @@ -114,23 +114,66 @@ inline var trace_quad_form(const Mat1& A, const Mat2& B) { check_square("trace_quad_form", "A", A); check_multiplicable("trace_quad_form", "A", A, "B", B); - auto arena_A = to_arena(to_var_value_if::value>(A)); - auto arena_B = to_arena(to_var_value_if::value>(B)); + var res; + + if(!is_constant::value && !is_constant::value) { + arena_t> arena_A = A; + arena_t> arena_B = B; - var res + res = (value_of(arena_B).transpose() * value_of(arena_A) * value_of(arena_B)) - .trace(); - - reverse_pass_callback([arena_A, arena_B, res]() mutable { - if (!is_constant::value) - forward_as>(arena_A).adj() - += res.adj() * value_of(arena_B) * value_of(arena_B).transpose(); - - if (!is_constant::value) - forward_as>(arena_B).adj() - += res.adj() * (value_of(arena_A) + value_of(arena_A).transpose()) - * value_of(arena_B); - }); + .trace(); + + reverse_pass_callback([arena_A, arena_B, res]() mutable { + if(is_var_matrix::value) { + arena_A.adj().noalias() + += res.adj() * value_of(arena_B) * value_of(arena_B).transpose(); + } else { + arena_A.adj() + += res.adj() * value_of(arena_B) * value_of(arena_B).transpose(); + } + + if(is_var_matrix::value) { + arena_B.adj().noalias() + += res.adj() * (value_of(arena_A) + value_of(arena_A).transpose()) + * value_of(arena_B); + } else { + arena_B.adj() + += res.adj() * (value_of(arena_A) + value_of(arena_A).transpose()) + * value_of(arena_B); + } + }); + } else if(!is_constant::value) { + arena_t> arena_A = value_of(A); + arena_t> arena_B = B; + + res + = (value_of(arena_B).transpose() * value_of(arena_A) * value_of(arena_B)) + .trace(); + + reverse_pass_callback([arena_A, arena_B, res]() mutable { + if(is_var_matrix::value) { + arena_B.adj().noalias() + += res.adj() * (arena_A + arena_A.transpose()) * value_of(arena_B); + } else { + arena_B.adj() + += res.adj() * (arena_A + arena_A.transpose()) * value_of(arena_B); + } + }); + } else { + arena_t> arena_A = A; + arena_t> arena_B = value_of(B); + + res = (arena_B.transpose() * value_of(arena_A) * arena_B).trace(); + + reverse_pass_callback([arena_A, arena_B, res]() mutable { + if(is_var_matrix::value) { + arena_A.adj().noalias() += res.adj() * arena_B * arena_B.transpose(); + } else { + arena_A.adj() += res.adj() * arena_B * arena_B.transpose(); + } + }); + } return res; } diff --git a/test/unit/math/rev/fun/to_var_value_test.cpp b/test/unit/math/rev/fun/to_var_value_test.cpp index 9fa2c766944..a940588efd1 100644 --- a/test/unit/math/rev/fun/to_var_value_test.cpp +++ b/test/unit/math/rev/fun/to_var_value_test.cpp @@ -45,84 +45,6 @@ TEST(AgradRevMatrix, to_var_value_types) { stan::math::recover_memory(); } -TEST(AgradRevMatrix, to_var_value_if_types) { - using stan::math::to_var_value_if; - using stan::math::var; - using stan::math::var_value; - - using mat_var = Eigen::Matrix; - using vec_var = Eigen::Matrix; - using row_vec_var = Eigen::Matrix; - - using var_mat = var_value; - using var_vec = var_value; - using var_row_vec = var_value; - - var a = 2.0; - mat_var b = Eigen::MatrixXd(2, 2); - vec_var c = Eigen::VectorXd(2); - row_vec_var d = Eigen::RowVectorXd(2); - - var_mat e = Eigen::MatrixXd(2, 2); - var_vec f = Eigen::VectorXd(2); - var_row_vec g = Eigen::RowVectorXd(2); - - double h = 4.0; - Eigen::MatrixXd i(2, 2); - Eigen::VectorXd j(2); - Eigen::RowVectorXd k(2); - - { - auto av = to_var_value_if(a); - auto bv = to_var_value_if(b); - auto cv = to_var_value_if(c); - auto dv = to_var_value_if(d); - - auto ev = to_var_value_if(e); - auto fv = to_var_value_if(f); - auto gv = to_var_value_if(g); - - test::expect_same_type(); - test::expect_same_type(); - test::expect_same_type(); - test::expect_same_type(); - - test::expect_same_type(); - test::expect_same_type(); - test::expect_same_type(); - } - - { - auto av = to_var_value_if(a); - auto bv = to_var_value_if(b); - auto cv = to_var_value_if(c); - auto dv = to_var_value_if(d); - - auto ev = to_var_value_if(e); - auto fv = to_var_value_if(f); - auto gv = to_var_value_if(g); - - auto hv = to_var_value_if(h); - auto iv = to_var_value_if(i); - auto jv = to_var_value_if(j); - - test::expect_same_type(); - test::expect_same_type(); - test::expect_same_type(); - test::expect_same_type(); - - test::expect_same_type(); - test::expect_same_type(); - test::expect_same_type(); - - test::expect_same_type(); - test::expect_same_type(); - test::expect_same_type(); - } - - stan::math::recover_memory(); -} - TEST(AgradRevMatrix, to_var_value_matrix_test) { Eigen::MatrixXd val(2, 3); val << 1, 2, 3, 4, 5, 6; From 421af8d3e1324e84fec3026b2ff22dee41a9f564 Mon Sep 17 00:00:00 2001 From: Stan Jenkins Date: Thu, 12 Nov 2020 16:23:23 +0000 Subject: [PATCH 10/12] [Jenkins] auto-formatting by clang-format version 6.0.0-1ubuntu2~16.04.1 (tags/RELEASE_600/final) --- stan/math/rev/fun/quad_form.hpp | 58 +++++++++++++-------------- stan/math/rev/fun/trace_quad_form.hpp | 58 +++++++++++++-------------- 2 files changed, 58 insertions(+), 58 deletions(-) diff --git a/stan/math/rev/fun/quad_form.hpp b/stan/math/rev/fun/quad_form.hpp index 16fbc1c7b5e..363286eb42f 100644 --- a/stan/math/rev/fun/quad_form.hpp +++ b/stan/math/rev/fun/quad_form.hpp @@ -119,11 +119,11 @@ inline auto quad_form_impl(const Mat1& A, const Mat2& B, bool symmetric) { check_multiplicable("quad_form", "A", A, "B", B); using return_t - = promote_var_matrix_t; + = promote_var_matrix_t; - if(!is_constant::value && !is_constant::value) { + if (!is_constant::value && !is_constant::value) { arena_t> arena_A = A; arena_t> arena_B = B; @@ -131,7 +131,7 @@ inline auto quad_form_impl(const Mat1& A, const Mat2& B, bool symmetric) { check_not_nan("multiply", "B", value_of(arena_B)); auto arena_res = to_arena(value_of(arena_B).transpose() * value_of(arena_A) - * value_of(arena_B)); + * value_of(arena_B)); if (symmetric) { arena_res += arena_res.transpose().eval(); @@ -142,33 +142,33 @@ inline auto quad_form_impl(const Mat1& A, const Mat2& B, bool symmetric) { reverse_pass_callback([arena_A, arena_B, res]() mutable { auto C_adj_B_t = (res.adj() * value_of(arena_B).transpose()).eval(); - if(is_var_matrix::value) { - arena_A.adj().noalias() += value_of(arena_B) * C_adj_B_t; + if (is_var_matrix::value) { + arena_A.adj().noalias() += value_of(arena_B) * C_adj_B_t; } else { - arena_A.adj() += value_of(arena_B) * C_adj_B_t; + arena_A.adj() += value_of(arena_B) * C_adj_B_t; } - if(is_var_matrix::value) { - arena_B.adj().noalias() - += value_of(arena_A) * C_adj_B_t.transpose() - + value_of(arena_A).transpose() * value_of(arena_B) * res.adj(); + if (is_var_matrix::value) { + arena_B.adj().noalias() + += value_of(arena_A) * C_adj_B_t.transpose() + + value_of(arena_A).transpose() * value_of(arena_B) * res.adj(); } else { - arena_B.adj() - += value_of(arena_A) * C_adj_B_t.transpose() - + value_of(arena_A).transpose() * value_of(arena_B) * res.adj(); + arena_B.adj() + += value_of(arena_A) * C_adj_B_t.transpose() + + value_of(arena_A).transpose() * value_of(arena_B) * res.adj(); } }); return res; - } else if(!is_constant::value) { + } else if (!is_constant::value) { arena_t> arena_A = value_of(A); arena_t> arena_B = B; check_not_nan("multiply", "A", arena_A); check_not_nan("multiply", "B", arena_B.val()); - auto arena_res = to_arena(value_of(arena_B).transpose() * arena_A - * value_of(arena_B)); + auto arena_res + = to_arena(value_of(arena_B).transpose() * arena_A * value_of(arena_B)); if (symmetric) { arena_res += arena_res.transpose().eval(); @@ -179,14 +179,13 @@ inline auto quad_form_impl(const Mat1& A, const Mat2& B, bool symmetric) { reverse_pass_callback([arena_A, arena_B, res]() mutable { auto C_adj_B_t = (res.adj() * value_of(arena_B).transpose()); - if(is_var_matrix::value) { - arena_B.adj().noalias() - += arena_A * C_adj_B_t.transpose() - + arena_A.transpose() * value_of(arena_B) * res.adj(); + if (is_var_matrix::value) { + arena_B.adj().noalias() + += arena_A * C_adj_B_t.transpose() + + arena_A.transpose() * value_of(arena_B) * res.adj(); } else { - arena_B.adj() - += arena_A * C_adj_B_t.transpose() - + arena_A.transpose() * value_of(arena_B) * res.adj(); + arena_B.adj() += arena_A * C_adj_B_t.transpose() + + arena_A.transpose() * value_of(arena_B) * res.adj(); } }); @@ -198,7 +197,8 @@ inline auto quad_form_impl(const Mat1& A, const Mat2& B, bool symmetric) { check_not_nan("multiply", "A", value_of(arena_A)); check_not_nan("multiply", "B", arena_B); - auto arena_res = to_arena(arena_B.transpose() * value_of(arena_A) * arena_B); + auto arena_res + = to_arena(arena_B.transpose() * value_of(arena_A) * arena_B); if (symmetric) { arena_res += arena_res.transpose().eval(); @@ -209,10 +209,10 @@ inline auto quad_form_impl(const Mat1& A, const Mat2& B, bool symmetric) { reverse_pass_callback([arena_A, arena_B, res]() mutable { auto C_adj_B_t = (res.adj() * arena_B.transpose()); - if(is_var_matrix::value) { - arena_A.adj().noalias() += arena_B * C_adj_B_t; + if (is_var_matrix::value) { + arena_A.adj().noalias() += arena_B * C_adj_B_t; } else { - arena_A.adj() += arena_B * C_adj_B_t; + arena_A.adj() += arena_B * C_adj_B_t; } }); diff --git a/stan/math/rev/fun/trace_quad_form.hpp b/stan/math/rev/fun/trace_quad_form.hpp index 4af4d35b249..cd789050b32 100644 --- a/stan/math/rev/fun/trace_quad_form.hpp +++ b/stan/math/rev/fun/trace_quad_form.hpp @@ -115,49 +115,49 @@ inline var trace_quad_form(const Mat1& A, const Mat2& B) { check_multiplicable("trace_quad_form", "A", A, "B", B); var res; - - if(!is_constant::value && !is_constant::value) { + + if (!is_constant::value && !is_constant::value) { arena_t> arena_A = A; arena_t> arena_B = B; - res - = (value_of(arena_B).transpose() * value_of(arena_A) * value_of(arena_B)) - .trace(); + res = (value_of(arena_B).transpose() * value_of(arena_A) + * value_of(arena_B)) + .trace(); reverse_pass_callback([arena_A, arena_B, res]() mutable { - if(is_var_matrix::value) { - arena_A.adj().noalias() - += res.adj() * value_of(arena_B) * value_of(arena_B).transpose(); + if (is_var_matrix::value) { + arena_A.adj().noalias() + += res.adj() * value_of(arena_B) * value_of(arena_B).transpose(); } else { - arena_A.adj() - += res.adj() * value_of(arena_B) * value_of(arena_B).transpose(); + arena_A.adj() + += res.adj() * value_of(arena_B) * value_of(arena_B).transpose(); } - if(is_var_matrix::value) { - arena_B.adj().noalias() - += res.adj() * (value_of(arena_A) + value_of(arena_A).transpose()) - * value_of(arena_B); + if (is_var_matrix::value) { + arena_B.adj().noalias() + += res.adj() * (value_of(arena_A) + value_of(arena_A).transpose()) + * value_of(arena_B); } else { - arena_B.adj() - += res.adj() * (value_of(arena_A) + value_of(arena_A).transpose()) - * value_of(arena_B); + arena_B.adj() += res.adj() + * (value_of(arena_A) + value_of(arena_A).transpose()) + * value_of(arena_B); } }); - } else if(!is_constant::value) { + } else if (!is_constant::value) { arena_t> arena_A = value_of(A); arena_t> arena_B = B; - res - = (value_of(arena_B).transpose() * value_of(arena_A) * value_of(arena_B)) - .trace(); + res = (value_of(arena_B).transpose() * value_of(arena_A) + * value_of(arena_B)) + .trace(); reverse_pass_callback([arena_A, arena_B, res]() mutable { - if(is_var_matrix::value) { - arena_B.adj().noalias() - += res.adj() * (arena_A + arena_A.transpose()) * value_of(arena_B); + if (is_var_matrix::value) { + arena_B.adj().noalias() + += res.adj() * (arena_A + arena_A.transpose()) * value_of(arena_B); } else { - arena_B.adj() - += res.adj() * (arena_A + arena_A.transpose()) * value_of(arena_B); + arena_B.adj() + += res.adj() * (arena_A + arena_A.transpose()) * value_of(arena_B); } }); } else { @@ -167,10 +167,10 @@ inline var trace_quad_form(const Mat1& A, const Mat2& B) { res = (arena_B.transpose() * value_of(arena_A) * arena_B).trace(); reverse_pass_callback([arena_A, arena_B, res]() mutable { - if(is_var_matrix::value) { - arena_A.adj().noalias() += res.adj() * arena_B * arena_B.transpose(); + if (is_var_matrix::value) { + arena_A.adj().noalias() += res.adj() * arena_B * arena_B.transpose(); } else { - arena_A.adj() += res.adj() * arena_B * arena_B.transpose(); + arena_A.adj() += res.adj() * arena_B * arena_B.transpose(); } }); } From e986107fb9db861926a69d1b5223de7e625aed33 Mon Sep 17 00:00:00 2001 From: Ben Date: Thu, 12 Nov 2020 12:33:18 -0500 Subject: [PATCH 11/12] Removed unnecessary deallocate (Issue #2101) --- test/unit/math/test_ad.hpp | 2 -- 1 file changed, 2 deletions(-) diff --git a/test/unit/math/test_ad.hpp b/test/unit/math/test_ad.hpp index 627844d0a96..0f2380ec32b 100644 --- a/test/unit/math/test_ad.hpp +++ b/test/unit/math/test_ad.hpp @@ -1933,8 +1933,6 @@ template void expect_ad_matvar(const F& f, const EigMat1& x, const EigMat2& y) { ad_tolerances tols; expect_ad_matvar(tols, f, x, y); - - stan::math::recover_memory(); } /** From dd8418ac06a5443f56a9c21c5ff2ae7b11469500 Mon Sep 17 00:00:00 2001 From: Ben Date: Mon, 16 Nov 2020 09:29:16 -0500 Subject: [PATCH 12/12] Added more template docs (Issue #2101) --- stan/math/rev/fun/quad_form.hpp | 10 ++++++++++ stan/math/rev/fun/trace_quad_form.hpp | 5 +++++ 2 files changed, 15 insertions(+) diff --git a/stan/math/rev/fun/quad_form.hpp b/stan/math/rev/fun/quad_form.hpp index 363286eb42f..f54d50e2ad1 100644 --- a/stan/math/rev/fun/quad_form.hpp +++ b/stan/math/rev/fun/quad_form.hpp @@ -289,6 +289,11 @@ inline var quad_form(const EigMat& A, const ColVec& B, bool symmetric = false) { * Symmetry of the resulting matrix is not guaranteed due to numerical * precision. * + * This overload handles arguments where one of Mat1 or Mat2 are + * `var_value` where `T` is an Eigen type. The other type can + * also be a `var_value` or it can be a matrix type that inherits + * from EigenBase + * * @tparam Mat1 type of the first (square) matrix * @tparam Mat2 type of the second matrix * @@ -310,6 +315,11 @@ inline auto quad_form(const Mat1& A, const Mat2& B, bool symmetric = false) { /** * Return the quadratic form \f$ B^T A B \f$. * + * This overload handles arguments where one of Mat or Vec are + * `var_value` where `T` is an Eigen type. The other type can + * also be a `var_value`, or it can be a type that inherits from + * EigenBase + * * @tparam Mat type of the matrix * @tparam Vec type of the vector * diff --git a/stan/math/rev/fun/trace_quad_form.hpp b/stan/math/rev/fun/trace_quad_form.hpp index cd789050b32..6b9ceb8ce0a 100644 --- a/stan/math/rev/fun/trace_quad_form.hpp +++ b/stan/math/rev/fun/trace_quad_form.hpp @@ -98,6 +98,11 @@ inline return_type_t trace_quad_form(const EigMat1& A, /** * Compute trace(B^T A B). * + * This overload handles arguments where one of Mat1 or Mat2 are + * `var_value` where `T` is an Eigen type. The other type can + * also be a `var_value` or it can be a type that inherits + * from EigenBase + * * @tparam Mat1 type of the first matrix * @tparam Mat2 type of the second matrix *