semiflow-py

PyO3 Python bindings for semiflow — Chernoff approximations of operator semigroups (Remizov 2025, Theorem 6).

Built on the semiflow core crate (ADR-0154, 2026-06-10). The Python surface has parity with all core kernel families via ADR-0111 Waves P1–P7 plus the 0.9.0-beta addition of ReverseHeat1D (reverse-mode AD, math §51, ADR-0156): 26 binding classes + 1 free function. Pyright errors: 0. Complete __init__.pyi stubs; py.typed marker; GIL released in all evolve paths (ADR-0031).

Tensor-train and gridless carriers (TtEvolver, TtState, TtCoupledEvolver, VarCoefTtEvolver, GridlessEvolver, MeasureState) are fully exposed via PyO3 as of 0.12.0-beta (Rust-only at 0.9.0-beta; bindings added in ADR-0171 / ADR-0178).

Installation

pip install semiflow-pde

Note: semiflow-pde is published on PyPI. Wheels for common platforms are also available via GitHub releases if you need a specific build: pip install semiflow_pde-*.whl.

Or build from source (requires Rust toolchain + maturin):

pip install maturin
maturin develop --profile release-ffi -m crates/semiflow-py/Cargo.toml

Array I/O conventions

• All real-valued state arrays are numpy.float64 (np.float64).
• Schrödinger and SchrodingerComplex1D state arrays are numpy.complex128.
• 2D state is flat float64 in row-major x-fastest order: index j*nx + i corresponds to u(x_i, y_j).
• 3D state is flat x-fastest: index k*nx*ny + j*nx + i.
• values() always returns a copy of the internal Rust state; mutations to the returned array do not affect the object.
• Inputs are validated for NaN/Inf at construction and before evolve; non-finite inputs raise SemiflowError(kind='NanInf').
• All finite-check and grid-size errors raise SemiflowError.

Error model

All semiflow-py operations raise a single exception type:

from semiflow import SemiflowError

The .kind attribute (a string) identifies the error category:

Boundary policies

All 1D/2D/3D kernels accept a keyword argument boundary:

Far-field / Dirichlet-like boundaries for finance

Users coming from classical finite-difference option pricers often look for an inhomogeneous Dirichlet far-field ("set V = S − K e^{−rτ} at S_max") and, finding only "zero", conclude the library cannot price calls. In fact boundary="linear" is the correct and recommended closure for European call pricing.

Why "linear" works for calls: The Black-Scholes / CEV solution satisfies V(S, τ) ≈ S − K e^{−rτ} as S → ∞. This asymptote is linear in S, so V_SS → 0. Linear extrapolation from the two outermost grid nodes reproduces an affine far-field exactly (to machine precision), introducing no spurious curvature or kink. Validated on Shift1D.with_arrays with S ∈ [0, 4K], n = 1025: ATM relative error ≈ 8.5e-5 with no boundary artifact.

Puts in log-price space: When a European put is priced in log-price coordinates (z = ln S), the solution satisfies u_z → 0 at both ends of a sufficiently wide domain. Either "reflect" (zero-flux Neumann) or "linear" works; "reflect" is marginally preferable because it imposes the correct zero-derivative condition exactly.

Example — Black-Scholes call via Shift1D:

import numpy as np
import semiflow as rp

# Black-Scholes PDE: dV/dt = (sigma^2 / 2) S^2 V_SS + r S V_S - r V
# Rewrite as L V = a(S) V_SS + b(S) V_S + c(S) V with:
#   a(S) = 0.5 * sigma^2 * S^2,  b(S) = r * S,  c(S) = -r

K, T, r, sigma = 100.0, 1.0, 0.05, 0.20
S_max = 4.0 * K      # wide enough that far-field is linear
n = 1025
S_nodes = np.linspace(0.0, S_max, n)

a_arr = 0.5 * sigma**2 * S_nodes**2
b_arr = r * S_nodes
c_arr = np.full(n, -r)

# European call payoff at maturity
u0 = np.maximum(S_nodes - K, 0.0)

state = rp.Shift1D.with_arrays(
    0.0, S_max, n,
    a_arr, b_arr, c_arr,
    c_norm_bound=r,  # upper bound on |c(x)|
    u0=u0,
    boundary="linear",   # <-- correct far-field closure for calls
)
state.evolve(t=T, n_steps=200)

V = state.values()
S_atm_idx = np.argmin(np.abs(S_nodes - K))
print(f"ATM call price (semiflow): {V[S_atm_idx]:.6f}")
# rel. error vs Black-Scholes closed form ≈ 8.5e-5 at n=1025, n_steps=200

Inhomogeneous Dirichlet (future work): The Rust core provides BoundaryPolicy::Dirichlet { value } (a constant ghost-node extension), but it is not yet wired into the Python boundary= string parser. If you need to pin u = g_lo at the left edge and u = g_hi at the right edge with distinct values, use "linear" as the near-exact idiom for asymptotically-linear payoffs. Support for boundary=("dirichlet", value) is tracked in issue #20.

Usage examples

1. Unit-diffusion 1D heat

Solve ∂_t u = ∂²_x u on [-10, 10] with a Gaussian initial condition:

import numpy as np
import semiflow as rp

n = 1000
xs = np.linspace(-10.0, 10.0, n)
u0 = np.exp(-(xs - 0.5)**2 / 0.01)   # narrow Gaussian at x=0.5

state = rp.Heat1D(-10.0, 10.0, n, u0)
state.evolve(t=1.0, n_steps=100)

u = state.values()    # float64 ndarray, shape (n,)
print(f"n={len(state)}, max={u.max():.6f}")

The GIL is released during evolve (ADR-0031); concurrent Python threads make progress during long calls.

2. SchrodingerComplex1D — native complex128 wavefunction

Solve i ψ_t = (−½∂²_x + V) ψ and verify unitarity:

import numpy as np
import semiflow as rp

n = 512
xs = np.linspace(-10.0, 10.0, n)
psi0 = np.exp(-xs**2 / 2.0).astype(np.complex128)  # normalised Gaussian
psi0 /= np.sqrt(np.trapz(np.abs(psi0)**2, xs))      # L2-normalise

sch = rp.SchrodingerComplex1D(-10.0, 10.0, n, psi0)
norm0 = sch.norm_squared()

sch.evolve(t=0.5, n_steps=200)

psi_t = sch.values()    # complex128 ndarray
assert abs(sch.norm_squared() / norm0 - 1.0) < 1e-12, "unitarity violated"
print(f"norm ratio = {sch.norm_squared() / norm0:.15f}")

3. Manifold2D — Riemannian manifold heat kernel

Solve ∂_t u = Δ_{S²} u on the 2-sphere via MMRS 2023 Chernoff formula:

import numpy as np
import semiflow as rp

nx, ny = 32, 64
u0 = np.zeros(nx * ny, dtype=np.float64)
u0[nx * (ny // 2) + nx // 2] = 1.0  # delta-like at chart centre

sphere = rp.Manifold2D(
    0.1, np.pi - 0.1, nx,    # theta axis
    0.0, 2 * np.pi,   ny,    # phi axis
    u0,
    manifold="sphere2",
    radius=1.0,
    curvature_correction=True,  # enables R/12 correction -> order 2
)
sphere.evolve(t=0.02, n_steps=50)

u_t = sphere.values()    # float64 ndarray, length nx*ny (row-major theta-fastest)
print(f"integral ≈ {u_t.sum() * (np.pi / nx) * (2 * np.pi / ny):.4f}")

Available manifolds: "torus" (flat T²), "sphere2" (S²(r)), "hyperbolic2" (Poincaré disk H²(s)). The radius parameter sets r or s.

Class reference

Classes are grouped by kernel family. All stateful classes expose at least evolve(t, n_steps=100) (mutates in-place, GIL released) and values() → NDArray[np.float64] (copy). See __init__.pyi for complete signatures.

1D diffusion family

Operator splitting — multi-dimensional

Schrödinger

Both support .with_potential(v_array) and .norm_squared().

Boundary-condition kernels

Time-dependent and subordinated

Geometry and hypoelliptic operators

Graph PDE

Matrix and point-eval kernels

Anisotropic multi-D

Adjoint and adaptive wrappers

Reverse-mode AD (0.9.0-beta, ADR-0156)

Constructor parameters:

.value_and_grad(tau, u0, target) -> (float, float):

import numpy as np
import semiflow as rp

n_grid = 24
xs = np.linspace(-4.0, 4.0, n_grid)

rc = rp.ReverseHeat1D(theta=0.4, xmin=-4.0, xmax=4.0, n_grid=n_grid, n_steps=8)
u0     = np.exp(-xs**2)
target = np.zeros(n_grid)

value, grad = rc.value_and_grad(tau=0.05, u0=u0, target=target)
print(f"loss={value:.6e}  ∂J/∂θ={grad:.6e}")

Raises SemiflowError with .kind in {'OutOfDomain', 'GridMismatch', 'NanInf'}.

NARROW scope (§51.5): constant-a DiffusionChernoff only; θ is the uniform diffusivity. Variable-coefficient and nonlinear kernels are out of scope for ReverseHeat1D.

Tensor-train carriers (ADR-0171 / ADR-0178)

Curse-escaped O(d·n·r²) storage for separable diagonal-A diffusion on ℝᵈ. State and evolvers share the TtState carrier; TtEvolver / TtCoupledEvolver cover constant-coefficient cases, VarCoefTtEvolver covers variable-coefficient additive-separable operators (issue #2, ADR-0178).

Gridless / particle carriers (ADR-0171)

Sparse signed-measure ensemble on ℝ (D=1). Curse-escaped: the 3ᴰ dense tree is never materialised; only sparse marginals and scalar observables cross the Python boundary.

v3 Evolver surface

Greeks and hyper-dual AD

Adjoint Fokker-Planck (particle pushforward)

Killing and absorbing boundary conditions

Obstacle / variational-inequality family

Resolvent time-jump family (TWS contour)

Matrix diffusion (coupled 2-component)

Wentzell dynamic boundary conditions

Additional 1D kernels

Sparse-grid and high-dimensional kernels

Graph adjoint and Krylov families

Symmetric operator and FEM utilities

Performance

GIL release follows the three-phase py.detach pattern (ADR-0031): acquire → snapshot inputs → detach → Rust compute → reacquire. Send + Sync is verified at compile time with static_assertions.

Indicative timings on i7-12700K (1000 nodes, 100 steps, Heat1D):

For large grids or many time steps, prefer .with_a_array over .with_a_function: the array path uses a pure-Rust Arc<Vec<f64>> Catmull-Rom interpolant and never re-acquires the GIL during evolve.

Type stubs

__init__.pyi and the py.typed marker ship with every wheel. Static type checkers (mypy, pyright, pylance) pick them up automatically. The pyrightconfig.json at the repo root adds crates/semiflow-py/python to extraPaths so local development also resolves the stubs correctly (0 reportAttributeAccessIssue errors).

Mathematical reference

I. D. Remizov, Vladikavkaz Math. J. 27(4) (2025) 124–135. DOI 10.46698/a3908-1212-5385-q

Changelog / Release notes

• CHANGELOG.md — full version history
• GitHub Releases — tagged release notes and wheel downloads

License

MIT OR Apache-2.0 — same as semiflow.