Merge pull request #107 from California-Planet-Search/next-release

Next Release (v1.0.3)

README.md

@@ -1,3 +1,4 @@
+
 # RadVel
 
 General Toolkit for Modeling Radial Velocities.
@@ -16,29 +17,44 @@ Please cite the following DOI if you wish to make use of this software in any pu
 
 ## Documentation
 
-Documentation is available on ReadTheDocs.org: http://radvel.readthedocs.io
+Documentation is available [here](http://radvel.readthedocs.io/)
 
 ## Features
 
-- Object-oriented (i.e. models, likelihoods, priors, and posteriors are defined as objects)
-- Extensible (i.e. naturally define new likelihoods, priors, parameterizatoins)
-- Convenient API to fix/float parameters
-- Easily plugs in to the the suite of `scipy.optimize` routines for max-likelihood fitting 
-- Works with `emcee` MCMC
-- parameters are represented as dicts not arrays
-- Can handle data from multiple telescopes
-- Easily convert between different parameterizations
-- Computation of Kepler's equation (numerically intensive) written in C
 
-## Future Improvements...
+With RadVel you can
+
+
+- *Optimize*
+  - leverages the suite of minimizers in [scipy.optimize](https://docs.scipy.org/doc/scipy/reference/optimize.html)
+- *Run MCMC*
+  - leverages [emcee](http://dfm.io/emcee/) package for MCMC exploration of posterior
+- *Visualize*
+  - creates quicklook summary plots and statistics
+
+RadVel is
 
-- PERF: Optimizations for low eccentricity orbits
-- Add Gaussian Process (GP) functionality
+- *Flexible*
+  - fix/float parameters that are indexed as strings (emulates [lmfit](https://github.com/lmfit/lmfit-py/) API)
+  - convert between different parameterizations e.g. `e omega <-> sqrtecosw sqrtesinw`
+  - incorporate RVs from multiple telescopes
+- *Extensible* 
+  - Object-oriented programing makes adding new likelihoods, priors, etc. easy
+- *Scriptable*
+  - Code can be run through convenient Command-line Interface (CLI) 
+- *Fast*
+   - Kepler's equation solved in C (slower pure python solver available)
+   - MCMC is multi-threaded
+
+## Future Improvements
+
+- Gaussian Process (GP) functionality
 
 ## Tutorials 
 
 Follow examples in
 
 - `radvel/tests/SyntheticData.ipynb`
-- `radvel/tests/EPIC-2037_Fitting+MCMC.ipynb`
+- `radvel/tests/K2-24_Fitting+MCMC.ipynb`
 - `radvel/tests/164922_Fitting+MCMC.ipynb`
+

radvel/__init__.py

@@ -13,7 +13,7 @@
 __all__=['model', 'likelihood', 'posterior', 'mcmc', 'prior', 'utils',
          'fitting', 'report', 'cli', 'driver']
 
-__version__ = '1.0.2'
+__version__ = '1.0.3'
 
 MODULEDIR, filename = os.path.split(__file__)
 DATADIR = os.path.join(sys.prefix, 'radvel_example_data')

radvel/basis.py

@@ -8,7 +8,9 @@
                'per tc secosw sesinw logk',
                'per tc secosw sesinw k',
                'per tc ecosw esinw k',
-               'per tc e w k']
+               'per tc e w k',
+               'logper tc secosw sesinw k',
+               'logper tc secosw sesinw logk']
 
 
 def _print_valid_basis():
@@ -45,6 +47,7 @@ class Basis(object):
         'per tc secosw sesinw k'  \n
         'per tc ecosw esinw k'  \n
         'per tc e w k' \n
+        'logper tc secosw sesinw k'\n
         'logper tc secosw sesinw logk'
     """
     cps_params = 'per tp e w k'.split()
@@ -193,6 +196,21 @@ def _setpar(key, new_value):
                 w = np.arctan2(sesinw, secosw)
                 tp = timetrans_to_timeperi(tc, per, e, w)
 
+            if basis_name == 'logper tc secosw sesinw logk':
+                # pull out parameters
+                logper = _getpar('logper')
+                tc = _getpar('tc')
+                secosw = _getpar('secosw')
+                sesinw = _getpar('sesinw')
+                k = _getpar('logk')
+
+                # transform into CPS basis
+                per = np.exp(logper)
+                e = secosw ** 2 + sesinw ** 2
+                k = np.exp(k)
+                w = np.arctan2(sesinw, secosw)
+                tp = timetrans_to_timeperi(tc, per, e, w)
+
             if basis_name == 'per tc ecosw esinw k':
                 # pull out parameters
                 per = _getpar('per')
@@ -362,6 +380,32 @@ def _delpar(key):
                 self.name = newbasis
                 self.params = newbasis.split()
 
+            if newbasis == 'logper tc secosw sesinw logk':
+                per = _getpar('per')
+                e = _getpar('e')
+                w = _getpar('w')
+                k = _getpar('k')
+                try:
+                    tp = _getpar('tp')
+                except KeyError:
+                    tp = timetrans_to_timeperi(_getpar('tc'), per, e, w)
+                    _setpar('tp', tp)
+                _setpar('logper', np.log(per))
+                _setpar('secosw', np.sqrt(e)*np.cos(w))
+                _setpar('sesinw', np.sqrt(e)*np.sin(w))
+                _setpar('logk', np.log(k))
+                _setpar('tc', timeperi_to_timetrans(tp, per, e, w))
+
+                if not kwargs.get('keep', True):
+                    _delpar('per')
+                    _delpar('tp')
+                    _delpar('e')
+                    _delpar('w')
+                    _delpar('k')
+
+                self.name = newbasis
+                self.params = newbasis.split()
+
             if newbasis == 'per tc ecosw esinw k':
                 per = _getpar('per')
                 e = _getpar('e')

radvel/kepler.py

@@ -74,8 +74,8 @@ def kepler(inbigM, inecc):
     Earr = Marr + np.sign(np.sin(Marr)) * k * eccarr  # first guess at E
     # fiarr should go to zero when converges
     fiarr = ( Earr - eccarr * np.sin(Earr) - Marr)  
-    convd = np.abs(fiarr) > conv  # which indices have not converged
-    nd = np.sum(convd is True)  # number of converged elements
+    convd = np.where(np.abs(fiarr) > conv)[0]  # which indices have not converged
+    nd = len(convd)  # number of unconverged elements
     count = 0
 
     while nd > 0:  # while unconverged elements exist

radvel/plotting.py

@@ -480,7 +480,7 @@ def corner_plot_derived_pars(chains, P, saveplot=None):
     if 'planet_letters' in dir(P):
         planet_letters = P.planet_letters
     else:
-        planet_letters = ['a', 'b', 'c', 'd', 'e', 'f', 'g', 'h', 'i', 'j', 'k']
+        planet_letters = {1:'b', 2:'c', 3:'d', 4:'e', 5:'f', 6:'g', 7:'h', 8:'i', 9:'j', 10:'k'}
 
     # Determine which columns to include in corner plot
     labels = []

radvel/prior.py

@@ -1,10 +1,16 @@
+
 import numpy as np
+
 from radvel import model
+from radvel import orbit
+from radvel import utils
+
 
 class Prior(object):
     def __repr__(self):
         return "Generic Prior"
 
+
 class Gaussian(Prior):
     """Gaussian prior
 
@@ -20,17 +26,19 @@ def __init__(self, param, mu, sigma):
         self.mu = mu
         self.sigma = sigma
         self.param = param
+
     def __call__(self, params):
         x = params[self.param].value
-        return -0.5 * ( ( x - self.mu ) / self.sigma )**2
+        return -0.5 * ((x - self.mu) / self.sigma)**2
+
     def __repr__(self):
         s = "Gaussian prior on {}, mu={}, sigma={}".format(
             self.param, self. mu, self.sigma
             )
         return s
+
     def __str__(self):
         try:
-            d = {self.param: self.mu}
             tex = model.Parameters(9).tex_labels(param_list=[self.param])[self.param]
 
             s = "Gaussian prior on {}: ${} \\pm {}$ \\\\".format(tex, self. mu, self.sigma)
@@ -39,6 +47,7 @@ def __str__(self):
 
         return s
 
+
 class EccentricityPrior(Prior):
     """Physical eccentricities
 
@@ -55,10 +64,9 @@ class EccentricityPrior(Prior):
             each planet can have a different eccentricity upper limit.
     """
 
-
     def __repr__(self):
         msg = ""
-        for i,num_planet in enumerate(self.planet_list):
+        for i, num_planet in enumerate(self.planet_list):
             msg += "e{} constrained to be < {}\n".format(num_planet, self.upperlims[i])
 
         return msg[:-1]
@@ -67,7 +75,7 @@ def __str__(self):
         tex = model.Parameters(9, basis='per tc e w k').tex_labels()
 
         msg = ""
-        for i,num_planet in enumerate(self.planet_list):
+        for i, num_planet in enumerate(self.planet_list):
             par = "e{}".format(num_planet)
             label = tex[par]
             msg += "{} constrained to be $<{}$ \\\\\\\\\n".format(label, self.upperlims[i])
@@ -77,7 +85,7 @@ def __str__(self):
     def __init__(self, num_planets, upperlims=0.99):
 
         if type(num_planets) == int:
-            self.planet_list = range(1,num_planets+1)
+            self.planet_list = range(1, num_planets+1)
             npl = len(self.planet_list)
         else:
             self.planet_list = num_planets
@@ -92,28 +100,28 @@ def __init__(self, num_planets, upperlims=0.99):
 
     def __call__(self, params):
         def _getpar(key, num_planet):
-            return params['{}{}'.format(key,num_planet)].value
+            return params['{}{}'.format(key, num_planet)].value
 
         parnames = params.basis.name.split()
 
-        for i,num_planet in enumerate(self.planet_list):
+        for i, num_planet in enumerate(self.planet_list):
             if 'e' in parnames:
                 ecc = _getpar('e', num_planet)
             elif 'secosw' in parnames:
-                secosw = _getpar('secosw',num_planet)
-                sesinw = _getpar('sesinw',num_planet)
+                secosw = _getpar('secosw', num_planet)
+                sesinw = _getpar('sesinw', num_planet)
                 ecc = secosw**2 + sesinw**2 
             elif 'ecosw' in parnames:
-                ecosw = _getpar('ecosw',num_planet)
-                esinw = _getpar('esinw',num_planet)
+                ecosw = _getpar('ecosw', num_planet)
+                esinw = _getpar('esinw', num_planet)
                 ecc = np.sqrt(ecosw**2 + esinw**2)
 
-
             if ecc > self.upperlims[i] or ecc < 0.0:
                 return -np.inf
 
         return 0
-
+
+
 class PositiveKPrior(Prior):
     """K must be positive
 
@@ -128,6 +136,7 @@ class PositiveKPrior(Prior):
 
     def __repr__(self):
         return "K constrained to be > 0"
+
     def __str__(self):
         return "$K$ constrained to be $>0$"
 
@@ -136,18 +145,19 @@ def __init__(self, num_planets):
 
     def __call__(self, params):
         def _getpar(key, num_planet):
-            return params['{}{}'.format(key,num_planet)].value
+            return params['{}{}'.format(key, num_planet)].value
 
-        for num_planet in range(1,self.num_planets+1):
+        for num_planet in range(1, self.num_planets+1):
             try:
                 k = _getpar('k', num_planet)
             except KeyError:
-                k = np.exp(_getpar('logk',num_planet))
+                k = np.exp(_getpar('logk', num_planet))
 
             if k < 0.0:    
                 return -np.inf
         return 0
 
+
 class HardBounds(Prior):
     """Prior for hard boundaries
 
@@ -164,26 +174,79 @@ def __init__(self, param, minval, maxval):
         self.minval = minval
         self.maxval = maxval
         self.param = param
+
     def __call__(self, params):
         x = params[self.param].value
         if x < self.minval or x > self.maxval:
             return -np.inf
         else:
             return 0.0
+
     def __repr__(self):
         s = "Bounded prior on {}, min={}, max={}".format(
             self.param, self.minval, self.maxval
             )
         return s
+
     def __str__(self):
         try:
-            d = {self.param: self.minval}
             tex = model.Parameters(9).tex_labels(param_list=[self.param])[self.param]
 
             s = "Bounded prior: ${} < {} < {}$".format(self.minval,
-                                                       tex.replace('$',''),
+                                                       tex.replace('$', ''),
                                                        self.maxval)
         except KeyError:
             s = self.__repr__()
 
         return s
+
+
+class SecondaryEclipsePrior(Prior):
+    """Secondary eclipse prior
+
+    Implied prior on eccentricity and omega by specifying measured secondary eclipse time
+
+    Args:
+        planet_num (int): Number of planet with measured secondary eclipse
+        ts (float): Secondary eclipse midpoint time.
+            Should be in the same units as the timestamps of your data.
+        ts_err (float): Uncertainty on secondary eclipse time
+    """
+
+    def __repr__(self):
+        msg = ""
+        msg += "secondary eclipse constraint: {} +/- {}\n".format(self.ts, self.ts_err)
+
+        return msg[:-1]
+
+    def __str__(self):
+        msg = "secondary eclipse prior: ${} \pm {}$ \\\\\\\\\n".format(self.ts, self.ts_err)
+
+        return msg[:-5]
+
+    def __init__(self, planet_num, ts, ts_err):
+
+        self.planet_num = planet_num
+        self.ts = ts
+        self.ts_err = ts_err
+
+    def __call__(self, params):
+        def _getpar(key):
+            return cps_params['{}{}'.format(key, self.planet_num)].value
+
+        cps_params = params.basis.to_cps(params)
+
+        tp = _getpar('tp')
+        per = _getpar('per')
+        ecc = _getpar('e')
+        omega = _getpar('w')
+
+        ts = orbit.timeperi_to_timetrans(tp, per, ecc, omega, secondary=True)
+        ts_phase = utils.t_to_phase(cps_params, ts, self.planet_num)
+
+        pts = utils.t_to_phase(cps_params, self.ts, self.planet_num)
+        epts = self.ts_err / per
+
+        penalty = -0.5 * ((ts_phase - pts) / epts)**2
+
+        return penalty