Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Generalize operator-like functions #1628

Merged
merged 34 commits into from
Feb 14, 2020
Merged

Generalize operator-like functions #1628

Show file tree
Hide file tree
Changes from 4 commits
Commits
Show all changes
34 commits
Select commit Hold shift + click to select a range
84b5ab6
generalized operator-like functions
t4c1 Jan 20, 2020
7d7779e
optimize fwd multiply
t4c1 Jan 20, 2020
4b50381
Merge commit '10cc6ba675743f09832c6749fbb1a92d74888bd2' into HEAD
yashikno Jan 20, 2020
3c23987
[Jenkins] auto-formatting by clang-format version 6.0.0 (tags/google/…
stan-buildbot Jan 20, 2020
7cd74ef
fixed call sites.
t4c1 Jan 20, 2020
8263fdb
fixed a test.
t4c1 Jan 20, 2020
fb2de8a
generalized transpose.hpp() input
t4c1 Jan 21, 2020
1fedd01
fixed transpose
t4c1 Jan 21, 2020
9029b14
[Jenkins] auto-formatting by clang-format version 5.0.2-svn328729-1~e…
stan-buildbot Jan 21, 2020
d73d0c8
removed compile time multiplicable requirements
t4c1 Jan 21, 2020
45319be
Merge commit '12cbb63adfcfb64440ade29e97a25db30c5ac8e2' into HEAD
yashikno Jan 21, 2020
f989b11
fix unnecessary vector to matrix conversions
t4c1 Jan 22, 2020
01a87f4
Merge commit '81248fdcd9bbd107f63cf8211333c2acbdbe5ed8' into HEAD
yashikno Jan 22, 2020
d34d6fe
[Jenkins] auto-formatting by clang-format version 6.0.0 (tags/google/…
stan-buildbot Jan 22, 2020
c68504a
addressed review comments
t4c1 Jan 22, 2020
aaa8c17
fix cpplint
t4c1 Jan 22, 2020
b98a5a0
[Jenkins] auto-formatting by clang-format version 6.0.0 (tags/google/…
stan-buildbot Jan 22, 2020
9445e07
Merge branch 'generalize_operators' of https://github.com/bstatcomp/m…
t4c1 Jan 22, 2020
9180fc6
generalized matrix_exp
t4c1 Jan 22, 2020
ed9b73a
Merge branch 'generalize_operators' of https://github.com/bstatcomp/m…
t4c1 Jan 22, 2020
f2b4e13
[Jenkins] auto-formatting by clang-format version 6.0.0 (tags/google/…
stan-buildbot Jan 22, 2020
d989040
replaced expression return values with matrices.
t4c1 Jan 27, 2020
ca27b44
fix dot_product
t4c1 Jan 28, 2020
e7d4d53
addressed review comments and generalized dot_product
t4c1 Feb 3, 2020
87d205c
Merge commit 'b3dde68ece20db22da4979175c07a4d42ccdc050' into HEAD
yashikno Feb 3, 2020
a9490b3
[Jenkins] auto-formatting by clang-format version 6.0.0 (tags/google/…
stan-buildbot Feb 3, 2020
c7162fd
added missing include
t4c1 Feb 3, 2020
b160720
one more header
t4c1 Feb 3, 2020
f88d89a
Merge branch 'develop' into generalize_operators
t4c1 Feb 10, 2020
d9284e8
Merge branch 'develop' into generalize_operators
t4c1 Feb 11, 2020
cf2b72d
fixed size zero multiply
t4c1 Feb 12, 2020
fb619c0
Merge branch 'develop' into generalize_operators
t4c1 Feb 12, 2020
5290f6b
fixed multiply templates
t4c1 Feb 12, 2020
9529ccc
addressed remaining review comments
t4c1 Feb 12, 2020
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
42 changes: 0 additions & 42 deletions stan/math/fwd/fun/dot_product.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -9,48 +9,6 @@
namespace stan {
namespace math {

template <typename T, int R1, int C1, int R2, int C2>
inline fvar<T> dot_product(const Eigen::Matrix<fvar<T>, R1, C1>& v1,
const Eigen::Matrix<fvar<T>, R2, C2>& v2) {
check_vector("dot_product", "v1", v1);
check_vector("dot_product", "v2", v2);
check_matching_sizes("dot_product", "v1", v1, "v2", v2);

fvar<T> ret(0, 0);
for (size_type i = 0; i < v1.size(); i++) {
ret += v1(i) * v2(i);
}
return ret;
}

template <typename T, int R1, int C1, int R2, int C2>
inline fvar<T> dot_product(const Eigen::Matrix<fvar<T>, R1, C1>& v1,
const Eigen::Matrix<double, R2, C2>& v2) {
check_vector("dot_product", "v1", v1);
check_vector("dot_product", "v2", v2);
check_matching_sizes("dot_product", "v1", v1, "v2", v2);

fvar<T> ret(0, 0);
for (size_type i = 0; i < v1.size(); i++) {
ret += v1(i) * v2(i);
}
return ret;
}

template <typename T, int R1, int C1, int R2, int C2>
inline fvar<T> dot_product(const Eigen::Matrix<double, R1, C1>& v1,
const Eigen::Matrix<fvar<T>, R2, C2>& v2) {
check_vector("dot_product", "v1", v1);
check_vector("dot_product", "v2", v2);
check_matching_sizes("dot_product", "v1", v1, "v2", v2);

fvar<T> ret(0, 0);
for (size_type i = 0; i < v1.size(); i++) {
ret += v1(i) * v2(i);
}
return ret;
}

template <typename T, int R1, int C1, int R2, int C2>
inline fvar<T> dot_product(const Eigen::Matrix<fvar<T>, R1, C1>& v1,
const Eigen::Matrix<fvar<T>, R2, C2>& v2,
Expand Down
131 changes: 38 additions & 93 deletions stan/math/fwd/fun/multiply.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -1,6 +1,7 @@
#ifndef STAN_MATH_FWD_FUN_MULTIPLY_HPP
#define STAN_MATH_FWD_FUN_MULTIPLY_HPP

#include <stan/math/prim/meta.hpp>
#include <stan/math/prim/err.hpp>
#include <stan/math/prim/fun/Eigen.hpp>
#include <stan/math/fwd/core.hpp>
Expand All @@ -10,120 +11,64 @@
namespace stan {
namespace math {

template <typename T, int R1, int C1>
inline Eigen::Matrix<fvar<T>, R1, C1> multiply(
const Eigen::Matrix<fvar<T>, R1, C1>& m, const fvar<T>& c) {
Eigen::Matrix<fvar<T>, R1, C1> res(m.rows(), m.cols());
for (int i = 0; i < m.size(); i++)
res(i) = c * m(i);
return res;
}

template <typename T, int R2, int C2>
inline Eigen::Matrix<fvar<T>, R2, C2> multiply(
const Eigen::Matrix<fvar<T>, R2, C2>& m, double c) {
Eigen::Matrix<fvar<T>, R2, C2> res(m.rows(), m.cols());
for (int i = 0; i < m.size(); i++)
res(i) = c * m(i);
return res;
}

template <typename T, int R1, int C1>
inline Eigen::Matrix<fvar<T>, R1, C1> multiply(
const Eigen::Matrix<double, R1, C1>& m, const fvar<T>& c) {
Eigen::Matrix<fvar<T>, R1, C1> res(m.rows(), m.cols());
for (int i = 0; i < m.size(); i++)
res(i) = c * m(i);
return res;
}

template <typename T, int R1, int C1>
inline Eigen::Matrix<fvar<T>, R1, C1> multiply(
const fvar<T>& c, const Eigen::Matrix<fvar<T>, R1, C1>& m) {
return multiply(m, c);
}

template <typename T, int R1, int C1>
inline Eigen::Matrix<fvar<T>, R1, C1> multiply(
double c, const Eigen::Matrix<fvar<T>, R1, C1>& m) {
return multiply(m, c);
}

template <typename T, int R1, int C1>
inline Eigen::Matrix<fvar<T>, R1, C1> multiply(
const fvar<T>& c, const Eigen::Matrix<double, R1, C1>& m) {
return multiply(m, c);
}

template <typename T, int R1, int C1, int R2, int C2>
inline Eigen::Matrix<fvar<T>, R1, C2> multiply(
const Eigen::Matrix<fvar<T>, R1, C1>& m1,
const Eigen::Matrix<fvar<T>, R2, C2>& m2) {
template <
typename T1, typename T2, typename = require_all_eigen_vt<is_fvar, T1, T2>,
typename = require_same_t<typename value_type_t<T1>::Scalar,
typename value_type_t<T2>::Scalar>,
typename
= require_t<bool_constant<T1::ColsAtCompileTime == T2::RowsAtCompileTime>>,
typename = require_not_t<
conjunction<is_eigen_row_vector<T1>, is_eigen_col_vector<T2>>>,
int = 0>
Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This int template parameter is required to make compiler treat these as different overloads (if they have same arguments and template parameters it treats them as same overload and gives a lot of errors about redefining default template parameter values.

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

! O that's how you do that. I was fighting with that compile error like a week ago. Thank you. :)

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

To me this solution seems a bit hacky, but I don't know any better one.

Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I'm not a fan of this approach. It's feels like it fights against the language. I think it would be good to post a stackoverflow Q about this

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I found a nice solution here: https://stackoverflow.com/questions/38305222/default-template-argument-when-using-stdenable-if-as-templ-param-why-ok-wit

The trick is to use requires as type of template parameters instead of defaults to typename parameters - require_t<...>* = nullptr instead of typename = require_t<...>

inline auto multiply(const T1& m1, const T2& m2) {
check_multiplicable("multiply", "m1", m1, "m2", m2);
Eigen::Matrix<fvar<T>, R1, C2> result(m1.rows(), m2.cols());
for (size_type i = 0; i < m1.rows(); i++) {
Eigen::Matrix<fvar<T>, 1, C1> crow = m1.row(i);
for (size_type j = 0; j < m2.cols(); j++) {
Eigen::Matrix<fvar<T>, R2, 1> ccol = m2.col(j);
result(i, j) = dot_product(crow, ccol);
}
}
return result;
return m1 * m2;
SteveBronder marked this conversation as resolved.
Show resolved Hide resolved
}

template <typename T, int R1, int C1, int R2, int C2>
inline Eigen::Matrix<fvar<T>, R1, C2> multiply(
const Eigen::Matrix<fvar<T>, R1, C1>& m1,
const Eigen::Matrix<double, R2, C2>& m2) {
template <typename T1, typename T2, typename = require_all_eigen_t<T1, T2>,
typename = require_fvar_t<value_type_t<T1>>,
typename = require_same_vt<double, T2>,
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Could shorten this up with

template <typename Mat1, typename Mat2, typename = require_eigen_vt<is_fvar, Mat1>,
  typename = require_eigen_vt<std::is_floating_point, Mat2>

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

done

typename = require_t<
bool_constant<T1::ColsAtCompileTime == T2::RowsAtCompileTime>>,
typename = require_not_t<
conjunction<is_eigen_row_vector<T1>, is_eigen_col_vector<T2>>>,
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Maybe just a meta function for is_dot_product for these instead of a requires specialization

int = 0>
inline auto multiply(const T1& m1, const T2& m2) {
check_multiplicable("multiply", "m1", m1, "m2", m2);
Eigen::Matrix<fvar<T>, R1, C2> result(m1.rows(), m2.cols());
Eigen::Matrix<value_type_t<T1>, T1::RowsAtCompileTime, T2::ColsAtCompileTime>
result(m1.rows(), m2.cols());
for (size_type i = 0; i < m1.rows(); i++) {
Eigen::Matrix<fvar<T>, 1, C1> crow = m1.row(i);
Eigen::Matrix<value_type_t<T1>, 1, T1::ColsAtCompileTime> crow = m1.row(i);
SteveBronder marked this conversation as resolved.
Show resolved Hide resolved
for (size_type j = 0; j < m2.cols(); j++) {
Eigen::Matrix<double, R2, 1> ccol = m2.col(j);
auto ccol = m2.col(j);
result(i, j) = dot_product(crow, ccol);
}
}
return result;
}

template <typename T, int R1, int C1, int R2, int C2>
inline Eigen::Matrix<fvar<T>, R1, C2> multiply(
const Eigen::Matrix<double, R1, C1>& m1,
const Eigen::Matrix<fvar<T>, R2, C2>& m2) {
template <typename T1, typename T2, typename = require_all_eigen_t<T1, T2>,
typename = require_same_vt<double, T1>,
typename = require_fvar_t<value_type_t<T2>>,
typename = require_t<
bool_constant<T1::ColsAtCompileTime == T2::RowsAtCompileTime>>,
typename = require_not_t<
conjunction<is_eigen_row_vector<T1>, is_eigen_col_vector<T2>>>,
char = 0>
inline auto multiply(const T1& m1, const T2& m2) {
check_multiplicable("multiply", "m1", m1, "m2", m2);
Eigen::Matrix<fvar<T>, R1, C2> result(m1.rows(), m2.cols());
Eigen::Matrix<value_type_t<T2>, T1::RowsAtCompileTime, T2::ColsAtCompileTime>
result(m1.rows(), m2.cols());
for (size_type i = 0; i < m1.rows(); i++) {
Eigen::Matrix<double, 1, C1> crow = m1.row(i);
Eigen::Matrix<double, 1, T1::ColsAtCompileTime> crow = m1.row(i);
for (size_type j = 0; j < m2.cols(); j++) {
Eigen::Matrix<fvar<T>, R2, 1> ccol = m2.col(j);
auto ccol = m2.col(j);
result(i, j) = dot_product(crow, ccol);
}
}
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Is there a reason to not put this as

for (size_type i = 0; i < m1.rows(); i++) {
  for (size_type j = 0; j < m2.cols(); j++) {
    result(i, j) = dot_product(m1.row(i), m2.col(j));
  }
}

Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I'm not a great assembly guy, but my usual heuristic of 'less is better' says the above looks a bit nicer

https://godbolt.org/z/hRS-jV

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This is a case where less is not better.

Doing Eigen::Matrix<double, 1, T1::ColsAtCompileTime> crow = m1.row(i); in outer loop benefits performance for large matrices. Consecutive elements in same row of the matrix are not in consecutive memory locations, so accessing them in order uses cache poorly. Copying them to a vector makes them consecutive in memory and accessing them faster.

Since we go over them in every iteration of inner loop it is better to inefficiently access them just once and make remaining accesses efficient.

If you want I can make line 51 na 52 into a single line, but it won't affect performance in any way.

Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Doing Eigen::Matrix<double, 1, T1::ColsAtCompileTime> crow = m1.row(i); in outer loop benefits performance for large matrices. Consecutive elements in same row of the matrix are not in consecutive memory locations, so accessing them in order uses cache poorly. Copying them to a vector makes them consecutive in memory and accessing them faster.

Since we go over them in every iteration of inner loop it is better to inefficiently access them just once and make remaining accesses efficient.

Word that's sensible to me.

If you want I can make line 51 na 52 into a single line, but it won't affect performance in any way.

idt it will effect performance though I think zipping 51 into 52 would be nice

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

done

return result;
}

template <typename T, int C1, int R2>
inline fvar<T> multiply(const Eigen::Matrix<fvar<T>, 1, C1>& rv,
const Eigen::Matrix<fvar<T>, R2, 1>& v) {
check_multiplicable("multiply", "rv", rv, "v", v);
return dot_product(rv, v);
}

template <typename T, int C1, int R2>
inline fvar<T> multiply(const Eigen::Matrix<fvar<T>, 1, C1>& rv,
const Eigen::Matrix<double, R2, 1>& v) {
check_multiplicable("multiply", "rv", rv, "v", v);
return dot_product(rv, v);
}

template <typename T, int C1, int R2>
inline fvar<T> multiply(const Eigen::Matrix<double, 1, C1>& rv,
const Eigen::Matrix<fvar<T>, R2, 1>& v) {
check_multiplicable("multiply", "rv", rv, "v", v);
return dot_product(rv, v);
}

} // namespace math
} // namespace stan
#endif
32 changes: 12 additions & 20 deletions stan/math/prim/err/check_matching_dims.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -14,12 +14,8 @@ namespace math {
/**
* Check if the two matrices are of the same size.
* This function checks the runtime sizes only.
* @tparam T1 scalar type of the first matrix
* @tparam T2 scalar type of the second matrix
* @tparam R1 number of rows in the first matrix, can be Eigen::Dynamic
* @tparam C1 number of columns in the first matrix, can be Eigen::Dynamic
* @tparam R2 number of rows in the second matrix, can be Eigen::Dynamic
* @tparam C2 number of columns in the second matrix, can be Eigen::Dynamic
* @tparam T1 type of the first matrix
* @tparam T2 type of the second matrix
* @param function name of function (for error messages)
* @param name1 variable name for the first matrix (for error messages)
* @param y1 first matrix to test
Expand All @@ -28,11 +24,9 @@ namespace math {
* @throw <code>std::invalid_argument</code> if the dimensions of the
* matrices do not match
*/
template <typename T1, typename T2, int R1, int C1, int R2, int C2>
template <typename T1, typename T2, typename = require_all_eigen_t<T1, T2>>
inline void check_matching_dims(const char* function, const char* name1,
const Eigen::Matrix<T1, R1, C1>& y1,
const char* name2,
const Eigen::Matrix<T2, R2, C2>& y2) {
const T1& y1, const char* name2, const T2& y2) {
check_size_match(function, "Rows of ", name1, y1.rows(), "rows of ", name2,
y2.rows());
check_size_match(function, "Columns of ", name1, y1.cols(), "columns of ",
Expand All @@ -47,10 +41,6 @@ inline void check_matching_dims(const char* function, const char* name1,
* @tparam check_compile Whether to check the static sizes
* @tparam T1 scalar type of the first matrix
* @tparam T2 scalar type of the second matrix
* @tparam R1 number of rows in the first matrix, can be Eigen::Dynamic
* @tparam C1 number of columns in the first matrix, can be Eigen::Dynamic
* @tparam R2 number of rows in the second matrix, can be Eigen::Dynamic
* @tparam C2 number of columns in the second matrix, can be Eigen::Dynamic
* @param function name of function (for error messages)
* @param name1 variable name for the first matrix (for error messages)
* @param y1 first matrix to test
Expand All @@ -59,13 +49,15 @@ inline void check_matching_dims(const char* function, const char* name1,
* @throw <code>std::invalid_argument</code> if the dimensions of the matrices
* do not match
*/
template <bool check_compile, typename T1, typename T2, int R1, int C1, int R2,
int C2>
template <bool check_compile, typename T1, typename T2,
typename = require_all_eigen_t<T1, T2>>
inline void check_matching_dims(const char* function, const char* name1,
const Eigen::Matrix<T1, R1, C1>& y1,
const char* name2,
const Eigen::Matrix<T2, R2, C2>& y2) {
if (check_compile && (R1 != R2 || C1 != C2)) {
const T1& y1, const char* name2, const T2& y2) {
if (check_compile
&& (static_cast<int>(T1::RowsAtCompileTime)
!= static_cast<int>(T2::RowsAtCompileTime)
|| static_cast<int>(T1::ColsAtCompileTime)
!= static_cast<int>(T2::ColsAtCompileTime))) {
std::ostringstream msg;
msg << "Static rows and cols of " << name1 << " and " << name2
<< " must match in size.";
Expand Down
46 changes: 20 additions & 26 deletions stan/math/prim/fun/add.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -1,6 +1,7 @@
#ifndef STAN_MATH_PRIM_FUN_ADD_HPP
#define STAN_MATH_PRIM_FUN_ADD_HPP

#include <stan/math/prim/meta.hpp>
#include <stan/math/prim/err.hpp>
#include <stan/math/prim/fun/Eigen.hpp>

Expand All @@ -11,56 +12,49 @@ namespace math {
* Return the sum of the specified matrices. The two matrices
* must have the same dimensions.
*
* @tparam T1 type of elements in the first matrix
* @tparam T2 type of elements in the second matrix
* @tparam R number of rows, can be Eigen::Dynamic
* @tparam C number of columns, can be Eigen::Dynamic
* @tparam T1 type of the first matrix or expression
* @tparam T2 type of the second matrix or expression
*
* @param m1 First matrix.
* @param m2 Second matrix.
* @param m1 First matrix or expression.
* @param m2 Second matrix or expression.
* @return Sum of the matrices.
* @throw std::invalid_argument if m1 and m2 do not have the same
* dimensions.
*/
template <typename T1, typename T2, int R, int C>
inline Eigen::Matrix<return_type_t<T1, T2>, R, C> add(
const Eigen::Matrix<T1, R, C>& m1, const Eigen::Matrix<T2, R, C>& m2) {
template <typename T1, typename T2, typename = require_all_eigen_t<T1, T2>>
inline auto add(const T1& m1, const T2& m2) {
check_matching_dims("add", "m1", m1, "m2", m2);
return m1 + m2;
}

/**
* Return the sum of the specified matrix and specified scalar.
*
* @tparam T1 type of elements in the matrix
* @tparam T2 type of scalar
* @tparam R number of rows, can be Eigen::Dynamic
* @tparam C number of columns, can be Eigen::Dynamic
* @param m Matrix.
* @tparam T1 type of the matrix or expression
* @tparam T2 type of the scalar
* @param m Matrix or expression.
* @param c Scalar.
* @return The matrix plus the scalar.
*/
template <typename T1, typename T2, int R, int C>
inline Eigen::Matrix<return_type_t<T1, T2>, R, C> add(
const Eigen::Matrix<T1, R, C>& m, const T2& c) {
return m.array() + c;
template <typename T1, typename T2, typename = require_eigen_t<T1>,
typename = require_stan_scalar_t<T2>>
inline auto add(const T1& m, const T2 c) {
return (m.array() + c).matrix();
}

/**
* Return the sum of the specified scalar and specified matrix.
*
* @tparam T1 type of scalar
* @tparam T2 type of elements in the matrix
* @tparam R number of rows, can be Eigen::Dynamic
* @tparam C number of columns, can be Eigen::Dynamic
* @tparam T1 type of the scalar
* @tparam T2 type of the matrix or expression
* @param c Scalar.
* @param m Matrix.
* @return The scalar plus the matrix.
*/
template <typename T1, typename T2, int R, int C>
inline Eigen::Matrix<return_type_t<T1, T2>, R, C> add(
const T1& c, const Eigen::Matrix<T2, R, C>& m) {
return c + m.array();
template <typename T1, typename T2, typename = require_stan_scalar_t<T1>,
typename = require_eigen_t<T2>>
inline auto add(const T1 c, const T2& m) {
return (c + m.array()).matrix();
}

} // namespace math
Expand Down
Loading