qlroms.operators

Common operator layer for all qlROM model families.

Every model family exposes the same forecast form on augmented chart-local states apod (r+1, 1) with the cluster id in the last row:

apod^{n+1} = model.step(apod^n)          # discrete map, all families

so the outer forecast workflow never depends on model-family internals. Continuous families (qlGalerkin) realize step() by an ETDRK4 integration of their da/dt = g(a).

This module owns the case-agnostic operator machinery: the contour-integration ETDRK4 coefficient builder (build_etdrk4_coeffs) and the dimension-agnostic quadratic reduced-space steppers (quadratic_etdrk4_step, quadratic_map_step).

build_etdrk4_coeffs(L, dt, M=32, R=15.0)

Unified ETDRK4 coefficients builder via contour integration.

Computes (E, E2, Q, f1, f2, f3) from an unscaled linear operator L via scipy.linalg.expm for exponentials and contour integration for resolvent integrals. Works for both reduced (r x r) and full-order (Nx x Nx) operators.

Parameters:

Name	Type	Description	Default
L		Unscaled linear operator matrix (n x n), typically a spectral matrix or reduced operator.	required
dt		Time step size (used to scale L internally).	required
M		Number of contour quadrature points (default 32).	32
R		Contour radius for resolvent integral (default 15.0).	15.0

Returns:

Type	Description
	Tuple (E, E2, Q, f1, f2, f3), each (n x n) real numpy arrays.

Source code in qlroms/operators.py

def build_etdrk4_coeffs(L, dt, M=32, R=15.0):
    """Unified ETDRK4 coefficients builder via contour integration.

    Computes (E, E2, Q, f1, f2, f3) from an unscaled linear operator L via
    scipy.linalg.expm for exponentials and contour integration for resolvent integrals.
    Works for both reduced (r x r) and full-order (Nx x Nx) operators.

    Args:
        L: Unscaled linear operator matrix (n x n), typically a spectral matrix or reduced operator.
        dt: Time step size (used to scale L internally).
        M: Number of contour quadrature points (default 32).
        R: Contour radius for resolvent integral (default 15.0).

    Returns:
        Tuple (E, E2, Q, f1, f2, f3), each (n x n) real numpy arrays.
    """
    # Convert to numpy if needed
    L_np = _to_numpy(L)
    rdim = L_np.shape[0]

    # Scale by dt to get the matrix exponential argument
    A = dt * L_np.astype(complex)

    # Exponentials via scipy
    E = expm(A)
    E2 = expm(A / 2.0)

    # Resolvent integrals via contour integration
    I = np.eye(rdim, dtype=complex)
    Q = np.zeros((rdim, rdim), dtype=complex)
    f1 = np.zeros((rdim, rdim), dtype=complex)
    f2 = np.zeros((rdim, rdim), dtype=complex)
    f3 = np.zeros((rdim, rdim), dtype=complex)

    m = np.arange(1, M + 1)
    Rvals = R * np.exp(1j * np.pi * (m - 0.5) / M)

    for z in Rvals:
        zIA = np.linalg.inv(z * I - A)
        ez2 = np.exp(z / 2.0)
        ez = np.exp(z)

        Q += dt * zIA * (ez2 - 1.0)
        f1 += dt * zIA * ((-4.0 - z + ez * (4.0 - 3.0 * z + z * z)) / (z * z))
        f2 += dt * zIA * ((2.0 + z + ez * (z - 2.0)) / (z * z))
        f3 += dt * zIA * ((-4.0 - 3.0 * z - z * z + ez * (4.0 - z)) / (z * z))

    Q = np.real(Q / M)
    f1 = np.real(f1 / M)
    f2 = np.real(f2 / M)
    f3 = np.real(f3 / M)
    E = np.real(E)
    E2 = np.real(E2)

    return E, E2, Q, f1, f2, f3

quadratic_etdrk4_step(a, b, B, etdrk4_coeffs)

One ETDRK4 step of da/dt = b + B(a, a) (A already folded into etdrk4_coeffs).

Dimension-agnostic: only ever reads (r,)-sized tensors, so the same function steps quadratic 1D and 2D clusters. etdrk4_coeffs are the (E, E2, Q, f1, f2, f3) coefficients built from the linear operator A (see qlROMBase.get_stepping_coeffs); b, B are that cluster's affine/quadratic parts. a is a pure (r, 1) reduced-coordinate column (no cluster-id row).

Source code in qlroms/operators.py

def quadratic_etdrk4_step(a: torch.Tensor, b: torch.Tensor, B: torch.Tensor, etdrk4_coeffs: tuple) -> torch.Tensor:
    """One ETDRK4 step of  da/dt = b + B(a, a)  (A already folded into etdrk4_coeffs).

    Dimension-agnostic: only ever reads (r,)-sized tensors, so the same function
    steps quadratic 1D and 2D clusters. `etdrk4_coeffs` are the (E, E2, Q, f1,
    f2, f3) coefficients built from the linear operator A (see
    qlROMBase.get_stepping_coeffs); b, B are that cluster's affine/quadratic
    parts. `a` is a pure (r, 1) reduced-coordinate column (no cluster-id row).
    """
    E, E2, Q, f1, f2, f3 = etdrk4_coeffs

    def nfun(av: torch.Tensor) -> torch.Tensor:
        return b[:, None] + torch.einsum("ijk,jm,km->im", B, av, av)

    Nu = nfun(a)
    aa = E2 @ a + Q @ Nu
    Na = nfun(aa)
    bb = E2 @ a + Q @ Na
    Nb = nfun(bb)
    cc = E2 @ aa + Q @ (2.0 * Nb - Nu)
    Nc = nfun(cc)
    return E @ a + f1 @ Nu + 2.0 * f2 @ (Na + Nb) + f3 @ Nc

quadratic_map_step(a, b, A, B=None)

One DISCRETE step of the map a^{n+1} = b + A a^n + B(a^n, a^n).

The direct one-step analog of quadratic_etdrk4_step (no ODE, no integrator): the operators ARE the map, applied as-is (QuadROM with discrete=True). a is a pure (r, 1) reduced-coordinate column. B=None (or empty) gives the affine/linear map b + A a -- classic DMD.

Source code in qlroms/operators.py

def quadratic_map_step(a: torch.Tensor, b: torch.Tensor, A: torch.Tensor,
                       B: torch.Tensor | None = None) -> torch.Tensor:
    """One DISCRETE step of the map  a^{n+1} = b + A a^n + B(a^n, a^n).

    The direct one-step analog of quadratic_etdrk4_step (no ODE, no integrator): the
    operators ARE the map, applied as-is (QuadROM with discrete=True). `a` is a pure
    (r, 1) reduced-coordinate column. B=None (or empty) gives the affine/linear map
    b + A a -- classic DMD.
    """
    out = b[:, None] + A @ a
    if B is not None and B.numel():
        out = out + torch.einsum("ijk,jm,km->im", B, a, a)
    return out