Question about optimisation of interpolated potentials

[smandhai]

Hi again!

I'm having an issue of sorts where the orbits I run using an interpolated potential are running extremely slow (at least for my purposes). 10 orbits within the following potential for example take around 5 minutes. I was wondering if there are any tips or optimisations that I can make/if I'm doing this correctly? I'm currently using a leapfrog integrator to process the orbits. Further context below. Currently using galpy==1.8.0 (1.8.1 has an issue with loading in the unpickled potentials, but this seemingly goes away with 1.8.2 - which I can't install on my current machine as my local pip installation can't find it. This is probably a local issue with my laptop as the issue goes away when I use a different machine.)

spi= pot.InterpSnapshotRZPotential(data_cut.s,rgrid=(np.log10(1e-4),np.log10(1e3),101),logR=True,zgrid=(-60.,60.,101),interpepifreq=True, interpverticalfreq=True,interpPot=True,zsym=False,enable_c=True,numcores=4)

Is the convention with the log grid correct here? I suspect as custom interpolated potentials are generally more detailed than fully analytical potentials, they require longer times to integrate orbits over - is this the case?

#Context:

I am aiming to distil a galaxy simulation (in this case AURIGA - built on GADGET) into a workable potential to calculate orbits via galpy. I have loaded in the galaxy box into pynbody, pre-processed it (cutting around the region of interest), and used galpy to produce a potential for the stellar disc using the snippet above (where data_cut.s is the trimmed stellar component of the trimmed data set). This works fine and I've pickled it using the method from the discussion in #507. I have since unpickled it and loaded it into a separate script for the orbits. The orbits themselves aren't anything fancy, as I just wanted to get the basic foundation functional. I compared the included MWPotential2014 against "spi" potential produced above. The former runs in <1s, whilst the latter can take up to 20-30s per orbit - with the same input parameters. Producing the rotation curve alone can also take some time.

For the DM and Gas components of the galaxy, I originally used a SCF profile for their potentials. However, when compared to the rotation curves produced directly in pynbody, they were off - which I guess was expected as the SCF reflected a Hernquist potential whilst the galaxy's DM for example was more akin to a NFW profile. I have since switched these two similar interpolations as the stellar case above. Stacking the potentials doesn't seem to change the speed of the orbits. The rotation curves are more akin to what pynbody reports too.

For reference: Each component has few x 10-100 x 10^6 particles that are being passed through the interpolation script.

[jobovy]

Hi!

I think I'll need some more info to be able to help.

• How exactly are you integrating the orbits (i.e., what method= are you using? If you are using method='leapfrog' then that will be very slow, even leapfrog_c might be very slow; I don't think there is any reason not to use one of the fast and accurate integrators like dop853_c)
• Could you try setting up an interpRZPotential for MWPotentia2014 using the same grids and then doing the orbit integration with that? Once the potential is interpolated, they should take equally long no matter what the original potential was, so that might help me to locally check the issue.
• I don't think this matters here, but if you're just doing orbit integration, you don't need to set interpepifreq=True, interpverticalfreq=True,interpPot=True.
• Weird about the pickling not working in 1.8.1, but if it works again in 1.8.2 then I guess there is nothing to fix! Maybe the reason you can't upgrade is that your Python version is < 3.8? galpy 1.8.2 requires Python >= 3.8.

could you try just setting up an interpRZPotential with MWPotential2014

[smandhai]

The combined interpolated potential works faster - it takes 10 orbits just under two minutes to finish. Using timesteps with a step size of 100, it takes around 2 seconds. I'll continue to play around but this is somewhat workable for now, as I won't be orbiting a large number of systems just yet.

As I've reduced the gridsize, I'll work in the partial python orbit integration that you suggested! I'll need to figure out the most optimal way to do this, especially if the systems I'm working with leave and re-enter the grid several times.

[jobovy]

That pickling error is very strange, according to the code I don't see how _mass would not be available (6d866af). If you can create a simple reproducible example of the bug, feel free to open an issue. Of course, if you can't reliably reproduce it yourself, then it might be an issue with how you are using the code (perhaps you are using a pickle created after removing the _s?).

If z-symmetry is important to you, you might want to open an issue about that. Should be an easy enough fix, although we'll need to come up with a good test (could add a small z offset bulge so there's a small z-asymmetry that would then be ignored by the C code).

In terms of leaving/re-entering the grid multiple times: I would just integrate any orbit that ever leaves the grid fully in Python, so not worry about resolving the multiple re-enters. Or is that too slow?

[smandhai]

Thank you for the response, Jo! For the pickling issue, I think it's a local issue on my laptop as the exact same script with the exact same pickled objects run fine on the work machines I'm using. I'll tinker with this and find a solution, I'll let you know if there's a fix in case someone else stumbles into the same issue.

I haven't tested the Python integration but I suspect that might be fine to be honest, especially as I assume most objects will take a while to re-enter the grid once they exit.

I'll revisit the z-symmetry issue at a later point as I'm not at a point where I can determine if its impact is huge. I'll think of a test nonetheless and compare the offsets to the results I get with it as it currently is.

For now, I think the adaptive timesteps are working well and the combined potentials (so I think as far as the optimisation in the original question go, this has been resolved!).

I'll continue working on the different parts and add to this thread when I've tested the ideas.

[smandhai]

Related to the _mass error, some additional notes:

Clean anaconda reinstall didn't solve the issue. I realised that the pickled potential was created in 1.8.0, and therefore needed to have _s removed before pickling to get around the Threadlock issue in #509. I tried to create the pickles for the potentials in 1.8.2 but this returns a AttributeError: 'NoneType' object has no attribute 'samples_to_coefficients' (probably owing to the fact that _s and _origPot are removed. I'll try adding these back in and see what happens).

It's still baffling why the scripts and pickles run perfectly fine on a venv environment with galpy 1.8.2 installed on the work machine. The only other explanation I can think of is that the work machine is defaulting to galpy 1.8.0 and ignoring the venv installation (I've checked both via pip to confirm their versions).

Edit: Suspicions confirmed - I loaded ipython within the virtual environment and the galpy version reported (galpy.__version__) is 1.7.2 (version on the base installation of python on my work machine) despite pip reporting 1.8.2. Not sure why this is but it would explain why there's a discontinuity between different machines. I suspect my laptop had a true 1.8.2 version! So I suspect the issue in the end was an incompatibility with pickling pre-1.8.2 snapshots and loading them in via 1.8.2. I'll do the final test where I add _s and _origPot back in and see if I can get everything working on 1.8.2.

Edit 2: Indeed, the _mass error does occur on the second machine when I force the usage of the usage of version 1.8.2. So that reaffirms the incompatibility issue!

[jobovy]

Okay, thanks for checking! So it seems like there is no intrinsic issue with the pickling then 👍