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

Version 1.3.3

radvel/__init__.py

@@ -25,7 +25,7 @@ def _custom_warningfmt(msg, *a, **b):
 __all__ = ['model', 'likelihood', 'posterior', 'mcmc', 'prior', 'utils',
          'fitting', 'report', 'cli', 'driver', 'gp']
 
-__version__ = '1.3.2'
+__version__ = '1.3.3'
 __spec__ = __name__
 __package__ = __path__[0]
 

radvel/driver.py

@@ -53,6 +53,7 @@ def plots(args):
 
         if ptype == 'rv':
             args.plotkw['uparams'] = post.uparams
+            args.plotkw['status'] = status
             saveto = os.path.join(
                 args.outputdir,conf_base+'_rv_multipanel.pdf'
             )
@@ -521,12 +522,25 @@ def _get_colname(key):
         except (AttributeError, KeyError):
             pass
 
-    print("Derived parameters:", outcols)
+    # Get quantiles and update posterior object to median
+    # values returned by MCMC chains
+    quantiles = chains.quantile([0.159, 0.5, 0.841])
+    csvfn = os.path.join(args.outputdir, conf_base+'_derived_quantiles.csv')
+    quantiles.to_csv(csvfn, columns=outcols)
+
+    # saved derived paramters to posterior file
+    postfile = os.path.join(args.outputdir,
+                            '{}_post_obj.pkl'.format(conf_base))
+    post.derived = quantiles[outcols]
+    post.writeto(postfile)
+    savestate['quantfile'] = os.path.relpath(csvfn)
 
     csvfn = os.path.join(args.outputdir, conf_base+'_derived.csv.tar.bz2')
     chains.to_csv(csvfn, columns=outcols, compression='bz2')
     savestate['chainfile'] = os.path.relpath(csvfn)
 
+    print("Derived parameters:", outcols)
+
     save_status(statfile, 'derive', savestate)
 
 

radvel/fitting.py

@@ -102,11 +102,11 @@ def model_comp(post, params=[], mc_list=[], verbose=False):
             print("BIC = %4.2f" % fitpost.likelihood.bic())
             print("AIC = %4.2f" % fitpost.likelihood.aic())
 
-        comparison_parameters = ['Free Params', '$N_{\\rm free}$', '$N_{\\rm data}$',
-                                 'RMS', '$\\ln{\\mathcal{L}}$', 'BIC', 'AICc']
+        comparison_parameters = ['Free Params', r'$N_{\rm free}$', r'$N_{\rm data}$',
+                                 'RMS', r'$\ln{\mathcal{L}}$', 'BIC', 'AICc']
         pdict = collections.OrderedDict.fromkeys(comparison_parameters)
-        pdict['$N_{\\rm data}$'] = (ndata, 'number of measurements')
-        pdict['$N_{\\rm free}$'] = (nfree, 'number of free parameters')
+        pdict[r'$N_{\rm data}$'] = (ndata, 'number of measurements')
+        pdict[r'$N_{\rm free}$'] = (nfree, 'number of free parameters')
         pdict['RMS'] = (
             np.round(np.std(fitpost.likelihood.residuals()), 2),
             'RMS of residuals in m s$^{-1}$'
@@ -115,7 +115,7 @@ def model_comp(post, params=[], mc_list=[], verbose=False):
         # pdict['$\\chi^{2}_{\\nu}$'] = (
         #    np.round(chi_red,2), "jitter fixed"
         # )
-        pdict['$\\ln{\\mathcal{L}}$'] = (
+        pdict[r'$\ln{\mathcal{L}}$'] = (
             np.round(fitpost.logprob(), 2), "natural log of the likelihood"
         )
         pdict['BIC'] = (

radvel/plot/mcmc_plots.py

@@ -186,6 +186,7 @@ def __init__(self, chains, P, saveplot=None):
         # Determine which columns to include in corner plot
         self.labels = []
         self.texlabels = []
+        self.units = []
         for i in np.arange(1, P.nplanets + 1, 1):
             letter = planet_letters[i]
 
@@ -211,15 +212,14 @@ def __init__(self, chains, P, saveplot=None):
                     if np.median(self.chains[label]) > 100:
                         unit = "M$_{\\odot}$"
                         self.chains[label] *= 0.000954265748
-
-                    tl += " [%s]" % unit
                 elif key == 'rhop':
-                    tl += " [g cm$^{-3}$]"
+                    unit = " g cm$^{-3}$"
                 elif key == 'a':
-                    tl += " [AU]"
+                    unit = " AU"
                 else:
-                    tl += " "
+                    unit = " "
 
+                self.units.append(unit)
                 self.labels.append(label)
                 self.texlabels.append(tl)
 
@@ -231,8 +231,13 @@ def plot(self):
         f = rcParams['font.size']
         rcParams['font.size'] = 12
 
+        plot_labels = []
+        for t, u in zip(self.texlabels, self.units):
+            label = '{} [{}]'.format(t, u)
+            plot_labels.append(label)
+
         _ = corner.corner(
-            self.chains[self.labels], labels=self.texlabels, label_kwargs={"fontsize": 14}, 
+            self.chains[self.labels], labels=plot_labels, label_kwargs={"fontsize": 14},
             plot_datapoints=False, bins=30, quantiles=[0.16, 0.50, 0.84],
             show_titles=True, title_kwargs={"fontsize": 14}, smooth=True
         )

radvel/plot/orbit_plots.py

@@ -1,11 +1,13 @@
 import numpy as np
+import pandas as pd
 from matplotlib import rcParams, gridspec
 from matplotlib import pyplot as pl
 from matplotlib.ticker import MaxNLocator
 from astropy.time import Time
 
 import radvel
 from radvel import plot
+from radvel.plot import mcmc_plots
 from radvel.utils import t_to_phase, fastbin, sigfig
 
 
@@ -55,12 +57,13 @@ class MultipanelPlot(object):
         highlight_last (bool): make the most recent measurement much larger in all panels
         show_rms (bool): show RMS of the residuals by instrument in the legend
         legend_kwargs (dict): dict of options to pass to legend (plotted in top panel)
+        status (ConfigParser): (optional) result of radvel.driver.load_status on the .stat status file
     """
     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, text_size=9, highlight_last=False,
-                 show_rms=False, legend_kwargs=dict(loc='best')):
+                 show_rms=False, legend_kwargs=dict(loc='best'), status=None):
 
         self.post = post
         self.saveplot = saveplot
@@ -86,6 +89,9 @@ def __init__(self, post, saveplot=None, epoch=2450000, yscale_auto=False, yscale
         self.legend_kwargs = legend_kwargs
         rcParams['font.size'] = text_size
 
+        if status is not None:
+            self.status = status
+
         if isinstance(self.post.likelihood, radvel.likelihood.CompositeLikelihood):
             self.like_list = self.post.likelihood.like_list
         else:
@@ -326,7 +332,7 @@ def plot_phasefold(self, pltletter, pnum):
             val = self.post.params["%s%d" % (print_params[l], pnum)].value
 
             if self.uparams is None:
-                _anotext = '$\\mathregular{%s}$ = %4.2f %s' % (labels[l].replace("$", ""), val, units[p])
+                _anotext = r'$\mathregular{%s}$ = %4.2f %s' % (labels[l].replace("$", ""), val, units[p])
             else:
                 if hasattr(self.post, 'medparams'):
                     val = self.post.medparams["%s%d" % (print_params[l], pnum)]
@@ -338,12 +344,40 @@ def plot_phasefold(self, pltletter, pnum):
                 err = self.uparams["%s%d" % (print_params[l], pnum)]
                 if err > 1e-15:
                     val, err, errlow = sigfig(val, err)
-                    _anotext = '$\\mathregular{%s}$ = %s $\\mathregular{\\pm}$ %s %s' \
+                    _anotext = r'$\mathregular{%s}$ = %s $\mathregular{\pm}$ %s %s' \
                                % (labels[l].replace("$", ""), val, err, units[p])
                 else:
-                    _anotext = '$\\mathregular{%s}$ = %4.2f %s' % (labels[l].replace("$", ""), val, units[p])
-
-            anotext += [_anotext] 
+                    _anotext = r'$\mathregular{%s}$ = %4.2f %s' % (labels[l].replace("$", ""), val, units[p])
+
+            anotext += [_anotext]
+
+        if hasattr(self.post, 'derived'):
+            chains = pd.read_csv(self.status['derive']['chainfile'])
+            self.post.nplanets = self.num_planets
+            dp = mcmc_plots.DerivedPlot(chains, self.post)
+            labels = dp.labels
+            texlabels = dp.texlabels
+            units = dp.units
+            derived_params = ['mpsini']
+            for l, par in enumerate(derived_params):
+                par_label = par + str(pnum)
+                if par_label in self.post.derived.columns:
+                    val = self.post.derived["%s%d" % (derived_params[l], pnum)].loc[0.500]
+                    low = self.post.derived["%s%d" % (derived_params[l], pnum)].loc[0.159]
+                    high = self.post.derived["%s%d" % (derived_params[l], pnum)].loc[0.841]
+                    err_low = val - low
+                    err_high = high - val
+                    err = np.mean([err_low, err_high])
+                    err = radvel.utils.round_sig(err)
+                    index = np.where(np.array(labels) == par_label)[0][0]
+                    if err > 1e-15:
+                        val, err, errlow = sigfig(val, err)
+                        _anotext = r'$\mathregular{%s}$ = %s $\mathregular{\pm}$ %s %s' \
+                                   % (texlabels[index].replace("$", ""), val, err, units[index])
+                    else:
+                        _anotext = r'$\mathregular{%s}$ = %4.2f %s' % (texlabels[index].replace("$", ""), val, units[index])
+
+                    anotext += [_anotext]
 
         anotext = '\n'.join(anotext)
         plot.add_anchored(
@@ -452,15 +486,16 @@ class GPMultipanelPlot(MultipanelPlot):
         subtract_orbit_model (bool, optional): if True, subtract the best-fit
             orbit model from the data and the model when plotting 
             the results. Useful for seeing the structure of correlated
-            noise in the data. Default: False. 
+            noise in the data. Default: False.
+        status (ConfigParser): (optional) result of radvel.driver.load_status on the .stat status file
 
     """
     def __init__(self, post, saveplot=None, epoch=2450000, yscale_auto=False, yscale_sigma=3.0,
                  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, text_size=9, legend_kwargs=dict(loc='best'), subtract_gp_mean_model=False,
-                 plot_likelihoods_separately=False, subtract_orbit_model=False):
+                 plot_likelihoods_separately=False, subtract_orbit_model=False, status=None):
 
         super(GPMultipanelPlot, self).__init__(
             post, saveplot=saveplot, epoch=epoch, yscale_auto=yscale_auto,
@@ -474,6 +509,8 @@ def __init__(self, post, saveplot=None, epoch=2450000, yscale_auto=False, yscale
         self.subtract_gp_mean_model = subtract_gp_mean_model
         self.plot_likelihoods_separately = plot_likelihoods_separately
         self.subtract_orbit_model = subtract_orbit_model
+        if status is not None:
+            self.status = status
 
         is_gp = False
         for like in self.like_list: