Allow constrain and free functions to work on Eigen types and std vectors #1766

SteveBronder · 2020-03-09T04:08:55Z

Summary

I was working on a PR for sparse matrices up in the Stan library and found a few places where we can use vectorized ops in the reader and writer. This PR adds vectorized versions of the *_constrain and *_free functions so they can take in eigen types and standard vectors. This also adds testing for the *_free functions which idt existed before.

Tests

Tests can be run with

./runTests.py ./test/unit/math/mix/fun/ -f constrain_test

Side Effects

The *_constrain and *_free can now take in vector like types.

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
- unit tests pass (to run, use: ./runTests.py test/unit)
- header checks pass, (make test-headers)
- dependencies checks pass, (make test-math-dependencies)
- docs build, (make doxygen)
- code passes the built in C++ standards checks (make cpplint)
the code is written in idiomatic C++ and changes are documented in the doxygen
the new changes are tested

…stable/2017-11-14)

stan-buildbot · 2020-03-09T12:15:40Z

Name	Old Result	New Result	Ratio	Performance change( 1 - new / old )
gp_pois_regr/gp_pois_regr.stan	4.86	4.8	1.01	1.32% faster
low_dim_corr_gauss/low_dim_corr_gauss.stan	0.02	0.02	0.97	-3.37% slower
eight_schools/eight_schools.stan	0.09	0.09	0.97	-2.71% slower
gp_regr/gp_regr.stan	0.22	0.22	1.01	1.46% faster
irt_2pl/irt_2pl.stan	6.48	6.45	1.01	0.52% faster
performance.compilation	87.35	86.22	1.01	1.29% faster
low_dim_gauss_mix_collapse/low_dim_gauss_mix_collapse.stan	7.54	7.58	1.0	-0.48% slower
pkpd/one_comp_mm_elim_abs.stan	20.83	19.99	1.04	4.02% faster
sir/sir.stan	93.9	91.58	1.03	2.48% faster
gp_regr/gen_gp_data.stan	0.05	0.05	0.97	-2.75% slower
low_dim_gauss_mix/low_dim_gauss_mix.stan	2.95	2.98	0.99	-1.04% slower
pkpd/sim_one_comp_mm_elim_abs.stan	0.31	0.33	0.94	-6.78% slower
arK/arK.stan	1.9	1.75	1.09	7.99% faster
arma/arma.stan	0.69	0.66	1.04	4.23% faster
garch/garch.stan	0.52	0.62	0.84	-19.18% slower
Mean result: 0.994623494379

Jenkins Console Log
Blue Ocean
Commit hash: 72c877f

Machine information

ProductName: Mac OS X ProductVersion: 10.11.6 BuildVersion: 15G22010

CPU:
Intel(R) Xeon(R) CPU E5-1680 v2 @ 3.00GHz

G++:
Configured with: --prefix=/Applications/Xcode.app/Contents/Developer/usr --with-gxx-include-dir=/usr/include/c++/4.2.1
Apple LLVM version 7.0.2 (clang-700.1.81)
Target: x86_64-apple-darwin15.6.0
Thread model: posix

Clang:
Apple LLVM version 7.0.2 (clang-700.1.81)
Target: x86_64-apple-darwin15.6.0
Thread model: posix

SteveBronder · 2020-03-09T19:08:15Z

This is ready for review! If the reviewer wants me to break this up into multiple PRs that's fine, though the changes are mostly the same across a lot of the files

t4c1

This PR is indeed quite large. It is managable, but next time plan for smaller ones. I tried not to repeat comments, so most of them apply to multiple places in code.

I think it only makes sense to use auto to replace types or metaprograms that are too long or confusing. it is subjective, where exactely is the threshold is, but built-in types (such as int) should never be replaced by auto.

t4c1 · 2020-03-13T11:33:28Z

stan/math/prim/err/check_ordered.hpp

-                   const std::vector<T_y>& y) {
-  for (size_t n = 1; n < y.size(); n++) {
+template <typename Vec, require_vector_like_t<Vec>* = nullptr>
+inline void check_ordered(const char* function, const char* name, Vec&& y) {


Arguments should be passed into functions as const references unless there is a good reason to pass them in some other way.

What's the reason behind prefering const&? I kind of like universal references since they preserve value types

It is simply the matter of using the simplest tool that does the job. It makes code easier to understand (especially to someone new to the codebase). By specifying const& anyone reading the code immediately knows that the function does not modify the argument.

What do you mean by preserving value types? In what cases is it beneficial?

I'll second @t4c1's comment. This is purely for readability. Yes, you can use universal references, but we don't have that propagated throughout the Math library. This PR doesn't need this particular change and if you wanted to change that all at once, we can do that (in a separate issue + PR).

t4c1 · 2020-03-13T11:34:49Z

stan/math/prim/err/check_ordered.hpp

-  for (size_t n = 1; n < y.size(); n++) {
+template <typename Vec, require_vector_like_t<Vec>* = nullptr>
+inline void check_ordered(const char* function, const char* name, Vec&& y) {
+  for (auto n = 1; n < y.size(); n++) {


No reason to use auto here. int is both shorter and clearer.

t4c1 · 2020-03-13T11:42:57Z

stan/math/prim/fun/divide.hpp

+inline auto divide(const Vec& x, Scal c) {
+  std::vector<value_type_t<Vec>> ret_x(x.size());
+  std::transform(x.begin(), x.end(), ret_x.begin(),
+                 [&c](auto&& x_iter) { return x_iter / c; });


This could use apply_vector_unary and be combined with previous overload.

Which overload are you referencing here?

template <typename Mat, typename Scal, typename = require_eigen_t<Mat>, typename = require_stan_scalar_t<Scal>, typename = require_all_not_var_t<scalar_type_t<Mat>, Scal>> inline auto divide(const Mat& m, Scal c) { return (m / c).eval(); }

@SteveBronder, did you see @t4c1's last comment here?

t4c1 · 2020-03-13T11:45:26Z

stan/math/prim/fun/dot_self.hpp

+ * @throw std::domain_error If v is not vector dimensioned.
+ */
+template <typename StdVec, require_std_vector_t<StdVec>* = nullptr>
+inline auto dot_self(StdVec&& x) {


[optional] In my opinion

template <Scalar> inline auto dot_self(const std::vector<Scalar>& x)

is easier to understand than using require, while producing the same result.

Also you can use apply_vector_unary to combine this overload with the next one.

I saw you changed this in your PR so I'll wait on these till your PR is merged

t4c1 · 2020-03-13T11:48:28Z

stan/math/prim/fun/dot_self.hpp

+template <typename StdVec, require_std_vector_t<StdVec>* = nullptr>
+inline auto dot_self(StdVec&& x) {
+  value_type_t<StdVec> sum = 0.0;
+  for (auto&& i : x) {


I would prefere value_type_t<StdVec> here. In this instance I might be ok with auto, but since i is scalar we don't need reference (especially not universal reference).

t4c1 · 2020-03-13T12:05:09Z

stan/math/prim/fun/ordered_constrain.hpp


-  size_type k = x.size();
-  Matrix<T, Dynamic, 1> y(k);
+  auto k = x.size();


No need for auto here.

I thought so too but it looks like users can change this and it's default is

EIGEN_DEFAULT_DENSE_INDEX_TYPE - the type for column and row indices in matrices, vectors and array (DenseBase::Index). Set to std::ptrdiff_t by default.

So I know it sounds weird but I kind of like auto here. Though casting to int is no biggie (idt) so I'm fine with just calling it int

https://eigen.tuxfamily.org/dox/TopicPreprocessorDirectives.html

In that case you can use Eigen::Index to avoid the cast.

t4c1 · 2020-03-13T12:10:59Z

stan/math/prim/fun/positive_ordered_constrain.hpp

-  size_type k = x.size();
-  Matrix<T, Dynamic, 1> y(k);
+  auto k = x.size();
+  plain_type_t<Vec> y(k);
  if (k == 0) {
    return y;
  }
  y[0] = exp(x[0]);
-  for (size_type i = 1; i < k; ++i) {
+  for (auto i = 1; i < k; ++i) {
    y[i] = y[i - 1] + exp(x[i]);
  }
  return y;


return cumulative_sum(exp(x))?

I need to look into this but that was giving me the test failure

test/unit/math/prim/fun/positive_ordered_transform_test.cpp:25: Failure Expected equality of these values: exp(1.0) + exp(-2.0) + exp(-5.0) Which is: 2.86036 y[2] Which is: 2.86036

Well the test uses EXPECT_EQ for floating point numbers. It would fail for anything that is not exactely the same expression. The test should be modified to have some tolerance.

t4c1 · 2020-03-13T12:12:49Z

stan/math/prim/fun/simplex_constrain.hpp

@@ -20,28 +20,24 @@ namespace math {
 *
 * The transform is based on a centered stick-breaking process.
 *
- * @tparam T type of elements in the vector
+ * @tparam Vec type deriving from `Eigen::MatgrixBase` with rows or columns


Suggested change

* @tparam Vec type deriving from `Eigen::MatgrixBase` with rows or columns

* @tparam Vec type deriving from `Eigen::MatrixBase` with rows or columns

t4c1 · 2020-03-13T12:26:20Z

stan/math/prim/fun/vec_concat.hpp

+  std::vector<value_type_t<Vec>> vec
+      = vec_concat(std::forward<VecArgs>(args)...);


[optional] This is a bit inefficient. The result vector vec should reserve its final size in advance.

t4c1 · 2020-03-13T12:29:28Z

test/unit/math/mix/fun/identity_constrain_test.cpp

+auto g3(const T& x) {
+  stan::value_type_t<T> lp = 0;
+  auto x_cons = stan::math::identity_constrain(x, lp);
+  auto x_free = stan::math::identity_free(x_cons);


You are only testing identity_constrain and identity_free together. Each should also be tested individually.

bob-carpenter · 2020-03-13T14:06:33Z

stan/math/prim/fun/dot_self.hpp

-  check_vector("dot_self", "v", v);
+template <typename EigMat,
+          require_eigen_vt<std::is_arithmetic, EigMat>* = nullptr>
+inline double dot_self(EigMat&& v) {


I'm not reviewing, but would very much appreciate seeing this renamed to EigVec (what I'd really like is V, but that's a bigger discussion we need to take up).

It is not an Eigen matrix, it's an Eigen vector. Eigen::Matrix<T, 1, -1> and Eigen::Matrix<T, -1, 1> and Eigen::Matrix<T, -1, -1> are three completely independent specializations of Eigen::Matrix<T, R, C>. Only the last of the tree is a matrix. The other two are vectors and do not even support T& operator()(ptrdiff_t, ptrdiff_t);.

I find it weird, but there are even some test checking that dot_self works with (non-vector) matrices.

bob-carpenter · 2020-03-13T14:11:00Z

stan/math/prim/fun/identity_free.hpp

@@ -13,13 +13,32 @@ namespace math {
 * <p>This function is a no-op and mainly useful as a placeholder
 * in auto-generated code.
 *
- * @tparam T type of value
- * @param[in] y value
+ * @tparam Scalar type of value


This is not a scalar any more, because now scalars include complex numbers. Going back to T works. Or if you insist on verbosely naming types, it should be Real.

Why not? Complex scalars are scalars as well.

The problem's the other way around. These functions won't work for complex inputs, so the required type should be called Real, not Scalar.

… doing constraints

…xp1~20180509124008.99 (branches/release_50)

…eanup/constrains

…xp1~20180509124008.99 (branches/release_50)

…eanup/constrains

…gs/RELEASE_500/final)

stan-buildbot · 2020-03-18T01:24:35Z

Name	Old Result	New Result	Ratio	Performance change( 1 - new / old )
gp_pois_regr/gp_pois_regr.stan	4.96	4.87	1.02	1.95% faster
low_dim_corr_gauss/low_dim_corr_gauss.stan	0.02	0.02	0.97	-2.71% slower
eight_schools/eight_schools.stan	0.09	0.09	1.01	0.59% faster
gp_regr/gp_regr.stan	0.22	0.22	1.0	0.03% faster
irt_2pl/irt_2pl.stan	6.5	6.44	1.01	0.84% faster
performance.compilation	88.52	86.1	1.03	2.73% faster
low_dim_gauss_mix_collapse/low_dim_gauss_mix_collapse.stan	7.53	7.53	1.0	0.07% faster
pkpd/one_comp_mm_elim_abs.stan	21.06	21.29	0.99	-1.09% slower
sir/sir.stan	93.3	90.74	1.03	2.74% faster
gp_regr/gen_gp_data.stan	0.05	0.05	0.99	-1.09% slower
low_dim_gauss_mix/low_dim_gauss_mix.stan	2.95	2.99	0.98	-1.54% slower
pkpd/sim_one_comp_mm_elim_abs.stan	0.31	0.31	1.0	-0.15% slower
arK/arK.stan	1.74	1.73	1.0	0.48% faster
arma/arma.stan	0.66	0.66	1.0	-0.41% slower
garch/garch.stan	0.52	0.62	0.84	-19.17% slower
Mean result: 0.991128698371

Jenkins Console Log
Blue Ocean
Commit hash: eb24ce2

Machine information

ProductName: Mac OS X ProductVersion: 10.11.6 BuildVersion: 15G22010

CPU:
Intel(R) Xeon(R) CPU E5-1680 v2 @ 3.00GHz

G++:
Configured with: --prefix=/Applications/Xcode.app/Contents/Developer/usr --with-gxx-include-dir=/usr/include/c++/4.2.1
Apple LLVM version 7.0.2 (clang-700.1.81)
Target: x86_64-apple-darwin15.6.0
Thread model: posix

Clang:
Apple LLVM version 7.0.2 (clang-700.1.81)
Target: x86_64-apple-darwin15.6.0
Thread model: posix

syclik

@SteveBronder, it looks like there's work to do on your part for this PR.

syclik · 2020-04-28T03:44:02Z

stan/math/prim/err/check_ordered.hpp

-                   const std::vector<T_y>& y) {
-  for (size_t n = 1; n < y.size(); n++) {
+template <typename Vec, require_vector_like_t<Vec>* = nullptr>
+inline void check_ordered(const char* function, const char* name, Vec&& y) {


I'll second @t4c1's comment. This is purely for readability. Yes, you can use universal references, but we don't have that propagated throughout the Math library. This PR doesn't need this particular change and if you wanted to change that all at once, we can do that (in a separate issue + PR).

syclik · 2020-04-28T03:45:03Z

stan/math/prim/err/check_positive_ordered.hpp

-
+template <typename Vec, require_vector_like_t<Vec>* = nullptr>
+inline void check_positive_ordered(const char* function, const char* name,
+                                   Vec&& y) {


Same comment about using the "universal reference."

syclik · 2020-04-28T03:47:07Z

stan/math/prim/fun/divide.hpp

+inline auto divide(const Vec& x, Scal c) {
+  std::vector<value_type_t<Vec>> ret_x(x.size());
+  std::transform(x.begin(), x.end(), ret_x.begin(),
+                 [&c](auto&& x_iter) { return x_iter / c; });


@SteveBronder, did you see @t4c1's last comment here?

syclik · 2020-04-28T03:48:40Z

stan/math/prim/fun/identity_constrain.hpp

-inline T identity_constrain(const T& x) {
-  return x;
+template <typename Scalar, require_all_stan_scalar_t<Scalar>* = nullptr>
+inline decltype(auto) identity_constrain(Scalar&& x) {


Any results from benchmarking?

…4.1 (tags/RELEASE_600/final)

…eanup/constrains

…4.1 (tags/RELEASE_600/final)

…eanup/constrains

…4.1 (tags/RELEASE_600/final)

SteveBronder · 2020-05-11T18:23:45Z

I got this to work on an upstream stan branch and didn't see any performance improvements over the basic benchmarks and a few in examples so going to close

SteveBronder added 3 commits March 8, 2020 22:23

make constrained types vectorized and accept any vector type

80cb1f2

restrain simplex to work with only eigen vectors

858315f

update docs

35474e2

SteveBronder changed the title ~~Cleanup/constrains~~ Allow constrain and free functions to work in Eigen types and std vectors Mar 9, 2020

SteveBronder changed the title ~~Allow constrain and free functions to work in Eigen types and std vectors~~ Allow constrain and free functions to work on Eigen types and std vectors Mar 9, 2020

stan-buildbot and others added 5 commits March 9, 2020 00:09

[Jenkins] auto-formatting by clang-format version 6.0.0 (tags/google/…

099e60d

…stable/2017-11-14)

Include what you use

c2b457d

test-headers

1ee05ab

Fix dot_self template

8f63917

use eigen_vt instead of eigen_vector_vt for dot_self in prim

72c877f

t4c1 requested changes Mar 13, 2020

View reviewed changes

bob-carpenter reviewed Mar 13, 2020

View reviewed changes

SteveBronder and others added 11 commits March 14, 2020 17:31

Merge remote-tracking branch 'origin/develop' into cleanup/constrains

cd23b53

make the vector and eigen constrain functions check the bounds before…

2712a1b

… doing constraints

[Jenkins] auto-formatting by clang-format version 5.0.2-svn328729-1~e…

ce83561

…xp1~20180509124008.99 (branches/release_50)

cpplint

edf5df3

Merge branch 'cleanup/constrains' of github.com:stan-dev/math into cl…

439223b

…eanup/constrains

eigen and std vec for lub_free

7ef27ae

[Jenkins] auto-formatting by clang-format version 5.0.2-svn328729-1~e…

4ef0aab

…xp1~20180509124008.99 (branches/release_50)

update docs

2709fd6

Merge branch 'cleanup/constrains' of github.com:stan-dev/math into cl…

1a06e87

…eanup/constrains

[Jenkins] auto-formatting by clang-format version 5.0.0-3~16.04.1 (ta…

88e7366

…gs/RELEASE_500/final)

cpplint header fix

eb24ce2

syclik reviewed Apr 28, 2020

View reviewed changes

SteveBronder added 2 commits May 11, 2020 12:20

merge to develop

5c6e7e2

remove some const refs

95365e8

stan-buildbot and others added 8 commits May 11, 2020 16:27

[Jenkins] auto-formatting by clang-format version 6.0.0-1ubuntu2~16.0…

0659770

…4.1 (tags/RELEASE_600/final)

merge with current check_simplex

ccb27ac

Merge branch 'cleanup/constrains' of github.com:stan-dev/math into cl…

15c1cd9

…eanup/constrains

remove added tests for checking against vector for pos_ordered_constrain

6842e6b

[Jenkins] auto-formatting by clang-format version 6.0.0-1ubuntu2~16.0…

e5aecfd

…4.1 (tags/RELEASE_600/final)

allow for fma to take in an eigen type

9bc422f

Merge branch 'cleanup/constrains' of github.com:stan-dev/math into cl…

8c44eaa

…eanup/constrains

[Jenkins] auto-formatting by clang-format version 6.0.0-1ubuntu2~16.0…

f683957

…4.1 (tags/RELEASE_600/final)

SteveBronder closed this May 11, 2020

	* @tparam Vec type deriving from `Eigen::MatgrixBase` with rows or columns
	* @tparam Vec type deriving from `Eigen::MatrixBase` with rows or columns

		std::vector<value_type_t<Vec>> vec
		= vec_concat(std::forward<VecArgs>(args)...);

Allow constrain and free functions to work on Eigen types and std vectors #1766

Allow constrain and free functions to work on Eigen types and std vectors #1766

Conversation

SteveBronder commented Mar 9, 2020

Summary

Tests

Side Effects

Checklist

stan-buildbot commented Mar 9, 2020

SteveBronder commented Mar 9, 2020

t4c1 left a comment • edited Loading

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

stan-buildbot commented Mar 18, 2020

syclik left a comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

SteveBronder commented May 11, 2020

t4c1 left a comment •

edited

Loading