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

Version 1.1.5

.readthedocs.yml

@@ -0,0 +1,5 @@
+requirements_file:
+    docs/requirements.txt
+python:
+    setup_py_install: true
+    version: 3.5

.travis.yml

@@ -26,7 +26,7 @@ before_install:
 install:
     - conda install --yes pip numpy scipy matplotlib cython pandas nose nbformat nbconvert jupyter_client ipykernel
     - if [ ${py} == "2" ]; then conda install --yes configparser==3.5.0b2; fi
-    - pip install coveralls
+    - pip install coveralls celerite
     - pip install .
     - python setup.py build_ext -i
 

docs/conf.py

@@ -18,7 +18,7 @@
 import matplotlib
 matplotlib.use('agg')
 
-autodoc_mock_imports = ['_tkinter', 'pandas', 'scipy.optimize']
+autodoc_mock_imports = ['_tkinter', 'pandas', 'scipy.optimize', 'celerite', 'celerite.solver']
 for mod_name in autodoc_mock_imports:
     sys.modules[mod_name] = mock.Mock()
 

docs/requirements.txt

@@ -0,0 +1,12 @@
+numpy>=1.11
+scipy>=0.17
+matplotlib>=1.5.3
+corner>=2.0
+cython>=0.22
+astropy>=1.1.2
+emcee>=2.2.1
+pandas>=0.20.0
+nbsphinx
+jupyter_client
+ipykernel
+pybind11

example_planets/k2-131.py

@@ -3,12 +3,15 @@
 import os
 import radvel
 
-# IMPORTANT: This config file uses Gaussian Process (GP) regression to model stellar activity.
-#            Don't start here. Make sure you understand how RadVel works using the basic configuration
-#            files (radvel/example_planets/epic203771098.py and radvel/example_planets/HD164922.py)
-#            before coming back to this configuration file.
+"""
+This config file uses Gaussian Process (GP) regression to model stellar activity.
 
-# Data from Dai+ 2017
+Don't start here. Make sure you understand how RadVel works using the basic configuration
+files (radvel/example_planets/epic203771098.py and radvel/example_planets/HD164922.py)
+before coming back to this configuration file.
+"""
+
+# Data from Dai et al. (2017)
 instnames = ['harps-n','pfs'] 
 data = pd.read_csv(os.path.join(radvel.DATADIR,'k2-131.txt'), sep=' ')
 t = np.array(data['time']) 
@@ -24,17 +27,19 @@
 nplanets = 1
 fitting_basis = 'per tc secosw sesinw k'
 
-# Numbers for priors used in Dai+ 2017
-gp_explength_mean = 9.5*np.sqrt(2.) # sqrt(2)*tau in Dai+ 2017 [days]
+# Numbers for priors used in Dai et al. (2017)
+gp_explength_mean = 9.5*np.sqrt(2.) # sqrt(2)*tau in Dai et al. (2017) [days]
 gp_explength_unc = 1.0*np.sqrt(2.)
-gp_perlength_mean = np.sqrt(1./(2*3.32)) # sqrt(1/2*gamma) in Dai+ 2017
+gp_perlength_mean = np.sqrt(1./(2*3.32)) # sqrt(1/2*gamma) in Dai et al. (2017)
 gp_perlength_unc = 0.04
-# NOTE: this prior isn't equivalent to the one Dai+ 2017 use. However,
-#    our formulation of the quasi-periodic kernel explicitly keeps the covariance
-#    matrix postitive semi-definite, so we use this instead. The orbit model
-#    results aren't affected.
-
-gp_per_mean = 9.64 # T_bar in Dai+ 2017 [days]
+""" 
+NOTE: this prior isn't equivalent to the one Dai et al. (2017) use. However,
+our formulation of the quasi-periodic kernel explicitly keeps the covariance
+matrix postitive semi-definite, so we use this instead. The orbit model
+results aren't affected.
+"""
+
+gp_per_mean = 9.64 # T_bar in Dai et al. (2017) [days]
 gp_per_unc = 0.12
 Porb = 0.3693038 # [days]
 Porb_unc = 0.0000091 
@@ -54,32 +59,37 @@
 time_base = np.median(t)
 
 # Define GP hyperparameters as Parameter objects.
-params['gp_amp'] = radvel.Parameter(value=26.0)
+params['gp_amp_j'] = radvel.Parameter(value=26.0)
+params['gp_amp_h'] = radvel.Parameter(value=26.0)
 params['gp_explength'] = radvel.Parameter(value=gp_explength_mean) 
 params['gp_per'] = radvel.Parameter(value=gp_per_mean) 
 params['gp_perlength'] = radvel.Parameter(value=gp_perlength_mean)
 
 
-# Define a dictionary, 'hnames', specifying the names
-# of the GP hyperparameters corresponding to a particular
-# data set. Use the strings in 'instnames' to tell radvel
-# which data set you're talking about.
+"""
+Define a dictionary, 'hnames', specifying the names
+of the GP hyperparameters corresponding to a particular
+data set. Use the strings in 'instnames' to tell radvel
+which data set you're talking about.
+"""
 hnames = {
-  'harps-n': ['gp_amp', # GP variability amplitude
+  'harps-n': ['gp_amp_h', # GP variability amplitude
               'gp_per', # GP variability period
               'gp_explength', # GP non-periodic characteristic length
               'gp_perlength'], # GP periodic characteristic length
-  'pfs': ['gp_amp', 
+  'pfs': ['gp_amp_j', 
           'gp_per', 
           'gp_explength', 
           'gp_perlength']
 }
 
 kernel_name = {'harps-n':"QuasiPer", 
                'pfs':"QuasiPer"}
-# If all kernels are quasi-periodic, you don't need to include the
-# 'kernel_name' lines. I included it to show you how to use different kernel
-# types in different likelihoods.
+"""
+NOTE: If all kernels are quasi-periodic, you don't need to include the
+'kernel_name' lines. I included it to show you how to use different kernel
+types in different likelihoods.
+"""
 
 jit_guesses = {'harps-n':2.0, 'pfs':5.3}
 
@@ -94,11 +104,12 @@ def initialize_instparams(tel_suffix):
 for tel in instnames:
     initialize_instparams(tel)
 
-# Add in priors (Dai+ Section 7.2.1)
+# Add in priors (Dai et al. 2017 Section 7.2.1)
 priors = [radvel.prior.Gaussian('per1', Porb, Porb_unc),
           radvel.prior.Gaussian('tc1', Tc, Tc_unc),
           radvel.prior.Jeffreys('k1', 0.01, 10.), # min and max for Jeffrey's priors estimated by Sarah
-          radvel.prior.Jeffreys('gp_amp', 0.01, 100.),
+          radvel.prior.Jeffreys('gp_amp_h', 0.01, 100.),
+          radvel.prior.Jeffreys('gp_amp_j', 0.01, 100.),
           radvel.prior.Jeffreys('jit_pfs', 0.01, 10.),
           radvel.prior.Jeffreys('jit_harps-n', 0.01,10.),
           radvel.prior.Gaussian('gp_explength', gp_explength_mean, gp_explength_unc),

radvel/__init__.py

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

radvel/gp.py

@@ -8,10 +8,12 @@
 from celerite.solver import CholeskySolver
 
 # implemented kernels & examples of possible names for their associated hyperparameters
-KERNELS = {"SqExp": ['gp_length','gp_amp'],
-           "Per": ['gp_per','gp_length','gp_amp'],
-           "QuasiPer": ['gp_per','gp_perlength','gp_explength','gp_amp'],
-           "Celerite": ['1_logA','1_logB','1_logC','1_logD']}
+KERNELS = {
+    "SqExp": ['gp_length','gp_amp'],
+    "Per": ['gp_per','gp_length','gp_amp'],
+    "QuasiPer": ['gp_per','gp_perlength','gp_explength','gp_amp'],
+    "Celerite": ['1_logA','1_logB','1_logC','1_logD']
+}
 
 if sys.version_info[0] < 3:
     ABC = abc.ABCMeta('ABC', (), {})
@@ -21,7 +23,12 @@
 class Kernel(ABC):
     """
     Abstract base class to store kernel info and compute covariance matrix.
-    All kernel objects inherit from this class
+    All kernel objects inherit from this class.
+
+    Note:
+        To implement your own kernel, create a class that inherits
+        from this class. It should have hyperparameters that follow
+        the name scheme 'gp_NAME_SUFFIX'.
 
     """
 

radvel/likelihood.py

@@ -307,8 +307,10 @@ def update_kernel_params(self):
         """ Update the Kernel object with new values of the hyperparameters 
         """
         for key in self.list_vary_params():
-            if key in self.hnames:  
-                self.kernel.hparams[key].value = self.params[key].value
+            if key in self.hnames: 
+                hparams_key = key.split('_')
+                hparams_key = hparams_key[0] + '_' + hparams_key[1]
+                self.kernel.hparams[hparams_key].value = self.params[key].value
 
     def _resids(self):
         """Residuals for internal GP calculations

radvel/prior.py

@@ -1,5 +1,6 @@
 
 import numpy as np
+from scipy.stats import gaussian_kde
 
 from radvel import model
 from radvel import orbit
@@ -343,3 +344,123 @@ def __str__(self):
             s = self.__repr__()
 
         return s
+
+class NumericalPrior(Prior):
+    """Prior defined by an input array of values
+
+    Wrapper for scipy.stats.gaussian_kde.
+
+    This prior uses Gaussian Kernel Density Estimation to
+    estimate the probability density function from which
+    a set of values are randomly drawn.
+
+    Useful for defining a prior given a posterior obtained
+    from a complementary fitting process. For example, you
+    might use transit data to obtain constraints on secosw and
+    sesinw, then use the posterior on secosw as a prior for
+    a RadVel fit.
+
+    Args:
+        param_list (list of str): list of parameter label(s). 
+        values (numpy array of float): values of ``param`` you
+            wish to use to define this prior. For example, this 
+            might be a posterior array of values of secosw 
+            derived from transit data. In case of univariate data 
+            this is a 1-D array, otherwise a 2-D array with shape 
+            (# of elements in param_list, # of data points).
+        bw_method (str, scalar, or callable [optional]): see 
+            scipy.stats.gaussian_kde
+
+    Note: the larger the input array of values, the longer it will
+    take for calls to this prior to be evaluated. Consider thinning
+    large input arrays to speed up performance.
+
+    """
+
+    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)
+
+    def __call__(self, params):
+        x = []
+        for param in self.param_list:
+            x.append(params[param].value)
+        return self.pdf_estimate(x)
+
+    def __repr__(self):
+        s = "Numerical prior on {}".format(
+            self.param_list
+            )
+        return s
+    def __str__(self):
+        try:
+            tex = model.Parameters(9).tex_labels(param_list=self.param_list)
+            t=[tex[key] for key in tex.keys()]
+            if len(self.param_list) == 1:
+                str2print = '{0}'.format(*t)
+            elif len(self.param_list) == 2:
+                str2print = '{} and {}'.format(*t)
+            else:
+                str2print = ''
+                for el in np.arange(len(self.param_list) - 1):
+                    str2print += '{}, '.format(t[el])
+                str2print += 'and {}'.format(t[el+1])
+            s = "Numerical prior on " + str2print + \
+                ", defined using Gaussian kernel density estimation."
+        except KeyError:
+            s = self.__repr__()
+
+        return s
+
+
+class UserDefinedPrior(Prior):
+    """Interface for user to define a prior 
+       with an arbitrary functional form. 
+
+    Args:
+        param_list (list of str): list of parameter label(s). 
+        func (function): a Python function that takes in  a list
+            of values (ordered as in ``param_list``), and returns
+            the corresponding log-value of a pdf. 
+        tex_rep (str): TeX-readable string representation of
+            this prior, to be passed into radvel report and 
+            plotting code.
+
+    Example:
+        >>> def myPriorFunc(inp_list):
+        ...     if inp_array == [0.]:
+        ...         return 0.
+        ...     else:
+        ...         return -np.inf
+        >>> myTexString = 'Delta Function Prior on $\sqrt{e}$'
+        >>> myPrior = radvel.prior.UserDefinedPrior(['se'], myPriorFunc, myTexString)
+
+    Note:
+        ``func`` must be properly normalized; i.e. integrating over the
+        entire parameter space must give probability 1. 
+    """
+
+    def __init__(self, param_list, func, tex_rep):
+        self.param_list = param_list
+        self.func = func
+        self.tex_rep = tex_rep
+
+    def __call__(self, params):
+        x = []
+        for param in self.param_list:
+            x.append(params[param].value)
+        return self.func(x)
+
+    def __repr__(self):
+        s = "User-defined prior on {}".format(
+            self.param_list
+            )
+        return s
+
+    def __str__(self):
+        s = self.tex_rep
+        return s
+