From d9345ea216137953aebe58d2ae2e8be0019f1734 Mon Sep 17 00:00:00 2001 From: Pablo Corcho-Caballero Date: Thu, 11 Jun 2026 12:07:48 +0200 Subject: [PATCH 1/9] store edges instead of bins --- src/besta/postprocess.py | 47 +++++++++++++++++++++++++++++++--------- 1 file changed, 37 insertions(+), 10 deletions(-) diff --git a/src/besta/postprocess.py b/src/besta/postprocess.py index 6b5bfe6b..2ce87818 100644 --- a/src/besta/postprocess.py +++ b/src/besta/postprocess.py @@ -350,24 +350,24 @@ def histogram_pdf_1d( # Convert probability per bin -> density with np.errstate(divide="ignore", invalid="ignore"): pdf = hist / dx - centers = 0.5 * (edges[:-1] + edges[1:]) # Ensure integrates to 1 (numerical) integral = np.nansum(pdf * dx) if integral > 0: pdf /= integral - return centers, pdf + return edges, pdf def kde_pdf_1d( x: np.ndarray, weights: np.ndarray, - grid: np.ndarray, + bins: np.ndarray, ) -> np.ndarray: """ KDE PDF on a provided grid; returns NaNs on failure. """ x = _as_float_array(x).ravel() w = _as_float_array(weights).ravel() - g = _as_float_array(grid).ravel() + edges = _as_float_array(bins).ravel() + g = 0.5 * (edges[:-1] + edges[1:]) mask = np.isfinite(x) & np.isfinite(w) x = x[mask] w = w[mask] @@ -450,6 +450,33 @@ def kde_or_hist_pdf_2d( Z /= integral return xc, yc, Z +def pit_from_pdf(x_edges, pdf, x_true): + """ + Compute PIT value for a true value given a PDF defined by edges and values. + + Parameters + ---------- + x_edges : array shape (N+1,) bin edges + pdf : array shape (N,) PDF values for each bin, normalized to integrate to 1 + x_true : scalar true value + + Returns + ------- + pit : scalar in [0,1] representing the cumulative probability up to ``x_true`` + """ + x_edges = _as_float_array(x_edges).ravel() + pdf = _as_float_array(pdf).ravel() + if x_edges.size != pdf.size + 1: + raise ValueError("x_edges must have one more element than pdf.") + if not np.isfinite(x_true): + raise ValueError("x_true must be finite.") + dx = np.diff(x_edges) + cdf = np.cumsum(pdf * dx) + if not np.isclose(cdf[-1], 1.0): + raise ValueError("pdf must be normalized to integrate to 1.") + + return np.interp(x_true, x_edges, np.r_[0.0, cdf]) + # ----------------------------------------------------------------------------- # Autocorrelation # ----------------------------------------------------------------------------- @@ -939,12 +966,12 @@ def to_fits(self) -> fits.HDUList: hdus.append(fits.BinTableHDU(t_pct, name="PERCENTILES", header=pct_hdr)) - # PDF1D table: store grid/pdf/kde per parameter as separate columns + # PDF1D table: store edges/pdf/kde per parameter as separate columns t_pdf1 = Table() for name, d in self.pdf_1d.items(): - t_pdf1[f"{name}_x"] = _as_float_array(d["grid"]) + t_pdf1[f"{name}_x"] = _as_float_array(d["edges"]) t_pdf1[f"{name}_pdf"] = _as_float_array(d["hist_pdf"]) - t_pdf1[f"{name}_kde"] = _as_float_array(d.get("kde_pdf", np.full_like(d["grid"], np.nan))) + t_pdf1[f"{name}_kde"] = _as_float_array(d.get("kde_pdf", np.full_like(d["edges"], np.nan))) if len(t_pdf1.colnames) > 0: hdus.append(fits.BinTableHDU(t_pdf1, name="PDF1D")) @@ -1274,10 +1301,10 @@ def summarize_results( pdf1d: Dict[str, Dict[str, np.ndarray]] = {} if compute_1d: for i, nm in enumerate(names): - grid, hist_pdf = histogram_pdf_1d(samples[i, :], w, bins=pdf_bins_1d) - d = {"grid": grid, "hist_pdf": hist_pdf} + edges, hist_pdf = histogram_pdf_1d(samples[i, :], w, bins=pdf_bins_1d) + d = {"edges": edges, "hist_pdf": hist_pdf} if kde_1d: - d["kde_pdf"] = kde_pdf_1d(samples[i, :], w, grid) + d["kde_pdf"] = kde_pdf_1d(samples[i, :], w, edges) pdf1d[nm] = d # 2D PDFs From 6ef6e9339382cf8285ac15c730b433d01c05943a Mon Sep 17 00:00:00 2001 From: Pablo Corcho-Caballero Date: Thu, 11 Jun 2026 12:09:31 +0200 Subject: [PATCH 2/9] update notebooks --- tutorials/fit_redshift/fit_jwst_redshift.ipynb | 4 ++-- tutorials/fit_sdss_spectra/sdss_spectra.ipynb | 10 ++++++---- 2 files changed, 8 insertions(+), 6 deletions(-) diff --git a/tutorials/fit_redshift/fit_jwst_redshift.ipynb b/tutorials/fit_redshift/fit_jwst_redshift.ipynb index 2f9e5fad..5c90498f 100644 --- a/tutorials/fit_redshift/fit_jwst_redshift.ipynb +++ b/tutorials/fit_redshift/fit_jwst_redshift.ipynb @@ -521,7 +521,7 @@ ], "metadata": { "kernelspec": { - "display_name": "cosmo-env", + "display_name": "besta", "language": "python", "name": "python3" }, @@ -535,7 +535,7 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.11.10" + "version": "3.14.4" } }, "nbformat": 4, diff --git a/tutorials/fit_sdss_spectra/sdss_spectra.ipynb b/tutorials/fit_sdss_spectra/sdss_spectra.ipynb index 8c9f34d0..ca87a29f 100644 --- a/tutorials/fit_sdss_spectra/sdss_spectra.ipynb +++ b/tutorials/fit_sdss_spectra/sdss_spectra.ipynb @@ -191,9 +191,11 @@ "metadata": {}, "outputs": [], "source": [ - "sdss_lsf_file = \"lsf_SDSS.dat\"\n", + "lsf_pixels = hdul[1].data[\"wdisp\"]\n", + "sdss_lsf_wl = lsf_pixels * np.gradient(wl.value)\n", + "sdss_lsf_fwhm = wl\n", + "# convert to FWHM in Angstrom at each wavelength accounting for each pixels varying dispersion\n", "emiles_lsf_file = \"e-miles_spectral_resolution.dat\"\n", - "sdss_lsf_wl, sdss_lsf_fwhm = np.loadtxt(sdss_lsf_file, unpack=True)\n", "emiles_lsf_wl, emiles_lsf_fwhm = np.loadtxt(emiles_lsf_file, unpack=True, usecols=(0, 1))\n", "\n", "wl_min = 3850\n", @@ -969,7 +971,7 @@ ], "metadata": { "kernelspec": { - "display_name": "Python 3", + "display_name": "besta", "language": "python", "name": "python3" }, @@ -983,7 +985,7 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.11.10" + "version": "3.14.4" } }, "nbformat": 4, From b1158cd6025ec40e6b8265ea9c65b7ffe306ad64 Mon Sep 17 00:00:00 2001 From: Pablo Corcho-Caballero Date: Sat, 13 Jun 2026 12:05:04 +0200 Subject: [PATCH 3/9] remove stales imports --- src/besta/grid/grid.py | 3 --- 1 file changed, 3 deletions(-) diff --git a/src/besta/grid/grid.py b/src/besta/grid/grid.py index fbaa31c6..766e154c 100644 --- a/src/besta/grid/grid.py +++ b/src/besta/grid/grid.py @@ -35,13 +35,10 @@ from besta.postprocess import ( enclosed_fraction_map, - pit_from_discrete_posterior, pdf_stats, - photoz_metrics, weighted_quantiles, ) -from besta.utils import available_memory_bytes from besta.logging import get_logger os.environ.setdefault("OMP_NUM_THREADS", "1") From a40124fdb3a4a1d4625b3a79a54da3c2a0003335 Mon Sep 17 00:00:00 2001 From: Pablo Corcho-Caballero Date: Sat, 13 Jun 2026 15:30:20 +0200 Subject: [PATCH 4/9] add test feature masking --- tests/test_run.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/tests/test_run.py b/tests/test_run.py index 3c2d8fb9..27bcd56d 100644 --- a/tests/test_run.py +++ b/tests/test_run.py @@ -111,6 +111,7 @@ def test_fit(self): "SFHModel": "ExponentialSFH", "velscale": 50.0, "ExtinctionLaw": "ccm89", + "use_features": "T", }} t0 = time() @@ -140,4 +141,4 @@ def test_postprocess(self): results = summarize_results(results.results_table) if __name__ == "__main__": - unittest.main() \ No newline at end of file + unittest.main() From efb3b3c4ac1262bc4ddbac62294899efcdcae81e Mon Sep 17 00:00:00 2001 From: Pablo Corcho-Caballero Date: Sat, 13 Jun 2026 15:31:09 +0200 Subject: [PATCH 5/9] allow for feature-based weights --- src/besta/pipeline_modules/base_module.py | 18 +++++++++++++----- 1 file changed, 13 insertions(+), 5 deletions(-) diff --git a/src/besta/pipeline_modules/base_module.py b/src/besta/pipeline_modules/base_module.py index b7b0683a..5e29fe08 100644 --- a/src/besta/pipeline_modules/base_module.py +++ b/src/besta/pipeline_modules/base_module.py @@ -655,6 +655,9 @@ def prepare_observed_spectra( if not (instrumental_lsf == 0).all(): self.config["lsf"] = instrumental_lsf + if options.get_bool("use_features", default=True): + self.get_feature_weights(options) + self.config["weights"] *= self.config["feature_weights"] _log("Configuration done.") def prepare_galaxy(self, options): @@ -804,6 +807,7 @@ def prepare_losvd_kernel(self, options): _log("Configuration done") def get_feature_weights(self, options): + """TODO""" logger.info("Computing feature weights from input spectra") # Estimate the continuum continuum, continuum_err = spectrum.estimate_continuum( @@ -817,12 +821,10 @@ def get_feature_weights(self, options): self.config["continuum"] = continuum self.config["continuum_err"] = continuum_err # Favour features over/under continuum - w = (np.abs(self.config["flux"] - continuum) / continuum_err)**2 + weight_powlaw = options.get_double("feature_weight_powlaw", 2.0) + w = (np.abs(self.config["flux"] - continuum) / continuum_err)**weight_powlaw w = np.where(np.isfinite(w), w, 0.0) - w_sum = np.nansum(w) - if w_sum <= 0: - raise ValueError("Feature-based weights sum to zero; please check the input data or disable feature-based weighting.") - w /= w_sum + w /= w.max() self.config["feature_weights"] = w def measure_emission_lines(self, solution: DataBlock, **kwargs): @@ -1047,6 +1049,12 @@ def plot_solution(self, solution: DataBlock, figname=None, plot_lines=True): ax.set_yscale("symlog", linthresh=1.0) ax.set_xlabel("Wavelength (AA)") + twax = ax.twinx() + twax.fill_between( + np.array(self.config["wavelength"]), 0.0, weights, + color="lime", alpha=0.2, label="Weights") + twax.set_ylabel("Weight") + ax = axs[1, 1] ax.hist( chi2, From ba792882ba7a9879be11b97efcdca9b959e1d38b Mon Sep 17 00:00:00 2001 From: Pablo Corcho-Caballero Date: Sat, 13 Jun 2026 15:31:59 +0200 Subject: [PATCH 6/9] write ini values file when use_transforms=True --- src/besta/sfh.py | 20 +++++++++++++++----- 1 file changed, 15 insertions(+), 5 deletions(-) diff --git a/src/besta/sfh.py b/src/besta/sfh.py index a7830045..43749122 100644 --- a/src/besta/sfh.py +++ b/src/besta/sfh.py @@ -97,9 +97,21 @@ def make_ini(self, ini_file): Path to the output .ini file. """ logger.info("Making ini file: %s", ini_file) + + free_params = self.free_params.copy() + + if self.use_transforms: + logger.info("transforming default SFH values into latent variables") + if getattr(self, "sfh_bin_keys", None): + sfh_values = self.to_latent([free_params[k][1] for k in self.sfh_bin_keys]) + for key, val in zip(self.sfh_bin_keys, sfh_values): + free_params[key] = [-5, val, 5] + with open(ini_file, "w", encoding="utf-8") as file: + file.write(f"; Default prior file for SFH model: {str(self.__class__)}\n") + file.write(f"; use_transforms: {str(self.use_transforms)}\n") file.write(f"[{self.sect_name}]\n") - for key, val in self.free_params.items(): + for key, val in free_params.items(): if len(val) > 1: file.write(f"{key} = {val[0]} {val[1]} {val[2]}\n") else: @@ -169,15 +181,13 @@ class FixedTimeSFH(ZPowerLawMixin, SFHBase, PieceWiseSFHMixin): The SFH of a galaxy is modelled as a stepwise function where the free parameters correspond to the mass fraction formed on each bin. - Upon initializaiton, the free parameters are set between -8 to 0 in terms - of log(M/Msun). The starting point corresponds to the mass fraction formed - assuming a constant star formation history. + The default starting point for uniform priors corresponds to the mass + fraction formed assuming a constant star formation history. Attributes ---------- lookback_time : astropy.units.Quantity Lookback time bin edges. - """ def __init__(self, lookback_time_bins, *args, **kwargs): From 3435e268c45da126fff3572eca631c41f72325c1 Mon Sep 17 00:00:00 2001 From: Pablo Corcho-Caballero Date: Mon, 15 Jun 2026 12:08:49 +0200 Subject: [PATCH 7/9] re-organise methods --- src/besta/postprocess.py | 182 ++++++++++++++++----------------------- 1 file changed, 74 insertions(+), 108 deletions(-) diff --git a/src/besta/postprocess.py b/src/besta/postprocess.py index 2ce87818..d4f90e44 100644 --- a/src/besta/postprocess.py +++ b/src/besta/postprocess.py @@ -477,6 +477,65 @@ def pit_from_pdf(x_edges, pdf, x_true): return np.interp(x_true, x_edges, np.r_[0.0, cdf]) +def pdf_stats(edges: np.ndarray, pdf: np.ndarray, + quantiles=None, + find_multimodal: bool = False) -> dict: + """ + Compute summary stats from a discrete pdf over bin centers. + + Parameters + ---------- + edges : ndarray, shape (K,) + Bin edges. + pdf : ndarray, shape (K,) + PDF defined by the edges. + quantiles : tuple of float, optional + Quantiles to report (default 16, 50, 84 percent). + find_multimodal : bool, optional + If True, return rough modes by local-maximum search. + + Returns + ------- + stats : dict + Keys: mean, std, map, q, lo68, hi68, modes (optional). + """ + # Ensure normalization + norm = np.sum(pdf * np.diff(edges)) + pdf /= norm if norm > 0 else 1.0 + centers = 0.5 * (edges[:-1] + edges[1:]) + # Trapecium integrals + mean = np.sum(pdf * centers * np.diff(edges)) + var = np.sum(pdf * (centers - mean)**2 * np.diff(edges)) + std = var ** 0.5 + k_map = np.argmax(pdf * np.diff(edges)) + v_map = centers[k_map] + + cdf = np.cumsum(pdf * np.diff(edges)) + cdf = np.insert(cdf, 0, 0) + + if quantiles is not None: + qs = np.array(quantiles, float) + qvals = np.interp(qs, cdf, edges, left=edges[0], right=edges[-1]) + else: + qvals = None + + median = np.interp(0.5, cdf, edges, left=edges[0], right=edges[-1]) + lo68 = np.interp(0.16, cdf, edges, left=edges[0], right=edges[-1]) + hi68 = np.interp(0.84, cdf, edges, left=edges[0], right=edges[-1]) + + out = {"mean": mean, "std": std, "map": v_map, "q": qvals, + "median": median, "lo68": lo68, "hi68": hi68} + + if find_multimodal: + modes = [] + for i in range(1, len(pdf) - 1): + if pdf[i] > pdf[i - 1] and pdf[i] > pdf[i + 1]: + modes.append(centers[i]) + if not modes: + modes = [v_map] + out["modes"] = np.asarray(modes) + return out + # ----------------------------------------------------------------------------- # Autocorrelation # ----------------------------------------------------------------------------- @@ -890,6 +949,11 @@ def _tolist(a): } return out + @classmethod + def from_json(cls): + #TODO + pass + def write_json(self, path: str, *, overwrite: bool = True, indent: int = 2) -> str: """Write a JSON summary file.""" if (not overwrite) and os.path.exists(path): @@ -994,9 +1058,10 @@ def write_fits(self, path: str, *, overwrite: bool = True) -> str: hdul.writeto(path, overwrite=overwrite) return path - # ------------------------------------------------------------------------- - # Plot helpers (optional) - # ------------------------------------------------------------------------- + @classmethod + def from_fits(cls): + #TODO + pass def plot_1d_pdfs( self, @@ -1090,15 +1155,7 @@ def corner_plot( dpi: int = 200, show: bool = False, ) -> str: - """ - Lightweight corner plot without external dependencies. - - Diagonal: 1D hist PDFs (weighted) - Off-diagonal: scatter of (subsampled) points colored by weight rank (simple) - - For serious usage, consider adding an optional dependency later (corner/arviz), - but this is a decent built-in baseline. - """ + """Build a corner plot""" npar = len(self.parameter_names) if npar == 0 or self.samples.size == 0: raise ValueError("No samples available to plot.") @@ -1114,7 +1171,10 @@ def corner_plot( S = self.samples[:, idx] w = self.weights[idx] - fig, axes = plt.subplots(npar, npar, figsize=(2.2 * npar, 2.2 * npar), constrained_layout=True) + fig, axes = plt.subplots( + npar, npar, figsize=(2.2 * npar, 2.2 * npar), + sharex="col", + constrained_layout=True) for i in range(npar): for j in range(npar): @@ -1128,7 +1188,7 @@ def corner_plot( ax.axvline(self.map[i], lw=1.0, alpha=0.8) ax.set_yticks([]) elif i > j: - ax.scatter(S[j, :], S[i, :], s=2, alpha=0.25) + ax.hist2d(S[j, :], S[i, :], weights=w) else: ax.axis("off") @@ -1503,100 +1563,6 @@ def weighted_quantiles(x: np.ndarray, w: np.ndarray, cdf = cdf / cdf[-1] return np.interp(qs, cdf, x) -def pdf_stats(edges: np.ndarray, pdf: np.ndarray, - quantiles=None, - find_multimodal: bool = False) -> dict: - """ - Compute summary stats from a discrete pdf over bin centers. - - Parameters - ---------- - edges : ndarray, shape (K,) - Bin edges. - pdf : ndarray, shape (K,) - PDF defined by the edges. - quantiles : tuple of float, optional - Quantiles to report (default 16, 50, 84 percent). - find_multimodal : bool, optional - If True, return rough modes by local-maximum search. - - Returns - ------- - stats : dict - Keys: mean, std, map, q, lo68, hi68, modes (optional). - """ - # Ensure normalization - norm = np.sum(pdf * np.diff(edges)) - pdf /= norm if norm > 0 else 1.0 - centers = 0.5 * (edges[:-1] + edges[1:]) - # Trapecium integrals - mean = np.sum(pdf * centers * np.diff(edges)) - var = np.sum(pdf * (centers - mean)**2 * np.diff(edges)) - std = var ** 0.5 - k_map = np.argmax(pdf * np.diff(edges)) - v_map = centers[k_map] - - cdf = np.cumsum(pdf * np.diff(edges)) - cdf = np.insert(cdf, 0, 0) - - if quantiles is not None: - qs = np.array(quantiles, float) - qvals = np.interp(qs, cdf, edges, left=edges[0], right=edges[-1]) - else: - qvals = None - - median = np.interp(0.5, cdf, edges, left=edges[0], right=edges[-1]) - lo68 = np.interp(0.16, cdf, edges, left=edges[0], right=edges[-1]) - hi68 = np.interp(0.84, cdf, edges, left=edges[0], right=edges[-1]) - - out = {"mean": mean, "std": std, "map": v_map, "q": qvals, - "median": median, "lo68": lo68, "hi68": hi68} - - if find_multimodal: - modes = [] - for i in range(1, len(pdf) - 1): - if pdf[i] > pdf[i - 1] and pdf[i] > pdf[i + 1]: - modes.append(centers[i]) - if not modes: - modes = [v_map] - out["modes"] = np.asarray(modes) - return out - - -def pit_from_discrete_posterior(z_true: np.ndarray, - posts: np.ndarray, - z_edges: np.ndarray) -> np.ndarray: - """ - Probability Integral Transform for discrete posteriors on bins. - - Assumes uniform density within each bin for within-bin interpolation. - - Parameters - ---------- - z_true : ndarray, shape (N,) - True values. - posts : ndarray, shape (N, K) - Row-normalised posteriors over K bins. - z_edges : ndarray, shape (K+1,) - Bin edges. - - Returns - ------- - pit : ndarray, shape (N,) - PIT values in [0, 1]. - """ - N, K = posts.shape - assert K == len(z_edges) - 1 - cdf_bins = np.cumsum(posts, axis=1) - j = np.clip(np.digitize(z_true, z_edges) - 1, 0, K - 1) - idx = np.arange(N) - below = np.where(j > 0, cdf_bins[idx, j - 1], 0.0) - widths = z_edges[1:] - z_edges[:-1] - frac = np.clip((z_true - z_edges[j]) / widths[j], 0.0, 1.0) - pit = below + posts[idx, j] * frac - return np.clip(pit, 0.0, 1.0) - - def photoz_metrics(z_true: np.ndarray, z_est: np.ndarray) -> dict: """ Standard photo-z metrics using delta z over 1+z. From bb4c3e6195b90296ce04606acebe24c5d8e517a6 Mon Sep 17 00:00:00 2001 From: Pablo Corcho-Caballero Date: Mon, 15 Jun 2026 14:20:57 +0200 Subject: [PATCH 8/9] rename test --- tests/_test_pymc_sampler.py | 250 ++++++++++++++++++++++++++++++++++++ 1 file changed, 250 insertions(+) create mode 100644 tests/_test_pymc_sampler.py diff --git a/tests/_test_pymc_sampler.py b/tests/_test_pymc_sampler.py new file mode 100644 index 00000000..9a0b04f4 --- /dev/null +++ b/tests/_test_pymc_sampler.py @@ -0,0 +1,250 @@ +from __future__ import annotations + +import shutil +from pathlib import Path +from types import SimpleNamespace + +import numpy as np +import pytest + +from cosmosis.output import InMemoryOutput + +from besta import io +from besta.pipeline import MainPipeline +from besta.samplers.pymc_sampler import PymcSampler + + +def _write_ini(path: Path, text: str) -> Path: + path.write_text(text, encoding="utf-8") + return path + + +def test_pymc_sampler_smoke_with_dummy_pipeline(tmp_path): + class DummyPipeline: + def __init__(self): + self.varied_params = [object(), object()] + self.modules = [] + + def output_names(self): + # Varied params + one extra output. + return ["p0", "p1", "extra0"] + + def denormalize_vector(self, u): + # Map unit-cube to a bounded physical space. + return 4.0 * np.asarray(u, dtype=float) - 2.0 + + def run_results(self, p): + p = np.asarray(p, dtype=float) + post = -0.5 * float(np.dot(p, p)) + prior = -0.1 * float(np.dot(p, p)) + like = post - prior + extra = np.array([float(p.sum())], dtype=float) + return SimpleNamespace(post=post, prior=prior, like=like, extra=extra) + + ini_path = _write_ini( + tmp_path / "pymc_smoke.ini", + """ +[runtime] +sampler = pymc + +[pymc] +samples = 24 +nsteps = 8 +burn_fraction = 0.0 +chains = 1 +progressbar = F +seed = 11 +""".strip(), + ) + + output = InMemoryOutput() + sampler = PymcSampler(str(ini_path), DummyPipeline(), output) + sampler.config() + sampler.execute() + + assert sampler.is_converged() + assert sampler.num_samples > 0 + assert len(output.rows) == sampler.num_samples + + col_names = [c[0] for c in output.columns] + assert "prior" in col_names + assert "post" in col_names + assert "like" in col_names + + +def test_pymc_sampler_accepts_slice_step_method(tmp_path): + class DummyPipeline: + def __init__(self): + self.varied_params = [object(), object()] + self.modules = [] + + def output_names(self): + return ["p0", "p1", "extra0"] + + def denormalize_vector(self, u): + return 4.0 * np.asarray(u, dtype=float) - 2.0 + + def run_results(self, p): + p = np.asarray(p, dtype=float) + post = -0.5 * float(np.dot(p, p)) + prior = -0.1 * float(np.dot(p, p)) + like = post - prior + extra = np.array([float(p.sum())], dtype=float) + return SimpleNamespace(post=post, prior=prior, like=like, extra=extra) + + ini_path = _write_ini( + tmp_path / "pymc_slice.ini", + """ +[runtime] +sampler = pymc + +[pymc] +samples = 12 +nsteps = 4 +burn_fraction = 0.0 +chains = 1 +progressbar = F +seed = 5 +step_method = slice +""".strip(), + ) + + output = InMemoryOutput() + sampler = PymcSampler(str(ini_path), DummyPipeline(), output) + sampler.config() + sampler.execute() + + assert sampler.is_converged() + assert sampler.num_samples > 0 + + +def test_pymc_sampler_rejects_gradient_step_methods(tmp_path): + class DummyPipeline: + def __init__(self): + self.varied_params = [object(), object()] + self.modules = [] + + def output_names(self): + return ["p0", "p1", "extra0"] + + def denormalize_vector(self, u): + return 4.0 * np.asarray(u, dtype=float) - 2.0 + + def run_results(self, p): + p = np.asarray(p, dtype=float) + post = -0.5 * float(np.dot(p, p)) + prior = -0.1 * float(np.dot(p, p)) + like = post - prior + extra = np.array([float(p.sum())], dtype=float) + return SimpleNamespace(post=post, prior=prior, like=like, extra=extra) + + ini_path = _write_ini( + tmp_path / "pymc_nuts.ini", + """ +[runtime] +sampler = pymc + +[pymc] +samples = 6 +nsteps = 3 +burn_fraction = 0.0 +chains = 1 +progressbar = F +seed = 1 +step_method = nuts +""".strip(), + ) + + output = InMemoryOutput() + sampler = PymcSampler(str(ini_path), DummyPipeline(), output) + sampler.config() + with pytest.raises(ValueError, match="requires gradients"): + sampler.execute() + + +def test_pymc_sampler_full_besta_run(tmp_path): + cosmosis_exe = shutil.which("cosmosis") + if cosmosis_exe is None: + pytest.skip("cosmosis executable not available in PATH") + + module_path = tmp_path / "dummy_like_module.py" + module_path.write_text( + """ +def setup(options): + return {} + + +def execute(block, config): + p1 = block['parameters', 'p1'] + p2 = block['parameters', 'p2'] + like = -0.5 * (p1**2 + p2**2) + block['extra', 'sum_p'] = p1 + p2 + block['likelihoods', 'dummy_like'] = like + return 0 + + +def cleanup(config): + return 0 +""".strip(), + encoding="utf-8", + ) + + output_root = tmp_path / "pymc_full_run" + values_path = tmp_path / "values.ini" + values_path.write_text( + """ +[parameters] +p1 = -2.0 0.0 2.0 +p2 = -2.0 0.0 2.0 +""".strip(), + encoding="utf-8", + ) + + sampler_path = Path(__file__).resolve().parents[1] / "src" / "besta" / "samplers" / "pymc_sampler.py" + + config = { + "runtime": { + "sampler": "pymc", + "import_samplers": str(sampler_path), + }, + "pymc": { + "samples": 16, + "nsteps": 8, + "burn_fraction": 0.0, + "chains": 1, + "progressbar": False, + "seed": 3, + }, + "output": { + "filename": str(output_root), + "format": "text", + }, + "pipeline": { + "modules": "DummyLike", + "values": str(values_path), + "likelihoods": "dummy", + "quiet": "T", + "debug": "T", + "extra_output": "extra/sum_p", + }, + "DummyLike": { + "file": str(module_path), + }, + } + + runner = MainPipeline( + [config], + n_cores_list=[1], + ini_values_files=[str(values_path)], + ) + status = runner.execute_all(plot_result=False) + assert status == 0 + + results_file = Path(str(output_root) + ".txt") + assert results_file.exists() + + table = io.read_results_file(str(results_file)) + assert len(table) > 0 + assert "post" in table.colnames + assert "prior" in table.colnames + assert "like" in table.colnames From 75f39f5fd0cae38bdf6617f4fa8a1b3d16484f16 Mon Sep 17 00:00:00 2001 From: Pablo Corcho-Caballero Date: Mon, 15 Jun 2026 14:19:04 +0200 Subject: [PATCH 9/9] bugfix --- src/besta/postprocess.py | 24 +++++++++++++++++++----- 1 file changed, 19 insertions(+), 5 deletions(-) diff --git a/src/besta/postprocess.py b/src/besta/postprocess.py index d4f90e44..125b94da 100644 --- a/src/besta/postprocess.py +++ b/src/besta/postprocess.py @@ -1032,13 +1032,25 @@ def to_fits(self) -> fits.HDUList: # PDF1D table: store edges/pdf/kde per parameter as separate columns t_pdf1 = Table() + t_pdf1_edges = Table() for name, d in self.pdf_1d.items(): - t_pdf1[f"{name}_x"] = _as_float_array(d["edges"]) - t_pdf1[f"{name}_pdf"] = _as_float_array(d["hist_pdf"]) - t_pdf1[f"{name}_kde"] = _as_float_array(d.get("kde_pdf", np.full_like(d["edges"], np.nan))) + edges = _as_float_array(d["edges"]) + centers = _as_float_array(d.get("grid", 0.5 * (edges[:-1] + edges[1:]))) + hist_pdf = _as_float_array(d["hist_pdf"]) + kde_pdf = _as_float_array(d.get("kde_pdf", np.full_like(hist_pdf, np.nan))) + + # FITS bin table columns must have consistent lengths across rows. + # Store centers/PDF/KDE in PDF1D (all length = n_bins). + t_pdf1[f"{name}_x"] = centers + t_pdf1[f"{name}_pdf"] = hist_pdf + t_pdf1[f"{name}_kde"] = kde_pdf + # Store raw histogram edges separately (length = n_bins + 1). + t_pdf1_edges[f"{name}_edges"] = edges if len(t_pdf1.colnames) > 0: hdus.append(fits.BinTableHDU(t_pdf1, name="PDF1D")) + if len(t_pdf1_edges.colnames) > 0: + hdus.append(fits.BinTableHDU(t_pdf1_edges, name="PDF1D_EDGES")) # PDF2D images for (n0, n1), d in self.pdf_2d.items(): @@ -1181,7 +1193,8 @@ def corner_plot( ax = axes[i, j] if i == j: x = self.samples[i, :] - xc, pdf = histogram_pdf_1d(x, self.weights, bins=bins) + edges, pdf = histogram_pdf_1d(x, self.weights, bins=bins) + xc = 0.5 * (edges[:-1] + edges[1:]) ax.plot(xc, pdf, lw=1.2) # mark mean/MAP ax.axvline(self.mean[i], lw=1.0, alpha=0.8) @@ -1362,7 +1375,8 @@ def summarize_results( if compute_1d: for i, nm in enumerate(names): edges, hist_pdf = histogram_pdf_1d(samples[i, :], w, bins=pdf_bins_1d) - d = {"edges": edges, "hist_pdf": hist_pdf} + centers = 0.5 * (edges[:-1] + edges[1:]) + d = {"grid": centers, "edges": edges, "hist_pdf": hist_pdf} if kde_1d: d["kde_pdf"] = kde_pdf_1d(samples[i, :], w, edges) pdf1d[nm] = d