derivative of cel

[rodluger]

@ericagol Is it a simple function of the elliptic integrals, like the derivatives of E and K? If so, we should code it up to avoid numerical instabilities in autodiff.

[ericagol]

Here is the definition from Bulirsch (1969):

[ericagol]

For the derivative with respect to k_c:

Not sure what happens when k_c^2 = p...

\frac{\partial \mathrm{cel}(k_c,p,a,b)}{\partial k_c} = \frac{-k_c}{p-k_c^2}\left[\mathrm{cel}(k_c,k_c^2,a,b)-\mathrm{cel}(k_c,p,a,b)\right]

[ericagol]

The derivatives with respect to a and b are straightforward:

[rodluger]

Nice, thanks!

On Thu, Jun 7, 2018 at 11:57 AM Eric Agol notifications@github.com wrote:

The derivatives with respect to a and b are straightforward:

[image: image] https://user-images.githubusercontent.com/243664/41120273-e26ee3aa-6a49-11e8-8977-5f43f9c68bd6.jpeg

[image: image] https://user-images.githubusercontent.com/243664/41120284-ead4bc7c-6a49-11e8-8b09-6f1601c9054f.jpeg

— You are receiving this because you authored the thread. Reply to this email directly, view it on GitHub https://github.com/rodluger/limbdark/issues/11#issuecomment-395528454, or mute the thread https://github.com/notifications/unsubscribe-auth/AI5FK7FDNbS2pTBfvKha976Ba8NKuDX2ks5t6XeMgaJpZM4Ue6bq .

[ericagol]

@rodluger The p derivative is trickier...

[ericagol]

The p derivative is not straightforward... may be easiest to handle with autodiff.

[rodluger]

OK, that should work. Thanks!

[ericagol]

@rodluger I figured out the derivative of cel with respect to p:

\frac{\partial \mathrm{cel}(k_c,p,a,b)}{\partial p} &= \frac{2pb-pa-b}{2p(1-p)}\mathrm{cel}(k_c,p,0,1)\\ &+\frac{b-ap}{2p(1-p)(p-k_c^2)}\mathrm{cel}(k_c,1,1,k_c^2)\\ &- \frac{b-ap}{2(1-p)(p-k_c^2)}\mathrm{cel}(k_c,p,1,1)

[ericagol]

I've added some julia code, cel_gradient.jl which computes the gradient analytically, with autodiff, and with numerical derivatives of BigFloat, and all agree. I just corrected the expression above for d cel/ d k_c which was missing a factor of -1/2.

[ericagol]

@rodluger Here's a version with only two cel evaluations, and which I think ought to be more stable due to only evaluating the difference \Pi-K, not \Pi alone.

\frac{\partial \mathrm{cel}(k_c,p,a,b)}{\partial p} &= \frac{\mathrm{cel}(k_c,p,0,\lambda) +(b-ap)\mathrm{cel}(k_c,1,1-p,k_c^2-p)}{2p(1-p)(p-k_c^2)},\\ \lambda & = k_c^2(b+ap-2bp)+p(3bp-ap^2-2b).

[rodluger]

@ericagol This worked beautifully. One note, though. When p = 1 or p = 0, the expression for dcel/dp diverges because the denominator blows up. This is actually a case we encounter all the time, since the expression for Lambda1 in the k^2 < 1 case calls cel with p = 1. In practice, it doesn't matter, because in that case p is a constant and the derivative should just be zero, so we can skip the calculation. As you mentioned above, though, there might be a singularity when p = kc^2.

[ericagol]

Yes, p=1 yields a complete elliptic integral of the first or second kind (or a combination thereof); I originally had a separate routine for evaluating this case, but it seemed sensible to reduce the number of routines to call. We can write a wrapper routine that just computes the derivatives with respect to k_c, a and b in these cases. I don't think we ever call p=0. For p=k_c^2, this occurs for k^2 < 1 when b-r=1 which corresponds to first/last contact, at which point the star/planet is unocculted, so all derivatives should be zero. For the k^2 >1 case, k_c^2=p when b=0 or r=0. The r=0 case should have all derivatives of zero, so that leaves the b=0 case, which is a special case, and so the derivatives need to to be handled carefully.... which I still haven't figured out yet.

[rodluger]

It is possible that autodiffing the iterative solution for cel is enough in these edge cases. It could also be more stable than the analytic expression in the vicinity of these singular points.

On Fri, Jun 15, 2018 at 10:32 AM Eric Agol notifications@github.com wrote:

Yes, p=1 yields a complete elliptic integral of the first or second kind (or a combination thereof); I originally had a separate routine for evaluating this case, but it seemed sensible to reduce the number of routines to call. We can write a wrapper routine that just computes the derivatives with respect to k_c, a and b in these cases. I don't think we ever call p=0. For p=k_c^2, this occurs for k^2 < 1 when b-r=1 which corresponds to first/last contact, at which point the star/planet is unocculted, so all derivatives should be zero. For the k^2 >1 case, k_c^2=p when b=0 or r=0. The r=0 case should have all derivatives of zero, so that leaves the b=0 case, which is a special case, and so the derivatives need to to be handled carefully.... which I still haven't figured out yet.

— You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub https://github.com/rodluger/limbdark/issues/11#issuecomment-397691099, or mute the thread https://github.com/notifications/unsubscribe-auth/AI5FKyDGIUsXJU6sp9QdOfq-BZ51LCj0ks5t8--1gaJpZM4Ue6bq .