Flux bounding: enumeration and stochastic sampling

Flux bounding: enumeration and stochastic sampling#

What’s in this notebook? This notebook introduces the bounded_fluxes class from jaxvacua.flux_bounding, which implements the systematic flux-bounding and enumeration algorithm of arXiv:2501.03984. The goal here is to understand the method and the two public workflows — complete enumeration and stochastic bounded sampling — on laptop-safe examples.

In this notebook, you will learn:

How eigenvalue bounds turn an infinite flux search into a finite NS-NS flux box.

How enumerate_fluxes and sample_bounded_fluxes differ.

Which parameters control runtime, memory and feasibility.

How to recognise when a scan should move from a tutorial notebook to an advanced case study.

Non-goals: full literature reproduction and Dataset B recovery now live in NB17 — Flux-bounding case studies. Solver and sampler benchmarks live in NB16 — Sampling benchmarks.

Prerequisites: NB05 — Finding flux vacua, NB06 — ISD principle, and NB07 — Practical ISD sampling.

Recommended paths through this notebook#

Flux bounding contains both a mathematical construction and practical search APIs. For a first pass, read:

Quick start — the smallest complete bounded-search example.
Algorithm overview — why the eigenvalue bounds make the search finite.
Step-by-step walkthrough — how candidates are filtered and refined.
Enumeration vs stochastic sampling — which public entry point to use.

The full mirror-octic and Dataset B reproductions use the same ideas, but they are intentionally separated into NB17 because they can require longer runtimes and more memory.

Setup#

import warnings, time, math
import numpy as np
from tqdm.auto import tqdm
from scipy.optimize import root

import jax
import jax.numpy as jnp
jax.config.update("jax_enable_x64", True)

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

import jaxvacua as jvc
from jaxvacua.flux_bounding import bounded_fluxes

warnings.filterwarnings("ignore")
# Keep the default notebook execution short. Turn these on only for local
# experiments where longer runtime and higher memory use are acceptable.
RUN_LARGE_NMAX_DEMO = False
RUN_NEWTON_REFINEMENT_DEMO = False

Canonical fixture. The main walkthrough uses the degree-18 hypersurface in \(\mathbb{CP}^{1,1,1,6,9}\) (\(h^{2,1}=2\), model_ID=1), the same compact test geometry used throughout the vacuum-finding tutorials. The sampler restricts \(\operatorname{Im}(z) \in [2, 3]\) and \(s = \operatorname{Im}(\tau) \in [2, 5]\) so that the tutorial examples stay small enough for a laptop.

h12 = 2
model = jvc.FluxVacuaFinder(h12=h12, model_ID=1)
print(model)

sampler = jvc.data_sampler(
    model,
    moduli_bounds=(2., 3.),
    axion_bounds=(-0.5, 0.5),
    dilaton_bounds=(2., 5.),
)
print(sampler)

Quick start: the 1D mirror octic#

Before the two-modulus walkthrough below, here is the entire flux-bounding workflow on the simplest geometry — the mirror octic, the one-parameter (\(h^{2,1}=1\)) mirror of the degree-8 hypersurface in \(\mathbb{WP}^4_{1,1,1,1,4}\) (\(\kappa=2\), \(c_2\cdot D=44\), \(\chi=-296\)).

The algorithm in one sentence: from the extremal eigenvalues of the ISD matrix \(\mathcal{M}\) and the gauge-kinetic matrix \(\mathcal{N}\) it builds a finite box of integer NS-NS fluxes \(h\), ISD-completes each to a full flux \([f\,|\,h]\), keeps those passing the tadpole and eigenvalue filters, and (optionally) Newton-refines to SUSY vacua. The call sequence is always the same:

model = FluxVacuaFinder(...)                    # the 4d EFT
bf    = bounded_fluxes(model, sampler=..., Nmax=...)
vacua = bf.enumerate_fluxes(...)                # or bf.sample_bounded_fluxes(...)

# Mirror octic: one-parameter model (h12 = 1), stored as KS model_ID = 1.
m_oct = jvc.FluxVacuaFinder(h12=1, model_ID=1)
print(m_oct)

sampler_oct = jvc.data_sampler(
    m_oct,
    moduli_bounds=(1.5, 7.0),     # Im(z) window in the LCS region
    axion_bounds=(-0.5, 0.5),
    dilaton_bounds=(1.0, 10.0),
    use_jax=True,
    seed=0,
)
bf_oct = bounded_fluxes(m_oct, sampler=sampler_oct, Nmax=4)

# Complete enumeration inside the eigenvalue box (provably complete at small Nmax).
vac_oct = bf_oct.enumerate_fluxes(
    n_sample=200,
    n_isd_per_h=20,
    refine=True,
    verbose=False,
    confirm_streaming=False,
)
print(f"Found {len(vac_oct)} octic SUSY flux vacua with N_flux <= 4.")
for v in vac_oct[:5]:
    nf = float(m_oct.tadpole(jnp.array(v["flux"])).real)
    print(f"  flux={np.round(np.array(v['flux'])).astype(int).tolist()}  "
          f"z={complex(v['moduli'][0]):.3f}  N_flux={nf:.0f}")

Literature cross-checks live in NB17#

The mirror-octic dataset from Plauschinn-Schlechter, arXiv:2310.06040, is now included locally as data_mirror_octic.csv. The full convention map and reproduction check are in NB17 — Flux-bounding case studies. NB08 keeps only the small bounded-search quickstart so the method notebook remains fast.

Algorithm overview#

The flux-bounding algorithm derives rigorous upper bounds on the integer flux quanta \((f, h)\) from the eigenvalues of two period-matrix-dependent objects:

The ISD matrix \(\mathcal{M}(z,\bar{z})\): its largest eigenvalue \(\lambda_{\max}\) controls the size of the full flux vector.
The gauge kinetic matrix \(\mathcal{N}(z,\bar{z})\): the extreme eigenvalues of \(-\mathrm{Im}(\mathcal{N})\) and \(\mathrm{Im}(\mathcal{N}^{-1})\) bound the NSNS-flux sub-vectors \(h_1, h_2\) individually.

Given a sample of moduli points, global extrema of these eigenvalues are computed and used to construct a bounding box for the integer NSNS-flux vectors \(h = (h_1, h_2)\). A complete set of candidate \(h\) vectors is then enumerated inside this box; for each candidate, an ISD-projected RR-flux \(f\) is constructed via

(1)#\[ f \approx \bigl(s\,\mathcal{M}\,\Sigma + c_0\bigr)\,h,\]

and the resulting flux configuration is tested against the D3-tadpole constraint and all local eigenvalue bounds.

Eigenvalue bounds#

Initialising `bounded_fluxes`#

The bounded_fluxes class is the central object. It takes a model, an optional sampler, a D3-tadpole bound Nmax, and an optional lower bound dil_min on \(s = \mathrm{Im}(\tau)\).

In the examples below we pass Nmax explicitly. This is the clearest tutorial style and avoids ambiguity about which tadpole bound is intended. If dil_min is not supplied it defaults to \(\sqrt{3}/2\), the boundary of the \(\mathrm{SL}(2,\mathbb{Z})\) fundamental domain.

bf = bounded_fluxes(model, sampler=sampler, Nmax=50)
print(bf)
print(f"n_fluxes     = {bf.n_fluxes}   (= 2*(h12+1))")
print(f"dimension_H3 = {bf.dimension_H3}   (= h12+1)")
print(f"Nmax         = {bf.Nmax}")
print(f"dil_min      = {bf.dil_min:.4f}  (= sqrt(3)/2)")

Eigenvalue quantities at a single moduli point#

The method compute_evs(moduli) computes the five key eigenvalue quantities at a single point in moduli space:

Quantity	Definition
\(\lambda_{\max}\)	Largest eigenvalue of \(\mathcal{M}\) (ISD matrix)
\(\mu_{\min}\)	Smallest eigenvalue of \(-\mathrm{Im}(\mathcal{N})\)
\(\mu_{\max}\)	Largest eigenvalue of \(-\mathrm{Im}(\mathcal{N})\)
\(\tilde\mu_{\min}\)	Smallest eigenvalue of \(\mathrm{Im}(\mathcal{N}^{-1})\)
\(\tilde\mu_{\max}\)	Largest eigenvalue of \(\mathrm{Im}(\mathcal{N}^{-1})\)

This method is JIT-compiled and fast to evaluate repeatedly:

# A test point inside the LCS region
moduli_test = jnp.array([0.1 + 2.5j, -0.1 + 3.0j])

lam_max, mu_min, mu_max, tmu_min, tmu_max = bf.compute_evs(moduli_test)

print(f"lambda_max   = {lam_max:.6f}")
print(f"mu_min       = {mu_min:.6f}")
print(f"mu_max       = {mu_max:.6f}")
print(f"tmu_min      = {tmu_min:.6f}")
print(f"tmu_max      = {tmu_max:.6f}")

The batched version compute_evs_vmap(moduli_batch) evaluates these quantities over a whole array of moduli in one vectorised call:

moduli_batch = sampler.get_complex_moduli(200)
evs_batch = bf.compute_evs_vmap(jnp.array(moduli_batch))
lam_arr, mu_min_arr, mu_max_arr, tmu_min_arr, tmu_max_arr = evs_batch

print(f"lambda_max range : [{float(lam_arr.min()):.4f}, {float(lam_arr.max()):.4f}]")
print(f"mu_min    range  : [{float(mu_min_arr.min()):.4f}, {float(mu_min_arr.max()):.4f}]")
print(f"tmu_min   range  : [{float(tmu_min_arr.min()):.4f}, {float(tmu_min_arr.max()):.4f}]")

Computing the global bounding box#

compute_bounding_box(moduli_sample) iterates over the sample, updating the global eigenvalue extrema, and returns the \(L^2\) radii of the bounding box for the NSNS-flux vector:

(2)#\[ \|h_1\|^2 \leq \frac{N_{\max}}{s_{\min}\,\tilde\mu_{\min}^{\rm gl}}, \quad \|h_2\|^2 \leq \frac{N_{\max}}{s_{\min}\,\mu_{\min}^{\rm gl}}, \quad \|h\|^2 \leq \frac{2\,\lambda_{\max}^{\rm gl}\,N_{\max}}{s_{\min}}.\]

The global extrema are also accessible as attributes of bf after the call.

moduli_sample = sampler.get_complex_moduli(500)
h1_box, h2_box, h_box = bf.compute_bounding_box(moduli_sample)

print("Global eigenvalue extrema:")
print(f"  lambda_max_gl   = {bf.lambda_max_gl:.4f}")
print(f"  mu_min_gl       = {bf.mu_min_gl:.4f}")
print(f"  tilde_mu_min_gl = {bf.tilde_mu_min_gl:.4f}")
print()
print("Bounding box radii (L2 norms):")
print(f"  h1_box = {h1_box:.3f}   (|h1|^2 <= {h1_box**2:.2f})")
print(f"  h2_box = {h2_box:.3f}   (|h2|^2 <= {h2_box**2:.2f})")
print(f"  h_box  = {h_box:.3f}   (|h|^2  <= {h_box**2:.2f})")
print()
print(f"Dilaton upper bound: dil_max = {bf.dil_max:.4f}")

Pre-computing eigenvalue bounds#

For large-scale scans, computing eigenvalue bounds separately (once, with many moduli points) is more efficient than recomputing them inside each enumerate_fluxes or sample_bounded_fluxes call. The compute_eigenvalue_bounds method does this and stores the results as class attributes:

# Pre-compute once and reuse the cached bounds below.
# Increase n_sample for production scans; the tutorial default is deliberately small.
bf.reset_eigenvalue_bounds()  # clear any previous bounds
h1_box, h2_box, h_box = bf.compute_eigenvalue_bounds(n_sample=2_000)
print(f"\nbounds_initialized = {bf.bounds_initialized}")

The box dimensions can always be retrieved later:

h1_box, h2_box, h_box = bf.get_h_box()
print(h1_box, h2_box, h_box)

Checking bounds for a given flux#

check_bounds(moduli, tau, flux) evaluates all eigenvalue inequalities of arXiv:2501.03984 for a specific flux configuration. It first calls update_local to compute the local eigenvalues and flux sub-vector norms, then evaluates every bound_* method.

Let us demonstrate this with a sample flux configuration:

# Draw a random initial guess from the sampler
seed = 42
rns_key = jvc.PRNGSequence(seed)

sampler.update_interior_points(num_pts=5000)
moduli_0, tau_0, flux_0 = sampler.initial_guesses(1, rns_key=rns_key)
moduli_0 = moduli_0[0]
tau_0    = complex(tau_0[0])
flux_0   = flux_0[0]

print("moduli =", moduli_0)
print("tau    =", tau_0)
print("flux   =", flux_0)

results = bf.check_bounds(moduli_0, tau_0, flux_0)

print(f"{'Bound':<20} {'Result'}")
print("-" * 40)
for result, label in results:
    print(f"{label:<20} {result}")

The local eigenvalue state is now stored in bf and can be inspected directly:

print(f"s          = {bf.s:.4f}")
print(f"c0         = {bf.c0:.4f}")
print(f"N_flux     = {bf.Nflux:.4f}")
print(f"lambda_max = {bf.lambda_max:.4f}")
print(f"mu_min     = {bf.mu_min:.4f}")
print(f"tmu_min    = {bf.tilde_mu_min:.4f}")
print()
print(f"||h||^2  = {bf.hnorm:.4f}")
print(f"||f||^2  = {bf.fnorm:.4f}")
print(f"||h1||^2 = {bf.h1norm:.4f}")
print(f"||h2||^2 = {bf.h2norm:.4f}")

check_bounds_flat() aggregates all bounds into a single pass/fail flag (no arguments — uses the current local state):

all_pass, detailed = bf.check_bounds_flat()
print(f"All bounds satisfied: {all_pass}")

Step-by-step walkthrough#

This section opens sample_bounded_fluxes at high resolution: each of the six internal stages — eigenvalue bounds, NS-NS box, ISD completion, tadpole + eigenvalue filters, ISD-condition sanity check, Newton/scipy refinement — gets its own code cell and a short prose introduction. The walkthrough uses a wider moduli/dilaton region than §4 so the box stays enumerable for cross-checks. The shadowed sampler and bf here are walkthrough-local; later sections build their own.

Model and sampler setup#

We use the same \(h^{1,2}=2\) Kreuzer–Skarke geometry as in the earlier flux-bounding examples. The moduli region is:

(3)#\[ \mathrm{Im}(z_i) \in [2, 3],\quad s = \mathrm{Im}(\tau) \in [2, 5],\quad c_0 = \mathrm{Re}(\tau) \in [-0.5, 0.5].\]

We set \(N_{\max} = 5\) so that the bounding box is small enough to enumerate exhaustively as a cross-check.

DILATON_BOUNDS  = [2, 5]
MODULI_BOUNDS   = [2, 3]
AXION_BOUNDS    = [-0.5, 0.5]
NMAX            = 5

sampler = jvc.data_sampler(
    model,
    flux_bounds=[-5, 5],
    moduli_bounds=MODULI_BOUNDS,
    axion_bounds=AXION_BOUNDS,
    dilaton_bounds=DILATON_BOUNDS,
    seed=0,
)

print(f"s   ∈ [{sampler.s_lower}, {sampler.s_upper}]")
print(f"Im(z) ∈ [{sampler.moduli_lower}, {sampler.moduli_upper}]")
print(f"c₀  ∈ [{sampler.axion_lower}, {sampler.axion_upper}]")

Step 1 — Global eigenvalue bounds#

sample_bounded_fluxes first samples \(n_{\mathrm{sample}}\) moduli points, computes three eigenvalue quantities at each point, and takes global extrema.

Quantity	Matrix	Role
\(\lambda_{\max}\)	ISD matrix \(\mathcal{M}\)	bounds \(\|h\|^2\) and \(\|f\|^2\)
\(\mu_{\min/\max}\)	\(-\mathrm{Im}(\mathcal{N})\)	bounds \(\|h_2\|^2\) (B-cycle NSNS sector)
\(\tilde{\mu}_{\min/\max}\)	\(\mathrm{Im}(\mathcal{N}^{-1})\)	bounds \(\|h_1\|^2\) (A-cycle NSNS sector)

Below we reproduce this step manually so we can inspect every number.

N_SAMPLE = 200_000
np.random.seed(42)

moduli_sample = sampler.get_complex_moduli(N_SAMPLE)
tau_sample    = sampler.get_complex_tau(N_SAMPLE)

flag = np.all(moduli_sample.imag>MODULI_BOUNDS[0],axis=1)
print(np.where(flag)[0].shape[0], "points with Im(z) >", MODULI_BOUNDS[0])
moduli_sample = moduli_sample[flag]
tau_sample    = tau_sample[flag]

print(f"Sampled {len(moduli_sample)} moduli points.")
print(f"Im(z) range: [{moduli_sample.imag.min():.3f}, {moduli_sample.imag.max():.3f}]")
print(f"s range:     [{tau_sample.imag.min():.3f}, {tau_sample.imag.max():.3f}]")

# Compute eigenvalues at a single representative moduli point.
z0  = moduli_sample[0]
tau0 = tau_sample[0]
z0c = jnp.conj(z0)

# Gauge kinetic matrix N(z, z_bar)
N_mat = model.gauge_kinetic_matrix(z0, z0c)
print("N matrix (gauge kinetic):")
print(N_mat)

# Eigenvalues of -Im(N)
mu_evs = jnp.linalg.eigvalsh(-N_mat.imag)
print(f"\nEigenvalues of -Im(N): {np.array(mu_evs)}")
print(f"  μ_min = {float(mu_evs.min()):.6f},  μ_max = {float(mu_evs.max()):.6f}")

# Eigenvalues of Im(N^{-1})
tmu_evs = jnp.linalg.eigvalsh(jnp.linalg.inv(N_mat).imag)
print(f"\nEigenvalues of Im(N⁻¹): {np.array(tmu_evs)}")
print(f"  μ̃_min = {float(tmu_evs.min()):.6f},  μ̃_max = {float(tmu_evs.max()):.6f}")

# ISD matrix M(z, z_bar) — formula from periods.py:
#   Block1 = [Im(N) + Re(N) @ Im(N)^{-1} @ Re(N)  ,  Re(N) @ Im(N)^{-1}]
#   Block2 = [Im(N)^{-1} @ Re(N)                   ,  Im(N)^{-1}         ]
#   M = -[Block1; Block2]

M0 = model.ISD_matrix(z0, z0c)
lm_evs = jnp.linalg.eigvalsh(M0)
print("ISD matrix M shape:", M0.shape)
print(f"ISD matrix eigenvalues: {np.array(lm_evs)}")
print(f"  λ_max = {float(lm_evs.max()):.6f}")
print()
# Verify M is positive definite
print(f"All eigenvalues positive? {bool(jnp.all(lm_evs > 0))}")

# Now compute global extrema over the full moduli sample using the JIT-compiled vmap.
# This is what compute_bounding_box does internally.
bf = bounded_fluxes(model, sampler=sampler, Nmax=NMAX)

# The dil_min should tighten from sqrt(3)/2 to sampler.s_lower=2
# (done inside sample_bounded_fluxes / enumerate_fluxes automatically)
bf.dil_min = float(sampler.s_lower)  # manual tightening for this walkthrough
print(f"dil_min (used for bounds) = {bf.dil_min}  (= sampler.s_lower)")

# JIT-vmapped eigenvalue computation over all moduli points
evs_all, M_all = bf._compute_evs_and_M_vmap(jnp.array(moduli_sample, dtype=complex))
lm_arr, mu_min_arr, mu_max_arr, tmu_min_arr, tmu_max_arr = evs_all

lambda_max_gl   = float(jnp.max(lm_arr))
mu_min_gl       = float(jnp.min(mu_min_arr))
mu_max_gl       = float(jnp.max(mu_max_arr))
tmu_min_gl      = float(jnp.min(tmu_min_arr))
tmu_max_gl      = float(jnp.max(tmu_max_arr))

print(f"\nGlobal eigenvalue extrema over {N_SAMPLE} moduli points:")
print(f"  λ_max   = {lambda_max_gl:.4f}")
print(f"  μ_min   = {mu_min_gl:.4f}")
print(f"  μ_max   = {mu_max_gl:.4f}")
print(f"  μ̃_min  = {tmu_min_gl:.6f}")
print(f"  μ̃_max  = {tmu_max_gl:.4f}")

fig, axes = plt.subplots(1, 3, dpi=130, figsize=(13, 3.2))

axes[0].hist(np.array(lm_arr), bins=40, edgecolor='k', lw=0.4)
axes[0].axvline(lambda_max_gl, color='r', lw=1.5, label=f'λ_max={lambda_max_gl:.1f}')
axes[0].set_xlabel(r'$\lambda_{\max}(\mathcal{M})$', fontsize=11)
axes[0].set_ylabel('Count', fontsize=11)
axes[0].set_title('ISD matrix max eigenvalue', fontsize=10)
axes[0].legend(fontsize=8)
axes[0].set_yscale('log')

axes[1].hist(np.array(mu_min_arr), bins=40, edgecolor='k', lw=0.4)
axes[1].axvline(mu_min_gl, color='r', lw=1.5, label=f'μ_min={mu_min_gl:.3f}')
axes[1].set_xlabel(r'$\mu_{\min}(-\mathrm{Im}\,\mathcal{N})$', fontsize=11)
axes[1].set_title(r'Gauge kinetic $\mu_{\min}$', fontsize=10)
axes[1].legend(fontsize=8)
axes[1].set_yscale('log')

axes[2].hist(np.array(tmu_min_arr), bins=40, edgecolor='k', lw=0.4)
axes[2].axvline(tmu_min_gl, color='r', lw=1.5, label=f'μ̃_min={tmu_min_gl:.5f}')
axes[2].set_xlabel(r'$\tilde\mu_{\min}(\mathrm{Im}\,\mathcal{N}^{-1})$', fontsize=11)
axes[2].set_title(r'Inverse gauge kinetic $\tilde\mu_{\min}$', fontsize=10)
axes[2].legend(fontsize=8)
axes[2].set_yscale('log')

plt.tight_layout()
plt.show()

Step 2 — Bounding box for NSNS-flux \(h\)#

From the inequalities in arXiv:2501.03984 the allowed integer NSNS-fluxes \(h = (h_1, h_2)\) must satisfy:

(4)#\[ \tilde\mu_{\min}\,\|h_1\|^2 \leq \frac{N_{\max}}{s_{\min}},\qquad \mu_{\min}\,\|h_2\|^2 \leq \frac{N_{\max}}{s_{\min}},\qquad \|h\|^2 \leq \frac{2\,\lambda_{\max}\,N_{\max}}{s_{\min}}.\]

Here \(s_{\min}\) is dil_min, which we set to sampler.s_lower = 2.

The key insight is that \(s_{\min}\) must reflect the actual lower bound of the sampled region, not the SL(2,Z) fundamental domain floor \(\sqrt{3}/2 \approx 0.866\). Using the correct \(s_{\min} = 2\) tightens the box radii by a factor \(\sqrt{2/0.866} \approx 1.5\), reducing box volume by \(\sim 3.4\times\) for each component.

s_min = bf.dil_min  # = 2.0 after tightening

h1_box = np.sqrt(NMAX / (s_min * tmu_min_gl))
h2_box = np.sqrt(NMAX / (s_min * mu_min_gl))
h_box  = np.sqrt(2.0 * lambda_max_gl * NMAX / s_min)

print(f"Bounding box (with s_min = {s_min}):")
print(f"  h1_box = sqrt(Nmax / (s_min * μ̃_min)) = sqrt({NMAX} / ({s_min} * {tmu_min_gl:.5f})) = {h1_box:.3f}")
print(f"  h2_box = sqrt(Nmax / (s_min * μ_min))  = sqrt({NMAX} / ({s_min} * {mu_min_gl:.5f})) = {h2_box:.3f}")
print(f"  h_box  = sqrt(2*λ_max*Nmax / s_min)    = sqrt(2*{lambda_max_gl:.3f}*{NMAX} / {s_min})  = {h_box:.3f}")

# Cross-check against compute_bounding_box
h1_box_cc, h2_box_cc, h_box_cc = bf.compute_bounding_box(
    jnp.array(moduli_sample, dtype=complex)
)
print(f"\nFrom compute_bounding_box:")
print(f"  h1_box={h1_box_cc:.3f}, h2_box={h2_box_cc:.3f}, h_box={h_box_cc:.3f}")
print(f"  (matches: {np.allclose([h1_box, h2_box, h_box], [h1_box_cc, h2_box_cc, h_box_cc], rtol=1e-4)})")

# What would happen if we forgot to tighten dil_min?
s_min_wrong = np.sqrt(3) / 2  # SL(2,Z) fundamental domain floor
h1_box_wrong = np.sqrt(NMAX / (s_min_wrong * tmu_min_gl))
h2_box_wrong = np.sqrt(NMAX / (s_min_wrong * mu_min_gl))
h_box_wrong  = np.sqrt(2.0 * lambda_max_gl * NMAX / s_min_wrong)

print(f"Comparison of bounding box radii:")
print(f"  s_min = {s_min_wrong:.4f} (wrong, SL(2,Z) floor):  h_box = {h_box_wrong:.2f}")
print(f"  s_min = {s_min:.4f} (correct, sampler floor): h_box = {h_box:.2f}")
print(f"  Volume ratio: ({h_box_wrong:.2f}/{h_box:.2f})^4 ≈ {(h_box_wrong/h_box)**4:.1f}x more candidates with wrong s_min")

# Enumerate all integer h vectors inside the bounding box (for small Nmax=2 this is feasible)

h_candidates = bf.get_h_candidates()
print(f"Total integer h candidates inside bounding box: {len(h_candidates)}")
print(f"h vector shape: {h_candidates.shape}  (n_candidates x n_fluxes)")
print(f"n_fluxes = {model.n_fluxes}, dimension_H3 = {model.dimension_H3}")
print(f"  → h = [h₁ | h₂] where h₁,h₂ each have {model.dimension_H3} components")

Step 3 — ISD completion: \(h \mapsto f\)#

For each NSNS-flux \(h\) and each sampled moduli point \((z, \tau)\), we compute the ISD-projected RR-flux:

(5)#\[ f = s\,\mathcal{M}(z,\bar z)\,\Sigma\,h + c_0\,h, \qquad \tau = c_0 + \mathrm{i}s.\]

This is mode "H" in sampling.py’s _ISD_sampling_FH. The result \(f\) is continuous (not integer); we round to the nearest integer and later check the tadpole constraint.

Convention verification#

We verify that flux_bounding.py’s formula matches sampling.py’s formula exactly.

# Take the first moduli point as our working example
z_ex  = jnp.array(moduli_sample[0], dtype=complex)
tau_ex = jnp.array(tau_sample[0], dtype=complex)
z_ex_c = jnp.conj(z_ex)
tau_ex_c = jnp.conj(tau_ex)

c0_ex = float(jnp.real(tau_ex))
s_ex  = float(jnp.imag(tau_ex))

print(f"Working moduli point:")
print(f"  z   = {np.array(z_ex)}")
print(f"  τ   = {complex(tau_ex):.6f}  (c₀={c0_ex:.4f}, s={s_ex:.4f})")

# Symplectic form sigma from model
sigma = model.periods.sigma
print(f"\nSymplectic form σ:")
print(np.array(sigma).astype(int))

# Pick a sample h vector
h_test = jnp.array([0, 1, 2, 0, 0, 1], dtype=float)  # [h1_1, h1_2, h2_1, h2_2]

# Method A: flux_bounding.py formula
M0_ex = model.ISD_matrix(z_ex, z_ex_c)
M0_sigma = M0_ex @ sigma
f_bounding = s_ex * (M0_sigma @ h_test) + c0_ex * h_test

# Method B: sampling.py _ISD_sampling_FH mode="H" formula
SigmaFlux = sigma @ h_test              # σ @ h
M0_SigmaFlux = M0_ex @ SigmaFlux        # M₀ @ σ @ h
f_sampling = M0_SigmaFlux * s_ex + h_test * c0_ex

# Method C: sampling.py's ISD_sampling function (the official API)
#flux_full_test = jnp.concatenate([jnp.zeros(model.n_fluxes), h_test])  # [f=0 | h]
flux_full_test = h_test
f_official = sampler.ISD_sampling(
    z_ex, z_ex_c, tau_ex, tau_ex_c,
    flux_full_test,
    mode="H", output="half", return_integer_flux=False,
)

print("ISD completion: f = s * M₀ @ (σ @ h) + c₀ * h")
print(f"  f (flux_bounding formula) = {np.array(f_bounding)}")
print(f"  f (sampling.py formula)   = {np.array(f_sampling)}")
print(f"  f (ISD_sampling API)      = {np.array(f_official)}")
print()
print(f"  All agree: {np.allclose(f_bounding, f_sampling) and np.allclose(f_bounding, f_official)}")

# Batch ISD completion: apply to all h candidates in parallel
# This is the inner loop of _process_h_at_modulus_jit

h_chunk = jnp.array(h_candidates[:200], dtype=float)  # first 200 for display

# Batched formula: (M0_sigma @ h_chunk.T).T  applies M0_sigma to each h as a column
f_chunk = s_ex * (M0_sigma @ h_chunk.T).T + c0_ex * h_chunk
f_int_chunk = jnp.round(f_chunk)  # round to nearest integer

print(f"ISD-completed f for first 5 h candidates (before rounding):")
for i in range(5):
    print(f"  h = {np.array(h_candidates[i])}, f_cont = {np.array(f_chunk[i]).round(3)}, f_int = {np.array(f_int_chunk[i]).astype(int)}")

Step 4 — Tadpole and eigenvalue bound filters#

After ISD completion we apply two filters to each candidate flux \([f | h]\):

Tadpole constraint: \(N_{\mathrm{flux}} = f^T \Sigma\, h \in (0, N_{\max}]\)
Eigenvalue bounds: a set of local and global inequalities from arXiv:2501.03984 involving \(\|h\|^2, \|f\|^2\), and the eigenvalues at this moduli point.

We verify the tadpole formula first.

Sign convention. We use the signed tadpole \(N_{\mathrm{flux}} = f^T \Sigma\, h\) with \(\Sigma = \begin{pmatrix}0 & \mathbb{1}\\ -\mathbb{1} & 0\end{pmatrix}= \) model.periods.sigma. This equals model.tadpole(flux) exactly. Since \(\Sigma\) is antisymmetric, \(f^T\Sigma h = -\,h^T\Sigma f\); for physical (ISD-completed) fluxes \(N_{\mathrm{flux}}>0\), so no absolute value is needed. The same convention is cross-checked against arXiv:2501.03984 in the Dataset B case study (§6).

# D3-tadpole formula: N_flux = f @ sigma @ h  (antisymmetric bilinear form)
# For a full flux vector [f | h] of length 2*n_fluxes:
#   N_flux = flux[:n_fl] @ sigma @ flux[n_fl:]

n_fl = model.n_fluxes
dim  = model.dimension_H3

# Use the first ISD-completed example
h_ex = jnp.array(h_candidates[10], dtype=float)
f_ex = jnp.round(s_ex * (M0_sigma @ h_ex) + c0_ex * h_ex)
flux_ex = jnp.concatenate([f_ex, h_ex])

# Compute tadpole three ways:
N1 = float((f_ex @ sigma @ h_ex).real)
N2 = float((flux_ex[:n_fl] @ jnp.dot(sigma, flux_ex[n_fl:])).real)
N3 = float(model.tadpole(flux_ex).real)

print(f"Flux example: f = {np.array(f_ex).astype(int)}, h = {np.array(h_ex).astype(int)}")
print(f"Tadpole N_flux = f @ σ @ h:")
print(f"  Method 1 (direct):          {N1:.6f}")
print(f"  Method 2 (via full vector): {N2:.6f}")
print(f"  Method 3 (model.tadpole):   {N3:.6f}")
print(f"  All agree: {np.isclose(N1, N2) and np.isclose(N1, N3)}")

# For an exactly ISD flux (not rounded), the tadpole has a clean form:
#   N_flux = s * h^T @ M^{-1} @ h  (positive for positive-definite M)

h_ex_f = jnp.array(h_candidates[10], dtype=float)
f_ex_cont = s_ex * (M0_sigma @ h_ex_f) + c0_ex * h_ex_f  # continuous (not rounded)

N_continuous = float((f_ex_cont @ sigma @ h_ex_f).real)
N_formula    = (s_ex * h_ex_f @ jnp.linalg.inv(M0_ex) @ h_ex_f).real

print("For the CONTINUOUS (pre-rounding) ISD flux:")
print(f"  N_flux = f_cont @ σ @ h        = {N_continuous:.8f}")
print(f"  N_flux = s * h^T M^{{-1}} h     = {N_formula:.8f}")
print(f"  Match (ISD identity):          {np.isclose(N_continuous, N_formula)}")
print()
print("Note: after rounding f to integers, N_flux changes slightly.")

# Run the JIT-compiled filter on h candidates at this one moduli point.
# We process the candidates in chunks: this avoids a large one-shot compile/allocation
# while still exercising the same kernel used by the production bounded-flux search.
from jaxvacua.flux_bounding import _process_h_at_modulus_jit

h_candidates_arr = jnp.array(h_candidates, dtype=float)
evs_ex = (
    float(jnp.max(jnp.linalg.eigvalsh(M0_ex))),
    float(jnp.min(jnp.linalg.eigvalsh(-N_mat.imag))),
    float(jnp.max(jnp.linalg.eigvalsh(-N_mat.imag))),
    float(jnp.min(jnp.linalg.eigvalsh(jnp.linalg.inv(N_mat).imag))),
    float(jnp.max(jnp.linalg.eigvalsh(jnp.linalg.inv(N_mat).imag))),
)
dil_max = float(lambda_max_gl * NMAX)

chunk_size = min(1024, len(h_candidates_arr))
flux_chunks = []
valid_chunks = []
n_chunks = 0
for start in tqdm(range(0, len(h_candidates_arr), chunk_size)):
    h_chunk = h_candidates_arr[start:start + chunk_size]
    flux_chunk, valid_chunk = _process_h_at_modulus_jit(
        h_chunk,
        M0_sigma,
        n_fl, dim,
        jnp.array(s_ex),
        jnp.array(c0_ex),
        evs_ex,
        tau_ex,
        sigma,
        lambda_max_gl, mu_min_gl, mu_max_gl, tmu_min_gl, tmu_max_gl,
        float(s_min), dil_max, float(NMAX),
    )
    flux_chunks.append(flux_chunk)
    valid_chunks.append(valid_chunk)
    n_chunks += len(flux_chunk)
    
print(f"Processed {n_chunks} flux candidates.")
flux_all = jnp.concatenate(flux_chunks, axis=0)
valid_all = jnp.concatenate(valid_chunks, axis=0)
n_valid = int(jnp.sum(valid_all))
print(f"Candidates passing all filters (at moduli point 0): {n_valid} / {len(h_candidates)}")
print(f"Pass rate: {100*n_valid/len(h_candidates):.2f}%")

# Manually check the tadpole for all passing candidates.
valid_idx = np.where(np.array(valid_all))[0]
flux_valid = np.array(flux_all)[valid_idx]

if len(flux_valid) == 0:
    print("No candidates passed all filters at this illustrative moduli point.")
else:
    tadpoles = np.array([float(model.tadpole(jnp.array(fv)).real) for fv in flux_valid])
    print("Tadpoles of passing candidates:")
    print(f"  min={tadpoles.min():.1f}, max={tadpoles.max():.1f}")
    print(f"  All in (0, {NMAX}]: {bool(np.all((tadpoles > 0) & (tadpoles <= NMAX)))}")
    print()
    print("First few passing fluxes [f | h] (integer):")
    for i, (fv, tad) in enumerate(zip(flux_valid[:5], tadpoles[:5])):
        print(f"  [{i}] {fv.astype(int).tolist()}  N_flux={tad:.0f}")

Step 5 — ISD convention consistency check#

We perform a direct cross-check: take a single ISD flux constructed by sampling.py’s initial_guesses_ISD, and verify that:

The ISD condition \(\star G_3 = \mathrm{i}\,G_3\) is approximately satisfied (up to the integrality rounding)
The tadpole \(N_{\mathrm{flux}} = f^T \Sigma h\) is positive and within bounds
The flux_bounding.py formula gives the same \(f\) for the same \((h, z, \tau)\)

# ISD condition: f = s * M(z) * sigma * h + c0 * h
# Equivalently: (f - c0*h) = s * M(z) * sigma * h
# Or: ftilde = f - c0*h satisfies ftilde = s * M * sigma * h
#
# Check the degree to which the integer-rounded ISD flux satisfies this.

f_int = jnp.round(f_ex_cont)  # integer-rounded ISD flux
ftilde = f_int - c0_ex * h_ex_f
rhs    = s_ex * (M0_sigma @ h_ex_f)

print("ISD residual check for integer-rounded flux:")
print(f"  f_tilde  = f - c₀ h = {np.array(ftilde).round(4)}")
print(f"  s M σ h             = {np.array(rhs).round(4)}")
print(f"  ||f_tilde - s M σ h|| = {float(jnp.linalg.norm(ftilde - rhs)):.4f}  (rounding error)")

# Full ISD condition in terms of G3 = f - tau*h (complex flux)
# Self-duality: G3 = (f - tau*h) should satisfy M * sigma * G3 = -G3
# (M is the ISD projector with eigenvalue -1 for ISD fluxes)

# For integer flux we check how well this is satisfied
G3 = f_int - tau_ex * h_ex_f
sigma_G3 = sigma @ G3
M_sigma_G3 = M0_ex @ sigma_G3

print("Self-duality check: M σ G₃ ≈ -G₃ ?")
print(f"  G₃         = {np.array(G3.real).round(3)} + i * {np.array(G3.imag).round(3)}")
print(f"  M σ G₃     = {np.array(M_sigma_G3.real).round(3)} + i * {np.array(M_sigma_G3.imag).round(3)}")
print(f"  -G₃        = {np.array(-G3.real).round(3)} + i * {np.array(-G3.imag).round(3)}")
print()
# For the continuous flux this should be exact
G3_cont = f_ex_cont - tau_ex * h_ex_f
sigma_G3_cont = sigma @ G3_cont
M_sigma_G3_cont = M0_ex @ sigma_G3_cont
print(f"  ||M σ G₃_cont - (-G₃_cont)|| = {float(jnp.linalg.norm(M_sigma_G3_cont + G3_cont)):.2e}  (continuous, should be ~0)")

Initial guess strategy#

The moduli point at which ISD completion was performed is an excellent initial guess: the rounded integer flux is close to ISD at that point, so \(D_I W\) is small there.

# Run sample_bounded_fluxes with refine=False to get a small pool of initial guesses.
# The production settings used in the paper are much larger; NB17 keeps those as opt-in case studies.
bf2 = bounded_fluxes(model, sampler=sampler, Nmax=NMAX)

candidates = bf2.sample_bounded_fluxes(
    n_target=100,
    n_batch=20_000,
    n_sample=500,
    n_mod=100,
    max_batches=5,
    verbose=True,
    seed=42,
    refine=False,
    return_moduli=True,
)
print(f"\nFound {len(candidates)} initial candidates.")

# Inspect the F-term residuals at the initial guess moduli points
print("F-term residuals |DW| at initial guess moduli:")
print(f"{'Flux':>12s}  {'|DW| pre-refinement':>52s}  N_flux")
print("-" * 60)
for r in candidates:
    z   = jnp.array(r["moduli"], dtype=complex)
    tau = jnp.array(r["tau"], dtype=complex)
    fl  = jnp.array(r["flux"])
    DW  = model.DW(z, jnp.conj(z), tau, jnp.conj(tau), fl)
    nfl = float(model.tadpole(fl).real)
    residual = float(jnp.linalg.norm(DW))
    if residual < 2:  # only show those with reasonably small residuals
        print(f"  {str(r['flux'].astype(int).tolist()).replace(" ",""):>12s}  {residual:>22.6e}  {nfl:.0f}")

# Define the root-finding problem for scipy.optimize.root
# Pack (Re(z1), Im(z1), Re(z2), Im(z2), Re(tau), Im(tau)) into a real 6-vector x

h12 = model.h12  # number of complex structure moduli

def pack(z, tau):
    """Complex moduli + tau -> real vector."""
    z_np = np.array(z)
    return np.concatenate([z_np.real, z_np.imag, [tau.real, tau.imag]])

def unpack(x):
    """Real vector -> (z_jax, tau_jax)."""
    z_re  = jnp.array(x[:h12])
    z_im  = jnp.array(x[h12:2*h12])
    z     = z_re + 1j * z_im
    tau   = x[2*h12] + 1j * x[2*h12 + 1]
    return z, tau

def residual_fn(x, flux):
    """F-term residual Re(DW) concatenated with Im(DW), as a real vector."""
    z, tau = unpack(x)
    DW = model.DW(z, jnp.conj(z), jnp.array(tau), jnp.conj(jnp.array(tau)), flux)
    return np.concatenate([np.array(DW.real), np.array(DW.imag)])

# Quick sanity check on a known initial guess
r0   = candidates[0]
fl0  = jnp.array(r0["flux"])
x0   = pack(r0["moduli"], complex(r0["tau"]))
res0 = residual_fn(x0, fl0)
print(f"Residual vector at initial guess: {res0.round(4)}")
print(f"||residual|| = {np.linalg.norm(res0):.4e}")

# Solve DW = 0 for a small subset of candidates using scipy.optimize.root.
# Keeping this subset small is important for laptop-safe execution.

TOL_SCIPY = 1e-10
MAX_SCIPY_CANDIDATES = 20

scipy_results = []
for r in tqdm(candidates[:MAX_SCIPY_CANDIDATES]):
    fl   = jnp.array(r["flux"])
    x0   = pack(r["moduli"], complex(r["tau"]))

    sol  = root(residual_fn, x0, args=(fl,), method="hybr", tol=TOL_SCIPY)
    z_sol, tau_sol = unpack(sol.x)

    flag = np.any(z_sol.imag > MODULI_BOUNDS[1]) or np.any(z_sol.imag < MODULI_BOUNDS[0])
    if flag:
        continue

    X = model._convert_complex_to_real_nondif(z_sol, tau_sol)
    model.is_physical(X)

    if not sol.success:
        continue
    print(f"Root found for {r['flux'].astype(int).tolist()}: {sol.message}")

    DW_final = residual_fn(sol.x, fl)
    res_final = np.linalg.norm(DW_final)
    converged = sol.success and res_final < TOL_SCIPY * 1e3

    z_sol, tau_sol, fl = model.map_to_fd(z_sol, tau_sol, fl)

    scipy_results.append({
        "flux":       np.array(fl),
        "moduli":     np.array(z_sol),
        "tau":        complex(tau_sol),
        "residual":   res_final,
        "converged":  converged,
        "scipy_msg":  sol.message,
    })

n_conv = sum(r["converged"] for r in scipy_results)
print(f"scipy.optimize.root: {n_conv}/{len(scipy_results)} converged from {min(len(candidates), MAX_SCIPY_CANDIDATES)} tested candidates (tol={TOL_SCIPY:.0e})")

# Filter converged solutions and check whether they lie inside the sampler's moduli patch.
# Patch conditions: Im(z_i) ∈ [moduli_lower, moduli_upper], s ∈ [s_lower, s_upper], c₀ ∈ [axion_lower, axion_upper]

def in_patch(z, tau):
    im_z = np.imag(z)
    s    = np.imag(tau)
    c0   = np.real(tau)
    return (
        np.all(im_z >= sampler.moduli_lower) and
        np.all(im_z <= sampler.moduli_upper) and
        s >= sampler.s_lower and s <= sampler.s_upper and
        c0 >= sampler.axion_lower and c0 <= sampler.axion_upper
    )

print(f"{'Flux':>45s}  {'|DW|':>12s}  {'in patch':>10s}  {'N_flux':>8s}")
print("-" * 90)
n_in_patch = 0
vacua = []
seen  = set()
for r in scipy_results:
    if not r["converged"]:
        continue
    key = tuple(int(x) for x in r["flux"])
    if key in seen:
        continue
    seen.add(key)
    patch = in_patch(r["moduli"], r["tau"])
    if patch:
        n_in_patch += 1
        vacua.append(r)
    fl_str = str(r["flux"].astype(int).tolist())
    nfl    = float(model.tadpole(jnp.array(r["flux"])).real)
    print(f"{fl_str:>45s}  {r['residual']:>12.2e}  {'yes' if patch else 'no':>10s}  {nfl:>8.0f}")

print(f"\nConverged and in-patch: {n_in_patch} unique vacua")

# Print the refined vacua in detail
if vacua:
    print(f"{'='*80}")
    print(f"Found {len(vacua)} refined SUSY vacua:")
    print(f"{'='*80}")
    for i, r in enumerate(vacua):
        z   = jnp.array(r["moduli"], dtype=complex)
        tau = jnp.array(r["tau"], dtype=complex)
        fl  = jnp.array(r["flux"])
        DW  = model.DW(z, jnp.conj(z), tau, jnp.conj(tau), fl)
        nfl = float(model.tadpole(fl).real)

        print(f"\nVacuum {i+1}:")
        print(f"  flux    = {r['flux'].astype(int).tolist()}")
        print(f"  moduli  = {np.array(r['moduli'])}")
        print(f"  tau     = {r['tau']:.8f}")
        print(f"  N_flux  = {nfl:.0f}")
        print(f"  |DW|    = {r['residual']:.2e}")
        print(f"  DW_vec  = {np.array(DW).round(12)}")
else:
    print("No converged in-patch vacua found in this run. Try increasing n_target or n_batch.")

Comparison: scipy vs Newton#

Let us verify that scipy.optimize.root and the built-in Newton method converge to the same vacua.

# Re-run sample_bounded_fluxes with Newton refinement on a small target.
bf3 = bounded_fluxes(model, sampler=sampler, Nmax=NMAX)

newton_vacua = bf3.sample_bounded_fluxes(
    n_target=20,
    n_batch=20_000,
    n_sample=500,
    n_mod=50,
    max_batches=5,
    verbose=True,
    seed=42,
    refine=True,
    newton_step_size=1,
    newton_tol=1e-10,
    newton_max_iters=100,
)
print(f"Newton refinement found {len(newton_vacua)} vacua.")
print(f"scipy.optimize.root found {len(vacua)} vacua.")

# Compare flux vectors found by both methods
scipy_flux_set  = {tuple(int(x) for x in r["flux"]) for r in vacua}
newton_flux_set = {tuple(int(x) for x in r["flux"]) for r in newton_vacua}

print(f"Flux vectors in scipy result:  {len(scipy_flux_set)}")
print(f"Flux vectors in Newton result: {len(newton_flux_set)}")
print(f"Intersection: {len(scipy_flux_set & newton_flux_set)}")
print(f"Only in scipy:  {scipy_flux_set - newton_flux_set}")
print(f"Only in Newton: {newton_flux_set - scipy_flux_set}")

# For fluxes found by both, compare the converged moduli
if scipy_flux_set & newton_flux_set:
    common_key = next(iter(scipy_flux_set & newton_flux_set))

    r_scipy  = next(r for r in vacua       if tuple(int(x) for x in r["flux"]) == common_key)
    r_newton = next(r for r in newton_vacua if tuple(int(x) for x in r["flux"]) == common_key)

    print(f"Flux: {list(common_key)}")
    print(f"  scipy  moduli = {r_scipy['moduli']}")
    print(f"  Newton moduli = {r_newton['moduli']}")
    print(f"  scipy  tau    = {r_scipy['tau']:.10f}")
    print(f"  Newton tau    = {r_newton['tau']:.10f}")
    print(f"  ||Δmoduli||   = {np.linalg.norm(r_scipy['moduli'] - r_newton['moduli']):.2e}")
    print(f"  |Δtau|        = {abs(r_scipy['tau'] - r_newton['tau']):.2e}")

Summary and consistency checks#

We now perform a final self-consistency check on each refined vacuum.

def check_vacuum(r, model, sampler, label=""):
    """Print a consistency table for one refined vacuum."""
    z   = jnp.array(r["moduli"], dtype=complex)
    tau = jnp.array(r["tau"], dtype=complex)
    fl  = jnp.array(r["flux"])

    DW   = model.DW(z, jnp.conj(z), tau, jnp.conj(tau), fl)
    nfl  = float(model.tadpole(fl).real)
    M0_v = model.ISD_matrix(z, jnp.conj(z))
    G3   = fl[:model.n_fluxes] - tau * fl[model.n_fluxes:]
    isd_res = float(jnp.linalg.norm(M0_v @ (model.periods.sigma @ G3) + G3))

    patch = in_patch(np.array(z), complex(tau))

    header = f"--- Vacuum {label} ---"
    print(header)
    print(f"  flux            = {np.array(fl).astype(int).tolist()}")
    print(f"  moduli Im(z)    = {np.imag(np.array(z))}")
    print(f"  moduli Re(z)    = {np.real(np.array(z))}")
    print(f"  tau             = {complex(tau):.8f}")
    print(f"  N_flux          = {nfl:.2f}  (should be integer in (0,{NMAX}])")
    print(f"  ||DW||          = {float(jnp.linalg.norm(DW)):.2e}  (should be < 1e-8)")
    print(f"  ||M σ G3 + G3|| = {isd_res:.2e}  (ISD residual; 0 = exact ISD)")
    print(f"  in sampler patch= {patch}")
    print()

if vacua:
    for i, r in enumerate(vacua[:3]):
        check_vacuum(r, model, sampler, label=str(i+1))
else:
    print("No vacua to check. Increase n_target or n_batch.")

# Tadpole and moduli distribution of all refined vacua
all_refined = vacua + newton_vacua  # combine both refinement methods
# Deduplicate by flux
seen_all = set()
all_refined_dedup = []
for r in all_refined:
    key = tuple(int(x) for x in r["flux"])
    if key not in seen_all:
        seen_all.add(key)
        all_refined_dedup.append(r)

print(f"Total unique refined vacua: {len(all_refined_dedup)}")

if len(all_refined_dedup) >= 2:
    tads = np.array([float(model.tadpole(jnp.array(r["flux"])).real) for r in all_refined_dedup])
    mods = np.array([r["moduli"] for r in all_refined_dedup])
    taus = np.array([r["tau"] for r in all_refined_dedup])

    fig, axes = plt.subplots(1, 3, dpi=130, figsize=(13, 3.5))

    axes[0].hist(tads, bins=range(0, NMAX + 2), edgecolor='k', lw=0.5)
    axes[0].set_xlabel(r'$N_{\rm flux}$', fontsize=12)
    axes[0].set_ylabel('Count', fontsize=12)
    axes[0].set_title('Tadpole distribution', fontsize=11)

    sc = axes[1].scatter(
        mods[:, 0].imag, mods[:, 1].imag,
        c=tads, cmap='viridis', s=20, alpha=0.8
    )
    plt.colorbar(sc, ax=axes[1], label=r'$N_{\rm flux}$')
    axes[1].set_xlabel(r'$\mathrm{Im}(z_1)$', fontsize=12)
    axes[1].set_ylabel(r'$\mathrm{Im}(z_2)$', fontsize=12)
    axes[1].set_title('Moduli distribution', fontsize=11)

    axes[2].scatter(taus.real, taus.imag, c=tads, cmap='viridis', s=20, alpha=0.8)
    axes[2].set_xlabel(r'$\mathrm{Re}(\tau)$', fontsize=12)
    axes[2].set_ylabel(r'$\mathrm{Im}(\tau) = s$', fontsize=12)
    axes[2].set_title(r'Axio-dilaton $\tau$', fontsize=11)

    plt.tight_layout()
    plt.show()
else:
    print("Not enough vacua to plot — try larger n_target.")

Enumeration vs stochastic sampling#

This section compares two complementary approaches to populating flux ensembles inside the eigenvalue-bounded box:

Systematic enumeration (bounded_fluxes.enumerate_fluxes): iterates every integer \(h\) in the bounding box, ISD-completes each, and applies tadpole + eigenvalue filters. Provably complete but scales as \(\sim N_{\max}^3\), so practical only for \(N_{\max} \lesssim 30\).
Stochastic Gaussian-prior sampling (bounded_fluxes.sample_bounded_fluxes): draws \(h \sim \mathcal{N}(0, \sigma^2 M)\) from the prior tuned to the ISD tadpole, then runs the same ISD/filter pipeline. Sublinear in \(N_{\max}\), gives up completeness in exchange. The Gaussian prior outperforms a uniform-box prior by ~50× yield (see NB07 §5).

Flux utility helpers#

The class provides several helpers to decompose flux vectors into their sub-components.

# Split flux = [f | h] into f and h (each of length n_fluxes = 2*(h12+1))
f, h = bf.get_fh(flux_0)
print("f =", f)
print("h =", h)

# Further split f = [f1 | f2] and h = [h1 | h2] (each of length dimension_H3 = h12+1)
h1, h2 = bf.get_flux_split(h)
f1, f2 = bf.get_flux_split(f)
print("h1 =", h1, "  h2 =", h2)
print("f1 =", f1, "  f2 =", f2)

# Or all at once:
h1, h2, f1, f2 = bf.get_subvector(flux_0)

# D3-tadpole charge
print(f"N_flux = {bf.get_nflux(flux_0):.4f}")

Full enumeration algorithm#

enumerate_fluxes() is the main entry point that combines all steps above. It:

Samples n_sample moduli points from the sampler and computes the global bounding box. If compute_eigenvalue_bounds() was called beforehand, this step is skipped.
Enumerates all integer \(h\) candidates inside the box.
For each \(h\), computes an ISD-projected \(f\) at up to n_isd_per_h moduli points.
Retains flux configurations that satisfy the D3-tadpole constraint and all local eigenvalue bounds.

New features (since v2):

compute_eigenvalue_bounds(n_sample) — pre-compute and cache eigenvalue bounds.
moduli_regions — stratified ISD starting points via Cartesian product across moduli dimensions.
use_linearised_shifts=True — iterative ISD refinement with flag-based early stopping (requires FluxVacuaFinder).
constraints — user-supplied constraint function (moduli, tau, flux) → bool.
chunk_size — override the default streaming chunk size.
KeyboardInterrupt — gracefully returns partial results found so far.

This requires a sampler to have been passed at initialisation.

# Re-initialise with a very tight tadpole to keep complete enumeration fast.
bf_enum = bounded_fluxes(model, sampler=sampler, Nmax=2)

# Pre-compute eigenvalue bounds (reused inside enumerate_fluxes).
bf_enum.compute_eigenvalue_bounds(n_sample=1_000)

valid_fluxes = bf_enum.enumerate_fluxes(
    n_sample=50,
    n_isd_per_h=5,
    verbose=True,
    confirm_streaming=False,
)

# Package as dicts (flux only) for the comparison / tadpole cells below.
enum_results = [{"flux": np.asarray(fv)} for fv in valid_fluxes]

if len(valid_fluxes) > 0:
    flux_mat = np.stack(valid_fluxes)   # shape (N_valid, 2*n_fluxes)

    # D3-tadpole for each valid flux
    tadpoles = np.array([
        bf_enum.get_nflux(jnp.array(fv)) for fv in valid_fluxes
    ])

    print(f"Valid flux configurations found: {len(valid_fluxes)}")
    print(f"N_flux range : [{tadpoles.min():.0f}, {tadpoles.max():.0f}]")
    print(f"N_flux mean  : {tadpoles.mean():.2f}")
    print()
    print("First three valid flux vectors [f | h]:")
    for fv in valid_fluxes[:3]:
        print(" ", fv.astype(int))
else:
    print("No valid flux vectors found — try increasing Nmax or the moduli sample size.")

import matplotlib.pyplot as plt

if len(valid_fluxes) > 0:
    fig, axes = plt.subplots(1, 2, dpi=150, figsize=(8, 3))

    axes[0].hist(tadpoles, bins=range(0, int(tadpoles.max()) + 2), edgecolor='k', linewidth=0.5)
    axes[0].set_xlabel(r"$N_{\rm flux}$", fontsize=12)
    axes[0].set_ylabel("Count", fontsize=12)
    axes[0].set_title("D3-tadpole distribution", fontsize=11)

    # h1 vs h2 sub-vector norms
    h1norms = np.array([
        bf_enum.compute_norm(bf_enum.get_flux_split(bf_enum.get_fh(jnp.array(fv))[1])[0])
        for fv in valid_fluxes
    ])
    h2norms = np.array([
        bf_enum.compute_norm(bf_enum.get_flux_split(bf_enum.get_fh(jnp.array(fv))[1])[1])
        for fv in valid_fluxes
    ])

    axes[1].scatter(h1norms, h2norms, s=15, alpha=0.7)
    axes[1].set_xlabel(r"$\|h_1\|^2$", fontsize=12)
    axes[1].set_ylabel(r"$\|h_2\|^2$", fontsize=12)
    axes[1].set_title(r"NSNS sub-vector norms", fontsize=11)
    h1b, h2b, _ = bf_enum.get_h_box()
    axes[1].axvline(h1b**2, color='r', linestyle='--', linewidth=0.8, label=r"$h_1$ box")
    axes[1].axhline(h2b**2, color='b', linestyle='--', linewidth=0.8, label=r"$h_2$ box")
    axes[1].legend(fontsize=9)

    plt.tight_layout()
    plt.show()

Stochastic alternative#

sample_bounded_fluxes(n_target, n_batch, ...) Monte-Carlos the Gaussian prior described above, ISD-completes, and applies the filters until n_target vacua are found or max_batches is exhausted. Below we run it at \(N_{\max}=20\) — where systematic enumeration starts to slow down — and compare yields.

bf_stoch = bounded_fluxes(model, sampler=sampler, Nmax=20)
bf_stoch.compute_eigenvalue_bounds(n_sample=2_000)

stoch_results = bf_stoch.sample_bounded_fluxes(
    n_target=100,
    n_batch=10_000,
    n_sample=500,
    n_mod=10,
    max_batches=10,
    refine=True,
    newton_step_size=1.,
    verbose=True,
    seed=42,
    return_moduli=True,
    moduli_regions=[(2., 2.5), (2.5, 3.)],
)
print(f"\nStochastic search found {len(stoch_results)} unique valid flux candidates.")

# Each result is a dict with "flux", "moduli", "tau"
if stoch_results:
    r = stoch_results[0]
    print("flux   =", r["flux"].astype(int))
    print("moduli =", r["moduli"])
    print("tau    =", r["tau"])
    z = r["moduli"]
    t = r["tau"]
    flux = r["flux"].astype(int)
    DW = model.DW(z,jnp.conj(z),t,jnp.conj(t),flux)
    print("DW    =", DW)

Comparison#

We check how many of the stochastic results also appear in the exhaustive enumeration.

# Build set of flux tuples from enumeration
enum_set = set()
for r in enum_results:
    enum_set.add(tuple(int(x) for x in r["flux"]))

# Check overlap
n_overlap = 0
for r in stoch_results:
    key = tuple(int(x) for x in r["flux"])
    if key in enum_set:
        n_overlap += 1

print(f"Enumeration:  {len(enum_results)} flux candidates")
print(f"Stochastic:   {len(stoch_results)} flux candidates")
print(f"Overlap:      {n_overlap} / {len(stoch_results)} stochastic results also found by enumeration")
if len(enum_results) > 0:
    print(f"Coverage:     {n_overlap}/{len(enum_results)} = {100*n_overlap/len(enum_results):.1f}% of enumerated vacua recovered")

Tadpole distributions#

def compute_tadpoles(bf_obj, results):
    return np.array([bf_obj.get_nflux(jnp.array(r["flux"])) for r in results])

if len(enum_results) > 0 and len(stoch_results) > 0:
    tad_enum  = compute_tadpoles(bf, enum_results)
    tad_stoch = compute_tadpoles(bf_stoch, stoch_results)

    fig, axes = plt.subplots(1, 3, dpi=150, figsize=(13, 3.5))

    bins = range(0, int(max(tad_enum.max(), tad_stoch.max())) + 2)

    axes[0].hist(tad_enum, bins=bins, edgecolor='k', linewidth=0.5, alpha=0.7, label='Enumeration')
    axes[0].hist(tad_stoch, bins=bins, edgecolor='k', linewidth=0.5, alpha=0.5, label='Stochastic')
    axes[0].set_xlabel(r"$N_{\rm flux}$", fontsize=12)
    axes[0].set_ylabel("Count", fontsize=12)
    axes[0].set_title("D3-tadpole distribution", fontsize=11)
    axes[0].legend(fontsize=9)

    # h-norm distributions
    def h_norms(bf_obj, results):
        return np.array([bf_obj.compute_norm(bf_obj.get_fh(jnp.array(r["flux"]))[1]) for r in results])

    hn_enum  = h_norms(bf, enum_results)
    hn_stoch = h_norms(bf_stoch, stoch_results)

    axes[1].hist(hn_enum, bins=20, edgecolor='k', linewidth=0.5, alpha=0.7, label='Enumeration')
    axes[1].hist(hn_stoch, bins=20, edgecolor='k', linewidth=0.5, alpha=0.5, label='Stochastic')
    axes[1].set_xlabel(r"$\|h\|^2$", fontsize=12)
    axes[1].set_ylabel("Count", fontsize=12)
    axes[1].set_title(r"NSNS-flux norm distribution", fontsize=11)
    axes[1].legend(fontsize=9)

    # Moduli scatter: Im(z_1) vs Im(z_2) at which each stochastic flux was found valid
    mod_arr = np.array([r["moduli"] for r in stoch_results])
    sc = axes[2].scatter(mod_arr[:, 0].imag, mod_arr[:, 1].imag,
                         c=tad_stoch, cmap="viridis", s=15, alpha=0.7)
    plt.colorbar(sc, ax=axes[2], label=r"$N_{\rm flux}$")
    axes[2].set_xlabel(r"$\mathrm{Im}(z_1)$", fontsize=12)
    axes[2].set_ylabel(r"$\mathrm{Im}(z_2)$", fontsize=12)
    axes[2].set_title("Moduli of stochastic vacua", fontsize=11)

    plt.tight_layout()
    plt.show()
else:
    print("Not enough data to plot.")

Scaling to large \(N_{\max}\)#

At \(N_{\max} = 200\), the bounding box contains far too many \(h\)-vectors for exhaustive enumeration. The stochastic sample_bounded_fluxes path is the relevant API at this scale, but even stochastic searches should be sized deliberately on a laptop.

This cell estimates how the box grows at large \(N_{\max}\). The actual stochastic search is useful, but it is disabled by default via RUN_LARGE_NMAX_DEMO = False because it can draw many batches and allocate large arrays.

bf_large = bounded_fluxes(model, sampler=sampler, Nmax=200)
bf_large.compute_eigenvalue_bounds(n_sample=2_000)

dim = bf_large.dimension_H3
h1_box, h2_box, h_box = bf_large._h1_box, bf_large._h2_box, bf_large._h_box
h1_max = int(np.ceil(h1_box))
h2_max = int(np.ceil(h2_box))
n_box = (2 * h1_max + 1) ** dim * (2 * h2_max + 1) ** dim

print(f"Nmax = {bf_large.Nmax}")
print(f"Bounding box: h1_box={h1_box:.2f}, h2_box={h2_box:.2f}, h_box={h_box:.2f}")
print(f"Unfiltered box size: {n_box:,} h-candidates")
print(f"  → Full enumeration is {'feasible' if n_box < 1_000_000 else 'INFEASIBLE'}")
print("The stochastic Nmax=200 scan below is disabled unless RUN_LARGE_NMAX_DEMO = True.")

if RUN_LARGE_NMAX_DEMO:
    large_results = bf_large.sample_bounded_fluxes(
        n_target=200,
        n_batch=50_000,
        n_sample=300,
        n_mod=10,
        max_batches=10,
        verbose=True,
        seed=123,
        return_moduli=True,
    )
    print(f"\nFound {len(large_results)} unique flux candidates at Nmax={bf_large.Nmax}.")
else:
    large_results = []
    print("RUN_LARGE_NMAX_DEMO = False — skipping the larger stochastic scan.")

if len(large_results) > 0:
    tad_large = compute_tadpoles(bf_large, large_results)
    mod_large = np.array([r["moduli"] for r in large_results])

    fig, axes = plt.subplots(1, 2, dpi=150, figsize=(9, 3.5))

    axes[0].hist(tad_large, bins=30, edgecolor='k', linewidth=0.5)
    axes[0].set_xlabel(r"$N_{\rm flux}$", fontsize=12)
    axes[0].set_ylabel("Count", fontsize=12)
    axes[0].set_title(f"Tadpole distribution ($N_{{\\max}} = {bf_large.Nmax}$)", fontsize=11)

    sc = axes[1].scatter(mod_large[:, 0].imag, mod_large[:, 1].imag,
                         c=tad_large, cmap="viridis", s=15, alpha=0.7)
    plt.colorbar(sc, ax=axes[1], label=r"$N_{\rm flux}$")
    axes[1].set_xlabel(r"$\mathrm{Im}(z_1)$", fontsize=12)
    axes[1].set_ylabel(r"$\mathrm{Im}(z_2)$", fontsize=12)
    axes[1].set_title("Moduli locations", fontsize=11)

    plt.tight_layout()
    plt.show()

    print(f"N_flux range: [{tad_large.min():.0f}, {tad_large.max():.0f}]")
    print(f"N_flux mean:  {tad_large.mean():.1f}")
else:
    print("No flux candidates found.")

Take-aways#

Eigenvalue bounds give a finite search region. Extremal eigenvalues of the ISD matrix and gauge-kinetic matrix bound the NS-NS flux components once a moduli region, dilaton floor and tadpole bound are fixed.
bounded_fluxes packages the workflow: eigenvalue bounds, NS-NS box construction, ISD completion, tadpole/eigenvalue filters and Newton/scipy refinement.
Use enumerate_fluxes only when the box is small enough. It is complete, but the Cartesian product grows rapidly with \(N_{\max}\) and dimension.
Use sample_bounded_fluxes for larger searches. It trades completeness for feasibility by sampling candidates from a physically motivated bounded distribution.
Precompute bounds when running several scans. Reusing compute_eigenvalue_bounds(...) avoids repeated expensive setup.
Full reproductions belong in NB17. The same algorithm validates the mirror-octic and Dataset B literature datasets, but those checks are intentionally separated from the laptop-safe tutorial.

Flux bounding: enumeration and stochastic sampling

Contents

Flux bounding: enumeration and stochastic sampling#

Recommended paths through this notebook#

Outline#

Setup#

Quick start: the 1D mirror octic#

Literature cross-checks live in NB17#

Algorithm overview#

Eigenvalue bounds#

Initialising `bounded_fluxes`#

Eigenvalue quantities at a single moduli point#

Computing the global bounding box#

Pre-computing eigenvalue bounds#

Checking bounds for a given flux#

Step-by-step walkthrough#

Model and sampler setup#

Step 1 — Global eigenvalue bounds#

Step 2 — Bounding box for NSNS-flux \(h\)#

Step 3 — ISD completion: \(h \mapsto f\)#

Convention verification#

Step 4 — Tadpole and eigenvalue bound filters#

Step 5 — ISD convention consistency check#

Step 6 — Refinement with `scipy.optimize.root`#

Initial guess strategy#

Comparison: scipy vs Newton#

Summary and consistency checks#

Enumeration vs stochastic sampling#

Flux utility helpers#

Full enumeration algorithm#

Stochastic alternative#

Comparison#

Tadpole distributions#

Scaling to large \(N_{\max}\)#

Newton refinement#

Take-aways#

Further reading#

Flux bounding: enumeration and stochastic sampling

Contents

Flux bounding: enumeration and stochastic sampling#

Recommended paths through this notebook#

Outline#

Setup#

Quick start: the 1D mirror octic#

Literature cross-checks live in NB17#

Algorithm overview#

Eigenvalue bounds#

Initialising bounded_fluxes#

Eigenvalue quantities at a single moduli point#

Computing the global bounding box#

Pre-computing eigenvalue bounds#

Checking bounds for a given flux#

Step-by-step walkthrough#

Model and sampler setup#

Step 1 — Global eigenvalue bounds#

Step 2 — Bounding box for NSNS-flux \(h\)#

Step 3 — ISD completion: \(h \mapsto f\)#

Convention verification#

Step 4 — Tadpole and eigenvalue bound filters#

Step 5 — ISD convention consistency check#

Step 6 — Refinement with scipy.optimize.root#

Initial guess strategy#

Comparison: scipy vs Newton#

Summary and consistency checks#

Enumeration vs stochastic sampling#

Flux utility helpers#

Full enumeration algorithm#

Stochastic alternative#

Comparison#

Tadpole distributions#

Scaling to large \(N_{\max}\)#

Newton refinement#

Take-aways#

Further reading#

Initialising `bounded_fluxes`#

Step 6 — Refinement with `scipy.optimize.root`#