exoplanet-dev / exoplanet

Fast & scalable MCMC for all your exoplanet needs!
https://docs.exoplanet.codes
MIT License
206 stars 52 forks source link

Can KeplerianOrbit functions accept theano variables? #58

Closed hposborn closed 5 years ago

hposborn commented 5 years ago

Ok, I'm slowly understanding the whole Theano/PyMc3 back-end, so this is probably a stupid question, but can the functions in KeplerianOrbit accept theano variables?

For example, I would like to include the relative velocity as a deterministic variable in my transit model. This let's me check the transit duration and compare with my past models (which used vel_R to fit the transit). This should be simple enough as the KeplerianOrbit has "get_relative_velocity" which effectively calculates this if you input the time of transit, t0, as the variable.

However, just including something like this: t0 = pm.Normal("t0", mu=list_t0s, sd=1.0, shape=2, testval=list_t0s) orbit = xo.orbits.KeplerianOrbit(period=period, t0=t0, b=b) vel_R = pm.Deterministic("vel_R", orbit.get_relative_velocity(t0)) fails:

`AttributeError Traceback (most recent call last)

in () 95 #print(orbit._get_true_anomaly(t0)[1].shape) 96 #vel_R = pm.Deterministic("vel_R", ((2*np.pi*sma*c.au.value/P)*(1 + ecc*orbit.cosf)/np.sqrt(1-ecc**2))/(Rs*c.R_sun.value)) ---> 97 vel_R = pm.Deterministic("vel_R", orbit.get_relative_velocity(t0)) 98 print(vel_R) 99 Tdur = pm.Deterministic("Tdur", (2*(1+RpRs)*np.sqrt(1-(b/(1+RpRs))**2))/vel_R) /Users/hosborn/anaconda3/lib/python3.6/site-packages/pymc3/model.py in Deterministic(name, var, model) 1496 """ 1497 model = modelcontext(model) -> 1498 var = var.copy(model.name_for(name)) 1499 model.deterministics.append(var) 1500 model.add_random_variable(var) AttributeError: 'tuple' object has no attribute 'copy'`

In this case, t0 has two subtensors (i.e. there are two planets with two t0s), so maybe there is something I can do with the t0 tensor to allow KeplerianOrbit.get_relative_velocity to process the input.

I have also tried using:orbit._get_true_anomaly(t0) to give me cosf and do the calculation myself (e.g. using Eq 12 in Barnes, 2008 https://arxiv.org/pdf/0708.0243.pdf), but this also fails...

Thanks!

setup:

dfm commented 5 years ago

@hposborn: Theano is a whole can of worms :), but here the issue is that get_relative_velocity returns three elements: the vx, vy, and vy. Each of these should have shape (ntimes, nplanets), I believe. So, I expect that what you actually want is the speed in the plane of the sky which you should be able to compute as:

vx, vy, vz = orbit.get_relative_velocity(t0)
vsky = tt.sqrt(vx**2 + vy**2)

but then there's one more catch. This will be a 2x2 matrix with the speed of each planet at each time so you'll want to take the diagonal of that (tt.diag(vsky) should work, I think) to get the right velocities (in R_sun/day). Try that and see if it works!

hposborn commented 5 years ago

Ah thanks. I thought that might be the case, but I couldn't get orbit.get_relative_velocity(t0)[0] (for example) to produce anything either (but I guess I shouldn't treat Theano variables like tuples). Works great now, thanks!