diff --git a/rrompy/hfengines/base/marginal_proxy_engine.py b/rrompy/hfengines/base/marginal_proxy_engine.py
index 29ac89b..f4fc5ae 100644
--- a/rrompy/hfengines/base/marginal_proxy_engine.py
+++ b/rrompy/hfengines/base/marginal_proxy_engine.py
@@ -1,380 +1,387 @@
# 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, List, paramVal,
- paramList, HFEng, sampList)
+from rrompy.utilities.base.types import (Np1D, Np2D, ScOp, List, ListAny,
+ 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] == fp])
@property
def fixedLocations(self):
return tuple([x for x in range(self.nparBase) if self.marg[x] != fp])
@property
def muFixed(self):
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):
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]:
+ def affineLinearSystemA(self, mu : paramVal = [],
+ sclF : ListAny = None) -> List[Np2D]:
"""
Assemble affine blocks of operator of linear system (just linear
blocks).
"""
+ if sclF is None: sclF = [1.] * self.npar
As = [None] * self.nAs
for j in range(self.nAs):
- As[j] = self.A(mu, hashI(j, self.npar))
+ derIdx = hashI(j, self.npar)
+ mult = np.prod([s ** d for (s, d) in zip(sclF, derIdx)])
+ As[j] = mult * self.A(mu, derIdx)
return As
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):
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 = [],
+ def affineLinearSystemb(self, mu : paramVal = [], sclF : ListAny = None,
homogeneized : bool = False) -> List[Np1D]:
"""
Assemble affine blocks of RHS of linear system (just linear blocks).
"""
+ if sclF is None: sclF = [1.] * self.npar
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)
+ derIdx = hashI(j, self.npar)
+ mult = np.prod([s ** d for (s, d) in zip(sclF, derIdx)])
+ bs[j] = mult * self.b(mu, derIdx, homogeneized)
return bs
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 46e0604..5effc6c 100644
--- a/rrompy/hfengines/base/matrix_engine_base.py
+++ b/rrompy/hfengines/base/matrix_engine_base.py
@@ -1,551 +1,517 @@
# 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,
+ ListAny, DictAny, paramVal, paramList,
sampList)
from rrompy.utilities.base import (purgeList, getNewFilename,
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]:
+ def affineLinearSystemA(self, mu : paramVal = [],
+ sclF : ListAny = None) -> List[Np2D]:
"""
Assemble affine blocks of operator of linear system (just linear
blocks).
"""
+ if sclF is None: sclF = [1.] * self.npar
As = [None] * self.nAs
for j in range(self.nAs):
- As[j] = self.A(mu, hashI(j, self.npar))
+ derIdx = hashI(j, self.npar)
+ mult = np.prod([s ** d for (s, d) in zip(sclF, derIdx)])
+ As[j] = mult * self.A(mu, derIdx)
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
@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 = [],
+ def affineLinearSystemb(self, mu : paramVal = [], sclF : ListAny = None,
homogeneized : bool = False) -> List[Np1D]:
"""
Assemble affine blocks of RHS of linear system (just linear blocks).
"""
+ if sclF is None: sclF = [1.] * self.npar
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)
+ derIdx = hashI(j, self.npar)
+ mult = np.prod([s ** d for (s, d) in zip(sclF, derIdx)])
+ bs[j] = mult * self.b(mu, derIdx, 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 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/reduction_methods/base/generic_approximant.py b/rrompy/reduction_methods/base/generic_approximant.py
index 858c97d..53ca38f 100644
--- a/rrompy/reduction_methods/base/generic_approximant.py
+++ b/rrompy/reduction_methods/base/generic_approximant.py
@@ -1,869 +1,868 @@
# 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
from itertools import product as iterprod
from copy import deepcopy as copy
from os import remove as osrm
from rrompy.sampling.standard import (SamplingEngineStandard,
SamplingEngineStandardPOD)
from rrompy.utilities.base.types import (Np1D, DictAny, HFEng, List, Tuple,
ListAny, strLst, paramVal, paramList,
sampList)
from rrompy.utilities.base import (purgeDict, verbosityManager as vbMng,
getNewFilename)
from rrompy.utilities.exception_manager import (RROMPyException, RROMPyAssert,
RROMPy_READY, RROMPy_FRAGILE)
from rrompy.utilities.base import pickleDump, pickleLoad
from rrompy.parameter import (emptyParameterList, checkParameter,
checkParameterList)
from rrompy.sampling import sampleList, emptySampleList
__all__ = ['GenericApproximant']
def addNormFieldToClass(self, fieldName):
def objFunc(self, mu:paramList, homogeneized : bool = False) -> Np1D:
uV = getattr(self.__class__, "get" + fieldName)(self, mu, homogeneized)
val = self.HFEngine.norm(uV)
return val
setattr(self.__class__, "norm" + fieldName, objFunc)
def addPlotFieldToClass(self, fieldName):
def objFunc(self, mu:paramList, *args, homogeneized : bool = False,
**kwargs):
uV = getattr(self.__class__, "get" + fieldName)(self, mu, homogeneized)
kwargsCopy = copy(kwargs)
for j, u in enumerate(uV):
if "name" in kwargs.keys():
kwargsCopy["name"] = kwargs["name"] + str(j)
self.HFEngine.plot(u, *args, **kwargs)
setattr(self.__class__, "plot" + fieldName, objFunc)
def addOutParaviewFieldToClass(self, fieldName):
def objFunc(self, mu:paramVal, *args, homogeneized : bool = False,
**kwargs):
if not hasattr(self.HFEngine, "outParaview"):
raise RROMPyException(("High fidelity engine cannot output to "
"Paraview."))
uV = getattr(self.__class__, "get" + fieldName)(self, mu, homogeneized)
kwargsCopy = copy(kwargs)
for j, u in enumerate(uV):
if "name" in kwargs.keys():
kwargsCopy["name"] = kwargs["name"] + str(j)
self.HFEngine.outParaview(u, *args, **kwargsCopy)
setattr(self.__class__, "outParaview" + fieldName, objFunc)
def addOutParaviewTimeDomainFieldToClass(self, fieldName):
def objFunc(self, mu:paramVal, *args,
homogeneized : bool = False, **kwargs):
if not hasattr(self.HFEngine, "outParaviewTimeDomain"):
raise RROMPyException(("High fidelity engine cannot output to "
"Paraview."))
uV = getattr(self.__class__, "get" + fieldName)(self, mu, homogeneized)
omega = args.pop(0) if len(args) > 0 else np.real(mu)
kwargsCopy = copy(kwargs)
for j, u in enumerate(uV):
if "name" in kwargs.keys():
kwargsCopy["name"] = kwargs["name"] + str(j)
self.HFEngine.outParaviewTimeDomain(u, omega, *args,
**kwargsCopy)
setattr(self.__class__, "outParaviewTimeDomain" + fieldName, objFunc)
class GenericApproximant:
"""
ABSTRACT
ROM 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.
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.
trainedModel: Trained model evaluator.
mu0: Default parameter.
homogeneized: Whether to homogeneize Dirichlet BCs.
approxParameters: Dictionary containing values for main parameters of
approximant. Recognized keys are in parameterList{Soft,Critical}.
parameterListSoft: Recognized keys of soft approximant parameters:
- 'POD': whether to compute POD of snapshots.
parameterListCritical: Recognized keys of critical approximant
parameters:
- 'S': total number of samples current approximant relies upon.
verbosity: Verbosity level.
POD: Whether to compute POD of snapshots.
S: Number of solution snapshots over which current approximant is
based upon.
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.
"""
__all__ += [ftype + dtype for ftype, dtype in iterprod(
["norm", "plot", "outParaview", "outParaviewTimeDomain"],
["HF", "RHS", "Approx", "Res", "Err"])]
def __init__(self, HFEngine:HFEng, mu0 : paramVal = None,
approxParameters : DictAny = {}, homogeneized : bool = False,
verbosity : int = 10, timestamp : bool = True):
self._preInit()
self._mode = RROMPy_READY
self.verbosity = verbosity
self.timestamp = timestamp
vbMng(self, "INIT",
"Initializing engine of type {}.".format(self.name()), 10)
self._HFEngine = HFEngine
self.trainedModel = None
self.lastSolvedHF = emptyParameterList()
self.uHF = emptySampleList()
self._addParametersToList(["POD"], [True], ["S"], [[1]])
if mu0 is None:
if hasattr(self.HFEngine, "mu0"):
self.mu0 = checkParameter(self.HFEngine.mu0)
else:
raise RROMPyException(("Center of approximation cannot be "
"inferred from HF engine. Parameter "
"required"))
else:
self.mu0 = checkParameter(mu0, self.HFEngine.npar)
self.resetSamples()
self.homogeneized = homogeneized
self.approxParameters = approxParameters
self._postInit()
### add norm{HF,RHS,Approx,Res,Err} methods
"""
Compute norm of * at arbitrary parameter.
Args:
mu: Target parameter.
homogeneized(optional): Whether to remove Dirichlet BC. Defaults to
False.
Returns:
Target norm of *.
"""
for objName in ["HF", "RHS", "Err"]:
addNormFieldToClass(self, objName)
def objFunc(self, mu:paramList, homogeneized : bool = False) -> Np1D:
# uV = getattr(self.__class__, "getRes")(self, mu, homogeneized,
# duality = False)
uV = self.getRes(mu, homogeneized, duality = False)
val = self.HFEngine.norm(uV, dual = True, duality = False)
return val
setattr(self.__class__, "normRes", objFunc)
if not hasattr(self, "normApprox"):
addNormFieldToClass(self, "Approx")
### add plot{HF,RHS,Approx,Res,Err} methods
"""
Do some nice plots of * at arbitrary parameter.
Args:
mu: Target parameter.
name(optional): Name to be shown as title of the plots. Defaults to
'u'.
what(optional): Which plots to do. If list, can contain 'ABS',
'PHASE', 'REAL', 'IMAG'. If str, same plus wildcard 'ALL'.
Defaults to 'ALL'.
save(optional): Where to save plot(s). Defaults to None, i.e. no
saving.
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.
homogeneized(optional): Whether to remove Dirichlet BC. Defaults to
False.
figspecs(optional key args): Optional arguments for matplotlib
figure creation.
"""
for objName in ["HF", "RHS", "Approx", "Res", "Err"]:
addPlotFieldToClass(self, objName)
### add outParaview{HF,RHS,Approx,Res,Err} methods
"""
Output * to ParaView file.
Args:
mu: Target parameter.
name(optional): Base name to be used for data output.
filename(optional): Name of output file.
time(optional): Timestamp.
what(optional): Which plots to do. If list, can contain 'MESH',
'ABS', 'PHASE', 'REAL', 'IMAG'. If str, same plus wildcard
'ALL'. Defaults to 'ALL'.
forceNewFile(optional): Whether to create new output file.
filePW(optional): Fenics File entity (for time series).
homogeneized(optional): Whether to remove Dirichlet BC. Defaults to
False.
"""
for objName in ["HF", "RHS", "Approx", "Res", "Err"]:
addOutParaviewFieldToClass(self, objName)
### add outParaviewTimeDomain{HF,RHS,Approx,Res,Err} methods
"""
Output * to ParaView file, converted to time domain.
Args:
mu: Target parameter.
omega(optional): frequency.
timeFinal(optional): final time of simulation.
periodResolution(optional): number of time steps per period.
name(optional): Base name to be used for data output.
filename(optional): Name of output file.
forceNewFile(optional): Whether to create new output file.
homogeneized(optional): Whether to remove Dirichlet BC. Defaults to
False.
"""
for objName in ["HF", "RHS", "Approx", "Res", "Err"]:
addOutParaviewTimeDomainFieldToClass(self, objName)
def _preInit(self):
if not hasattr(self, "depth"): self.depth = 0
else: self.depth += 1
@property
def parameterList(self):
"""Value of parameterListSoft + parameterListCritical."""
return self.parameterListSoft + self.parameterListCritical
def _addParametersToList(self, whatSoft:strLst, defaultSoft:ListAny,
whatCritical : strLst = [],
defaultCritical : ListAny = [],
toBeExcluded : strLst = []):
if not hasattr(self, "parameterToBeExcluded"):
self.parameterToBeExcluded = []
self.parameterToBeExcluded += toBeExcluded
if not hasattr(self, "parameterListSoft"):
self.parameterListSoft = []
if not hasattr(self, "parameterDefaultSoft"):
self.parameterDefaultSoft = {}
if not hasattr(self, "parameterListCritical"):
self.parameterListCritical = []
if not hasattr(self, "parameterDefaultCritical"):
self.parameterDefaultCritical = {}
for j, what in enumerate(whatSoft):
if what not in self.parameterToBeExcluded:
self.parameterListSoft += [what]
self.parameterDefaultSoft[what] = defaultSoft[j]
for j, what in enumerate(whatCritical):
if what not in self.parameterToBeExcluded:
self.parameterListCritical += [what]
self.parameterDefaultCritical[what] = defaultCritical[j]
def _postInit(self):
if self.depth == 0:
vbMng(self, "DEL", "Done initializing.", 10)
del self.depth
else: self.depth -= 1
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 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 = SamplingEngineStandardPOD
else:
SamplingEngine = SamplingEngineStandard
self.samplingEngine = SamplingEngine(self.HFEngine,
verbosity = self.verbosity,
allowRepeatedSamples = True)
@property
def HFEngine(self):
"""Value of HFEngine."""
return self._HFEngine
@HFEngine.setter
def HFEngine(self, HFEngine):
raise RROMPyException("Cannot change HFEngine.")
@property
def mu0(self):
"""Value of mu0."""
return self._mu0
@mu0.setter
def mu0(self, mu0):
mu0 = checkParameter(mu0)
if not hasattr(self, "_mu0") or mu0 != self.mu0:
self.resetSamples()
self._mu0 = mu0
@property
def npar(self):
"""Number of parameters."""
return self.mu0.shape[1]
@property
def approxParameters(self):
"""Value of approximant parameters."""
return self._approxParameters
@approxParameters.setter
def approxParameters(self, approxParams):
if not hasattr(self, "approxParameters"):
self._approxParameters = {}
approxParameters = purgeDict(approxParams, self.parameterList,
dictname = self.name() + ".approxParameters",
baselevel = 1)
keyList = list(approxParameters.keys())
for key in self.parameterListCritical:
if key in keyList:
setattr(self, "_" + key, self.parameterDefaultCritical[key])
for key in self.parameterListSoft:
if key in keyList:
setattr(self, "_" + key, self.parameterDefaultSoft[key])
fragile = False
for key in self.parameterListCritical:
if key in keyList:
val = approxParameters[key]
else:
val = getattr(self, "_" + key, None)
if val is None:
val = self.parameterDefaultCritical[key]
getattr(self.__class__, key, None).fset(self, val)
fragile = fragile or val is None
for key in self.parameterListSoft:
if key in keyList:
val = approxParameters[key]
else:
val = getattr(self, "_" + key, None)
if val is None:
val = self.parameterDefaultSoft[key]
getattr(self.__class__, key, None).fset(self, val)
if fragile:
self._mode = RROMPy_FRAGILE
@property
def POD(self):
"""Value of POD."""
return self._POD
@POD.setter
def POD(self, POD):
if hasattr(self, "_POD"): PODold = self.POD
else: PODold = -1
self._POD = POD
self._approxParameters["POD"] = self.POD
if PODold != self.POD:
self.samplingEngine = None
self.resetSamples()
@property
def S(self):
"""Value of S."""
return self._S
@S.setter
def S(self, S):
if not hasattr(S, "__len__"): S = [S]
if any([s <= 0 for s in S]):
raise RROMPyException("S must be positive.")
if hasattr(self, "_S") and self._S is not None: Sold = tuple(self.S)
else: Sold = -1
self._S = S
self._approxParameters["S"] = self.S
if Sold != tuple(self.S):
self.resetSamples()
@property
def homogeneized(self):
"""Value of homogeneized."""
return self._homogeneized
@homogeneized.setter
def homogeneized(self, homogeneized):
if not hasattr(self, "_homogeneized"):
self._homogeneized = None
if homogeneized != self.homogeneized:
self._homogeneized = homogeneized
self.resetSamples()
@property
def trainedModel(self):
"""Value of trainedModel."""
return self._trainedModel
@trainedModel.setter
def trainedModel(self, trainedModel):
self._trainedModel = trainedModel
if self._trainedModel is not None:
self._trainedModel.lastSolvedApproxReduced = emptyParameterList()
self._trainedModel.lastSolvedApprox = emptyParameterList()
self.lastSolvedApproxReduced = emptyParameterList()
self.lastSolvedApprox = emptyParameterList()
self.uApproxReduced = emptySampleList()
self.uApprox = emptySampleList()
def resetSamples(self):
if hasattr(self, "samplingEngine") and self.samplingEngine is not None:
self.samplingEngine.resetHistory()
else:
self.setupSampling()
self._mode = RROMPy_READY
def plotSamples(self, warping : List[callable] = None, name : str = "u",
save : str = None, what : strLst = 'all',
saveFormat : str = "eps", saveDPI : int = 100,
show : bool = True, plotArgs : dict = {}, **figspecs):
"""
Do some nice plots of the samples.
Args:
warping(optional): Domain warping functions.
name(optional): Name to be shown as title of the plots. Defaults to
'u'.
what(optional): Which plots to do. If list, can contain 'ABS',
'PHASE', 'REAL', 'IMAG'. If str, same plus wildcard 'ALL'.
Defaults to 'ALL'.
save(optional): Where to save plot(s). Defaults to None, i.e. no
saving.
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.
plotArgs(optional): Optional arguments for fen/pyplot.
figspecs(optional key args): Optional arguments for matplotlib
figure creation.
"""
RROMPyAssert(self._mode, message = "Cannot plot samples.")
self.samplingEngine.plotSamples(warping, name, save, what, saveFormat,
saveDPI, show, plotArgs, **figspecs)
def outParaviewSamples(self, name : str = "u", filename : str = "out",
times : Np1D = None, what : strLst = 'all',
forceNewFile : bool = True, folders : bool = False,
filePW = None):
"""
Output samples to ParaView file.
Args:
name(optional): Base name to be used for data output.
filename(optional): Name of output file.
times(optional): Timestamps.
what(optional): Which plots to do. If list, can contain 'MESH',
'ABS', 'PHASE', 'REAL', 'IMAG'. If str, same plus wildcard
'ALL'. Defaults to 'ALL'.
forceNewFile(optional): Whether to create new output file.
folders(optional): Whether to split output in folders.
filePW(optional): Fenics File entity (for time series).
"""
RROMPyAssert(self._mode, message = "Cannot output samples.")
self.samplingEngine.outParaviewSamples(name = name,
filename = filename,
times = times, what = what,
forceNewFile = forceNewFile,
folders = folders,
filePW = filePW)
def outParaviewTimeDomainSamples(self, omegas : Np1D = None,
timeFinal : Np1D = None,
periodResolution : int = 20,
name : str = "u",
filename : str = "out",
forceNewFile : bool = True,
folders : bool = False):
"""
Output samples to ParaView file, converted to time domain.
Args:
omegas(optional): frequencies.
timeFinal(optional): final time of simulation.
periodResolution(optional): number of time steps per period.
name(optional): Base name to be used for data output.
filename(optional): Name of output file.
forceNewFile(optional): Whether to create new output file.
folders(optional): Whether to split output in folders.
"""
RROMPyAssert(self._mode, message = "Cannot output samples.")
self.samplingEngine.outParaviewTimeDomainSamples(omegas = omegas,
timeFinal = timeFinal,
periodResolution = periodResolution,
name = name, filename = filename,
forceNewFile = forceNewFile,
folders = folders)
def setSamples(self, samplingEngine):
"""Copy samplingEngine and samples."""
vbMng(self, "INIT", "Transfering samples.", 10)
self.samplingEngine = copy(samplingEngine)
vbMng(self, "DEL", "Done transfering samples.", 10)
def setTrainedModel(self, model):
"""Deepcopy approximation from trained model."""
if hasattr(model, "storeTrainedModel"):
verb = model.verbosity
model.verbosity = 0
fileOut = model.storeTrainedModel()
model.verbosity = verb
else:
try:
fileOut = getNewFilename("trained_model", "pkl")
pickleDump(model.data.__dict__, fileOut)
except:
raise RROMPyException(("Failed to store model data. Parameter "
"model must have either "
"storeTrainedModel or "
"data.__dict__ properties."))
self.loadTrainedModel(fileOut)
osrm(fileOut)
@abstractmethod
def setupApprox(self):
"""
Setup approximant. (ABSTRACT)
Any specialization should include something like
if self.checkComputedApprox():
return
RROMPyAssert(self._mode, message = "Cannot setup approximant.")
...
self.trainedModel = ...
self.trainedModel.data = ...
self.trainedModel.data.approxParameters = copy(
self.approxParameters)
"""
pass
def checkComputedApprox(self) -> bool:
"""
Check if setup of new approximant is not needed.
Returns:
True if new setup is not needed. False otherwise.
"""
return self._mode == RROMPy_FRAGILE or (self.trainedModel is not None
and self.trainedModel.data.approxParameters == self.approxParameters)
def _pruneBeforeEval(self, mu:paramList, field:str, append:bool,
prune:bool) -> Tuple[paramList, Np1D, Np1D, bool]:
mu = checkParameterList(mu, self.npar)[0]
idx = np.empty(len(mu), dtype = np.int)
if prune:
jExtra = np.zeros(len(mu), dtype = bool)
muExtra = emptyParameterList()
lastSolvedMus = getattr(self, "lastSolved" + field)
if (len(mu) > 0 and len(mu) == len(lastSolvedMus)
and mu == lastSolvedMus):
idx = np.arange(len(mu), dtype = np.int)
return muExtra, jExtra, idx, True
muKeep = copy(muExtra)
for j in range(len(mu)):
jPos = lastSolvedMus.find(mu[j])
if jPos is not None:
idx[j] = jPos
muKeep.append(mu[j])
else:
jExtra[j] = True
muExtra.append(mu[j])
if len(muKeep) > 0 and not append:
lastSolvedu = getattr(self, "u" + field)
idx[~jExtra] = getattr(self.__class__, "set" + field)(self,
muKeep, lastSolvedu[idx[~jExtra]], append)
append = True
else:
jExtra = np.ones(len(mu), dtype = bool)
muExtra = mu
return muExtra, jExtra, idx, append
def _setObject(self, mu:paramList, field:str, object:sampList,
append:bool) -> List[int]:
newMus = checkParameterList(mu, self.npar)[0]
newObj = sampleList(object)
if append:
getattr(self, "lastSolved" + field).append(newMus)
getattr(self, "u" + field).append(newObj)
Ltot = len(getattr(self, "u" + field))
return list(range(Ltot - len(newObj), Ltot))
setattr(self, "lastSolved" + field, copy(newMus))
setattr(self, "u" + field, copy(newObj))
return list(range(len(getattr(self, "u" + field))))
def setHF(self, muHF:paramList, uHF:sampleList,
append : bool = False) -> List[int]:
"""Assign high fidelity solution."""
return self._setObject(muHF, "HF", uHF, append)
def evalHF(self, mu:paramList, append : bool = False,
prune : bool = True) -> List[int]:
"""
Find high fidelity solution with original parameters and arbitrary
parameter.
Args:
mu: Target parameter.
append(optional): Whether to append new HF solutions to old ones.
prune(optional): Whether to remove duplicates of already appearing
HF solutions.
"""
muExtra, jExtra, idx, append = self._pruneBeforeEval(mu, "HF", append,
prune)
if len(muExtra) > 0:
vbMng(self, "INIT", "Solving HF model for mu = {}.".format(mu),
15)
newuHFs = self.HFEngine.solve(muExtra,
homogeneized = self.homogeneized)
vbMng(self, "DEL", "Done solving HF model.", 15)
idx[jExtra] = self.setHF(muExtra, newuHFs, append)
return list(idx)
def setApproxReduced(self, muApproxR:paramList, uApproxR:sampleList,
append : bool = False) -> List[int]:
"""Assign high fidelity solution."""
return self._setObject(muApproxR, "ApproxReduced", uApproxR, append)
def evalApproxReduced(self, mu:paramList, append : bool = False,
prune : bool = True) -> List[int]:
"""
Evaluate reduced representation of approximant at arbitrary parameter.
Args:
mu: Target parameter.
append(optional): Whether to append new HF solutions to old ones.
prune(optional): Whether to remove duplicates of already appearing
HF solutions.
"""
self.setupApprox()
muExtra, jExtra, idx, append = self._pruneBeforeEval(mu,
"ApproxReduced",
append, prune)
if len(muExtra) > 0:
newuApproxs = self.trainedModel.getApproxReduced(muExtra)
idx[jExtra] = self.setApproxReduced(muExtra, newuApproxs, append)
return list(idx)
def setApprox(self, muApprox:paramList, uApprox:sampleList,
append : bool = False) -> List[int]:
"""Assign high fidelity solution."""
return self._setObject(muApprox, "Approx", uApprox, append)
def evalApprox(self, mu:paramList, append : bool = False,
prune : bool = True) -> List[int]:
"""
Evaluate approximant at arbitrary parameter.
Args:
mu: Target parameter.
append(optional): Whether to append new HF solutions to old ones.
prune(optional): Whether to remove duplicates of already appearing
HF solutions.
"""
self.setupApprox()
muExtra, jExtra, idx, append = self._pruneBeforeEval(mu, "Approx",
append, prune)
if len(muExtra) > 0:
newuApproxs = self.trainedModel.getApprox(muExtra)
idx[jExtra] = self.setApprox(muExtra, newuApproxs, append)
return list(idx)
def getHF(self, mu:paramList, homogeneized : bool = False,
append : bool = False, prune : bool = True) -> sampList:
"""
Get HF solution at arbitrary parameter.
Args:
mu: Target parameter.
homogeneized(optional): Whether to remove Dirichlet BC. Defaults to
False.
Returns:
HFsolution.
"""
mu = checkParameterList(mu, self.npar)[0]
idx = self.evalHF(mu, append = append, prune = prune)
uHFs = self.uHF(idx)
if self.homogeneized and not homogeneized:
for j, m in enumerate(mu):
uHFs[j] += self.HFEngine.liftDirichletData(m)
if not self.homogeneized and homogeneized:
for j, m in enumerate(mu):
uHFs[j] -= self.HFEngine.liftDirichletData(m)
return uHFs
def getRHS(self, mu:paramList, homogeneized : bool = False,
duality : bool = True) -> sampList:
"""
Get linear system RHS at arbitrary parameter.
Args:
mu: Target parameter.
homogeneized(optional): Whether to remove Dirichlet BC. Defaults to
False.
Returns:
Linear system RHS.
"""
return self.HFEngine.residual(None, mu, homogeneized = homogeneized,
duality = duality)
def getApproxReduced(self, mu:paramList, append : bool = False,
prune : bool = True) -> sampList:
"""
Get approximant at arbitrary parameter.
Args:
mu: Target parameter.
Returns:
Reduced approximant.
"""
mu = checkParameterList(mu, self.npar)[0]
idx = self.evalApproxReduced(mu, append = append, prune = prune)
uApproxRs = self.uApproxReduced(idx)
return uApproxRs
def getApprox(self, mu:paramList, homogeneized : bool = False,
append : bool = False, prune : bool = True) -> sampList:
"""
Get approximant at arbitrary parameter.
Args:
mu: Target parameter.
homogeneized(optional): Whether to remove Dirichlet BC. Defaults to
False.
Returns:
Approximant.
"""
mu = checkParameterList(mu, self.npar)[0]
idx = self.evalApprox(mu, append = append, prune = prune)
uApproxs = self.uApprox(idx)
if self.homogeneized and not homogeneized:
for j, m in enumerate(mu):
uApproxs[j] += self.HFEngine.liftDirichletData(m)
if not self.homogeneized and homogeneized:
for j, m in enumerate(mu):
uApproxs[j] -= self.HFEngine.liftDirichletData(m)
return uApproxs
def getRes(self, mu:paramList, homogeneized : bool = False,
duality : bool = True) -> sampList:
"""
Get residual at arbitrary parameter.
Args:
mu: Target parameter.
homogeneized(optional): Whether to remove Dirichlet BC. Defaults to
False.
Returns:
Approximant residual.
"""
return self.HFEngine.residual(self.getApprox(mu, homogeneized), mu,
homogeneized = homogeneized,
duality = duality)
def getErr(self, mu:paramList, homogeneized : bool = False,
append : bool = False, prune : bool = True) -> sampList:
"""
Get error at arbitrary parameter.
Args:
mu: Target parameter.
homogeneized(optional): Whether to remove Dirichlet BC. Defaults to
False.
Returns:
Approximant error.
"""
return (self.getApprox(mu, homogeneized, append = append, prune =prune)
- self.getHF(mu, homogeneized, append = append, prune = prune))
def getPoles(self, *args, **kwargs) -> Np1D:
"""
Obtain approximant poles.
Returns:
Numpy complex vector of poles.
"""
self.setupApprox()
vbMng(self, "INIT", "Computing poles of model.", 20)
poles = self.trainedModel.getPoles(*args, **kwargs)
vbMng(self, "DEL", "Done computing poles.", 20)
return poles
def storeTrainedModel(self, filenameBase : str = "trained_model",
forceNewFile : bool = True) -> str:
"""Store trained reduced model to file."""
self.setupApprox()
vbMng(self, "INIT", "Storing trained model to file.", 20)
if forceNewFile:
filename = getNewFilename(filenameBase, "pkl")
else:
filename = "{}.pkl".format(filenameBase)
pickleDump(self.trainedModel.data.__dict__, filename)
vbMng(self, "DEL", "Done storing trained model.", 20)
return filename
def loadTrainedModel(self, filename:str):
"""Load trained reduced model from file."""
vbMng(self, "INIT", "Loading pre-trained model from file.", 20)
datadict = pickleLoad(filename)
name = datadict.pop("name")
if name == "TrainedModelRational":
from rrompy.reduction_methods.trained_model import \
TrainedModelRational as tModel
elif name == "TrainedModelRB":
from rrompy.reduction_methods.trained_model import \
TrainedModelRB as tModel
else:
raise RROMPyException(("Trained model name not recognized. "
"Loading failed."))
self.mu0 = datadict.pop("mu0")
from rrompy.reduction_methods.trained_model import TrainedModelData
trainedModel = tModel()
trainedModel.verbosity = self.verbosity
trainedModel.timestamp = self.timestamp
+ self.scaleFactor = datadict.pop("scaleFactor")
data = TrainedModelData(name, self.mu0, datadict.pop("projMat"),
- datadict.pop("rescalingExp"))
+ self.scaleFactor, datadict.pop("rescalingExp"))
if "mus" in datadict:
data.mus = datadict.pop("mus")
approxParameters = datadict.pop("approxParameters")
data.approxParameters = copy(approxParameters)
if "sampler" in approxParameters:
self._approxParameters["sampler"] = approxParameters.pop("sampler")
self.approxParameters = copy(approxParameters)
if "mus" in data.__dict__:
self.mus = copy(data.mus)
- self.scaleFactor = datadict.pop("scaleFactor")
- data.scaleFactor = self.scaleFactor
for key in datadict:
setattr(data, key, datadict[key])
trainedModel.data = data
self.trainedModel = trainedModel
self._mode = RROMPy_FRAGILE
vbMng(self, "DEL", "Done loading pre-trained model.", 20)
diff --git a/rrompy/reduction_methods/base/rb_utils.py b/rrompy/reduction_methods/base/rb_utils.py
index 78b37cc..fde59e0 100644
--- a/rrompy/reduction_methods/base/rb_utils.py
+++ b/rrompy/reduction_methods/base/rb_utils.py
@@ -1,62 +1,68 @@
# 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 Np1D, Np2D, Tuple, List, sampList
+from rrompy.utilities.base.types import (Np1D, Np2D, Tuple, List, ListAny,
+ sampList)
+from rrompy.utilities.poly_fitting.polynomial import (
+ hashIdxToDerivative as hashI)
from rrompy.utilities.exception_manager import RROMPyAssert
from rrompy.sampling import sampleList
__all__ = ['projectAffineDecomposition']
def projectAffineDecomposition(As:List[Np2D], bs:List[Np1D], pMat:sampList,
ARBsOld : List[Np2D] = None,
bRBsOld : List[Np1D] = None,
pMatOld : sampList = None)\
-> Tuple[List[Np2D], List[Np1D]]:
"""Project affine decomposition of linear system onto basis."""
RROMPyAssert((ARBsOld is None, bRBsOld is None),
(pMatOld is None, pMatOld is None),
"Old affine projected terms")
if isinstance(pMat, (sampleList,)): pMat = pMat.data
pMatH = pMat.T.conj()
ARBs = [None] * len(As)
bRBs = [None] * len(bs)
if pMatOld is None:
- for j in range(len(As)):
- ARBs[j] = pMatH.dot(As[j].dot(pMat))
- for j in range(len(bs)):
- bRBs[j] = pMatH.dot(bs[j])
+ for j in range(max(len(As), len(bs))):
+ if j < len(As):
+ ARBs[j] = pMatH.dot(As[j].dot(pMat))
+ if j < len(bs):
+ bRBs[j] = pMatH.dot(bs[j])
else:
RROMPyAssert((len(ARBsOld), len(bRBsOld)), (len(As), len(bs)),
"Old affine projected terms")
if isinstance(pMatOld, (sampleList,)): pMatOld = pMatOld.data
pMatOldH = pMatOld.T.conj()
Sold = pMatOld.shape[1]
Snew = pMat.shape[1]
- for j in range(len(As)):
- ARBs[j] = np.empty((Sold + Snew, Sold + Snew), dtype = np.complex)
- ARBs[j][: Sold, : Sold] = ARBsOld[j]
- ARBs[j][: Sold, Sold :] = pMatOldH.dot(As[j].dot(pMat))
- ARBs[j][Sold :, : Sold] = pMatH.dot(As[j].dot(pMatOld))
- ARBs[j][Sold :, Sold :] = pMatH.dot(As[j].dot(pMat))
- for j in range(len(bs)):
- bRBs[j] = np.empty((Sold + Snew), dtype = np.complex)
- bRBs[j][: Sold] = bRBsOld[j]
- bRBs[j][Sold :] = pMatH.dot(bs[j])
+ for j in range(max(len(As), len(bs))):
+ if j < len(As):
+ ARBs[j] = np.empty((Sold + Snew, Sold + Snew),
+ dtype = np.complex)
+ ARBs[j][: Sold, : Sold] = ARBsOld[j]
+ ARBs[j][: Sold, Sold :] = pMatOldH.dot(As[j].dot(pMat))
+ ARBs[j][Sold :, : Sold] = pMatH.dot(As[j].dot(pMatOld))
+ ARBs[j][Sold :, Sold :] = pMatH.dot(As[j].dot(pMat))
+ if j < len(bs):
+ bRBs[j] = np.empty((Sold + Snew), dtype = np.complex)
+ bRBs[j][: Sold] = bRBsOld[j]
+ bRBs[j][Sold :] = pMatH.dot(bs[j])
return ARBs, bRBs
diff --git a/rrompy/reduction_methods/distributed/rational_interpolant.py b/rrompy/reduction_methods/distributed/rational_interpolant.py
index 7d69bf0..7439332 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.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
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.scaleFactor,
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 ea97567..3f142ed 100644
--- a/rrompy/reduction_methods/distributed/rb_distributed.py
+++ b/rrompy/reduction_methods/distributed/rb_distributed.py
@@ -1,203 +1,207 @@
# 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.types import (Np1D, Np2D, List, ListAny, 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.
bs: List of numpy vectors representing coefficients of linear system
RHS.
ARBs: List of sparse matrices (in CSC format) representing coefficients
of compressed linear system matrix.
bRBs: List of numpy vectors representing coefficients of compressed
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.computeScaleFactor()
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)
+ pMat, self.scaleFactor,
+ self.HFEngine.rescalingExp)
data.mus = copy(self.mus)
self.trainedModel.data = data
else:
self.trainedModel = self.trainedModel
self.trainedModel.data.projMat = copy(pMat)
+ vbMng(self, "INIT", "Computing affine blocks of system.", 10)
+ self.As = self.HFEngine.affineLinearSystemA(self.mu0, self.scaleFactor)
+ self.bs = self.HFEngine.affineLinearSystemb(self.mu0, self.scaleFactor,
+ self.homogeneized)
+ vbMng(self, "DEL", "Done computing affine blocks.", 10)
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/rational_interpolant_greedy.py b/rrompy/reduction_methods/distributed_greedy/rational_interpolant_greedy.py
index 9b47c7a..b1f06b1 100644
--- a/rrompy/reduction_methods/distributed_greedy/rational_interpolant_greedy.py
+++ b/rrompy/reduction_methods/distributed_greedy/rational_interpolant_greedy.py
@@ -1,412 +1,412 @@
# 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_greedy_approximant import \
GenericDistributedGreedyApproximant
from rrompy.utilities.poly_fitting.polynomial import polybases
from rrompy.reduction_methods.distributed import RationalInterpolant
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, DictAny, HFEng, paramVal
from rrompy.utilities.base import verbosityManager as vbMng
from rrompy.utilities.exception_manager import RROMPyWarning
from rrompy.utilities.exception_manager import RROMPyException, RROMPyAssert
__all__ = ['RationalInterpolantGreedy']
class RationalInterpolantGreedy(GenericDistributedGreedyApproximant,
RationalInterpolant):
"""
ROM greedy 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': number of starting training points;
- 'sampler': sample point generator;
- '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;
- '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 training points to be exhausted
before training 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;
- 'polybasis': type of basis for interpolation; allowed values
include 'MONOMIAL', 'CHEBYSHEV' and 'LEGENDRE'; defaults to
'MONOMIAL';
- 'Delta': difference between M and N in rational approximant;
defaults to 0;
- 'errorEstimatorKind': kind of error estimator; available values
include 'EXACT', 'BASIC', and 'BARE'; defaults to 'EXACT';
- '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.
- '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;
- '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;
- 'Delta': difference between M and N in rational approximant;
- 'errorEstimatorKind': kind of error estimator;
- 'interpRcond': tolerance for interpolation;
- '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 test points.
sampler: Sample point generator.
radialBasis: Radial basis family for interpolant numerator.
radialBasisWeights: Radial basis weights for interpolant numerator.
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.
robustTol: tolerance for robust rational denominator management.
Delta: difference between M and N in rational approximant.
errorEstimatorKind: kind of error estimator.
interpRcond: tolerance for interpolation.
robustTol: tolerance for robust rational denominator management.
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.
"""
_allowedEstimatorKinds = ["EXACT", "BASIC", "BARE"]
def __init__(self, HFEngine:HFEng, mu0 : paramVal = None,
approxParameters : DictAny = {}, homogeneized : bool = False,
verbosity : int = 10, timestamp : bool = True):
self._preInit()
self._addParametersToList(["Delta", "polybasis", "errorEstimatorKind",
"interpRcond", "robustTol"],
[0, "MONOMIAL", "EXACT", -1, 0],
toBeExcluded = ["E"])
super().__init__(HFEngine = HFEngine, mu0 = mu0,
approxParameters = approxParameters,
homogeneized = homogeneized,
verbosity = verbosity, timestamp = timestamp)
vbMng(self, "INIT", "Computing Taylor blocks of system.", 7)
nAs = self.HFEngine.nAs - 1
nbs = max(self.HFEngine.nbs, (nAs + 1) * self.homogeneized)
self.As = [self.HFEngine.A(self.mu0, j + 1) for j in range(nAs)]
self.bs = [self.HFEngine.b(self.mu0, j, self.homogeneized)
for j in range(nbs)]
vbMng(self, "DEL", "Done computing Taylor blocks.", 7)
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("Sample type not recognized.")
self._polybasis = polybasis
except:
RROMPyWarning(("Prescribed polybasis not recognized. Overriding "
"to 'MONOMIAL'."))
self._polybasis = "MONOMIAL"
self._approxParameters["polybasis"] = self.polybasis
@property
def Delta(self):
"""Value of Delta."""
return self._Delta
@Delta.setter
def Delta(self, Delta):
if not np.isclose(Delta, np.floor(Delta)):
raise RROMPyException("Delta must be an integer.")
if Delta < 0:
RROMPyWarning(("Error estimator unreliable for Delta < 0. "
"Overloading of errorEstimator is suggested."))
else:
Deltamin = (max(self.HFEngine.nbs,
self.HFEngine.nAs * self.homogeneized)
- 1 - 1 * (self.HFEngine.nAs > 1))
if Delta < Deltamin:
RROMPyWarning(("Method may be unreliable for selected Delta. "
"Suggested minimal value of Delta: {}.").format(
Deltamin))
self._Delta = Delta
self._approxParameters["Delta"] = self.Delta
@property
def errorEstimatorKind(self):
"""Value of errorEstimatorKind."""
return self._errorEstimatorKind
@errorEstimatorKind.setter
def errorEstimatorKind(self, errorEstimatorKind):
errorEstimatorKind = errorEstimatorKind.upper()
if errorEstimatorKind not in self._allowedEstimatorKinds:
RROMPyWarning(("Error estimator kind not recognized. Overriding "
"to 'EXACT'."))
errorEstimatorKind = "EXACT"
self._errorEstimatorKind = errorEstimatorKind
self._approxParameters["errorEstimatorKind"] = self.errorEstimatorKind
@property
def nTestPoints(self):
"""Value of nTestPoints."""
return self._nTestPoints
@nTestPoints.setter
def nTestPoints(self, nTestPoints):
if nTestPoints <= np.abs(self.Delta):
RROMPyWarning(("nTestPoints must be at least abs(Delta) + 1. "
"Increasing value to abs(Delta) + 1."))
nTestPoints = np.abs(self.Delta) + 1
if not np.isclose(nTestPoints, np.int(nTestPoints)):
raise RROMPyException("nTestPoints must be an integer.")
nTestPoints = np.int(nTestPoints)
if hasattr(self, "_nTestPoints") and self.nTestPoints is not None:
nTestPointsold = self.nTestPoints
else: nTestPointsold = -1
self._nTestPoints = nTestPoints
self._approxParameters["nTestPoints"] = self.nTestPoints
if nTestPointsold != self.nTestPoints:
self.resetSamples()
def _errorSamplingRatio(self, mus:Np1D, muRTest:Np1D,
muRTrain:Np1D) -> Np1D:
"""Scalar ratio in explicit error estimator."""
self.setupApprox()
testTile = np.tile(np.reshape(muRTest, (-1, 1)), [1, len(muRTrain)])
nodalVals = np.prod(testTile - np.reshape(muRTrain, (1, -1)), axis = 1)
denVals = self.trainedModel.getQVal(mus)
return np.abs(nodalVals / denVals)
def _RHSNorms(self, radiusb0:Np2D) -> Np1D:
"""High fidelity system RHS norms."""
self.assembleReducedResidualBlocks(full = False)
# 'ij,jk,ik->k', resbb, radiusb0, radiusb0.conj()
b0resb0 = np.sum(self.trainedModel.data.resbb.dot(radiusb0)
* radiusb0.conj(), axis = 0)
RHSnorms = np.power(np.abs(b0resb0), .5)
return RHSnorms
def _errorEstimatorBare(self) -> Np1D:
"""Bare residual-based error estimator."""
self.setupApprox()
self.assembleReducedResidualGramian(self.trainedModel.data.projMat)
pDom = self.trainedModel.data.P.domCoeff
LL = pDom.conj().dot(self.trainedModel.data.gramian.dot(pDom))
Adiag = self.As[0].diagonal()
LL = ((self.scaleFactor[0] * np.linalg.norm(Adiag)) ** 2. * LL
/ np.size(Adiag))
return np.power(np.abs(LL), .5)
def _errorEstimatorBasic(self, muTest:paramVal, ratioTest:complex) -> Np1D:
"""Basic residual-based error estimator."""
resmu = self.HFEngine.residual(self.trainedModel.getApprox(muTest),
muTest, self.homogeneized,
duality = False)[0]
return np.abs(self.estimatorNormEngine.norm(resmu) / ratioTest)
def _errorEstimatorExact(self, muRTrain:Np1D, vanderBase:Np2D) -> Np1D:
"""Exact residual-based error estimator."""
self.setupApprox()
self.assembleReducedResidualBlocks(full = True)
nAs = self.HFEngine.nAs - 1
nbs = max(self.HFEngine.nbs - 1, nAs * self.homogeneized)
delta = len(self.mus) - self.trainedModel.data.Q.deg[0]
nbsEff = max(0, nbs - delta)
momentQ = np.zeros(nbsEff, dtype = np.complex)
momentQu = np.zeros((len(self.mus), nAs), dtype = np.complex)
radiusbTen = np.zeros((nbsEff, nbsEff, vanderBase.shape[1]),
dtype = np.complex)
radiusATen = np.zeros((nAs, nAs, vanderBase.shape[1]),
dtype = np.complex)
if nbsEff > 0:
momentQ[0] = self.trainedModel.data.Q.domCoeff
radiusbTen[0, :, :] = vanderBase[: nbsEff, :]
momentQu[:, 0] = self.trainedModel.data.P.domCoeff
radiusATen[0, :, :] = vanderBase[: nAs, :]
Qvals = self.trainedModel.getQVal(self.mus)
for k in range(1, max(nAs, nbs * (nbsEff > 0))):
Qvals = Qvals * muRTrain
if k > delta and k < nbs:
momentQ[k - delta] = self._fitinv.dot(Qvals)
radiusbTen[k - delta, k :, :] = (
radiusbTen[0, : delta - k, :])
if k < nAs:
momentQu[:, k] = Qvals * self._fitinv
radiusATen[k, k :, :] = radiusATen[0, : - k, :]
if self.POD and nAs > 1:
momentQu[:, 1 :] = self.samplingEngine.RPOD.dot(
momentQu[:, 1 :])
radiusA = np.tensordot(momentQu, radiusATen, 1)
if nbsEff > 0:
radiusb = np.tensordot(momentQ, radiusbTen, 1)
# 'ij,jk,ik->k', resbb, radiusb, radiusb.conj()
ff = np.sum(self.trainedModel.data.resbb[delta + 1 :, delta + 1 :]\
.dot(radiusb) * radiusb.conj(), axis = 0)
# 'ijk,jkl,il->l', resAb, radiusA, radiusb.conj()
Lf = np.sum(np.tensordot(
self.trainedModel.data.resAb[delta :, :, :], radiusA, 2)
* radiusb.conj(), axis = 0)
else:
ff, Lf = 0., 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.power(np.abs(ff - 2. * np.real(Lf) + LL), .5)
def errorEstimator(self, mus:Np1D) -> Np1D:
"""Standard residual-based error estimator."""
self.setupApprox()
if (np.any(np.isnan(self.trainedModel.data.P.domCoeff))
or np.any(np.isinf(self.trainedModel.data.P.domCoeff))):
err = np.empty(len(mus))
err[:] = np.inf
return err
nAs = self.HFEngine.nAs - 1
nbs = max(self.HFEngine.nbs - 1, nAs * self.homogeneized)
muRTest = self.centerNormalize(mus)(0)
muRTrain = self.centerNormalize(self.mus)(0)
samplingRatio = self._errorSamplingRatio(mus, muRTest, muRTrain)
vanderBase = np.polynomial.polynomial.polyvander(muRTest,
max(nAs, nbs)).T
RHSnorms = self._RHSNorms(vanderBase[: nbs + 1, :])
if self.errorEstimatorKind == "BARE":
jOpt = self._errorEstimatorBare()
elif self.errorEstimatorKind == "BASIC":
idx_muTestSample = np.argmax(samplingRatio)
jOpt = self._errorEstimatorBasic(mus[idx_muTestSample],
samplingRatio[idx_muTestSample])
else: #if self.errorEstimatorKind == "EXACT":
jOpt = self._errorEstimatorExact(muRTrain, vanderBase[: -1, :])
return jOpt * samplingRatio / RHSnorms
def setupApprox(self, plotEst : bool = False):
"""
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.greedy(plotEst)
self._S = len(self.mus)
self._N, self._M, self._E = [self._S - 1] * 3
if self.Delta < 0:
self._M += self.Delta
else:
self._N -= self.Delta
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.scaleFactor,
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)
self.trainedModel.data.mus = copy(self.mus)
if min(self.M, self.N) < 0:
vbMng(self, "MAIN", "Minimal sample size not achieved.", 5)
Q = np.empty(tuple([max(self.N, 0) + 1] * self.npar),
dtype = np.complex)
P = np.empty(tuple([max(self.M, 0) + 1] * self.npar)
+ (len(self.mus),), dtype = np.complex)
Q[:] = np.nan
P[:] = np.nan
self.trainedModel.data.Q = copy(Q)
self.trainedModel.data.P = copy(P)
self.trainedModel.data.approxParameters = copy(
self.approxParameters)
vbMng(self, "DEL", "Aborting computation of approximant.", 5)
return
if self.N > 0:
Q = self._setupDenominator()
else:
Q = np.ones(tuple([1] * self.npar), dtype = np.complex)
self.trainedModel.data.Q = copy(Q)
self.trainedModel.data.P = copy(self._setupNumerator())
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 reduced linear system through projections."""
scaling = self.trainedModel.data.scaleFactor[0]
self.assembleReducedResidualBlocksbb(self.bs, scaling)
if full:
pMat = self.trainedModel.data.projMat
self.assembleReducedResidualBlocksAb(self.As, self.bs[1 :],
pMat, scaling)
self.assembleReducedResidualBlocksAA(self.As, pMat, scaling)
diff --git a/rrompy/reduction_methods/distributed_greedy/rb_distributed_greedy.py b/rrompy/reduction_methods/distributed_greedy/rb_distributed_greedy.py
index 1529dc1..1bb837c 100644
--- a/rrompy/reduction_methods/distributed_greedy/rb_distributed_greedy.py
+++ b/rrompy/reduction_methods/distributed_greedy/rb_distributed_greedy.py
@@ -1,237 +1,238 @@
# 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.
bs: List of numpy vectors representing coefficients of linear system
RHS.
ARBs: List of sparse matrices (in CSC format) representing coefficients
of compressed linear system matrix.
bRBs: List of numpy vectors representing coefficients of compressed
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]
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):
uApproxR = uApproxRs[j]
for i in range(max(nAs, nbs)):
- thetai = (muPL[0] ** self.HFEngine.rescalingExp[0]
- - mu0Eff) ** i
+ thetai = ((muPL[0] ** self.HFEngine.rescalingExp[0]
+ - mu0Eff) / self.scaleFactor[0]) ** 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.computeScaleFactor()
self.greedy(plotEst)
vbMng(self, "INIT", "Computing projection matrix.", 7)
pMat = self.samplingEngine.samples
-
+ vbMng(self, "INIT", "Computing affine blocks of system.", 10)
+ self.As = self.HFEngine.affineLinearSystemA(self.mu0, self.scaleFactor)
+ self.bs = self.HFEngine.affineLinearSystemb(self.mu0, self.scaleFactor,
+ self.homogeneized)
+ vbMng(self, "DEL", "Done computing affine blocks.", 10)
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)
+ pMat, self.scaleFactor,
+ self.HFEngine.rescalingExp)
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 004b3fa..9830125 100644
--- a/rrompy/reduction_methods/pivoting/__init__.py
+++ b/rrompy/reduction_methods/pivoting/__init__.py
@@ -1,35 +1,29 @@
# 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',
- 'RationalInterpolantPivotedPoleMatching',
- 'RBDistributedPivotedPoleMatching'
+ 'RBDistributedPivoted'
]
diff --git a/rrompy/reduction_methods/pivoting/rational_interpolant_pivoted.py b/rrompy/reduction_methods/pivoting/rational_interpolant_pivoted.py
index d12c423..15a07ea 100644
--- a/rrompy/reduction_methods/pivoting/rational_interpolant_pivoted.py
+++ b/rrompy/reduction_methods/pivoting/rational_interpolant_pivoted.py
@@ -1,605 +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.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, 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.scaleFactor,
self.HFEngine.rescalingExp,
- self.directionPivot,
- self.scaleFactor)
+ self.directionPivot)
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/rb_distributed_pivoted.py b/rrompy/reduction_methods/pivoting/rb_distributed_pivoted.py
index f18460b..4fabd41 100644
--- a/rrompy/reduction_methods/pivoting/rb_distributed_pivoted.py
+++ b/rrompy/reduction_methods/pivoting/rb_distributed_pivoted.py
@@ -1,282 +1,284 @@
# 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.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.
bs: List of numpy vectors representing coefficients of linear system
RHS.
ARBs: List of sparse matrices (in CSC format) representing coefficients
of compressed linear system matrix.
bRBs: List of numpy vectors representing coefficients of compressed
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)
+ self.AsList[k] = MEnginek.affineLinearSystemA(mu0Eff,
+ self.scaleFactorPivot)
if not hasbs:
- self.bsList[k] = MEnginek.affineLinearSystemb(mu0Eff,
- self.homogeneized)
+ self.bsList[k] = MEnginek.affineLinearSystemb(mu0Eff,
+ self.scaleFactorPivot,
+ 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)
+ pList, self.scaleFactor,
+ self.HFEngine.rescalingExp,
+ self.directionPivot)
data.musPivot = copy(self.musPivot)
data.musMarginal = copy(self.musMarginal)
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/trained_model/trained_model_data.py b/rrompy/reduction_methods/pole_matching/__init__.py
similarity index 53%
copy from rrompy/reduction_methods/trained_model/trained_model_data.py
copy to rrompy/reduction_methods/pole_matching/__init__.py
index 63620a7..ae1e7db 100644
--- a/rrompy/reduction_methods/trained_model/trained_model_data.py
+++ b/rrompy/reduction_methods/pole_matching/__init__.py
@@ -1,35 +1,29 @@
# 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
-from rrompy.utilities.base.types import Np2D, List, paramVal
-from rrompy.utilities.exception_manager import RROMPyAssert
+from .generic_pole_matching_approximant import GenericPoleMatchingApproximant
+from .rational_interpolant_pole_matching import RationalInterpolantPoleMatching
+from .rb_distributed_pole_matching import RBDistributedPoleMatching
-__all__ = ['TrainedModelData']
+__all__ = [
+ 'GenericPoleMatchingApproximant',
+ 'RationalInterpolantPoleMatching',
+ 'RBDistributedPoleMatching'
+ ]
-class TrainedModelData:
- """ROM approximant evaluation data (must be pickle-able)."""
- def __init__(self, name:str, mu0:paramVal, projMat:Np2D,
- rescalingExp : List[float] = [1.]):
- self.npar = len(rescalingExp)
- RROMPyAssert(mu0.shape[1], self.npar, "Number of parameters")
- self.name = name
- self.mu0 = mu0
- self.projMat = copy(projMat)
- self.rescalingExp = rescalingExp
diff --git a/rrompy/reduction_methods/pivoting/rb_distributed_pivoted_pole_matching.py b/rrompy/reduction_methods/pole_matching/generic_pole_matching_approximant.py
similarity index 65%
copy from rrompy/reduction_methods/pivoting/rb_distributed_pivoted_pole_matching.py
copy to rrompy/reduction_methods/pole_matching/generic_pole_matching_approximant.py
index 19afb42..f831d2a 100644
--- a/rrompy/reduction_methods/pivoting/rb_distributed_pivoted_pole_matching.py
+++ b/rrompy/reduction_methods/pole_matching/generic_pole_matching_approximant.py
@@ -1,198 +1,179 @@
# 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
-from .rb_distributed_pivoted import RBDistributedPivoted
-from rrompy.reduction_methods.trained_model import (TrainedModelPivotedData,
- 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)
+from numpy import inf
+from rrompy.reduction_methods.pivoting.generic_pivoted_approximant import (
+ GenericPivotedApproximant)
+from rrompy.utilities.base.types import ListAny, DictAny, HFEng, paramVal
+from rrompy.utilities.exception_manager import RROMPyException, RROMPyWarning
-__all__ = ['RBDistributedPivotedPoleMatching']
+__all__ = ['GenericPoleMatchingApproximant']
-class RBDistributedPivotedPoleMatching(RBDistributedPivoted):
+class GenericPoleMatchingApproximant(GenericPivotedApproximant):
"""
- ROM pivoted RB approximant computation for parametric problems with pole
- matching.
+ ROM pole matching 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;
- 'matchingWeight': weight for pole matching optimization; defaults
to 1;
+ - 'cutOffTolerance': tolerance for ignoring parasitic poles;
+ defaults to np.inf;
+ - 'cutOffType': rule for tolerance computation for parasitic poles;
+ defaults to 'MAGNITUDE';
- '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.
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;
+ - 'cutOffTolerance': tolerance for ignoring parasitic poles;
+ - 'cutOffType': rule for tolerance computation for parasitic poles;
- '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.
+ cutOffTolerance: Tolerance for ignoring parasitic poles.
+ cutOffType: Rule for tolerance computation for parasitic poles.
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.
- bs: List of numpy vectors representing coefficients of linear system
- RHS.
- ARBs: List of sparse matrices (in CSC format) representing coefficients
- of compressed linear system matrix.
- bRBs: List of numpy vectors representing coefficients of compressed
- 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()
if len(directionPivot) > 1:
raise RROMPyException(("Exactly 1 pivot parameter allowed in pole "
"matching."))
- self._addParametersToList(["matchingWeight"], [1])
+ self._addParametersToList(["matchingWeight", "cutOffTolerance",
+ "cutOffType"], [1, inf, "MAGNITUDE"])
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 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)
- 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)
+ @property
+ def cutOffTolerance(self):
+ """Value of cutOffTolerance."""
+ return self._cutOffTolerance
+ @cutOffTolerance.setter
+ def cutOffTolerance(self, cutOffTolerance):
+ self._cutOffTolerance = cutOffTolerance
+ self._approxParameters["cutOffTolerance"] = self.cutOffTolerance
+
+ @property
+ def cutOffType(self):
+ """Value of cutOffType."""
+ return self._cutOffType
+ @cutOffType.setter
+ def cutOffType(self, cutOffType):
+ try:
+ cutOffType = cutOffType.upper().strip().replace(" ","")
+ if cutOffType not in ["MAGNITUDE", "POTENTIAL"]:
+ raise RROMPyException("Prescribed cutOffType not recognized.")
+ self._cutOffType = cutOffType
+ except:
+ RROMPyWarning(("Prescribed cutOffType not recognized. Overriding "
+ "to 'MAGNITUDE'."))
+ self._cutOffType = "MAGNITUDE"
+ self._approxParameters["cutOffType"] = self.cutOffType
diff --git a/rrompy/reduction_methods/pivoting/rational_interpolant_pivoted_pole_matching.py b/rrompy/reduction_methods/pole_matching/rational_interpolant_pole_matching.py
similarity index 86%
rename from rrompy/reduction_methods/pivoting/rational_interpolant_pivoted_pole_matching.py
rename to rrompy/reduction_methods/pole_matching/rational_interpolant_pole_matching.py
index 5c465f6..bc8eeec 100644
--- a/rrompy/reduction_methods/pivoting/rational_interpolant_pivoted_pole_matching.py
+++ b/rrompy/reduction_methods/pole_matching/rational_interpolant_pole_matching.py
@@ -1,231 +1,222 @@
# 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.reduction_methods.pivoting.rational_interpolant_pivoted import \
+ RationalInterpolantPivoted
+from .generic_pole_matching_approximant import GenericPoleMatchingApproximant
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)
+from rrompy.utilities.exception_manager import RROMPyAssert
-__all__ = ['RationalInterpolantPivotedPoleMatching']
+__all__ = ['RationalInterpolantPoleMatching']
-class RationalInterpolantPivotedPoleMatching(RationalInterpolantPivoted):
+class RationalInterpolantPoleMatching(GenericPoleMatchingApproximant,
+ 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;
+ - 'cutOffTolerance': tolerance for ignoring parasitic poles;
+ defaults to np.inf;
+ - 'cutOffType': rule for tolerance computation for parasitic poles;
+ defaults to 'MAGNITUDE';
- '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;
+ - 'cutOffTolerance': tolerance for ignoring parasitic poles;
+ - 'cutOffType': rule for tolerance computation for parasitic poles;
- '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.
+ cutOffTolerance: Tolerance for ignoring parasitic poles.
+ cutOffType: Rule for tolerance computation for parasitic poles.
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.scaleFactor,
self.HFEngine.rescalingExp,
- self.directionPivot,
- self.scaleFactor)
+ self.directionPivot)
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)
+ if not np.isinf(self.cutOffTolerance):
+ vbMng(self, "INIT", "Recompressing by cut-off.", 10)
+ msg = self.trainedModel.recompressByCutOff([-1., 1.],
+ self.cutOffTolerance,
+ self.cutOffType)
+ vbMng(self, "DEL", "Done recompressing." + msg, 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_pole_matching.py b/rrompy/reduction_methods/pole_matching/rb_distributed_pole_matching.py
similarity index 83%
rename from rrompy/reduction_methods/pivoting/rb_distributed_pivoted_pole_matching.py
rename to rrompy/reduction_methods/pole_matching/rb_distributed_pole_matching.py
index 19afb42..d211dd3 100644
--- a/rrompy/reduction_methods/pivoting/rb_distributed_pivoted_pole_matching.py
+++ b/rrompy/reduction_methods/pole_matching/rb_distributed_pole_matching.py
@@ -1,198 +1,190 @@
# 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 isinf
from copy import deepcopy as copy
-from .rb_distributed_pivoted import RBDistributedPivoted
+from rrompy.reduction_methods.pivoting.rb_distributed_pivoted import \
+ RBDistributedPivoted
+from .generic_pole_matching_approximant import GenericPoleMatchingApproximant
from rrompy.reduction_methods.trained_model import (TrainedModelPivotedData,
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)
+from rrompy.utilities.exception_manager import RROMPyAssert
-__all__ = ['RBDistributedPivotedPoleMatching']
+__all__ = ['RBDistributedPoleMatching']
-class RBDistributedPivotedPoleMatching(RBDistributedPivoted):
+class RBDistributedPoleMatching(GenericPoleMatchingApproximant,
+ RBDistributedPivoted):
"""
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;
+ - 'cutOffTolerance': tolerance for ignoring parasitic poles;
+ defaults to np.inf;
+ - 'cutOffType': rule for tolerance computation for parasitic poles;
+ defaults to 'MAGNITUDE';
- '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.
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;
+ - 'cutOffTolerance': tolerance for ignoring parasitic poles;
+ - 'cutOffType': rule for tolerance computation for parasitic poles;
- '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.
+ cutOffTolerance: Tolerance for ignoring parasitic poles.
+ cutOffType: Rule for tolerance computation for parasitic poles.
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.
bs: List of numpy vectors representing coefficients of linear system
RHS.
ARBs: List of sparse matrices (in CSC format) representing coefficients
of compressed linear system matrix.
bRBs: List of numpy vectors representing coefficients of compressed
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()
- 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 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)
+ pList, self.scaleFactor,
+ self.HFEngine.rescalingExp,
+ self.directionPivot)
data.musPivot = copy(self.musPivot)
data.musMarginal = copy(self.musMarginal)
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)
+ if not isinf(self.cutOffTolerance):
+ vbMng(self, "INIT", "Recompressing by cut-off.", 10)
+ msg = self.trainedModel.recompressByCutOff([- 1., 1.],
+ self.cutOffTolerance,
+ self.cutOffType)
+ vbMng(self, "DEL", "Done recompressing." + msg, 10)
self.trainedModel.data.approxParameters = copy(self.approxParameters)
vbMng(self, "DEL", "Done setting up approximant.", 5)
diff --git a/rrompy/reduction_methods/trained_model/trained_model_data.py b/rrompy/reduction_methods/trained_model/trained_model_data.py
index 63620a7..766f236 100644
--- a/rrompy/reduction_methods/trained_model/trained_model_data.py
+++ b/rrompy/reduction_methods/trained_model/trained_model_data.py
@@ -1,35 +1,37 @@
# 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
from rrompy.utilities.base.types import Np2D, List, paramVal
from rrompy.utilities.exception_manager import RROMPyAssert
__all__ = ['TrainedModelData']
class TrainedModelData:
"""ROM approximant evaluation data (must be pickle-able)."""
def __init__(self, name:str, mu0:paramVal, projMat:Np2D,
+ scaleFactor : List[float] = [1.],
rescalingExp : List[float] = [1.]):
self.npar = len(rescalingExp)
RROMPyAssert(mu0.shape[1], self.npar, "Number of parameters")
self.name = name
self.mu0 = mu0
self.projMat = copy(projMat)
+ self.scaleFactor = scaleFactor
self.rescalingExp = rescalingExp
diff --git a/rrompy/reduction_methods/trained_model/trained_model_pivoted_data.py b/rrompy/reduction_methods/trained_model/trained_model_pivoted_data.py
index d5e6539..d8ba156 100644
--- a/rrompy/reduction_methods/trained_model/trained_model_pivoted_data.py
+++ b/rrompy/reduction_methods/trained_model/trained_model_pivoted_data.py
@@ -1,71 +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 .
#
from .trained_model_data import TrainedModelData
from rrompy.utilities.base.types import Np2D, List, ListAny, paramVal
from rrompy.parameter import checkParameterList
__all__ = ['TrainedModelPivotedData']
class TrainedModelPivotedData(TrainedModelData):
"""ROM approximant evaluation data (must be pickle-able)."""
def __init__(self, name:str, mu0:paramVal, projMat:Np2D,
+ scaleFactor : ListAny = [1.],
rescalingExp : List[float] = [1.],
- directionPivot : ListAny = [0],
- scaleFactor : ListAny = [1.]):
- super().__init__(name, mu0, projMat, rescalingExp)
+ directionPivot : ListAny = [0]):
+ super().__init__(name, mu0, projMat, scaleFactor, rescalingExp)
self.directionPivot = directionPivot
- self.scaleFactor = scaleFactor
@property
def directionMarginal(self):
return tuple([x for x in range(self.npar) \
if x not in self.directionPivot])
@property
def mu0Pivot(self):
return checkParameterList(self.mu0(0, self.directionPivot))[0]
@property
def mu0Marginal(self):
return checkParameterList(self.mu0(0, self.directionMarginal))[0]
@property
def nparPivot(self):
return len(self.directionPivot)
@property
def nparMarginal(self):
return self.npar - self.nparPivot
@property
def rescalingExpPivot(self):
return [self.rescalingExp[x] for x in self.directionPivot]
@property
def rescalingExpMarginal(self):
return [self.rescalingExp[x] for x in self.directionMarginal]
@property
def scaleFactorPivot(self):
return [self.scaleFactor[x] for x in self.directionPivot]
@property
def scaleFactorMarginal(self):
return [self.scaleFactor[x] for x in self.directionMarginal]
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
index 8e71380..fc1720c 100644
--- 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
@@ -1,196 +1,244 @@
# 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.types import (Np1D, Tuple, 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))
+ explored = [0]
+ unexplored = list(range(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)
+ minIex = explored[minIdx // len(unexplored)]
+ minIunex = unexplored[minIdx % len(unexplored)]
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]
+ explored += [minIunex]
+ unexplored.remove(minIunex)
+# 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 recompressByCutOff(self, murange : Tuple[float, float] = [- 1., 1.],
+ tol : float = np.inf, rtype : str = "MAGNITUDE"):
+ if np.isinf(tol): return " No poles erased."
+ N = len(self.data.HIs[0].poles)
+ mu0 = np.mean(murange)
+ musig = murange[0] - mu0
+ if np.isclose(musig, 0.):
+ radius = lambda x: np.abs(x - mu0)
+ else:
+ if rtype == "MAGNITUDE":
+ murdir = (murange[0] - mu0) / np.abs(musig)
+ def radius(x):
+ scalprod = (x - mu0) * murdir.conj() / np.abs(musig)
+ rescalepar = np.abs(np.real(scalprod))
+ rescaleort = np.abs(np.imag(scalprod))
+ return ((rescalepar - 1.) ** 2. * (rescalepar > 1.)
+ + rescaleort ** 2.) ** .5
+ else:#if rtype == "POTENTIAL":
+ def radius(x):
+ rescale = (x - mu0) / musig
+ return np.max(np.abs(rescale * np.array([-1., 1.])
+ + (rescale ** 2. - 1) ** .5)) - 1.
+ keepPole, removePole = [], []
+ for j in range(N):
+ for hi in self.data.HIs:
+ if radius(hi.poles[j]) <= tol:
+ keepPole += [j]
+ break
+ if len(keepPole) == 0 or keepPole[-1] != j: removePole += [j]
+ if len(keepPole) == N: return " No poles erased."
+ keepCoeff = keepPole + [N]
+ keepCoeff = keepCoeff + list(range(N + 1,len(self.data.HIs[0].coeffs)))
+ for hi in self.data.HIs:
+ polyCorrection = np.zeros_like(hi.coeffs[0, :])
+ for j in removePole:
+ polyCorrection += hi.coeffs[j, :] / (mu0 - hi.poles[j])
+ if len(hi.coeffs) == N:
+ hi.coeffs = np.vstack((hi.coeffs, polyCorrection))
+ else:
+ hi.coeffs[N, :] += polyCorrection
+ hi.poles = hi.poles[keepPole]
+ hi.coeffs = hi.coeffs[keepCoeff, :]
+ return (" Erased {} pole".format(len(removePole))
+ + "s" * (len(removePole) > 1) + ".")
+
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])
+ if isinstance(cs, (list, tuple,)): cs = np.array(cs)
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)
+ "Assembling reduced model for mu = {}.".format(muPL), 87)
hsL = self.interpolateMarginalInterpolator(muM)
- vbMng(self, "DEL", "Done assembling reduced model.", 17)
+ vbMng(self, "DEL", "Done assembling reduced model.", 87)
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 36bdbd4..2417a25 100644
--- a/rrompy/reduction_methods/trained_model/trained_model_pivoted_rational.py
+++ b/rrompy/reduction_methods/trained_model/trained_model_pivoted_rational.py
@@ -1,224 +1,225 @@
# 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)
+ 95)
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.reshape(len(sEff[0]), - 1) * mIj(muC, der, scl)
- vbMng(self, "DEL", "Done interpolating marginal.", 35)
+ if len(sEff[0]) > 0:
+ for mIj, spj in zip(self.data.marginalInterp, sEff):
+ p = p + spj.reshape(len(sEff[0]), - 1) * mIj(muC, der, scl)
+ vbMng(self, "DEL", "Done interpolating marginal.", 95)
if not sList:
p = p.data.flatten()
return p
def getPValsPivot(self, mu : paramList = []) -> sampList:
"""
Evaluate rational numerators at arbitrary pivot parameter.
Args:
mu: Target parameter.
"""
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)) for P in self.data.Ps]
vbMng(self, "DEL", "Done evaluating numerator.", 17)
return pP
def getPVal(self, mu : paramList = []) -> sampList:
"""
Evaluate rational numerator at arbitrary parameter.
Args:
mu: Target parameter.
"""
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]
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_rb.py b/rrompy/reduction_methods/trained_model/trained_model_pivoted_rb.py
index 7a5b260..33c669d 100644
--- a/rrompy/reduction_methods/trained_model/trained_model_pivoted_rb.py
+++ b/rrompy/reduction_methods/trained_model/trained_model_pivoted_rb.py
@@ -1,205 +1,202 @@
# 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_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.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)
+ if mu0 is None: mu0 = self.data.mu0Pivot
rad = ((mu ** self.data.rescalingExpPivot
- - mu0 ** 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)
+ if mu0 is None: mu0 = self.data.mu0Marginal
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]
+ sEff[j] = samples[j].data
else:
sEff[j] = samples[j]
vbMng(self, "INIT", "Interpolating marginal at mu = {}.".format(mu),
- 35)
+ 95)
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)
+ p = np.zeros_like(sEff[0], dtype = sEff[0].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)
+ p += mIj(muC) * spj
+ vbMng(self, "DEL", "Done interpolating marginal.", 95)
+# if sList: p = sampleList(p)
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]
if (not hasattr(self, "lastSolvedApproxReduced")
or self.lastSolvedApproxReduced != mu):
vbMng(self, "INIT",
"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)
- 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)
+ "Assembling reduced model for mu = {}.".format(muPL), 87)
ARBs = self.interpolateMarginal(muM, self.data.ARBsList)
bRBs = self.interpolateMarginal(muM, self.data.bRBsList)
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))
+ theta = np.prod(self.centerNormalizePivot(muP)
+ ** 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, "DEL", "Done assembling reduced model.", 87)
vbMng(self, "INIT",
- "Solving reduced model for mu = {}.".format(muPL), 15)
+ "Solving reduced model for mu = {}.".format(muPL), 75)
self.uApproxReduced[i] = np.linalg.solve(ARBmu, bRBmu)
- vbMng(self, "DEL", "Done solving reduced model.", 15)
+ vbMng(self, "DEL", "Done solving reduced model.", 75)
vbMng(self, "DEL", "Done computing RB solution.", 12)
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], jSupp : int = 1,
**kwargs) -> Np1D:
"""
Obtain approximant poles.
Returns:
Numpy complex vector of poles.
"""
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."))
if rDim in self.data.directionMarginal:
raise RROMPyException(("'freepar' entry in marginalVals must be "
"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)
if self.data.nparPivot > 1:
mValsP[rDimP] = self.data.mu0Pivot(rDimP)
- mValsP = (checkParameter(mValsP) ** self.data.rescalingExpPivot
- - self.data.mu0Pivot ** self.data.rescalingExpPivot)
+ mValsP = self.centerNormalizePivot(mValsP)
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])
+ + ev * self.scaleFactor[rDim],
+ 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
index 608bdf4..f2cf22d 100644
--- 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
@@ -1,47 +1,55 @@
# 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_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):
+ def initializeFromAffine(self, HFEngine:HFEng, matchingWeight : float = 1.,
+ POD : bool = True, jSupp : int = 1):
"""Initialize Heaviside representation."""
RROMPyAssert(self.data.nparPivot, 1, "Number of pivot parameters")
poles, coeffs = [], []
+ nbadpls = 0
for As, bs in zip(self.data.ARBsList, self.data.bRBsList):
- cfs, pls, basis = affine2heaviside(As, bs, self.data.mu0Pivot)
+ cfs, pls, basis = affine2heaviside(As, bs, jSupp)
poles += [pls]
coeffs += [cfs]
+ nbadpls = max(nbadpls, np.sum(np.isinf(pls)))
+ if nbadpls > 0:
+ for j in range(len(poles)):
+ plsgood = np.argsort(np.abs(pls))[: - nbadpls]
+ poles[j] = poles[j][plsgood]
+ coeffs[j] = coeffs[j][plsgood, :]
self.initializeFromLists(poles, coeffs, basis, HFEngine,
matchingWeight, POD)
diff --git a/rrompy/reduction_methods/trained_model/trained_model_rb.py b/rrompy/reduction_methods/trained_model/trained_model_rb.py
index e656919..bd5ad7f 100644
--- a/rrompy/reduction_methods/trained_model/trained_model_rb.py
+++ b/rrompy/reduction_methods/trained_model/trained_model_rb.py
@@ -1,106 +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
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.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]
if (not hasattr(self, "lastSolvedApproxReduced")
or self.lastSolvedApproxReduced != mu):
vbMng(self, "INIT",
"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)
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))
+ derIdx = hashI(j, self.data.npar)
+ theta = np.prod(((mu[i] ** self.data.rescalingExp - mu0Eff)
+ / self.data.scaleFactor) ** derIdx)
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[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
return self.uApproxReduced
def getPoles(self, marginalVals : ListAny = [fp], jSupp : int = 1,
**kwargs) -> Np1D:
"""
Obtain approximant poles.
Returns:
Numpy complex vector of poles.
"""
if not hasattr(marginalVals, "__len__"): marginalVals = [marginalVals]
mVals = list(marginalVals)
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."))
ARBs = self.data.ARBs
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]
- + ev, 1. / self.data.rescalingExp[rDim])
+ + ev * self.data.scaleFactor[rDim],
+ 1. / self.data.rescalingExp[rDim])
diff --git a/rrompy/utilities/numerical/point_matching.py b/rrompy/utilities/numerical/point_matching.py
index b3c9cfa..9029e02 100644
--- a/rrompy/utilities/numerical/point_matching.py
+++ b/rrompy/utilities/numerical/point_matching.py
@@ -1,156 +1,165 @@
# 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:
+ max_iter : int = 10, verbObj = None,
+ expl_threshold : int = 8) -> 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)
+ if len(c1) > expl_threshold:
+ optP = np.random.permutation(c1)[: len(c2)]
+ max_iter += int(np.ceil(np.log2(len(c1))))
+ else:
+ optP = matchExplicit(distanceMatrixEff.T, c2, c1)
for i, k in enumerate(optP):
matches[k] = c2[i]
else:
- optP = matchExplicit(distanceMatrixEff, c1, c2)
+ if len(c2) > expl_threshold:
+ optP = np.random.permutation(c2)[: len(c1)]
+ max_iter += int(np.ceil(np.log2(len(c2))))
+ 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)
+ verbObj, expl_threshold)
return matchExplicit(distanceMatrix, idx1, idx2, verbObj)
diff --git a/rrompy/utilities/poly_fitting/heaviside/heaviside_to_from_affine.py b/rrompy/utilities/poly_fitting/heaviside/heaviside_to_from_affine.py
index 13527a1..11cbc3b 100644
--- a/rrompy/utilities/poly_fitting/heaviside/heaviside_to_from_affine.py
+++ b/rrompy/utilities/poly_fitting/heaviside/heaviside_to_from_affine.py
@@ -1,94 +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, List, ListAny, Tuple, paramVal
+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', 'affine2heaviside']
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)))],
[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)
return c.reshape(c.shape[0], -1).T, As, bs
-def affine2heaviside(As:ListAny, bs:ListAny, mu : paramVal = [],
+def affine2heaviside(As:ListAny, bs:ListAny,
jSupp : int = 1) -> Tuple[Np2D, Np1D, str]:
- mu = checkParameter(mu, 1)(0, 0)
+ if jSupp != 1 and not (isinstance(jSupp, (int,))
+ and jSupp.upper() == "COMPANION"):
+ raise RROMPyException(("Affine to heaviside conversion does not allow "
+ "nonlinear eigenproblem support outside first "
+ "block row."))
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]
+ Q = Q[:, : 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