Quick Start

This guide walks through the core hydra_tod workflow: simulating time-ordered data, running the full Gibbs sampler, assessing convergence, and using individual sub-samplers in isolation.

For complete annotated examples, see the Jupyter notebooks in the examples/ directory of the repository.

1. Simulating a Single TOD

TODSimulation generates a synthetic MeerKAT-like observation — scan pattern, beam projection, gain variations, noise-diode injections, and correlated \(1/f\) noise — in one call.

from hydra_tod.simulation import TODSimulation

sim = TODSimulation(
    nside=64,              # HEALPix resolution
    elevation=41.5,        # Telescope elevation (degrees)
    freq=750,              # Observing frequency (MHz)
    alpha=2.0,             # Flicker noise spectral index
    logf0=-4.87,           # Log10 of knee frequency (rad/s)
    ptsrc_path="gleam_nside512_K_allsky_408MHz.npy",
)

tod    = sim.TOD_setting       # Observed data vector
sky    = sim.sky_params        # True sky pixel temperatures
gains  = sim.gains_setting     # True gain time series
noise  = sim.noise_setting     # True 1/f noise realisation

The simulation also creates the linear operators needed by the sampler:

gain_proj    = sim.gain_proj       # Legendre gain projection matrix
Tsky_op      = sim.Tsky_operator   # Beam projection (sky → TOD)
Tloc_op      = sim.Tloc_operator   # Local-temp projection (rec + nd)
t_list       = sim.time_list       # Observation times (seconds)
logfc        = sim.logfc           # Log cutoff frequency

2. Running the Full Gibbs Sampler

TOD_Gibbs_sampler() takes the operators from the simulation and returns posterior samples for all parameter blocks.

import numpy as np
from hydra_tod.full_Gibbs_sampler import TOD_Gibbs_sampler

# --- priors ---
n_sky   = sim.sky_params.shape[0]
n_loc   = sim.Tloc_operator.shape[1]
n_gain  = sim.gain_proj.shape[1]

# 20 % Gaussian prior on sky pixels
prior_sky_cov_inv  = (1.0 / (0.2 * sim.sky_params)) ** 2
prior_sky_mean     = sim.sky_params          # prior centred on truth

# 10 % prior on local temperature and gain coefficients
prior_loc_cov_inv  = np.ones(n_loc) * (10.0) ** 2
prior_loc_mean     = np.zeros(n_loc)
prior_gain_cov_inv = np.ones(n_gain) * (10.0) ** 2
prior_gain_mean    = np.zeros(n_gain)

# --- initial conditions ---
init_sky  = sim.sky_params * 1.05           # slight offset from truth
init_loc  = np.zeros(n_loc)
init_noise = [sim.logf0, sim.alpha]

# --- run sampler ---
Tsky_samples, gain_samples, noise_samples, Tloc_samples = TOD_Gibbs_sampler(
    local_TOD_list=[tod],
    local_t_lists=[t_list],
    local_gain_operator_list=[gain_proj],
    local_Tsky_operator_list=[Tsky_op],
    local_Tloc_operator_list=[Tloc_op],
    init_Tsky_params=init_sky,
    init_Tloc_params_list=[init_loc],
    init_noise_params_list=[init_noise],
    local_logfc_list=[logfc],
    prior_cov_inv_Tsky=prior_sky_cov_inv,
    prior_mean_Tsky=prior_sky_mean,
    prior_cov_inv_Tloc_list=[prior_loc_cov_inv],
    prior_mean_Tloc_list=[prior_loc_mean],
    prior_cov_inv_gain_list=[prior_gain_cov_inv],
    prior_mean_gain_list=[prior_gain_mean],
    n_samples=2000,
    gain_model="linear",
    sampler="emcee",
    jeffreys=True,
)

# Tsky_samples  shape: (2000, n_sky)
# gain_samples  shape: (2000, n_gain)
# noise_samples shape: (2000, 2)   — (logf0, alpha) per iteration
# Tloc_samples  shape: (2000, n_loc)

Note

For multi-TOD analysis with MPI, launch the script with mpiexec -n N python script.py. Each rank receives a slice of local_TOD_list and the sampler handles synchronisation automatically.

3. Inspecting Posterior Samples

Use the built-in diagnostics to check convergence before interpreting results.

from hydra_tod.mcmc_diagnostics import diagnostics
import matplotlib
matplotlib.use("Agg")

# Wrap single chain as shape (1, n_iter, n_params)
noise_chains = noise_samples[np.newaxis, :, :]   # (1, 2000, 2)
summary = diagnostics(
    noise_chains,
    param_names=["logf0", "alpha"],
    max_plots=2,      # produce trace + ACF plots
)

print(summary["logf0"])
# {"ESS_min": 342.1, "ESS_median": 342.1, "Rhat_split": 1.002}

A well-converged chain should have ESS > ~100 and split-\(\hat{R}\) close to 1.0.

Posterior summaries for sky maps:

sky_mean = Tsky_samples.mean(axis=0)          # posterior mean
sky_std  = Tsky_samples.std(axis=0)           # posterior std
residual = sky_mean - sim.sky_params           # bias map
z_score  = residual / sky_std                 # should be ~ N(0,1)

4. Using Sub-samplers Directly

Each Gibbs step can be called in isolation, which is useful for debugging or for building custom inference pipelines.

Sample gains only

from hydra_tod.gain_sampler import gain_sampler
from hydra_tod.flicker_model import flicker_cov
from scipy.linalg import toeplitz
import numpy as np

# Build noise covariance inverse from current noise params
from hydra_tod.utils import lag_list, cho_compute_mat_inv
from hydra_tod.flicker_model import flicker_cov_vec

lags     = lag_list(t_list)
cov_row  = flicker_cov_vec(lags, 10**sim.logf0, 10**sim.logfc, sim.alpha)
Ncov_inv = cho_compute_mat_inv(toeplitz(cov_row))

gain_coeffs, gains = gain_sampler(
    data=tod,
    t_list=t_list,
    gain_proj=gain_proj,
    Tsys=sky[..., np.newaxis] * np.ones_like(tod),  # current Tsys estimate
    noise_params=(sim.logf0, sim.alpha),
    logfc=sim.logfc,
    model="linear",
    prior_cov_inv=prior_gain_cov_inv,
    prior_mean=prior_gain_mean,
)

Sample noise parameters only

from hydra_tod.noise_sampler_fixed_fc import flicker_sampler

noise_sample = flicker_sampler(
    TOD=tod,
    gains=gains,          # current gain estimate
    Tsys=tsys_estimate,   # current Tsys estimate
    init_params=[sim.logf0, sim.alpha],
    n_samples=1,
    jeffreys=True,
    sampler="emcee",
    consecutive=True,     # set False for flagged/non-contiguous data
)
logf0_new, alpha_new = noise_sample