-
-
Notifications
You must be signed in to change notification settings - Fork 189
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.
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
Add simplex row and column matrix constraints #2992
Conversation
@andrjohns could you take a look at this? |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Functionality makes sense. Putting in comments to help make the code cleaner going into the code base.
*/ | ||
template <typename Mat, require_eigen_matrix_dynamic_t<Mat>* = nullptr, | ||
require_not_st_var<Mat>* = nullptr> | ||
inline auto simplex_column_constrain(const Mat& y) { |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Can we make the return type more explicit? (Instead of auto
)
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Same below.
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Yeah we can make it plain_type_t<Mat>
and that should be fine
*/ | ||
template <typename Mat, require_eigen_matrix_dynamic_t<Mat>* = nullptr, | ||
require_not_st_var<Mat>* = nullptr> | ||
inline auto simplex_row_constrain(const Mat& y) { |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Same here for being more explicit with types instead of auto
.
Jenkins Console Log Machine informationNo LSB modules are available. Distributor ID: Ubuntu Description: Ubuntu 20.04.3 LTS Release: 20.04 Codename: focalCPU: G++: Clang: |
@syclik could you give this anotherl look? |
@SteveBronder, sorry, I just saw a comment that I missed. Thought I'd put it here so you can track it. (It's more of a question about the PR than a single line of code.)
The question I was asking was whether both stan/math/prim/fun/simplex_row_constrain.hpp
stan/math/rev/fun/simplex_column_constrain.hpp
The version with the jacobian also repeats. I guess I was wondering if it's better to have a C-like function that takes the pointer to data and length and consolidates the computation. We'd still have the same functions from the outside to handle types and everything. It's a question, not even a suggestion. I was just reading the code and found the logic to be complicated enough and duplicated where it made me wonder. (Repeating code twice isn't that big of a deal in the long run.) |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Thanks for updating the PR with the doc and return types.
I had a couple minor comments. I think one of the function docs had different indexes than everything else; it'd be nice to change that. Almost everything else is just using slightly more explicit return types. (It just makes it clear that the input type is the output type, which cuts out scanning the function for the return type.)
* dimensionality (K, M). | ||
*/ | ||
template <bool Jacobian, typename T, require_std_vector_t<T>* = nullptr> | ||
inline auto simplex_column_constrain(const T& y, return_type_t<T>& lp) { |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
inline auto simplex_column_constrain(const T& y, return_type_t<T>& lp) { | |
inline plain_type_t<T> simplex_column_constrain(const T& y, return_type_t<T>& lp) { |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
This one we want to keep auto
because it can return a vector with any amount of vectors of matrices within it.
@@ -24,7 +24,7 @@ namespace math { | |||
* @param y Free vector input of dimensionality K - 1. | |||
* @return Simplex of dimensionality K. | |||
*/ | |||
template <typename Vec, require_eigen_col_vector_t<Vec>* = nullptr, | |||
template <typename Vec, require_eigen_vector_t<Vec>* = nullptr, | |||
require_not_st_var<Vec>* = nullptr> | |||
inline auto simplex_constrain(const Vec& y) { |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
inline auto simplex_constrain(const Vec& y) { | |
inline plain_type_t<Vec> simplex_constrain(const Vec& y) { |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
I changed all of these to plain_type_t
@@ -56,7 +56,7 @@ inline auto simplex_constrain(const Vec& y) { | |||
* @param lp Log probability reference to increment. | |||
* @return Simplex of dimensionality K. | |||
*/ | |||
template <typename Vec, require_eigen_col_vector_t<Vec>* = nullptr, | |||
template <typename Vec, require_eigen_vector_t<Vec>* = nullptr, | |||
require_not_st_var<Vec>* = nullptr> | |||
inline auto simplex_constrain(const Vec& y, value_type_t<Vec>& lp) { |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
inline auto simplex_constrain(const Vec& y, value_type_t<Vec>& lp) { | |
inline plain_type_t<Vec> simplex_constrain(const Vec& y, value_type_t<Vec>& lp) { |
* @return vector of matrices with simplex rows of dimensionality (N, K) | ||
*/ | ||
template <bool Jacobian, typename T, require_std_vector_t<T>* = nullptr> | ||
inline auto simplex_row_constrain(const T& y, return_type_t<T>& lp) { |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
inline auto simplex_row_constrain(const T& y, return_type_t<T>& lp) { | |
inline plain_type_t<T> simplex_row_constrain(const T& y, return_type_t<T>& lp) { |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Same as above for std vectors of std vectors of matrices. idt plain_type_t here gives any more information
namespace math { | ||
|
||
namespace internal { | ||
template <typename Mat, int NewOptions, typename = void> |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Would this code be useful in stan/math/fun/Eigen.hpp
? If we think we could use it again, I'd suggest moving it there now.
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
This requires some things in meta so I'm going to move it here
* @return Matrix of simplex columns of dimensionality (N + 1, M). | ||
*/ | ||
template <typename T, require_rev_matrix_t<T>* = nullptr> | ||
auto simplex_column_constrain(const T& y, scalar_type_t<T>& lp) { |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
auto simplex_column_constrain(const T& y, scalar_type_t<T>& lp) { | |
inline plain_type_t<T> simplex_column_constrain(const T& y, scalar_type_t<T>& lp) { |
* @return Matrix with Simplexes along the rows of dimensionality (N, K) | ||
*/ | ||
template <typename T, require_rev_matrix_t<T>* = nullptr> | ||
inline auto simplex_row_constrain(const T& y) { |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
inline auto simplex_row_constrain(const T& y) { | |
inline plain_type_t<T> simplex_row_constrain(const T& y) { |
* stick-breaking process. | ||
* | ||
* @tparam T type of the matrix to constrain | ||
* @param y Free matrix input of dimensionality (N, M). |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
All the other doc says (N, K-1)
. Can we change this to be consistent?
* @tparam T type of the matrix to constrain | ||
* @param y Free matrix input of dimensionality (N, M). | ||
* @param lp Log probability reference to increment. | ||
* @return Matrix with simplexes along the rows of dimensionality (N, M + 1). |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
And also change this to (N, K)
.
* @return Matrix with simplexes along the rows of dimensionality (N, M + 1). | ||
*/ | ||
template <typename T, require_rev_matrix_t<T>* = nullptr> | ||
inline auto simplex_row_constrain(const T& y, scalar_type_t<T>& lp) { |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
inline auto simplex_row_constrain(const T& y, scalar_type_t<T>& lp) { | |
inline plain_type_t<T> simplex_row_constrain(const T& y, scalar_type_t<T>& lp) { |
@SteveBronder, minor changes requested. Once addressed, we can merge! |
Thanks @syclik! Just updated everything |
Actually @syclik do you mind if I just add the free functions to this PR? I think I'm just going to have them call the current |
@SteveBronder, go for it! I think it makes sense. |
Jenkins Console Log Machine informationNo LSB modules are available. Distributor ID: Ubuntu Description: Ubuntu 20.04.3 LTS Release: 20.04 Codename: focalCPU: G++: Clang: |
Jenkins Console Log Machine informationNo LSB modules are available. Distributor ID: Ubuntu Description: Ubuntu 20.04.3 LTS Release: 20.04 Codename: focalCPU: G++: Clang: |
@syclik I think this is ready to merge! |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
All requests addressed, approving and merging
Summary
This adds the functions
simplex_row_constrain
andsimplex_column_constrain
to be used by the stan language for matrices with simplexes along the rows or columns of the matrix. This is just the constrain version which is already a lot of code. I'm going to add the free functions in a separate PRTests
Tests were added in mix and prim and can be run with
A prim test for
simplex_column_constrain
was not written because under the hoodsimplex_column_constrain
for doubles just callssimplex_constrain
Side Effects
There is one interesting thing to point out in the reverse mode versio nof
simplex_column_constrain
. I think this is one of the only places in the Stan math library we useEigen::RowMajor
to have row major underlying matrices.For each of these algorithms we want to access memory in a contiguous manner. For rows of simplex matrices we start with a (N, M) and get back a (N, M + 1) where we operate on each column of the matrix from left to right. For columns of simplex matrices we start with an (N, M) and get back a (N + 1, M) where we need to operate on the row of the matrix from top to bottom. So for the row simplex everything is nice with column major order because that's how our access happens. But for columns of simplex vectors for reverse mode, we make a hard copy of the input changing it from column major to row major so that the rest of our memory access is contiguous. That means that for the rest of the reverse mode code we access things nicely, but then pay the cost of different storage orders when copy the input to the arena matrix in forward mode and then when we copy the adjoint from the reverse mode into the input's adjoint. But reading from different storage orders twice I think is better than reading from the noncontiguous memory for the entire algorithm.
Release notes
Adds constrain function for creating matrices holding row or column simplexes
Checklist
Copyright holder: Simons foundation
The copyright holder is typically you or your assignee, such as a university or company. By submitting this pull request, the copyright holder is agreeing to the license the submitted work under the following licenses:
- Code: BSD 3-clause (https://opensource.org/licenses/BSD-3-Clause)
- Documentation: CC-BY 4.0 (https://creativecommons.org/licenses/by/4.0/)
the basic tests are passing
./runTests.py test/unit
)make test-headers
)make test-math-dependencies
)make doxygen
)make cpplint
)the code is written in idiomatic C++ and changes are documented in the doxygen
the new changes are tested