Homec4science
Diffusion RationalROMPy (develop) Browse

RROMPy/rrompy/solver/norm_utilities.py
4fa12654f76edevelop

Last Change
Raw File

norm_utilities.py
View Options

# Copyright (C) 2018-2020 by the RROMPy authors
#
# This file is part of RROMPy.
#
# RROMPy is free software: you can redistribute it and/or modify
# it under the terms of the GNU Lesser General Public License as published by
# the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
#
# RROMPy is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU Lesser General Public License for more details.
#
# You should have received a copy of the GNU Lesser General Public License
# along with RROMPy. If not, see <http://www.gnu.org/licenses/>.
#
from abc import abstractmethod
import numpy as np
from numbers import Number
from copy import deepcopy as copy
from rrompy.utilities.base.types import Np1D, Np2D, DictAny
from rrompy.utilities.numerical import dot as tdot, solve as tsolve
from rrompy.sampling.sample_list import sampleList
from rrompy.parameter.parameter_list import parameterList
from .linear_solver import setupSolver
from rrompy.utilities.exception_manager import RROMPyException
__all__ = ['Np2DLike', 'Np2DLikeGramian', 'Np2DLikeInv', 'Np2DLikeInvLowRank',
'normEngine']
class Np2DLike:
@abstractmethod
def dot(self, u:Np2D) -> Np2D:
pass
class Np2DLikeGramian(Np2DLike):
def __init__(self, L : Np2D = None, R : Np2D = None):
if L is None and R is None:
raise RROMPyException(("Must specify at least one of low-rank "
"factors."))
self.L = R.T.conj() if L is None else L
self.R = L.T.conj() if R is None else R
def dot(self, u:Np2D) -> Np2D:
return tdot(self.L, tdot(self.R, u)).reshape(u.shape)
class Np2DLikeInv(Np2DLike):
def __init__(self, K:Np2D, M:Np2D, solverType:str, solverArgs:DictAny):
self.K, self.M = K, M
self.MH = np.conj(M) if isinstance(self.M, Number) else M.T.conj()
try:
self.solver, self.solverArgs = setupSolver(solverType, solverArgs)
except:
self.solver, self.solverArgs = solverType, solverArgs
def dot(self, u:Np2D) -> Np2D:
return tdot(self.MH, tsolve(self.K, tdot(self.M, u), self.solver,
self.solverArgs)).reshape(u.shape)
@property
def shape(self):
if isinstance(self.M, Number): return self.K.shape
return (self.MH.shape[0], self.M.shape[1])
class Np2DLikeInvLowRank(Np2DLike):
def __init__(self, K:Np2D, M:Np2D, solverType:str, solverArgs:DictAny,
rank:int, oversampling : int = 10, seed : int = 420):
sizeO = K.shape[1] if hasattr(K, "shape") else M.shape[1]
if rank > sizeO:
raise RROMPyException(("Cannot select compressed rank larger than "
"original size."))
if oversampling < 0:
raise RROMPyException("Oversampling parameter must be positive.")
HF = Np2DLikeInv(K, M, solverType, solverArgs)
np.random.seed(seed)
xs = np.random.randn(sizeO, rank + oversampling)
samples = HF.dot(xs)
Q, _ = np.linalg.qr(samples, mode = "reduced")
R = HF.dot(Q).T.conj() # assuming HF (i.e. K) hermitian...
U, s, Vh = np.linalg.svd(R, full_matrices = False)
self.L = Q.dot(U[:, : rank]) * s[: rank]
self.R = Vh[: rank, :]
def dot(self, u:Np2D) -> Np2D:
return tdot(self.L, tdot(self.R, u)).reshape(u.shape)
@property
def shape(self):
return (self.L.shape[0], self.R.shape[1])
class normEngine:
def __init__(self, normMatrix:Np2D):
self.normMatrix = copy(normMatrix)
def innerProduct(self, u:Np2D, v:Np2D, onlyDiag : bool = False) -> Np2D:
if isinstance(u, (parameterList, sampleList)): u = u.data
if isinstance(v, (parameterList, sampleList)): v = v.data
if onlyDiag:
return np.sum(tdot(self.normMatrix, u) * v.conj(), axis = 0)
return tdot(tdot(self.normMatrix, u).T, v.conj()).T
def norm(self, u:Np2D) -> Np1D:
return np.power(np.abs(self.innerProduct(u, u, onlyDiag = True)), .5)
c4science · Help