.. DO NOT EDIT. .. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY. .. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE: .. "gallery/abp_to_spde_demo.py" .. LINE NUMBERS ARE GIVEN BELOW. .. only:: html .. note:: :class: sphx-glr-download-link-note :ref:`Go to the end ` to download the full example code. .. rst-class:: sphx-glr-example-title .. _sphx_glr_gallery_abp_to_spde_demo.py: Discovering Toner–Tu hydrodynamics from agent-based flocking ============================================================== .. note:: Uses the **experimental** SPDE toolbox — see :doc:`/spde/index`. A polar active-matter "round trip": simulate 10 000 self-propelled particles with Vicsek-like alignment, **coarse-grain** their positions and orientations into continuum fields ``(ρ, m_x, m_y)``, and let SFI + PASTIS **discover** the governing SPDE — without ever telling it which terms to expect. Pipeline: 1. **Microscopic ABPs** (`SFI.langevin.OverdampedProcess` + chunked neighbor-list rebuilds): self-propulsion + soft repulsion + short-range polar alignment, tuned to the *banded* flocking regime. 2. **Coarse-graining** (Gaussian-smoothed bilinear deposition, :func:`examples._gallery_utils.coarse_grain.coarse_grain_polar`): particles → density ``ρ`` and polar momentum ``m = ρ⟨\hat e_θ⟩`` on a periodic 64×64 grid. 3. **Overcomplete SPDE basis** (``GridLayout`` with one ``ScalarSector`` for ρ and one spatial ``VectorSector`` for m): ~30 candidate terms covering continuity, polar Landau ordering, Frank elasticity, advection, and pressure-like density gradients. 4. **Linear inference + PASTIS** sparsification recovers the dominant Toner–Tu terms. 5. **Bootstrap** the inferred SPDE from random initial conditions and compare its bands to the agent-based ones. .. rubric:: Tags synthetic · overdamped · multi-particle · linear · spde · pastis · interactions .. GENERATED FROM PYTHON SOURCE LINES 38-40 .. code-block:: Python :dedent: 1 .. GENERATED FROM PYTHON SOURCE LINES 75-94 Microscopic flocking model ---------------------------- Each agent carries position :math:`(x_i, y_i)` and heading :math:`\theta_i`. The force on the 3-vector state :math:`(x, y, \theta)` is .. math:: \mathbf F_i = \underbrace{c_0\, \hat e_{\theta_i}}_{\text{propulsion}} \;+\; \underbrace{\sum_j (-\varepsilon)\,e^{-r_{ij}/R_0}\, \hat r_{ij}}_{\text{soft repulsion}} \;+\; \underbrace{\hat z\sum_j A\,e^{-r_{ij}/L_a}\, \sin(\theta_j - \theta_i)}_{\text{Vicsek alignment}} acting on :math:`(x, y, \theta)` respectively, plus rotational noise. We work in the **banded** flocking regime: high alignment, low rotational noise — propagating density bands aligned along the mean heading direction emerge spontaneously. .. GENERATED FROM PYTHON SOURCE LINES 94-140 .. code-block:: Python from SFI.bases.pairs import ( angle_coupling, heading_vector, pair_direction, parametric_radial_kernel, ) from SFI.langevin import OverdampedProcess from SFI.langevin.chunked import simulate_chunked from SFI.statefunc import Basis from SFI.utils.neighbors import make_neighbor_extras # ── physical parameters ── N_particles = 10_000 Lx = Ly = 64.0 Nsteps = 1500 dt_sim = 0.05 seed = 7 # Particle force parameters (tuned for banded regime) c0_true = 0.50 # self-propulsion speed eps_true = 2.00 # repulsion strength R0_true = 0.50 # repulsion length A_true = 1.20 # alignment torque amplitude La_true = 1.50 # alignment kernel length # Anisotropic noise: low translational, moderate rotational D_xy = 0.005 D_theta = 0.10 D_matrix = jnp.diag(jnp.array([D_xy, D_xy, D_theta])) # Neighbor-list parameters # cutoff = 3.0 captures repulsion (R0=0.5, exp(-6) ~ 2e-3) and alignment # (La=1.5, exp(-2) ~ 0.13) tails to good accuracy. cutoff = 3.0 skin = 1.5 rebuild_every = 5 box = jnp.array([Lx, Ly]) box_np = np.array([Lx, Ly]) theta_F_exact = dict(c0=c0_true, eps=eps_true, R0=R0_true, A=A_true, La=La_true) print(f"N = {N_particles}, box = {Lx:.0f}×{Ly:.0f}, ρ̄ = {N_particles/(Lx*Ly):.2f}") .. rst-class:: sphx-glr-script-out .. code-block:: none N = 10000, box = 64×64, ρ̄ = 2.44 .. GENERATED FROM PYTHON SOURCE LINES 141-143 Building the parametric simulation force ------------------------------------------ .. GENERATED FROM PYTHON SOURCE LINES 143-172 .. code-block:: Python dim = 3 # (x, y, θ) per particle B_heading = heading_vector(dim=dim, angle_index=2) e_ij = pair_direction( dim=dim, box="extras", spatial_dims=slice(0, 2), embed_dim=dim, embed_axes=[0, 1], ) g_align = angle_coupling(jnp.sin, dim=dim, angle_index=2) k_repel = parametric_radial_kernel( lambda r, p: -p["eps"] * jnp.exp(-r / p["R0"]), params={"eps": (), "R0": ()}, dim=dim, box="extras", spatial_dims=slice(0, 2), ) k_align = parametric_radial_kernel( lambda r, p: p["A"] * jnp.exp(-r / p["La"]), params={"A": (), "La": ()}, dim=dim, box="extras", spatial_dims=slice(0, 2), ) csr_kw = dict(indptr_key="indptr", indices_key="indices") F_sim = ( B_heading.to_psf(coeff_key="c0") + (k_repel * e_ij).dispatch_pairs_from_extras(**csr_kw, return_as="psf") + (k_align * g_align).dispatch_pairs_from_extras(**csr_kw, return_as="psf") ) .. GENERATED FROM PYTHON SOURCE LINES 173-179 Chunked simulation in the banded flocking regime -------------------------------------------------- Initial conditions: uniform positions, uniform random headings. We let the system burn in (the first ~20 time units relax to the polar-ordered manifold; bands sharpen later). .. GENERATED FROM PYTHON SOURCE LINES 179-286 .. code-block:: Python key = random.PRNGKey(seed) key, kx, kth = random.split(key, 3) X0_xy = random.uniform(kx, (N_particles, 2)) * box TH0 = random.uniform(kth, (N_particles,), minval=-jnp.pi, maxval=jnp.pi) x0 = jnp.concatenate([X0_xy, TH0[:, None]], axis=1) print("Building initial neighbor list ...") t0 = time.perf_counter() nbr0 = make_neighbor_extras(np.asarray(x0[:, :2]), cutoff + skin, box_np) print(f" nnz = {len(nbr0['indices'])}, " f"⟨neighbors⟩ = {len(nbr0['indices']) / N_particles:.1f} " f"({time.perf_counter() - t0:.2f}s)") extras0 = {"box": box} extras0.update(nbr0) proc = OverdampedProcess(F_sim, D=D_matrix, extras_global=extras0) proc.set_params(theta_F=theta_F_exact) proc.initialize(x0) # Optional disk cache for the microscopic trajectory so re-running the # demo (e.g. to iterate on plots or inference) does not redo the # expensive ABP integration. ``__file__`` is undefined when running # under sphinx-gallery's ``exec``; we search upward for the canonical # ``examples/gallery/_cache`` directory. def _find_cache_dir() -> str: try: start = os.path.dirname(os.path.abspath(__file__)) except NameError: start = os.getcwd() cur = start for _ in range(8): cand = os.path.join(cur, "examples", "gallery", "_cache") if os.path.isdir(cand): return cand if os.path.basename(cur) == "gallery" and \ os.path.basename(os.path.dirname(cur)) == "examples": cand = os.path.join(cur, "_cache") if os.path.isdir(cand): return cand nxt = os.path.dirname(cur) if nxt == cur: break cur = nxt return os.path.join(start, "_cache") _CACHE_DIR = _find_cache_dir() os.makedirs(_CACHE_DIR, exist_ok=True) _cache_tag = (f"abp_to_spde_N{N_particles}_L{int(Lx)}" f"_S{Nsteps}_dt{dt_sim}_seed{seed}.npy") _cache_path = os.path.join(_CACHE_DIR, _cache_tag) if os.path.exists(_cache_path): X_micro = np.load(_cache_path) print(f"Loaded cached microscopic trajectory: {X_micro.shape} " f"(from {_cache_path})") else: print(f"Simulating {Nsteps} steps with neighbor rebuild every {rebuild_every} step(s) ...") t0 = time.perf_counter() key, sub = random.split(key) coll_micro = simulate_chunked( proc, dt=dt_sim, Nsteps=Nsteps, key=sub, cutoff=cutoff, box=box_np, skin=skin, rebuild_every=rebuild_every, save_every=10, spatial_dims=slice(0, 2), nnz_safety=3.0, verbose=False, ) sim_time = time.perf_counter() - t0 n_chunks = len(coll_micro.datasets) print(f"Simulation done in {sim_time:.0f}s ({n_chunks} chunks)") # Concatenate frames across chunks: shape (T, N, 3) X_micro = coll_micro.to_array(axis="time") np.save(_cache_path, X_micro) print(f"Cached microscopic trajectory → {_cache_path}") T_total = X_micro.shape[0] print(f"Recorded {T_total} frames at Δt = {dt_sim * 10:.2f}") # Restrict analysis to the **banded transient** window. At early # times the system has not yet polarized; at late times the bands # dissolve into a uniform polar phase. Empirically (see the σ/μ # diagnostic printed below) the propagating bands are sharpest in # the middle third of the trajectory. T_lo = int(0.20 * T_total) T_hi = int(0.60 * T_total) X_use = X_micro[T_lo:T_hi] print(f"Banded window: frames [{T_lo}:{T_hi}] → {X_use.shape[0]} frames") # Quick banding diagnostic on the chosen window: density variance # coefficient σ/μ on a 32×32 grid (>0.3 ↔ visibly banded). _diag_sigs = [] for _ti in range(0, X_use.shape[0], max(1, X_use.shape[0] // 20)): _pos = np.asarray(X_use[_ti, :, :2]) % Lx _H, _, _ = np.histogram2d(_pos[:, 0], _pos[:, 1], bins=32, range=[[0, Lx], [0, Ly]]) _rho = _H / (Lx / 32) ** 2 _diag_sigs.append(_rho.std() / _rho.mean()) _phi_win = float(np.sqrt( np.cos(X_use[..., 2]).mean(axis=1) ** 2 + np.sin(X_use[..., 2]).mean(axis=1) ** 2 ).mean()) print(f"Banding diagnostic on window: ⟨σ_ρ/μ_ρ⟩ = {np.mean(_diag_sigs):.3f} " f"(>0.3 ⇒ banded); polar order ⟨φ⟩ = {_phi_win:.3f}") .. rst-class:: sphx-glr-script-out .. code-block:: none Building initial neighbor list ... nnz = 1554306, ⟨neighbors⟩ = 155.4 (0.24s) Loaded cached microscopic trajectory: (1500, 10000, 3) (from examples/gallery/_cache/abp_to_spde_N10000_L64_S1500_dt0.05_seed7.npy) Recorded 1500 frames at Δt = 0.50 Banded window: frames [300:900] → 600 frames Banding diagnostic on window: ⟨σ_ρ/μ_ρ⟩ = 0.306 (>0.3 ⇒ banded); polar order ⟨φ⟩ = 0.900 .. GENERATED FROM PYTHON SOURCE LINES 287-295 Particle snapshots -------------------- Three snapshots inside the **banded transient window** (frames ``T_lo``, midpoint, ``T_hi``). Heading is HSV-colour-coded; the narrow stripes of common colour, separated by lower-density gaps, are the propagating polar bands characteristic of Toner–Tu / Vicsek flocking. .. GENERATED FROM PYTHON SOURCE LINES 295-319 .. code-block:: Python from SFI.trajectory import TrajectoryCollection # Wrap the (cached or freshly simulated) microscopic trajectory in a # collection so the canonical particle/SPDE plotters can read frames # directly (heading lives in state dimension 2). coll_micro = TrajectoryCollection.from_arrays(X=X_micro, dt=dt_sim * 10) snap_idx = [T_lo, (T_lo + T_hi) // 2, T_hi - 1] fig_snap, axes_snap = plt.subplots(1, 3, figsize=(12, 4.2)) for ax, ti in zip(axes_snap, snap_idx): plot_particles( coll_micro, t_index=ti, color_dim=2, cmap="hsv", vmin=-np.pi, vmax=np.pi, box=box_np, s=2.0, alpha=0.85, edgecolors="none", ax=ax, ) ax.set_xlim(0, Lx); ax.set_ylim(0, Ly) ax.set_title(f"t = {ti * dt_sim * 10:.1f}") fig_snap.suptitle( "Microscopic flock — 10 000 self-propelled particles (colour = heading angle)", fontsize=11, ) plt.show() .. image-sg:: /gallery/images/sphx_glr_abp_to_spde_demo_001.png :alt: Microscopic flock — 10 000 self-propelled particles (colour = heading angle), t = 150.0, t = 300.0, t = 449.5 :srcset: /gallery/images/sphx_glr_abp_to_spde_demo_001.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 320-327 Animated banded flock ----------------------- Sub-sample particles for clarity, then animate every few frames through the banded window. The HSV colour wheel encodes heading angle, so density stripes of a single colour reveal coherently moving polar bands sweeping across the periodic box. .. GENERATED FROM PYTHON SOURCE LINES 327-337 .. code-block:: Python skip = max(1, (T_hi - T_lo) // 200) fig_anim, ax_anim = plt.subplots(figsize=(5.4, 5.2)) anim = animate_particles( coll_micro, color_dim=2, cmap="hsv", vmin=-np.pi, vmax=np.pi, box=box_np, skip=skip, s=2.5, ax=ax_anim, ) plt.show() .. video:: /gallery/images/sphx_glr_abp_to_spde_demo_002.mp4 :class: sphx-glr-single-img :height: 780 :width: 810 :autoplay: :loop: .. GENERATED FROM PYTHON SOURCE LINES 338-346 Coarse-graining: particles → continuum fields ---------------------------------------------- At each frame we deposit ``(1, cosθ_i, sinθ_i)`` for each particle bilinearly onto a periodic 64×64 grid, then convolve with a wrapped Gaussian of standard deviation ``σ = 1.5`` cells. The result is three smooth fields per cell — the local density ``ρ`` and the polar momentum density ``m = ρ\\langle\\hat e_θ\\rangle``. .. GENERATED FROM PYTHON SOURCE LINES 346-377 .. code-block:: Python GRID = (64, 64) DX_GRID = Lx / GRID[0] SIGMA_CELLS = 1.5 t0 = time.perf_counter() fields = coarse_grain_polar( jnp.asarray(X_use), box=box, grid_shape=GRID, sigma_cells=SIGMA_CELLS, angle_index=2, ) print(f"Coarse-graining: {fields.shape} in {time.perf_counter() - t0:.1f}s") # Sanity: integrated density should equal particle count total_mass = float(fields[0, :, 0].sum() * DX_GRID * DX_GRID) print(f"∫ρ dA = {total_mass:.1f} (expected {N_particles})") # Wrap as a TrajectoryCollection on the grid; dt is the frame stride # of the saved microscopic trajectory. from SFI.bases.spde import square_grid_extras from SFI.trajectory import TrajectoryCollection DT_FRAME = dt_sim * 10 # save_every=10 box_extras_grid = square_grid_extras(grid_shape=GRID, dx=DX_GRID) coll_fields = TrajectoryCollection.from_arrays( X=np.asarray(fields, dtype=np.float32), dt=DT_FRAME, extras_global=box_extras_grid, ) print(f"Field trajectory: ({coll_fields.T}, {coll_fields.N}, {coll_fields.d})") .. rst-class:: sphx-glr-script-out .. code-block:: none Coarse-graining: (600, 4096, 3) in 2.3s ∫ρ dA = 10000.0 (expected 10000) Field trajectory: (600, 4096, 3) .. GENERATED FROM PYTHON SOURCE LINES 378-384 Coarse-grained snapshots -------------------------- Density ``ρ`` (heat map) and polar momentum ``m`` (streamlines) at the same physical times as the microscopic snapshots above — minus the burn-in offset. .. GENERATED FROM PYTHON SOURCE LINES 384-403 .. code-block:: Python cg_idx = [0, coll_fields.T // 2, coll_fields.T - 1] vmax_rho = float(np.percentile(fields[..., 0], 99.5)) stream_kw = dict(density=0.9, color="white", linewidth=0.7, arrowsize=0.7) fig_cg, axes_cg = plt.subplots(1, 3, figsize=(12, 4.2)) plot_spde_snapshot( coll_fields, cg_idx, scalar_channel=0, vector_channels=(1, 2), grid_shape=GRID, dx=DX_GRID, render="streamplot", vmin=0.0, vmax=vmax_rho, streamplot_kw=stream_kw, axes=axes_cg, ) for ax, ti in zip(axes_cg, cg_idx): ax.set_xlim(0, Lx); ax.set_ylim(0, Ly); ax.set_aspect("equal") ax.set_title(f"frame {ti}") fig_cg.colorbar(axes_cg[0].images[0], ax=axes_cg.tolist(), shrink=0.78, label=r"density $\rho$") fig_cg.suptitle("Coarse-grained fields: ρ (colour) + m (streamlines)", fontsize=11) plt.show() .. image-sg:: /gallery/images/sphx_glr_abp_to_spde_demo_003.png :alt: Coarse-grained fields: ρ (colour) + m (streamlines), frame 0, frame 300, frame 599 :srcset: /gallery/images/sphx_glr_abp_to_spde_demo_003.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 404-415 Overcomplete SPDE basis: candidate Toner–Tu terms --------------------------------------------------- We assemble ~30 hydrodynamic features with :class:`~SFI.statefunc.layout.GridLayout`: - **ρ-equation** (scalar): mass-flux divergence, gradient pressures, and a few pointwise nonlinearities. - **m-equation** (vector): linear, cubic, and quintic Landau terms; Frank-elastic ``∇²m``; pressure-like ``∇ρ`` and ``∇|m|²``; advective ``(m·∇)m``. .. GENERATED FROM PYTHON SOURCE LINES 415-469 .. code-block:: Python from SFI.statefunc.layout import GridLayout, ScalarSector, VectorSector from SFI.statefunc.structexpr import StructuredExpr layout = GridLayout( rho=ScalarSector([0]), m=VectorSector([1, 2], sdim=2, spatial=True), dim=3, ndim=2, bc="pbc", ) rho = layout.rho m = layout.m # Squared magnitude |m|² (scalar) via einsum contraction m_sq = StructuredExpr.einsum("i,i->", m, m).with_label("|m|²") # --- ρ-equation candidates (scalar output per cell) --- rho_terms = [ layout.div(m).with_label("∇·m"), layout.div(rho * m).with_label("∇·(ρm)"), layout.lap(rho), (layout.lap(rho * rho)).with_label("∇²ρ²"), layout.lap(m_sq).with_label("∇²|m|²"), layout.advection_by(m, rho), rho.with_label("ρ"), (rho * rho).with_label("ρ²"), m_sq, ] rho_basis = rho_terms[0] for t in rho_terms[1:]: rho_basis = rho_basis & t # --- m-equation candidates (vector output per cell) --- m_terms = [ m.with_label("m"), (m_sq * m).with_label("|m|²m"), ((m_sq * m_sq) * m).with_label("|m|⁴m"), layout.lap(m), layout.grad(layout.div(m)).with_label("∇(∇·m)"), layout.grad(rho).with_label("∇ρ"), (layout.grad(rho * rho)).with_label("∇ρ²"), layout.grad(m_sq).with_label("∇|m|²"), layout.advection_by(m, m), (m * layout.div(m)).with_label("m(∇·m)"), (rho * m).with_label("ρm"), ] m_basis = m_terms[0] for t in m_terms[1:]: m_basis = m_basis & t BASIS = layout.embed(rank=1, rho=rho_basis, m=m_basis) n_feat = BASIS.n_features print(f"Overcomplete basis: {n_feat} candidate features " f"({len(rho_terms)} for ρ + {len(m_terms)} for m)") .. rst-class:: sphx-glr-script-out .. code-block:: none Overcomplete basis: 20 candidate features (9 for ρ + 11 for m) .. GENERATED FROM PYTHON SOURCE LINES 470-472 Linear inference + PASTIS sparsification ------------------------------------------ .. GENERATED FROM PYTHON SOURCE LINES 472-488 .. code-block:: Python from SFI import OverdampedLangevinInference inf = OverdampedLangevinInference(coll_fields) inf.compute_diffusion_constant(method="WeakNoise") inf.infer_force_linear(BASIS, M_mode="Ito") inf.compute_force_error() inf.sparsify_force(criterion="PASTIS", p=0.001) inf.compute_force_error() k_sel, support_sel, _, coeffs_sel = \ inf.force_sparsity_result.select_by_ic("PASTIS") inf.print_report() .. rst-class:: sphx-glr-script-out .. code-block:: none --- StochasticForceInference Report --- Average diffusion tensor: [[ 1.7269813e-05 6.2673080e-06 -1.4041916e-05] [ 6.2673080e-06 4.1185002e-04 1.4832322e-04] [-1.4041916e-05 1.4832322e-04 1.3275006e-04]] Measurement noise tensor: [[-2.8831907e-06 -8.2338357e-07 9.7881082e-07] [-8.2338369e-07 2.5146157e-05 1.8341721e-05] [ 9.7881082e-07 1.8341721e-05 -1.5236627e-05]] Force estimated information: 684290.3125 Force: estimated normalized mean squared error (sampling only): 1.4613569012135748e-05 Force Coefficient Table ───────────────────────────────────────────────────────────────── # Label Coefficient Std.Err SNR Sig ───────────────────────────────────────────────────────────────── 0 ∇·m -6.42026e-02 2.15407e-04 298.1 *** 1 ∇·(ρm) 8.27530e-03 9.87390e-05 83.8 ** 2 ∇²rho 6.85043e-02 8.50566e-04 80.5 ** 3 ∇²ρ² 4.55006e-02 2.37074e-04 191.9 *** 4 ∇²|m|² -2.23570e-03 1.50703e-04 14.8 ** 5 (m·∇)rho -1.19461e-02 1.16922e-04 102.2 *** 6 ρ 1.30153e-03 1.75887e-05 74.0 ** 7 ρ² -3.03218e-05 8.01770e-06 3.8 * 8 |m|² -4.79972e-04 1.16139e-05 41.3 ** 9 m 3.65620e-03 6.22864e-05 58.7 ** 10 |m|²m -4.71209e-04 1.20343e-05 39.2 ** 11 |m|⁴m 3.54635e-05 9.29253e-07 38.2 ** 12 ∇²m 3.24034e-01 7.49543e-04 432.3 *** 13 ∇(∇·m) -2.26157e-01 1.15357e-03 196.1 *** 14 ∇ρ -8.33872e-02 7.32461e-04 113.8 *** 15 ∇ρ² 9.73189e-03 1.14047e-04 85.3 ** 16 ∇|m|² -2.23048e-02 1.29393e-04 172.4 *** 17 (m·∇)m 4.63086e-02 9.44763e-05 490.2 *** 18 m(∇·m) -3.17987e-03 1.03356e-04 30.8 ** 19 ρm -5.32718e-04 4.02552e-05 13.2 ** ───────────────────────────────────────────────────────────────── 20/20 basis functions in support, sig: 20* / 19** / 8*** (|SNR| ≥ 2 / 10 / 100) (Std.err. reflects sampling error only; discretization bias is not included.) .. GENERATED FROM PYTHON SOURCE LINES 489-500 Inferred SPDE: term coefficients ---------------------------------- We sort the inferred basis by the magnitude of its coefficient. The few largest terms — Frank elasticity ``∇²m``, the density Laplacian / divergence-of-flux pair, advective ``(m·∇)m`` and the pressure-like ``∇ρ`` — are precisely the canonical Toner–Tu hydrodynamic terms. Sub-dominant nonlinearities span four orders of magnitude below them; they remain *statistically* significant at this enormous sample size (1050 frames × 4096 cells), so PASTIS keeps them, but they are dynamically negligible. .. GENERATED FROM PYTHON SOURCE LINES 500-523 .. code-block:: Python labels = list(BASIS.labels) if BASIS.labels else [ f"f{i}" for i in range(n_feat) ] if len(labels) != n_feat: labels = [f"f{i}" for i in range(n_feat)] # Magnitudes of the PASTIS-selected coefficients, sorted descending on a # log axis; ``show_pruned`` appends faded zero-bars for the rejected # candidate terms so the full overcomplete basis stays visible. fig_bar, ax_bar = plt.subplots(figsize=(13, 4.6)) plot_recovery_bar( np.abs(np.asarray(coeffs_sel)), support_sel, labels=labels, yscale="log", sort=True, show_pruned=True, ax=ax_bar, ) ax_bar.set_ylabel(r"$|\,\mathrm{coefficient}\,|$ (log)") ax_bar.set_title( f"Inferred SPDE — coefficients sorted by magnitude " f"(PASTIS keeps {k_sel}/{n_feat})" ) fig_bar.tight_layout() plt.show() .. image-sg:: /gallery/images/sphx_glr_abp_to_spde_demo_004.png :alt: Inferred SPDE — coefficients sorted by magnitude (PASTIS keeps 20/20) :srcset: /gallery/images/sphx_glr_abp_to_spde_demo_004.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-script-out .. code-block:: none examples/gallery/abp_to_spde_demo.py:520: UserWarning: The figure layout has changed to tight fig_bar.tight_layout() .. GENERATED FROM PYTHON SOURCE LINES 524-531 Bootstrap: simulate the *inferred* SPDE ----------------------------------------- Initialise from random ρ, m fluctuations and integrate the sparse inferred SPDE. If the dynamics are correct, propagating bands should re-emerge spontaneously, just like in the agent-based simulation. .. GENERATED FROM PYTHON SOURCE LINES 531-539 .. code-block:: Python key, bkey = random.split(key) coll_boot, _ = inf.simulate_bootstrapped_trajectory( bkey, oversampling=4, simulate=True, ) T_boot = coll_boot.T print(f"Bootstrap SPDE: {T_boot} frames") .. rst-class:: sphx-glr-script-out .. code-block:: none Bootstrap SPDE: 600 frames .. GENERATED FROM PYTHON SOURCE LINES 540-546 Bootstrap visualisation ------------------------- Side-by-side: agent-based bands (left, coarse-grained) vs the bands generated by the inferred SPDE (right) — both at their respective late times. .. GENERATED FROM PYTHON SOURCE LINES 546-576 .. code-block:: Python fig_bs, axes_bs = plt.subplots(1, 2, figsize=(9, 4.6)) vmax_rho = float(np.percentile(fields[..., 0], 99.5)) stream_kw = dict(density=0.9, color="white", linewidth=0.7, arrowsize=0.7) # Microscopic CG (last frame) vs inferred SPDE (last frame), both as # ρ heat-map + m streamlines. plot_spde_snapshot( coll_fields, coll_fields.T - 1, scalar_channel=0, vector_channels=(1, 2), grid_shape=GRID, dx=DX_GRID, render="streamplot", vmin=0.0, vmax=vmax_rho, streamplot_kw=stream_kw, axes=axes_bs[0], ) axes_bs[0].set_title("Agent-based (coarse-grained)") axes_bs[0].set_aspect("equal") plot_spde_snapshot( coll_boot, T_boot - 1, scalar_channel=0, vector_channels=(1, 2), grid_shape=GRID, dx=DX_GRID, render="streamplot", vmin=0.0, vmax=vmax_rho, streamplot_kw=stream_kw, axes=axes_bs[1], ) axes_bs[1].set_title("Inferred SPDE (PASTIS-sparse)") axes_bs[1].set_aspect("equal") fig_bs.suptitle( "Bootstrap test — inferred SPDE reproduces the agent-based flock", fontsize=11, ) plt.show() .. image-sg:: /gallery/images/sphx_glr_abp_to_spde_demo_005.png :alt: Bootstrap test — inferred SPDE reproduces the agent-based flock, Agent-based (coarse-grained), Inferred SPDE (PASTIS-sparse) :srcset: /gallery/images/sphx_glr_abp_to_spde_demo_005.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 577-586 Side-by-side animation: agent-based vs inferred SPDE ------------------------------------------------------ The climax: animate the coarse-grained agent-based fields (left) next to the freely-running inferred SPDE (right). Both initialised independently (the SPDE from random fluctuations), and both develop propagating polar bands of comparable amplitude and wavelength — evidence that PASTIS has recovered the *dynamics*, not just the instantaneous structure. .. GENERATED FROM PYTHON SOURCE LINES 586-600 .. code-block:: Python T_anim = min(coll_fields.T, T_boot) n_anim_frames = 150 stride = max(1, T_anim // n_anim_frames) vmax_rho = float(np.percentile(fields[..., 0], 99.5)) anim_dual = animate_spde_comparison( coll_fields, coll_boot, grid_shape=GRID, field_component=0, skip=stride, vmin=0.0, vmax=vmax_rho, interval=60, titles=("Agent-based (coarse-grained)", "Inferred SPDE (PASTIS-sparse)"), ) plt.show() .. video:: /gallery/images/sphx_glr_abp_to_spde_demo_006.mp4 :class: sphx-glr-single-img :height: 675 :width: 1350 :autoplay: :loop: .. rst-class:: sphx-glr-script-out .. code-block:: none [Generated: 2026-06-30 12:57] .. rst-class:: sphx-glr-timing **Total running time of the script:** (6 minutes 27.927 seconds) .. rst-class:: sphx-glr-example-tags 🏷 Tags: synthetic, overdamped, multi-particle, linear, spde, experimental, pastis, interactions .. _sphx_glr_download_gallery_abp_to_spde_demo.py: .. only:: html .. container:: sphx-glr-footer sphx-glr-footer-example .. container:: sphx-glr-download sphx-glr-download-jupyter :download:`Download Jupyter notebook: abp_to_spde_demo.ipynb ` .. container:: sphx-glr-download sphx-glr-download-python :download:`Download Python source code: abp_to_spde_demo.py ` .. container:: sphx-glr-download sphx-glr-download-zip :download:`Download zipped: abp_to_spde_demo.zip ` .. only:: html .. rst-class:: sphx-glr-signature `Gallery generated by Sphinx-Gallery `_