Secure your code as it's written. Use Snyk Code to scan source code in minutes - no build needed - and fix issues immediately.
with model:
if 'K' not in pars:
if sigma_K0 is None or P0 is None:
raise ValueError("If using the default prior form on K, you "
"must pass in a variance scale (sigma_K0) "
"and a reference period (P0)")
# Default prior on semi-amplitude: scales with period and
# eccentricity such that it is flat with companion mass
v_unit = sigma_K0.unit
out_pars['K'] = xu.with_unit(FixedCompanionMass('K', P=P, e=e,
sigma_K0=sigma_K0,
P0=P0),
v_unit)
else:
v_unit = getattr(pars['K'], xu.UNIT_ATTR_NAME, u.one)
for i, name in enumerate(v_names):
if name not in pars:
# Default priors are independent gaussians
# FIXME: make mean, mu_v, customizable
out_pars[name] = xu.with_unit(
pm.Normal(name, 0.,
sigma_v[name].value),
sigma_v[name].unit)
for k in pars.keys():
out_pars[k] = pars[k]
return out_pars
def __init__(self, P, e, sigma_K0, P0, mu=0., max_K=500*u.km/u.s,
K_unit=None, **kwargs):
self._sigma_K0 = sigma_K0
self._P0 = P0
self._max_K = max_K
if K_unit is not None:
self._sigma_K0 = self.sigma_K0.to(K_unit)
self._max_K = self._max_K.to(self._sigma_K0.unit)
if hasattr(P, xu.UNIT_ATTR_NAME):
self._P0 = self._P0.to(getattr(P, xu.UNIT_ATTR_NAME))
sigma_K0 = self._sigma_K0.value
P0 = self._P0.value
sigma = tt.min([self._max_K.value,
sigma_K0 * (P/P0)**(-1/3) / np.sqrt(1-e**2)])
super().__init__(mu=mu, sigma=sigma)
self._nonlinear_equiv_units = get_nonlinear_equiv_units()
self._linear_equiv_units = get_linear_equiv_units(self.poly_trend)
self._v0_offsets_equiv_units = get_v0_offsets_equiv_units(self.n_offsets)
self._all_par_unit_equiv = {**self._nonlinear_equiv_units,
**self._linear_equiv_units,
**self._v0_offsets_equiv_units}
# At this point, pars must be a dictionary: validate that all
# parameters are specified and that they all have units
for name in self.par_names:
if name not in pars:
raise ValueError(f"Missing prior for parameter '{name}': "
"you must specify a prior distribution for "
"all parameters.")
if not hasattr(pars[name], xu.UNIT_ATTR_NAME):
raise ValueError(f"Parameter '{name}' does not have associated "
"units: Use exoplanet.units to specify units "
"for your pymc3 variables. See the "
"documentation for examples: thejoker.rtfd.io")
equiv_unit = self._all_par_unit_equiv[name]
if not getattr(pars[name],
xu.UNIT_ATTR_NAME).is_equivalent(equiv_unit):
raise ValueError(f"Parameter '{name}' has an invalid unit: "
f"The units for this parameter must be "
f"transformable to '{equiv_unit}'")
# Enforce that the priors on all linear parameters are Normal (or a
# subclass of Normal)
for name in (list(self._linear_equiv_units.keys())
+ list(self._v0_offsets_equiv_units.keys())):
df = pm.trace_to_dataframe(trace)
data, *_ = validate_prepare_data(data,
self.prior.poly_trend,
self.prior.n_offsets)
samples = JokerSamples(poly_trend=self.prior.poly_trend,
n_offsets=self.prior.n_offsets,
t0=data.t0)
if names is None:
names = self.prior.par_names
for name in names:
par = self.prior.pars[name]
unit = getattr(par, xu.UNIT_ATTR_NAME)
samples[name] = df[name].values * unit
return samples
# First, prepare the joker_samples:
if not isinstance(joker_samples, JokerSamples):
raise TypeError("You must pass in a JokerSamples instance to the "
"joker_samples argument.")
if len(joker_samples) > 1:
# check if unimodal in P, if not, warn
if not is_P_unimodal(joker_samples, data):
logger.warn("TODO: samples ain't unimodal")
joker_samples = joker_samples.median()
mcmc_init = dict()
for name in self.prior.par_names:
unit = getattr(self.prior.pars[name], xu.UNIT_ATTR_NAME)
if joker_samples[name].shape == ():
mcmc_init[name] = joker_samples[name].to_value(unit)
else:
mcmc_init[name] = joker_samples[name].to_value(unit)[0]
if 't_peri' in model.named_vars and use_cached:
logger.debug('pymc3 model has already been setup for running MCMC '
'- using the previously setup model parameters.')
return mcmc_init
# remove previously-added parameters:
names = ['t_peri', 'obs']
to_remove = []
for name in names:
for i, par in enumerate(model.vars):
if par.name == name:
def __init__(self, P, e, sigma_K0, P0, mu=0., max_K=500*u.km/u.s,
K_unit=None, **kwargs):
self._sigma_K0 = sigma_K0
self._P0 = P0
self._max_K = max_K
if K_unit is not None:
self._sigma_K0 = self.sigma_K0.to(K_unit)
self._max_K = self._max_K.to(self._sigma_K0.unit)
if hasattr(P, xu.UNIT_ATTR_NAME):
self._P0 = self._P0.to(getattr(P, xu.UNIT_ATTR_NAME))
sigma_K0 = self._sigma_K0.value
P0 = self._P0.value
sigma = tt.min([self._max_K.value,
sigma_K0 * (P/P0)**(-1/3) / np.sqrt(1-e**2)])
super().__init__(mu=mu, sigma=sigma)
# parameters are specified and that they all have units
for name in self.par_names:
if name not in pars:
raise ValueError(f"Missing prior for parameter '{name}': "
"you must specify a prior distribution for "
"all parameters.")
if not hasattr(pars[name], xu.UNIT_ATTR_NAME):
raise ValueError(f"Parameter '{name}' does not have associated "
"units: Use exoplanet.units to specify units "
"for your pymc3 variables. See the "
"documentation for examples: thejoker.rtfd.io")
equiv_unit = self._all_par_unit_equiv[name]
if not getattr(pars[name],
xu.UNIT_ATTR_NAME).is_equivalent(equiv_unit):
raise ValueError(f"Parameter '{name}' has an invalid unit: "
f"The units for this parameter must be "
f"transformable to '{equiv_unit}'")
# Enforce that the priors on all linear parameters are Normal (or a
# subclass of Normal)
for name in (list(self._linear_equiv_units.keys())
+ list(self._v0_offsets_equiv_units.keys())):
if not isinstance(pars[name].distribution, pm.Normal):
raise ValueError("Priors on the linear parameters (K, v0, "
"etc.) must be independent Normal "
"distributions, not '{}'"
.format(type(pars[name].distribution)))
self.pars = pars