Parametric SLM: emulating the TOV solver across EOS parameters
This tutorial shows how to train the SLM emulator of the TOV solver
(slmemulator.TOVEmulator) on a grid of equation-of-state (EOS) parameters and
then predict mass–radius and tidal curves at new parameter values in a
fraction of a millisecond — no TOV solve required at prediction time.
How it works. For each training parameter point the emulator generates the EOS, solves the TOV equations, and stores the log-transformed curves \([\log R,\ \log P_c,\ \log M,\ \log k_2]\) (the star log-extended structure of the SLM method, arXiv:2411.10556). The curves are then interpolated across min–max-normalized parameter space with the Banach Greedy Recombination Interpolation Method (Banach GRIM, arXiv:2205.07495): the dense kernel interpolant is thinned by recombination against a greedily grown set of evaluation functionals, which is exact at the matched functionals and smooth in between.
For a scripted version of this workflow (with accuracy/timing tables), see
tests/testpSLM.py.
1. Train the emulator
We train on a \(6\times3\) grid of the Quarkyonia EOS parameters
\((\Lambda, \kappa)\) — 18 parameter sets, as in the paper. The TOV curves vary
much more strongly with \(\Lambda\) (the central pressure at maximum mass roughly
triples from \(\Lambda=300\) to \(500\)) than with \(\kappa\), so \(\Lambda\) gets the
denser grid; spending the same 18 points \(\kappa\)-dense (\(3\times6\)) is no more
accurate than a \(3\times3\) grid. For tutorial speed we lower the TOV
integration steps to sol_pts=1000; use the default sol_pts=4000 for
production runs.
import numpy as np
import matplotlib.pyplot as plt
from time import perf_counter
from slmemulator import TOVEmulator, quarkyonia_generator
# training grid: 6 x 3 in (lambda, kappa) — 18 points, lambda-dense because
# the curves vary far more strongly with lambda than with kappa
lam_grid = np.linspace(300.0, 500.0, 6)
kappa_grid = np.linspace(0.1, 0.3, 3)
train_params = np.array([(l, k) for l in lam_grid for k in kappa_grid])
emulator = TOVEmulator(
quarkyonia_generator, # writes an EOS file for given (lambda, kappa)
work_dir="EOS_files/QUARKYONIA_TUTORIAL",
param_names=["lambda", "kappa"],
tidal=True,
sol_pts=1000, # tutorial speed; default is 4000
)
emulator.train(train_params)
print(f"training data: {emulator.train_data_time_:.1f} s "
f"({len(train_params)} TOV solves), fit: {emulator.fit_time_*1e3:.0f} ms")
2. Predict at new parameters and validate
predict maps parameters directly to the curves (radius in km, central
pressure in MeV/fm\(^3\), mass in \(M_\odot\), and the Love number \(k_2\)). We
validate against a direct TOV solve at a held-out point.
theta = [370.0, 0.17] # held-out: not on the training grid
curves_emulated = emulator.predict(theta) # instant
curves_true = emulator.solve(theta, tag="_test") # full physics (EOS + TOV)
names = emulator.quantity_names
rel_err = np.max(np.abs(curves_emulated - curves_true) / np.abs(curves_true), axis=1)
for name, err in zip(names, rel_err):
print(f"max relative error, {name:<18}: {err:.2e}")
The single most-quoted summary of a TOV solution is the maximum-mass configuration. Comparing it between the solver and the emulator gives an immediate feel for the accuracy:
R_emu, Pc_emu, M_emu, k2_emu = curves_emulated
R_true, Pc_true, M_true, k2_true = curves_true
# maximum-mass star: TOV solver vs pSLM emulator
i_tov, i_emu = np.argmax(M_true), np.argmax(M_emu)
print(f"{'':15}{'M_max [M_sun]':>15}{'R(M_max) [km]':>15}{'P_c(M_max) [MeV/fm^3]':>23}")
print(f"{'TOV solver':<15}{M_true[i_tov]:>15.4f}{R_true[i_tov]:>15.3f}{Pc_true[i_tov]:>23.2f}")
print(f"{'pSLM emulator':<15}{M_emu[i_emu]:>15.4f}{R_emu[i_emu]:>15.3f}{Pc_emu[i_emu]:>23.2f}")
print(f"\nmax-mass point differs by "
f"{abs(M_emu[i_emu] - M_true[i_tov]) / M_true[i_tov] * 100:.2f}% in mass, "
f"{abs(R_emu[i_emu] - R_true[i_tov]) / R_true[i_tov] * 100:.2f}% in radius")
fig, axes = plt.subplots(1, 2, figsize=(11, 4))
axes[0].plot(R_true, M_true, ".", label="TOV solver")
axes[0].plot(R_emu, M_emu, "-", label="pSLM emulator")
axes[0].plot(R_true[i_tov], M_true[i_tov], "k*", ms=12, label="TOV max mass")
axes[0].plot(R_emu[i_emu], M_emu[i_emu], "r+", ms=12, mew=2, label="pSLM max mass")
axes[0].set_xlabel("Radius [km]")
axes[0].set_ylabel(r"Mass [$M_\odot$]")
axes[0].legend()
axes[1].semilogy(M_true, k2_true, ".", label="TOV solver")
axes[1].semilogy(M_emu, k2_emu, "-", label="pSLM emulator")
axes[1].set_xlabel(r"Mass [$M_\odot$]")
axes[1].set_ylabel(r"$k_2$")
axes[1].legend()
fig.suptitle(rf"$\Lambda = {theta[0]:g}$, $\kappa = {theta[1]:g}$ (held out)")
fig.tight_layout()
plt.show()
3. Speedup
The one-time training cost is a TOV solve per grid point; after that each prediction is a kernel evaluation.
t0 = perf_counter()
emulator.solve(theta, tag="_timing")
t_tov = perf_counter() - t0
t0 = perf_counter()
for _ in range(100):
emulator.predict(theta)
t_emu = (perf_counter() - t0) / 100
print(f"TOV solver : {t_tov:.2f} s per parameter point")
print(f"emulator : {t_emu*1e3:.3f} ms per parameter point "
f"(~{t_tov/t_emu:,.0f}x speedup)")
4. Persist the trained emulator
Training is the expensive part, so save it; load restores the emulator
without re-solving anything.
5. SLM/DMD modes of an emulated curve
The emulated curve can be fed straight into the SLM decomposition (log-space DMD with quadratic augmentation) to obtain modes, eigenvalues and a reduced reconstruction — the parametric SLM (pSLM) pipeline end to end.
6. File-based workflow: ParametricSLM
If you already have TOV data files (columns: radius, central pressure, mass,
\(k_2\)) computed at known parameters, ParametricSLM fits the same Banach-GRIM
interpolant directly on the files. Its predict returns the familiar DMD
component tuple, computed from the interpolated curve.
from pathlib import Path
from slmemulator import ParametricSLM
# write the training curves as MR files (as tests/testSLM.py produces)
mr_dir = Path("TOV_data/QUARKYONIA_TUTORIAL")
mr_dir.mkdir(parents=True, exist_ok=True)
file_list = []
for p, X_log in zip(emulator.train_params, emulator.train_curves_log):
f = mr_dir / f"MR_{p[0]:.2e}_{p[1]:.2e}.txt"
np.savetxt(f, np.exp(X_log).T)
file_list.append(f.name)
pslm = ParametricSLM(file_list, filePath=mr_dir, tidal=True) # params parsed
pslm.fit() # from file names
Phi, omega, eigs, b, Xdmd, t = pslm.predict(theta)
print("predicted curve matrix:", Xdmd.shape, "| physical curves: np.exp(Xdmd.real)")
Practical notes
- Grid density matters — put points where the curves vary. Two points per
parameter (grid corners) cannot emulate the interior when the curves vary
strongly. Just as important is where the points go: for Quarkyonia the
curves change far more with \(\Lambda\) than with \(\kappa\), so the 18-point
\(6\times3\) (\(\Lambda\)-dense) grid used above reaches sub-percent errors in
\(R\), \(P_c\) and \(M\) (a few percent for \(k_2\), the hardest quantity), while
the same 18 points spent \(\kappa\)-dense (\(3\times6\)) are no better than a
\(3\times3\) grid. Prediction cost is one kernel-row evaluation, so it is
essentially independent of the training-set size — a denser grid only adds
one-time training solves. Validate with held-out points, e.g. via
python tests/testpSLM.py. - Centering. The interpolant works on deviations from the mean training
curve (
center=Trueinternally), so quantities that barely vary across the training set are reproduced (nearly) exactly. - Interpolation method.
BanachGRIMInterpolator.fituses the Banach GRIM algorithm (method="grim") by default;method="pgreedy"gives f·P-greedy kernel interpolation, which is faster for very large snapshot sets. At tight tolerance both reproduce the dense kernel interpolant. - Tolerances. For the emulator,
tolis the maximum entry-wise training residual in log space; the default1e-12keeps full interpolation accuracy.