-
-
Notifications
You must be signed in to change notification settings - Fork 187
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
Create Jacobian template parameter overload for constrain functions #2559
Conversation
…e supplied directly to the constrain function
…4.1 (tags/RELEASE_600/final)
Jenkins Console Log Machine informationProductName: Mac OS X ProductVersion: 10.11.6 BuildVersion: 15G22010CPU: G++: Clang: |
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 think you may fire me as reviewer after all these change requests. Most of them are around doc and I'm happy to help with that if you'd like---just let me know.
*/ | ||
template <bool Jacobian, typename T> | ||
inline auto cholesky_corr_constrain(const T& y, int K, scalar_type_t<T>& lp) { | ||
if (Jacobian) { |
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.
[question]
I'm unclear on how this function gets resolved since it seems to call itself inside the if (Jacobian)
block. Is this only resolved because the nested call doesn't provide a template parameter?
[comment]
I opened an issue to refactor the implementation: #2561
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.
Yes since there's no default template value for Jacobian
it diverts to the signatures without the Jacobian
template parameter
Eigen::Matrix<value_type_t<T>, Eigen::Dynamic, Eigen::Dynamic> | ||
cholesky_factor_constrain(const T& x, int M, int N, value_type_t<T>& lp) { | ||
inline Eigen::Matrix<value_type_t<T>, Eigen::Dynamic, Eigen::Dynamic> | ||
cholesky_factor_constrain(const T& x, int M, int N, scalar_type_t<T>& lp) { | ||
check_size_match("cholesky_factor_constrain", "x.size()", x.size(), |
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.
[optional]
Ideally the error messages will be phrased in a way that a Stan user can understand them from the Stan program. So maybe "constrained size" would be better than "x.size()". I don't think this can fail for Stan programs because we manage the sizes, and it's not part of this issue, so I'm marking "optional".
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.
Yep just fixed these up to `constrained size'.
template <bool Jacobian, typename T> | ||
inline auto cholesky_factor_constrain(const T& x, int M, int N, | ||
scalar_type_t<T>& lp) { | ||
if (Jacobian) { |
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.
[optional]
I generally prefer a ternary operator when appropriate, so this would be
return Jacobian ? cholesky_factor_constrain(x, M, N, lp) : cholesky_factor_constrain(X, M, N);
if it fits or
return Jacobian
? cholesky_factor_constrain(x, M, N, lp)
: cholesky_factor_constrain(X, M, N);
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.
imo I like the if (Jacobian)
for these just because it future proofs us a bit so that if we ever move to C++17 we can pretty much just do a copy replace for if (Jacobian)
for if constexpr (Jacobian)
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.
Interesting. Does the ternary operator not get reduced at compile time if the condition is constexpr
? Obviously the template parameters are all constexpr
!
* | ||
* <p>\f$\log | \frac{d}{dx} \tanh x | = \log (1 - \tanh^2 x)\f$. | ||
* | ||
* @tparam Jacobian If true, incremented `lp` with the log Jacobian |
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.
[optional]
This should just be "increment" and it's "log absolute Jacobian determinant" if you want to be precise. Or in this case, "log absolute derivative" or if you want to be verbose "log absolute derivative of the constraining transform".
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 went with the bit of doc you wrote for these in the comment below
@@ -49,6 +49,31 @@ inline auto corr_constrain(const T_x& x, T_lp& lp) { | |||
return tanh_x; | |||
} | |||
|
|||
/** | |||
* Return the result of transforming the specified scalar or container of values |
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.
The doc for all of these needs to mention that the Jacobian is incremented with the log absolute Jacobian determinant of the transform if the template parameter Jacobian
is set to true
. The goal is that the first sentence of the doc should give as thorough and precise an explanation of the function's arguments and returns as possible. If you don't think it fits in one sentence, add a second sentence on the Jacobian.
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.
Yes I appended the doc you mention below with all of these so should be good 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.
Looks good. Thanks.
|
||
/** | ||
* Return the unit length vector corresponding to the free vector y. | ||
* See https://en.wikipedia.org/wiki/N-sphere#Generating_random_points |
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 comments as other doc on return, etc.
- refer to our reference manual, which in turn has references into the literature
typename stan::scalar_type<T>::type lp = 0; | ||
stan::math::cholesky_corr_constrain(x, inv_size(x), lp); | ||
auto g3(const T& x) { | ||
stan::scalar_type_t<T> lp = 0; |
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.
[comment]
thanks for tidying up!
@@ -3,7 +3,8 @@ | |||
TEST(mathMixMatFun, cholesky_factor_constrain) { | |||
auto f = [](int M, int N) { | |||
return [M, N](const auto& x1) { | |||
return stan::math::cholesky_factor_constrain(x1, M, N); | |||
stan::scalar_type_t<std::decay_t<decltype(x1)>> lp = 0.0; |
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.
[question]
Doesn't scalar_type_t
do all the decaying for you? If not, what's this for?
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.
Yep, fixed!
}; | ||
}; | ||
|
||
auto f2 = [](int M, int N) { | ||
return [M, N](const auto& x1) { | ||
stan::scalar_type_t<std::decay_t<decltype(x1)>> lp = 0.0; | ||
stan::math::cholesky_factor_constrain(x1, M, N, lp); | ||
stan::math::cholesky_factor_constrain<true>(x1, M, N, 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.
[question]
Are we still testing cholesky_factor_constrain
? Given that it's still in the API, it should still be tested. Especially since I've filed the issue to refactor the implementations.
Same question for all of these tests.
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.
Yep! See my comment here
}; | ||
auto f3 = [](const auto& x, const auto& lb, const auto& ub) { | ||
stan::return_type_t<decltype(x), decltype(lb), decltype(ub)> lp = 0; | ||
stan::math::lub_constrain(x, lb, ub, lp); | ||
stan::math::lub_constrain<true>(x, lb, ub, 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.
Isn't this redundant with the test on line 12? Are we still testing the old lub_constrain
? We should be.
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.
Yes so for each of these tests corresponds to the constraint without lp and is for the constrain with lp. In order to test the API I thought it would just be easier to change all of these to use it, though eod it still differs to the original functions the test was checking
…emplate-constrains
…4.1 (tags/RELEASE_600/final)
…emplate-constrains
Jenkins Console Log Machine informationProductName: Mac OS X ProductVersion: 10.11.6 BuildVersion: 15G22010CPU: G++: Clang: |
Is this ready for a final review? I should be able to get to it tomorrow (Wednesday). |
Yep! Though I had two doc qs for the below |
Jenkins Console Log Machine informationProductName: Mac OS X ProductVersion: 10.11.6 BuildVersion: 15G22010CPU: G++: Clang: |
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.
Looks good. Thanks! I left some residual comments, but no to-do items.
template <bool Jacobian, typename T> | ||
inline auto cholesky_factor_constrain(const T& x, int M, int N, | ||
scalar_type_t<T>& lp) { | ||
if (Jacobian) { |
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.
Interesting. Does the ternary operator not get reduced at compile time if the condition is constexpr
? Obviously the template parameters are all constexpr
!
@@ -49,6 +49,31 @@ inline auto corr_constrain(const T_x& x, T_lp& lp) { | |||
return tanh_x; | |||
} | |||
|
|||
/** | |||
* Return the result of transforming the specified scalar or container of values |
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.
Looks good. Thanks.
* scalar type. | ||
* @param y vector of K unrestricted variables | ||
* @return Unit length vector of dimension K | ||
*/ | ||
template <typename T, require_eigen_col_vector_t<T>* = nullptr, | ||
require_not_vt_autodiff<T>* = nullptr> | ||
inline auto unit_vector_constrain(const T& y) { | ||
inline plain_type_t<T> unit_vector_constrain(const T& y) { | ||
using std::sqrt; | ||
check_nonzero_size("unit_vector_constrain", "y", 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.
No. The idea is that y[N] = sqrt(1 - sum(y[1:N-1]))
. So if N = 1, the sum is zero, and we're left with sqrt(1)
. That's good, because [1]'
is a unit vector, but [0]'
is not!
I opened the new issue here: #2569
Summary
Related to stan-dev/stanc3#947 this PR allows stanc3 to pass a bool template parameter to the constraints to decide whether or not they should accumulate the jacobian calculation into the
lp
parameter. See the PR above to note see examples of how this is used.All I really did here was add a new overload to the
*_constrain
functions that has abool Jacobian
template parameter. If that is flipped totrue
then we do the jacobian calculation.Tests
I changed all of the current constrain tests to use the new
Jacobian
interface, which should still cover all the current test examples and the new API.Side Effects
Release notes
Adds an overload for the constrain functions on whether to accumulate jacobians into log probability argument.
Checklist
Math issue #(issue number)
Copyright holder: Steve Bronder
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