# Licensed under a 3-clause BSD style license - see LICENSE.rst
import platform
import numpy as np
import arviz_stats
import pymc as pm
import astropy.units as u
# conda-forge adds their own clang++ to pymc dependencies, but it conflicts with
# the xcode command-line tools libraries. enforce the use of apple's clang++.
if platform.system() == "Darwin":
import pytensor
pytensor.config.cxx = "/usr/bin/clang++"
__all__ = ["azel_fit", "best_fit_pars"]
DEG2RAD = np.pi / 180
AZEL_TERMS = ("ia", "ie", "an", "aw", "ca", "npae", "tf", "tx")
HADEC_TERMS = ("ih", "id", "np", "ch", "ma", "me", "tf")
AZIMUTH_FUNCS = {
"ia": lambda az, el: -1.0,
"an": lambda az, el: -1.0 * pm.math.sin(DEG2RAD * az) * pm.math.tan(DEG2RAD * el),
"aw": lambda az, el: -1.0 * pm.math.cos(DEG2RAD * az) * pm.math.tan(DEG2RAD * el),
"ca": lambda az, el: -1.0 / pm.math.cos(DEG2RAD * el),
"npae": lambda az, el: -1.0 * pm.math.tan(DEG2RAD * el),
}
ELEVATION_FUNCS = {
"ie": lambda az, el: 1.0,
"an": lambda az, el: -1.0 * pm.math.cos(DEG2RAD * az),
"aw": lambda az, el: pm.math.sin(DEG2RAD * az),
"tf": lambda az, el: -1.0 * pm.math.cos(DEG2RAD * el),
"tx": lambda az, el: -1.0 / pm.math.tan(DEG2RAD * el),
}
[docs]
def azel_fit(
coo_ref,
coo_meas,
nsamp=2000,
ntune=500,
target_accept=0.95,
random_seed=8675309,
cores=None,
fit_terms=AZEL_TERMS,
fixed_terms=None,
init_pars=None,
prior_sigmas=None,
):
"""
Fit full az/el pointing model using PyMC. The terms are analogous to those used by
TPOINT(tm). This fit includes the eight normal terms used and described in
`~pytelpoint.transform.azel` with additional terms, az_sigma and el_sigma, that
describe the intrinsic/observational scatter.
Parameters
----------
coo_ref : `~astropy.coordinates.SkyCoord` instance
Reference coordinates
coo_meas : `~astropy.coordinates.SkyCoord` instance
Measured coordinates
nsamp : int (default: 2000)
Number of inference samples per chain
ntune : int (default: 500)
Number of burn-in samples per chain
target_accept : float (default: 0.95)
Sets acceptance probability target for determining step size
random_seed : int (default: 8675309)
Seed number for random number generator
cores : int (default: None)
Number of cores to use for parallel chains. The default of None
will use the number of available cores, but no more than 4.
fit_terms : list-like (default: AZEL_TERMS)
Model terms to include in the fit.
fixed_terms : dict (default: {})
Dict of terms to fix to a specified value.
init_pars : dict (default: None -> {'ia': 1200.})
Initial guesses for the fit parameters. Keys are the same those provided by
`~pytelpoint.fitting.best_fit_pars` and described in
`~pytelpoint.transform.azel`: 'ia', 'ie', 'an', 'aw', 'ca', 'npae', 'tf', 'tx',
'az_sigma', 'el_sigma'. The default for 'ia' is appropriate for the MMTO.
If not specified, then the initial guess for a parameter is assumed to be 0.
prior_sigmas : dict (default: None -> {'ia': 100., 'ie': 50.})
The priors for the fit parameters are assumed to be `~pymc.Normal`
distributions. The sigmas for these can be specified here. The index
parameters, 'ia' and 'ie', have default sigma values of 100 and
50, respectively. The rest default to 25 if not specified.
Returns
-------
idata : `~arviz.InferenceData`
Inference data from the pointing model
"""
if fixed_terms is None:
fixed_terms = {}
if init_pars is None:
init_pars = {"ia": 1200.0}
if prior_sigmas is None:
prior_sigmas = {"ia": 100.0, "ie": 50.0}
pointing_model = pm.Model()
with pointing_model:
# az/el are the astrometric reference values. az_raw/el_raw are the observed
# encoder values. they should be in degrees, but are converted here
# just in case.
az = pm.Data("az", coo_ref.az.to(u.deg).value)
el = pm.Data("el", coo_ref.alt.to(u.deg).value)
az_raw = pm.Data("az_raw", coo_meas.az.to(u.deg).value)
el_raw = pm.Data("el_raw", coo_meas.alt.to(u.deg).value)
terms = {}
combined_terms = list(fit_terms) + list(fixed_terms.keys())
for term in combined_terms:
if term not in AZEL_TERMS:
raise ValueError(f"Invalid az/el fitting term, {term}.")
if term in fixed_terms:
terms[term] = fixed_terms[term]
else:
terms[term] = pm.Normal(
term, init_pars.get(term, 0.0), prior_sigmas.get(term, 25.0)
)
az_sigma = pm.HalfNormal("az_sigma", sigma=init_pars.get("az_sigma", 1.0))
el_sigma = pm.HalfNormal("el_sigma", sigma=init_pars.get("el_sigma", 1.0))
daz = 0.0
for k, f in AZIMUTH_FUNCS.items():
if k in combined_terms:
daz += terms[k] * f(az, el)
# using 'dalt' because 'del' is a python built-in
dalt = 0.0
for k, f in ELEVATION_FUNCS.items():
if k in combined_terms:
dalt += terms[k] * f(az, el)
# models are the raw encoder values plus pointing model; observed are
# the actual az/el
mu_az = az_raw + daz / 3600.0
mu_el = el_raw + dalt / 3600.0
_ = pm.Normal("azerr", mu=mu_az, sigma=az_sigma / 3600, observed=az)
_ = pm.Normal("elerr", mu=mu_el, sigma=el_sigma / 3600, observed=el)
idata = pm.sample(
nsamp,
tune=ntune,
target_accept=target_accept,
return_inferencedata=True,
random_seed=random_seed,
cores=cores,
)
return idata
[docs]
def best_fit_pars(idata):
"""
Pull out the best fit parameters from a pointing model fit and
return them as a dict.
Parameters
----------
idata : `~arviz.InferenceData`
Inference data from pointing model.
Returns
-------
pointing_pars : dict
Best-fit pointing parameters
"""
t_fit = arviz_stats.summary(idata, round_to=8)
pointing_pars = {}
for p in t_fit.index:
pointing_pars[p] = t_fit.loc[p, "mean"]
return pointing_pars