1#
  2# ************************************************************************
  3#
  4#                         ROM Tools and Workflows
  5# Copyright 2019 National Technology & Engineering Solutions of Sandia,LLC
  6#                              (NTESS)
  7#
  8# Under the terms of Contract DE-NA0003525 with NTESS, the
  9# U.S. Government retains certain rights in this software.
 10#
 11# ROM Tools and Workflows is licensed under BSD-3-Clause terms of use:
 12#
 13# Redistribution and use in source and binary forms, with or without
 14# modification, are permitted provided that the following conditions
 15# are met:
 16#
 17# 1. Redistributions of source code must retain the above copyright
 18# notice, this list of conditions and the following disclaimer.
 19#
 20# 2. Redistributions in binary form must reproduce the above copyright
 21# notice, this list of conditions and the following disclaimer in the
 22# documentation and/or other materials provided with the distribution.
 23#
 24# 3. Neither the name of the copyright holder nor the names of its
 25# contributors may be used to endorse or promote products derived
 26# from this software without specific prior written permission.
 27#
 28# THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
 29# "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
 30# LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS
 31# FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE
 32# COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT,
 33# INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
 34# (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
 35# SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
 36# HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT,
 37# STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING
 38# IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
 39# POSSIBILITY OF SUCH DAMAGE.
 40#
 41# Questions? Contact Eric Parish (ejparis@sandia.gov)
 42#
 43# ************************************************************************
 44#
 45from typing import Any, Tuple
 46import numpy as np
 47import romtools.linalg.linalg as la
 48
 49
 50class SvdMethodOfSnapshots:
 51    '''
 52    #Parallel implementation of the method of snapshots to mimic the SVD for basis construction
 53    Sample usage:
 54
 55                  mySvd = SvdMethodOfSnapshots(comm)
 56                  U,s,_ = mySvd(snapshots)
 57
 58    where snapshots is the local portion of a distributed memory array.
 59
 60    The standard reduced-basis problem requires solving the optimization problem
 61    $$ \\boldsymbol \\Phi = \\underset{ \\boldsymbol \\Phi_{\\*} \\in \\mathbb{R}^{N \\times K} | \\boldsymbol
 62    \\Phi_{\\*}^T \\boldsymbol \\Phi_{\\*} = \\mathbf{I}}{ \\mathrm{arg \\; min} } \\| \\Phi_{\\*} \\Phi_{\\*}^T
 63    \\mathbf{S} - \\mathbf{S} \\|_2,$$
 64    where $\\mathbf{S} \\in \\mathbb{R}^{N \\times N_s}$, with $N_s$ being the number of snapshots.
 65    The standard way to solve this is with the thin SVD. An alternative approach is to use the method of
 66    snapshts/kernel trick, see, e.g., https://web.stanford.edu/group/frg/course_work/CME345/CA-CME345-Ch4.pdf.
 67    Here, we instead solve the eigenvalue probelm
 68    $$ \\mathbf{S}^T \\mathbf{S} \\boldsymbol \\psi_i = \\lambda_i \\boldsymbol \\psi_i$$
 69    for $i = 1,\\ldots,N_s$. It can be shown that the left singular vectors from the SVD of $\\mathbf{S}$ are
 70    related to the eigen-vectors of the above by
 71    $$ \\mathbf{u}_i = \\frac{1}{\\sqrt{\\lambda_i}} \\mathbf{S} \\boldsymbol \\psi_i.$$
 72
 73    An advantage of the method of snapshots is that it can be easily parallelized and is efficient if we don't
 74    have many snapshots. We compute $\\mathbf{S}^T \\mathbf{S}$ in parallel, and then solve the (typically small)
 75    eigenvalue problem in serial.
 76    '''
 77
 78    def __init__(self, comm) -> None:
 79        self._comm = comm
 80
 81    def __call__(self, snapshots: np.ndarray,
 82                 full_matrices: bool = False,
 83                 compute_uv: bool = False,
 84                 hermitian: bool = False) -> Tuple[np.ndarray, np.ndarray, Any]:
 85        U, s = la.thin_svd(snapshots, self._comm, method='method_of_snapshots')
 86        return U, s, 'not_computed_in_method_of_snapshots'
 87
 88
 89class SvdMethodOfSnapshotsForQr:
 90    '''
 91    Same as SvdMethodOfSnapshots, but call only returns two arguments to be
 92    compatible with QR routine.
 93    '''
 94    def __init__(self, comm) -> None:
 95        self._comm = comm
 96
 97    def __call__(self, snapshots: np.ndarray,
 98                 full_matrices: bool = False,
 99                 compute_uv: bool = False,
100                 hermitian: bool = False) -> Tuple[np.ndarray, Any]:
101        U, _ = la.thin_svd(snapshots, self._comm, method='method_of_snapshots')
102        return U, 'not_computed_in_method_of_snapshots'