Spherical harmonics functions

I’m looking at using spherical harmonic coefficients in a Stan model. For now, they wouldn’t depend on parameters, so I could precalculate them and load them via data but it’d be nice to compute them in a transformed data section. I noticed Boost has an implementation with real and imaginary parts separated. Would this be a useful thing to contribute to Stan?

2 Likes

This is on my list as part of plan to support PDE solution. It’d be great if someone is also ineterested and can contribute, so we can more or less cross special-domain Laplace equation from the list.

1 Like

We’re interested in SHT of neural field equations involving delays. I’m not an expert in PDEs but happy to help.

Looking into SH definitions, it seems prudent to implement in Stan instead of reusing Boost as we’d get AD for free.

What else is on your list? Sparse linear solvers? I’ve noticed PDEs are a numerical analysis rabbit hole so I’m wary of going too deep.

You’re absolutely right. It’s not realistic to reinvent the wheel by putting a PDE solver inside Stan. So my plan is to link a library with adjoint equation capability. In general this would be a mostly a FEM solver. But for special treatment like Laplace equation and wave-like equation we can certainly use SH or any other solutions for better accuracy.

did we get SH yet? I would like to start a prototype with it. Otherwise do you think Boost is the best place to start?

No I haven’t got time working on this, or adjoint solver, and I agree boost can be a good start for this.

OK. I saw elsewhere the complex number support is being worked on in Stan, I might wait a bit for that since it would simplify handling output of these functions.

yeah. I’m also waiting for that PR :)

The complex number thing is hung up on Eigen traits. See this Discourse topic.

It works now other than supporting multiplication of dynamically sized mixed type matrices, like Matrix<complex<double>, -1, -1> by Matrix<var, -1, -1> or Matrix<complex<var>, -1, -1>.

Thanks for the info! I was already taking a shot at a userfunc returning just real or imaginary parts, but it doesn’t find the Boost header for this. Does Stan ship with all of boost::math or just a subset?

it should be all of boost.

The boost lib is pruned of the following files:
doc/, libs/*/test/, libs/*/example/. libs/*/doc/`

and some other files in the main folder: boost.css, boost.png, index.htm, index.html, INSTALL, rst.css

see https://github.com/stan-dev/math/blob/develop/lib/upgrade-boost.sh#L114 and the line below.

OK indeed I found it

prior@posterior:~/Downloads/cmdstan-2.21.0$ find . -name '*spher*harm*.hpp'
./stan/lib/stan_math/lib/boost_1.69.0/boost/math/special_functions/spherical_harmonic.hpp

but compiler can’t find the header (here my model with user defined function is in sh.stan and I’ve generated sh.hpp already)

prior@posterior:~/Downloads/cmdstan-2.21.0$ make sh

--- Compiling, linking C++ code ---
g++ -std=c++1y -pthread -D_REENTRANT -Wno-sign-compare     -I stan/lib/stan_math/lib/tbb_2019_U8/include -O3 -I src -I stan/src -I lib/rapidjson_1.1.0/ -I stan/lib/stan_math/ -I stan/lib/stan_math/lib/eigen_3.3.3 -I stan/lib/stan_math/lib/boost_1.69.0 -I stan/lib/stan_math/lib/sundials_4.1.0/include    -DBOOST_DISABLE_ASSERTS      -c  -x c++ -o sh.o sh.hpp
sh.hpp:25:10: fatal error: boost/math/special_functions/spheric_harmonic.hpp: No such file or directory
 #include "boost/math/special_functions/spheric_harmonic.hpp"
          ^~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
compilation terminated.
make/program:37: recipe for target 'sh' failed
make: *** [sh] Error 1

The folder with Boost is included in the flags to g++ and then I used the relative path as indicated in the error… am I missing something rly basic? (sorry in advance)

Just replacing " with < should do it I think.

We use some of the files in /special_functions, see example here https://github.com/stan-dev/math/blob/develop/stan/math/prim/scal/fun/gamma_p.hpp

Oops, 'twas a typo in the Boost docs linked above, i.e. the header is spherical_harmonic.hpp i.e. with al, which the find command says above but I can’t read. Now I have a boost metaprogramming error traceback to work with.

I’ve the include working and a C++ function

template <typename T2__, typename T3__>
typename boost::math::tools::promote_args<T2__, T3__>::type
spherical_harmonic_real(const int& n,
                            const int& m,
                            const T2__& theta,
                            const T3__& phi, std::ostream* pstream__)
{
        return boost::math::spherical_harmonic_r(n, m, theta, phi);
}

where Boost’s function has the signature

template <class T1, class T2>
calculated-result-type spherical_harmonic_r(unsigned n, int m, T1 theta, T2 phi);

but this ends in tragedy with boost trying to truncate a stan::math::var to an int in tgamma_delta_ratio

template <class T, class Policy>
T tgamma_delta_ratio_imp(T z, T delta, const Policy& pol)
{
...
  return unchecked_factorial<T>((unsigned)itrunc(z, pol) - 1) / unchecked_factorial<T>((unsigned)itrunc(T(z + delta), pol) - 1);

but wait,

         // Both z and delta are integers, see if we can just use table lookup
         // of the factorials to get the result:
         //
         if((z <= max_factorial<T>::value) && (z + delta <= max_factorial<T>::value))

It seems that when Boost chooses a table lookup for what appear to be integers, the cast to int fails for stan::math::var type.

I wonder if there is a Boost policy to tell it not to try this? Is it possible to wrap a var before passing it to Boost to avoid a itrunc failure? Does this mean a hand code derivative is required?

Experience on how to proceed is welcome…

1 Like

It turns out there are known recurrence relations for derivatives of the associated Legendre polynomial of which SH is simple function, so I’ll try the appropriate specializations of this next.

1 Like

Try

return boost::math::spherical_harmonic_r(n, m, value_of_rec(theta), value_of_rec(phi));

which always passes the double wrapped in the var (or var) instead of the var itself. You’ll need to include <stan/math/prim/mat/fun/value_of_rec.hpp>.

Thanks, I will do that for one of the specializations to doubles. From the above link I found the derivatives so will also code that for the Var specialization.

I haven’t found the time to code this up, but the derivatives I found are as follows
\frac{d}{d\phi} Y^m_n(\theta, \phi)=\frac{m-n+1}{1-\cos^2\phi}Y^m+{n+1}(\theta,\phi)-n(\cos \phi)Y^m_n(\theta,\phi)
and
\frac{d}{d\theta} Y^m_n(\theta, \phi)=i m\theta Y^m_n(\theta,\phi)