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

Version 1.1.8

docs/quickstartcli.rst

@@ -25,6 +25,9 @@ If you are running OSX, and want to perform Gaussian Process likelihood
 computations in parallel, you may need to perform some additional
 installation steps. See :ref:`OSX-multiprocessing`.
 
+If you wish to use the celerite GP kernels you will also need to install celerite.
+See the `celerite install instructions <http://celerite.readthedocs.io/en/stable/python/install/#using-pip>`_.
+
 
 Example Fit
 +++++++++++

example_planets/k2-131_celerite.py

@@ -47,16 +47,16 @@
 # First Celerite Term. AC>=BD must be true to ensure 
 # positive-definiteness. (A == params['1_logA'].value, etc.)
 params['1_logA'] = radvel.Parameter(value=np.log(26.))
-params['1_logB'] = radvel.Parameter(value=np.log(.0005)) 
+params['1_logB'] = radvel.Parameter(value=np.log(.0005))
 params['1_logC'] = radvel.Parameter(value=np.log(.5))
 params['1_logD'] = radvel.Parameter(value=np.log(.0005))
 
 # Second Celerite Term (real). Setting vary=False for 1_logB and 
 # 1_logD ensures that this term remains real throughout the fitting process.
 params['2_logA'] = radvel.Parameter(value=np.log(26.)) 
-params['2_logB'] = radvel.Parameter(value=np.log(0.),vary=False)
+params['2_logB'] = radvel.Parameter(value=np.log(0.005), vary=False)
 params['2_logC'] = radvel.Parameter(value=np.log(0.1)) 
-params['2_logD'] = radvel.Parameter(value=np.log(0.),vary=False)
+params['2_logD'] = radvel.Parameter(value=np.log(0.005), vary=False)
 
 hnames = {}
 for tel in instnames:
@@ -84,6 +84,7 @@ def initialize_instparams(tel_suffix):
           radvel.prior.Gaussian('tc1', Tc, Tc_unc),
           radvel.prior.Jeffreys('k1', 0.01, 10.),
           radvel.prior.Jeffreys('jit_pfs', 0.01, 10.),
-          radvel.prior.Jeffreys('jit_harps-n', 0.01,10.)
+          radvel.prior.Jeffreys('jit_harps-n', 0.01,10.),
+          radvel.prior.Jeffreys('1_logA', 0.005, 75.0)
           ]
 

radvel/__init__.py

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

radvel/driver.py

@@ -282,8 +282,12 @@ def tables(args):
                 P, post, chains, compstats=compstats
             )
             tex = tabletex.tab_comparison()
-        else:
-            tex = tabletex.tab_comparison()
+        elif tabtype == 'rv':
+            tex = tabletex.tab_rv()
+        elif tabtype == 'params':
+            tex = tabletex.tab_params()
+        elif tabtype == 'priors':
+            tex = tabletex.tab_prior_summary()
 
         saveto = os.path.join(
             args.outputdir, '{}_{}_.tex'.format(conf_base,tabtype)

radvel/fitting.py

@@ -4,37 +4,37 @@
 import collections
 
 
-def maxlike_fitting(post, verbose=True):
+def maxlike_fitting(post, verbose=True, method='Powell'):
     """Maximum Likelihood Fitting
 
     Perform a maximum likelihood fit.
 
     Args:
         post (radvel.Posterior): Posterior object with initial guesses
         verbose (bool [optional]): Print messages and fitted values?
+        method (string [optional]): Minimization method. See documentation for `scipy.optimize.minimize` for available
+            options.
 
     Returns: 
         radvel.Posterior : Posterior object with parameters
         updated to their maximum likelihood values
 
     """
 
-    post0 = copy.copy(post)
     if verbose:
-        print("Initial loglikelihood = %f" % post0.logprob())
+        print("Initial loglikelihood = %f" % post.logprob())
         print("Performing maximum likelihood fit...")
     res = scipy.optimize.minimize(
-        post.neglogprob_array, post.get_vary_params(), method='Nelder-Mead',
+        post.neglogprob_array, post.get_vary_params(), method=method,
         options=dict(xatol=1e-8, maxiter=200, maxfev=100000)
     )
-    synthpost = copy.copy(post)
     synthparams = post.params.basis.to_synth(post.params, noVary = True) # setting "noVary" assigns each new parameter a vary attribute
-    synthpost.params.update(synthparams)                                 # of '', for printing purposes
+    post.params.update(synthparams)                                 # of '', for printing purposes
 
     if verbose:
         print("Final loglikelihood = %f" % post.logprob())
         print("Best-fit parameters:")
-        print(synthpost)
+        print(post)
 
     return post
 
@@ -56,7 +56,7 @@ def model_comp(post, verbose=False):
             dictionary is a tuple with the statistic value as the first 
             element and a description of that statistic in the second element.
     """
-    ipost = copy.deepcopy(post)
+    iparams = copy.deepcopy(post.params)
 
     num_planets = post.likelihood.model.num_planets
 
@@ -75,8 +75,8 @@ def model_comp(post, verbose=False):
                         post.params[par].value = 0.0
                     post.params[par].vary = False
                 else:
-                    post.params[par].value = ipost.params[par].value
-                    post.params[par].vary = ipost.params[par].vary
+                    post.params[par].value = iparams[par].value
+                    post.params[par].vary = iparams[par].vary
             except (ValueError, KeyError):
                 pass
 

radvel/gp.py

@@ -4,8 +4,6 @@
 from scipy import spatial
 import abc
 import numpy as np
-import celerite
-from celerite.solver import CholeskySolver
 
 # implemented kernels & examples of possible names for their associated hyperparameters
 KERNELS = {
@@ -20,6 +18,23 @@
 else:
     ABC = abc.ABC
 
+
+# celerite is an optional dependency
+def _try_celerite():
+    try:
+        import celerite
+        from celerite.solver import CholeskySolver
+        return True
+    except ImportError:
+        print("WARNING: celerite not installed. GP kernals using celerite will not work.")
+        print("Try installing celerite using 'pip install celerite'")
+        return False
+
+_has_celerite = _try_celerite()
+if _has_celerite:
+    import celerite
+    from celerite.solver import CholeskySolver
+
 class Kernel(ABC):
     """
     Abstract base class to store kernel info and compute covariance matrix.

radvel/likelihood.py

@@ -3,8 +3,10 @@
 from radvel import gp
 from scipy.linalg import cho_factor, cho_solve
 from scipy import matrix
-import celerite
 
+_has_celerite = gp._try_celerite()
+if _has_celerite:
+    import celerite
 
 class Likelihood(object):
     """
@@ -487,6 +489,7 @@ def logprob(self):
         else:
            print("WARNING: CeleriteKernel has encountered a non-positive-definite"
                 +" kernel. Ensure that ac>=bd for all complex kernel terms.")
+
            return -np.inf
 
     def predict(self,xpred):

radvel/mcmc.py

@@ -1,6 +1,5 @@
 import sys
 import time
-import copy
 
 import multiprocessing as mp
 
@@ -160,17 +159,16 @@ def mcmc(post, nwalkers=50, nrun=10000, ensembles=8, checkinterval=50, burnGR=1.
             pscale = post.params[par].mcmcscale
         pscales.append(pscale)
     pscales = np.array(pscales)
-    
+
     statevars.samplers = []
     statevars.initial_positions = []
     for e in range(ensembles):
-        pcopy = copy.deepcopy(post)
-        pi = pcopy.get_vary_params()
+        pi = post.get_vary_params()
         p0 = np.vstack([pi]*statevars.nwalkers)
         p0 += [np.random.rand(statevars.ndim)*pscales for i in range(statevars.nwalkers)]
         statevars.initial_positions.append(p0)
         statevars.samplers.append(emcee.EnsembleSampler( 
-            statevars.nwalkers, statevars.ndim, pcopy.logprob_array, threads=1))
+            statevars.nwalkers, statevars.ndim, post.logprob_array, threads=1))
 
     num_run = int(np.round(nrun / checkinterval))
     statevars.totsteps = nrun*statevars.nwalkers*statevars.ensembles
@@ -256,7 +254,7 @@ def mcmc(post, nwalkers=50, nrun=10000, ensembles=8, checkinterval=50, burnGR=1.
 
     df = pd.DataFrame(
         statevars.tchains.reshape(statevars.ndim,statevars.tchains.shape[1]*statevars.tchains.shape[2]).transpose(),
-        columns=pcopy.list_vary_params())
+        columns=post.list_vary_params())
     df['lnprobability'] = np.hstack(statevars.lnprob)
 
     df = df.iloc[::thin]

radvel/prior.py

@@ -379,10 +379,7 @@ class NumericalPrior(Prior):
 
     def __init__(self, param_list, values, bw_method=None):
         self.param_list = param_list
-        self.values = values
-        self.bw_method = bw_method
-
-        self.pdf_estimate = gaussian_kde(self.values, bw_method=self.bw_method)
+        self.pdf_estimate = gaussian_kde(values, bw_method=bw_method)
 
     def __call__(self, params):
         x = []

radvel/report.py

@@ -52,22 +52,20 @@ def __init__(self, planet, post, chains, compstats=None):
         self.starname_tex = planet.starname.replace('_', '\\_')
         self.runname = self.starname_tex
 
-        printpost = copy.deepcopy(post)
-        printpost.params = printpost.params.basis.to_synth(printpost.params)
-        printpost.params = printpost.params.basis.from_synth(
-            printpost.params, print_basis
+        post.params = post.params.basis.to_synth(post.params)
+        post.params = post.params.basis.from_synth(
+            post.params, print_basis
         )
-        self.latex_dict = printpost.params.tex_labels()
+        self.latex_dict = post.params.tex_labels()
 
-        printchains = copy.copy(chains)
         for p in post.params.keys():
             if p not in chains.columns:
                 chains[p] = post.params[p].value
 
-        self.chains = printpost.params.basis.to_synth(
+        self.chains = post.params.basis.to_synth(
             chains, basis_name=planet.fitting_basis
         )
-        self.chains = printpost.params.basis.from_synth(
+        self.chains = post.params.basis.from_synth(
             self.chains, print_basis
         )
         self.quantiles = self.chains.quantile([0.159, 0.5, 0.841])

radvel/tests/test_api.py

@@ -255,4 +255,4 @@ def test_kepler():
 
 
 if __name__ == '__main__':
-    test_k2(setupfn='/Users/petigura/code/radvel/example_planets/epic203771098.py')
+    test_celerite()

radvel/utils.py

@@ -421,12 +421,8 @@ def semi_amplitude(Msini, P, Mtotal, e, Msini_units='jupiter'):
         P (float): Orbital period [days]
         Mtotal (float): Mass of star + mass of planet [Msun]
         e (float): eccentricity
-<<<<<<< HEAD
         Msini_units (Optional[str]): Units of Msini {'earth','jupiter'}
             default: 'jupiter'
-=======
-        Msini_units (string): (optional) Units of Msini {'earth','jupiter'} (default = 'jupiter')
->>>>>>> c0b22c0144d98ac2d25d2d5b135eee4ed08fbd42
 
     Returns:
         Doppler semi-amplitude [m/s]
@@ -517,17 +513,10 @@ def density(mass, radius, MR_units='earth'):
     Args:
         mass (float): mass [MR_units]
         radius (float): radius [MR_units]
-<<<<<<< HEAD
         MR_units (string): (optional) units of mass and radius. Must be 'earth', or 'jupiter' (default 'earth').
 
-        float: density in g/cc
-=======
-        MR_units (Optional[str]): Units of Msini {'earth','jupiter'} 
-            default: 'earth'
-
     Returns:
-        float: density [g/cc
->>>>>>> c0b22c0144d98ac2d25d2d5b135eee4ed08fbd42
+        float: density in g/cc
     """
 
     mass = np.array(mass)

requirements.txt

@@ -12,4 +12,3 @@ ipykernel
 pybind11
 celerite>=0.3.0
 jinja2
-tornado==4.5.3