Swap gsl_gamma_inc with gsl_gamma_inc_P for mass calculation

[jamesgrimmett]

Hi @jobovy, here is the change that we discussed in (i) of #701

For test cases, I've used integrations of the Sun over 1000Gyr (in order to get a decent sample size for profiling) and MW Globular clusters over 10Gyr;

        ts = np.linspace(0,1000,100001)*u.Gyr
        sun = Orbit()
        t = time.time()
        sun.integrate(ts, mw_pot, method="symplec4_c")
        print(f"Duration Sun: {time.time() - t}")

        ts = np.linspace(0,10,1001)*u.Gyr
        o = Orbit.from_name('MW globular clusters')
        t = time.time()
        o.integrate(ts, mw_pot, method="symplec4_c")
        print(f"Duration MW globular clusters: {time.time() - t}")

With the current main branch (using gsl_sf_gamma_inc), I find

    Duration Sun: 9.844529867172241
    Duration MW globular clusters: 6.79642391204834

and with this change (using gsl_sf_gamma_inc_P), I find

    Duration Sun: 7.059226989746094
    Duration MW globular clusters: 1.4422881603240967

So, just less than a 30% improvement for the Sun, and 80% for the collection of MW Globular Clusters, on my machine.

I also took a look at trying out gsl_sf_gamma * (1 - gsl_sf_gamma_inc_Q), but found very similar runtime to with current main;

    Duration Sun: 9.121899127960205
    Duration MW globular clusters: 6.627019882202148

It looks like many of the calls to the gsl P function are diverted to a variation of the Q function (gamma_inc_Q_large_x), but if I call the gsl_sf_gamma_inc_Q directly we mostly end up in gamma_inc_Q_CF => gamma_inc_F_CF (the same place we end up via gsl_sf_gamma_inc), which is interesting... Looking at the gsl code, to filter between gamma_inc_Q_CF and gamma_inc_Q_large_x, the P function considers "large_x" as x => 5*a, whereas the Q function uses x > 1.0e+06. Anyway, it looks like we end up in the most efficient function for typical orbits by calling the P function.
For reference, I've included some of the callmaps obtained by running with gperftools profiler below.

gsl_sf_gamma - gsl_sf_gamma_inc

gsl_sf_gamma * gsl_sf_gamma_inc_P

gsl_sf_gamma * (1 - gsl_sf_gamma_inc_Q)

[codecov]

Codecov Report

All modified and coverable lines are covered by tests ✅

Project coverage is 99.90%. Comparing base (ef6eaf0) to head (0a8c448).
Report is 1 commits behind head on main.

Additional details and impacted files

@@           Coverage Diff           @@
##             main     #706   +/-   ##
=======================================
  Coverage   99.90%   99.90%           
=======================================
  Files         200      200           
  Lines       29433    29433           
  Branches      563      563           
=======================================
  Hits        29406    29406           
  Misses         27       27           

☔ View full report in Codecov by Sentry.
📢 Have feedback on the report? Share it here.

[jobovy]

Thanks @jamesgrimmett, this is fantastic!

I experimented with this a bit myself and get consistent results, a 35% speed-up for the orbit of the Sun and almost 5x speed-up for the globular clusters. To try to figure out where the difference really happens, I also did the MW satellite galaxies and got only an 18% speed-up and then I did an orbit within the bulge region (o = Orbit([0.1,0.1,0.6,0.1,0.1,0.1])) and got an ~11x speed-up! So it seems like the difference is biggest for orbits close to the bulge, which is probably because of the different branching that you mention.

Anyway, everything is generally faster, so I'll merge this. Thanks again for looking into this!