diff --git a/stan/math/opencl/copy.hpp b/stan/math/opencl/copy.hpp index 6a4b8657c77..e795f7c1b04 100644 --- a/stan/math/opencl/copy.hpp +++ b/stan/math/opencl/copy.hpp @@ -11,6 +11,7 @@ #include #include #include +#include #include #include #include @@ -101,6 +102,23 @@ inline Eigen::Matrix from_matrix_cl(const matrix_cl& src) { return dst; } +/** \ingroup opencl + * Copies result of a kernel generator expression to the + * destination Eigen matrix. + * + * @tparam R rows type of the destination + * @tparam C cols type of the destination + * @tparam T type of expression + * @param src source expression + * @return Eigen matrix with a copy of the data in the source matrix + */ +template * = nullptr, + require_not_matrix_cl_t* = nullptr> +inline Eigen::Matrix, R, C> from_matrix_cl(const T& src) { + return from_matrix_cl(src.eval()); +} + /** \ingroup opencl * Packs the flat triangular matrix on the OpenCL device and * copies it to the std::vector. diff --git a/stan/math/opencl/kernel_generator.hpp b/stan/math/opencl/kernel_generator.hpp index 8fe240391c9..eedd39d9b03 100644 --- a/stan/math/opencl/kernel_generator.hpp +++ b/stan/math/opencl/kernel_generator.hpp @@ -114,6 +114,7 @@ #include #include +#include #include #include #include diff --git a/stan/math/opencl/kernel_generator/constant.hpp b/stan/math/opencl/kernel_generator/constant.hpp new file mode 100644 index 00000000000..ff9e48f07e7 --- /dev/null +++ b/stan/math/opencl/kernel_generator/constant.hpp @@ -0,0 +1,132 @@ +#ifndef STAN_MATH_OPENCL_KERNEL_GENERATOR_CONSTANT_HPP +#define STAN_MATH_OPENCL_KERNEL_GENERATOR_CONSTANT_HPP +#ifdef STAN_OPENCL + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +namespace stan { +namespace math { + +/** \addtogroup opencl_kernel_generator + * @{ + */ + +/** + * Represents a matrix of single repeated value in kernel generator expressions. + * @tparam T type of the scalar + */ +template +class constant_ : public operation_cl, T> { + T a_; + int rows_; + int cols_; + + public: + static_assert(std::is_arithmetic::value, + "class scalar_: std::is_arithmetic must be true!"); + using Scalar = T; + using base = operation_cl, T>; + using base::var_name_; + + /** + * Constructor + * @param a scalar value + * @param rows number of rows of the matrix + * @param cols number of columns of the matrix + */ + explicit constant_(const T a, int rows, int cols) + : a_(a), rows_(rows), cols_(cols) {} + + /** + * Creates a deep copy of this expression. + * @return copy of \c *this + */ + inline constant_ deep_copy() const { + return constant_(a_, rows_, cols_); + } + + /** + * Generates kernel code for this expression. + * @param row_index_name row index variable name + * @param col_index_name column index variable name + * @param view_handled whether whether caller already handled matrix view + * @return part of kernel with code for this expression + */ + inline kernel_parts generate(const std::string& row_index_name, + const std::string& col_index_name, + const bool view_handled) const { + kernel_parts res{}; + res.args = type_str() + " " + var_name_ + ", "; + return res; + } + + /** + * Sets kernel arguments for this and nested expressions. + * @param[in,out] generated set of expressions that already set their kernel + * arguments + * @param kernel kernel to set arguments on + * @param[in,out] arg_num consecutive number of the first argument to set. + * This is incremented for each argument set by this function. + */ + inline void set_args(std::set& generated, + cl::Kernel& kernel, int& arg_num) const { + kernel.setArg(arg_num++, a_); + } + + /** + * Number of rows of a matrix that would be the result of evaluating this + * expression. + * @return number of rows + */ + inline int rows() const { return rows_; } + + /** + * Number of columns of a matrix that would be the result of evaluating this + * expression. + * @return number of columns + */ + inline int cols() const { return cols_; } + + /** + * Determine indices of extreme sub- and superdiagonals written. + * @return pair of indices - bottom and top diagonal + */ + inline std::pair extreme_diagonals() const { + return {std::numeric_limits::min(), std::numeric_limits::max()}; + } +}; + +/** + * Matrix of repeated values in kernel generator expressions. + * + * In most cases scalars should be directly used instead of this. This is, + * however, useful for initializing some expression to specific value if that + * expresssion could also be plain `matrix_cl`. + * + * @tparam T type of argument + * @param a input argument + * @param rows number of rows + * @param cols number of columns + * @return Block of given expression + */ +template > +inline auto constant(const T a, int rows, int cols) { + return constant_(a, rows, cols); +} + +/** @}*/ +} // namespace math +} // namespace stan + +#endif +#endif diff --git a/stan/math/opencl/kernel_generator/is_kernel_expression.hpp b/stan/math/opencl/kernel_generator/is_kernel_expression.hpp index 5bb54c818d2..186f3f84ba1 100644 --- a/stan/math/opencl/kernel_generator/is_kernel_expression.hpp +++ b/stan/math/opencl/kernel_generator/is_kernel_expression.hpp @@ -7,16 +7,21 @@ #include namespace stan { -namespace math { /** \addtogroup opencl_kernel_generator * @{ */ + /** - * Non-templated base of \c operation is needed for easy checking if something - * is a subclass of \c operation. + * Non-templated base of `operation_cl` is needed for easy checking if something + * is a subclass of `operation_cl`. */ class operation_cl_base {}; +/** + * Non-templated base of `operation_cl_lhs` is needed for easy checking if + * something is a subclass of `operation_cl_lhs`. + */ +class operation_cl_lhs_base {}; /** * Determines whether a type is non-scalar type that is a valid kernel generator @@ -26,7 +31,6 @@ template struct is_kernel_expression_and_not_scalar : bool_constant>::value> {}; - template struct is_kernel_expression_and_not_scalar> : std::true_type {}; @@ -57,8 +61,21 @@ using require_all_kernel_expressions_and_none_scalar_t template using require_all_kernel_expressions_t = require_all_t...>; + +/** + * Determines whether a type is an assignable kernel generator + * expression. + */ +template +struct is_kernel_expression_lhs + : bool_constant>::value> {}; +template +struct is_kernel_expression_lhs> : std::true_type {}; + /** @}*/ -} // namespace math +STAN_ADD_REQUIRE_UNARY(kernel_expression_lhs, is_kernel_expression_lhs, + opencl_kernel_generator); } // namespace stan #endif diff --git a/stan/math/opencl/kernel_generator/operation_cl.hpp b/stan/math/opencl/kernel_generator/operation_cl.hpp index 846e4403606..a20478a793b 100644 --- a/stan/math/opencl/kernel_generator/operation_cl.hpp +++ b/stan/math/opencl/kernel_generator/operation_cl.hpp @@ -314,6 +314,8 @@ class operation_cl : public operation_cl_base { * @return number of rows */ inline int rows() const { + static_assert( + N > 0, "default rows does not work on expressions with no arguments!"); return index_apply([&](auto... Is) { // assuming all non-dynamic sizes match return std::max({this->get_arg().rows()...}); @@ -326,6 +328,8 @@ class operation_cl : public operation_cl_base { * @return number of columns */ inline int cols() const { + static_assert( + N > 0, "default cols does not work on expressions with no arguments!"); return index_apply([&](auto... Is) { // assuming all non-dynamic sizes match return std::max({this->get_arg().cols()...}); @@ -352,6 +356,9 @@ class operation_cl : public operation_cl_base { * @return pair of indices - bottom and top diagonal */ inline std::pair extreme_diagonals() const { + static_assert(N > 0, + "default extreme_diagonals does not work on expressions with " + "no arguments!"); return index_apply([&](auto... Is) { auto arg_diags = std::make_tuple(this->get_arg().extreme_diagonals()...); diff --git a/stan/math/opencl/kernel_generator/operation_cl_lhs.hpp b/stan/math/opencl/kernel_generator/operation_cl_lhs.hpp index ef2c73b812d..c7b47d8bb98 100644 --- a/stan/math/opencl/kernel_generator/operation_cl_lhs.hpp +++ b/stan/math/opencl/kernel_generator/operation_cl_lhs.hpp @@ -24,7 +24,8 @@ namespace math { * @tparam Args types of arguments to this operation */ template -class operation_cl_lhs : public operation_cl { +class operation_cl_lhs : public operation_cl, + public operation_cl_lhs_base { protected: using base = operation_cl; static constexpr int N = sizeof...(Args); diff --git a/stan/math/opencl/kernel_generator/scalar.hpp b/stan/math/opencl/kernel_generator/scalar.hpp index 2ddd13b5d01..5282b1f3c4f 100644 --- a/stan/math/opencl/kernel_generator/scalar.hpp +++ b/stan/math/opencl/kernel_generator/scalar.hpp @@ -27,7 +27,7 @@ namespace math { */ template class scalar_ : public operation_cl, T> { - private: + protected: T a_; public: diff --git a/stan/math/opencl/rev/copy.hpp b/stan/math/opencl/rev/copy.hpp index 2dda8fe80d5..b343a4a06f9 100644 --- a/stan/math/opencl/rev/copy.hpp +++ b/stan/math/opencl/rev/copy.hpp @@ -16,6 +16,76 @@ namespace stan { namespace math { +namespace internal { +template * = nullptr> +class op_copy_to_cl_vari final + : public vari_value>> { + T_arg_adj arg_adj_; + + public: + template * = nullptr, + require_vt_same* = nullptr> + explicit op_copy_to_cl_vari(const T_arg_val& val, T_arg_adj adj) + : vari_value>>(to_matrix_cl(val)), + arg_adj_(adj) {} + + virtual void chain() { + arg_adj_ += from_matrix_cl::RowsAtCompileTime, + std::decay_t::ColsAtCompileTime>( + this->adj_); + } +}; +} // namespace internal + +/** \ingroup opencl + * Copies the source var containing Eigen matrices to destination var that has + * data stored on the OpenCL device. + * + * @tparam T type of the Eigen matrix + * @param a source Eigen matrix + * @return var with a copy of the data on the OpenCL device + */ +template +inline var_value>> to_matrix_cl( + const var_value& a) { + return new internal::op_copy_to_cl_variadj_)>(a.val(), + a.vi_->adj_); +} + +namespace internal { +template * = nullptr> +class op_copy_from_cl_vari final + : public vari_value, Rows, Cols>> { + vari_value& a_; + + public: + explicit op_copy_from_cl_vari(vari_value& a) + : vari_value, Rows, Cols>>( + from_matrix_cl(a.val_)), + a_(a) {} + + virtual void chain() { a_.adj_ = a_.adj_ + to_matrix_cl(this->adj_); } +}; +} // namespace internal + +/** \ingroup opencl + * Copies the source var that has data stored on the OpenCL device to + * destination var containing Eigen matrices. + * + * @tparam Rows number of compile time rows of the destination matrix + * @tparam Rows number of compile time columns of the destination matrix + * @tparam T type of the matrix or expression on the OpenCL device + * @param a source matrix_cl or expression + * @return var with a copy of the data on the host + */ +template * = nullptr> +inline var_value, Rows, Cols>> from_matrix_cl( + const var_value& a) { + return new internal::op_copy_from_cl_vari(*a.vi_); +} + /** \ingroup opencl * Copies the source Eigen matrix of vars to * the destination matrix that is stored @@ -26,19 +96,17 @@ namespace math { * @param src source Eigen matrix * @return matrix_cl with a copy of the data in the source matrix */ -template * = nullptr> -inline matrix_cl to_matrix_cl(const Eigen::Matrix& src) try { - matrix_cl dst(src.rows(), src.cols()); - if (src.size() == 0) { - return dst; - } - dst.val() = to_matrix_cl(src.val().eval()); - dst.adj() = to_matrix_cl(src.adj().eval()); - return dst; -} catch (const cl::Error& e) { - check_opencl_error("copy var Eigen->(OpenCL)", e); - matrix_cl dst(src.rows(), src.cols()); - return dst; +template * = nullptr> +inline var_value>>> to_matrix_cl( + const T& src) { + // the matrix can go out of scope before chain() is called. So we store a map + // to the data + var* src_array + = ChainableStack::instance_->memalloc_.alloc_array(src.size()); + Eigen::Map> src_stacked(src_array, src.rows(), src.cols()); + src_stacked = src; + return new internal::op_copy_to_cl_vari( + src_stacked.val(), src_stacked.adj()); } /** \ingroup opencl diff --git a/stan/math/opencl/rev/vari.hpp b/stan/math/opencl/rev/vari.hpp new file mode 100644 index 00000000000..b02edc27d19 --- /dev/null +++ b/stan/math/opencl/rev/vari.hpp @@ -0,0 +1,204 @@ +#ifndef STAN_MATH_OPENCL_REV_VARI_HPP +#define STAN_MATH_OPENCL_REV_VARI_HPP +#ifdef STAN_OPENCL + +#include +#include +#include + +namespace stan { +namespace math { +/** + * The variable implementation for `matrix_cl`. + * + * This class is complete (not abstract) and may be used for + * constants. + * + * A variable implementation is constructed with a constant + * value. It also stores the adjoint for storing the partial + * derivative with respect to the root of the derivative tree. + * + */ +template +class vari_value> + : public vari_base, public chainable_alloc { + public: + /** + * The adjoint of this variable, which is the partial derivative + * of this variable with respect to the root variable. + */ + T adj_; + + /** + * The value of this variable. + */ + T val_; + + /** + * Construct a matrix_cl variable implementation from a value. The + * adjoint is initialized to zero. + * + * All constructed variables are added to the stack. Variables + * should be constructed before variables on which they depend + * to insure proper partial derivative propagation. During + * derivative propagation, the chain() method of each variable + * will be called in the reverse order of construction. + * + * @tparam S A `matrix_cl` or kernel generator expression type that is + * convertible to `T` + * @param x Value of the constructed variable. + */ + template * = nullptr> + explicit vari_value(S&& x) + : chainable_alloc(), + adj_(constant(0, x.rows(), x.cols())), + val_(std::forward(x)) { + ChainableStack::instance_->var_stack_.push_back(this); + } + + /** + * Construct a matrix_cl variable implementation from an Eigen value. The + * adjoint is initialized to zero. + * + * All constructed variables are added to the stack. Variables + * should be constructed before variables on which they depend + * to insure proper partial derivative propagation. During + * derivative propagation, the chain() method of each variable + * will be called in the reverse order of construction. + * + * @tparam S A dense Eigen value or expression that has same scalar type as T + * @param x Value of the constructed variable. + */ + template * = nullptr, + require_vt_same* = nullptr> + explicit vari_value(const S& x) + : chainable_alloc(), adj_(constant(0, x.rows(), x.cols())), val_(x) { + ChainableStack::instance_->var_stack_.push_back(this); + } + + /** + * Construct an matrix_cl variable implementation from a value. The + * adjoint is initialized to zero and if `stacked` is `false` this vari + * will be not be put on the var_stack. Instead it will only be put on + * a stack to keep track of whether the adjoint needs to be set to zero. + * + * All constructed variables are added to a stack. Variables + * should be constructed before variables on which they depend + * to insure proper partial derivative propagation. During + * derivative propagation, the chain() method of each variable + * will be called in the reverse order of construction. + * + * @tparam S A `matrix_cl` or kernel generator expression type that is + * convertible to `T` + * @param x Value of the constructed variable. + * @param stacked If false will put this this vari on the nochain stack so + * that its `chain()` method is not called. + */ + template * = nullptr> + vari_value(S&& x, bool stacked) + : chainable_alloc(), + adj_(constant(0, x.rows(), x.cols())), + val_(std::forward(x)) { + if (stacked) { + ChainableStack::instance_->var_stack_.push_back(this); + } else { + ChainableStack::instance_->var_nochain_stack_.push_back(this); + } + } + + /** + * Returns a view into a block of matrix. + * @param row starting row of the block + * @param col starting column of the block + * @param rows number of rows in the block + * @param cols number of columns in the block + * @return block + */ + auto block(int row, int col, int rows, int cols) { + auto&& val_block = stan::math::block(val_, row, col, rows, cols); + auto&& adj_block = stan::math::block(adj_, row, col, rows, cols); + return vari_value>(std::move(val_block), + std::move(adj_block)); + } + + /** + * Return the number of rows for this class's `val_` member + */ + const Eigen::Index rows() const { return val_.rows(); } + /** + * Return the number of columns for this class's `val_` member + */ + const Eigen::Index cols() const { return val_.cols(); } + /** + * Return the size of this class's `val_` member + */ + const Eigen::Index size() const { return val_.size(); } + + virtual void chain() {} + + /** + * Set the adjoint value of this variable to 0. This is used to + * reset adjoints before propagating derivatives again (for + * example in a Jacobian calculation). + */ + inline void set_zero_adjoint() final { adj_ = constant(0, rows(), cols()); } + + /** + * Allocate memory from the underlying memory pool. This memory is + * is managed as a whole externally. + * + * Warning: Classes should not be allocated with this operator + * if they have non-trivial destructors. + * + * @param nbytes Number of bytes to allocate. + * @return Pointer to allocated bytes. + */ + static inline void* operator new(size_t nbytes) { + return ChainableStack::instance_->memalloc_.alloc(nbytes); + } + + /** + * Delete a pointer from the underlying memory pool. + * + * This no-op implementation enables a subclass to throw + * exceptions in its constructor. An exception thrown in the + * constructor of a subclass will result in an error being + * raised, which is in turn caught and calls delete(). + * + * See the discussion of "plugging the memory leak" in: + * http://www.parashift.com/c++-faq/memory-pools.html + */ + static inline void operator delete(void* /* ignore arg */) { /* no op */ + } + + protected: + // to allow access to this constructor from instantinations with different + // template parameters + template + friend class vari_value; + + /** + * Construct a matrix_cl variable implementation from a value and + * adjoint. + * + * This constructor does not add the matrix to any stack! It is intended only + * for views into matrices already on stack. + * + * @param val Value of the constructed variable. + * @param adj Adjoint of the constructed variable. + */ + vari_value(T&& val, T&& adj) + : chainable_alloc(), + adj_(std::forward(adj)), + val_(std::forward(val)) {} + + private: + template + friend class var_value; +}; + +} // namespace math +} // namespace stan + +#endif +#endif diff --git a/stan/math/opencl/value_type.hpp b/stan/math/opencl/value_type.hpp index dbbc7750655..2effcf27184 100644 --- a/stan/math/opencl/value_type.hpp +++ b/stan/math/opencl/value_type.hpp @@ -12,8 +12,7 @@ namespace stan { * Return the value type of an OpenCL matrix. */ template -struct value_type> { +struct value_type> { using type = typename std::decay_t::Scalar; }; } // namespace stan diff --git a/stan/math/prim/fun/MatrixExponential.h b/stan/math/prim/fun/MatrixExponential.h index 85e9f8274fa..74ae2a47ee2 100644 --- a/stan/math/prim/fun/MatrixExponential.h +++ b/stan/math/prim/fun/MatrixExponential.h @@ -185,7 +185,7 @@ namespace Eigen { { template static void run(const MatrixType& arg, MatrixType& U, MatrixType& V, int& squarings, - T scalar_type) + T scalar_type) { using std::frexp; using std::pow; diff --git a/stan/math/prim/meta.hpp b/stan/math/prim/meta.hpp index 4542d65cacc..ef1128d36b4 100644 --- a/stan/math/prim/meta.hpp +++ b/stan/math/prim/meta.hpp @@ -187,9 +187,12 @@ #include #include #include +#include #include #include +#include #include +#include #include #include #include diff --git a/stan/math/prim/meta/is_eigen_dense_base.hpp b/stan/math/prim/meta/is_eigen_dense_base.hpp new file mode 100644 index 00000000000..3ce1dbe8ebc --- /dev/null +++ b/stan/math/prim/meta/is_eigen_dense_base.hpp @@ -0,0 +1,31 @@ +#ifndef STAN_MATH_PRIM_META_IS_EIGEN_DENSE_BASE_HPP +#define STAN_MATH_PRIM_META_IS_EIGEN_DENSE_BASE_HPP + +#include +#include +#include +#include +#include + +namespace stan { + +/** + * Checks whether type T is derived from Eigen::DenseBase. + * If true this will have a static member function named value with a type + * of true, else value is false. + * @tparam T Type to check if it is derived from `DenseBase` + * @tparam Enable used for SFINAE deduction. + * @ingroup type_trait + */ +template +struct is_eigen_dense_base + : bool_constant::value> {}; + +STAN_ADD_REQUIRE_UNARY(eigen_dense_base, is_eigen_dense_base, + require_eigens_types); +STAN_ADD_REQUIRE_CONTAINER(eigen_dense_base, is_eigen_dense_base, + require_eigens_types); + +} // namespace stan + +#endif diff --git a/stan/math/prim/meta/is_eigen_matrix_base.hpp b/stan/math/prim/meta/is_eigen_matrix_base.hpp index f6d47c0b2e1..9e4e6e9465e 100644 --- a/stan/math/prim/meta/is_eigen_matrix_base.hpp +++ b/stan/math/prim/meta/is_eigen_matrix_base.hpp @@ -13,7 +13,7 @@ namespace stan { * Checks whether type T is derived from Eigen::MatrixBase. * If true this will have a static member function named value with a type * of true, else value is false. - * @tparam T Type to check if it is derived from `EigenBase` + * @tparam T Type to check if it is derived from `MatrixBase` * @tparam Enable used for SFINAE deduction. * @ingroup type_trait */ diff --git a/stan/math/prim/meta/is_eigen_sparse_base.hpp b/stan/math/prim/meta/is_eigen_sparse_base.hpp new file mode 100644 index 00000000000..2c1028730bf --- /dev/null +++ b/stan/math/prim/meta/is_eigen_sparse_base.hpp @@ -0,0 +1,32 @@ +#ifndef STAN_MATH_PRIM_META_IS_EIGEN_SPARSE_BASE_HPP +#define STAN_MATH_PRIM_META_IS_EIGEN_SPARSE_BASE_HPP + +#include +#include +#include +#include +#include + +namespace stan { + +/** + * Checks whether type T is derived from Eigen::SparseMatrixBase. + * If true this will have a static member function named value with a type + * of true, else value is false. + * @tparam T Type to check if it is derived from `SparseMatrixBase` + * @tparam Enable used for SFINAE deduction. + * @ingroup type_trait + */ +template +struct is_eigen_sparse_base + : bool_constant< + is_base_pointer_convertible::value> {}; + +STAN_ADD_REQUIRE_UNARY(eigen_sparse_base, is_eigen_sparse_base, + require_eigens_types); +STAN_ADD_REQUIRE_CONTAINER(eigen_sparse_base, is_eigen_sparse_base, + require_eigens_types); + +} // namespace stan + +#endif diff --git a/stan/math/prim/meta/is_plain_type.hpp b/stan/math/prim/meta/is_plain_type.hpp new file mode 100644 index 00000000000..4588d08a7d6 --- /dev/null +++ b/stan/math/prim/meta/is_plain_type.hpp @@ -0,0 +1,20 @@ +#ifndef STAN_MATH_PRIM_META_IS_PLAIN_TYPE_HPP +#define STAN_MATH_PRIM_META_IS_PLAIN_TYPE_HPP + +#include +#include +#include + +namespace stan { +/** \ingroup type_trait + * Checks whether the template type `T` is an assignable type. This is used + * to detect whether a type is an Eigen matrix expression. + */ +template +using is_plain_type = std::is_same, plain_type_t>; + +STAN_ADD_REQUIRE_UNARY(plain_type, is_plain_type, require_eigens_types); +STAN_ADD_REQUIRE_UNARY_INNER(plain_type, is_plain_type, require_eigens_types); + +} // namespace stan +#endif diff --git a/stan/math/prim/meta/ref_type.hpp b/stan/math/prim/meta/ref_type.hpp index 1660f9eb183..03a83a76607 100644 --- a/stan/math/prim/meta/ref_type.hpp +++ b/stan/math/prim/meta/ref_type.hpp @@ -73,12 +73,13 @@ struct ref_type_for_opencl { = std::conditional_t::value, T_val, const T&>; // Setting Outer stride of Ref to 0 (default) won't actually check that // expression has contiguous outer stride. Instead we need to check that - // evaluator flags contain LinearAccessBit. + // evaluator flags contain LinearAccessBit and PacketAccessBit. using type = std::conditional_t< Eigen::internal::traits>>:: template match::MatchAtCompileTime + && (Eigen::internal::evaluator::Flags & Eigen::LinearAccessBit) && (Eigen::internal::evaluator::Flags - & Eigen::LinearAccessBit), + & Eigen::PacketAccessBit), T_optionally_ref, T_plain_col_major>; }; diff --git a/stan/math/rev/core/chainablestack.hpp b/stan/math/rev/core/chainablestack.hpp index 4785afe1146..bb28eb1772d 100644 --- a/stan/math/rev/core/chainablestack.hpp +++ b/stan/math/rev/core/chainablestack.hpp @@ -4,11 +4,8 @@ #include namespace stan { namespace math { - -template -class vari_value; -class vari_base; class chainable_alloc; +class vari_base; using ChainableStack = AutodiffStackSingleton; } // namespace math diff --git a/stan/math/rev/core/grad.hpp b/stan/math/rev/core/grad.hpp index efe8fee1e87..570b86d18d6 100644 --- a/stan/math/rev/core/grad.hpp +++ b/stan/math/rev/core/grad.hpp @@ -11,6 +11,26 @@ namespace stan { namespace math { +/** + * Compute the gradient for all variables starting from the end of the AD tape. + * This function does not recover memory. The chain + * rule is applied working down the stack from the last vari created on the + * AD tape and then calling each vari's `chain()` method in turn. + * + *

This function computes a nested gradient only going back as far + * as the last nesting. + * + *

This function does not recover any memory from the computation. + * + */ +static void grad() { + size_t end = ChainableStack::instance_->var_stack_.size(); + size_t beginning = empty_nested() ? 0 : end - nested_size(); + for (size_t i = end; i-- > beginning;) { + ChainableStack::instance_->var_stack_[i]->chain(); + } +} + /** * Compute the gradient for all variables starting from the * specified root variable implementation. Does not recover @@ -30,11 +50,7 @@ namespace math { template static void grad(Vari* vi) { vi->init_dependent(); - size_t end = ChainableStack::instance_->var_stack_.size(); - size_t beginning = empty_nested() ? 0 : end - nested_size(); - for (size_t i = end; i-- > beginning;) { - ChainableStack::instance_->var_stack_[i]->chain(); - } + grad(); } } // namespace math diff --git a/stan/math/rev/core/operator_divide_equal.hpp b/stan/math/rev/core/operator_divide_equal.hpp index 23687c7011e..5374ad3467b 100644 --- a/stan/math/rev/core/operator_divide_equal.hpp +++ b/stan/math/rev/core/operator_divide_equal.hpp @@ -8,15 +8,13 @@ namespace stan { namespace math { template -inline var_value& var_value>::operator/=( - const var_value& b) { +inline var_value& var_value::operator/=(const var_value& b) { vi_ = new internal::divide_vv_vari(vi_, b.vi_); return *this; } template -inline var_value& var_value>::operator/=( - T b) { +inline var_value& var_value::operator/=(T b) { if (b == 1.0) { return *this; } diff --git a/stan/math/rev/core/operator_minus_equal.hpp b/stan/math/rev/core/operator_minus_equal.hpp index f72f7d7412d..2a9d33dea2b 100644 --- a/stan/math/rev/core/operator_minus_equal.hpp +++ b/stan/math/rev/core/operator_minus_equal.hpp @@ -9,15 +9,13 @@ namespace stan { namespace math { template -inline var_value& var_value>::operator-=( - const var_value& b) { +inline var_value& var_value::operator-=(const var_value& b) { vi_ = new internal::subtract_vv_vari(vi_, b.vi_); return *this; } template -inline var_value& var_value>::operator-=( - T b) { +inline var_value& var_value::operator-=(T b) { if (b == 0.0) { return *this; } diff --git a/stan/math/rev/core/operator_multiply_equal.hpp b/stan/math/rev/core/operator_multiply_equal.hpp index 61c13263c19..45035a6acfe 100644 --- a/stan/math/rev/core/operator_multiply_equal.hpp +++ b/stan/math/rev/core/operator_multiply_equal.hpp @@ -9,15 +9,13 @@ namespace stan { namespace math { template -inline var_value& var_value>::operator*=( - const var_value& b) { +inline var_value& var_value::operator*=(const var_value& b) { vi_ = new internal::multiply_vv_vari(vi_, b.vi_); return *this; } template -inline var_value& var_value>::operator*=( - T b) { +inline var_value& var_value::operator*=(T b) { if (b == 1.0) { return *this; } diff --git a/stan/math/rev/core/operator_plus_equal.hpp b/stan/math/rev/core/operator_plus_equal.hpp index 32d41b30772..769a96d0db1 100644 --- a/stan/math/rev/core/operator_plus_equal.hpp +++ b/stan/math/rev/core/operator_plus_equal.hpp @@ -9,15 +9,13 @@ namespace stan { namespace math { template -inline var_value& var_value>::operator+=( - const var_value& b) { +inline var_value& var_value::operator+=(const var_value& b) { vi_ = new internal::add_vv_vari(vi_, b.vi_); return *this; } template -inline var_value& var_value>::operator+=( - T b) { +inline var_value& var_value::operator+=(T b) { if (b == 0.0) { return *this; } diff --git a/stan/math/rev/core/set_zero_all_adjoints_nested.hpp b/stan/math/rev/core/set_zero_all_adjoints_nested.hpp index 34dc367556c..22c74b5ccc8 100644 --- a/stan/math/rev/core/set_zero_all_adjoints_nested.hpp +++ b/stan/math/rev/core/set_zero_all_adjoints_nested.hpp @@ -23,14 +23,15 @@ static void set_zero_all_adjoints_nested() { "empty_nested() must be false before calling" " set_zero_all_adjoints_nested()"); } - size_t start1 = ChainableStack::instance_->nested_var_stack_sizes_.back(); + const size_t start1 + = ChainableStack::instance_->nested_var_stack_sizes_.back(); // avoid wrap with unsigned when start1 == 0 for (size_t i = (start1 == 0U) ? 0U : (start1 - 1); i < ChainableStack::instance_->var_stack_.size(); ++i) { ChainableStack::instance_->var_stack_[i]->set_zero_adjoint(); } - size_t start2 + const size_t start2 = ChainableStack::instance_->nested_var_nochain_stack_sizes_.back(); for (size_t i = (start2 == 0U) ? 0U : (start2 - 1); i < ChainableStack::instance_->var_nochain_stack_.size(); ++i) { diff --git a/stan/math/rev/core/var.hpp b/stan/math/rev/core/var.hpp index a7dbaf7b149..b7b453d7a28 100644 --- a/stan/math/rev/core/var.hpp +++ b/stan/math/rev/core/var.hpp @@ -8,6 +8,9 @@ #include #include #include +#ifdef STAN_OPENCL +#include +#endif namespace stan { namespace math { @@ -30,13 +33,20 @@ static void grad(Vari* vi); * var values objects. * @tparam T An Floating point type. */ -template -class var_value {}; - template -class var_value> { +class var_value { + static_assert( + is_plain_type::value, + "The template for this var is an" + " expression but a var_value's inner type must be assignable such as" + " a double, Eigen::Matrix, or Eigen::Array"); + static_assert( + std::is_floating_point>::value, + "The template for must be a floating point or a container holding" + " floating point types"); + public: - using value_type = std::decay_t; // Numeric type in vari_value. + using value_type = std::decay_t; // type in vari_value. using vari_type = vari_value; // Type of underlying vari impl. /** @@ -57,7 +67,7 @@ class var_value> { * @return true if this variable does not yet have * a defined variable. */ - bool is_uninitialized() { return (vi_ == nullptr); } + inline bool is_uninitialized() { return (vi_ == nullptr); } /** * Construct a variable for later assignment. @@ -76,14 +86,13 @@ class var_value> { * @param x Value of the variable. */ template * = nullptr> - var_value(const S& x) : vi_(new vari_type(x, false)) {} // NOLINT + var_value(S&& x) : vi_(new vari_type(std::forward(x), false)) {} // NOLINT /** * Construct a variable from a pointer to a variable implementation. * @param vi A vari_value pointer. */ - var_value(vari_value* vi) // NOLINT - : vi_(vi) {} + var_value(vari_type* vi) : vi_(vi) {} // NOLINT /** * Return a constant reference to the value of this variable. @@ -120,10 +129,15 @@ class var_value> { * The grad() function does not recover memory. In Stan * 2.4 and earlier, this function did recover memory. * + * @tparam CheckContainer Not set by user. The default value of value_type + * is used to require that grad is only available for scalar `var_value` + * types. * @param x Vector of independent variables. * @param g Gradient vector of partial derivatives of this * variable with respect to x. */ + template * = nullptr> inline void grad(std::vector>& x, std::vector& g) { stan::math::grad(vi_); g.resize(x.size()); @@ -136,9 +150,16 @@ class var_value> { * Compute the gradient of this (dependent) variable with respect * to all (independent) variables. * + * @tparam CheckContainer Not set by user. The default value of value_type + * is used to require that grad is only available for scalar `var_value` + * types. * The grad() function does not recover memory. */ - void grad() { stan::math::grad(vi_); } + template * = nullptr> + void grad() { + stan::math::grad(vi_); + } // POINTER OVERRIDES diff --git a/stan/math/rev/core/vari.hpp b/stan/math/rev/core/vari.hpp index 6515174ae8b..b17e532870c 100644 --- a/stan/math/rev/core/vari.hpp +++ b/stan/math/rev/core/vari.hpp @@ -11,7 +11,7 @@ namespace stan { namespace math { // forward declaration of var -template +template class var_value; /** * Abstract base class that all `vari_value` and it's derived classes inherit. @@ -33,8 +33,10 @@ class vari_base { virtual ~vari_base() noexcept {} }; +template +class vari_value; /** - * The variable implementation base class. + * The variable implementation for floating point types. * * This class is complete (not abstract) and may be used for * constants. @@ -45,24 +47,18 @@ class vari_base { * */ template -class vari_value::value>> - : public vari_base { - private: - template - friend class var_value; - +class vari_value> : public vari_base { public: - using Scalar = T; - using value_type = Scalar; + using value_type = std::decay_t; /** * The value of this variable. */ - const Scalar val_; + const value_type val_; /** * The adjoint of this variable, which is the partial derivative * of this variable with respect to the root variable. */ - Scalar adj_; + value_type adj_; /** * Construct a variable implementation from a value. The @@ -77,8 +73,7 @@ class vari_value::value>> * @tparam S a floating point type. * @param x Value of the constructed variable. */ - template ::value>* = nullptr> + template * = nullptr> vari_value(S x) noexcept : val_(x), adj_(0.0) { // NOLINT ChainableStack::instance_->var_stack_.emplace_back(this); } @@ -88,9 +83,7 @@ class vari_value::value>> * adjoint is initialized to zero and if `stacked` is `false` this vari * will be not be put on the var_stack. Instead it will only be put on * a stack to keep track of whether the adjoint needs to be set to zero. - * - * All constructed variables are added to a stack. Variables - * should be constructed before variables on which they depend + * Variables should be constructed before variables on which they depend * to insure proper partial derivative propagation. During * derivative propagation, the chain() method of each variable * will be called in the reverse order of construction. @@ -100,8 +93,7 @@ class vari_value::value>> * @param stacked If false will put this this vari on the nochain stack so * that its `chain()` method is not called. */ - template ::value>* = nullptr> + template * = nullptr> vari_value(S x, bool stacked) noexcept : val_(x), adj_(0.0) { if (stacked) { ChainableStack::instance_->var_stack_.emplace_back(this); @@ -110,15 +102,6 @@ class vari_value::value>> } } - /** - * Constructor from vari_value - * @tparam S A floating point type - * @param x A vari_value - */ - vari_value(const vari_value& x) noexcept : val_(x.val_), adj_(x.adj_) { - ChainableStack::instance_->var_stack_.emplace_back(this); - } - ~vari_value() = default; inline void chain() {} @@ -179,6 +162,372 @@ class vari_value::value>> static inline void operator delete( void* /* ignore arg */) noexcept { /* no op */ } + + private: + template + friend class var_value; +}; + +// For backwards compatability the default is double +using vari = vari_value; + +/** + * The variable implementation for Eigen dense matrix types. + * + * This class is complete (not abstract) and may be used for + * constants. + * + * A variable implementation is constructed with a constant + * value. It also stores the adjoint for storing the partial + * derivative with respect to the root of the derivative tree. + * + */ +template +class vari_value> : public vari_base { + public: + static_assert( + is_plain_type::value, + "The template for this var is an" + " expression but a var_value's inner type must be assignable such as" + " a double, Eigen::Matrix, or Eigen::Array"); + + /** + * `PlainObject` represents a user constructible type such as Matrix or Array + */ + using PlainObject = plain_type_t; + using value_type = PlainObject; // The underlying type for this class + using eigen_scalar = value_type_t; // A floating point type + /** + * Maps for adj_ and val_ + */ + using eigen_map = Eigen::Map; + using vari_type = vari_value>; + eigen_scalar* val_mem_; // Pointer to memory allocated on the stack for val_ + eigen_scalar* adj_mem_; // Pointer to memory allocated on the stack for adj_ + /** + * Number of rows known at compile time + */ + static constexpr int RowsAtCompileTime = PlainObject::RowsAtCompileTime; + /** + * Number of columns known at compile time + */ + static constexpr int ColsAtCompileTime = PlainObject::ColsAtCompileTime; + + /** + * The value of this variable. + */ + const eigen_map val_; + + /** + * The adjoint of this variable, which is the partial derivative + * of this variable with respect to the root variable. + */ + eigen_map adj_; + + /** + * Construct a dense Eigen variable implementation from a value. The + * adjoint is initialized to zero. + * + * All constructed variables are added to the stack. Variables + * should be constructed before variables on which they depend + * to insure proper partial derivative propagation. During + * derivative propagation, the chain() method of each variable + * will be called in the reverse order of construction. + * + * @tparam S A dense Eigen type that is convertible to `value_type` + * @param x Value of the constructed variable. + */ + template * = nullptr> + explicit vari_value(const S& x) + : val_mem_(ChainableStack::instance_->memalloc_.alloc_array( + x.size())), + adj_mem_(ChainableStack::instance_->memalloc_.alloc_array( + x.size())), + val_(eigen_map(val_mem_, x.rows(), x.cols()) = x), + adj_(eigen_map(adj_mem_, x.rows(), x.cols()).setZero()) { + ChainableStack::instance_->var_stack_.push_back(this); + } + + /** + * Construct a dense Eigen variable implementation from a value. The + * adjoint is initialized to zero and if `stacked` is `false` this vari + * will be not be put on the var_stack. Instead it will only be put on + * a stack to keep track of whether the adjoint needs to be set to zero. + * Variables should be constructed before variables on which they depend + * to insure proper partial derivative propagation. During + * derivative propagation, the chain() method of each variable + * will be called in the reverse order of construction. + * + * @tparam S A dense Eigen type that is convertible to `value_type` + * @param x Value of the constructed variable. + * @param stacked If false will put this this vari on the nochain stack so + * that its `chain()` method is not called. + */ + template * = nullptr> + vari_value(const S& x, bool stacked) + : val_mem_(ChainableStack::instance_->memalloc_.alloc_array( + x.size())), + adj_mem_(ChainableStack::instance_->memalloc_.alloc_array( + x.size())), + val_(eigen_map(val_mem_, x.rows(), x.cols()) = x), + adj_(eigen_map(adj_mem_, x.rows(), x.cols()).setZero()) { + if (stacked) { + ChainableStack::instance_->var_stack_.emplace_back(this); + } else { + ChainableStack::instance_->var_nochain_stack_.emplace_back(this); + } + } + + /** + * Return the number of rows for this class's `val_` member + */ + const Eigen::Index rows() const { return val_.rows(); } + /** + * Return the number of columns for this class's `val_` member + */ + const Eigen::Index cols() const { return val_.rows(); } + /** + * Return the size of this class's `val_` member + */ + const Eigen::Index size() const { return val_.size(); } + + virtual void chain() {} + /** + * Initialize the adjoint for this (dependent) variable to 1. + * This operation is applied to the dependent variable before + * propagating derivatives, setting the derivative of the + * result with respect to itself to be 1. + */ + inline void init_dependent() { adj_.setOnes(); } + + /** + * Set the adjoint value of this variable to 0. This is used to + * reset adjoints before propagating derivatives again (for + * example in a Jacobian calculation). + */ + inline void set_zero_adjoint() final { adj_.setZero(); } + + /** + * Insertion operator for vari. Prints the current value and + * the adjoint value. + * + * @param os [in, out] ostream to modify + * @param v [in] vari object to print. + * + * @return The modified ostream. + */ + friend std::ostream& operator<<(std::ostream& os, const vari_value* v) { + return os << "val: \n" << v->val_ << " \nadj: \n" << v->adj_; + } + + /** + * Allocate memory from the underlying memory pool. This memory is + * is managed as a whole externally. + * + * Warning: Classes should not be allocated with this operator + * if they have non-trivial destructors. + * + * @param nbytes Number of bytes to allocate. + * @return Pointer to allocated bytes. + */ + static inline void* operator new(size_t nbytes) { + return ChainableStack::instance_->memalloc_.alloc(nbytes); + } + + /** + * Delete a pointer from the underlying memory pool. + * + * This no-op implementation enables a subclass to throw + * exceptions in its constructor. An exception thrown in the + * constructor of a subclass will result in an error being + * raised, which is in turn caught and calls delete(). + * + * See the discussion of "plugging the memory leak" in: + * http://www.parashift.com/c++-faq/memory-pools.html + */ + static inline void operator delete(void* /* ignore arg */) { /* no op */ + } + + private: + template + friend class var_value; + template + friend class vari_value; +}; + +/** + * The variable implementation for Eigen sparse matrix types. + * + * This class is complete (not abstract) and may be used for + * constants. + * + * A variable implementation is constructed with a constant + * value. It also stores the adjoint for storing the partial + * derivative with respect to the root of the derivative tree. + * + */ +template +class vari_value> : public vari_base, + chainable_alloc { + public: + using PlainObject = plain_type_t; // Base type of Eigen class + using value_type = PlainObject; // vari's adj_ and val_ member type + /** + * Rows at compile time + */ + static constexpr int RowsAtCompileTime = T::RowsAtCompileTime; + /** + * Columns at compile time + */ + static constexpr int ColsAtCompileTime = T::ColsAtCompileTime; + + /** + * The adjoint of this variable, which is the partial derivative + * of this variable with respect to the root variable. + */ + PlainObject adj_; + + /** + * The value of this variable. + */ + const PlainObject val_; + + /** + * Construct a variable implementation from a value. The + * adjoint is initialized to zero. + * + * All constructed variables are added to the stack. For a sparse eigen matrix + * this includes the nozero values as well the inner and outer indices. + * Variables should be constructed before variables on which they depend + * to insure proper partial derivative propagation. During + * derivative propagation, the chain() method of each variable + * will be called in the reverse order of construction. + * + * @tparam S A sparse Eigen type that is convertible to `value_type` + * @param x Value of the constructed variable. + */ + template * = nullptr> + explicit vari_value(S&& x) + : adj_(x), val_(std::forward(x)), chainable_alloc() { + this->set_zero_adjoint(); + ChainableStack::instance_->var_stack_.push_back(this); + } + /** + * Construct an sparse Eigen variable implementation from a value. The + * adjoint is initialized to zero and if `stacked` is `false` this vari + * will be not be put on the var_stack. Instead it will only be put on + * a stack to keep track of whether the adjoint needs to be set to zero. + * + * All constructed variables are added to a stack. Variables + * should be constructed before variables on which they depend + * to insure proper partial derivative propagation. During + * derivative propagation, the chain() method of each variable + * will be called in the reverse order of construction. + * + * @tparam S A sparse Eigen type that is convertible to `value_type` + * @param x Value of the constructed variable. + * @param stacked If false will put this this vari on the nochain stack so + * that its `chain()` method is not called. + */ + template * = nullptr> + vari_value(S&& x, bool stacked) + : adj_(x), val_(std::forward(x)), chainable_alloc() { + this->set_zero_adjoint(); + if (stacked) { + ChainableStack::instance_->var_stack_.push_back(this); + } else { + ChainableStack::instance_->var_nochain_stack_.push_back(this); + } + } + + /** + * Return the number of rows for this class's `val_` member + */ + const Eigen::Index rows() const { return val_.rows(); } + /** + * Return the number of columns for this class's `val_` member + */ + const Eigen::Index cols() const { return val_.cols(); } + /** + * Return the size of this class's `val_` member + */ + const Eigen::Index size() const { return val_.size(); } + + void chain() {} + /** + * Initialize the adjoint for this (dependent) variable to 1. + * This operation is applied to the dependent variable before + * propagating derivatives, setting the derivative of the + * result with respect to itself to be 1. + */ + inline void init_dependent() { + for (int k = 0; k < adj_.outerSize(); ++k) { + for (typename PlainObject::InnerIterator it(adj_, k); it; ++it) { + it.valueRef() = 1.0; + } + } + } + + /** + * Set the adjoint value of this variable to 0. This is used to + * reset adjoints before propagating derivatives again (for + * example in a Jacobian calculation). + */ + inline void set_zero_adjoint() noexcept final { + for (int k = 0; k < adj_.outerSize(); ++k) { + for (typename PlainObject::InnerIterator it(adj_, k); it; ++it) { + it.valueRef() = 0.0; + } + } + } + + /** + * Insertion operator for vari. Prints the current value and + * the adjoint value. + * + * @param os [in, out] ostream to modify + * @param v [in] vari object to print. + * + * @return The modified ostream. + */ + friend std::ostream& operator<<(std::ostream& os, const vari_value* v) { + return os << "val: \n" << v->val_ << " \nadj: \n" << v->adj_; + } + + /** + * Allocate memory from the underlying memory pool. This memory is + * is managed as a whole externally. + * + * Warning: Classes should not be allocated with this operator + * if they have non-trivial destructors. + * + * @param nbytes Number of bytes to allocate. + * @return Pointer to allocated bytes. + */ + static inline void* operator new(size_t nbytes) noexcept { + return ChainableStack::instance_->memalloc_.alloc(nbytes); + } + + /** + * Delete a pointer from the underlying memory pool. + * + * This no-op implementation enables a subclass to throw + * exceptions in its constructor. An exception thrown in the + * constructor of a subclass will result in an error being + * raised, which is in turn caught and calls delete(). + * + * See the discussion of "plugging the memory leak" in: + * http://www.parashift.com/c++-faq/memory-pools.html + */ + static inline void operator delete( + void* /* ignore arg */) noexcept { /* no op */ + } + + private: + template + friend class var_value; + template + friend class vari_value; }; // For backwards compatability the default is double diff --git a/stan/math/rev/fun/lambert_w.hpp b/stan/math/rev/fun/lambert_w.hpp index 0f2455f3cc6..b1c845585cf 100644 --- a/stan/math/rev/fun/lambert_w.hpp +++ b/stan/math/rev/fun/lambert_w.hpp @@ -5,6 +5,7 @@ #include #include #include +#include #include namespace stan { diff --git a/test/unit/math/mix/meta/is_eigen_sparse_base_test.cpp b/test/unit/math/mix/meta/is_eigen_sparse_base_test.cpp new file mode 100644 index 00000000000..cc513e83160 --- /dev/null +++ b/test/unit/math/mix/meta/is_eigen_sparse_base_test.cpp @@ -0,0 +1,58 @@ +#include +#include +#include +#include + +TEST(MathMetaPrim, is_eigen_sparse_base_hierarchy_tests) { + using Eigen::Array; + using Eigen::Matrix; + using stan::is_eigen_sparse_base; + using stan::math::test::all_eigen_dense; + all_eigen_dense(); + all_eigen_dense(); + all_eigen_dense(); + all_eigen_dense(); + all_eigen_dense(); + all_eigen_dense(); +} + +TEST(MathMetaPrim, is_eigen_sparse_base_sparse_tests) { + using Eigen::Array; + using Eigen::Matrix; + using stan::is_eigen_sparse_base; + using stan::math::test::all_eigen_sparse; + all_eigen_sparse(); +} + +TEST(MathMetaPrim, is_eigen_sparse_base_decomp_tests) { + using Eigen::Array; + using Eigen::Matrix; + using stan::is_eigen_sparse_base; + using stan::math::test::all_eigen_dense_decomp; + all_eigen_dense_decomp(); +} + +TEST(MathMetaPrim, is_eigen_sparse_base_expr_tests) { + using Eigen::Array; + using Eigen::Matrix; + using stan::is_eigen_sparse_base; + using stan::math::test::all_eigen_dense_exprs; + all_eigen_dense_exprs(); + all_eigen_dense_exprs(); + all_eigen_dense_exprs(); + all_eigen_dense_exprs(); + all_eigen_dense_exprs(); + all_eigen_dense_exprs(); +} diff --git a/test/unit/math/opencl/kernel_generator/constant_test.cpp b/test/unit/math/opencl/kernel_generator/constant_test.cpp new file mode 100644 index 00000000000..a14dbabff5b --- /dev/null +++ b/test/unit/math/opencl/kernel_generator/constant_test.cpp @@ -0,0 +1,17 @@ +#ifdef STAN_OPENCL + +#include +#include +#include +#include +#include +#include + +TEST(KernelGenerator, constant_test) { + stan::math::matrix_cl m1_cl(stan::math::constant(1.2, 3, 4)); + + Eigen::MatrixXd res = stan::math::from_matrix_cl(m1_cl); + EXPECT_MATRIX_EQ(res, Eigen::MatrixXd::Constant(3, 4, 1.2)); +} + +#endif diff --git a/test/unit/math/opencl/rev/copy_test.cpp b/test/unit/math/opencl/rev/copy_test.cpp index c5e43d49db3..b55432ad6bc 100644 --- a/test/unit/math/opencl/rev/copy_test.cpp +++ b/test/unit/math/opencl/rev/copy_test.cpp @@ -1,5 +1,6 @@ #ifdef STAN_OPENCL #include +#include #include #include #include @@ -14,7 +15,8 @@ TEST(MathMatrixRevCL, matrix_cl_vector_copy) { // vector stan::math::matrix_cl d11_cl(3, 1); stan::math::matrix_cl d111_cl(3, 1); - EXPECT_NO_THROW(d11_cl = stan::math::to_matrix_cl(d1_cpu)); + d11_cl.val() = stan::math::to_matrix_cl(d1_cpu.val()); + d11_cl.adj() = stan::math::to_matrix_cl(d1_cpu.adj()); EXPECT_NO_THROW(d11_cl.adj() = stan::math::copy_cl(d11_cl.val())); EXPECT_NO_THROW(d111_cl = stan::math::copy_cl(d11_cl)); EXPECT_NO_THROW(d1_a_cpu = stan::math::from_matrix_cl(d11_cl)); @@ -39,7 +41,8 @@ TEST(MathMatrixRevCL, matrix_cl_matrix_copy) { stan::math::matrix_cl d000_cl; stan::math::matrix_cl d22_cl(2, 3); stan::math::matrix_cl d222_cl(2, 3); - EXPECT_NO_THROW(d22_cl = stan::math::to_matrix_cl(d2_cpu)); + d22_cl.val() = stan::math::to_matrix_cl(d2_cpu.val()); + d22_cl.adj() = stan::math::to_matrix_cl(d2_cpu.adj()); EXPECT_NO_THROW(d222_cl = stan::math::copy_cl(d22_cl)); EXPECT_NO_THROW(d2_a_cpu = stan::math::from_matrix_cl(d22_cl.val())); EXPECT_NO_THROW(d2_b_cpu = stan::math::from_matrix_cl(d222_cl.val())); @@ -55,8 +58,8 @@ TEST(MathMatrixRevCL, matrix_cl_matrix_copy) { EXPECT_EQ(4, d2_b_cpu(1, 0)); EXPECT_EQ(5, d2_b_cpu(1, 1)); EXPECT_EQ(6, d2_b_cpu(1, 2)); - d00_cl = stan::math::to_matrix_cl(d0_cpu); - EXPECT_NO_THROW(d00_cl = stan::math::to_matrix_cl(d0_cpu)); + d00_cl.val() = stan::math::to_matrix_cl(d0_cpu.val()); + d00_cl.adj() = stan::math::to_matrix_cl(d0_cpu.adj()); EXPECT_NO_THROW(d0_cpu = stan::math::from_matrix_cl(d00_cl.val())); EXPECT_NO_THROW(d000_cl = stan::math::copy_cl(d00_cl)); } @@ -150,4 +153,44 @@ TEST(MathMatrixRevCL, matrix_cl_pack_unpack_copy_exception) { std::invalid_argument); */ } + +TEST(VariCL, var_matrix_to_matrix_cl) { + using stan::math::var_value; + Eigen::MatrixXd vals(2, 3); + vals << 1, 2, 3, 4, 5, 6; + var_value a(vals); + var_value> a_cl = stan::math::to_matrix_cl(a); + EXPECT_MATRIX_EQ(from_matrix_cl(a_cl.val()), vals); + a_cl.adj() = stan::math::constant(1, 2, 3); + a_cl.vi_->chain(); + EXPECT_MATRIX_EQ(a.adj(), Eigen::MatrixXd::Constant(2, 3, 1)); +} + +TEST(VariCL, matrix_var_to_matrix_cl) { + using stan::math::var_value; + Eigen::MatrixXd vals(2, 3); + vals << 1, 2, 3, 4, 5, 6; + stan::math::matrix_v vars = vals; + var_value> a_cl + = stan::math::to_matrix_cl(vars); + EXPECT_MATRIX_EQ(from_matrix_cl(a_cl.val()), vals); + a_cl.adj() = stan::math::constant(1, 2, 3); + a_cl.vi_->chain(); + EXPECT_MATRIX_EQ(vars.adj(), Eigen::MatrixXd::Constant(2, 3, 1)); +} + +TEST(VariCL, from_matrix_cl) { + using stan::math::var_value; + Eigen::MatrixXd vals(2, 3); + vals << 1, 2, 3, 4, 5, 6; + stan::math::matrix_cl vals_cl(vals); + var_value> a_cl(vals_cl); + var_value a = stan::math::from_matrix_cl(a_cl); + EXPECT_MATRIX_EQ(a.val(), vals); + a.adj() = Eigen::MatrixXd::Constant(2, 3, 1); + a.vi_->chain(); + EXPECT_MATRIX_EQ(from_matrix_cl(a_cl.adj()), + Eigen::MatrixXd::Constant(2, 3, 1)); +} + #endif diff --git a/test/unit/math/opencl/rev/sub_block_test.cpp b/test/unit/math/opencl/rev/sub_block_test.cpp index f492b26de63..32ccc28cded 100644 --- a/test/unit/math/opencl/rev/sub_block_test.cpp +++ b/test/unit/math/opencl/rev/sub_block_test.cpp @@ -64,7 +64,8 @@ TEST(MathMatrixRevCL, sub_block_triangular) { stan::math::matrix_cl b_cl(b); Eigen::Matrix c; - a_cl = stan::math::to_matrix_cl(a); + a_cl.val() = stan::math::to_matrix_cl(a.val()); + a_cl.adj() = stan::math::to_matrix_cl(a.adj()); a_cl.view(stan::math::matrix_cl_view::Lower); b_cl.view(stan::math::matrix_cl_view::Lower); a_cl.sub_block(b_cl, 0, 1, 0, 1, 2, 2); @@ -80,7 +81,8 @@ TEST(MathMatrixRevCL, sub_block_triangular) { EXPECT_EQ(c(2, 1), 0); EXPECT_EQ(c(2, 2), 0); - a_cl = stan::math::to_matrix_cl(a); + a_cl.val() = stan::math::to_matrix_cl(a.val()); + a_cl.adj() = stan::math::to_matrix_cl(a.adj()); a_cl.view(stan::math::matrix_cl_view::Lower); b_cl.view(stan::math::matrix_cl_view::Lower); a_cl.sub_block(b_cl, 1, 0, 1, 1, 2, 2); @@ -96,7 +98,8 @@ TEST(MathMatrixRevCL, sub_block_triangular) { EXPECT_EQ(c(2, 1), 1); EXPECT_EQ(c(2, 2), 1); - a_cl = stan::math::to_matrix_cl(a); + a_cl.val() = stan::math::to_matrix_cl(a.val()); + a_cl.adj() = stan::math::to_matrix_cl(a.adj()); a_cl.view(stan::math::matrix_cl_view::Lower); b_cl.view(stan::math::matrix_cl_view::Upper); a_cl.sub_block(b_cl, 0, 0, 1, 0, 2, 2); @@ -112,7 +115,8 @@ TEST(MathMatrixRevCL, sub_block_triangular) { EXPECT_EQ(c(2, 1), 1); EXPECT_EQ(c(2, 2), 0); - a_cl = stan::math::to_matrix_cl(a); + a_cl.val() = stan::math::to_matrix_cl(a.val()); + a_cl.adj() = stan::math::to_matrix_cl(a.adj()); a_cl.view(stan::math::matrix_cl_view::Lower); b_cl.view(stan::math::matrix_cl_view::Upper); a_cl.sub_block(b_cl, 0, 0, 0, 0, 2, 2); @@ -128,7 +132,8 @@ TEST(MathMatrixRevCL, sub_block_triangular) { EXPECT_EQ(c(2, 1), 0); EXPECT_EQ(c(2, 2), 0); - a_cl = stan::math::to_matrix_cl(a); + a_cl.val() = stan::math::to_matrix_cl(a.val()); + a_cl.adj() = stan::math::to_matrix_cl(a.adj()); a_cl.view(stan::math::matrix_cl_view::Lower); b_cl.view(stan::math::matrix_cl_view::Upper); a_cl.sub_block(b_cl, 1, 0, 1, 0, 2, 2); @@ -144,7 +149,8 @@ TEST(MathMatrixRevCL, sub_block_triangular) { EXPECT_EQ(c(2, 1), 0); EXPECT_EQ(c(2, 2), 0); - a_cl = stan::math::to_matrix_cl(a); + a_cl.val() = stan::math::to_matrix_cl(a.val()); + a_cl.adj() = stan::math::to_matrix_cl(a.adj()); a_cl.view(stan::math::matrix_cl_view::Lower); b_cl.view(stan::math::matrix_cl_view::Upper); a_cl.sub_block(b_cl, 1, 0, 1, 0, 2, 3); @@ -160,7 +166,8 @@ TEST(MathMatrixRevCL, sub_block_triangular) { EXPECT_EQ(c(2, 1), 0); EXPECT_EQ(c(2, 2), 1); - a_cl = stan::math::to_matrix_cl(a); + a_cl.val() = stan::math::to_matrix_cl(a.val()); + a_cl.adj() = stan::math::to_matrix_cl(a.adj()); a_cl.view(stan::math::matrix_cl_view::Upper); b_cl.view(stan::math::matrix_cl_view::Upper); a_cl.sub_block(b_cl, 1, 0, 1, 0, 2, 2); @@ -176,7 +183,8 @@ TEST(MathMatrixRevCL, sub_block_triangular) { EXPECT_EQ(c(2, 1), 0); EXPECT_EQ(c(2, 2), 0); - a_cl = stan::math::to_matrix_cl(a); + a_cl.val() = stan::math::to_matrix_cl(a.val()); + a_cl.adj() = stan::math::to_matrix_cl(a.adj()); a_cl.view(stan::math::matrix_cl_view::Upper); b_cl.view(stan::math::matrix_cl_view::Upper); a_cl.sub_block(b_cl, 0, 1, 1, 1, 2, 2); @@ -192,7 +200,8 @@ TEST(MathMatrixRevCL, sub_block_triangular) { EXPECT_EQ(c(2, 1), 1); EXPECT_EQ(c(2, 2), 1); - a_cl = stan::math::to_matrix_cl(a); + a_cl.val() = stan::math::to_matrix_cl(a.val()); + a_cl.adj() = stan::math::to_matrix_cl(a.adj()); a_cl.view(stan::math::matrix_cl_view::Upper); b_cl.view(stan::math::matrix_cl_view::Lower); a_cl.sub_block(b_cl, 0, 0, 0, 1, 2, 2); @@ -208,7 +217,8 @@ TEST(MathMatrixRevCL, sub_block_triangular) { EXPECT_EQ(c(2, 1), 0); EXPECT_EQ(c(2, 2), 0); - a_cl = stan::math::to_matrix_cl(a); + a_cl.val() = stan::math::to_matrix_cl(a.val()); + a_cl.adj() = stan::math::to_matrix_cl(a.adj()); a_cl.view(stan::math::matrix_cl_view::Upper); b_cl.view(stan::math::matrix_cl_view::Lower); a_cl.sub_block(b_cl, 0, 0, 0, 0, 2, 2); @@ -224,7 +234,8 @@ TEST(MathMatrixRevCL, sub_block_triangular) { EXPECT_EQ(c(2, 1), 0); EXPECT_EQ(c(2, 2), 0); - a_cl = stan::math::to_matrix_cl(a); + a_cl.val() = stan::math::to_matrix_cl(a.val()); + a_cl.adj() = stan::math::to_matrix_cl(a.adj()); a_cl.view(stan::math::matrix_cl_view::Upper); b_cl.view(stan::math::matrix_cl_view::Lower); a_cl.sub_block(b_cl, 0, 1, 0, 1, 2, 2); @@ -240,7 +251,8 @@ TEST(MathMatrixRevCL, sub_block_triangular) { EXPECT_EQ(c(2, 1), 0); EXPECT_EQ(c(2, 2), 0); - a_cl = stan::math::to_matrix_cl(a); + a_cl.val() = stan::math::to_matrix_cl(a.val()); + a_cl.adj() = stan::math::to_matrix_cl(a.adj()); a_cl.view(stan::math::matrix_cl_view::Upper); b_cl.view(stan::math::matrix_cl_view::Lower); a_cl.sub_block(b_cl, 0, 1, 0, 1, 3, 2); @@ -256,7 +268,8 @@ TEST(MathMatrixRevCL, sub_block_triangular) { EXPECT_EQ(c(2, 1), 1); EXPECT_EQ(c(2, 2), 1); - a_cl = stan::math::to_matrix_cl(a); + a_cl.val() = stan::math::to_matrix_cl(a.val()); + a_cl.adj() = stan::math::to_matrix_cl(a.adj()); a_cl.view(stan::math::matrix_cl_view::Upper); b_cl.view(stan::math::matrix_cl_view::Lower); a_cl.sub_block(b_cl, 0, 0, 0, 0, 3, 3); diff --git a/test/unit/math/opencl/rev/vari_test.cpp b/test/unit/math/opencl/rev/vari_test.cpp new file mode 100644 index 00000000000..4960f887324 --- /dev/null +++ b/test/unit/math/opencl/rev/vari_test.cpp @@ -0,0 +1,25 @@ +#ifdef STAN_OPENCL +#include +#include +#include +#include + +TEST(AgradRev, matrix_cl_vari_block) { + using stan::math::vari_value; + Eigen::MatrixXd a = Eigen::MatrixXd::Random(3, 3); + Eigen::MatrixXd b = Eigen::MatrixXd::Random(3, 3); + stan::math::matrix_cl a_cl(a); + stan::math::matrix_cl b_cl(b); + + vari_value> A(a); + EXPECT_MATRIX_EQ(a.block(0, 1, 2, 2), + stan::math::from_matrix_cl(A.block(0, 1, 2, 2).val_)); + vari_value> B(a_cl); + B.adj_ = b_cl; + EXPECT_MATRIX_EQ(a.block(0, 1, 2, 2), + stan::math::from_matrix_cl(B.block(0, 1, 2, 2).val_)); + EXPECT_MATRIX_EQ(b.block(0, 1, 2, 2), + stan::math::from_matrix_cl(B.block(0, 1, 2, 2).adj_)); +} + +#endif diff --git a/test/unit/math/prim/meta/is_plain_type_test.cpp b/test/unit/math/prim/meta/is_plain_type_test.cpp new file mode 100644 index 00000000000..a845657ede0 --- /dev/null +++ b/test/unit/math/prim/meta/is_plain_type_test.cpp @@ -0,0 +1,28 @@ +#include +#include +#include + +TEST(MathMetaPrim, is_plain_type) { + using stan::is_plain_type; + EXPECT_TRUE((is_plain_type::value)); + EXPECT_TRUE((is_plain_type::value)); + EXPECT_TRUE((is_plain_type::value)); + + EXPECT_TRUE((is_plain_type>::value)); + EXPECT_TRUE((is_plain_type>::value)); + EXPECT_TRUE((is_plain_type&>::value)); + + EXPECT_TRUE((is_plain_type>::value)); + EXPECT_TRUE((is_plain_type>::value)); + + EXPECT_TRUE((is_plain_type>::value)); + EXPECT_TRUE((is_plain_type&>::value)); + Eigen::Matrix a; + Eigen::Matrix b; + + EXPECT_FALSE((is_plain_type::value)); + EXPECT_FALSE((is_plain_type::value)); + EXPECT_FALSE((is_plain_type< + Eigen::MatrixBase>&&>::value)); + EXPECT_FALSE((is_plain_type>::value)); +} diff --git a/test/unit/math/rev/core/var_test.cpp b/test/unit/math/rev/core/var_test.cpp index c9bb2dfb48b..9d9d42774cc 100644 --- a/test/unit/math/rev/core/var_test.cpp +++ b/test/unit/math/rev/core/var_test.cpp @@ -1,5 +1,6 @@ #include #include +#include #include #include #include @@ -14,8 +15,10 @@ struct AgradRev : public testing::Test { } }; +namespace stan { +namespace test { template -void ctor_overloads_impl() { +void ctor_overloads_float_impl() { using stan::math::var_value; using stan::math::vari_value; using stan::math::test::type_name; @@ -39,27 +42,122 @@ void ctor_overloads_impl() { } template -void ctor_overloads() { - ctor_overloads_impl(); - ctor_overloads_impl(); - ctor_overloads_impl(); - ctor_overloads_impl(); - ctor_overloads_impl(); - ctor_overloads_impl(); - ctor_overloads_impl(); - ctor_overloads_impl(); - ctor_overloads_impl(); - ctor_overloads_impl(); - ctor_overloads_impl(); - ctor_overloads_impl(); - ctor_overloads_impl(); - ctor_overloads_impl(); +void ctor_overloads_float() { + ctor_overloads_float_impl(); + ctor_overloads_float_impl(); + ctor_overloads_float_impl(); + ctor_overloads_float_impl(); + ctor_overloads_float_impl(); + ctor_overloads_float_impl(); + ctor_overloads_float_impl(); + ctor_overloads_float_impl(); + ctor_overloads_float_impl(); + ctor_overloads_float_impl(); + ctor_overloads_float_impl(); + ctor_overloads_float_impl(); + ctor_overloads_float_impl(); + ctor_overloads_float_impl(); } -TEST_F(AgradRev, ctorOverloads) { - ctor_overloads(); - ctor_overloads(); - ctor_overloads(); +template +void ctor_overloads_matrix(EigenMat&& xx) { + using stan::math::var_value; + using stan::math::vari_value; + using stan::math::test::type_name; + using eigen_plain = std::decay_t>; + + eigen_plain x = xx; + // standard constructor + EXPECT_MATRIX_FLOAT_EQ((x + x).eval(), var_value(x + x).val()); + // make sure copy ctor is used rather than casting vari* to unsigned int + EXPECT_MATRIX_FLOAT_EQ( + x, var_value(new vari_value(x)).val()); + // make sure rvalue var_value can be accepted + EXPECT_MATRIX_FLOAT_EQ( + x, var_value(var_value(x)).val()); + // test init_dependent for adj + auto test_var_x = var_value(var_value(x)); + test_var_x.vi_->init_dependent(); + EXPECT_MATRIX_FLOAT_EQ(eigen_plain::Ones(x.rows(), x.cols()), + test_var_x.adj()); +} + +template +void ctor_overloads_sparse_matrix(EigenMat&& x) { + using stan::math::var_value; + using stan::math::vari_value; + using stan::math::test::type_name; + using eigen_plain = std::decay_t>; + using inner_iterator = typename eigen_plain::InnerIterator; + // standard constructor with eigen expression + eigen_plain matmul_x = x * x; + eigen_plain matmul_xx = var_value(x * x).val(); + for (int k = 0; k < matmul_x.outerSize(); ++k) { + for (inner_iterator it(matmul_x, k), iz(matmul_xx, k); it; ++it, ++iz) { + EXPECT_FLOAT_EQ(iz.value(), it.value()); + } + } + const eigen_plain const_matmul_x = x * x; + eigen_plain const_matmul_xx = var_value(const_matmul_x).val(); + for (int k = 0; k < matmul_x.outerSize(); ++k) { + for (inner_iterator it(const_matmul_x, k), iz(const_matmul_xx, k); it; + ++it, ++iz) { + EXPECT_FLOAT_EQ(iz.value(), it.value()); + } + } + + // make sure rvalue var_value can be accepted + eigen_plain x_rv = var_value(var_value(x)).val(); + for (int k = 0; k < x.outerSize(); ++k) { + for (inner_iterator it(x, k), iz(x_rv, k); it; ++it, ++iz) { + EXPECT_FLOAT_EQ(iz.value(), it.value()); + } + } + + // from a vari_value with sparse eigen expression + eigen_plain x_from_vari + = var_value(new vari_value(x * x)).val(); + for (int k = 0; k < matmul_x.outerSize(); ++k) { + for (inner_iterator it(matmul_x, k), iz(x_from_vari, k); it; ++it, ++iz) { + EXPECT_FLOAT_EQ(iz.value(), it.value()); + } + } + // test inplace addition works + auto inplace_add_var = var_value(new vari_value(x)); + eigen_plain test_y = make_sparse_matrix_random(10, 10); + inplace_add_var.vi_->init_dependent(); + inplace_add_var.adj() += test_y; + // adjoints sparsity pattern will be pattern of x and test_y for addition + for (int k = 0; k < x.outerSize(); ++k) { + for (inner_iterator it(test_y, k), iz(inplace_add_var.adj(), k); iz; ++iz) { + if (iz.row() == it.row() && iz.col() == it.col()) { + EXPECT_FLOAT_EQ(iz.value() - 1, it.value()); + ++it; + } else { + EXPECT_FLOAT_EQ(iz.value(), 1.0); + } + } + } +} + +} // namespace test +} // namespace stan +TEST_F(AgradRev, ctorfloatOverloads) { + stan::test::ctor_overloads_float(); + stan::test::ctor_overloads_float(); + stan::test::ctor_overloads_float(); +} + +TEST_F(AgradRev, ctormatrixOverloads) { + using dense_mat = Eigen::Matrix; + using sparse_mat = Eigen::SparseMatrix; + stan::test::ctor_overloads_matrix(dense_mat::Random(10, 10)); + using dense_vec = Eigen::Matrix; + stan::test::ctor_overloads_matrix(dense_vec::Random(10)); + using dense_row_vec = Eigen::Matrix; + stan::test::ctor_overloads_matrix(dense_row_vec::Random(10)); + sparse_mat sparse_x = stan::test::make_sparse_matrix_random(10, 10); + stan::test::ctor_overloads_sparse_matrix(sparse_x); } TEST_F(AgradRev, a_eq_x) { diff --git a/test/unit/math/rev/core/vari_test.cpp b/test/unit/math/rev/core/vari_test.cpp index fc5a6e34ab5..4ad4497c010 100644 --- a/test/unit/math/rev/core/vari_test.cpp +++ b/test/unit/math/rev/core/vari_test.cpp @@ -1,4 +1,5 @@ #include +#include #include #include @@ -15,3 +16,44 @@ TEST(AgradRev, long_double_test) { ss << &v; EXPECT_EQ("5:0", ss.str()); } + +TEST(AgradRev, dense_matrix_vari) { + using stan::math::vari_value; + using eig_mat = Eigen::MatrixXd; + vari_value A_vari(eig_mat::Random(3, 3)); + eig_mat B(eig_mat::Random(3, 3)); + vari_value B_vari(B); + EXPECT_MATRIX_FLOAT_EQ(B, B_vari.val_); +} +TEST(AgradRev, dense_vector_vari) { + using stan::math::vari_value; + using eig_vec = Eigen::Matrix; + vari_value A_vari(eig_vec::Random(3)); + eig_vec B(eig_vec::Random(3)); + vari_value B_vari(B); + EXPECT_MATRIX_FLOAT_EQ(B, B_vari.val_); +} + +TEST(AgradRev, dense_row_vector_vari) { + using stan::math::vari_value; + using eig_row_vec = Eigen::Matrix; + vari_value A_vari(eig_row_vec::Random(3)); + eig_row_vec B(eig_row_vec::Random(3)); + vari_value B_vari(B); + EXPECT_MATRIX_FLOAT_EQ(B, B_vari.val_); +} + +TEST(AgradRev, sparse_matrix_vari) { + using stan::math::vari_value; + using eig_mat = Eigen::SparseMatrix; + using inner_iterator = typename eig_mat::InnerIterator; + using stan::test::make_sparse_matrix_random; + vari_value A_vari(make_sparse_matrix_random(10, 10)); + eig_mat B = make_sparse_matrix_random(10, 10); + vari_value B_vari(B); + for (int k = 0; k < B.outerSize(); ++k) { + for (inner_iterator it(B, k), iz(B_vari.val_, k); it; ++it, ++iz) { + EXPECT_FLOAT_EQ(iz.value(), it.value()); + } + } +} diff --git a/test/unit/math/rev/meta/is_fvar_test.cpp b/test/unit/math/rev/meta/is_fvar_test.cpp index 4a159531f99..46e38837ca7 100644 --- a/test/unit/math/rev/meta/is_fvar_test.cpp +++ b/test/unit/math/rev/meta/is_fvar_test.cpp @@ -1,4 +1,5 @@ #include +#include #include TEST(MetaTraitsRevScal, is_fvar) { diff --git a/test/unit/math/rev/meta/is_var_or_arithmetic_test.cpp b/test/unit/math/rev/meta/is_var_or_arithmetic_test.cpp index ff1188df1bc..4778049da7a 100644 --- a/test/unit/math/rev/meta/is_var_or_arithmetic_test.cpp +++ b/test/unit/math/rev/meta/is_var_or_arithmetic_test.cpp @@ -1,4 +1,5 @@ #include +#include #include #include #include diff --git a/test/unit/math/rev/meta/is_var_test.cpp b/test/unit/math/rev/meta/is_var_test.cpp index cb225ff7d1c..ec276e7ab7a 100644 --- a/test/unit/math/rev/meta/is_var_test.cpp +++ b/test/unit/math/rev/meta/is_var_test.cpp @@ -1,4 +1,5 @@ #include +#include #include TEST(MetaTraitsRevScal, is_var) { diff --git a/test/unit/util.hpp b/test/unit/util.hpp index 4f0bf854254..b8746686e56 100644 --- a/test/unit/util.hpp +++ b/test/unit/util.hpp @@ -2,6 +2,8 @@ #define TEST_UNIT_UTIL_HPP #include +#include +#include #include #include #include @@ -195,6 +197,25 @@ void expect_type_matrix(const Eigen::Matrix& x) { EXPECT_EQ(Eigen::Dynamic, C); EXPECT_EQ(Eigen::Dynamic, R); } + +auto make_sparse_matrix_random(int rows, int cols) { + using eigen_triplet = Eigen::Triplet; + boost::mt19937 gen; + boost::random::uniform_real_distribution dist(0.0, 1.0); + std::vector tripletList; + for (int i = 0; i < rows; ++i) { + for (int j = 0; j < cols; ++j) { + auto v_ij = dist(gen); + if (v_ij < 0.1) { + tripletList.push_back(eigen_triplet(i, j, v_ij)); + } + } + } + Eigen::SparseMatrix mat(rows, cols); + mat.setFromTriplets(tripletList.begin(), tripletList.end()); + return mat; +} + } // namespace test } // namespace stan