17447393926

Committed 03 Sep 2025 10:11PM UTC coverage: 87.165%. First build

Build # 17447393926

Build Type

Pull #395

github

Committed by

bjfultn

Commit Message

Fix Python 3.11 compatibility: replace old Python 2 raise syntax and bare except clauses

Pull Request Pull Request #395: testing new ci

Run Details

170 of 241 new or added lines in 8 files covered. (70.54%)

3735 of 4285 relevant lines covered (87.16%)

0.87 hits per line

Source File
Press 'n' to go to next uncovered line, 'b' for previous

94.41

/radvel/report.py

import numpy as np
import pandas as pd

import subprocess
import os
import tempfile
import shutil
from operator import itemgetter
from jinja2 import Environment, PackageLoader, select_autoescape

import radvel

env = Environment(loader=PackageLoader('radvel', 'templates'))
print_basis = 'per tc tp e w k'
units = {
    'per': 'days',
    'tp': 'JD',
    'tc': 'JD',
    'e': '',
    'w': 'radians',
    'k': 'm s$^{-1}$',
    'logk': r'$\ln{(\rm m\ s^{-1})}$',
    'secosw': '',
    'sesinw': '',
    'gamma': 'm s$-1$',
    'jitter': 'm s$-1$',
    'logjit': r'$\ln{(\rm m\ s^{-1})}$',
    'jit': r'$\rm m\ s^{-1}$',
    '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': '',
    # 'mpsini': '$M_\earth$',
    # 'rp': '$R_\earth$',
    # 'rhop': 'g cm$^{-3}$',
}


class RadvelReport(object):
    """Radvel report

    Class to handle the creation of the radvel summary PDF

    Args:
        planet (planet object): planet configuration object loaded in
            `kepfit.py` using `imp.load_source`
        post (radvel.posterior):
            radvel.posterior object containing the best-fit parameters in
                post.params
        compstats (dict): dictionary of model comparison results from `radvel ic`
        derived (bool): included table of derived parameters
        chains (DataFrame): output DataFrame from a `radvel.mcmc` run
        criterion (DataFrame): output DataFrame from a 'radvel.mcmc' run
    """

    def __init__(self, planet, post, chains, minafactor, maxarchange, maxgr, mintz, compstats=None,
                 derived=False):
        self.planet = planet
        self.post = post
        self.starname = planet.starname
        self.starname_tex = planet.starname.replace('_', '\\_')
        self.runname = self.starname_tex
        self.derived = derived

        post.params = post.params.basis.to_synth(post.params)
        post.params = post.params.basis.from_synth(
            post.params, print_basis
        )
        self.latex_dict = post.params.tex_labels()

        for p in post.params.keys():
            if p not in chains.columns:
                chains[p] = post.params[p].value

        self.chains = post.params.basis.to_synth(
            chains, basis_name=planet.fitting_basis
        )
        self.chains = post.params.basis.from_synth(
            self.chains, print_basis
        )
        self.minafactor = minafactor
        self.maxarchange = maxarchange
        self.maxgr = maxgr
        self.mintz = mintz
        self.quantiles = self.chains.quantile([0.159, 0.5, 0.841])
        self.compstats = compstats
        self.num_planets = self.post.params.num_planets

    def texdoc(self):
        """TeX for entire document

        Returns:
            str: TeX code for report
        """
        reportkw = {}
        reportkw['version'] = radvel.__version__

        # Render TeX for figures
        figtypes = ['rv_multipanel', 'corner', 'corner_derived_pars']
        for figtype in figtypes:
            infile = "{}_{}.pdf".format(self.runname, figtype)
            tmpfile = 'fig_{}.tex'.format(figtype)
            key = 'fig_{}'.format(figtype)
            if os.path.exists(infile):
                t = env.get_template(tmpfile)
                reportkw[key] = t.render(report=self, infile=infile)

        # Render TeX for tables
        textable = TexTable(self)
        reportkw['tab_crit'] = textable.tab_crit()
        reportkw['tab_rv'] = textable.tab_rv()
        reportkw['tab_params'] = textable.tab_params()
        reportkw['tab_derived'] = textable.tab_derived()
        reportkw['tab_prior_summary'] = textable.tab_prior_summary()

        if self.compstats is not None:
            reportkw['tab_comparison'] = textable.tab_comparison()

        t = env.get_template('report.tex')
        out = t.render(report=self, **reportkw)
        return out

    def compile(self, pdfname, latex_compiler='pdflatex', depfiles=[]):
        """Compile radvel report

        Compile the radvel report from a string containing TeX code
        and save the resulting PDF to a file.

        Args:
            pdfname (str): name of the output PDF file
            latex_compiler (str): path to latex compiler
            depfiles (list): list of file names of dependencies needed for
                LaTex compilation (e.g. figure files)

        """
        texname = os.path.basename(pdfname).split('.')[0] + '.tex'
        current = os.getcwd()
        temp = tempfile.mkdtemp()
        for fname in depfiles:
            shutil.copy2(
                os.path.join(current, fname), os.path.join(temp, fname)
            )

        os.chdir(temp)

        f = open(texname, 'w')
        f.write(self.texdoc())
        f.close()

        shutil.copy(texname, current)

        try:
            for i in range(3):
                # LaTex likes to be compiled a few times
                # to get the table widths correct
                proc = subprocess.Popen(
                    [latex_compiler, texname], stdout=subprocess.PIPE,
                )
                proc.communicate()  # Let the subprocess complete
        except (OSError):
            msg = """
WARNING: REPORT: could not run %s. Ensure that %s is in your PATH
or pass in the path as an argument
""" % (latex_compiler, latex_compiler)
            print(msg)
            return

        shutil.copy(pdfname, current)

        shutil.rmtree(temp)
        os.chdir(current)


class TexTable(RadvelReport):
    """LaTeX table

    Class to handle generation of the LaTeX tables within the summary PDF.

    Args:
        report (radvel.report.RadvelReport): radvel report object
        full (bool): get full-length RV table [default: True]
    """

    def __init__(self, report):
        self.report = report
        self.post = report.post
        self.quantiles = report.quantiles
        self.fitting_basis = report.post.params.basis.name
        self.minafactor = report.minafactor
        self.maxarchange = report.maxarchange
        self.maxgr = report.maxgr
        self.mintz = report.mintz

    def _row(self, param, unit):
        """
        Helper function to output the rows in the parameter table
        """
        if param not in self.quantiles.keys():
            param = param[:-1]

        if unit == 'radians':
            par = radvel.utils.geterr(self.report.chains[param], angular=True)
            med, low, high = par
        else:
            med = self.quantiles[param][0.5]
            low = self.quantiles[param][0.5] - self.quantiles[param][0.159]
            high = self.quantiles[param][0.841] - self.quantiles[param][0.5]

        maxlike = self.post.maxparams[param]
        tex = self.report.latex_dict[param]

        low = radvel.utils.round_sig(low)
        high = radvel.utils.round_sig(high)
        maxlike, errlow, errhigh = radvel.utils.sigfig(maxlike, low, high)
        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)
            errfmt = ''
        else:
            if errhigh == errlow:
                errfmt = r'\pm %s' % errhigh
            else:
                errfmt = '^{+%s}_{-%s}' % (errhigh, errlow)

        row = "%s & $%s%s$ & $%s$ & %s" % (tex, med, errfmt, maxlike, unit)
        return row

    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 = []

        nloop = self.report.planet.nplanets+1
        if dontloop:
            nloop=2

        for n in range(1, nloop):
            for p in basis.split():  # loop over variables
                par = p+str(n)

                # get unit for parameter
                unit = units.get(p, '')

                # try to remove suffix
                if unit == '' and par not in units.keys():
                    for s in suffixes:
                        if s in p:
                            unit = units.get(p.replace(s, ''), '')
                            break

                # if still can't find units get it by the parameter itself
                # (derived parameters)
                if unit == '':
                    unit = units.get(par, '')

                try:
                    row = self._row(par, unit)
                except KeyError:
                    continue

                rows.append(row)

        return rows

    def tab_prior_summary(self, name_in_title=False):
        """Summary of priors

        Args:
            name_in_title (Bool [optional]): if True, include
                the name of the star in the table title
        """
        texdict = self.post.likelihood.params.tex_labels()
        prior_list = self.post.priors
        rows = []
        for prior in prior_list:
            row = prior.__str__()
            if not row.endswith("\\\\"):
                row = row + "\\\\"
            rows.append(row)

        kw = {}
        if name_in_title:
            kw['title'] = "{} Summary of Priors".format(self.report.starname)
        else:
            kw['title'] = "Summary of Priors"
        tmpfile = 'tab_prior_summary.tex'
        t = env.get_template(tmpfile)
        out = t.render(rows=rows, **kw)
        return out

    def tab_rv(self, name_in_title=False, max_lines=50):
        """Table of input velocities

        Args:
            name_in_title (Bool [optional]): if True, include
                the name of the star in the table title
        """

        kw = {}
        nvels = len(self.post.likelihood.x)

        if max_lines is None:
            iters = range(nvels)
            kw['notes'] = ''
        else:
            max_lines = int(np.round(max_lines))
            iters = range(nvels)[:max_lines]
            kw['notes'] = r"""Only the first %d of %d RVs are displayed in this table. \
Use \texttt{radvel table -t rv} to save the full \LaTeX\ table as a separate file.""" % (max_lines, nvels)

        rows = []
        for i in iters:
            t = self.post.likelihood.x[i]
            v = self.post.likelihood.y[i]
            e = self.post.likelihood.yerr[i]
            inst = self.post.likelihood.telvec[i]
            if '_' in inst:
                inst = inst.replace('_', r'$_{\rm ')
                inst += '}$'
            row = "{:.5f} & {:.2f} & {:.2f} & {:s}".format(t, v, e, inst)
            rows.append(row)

        if name_in_title:
            kw['title'] = "{} Radial Velocities".format(self.report.starname)
        else:
            kw['title'] = "Radial Velocities"
        tmpfile = 'tab_rv.tex'
        t = env.get_template(tmpfile)
        out = t.render(rows=rows, **kw)
        return out

    def tab_params(self, name_in_title=False):
        """ Table of final parameter values
        Args:
            name_in_title (Bool [optional]): if True, include
                the name of the star in the table title
        """
        # Sort extra params
        ep = []
        order = ['gamma', 'dvdt', 'curv', 'jit']
        for o in order:
            op = []
            for p in self.post.likelihood.extra_params:
                if o in p:
                    op.append(p)
            if len(op)==0:
                op = [o]
            [ep.append(i) for i in sorted(op)[::-1]]

        # Add GP parameters
        for par in self.report.post.likelihood.vector.names:
            if par.startswith('gp_'):
                ep.append(par)

        ep = ' '.join(ep)

        kw = {}
        kw['fitting_basis_rows'] = self._data(self.fitting_basis)
        kw['print_basis_rows'] = self._data(print_basis)
        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:
            kw['title'] = "{} Posteriors".format(self.report.starname)
        else:
            kw['title'] = "Posteriors"
        tmpfile = 'tab_params.tex'
        t = env.get_template(tmpfile)
        out = t.render(**kw)
        return out

    def tab_crit(self, name_in_title=False):
        """Table of final convergence criterion values
        Args:
            name_in_title (Bool [optional]): if True, include
                the name of the star in the table title
        """

        names = ['minAfactor', 'maxArchange', 'maxGR', 'minTz']
        values = [self.minafactor, self.maxarchange, self.maxgr, self.mintz]
        rows = []
        for i in range(0,len(names)):
            val = values[i]
            val_str = '${:7.3f}$'.format(float(val)) if val is not None else "None"
            rows.append(r"%s & %s" % (names[i], val_str))

        kw = dict()
        kw['rows'] = rows
        if name_in_title:
            kw['title'] = "{} Final Convergence Criterion".format(self.report.starname)
        else:
            kw['title'] = "Final Convergence Criterion"

        tmpfile = 'tab_crit.tex'
        t = env.get_template(tmpfile)
        out = t.render(**kw)

        return out

    def tab_derived(self, name_in_title=False):
        """ Table of derived parameter values
        Args:
            name_in_title (Bool [optional]): if True, include
                the name of the star in the table title
        """

        if not self.report.derived:
            return ""

        dpl = radvel.plot.mcmc_plots.DerivedPlot(self.report.chains, self.report.planet)
        derived_params = dpl.labels
        derived_basis = ' '.join(set([s[:-1] for s in derived_params]))
        derived_tex = dpl.texlabels
        derived_units = dpl.units

        self.report.latex_dict.update(dict(zip(derived_params, derived_tex)))
        units.update(dict(zip(derived_params, derived_units)))

        self.quantiles = dpl.chains.quantile([0.159, 0.5, 0.841])

        for par in derived_params:
            # self.report.post.maxparams[par] = self.report.chains[par].iloc[
            #     self.report.chains['lnprobability'].argmax]
            self.post.maxparams[par] = dpl.chains.loc[dpl.chains['lnprobability'].idxmax(), par]

        kw = dict()
        kw['derived_rows'] = self._data(derived_basis)
        if name_in_title:
            kw['title'] = "{} Derived Posteriors".format(self.report.starname)
        else:
            kw['title'] = "Derived Posteriors"
        tmpfile = 'tab_derived.tex'
        t = env.get_template(tmpfile)
        out = t.render(**kw)

        return out

    def tab_comparison(self):
        """Model comparisons
        """
        statsdict = self.report.compstats

        if statsdict is None or len(statsdict) < 1:
            return ""

        statsdict_sorted = sorted(statsdict, key=itemgetter('AICc'), reverse=False)

        n_test = len(statsdict_sorted)
        if n_test > 50:
            print("Warning, the number of model comparisons is very"
                  + " large. Printing 50 best models.\nConsider using"
                  + " the --unmixed flag when performing ic comparisons")
            n_test = 50
            # statsdict_sorted = statsdict_sorted[:50]

        statskeys = statsdict_sorted[0].keys()
        coldefs = r"\begin{deluxetable*}{%s}" % ('l'+'l'+'r'*(len(statskeys)-1) + 'r')
        head = r"\tablehead{"
        head += r"\colhead{AICc Qualitative Comparison}"
        head += r" & \colhead{Free Parameters}"
        for s in statskeys:
            if s == 'Free Params':
                pass
            else:
                head += r" & \colhead{%s}" % s
        head += r" & \colhead{$\Delta$AICc}"
        head += r"}"

        minAIC = statsdict_sorted[0]['AICc'][0]
        # See Burnham + Anderson 2004
        deltaAIClevels = [0., 2., 4., 10.]
        deltaAICmessages = ["Nearly Indistinguishable", "Somewhat Disfavored", \
            "Strongly Disfavored", "Ruled Out"]
        deltaAICtrigger = 0
        maxtrigger = len(deltaAIClevels)

        rows = []
        for i in range(n_test):
            row = ""
            for s in statsdict_sorted[i].keys():
                val = statsdict_sorted[i][s][0]
                if type(val) is int:
                    row += " & %s" % str(val)
                elif type(val) is float:
                    row += " & %.2f" % val
                elif type(val) is str:
                    row += " & %s" % val
                elif type(val) is list:
                    row += " &"
                    for item in val:
                        row += " %s," %item
                    row += r" {$\gamma$}"
                    #row = row[:-1]
                else:
                    raise ValueError("Failed to format values for LaTeX: {}  {}".format(s, val))
            row += " & %.2f" % (statsdict_sorted[i]['AICc'][0] - minAIC)
            # row = row[3:]
            if i == 0:
                # row = "{\\bf" + row + "}"
                row = "AICc Favored Model"+ row
            appendhline = False
            if (deltaAICtrigger < maxtrigger) and ((statsdict_sorted[i]['AICc'][0] - minAIC) > deltaAIClevels[deltaAICtrigger]):
                deltaAICtrigger += 1
                while (deltaAICtrigger < maxtrigger) and ((statsdict_sorted[i]['AICc'][0] - minAIC) > deltaAIClevels[deltaAICtrigger]):
                    deltaAICtrigger += 1
                row = deltaAICmessages[deltaAICtrigger-1] + row
                appendhline = True
            if appendhline or (i == 1):
                rows.append(r"\hline")
            rows.append(row)

        t = env.get_template('tab_comparison.tex')
        out = t.render(coldefs=coldefs, head=head, rows=rows)
        return out

1	import numpy as np	1✔
2	import pandas as pd	1✔
3
4	import subprocess	1✔
5	import os	1✔
6	import tempfile	1✔
7	import shutil	1✔
8	from operator import itemgetter	1✔
9	from jinja2 import Environment, PackageLoader, select_autoescape	1✔
10
11	import radvel	1✔
12
13	env = Environment(loader=PackageLoader('radvel', 'templates'))	1✔
14	print_basis = 'per tc tp e w k'	1✔
15	units = {	1✔
16	'per': 'days',
17	'tp': 'JD',
18	'tc': 'JD',
19	'e': '',
20	'w': 'radians',
21	'k': 'm s$^{-1}$',
22	'logk': r'$\ln{(\rm m\ s^{-1})}$',
23	'secosw': '',
24	'sesinw': '',
25	'gamma': 'm s$-1$',
26	'jitter': 'm s$-1$',
27	'logjit': r'$\ln{(\rm m\ s^{-1})}$',
28	'jit': r'$\rm m\ s^{-1}$',
29	'dvdt': 'm s$^{-1}$ d$^{-1}$',
30	'curv': 'm s$^{-1}$ d$^{-2}$',
31	'gp_amp': 'm s$-1$',
32	'gp_explength': 'days',
33	'gp_per': 'days',
34	'gp_perlength': '',
35	# 'mpsini': '$M_\earth$',
36	# 'rp': '$R_\earth$',
37	# 'rhop': 'g cm$^{-3}$',
38	}
39
40
41	class RadvelReport(object):	1✔
42	"""Radvel report
43
44	Class to handle the creation of the radvel summary PDF
45
46	Args:
47	planet (planet object): planet configuration object loaded in
48	`kepfit.py` using `imp.load_source`
49	post (radvel.posterior):
50	radvel.posterior object containing the best-fit parameters in
51	post.params
52	compstats (dict): dictionary of model comparison results from `radvel ic`
53	derived (bool): included table of derived parameters
54	chains (DataFrame): output DataFrame from a `radvel.mcmc` run
55	criterion (DataFrame): output DataFrame from a 'radvel.mcmc' run
56	"""
57
58	def __init__(self, planet, post, chains, minafactor, maxarchange, maxgr, mintz, compstats=None,	1✔
59	derived=False):
60	self.planet = planet	1✔
61	self.post = post	1✔
62	self.starname = planet.starname	1✔
63	self.starname_tex = planet.starname.replace('_', '\\_')	1✔
64	self.runname = self.starname_tex	1✔
65	self.derived = derived	1✔
66
67	post.params = post.params.basis.to_synth(post.params)	1✔
68	post.params = post.params.basis.from_synth(	1✔
69	post.params, print_basis
70	)
71	self.latex_dict = post.params.tex_labels()	1✔
72
73	for p in post.params.keys():	1✔
74	if p not in chains.columns:	1✔
75	chains[p] = post.params[p].value	1✔
76
77	self.chains = post.params.basis.to_synth(	1✔
78	chains, basis_name=planet.fitting_basis
79	)
80	self.chains = post.params.basis.from_synth(	1✔
81	self.chains, print_basis
82	)
83	self.minafactor = minafactor	1✔
84	self.maxarchange = maxarchange	1✔
85	self.maxgr = maxgr	1✔
86	self.mintz = mintz	1✔
87	self.quantiles = self.chains.quantile([0.159, 0.5, 0.841])	1✔
88	self.compstats = compstats	1✔
89	self.num_planets = self.post.params.num_planets	1✔
90
91	def texdoc(self):	1✔
92	"""TeX for entire document
93
94	Returns:
95	str: TeX code for report
96	"""
97	reportkw = {}	1✔
98	reportkw['version'] = radvel.__version__	1✔
99
100	# Render TeX for figures
101	figtypes = ['rv_multipanel', 'corner', 'corner_derived_pars']	1✔
102	for figtype in figtypes:	1✔
103	infile = "{}_{}.pdf".format(self.runname, figtype)	1✔
104	tmpfile = 'fig_{}.tex'.format(figtype)	1✔
105	key = 'fig_{}'.format(figtype)	1✔
106	if os.path.exists(infile):	1✔
107	t = env.get_template(tmpfile)	1✔
108	reportkw[key] = t.render(report=self, infile=infile)	1✔
109
110	# Render TeX for tables
111	textable = TexTable(self)	1✔
112	reportkw['tab_crit'] = textable.tab_crit()	1✔
113	reportkw['tab_rv'] = textable.tab_rv()	1✔
114	reportkw['tab_params'] = textable.tab_params()	1✔
115	reportkw['tab_derived'] = textable.tab_derived()	1✔
116	reportkw['tab_prior_summary'] = textable.tab_prior_summary()	1✔
117
118	if self.compstats is not None:	1✔
119	reportkw['tab_comparison'] = textable.tab_comparison()	1✔
120
121	t = env.get_template('report.tex')	1✔
122	out = t.render(report=self, **reportkw)	1✔
123	return out	1✔
124
125	def compile(self, pdfname, latex_compiler='pdflatex', depfiles=[]):	1✔
126	"""Compile radvel report
127
128	Compile the radvel report from a string containing TeX code
129	and save the resulting PDF to a file.
130
131	Args:
132	pdfname (str): name of the output PDF file
133	latex_compiler (str): path to latex compiler
134	depfiles (list): list of file names of dependencies needed for
135	LaTex compilation (e.g. figure files)
136
137	"""
138	texname = os.path.basename(pdfname).split('.')[0] + '.tex'	1✔
139	current = os.getcwd()	1✔
140	temp = tempfile.mkdtemp()	1✔
141	for fname in depfiles:	1✔
142	shutil.copy2(	1✔
143	os.path.join(current, fname), os.path.join(temp, fname)
144	)
145
146	os.chdir(temp)	1✔
147
148	f = open(texname, 'w')	1✔
149	f.write(self.texdoc())	1✔
150	f.close()	1✔
151
152	shutil.copy(texname, current)	1✔
153
154	try:	1✔
155	for i in range(3):	1✔
156	# LaTex likes to be compiled a few times
157	# to get the table widths correct
158	proc = subprocess.Popen(	1✔
159	[latex_compiler, texname], stdout=subprocess.PIPE,
160	)
161	proc.communicate() # Let the subprocess complete	×
162	except (OSError):	1✔
163	msg = """	1✔
164	WARNING: REPORT: could not run %s. Ensure that %s is in your PATH
165	or pass in the path as an argument
166	""" % (latex_compiler, latex_compiler)
167	print(msg)	1✔
168	return	1✔
169
170	shutil.copy(pdfname, current)	×
171
172	shutil.rmtree(temp)	×
173	os.chdir(current)	×
174
175
176	class TexTable(RadvelReport):	1✔
177	"""LaTeX table
178
179	Class to handle generation of the LaTeX tables within the summary PDF.
180
181	Args:
182	report (radvel.report.RadvelReport): radvel report object
183	full (bool): get full-length RV table [default: True]
184	"""
185
186	def __init__(self, report):	1✔
187	self.report = report	1✔
188	self.post = report.post	1✔
189	self.quantiles = report.quantiles	1✔
190	self.fitting_basis = report.post.params.basis.name	1✔
191	self.minafactor = report.minafactor	1✔
192	self.maxarchange = report.maxarchange	1✔
193	self.maxgr = report.maxgr	1✔
194	self.mintz = report.mintz	1✔
195
196	def _row(self, param, unit):	1✔
197	"""
198	Helper function to output the rows in the parameter table
199	"""
200	if param not in self.quantiles.keys():	1✔
201	param = param[:-1]	1✔
202
203	if unit == 'radians':	1✔
204	par = radvel.utils.geterr(self.report.chains[param], angular=True)	1✔
205	med, low, high = par	1✔
206	else:
207	med = self.quantiles[param][0.5]	1✔
208	low = self.quantiles[param][0.5] - self.quantiles[param][0.159]	1✔
209	high = self.quantiles[param][0.841] - self.quantiles[param][0.5]	1✔
210
211	maxlike = self.post.maxparams[param]	1✔
212	tex = self.report.latex_dict[param]	1✔
213
214	low = radvel.utils.round_sig(low)	1✔
215	high = radvel.utils.round_sig(high)	1✔
216	maxlike, errlow, errhigh = radvel.utils.sigfig(maxlike, low, high)	1✔
217	med, errlow, errhigh = radvel.utils.sigfig(med, low, high)	1✔
218
219	if min(errlow,errhigh) <= 1e-12:	1✔
220	med = maxlike = r"\equiv%s" % round(self.quantiles[param][0.5], 4)	1✔
221	errfmt = ''	1✔
222	else:
223	if errhigh == errlow:	1✔
224	errfmt = r'\pm %s' % errhigh	1✔
225	else:
226	errfmt = '^{+%s}_{-%s}' % (errhigh, errlow)	1✔
227
228	row = "%s & $%s%s$ & $%s$ & %s" % (tex, med, errfmt, maxlike, unit)	1✔
229	return row	1✔
230
231	def _data(self, basis, dontloop=False):	1✔
232	"""
233	Helper function to output the rows in the parameter table
234
235	Args:
236	basis (str): name of Basis object (see basis.py) to be printed
237	dontloop (Bool): if True, don't loop over number of planets (useful for
238	printing out gamma, dvdt, jitter, curv)
239	"""
240	suffixes = ['_'+j for j in self.report.post.likelihood.suffixes]	1✔
241	rows = []	1✔
242
243	nloop = self.report.planet.nplanets+1	1✔
244	if dontloop:	1✔
245	nloop=2	1✔
246
247	for n in range(1, nloop):	1✔
248	for p in basis.split(): # loop over variables	1✔
249	par = p+str(n)	1✔
250
251	# get unit for parameter
252	unit = units.get(p, '')	1✔
253
254	# try to remove suffix
255	if unit == '' and par not in units.keys():	1✔
256	for s in suffixes:	1✔
257	if s in p:	1✔
258	unit = units.get(p.replace(s, ''), '')	1✔
259	break	1✔
260
261	# if still can't find units get it by the parameter itself
262	# (derived parameters)
263	if unit == '':	1✔
264	unit = units.get(par, '')	1✔
265
266	try:	1✔
267	row = self._row(par, unit)	1✔
268	except KeyError:	×
269	continue	×
270
271	rows.append(row)	1✔
272
273	return rows	1✔
274
275	def tab_prior_summary(self, name_in_title=False):	1✔
276	"""Summary of priors
277
278	Args:
279	name_in_title (Bool [optional]): if True, include
280	the name of the star in the table title
281	"""
282	texdict = self.post.likelihood.params.tex_labels()	1✔
283	prior_list = self.post.priors	1✔
284	rows = []	1✔
285	for prior in prior_list:	1✔
286	row = prior.__str__()	1✔
287	if not row.endswith("\\\\"):	1✔
288	row = row + "\\\\"	1✔
289	rows.append(row)	1✔
290
291	kw = {}	1✔
292	if name_in_title:	1✔
293	kw['title'] = "{} Summary of Priors".format(self.report.starname)	×
294	else:
295	kw['title'] = "Summary of Priors"	1✔
296	tmpfile = 'tab_prior_summary.tex'	1✔
297	t = env.get_template(tmpfile)	1✔
298	out = t.render(rows=rows, **kw)	1✔
299	return out	1✔
300
301	def tab_rv(self, name_in_title=False, max_lines=50):	1✔
302	"""Table of input velocities
303
304	Args:
305	name_in_title (Bool [optional]): if True, include
306	the name of the star in the table title
307	"""
308
309	kw = {}	1✔
310	nvels = len(self.post.likelihood.x)	1✔
311
312	if max_lines is None:	1✔
313	iters = range(nvels)	1✔
314	kw['notes'] = ''	1✔
315	else:
316	max_lines = int(np.round(max_lines))	1✔
317	iters = range(nvels)[:max_lines]	1✔
318	kw['notes'] = r"""Only the first %d of %d RVs are displayed in this table. \	1✔
319	Use \texttt{radvel table -t rv} to save the full \LaTeX\ table as a separate file.""" % (max_lines, nvels)
320
321	rows = []	1✔
322	for i in iters:	1✔
323	t = self.post.likelihood.x[i]	1✔
324	v = self.post.likelihood.y[i]	1✔
325	e = self.post.likelihood.yerr[i]	1✔
326	inst = self.post.likelihood.telvec[i]	1✔
327	if '_' in inst:	1✔
328	inst = inst.replace('_', r'$_{\rm ')	×
329	inst += '}$'	×
330	row = "{:.5f} & {:.2f} & {:.2f} & {:s}".format(t, v, e, inst)	1✔
331	rows.append(row)	1✔
332
333	if name_in_title:	1✔
334	kw['title'] = "{} Radial Velocities".format(self.report.starname)	×
335	else:
336	kw['title'] = "Radial Velocities"	1✔
337	tmpfile = 'tab_rv.tex'	1✔
338	t = env.get_template(tmpfile)	1✔
339	out = t.render(rows=rows, **kw)	1✔
340	return out	1✔
341
342	def tab_params(self, name_in_title=False):	1✔
343	""" Table of final parameter values
344	Args:
345	name_in_title (Bool [optional]): if True, include
346	the name of the star in the table title
347	"""
348	# Sort extra params
349	ep = []	1✔
350	order = ['gamma', 'dvdt', 'curv', 'jit']	1✔
351	for o in order:	1✔
352	op = []	1✔
353	for p in self.post.likelihood.extra_params:	1✔
354	if o in p:	1✔
355	op.append(p)	1✔
356	if len(op)==0:	1✔
357	op = [o]	1✔
358	[ep.append(i) for i in sorted(op)[::-1]]	1✔
359
360	# Add GP parameters
361	for par in self.report.post.likelihood.vector.names:	1✔
362	if par.startswith('gp_'):	1✔
363	ep.append(par)	×
364
365	ep = ' '.join(ep)	1✔
366
367	kw = {}	1✔
368	kw['fitting_basis_rows'] = self._data(self.fitting_basis)	1✔
369	kw['print_basis_rows'] = self._data(print_basis)	1✔
370	kw['ep_rows'] = self._data(ep, dontloop=True)	1✔
371	kw['nlinks'] = len(self.report.chains)	1✔
372	kw['time_base'] = self.report.post.likelihood.model.time_base	1✔
373	if name_in_title:	1✔
374	kw['title'] = "{} Posteriors".format(self.report.starname)	×
375	else:
376	kw['title'] = "Posteriors"	1✔
377	tmpfile = 'tab_params.tex'	1✔
378	t = env.get_template(tmpfile)	1✔
379	out = t.render(**kw)	1✔
380	return out	1✔
381
382	def tab_crit(self, name_in_title=False):	1✔
383	"""Table of final convergence criterion values
384	Args:
385	name_in_title (Bool [optional]): if True, include
386	the name of the star in the table title
387	"""
388
389	names = ['minAfactor', 'maxArchange', 'maxGR', 'minTz']	1✔
390	values = [self.minafactor, self.maxarchange, self.maxgr, self.mintz]	1✔
391	rows = []	1✔
392	for i in range(0,len(names)):	1✔
393	val = values[i]	1✔
394	val_str = '${:7.3f}$'.format(float(val)) if val is not None else "None"	1✔
395	rows.append(r"%s & %s" % (names[i], val_str))	1✔
396
397	kw = dict()	1✔
398	kw['rows'] = rows	1✔
399	if name_in_title:	1✔
400	kw['title'] = "{} Final Convergence Criterion".format(self.report.starname)	×
401	else:
402	kw['title'] = "Final Convergence Criterion"	1✔
403
404	tmpfile = 'tab_crit.tex'	1✔
405	t = env.get_template(tmpfile)	1✔
406	out = t.render(**kw)	1✔
407
408	return out	1✔
409
410	def tab_derived(self, name_in_title=False):	1✔
411	""" Table of derived parameter values
412	Args:
413	name_in_title (Bool [optional]): if True, include
414	the name of the star in the table title
415	"""
416
417	if not self.report.derived:	1✔
418	return ""	1✔
419
420	dpl = radvel.plot.mcmc_plots.DerivedPlot(self.report.chains, self.report.planet)	1✔
421	derived_params = dpl.labels	1✔
422	derived_basis = ' '.join(set([s[:-1] for s in derived_params]))	1✔
423	derived_tex = dpl.texlabels	1✔
424	derived_units = dpl.units	1✔
425
426	self.report.latex_dict.update(dict(zip(derived_params, derived_tex)))	1✔
427	units.update(dict(zip(derived_params, derived_units)))	1✔
428
429	self.quantiles = dpl.chains.quantile([0.159, 0.5, 0.841])	1✔
430
431	for par in derived_params:	1✔
432	# self.report.post.maxparams[par] = self.report.chains[par].iloc[
433	# self.report.chains['lnprobability'].argmax]
434	self.post.maxparams[par] = dpl.chains.loc[dpl.chains['lnprobability'].idxmax(), par]	1✔
435
436	kw = dict()	1✔
437	kw['derived_rows'] = self._data(derived_basis)	1✔
438	if name_in_title:	1✔
439	kw['title'] = "{} Derived Posteriors".format(self.report.starname)	×
440	else:
441	kw['title'] = "Derived Posteriors"	1✔
442	tmpfile = 'tab_derived.tex'	1✔
443	t = env.get_template(tmpfile)	1✔
444	out = t.render(**kw)	1✔
445
446	return out	1✔
447
448	def tab_comparison(self):	1✔
449	"""Model comparisons
450	"""
451	statsdict = self.report.compstats	1✔
452
453	if statsdict is None or len(statsdict) < 1:	1✔
454	return ""	×
455
456	statsdict_sorted = sorted(statsdict, key=itemgetter('AICc'), reverse=False)	1✔
457
458	n_test = len(statsdict_sorted)	1✔
459	if n_test > 50:	1✔
460	print("Warning, the number of model comparisons is very"	1✔
461	+ " large. Printing 50 best models.\nConsider using"
462	+ " the --unmixed flag when performing ic comparisons")
463	n_test = 50	1✔
464	# statsdict_sorted = statsdict_sorted[:50]
465
466	statskeys = statsdict_sorted[0].keys()	1✔
467	coldefs = r"\begin{deluxetable}{%s}" % ('l'+'l'+'r'(len(statskeys)-1) + 'r')	1✔
468	head = r"\tablehead{"	1✔
469	head += r"\colhead{AICc Qualitative Comparison}"	1✔
470	head += r" & \colhead{Free Parameters}"	1✔
471	for s in statskeys:	1✔
472	if s == 'Free Params':	1✔
473	pass	1✔
474	else:
475	head += r" & \colhead{%s}" % s	1✔
476	head += r" & \colhead{$\Delta$AICc}"	1✔
477	head += r"}"	1✔
478
479	minAIC = statsdict_sorted[0]['AICc'][0]	1✔
480	# See Burnham + Anderson 2004
481	deltaAIClevels = [0., 2., 4., 10.]	1✔
482	deltaAICmessages = ["Nearly Indistinguishable", "Somewhat Disfavored", \	1✔
483	"Strongly Disfavored", "Ruled Out"]
484	deltaAICtrigger = 0	1✔
485	maxtrigger = len(deltaAIClevels)	1✔
486
487	rows = []	1✔
488	for i in range(n_test):	1✔
489	row = ""	1✔
490	for s in statsdict_sorted[i].keys():	1✔
491	val = statsdict_sorted[i][s][0]	1✔
492	if type(val) is int:	1✔
493	row += " & %s" % str(val)	1✔
494	elif type(val) is float:	1✔
495	row += " & %.2f" % val	1✔
496	elif type(val) is str:	1✔
497	row += " & %s" % val	×
498	elif type(val) is list:	1✔
499	row += " &"	1✔
500	for item in val:	1✔
501	row += " %s," %item	1✔
502	row += r" {$\gamma$}"	1✔
503	#row = row[:-1]
504	else:
NEW 505	raise ValueError("Failed to format values for LaTeX: {} {}".format(s, val))	×
506	row += " & %.2f" % (statsdict_sorted[i]['AICc'][0] - minAIC)	1✔
507	# row = row[3:]
508	if i == 0:	1✔
509	# row = "{\\bf" + row + "}"
510	row = "AICc Favored Model"+ row	1✔
511	appendhline = False	1✔
512	if (deltaAICtrigger < maxtrigger) and ((statsdict_sorted[i]['AICc'][0] - minAIC) > deltaAIClevels[deltaAICtrigger]):	1✔
513	deltaAICtrigger += 1	1✔
514	while (deltaAICtrigger < maxtrigger) and ((statsdict_sorted[i]['AICc'][0] - minAIC) > deltaAIClevels[deltaAICtrigger]):	1✔
515	deltaAICtrigger += 1	1✔
516	row = deltaAICmessages[deltaAICtrigger-1] + row	1✔
517	appendhline = True	1✔
518	if appendhline or (i == 1):	1✔
519	rows.append(r"\hline")	1✔
520	rows.append(row)	1✔
521
522	t = env.get_template('tab_comparison.tex')	1✔
523	out = t.render(coldefs=coldefs, head=head, rows=rows)	1✔
524	return out	1✔

California-Planet-Search / radvel / 17447393926

Source File Press 'n' to go to next uncovered line, 'b' for previous

Source File
Press 'n' to go to next uncovered line, 'b' for previous