Skip to content

Commit

Permalink
Merge pull request #313 from California-Planet-Search/next-release
Browse files Browse the repository at this point in the history 
Version 1.4.0
  • Loading branch information
@bjfultn
bjfultn authored May 7, 2020
2 parents cdcaa14 + f574808 commit 75fdf4e
Show file tree
Hide file tree
Showing 19 changed files with 1,518 additions and 583 deletions.
2 changes: 2 additions & 0 deletions docs/tutorials.rst
View file
Original file line number Diff line number Diff line change
Expand Up @@ -13,3 +13,5 @@ in the `docs/tutorials` subdirectory of the radvel package.

tutorials/K2-24_Fitting+MCMC.ipynb
tutorials/GaussianProcess-tutorial.ipynb
tutorials/CustomModel-tutorial.ipynb

159 changes: 59 additions & 100 deletions docs/tutorials/164922_Fitting+MCMC.ipynb
View file

Large diffs are not rendered by default.

725 changes: 725 additions & 0 deletions docs/tutorials/CustomModel-tutorial.ipynb
View file

Large diffs are not rendered by default.

104 changes: 38 additions & 66 deletions docs/tutorials/GaussianProcess-tutorial.ipynb
View file

Large diffs are not rendered by default.

88 changes: 36 additions & 52 deletions docs/tutorials/K2-24_Fitting+MCMC.ipynb
View file

Large diffs are not rendered by default.

193 changes: 64 additions & 129 deletions docs/tutorials/SyntheticData.ipynb
View file

Large diffs are not rendered by default.

10 changes: 7 additions & 3 deletions radvel/__init__.py
View file
Original file line number Diff line number Diff line change
@@ -1,5 +1,11 @@
from __future__ import absolute_import

# turn off numpy multithreading
import os
os.environ['MKL_NUM_THREADS'] = '1'
os.environ['NUMEXPR_NUM_THREADS'] = '1'
os.environ['OMP_NUM_THREADS'] = '1'

# pre-import the big packages to avoid some warnings
import emcee # producing ABC warning
import nbsphinx # producing ABC warning
Expand All @@ -17,15 +23,13 @@
import warnings
warnings.filterwarnings("ignore")


def _custom_warningfmt(msg, *a, **b):
return "WARNING:", str(msg) + '\n'


__all__ = ['model', 'likelihood', 'posterior', 'mcmc', 'prior', 'utils',
'fitting', 'report', 'cli', 'driver', 'gp']

__version__ = '1.3.8'
__version__ = '1.4.0'
__spec__ = __name__
__package__ = __path__[0]

Expand Down
198 changes: 194 additions & 4 deletions radvel/basis.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -115,21 +115,124 @@ def to_any_basis(self, params_in, newbasis):
arbbasis_params = self.from_synth(synth_params, newbasis, keep=False)
return arbbasis_params

def v_to_any_basis(self, params_in, newbasis):
synth_vector = self.v_to_synth(params_in)
arbbasis_vector = self.v_from_synth(synth_vector, newbasis)
return arbbasis_vector

def v_to_synth(self, params_in, **kwargs):
basis_name = kwargs.setdefault('basis_name', self.name)
if isinstance(params_in, radvel.Vector):
vector = params_in.vector
else:
vector = params_in

vector2 = vector.copy()

def setvary(indices):
if kwargs.get('noVary', True):
for i in indices:
vector2[i][1] = False
vector2[i][2] = None
else:
for i in indices:
vector2[i][1] = True
vector2[i][2] = None

for num_planet in range(self.num_planets):

#per: (5*num_planet)
#tp: 1 + (5*num_planet)
#e: 2 + (5*num_planet)
#w: 3 + (5*num_planet)
#k: 4 + (5*num_planet)

if basis_name == 'per tp e w k':
return vector2

if basis_name == 'per tc e w k':
vector2[1+(5*num_planet)][0] = timetrans_to_timeperi(vector[1+(5*num_planet)][0],
vector[(5*num_planet)][0],
vector[2+(5*num_planet)][0],
vector[3+(5*num_planet)][0])
setvary([1+(5*num_planet)])

if basis_name == 'per tc se w k':
vector2[2+(5 * num_planet)][0] = vector[2+(5 * num_planet)][0] ** 2
vector2[1+(5*num_planet)][0] = timetrans_to_timeperi(vector[1+(5*num_planet)][0],
vector[(5*num_planet)][0],
vector2[2+(5*num_planet)][0],
vector[3+(5*num_planet)][0])
setvary([2+(5*num_planet), 1+(5*num_planet)])

if basis_name == 'per tc secosw sesinw logk':
vector2[4+(5*num_planet)][0] = np.exp(vector[4+(5*num_planet)][0])
vector2[2+(5*num_planet)][0] = vector[2+(5*num_planet)][0]**2 + vector[3+(5*num_planet)][0]**2
vector2[3+(5 * num_planet)][0] = np.arctan2(vector[3+(5 * num_planet)][0], vector[2+(5 * num_planet)][0])
vector2[1+(5*num_planet)][0] = timetrans_to_timeperi(vector[1+(5*num_planet)][0],
vector[(5*num_planet)][0],
vector2[2+(5*num_planet)][0],
vector2[3+(5*num_planet)][0])
setvary([4+(5*num_planet),2+(5 * num_planet),3+(5 * num_planet),1+(5 * num_planet)])

if basis_name == 'per tc secosw sesinw k':
vector2[2 + (5 * num_planet)][0] = vector[2+(5 * num_planet)][0] ** 2 + vector[3+(5 * num_planet)][0] ** 2
vector2[3 + (5 * num_planet)][0] = np.arctan2(vector[3+(5 * num_planet)][0], vector[2+(5 * num_planet)][0])
vector2[1 + (5 * num_planet)][0] = timetrans_to_timeperi(vector[1+(5 * num_planet)][0],
vector[(5 * num_planet)][0],
vector2[2+(5 * num_planet)][0],
vector2[3+(5 * num_planet)][0])
setvary([2+(5 * num_planet),3+(5 * num_planet),1+(5 * num_planet)])

if basis_name == 'logper tc secosw sesinw k':
vector2[2+(5 * num_planet)][0] = vector[2+(5 * num_planet)][0] ** 2 + vector[3+(5 * num_planet)][0] ** 2
vector2[3+(5 * num_planet)][0] = np.arctan2(vector[3+(5 * num_planet)][0], vector[2+(5 * num_planet)][0])
vector2[(5 * num_planet)][0] = np.exp(vector[(5 * num_planet)][0])
vector2[1 + (5 * num_planet)][0] = timetrans_to_timeperi(vector[1 + (5 * num_planet)][0],
vector2[(5 * num_planet)][0],
vector2[2+(5 * num_planet)][0],
vector2[3 + (5 * num_planet)][0])
setvary([2+(5 * num_planet),3+(5 * num_planet),(5 * num_planet),1 + (5 * num_planet)])

if basis_name == 'logper tc secosw sesinw logk':
vector2[2+(5 * num_planet)][0] = vector[2+(5 * num_planet)][0] ** 2 + vector[3+(5 * num_planet)][0] ** 2
vector2[3+(5 * num_planet)][0] = np.arctan2(vector[3+(5 * num_planet)][0], vector[2+(5 * num_planet)][0])
vector2[(5 * num_planet)][0] = np.exp(vector[(5 * num_planet)][0])
vector2[4 + (5 * num_planet)][0] = np.exp(vector[4 + (5 * num_planet)][0])
vector2[1 + (5 * num_planet)][0] = timetrans_to_timeperi(vector[1 + (5 * num_planet)][0],
vector2[(5 * num_planet)][0],
vector2[2+(5 * num_planet)][0],
vector2[3+(5 * num_planet)][0])
setvary([2+(5 * num_planet),3 + (5 * num_planet),(5 * num_planet),4 + (5 * num_planet),1 + (5 * num_planet)])

if basis_name == 'per tc ecosw esinw k':
vector2[2 + (5 * num_planet)][0] = np.sqrt(vector[2 + (5 * num_planet)][0] ** 2 + vector[3 + (5 * num_planet)][0] ** 2)
vector2[3+(5 * num_planet)][0] = np.arctan2(vector[3+(5 * num_planet)][0], vector[2+(5 * num_planet)][0])
vector2[1 + (5 * num_planet)][0] = timetrans_to_timeperi(vector[1 + (5 * num_planet)][0],
vector[(5 * num_planet)][0],
vector2[2+(5 * num_planet)][0],
vector2[3+(5 * num_planet)][0])
setvary([3+(5 * num_planet),2+(5 * num_planet),1 + (5 * num_planet)])

if basis_name == 'logper tp e w logk':
vector2[4 + (5 * num_planet)][0] = np.exp(vector[4 + (5 * num_planet)][0])
vector2[(5 * num_planet)][0] = np.exp(vector[(5 * num_planet)][0])
setvary([4+(5 * num_planet),(5 * num_planet)])

return vector2


def to_synth(self, params_in, **kwargs):
"""Convert to synth basis
Convert Parameters object with parameters of a given basis into the
synth basis
Args:
params_in (radvel.Parameters or pandas.DataFrame): radvel.Parameters object or pandas.Dataframe containing
orbital parameters expressed in current basis
noVary (bool [optional]): if True, set the 'vary' attribute of the returned Parameter objects
to '' (used for displaying best fit parameters)
Returns:
Parameters or DataFrame: parameters expressed in the synth basis
"""
basis_name = kwargs.setdefault('basis_name', self.name)

Expand Down Expand Up @@ -292,6 +395,93 @@ def _setpar(key, new_value):
params_out.basis = Basis('per tp e w k', self.num_planets)
return params_out

def v_from_synth(self, params_in, newbasis):
if isinstance(params_in, radvel.Vector):
vector = params_in.vector
else:
vector = params_in

vector2 = vector.copy()

def setvary(indices):
for i in indices:
vector2[i][1] = True
vector2[i][2] = None

for num_planet in range(self.num_planets):

if newbasis == 'per tc e w k':
vector2[1+(5*num_planet)][0] = timeperi_to_timetrans(vector[1+(5*num_planet)][0],
vector[(5*num_planet)][0],
vector[2+(5*num_planet)][0],
vector[3+(5*num_planet)][0])
setvary([3+(5*num_planet)])

if newbasis == 'per tc se w k':
vector2[1 + (5 * num_planet)][0] = timeperi_to_timetrans(vector[1 + (5 * num_planet)][0],
vector[(5 * num_planet)][0],
vector[2 + (5 * num_planet)][0],
vector[3 + (5 * num_planet)][0])
vector2[2+(5*num_planet)][0] = np.sqrt(vector[2+(5*num_planet)][0])
setvary([1 + (5 * num_planet),2+(5*num_planet)])

if newbasis == 'per tc secosw sesinw logk':
vector2[1 + (5 * num_planet)][0] = timeperi_to_timetrans(vector[1 + (5 * num_planet)][0],
vector[(5 * num_planet)][0],
vector[2 + (5 * num_planet)][0],
vector[3 + (5 * num_planet)][0])
vector2[2+(5*num_planet)][0] = np.sqrt(vector[2+(5*num_planet)][0])*np.cos(vector[3+(5*num_planet)][0])
vector2[3+(5 * num_planet)][0] = np.sqrt(vector[2+(5 * num_planet)][0]) * np.sin(vector[3 + (5 * num_planet)][0])
vector2[4+(5*num_planet)][0] = np.log(vector[4+(5*num_planet)][0])
setvary([1 + (5 * num_planet),2+(5*num_planet),3+(5 * num_planet),4+(5*num_planet)])

if newbasis == 'per tc secosw sesinw k':
vector2[1 + (5 * num_planet)][0] = timeperi_to_timetrans(vector[1 + (5 * num_planet)][0],
vector[(5 * num_planet)][0],
vector[2 + (5 * num_planet)][0],
vector[3 + (5 * num_planet)][0])
vector2[2 + (5 * num_planet)][0] = np.sqrt(vector[2 + (5 * num_planet)][0]) * np.cos(vector[3 + (5 * num_planet)][0])
vector2[3 + (5 * num_planet)][0] = np.sqrt(vector[2 + (5 * num_planet)][0]) * np.sin(vector[3 + (5 * num_planet)][0])
setvary([1 + (5 * num_planet), 2+(5 * num_planet), 3 + (5 * num_planet)])

if newbasis == 'logper tc secosw sesinw k':
vector2[1 + (5 * num_planet)][0] = timeperi_to_timetrans(vector[1 + (5 * num_planet)][0],
vector[(5 * num_planet)][0],
vector[2 + (5 * num_planet)][0],
vector[3 + (5 * num_planet)][0])
vector2[2 + (5 * num_planet)][0] = np.sqrt(vector[2 + (5 * num_planet)][0]) * np.cos(vector[3 + (5 * num_planet)][0])
vector2[3 + (5 * num_planet)][0] = np.sqrt(vector[2 + (5 * num_planet)][0]) * np.sin(vector[3 + (5 * num_planet)][0])
vector2[(5*num_planet)][0] = np.log(vector[(5*num_planet)][0])
setvary([1 + (5 * num_planet), 2+(5 * num_planet),3 + (5 * num_planet),(5*num_planet)])

if newbasis == 'logper tc secosw sesinw logk':
vector2[1 + (5 * num_planet)][0] = timeperi_to_timetrans(vector[1 + (5 * num_planet)][0],
vector[(5 * num_planet)][0],
vector[2 + (5 * num_planet)][0],
vector[3 + (5 * num_planet)][0])
vector2[2 + (5 * num_planet)][0] = np.sqrt(vector[2 + (5 * num_planet)][0]) * np.cos(vector[3 + (5 * num_planet)][0])
vector2[3 + (5 * num_planet)][0] = np.sqrt(vector[2 + (5 * num_planet)][0]) * np.sin(vector[3 + (5 * num_planet)][0])
vector2[(5 * num_planet)][0] = np.log(vector[(5 * num_planet)][0])
vector2[4+(5*num_planet)][0] = np.log(vector[4+(5*num_planet)][0])
setvary([1 + (5 * num_planet), 2+(5 * num_planet), 3+ (5 * num_planet), (5 * num_planet),4+(5*num_planet)])

if newbasis == 'per tc ecosw esinw k':
vector2[1 + (5 * num_planet)][0] = timeperi_to_timetrans(vector[1 + (5 * num_planet)][0],
vector[(5 * num_planet)][0],
vector[2 + (5 * num_planet)][0],
vector[3 + (5 * num_planet)][0])
vector2[2+(5*num_planet)][0] = vector[2+(5*num_planet)][0]*np.cos(vector[3+(5*num_planet)][0])
vector2[3+(5 * num_planet)][0] = vector[2+(5 * num_planet)][0] * np.sin(vector[3 + (5 * num_planet)][0])
setvary([1 + (5 * num_planet),2+(5*num_planet),3+(5 * num_planet)])

if newbasis == 'logper tp e w logk':
vector2[(5 * num_planet)][0] = np.log(vector[(5 * num_planet)][0])
vector2[4 + (5 * num_planet)][0] = np.log(vector[4 + (5 * num_planet)][0])
setvary([(5 * num_planet),4 + (5 * num_planet)])

return vector2


def from_synth(self, params_in, newbasis, **kwargs):
"""Convert from synth basis into another basis
Expand Down
16 changes: 9 additions & 7 deletions radvel/driver.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -205,8 +205,8 @@ def mcmc(args):
# Convert chains into synth basis
synthchains = chains.copy()
for par in post.params.keys():
if not post.params[par].vary:
synthchains[par] = post.params[par].value
if not post.vector.vector[post.vector.indices[par]][1]:
synthchains[par] = post.vector.vector[post.vector.indices[par]][0]

synthchains = post.params.basis.to_synth(synthchains)
synth_quantile = synthchains.quantile([0.159, 0.5, 0.841])
Expand All @@ -215,9 +215,11 @@ def mcmc(args):
# values returned by MCMC chains
post_summary = chains.quantile([0.159, 0.5, 0.841])

for k in chains.keys():
for k in chains.columns:
if k in post.params.keys():
post.params[k].value = post_summary[k][0.5]
post.vector.vector[post.vector.indices[k]][0] = post_summary[k][0.5]

post.vector.vector_to_dict()

print("Performing post-MCMC maximum likelihood fit...")
post = radvel.fitting.maxlike_fitting(post, verbose=False)
Expand All @@ -227,15 +229,15 @@ def mcmc(args):
final_chisq = np.sum(post.likelihood.residuals()**2 / (post.likelihood.errorbars()**2))
deg_of_freedom = len(post.likelihood.y) - len(post.likelihood.get_vary_params())
final_chisq_reduced = final_chisq / deg_of_freedom
post.vector.vector_to_dict()
synthparams = post.params.basis.to_synth(post.params)
post.params.update(synthparams)

print("Calculating uncertainties...")
post.uparams = {}
post.medparams = {}
post.maxparams = {}
for par in post.params.keys():
maxlike = post.params[par].value
for par in synthparams.keys():
maxlike = synthparams[par].value
med = synth_quantile[par][0.5]
high = synth_quantile[par][0.841] - med
low = med - synth_quantile[par][0.159]
Expand Down
Loading

0 comments on commit 75fdf4e

Please sign in to comment.