diff --git a/rrompy/hfengines/base/marginal_proxy_engine.py b/rrompy/hfengines/base/marginal_proxy_engine.py index de301b4..29ac89b 100644 --- a/rrompy/hfengines/base/marginal_proxy_engine.py +++ b/rrompy/hfengines/base/marginal_proxy_engine.py @@ -1,440 +1,380 @@ # 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 . # import numpy as np import scipy.sparse as scsp from copy import deepcopy as copy, copy as softcopy -from rrompy.utilities.base.types import (Np1D, Np2D, ScOp, Tuple, List, - paramVal, paramList, HFEng, sampList) -from rrompy.utilities.base import multibinom, freepar as fp +from rrompy.utilities.base.types import (Np1D, Np2D, ScOp, List, paramVal, + paramList, HFEng, sampList) +from rrompy.utilities.base import freepar as fp +from rrompy.utilities.numerical import multibinom, marginalizePolyList from rrompy.utilities.poly_fitting.polynomial import ( hashDerivativeToIdx as hashD, hashIdxToDerivative as hashI) from rrompy.utilities.exception_manager import RROMPyAssert from rrompy.parameter import checkParameter, checkParameterList from rrompy.sampling import sampleList, emptySampleList from rrompy.solver.scipy import tensorizeLS, detensorizeLS __all__ = ['MarginalProxyEngine'] class MarginalProxyEngine: """ Marginalized should prescribe fixed value for the marginalized parameters and leave freepar/None elsewhere. """ def __init__(self, HFEngine:HFEng, marginalized:Np1D): self.baseHF = HFEngine self.marg = marginalized self._setupAs() self._setupbs() @property def freeLocations(self): - return tuple([x for x in range(self.nparBase) if self.marg[x] is fp]) + return tuple([x for x in range(self.nparBase) if self.marg[x] == fp]) @property def fixedLocations(self): - return tuple([x for x in range(self.nparBase) \ - if self.marg[x] is not fp]) + return tuple([x for x in range(self.nparBase) if self.marg[x] != fp]) @property def muFixed(self): - muF = checkParameter([self.marg[x] for x in self.fixedLocations]) + muF = checkParameter([m for m in self.marg if m != fp]) if self.nparBase - self.npar > 0: muF = muF[0] return muF @property def rescalingExp(self): return [self.baseHF.rescalingExp[x] for x in self.freeLocations] @property def rescalingExpFixed(self): return [self.baseHF.rescalingExp[x] for x in self.fixedLocations] @property def V(self): return self.baseHF.V def name(self) -> str: return "{}-proxy for {}".format(self.freeLocations, self.baseHF.name()) def __str__(self) -> str: return self.name() def __repr__(self) -> str: return self.__str__() + " at " + hex(id(self)) def __dir_base__(self): return [x for x in self.__dir__() if x[:2] != "__"] def __deepcopy__(self, memo): return softcopy(self) @property def npar(self): """Value of npar.""" return len(self.freeLocations) @property def nparBase(self): """Value of npar.""" return self.baseHF.npar @property def nAs(self): """Value of nAs.""" return self._nAs @property def nbs(self): """Value of nbs.""" return self._nbs @property def nbsH(self) -> int: return max(self.nbs, self.nAs) def spacedim(self): return self.As[0].shape[1] def checkParameter(self, mu:paramList): return checkParameter(mu, self.npar) def checkParameterList(self, mu:paramList): return checkParameterList(mu, self.npar) def buildEnergyNormForm(self): self.baseHF.buildEnergyNormForm() self.energyNormMatrix = self.baseHF.energyNormMatrix def buildEnergyNormDualForm(self): self.baseHF.buildEnergyNormDualForm() self.energyNormDualMatrix = self.baseHF.energyNormDualMatrix def buildDualityPairingForm(self): self.baseHF.buildDualityPairingForm() self.dualityMatrix = self.baseHF.dualityMatrix def buildEnergyNormPartialDualForm(self): self.baseHF.buildEnergyNormPartialDualForm() self.energyNormPartialDualMatrix = ( self.baseHF.energyNormPartialDualMatrix) def innerProduct(self, u:Np2D, v:Np2D, onlyDiag : bool = False, dual : bool = False, duality : bool = True) -> Np2D: return self.baseHF.innerProduct(u, v, onlyDiag, dual, duality) def norm(self, u:Np2D, dual : bool = False, duality : bool = True) -> Np1D: return self.baseHF.norm(u, dual, duality) def checkAInBounds(self, derI : int = 0): """Check if derivative index is oob for operator of linear system.""" if derI < 0 or derI >= self.nAs: d = self.spacedim() return scsp.csr_matrix((np.zeros(0), np.zeros(0), np.zeros(d + 1)), shape = (d, d), dtype = np.complex) def checkbInBounds(self, derI : int = 0, homogeneized : bool = False): """Check if derivative index is oob for RHS of linear system.""" nbs = self.nbsH if homogeneized else self.nbs if derI < 0 or derI >= nbs: return np.zeros(self.spacedim(), dtype = np.complex) def _assembleA(self, mu : paramVal = [], der : List[int] = 0, derI : int = None) -> ScOp: """Assemble (derivative of) operator of linear system.""" mu = self.checkParameter(mu) if self.npar > 0: mu = mu[0] if not hasattr(der, "__len__"): der = [der] * self.npar if derI is None: derI = hashD(der) Anull = self.checkAInBounds(derI) if Anull is not None: return Anull rExp = self.rescalingExp A = copy(self.As[derI]) for j in range(derI + 1, self.nAs): derIdx = hashI(j, self.npar) diffIdx = [x - y for (x, y) in zip(derIdx, der)] if np.all([x >= 0 for x in diffIdx]): A = A + (np.prod((mu ** rExp) ** diffIdx) * multibinom(derIdx, diffIdx) * self.As[j]) return A def A(self, mu : paramVal = [], der : List[int] = 0) -> ScOp: """ Assemble terms of operator of linear system and return it (or its derivative) at a given parameter. """ mu = self.checkParameter(mu) if not hasattr(der, "__len__"): der = [der] * self.npar derI = hashD(der) return self._assembleA(mu, der, derI) def _setupAs(self): - self.As = [] - self._nAs = 0 - for j in range(self.baseHF.nAs): - derjBase = hashI(j, self.nparBase) - jNew = hashD([derjBase[i] for i in self.freeLocations]) - derjFixed = [derjBase[i] for i in self.fixedLocations] - A = (np.prod((self.muFixed ** self.rescalingExpFixed) ** derjFixed) - * self.baseHF.As[j]) - if jNew >= self._nAs: - for _ in range(self._nAs, jNew): - self.As += [self.checkAInBounds(-1)] - self.As += [A] - self._nAs = jNew + 1 - else: - self.As[jNew] = self.As[jNew] + A + mVals = [fp] * self.nparBase + muFixedEff = self.muFixed ** self.rescalingExpFixed + for i, mFE in zip(self.fixedLocations, muFixedEff): + mVals[i] = mFE + self.As = marginalizePolyList(self.baseHF.As, mVals, "auto") + self._nAs = len(self.As) def affineLinearSystemA(self, mu : paramVal = []) -> List[Np2D]: """ Assemble affine blocks of operator of linear system (just linear blocks). """ As = [None] * self.nAs for j in range(self.nAs): As[j] = self.A(mu, hashI(j, self.npar)) return As - def affineWeightsA(self, mu : paramVal = []) -> List[str]: - """ - Assemble affine blocks of operator of linear system (just affine - weights). Stored as strings for the sake of pickling. - """ - mu = self.checkParameter(mu) - lambdasA = ["1."] - mu0Eff = mu ** self.rescalingExp - for j in range(1, self.nAs): - lambdasA += ["np.prod((mu ** ({1}) - [{0}]) ** ({2}))".format( - ','.join([str(x) for x in mu0Eff[0]]), - self.rescalingExp, hashI(j, self.npar))] - return lambdasA - - def affineBlocksA(self, mu : paramVal = [])\ - -> Tuple[List[Np2D], List[str]]: - """Assemble affine blocks of operator of linear system.""" - return self.affineLinearSystemA(mu), self.affineWeightsA(mu) - def _assembleb(self, mu : paramVal = [], der : List[int] = 0, derI : int = None, homogeneized : bool = False) -> ScOp: """Assemble (derivative of) (homogeneized) RHS of linear system.""" mu = self.checkParameter(mu) if self.npar > 0: mu = mu[0] if not hasattr(der, "__len__"): der = [der] * self.npar if derI is None: derI = hashD(der) bnull = self.checkbInBounds(derI, homogeneized) if bnull is not None: return bnull bs = self.bsH if homogeneized else self.bs rExp = self.rescalingExp b = copy(bs[derI]) for j in range(derI + 1, len(bs)): derIdx = hashI(j, self.npar) diffIdx = [x - y for (x, y) in zip(derIdx, der)] if np.all([x >= 0 for x in diffIdx]): b = b + (np.prod((mu ** rExp) ** diffIdx) * multibinom(derIdx, diffIdx) * bs[j]) return b def b(self, mu : paramVal = [], der : List[int] = 0, homogeneized : bool = False) -> Np1D: """ Assemble terms of (homogeneized) RHS of linear system and return it (or its derivative) at a given parameter. """ mu = self.checkParameter(mu) if not hasattr(der, "__len__"): der = [der] * self.npar derI = hashD(der) return self._assembleb(mu, der, derI, homogeneized) def _setupbs(self): - self.bs = [] - self._nbs = 0 - for j in range(self.baseHF.nbs): - derjBase = hashI(j, self.nparBase) - jNew = hashD([derjBase[i] for i in self.freeLocations]) - derjFixed = [derjBase[i] for i in self.fixedLocations] - b = (np.prod((self.muFixed ** self.rescalingExpFixed) ** derjFixed) - * self.baseHF.bs[j]) - if jNew >= self._nbs: - for _ in range(self._nbs, jNew): - self.bs += [self.checkbInBounds(-1)] - self.bs += [b] - self._nbs = jNew + 1 - else: - self.bs[jNew] = self.bs[jNew] + b - - + mVals = [fp] * self.nparBase + muFixedEff = self.muFixed ** self.rescalingExpFixed + for i, mFE in zip(self.fixedLocations, muFixedEff): + mVals[i] = mFE + self.bs = marginalizePolyList(self.baseHF.bs, mVals, "auto") + self._nbs = len(self.bs) + def affineLinearSystemb(self, mu : paramVal = [], homogeneized : bool = False) -> List[Np1D]: """ Assemble affine blocks of RHS of linear system (just linear blocks). """ nbs = self.nbsH if homogeneized else self.nbs bs = [None] * nbs for j in range(nbs): bs[j] = self.b(mu, hashI(j, self.npar), homogeneized) return bs - def affineWeightsb(self, mu : paramVal = [], - homogeneized : bool = False) -> List[str]: - """ - Assemble affine blocks of RHS of linear system (just affine weights). - Stored as strings for the sake of pickling. - """ - mu = self.checkParameter(mu) - nbs = self.nbsH if homogeneized else self.nbs - lambdasb = ["1."] - mu0Eff = mu ** self.rescalingExp - for j in range(1, nbs): - lambdasb += ["np.prod((mu ** ({1}) - [{0}]) ** ({2}))".format( - ','.join([str(x) for x in mu0Eff[0]]), - self.rescalingExp, hashI(j, self.npar))] - return lambdasb - - def affineBlocksb(self, mu : paramVal = [], homogeneized : bool = False)\ - -> Tuple[List[Np1D], List[str]]: - """Assemble affine blocks of RHS of linear system.""" - return (self.affineLinearSystemb(mu, homogeneized), - self.affineWeightsb(mu, homogeneized)) - def solve(self, mu : paramList = [], RHS : sampList = None, homogeneized : bool = False) -> sampList: """ Find solution of linear system. Args: mu: parameter value. RHS: RHS of linear system. If None, defaults to that of parametric system. Defaults to None. """ mu = self.checkParameterList(mu)[0] if mu.shape[0] == 0: mu.reset((1, self.npar), mu.dtype) sol = emptySampleList() if len(mu) > 0: if RHS is None: RHS = [self.b(m, homogeneized = homogeneized) for m in mu] RHS = sampleList(RHS) mult = 0 if len(RHS) == 1 else 1 RROMPyAssert(mult * (len(mu) - 1) + 1, len(RHS), "Sample size") u = self.baseHF._solver(self.A(mu[0]), RHS[0], self.baseHF._solverArgs) sol.reset((len(u), len(mu)), dtype = u.dtype) sol[0] = u for j in range(1, len(mu), self.baseHF._solveBatchSize): if self.baseHF._solveBatchSize != 1: iRange = list(range(j, min(j + self.baseHF._solveBatchSize, len(mu)))) As = [self.A(mu[i]) for i in iRange] bs = [RHS[mult * i] for i in iRange] A, b = tensorizeLS(As, bs) else: A, b = self.A(mu[j]), RHS[mult * j] solStack = self.baseHF._solver(A, b, self.baseHF._solverArgs) if self.baseHF._solveBatchSize != 1: sol[iRange] = detensorizeLS(solStack, len(iRange)) else: sol[j] = solStack return sol def residual(self, u:sampList, mu : paramList = [], homogeneized : bool = False, duality : bool = True) -> sampList: """ Find residual of linear system for given approximate solution. Args: u: numpy complex array with function dofs. If None, set to 0. mu: parameter value. """ mu = self.checkParameterList(mu)[0] if mu.shape[0] == 0: mu.reset((1, self.npar), mu.dtype) if u is not None: u = sampleList(u) mult = 0 if len(u) == 1 else 1 RROMPyAssert(mult * (len(mu) - 1) + 1, len(u), "Sample size") res = emptySampleList() if duality and not hasattr(self, "dualityMatrix"): self.buildDualityPairingForm() for j in range(len(mu)): b = self.b(mu[j], homogeneized = homogeneized) if u is None: r = b else: r = b - self.A(mu[j]).dot(u[mult * j]) if j == 0: res.reset((len(r), len(mu)), dtype = r.dtype) if duality: r = self.dualityMatrix.dot(r) res[j] = r return res def _rayleighQuotient(self, *args, **kwargs) -> float: return self.baseHF._rayleighQuotient(*args, **kwargs) def stabilityFactor(self, u:sampList, mu : paramList = [], nIterP : int = 10, nIterR : int = 10) -> sampList: """ Find stability factor of matrix of linear system using iterative inverse power iteration- and Rayleigh quotient-based procedure. Args: u: numpy complex arrays with function dofs. mu: parameter values. nIterP: number of iterations of power method. nIterR: number of iterations of Rayleigh quotient method. """ mu = self.checkParameterList(mu)[0] if mu.shape[0] == 0: mu.reset((1, self.npar), mu.dtype) u = sampleList(u) mult = 0 if len(u) == 1 else 1 RROMPyAssert(mult * (len(mu) - 1) + 1, len(u), "Sample size") stabFact = np.empty(len(mu), dtype = float) if not hasattr(self, "energyNormMatrix"): self.buildEnergyNormForm() for j in range(len(mu)): stabFact[j] = self._rayleighQuotient(self.A(mu[j]), u[mult * j], self.energyNormMatrix, 0., nIterP, nIterR) return stabFact def plot(self, *args, **kwargs): """ Do some nice plots of the complex-valued function with given dofs. """ self.baseHF.plot(*args, **kwargs) def plotmesh(self, *args, **kwargs): """ Do a nice plot of the mesh. """ self.baseHF.plotmesh(*args, **kwargs) def outParaview(self, *args, **kwargs): """ Output complex-valued function with given dofs to ParaView file. """ self.baseHF.outParaview(*args, **kwargs) def outParaviewTimeDomain(self, *args, **kwargs): """ Output complex-valued function with given dofs to ParaView file, converted to time domain. """ self.baseHF.outParaviewTimeDomain(*args, **kwargs) diff --git a/rrompy/hfengines/base/matrix_engine_base.py b/rrompy/hfengines/base/matrix_engine_base.py index 54fa625..46e0604 100644 --- a/rrompy/hfengines/base/matrix_engine_base.py +++ b/rrompy/hfengines/base/matrix_engine_base.py @@ -1,550 +1,551 @@ # 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 . # from abc import abstractmethod import numpy as np import scipy.sparse as scsp from matplotlib import pyplot as plt from copy import deepcopy as copy, copy as softcopy from rrompy.utilities.base.types import (Np1D, Np2D, ScOp, strLst, Tuple, List, DictAny, paramVal, paramList, sampList) from rrompy.utilities.base import (purgeList, getNewFilename, - verbosityManager as vbMng, multibinom) + verbosityManager as vbMng) +from rrompy.utilities.numerical import multibinom from rrompy.utilities.poly_fitting.polynomial import ( hashDerivativeToIdx as hashD, hashIdxToDerivative as hashI) from rrompy.utilities.exception_manager import RROMPyException, RROMPyAssert from rrompy.parameter import checkParameter, checkParameterList from rrompy.sampling import sampleList, emptySampleList from rrompy.solver import setupSolver, Np2DLikeEye from rrompy.solver.scipy import tensorizeLS, detensorizeLS __all__ = ['MatrixEngineBase'] class MatrixEngineBase: """ Generic solver for parametric matrix problems. Attributes: verbosity: Verbosity level. As: Scipy sparse array representation (in CSC format) of As. bs: Numpy array representation of bs. bsH: Numpy array representation of homogeneized bs. energyNormMatrix: Scipy sparse matrix representing inner product. energyNormDualMatrix: Scipy sparse matrix representing dual inner product. dualityMatrix: Scipy sparse matrix representing duality. energyNormPartialDualMatrix: Scipy sparse matrix representing dual inner product without duality. """ _solveBatchSize = 1 def __init__(self, verbosity : int = 10, timestamp : bool = True): self.verbosity = verbosity self.timestamp = timestamp self.nAs, self.nbs = 1, 1 self.setSolver("SPSOLVE", {"use_umfpack" : False}) self.npar = 0 def name(self) -> str: return self.__class__.__name__ def __str__(self) -> str: return self.name() def __repr__(self) -> str: return self.__str__() + " at " + hex(id(self)) def __dir_base__(self): return [x for x in self.__dir__() if x[:2] != "__"] def __deepcopy__(self, memo): return softcopy(self) @property def npar(self): """Value of npar.""" return self._npar @npar.setter def npar(self, npar): nparOld = self._npar if hasattr(self, "_npar") else -1 if npar != nparOld: self.rescalingExp = [1.] * npar self._npar = npar @property def nAs(self): """Value of nAs.""" return self._nAs @nAs.setter def nAs(self, nAs): self._nAs = nAs self.resetAs() @property def nbs(self): """Value of nbs.""" return self._nbs @nbs.setter def nbs(self, nbs): self._nbs = nbs self.resetbs() @property def nbsH(self) -> int: return max(self.nbs, self.nAs) def spacedim(self): return self.As[0].shape[1] def checkParameter(self, mu:paramList): return checkParameter(mu, self.npar) def checkParameterList(self, mu:paramList): return checkParameterList(mu, self.npar) def buildEnergyNormForm(self): # eye """ Build sparse matrix (in CSR format) representative of scalar product. """ vbMng(self, "INIT", "Assembling energy matrix.", 20) self.energyNormMatrix = Np2DLikeEye() vbMng(self, "DEL", "Done assembling energy matrix.", 20) def buildEnergyNormDualForm(self): """ Build sparse matrix (in CSR format) representative of dual scalar product. """ if not hasattr(self, "energyNormMatrix"): self.buildEnergyNormForm() vbMng(self, "INIT", "Assembling energy dual matrix.", 20) self.energyNormDualMatrix = self.energyNormMatrix vbMng(self, "DEL", "Done assembling energy dual matrix.", 20) def buildDualityPairingForm(self): """Build sparse matrix (in CSR format) representative of duality.""" if not hasattr(self, "energyNormMatrix"): self.buildEnergyNormForm() vbMng(self, "INIT", "Assembling duality matrix.", 20) self.dualityMatrix = self.energyNormMatrix vbMng(self, "DEL", "Done assembling duality matrix.", 20) def buildEnergyNormPartialDualForm(self): """ Build sparse matrix (in CSR format) representative of dual scalar product without duality. """ if not hasattr(self, "energyNormDualMatrix"): self.buildEnergyNormDualForm() vbMng(self, "INIT", "Assembling energy partial dual matrix.", 20) self.energyNormPartialDualMatrix = self.energyNormDualMatrix vbMng(self, "DEL", "Done assembling energy partial dual matrix.", 20) def innerProduct(self, u:Np2D, v:Np2D, onlyDiag : bool = False, dual : bool = False, duality : bool = True) -> Np2D: """Scalar product.""" if dual: if duality: if not hasattr(self, "energyNormDualMatrix"): self.buildEnergyNormDualForm() energyMat = self.energyNormDualMatrix else: if not hasattr(self, "energyNormPartialDualMatrix"): self.buildEnergyNormPartialDualForm() energyMat = self.energyNormPartialDualMatrix else: if not hasattr(self, "energyNormMatrix"): self.buildEnergyNormForm() energyMat = self.energyNormMatrix if not isinstance(u, (np.ndarray,)): u = u.data if not isinstance(v, (np.ndarray,)): v = v.data if onlyDiag: return np.sum(energyMat.dot(u) * v.conj(), axis = 0) return v.T.conj().dot(energyMat.dot(u)) def norm(self, u:Np2D, dual : bool = False, duality : bool = True) -> Np1D: return np.abs(self.innerProduct(u, u, onlyDiag = True, dual = dual, duality = duality)) ** .5 def checkAInBounds(self, derI : int = 0): """Check if derivative index is oob for operator of linear system.""" if derI < 0 or derI >= self.nAs: d = self.spacedim() return scsp.csr_matrix((np.zeros(0), np.zeros(0), np.zeros(d + 1)), shape = (d, d), dtype = np.complex) def checkbInBounds(self, derI : int = 0, homogeneized : bool = False): """Check if derivative index is oob for RHS of linear system.""" nbs = self.nbsH if homogeneized else self.nbs if derI < 0 or derI >= nbs: return np.zeros(self.spacedim(), dtype = np.complex) def resetAs(self): """Reset (derivatives of) operator of linear system.""" self.setAs([None] * self.nAs) if hasattr(self, "_nbs"): self.resetbsH() def resetbs(self): """Reset (derivatives of) RHS of linear system.""" self.setbs([None] * self.nbs) if hasattr(self, "_nAs"): self.resetbsH() def resetbsH(self): """Reset (derivatives of) homogeneized RHS of linear system.""" self.setbsH([None] * self.nbsH) def setAs(self, As:List[Np2D]): """Assign terms of operator of linear system.""" if len(As) != self.nAs: raise RROMPyException(("Expected number {} of terms of As not " "matching given list length {}.").format(self.nAs, len(As))) self.As = [copy(A) for A in As] def setbs(self, bs:List[Np1D]): """Assign terms of RHS of linear system.""" if len(bs) != self.nbs: raise RROMPyException(("Expected number {} of terms of bs not " "matching given list length {}.").format(self.nbs, len(bs))) self.bs = [copy(b) for b in bs] def setbsH(self, bsH:List[Np1D]): """Assign terms of homogeneized RHS of linear system.""" if len(bsH) != self.nbsH: raise RROMPyException(("Expected number {} of terms of bsH not " "matching given list length {}.").format(self.nbsH, len(bsH))) self.bsH = [copy(bH) for bH in bsH] def _assembleA(self, mu : paramVal = [], der : List[int] = 0, derI : int = None) -> ScOp: """Assemble (derivative of) operator of linear system.""" mu = self.checkParameter(mu) if self.npar > 0: mu = mu[0] if not hasattr(der, "__len__"): der = [der] * self.npar if derI is None: derI = hashD(der) Anull = self.checkAInBounds(derI) if Anull is not None: return Anull rExp = self.rescalingExp A = copy(self.As[derI]) for j in range(derI + 1, self.nAs): derIdx = hashI(j, self.npar) diffIdx = [x - y for (x, y) in zip(derIdx, der)] if np.all([x >= 0 for x in diffIdx]): A = A + (np.prod((mu ** rExp) ** diffIdx) * multibinom(derIdx, diffIdx) * self.As[j]) return A @abstractmethod def A(self, mu : paramVal = [], der : List[int] = 0) -> ScOp: """ Assemble terms of operator of linear system and return it (or its derivative) at a given parameter. """ if not hasattr(der, "__len__"): der = [der] * self.npar derI = hashD(der) for j in range(derI, self.nAs): if self.As[j] is None: self.As[j] = self.checkAInBounds(-1) return self._assembleA(mu, der, derI) def affineLinearSystemA(self, mu : paramVal = []) -> List[Np2D]: """ Assemble affine blocks of operator of linear system (just linear blocks). """ As = [None] * self.nAs for j in range(self.nAs): As[j] = self.A(mu, hashI(j, self.npar)) return As - def affineWeightsA(self, mu : paramVal = []) -> List[str]: - """ - Assemble affine blocks of operator of linear system (just affine - weights). Stored as strings for the sake of pickling. - """ - mu = self.checkParameter(mu) - lambdasA = ["1."] - mu0Eff = mu ** self.rescalingExp - for j in range(1, self.nAs): - lambdasA += ["np.prod((mu ** ({1}) - [{0}]) ** ({2}))".format( - ','.join([str(x) for x in mu0Eff[0]]), - self.rescalingExp, hashI(j, self.npar))] - return lambdasA - - def affineBlocksA(self, mu : paramVal = [])\ - -> Tuple[List[Np2D], List[str]]: - """Assemble affine blocks of operator of linear system.""" - return self.affineLinearSystemA(mu), self.affineWeightsA(mu) +# def affineWeightsA(self, mu : paramVal = []) -> List[str]: +# """ +# Assemble affine blocks of operator of linear system (just affine +# weights). Stored as strings for the sake of pickling. +# """ +# mu = self.checkParameter(mu) +# lambdasA = ["1."] +# mu0Eff = mu ** self.rescalingExp +# for j in range(1, self.nAs): +# lambdasA += ["np.prod((mu ** ({1}) - [{0}]) ** ({2}))".format( +# ','.join([str(x) for x in mu0Eff[0]]), +# self.rescalingExp, hashI(j, self.npar))] +# return lambdasA + +# def affineBlocksA(self, mu : paramVal = [])\ +# -> Tuple[List[Np2D], List[str]]: +# """Assemble affine blocks of operator of linear system.""" +# return self.affineLinearSystemA(mu), self.affineWeightsA(mu) def _assembleb(self, mu : paramVal = [], der : List[int] = 0, derI : int = None, homogeneized : bool = False) -> ScOp: """Assemble (derivative of) (homogeneized) RHS of linear system.""" mu = self.checkParameter(mu) if self.npar > 0: mu = mu[0] if not hasattr(der, "__len__"): der = [der] * self.npar if derI is None: derI = hashD(der) bnull = self.checkbInBounds(derI, homogeneized) if bnull is not None: return bnull bs = self.bsH if homogeneized else self.bs rExp = self.rescalingExp b = copy(bs[derI]) for j in range(derI + 1, len(bs)): derIdx = hashI(j, self.npar) diffIdx = [x - y for (x, y) in zip(derIdx, der)] if np.all([x >= 0 for x in diffIdx]): b = b + (np.prod((mu ** rExp) ** diffIdx) * multibinom(derIdx, diffIdx) * bs[j]) return b @abstractmethod def b(self, mu : paramVal = [], der : List[int] = 0, homogeneized : bool = False) -> Np1D: """ Assemble terms of (homogeneized) RHS of linear system and return it (or its derivative) at a given parameter. """ if not hasattr(der, "__len__"): der = [der] * self.npar derI = hashD(der) if homogeneized: for j in range(derI, self.nbsH): if self.bsH[j] is None: self.bsH[j] = self.checkbInBounds(-1) else: for j in range(derI, self.nbs): if self.bs[j] is None: self.bs[j] = self.checkbInBounds(-1) return self._assembleb(mu, der, derI, homogeneized) def affineLinearSystemb(self, mu : paramVal = [], homogeneized : bool = False) -> List[Np1D]: """ Assemble affine blocks of RHS of linear system (just linear blocks). """ nbs = self.nbsH if homogeneized else self.nbs bs = [None] * nbs for j in range(nbs): bs[j] = self.b(mu, hashI(j, self.npar), homogeneized) return bs - def affineWeightsb(self, mu : paramVal = [], - homogeneized : bool = False) -> List[str]: - """ - Assemble affine blocks of RHS of linear system (just affine weights). - Stored as strings for the sake of pickling. - """ - mu = self.checkParameter(mu) - nbs = self.nbsH if homogeneized else self.nbs - lambdasb = ["1."] - mu0Eff = mu ** self.rescalingExp - for j in range(1, nbs): - lambdasb += ["np.prod((mu ** ({1}) - [{0}]) ** ({2}))".format( - ','.join([str(x) for x in mu0Eff[0]]), - self.rescalingExp, hashI(j, self.npar))] - return lambdasb - - def affineBlocksb(self, mu : paramVal = [], homogeneized : bool = False)\ - -> Tuple[List[Np1D], List[str]]: - """Assemble affine blocks of RHS of linear system.""" - return (self.affineLinearSystemb(mu, homogeneized), - self.affineWeightsb(mu, homogeneized)) +# def affineWeightsb(self, mu : paramVal = [], +# homogeneized : bool = False) -> List[str]: +# """ +# Assemble affine blocks of RHS of linear system (just affine weights). +# Stored as strings for the sake of pickling. +# """ +# mu = self.checkParameter(mu) +# nbs = self.nbsH if homogeneized else self.nbs +# lambdasb = ["1."] +# mu0Eff = mu ** self.rescalingExp +# for j in range(1, nbs): +# lambdasb += ["np.prod((mu ** ({1}) - [{0}]) ** ({2}))".format( +# ','.join([str(x) for x in mu0Eff[0]]), +# self.rescalingExp, hashI(j, self.npar))] +# return lambdasb + +# def affineBlocksb(self, mu : paramVal = [], homogeneized : bool = False)\ +# -> Tuple[List[Np1D], List[str]]: +# """Assemble affine blocks of RHS of linear system.""" +# return (self.affineLinearSystemb(mu, homogeneized), +# self.affineWeightsb(mu, homogeneized)) def setSolver(self, solverType:str, solverArgs : DictAny = {}): """Choose solver type and parameters.""" self._solver, self._solverArgs = setupSolver(solverType, solverArgs) def solve(self, mu : paramList = [], RHS : sampList = None, homogeneized : bool = False) -> sampList: """ Find solution of linear system. Args: mu: parameter value. RHS: RHS of linear system. If None, defaults to that of parametric system. Defaults to None. """ mu = self.checkParameterList(mu)[0] if mu.shape[0] == 0: mu.reset((1, self.npar), mu.dtype) sol = emptySampleList() if len(mu) > 0: if RHS is None: RHS = [self.b(m, homogeneized = homogeneized) for m in mu] RHS = sampleList(RHS) mult = 0 if len(RHS) == 1 else 1 RROMPyAssert(mult * (len(mu) - 1) + 1, len(RHS), "Sample size") u = self._solver(self.A(mu[0]), RHS[0], self._solverArgs) sol.reset((len(u), len(mu)), dtype = u.dtype) sol[0] = u for j in range(1, len(mu), self._solveBatchSize): if self._solveBatchSize != 1: iRange = list(range(j, min(j + self._solveBatchSize, len(mu)))) As = [self.A(mu[i]) for i in iRange] bs = [RHS[mult * i] for i in iRange] A, b = tensorizeLS(As, bs) else: A, b = self.A(mu[j]), RHS[mult * j] solStack = self._solver(A, b, self._solverArgs) if self._solveBatchSize != 1: sol[iRange] = detensorizeLS(solStack, len(iRange)) else: sol[j] = solStack return sol def residual(self, u:sampList, mu : paramList = [], homogeneized : bool = False, duality : bool = True) -> sampList: """ Find residual of linear system for given approximate solution. Args: u: numpy complex array with function dofs. If None, set to 0. mu: parameter value. """ mu = self.checkParameterList(mu)[0] if mu.shape[0] == 0: mu.reset((1, self.npar), mu.dtype) if u is not None: u = sampleList(u) mult = 0 if len(u) == 1 else 1 RROMPyAssert(mult * (len(mu) - 1) + 1, len(u), "Sample size") res = emptySampleList() if duality and not hasattr(self, "dualityMatrix"): self.buildDualityPairingForm() for j in range(len(mu)): b = self.b(mu[j], homogeneized = homogeneized) if u is None: r = b else: r = b - self.A(mu[j]).dot(u[mult * j]) if j == 0: res.reset((len(r), len(mu)), dtype = r.dtype) if duality: r = self.dualityMatrix.dot(r) res[j] = r return res def _rayleighQuotient(self, A:Np2D, v0:Np1D, M:Np2D, sigma : float = 0., nIterP : int = 10, nIterR : int = 10) -> float: nIterP = min(nIterP, len(v0) // 2) nIterR = min(nIterR, (len(v0) + 1) // 2) v0 /= v0.T.conj().dot(M.dot(v0)) ** .5 for j in range(nIterP): v0 = self._solver(A - sigma * M, M.dot(v0), self._solverArgs) v0 /= v0.T.conj().dot(M.dot(v0)) ** .5 l0 = v0.T.conj().dot(A.dot(v0)) for j in range(nIterR): v0 = self._solver(A - l0 * M, M.dot(v0), self._solverArgs) v0 /= v0.T.conj().dot(M.dot(v0)) ** .5 l0 = v0.T.conj().dot(A.dot(v0)) if np.isnan(l0): l0 = np.finfo(float).eps return np.abs(l0) def stabilityFactor(self, u:sampList, mu : paramList = [], nIterP : int = 10, nIterR : int = 10) -> sampList: """ Find stability factor of matrix of linear system using iterative inverse power iteration- and Rayleigh quotient-based procedure. Args: u: numpy complex arrays with function dofs. mu: parameter values. nIterP: number of iterations of power method. nIterR: number of iterations of Rayleigh quotient method. """ mu = self.checkParameterList(mu)[0] if mu.shape[0] == 0: mu.reset((1, self.npar), mu.dtype) u = sampleList(u) mult = 0 if len(u) == 1 else 1 RROMPyAssert(mult * (len(mu) - 1) + 1, len(u), "Sample size") stabFact = np.empty(len(mu), dtype = float) if not hasattr(self, "energyNormMatrix"): self.buildEnergyNormForm() for j in range(len(mu)): stabFact[j] = self._rayleighQuotient(self.A(mu[j]), u[mult * j], self.energyNormMatrix, 0., nIterP, nIterR) return stabFact def plot(self, u:Np1D, warping : List[callable] = None, name : str = "u", save : str = None, what : strLst = 'all', saveFormat : str = "eps", saveDPI : int = 100, show : bool = True, pyplotArgs : dict = {}, **figspecs): """ Do some nice plots of the complex-valued function with given dofs. Args: u: numpy complex array with function dofs. name(optional): Name to be shown as title of the plots. Defaults to 'u'. save(optional): Where to save plot(s). Defaults to None, i.e. no saving. what(optional): Which plots to do. If list, can contain 'ABS', 'PHASE', 'REAL', 'IMAG'. If str, same plus wildcard 'ALL'. Defaults to 'ALL'. saveFormat(optional): Format for saved plot(s). Defaults to "eps". saveDPI(optional): DPI for saved plot(s). Defaults to 100. show(optional): Whether to show figure. Defaults to True. pyplotArgs(optional): Optional arguments for pyplot. figspecs(optional key args): Optional arguments for matplotlib figure creation. """ if isinstance(what, (str,)): if what.upper() == 'ALL': what = ['ABS', 'PHASE', 'REAL', 'IMAG'] else: what = [what] what = purgeList(what, ['ABS', 'PHASE', 'REAL', 'IMAG'], listname = self.name() + ".what", baselevel = 1) if len(what) == 0: return if 'figsize' not in figspecs.keys(): figspecs['figsize'] = (13. * len(what) / 4, 3) subplotcode = 100 + len(what) * 10 idxs = np.arange(self.spacedim()) if warping is not None: idxs = warping[0](np.arange(self.spacedim())) plt.figure(**figspecs) plt.jet() if 'ABS' in what: subplotcode = subplotcode + 1 plt.subplot(subplotcode) plt.plot(idxs, np.abs(u).flatten(), **pyplotArgs) plt.title("|{0}|".format(name)) if 'PHASE' in what: subplotcode = subplotcode + 1 plt.subplot(subplotcode) plt.plot(idxs, np.angle(u).flatten(), **pyplotArgs) plt.title("phase({0})".format(name)) if 'REAL' in what: subplotcode = subplotcode + 1 plt.subplot(subplotcode) plt.plot(idxs, np.real(u).flatten(), **pyplotArgs) plt.title("Re({0})".format(name)) if 'IMAG' in what: subplotcode = subplotcode + 1 plt.subplot(subplotcode) plt.plot(idxs, np.imag(u).flatten(), **pyplotArgs) plt.title("Im({0})".format(name)) if save is not None: save = save.strip() plt.savefig(getNewFilename("{}_fig_".format(save), saveFormat), format = saveFormat, dpi = saveDPI) if show: plt.show() plt.close() diff --git a/rrompy/hfengines/linear_problem/bidimensional/helmholtz_square_domain_problem_engine.py b/rrompy/hfengines/linear_problem/bidimensional/helmholtz_square_domain_problem_engine.py index aa3c9d8..343e450 100644 --- a/rrompy/hfengines/linear_problem/bidimensional/helmholtz_square_domain_problem_engine.py +++ b/rrompy/hfengines/linear_problem/bidimensional/helmholtz_square_domain_problem_engine.py @@ -1,162 +1,159 @@ # 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 . # import numpy as np from numpy.polynomial.polynomial import polyfit as fit import fenics as fen from rrompy.utilities.base.types import (ScOp, List, Tuple, paramVal, Np1D, FenExpr) from rrompy.solver.fenics import fenZERO from rrompy.hfengines.linear_problem.helmholtz_problem_engine import ( HelmholtzProblemEngine) from rrompy.utilities.base import verbosityManager as vbMng from rrompy.utilities.poly_fitting.polynomial import ( hashDerivativeToIdx as hashD, hashIdxToDerivative as hashI) from rrompy.solver.fenics import fenics2Sparse, fenics2Vector __all__ = ['HelmholtzSquareDomainProblemEngine'] class HelmholtzSquareDomainProblemEngine(HelmholtzProblemEngine): """ Solver for square Helmholtz problems with parametric laplacian. - \dxx u - mu_2**2 \dyy u - mu_1**2 * u = f(mu_2) in \Omega = [0,\pi]^2 u = 0 on \partial\Omega """ def __init__(self, kappa:float, theta:float, n:int, mu0 : paramVal = [12. ** .5, 1.], degree_threshold : int = np.inf, verbosity : int = 10, timestamp : bool = True): super().__init__(mu0 = mu0, degree_threshold = degree_threshold, verbosity = verbosity, timestamp = timestamp) self.nAs, self.nbs = 6, 11 * 12 // 2 self.npar = 2 self.rescalingExp = [2., 1.] self._theta = theta self._kappa = kappa pi = np.pi mesh = fen.RectangleMesh(fen.Point(0, 0), fen.Point(pi, pi), 3 * n, 3 * n) self.V = fen.FunctionSpace(mesh, "P", 1) def getForcingTerm(self, mu : paramVal = []) -> Tuple[FenExpr, FenExpr]: """Compute forcing term.""" - mu = self.checkParameter(mu) vbMng(self, "INIT", ("Assembling base expression for forcing term " "at {}.").format(mu), 25) - mu = mu(0, 1) c, s = np.cos(self._theta), np.sin(self._theta) x, y = fen.SpatialCoordinate(self.V.mesh())[:] forcingTerm = [fen.cos(self._kappa * (c * x + s / mu * y)), fen.sin(self._kappa * (c * x + s / mu * y))] vbMng(self, "DEL", "Done assembling base expression.", 25) return forcingTerm def A(self, mu : paramVal = [], der : List[int] = 0) -> ScOp: """Assemble (derivative of) operator of linear system.""" mu = self.checkParameter(mu) if not hasattr(der, "__len__"): der = [der] * self.npar derI = hashD(der) self.autoSetDS() for j in range(2, 5): if derI <= j and self.As[j] is None: self.As[j] = self.checkAInBounds(-1) if derI <= 0 and self.As[0] is None: vbMng(self, "INIT", "Assembling operator term A0.", 20) DirichletBC0 = fen.DirichletBC(self.V, fenZERO, self.DirichletBoundary) a0Re = fen.dot(self.u.dx(0), self.v.dx(0)) * fen.dx self.As[0] = fenics2Sparse(a0Re, {}, DirichletBC0, 1) vbMng(self, "DEL", "Done assembling operator term.", 20) if derI <= 1 and self.As[1] is None: vbMng(self, "INIT", "Assembling operator term A1.", 20) DirichletBC0 = fen.DirichletBC(self.V, fenZERO, self.DirichletBoundary) nRe, nIm = self.refractionIndex n2Re, n2Im = nRe * nRe - nIm * nIm, 2 * nRe * nIm parsRe = self.iterReduceQuadratureDegree(zip([n2Re], ["refractionIndexSquaredReal"])) parsIm = self.iterReduceQuadratureDegree(zip([n2Im], ["refractionIndexSquaredImag"])) a1Re = - n2Re * fen.dot(self.u, self.v) * fen.dx a1Im = - n2Im * fen.dot(self.u, self.v) * fen.dx self.As[1] = (fenics2Sparse(a1Re, parsRe, DirichletBC0, 0) + 1.j * fenics2Sparse(a1Im, parsIm, DirichletBC0, 0)) vbMng(self, "DEL", "Done assembling operator term.", 20) if derI <= 5 and self.As[5] is None: vbMng(self, "INIT", "Assembling operator term A5.", 20) DirichletBC0 = fen.DirichletBC(self.V, fenZERO, self.DirichletBoundary) a5Re = fen.dot(self.u.dx(1), self.v.dx(1)) * fen.dx self.As[5] = fenics2Sparse(a5Re, {}, DirichletBC0, 0) vbMng(self, "DEL", "Done assembling operator term.", 20) return self._assembleA(mu, der, derI) def b(self, mu : paramVal = [], der : List[int] = 0, homogeneized : bool = False) -> Np1D: """Assemble (derivative of) RHS of linear system.""" mu = self.checkParameter(mu) if not hasattr(der, "__len__"): der = [der] * self.npar derI = hashD(der) nbsTot = self.nbsH if homogeneized else self.nbs bs = self.bsH if homogeneized else self.bs if homogeneized and self.mu0 != self.mu0BC: self.liftDirichletData(self.mu0) bDEIMCoeffs = None for j in range(derI, nbsTot): derH = hashI(j, self.npar) if bs[j] is None: if np.sum(derH) != derH[-1]: if homogeneized: self.bsH[j] = self.checkbInBounds(-1) Ader = self.A(0, derH) self.bsH[j] -= Ader.dot(self.liftedDirichletDatum) else: self.bs[j] = self.checkbInBounds(-1) continue vbMng(self, "INIT", "Assembling forcing term b{}.".format(j), 20) if bDEIMCoeffs is None: - bDEIM = np.empty((self.nbs, self.spacedim()), - dtype = np.complex) - muDEIM = np.linspace(.5, 4., - np.sum(hashI(nbsTot - 1, self.npar))) + bDEIM = np.empty((np.sum(hashI(nbsTot, self.npar)), + self.spacedim()), dtype = np.complex) + muDEIM = np.linspace(.5, 4., bDEIM.shape[0]) for jj, muD in enumerate(muDEIM): fRe, fIm = self.getForcingTerm(muD) parsRe = self.iterReduceQuadratureDegree(zip([fRe], ["forcingTerm{}Real".format(jj)])) parsIm = self.iterReduceQuadratureDegree(zip([fIm], ["forcingTerm{}Imag".format(jj)])) LR = fen.dot(fRe, self.v) * fen.dx LI = fen.dot(fIm, self.v) * fen.dx DBC0 = fen.DirichletBC(self.V, fenZERO, self.DirichletBoundary) bDEIM[jj] = (fenics2Vector(LR, parsRe, DBC0, 1) + 1.j * fenics2Vector(LI, parsIm, DBC0, 1)) bDEIMCoeffs = fit(muDEIM, bDEIM, len(muDEIM) - 1) b = bDEIMCoeffs[derH[-1]] if homogeneized: Ader = self.A(0, derH) b -= Ader.dot(self.liftedDirichletDatum) if homogeneized: self.bsH[j] = b else: self.bs[j] = b vbMng(self, "DEL", "Done assembling forcing term.", 20) return self._assembleb(mu, der, derI, homogeneized) diff --git a/rrompy/hfengines/linear_problem/helmholtz_square_bubble_domain_problem_engine.py b/rrompy/hfengines/linear_problem/helmholtz_square_bubble_domain_problem_engine.py index a50171c..a455d2b 100644 --- a/rrompy/hfengines/linear_problem/helmholtz_square_bubble_domain_problem_engine.py +++ b/rrompy/hfengines/linear_problem/helmholtz_square_bubble_domain_problem_engine.py @@ -1,162 +1,160 @@ # 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 . # import numpy as np from numpy.polynomial.polynomial import polyfit as fit import fenics as fen from rrompy.utilities.base.types import (Np1D, ScOp, Tuple, List, FenExpr, paramVal) from rrompy.solver.fenics import fenZERO from .helmholtz_problem_engine import HelmholtzProblemEngine from rrompy.utilities.base import verbosityManager as vbMng from rrompy.utilities.poly_fitting.polynomial import ( hashDerivativeToIdx as hashD, hashIdxToDerivative as hashI) from rrompy.solver.fenics import fenics2Sparse, fenics2Vector __all__ = ['HelmholtzSquareBubbleDomainProblemEngine'] class HelmholtzSquareBubbleDomainProblemEngine(HelmholtzProblemEngine): """ Solver for square bubble Helmholtz problems with parametric domain heigth. - \Delta u - kappa^2 * u = f in \Omega_mu = [0,\pi] x [0,\mu\pi] u = 0 on \Gamma_mu = \partial\Omega_mu with exact solution square bubble times plane wave. """ def __init__(self, kappa:float, theta:float, n:int, mu0 : paramVal = [1.], degree_threshold : int = np.inf, verbosity : int = 10, timestamp : bool = True): super().__init__(mu0 = mu0, degree_threshold = degree_threshold, verbosity = verbosity, timestamp = timestamp) self.nAs, self.nbs = 3, 15 self.kappa = kappa self.theta = theta self.forcingTermMu = np.nan mesh = fen.RectangleMesh(fen.Point(0, 0), fen.Point(np.pi, np.pi), 3 * n, 3 * n) self.V = fen.FunctionSpace(mesh, "P", 1) self.rescalingExp = [1.] def getForcingTerm(self, mu : paramVal = []) -> Tuple[FenExpr, FenExpr]: """Compute forcing term.""" - mu = self.checkParameter(mu) vbMng(self, "INIT", ("Assembling base expression for forcing term " "at {}.").format(mu), 25) pi = np.pi c, s = np.cos(self.theta), np.sin(self.theta) x, y = fen.SpatialCoordinate(self.V.mesh())[:] - muR, muI = np.real(mu(0, 0)), np.imag(mu(0, 0)) - mu2R, mu2I = np.real(mu(0, 0) ** 2.), np.imag(mu(0, 0) ** 2.) + muR, muI = np.real(mu), np.imag(mu) + mu2R, mu2I = np.real(mu ** 2.), np.imag(mu ** 2.) C = 16. / pi ** 4. bR = C * (2 * (x * (pi - x) + y * (pi - y)) + (self.kappa * s) ** 2. * (mu2R - 1.) * x * (pi - x) * y * (pi - y)) bI = C * (2 * self.kappa * (c * (pi - 2 * x) * y * (pi - y) + s * x * (pi - x) * (pi - 2 * y)) + (self.kappa * s) ** 2. * mu2I * x * (pi - x) * y * (pi - y)) wR = (fen.cos(self.kappa * (c * x + s * muR * y)) * fen.exp(self.kappa * s * muI * y)) wI = (fen.sin(self.kappa * (c * x + s * muR * y)) * fen.exp(self.kappa * s * muI * y)) fRe, fIm = bR * wR + bI * wI, bI * wR - bR * wI - cRe, cIm = np.real(mu(0, 0) ** 2.), np.imag(mu(0, 0) ** 2.) - forcingTerm = [cRe * fRe - cIm * fIm + fenZERO, - cRe * fIm + cIm * fRe + fenZERO] + forcingTerm = [mu2R * fRe - mu2I * fIm + fenZERO, + mu2R * fIm + mu2I * fRe + fenZERO] vbMng(self, "DEL", "Done assembling base expression.", 25) return forcingTerm def A(self, mu : paramVal = [], der : List[int] = 0) -> ScOp: """Assemble (derivative of) operator of linear system.""" mu = self.checkParameter(mu) if not hasattr(der, "__len__"): der = [der] * self.npar derI = hashD(der) self.autoSetDS() if derI <= 0 and self.As[0] is None: vbMng(self, "INIT", "Assembling operator term A0.", 20) DirichletBC0 = fen.DirichletBC(self.V, fenZERO, self.DirichletBoundary) a0Re = fen.dot(self.u.dx(1), self.v.dx(1)) * fen.dx self.As[0] = fenics2Sparse(a0Re, {}, DirichletBC0, 1) vbMng(self, "DEL", "Done assembling operator term.", 20) if derI <= 1 and self.As[1] is None: self.As[1] = self.checkAInBounds(-1) if derI <= 2 and self.As[2] is None: vbMng(self, "INIT", "Assembling operator term A2.", 20) DirichletBC0 = fen.DirichletBC(self.V, fenZERO, self.DirichletBoundary) nRe, nIm = self.refractionIndex n2Re, n2Im = nRe * nRe - nIm * nIm, 2 * nRe * nIm k2Re, k2Im = np.real(self.omega ** 2), np.imag(self.omega ** 2) k2n2Re = k2Re * n2Re - k2Im * n2Im k2n2Im = k2Re * n2Im + k2Im * n2Re parsRe = self.iterReduceQuadratureDegree(zip([k2n2Re], ["kappaSquaredRefractionIndexSquaredReal"])) parsIm = self.iterReduceQuadratureDegree(zip([k2n2Im], ["kappaSquaredRefractionIndexSquaredImag"])) a2Re = (fen.dot(self.u.dx(0), self.v.dx(0)) - k2n2Re * fen.dot(self.u, self.v)) * fen.dx a2Im = - k2n2Im * fen.dot(self.u, self.v) * fen.dx self.As[2] = (fenics2Sparse(a2Re, parsRe, DirichletBC0, 0) + 1.j * fenics2Sparse(a2Im, parsIm, DirichletBC0, 0)) vbMng(self, "DEL", "Done assembling operator term.", 20) return self._assembleA(mu, der, derI) def b(self, mu : paramVal = [], der : List[int] = 0, homogeneized : bool = False) -> Np1D: """Assemble (derivative of) RHS of linear system.""" mu = self.checkParameter(mu) if not hasattr(der, "__len__"): der = [der] * self.npar derI = hashD(der) nbsTot = self.nbsH if homogeneized else self.nbs bs = self.bsH if homogeneized else self.bs if homogeneized and self.mu0 != self.mu0BC: self.liftDirichletData(self.mu0) bDEIMCoeffs = None for j in range(derI, nbsTot): if bs[j] is None: vbMng(self, "INIT", "Assembling forcing term b{}.".format(j), 20) if bDEIMCoeffs is None: bDEIM = np.empty((self.nbs, self.spacedim()), dtype = np.complex) muDEIM = np.linspace(.5, 4., self.nbs) for jj, muD in enumerate(muDEIM): fRe, fIm = self.getForcingTerm(muD) parsRe = self.iterReduceQuadratureDegree(zip([fRe], ["forcingTerm{}Real".format(jj)])) parsIm = self.iterReduceQuadratureDegree(zip([fIm], ["forcingTerm{}Imag".format(jj)])) LR = fen.dot(fRe, self.v) * fen.dx LI = fen.dot(fIm, self.v) * fen.dx DBC0 = fen.DirichletBC(self.V, fenZERO, self.DirichletBoundary) bDEIM[jj] = (fenics2Vector(LR, parsRe, DBC0, 1) + 1.j * fenics2Vector(LI, parsIm, DBC0, 1)) bDEIMCoeffs = fit(muDEIM, bDEIM, self.nbs - 1) b = bDEIMCoeffs[j] if homogeneized: Ader = self.A(0, hashI(j, self.npar)) b -= Ader.dot(self.liftedDirichletDatum) if homogeneized: self.bsH[j] = b else: self.bs[j] = b vbMng(self, "DEL", "Done assembling forcing term.", 20) return self._assembleb(mu, der, derI, homogeneized) diff --git a/rrompy/parameter/parameter_sampling/fft_sampler.py b/rrompy/parameter/parameter_sampling/fft_sampler.py index d55a63f..5ea7e20 100644 --- a/rrompy/parameter/parameter_sampling/fft_sampler.py +++ b/rrompy/parameter/parameter_sampling/fft_sampler.py @@ -1,50 +1,50 @@ # 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 . # import numpy as np from .generic_sampler import GenericSampler from rrompy.utilities.base.types import List, paramList -from rrompy.utilities.base import lowDiscrepancy, kroneckerer +from rrompy.utilities.numerical import lowDiscrepancy, kroneckerer from rrompy.parameter import checkParameterList __all__ = ['FFTSampler'] class FFTSampler(GenericSampler): """Generator of FFT-type sample points on scaled roots of unity.""" def generatePoints(self, n:List[int], reorder : bool = True) -> paramList: """Array of sample points.""" if not hasattr(n, "__len__"): n = [n] super().generatePoints(n) nleft, nright = 1, np.prod(n) xmat = np.empty((nright, self.npar), dtype = np.complex) for d in range(self.npar): nright //= n[d] a, b = self.lims(0, d), self.lims(1, d) if self.scaling is not None: a, b = self.scaling[d](a), self.scaling[d](b) c, r = (a + b) / 2., np.abs(a - b) / 2. xd = c + r * np.exp(1.j * np.linspace(0, 2 * np.pi, n[d] + 1)[:-1]) if self.scalingInv is not None: xd = self.scalingInv[d](xd) xmat[:, d] = kroneckerer(xd, nleft, nright) nleft *= n[d] if reorder: xmat = xmat[lowDiscrepancy(np.prod(n)), :] x = checkParameterList(xmat, self.npar)[0] return x diff --git a/rrompy/parameter/parameter_sampling/quadrature_sampler.py b/rrompy/parameter/parameter_sampling/quadrature_sampler.py index 1ca7bda..5fb7ca9 100644 --- a/rrompy/parameter/parameter_sampling/quadrature_sampler.py +++ b/rrompy/parameter/parameter_sampling/quadrature_sampler.py @@ -1,87 +1,87 @@ # 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 . # import numpy as np from .generic_sampler import GenericSampler from rrompy.utilities.base.types import List, paramList from rrompy.utilities.exception_manager import RROMPyException -from rrompy.utilities.base import lowDiscrepancy, kroneckerer +from rrompy.utilities.numerical import lowDiscrepancy, kroneckerer from rrompy.parameter import checkParameterList __all__ = ['QuadratureSampler'] class QuadratureSampler(GenericSampler): """Generator of quadrature sample points.""" allowedKinds = ["UNIFORM", "CHEBYSHEV", "GAUSSLEGENDRE"] def __init__(self, lims:paramList, kind : str = "UNIFORM", scaling : List[callable] = None, scalingInv : List[callable] = None): super().__init__(lims = lims, scaling = scaling, scalingInv = scalingInv) self.kind = kind def __str__(self) -> str: return "{}_{}".format(super().__str__(), self.kind) def __repr__(self) -> str: return self.__str__() + " at " + hex(id(self)) @property def kind(self): """Value of kind.""" return self._kind @kind.setter def kind(self, kind): if kind.upper() not in self.allowedKinds: raise RROMPyException("Generator kind not recognized.") self._kind = kind.upper() def generatePoints(self, n:List[int], reorder : bool = True) -> paramList: """Array of sample points.""" if not hasattr(n, "__len__"): n = [n] super().generatePoints(n) nleft, nright = 1, np.prod(n) xmat = np.empty((nright, self.npar), dtype = self.lims.dtype) for d in range(self.npar): nright //= n[d] a, b = self.lims(0, d), self.lims(1, d) if self.scaling is not None: a, b = self.scaling[d](a), self.scaling[d](b) c, r = (a + b) / 2., (a - b) / 2. dAbs = 2. * np.abs(r) if self.kind == "UNIFORM": xd = np.linspace(a, b, n[d]) elif self.kind == "CHEBYSHEV": nodes, _ = np.polynomial.chebyshev.chebgauss(n[d]) xd = c + r * nodes elif self.kind == "GAUSSLEGENDRE": nodes, _ = np.polynomial.legendre.leggauss(n[d]) xd = c + r * nodes[::-1] if self.scalingInv is not None: xd = self.scalingInv[d](xd) xmat[:, d] = kroneckerer(xd, nleft, nright) nleft *= n[d] nright = np.prod(n) if nright > 1 and reorder: fejerOrdering = [nright - 1] + lowDiscrepancy(nright - 1) xmat = xmat[fejerOrdering, :] x = checkParameterList(xmat, self.npar)[0] return x diff --git a/rrompy/parameter/parameter_sampling/quadrature_sampler_total.py b/rrompy/parameter/parameter_sampling/quadrature_sampler_total.py index ca49ca6..99d5cf3 100644 --- a/rrompy/parameter/parameter_sampling/quadrature_sampler_total.py +++ b/rrompy/parameter/parameter_sampling/quadrature_sampler_total.py @@ -1,60 +1,60 @@ # 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 . # import numpy as np from scipy.special import binom, factorial as fact from .quadrature_sampler import QuadratureSampler from rrompy.utilities.base.types import paramList -from rrompy.utilities.base import lowDiscrepancy +from rrompy.utilities.numerical import lowDiscrepancy from rrompy.parameter import checkParameterList __all__ = ['QuadratureSamplerTotal'] class QuadratureSamplerTotal(QuadratureSampler): """ Generator of quadrature sample points for total degree polynomial computations. """ def generatePoints(self, n:int, reorder : bool = True) -> paramList: """Array of sample points.""" if hasattr(n, "__len__"): n = n[0] d = self.npar n1d = int((fact(d) * n) ** (1. / d)) while binom(n1d + d - 1, d) > n: n1d -= 1 x = super().generatePoints([n1d] * d, reorder = False) nTot = n1d ** d indicesBase = np.zeros(nTot, dtype = int) idxBase = [x + 1 for x in lowDiscrepancy(n1d - 1, inverse = True)] linearIdxs = np.array(idxBase + [0]) nleft, nright = 1, nTot for j in range(d): nright //= n1d kronIn = np.repeat(linearIdxs, nright) indicesBase += np.tile(kronIn, nleft) nleft *= n1d keepIdxs = np.zeros(nTot, dtype = bool) keepIdxs[indicesBase < n1d] = True xmat = x.data[keepIdxs, :] if reorder: fejerTot = np.array([nTot - 1] + list(lowDiscrepancy(nTot - 1))) xmat = xmat[np.argsort(np.argsort(fejerTot[keepIdxs])), :] x = checkParameterList(xmat, d)[0] return x diff --git a/rrompy/parameter/parameter_sampling/random_sampler.py b/rrompy/parameter/parameter_sampling/random_sampler.py index 6193a4a..460d4e8 100644 --- a/rrompy/parameter/parameter_sampling/random_sampler.py +++ b/rrompy/parameter/parameter_sampling/random_sampler.py @@ -1,77 +1,76 @@ # 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 . # import numpy as np from .generic_sampler import GenericSampler -from rrompy.utilities.base.halton import haltonGenerate -from rrompy.utilities.base.sobol import sobolGenerate +from rrompy.utilities.numerical import haltonGenerate, sobolGenerate from rrompy.utilities.base.types import List, paramList from rrompy.utilities.exception_manager import RROMPyException from rrompy.parameter import checkParameterList __all__ = ['RandomSampler'] class RandomSampler(GenericSampler): """Generator of quadrature sample points.""" allowedKinds = ["UNIFORM", "HALTON", "SOBOL"] def __init__(self, lims:paramList, kind : str = "UNIFORM", scaling : List[callable] = None, scalingInv : List[callable] = None, seed : int = 42): super().__init__(lims = lims, scaling = scaling, scalingInv = scalingInv) self.kind = kind self.seed = seed def __str__(self) -> str: return "{}_{}".format(super().__str__(), self.kind) def __repr__(self) -> str: return self.__str__() + " at " + hex(id(self)) @property def kind(self): """Value of kind.""" return self._kind @kind.setter def kind(self, kind): if kind.upper() not in self.allowedKinds: raise RROMPyException("Generator kind not recognized.") self._kind = kind.upper() def generatePoints(self, n:int) -> paramList: """Array of quadrature points.""" if hasattr(n, "__len__"): n = n[0] if self.kind == "UNIFORM": np.random.seed(self.seed) xmat = np.random.uniform(size = (n, self.npar)) elif self.kind == "HALTON": xmat = haltonGenerate(self.npar, n, self.seed) else: xmat = sobolGenerate(self.npar, n, self.seed) for d in range(self.npar): a, b = self.lims(0, d), self.lims(1, d) if self.scaling is not None: a, b = self.scaling[d](a), self.scaling[d](b) xmat[:, d] = a + (b - a) * xmat[:, d] if self.scalingInv is not None: xmat[:, d] = self.scalingInv[d](xmat[:, d]) x = checkParameterList(xmat, self.npar)[0] return x diff --git a/rrompy/reduction_methods/distributed/rational_interpolant.py b/rrompy/reduction_methods/distributed/rational_interpolant.py index 68306aa..7d69bf0 100644 --- a/rrompy/reduction_methods/distributed/rational_interpolant.py +++ b/rrompy/reduction_methods/distributed/rational_interpolant.py @@ -1,550 +1,550 @@ # 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 . # from copy import deepcopy as copy import numpy as np from rrompy.reduction_methods.base import checkRobustTolerance from .generic_distributed_approximant import GenericDistributedApproximant -from rrompy.utilities.poly_fitting import customPInv from rrompy.utilities.poly_fitting.polynomial import (polybases, polyfitname, nextDerivativeIndices, hashDerivativeToIdx as hashD, hashIdxToDerivative as hashI, homogeneizedpolyvander as hpvP, homogeneizedToFull, PolynomialInterpolator as PI) from rrompy.utilities.poly_fitting.radial_basis import (rbbases, RadialBasisInterpolator as RBI) from rrompy.reduction_methods.trained_model import ( TrainedModelRational as tModel) from rrompy.reduction_methods.trained_model import TrainedModelData from rrompy.utilities.base.types import (Np1D, Np2D, HFEng, DictAny, Tuple, List, paramVal, paramList, sampList) -from rrompy.utilities.base import verbosityManager as vbMng, multifactorial +from rrompy.utilities.base import verbosityManager as vbMng +from rrompy.utilities.numerical import multifactorial, customPInv from rrompy.utilities.exception_manager import (RROMPyException, RROMPyAssert, RROMPyWarning) __all__ = ['RationalInterpolant'] class RationalInterpolant(GenericDistributedApproximant): """ ROM rational interpolant computation for parametric problems. Args: HFEngine: HF problem solver. mu0(optional): Default parameter. Defaults to 0. approxParameters(optional): Dictionary containing values for main parameters of approximant. Recognized keys are: - 'POD': whether to compute POD of snapshots; defaults to True; - 'S': total number of samples current approximant relies upon; - 'sampler': sample point generator; - 'polybasis': type of polynomial basis for interpolation; allowed values include 'MONOMIAL', 'CHEBYSHEV' and 'LEGENDRE'; defaults to 'MONOMIAL'; - 'E': number of derivatives used to compute interpolant; defaults to 0; - 'M': degree of rational interpolant numerator; defaults to 0; - 'N': degree of rational interpolant denominator; defaults to 0; - 'radialBasis': radial basis family for interpolant numerator; defaults to 0, i.e. no radial basis; - 'radialBasisWeights': radial basis weights for interpolant numerator; defaults to 0, i.e. identity; - 'interpRcond': tolerance for interpolation; defaults to None; - 'robustTol': tolerance for robust rational denominator management; defaults to 0. Defaults to empty dict. homogeneized(optional): Whether to homogeneize Dirichlet BCs. Defaults to False. verbosity(optional): Verbosity level. Defaults to 10. Attributes: HFEngine: HF problem solver. mu0: Default parameter. mus: Array of snapshot parameters. homogeneized: Whether to homogeneize Dirichlet BCs. approxParameters: Dictionary containing values for main parameters of approximant. Recognized keys are in parameterList. parameterListSoft: Recognized keys of soft approximant parameters: - 'POD': whether to compute POD of snapshots; - 'polybasis': type of polynomial basis for interpolation; - 'E': number of derivatives used to compute interpolant; - 'M': degree of rational interpolant numerator; - 'N': degree of rational interpolant denominator; - 'radialBasis': radial basis family for interpolant numerator; - 'radialBasisWeights': radial basis weights for interpolant numerator; - 'interpRcond': tolerance for interpolation via numpy.polyfit; - 'robustTol': tolerance for robust rational denominator management. parameterListCritical: Recognized keys of critical approximant parameters: - 'S': total number of samples current approximant relies upon; - 'sampler': sample point generator. POD: Whether to compute POD of snapshots. S: Number of solution snapshots over which current approximant is based upon. sampler: Sample point generator. polybasis: type of polynomial basis for interpolation. M: Numerator degree of approximant. N: Denominator degree of approximant. radialBasis: Radial basis family for interpolant numerator. radialBasisWeights: Radial basis weights for interpolant numerator. interpRcond: Tolerance for interpolation via numpy.polyfit. robustTol: Tolerance for robust rational denominator management. E: Complete derivative depth level of S. muBounds: list of bounds for parameter values. samplingEngine: Sampling engine. uHF: High fidelity solution(s) with parameter(s) lastSolvedHF as sampleList. lastSolvedHF: Parameter(s) corresponding to last computed high fidelity solution(s) as parameterList. uApproxReduced: Reduced approximate solution(s) with parameter(s) lastSolvedApprox as sampleList. lastSolvedApproxReduced: Parameter(s) corresponding to last computed reduced approximate solution(s) as parameterList. uApprox: Approximate solution(s) with parameter(s) lastSolvedApprox as sampleList. lastSolvedApprox: Parameter(s) corresponding to last computed approximate solution(s) as parameterList. Q: Numpy 1D vector containing complex coefficients of approximant denominator. P: Numpy 2D vector whose columns are FE dofs of coefficients of approximant numerator. """ def __init__(self, HFEngine:HFEng, mu0 : paramVal = None, approxParameters : DictAny = {}, homogeneized : bool = False, verbosity : int = 10, timestamp : bool = True): self._preInit() self._addParametersToList(["polybasis", "E", "M", "N", "radialBasis", "radialBasisWeights", "interpRcond", "robustTol"], ["MONOMIAL", -1, 0, 0, 0, 1, -1, 0]) super().__init__(HFEngine = HFEngine, mu0 = mu0, approxParameters = approxParameters, homogeneized = homogeneized, verbosity = verbosity, timestamp = timestamp) self._postInit() @property def polybasis(self): """Value of polybasis.""" return self._polybasis @polybasis.setter def polybasis(self, polybasis): try: polybasis = polybasis.upper().strip().replace(" ","") if polybasis not in polybases: raise RROMPyException("Prescribed polybasis not recognized.") self._polybasis = polybasis except: RROMPyWarning(("Prescribed polybasis not recognized. Overriding " "to 'MONOMIAL'.")) self._polybasis = "MONOMIAL" self._approxParameters["polybasis"] = self.polybasis @property def radialBasis(self): """Value of radialBasis.""" return self._radialBasis @radialBasis.setter def radialBasis(self, radialBasis): try: if radialBasis != 0: radialBasis = radialBasis.upper().strip().replace(" ","") if radialBasis not in rbbases: raise RROMPyException(("Prescribed radialBasis not " "recognized.")) self._radialBasis = radialBasis except: RROMPyWarning(("Prescribed radialBasis not recognized. Overriding " "to 0.")) self._radialBasis = 0 self._approxParameters["radialBasis"] = self.radialBasis @property def polybasisP(self): if self.radialBasis == 0: return self._polybasis return self._polybasis + "_" + self.radialBasis @property def interpRcond(self): """Value of interpRcond.""" return self._interpRcond @interpRcond.setter def interpRcond(self, interpRcond): self._interpRcond = interpRcond self._approxParameters["interpRcond"] = self.interpRcond @property def radialBasisWeights(self): """Value of radialBasisWeights.""" return self._radialBasisWeights @radialBasisWeights.setter def radialBasisWeights(self, radialBasisWeights): self._radialBasisWeights = radialBasisWeights self._approxParameters["radialBasisWeights"] = self.radialBasisWeights @property def M(self): """Value of M. Its assignment may change S.""" return self._M @M.setter def M(self, M): if M < 0: raise RROMPyException("M must be non-negative.") self._M = M self._approxParameters["M"] = self.M if hasattr(self, "_E") and self.E >= 0 and self.E < self.M: RROMPyWarning("Prescribed S is too small. Decreasing M.") self.M = self.E @property def N(self): """Value of N. Its assignment may change S.""" return self._N @N.setter def N(self, N): if N < 0: raise RROMPyException("N must be non-negative.") self._N = N self._approxParameters["N"] = self.N if hasattr(self, "_E") and self.E >= 0 and self.E < self.N: RROMPyWarning("Prescribed S is too small. Decreasing N.") self.N = self.E @property def E(self): """Value of E.""" return self._E @E.setter def E(self, E): if E < 0: if not hasattr(self, "_S"): raise RROMPyException(("Value of E must be positive if S is " "not yet assigned.")) E = np.sum(hashI(np.prod(self.S), self.npar)) - 1 self._E = E self._approxParameters["E"] = self.E if (hasattr(self, "_S") and self.E >= np.sum(hashI(np.prod(self.S), self.npar))): RROMPyWarning("Prescribed S is too small. Decreasing E.") self.E = -1 if hasattr(self, "_M"): self.M = self.M if hasattr(self, "_N"): self.N = self.N @property def robustTol(self): """Value of tolerance for robust rational denominator management.""" return self._robustTol @robustTol.setter def robustTol(self, robustTol): if robustTol < 0.: RROMPyWarning(("Overriding prescribed negative robustness " "tolerance to 0.")) robustTol = 0. self._robustTol = robustTol self._approxParameters["robustTol"] = self.robustTol @property def S(self): """Value of S.""" return self._S @S.setter def S(self, S): GenericDistributedApproximant.S.fset(self, S) if hasattr(self, "_M"): self.M = self.M if hasattr(self, "_N"): self.N = self.N if hasattr(self, "_E"): self.E = self.E def resetSamples(self): """Reset samples.""" super().resetSamples() self._musUniqueCN = None self._derIdxs = None self._reorder = None def _setupInterpolationIndices(self): """Setup parameters for polyvander.""" RROMPyAssert(self._mode, message = "Cannot setup interpolation indices.") if self._musUniqueCN is None or len(self._reorder) != len(self.mus): self._musUniqueCN, musIdxsTo, musIdxs, musCount = ( self.centerNormalize(self.mus).unique(return_index = True, return_inverse = True, return_counts = True)) self._musUnique = self.mus[musIdxsTo] self._derIdxs = [None] * len(self._musUniqueCN) self._reorder = np.empty(len(musIdxs), dtype = int) filled = 0 for j, cnt in enumerate(musCount): self._derIdxs[j] = nextDerivativeIndices([], self.mus.shape[1], cnt) jIdx = np.nonzero(musIdxs == j)[0] self._reorder[jIdx] = np.arange(filled, filled + cnt) filled += cnt def _setupDenominator(self): """Compute rational denominator.""" RROMPyAssert(self._mode, message = "Cannot setup denominator.") vbMng(self, "INIT", "Starting computation of denominator.", 7) while self.N > 0: invD = self._computeInterpolantInverseBlocks() if self.POD: ev, eV = self.findeveVGQR(self.samplingEngine.RPOD, invD) else: ev, eV = self.findeveVGExplicit(self.samplingEngine.samples, invD) newParams = checkRobustTolerance(ev, self.N, self.robustTol) if not newParams: break self.approxParameters = newParams if self.N <= 0: self._N = 0 eV = np.ones((1, 1)) q = PI() q.npar = self.npar q.polybasis = self.polybasis q.coeffs = homogeneizedToFull(tuple([self.N + 1] * self.npar), self.npar, eV[:, 0]) vbMng(self, "DEL", "Done computing denominator.", 7) return q def _setupNumerator(self): """Compute rational numerator.""" RROMPyAssert(self._mode, message = "Cannot setup numerator.") vbMng(self, "INIT", "Starting computation of numerator.", 7) Qevaldiag = np.zeros((len(self.mus), len(self.mus)), dtype = np.complex) verb = self.trainedModel.verbosity self.trainedModel.verbosity = 0 self._setupInterpolationIndices() idxGlob = 0 for j, derIdxs in enumerate(self._derIdxs): nder = len(derIdxs) idxLoc = np.arange(len(self.mus))[(self._reorder >= idxGlob) * (self._reorder < idxGlob + nder)] idxGlob += nder Qval = [0] * nder for der in range(nder): derIdx = hashI(der, self.npar) Qval[der] = (self.trainedModel.getQVal( self._musUnique[j], derIdx, scl = np.power(self.scaleFactor, -1.)) / multifactorial(derIdx)) for derU, derUIdx in enumerate(derIdxs): for derQ, derQIdx in enumerate(derIdxs): diffIdx = [x - y for (x, y) in zip(derUIdx, derQIdx)] if all([x >= 0 for x in diffIdx]): diffj = hashD(diffIdx) Qevaldiag[idxLoc[derU], idxLoc[derQ]] = Qval[diffj] if self.POD: Qevaldiag = Qevaldiag.dot(self.samplingEngine.RPOD.T) self.trainedModel.verbosity = verb while self.M >= 0: if self.radialBasis == 0: p = PI() wellCond, msg = p.setupByInterpolation( self._musUniqueCN, Qevaldiag, self.M, self.polybasisP, self.verbosity >= 5, True, {"derIdxs": self._derIdxs, "reorder": self._reorder, "scl": np.power(self.scaleFactor, -1.)}, {"rcond": self.interpRcond}) else: p = RBI() wellCond, msg = p.setupByInterpolation( self._musUniqueCN, Qevaldiag, self.M, self.polybasisP, self.radialBasisWeights, self.verbosity >= 5, True, {"derIdxs": self._derIdxs, "reorder": self._reorder, "scl": np.power(self.scaleFactor, -1.)}, {"rcond": self.interpRcond}) vbMng(self, "MAIN", msg, 5) if wellCond: break RROMPyWarning("Polyfit is poorly conditioned. Reducing M by 1.") self.M -= 1 if self.M < 0: raise RROMPyException(("Instability in computation of numerator. " "Aborting.")) vbMng(self, "DEL", "Done computing numerator.", 7) return p def setupApprox(self): """ Compute rational interpolant. SVD-based robust eigenvalue management. """ if self.checkComputedApprox(): return RROMPyAssert(self._mode, message = "Cannot setup approximant.") vbMng(self, "INIT", "Setting up {}.". format(self.name()), 5) self.computeScaleFactor() self.computeSnapshots() if self.trainedModel is None: self.trainedModel = tModel() self.trainedModel.verbosity = self.verbosity self.trainedModel.timestamp = self.timestamp data = TrainedModelData(self.trainedModel.name(), self.mu0, self.samplingEngine.samples, self.HFEngine.rescalingExp) data.scaleFactor = self.scaleFactor data.mus = copy(self.mus) self.trainedModel.data = data else: self.trainedModel = self.trainedModel self.trainedModel.data.projMat = copy(self.samplingEngine.samples) if self.N > 0: Q = self._setupDenominator() else: Q = PI() Q.coeffs = np.ones(tuple([1] * self.npar), dtype = np.complex) Q.npar = self.npar Q.polybasis = self.polybasis self.trainedModel.data.Q = Q self.trainedModel.data.P = self._setupNumerator() self.trainedModel.data.approxParameters = copy(self.approxParameters) vbMng(self, "DEL", "Done setting up approximant.", 5) def _computeInterpolantInverseBlocks(self) -> List[Np2D]: """ Compute inverse factors for minimal interpolant target functional. """ RROMPyAssert(self._mode, message = "Cannot solve eigenvalue problem.") self._setupInterpolationIndices() while self.E >= 0: eWidth = (hashD([self.E + 1] + [0] * (self.npar - 1)) - hashD([self.E] + [0] * (self.npar - 1))) TE, _, argIdxs = hpvP(self._musUniqueCN, self.E, self.polybasis, self._derIdxs, self._reorder, scl = np.power(self.scaleFactor, -1.)) fitOut = customPInv(TE[:, argIdxs], rcond = self.interpRcond, full = True) vbMng(self, "MAIN", ("Fitting {} samples with degree {} through {}... " "Conditioning of pseudoinverse system: {:.4e}.").format( TE.shape[0], self.E, polyfitname(self.polybasis), fitOut[1][1][0] / fitOut[1][1][-1]), 5) if fitOut[1][0] == len(argIdxs): self._fitinv = fitOut[0][- eWidth : , :] break RROMPyWarning("Polyfit is poorly conditioned. Reducing E by 1.") self.E -= 1 if self.E < 0: raise RROMPyException(("Instability in computation of " "denominator. Aborting.")) TN, _, argIdxs = hpvP(self._musUniqueCN, self.N, self.polybasis, self._derIdxs, self._reorder, scl = np.power(self.scaleFactor, -1.)) TN = TN[:, argIdxs] invD = [None] * (eWidth) for k in range(eWidth): pseudoInv = np.diag(self._fitinv[k, :]) idxGlob = 0 for j, derIdxs in enumerate(self._derIdxs): nder = len(derIdxs) idxGlob += nder if nder > 1: idxLoc = np.arange(len(self.mus))[ (self._reorder >= idxGlob - nder) * (self._reorder < idxGlob)] invLoc = self._fitinv[k, idxLoc] pseudoInv[np.ix_(idxLoc, idxLoc)] = 0. for diffj, diffjIdx in enumerate(derIdxs): for derQ, derQIdx in enumerate(derIdxs): derUIdx = [x - y for (x, y) in zip(diffjIdx, derQIdx)] if all([x >= 0 for x in derUIdx]): derU = hashD(derUIdx) pseudoInv[idxLoc[derU], idxLoc[derQ]] = ( invLoc[diffj]) invD[k] = pseudoInv.dot(TN) return invD def findeveVGExplicit(self, sampleE:sampList, invD:List[Np2D]) -> Tuple[Np1D, Np2D]: """ Compute explicitly eigenvalues and eigenvectors of rational denominator matrix. """ RROMPyAssert(self._mode, message = "Cannot solve eigenvalue problem.") nEnd = invD[0].shape[1] eWidth = len(invD) vbMng(self, "INIT", "Building gramian matrix.", 10) gramian = self.HFEngine.innerProduct(sampleE, sampleE) G = np.zeros((nEnd, nEnd), dtype = np.complex) for k in range(eWidth): G += invD[k].T.conj().dot(gramian.dot(invD[k])) vbMng(self, "DEL", "Done building gramian.", 10) vbMng(self, "INIT", "Solving eigenvalue problem for gramian matrix.", 7) ev, eV = np.linalg.eigh(G) vbMng(self, "MAIN", ("Solved eigenvalue problem of size {} with condition number " "{:.4e}.").format(nEnd, ev[-1] / ev[0]), 5) vbMng(self, "DEL", "Done solving eigenvalue problem.", 7) return ev, eV def findeveVGQR(self, RPODE:Np2D, invD:List[Np2D]) -> Tuple[Np1D, Np2D]: """ Compute eigenvalues and eigenvectors of rational denominator matrix through SVD of R factor. """ RROMPyAssert(self._mode, message = "Cannot solve eigenvalue problem.") nEnd = invD[0].shape[1] S = RPODE.shape[0] eWidth = len(invD) vbMng(self, "INIT", "Building half-gramian matrix stack.", 10) Rstack = np.zeros((S * eWidth, nEnd), dtype = np.complex) for k in range(eWidth): Rstack[k * S : (k + 1) * S, :] = RPODE.dot(invD[k]) vbMng(self, "DEL", "Done building half-gramian.", 10) vbMng(self, "INIT", "Solving svd for square root of gramian matrix.", 7) _, s, eV = np.linalg.svd(Rstack, full_matrices = False) ev = s[::-1] eV = eV[::-1, :].T.conj() vbMng(self, "MAIN", ("Solved svd problem of size {} x {} with condition number " "{:.4e}.").format(*Rstack.shape, s[0] / s[-1]), 5) vbMng(self, "DEL", "Done solving svd.", 7) return ev, eV def centerNormalize(self, mu : paramList = [], mu0 : paramVal = None) -> paramList: """ Compute normalized parameter to be plugged into approximant. Args: mu: Parameter(s) 1. mu0: Parameter(s) 2. If None, set to self.mu0. Returns: Normalized parameter. """ return self.trainedModel.centerNormalize(mu, mu0) def getResidues(self, *args, **kwargs) -> Np1D: """ Obtain approximant residues. Returns: Matrix with residues as columns. """ return self.trainedModel.getResidues(*args, **kwargs) diff --git a/rrompy/reduction_methods/distributed/rb_distributed.py b/rrompy/reduction_methods/distributed/rb_distributed.py index aa64217..ea97567 100644 --- a/rrompy/reduction_methods/distributed/rb_distributed.py +++ b/rrompy/reduction_methods/distributed/rb_distributed.py @@ -1,210 +1,203 @@ # 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 . # from copy import deepcopy as copy import numpy as np from .generic_distributed_approximant import GenericDistributedApproximant from rrompy.reduction_methods.trained_model import TrainedModelRB as tModel from rrompy.reduction_methods.trained_model import TrainedModelData from rrompy.reduction_methods.base.rb_utils import projectAffineDecomposition from rrompy.utilities.base.types import (Np1D, Np2D, List, Tuple, DictAny, HFEng, paramVal, sampList) from rrompy.utilities.base import verbosityManager as vbMng from rrompy.utilities.exception_manager import (RROMPyWarning, RROMPyException, RROMPyAssert) __all__ = ['RBDistributed'] class RBDistributed(GenericDistributedApproximant): """ ROM RB approximant computation for parametric problems. Args: HFEngine: HF problem solver. mu0(optional): Default parameter. Defaults to 0. approxParameters(optional): Dictionary containing values for main parameters of approximant. Recognized keys are: - 'POD': whether to compute POD of snapshots; defaults to True; - 'S': total number of samples current approximant relies upon; - 'sampler': sample point generator; - 'R': rank for Galerkin projection; defaults to prod(S); - 'PODTolerance': tolerance for snapshots POD; defaults to -1. Defaults to empty dict. homogeneized(optional): Whether to homogeneize Dirichlet BCs. Defaults to False. verbosity(optional): Verbosity level. Defaults to 10. Attributes: HFEngine: HF problem solver. mu0: Default parameter. mus: Array of snapshot parameters. homogeneized: Whether to homogeneize Dirichlet BCs. approxRadius: Dummy radius of approximant (i.e. distance from mu0 to farthest sample point). approxParameters: Dictionary containing values for main parameters of approximant. Recognized keys are in parameterList. parameterListSoft: Recognized keys of soft approximant parameters: - 'POD': whether to compute POD of snapshots. - 'R': rank for Galerkin projection; - 'PODTolerance': tolerance for snapshots POD. parameterListCritical: Recognized keys of critical approximant parameters: - 'S': total number of samples current approximant relies upon; - 'sampler': sample point generator; POD: Whether to compute POD of snapshots. S: Number of solution snapshots over which current approximant is based upon. sampler: Sample point generator. R: Rank for Galerkin projection. muBounds: list of bounds for parameter values. samplingEngine: Sampling engine. uHF: High fidelity solution(s) with parameter(s) lastSolvedHF as sampleList. lastSolvedHF: Parameter(s) corresponding to last computed high fidelity solution(s) as parameterList. uApproxReduced: Reduced approximate solution(s) with parameter(s) lastSolvedApprox as sampleList. lastSolvedApproxReduced: Parameter(s) corresponding to last computed reduced approximate solution(s) as parameterList. uApprox: Approximate solution(s) with parameter(s) lastSolvedApprox as sampleList. lastSolvedApprox: Parameter(s) corresponding to last computed approximate solution(s) as parameterList. As: List of sparse matrices (in CSC format) representing coefficients - of linear system matrix wrt theta(mu). + of linear system matrix. bs: List of numpy vectors representing coefficients of linear system - RHS wrt theta(mu). - thetaAs: List of callables representing coefficients of linear system - matrix wrt mu. - thetabs: List of callables representing coefficients of linear system - RHS wrt mu. + RHS. ARBs: List of sparse matrices (in CSC format) representing coefficients - of compressed linear system matrix wrt theta(mu). + of compressed linear system matrix. bRBs: List of numpy vectors representing coefficients of compressed - linear system RHS wrt theta(mu). + linear system RHS. """ def __init__(self, HFEngine:HFEng, mu0 : paramVal = None, approxParameters : DictAny = {}, homogeneized : bool = False, verbosity : int = 10, timestamp : bool = True): self._preInit() self._addParametersToList(["R", "PODTolerance"], [1, -1]) super().__init__(HFEngine = HFEngine, mu0 = mu0, approxParameters = approxParameters, homogeneized = homogeneized, verbosity = verbosity, timestamp = timestamp) vbMng(self, "INIT", "Computing affine blocks of system.", 10) self.As = self.HFEngine.affineLinearSystemA(self.mu0) self.bs = self.HFEngine.affineLinearSystemb(self.mu0, self.homogeneized) vbMng(self, "DEL", "Done computing affine blocks.", 10) self._postInit() @property def S(self): """Value of S.""" return self._S @S.setter def S(self, S): GenericDistributedApproximant.S.fset(self, S) if not hasattr(self, "_R"): self._R = np.prod(self.S) @property def R(self): """Value of R. Its assignment may change S.""" return self._R @R.setter def R(self, R): if R < 0: raise RROMPyException("R must be non-negative.") self._R = R self._approxParameters["R"] = self.R if hasattr(self, "_S") and np.prod(self.S) < self.R: RROMPyWarning("Prescribed S is too small. Decreasing R.") self.R = np.prod(self.S) @property def PODTolerance(self): """Value of PODTolerance.""" return self._PODTolerance @PODTolerance.setter def PODTolerance(self, PODTolerance): self._PODTolerance = PODTolerance self._approxParameters["PODTolerance"] = self.PODTolerance def setupApprox(self): """Compute RB projection matrix.""" if self.checkComputedApprox(): return RROMPyAssert(self._mode, message = "Cannot setup approximant.") vbMng(self, "INIT", "Setting up {}.". format(self.name()), 5) self.computeSnapshots() vbMng(self, "INIT", "Computing projection matrix.", 7) if self.POD: U, s, _ = np.linalg.svd(self.samplingEngine.RPOD) s = s ** 2. else: Gramian = self.HFEngine.innerProduct(self.samplingEngine.samples, self.samplingEngine.samples) U, s, _ = np.linalg.svd(Gramian) nsamples = self.samplingEngine.nsamples snorm = np.cumsum(s[::-1]) / np.sum(s) nPODTrunc = min(nsamples - np.argmax(snorm > self.PODTolerance), self.R) pMat = self.samplingEngine.samples.dot(U[:, : nPODTrunc]) vbMng(self, "MAIN", ("Assembling {}x{} projection matrix from {} " "samples.").format(*(pMat.shape), nsamples), 5) if self.trainedModel is None: self.trainedModel = tModel() self.trainedModel.verbosity = self.verbosity self.trainedModel.timestamp = self.timestamp data = TrainedModelData(self.trainedModel.name(), self.mu0, pMat, self.HFEngine.rescalingExp) - data.thetaAs = self.HFEngine.affineWeightsA(self.mu0) - data.thetabs = self.HFEngine.affineWeightsb(self.mu0, - self.homogeneized) data.mus = copy(self.mus) self.trainedModel.data = data else: self.trainedModel = self.trainedModel self.trainedModel.data.projMat = copy(pMat) ARBs, bRBs = self.assembleReducedSystem(pMat) self.trainedModel.data.ARBs = ARBs self.trainedModel.data.bRBs = bRBs vbMng(self, "DEL", "Done computing projection matrix.", 7) self.trainedModel.data.approxParameters = copy(self.approxParameters) vbMng(self, "DEL", "Done setting up approximant.", 5) def assembleReducedSystem(self, pMat : sampList = None, pMatOld : sampList = None)\ -> Tuple[List[Np2D], List[Np1D]]: """Build affine blocks of RB linear system through projections.""" if pMat is None: self.setupApprox() ARBs = self.trainedModel.data.ARBs bRBs = self.trainedModel.data.bRBs else: vbMng(self, "INIT", "Projecting affine terms of HF model.", 10) ARBsOld = None if pMatOld is None else self.trainedModel.data.ARBs bRBsOld = None if pMatOld is None else self.trainedModel.data.bRBs ARBs, bRBs = projectAffineDecomposition(self.As, self.bs, pMat, ARBsOld, bRBsOld, pMatOld) vbMng(self, "DEL", "Done projecting affine terms.", 10) return ARBs, bRBs diff --git a/rrompy/reduction_methods/distributed_greedy/rb_distributed_greedy.py b/rrompy/reduction_methods/distributed_greedy/rb_distributed_greedy.py index a80695b..1529dc1 100644 --- a/rrompy/reduction_methods/distributed_greedy/rb_distributed_greedy.py +++ b/rrompy/reduction_methods/distributed_greedy/rb_distributed_greedy.py @@ -1,241 +1,237 @@ # 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 . # import numpy as np from copy import deepcopy as copy from .generic_distributed_greedy_approximant import \ GenericDistributedGreedyApproximant from rrompy.reduction_methods.distributed import RBDistributed from rrompy.reduction_methods.trained_model import TrainedModelRB as tModel from rrompy.reduction_methods.trained_model import TrainedModelData from rrompy.utilities.base.types import Np1D, DictAny, HFEng, paramVal from rrompy.utilities.base import verbosityManager as vbMng +from rrompy.utilities.poly_fitting.polynomial import ( + hashIdxToDerivative as hashI) from rrompy.utilities.exception_manager import RROMPyWarning, RROMPyAssert from rrompy.parameter import checkParameterList __all__ = ['RBDistributedGreedy'] class RBDistributedGreedy(GenericDistributedGreedyApproximant, RBDistributed): """ ROM greedy RB approximant computation for parametric problems. Args: HFEngine: HF problem solver. mu0(optional): Default parameter. Defaults to 0. approxParameters(optional): Dictionary containing values for main parameters of approximant. Recognized keys are: - 'POD': whether to compute POD of snapshots; defaults to True; - 'S': number of starting training points; - 'sampler': sample point generator; - 'greedyTol': uniform error tolerance for greedy algorithm; defaults to 1e-2; - 'interactive': whether to interactively terminate greedy algorithm; defaults to False; - 'maxIter': maximum number of greedy steps; defaults to 1e2; - 'refinementRatio': ratio of test points to be exhausted before test set refinement; defaults to 0.2; - 'nTestPoints': number of test points; defaults to 5e2; - 'trainSetGenerator': training sample points generator; defaults to Chebyshev sampler within muBounds. Defaults to empty dict. homogeneized(optional): Whether to homogeneize Dirichlet BCs. Defaults to False. verbosity(optional): Verbosity level. Defaults to 10. Attributes: HFEngine: HF problem solver. mu0: Default parameter. mus: Array of snapshot parameters. homogeneized: Whether to homogeneize Dirichlet BCs. approxParameters: Dictionary containing values for main parameters of approximant. Recognized keys are in parameterList. parameterListSoft: Recognized keys of soft approximant parameters: - 'POD': whether to compute POD of snapshots. - 'greedyTol': uniform error tolerance for greedy algorithm; - 'interactive': whether to interactively terminate greedy algorithm; - 'maxIter': maximum number of greedy steps; - 'refinementRatio': ratio of test points to be exhausted before test set refinement; - 'nTestPoints': number of test points; - 'trainSetGenerator': training sample points generator. parameterListCritical: Recognized keys of critical approximant parameters: - 'S': total number of samples current approximant relies upon; - 'sampler': sample point generator. POD: whether to compute POD of snapshots. S: number of test points. sampler: Sample point generator. greedyTol: uniform error tolerance for greedy algorithm. interactive: whether to interactively terminate greedy algorithm. maxIter: maximum number of greedy steps. refinementRatio: ratio of training points to be exhausted before training set refinement. nTestPoints: number of starting training points. trainSetGenerator: training sample points generator. muBounds: list of bounds for parameter values. samplingEngine: Sampling engine. estimatorNormEngine: Engine for estimator norm computation. uHF: High fidelity solution(s) with parameter(s) lastSolvedHF as sampleList. lastSolvedHF: Parameter(s) corresponding to last computed high fidelity solution(s) as parameterList. uApproxReduced: Reduced approximate solution(s) with parameter(s) lastSolvedApprox as sampleList. lastSolvedApproxReduced: Parameter(s) corresponding to last computed reduced approximate solution(s) as parameterList. uApprox: Approximate solution(s) with parameter(s) lastSolvedApprox as sampleList. lastSolvedApprox: Parameter(s) corresponding to last computed approximate solution(s) as parameterList. As: List of sparse matrices (in CSC format) representing coefficients - of linear system matrix wrt theta(mu). + of linear system matrix. bs: List of numpy vectors representing coefficients of linear system - RHS wrt theta(mu). - thetaAs: List of callables representing coefficients of linear system - matrix wrt mu. - thetabs: List of callables representing coefficients of linear system - RHS wrt mu. + RHS. ARBs: List of sparse matrices (in CSC format) representing coefficients - of compressed linear system matrix wrt theta(mu). + of compressed linear system matrix. bRBs: List of numpy vectors representing coefficients of compressed - linear system RHS wrt theta(mu). + linear system RHS. """ def __init__(self, HFEngine:HFEng, mu0 : paramVal = None, approxParameters : DictAny = {}, homogeneized : bool = False, verbosity : int = 10, timestamp : bool = True): self._preInit() super().__init__(HFEngine = HFEngine, mu0 = mu0, approxParameters = approxParameters, homogeneized = homogeneized, verbosity = verbosity, timestamp = timestamp) vbMng(self, "INIT", "Computing affine blocks of system.", 10) self.As = self.HFEngine.affineLinearSystemA(self.mu0) self.bs = self.HFEngine.affineLinearSystemb(self.mu0, self.homogeneized) vbMng(self, "DEL", "Done computing affine blocks.", 10) self._postInit() @property def R(self): """Value of R.""" self._R = np.prod(self._S) return self._R @R.setter def R(self, R): RROMPyWarning(("R is used just to simplify inheritance, and its value " "cannot be changed from that of prod(S).")) @property def PODTolerance(self): """Value of PODTolerance.""" self._PODTolerance = -1 return self._PODTolerance @PODTolerance.setter def PODTolerance(self, PODTolerance): RROMPyWarning(("PODTolerance is used just to simplify inheritance, " "and its value cannot be changed from -1.")) def errorEstimator(self, mus:Np1D) -> Np1D: """ Standard residual-based error estimator. Unreliable for unstable problems (inf-sup constant is missing). """ self.setupApprox() self.assembleReducedResidualBlocks(full = True) nmus = len(mus) nAs = self.trainedModel.data.resAA.shape[1] nbs = self.trainedModel.data.resbb.shape[0] - thetaAs = self.trainedModel.data.thetaAs - thetabs = self.trainedModel.data.thetabs radiusA = np.empty((len(self.mus), nAs, nmus), dtype = np.complex) radiusb = np.empty((nbs, nmus), dtype = np.complex) verb = self.trainedModel.verbosity self.trainedModel.verbosity = 0 if verb >= 5: mustr = mus if nmus > 2: mustr = "[{} ..({}).. {}]".format(mus[0], nmus - 2, mus[-1]) vbMng(self, "INIT", "Computing RB solution at mu = {}.".format(mustr), 0) parmus = checkParameterList(mus, self.npar)[0] uApproxRs = self.getApproxReduced(parmus) + mu0Eff = self.mu0(0, 0) ** self.HFEngine.rescalingExp[0] for j, muPL in enumerate(parmus): - mu = muPL[0] uApproxR = uApproxRs[j] - for i in range(nAs): - radiusA[:, i, j] = eval(thetaAs[i]) * uApproxR - for i in range(nbs): - radiusb[i, j] = eval(thetabs[i]) + for i in range(max(nAs, nbs)): + thetai = (muPL[0] ** self.HFEngine.rescalingExp[0] + - mu0Eff) ** i + if i < nAs: + radiusA[:, i, j] = thetai * uApproxR + if i < nbs: + radiusb[i, j] = thetai vbMng(self, "DEL", "Done computing RB solution.", 5) self.trainedModel.verbosity = verb # 'ij,jk,ik->k', resbb, radiusb, radiusb.conj() ff = np.sum(self.trainedModel.data.resbb.dot(radiusb) * radiusb.conj(), axis = 0) # 'ijk,jkl,il->l', resAb, radiusA, radiusb.conj() Lf = np.sum(np.tensordot(self.trainedModel.data.resAb, radiusA, 2) * radiusb.conj(), axis = 0) # 'ijkl,klt,ijt->t', resAA, radiusA, radiusA.conj() LL = np.sum(np.tensordot(self.trainedModel.data.resAA, radiusA, 2) * radiusA.conj(), axis = (0, 1)) return np.abs((LL - 2. * np.real(Lf) + ff) / ff) ** .5 def setupApprox(self, plotEst : bool = False): """Compute RB projection matrix.""" if self.checkComputedApprox(): return RROMPyAssert(self._mode, message = "Cannot setup approximant.") vbMng(self, "INIT", "Setting up {}.". format(self.name()), 5) self.greedy(plotEst) vbMng(self, "INIT", "Computing projection matrix.", 7) pMat = self.samplingEngine.samples if self.trainedModel is None: self.trainedModel = tModel() self.trainedModel.verbosity = self.verbosity self.trainedModel.timestamp = self.timestamp data = TrainedModelData(self.trainedModel.name(), self.mu0, pMat, self.HFEngine.rescalingExp) - data.thetaAs = self.HFEngine.affineWeightsA(self.mu0) - data.thetabs = self.HFEngine.affineWeightsb(self.mu0, - self.homogeneized) ARBs, bRBs = self.assembleReducedSystem(pMat) self.trainedModel.data = data else: self.trainedModel = self.trainedModel pMatOld = self.trainedModel.data.projMat Sold = pMatOld.shape[1] idxNew = list(range(Sold, pMat.shape[1])) ARBs, bRBs = self.assembleReducedSystem(pMat(idxNew), pMatOld) self.trainedModel.data.projMat = copy(pMat) self.trainedModel.data.mus = copy(self.mus) self.trainedModel.data.ARBs = ARBs self.trainedModel.data.bRBs = bRBs vbMng(self, "DEL", "Done computing projection matrix.", 7) self.trainedModel.data.approxParameters = copy(self.approxParameters) vbMng(self, "DEL", "Done setting up approximant.", 5) def assembleReducedResidualBlocks(self, full : bool = False): """Build affine blocks of RB linear system through projections.""" self.assembleReducedResidualBlocksbb(self.bs) if full: pMat = self.trainedModel.data.projMat self.assembleReducedResidualBlocksAb(self.As, self.bs, pMat) self.assembleReducedResidualBlocksAA(self.As, pMat) diff --git a/rrompy/reduction_methods/pivoting/__init__.py b/rrompy/reduction_methods/pivoting/__init__.py index 9830125..004b3fa 100644 --- a/rrompy/reduction_methods/pivoting/__init__.py +++ b/rrompy/reduction_methods/pivoting/__init__.py @@ -1,29 +1,35 @@ # 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 . # from .generic_pivoted_approximant import GenericPivotedApproximant from .rational_interpolant_pivoted import RationalInterpolantPivoted from .rb_distributed_pivoted import RBDistributedPivoted +from .rational_interpolant_pivoted_pole_matching import \ + RationalInterpolantPivotedPoleMatching +from .rb_distributed_pivoted_pole_matching import \ + RBDistributedPivotedPoleMatching __all__ = [ 'GenericPivotedApproximant', 'RationalInterpolantPivoted', - 'RBDistributedPivoted' + 'RBDistributedPivoted', + 'RationalInterpolantPivotedPoleMatching', + 'RBDistributedPivotedPoleMatching' ] diff --git a/rrompy/reduction_methods/pivoting/generic_pivoted_approximant.py b/rrompy/reduction_methods/pivoting/generic_pivoted_approximant.py index d9c2cc6..57576d0 100644 --- a/rrompy/reduction_methods/pivoting/generic_pivoted_approximant.py +++ b/rrompy/reduction_methods/pivoting/generic_pivoted_approximant.py @@ -1,481 +1,480 @@ # 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 . # import numpy as np from copy import deepcopy as copy from rrompy.reduction_methods.base.generic_approximant import ( GenericApproximant) -from rrompy.hfengines.base import MarginalProxyEngine from rrompy.utilities.poly_fitting.polynomial import (polybases, nextDerivativeIndices, PolynomialInterpolator as PI) from rrompy.utilities.poly_fitting.radial_basis import ( RadialBasisInterpolator as RBI) from rrompy.utilities.poly_fitting.radial_basis import rbbases from rrompy.sampling.pivoted import (SamplingEnginePivoted, SamplingEnginePivotedPOD) from rrompy.utilities.base.types import (ListAny, DictAny, HFEng, paramVal, paramList) -from rrompy.utilities.base import verbosityManager as vbMng, freepar as fp +from rrompy.utilities.base import verbosityManager as vbMng from rrompy.utilities.exception_manager import (RROMPyException, RROMPyAssert, RROMPyWarning) from rrompy.parameter import emptyParameterList __all__ = ['GenericPivotedApproximant'] class GenericPivotedApproximant(GenericApproximant): """ ROM pivoted approximant computation for parametric problems (ABSTRACT). Args: HFEngine: HF problem solver. mu0(optional): Default parameter. Defaults to 0. directionPivot(optional): Pivot components. Defaults to [0]. approxParameters(optional): Dictionary containing values for main parameters of approximant. Recognized keys are: - 'POD': whether to compute POD of snapshots; defaults to True; - 'S': total number of pivot samples current approximant relies upon; - 'samplerPivot': pivot sample point generator; - 'SMarginal': total number of marginal samples current approximant relies upon; - 'samplerMarginal': marginal sample point generator; - 'polybasisMarginal': type of polynomial basis for marginal interpolation; allowed values include 'MONOMIAL', 'CHEBYSHEV' and 'LEGENDRE'; defaults to 'MONOMIAL'; - 'MMarginal': degree of marginal interpolant; defaults to 0; - 'radialBasisMarginal': radial basis family for marginal interpolant; defaults to 0, i.e. no radial basis; - 'radialBasisWeightsMarginal': radial basis weights for marginal interpolant; defaults to 0, i.e. identity; - 'interpRcondMarginal': tolerance for marginal interpolation; defaults to None. Defaults to empty dict. homogeneized(optional): Whether to homogeneize Dirichlet BCs. Defaults to False. verbosity(optional): Verbosity level. Defaults to 10. Attributes: HFEngine: HF problem solver. mu0: Default parameter. directionPivot: Pivot components. mus: Array of snapshot parameters. musPivot: Array of pivot snapshot parameters. musMarginal: Array of marginal snapshot parameters. homogeneized: Whether to homogeneize Dirichlet BCs. approxParameters: Dictionary containing values for main parameters of approximant. Recognized keys are in parameterList. parameterListSoft: Recognized keys of soft approximant parameters: - 'POD': whether to compute POD of snapshots; - 'polybasisMarginal': type of polynomial basis for marginal interpolation; - 'MMarginal': degree of marginal interpolant; - 'radialBasisMarginal': radial basis family for marginal interpolant; - 'radialBasisWeightsMarginal': radial basis weights for marginal interpolant; - 'interpRcondMarginal': tolerance for marginal interpolation. parameterListCritical: Recognized keys of critical approximant parameters: - 'S': total number of pivot samples current approximant relies upon; - 'samplerPivot': pivot sample point generator; - 'SMarginal': total number of marginal samples current approximant relies upon; - 'samplerMarginal': marginal sample point generator. POD: Whether to compute POD of snapshots. S: Total number of pivot samples current approximant relies upon. samplerPivot: Pivot sample point generator. SMarginal: Total number of marginal samples current approximant relies upon. samplerMarginal: Marginal sample point generator. polybasisMarginal: Type of polynomial basis for marginal interpolation. MMarginal: Degree of marginal interpolant. radialBasisMarginal: Radial basis family for marginal interpolant. radialBasisWeightsMarginal: Radial basis weights for marginal interpolant. interpRcondMarginal: Tolerance for marginal interpolation. muBoundsPivot: list of bounds for pivot parameter values. muBoundsMarginal: list of bounds for marginal parameter values. samplingEngine: Sampling engine. uHF: High fidelity solution(s) with parameter(s) lastSolvedHF as sampleList. lastSolvedHF: Parameter(s) corresponding to last computed high fidelity solution(s) as parameterList. uApproxReduced: Reduced approximate solution(s) with parameter(s) lastSolvedApprox as sampleList. lastSolvedApproxReduced: Parameter(s) corresponding to last computed reduced approximate solution(s) as parameterList. uApprox: Approximate solution(s) with parameter(s) lastSolvedApprox as sampleList. lastSolvedApprox: Parameter(s) corresponding to last computed approximate solution(s) as parameterList. """ def __init__(self, HFEngine:HFEng, mu0 : paramVal = None, directionPivot : ListAny = [0], approxParameters : DictAny = {}, homogeneized : bool = False, verbosity : int = 10, timestamp : bool = True): self._preInit() from rrompy.parameter.parameter_sampling import QuadratureSampler as QS QSBase = QS([[0], [1]], "UNIFORM") self._addParametersToList( ["polybasisMarginal", "MMarginal", "radialBasisMarginal", "radialBasisWeightsMarginal", "interpRcondMarginal"], ["MONOMIAL", 0, 0, 1, -1], ["samplerPivot", "SMarginal", "samplerMarginal"], [QSBase, [1], QSBase]) del QS self._directionPivot = directionPivot super().__init__(HFEngine = HFEngine, mu0 = mu0, approxParameters = approxParameters, homogeneized = homogeneized, verbosity = verbosity, timestamp = timestamp) self._postInit() def setupSampling(self): """Setup sampling engine.""" RROMPyAssert(self._mode, message = "Cannot setup sampling engine.") if not hasattr(self, "_POD") or self._POD is None: return if self.POD: SamplingEngine = SamplingEnginePivotedPOD else: SamplingEngine = SamplingEnginePivoted self.samplingEngine = SamplingEngine(self.HFEngine, self.directionPivot, verbosity = self.verbosity, allowRepeatedSamples = True) @property def SMarginal(self): """Value of SMarginal.""" return self._SMarginal @SMarginal.setter def SMarginal(self, SMarginal): if not hasattr(SMarginal, "__len__"): SMarginal = [SMarginal] if any([s <= 0 for s in SMarginal]): raise RROMPyException("SMarginal must be positive.") if hasattr(self, "_SMarginal") and self._SMarginal is not None: Sold = tuple(self.SMarginal) else: Sold = -1 self._SMarginal = SMarginal self._approxParameters["SMarginal"] = self.SMarginal if Sold != tuple(self.SMarginal): self.resetSamples() @property def polybasisMarginal(self): """Value of polybasisMarginal.""" return self._polybasisMarginal @polybasisMarginal.setter def polybasisMarginal(self, polybasisMarginal): try: polybasisMarginal = polybasisMarginal.upper().strip().replace(" ", "") if polybasisMarginal not in polybases: raise RROMPyException( "Prescribed marginal polybasis not recognized.") self._polybasisMarginal = polybasisMarginal except: RROMPyWarning(("Prescribed marginal polybasis not recognized. " "Overriding to 'MONOMIAL'.")) self._polybasisMarginal = "MONOMIAL" self._approxParameters["polybasisMarginal"] = self.polybasisMarginal @property def MMarginal(self): """Value of MMarginal.""" return self._MMarginal @MMarginal.setter def MMarginal(self, MMarginal): if MMarginal < 0: raise RROMPyException("MMarginal must be non-negative.") self._MMarginal = MMarginal self._approxParameters["MMarginal"] = self.MMarginal @property def radialBasisMarginal(self): """Value of radialBasisMarginal.""" return self._radialBasisMarginal @radialBasisMarginal.setter def radialBasisMarginal(self, radialBasisMarginal): try: if radialBasisMarginal != 0: radialBasisMarginal = radialBasisMarginal.upper().strip()\ .replace(" ","") if radialBasisMarginal not in rbbases: raise RROMPyException(("Prescribed marginal radialBasis " "not recognized.")) self._radialBasisMarginal = radialBasisMarginal except: RROMPyWarning(("Prescribed marginal radialBasis not recognized. " "Overriding to 0.")) self._radialBasisMarginal = 0 self._approxParameters["radialBasisMarginal"] = ( self.radialBasisMarginal) @property def radialBasisWeightsMarginal(self): """Value of radialBasisWeightsMarginal.""" return self._radialBasisWeightsMarginal @radialBasisWeightsMarginal.setter def radialBasisWeightsMarginal(self, radialBasisWeightsMarginal): self._radialBasisWeightsMarginal = radialBasisWeightsMarginal self._approxParameters["radialBasisWeightsMarginal"] = ( self.radialBasisWeightsMarginal) @property def polybasisMarginalP(self): if self.radialBasisMarginal == 0: return self._polybasisMarginal return self._polybasisMarginal + "_" + self.radialBasisMarginal @property def interpRcondMarginal(self): """Value of interpRcondMarginal.""" return self._interpRcondMarginal @interpRcondMarginal.setter def interpRcondMarginal(self, interpRcondMarginal): self._interpRcondMarginal = interpRcondMarginal self._approxParameters["interpRcondMarginal"] = ( self.interpRcondMarginal) @property def directionPivot(self): """Value of directionPivot. Its assignment may reset snapshots.""" return self._directionPivot @directionPivot.setter def directionPivot(self, directionPivot): if hasattr(self, '_directionPivot'): directionPivotOld = copy(self.directionPivot) else: directionPivotOld = None if (directionPivotOld is None or len(directionPivot) != len(directionPivotOld) or not directionPivot == directionPivotOld): self.resetSamples() self._directionPivot = directionPivot @property def directionMarginal(self): return [x for x in range(self.HFEngine.npar) \ if x not in self.directionPivot] @property def rescalingExpPivot(self): return [self.HFEngine.rescalingExp[x] for x in self.directionPivot] @property def rescalingExpMarginal(self): return [self.HFEngine.rescalingExp[x] for x in self.directionMarginal] @property def muBoundsPivot(self): """Value of muBoundsPivot.""" return self.samplerPivot.lims @property def muBoundsMarginal(self): """Value of muBoundsMarginal.""" return self.samplerMarginal.lims @property def samplerPivot(self): """Value of samplerPivot.""" return self._samplerPivot @samplerPivot.setter def samplerPivot(self, samplerPivot): if 'generatePoints' not in dir(samplerPivot): raise RROMPyException("Pivot sampler type not recognized.") if hasattr(self, '_samplerPivot') and self._samplerPivot is not None: samplerOld = self.samplerPivot self._samplerPivot = samplerPivot self._approxParameters["samplerPivot"] = self.samplerPivot.__str__() if not 'samplerOld' in locals() or samplerOld != self.samplerPivot: self.resetSamples() @property def samplerMarginal(self): """Value of samplerMarginal.""" return self._samplerMarginal @samplerMarginal.setter def samplerMarginal(self, samplerMarginal): if 'generatePoints' not in dir(samplerMarginal): raise RROMPyException("Marginal sampler type not recognized.") if (hasattr(self, '_samplerMarginal') and self._samplerMarginal is not None): samplerOld = self.samplerMarginal self._samplerMarginal = samplerMarginal self._approxParameters["samplerMarginal"] = ( self.samplerMarginal.__str__()) if not 'samplerOld' in locals() or samplerOld != self.samplerMarginal: self.resetSamples() def resetSamples(self): """Reset samples.""" super().resetSamples() self._musMUniqueCN = None self._derMIdxs = None self._reorderM = None def setSamples(self, samplingEngine): """Copy samplingEngine and samples.""" super().setSamples(samplingEngine) self.mus = copy(self.samplingEngine[0].mus) for sEj in self.samplingEngine[1:]: self.mus.append(sEj.mus) def computeSnapshots(self): """Compute snapshots of solution map.""" RROMPyAssert(self._mode, message = "Cannot start snapshot computation.") Sp = np.prod(self.S) Seff = Sp * np.prod(self.SMarginal) if self.samplingEngine.nsamplesTot != Seff: self.resetSamples() vbMng(self, "INIT", "Starting computation of snapshots.", 5) if self.HFEngine.As[0] is None: self.HFEngine.A(self.mu0) if self.HFEngine.bs[0] is None: self.HFEngine.b(self.mu0) self.musPivot = self.samplerPivot.generatePoints(self.S) self.musMarginal = self.samplerMarginal.generatePoints( self.SMarginal) self.mus = emptyParameterList() self.mus.reset((Seff, self.HFEngine.npar)) self.samplingEngine.resetHistory(Seff // Sp) for j, muMarg in enumerate(self.musMarginal): for k in range(j * Sp, (j + 1) * Sp): self.mus.data[k, self.directionPivot] = ( self.musPivot[k - j * Sp].data) self.mus.data[k, self.directionMarginal] = muMarg.data self.samplingEngine.iterSample(self.musPivot, self.musMarginal, homogeneized = self.homogeneized) if self.POD: self.samplingEngine.coalesceSamples(self.interpRcondMarginal) else: self.samplingEngine.coalesceSamples() vbMng(self, "DEL", "Done computing snapshots.", 5) def _setupMarginalInterpolationIndices(self): """Setup parameters for polyvander.""" RROMPyAssert(self._mode, message = "Cannot setup interpolation indices.") if (self._musMUniqueCN is None or len(self._reorderM) != len(self.musMarginal)): self._musMUniqueCN, musMIdxsTo, musMIdxs, musMCount = ( self.centerNormalizeMarginal(self.musMarginal).unique( return_index = True, return_inverse = True, return_counts = True)) self._musMUnique = self.musMarginal[musMIdxsTo] self._derMIdxs = [None] * len(self._musMUniqueCN) self._reorderM = np.empty(len(musMIdxs), dtype = int) filled = 0 for j, cnt in enumerate(musMCount): self._derMIdxs[j] = nextDerivativeIndices([], self.musMarginal.shape[1], cnt) jIdx = np.nonzero(musMIdxs == j)[0] self._reorderM[jIdx] = np.arange(filled, filled + cnt) filled += cnt def _setupMarginalInterp(self): """Compute marginal interpolator.""" RROMPyAssert(self._mode, message = "Cannot setup numerator.") vbMng(self, "INIT", "Starting computation of marginal interpolator.", 7) self._setupMarginalInterpolationIndices() MMarginal0 = copy(self.MMarginal) mI = [] for j in range(len(self.musMarginal)): canonicalj = 1. * (np.arange(len(self.musMarginal)) == j) self._MMarginal = MMarginal0 while self.MMarginal >= 0: if self.radialBasisMarginal == 0: p = PI() wellCond, msg = p.setupByInterpolation( self._musMUniqueCN, canonicalj, self.MMarginal, self.polybasisMarginal, self.verbosity >= 5, True, {"derIdxs": self._derMIdxs, "reorder": self._reorderM, "scl": np.power(self.scaleFactorMarginal, -1.)}, {"rcond": self.interpRcondMarginal}) else: p = RBI() wellCond, msg = p.setupByInterpolation( self._musMUniqueCN, canonicalj, self.MMarginal, self.polybasisMarginalP, self.radialBasisWeightsMarginal, self.verbosity >= 5, True, {"derIdxs": self._derMIdxs, "reorder": self._reorderM, "scl": np.power(self.scaleFactorMarginal, -1.)}, {"rcond": self.interpRcondMarginal}) vbMng(self, "MAIN", msg, 5) if wellCond: break RROMPyWarning(("Polyfit is poorly conditioned. Reducing " "MMarginal by 1.")) self.MMarginal -= 1 if self.MMarginal < 0: raise RROMPyException(("Instability in computation of " "marginal interpolator. Aborting.")) mI = mI + [copy(p)] vbMng(self, "DEL", "Done computing marginal interpolator.", 7) return mI def normApprox(self, mu:paramList, homogeneized : bool = False) -> float: """ Compute norm of approximant at arbitrary parameter. Args: mu: Target parameter. homogeneized(optional): Whether to remove Dirichlet BC. Defaults to False. Returns: Target norm of approximant. """ if not self.POD or self.homogeneized != homogeneized: return super().normApprox(mu, homogeneized) return np.linalg.norm(self.getApproxReduced(mu).data, axis = 0) def computeScaleFactor(self): """Compute parameter rescaling factor.""" RROMPyAssert(self._mode, message = "Cannot compute rescaling factor.") self.scaleFactorPivot = .5 * np.abs( self.muBoundsPivot[0] ** self.rescalingExpPivot - self.muBoundsPivot[1] ** self.rescalingExpPivot) self.scaleFactorMarginal = .5 * np.abs( self.muBoundsMarginal[0] ** self.rescalingExpMarginal - self.muBoundsMarginal[1] ** self.rescalingExpMarginal) self.scaleFactor = np.empty(self.npar) self.scaleFactor[self.directionPivot] = self.scaleFactorPivot self.scaleFactor[self.directionMarginal] = self.scaleFactorMarginal def centerNormalizeMarginal(self, mu : paramList = [], mu0 : paramVal = None) -> paramList: """ Compute normalized parameter to be plugged into approximant. Args: mu: Parameter(s) 1. mu0: Parameter(s) 2. If None, set to self.mu0. Returns: Normalized parameter. """ return self.trainedModel.centerNormalizeMarginal(mu, mu0) diff --git a/rrompy/reduction_methods/pivoting/rational_interpolant_pivoted.py b/rrompy/reduction_methods/pivoting/rational_interpolant_pivoted.py index 32fb25b..d12c423 100644 --- a/rrompy/reduction_methods/pivoting/rational_interpolant_pivoted.py +++ b/rrompy/reduction_methods/pivoting/rational_interpolant_pivoted.py @@ -1,606 +1,605 @@ # 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 . # from copy import deepcopy as copy import numpy as np from rrompy.reduction_methods.base import checkRobustTolerance from .generic_pivoted_approximant import GenericPivotedApproximant from rrompy.reduction_methods.distributed.rational_interpolant import ( RationalInterpolant as RI) -from rrompy.utilities.poly_fitting import customPInv from rrompy.utilities.poly_fitting.polynomial import (polybases, polyfitname, nextDerivativeIndices, hashDerivativeToIdx as hashD, hashIdxToDerivative as hashI, homogeneizedpolyvander as hpvP, homogeneizedToFull, PolynomialInterpolator as PI) from rrompy.utilities.poly_fitting.radial_basis import (rbbases, RadialBasisInterpolator as RBI) from rrompy.reduction_methods.trained_model import (TrainedModelPivotedData, TrainedModelPivotedRational as tModel) from rrompy.utilities.base.types import (Np1D, Np2D, HFEng, DictAny, List, ListAny, paramVal, paramList) -from rrompy.utilities.base import (verbosityManager as vbMng, multifactorial, - freepar as fp) +from rrompy.utilities.base import verbosityManager as vbMng, freepar as fp +from rrompy.utilities.numerical import multifactorial, customPInv from rrompy.utilities.exception_manager import (RROMPyException, RROMPyAssert, RROMPyWarning) from rrompy.parameter import checkParameter __all__ = ['RationalInterpolantPivoted'] class RationalInterpolantPivoted(GenericPivotedApproximant): """ ROM pivoted rational interpolant computation for parametric problems. Args: HFEngine: HF problem solver. mu0(optional): Default parameter. Defaults to 0. directionPivot(optional): Pivot components. Defaults to [0]. approxParameters(optional): Dictionary containing values for main parameters of approximant. Recognized keys are: - 'POD': whether to compute POD of snapshots; defaults to True; - 'S': total number of pivot samples current approximant relies upon; - 'samplerPivot': pivot sample point generator; - 'SMarginal': total number of marginal samples current approximant relies upon; - 'samplerMarginal': marginal sample point generator; - 'polybasisPivot': type of polynomial basis for pivot interpolation; allowed values include 'MONOMIAL', 'CHEBYSHEV' and 'LEGENDRE'; defaults to 'MONOMIAL'; - 'polybasisMarginal': type of polynomial basis for marginal interpolation; allowed values include 'MONOMIAL', 'CHEBYSHEV' and 'LEGENDRE'; defaults to 'MONOMIAL'; - 'E': number of derivatives used to compute interpolant; defaults to 0; - 'M': degree of rational interpolant numerator; defaults to 0; - 'N': degree of rational interpolant denominator; defaults to 0; - 'MMarginal': degree of marginal interpolant; defaults to 0; - 'radialBasisPivot': radial basis family for pivot numerator; defaults to 0, i.e. no radial basis; - 'radialBasisMarginal': radial basis family for marginal interpolant; defaults to 0, i.e. no radial basis; - 'radialBasisWeightsPivot': radial basis weights for pivot numerator; defaults to 0, i.e. identity; - 'radialBasisWeightsMarginal': radial basis weights for marginal interpolant; defaults to 0, i.e. identity; - 'interpRcondPivot': tolerance for pivot interpolation; defaults to None; - 'interpRcondMarginal': tolerance for marginal interpolation; defaults to None; - 'robustTol': tolerance for robust rational denominator management; defaults to 0. Defaults to empty dict. homogeneized(optional): Whether to homogeneize Dirichlet BCs. Defaults to False. verbosity(optional): Verbosity level. Defaults to 10. Attributes: HFEngine: HF problem solver. mu0: Default parameter. directionPivot: Pivot components. mus: Array of snapshot parameters. musPivot: Array of pivot snapshot parameters. musMarginal: Array of marginal snapshot parameters. homogeneized: Whether to homogeneize Dirichlet BCs. approxParameters: Dictionary containing values for main parameters of approximant. Recognized keys are in parameterList. parameterListSoft: Recognized keys of soft approximant parameters: - 'POD': whether to compute POD of snapshots; - 'polybasisPivot': type of polynomial basis for pivot interpolation; - 'polybasisMarginal': type of polynomial basis for marginal interpolation; - 'E': number of derivatives used to compute interpolant; - 'M': degree of rational interpolant numerator; - 'N': degree of rational interpolant denominator; - 'MMarginal': degree of marginal interpolant; - 'radialBasisPivot': radial basis family for pivot numerator; - 'radialBasisMarginal': radial basis family for marginal interpolant; - 'radialBasisWeightsPivot': radial basis weights for pivot numerator; - 'radialBasisWeightsMarginal': radial basis weights for marginal interpolant; - 'interpRcondPivot': tolerance for pivot interpolation; - 'interpRcondMarginal': tolerance for marginal interpolation; - 'robustTol': tolerance for robust rational denominator management. parameterListCritical: Recognized keys of critical approximant parameters: - 'S': total number of pivot samples current approximant relies upon; - 'samplerPivot': pivot sample point generator; - 'SMarginal': total number of marginal samples current approximant relies upon; - 'samplerMarginal': marginal sample point generator. POD: Whether to compute POD of snapshots. S: Total number of pivot samples current approximant relies upon. sampler: Pivot sample point generator. polybasisPivot: Type of polynomial basis for pivot interpolation. polybasisMarginal: Type of polynomial basis for marginal interpolation. M: Numerator degree of approximant. N: Denominator degree of approximant. MMarginal: Degree of marginal interpolant. radialBasisPivot: Radial basis family for pivot numerator. radialBasisMarginal: Radial basis family for marginal interpolant. radialBasisWeightsPivot: Radial basis weights for pivot numerator. radialBasisWeightsMarginal: Radial basis weights for marginal interpolant. interpRcondPivot: Tolerance for pivot interpolation. interpRcondMarginal: Tolerance for marginal interpolation. robustTol: Tolerance for robust rational denominator management. E: Complete derivative depth level of S. muBoundsPivot: list of bounds for pivot parameter values. muBoundsMarginal: list of bounds for marginal parameter values. samplingEngine: Sampling engine. uHF: High fidelity solution(s) with parameter(s) lastSolvedHF as sampleList. lastSolvedHF: Parameter(s) corresponding to last computed high fidelity solution(s) as parameterList. uApproxReduced: Reduced approximate solution(s) with parameter(s) lastSolvedApprox as sampleList. lastSolvedApproxReduced: Parameter(s) corresponding to last computed reduced approximate solution(s) as parameterList. uApprox: Approximate solution(s) with parameter(s) lastSolvedApprox as sampleList. lastSolvedApprox: Parameter(s) corresponding to last computed approximate solution(s) as parameterList. Q: Numpy 1D vector containing complex coefficients of approximant denominator. P: Numpy 2D vector whose columns are FE dofs of coefficients of approximant numerator. """ def __init__(self, HFEngine:HFEng, mu0 : paramVal = None, directionPivot : ListAny = [0], approxParameters : DictAny = {}, homogeneized : bool = False, verbosity : int = 10, timestamp : bool = True): self._preInit() self._addParametersToList( ["polybasisPivot", "E", "M", "N", "radialBasisPivot", "radialBasisWeightsPivot", "interpRcondPivot", "robustTol"], ["MONOMIAL", -1, 0, 0, 0, 1, -1, 0]) super().__init__(HFEngine = HFEngine, mu0 = mu0, directionPivot = directionPivot, approxParameters = approxParameters, homogeneized = homogeneized, verbosity = verbosity, timestamp = timestamp) self._postInit() @property def polybasisPivot(self): """Value of polybasisPivot.""" return self._polybasisPivot @polybasisPivot.setter def polybasisPivot(self, polybasisPivot): try: polybasisPivot = polybasisPivot.upper().strip().replace(" ","") if polybasisPivot not in polybases: raise RROMPyException( "Prescribed pivot polybasis not recognized.") self._polybasisPivot = polybasisPivot except: RROMPyWarning(("Prescribed pivot polybasis not recognized. " "Overriding to 'MONOMIAL'.")) self._polybasisPivot = "MONOMIAL" self._approxParameters["polybasisPivot"] = self.polybasisPivot @property def radialBasisPivot(self): """Value of radialBasisPivot.""" return self._radialBasisPivot @radialBasisPivot.setter def radialBasisPivot(self, radialBasisPivot): try: if radialBasisPivot != 0: radialBasisPivot = radialBasisPivot.upper().strip().replace( " ","") if radialBasisPivot not in rbbases: raise RROMPyException(("Prescribed pivot radialBasis not " "recognized.")) self._radialBasisPivot = radialBasisPivot except: RROMPyWarning(("Prescribed pivot radialBasis not recognized. " "Overriding to 0.")) self._radialBasisPivot = 0 self._approxParameters["radialBasisPivot"] = self.radialBasisPivot @property def radialBasisWeightsPivot(self): """Value of radialBasisWeightsPivot.""" return self._radialBasisWeightsPivot @radialBasisWeightsPivot.setter def radialBasisWeightsPivot(self, radialBasisWeightsPivot): self._radialBasisWeightsPivot = radialBasisWeightsPivot self._approxParameters["radialBasisWeightsPivot"] = ( self.radialBasisWeightsPivot) @property def polybasisPivotP(self): if self.radialBasisPivot == 0: return self._polybasisPivot return self._polybasisPivot + "_" + self.radialBasisPivot @property def interpRcondPivot(self): """Value of interpRcondPivot.""" return self._interpRcondPivot @interpRcondPivot.setter def interpRcondPivot(self, interpRcondPivot): self._interpRcondPivot = interpRcondPivot self._approxParameters["interpRcondPivot"] = self.interpRcondPivot @property def M(self): """Value of M. Its assignment may change S.""" return self._M @M.setter def M(self, M): if M < 0: raise RROMPyException("M must be non-negative.") self._M = M self._approxParameters["M"] = self.M if hasattr(self, "_E") and self.E >= 0 and self.E < self.M: RROMPyWarning("Prescribed S is too small. Decreasing M.") self.M = self.E @property def N(self): """Value of N. Its assignment may change S.""" return self._N @N.setter def N(self, N): if N < 0: raise RROMPyException("N must be non-negative.") self._N = N self._approxParameters["N"] = self.N if hasattr(self, "_E") and self.E >= 0 and self.E < self.N: RROMPyWarning("Prescribed S is too small. Decreasing N.") self.N = self.E @property def E(self): """Value of E.""" return self._E @E.setter def E(self, E): if E < 0: if not hasattr(self, "_S"): raise RROMPyException(("Value of E must be positive if S is " "not yet assigned.")) E = np.sum(hashI(np.prod(self.S), len(self.directionPivot))) - 1 self._E = E self._approxParameters["E"] = self.E if (hasattr(self, "_S") and self.E >= np.sum(hashI(np.prod(self.S),len(self.directionPivot)))): RROMPyWarning("Prescribed S is too small. Decreasing E.") self.E = -1 if hasattr(self, "_M"): self.M = self.M if hasattr(self, "_N"): self.N = self.N @property def robustTol(self): """Value of tolerance for robust rational denominator management.""" return self._robustTol @robustTol.setter def robustTol(self, robustTol): if robustTol < 0.: RROMPyWarning(("Overriding prescribed negative robustness " "tolerance to 0.")) robustTol = 0. self._robustTol = robustTol self._approxParameters["robustTol"] = self.robustTol @property def S(self): """Value of S.""" return self._S @S.setter def S(self, S): GenericPivotedApproximant.S.fset(self, S) if hasattr(self, "_M"): self.M = self.M if hasattr(self, "_N"): self.N = self.N if hasattr(self, "_E"): self.E = self.E def resetSamples(self): """Reset samples.""" super().resetSamples() self._musPUniqueCN = None self._derPIdxs = None self._reorderP = None def _setupPivotInterpolationIndices(self): """Setup parameters for polyvander.""" RROMPyAssert(self._mode, message = "Cannot setup interpolation indices.") if (self._musPUniqueCN is None or len(self._reorderP) != len(self.musPivot)): self._musPUniqueCN, musPIdxsTo, musPIdxs, musPCount = ( self.centerNormalizePivot(self.musPivot).unique( return_index = True, return_inverse = True, return_counts = True)) self._musPUnique = self.mus[musPIdxsTo] self._derPIdxs = [None] * len(self._musPUniqueCN) self._reorderP = np.empty(len(musPIdxs), dtype = int) filled = 0 for j, cnt in enumerate(musPCount): self._derPIdxs[j] = nextDerivativeIndices([], self.musPivot.shape[1], cnt) jIdx = np.nonzero(musPIdxs == j)[0] self._reorderP[jIdx] = np.arange(filled, filled + cnt) filled += cnt def _setupDenominator(self): """Compute rational denominator.""" RROMPyAssert(self._mode, message = "Cannot setup denominator.") vbMng(self, "INIT", "Starting computation of denominator.", 7) NinvD = None N0 = copy(self.N) qs = [] self.verbosity -= 10 for j in range(len(self.musMarginal)): self._N = N0 while self.N > 0: if NinvD != self.N: invD = self._computeInterpolantInverseBlocks() NinvD = self.N if self.POD: ev, eV = RI.findeveVGQR(self, self.samplingEngine.RPOD[j], invD) else: ev, eV = RI.findeveVGExplicit(self, self.samplingEngine.samples[j], invD) newParams = checkRobustTolerance(ev, self.N, self.robustTol) if not newParams: break self.approxParameters = newParams if self.N <= 0: self._N = 0 eV = np.ones((1, 1)) q = PI() q.npar = self.musPivot.shape[1] q.polybasis = self.polybasisPivot q.coeffs = homogeneizedToFull(tuple([self.N + 1] * q.npar), q.npar, eV[:, 0]) qs = qs + [copy(q)] self.verbosity += 10 vbMng(self, "DEL", "Done computing denominator.", 7) return qs def _setupNumerator(self): """Compute rational numerator.""" RROMPyAssert(self._mode, message = "Cannot setup numerator.") vbMng(self, "INIT", "Starting computation of numerator.", 7) Qevaldiag = np.zeros((len(self.musPivot), len(self.musPivot)), dtype = np.complex) verb = self.trainedModel.verbosity self.trainedModel.verbosity = 0 self._setupPivotInterpolationIndices() M0 = copy(self.M) tensor_idx = 0 ps = [] for k, muM in enumerate(self.musMarginal): self._M = M0 idxGlob = 0 for j, derIdxs in enumerate(self._derPIdxs): mujEff = [fp] * self.npar for jj, kk in enumerate(self.directionPivot): mujEff[kk] = self._musPUnique[j, jj] for jj, kk in enumerate(self.directionMarginal): mujEff[kk] = muM(0, jj) mujEff = checkParameter(mujEff, self.npar) nder = len(derIdxs) idxLoc = np.arange(len(self.musPivot))[ (self._reorderP >= idxGlob) * (self._reorderP < idxGlob + nder)] idxGlob += nder Qval = [0] * nder for der in range(nder): derIdx = hashI(der, self.musPivot.shape[1]) derIdxEff = [0] * self.npar sclEff = [0] * self.npar for jj, kk in enumerate(self.directionPivot): derIdxEff[kk] = derIdx[jj] sclEff[kk] = self.scaleFactorPivot[jj] ** -1. Qval[der] = (self.trainedModel.getQVal(mujEff, derIdxEff, scl = sclEff) / multifactorial(derIdx)) for derU, derUIdx in enumerate(derIdxs): for derQ, derQIdx in enumerate(derIdxs): diffIdx = [x - y for (x, y) in zip(derUIdx, derQIdx)] if all([x >= 0 for x in diffIdx]): diffj = hashD(diffIdx) Qevaldiag[idxLoc[derU], idxLoc[derQ]] = Qval[diffj] while self.M >= 0: if self.radialBasisPivot == 0: p = PI() wellCond, msg = p.setupByInterpolation( self._musPUniqueCN, Qevaldiag, self.M, self.polybasisPivotP, self.verbosity >= 5, True, {"derIdxs": self._derPIdxs, "reorder": self._reorderP, "scl": np.power(self.scaleFactorPivot, -1.)}, {"rcond": self.interpRcondPivot}) else: p = RBI() wellCond, msg = p.setupByInterpolation( self._musPUniqueCN, Qevaldiag, self.M, self.polybasisPivotP, self.radialBasisWeightsPivot, self.verbosity >= 5, True, {"derIdxs": self._derPIdxs, "reorder": self._reorderP, "scl": np.power(self.scaleFactorPivot, -1.)}, {"rcond": self.interpRcondPivot}) vbMng(self, "MAIN", msg, 5) if wellCond: break RROMPyWarning(("Polyfit is poorly conditioned. " "Reducing M by 1.")) self.M -= 1 if self.M < 0: raise RROMPyException(("Instability in computation of " "numerator. Aborting.")) tensor_idx_new = tensor_idx + Qevaldiag.shape[1] if self.POD: p.postmultiplyTensorize(self.samplingEngine.RPODCoalesced.T[ tensor_idx : tensor_idx_new, :]) else: p.pad(tensor_idx, len(self.mus) - tensor_idx_new) tensor_idx = tensor_idx_new ps = ps + [copy(p)] self.trainedModel.verbosity = verb vbMng(self, "DEL", "Done computing numerator.", 7) return ps def setupApprox(self): """ Compute rational interpolant. SVD-based robust eigenvalue management. """ if self.checkComputedApprox(): return RROMPyAssert(self._mode, message = "Cannot setup approximant.") vbMng(self, "INIT", "Setting up {}.". format(self.name()), 5) self.computeScaleFactor() self.computeSnapshots() if self.trainedModel is None: self.trainedModel = tModel() self.trainedModel.verbosity = self.verbosity self.trainedModel.timestamp = self.timestamp data = TrainedModelPivotedData(self.trainedModel.name(), self.mu0, self.samplingEngine.samplesCoalesced, self.HFEngine.rescalingExp, self.directionPivot, self.scaleFactor) data.musPivot = copy(self.musPivot) data.musMarginal = copy(self.musMarginal) self.trainedModel.data = data else: self.trainedModel = self.trainedModel self.trainedModel.data.projMat = copy(self.samplingEngine.samples) self.trainedModel.data.marginalInterp = self._setupMarginalInterp() if self.N > 0: Qs = self._setupDenominator() else: Q = PI() Q.npar = self.musPivot.shape[1] Q.coeffs = np.ones(tuple([1] * Q.npar), dtype = np.complex) Q.polybasis = self.polybasisPivot Qs = [Q for _ in range(len(self.musMarginal))] self.trainedModel.data.Qs = Qs self.trainedModel.data.Ps = self._setupNumerator() self.trainedModel.data.approxParameters = copy(self.approxParameters) vbMng(self, "DEL", "Done setting up approximant.", 5) def _computeInterpolantInverseBlocks(self) -> List[Np2D]: """ Compute inverse factors for minimal interpolant target functional. """ RROMPyAssert(self._mode, message = "Cannot solve eigenvalue problem.") self._setupPivotInterpolationIndices() while self.E >= 0: eWidth = (hashD([self.E + 1] + [0] * (self.musPivot.shape[1] - 1)) - hashD([self.E] + [0] * (self.musPivot.shape[1] - 1))) TE, _, argIdxs = hpvP(self._musPUniqueCN, self.E, self.polybasisPivot, self._derPIdxs, self._reorderP, scl = np.power(self.scaleFactorPivot, -1.)) fitOut = customPInv(TE[:, argIdxs], rcond = self.interpRcondPivot, full = True) vbMng(self, "MAIN", ("Fitting {} samples with degree {} through {}... " "Conditioning of pseudoinverse system: {:.4e}.").format( TE.shape[0], self.E, polyfitname(self.polybasisPivot), fitOut[1][1][0] / fitOut[1][1][-1]), 5) if fitOut[1][0] == len(argIdxs): self._fitinvP = fitOut[0][- eWidth : , :] break RROMPyWarning("Polyfit is poorly conditioned. Reducing E by 1.") self.E -= 1 if self.E < 0: raise RROMPyException(("Instability in computation of " "denominator. Aborting.")) TN, _, argIdxs = hpvP(self._musPUniqueCN, self.N, self.polybasisPivot, self._derPIdxs, self._reorderP, scl = np.power(self.scaleFactorPivot, -1.)) TN = TN[:, argIdxs] invD = [None] * (eWidth) for k in range(eWidth): pseudoInv = np.diag(self._fitinvP[k, :]) idxGlob = 0 for j, derIdxs in enumerate(self._derPIdxs): nder = len(derIdxs) idxGlob += nder if nder > 1: idxLoc = np.arange(len(self.musPivot))[ (self._reorderP >= idxGlob - nder) * (self._reorderP < idxGlob)] invLoc = self._fitinvP[k, idxLoc] pseudoInv[np.ix_(idxLoc, idxLoc)] = 0. for diffj, diffjIdx in enumerate(derIdxs): for derQ, derQIdx in enumerate(derIdxs): derUIdx = [x - y for (x, y) in zip(diffjIdx, derQIdx)] if all([x >= 0 for x in derUIdx]): derU = hashD(derUIdx) pseudoInv[idxLoc[derU], idxLoc[derQ]] = ( invLoc[diffj]) invD[k] = pseudoInv.dot(TN) return invD def centerNormalize(self, mu : paramList = [], mu0 : paramVal = None) -> paramList: """ Compute normalized parameter to be plugged into approximant. Args: mu: Parameter(s) 1. mu0: Parameter(s) 2. If None, set to self.mu0. Returns: Normalized parameter. """ return self.trainedModel.centerNormalize(mu, mu0) def centerNormalizePivot(self, mu : paramList = [], mu0 : paramVal = None) -> paramList: """ Compute normalized parameter to be plugged into approximant. Args: mu: Parameter(s) 1. mu0: Parameter(s) 2. If None, set to self.mu0. Returns: Normalized parameter. """ return self.trainedModel.centerNormalizePivot(mu, mu0) def getResidues(self, *args, **kwargs) -> Np1D: """ Obtain approximant residues. Returns: Matrix with residues as columns. """ return self.trainedModel.getResidues(*args, **kwargs) diff --git a/rrompy/reduction_methods/pivoting/rational_interpolant_pivoted_pole_matching.py b/rrompy/reduction_methods/pivoting/rational_interpolant_pivoted_pole_matching.py new file mode 100644 index 0000000..5c465f6 --- /dev/null +++ b/rrompy/reduction_methods/pivoting/rational_interpolant_pivoted_pole_matching.py @@ -0,0 +1,231 @@ +# 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 . +# + +from copy import deepcopy as copy +import numpy as np +from .rational_interpolant_pivoted import RationalInterpolantPivoted +from rrompy.utilities.poly_fitting.polynomial import ( + PolynomialInterpolator as PI) +from rrompy.reduction_methods.trained_model import (TrainedModelPivotedData, + TrainedModelPivotedRationalPoleMatching as tModel) +from rrompy.utilities.base.types import HFEng, DictAny, ListAny, paramVal +from rrompy.utilities.base import verbosityManager as vbMng +from rrompy.utilities.exception_manager import (RROMPyException, RROMPyAssert) + +__all__ = ['RationalInterpolantPivotedPoleMatching'] + +class RationalInterpolantPivotedPoleMatching(RationalInterpolantPivoted): + """ + ROM pivoted rational interpolant computation for parametric problems with + pole matching. + + Args: + HFEngine: HF problem solver. + mu0(optional): Default parameter. Defaults to 0. + directionPivot(optional): Pivot components. Defaults to [0]. + approxParameters(optional): Dictionary containing values for main + parameters of approximant. Recognized keys are: + - 'POD': whether to compute POD of snapshots; defaults to True; + - 'matchingWeight': weight for pole matching optimization; defaults + to 1; + - 'S': total number of pivot samples current approximant relies + upon; + - 'samplerPivot': pivot sample point generator; + - 'SMarginal': total number of marginal samples current approximant + relies upon; + - 'samplerMarginal': marginal sample point generator; + - 'polybasisPivot': type of polynomial basis for pivot + interpolation; allowed values include 'MONOMIAL', 'CHEBYSHEV' + and 'LEGENDRE'; defaults to 'MONOMIAL'; + - 'polybasisMarginal': type of polynomial basis for marginal + interpolation; allowed values include 'MONOMIAL', 'CHEBYSHEV' + and 'LEGENDRE'; defaults to 'MONOMIAL'; + - 'E': number of derivatives used to compute interpolant; defaults + to 0; + - 'M': degree of rational interpolant numerator; defaults to 0; + - 'N': degree of rational interpolant denominator; defaults to 0; + - 'MMarginal': degree of marginal interpolant; defaults to 0; + - 'radialBasisPivot': radial basis family for pivot numerator; + defaults to 0, i.e. no radial basis; + - 'radialBasisMarginal': radial basis family for marginal + interpolant; defaults to 0, i.e. no radial basis; + - 'radialBasisWeightsPivot': radial basis weights for pivot + numerator; defaults to 0, i.e. identity; + - 'radialBasisWeightsMarginal': radial basis weights for marginal + interpolant; defaults to 0, i.e. identity; + - 'interpRcondPivot': tolerance for pivot interpolation; defaults + to None; + - 'interpRcondMarginal': tolerance for marginal interpolation; + defaults to None; + - 'robustTol': tolerance for robust rational denominator + management; defaults to 0. + Defaults to empty dict. + homogeneized(optional): Whether to homogeneize Dirichlet BCs. Defaults + to False. + verbosity(optional): Verbosity level. Defaults to 10. + + Attributes: + HFEngine: HF problem solver. + mu0: Default parameter. + directionPivot: Pivot components. + mus: Array of snapshot parameters. + musPivot: Array of pivot snapshot parameters. + musMarginal: Array of marginal snapshot parameters. + homogeneized: Whether to homogeneize Dirichlet BCs. + approxParameters: Dictionary containing values for main parameters of + approximant. Recognized keys are in parameterList. + parameterListSoft: Recognized keys of soft approximant parameters: + - 'POD': whether to compute POD of snapshots; + - 'matchingWeight': weight for pole matching optimization; + - 'polybasisPivot': type of polynomial basis for pivot + interpolation; + - 'polybasisMarginal': type of polynomial basis for marginal + interpolation; + - 'E': number of derivatives used to compute interpolant; + - 'M': degree of rational interpolant numerator; + - 'N': degree of rational interpolant denominator; + - 'MMarginal': degree of marginal interpolant; + - 'radialBasisPivot': radial basis family for pivot numerator; + - 'radialBasisMarginal': radial basis family for marginal + interpolant; + - 'radialBasisWeightsPivot': radial basis weights for pivot + numerator; + - 'radialBasisWeightsMarginal': radial basis weights for marginal + interpolant; + - 'interpRcondPivot': tolerance for pivot interpolation; + - 'interpRcondMarginal': tolerance for marginal interpolation; + - 'robustTol': tolerance for robust rational denominator + management. + parameterListCritical: Recognized keys of critical approximant + parameters: + - 'S': total number of pivot samples current approximant relies + upon; + - 'samplerPivot': pivot sample point generator; + - 'SMarginal': total number of marginal samples current approximant + relies upon; + - 'samplerMarginal': marginal sample point generator. + POD: Whether to compute POD of snapshots. + matchingWeight: Weight for pole matching optimization. + S: Total number of pivot samples current approximant relies upon. + sampler: Pivot sample point generator. + polybasisPivot: Type of polynomial basis for pivot interpolation. + polybasisMarginal: Type of polynomial basis for marginal interpolation. + M: Numerator degree of approximant. + N: Denominator degree of approximant. + MMarginal: Degree of marginal interpolant. + radialBasisPivot: Radial basis family for pivot numerator. + radialBasisMarginal: Radial basis family for marginal interpolant. + radialBasisWeightsPivot: Radial basis weights for pivot numerator. + radialBasisWeightsMarginal: Radial basis weights for marginal + interpolant. + interpRcondPivot: Tolerance for pivot interpolation. + interpRcondMarginal: Tolerance for marginal interpolation. + robustTol: Tolerance for robust rational denominator management. + E: Complete derivative depth level of S. + muBoundsPivot: list of bounds for pivot parameter values. + muBoundsMarginal: list of bounds for marginal parameter values. + samplingEngine: Sampling engine. + uHF: High fidelity solution(s) with parameter(s) lastSolvedHF as + sampleList. + lastSolvedHF: Parameter(s) corresponding to last computed high fidelity + solution(s) as parameterList. + uApproxReduced: Reduced approximate solution(s) with parameter(s) + lastSolvedApprox as sampleList. + lastSolvedApproxReduced: Parameter(s) corresponding to last computed + reduced approximate solution(s) as parameterList. + uApprox: Approximate solution(s) with parameter(s) lastSolvedApprox as + sampleList. + lastSolvedApprox: Parameter(s) corresponding to last computed + approximate solution(s) as parameterList. + Q: Numpy 1D vector containing complex coefficients of approximant + denominator. + P: Numpy 2D vector whose columns are FE dofs of coefficients of + approximant numerator. + """ + + def __init__(self, HFEngine:HFEng, mu0 : paramVal = None, + directionPivot : ListAny = [0], + approxParameters : DictAny = {}, homogeneized : bool = False, + verbosity : int = 10, timestamp : bool = True): + self._preInit() + if len(directionPivot) > 1: + raise RROMPyException(("Exactly 1 pivot parameter allowed in pole " + "matching.")) + self._addParametersToList(["matchingWeight"], [1]) + super().__init__(HFEngine = HFEngine, mu0 = mu0, + directionPivot = directionPivot, + approxParameters = approxParameters, + homogeneized = homogeneized, + verbosity = verbosity, timestamp = timestamp) + self._postInit() + + @property + def matchingWeight(self): + """Value of matchingWeight.""" + return self._matchingWeight + @matchingWeight.setter + def matchingWeight(self, matchingWeight): + self._matchingWeight = matchingWeight + self._approxParameters["matchingWeight"] = self.matchingWeight + + def setupApprox(self): + """ + Compute rational interpolant. + SVD-based robust eigenvalue management. + """ + if self.checkComputedApprox(): + return + RROMPyAssert(self._mode, message = "Cannot setup approximant.") + vbMng(self, "INIT", "Setting up {}.". format(self.name()), 5) + self.computeScaleFactor() + self.computeSnapshots() + if self.trainedModel is None: + self.trainedModel = tModel() + self.trainedModel.verbosity = self.verbosity + self.trainedModel.timestamp = self.timestamp + data = TrainedModelPivotedData(self.trainedModel.name(), self.mu0, + self.samplingEngine.samplesCoalesced, + self.HFEngine.rescalingExp, + self.directionPivot, + self.scaleFactor) + data.musPivot = copy(self.musPivot) + data.musMarginal = copy(self.musMarginal) + self.trainedModel.data = data + else: + self.trainedModel = self.trainedModel + self.trainedModel.data.projMat = copy( + self.samplingEngine.samplesCoalesced) + self.trainedModel.data.marginalInterp = self._setupMarginalInterp() + if self.N > 0: + Qs = self._setupDenominator() + else: + Q = PI() + Q.npar = self.musPivot.shape[1] + Q.coeffs = np.ones(tuple([1] * Q.npar), dtype = np.complex) + Q.polybasis = self.polybasisPivot + Qs = [Q for _ in range(len(self.musMarginal))] + self.trainedModel.data._temporary = True + self.trainedModel.data.Qs = Qs + self.trainedModel.data.Ps = self._setupNumerator() + vbMng(self, "INIT", "Matching poles.", 10) + self.trainedModel.initializeFromRational(self.HFEngine, + self.matchingWeight, self.POD) + vbMng(self, "DEL", "Done matching poles.", 10) + self.trainedModel.data.approxParameters = copy(self.approxParameters) + vbMng(self, "DEL", "Done setting up approximant.", 5) + \ No newline at end of file diff --git a/rrompy/reduction_methods/pivoting/rb_distributed_pivoted.py b/rrompy/reduction_methods/pivoting/rb_distributed_pivoted.py index af9c165..f18460b 100644 --- a/rrompy/reduction_methods/pivoting/rb_distributed_pivoted.py +++ b/rrompy/reduction_methods/pivoting/rb_distributed_pivoted.py @@ -1,290 +1,282 @@ # 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 . # from copy import deepcopy as copy import numpy as np from .generic_pivoted_approximant import GenericPivotedApproximant from rrompy.hfengines.base import MarginalProxyEngine from rrompy.reduction_methods.trained_model import (TrainedModelPivotedData, TrainedModelPivotedRB as tModel) from rrompy.reduction_methods.base.rb_utils import projectAffineDecomposition from rrompy.utilities.base.types import (Np1D, Np2D, List, ListAny, Tuple, DictAny, HFEng, paramVal, sampList) from rrompy.utilities.base import verbosityManager as vbMng, freepar as fp -from rrompy.utilities.poly_fitting import customPInv +from rrompy.utilities.numerical import customPInv from rrompy.utilities.exception_manager import (RROMPyWarning, RROMPyException, RROMPyAssert) __all__ = ['RBDistributedPivoted'] class RBDistributedPivoted(GenericPivotedApproximant): """ ROM pivoted RB approximant computation for parametric problems. Args: HFEngine: HF problem solver. mu0(optional): Default parameter. Defaults to 0. directionPivot(optional): Pivot components. Defaults to [0]. approxParameters(optional): Dictionary containing values for main parameters of approximant. Recognized keys are: - - 'POD': whether to compute POD of snapshots; defaults to True; - 'S': total number of pivot samples current approximant relies upon; - 'samplerPivot': pivot sample point generator; - 'SMarginal': total number of marginal samples current approximant relies upon; - 'samplerMarginal': marginal sample point generator; - 'R': rank for pivot Galerkin projection; defaults to prod(S); - 'PODTolerance': tolerance for pivot snapshots POD; defaults to -1; - 'polybasisMarginal': type of polynomial basis for marginal interpolation; allowed values include 'MONOMIAL', 'CHEBYSHEV' and 'LEGENDRE'; defaults to 'MONOMIAL'; - 'MMarginal': degree of marginal interpolant; defaults to 0; - 'radialBasisMarginal': radial basis family for marginal interpolant; defaults to 0, i.e. no radial basis; - 'radialBasisWeightsMarginal': radial basis weights for marginal interpolant; defaults to 0, i.e. identity; - 'interpRcondMarginal': tolerance for marginal interpolation; defaults to None. Defaults to empty dict. homogeneized(optional): Whether to homogeneize Dirichlet BCs. Defaults to False. verbosity(optional): Verbosity level. Defaults to 10. Attributes: HFEngine: HF problem solver. mu0: Default parameter. directionPivot: Pivot components. mus: Array of snapshot parameters. homogeneized: Whether to homogeneize Dirichlet BCs. approxRadius: Dummy radius of approximant (i.e. distance from mu0 to farthest sample point). approxParameters: Dictionary containing values for main parameters of approximant. Recognized keys are in parameterList. parameterListSoft: Recognized keys of soft approximant parameters: - 'POD': whether to compute POD of snapshots; - 'R': rank for Galerkin projection; - 'PODTolerance': tolerance for snapshots POD; - 'polybasisMarginal': type of polynomial basis for marginal interpolation; - 'MMarginal': degree of marginal interpolant; - 'radialBasisMarginal': radial basis family for marginal interpolant; - 'radialBasisWeightsMarginal': radial basis weights for marginal interpolant; - 'interpRcondMarginal': tolerance for marginal interpolation. parameterListCritical: Recognized keys of critical approximant parameters: - 'S': total number of pivot samples current approximant relies upon; - 'samplerPivot': pivot sample point generator; - 'SMarginal': total number of marginal samples current approximant relies upon; - 'samplerMarginal': marginal sample point generator. POD: Whether to compute POD of snapshots. S: Total number of pivot samples current approximant relies upon. samplerPivot: Pivot sample point generator. SMarginal: Total number of marginal samples current approximant relies upon. samplerMarginal: Marginal sample point generator. R: Rank for Galerkin projection. PODTolerance: Tolerance for pivot snapshots POD. polybasisMarginal: Type of polynomial basis for marginal interpolation. MMarginal: Degree of marginal interpolant. radialBasisMarginal: Radial basis family for marginal interpolant. radialBasisWeightsMarginal: Radial basis weights for marginal interpolant. interpRcondMarginal: Tolerance for marginal interpolation. muBoundsPivot: list of bounds for pivot parameter values. muBoundsMarginal: list of bounds for marginal parameter values. samplingEngine: Sampling engine. uHF: High fidelity solution(s) with parameter(s) lastSolvedHF as sampleList. lastSolvedHF: Parameter(s) corresponding to last computed high fidelity solution(s) as parameterList. uApproxReduced: Reduced approximate solution(s) with parameter(s) lastSolvedApprox as sampleList. lastSolvedApproxReduced: Parameter(s) corresponding to last computed reduced approximate solution(s) as parameterList. uApprox: Approximate solution(s) with parameter(s) lastSolvedApprox as sampleList. lastSolvedApprox: Parameter(s) corresponding to last computed approximate solution(s) as parameterList. As: List of sparse matrices (in CSC format) representing coefficients - of linear system matrix wrt theta(mu). + of linear system matrix. bs: List of numpy vectors representing coefficients of linear system - RHS wrt theta(mu). - thetaAs: List of callables representing coefficients of linear system - matrix wrt mu. - thetabs: List of callables representing coefficients of linear system - RHS wrt mu. + RHS. ARBs: List of sparse matrices (in CSC format) representing coefficients - of compressed linear system matrix wrt theta(mu). + of compressed linear system matrix. bRBs: List of numpy vectors representing coefficients of compressed - linear system RHS wrt theta(mu). + linear system RHS. """ def __init__(self, HFEngine:HFEng, mu0 : paramVal = None, directionPivot : ListAny = [0], approxParameters : DictAny = {}, homogeneized : bool = False, verbosity : int = 10, timestamp : bool = True): self._preInit() self._addParametersToList(["R", "PODTolerance"], [1, -1]) super().__init__(HFEngine = HFEngine, mu0 = mu0, directionPivot = directionPivot, approxParameters = approxParameters, homogeneized = homogeneized, verbosity = verbosity, timestamp = timestamp) self._postInit() @property def S(self): """Value of S.""" return self._S @S.setter def S(self, S): GenericPivotedApproximant.S.fset(self, S) if not hasattr(self, "_R"): self._R = np.prod(self.S) * np.prod(self.SMarginal) @property def R(self): """Value of R. Its assignment may change S.""" return self._R @R.setter def R(self, R): if R < 0: raise RROMPyException("R must be non-negative.") self._R = R self._approxParameters["R"] = self.R if (hasattr(self, "_S") and hasattr(self, "_SMarginal") and np.prod(self.S) * np.prod(self.SMarginal) < self.R): RROMPyWarning(("Prescribed S and SMarginal are too small. " "Decreasing R.")) self.R = np.prod(self.S) * np.prod(self.SMarginal) @property def PODTolerance(self): """Value of PODTolerance.""" return self._PODTolerance @PODTolerance.setter def PODTolerance(self, PODTolerance): self._PODTolerance = PODTolerance self._approxParameters["PODTolerance"] = self.PODTolerance def _setupProjectionMatrix(self): """Compute projection matrix.""" RROMPyAssert(self._mode, message = "Cannot setup numerator.") vbMng(self, "INIT", "Starting computation of projection matrix.", 7) if self.POD: U, s, _ = np.linalg.svd(self.samplingEngine.RPODCoalesced) s = s ** 2. else: Gramian = self.HFEngine.innerProduct( self.samplingEngine.samplesCoalesced, self.samplingEngine.samplesCoalesced) U, s, _ = np.linalg.svd(Gramian) nsamples = self.samplingEngine.samplesCoalesced.shape[1] snorm = np.cumsum(s[::-1]) / np.sum(s) nPODTrunc = min(nsamples - np.argmax(snorm > self.PODTolerance), self.R) pMat = self.samplingEngine.samplesCoalesced.dot(U[:, : nPODTrunc]) vbMng(self, "MAIN", "Assembling {}x{} projection matrix from {} samples.".format( *(pMat.shape), nsamples), 5) vbMng(self, "DEL", "Done computing projection matrix.", 7) return pMat def _setupAffineBlocks(self): """Compute list of marginalized affine blocks of system.""" hasAs = hasattr(self, "AsList") and self.AsList is not None hasbs = hasattr(self, "bsList") and self.bsList is not None if hasAs and hasbs: return vbMng(self, "INIT", "Computing affine blocks of system.", 10) mu0Eff = self.mu0.data[0, self.directionPivot] if not hasAs: self.AsList = [None] * len(self.musMarginal) if not hasbs: self.bsList = [None] * len(self.musMarginal) for k, muM in enumerate(self.musMarginal): muEff = [fp] * self.npar for jj, kk in enumerate(self.directionMarginal): muEff[kk] = muM(0, jj) MEnginek = MarginalProxyEngine(self.HFEngine, muEff) if not hasAs: self.AsList[k] = MEnginek.affineLinearSystemA(mu0Eff) if not hasbs: self.bsList[k] = MEnginek.affineLinearSystemb(mu0Eff, self.homogeneized) - self.thetaAs = MEnginek.affineWeightsA(mu0Eff) - self.thetabs = MEnginek.affineWeightsb(mu0Eff, self.homogeneized) vbMng(self, "DEL", "Done computing affine blocks.", 10) def setupApprox(self): """Compute RB projection matrix.""" if self.checkComputedApprox(): return RROMPyAssert(self._mode, message = "Cannot setup approximant.") vbMng(self, "INIT", "Setting up {}.".format(self.name()), 5) self.computeScaleFactor() self.computeSnapshots() self._setupAffineBlocks() pMat = self._setupProjectionMatrix() ARBsList, bRBsList, pList = self.assembleReducedSystemMarginal(pMat) if self.trainedModel is None: self.trainedModel = tModel() self.trainedModel.verbosity = self.verbosity self.trainedModel.timestamp = self.timestamp data = TrainedModelPivotedData(self.trainedModel.name(), self.mu0, pList, self.HFEngine.rescalingExp, self.directionPivot, self.scaleFactor) data.musPivot = copy(self.musPivot) data.musMarginal = copy(self.musMarginal) - data.thetaAs, data.thetabs = self.thetaAs, self.thetabs self.trainedModel.data = data else: self.trainedModel = self.trainedModel self.trainedModel.data.projMat = copy(pList) self.trainedModel.data.marginalInterp = self._setupMarginalInterp() self.trainedModel.data.ARBsList = ARBsList self.trainedModel.data.bRBsList = bRBsList self.trainedModel.data.approxParameters = copy(self.approxParameters) vbMng(self, "DEL", "Done setting up approximant.", 5) def assembleReducedSystemMarginal(self, pMat : sampList = None)\ -> Tuple[List[List[Np2D]], List[List[Np1D]], List[Np2D]]: """Build affine blocks of RB linear system through projections.""" if pMat is None: self.setupApprox() ARBsList = self.trainedModel.data.ARBsList bRBsList = self.trainedModel.data.bRBsList else: vbMng(self, "INIT", "Projecting affine terms of HF model.", 10) ARBsList = [None] * len(self.musMarginal) bRBsList = [None] * len(self.musMarginal) projList = [None] * len(self.musMarginal) for k, (As, bs, sample) in enumerate(zip(self.AsList, self.bsList, self.samplingEngine.samples)): compLocal = self.HFEngine.innerProduct(sample, pMat) projList[k] = sample.dot(customPInv(compLocal)) ARBsList[k], bRBsList[k] = projectAffineDecomposition(As, bs, projList[k]) vbMng(self, "DEL", "Done projecting affine terms.", 10) return ARBsList, bRBsList, projList diff --git a/rrompy/reduction_methods/pivoting/rb_distributed_pivoted.py b/rrompy/reduction_methods/pivoting/rb_distributed_pivoted_pole_matching.py similarity index 56% copy from rrompy/reduction_methods/pivoting/rb_distributed_pivoted.py copy to rrompy/reduction_methods/pivoting/rb_distributed_pivoted_pole_matching.py index af9c165..19afb42 100644 --- a/rrompy/reduction_methods/pivoting/rb_distributed_pivoted.py +++ b/rrompy/reduction_methods/pivoting/rb_distributed_pivoted_pole_matching.py @@ -1,290 +1,198 @@ # 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 . # from copy import deepcopy as copy -import numpy as np -from .generic_pivoted_approximant import GenericPivotedApproximant -from rrompy.hfengines.base import MarginalProxyEngine +from .rb_distributed_pivoted import RBDistributedPivoted from rrompy.reduction_methods.trained_model import (TrainedModelPivotedData, - TrainedModelPivotedRB as tModel) -from rrompy.reduction_methods.base.rb_utils import projectAffineDecomposition -from rrompy.utilities.base.types import (Np1D, Np2D, List, ListAny, Tuple, - DictAny, HFEng, paramVal, sampList) -from rrompy.utilities.base import verbosityManager as vbMng, freepar as fp -from rrompy.utilities.poly_fitting import customPInv -from rrompy.utilities.exception_manager import (RROMPyWarning, RROMPyException, - RROMPyAssert) + TrainedModelPivotedRBPoleMatching as tModel) +from rrompy.utilities.base.types import HFEng, DictAny, ListAny, paramVal +from rrompy.utilities.base import verbosityManager as vbMng +from rrompy.utilities.exception_manager import (RROMPyException, RROMPyAssert) -__all__ = ['RBDistributedPivoted'] +__all__ = ['RBDistributedPivotedPoleMatching'] -class RBDistributedPivoted(GenericPivotedApproximant): +class RBDistributedPivotedPoleMatching(RBDistributedPivoted): """ - ROM pivoted RB approximant computation for parametric problems. + ROM pivoted RB approximant computation for parametric problems with pole + matching. Args: HFEngine: HF problem solver. mu0(optional): Default parameter. Defaults to 0. directionPivot(optional): Pivot components. Defaults to [0]. approxParameters(optional): Dictionary containing values for main parameters of approximant. Recognized keys are: - - 'POD': whether to compute POD of snapshots; defaults to True; + - 'matchingWeight': weight for pole matching optimization; defaults + to 1; - 'S': total number of pivot samples current approximant relies upon; - 'samplerPivot': pivot sample point generator; - 'SMarginal': total number of marginal samples current approximant relies upon; - 'samplerMarginal': marginal sample point generator; - 'R': rank for pivot Galerkin projection; defaults to prod(S); - 'PODTolerance': tolerance for pivot snapshots POD; defaults to -1; - 'polybasisMarginal': type of polynomial basis for marginal interpolation; allowed values include 'MONOMIAL', 'CHEBYSHEV' and 'LEGENDRE'; defaults to 'MONOMIAL'; - 'MMarginal': degree of marginal interpolant; defaults to 0; - 'radialBasisMarginal': radial basis family for marginal interpolant; defaults to 0, i.e. no radial basis; - 'radialBasisWeightsMarginal': radial basis weights for marginal interpolant; defaults to 0, i.e. identity; - 'interpRcondMarginal': tolerance for marginal interpolation; defaults to None. Defaults to empty dict. homogeneized(optional): Whether to homogeneize Dirichlet BCs. Defaults to False. verbosity(optional): Verbosity level. Defaults to 10. Attributes: HFEngine: HF problem solver. mu0: Default parameter. directionPivot: Pivot components. mus: Array of snapshot parameters. + musPivot: Array of pivot snapshot parameters. + musMarginal: Array of marginal snapshot parameters. homogeneized: Whether to homogeneize Dirichlet BCs. - approxRadius: Dummy radius of approximant (i.e. distance from mu0 to - farthest sample point). approxParameters: Dictionary containing values for main parameters of approximant. Recognized keys are in parameterList. parameterListSoft: Recognized keys of soft approximant parameters: - 'POD': whether to compute POD of snapshots; + - 'matchingWeight': weight for pole matching optimization; - 'R': rank for Galerkin projection; - 'PODTolerance': tolerance for snapshots POD; - 'polybasisMarginal': type of polynomial basis for marginal interpolation; - 'MMarginal': degree of marginal interpolant; - 'radialBasisMarginal': radial basis family for marginal interpolant; - 'radialBasisWeightsMarginal': radial basis weights for marginal interpolant; - 'interpRcondMarginal': tolerance for marginal interpolation. parameterListCritical: Recognized keys of critical approximant parameters: - 'S': total number of pivot samples current approximant relies upon; - 'samplerPivot': pivot sample point generator; - 'SMarginal': total number of marginal samples current approximant relies upon; - 'samplerMarginal': marginal sample point generator. POD: Whether to compute POD of snapshots. + matchingWeight: Weight for pole matching optimization. S: Total number of pivot samples current approximant relies upon. samplerPivot: Pivot sample point generator. SMarginal: Total number of marginal samples current approximant relies upon. samplerMarginal: Marginal sample point generator. R: Rank for Galerkin projection. PODTolerance: Tolerance for pivot snapshots POD. polybasisMarginal: Type of polynomial basis for marginal interpolation. MMarginal: Degree of marginal interpolant. radialBasisMarginal: Radial basis family for marginal interpolant. radialBasisWeightsMarginal: Radial basis weights for marginal interpolant. interpRcondMarginal: Tolerance for marginal interpolation. muBoundsPivot: list of bounds for pivot parameter values. muBoundsMarginal: list of bounds for marginal parameter values. samplingEngine: Sampling engine. uHF: High fidelity solution(s) with parameter(s) lastSolvedHF as sampleList. lastSolvedHF: Parameter(s) corresponding to last computed high fidelity solution(s) as parameterList. + lastSolvedHF: Parameter(s) corresponding to last computed high fidelity + solution(s) as parameterList. uApproxReduced: Reduced approximate solution(s) with parameter(s) lastSolvedApprox as sampleList. lastSolvedApproxReduced: Parameter(s) corresponding to last computed reduced approximate solution(s) as parameterList. uApprox: Approximate solution(s) with parameter(s) lastSolvedApprox as sampleList. lastSolvedApprox: Parameter(s) corresponding to last computed approximate solution(s) as parameterList. As: List of sparse matrices (in CSC format) representing coefficients - of linear system matrix wrt theta(mu). + of linear system matrix. bs: List of numpy vectors representing coefficients of linear system - RHS wrt theta(mu). - thetaAs: List of callables representing coefficients of linear system - matrix wrt mu. - thetabs: List of callables representing coefficients of linear system - RHS wrt mu. + RHS. ARBs: List of sparse matrices (in CSC format) representing coefficients - of compressed linear system matrix wrt theta(mu). + of compressed linear system matrix. bRBs: List of numpy vectors representing coefficients of compressed - linear system RHS wrt theta(mu). + linear system RHS. """ def __init__(self, HFEngine:HFEng, mu0 : paramVal = None, directionPivot : ListAny = [0], approxParameters : DictAny = {}, homogeneized : bool = False, verbosity : int = 10, timestamp : bool = True): self._preInit() - self._addParametersToList(["R", "PODTolerance"], [1, -1]) + if len(directionPivot) > 1: + raise RROMPyException(("Exactly 1 pivot parameter allowed in pole " + "matching.")) + self._addParametersToList(["matchingWeight"], [1]) super().__init__(HFEngine = HFEngine, mu0 = mu0, directionPivot = directionPivot, approxParameters = approxParameters, homogeneized = homogeneized, verbosity = verbosity, timestamp = timestamp) self._postInit() @property - def S(self): - """Value of S.""" - return self._S - @S.setter - def S(self, S): - GenericPivotedApproximant.S.fset(self, S) - if not hasattr(self, "_R"): - self._R = np.prod(self.S) * np.prod(self.SMarginal) - - @property - def R(self): - """Value of R. Its assignment may change S.""" - return self._R - @R.setter - def R(self, R): - if R < 0: raise RROMPyException("R must be non-negative.") - self._R = R - self._approxParameters["R"] = self.R - if (hasattr(self, "_S") and hasattr(self, "_SMarginal") - and np.prod(self.S) * np.prod(self.SMarginal) < self.R): - RROMPyWarning(("Prescribed S and SMarginal are too small. " - "Decreasing R.")) - self.R = np.prod(self.S) * np.prod(self.SMarginal) - - @property - def PODTolerance(self): - """Value of PODTolerance.""" - return self._PODTolerance - @PODTolerance.setter - def PODTolerance(self, PODTolerance): - self._PODTolerance = PODTolerance - self._approxParameters["PODTolerance"] = self.PODTolerance + def matchingWeight(self): + """Value of matchingWeight.""" + return self._matchingWeight + @matchingWeight.setter + def matchingWeight(self, matchingWeight): + self._matchingWeight = matchingWeight + self._approxParameters["matchingWeight"] = self.matchingWeight - def _setupProjectionMatrix(self): - """Compute projection matrix.""" - RROMPyAssert(self._mode, message = "Cannot setup numerator.") - vbMng(self, "INIT", "Starting computation of projection matrix.", 7) - if self.POD: - U, s, _ = np.linalg.svd(self.samplingEngine.RPODCoalesced) - s = s ** 2. - else: - Gramian = self.HFEngine.innerProduct( - self.samplingEngine.samplesCoalesced, - self.samplingEngine.samplesCoalesced) - U, s, _ = np.linalg.svd(Gramian) - nsamples = self.samplingEngine.samplesCoalesced.shape[1] - snorm = np.cumsum(s[::-1]) / np.sum(s) - nPODTrunc = min(nsamples - np.argmax(snorm > self.PODTolerance), - self.R) - pMat = self.samplingEngine.samplesCoalesced.dot(U[:, : nPODTrunc]) - vbMng(self, "MAIN", - "Assembling {}x{} projection matrix from {} samples.".format( - *(pMat.shape), nsamples), 5) - vbMng(self, "DEL", "Done computing projection matrix.", 7) - return pMat - - def _setupAffineBlocks(self): - """Compute list of marginalized affine blocks of system.""" - hasAs = hasattr(self, "AsList") and self.AsList is not None - hasbs = hasattr(self, "bsList") and self.bsList is not None - if hasAs and hasbs: return - vbMng(self, "INIT", "Computing affine blocks of system.", 10) - mu0Eff = self.mu0.data[0, self.directionPivot] - if not hasAs: self.AsList = [None] * len(self.musMarginal) - if not hasbs: self.bsList = [None] * len(self.musMarginal) - for k, muM in enumerate(self.musMarginal): - muEff = [fp] * self.npar - for jj, kk in enumerate(self.directionMarginal): - muEff[kk] = muM(0, jj) - MEnginek = MarginalProxyEngine(self.HFEngine, muEff) - if not hasAs: - self.AsList[k] = MEnginek.affineLinearSystemA(mu0Eff) - if not hasbs: - self.bsList[k] = MEnginek.affineLinearSystemb(mu0Eff, - self.homogeneized) - self.thetaAs = MEnginek.affineWeightsA(mu0Eff) - self.thetabs = MEnginek.affineWeightsb(mu0Eff, self.homogeneized) - vbMng(self, "DEL", "Done computing affine blocks.", 10) - def setupApprox(self): """Compute RB projection matrix.""" if self.checkComputedApprox(): return RROMPyAssert(self._mode, message = "Cannot setup approximant.") vbMng(self, "INIT", "Setting up {}.".format(self.name()), 5) self.computeScaleFactor() self.computeSnapshots() self._setupAffineBlocks() pMat = self._setupProjectionMatrix() ARBsList, bRBsList, pList = self.assembleReducedSystemMarginal(pMat) if self.trainedModel is None: self.trainedModel = tModel() self.trainedModel.verbosity = self.verbosity self.trainedModel.timestamp = self.timestamp data = TrainedModelPivotedData(self.trainedModel.name(), self.mu0, pList, self.HFEngine.rescalingExp, self.directionPivot, self.scaleFactor) data.musPivot = copy(self.musPivot) data.musMarginal = copy(self.musMarginal) - data.thetaAs, data.thetabs = self.thetaAs, self.thetabs self.trainedModel.data = data else: self.trainedModel = self.trainedModel self.trainedModel.data.projMat = copy(pList) self.trainedModel.data.marginalInterp = self._setupMarginalInterp() self.trainedModel.data.ARBsList = ARBsList self.trainedModel.data.bRBsList = bRBsList + vbMng(self, "INIT", "Matching poles.", 10) + self.trainedModel.initializeFromAffine(self.HFEngine, + self.matchingWeight, self.POD) + vbMng(self, "DEL", "Done matching poles.", 10) self.trainedModel.data.approxParameters = copy(self.approxParameters) vbMng(self, "DEL", "Done setting up approximant.", 5) - - def assembleReducedSystemMarginal(self, pMat : sampList = None)\ - -> Tuple[List[List[Np2D]], - List[List[Np1D]], List[Np2D]]: - """Build affine blocks of RB linear system through projections.""" - if pMat is None: - self.setupApprox() - ARBsList = self.trainedModel.data.ARBsList - bRBsList = self.trainedModel.data.bRBsList - else: - vbMng(self, "INIT", "Projecting affine terms of HF model.", 10) - ARBsList = [None] * len(self.musMarginal) - bRBsList = [None] * len(self.musMarginal) - projList = [None] * len(self.musMarginal) - for k, (As, bs, sample) in enumerate(zip(self.AsList, self.bsList, - self.samplingEngine.samples)): - compLocal = self.HFEngine.innerProduct(sample, pMat) - projList[k] = sample.dot(customPInv(compLocal)) - ARBsList[k], bRBsList[k] = projectAffineDecomposition(As, bs, - projList[k]) - vbMng(self, "DEL", "Done projecting affine terms.", 10) - return ARBsList, bRBsList, projList - diff --git a/rrompy/reduction_methods/trained_model/__init__.py b/rrompy/reduction_methods/trained_model/__init__.py index e77d9fd..e16c89a 100644 --- a/rrompy/reduction_methods/trained_model/__init__.py +++ b/rrompy/reduction_methods/trained_model/__init__.py @@ -1,37 +1,45 @@ # 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 . # from .trained_model import TrainedModel from .trained_model_data import TrainedModelData from .trained_model_rational import TrainedModelRational from .trained_model_rb import TrainedModelRB from .trained_model_pivoted_data import TrainedModelPivotedData from .trained_model_pivoted_rational import TrainedModelPivotedRational from .trained_model_pivoted_rb import TrainedModelPivotedRB +from .trained_model_pivoted_general_pole_matching import \ + TrainedModelPivotedGeneralPoleMatching +from .trained_model_pivoted_rational_pole_matching import \ + TrainedModelPivotedRationalPoleMatching +from .trained_model_pivoted_rb_pole_matching import \ + TrainedModelPivotedRBPoleMatching __all__ = [ 'TrainedModel', 'TrainedModelData', 'TrainedModelRational', 'TrainedModelRB', 'TrainedModelPivotedData', 'TrainedModelPivotedRational', - 'TrainedModelPivotedRB' + 'TrainedModelPivotedRB', + 'TrainedModelPivotedRationalPoleMatching', + 'TrainedModelPivotedRBPoleMatching' ] diff --git a/rrompy/reduction_methods/trained_model/trained_model_pivoted_general_pole_matching.py b/rrompy/reduction_methods/trained_model/trained_model_pivoted_general_pole_matching.py new file mode 100644 index 0000000..8e71380 --- /dev/null +++ b/rrompy/reduction_methods/trained_model/trained_model_pivoted_general_pole_matching.py @@ -0,0 +1,196 @@ +# 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 . +# + +import numpy as np +from .trained_model import TrainedModel +from rrompy.utilities.base.types import (Np1D, ListAny, paramVal, paramList, + sampList, HFEng) +from rrompy.utilities.base import verbosityManager as vbMng, freepar as fp +from rrompy.utilities.numerical import pointMatching +from rrompy.utilities.poly_fitting.heaviside import HeavisideInterpolator as HI +from rrompy.utilities.exception_manager import RROMPyException, RROMPyAssert +from rrompy.parameter import checkParameter, checkParameterList +from rrompy.sampling import emptySampleList + +__all__ = ['TrainedModelPivotedGeneralPoleMatching'] + +class TrainedModelPivotedGeneralPoleMatching(TrainedModel): + """ + ROM approximant evaluation for pivoted approximants with pole matching. + + Attributes: + Data: dictionary with all that can be pickled. + """ + + def initializeFromLists(self, poles:ListAny, coeffs:ListAny, basis:str, + HFEngine:HFEng, matchingWeight : float = 1., + POD : bool = True): + """Initialize Heaviside representation.""" + musM = self.data.musMarginal + margAbsDist = np.sum(np.abs(np.tile(musM.data.reshape(len(musM),1,-1), + [1, len(musM), 1]) + - np.tile(musM.data.reshape(1,len(musM),-1), + [len(musM), 1, 1])), axis = 2) + N = len(poles[0]) + explored = np.zeros(1, dtype = int) + unexplored = np.arange(1, len(musM)) + for _ in range(1, len(musM)): + minIdx = np.argmin(np.concatenate([margAbsDist[ex, unexplored] \ + for ex in explored])) + minIex = minIdx // len(explored) + minIunex = minIdx % len(explored) + dist = np.abs(np.tile(poles[minIex].reshape(-1, 1), N) + - poles[minIunex].reshape(1, -1)) + if matchingWeight != 0: + resex = coeffs[minIex][: N] + resunex = coeffs[minIunex][: N] + if POD: + distR = resex.dot(resunex.T.conj()) + distR = (distR.T / np.linalg.norm(resex, axis = 1)).T + distR = distR / np.linalg.norm(resunex, axis = 1) + else: + resex = self.data.projMat.dot(resex.T) + resunex = self.data.projMat.dot(resunex.T) + distR = HFEngine.innerProduct(resex, resunex).T + distR = (distR.T / HFEngine.norm(resex)).T + distR = distR / HFEngine.norm(resunex) + distR = np.abs(distR) + distR[distR > 1.] = 1. + dist += 2. / np.pi * matchingWeight * np.arccos(distR) + reordering = pointMatching(dist, verbObj = self) + poles[minIunex] = poles[minIunex][reordering] + coeffs[minIunex][: N] = coeffs[minIunex][reordering] + explored = np.append(explored, minIunex) + unexplored = unexplored[unexplored != minIunex] + HIs = [] + for pls, cfs in zip(poles, coeffs): + hsi = HI() + hsi.poles = pls + hsi.coeffs = cfs + hsi.npar = 1 + hsi.polybasis = basis + HIs += [hsi] + self.data.HIs = HIs + + def interpolateMarginalInterpolator(self, mu : paramVal = []) -> Np1D: + """Obtain interpolated approximant interpolator.""" + mu = checkParameter(mu, self.data.nparMarginal)[0] + hsi = HI() + hsi.poles = self.interpolateMarginalPoles(mu) + hsi.coeffs = self.interpolateMarginalCoeffs(mu) + hsi.npar = 1 + hsi.polybasis = self.data.HIs[0].polybasis + return hsi + + def interpolateMarginalPoles(self, mu : paramList = []) -> Np1D: + """Obtain interpolated approximant poles.""" + mu = checkParameterList(mu, self.data.nparMarginal)[0] + return self.interpolateMarginal(mu, [hi.poles for hi in self.data.HIs]) + + def interpolateMarginalCoeffs(self, mu : paramList = []) -> Np1D: + """Obtain interpolated approximant coefficients.""" + mu = checkParameterList(mu, self.data.nparMarginal)[0] + cs = self.interpolateMarginal(mu, [hi.coeffs for hi in self.data.HIs]) + return cs.reshape(self.data.HIs[0].coeffs.shape) + + def getApproxReduced(self, mu : paramList = []) -> sampList: + """ + Evaluate reduced representation of approximant at arbitrary parameter. + + Args: + mu: Target parameter. + """ + RROMPyAssert(self.data.nparPivot, 1, "Number of pivot parameters") + mu = checkParameterList(mu, self.data.npar)[0] + if (not hasattr(self, "lastSolvedApproxReduced") + or self.lastSolvedApproxReduced != mu): + vbMng(self, "INIT", + "Evaluating approximant at mu = {}.".format(mu), 12) + self.uApproxReduced = emptySampleList() + self.uApproxReduced.reset((self.data.HIs[0].coeffs.shape[1], + len(mu)), self.data.projMat[0].dtype) + + for i, muPL in enumerate(mu): + muL = self.centerNormalizePivot([muPL(0, x) \ + for x in self.data.directionPivot]) + muM = [muPL(0, x) for x in self.data.directionMarginal] + vbMng(self, "INIT", + "Assembling reduced model for mu = {}.".format(muPL), 17) + hsL = self.interpolateMarginalInterpolator(muM) + vbMng(self, "DEL", "Done assembling reduced model.", 17) + self.uApproxReduced[i] = hsL(muL) + vbMng(self, "DEL", "Done evaluating approximant.", 12) + self.lastSolvedApproxReduced = mu + return self.uApproxReduced + + def getPoles(self, *args, **kwargs) -> Np1D: + """ + Obtain approximant poles. + + Returns: + Numpy complex vector of poles. + """ + RROMPyAssert(self.data.nparPivot, 1, "Number of pivot parameters") + if len(args) + len(kwargs) > 1: + raise RROMPyException(("Wrong number of parameters passed. " + "Only 1 available.")) + elif len(args) + len(kwargs) == 1: + if len(args) == 1: + mVals = args[0] + else: + mVals = kwargs["marginalVals"] + if not hasattr(mVals, "__len__"): mVals = [mVals] + mVals = list(mVals) + else: + mVals = [fp] + try: + rDim = mVals.index(fp) + if rDim < len(mVals) - 1 and fp in mVals[rDim + 1 :]: + raise + except: + raise RROMPyException(("Exactly 1 'freepar' entry in " + "marginalVals must be provided.")) + if rDim != self.data.directionPivot[0]: + raise RROMPyException(("'freepar' entry in marginalVals must " + "coincide with pivot direction.")) + mVals[rDim] = self.data.mu0(rDim) + mMarg = [mVals[j] for j in range(len(mVals)) if j != rDim] + roots = np.array(self.interpolateMarginalPoles(mMarg)) + return np.power(self.data.mu0(rDim) ** self.data.rescalingExp[rDim] + + self.data.scaleFactor[rDim] * roots, + 1. / self.data.rescalingExp[rDim]) + + def getResidues(self, *args, **kwargs) -> Np1D: + """ + Obtain approximant residues. + + Returns: + Numpy matrix with residues as columns. + """ + pls = self.getPoles(*args, **kwargs) + if len(args) == 1: + mVals = args[0] + else: + mVals = kwargs["marginalVals"] + if not hasattr(mVals, "__len__"): mVals = [mVals] + mVals = list(mVals) + rDim = mVals.index(fp) + mMarg = [mVals[j] for j in range(len(mVals)) if j != rDim] + residues = self.interpolateMarginalCoeffs(mMarg)[: len(pls)] + res = self.data.projMat.dot(residues) + return pls, res diff --git a/rrompy/reduction_methods/trained_model/trained_model_pivoted_rational.py b/rrompy/reduction_methods/trained_model/trained_model_pivoted_rational.py index 1ae9dec..36bdbd4 100644 --- a/rrompy/reduction_methods/trained_model/trained_model_pivoted_rational.py +++ b/rrompy/reduction_methods/trained_model/trained_model_pivoted_rational.py @@ -1,238 +1,224 @@ # 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 . # import numpy as np from .trained_model_rational import TrainedModelRational from rrompy.utilities.base.types import (Np1D, List, ListAny, paramVal, paramList, sampList) from rrompy.utilities.base import verbosityManager as vbMng, freepar as fp from rrompy.utilities.poly_fitting.polynomial import polyroots from rrompy.utilities.exception_manager import RROMPyException from rrompy.parameter import checkParameter, checkParameterList from rrompy.sampling import sampleList, emptySampleList __all__ = ['TrainedModelPivotedRational'] class TrainedModelPivotedRational(TrainedModelRational): """ ROM approximant evaluation for pivoted Rational approximant. Attributes: Data: dictionary with all that can be pickled. """ def centerNormalizePivot(self, mu : paramList = [], mu0 : paramVal = None) -> paramList: """ Compute normalized parameter to be plugged into approximant. Args: mu: Parameter(s) 1. mu0: Parameter(s) 2. If None, set to self.data.mu0. Returns: Normalized parameter. """ mu = checkParameterList(mu, self.data.nparPivot)[0] if mu0 is None: mu0 = self.data.mu0(self.data.directionPivot) rad = ((mu ** self.data.rescalingExpPivot - mu0 ** self.data.rescalingExpPivot) / self.data.scaleFactorPivot) return rad def centerNormalizeMarginal(self, mu : paramList = [], mu0 : paramVal = None) -> paramList: """ Compute normalized parameter to be plugged into approximant. Args: mu: Parameter(s) 1. mu0: Parameter(s) 2. If None, set to self.data.mu0. Returns: Normalized parameter. """ mu = checkParameterList(mu, self.data.nparMarginal)[0] if mu0 is None: mu0 = self.data.mu0(self.data.directionMarginal) rad = ((mu ** self.data.rescalingExpMarginal - mu0 ** self.data.rescalingExpMarginal) / self.data.scaleFactorMarginal) return rad def interpolateMarginal(self, mu : paramList = [], samples : ListAny = [], der : List[int] = None, scl : Np1D = None) -> sampList: """ Evaluate marginal interpolator at arbitrary marginal parameter. Args: mu: Target parameter. samples: Objects to interpolate. der(optional): Derivatives to take before evaluation. """ mu = checkParameterList(mu, self.data.nparMarginal)[0] sList = isinstance(samples[0], sampleList) sEff = [None] * len(samples) for j in range(len(samples)): if sList: sEff[j] = samples[j].data else: sEff[j] = samples[j] vbMng(self, "INIT", "Interpolating marginal at mu = {}.".format(mu), 35) muC = self.centerNormalizeMarginal(mu) p = emptySampleList() p.reset((len(sEff[0]), len(muC))) p.data[:] = 0. for mIj, spj in zip(self.data.marginalInterp, sEff): - p = p + spj * mIj(muC, der, scl) + p = p + spj.reshape(len(sEff[0]), - 1) * mIj(muC, der, scl) vbMng(self, "DEL", "Done interpolating marginal.", 35) if not sList: p = p.data.flatten() return p - def getPValsPivot(self, mu : paramList = [], der : List[int] = None, - scl : Np1D = None) -> sampList: + def getPValsPivot(self, mu : paramList = []) -> sampList: """ Evaluate rational numerators at arbitrary pivot parameter. Args: mu: Target parameter. - der(optional): Derivatives to take before evaluation. """ mu = checkParameterList(mu, self.data.nparPivot)[0] vbMng(self, "INIT", "Evaluating numerator at mu = {}.".format(mu), 17) muC = self.centerNormalizePivot(mu) - pP = [sampleList(P(muC, der, scl)) for P in self.data.Ps] + pP = [sampleList(P(muC)) for P in self.data.Ps] vbMng(self, "DEL", "Done evaluating numerator.", 17) return pP - def getPVal(self, mu : paramList = [], der : List[int] = None, - scl : Np1D = None) -> sampList: + def getPVal(self, mu : paramList = []) -> sampList: """ Evaluate rational numerator at arbitrary parameter. Args: mu: Target parameter. - der(optional): Derivatives to take before evaluation. """ mu = checkParameterList(mu, self.data.npar)[0] muP = checkParameterList(mu.data[:, self.data.directionPivot], self.data.nparPivot)[0] muM = checkParameterList(mu.data[:, self.data.directionMarginal], self.data.nparMarginal)[0] - if der is None: - derP, derM = None, None - else: - derP = [der[x] for x in self.data.directionPivot] - derM = [der[x] for x in self.data.directionMarginal] - if scl is None: - sclP, sclM = None, None - else: - sclP = [scl[x] for x in self.data.directionPivot] - sclM = [scl[x] for x in self.data.directionMarginal] - pP = self.getPValsPivot(muP, derP, sclP) - return self.interpolateMarginal(muM, pP, derM, sclM) + pP = self.getPValsPivot(muP) + return self.interpolateMarginal(muM, pP) def getQValsPivot(self, mu:Np1D, der : List[int] = None, scl : Np1D = None) -> Np1D: """ Evaluate rational denominators at arbitrary pivot parameter. Args: mu: Target parameter. der(optional): Derivatives to take before evaluation. """ mu = checkParameterList(mu, self.data.nparPivot)[0] vbMng(self, "INIT", "Evaluating denominator at mu = {}.".format(mu), 17) muC = self.centerNormalizePivot(mu) qP = [Q(muC, der, scl).reshape(1, -1) for Q in self.data.Qs] vbMng(self, "DEL", "Done evaluating denominator.", 17) return qP def getQVal(self, mu:Np1D, der : List[int] = None, scl : Np1D = None) -> Np1D: """ Evaluate rational denominator at arbitrary parameter. Args: mu: Target parameter. der(optional): Derivatives to take before evaluation. """ mu = checkParameterList(mu, self.data.npar)[0] muP = checkParameterList(mu.data[:, self.data.directionPivot], self.data.nparPivot)[0] muM = checkParameterList(mu.data[:, self.data.directionMarginal], self.data.nparMarginal)[0] if der is None: derP, derM = None, None else: derP = [der[x] for x in self.data.directionPivot] derM = [der[x] for x in self.data.directionMarginal] if scl is None: sclP, sclM = None, None else: sclP = [scl[x] for x in self.data.directionPivot] sclM = [scl[x] for x in self.data.directionMarginal] qP = self.getQValsPivot(muP, derP, sclP) return self.interpolateMarginal(muM, qP, derM, sclM) def getPoles(self, *args, **kwargs) -> Np1D: """ Obtain approximant poles. Returns: Numpy complex vector of poles. """ if len(args) + len(kwargs) > 1: raise RROMPyException(("Wrong number of parameters passed. " "Only 1 available.")) elif len(args) + len(kwargs) == 1: if len(args) == 1: mVals = args[0] else: mVals = kwargs["marginalVals"] if not hasattr(mVals, "__len__"): mVals = [mVals] mVals = list(mVals) else: mVals = [fp] mValsP = [mVals[x] for x in self.data.directionPivot] mValsM = [mVals[x] for x in self.data.directionMarginal] try: rDim = mValsP.index(fp) if rDim < len(mValsP) - 1 and fp in mValsP[rDim + 1 :]: raise except: raise RROMPyException(("Exactly 1 'freepar' entry in " "marginalVals must be provided.")) if fp in mValsM: raise RROMPyException(("'freepar' entry in marginalVals must be " "pivot dimension.")) rDimB = self.data.directionPivot[rDim] Qeff = self.interpolateMarginal(mValsM, [Q.coeffs[:, np.newaxis] \ for Q in self.data.Qs]) Qeff = Qeff.reshape(self.data.Qs[0].coeffs.shape) mValsP[rDim] = self.data.mu0(rDimB) mValsP = self.centerNormalizePivot(checkParameter(mValsP)) mValsP = list(mValsP.data.flatten()) mValsP[rDim] = fp return np.power(self.data.mu0(rDimB) ** self.data.rescalingExp[rDimB] + self.data.scaleFactor[rDimB] * polyroots(Qeff, self.data.Qs[0].polybasis, mValsP), 1. / self.data.rescalingExp[rDimB]) diff --git a/rrompy/reduction_methods/trained_model/trained_model_pivoted_rational_pole_matching.py b/rrompy/reduction_methods/trained_model/trained_model_pivoted_rational_pole_matching.py new file mode 100644 index 0000000..a7e9fba --- /dev/null +++ b/rrompy/reduction_methods/trained_model/trained_model_pivoted_rational_pole_matching.py @@ -0,0 +1,110 @@ +# 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 . +# + +import numpy as np +from scipy.special import factorial as fact +from itertools import combinations +from .trained_model_pivoted_general_pole_matching import \ + TrainedModelPivotedGeneralPoleMatching +from .trained_model_pivoted_rational import TrainedModelPivotedRational +from rrompy.utilities.base.types import Np1D, List, paramList, sampList, HFEng +from rrompy.utilities.poly_fitting.heaviside import rational2heaviside +from rrompy.utilities.exception_manager import RROMPyException, RROMPyAssert +from rrompy.parameter import checkParameterList +from rrompy.sampling import emptySampleList + +__all__ = ['TrainedModelPivotedRationalPoleMatching'] + +class TrainedModelPivotedRationalPoleMatching( + TrainedModelPivotedGeneralPoleMatching, + TrainedModelPivotedRational): + """ + ROM approximant evaluation for pivoted Rational approximant with pole + matching. + + Attributes: + Data: dictionary with all that can be pickled. + """ + + def initializeFromRational(self, HFEngine:HFEng, + matchingWeight : float = 1., POD : bool = True): + """Initialize Heaviside representation.""" + RROMPyAssert(self.data.nparPivot, 1, "Number of pivot parameters") + poles, coeffs = [], [] + for Q, P in zip(self.data.Qs, self.data.Ps): + cfs, pls, basis = rational2heaviside(P, Q) + poles += [pls] + coeffs += [cfs] + self.initializeFromLists(poles, coeffs, basis, HFEngine, + matchingWeight, POD) + self.data._temporary = False + + def getPVal(self, mu : paramList = []) -> sampList: + """ + Evaluate rational numerator at arbitrary parameter. + + Args: + mu: Target parameter. + """ + if self.data._temporary: return super().getPVal(mu) + RROMPyAssert(self.data.nparPivot, 1, "Number of pivot parameters") + mu = checkParameterList(mu, self.data.npar)[0] + p = emptySampleList() + p.reset((len(self.data.HIs[0].coeffs.shape[1]), len(mu))) + for i, muPL in enumerate(mu): + muL = self.centerNormalizePivot([muPL(0, x) \ + for x in self.data.directionPivot]) + muM = [muPL(0, x) for x in self.data.directionMarginal] + hsL = self.interpolateMarginalInterpolator(muM) + p[i] = hsL(muL) * np.prod(muL(0, 0) - hsL.poles) + return p + + def getQVal(self, mu:Np1D, der : List[int] = None, + scl : Np1D = None) -> Np1D: + """ + Evaluate rational denominator at arbitrary parameter. + + Args: + mu: Target parameter. + der(optional): Derivatives to take before evaluation. + """ + if self.data._temporary: return super().getQVal(mu, der, scl) + RROMPyAssert(self.data.nparPivot, 1, "Number of pivot parameters") + mu = checkParameterList(mu, self.data.npar)[0] + muP = self.centerNormalizePivot(checkParameterList( + mu.data[:, self.data.directionPivot], 1)[0]) + muM = checkParameterList(mu.data[:, self.data.directionMarginal], + self.data.npar - 1)[0] + if der is None: + derP, derM = 0, [0] + else: + derP = der[self.data.directionPivot[0]] + derM = [der[x] for x in self.data.directionMarginal] + if np.any(np.array(derM) != 0): + raise RROMPyException(("Derivatives of Q with respect to marginal " + "parameters not allowed.")) + sclP = 1 if scl is None else scl[self.data.directionPivot[0]] + derVal = np.zeros(len(mu), dtype = self.data.HIs[0].poles.dtype) + N = len(self.data.HIs[0].poles) + if derP == N: derVal[:] = 1. + elif derP >= 0 and derP < N: + pls = self.interpolateMarginalPoles(muM).reshape(-1, len(mu)).T + plsDist = muP.data.reshape(-1, 1) - pls + for terms in combinations(np.arange(N), N - derP): + derVal += np.prod(plsDist[:, list(terms)], axis = 1) + return sclP ** derP * fact(derP) * derVal diff --git a/rrompy/reduction_methods/trained_model/trained_model_pivoted_rb.py b/rrompy/reduction_methods/trained_model/trained_model_pivoted_rb.py index c5979a2..7a5b260 100644 --- a/rrompy/reduction_methods/trained_model/trained_model_pivoted_rb.py +++ b/rrompy/reduction_methods/trained_model/trained_model_pivoted_rb.py @@ -1,196 +1,205 @@ # 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 . # import numpy as np -from scipy.linalg import eigvals from .trained_model_rb import TrainedModelRB from rrompy.utilities.base.types import (Np1D, ListAny, paramVal, paramList, sampList) from rrompy.utilities.base import verbosityManager as vbMng, freepar as fp -from rrompy.utilities.exception_manager import (RROMPyException, RROMPyWarning, - RROMPyAssert) +from rrompy.utilities.poly_fitting.polynomial import ( + hashIdxToDerivative as hashI) +from rrompy.utilities.numerical import (eigvalsNonlinearDense, + marginalizePolyList) +from rrompy.utilities.exception_manager import RROMPyException from rrompy.parameter import checkParameter, checkParameterList from rrompy.sampling import emptySampleList, sampleList __all__ = ['TrainedModelPivotedRB'] class TrainedModelPivotedRB(TrainedModelRB): """ ROM approximant evaluation for pivoted RB approximant. Attributes: Data: dictionary with all that can be pickled. """ + def centerNormalizePivot(self, mu : paramList = [], + mu0 : paramVal = None) -> paramList: + """ + Compute normalized parameter to be plugged into approximant. + + Args: + mu: Parameter(s) 1. + mu0: Parameter(s) 2. If None, set to self.data.mu0. + + Returns: + Normalized parameter. + """ + mu = checkParameterList(mu, self.data.nparPivot)[0] + if mu0 is None: mu0 = self.data.mu0(self.data.directionPivot) + rad = ((mu ** self.data.rescalingExpPivot + - mu0 ** self.data.rescalingExpPivot)) + return rad + def centerNormalizeMarginal(self, mu : paramList = [], mu0 : paramVal = None) -> paramList: """ Compute normalized parameter to be plugged into approximant. Args: mu: Parameter(s) 1. mu0: Parameter(s) 2. If None, set to self.data.mu0. Returns: Normalized parameter. """ mu = checkParameterList(mu, self.data.nparMarginal)[0] if mu0 is None: mu0 = self.data.mu0(self.data.directionMarginal) rad = ((mu ** self.data.rescalingExpMarginal - mu0 ** self.data.rescalingExpMarginal) / self.data.scaleFactorMarginal) return rad def interpolateMarginal(self, mu : paramVal = [], samples : ListAny = []) -> ListAny: """ Evaluate marginal interpolator at arbitrary marginal parameter. Args: mu: Target parameter. samples: Objects to interpolate. """ mu = checkParameter(mu, self.data.nparMarginal) sList = isinstance(samples[0], sampleList) sEff = [None] * len(samples) for j in range(len(samples)): if sList: sEff[j] = [samples[j].data] else: sEff[j] = samples[j] vbMng(self, "INIT", "Interpolating marginal at mu = {}.".format(mu), 35) muC = self.centerNormalizeMarginal(mu) p = [None] * len(sEff[0]) for k in range(len(sEff[0])): p[k] = np.zeros_like(sEff[0][k], dtype = sEff[0][k].dtype) for mIj, spj in zip(self.data.marginalInterp, sEff): mIjEval = mIj(muC) for k in range(len(sEff[0])): p[k] = p[k] + mIjEval * spj[k] vbMng(self, "DEL", "Done interpolating marginal.", 35) return p def getApproxReduced(self, mu : paramList = []) -> sampList: """ Evaluate reduced representation of approximant at arbitrary parameter. Args: mu: Target parameter. """ - Mu = checkParameterList(mu, self.data.npar)[0] + mu = checkParameterList(mu, self.data.npar)[0] if (not hasattr(self, "lastSolvedApproxReduced") - or self.lastSolvedApproxReduced != Mu): + or self.lastSolvedApproxReduced != mu): vbMng(self, "INIT", - "Computing RB solution at mu = {}.".format(Mu), 12) - thetaAs, thetabs = self.data.thetaAs, self.data.thetabs + "Computing RB solution at mu = {}.".format(mu), 12) self.uApproxReduced = emptySampleList() self.uApproxReduced.reset((self.data.ARBsList[0][0].shape[0], len(mu)), self.data.projMat[0].dtype) - for i, muPL in enumerate(Mu): - mu = checkParameter([muPL(0, x) \ + mu0Eff = (self.data.mu0Pivot ** self.data.rescalingExpPivot).data + for i, muPL in enumerate(mu): + muP = checkParameter([muPL(0, x) \ for x in self.data.directionPivot]) muM = [muPL(0, x) for x in self.data.directionMarginal] vbMng(self, "INIT", "Assembling reduced model for mu = {}.".format(muPL), 17) ARBs = self.interpolateMarginal(muM, self.data.ARBsList) bRBs = self.interpolateMarginal(muM, self.data.bRBsList) - ARBmu = eval(thetaAs[0]) * ARBs[0] - bRBmu = eval(thetabs[0]) * bRBs[0] - for j in range(1, len(ARBs)): - ARBmu += eval(thetaAs[j]) * ARBs[j] - for j in range(1, len(bRBs)): - bRBmu += eval(thetabs[j]) * bRBs[j] + ARBmu = ARBs[0] + bRBmu = bRBs[0] + for j in range(1, max(len(ARBs), len(bRBs))): + theta = np.prod((muP ** self.data.rescalingExpPivot + - mu0Eff) ** hashI(j, self.data.nparPivot)) + if j < len(ARBs): + ARBmu += theta * ARBs[j] + if j < len(bRBs): + bRBmu += theta * bRBs[j] vbMng(self, "DEL", "Done assembling reduced model.", 17) vbMng(self, "INIT", "Solving reduced model for mu = {}.".format(muPL), 15) self.uApproxReduced[i] = np.linalg.solve(ARBmu, bRBmu) vbMng(self, "DEL", "Done solving reduced model.", 15) vbMng(self, "DEL", "Done computing RB solution.", 12) - self.lastSolvedApproxReduced = Mu + self.lastSolvedApproxReduced = mu return self.uApproxReduced def getApprox(self, mu : paramList = []) -> sampList: """ Evaluate approximant at arbitrary parameter. Args: mu: Target parameter. """ mu = checkParameterList(mu, self.data.npar)[0] if (not hasattr(self, "lastSolvedApprox") or self.lastSolvedApprox != mu): uApproxR = self.getApproxReduced(mu) self.uApprox = emptySampleList() self.uApprox.reset((self.data.projMat[0].shape[0], len(mu)), self.data.projMat[0].dtype) for i in range(len(mu)): muM = [mu(i, x) for x in self.data.directionMarginal] projLoc = self.interpolateMarginal(muM, self.data.projMat)[0] self.uApprox[i] = projLoc.dot(uApproxR[i]) self.lastSolvedApprox = mu return self.uApprox - def getPoles(self, marginalVals : ListAny = [fp], *args, **kwargs) -> Np1D: + def getPoles(self, marginalVals : ListAny = [fp], jSupp : int = 1, + **kwargs) -> Np1D: """ Obtain approximant poles. Returns: Numpy complex vector of poles. """ - RROMPyAssert(self.data.nparPivot, 1, "Number of parameters") - RROMPyWarning(("Impossible to compute poles in general affine " - "parameter dependence. Results subject to " - "interpretation/rescaling, or possibly completely " - "wrong.")) if not hasattr(marginalVals, "__len__"): marginalVals = [marginalVals] marginalVals = list(marginalVals) try: rDim = marginalVals.index(fp) if rDim < len(marginalVals) - 1 and fp in marginalVals[rDim + 1 :]: raise except: raise RROMPyException(("Exactly 1 'freepar' entry in " "marginalVals must be provided.")) - mValsM = [marginalVals[x] for x in self.data.directionMarginal] - try: - rDim = marginalVals.index(fp) - if rDim < len(marginalVals) - 1 and fp in marginalVals[rDim + 1 :]: - raise - except: - raise RROMPyException(("Exactly 1 'freepar' entry in " - "marginalVals must be provided.")) - if fp in mValsM: + if rDim in self.data.directionMarginal: raise RROMPyException(("'freepar' entry in marginalVals must be " - "pivot dimension.")) - rDimB = self.data.directionPivot[rDim] + "in pivot dimension.")) + mValsP = [marginalVals[x] for x in self.data.directionPivot] + rDimP = mValsP.index(fp) + mValsM = [marginalVals[x] for x in self.data.directionMarginal] ARBs = self.interpolateMarginal(mValsM, self.data.ARBsList) - R = ARBs[0].shape[0] - if len(ARBs) < 2: - return - A = np.eye(R * (len(ARBs) - 1), dtype = np.complex) - B = np.zeros_like(A) - A[: R, : R] = - ARBs[0] - for j in range(len(ARBs) - 1): - Aj = ARBs[j + 1] - B[: R, j * R : (j + 1) * R] = Aj - II = np.arange(R, R * (len(ARBs) - 1)) - B[II, II - R] = 1. - return np.power(self.data.mu0(rDimB) ** self.data.rescalingExp[rDimB] - + eigvals(A, B, *args, **kwargs), - 1. / self.data.rescalingExp[rDimB]) + if self.data.nparPivot > 1: + mValsP[rDimP] = self.data.mu0Pivot(rDimP) + mValsP = (checkParameter(mValsP) ** self.data.rescalingExpPivot + - self.data.mu0Pivot ** self.data.rescalingExpPivot) + mValsP = mValsP.data.flatten() + mValsP[rDimP] = fp + ARBs = marginalizePolyList(ARBs, mValsP, "auto") + ev = eigvalsNonlinearDense(ARBs, jSupp = jSupp, **kwargs) + return np.power(self.data.mu0(rDim) ** self.data.rescalingExp[rDim] + + ev, 1. / self.data.rescalingExp[rDim]) diff --git a/rrompy/reduction_methods/trained_model/trained_model_pivoted_rb_pole_matching.py b/rrompy/reduction_methods/trained_model/trained_model_pivoted_rb_pole_matching.py new file mode 100644 index 0000000..608bdf4 --- /dev/null +++ b/rrompy/reduction_methods/trained_model/trained_model_pivoted_rb_pole_matching.py @@ -0,0 +1,47 @@ +# 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 . +# + +from .trained_model_pivoted_general_pole_matching import \ + TrainedModelPivotedGeneralPoleMatching +from .trained_model_pivoted_rb import TrainedModelPivotedRB +from rrompy.utilities.base.types import HFEng +from rrompy.utilities.poly_fitting.heaviside import affine2heaviside +from rrompy.utilities.exception_manager import RROMPyAssert + +__all__ = ['TrainedModelPivotedRBPoleMatching'] + +class TrainedModelPivotedRBPoleMatching(TrainedModelPivotedGeneralPoleMatching, + TrainedModelPivotedRB): + """ + ROM approximant evaluation for pivoted RB approximant with pole matching. + + Attributes: + Data: dictionary with all that can be pickled. + """ + + def initializeFromAffine(self, HFEngine:HFEng, + matchingWeight : float = 1., POD : bool = True): + """Initialize Heaviside representation.""" + RROMPyAssert(self.data.nparPivot, 1, "Number of pivot parameters") + poles, coeffs = [], [] + for As, bs in zip(self.data.ARBsList, self.data.bRBsList): + cfs, pls, basis = affine2heaviside(As, bs, self.data.mu0Pivot) + poles += [pls] + coeffs += [cfs] + self.initializeFromLists(poles, coeffs, basis, HFEngine, + matchingWeight, POD) diff --git a/rrompy/reduction_methods/trained_model/trained_model_rational.py b/rrompy/reduction_methods/trained_model/trained_model_rational.py index 7d34787..8f5ec51 100644 --- a/rrompy/reduction_methods/trained_model/trained_model_rational.py +++ b/rrompy/reduction_methods/trained_model/trained_model_rational.py @@ -1,165 +1,163 @@ # 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 . # import numpy as np from . import TrainedModel from rrompy.utilities.base.types import (Np1D, List, paramVal, paramList, sampList) from rrompy.utilities.base import verbosityManager as vbMng, freepar as fp from rrompy.utilities.exception_manager import RROMPyException from rrompy.parameter import (checkParameter, checkParameterList, emptyParameterList) from rrompy.sampling import sampleList __all__ = ['TrainedModelRational'] class TrainedModelRational(TrainedModel): """ ROM approximant evaluation for rational approximant. Attributes: Data: dictionary with all that can be pickled. """ def centerNormalize(self, mu : paramList = [], mu0 : paramVal = None) -> paramList: """ Compute normalized parameter to be plugged into approximant. Args: mu: Parameter(s) 1. mu0: Parameter(s) 2. If None, set to self.data.mu0. Returns: Normalized parameter. """ mu = checkParameterList(mu, self.data.npar)[0] if mu0 is None: mu0 = self.data.mu0 rad = ((mu ** self.data.rescalingExp - mu0 ** self.data.rescalingExp) / self.data.scaleFactor) return rad - def getPVal(self, mu : paramList = [], der : List[int] = None, - scl : Np1D = None) -> sampList: + def getPVal(self, mu : paramList = []) -> sampList: """ Evaluate rational numerator at arbitrary parameter. Args: mu: Target parameter. - der(optional): Derivatives to take before evaluation. """ mu = checkParameterList(mu, self.data.npar)[0] vbMng(self, "INIT", "Evaluating numerator at mu = {}.".format(mu), 17) muCenter = self.centerNormalize(mu) - p = sampleList(self.data.P(muCenter, der, scl)) + p = sampleList(self.data.P(muCenter)) vbMng(self, "DEL", "Done evaluating numerator.", 17) return p def getQVal(self, mu:Np1D, der : List[int] = None, scl : Np1D = None) -> Np1D: """ Evaluate rational denominator at arbitrary parameter. Args: mu: Target parameter. der(optional): Derivatives to take before evaluation. """ mu = checkParameterList(mu, self.data.npar)[0] vbMng(self, "INIT", "Evaluating denominator at mu = {}.".format(mu), 17) muCenter = self.centerNormalize(mu) q = self.data.Q(muCenter, der, scl) vbMng(self, "DEL", "Done evaluating denominator.", 17) return q def getApproxReduced(self, mu : paramList = []) -> sampList: """ Evaluate reduced representation of approximant at arbitrary parameter. Args: mu: Target parameter. """ mu = checkParameterList(mu, self.data.npar)[0] if (not hasattr(self, "lastSolvedApproxReduced") or self.lastSolvedApproxReduced != mu): vbMng(self, "INIT", "Evaluating approximant at mu = {}.".format(mu), 12) self.uApproxReduced = self.getPVal(mu) / self.getQVal(mu) vbMng(self, "DEL", "Done evaluating approximant.", 12) self.lastSolvedApproxReduced = mu return self.uApproxReduced def getPoles(self, *args, **kwargs) -> Np1D: """ Obtain approximant poles. Returns: Numpy complex vector of poles. """ if len(args) + len(kwargs) > 1: raise RROMPyException(("Wrong number of parameters passed. " "Only 1 available.")) elif len(args) + len(kwargs) == 1: if len(args) == 1: mVals = args[0] else: mVals = kwargs["marginalVals"] if not hasattr(mVals, "__len__"): mVals = [mVals] mVals = list(mVals) else: mVals = [fp] try: rDim = mVals.index(fp) if rDim < len(mVals) - 1 and fp in mVals[rDim + 1 :]: raise except: raise RROMPyException(("Exactly 1 'freepar' entry in " "marginalVals must be provided.")) mVals[rDim] = self.data.mu0(rDim) mVals = self.centerNormalize(checkParameter(mVals, len(mVals))) mVals = list(mVals.data.flatten()) mVals[rDim] = fp return np.power(self.data.mu0(rDim) ** self.data.rescalingExp[rDim] + self.data.scaleFactor[rDim] * self.data.Q.roots(mVals), 1. / self.data.rescalingExp[rDim]) def getResidues(self, *args, **kwargs) -> Np1D: """ Obtain approximant residues. Returns: Numpy matrix with residues as columns. """ pls = self.getPoles(*args, **kwargs) if len(args) == 1: mVals = args[0] else: mVals = kwargs["marginalVals"] if not hasattr(mVals, "__len__"): mVals = [mVals] mVals = list(mVals) rDim = mVals.index(fp) poles = emptyParameterList() poles.reset((len(pls), self.data.npar), dtype = pls.dtype) for k, pl in enumerate(pls): poles[k] = mVals poles.data[k, rDim] = pl res = (self.data.projMat.dot(self.getPVal(poles).data) / self.getQVal(poles, list(1 * (np.arange(self.data.npar) == rDim)))) return pls, res diff --git a/rrompy/reduction_methods/trained_model/trained_model_rb.py b/rrompy/reduction_methods/trained_model/trained_model_rb.py index f508b3e..e656919 100644 --- a/rrompy/reduction_methods/trained_model/trained_model_rb.py +++ b/rrompy/reduction_methods/trained_model/trained_model_rb.py @@ -1,111 +1,106 @@ # 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 . # import numpy as np -from scipy.linalg import eigvals from .trained_model import TrainedModel from rrompy.utilities.base.types import Np1D, ListAny, paramList, sampList from rrompy.utilities.base import verbosityManager as vbMng, freepar as fp -from rrompy.utilities.exception_manager import (RROMPyException, RROMPyWarning, - RROMPyAssert) -from rrompy.parameter import checkParameterList +from rrompy.utilities.poly_fitting.polynomial import ( + hashIdxToDerivative as hashI) +from rrompy.utilities.numerical import (eigvalsNonlinearDense, + marginalizePolyList) +from rrompy.utilities.exception_manager import RROMPyException +from rrompy.parameter import checkParameter, checkParameterList from rrompy.sampling import emptySampleList __all__ = ['TrainedModelRB'] class TrainedModelRB(TrainedModel): """ ROM approximant evaluation for RB approximant. Attributes: Data: dictionary with all that can be pickled. """ def getApproxReduced(self, mu : paramList = []) -> sampList: """ Evaluate reduced representation of approximant at arbitrary parameter. Args: mu: Target parameter. """ - Mu = checkParameterList(mu, self.data.npar)[0] + mu = checkParameterList(mu, self.data.npar)[0] if (not hasattr(self, "lastSolvedApproxReduced") - or self.lastSolvedApproxReduced != Mu): + or self.lastSolvedApproxReduced != mu): vbMng(self, "INIT", - "Computing RB solution at mu = {}.".format(Mu), 12) - thetaAs, thetabs = self.data.thetaAs, self.data.thetabs + "Computing RB solution at mu = {}.".format(mu), 12) ARBs, bRBs = self.data.ARBs, self.data.bRBs self.uApproxReduced = emptySampleList() self.uApproxReduced.reset((ARBs[0].shape[0], len(mu)), self.data.projMat.dtype) - for i, muPL in enumerate(Mu): - mu = muPL[0] - vbMng(self, "INIT", - "Assembling reduced model for mu = {}.".format(mu), 17) - ARBmu = eval(thetaAs[0]) * ARBs[0] - bRBmu = eval(thetabs[0]) * bRBs[0] - for j in range(1, len(ARBs)): - ARBmu += eval(thetaAs[j]) * ARBs[j] - for j in range(1, len(bRBs)): - bRBmu += eval(thetabs[j]) * bRBs[j] + mu0Eff = (self.data.mu0 ** self.data.rescalingExp).data + for i in range(len(mu)): + vbMng(self, "INIT", "Assembling reduced model for mu = {}."\ + .format(mu[i]), 17) + ARBmu = 1. * ARBs[0] + bRBmu = 1. * bRBs[0] + for j in range(1, max(len(ARBs), len(bRBs))): + theta = np.prod((mu[i] ** self.data.rescalingExp + - mu0Eff) ** hashI(j, self.data.npar)) + if j < len(ARBs): + ARBmu += theta * ARBs[j] + if j < len(bRBs): + bRBmu += theta * bRBs[j] vbMng(self, "DEL", "Done assembling reduced model.", 17) vbMng(self, "INIT", - "Solving reduced model for mu = {}.".format(mu), 15) + "Solving reduced model for mu = {}.".format(mu[i]), 15) self.uApproxReduced[i] = np.linalg.solve(ARBmu, bRBmu) vbMng(self, "DEL", "Done solving reduced model.", 15) vbMng(self, "DEL", "Done computing RB solution.", 12) - self.lastSolvedApproxReduced = Mu + self.lastSolvedApproxReduced = mu return self.uApproxReduced - def getPoles(self, marginalVals : ListAny = [fp], *args, **kwargs) -> Np1D: + def getPoles(self, marginalVals : ListAny = [fp], jSupp : int = 1, + **kwargs) -> Np1D: """ Obtain approximant poles. Returns: Numpy complex vector of poles. """ - RROMPyAssert(self.data.npar, 1, "Number of parameters") - RROMPyWarning(("Impossible to compute poles in general affine " - "parameter dependence. Results subject to " - "interpretation/rescaling, or possibly completely " - "wrong.")) if not hasattr(marginalVals, "__len__"): marginalVals = [marginalVals] - marginalVals = list(marginalVals) + mVals = list(marginalVals) try: - rDim = marginalVals.index(fp) - if rDim < len(marginalVals) - 1 and fp in marginalVals[rDim + 1 :]: + rDim = mVals.index(fp) + if rDim < len(mVals) - 1 and fp in mVals[rDim + 1 :]: raise except: raise RROMPyException(("Exactly 1 'freepar' entry in " "marginalVals must be provided.")) ARBs = self.data.ARBs - R = ARBs[0].shape[0] - if len(ARBs) < 2: - return - A = np.eye(R * (len(ARBs) - 1), dtype = np.complex) - B = np.zeros_like(A) - A[: R, : R] = - ARBs[0] - for j in range(len(ARBs) - 1): - Aj = ARBs[j + 1] - B[: R, j * R : (j + 1) * R] = Aj - II = np.arange(R, R * (len(ARBs) - 1)) - B[II, II - R] = 1. + if self.data.npar > 1: + mVals[rDim] = self.data.mu0(rDim) + mVals = (checkParameter(mVals) ** self.data.rescalingExp + - self.data.mu0 ** self.data.rescalingExp) + mVals = mVals.data.flatten() + mVals[rDim] = fp + ARBs = marginalizePolyList(ARBs, mVals, "auto") + ev = eigvalsNonlinearDense(ARBs, jSupp = jSupp, **kwargs) return np.power(self.data.mu0(rDim) ** self.data.rescalingExp[rDim] - + eigvals(A, B, *args, **kwargs), - 1. / self.data.rescalingExp[rDim]) - + + ev, 1. / self.data.rescalingExp[rDim]) diff --git a/rrompy/utilities/base/__init__.py b/rrompy/utilities/base/__init__.py index 99a50e7..c2be650 100644 --- a/rrompy/utilities/base/__init__.py +++ b/rrompy/utilities/base/__init__.py @@ -1,59 +1,42 @@ # 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 . # from .find_dict_str_key import findDictStrKey from .get_new_filename import getNewFilename -from .kroneckerer import kroneckerer -from .factorials import multibinom, multifactorial from .pickle_utilities import pickleDump, pickleLoad from .purge_dict import purgeDict from .purge_list import purgeList -from .number_theory import (squareResonances, primeFactorize, - getLowestPrimeFactor) -from .halton import haltonGenerate -from .sobol import sobolGenerate -from .low_discrepancy import vanderCorput, lowDiscrepancy from . import types as Types from .verbosity_depth import verbosityDepth, verbosityManager freepar = None __all__ = [ 'findDictStrKey', 'getNewFilename', - 'kroneckerer', - 'multibinom', - 'multifactorial', 'pickleDump', 'pickleLoad', 'purgeDict', 'purgeList', - 'squareResonances', - 'primeFactorize', - 'getLowestPrimeFactor', - 'haltonGenerate', - 'sobolGenerate', - 'vanderCorput', - 'lowDiscrepancy', 'Types', 'verbosityDepth', 'verbosityManager', 'freepar' ] diff --git a/rrompy/utilities/base/__init__.py b/rrompy/utilities/numerical/__init__.py similarity index 56% copy from rrompy/utilities/base/__init__.py copy to rrompy/utilities/numerical/__init__.py index 99a50e7..6af5afa 100644 --- a/rrompy/utilities/base/__init__.py +++ b/rrompy/utilities/numerical/__init__.py @@ -1,59 +1,49 @@ # 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 . # -from .find_dict_str_key import findDictStrKey -from .get_new_filename import getNewFilename -from .kroneckerer import kroneckerer +from .custom_pinv import customPInv from .factorials import multibinom, multifactorial -from .pickle_utilities import pickleDump, pickleLoad -from .purge_dict import purgeDict -from .purge_list import purgeList -from .number_theory import (squareResonances, primeFactorize, - getLowestPrimeFactor) from .halton import haltonGenerate +from .kroneckerer import kroneckerer +from .low_discrepancy import lowDiscrepancy +from .marginalize_poly_list import marginalizePolyList +from .nonlinear_eigenproblem import (linearizeDense, eigNonlinearDense, + eigvalsNonlinearDense) +from .number_theory import squareResonances +from .point_matching import pointMatching from .sobol import sobolGenerate -from .low_discrepancy import vanderCorput, lowDiscrepancy -from . import types as Types -from .verbosity_depth import verbosityDepth, verbosityManager freepar = None __all__ = [ - 'findDictStrKey', - 'getNewFilename', - 'kroneckerer', + 'customPInv', 'multibinom', 'multifactorial', - 'pickleDump', - 'pickleLoad', - 'purgeDict', - 'purgeList', - 'squareResonances', - 'primeFactorize', - 'getLowestPrimeFactor', 'haltonGenerate', - 'sobolGenerate', - 'vanderCorput', + 'kroneckerer', 'lowDiscrepancy', - 'Types', - 'verbosityDepth', - 'verbosityManager', - 'freepar' + 'marginalizePolyList', + 'linearizeDense', + 'eigNonlinearDense', + 'eigvalsNonlinearDense', + 'squareResonances', + 'pointMatching', + 'sobolGenerate' ] diff --git a/rrompy/utilities/poly_fitting/custom_pinv.py b/rrompy/utilities/numerical/custom_pinv.py similarity index 100% rename from rrompy/utilities/poly_fitting/custom_pinv.py rename to rrompy/utilities/numerical/custom_pinv.py diff --git a/rrompy/utilities/base/factorials.py b/rrompy/utilities/numerical/factorials.py similarity index 100% rename from rrompy/utilities/base/factorials.py rename to rrompy/utilities/numerical/factorials.py diff --git a/rrompy/utilities/base/halton.py b/rrompy/utilities/numerical/halton.py similarity index 100% rename from rrompy/utilities/base/halton.py rename to rrompy/utilities/numerical/halton.py diff --git a/rrompy/utilities/base/kroneckerer.py b/rrompy/utilities/numerical/kroneckerer.py similarity index 100% rename from rrompy/utilities/base/kroneckerer.py rename to rrompy/utilities/numerical/kroneckerer.py diff --git a/rrompy/utilities/base/low_discrepancy.py b/rrompy/utilities/numerical/low_discrepancy.py similarity index 96% rename from rrompy/utilities/base/low_discrepancy.py rename to rrompy/utilities/numerical/low_discrepancy.py index 0dc73eb..e765fee 100644 --- a/rrompy/utilities/base/low_discrepancy.py +++ b/rrompy/utilities/numerical/low_discrepancy.py @@ -1,42 +1,42 @@ # 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 . # import numpy as np from rrompy.utilities.base.types import List from rrompy.utilities.exception_manager import RROMPyException -__all__ = ['vanderCorput', 'lowDiscrepancy'] +__all__ = ['lowDiscrepancy'] def vanderCorput(n:int) -> List[int]: if n <= 0: raise RROMPyException("Only positive integers allowed.") x = [0] * n ln = int(np.ceil(np.log2(n))) for j in range(n): x[j] = int(np.binary_repr(j, width = ln)[::-1], 2) return x def lowDiscrepancy(n:int, inverse : bool = False) -> List[int]: if n <= 0: raise RROMPyException("Only positive integers allowed.") max2Fac = 2 ** int(np.ceil(np.log2(n))) xBase = np.array(vanderCorput(max2Fac), dtype = np.int) x = list(xBase[xBase < n]) if inverse: x = list(np.argsort(x)) return x diff --git a/rrompy/utilities/numerical/marginalize_poly_list.py b/rrompy/utilities/numerical/marginalize_poly_list.py new file mode 100644 index 0000000..9780834 --- /dev/null +++ b/rrompy/utilities/numerical/marginalize_poly_list.py @@ -0,0 +1,79 @@ +# 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 . +# + +import numpy as np +from scipy.sparse import csr +from rrompy.utilities.base.types import Np1D, Np2D, ListAny +from rrompy.utilities.base import freepar as fp +from rrompy.utilities.poly_fitting.polynomial import ( + hashDerivativeToIdx as hashD, hashIdxToDerivative as hashI) +from rrompy.parameter import checkParameter + +__all__ = ['marginalizePolyList'] + +def marginalizePolyList(objs:ListAny, marginalVals : Np1D = [fp], + zeroObj : Np2D = 0., + recompress : bool = True) -> ListAny: + res = [] + freeLocations = [] + fixedLocations = [] + muFixed = [] + if not hasattr(marginalVals, "__len__"): marginalVals = [marginalVals] + for i, m in enumerate(marginalVals): + if m == fp: + freeLocations += [i] + else: + fixedLocations += [i] + muFixed += [m] + muFixed = checkParameter(muFixed, len(fixedLocations)) + if zeroObj == "auto": + if isinstance(objs[0], np.ndarray): + zeroObj = np.zeros_like(objs[0]) + elif isinstance(objs[0], csr.csr_matrix): + d = objs[0].shape[0] + zeroObj = csr.csr_matrix(([], [], np.zeros(d + 1)), + shape = objs[0].shape, + dtype = objs[0].dtype) + else: + zeroObj = 0. + for j, obj in enumerate(objs): + derjBase = hashI(j, len(marginalVals)) + jNew = hashD([derjBase[i] for i in freeLocations]) + derjFixed = [derjBase[i] for i in fixedLocations] + obj = np.prod((muFixed ** derjFixed).data) * obj + if jNew >= len(res): + for _ in range(len(res), jNew): + res += [zeroObj] + res += [obj] + else: + res[jNew] = res[jNew] + obj + if recompress: + for re in res[::-1]: + try: + if isinstance(re, np.ndarray): + iszero = np.allclose(re, zeroObj, + atol = 2 * np.finfo(re.dtype).eps) + elif isinstance(re, csr.csr_matrix): + iszero = re.nnz == 0 + else: + break + if not iszero: break + except: break + res.pop() + return res + \ No newline at end of file diff --git a/rrompy/utilities/numerical/nonlinear_eigenproblem.py b/rrompy/utilities/numerical/nonlinear_eigenproblem.py new file mode 100644 index 0000000..e43ba93 --- /dev/null +++ b/rrompy/utilities/numerical/nonlinear_eigenproblem.py @@ -0,0 +1,108 @@ +# 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 . +# + +import numpy as np +import scipy.linalg as scla +#import scipy.sparse as scsp +from rrompy.utilities.base.types import Tuple, List, Np1D, Np2D +from .custom_pinv import customPInv + +__all__ = ['linearizeDense', 'eigNonlinearDense', 'eigvalsNonlinearDense'] + +def linearizeDense(As:List[Np2D], jSupp : int = 1) -> Tuple[Np2D, Np2D]: + N = len(As) + n = As[0].shape[0] + stiff = np.zeros(((N - 1) * n, (N - 1) * n), dtype = As[0].dtype) + mass = np.zeros(((N - 1) * n, (N - 1) * n), dtype = As[0].dtype) + if isinstance(jSupp, str) and jSupp.upper() == "COMPANION": + II = np.arange(n, (N - 1) * n) + stiff = np.pad(- np.hstack(As[-2 :: -1]), + [[0, (N - 2) * n], [0, 0]], "constant") + stiff[II, II - n] = 1. + mass = np.pad(As[-1], [0, (N - 2) * n], "constant") + mass[II, II] = 1. + else: + for j in range(jSupp): + for k in range(jSupp - j - 1, jSupp): + mass[n * j : n * (j + 1), k * n : (k + 1) * n] = \ + As[N - 2 + jSupp - k - j] + for j in range(jSupp - 1, N - 1): + for k in range(jSupp, N - 1 + jSupp - j): + stiff[n * j : n * (j + 1), (k - 1) * n : k * n] = \ + - As[jSupp - k + N - 2 - j] + stiff[: n * (jSupp - 1), : n * (jSupp - 1)] = mass[: n * (jSupp - 1), + n : n * jSupp] + mass[n * jSupp :, n * jSupp :] = stiff[n * (jSupp - 1) : - n, + n * jSupp :] + return stiff, mass + +def eigNonlinearDense(As:List[Np2D], jSupp : int = 1, + return_inverse : bool = False, + **kwargs_eig) -> Tuple[Np1D, Np2D]: + stiff, mass = linearizeDense(As, jSupp) + ev, P = scla.eig(stiff, mass, overwrite_a = True, overwrite_b = True, + **kwargs_eig) + if not return_inverse: return ev, P + Pinv = customPInv(P) + return ev, P, Pinv + +def eigvalsNonlinearDense(As:List[Np2D], jSupp : int = 1, + **kwargs_eigvals) -> Np1D: + stiff, mass = linearizeDense(As, jSupp) + return scla.eigvals(stiff, mass, overwrite_a = True, **kwargs_eigvals) + +#def linearizeSparse(As:List[Np2D], jSupp : int = 1) -> Tuple[Np2D, Np2D]: +# N = len(As) +# n = As[0].shape[0] +# if isinstance(jSupp, str) and jSupp.upper() == "COMPANION": +# II = np.arange(n, (N - 1) * n) +# III = np.arange((N - 2) * n + 1) +# IIII = np.arange(0, n ** 2, n) +# improve management of sparse As... +# Alist = - np.hstack([A.todense() for A in As[-2 :: -1]]) +# stiffD = np.concatenate((Alist.flatten(), np.ones((N - 2) * n))) +# stiffP = np.concatenate(((N - 1) * IIII, (N - 1) * n ** 2 + III)) +# stiffI = np.concatenate((np.tile(np.arange((N - 1) * n), n), II - n)) +# massD = np.concatenate((As[-1].todense().flatten(), +# np.ones((N - 2) * n))) +# massP = np.concatenate((IIII, n ** 2 + III)) +# massI = np.concatenate((np.tile(np.arange(n), n), II)) +# else: +# compute stiffD, stiffP, stiffI depending on jSupp +# compute massD, massP, massI depending on jSupp +# stiff = scsp.csr_matrix((stiffD, stiffI, stiffP), +# shape = ((N - 1) * n, (N - 1) * n)) +# mass = scsp.csr_matrix((massD, massI, massP), +# shape = ((N - 1) * n, (N - 1) * n)) +# return stiff, mass +# +#def eigNonlinearSparse(As:List[Np2D], jSupp : int = 1, +# return_inverse : bool = False, +# **kwargs_eig) -> Tuple[Np1D, Np2D]: +# stiff, mass = linearizeSparse(As, jSupp) +# ev, P = scsp.linalg.eig(stiff, M = mass, return_eigenvectors = True, +# **kwargs_eig) +# if not return_inverse: return ev, P +# Pinv = customPInv(P) +# return ev, P, Pinv +# +#def eigvalsNonlinearSparse(As:List[Np2D], jSupp : int = 1, +# **kwargs_eigvals) -> Np1D: +# stiff, mass = linearizeSparse(As, jSupp) +# return scsp.linalg.eig(stiff, M = mass, return_eigenvectors = False, +# **kwargs_eigvals) diff --git a/rrompy/utilities/base/number_theory.py b/rrompy/utilities/numerical/number_theory.py similarity index 96% rename from rrompy/utilities/base/number_theory.py rename to rrompy/utilities/numerical/number_theory.py index c71c9df..cc728ab 100644 --- a/rrompy/utilities/base/number_theory.py +++ b/rrompy/utilities/numerical/number_theory.py @@ -1,70 +1,70 @@ # 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 . # import numpy as np -__all__ = ['squareResonances','getLowestPrimeFactor','primeFactorize'] +__all__ = ['squareResonances'] prime_v = [] def squareResonances(a:int, b:int, zero_terms : bool = True): spectrum = [] a = max(int(np.floor(a)), 0) b = max(int(np.ceil(b)), 0) global prime_v if len(prime_v) == 0: prime_v = [2, 3] if a > prime_v[-1]: for i in range(prime_v[-1], a, 2): getLowestPrimeFactor(i) for i in range(a, b + 1): spectrum = spectrum + [i] * countSquareSums(i, zero_terms) return np.array(spectrum) def getLowestPrimeFactor(n:int): global prime_v for x in prime_v: if n % x == 0: return x if prime_v[-1] < n: prime_v = prime_v + [n] return n def primeFactorize(n:int): factors = [] number = n while number > 1: factor = getLowestPrimeFactor(number) factors.append(factor) number = int(number / factor) if n < -1: factors[0] = -factors[0] return list(factors) def countSquareSums(n:int, zero_terms : bool = True): if n < 2: return (n + 1) * zero_terms factors = primeFactorize(n) funique, fcounts = np.unique(factors, return_counts = True) count = 1 for fac, con in zip(funique, fcounts): #using number theory magic if fac % 4 == 1: count = count * (con + 1) elif fac % 4 == 3 and con % 2 == 1: return 0 return count + (2 * zero_terms - 1) * (int(n ** .5) ** 2 == n) diff --git a/rrompy/utilities/numerical/point_matching.py b/rrompy/utilities/numerical/point_matching.py new file mode 100644 index 0000000..b3c9cfa --- /dev/null +++ b/rrompy/utilities/numerical/point_matching.py @@ -0,0 +1,156 @@ +# 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 . +# + +import numpy as np +from itertools import permutations +from rrompy.utilities.base.types import Np1D, Np2D +from rrompy.utilities.base import verbosityManager as vbMng +from rrompy.utilities.exception_manager import RROMPyException + +__all__ = ['pointMatching'] + +def matchExplicit(distanceMatrix:Np2D, idx1:Np1D, idx2:Np1D, + verbObj = None) -> Np1D: + if verbObj is not None: + vbMng(verbObj, "INIT", ("Starting exact point matching of {} " + "points.").format(len(distanceMatrix)), 25) + L1 = len(idx1) + matches, optV = None, None + idxLead = 0 + for idx2p in permutations(idx2, L1): + if verbObj is not None and idxLead != idx2p[0]: + idxLead += 1 + explRatio = 100 * idxLead / len(idx2) + vbMng(verbObj, "MAIN", + "Explored {}% of permutations.".format(int(explRatio)), 65) + val = np.sum(distanceMatrix[idx1, idx2p]) + if optV is None or val < optV: + optV, matches = val, idx2p + if verbObj is not None: + vbMng(verbObj, "DEL", "Done point matching.", 25) + return np.array(matches) + +def findClusterBipartiteRecursive(i, edges, clusters, explored, right): + if explored[right][i] == 0: return clusters, explored + explored[right][i] = 0 + if edges[right][i] not in clusters[1 - right]: + clusters[1 - right] += [edges[right][i]] + clusters, explored = findClusterBipartiteRecursive(edges[right][i], edges, + clusters, explored, + 1 - right) + for k, i2 in enumerate(edges[1 - right]): + if i2 == i and k not in clusters[1 - right]: + clusters[1 - right] += [k] + clusters, explored = findClusterBipartiteRecursive(k, edges, + clusters, explored, + 1 - right) + return clusters, explored + +def pointMatchingHeuristic(distanceMatrix:Np2D, idx1:Np1D, idx2:Np1D, + max_iter : int = 10, verbObj = None) -> Np1D: + if verbObj is not None: + vbMng(verbObj, "INIT", ("Starting heuristic point matching of {} " + "points.").format(len(idx1)), 25) + distanceMatrixEff = distanceMatrix[idx1, :] + distanceMatrixEff = distanceMatrixEff[:, idx2] + N = len(distanceMatrixEff) + matches = - np.ones(N, dtype = int) + fBest = np.argmin(distanceMatrixEff, axis = 1) + bBest = np.argmin(distanceMatrixEff, axis = 0) + clusters1, clusters2 = [], [] + fIdxs = [np.ones(N, dtype = bool), np.ones(N, dtype = bool)] + if verbObj is not None: vbMng(verbObj, "INIT", "Starting clustering.", 65) + for i in range(N): + if fIdxs[0][i]: + cloc, fIdxs = findClusterBipartiteRecursive(i, [fBest, bBest], + [[], []], fIdxs, 0) + clusters1 += [cloc[0]] + clusters2 += [cloc[1]] + if verbObj is not None: vbMng(verbObj, "DEL", "Done clustering.", 65) + if verbObj is not None: + vbMng(verbObj, "INIT", "Starting optimization of clustered points.", + 65) + for c1, c2 in zip(clusters1, clusters2): + if len(c1) > len(c2): + optP = matchExplicit(distanceMatrixEff.T, c2, c1) + for i, k in enumerate(optP): + matches[k] = c2[i] + else: + optP = matchExplicit(distanceMatrixEff, c1, c2) + for i, k in enumerate(c1): + matches[k] = optP[i] + front1 = np.where(matches == -1)[0] + if len(front1) == 0: + if verbObj is not None: vbMng(verbObj, "DEL", "Done optimizing.", 65) + if verbObj is not None: + vbMng(verbObj, "DEL", "Done point matching.", 25) + return idx2[matches] + if len(front1) == N: + raise RROMPyException(("Heuristic point matching algorithm not " + "converging. Must increase threshold.")) + optP = np.array([i for i in range(N) if i not in matches], dtype = int) + if verbObj is not None: vbMng(verbObj, "DEL", "Done optimizing.", 65) + for _ in range(max_iter): + bulk1 = [i for i in range(N) if i not in front1] + if verbObj is not None: + vbMng(verbObj, "INIT", ("Starting optimization of {} unclustered " + "points.").format(len(front1)), 65) + optP = pointMatching(distanceMatrixEff, idx1 = front1, idx2 = optP) + for i, k in enumerate(front1): matches[k] = optP[i] + if verbObj is not None: vbMng(verbObj, "DEL", "Done optimizing.", 65) + if verbObj is not None: + vbMng(verbObj, "INIT", + ("Starting optimization of {}x{} frontier " + "points.").format(len(front1), len(bulk1)), 65) + keepfront1, addfront1, addoptP = [], [], [] + for i, il1 in enumerate(front1): + il2 = optP[i] + change = False + for ib1 in bulk1: + ib2 = matches[ib1] + val = ( + distanceMatrixEff[il1, ib2] + distanceMatrixEff[ib1, il2] + - distanceMatrixEff[il1, il2] - distanceMatrixEff[ib1, ib2]) + if val < 0: + addfront1 += [ib1] + addoptP += [il2] + matches[ib1] = il2 + matches[il1] = ib2 + il2 = ib2 + optP[i] = ib2 + change = True + if change: keepfront1 += [i] + front1, optP = front1[keepfront1], optP[keepfront1] + for (addb, addo) in zip(addfront1[::-1], addoptP[::-1]): + if addb not in front1: + front1, optP = np.append(front1, addb), np.append(optP, addo) + if verbObj is not None: vbMng(verbObj, "DEL", "Done optimizing.", 65) + if len(front1) == 0: break + if verbObj is not None: vbMng(verbObj, "DEL", "Done point matching.", 25) + return idx2[matches] + +def pointMatching(distanceMatrix:Np2D, expl_threshold : int = 8, + max_iter : int = 10, verbObj = None, idx1 : Np1D = None, + idx2 : Np1D = None) -> Np1D: + N = len(distanceMatrix) + if idx1 is None: idx1 = np.arange(N) + if idx2 is None: idx2 = np.arange(N) + if len(idx1) > expl_threshold: + return pointMatchingHeuristic(distanceMatrix, idx1, idx2, max_iter, + verbObj) + return matchExplicit(distanceMatrix, idx1, idx2, verbObj) diff --git a/rrompy/utilities/base/sobol.py b/rrompy/utilities/numerical/sobol.py similarity index 100% rename from rrompy/utilities/base/sobol.py rename to rrompy/utilities/numerical/sobol.py diff --git a/rrompy/utilities/poly_fitting/__init__.py b/rrompy/utilities/poly_fitting/__init__.py index 0e297b8..7655844 100644 --- a/rrompy/utilities/poly_fitting/__init__.py +++ b/rrompy/utilities/poly_fitting/__init__.py @@ -1,29 +1,27 @@ # 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 . # from .custom_fit import customFit -from .custom_pinv import customPInv from .interpolator import GenericInterpolator __all__ = [ 'customFit', - 'customPInv', 'GenericInterpolator' ] diff --git a/rrompy/utilities/poly_fitting/heaviside/__init__.py b/rrompy/utilities/poly_fitting/heaviside/__init__.py index d815ffa..4074d47 100644 --- a/rrompy/utilities/poly_fitting/heaviside/__init__.py +++ b/rrompy/utilities/poly_fitting/heaviside/__init__.py @@ -1,40 +1,41 @@ # 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 . # from .base import polybases, polyfitname, polydomcoeff from .val import polyval from .vander import polyvander from .homogeneization import homogeneizedpolyvander from .heaviside_to_from_rational import heaviside2rational, rational2heaviside -from .heaviside_to_affine import heaviside2affine +from .heaviside_to_from_affine import heaviside2affine, affine2heaviside from .heaviside_interpolator import HeavisideInterpolator __all__ = [ 'polybases', 'polyfitname', 'polydomcoeff', 'polyval', 'polyvander', 'heaviside2rational', 'rational2heaviside', 'heaviside2affine', + 'affine2heaviside', 'homogeneizedpolyvander', 'HeavisideInterpolator' ] diff --git a/rrompy/utilities/poly_fitting/heaviside/heaviside_interpolator.py b/rrompy/utilities/poly_fitting/heaviside/heaviside_interpolator.py index 459b102..79334ef 100644 --- a/rrompy/utilities/poly_fitting/heaviside/heaviside_interpolator.py +++ b/rrompy/utilities/poly_fitting/heaviside/heaviside_interpolator.py @@ -1,77 +1,78 @@ # 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 . # import numpy as np from copy import deepcopy as copy from rrompy.utilities.base.types import (List, ListAny, Np1D, paramList, interpEng) +from rrompy.utilities.base import freepar as fp from rrompy.utilities.poly_fitting.polynomial.polynomial_interpolator import ( PolynomialInterpolator) from rrompy.utilities.poly_fitting.polynomial.roots import polyroots from rrompy.utilities.poly_fitting.heaviside.val import polyval from rrompy.utilities.poly_fitting.heaviside.heaviside_to_from_rational import\ heaviside2rational, rational2heaviside from rrompy.utilities.exception_manager import RROMPyAssert __all__ = ['HeavisideInterpolator'] class HeavisideInterpolator(PolynomialInterpolator): """HERE""" def __init__(self, other = None): if other is None: return self.poles = other.poles super().__init__(other) def __call__(self, mu:paramList, der : List[int] = None, scl : Np1D = None): - return polyval(mu, self.coeffs, self.polybasis, der, scl) + return polyval(mu, self.coeffs, self.poles, self.polybasis) def __copy__(self): return HeavisideInterpolator(self) def __deepcopy__(self, memo): other = HeavisideInterpolator() other.poles, other.coeffs, other.npar, other.polybasis = copy( (self.poles, self.coeffs, self.npar, self.polybasis), memo) return other def pad(self, nleft : List[int] = None, nright : List[int] = None): if nleft is None: nleft = [0] * len(self.shape) if nright is None: nright = [0] * len(self.shape) if not hasattr(nleft, "__len__"): nleft = [nleft] if not hasattr(nright, "__len__"): nright = [nright] RROMPyAssert(nleft[0], 0, "Padding in free direction") super().pad(nleft, nright) def setupFromRational(self, num:interpEng, den:interpEng, murange : Np1D = np.array([-1., 1.]), scl : Np1D = None, scaling : List[callable] = None, scalingInv : List[callable] = None): self.coeffs, self.poles, self.polybasis = rational2heaviside(num, den, murange, scl, scaling, scalingInv) - def roots(self, marginalVals : ListAny = [None], murange : Np1D = None, + def roots(self, marginalVals : ListAny = [fp], murange : Np1D = None, scaling : List[callable] = None, scalingInv : List[callable] = None): RROMPyAssert(self.shape, (1,), "Shape of output") - RROMPyAssert(marginalVals, [None], "Marginal values") + RROMPyAssert(marginalVals, [fp], "Marginal values") basisN = self.polybasis.split("_")[0] coeffsN = heaviside2rational(self.coeffs, self.poles, murange, basisN, scaling = scaling, scalingInv = scalingInv)[0] return polyroots(coeffsN, basisN) diff --git a/rrompy/utilities/poly_fitting/heaviside/heaviside_to_affine.py b/rrompy/utilities/poly_fitting/heaviside/heaviside_to_from_affine.py similarity index 51% rename from rrompy/utilities/poly_fitting/heaviside/heaviside_to_affine.py rename to rrompy/utilities/poly_fitting/heaviside/heaviside_to_from_affine.py index 2cb7029..13527a1 100644 --- a/rrompy/utilities/poly_fitting/heaviside/heaviside_to_affine.py +++ b/rrompy/utilities/poly_fitting/heaviside/heaviside_to_from_affine.py @@ -1,63 +1,94 @@ # 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 . # import numpy as np +from scipy.special import binom import scipy.sparse as sp -from rrompy.utilities.base.types import Np1D, Np2D, Tuple, paramVal, interpEng +from rrompy.utilities.base.types import Np1D, Np2D, List, ListAny, Tuple, paramVal +from rrompy.utilities.numerical import eigNonlinearDense from rrompy.utilities.exception_manager import RROMPyException from rrompy.parameter import checkParameter -__all__ = ['heaviside2affine'] +__all__ = ['heaviside2affine', 'affine2heaviside'] -def heaviside2affine(c:Np2D, poles:Np1D, murange : Np1D = None, - basis : str = "MONOMIAL_HEAVISIDE", mu : paramVal = [], - rescalingExp : int = 1, - sparse : bool = False) -> Tuple[interpEng, interpEng]: +def heaviside2affine(c:Np2D, poles:Np1D, mu : paramVal = [], + basis : str = "MONOMIAL_HEAVISIDE", + sparse : bool = False) \ + -> Tuple[Np2D, List[Np2D], List[Np1D]]: mu = checkParameter(mu, 1)(0, 0) n, d = len(poles), len(c) - len(poles) basisN = basis.split("_")[0] if basisN not in ["MONOMIAL", "CHEBYSHEV", "LEGENDRE"]: raise RROMPyException("Polynomial basis not recognized.") if sparse: - A0 = sp.spdiags([np.concatenate((mu - poles, np.ones(d)))], + A0 = sp.spdiags([np.concatenate((- mu - poles, np.ones(d)))], [0], n + d, n + d) A1 = sp.spdiags([np.concatenate((np.ones(n), np.zeros(d)))], [0], n + d, n + d) else: A0 = np.diag(np.concatenate((mu - poles, np.ones(d)))) A1 = np.diag(np.concatenate((np.ones(n), np.zeros(d)))) As = [A0, A1] bs = np.zeros((d, n + d), dtype = poles.dtype) bs[0, :] = 1. if d > 0: bs[0, n + 1 :] = 0. if d > 1: bs[1, n + 1] = 1. for j in range(2, d): if basisN == "MONOMIAL": bs[j, n + j] = 1. else: alpha = - 1. if basisN == "CHEBYSHEV" else 1. / j - 1. bs[:, n + j] = alpha * bs[:, n + j - 2] bs[1 :, n + j] += (1. - alpha) * bs[: -1, n + j - 1] bs = list(bs) - thetaAs = ["1.", "(mu ** ({}) - {})".format(rescalingExp, mu)] - thetabs = ["1."] - for j in range(1, d): - thetabs += ["((mu ** ({})- {}) ** ({}))".format(rescalingExp, mu, j)] - return c.reshape(c.shape[0], -1).T, As, bs, thetaAs, thetabs + return c.reshape(c.shape[0], -1).T, As, bs + +def affine2heaviside(As:ListAny, bs:ListAny, mu : paramVal = [], + jSupp : int = 1) -> Tuple[Np2D, Np1D, str]: + mu = checkParameter(mu, 1)(0, 0) + N = len(As) + M = len(bs) + n = As[0].shape[0] + if N == 1: + poles = np.empty(0, dtype = np.complex) + Q = np.eye(n) + else: + basis = "MONOMIAL_HEAVISIDE" + poles, P, Q = eigNonlinearDense(As, jSupp = jSupp, + return_inverse = True) + poles = poles + mu + P = P[- n :, :] + Q = Q[:, (jSupp - 1) * n : jSupp * n] + bEffs = np.array([Q.dot(np.linalg.solve(As[-1], b)) for b in bs]) + if N == 1: + c = bEffs + else: + c = np.zeros((len(poles) + M - 1, As[0].shape[1]), dtype = As[0].dtype) + for l, pl in enumerate(poles): + for i in range(M): + c[l, :] = pl ** i * bEffs[i, l] * P[:, l] + for l in range(M - 1): + for i in range(l + 1, M): + c[len(poles) + l, :] = P.dot(poles ** (i- 1 - l) * bEffs[i, :]) + for l in range(M - 1): + for i in range(l + 1, M - 1): + c[len(poles) + l, :] += (binom(i, l) * (- mu) ** (i - l) + * c[len(poles) + i, :]) + return c, poles, basis diff --git a/rrompy/utilities/poly_fitting/heaviside/heaviside_to_from_rational.py b/rrompy/utilities/poly_fitting/heaviside/heaviside_to_from_rational.py index 37104eb..697ac43 100644 --- a/rrompy/utilities/poly_fitting/heaviside/heaviside_to_from_rational.py +++ b/rrompy/utilities/poly_fitting/heaviside/heaviside_to_from_rational.py @@ -1,99 +1,99 @@ # 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 . # import numpy as np from rrompy.utilities.poly_fitting.polynomial.polynomial_interpolator import \ PolynomialInterpolator from rrompy.utilities.poly_fitting.polynomial.vander import polyvander from rrompy.utilities.poly_fitting.custom_fit import customFit from rrompy.utilities.base.types import Np1D, Np2D, List, Tuple, interpEng from rrompy.parameter.parameter_sampling import (RandomSampler as RS, QuadratureSampler as QS) from .val import polyval from rrompy.utilities.exception_manager import RROMPyException, RROMPyAssert __all__ = ['heaviside2rational', 'rational2heaviside'] def heaviside2rational(c:Np2D, poles:Np1D, murange : Np1D = None, basis : str = "MONOMIAL_HEAVISIDE", basisD : str = None, scaling : List[callable] = None, scalingInv : List[callable] = None) \ -> Tuple[interpEng, interpEng]: num = PolynomialInterpolator() den = PolynomialInterpolator() basisN = basis.split("_")[0] if basisD is None: basisD = basisN if murange is None: multiplier = [1., 1.j] avgs = [np.mean(np.real(poles)), np.mean(np.imag(poles))] ranges = np.array([[np.min(np.real(poles)), np.max(np.real(poles))], [np.min(np.imag(poles)), np.max(np.imag(poles))]]) domIdx = np.argmax(ranges[:, 1] - ranges[:, 0]) murange = [multiplier[domIdx] * x + multiplier[1 - domIdx] * avgs[1 - domIdx] for x in ranges[domIdx, :]] extraPt = None while extraPt is None or np.any(np.isclose(extraPt, poles)): extraPt = murange[0] + (murange[1] - murange[0]) * np.random.rand(1) denAuxPts = np.concatenate((poles, extraPt)) denAuxVals = np.concatenate((np.zeros(len(poles)), [1.])) den.setupByInterpolation(denAuxPts, denAuxVals, len(poles), basisD) den.coeffs /= np.linalg.norm(den.coeffs) if basis == "CHEBYSHEV": - sampler = QS(murange, "CHEBYSHEV") + sampler = QS(murange, "CHEBYSHEV", scaling, scalingInv) elif basis == "LEGENDRE": - sampler = QS(murange, "GAUSSLEGENDRE") + sampler = QS(murange, "GAUSSLEGENDRE", scaling, scalingInv) else: - sampler = RS(murange, "HALTON") + sampler = RS(murange, "HALTON", scaling, scalingInv) xAux = sampler.generatePoints(len(c)) valsAux = den(xAux) * polyval(xAux, c, poles, basis) num.setupByInterpolation(xAux, valsAux, len(c) - 1, basisN) return num, den def rational2heaviside(num:interpEng, den:interpEng, murange : Np1D = np.array([-1., 1.]), scl : Np1D = None, scaling : List[callable] = None, scalingInv : List[callable] = None) \ -> Tuple[Np2D, Np1D, str]: if (not isinstance(num, PolynomialInterpolator) or not isinstance(den, PolynomialInterpolator)): raise RROMPyException(("Rational numerator and denominator must be " "polynomial interpolators.")) RROMPyAssert(num.npar, 1, "Number of parameters") RROMPyAssert(den.npar, 1, "Number of parameters") basis = num.polybasis + "_HEAVISIDE" c = np.zeros_like(num.coeffs) poles = den.roots() Psp = num(poles) Qspder = den(poles, 1, scl) - c[: len(poles)] = Psp / Qspder + c[: len(poles)] = (Psp / Qspder).T if len(c) > len(poles): from rrompy.parameter.parameter_sampling import (RandomSampler as RS, QuadratureSampler as QS) if num.polybasis == "CHEBYSHEV": - sampler = QS(murange, "CHEBYSHEV") + sampler = QS(murange, "CHEBYSHEV", scaling, scalingInv) elif num.polybasis == "LEGENDRE": - sampler = QS(murange, "GAUSSLEGENDRE") + sampler = QS(murange, "GAUSSLEGENDRE", scaling, scalingInv) else: - sampler = RS(murange, "HALTON") + sampler = RS(murange, "HALTON", scaling, scalingInv) xAux = sampler.generatePoints(len(c)) valsAux = (num(xAux) / den(xAux) - - polyval(xAux, c, poles, basis)) + - polyval(xAux, c, poles, basis)).T VanAux = polyvander(xAux, [len(c) - len(poles) - 1], num.polybasis) c[len(poles) :] = customFit(VanAux, valsAux) return c, poles, basis diff --git a/rrompy/utilities/poly_fitting/heaviside/val.py b/rrompy/utilities/poly_fitting/heaviside/val.py index 5ee8433..2a8fa45 100644 --- a/rrompy/utilities/poly_fitting/heaviside/val.py +++ b/rrompy/utilities/poly_fitting/heaviside/val.py @@ -1,46 +1,46 @@ # 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 . # import numpy as np from copy import deepcopy as copy from rrompy.utilities.base.types import Np1D, Np2D, paramList from rrompy.parameter import checkParameterList from rrompy.utilities.poly_fitting.polynomial import polyval as pvP __all__ = ['polyval'] def polyval(x:paramList, c:Np2D, poles:Np1D, basis : str = "MONOMIAL_HEAVISIDE") -> Np2D: x = checkParameterList(x, 1)[0] poles = poles.flatten() if len(c) > len(poles): basisp = basis.split("_")[0] c0 = pvP(x, c[len(poles) :], basisp) else: - c0 = np.zeros((len(x),) + c.shape[1:], dtype = c.dtype) + c0 = np.zeros(c.shape[1:] + (len(x),), dtype = c.dtype) csh = copy(c0.shape) if len(csh) == 1: c0 = c0.reshape(1, -1) for j in range(len(x)): muDiff = 1. / (x[j] - poles) val = muDiff.dot(c[: len(poles)]) try: c0[..., j] += val except: c0[..., j] += val.flatten() if len(csh) == 1: c0 = c0.flatten() return c0 diff --git a/rrompy/utilities/poly_fitting/polynomial/marginalize.py b/rrompy/utilities/poly_fitting/polynomial/marginalize.py index a8d1702..67ce917 100644 --- a/rrompy/utilities/poly_fitting/polynomial/marginalize.py +++ b/rrompy/utilities/poly_fitting/polynomial/marginalize.py @@ -1,57 +1,58 @@ # 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 . # from numpy import array, polynomial as po from copy import deepcopy as copy from rrompy.utilities.base.types import Np1D, Np2D +from rrompy.utilities.base import freepar as fp from rrompy.utilities.exception_manager import RROMPyAssert, RROMPyException __all__ = ['polymarginalize'] -def polymarginalize(c:Np2D, basis:str, marginalVals : Np1D = [None], +def polymarginalize(c:Np2D, basis:str, marginalVals : Np1D = [fp], nMarginal : int = None) -> Np1D: if not hasattr(c, "ndim"): c = array(c) ndim = c.ndim if not hasattr(marginalVals, "__len__"): marginalVals = [marginalVals] marginalVals = list(marginalVals) try: polyvalbase = {"CHEBYSHEV" : po.chebyshev.chebval, "LEGENDRE" : po.legendre.legval, "MONOMIAL" : po.polynomial.polyval}[basis.upper()] except: raise RROMPyException("Polynomial basis not recognized.") RROMPyAssert(ndim, len(marginalVals), "Marginalized variables") marginalDims = [] for j in range(len(marginalVals)): - if marginalVals[j] is None: + if marginalVals[j] == fp: marginalDims += [c.shape[j]] if nMarginal is not None and len(marginalDims) != nMarginal: - raise RROMPyException(("Exactly {} 'None' entries in marginalVals " + raise RROMPyException(("Exactly {} 'freepar' entries in marginalVals " "must be provided.").format(nMarginal)) cEff = [copy(c)] for d in range(ndim): - if marginalVals[d] is not None: + if marginalVals[d] != fp: for dj in range(len(cEff)): cEff[dj] = polyvalbase(marginalVals[d], cEff[dj], tensor = False) else: cEff2 = [] for dj in range(len(cEff)): cEff2 += list(cEff[dj]) cEff = copy(cEff2) return array(cEff).reshape(tuple(marginalDims)) diff --git a/rrompy/utilities/poly_fitting/polynomial/polynomial_interpolator.py b/rrompy/utilities/poly_fitting/polynomial/polynomial_interpolator.py index 26fac46..90f8472 100644 --- a/rrompy/utilities/poly_fitting/polynomial/polynomial_interpolator.py +++ b/rrompy/utilities/poly_fitting/polynomial/polynomial_interpolator.py @@ -1,141 +1,141 @@ # 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 . # import numpy as np from copy import deepcopy as copy from rrompy.utilities.base.types import (List, ListAny, DictAny, Np1D, Np2D, paramList) +from rrompy.utilities.base import freepar as fp from rrompy.utilities.poly_fitting.interpolator import GenericInterpolator from rrompy.utilities.poly_fitting.custom_fit import customFit from rrompy.utilities.poly_fitting.polynomial.base import (polyfitname, polydomcoeff) from rrompy.utilities.poly_fitting.polynomial.val import polyval from rrompy.utilities.poly_fitting.polynomial.roots import polyroots from rrompy.utilities.poly_fitting.polynomial.homogeneization import ( homogeneizedpolyvander as hpv, homogeneizedToFull) from rrompy.utilities.poly_fitting.polynomial.vander import polyvander as pv from rrompy.utilities.exception_manager import RROMPyAssert, RROMPyException from rrompy.parameter import checkParameterList __all__ = ['PolynomialInterpolator'] class PolynomialInterpolator(GenericInterpolator): """HERE""" def __init__(self, other = None): if other is None: return self.coeffs = other.coeffs self.npar = other.npar self.polybasis = other.polybasis @property def shape(self): if self.coeffs.ndim > self.npar: sh = self.coeffs.shape[self.npar :] else: sh = tuple([1]) return sh @property def deg(self): return [x - 1 for x in self.coeffs.shape[: self.npar]] def __getitem__(self, key): return self.coeffs[key] def __call__(self, mu:paramList, der : List[int] = None, scl : Np1D = None): return polyval(mu, self.coeffs, self.polybasis, der, scl) def __copy__(self): return PolynomialInterpolator(self) def __deepcopy__(self, memo): other = PolynomialInterpolator() other.coeffs, other.npar, other.polybasis = copy( (self.coeffs, self.npar, self.polybasis), memo) return other @property def domCoeff(self): RROMPyAssert(self.npar, 1, "Number of parameters") return polydomcoeff(self.deg, self.polybasis) * self[-1] def pad(self, nleft : List[int] = None, nright : List[int] = None): if nleft is None: nleft = [0] * len(self.shape) if nright is None: nright = [0] * len(self.shape) if not hasattr(nleft, "__len__"): nleft = [nleft] if not hasattr(nright, "__len__"): nright = [nright] RROMPyAssert(len(self.shape), len(nleft), "Shape of output") RROMPyAssert(len(self.shape), len(nright), "Shape of output") padwidth = [(0, 0)] * self.npar padwidth = padwidth + [(l, r) for l, r in zip(nleft, nright)] self.coeffs = np.pad(self.coeffs, padwidth, "constant", constant_values = (0., 0.)) def postmultiplyTensorize(self, A:Np2D): RROMPyAssert(A.shape[0], self.shape[-1], "Shape of output") self.coeffs = np.tensordot(self.coeffs, A, axes = (-1, 0)) def setupByInterpolation(self, support:paramList, values:ListAny, deg:int, polybasis : str = "MONOMIAL", verbose : bool = True, homogeneized : bool = True, vanderCoeffs : DictAny = {}, fitCoeffs : DictAny = {}): support = checkParameterList(support)[0] self.npar = support.shape[1] self.polybasis = polybasis if homogeneized: vander, _, reorder = hpv(support, deg = deg, basis = polybasis, **vanderCoeffs) vander = vander[:, reorder] else: if not hasattr(deg, "__len__"): deg = [deg] * self.npar vander = pv(support, deg = deg, basis = polybasis, **vanderCoeffs) outDim = values.shape[1:] values = values.reshape(values.shape[0], -1) fitOut = customFit(vander, values, full = True, **fitCoeffs) P = fitOut[0] if verbose: msg = ("Fitting {} samples with degree {} through {}... " "Conditioning of LS system: {:.4e}.").format( len(vander), deg, polyfitname(self.polybasis), fitOut[1][2][0] / fitOut[1][2][-1]) else: msg = None if homogeneized: self.coeffs = homogeneizedToFull( tuple([deg + 1] * self.npar) + outDim, self.npar, P) else: self.coeffs = P.reshape(tuple([d + 1 for d in deg]) + outDim) return fitOut[1][1] == vander.shape[1], msg - def roots(self, marginalVals : ListAny = [None]): + def roots(self, marginalVals : ListAny = [fp]): RROMPyAssert(self.shape, (1,), "Shape of output") RROMPyAssert(len(marginalVals), self.npar, "Number of parameters") try: - rDim = marginalVals.index(None) - if (rDim < len(marginalVals) - 1 - and None in marginalVals[rDim + 1 :]): + rDim = marginalVals.index(fp) + if rDim < len(marginalVals) - 1 and fp in marginalVals[rDim + 1 :]: raise except: - raise RROMPyException(("Exactly 1 'None' entry in " + raise RROMPyException(("Exactly 1 'freepar' entry in " "marginalVals must be provided.")) return polyroots(self.coeffs, self.polybasis, marginalVals) diff --git a/rrompy/utilities/poly_fitting/polynomial/roots.py b/rrompy/utilities/poly_fitting/polynomial/roots.py index 40563f9..f392eb7 100644 --- a/rrompy/utilities/poly_fitting/polynomial/roots.py +++ b/rrompy/utilities/poly_fitting/polynomial/roots.py @@ -1,33 +1,34 @@ # 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 . # from numpy import polynomial as po from rrompy.utilities.base.types import Np1D, Np2D +from rrompy.utilities.base import freepar as fp from .marginalize import polymarginalize from rrompy.utilities.exception_manager import RROMPyException __all__ = ['polyroots'] -def polyroots(c:Np2D, basis:str, marginalVals : Np1D = [None]) -> Np1D: +def polyroots(c:Np2D, basis:str, marginalVals : Np1D = [fp]) -> Np1D: try: rootsbase = {"CHEBYSHEV" : po.chebyshev.chebroots, "LEGENDRE" : po.legendre.legroots, "MONOMIAL" : po.polynomial.polyroots}[basis.upper()] except: raise RROMPyException("Polynomial basis not recognized.") return rootsbase(polymarginalize(c, basis, marginalVals, 1)) diff --git a/rrompy/utilities/poly_fitting/radial_basis/radial_basis_interpolator.py b/rrompy/utilities/poly_fitting/radial_basis/radial_basis_interpolator.py index 6946292..cd2584c 100644 --- a/rrompy/utilities/poly_fitting/radial_basis/radial_basis_interpolator.py +++ b/rrompy/utilities/poly_fitting/radial_basis/radial_basis_interpolator.py @@ -1,140 +1,142 @@ # 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 . # import numpy as np from copy import deepcopy as copy from rrompy.utilities.base.types import (List, ListAny, DictAny, Np1D, Np2D, paramList) from rrompy.utilities.poly_fitting.interpolator import GenericInterpolator from rrompy.utilities.poly_fitting.custom_fit import customFit from rrompy.utilities.poly_fitting.radial_basis.base import polyfitname from rrompy.utilities.poly_fitting.radial_basis.val import polyval from rrompy.utilities.poly_fitting.radial_basis.homogeneization import ( homogeneizedpolyvander as hpv) from rrompy.utilities.poly_fitting.radial_basis.vander import polyvander as pv from rrompy.utilities.poly_fitting.polynomial.homogeneization import ( homogeneizedToFull) from rrompy.utilities.exception_manager import RROMPyException, RROMPyAssert from rrompy.parameter import checkParameterList __all__ = ['RadialBasisInterpolator'] class RadialBasisInterpolator(GenericInterpolator): """HERE""" def __init__(self, other = None): if other is None: return self.support = other.support self.coeffsGlobal = other.coeffsGlobal self.coeffsLocal = other.coeffsLocal self.directionalWeights = other.directionalWeights self.npar = other.npar self.polybasis = other.polybasis @property def shape(self): sh = self.coeffsLocal.shape[1 :] if self.coeffsLocal.ndim > 1 else 1 return sh @property def deg(self): return [x - 1 for x in self.coeffsGlobal.shape[: self.npar]] def __call__(self, mu:paramList, der : List[int] = None, scl : Np1D = None): if der is not None and np.sum(der) > 0: raise RROMPyException(("Cannot take derivatives of radial basis " "function.")) return polyval(mu, self.coeffsGlobal, self.coeffsLocal, self.support, self.directionalWeights, self.polybasis) def __copy__(self): return RadialBasisInterpolator(self) def __deepcopy__(self, memo): other = RadialBasisInterpolator() (other.support, other.coeffsGlobal, other.coeffsLocal, other.directionalWeights, other.npar, other.polybasis) = copy( (self.support, self.coeffsGlobal, self.coeffsLocal, self.directionalWeights, self.npar, self.polybasis), memo) return other def postmultiplyTensorize(self, A:Np2D): RROMPyAssert(A.shape[0], self.shape[-1], "Shape of output") self.coeffsLocal = np.tensordot(self.coeffsLocal, A, axes = (-1, 0)) self.coeffsGlobal = np.tensordot(self.coeffsGlobal, A, axes = (-1, 0)) def pad(self, nleft : List[int] = None, nright : List[int] = None): if nleft is None: nleft = [0] * len(self.shape) if nright is None: nright = [0] * len(self.shape) if not hasattr(nleft, "__len__"): nleft = [nleft] if not hasattr(nright, "__len__"): nright = [nright] RROMPyAssert(len(self.shape), len(nleft), "Shape of output") RROMPyAssert(len(self.shape), len(nright), "Shape of output") padwidth = [(0, 0)] + [(l, r) for l, r in zip(nleft, nright)] self.coeffsLocal = np.pad(self.coeffsLocal, padwidth, "constant", constant_values = (0., 0.)) padwidth = [(0, 0)] * (self.npar - 1) + padwidth self.coeffsGlobal = np.pad(self.coeffsGlobal, padwidth, "constant", constant_values = (0., 0.)) def setupByInterpolation(self, support:paramList, values:ListAny, deg:int, polybasis : str = "MONOMIAL_GAUSSIAN", directionalWeights : Np1D = None, verbose : bool = True, homogeneized : bool = True, vanderCoeffs : DictAny = {}, fitCoeffs : DictAny = {}): support = checkParameterList(support)[0] self.support = copy(support) + if "reorder" in vanderCoeffs.keys(): + self.support = self.support[vanderCoeffs["reorder"]] self.npar = support.shape[1] if directionalWeights is None: directionalWeights = np.ones(self.npar) self.directionalWeights = directionalWeights self.polybasis = polybasis if homogeneized: vander, _, reorder = hpv(support, deg = deg, basis = polybasis, directionalWeights = self.directionalWeights, **vanderCoeffs) vander = vander[reorder] vander = vander[:, reorder] else: if not hasattr(deg, "__len__"): deg = [deg] * self.npar vander = pv(support, deg = deg, basis = polybasis, **vanderCoeffs) outDim = values.shape[1:] values = values.reshape(values.shape[0], -1) values = np.pad(values, ((0, len(vander) - len(values)), (0, 0)), "constant") fitOut = customFit(vander, values, full = True, **fitCoeffs) P = fitOut[0][len(support) :] if verbose: msg = ("Fitting {}+{} samples with degree {} through {}... " "Conditioning of LS system: {:.4e}.").format( len(support), len(vander) - len(support), deg, polyfitname(self.polybasis), fitOut[1][2][0] / fitOut[1][2][-1]) else: msg = None self.coeffsLocal = fitOut[0][: len(support)] if homogeneized: self.coeffsGlobal = homogeneizedToFull( tuple([deg + 1] * self.npar) + outDim, self.npar, P) else: self.coeffsGlobal = P.reshape(tuple([d + 1 for d in deg]) + outDim) return fitOut[1][1] == vander.shape[1], msg diff --git a/tests/utilities/basic_routines.py b/tests/utilities/basic_routines.py index dfb251e..55d620b 100644 --- a/tests/utilities/basic_routines.py +++ b/tests/utilities/basic_routines.py @@ -1,66 +1,67 @@ # 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 . # import pytest import numpy as np from rrompy.utilities.exception_manager import RROMPyException, RROMPyWarning from rrompy.utilities.base import (findDictStrKey, purgeDict, purgeList, - lowDiscrepancy, verbosityDepth) + verbosityDepth) +from rrompy.utilities.numerical import lowDiscrepancy def test_exception(): with pytest.raises(Exception): raise RROMPyException("This is an exception!") def test_warning(capsys): RROMPyWarning("This is a warning.") out, err = capsys.readouterr() assert "This is a warning." in out assert len(err) == 0 def test_dict_list(): dictBase = {str(x): str(x**2) for x in range(10)} assert findDictStrKey("-1", dictBase.keys()) is None assert findDictStrKey("5", dictBase.keys()) == "5" dictPurged = purgeDict(dictBase, [str(x) for x in range(4)], silent = True) assert len(dictPurged.keys()) == 4 listBase = list(dictBase.values()) listPurged = purgeList(listBase, [str(x**2) for x in range(4,6)], silent = True) assert len(listPurged) == 2 def test_vandercorput(): idxs = lowDiscrepancy(8) assert np.all(idxs == np.array([0,4,2,6,1,5,3,7])) def test_verbositydepth(capsys): verbosityDepth("INIT", "Start message", timestamp = False) out, err = capsys.readouterr() assert out == "┌Start message\n" assert len(err) == 0 verbosityDepth("MAIN", "Main message", end = "\n\n", timestamp = False) out, err = capsys.readouterr() assert out == "├Main message\n\n" assert len(err) == 0 verbosityDepth("INIT", "Start message2") verbosityDepth("DEL", "Delete message2") verbosityDepth("DEL", "Delete message") with pytest.raises(Exception): verbosityDepth("DEL", "Wrong deletion")