Inverse Model Demo

demos/demo_inverse.py synthesises observations from known solar parameters (with 0.02 mag Gaussian noise per band), runs the MCMC inverse model on them, prints the posterior summary, saves the chain to inverse_posterior.npz, and produces two diagnostic figures: a corner-style plot (inverse_corner.png, with the truth overplotted in red) and per-parameter chain traces (inverse_chains.png).

To fit real data instead, set SYNTHESISE = False and fill in the OBS_* arrays — the rest of the script is unchanged.

Concepts demonstrated: the inverse model, walker initialisation, posterior persistence and summaries.

Source

"""
demo_inverse.py
===============
Demonstration of the SED_Model inverse model.

Given observed magnitudes (with uncertainties) in a set of filters,
recovers the posterior distribution on (Teff, logg, [M/H]) via MCMC.

Usage
-----
    python demo_inverse.py

Edit the OBSERVATIONS and DATA PATHS sections below to suit your setup.
The demo generates synthetic observations from known parameters so it
works out of the box.  Replace the OBS_* arrays with your real data.
"""

import os
import sys
from pathlib import Path

import numpy as np
import matplotlib.pyplot as plt

from sed_model import (
    load_grid, load_filters_from_instrument_dir,
    run_forward, run_inverse, InverseResult,
)

sys.path.insert(0, str(Path(__file__).parent))
from plot_utils import plot_posterior

# =============================================================================
# DATA PATHS  —  edit these
# =============================================================================

GRID_DIR    = "~/SED_Tools/data/stellar_models/Kurucz2003all"
FILTER_DIR  = "~/SED_Tools/data/filters/Generic/Johnson"
VEGA_SED    = "~/SED_Tools/data/stellar_models/vega_flam.csv"

# =============================================================================
# FIXED PARAMETERS
# =============================================================================

R      = 6.957e10    # Stellar radius (cm)  —  must be supplied by the user
D      = 3.0857e19   # Distance (cm)        —  10 pc

MAG_SYSTEM = "AB"

# =============================================================================
# OBSERVATIONS  —  replace with your real data
#
# If you have real magnitudes, set:
#   OBS_FILTER_NAMES  = ["B", "V", "R"]  (must match filter file stems)
#   OBS_MAGNITUDES    = [12.3, 11.8, 11.4]
#   OBS_UNCERTAINTIES = [0.02, 0.02, 0.03]
#
# Below we synthesise observations from known parameters for demonstration.
# =============================================================================

SYNTHESISE = True   # set False and fill arrays below to use real data

OBS_FILTER_NAMES  = None   # e.g. ["B", "V", "R", "I"]
OBS_MAGNITUDES    = None   # e.g. [12.3, 11.8, 11.4, 11.0]
OBS_UNCERTAINTIES = None   # e.g. [0.02, 0.02, 0.03, 0.03]

# =============================================================================
# MCMC SETTINGS
# =============================================================================

N_WALKERS = 32
N_STEPS   = 1000
N_BURN    = 300
N_THIN    = 2

# =============================================================================
# LOAD DATA
# =============================================================================

GRID_DIR   = os.path.expanduser(GRID_DIR)
FILTER_DIR = os.path.expanduser(FILTER_DIR)
VEGA_SED   = os.path.expanduser(VEGA_SED)

print("Loading atmosphere grid ...")
grid = load_grid(GRID_DIR)
print(f"  {grid}")

print("\nLoading filters ...")
vega_path = VEGA_SED if os.path.isfile(VEGA_SED) else None
filters = load_filters_from_instrument_dir(FILTER_DIR, vega_sed_path=vega_path)
print(f"  {len(filters)} filters: {[f.name for f in filters]}")

# =============================================================================
# SYNTHESISE OBSERVATIONS (demo only)
# =============================================================================

TRUE_TEFF = 5778.0
TRUE_LOGG = 4.44
TRUE_META = 0.0

if SYNTHESISE:
    print(f"\nSynthesising observations from "
          f"Teff={TRUE_TEFF:.0f} K, logg={TRUE_LOGG:.2f}, [M/H]={TRUE_META:.2f} ...")
    fwd = run_forward(
        teff=TRUE_TEFF, logg=TRUE_LOGG, meta=TRUE_META,
        R=R, d=D, grid=grid, filters=filters, mag_system=MAG_SYSTEM,
    )
    rng = np.random.default_rng(42)
    OBS_FILTER_NAMES  = [f.name for f in filters]
    OBS_UNCERTAINTIES = np.full(len(filters), 0.02)
    OBS_MAGNITUDES    = (np.array([fwd.magnitudes[n] for n in OBS_FILTER_NAMES])
                         + rng.normal(0, OBS_UNCERTAINTIES))
    print("  Synthesised magnitudes (with 0.02 mag noise):")
    for name, mag, err in zip(OBS_FILTER_NAMES, OBS_MAGNITUDES, OBS_UNCERTAINTIES):
        print(f"    {name:>6s} : {mag:.4f} ± {err:.4f}")

# =============================================================================
# RUN INVERSE MODEL
# =============================================================================

print(f"\nRunning MCMC inference ...")
print(f"  Walkers={N_WALKERS}   Steps={N_STEPS}   Burn={N_BURN}   Thin={N_THIN}")
print(f"  This will take ~{N_WALKERS * N_STEPS / 1000:.0f}k forward evaluations ...")

posterior = run_inverse(
    obs_magnitudes=OBS_MAGNITUDES,
    obs_uncertainties=OBS_UNCERTAINTIES,
    filter_names=OBS_FILTER_NAMES,
    R=R, d=D,
    grid=grid,
    filters=filters,
    mag_system=MAG_SYSTEM,
    n_walkers=N_WALKERS,
    n_steps=N_STEPS,
    n_burn=N_BURN,
    n_thin=N_THIN,
    progress=True,
    seed=42,
)

# =============================================================================
# PRINT SUMMARY
# =============================================================================

posterior.print_summary()

# =============================================================================
# SAVE POSTERIOR
# =============================================================================

posterior.save("inverse_posterior.npz")
print("Saved: inverse_posterior.npz")
print("  Reload with: posterior = InverseResult.load('inverse_posterior.npz')")

# =============================================================================
# PLOT — corner plot + chain convergence
# =============================================================================

truths = (
    {'teff': TRUE_TEFF, 'logg': TRUE_LOGG, 'meta': TRUE_META}
    if SYNTHESISE else None
)
plot_posterior(posterior, "inverse_corner.png", "inverse_chains.png", truths=truths)