From ec9134236271169665425184805f0af981e0750c Mon Sep 17 00:00:00 2001 From: dpsimpson Date: Tue, 31 Mar 2020 02:46:04 -0400 Subject: [PATCH 1/2] Reverse mode specializaiton for matrix add and subtract --- stan/math/prim/fun/add.hpp | 18 +- stan/math/prim/fun/subtract.hpp | 18 +- stan/math/rev/fun.hpp | 2 + stan/math/rev/fun/add.hpp | 548 +++++++++++++++++++++ stan/math/rev/fun/subtract.hpp | 601 +++++++++++++++++++++++ test/unit/math/rev/fun/add_test.cpp | 80 +++ test/unit/math/rev/fun/subtract_test.cpp | 108 ++++ 7 files changed, 1365 insertions(+), 10 deletions(-) create mode 100644 stan/math/rev/fun/add.hpp create mode 100644 stan/math/rev/fun/subtract.hpp create mode 100644 test/unit/math/rev/fun/add_test.cpp create mode 100644 test/unit/math/rev/fun/subtract_test.cpp diff --git a/stan/math/prim/fun/add.hpp b/stan/math/prim/fun/add.hpp index 218fa19ad6b..a58adb5df03 100644 --- a/stan/math/prim/fun/add.hpp +++ b/stan/math/prim/fun/add.hpp @@ -4,6 +4,7 @@ #include #include #include +#include namespace stan { namespace math { @@ -22,7 +23,8 @@ namespace math { * dimensions. */ template > + require_all_eigen_t* = nullptr, + require_all_not_eigen_vt* = nullptr> inline auto add(const Mat1& m1, const Mat2& m2) { check_matching_dims("add", "m1", m1, "m2", m2); return (m1 + m2).eval(); @@ -37,8 +39,11 @@ inline auto add(const Mat1& m1, const Mat2& m2) { * @param c Scalar. * @return The matrix plus the scalar. */ -template , - typename = require_stan_scalar_t> +template * = nullptr, + require_stan_scalar_t* = nullptr, + require_not_var_t* = nullptr, + require_not_eigen_vt* = nullptr> inline auto add(const Mat& m, const Scal c) { return (m.array() + c).matrix().eval(); } @@ -52,8 +57,11 @@ inline auto add(const Mat& m, const Scal c) { * @param m Matrix. * @return The scalar plus the matrix. */ -template , - typename = require_eigen_t> + template * = nullptr, + require_stan_scalar_t* = nullptr, + require_not_var_t* = nullptr, + require_not_eigen_vt* = nullptr> inline auto add(const Scal c, const Mat& m) { return (c + m.array()).matrix().eval(); } diff --git a/stan/math/prim/fun/subtract.hpp b/stan/math/prim/fun/subtract.hpp index 5ba3459d57d..1230a9f06a6 100644 --- a/stan/math/prim/fun/subtract.hpp +++ b/stan/math/prim/fun/subtract.hpp @@ -4,6 +4,7 @@ #include #include #include +#include namespace stan { namespace math { @@ -20,7 +21,8 @@ namespace math { * @return Difference between first matrix and second matrix. */ template > + require_all_eigen_t* = nullptr, + require_all_not_eigen_vt* = nullptr> inline auto subtract(const Mat1& m1, const Mat2& m2) { check_matching_dims("subtract", "m1", m1, "m2", m2); return (m1 - m2).eval(); @@ -36,8 +38,11 @@ inline auto subtract(const Mat1& m1, const Mat2& m2) { * @param m Matrix or expression. * @return The scalar minus the matrix. */ -template , - typename = require_eigen_t> + template * = nullptr, + require_stan_scalar_t* = nullptr, + require_not_var_t* = nullptr, + require_not_eigen_vt* = nullptr> inline auto subtract(const Scal c, const Mat& m) { return (c - m.array()).matrix().eval(); } @@ -52,8 +57,11 @@ inline auto subtract(const Scal c, const Mat& m) { * @param c Scalar. * @return The matrix minus the scalar. */ -template , - typename = require_stan_scalar_t> + template * = nullptr, + require_stan_scalar_t* = nullptr, + require_not_var_t* = nullptr, + require_not_eigen_vt* = nullptr> inline auto subtract(const Mat& m, const Scal c) { return (m.array() - c).matrix().eval(); } diff --git a/stan/math/rev/fun.hpp b/stan/math/rev/fun.hpp index a6025090f21..c0c31b3b5ae 100644 --- a/stan/math/rev/fun.hpp +++ b/stan/math/rev/fun.hpp @@ -12,6 +12,7 @@ #include #include #include +#include #include #include #include @@ -127,6 +128,7 @@ #include #include #include +#include #include #include #include diff --git a/stan/math/rev/fun/add.hpp b/stan/math/rev/fun/add.hpp new file mode 100644 index 00000000000..e46f064e95d --- /dev/null +++ b/stan/math/rev/fun/add.hpp @@ -0,0 +1,548 @@ +#ifndef STAN_MATH_REV_FUN_MATRIX_ADD_HPP +#define STAN_MATH_REV_FUN_MATRIX_ADD_HPP + +#include +#include +#include +#include +#include + +namespace stan { +namespace math { + +/** + * This is a subclass of the vari class for matrix + * addition A + B, where A and B are N by M matrices. + * + * The class stores the structure of each matrix, + * the double values of A and B, and pointers to + * the varis for A and B if A or B is a var. It + * also instatiates and stores pointers to varis + * for all elements of A + B. + * + * @tparam Ta type of elements in matrix A + * @tparam Tb type of elements in matrix B + * @tparam R number of rows (common to A and B) + * @tparam C number of columns (common to A and B) + */ + + template + class add_mat_vari : public vari { + public: + int rows_; + int cols_; + int size_; + double* Ad_; + double* Bd_; + vari** variRefA_; + vari** variRefB_; + vari** variRefAplusB_; + + /** + * Constructor for add_mat_vari. + * + * All memeory is allocated in ChainableStack's + * stack_alloc arena. + * + * It is critical for the efficiency of this object + * that the constructor create new varis that aren't + * popped onto the var_stack_, but rather are + * popped onto the var_nochain_stack_. This is + * controlled by the second argument to + * vari's constructor. + * + * @param A matrix + * @param B matrix + **/ + add_mat_vari(const Eigen::Matrix& A, + const Eigen::Matrix& B) + : vari(0.0), + rows_(A.rows()), + cols_(B.cols()), + size_(A.size()), + Ad_(ChainableStack::instance_->memalloc_.alloc_array(size_)), + Bd_(ChainableStack::instance_->memalloc_.alloc_array(size_)), + variRefA_( + ChainableStack::instance_->memalloc_.alloc_array(size_)), + variRefB_( + ChainableStack::instance_->memalloc_.alloc_array(size_)), + variRefAplusB_( + ChainableStack::instance_->memalloc_.alloc_array(size_)) { + using Eigen::Map; + Map(variRefA_, rows_, cols_) = A.vi(); + Map(variRefB_, rows_, cols_) = B.vi(); + Map Ad(Ad_, rows_, cols_); + Map Bd(Bd_, rows_, cols_); + Ad = A.val(); + Bd = B.val(); + Map(variRefAplusB_, rows_, cols_) + = (Ad + Bd).unaryExpr([](double x) { return new vari(x, false); }); + + } + + virtual void chain() { + using Eigen::Map; + matrix_d adjAplusB(rows_, cols_); + adjAplusB = Map(variRefAplusB_, rows_, cols_).adj(); + Map(variRefA_, rows_, cols_).adj() += adjAplusB; + Map(variRefB_, rows_, cols_).adj() += adjAplusB; + } + }; + + /** + * This is a subclass of the vari class for matrix + * addition A + B, where A is an N by M matrix of + * doubles and and B is an N by M matrix of vars. + * + * The class stores the structure of each matrix, + * the double values of A and B, and pointers to + * the varis for A and B if A or B is a var. It + * also instatiates and stores pointers to varis + * for all elements of A + B. + * + * @tparam Tb type of elements in matrix B + * @tparam R number of rows (common to A and B) + * @tparam C number of columns (common to A and B) + */ + + template + class add_mat_vari : public vari { + public: + int rows_; + int cols_; + int size_; + double* Ad_; + double* Bd_; + vari** variRefB_; + vari** variRefAplusB_; + + /** + * Constructor for add_mat_vari. + * + * All memeory is allocated in ChainableStack's + * stack_alloc arena. + * + * It is critical for the efficiency of this object + * that the constructor create new varis that aren't + * popped onto the var_stack_, but rather are + * popped onto the var_nochain_stack_. This is + * controlled by the second argument to + * vari's constructor. + * + * @param A matrix + * @param B matrix + **/ + add_mat_vari(const Eigen::Matrix& A, + const Eigen::Matrix& B) + : vari(0.0), + rows_(A.rows()), + cols_(B.cols()), + size_(A.size()), + Ad_(ChainableStack::instance_->memalloc_.alloc_array(size_)), + Bd_(ChainableStack::instance_->memalloc_.alloc_array(size_)), + variRefB_( + ChainableStack::instance_->memalloc_.alloc_array(size_)), + variRefAplusB_( + ChainableStack::instance_->memalloc_.alloc_array(size_)) { + using Eigen::Map; + Map(variRefB_, rows_, cols_) = B.vi(); + Map Ad(Ad_, rows_, cols_); + Map Bd(Bd_, rows_, cols_); + Ad = A; + Bd = B.val(); + Map(variRefAplusB_, rows_, cols_) + = (Ad + Bd).unaryExpr([](double x) { return new vari(x, false); }); + + } + + virtual void chain() { + using Eigen::Map; + matrix_d adjAplusB(rows_, cols_); + adjAplusB = Map(variRefAplusB_, rows_, cols_).adj(); + Map(variRefB_, rows_, cols_).adj() += adjAplusB; + } + }; + + /** + * This is a subclass of the vari class for matrix + * addition A + B, where A is an N by M matrix of vars + * and B is an N by M matrix of doubles. + * + * The class stores the structure of each matrix, + * the double values of A and B, and pointers to + * the varis for A and B if A or B is a var. It + * also instatiates and stores pointers to varis + * for all elements of A + B. + * + * @tparam Ta type of elements in matrix A + * @tparam R number of rows (common to A and B) + * @tparam C number of columns (common to A and B) + */ + + template + class add_mat_vari : public vari { + public: + int rows_; + int cols_; + int size_; + double* Ad_; + double* Bd_; + vari** variRefA_; + vari** variRefAplusB_; + + /** + * Constructor for add_mat_vari. + * + * All memeory is allocated in ChainableStack's + * stack_alloc arena. + * + * It is critical for the efficiency of this object + * that the constructor create new varis that aren't + * popped onto the var_stack_, but rather are + * popped onto the var_nochain_stack_. This is + * controlled by the second argument to + * vari's constructor. + * + * @param A matrix + * @param B matrix + **/ + add_mat_vari(const Eigen::Matrix& A, + const Eigen::Matrix& B) + : vari(0.0), + rows_(A.rows()), + cols_(B.cols()), + size_(A.size()), + Ad_(ChainableStack::instance_->memalloc_.alloc_array(size_)), + Bd_(ChainableStack::instance_->memalloc_.alloc_array(size_)), + variRefA_( + ChainableStack::instance_->memalloc_.alloc_array(size_)), + variRefAplusB_( + ChainableStack::instance_->memalloc_.alloc_array(size_)) { + using Eigen::Map; + Map(variRefA_, rows_, cols_) = A.vi(); + Map Ad(Ad_, rows_, cols_); + Map Bd(Bd_, rows_, cols_); + Ad = A.val(); + Bd = B; + Map(variRefAplusB_, rows_, cols_) + = (Ad + Bd).unaryExpr([](double x) { return new vari(x, false); }); + + } + + virtual void chain() { + using Eigen::Map; + matrix_d adjAplusB(rows_, cols_); + adjAplusB = Map(variRefAplusB_, rows_, cols_).adj(); + Map(variRefA_, rows_, cols_).adj() += adjAplusB; + } + }; + + /** + * Return the sum of two matrices + * @tparam Mat1 type of first matrix + * @tparam Mat2 type of second matrix + * + * @param[in] A matrix + * @param[in] B matrix + * @return Sum of two matrices + */ + template * = nullptr, + require_any_eigen_vt* = nullptr> + auto add(const Mat1& A, const Mat2& B) { + using Ta = value_type_t; + using Tb = value_type_t; + constexpr int R = Mat1::RowsAtCompileTime; + constexpr int C = Mat1::ColsAtCompileTime; + check_matching_dims("add", "A", A, "B", B); + check_not_nan("add", "A", A); + check_not_nan("add", "B", B); + + // Manage memory in arena allocator + add_mat_vari* baseVari + = new add_mat_vari(A, B); + Eigen::Matrix AplusB_v(A.rows(), A.cols()); + AplusB_v.vi() = Eigen::Map(&baseVari->variRefAplusB_[0], + A.rows(), A.cols()); + + return AplusB_v; + } + + + /** + * This is a subclass of the vari class for matrix + * addition A + c, where A is a N by M matrix and + * c is a scalar. + * + * The class stores the structure of each matrix, + * the double values of A and c, and pointers to + * the varis for A and c if A or c is a var. It + * also instatiates and stores pointers to varis + * for all elements of A + c. + * + * @tparam Ta type of elements in matrix A + * @tparam Tc type of elements in matrix B + * @tparam R number of rows (common to A and B) + * @tparam C number of columns (common to A and B) + */ + + template + class add_mat_scal_vari : public vari { + public: + int rows_; + int cols_; + int size_; + double* Ad_; + double cd_; + vari** variRefA_; + vari* variRefc_; + vari** variRefAplusc_; + + /** + * Constructor for add_mat_vari. + * + * All memory is allocated in ChainableStack's + * stack_alloc arena. + * + * It is critical for the efficiency of this object + * that the constructor create new varis that aren't + * popped onto the var_stack_, but rather are + * popped onto the var_nochain_stack_. This is + * controlled by the second argument to + * vari's constructor. + * + * @param A matrix + * @param c scalar + **/ + add_mat_scal_vari(const Eigen::Matrix& A, + const Tc c) + : vari(0.0), + rows_(A.rows()), + cols_(A.cols()), + size_(A.size()), + Ad_(ChainableStack::instance_->memalloc_.alloc_array(size_)), + cd_(c.val()), + variRefA_( + ChainableStack::instance_->memalloc_.alloc_array(size_)), + variRefc_(c.vi_), + variRefAplusc_( + ChainableStack::instance_->memalloc_.alloc_array(size_)) { + using Eigen::Map; + Map(variRefA_, rows_, cols_) = A.vi(); + Map Ad(Ad_, rows_, cols_); + double cd = c.val(); + Ad = A.val(); + Map(variRefAplusc_, rows_, cols_) + = (Ad.array() + cd).matrix().eval().unaryExpr([](double x) { return new vari(x, false); }); + + } + + virtual void chain() { + using Eigen::Map; + matrix_d adjAplusc(rows_, cols_); + adjAplusc = Map(variRefAplusc_, rows_, cols_).adj(); + Map(variRefA_, rows_, cols_).adj() += adjAplusc; + variRefc_->adj_ += adjAplusc.sum(); + + } + }; + + /** + * This is a subclass of the vari class for matrix + * addition A + c, where A is a N by M matrix and + * c is a double. + * + * The class stores the structure of each matrix, + * the double values of A and c, and pointers to + * the varis for A and c if A or c is a var. It + * also instatiates and stores pointers to varis + * for all elements of A + c. + * + * @tparam Ta type of elements in matrix A + * @tparam Tc type of elements in matrix B + * @tparam R number of rows (common to A and B) + * @tparam C number of columns (common to A and B) + */ + + template + class add_mat_scal_vari : public vari { + public: + int rows_; + int cols_; + int size_; + double* Ad_; + double cd_; + vari** variRefA_; + vari** variRefAplusc_; + + /** + * Constructor for add_mat_vari. + * + * All memory is allocated in ChainableStack's + * stack_alloc arena. + * + * It is critical for the efficiency of this object + * that the constructor create new varis that aren't + * popped onto the var_stack_, but rather are + * popped onto the var_nochain_stack_. This is + * controlled by the second argument to + * vari's constructor. + * + * @param A matrix + * @param c scalar + **/ + add_mat_scal_vari(const Eigen::Matrix& A, + const double c) + : vari(0.0), + rows_(A.rows()), + cols_(A.cols()), + size_(A.size()), + Ad_(ChainableStack::instance_->memalloc_.alloc_array(size_)), + cd_(c), + variRefA_( + ChainableStack::instance_->memalloc_.alloc_array(size_)), + variRefAplusc_( + ChainableStack::instance_->memalloc_.alloc_array(size_)) { + using Eigen::Map; + Map(variRefA_, rows_, cols_) = A.vi(); + Map Ad(Ad_, rows_, cols_); + double cd = c; + Ad = A.val(); + Map(variRefAplusc_, rows_, cols_) + = (Ad.array() + cd).matrix().eval().unaryExpr([](double x) { return new vari(x, false); }); + + } + + virtual void chain() { + using Eigen::Map; + matrix_d adjAplusc(rows_, cols_); + adjAplusc = Map(variRefAplusc_, rows_, cols_).adj(); + Map(variRefA_, rows_, cols_).adj() += adjAplusc; + } + }; + + + /** + * This is a subclass of the vari class for matrix + * addition A + c, where A is a N by M matrix of + * doubles c is a scalar var. + * + * The class stores the structure of each matrix, + * the double values of A and c, and pointers to + * the varis for A and c if A or c is a var. It + * also instatiates and stores pointers to varis + * for all elements of A + c. + * + * @tparam Ta type of elements in matrix A + * @tparam Tc type of elements in matrix B + * @tparam R number of rows (common to A and B) + * @tparam C number of columns (common to A and B) + */ + + template + class add_mat_scal_vari : public vari { + public: + int rows_; + int cols_; + int size_; + double* Ad_; + double cd_; + vari* variRefc_; + vari** variRefAplusc_; + + /** + * Constructor for add_mat_vari. + * + * All memory is allocated in ChainableStack's + * stack_alloc arena. + * + * It is critical for the efficiency of this object + * that the constructor create new varis that aren't + * popped onto the var_stack_, but rather are + * popped onto the var_nochain_stack_. This is + * controlled by the second argument to + * vari's constructor. + * + * @param A matrix + * @param c scalar + **/ + add_mat_scal_vari(const Eigen::Matrix& A, + const Tc c) + : vari(0.0), + rows_(A.rows()), + cols_(A.cols()), + size_(A.size()), + Ad_(ChainableStack::instance_->memalloc_.alloc_array(size_)), + cd_(c.val()), + variRefc_(c.vi_), + variRefAplusc_( + ChainableStack::instance_->memalloc_.alloc_array(size_)) { + using Eigen::Map; + + double cd = c.val(); + Map Ad(Ad_, rows_, cols_); + Ad = A.val(); + Map(variRefAplusc_, rows_, cols_) + = (Ad.array() + cd).matrix().unaryExpr([](double x) { return new vari(x, false); }); + + } + + virtual void chain() { + using Eigen::Map; + matrix_d adjAplusc(rows_, cols_); + adjAplusc = Map(variRefAplusc_, rows_, cols_).adj(); + variRefc_->adj_ += adjAplusc.sum(); + } + }; + + + /** + * Return the sum of the specified matrix and specified scalar. + * + * @tparam Mat type of the matrix or expression + * @tparam Scal type of the scalar + * @param m Matrix or expression. + * @param c Scalar. + * @return The matrix plus the scalar. + */ + template * = nullptr, + require_stan_scalar_t* = nullptr, + require_any_var_t::type, Scal>* = nullptr> + auto add(const Mat& A, const Scal c) { + using Ta = value_type_t; + constexpr int R = Mat::RowsAtCompileTime; + constexpr int C = Mat::ColsAtCompileTime; + check_not_nan("add matrix-scalar A + c", "A", A); + check_not_nan("add matrix-scalar A + c", "c", c); + + add_mat_scal_vari* baseVari + = new add_mat_scal_vari(A, c); + Eigen::Matrix Aplusc_v(A.rows(), A.cols()); + Aplusc_v.vi() = Eigen::Map( + &baseVari->variRefAplusc_[0], A.rows(), A.cols()); + + return Aplusc_v; + } + + /** + * Return the sum of the specified scalar and specified matrix. + * + * @tparam Scal type of the scalar + * @tparam Mat type of the matrix or expression + * @param c Scalar. + * @param m Matrix. + * @return The scalar plus the matrix. + */ + template * = nullptr, + require_stan_scalar_t* = nullptr, + require_any_var_t::type, Scal>* = nullptr> + inline auto add(const Scal c, const Mat& A) { + return add(A, c); + } + +} // namespace math +} // namespace stan + +#endif diff --git a/stan/math/rev/fun/subtract.hpp b/stan/math/rev/fun/subtract.hpp new file mode 100644 index 00000000000..837f950d946 --- /dev/null +++ b/stan/math/rev/fun/subtract.hpp @@ -0,0 +1,601 @@ +#ifndef STAN_MATH_REV_FUN_MATRIX_SUBTRACT_HPP +#define STAN_MATH_REV_FUN_MATRIX_SUBTRACT_HPP + +#include +#include +#include +#include + +namespace stan { +namespace math { + +/** + * This is a subclass of the vari class for matrix + * subtraction A - B, whre A and B are N by M matrices. + * + * The class stores the structure of each matrix, + * the double values of A and B, and pointers to + * the varis for A and B if A or B is a var. It + * also instatiates and stores pointers to varis + * for all elements of A - B. + * + * @tparam Ta type of elements in matrix A + * @tparam Tb type of elements in matrix B + * @tparam R number of rows (common to A and B) + * @tparam C number of columns (common to A and B) + */ + + template + class subtract_mat_vari : public vari { + public: + int rows_; + int cols_; + int size_; + double* Ad_; + double* Bd_; + vari** variRefA_; + vari** variRefB_; + vari** variRefAminusB_; + + /** + * Constructor for subtract_mat_vari. + * + * All memeory is allocated in ChainableStack's + * stack_alloc arena. + * + * It is critical for the efficiency of this object + * that the constructor create new varis that aren't + * popped onto the var_stack_, but rather are + * popped onto the var_nochain_stack_. This is + * controlled by the second argument to + * vari's constructor. + * + * @param A matrix + * @param B matrix + **/ + subtract_mat_vari(const Eigen::Matrix& A, + const Eigen::Matrix& B) + : vari(0.0), + rows_(A.rows()), + cols_(B.cols()), + size_(A.size()), + Ad_(ChainableStack::instance_->memalloc_.alloc_array(size_)), + Bd_(ChainableStack::instance_->memalloc_.alloc_array(size_)), + variRefA_( + ChainableStack::instance_->memalloc_.alloc_array(size_)), + variRefB_( + ChainableStack::instance_->memalloc_.alloc_array(size_)), + variRefAminusB_( + ChainableStack::instance_->memalloc_.alloc_array(size_)) { + using Eigen::Map; + Map(variRefA_, rows_, cols_) = A.vi(); + Map(variRefB_, rows_, cols_) = B.vi(); + Map Ad(Ad_, rows_, cols_); + Map Bd(Bd_, rows_, cols_); + Ad = A.val(); + Bd = B.val(); + Map(variRefAminusB_, rows_, cols_) + = (Ad - Bd).unaryExpr([](double x) { return new vari(x, false); }); + + } + + virtual void chain() { + using Eigen::Map; + matrix_d adjAminusB(rows_, cols_); + adjAminusB = Map(variRefAminusB_, rows_, cols_).adj(); + Map(variRefA_, rows_, cols_).adj() += adjAminusB; + Map(variRefB_, rows_, cols_).adj() -= adjAminusB; + } + }; + + /** + * This is a subclass of the vari class for matrix + * subtraction A - B, where A is an N by M matrix of + * doubles and and B is an N by M matrix of vars. + * + * The class stores the structure of each matrix, + * the double values of A and B, and pointers to + * the varis for A and B if A or B is a var. It + * also instatiates and stores pointers to varis + * for all elements of A - B. + * + * @tparam Tb type of elements in matrix B + * @tparam R number of rows (common to A and B) + * @tparam C number of columns (common to A and B) + */ + + template + class subtract_mat_vari : public vari { + public: + int rows_; + int cols_; + int size_; + double* Ad_; + double* Bd_; + vari** variRefB_; + vari** variRefAminusB_; + + /** + * Constructor for subtract_mat_vari. + * + * All memeory is allocated in ChainableStack's + * stack_alloc arena. + * + * It is critical for the efficiency of this object + * that the constructor create new varis that aren't + * popped onto the var_stack_, but rather are + * popped onto the var_nochain_stack_. This is + * controlled by the second argument to + * vari's constructor. + * + * @param A matrix + * @param B matrix + **/ + subtract_mat_vari(const Eigen::Matrix& A, + const Eigen::Matrix& B) + : vari(0.0), + rows_(A.rows()), + cols_(B.cols()), + size_(A.size()), + Ad_(ChainableStack::instance_->memalloc_.alloc_array(size_)), + Bd_(ChainableStack::instance_->memalloc_.alloc_array(size_)), + variRefB_( + ChainableStack::instance_->memalloc_.alloc_array(size_)), + variRefAminusB_( + ChainableStack::instance_->memalloc_.alloc_array(size_)) { + using Eigen::Map; + Map(variRefB_, rows_, cols_) = B.vi(); + Map Ad(Ad_, rows_, cols_); + Map Bd(Bd_, rows_, cols_); + Ad = A; + Bd = B.val(); + Map(variRefAminusB_, rows_, cols_) + = (Ad - Bd).unaryExpr([](double x) { return new vari(x, false); }); + + } + + virtual void chain() { + using Eigen::Map; + matrix_d adjAminusB(rows_, cols_); + adjAminusB = Map(variRefAminusB_, rows_, cols_).adj(); + Map(variRefB_, rows_, cols_).adj() -= adjAminusB; + } + }; + + /** + * This is a subclass of the vari class for matrix + * subtraction A - B, whre A is an N by M matrix of vars + * and B is an N by M matrix of doubles. + * + * The class stores the structure of each matrix, + * the double values of A and B, and pointers to + * the varis for A and B if A or B is a var. It + * also instatiates and stores pointers to varis + * for all elements of A - B. + * + * @tparam Ta type of elements in matrix A + * @tparam R number of rows (common to A and B) + * @tparam C number of columns (common to A and B) + */ + + template + class subtract_mat_vari : public vari { + public: + int rows_; + int cols_; + int size_; + double* Ad_; + double* Bd_; + vari** variRefA_; + vari** variRefAminusB_; + + /** + * Constructor for subtract_mat_vari. + * + * All memeory is allocated in ChainableStack's + * stack_alloc arena. + * + * It is critical for the efficiency of this object + * that the constructor create new varis that aren't + * popped onto the var_stack_, but rather are + * popped onto the var_nochain_stack_. This is + * controlled by the second argument to + * vari's constructor. + * + * @param A matrix + * @param B matrix + **/ + subtract_mat_vari(const Eigen::Matrix& A, + const Eigen::Matrix& B) + : vari(0.0), + rows_(A.rows()), + cols_(B.cols()), + size_(A.size()), + Ad_(ChainableStack::instance_->memalloc_.alloc_array(size_)), + Bd_(ChainableStack::instance_->memalloc_.alloc_array(size_)), + variRefA_( + ChainableStack::instance_->memalloc_.alloc_array(size_)), + variRefAminusB_( + ChainableStack::instance_->memalloc_.alloc_array(size_)) { + using Eigen::Map; + Map(variRefA_, rows_, cols_) = A.vi(); + Map Ad(Ad_, rows_, cols_); + Map Bd(Bd_, rows_, cols_); + Ad = A.val(); + Bd = B; + Map(variRefAminusB_, rows_, cols_) + = (Ad - Bd).unaryExpr([](double x) { return new vari(x, false); }); + + } + + virtual void chain() { + using Eigen::Map; + matrix_d adjAminusB(rows_, cols_); + adjAminusB = Map(variRefAminusB_, rows_, cols_).adj(); + Map(variRefA_, rows_, cols_).adj() += adjAminusB; + } + }; + + /** + * Return the sum of two matrices + * @tparam Mat1 type of first matrix + * @tparam Mat2 type of second matrix + * + * @param[in] A matrix + * @param[in] B matrix + * @return Sum of two matrices + */ + template * = nullptr, + require_any_eigen_vt* = nullptr> + auto subtract(const Mat1& A, const Mat2& B) { + using Ta = value_type_t; + using Tb = value_type_t; + constexpr int R = Mat1::RowsAtCompileTime; + constexpr int C = Mat1::ColsAtCompileTime; + check_matching_dims("subtract", "A", A, "B", B); + check_not_nan("subtract", "A", A); + check_not_nan("subtract", "B", B); + + // Manage memory in arena allocator + subtract_mat_vari* baseVari + = new subtract_mat_vari(A, B); + Eigen::Matrix AminusB_v(A.rows(), A.cols()); + AminusB_v.vi() = Eigen::Map(&baseVari->variRefAminusB_[0], + A.rows(), A.cols()); + + return AminusB_v; + } + + /** + * This is a subclass of the vari class for matrix + * subtraction A + c, where A is a N by M matrix and + * c is a scalar. + * + * The class stores the structure of each matrix, + * the double values of A and c, and pointers to + * the varis for A and c if A or c is a var. It + * also instatiates and stores pointers to varis + * for all elements of A + c. + * + * @tparam Ta type of elements in matrix A + * @tparam Tc type of elements in matrix B + * @tparam R number of rows (common to A and B) + * @tparam C number of columns (common to A and B) + * @tparam MatMinusScal: true if computing A - c, false if + * computing c-A + */ + + template + class subtract_mat_scal_vari : public vari { + public: + int rows_; + int cols_; + int size_; + double* Ad_; + double cd_; + vari** variRefA_; + vari* variRefc_; + vari** variRefAminusc_; + + /** + * Constructor for add_mat_vari. + * + * All memory is allocated in ChainableStack's + * stack_alloc arena. + * + * It is critical for the efficiency of this object + * that the constructor create new varis that aren't + * popped onto the var_stack_, but rather are + * popped onto the var_nochain_stack_. This is + * controlled by the second argument to + * vari's constructor. + * + * @param A matrix + * @param c scalar + **/ + subtract_mat_scal_vari(const Eigen::Matrix& A, + const Tc c) + : vari(0.0), + rows_(A.rows()), + cols_(A.cols()), + size_(A.size()), + Ad_(ChainableStack::instance_->memalloc_.alloc_array(size_)), + cd_(c.val()), + variRefA_( + ChainableStack::instance_->memalloc_.alloc_array(size_)), + variRefc_(c.vi_), + variRefAminusc_( + ChainableStack::instance_->memalloc_.alloc_array(size_)) { + using Eigen::Map; + Map(variRefA_, rows_, cols_) = A.vi(); + Map Ad(Ad_, rows_, cols_); + double cd = c.val(); + Ad = A.val(); + if constexpr (MatMinusScal){ + Map(variRefAminusc_, rows_, cols_) + = (Ad.array() - cd).matrix().eval().unaryExpr([](double x) + { return new vari(x, false); }); + } else { + Map(variRefAminusc_, rows_, cols_) + = (cd - Ad.array()).matrix().eval().unaryExpr([](double x) + { return new vari(x, false); }); + } + + } + + virtual void chain() { + using Eigen::Map; + matrix_d adjAminusc(rows_, cols_); + adjAminusc = Map(variRefAminusc_, rows_, cols_).adj(); + if constexpr (MatMinusScal) { + Map(variRefA_, rows_, cols_).adj() += adjAminusc; + variRefc_->adj_ -= adjAminusc.sum(); + } else { + Map(variRefA_, rows_, cols_).adj() -= adjAminusc; + variRefc_->adj_ += adjAminusc.sum(); + } + + } + }; + + /** + * This is a subclass of the vari class for matrix + * subtraction A + c, where A is a N by M matrix and + * c is a double. + * + * The class stores the structure of each matrix, + * the double values of A and c, and pointers to + * the varis for A and c if A or c is a var. It + * also instatiates and stores pointers to varis + * for all elements of A + c. + * + * @tparam Ta type of elements in matrix A + * @tparam Tc type of elements in matrix B + * @tparam R number of rows (common to A and B) + * @tparam C number of columns (common to A and B) + * @tparam MatMinusScal: true if computing A - c, false if + * computing c-A + */ + + template + class subtract_mat_scal_vari + : public vari { + public: + int rows_; + int cols_; + int size_; + double* Ad_; + double cd_; + vari** variRefA_; + vari** variRefAminusc_; + + /** + * Constructor for add_mat_vari. + * + * All memory is allocated in ChainableStack's + * stack_alloc arena. + * + * It is critical for the efficiency of this object + * that the constructor create new varis that aren't + * popped onto the var_stack_, but rather are + * popped onto the var_nochain_stack_. This is + * controlled by the second argument to + * vari's constructor. + * + * @param A matrix + * @param c scalar + **/ + subtract_mat_scal_vari(const Eigen::Matrix& A, + const double c) + : vari(0.0), + rows_(A.rows()), + cols_(A.cols()), + size_(A.size()), + Ad_(ChainableStack::instance_->memalloc_.alloc_array(size_)), + cd_(c), + variRefA_( + ChainableStack::instance_->memalloc_.alloc_array(size_)), + variRefAminusc_( + ChainableStack::instance_->memalloc_.alloc_array(size_)) { + using Eigen::Map; + Map(variRefA_, rows_, cols_) = A.vi(); + Map Ad(Ad_, rows_, cols_); + double cd = c; + Ad = A.val(); + if constexpr (MatMinusScal){ + Map(variRefAminusc_, rows_, cols_) + = (Ad.array() - cd).matrix().eval().unaryExpr([](double x) + { return new vari(x, false); }); + } else { + Map(variRefAminusc_, rows_, cols_) + = (cd - Ad.array()).matrix().eval().unaryExpr([](double x) + { return new vari(x, false); }); + } + + } + + virtual void chain() { + using Eigen::Map; + matrix_d adjAminusc(rows_, cols_); + adjAminusc = Map(variRefAminusc_, rows_, cols_).adj(); + + if constexpr (MatMinusScal) { + Map(variRefA_, rows_, cols_).adj() += adjAminusc; + } else { + Map(variRefA_, rows_, cols_).adj() -= adjAminusc; + } + } + }; + + + /** + * This is a subclass of the vari class for matrix + * subtraction A + c, where A is a N by M matrix of + * doubles c is a scalar var. + * + * The class stores the structure of each matrix, + * the double values of A and c, and pointers to + * the varis for A and c if A or c is a var. It + * also instatiates and stores pointers to varis + * for all elements of A + c. + * + * @tparam Ta type of elements in matrix A + * @tparam Tc type of elements in matrix B + * @tparam R number of rows (common to A and B) + * @tparam C number of columns (common to A and B) + * @tparam MatMinusScal: true if computing A - c, false if + * computing c-A + */ + + template + class subtract_mat_scal_vari + : public vari { + public: + int rows_; + int cols_; + int size_; + double* Ad_; + double cd_; + vari* variRefc_; + vari** variRefAminusc_; + + /** + * Constructor for add_mat_vari. + * + * All memory is allocated in ChainableStack's + * stack_alloc arena. + * + * It is critical for the efficiency of this object + * that the constructor create new varis that aren't + * popped onto the var_stack_, but rather are + * popped onto the var_nochain_stack_. This is + * controlled by the second argument to + * vari's constructor. + * + * @param A matrix + * @param c scalar + **/ + subtract_mat_scal_vari(const Eigen::Matrix& A, + const Tc c) + : vari(0.0), + rows_(A.rows()), + cols_(A.cols()), + size_(A.size()), + Ad_(ChainableStack::instance_->memalloc_.alloc_array(size_)), + cd_(c.val()), + variRefc_(c.vi_), + variRefAminusc_( + ChainableStack::instance_->memalloc_.alloc_array(size_)) { + using Eigen::Map; + + double cd = c.val(); + Map Ad(Ad_, rows_, cols_); + Ad = A.val(); + if constexpr (MatMinusScal){ + Map(variRefAminusc_, rows_, cols_) + = (Ad.array() - cd).matrix().eval().unaryExpr([](double x) + { return new vari(x, false); }); + } else { + Map(variRefAminusc_, rows_, cols_) + = (cd - Ad.array()).matrix().eval().unaryExpr([](double x) + { return new vari(x, false); }); + } + + } + + virtual void chain() { + using Eigen::Map; + matrix_d adjAminusc(rows_, cols_); + adjAminusc = Map(variRefAminusc_, rows_, cols_).adj(); + if constexpr (MatMinusScal) { + variRefc_->adj_ -= adjAminusc.sum(); + } else { + variRefc_->adj_ += adjAminusc.sum(); + } + } + }; + + + /** + * Return the difference A - c of a matrix A and a scalar c. + * + * @tparam Mat type of the matrix or expression + * @tparam Scal type of the scalar + * @param m Matrix or expression. + * @param c Scalar. + * @return The matrix minus the scalar. + */ + template * = nullptr, + require_stan_scalar_t* = nullptr, + require_any_var_t::type, Scal>* = nullptr> + auto subtract(const Mat& A, const Scal c) { + using Ta = value_type_t; + constexpr int R = Mat::RowsAtCompileTime; + constexpr int C = Mat::ColsAtCompileTime; + check_not_nan("subtract matrix-scalar A - c", "A", A); + check_not_nan("subtract matrix-scalar A - c", "c", c); + + subtract_mat_scal_vari* baseVari + = new subtract_mat_scal_vari(A, c); + Eigen::Matrix Aminusc_v(A.rows(), A.cols()); + Aminusc_v.vi() = Eigen::Map( + &baseVari->variRefAminusc_[0], A.rows(), A.cols()); + + return Aminusc_v; + } + + /** + * Return the difference c - A of a scalar c and a matrix A. + * + * @tparam Scal type of the scalar + * @tparam Mat type of the matrix or expression + * @param c Scalar. + * @param m Matrix. + * @return The scalar minus the matrix. + */ + template * = nullptr, + require_stan_scalar_t* = nullptr, + require_any_var_t::type, Scal>* = nullptr> + auto subtract(const Scal c, const Mat& A) { + using Ta = value_type_t; + constexpr int R = Mat::RowsAtCompileTime; + constexpr int C = Mat::ColsAtCompileTime; + check_not_nan("subtract matrix-scalar A - c", "A", A); + check_not_nan("subtract matrix-scalar A - c", "c", c); + + subtract_mat_scal_vari* baseVari + = new subtract_mat_scal_vari(A, c); + Eigen::Matrix Aminusc_v(A.rows(), A.cols()); + Aminusc_v.vi() = Eigen::Map( + &baseVari->variRefAminusc_[0], A.rows(), A.cols()); + + return Aminusc_v; + } + +} // namespace math +} // namespace stan + +#endif diff --git a/test/unit/math/rev/fun/add_test.cpp b/test/unit/math/rev/fun/add_test.cpp new file mode 100644 index 00000000000..1868881f48f --- /dev/null +++ b/test/unit/math/rev/fun/add_test.cpp @@ -0,0 +1,80 @@ +#include +#include + + +TEST(MathRevAdd, add_v_d) { + using stan::math::matrix_d; + using stan::math::matrix_v; + using stan::math::var; + using stan::math::add; + matrix_d Ad = Eigen::MatrixXd::Random(4,6); + matrix_d Bd = Eigen::MatrixXd::Random(4,6); + var c = 3.2; + var d = -2.73; + matrix_v Av = c * Ad; + matrix_v Bv = d * Bd; + matrix_v Cvd = add(Av, Bd); + matrix_v Cdv = add(Ad, Bv); + matrix_v Cvv = add(Av, Bv); + matrix_v Cdvd = add(Ad, + add(Bv, Eigen::MatrixXd::Random(4,6))); + matrix_v Adcv = add(c, Ad); + matrix_v Bvcv = add(c, Bv); + matrix_v Bvscal = add(0.7, Bv); + Cvd(0,1).grad(); + EXPECT_FLOAT_EQ(c.adj(), Ad(0,1)); + EXPECT_FLOAT_EQ(d.adj(), 0.0); + + stan::math::set_zero_all_adjoints(); + Cdv(2,3).grad(); + EXPECT_FLOAT_EQ(d.adj(), Bd(2,3)); + EXPECT_FLOAT_EQ(c.adj(), 0.0); + + stan::math::set_zero_all_adjoints(); + Cvv(3,2).grad(); + EXPECT_FLOAT_EQ(d.adj(), Bd(3,2)); + EXPECT_FLOAT_EQ(c.adj(), Ad(3,2)); + + stan::math::set_zero_all_adjoints(); + Cdvd(3,4).grad(); + EXPECT_FLOAT_EQ(d.adj(), Bd(3,4)); + EXPECT_FLOAT_EQ(c.adj(), 0.0); + + stan::math::set_zero_all_adjoints(); + Adcv(1,1).grad(); + EXPECT_FLOAT_EQ(c.adj(), 1.0); + EXPECT_FLOAT_EQ(d.adj(), 0.0); + + stan::math::set_zero_all_adjoints(); + Bvcv(2,1).grad(); + EXPECT_FLOAT_EQ(c.adj(), 1.0); + EXPECT_FLOAT_EQ(d.adj(), Bd(2,1)); + + stan::math::set_zero_all_adjoints(); + Bvscal(2,0).grad(); + EXPECT_FLOAT_EQ(d.adj(), Bd(2,0)); + EXPECT_FLOAT_EQ(c.adj(), 0.0); +} + +TEST(MathRevAdd, add_check_dim) { + using stan::math::matrix_d; + using stan::math::matrix_v; + using stan::math::add; + matrix_v A, B; + + A.resize(3, 2); + B.resize(3, 2); + EXPECT_NO_THROW(add(A,B)); + + A.resize(1, 0); + B.resize(1, 0); + EXPECT_NO_THROW(add(A,B)); + + A.resize(0, 0); + B.resize(0, 0); + EXPECT_NO_THROW(add(A,B)); + + A.resize(2, 3); + B.resize(3, 2); + EXPECT_THROW(add(A,B), std::invalid_argument); +}; diff --git a/test/unit/math/rev/fun/subtract_test.cpp b/test/unit/math/rev/fun/subtract_test.cpp new file mode 100644 index 00000000000..21ea614459f --- /dev/null +++ b/test/unit/math/rev/fun/subtract_test.cpp @@ -0,0 +1,108 @@ +#include +#include + + +TEST(MathRevSubtract, subtract_v_d) { + using stan::math::matrix_d; + using stan::math::matrix_v; + using stan::math::var; + using stan::math::subtract; + matrix_d Ad = Eigen::MatrixXd::Random(4,6); + matrix_d Bd = Eigen::MatrixXd::Random(4,6); + var c(3.2); + var d(-2.73); + matrix_v Av = c * Ad; + matrix_v Bv = d * Bd; + matrix_v Cvd = subtract(Av, Bd); + matrix_v Cdv = subtract(Ad, Bv); + matrix_v Cvv = subtract(Av, Bv); + matrix_v Cdvd = subtract(Ad, subtract(Bv, Eigen::MatrixXd::Random(4,6))); + matrix_v Cddv = subtract(Ad, subtract(Eigen::MatrixXd::Random(4,6), Bv)); + matrix_v Adcv = subtract(Ad, c); + matrix_v cvAd = subtract(c, Ad); + matrix_v Bvcv = subtract(Bv, c); + matrix_v cvBv = subtract(c, Bv); + matrix_v Bvscal = subtract(Bv, 0.7); + matrix_v scalBv = subtract(0.7, Bv); + + Cvd(0,1).grad(); + EXPECT_FLOAT_EQ(c.adj(), Ad(0,1)); + + stan::math::set_zero_all_adjoints(); + Cdv(2,3).grad(); + EXPECT_FLOAT_EQ(d.adj(), -Bd(2,3)); + + stan::math::set_zero_all_adjoints(); + Cvv(3,2).grad(); + EXPECT_FLOAT_EQ(d.adj(), -Bd(3,2)); + EXPECT_FLOAT_EQ(c.adj(), Ad(3,2)); + + stan::math::set_zero_all_adjoints(); + Cdvd(3,4).grad(); + EXPECT_FLOAT_EQ(d.adj(), -Bd(3,4)); + + stan::math::set_zero_all_adjoints(); + Cddv(2,5).grad(); + EXPECT_FLOAT_EQ(d.adj(), Bd(2,5)); + + + + + + + + stan::math::set_zero_all_adjoints(); + Adcv(1,1).grad(); + EXPECT_FLOAT_EQ(c.adj(), -1.0); + EXPECT_FLOAT_EQ(d.adj(), 0.0); + + stan::math::set_zero_all_adjoints(); + cvAd(1,1).grad(); + EXPECT_FLOAT_EQ(c.adj(), 1.0); + EXPECT_FLOAT_EQ(d.adj(), 0.0); + + stan::math::set_zero_all_adjoints(); + Bvcv(2,1).grad(); + EXPECT_FLOAT_EQ(c.adj(), -1.0); + EXPECT_FLOAT_EQ(d.adj(), Bd(2,1)); + + stan::math::set_zero_all_adjoints(); + cvBv(2,1).grad(); + EXPECT_FLOAT_EQ(c.adj(), 1.0); + EXPECT_FLOAT_EQ(d.adj(), -Bd(2,1)); + + stan::math::set_zero_all_adjoints(); + Bvscal(2,0).grad(); + EXPECT_FLOAT_EQ(d.adj(), Bd(2,0)); + EXPECT_FLOAT_EQ(c.adj(), 0.0); + + stan::math::set_zero_all_adjoints(); + scalBv(2,0).grad(); + EXPECT_FLOAT_EQ(d.adj(), -Bd(2,0)); + EXPECT_FLOAT_EQ(c.adj(), 0.0); +} + + + +TEST(MathRevSubtract, subtract_check_dim) { + using stan::math::matrix_d; + using stan::math::matrix_v; + using stan::math::subtract; + matrix_v A, B; + + A.resize(3, 2); + B.resize(3, 2); + EXPECT_NO_THROW(subtract(A,B)); + + A.resize(1, 0); + B.resize(1, 0); + EXPECT_NO_THROW(subtract(A,B)); + + A.resize(0, 0); + B.resize(0, 0); + EXPECT_NO_THROW(subtract(A,B)); + + A.resize(2, 3); + B.resize(3, 2); + EXPECT_THROW(subtract(A,B), std::invalid_argument); +}; From 7c71f39832bfc92f480be9691d4022b13dff1abd Mon Sep 17 00:00:00 2001 From: Stan Jenkins Date: Tue, 31 Mar 2020 11:04:27 -0400 Subject: [PATCH 2/2] [Jenkins] auto-formatting by clang-format version 6.0.0 (tags/google/stable/2017-11-14) --- stan/math/opencl/kernel_generator/load.hpp | 2 +- stan/math/prim/fun/add.hpp | 16 +- stan/math/prim/fun/subtract.hpp | 18 +- stan/math/rev/fun/add.hpp | 972 +++++++++--------- stan/math/rev/fun/subtract.hpp | 1072 ++++++++++---------- test/unit/math/rev/fun/add_test.cpp | 50 +- test/unit/math/rev/fun/subtract_test.cpp | 69 +- 7 files changed, 1082 insertions(+), 1117 deletions(-) diff --git a/stan/math/opencl/kernel_generator/load.hpp b/stan/math/opencl/kernel_generator/load.hpp index 799cdb551a0..ae0f7d1eb6d 100644 --- a/stan/math/opencl/kernel_generator/load.hpp +++ b/stan/math/opencl/kernel_generator/load.hpp @@ -49,7 +49,7 @@ class load_ * Creates a deep copy of this expression. * @return copy of \c *this */ - inline load_ deep_copy() const & { return load_(a_); } + inline load_ deep_copy() const& { return load_(a_); } inline load_ deep_copy() && { return load_(std::forward(a_)); } /** diff --git a/stan/math/prim/fun/add.hpp b/stan/math/prim/fun/add.hpp index a58adb5df03..d20b0d2fa8a 100644 --- a/stan/math/prim/fun/add.hpp +++ b/stan/math/prim/fun/add.hpp @@ -23,8 +23,8 @@ namespace math { * dimensions. */ template * = nullptr, - require_all_not_eigen_vt* = nullptr> + require_all_eigen_t* = nullptr, + require_all_not_eigen_vt* = nullptr> inline auto add(const Mat1& m1, const Mat2& m2) { check_matching_dims("add", "m1", m1, "m2", m2); return (m1 + m2).eval(); @@ -39,8 +39,7 @@ inline auto add(const Mat1& m1, const Mat2& m2) { * @param c Scalar. * @return The matrix plus the scalar. */ -template * = nullptr, +template * = nullptr, require_stan_scalar_t* = nullptr, require_not_var_t* = nullptr, require_not_eigen_vt* = nullptr> @@ -57,11 +56,10 @@ inline auto add(const Mat& m, const Scal c) { * @param m Matrix. * @return The scalar plus the matrix. */ - template * = nullptr, - require_stan_scalar_t* = nullptr, - require_not_var_t* = nullptr, - require_not_eigen_vt* = nullptr> +template * = nullptr, + require_stan_scalar_t* = nullptr, + require_not_var_t* = nullptr, + require_not_eigen_vt* = nullptr> inline auto add(const Scal c, const Mat& m) { return (c + m.array()).matrix().eval(); } diff --git a/stan/math/prim/fun/subtract.hpp b/stan/math/prim/fun/subtract.hpp index 1230a9f06a6..fe5b44906cd 100644 --- a/stan/math/prim/fun/subtract.hpp +++ b/stan/math/prim/fun/subtract.hpp @@ -38,11 +38,10 @@ inline auto subtract(const Mat1& m1, const Mat2& m2) { * @param m Matrix or expression. * @return The scalar minus the matrix. */ - template * = nullptr, - require_stan_scalar_t* = nullptr, - require_not_var_t* = nullptr, - require_not_eigen_vt* = nullptr> +template * = nullptr, + require_stan_scalar_t* = nullptr, + require_not_var_t* = nullptr, + require_not_eigen_vt* = nullptr> inline auto subtract(const Scal c, const Mat& m) { return (c - m.array()).matrix().eval(); } @@ -57,11 +56,10 @@ inline auto subtract(const Scal c, const Mat& m) { * @param c Scalar. * @return The matrix minus the scalar. */ - template * = nullptr, - require_stan_scalar_t* = nullptr, - require_not_var_t* = nullptr, - require_not_eigen_vt* = nullptr> +template * = nullptr, + require_stan_scalar_t* = nullptr, + require_not_var_t* = nullptr, + require_not_eigen_vt* = nullptr> inline auto subtract(const Mat& m, const Scal c) { return (m.array() - c).matrix().eval(); } diff --git a/stan/math/rev/fun/add.hpp b/stan/math/rev/fun/add.hpp index e46f064e95d..8c892d0dc56 100644 --- a/stan/math/rev/fun/add.hpp +++ b/stan/math/rev/fun/add.hpp @@ -26,19 +26,96 @@ namespace math { * @tparam C number of columns (common to A and B) */ - template - class add_mat_vari : public vari { +template +class add_mat_vari : public vari { public: - int rows_; - int cols_; - int size_; - double* Ad_; - double* Bd_; - vari** variRefA_; - vari** variRefB_; - vari** variRefAplusB_; - - /** + int rows_; + int cols_; + int size_; + double* Ad_; + double* Bd_; + vari** variRefA_; + vari** variRefB_; + vari** variRefAplusB_; + + /** + * Constructor for add_mat_vari. + * + * All memeory is allocated in ChainableStack's + * stack_alloc arena. + * + * It is critical for the efficiency of this object + * that the constructor create new varis that aren't + * popped onto the var_stack_, but rather are + * popped onto the var_nochain_stack_. This is + * controlled by the second argument to + * vari's constructor. + * + * @param A matrix + * @param B matrix + **/ + add_mat_vari(const Eigen::Matrix& A, + const Eigen::Matrix& B) + : vari(0.0), + rows_(A.rows()), + cols_(B.cols()), + size_(A.size()), + Ad_(ChainableStack::instance_->memalloc_.alloc_array(size_)), + Bd_(ChainableStack::instance_->memalloc_.alloc_array(size_)), + variRefA_( + ChainableStack::instance_->memalloc_.alloc_array(size_)), + variRefB_( + ChainableStack::instance_->memalloc_.alloc_array(size_)), + variRefAplusB_( + ChainableStack::instance_->memalloc_.alloc_array(size_)) { + using Eigen::Map; + Map(variRefA_, rows_, cols_) = A.vi(); + Map(variRefB_, rows_, cols_) = B.vi(); + Map Ad(Ad_, rows_, cols_); + Map Bd(Bd_, rows_, cols_); + Ad = A.val(); + Bd = B.val(); + Map(variRefAplusB_, rows_, cols_) + = (Ad + Bd).unaryExpr([](double x) { return new vari(x, false); }); + } + + virtual void chain() { + using Eigen::Map; + matrix_d adjAplusB(rows_, cols_); + adjAplusB = Map(variRefAplusB_, rows_, cols_).adj(); + Map(variRefA_, rows_, cols_).adj() += adjAplusB; + Map(variRefB_, rows_, cols_).adj() += adjAplusB; + } +}; + +/** + * This is a subclass of the vari class for matrix + * addition A + B, where A is an N by M matrix of + * doubles and and B is an N by M matrix of vars. + * + * The class stores the structure of each matrix, + * the double values of A and B, and pointers to + * the varis for A and B if A or B is a var. It + * also instatiates and stores pointers to varis + * for all elements of A + B. + * + * @tparam Tb type of elements in matrix B + * @tparam R number of rows (common to A and B) + * @tparam C number of columns (common to A and B) + */ + +template +class add_mat_vari : public vari { + public: + int rows_; + int cols_; + int size_; + double* Ad_; + double* Bd_; + vari** variRefB_; + vari** variRefAplusB_; + + /** * Constructor for add_mat_vari. * * All memeory is allocated in ChainableStack's @@ -54,493 +131,406 @@ namespace math { * @param A matrix * @param B matrix **/ - add_mat_vari(const Eigen::Matrix& A, - const Eigen::Matrix& B) - : vari(0.0), - rows_(A.rows()), - cols_(B.cols()), - size_(A.size()), - Ad_(ChainableStack::instance_->memalloc_.alloc_array(size_)), - Bd_(ChainableStack::instance_->memalloc_.alloc_array(size_)), - variRefA_( - ChainableStack::instance_->memalloc_.alloc_array(size_)), - variRefB_( - ChainableStack::instance_->memalloc_.alloc_array(size_)), - variRefAplusB_( - ChainableStack::instance_->memalloc_.alloc_array(size_)) { - using Eigen::Map; - Map(variRefA_, rows_, cols_) = A.vi(); - Map(variRefB_, rows_, cols_) = B.vi(); - Map Ad(Ad_, rows_, cols_); - Map Bd(Bd_, rows_, cols_); - Ad = A.val(); - Bd = B.val(); - Map(variRefAplusB_, rows_, cols_) + add_mat_vari(const Eigen::Matrix& A, + const Eigen::Matrix& B) + : vari(0.0), + rows_(A.rows()), + cols_(B.cols()), + size_(A.size()), + Ad_(ChainableStack::instance_->memalloc_.alloc_array(size_)), + Bd_(ChainableStack::instance_->memalloc_.alloc_array(size_)), + variRefB_( + ChainableStack::instance_->memalloc_.alloc_array(size_)), + variRefAplusB_( + ChainableStack::instance_->memalloc_.alloc_array(size_)) { + using Eigen::Map; + Map(variRefB_, rows_, cols_) = B.vi(); + Map Ad(Ad_, rows_, cols_); + Map Bd(Bd_, rows_, cols_); + Ad = A; + Bd = B.val(); + Map(variRefAplusB_, rows_, cols_) = (Ad + Bd).unaryExpr([](double x) { return new vari(x, false); }); + } + + virtual void chain() { + using Eigen::Map; + matrix_d adjAplusB(rows_, cols_); + adjAplusB = Map(variRefAplusB_, rows_, cols_).adj(); + Map(variRefB_, rows_, cols_).adj() += adjAplusB; + } +}; + +/** + * This is a subclass of the vari class for matrix + * addition A + B, where A is an N by M matrix of vars + * and B is an N by M matrix of doubles. + * + * The class stores the structure of each matrix, + * the double values of A and B, and pointers to + * the varis for A and B if A or B is a var. It + * also instatiates and stores pointers to varis + * for all elements of A + B. + * + * @tparam Ta type of elements in matrix A + * @tparam R number of rows (common to A and B) + * @tparam C number of columns (common to A and B) + */ + +template +class add_mat_vari : public vari { + public: + int rows_; + int cols_; + int size_; + double* Ad_; + double* Bd_; + vari** variRefA_; + vari** variRefAplusB_; + + /** + * Constructor for add_mat_vari. + * + * All memeory is allocated in ChainableStack's + * stack_alloc arena. + * + * It is critical for the efficiency of this object + * that the constructor create new varis that aren't + * popped onto the var_stack_, but rather are + * popped onto the var_nochain_stack_. This is + * controlled by the second argument to + * vari's constructor. + * + * @param A matrix + * @param B matrix + **/ + add_mat_vari(const Eigen::Matrix& A, + const Eigen::Matrix& B) + : vari(0.0), + rows_(A.rows()), + cols_(B.cols()), + size_(A.size()), + Ad_(ChainableStack::instance_->memalloc_.alloc_array(size_)), + Bd_(ChainableStack::instance_->memalloc_.alloc_array(size_)), + variRefA_( + ChainableStack::instance_->memalloc_.alloc_array(size_)), + variRefAplusB_( + ChainableStack::instance_->memalloc_.alloc_array(size_)) { + using Eigen::Map; + Map(variRefA_, rows_, cols_) = A.vi(); + Map Ad(Ad_, rows_, cols_); + Map Bd(Bd_, rows_, cols_); + Ad = A.val(); + Bd = B; + Map(variRefAplusB_, rows_, cols_) + = (Ad + Bd).unaryExpr([](double x) { return new vari(x, false); }); + } + + virtual void chain() { + using Eigen::Map; + matrix_d adjAplusB(rows_, cols_); + adjAplusB = Map(variRefAplusB_, rows_, cols_).adj(); + Map(variRefA_, rows_, cols_).adj() += adjAplusB; + } +}; + +/** + * Return the sum of two matrices + * @tparam Mat1 type of first matrix + * @tparam Mat2 type of second matrix + * + * @param[in] A matrix + * @param[in] B matrix + * @return Sum of two matrices + */ +template * = nullptr, + require_any_eigen_vt* = nullptr> +auto add(const Mat1& A, const Mat2& B) { + using Ta = value_type_t; + using Tb = value_type_t; + constexpr int R = Mat1::RowsAtCompileTime; + constexpr int C = Mat1::ColsAtCompileTime; + check_matching_dims("add", "A", A, "B", B); + check_not_nan("add", "A", A); + check_not_nan("add", "B", B); + + // Manage memory in arena allocator + add_mat_vari* baseVari = new add_mat_vari(A, B); + Eigen::Matrix AplusB_v(A.rows(), A.cols()); + AplusB_v.vi() + = Eigen::Map(&baseVari->variRefAplusB_[0], A.rows(), A.cols()); + + return AplusB_v; +} + +/** + * This is a subclass of the vari class for matrix + * addition A + c, where A is a N by M matrix and + * c is a scalar. + * + * The class stores the structure of each matrix, + * the double values of A and c, and pointers to + * the varis for A and c if A or c is a var. It + * also instatiates and stores pointers to varis + * for all elements of A + c. + * + * @tparam Ta type of elements in matrix A + * @tparam Tc type of elements in matrix B + * @tparam R number of rows (common to A and B) + * @tparam C number of columns (common to A and B) + */ - } - - virtual void chain() { - using Eigen::Map; - matrix_d adjAplusB(rows_, cols_); - adjAplusB = Map(variRefAplusB_, rows_, cols_).adj(); - Map(variRefA_, rows_, cols_).adj() += adjAplusB; - Map(variRefB_, rows_, cols_).adj() += adjAplusB; - } - }; - - /** - * This is a subclass of the vari class for matrix - * addition A + B, where A is an N by M matrix of - * doubles and and B is an N by M matrix of vars. - * - * The class stores the structure of each matrix, - * the double values of A and B, and pointers to - * the varis for A and B if A or B is a var. It - * also instatiates and stores pointers to varis - * for all elements of A + B. - * - * @tparam Tb type of elements in matrix B - * @tparam R number of rows (common to A and B) - * @tparam C number of columns (common to A and B) - */ - - template - class add_mat_vari : public vari { - public: - int rows_; - int cols_; - int size_; - double* Ad_; - double* Bd_; - vari** variRefB_; - vari** variRefAplusB_; - - /** - * Constructor for add_mat_vari. - * - * All memeory is allocated in ChainableStack's - * stack_alloc arena. - * - * It is critical for the efficiency of this object - * that the constructor create new varis that aren't - * popped onto the var_stack_, but rather are - * popped onto the var_nochain_stack_. This is - * controlled by the second argument to - * vari's constructor. - * - * @param A matrix - * @param B matrix - **/ - add_mat_vari(const Eigen::Matrix& A, - const Eigen::Matrix& B) - : vari(0.0), - rows_(A.rows()), - cols_(B.cols()), - size_(A.size()), - Ad_(ChainableStack::instance_->memalloc_.alloc_array(size_)), - Bd_(ChainableStack::instance_->memalloc_.alloc_array(size_)), - variRefB_( +template +class add_mat_scal_vari : public vari { + public: + int rows_; + int cols_; + int size_; + double* Ad_; + double cd_; + vari** variRefA_; + vari* variRefc_; + vari** variRefAplusc_; + + /** + * Constructor for add_mat_vari. + * + * All memory is allocated in ChainableStack's + * stack_alloc arena. + * + * It is critical for the efficiency of this object + * that the constructor create new varis that aren't + * popped onto the var_stack_, but rather are + * popped onto the var_nochain_stack_. This is + * controlled by the second argument to + * vari's constructor. + * + * @param A matrix + * @param c scalar + **/ + add_mat_scal_vari(const Eigen::Matrix& A, const Tc c) + : vari(0.0), + rows_(A.rows()), + cols_(A.cols()), + size_(A.size()), + Ad_(ChainableStack::instance_->memalloc_.alloc_array(size_)), + cd_(c.val()), + variRefA_( ChainableStack::instance_->memalloc_.alloc_array(size_)), - variRefAplusB_( + variRefc_(c.vi_), + variRefAplusc_( ChainableStack::instance_->memalloc_.alloc_array(size_)) { - using Eigen::Map; - Map(variRefB_, rows_, cols_) = B.vi(); - Map Ad(Ad_, rows_, cols_); - Map Bd(Bd_, rows_, cols_); - Ad = A; - Bd = B.val(); - Map(variRefAplusB_, rows_, cols_) - = (Ad + Bd).unaryExpr([](double x) { return new vari(x, false); }); - - } - - virtual void chain() { - using Eigen::Map; - matrix_d adjAplusB(rows_, cols_); - adjAplusB = Map(variRefAplusB_, rows_, cols_).adj(); - Map(variRefB_, rows_, cols_).adj() += adjAplusB; - } - }; + using Eigen::Map; + Map(variRefA_, rows_, cols_) = A.vi(); + Map Ad(Ad_, rows_, cols_); + double cd = c.val(); + Ad = A.val(); + Map(variRefAplusc_, rows_, cols_) + = (Ad.array() + cd).matrix().eval().unaryExpr([](double x) { + return new vari(x, false); + }); + } + + virtual void chain() { + using Eigen::Map; + matrix_d adjAplusc(rows_, cols_); + adjAplusc = Map(variRefAplusc_, rows_, cols_).adj(); + Map(variRefA_, rows_, cols_).adj() += adjAplusc; + variRefc_->adj_ += adjAplusc.sum(); + } +}; + +/** + * This is a subclass of the vari class for matrix + * addition A + c, where A is a N by M matrix and + * c is a double. + * + * The class stores the structure of each matrix, + * the double values of A and c, and pointers to + * the varis for A and c if A or c is a var. It + * also instatiates and stores pointers to varis + * for all elements of A + c. + * + * @tparam Ta type of elements in matrix A + * @tparam Tc type of elements in matrix B + * @tparam R number of rows (common to A and B) + * @tparam C number of columns (common to A and B) + */ + +template +class add_mat_scal_vari : public vari { + public: + int rows_; + int cols_; + int size_; + double* Ad_; + double cd_; + vari** variRefA_; + vari** variRefAplusc_; /** - * This is a subclass of the vari class for matrix - * addition A + B, where A is an N by M matrix of vars - * and B is an N by M matrix of doubles. + * Constructor for add_mat_vari. * - * The class stores the structure of each matrix, - * the double values of A and B, and pointers to - * the varis for A and B if A or B is a var. It - * also instatiates and stores pointers to varis - * for all elements of A + B. + * All memory is allocated in ChainableStack's + * stack_alloc arena. * - * @tparam Ta type of elements in matrix A - * @tparam R number of rows (common to A and B) - * @tparam C number of columns (common to A and B) - */ - - template - class add_mat_vari : public vari { - public: - int rows_; - int cols_; - int size_; - double* Ad_; - double* Bd_; - vari** variRefA_; - vari** variRefAplusB_; - - /** - * Constructor for add_mat_vari. - * - * All memeory is allocated in ChainableStack's - * stack_alloc arena. - * - * It is critical for the efficiency of this object - * that the constructor create new varis that aren't - * popped onto the var_stack_, but rather are - * popped onto the var_nochain_stack_. This is - * controlled by the second argument to - * vari's constructor. - * - * @param A matrix - * @param B matrix - **/ - add_mat_vari(const Eigen::Matrix& A, - const Eigen::Matrix& B) - : vari(0.0), - rows_(A.rows()), - cols_(B.cols()), - size_(A.size()), - Ad_(ChainableStack::instance_->memalloc_.alloc_array(size_)), - Bd_(ChainableStack::instance_->memalloc_.alloc_array(size_)), - variRefA_( - ChainableStack::instance_->memalloc_.alloc_array(size_)), - variRefAplusB_( - ChainableStack::instance_->memalloc_.alloc_array(size_)) { - using Eigen::Map; - Map(variRefA_, rows_, cols_) = A.vi(); - Map Ad(Ad_, rows_, cols_); - Map Bd(Bd_, rows_, cols_); - Ad = A.val(); - Bd = B; - Map(variRefAplusB_, rows_, cols_) - = (Ad + Bd).unaryExpr([](double x) { return new vari(x, false); }); - - } - - virtual void chain() { - using Eigen::Map; - matrix_d adjAplusB(rows_, cols_); - adjAplusB = Map(variRefAplusB_, rows_, cols_).adj(); - Map(variRefA_, rows_, cols_).adj() += adjAplusB; - } - }; - - /** - * Return the sum of two matrices - * @tparam Mat1 type of first matrix - * @tparam Mat2 type of second matrix - * - * @param[in] A matrix - * @param[in] B matrix - * @return Sum of two matrices - */ - template * = nullptr, - require_any_eigen_vt* = nullptr> - auto add(const Mat1& A, const Mat2& B) { - using Ta = value_type_t; - using Tb = value_type_t; - constexpr int R = Mat1::RowsAtCompileTime; - constexpr int C = Mat1::ColsAtCompileTime; - check_matching_dims("add", "A", A, "B", B); - check_not_nan("add", "A", A); - check_not_nan("add", "B", B); - - // Manage memory in arena allocator - add_mat_vari* baseVari - = new add_mat_vari(A, B); - Eigen::Matrix AplusB_v(A.rows(), A.cols()); - AplusB_v.vi() = Eigen::Map(&baseVari->variRefAplusB_[0], - A.rows(), A.cols()); - - return AplusB_v; - } - - - /** - * This is a subclass of the vari class for matrix - * addition A + c, where A is a N by M matrix and - * c is a scalar. - * - * The class stores the structure of each matrix, - * the double values of A and c, and pointers to - * the varis for A and c if A or c is a var. It - * also instatiates and stores pointers to varis - * for all elements of A + c. - * - * @tparam Ta type of elements in matrix A - * @tparam Tc type of elements in matrix B - * @tparam R number of rows (common to A and B) - * @tparam C number of columns (common to A and B) - */ - - template - class add_mat_scal_vari : public vari { - public: - int rows_; - int cols_; - int size_; - double* Ad_; - double cd_; - vari** variRefA_; - vari* variRefc_; - vari** variRefAplusc_; - - /** - * Constructor for add_mat_vari. - * - * All memory is allocated in ChainableStack's - * stack_alloc arena. - * - * It is critical for the efficiency of this object - * that the constructor create new varis that aren't - * popped onto the var_stack_, but rather are - * popped onto the var_nochain_stack_. This is - * controlled by the second argument to - * vari's constructor. - * - * @param A matrix - * @param c scalar - **/ - add_mat_scal_vari(const Eigen::Matrix& A, - const Tc c) - : vari(0.0), - rows_(A.rows()), - cols_(A.cols()), - size_(A.size()), - Ad_(ChainableStack::instance_->memalloc_.alloc_array(size_)), - cd_(c.val()), - variRefA_( - ChainableStack::instance_->memalloc_.alloc_array(size_)), - variRefc_(c.vi_), - variRefAplusc_( - ChainableStack::instance_->memalloc_.alloc_array(size_)) { - using Eigen::Map; - Map(variRefA_, rows_, cols_) = A.vi(); - Map Ad(Ad_, rows_, cols_); - double cd = c.val(); - Ad = A.val(); - Map(variRefAplusc_, rows_, cols_) - = (Ad.array() + cd).matrix().eval().unaryExpr([](double x) { return new vari(x, false); }); - - } - - virtual void chain() { - using Eigen::Map; - matrix_d adjAplusc(rows_, cols_); - adjAplusc = Map(variRefAplusc_, rows_, cols_).adj(); - Map(variRefA_, rows_, cols_).adj() += adjAplusc; - variRefc_->adj_ += adjAplusc.sum(); - - } - }; - - /** - * This is a subclass of the vari class for matrix - * addition A + c, where A is a N by M matrix and - * c is a double. - * - * The class stores the structure of each matrix, - * the double values of A and c, and pointers to - * the varis for A and c if A or c is a var. It - * also instatiates and stores pointers to varis - * for all elements of A + c. - * - * @tparam Ta type of elements in matrix A - * @tparam Tc type of elements in matrix B - * @tparam R number of rows (common to A and B) - * @tparam C number of columns (common to A and B) - */ - - template - class add_mat_scal_vari : public vari { - public: - int rows_; - int cols_; - int size_; - double* Ad_; - double cd_; - vari** variRefA_; - vari** variRefAplusc_; - - /** - * Constructor for add_mat_vari. - * - * All memory is allocated in ChainableStack's - * stack_alloc arena. - * - * It is critical for the efficiency of this object - * that the constructor create new varis that aren't - * popped onto the var_stack_, but rather are - * popped onto the var_nochain_stack_. This is - * controlled by the second argument to - * vari's constructor. - * - * @param A matrix - * @param c scalar - **/ - add_mat_scal_vari(const Eigen::Matrix& A, - const double c) - : vari(0.0), - rows_(A.rows()), - cols_(A.cols()), - size_(A.size()), - Ad_(ChainableStack::instance_->memalloc_.alloc_array(size_)), - cd_(c), - variRefA_( - ChainableStack::instance_->memalloc_.alloc_array(size_)), - variRefAplusc_( - ChainableStack::instance_->memalloc_.alloc_array(size_)) { - using Eigen::Map; - Map(variRefA_, rows_, cols_) = A.vi(); - Map Ad(Ad_, rows_, cols_); - double cd = c; - Ad = A.val(); - Map(variRefAplusc_, rows_, cols_) - = (Ad.array() + cd).matrix().eval().unaryExpr([](double x) { return new vari(x, false); }); - - } - - virtual void chain() { - using Eigen::Map; - matrix_d adjAplusc(rows_, cols_); - adjAplusc = Map(variRefAplusc_, rows_, cols_).adj(); - Map(variRefA_, rows_, cols_).adj() += adjAplusc; - } - }; - - - /** - * This is a subclass of the vari class for matrix - * addition A + c, where A is a N by M matrix of - * doubles c is a scalar var. - * - * The class stores the structure of each matrix, - * the double values of A and c, and pointers to - * the varis for A and c if A or c is a var. It - * also instatiates and stores pointers to varis - * for all elements of A + c. - * - * @tparam Ta type of elements in matrix A - * @tparam Tc type of elements in matrix B - * @tparam R number of rows (common to A and B) - * @tparam C number of columns (common to A and B) - */ - - template - class add_mat_scal_vari : public vari { - public: - int rows_; - int cols_; - int size_; - double* Ad_; - double cd_; - vari* variRefc_; - vari** variRefAplusc_; - - /** - * Constructor for add_mat_vari. - * - * All memory is allocated in ChainableStack's - * stack_alloc arena. - * - * It is critical for the efficiency of this object - * that the constructor create new varis that aren't - * popped onto the var_stack_, but rather are - * popped onto the var_nochain_stack_. This is - * controlled by the second argument to - * vari's constructor. - * - * @param A matrix - * @param c scalar - **/ - add_mat_scal_vari(const Eigen::Matrix& A, - const Tc c) - : vari(0.0), - rows_(A.rows()), - cols_(A.cols()), - size_(A.size()), - Ad_(ChainableStack::instance_->memalloc_.alloc_array(size_)), - cd_(c.val()), - variRefc_(c.vi_), - variRefAplusc_( - ChainableStack::instance_->memalloc_.alloc_array(size_)) { - using Eigen::Map; - - double cd = c.val(); - Map Ad(Ad_, rows_, cols_); - Ad = A.val(); - Map(variRefAplusc_, rows_, cols_) - = (Ad.array() + cd).matrix().unaryExpr([](double x) { return new vari(x, false); }); - - } - - virtual void chain() { - using Eigen::Map; - matrix_d adjAplusc(rows_, cols_); - adjAplusc = Map(variRefAplusc_, rows_, cols_).adj(); - variRefc_->adj_ += adjAplusc.sum(); - } - }; - - - /** - * Return the sum of the specified matrix and specified scalar. - * - * @tparam Mat type of the matrix or expression - * @tparam Scal type of the scalar - * @param m Matrix or expression. - * @param c Scalar. - * @return The matrix plus the scalar. - */ - template * = nullptr, - require_stan_scalar_t* = nullptr, - require_any_var_t::type, Scal>* = nullptr> - auto add(const Mat& A, const Scal c) { - using Ta = value_type_t; - constexpr int R = Mat::RowsAtCompileTime; - constexpr int C = Mat::ColsAtCompileTime; - check_not_nan("add matrix-scalar A + c", "A", A); - check_not_nan("add matrix-scalar A + c", "c", c); - - add_mat_scal_vari* baseVari - = new add_mat_scal_vari(A, c); - Eigen::Matrix Aplusc_v(A.rows(), A.cols()); - Aplusc_v.vi() = Eigen::Map( - &baseVari->variRefAplusc_[0], A.rows(), A.cols()); - - return Aplusc_v; - } - - /** - * Return the sum of the specified scalar and specified matrix. - * - * @tparam Scal type of the scalar - * @tparam Mat type of the matrix or expression - * @param c Scalar. - * @param m Matrix. - * @return The scalar plus the matrix. - */ - template * = nullptr, - require_stan_scalar_t* = nullptr, - require_any_var_t::type, Scal>* = nullptr> - inline auto add(const Scal c, const Mat& A) { - return add(A, c); - } + * It is critical for the efficiency of this object + * that the constructor create new varis that aren't + * popped onto the var_stack_, but rather are + * popped onto the var_nochain_stack_. This is + * controlled by the second argument to + * vari's constructor. + * + * @param A matrix + * @param c scalar + **/ + add_mat_scal_vari(const Eigen::Matrix& A, const double c) + : vari(0.0), + rows_(A.rows()), + cols_(A.cols()), + size_(A.size()), + Ad_(ChainableStack::instance_->memalloc_.alloc_array(size_)), + cd_(c), + variRefA_( + ChainableStack::instance_->memalloc_.alloc_array(size_)), + variRefAplusc_( + ChainableStack::instance_->memalloc_.alloc_array(size_)) { + using Eigen::Map; + Map(variRefA_, rows_, cols_) = A.vi(); + Map Ad(Ad_, rows_, cols_); + double cd = c; + Ad = A.val(); + Map(variRefAplusc_, rows_, cols_) + = (Ad.array() + cd).matrix().eval().unaryExpr([](double x) { + return new vari(x, false); + }); + } + + virtual void chain() { + using Eigen::Map; + matrix_d adjAplusc(rows_, cols_); + adjAplusc = Map(variRefAplusc_, rows_, cols_).adj(); + Map(variRefA_, rows_, cols_).adj() += adjAplusc; + } +}; + +/** + * This is a subclass of the vari class for matrix + * addition A + c, where A is a N by M matrix of + * doubles c is a scalar var. + * + * The class stores the structure of each matrix, + * the double values of A and c, and pointers to + * the varis for A and c if A or c is a var. It + * also instatiates and stores pointers to varis + * for all elements of A + c. + * + * @tparam Ta type of elements in matrix A + * @tparam Tc type of elements in matrix B + * @tparam R number of rows (common to A and B) + * @tparam C number of columns (common to A and B) + */ + +template +class add_mat_scal_vari : public vari { + public: + int rows_; + int cols_; + int size_; + double* Ad_; + double cd_; + vari* variRefc_; + vari** variRefAplusc_; + + /** + * Constructor for add_mat_vari. + * + * All memory is allocated in ChainableStack's + * stack_alloc arena. + * + * It is critical for the efficiency of this object + * that the constructor create new varis that aren't + * popped onto the var_stack_, but rather are + * popped onto the var_nochain_stack_. This is + * controlled by the second argument to + * vari's constructor. + * + * @param A matrix + * @param c scalar + **/ + add_mat_scal_vari(const Eigen::Matrix& A, const Tc c) + : vari(0.0), + rows_(A.rows()), + cols_(A.cols()), + size_(A.size()), + Ad_(ChainableStack::instance_->memalloc_.alloc_array(size_)), + cd_(c.val()), + variRefc_(c.vi_), + variRefAplusc_( + ChainableStack::instance_->memalloc_.alloc_array(size_)) { + using Eigen::Map; + + double cd = c.val(); + Map Ad(Ad_, rows_, cols_); + Ad = A.val(); + Map(variRefAplusc_, rows_, cols_) + = (Ad.array() + cd).matrix().unaryExpr([](double x) { + return new vari(x, false); + }); + } + + virtual void chain() { + using Eigen::Map; + matrix_d adjAplusc(rows_, cols_); + adjAplusc = Map(variRefAplusc_, rows_, cols_).adj(); + variRefc_->adj_ += adjAplusc.sum(); + } +}; + +/** + * Return the sum of the specified matrix and specified scalar. + * + * @tparam Mat type of the matrix or expression + * @tparam Scal type of the scalar + * @param m Matrix or expression. + * @param c Scalar. + * @return The matrix plus the scalar. + */ +template * = nullptr, + require_stan_scalar_t* = nullptr, + require_any_var_t::type, Scal>* = nullptr> +auto add(const Mat& A, const Scal c) { + using Ta = value_type_t; + constexpr int R = Mat::RowsAtCompileTime; + constexpr int C = Mat::ColsAtCompileTime; + check_not_nan("add matrix-scalar A + c", "A", A); + check_not_nan("add matrix-scalar A + c", "c", c); + + add_mat_scal_vari* baseVari + = new add_mat_scal_vari(A, c); + Eigen::Matrix Aplusc_v(A.rows(), A.cols()); + Aplusc_v.vi() + = Eigen::Map(&baseVari->variRefAplusc_[0], A.rows(), A.cols()); + + return Aplusc_v; +} + +/** + * Return the sum of the specified scalar and specified matrix. + * + * @tparam Scal type of the scalar + * @tparam Mat type of the matrix or expression + * @param c Scalar. + * @param m Matrix. + * @return The scalar plus the matrix. + */ +template * = nullptr, + require_stan_scalar_t* = nullptr, + require_any_var_t::type, Scal>* = nullptr> +inline auto add(const Scal c, const Mat& A) { + return add(A, c); +} } // namespace math } // namespace stan diff --git a/stan/math/rev/fun/subtract.hpp b/stan/math/rev/fun/subtract.hpp index 837f950d946..98821c9d80a 100644 --- a/stan/math/rev/fun/subtract.hpp +++ b/stan/math/rev/fun/subtract.hpp @@ -25,19 +25,19 @@ namespace math { * @tparam C number of columns (common to A and B) */ - template - class subtract_mat_vari : public vari { +template +class subtract_mat_vari : public vari { public: - int rows_; - int cols_; - int size_; - double* Ad_; - double* Bd_; - vari** variRefA_; - vari** variRefB_; - vari** variRefAminusB_; - - /** + int rows_; + int cols_; + int size_; + double* Ad_; + double* Bd_; + vari** variRefA_; + vari** variRefB_; + vari** variRefAminusB_; + + /** * Constructor for subtract_mat_vari. * * All memeory is allocated in ChainableStack's @@ -53,547 +53,537 @@ namespace math { * @param A matrix * @param B matrix **/ - subtract_mat_vari(const Eigen::Matrix& A, - const Eigen::Matrix& B) - : vari(0.0), - rows_(A.rows()), - cols_(B.cols()), - size_(A.size()), - Ad_(ChainableStack::instance_->memalloc_.alloc_array(size_)), - Bd_(ChainableStack::instance_->memalloc_.alloc_array(size_)), - variRefA_( - ChainableStack::instance_->memalloc_.alloc_array(size_)), - variRefB_( - ChainableStack::instance_->memalloc_.alloc_array(size_)), - variRefAminusB_( - ChainableStack::instance_->memalloc_.alloc_array(size_)) { - using Eigen::Map; - Map(variRefA_, rows_, cols_) = A.vi(); - Map(variRefB_, rows_, cols_) = B.vi(); - Map Ad(Ad_, rows_, cols_); - Map Bd(Bd_, rows_, cols_); - Ad = A.val(); - Bd = B.val(); - Map(variRefAminusB_, rows_, cols_) + subtract_mat_vari(const Eigen::Matrix& A, + const Eigen::Matrix& B) + : vari(0.0), + rows_(A.rows()), + cols_(B.cols()), + size_(A.size()), + Ad_(ChainableStack::instance_->memalloc_.alloc_array(size_)), + Bd_(ChainableStack::instance_->memalloc_.alloc_array(size_)), + variRefA_( + ChainableStack::instance_->memalloc_.alloc_array(size_)), + variRefB_( + ChainableStack::instance_->memalloc_.alloc_array(size_)), + variRefAminusB_( + ChainableStack::instance_->memalloc_.alloc_array(size_)) { + using Eigen::Map; + Map(variRefA_, rows_, cols_) = A.vi(); + Map(variRefB_, rows_, cols_) = B.vi(); + Map Ad(Ad_, rows_, cols_); + Map Bd(Bd_, rows_, cols_); + Ad = A.val(); + Bd = B.val(); + Map(variRefAminusB_, rows_, cols_) = (Ad - Bd).unaryExpr([](double x) { return new vari(x, false); }); + } - } + virtual void chain() { + using Eigen::Map; + matrix_d adjAminusB(rows_, cols_); + adjAminusB = Map(variRefAminusB_, rows_, cols_).adj(); + Map(variRefA_, rows_, cols_).adj() += adjAminusB; + Map(variRefB_, rows_, cols_).adj() -= adjAminusB; + } +}; + +/** + * This is a subclass of the vari class for matrix + * subtraction A - B, where A is an N by M matrix of + * doubles and and B is an N by M matrix of vars. + * + * The class stores the structure of each matrix, + * the double values of A and B, and pointers to + * the varis for A and B if A or B is a var. It + * also instatiates and stores pointers to varis + * for all elements of A - B. + * + * @tparam Tb type of elements in matrix B + * @tparam R number of rows (common to A and B) + * @tparam C number of columns (common to A and B) + */ + +template +class subtract_mat_vari : public vari { + public: + int rows_; + int cols_; + int size_; + double* Ad_; + double* Bd_; + vari** variRefB_; + vari** variRefAminusB_; + + /** + * Constructor for subtract_mat_vari. + * + * All memeory is allocated in ChainableStack's + * stack_alloc arena. + * + * It is critical for the efficiency of this object + * that the constructor create new varis that aren't + * popped onto the var_stack_, but rather are + * popped onto the var_nochain_stack_. This is + * controlled by the second argument to + * vari's constructor. + * + * @param A matrix + * @param B matrix + **/ + subtract_mat_vari(const Eigen::Matrix& A, + const Eigen::Matrix& B) + : vari(0.0), + rows_(A.rows()), + cols_(B.cols()), + size_(A.size()), + Ad_(ChainableStack::instance_->memalloc_.alloc_array(size_)), + Bd_(ChainableStack::instance_->memalloc_.alloc_array(size_)), + variRefB_( + ChainableStack::instance_->memalloc_.alloc_array(size_)), + variRefAminusB_( + ChainableStack::instance_->memalloc_.alloc_array(size_)) { + using Eigen::Map; + Map(variRefB_, rows_, cols_) = B.vi(); + Map Ad(Ad_, rows_, cols_); + Map Bd(Bd_, rows_, cols_); + Ad = A; + Bd = B.val(); + Map(variRefAminusB_, rows_, cols_) + = (Ad - Bd).unaryExpr([](double x) { return new vari(x, false); }); + } + + virtual void chain() { + using Eigen::Map; + matrix_d adjAminusB(rows_, cols_); + adjAminusB = Map(variRefAminusB_, rows_, cols_).adj(); + Map(variRefB_, rows_, cols_).adj() -= adjAminusB; + } +}; + +/** + * This is a subclass of the vari class for matrix + * subtraction A - B, whre A is an N by M matrix of vars + * and B is an N by M matrix of doubles. + * + * The class stores the structure of each matrix, + * the double values of A and B, and pointers to + * the varis for A and B if A or B is a var. It + * also instatiates and stores pointers to varis + * for all elements of A - B. + * + * @tparam Ta type of elements in matrix A + * @tparam R number of rows (common to A and B) + * @tparam C number of columns (common to A and B) + */ + +template +class subtract_mat_vari : public vari { + public: + int rows_; + int cols_; + int size_; + double* Ad_; + double* Bd_; + vari** variRefA_; + vari** variRefAminusB_; + + /** + * Constructor for subtract_mat_vari. + * + * All memeory is allocated in ChainableStack's + * stack_alloc arena. + * + * It is critical for the efficiency of this object + * that the constructor create new varis that aren't + * popped onto the var_stack_, but rather are + * popped onto the var_nochain_stack_. This is + * controlled by the second argument to + * vari's constructor. + * + * @param A matrix + * @param B matrix + **/ + subtract_mat_vari(const Eigen::Matrix& A, + const Eigen::Matrix& B) + : vari(0.0), + rows_(A.rows()), + cols_(B.cols()), + size_(A.size()), + Ad_(ChainableStack::instance_->memalloc_.alloc_array(size_)), + Bd_(ChainableStack::instance_->memalloc_.alloc_array(size_)), + variRefA_( + ChainableStack::instance_->memalloc_.alloc_array(size_)), + variRefAminusB_( + ChainableStack::instance_->memalloc_.alloc_array(size_)) { + using Eigen::Map; + Map(variRefA_, rows_, cols_) = A.vi(); + Map Ad(Ad_, rows_, cols_); + Map Bd(Bd_, rows_, cols_); + Ad = A.val(); + Bd = B; + Map(variRefAminusB_, rows_, cols_) + = (Ad - Bd).unaryExpr([](double x) { return new vari(x, false); }); + } + + virtual void chain() { + using Eigen::Map; + matrix_d adjAminusB(rows_, cols_); + adjAminusB = Map(variRefAminusB_, rows_, cols_).adj(); + Map(variRefA_, rows_, cols_).adj() += adjAminusB; + } +}; + +/** + * Return the sum of two matrices + * @tparam Mat1 type of first matrix + * @tparam Mat2 type of second matrix + * + * @param[in] A matrix + * @param[in] B matrix + * @return Sum of two matrices + */ +template * = nullptr, + require_any_eigen_vt* = nullptr> +auto subtract(const Mat1& A, const Mat2& B) { + using Ta = value_type_t; + using Tb = value_type_t; + constexpr int R = Mat1::RowsAtCompileTime; + constexpr int C = Mat1::ColsAtCompileTime; + check_matching_dims("subtract", "A", A, "B", B); + check_not_nan("subtract", "A", A); + check_not_nan("subtract", "B", B); + + // Manage memory in arena allocator + subtract_mat_vari* baseVari + = new subtract_mat_vari(A, B); + Eigen::Matrix AminusB_v(A.rows(), A.cols()); + AminusB_v.vi() = Eigen::Map(&baseVari->variRefAminusB_[0], + A.rows(), A.cols()); + + return AminusB_v; +} + +/** + * This is a subclass of the vari class for matrix + * subtraction A + c, where A is a N by M matrix and + * c is a scalar. + * + * The class stores the structure of each matrix, + * the double values of A and c, and pointers to + * the varis for A and c if A or c is a var. It + * also instatiates and stores pointers to varis + * for all elements of A + c. + * + * @tparam Ta type of elements in matrix A + * @tparam Tc type of elements in matrix B + * @tparam R number of rows (common to A and B) + * @tparam C number of columns (common to A and B) + * @tparam MatMinusScal: true if computing A - c, false if + * computing c-A + */ + +template +class subtract_mat_scal_vari : public vari { + public: + int rows_; + int cols_; + int size_; + double* Ad_; + double cd_; + vari** variRefA_; + vari* variRefc_; + vari** variRefAminusc_; - virtual void chain() { - using Eigen::Map; - matrix_d adjAminusB(rows_, cols_); - adjAminusB = Map(variRefAminusB_, rows_, cols_).adj(); - Map(variRefA_, rows_, cols_).adj() += adjAminusB; - Map(variRefB_, rows_, cols_).adj() -= adjAminusB; + /** + * Constructor for add_mat_vari. + * + * All memory is allocated in ChainableStack's + * stack_alloc arena. + * + * It is critical for the efficiency of this object + * that the constructor create new varis that aren't + * popped onto the var_stack_, but rather are + * popped onto the var_nochain_stack_. This is + * controlled by the second argument to + * vari's constructor. + * + * @param A matrix + * @param c scalar + **/ + subtract_mat_scal_vari(const Eigen::Matrix& A, const Tc c) + : vari(0.0), + rows_(A.rows()), + cols_(A.cols()), + size_(A.size()), + Ad_(ChainableStack::instance_->memalloc_.alloc_array(size_)), + cd_(c.val()), + variRefA_( + ChainableStack::instance_->memalloc_.alloc_array(size_)), + variRefc_(c.vi_), + variRefAminusc_( + ChainableStack::instance_->memalloc_.alloc_array(size_)) { + using Eigen::Map; + Map(variRefA_, rows_, cols_) = A.vi(); + Map Ad(Ad_, rows_, cols_); + double cd = c.val(); + Ad = A.val(); + if constexpr (MatMinusScal) { + Map(variRefAminusc_, rows_, cols_) + = (Ad.array() - cd).matrix().eval().unaryExpr([](double x) { + return new vari(x, false); + }); + } else { + Map(variRefAminusc_, rows_, cols_) + = (cd - Ad.array()).matrix().eval().unaryExpr([](double x) { + return new vari(x, false); + }); + } + } + + virtual void chain() { + using Eigen::Map; + matrix_d adjAminusc(rows_, cols_); + adjAminusc = Map(variRefAminusc_, rows_, cols_).adj(); + if constexpr (MatMinusScal) { + Map(variRefA_, rows_, cols_).adj() += adjAminusc; + variRefc_->adj_ -= adjAminusc.sum(); + } else { + Map(variRefA_, rows_, cols_).adj() -= adjAminusc; + variRefc_->adj_ += adjAminusc.sum(); } - }; - - /** - * This is a subclass of the vari class for matrix - * subtraction A - B, where A is an N by M matrix of - * doubles and and B is an N by M matrix of vars. - * - * The class stores the structure of each matrix, - * the double values of A and B, and pointers to - * the varis for A and B if A or B is a var. It - * also instatiates and stores pointers to varis - * for all elements of A - B. - * - * @tparam Tb type of elements in matrix B - * @tparam R number of rows (common to A and B) - * @tparam C number of columns (common to A and B) - */ - - template - class subtract_mat_vari : public vari { - public: - int rows_; - int cols_; - int size_; - double* Ad_; - double* Bd_; - vari** variRefB_; - vari** variRefAminusB_; - - /** - * Constructor for subtract_mat_vari. - * - * All memeory is allocated in ChainableStack's - * stack_alloc arena. - * - * It is critical for the efficiency of this object - * that the constructor create new varis that aren't - * popped onto the var_stack_, but rather are - * popped onto the var_nochain_stack_. This is - * controlled by the second argument to - * vari's constructor. - * - * @param A matrix - * @param B matrix - **/ - subtract_mat_vari(const Eigen::Matrix& A, - const Eigen::Matrix& B) - : vari(0.0), - rows_(A.rows()), - cols_(B.cols()), - size_(A.size()), - Ad_(ChainableStack::instance_->memalloc_.alloc_array(size_)), - Bd_(ChainableStack::instance_->memalloc_.alloc_array(size_)), - variRefB_( + } +}; + +/** + * This is a subclass of the vari class for matrix + * subtraction A + c, where A is a N by M matrix and + * c is a double. + * + * The class stores the structure of each matrix, + * the double values of A and c, and pointers to + * the varis for A and c if A or c is a var. It + * also instatiates and stores pointers to varis + * for all elements of A + c. + * + * @tparam Ta type of elements in matrix A + * @tparam Tc type of elements in matrix B + * @tparam R number of rows (common to A and B) + * @tparam C number of columns (common to A and B) + * @tparam MatMinusScal: true if computing A - c, false if + * computing c-A + */ + +template +class subtract_mat_scal_vari : public vari { + public: + int rows_; + int cols_; + int size_; + double* Ad_; + double cd_; + vari** variRefA_; + vari** variRefAminusc_; + + /** + * Constructor for add_mat_vari. + * + * All memory is allocated in ChainableStack's + * stack_alloc arena. + * + * It is critical for the efficiency of this object + * that the constructor create new varis that aren't + * popped onto the var_stack_, but rather are + * popped onto the var_nochain_stack_. This is + * controlled by the second argument to + * vari's constructor. + * + * @param A matrix + * @param c scalar + **/ + subtract_mat_scal_vari(const Eigen::Matrix& A, const double c) + : vari(0.0), + rows_(A.rows()), + cols_(A.cols()), + size_(A.size()), + Ad_(ChainableStack::instance_->memalloc_.alloc_array(size_)), + cd_(c), + variRefA_( ChainableStack::instance_->memalloc_.alloc_array(size_)), - variRefAminusB_( + variRefAminusc_( ChainableStack::instance_->memalloc_.alloc_array(size_)) { - using Eigen::Map; - Map(variRefB_, rows_, cols_) = B.vi(); - Map Ad(Ad_, rows_, cols_); - Map Bd(Bd_, rows_, cols_); - Ad = A; - Bd = B.val(); - Map(variRefAminusB_, rows_, cols_) - = (Ad - Bd).unaryExpr([](double x) { return new vari(x, false); }); - - } - - virtual void chain() { - using Eigen::Map; - matrix_d adjAminusB(rows_, cols_); - adjAminusB = Map(variRefAminusB_, rows_, cols_).adj(); - Map(variRefB_, rows_, cols_).adj() -= adjAminusB; - } - }; + using Eigen::Map; + Map(variRefA_, rows_, cols_) = A.vi(); + Map Ad(Ad_, rows_, cols_); + double cd = c; + Ad = A.val(); + if constexpr (MatMinusScal) { + Map(variRefAminusc_, rows_, cols_) + = (Ad.array() - cd).matrix().eval().unaryExpr([](double x) { + return new vari(x, false); + }); + } else { + Map(variRefAminusc_, rows_, cols_) + = (cd - Ad.array()).matrix().eval().unaryExpr([](double x) { + return new vari(x, false); + }); + } + } + + virtual void chain() { + using Eigen::Map; + matrix_d adjAminusc(rows_, cols_); + adjAminusc = Map(variRefAminusc_, rows_, cols_).adj(); + + if constexpr (MatMinusScal) { + Map(variRefA_, rows_, cols_).adj() += adjAminusc; + } else { + Map(variRefA_, rows_, cols_).adj() -= adjAminusc; + } + } +}; + +/** + * This is a subclass of the vari class for matrix + * subtraction A + c, where A is a N by M matrix of + * doubles c is a scalar var. + * + * The class stores the structure of each matrix, + * the double values of A and c, and pointers to + * the varis for A and c if A or c is a var. It + * also instatiates and stores pointers to varis + * for all elements of A + c. + * + * @tparam Ta type of elements in matrix A + * @tparam Tc type of elements in matrix B + * @tparam R number of rows (common to A and B) + * @tparam C number of columns (common to A and B) + * @tparam MatMinusScal: true if computing A - c, false if + * computing c-A + */ + +template +class subtract_mat_scal_vari : public vari { + public: + int rows_; + int cols_; + int size_; + double* Ad_; + double cd_; + vari* variRefc_; + vari** variRefAminusc_; /** - * This is a subclass of the vari class for matrix - * subtraction A - B, whre A is an N by M matrix of vars - * and B is an N by M matrix of doubles. + * Constructor for add_mat_vari. * - * The class stores the structure of each matrix, - * the double values of A and B, and pointers to - * the varis for A and B if A or B is a var. It - * also instatiates and stores pointers to varis - * for all elements of A - B. + * All memory is allocated in ChainableStack's + * stack_alloc arena. + * + * It is critical for the efficiency of this object + * that the constructor create new varis that aren't + * popped onto the var_stack_, but rather are + * popped onto the var_nochain_stack_. This is + * controlled by the second argument to + * vari's constructor. * - * @tparam Ta type of elements in matrix A - * @tparam R number of rows (common to A and B) - * @tparam C number of columns (common to A and B) - */ - - template - class subtract_mat_vari : public vari { - public: - int rows_; - int cols_; - int size_; - double* Ad_; - double* Bd_; - vari** variRefA_; - vari** variRefAminusB_; - - /** - * Constructor for subtract_mat_vari. - * - * All memeory is allocated in ChainableStack's - * stack_alloc arena. - * - * It is critical for the efficiency of this object - * that the constructor create new varis that aren't - * popped onto the var_stack_, but rather are - * popped onto the var_nochain_stack_. This is - * controlled by the second argument to - * vari's constructor. - * - * @param A matrix - * @param B matrix - **/ - subtract_mat_vari(const Eigen::Matrix& A, - const Eigen::Matrix& B) - : vari(0.0), - rows_(A.rows()), - cols_(B.cols()), - size_(A.size()), - Ad_(ChainableStack::instance_->memalloc_.alloc_array(size_)), - Bd_(ChainableStack::instance_->memalloc_.alloc_array(size_)), - variRefA_( - ChainableStack::instance_->memalloc_.alloc_array(size_)), - variRefAminusB_( - ChainableStack::instance_->memalloc_.alloc_array(size_)) { - using Eigen::Map; - Map(variRefA_, rows_, cols_) = A.vi(); - Map Ad(Ad_, rows_, cols_); - Map Bd(Bd_, rows_, cols_); - Ad = A.val(); - Bd = B; - Map(variRefAminusB_, rows_, cols_) - = (Ad - Bd).unaryExpr([](double x) { return new vari(x, false); }); - - } - - virtual void chain() { - using Eigen::Map; - matrix_d adjAminusB(rows_, cols_); - adjAminusB = Map(variRefAminusB_, rows_, cols_).adj(); - Map(variRefA_, rows_, cols_).adj() += adjAminusB; - } - }; - - /** - * Return the sum of two matrices - * @tparam Mat1 type of first matrix - * @tparam Mat2 type of second matrix - * - * @param[in] A matrix - * @param[in] B matrix - * @return Sum of two matrices - */ - template * = nullptr, - require_any_eigen_vt* = nullptr> - auto subtract(const Mat1& A, const Mat2& B) { - using Ta = value_type_t; - using Tb = value_type_t; - constexpr int R = Mat1::RowsAtCompileTime; - constexpr int C = Mat1::ColsAtCompileTime; - check_matching_dims("subtract", "A", A, "B", B); - check_not_nan("subtract", "A", A); - check_not_nan("subtract", "B", B); - - // Manage memory in arena allocator - subtract_mat_vari* baseVari - = new subtract_mat_vari(A, B); - Eigen::Matrix AminusB_v(A.rows(), A.cols()); - AminusB_v.vi() = Eigen::Map(&baseVari->variRefAminusB_[0], - A.rows(), A.cols()); - - return AminusB_v; + * @param A matrix + * @param c scalar + **/ + subtract_mat_scal_vari(const Eigen::Matrix& A, const Tc c) + : vari(0.0), + rows_(A.rows()), + cols_(A.cols()), + size_(A.size()), + Ad_(ChainableStack::instance_->memalloc_.alloc_array(size_)), + cd_(c.val()), + variRefc_(c.vi_), + variRefAminusc_( + ChainableStack::instance_->memalloc_.alloc_array(size_)) { + using Eigen::Map; + + double cd = c.val(); + Map Ad(Ad_, rows_, cols_); + Ad = A.val(); + if constexpr (MatMinusScal) { + Map(variRefAminusc_, rows_, cols_) + = (Ad.array() - cd).matrix().eval().unaryExpr([](double x) { + return new vari(x, false); + }); + } else { + Map(variRefAminusc_, rows_, cols_) + = (cd - Ad.array()).matrix().eval().unaryExpr([](double x) { + return new vari(x, false); + }); } + } + + virtual void chain() { + using Eigen::Map; + matrix_d adjAminusc(rows_, cols_); + adjAminusc = Map(variRefAminusc_, rows_, cols_).adj(); + if constexpr (MatMinusScal) { + variRefc_->adj_ -= adjAminusc.sum(); + } else { + variRefc_->adj_ += adjAminusc.sum(); + } + } +}; - /** - * This is a subclass of the vari class for matrix - * subtraction A + c, where A is a N by M matrix and - * c is a scalar. - * - * The class stores the structure of each matrix, - * the double values of A and c, and pointers to - * the varis for A and c if A or c is a var. It - * also instatiates and stores pointers to varis - * for all elements of A + c. - * - * @tparam Ta type of elements in matrix A - * @tparam Tc type of elements in matrix B - * @tparam R number of rows (common to A and B) - * @tparam C number of columns (common to A and B) - * @tparam MatMinusScal: true if computing A - c, false if - * computing c-A - */ - - template - class subtract_mat_scal_vari : public vari { - public: - int rows_; - int cols_; - int size_; - double* Ad_; - double cd_; - vari** variRefA_; - vari* variRefc_; - vari** variRefAminusc_; - - /** - * Constructor for add_mat_vari. - * - * All memory is allocated in ChainableStack's - * stack_alloc arena. - * - * It is critical for the efficiency of this object - * that the constructor create new varis that aren't - * popped onto the var_stack_, but rather are - * popped onto the var_nochain_stack_. This is - * controlled by the second argument to - * vari's constructor. - * - * @param A matrix - * @param c scalar - **/ - subtract_mat_scal_vari(const Eigen::Matrix& A, - const Tc c) - : vari(0.0), - rows_(A.rows()), - cols_(A.cols()), - size_(A.size()), - Ad_(ChainableStack::instance_->memalloc_.alloc_array(size_)), - cd_(c.val()), - variRefA_( - ChainableStack::instance_->memalloc_.alloc_array(size_)), - variRefc_(c.vi_), - variRefAminusc_( - ChainableStack::instance_->memalloc_.alloc_array(size_)) { - using Eigen::Map; - Map(variRefA_, rows_, cols_) = A.vi(); - Map Ad(Ad_, rows_, cols_); - double cd = c.val(); - Ad = A.val(); - if constexpr (MatMinusScal){ - Map(variRefAminusc_, rows_, cols_) - = (Ad.array() - cd).matrix().eval().unaryExpr([](double x) - { return new vari(x, false); }); - } else { - Map(variRefAminusc_, rows_, cols_) - = (cd - Ad.array()).matrix().eval().unaryExpr([](double x) - { return new vari(x, false); }); - } - - } - - virtual void chain() { - using Eigen::Map; - matrix_d adjAminusc(rows_, cols_); - adjAminusc = Map(variRefAminusc_, rows_, cols_).adj(); - if constexpr (MatMinusScal) { - Map(variRefA_, rows_, cols_).adj() += adjAminusc; - variRefc_->adj_ -= adjAminusc.sum(); - } else { - Map(variRefA_, rows_, cols_).adj() -= adjAminusc; - variRefc_->adj_ += adjAminusc.sum(); - } - - } - }; - - /** - * This is a subclass of the vari class for matrix - * subtraction A + c, where A is a N by M matrix and - * c is a double. - * - * The class stores the structure of each matrix, - * the double values of A and c, and pointers to - * the varis for A and c if A or c is a var. It - * also instatiates and stores pointers to varis - * for all elements of A + c. - * - * @tparam Ta type of elements in matrix A - * @tparam Tc type of elements in matrix B - * @tparam R number of rows (common to A and B) - * @tparam C number of columns (common to A and B) - * @tparam MatMinusScal: true if computing A - c, false if - * computing c-A - */ - - template - class subtract_mat_scal_vari - : public vari { - public: - int rows_; - int cols_; - int size_; - double* Ad_; - double cd_; - vari** variRefA_; - vari** variRefAminusc_; - - /** - * Constructor for add_mat_vari. - * - * All memory is allocated in ChainableStack's - * stack_alloc arena. - * - * It is critical for the efficiency of this object - * that the constructor create new varis that aren't - * popped onto the var_stack_, but rather are - * popped onto the var_nochain_stack_. This is - * controlled by the second argument to - * vari's constructor. - * - * @param A matrix - * @param c scalar - **/ - subtract_mat_scal_vari(const Eigen::Matrix& A, - const double c) - : vari(0.0), - rows_(A.rows()), - cols_(A.cols()), - size_(A.size()), - Ad_(ChainableStack::instance_->memalloc_.alloc_array(size_)), - cd_(c), - variRefA_( - ChainableStack::instance_->memalloc_.alloc_array(size_)), - variRefAminusc_( - ChainableStack::instance_->memalloc_.alloc_array(size_)) { - using Eigen::Map; - Map(variRefA_, rows_, cols_) = A.vi(); - Map Ad(Ad_, rows_, cols_); - double cd = c; - Ad = A.val(); - if constexpr (MatMinusScal){ - Map(variRefAminusc_, rows_, cols_) - = (Ad.array() - cd).matrix().eval().unaryExpr([](double x) - { return new vari(x, false); }); - } else { - Map(variRefAminusc_, rows_, cols_) - = (cd - Ad.array()).matrix().eval().unaryExpr([](double x) - { return new vari(x, false); }); - } - - } - - virtual void chain() { - using Eigen::Map; - matrix_d adjAminusc(rows_, cols_); - adjAminusc = Map(variRefAminusc_, rows_, cols_).adj(); - - if constexpr (MatMinusScal) { - Map(variRefA_, rows_, cols_).adj() += adjAminusc; - } else { - Map(variRefA_, rows_, cols_).adj() -= adjAminusc; - } - } - }; - - - /** - * This is a subclass of the vari class for matrix - * subtraction A + c, where A is a N by M matrix of - * doubles c is a scalar var. - * - * The class stores the structure of each matrix, - * the double values of A and c, and pointers to - * the varis for A and c if A or c is a var. It - * also instatiates and stores pointers to varis - * for all elements of A + c. - * - * @tparam Ta type of elements in matrix A - * @tparam Tc type of elements in matrix B - * @tparam R number of rows (common to A and B) - * @tparam C number of columns (common to A and B) - * @tparam MatMinusScal: true if computing A - c, false if - * computing c-A - */ - - template - class subtract_mat_scal_vari - : public vari { - public: - int rows_; - int cols_; - int size_; - double* Ad_; - double cd_; - vari* variRefc_; - vari** variRefAminusc_; - - /** - * Constructor for add_mat_vari. - * - * All memory is allocated in ChainableStack's - * stack_alloc arena. - * - * It is critical for the efficiency of this object - * that the constructor create new varis that aren't - * popped onto the var_stack_, but rather are - * popped onto the var_nochain_stack_. This is - * controlled by the second argument to - * vari's constructor. - * - * @param A matrix - * @param c scalar - **/ - subtract_mat_scal_vari(const Eigen::Matrix& A, - const Tc c) - : vari(0.0), - rows_(A.rows()), - cols_(A.cols()), - size_(A.size()), - Ad_(ChainableStack::instance_->memalloc_.alloc_array(size_)), - cd_(c.val()), - variRefc_(c.vi_), - variRefAminusc_( - ChainableStack::instance_->memalloc_.alloc_array(size_)) { - using Eigen::Map; - - double cd = c.val(); - Map Ad(Ad_, rows_, cols_); - Ad = A.val(); - if constexpr (MatMinusScal){ - Map(variRefAminusc_, rows_, cols_) - = (Ad.array() - cd).matrix().eval().unaryExpr([](double x) - { return new vari(x, false); }); - } else { - Map(variRefAminusc_, rows_, cols_) - = (cd - Ad.array()).matrix().eval().unaryExpr([](double x) - { return new vari(x, false); }); - } - - } - - virtual void chain() { - using Eigen::Map; - matrix_d adjAminusc(rows_, cols_); - adjAminusc = Map(variRefAminusc_, rows_, cols_).adj(); - if constexpr (MatMinusScal) { - variRefc_->adj_ -= adjAminusc.sum(); - } else { - variRefc_->adj_ += adjAminusc.sum(); - } - } - }; - - - /** - * Return the difference A - c of a matrix A and a scalar c. - * - * @tparam Mat type of the matrix or expression - * @tparam Scal type of the scalar - * @param m Matrix or expression. - * @param c Scalar. - * @return The matrix minus the scalar. - */ - template * = nullptr, - require_stan_scalar_t* = nullptr, - require_any_var_t::type, Scal>* = nullptr> - auto subtract(const Mat& A, const Scal c) { - using Ta = value_type_t; - constexpr int R = Mat::RowsAtCompileTime; - constexpr int C = Mat::ColsAtCompileTime; - check_not_nan("subtract matrix-scalar A - c", "A", A); - check_not_nan("subtract matrix-scalar A - c", "c", c); - - subtract_mat_scal_vari* baseVari - = new subtract_mat_scal_vari(A, c); - Eigen::Matrix Aminusc_v(A.rows(), A.cols()); - Aminusc_v.vi() = Eigen::Map( - &baseVari->variRefAminusc_[0], A.rows(), A.cols()); - - return Aminusc_v; - } - - /** - * Return the difference c - A of a scalar c and a matrix A. - * - * @tparam Scal type of the scalar - * @tparam Mat type of the matrix or expression - * @param c Scalar. - * @param m Matrix. - * @return The scalar minus the matrix. - */ - template * = nullptr, - require_stan_scalar_t* = nullptr, - require_any_var_t::type, Scal>* = nullptr> - auto subtract(const Scal c, const Mat& A) { - using Ta = value_type_t; - constexpr int R = Mat::RowsAtCompileTime; - constexpr int C = Mat::ColsAtCompileTime; - check_not_nan("subtract matrix-scalar A - c", "A", A); - check_not_nan("subtract matrix-scalar A - c", "c", c); - - subtract_mat_scal_vari* baseVari - = new subtract_mat_scal_vari(A, c); - Eigen::Matrix Aminusc_v(A.rows(), A.cols()); - Aminusc_v.vi() = Eigen::Map( - &baseVari->variRefAminusc_[0], A.rows(), A.cols()); - - return Aminusc_v; - } +/** + * Return the difference A - c of a matrix A and a scalar c. + * + * @tparam Mat type of the matrix or expression + * @tparam Scal type of the scalar + * @param m Matrix or expression. + * @param c Scalar. + * @return The matrix minus the scalar. + */ +template * = nullptr, + require_stan_scalar_t* = nullptr, + require_any_var_t::type, Scal>* = nullptr> +auto subtract(const Mat& A, const Scal c) { + using Ta = value_type_t; + constexpr int R = Mat::RowsAtCompileTime; + constexpr int C = Mat::ColsAtCompileTime; + check_not_nan("subtract matrix-scalar A - c", "A", A); + check_not_nan("subtract matrix-scalar A - c", "c", c); + + subtract_mat_scal_vari* baseVari + = new subtract_mat_scal_vari(A, c); + Eigen::Matrix Aminusc_v(A.rows(), A.cols()); + Aminusc_v.vi() = Eigen::Map(&baseVari->variRefAminusc_[0], + A.rows(), A.cols()); + + return Aminusc_v; +} + +/** + * Return the difference c - A of a scalar c and a matrix A. + * + * @tparam Scal type of the scalar + * @tparam Mat type of the matrix or expression + * @param c Scalar. + * @param m Matrix. + * @return The scalar minus the matrix. + */ +template * = nullptr, + require_stan_scalar_t* = nullptr, + require_any_var_t::type, Scal>* = nullptr> +auto subtract(const Scal c, const Mat& A) { + using Ta = value_type_t; + constexpr int R = Mat::RowsAtCompileTime; + constexpr int C = Mat::ColsAtCompileTime; + check_not_nan("subtract matrix-scalar A - c", "A", A); + check_not_nan("subtract matrix-scalar A - c", "c", c); + + subtract_mat_scal_vari* baseVari + = new subtract_mat_scal_vari(A, c); + Eigen::Matrix Aminusc_v(A.rows(), A.cols()); + Aminusc_v.vi() = Eigen::Map(&baseVari->variRefAminusc_[0], + A.rows(), A.cols()); + + return Aminusc_v; +} } // namespace math } // namespace stan diff --git a/test/unit/math/rev/fun/add_test.cpp b/test/unit/math/rev/fun/add_test.cpp index 1868881f48f..a27ed8143f5 100644 --- a/test/unit/math/rev/fun/add_test.cpp +++ b/test/unit/math/rev/fun/add_test.cpp @@ -1,14 +1,13 @@ #include #include - TEST(MathRevAdd, add_v_d) { + using stan::math::add; using stan::math::matrix_d; using stan::math::matrix_v; using stan::math::var; - using stan::math::add; - matrix_d Ad = Eigen::MatrixXd::Random(4,6); - matrix_d Bd = Eigen::MatrixXd::Random(4,6); + matrix_d Ad = Eigen::MatrixXd::Random(4, 6); + matrix_d Bd = Eigen::MatrixXd::Random(4, 6); var c = 3.2; var d = -2.73; matrix_v Av = c * Ad; @@ -16,65 +15,64 @@ TEST(MathRevAdd, add_v_d) { matrix_v Cvd = add(Av, Bd); matrix_v Cdv = add(Ad, Bv); matrix_v Cvv = add(Av, Bv); - matrix_v Cdvd = add(Ad, - add(Bv, Eigen::MatrixXd::Random(4,6))); + matrix_v Cdvd = add(Ad, add(Bv, Eigen::MatrixXd::Random(4, 6))); matrix_v Adcv = add(c, Ad); matrix_v Bvcv = add(c, Bv); matrix_v Bvscal = add(0.7, Bv); - Cvd(0,1).grad(); - EXPECT_FLOAT_EQ(c.adj(), Ad(0,1)); + Cvd(0, 1).grad(); + EXPECT_FLOAT_EQ(c.adj(), Ad(0, 1)); EXPECT_FLOAT_EQ(d.adj(), 0.0); - stan::math::set_zero_all_adjoints(); - Cdv(2,3).grad(); - EXPECT_FLOAT_EQ(d.adj(), Bd(2,3)); + stan::math::set_zero_all_adjoints(); + Cdv(2, 3).grad(); + EXPECT_FLOAT_EQ(d.adj(), Bd(2, 3)); EXPECT_FLOAT_EQ(c.adj(), 0.0); stan::math::set_zero_all_adjoints(); - Cvv(3,2).grad(); - EXPECT_FLOAT_EQ(d.adj(), Bd(3,2)); - EXPECT_FLOAT_EQ(c.adj(), Ad(3,2)); + Cvv(3, 2).grad(); + EXPECT_FLOAT_EQ(d.adj(), Bd(3, 2)); + EXPECT_FLOAT_EQ(c.adj(), Ad(3, 2)); stan::math::set_zero_all_adjoints(); - Cdvd(3,4).grad(); - EXPECT_FLOAT_EQ(d.adj(), Bd(3,4)); + Cdvd(3, 4).grad(); + EXPECT_FLOAT_EQ(d.adj(), Bd(3, 4)); EXPECT_FLOAT_EQ(c.adj(), 0.0); stan::math::set_zero_all_adjoints(); - Adcv(1,1).grad(); + Adcv(1, 1).grad(); EXPECT_FLOAT_EQ(c.adj(), 1.0); EXPECT_FLOAT_EQ(d.adj(), 0.0); stan::math::set_zero_all_adjoints(); - Bvcv(2,1).grad(); + Bvcv(2, 1).grad(); EXPECT_FLOAT_EQ(c.adj(), 1.0); - EXPECT_FLOAT_EQ(d.adj(), Bd(2,1)); + EXPECT_FLOAT_EQ(d.adj(), Bd(2, 1)); stan::math::set_zero_all_adjoints(); - Bvscal(2,0).grad(); - EXPECT_FLOAT_EQ(d.adj(), Bd(2,0)); + Bvscal(2, 0).grad(); + EXPECT_FLOAT_EQ(d.adj(), Bd(2, 0)); EXPECT_FLOAT_EQ(c.adj(), 0.0); } TEST(MathRevAdd, add_check_dim) { + using stan::math::add; using stan::math::matrix_d; using stan::math::matrix_v; - using stan::math::add; matrix_v A, B; A.resize(3, 2); B.resize(3, 2); - EXPECT_NO_THROW(add(A,B)); + EXPECT_NO_THROW(add(A, B)); A.resize(1, 0); B.resize(1, 0); - EXPECT_NO_THROW(add(A,B)); + EXPECT_NO_THROW(add(A, B)); A.resize(0, 0); B.resize(0, 0); - EXPECT_NO_THROW(add(A,B)); + EXPECT_NO_THROW(add(A, B)); A.resize(2, 3); B.resize(3, 2); - EXPECT_THROW(add(A,B), std::invalid_argument); + EXPECT_THROW(add(A, B), std::invalid_argument); }; diff --git a/test/unit/math/rev/fun/subtract_test.cpp b/test/unit/math/rev/fun/subtract_test.cpp index 21ea614459f..9e73d317381 100644 --- a/test/unit/math/rev/fun/subtract_test.cpp +++ b/test/unit/math/rev/fun/subtract_test.cpp @@ -1,14 +1,13 @@ #include #include - TEST(MathRevSubtract, subtract_v_d) { using stan::math::matrix_d; using stan::math::matrix_v; - using stan::math::var; using stan::math::subtract; - matrix_d Ad = Eigen::MatrixXd::Random(4,6); - matrix_d Bd = Eigen::MatrixXd::Random(4,6); + using stan::math::var; + matrix_d Ad = Eigen::MatrixXd::Random(4, 6); + matrix_d Bd = Eigen::MatrixXd::Random(4, 6); var c(3.2); var d(-2.73); matrix_v Av = c * Ad; @@ -16,8 +15,8 @@ TEST(MathRevSubtract, subtract_v_d) { matrix_v Cvd = subtract(Av, Bd); matrix_v Cdv = subtract(Ad, Bv); matrix_v Cvv = subtract(Av, Bv); - matrix_v Cdvd = subtract(Ad, subtract(Bv, Eigen::MatrixXd::Random(4,6))); - matrix_v Cddv = subtract(Ad, subtract(Eigen::MatrixXd::Random(4,6), Bv)); + matrix_v Cdvd = subtract(Ad, subtract(Bv, Eigen::MatrixXd::Random(4, 6))); + matrix_v Cddv = subtract(Ad, subtract(Eigen::MatrixXd::Random(4, 6), Bv)); matrix_v Adcv = subtract(Ad, c); matrix_v cvAd = subtract(c, Ad); matrix_v Bvcv = subtract(Bv, c); @@ -25,65 +24,57 @@ TEST(MathRevSubtract, subtract_v_d) { matrix_v Bvscal = subtract(Bv, 0.7); matrix_v scalBv = subtract(0.7, Bv); - Cvd(0,1).grad(); - EXPECT_FLOAT_EQ(c.adj(), Ad(0,1)); + Cvd(0, 1).grad(); + EXPECT_FLOAT_EQ(c.adj(), Ad(0, 1)); stan::math::set_zero_all_adjoints(); - Cdv(2,3).grad(); - EXPECT_FLOAT_EQ(d.adj(), -Bd(2,3)); + Cdv(2, 3).grad(); + EXPECT_FLOAT_EQ(d.adj(), -Bd(2, 3)); stan::math::set_zero_all_adjoints(); - Cvv(3,2).grad(); - EXPECT_FLOAT_EQ(d.adj(), -Bd(3,2)); - EXPECT_FLOAT_EQ(c.adj(), Ad(3,2)); + Cvv(3, 2).grad(); + EXPECT_FLOAT_EQ(d.adj(), -Bd(3, 2)); + EXPECT_FLOAT_EQ(c.adj(), Ad(3, 2)); stan::math::set_zero_all_adjoints(); - Cdvd(3,4).grad(); - EXPECT_FLOAT_EQ(d.adj(), -Bd(3,4)); + Cdvd(3, 4).grad(); + EXPECT_FLOAT_EQ(d.adj(), -Bd(3, 4)); stan::math::set_zero_all_adjoints(); - Cddv(2,5).grad(); - EXPECT_FLOAT_EQ(d.adj(), Bd(2,5)); - - - - - - + Cddv(2, 5).grad(); + EXPECT_FLOAT_EQ(d.adj(), Bd(2, 5)); stan::math::set_zero_all_adjoints(); - Adcv(1,1).grad(); + Adcv(1, 1).grad(); EXPECT_FLOAT_EQ(c.adj(), -1.0); EXPECT_FLOAT_EQ(d.adj(), 0.0); stan::math::set_zero_all_adjoints(); - cvAd(1,1).grad(); + cvAd(1, 1).grad(); EXPECT_FLOAT_EQ(c.adj(), 1.0); EXPECT_FLOAT_EQ(d.adj(), 0.0); stan::math::set_zero_all_adjoints(); - Bvcv(2,1).grad(); + Bvcv(2, 1).grad(); EXPECT_FLOAT_EQ(c.adj(), -1.0); - EXPECT_FLOAT_EQ(d.adj(), Bd(2,1)); + EXPECT_FLOAT_EQ(d.adj(), Bd(2, 1)); stan::math::set_zero_all_adjoints(); - cvBv(2,1).grad(); + cvBv(2, 1).grad(); EXPECT_FLOAT_EQ(c.adj(), 1.0); - EXPECT_FLOAT_EQ(d.adj(), -Bd(2,1)); + EXPECT_FLOAT_EQ(d.adj(), -Bd(2, 1)); stan::math::set_zero_all_adjoints(); - Bvscal(2,0).grad(); - EXPECT_FLOAT_EQ(d.adj(), Bd(2,0)); + Bvscal(2, 0).grad(); + EXPECT_FLOAT_EQ(d.adj(), Bd(2, 0)); EXPECT_FLOAT_EQ(c.adj(), 0.0); stan::math::set_zero_all_adjoints(); - scalBv(2,0).grad(); - EXPECT_FLOAT_EQ(d.adj(), -Bd(2,0)); + scalBv(2, 0).grad(); + EXPECT_FLOAT_EQ(d.adj(), -Bd(2, 0)); EXPECT_FLOAT_EQ(c.adj(), 0.0); } - - TEST(MathRevSubtract, subtract_check_dim) { using stan::math::matrix_d; using stan::math::matrix_v; @@ -92,17 +83,17 @@ TEST(MathRevSubtract, subtract_check_dim) { A.resize(3, 2); B.resize(3, 2); - EXPECT_NO_THROW(subtract(A,B)); + EXPECT_NO_THROW(subtract(A, B)); A.resize(1, 0); B.resize(1, 0); - EXPECT_NO_THROW(subtract(A,B)); + EXPECT_NO_THROW(subtract(A, B)); A.resize(0, 0); B.resize(0, 0); - EXPECT_NO_THROW(subtract(A,B)); + EXPECT_NO_THROW(subtract(A, B)); A.resize(2, 3); B.resize(3, 2); - EXPECT_THROW(subtract(A,B), std::invalid_argument); + EXPECT_THROW(subtract(A, B), std::invalid_argument); };