Legendre polynomials

Interfaces RStan
1

Hello everyone,
I’m having trouble using the boost::math::legendre_p function, despite the fact that boost::math::cyl_bessel_k as in the C++ vignette works fine in my code. I’m replacing

double besselK(const double& v, const double& z, std::ostream* pstream__) {
return boost::math::cyl_bessel_k(v, z);
}

in the hpp file with

double legendreP(const int& v, const double& z, std::ostream* pstream__) {
return boost::math::legendre_p(v, z);
}

and adding

real legendreP(int n, real v)

in the functions block. A call to stan_model() with the appropriate include, results in the following error,

Error in dyn.load(libLFile) :
unable to load shared object ‘/var/folders/h_/mp3xxlhd1w3dq4vjwk1f8kcm0000gn/T//RtmpDZF35C/file4816dc396ec.so’:
dlopen(/var/folders/h_/mp3xxlhd1w3dq4vjwk1f8kcm0000gn/T//RtmpDZF35C/file4816dc396ec.so, 6): Symbol not found: ZN5boost4math10legendre_pIdEENS0_5tools12promote_argsIT_fffffE4typeEiS4
Referenced from: /var/folders/h/mp3xxlhd1w3dq4vjwk1f8kcm0000gn/T//RtmpDZF35C/file4816dc396ec.so
Expected in: flat namespace
in /var/folders/h_/mp3xxlhd1w3dq4vjwk1f8kcm0000gn/T//RtmpDZF35C/file4816dc396ec.so

Any insights would be greatly appreciated!

2

My guess is that you are calling it with v as a parameter, in which case you are going to have to do the templating. If you just generate the C++ code, it should tell you what the signature should be but I think it would be something like

template<typename T>
auto legendreP(const int& v, const T&z, std::ostream* pstream__) -> decltype(v + z) {
  return boost::math::legendre_p(v, z);
}
4

No such luck, even tried the Airy function,

template <typename T1__>
typename boost::math::tools::promote_args<T1__>::type
airy(const T1__& z, std::ostream* pstream__) {
return boost::math::airy_ai(z);
}

and got back the same error type. And yet the cyl_bessel_k works fine; very odd.

5

Did you #include the appropriate headers from Boost?

6

Hi Ben,
After trying to add an include within stan_model() without any breakthrough, I decided to modify the stan/model/model_header.hpp file, adding the line

#include <boost/math/special_functions/legendre.hpp>

below the other boost includes, and it worked fine, thanks for your time!

7

Hi Ben,
As a follow up to this discussion, in the preparation of the C++ vignette, was the cyl_bessel_k function ever tested with v or z as parameters? I’ve been trying to use the modified bessel function within a transformed parameter block, and as long as the arguments are fixed, the code compiles, but once I make any one of the arguments a parameter, I start getting a string of 8 errors, starting with

/Library/Frameworks/R.framework/Versions/3.3/Resources/library/BH/include/boost/math/special_functions/trunc.hpp:65:11: error: cannot convert ‘result_type’ (aka ‘stan::math::var’) to ‘int’ without a conversion operator
return static_cast<int>(r);
^~~~~~~~~~~~~~~~~~~
/Library/Frameworks/R.framework/Versions/3.3/Resources/library/BH/include/boost/math/special_functions/trunc.hpp:70:11: note: in instantiation of function template specialization ‘boost::math::itrunc<stan::math::var, boost::math::policies::policy<boost::math::policies::default_policy, boost::math::policies::default_policy, boost::math::policies::default_policy, boost::math::policies::default_policy, boost::math::policies::default_policy, boost::math::policies::default_policy, boost::math::policies::default_policy, boost::math::policies::default_policy, boost::math::policies::default_policy, boost::math::policies::default_policy, boost::math::policies::default_policy> >’ requested here
return itrunc(v, policies::policy<>());
^
/Library/Frameworks/R.framework/Versions/3.3/Resources/library/BH/include/boost/math/special_functions/bessel.hpp:256:24: note: in instantiation of function template specialization ‘boost::math::itrunc<stan::math::var>’ requested here
return bessel_kn(itrunc(v), x, pol);
^
/Library/Frameworks/R.framework/Versions/3.3/Resources/library/BH/include/boost/math/special_functions/bessel.hpp:583:73: note: in instantiation of function template specialization ‘boost::math::detail::cyl_bessel_k_imp<stan::math::var, boost::math::policies::policy<detail::forwarding_arg1, detail::forwarding_arg2, boost::math::policies::default_policy, boost::math::policies::default_policy, boost::math::policies::default_policy, boost::math::policies::default_policy, boost::math::policies::default_policy, boost::math::policies::default_policy, boost::math::policies::default_policy, boost::math::policies::default_policy, boost::math::policies::default_policy> >’ requested here
return policies::checked_narrowing_cast<result_type, Policy>(detail::cyl_bessel_k_imp<value_type>(v, static_cast<value_type>(x), tag_type(), forwarding_policy()), “boost::math::cyl_bessel_k<%1%>(%1%,%1%)”);
^
/Library/Frameworks/R.framework/Versions/3.3/Resources/library/BH/include/boost/math/special_functions/bessel.hpp:589:11: note: in instantiation of function template specialization ‘boost::math::cyl_bessel_k<double, stan::math::var, boost::math::policies::policy<boost::math::policies::default_policy, boost::math::policies::default_policy, boost::math::policies::default_policy, boost::math::policies::default_policy, boost::math::policies::default_policy, boost::math::policies::default_policy, boost::math::policies::default_policy, boost::math::policies::default_policy, boost::math::policies::default_policy, boost::math::policies::default_policy, boost::math::policies::default_policy> >’ requested here
return cyl_bessel_k(v, x, policies::policy<>());

thanks.

8

No. I changed the example for the in CRAN purgatory 2.18 release. You can view it on GitHub

9

I see, thanks. So I assume that this new signature, that’s not based on boost promotion type, runs only on v. 2.18. I just installed v. 2.17.4 and stanc still spits out the boost signature.

10

If I understand you correctly, that is what is supposed to happen. The stanc program generates a C++ file from the Stan file using a Boost metaprogram to calculate the return type of your function (in most cases) and then you need to write C++ implementation of that function for whatever that Boost metaprogram evaluates to.

11

I was referring to the fact that sinc.hpp in the new vignette doesn’t compile under 2.17.4. Only after changing the signature to the boost-related signature does it do so. And even then, when attempting to compile a model that assumes that the argument of the sinc function is unknown – thus requiring the second definition – it again fails to compile.

12

Hi @bgoodri, I installed rstan 2.18.0 hoping that that would take care of the sinc.hpp issues, but they still seem to persist. For example, a model such as

sinc_sim ←
’
functions { real sinc(real x); }
parameters{ real y;}
model {real x = sinc(3.0);
y~normal(x,1);}
’

compiles, but this one doesn’t

sinc_fit ←
’
functions { real sinc(real x); }
data{ int N;
real y[N];}
parameters{ real<lower=0> x;}
model {
x~normal(2,1);
y~normal(sinc(x),1);}
’

thanks.

13

@Bob_Carpenter One of the examples at

#include "my_foo.hpp"
#include <stan/math/rev.hpp>

namespace stan {
  namespace math {
    double foo(double x) {
      return bar::foo(x);
    }

    var foo(const var& x) {
      double a = x.val();
      double fa = bar::foo(x_d);
      double dfa_da = bar::d_foo(a);
      return precomp_v_vari(fa, x.vi_, dfa_da);
    }
  }
}

appears not to be viable because precomp_v_vari returns a vari and it can’t convert a vari to a var. What is it supposed to say?

1 Like
14

Indeed, the output of stan_model when attempting compilation of the model sinc_fit contains the line

[...]/StanHeaders/include/stan/math/prim/mat/fun/assign.hpp:48:5: error: no match for ‘operator=’ (operand types are ‘stan::math::var’ and ‘const stan::math::vari’)

15

@Bob_Carpenter What is supposed to happen with precomp_v_vari? I think the vignette is lying to people like @Phil.

16

Apparently sinc.hpp should have said

double
sinc(const double& x, std::ostream* pstream__) {
  return x != 0.0 ? sin(x) / x : 1.0;
}

stan::math::var
sinc(const stan::math::var& x, std::ostream* pstream__) {
  double x_ = x.val();
  double f = x_ != 0.0 ? sin(x_) / x_ : 1.0;
  double dfdx_ = x_ != 0.0 ? (cos(x_) - sin(x_)) / x_ : 0.0;
  return stan::math::precomputed_gradients(f, {x.vi_}, {dfdx_});
}

which calls precomputed_gradients at the end rather than precomputed_gradients_vari as suggested at

17

Thanks @bgoodri, it’s working perfectly now!

18

It’s super simple, here’s the source:

class precomp_v_vari : public op_v_vari {
 protected:
  double da_;

 public:
  precomp_v_vari(double val, vari* avi, double da)
      : op_v_vari(val, avi), da_(da) {}
  void chain() { avi_->adj_ += adj_ * da_; }
};

It stores one operand pointer (avi) and one partial (da) and applies it as you’d expect. The op_v_vari superclass just stores the value and the operand pointer and makes the operand pointer available as avi_.

19

We figured it out, but the one that returns a vari is not what @Phil needed in order to declare but not define a function in functions and the implement it in C++. He and I needed to be using precomputed_gradients().

20

This is just a more efficient precomputed gradients when there’s only a single argument.

21

But the generated Stan code was expecting a var and said it couldn’t construct a var from the vari.