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

Version 1.1.11

docs/requirements.txt

@@ -2,7 +2,7 @@ numpy>=1.11
 scipy>=0.17
 matplotlib>=1.5.3
 corner>=2.0
-cython>=0.22
+cython>=0.23
 astropy>=1.1.2
 emcee>=2.2.1
 pandas>=0.20.0

radvel/__init__.py

@@ -19,7 +19,7 @@ def _custom_warningfmt(msg, *a, **b):
 __all__=['model', 'likelihood', 'posterior', 'mcmc', 'prior', 'utils',
          'fitting', 'report', 'cli', 'driver', 'gp']
 
-__version__ = '1.1.10'
+__version__ = '1.1.11'
 
 MODULEDIR, filename = os.path.split(__file__)
 DATADIR = os.path.join(sys.prefix, 'radvel_example_data')

radvel/driver.py

@@ -20,6 +20,8 @@
 import radvel
 from radvel.plot import orbit_plots, mcmc_plots
 from astropy import constants as c
+from numpy import inf
+
 
 def plots(args):
     """

radvel/likelihood.py

@@ -120,7 +120,7 @@ def __init__(self, like_list):
         params = like0.params
         self.model = like0.model
         self.x = like0.x
-        self.y = like0.y - params[like0.gamma_param].value
+        self.y = like0.y #- params[like0.gamma_param].value
         self.yerr = like0.yerr
         self.telvec = like0.telvec
         self.extra_params = like0.extra_params

radvel/plot/orbit_plots.py

@@ -50,11 +50,14 @@ class MultipanelPlot(object):
             plots and residuals plots.
         set_xlim (list of float): limits to use for x-axes of the timeseries and residuals plots, in
             JD - `epoch`. Ex: [7000., 70005.]
+        text_size (int): set matplotlib.rcParams['font.size'] (default: 9)
+        legend_kwargs (dict): dict of options to pass to legend (plotted in top panel)
     """
     def __init__(self, post, saveplot=None, epoch=2450000, yscale_auto=False, yscale_sigma=3.0,
                 phase_nrows=None, phase_ncols=None, uparams=None, telfmts={},legend=True,
                 phase_limits=[], nobin=False, phasetext_size='large', rv_phase_space=0.08, 
-                figwidth=7.5, fit_linewidth=2.0, set_xlim=None):
+                figwidth=7.5, fit_linewidth=2.0, set_xlim=None, text_size=9,
+                legend_kwargs=dict(loc='best')):
 
         self.post = post
         self.saveplot=saveplot
@@ -75,6 +78,8 @@ def __init__(self, post, saveplot=None, epoch=2450000, yscale_auto=False, yscale
         self.figwidth = figwidth
         self.fit_linewidth = fit_linewidth
         self.set_xlim = set_xlim
+        self.legend_kwargs = legend_kwargs
+        rcParams['font.size'] = text_size
 
         if isinstance(self.post.likelihood, radvel.likelihood.CompositeLikelihood):
             self.like_list = self.post.likelihood.like_list
@@ -180,7 +185,7 @@ def plot_timeseries(self):
 
         # legend
         if self.legend:
-            ax.legend(numpoints=1, loc='best')
+            ax.legend(numpoints=1, **self.legend_kwargs)
 
         # years on upper axis
         axyrs = ax.twiny()
@@ -190,6 +195,7 @@ def plot_timeseries(self):
         axyrs.get_xaxis().get_major_formatter().set_useOffset(False)
         axyrs.set_xlim(*decimalyear)
         axyrs.set_xlabel('Year', fontweight='bold')
+        pl.locator_params(axis='x', nbins=5)
 
         if not self.yscale_auto: 
             scale = np.std(self.rawresid+self.rvmod)
@@ -329,13 +335,17 @@ def plot_phasefold(self, pltletter, pnum):
         )
 
 
-    def plot_multipanel(self, nophase=False):
+    def plot_multipanel(self, nophase=False, letter_labels=True):
         """
         Provision and plot an RV multipanel plot
 
         Args:
             nophase (bool, optional): if True, don't
                 include phase plots. Default: False.
+            letter_labels (bool, optional): if True, include 
+                letter labels on orbit and residual plots.
+                Default: True.
+
         Returns:
             tuple containing:
                 - current matplotlib Figure object
@@ -366,18 +376,20 @@ def plot_multipanel(self, nophase=False):
 
         pl.sca(ax_rv)
         self.plot_timeseries()
-        pltletter = ord('a')
-        plot.labelfig(pltletter)
-        pltletter += 1
+        if letter_labels:
+            pltletter = ord('a')
+            plot.labelfig(pltletter)
+            pltletter += 1
 
          # residuals
         ax_resid = pl.subplot(gs_rv[1, 0])
         self.ax_list += [ax_resid]
 
         pl.sca(ax_resid)
         self.plot_residuals()
-        plot.labelfig(pltletter)
-        pltletter += 1
+        if letter_labels:
+            plot.labelfig(pltletter)
+            pltletter += 1
 
 
         # phase-folded plots
@@ -435,15 +447,16 @@ def __init__(self, post, saveplot=None, epoch=2450000, yscale_auto=False, yscale
                 phase_nrows=None, phase_ncols=None, uparams=None, rv_phase_space=0.08, telfmts={},
                 legend=True,
                 phase_limits=[], nobin=False, phasetext_size='large',  figwidth=7.5, fit_linewidth=2.0,
-                set_xlim=None, subtract_gp_mean_model=False,
+                set_xlim=None, text_size=9, legend_kwargs=dict(loc='best'), subtract_gp_mean_model=False,
                 plot_likelihoods_separately=False, subtract_orbit_model=False):
 
         super(GPMultipanelPlot, self).__init__(
             post,saveplot=saveplot, epoch=epoch, yscale_auto=yscale_auto,
             yscale_sigma=yscale_sigma,phase_nrows=phase_nrows, phase_ncols=phase_ncols,
             uparams=uparams, rv_phase_space=rv_phase_space, telfmts=telfmts, legend=legend,
             phase_limits=phase_limits, nobin=nobin, phasetext_size=phasetext_size, 
-            figwidth=figwidth, fit_linewidth=fit_linewidth, set_xlim=set_xlim
+            figwidth=figwidth, fit_linewidth=fit_linewidth, set_xlim=set_xlim, text_size=text_size,
+            legend_kwargs=legend_kwargs
         )
 
         self.subtract_gp_mean_model=subtract_gp_mean_model
@@ -563,7 +576,7 @@ def plot_timeseries(self):
 
         # legend
         if self.legend:
-            ax.legend(numpoints=1, loc='best')
+            ax.legend(numpoints=1, **self.legend_kwargs)
 
         # years on upper axis
         axyrs = ax.twiny()
@@ -572,6 +585,7 @@ def plot_timeseries(self):
         axyrs.plot(decimalyear, decimalyear)
         axyrs.get_xaxis().get_major_formatter().set_useOffset(False)
         axyrs.set_xlim(*decimalyear)
+        pl.locator_params(axis='x', nbins=5)
         axyrs.set_xlabel('Year', fontweight='bold')
 
 
@@ -647,7 +661,7 @@ def plot_multipanel(self, nophase=False):
 
                 # legend
                 if self.legend and i==1:
-                    ax.legend(numpoints=1, loc='best')
+                    ax.legend(numpoints=1, **self.legend_kwargs)
 
                 # years on upper axis
                 if i==1:

radvel/posterior.py

@@ -77,7 +77,16 @@ def aic(self):
         n = len(self.likelihood.y)
         k = len(self.likelihood.get_vary_params())
         aic = - 2.0 * self.logprob() + 2.0 * k
-        _aicc = aic + (2.0 * k * (k + 1.0)) / (n - k - 1.0)
+        # Small sample correction
+        _aicc = aic
+        denom = (n - k - 1.0)
+        if denom > 0:
+            _aicc += (2.0 * k * (k + 1.0)) / denom 
+        else:
+            print("Warning: The number of free parameters is greater than or equal to") 
+            print("         the number of data points. The AICc comparison calculations") 
+            print("         will fail in this case.") 
+            _aicc = np.inf
         return _aicc
 
     def logprob_array(self, param_values_array):

radvel/report.py

@@ -24,15 +24,18 @@
     'jitter': 'm s$-1$',
     'logjit': '$\\ln{(\\rm m\\ s^{-1})}$',
     'jit': '$\\rm m\\ s^{-1}$',
-    'dvdt': 'm s$^{-1}$ day$^{-1}$',
-    'curv': 'm s$^{-1}$ day$^{-2}$',
+    'dvdt': 'm s$^{-1}$ d$^{-1}$',
+    'curv': 'm s$^{-1}$ d$^{-2}$',
     'gp_amp': 'm s$-1$',
     'gp_explength': 'days',
     'gp_per': 'days',
-    'gp_perlength': ''
+    'gp_perlength': '',
+    'mpsini': '$M_\earth$',
+    'rp': '$R_\earth$',
+    'rhop': 'g cm$^{-3}$',
 }
 
-class RadvelReport():
+class RadvelReport(object):
     """Radvel report
 
     Class to handle the creation of the radvel summary PDF
@@ -191,22 +194,32 @@ def _row(self, param, unit):
         med, errlow, errhigh = radvel.utils.sigfig(med, low, high)
 
         if min(errlow,errhigh) <= 1e-12:
-            med = maxlike = r"$\equiv$ %s" % round(self.quantiles[param][0.5],4)
+            med = maxlike = r"\equiv%s" % round(self.quantiles[param][0.5],4)
             errfmt = ''
         else:
-            if errhigh == errlow: errfmt = '$\pm %s$' % (errhigh)    
-            else: errfmt = '$^{+%s}_{-%s}$' % (errhigh,errlow)
+            if errhigh == errlow: errfmt = '\pm %s' % (errhigh)    
+            else: errfmt = '^{+%s}_{-%s}' % (errhigh,errlow)
 
-        row = "%s & %s %s & %s & %s" % (tex,med,errfmt,maxlike,unit)
+        row = "%s & $%s%s$ & $%s$ & %s" % (tex,med,errfmt,maxlike,unit)
         return row
 
-    def _data(self, basis):
+    def _data(self, basis, dontloop=False):
         """
         Helper function to output the rows in the parameter table
+
+        Args:
+            basis (str): name of Basis object (see basis.py) to be printed
+            dontloop (Bool): if True, don't loop over number of planets (useful for
+                printing out gamma, dvdt, jitter, curv)
         """
         suffixes = ['_'+j for j in self.report.post.likelihood.suffixes]
         rows = []
-        for n in range(1,self.report.planet.nplanets+1):
+
+        nloop = self.report.planet.nplanets+1
+        if dontloop:
+            nloop=2
+
+        for n in range(1,nloop):
             for p in basis.split(): # loop over variables
                 unit = units.get(p, '')
                 if unit == '':
@@ -299,7 +312,7 @@ def tab_params(self, name_in_title=False):
         kw = {}
         kw['fitting_basis_rows'] = self._data(self.fitting_basis)
         kw['print_basis_rows'] = self._data(print_basis)
-        kw['ep_rows'] = self._data(ep)
+        kw['ep_rows'] = self._data(ep, dontloop=True)
         kw['nlinks'] = len(self.report.chains)
         kw['time_base'] = self.report.post.likelihood.model.time_base
         if name_in_title:

requirements.txt

@@ -2,7 +2,7 @@ numpy>=1.11
 scipy>=0.17
 matplotlib>=1.5.3
 corner>=2.0
-cython>=0.22
+cython>=0.23
 astropy>=1.1.2
 emcee>=2.2.1
 pandas>=0.20.0