sblunt / orbitize Goto Github PK

This should be a field in the System class that lists the labels (e.g., 'a', 'e') for each parameter so we (humans) know which parameter each is. We could consider using 'a1', 'e2' to help separate multiple planets.

Make sure we can install and compile it.

Allow multiple planets to be fit at once (useful for stellar RVs & astrometry)

From Rob: A nice functionality might be to specify what parameter you want on the color axis. while date can look cool, so can things like eccentricity (as jason’s nice plot of 51 eri shows)

The API (written by @sarah Blunt) specifies that the posterior sample array (or orbits) is a 2D array with shape (num_parameters, num_samples)

However, for MCMC output and input to corner etc. the 2D posterior sample array usually has shape (num_samples,num_parameters).

For now, I will proceed with the MCMC and corner.py standard but if OFTI is already implemented the opposite way, we can switch or discuss which way we should go.

One other reason to choose the MCMC/corner.py way (num_samples,num_parameters) is that in numpy, a 2-D array with shape (M,N) has M rows and N columns. So basically, it is a set of M arrays that each of length N. If M = num_samples and N=num_parameters, then each of the M arrays is a complete orbit on its own, whether computed by OFTI or MCMC. So MxN means M orbits each with N parameters, which makes more intuitive sense to me!

I have the code. I need to put it in and make it usable. A stretch goal is to add Rob's eccentricity solver

• turn for loops into vector array operations

This bug was causing a bunch of my unit tests on Driver to fail when I was testing for the ability to accept either a filename or an astropy.table.Table object.

When a Driver object is created with the OFTI algorithm, initializing the Sampler.OFTI object will convert all of the ra/dec input to sep/pa. This changes the quant1, quant1_err, quant2, and quant2_err values. However, it does not change the last column, which is quant_type. That entry remains radec.

This causes problems when I (or another user) tries to access the data table through the sampler.data_table attribute because the wrong information was shown.

There three potential fixes:

1. Add one line to the for-loop where the ra/dec values are being converted to also change this last column (from 'radec' to 'seppa'). So, accessing the data_table at any time will provide consistent values, however, it will not be the same as the input (i.e. you get the table used to by OFTI but not the input table for troubleshooting).

2. Change lines 97-100 of sampler.py (where astropy.table columns are being assigned to local attributes) to be copying the values instead of just assignment. Since Table objects are passed by reference, not by value, when the variables to which these columns were assigned are changed later, it changes the original too. In addition, with this solution, we would have two separate tables stored for the sampler object: input_table showing the original input and data_table for showing the values used by OFTI to calculate.

3. Do the copying thing as in option 2, but do not expose the data_table at all to users to keep it behind the scenes.

I personally like option 2 best but option 3 also makes sense if you want to keep it under the hood. I can implement any of these but it would be good to discuss ideas.

Referenced in PR #51.

Sometimes we have measurements that don't fit the standard orbitize format yet (e.g., Gaia/Hipparcos astrometry), and it would be nice to fit it as a one-off before making it a full feature (in case the format is supposed to change). It would be nice to be able to pass into the system(?) object, an extra likelihood function that is custom by the user. It should take a 2-D array of fit parameters, and return an array of likelihood values.

Especially for OFTI, it'd be great to pass in arrays of orbital parameters into kepler.calc_orbit.

Basically the epochs argument can be of length M and the sma, ecc, etc arguments can be of length N, and the return values (e.g., raoff) should be of dimension NxM.

This will requiring formatting the code to angle scalars/1-D vectors at the same time, as well as 1-D/2-D vectors at the same time. It would be great to avoid code duplication. I think this is possible.

We should make system.py and sampler.py use the same ordering of orbital elements, and make it less arbitrary. Henry suggested the following order that makes more sense:

a (semi-major axis)
e (eccentricity)
i (inclination)
omega or w (arg of pericentre)
Omega or O (long. of ascending node)
tau (or Mean anomaly or whatever equivalent coordinate one uses to denote a specific time from a specific point in orbit)

Code for both samplers as well as system.py need to be modified to support this.

Add policy from slack to contributor guidelines on how to collaborate on a single branch.

From #35

Add a function in the MCMC sampler to inspect the MCMC chains (either visually or by some metrics) and allow the user to chop off additional steps as "burn-in" if part of the MCMC chain hasn't converged yet. Update the results object with the updated chain.

To be completed after multi-threaded MCMC is complete (#30)

Without RVs, imaging astrometry has argument of periastron and longitude of ascending node doubly peaked, with the peaks 180 degrees apart. The distribution is identical between [0,180) and [180,360). Sampling multimodal distributions with MCMC is annoying, so it'd be nice if we had the option to only sample in the [0,180) degree range when we want to be lazy.

Technically, I think this just requires an optional argument and changing the uniform bounds on the angles, so it should not be too hard.

A stretch goal can be to estimate quickly (via maybe LSQR, or the burn in from MCMC) where the peaks are, and center the sampling around the peak (e.g., [peak-90, peak+90) degrees). This would prevent the edge case that if the peak is centered about zero, you'd still have a multimodal distribution by construction of the prior bounds.

The issue is that MCMC doesn't need to run a bunch of orbits at once normally, whereas OFTI does. However, if we put "number of models" as the leading dimension, then everything needs to be index it explicitly (e.g., raoff[model_number, epoch_number]). It would make things easier if we put "number of models" as the last dimension so that we don't have to think about it when not programming explicitly OFTI-related things (e.g., raoff[epoch_number] could be a number or an array).

lnlike, as per the docstring, should run with an extra dimension being different orbital parameter sets. @isabelangelo says this is not working right now so we should fix that. She will update her branch, and we can reproduce the bug running OFTI from her branch.

Somewhat related to #61. We should have an option for a user defined start date for tau (currently it's starting from MJD=0). There is possible evidence that having tau near the date of current observations helps with MCMC convergence.

chi2_lnlike function should compute the difference between an model (ang1) and observed (ang2) position angle as

ang1-ang2 = np.arctan2(np.sin(ang1 - ang2), np.cos(ang1 - ang2))

so the function will have to be know what the unit of the data is. If it's of type u.rad/u.deg, then it has to use this rather than the simple data-model.

Would be great to have more than one prior! Specifically, adding a uniform bounded prior and a linear prior (with parameters m, b of y=mx+b) (and any others you think might be useful) would be very helpful. For drawing samples from an arbitrary distribution, see: en.wikipedia.org/wiki/Inverse_transform_sampling.

Because dynamical masses are nice. Should be straightforward, but making an issue to keep track of it.

I am on an CentOs 7 machine with the most recent anaconda. numpy and cathon are upgraded.
Installing orbitize (no matter whther fron pip or source) gives the following error:

pip install orbitize
Collecting orbitize
Downloading https://files.pythonhosted.org/packages/32/14/e0da52c27d509cea675794d41cbc6b252c5ea3953d8a9427b386baac04cb/orbitize-1.0.0.tar.gz (76kB)
100% |████████████████████████████████| 81kB 10.3MB/s
Complete output from command python setup.py egg_info:
Traceback (most recent call last):
File "", line 1, in
File "/tmp/pip-install-cwMWKA/orbitize/setup.py", line 41, in
ext_modules=get_extensions(),
File "/tmp/pip-install-cwMWKA/orbitize/setup.py", line 29, in get_extensions
extensions = cythonize([Extension("orbitize._kepler", ["orbitize/_kepler.pyx"])])
File "/usr/lib64/python2.7/site-packages/Cython/Build/Dependencies.py", line 956, in cythonize
aliases=aliases)
File "/usr/lib64/python2.7/site-packages/Cython/Build/Dependencies.py", line 801, in create_extension_list
for file in nonempty(sorted(extended_iglob(filepattern)), "'%s' doesn't match any files" % filepattern):
File "/usr/lib64/python2.7/site-packages/Cython/Build/Dependencies.py", line 111, in nonempty
raise ValueError(error_msg)
ValueError: 'orbitize/_kepler.pyx' doesn't match any files

What could be the reason?
Mny thanks in advance for the help.
-ulrich

It looks like we can initialize system.System with a mass_err and plx_err but they aren't saved as attributes in the class definition so there's no way to access them later (e.g. I'm trying to do so when creating a results.Results object).

On a similar note, it doesn't look like system.System has a mass or plx attribute either, but instead has abs_system_mass and abs_plx? Is this intentional to have different naming? Are these attributes meant to be "public" and accessible outside of the class definition? If so, I think the different names should be mentioned in the docstring to avoid confusion!

See https://arxiv.org/pdf/1809.05490.pdf, awesome work by K. Kosmo O'Neil. Would be great to have these priors as options in orbitize!

After we get the OFTI sampler working (after Issue #4 )

print how many orbits have been accepted out of the desired number of output orbits (standardize this for OFTI and MCMC

At least support CSV for now. Include some documentation on what the format should be.

Get it working for a 1 planet, imaging-only dataset

The RadVel folks had a big discussion about how argp was defined a few months ago, and I wanted to make sure I know how we're defining it. Lauren Weiss suggested I make sure I understand this.

See Figure 1 of BJ's 2017 RadVel paper. Are we defining argp the same way they do, with a 180deg offset, or some other way?

We should be able to marginalize over plate scale/north angle uncertainties.

When running OFTI in Jupyter Lab, the progress bar creates many empty spaces instead of augmenting numerically.

(occurs for the following lines:)

s = myDriver.sampler
orbits = s.run_sampler(1000)

The end result of a bout of orbit-fitting will be an mxn array of orbital parameters, where m is the number of fitted parameters (e.g. semi-major axis), and n is the number of steps (in the case of MCMC) or accepted orbits (in the case of OFTI). Create a new file in orbitize/orbitize that will handle plot creation. Create two classes, one for triangle plots (use the corner Python package for this) and one for orbit plots (see appendix of http://adsabs.harvard.edu/abs/2017AJ....153..229B for examples). Each class should have a main method that creates the appropriate plot given an array of orbital parameters.

This is a relatively big job, and the orbit plots will require using a bunch of orbitize functionality. Maybe a good project for two people to collaborate on?

E.g., travis. Also means we need some tests.

add this after I push MCMC changes to master

We can probably use cython to make the iterative solver for eccentric anomaly faster. It should check though if cython is supported on the user machine, and default back to the python solver if needed.

This might be something @logan-pearce is interested in working on. In particular, I think we should implement the following:

• use a scipy optimization routine to estimate the minimum chi square for a given dataset, then use that to calculate a modified chi square for each tested orbit.
• adaptively modify priors as regions of parameter space are ruled out.

Tau is great for sampling, but somethings annoying to work with. Instead of having each person re-figure out tau every time, it'd be great to have some helper functions in orbitize. I was thinking of the following

• Convert tau to T0 (actual MJD date of periastron passage). Optionally, T0 can be the first periastron passage after a certain date.
• Convert tau to a different tau (we base tau off of MJD=0, but having the flexibility of switching to a different date is also good).

I think the first one is the most important, but the second one could also be useful while we are attacking this problem.

Currently, Driver expects a csv. We should make it flexible enough to accept either a csv or an astropy.table.Table. Thanks @logan-pearce for the suggestion!

create a dictionary in the system class that contains the orbital parameters and their corresponding indices that are output by OFTI/MCMC (useful for plotting)

driver.py takes parallax in arcseconds, while kepler.calc_orbit takes parallax input in mas. Which one should be the standard?

@semaphoreP (me): Make a input data file with the Beta pic b astrometry

Then, @henry-ngo and I can run it, and do some end-to-end validation.

Get it working for 1 planet with only imaging data.

Currently, the tutorial points to the documentation, which includes some examples but a tutorial can better capture the variety of input possible.

Goal for now is one for basic fitting with OFTI and one for basic fitting with MCMC. Add more as code capabilities expand.

Related to issue #37. We should have some approved tools for assessing convergence of MCMC walkers, since it is probably the trickiest part of MCMC orbit fitting. Similarly, we should expand the MCMC tutorial.

• function that returns the autocorrelation time of the chains
• documentation in the MCMC tutorial on how to use the above feature to assess convergence

Instead of defining separate variables for mass and parallax when they are fixed, let's just add them to system.sys_priors as floats instead of Prior objects

Set up a readthedocs site. Doesn't have to have anything.

I had some issues getting cython + Windows + Python 3 to work together, but supposedly according to the Cython docs it is possible: https://github.com/cython/cython/wiki/CythonExtensionsOnWindows

Regardless, it seems like having an option to disable cython featuers on install if the user is running into issues seems like a good idea. Basically add a flag to setup.py that disables it (E.g., python setup.py install --disable-cython)