API reference¶
Auto-generated from the source docstrings. The public API is the symbols below; treat everything else as internal.
Stellar populations¶
ceridwen.ssps.SSPData
dataclass
¶
SSPData(
ssp_lgmet,
ssp_lg_age_gyr,
ssp_wave,
ssp_flux,
isoc_type=None,
spec_library=None,
imf_type=None,
fsps_version=None,
fsps_kwargs=dict(),
wave_min=None,
wave_max=None,
schema_version=None,
)
Immutable container for the SSP interpolation grids (+ provenance).
Attributes:
| Name | Type | Description |
|---|---|---|
ssp_lgmet |
(ndarray, shape(n_met))
|
|
ssp_lg_age_gyr |
(ndarray, shape(n_ages))
|
|
ssp_wave |
(ndarray, shape(n_wave))
|
Wavelength grid in Angstroms. |
ssp_flux |
(ndarray, shape(n_met, n_ages, n_wave))
|
SSP flux density in |
isoc_type |
str or None
|
Isochrone library the grid was built with (e.g. |
spec_library |
str or None
|
Spectral library (e.g. |
imf_type |
int or None
|
FSPS IMF selector the grid was built with. |
fsps_version |
str or None
|
|
fsps_kwargs |
dict
|
The (whitelisted) FSPS build kwargs actually used. |
wave_min, wave_max |
float or None
|
Wavelength range (Å) of |
schema_version |
str or None
|
On-disk metadata schema tag. |
Notes
The provenance fields are ordinary Python objects (str / int / dict /
None); they are never JAX arrays and never enter a @jit kernel.
They are excluded from equality/hashing (compare=False).
No log_qq table is stored; the nebular model computes the ionising
photon rate internally. HDF5 files that contain a log_qq dataset
are loaded transparently — the field is simply ignored.
__post_init__ ¶
Validate grid consistency.
Source code in ceridwen/ssps/ssp_data.py
display ¶
Print a summary of the grid and its provenance.
Intended as a sanity check: call it right after from_fsps or
load to confirm the isochrone set, spectral library, IMF, and
grid coverage are what you expect before you build a CSPBasis.
Purely diagnostic — none of this is touched by the forward model.
Parameters:
| Name | Type | Description | Default |
|---|---|---|---|
return_str
|
bool
|
Return the formatted string instead of printing it. Default False. |
False
|
file
|
file - like
|
Destination for the print (default |
None
|
Returns:
| Type | Description |
|---|---|
str or None
|
The formatted string if |
Source code in ceridwen/ssps/ssp_data.py
save ¶
Serialise the SSP grids (and provenance metadata) to HDF5.
Provenance fields are written to the file attrs; any that are
None are simply omitted, so re-saving a legacy grid does not
invent metadata. fsps_kwargs is stored as a JSON string.
Parameters:
| Name | Type | Description | Default |
|---|---|---|---|
filename
|
str or Path
|
Output file path. Will be overwritten if it exists. |
required |
Source code in ceridwen/ssps/ssp_data.py
load
classmethod
¶
Load an :class:SSPData from an HDF5 file.
Backward compatible: files written before provenance tracking carry
none of the metadata attrs, so every provenance field is populated
as None / {} rather than raising. Any legacy log_qq
dataset is silently ignored — the nebular model computes its own
ionising-photon rate from ssp_flux.
Source code in ceridwen/ssps/ssp_data.py
from_fsps
classmethod
¶
Build an :class:SSPData directly from FSPS, recording provenance.
Only kwargs that define the stellar library / IMF are accepted —
the things that legitimately belong at SSP-build time. Anything the
CSP forward model applies itself (star-formation history, dust,
nebular emission, IGM, redshift, LOSVD smoothing, or a fixed
metallicity) raises :class:ValueError, so the grid can never be
silently double-processed or made inconsistent with
:class:ceridwen.csp.CSPBasis. zcontinuous and sfh are
fixed internally (the grid is built on FSPS's discrete zlegend
metallicity points).
Parameters:
| Name | Type | Description | Default |
|---|---|---|---|
save_to
|
str or Path
|
If given, the result is also persisted to this path via
:meth: |
None
|
**fsps_kwargs
|
Forwarded to :class: |
{}
|
Raises:
| Type | Description |
|---|---|
ValueError
|
If any kwarg is not a library/IMF-defining parameter. |
Source code in ceridwen/ssps/ssp_data.py
Composite stellar population (forward model)¶
ceridwen.csp.CSPBasis ¶
CSPBasis(
SSPData,
theta=None,
tuniv=13.8,
tiny_logt=-70,
zh_const=False,
add_neb=True,
init_neb_params=None,
nebemlineinspec=False,
add_dust=True,
add_diffuse_dust=True,
add_dust_emission=False,
add_igm=False,
igm_model="madau1995",
igm_factor=1.0,
sps_home=None,
init_dust_params=None,
diffuse_law="kriek_conroy",
verbose=True,
sfh_interp="step",
sigma_losvd_kms=300.0,
track_zred_age=False,
lookback_time=None,
sfh_per_bin=False,
**kwargs
)
Composite Stellar Population basis using a dict-valued theta.
The public interface is identical to csp.CSPBasis except that
predict(theta) now expects (and theta_init now is) a
dict[str, Array] rather than a flat 1-D jnp.ndarray.
Parameters:
| Name | Type | Description | Default |
|---|---|---|---|
SSPData
|
SSPData
|
Frozen dataclass with SSP grids (wave, flux, ages, zmet). The
nebular model computes its own ionising-photon rate from
|
required |
theta
|
dict
|
Initial parameter values. Must contain
The construction-time grid is a default, not a straitjacket: an
explicit Lookback-time convention
|
None
|
tuniv
|
float
|
Age of the Universe in Gyr. Default 13.8. |
13.8
|
zh_const
|
bool
|
If True, use constant metallicity (requires key |
False
|
add_neb
|
bool
|
Physics switches. |
True
|
add_dust
|
bool
|
Physics switches. |
True
|
add_diffuse_dust
|
bool
|
Physics switches. |
True
|
add_dust_emission
|
bool
|
Physics switches. |
True
|
sps_home
|
str
|
Path to the FSPS data directory (needed for nebular and dust-emission
grid loading). Defaults to the |
None
|
init_neb_params
|
dict
|
Keyword arguments forwarded to |
None
|
init_dust_params
|
dict
|
Keyword arguments forwarded to |
None
|
diffuse_law
|
str
|
Attenuation law name for the diffuse dust component. |
'kriek_conroy'
|
verbose
|
bool
|
Print parameter summary after initialization. |
True
|
lookback_time
|
array - like
|
Shortcut alternative to Example:: |
None
|
sfh_per_bin
|
bool
|
Only used with the |
False
|
sfh_interp : {'step', 'linear'} Controls the SFH integration scheme used when computing SSP weights.
``'step'`` (default) — piecewise-constant (FastStepBasis-style).
The SFR is held at the mean of the two endpoint values within
each SFH time bin. The weight of each SSP age bin is the
product of that constant SFR and the linear-time overlap between
the SFH bin and the SSP age bin. Weights are non-negative by
construction — no clipping is ever needed.
``'linear'`` — piecewise-linear.
Analytically integrates a linearly-interpolated SFH against the
SSP age bins in log-age space (``intsfwght``). Higher-order
accurate, but can produce small negative weights for steep SFH
gradients, which are then clipped.
To switch at runtime::
csp.calculate_ssp_weights = csp.calculate_ssp_weights_const_zh_step
# or
csp.calculate_ssp_weights = csp.calculate_ssp_weights_const_zh
Source code in ceridwen/csp/csp.py
250 251 252 253 254 255 256 257 258 259 260 261 262 263 264 265 266 267 268 269 270 271 272 273 274 275 276 277 278 279 280 281 282 283 284 285 286 287 288 289 290 291 292 293 294 295 296 297 298 299 300 301 302 303 304 305 306 307 308 309 310 311 312 313 314 315 316 317 318 319 320 321 322 323 324 325 326 327 328 329 330 331 332 333 334 335 336 337 338 339 340 341 342 343 344 345 346 347 348 349 350 351 352 353 354 355 356 357 358 359 360 361 362 363 364 365 366 367 368 369 370 371 372 373 374 375 376 377 378 379 380 381 382 383 384 385 386 387 388 389 390 391 392 393 394 395 396 397 398 399 400 401 402 403 404 405 406 407 408 409 410 411 412 413 414 415 416 417 418 419 420 421 422 423 424 425 426 427 428 429 430 431 432 433 434 435 436 437 438 439 440 441 442 443 444 445 446 447 448 449 450 451 452 453 454 455 456 457 458 459 460 461 462 463 464 465 466 467 468 469 470 471 472 473 474 475 476 477 478 479 480 481 482 483 484 485 486 487 488 489 490 491 492 493 494 495 496 497 498 499 500 | |
initialize_model_structure ¶
Validate the incoming theta and store self.theta_init.
The dict is validated, converted to JAX arrays, and stored directly.
Required keys
"sfh" : shape (n_time,)
"lookback_time" : shape (n_time,)
Either "Z" (scalar, constant metallicity) or "zh" (shape
(n_time,), time-varying metallicity) must be present, depending on
the zh_const flag set during __init__.
Source code in ceridwen/csp/csp.py
504 505 506 507 508 509 510 511 512 513 514 515 516 517 518 519 520 521 522 523 524 525 526 527 528 529 530 531 532 533 534 535 536 537 538 539 540 541 542 543 544 545 546 547 548 549 550 551 552 553 554 555 556 557 558 559 560 561 562 563 564 565 566 567 568 569 570 571 572 573 574 575 576 577 578 579 580 581 582 583 584 585 586 587 588 589 590 591 592 593 594 595 596 597 598 599 600 601 602 603 604 605 606 607 608 609 610 611 612 613 614 615 616 617 618 619 620 621 622 623 624 625 626 627 628 629 630 631 632 633 634 635 636 637 638 639 640 641 642 643 644 645 646 647 648 649 650 651 652 653 654 655 656 657 658 659 660 661 662 663 664 665 666 667 668 669 670 671 672 673 674 675 676 677 678 679 680 681 682 683 684 685 686 687 | |
register_known_theta_keys ¶
Register additional recognized theta keys so they are not mis-flagged
as typos by :meth:_warn_unknown_theta_keys.
SedModel calls this with its model-level free parameters (e.g.
logsfr_ratios, which the sfh transform consumes): those keys are
forwarded through to predict in the full theta dict but are not CSP
parameters, so without this they would trigger a spurious typo warning.
Source code in ceridwen/csp/csp.py
check_param_ranges ¶
Diagnostic (NON-jitted): list parameters that fall outside the interpolation grids, where the model silently clamps to the nearest grid edge and thus hides extrapolation.
Intended to be called once on your theta (or theta bounds) before a
fit; it is never invoked from the hot path. Returns the list of
human-readable messages (and emits them as warnings when warn).
Source code in ceridwen/csp/csp.py
set_attenuation_function ¶
Build and assign self.attenuate_dust(wave, theta) → (attn, attn_diffuse).
With the dict theta, each dust model simply reads the keys it knows about from the shared theta dict. No NamedTuple construction needed.
Source code in ceridwen/csp/csp.py
get_spectrum_components ¶
Return the canonical (continuum, lines) line decomposition.
Both arrays are on the rest-frame model grid self.wave and are
unscaled -- mass / redshift / IGM factors are applied downstream
by predict and get_line_spec, exactly as for get_spectrum.
continuum-- the line-free spectrum (stellar continuum + nebular continuum), dust-attenuated. Identical toget_spectrum(theta, include_lines=False).lines-- the broadened nebular emission-line component alone, carried through the same dust attenuation, so the full SED is recovered ascontinuum + lines.
This is the single source of truth for "spectrum with vs. without
emission lines": predict builds the photometry and slit spectra
from it and get_line_spec returns its lines term.
nebemlineinspec does not affect this method -- it only sets the
default of the single-array public get_spectrum. With
add_neb=False there is no nebular module and lines is
identically zero.
Source code in ceridwen/csp/csp.py
predict ¶
Compute the CSP spectrum and project it onto every observation.
This method is the primary hot-path entry point for the sampler.
It is designed to be fully JAX JIT-compatible with zero Python
if / isinstance branches in the traced code path:
get_spectrum(theta)is pure JAX.- The Python
forloop overobservationsis unrolled at trace time becauseobservationsis a static Python list (part of the closure, not a traced argument). obs.predict(spectrum, self.wave)dispatches through Python's method resolution order (static at trace time) to the appropriate subclass implementation — either a dense matrix–vector multiply (Spectrum,Lines) or a filter-set convolution (Photometry). The XLA kernel contains no conditional branches.
Pre-condition: every Observation in observations must have
had obs.setup_for_model(self.wave) called before the first JIT
trace. SedModel.__init__ does this automatically.
For a raw model spectrum without projection, use get_spectrum(theta)
directly; for the separate line-free continuum and emission-line
component, use get_spectrum_components(theta).
Parameters:
| Name | Type | Description | Default |
|---|---|---|---|
theta
|
dict[str, Array]
|
Free-parameter dict. Must contain at minimum |
required |
observations
|
list of Observation
|
Observations to project onto. Must be the same Python objects (same list structure, same types) on every call — changing the list forces a retrace. |
required |
Returns:
| Name | Type | Description |
|---|---|---|
predictions |
dict[str, Array]
|
Keyed by
|
.. warning::
|
The outputs are observed-frame AB maggies only if |
Source code in ceridwen/csp/csp.py
1092 1093 1094 1095 1096 1097 1098 1099 1100 1101 1102 1103 1104 1105 1106 1107 1108 1109 1110 1111 1112 1113 1114 1115 1116 1117 1118 1119 1120 1121 1122 1123 1124 1125 1126 1127 1128 1129 1130 1131 1132 1133 1134 1135 1136 1137 1138 1139 1140 1141 1142 1143 1144 1145 1146 1147 1148 1149 1150 1151 1152 1153 1154 1155 | |
get_line_spec ¶
Return the broadened-emission-line component of the model spectrum.
Computed as get_spectrum(include_lines=True) - get_spectrum(include_lines=False),
which gives the contribution of the nebular lines alone -- the
prospector-style "line-only" spectrum. Mass + redshift + IGM
scaling are applied identically to CSPBasis.predict so the
output is at the same physical scale as the observation arrays.
add_neb=False makes this return zero (no nebular module).
Source code in ceridwen/csp/csp.py
display_sfh ¶
display_sfh(
theta=None,
ax=None,
*,
overlay_nodes=True,
show_bin_edges=False,
units="Gyr",
**plot_kwargs
)
Plot the SFH against lookback time, rendered identically to the
interpretation used by :meth:_ssp_weights.
For sfh_interp == "step" this draws a piecewise-constant function
with one horizontal segment per bin [T_{i+1}, T_i] at height
:math:\bar\psi_i (the per-bin SFR consumed by
calculate_ssp_weights_*_step). For sfh_interp == "linear" it
draws the piecewise-linear interpolant between per-node SFR values --
the same function whose analytic integral against the SSP age grid
is computed by intsfwght.
Lookback time is read from theta["lookback_time"] if supplied
(units: Gyr) and otherwise falls back to self.sfh_times (which
is stored in years and converted back to Gyr here). theta_init
intentionally does NOT carry lookback_time -- it is a static
grid, not a free parameter -- so the default fallback path is the
common case.
The x-axis runs left-to-right in increasing lookback time: present
day (T = 0) sits at the origin on the left, and the oldest sampled
node sits on the right. This matches the natural index order of
theta["lookback_time"].
Parameters:
| Name | Type | Description | Default |
|---|---|---|---|
theta
|
dict
|
Parameter dict to display. Defaults to |
None
|
ax
|
Axes
|
Axes to draw into. If None, a new figure is created. |
None
|
overlay_nodes
|
bool
|
If True, mark per-bin SFR values at bin midpoints (step mode) or per-node SFR values at lookback nodes (linear mode). |
True
|
show_bin_edges
|
bool
|
If True, draw vertical dotted lines at every node |
False
|
units
|
(Gyr, yr, Myr)
|
X-axis units for the lookback-time axis. The SFR axis is always [M_sun / yr]. |
"Gyr"
|
**plot_kwargs
|
Forwarded to the per-segment |
{}
|
Returns:
| Name | Type | Description |
|---|---|---|
ax |
Axes
|
The axes containing the plot. |
Raises:
| Type | Description |
|---|---|
AssertionError
|
If the per-bin integral of the displayed SFR disagrees with the
per-bin mass |
Source code in ceridwen/csp/csp.py
1350 1351 1352 1353 1354 1355 1356 1357 1358 1359 1360 1361 1362 1363 1364 1365 1366 1367 1368 1369 1370 1371 1372 1373 1374 1375 1376 1377 1378 1379 1380 1381 1382 1383 1384 1385 1386 1387 1388 1389 1390 1391 1392 1393 1394 1395 1396 1397 1398 1399 1400 1401 1402 1403 1404 1405 1406 1407 1408 1409 1410 1411 1412 1413 1414 1415 1416 1417 1418 1419 1420 1421 1422 1423 1424 1425 1426 1427 1428 1429 1430 1431 1432 1433 1434 1435 1436 1437 1438 1439 1440 1441 1442 1443 1444 1445 1446 1447 1448 1449 1450 1451 1452 1453 1454 1455 1456 1457 1458 1459 1460 1461 1462 1463 1464 1465 1466 1467 1468 1469 1470 1471 1472 1473 1474 1475 1476 1477 1478 1479 1480 1481 1482 1483 1484 1485 1486 1487 1488 1489 1490 1491 1492 1493 1494 1495 1496 1497 1498 1499 1500 1501 1502 1503 1504 1505 1506 1507 1508 1509 1510 1511 1512 1513 1514 1515 1516 1517 1518 1519 1520 1521 1522 1523 1524 1525 1526 1527 1528 1529 1530 1531 1532 1533 1534 1535 1536 1537 1538 1539 1540 1541 1542 1543 1544 1545 1546 1547 1548 1549 | |
calculate_ssp_weights_const_zh ¶
Constant-metallicity, piecewise-linear SFH weights.
Thin wrapper over :meth:_ssp_weights; reads theta["sfh"]
(shape (n_time,), linear SFR) and theta["Z"] (shape (1,),
log10 absolute metallicity on the self.zmet / ssp_lgmet grid —
NOT log10 Z/Zsun). Same units as the var-zh variants' theta["zh"].
Source code in ceridwen/csp/csp.py
calculate_ssp_weights_const_zh_step ¶
Constant-metallicity, piecewise-constant (FastStepBasis-style) SFH
weights. Thin wrapper over :meth:_ssp_weights. Reads
theta["sfh"] and theta["Z"].
Source code in ceridwen/csp/csp.py
calculate_ssp_weights_var_zh ¶
Time-varying-metallicity, piecewise-linear SFH weights.
Thin wrapper over :meth:_ssp_weights; reads theta["sfh"]
(shape (n_time,), linear SFR) and theta["zh"] (shape
(n_time,), log10 absolute metallicity at each lookback time, on the
self.zmet / ssp_lgmet grid — NOT log10 Z/Zsun). Identical units
to the const-zh variants' theta["Z"].
Source code in ceridwen/csp/csp.py
calculate_ssp_weights_var_zh_step ¶
Time-varying-metallicity, piecewise-constant (FastStepBasis-style)
SFH weights. Thin wrapper over :meth:_ssp_weights. Reads
theta["sfh"] and theta["zh"].
Source code in ceridwen/csp/csp.py
get_spectrum_dattn_nodem_neb ¶
Dust attenuation + nebular emission, no dust emission.
include_lines:
None (default) -> use self.nebemlineinspec.
True / False -> override (csp.predict always passes True).
Source code in ceridwen/csp/csp.py
get_spectrum_dattn_dem_neb ¶
Dust attenuation + nebular emission + dust emission.
Source code in ceridwen/csp/csp.py
get_spectrum_dattn_nodem_noneb ¶
Dust attenuation, no nebular, no dust emission.
include_lines is accepted but ignored: there is no nebular
emission to include / exclude when add_neb=False.
Source code in ceridwen/csp/csp.py
get_spectrum_dattn_dem_noneb ¶
Dust attenuation + dust emission, no nebular. include_lines accepted but ignored.
Source code in ceridwen/csp/csp.py
get_spectrum_nodattn_nodem_noneb ¶
Stellar continuum only — no dust, no nebular. include_lines ignored.
Source code in ceridwen/csp/csp.py
get_spectrum_nodattn_nodem_neb ¶
Nebular emission only, no dust.
Source code in ceridwen/csp/csp.py
Observations¶
ceridwen.observation.Photometry ¶
Bases: Observation
Broadband photometric observation in AB maggies.
Flux and uncertainty are stored in maggies (linear AB flux units;
1 maggie = 3631 Jy). Filter information is held in a
sedpy_jax.observate.FilterSet.
Parameters:
| Name | Type | Description | Default |
|---|---|---|---|
filters
|
list of str or list of Filter objects
|
Filters to include. Strings are resolved to |
[]
|
flux
|
(array - like, shape(n_filters))
|
Observed maggies. |
required |
uncertainty
|
(array - like, shape(n_filters))
|
1-sigma uncertainties in maggies. |
required |
mask
|
array-like of bool, shape (n_filters,)
|
True for filters that should be included in the fit. |
required |
Examples:
>>> phot = Photometry(
... filters=["sdss_u0", "sdss_g0", "sdss_r0", "sdss_i0", "sdss_z0"],
... flux=obs_maggies,
... uncertainty=obs_maggies_unc,
... )
>>> model_maggies = phot.get_maggies(model_wave, model_fnu)
>>> chi2 = phot.chi_sq(model_maggies)
Parameters:
| Name | Type | Description | Default |
|---|---|---|---|
upper_limit
|
array-like of bool, shape (n_filters,)
|
Per-band non-detection flags. If True for band .. math:: This mirrors the convention already used in
:class: |
None
|
Source code in ceridwen/observation/photometry.py
set_filters ¶
Set the filter list. filters may be a list of filter-name
strings or of sedpy_jax Filter objects.
Source code in ceridwen/observation/photometry.py
get_maggies ¶
Project a model spectrum onto the filters and return synthetic maggies.
The model spectrum is expected in F_nu units (e.g. L_sun Hz^{-1}
M_sun^{-1} as returned by CSPBasis.get_spectrum). Internally the
spectrum is converted to F_lambda before being projected through the
AB-normalised FilterSet transmission matrix, so the output has the
same relative normalisation as a standard AB photometric integral.
Parameters:
| Name | Type | Description | Default |
|---|---|---|---|
model_wave
|
(array - like, shape(n_wave))
|
Wavelength grid [Å]. |
required |
model_fnu
|
(array - like, shape(n_wave))
|
Model spectrum in F_nu units (L_sun/Hz/M_sun or erg/s/Hz/cm^2). |
required |
Returns:
| Name | Type | Description |
|---|---|---|
maggies |
(ndarray, shape(n_filters))
|
Synthetic photometry with the same relative normalisation as the input flux. |
Notes
The AB normalisation constant in sedpy_jax cancels dimensionally when both the Ceridwen and FSPS spectra are expressed in the same units, making model/data comparisons unit-independent.
Source code in ceridwen/observation/photometry.py
setup_for_model ¶
Precompute a (n_filters, n_wave) projection matrix _T so that
predict reduces to a single GEMV: maggies = _T @ F_nu.
The matrix folds together three operations that
FilterSet.get_sed_maggies does per call:
- F_nu -> F_lambda conversion:
F_lam = F_nu * c / lam^2 - Interpolation from the model wavelength grid onto the
FilterSet's internal grid (
interp_source) - Dot product with the precomputed
FilterSet.transmatrix
By composing these into a single static matrix _T of shape
(n_filters, n_wave_model), all three steps collapse into one
GEMV at predict time.
Must be called once before predict (and before JIT compilation).
SedModel.__init__ calls this automatically.
Parameters:
| Name | Type | Description | Default |
|---|---|---|---|
wave_model
|
array - like
|
Rest-frame wavelength grid of the model spectrum [Å]. |
required |
zred
|
float
|
Fixed redshift at which to precompute the filter projection.
Defaults to 0 (rest-frame). For a non-zero |
0.0
|
Source code in ceridwen/observation/photometry.py
158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200 201 202 203 204 205 206 207 208 209 210 211 212 213 214 215 216 217 218 219 220 221 222 223 224 225 226 227 228 229 230 231 232 233 234 235 236 237 238 239 240 241 242 | |
predict ¶
Project a model F_nu spectrum onto the filters.
If setup_for_model has been called, this is a single GEMV
(_T @ spectrum). Otherwise falls back to get_maggies.
Parameters:
| Name | Type | Description | Default |
|---|---|---|---|
spectrum
|
(Array, shape(n_wave))
|
Model spectrum in F_nu units. |
required |
wave_model
|
(Array, shape(n_wave))
|
Model wavelength grid [Å]. |
required |
Returns:
| Type | Description |
|---|---|
(Array, shape(n_filters))
|
Synthetic AB maggies. |
Source code in ceridwen/observation/photometry.py
predict_at_redshift ¶
Project an observer-frame F_nu spectrum through the filters when the redshift is a traced (sampled) JAX scalar.
This is the free-redshift counterpart of :meth:predict. The
GEMV fast path baked by :meth:setup_for_model assumes a single
Python-scalar zred was known at trace time and bakes the
observed-frame wavelength grid into the projection matrix
_T; that path cannot be used for sampling. Here the
observed-frame wavelength grid is reconstructed per-sample as
wave_obs = (1 + zred) * wave_rest and the spectrum is
projected via :meth:FilterSet.get_sed_maggies with the
traced sourcewave.
Pre-condition: spectrum_fnu_observed is the observer-frame
F_nu, i.e. CSPBasis.predict has already multiplied by
flux_factor_maggies(zred) and (optionally) by the IGM
transmission. This method only handles the wavelength-grid
bookkeeping and the filter integral.
Parameters:
| Name | Type | Description | Default |
|---|---|---|---|
spectrum_fnu_observed
|
(Array, shape(n_wave))
|
Observer-frame F_nu on the rest-frame model wavelength grid (the standard ceridwen output of get_spectrum + mass + flux-factor + IGM). |
required |
wave_rest
|
(Array, shape(n_wave))
|
Rest-frame model wavelength grid [Å] (typically |
required |
zred
|
(Array, scalar)
|
Sampled redshift. May be a traced array; the entire path
below is JIT-compatible and differentiable in |
required |
Returns:
| Type | Description |
|---|---|
(Array, shape(n_filters))
|
Synthetic AB maggies in observer frame. |
Notes
- Cost is one filter interpolation + one trans-matrix dot per sample, vs the single GEMV of the fixed-z path. For a 14-d NUTS / 4000-particle NS / 20 000-step SVI run on a 40 GB A100 this is ~10-20x slower than the GEMV but still saturates the GPU.
- Numerically equivalent to ``setup_for_model(wave_rest, zred=z)
- predict(spectrum, wave_rest)`` evaluated at the same z, to float32 precision.
- Works regardless of whether
setup_for_modelhas been called. When both paths are wired (e.g. for compare-mode plots), prefer the GEMV path for any fixed-z observation and this method for any free-z observation.
Source code in ceridwen/observation/photometry.py
chi_sq ¶
Chi-squared contribution from this photometric observation.
For bands flagged as upper limits (self.upper_limit[i] = True),
the contribution is one-sided: a penalty is applied only when the
model flux exceeds the observed upper-limit value. Matches the
convention used in :class:Lines and Prospector's recommended
treatment of non-detections.
Parameters:
| Name | Type | Description | Default |
|---|---|---|---|
model_maggies
|
(array - like, shape(n_filters))
|
Synthetic photometry on the same filter set. |
required |
Returns:
| Name | Type | Description |
|---|---|---|
chi2 |
float
|
|
Source code in ceridwen/observation/photometry.py
residuals ¶
Per-filter (data − model) / sigma. Masked filters are set to NaN.
For bands flagged as upper limits, residuals are clipped to 0 when
the model is safely below the limit (positive residual), so the
returned vector matches what enters chi_sq band-by-band.
Returns:
| Name | Type | Description |
|---|---|---|
res |
(ndarray, shape(n_filters))
|
|
Source code in ceridwen/observation/photometry.py
ceridwen.observation.Spectrum ¶
Spectrum(
wavelength=None,
flux=None,
uncertainty=None,
mask=slice(None),
noise=None,
name=None,
resolution=None,
calibration=None,
logify_spectrum=False,
smoothtype=None,
inres=0.0,
sky=None,
noise_floor=0.0,
sigma_losvd=None,
fit_sigma_smooth=False,
**kwargs
)
Bases: Observation
Spectroscopic observation.
Stores a densely-sampled spectrum together with optional resolution and multiplicative flux-calibration arrays. Provides helpers for masking spectral regions, computing chi-squared residuals, and projecting the spectrum onto broadband filters.
Parameters:
| Name | Type | Description | Default |
|---|---|---|---|
wavelength
|
(array - like, shape(n_pix))
|
Wavelength grid [Å], vacuum, observed frame (as delivered by the
instrument). |
None
|
flux
|
(array - like, shape(n_pix))
|
Observed flux. Units must be consistent with |
None
|
uncertainty
|
(array - like, shape(n_pix))
|
1-sigma uncertainty, same units as |
None
|
mask
|
array-like of bool or slice, shape (n_pix,)
|
True for pixels that are used (not masked). |
slice(None)
|
resolution
|
float or array - like
|
Instrumental smoothing width. Interpretation depends on
|
None
|
smoothtype
|
(vel, R, 'lambda', lsf)
|
Type of instrumental broadening to apply in |
"vel"
|
inres
|
float
|
Intrinsic resolution of the model library, subtracted in quadrature
before applying the target smoothing. Units match |
0.0
|
calibration
|
(array - like, shape(n_pix))
|
Multiplicative flux-calibration vector (model × calibration ≈ data).
When set, it is applied as a per-pixel multiplicative correction to
the model inside |
None
|
logify_spectrum
|
bool
|
If True, |
False
|
sky
|
(array - like, shape(n_pix))
|
Observed sky background spectrum, same units and pixel grid as
|
None
|
noise_floor
|
float
|
Fractional uncertainty floor applied to the model flux. The effective per-pixel sigma used in chi-squared becomes: .. math:: where :math: |
0.0
|
sigma_losvd
|
float or None
|
Galaxy line-of-sight velocity dispersion [km/s]. When set, an
additional velocity-broadening step is applied to the model
spectrum before any instrumental smoothing specified by
|
None
|
noise
|
GaussianProcess or None
|
If a |
None
|
Examples:
>>> spec = Spectrum(
... wavelength=wave_aa,
... flux=obs_fnu,
... uncertainty=obs_fnu_unc,
... resolution=100.0,
... smoothtype="vel", # 100 km/s instrumental broadening
... noise_floor=0.01, # 1% calibration floor
... sigma_losvd=150.0, # 150 km/s galaxy velocity dispersion
... )
>>> spec.setup_for_model(model_wave)
>>> predicted = spec.predict(model_spectrum, model_wave)
>>> spec.mask_lines([6563., 4861.], dv=500.) # mask Hα, Hβ
>>> chi2 = spec.chi_sq(model_fnu)
>>> coeffs, cal_model = spec.fit_polynomial_calibration(predicted, order=4)
>>> phot = spec.synthetic_photometry(filterset)
Parameters:
| Name | Type | Description | Default |
|---|---|---|---|
smoothtype
|
(vel, R, 'lambda', lsf)
|
Which instrumental smoothing kernel to apply in
|
"vel"
|
inres
|
float
|
Intrinsic (library) resolution of the input model spectrum,
subtracted in quadrature before applying the target smoothing.
Units must match |
0.0
|
Source code in ceridwen/observation/spectrum.py
133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200 201 202 203 204 205 206 207 208 209 210 211 212 213 214 215 216 217 218 219 220 221 222 223 224 | |
setup_for_model ¶
Precompute projection matrices and/or smoothing kernels, then build
_predict_fn — the single callable used by predict.
Must be called once (Python-level, outside JIT) after the model
wavelength grid is known. SedModel.__init__ calls this
automatically.
Two behaviours depending on self.smoothtype:
No smoothing (smoothtype=None)
Builds the dense (n_pix, n_wave) linear-interpolation matrix
_H and sets _predict_fn(spec) = _H @ spec.
With instrumental smoothing (smoothtype in
{"vel", "R", "lambda", "lsf"})
Uses a factory function from sedpy_jax.smoothing to
precompute all FFT grid transforms. The returned closure is
fully JAX-JIT-compilable with respect to the spectrum.
_predict_fn(spec) applies smoothing and interpolation to
the observed pixel grid in one call. _H is also built (used
only by the no-smoothing fast path).
Parameters:
| Name | Type | Description | Default |
|---|---|---|---|
wave_model
|
(array - like, shape(n_wave))
|
Rest-frame model wavelength grid [Å], strictly increasing. |
required |
zred
|
float
|
Fixed redshift at which to precompute the spectral projection.
Default 0 (rest-frame). For |
0.0
|
Source code in ceridwen/observation/spectrum.py
243 244 245 246 247 248 249 250 251 252 253 254 255 256 257 258 259 260 261 262 263 264 265 266 267 268 269 270 271 272 273 274 275 276 277 278 279 280 281 282 283 284 285 286 287 288 289 290 291 292 293 294 295 296 297 298 299 300 301 302 303 304 305 306 307 308 309 310 311 312 313 314 315 316 317 318 319 320 321 322 323 324 325 326 327 328 329 330 331 332 333 334 335 336 337 338 339 340 341 342 343 344 345 346 347 348 349 350 351 352 353 354 355 356 357 358 359 360 361 362 363 364 365 366 367 368 369 370 371 372 373 374 375 376 377 378 379 380 381 382 383 384 385 386 387 388 389 390 391 392 393 394 395 396 397 398 399 400 401 402 403 404 405 406 407 408 409 410 411 412 413 414 415 416 417 418 419 420 421 422 423 424 425 426 427 428 429 430 431 432 433 434 435 436 437 438 439 440 441 442 443 444 445 446 447 448 449 450 451 452 453 454 455 456 457 458 459 460 461 462 463 464 465 466 467 468 469 470 471 472 473 474 475 476 477 478 479 480 481 482 483 484 485 486 487 488 489 490 491 | |
predict ¶
Project the model spectrum onto the observed pixel grid, applying instrumental smoothing if configured.
Calls _predict_fn(spectrum[, sigma_smooth]) which was constructed
by setup_for_model. Depending on self.smoothtype:
None— pure linear interpolation (_H @ spectrum)."vel"/"R"— constant-velocity FFT broadening then interpolation to observed pixels."lambda"— constant-wavelength FFT broadening then interpolation."lsf"— wavelength-dependent LSF broadening (CDF-transform FFT) then interpolation.
In all smoothing cases, the full smooth→interpolate pipeline is a
single closure that is JAX-JIT-compilable with respect to spectrum.
Must call setup_for_model(wave_model) before this method.
Parameters:
| Name | Type | Description | Default |
|---|---|---|---|
spectrum
|
(Array, shape(n_wave))
|
Model spectrum in F_nu units on the model wavelength grid. |
required |
wave_model
|
(Array, shape(n_wave))
|
Model wavelength grid [Å] (accepted for interface consistency;
the grid mapping was precomputed by |
required |
sigma_smooth
|
jax.Array scalar
|
Runtime galaxy LOSVD [km/s] -- the Prospector |
None
|
Returns:
| Type | Description |
|---|---|
(Array, shape(n_pix))
|
Model F_nu (smoothed and) interpolated onto |
Source code in ceridwen/observation/spectrum.py
synthetic_photometry ¶
Project this spectrum onto a FilterSet and return synthetic maggies.
The spectrum is assumed to be in F_nu units (e.g. L_sun/Hz/M_sun) and is converted to F_lambda before the AB-normalised projection.
Parameters:
| Name | Type | Description | Default |
|---|---|---|---|
filterset
|
FilterSet
|
|
required |
Returns:
| Name | Type | Description |
|---|---|---|
maggies |
(ndarray, shape(n_filters))
|
Returns |
Source code in ceridwen/observation/spectrum.py
mask_wavelength_range ¶
Mask pixels with wavelengths in [wave_min, wave_max] Å (inclusive).
Sets self.mask[i] = False for all pixels whose wavelength falls
inside the specified range.
Parameters:
| Name | Type | Description | Default |
|---|---|---|---|
wave_min
|
float
|
Wavelength bounds [Å]. |
required |
wave_max
|
float
|
Wavelength bounds [Å]. |
required |
Source code in ceridwen/observation/spectrum.py
mask_lines ¶
Mask spectral lines by zeroing the mask within ±dv km/s of each line.
Parameters:
| Name | Type | Description | Default |
|---|---|---|---|
line_waves
|
array - like
|
Rest-frame central wavelengths [Å]; redshifted internally by
(1 + |
required |
dv
|
float
|
Half-width to mask on each side [km/s]. Default 1000 km/s. |
1000.0
|
zred
|
float
|
Redshift to apply to |
0.0
|
Source code in ceridwen/observation/spectrum.py
chi_sq ¶
Chi-squared contribution from this spectrum.
Accounts for sky subtraction (self.sky) and a fractional noise
floor (self.noise_floor). If self.logify_spectrum is True,
residuals are computed in log-flux space: Δln f / (σ_eff / f).
Parameters:
| Name | Type | Description | Default |
|---|---|---|---|
model_flux
|
(array - like, shape(n_pix))
|
Model flux on the observed pixel grid (output of |
required |
Returns:
| Name | Type | Description |
|---|---|---|
chi2 |
float
|
Sum of squared normalised residuals over unmasked pixels. |
Source code in ceridwen/observation/spectrum.py
residuals ¶
Per-pixel (sky-corrected data − model) / sigma_eff. Masked pixels are set to NaN.
Returns:
| Name | Type | Description |
|---|---|---|
res |
(ndarray, shape(n_pix))
|
|
Source code in ceridwen/observation/spectrum.py
log_likelihood ¶
Full log-likelihood for this spectrum.
Combines the standard pixel-independent Gaussian log-likelihood with
an optional Gaussian Process (GP) correction for correlated residuals
if self.noise is a GaussianProcess instance.
.. math::
\log\mathcal{L} =
-\tfrac{1}{2}\sum_{i\,\in\,\rm mask} r_i^2
+ \log p_{\rm GP}(\mathbf{r} \mid \mathrm{GP})
where :math:r_i = ((d_i - s_i) - c_i m_i)/\sigma_{\rm eff,i}
(sky :math:s_i and calibration :math:c_i are optional; the
noise floor in :math:\sigma_{\rm eff,i} is scaled against the
calibrated model :math:c_i m_i when calibration is set), and
the GP term is zero if no noise model is set.
Parameters:
| Name | Type | Description | Default |
|---|---|---|---|
model_flux
|
(array - like, shape(n_pix))
|
Model flux on the observed pixel grid. |
required |
Returns:
| Type | Description |
|---|---|
float
|
Log-likelihood (larger is better). |
Source code in ceridwen/observation/spectrum.py
fit_polynomial_calibration ¶
Analytically fit a Chebyshev multiplicative calibration polynomial
P(λ) such that data ≈ P(λ) × model_flux.
The polynomial coefficients are solved at each call via weighted linear least squares, making this suitable for marginalising out the calibration at every likelihood evaluation without a parameter-space penalty.
The polynomial is evaluated in a normalised wavelength coordinate
:math:x \in [-1, 1] using Chebyshev basis functions
:math:T_n(x), which are numerically stable for high orders.
Parameters:
| Name | Type | Description | Default |
|---|---|---|---|
model_flux
|
(array - like, shape(n_pix))
|
Model flux on the observed pixel grid (output of |
required |
order
|
int
|
Polynomial order. 0 = constant, 1 = linear, etc. Default 3. |
3
|
Returns:
| Name | Type | Description |
|---|---|---|
coeffs |
(ndarray, shape(order + 1))
|
Chebyshev polynomial coefficients. |
calibrated_flux |
(ndarray, shape(n_pix))
|
|
Notes
Only unmasked pixels enter the least-squares fit. The returned
calibrated_flux is evaluated over the full pixel grid.
Source code in ceridwen/observation/spectrum.py
ceridwen.observation.Lines ¶
Bases: Observation
Observed nebular emission-line fluxes.
Stores a set of emission-line fluxes together with their FSPS line-array
indices, vacuum rest-frame wavelengths, and per-line 1-sigma uncertainties.
The interface is deliberately compatible with
prospect.observation.Lines: the line_ind attribute holds integer
indices into the FSPS emline_luminosity array, and the alias
mapping exposes "line_inds" as an alias for line_ind so that
existing Prospector model code can address this object without modification.
Beyond Prospector, this class adds JAX-native chi_sq / residuals
(JIT-compilable through a fitter), and mask_by_name / select_by_name
helpers that operate on human-readable line names.
Parameters:
| Name | Type | Description | Default |
|---|---|---|---|
line_ind
|
array-like of int
|
Indices of the observed lines in the FSPS emission-line array
( |
required |
line_names
|
list of str
|
Human-readable names, one per line (e.g. |
None
|
wavelength
|
array-like of float
|
Vacuum rest-frame wavelengths [Å], length = |
None
|
flux
|
array-like of float
|
Observed line fluxes. Units should be consistent with any model
prediction passed to |
required |
uncertainty
|
array-like of float
|
1-sigma line-flux uncertainties, same units as |
required |
mask
|
array-like of bool
|
True for lines to include in chi-squared. Defaults to all-True. |
required |
upper_limit
|
array-like of bool, shape (n_lines,)
|
If True for a given line, that line is treated as a non-detection upper limit rather than a positive detection. The chi-squared contribution for such lines is one-sided: a penalty is applied only when the model flux exceeds the observed value (i.e., the model predicts more emission than the upper limit allows): .. math:: where :math: |
None
|
Examples:
>>> lines = Lines(
... line_ind = [59, 63, 71],
... line_names = ["Hbeta", "[OIII]5007", "Halpha"],
... wavelength = [4861., 5007., 6563.],
... flux = obs_fluxes,
... uncertainty= obs_unc,
... upper_limit= [False, True, False], # [OIII]5007 is a non-detection
... )
>>> lines.mask_by_name(["[OIII]5007"]) # exclude one line
>>> chi2 = lines.chi_sq(model_fluxes)
>>> subset = lines.select_by_name(["Hbeta", "Halpha"])
Source code in ceridwen/observation/lines.py
setup_for_model ¶
Precompute the (n_lines, n_wave) Gaussian-aperture weight matrix
_W that extracts line fluxes from a model spectrum via a single
matrix–vector multiply.
Must be called once before predict (and before JIT-compiling any
function containing predict). SedModel.__init__ calls this
automatically.
Physical description
For each emission line centred at wavelength :math:\lambda_k, the
integrated line flux is estimated as a Gaussian-weighted integral over
the model spectrum:
.. math::
F_k = \int w_k(\lambda)\, f_\nu(\lambda)\, \mathrm{d}\lambda
where
.. math::
w_k(\lambda) = \exp\!\left[-\frac{1}{2}
\left(\frac{\lambda - \lambda_k}{\sigma_k}\right)^2\right],
\quad \sigma_k = \lambda_k\, \frac{\sigma_v}{c}
Discretised with the trapezoidal rule on the model wavelength grid,
this becomes _W @ spectrum where
_W[k, j] = w_k(wave_j) * dlambda_j and dlambda_j are the
trapezoidal quadrature weights.
Parameters:
| Name | Type | Description | Default |
|---|---|---|---|
wave_model
|
(array - like, shape(n_wave))
|
Model wavelength grid [Å], strictly increasing. |
required |
sigma_v
|
float
|
1-sigma Gaussian aperture width [km/s]. Default 200 km/s.
Sufficient to capture narrow nebular lines as generated by
|
200.0
|
Source code in ceridwen/observation/lines.py
147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200 201 202 203 204 205 206 207 208 209 210 211 212 213 214 215 216 217 218 219 220 221 222 223 224 225 226 227 228 229 230 231 232 233 234 235 236 237 238 239 240 241 242 243 244 245 246 247 248 249 250 251 252 | |
predict ¶
Extract emission-line fluxes from the model spectrum via Gaussian-
aperture integration: computes _W @ spectrum where _W was
precomputed once in setup_for_model. On GPU this is a single
GEMV; XLA constant-folds _W into the compiled graph.
Must call setup_for_model(wave_model, sigma_v=...) first.
Parameters:
| Name | Type | Description | Default |
|---|---|---|---|
spectrum
|
(Array, shape(n_wave))
|
Model spectrum in F_nu units. |
required |
wave_model
|
(Array, shape(n_wave))
|
Accepted for interface consistency; not used inside this method. |
required |
Returns:
| Type | Description |
|---|---|
(Array, shape(n_lines))
|
Gaussian-aperture integrated flux for each line. |
Source code in ceridwen/observation/lines.py
chi_sq ¶
Chi-squared contribution from the observed line fluxes.
For lines flagged as upper limits (self.upper_limit[k] = True),
the contribution is one-sided: a penalty is applied only when the
model flux exceeds the observed upper-limit value.
Parameters:
| Name | Type | Description | Default |
|---|---|---|---|
model_fluxes
|
(array - like, shape(n_lines))
|
Predicted line fluxes, same units as |
required |
Returns:
| Name | Type | Description |
|---|---|---|
chi2 |
float
|
|
Source code in ceridwen/observation/lines.py
residuals ¶
Per-line (data − model) / sigma. Masked lines are set to NaN.
Upper-limit lines where the model does not exceed the limit are set
to zero (no tension) rather than showing a negative residual.
Returns:
| Name | Type | Description |
|---|---|---|
res |
(ndarray, shape(n_lines))
|
|
Source code in ceridwen/observation/lines.py
mask_by_name ¶
Exclude lines whose name appears in names from chi-squared.
Sets self.mask[i] = False for all lines whose entry in
self.line_names matches any element of names. A no-op if
self.line_names is not set.
Parameters:
| Name | Type | Description | Default |
|---|---|---|---|
names
|
list of str
|
|
required |
Source code in ceridwen/observation/lines.py
select_by_name ¶
Return a new Lines instance containing only the named lines.
Parameters:
| Name | Type | Description | Default |
|---|---|---|---|
names
|
list of str
|
Must all be present in |
required |
Returns:
| Type | Description |
|---|---|
Lines
|
|
Raises:
| Type | Description |
|---|---|
ValueError
|
If |
KeyError
|
If any element of |
Source code in ceridwen/observation/lines.py
Model¶
ceridwen.model.SedModel ¶
Parameter manager + prediction layer for Ceridwen SED fitting.
Parameters:
| Name | Type | Description | Default |
|---|---|---|---|
csp
|
CSPBasis
|
Initialised composite stellar population model. Must expose
|
required |
observations
|
list of Observation
|
Data containers (Photometry, Spectrum, Lines). Each must have a
unique |
required |
priors
|
dict[str, Prior]
|
Mapping from free-parameter name to a prior object implementing
|
None
|
transforms
|
dict[str, callable]
|
Mapping from derived (CSP) parameter name to a callable that computes its value from the free-parameter dict:: Derived parameters listed here are removed from the free-parameter
list and replaced by the new free parameters supplied via
Example — fitting log-ratios of SFR bins instead of raw SFH:: |
None
|
free_param_init
|
dict[str, Array]
|
Initial values for free parameters that replace derived ones.
Keys in this dict are added to |
None
|
Attributes:
| Name | Type | Description |
|---|---|---|
theta_init |
dict[str, Array]
|
Initial values for the free parameters only (derived params
are absent; their replacements from |
param_names |
list[str]
|
Ordered list of free-parameter names. |
transforms |
dict[str, callable]
|
Registered transforms (empty dict if none). |
obs_dict |
dict[str, Observation]
|
Observations keyed by |
wave |
(Array, shape(n_wave))
|
Model wavelength grid [Å]. |
Source code in ceridwen/model/model.py
141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200 201 202 203 204 205 206 207 208 209 210 211 212 213 214 215 216 217 218 219 220 221 222 223 224 225 226 227 228 229 230 231 232 233 234 235 236 237 238 239 240 241 242 243 244 245 246 247 248 249 250 251 252 253 254 255 256 257 | |
obs_dict
property
¶
Observations as a dict keyed by obs.name.
Pass this to MultiObservationLikelihood.make_lnprobfn as the
observations argument::
lnprobfn = multi_lhood.make_lnprobfn(model.obs_dict, model, model)
apply_transforms ¶
Apply all registered transforms to produce a CSP-compatible model_theta.
Starts from a shallow copy of free_theta and computes each
derived parameter by calling the corresponding transform callable::
model_theta[derived] = transform_fn(free_theta)
The free-parameter keys (e.g. "logsfr_ratios") are kept in
model_theta alongside the derived ones; the CSP simply ignores
any keys it does not recognise.
Parameters:
| Name | Type | Description | Default |
|---|---|---|---|
free_theta
|
dict[str, Array]
|
Free-parameter dict as used by the sampler. |
required |
Returns:
| Name | Type | Description |
|---|---|---|
model_theta |
dict[str, Array]
|
Extended dict suitable for |
Source code in ceridwen/model/model.py
predict ¶
Project the CSP model spectrum onto all observations.
If transforms are registered, theta is treated as the
free-parameter dict; apply_transforms is called first to
obtain the CSP-compatible model_theta before forwarding to
csp.predict.
Internally calls csp.predict(model_theta, observations) which
computes the spectrum once and projects it onto each observation:
Photometry→ synthetic AB maggies via filter convolutionSpectrum→ model F_ν interpolated onto observed wavelength gridLines→ Gaussian-aperture integrated line fluxes
Mass scaling — if "logmass" is present in theta, the
spectrum is multiplied by 10 ** logmass inside csp.predict()
before projection. The logsfr_ratios_to_sfh transform
normalises the SFH so that the trapezoidal integral of SFR over
the lookback grid equals 1 M⊙ (Prospector / FSPS convention), so
this factor sets the physical amplitude for a galaxy with stellar
mass M = 10^logmass M⊙. Scaling once before projection is
more efficient than scaling each observation separately.
Parameters:
| Name | Type | Description | Default |
|---|---|---|---|
theta
|
dict[str, Array]
|
Free-parameter dict (before any transforms). May optionally
include |
required |
Returns:
| Type | Description |
|---|---|
dict[str, Array]
|
Keyed by |
Source code in ceridwen/model/model.py
predict_jit ¶
JIT-compiled version of :meth:predict.
Identical semantics, but the first call triggers XLA compilation
and subsequent calls with the same dict structure hit the compiled
cache. Use this for interactive evaluation (sanity checks,
posterior predictive checks) outside the sampler hot path, where
run_sampler already wraps the full log-posterior in @jax.jit.
For vectorised evaluation over many parameter draws, prefer
:meth:predict_vmap.
Source code in ceridwen/model/model.py
predict_vmap ¶
Vectorised prediction over a batch of parameter dicts.
Parameters:
| Name | Type | Description | Default |
|---|---|---|---|
theta_batch
|
dict[str, Array]
|
Each value has a leading batch dimension, e.g.
|
required |
Returns:
| Type | Description |
|---|---|
dict[str, Array]
|
Each value has a leading batch dimension, e.g.
|
Source code in ceridwen/model/model.py
ln_prior ¶
Evaluate the log-prior for all registered free parameters.
For each parameter p in self.priors, computes
sum(prior.logpdf(theta[p])) (the sum handles vector-valued
parameters such as a non-parametric SFH) and accumulates the total.
Parameters absent from self.priors contribute 0 (flat prior).
Parameters:
| Name | Type | Description | Default |
|---|---|---|---|
theta
|
dict[str, Array]
|
|
required |
Returns:
| Name | Type | Description |
|---|---|---|
lnp |
(Array, scalar)
|
|
Source code in ceridwen/model/model.py
log_prob ¶
Alias for ln_prior.
Required by DiagonalGaussianLikelihood.make_lnprobfn and
MultiObservationLikelihood.make_lnprobfn, which call
prior.log_prob(theta).
Parameters:
| Name | Type | Description | Default |
|---|---|---|---|
theta
|
dict[str, Array]
|
|
required |
Returns:
| Type | Description |
|---|---|
(Array, scalar)
|
|
Source code in ceridwen/model/model.py
summary ¶
Return a multi-line human-readable summary of the model configuration.
Covers: registered free parameters (with shapes and prior types), active transforms (free → derived param mappings), all observation objects, and the CSP physics switch configuration.
Source code in ceridwen/model/model.py
display ¶
Draw a publication-quality probabilistic graphical model (PGM) diagram.
The diagram follows standard PGM conventions:
- Open circles — stochastic latent variables (free parameters θᵢ)
- Stacked circles — vector-valued parameters (e.g. SFH weight vector)
- Double-bordered rectangle — deterministic SED computation f_ν(λ)
- Coloured rectangles — observation projection operators
- Filled circles — observed data (shaded = conditioned upon)
- Dashed arrows — prior ↦ parameter dependency (ε ≡ stochastic edge)
- Solid arrows — deterministic dependency (θ → f_SED → ŷ → y)
The figure adapts dynamically to however many parameters and observations are registered, making it immediately suitable for inclusion in a paper.
Parameters:
| Name | Type | Description | Default |
|---|---|---|---|
ax
|
Axes
|
Axes to draw into. If None a new figure is created. |
None
|
figsize
|
(float, float)
|
Figure size in inches |
None
|
return_fig
|
bool
|
If True return |
False
|
Returns:
| Type | Description |
|---|---|
(fig, ax) or None
|
Only returned when |
Examples:
Source code in ceridwen/model/model.py
524 525 526 527 528 529 530 531 532 533 534 535 536 537 538 539 540 541 542 543 544 545 546 547 548 549 550 551 552 553 554 555 556 557 558 559 560 561 562 563 564 565 566 567 568 569 570 571 572 573 574 575 576 577 578 579 580 581 582 583 584 585 586 587 588 589 590 591 592 593 594 595 596 597 598 599 600 601 602 603 604 605 606 607 608 609 610 611 612 613 614 615 616 617 618 619 620 621 622 623 624 625 626 627 628 629 630 631 632 633 634 635 636 637 638 639 640 641 642 643 644 645 646 647 648 649 650 651 652 653 654 655 656 657 658 659 660 661 662 663 664 665 666 667 668 669 670 671 672 673 674 675 676 677 678 679 680 681 682 683 684 685 686 687 688 689 690 691 692 693 694 695 696 697 698 699 700 701 702 703 704 705 706 707 708 709 710 711 712 713 714 715 716 717 718 719 720 721 722 723 724 725 726 727 728 729 730 731 732 733 734 735 736 737 738 739 740 741 742 743 744 745 746 747 748 749 750 751 752 753 754 755 756 757 758 759 760 761 762 763 764 765 766 767 768 769 770 771 772 773 774 775 776 777 778 779 780 781 782 783 784 785 786 787 788 789 790 791 792 793 794 795 796 797 798 799 800 801 802 803 804 805 806 807 808 809 810 811 812 813 814 815 816 817 818 819 820 821 822 823 824 825 826 827 828 829 830 831 832 833 834 835 836 837 838 839 840 841 842 843 844 845 846 847 848 849 850 851 852 853 854 855 856 857 858 859 860 861 862 863 864 865 866 867 868 869 870 871 872 873 874 875 876 877 878 879 880 881 882 883 884 885 886 887 888 889 890 891 892 893 894 895 896 897 898 899 900 901 902 903 904 905 906 907 908 909 910 911 912 913 914 915 916 917 918 919 920 921 922 923 924 925 926 927 928 929 930 931 932 933 934 935 936 937 938 939 940 941 942 943 944 945 946 947 948 949 950 951 952 953 954 955 956 957 958 959 960 961 962 963 964 965 966 967 968 969 970 971 972 973 974 975 976 977 978 979 980 981 982 983 984 985 986 987 988 989 990 991 992 993 994 995 996 997 998 999 1000 1001 1002 1003 1004 1005 1006 1007 1008 1009 1010 1011 1012 1013 1014 1015 1016 1017 1018 1019 1020 1021 1022 1023 1024 1025 1026 1027 1028 1029 1030 1031 1032 1033 1034 1035 1036 1037 1038 1039 1040 1041 | |
ceridwen.model.transforms.logsfr_ratios_to_sfh ¶
Convert log-ratios of consecutive SFR bins to a unit-mass SFH weight vector.
Parameters:
| Name | Type | Description | Default |
|---|---|---|---|
logsfr_ratios
|
(array_like, shape(n - 1))
|
Log10 ratios of consecutive SFR bins.
Under the lookback-time convention (index 0 =
today, last index = oldest), |
required |
sfh_times_yr
|
(array_like, shape(n))
|
Lookback-time grid in years (same as |
None
|
Returns:
| Name | Type | Description |
|---|---|---|
sfh |
(ndarray, shape(n))
|
Unit-mass SFH weight vector suitable for |
Notes
The normalisation enforces sum(sfh * w) = 1 where w are
the standard trapezoidal quadrature weights (half-width at
boundaries), so that the trapezoidal integral of the SFH over the
lookback-time grid equals 1 Msun:
.. math:: \int_0^{t_{\rm univ}} \mathrm{SFR}(t)\,dt = 1\;\mathrm{M_\odot}.
The total stellar mass of the model is then set by
theta["logmass"] (applied as a multiplicative 10**logmass
inside CSPBasis.predict), matching the Prospector / FSPS
convention.
This is consistent with the per-bin mass integral
m2 = sfh_mid * dt computed inside
CSPBasis.calculate_ssp_weights_const_zh_step and the
piecewise-linear integral used by
calculate_ssp_weights_const_zh: on a shared node grid the
trapezoid sum here and the midpoint sum there are algebraically
identical, so no compensating rescale is needed inside the CSP.
JAX compatibility
The function is fully JIT-compatible: every operation is a
jnp primitive on traced arrays, the sfh_times_yr is not
None branch resolves at trace time (it is a Python-level check
on the closure argument, not a runtime decision on a traced
value), and there are no data-dependent shapes.
Source code in ceridwen/model/transforms.py
79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 | |
Priors¶
ceridwen.priors ¶
Prior distributions for CERIDWEN models.
A clean, discoverable import path for the priors (implemented in
ceridwen.sampler.priors)::
from ceridwen.priors import Uniform, Normal, ClippedNormal, LogNormal, StudentT
Each prior exposes log_prob, sample, unit_transform (inverse-CDF,
used by nested sampling) and inverse_unit_transform (CDF), and is a JAX
PyTree so it can flow through jit/grad.
Prior
dataclass
¶
Bases: ABC
JAX-friendly prior base class that delegates all probability operations to
a TFP-JAX distribution. Subclasses must implement tfp_dist().
Source code in ceridwen/sampler/priors.py
tfp_dist
abstractmethod
¶
Return a TFP-JAX distribution object built from self.params.
Must be implemented by subclasses.
Uniform
dataclass
¶
Bases: Prior
Uniform distribution on [low, high].
Source code in ceridwen/sampler/priors.py
TopHat
dataclass
¶
Bases: Uniform
Uniform distribution between two bounds, renamed for backwards compatibility :param low: Minimum of the distribution
:param high: Maximum of the distribution
Source code in ceridwen/sampler/priors.py
Normal
dataclass
¶
Bases: Prior
A simple gaussian prior.
:param mean: Mean of the distribution
:param sigma: Standard deviation of the distribution
Source code in ceridwen/sampler/priors.py
ClippedNormal
dataclass
¶
Bases: Prior
A Gaussian prior clipped to some range.
:param mean: Mean of the normal distribution
:param sigma: Standard deviation of the normal distribution
:param low: Minimum of the distribution
:param high: Maximum of the distribution
Source code in ceridwen/sampler/priors.py
LogNormal
dataclass
¶
Bases: Prior
A log-normal prior, where the natural log of the variable is distributed normally. Useful for parameters that cannot be less than zero.
Note that LogNormal(np.exp(mode) / f) == LogNormal(np.exp(mode) * f)
and f = np.exp(sigma) corresponds to "one sigma" from the peak.
:param mode: Natural log of the variable value at which the probability density is highest.
:param sigma: Standard deviation of the distribution of the natural log of the variable.
Source code in ceridwen/sampler/priors.py
StudentT
dataclass
¶
Bases: Prior
A Student's T distribution
:param mean: Mean of the distribution
:param scale: Size of the distribution, analogous to the standard deviation
:param df: Number of degrees of freedom
Source code in ceridwen/sampler/priors.py
Likelihood¶
ceridwen.likelihood.DiagonalGaussianLikelihood
dataclass
¶
Bases: LikelihoodBase
Gaussian log-likelihood with an independent (diagonal) noise model.
This class sequences three operations:
- Call
noise_model.compute(sigma_obs, mu, mask, params)to obtain per-datuminv_varandlog_det. - Call
lnlike_diag_gaussian(y, mu, inv_var, log_det, mask)to get the scalar log-likelihood and diagnostics.
Because both steps are pure JAX functions, the entire __call__ is
JIT-compilable and differentiable end-to-end.
Parameters:
| Name | Type | Description | Default |
|---|---|---|---|
noise_model
|
DiagonalNoiseModel
|
Noise model instance. Defaults to a plain observational-uncertainty- only model (no jitter, no fractional error). |
DiagonalNoiseModel()
|
Examples:
Basic usage with no nuisance parameters::
lhood = DiagonalGaussianLikelihood()
lnl, aux = lhood(y, mu, sigma_obs, mask)
# aux.chi -- normalised residuals
# aux.lnl_pointwise -- per-datum contributions
With jitter as a sampled parameter::
lhood = DiagonalGaussianLikelihood(
noise_model=DiagonalNoiseModel(use_jitter=True)
)
lnl, aux = lhood(y, mu, sigma_obs, mask,
params={"log_jitter": theta["log_jitter"]})
Gradient of log-likelihood w.r.t. theta (for HMC/gradient-based samplers)::
def lnl_fn(theta):
mu = model.predict(theta)
return lhood(y, mu, sigma_obs, mask, params=theta)[0]
grad = jax.grad(lnl_fn)(theta)
Using has_aux=True to get diagnostics and gradient in one pass::
(lnl, aux), grad = jax.value_and_grad(
lambda t: lhood(y, model.predict(t), sigma_obs, mask, t),
has_aux=True,
)(theta)
__call__ ¶
Evaluate log-likelihood and return diagnostics.
This method is the intended hot path for sampling. It is safe inside
jax.jit, jax.grad, and jax.vmap.
Parameters:
| Name | Type | Description | Default |
|---|---|---|---|
y
|
(Array, shape(n_data))
|
|
required |
mu
|
(Array, shape(n_data))
|
|
required |
sigma_obs
|
(Array, shape(n_data))
|
|
required |
mask
|
Array of bool, shape (n_data,)
|
|
required |
params
|
dict
|
Nuisance parameters expected by the noise model. |
None
|
Returns:
| Name | Type | Description |
|---|---|---|
lnl_total |
(Array, scalar)
|
|
aux |
LikelihoodOutput
|
|
Source code in ceridwen/likelihood/likelihood.py
make_lnprobfn ¶
Build a JIT-compiled log-posterior for a single observation type.
observations must have attributes .flux, .uncertainty, and
.mask (the interface defined by ceridwen.observation.Observation).
The returned function has the exact signature expected by blackjax::
kernel = blackjax.nuts(lnprobfn, step_size)
state = kernel.init(initial_theta)
Parameters:
| Name | Type | Description | Default |
|---|---|---|---|
observations
|
Observation
|
A single observation object (photometry or spectroscopy or
lines). For multiple modalities use |
required |
model
|
object
|
Must implement |
required |
prior
|
object
|
Must implement |
required |
Returns:
| Type | Description |
|---|---|
Callable[[theta], Array]
|
JIT-compiled log-posterior. |
Notes
observations, model, prior, and self are all closed
over at factory time and become compile-time constants. Only theta
is traced. This means the first call incurs a one-time XLA compilation
cost; subsequent calls are fully compiled and have minimal Python
overhead.
Source code in ceridwen/likelihood/likelihood.py
ceridwen.likelihood.MultiObservationLikelihood
dataclass
¶
Bases: LikelihoodBase
Combine log-likelihoods across multiple independent data modalities.
In SED fitting, the posterior typically conditions on several observation types simultaneously -- broadband photometry (in maggies), a spectrum (in F_lambda), and emission-line fluxes. Each has different units, different noise properties, and potentially different noise model nuisance parameters.
This class owns a static mapping from string keys to likelihood objects. At trace time the Python loop over keys is unrolled by XLA; there is no runtime dispatch overhead.
Parameters:
| Name | Type | Description | Default |
|---|---|---|---|
keys
|
tuple of str
|
Ordered observation keys, e.g. |
tuple()
|
likelihoods
|
tuple of LikelihoodBase
|
One likelihood per key, in the same order. Different elements may have different noise models (e.g. photometry with fractional error, spectroscopy with jitter). |
tuple()
|
Examples:
Combine photometry and spectroscopy::
multi = MultiObservationLikelihood(
keys=("phot", "spec"),
likelihoods=(
DiagonalGaussianLikelihood(
DiagonalNoiseModel(use_fractional=True)
),
DiagonalGaussianLikelihood(
DiagonalNoiseModel(use_jitter=True)
),
),
)
lnprobfn = multi.make_lnprobfn(observations, model, prior)
where observations is a dict {"phot": phot_obs, "spec": spec_obs}
and model.predict(theta) returns a dict with the same keys.
__call__ ¶
Evaluate total log-likelihood across all modalities.
Parameters:
| Name | Type | Description | Default |
|---|---|---|---|
y
|
dict[str, Array]
|
Per-modality data arrays, keyed by the same strings as |
required |
mu
|
dict[str, Array]
|
Per-modality data arrays, keyed by the same strings as |
required |
sigma_obs
|
dict[str, Array]
|
Per-modality data arrays, keyed by the same strings as |
required |
mask
|
dict[str, Array]
|
Per-modality data arrays, keyed by the same strings as |
required |
params
|
dict[str, Array]
|
All nuisance parameters (shared across modalities). Each likelihood/noise model extracts only the keys it needs. |
None
|
Returns:
| Name | Type | Description |
|---|---|---|
lnl_total |
(Array, scalar)
|
Sum of log-likelihoods across all modalities. |
aux |
dict[str, LikelihoodOutput]
|
Per-modality diagnostics. |
Source code in ceridwen/likelihood/likelihood.py
make_lnprobfn ¶
Build a JIT-compiled log-posterior for multiple observation types.
Parameters:
| Name | Type | Description | Default |
|---|---|---|---|
observations
|
dict[str, Observation]
|
Keyed by the same strings as |
required |
model
|
object
|
Must implement |
required |
prior
|
object
|
Must implement |
required |
Returns:
| Type | Description |
|---|---|
Callable[[theta], Array]
|
JIT-compiled log-posterior, suitable for blackjax. |
Source code in ceridwen/likelihood/likelihood.py
Fitting¶
ceridwen.fit.fitSED ¶
fitSED(
model,
observations=None,
output_dir=".",
*,
sampler="nested",
rng_key=None,
sampler_kwargs=None,
vi=None,
vi_kwargs=None,
filename="ceridwen_result.h5",
overwrite=True,
verbose=True
)
Fit an SED model to observations and write results to HDF5.
This is the main user-facing function. It constructs the likelihood, configures and runs the requested sampler, and writes the posterior samples, model metadata, and observation data to disk.
Parameters:
| Name | Type | Description | Default |
|---|---|---|---|
model
|
SedModel
|
A fully configured |
required |
observations
|
list[Observation]
|
Observations to fit. If None, uses |
None
|
output_dir
|
str or Path
|
Directory for the output HDF5 file. Created if it does not exist. |
'.'
|
sampler
|
str
|
|
'nested'
|
rng_key
|
Array
|
JAX PRNG key. Defaults to |
None
|
sampler_kwargs
|
dict
|
Extra keyword arguments forwarded to the sampler adapter
constructor (e.g. |
None
|
vi
|
None, str, VariationalMap, or TrainedMap
|
Variational-inference preconditioning (only used when
When a map is supplied, NUTS samples in the whitened z-space, giving dramatically shorter warmup. See Hoffman et al. 2019, arXiv:1903.03704. |
None
|
vi_kwargs
|
dict
|
Forwarded to :func: |
None
|
filename
|
str
|
Name of the output HDF5 file. Default |
'ceridwen_result.h5'
|
overwrite
|
bool
|
If True (default), overwrite an existing file. |
True
|
verbose
|
bool
|
Print progress to the console. Default True. |
True
|
Returns:
| Type | Description |
|---|---|
SamplingResult
|
The sampling result object (same as returned by |
Source code in ceridwen/fit.py
67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200 201 202 203 204 205 206 207 208 209 210 211 212 213 214 215 216 217 218 219 220 221 222 223 224 225 226 227 228 229 230 231 232 233 234 235 236 237 238 239 240 241 242 243 244 245 246 247 248 249 250 251 252 253 254 255 256 257 258 259 260 261 262 263 264 265 266 267 268 269 270 271 272 273 274 275 276 277 278 279 280 281 282 | |
ceridwen.sampler.run_sampler ¶
Build JAX callables from a Ceridwen model and dispatch to a sampler.
This is the primary user-facing entry point for posterior sampling in Ceridwen. It is intentionally thin:
- Constructs a JIT-compiled
loglike_fnthat sums log-likelihoods over all registered observations (no prior contribution). - Wraps
model.ln_prioras a JIT-compiledlogprior_fn. - Delegates to
adapter.run(loglike_fn, logprior_fn, theta_init, key).
Parameters:
| Name | Type | Description | Default |
|---|---|---|---|
model
|
SedModel
|
Initialised model with observations set up and priors registered.
|
required |
likelihood
|
MultiObservationLikelihood
|
Likelihood object for the registered observations. Its |
required |
adapter
|
SamplerAdapter
|
Concrete backend (e.g. |
required |
rng_key
|
Array
|
JAX PRNGKey for reproducibility. |
required |
Returns:
| Type | Description |
|---|---|
SamplingResult
|
Call |
Examples:
Nested sampling with BlackJAX NSS::
from ceridwen.sampler import run_sampler
from ceridwen.sampler.nested import BlackJAXNestedSamplerAdapter
adapter = BlackJAXNestedSamplerAdapter(
priors = model.priors,
num_live = 500,
num_inner_steps = len(model.param_names) * 5,
)
result = run_sampler(model, multi_likelihood, adapter,
jax.random.PRNGKey(42))
ns = result.to_anesthetic(labels={"Z": r"$\log Z/Z_\odot$"})
print(result.summary())
Future HMC adapter (same call signature)::
adapter = BlackJAXHMCAdapter(step_size=0.01, n_warmup=500)
result = run_sampler(model, multi_likelihood, adapter, key)
Source code in ceridwen/sampler/runner.py
275 276 277 278 279 280 281 282 283 284 285 286 287 288 289 290 291 292 293 294 295 296 297 298 299 300 301 302 303 304 305 306 307 308 309 310 311 312 313 314 315 316 317 318 319 320 321 322 323 324 325 326 327 328 329 330 331 332 333 334 335 336 337 338 339 340 341 342 343 344 345 346 347 348 349 350 351 352 353 354 355 356 357 358 359 360 361 362 363 364 365 366 | |
ceridwen.sampler.nested.BlackJAXNestedSamplerAdapter ¶
BlackJAXNestedSamplerAdapter(
priors,
num_live=500,
num_inner_steps=None,
num_delete=None,
logZ_tol=-3.0,
verbose=True,
checkpoint_interval_s=1200.0,
checkpoint_dir=None,
)
Bases: SamplerAdapter
Adapter wrapping blackjax.nss for use with any Ceridwen SedModel.
Handles three concerns that are specific to nested sampling in Ceridwen:
-
Live-point initialisation — samples each free parameter independently from its registered prior using the
Prior.samplemethod fromceridwen.sampler.priors. -
Shape reconciliation — the BlackJAX NSS live-point dict has shape
{name: (num_live, *param_shape)}. The step function vmaps over axis 0, delivering single-particle slices of shape(*param_shape,)tologlike_fn/logprior_fn. This matches the Ceridwen dict-theta convention exactly. -
Evidence extraction — optionally uses
anestheticfor a more accurate :math:\ln Zestimate with uncertainty.
Parameters:
| Name | Type | Description | Default |
|---|---|---|---|
priors
|
dict[str, Prior]
|
Mapping from free-parameter name to a Ceridwen |
required |
num_live
|
int
|
Number of live points. Default 500. |
500
|
num_inner_steps
|
int
|
Inner MCMC steps per NS iteration. Default |
None
|
num_delete
|
int
|
Live points discarded per iteration. Default |
None
|
logZ_tol
|
float
|
Convergence threshold on :math: |
-3.0
|
verbose
|
bool
|
Print a |
True
|
Source code in ceridwen/sampler/nested.py
load_checkpoint
staticmethod
¶
Load a checkpoint / rescue pickle written by this adapter.
Returns the dict {positions, loglikelihood, loglikelihood_birth,
logZ, n_dead, partial}. A partial=True snapshot is a usable
(under-converged) posterior from a run killed before convergence ---
feed positions + loglikelihood + loglikelihood_birth to
anesthetic.NestedSamples exactly as the end-of-run path does.
Source code in ceridwen/sampler/nested.py
run ¶
Run BlackJAX NSS and return a SamplingResult.
Parameters:
| Name | Type | Description | Default |
|---|---|---|---|
loglike_fn
|
callable
|
JIT-compiled log-likelihood (no prior). |
required |
logprior_fn
|
callable
|
JIT-compiled log-prior (must be proper). |
required |
theta_init
|
dict[str, Array]
|
Reference parameter dict (shapes / dtypes). |
required |
rng_key
|
Array
|
|
required |
Returns:
| Type | Description |
|---|---|
SamplingResult
|
|
Source code in ceridwen/sampler/nested.py
290 291 292 293 294 295 296 297 298 299 300 301 302 303 304 305 306 307 308 309 310 311 312 313 314 315 316 317 318 319 320 321 322 323 324 325 326 327 328 329 330 331 332 333 334 335 336 337 338 339 340 341 342 343 344 345 346 347 348 349 350 351 352 353 354 355 356 357 358 359 360 361 362 363 364 365 366 367 368 369 370 371 372 373 374 375 376 377 378 379 380 381 382 383 384 385 386 387 388 389 390 391 392 393 394 395 396 397 398 399 400 401 402 403 404 405 406 407 408 409 410 411 412 413 414 415 416 417 418 419 420 421 422 423 424 425 426 427 428 429 430 431 432 433 434 435 436 437 438 439 440 441 442 443 444 445 446 447 448 449 450 451 452 453 454 455 456 457 458 459 460 461 462 463 464 465 466 467 468 469 470 471 472 473 474 475 476 477 478 479 480 481 482 483 484 485 486 487 488 489 490 491 492 493 494 495 496 497 498 499 500 501 502 503 504 505 506 507 508 509 510 511 512 513 514 515 516 517 518 519 520 521 522 523 524 525 526 527 528 529 530 531 532 533 534 535 536 537 538 539 540 541 542 543 544 545 546 547 548 549 550 551 552 553 554 555 556 557 558 559 560 561 562 563 564 565 566 567 568 569 570 571 572 573 574 575 576 577 578 579 580 581 | |
ceridwen.sampler.nuts.BlackJAXNUTSAdapter ¶
BlackJAXNUTSAdapter(
num_warmup=None,
num_samples=2000,
num_chains=4,
initial_step_size=None,
target_acceptance=0.95,
max_num_doublings=10,
dense_mass=None,
bounds=None,
vi=None,
vi_kwargs=None,
verbose=True,
)
Bases: SamplerAdapter
Adapter wrapping blackjax.nuts with window adaptation for Ceridwen.
Bounded (uniform-prior) parameters are automatically reparameterised onto an unconstrained space via sigmoid/logit, eliminating the hard boundary walls that cause divergent transitions.
Optionally, a variational transport map (see :mod:ceridwen.sampler.vi)
may be supplied via vi. When set, the adapter trains the map
against the unconstrained posterior and then runs NUTS on the
whitened target :math:\log p(f(z)) + \log|\partial f/\partial z|
(Hoffman et al. 2019, arXiv:1903.03704). In whitened space the
target is approximately :math:\mathcal{N}(0, I) so identity mass
matrix and step size :math:\mathcal{O}(1) are near-optimal,
giving dramatically shorter warmup.
Parameters:
| Name | Type | Description | Default |
|---|---|---|---|
num_warmup
|
int
|
Number of warmup (adaptation) steps per chain. Default 1500
for native NUTS. When |
None
|
num_samples
|
int
|
Number of post-warmup posterior draws per chain. Default 2000. |
2000
|
num_chains
|
int
|
Number of independent chains. Default 4. |
4
|
initial_step_size
|
float
|
Starting step size for the leapfrog integrator before adaptation.
Default 0.01 for native NUTS. When |
None
|
target_acceptance
|
float
|
Target acceptance probability for dual averaging. Default 0.95; higher values reduce divergences in VI-preconditioned NUTS. |
0.95
|
max_num_doublings
|
int
|
Maximum tree depth (2^max_num_doublings leapfrog steps). Default 10. |
10
|
dense_mass
|
bool
|
Use a dense (full) inverse mass matrix. Default True for
native NUTS. When |
None
|
bounds
|
dict
|
Maps parameter names to (low, high) tuples for bounded params.
If None (default), bounds are auto-detected from the model priors
passed through |
None
|
vi
|
None, str, VariationalMap, or TrainedMap
|
Variational preconditioning mode.
|
None
|
vi_kwargs
|
dict
|
Forwarded either to the VI map constructor (when |
None
|
verbose
|
bool
|
Print progress information. Default True. |
True
|
Source code in ceridwen/sampler/nuts.py
run ¶
Run BlackJAX NUTS with window adaptation.
Strategy (GPU-optimised): 1. Run ONE warmup to adapt step size + mass matrix. 2. Build a single NUTS kernel from the adapted parameters. 3. vmap the sampling across all chains in parallel. This compiles ONE XLA program and runs all chains simultaneously, fully utilising GPU parallelism.
Parameters:
| Name | Type | Description | Default |
|---|---|---|---|
loglike_fn
|
callable
|
JIT-compiled log-likelihood (no prior). |
required |
logprior_fn
|
callable
|
JIT-compiled log-prior. |
required |
theta_init
|
dict[str, Array]
|
Initial parameter values with correct shapes. |
required |
rng_key
|
Array
|
|
required |
Returns:
| Type | Description |
|---|---|
SamplingResult
|
|
Source code in ceridwen/sampler/nuts.py
330 331 332 333 334 335 336 337 338 339 340 341 342 343 344 345 346 347 348 349 350 351 352 353 354 355 356 357 358 359 360 361 362 363 364 365 366 367 368 369 370 371 372 373 374 375 376 377 378 379 380 381 382 383 384 385 386 387 388 389 390 391 392 393 394 395 396 397 398 399 400 401 402 403 404 405 406 407 408 409 410 411 412 413 414 415 416 417 418 419 420 421 422 423 424 425 426 427 428 429 430 431 432 433 434 435 436 437 438 439 440 441 442 443 444 445 446 447 448 449 450 451 452 453 454 455 456 457 458 459 460 461 462 463 464 465 466 467 468 469 470 471 472 473 474 475 476 477 478 479 480 481 482 483 484 485 486 487 488 489 490 491 492 493 494 495 496 497 498 499 500 501 502 503 504 505 506 507 508 509 510 511 512 513 514 515 516 517 518 519 520 521 522 523 524 525 526 527 528 529 530 531 532 533 534 535 536 537 538 539 540 541 542 543 544 545 546 547 548 549 550 551 552 553 554 555 556 557 558 559 560 561 562 563 564 565 566 567 568 569 570 571 572 573 574 575 576 577 578 579 580 581 582 583 584 585 586 587 588 589 590 591 592 593 594 595 596 597 598 599 600 601 602 603 604 605 606 607 608 609 610 611 612 613 614 615 616 617 618 619 620 621 622 623 624 625 626 627 628 629 630 631 632 633 634 635 636 637 638 639 640 641 642 643 644 645 646 647 648 649 650 651 652 653 654 655 656 657 658 659 660 661 662 663 664 665 666 667 668 669 670 671 672 673 674 675 676 677 678 679 680 681 682 683 684 685 686 687 688 689 690 691 692 693 694 695 696 697 698 699 700 701 702 703 704 705 706 707 708 709 710 711 712 713 714 715 716 717 718 719 720 721 722 723 724 725 726 727 728 729 730 731 | |
Dust, nebular, IGM, cosmology¶
ceridwen.dust.DustModel ¶
Dust ¶
JAX-compatible modular dust model that supports multiple attenuation laws per bin.
Parameters are passed as plain dicts (dict[str, Array]).
Call Dust.describe_attenuation_laws() to list all available models.
Source code in ceridwen/dust/DustModel.py
compute_attenuation ¶
Compute bin-wise attenuation curves.
Parameters:
| Name | Type | Description | Default |
|---|---|---|---|
wave
|
ndarray
|
Wavelength array in Angstroms. |
required |
fit_params
|
dict[str, Array]
|
Parameter dict. Each law wrapper extracts only the keys it needs. |
required |
Returns:
| Type | Description |
|---|---|
(ndarray, shape(num_bins, len(wave)))
|
|
Source code in ceridwen/dust/DustModel.py
get_default_fit_params ¶
Return a plain dict of default fit parameters.
Keys are the parameter names used in the active dust laws; values are JAX scalars. Returned as a dict so it merges directly into the global theta dict.
Source code in ceridwen/dust/DustModel.py
DiffuseDust ¶
Bases: Dust
Single-bin dust model (one law covering all ages) with diffuse_ prefixed
parameter names to avoid collisions with birth-cloud parameters in a shared
theta dict.
Source code in ceridwen/dust/DustModel.py
get_default_params ¶
Return a plain dict of default diffuse-dust parameters (diffuse_* keys).
Source code in ceridwen/dust/DustModel.py
compute_attenuation ¶
Compute the diffuse attenuation curve (single bin).
Parameters:
| Name | Type | Description | Default |
|---|---|---|---|
wave
|
ndarray
|
|
required |
fit_params
|
dict[str, Array]
|
|
required |
Returns:
| Type | Description |
|---|---|
(ndarray, shape(len(wave)))
|
|
Source code in ceridwen/dust/DustModel.py
make_law_wrapper ¶
Return a JAX-traceable wrapper that extracts named params from a dict.
Parameters missing from fit_params fall back to defaults (the
registry defaults), so a law parameter can be pinned via theta without
breaking older suites whose theta dicts never carried it. The
membership test is a static Python check resolved at trace time.
Source code in ceridwen/dust/DustModel.py
modify_function ¶
Return a parameter-renamed wrapper of func (no exec).
When the same attenuation law is used for multiple age bins, its parameters
are suffixed with the bin number so they do not collide in the shared
theta dict (e.g. tau_pow → tau_pow{number}). The returned wrapper
is invoked positionally by :func:make_law_wrapper
(func(wave, *args)); the renaming is purely so that
inspect.signature reports the suffixed names (which :class:Dust reads
to build the per-bin parameter-extraction list). We therefore forward the
positional arguments unchanged and attach a renamed __signature__.
The single-law path never calls this function and is unaffected.
Source code in ceridwen/dust/DustModel.py
ceridwen.dust.DustEmission ¶
DustEmission ¶
Initialize the DustEmission object with parameters for dust emission modeling.
Parameters:
| Name | Type | Description | Default |
|---|---|---|---|
duste_model
|
str
|
Dust emission model to use: 'DL07' or 'THEMIS'. |
'DL07'
|
dust_file
|
str
|
Path to the dust emission file (required). |
None
|
spec_lambda
|
ndarray
|
Wavelength grid over which dust emission will be evaluated (required). |
None
|
kwargs
|
dict
|
Optional keyword arguments to override default dust parameters. Supported: duste_qpah, duste_umin, duste_gamma |
{}
|
Source code in ceridwen/dust/DustEmission.py
__repr__ ¶
Custom string representation of the DustEmission object. Provides a readable summary of model settings and parameters.
Source code in ceridwen/dust/DustEmission.py
get_default_params ¶
Return a plain dict of default dust-emission fit parameters.
Returned as a dict so it merges directly into the global theta dict.
Source code in ceridwen/dust/DustEmission.py
update_dust_params ¶
Update dust parameters for emission calculations.
Parameters: duste_qpah (float, optional): New PAH fraction. duste_umin (float, optional): New minimum U radiation field. duste_gamma (float, optional): New fraction of high U component.
Source code in ceridwen/dust/DustEmission.py
compute_dust_emission ¶
compute_dust_emission(
spec_attn,
spec_dustfree,
spec_lambda,
diffuse_curve,
duste_qpah,
duste_umin,
duste_gamma,
)
Compute dust emission using precomputed trapezoidal weights for fast integration (bolometric luminosities via jnp.dot against the static weight vector).
Parameters: spec_attn (jnp.ndarray): Attenuated spectrum after dust absorption. spec_dustfree (jnp.ndarray): Stellar spectrum before attenuation. spec_lambda (jnp.ndarray): Wavelength array in Angstroms. diffuse_curve (jnp.ndarray): exp(-tau_diffuse), shape (1, n_wave) or (n_wave,). duste_qpah (float): PAH fraction. duste_umin (float): Minimum U radiation field. duste_gamma (float): Fraction of high U component.
Returns: tuple: (spectrum with dust emission, dust mass, total dust emission)
Source code in ceridwen/dust/DustEmission.py
208 209 210 211 212 213 214 215 216 217 218 219 220 221 222 223 224 225 226 227 228 229 230 231 232 233 234 235 236 237 238 239 240 241 242 243 244 245 246 247 248 249 250 251 252 253 254 255 256 257 258 259 260 261 262 263 264 265 266 267 268 269 270 271 272 273 274 275 276 277 278 279 280 281 282 283 284 285 286 287 | |
ceridwen.neb.NebularModel ¶
NebularModel(
cloudy_dust,
sps_home,
csp_lambda,
ssp_flux=None,
ssp_ages_lgyr=None,
isoc_type="mist",
nebnz=11,
nebnage=10,
nebnip=7,
smooth_velocity=True,
sigma_smooth=100.0,
res_floor_factor=2.0,
nebular_smooth_init=None,
)
CLOUDY-grid nebular emission model — physically strict variant.
The continuum cube and the line cube are each interpolated against
their own (logZ, age, logU) axes, ensuring that the returned
nebular spectrum corresponds to the parameters CLOUDY was actually
run at. This differs from FSPS's run-time convention; see the
module docstring and
:class:ceridwen.neb.NebularGridModel_fsps_match.NebularModelFSPSMatch
for the FSPS-matching variant.
Parameters:
| Name | Type | Description | Default |
|---|---|---|---|
cloudy_dust
|
bool
|
|
required |
sps_home
|
str | Path
|
FSPS root directory; expects |
required |
csp_lambda
|
(array, shape(nspec))
|
Wavelength grid (Å) the nebular spectrum will be projected onto. |
required |
ssp_flux
|
(array, shape(n_z, n_age, n_wave) or None)
|
Optional. If provided, |
None
|
ssp_ages_lgyr
|
(array, shape(n_age) or None)
|
Optional. |
None
|
isoc_type
|
(mist, pdva, prsc, bpss)
|
Isochrone tag identifying the ZAU file suffix. |
'mist'
|
nebnz
|
int
|
Dimensions of each CLOUDY grid (default 11, 10, 7 for MIST). |
11
|
nebnage
|
int
|
Dimensions of each CLOUDY grid (default 11, 10, 7 for MIST). |
11
|
nebnip
|
int
|
Dimensions of each CLOUDY grid (default 11, 10, 7 for MIST). |
11
|
smooth_velocity
|
bool
|
|
True
|
sigma_smooth
|
float
|
Line-broadening σ. Default 0.0 (matches FSPS). |
100.0
|
nebular_smooth_init
|
float or None
|
Backwards-compatible alias for |
None
|
Attributes:
| Name | Type | Description |
|---|---|---|
nebem_cont |
(nspec, nebnz, nebnage, nebnip)
|
|
nebem_line |
(nemline, nebnz, nebnage, nebnip)
|
|
nebem_cont_logz, nebem_cont_age, nebem_cont_logu |
1 - D
|
The |
nebem_line_logz, nebem_line_age, nebem_line_logu |
1 - D
|
The |
nebem_logz, nebem_age, nebem_logu |
1 - D
|
Legacy aliases — point at the line axes for compatibility
with code (e.g. |
nebem_line_pos |
(nemline,)
|
Rest-frame emission-line centroids in Å. |
gaussnebarr |
(nspec, nemline)
|
Pre-computed Gaussian profiles including the FSPS |
log_qq |
(n_z, n_age) or None
|
Self-consistent |
young_mask, young_idx |
Boolean mask / index array over |
Source code in ceridwen/neb/NebularGridModel.py
compute_log_qq ¶
log10(Q) for every (Z, age) SSP, matching FSPS's run-time
formula::
qq = ∫ L_nu / lambda dλ (0 < λ < 912 Å)
Q = (L_sun_erg / h_erg_s) × qq (photons / s)
Forced to float64 because the deep-UV fluxes are tiny.
Source code in ceridwen/neb/NebularGridModel.py
evaluate ¶
Single-point evaluation — returns (cont, lines) in Lsun/Hz.
The continuum cube is interpolated against the cont axes, the line cube against the line axes (each cube is therefore evaluated at the physical point CLOUDY was actually run at).
Source code in ceridwen/neb/NebularGridModel.py
evaluate_batch ¶
Vectorised evaluation for all (Z_ssp, age_young) pairs at a
single (logZ_gas, logU). Each cube is bilinearly collapsed
in (Z, U) on its OWN axes, then linearly interpolated in age
across the young-SSP set.
Parameters:
| Name | Type | Description | Default |
|---|---|---|---|
return_components
|
bool
|
When |
False
|
Source code in ceridwen/neb/NebularGridModel.py
ceridwen.igm ¶
ceridwen/igm.py¶
Intergalactic-medium (IGM) absorption models for redshift-aware SED fitting. Applied as a wavelength-dependent transmission curve to the rest-frame source spectrum as a function of the source redshift; the returned factor can simply be multiplied into the spectrum before the observation projection.
Interface
:class:IGMModel is the abstract base. Each concrete subclass
implements tau(lam_rest, zred) returning the optical depth on the
rest-frame wavelength grid at the given source redshift. The public
:meth:IGMModel.attenuation method returns exp(-tau * factor) —
the transmission curve — with an optional scalar factor matching
FSPS's igm_factor fudge (so the IGM strength can be fit).
Models
:class:Madau1995
Line-by-line port of FSPS's igm_absorb.f90: 17 Lyman-series
transitions (Ly-α through Ly-17) with blanketing index 3.46,
metal blanketing on Ly-α, and the analytic Lyman-continuum
approximation (Eq. 16 of Madau 1995). τ is capped at its peak
short-wavelength value so the fitting-function breakdown below
~100 Å does not produce a non-monotonic curve.
:class:NoIGM
Identity (τ ≡ 0). Useful as a null model so the forward
pass does not need a Python-level branch on self.igm_model is
None.
Extensibility
Adding e.g. Inoue et al. (2014) is a matter of subclassing
:class:IGMModel, implementing tau, and registering in
:data:_MODEL_REGISTRY. The public constructor
:func:make_igm_model accepts either a string name or a pre-built
instance:
from ceridwen.igm import make_igm_model, Madau1995 igm = make_igm_model("madau1995") igm = make_igm_model(Madau1995()) igm = make_igm_model(None) # NoIGM
Tied to CSPBasis
CSPBasis accepts add_igm and igm_model kwargs. When
add_igm=True, CSPBasis.predict multiplies the post-mass-
post-flux-factor spectrum by the transmission curve at
theta["zred"] (and, optionally, theta["igm_factor"]) — so IGM
strength can be fit.
References
Madau, P. (1995), ApJ, 441, 18. "Radiative Transfer in a Clumpy Universe: The Colors of High-Redshift Galaxies".
IGMModel ¶
Bases: ABC
Abstract IGM attenuation model.
Concrete subclasses implement tau(lam_rest, zred). The
:meth:attenuation method returns exp(-tau * factor) and is
what callers (e.g. CSPBasis.predict) multiply into the
spectrum.
Madau1995 ¶
Bases: IGMModel
Madau (1995) IGM attenuation.
Direct port of FSPS's igm_absorb.f90. Three components:
- Lyman-series line blanketing (17 transitions from Ly-α down to Ly-17):
.. math:: \tau_\mathrm{Ly\,series}(\lambda_\mathrm{rest}) = \sum_{i=1}^{17} A_i \left(\frac{\lambda_\mathrm{rest}(1+z)}{\lambda_i}\right)^{3.46},
applied only where :math:\lambda_\mathrm{rest} < \lambda_i.
- Metal blanketing on Ly-α with a softer index 1.68:
.. math:: \tau_\mathrm{metal} = 0.0017 \bigl(\lambda_\mathrm{obs}/1215.67\bigr)^{1.68}.
- Lyman-continuum photoelectric absorption (Madau 1995 Eq. 16
analytic approximation) for :math:
\lambda_\mathrm{rest} < 911.75\,\text{Å}.
A small post-processing step caps τ at its peak short-wavelength value, suppressing the fitting-function breakdown at ~100 Å.
At zred = 0 τ is forced to zero — the Madau 1995 approximation
is not physically valid for a line-of-sight of zero length and
would otherwise give spurious residual attenuation to local-
universe sources. This matches FSPS's own usage, which only
invokes IGM_ABSORB when zred > 0.
make_igm_model ¶
Build an :class:IGMModel from a string name, an instance, or
None.
Parameters:
| Name | Type | Description | Default |
|---|---|---|---|
name_or_model
|
str, IGMModel, or None
|
String name (key of :data: |
required |
Examples:
>>> make_igm_model("madau1995") # -> Madau1995()
>>> make_igm_model(None) # -> NoIGM()
>>> make_igm_model(Madau1995()) # -> same instance
Source code in ceridwen/igm.py
ceridwen.cosmology ¶
ceridwen/cosmology.py¶
Redshift helpers for ceridwen, with two back-ends:
-
JAX native (default, always available) — flat LambdaCDM, matter only, Planck-2018-like parameters, Simpson integrator on a fixed 128-node grid. JIT-compilable, differentiable in z — required whenever zred is a sampled free parameter (NUTS needs gradients).
-
astropy (optional, auto-detected) — uses
astropy.cosmology.Planck18.luminosity_distanceand therefore includes radiation + massive neutrinos, closing the ~0.5-1% gap between the native integrator and published Planck18 tables. astropy's code is NumPy-based, not JAX — so it can only be used whenzis a concrete Python scalar (not a traced value), typically for the one-off precompute step atSedModel(csp, obs, priors, zred=0.5). Setbackend='astropy'in the top-level helpers to opt in explicitly.
References
Hogg, "Distance measures in cosmology" (astro-ph/9905116). Planck Collaboration 2020, A&A 641, A6 (Planck 2018 VI). astropy Cosmology docs: https://docs.astropy.org/en/stable/cosmology/
Conventions
- distances in megaparsec (Mpc); luminosity distance
D_L(z)in Mpc. - Native back-end integrates :math:
\int_0^z dz'/E(z')with Simpson's rule on a fixedn_nodesgrid — fully differentiable, constant XLA graph cost independent ofz.
Use
from ceridwen.cosmology import luminosity_distance_mpc luminosity_distance_mpc(0.5) # JAX, any z luminosity_distance_mpc(0.5, backend='astropy') # Planck18 exact, scalar z only
Cosmology
dataclass
¶
Flat LambdaCDM cosmology with photons + effective neutrinos.
This is a deliberately-minimal JAX-friendly port of astropy's
FlatLambdaCDM. The massive-neutrino contribution is folded into
Om0 as cold matter (the late-time approximation; valid for
m_nu <= 1 eV at z <= 10 to <0.1% on D_L), so we do NOT need the
full relativistic transition that astropy evaluates point-by-point.
All properties are Python scalars and constant-fold into the XLA graph when this object is passed as a closure to a jitted function.
Onu0_massive_as_matter
property
¶
Massive-neutrino density today treated as cold matter
(Mangano et al. 2005 relation :math:\Omega_\nu h^2 = \sum m_\nu/93.14\,\mathrm{eV}).
Onu0_relativistic
property
¶
Relativistic (massless-equivalent) neutrino density today.
For Planck18 with one 0.06 eV neutrino, Neff = 3.046 counts the total at CMB epoch but one species becomes non-relativistic at late times. Subtract one for the late-time relativistic count.
E_of_z ¶
Dimensionless expansion rate
.. math:: E(z) = \sqrt{\Omega_m^\mathrm{eff}(1+z)^3 + \Omega_r (1+z)^4 + \Omega_\Lambda}.
Includes radiation (photons + massless neutrinos) exactly; the
single massive neutrino is folded into :math:\Omega_m^\mathrm{eff}
as cold matter. This matches
:class:astropy.cosmology.FlatLambdaCDM with the same cosmological
parameters to < 0.1% on :math:D_L(z) for :math:z \le 10.
The input is clamped to z >= 0 before use. VI ELBO training
and NUTS leapfrog steps can transiently explore z < 0 while
the proposal is being built, and at z < -1 the radicand
:math:\Omega_m(1+z)^3 + \Omega_r(1+z)^4 + \Omega_\Lambda can go
negative, giving sqrt(negative) = NaN. A single NaN in the
ELBO gradient poisons the entire TriL map and causes 100%
divergences downstream. Clamping at 0 keeps the trace finite;
the proposal is still rejected by whatever prior / transform is
in play.
The formula is identical to astropy's efunc implementation
modulo the massive-neutrino approximation; structurally
.. math:: E(z)^2 = (1+z)^3\bigl[\Omega_r(1+z) + \Omega_m^\mathrm{eff}\bigr] + \Omega_\Lambda,
which is the form astropy uses to factor out one (1+z)^3 for
numerical stability at high z.
Source code in ceridwen/cosmology.py
comoving_distance_mpc ¶
Comoving (line-of-sight) distance
:math:D_C(z) = D_H \int_0^z dz'/E(z').
Source code in ceridwen/cosmology.py
age_gyr ¶
Age of the Universe at redshift z in Gyr (JAX-native).
.. math:: t(z) = \frac{1}{H_0}\int_0^{a(z)}\frac{da'}{a'\,E(z'(a'))}, \qquad a = \frac{1}{1+z},\; z'(a') = 1/a' - 1.
Simpson's rule on a fixed scale-factor grid (traced shape static), so
the result is differentiable in z and usable inside a JIT'd model
(e.g. to recompute SFH age-bins from a sampled redshift). The
integrand vanishes analytically as :math:a'\to 0 (radiation era),
so the a'=0 node is set to zero to avoid 0/0. Matches
:meth:astropy.cosmology.FlatLambdaCDM.age to <~0.3% over
:math:0 \le z \le 20 (residual dominated by the massive-neutrino
approximation in :func:E_of_z, not the integrator).
Source code in ceridwen/cosmology.py
luminosity_distance_mpc ¶
Luminosity distance :math:D_L(z) = (1+z) D_C(z) in Mpc.
Parameters:
| Name | Type | Description | Default |
|---|---|---|---|
z
|
Array or float
|
Redshift. JAX arrays are supported by the native backend; scalars only for the astropy backend. |
required |
cosmo
|
Cosmology
|
Only used by the native backend (astropy uses Planck18 internally). |
DEFAULT_COSMO
|
n_nodes
|
int
|
Simpson-integrator node count (native only). 128 is accurate to ~1e-6 for the integral itself; the residual error vs astropy is dominated by missing radiation + neutrinos, not the integrator. |
128
|
backend
|
('native', 'astropy')
|
|
'native'
|
Source code in ceridwen/cosmology.py
flux_factor ¶
Cosmological flux factor to turn rest-frame per-Hz luminosity into
observed-frame flux density, including the :math:(1+z) term for
:math:F_\nu.
Specifically, for a rest-frame luminosity per unit frequency
:math:L_\nu^\mathrm{rest} in erg s^{-1} Hz^{-1} one has
.. math:: F_\nu^\mathrm{obs}(\nu_\mathrm{obs}) = \frac{(1+z)\, L_\nu^\mathrm{rest}(\nu_\mathrm{rest})} {4 \pi\, D_L(z)^2}.
This helper returns :math:(1+z) / (4\pi D_L^2) in CGS-compatible units
such that multiplying a spectrum in erg s^{-1} Hz^{-1} yields the
observed-frame :math:F_\nu in erg s^{-1} cm^{-2} Hz^{-1}.
For the ceridwen "maggies"-style convention where the rest-frame
spectrum is in L_sun Hz^{-1} M_sun^{-1} (as produced by FSPS / SSP
tables), this factor should be composed with any unit conversion
the forward model already does. See
:func:flux_factor_maggies for the maggies-ready version.
Source code in ceridwen/cosmology.py
flux_factor_maggies ¶
Redshift scaling relative to the CSP's fiducial d = 10 pc.
Multiplies a CSPBasis-computed (rest-frame, per-source) spectrum to
produce the observed-frame spectrum at the luminosity distance for
redshift z. At z = 0 the CSP normalisation already corresponds
to 10 pc, so :math:\mathrm{ff}(0) = 1. For z > 0:
.. math:: \mathrm{ff}(z) = \frac{(1 + z)\, D_L(z=0)^2}{D_L(z)^2} = \frac{(1 + z)\,(10\,\mathrm{pc})^2}{D_L(z)^2}.
This assumes the observation flux unit is maggies relative to the same absolute magnitude convention.
Parameters:
| Name | Type | Description | Default |
|---|---|---|---|
z
|
Array or float
|
Redshift. |
required |
backend
|
('native', 'astropy')
|
See :func: |
'native'
|
Source code in ceridwen/cosmology.py
have_astropy ¶
True if astropy is importable. Used by SedModel to decide whether to use the higher-accuracy Planck18 backend for fixed-z precompute.