Custom c++ function may break when called with int parameters

[martinmodrak]

Summary:

Custom C++ functions may return slightly different results when called with seemingly equivalent inputs from Stan. The problem lies in not promoting int to real.

Description:

When using the standard templated header for a user-defined C++ function within a Stan model, it is possible that some of the types resolve to int, despite the parameters being declared as real, resulting in unexpected rounding of values. This may propagate to functions that are called and subtly change behavior at unexpected places.

When the function has only one parameter, passing an int results in a compile-time error (similarly to #2484)

Reproducible Steps:

This code for RStan shows the issue:

library(rstan)

hpp = "
template <typename T0__, typename T1__>
typename boost::math::tools::promote_args<T0__, T1__>::type
test_conversion(const T0__& v, const T1__&z, std::ostream *msgs) {
  T0__ q = log(5);
  return(q);
}
"

stan_code = paste0("
functions {
  real test_conversion(real x, real y);
}


generated quantities {
  print(test_conversion(1, 5));
  print(test_conversion(1.0, 5.0));
}
")

compiled_as_model <- stan_model(model_code = stan_code, allow_undefined = TRUE, includes = hpp)
sampling(compiled_as_model, chains = 1, iter = 1, algorithm = "Fixed_param")

Current Output:

SAMPLING FOR MODEL '0e50413da2d3f8801d6a2a2bbebd1475' NOW (CHAIN 1).
Chain 1: Iteration: 1 / 1 [100%]  (Sampling)
Chain 1: 1
1.60944

Expected Output:

Both results equal, e.g.

Chain 1: Iteration: 1 / 1 [100%]  (Sampling)
Chain 1: 1.60944
1.60944

Current Version:

v2.18.2

I am on Windows & g++

g++ (x86_64-posix-seh, Built by MinGW-W64 project) 4.9.3

[bob-carpenter]

Thanks for reporting. Here's a simpler example that doesn't rely on an external function.

functions {

  real foo(real x, real y) {
    return x / y;
  }

}
transformed data {
  real x = foo(1, 2);
  print("x = ", x);
}

and then running

stan("foo.stan", algorithm = "Fixed_param", chains = 1, iter = 2)

leads to

x = 0

We need to automatically cast int to double when the argument is declared real.

[bob-carpenter]

We can do this by wrapping all the real arguments in a to_real() function, which should cast up integers and get compiled away everywhere else. Converting to pass by value is no big deal here as they're integers to begin with. The first argument could even be (int x), but I think it's clearer as written below.

inline double to_real(const int& x) {
  return x;
}

inline const T& to_real(const T& x) {
  return x;
}