Skip to content

changes_tracked.md

Latest commit

 

History

History
1399 lines (1152 loc) · 73.2 KB

changes_tracked.md

File metadata and controls

1399 lines (1152 loc) · 73.2 KB
Raw

XRB Lightcurve Project — Change Log

Tracks the evolution of the IC 10 X-1 X-ray binary lightcurve simulation, fitting, and inference stack since the original R port.

Table of Contents

  1. Project Overview
  2. Phase 1 — R → Python Migration
  3. Phase 2 — Flux Integration & XSPEC Support
  4. Phase 3 — Light-Curve Data Pipeline
  5. Phase 4 — Spectral Model Analysis
  6. Phase 5 — Phase Folding & χ² Fitting
  7. Phase 6 — Eclipse Geometry & LOS Cutoff Fixes
  8. Phase 7 — Unified Wind Model
  9. Phase 8 — MCMC Pipeline & Performance
  10. Phase 9 — Convergence Improvements (Reparameterization)
  11. Phase 10 — Wind-Shape MCMC Parameters
  12. Phase 11 — N-D Precomputed Grid for Shape-Fit MCMC (later removed)
  13. Phase 12 — flux_t Errors & First-Class Unbinned MCMC
  14. Phase 13 — Wind Normalization Constants
  15. Phase 14 — Pooled Direct-Model MCMC
  16. Phase 15 — Speed/Memory Pass & BIC Replaces WAIC/LOO
  17. Phase 16 — Frozen Parameters & Kepler Mode (ParamSpec)
  18. Phase 17 — Per-Sample Phase-Shift Alignment
  19. Phase 18 — Adaptive Constant-SNR Binning & Grid Removal
  20. Phase 19 — Gaussian Phase Smoothing, Scattered Flux & Residual Panels
  21. Phase 20 — Remove Multiplicative Flux Scale from the Single-Model χ² Fit
  22. Phase 21 — MCMC Scatter-Path Audit
  23. Phase 22 — Adaptive Binning in the Single-Model CLI
  24. Phase 23 — utils/ Extraction and a Single Plotting Routine
  25. Phase 24 — Run-Config Persistence and a Replot-Mode Fix
  26. Phase 25 — MCMC Script Slimming
  27. Side Investigation — Reference Epoch Recalibration
  28. Current File Inventory
  29. Current Status & Quick Commands

Project Overview

Target: IC 10 X-1 — eclipsing X-ray binary in the Local Group galaxy IC 10.

System (working values): WR companion R ≈ 2 R☉, accretion-disk r ≈ 0.001 R☉, inclination i₀ ≈ 26°, separation d ≈ 19 R☉.

Best spectral model: TBabs × powerlaw — nH ≈ 0.75×10²² cm⁻², Γ ≈ 1.86, χ²_red ≈ 1.52 (preferred over phabs by Δχ² ≈ 8.5).

Phase 1 — R → Python Migration

Ported new11.R, grid4.R, wind_los2.R, density_fnc.R (≈260 R lines) to a single xrb_lightcurve.py with full type hints, an argparse CLI, and fully vectorized NumPy core (10–100× faster than R). Public surface: simulate_lightcurve, create_grid, wind_los_integral, density_function. Companion files: example_usage.py, plot_results.py, requirements.txt.

Phase 2 — Flux Integration & XSPEC Support

Added three flux-conversion paths selectable via --flux_method:

  • legacy — hardcoded exponential fits.
  • interpolate — log-log interpolation from an XSPEC flux vs nH CSV.
  • refit — re-fit exponentials to the XSPEC table.

Helpers: load_flux_vs_nh_csv, interpolate_flux_from_nh, fit_exponential_to_csv. CLI gained --flux_csv and --lam (target mean nH; --lam2 retired in Phase 7). XSPEC table is generated by compute_flux_vs_nH.py --specdir … --out_csv data_flux_vs_nH.csv.

Units convention: flx (atoms / R☉⁴, raw column-density integral), fl (10²² cm⁻²), nfl_{band} (photons/cm²/s).

Phase 3 — Light-Curve Data Pipeline

Chandra .txt light curves were actually FITS binary; built a small toolchain to (a) convert FITS→TXT and (b) attach a calibrated FLUX column. Tools (in utils/): convert_fits_to_txt.py and the heasoft wrapper convert_fits_to_txt_heasoft.sh, add_flux_simple.py, add_flux_to_lightcurves.py, compute_count_to_flux_factor.py, get_average_count_rates.py. Output layout: data/IC_10_X1_LC/{Broad,Soft,Hard}{_converted,_with_flux}/.

Time-averaged count rates (cts/s): broad 0.1132, soft 0.0635, hard 0.0497.

Phase 4 — Spectral Model Analysis

XSPEC tooling: xspec_get_conversion_factors_tbabs.xcm, get_conversion_factors.sh, get_xspec_nH.py, compare_absorption_models.xcm, compare_models.sh. Comparison run picked TBabs × powerlaw over phabs × powerlaw (Δχ² = 8.54 in favor of TBabs; nH 0.75 vs 0.78). 0.5–7 keV model fluxes: TBabs 1.032×10⁻¹², phabs 1.031×10⁻¹² erg/cm²/s.

Phase 5 — Phase Folding & χ² Fitting

chandra_phase_analysis.py (≈1.0k lines) folds observations onto orbital phase using the ephemeris (REF_EPOCH = 278801348 s, P = 125431 s) and fits simulation models via χ² minimization with optional phase-shift / scale rescaling. Auto-detects flux columns (nfl_*, pho_count_*), handles both standard whitespace and CIAO #Columns: formats, supports multi-column simultaneous fitting.

Phase 6 — Eclipse Geometry & LOS Cutoff Fixes

Two physics fixes plus a runtime knob in xrb_lightcurve.py:

  1. Front-emitter spurious eclipse — eclipse gating moved from gma < π to sin(gma) > 0 so emitter-in-front geometries are not incorrectly occulted.
  2. Total-eclipse flux — is_eclipsed flag added; during eclipse all nfl_* and pho_count_* columns are forced to 0 (previously they collapsed to maximum flux because flx=0 → e⁰=1).
  3. --Rmax and --converge-rmax — configurable LOS cutoff. Converged mode integrates to −∞ using a closed form for the constant-velocity term plus a cached quadrature + asymptotic tail for the accelerating-wind term, comparable in speed to the legacy Rmax = 2d heuristic. (The cached _ACCEL_* table was retired in Phase 7 once a single Numba kernel covered all profiles.)

Phase 7 — Unified Wind Model

Plan: mcmc_wind_shape_params_8b9c89d2.plan.md (precursor), implementation under unified_wind_model_77726ced.plan.md.

Goal was to drop the hardcoded AV / CV wind duality (flx2, lam2, *_cv columns) in favour of a pluggable density profile selected at runtime.

New profile registry in xrb_lightcurve.py (each exposes a dimensionless g(r) consumed inline by a single Numba kernel):

wind_model id Free shape params Notes
broken_pl 0 Rb, p Piecewise PL (kept for back-compat).
smooth_pl 1 Rb, p, Delta Default; smoothly broken PL.
beta_law 2 R_star, beta, H CAK velocity profile, g = 1/(r²·v(r)).
confinement 3 R_star, fconf, ell 1/r² with inner exp. compression.

Kernel rewrite — _wind_los_profile_numba (and a Gauss-Legendre variant _los_gl_quadrature) replaces the old converged + fixed-Rmax pair, and _simulate_phases_numba runs the full per-phase sweep under @njit(parallel=True, prange). Cached _ACCEL_U_GRID/_ACCEL_F_GRID deleted.

API changes:

  • simulate_lightcurve(...) lost lam2 and gained wind_model: str = "smooth_pl" and wind_params: dict | None.
  • Output collapsed to a single flx/fl plus one nfl_{band} column per band (no more _av/_cv suffix).
  • New helpers: pack_wind_params, default_wind_params, evaluate_g_profile, compute_surface_density(sim_df, lam, R_star, wind_model, wind_params), and wind_density_posterior(...) for converting MCMC lam posteriors into surface-density n₀ posteriors.
  • CLI: removed --lam2; added --wind-model {broken_pl, smooth_pl, beta_law, confinement} plus per-model --Rb / --p / --Delta / --beta / --H / --fconf / --ell.

Downstream cleanup — plot_results.py, chandra_phase_analysis.py, chandra_analysis_combined_flux.py, utils/test_flux_methods.py all updated to look for nfl_{band} (no _av/_cv filter).

Phase 8 — MCMC Pipeline & Performance

Plan: mcmc_performance_and_statistics_8989cd39.plan.md.

mcmc_lightcurve_fit.py (≈3.1k lines) wraps xrb_lightcurve.py in a full emcee/zeus pipeline. Headline pieces:

  • Forward-model paths
    • PrecomputedModelGrid — geometry-only 6-D grid (d1, d2, r, R, i0, phase) with vectorized RegularGridInterpolator lookup. Multi-process build via _compute_single_model + Pool. Supports --save-grid / --load-grid.
    • DirectLightCurveModel — wraps simulate_lightcurve directly; required when shape params are sampled (Phase 10).
  • Samplers — emcee (default, stretch move) and zeus (ensemble slice) via the same harness.
  • Likelihoods (--likelihood) — chi2 (Gaussian, default), jitter (Gaussian + free log_f systematic term), studentt (heavy-tailed, configurable --studentt-nu). Cash/Poisson was scoped out. (PrecomputedModelGrid was deleted in Phase 18, along with studentt; only chi2 and jitter remain.)
  • Diagnostics — ArviZ summary (r_hat, mcse_*, eti89_*), autocorrelation, optional WAIC/LOO via --compute-waic, corner plots, best-fit overlay with reduced χ². (WAIC/LOO replaced by BIC in Phase 15; --compute-waic removed in Phase 18.)

Performance work done in xrb_lightcurve.py:

  • Pre-optimization state (post Phase 7). Introducing the unified wind profile registry dropped per-LOS Python dispatch into the hot loop: the LOS integral called back into Python for g(r) at every angular cell × every step on the dz grid × every phase. Combined with the fixed-step trapezoid quadrature and per-cell eclipse checks, a single simulate_lightcurve call had ballooned to several seconds (≈3–8 s depending on dz / d2h), which made Phase 8-style MCMC (10⁴–10⁵ model evaluations) completely infeasible without a precomputed grid — and even the grid build was painfully slow.
  • Fixes (together ≈ 2 orders of magnitude):
    1. @njit(cache=True[, parallel=True]) on _g_profile, the LOS integrand, _los_gl_quadrature, _wind_los_profile_numba, and the full per-phase sweep _simulate_phases_numba (using prange).
    2. Gauss-Legendre quadrature replaces the old fixed-step trapezoid; far fewer integrand evaluations for the same accuracy, and the GL nodes/weights are precomputed once.
    3. Inlined eclipse test — eclipse gating collapsed into the kernel so eclipsed phases early-exit without a Python round-trip.
    4. All wind_params / sim params flattened to scalar JIT arguments at the Python/Numba boundary; no dict lookups in the hot path.
  • Result: a single simulate_lightcurve call is now < 1 s (≈ 63 ms) on a laptop, i.e. ~50–100× faster than the post-Phase-7 regression and 10⁵–10⁶× the sustained rate needed for MCMC. The direct evaluator became viable for MCMC without a precomputed grid, and the grid build itself dropped from minutes to seconds.

GPU acceleration was evaluated and rejected (Intel Iris ≠ CUDA / JAX-Metal / PyTorch-MPS targets). NumPyro/NUTS port was scoped out for the same reason plus the non-differentiable interpolator in the fast path.

Phase 9 — Convergence Improvements (Reparameterization)

Plan: mcmc_convergence_improvements_c648afb3.plan.md.

(d1, d2) are strongly correlated because the wind absorption sees only their sum. --reparam swaps the sampling space to:

  • a = d1 + d2 (separation, well-constrained)
  • q = d1 / (d1 + d2) (ratio, weakly constrained)

Implementation lives entirely in mcmc_lightcurve_fit.py: REPARAM_PRIORS, get_param_config(reparam=True) returns ['a','q','r','R','i0'], _evaluate_model inverts (a,q) → (d1,d2), log_prior applies priors in (a,q) space with the +log(a) Jacobian, walker init / corner / summaries report derived d1, d2. Grid build itself is unchanged (still indexed in physical d1, d2).

HMC / NUTS was assessed and not pursued — piecewise-linear RegularGridInterpolator gradients, hard eclipse branches, and Numba kernels all break autodiff. Future paths (smooth GP/NN emulator, cubic grid

  • finite differences, JAX rewrite) are documented in the plan.

Phase 10 — Wind-Shape MCMC Parameters

Plan: mcmc_wind_shape_params_8b9c89d2.plan.md.

After Phase 8 brought the direct evaluator to ~63 ms/LC, wind-shape parameters are now first-class MCMC dimensions, gated by --fit-wind-shape.

Per-model active set (in WIND_SHAPE_FIT):

--wind-model Free Fixed Tied to geometry
smooth_pl Rb, p Delta = 2.0 —
beta_law beta H = 1.0 R_star = R
confinement fconf, ell — R_star = R

broken_pl is intentionally skipped (smooth_pl generalizes it). lam (overall normalization) stays fixed from spectral fits — shape params are constrained only by the LC shape.

Mechanics in mcmc_lightcurve_fit.py:

  • New registries: WIND_MODELS, WIND_SHAPE_FIT, WIND_SHAPE_FIXED, WIND_SHAPE_LABELS, WIND_SHAPE_PRIORS, plus helpers get_active_priors(...) and _to_wind_params(theta, active_names, wind_model, R_value, fit_wind_shape).
  • get_param_config and walker init are now fully dynamic in the active_names list (no more hardcoded n_phys=5 / theta[5] / pos[:,5]; the jitter log_f index is looked up by name).
  • log_prior, log_likelihood_*, _evaluate_model, compute_chi2_for_samples, compute_pointwise_loglik, plot_best_fit, run_arviz_diagnostics, run_single_fit, and replot_from_existing all accept wind_model, fit_wind_shape, active_names and route the wind-shape sample correctly.
  • load_existing_results accepts any column set (geometry columns required, extras kept) and returns (samples, stats, loaded_names) so re-plots work for both geometry-only and shape-fit chains.
  • all_results keys switched from f"{band}_{wind_model}" to (band, wind_model) tuples (the underscore in smooth_pl etc. broke the previous rsplit('_', 1) summary writer).

CLI additions:

  • --wind-model {smooth_pl, beta_law, confinement} (default smooth_pl). The legacy av/cv/both choices and the 'both' loop are gone.
  • --fit-wind-shape — adds the wind model's active shape params to the MCMC vector. Works with either the precomputed grid (see Phase 11) or --no-grid (direct evaluator).
  • Per-shape prior overrides: --prior-Rb / -p / -beta / -fconf / -ell using the same mean,std,min,max format as the geometry priors.
  • --lam2 removed everywhere.

Smoke verified — 50-step chains for smooth_pl (geom-only and +Rb,p), beta_law (+beta), confinement (+fconf, ell), plus a grid-path geometry-only run. All produce valid corner data, ArviZ summaries, samples CSV, and mcmc_summary.txt.

Phase 11 — N-D Precomputed Grid for Shape-Fit MCMC

Superseded: the entire precomputed-grid path described below was deleted in Phase 18. MCMC now always uses the direct evaluator. Kept here for history.

After Phase 10 a short (~1k-step) shape-fit chain still took ~2 h because --fit-wind-shape auto-forced the direct evaluator (≈ 63 ms × 32 walkers × 1000 steps × all phases ≈ hours). The grid class was extended to cover wind-shape axes dynamically instead of forcing a slow path.

  • PrecomputedModelGrid is now N-D. self.axis_names holds the ordered list of axes (geometry first, then the active shape axes for the chosen wind_model). self.param_grids and self.flux_grid extend correspondingly (geometry-only grids are unchanged shape-wise, so the refactor is zero-cost for existing workflows).
  • Worker path. _precompute_models now iterates via np.ndindex over the full axis set, resolves a per-combo wind_params dict (fixed template + per-point shape values + R_star = R when the model ties it), and hands it to _compute_single_model. The r >= R geometry filter still applies.
  • Evaluation. A single RegularGridInterpolator spans (axes…, phase); evaluate() looks shape values up from the caller's wind_params dict when shape axes are present, and falls back to the fixed template for geometry-only grids.
  • I/O — .npz files persist axis_names, shape_axes, fit_wind_shape, and one <name>_grid array per axis. Legacy geometry-only .npz files (missing axis_names) still load unchanged via a backward-compat branch in _load_grid.
  • CLI — --fit-wind-shape no longer force-disables the grid; instead:
    • Grid size expands by (shape_grid_points)^k for k shape axes, with the new --shape-grid-points N knob (default 5). A [notice] prints the expected total grid size up front.
    • --save-grid / --load-grid work as usual and are now actually useful for shape-fit MCMC (build once, run many chains).
    • --no-grid stays as the escape hatch for short debugging chains.
  • Logging cleanup. load_flux_vs_nh_csv(..., verbose=verbose) and interpolate_flux_from_nh(..., warn_extrapolation=verbose) now honor the verbose flag that simulate_lightcurve already threads through. MCMC calls with verbose=False, so the per-call "Detected energy bands in CSV: …" print and repeated "nH outside CSV range / Extrapolation will be used" UserWarning are gone (they only fire in notebooks / interactive runs now). A belt-and-suspenders warnings.filterwarnings(...) in mcmc_lightcurve_fit.py keeps the extrapolation warning quiet even if a future callsite forgets to pass verbose=False.
  • Smoke verified — N-D grid build + save/load round-trip (bit-exact), end-to-end MCMC with --fit-wind-shape on smooth_pl (+ Rb, p) and beta_law (+ beta), --no-grid shape-fit path, and legacy geometry-only grid .npz backward compatibility all pass.

Phase 12 — flux_t Errors & First-Class Unbinned MCMC

Plan: flux_t_error_and_unbinned_mcmc_1bf81aa9.plan.md. Commit fa672e3.

CIAO files carry flux_t but no flux_t_err, so the MCMC's error_column = obs_column + "_ERR" guess (flux_t_ERR) never matched and silently fell back to 0.1·|flux|.

  • chandra_phase_analysis.py — new _derive_err_from_rate_err(df, obs_col), invoked only when no error column is matched. Derives per-row err = rate_err · (obs/rate) for rate > 0, falling back to a file-level cf = median(obs/rate) for the rest (verified constant per file, e.g. cf ≈ 4.56e-12 for soft 11080).
  • Master-file code removed — verify_master_contains_individual(), the master_file parameter in load_data/read_observation, and the --master-file / --verify-master CLI args are gone.
  • mcmc_lightcurve_fit.py — passes obs_error_column=None when the user didn't set it explicitly so auto-derivation fires; keeps the flux > 0 drop in both binned and unbinned modes (zero-flux rows are GTI gaps, and with jitter σ²_eff = σ_obs² + (f·model)² ≈ 0 would blow up log σ²), and prints the drop count. The invalid-error patch tightened from 0.1·|flux| to max(0.1·|flux|, median(valid obs_err)).
  • is_binned threaded from main() through run_single_fit / replot_from_existing / plot_best_fit so labels and marker style match reality ("Observed (phase-binned)" with error bars vs "Observed (raw 100s)" as translucent scatter).
  • Jitter-aware χ² — plot_best_fit and compute_chi2_for_samples report the effective-variance χ² alongside the classical measurement-error χ² when --likelihood jitter is active.
  • --no-phase-bin help now recommends pairing with --likelihood jitter, and the module docstring gained an unbinned flux_t example.

Phase 13 — Wind Normalization Constants

Plan: wind_normalization_constants_76f1ef51.plan.md. Commit fa672e3.

Since every profile is coded as a dimensionless g(r), the physical amplitude has to be recovered after the fit. Added to xrb_lightcurve.py:

  • Constants M_H_G, M_SUN_G, KM_TO_CM alongside the existing R_SUN_CM.
  • compute_wind_normalization_constants(lam, flx_mean, wind_model, wind_params, v_inf=None, mu=1.4) — computes n0 = lam·1e22/(R_sun·flx_mean) then, per model: smooth_pl → g_break, n_break_cm3, rho_b_g_cm3; beta_law/confinement → n_surface_cm3, rho_surface_g_cm3, mdot_over_vinf_g_per_cm, plus mdot_g_s / mdot_msun_yr when v_inf (km/s) is supplied. v_inf is not fitted — it was absorbed into n0 by the lam normalization — so the caller must provide it for a mass-loss rate.
  • wind_normalization_constants_posterior(...) — mirrors wind_density_posterior, looping the point estimator over posterior samples and returning {samples, median, p16, p84} per constant.

Phase 14 — Pooled Direct-Model MCMC

Commit 36c8c8b.

The direct evaluator's simulate_lightcurve is already Numba parallel=True, so naively adding --n-threads worker processes caused severe thread oversubscription (workers × numba threads ≫ cores), while pinning workers to one thread made each LC so slow that process parallelism barely beat serial.

  • _init_numba_worker(max_numba_threads) is used as the Pool initializer and calls numba.set_num_threads(...) inside each worker.
  • --numba-threads-per-worker exposes the knob; the default auto is max(1, cpu_count // n_threads), so workers collectively use ≈ one thread per logical CPU.
  • The pool uses the spawn context (safest cross-platform, avoids inheriting heavy state) and is only enabled for DirectLightCurveModel; requesting --n-threads > 1 with any other model type prints a [notice] and runs serial.
  • Also fixed an indentation bug in xrb_lightcurve.py.

Phase 15 — Speed/Memory Pass & BIC Replaces WAIC/LOO

Plan: mcmc_speed_memory_optimization_68fc497a.plan.md. Commit 7116302.

A staged, accuracy-preserving optimization program with a reproducible benchmark harness (utils/benchmark_mcmc_performance.py, documented in PERFORMANCE_VALIDATION_REPORT.md; acceptance thresholds: parameter median drift < 0.1σ, reduced-χ² change ≤ 2%, BIC ranking consistency).

Runtime (xrb_lightcurve.py):

  • Module-level _FLUX_CACHE keyed by (abs csv path, flux_type) caching cleaned arrays, prebuilt log-log interp1d objects, and per-band exponential fits (_build_flux_context, _interpolate_flux_from_context). Previously every simulate_lightcurve call re-read and re-sorted the CSV.
  • Mega-kernel results assembled straight from NumPy arrays into the DataFrame (no tolist()/zip() round-trips).
  • Theta-ring trig tables precomputed once per _simulate_phases_numba call instead of per phase.

Memory (mcmc_lightcurve_fit.py): grid precompute streamed worker results directly into flux_grid instead of materializing a full list(...), and interpolator setup avoided whole-grid copies during NaN cleanup. (Both became moot once the grid path was deleted in Phase 18.)

Hot path: likelihood invariants precomputed once in run_mcmc and passed as a like_terms dict (obs_err2, jitter_logf_index), replacing repeated active_names.index('log_f') lookups and per-call obs_err**2 allocations. _interp_periodic_phases gained a monotonic fast path that skips the sort.

Output: save_samples_csv_chunked writes samples in configurable chunks (--csv-chunk-size, default 50000); --compact-output adds an NPZ companion; --no-csv-output skips the large CSV entirely.

Model comparison: WAIC/LOO removed in favor of BIC. compute_bic_metrics computes BIC = k·ln n - 2 ln L̂ with k = len(active_names), n = len(obs_flux) after all binning/filtering, and L̂ obtained by calling the run's own likelihood at the max-log-prob sample (not by subtracting priors from log_prob, avoiding Jacobian bookkeeping). Reports bic, logL_hat, k_params, n_obs, theta_source (map_log_prob | median_fallback) to the console, mcmc_summary.txt, and *_model_metrics.csv; ΔBIC is computed against the best model in the run. ArviZ is now used for convergence diagnostics only.

Phase 16 — Frozen Parameters & Kepler Mode (ParamSpec)

Plan: freeze_params_and_kepler_3368d554.plan.md. Commit be56359.

Two features plus the refactor that made both tractable.

ParamSpec dataclass (built once in main() by build_param_spec(...) and threaded through every consumer) replaces ad-hoc positional theta indexing: mode, active_names, active_labels, frozen, fit_wind_shape, fit_scatter, wind_model, likelihood, orbital_period_s, K_kepler. Central resolvers by name, not index:

  • _resolve_geom(theta, spec) (replaces _to_physical, kept as a shim),
  • _resolve_shape(theta, spec) (generalizes _to_wind_params),
  • _theta_value(theta, name, names, frozen).

--freeze NAME=VAL[,NAME=VAL,…] — pins parameters and drops them from the chain. Valid: d1, d2, a, q, r, R, i0, M_X, M_RH, f_scatter, Rb, p, beta, fconf, ell. Shape params can be frozen even without --fit-wind-shape. log_f cannot be frozen. Unknown names are rejected with the allowed list; values outside the prior box warn but proceed; Rb < R with both frozen fails fast before sampling. get_active_priors drops frozen entries.

--kepler — samples (M_X, M_RH) in M☉ with KEPLER_PARAM_NAMES / KEPLER_PARAM_LABELS / KEPLER_PRIORS, then derives a = K·M_tot^{1/3} (with K = (G·M☉·P²/4π²)^{1/3}/R☉ precomputed by _compute_kepler_prefactor) and q = M_RH/M_tot from the lever arm d1·M_X = d2·M_RH. --orbital-period sets P (default ORBITAL_PERIOD); --prior-MX / --prior-MRH override the mass priors. Mutually exclusive with --reparam. Composes with freezing (e.g. --kepler --freeze M_RH=20).

Prior rewrite — log_prior iterates spec.active_names, applies box + Gaussian per dim, then enforces constraints on resolved values so they hold under freezing and Kepler mode: r < R (was a positional theta[2] >= theta[3]) and Rb ≥ R for smooth_pl. _log_jacobian contributes +log(a) only in reparam mode (Kepler priors are already in mass space).

Stats / persistence — compute_statistics derives (a, q, d1, d2) in Kepler mode and (d1, d2) in reparam mode, and records a MAP entry (max-log-prob single sample) which — unlike the marginal medians — exactly satisfies d1+d2 = a and d1/(d1+d2) = q. *_chain.npz gained mode, frozen_names, frozen_values, orbital_period_s; replot_from_existing reads them back to rebuild the spec, with old chains defaulting to previous behavior. plot_best_fit gained _value_from_stats_or_frozen so frozen parameters no longer break the overlay.

Also in this commit: xrb_lightcurve.py --Delta default changed from 2.0 to 1.0. Note this was not mirrored in default_wind_params or WIND_SHAPE_FIXED, which still use 2.0 — see PROJECT.md "Known rough edges".

Phase 17 — Per-Sample Phase-Shift Alignment

Commit 67448c3.

Rather than trusting the ephemeris to align model and data, every likelihood call now minimizes weighted χ² over a phase shift, so the fit is insensitive to residual epoch error.

  • _build_phase_shift_terms(enabled, obs_phase, grid_size, eval_points, refine_points) precomputes the coarse shift grid, the dense model evaluation grid, and the shifted observation-phase matrix once per run.
  • _apply_best_phase_shift(...) does a two-stage search: a coarse uniform grid over [0,1), then a local refinement across ±1 coarse step around the best shift — near-fine-grid accuracy at a fraction of the cost. The model is evaluated once on the dense grid and re-interpolated per trial shift.
  • Defaults: DEFAULT_PHASE_SHIFT_GRID_SIZE = 25, DEFAULT_PHASE_SHIFT_EVAL_POINTS = 240, DEFAULT_PHASE_SHIFT_REFINE_POINTS = 9. CLI: --no-fit-phase-shift, --phase-shift-grid-size, --phase-shift-eval-points.
  • Because the shift is a per-sample nuisance minimization rather than a sampled parameter, it is applied consistently in log_likelihood_chi2, log_likelihood_jitter, compute_chi2_for_samples, compute_pointwise_loglik, compute_bic_metrics, and plot_best_fit (which also draws the model shifted by the best-fit value and prints it in the annotation box).
  • mcmc_summary.txt gained a "Run configuration" block recording fit_phase_shift, grid size, eval points, and wall time (stats['_run_meta']), plus a chain-diagnostics block from print_diagnostics (stats['_diagnostics']).

Phase 18 — Adaptive Constant-SNR Binning & Grid Removal

Plan: adaptive_constant-snr_binning_249c8cba.plan.md. Commit b376585.

Constant-counts binning. counts is now carried through read_observation / load_data (counts_column='counts') and load_observed_lightcurves. New phase_bin_data_snr(df, counts_per_bin=100, …) in chandra_phase_analysis.py (with a flux-naming wrapper in the MCMC module) sorts by phase and greedily accumulates counts until each bin reaches the target, so every binned point carries roughly equal Poisson weight (100 counts ⇒ SNR ≈ 10) and low-signal eclipse troughs merge into wide bins instead of many noisy narrow ones. Per bin it returns the counts-weighted phase center, inverse-variance weighted flux, error = √(1/Σw), n_points, total_counts, and phase_lo/phase_hi/width; a trailing under-target bin is merged into its predecessor.

Mode is selected by argument presence, not a --bin-mode flag: --no-phase-bin (raw) > --counts-per-bin N (adaptive) > --n-phase-bins N (fixed) > neither (50 fixed bins, backward compatible). --n-phase-bins default changed from 50 to None; supplying both binners is an error. Bin widths are threaded as obs_phase_width into plot_best_fit and drawn as horizontal error bars.

Precomputed grid deleted. PrecomputedModelGrid and its _precompute_models / _compute_single_model / _setup_interpolators / _save_grid / _load_grid machinery are gone, along with --save-grid, --load-grid, --no-grid, --grid-points, and --shape-grid-points. MCMC always uses DirectLightCurveModel: at ~60 ms/LC it is fast enough, it avoids grid interpolation artifacts, and it is the only path that keeps the likelihood physically faithful for every sample (including per-step wind shape).

Likelihood CLI simplified. Student-t and the deprecated WAIC shim removed: log_likelihood_studentt(), 'studentt' from LIKELIHOOD_TYPES and --likelihood choices, --studentt-nu, all studentt_nu plumbing through log_probability / run_mcmc / compute_chi2_for_samples / BIC / run_single_fit / save / replot, the scipy.special.gammaln import, and --compute-waic plus its shim mapping to --compute-bic. Only chi2 and jitter remain; --compute-bic is the model-comparison flag.

Phase 19 — Gaussian Phase Smoothing, Scattered Flux & Residual Panels

Plan: gaussian_phase_smoothing_reference_d1a46172.plan.md. Status: implemented but uncommitted on branch add_generic_wind.

Three reusable primitives live in chandra_phase_analysis.py so both the single-model and MCMC plot paths share one source of truth:

  1. smooth_lightcurve(phase, flux, flux_err, sigma=0.01, n_eval=300, n_mc=2000, random_state=None) — periodic Gaussian-kernel phase smoother. Periodic distance d = |((φ_i-φ_eval+0.5) mod 1) - 0.5|, weights exp(-½(d/σ)²), so it is continuous across phase = 0/1. Generalizes the MATLAB reference in temp/LC_MC/*.m from index windows to a phase-distance kernel, so it is correct for fixed-width bins, constant-SNR bins, and raw unbinned data. The 1σ band is a vectorized Monte Carlo (perturb all points at once, one matmul, std over realizations) rather than a Python loop. Kernel is phase-distance only — no inverse-variance weighting — matching the reference. σ = 0.01 sits well below the ~0.1–0.25 phase scale of real features and above the ~0.0002 raw sampling.
  2. estimate_scattered_flux(phase, flux, window=(0.4, 0.6)) — mean observed flux in the mid-eclipse window (fallback 0.1 × median, clamped ≥ 0).
  3. add_residual_panel(ax, phase, obs, model, err, xerr=None) — normalized pulls (O-M)/σ with 0/±1 reference lines.

xrb_lightcurve.py — simulate_lightcurve(scattered_flux=0.0) adds a constant, phase-independent offset to every nfl_* column after eclipse handling, so notebooks can bake a scattered-light floor into a directly generated model. The fit paths deliberately add scatter at overlay/evaluation time instead, so fit_simulation's multiplicative scale and MCMC's per-step scaling don't rescale an additive constant.

chandra_phase_analysis.py single-model path — fit_simulation(scatter=…) adds the constant after scaling inside the inner χ²; plot_phase becomes a 2-panel figure (3:1 heights, shared x) with a residual panel whenever a model overlay is present, plus the optional dashed-green smoothed curve and MC band. New CLI: --smooth, --smooth-sigma, --smooth-n-mc, --smooth-seed, --scatter, --scatter-eclipse-phase. When --fit runs without --scatter, the value is estimated from the eclipse window.

mcmc_lightcurve_fit.py — imports the three primitives. Same --smooth* flags, computed once per band and threaded into plot_best_fit, which is now a 2-panel figure with the residual panel clipped to ±5σ. New --fit-scatter promotes f_scatter to a free MCMC parameter (mirroring log_f): added by build_param_spec / ParamSpec.fit_scatter, resolved by _resolve_scatter and applied in _evaluate_model so all likelihoods, BIC, and pointwise log-lik inherit it. Its prior is centered on estimate_scattered_flux(...) with min = 0 and max = nanmax(obs_flux). Being phase-invariant it is unaffected by the phase-shift search; being physical it counts toward dof automatically (n_phys subtracts only log_f). It is saved as a normal active column, excluded from wind-shape extra-dim detection on replot, and can be pinned via --fit-scatter --freeze f_scatter=<v>.

Phase 20 — Remove Multiplicative Flux Scale from the Single-Model χ² Fit

chandra_phase_analysis.py. Status: uncommitted on branch add_generic_wind.

fit_simulation used to fit two parameters — a phase shift and a multiplicative flux scale — which meant reduced χ² never tested the model's absolute normalization. That normalization is not free: it is pinned by lam (the orbit-averaged nH from the spectral fit) together with the XSPEC flux vs nH table. A free y-scale therefore absorbed any normalization error instead of exposing it. The MCMC path never had such a parameter — it fits a per-sample phase shift and an additive f_scatter, nothing multiplicative — so the two χ² paths were also inconsistent with each other.

Note the additive scattered-flux floor is not a substitute for the multiplicative scale (they are different degrees of freedom); the justification for dropping scale is that the normalization is externally fixed, and the eclipse-floor offset is what the additive term legitimately covers.

  • fit_simulation(obs_df, sim_df, sim_column, fit_phase_shift=False, scatter=0.0, n_shift_grid=1000) now returns (shift, reduced_chi2) — the scale element is gone from both the fit and the return tuple. The model is interp((φ_obs - shift) mod 1) + scatter.
  • Robust shift search. χ²(shift) is periodic and strongly multi-modal because of the eclipse, and the old Nelder-Mead started at shift = 0 routinely settled in the wrong basin — a latent defect that mattered more once the shift became the only fitted parameter. Replaced with a vectorized coarse scan over the full period (default 1000 nodes, one np.interp over all (shift, φ_obs) pairs) followed by a bounded minimize_scalar refinement within one coarse step. Verified to recover injected shifts of 0.02, 0.37, 0.51, and 0.95 to < 1e-3.
  • dof corrected. Was N - 2 in both branches (wrong even before: the no-rescale branch fitted nothing). Now N - 1 when the shift is fitted and N when it is not.
  • plot_phase lost its scale parameter; the overlay is drawn at native normalization plus scatter. rescaled= renamed to shift_fitted=, and the title annotation now reports the phase shift instead of a rescaled/not label.
  • plot_multi_column_fits — fit_results entries are (shift, chi2); rescaled= → shift_fitted=.
  • CLI — --rescale renamed to --fit-phase-shift, with --rescale kept as a deprecated alias on the same dest so existing commands keep working.
  • Unused import removed — scipy.optimize.minimize → minimize_scalar.
  • Notebook call sites migrated (10 cells across notebooks/xrb_toy_wind_models.ipynb, xrb_model_analysis.ipynb, xrb_model_analysis_single_15803.ipynb); .bak copies left alongside.

Effect on reported fit quality. On the soft band (12 obs, 100 fixed bins, sim_flux_tbabs_15803_smooth_pl_soft.csv) the old two-parameter fit returned scale = 0.413 with χ²/dof = 2.06; with the scale removed the same data/model give χ²/dof = 20.97. The free scale had been hiding a ≈ 59 % flux-normalization deficit. Expect previously "acceptable" reduced χ² values to rise across the board — that is the intended behavior, and a coherent non-zero-centered residual band is now the diagnostic signature of a normalization mismatch (as opposed to a shape mismatch).

Not changed: chandra_analysis_combined_flux.py carries its own older duplicate of fit_simulation / plot_phase / plot_multi_column_fits (no scatter support at all) and still fits a multiplicative scale.

Follow-up — model/χ²/residual consistency

Two ways the plot could disagree with the number printed on it, both caused by "model evaluated at the observed phases" being written out three separate times with three subtly different formulas:

  1. Overlay could omit the scatter floor. plot_phase takes its own scatter argument, so a caller that passed scatter= to fit_simulation but not to plot_phase got a curve drawn scatter too low while the title showed the χ² that included it. main() always threaded both, but the notebook call sites did not.
  2. The residual panel ignored the phase shift entirely. It built its interpolation grid from the unshifted model (model_interp_phase = phase_sorted) and evaluated np.interp(obs_phase, …), whereas χ² uses interp((obs_phase - shift) mod 1, …). The plotted line was shifted correctly, so only the pulls were wrong — and badly: on a synthetic shift = 0.30 case with an otherwise perfect fit (χ²/dof = 0), the panel showed residuals implying χ²/dof ≈ 18050, with 50/120 points disagreeing with the χ² model.

Fixed by collapsing all three call sites onto one definition:

  • _prepare_model_interpolator(sim_df, sim_column) — builds the wrap-around (phase_wrap, flux_wrap) arrays, accepting either a phase or a deg column, and no longer mutates the caller's sim_df.
  • _model_from_wrap(phase_wrap, flux_wrap, phases, shift, scatter) — the single definition of the model. Accepts an array-valued shift so the coarse-scan batching in fit_simulation uses the same expression as a single evaluation.
  • evaluate_model_at_phases(sim_df, sim_column, phases, shift, scatter) — public one-shot wrapper (useful from notebooks).
  • _obs_errors(obs_df) — shared uncertainty extraction (provided errors else sqrt(|rate|), with zero/negative/non-finite floored to 1e-3), so the fit and the residual panel weight points identically. Non-finite errors are now floored too; previously a NaN error propagated into a NaN χ².
  • plot_phase draws the overlay on a dense 721-point grid via the shared evaluator instead of shifting the model's own sample points, and computes residuals with the same shift and scatter.

Added a self-check. When plot_phase is given a chi2 to display, it recomputes reduced χ² from the curve it actually drew (matching dof via shift_fitted) and warns if the two differ by more than 1 %. This catches a mismatched scatter, shift, or sim_column at the point of display rather than leaving a plausible-looking number over the wrong curve — it reproduces both bugs above as warnings. Verified silent on correct calls.

Notebook fix. One genuine mismatch existed — notebooks/xrb_toy_wind_models.ipynb cell 8 fitted with scatter=1.92174e-13 but called plot_phase without it; now passes it. An audit of all fit_simulation → plot_phase/plot_multi_column_fits pairs across the four notebooks found no others.

Phase 21 — MCMC Scatter-Path Audit

mcmc_lightcurve_fit.py. Status: uncommitted on branch add_generic_wind. All tests run under the henv conda env (Python 3.13 / numpy 2.2 / emcee 3.1.6 / arviz 1.0).

Audit of whether the Phase 20 fixes have analogues in the MCMC path. The MCMC never had a multiplicative flux scale, so nothing to remove there, and its per-sample phase-shift search (Phase 17) was already the model for the chandra_phase_analysis one. plot_best_fit was already fully consistent — it applies both the phase shift and f_scatter to the overlay curve, the shift-search model, and the residual basis (verified to agree with the likelihood to 1.7e-16 relative). Three genuine problems in the --fit-scatter path were found and fixed.

1. compute_chi2_for_samples silently dropped f_scatter. It duplicated the geometry/shape resolution and called model.evaluate directly, never invoking _resolve_scatter. So every per-sample χ² written to *_chi2.csv.gz was computed against a model missing the additive floor whenever --fit-scatter was active. Measured bias on a 40-bin broad-band fit:

f_scatter correct χ² χ² with floor dropped error
0 6228.11 6228.11 0.0 %
1e-13 7697.40 6228.11 −19.1 %
3e-13 11120.15 6228.11 −44.0 %
1e-12 28183.54 6228.11 −77.9 %

Fixed by routing it through _evaluate_model (the same entry point the likelihood uses) instead of re-implementing the evaluation, which also picks up wind-shape resolution and any future additions automatically. The now-dead local frozen binding was removed. Verified: all five reporting paths — log_likelihood_chi2, compute_chi2_for_samples, plot_best_fit, compute_bic_metrics, compute_pointwise_loglik — agree to machine precision with f_scatter active.

2. --fit-scatter could not start at all under emcee. Walker initialization clipped each dimension to min + 0.01*|min| + 1e-12 … max - 0.01*|max| - 1e-12. That absolute 1e-12 epsilon is meaningless for a parameter whose natural scale is ~1e-13: with f_scatter prior min = 0, the lower bound became exactly 1e-12, roughly 6× larger than the entire plausible range, so every walker was clipped to the same value. The column had zero variance, emcee.walkers_independent returned False, and the run died with ValueError: Initial state has a large condition number. The inset is now relative to each parameter's own prior span (pad = 1e-9 * (max - min)), with a fallback to the raw box if the inset inverts. Regression-checked: every existing parameter's clip bounds move by ~1e-9 relative (numerically identical); only f_scatter changes. Added a final guard that re-spreads any dimension that still collapses to a constant, with a warning naming the parameter, so a badly scaled prior degrades loudly instead of aborting inside emcee.

3. mcmc_summary.txt reported f_scatter: 0.000000. The summary and console writers used fixed-point %.6f, which rounds a ~1e-14 flux floor to zero — it read as "not fitted" even though the posterior was well constrained. Added _fmt_val(value, width=0), which switches to %.6e for non-zero magnitudes below 1e-4, and applied it to the marginal-posterior, derived-parameter and MAP blocks plus the print_results console table.

Verified end-to-end under henv: --fit-scatter with both chi2 and jitter likelihoods, --save-chi2 --compute-bic --smooth, produces all artifacts; f_scatter is genuinely sampled (230 unique values spanning 8.2e-17 to 1.7e-13, no longer pinned); the χ² table is 100 % finite; and plot_best_fit's reduced χ² at the MAP matches the χ² table row for that same sample (239.2630 vs 239.2631).

4. plot_best_fit duplicated the evaluation logic. It called eval_fn = getattr(model, 'evaluate_direct', model.evaluate) in three places and added f_scatter by hand to each, and it re-implemented geometry resolution (its own if reparam: stats['d1']… branch). Numerically correct, but it was the same duplication pattern that caused problem 1, and the evaluate_direct fallback was dead — it referenced PrecomputedModelGrid, removed in Phase 18. Now it reconstructs the point-estimate theta in active-parameter order and calls _evaluate_model through a local _eval_at(phases) helper, so geometry (phys / reparam / kepler / frozen), wind shape, and f_scatter all resolve through one implementation. f_scatter_best is obtained from _resolve_scatter and used for display only. Verified behaviour-preserving: plot_best_fit's χ² matches the likelihood exactly across eight configurations — phys, phys + f_scatter, reparam, kepler, jitter + f_scatter, smooth_pl shape fit, beta_law shape fit (with R_star tied to R), and frozen R + f_scatter.

Phase 22 — Adaptive Binning in the Single-Model CLI

chandra_phase_analysis.py. Status: uncommitted on branch add_generic_wind.

Phase 18 added phase_bin_data_snr to chandra_phase_analysis.py but wired the CLI flag only into mcmc_lightcurve_fit.py, so the single-model script could not use adaptive constant-counts bins from the command line. Added --counts-per-bin with the same argument-presence semantics as the MCMC script:

--no-phase-bin > --counts-per-bin N > --n-phase-bins N > 50 fixed-width bins.

  • --n-phase-bins default changed from 50 to None so mode selection is unambiguous; the effective fallback is still 50, so existing commands behave identically.
  • Supplying both binners is a parser.error, as are non-positive values.
  • --counts-per-bin on data with no counts column fails fast with a message pointing at --n-phase-bins, rather than raising from inside the binner.
  • --min-points-per-bin help now notes it applies to fixed-width binning only.
  • Variable bin widths already flowed through to horizontal error bars, since plot_phase picks up the width column when is_binned — no plotting change needed.

Verified on ObsID 15803 soft: --counts-per-bin 100 gives 87 bins averaging 104.3 counts (vs 60 fixed-width bins), and χ²/dof drops from 8.501 to 6.385 as the noisy narrow eclipse-trough bins merge.

Phase 23 — utils/ Extraction and a Single Plotting Routine

chandra_phase_analysis.py, mcmc_lightcurve_fit.py, new utils/utils.py, utils/plot_utils.py, utils/__init__.py. Status: uncommitted on branch add_generic_wind.

Both analysis scripts carried their own plotting functions, and mcmc_lightcurve_fit.py imported its data-layer helpers from chandra_phase_analysis.py — so the CLI script was simultaneously a library and a front end, and mcmc depended on it for reasons unrelated to Chandra data. Everything shared now lives in a utils package and neither script imports the other.

New layout

Module Contents
utils/utils.py (1071 lines) REF_EPOCH, ORBITAL_PERIOD, frac, fmt_val, band_label_from_column, detect_flux_columns, validate_sim_columns, read_observation, load_data, phase_bin_data, phase_bin_data_snr, smooth_lightcurve, estimate_scattered_flux, prepare_model_interpolator, model_from_wrap, evaluate_model_at_phases, interp_periodic_phases, obs_errors, fit_simulation. numpy/pandas/scipy only.
utils/plot_utils.py (597 lines) plot_lightcurve_fit (the one drawing routine), plot_phase, plot_multi_column_fits, plot_corner, plot_trace, add_residual_panel, build_fit_title, format_reduced_chi2, half_widths.

chandra_phase_analysis.py drops from 1639 to 457 lines and is now only the argparse CLI plus an __all__ re-export block, so from chandra_phase_analysis import * — the notebooks' import style — is unchanged. mcmc_lightcurve_fit.py drops from 4124 to 3952 lines.

One plotting function

plot_lightcurve_fit is the drawing code that used to be inlined in mcmc_lightcurve_fit.plot_best_fit, generalized so both paths reach it. It draws only what it is handed — observed arrays, an already-shifted overlay curve, and the model evaluated at the observed phases — which is what makes it usable from both:

  • plot_best_fit still owns the MCMC-specific work (MAP-vs-median point estimate, _evaluate_model, _apply_best_phase_shift) and then delegates.
  • plot_phase is now a thin adapter: interpolate sim_df at the given shift and additive scatter, then delegate. Its displayed-χ² self-check warning (Phase 20) moved with it.

Feature parity required generalizing three things: an optional obs_group array (so Chandra's per-observation series still get their own colors and legend entries), an ax/ax_res pair (so plot_multi_column_fits can still draw into a grid cell, residual-panel-free), and a title escape hatch for the no-model data plot. Two incidental fixes fell out: a model overlay with no error column no longer produces an empty residual panel, and --counts-per-bin widths reach the x-error bars through one code path instead of two.

Removed from the figure

Per request, plot_best_fit's wheat-colored annotation box is gone. The title is now only the energy band and χ²/dof:

SOFT band  —  χ²/dof = 13.018

The information it carried is not lost — point estimate, phase_shift, f, chi2_eff/dof and f_scatter are printed to stdout next to the existing parameter table, which is also written to the run summary. plot_phase titles gained the band label the same way: nfl_soft -> "SOFT band" via band_label_from_column, so grid panels remain identifiable.

Deduplicated

  • _interp_periodic_phases (mcmc) and _model_from_wrap (chandra) were two spellings of periodic model interpolation; both now live in utils.utils as interp_periodic_phases / model_from_wrap, documented as the array-in and prepared-interpolator forms of the same operation.
  • plot_corner, plot_trace, add_residual_panel, fmt_val moved out of mcmc_lightcurve_fit.py verbatim.
  • mcmc_lightcurve_fit.py no longer imports matplotlib.pyplot or corner at all — all figure work is behind utils.plot_utils.

Verification (conda env henv, Python 3.13.5 / numpy 2.2.6)

  • py_compile on all five files; import of both scripts; star-import exports 23 names with nothing missing.
  • Chandra CLI, four modes on ObsID soft data: adaptive bins + fitted shift + smoothing (χ²/dof 30.605, no self-check warning), 3-column grid (98.007 / 15.071 / 46.100, one title per band), fixed-width no-fit plot, raw unbinned fit.
  • MCMC: 24 walkers × 60 steps with --fit-scatter --smooth produced χ²/dof 13.0179; --replot from the saved chain reproduced 13.0179 exactly.
  • Consistency suite, 15 assertions: across phys / phys, shift fixed / phys + f_scatter, the χ²/dof plot_best_fit reports equals the χ² of the arrays it hands the plotter to rtol=1e-12, and the drawn overlay reinterpolated onto the observed phases matches the residual basis to ≤5.0e-4 relative (the Phase 20 regression would blow this up by orders of magnitude). fit_simulation recovers an injected shift=0.30 to 5 decimals with χ²/dof = 0; the self-check guard stays silent when shift/scatter agree and fires when scatter is wrong.
  • All five call forms the notebooks use for plot_phase / plot_multi_column_fits run against real data with no self-check warnings. (End-to-end nbconvert execution still stops at the notebooks' import xspec cell — a pre-existing henv limitation, unrelated to this change.)

Not migrated: chandra_analysis_combined_flux.py keeps its own older copies of fit_simulation / plot_phase / plot_multi_column_fits and still fits a multiplicative flux scale. See Known rough edges.

Phase 24 — Run-Config Persistence and a Replot-Mode Fix

mcmc_lightcurve_fit.py. Status: uncommitted on branch add_generic_wind.

Bug: --replot could not read a kepler- or reparam-mode chain

load_existing_results validated the samples CSV against REPARAM_PARAM_NAMES if reparam else PARAM_NAMES — the kepler argument it accepted was never used. Replotting a --kepler run therefore failed with

Error: Samples file missing required geometry columns: ['d1', 'd2']

even though the file correctly held M_X, M_RH, r, R, i0, .... Fixed to select the geometry block by mode, and two related gaps closed:

  • Frozen parameters are not sampled and so are legitimately absent from the CSV; they are now excluded from the required set (--freeze R=2.0 runs were unreplottable for the same reason).
  • The error now lists the columns actually present and, when they match another parameterization, names the flag to use.

Feature: the CLI of every fit is saved and restored

--replot previously recovered only mode, frozen values, orbital_period_s and likelihood from the chain .npz. Everything else fell back to argparse defaults — including --data-dir, --obs-column, --time-column, --counts-per-bin/--n-phase-bins, --lam, --dth, --d2h and all --prior-*. Those options determine the observed arrays, so a replot that missed them reported a χ²/dof for a dataset the posterior had never seen. On the existing broad kepler run, replotting with default --lam/--flux-csv gave χ²/dof = 0.916 against the true 1.064.

Every fit now writes <band>_<wind_model>_run_config.json into --output-dir:

{
  "created": "2026-08-26T16:53:02-0400",
  "command": "mcmc_lightcurve_fit.py --band broad --kepler ...",
  "band": "broad",
  "wind_model": "smooth_pl",
  "args": { "...": "every argparse dest" }
}

It is written before sampling starts, so it survives an interrupted run. On --replot, apply_saved_run_config fills in every option the user did not type:

  • Explicit flags always win. Which options were typed is determined by scanning sys.argv against the parser's option strings (including unambiguous prefixes), not by comparing against defaults — a user who explicitly passes the default value still overrides the saved config.
  • replot and output_dir are never restored: the first would cancel the replot (a saved fit always recorded replot=False), and the second is defined by where the config was found.
  • --band and --flux-csv changed from required=True to being validated after the restore, so --replot alone is a complete command. They are still required for a real fit.
  • Ambiguity is handled: with several configs in one directory, ones differing only by band (the --band all case) restore identically and the first is used; genuinely different configs produce an error listing them.
  • Self-healing: a --replot that finds no config writes one from the options it just used, so pre-existing result directories become self-sufficient after one full-CLI replot.

As an independent backstop, replot_from_existing now compares the observed point count against n_obs in the chain metadata and warns on a mismatch. n_obs is also now stored unconditionally rather than only when --compute-bic ran, so the check works for every run.

Verification (conda env henv)

  • The existing mcmc_results/broad_smooth_pl kepler + wind-shape + f_scatter run replots without error. With the original CLI it reproduces the saved numbers exactly: χ²/dof 1.06423 (summary recorded 1.064) and BIC 195.673 (logL_hat=-77.689, k=8, n=154).
  • A second --replot --output-dir mcmc_results with no other arguments reproduces the same 1.06423 / 195.673 from the auto-written config.
  • Fresh 24×40-step soft fit: config written, then bare --replot round-tripped χ²/dof = 13.2843 identically.
  • --replot --lam 0.9 reports kept from the command line: --lam, confirming precedence.
  • Mismatch guard fires as intended: forcing --counts-per-bin 300 on the broad run warns "Replot is using 53 observed points but the saved fit used 154".
  • --band/--flux-csv still error out when missing on a non-replot run.

Phase 25 — MCMC Script Slimming

mcmc_lightcurve_fit.py, utils/utils.py. Status: uncommitted on branch add_generic_wind. Behaviour-preserving throughout — every regression number below is bit-identical to the pre-cleanup run.

mcmc_lightcurve_fit.py: 4246 -> 3553 lines (-693, -16%).

Moved into utils/utils.py (generic, not MCMC-specific)

Moved Lines
resolve_band_directory, load_observed_lightcurves 90
build_phase_shift_terms, apply_best_phase_shift + the DEFAULT_PHASE_SHIFT_* constants 95
run-config persistence: save_run_config, find_run_configs, apply_saved_run_config, _explicit_cli_dests, _jsonable, run_config_path 165
save_samples_csv_chunked 23

The phase-shift search now sits next to fit_simulation, which runs the same coarse-scan-then-refine algorithm on a tabulated model — both spellings of the idea are in one place. utils/utils.py gained a stdlib-only dependency footprint (argparse, csv, json, shlex, time) and stays free of any import from either analysis script.

Deleted by making the shared binners honest

phase_bin_data / phase_bin_data_snr hardcoded rate/error as their output column names, so the MCMC script carried two 25-50 line wrappers that renamed flux/flux_err in and back out again. The binners now name their value columns after the rate_column/error_column arguments they were given, so the MCMC path calls them directly:

bin_cols = dict(rate_column='flux', error_column='flux_err')
obs_df = phase_bin_data_snr(obs_df, counts_per_bin=..., **bin_cols)

chandra_phase_analysis passes rate/error and is unaffected. -73 lines.

Deduplicated

  • Likelihood front half. log_likelihood_chi2 and log_likelihood_jitter shared ~23 lines of identical model-evaluation + phase-shift-alignment preamble. Extracted as _aligned_model_flux; each likelihood is now its own formula plus a call. f = exp(log_f) moved after the alignment (it does not depend on the model, so the result is unchanged).
  • CLI plumbing. _phase_shift_opts(args) replaces 8 copies of the same three-line fit_phase_shift= / phase_shift_grid_size= / phase_shift_eval_points= getattr triple. _smooth_plot_kwargs(smoothed, args) replaces two 13-line blocks of smoothed[...] if smoothed is not None else None.
  • --prior-* definitions. Nine near-identical 8-line add_argument calls became a table-driven loop reading defaults straight from DEFAULT_PRIORS / REPARAM_PRIORS / KEPLER_PRIORS, so the help text cannot drift from the values actually used. It already had: --prior-r advertised max=0.01 while the code used 0.1. -50 lines, one stale-help bug fixed.
  • Prior parsing. The geometry and wind-shape override loops were the same 20 lines twice; now one _parse_prior_overrides(parser, args, names).
  • replot_from_existing opened *_chain.npz twice, in two separate try/except blocks, to read mode/frozen/n_obs and then likelihood. Merged into one read.
  • Label construction. The 17-line "name -> corner label" if/elif chain in the replot path became _labels_for_names(names, mode), built on the existing get_mode_name_label.
  • _default_priors(reparam, kepler) replaces three copies of the if kepler / elif reparam / else prior-selection block (and upgrades two of them from .copy() to deepcopy, so a prior override can no longer reach the module-level dict).
  • _maybe_save_chi2(...) replaces two 15-line --save-chi2 blocks.

Dead code removed

  • compute_pointwise_loglik (84 lines) — the per-observation log-likelihood matrix that WAIC/LOO consumed. Traced before removing: it was added with its caller in b867ff1 and called for six commits from inside run_arviz_diagnostics, gated on --compute-waic, feeding log_lik_dict = {"obs": ll[np.newaxis, :, :]} into _build_inference_data and thence az.waic / az.loo. Phase 15 (7116302) deleted that call when BIC replaced WAIC/LOO — L̂ from the run's own likelihood at the max-log-prob sample needs no pointwise terms — leaving --compute-waic as a shim; Phase 18 (b376585) deleted the flag and shim outright. Unreachable ever since, with the only remaining trace being _build_inference_data's optional log_lik_dict parameter, whose sole caller passes None. Note it was maintained dead code (it picked up the Phase 17 phase-shift alignment and lost its studentt branch with Phase 18), so it is a working starting point if LOO is ever wanted back: git show 19dc637:mcmc_lightcurve_fit.py.
  • _to_physical — a "backward-compatible wrapper for legacy callers" with no callers anywhere in the repo.
  • A dead import gzip inside compute_chi2_for_samples (the gzip write goes through to_csv(compression='gzip')), and now-unused module imports (csv, json, shlex, sys, glob, pickle).

Comments trimmed

Thirteen multi-line rationale comments (7-10 lines each) cut to 1-3 lines, keeping the why and dropping the retold debugging narrative — e.g. the walker-clip block went from eight lines to three while still naming the failure it prevents (f_scatter collapsing and aborting emcee on the condition number). -37 lines.

Verification (conda env henv) — all numbers unchanged

  • Bare --replot --output-dir mcmc_results on the real broad kepler + wind-shape + f_scatter run: χ²/dof 1.06423, BIC 195.673 (logL_hat=-77.689, k=8, n=154) — identical to Phase 24.
  • Consistency suite: 18/18 pass, including three plot_best_fit configurations agreeing with the drawn arrays to rtol=1e-12 at the same values as before (1.117909384 / 1.118796619 / 1.392422282), and new checks that the moved functions resolve to utils.utils and the binners preserve column names.
  • Fresh fits round-trip through --replot exactly in every mode: chi2 + f_scatter (12.9794), jitter + f_scatter + --save-chi2 + --compute-bic (16.698, BIC -2011.323), --reparam (5.71919), --kepler --fit-wind-shape --freeze R=13.0 (9.82407).
  • Chandra CLI, all four modes: 30.605 / (98.007, 15.071, 46.100) / 40 bins / 3.129 — unchanged.
  • All five notebook plot_phase / plot_multi_column_fits call forms run with no self-check warnings. ArviZ, corner and trace paths exercised.
  • No unused imports remain in any of the four files.

Side Investigation — Reference Epoch Recalibration

Plan: reference_epoch_recalibration_ae1cf98a.plan.md.

Per Laycock et al. 2015 (stu2151.pdf, §4), T0 = 278801348 s = MJD 54040.87 is the mid-eclipse time of ObsID 07082, defined to lie at phase 0.5, with a full eclipse width ~0.2 in phase (~7 h); the paper folds with φ = (t - T0 - 100000P)/P. The codebase formula frac((t - T0)/P) instead puts that reference mid-eclipse at phase 0.0, and the precomputed phase column stored in the CIAO files reproduces neither exactly.

A read-only study script (find_reference_epoch.py) located mid-eclipse two ways — anchored on ObsID 15803 (the only observation with a full clean eclipse, ~1.73 d > one 1.45 d orbit) via sliding-window / threshold / trapezoid-fit, and via a joint all-observation scan over trial phase offsets maximizing eclipse contrast — with the period held fixed at 125431 s.

Outcome: a corrected epoch of 278800407.267 s, which sits commented out next to REF_EPOCH in chandra_phase_analysis.py; the active value is still 278801348. In practice the Phase 17 per-sample phase-shift search absorbs the offset, so this now mainly affects the interpretability of plotted phases. The study script itself is not present in the working tree.

Current File Inventory

Core simulation / inference

File Lines Status Description
xrb_lightcurve.py 2115 active Forward model: profiles, Numba GL LOS kernels, simulate_lightcurve, physical back-calculation helpers.
mcmc_lightcurve_fit.py 3553 active emcee/zeus MCMC: ParamSpec (phys/reparam/kepler), freeze, wind-shape, f_scatter, phase-shift search, BIC, ArviZ, replot. Direct evaluator only.
chandra_phase_analysis.py 457 active CLI front end for the single-model χ² fit; re-exports the shared utils/ API for the notebooks.
utils/utils.py 1474 active Shared layer: ephemeris, loading, both binners, smoothing, periodic model interpolation + phase-shift search, fit_simulation, run-config persistence.
utils/plot_utils.py 597 active All plotting, built on the single plot_lightcurve_fit drawing routine.
compute_flux_vs_nH.py 934 active XSPEC table generator (flux vs nH).
xspec_fit_mcmc.py 702 active XSPEC-side MCMC for spectral fits.
chandra_analysis_combined_flux.py 539 active Pre-folded combined-flux phase analysis.
plot_results.py 272 active Standalone plotting of simulation CSVs.
compute_count_to_flux_factor.py 147 active Count-rate → flux conversion factor.
example_usage.py 96 active Programmatic simulate_lightcurve examples.

XSPEC / shell helpers

compare_models.sh, convert_fits_to_txt_heasoft.sh, xspec_tbabs_fit_results.xcm, rkp_run_w_mcmc_cmds.sh (command scrapbook — contains dead flags from earlier phases).

Missing from the working tree despite being referenced: compare_absorption_models.xcm (invoked by compare_models.sh), xspec_get_conversion_factors_tbabs.xcm, get_xspec_nH.py, utils/get_conversion_factors.sh, find_reference_epoch.py.

Utilities (utils/)

Now a package (__init__.py). utils.py and plot_utils.py are library code imported by both analysis scripts (see Core above). Standalone data-prep scripts, not part of the package API: add_flux_simple.py, add_flux_to_lightcurves.py, convert_fits_to_txt.py, get_average_count_rates.py, test_flux_methods.py, benchmark_mcmc_performance.py.

Notebooks

xrb_toy_wind_models.ipynb (active wind-profile exploration), xrb_model_analysis.ipynb, xrb_model_analysis_single_15803.ipynb, xrb_flux_nH_abs.ipynb.

Documentation

PROJECT.md (current-state reference), changes_tracked.md (this file), mcmc_chi2_jitter_explanation.md, PERFORMANCE_VALIDATION_REPORT.md, FLUX_INTEGRATION_SUMMARY.md, FLUX_METHODS_QUICKREF.md, XSPEC_CONVERSION_GUIDE.md, FITS_CONVERSION_README.md, QUICK_START_FLUX_CONVERSION.md, CONVERSION_WORKFLOW.md, README_CONVERSION_TOOLS.md. Reference PDFs: Wind_Density.pdf (profile equations), stu2151.pdf (Laycock et al. 2015 ephemeris). Stale: README.md, MIGRATION_SUMMARY.md — still describe the original R→Python port and removed API (--lam2, flx2/fl2, nfl_*_av/_cv, pho_count_*).

Plans (.cursor/plans/)

unified_wind_model_77726ced (Phase 7), mcmc_performance_and_statistics_8989cd39 (Phase 8), mcmc_convergence_improvements_c648afb3 (Phase 9), mcmc_wind_shape_params_8b9c89d2 (Phase 10), flux_t_error_and_unbinned_mcmc_1bf81aa9 (Phase 12), wind_normalization_constants_76f1ef51 (Phase 13), mcmc_speed_memory_optimization_68fc497a (Phase 15), freeze_params_and_kepler_3368d554 (Phase 16), adaptive_constant-snr_binning_249c8cba (Phase 18), gaussian_phase_smoothing_reference_d1a46172 (Phase 19), reference_epoch_recalibration_ae1cf98a (side investigation).

Legacy (preserved)

legacy_r_code/ (new11.R, grid4.R, wind_los2.R, density_fnc.R), light_curve_model_opt_bw.R.

Current Status & Quick Commands

See PROJECT.md for the full current-state reference (module-level API, data layout, outputs, and known rough edges). Summary:

Environment: henv conda env (heasoft/XSPEC + python deps including numba, emcee, zeus-mcmc, arviz, corner, astropy, scipy>=1.12). requirements.txt is incomplete — it predates the numba/arviz/zeus additions.

Working spectral model: TBabs × powerlaw, nH = 0.75×10²² cm⁻², Γ = 1.86, χ²_red = 1.52.

Default forward model: smooth_pl wind, single nfl_{band} flux column, Gauss-Legendre mega-kernel (~60 ms per light curve). Note Rb=5, p=4 throughout, but Delta defaults to 1.0 on the xrb_lightcurve.py CLI and 2.0 in default_wind_params / WIND_SHAPE_FIXED.

MCMC defaults: phys mode (d1, d2, r, R, i0), chi2 likelihood, emcee, 50 fixed-width phase bins, per-sample phase-shift alignment on, direct evaluator (no grid path exists).

# Build/refresh the XSPEC flux-vs-nH table
python compute_flux_vs_nH.py --specdir ./data/IC10X1_spec --model tbabs \
    --bands broad soft medium hard --out_csv flux_vs_nH_tbabs_broad.csv

# FITS pipeline
./convert_fits_to_txt_heasoft.sh
python utils/add_flux_simple.py \
    data/IC_10_X1_LC/Broad_converted/ data/IC_10_X1_LC/Broad_with_flux/ 1.500509e-11

# Simulate one light curve
python xrb_lightcurve.py --flux_method interpolate \
    --flux_csv flux_vs_nH_tbabs_broad.csv \
    --wind-model smooth_pl --Rb 5 --p 4 --Delta 1 \
    --i0 12.0 --lam 0.572385 --output sim_broad.csv

# Single-model χ² fit + smoothed overlay + residual panel
python chandra_phase_analysis.py --data-dir data/IC_10_X1_LC_CIAO/broad \
    --obs-column flux_t --time-column t_raw \
    --fit --sim-file sim_broad.csv --rescale --smooth \
    --n-phase-bins 100 --output fit_broad.png

# MCMC: geometry only, reparameterized, adaptive constant-SNR bins
python mcmc_lightcurve_fit.py --band broad \
    --flux-csv flux_vs_nH_tbabs_broad.csv --data-dir data/IC_10_X1_LC_CIAO \
    --obs-column flux_t --time-column t_raw \
    --wind-model smooth_pl --reparam --likelihood chi2 \
    --counts-per-bin 100 --sampler zeus --dth 4.0 \
    --n-walkers 24 --n-steps 20000 --n-burn 2000 \
    --compute-bic --smooth --output-dir mcmc_results/broad/smooth_pl/geom

# MCMC: Kepler masses + wind shape + free scattered-flux floor
python mcmc_lightcurve_fit.py --band broad \
    --flux-csv flux_vs_nH_tbabs_broad.csv --data-dir data/IC_10_X1_LC_CIAO \
    --obs-column flux_t --time-column t_raw \
    --wind-model smooth_pl --fit-wind-shape --kepler --fit-scatter \
    --likelihood jitter --counts-per-bin 100 \
    --sampler zeus --n-walkers 24 --n-steps 21000 --n-burn 2000 \
    --n-threads 4 --dth 4.0 \
    --prior-MX 30,10,1,100 --prior-MRH 20,10,1,100 \
    --prior-Rb 6,3,3,80 --prior-p 4,2,2,8 \
    --compute-bic --output-dir mcmc_results/broad/smooth_pl/kepler_shape

# MCMC: raw unbinned 100s data (pair with jitter)
python mcmc_lightcurve_fit.py --band soft \
    --flux-csv flux_vs_nH_tbabs_soft.csv --data-dir data/IC_10_X1_LC_CIAO \
    --obs-column flux_t --time-column t_raw \
    --no-phase-bin --likelihood jitter --output-dir mcmc_results/soft/raw_jitter

# Freeze parameters out of the chain
python mcmc_lightcurve_fit.py --band broad --flux-csv flux_vs_nH_tbabs_broad.csv \
    --reparam --freeze q=0.5,Rb=6.0 --n-steps 2000 \
    --output-dir mcmc_results/broad/frozen

# Re-plot / recompute BIC from saved results (pass the same data/binning flags)
python mcmc_lightcurve_fit.py --band broad --flux-csv flux_vs_nH_tbabs_broad.csv \
    --data-dir data/IC_10_X1_LC_CIAO --obs-column flux_t --time-column t_raw \
    --wind-model smooth_pl --counts-per-bin 100 \
    --replot --compute-bic --output-dir mcmc_results/broad/smooth_pl/geom

Uncommitted work in progress (branch add_generic_wind): Phase 19 — Gaussian smoothing, f_scatter, and residual panels — is implemented in chandra_phase_analysis.py, mcmc_lightcurve_fit.py, xrb_lightcurve.py, and xrb_toy_wind_models.ipynb but not yet committed.

Last Updated: August 12, 2026
Maintainer: R. Panchal