"""
MAP (Maximum A Posteriori) optimizer for DE-MCMC warm-start.
Runs a bounded global optimization of the negative log-posterior (likelihood +
prior) produced by ModelData.lh_function_gib, then refines with a local
polish. The optimum replaces bestpars_data.list_bestpars_initial_value so
that DE-MCMC chains are initialized around the posterior mode instead of a
user-provided guess. This mitigates the DE-MCMC feedback loop where
unconstrained dimensions diffuse, (r1 - r2) grows unbounded, and joint
proposal size blows up causing acceptance collapse.
Intended to be called from Retrieval.main exactly once, before
run_mcmc_sampling, and only when the user enabled the GUI checkbox.
"""
from __future__ import annotations
import sys
import time
from multiprocessing.managers import SyncManager
from typing import List, Tuple
import numpy as np
from scipy.optimize import differential_evolution, minimize
from GUIBRUSHR.Retrieval.ModelCalculation.Classes.ParamForModel import RANGE_SENTINEL
from GUIBRUSHR.Retrieval.ExofastMCMC.process_safety import (
install_parent_death_signal as _install_parent_death_signal,
mp_context as _mp_context,
)
_INVALID_COST = 1.0e30
def _build_bounds_aligned_with_chain(model_obj) -> List[Tuple[float, float]]:
"""
Build per-chain-index (range_min, range_max) bounds.
The order matches the iteration used by
``ParameterHandler.create_param_full``: each chain entry maps to a
base parameter slot via the explicit ``chain_to_base`` routing
(built from ``list_bestpars`` + ``list_multiple_param``), so we no
longer rely on substring matching to attribute consecutive chain
entries to the right multi slot.
"""
from re import sub
param_handler = model_obj.param_handler
bestpars_data = model_obj.bestpars_data
list_bestpars = bestpars_data.list_bestpars
list_multiple_param = param_handler.list_multiple_param
params_list = param_handler.params_list
initial_param_array = param_handler.initial_param_array
bounds: List[Tuple[float, float]] = []
for chain_name in list_bestpars:
if chain_name in params_list:
base = chain_name
else:
base = sub(r"\d+", "", chain_name)
if base not in list_multiple_param:
raise RuntimeError(
f"MAP optimizer: chain entry {chain_name!r} stripped to "
f"{base!r}, which is not in list_multiple_param."
)
slot = initial_param_array[params_list.index(base)]
if slot is None:
raise RuntimeError(
f"MAP optimizer: no slot for base {base!r} (chain {chain_name!r})."
)
if not slot.status:
raise RuntimeError(
f"MAP optimizer: slot for base {base!r} (chain {chain_name!r}) "
"is not active (status=False); cannot supply MCMC bounds for a "
"fixed parameter."
)
lo = float(slot.range_min)
hi = float(slot.range_max)
# Unbounded ranges (+-inf, from the +-999999 sentinels) must be
# clamped back to a finite search box: differential_evolution
# requires finite bounds. The sentinel keeps the legacy window.
if not np.isfinite(lo):
lo = -float(RANGE_SENTINEL)
if not np.isfinite(hi):
hi = float(RANGE_SENTINEL)
bounds.append((lo, hi))
if len(bounds) != bestpars_data.nfit:
raise RuntimeError(
f"MAP optimizer: bounds length {len(bounds)} does not match "
f"nfit={bestpars_data.nfit}. Parameter alignment is broken."
)
return bounds
def _neg_log_post(theta: np.ndarray, model_obj) -> float:
"""Negative log-posterior objective, finite even for invalid points."""
try:
param_full = model_obj.create_param_full(np.asarray(theta, dtype=float))
except (AssertionError, ValueError, KeyError, IndexError):
return _INVALID_COST
for p in param_full:
if p is not None and not p.boundaries_check():
return _INVALID_COST
try:
lhood, log_prior, _, _ = model_obj.lh_function_gib(param_full)
except (ValueError, FloatingPointError, RuntimeError):
return _INVALID_COST
value = -(float(lhood) + float(log_prior))
if not np.isfinite(value):
return _INVALID_COST
return value
def _clip_to_bounds(
theta: np.ndarray, bounds: List[Tuple[float, float]]
) -> np.ndarray:
lows = np.array([b[0] for b in bounds], dtype=float)
highs = np.array([b[1] for b in bounds], dtype=float)
return np.minimum(np.maximum(theta, lows), highs)
def _chunk_worker(chunk_thetas, idx_offset, return_dict, model_obj, worker_seed):
"""Per-call worker: evaluate a static slice of trial vectors.
Direct mirror of ``ModelData.parallel_chain``: each Process receives its
own slice of the work (here a chunk of theta vectors instead of a chain
range), iterates through it autonomously, and writes results into a
``Manager.dict`` keyed by the global trial index. No IPC during execution;
parent collects after ``join()``.
"""
# SAFETY (requirement 1): forked worker must die if the parent crashes
# (segfault/OOM) instead of becoming a CPU-burning orphan. Linux-only.
_install_parent_death_signal()
if worker_seed is not None and getattr(model_obj, "random_obj", None) is not None:
# Same safeguard as parallel_chain (ModelData.py:714-720): replace the
# fork-inherited RNG with an independent stream from a SeedSequence
# child. Kept even though the current forward model is RNG-free.
model_obj.random_obj.rng = np.random.default_rng(worker_seed)
for k, theta in enumerate(chunk_thetas):
try:
cost = _neg_log_post(np.asarray(theta, dtype=float), model_obj)
except Exception:
cost = _INVALID_COST
return_dict[idx_offset + k] = float(cost)
def _parallel_evaluate(model_obj, n_workers: int, thetas: List[np.ndarray]) -> List[float]:
"""Evaluate ``thetas`` across ``n_workers`` processes, parallel_chain-style.
Splits the trial vectors into ``n_workers`` contiguous chunks, spawns one
``multiprocessing.Process`` per chunk, joins, then reconstructs the cost
list in original order. Lifecycle (spawn -> start -> join -> terminate) is
a 1:1 transcription of ``run_multiple_processes_manager_version`` in
``ModelData.py:1080-1146`` - same Manager.dict result protocol, same
``seed_seq.spawn`` for child seeds, same per-call spawn cost on macOS.
"""
n_trials = len(thetas)
if n_trials == 0:
return []
# Cap workers to actual work (avoid empty chunks creating idle processes).
n_workers = max(1, min(int(n_workers), n_trials))
chunk_sizes = [
len(c) for c in np.array_split(np.arange(n_trials), n_workers)
]
offsets = [0]
for s in chunk_sizes[:-1]:
offsets.append(offsets[-1] + s)
child_seeds: List = [None] * n_workers
if getattr(model_obj, "random_obj", None) is not None:
try:
child_seeds = list(model_obj.random_obj.seed_seq.spawn(n_workers))
except Exception:
child_seeds = [None] * n_workers
# Manager shut down in ``finally`` so its server process is reaped even if
# an exception escapes mid-loop. Without this, every DE iteration would
# leak a Manager process (40 iter * Manager = zombie accumulation).
# Context picked via _mp_context() to enforce ``fork`` on macOS as well
# (avoids B8: spawn pickle of pRT atmosphere = 53min-5.3hr overhead).
# Started with a parent-death initializer (requirement 1) so the Manager
# server also dies if the parent crashes instead of orphaning.
ctx = _mp_context()
proc: List = []
manager = SyncManager(ctx=ctx)
manager.start(_install_parent_death_signal)
try:
return_dict = manager.dict()
try:
for i in range(n_workers):
start = offsets[i]
end = start + chunk_sizes[i]
chunk = thetas[start:end]
proc.append(
ctx.Process(
target=_chunk_worker,
args=(chunk, start, return_dict, model_obj, child_seeds[i]),
)
)
proc[i].start()
for i in range(n_workers):
proc[i].join()
# Diagnostic: worker that died via SIGSEGV/SIGKILL/abnormal
# exit leaves no Python exception. The missing return_dict
# entries are silently substituted with _INVALID_COST below;
# log a WARNING so the user knows a crash happened (audit
# C5/C6/A8: silent worker death = invisible degradation).
ec = proc[i].exitcode
if ec is not None and ec != 0:
print(
f"[MAP optimizer] WARNING: worker {i} "
f"(chunk[{offsets[i]}:{offsets[i] + chunk_sizes[i]}]) "
f"exited abnormally with code {ec} "
f"(SIGSEGV=-11, SIGKILL/OOM=-9); "
f"affected trials get _INVALID_COST",
flush=True,
)
proc[i].terminate()
finally:
# Defensive cleanup on any exception escape (KeyboardInterrupt,
# MemoryError, etc.): terminate any worker still alive and reap
# it. Without this, SIGINT during _parallel_evaluate left up to
# n_workers zombie processes (audit C4 confirmed empirically).
for p in proc:
if p.is_alive():
p.terminate()
p.join(timeout=5.0)
if p.is_alive():
p.kill()
p.join(timeout=2.0)
# Reconstruct ordered list; missing entries (worker died) get _INVALID_COST.
results = [float(return_dict.get(i, _INVALID_COST)) for i in range(n_trials)]
finally:
manager.shutdown()
return results
[docs]
def compute_hessian_at_map(
model_obj,
x_map: np.ndarray,
bounds,
n_workers: int = 1,
h_rel: float = 1.0e-2,
eps_floor: float = 1.0e-10,
max_ridge_iter: int = 20,
full_matrix: bool = True,
):
"""Compute the Hessian of ``-log_post`` at ``x_map`` via central FD.
When ``full_matrix=True`` the off-diagonal terms are computed too
(``1 + 2N + 4·C(N,2)`` evaluations); when ``False`` only the diagonal is
populated (``1 + 2N`` evaluations) - same finite-difference machinery,
just skip the cross-term loop. This matches the user-facing
``init_mode=="diagonal"`` choice that explicitly opts out of off-diagonal
information.
Returns
-------
(H, C, sigma_p, L) or None
H: regularized Hessian (PSD), shape (d, d).
C: H^-1 (covariance approximation), shape (d, d).
sigma_p: sqrt(diag(C)), capped at ``0.5 * (range_max-range_min)``.
L: Cholesky factor of C (lower triangular), shape (d, d) or None
when ``full_matrix=False``.
Returns ``None`` on any non-recoverable failure (caller should fall
back to the legacy isotropic init and emit a JSONL fallback event).
"""
d = x_map.size
scale_vec = np.asarray(model_obj.retrieval_data.scale_vector_params, dtype=float)
if scale_vec.shape[0] != d:
# Aliased fallback: use prior widths if scale shape mismatches (rare).
lows_b = np.array([b[0] for b in bounds], dtype=float)
highs_b = np.array([b[1] for b in bounds], dtype=float)
scale_vec = 0.01 * (highs_b - lows_b)
h = h_rel * np.maximum(np.abs(scale_vec), 1.0)
# Shrink h if any stencil point would cross a prior bound.
lows = np.array([b[0] for b in bounds], dtype=float)
highs = np.array([b[1] for b in bounds], dtype=float)
for i in range(d):
room = min(x_map[i] - lows[i], highs[i] - x_map[i])
if room <= 0:
# MAP sits on a bound - Hessian is undefined here. Bail out.
return None
if h[i] >= room:
h[i] = 0.45 * room
# Build the ordered evaluation list deterministically.
thetas = [x_map.copy()]
diag_idx = {}
for i in range(d):
for s in (+1, -1):
t = x_map.copy()
t[i] += s * h[i]
diag_idx[(i, s)] = len(thetas)
thetas.append(t)
cross_idx = {}
if full_matrix:
for i in range(d):
for j in range(i + 1, d):
for si in (+1, -1):
for sj in (+1, -1):
t = x_map.copy()
t[i] += si * h[i]
t[j] += sj * h[j]
cross_idx[(i, j, si, sj)] = len(thetas)
thetas.append(t)
costs = _parallel_evaluate(model_obj, max(1, int(n_workers)), thetas)
f = -np.array(costs, dtype=float) # convert -log_post → log_post
if not np.all(np.isfinite(f)):
# Stencil hit a -inf region (opacity grid edge or invalid params);
# the resulting Hessian row would be poisoned. Bail out cleanly.
return None
f0 = f[0]
H = np.zeros((d, d), dtype=float)
for i in range(d):
fp = f[diag_idx[(i, +1)]]
fm = f[diag_idx[(i, -1)]]
H[i, i] = -(fp - 2 * f0 + fm) / (h[i] ** 2)
if full_matrix:
for i in range(d):
for j in range(i + 1, d):
fpp = f[cross_idx[(i, j, +1, +1)]]
fpm = f[cross_idx[(i, j, +1, -1)]]
fmp = f[cross_idx[(i, j, -1, +1)]]
fmm = f[cross_idx[(i, j, -1, -1)]]
H[i, j] = -(fpp - fpm - fmp + fmm) / (4.0 * h[i] * h[j])
H[j, i] = H[i, j]
H = 0.5 * (H + H.T) # symmetrize against round-off
# Regularize toward PSD via ridge bumping until Cholesky succeeds.
eig_min = np.linalg.eigvalsh(H).min() if full_matrix else float(np.diag(H).min())
eps = max(eps_floor, abs(eig_min) * 1.1) if eig_min <= 0 else 0.0
H_reg = H + eps * np.eye(d)
if full_matrix:
for _ in range(max_ridge_iter):
try:
np.linalg.cholesky(H_reg)
break
except np.linalg.LinAlgError:
eps = max(eps * 10.0, eps_floor)
H_reg = H + eps * np.eye(d)
else:
return None # could not regularize to PSD
# Diagonal-only path: H is (intentionally) only populated on the diag, so
# C is also diagonal. Avoids the inversion of a sparse matrix.
if not full_matrix:
diag_H = np.diag(H_reg)
if np.any(diag_H <= 0):
return None
sigma_p = 1.0 / np.sqrt(diag_H)
# Cap at half the prior width so orphan dimensions (H_ii ≈ 0 → σ_p
# blows up after ridge) don't blow chains past the bounds.
rmin_v = np.asarray(model_obj.retrieval_data.range_min_vector_params, dtype=float)
rmax_v = np.asarray(model_obj.retrieval_data.range_max_vector_params, dtype=float)
if rmin_v.shape[0] == d and rmax_v.shape[0] == d:
sigma_p = np.minimum(sigma_p, 0.5 * (rmax_v - rmin_v))
C = np.diag(sigma_p ** 2)
return H_reg, C, sigma_p, None
# Full-matrix path: invert, extract σ_p (capped), and compute Cholesky.
try:
C = np.linalg.inv(H_reg)
except np.linalg.LinAlgError:
return None
sigma_p = np.sqrt(np.maximum(np.diag(C), 0.0))
rmin_v = np.asarray(model_obj.retrieval_data.range_min_vector_params, dtype=float)
rmax_v = np.asarray(model_obj.retrieval_data.range_max_vector_params, dtype=float)
if rmin_v.shape[0] == d and rmax_v.shape[0] == d:
sigma_p = np.minimum(sigma_p, 0.5 * (rmax_v - rmin_v))
# Rebuild C consistent with capped σ_p (preserves correlations but
# rescales rows/columns of the affected dimensions). This keeps the
# Cholesky downstream consistent with the per-axis cap.
scale_diag = sigma_p / np.sqrt(np.maximum(np.diag(C), 1.0e-300))
C = (C * scale_diag[:, None]) * scale_diag[None, :]
try:
L = np.linalg.cholesky(C + 1.0e-12 * np.eye(d))
except np.linalg.LinAlgError:
return None
return H_reg, C, sigma_p, L
[docs]
def find_map_start(
model_obj,
maxiter_global: int = 40,
popsize: int = 10,
tol_global: float = 1.0e-2,
maxiter_local: int = 150,
seed: int | None = None,
n_workers: int = 1,
init_mode: str = "isotropic",
) -> np.ndarray:
"""
Find the MAP estimate and overwrite bestpars_data.list_bestpars_initial_value.
Uses scipy.optimize.differential_evolution for global exploration followed
by a Nelder-Mead polish. Chain bounds and the current initial value are
read from model_obj. On success, list_bestpars_initial_value is replaced
in place with the MAP point so that subsequent DE-MCMC chain
initialization uses the posterior mode as its center.
Parameters
----------
model_obj :
Fully initialized ModelData object (bestpars_data, retrieval_data,
param_handler available).
maxiter_global : int
Maximum iterations for differential_evolution. Default 40.
popsize : int
Population size multiplier for differential_evolution. Default 10.
tol_global : float
Relative tolerance for differential_evolution convergence. Default 1e-2.
maxiter_local : int
Maximum iterations for Nelder-Mead polish. Set to 0 to skip the
polish stage entirely. Default 150.
seed : int or None
RNG seed for reproducibility; None pulls from OS entropy.
n_workers : int
Number of multiprocessing.Process workers used to evaluate the DE
population in parallel. ``n_workers <= 1`` keeps the historical
serial path. With ``n_workers > 1`` the function spawns a persistent
pool mirroring ``ModelData.parallel_chain`` (mp.Process + child seeds
from ``model_obj.random_obj.seed_seq``). On any pool failure the
function falls back to the serial path with the same init population.
Returns
-------
numpy.ndarray
The MAP parameter vector (length nfit).
"""
bestpars_data = model_obj.bestpars_data
table_output_file = model_obj.retrieval_data.table_output_file
bounds = _build_bounds_aligned_with_chain(model_obj)
x0 = np.array(bestpars_data.list_bestpars_initial_value, dtype=float)
x0 = _clip_to_bounds(x0, bounds)
cost_initial = _neg_log_post(x0, model_obj)
def _log(msg: str) -> None:
print(msg)
if table_output_file is not None:
try:
# Append (not truncate): MAP runs many iterations and the
# historical log on disk is the only post-mortem artifact -
# E7 confirmed `"w"` was destroying all but the last line.
with open(table_output_file, "a") as f:
f.write(msg + "\n")
except OSError:
pass
# Display convention: we minimize -log_post internally (scipy needs a
# cost to minimize), but we log log_post directly so the numbers behave
# bayesian-naturally (higher = better fit + prior). Δlog_post reports
# the gain of each stage relative to its reference point, which is the
# only quantity that matters for judging whether the optimizer worked.
_log("[MAP optimizer] starting warm-start search")
_log(f"[MAP optimizer] nfit = {bestpars_data.nfit}")
_log(f"[MAP optimizer] initial log_post = {-cost_initial:.12g}")
# Shared progress state. nfev and best are cumulative across DE +
# polish; iter counters are stage-local.
state = {"nfev": 0, "best": float("inf"), "de_iter": 0, "polish_iter": 0}
log_every_eval = max(50, popsize * bestpars_data.nfit // 2)
def _objective(theta: np.ndarray) -> float:
value = _neg_log_post(theta, model_obj)
state["nfev"] += 1
if value < state["best"]:
state["best"] = value
if state["nfev"] % log_every_eval == 0:
_log(
f"[MAP optimizer] nfev = {state['nfev']}, "
f"best log_post so far = {-state['best']:.6g}"
)
return value
def _de_callback(x: np.ndarray, convergence: float) -> bool:
state["de_iter"] += 1
_log(
f"[MAP optimizer] DE iter {state['de_iter']}: "
f"best log_post = {-state['best']:.6g}, "
f"convergence = {convergence:.4g}, nfev = {state['nfev']}"
)
return False
def _polish_callback(xk: np.ndarray) -> bool:
state["polish_iter"] += 1
if state["polish_iter"] % 20 == 0:
_log(
f"[MAP optimizer] polish iter {state['polish_iter']}: "
f"best log_post = {-state['best']:.6g}, "
f"nfev = {state['nfev']}"
)
return False
init_population = np.tile(x0, (max(popsize * bestpars_data.nfit, 5), 1))
rng = np.random.default_rng(seed)
lows = np.array([b[0] for b in bounds], dtype=float)
highs = np.array([b[1] for b in bounds], dtype=float)
spread = 0.1 * (highs - lows)
init_population = init_population + rng.normal(
scale=spread, size=init_population.shape
)
init_population = np.clip(init_population, lows, highs)
init_population[0] = x0
# scipy >= 1.15 deprecates ``seed=`` in favor of ``rng=`` (removal in 1.17).
# Both accept int | Generator | SeedSequence; ``rng=`` is the supported form.
de_kwargs = dict(
bounds=bounds,
maxiter=maxiter_global,
popsize=popsize,
tol=tol_global,
init=init_population,
polish=False,
rng=seed,
updating="deferred",
callback=_de_callback,
)
n_workers = max(1, int(n_workers or 1))
_log(
f"[MAP optimizer] DE: n_workers={n_workers}, "
f"start_method={_mp_context().get_start_method()} "
f"(platform={sys.platform})"
)
# parallel_map: scipy hands us a population per DE iteration; we split it
# into n_workers chunks and spawn N Process per call (parallel_chain
# pattern: spawn -> start -> join -> terminate, Manager.dict for
# results). ``func`` is scipy's wrapped objective; we ignore it and call
# ``_neg_log_post`` directly inside the worker.
#
# State (nfev, best) MUST be updated parent-side here: with fork, worker
# processes mutate their own copy of ``state`` and changes never propagate
# back. The serial path goes through ``_objective`` which updates
# ``state`` directly; the parallel path bypasses it.
def parallel_map(func, iterable): # noqa: ARG001
items = list(iterable)
costs = _parallel_evaluate(model_obj, n_workers, items)
for c in costs:
state["nfev"] += 1
if c < state["best"]:
state["best"] = c
if state["nfev"] % log_every_eval == 0:
_log(
f"[MAP optimizer] nfev = {state['nfev']}, "
f"best log_post so far = {-state['best']:.6g}"
)
return costs
workers_arg = parallel_map if n_workers > 1 else 1
de_result = None # ensure defined even if both attempts raise
de_t0 = time.perf_counter()
try:
try:
de_result = differential_evolution(
_objective, workers=workers_arg, **de_kwargs
)
except Exception as exc_parallel: # pragma: no cover - safety net
if workers_arg == 1:
raise
_log(
f"[MAP optimizer] parallel evaluation failed "
f"({type(exc_parallel).__name__}: {exc_parallel}); "
"falling back to serial DE with the same init population"
)
try:
de_result = differential_evolution(
_objective, workers=1, **de_kwargs
)
except Exception as exc_fallback:
_log(
f"[MAP optimizer] serial fallback also failed "
f"({type(exc_fallback).__name__}: {exc_fallback})"
)
raise RuntimeError(
"MAP optimizer: both parallel and serial DE failed"
) from exc_parallel
finally:
de_wall = time.perf_counter() - de_t0
_log(f"[MAP optimizer] DE wall-clock = {de_wall:.1f}s")
x_de = _clip_to_bounds(np.asarray(de_result.x, dtype=float), bounds)
cost_de = _neg_log_post(x_de, model_obj)
delta_de = cost_initial - cost_de
_log(
f"[MAP optimizer] DE finished: log_post = {-cost_de:.12g}, "
f"Δlog_post vs initial = {delta_de:+.6g} "
f"(nit={de_result.nit}, nfev={de_result.nfev})"
)
if maxiter_local > 0:
polish_result = minimize(
_objective,
x_de,
method="Nelder-Mead",
callback=_polish_callback,
options={
"xatol": 1.0e-4,
"fatol": 1.0e-3,
"maxiter": maxiter_local,
"adaptive": True,
},
)
x_polished = _clip_to_bounds(np.asarray(polish_result.x, dtype=float), bounds)
cost_polished = _neg_log_post(x_polished, model_obj)
delta_polish_vs_de = cost_de - cost_polished
delta_polish_vs_initial = cost_initial - cost_polished
_log(
f"[MAP optimizer] polish finished: log_post = {-cost_polished:.12g}, "
f"Δlog_post vs DE = {delta_polish_vs_de:+.6g}, "
f"Δlog_post vs initial = {delta_polish_vs_initial:+.6g} "
f"(nit={polish_result.nit}, nfev={polish_result.nfev})"
)
candidates = [(cost_initial, x0), (cost_de, x_de), (cost_polished, x_polished)]
else:
_log("[MAP optimizer] polish skipped (maxiter_local <= 0)")
candidates = [(cost_initial, x0), (cost_de, x_de)]
cost_best, x_best = min(candidates, key=lambda t: t[0])
delta_best = cost_initial - cost_best
_log(
f"[MAP optimizer] best log_post = {-cost_best:.12g}, "
f"total Δlog_post = {delta_best:+.6g}"
)
if cost_best >= _INVALID_COST:
_log(
"[MAP optimizer] WARNING: no valid point found. "
"Keeping original list_bestpars_initial_value."
)
return x0
if not np.all(np.isfinite(x_best)):
# MCMC downstream cannot tolerate NaN/inf in the starting point -
# better to keep the user's guess than to poison the chains silently.
_log(
"[MAP optimizer] WARNING: best candidate contains non-finite "
"values; keeping original list_bestpars_initial_value."
)
return x0
bestpars_data.list_bestpars_initial_value = x_best.tolist()
_log("[MAP optimizer] list_bestpars_initial_value replaced with MAP estimate")
# Optional Hessian extraction at the MAP, gated by the user-facing
# ``init_mode``. ``"isotropic"`` keeps the byte-identical legacy path
# (chains scattered with the YAML ``scale``); ``"diagonal"`` and
# ``"correlated"`` pay the FD cost to populate Bestpars attributes that
# exofast_demc's chain-init dispatch consumes. Failure here is non-fatal
# - the dispatch falls back to isotropic and emits a JSONL event.
bestpars_data.init_scale_per_param = None
bestpars_data.init_cholesky = None
mode = (init_mode or "isotropic").strip().lower()
if mode in ("diagonal", "correlated"):
full = (mode == "correlated")
d = x_best.size
n_evals = 1 + 2 * d + (4 * d * (d - 1) // 2 if full else 0)
_log(f"[MAP optimizer] Hessian (mode={mode}): "
f"~{n_evals} forward-model evals, {n_workers} workers")
h_t0 = time.perf_counter()
try:
result = compute_hessian_at_map(
model_obj, x_best, bounds,
n_workers=n_workers, full_matrix=full,
)
except Exception as exc: # pragma: no cover - safety net
result = None
_log(f"[MAP optimizer] Hessian extraction raised "
f"{type(exc).__name__}: {exc}")
h_wall = time.perf_counter() - h_t0
if result is None:
_log(f"[MAP optimizer] Hessian extraction FAILED in {h_wall:.1f}s "
f"- chain init will fall back to isotropic with a warning")
else:
H_reg, C, sigma_p, L = result
bestpars_data.init_scale_per_param = np.asarray(sigma_p, dtype=float)
if full and L is not None:
bestpars_data.init_cholesky = np.asarray(L, dtype=float)
try:
cond = float(np.linalg.cond(H_reg))
except np.linalg.LinAlgError:
cond = float("nan")
_log(f"[MAP optimizer] Hessian OK in {h_wall:.1f}s - "
f"cond(H)={cond:.3e}, σ_p sample={np.array2string(sigma_p, precision=3)}")
# JSONL diagnostic for the analyzer's
# init_sigma_vs_observed_sigma plot. Best-effort: never raise.
try:
from GUIBRUSHR.Retrieval.debug_log import emit_event, get_log_path
log_path = get_log_path(model_obj)
emit_event(log_path, {
"event": "init_covariance",
"mode": mode,
"sigma_p": [float(v) for v in sigma_p],
"cond_H": cond,
"n_evals": int(n_evals),
})
except Exception:
pass
return x_best