Problem implementing Zipf (power law) distribution -- external C++


#21

That is supposed to work. If you add LDFLAGS += --no-undefined to ~/.R/Makevars you get a slightly more informative error message (an a linker error rather than an error from dyn.load) like:

file5e8d2fcaf228.o: In function `boost::math::tools::promote_args<double, float, float, float, float, float>::type model5e8d11fab598_discrete_power_law_namespace::power_law<double>(int const&, double const&, std::ostream*)':
file5e8d2fcaf228.cpp:(.text._ZN46model5e8d11fab598_discrete_power_law_namespace9power_lawIdEEN5boost4math5tools12promote_argsIT_fffffE4typeERKiRKS5_PSo[_ZN46model5e8d11fab598_discrete_power_law_namespace9power_lawIdEEN5boost4math5tools12promote_argsIT_fffffE4typeERKiRKS5_PSo]+0x2f): undefined reference to `boost::math::tools::promote_args<double, float, float, float, float, float>::type model5e8d11fab598_discrete_power_law_namespace::zeta<double>(double const&, std::ostream*)'

I don’t understand why it thinks it needs a definition of
boost::math::tools::promote_args<double, float, float, float, float, float>::type blah_blah_namespace::zeta<double>(double const&, std::ostream*)
as opposed to what is defined which is simply

double zeta(const double& s, std::ostream* pstream__) {
  return boost::math::zeta(s);
}

And then it does it again for the var specialization

file5e8d2fcaf228.o: In function `boost::math::tools::promote_args<stan::math::var, float, float, float, float, float>::type model5e8d11fab598_discrete_power_law_namespace::power_law<stan::math::var>(int const&, stan::math::var const&, std::ostream*)':
file5e8d2fcaf228.cpp:(.text._ZN46model5e8d11fab598_discrete_power_law_namespace9power_lawIN4stan4math3varEEEN5boost4math5tools12promote_argsIT_fffffE4typeERKiRKS8_PSo[_ZN46model5e8d11fab598_discrete_power_law_namespace9power_lawIN4stan4math3varEEEN5boost4math5tools12promote_argsIT_fffffE4typeERKiRKS8_PSo]+0x109): undefined reference to `boost::math::tools::promote_args<stan::math::var, float, float, float, float, float>::type model5e8d11fab598_discrete_power_law_namespace::zeta<stan::math::var>(stan::math::var const&, std::ostream*)'

#22

OK. I have the problem surrounded. This works, i.e. just declare real zeta(real alpha); but don’t define a power_law function that calls it:

functions{
    real zeta(real s);
    /*TODO: implement rng for prior and posterior predictive checks*/
}
data{
  int<lower=0> K; // number of unique values
  int values[K];
  int<lower=0> frequencies[K];
}
parameters{
  real <lower=1> alpha;
}
model{
  real constant = log(zeta(alpha));
  for (k in 1:K) {
    target += frequencies[k] * (-alpha * log(values[k]) - constant);
  }
} 

That is actually better (iff alpha is a scalar) because it only evaluates the zeta function once. But it does not answer the question (for @Bob_Carpenter or @seantalts) of why you can declare a zeta function, define it in an external C++ file, but not call it from another Stan function without incurring a linker error. In other words, the original Stan program fails to link, even though it should inline to the same thing

functions{
    real zeta(real s);
    real power_law_lpmf(int x, real alpha){
     return (-alpha * log(x) - log(zeta(alpha)) );
    }
    /*TODO: implement rng for prior and posterior predictive checks*/
}
data{
  int<lower=0> K; // number of unique values  
  int values[K];
  int<lower=0> frequencies[K]; 
}
parameters{
  real <lower=1> alpha;
}
model{
  for (k in 1:K) target += frequencies[k] * power_law_lpmf(values[k] | alpha);
}

with errors like


#23

My goal with tagging you and Bob and Ben wasn’t just to ask for help. For example, Bob was the one who initially suggested the Zipf stuff be added to Stan [even though the issue was posted by Daniel Lee] and I wanted him to see the - shy, admittedly - progress on the maths side of things.
I’m hoping that whatever I do here for my own stuff, which will be using Rstan, can then with a bit of work be added to the math library.

By the way, I tried CmdStan and it wasn’t much more intuitive or gave much clearer [to me] error messages. As it unfolded, I would never guessed what the root cause was.


#24

Thanks so much for all the effort you put into this, Ben. This compiles at least. Still some stuff to iron out, but progress in sight.

Yup, nice catch.


#25

I appreciate your sentiment here - someone may indeed come along and find that your code helps them with the full implementation. Cross-posting a reference to this discussion or your code on the github issue would help with that. The core maintainers are a bit overworked so we’d be unlikely to just pick up the issue, but if you (or anyone else reading this) wanted to learn a bit more about contributing to the Math library we’d be happy to help you turn your prototype into production code.

I didn’t mean to imply that anything about C++ development would be intuitive, just that CmdStan has fewer layers of complexity to deconstruct when debugging. So even if you didn’t understand its output, posting that here would have helped the people you were tagging. I actually didn’t know you could even include external C++ files with RStan until reading this thread and this vignette.

You’re sure the external C++ file was included properly and had all the proper overloads defined? Putting up the full model file and include file as a gist on github with the compiler invocations might shed some light on it.


#26

I should have added that “with a bit of work” meant “with a bit of work from me”. I do intend to contribute this to math. But I needed to solve my problem first in order to get on with my research, and that was writing a .hpp that I could call from RStan.
I’ve read the materials you linked and they’re helpful. Due to my C++ illiteracy and inherent thickness, it’ll probably take a while until I can produce anything close to semi-production code for you guys to review.
My main doubt, which I guess applies to both yourself and @Bob_Carpenter is: would it be OK to compute the gradient using the integral representation (and hence an exp_sinh quadrature call)? My fear is that this might be rather too slow and/or inaccurate for the math library.

Let me know if this helps. I’ll edit the gist to add any info you might judge helpful.


#27

I think Ben or another RStan person would know more here - what you’ve posted doesn’t seem to show zeta.hpp getting included in the compiler invocations or the generated hpp. There must be more happening during the RStan invocation to link those up.


#28

Well, there is a linker error if you try to call the C++ function from within another Stan function, so I must have overlooked something but I don’t know what that is.

Here is the C++ with var and a double specialization of the declared but not defined zeta function in a Stan program.

stan::math::var zeta(stan::math::var const& s, std::ostream* pstream__) {
  double s_ = value_of(s);
  double f = boost::math::zeta(s_);
  double d = 0;
  double k = 2;
  double term = 1;
  // http://functions.wolfram.com/ZetaFunctionsandPolylogarithms/Zeta/20/01/01/0001/MainEq1.L.gif
  while (term > 1e-16) {
    term = log(k) / pow(k, s_);
    d -= term;
    k++;
  }
  return var(new precomp_v_vari(f, s.vi_, d));
}

double zeta(const double& s, std::ostream* pstream__) {
  return boost::math::zeta(s);
}

This Stan program works

functions{
    real zeta(real s);
    /*TODO: implement rng for prior and posterior predictive checks*/
}
data{
  int<lower=0> K; // number of unique values  
  int values[K];
  int<lower=0> frequencies[K]; 
}
parameters{
  real <lower=1> alpha;
}
model{
  real constant = log(zeta(alpha));
  for (k in 1:K) {
    target += frequencies[k] * (-alpha * log(values[k]) - constant);
  }
}

with the usual

stan_model(file = "discrete_power_law.stan",  allow_undefined = TRUE,
           includes = paste0('\n#include "', file.path(getwd(), 'zeta.hpp'), '"\n')

but the original Stan model has a linker error

functions{
    real zeta(real s);
    real power_law_lpmf(int x, real alpha){
     return (-alpha * log(x) - log(zeta(alpha)) );
    }
    /*TODO: implement rng for prior and posterior predictive checks*/
}
data{
  int<lower=0> K; // number of unique values  
  int values[K];
  int<lower=0> frequencies[K]; 
}
parameters{
  real <lower=1> alpha;
}
model{
  for (k in 1:K) target += frequencies[k] * power_law_lpmf(values[k] | alpha);
}

So, the only difference is that in the working case, the zeta function is called in the model block directly, whereas in the non-working case, the model block calls power_law_lpmf which calls zeta. The compiler should inline the power_law_lpmf function anyway, in which case, I have no idea why one links correctly and the other has undefined symbols.