Page MenuHomec4science
Files F61472253

approximant_taylor_pade.py
No OneTemporary
Actions

File Metadata

Created
Mon, May 6, 21:54

approximant_taylor_pade.py
View Options

# Copyright (C) 2018 by the RROMPy authors
#
# This file is part of RROMPy.
#
# RROMPy is free software: you can redistribute it and/or modify
# it under the terms of the GNU Lesser General Public License as published by
# the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
#
# RROMPy is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU Lesser General Public License for more details.
#
# You should have received a copy of the GNU Lesser General Public License
# along with RROMPy. If not, see <http://www.gnu.org/licenses/>.
#
from copy import copy
import numpy as np
from rrompy.reduction_methods.base import checkRobustTolerance
from .generic_approximant_taylor import GenericApproximantTaylor
from rrompy.sampling.base.pod_engine import PODEngine
from rrompy.utilities.base.types import Np1D, Np2D, List, Tuple, DictAny
from rrompy.utilities.base.types import HFEng
from rrompy.utilities.base import purgeDict, verbosityDepth
from rrompy.utilities.warning_manager import warn
__all__ = ['ApproximantTaylorPade']
class ApproximantTaylorPade(GenericApproximantTaylor):
"""
ROM single-point fast Pade' 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;
- 'rho': weight for computation of original Pade' approximant;
defaults to np.inf, i.e. fast approximant;
- 'M': degree of Pade' approximant numerator; defaults to 0;
- 'N': degree of Pade' approximant denominator; defaults to 0;
- 'E': total number of derivatives current approximant relies upon;
defaults to Emax;
- 'Emax': total number of derivatives of solution map to be
computed; defaults to E;
- 'robustTol': tolerance for robust Pade' denominator management;
defaults to 0;
- 'sampleType': label of sampling type; available values are:
- 'ARNOLDI': orthogonalization of solution derivatives through
Arnoldi algorithm;
- 'KRYLOV': standard computation of solution derivatives.
Defaults to 'KRYLOV'.
Defaults to empty dict.
homogeneized: Whether to homogeneize Dirichlet BCs. Defaults to False.
verbosity(optional): Verbosity level. Defaults to 10.
Attributes:
HFEngine: HF problem solver.
mu0: Default parameter.
homogeneized: Whether to homogeneize Dirichlet BCs.
approxParameters: Dictionary containing values for main parameters of
approximant. Recognized keys are in parameterList.
parameterList: Recognized keys of approximant parameters:
- 'POD': whether to compute POD of snapshots;
- 'rho': weight for computation of original Pade' approximant;
- 'M': degree of Pade' approximant numerator;
- 'N': degree of Pade' approximant denominator;
- 'E': total number of derivatives current approximant relies upon;
- 'Emax': total number of derivatives of solution map to be
computed;
- 'robustTol': tolerance for robust Pade' denominator management;
- 'sampleType': label of sampling type.
POD: Whether to compute QR factorization of derivatives.
rho: Weight of approximant.
M: Numerator degree of approximant.
N: Denominator degree of approximant.
E: Number of solution derivatives over which current approximant is
based upon.
Emax: Total number of solution derivatives to be computed.
robustTol: Tolerance for robust Pade' denominator management.
sampleType: Label of sampling type.
initialHFData: HF problem initial data.
uHF: High fidelity solution with wavenumber lastSolvedHF as numpy
complex vector.
lastSolvedHF: Wavenumber corresponding to last computed high fidelity
solution.
G: Square Numpy 2D vector of size (N+1) corresponding to Pade'
denominator matrix (see paper).
Q: Numpy 1D vector containing complex coefficients of approximant
denominator.
P: Numpy 2D vector whose columns are FE dofs of coefficients of
approximant numerator.
uApp: Last evaluated approximant as numpy complex vector.
lastApproxParameters: List of parameters corresponding to last
computed approximant.
"""
def __init__(self, HFEngine:HFEng, mu0 : complex = 0,
approxParameters : DictAny = {}, homogeneized : bool = False,
verbosity : int = 10):
self._preInit()
self._addParametersToList(["M", "N", "robustTol", "rho"])
super().__init__(HFEngine = HFEngine, mu0 = mu0,
approxParameters = approxParameters,
homogeneized = homogeneized,
verbosity = verbosity)
self._postInit()
@property
def approxParameters(self):
"""Value of approximant parameters."""
return self._approxParameters
@approxParameters.setter
def approxParameters(self, approxParams):
approxParameters = purgeDict(approxParams, self.parameterList,
dictname = self.name() + ".approxParameters",
baselevel = 1)
approxParametersCopy = purgeDict(approxParameters,
["M", "N", "robustTol", "rho"],
True, True, baselevel = 1)
keyList = list(approxParameters.keys())
if "rho" in keyList:
self._rho = approxParameters["rho"]
elif hasattr(self, "rho"):
self._rho = self.rho
else:
self._rho = np.inf
GenericApproximantTaylor.approxParameters.fset(self,
approxParametersCopy)
self.rho = self._rho
if "robustTol" in keyList:
self.robustTol = approxParameters["robustTol"]
elif hasattr(self, "robustTol"):
self.robustTol = self.robustTol
else:
self.robustTol = 0
self._ignoreParWarnings = True
if "M" in keyList:
self.M = approxParameters["M"]
elif hasattr(self, "M"):
self.M = self.M
else:
self.M = 0
del self._ignoreParWarnings
if "N" in keyList:
self.N = approxParameters["N"]
elif hasattr(self, "N"):
self.N = self.N
else:
self.N = 0
@property
def rho(self):
"""Value of rho."""
return self._rho
@rho.setter
def rho(self, rho):
self._rho = np.abs(rho)
self._approxParameters["rho"] = self.rho
@property
def M(self):
"""Value of M. Its assignment may change Emax and E."""
return self._M
@M.setter
def M(self, M):
if M < 0: raise ArithmeticError("M must be non-negative.")
self._M = M
self._approxParameters["M"] = self.M
if not hasattr(self, "_ignoreParWarnings"):
self.checkMNEEmax()
@property
def N(self):
"""Value of N. Its assignment may change Emax."""
return self._N
@N.setter
def N(self, N):
if N < 0: raise ArithmeticError("N must be non-negative.")
self._N = N
self._approxParameters["N"] = self.N
self.checkMNEEmax()
def checkMNEEmax(self):
"""Check consistency of M, N, E, and Emax."""
M = self.M if hasattr(self, "_M") else 0
N = self.N if hasattr(self, "_N") else 0
E = self.E if hasattr(self, "_E") else M + N
Emax = self.Emax if hasattr(self, "_Emax") else M + N
if self.rho == np.inf:
if Emax < max(N, M):
warn(("Prescribed Emax is too small. Updating Emax to "
"max(M, N)."))
self.Emax = max(N, M)
if E < max(N, M):
warn("Prescribed E is too small. Updating E to max(M, N).")
self.E = max(N, M)
else:
if Emax < N + M:
warn("Prescribed Emax is too small. Updating Emax to M + N.")
self.Emax = self.N + M
if E < N + M:
warn("Prescribed E is too small. Updating E to M + N.")
self.E = self.N + M
@property
def robustTol(self):
"""Value of tolerance for robust Pade' denominator management."""
return self._robustTol
@robustTol.setter
def robustTol(self, robustTol):
if robustTol < 0.:
warn("Overriding prescribed negative robustness tolerance to 0.")
robustTol = 0.
self._robustTol = robustTol
self._approxParameters["robustTol"] = self.robustTol
@property
def E(self):
"""Value of E. Its assignment may change Emax."""
return self._E
@E.setter
def E(self, E):
if E < 0: raise ArithmeticError("E must be non-negative.")
self._E = E
self.checkMNEEmax()
self._approxParameters["E"] = self.E
if hasattr(self, "Emax") and self.Emax < self.E:
warn("Prescribed Emax is too small. Updating Emax to E.")
self.Emax = self.E
@property
def Emax(self):
"""Value of Emax. Its assignment may reset computed derivatives."""
return self._Emax
@Emax.setter
def Emax(self, Emax):
if Emax < 0: raise ArithmeticError("Emax must be non-negative.")
if hasattr(self, "Emax"): EmaxOld = self.Emax
else: EmaxOld = -1
if hasattr(self, "_N"): N = self.N
else: N = 0
if hasattr(self, "_M"): M = self.M
else: M = 0
if hasattr(self, "_E"): E = self.E
else: E = 0
if self.rho == np.inf:
if max(N, M, E) > Emax:
warn(("Prescribed Emax is too small. Updating Emax to "
"max(N, M, E)."))
Emax = max(N, M, E)
else:
if max(N + M, E) > Emax:
warn(("Prescribed Emax is too small. Updating Emax to "
"max(N + M, E)."))
Emax = max(N + M, E)
self._Emax = Emax
self._approxParameters["Emax"] = Emax
if EmaxOld >= self.Emax and self.samplingEngine.samples is not None:
self.samplingEngine.samples = self.samplingEngine.samples[:,
: self.Emax + 1]
if (self.sampleType == "ARNOLDI"
and self.samplingEngine.HArnoldi is not None):
self.samplingEngine.HArnoldi = self.samplingEngine.HArnoldi[
: self.Emax + 1,
: self.Emax + 1]
self.samplingEngine.RArnoldi = self.samplingEngine.RArnoldi[
: self.Emax + 1,
: self.Emax + 1]
def setupApprox(self):
"""
Compute Pade' approximant. SVD-based robust eigenvalue management.
"""
if not self.checkComputedApprox():
if self.verbosity >= 5:
verbosityDepth("INIT", "Setting up {}.". format(self.name()))
self.computeDerivatives()
if self.N > 0:
if self.verbosity >= 7:
verbosityDepth("INIT", ("Starting computation of "
"denominator."))
while self.N > 0:
if self.POD:
ev, eV = self.findeveVGQR()
else:
ev, eV = self.findeveVGExplicit()
newParameters = checkRobustTolerance(ev, self.E,
self.robustTol)
if not newParameters:
break
self.approxParameters = newParameters
if self.N <= 0:
eV = np.ones((1, 1))
self.Q = np.poly1d(eV[:, 0])
if self.verbosity >= 7:
verbosityDepth("DEL", "Done computing denominator.")
else:
self.Q = np.poly1d([1])
if self.verbosity >= 7:
verbosityDepth("INIT", "Starting computation of numerator.")
self.P = np.zeros((self.Emax + 1, self.M + 1), dtype = np.complex)
for i in range(self.Q.order):
rng = np.arange(self.M + 1 - i)
self.P[rng, - 1 - rng - i] = self.Q[i]
if self.sampleType == "ARNOLDI":
self.P = self.samplingEngine.RArnoldi.dot(self.P)
if self.verbosity >= 7:
verbosityDepth("DEL", "Done computing numerator.")
self.lastApproxParameters = copy(self.approxParameters)
if hasattr(self, "lastSolvedApp"): del self.lastSolvedApp
if self.verbosity >= 5:
verbosityDepth("DEL", "Done setting up approximant.\n")
def rescaleParameter(self, R:Np2D, A : Np2D = None,
exponent : float = 1.) -> Np2D:
"""
Prepare parameter rescaling.
Args:
R: Matrix whose columns need rescaling.
A(optional): Matrix whose diagonal defines scaling factor. If None,
previous value of scaleFactor is used. Defaults to None.
exponent(optional): Exponent of scaling factor in matrix diagonal.
Defaults to 1.
Returns:
Rescaled matrix.
"""
if A is not None:
aDiag = np.diag(A)
scaleCoeffs = np.polyfit(np.arange(A.shape[1]),
np.log(aDiag), 1)
self.scaleFactor = np.exp(- scaleCoeffs[0] / exponent)
return np.multiply(R, np.power(self.scaleFactor,np.arange(R.shape[1])))
def buildG(self):
"""Assemble Pade' denominator matrix."""
self.computeDerivatives()
if self.rho == np.inf:
Nmin = self.E - self.N
else:
Nmin = self.M - self.N + 1
if self.sampleType == "KRYLOV":
DerE = self.samplingEngine.samples[:, Nmin : self.E + 1]
G = self.HFEngine.innerProduct(DerE, DerE)
DerE = self.rescaleParameter(DerE, G, 2.)
G = self.HFEngine.innerProduct(DerE, DerE)
else:
RArnE = self.samplingEngine.RArnoldi[: self.E + 1,
Nmin : self.E + 1]
RArnE = self.rescaleParameter(RArnE, RArnE[Nmin :, :])
G = RArnE.conj().T.dot(RArnE)
if self.rho == np.inf:
self.G = G
else:
Gbig = G
self.G = np.zeros((self.N + 1, self.N + 1), dtype = np.complex)
for k in range(self.E - self.M):
self.G += self.rho ** (2 * k) * Gbig[k : k + self.N + 1,
k : k + self.N + 1]
def findeveVGExplicit(self) -> Tuple[Np1D, Np2D]:
"""
Compute explicitly eigenvalues and eigenvectors of Pade' denominator
matrix.
"""
if self.verbosity >= 10:
verbosityDepth("INIT", "Building gramian matrix.")
self.buildG()
if self.verbosity >= 10:
verbosityDepth("DEL", "Done building gramian.")
if self.verbosity >= 7:
verbosityDepth("INIT",
"Solving eigenvalue problem for gramian matrix.")
ev, eV = np.linalg.eigh(self.G)
eV = self.rescaleParameter(eV.T).T
if self.verbosity >= 5:
try: condev = ev[-1] / ev[0]
except: condev = np.inf
verbosityDepth("MAIN", ("Solved eigenvalue problem of size {} "
"with condition number {:.4e}.").format(
self.N + 1, condev))
if self.verbosity >= 7:
verbosityDepth("DEL", "Done solving eigenvalue problem.")
return ev, eV
def findeveVGQR(self) -> Tuple[Np1D, Np2D]:
"""
Compute eigenvalues and eigenvectors of Pade' denominator matrix
through SVD of R factor. See ``Householder triangularization of a
quasimatrix'', L.Trefethen, 2008 for QR algorithm.
Returns:
Eigenvalues in ascending order and corresponding eigenvector
matrix.
"""
self.computeDerivatives()
if self.rho == np.inf:
Nmin = self.E - self.N
else:
Nmin = self.M - self.N + 1
if self.sampleType == "KRYLOV":
A = copy(self.samplingEngine.samples[:, Nmin : self.E + 1])
self.PODEngine = PODEngine(self.HFEngine)
if self.verbosity >= 10:
verbosityDepth("INIT", "Orthogonalizing samples.", end = "")
R = self.PODEngine.QRHouseholder(A, only_R = True)
if self.verbosity >= 10:
verbosityDepth("DEL", " Done.", inline = True)
else:
R = self.samplingEngine.RArnoldi[: self.E + 1, Nmin : self.E + 1]
R = self.rescaleParameter(R, R[Nmin :, :])
if self.rho == np.inf:
if self.verbosity >= 10:
verbosityDepth("INIT", ("Solving svd for square root of "
"gramian matrix."), end = "")
sizeI = R.shape[0]
_, s, V = np.linalg.svd(R, full_matrices = False)
else:
if self.verbosity >= 10:
verbosityDepth("INIT", ("Building matrix stack for square "
"root of gramian."), end = "")
Rtower = np.zeros((R.shape[0] * (self.E - self.M), self.N + 1),
dtype = np.complex)
for k in range(self.E - self.M):
Rtower[k * R.shape[0] : (k + 1) * R.shape[0], :] = (
self.rho ** k
* R[:, self.M - self.N + 1 + k : self.M + 1 + k])
if self.verbosity >= 10:
verbosityDepth("DEL", " Done.", inline = True)
if self.verbosity >= 7:
verbosityDepth("INIT", ("Solving svd for square root of "
"gramian matrix."), end = "")
sizeI = Rtower.shape[0]
_, s, V = np.linalg.svd(Rtower, full_matrices = False)
eV = V.conj().T[:, ::-1]
eV = self.rescaleParameter(eV.T).T
if self.verbosity >= 5:
try: condev = s[0] / s[-1]
except: condev = np.inf
verbosityDepth("MAIN", ("Solved svd problem of size {} x {} with "
"condition number {:.4e}.").format(
sizeI, self.N + 1, condev))
if self.verbosity >= 7:
verbosityDepth("DEL", " Done.", inline = True)
return s[::-1], eV
def radiusPade(self, mu:Np1D, mu0 : float = None) -> float:
"""
Compute translated radius to be plugged into Pade' approximant.
Args:
mu: Parameter(s) 1.
mu0: Parameter(s) 2. If None, set to self.mu0.
Returns:
Translated radius to be plugged into Pade' approximant.
"""
if mu0 is None: mu0 = self.mu0
return self.HFEngine.rescaling(mu) - self.HFEngine.rescaling(mu0)
def getPVal(self, mu:List[complex]):
"""
Evaluate Pade' numerator at arbitrary parameter.
Args:
mu: Target parameter.
"""
self.setupApprox()
if self.verbosity >= 10:
verbosityDepth("INIT",
"Evaluating numerator at mu = {}.".format(mu),
end = "")
try:
len(mu)
except:
mu = [mu]
powerlist = np.vander(self.radiusPade(mu), self.M + 1).T
p = self.P.dot(powerlist)
if len(mu) == 1:
p = p.flatten()
if self.verbosity >= 10:
verbosityDepth("DEL", " Done.", inline = True)
return p
def getQVal(self, mu:List[complex]):
"""
Evaluate Pade' denominator at arbitrary parameter.
Args:
mu: Target parameter.
"""
self.setupApprox()
if self.verbosity >= 10:
verbosityDepth("INIT",
"Evaluating denominator at mu = {}.".format(mu),
end = "")
q = self.Q(self.radiusPade(mu))
if self.verbosity >= 10:
verbosityDepth("DEL", " Done.", inline = True)
return q
def evalApproxReduced(self, mu:complex):
"""
Evaluate Pade' approximant at arbitrary parameter.
Args:
mu: Target parameter.
"""
self.setupApprox()
if (not hasattr(self, "lastSolvedApp")
or not np.isclose(self.lastSolvedApp, mu)):
if self.verbosity >= 5:
verbosityDepth("INIT",
"Evaluating approximant at mu = {}.".format(mu))
self.uAppReduced = self.getPVal(mu) / self.getQVal(mu)
self.lastSolvedApp = mu
if self.verbosity >= 5:
verbosityDepth("DEL", "Done evaluating approximant.")
def evalApprox(self, mu:complex):
"""
Evaluate approximant at arbitrary parameter.
Args:
mu: Target parameter.
"""
self.evalApproxReduced(mu)
self.uApp = self.samplingEngine.samples.dot(self.uAppReduced)
def getPoles(self) -> Np1D:
"""
Obtain approximant poles.
Returns:
Numpy complex vector of poles.
"""
self.setupApprox()
return self.HFEngine.rescalingInv(self.Q.r
+ self.HFEngine.rescaling(self.mu0))

Event Timeline

c4science · Help