Sampling flux vacua

Sampling flux vacua#

What’s in this notebook? This notebook collects the workflows for sampling SUSY and non-SUSY flux vacua at scale, built on top of the ISD principle introduced in NB06. We walk through three complementary workflows (canonical SUSY from input fluxes, the wrapper API, the non-SUSY variation), compare the available sampling methods on a single tadpole-bounded benchmark, and end with an advanced performance-tuning section.

In this notebook, you will learn:

How to refine an integer-rounded ISD flux to a true SUSY vacuum via newton_method_flux_vacua (Workflow A).

How to drive the same pipeline through the high-level sample_SUSY_flux_vacua wrapper (Workflow B).

How the non-SUSY variation differs in mode flag, sign conventions, and result analysis (Workflow C).

How the four ISD-shift methods compare on tadpole regimes and what the bounding-box scaling looks like (Method comparison).

How to tune sampling regions, \(N_{\max}\) knobs, and hybrid solvers for production runs (Advanced).

Prerequisites: NB04 — Sampling module, NB06 — ISD sampling principle, NB05 — Finding flux vacua.

Recommended paths through this notebook#

This tutorial is deliberately broad because ISD sampling is used both for small demonstrations and for larger searches. The most common routes are:

If you want to…	Read
understand the minimal SUSY pipeline	Workflow A, then the take-aways
use the high-level production wrapper	Workflow B
compare SUSY and non-SUSY sampling	Workflows A and C
tune performance or method choices	Method comparison and advanced sections

A shorter pipeline-building version of this material would keep only Workflow A, Workflow B and the final take-aways. The remaining sections are kept here as a reference for method choices and diagnostics.

Outline#

Setup
Workflow A: SUSY ISD sampling from input fluxes
Workflow B: SUSY ISD sampling via the wrapper API
Workflow C: Non-SUSY ISD sampling
Method comparison and benchmarks
Advanced: performance tuning
Take-aways
Further reading

Setup#

# General imports
import warnings, sys, time
import numpy as np
from tqdm.auto import tqdm
from scipy.optimize import root

# JAX imports
import jax
from jax import jit, vmap
import jax.numpy as jnp
jax.config.update("jax_enable_x64", True)

# Plotting
import seaborn as sn
import matplotlib.pyplot as plt
cmap = sn.color_palette("viridis", as_cmap=True)

# JAXVacua
import jaxvacua as jvc
from jaxvacua.util import vmapping_func

warnings.filterwarnings("ignore")

Canonical fixture. Throughout this notebook we use the degree-18 hypersurface in \(\mathbb{CP}^{1,1,1,6,9}\) (\(h^{2,1}=2\), \(Q_{D3}=276\), model_ID=1) as the running example. The same fixture is reused in every workflow and benchmark so that results are directly comparable. The h12=3 mini-examples at the end of §A and §C illustrate that the workflows generalise without code changes.

h12 = 2
model = jvc.FluxVacuaFinder(
    h12=h12,
    Q=276,
    model_ID=1,
    maximum_degree=2,
    limit="LCS",
    model_type="KS",
    map_to_fd=True,   # map returned tau to the SL(2,Z) fundamental domain
)

We build a data_sampler with mildly restrictive flux and moduli bounds that we then reuse across §A, §B, §C, and the benchmark in §5. The construction is identical to that of NB04 — the sampler is purely infrastructure here.

sampler = jvc.data_sampler(
    model,
    flux_bounds=[-5, 5],
    moduli_bounds=[0, 10],
    axion_bounds=[-0.5, 0.5],
)

For convenient ensemble evaluation we vectorise the F-term function model.DW once and reuse it throughout. The cross-link to the ISD principle: the algorithmic intro to ISD sampling lives in NB06; here we focus on the workflow of producing vacua at scale.

DW_v = vmap(model.DW)

Workflow A: SUSY ISD sampling from input fluxes#

Workflow A starts from a fixed integer flux choice and searches for all moduli configurations that yield a SUSY vacuum compatible with it. This is the complement of Workflow B (moduli-first); together they cover the two natural ways of populating an ensemble. The algorithmic primitive is the same in both — the ISD-rounded shift from NB06 ISD principle — but the order in which fluxes and moduli are sampled differs.

The algorithm#

The pipeline has two stages, each vmap’d twice (once over moduli starting points and once over candidate flux vectors):

find_solution_init_vmap — apply model.linearised_shifts to round the input flux into an ISD-compatible integer flux and produce near-vacuum starting moduli.
find_solution_steps_vmap — iterate model.linearised_shifts to drive the residual \(|D_I W|\) to zero (Newton-like).

Below we wire up both stages, run a few iterations on two flux samples, and confirm the residual shrinks at each step.

# Common kwargs for the optimiser. "Fflux" fixes the RR-fluxes as input;
# the linearised_shifts call then solves for the NS-NS half.
mode = "Fflux"
kwargs = {"mode": mode, "Q": model.D3_tadpole,
          "constraints": None, "remove_NANs": True}

# Double vmap of model.DW — once over moduli starting points, once over fluxes.
DW_vv = jax.vmap(jax.vmap(model.DW))

The two doubly-vmapped optimisers — one for the ISD initialisation step (find_solution_init_vmap) and one for the Newton refinement step (find_solution_steps_vmap) — are built with the vmapping_func helper imported in §1:

# Step 1: ISD initialisation. Outer vmap: moduli starting points. Inner: same flux.
find_solution_init = vmapping_func(
    model.linearised_shifts, in_axes=(0, 0, None), return_flag=False, **kwargs)
find_solution_init_vmap = vmapping_func(
    find_solution_init, in_axes=(None, None, 0))

# Step 2: Newton refinement. Both axes vmapped.
find_solution_steps = vmapping_func(
    model.linearised_shifts, in_axes=(0, 0, 0), return_flag=True, **kwargs)
find_solution_steps_vmap = vmapping_func(
    find_solution_steps, in_axes=(0, 0, 0))

We supply two starting points and two candidate flux vectors. The doubly-vmapped call evaluates every (moduli, flux) pair in parallel.

moduli0 = np.array([
    [ 0.22924045 + 3.00550378j, -0.46016314 + 2.28651085j],
    [-0.04468382 + 0.63834065j,  0.25839392 + 2.36493542j],
])
tau0    = np.array([0.18254317 + 6.81361114j, 0.49828078 + 8.33716844j])
fluxes0 = np.array([[ 1,  5,  1,  2,  4, -1],
                    [ 5, -3,  0, -4, -4, -3]])

# Step 1: ISD initialisation rounds the fluxes and shifts the moduli.
moduli, tau, fluxes = find_solution_init_vmap(moduli0, tau0, fluxes0)

Iterating the Newton-like step find_solution_steps_vmap then drives \(|D_I W|\) toward zero. We track the per-flux residual at each step:

moduli1, tau1, fluxes1, _ = find_solution_steps_vmap(moduli, tau, fluxes)

for i in range(25):
    moduli1, tau1, fluxes1, _ = find_solution_steps_vmap(moduli1, tau1, fluxes1)
    DW = DW_vv(moduli1, jnp.conj(moduli1), tau1, jnp.conj(tau1), fluxes1)
    res = np.sum(np.abs(DW), axis=2)
    print(f"step {i:2d}   |DW| flux1: {res[0]}   flux2: {res[1]}")

Large-scale workflow#

All of the above is packaged in a single wrapper, model.sample_SUSY_vacua_from_fluxes, which combines the initialisation, the Newton refinement, and the bookkeeping. Below we sample 2,000 starting configurations and refine them in one call (takes a couple of minutes):

rns_key = jax.random.PRNGKey(42)

moduli, tau, fluxes, residual = model.sample_SUSY_vacua_from_fluxes(
    fluxes_init=fluxes0,
    N=10**3,
    sampler=sampler,
    rns_key=rns_key,
    moduli_sampling_mode="cone",
    max_iters=10,
    objective_fct=DW_vv,
    optimiser_init=find_solution_init_vmap,
    optimiser_steps=find_solution_steps_vmap,
    vmap_dim_pts=100,
    mode=mode,
    print_progress=True,
)

DW = jax.vmap(model.DW)(moduli, jnp.conj(moduli), tau, jnp.conj(tau), fluxes)
print(f"all |DW| residuals below 1e-10: "
      f"{np.all(np.sum(np.abs(DW), axis=1) < 1e-10)}")

Improved constrained sampling — input fluxes across a moduli region#

By passing fluxes_init=None we instead let the sampler draw fluxes internally, while still constraining the moduli to a chosen sub-region of the Kähler cone. This is the flux-first analog of moduli-first sampling: instead of fixing fluxes and searching moduli, we fix a region of moduli space and let the sampler explore fluxes consistent with it.

moduli, tau, fluxes, residuals = model.sample_SUSY_vacua_from_fluxes(
    fluxes_init=None,
    N=100,
    sampler=sampler,
    rns_key=rns_key,
    max_iters=10,
    moduli_sampling_mode="cone",
    max_tadpole=276,
    objective_fct=DW_vv,
    optimiser_init=find_solution_init_vmap,
    optimiser_steps=find_solution_steps_vmap,
    vmap_dim_flux=100,
    vmap_dim_pts=100,
    mode=mode,
    print_progress=True,
)

Higher-\(h^{1,2}\) example: SUSY workflow at \(h^{1,2}=3\)#

The workflow above generalises without code changes to higher \(h^{1,2}\). As an illustration we run the same Workflow-A pipeline on a \(h^{1,2}=3\) model (model_ID=4), reusing the doubly-vmapped optimisers built above — only the model and sampler need refreshing.

# h12=3 example (model_ID=4 has hyperplanes + sensible sampling region).
model_3 = jvc.FluxVacuaFinder(h12=3, model_ID=4, model_type="KS",
                              maximum_degree=2, limit="LCS")
sampler_3 = jvc.data_sampler(
    model_3,
    flux_bounds=[-5, 5],
    moduli_bounds=(2., 8.),
    dilaton_bounds=(np.sqrt(3) / 2, 10.),
    axion_bounds=(-0.5, 0.5),
    seed=42,
)

# Rebuild the optimisers with the new model.
kwargs_3 = {**kwargs, "Q": model_3.D3_tadpole}
_init_3 = vmapping_func(model_3.linearised_shifts, in_axes=(0, 0, None),
                       return_flag=False, **kwargs_3)
init_3_vmap = vmapping_func(_init_3, in_axes=(None, None, 0))
_step_3 = vmapping_func(model_3.linearised_shifts, in_axes=(0, 0, 0),
                       return_flag=True, **kwargs_3)
step_3_vmap = vmapping_func(_step_3, in_axes=(0, 0, 0))
DW_vv_3 = jax.vmap(jax.vmap(model_3.DW))

moduli3, tau3, fluxes3, res3 = model_3.sample_SUSY_vacua_from_fluxes(
    fluxes_init=None,
    N=500,
    sampler=sampler_3,
    rns_key=jax.random.PRNGKey(42),
    max_iters=10,
    moduli_sampling_mode="cone",
    max_tadpole=276,
    objective_fct=DW_vv_3,
    optimiser_init=init_3_vmap,
    optimiser_steps=step_3_vmap,
    vmap_dim_flux=50,
    vmap_dim_pts=10,
    mode=mode,
    print_progress=True,
)
DW3 = jax.vmap(model_3.DW)(moduli3, jnp.conj(moduli3),
                          tau3, jnp.conj(tau3), fluxes3)
print(f"h12=3 — {len(moduli3)} vacua; "
      f"all |DW| < 1e-10: {np.all(np.sum(np.abs(DW3), axis=1) < 1e-10)}")

Workflow B: SUSY ISD sampling via the wrapper API#

What’s new vs §A. §A wired the two-stage pipeline manually: build two doubly-vmapped optimisers (find_solution_init_vmap and find_solution_steps_vmap), call them in sequence, iterate. §B wraps the same pipeline behind a single high-level driver, model.sample_SUSY_flux_vacua, that handles the moduli draws, ISD rounding, and Newton refinement in one call. The trade-off: less fine-grained control over the inner loop, but a much shorter user-facing snippet that suffices for routine production runs.

ISD sampling with a single optimiser#

The simplest usage passes one optimiser — here linearised_shifts_F (constructed once with vmapping_func) — together with the §4 sampler. vmap_dim controls the batch size of the inner XLA kernel.

rns_key = jax.random.PRNGKey(42)

kwargs = {"Q": model.D3_tadpole, "return_flag": True,
          "constraints": None, "remove_NANs": True,
          "in_axes": (0, 0, 0)}

linearised_shifts_F = vmapping_func(
    model.linearised_shifts, mode="Fflux", **kwargs)

# Single-optimiser run targeting N=1000 SUSY vacua.
moduli, tau, fluxes, residuals = model.sample_SUSY_flux_vacua(
    N=10**3,
    rns_key=rns_key,
    sampler=sampler,
    optimiser=linearised_shifts_F,
    objective_fct=DW_v,
    vmap_dim=10**2,
    print_progress=True,
)

Running with multiple optimisers#

Different ISD modes ("ISD", "Hflux", "Fflux") cover different halves of the flux vector. Passing a list as optimisers=... lets the wrapper round-robin between them and pool the resulting vacua — handy when one mode alone yields too few convergent points.

linearised_shifts_ISD = vmapping_func(
    model.linearised_shifts, mode="ISD",   **kwargs)
linearised_shifts_H   = vmapping_func(
    model.linearised_shifts, mode="Hflux", **kwargs)
linearised_shifts_F   = vmapping_func(
    model.linearised_shifts, mode="Fflux", **kwargs)

optimisers = [linearised_shifts_ISD, linearised_shifts_H, linearised_shifts_F]

moduli, tau, fluxes, residuals = model.sample_SUSY_flux_vacua(
    N=10**2,
    rns_key=rns_key,
    sampler=sampler,
    optimisers=optimisers,
    objective_fct=DW_v,
    vmap_dim=10**2,
    print_progress=True,
)

Constrained sampling — moduli-region restrictions#

The wrapper accepts a JIT-compiled boolean constraints(moduli, tau, flux) callback that filters vacua at every batch. As a proof of concept we restrict the moduli to \(2<\operatorname{Im}(z^i)<3\) and the axio-dilaton to \(\operatorname{Im}(\tau)>4\), and rebuild the sampler with matching bounds so the initial guesses sit in the feasible region.

@jit
def constraints_model(moduli, tau, flux):
    return (jnp.all(jnp.imag(moduli) < 3., axis=0)
            & jnp.all(jnp.imag(moduli) > 2., axis=0)
            & (jnp.imag(tau) > 4))

kwargs_C = {"Q": model.D3_tadpole, "return_flag": True,
            "constraints": constraints_model, "remove_NANs": True,
            "in_axes": (0, 0, 0)}
linearised_shifts_F_v = vmapping_func(
    model.linearised_shifts, mode="Fflux", **kwargs_C)

# Tighter sampler matching the constraint region.
sampler_constrained = jvc.data_sampler(
    model, flux_bounds=[-5, 5],
    moduli_bounds=[2, 3], dilaton_bounds=[4, 10])

moduli, tau, fluxes, residuals = model.sample_SUSY_flux_vacua(
    N=10**2,
    rns_key=rns_key,
    sampler=sampler_constrained,
    optimiser=linearised_shifts_F_v,
    objective_fct=DW_v,
    vmap_dim=10**2,
    print_progress=True,
)

Workflow C: Non-SUSY ISD sampling#

What’s new vs §A/§B. Workflows A and B solve the F-term conditions \(D_I W = 0\) (SUSY vacua) using FluxVacuaFinder’s Newton + ISD pipeline. Workflow C targets the broader class of critical points of the scalar potential, \(\partial_I V = 0\), which includes non-SUSY vacua where \(D_I W \neq 0\). The non-SUSY workflow lives directly on FluxVacuaFinder via :func:sample_critical_points; the inner solve targets \(\partial V = 0\) instead of \(DW = 0\).

Background — non-SUSY critical points#

SUSY flux vacua satisfy \(D_I W = 0\) (F-flatness), which implies \(\partial_I V = 0\) automatically. But \(V\) also has critical points where \(\partial_I V = 0\) while \(D_I W \neq 0\) — these are non-SUSY vacua. Two scalar-potential variants are commonly used:

No-scale \(V_{\rm ns} = e^K |DW|^2 \geq 0\): all critical points are either SUSY minima (\(DW = 0\), \(V = 0\)) or non-SUSY critical points (\(DW \neq 0\), \(V > 0\)). Empirically ~85% of the critical points found are minima (not saddles).
Full SUGRA \(V = e^K(|DW|^2 - 3|W|^2)\): can be negative, producing AdS minima and more saddles.

sample_critical_points defaults to the no-scale potential (noscale=True); switch to the full SUGRA potential via noscale=False.

Quick start: finding non-SUSY critical points#

Non-SUSY sampling needs tighter moduli bounds than the §4 canonical sampler — we restrict \(\operatorname{Im}(z)\in(2,5)\) and \(s \in (\sqrt 3/2, 10)\) to stay in a well-controlled region. The ISD-flux completion mode that empirically yields the most non-SUSY points is "ISD-".

sampler_nonSUSY = jvc.data_sampler(
    model,
    moduli_bounds=(2., 5.),
    dilaton_bounds=(np.sqrt(3) / 2, 10.),
    axion_bounds=(-0.5, 0.5),
    seed=42,
)

Nmax = 200  # tadpole bound for sampled fluxes

# moduli_max=100 (default) filters runaway solutions with |z_i| > 100.
finder = model    # post-merge: FluxVacuaFinder is the finder

results = finder.sample_critical_points(
    Nmax=Nmax, noscale=True, sampler=sampler_nonSUSY, n_target=50,
    n_batch=200,
    max_batches=5,
    isd_mode="ISD-",
    solver="scipy",
    verbose=True,
)

Analysing the results#

Each result is a dict with keys moduli, tau, flux, V, |DW|, is_susy, is_minimum, Nflux, etc. The lines below classify the batch into SUSY-vs-non-SUSY and minima-vs-saddles, then print the first ten entries.

n_susy = sum(1 for r in results if r['is_susy'])
n_nonsusy = len(results) - n_susy
n_min = sum(1 for r in results if r['is_minimum'])
n_sad = len(results) - n_min

print(f"Total critical points: {len(results)}")
print(f"  SUSY:     {n_susy} "
      f"({100*n_susy/max(len(results),1):.0f}%)")
print(f"  non-SUSY: {n_nonsusy} "
      f"({100*n_nonsusy/max(len(results),1):.0f}%)")
print(f"  minima:   {n_min} "
      f"({100*n_min/max(len(results),1):.0f}%)")
print(f"  saddle:   {n_sad} "
      f"({100*n_sad/max(len(results),1):.0f}%)")

if results:
    print(f"\n{'#':>3} {'V':>12} {'|DW|':>12} {'min_eig':>10} "
          f"{'type':>12} {'Nflux':>6}")
    print("-" * 60)
    for i, r in enumerate(results[:10]):
        label = "SUSY" if r['is_susy'] else "non-SUSY"
        kind = "minimum" if r['is_minimum'] else "saddle"
        print(f"{i:>3} {r['V']:>12.3e} {r['|DW|']:>12.3e} "
              f"{r.get('min_eig', float('nan')):>10.2e} "
              f"{label+'/'+kind:>12} "
              f"{r.get('Nflux', 0):>6.0f}")

Higher-\(h^{1,2}\) example: non-SUSY workflow at \(h^{1,2}=3\)#

As in §A, the workflow generalises without code changes. For higher \(h^{1,2}\), runaway solutions (\(|z_i|\) or \(|\tau|\) in the \(10^{8}-10^{30}\) range) become common; the moduli_max parameter filters them efficiently and Newton’s early-escape short-circuits iterations as soon as the trajectory leaves the box (≈4× speed-up). For this \(h^{1,2}=3\) model (model_ID=4), physical solutions sit around \(|z|\sim 200-500\), so we widen moduli_max=500.

# h12=3 example (model_ID=4 has hyperplanes + solutions near the sampling region).
model_3_nonSUSY = jvc.FluxVacuaFinder(
    h12=3, model_ID=4, model_type="KS",
    maximum_degree=2, limit="LCS", map_to_fd=True,
)
sampler_3_nonSUSY = jvc.data_sampler(
    model_3_nonSUSY,
    moduli_bounds=(2., 8.),
    dilaton_bounds=(np.sqrt(3) / 2, 10.),
    axion_bounds=(-0.5, 0.5),
    seed=42,
)

finder_3 = model    # post-merge: FluxVacuaFinder is the finder

results_3 = finder_3.sample_critical_points(
    Nmax=50, noscale=True, moduli_max=500, sampler=sampler_3_nonSUSY, n_target=20,
    n_batch=500,
    max_batches=2,
    isd_mode="ISD-",
    solver="newton",
    step_size=0.5,
    newton_max_iters=500,
    verbose=True,
)

n3 = len(results_3)
n3_ns = sum(1 for r in results_3 if not r['is_susy'])
n3_min = sum(1 for r in results_3 if r['is_minimum'])
print(f"\nh12=3 non-SUSY results: {n3} critical points "
      f"({n3_ns} non-SUSY, {n3_min} minima)")

Method comparison and benchmarks#

This section benchmarks the methods exposed in §A/§B against each other (and against bounded_fluxes.enumerate_fluxes / sample_bounded_fluxes from NB08) on the §4 canonical fixture. We use a tighter region \(\operatorname{Im}(z)\in(2,5)\) and \(s\in(\sqrt 3/2, 10)\) for a fair head-to-head.

Bounding-box scaling with \(N_{\max}\)#

The eigenvalue-based bounding box from arXiv:2501.03984 defines the search region for NS-NS fluxes. The Cartesian product over the full \(h\)-vector grows as \(\sim N_{\max}^3\), making systematic enumeration infeasible beyond \(N_{\max} \sim 30\).

# Show how bounding box scales with Nmax
dil_bounds = (math.sqrt(3)/2, 10.)
moduli_bounds = (2., 5.)
axion_bounds = (-0.5, 0.5)

Nmax_values = [5, 10, 20, 50, 100, 200, 500]
box_sizes = []

for Nm in Nmax_values:
    s = jvc.data_sampler(model, moduli_bounds=moduli_bounds,
        dilaton_bounds=dil_bounds, axion_bounds=axion_bounds, seed=42)
    bf = bounded_fluxes(model, sampler=s, Nmax=Nm)
    bf.compute_eigenvalue_bounds(200, verbose=False)
    dim = bf.dimension_H3
    h1m = int(np.ceil(bf._h1_box))
    h2m = int(np.ceil(bf._h2_box))
    n_box = (2*h1m+1)**dim * (2*h2m+1)**dim
    box_sizes.append(n_box)
    feasible = "✓ enumerate" if n_box < 1e8 else ("⚠ streaming" if n_box < 1e10 else "✗ sample only")
    print(f"Nmax={Nm:>3}: h1_max={h1m:>3}, h2_max={h2m:>2}, "
          f"box={n_box:>18,}  {feasible}")

fig, ax = plt.subplots(figsize=(8, 4),dpi=150)
ax.semilogy(Nmax_values, box_sizes, 'o-', color='steelblue', linewidth=2)
ax.axhline(1e8, color='green', linestyle='--', alpha=0.5, label='enumerate feasible (<10⁸)')
ax.axhline(1e10, color='orange', linestyle='--', alpha=0.5, label='streaming feasible (<10¹⁰)')
ax.set_xlabel(r'$N_{\max}$')
ax.set_ylabel('Bounding box size')
ax.set_title('Search space growth with tadpole bound')
ax.legend()
plt.tight_layout()
plt.show()

The four sampling methods#

Method	Approach	Completeness	Best regime
Manual ISD	Random \((z, \tau, h)\) → ISD complete \(f\) → Newton	None	Baseline only
`enumerate_fluxes`	Systematic: all \(h\) in bounding box, stream + filter	Complete	\(N_{\max} \lesssim 30\)
`sample_bounded_fluxes` (Gaussian)	\(h \sim \mathcal{N}(0, \sigma^2 M)\) → ISD → Newton	Stochastic	\(N_{\max} \gtrsim 50\)
`sample_bounded_fluxes` (Gaussian + LS)	Same prior + `linearised_shifts` (moves moduli)	Stochastic	Medium \(N_{\max}\)

Why the Gaussian prior \(h \sim \mathcal{N}(0, \sigma^2 M)\)?#

The continuous ISD tadpole is \(N_{\rm flux}^{\rm cont} = s \cdot h^T M^{-1} h\). For \(h\) drawn from \(\mathcal{N}(0, \sigma^2 M)\):

\[\mathbb{E}\bigl[s_{\min} \cdot h^T M^{-1} h\bigr] = s_{\min} \cdot \sigma^2 \cdot \mathrm{tr}(I) = s_{\min} \cdot \sigma^2 \cdot d\]

Choosing \(\sigma^2 = N_{\max} / (d \cdot s_{\min})\) ensures most samples satisfy \(N_{\rm flux} \lesssim N_{\max}\). The \(M\)-weighting concentrates samples along eigendirections where \(M^{-1}\) has small eigenvalues (low tadpole cost), matching the physical vacuum distribution.

Benchmark helper#

We define a helper function that runs each method with consistent parameters and collects results.

def run_benchmark(Nmax, moduli_bounds=(2., 5.), dil_bounds=(math.sqrt(3)/2, 10.),
                  axion_bounds=(-0.5, 0.5), n_manual=5000, max_batches_sample=5,
                  n_batch_sample=10_000, skip_enumerate=False):
    """Run all four methods and return a results dict."""
    
    sampler = jvc.data_sampler(model, moduli_bounds=moduli_bounds,
        dilaton_bounds=dil_bounds, axion_bounds=axion_bounds, seed=42)
    results = {}
    
    # ── Method 1: Manual ISD ──
    rng = np.random.default_rng(42)
    t0 = time.perf_counter()
    mod_pts = jnp.array(sampler.get_complex_moduli(n_manual), dtype=complex)
    tau_pts = jnp.array(sampler.get_complex_tau(n_manual), dtype=complex)
    h_all = rng.integers(-5, 6, (n_manual, n_fl))
    
    good = []
    for i in range(n_manual):
        h_i = jnp.array(h_all[i], dtype=float)
        if jnp.all(h_i == 0): continue
        try:
            fl = sampler.ISD_sampling(mod_pts[i], jnp.conj(mod_pts[i]),
                                       tau_pts[i], jnp.conj(tau_pts[i]),
                                       h_i, mode="H")
            if fl is not None and jnp.all(jnp.isfinite(fl)):
                fl_r = jnp.array(fl).real
                fl_r = fl_r.at[:n_fl].set(jnp.round(fl_r[:n_fl]))
                tad = abs(float(jnp.real(model.tadpole(fl_r))))
                if 0 < tad <= Nmax:
                    good.append((mod_pts[i], tau_pts[i], fl_r))
        except: pass
    
    n_conv = 0
    if good:
        bf_tmp = bounded_fluxes(model, sampler=sampler, Nmax=Nmax)
        m_arr = jnp.array([x[0] for x in good], dtype=complex)
        t_arr = jnp.array([x[1] for x in good], dtype=complex)
        f_arr = jnp.array([x[2] for x in good])
        m_o, t_o, r_o = bf_tmp.newton_refine_batch(m_arr, t_arr, f_arr,
            step_size=1.0, tol=1e-10, max_iters=100)
        r_np = np.asarray(r_o)
        in_p = np.asarray(bf_tmp.in_patch_batch(m_o, t_o))
        n_conv = int(np.sum((r_np < 1e-10) & in_p))
    dt = time.perf_counter() - t0
    results['Manual ISD'] = {'vacua': n_conv, 'time': dt,
                              'isd_yield': f"{len(good)}/{n_manual}"}
    
    # ── Method 2: enumerate_fluxes ──
    if not skip_enumerate:
        bf_e = bounded_fluxes(model, sampler=sampler, Nmax=Nmax)
        t0 = time.perf_counter()
        res_e = bf_e.enumerate_fluxes(
            n_sample=200, n_isd_per_h=10, max_h_candidates=5_000_000,
            refine=True, return_moduli=True, newton_tol=1e-10,
            confirm_streaming=False, verbose=False, n_moduli_batches=1)
        dt = time.perf_counter() - t0
        results['enumerate_fluxes'] = {'vacua': len(res_e), 'time': dt}
    
    # ── Method 3: sample_bounded (Gaussian, no LS) ──
    bf_s = bounded_fluxes(model, sampler=sampler, Nmax=Nmax)
    t0 = time.perf_counter()
    res_s = bf_s.sample_bounded_fluxes(
        n_target=500, n_batch=n_batch_sample, n_sample=100, n_mod=10,
        max_batches=max_batches_sample, refine=True, return_moduli=True,
        newton_tol=1e-10, verbose=False)
    dt = time.perf_counter() - t0
    results['sample (Gaussian)'] = {'vacua': len(res_s), 'time': dt}
    
    # ── Method 4: sample_bounded (Gaussian + linearised_shifts) ──
    bf_ls = bounded_fluxes(model, sampler=sampler, Nmax=Nmax)
    t0 = time.perf_counter()
    res_ls = bf_ls.sample_bounded_fluxes(
        n_target=500, n_batch=min(n_batch_sample, 5000), n_sample=100, n_mod=10,
        max_batches=max_batches_sample, refine=True, return_moduli=True,
        newton_tol=1e-10, verbose=False,
        use_linearised_shifts=True, n_isd_iters=5, n_moduli_batches=1)
    dt = time.perf_counter() - t0
    results['sample (Gauss+LS)'] = {'vacua': len(res_ls), 'time': dt}
    
    return results


def print_results(results, Nmax):
    """Pretty-print benchmark results."""
    print(f"\n{'='*65}")
    print(f"  Nmax = {Nmax}")
    print(f"{'='*65}")
    print(f"{'Method':<28} {'Vacua':>7} {'Time':>9} {'Rate':>10}")
    print(f"{'-'*65}")
    for name, r in results.items():
        n = r['vacua']; dt = r['time']
        rate = n / max(dt, 0.1)
        if dt < 120: ts = f"{dt:.1f}s"
        elif dt < 3600: ts = f"{dt/60:.1f}m"
        else: ts = f"{dt/3600:.1f}h"
        extra = f"  (ISD: {r['isd_yield']})" if 'isd_yield' in r else ""
        print(f"{name:<28} {n:>7} {ts:>9} {rate:>9.2f}/s{extra}")
    print(f"{'-'*65}")

Benchmark: \(N_{\max} = 10\) (small tadpole)#

At small \(N_{\max}\), the bounding box is manageable and enumerate_fluxes can systematically check all candidates. The tadpole pre-filter (continuous quadratic form \(s_{\min} h^T M^{-1} h \leq N_{\max}\)) eliminates >99% of the box, making enumeration fast.

All reported vacua are Newton-refined to \(\sum_I |D_I W| < 10^{-10}\) and verified to lie inside the sampler’s moduli patch.

res_10 = run_benchmark(Nmax=5, n_manual=5000, max_batches_sample=5, n_batch_sample=10_000)
print_results(res_10, 10)

Benchmark: \(N_{\max} = 200\) (large tadpole)#

At large \(N_{\max}\), the bounding box is too large for enumeration (~\(10^{13}\) candidates). Stochastic sampling with the Gaussian prior becomes the method of choice.

res_200 = run_benchmark(Nmax=200, skip_enumerate=True,
                        n_manual=5000, max_batches_sample=5, n_batch_sample=10_000)
print_results(res_200, 200)

Gaussian prior analysis#

The Gaussian prior \(h \sim \mathcal{N}(0, \sigma^2 M)\) dramatically outperforms uniform sampling because it concentrates samples where the tadpole is small. Let’s visualise this.

# Compare ||h|| distributions: Gaussian prior vs Dataset B
Nmax_plot = 10
s_min_plot = math.sqrt(3) / 2
sampler_plot = jvc.data_sampler(model, moduli_bounds=(2.,5.),
    dilaton_bounds=(s_min_plot, 10.), axion_bounds=(-0.5, 0.5), seed=42)

# Representative ISD matrix
z_rep = jnp.array(sampler_plot.get_complex_moduli(1)[0], dtype=complex)
M_rep = np.asarray(model.ISD_matrix(z_rep, jnp.conj(z_rep)).real)
M_inv_rep = np.linalg.inv(M_rep)

# Gaussian samples
eigvals_M, eigvecs_M = np.linalg.eigh(M_rep)
sigma_sq = Nmax_plot / (n_fl * s_min_plot)
scales = np.sqrt(sigma_sq * np.abs(eigvals_M))
rng_plot = np.random.default_rng(0)
N_samp = 50_000
z_std = rng_plot.standard_normal((N_samp, n_fl))
h_gauss = np.round((z_std * scales[None, :]) @ eigvecs_M.T).astype(int)
h_gauss_norms = np.sqrt(np.sum(h_gauss**2, axis=1))
h_gauss = h_gauss[h_gauss_norms > 0]
h_gauss_norms = h_gauss_norms[h_gauss_norms > 0]

# Tadpole of Gaussian samples
hMh_gauss = np.abs(s_min_plot * np.einsum('hi,ij,hj->h', h_gauss.astype(float), M_inv_rep, h_gauss.astype(float)))

# Uniform-from-box samples (for comparison)
bf_plot = bounded_fluxes(model, sampler=sampler_plot, Nmax=Nmax_plot)
bf_plot.compute_eigenvalue_bounds(200, verbose=False)
h_unif = bf_plot._sample_h_vectors(min(N_samp, 50_000), rng=np.random.default_rng(0))
h_unif_norms = np.sqrt(np.sum(h_unif**2, axis=1))
hMh_unif = np.abs(s_min_plot * np.einsum('hi,ij,hj->h', h_unif.astype(float), M_inv_rep, h_unif.astype(float)))

# Dataset B (if available)
try:
    datsB = jvc.util.load_zipped_pickle("../../../notebooks/dataset_B.p")
    h_B_norms = np.sqrt(np.sum(datsB[:, 12:18]**2, axis=1))
    has_datsB = True
except:
    has_datsB = False

fig, axes = plt.subplots(1, 3, figsize=(16, 4.5))

# Panel 1: ||h|| distribution
bins_h = np.arange(0, 60, 1)
axes[0].hist(h_unif_norms, bins=bins_h, alpha=0.5, density=True, label='Uniform (box)', color='grey')
axes[0].hist(h_gauss_norms, bins=bins_h, alpha=0.6, density=True, label=r'Gaussian $\mathcal{N}(0,\sigma^2 M)$', color='coral')
if has_datsB:
    axes[0].hist(h_B_norms, bins=bins_h, alpha=0.5, density=True, label='Dataset B', color='steelblue')
axes[0].set_xlabel(r'$\|h\|$')
axes[0].set_ylabel('density')
axes[0].set_title(r'$\|h\|$ distribution')
axes[0].legend(fontsize=9)
axes[0].set_xlim(0, 50)

# Panel 2: Continuous tadpole distribution
bins_t = np.linspace(0, 100, 50)
axes[1].hist(hMh_unif, bins=bins_t, alpha=0.5, density=True, label='Uniform', color='grey')
axes[1].hist(hMh_gauss, bins=bins_t, alpha=0.6, density=True, label='Gaussian', color='coral')
axes[1].axvline(Nmax_plot, color='red', linestyle='--', linewidth=2, label=f'$N_{{\\max}}={Nmax_plot}$')
axes[1].set_xlabel(r'$s_{\min} \cdot h^T M^{-1} h$')
axes[1].set_ylabel('density')
axes[1].set_title('Continuous tadpole')
axes[1].legend(fontsize=9)
axes[1].set_xlim(0, 100)

# Panel 3: Tadpole pass rate vs Nmax
Nmax_range = np.arange(1, 201)
pass_gauss = [np.sum(hMh_gauss <= nm) / len(hMh_gauss) * 100 for nm in Nmax_range]
pass_unif = [np.sum(hMh_unif <= nm) / len(hMh_unif) * 100 for nm in Nmax_range]
axes[2].plot(Nmax_range, pass_gauss, color='coral', linewidth=2, label='Gaussian')
axes[2].plot(Nmax_range, pass_unif, color='grey', linewidth=2, label='Uniform (box)')
axes[2].set_xlabel(r'$N_{\max}$')
axes[2].set_ylabel('% passing tadpole')
axes[2].set_title('Tadpole acceptance rate')
axes[2].legend()

plt.suptitle(f'Gaussian prior vs uniform box sampling (Nmax={Nmax_plot})', y=1.02, fontsize=13)
plt.tight_layout()
plt.show()

# Print key numbers
print(f"Tadpole acceptance at Nmax={Nmax_plot}:")
print(f"  Gaussian: {np.sum(hMh_gauss <= Nmax_plot)/len(hMh_gauss)*100:.1f}%")
print(f"  Uniform:  {np.sum(hMh_unif <= Nmax_plot)/len(hMh_unif)*100:.1f}%")
print(f"  → Gaussian is {np.sum(hMh_gauss<=Nmax_plot)/max(np.sum(hMh_unif<=Nmax_plot),1):.0f}× more efficient")

Legacy: uniform sampling of fluxes from a box#

Before ISD-based sampling, the common practice was to draw integer fluxes uniformly from a hypercube and Newton-refine every survivor. We include this here as a baseline so the speed-up from the four ISD-based methods above is concrete.

# Legacy uniform-box sampling for comparison.
kwargs_box = {"Q": model.D3_tadpole, "return_flag": True,
              "constraints": None, "remove_NANs": True,
              "in_axes": (0, 0, 0)}
linearised_shifts_box = vmapping_func(
    model.linearised_shifts, mode="box", **kwargs_box)

t0 = time.perf_counter()
moduli_box, tau_box, fluxes_box, residuals_box = (
    model.sample_SUSY_flux_vacua(
        N=2 * 10**2,
        rns_key=rns_key,
        sampler=sampler,
        optimiser=linearised_shifts_box,
        objective_fct=DW_v,
        vmap_dim=10**2,
        print_progress=True,
    )
)
dt = time.perf_counter() - t0
print(f"Legacy uniform-box: {len(moduli_box)} vacua in {dt:.1f}s")

ISD mode comparison#

The sample_critical_points supports four ISD sampling modes for generating flux candidates:

Mode	Input	Completes via	Best for
F	random \(f\)	\(h = M^{-1}\Sigma(f - c_0 h)\)	Mixed SUSY + non-SUSY (85% minima rate)
H	random \(h\)	\(f = sM\Sigma h + c_0 h\)	SUSY-dominated
ISD+	random \([f_2 \mid h_2]\)	compute \([f_1 \mid h_1]\) via \(\mathcal{N}\)	Almost all SUSY
ISD-	random \([f_1 \mid h_1]\)	compute \([f_2 \mid h_2]\) via \(\mathcal{N}^{-1}\)	Most non-SUSY (78% rate)

Let’s compare all four modes on 500 input candidates each.

modes = ["F", "H", "ISD+", "ISD-"]
mode_results = {}

for mode in modes:
    print(f"\n{'='*60}")
    print(f"  Mode: {mode}")
    print(f"{'='*60}")
    finder_m = model    # post-merge: FluxVacuaFinder is the finder
    res = finder_m.sample_critical_points(
Nmax=Nmax, noscale=True, sampler=sampler,         n_target=500,       # large target so we don't stop early
        n_batch=500,
        max_batches=1,      # single batch of 500 candidates
        isd_mode=mode,
        solver="scipy",
        verbose=True,
    )
    mode_results[mode] = res

# Summary table
print(f"{'Mode':<8} {'Found':>6} {'SUSY':>6} {'non-SUSY':>9} {'minima':>7} {'saddle':>7} {'min rate':>9} {'non-SUSY%':>10}")
print("-" * 72)
for mode in modes:
    res = mode_results[mode]
    n = len(res)
    n_s = sum(1 for r in res if r.get('is_susy', False))
    n_ns = n - n_s
    n_m = sum(1 for r in res if r.get('is_minimum', False))
    n_sad = n - n_m
    mr = f"{100*n_m/n:.0f}%" if n > 0 else "—"
    nsr = f"{100*n_ns/n:.0f}%" if n > 0 else "—"
    print(f"{mode:<8} {n:>6} {n_s:>6} {n_ns:>9} {n_m:>7} {n_sad:>7} {mr:>9} {nsr:>10}")

Solver comparison: Newton vs scipy#

The sample_critical_points supports two main solver backends:

Newton (solver="newton"): JIT-compiled Newton’s method using model.dV_x and model.ddV_x. Converges quadratically when close to a solution.
scipy (solver="scipy"): Uses scipy.optimize.root with the hybrid Powell method. Generally faster per-candidate due to adaptive step control.

Let’s compare them on the same ISD- flux candidates.

solvers = ["newton", "scipy"]
solver_results = {}

for slv in solvers:
    print(f"\n{'='*60}")
    print(f"  Solver: {slv}")
    print(f"{'='*60}")
    finder_s = model    # post-merge: FluxVacuaFinder is the finder
    t_start = time.perf_counter()
    res = finder_s.sample_critical_points(
Nmax=Nmax, noscale=True, sampler=sampler,         n_target=500,
        n_batch=500,
        max_batches=1,
        isd_mode="ISD-",
        solver=slv,
        verbose=True,
    )
    t_end = time.perf_counter()
    solver_results[slv] = {'results': res, 'time': t_end - t_start}

# Comparison
print(f"\n{'Solver':<10} {'Found':>6} {'SUSY':>6} {'non-SUSY':>9} {'minima':>7} {'time (s)':>9}")
print("-" * 52)
for slv in solvers:
    res = solver_results[slv]['results']
    n = len(res)
    n_s = sum(1 for r in res if r.get('is_susy', False))
    n_ns = n - n_s
    n_m = sum(1 for r in res if r.get('is_minimum', False))
    t = solver_results[slv]['time']
    print(f"{slv:<10} {n:>6} {n_s:>6} {n_ns:>9} {n_m:>7} {t:>8.1f}s")

Advanced: performance tuning#

This section addresses performance and quality tuning for non-SUSY sampling at production scale. We compare the vectorised Adam and the hybrid solver against scipy/Newton, then walk through the seven ablation knobs from NB19 §10 (moduli_max, sampling region, ISD mode, solution magnitudes, \(N_{\max}\), and a production-tuned run).

Vectorised optax solver and hybrid approach#

The per-candidate loop in the Newton and scipy solvers becomes a bottleneck for large batches. sample_critical_points includes a vectorised optax Adam solver (solver="adam_v") that uses lax.scan for the step loop and vmap over the batch — compiled into a single XLA kernel.

Even better is the hybrid solver (solver="hybrid"): it runs a short vectorised Adam warm-start to move all candidates into Newton’s convergence basin, then finishes with Newton’s quadratic convergence. This typically finds 40% more critical points than scipy alone.

Parameter analysis#

The key hyperparameters for the vectorised Adam solver are:

Learning rate schedule: cosine annealing with init_lr=0.5 works best
Gradient clipping: clip_by_global_norm(10.0) prevents divergence
Step count: 1000 is the optimal warm-start — more steps have diminishing returns
Objective function: |dV|² outperforms log(|dV|²+ε) and V directly

import optax

# Generate a shared set of ISD- candidates for fair comparison
finder_bench = model    # post-merge: FluxVacuaFinder is the finder
mod_pts = jnp.array(sampler.get_complex_moduli(500), dtype=complex)
tau_pts = jnp.array(sampler.get_complex_tau(500), dtype=complex)
x0_all, fl_all, _ = finder_bench._generate_flux_candidates(
    500, mod_pts, tau_pts, isd_mode="ISD-")
n_bench = min(200, len(x0_all))
print(f"Benchmarking on {n_bench} ISD- candidates\n")

tol = 1e-8

# --- Parameter sweep for Adam_v ---
print("Learning rate schedule comparison (Adam, cosine annealing, 5000 steps):")
print(f"{'lr':>6} {'conv':>6} {'min_res':>10} {'p10_res':>10} {'med_res':>10} {'time':>7}")
print("-" * 55)
for lr in [0.1, 0.3, 0.5, 0.7, 1.0]:
    sched = optax.cosine_decay_schedule(init_value=lr, decay_steps=5000)
    opt = optax.chain(optax.clip_by_global_norm(10.0), optax.adam(learning_rate=sched))
    t0 = time.perf_counter()
    _, res, conv = finder_bench._solve_dV_optax_batch(
        x0_all[:n_bench], fl_all[:n_bench], optimiser=opt,
        n_steps=5000, objective="dV2", tol=tol, sub_batch_size=50)
    dt = time.perf_counter() - t0
    sr = np.sort(res)
    print(f"{lr:>6.1f} {int(conv.sum()):>6} {sr[0]:>10.2e} {sr[int(0.1*len(sr))]:>10.2e} "
          f"{np.median(res):>10.2e} {dt:>6.1f}s")

# --- Objective function comparison ---
print("Objective function comparison (Adam cosine lr=0.5, 5000 steps):")
print(f"{'objective':>10} {'conv':>6} {'min_res':>10} {'med_res':>10}")
print("-" * 40)
for obj in ["dV2", "log_dV2", "V"]:
    sched = optax.cosine_decay_schedule(init_value=0.5, decay_steps=5000)
    opt = optax.chain(optax.clip_by_global_norm(10.0), optax.adam(learning_rate=sched))
    _, res, conv = finder_bench._solve_dV_optax_batch(
        x0_all[:n_bench], fl_all[:n_bench], optimiser=opt,
        n_steps=5000, objective=obj, tol=tol, sub_batch_size=50)
    print(f"{obj:>10} {int(conv.sum()):>6} {np.min(res):>10.2e} {np.median(res):>10.2e}")

# --- Full solver comparison: scipy vs Newton vs Adam_v vs Hybrid ---
from scipy.optimize import root as scipy_root

solver_bench = {}

# 1. scipy
t0 = time.perf_counter()
scipy_conv = 0
for j in range(n_bench):
    fl_np = fl_all[j].copy()
    def f(x, _fl=fl_np): return np.asarray(
        model.dV_x(jnp.array(x), jnp.array(_fl), noscale=True))
    def jac(x, _fl=fl_np): return np.asarray(
        model.ddV_x(jnp.array(x), jnp.array(_fl), noscale=True))
    r = scipy_root(f, x0=x0_all[j], method='hybr', jac=jac)
    if r.success and max(abs(r.fun)) < tol: scipy_conv += 1
solver_bench['scipy'] = (scipy_conv, time.perf_counter() - t0)

# 2. Newton
t0 = time.perf_counter()
newton_conv = 0
for j in range(n_bench):
    _, _, c = finder_bench._solve_dV_newton_single(
        x0_all[j], fl_all[j], tol=tol, max_iters=300)
    if c: newton_conv += 1
solver_bench['Newton'] = (newton_conv, time.perf_counter() - t0)

# 3. Hybrid (Adam 1k + Newton)
sched = optax.cosine_decay_schedule(init_value=0.5, decay_steps=1000)
opt = optax.chain(optax.clip_by_global_norm(10.0), optax.adam(learning_rate=sched))
t0 = time.perf_counter()
x_warm, res_warm, _ = finder_bench._solve_dV_optax_batch(
    x0_all[:n_bench], fl_all[:n_bench], optimiser=opt,
    n_steps=1000, objective="dV2", tol=tol, sub_batch_size=50)
hybrid_conv = int(np.sum(res_warm < tol))
mask = (res_warm < 1.0) & (res_warm >= tol)
for j in np.where(mask)[0]:
    _, res_n, conv_n = finder_bench._solve_dV_newton_single(
        x_warm[j], fl_all[j], tol=tol, max_iters=300)
    if conv_n: hybrid_conv += 1
solver_bench['Hybrid 1k'] = (hybrid_conv, time.perf_counter() - t0)

# Print comparison
print(f"\n{'Solver':<14} {'Converged':>10} {'Rate':>8} {'Time':>8}")
print("-" * 44)
for name in ['scipy', 'Newton', 'Hybrid 1k']:
    nc, t = solver_bench[name]
    print(f"{name:<14} {nc:>10} {100*nc/n_bench:>7.1f}% {t:>7.1f}s")

Using the hybrid solver via `sample_critical_points`#

The hybrid solver is the recommended choice when you want to maximise the number of critical points found:

# The hybrid solver: vectorised Adam warm-start (1k steps) + Newton refinement
finder_hybrid = model    # post-merge: FluxVacuaFinder is the finder

results_hybrid = finder_hybrid.sample_critical_points(
Nmax=Nmax, noscale=True, sampler=sampler,     n_target=50,
    n_batch=200,
    max_batches=3,
    isd_mode="ISD-",
    solver="hybrid",
    optax_steps=1000,           # Adam warm-start steps (1k is optimal)
    newton_max_iters=300,       # Newton refinement iterations
    verbose=True,
)

n_h = len(results_hybrid)
n_ns_h = sum(1 for r in results_hybrid if not r.get('is_susy', False))
n_min_h = sum(1 for r in results_hybrid if r.get('is_minimum', False))
print(f"\nHybrid results: {n_h} critical points "
      f"({n_ns_h} non-SUSY, {n_min_h} minima)")

Performance analysis: moduli_max, sampling regions, and Nmax#

The moduli_max parameter, the moduli/dilaton sampling region, and \(N_{\max}\) all affect how many physical critical points are found vs. how many are runaways. This section benchmarks their interplay.

10a. Effect of `moduli_max`#

The moduli_max filter rejects solutions where any real coordinate component exceeds the threshold. Setting it too tight loses physical solutions; too loose lets runaways through.

# 10a. moduli_max sweep (ISD-, scipy, default region)
print(f"{'moduli_max':>12} {'CritPts':>8} {'SUSY':>5} {'nonSUSY':>8} "
      f"{'minima':>7} {'saddle':>7} {'Time':>7}")
print("-" * 62)

sampler_default = jvc.data_sampler(model, moduli_bounds=(2., 5.),
    dilaton_bounds=(math.sqrt(3)/2, 10.), axion_bounds=(-0.5, 0.5), seed=42)

for mmax in [10, 20, 50, 100, 200, 500, None]:
    f = model    # post-merge: FluxVacuaFinder is the finder
    t0 = time.perf_counter()
    res = f.sample_critical_points(
        Nmax=Nmax, noscale=True, moduli_max=mmax, sampler=sampler_default, n_target=500, n_batch=500, max_batches=1,
        isd_mode="ISD-", solver="scipy", verbose=False)
    dt = time.perf_counter() - t0
    n_s = sum(1 for r in res if r.get('is_susy', False))
    n_m = sum(1 for r in res if r.get('is_minimum', False))
    label = str(mmax) if mmax is not None else "None"
    print(f"{label:>12} {len(res):>8} {n_s:>5} {len(res)-n_s:>8} "
          f"{n_m:>7} {len(res)-n_m:>7} {dt:>6.1f}s")

Effect of sampling region#

The moduli/dilaton sampling bounds determine where starting points are drawn. Tighter regions (closer to the origin) produce more physical solutions because the solver has less distance to travel.

# 10b. Sampling region sweep (ISD-, scipy, moduli_max=100)
regions = [
    ((1., 2.), (math.sqrt(3)/2, 3.), "very tight: Im(z)∈(1,2), s∈(√3/2,3)"),
    ((1., 3.), (math.sqrt(3)/2, 5.), "tight: Im(z)∈(1,3), s∈(√3/2,5)"),
    ((2., 5.), (math.sqrt(3)/2, 10.), "default: Im(z)∈(2,5), s∈(√3/2,10)"),
    ((1., 10.), (math.sqrt(3)/2, 10.), "wide: Im(z)∈(1,10), s∈(√3/2,10)"),
    ((3., 8.), (1., 15.), "moderate: Im(z)∈(3,8), s∈(1,15)"),
    ((5., 15.), (2., 20.), "deep LCS: Im(z)∈(5,15), s∈(2,20)"),
]

print(f"{'Region':>48} {'CritPts':>8} {'SUSY':>5} {'nonSUSY':>8} "
      f"{'minima':>7} {'Time':>7}")
print("-" * 90)

for mod_b, dil_b, label in regions:
    s = jvc.data_sampler(model, moduli_bounds=mod_b,
        dilaton_bounds=dil_b, axion_bounds=(-0.5, 0.5), seed=42)
    f = model    # post-merge: FluxVacuaFinder is the finder
    t0 = time.perf_counter()
    res = f.sample_critical_points(
        Nmax=Nmax, noscale=True, moduli_max=100, sampler=s, n_target=500, n_batch=500, max_batches=1,
        isd_mode="ISD-", solver="scipy", verbose=False)
    dt = time.perf_counter() - t0
    n_s = sum(1 for r in res if r.get('is_susy', False))
    n_m = sum(1 for r in res if r.get('is_minimum', False))
    print(f"{label:>48} {len(res):>8} {n_s:>5} {len(res)-n_s:>8} "
          f"{n_m:>7} {dt:>6.1f}s")

ISD mode comparison at tight region#

With moduli_max=100 and a tight sampling region, the relative performance of ISD modes changes significantly — modes H and ISD+ now outperform ISD- because they produce solutions closer to the starting points.

# 10c. ISD mode comparison at tight region, moduli_max=100
sampler_tight = jvc.data_sampler(model, moduli_bounds=(1., 2.),
    dilaton_bounds=(math.sqrt(3)/2, 3.), axion_bounds=(-0.5, 0.5), seed=42)

print(f"{'Mode':>6} {'Solver':>8} {'CritPts':>8} {'SUSY':>5} {'nonSUSY':>8} "
      f"{'minima':>7} {'Time':>7}")
print("-" * 58)

for mode in ['ISD-', 'F', 'H', 'ISD+']:
    for solver in ['scipy', 'hybrid']:
        f = model    # post-merge: FluxVacuaFinder is the finder
        t0 = time.perf_counter()
        kw = {'optax_steps': 1000} if solver == 'hybrid' else {}
        res = f.sample_critical_points(
            Nmax=Nmax, noscale=True, moduli_max=100, sampler=sampler_tight, n_target=500, n_batch=500, max_batches=1,
            isd_mode=mode, solver=solver, verbose=False, **kw)
        dt = time.perf_counter() - t0
        n_s = sum(1 for r in res if r.get('is_susy', False))
        n_m = sum(1 for r in res if r.get('is_minimum', False))
        print(f"{mode:>6} {solver:>8} {len(res):>8} {n_s:>5} {len(res)-n_s:>8} "
              f"{n_m:>7} {dt:>6.1f}s")

Solution magnitude distribution#

Without moduli_max, the vast majority of solutions are runaways. The tight region has the highest fraction of physical solutions.

# 10d. Solution magnitude distributions (ISD-, scipy, no moduli_max)
dist_regions = [
    ((1., 2.), (math.sqrt(3)/2, 3.), "very tight"),
    ((2., 5.), (math.sqrt(3)/2, 10.), "default"),
    ((5., 15.), (2., 20.), "deep LCS"),
]

for mod_b, dil_b, label in dist_regions:
    s = jvc.data_sampler(model, moduli_bounds=mod_b,
        dilaton_bounds=dil_b, axion_bounds=(-0.5, 0.5), seed=42)
    f = model    # post-merge: FluxVacuaFinder is the finder
    res = f.sample_critical_points(
        Nmax=Nmax, noscale=True, moduli_max=None, sampler=s, n_target=500, n_batch=500, max_batches=1,
        isd_mode="ISD-", solver="scipy", verbose=False)

    if not res:
        print(f"{label}: 0 solutions")
        continue

    max_abs = np.array([float(np.max(np.abs(np.append(r['moduli'], r['tau']))))
                        for r in res])
    print(f"{label} ({len(res)} pts):")
    print(f"  max(|z|,|tau|): min={max_abs.min():.1f}, "
          f"median={np.median(max_abs):.1e}, max={max_abs.max():.1e}")
    for hi in [5, 10, 50, 100]:
        n_in = int(np.sum(max_abs <= hi))
        print(f"    ≤{hi:>5}: {n_in:>4}/{len(res)} ({100*n_in/len(res):>5.1f}%)")
    print()

\(N_{\max}\) sensitivity#

The tadpole constraint controls flux magnitudes. Larger \(N_{\max}\) allows more diverse fluxes but also more runaways.

# 10e. Nmax sweep (tight region, ISD-, scipy, moduli_max=100)
print(f"{'Nmax':>6} {'CritPts':>8} {'SUSY':>5} {'nonSUSY':>8} "
      f"{'minima':>7} {'Time':>7}")
print("-" * 48)

for nmax in [10, 20, 50, 100, 200, 500]:
    f = model    # post-merge: FluxVacuaFinder is the finder
    t0 = time.perf_counter()
    res = f.sample_critical_points(
        Nmax=nmax, noscale=True, moduli_max=100, sampler=sampler_tight, n_target=500, n_batch=500, max_batches=1,
        isd_mode="ISD-", solver="scipy", verbose=False)
    dt = time.perf_counter() - t0
    n_s = sum(1 for r in res if r.get('is_susy', False))
    n_m = sum(1 for r in res if r.get('is_minimum', False))
    print(f"{nmax:>6} {len(res):>8} {n_s:>5} {len(res)-n_s:>8} "
          f"{n_m:>7} {dt:>6.1f}s")

Production run: tight region + hybrid + moduli_max=100#

Combining the best settings: tight sampling region near the origin, hybrid solver, and moduli_max=100. All solutions are guaranteed to have \(|z_i|, |\tau| \leq 100\).

# 10f. Production run: tight region + hybrid + moduli_max=100
finder_prod = model    # post-merge: FluxVacuaFinder is the finder

results_prod = finder_prod.sample_critical_points(
    Nmax=Nmax, noscale=True, moduli_max=100, sampler=sampler_tight, n_target=200, n_batch=500, max_batches=5,
    isd_mode="ISD-", solver="hybrid", optax_steps=1000, verbose=True)

n_p = len(results_prod)
if n_p > 0:
    max_abs = np.array([float(np.max(np.abs(np.append(r['moduli'], r['tau']))))
                        for r in results_prod])
    Vs = np.array([r['V'] for r in results_prod])
    DWs = np.array([r['|DW|'] for r in results_prod])
    Nfluxes = np.array([r['Nflux'] for r in results_prod])

    print(f"\n{n_p} critical points found:")
    print(f"  max(|z|,|tau|): min={max_abs.min():.1f}, "
          f"median={np.median(max_abs):.1f}, max={max_abs.max():.1f}")
    print(f"  V: min={Vs.min():.2e}, median={np.median(Vs):.2e}")
    print(f"  |DW|: min={DWs.min():.2e}, median={np.median(DWs):.2e}")
    print(f"  Nflux: min={Nfluxes.min():.0f}, median={np.median(Nfluxes):.0f}, "
          f"max={Nfluxes.max():.0f}")

    # All solutions are within moduli_max
    print(f"\n  100% of solutions have max(|z|,|tau|) ≤ {max_abs.max():.1f} "
          f"(moduli_max={finder_prod.moduli_max})")

Summary of performance analysis#

Key findings for this model (h12=2, CP\(^4\)[1,1,1,6,9][18]):

moduli_max effect: Without filtering, ~85% of scipy solutions are runaways (\(|z| > 100\)). Setting moduli_max=100 retains only physical solutions with a modest reduction in count (~30-60 vs ~160 unfiltered). The early Newton escape provides a ~2× speedup.
Sampling region matters more than moduli_max: The very tight region Im(z) \(\in\) (1,2), \(s \in\) (\(\sqrt{3}/2\), 3) produces the most physical solutions — 59 with scipy vs 32 at the default region. Starting closer to the origin helps because the solver is less likely to escape.
ISD mode ranking changes with moduli_max: Without filtering, ISD- dominates (most solutions overall). With filtering, ISD+ (163 pts) and H (117 pts) outperform ISD- (54 pts) at the tight region — they produce solutions closer to the starting points.
\(N_{\max}\) has diminishing returns with moduli_max: Going from \(N_{\max}=50\) to 200 adds only ~5 extra physical solutions (58 vs 53). Most of the additional flux diversity at higher \(N_{\max}\) produces runaways.
Best production config: tight region + solver="hybrid" + moduli_max=100 → ~100 physical critical points per 500 candidates, all with \(|z|, |\tau| \leq 26\).

Take-aways#

Three workflows, one fixture. Workflow A (manual, doubly-vmapped), Workflow B (model.sample_SUSY_flux_vacua wrapper), and Workflow C (FluxVacuaFinder.sample_critical_points) all consume the same §4 model + sampler. The algorithmic differences are the mode flag (SUSY vs non-SUSY), the ISD-completion direction ("F" / "H" / "ISD±"), and whether the inner loop is hand-rolled or wrapped.
Use the wrapper for routine work. §B’s sample_SUSY_flux_vacua(..., optimiser=...) is the shortest path to a batch of refined SUSY vacua. Drop to §A’s doubly-vmapped pipeline only when you need to inspect intermediate ISD shifts or interleave custom logic between rounding and Newton refinement.
For non-SUSY vacua, use FluxVacuaFinder.sample_critical_points with noscale=True and isd_mode="ISD-". This combination empirically yields ~85% minima (not saddles), the highest non-SUSY yield among ISD modes, and no \(V \to -\infty\) runaways from negative-cost SUGRA terms.
Tadpole regime dictates the method. At \(N_{\max} \lesssim 30\) use enumerate_fluxes (complete, fast). At \(N_{\max} \gtrsim 50\) switch to sample_bounded_fluxes with the Gaussian prior \(h \sim \mathcal{N}(0, \sigma^2 M)\) — it concentrates samples along low-cost eigendirections and out-performs uniform-box sampling by ~50×.
Performance tuning has six knobs (§Advanced): moduli_max for runaway filtering, the moduli/dilaton sampling region, the ISD mode, \(N_{\max}\), the solver (scipy for raw speed, hybrid for convergence basin extension, newton for tight quadratic refinement near a vacuum), and n_target/n_batch. Tighter regions + moduli_max=100 + solver="hybrid" is the recommended starting point for production runs.
The workflows scale to higher \(h^{1,2}\) without code changes (see the §A/§C h^{1,2}=3 mini-examples) — only the model, sampler, and moduli_max/Nmax knobs need adjustment.
The legacy uniform-box approach (§6) still works but produces an order of magnitude fewer vacua per unit wall-time. Useful only as a sanity-check baseline.