#
# ************************************************************************
#
# 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)
#
# ************************************************************************
#
'''
Notes
-----
The vector defining the affine offset for a linear subspace is viewed as a matrix of shape
.. math::
\\mathbf{u}_{\\mathrm{shift}} \\in \\mathbb{R}^{N_{\\mathrm{vars}} \\times N_{\\mathrm{x}} }
Theory
------
*What is a shift vector, and why would I use it?* In ROMs, we restrict a state to belong to a low-dimensional affine
vector space,
.. math::
\\mathbf{u} \\approx \\tilde{\\mathbf{u}} \\in \\mathcal{V} + \\mathbf{u}_{\\mathrm{shift}}
where
:math:`\\mathcal{V} \\equiv \\mathrm{range}(\\boldsymbol \\Phi)`. Here :math:`\\mathbf{u}_{\\mathrm{shift}}` defines an
affine offset.
Affine offsets can be useful for a variety of reasons, including satisfying boundary conditions, and satisfying initial
conditions.
The :class:`Shifter` protocol encapsulates the affine offset. The
:class:`StreamingShifter` protocol additionally supports initializing a
data-derived shift vector from snapshot blocks loaded on demand.
API
---
'''
import sys
from numbers import Number
from typing import Protocol
import numpy as np
import romtools.linalg.linalg as la
from romtools.vector_space.utils.snapshot_loader import SnapshotLoader
[docs]
class Shifter(Protocol):
'''Interface for the Shifter class.'''
[docs]
def apply_shift(self, my_array: np.ndarray) -> None:
'''Shifts the snapshot matrix by subtracting a vector generated by the public-facing free functions.'''
...
[docs]
def apply_inverse_shift(self, my_array: np.ndarray) -> None:
'''Shifts the snapshot matrix by adding a vector generated by the public-facing free functions.'''
...
[docs]
def get_shift_vector(self) -> np.ndarray:
'''Returns the vector used to shift the data.'''
...
[docs]
class StreamingShifter(Shifter, Protocol):
'''Shifter interface required by streaming POD vector spaces.'''
[docs]
def initialize_shift_vector_from_loader(
self,
snapshot_loader: SnapshotLoader,
block_size: int,
n_snapshots: int,
comm=None) -> None:
'''Initialize a shift vector from snapshot blocks when required.'''
...
class _Shifter():
'''
Shifts the data by a vector generated by the public-facing free functions.
This class conforms to the :class:`Shifter` protocol.
'''
def __init__(self, shift_vector: np.ndarray) -> None:
'''
Constructor
Args:
shift_vector (np.ndarray): The vector to shift the data by.
'''
self.__shift_vector = shift_vector.copy()
def apply_shift(self, my_array: np.ndarray) -> None:
'''Shifts the input array in place by subtracting the provided shift vector.'''
my_array -= self.__shift_vector[..., None]
def initialize_shift_vector_from_loader(
self,
snapshot_loader: SnapshotLoader,
block_size: int,
n_snapshots: int,
comm=None) -> None:
'''No-op initialization for a shifter with a fixed shift vector.'''
_ = snapshot_loader, block_size, n_snapshots, comm
def apply_inverse_shift(self, my_array: np.ndarray) -> None:
'''Shifts the input array in place by adding the provided shift vector.'''
my_array += self.__shift_vector[..., None]
def get_shift_vector(self) -> np.ndarray:
'''Returns the shift vector.'''
return self.__shift_vector
class _StreamingDataDerivedShifter:
'''Base implementation for shifters initialized from snapshot blocks.'''
def __init__(self) -> None:
self._shift_vector = None
def apply_shift(self, my_array: np.ndarray) -> None:
'''Shift an input block in place.'''
self._check_initialized()
my_array -= self._shift_vector[..., None]
def apply_inverse_shift(self, my_array: np.ndarray) -> None:
'''Inverse-shift an input block in place.'''
self._check_initialized()
my_array += self._shift_vector[..., None]
def get_shift_vector(self) -> np.ndarray:
'''Return the initialized shift vector.'''
self._check_initialized()
return self._shift_vector
def _check_initialized(self) -> None:
if self._shift_vector is None:
raise RuntimeError(
"Streaming shifter has not been initialized from a snapshot loader"
)
class _StreamingAverageShifter(_StreamingDataDerivedShifter):
'''Streaming temporal-average shifter implementation.'''
def initialize_shift_vector_from_loader(
self,
snapshot_loader: SnapshotLoader,
block_size: int,
n_snapshots: int,
comm=None) -> None:
_validate_streaming_shifter_parameters(block_size, n_snapshots)
state_shape = None
snapshot_sum = None
for start in range(0, n_snapshots, block_size):
end = min(start + block_size, n_snapshots)
block = np.asarray(snapshot_loader(start, end))
state_shape = _validate_streaming_snapshot_block(
block, start, end, state_shape, comm
)
if snapshot_sum is None:
accumulation_dtype = np.result_type(block.dtype, np.float64)
snapshot_sum = np.zeros(state_shape, dtype=accumulation_dtype)
snapshot_sum += np.sum(block, axis=-1)
self._shift_vector = snapshot_sum / n_snapshots
class _StreamingFirstVectorShifter(_StreamingDataDerivedShifter):
'''Streaming first-snapshot shifter implementation.'''
def initialize_shift_vector_from_loader(
self,
snapshot_loader: SnapshotLoader,
block_size: int,
n_snapshots: int,
comm=None) -> None:
_validate_streaming_shifter_parameters(block_size, n_snapshots)
block = np.asarray(snapshot_loader(0, 1))
_validate_streaming_snapshot_block(block, 0, 1, None, comm)
self._shift_vector = np.array(block[..., 0], copy=True)
def _validate_streaming_shifter_parameters(
block_size: int, n_snapshots: int) -> None:
if block_size <= 0:
raise ValueError("block_size must be positive")
if n_snapshots <= 0:
raise ValueError("n_snapshots must be positive")
def _validate_streaming_snapshot_block(
block: np.ndarray,
start: int,
end: int,
state_shape,
comm=None):
if block.ndim != 3:
raise ValueError("snapshot loader must return three-dimensional blocks")
if block.shape[-1] != end - start:
raise ValueError("snapshot loader returned an incorrect number of snapshots")
if state_shape is not None and block.shape[:-1] != state_shape:
raise ValueError("snapshot loader returned inconsistent state dimensions")
if state_shape is None and comm is not None and comm.Get_size() > 1:
from mpi4py import MPI
variable_counts = comm.allgather(block.shape[0])
if any(count != variable_counts[0] for count in variable_counts):
raise ValueError("variable count must match across MPI ranks")
minimum_local_dofs = comm.allreduce(block.shape[1], op=MPI.MIN)
if minimum_local_dofs <= 0:
raise ValueError("each MPI rank must own at least one spatial DOF")
return block.shape[:-1]
[docs]
def create_noop_shifter(my_array: np.ndarray) -> StreamingShifter:
'''No op implementation.'''
shift_vector = np.zeros((my_array.shape[0], my_array.shape[1]))
shifter = _Shifter(shift_vector)
return shifter
[docs]
def create_constant_shifter(
shift_value, my_array: np.ndarray) -> StreamingShifter:
'''Shifts the data by a constant value.'''
if isinstance(shift_value, np.ndarray):
shift_vector = np.empty((my_array.shape[0], my_array.shape[1],))
assert my_array.shape[0] == shift_value.size
for i in range(0, my_array.shape[0]):
shift_vector[i] = shift_value[i]
elif isinstance(shift_value, Number):
shift_vector = np.full((my_array.shape[0], my_array.shape[1],), shift_value)
else:
sys.exit("Error: shift_value must be either a number or np.ndarray.")
shifter = _Shifter(shift_vector)
return shifter
[docs]
def create_vector_shifter(shift_vector: np.ndarray) -> StreamingShifter:
'''Shifts the data by a user-input vector.'''
shifter = _Shifter(shift_vector)
return shifter
[docs]
def create_average_shifter(my_array: np.ndarray) -> StreamingShifter:
'''Shifts the data by the average of a data matrix.'''
shift_vector = la.mean(my_array, axis=2)
return _Shifter(shift_vector)
[docs]
def create_firstvec_shifter(my_array: np.ndarray) -> StreamingShifter:
'''Shifts the data by the first vector of a data matrix.'''
shift_vector = my_array[:, :, 0]
shifter = _Shifter(shift_vector)
return shifter
[docs]
def create_streaming_average_shifter() -> StreamingShifter:
'''Create a shifter initialized from the temporal average of streamed data.'''
return _StreamingAverageShifter()
[docs]
def create_streaming_firstvec_shifter() -> StreamingShifter:
'''Create a shifter initialized from the first streamed snapshot.'''
return _StreamingFirstVectorShifter()