lkreidberg / batman

Fast transit light curves models in Python.
http://lkreidberg.github.io/batman
GNU General Public License v3.0
83 stars 52 forks source link

Slow computation of light curves with nonlinear limb darkening models. #65

Open NTSchragal opened 6 months ago

NTSchragal commented 6 months ago

Hello! I am currently attempting to use BATMAN for light curve fitting using nonlinear limb darkening coefficients. In pursuit of this, I am doing an MCMC (thanks emcee) analysis which, of course, relies greatly on computational efficiency. I have discovered that light curve evaluation with nonlinear limb darkening is significantly (~25x) slower than evaluation with quadratic limb darkening. I think I have found why; my quarrel is with this chunk of c_src/_nonlinear_ld.c:

inline double intensity(double x, double* args)
{
    double c1=args[0], c2=args[1], c3=args[2], c4=args[3], norm=args[4];
    if(x > 0.99995) x = 0.99995;
    double sqrtmu = pow(1. - x*x,0.25);
    return (1. - c1*(1. - sqrtmu) - c2*(1. - pow(sqrtmu,2.)) - c3*(1. - pow(sqrtmu, 3.)) - c4*(1. - pow(sqrtmu,4.)))/norm; 
}

I admit I am no expert, but I believe the repeated calls to pow() are slowing things down.

I don't have the time right now to test this myself (hence opening the issue), but could this be accelerated by, in the 3rd line of the function, calling sqrt() twice instead of pow(__, 0.25) (since sqrt() is a less general function than pow()) and, in the 4th line, explicitly writing out the powers as multiplication. E.g.:

inline double intensity(double x, double* args)
{
    double c1=args[0], c2=args[1], c3=args[2], c4=args[3], norm=args[4];
    if(x > 0.99995) x = 0.99995;
    double sqrtmu = sqrt(sqrt(1. - x*x));
    return (1. - c1*(1. - sqrtmu) - c2*(1. - sqrtmu*sqrtmu) - c3*(1. - sqrtmu*sqrtmu*sqrtmu) - c4*(1. - sqrtmu*sqrtmu*sqrtmu*sqrtmu))/norm; 
}

(Yes, that last bit of the 4th line looks horrible but c'est la vie.)

Thanks in advance for taking a look at this. If I get time I will try to test this solution myself, but for now here is the issue for posterity.

lkreidberg commented 5 months ago

hi there, thanks for the message! I think you're right the POW is slower than SQRT, but what is really slowing thing down is that the nonlinear LD law requires an integral over the obscured part of the star, whereas the quadratic law is an analytic function.

There are some new and improved analytic models from Agol et al. 2020 for more general cases, but I haven't gotten around to implementing these. Leaving this issue open in case anyone is inspired to work on it :)