legendre_assoc(): overflow

[CoChrists]

Using thor.scatter.sph_hrm_coefficients() with q_magnitudes = (1.0, 2.0, 3.0) and num_coefficients = 32 I get this warning:

/usr/local/lib/python2.7/dist-packages/thor/math2.py:375: RuntimeWarning: overflow encountered in int_scalars

Sounds serious to me...

[tjlane]

Can you do me a favor and track down exactly where this is coming from and/or provide me a test case?

[tjlane]

I'm just going to keep encouraging you to get your hands dirty and start writing/fixing code :).

[CoChrists]

here is my case:

#!/usr/bin/python -Wall

import numpy
import thor
import mdtraj.formats.pdb.pdbfile
import sys
import posix

def main(argv):
    pdbfilename = "4IZQ.pdb"
    traj = mdtraj.formats.pdb.pdbfile.load_pdb(pdbfilename)

    if traj is None:
        print "error reading %s" % pdbfilename
        return posix.EX_IOERR

    numpy.random.seed(10101)

    weights          = None
    q_magnitudes     = (1.0, 2.0, 3.0)
    num_coefficients = 32
    coeffsh_even = thor.scatter.sph_hrm_coefficients(traj, q_magnitudes,
                                                     weights,
                                                     num_coefficients)

    return posix.EX_OK

exit(main(sys.argv))

[CoChrists]

Could you provide a reference for your implementation, please? It might help to identify possible sources of numerical instability... Thanks!

[tjlane]

Mmmm there is no real reference besides numpy's docs:
http://docs.scipy.org/doc/numpy/reference/generated/numpy.polynomial.legendre.Legendre.html

[CoChrists]

@tjlane: I don't understand your answer. I am looking for a reference for the formulae in lines https://github.com/tjlane/thor/blob/master/src/python/math2.py#L343 and following.
If I am right, this is your implementation of Associated Legendre functions. Locally, I have implemented another (recursive) version, which I took from Wolfram (see comment #1 (comment) ). It does not show the numerical instabilities but it is very slow (probably due to memory management problems). So I would like to keep your version, but it has to be stable of course. Could you show me which formulae you use in that implementation, please?

[tjlane]

aha yes, I gotcha now!

I stole it from here:

http://mpmath.googlecode.com/svn/trunk/doc/build/functions/orthogonal.html?highlight=legendre#legenp

There is likely a better way. Note that as Robert mentioned previously, scipy is working on improving their spherical harmonic/legendre implementations this summer, so we should team up with them if possible.

[CoChrists]

Ok. Many thanks! What do you mean by "I stole it"? Did you copy their source code or did you just write down their description at
http://mpmath.googlecode.com/svn/trunk/doc/build/functions/orthogonal.html?highlight=legendre#legenp
There are certain differences in your implementation compared to their description. As I don't understand all of them, maybe it would be more efficient for me to meet you somewhere and you tell me the reasoning behind?

[tjlane]

Ah I just did my own implementation of their math equation. Didn't copy any code -- though that might have been a good idea. TBH I was looking for the fastest thing (in TJ time, not computer time) I could get working!

[CoChrists]

Thank you for your help.

[tjlane]

@CoChrists going through old issues and this one came up. Should we just wait for scipy to improve their associated legendre implementation? Or would it be helpful for you if we tackled this?

[CoChrists]

I have not looked into this for a while, but yes, it would be great if these functions would work reliably. We can start investigating again, by trying to reproduce the issue and then go on from there?