Vectorize everything with broadcasting

[bob-carpenter]

Vectorize all functions for matching shapes and allow broadcasting of lower-dimensional types.

This will require

• the function to be vectorized define a core functor-like static class with templated operator and appropriate typedefs
• a template function that can accept the core functor class as a template argument and implement the appropriate functions with template programs

There is an example for the function inv() in stan/math/prim/mat/fun/inv.hpp in the branch feature/issue-202-vectorize-all.

This issue involves the first step above for all of the Stan functions, so it may need to be broken down into stages, starting with unary functions.

For the stan-dev/stan issue: stan-dev/stan#1683

For each of thee functions, we need (1) functor, (2) general template definition, (3) doc in manual, (4) signatures [need function for this], (5) signature tests for all instantiations.

When I did the example for inv(), I had to remove the templating for inv() itself and replace it with double to remove the ambiguity. For functions that only have a template version (i.e., there is not also a version for fwd and rev), I believe there will need to be an enable_if on the base definition that restricts application to the base class.

We'll also need a testing framework. This should be easy to do following the way apply_scalar_unary is defined (across prim, rev, and fwd). It will require a simple templated testing functor for each function being tested.

(int):int

• abs
• int_step
• operator!
• operator+
• operator-

(int,int):int

• max
• min
• operator*
• operator+
• operator-
• operator/
• operator<
• operator<=
• operator>
• operator>=
• operator==
• operator&&
• operator||

(int,real):real

• bessel_first_kind
• bessel_second_kind
• modified_bessel_first_kind
• modified_bessel_second_kind

(real):int

• int_step
• is_inf
• is_nan
• operator!

(real):real

• abs
• acos
• acosh
• asin
• asinh
• atan
• atanh
• cbrt
• ceil
• cos
• cosh
• digamma
• erf
• erfc
• exp
• exp2
• expm1
• fabs
• floor
• inv
• inv_cloglog
• inv_logit
• inv_Phi
• inv_sqrt
• inv_square
• lgamma
• log
• log10
• log1m
• log1m_exp
• log1m_inv_logit
• log1p
• log1p_exp
• log2
• log_inv_logit
• logit
• Phi
• Phi_approx
• round
• sin
• sinh
• sqrt
• square
• step
• tan
• tanh
• tgamma
• trigamma
• trunc

(real,real):int

• operator<
• operator<=
• operator>
• operator>=
• operator==
• operator&&
• operator||

(real,real):real

• atan2
• binomial_coefficient_log
• falling_factorial
• fdim
• fmax
• fmin
• fmod
• gamma_p
• gamma_q
• hypot
• lbeta
• log_diff_exp
• log_falling_factorial
• log_rising_factorial
• log_sum_exp
• max
• min
• multiply_log
• operator*
• operator+
• operator-
• operator/
• operator^
• owens_t
• pow
• rising_factorial

(int, real):real

• binary_log_loss
• lmgamma

(real,real,real):real

• fma
• log_mix

(int,real,real):real

• if_else

[rayleigh]

I'm reopening this issue because the pull request that was merged in only provides the infrastructure to test these vectorize functions. However, if you think that we should close this issue and open a new one, I can do that.

[bob-carpenter]

Either way is OK. You need to make the checklist of all
unary functions to which this applies. You can probably knock
them off in multiple issues.

On Mar 15, 2016, at 4:52 PM, rayleigh [email protected] wrote:

I'm reopening this issue because the pull request that was merged in only provides the infrastructure to test these vectorize functions. However, if you think that we should close this issue and open a new one, I can do that.

—
You are receiving this because you were assigned.
Reply to this email directly or view it on GitHub

[rayleigh]

Okay. Copying over your response to a question of where to put files:

1. The function definitions so should go here:

stan/math/prim/mat/fun/<function-name>.hpp

1. For functions already partly vectorized, you need
   to remove the original vectorized definitions (just get
   rid of them and their associated tests).
2. The function signatures class is going to need a utility
   method to declare these, presumably something like
   ``add_unary_vectorized("exp") and so on.

This should declare all the basic vectorizations and
perhaps up to 4 deep for arrays of all types.

1. You're going to need to write (or have me write)
   a. a new intro section to the functions guide in the manual,
   with some new syntax for this vectorization, and
   b. update the doc for each of the added functions

Also, where should the test files go? Should they go in /test/unit/math/mix/mat/fun/<function-name>_test.hpp?

[bob-carpenter]

Thanks and yes, that's the right place for the test files, because
they depend on mix.

On Mar 17, 2016, at 10:30 AM, rayleigh [email protected] wrote:

Okay. Copying over your response to a question of where to put files:

• The function definitions so should go here:
stan/math/prim/mat/fun/.hpp

• For functions already partly vectorized, you need
to remove the original vectorized definitions (just get
rid of them and their associated tests).

• The function signatures class is going to need a utility
method to declare these, presumably something like
``add_unary_vectorized("exp") and so on.

This should declare all the basic vectorizations and
perhaps up to 4 deep for arrays of all types.

• You're going to need to write (or have me write) a. a new intro section to the functions guide in the manual, with some new syntax for this vectorization, and b. update the doc for each of the added functions
Also, where should the test files go? Should they go in /test/unit/math/mix/mat/fun/_test.hpp?

—
You are receiving this because you were assigned.
Reply to this email directly or view it on GitHub

[rayleigh]

Thanks. Looking through the 2.8.0 Stan manual, I noticed that abs is to be deprecated. Should I vectorize it or should I only vectorize fabs?

[bob-carpenter]

Ack, I think we undeprecated abs and are going to
deprecate fabs (to make it more like math and less
like C++). So go ahead and do both and we can see where
the chips land.

On Mar 21, 2016, at 11:45 AM, rayleigh [email protected] wrote:

Thanks. Looking through the 2.8.0 Stan manual, I noticed that abs is to be deprecated. Should I vectorize it or should I only vectorize fabs?

—
You are receiving this because you were assigned.
Reply to this email directly or view it on GitHub

[rayleigh]

Sounds good; I vectorized both. I've been trying to vectorize the function cbrt. It's defined for 0, but its derivatives aren't (it returns NaN). From this, I realized that the testing framework assumed that if a function is defined for a value, its derivatives are as well. To handle this, does the vectorize function throw an error if an user enters 0? Or, do I write separate code to test this case?

[syclik]

I think we should test derivatives.

We had this discussion in the past regarding defined values and undefined
derivatives. I think the functions should trap that error and fail early if
possible, even if it conflicts with other implementations of the function.
If we're using these functions primarily where they need derivatives, it
makes no sense to knowingly propagate NaNs.

On Mon, Mar 21, 2016 at 7:00 PM, rayleigh [email protected] wrote:

Sounds good; I vectorized both. I've been trying to vectorize the function
cbrt. It's defined for 0, but its derivatives aren't (it returns NaN).
From this, I realized that the testing framework assumed that if a function
is defined for a value, its derivatives are as well. To handle this, does
the vectorize function throw an error if an user enters 0? Or, do I write
separate code to test this case?

—
You are receiving this because you modified the open/close state.
Reply to this email directly or view it on GitHub
#202 (comment)

[rayleigh]

Okay. Just to clarify, should the vectorize-all function or the individual scalar functions include the error checking? Or, both?

[syclik]

That was a general comment.

What test am I supposed to run? Can you provide the full path?

On Tue, Mar 22, 2016 at 11:34 AM, rayleigh [email protected] wrote:

Okay. Just to clarify, should the vectorize-all function or the individual
scalar functions include the error checking? Or, both?

—
You are receiving this because you modified the open/close state.
Reply to this email directly or view it on GitHub
#202 (comment)

[rayleigh]

I pushed more changes to the branch feature/issue-202-vectorize-all. However, I only wanted to commit things that worked. So, if you pull that branch, you'll need to change test/unit/math/mix/mat/fun/cbrt_test.cpp before running it. In the function valid_inputs, change any one of the values to 0.

[rayleigh]

I pushed the change to exp.hpp to feature/issue-202-vectorize-all. The test to run is test/unit/math/mix/mat/vectorize/foo_test.cpp.

[bob-carpenter]

Is it only exp() that's causing the issue?

• Bob

On Mar 22, 2016, at 2:28 PM, rayleigh [email protected] wrote:

I pushed the change to exp.hpp to feature/issue-202-vectorize-all. The test to run is test/unit/math/mix/mat/vectorize/foo_test.cpp.

—
You are receiving this because you were assigned.
Reply to this email directly or view it on GitHub

[syclik]

No. This isn't an include-order issue.

Here's where I'm at so far with exp:

1. There's a definition of a struct exp_fun in
   stan/math/prim/scal/mat/fun/exp.hpp.
2. There is a definition of a templated function exp that should handle the
   vectorization in stan/math/prim/mat/fun/exp.hpp.
3. There is a definition of a stan::math::var exp(const stan::math::var& a) in stan/math/rev/scal/fun/exp.hpp.
4. When calling foo_test.cpp, which should really be changed to
   exp_test.cpp, with the appropriate changes, this is the compiler error
   message that crops up:

In file included from test/unit/math/mix/mat/vectorize/foo_test.cpp:1:
In file included from ./stan/math/mix/mat.hpp:12:
In file included from ./stan/math/prim/mat.hpp:91:
In file included from ./stan/math/prim/mat/fun/exp.hpp:4:
./stan/math/prim/mat/vectorize/apply_scalar_unary.hpp:40:41: error:
implicit instantiation of
      undefined template 'Eigen::internal::traits<stan::math::var>'
      typedef typename Eigen::internal::traits<T>::Scalar scalar_t;
                                        ^
./stan/math/prim/mat/fun/exp.hpp:23:21: note: in instantiation of template
class
      'stan::math::apply_scalar_unary<stan::math::exp_fun,
stan::math::var>' requested here
    inline typename apply_scalar_unary<exp_fun, T>::return_t
                    ^
./stan/math/rev/scal/fun/grad_inc_beta.hpp:28:16: note: while substituting
deduced template arguments
      into function template 'exp' [with T = stan::math::var]
      var c3 = exp(lbeta(a, b)) * inc_beta(a, b, z);
               ^
lib/eigen_3.2.4/Eigen/src/Core/util/ForwardDeclarations.h:17:29: note:
template is declared here
template<typename T> struct traits;
                            ^
In file included from test/unit/math/mix/mat/vectorize/foo_test.cpp:1:
In file included from ./stan/math/mix/mat.hpp:12:
In file included from ./stan/math/prim/mat.hpp:91:
In file included from ./stan/math/prim/mat/fun/exp.hpp:4:

1. That compiler error is coming from the template specialization of
   apply_scalar_unary() not being used from
   stan/math/rev/mat/vectorize/apply_scalar_unary.hpp.

I think if I can figure out why that specialization is not being used, it
would go a long way.

Daniel

On Tue, Mar 22, 2016 at 11:09 PM, Bob Carpenter [email protected]
wrote:

Is it only exp() that's causing the issue?

• Bob

On Mar 22, 2016, at 2:28 PM, rayleigh [email protected] wrote:

I pushed the change to exp.hpp to feature/issue-202-vectorize-all. The
test to run is test/unit/math/mix/mat/vectorize/foo_test.cpp.

—
You are receiving this because you were assigned.
Reply to this email directly or view it on GitHub

—
You are receiving this because you modified the open/close state.
Reply to this email directly or view it on GitHub
#202 (comment)

[syclik]

Ok. Those headers aren't included in stan/math/prim/mat.hpp or rev/mat.hpp
or fwd/mat.hpp.

I'll fix that now.

On Tue, Mar 22, 2016 at 11:22 PM, Daniel Lee [email protected] wrote:

No. This isn't an include-order issue.

Here's where I'm at so far with exp:

1. There's a definition of a struct exp_fun in
   stan/math/prim/scal/mat/fun/exp.hpp.
2. There is a definition of a templated function exp that should handle
   the vectorization in stan/math/prim/mat/fun/exp.hpp.
3. There is a definition of a stan::math::var exp(const stan::math::var& a) in stan/math/rev/scal/fun/exp.hpp.
4. When calling foo_test.cpp, which should really be changed to
   exp_test.cpp, with the appropriate changes, this is the compiler error
   message that crops up:

In file included from test/unit/math/mix/mat/vectorize/foo_test.cpp:1:
In file included from ./stan/math/mix/mat.hpp:12:
In file included from ./stan/math/prim/mat.hpp:91:
In file included from ./stan/math/prim/mat/fun/exp.hpp:4:
./stan/math/prim/mat/vectorize/apply_scalar_unary.hpp:40:41: error:
implicit instantiation of
      undefined template 'Eigen::internal::traits<stan::math::var>'
      typedef typename Eigen::internal::traits<T>::Scalar scalar_t;
                                        ^
./stan/math/prim/mat/fun/exp.hpp:23:21: note: in instantiation of template
class
      'stan::math::apply_scalar_unary<stan::math::exp_fun,
stan::math::var>' requested here
    inline typename apply_scalar_unary<exp_fun, T>::return_t
                    ^
./stan/math/rev/scal/fun/grad_inc_beta.hpp:28:16: note: while substituting
deduced template arguments
      into function template 'exp' [with T = stan::math::var]
      var c3 = exp(lbeta(a, b)) * inc_beta(a, b, z);
               ^
lib/eigen_3.2.4/Eigen/src/Core/util/ForwardDeclarations.h:17:29: note:
template is declared here
template<typename T> struct traits;
                            ^
In file included from test/unit/math/mix/mat/vectorize/foo_test.cpp:1:
In file included from ./stan/math/mix/mat.hpp:12:
In file included from ./stan/math/prim/mat.hpp:91:
In file included from ./stan/math/prim/mat/fun/exp.hpp:4:

1. That compiler error is coming from the template specialization of
   apply_scalar_unary() not being used from
   stan/math/rev/mat/vectorize/apply_scalar_unary.hpp.

I think if I can figure out why that specialization is not being used, it
would go a long way.

Daniel

On Tue, Mar 22, 2016 at 11:09 PM, Bob Carpenter [email protected]
wrote:

Is it only exp() that's causing the issue?

• Bob

On Mar 22, 2016, at 2:28 PM, rayleigh [email protected] wrote:

I pushed the change to exp.hpp to feature/issue-202-vectorize-all. The
test to run is test/unit/math/mix/mat/vectorize/foo_test.cpp.

—
You are receiving this because you were assigned.
Reply to this email directly or view it on GitHub

—
You are receiving this because you modified the open/close state.
Reply to this email directly or view it on GitHub
#202 (comment)

[syclik]

I just pushed, but it's not working yet. Bob, I might need your help
reasoning through this.

apply_scalar_unary() in prim is being called first. Even when the mix
header is included. That gets instantiated with stan::math::var as its
second argument instead of being found with argument dependent lookup. If
you want more detail as to the include order, try running the test (it's ok
if it doesn't compile), then open up
test/unit/math/mix/vectorize/foo_test.d. That'll show you the order in
which files get opened.

If you include the rev version as the first thing in the test file, it
works.

I'm not sure how to get past this problem. I'm sure we've figured this out
in the past.

On Tue, Mar 22, 2016 at 11:24 PM, Daniel Lee [email protected] wrote:

Ok. Those headers aren't included in stan/math/prim/mat.hpp or rev/mat.hpp
or fwd/mat.hpp.

I'll fix that now.

On Tue, Mar 22, 2016 at 11:22 PM, Daniel Lee [email protected] wrote:

No. This isn't an include-order issue.

Here's where I'm at so far with exp:

1. There's a definition of a struct exp_fun in
   stan/math/prim/scal/mat/fun/exp.hpp.
2. There is a definition of a templated function exp that should handle
   the vectorization in stan/math/prim/mat/fun/exp.hpp.
3. There is a definition of a stan::math::var exp(const stan::math::var& a) in stan/math/rev/scal/fun/exp.hpp.
4. When calling foo_test.cpp, which should really be changed to
   exp_test.cpp, with the appropriate changes, this is the compiler error
   message that crops up:

In file included from test/unit/math/mix/mat/vectorize/foo_test.cpp:1:
In file included from ./stan/math/mix/mat.hpp:12:
In file included from ./stan/math/prim/mat.hpp:91:
In file included from ./stan/math/prim/mat/fun/exp.hpp:4:
./stan/math/prim/mat/vectorize/apply_scalar_unary.hpp:40:41: error:
implicit instantiation of
      undefined template 'Eigen::internal::traits<stan::math::var>'
      typedef typename Eigen::internal::traits<T>::Scalar scalar_t;
                                        ^
./stan/math/prim/mat/fun/exp.hpp:23:21: note: in instantiation of
template class
      'stan::math::apply_scalar_unary<stan::math::exp_fun,
stan::math::var>' requested here
    inline typename apply_scalar_unary<exp_fun, T>::return_t
                    ^
./stan/math/rev/scal/fun/grad_inc_beta.hpp:28:16: note: while
substituting deduced template arguments
      into function template 'exp' [with T = stan::math::var]
      var c3 = exp(lbeta(a, b)) * inc_beta(a, b, z);
               ^
lib/eigen_3.2.4/Eigen/src/Core/util/ForwardDeclarations.h:17:29: note:
template is declared here
template<typename T> struct traits;
                            ^
In file included from test/unit/math/mix/mat/vectorize/foo_test.cpp:1:
In file included from ./stan/math/mix/mat.hpp:12:
In file included from ./stan/math/prim/mat.hpp:91:
In file included from ./stan/math/prim/mat/fun/exp.hpp:4:

1. That compiler error is coming from the template specialization of
   apply_scalar_unary() not being used from
   stan/math/rev/mat/vectorize/apply_scalar_unary.hpp.

I think if I can figure out why that specialization is not being used, it
would go a long way.

Daniel

On Tue, Mar 22, 2016 at 11:09 PM, Bob Carpenter <[email protected]

wrote:

Is it only exp() that's causing the issue?

• Bob

On Mar 22, 2016, at 2:28 PM, rayleigh [email protected]
wrote:

I pushed the change to exp.hpp to feature/issue-202-vectorize-all. The
test to run is test/unit/math/mix/mat/vectorize/foo_test.cpp.

—
You are receiving this because you were assigned.
Reply to this email directly or view it on GitHub

—
You are receiving this because you modified the open/close state.
Reply to this email directly or view it on GitHub
#202 (comment)

[rayleigh]

Any luck getting past this problem?

[rayleigh]

Copying this over from an email chain because it might have the solution (I haven't tried yet):

If you edit out the very first include line in:

/test/unit/math/mix/mat/vectorize/foo_test.cpp,

everything passes. So the new file looks like:

---------------------
// #include <stan/math/mix/mat.hpp>
#include <gtest/gtest.h>
#include <test/unit/math/prim/mat/vectorize/prim_scalar_unary_test.hpp>
#include <test/unit/math/rev/mat/vectorize/rev_scalar_unary_test.hpp>
#include <test/unit/math/fwd/mat/vectorize/fwd_scalar_unary_test.hpp>
#include <test/unit/math/mix/mat/vectorize/mix_scalar_unary_test.hpp>
#include <test/unit/math/prim/mat/vectorize/foo_fun.hpp>
#include <test/unit/math/prim/mat/vectorize/vector_builder.hpp>

/**
 * This is the structure for testing mock function foo (defined in the
 * testing framework).  See README.txt for more instructions.
 */
struct foo_test {
-------------------------------

This makes me think there may be a problem in the include order in mix/mat.hpp.

You might want to read this:

http://stackoverflow.com/questions/14626033/implicit-instantiation-of-undefined-template-class

It's all about steps 1, 2, and 3 in the answer.

[rayleigh]

Other questions:

1. I tried vectorizing trunc.hpp, but the matrix tests fail to run for it. When I looked at boost/math/special_functions/trunc.hpp, I saw the following comment:

// The following functions will not compile unless T has an
// implicit convertion to the integer types. 

So, is trunc.hpp designed to work on matrix?

2. According to the manual, the function step returns a real, but when I look at the prim/scal/fun/step.hpp, it returns an int. Should I change the function to return a real?

3. For cbrt, I know that I asked about it above, but I want to clarify what I should do about it. In the tests, should I allow the code to pass if both the expected and tested derivatives are NaN?

[rayleigh]

Listing why certain unary functions weren't vectorized:
Unvectorized acosh and atanh because error handling isn't the same across all modes. Asinh looks to have similar coding problems as acosh and atanh, but it's defined for all values so I'm leaving it vectorized.

Unvectorized exp2, log2, and log1p because of an ambiguity error when making Stan model's C++ code for untemplated prim versions.

trunc and step aren't vectorized because both aren't compatible with matrices.

Other (real):real functions aren't vectorized because they're missing a mode and their prim version can't be untemplated.

[bob-carpenter]

Let's just get what you have working in and then create new issues
for these other cases. The other four cases with solutions are
as follows:

1. acosh etc (varying failure conditions in prim/rev/fwd)
   generalize testing or change error handling.
2. exp2 etc (defined only in ::)
   just need to resolve instantiation ambiguity as we've done for other cases,
   perhaps with enable_if, perhaps with appropriate qualification in generated
   code (see existing case of max() in generator).
3. step etc. (integer returns)
   just need appropriate enable_if to remove matrix arguments;
   or they can be defined and just not included in function_signatures in
   the Stan language; there's no API problem with allowing Matrix<int,...>.
4. others (?) (missing fwd or rev, so need templated prim case)
   easy fix to define prim case for double and then directly define fwd
   and rev without templates (well, fwd needs a template, but the whole arg's
   not templated)
   a more general fix requiring less code dup would be to disable_if
   the templated prim case to only apply to scalars (which is easy in theory,
   but may be tricky in practice to deal with our include orders among
   prim/fwd/rev and scal/arr/mat.

• Bob

On Apr 11, 2016, at 2:50 PM, rayleigh [email protected] wrote:

Listing why certain unary functions weren't vectorized:
Unvectorized acosh and atanh because error handling isn't the same across all modes. Asinh looks to have similar coding problems as acosh and atanh, but it's defined for all values so I'm leaving it vectorized.

Unvectorized exp2, log2, and log1p because of an ambiguity error when making Stan model's C++ code for untemplated prim versions.

trunc and step aren't vectorized because both aren't compatible with matrices.

Other (real):real functions aren't vectorized because they're missing a mode and their prim version can't be untemplated.

—
You are receiving this because you were assigned.
Reply to this email directly or view it on GitHub

[rayleigh]

Sounds good though I disagree about the solution to 1:

1. In the rev version of acosh.hpp, you have:

namespace {
      class acosh_vari : public op_v_vari {
      public:
        acosh_vari(double val, vari* avi) :
          op_v_vari(val, avi) {
        }
        void chain() {
          avi_->adj_ += adj_ / std::sqrt(avi_->val_ * avi_->val_ - 1.0);
        }
      };
    }

In the rev version of digamma.hpp (which also uses code from boost/math/special_functions), you have

namespace {
      class digamma_vari : public op_v_vari {
      public:
        explicit digamma_vari(vari* avi) :
          op_v_vari(digamma(avi->val_), avi) {
        }
        void chain() {
          avi_->adj_ += adj_ * trigamma(avi_->val_);
        }
      };
    }

The vectorized version of digamma handles error consistently across all versions so I think the rev and fwd mode code needs to be fixed for acosh, atanh, and asinh.

[bob-carpenter]

Here's the actual function for acosh:

inline var acosh(const var& a) {
if (boost::math::isinf(a.val()) && a > 0.0)
return var(new acosh_vari(a.val(), a.vi_));
return var(new acosh_vari(::acosh(a.val()), a.vi_));
}

We could move the call to ::acosh() into acosh_vari constructor,
but that wouldn't change the error handling. Or is it the
conditioning on isinf that's the issue?

So I'm not sure what you're suggesting or where there is
currently an inconsistency.

• Bob

[rayleigh]

Thanks for pointing that out. I took another look and I think the includes might be the issue because for acosh.hpp, the includes are:

#include <math.h>
#include <stan/math/rev/core.hpp>
#include <boost/math/special_functions/fpclassify.hpp>
#include <cmath>

#ifdef _MSC_VER
#include <boost/math/special_functions/acosh.hpp>
using boost::math::acosh;
#endif

So, unless it's being compiled on a Visual C++ compiler, I don't think it'll use acosh.hpp from boost/math/special_functions, which is providing the error handling. Because Stan doesn't use C++11 and acosh.hpp is only available in C++11's cmath library, the vectorized version of acosh.hpp uses acosh.hpp from boost/math/special_functions as the base function. I think whether boost/math/special_functions/acosh.hpp is included explains the inconsistency in error handling that I'm seeing.

[bob-carpenter]

That conditional include is nasty. Maybe we should just use
boost::math::acosh independently of platform.

Is there ever a reason to include both <math.h> and ?
And doesn't that <math.h> include have to go last?

The error handling can be fixed by just adding it to the
code, even if we do use Boost's definition for computing
non-error inputs.

On Apr 11, 2016, at 5:40 PM, rayleigh [email protected] wrote:

Thanks for pointing that out. I took another look and I think the includes might be the issue because for acosh.hpp, the includes are:

#include <math.h>
#include <stan/math/rev/core.hpp>
#include <boost/math/special_functions/fpclassify.hpp>
#include

#ifdef _MSC_VER
#include <boost/math/special_functions/acosh.hpp>
using boost::math::acosh;
#endif

So, unless it's being compiled on a Visual C++ compiler, I don't think it'll use acosh.hpp from boost/math/special_functions, which is providing the error handling. Because Stan doesn't use C++11 and acosh.hpp is only available in C++11's cmath library, the vectorized version of acosh.hpp uses acosh.hpp from boost/math/special_functions as the base function. I think whether boost/math/special_functions/acosh.hpp is included explains the inconsistency in error handling that I'm seeing.

—
You are receiving this because you were assigned.
Reply to this email directly or view it on GitHub

[bob-carpenter]

@rayleigh Is the current checklist above up to date? Is there a reason the dozen or so (real): real signature ones aren't implemented?

[bob-carpenter]

Closing in favor of newer issue #347 given that it's partially done.

[jessexknight]

Can we reopen this? The logical operators (at least) seem to remain un-vectorized

Binary infix operator == with functional interpretation logical_eq requires arguments of primitive type (int or real), found left type=int[ ], right arg type=int.

Thanks,

[bob-carpenter]

@jessexknight Probably better to start a new issue specifically or logical operations, because it's not immediately obvious how they should behave. This issue was specifically for real-valued operations and used a specific vectorization assuming real valued functions with autodiff---here there are no derivatives because we get integer values out.

I can imagine two alternatives for the logical operators:

1. int logical_eq(reals x, reals y); which returns a single truth value conjoining the element wise result.

2. int[] logical_eq(reals x, reals y); which returns a container of results element wise.

Both approaches can broadcast scalars to containers. If we go with (2), then we probably want functions all() and any() like in NumPy that return the conjunction and disjunction of a container.

This is related to the issue of comparing containers to elements to return booleans, the way that R allows. That'd give us something like this

int[] xs = { 1, 0, 2, 3, 0, -1};
int[] ys = (xs == 0);   // after eval, ys = { 2, 5 }

[jessexknight]

Thanks, I've evidently added an issue...

Just want to add that I would strongly prefer (2) to allow more flexibility.

Not sure I follow the last issue you mentioned -- I'm not familiar with this behaviour in R, and personally wouldn't find this a priority vs the vectorization of the base logical functions.

[bob-carpenter]

Thanks for opening the issue. I edited it slightly to make it easier for our devs.

In R, you can do this:

> sex = c(0, 1, 1, 0, 0, 0, 1)
> age = c(23, 29, 30, 12, 15, 18)
> age[sex == 0]
[1] 23 12 15 18
> age[sex == 1]
[1] 29 30 NA
> age = c(23, 29, 30, 22, 25, 18, 31)
> sex = c(0,   1,  1 , 0,  0 , 0,  1)

> sex == 0
[1]  TRUE FALSE FALSE  TRUE  TRUE  TRUE FALSE

> age[sex == 0]
[1] 23 22 25 18

As you can see, it's useful for picking subgroups out of parallel sequences. But it relies on really odd behavior where if you give R a list of boolean arguments, it'll include the ones that are TRUE in the result. This doesn't make much sense, as TRUE evaluates to 1 and FALSE evaluates to 0, but you get very different results if you replace the booleans with integers here.

> age[c(1, 0, 0, 1, 1, 1, 0)]
[1] 23 23 23 23

This is a very R result in that it seems to just ignore the out of range 0 inputs!

[jessexknight]

Oh, I see. I think this logical indexing is available in Numpy and Matlab too. In fact, I think logical indexing can be faster in Numpy, besides the fact that 0 is not out of bounds ;)

I suppose this raises the idea of a logical data type in Stan, but I think this type of indexing might be one of the only use cases ...

[bob-carpenter]

logical indexing can be faster in Numpy

It's always going to be bound by having to evaluate the condition for every element of the container. You could potentially do it without constructing the intermediate sex == 0 array by evaluating it lazily with an expression template, which would be more efficient.

I suppose this raises the idea of a logical data type in Stan

We've already committed to the C-style coding of 0 for false and everything else being true, so if we did go down the boolean route, we'd at least need the type to be promotable to int.

[andrjohns]

As a bit of background for this:

• Vectorised any/all/select functions in progress
• Discussion on bool data type in stanc3