#
# ************************************************************************
#
# ROM Tools and Workflows
# Copyright 2019 National Technology & Engineering Solutions of Sandia,LLC
# (NTESS)
#
# Under the terms of Contract DE-NA0003525 with NTESS, the
# U.S. Government retains certain rights in this software.
#
# ROM Tools and Workflows is licensed under BSD-3-Clause terms of use:
#
# Redistribution and use in source and binary forms, with or without
# modification, are permitted provided that the following conditions
# are met:
#
# 1. Redistributions of source code must retain the above copyright
# notice, this list of conditions and the following disclaimer.
#
# 2. Redistributions in binary form must reproduce the above copyright
# notice, this list of conditions and the following disclaimer in the
# documentation and/or other materials provided with the distribution.
#
# 3. Neither the name of the copyright holder nor the names of its
# contributors may be used to endorse or promote products derived
# from this software without specific prior written permission.
#
# THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
# "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
# LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS
# FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE
# COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT,
# INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
# (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
# SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
# HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT,
# STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING
# IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
# POSSIBILITY OF SUCH DAMAGE.
#
# Questions? Contact Eric Parish (ejparis@sandia.gov)
#
# ************************************************************************
#
'''
This module defines the API to work with a vector subspace.
A vector subspace is foundational to reduced-order models.
In a ROM, a high-dimensional state is restricted to live within a low-dimensional vector space, known as a trial space.
Mathematically, given a "FOM" vector :math:`\\mathbf{u} \\in \\mathbb{R}^{N_{\\mathrm{vars}} N_{\\mathrm{x}}}`,
we can write
.. math::
\\mathbf{u} \\approx \\tilde{\\mathbf{u}} \\in \\mathcal{V} + \\mathbf{u}_{\\mathrm{shift}}
where
- :math:`\\mathcal{V}` with :math:`\\text{dim}(\\mathcal{V}) = K \\le N_{\\mathrm{vars}} N_{\\mathrm{x}}`
is the trial space
- :math:`N_{\\mathrm{vars}}` is the number of PDE variables (e.g., 5 for the compressible Navier-Stokes equations in 3D)
- :math:`N_{\\mathrm{x}}` is the number of spatial DOFs
Formally, we can describe this low-dimensional representation with a basis and an affine offset,
.. math::
\\tilde{\\mathbf{u}} = \\boldsymbol \\Phi \\hat{\\mathbf{u}} + \\mathbf{u}_{\\mathrm{shift}}
where :math:`\\boldsymbol \\Phi \\in \\mathbb{R}^{ N_{\\mathrm{vars}} N_{\\mathrm{x}} \\times K}` is the basis matrix
(:math:`K` is the number of basis), :math:`\\hat{\\mathbf{u}} \\in \\mathbb{R}^{K}` are the reduced, or generalized
coordinates, :math:`\\mathbf{u}_{\\mathrm{shift}} \\in \\mathbb{R}^{ N_{\\mathrm{vars}} N_{\\mathrm{x}}}` is the
shift vector (or affine offset), and, by definition, :math:`\\mathcal{V} \\equiv \\mathrm{range}(\\boldsymbol \\Phi)`.
The `VectorSpace` abstract class defined below encapsulates the information of an affine vector space,
:math:`\\mathcal{V}`, by virtue of providing access to a basis matrix, a shift vector, and the dimensionality of the
vector space, while decoupling this representation from *how* it is computed.
Tensor representation
---------------------
Our representation of the basis and the affine offset for a vector space is based on tensors
.. math::
\\mathcal{\\Phi} \\in \\mathbb{R}^{ N_{\\mathrm{vars}} \\times N_{\\mathrm{x}} \\times K}
.. math::
\\mathcal{u}_{\\mathrm{shift}} \\in \\mathbb{R}^{ N_{\\mathrm{vars}} \\times N_{\\mathrm{x}}}
Internally, we remark that all tensors are reshaped into 2D matrices, e.g., when performing SVD.
Content
-------
We currently provide the following concrete classes:
- `DictionaryVectorSpace`: construct a vector space from a matrix without truncation.
- `VectorSpaceFromPOD`: construct a vector subspace computed via SVD.
- `VectorSpaceFromStreamingPOD`: construct a POD vector space from snapshot blocks.
which derive from the abstract class `VectorSpace`.
API
---
'''
from typing import Tuple, Protocol, Callable
import numpy as np
from romtools.vector_space.utils.truncater import LeftSingularVectorTruncater, NoOpTruncater
from romtools.vector_space.utils.shifter import (
Shifter, StreamingShifter, create_noop_shifter
)
from romtools.vector_space.utils.scaler import Scaler, StreamingScaler, NoOpScaler
from romtools.vector_space.utils.orthogonalizer import Orthogonalizer, NoOpOrthogonalizer
from romtools.vector_space.utils.snapshot_loader import SnapshotLoader
import romtools.linalg.linalg as la
[docs]
class VectorSpace(Protocol):
'''
Abstract base class for vector space implementations.
Methods:
'''
[docs]
def extents(self) -> Tuple[int, int, int]:
'''
Retrieves the dimension of the vector space
Returns:
A Tuple with the dimensions of the vector space (n_var,nx,K).
'''
...
[docs]
def get_shift_vector(self) -> np.ndarray:
'''
Retrieves the shift vector of the vector space.
Returns:
`np.ndarray`: The shift vector in tensor form.
'''
...
[docs]
def get_basis(self) -> np.ndarray:
'''
Retrieves the basis vectors of the vector space.
Returns:
`np.ndarray`: The basis of the vector space in tensor form.
'''
...
[docs]
class DictionaryVectorSpace():
'''
Reduced basis vector space (no truncation).
This class conforms to `VectorSpace` protocol.
Given a snapshot matrix :math:`\\mathbf{S}`, we set the basis to be
.. math::
\\boldsymbol \\Phi = \\mathrm{orthogonalize}(\\mathbf{S} - \\mathbf{u}_{\\mathrm{shift}})
where the orthogonalization and shifts are defined by their
respective classes
'''
def __init__(self,
snapshots,
shifter: Shifter = None,
orthogonalizer: Orthogonalizer = NoOpOrthogonalizer()) -> None:
'''
Constructor.
Args:
snapshots (np.ndarray): Snapshot data in tensor form
:math:`\\in \\mathbb{R}^{ N_{\\mathrm{vars}} \\times N_{\\mathrm{x}} \\times N_{samples}}`
shifter: Class that shifts the basis.
orthogonalizer: Class that orthogonalizes the basis.
This constructor initializes a vector space by performing basis
manipulation operations on the provided snapshot data.
'''
# Create noop shifter if not provided
if shifter is None:
shifter = create_noop_shifter(snapshots)
# compute basis
n_var = snapshots.shape[0]
shifter.apply_shift(snapshots)
snapshot_matrix = _tensor_to_matrix(snapshots)
self.__basis = orthogonalizer.orthogonalize(snapshot_matrix)
self.__basis = _matrix_to_tensor(n_var, self.__basis)
self.__shift_vector = shifter.get_shift_vector()
[docs]
def get_shift_vector(self) -> np.ndarray:
'''
Concrete implementation of `VectorSpace.get_shift_vector()`
'''
return self.__shift_vector
[docs]
def get_basis(self) -> np.ndarray:
'''
Concrete implementation of `VectorSpace.get_basis()`
'''
return self.__basis
def extents(self) -> Tuple[int, int, int]:
return self.__basis.shape
[docs]
class VectorSpaceFromPOD():
'''
POD vector space (constructed via SVD).
This class conforms to `VectorSpace` protocol.
Given a snapshot matrix :math:`\\mathbf{S}`, we compute the basis :math:`\\boldsymbol \\Phi` as
.. math::
\\boldsymbol U = \\mathrm{SVD}(\\mathrm{prescale}(\\mathbf{S} - \\mathbf{u}_{\\mathrm{shift}}))
.. math::
\\boldsymbol \\Phi = \\mathrm{orthogonalize}(\\mathrm{postscale}(\\mathrm{truncate}( \\boldsymbol U )))
where :math:`\\boldsymbol U` are the left singular vectors and the
orthogonalization, truncation, scaling, and shifts are defined by their respective classes.
For truncation, we enable truncation based on a fixed dimension or the
decay of singular values; please refer to the documentation for the truncater.
'''
def __init__(self,
snapshots,
truncater: LeftSingularVectorTruncater = NoOpTruncater(),
shifter: Shifter = None,
orthogonalizer: Orthogonalizer = NoOpOrthogonalizer(),
scaler: Scaler = NoOpScaler(),
svdFnc: Callable = None) -> None:
'''
Constructor.
Args:
snapshots (np.ndarray): Snapshot data in tensor form
:math:`\\in \\mathbb{R}^{ N_{\\mathrm{vars}} \\times N_{\\mathrm{x}} \\times N_{samples}}`
truncater (Truncater): Concrete implementation for truncating the basis.
shifter (Shifter): Concrete implementation responsible for shifting the basis.
orthogonalizer (Orthogonalizer): Concrete implementation that orthogonalizes the basis.
scaler: Concrete implementation that scales the basis.
svdFnc: a callable to use for computing the SVD on the snapshots data.
IMPORTANT: must conform to the API of
[np.linalg.svd](https://numpy.org/doc/stable/reference/generated/numpy.linalg.svd.html#numpy-linalg-svd).
If `None`, internally we use `np.linalg.svd`.
Note: this is useful when you want to use a custom svd, for
example when your snapshots are distributed with MPI, or
maybe you have a fancy svd function that you can use.
This constructor initializes a POD vector space by performing SVD on the provided snapshot
data and applying various basis manipulation operations, including truncation, shifting, scaling,
and orthogonalization.
'''
if shifter is None:
shifter = create_noop_shifter(snapshots)
n_var = snapshots.shape[0]
shifter.apply_shift(snapshots)
scaler.pre_scale(snapshots)
snapshot_matrix = _tensor_to_matrix(snapshots)
svd_picked = np.linalg.svd if svdFnc is None else svdFnc
lsv, svals, _ = svd_picked(snapshot_matrix, full_matrices=False,
compute_uv=True, hermitian=False)
self.__svals = svals
self.__basis = truncater.truncate(lsv, svals)
self.__basis = _matrix_to_tensor(n_var, self.__basis)
scaler.post_scale(self.__basis)
self.__basis = _tensor_to_matrix(self.__basis)
self.__basis = orthogonalizer.orthogonalize(self.__basis)
self.__basis = _matrix_to_tensor(n_var, self.__basis)
self.__shift_vector = shifter.get_shift_vector()
[docs]
def get_singular_values(self) -> np.ndarray:
'''
Returns array of singular values
'''
return self.__svals
[docs]
def get_shift_vector(self) -> np.ndarray:
'''
Concrete implementation of `VectorSpace.get_shift_vector()`
'''
return self.__shift_vector
[docs]
def get_basis(self) -> np.ndarray:
'''
Concrete implementation of `VectorSpace.get_basis()`
'''
return self.__basis
def extents(self) -> Tuple[int, int, int]:
return self.__basis.shape
[docs]
class VectorSpaceFromStreamingPOD():
'''
POD vector space constructed from snapshot blocks loaded on demand.
This class conforms to the VectorSpace protocol. Snapshot data is loaded
twice by the streaming POD algorithm, and only one block is held at a time.
Data-derived shifters and scalers perform initialization before those two
passes. The randomized algorithm first computes up to
``max_basis_dimension`` candidate modes, after which the supplied truncater
selects the retained basis. Energy truncation is measured against the
complete shifted and pre-scaled snapshot matrix.
'''
def __init__(self,
snapshot_loader: SnapshotLoader,
block_size: int,
n_snapshots: int,
max_basis_dimension: int,
truncater: LeftSingularVectorTruncater = NoOpTruncater(),
shifter: StreamingShifter = None,
orthogonalizer: Orthogonalizer = NoOpOrthogonalizer(),
scaler: StreamingScaler = NoOpScaler(),
svdFnc: Callable = None,
comm=None) -> None:
'''
Constructor.
Args:
snapshot_loader: Snapshot loader returning a
tensor of shape (n_var, n_dofs, end-start).
block_size: Maximum number of snapshots loaded at once.
n_snapshots: Total number of snapshots available to the loader.
max_basis_dimension: Number of randomized candidate POD modes.
truncater: Truncater selecting the final basis from the candidate
modes.
shifter: Optional shifter supporting loader-based initialization.
orthogonalizer: Optional final basis orthogonalizer.
scaler: Scaler supporting streaming initialization.
svdFnc: NumPy-compatible range-sketch SVD callable.
comm: Optional communicator for row-distributed streaming POD.
'''
if shifter is not None:
shifter.initialize_shift_vector_from_loader(
snapshot_loader, block_size, n_snapshots, comm
)
def shifted_snapshot_loader(start: int, end: int) -> np.ndarray:
snapshot_block = np.array(snapshot_loader(start, end), copy=True)
if shifter is not None:
shifter.apply_shift(snapshot_block)
return snapshot_block
scaler.initialize_scalings_from_loader(
shifted_snapshot_loader, block_size, n_snapshots, comm
)
def shifted_and_scaled_snapshot_loader(
start: int, end: int) -> np.ndarray:
snapshot_block = shifted_snapshot_loader(start, end)
scaler.pre_scale(snapshot_block)
return snapshot_block
basis, svals, _, total_energy = la.streaming_pod(
snapshot_loader=shifted_and_scaled_snapshot_loader,
block_size=block_size,
n_snapshots=n_snapshots,
max_basis_dimension=max_basis_dimension,
svdFnc=svdFnc,
comm=comm,
)
if basis.ndim != 3:
raise ValueError(
"VectorSpaceFromStreamingPOD requires three-dimensional snapshot blocks"
)
self.__svals = svals
n_var = basis.shape[0]
basis_matrix = _tensor_to_matrix(basis)
basis_matrix = truncater.truncate(
basis_matrix, svals, total_energy=total_energy
)
basis = _matrix_to_tensor(n_var, basis_matrix)
scaler.post_scale(basis)
basis_matrix = _tensor_to_matrix(basis)
basis_matrix = orthogonalizer.orthogonalize(basis_matrix)
self.__basis = _matrix_to_tensor(n_var, basis_matrix)
if shifter is None:
self.__shift_vector = np.zeros(self.__basis.shape[:2])
else:
self.__shift_vector = shifter.get_shift_vector()
[docs]
def get_singular_values(self) -> np.ndarray:
'''Returns the approximate singular values of the candidate basis.'''
return self.__svals
[docs]
def get_shift_vector(self) -> np.ndarray:
'''Concrete implementation of VectorSpace.get_shift_vector().'''
return self.__shift_vector
[docs]
def get_basis(self) -> np.ndarray:
'''Concrete implementation of VectorSpace.get_basis().'''
return self.__basis
def extents(self) -> Tuple[int, int, int]:
return self.__basis.shape
def _tensor_to_matrix(tensor_input: np.ndarray) -> np.ndarray:
'''
Converts a tensor with shape `[N, M, P]` to a matrix representation
in which the first two dimension are collapsed `[N M, P]`.
'''
output_tensor = tensor_input.reshape(tensor_input.shape[0]*tensor_input.shape[1],
tensor_input.shape[2])
return output_tensor
def _matrix_to_tensor(n_var: int, matrix_input: np.ndarray) -> np.ndarray:
'''
Inverse operation of `_tensor_to_matrix`
'''
d1 = int(matrix_input.shape[0] / n_var)
d2 = matrix_input.shape[1]
output_matrix = matrix_input.reshape(n_var, d1, d2)
return output_matrix