qlroms.charts

Local charts: clustering, per-cluster POD, and the Chart building block.

Implements section 1 of docs/streamlined_vision_building_blocks.md: - fit_clusters: K-means++ clustering ("full" | "minibatch" | "warmstart"), optional overlapping assignment with a tolerance. - clustering_features: which representation of the snapshots k-means sees ("physical" snapshots, lossless global POD, reduced global POD, or a custom callable). - compute_pod_basis: local basis construction (snapshot POD; exact or randomized SVD; Mw-weighted method of snapshots when a mass matrix is given). - Chart: one cluster's physical <-> reduced round trip under the inner product = u^T Mw v (Chart.from_fenics + mass_to_torch_sparse adapt dolfinx-side POD output without qlroms ever importing dolfinx/petsc4py).

The Atlas container and the one-call build (fit_charts) live in qlroms.atlas.

Chart

One cluster's local chart: the physical(N-dim) <-> reduced(r-dim) round trip for a SINGLE chart (its own Phi, its own centroid), under the inner product = u^T Mw v -- Mw is None for the plain M = I case (default), a (N,) diagonal (e.g. ks2d's quadrature weights), or a dense/sparse (N, N) matrix (e.g. fenics' assembled FEM mass matrix). Phi is assumed Mw-orthonormal, so projection is a = Phi^T Mw (u - c) and recovery is c + Phi a.

This is both the building block Atlas combines pairwise into the cached tmap/tshift, and the parent of every case single-cluster ROM class (ks1d.rom.ROM, ks2d.rom.ROM, pinball.rom.ROM): those add the cluster's dynamics operators + step_reduced on top of this geometry.

Methods only read self.Phi / self.centroid / self.Mw / self.device / self.rdtype, so dataclass subclasses that define those as fields inherit them without calling init here.

Source code in qlroms/charts.py

class Chart:
    """One cluster's local chart: the physical(N-dim) <-> reduced(r-dim) round trip
    for a SINGLE chart (its own Phi, its own centroid), under the inner product
    <u, v> = u^T Mw v -- Mw is None for the plain M = I case (default), a (N,)
    diagonal (e.g. ks2d's quadrature weights), or a dense/sparse (N, N) matrix
    (e.g. fenics' assembled FEM mass matrix). Phi is assumed Mw-orthonormal, so
    projection is a = Phi^T Mw (u - c) and recovery is c + Phi a.

    This is both the building block Atlas combines pairwise into the cached
    tmap/tshift, and the parent of every case single-cluster ROM class
    (ks1d.rom.ROM, ks2d.rom.ROM, pinball.rom.ROM): those add
    the cluster's dynamics operators + `step_reduced` on top of this geometry.

    Methods only read self.Phi / self.centroid / self.Mw / self.device /
    self.rdtype, so dataclass subclasses that define those as fields inherit them
    without calling __init__ here."""

    Mw = None   # default inner-product weight (identity); fields/instances may override

    def __init__(self, Phi: torch.Tensor, centroid: torch.Tensor, r: int | None = None,
                 Mw: torch.Tensor | None = None):
        self.Phi = Phi                    # (N, r)
        self.centroid = centroid          # (N,)
        self.r = int(Phi.shape[1]) if r is None else r
        self.Mw = Mw
        self.device, self.rdtype = Phi.device, Phi.dtype

    def _phi_cent(self) -> tuple[torch.Tensor, torch.Tensor]:
        """Phi and column centroid on the configured device/dtype (no-op for plain
        Charts; dataclass subclasses may carry config device != tensor device)."""
        phi = self.Phi.to(dtype=self.rdtype, device=self.device)
        cent = self.centroid.reshape(-1, 1).to(dtype=self.rdtype, device=self.device)
        return phi, cent

    def project_state(self, state) -> torch.Tensor:
        """(N,) or (N, m) physical -> (r, m) reduced: a = Phi^T Mw (u - c)."""
        s = torch.as_tensor(state, dtype=self.rdtype, device=self.device)
        if s.ndim == 1:
            s = s[:, None]
        phi, cent = self._phi_cent()
        return phi.T @ apply_weight(self.Mw, s - cent)

    def recover_state(self, a) -> torch.Tensor:
        """(r,) reduced -> (N,), or (r, m) -> (N, m) physical."""
        a = torch.as_tensor(a, dtype=self.rdtype, device=self.device)
        squeezed = a.ndim == 1
        if squeezed:
            a = a[:, None]
        phi, cent = self._phi_cent()
        out = cent + phi @ a
        return out.squeeze(-1) if squeezed else out

    @classmethod
    def from_fenics(cls, modes, mean, mass=None) -> "Chart":
        """Build a Chart from dolfinx-side POD output (duck-typed: only
        Function.x.array and Mat.getValuesCSR/getSize are read, so qlroms never
        imports dolfinx or petsc4py). An atlas of these is
        qlroms.atlas.Atlas.from_charts(charts) -- Atlas is Mw-aware.

        Args:
            modes: sequence of dolfinx Functions -- the cluster's mass-orthonormal
                POD modes (columns of Phi).
            mean: dolfinx Function -- the cluster centroid (e.g. velocity mean).
            mass: assembled PETSc mass matrix, an already-converted torch tensor,
                or None (M = I).
        """
        Phi = torch.from_numpy(np.stack([np.asarray(m.x.array, dtype=float).copy()
                                         for m in modes], axis=1))
        centroid = torch.from_numpy(np.asarray(mean.x.array, dtype=float).copy())
        if mass is not None and not torch.is_tensor(mass):
            mass = mass_to_torch_sparse(mass)
        return cls(Phi, centroid, Mw=mass)

project_state(state)

(N,) or (N, m) physical -> (r, m) reduced: a = Phi^T Mw (u - c).

Source code in qlroms/charts.py

def project_state(self, state) -> torch.Tensor:
    """(N,) or (N, m) physical -> (r, m) reduced: a = Phi^T Mw (u - c)."""
    s = torch.as_tensor(state, dtype=self.rdtype, device=self.device)
    if s.ndim == 1:
        s = s[:, None]
    phi, cent = self._phi_cent()
    return phi.T @ apply_weight(self.Mw, s - cent)

recover_state(a)

(r,) reduced -> (N,), or (r, m) -> (N, m) physical.

Source code in qlroms/charts.py

def recover_state(self, a) -> torch.Tensor:
    """(r,) reduced -> (N,), or (r, m) -> (N, m) physical."""
    a = torch.as_tensor(a, dtype=self.rdtype, device=self.device)
    squeezed = a.ndim == 1
    if squeezed:
        a = a[:, None]
    phi, cent = self._phi_cent()
    out = cent + phi @ a
    return out.squeeze(-1) if squeezed else out

from_fenics(modes, mean, mass=None) classmethod

Build a Chart from dolfinx-side POD output (duck-typed: only Function.x.array and Mat.getValuesCSR/getSize are read, so qlroms never imports dolfinx or petsc4py). An atlas of these is qlroms.atlas.Atlas.from_charts(charts) -- Atlas is Mw-aware.

Parameters:

Name	Type	Description	Default
modes		sequence of dolfinx Functions -- the cluster's mass-orthonormal POD modes (columns of Phi).	required
mean		dolfinx Function -- the cluster centroid (e.g. velocity mean).	required
mass		assembled PETSc mass matrix, an already-converted torch tensor, or None (M = I).	None

Source code in qlroms/charts.py

@classmethod
def from_fenics(cls, modes, mean, mass=None) -> "Chart":
    """Build a Chart from dolfinx-side POD output (duck-typed: only
    Function.x.array and Mat.getValuesCSR/getSize are read, so qlroms never
    imports dolfinx or petsc4py). An atlas of these is
    qlroms.atlas.Atlas.from_charts(charts) -- Atlas is Mw-aware.

    Args:
        modes: sequence of dolfinx Functions -- the cluster's mass-orthonormal
            POD modes (columns of Phi).
        mean: dolfinx Function -- the cluster centroid (e.g. velocity mean).
        mass: assembled PETSc mass matrix, an already-converted torch tensor,
            or None (M = I).
    """
    Phi = torch.from_numpy(np.stack([np.asarray(m.x.array, dtype=float).copy()
                                     for m in modes], axis=1))
    centroid = torch.from_numpy(np.asarray(mean.x.array, dtype=float).copy())
    if mass is not None and not torch.is_tensor(mass):
        mass = mass_to_torch_sparse(mass)
    return cls(Phi, centroid, Mw=mass)

QuadROM

Bases: Chart

The family-neutral quadratic single-cluster ROM member: chart geometry plus reduced (b, A, B). Shared by BOTH model families -- once the operators are reduced, intrusive (Galerkin-projected) and non-intrusive (OpInf-fitted) members are the same compute, the ODE da/dt = b + A a + B(a, a) stepped by ETDRK4 (or, discrete=True, the one-step map applied directly; B may be None for a linear/affine map). Where the operators CAME from is recorded by the thin family subclasses (qlroms.data_driven_qlroms.opinf.OpInfROM, qlroms.intrusive_qlroms.galerkin.GalerkinROM) and by the compilation class (qlOpinf vs qlGalerkin); this class carries no provenance. Operators are None until seeded via set_operators.

Source code in qlroms/charts.py

class QuadROM(Chart):
    """The family-neutral quadratic single-cluster ROM member: chart geometry plus
    reduced (b, A, B). Shared by BOTH model families -- once the operators are
    reduced, intrusive (Galerkin-projected) and non-intrusive (OpInf-fitted) members
    are the same compute, the ODE da/dt = b + A a + B(a, a) stepped by ETDRK4 (or,
    `discrete=True`, the one-step map applied directly; B may be None for a
    linear/affine map). Where the operators CAME from is recorded by the thin family
    subclasses (qlroms.data_driven_qlroms.opinf.OpInfROM, qlroms.intrusive_qlroms.galerkin.GalerkinROM) and
    by the compilation class (qlOpinf vs qlGalerkin); this class carries no
    provenance. Operators are None until seeded via set_operators."""

    discrete = False   # class-level default so pre-`discrete` pickled caches still load

    def __init__(self, Phi, centroid, dt: float, Mw=None,
                 Mcontour: int = 32, Rcontour: float = 15.0, discrete: bool = False):
        super().__init__(Phi, centroid, Mw=Mw)
        self.dt = float(dt)   # for discrete maps this is metadata only
        self.Mcontour, self.Rcontour = int(Mcontour), float(Rcontour)
        self.discrete = bool(discrete)
        self.b = self.A = self.B = None

    def set_operators(self, b, A, B) -> None:
        self.b, self.A, self.B = b, A, B
        self.__dict__.pop("etdrk4_rom", None)   # invalidate the cached coefficients

    @cached_property
    def etdrk4_rom(self) -> tuple:
        """(E, E2, Q, f1, f2, f3) ETDRK4 coefficients for the linear part A."""
        from .operators import build_etdrk4_coeffs
        if self.A is None:
            raise RuntimeError("Operators not set; call set_operators (or the owning "
                               "compilation's fit) first.")
        coeffs = build_etdrk4_coeffs(self.A.detach().cpu().numpy(), self.dt,
                                     M=self.Mcontour, R=self.Rcontour)
        return tuple(torch.from_numpy(c).to(dtype=self.rdtype, device=self.device) for c in coeffs)

    def step_reduced(self, a: torch.Tensor) -> torch.Tensor:
        """One step of the pure (r, 1) reduced state: a direct map application if
        discrete, an ETDRK4 integration of the reduced ODE otherwise."""
        from .operators import quadratic_etdrk4_step, quadratic_map_step
        if self.A is None:
            raise RuntimeError("Operators not set; call set_operators (or the owning "
                               "compilation's fit) first.")
        if self.discrete:
            return quadratic_map_step(a, self.b, self.A, self.B)
        return quadratic_etdrk4_step(a, self.b, self.B, self.etdrk4_rom)

etdrk4_rom cached property

(E, E2, Q, f1, f2, f3) ETDRK4 coefficients for the linear part A.

step_reduced(a)

One step of the pure (r, 1) reduced state: a direct map application if discrete, an ETDRK4 integration of the reduced ODE otherwise.

Source code in qlroms/charts.py

def step_reduced(self, a: torch.Tensor) -> torch.Tensor:
    """One step of the pure (r, 1) reduced state: a direct map application if
    discrete, an ETDRK4 integration of the reduced ODE otherwise."""
    from .operators import quadratic_etdrk4_step, quadratic_map_step
    if self.A is None:
        raise RuntimeError("Operators not set; call set_operators (or the owning "
                           "compilation's fit) first.")
    if self.discrete:
        return quadratic_map_step(a, self.b, self.A, self.B)
    return quadratic_etdrk4_step(a, self.b, self.B, self.etdrk4_rom)

fit_clusters(Xtrain, K, random_state=1, kmeans_method='minibatch', kmeans_n_init=KMEANS_N_INIT, kmeans_max_iter=KMEANS_MAX_ITER, assign_overlapping=False, overlap_tolerance=1.1)

Cluster training snapshots with KMeans.

Parameters:

Name	Type	Description	Default
Xtrain	Tensor	Snapshot matrix (Ndof, Ntrain). Real or complex; only real part is used.	required
K	int	Number of clusters.	required
random_state	int	KMeans random seed.	1
assign_overlapping	bool	Whether to assign points to multiple clusters if they are equidistant.	False
overlap_tolerance	float	Tolerance for overlapping assignments.	1.1

kmeans_method selects the fitting strategy, all k-means++ seeded: "full": KMeans from scratch, n_init random restarts. "minibatch": MiniBatchKMeans only -- fast, approximate (the default). "warmstart": MiniBatchKMeans first for cheap centroids, then one full KMeans pass (n_init=1) initialized from them for full-quality convergence.

Returns:

Name	Type	Description
centroids	Tensor	(K, Ndof) tensor on the same device/dtype as Xtrain.
labels	Tensor	(Ntrain,) long tensor of cluster assignments.
cluster_sizes	list[int]	list of int cluster populations.
Xtrain_augment	Tensor	Augmented training data tensor. (same as Xtrain if assign_overlapping=False, else repeated columns for overlapping assignments)
aug_idx	Tensor	(Naug,) long tensor mapping each augmented column back to its original Xtrain column (arange(Ntrain) if assign_overlapping=False) -- lets callers re-express the augmentation in another space (fit_charts clusters in feature space but needs physical snapshots for the POD).

Source code in qlroms/charts.py

def fit_clusters(
    Xtrain: torch.Tensor,
    K: int,
    random_state: int = 1,
    kmeans_method: KMEANS_METHOD = "minibatch",
    kmeans_n_init: int = KMEANS_N_INIT,
    kmeans_max_iter: int = KMEANS_MAX_ITER,
    assign_overlapping: bool = False,
    overlap_tolerance: float = 1.1,
) -> tuple[torch.Tensor, torch.Tensor, list[int], torch.Tensor, torch.Tensor]:
    """Cluster training snapshots with KMeans.

    Args:
        Xtrain: Snapshot matrix (Ndof, Ntrain). Real or complex; only real part is used.
        K: Number of clusters.
        random_state: KMeans random seed.
        assign_overlapping: Whether to assign points to multiple clusters if they are equidistant.
        overlap_tolerance: Tolerance for overlapping assignments.

    kmeans_method selects the fitting strategy, all k-means++ seeded:
        "full":      KMeans from scratch, n_init random restarts.
        "minibatch": MiniBatchKMeans only -- fast, approximate (the default).
        "warmstart": MiniBatchKMeans first for cheap centroids, then one full KMeans
                     pass (n_init=1) initialized from them for full-quality convergence.

    Returns:
        centroids: (K, Ndof) tensor on the same device/dtype as Xtrain.
        labels: (Ntrain,) long tensor of cluster assignments.
        cluster_sizes: list of int cluster populations.
        Xtrain_augment: Augmented training data tensor. (same as Xtrain if assign_overlapping=False, else repeated columns for overlapping assignments)
        aug_idx: (Naug,) long tensor mapping each augmented column back to its
            original Xtrain column (arange(Ntrain) if assign_overlapping=False) --
            lets callers re-express the augmentation in another space (fit_charts
            clusters in feature space but needs physical snapshots for the POD).
    """
    if K > 1:
        xtrain_np = Xtrain.real.cpu().numpy().T

        if kmeans_method == "minibatch":
            kmeans = MiniBatchKMeans(
                n_clusters=K,
                n_init=kmeans_n_init,
                max_iter=kmeans_max_iter,
                random_state=random_state,
            )
        else:
            if kmeans_method == "warmstart":
                mbk = MiniBatchKMeans(
                    n_clusters=K,
                    n_init=kmeans_n_init,
                    max_iter=kmeans_max_iter,
                    random_state=random_state,
                )
                mbk.fit(xtrain_np)
                init = mbk.cluster_centers_
                n_init = 1  # already seeded from a converged solution; restarts would waste it
            else:
                init = "k-means++"
                n_init = kmeans_n_init
            kmeans = KMeans(
                n_clusters=K,
                init=init,
                n_init=n_init,
                max_iter=kmeans_max_iter,
                random_state=random_state,
                algorithm="elkan",
            )

        kmeans.fit(xtrain_np)
        centers_np = kmeans.cluster_centers_
        if assign_overlapping:
            distances = kmeans.transform(xtrain_np)  # shape (Ntrain, K)
            # assign points to clusters within overlap_tolerance x nearest distance;
            # multiplicative form (vs. distances/min_distances <= tolerance) avoids
            # a 0/0 when a point sits exactly on its nearest centroid.
            min_distances = distances.min(axis=1, keepdims=True)
            threshold = overlap_tolerance * min_distances
            labels_np = [
                [k for k in range(K) if distances[i, k] <= threshold[i]]
                for i in range(xtrain_np.shape[0])
            ]
            aug_idx_np = np.repeat(np.arange(len(labels_np)), [len(ks) for ks in labels_np])
            Xtrain_augment = Xtrain[:, torch.from_numpy(aug_idx_np).to(Xtrain.device)]
            labels_np = np.array([k for sublist in labels_np for k in sublist])  # flatten the list of lists

        else:
            labels_np = kmeans.labels_
            aug_idx_np = np.arange(Xtrain.shape[1])
            Xtrain_augment = Xtrain

        labels = torch.from_numpy(labels_np).to(device=Xtrain.device, dtype=torch.long)
        aug_idx = torch.from_numpy(aug_idx_np).to(device=Xtrain.device, dtype=torch.long)
        centroids = torch.from_numpy(centers_np).to(
            device=Xtrain.device, dtype=Xtrain.real.dtype
        )
    else:
        labels = torch.zeros(Xtrain.shape[1], dtype=torch.long, device=Xtrain.device)
        centroids = Xtrain.real.mean(dim=1, keepdim=True).T.contiguous()
        Xtrain_augment = Xtrain
        aug_idx = torch.arange(Xtrain.shape[1], device=Xtrain.device)

    cluster_sizes = []
    for k in range(K):
        nk = int((labels == k).sum().item())
        if nk < 2:
            raise ValueError(f"Cluster {k} has only {nk} snapshots.")
        cluster_sizes.append(nk)

    return centroids, labels, cluster_sizes, Xtrain_augment, aug_idx

compute_pod_basis(X, r, method='exact', oversampling=10, n_iter=2, block_size=256, random_state=0, return_singular_values=False, Mw=None)

Leading r left singular vectors (POD spatial modes) of an already-centered column-snapshot matrix X (N, n).

Parameters:

Name	Type	Description	Default
Mw		inner-product weight (None = identity). When given (a (N,) diagonal or a dense/sparse (N, N) mass matrix), the basis is computed by the method of snapshots under = u^T Mw v -- the returned modes are Mw-orthonormal (Phi^T Mw Phi = I), as a FEM mass-matrix POD requires. method and the randomized knobs are ignored in that case (the (n, n) snapshot Gram is cheap whenever n << N, the regime where a weighted POD matters).	None
method	POD_METHODS	"exact": full economy SVD (torch.linalg.svd) -- exact, O(N * n * min(N, n)). "randomized": torch.svd_lowrank, PyTorch's native randomized range-finder + power iteration -- the same approach ks2d/config.py already uses for its larger spectral grids; O(N * n * (r + oversampling)). "randomized_blocked": manual block-processed randomized SVD (see _randomized_pod_basis_blocked); same asymptotic cost as "randomized" but with an explicit, inspectable oversampling/n_iter/block_size and no dependency on torch's internal implementation.	'exact'
oversampling/n_iter		extra sampling dimension and power-iteration count for the two randomized methods (ignored for "exact").	required
block_size	int	"randomized_blocked" only -- columns processed per matmul chunk.	256
random_state	int	seed for the randomized methods' Gaussian sketch (ignored for "exact").	0
return_singular_values	bool	if True, also return the corresponding leading r singular values -- e.g. for a POD spectrum plot without a second, separate (and possibly inconsistent) SVD call.	False

Returns: Phi: (N, r) tensor (r may be less than requested if min(N, n) < r), or (Phi, svals) if return_singular_values is True.

Source code in qlroms/charts.py

def compute_pod_basis(
    X: torch.Tensor,
    r: int,
    method: POD_METHODS = "exact",
    oversampling: int = 10,
    n_iter: int = 2,
    block_size: int = 256,
    random_state: int = 0,
    return_singular_values: bool = False,
    Mw=None,
) -> torch.Tensor | tuple[torch.Tensor, torch.Tensor]:
    """Leading r left singular vectors (POD spatial modes) of an already-centered
    column-snapshot matrix X (N, n).

    Args:
        Mw: inner-product weight (None = identity). When given (a (N,) diagonal or a
            dense/sparse (N, N) mass matrix), the basis is computed by the method of
            snapshots under <u, v> = u^T Mw v -- the returned modes are Mw-orthonormal
            (Phi^T Mw Phi = I), as a FEM mass-matrix POD requires. `method` and the
            randomized knobs are ignored in that case (the (n, n) snapshot Gram is
            cheap whenever n << N, the regime where a weighted POD matters).
        method:
            "exact": full economy SVD (torch.linalg.svd) -- exact, O(N * n * min(N, n)).
            "randomized": torch.svd_lowrank, PyTorch's native randomized range-finder +
                power iteration -- the same approach ks2d/config.py already uses for its
                larger spectral grids; O(N * n * (r + oversampling)).
            "randomized_blocked": manual block-processed randomized SVD (see
                _randomized_pod_basis_blocked); same asymptotic cost as "randomized" but
                with an explicit, inspectable oversampling/n_iter/block_size and no
                dependency on torch's internal implementation.
        oversampling/n_iter: extra sampling dimension and power-iteration count for the
            two randomized methods (ignored for "exact").
        block_size: "randomized_blocked" only -- columns processed per matmul chunk.
        random_state: seed for the randomized methods' Gaussian sketch (ignored for
            "exact").
        return_singular_values: if True, also return the corresponding leading r
            singular values -- e.g. for a POD spectrum plot without a second, separate
            (and possibly inconsistent) SVD call.
    Returns:
        Phi: (N, r) tensor (r may be less than requested if min(N, n) < r), or
        (Phi, svals) if return_singular_values is True.
    """
    N, n = X.shape
    r_eff = min(r, n, N)
    if Mw is not None and not torch.is_tensor(Mw):
        Mw = torch.as_tensor(Mw, dtype=X.dtype, device=X.device)
    if Mw is not None and Mw.ndim == 1:
        # Diagonal weight: reduce to an unweighted SVD of sqrt(w) X (honoring
        # `method`), then unscale -- Phi = Q / sqrt(w) satisfies Phi^T diag(w)
        # Phi = I. O(N n r) instead of the (n, n) snapshot Gram below, which is
        # the wrong regime when n >> N (e.g. many snapshots per cluster).
        sqrtw = torch.sqrt(Mw)[:, None]
        out = compute_pod_basis(sqrtw * X, r, method=method, oversampling=oversampling,
                                n_iter=n_iter, block_size=block_size,
                                random_state=random_state,
                                return_singular_values=return_singular_values, Mw=None)
        if return_singular_values:
            return out[0] / sqrtw, out[1]
        return out / sqrtw
    if Mw is not None:
        # Method of snapshots under the Mw inner product: C = X^T Mw X (n, n),
        # modes = X V / sqrt(lambda) are Mw-orthonormal.
        C = X.T @ apply_weight(Mw, X)
        evals, evecs = torch.linalg.eigh(C)                  # ascending
        evals, evecs = evals.flip(0), evecs.flip(1)          # descending
        positive = evals > 0
        r_eff = min(r_eff, int(positive.sum().item()))
        Phi = X @ (evecs[:, :r_eff] / torch.sqrt(evals[:r_eff]))
        svals = torch.sqrt(evals[:r_eff])
        return (Phi, svals) if return_singular_values else Phi
    if method == "exact":
        U, S, _ = torch.linalg.svd(X, full_matrices=False)
        Phi, svals = U[:, :r_eff], S[:r_eff]
    elif method == "randomized":
        torch.manual_seed(random_state)  # torch.svd_lowrank has no generator argument
        q = min(r_eff + oversampling, N, n)
        Uk, Sk, _ = torch.svd_lowrank(X, q=q, niter=n_iter)
        Phi, svals = Uk[:, :r_eff], Sk[:r_eff]
    elif method == "randomized_blocked":
        Phi, svals = _randomized_pod_basis_blocked(
            X, r_eff, oversampling=oversampling, n_iter=n_iter,
            block_size=block_size, random_state=random_state,
        )
    else:
        raise ValueError(
            f"Unknown POD method '{method}'; choose 'exact', 'randomized', or 'randomized_blocked'."
        )
    return (Phi, svals) if return_singular_values else Phi

clustering_features(X, r=None, cluster_space='physical', pod_method='exact', random_state=1, **pod_kwargs)

Representation of the snapshots that k-means actually clusters on.

The clustering feature options. Only the distance metric used to assign labels changes; fit_charts always recomputes physical-space centroids afterwards when a non-physical space is used.

Parameters:

Name	Type	Description	Default
X	Tensor	(Ndof, Ntrain) snapshot matrix, columns = snapshots.	required
cluster_space	CLUSTER_SPACES	"physical": cluster directly on the raw snapshots -- correct but the most expensive k-means call when Ndof >> Ntrain. "pod_lossless": cluster on coefficients in a full-rank (rank = min(Ndof, Ntrain - 1)) global POD basis -- an exact, distance- preserving rotation of physical space, so the partition is identical to "physical" (up to k-means++'s random init) but far cheaper when Ndof >> Ntrain. "pod_r": cluster on coefficients in the same r-mode (lossy) truncation the local ROMs will use -- cheapest, but the partition can genuinely differ since truncated energy is invisible to k-means. callable: any custom feature map X -> (Nfeat, Ntrain).	'physical'
r	int | None	number of modes for "pod_r" (required for that space).	None
pod_method/pod_kwargs		forwarded to compute_pod_basis for "pod_r".	required

Returns: (Nfeat, Ntrain) feature matrix whose COLUMNS k-means clusters.

Source code in qlroms/charts.py

def clustering_features(
    X: torch.Tensor,
    r: int | None = None,
    cluster_space: CLUSTER_SPACES = "physical",
    pod_method: POD_METHODS = "exact",
    random_state: int = 1,
    **pod_kwargs,
) -> torch.Tensor:
    """Representation of the snapshots that k-means actually clusters on.

    The clustering feature options. Only the distance metric
    used to assign labels changes; fit_charts always recomputes physical-space
    centroids afterwards when a non-physical space is used.

    Args:
        X: (Ndof, Ntrain) snapshot matrix, columns = snapshots.
        cluster_space:
            "physical": cluster directly on the raw snapshots -- correct but the most
                expensive k-means call when Ndof >> Ntrain.
            "pod_lossless": cluster on coefficients in a full-rank (rank =
                min(Ndof, Ntrain - 1)) global POD basis -- an exact, distance-
                preserving rotation of physical space, so the partition is identical
                to "physical" (up to k-means++'s random init) but far cheaper when
                Ndof >> Ntrain.
            "pod_r": cluster on coefficients in the same r-mode (lossy) truncation the
                local ROMs will use -- cheapest, but the partition can genuinely differ
                since truncated energy is invisible to k-means.
            callable: any custom feature map X -> (Nfeat, Ntrain).
        r: number of modes for "pod_r" (required for that space).
        pod_method/pod_kwargs: forwarded to compute_pod_basis for "pod_r".
    Returns:
        (Nfeat, Ntrain) feature matrix whose COLUMNS k-means clusters.
    """
    if callable(cluster_space):
        return cluster_space(X)
    if cluster_space == "physical":
        return X
    Xc = X - X.mean(dim=1, keepdim=True)
    if cluster_space == "pod_lossless":
        r_lossless = min(X.shape[0], X.shape[1] - 1)
        Phi_c = compute_pod_basis(Xc, r_lossless, method="exact")
    elif cluster_space == "pod_r":
        if r is None:
            raise ValueError("cluster_space='pod_r' needs r.")
        Phi_c = compute_pod_basis(Xc, r, method=pod_method, random_state=random_state, **pod_kwargs)
    else:
        raise ValueError(
            f"Unknown cluster_space '{cluster_space}'; choose 'physical', 'pod_lossless', "
            f"'pod_r', or pass a callable."
        )
    return Phi_c.T @ Xc

apply_weight(Mw, x)

M @ x for the chart inner-product weight: identity when Mw is None, elementwise for a (N,) diagonal, matmul for a dense/sparse (N, N) matrix.

Source code in qlroms/charts.py

def apply_weight(Mw, x: torch.Tensor) -> torch.Tensor:
    """M @ x for the chart inner-product weight: identity when Mw is None,
    elementwise for a (N,) diagonal, matmul for a dense/sparse (N, N) matrix."""
    if Mw is None:
        return x
    if Mw.ndim == 1:
        return Mw[:, None] * x
    return torch.sparse.mm(Mw, x) if Mw.is_sparse else Mw @ x

mass_orthonormalize(X, Mw=None, tol=1e-12)

Modified Gram-Schmidt of the columns of X (N, m) under = u^T Mw v.

Columns whose residual Mw-norm falls below tol are dropped (near-linear dependence), so the result may have fewer columns than X. Used e.g. to re-orthonormalize a POD basis after supremizer enrichment (fenics intrusive build); the output satisfies Q^T Mw Q = I.

Source code in qlroms/charts.py

def mass_orthonormalize(X: torch.Tensor, Mw=None, tol: float = 1e-12) -> torch.Tensor:
    """Modified Gram-Schmidt of the columns of X (N, m) under <u, v> = u^T Mw v.

    Columns whose residual Mw-norm falls below `tol` are dropped (near-linear
    dependence), so the result may have fewer columns than X. Used e.g. to
    re-orthonormalize a POD basis after supremizer enrichment (fenics intrusive
    build); the output satisfies Q^T Mw Q = I."""
    cols: list[torch.Tensor] = []
    for j in range(X.shape[1]):
        v = X[:, [j]].clone()
        for q in cols:
            v = v - q @ (q.T @ apply_weight(Mw, v))
        norm = torch.sqrt((v.T @ apply_weight(Mw, v)).clamp_min(0.0)).item()
        if norm < tol:
            continue
        cols.append(v / norm)
    return torch.cat(cols, dim=1) if cols else X[:, :0]

mass_to_torch_sparse(mass)

PETSc AIJ mass matrix -> torch sparse COO (the chart inner-product weight Mw). Duck-typed (getValuesCSR / getSize only), so qlroms never imports petsc4py.

Source code in qlroms/charts.py

def mass_to_torch_sparse(mass) -> torch.Tensor:
    """PETSc AIJ mass matrix -> torch sparse COO (the chart inner-product weight Mw).
    Duck-typed (getValuesCSR / getSize only), so qlroms never imports petsc4py."""
    ia, ja, av = mass.getValuesCSR()
    n = mass.getSize()[0]
    return torch.sparse_csr_tensor(
        torch.from_numpy(ia.astype(np.int64)), torch.from_numpy(ja.astype(np.int64)),
        torch.from_numpy(np.asarray(av, dtype=float)), size=(n, n), dtype=torch.float64,
    ).to_sparse_coo()