Coder Social home page Coder Social logo

exoplanet-dev / exoplanet Goto Github PK

View Code? Open in Web Editor NEW
202.0 18.0 54.0 15.13 MB

Fast & scalable MCMC for all your exoplanet needs!

Home Page: https://docs.exoplanet.codes

License: MIT License

Python 80.98% TeX 19.02%
exoplanet exoplanet-transits exoplanet-radial-velocities mcmc gaussian-processes bayesian-inference pymc3 time-series time-series-analysis astronomy

exoplanet's Introduction



Tests Docs Coverage

exoplanet

Fast & scalable MCMC for all your exoplanet needs! exoplanet is a toolkit for probabilistic modeling of time series data in astronomy with a focus on observations of exoplanets, using PyMC. PyMC is a flexible and high-performance model-building language and inference engine that scales well to problems with a large number of parameters. exoplanet extends PyMC's language to support many of the custom functions and distributions required when fitting exoplanet datasets.

Read the full documentation at docs.exoplanet.codes.

Installation

The quickest way to get started is to use pip:

python -m pip install "exoplanet[pymc]"

Note that you will need Python (>=3.6) installed for this to work, but then this will install all the required dependencies.

Check out the main installation documentation page for more options.

Usage

Check out the tutorials and API docs on the docs page for example usage and much more info. You can also find more in-depth examples on the exoplanet case studies page.

Contributing

exoplanet is an open-source project, and we would love it if you wanted to contribute. Check out the developer documentation for more info about getting started.

exoplanet's People

Contributors

adrn avatar arjunsavel avatar avivajpeyi avatar christinahedges avatar dependabot[bot] avatar dfm avatar emilygilbert avatar ericagol avatar gjgilbert avatar iancze avatar lgbouma avatar mrtommyb avatar pre-commit-ci[bot] avatar rodluger avatar t-brandt avatar vandalt avatar

Stargazers

 avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar

Watchers

 avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar

exoplanet's Issues

Wavelet Log Likelihood (Carter & Winn 2009)

Is your feature request related to a problem? Please describe.
I'm using exoplanet to integrate PyMC3 with efficient transit modeling. I have the base system running smoothly -- thanks to great tutorials. Thank you.

I was requested by a collaborator to include the Carter & Winn (2009) Daubechies wavelet likelihood.

(1) I am able to generate log_likelihood as a float value. So that first feature request that I would ask is how to add that into the xo.optimize or pm.sample calls?

(2) If we can bypass wavelets altogether by using a GP to estimate sigma_r and sigma_w efficiently, then I would be more than happy to do that. We've previously discussed using a GP to model the Daubechies wavelet likelihood; but I would need for help in selecting the GP+Kernels combinations to do it.

As a secondary feature request, I would ask for a tutorial to estimate the residual red noise in a light curve the way that our community to put in tables.

Describe the solution you'd like
I see from #70, that adding in a likelihood function is feasible; but I get the impression that there is more to do than a simple xo_logp + wavelet_logp.

Is there a way to take a function (with 2 inputs + 3 PyM3 parameters) and "add" it to the likelihood computed inside exoplanet?

Describe alternatives you've considered
I grabbed the necessary bits of code form @pcubillos's package mc3
https://github.com/pcubillos/mc3/blob/master/mc3/stats/stats.py [line 209]

I can compute the wavelet likelihood (outside of PyMC3+XO) as a float value; but how to wrap this into something that xo can interpret as an addition to the log-likelihood built in?

Additional context
If there is a much simpler way to use the GP module to estimate sigma_r and sigma_w, then I could bypass this whole issue.

At the same time, others may still want a wavelet likelihood added onto their light curve posteriors. As such, may I suggest this feature overall.

Consider whether exoplanet should use DensityDist rather than Potential

PyMC3 recommended way to include a custom likelihood function is with pm.DensityDist, e.g.

loglike = pm.DensityDist('loglike', logp, observed=y)

However in the exoplanet examples, the the GP likelihood is implemented as pm.Potential. The results are the same but it might be sensible to implement the samples using the standard PyMC3 way of doing it. Advantages include clear separating out of the observables. Additionally, some of the nice plotting functions like model_to_graphviz ignore potentials.

Here's an example of what it would look like

pm.DensityDist("loglike", gp.log_likelihood, observed=(y - light_curve - mean))

Support for transit duration variations

A feature request from @gjgilbert.

A more general interface could be to allow any of the orbital elements to be a function of transit number. Right now, transit time is hard coded as the only one to vary, but it shouldn't be too bad to add impact parameter (or inclination? or both? things start to get a little tricky) as another parameter that varies. I'm inclined to think that the more general interface might be easier to maintain in the long run, but I'm open to suggestions!

get_true_anomaly when mean anomaly is pi

I've been using exoplanet quite a bit the past few days and noticed some weird behaviour in exoplanet.orbits.get_true_anomaly. When the mean anomaly is exactly pi, get_true_anomaly returns 0, whereas my implementation (and some others I've tried) return pi. You can see this in action below.

from exoplanet.utils import eval_in_model
from exoplanet.orbits import get_true_anomaly
import pymc3 as pm
import theano.tensor as tt

period = 10.
eccen = 0.5
asini = 100.
varpi = 1.

# TESS time sampling
time = np.arange(0, 27, 1.0 / (24 * 30))

with pm.Model() as model:
    M = tt.zeros_like(tt.constant(time)) + 2.0 * np.pi * (tt.constant(time)) / period
    f = get_true_anomaly(M, eccen + tt.zeros_like(M))
    tau_tens = (- (1 - tt.square(eccen)) * tt.sin(f+varpi) / (1 + eccen*tt.cos(f))) * (asini / 86400.)
    
    tau_val = eval_in_model(tau_tens)
    true_anom_val = eval_in_model(f)
    
mean_anom = np.zeros_like(time) + 2.0 * np.pi * (time) / period
# kepler() is a basic Newton's method solver I implemented
ecc_anom = kepler(mean_anom, eccen + np.zeros_like(mean_anom))
true_anom_python = 2.0 * np.arctan2(np.sqrt(1.0+eccen)*np.tan(0.5*ecc_anom),np.sqrt(1.0-eccen) + np.zeros_like(time))
tau_python = -(asini/86400) * (1.0 - np.square(eccen)) * np.sin(true_anom_python + varpi) / (1.0 + eccen*np.cos(true_anom_python))

residual

This flips the sign on the cos term in the tau variable, and makes the plot look a bit funky:

tau_compare

Write some TESS tutorials

  • Figure out how to fit the TESS alerts
  • write a tutorial with demos with initial TESS discoveries
  • resample a TESS light curve as if it is an FFI and fit that

Help with program

Hi everyone,

My name is Alex and im an undergraduate student at Loughborough university, currently undergoing my final year project in which im using data from TESS to hunt for my own exoplanet candidates.

Im trying to use radial velocities to help calculate masses, in order to negate eclipsing binaries from my shortlist etc but il admit my coding is sub par and im not sure how to go about it, however i was pointed in the direction of your program!

Im using TESS target pixel files and have created lightcurves that include all flux and time metadata, but where do you find values for radial velocities such that i can create a radial velocity plot? Is it something that i calculate using my flux values or do you go to an archive somewhere that has dedicated RV values?

Sorry about the probably face palm stupid question! I appreciate any help that can be given to steer me in the right direction towards calculating orbiting masses.

Many thanks
Alex

Installation error in windows gcc exit status 1 with no sys/mman.h file

I run the following commands to install exoplanet in a new anaconda environment:

conda create -n exoplanet_env python=3.7
conda activate exoplanet_env
python -m pip install --upgrade pip
conda install numpy astropy pymc3
python -m pip install -U exoplanet

And then I get:

(exoplanet_env) C:\Users\Vital>python -m pip install -U exoplanet
Collecting exoplanet
  Using cached https://files.pythonhosted.org/packages/7e/c6/a59b918455922c915e71d040cb0f069d0c6d2653d33f38e12d27c875c6c8/exoplanet-0.2.3-py2.py3-none-any.whl
Requirement already satisfied, skipping upgrade: pymc3>=3.5 in d:\programs\anaconda3\envs\exoplanet_env\lib\site-packages (from exoplanet) (3.7)
Requirement already satisfied, skipping upgrade: theano>=1.0.4 in d:\programs\anaconda3\envs\exoplanet_env\lib\site-packages (from exoplanet) (1.0.4)
Requirement already satisfied, skipping upgrade: astropy>=3.1 in d:\programs\anaconda3\envs\exoplanet_env\lib\site-packages (from exoplanet) (3.2.3)
Requirement already satisfied, skipping upgrade: numpy>=1.13.0 in d:\programs\anaconda3\envs\exoplanet_env\lib\site-packages (from exoplanet) (1.17.4)
Processing c:\users\vital\appdata\local\pip\cache\wheels\62\39\dd\45d51f03ee2b1c044c89fcf3a08cb0496734605d7b9ef4e5b6\rebound_pymc3-0.0.3-cp37-none-any.whl
Requirement already satisfied, skipping upgrade: h5py>=2.7.0 in d:\programs\anaconda3\envs\exoplanet_env\lib\site-packages (from pymc3>=3.5->exoplanet) (2.9.0)
Requirement already satisfied, skipping upgrade: patsy>=0.4.0 in d:\programs\anaconda3\envs\exoplanet_env\lib\site-packages (from pymc3>=3.5->exoplanet) (0.5.1)
Requirement already satisfied, skipping upgrade: tqdm>=4.8.4 in d:\programs\anaconda3\envs\exoplanet_env\lib\site-packages (from pymc3>=3.5->exoplanet) (4.40.2)
Requirement already satisfied, skipping upgrade: scipy>=0.18.1 in d:\programs\anaconda3\envs\exoplanet_env\lib\site-packages (from pymc3>=3.5->exoplanet) (1.3.2)
Requirement already satisfied, skipping upgrade: pandas>=0.18.0 in d:\programs\anaconda3\envs\exoplanet_env\lib\site-packages (from pymc3>=3.5->exoplanet) (0.25.3)
Requirement already satisfied, skipping upgrade: six>=1.9.0 in d:\programs\anaconda3\envs\exoplanet_env\lib\site-packages (from theano>=1.0.4->exoplanet) (1.13.0)
Collecting rebound>=3.10.1
  Using cached https://files.pythonhosted.org/packages/fc/c5/0bab299d215e9da2d93607344957e726ab29b2f13a86c0cf09b637bd48a0/rebound-3.12.1.tar.gz
Requirement already satisfied, skipping upgrade: python-dateutil>=2.6.1 in d:\programs\anaconda3\envs\exoplanet_env\lib\site-packages (from pandas>=0.18.0->pymc3>=3.5->exoplanet) (2.8.1)
Requirement already satisfied, skipping upgrade: pytz>=2017.2 in d:\programs\anaconda3\envs\exoplanet_env\lib\site-packages (from pandas>=0.18.0->pymc3>=3.5->exoplanet) (2019.3)
Building wheels for collected packages: rebound
  Building wheel for rebound (setup.py) ... error
  ERROR: Command errored out with exit status 1:
   command: 'D:\Programs\Anaconda3\envs\exoplanet_env\python.exe' -u -c 'import sys, setuptools, tokenize; sys.argv[0] = '"'"'C:\\Users\\Vital\\AppData\\Local\\Temp\\pip-install-9vkxn_kd\\rebound\\setup.py'"'"'; __file__='"'"'C:\\Users\\Vital\\AppData\\Local\\Temp\\pip-install-9vkxn_kd\\rebound\\setup.py'"'"';f=getattr(tokenize, '"'"'open'"'"', open)(__file__);code=f.read().replace('"'"'\r\n'"'"', '"'"'\n'"'"');f.close();exec(compile(code, __file__, '"'"'exec'"'"'))' bdist_wheel -d 'C:\Users\Vital\AppData\Local\Temp\pip-wheel-u01p09vw' --python-tag cp37
       cwd: C:\Users\Vital\AppData\Local\Temp\pip-install-9vkxn_kd\rebound\
  Complete output (31 lines):
  fatal: not a git repository (or any of the parent directories): .git
  running bdist_wheel
  running build
  running build_py
  creating build
  creating build\lib.win-amd64-3.7
  creating build\lib.win-amd64-3.7\rebound
  copying rebound\data.py -> build\lib.win-amd64-3.7\rebound
  copying rebound\horizons.py -> build\lib.win-amd64-3.7\rebound
  copying rebound\interruptible_pool.py -> build\lib.win-amd64-3.7\rebound
  copying rebound\particle.py -> build\lib.win-amd64-3.7\rebound
  copying rebound\plotting.py -> build\lib.win-amd64-3.7\rebound
  copying rebound\simulation.py -> build\lib.win-amd64-3.7\rebound
  copying rebound\simulationarchive.py -> build\lib.win-amd64-3.7\rebound
  copying rebound\tools.py -> build\lib.win-amd64-3.7\rebound
  copying rebound\units.py -> build\lib.win-amd64-3.7\rebound
  copying rebound\widget.py -> build\lib.win-amd64-3.7\rebound
  copying rebound\__init__.py -> build\lib.win-amd64-3.7\rebound
  copying rebound\rebound.h -> build\lib.win-amd64-3.7\rebound
  running build_ext
  building 'librebound' extension
  creating build\temp.win-amd64-3.7
  creating build\temp.win-amd64-3.7\Release
  creating build\temp.win-amd64-3.7\Release\src
  D:\Programs\Anaconda3\envs\exoplanet_env\Library\mingw-w64\bin\gcc.exe -mdll -O -Wall -DMS_WIN64 -DLIBREBOUND -Isrc -ID:\Programs\Anaconda3\envs\exoplanet_env\include -ID:\Programs\Anaconda3\envs\exoplanet_env\include -c src/rebound.c -o build\temp.win-amd64-3.7\Release\src\rebound.o -fstrict-aliasing -O3 -std=c99 -Wno-unknown-pragmas -DGITHASH=d15889dce6327b47ace8111b0b376361aba81ff9 -DLIBREBOUND -D_GNU_SOURCE -fPIC
  src/rebound.c:1:0: warning: -fPIC ignored for target (all code is position independent)
   /**
   ^
  src/rebound.c:33:22: fatal error: sys/mman.h: No such file or directory
  compilation terminated.
  error: command 'D:\\Programs\\Anaconda3\\envs\\exoplanet_env\\Library\\mingw-w64\\bin\\gcc.exe' failed with exit status 1
  ----------------------------------------
  ERROR: Failed building wheel for rebound
  Running setup.py clean for rebound
Failed to build rebound
Installing collected packages: rebound, rebound-pymc3, exoplanet
    Running setup.py install for rebound ... error
    ERROR: Command errored out with exit status 1:
     command: 'D:\Programs\Anaconda3\envs\exoplanet_env\python.exe' -u -c 'import sys, setuptools, tokenize; sys.argv[0] = '"'"'C:\\Users\\Vital\\AppData\\Local\\Temp\\pip-install-9vkxn_kd\\rebound\\setup.py'"'"'; __file__='"'"'C:\\Users\\Vital\\AppData\\Local\\Temp\\pip-install-9vkxn_kd\\rebound\\setup.py'"'"';f=getattr(tokenize, '"'"'open'"'"', open)(__file__);code=f.read().replace('"'"'\r\n'"'"', '"'"'\n'"'"');f.close();exec(compile(code, __file__, '"'"'exec'"'"'))' install --record 'C:\Users\Vital\AppData\Local\Temp\pip-record-uxwu9jsr\install-record.txt' --single-version-externally-managed --compile
         cwd: C:\Users\Vital\AppData\Local\Temp\pip-install-9vkxn_kd\rebound\
    Complete output (31 lines):
    fatal: not a git repository (or any of the parent directories): .git
    running install
    running build
    running build_py
    creating build
    creating build\lib.win-amd64-3.7
    creating build\lib.win-amd64-3.7\rebound
    copying rebound\data.py -> build\lib.win-amd64-3.7\rebound
    copying rebound\horizons.py -> build\lib.win-amd64-3.7\rebound
    copying rebound\interruptible_pool.py -> build\lib.win-amd64-3.7\rebound
    copying rebound\particle.py -> build\lib.win-amd64-3.7\rebound
    copying rebound\plotting.py -> build\lib.win-amd64-3.7\rebound
    copying rebound\simulation.py -> build\lib.win-amd64-3.7\rebound
    copying rebound\simulationarchive.py -> build\lib.win-amd64-3.7\rebound
    copying rebound\tools.py -> build\lib.win-amd64-3.7\rebound
    copying rebound\units.py -> build\lib.win-amd64-3.7\rebound
    copying rebound\widget.py -> build\lib.win-amd64-3.7\rebound
    copying rebound\__init__.py -> build\lib.win-amd64-3.7\rebound
    copying rebound\rebound.h -> build\lib.win-amd64-3.7\rebound
    running build_ext
    building 'librebound' extension
    creating build\temp.win-amd64-3.7
    creating build\temp.win-amd64-3.7\Release
    creating build\temp.win-amd64-3.7\Release\src
    D:\Programs\Anaconda3\envs\exoplanet_env\Library\mingw-w64\bin\gcc.exe -mdll -O -Wall -DMS_WIN64 -DLIBREBOUND -Isrc -ID:\Programs\Anaconda3\envs\exoplanet_env\include -ID:\Programs\Anaconda3\envs\exoplanet_env\include -c src/rebound.c -o build\temp.win-amd64-3.7\Release\src\rebound.o -fstrict-aliasing -O3 -std=c99 -Wno-unknown-pragmas -DGITHASH=d15889dce6327b47ace8111b0b376361aba81ff9 -DLIBREBOUND -D_GNU_SOURCE -fPIC
    src/rebound.c:1:0: warning: -fPIC ignored for target (all code is position independent)
     /**
     ^
    src/rebound.c:33:22: fatal error: sys/mman.h: No such file or directory
    compilation terminated.
    error: command 'D:\\Programs\\Anaconda3\\envs\\exoplanet_env\\Library\\mingw-w64\\bin\\gcc.exe' failed with exit status 1
    ----------------------------------------
ERROR: Command errored out with exit status 1: 'D:\Programs\Anaconda3\envs\exoplanet_env\python.exe' -u -c 'import sys, setuptools, tokenize; sys.argv[0] = '"'"'C:\\Users\\Vital\\AppData\\Local\\Temp\\pip-install-9vkxn_kd\\rebound\\setup.py'"'"'; __file__='"'"'C:\\Users\\Vital\\AppData\\Local\\Temp\\pip-install-9vkxn_kd\\rebound\\setup.py'"'"';f=getattr(tokenize, '"'"'open'"'"', open)(__file__);code=f.read().replace('"'"'\r\n'"'"', '"'"'\n'"'"');f.close();exec(compile(code, __file__, '"'"'exec'"'"'))' install --record 'C:\Users\Vital\AppData\Local\Temp\pip-record-uxwu9jsr\install-record.txt' --single-version-externally-managed --compile Check the logs for full command output.

I should just install linux in this machine but if anyone has encountered this issue I would love to hear your advice (merry xmas :)

system/environment issue on supercomputer

Hi,
I am doing unit-test of exoplanet on supercomputer cluster. However I met a problem of system or environment issue about it. It seems about the flag of theano module. Here is the test feedback. Do you have any insight how to configure those possible wrong flag value in the theano backend?
Thanks,
Yangyang
bug_infos.log

jitter not used in example 51 Peg fit?

In the example of fitting 51 Peg data with exoplanet, there's a log-jitter variable defined, but I don't see where it comes into the fit.

My guess on how to fix this is that you should change from
pm.Normal("obs", mu=rvmodel, sd=rv_err, observed=rv)
to
pm.Normal("obs", mu=rvmodel, sd=tt.sqrt(tt.sqr(rv_err) + s2), observed=rv).
Is that ok, or does there need to be a penalty term in the likelihood function for large jitter explicitly defined also? Or am I misunderstanding what's happening here entirely??

Variable exposure time broken when `use_in_transit=True`

The following code

import numpy as np
import matplotlib.pyplot as plt
import exoplanet as xo

# Generate a starry light curve
orbit = xo.orbits.KeplerianOrbit(period=1)
t = np.linspace(-0.5, 0.5, 1000)
u = [0.3, 0.2]
starry_op = xo.StarryLightCurve(u)

# Compute the light curve
lc_int = starry_op.get_light_curve(orbit=orbit, r=0.1, t=t, texp=0.1 * np.ones_like(t)).eval()

throws a very long error starting with

ERROR (theano.gof.opt): Optimization failure due to: constant_folding
ERROR (theano.gof.opt): node: Elemwise{add,no_inplace}(TensorConstant{[[-0.09159..09159159]]}, TensorConstant{[[-0.05   ..05      ]]})
ERROR (theano.gof.opt): TRACEBACK:
ERROR (theano.gof.opt): Traceback (most recent call last):
  File "/anaconda3/lib/python3.7/site-packages/theano/gof/opt.py", line 2034, in process_node
    replacements = lopt.transform(node)
  File "/anaconda3/lib/python3.7/site-packages/theano/tensor/opt.py", line 6518, in constant_folding
    required = thunk()
  File "/anaconda3/lib/python3.7/site-packages/theano/gof/op.py", line 862, in rval
    thunk()
  File "/anaconda3/lib/python3.7/site-packages/theano/gof/cc.py", line 1739, in __call__
    reraise(exc_type, exc_value, exc_trace)
  File "/anaconda3/lib/python3.7/site-packages/six.py", line 693, in reraise
    raise value
ValueError: Input dimension mis-match. (input[0].shape[0] = 184, input[1].shape[0] = 1000)

The issue is on this line, which assumes that dt and t have the same shape. They don't because use_in_transit leads to t being cropped.

NotADirectoryError when importing exoplanet

I ran the tutorial as written here:
http://exoplanet.dfm.io/en/latest/tutorials/stellar-variability/

and get the following error. I'm using a git clone of the master branch

NotADirectoryError                        Traceback (most recent call last)
<ipython-input-2-19985ffa6e50> in <module>()
----> 1 import exoplanet as xo
      2 
      3 results = xo.estimators.lomb_scargle_estimator(
      4     x, y, max_peaks=1, min_period=5.0, max_period=100.0,
      5     samples_per_peak=50)

~/anaconda3/lib/python3.6/importlib/_bootstrap.py in _find_and_load(name, import_)

~/anaconda3/lib/python3.6/importlib/_bootstrap.py in _find_and_load_unlocked(name, import_)

~/anaconda3/lib/python3.6/importlib/_bootstrap.py in _load_unlocked(spec)

~/anaconda3/lib/python3.6/importlib/_bootstrap.py in _load_backward_compatible(spec)

~/anaconda3/lib/python3.6/site-packages/exoplanet-0.1.2.dev0-py3.6.egg/exoplanet/__init__.py in <module>()
     14     from .sampling import *  # NOQA
     15     from .estimators import *  # NOQA
---> 16     from .light_curves import *  # NOQA
     17 
     18     from . import distributions, gp, orbits

~/anaconda3/lib/python3.6/importlib/_bootstrap.py in _find_and_load(name, import_)

~/anaconda3/lib/python3.6/importlib/_bootstrap.py in _find_and_load_unlocked(name, import_)

~/anaconda3/lib/python3.6/importlib/_bootstrap.py in _load_unlocked(spec)

~/anaconda3/lib/python3.6/importlib/_bootstrap.py in _load_backward_compatible(spec)

~/anaconda3/lib/python3.6/site-packages/exoplanet-0.1.2.dev0-py3.6.egg/exoplanet/light_curves.py in <module>()
     13 from .theano_ops.starry.limbdark import LimbDarkOp
     14 
---> 15 get_cl = GetClOp()
     16 limbdark = LimbDarkOp()
     17 

~/anaconda3/lib/python3.6/site-packages/exoplanet-0.1.2.dev0-py3.6.egg/exoplanet/theano_ops/starry/get_cl.py in __init__(self)
     20 
     21     def __init__(self):
---> 22         self.grad_op = GetClRevOp()
     23         super(GetClOp, self).__init__()
     24 

~/anaconda3/lib/python3.6/site-packages/exoplanet-0.1.2.dev0-py3.6.egg/exoplanet/theano_ops/starry/base_op.py in __init__(self)
     19 
     20     def __init__(self):
---> 21         super(StarryBaseOp, self).__init__(self.func_file, self.func_name)
     22 
     23     # def c_code_cache_version(self):

~/anaconda3/lib/python3.6/site-packages/theano/gof/op.py in __init__(self, func_files, func_name)
   1306         # find the new path and new version of the file in Theano.
   1307         self.func_files = func_files
-> 1308         self.load_c_code(func_files)
   1309 
   1310         if len(self.code_sections) == 0:

~/anaconda3/lib/python3.6/site-packages/theano/gof/op.py in load_c_code(self, func_files)
   1329         for func_file in func_files:
   1330             # U (universal) will convert all new lines format to \n.
-> 1331             with _open_u(func_file) as f:
   1332                 self.func_codes.append(f.read())
   1333 

NotADirectoryError: [Errno 20] Not a directory: '/Users/tom/anaconda3/lib/python3.6/site-packages/exoplanet-0.1.2.dev0-py3.6.egg/exoplanet/theano_ops/starry/get_cl_rev.cc'

Exposure integration fails for large exposure times

When the exposure time is comparable to the duration of the transit (which can occur for ultra-short period planets), the integration algorithm in exoplanet produces spurious wiggles. The code below compares a light curve integrated with exoplanet (orange line in plot) and one using a simple moving average (green line in plot). This occurs (with varying degrees) for all three integration order values and is not related to the boundary effects: even though I'm plotting the light curve only in the vicinity of transit, I'm computing it on a much wider time grid. Note that setting use_in_transit = False does not fix this issue.

import numpy as np
import matplotlib.pyplot as plt
import exoplanet as xo


def moving_average(a, n):
    """
    Compute a moving average over a window
    of `n` points. Based on
    https://stackoverflow.com/a/14314054
    """
    if n % 2 != 0:
        n += 1
    ret = np.cumsum(a, dtype=float)
    ret[n:] = ret[n:] - ret[:-n]
    result = ret[n - 1:] / n
    return np.concatenate([np.ones(n // 2) * result[0],
                           result,
                           np.ones(n // 2 - 1) * result[-1]])

# Generate a starry light curve
orbit = xo.orbits.KeplerianOrbit(period=1)
t = np.linspace(-0.5, 0.5, 1000)
u = [0.3, 0.2]
starry_op = xo.StarryLightCurve(u)

# Compute the light curve
lc = starry_op.get_light_curve(orbit=orbit, r=0.1, t=t).eval()

# Compute the light curve with finite exposure time
lc_int = starry_op.get_light_curve(orbit=orbit, r=0.1, t=t, texp=0.1, order=2).eval()

# Same thing, but numerically using a moving average
lc_int_num = moving_average(lc, int(0.1 / (t[1] - t[0])))

# Plot them
plt.plot(t, lc, lw=2, label="Flux")
plt.plot(t, lc_int, lw=2, label="Integrated with exoplanet")
plt.plot(t, lc_int_num, lw=2, label="Integrated with moving average")
plt.ylabel("relative flux")
plt.xlabel("time [days]")
plt.xlim(-0.1, 0.1)
plt.margins(None, 0.3)
plt.legend()
plt.show()

image

Hacking exoplanet to fit two disparate transits without "islands" in the period PDF...

For a few transiting systems (and increasingly more with the gaps likely to occur in TESS) we have two individual transits that are separated by a long gap, meaning we have no real knowledge of the orbital period except for a few discrete integer period ratios. These create "islands" in the posteriors and make MCMC modelling extremely difficult... The most well known such case is HIP 41378 (https://arxiv.org/abs/1809.10688).
I'm trying to come up with a (probably hacky) way that exoplanet can model such systems without converging on only a single period island... Any ideas would be appreciated.

  • The least hacky method is to simply force a model to converge on each of the N period "islands" and compute an evidence (or BIC) for each... Which is obviously not simple with an MCMC compared to, say, PyMultiNest.

  • An easier way (and the way I probably have hope of coding) might be to model each transit separately and find the intersection of the posteriors, but I feel like we're losing information (and certainly efficiency) by treating each transit as independent.

  • Another relatively easy way would be to fit two t0s for each transit and offset the second dataset such that they overlap (i.e. the model sees two sets of transits at t0...), while leaving the period continuous. But this seems annoying to model with a GP and I'm not sure how exoplanet/PyMC works with a parameter directly modifying the input data..?

  • A final way might be to, instead of modelling a period, model some time difference plus a scalar integer period ratio (eg (t_1-t_0)/n where n={1,2,3,4,...,max}). This would seem to be the best solution. But I don't know if using non-continuous variables is implementable and/or differentiable as PyMC requires?

Any ideas which of these options seems most reasonable, and also how simple it might be to hack exoplanet to fit them?

Where should +Z point in the Cartesian XYZ frame convention?

Following #15, we should have a discussion about the pros and cons of various orientations of the XYZ Cartesian frame in observer coordinates. If we are sticking to right-handed systems, there are basically two options

  1. +Z points away from the observer. This is currently the way exoplanet is structured.
  2. +Z points toward the observer. This is the way Murray and Corriea seem to have things.

There are some pros and cons to each system.

System 1 yields v_z = v_radial, in that positive radial velocities are receeding from the observer. However, matching up astrometric orbits is a little bit tricky because after we go from the x,y,z in the perifocal frame to the XYZ in the sky frame, we need to associate X and Y with RA and DEC appropriately. The observational definition is that if the position angle of the ascending node (Omega) is 0, then the node is located 0 degrees east of north. This makes the X axis correspond to the north axis, and so to conform to the right hand rule, since we usually draw north up and east to the left, the Y axis must point to the right, or Y = - RA axis. This means that mentally mapping a simple orbit in the perifocal plane to the observer frame requires a reflection about the X axis. I always found this a little bit confusing because when plotting orbits we typically invert the RA axis to have it increasing to the left. This convention would also have us negating it as well.

System 2 does better with associating X and Y with cardinal directions. The X axis still corresponds to north, but now Y points to the left, so Y = + RA axis. So mapping from the perifocal frame to the observer frame only requires a 90 degree rotation. This is easy enough to plot an orbit as plt.plot(Y, X). This system can't have it all, though, since now Z points towards the observer and so v_radial = - v_z.

Thoughts on which would make our lives easiest across all of the package?

Second, following the choice of where +Z points, is a discussion of which node is the ascending node, i.e., the one where the secondary is coming towards the observer, or the one where the secondary is receeding from the observer.

A third observational convention is that inclination is [0, 90) degrees when the orbit is counter-clockwise on the sky and (90, 180] degrees when the orbit is clockwise on the sky. For this reason I tend to prefer System 2 because it means the inclination is simply the angle between the angular momentum vector of the orbit and the +Z axis (pointed towards the observer).

Once we make decisions on the first two conventions then I think the rotation matrices from the perifocal frame to the observer frame are uniquely specified.

TTVOrbit for StarryLightCurve

Thanks for developing the exoplanet package!

I wonder is the TTVOrbit option on? If so, does the StarryLightCurve work properly on it?

  • TTVOrbit with the argument transit_times gives me an error message, e.g.,
orbit = xo.orbits.TTVOrbit(transit_times=[1,2,3])
  • The ttvs argument works, but when I input the orbit to Starry, it gives an error message.

Here is what I tried.

orbit = xo.orbits.TTVOrbit(t0=[0],period=[10],ttvs=[1,2,3])
t = np.linspace(0, 20, 100)
u = [0.3, 0.2]
light_curve = xo.StarryLightCurve(u).get_light_curve(orbit=orbit, r=0.1, t=t, texp=0.02).eval()

Here is what I got.

---------------------------------------------------------------------------
IndexError                                Traceback (most recent call last)
<ipython-input-12-e5d5e061674b> in <module>
      2 t = np.linspace(0, 20, 100)
      3 u = [0.3, 0.2]
----> 4 light_curve = xo.StarryLightCurve(u).get_light_curve(orbit=orbit, r=0.1, t=t, texp=0.02).eval()

~/.conda/envs/env_WaLEs/lib/python3.7/site-packages/exoplanet-0.1.7.dev0-py3.7.egg/exoplanet/light_curves.py in get_light_curve(self, orbit, r, t, texp, oversample, order, use_in_transit)
     89             transit_model = tt.shape_padleft(tt.zeros_like(r), t.ndim) \
     90                 + tt.shape_padright(tt.zeros_like(t), r.ndim)
---> 91             inds = orbit.in_transit(t, r=r, texp=texp)
     92             t = t[inds]
     93 

~/.conda/envs/env_WaLEs/lib/python3.7/site-packages/exoplanet-0.1.7.dev0-py3.7.egg/exoplanet/orbits/keplerian.py in in_transit(self, t, r, texp)
    614         # Wrap the times into time since transit
    615         hp = 0.5 * self.period
--> 616         dt = tt.mod(self._warp_times(t) - self.t0 + hp, self.period) - hp
    617 
    618         if self.ecc is None:

~/.conda/envs/env_WaLEs/lib/python3.7/site-packages/exoplanet-0.1.7.dev0-py3.7.egg/exoplanet/orbits/ttv.py in _warp_times(self, t)
     78         ind = tt.cast(tt.floor(delta), "int64")
     79         dt = tt.stack([ttv[tt.clip(ind[i], 0, ttv.shape[0]-1)]
---> 80                        for i, ttv in enumerate(self.ttvs)], -1)
     81         return tt.shape_padright(t) + dt

~/.conda/envs/env_WaLEs/lib/python3.7/site-packages/exoplanet-0.1.7.dev0-py3.7.egg/exoplanet/orbits/ttv.py in <listcomp>(.0)
     78         ind = tt.cast(tt.floor(delta), "int64")
     79         dt = tt.stack([ttv[tt.clip(ind[i], 0, ttv.shape[0]-1)]
---> 80                        for i, ttv in enumerate(self.ttvs)], -1)
     81         return tt.shape_padright(t) + dt

~/.conda/envs/env_WaLEs/lib/python3.7/site-packages/theano/tensor/var.py in __getitem__(self, args)
    508         # Check if the number of dimensions isn't too large.
    509         if self.ndim < index_dim_count:
--> 510             raise IndexError('too many indices for array')
    511 
    512         # Convert an Ellipsis if provided into an appropriate number of

IndexError: too many indices for array

I may just misunderstand the functions. Is there any demo code I could check using TTVOrbit on StarryLightCurve?

Thanks a lot!

error when using a small number of samples in tune

When the number of samples in tune is small (lower than start/finish perhaps), I got an error

sampler = xo.PyMC3Sampler(window=100, start=200, finish=200)
with model:
    sampler.tune(tune=200, start=map_soln, step_kwargs=dict(target_accept=0.9))

/home/tom/anaconda3/lib/python3.6/site-packages/theano/tensor/subtensor.py:2197: FutureWarning: Using a non-tuple sequence for multidimensional indexing is deprecated; use `arr[tuple(seq)]` instead of `arr[seq]`. In the future this will be interpreted as an array index, `arr[np.array(seq)]`, which will result either in an error or a different result.
  rval = inputs[0].__getitem__(inputs[1:])
Sampling 4 chains: 100%|██████████| 808/808 [21:22<00:00,  4.35s/draws]
Sampling 4 chains: 100%|██████████| 808/808 [19:31<00:00,  4.24s/draws]

AttributeError                            Traceback (most recent call last)
<ipython-input-17-77b8d9feadf0> in <module>()
      1 sampler = xo.PyMC3Sampler(window=100, start=200, finish=200)
      2 with model:
----> 3     sampler.tune(tune=200, start=map_soln, step_kwargs=dict(target_accept=0.9))

~/gitcode/exoplanet/exoplanet/sampling.py in tune(self, tune, start, step_kwargs, **kwargs)
    180         self.extend_tune(start=start, step_kwargs=step_kwargs,
    181                          steps=self.finish, trace=trace, **kwargs)
--> 182         self._current_step.stop_tuning()
    183 
    184         return self._current_trace

AttributeError: 'NoneType' object has no attribute 'stop_tuning'

GP implementation of a precomputed covariance matrix

PyMC3 supports multiplying a GP covariance function with a precomputed covariance matrix, as seen with the example below copied from the PyMC3 docs (https://docs.pymc.io/notebooks/GP-MeansAndCovs.html). Is it possible to do something similar directly with exoplanet?

Thanks!

import matplotlib.pyplot as plt
import pymc3 as pm
import theano
import numpy as np

# data vector
X = np.linspace(0, 2, 200)[:,None]

# first evaluate a covariance function into a matrix
period = 0.2
cov_cos = pm.gp.cov.Cosine(1, period)
K_cos = theano.function([], cov_cos(X))()

# now multiply it with a covariance *function*
cov = pm.gp.cov.Matern32(1, 0.5) * K_cos

K = cov(X).eval()

plt.figure(figsize=(14,4))
plt.plot(X, pm.MvNormal.dist(mu=np.zeros(K.shape[0]), cov=K).random(size=3).T);
plt.title("Samples from the GP prior");
plt.ylabel("y");
plt.xlabel("X");

image

Can KeplerianOrbit functions accept theano variables?

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", ((2np.pismac.au.value/P)(1 + eccorbit.cosf)/np.sqrt(1-ecc**2))/(Rsc.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:

  • exoplanet 0.2.1.dev0
  • Operating system: MacOSX
  • Python version & installation method (pip, conda, etc.): pip+git

getting OSError: [Errno 12] Cannot allocate memory

This is really a question on sampling. I am following your examples but using a three planet system.
I ran the tuning step with

np.random.seed(42)
sampler = xo.PyMC3Sampler(window=100, start=500, finish=200)
with model:
    burnin = sampler.tune(tune=5000, start=map_soln, step_kwargs=dict(target_accept=0.9), cores=6)

/home/tom/anaconda3/lib/python3.6/site-packages/theano/tensor/subtensor.py:2197: FutureWarning: Using a non-tuple sequence for multidimensional indexing is deprecated; use `arr[tuple(seq)]` instead of `arr[seq]`. In the future this will be interpreted as an array index, `arr[np.array(seq)]`, which will result either in an error or a different result.
  rval = inputs[0].__getitem__(inputs[1:])
Sampling 6 chains: 100%|██████████| 3012/3012 [1:50:47<00:00,  6.27s/draws]  
Sampling 6 chains: 100%|██████████| 612/612 [07:34<00:00,  2.57s/draws]
Sampling 6 chains: 100%|██████████| 1212/1212 [15:29<00:00,  2.65s/draws]
Sampling 6 chains: 100%|██████████| 2412/2412 [1:18:08<00:00,  8.68s/draws]
Sampling 6 chains: 100%|██████████| 4812/4812 [2:54:50<00:00, 10.13s/draws]  
Sampling 6 chains: 100%|██████████| 18012/18012 [12:41:31<00:00, 10.79s/draws]  

When I try to do the sampling step I get

with model:
    trace = sampler.sample(draws=2000, cores=6)

Multiprocess sampling (6 chains in 6 jobs)
NUTS: [logw0, logS0, logs2, omega, ecc, rb, t0, logP, rho_star, r_star, u_star, mean]
/home/tom/anaconda3/lib/python3.6/site-packages/theano/tensor/subtensor.py:2197: FutureWarning: Using a non-tuple sequence for multidimensional indexing is deprecated; use `arr[tuple(seq)]` instead of `arr[seq]`. In the future this will be interpreted as an array index, `arr[np.array(seq)]`, which will result either in an error or a different result.
  rval = inputs[0].__getitem__(inputs[1:])

---------------------------------------------------------------------------
OSError                                   Traceback (most recent call last)
<ipython-input-29-82e710a8cf2e> in <module>()
      1 with model:
----> 2     trace = sampler.sample(draws=1000, cores=6)

~/anaconda3/lib/python3.6/site-packages/exoplanet/sampling.py in sample(self, trace, step, start, step_kwargs, **kwargs)
    191             start=start, step_kwargs=step_kwargs, trace=trace, step=step)
    192         kwargs["tune"] = self.finish
--> 193         return pm.sample(start=start, step=step, **kwargs)

~/anaconda3/lib/python3.6/site-packages/pymc3/sampling.py in sample(draws, step, init, n_init, start, trace, chain_idx, chains, cores, tune, nuts_kwargs, step_kwargs, progressbar, model, random_seed, live_plot, discard_tuned_samples, live_plot_kwargs, compute_convergence_checks, use_mmap, **kwargs)
    447             _print_step_hierarchy(step)
    448             try:
--> 449                 trace = _mp_sample(**sample_args)
    450             except pickle.PickleError:
    451                 _log.warning("Could not pickle model, sampling singlethreaded.")

~/anaconda3/lib/python3.6/site-packages/pymc3/sampling.py in _mp_sample(draws, tune, step, chains, cores, chain, random_seed, start, progressbar, trace, model, use_mmap, **kwargs)
    994         sampler = ps.ParallelSampler(
    995             draws, tune, chains, cores, random_seed, start, step,
--> 996             chain, progressbar)
    997         try:
    998             with sampler:

~/anaconda3/lib/python3.6/site-packages/pymc3/parallel_sampling.py in __init__(self, draws, tune, chains, cores, seeds, start_points, step_method, start_chain_num, progressbar)
    273             ProcessAdapter(draws, tune, step_method,
    274                            chain + start_chain_num, seed, start)
--> 275             for chain, seed, start in zip(range(chains), seeds, start_points)
    276         ]
    277 

~/anaconda3/lib/python3.6/site-packages/pymc3/parallel_sampling.py in <listcomp>(.0)
    273             ProcessAdapter(draws, tune, step_method,
    274                            chain + start_chain_num, seed, start)
--> 275             for chain, seed, start in zip(range(chains), seeds, start_points)
    276         ]
    277 

~/anaconda3/lib/python3.6/site-packages/pymc3/parallel_sampling.py in __init__(self, draws, tune, step_method, chain, seed, start)
    180             draws, tune, seed)
    181         # We fork right away, so that the main process can start tqdm threads
--> 182         self._process.start()
    183 
    184     @property

~/anaconda3/lib/python3.6/multiprocessing/process.py in start(self)
    103                'daemonic processes are not allowed to have children'
    104         _cleanup()
--> 105         self._popen = self._Popen(self)
    106         self._sentinel = self._popen.sentinel
    107         # Avoid a refcycle if the target function holds an indirect

~/anaconda3/lib/python3.6/multiprocessing/context.py in _Popen(process_obj)
    221     @staticmethod
    222     def _Popen(process_obj):
--> 223         return _default_context.get_context().Process._Popen(process_obj)
    224 
    225 class DefaultContext(BaseContext):

~/anaconda3/lib/python3.6/multiprocessing/context.py in _Popen(process_obj)
    275         def _Popen(process_obj):
    276             from .popen_fork import Popen
--> 277             return Popen(process_obj)
    278 
    279     class SpawnProcess(process.BaseProcess):

~/anaconda3/lib/python3.6/multiprocessing/popen_fork.py in __init__(self, process_obj)
     17         util._flush_std_streams()
     18         self.returncode = None
---> 19         self._launch(process_obj)
     20 
     21     def duplicate_for_child(self, fd):

~/anaconda3/lib/python3.6/multiprocessing/popen_fork.py in _launch(self, process_obj)
     64         code = 1
     65         parent_r, child_w = os.pipe()
---> 66         self.pid = os.fork()
     67         if self.pid == 0:
     68             try:

OSError: [Errno 12] Cannot allocate memory

expanding `_get_consistent_inputs` for astrometric + rv parameter sets

For an orbit fit to both double-lined RV and astrometric data of a binary star, you basically fully orient everything in 3D space and determine the individual stellar masses. (There is no information about the radii of the stars). So, we'd like to choose a basic set of (independent) orbital parameters from which every other property can be derived.

Here are a few choices that might work

a (physical)
period
kappa = a1/(a1 + a2)
-> Mtot is calculated from (a, period)
-> M2 is calculated from Mtot * kappa
-> M2 is calculated from Mtot - M2

Alternatively,

a (physical)
period
M2
-> Mtot is calculated from (a, period)
-> M1 is calculated from Mtot - M2

or

a (physical)
period
q = M2/M1
-> Mtot is calculated from (a, period)
-> M1 is calculated from Mtot and q
-> M2 is calculated from Mtot and q

As currently written, if a, period, and m_star are given, _get_consistent_inputs blocks us.

I think option #2 parameterization is going to be the easiest to merge with the current code. So, if we give a, period, and m_planet = M2, then exoplanet calculates rho_star assuming r_star = 1.0.

Then, this flows down the control and m_star is calculated from rho_star and r_star, but, because r_star = 1.0, this isn't the right stellar mass implied by a, period, and m_planet.

# conversion constant
au_to_R_sun = (constants.au / constants.R_sun).value

a_ang = 0.324 # arcsec
parallax = 24.05 # milliarcsec
dpc = 1e3 / parallax
a = a_ang * dpc # au

e = 0.798
i = 96.0 * deg # [rad]
omega = 251.6 * deg - np.pi # Pourbaix reports omega_2, but we want omega_1
Omega = 159.6 * deg
P = 28.8 * 365.25 # days

T0 = Time(1989.92, format="decimalyear")
T0.format = "jd"
T0 = T0.value # [Julian Date]

gamma = 47.97 # km/s; systemic velocity

# kappa = a1 / (a1 + a2)
kappa = 0.45

# calculate Mtot from a, P
Mtot = calc_Mtot(a, P)

M2 = kappa * Mtot
M1 = Mtot - M2

# calculate M1, M2 from kappa and Mtot
print("Masses from Kepler's third law M1:{:.2f}, M2:{:.2f}, Mtot:{:.2f}".format(M1, M2, Mtot))

# n.b. that we include an extra conversion for a, because exoplanet expects a in R_sun
orbit = xo.orbits.KeplerianOrbit(a=a * au_to_R_sun, t_periastron=T0, period=P,
                               incl=i, ecc=e, omega=omega, Omega=Omega, m_planet=M2)


print("m_planet: {:.2f}".format(orbit.m_planet.eval()))
print("m_star: {:.2f}".format(orbit.m_star.eval()))
print("m_tot: {:.2f}".format(orbit.m_total.eval()))

# Output
Correct masses from Kepler's third law M1:1.62, M2:1.33, Mtot:2.95
Calculated masses from _get_consistent_inputs
m_planet: 1.33
m_star: 0.27
m_tot: 1.60

I'm not super familiar with the minimum parameter sets required for transit fitting, or transit + RV fitting, so I'm a little hesitant to dive in and start disrupting control flow through this routine. I would be inclined to change things such that the user specifies a, P, and m_planet and instead of calculating rho_star directly, _get_consistent_inputs calculates m_total and m_star = m_total - m_planet. I could use a little guidance on whether this is advisable package-wide, however.

better eccentricity distribution?

Currently the orbit-fitting tutorials use xo.distributions.UnitUniform or Uniform as a prior for orbital eccentricity, but this can bias the resulting maximum-likelihood estimate to high e. Are there plans to add any built-in distributions that would make better priors on eccentricity? One option might be a beta-distribution (Kipping 2014).

Planetary coordinate system in multi-planet case

@dfm In the multi-planet case, in which coordinate system do you solve Kepler's equations? The ideal is to use Jacobi coordinates, although heliocentric and barycentric are easier to implement. However, both heliocentric and barycentric have drawbacks.

Luminosity of secondary

Using exoplanet is great! I'm now wondering, starry has the capability to make the secondary luminous, meaning that it will automatically add secondary eclipses, and phase curves. Is this functionality preserved in exoplanet? How can I add secondary eclipses to my fit exoplanet model?

Many thanks!
C

K2-24 tutorial does not converge

I did a fresh conda environment and installed exoplanet and dependencies. I ran the K2-24 tutorial in a downloaded notebook and I get the following from the sample

with model:

    trace = sampler.sample(draws=2000)

Sampling 4 chains: 100%|██████████| 8000/8000 [00:37<00:00, 215.91draws/s]
ERROR:pymc3:The chain contains only diverging samples. The model is probably misspecified.
WARNING:pymc3:The acceptance probability does not match the target. It is 0.0, but should be close to 0.9. Try to increase the number of tuning steps.
ERROR:pymc3:There were 1999 divergences after tuning. Increase `target_accept` or reparameterize.
WARNING:pymc3:The acceptance probability does not match the target. It is 9.77952587176284e-17, but should be close to 0.9. Try to increase the number of tuning steps.
ERROR:pymc3:There were 1907 divergences after tuning. Increase `target_accept` or reparameterize.
WARNING:pymc3:The acceptance probability does not match the target. It is 0.006233605808817455, but should be close to 0.9. Try to increase the number of tuning steps.
ERROR:pymc3:There were 1918 divergences after tuning. Increase `target_accept` or reparameterize.
WARNING:pymc3:The acceptance probability does not match the target. It is 6.736254608770937e-15, but should be close to 0.9. Try to increase the number of tuning steps.
ERROR:pymc3:The gelman-rubin statistic is larger than 1.4 for some parameters. The sampler did not converge.
ERROR:pymc3:The estimated number of effective samples is smaller than 200 for some parameters.

The summary output is similarly bad

  | mean | sd | mc_error | hpd_2.5 | hpd_97.5 | n_eff | Rhat
-- | -- | -- | -- | -- | -- | -- | --
42.363591 | 0.000301 | 0.000030 | 42.363218 | 42.364071 | NaN | 7.343693
20.885450 | 0.000205 | 0.000021 | 20.885138 | 20.885673 | 2.018893 | 12.642319
0.063015 | 0.000591 | 0.000059 | 0.062178 | 0.063959 | NaN | 6.166608
0.046068 | 0.000590 | 0.000059 | 0.045105 | 0.046691 | 2.028321 | 10.391917
29.863008 | 2.787906 | 0.278459 | 25.333686 | 33.197244 | 2.035109 | 8.949308
24.414847 | 2.274746 | 0.227115 | 22.034980 | 28.812748 | NaN | 3.937818
0.119697 | 0.069082 | 0.006903 | 0.004913 | 0.187045 | 2.018684 | 12.257648
0.228654 | 0.043500 | 0.004344 | 0.166656 | 0.279087 | NaN | 9.578997
0.575165 | 0.042205 | 0.004219 | 0.522428 | 0.645138 | NaN | 7.824198
0.607800 | 0.057467 | 0.005744 | 0.519003 | 0.679079 | 2.014646 | 13.941053

Aperiodic error when computing xo.orbits.KerplerianOrbit

Describe the bug
When creating an instance of the xo.orbits.KeplerianOrbit object, I sometimes incur an error that reports ValueError: Cannot compute test value ... missing default value. [see full traceback below].

I always provide the exact same values every time; but every ~5th try, Theano triggers this error.

The behaviour of xo.orbits.KeplerianOrbit is such that we can either enter the inclination, impact parameter, or transit duration; but only one a time. [this is built in; no problem]

The problem arises that I've tried this exact same code on 3 machines, and ~10 times per machine. Depending on the machine and time of day (i.e. randomness), the KeplerianOrbit object does not accept the form of one or more values in the set [inclination, impact parameter, duration] that I provide. Note that I always give it the exact same values every time [see example code below].

If I comment out the inclination and uncomment either the duration or impact parameter line [or vice versa], then the code will work for ~5 more attempts; at some point, it will fail again with the related error. At which point, I would permute through the list and KeplerianOrbit nominally works again.

Note that the exact form of the error changes related to the kwarg used [inc, b, dur]; but 99% of the error is identical to the traceback below.

To Reproduce

import exoplanet as xo
import numpy as np

def build_synthetic_model(min_phase=-0.1, max_phase=0.1, size=1000, u=[0.0]):
    ''' Create synthetic transit light curve by grabbing data from exo.mast
         and populating xo.orbits.KeplerianOrbit with that

        min_phase (float): minimum value for the transit light curve phase
        max_phase (float): maximum value for the transit light curve phase
        size (int): number of data points
        u (list): limb darkening coefficients; default: no limb darkening
    '''
    # Planetary orbital parameters
    t0 = 54278.93671400007
    period = 2.21857567
    Rp_Rs = 0.1546877455085016
    deg2rad = np.pi / 180

    orbit = xo.orbits.KeplerianOrbit(
        t0=t0,
        period=period,
        a=8.83602,
        # b=0.6631,
        # incl=1.4959217018843398,
        duration=0.0759722,
        ecc=0.0,
        omega=1.5707963267948966  # == pi/2
    )
    star = xo.LimbDarkLightCurve(u)
    phase = np.linspace(min_phase, max_phase, size)
    times = phase * period + t0
    lightcurve = star.get_light_curve(
        r=Rp_Rs,
        orbit=orbit,
        t=times)

    lightcurve = lightcurve.eval().flatten()

    return times, lightcurve, orbit

if __name__ == '__main__':
    times, lightcurve, orbit = build_synthetic_model()

Expected behavior
I expect to generate a transit model, via KeplerianOrbit, every time I run the code above.

Your setup (please complete the following information):

  • Version of exoplanet: 0.2.6.dev1+gcc717ef
  • Operating system: Ubuntu 18.04 (Linux 5.3.0-40-generic #32~18.04.1-Ubuntu x86_64)
  • Python version & installation method (pip, conda, etc.): Python 3.7.6
    Conda version: 4.8.2
    Method: exoplanet install through github repo; all dependencies through conda

Additional context
In writing this issue, I realized that I am using the word model inside the function; but that I call this function inside a PyMC3 wrapper which includes with pm.Model() as model. This could easily be the issue. I will report back after changing the variable name that I store the xo.LimbDarkLightCurve(u).get_light_curve instance.

Note that the tracback for this error points to the KeplerianOrbit; but I will change the instance name and respond if I do not receive the error again.

Full Traceback

---------------------------------------------------------------------------
AttributeError                            Traceback (most recent call last)
/usr/local/lib/python3.6/dist-packages/theano/gof/op.py in __call__(self, *inputs, **kwargs)
    624                 try:
--> 625                     storage_map[ins] = [self._get_test_value(ins)]
    626                     compute_map[ins] = [True]

7 frames
/usr/local/lib/python3.6/dist-packages/theano/gof/op.py in _get_test_value(cls, v)
    580         detailed_err_msg = utils.get_variable_trace_string(v)
--> 581         raise AttributeError('%s has no test value %s' % (v, detailed_err_msg))
    582 

AttributeError: Elemwise{mul,no_inplace}.0 has no test value  
Backtrace when that variable is created:

  File "/usr/local/lib/python3.6/dist-packages/ipykernel/ipkernel.py", line 196, in do_execute
    res = shell.run_cell(code, store_history=store_history, silent=silent)
  File "/usr/local/lib/python3.6/dist-packages/ipykernel/zmqshell.py", line 533, in run_cell
    return super(ZMQInteractiveShell, self).run_cell(*args, **kwargs)
  File "/usr/local/lib/python3.6/dist-packages/IPython/core/interactiveshell.py", line 2718, in run_cell
    interactivity=interactivity, compiler=compiler, result=result)
  File "/usr/local/lib/python3.6/dist-packages/IPython/core/interactiveshell.py", line 2822, in run_ast_nodes
    if self.run_code(code, result):
  File "/usr/local/lib/python3.6/dist-packages/IPython/core/interactiveshell.py", line 2882, in run_code
    exec(code_obj, self.user_global_ns, self.user_ns)
  File "<ipython-input-9-866db5884fd4>", line 2, in <module>
    times, synthetic_eclipse, orbit = build_synthetic_model(size=n_pts)
  File "<ipython-input-5-cb8f0767a3aa>", line 25, in build_synthetic_model
    omega=1.5707963267948966  # == pi/2
  File "/usr/local/lib/python3.6/dist-packages/exoplanet/orbits/keplerian.py", line 207, in __init__
    self.cos_incl = self.dcosidb * self.b


During handling of the above exception, another exception occurred:

ValueError                                Traceback (most recent call last)
<ipython-input-16-93c8d291cba6> in <module>()
      2                                                 log_Q=np.log(1 / np.sqrt(2)),
      3                                                 tune=tune, draws=draws,
----> 4                                                 target_accept=0.9)

<ipython-input-4-92c977bcc44c> in run_pymc3_with_gp(times, data, dataerr, orbit, log_Q, tune, draws, target_accept, u)
     30             return pm.math.sum(light_curves, axis=-1) + mean
     31 
---> 32         light_curves = lc_model(times)
     33 
     34         # The likelihood function assuming known Gaussian uncertainty

<ipython-input-4-92c977bcc44c> in lc_model(times)
     24             # Compute the model light curve using starry
     25             light_curves = xo.LimbDarkLightCurve(u).get_light_curve(
---> 26                 orbit=orbit, r=r, t=times
     27             )
     28 

/usr/local/lib/python3.6/dist-packages/exoplanet/light_curves.py in get_light_curve(self, orbit, r, t, texp, oversample, order, use_in_transit)
    103                 tt.zeros_like(r), t.ndim
    104             ) + tt.shape_padright(tt.zeros_like(t), r.ndim)
--> 105             inds = orbit.in_transit(t, r=r, texp=texp)
    106             t = t[inds]
    107 

/usr/local/lib/python3.6/dist-packages/exoplanet/orbits/keplerian.py in in_transit(self, t, r, texp)
    605                 self.cos_omega,
    606                 self.sin_omega,
--> 607                 self.cos_incl + z,
    608                 self.sin_incl + z,
    609                 R + r,

/usr/local/lib/python3.6/dist-packages/theano/tensor/var.py in __add__(self, other)
    126     def __add__(self, other):
    127         try:
--> 128             return theano.tensor.basic.add(self, other)
    129         # We should catch the minimum number of exception here.
    130         # Otherwise this will convert error when Theano flags

/usr/local/lib/python3.6/dist-packages/theano/gof/op.py in __call__(self, *inputs, **kwargs)
    637                         raise ValueError(
    638                             'Cannot compute test value: input %i (%s) of Op %s missing default value. %s' %
--> 639                             (i, ins, node, detailed_err_msg))
    640                     elif config.compute_test_value == 'ignore':
    641                         # silently skip test

ValueError: Cannot compute test value: input 0 (Elemwise{mul,no_inplace}.0) of Op Elemwise{add,no_inplace}(Elemwise{mul,no_inplace}.0, Elemwise{second,no_inplace}.0) missing default value.  
Backtrace when that variable is created:

  File "/usr/local/lib/python3.6/dist-packages/ipykernel/ipkernel.py", line 196, in do_execute
    res = shell.run_cell(code, store_history=store_history, silent=silent)
  File "/usr/local/lib/python3.6/dist-packages/ipykernel/zmqshell.py", line 533, in run_cell
    return super(ZMQInteractiveShell, self).run_cell(*args, **kwargs)
  File "/usr/local/lib/python3.6/dist-packages/IPython/core/interactiveshell.py", line 2718, in run_cell
    interactivity=interactivity, compiler=compiler, result=result)
  File "/usr/local/lib/python3.6/dist-packages/IPython/core/interactiveshell.py", line 2822, in run_ast_nodes
    if self.run_code(code, result):
  File "/usr/local/lib/python3.6/dist-packages/IPython/core/interactiveshell.py", line 2882, in run_code
    exec(code_obj, self.user_global_ns, self.user_ns)
  File "<ipython-input-9-866db5884fd4>", line 2, in <module>
    times, synthetic_eclipse, orbit = build_synthetic_model(size=n_pts)
  File "<ipython-input-5-cb8f0767a3aa>", line 25, in build_synthetic_model
    omega=1.5707963267948966  # == pi/2
  File "/usr/local/lib/python3.6/dist-packages/exoplanet/orbits/keplerian.py", line 207, in __init__
    self.cos_incl = self.dcosidb * self.b

Recommend Projects

  • React photo React

    A declarative, efficient, and flexible JavaScript library for building user interfaces.

  • Vue.js photo Vue.js

    🖖 Vue.js is a progressive, incrementally-adoptable JavaScript framework for building UI on the web.

  • Typescript photo Typescript

    TypeScript is a superset of JavaScript that compiles to clean JavaScript output.

  • TensorFlow photo TensorFlow

    An Open Source Machine Learning Framework for Everyone

  • Django photo Django

    The Web framework for perfectionists with deadlines.

  • D3 photo D3

    Bring data to life with SVG, Canvas and HTML. 📊📈🎉

Recommend Topics

  • javascript

    JavaScript (JS) is a lightweight interpreted programming language with first-class functions.

  • web

    Some thing interesting about web. New door for the world.

  • server

    A server is a program made to process requests and deliver data to clients.

  • Machine learning

    Machine learning is a way of modeling and interpreting data that allows a piece of software to respond intelligently.

  • Game

    Some thing interesting about game, make everyone happy.

Recommend Org

  • Facebook photo Facebook

    We are working to build community through open source technology. NB: members must have two-factor auth.

  • Microsoft photo Microsoft

    Open source projects and samples from Microsoft.

  • Google photo Google

    Google ❤️ Open Source for everyone.

  • D3 photo D3

    Data-Driven Documents codes.