Page MenuHomec4science
Files F60162778

norm_utilities.py
No OneTemporary
Actions

File Metadata

Created
Sat, Apr 27, 23:46

norm_utilities.py
View Options

# Copyright (C) 2018 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.solver.linear_solver import setupSolver
from rrompy.utilities.exception_manager import RROMPyException
__all__ = ['Np2DLike', 'Np2DLikeInv', 'Np2DLikeInvLowRank', 'normEngine']
@abstractmethod
class Np2DLike:
def dot(self, u:Np2D) -> Np2D:
pass
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 not isinstance(u, (np.ndarray,)): u = u.data
if not isinstance(v, (np.ndarray,)): 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)

Event Timeline

c4science · Help