Page MenuHomec4science
Files F60932487

rational_interpolant.py
No OneTemporary
Actions

File Metadata

Created
Fri, May 3, 11:37

rational_interpolant.py
View Options

# Copyright (C) 2018-2020 by the RROMPy authors
#
# This file is part of RROMPy.
#
# RROMPy is free software: you can redistribute it and/or modify
# it under the terms of the GNU Lesser General Public License as published by
# the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
#
# RROMPy is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU Lesser General Public License for more details.
#
# You should have received a copy of the GNU Lesser General Public License
# along with RROMPy. If not, see <http://www.gnu.org/licenses/>.
#
from copy import deepcopy as copy
import numpy as np
from scipy.sparse import coo_matrix
from scipy.linalg import eig, eigvals, khatri_rao
from .generic_standard_approximant import GenericStandardApproximant
from rrompy.utilities.poly_fitting.polynomial import (
polybases, polyfitname,
polyvander as pvP, polyTimes,
PolynomialInterpolator as PI,
PolynomialInterpolatorNodal as PIN)
from rrompy.utilities.poly_fitting.heaviside import rational2heaviside
from rrompy.utilities.base.types import (Np1D, Np2D, Tuple, List, paramList,
interpEng)
from rrompy.utilities.base import verbosityManager as vbMng
from rrompy.utilities.numerical import pseudoInverse, dot, baseDistanceMatrix
from rrompy.utilities.numerical.factorials import multifactorial
from rrompy.utilities.numerical.degree import (mapTotalToTensor, degreeToDOFs,
reduceDegreeMN, reduceDegreeN)
from rrompy.utilities.numerical.hash_derivative import (nextDerivativeIndices,
hashDerivativeToIdx as hashD,
hashIdxToDerivative as hashI)
from rrompy.utilities.exception_manager import (RROMPyException, RROMPyAssert,
RROMPyWarning)
__all__ = ['RationalInterpolant']
def polyTimesTable(P:interpEng, mus:Np1D, reorder:List[int],
derIdxs:List[List[List[int]]], scl : Np1D = None) -> Np2D:
"""Table of polynomial products."""
if not isinstance(P, PI):
raise RROMPyException(("Polynomial to evaluate must be a polynomial "
"interpolator."))
Pvals = [[0.] * len(derIdx) for derIdx in derIdxs]
for j, derIdx in enumerate(derIdxs):
nder = len(derIdx)
for der in range(nder):
derI = hashI(der, P.npar)
Pvals[j][der] = P([mus[j]], derI, scl) / multifactorial(derI)
return blockDiagDer(Pvals, reorder, derIdxs, format = "dense")
def vanderInvTable(vanInv:Np2D, reorder:List[int],
derIdxs:List[List[List[int]]]) -> Np2D:
"""Table of Vandermonde pseudo-inverse."""
S = len(reorder)
Ts = [None] * len(vanInv)
for k in range(len(vanInv)):
invLocs = [None] * len(derIdxs)
idxGlob = 0
for j, derIdx in enumerate(derIdxs):
nder = len(derIdx)
idxGlob += nder
idxLoc = np.arange(S)[(reorder >= idxGlob - nder)
* (reorder < idxGlob)]
invLocs[j] = vanInv[k, idxLoc]
Ts[k] = blockDiagDer(invLocs, reorder, derIdxs, [2, 1, 0])
return Ts
def blockDiagDer(vals:List[Np1D], reorder:List[int],
derIdxs:List[List[List[int]]],
permute : List[int] = None, format : str = "sparse") -> Np2D:
"""Table of derivative values for point confluence."""
S = len(reorder)
if format == "sparse":
data, idxI, idxJ = [], [], []
else: #format == "dense":
T = np.zeros((S, S), dtype = np.complex)
if permute is None: permute = [0, 1, 2]
idxGlob = 0
for j, derIdx in enumerate(derIdxs):
nder = len(derIdx)
idxGlob += nder
idxLoc = np.arange(S)[(reorder >= idxGlob - nder)
* (reorder < idxGlob)]
val = vals[j]
for derI, derIdxI in enumerate(derIdx):
for derJ, derIdxJ in enumerate(derIdx):
diffIdx = [x - y for (x, y) in zip(derIdxI, derIdxJ)]
if all([x >= 0 for x in diffIdx]):
diffj = hashD(diffIdx)
i1, i2, i3 = np.array([derI, derJ, diffj])[permute]
if format == "sparse":
data += [val[i3]]
idxI += [idxLoc[i1]]
idxJ += [idxLoc[i2]]
else:
T[idxLoc[i1], idxLoc[i2]] = val[i3]
if format == "sparse":
T = coo_matrix((data, (idxI, idxJ)), shape = (S, S)).tocsr()
return T
def barycentricRoots(coeffs:Np1D, supp:Np1D,
return_vectors : bool = False) -> Np1D:
m = len(supp)
if len(coeffs) == m: coeffs = np.append(0., coeffs)
arrow = np.diag(np.append(0., supp) + 0.j)
arrow[0], arrow[1 :, 0] = coeffs, 1.
active = np.diag(np.append(0., np.ones(m)))
if return_vectors:
roots, vectors = eig(arrow, active)
else:
roots = eigvals(arrow, active)
badRoots = np.isinf(roots) + np.isnan(roots) == False
if return_vectors: return roots[badRoots], vectors[:, badRoots]
return roots[badRoots]
class RationalInterpolant(GenericStandardApproximant):
"""
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': kind of snapshots orthogonalization; allowed values
include 0, 1/2, and 1; defaults to 1, i.e. POD;
- 'scaleFactorDer': scaling factors for derivative computation;
defaults to 'AUTO';
- 'S': total number of samples current approximant relies upon;
- 'sampler': sample point generator;
- 'rationalMode': mode of rational approximation; allowed values
include 'MINIMAL[_STATE]', 'MINIMAL_OUTPUT',
'BARYCENTRIC[_STATE]', 'BARYCENTRIC_OUTPUT',
'STANDARD[_STATE]', and 'STANDARD_OUTPUT'; defaults to
'MINIMAL';
- 'polybasis': type of polynomial basis for interpolation; defaults
to 'MONOMIAL';
- 'polydegreetype': type of polynomial degree; allowed values are
'*' or '*_*', with * either 'TENSOR' or 'TOTAL'; defaults to
'TOTAL';
- 'M': degree of rational interpolant numerator; defaults to
'AUTO', i.e. maximum allowed;
- 'N': degree of rational interpolant denominator; defaults to
'AUTO', i.e. maximum allowed;
- 'MAuxiliary': degree of rational interpolant numerator with
respect to non-pivot parameters; defaults to 0;
- 'NAuxiliary': degree of rational interpolant denominator with
respect to non-pivot parameters; defaults to 0;
- 'functionalSolve': strategy for minimization of denominator
functional; allowed values include 'NORM' and 'DOMINANT';
defaults to 'NORM';
- 'interpTol': tolerance for interpolation; defaults to None;
- 'forceQReal': force denominator to have real coefficients;
defaults to False;
- 'polyTruncateTol': tolerance for truncation of rational terms;
defaults to 0;
- 'QTol': tolerance for robust rational denominator management;
defaults to 0.
Defaults to empty dict.
verbosity(optional): Verbosity level. Defaults to 10.
Attributes:
HFEngine: HF problem solver.
mu0: Default parameter.
mus: Array of snapshot parameters.
approxParameters: Dictionary containing values for main parameters of
approximant. Recognized keys are in parameterList.
parameterListSoft: Recognized keys of soft approximant parameters:
- 'POD': kind of snapshots orthogonalization;
- 'scaleFactorDer': scaling factors for derivative computation;
- 'rationalMode': mode of rational approximation;
- 'polybasis': type of polynomial basis for interpolation;
- 'polydegreetype': type of polynomial degree;
- 'M': degree of rational interpolant numerator;
- 'N': degree of rational interpolant denominator;
- 'MAuxiliary': degree of rational interpolant numerator with
respect to non-pivot parameters; defaults to 0;
- 'NAuxiliary': degree of rational interpolant denominator with
respect to non-pivot parameters; defaults to 0;
- 'functionalSolve': strategy for minimization of denominator
functional;
- 'interpTol': tolerance for interpolation via numpy.polyfit;
- 'forceQReal': force denominator to have real coefficients;
- 'polyTruncateTol': tolerance for truncation of rational terms;
- 'QTol': 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.
verbosity: Verbosity level.
POD: Kind of snapshots orthogonalization.
scaleFactorDer: Scaling factors for derivative computation.
S: Number of solution snapshots over which current approximant is
based upon.
sampler: Sample point generator.
rationalMode: Mode of rational approximation.
polybasis: Type of polynomial basis for interpolation.
polydegreetype: Type of polynomial degree.
M: Numerator degree of approximant.
N: Denominator degree of approximant.
MAuxiliary: Degree of rational interpolant numerator with respect to
non-pivot parameters.
NAuxiliary: Degree of rational interpolant denominator with respect to
non-pivot parameters.
functionalSolve: Strategy for minimization of denominator functional.
interpTol: Tolerance for interpolation via numpy.polyfit.
forceQReal: Force denominator to have real coefficients.
polyTruncateTol: Tolerance for truncation of rational terms.
QTol: Tolerance for robust rational denominator management.
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, *args, **kwargs):
self._preInit()
self._addParametersToList(["rationalMode", "polybasis", "M", "N",
"polydegreetype", "MAuxiliary",
"NAuxiliary", "functionalSolve",
"interpTol", "forceQReal",
"polyTruncateTol", "QTol"],
["MINIMAL", "MONOMIAL", "AUTO", "AUTO",
"TOTAL", 0, 0, "NORM", -1, 0, 0., 0.])
super().__init__(*args, **kwargs)
self._postInit()
@property
def tModelType(self):
from .trained_model.trained_model_rational import TrainedModelRational
return TrainedModelRational
@property
def rationalMode(self):
"""Value of rationalMode."""
return self._rationalMode
@rationalMode.setter
def rationalMode(self, rationalMode):
try:
rationalMode = rationalMode.upper().strip().replace(" ","")
if rationalMode in ["MINIMAL", "BARYCENTRIC", "STANDARD"]:
rationalMode += "_STATE"
if rationalMode not in ["MINIMAL_STATE", "MINIMAL_OUTPUT",
"BARYCENTRIC_STATE", "BARYCENTRIC_OUTPUT",
"STANDARD_STATE", "STANDARD_OUTPUT"]:
raise RROMPyException("Prescribed mode not recognized.")
except Exception as E:
RROMPyWarning(str(E) + " Overriding to 'MINIMAL'.")
rationalMode = "MINIMAL_STATE"
self._rationalMode = rationalMode
self._approxParameters["rationalMode"] = self.rationalMode
@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.")
except Exception as E:
RROMPyWarning(str(E) + " Overriding to 'MONOMIAL'.")
polybasis = "MONOMIAL"
self._polybasis = polybasis
self._approxParameters["polybasis"] = self.polybasis
@property
def functionalSolve(self):
"""Value of functionalSolve."""
return self._functionalSolve
@functionalSolve.setter
def functionalSolve(self, functionalSolve):
try:
functionalSolve = functionalSolve.upper().strip().replace(" ","")
if functionalSolve not in ["NORM", "DOMINANT"]:
raise RROMPyException(("Prescribed functionalSolve not "
"recognized."))
except Exception as E:
RROMPyWarning(str(E) + " Overriding to 'NORM'.")
functionalSolve = "NORM"
self._functionalSolve = functionalSolve
self._approxParameters["functionalSolve"] = self.functionalSolve
@property
def interpTol(self):
"""Value of interpTol."""
return self._interpTol
@interpTol.setter
def interpTol(self, interpTol):
self._interpTol = interpTol
self._approxParameters["interpTol"] = self.interpTol
@property
def M(self):
"""Value of M."""
return self._M
@M.setter
def M(self, M):
if isinstance(M, str):
M = M.strip().replace(" ","")
if "-" not in M: M = M + "-0"
self._M_isauto, self._M_shift = True, int(M.split("-")[-1])
M = 0
if M < 0: raise RROMPyException("M must be non-negative.")
self._M = M
self._approxParameters["M"] = self.M
def _setMAuto(self):
if self.rationalMode != "STANDARD_OUTPUT":
SEff = self.S
else:
dim = self.HFEngine.outputdim
if hasattr(self, "_N_isauto") and self._N_isauto:
M, N = reduceDegreeMN(self.S - 2, self.S - 1, self.S,
self.npar, self.polydegreetype,
self.MAuxiliary, self.NAuxiliary, dim, 0)
self.M = max(0, M - self._M_shift)
self.N = max(0, N - self._N_shift)
vbMng(self, "MAIN",
"Automatically setting M to {} and N to {}.".format(
self.M, self.N), 25)
return
degN = degreeToDOFs([self.N] + [self.NAuxiliary] * (self.npar - 1),
self.polydegreetypeM) - 1
SEff = int(np.floor(self.S - degN / dim))
M = reduceDegreeN(SEff - 1, SEff, self.npar, self.polydegreetypeM,
self.MAuxiliary)
self.M = max(0, M - self._M_shift)
vbMng(self, "MAIN", "Automatically setting M to {}.".format(self.M),
25)
@property
def N(self):
"""Value of N."""
return self._N
@N.setter
def N(self, N):
if isinstance(N, str):
N = N.strip().replace(" ","")
if "-" not in N: N = N + "-0"
self._N_isauto, self._N_shift = True, int(N.split("-")[-1])
N = 0
if N < 0: raise RROMPyException("N must be non-negative.")
self._N = N
self._approxParameters["N"] = self.N
def _setNAuto(self):
if self.rationalMode[-6:] == "OUTPUT":
dim = self.HFEngine.outputdim
else:
dim = self.HFEngine.spacedim
if (self.rationalMode[:7] == "MINIMAL"
or self.rationalMode[:11] == "BARYCENTRIC"):
SEff = min(self.S, dim + 1)
else:
if hasattr(self, "_M_isauto") and self._M_isauto:
if self.rationalMode == "STANDARD_OUTPUT":
return self._setMAuto()
M = reduceDegreeN(self.S - 2, self.S - 1, self.npar,
self.polydegreetypeM, self.MAuxiliary)
else:
M = reduceDegreeN(self.M, self.S - 1, self.npar,
self.polydegreetypeM, self.MAuxiliary)
degM = degreeToDOFs([M] + [self.MAuxiliary] * (self.npar - 1),
self.polydegreetypeN)
SEff = min(self.S, dim * (self.S - degM) + 1)
N = reduceDegreeN(SEff - 1, SEff, self.npar, self.polydegreetypeN,
self.NAuxiliary)
self.N = max(0, N - self._N_shift)
vbMng(self, "MAIN", "Automatically setting N to {}.".format(self.N),
25)
@property
def MAuxiliary(self):
"""Value of MAuxiliary."""
return self._MAuxiliary
@MAuxiliary.setter
def MAuxiliary(self, MAuxiliary):
if MAuxiliary < 0:
raise RROMPyException("MAuxiliary must be non-negative.")
self._MAuxiliary = MAuxiliary
self._approxParameters["MAuxiliary"] = self.MAuxiliary
@property
def NAuxiliary(self):
"""Value of NAuxiliary."""
return self._NAuxiliary
@NAuxiliary.setter
def NAuxiliary(self, NAuxiliary):
if NAuxiliary < 0:
raise RROMPyException("NAuxiliary must be non-negative.")
self._NAuxiliary = NAuxiliary
self._approxParameters["NAuxiliary"] = self.NAuxiliary
@property
def polydegreetype(self):
"""Value of polydegreetype."""
return self._polydegreetype
@polydegreetype.setter
def polydegreetype(self, polydegreetype):
try:
polydegreetype = polydegreetype.upper().strip().replace(" ","")
if "_" not in polydegreetype:
polydegreetype = "{0}_{0}".format(polydegreetype)
polydegreetypes = polydegreetype.split("_")
if not np.all([pdt in ["TOTAL", "TENSOR"]
for pdt in polydegreetypes]):
raise RROMPyException(("Prescribed polydegreetype not "
"recognized."))
except Exception as E:
RROMPyWarning(str(E) + " Overriding to 'TENSOR_TOTAL'.")
polydegreetype = "TENSOR_TOTAL"
self._polydegreetype = polydegreetype
self._approxParameters["polydegreetype"] = self.polydegreetype
@property
def polydegreetypeM(self): return self.polydegreetype.split("_")[0]
@property
def polydegreetypeN(self): return self.polydegreetype.split("_")[1]
@property
def forceQReal(self):
"""Value of forceQReal."""
return self._forceQReal
@forceQReal.setter
def forceQReal(self, forceQReal):
self._forceQReal = forceQReal
self._approxParameters["forceQReal"] = self.forceQReal
@property
def polyTruncateTol(self):
"""Value of tolerance for truncation of rational terms."""
return self._polyTruncateTol
@polyTruncateTol.setter
def polyTruncateTol(self, polyTruncateTol):
if polyTruncateTol < 0.:
RROMPyWarning(("Overriding prescribed negative numerator "
"tolerance to 0."))
polyTruncateTol = 0.
self._polyTruncateTol = polyTruncateTol
self._approxParameters["polyTruncateTol"] = self.polyTruncateTol
@property
def QTol(self):
"""Value of tolerance for robust rational denominator management."""
return self._QTol
@QTol.setter
def QTol(self, QTol):
if QTol < 0.:
RROMPyWarning(("Overriding prescribed negative robustness "
"tolerance to 0."))
QTol = 0.
self._QTol = QTol
self._approxParameters["QTol"] = self.QTol
def resetSamples(self):
"""Reset samples."""
super().resetSamples()
self._musUniqueCN = None
self._derIdxs = None
self._reorder = None
def _setupInterpolationIndices(self):
"""Setup parameters for polyvander."""
if self._musUniqueCN is None or len(self._reorder) != len(self.mus):
self._musUniqueCN, musIdxsTo, musIdxs, musCount = (
self.trainedModel.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)
if self.functionalSolve != "NORM" and self.npar > 1:
RROMPyWarning(("Strategy for functional optimization must be "
"'NORM' for more than one parameter. Overriding to "
"'NORM'."))
self.functionalSolve = "NORM"
if (self.rationalMode[:11] == "BARYCENTRIC" and self.npar > 1):
RROMPyWarning(("Minimal barycentric strategy cannot be applied "
"with more than one parameter. Overriding to "
"'MINIMAL'."))
self.rationalMode = "MINIMAL" + self.rationalMode[11:]
if (self.rationalMode[:11] == "BARYCENTRIC"
and not hasattr(self, "_N_isauto") and self.N + 1 < self.S):
RROMPyWarning(("Minimal barycentric strategy cannot be applied "
"with Least Squares. Overriding to 'MINIMAL'."))
self.rationalMode = "MINIMAL" + self.rationalMode[11:]
if self.rationalMode[:11] == "BARYCENTRIC":
self._setupInterpolationIndices()
if len(self._musUnique) != self.S:
RROMPyWarning(("Barycentric functional optimization cannot be "
"applied to repeated samples. Overriding to "
"'NORM'."))
self.rationalMode = "MINIMAL" + self.rationalMode[11:]
if hasattr(self, "_N_isauto"): self._setNAuto()
if self.rationalMode[:8] == "STANDARD" and hasattr(self, "_M_isauto"):
self._setMAuto()
if self.rationalMode[-6:] == "OUTPUT":
dim = self.HFEngine.outputdim
else:
dim = self.HFEngine.spacedim
if (self.rationalMode[:7] == "MINIMAL"
or self.rationalMode[:11] == "BARYCENTRIC"):
N = reduceDegreeN(self.N, min(self.S, dim + 1), self.npar,
self.polydegreetypeN, self.NAuxiliary)
if N < self.N:
RROMPyWarning(("N too large compared to S. Reducing N by "
"{}.").format(self.N - N))
self.N = N
else:
N = reduceDegreeMN(self.M, self.N, self.S, self.npar,
self.polydegreetype, self.MAuxiliary,
self.NAuxiliary, dim)[1]
if N < self.N:
RROMPyWarning(("N too large compared to S. Reducing N by "
"{}.").format(self.N - N))
self.N = N
if self.rationalMode[:11] != "BARYCENTRIC":
invD, TN = self._computeInterpolantInverseBlocks()
Rscaling = None
if self.rationalMode[-6:] == "OUTPUT":
sampleE = self.samplingEngine.samples_output
elif self.POD == 1:
sampleE = self.samplingEngine.Rscale
elif self.POD == 1/2:
sampleE = self.samplingEngine.samples_normal
Rscaling = self.samplingEngine.Rscale
else:
sampleE = self.samplingEngine.samples
while self.N > 0:
if self.rationalMode[:11] != "BARYCENTRIC":
NEff = degreeToDOFs([self.N]
+ [self.NAuxiliary] * (self.npar - 1),
self.polydegreetypeN)
TN = TN[:, : NEff]
if self.rationalMode[:7] == "MINIMAL":
ev, eV = self.findeveVGMinimal(sampleE, invD, TN, Rscaling)
elif self.rationalMode[:11] == "BARYCENTRIC":
ev, eV = self.findeveVGBarycentric(sampleE, Rscaling)
break
else: #if self.rationalMode[:8] == "STANDARD":
ev, eV = self.findeveVGStandard(sampleE, invD, TN, Rscaling)
nevBad = np.sum(np.abs(ev / ev[-1]) < self.QTol)
if not nevBad: break
dN = int(np.ceil(nevBad / (1 + self.NAuxiliary)
** (self.npar - 1)))
vbMng(self, "MAIN",
("Smallest {} eigenvalue{} below tolerance. Reducing N by "
"{}.").format(nevBad, "s" * (nevBad > 1), dN), 10)
self.N = self.N - dN
if self.rationalMode[-5:] == "STATE" and self.POD != 1: del self._gram
if self.N <= 0:
self.N = 0
eV = np.zeros((1,) + (self.NAuxiliary + 1,) * (self.npar - 1),
dtype = np.complex)
eV[(0,) * self.npar] = 1.
if self.N > 0 and self.rationalMode[:11] == "BARYCENTRIC":
q = PIN()
q.polybasis, q.nodes = self.polybasis, eV
else:
q = PI()
q.npar, q.polybasis = self.npar, self.polybasis
deg = [self.N + 1] + [self.NAuxiliary + 1] * (self.npar - 1)
if self.polydegreetypeN == "TENSOR":
q.coeffs = eV.reshape(deg)
else: #if self.polydegreetypeN == "TOTAL":
q.coeffs = mapTotalToTensor(deg, self.npar, eV)
q.truncate(self.polyTruncateTol)
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)
self._setupInterpolationIndices()
Qevaldiag = polyTimesTable(self.trainedModel.data.Q, self._musUniqueCN,
self._reorder, self._derIdxs,
self.scaleFactorRel)
if self.rationalMode[-6:] == "OUTPUT":
Qevaldiag = Qevaldiag.dot(self.samplingEngine.samples_output.T)
else:
if self.POD == 1:
Qevaldiag = Qevaldiag.dot(self.samplingEngine.Rscale.T)
elif self.POD == 1/2:
Qevaldiag = Qevaldiag * self.samplingEngine.Rscale
if hasattr(self, "_M_isauto"): self._setMAuto()
M = reduceDegreeN(self.M, self.S, self.npar, self.polydegreetypeM,
self.MAuxiliary)
if M < self.M:
RROMPyWarning(("M too large compared to S. Reducing M by "
"{}").format(self.M - M))
self.M = M
while self.M >= 0:
deg = [self.M] + [self.MAuxiliary] * (self.npar - 1)
pParRest = [deg, self.polybasis, self.verbosity >= 5,
self.polydegreetypeM, self.polyTruncateTol,
{"derIdxs": self._derIdxs, "reorder": self._reorder,
"scl": self.scaleFactorRel,
"degreetype": self.polydegreetypeM},
{"rcond": self.interpTol}]
p = PI()
wellCond, msg = p.setupByInterpolation(self._musUniqueCN,
Qevaldiag, *pParRest)
vbMng(self, "MAIN", msg, 5)
if wellCond: break
vbMng(self, "MAIN", ("Polyfit is poorly conditioned. Reducing M "
"by 1."), 10)
self.M = 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) -> int:
"""Compute rational interpolant."""
if self.checkComputedApprox(): return -1
RROMPyAssert(self._mode, message = "Cannot setup approximant.")
vbMng(self, "INIT", "Setting up {}.".format(self.name()), 5)
self.computeSnapshots()
if self.rationalMode[-6:] == "OUTPUT":
vbMng(self, "INIT", "Extracting system output from state.", 35)
self.samplingEngine.samples_output = self.HFEngine.applyC(
self.samplingEngine.samples, self.mus)
_POD, self._POD = self.POD, 1
vbMng(self, "DEL", "Done extracting system output.", 35)
self._setupTrainedModel(self.samplingEngine.projectionMatrix,
collapsed = (self.rationalMode[-6:] == "OUTPUT"))
self._setupRational(self._setupDenominator())
self.trainedModel.data.approxParameters = copy(self.approxParameters)
if self.rationalMode[-6:] == "OUTPUT": self._POD = _POD
vbMng(self, "DEL", "Done setting up approximant.", 5)
return 0
def _setupRational(self, Q:interpEng, P : interpEng = None):
vbMng(self, "INIT", "Starting approximant finalization.", 10)
self.trainedModel.data.Q = Q
if P is None: P = self._setupNumerator()
while self.N > 0 and self.npar == 1:
if self.HFEngine._ignoreResidues:
pls = Q.roots()
cfs, projMat = None, None
else:
cfs, pls, _ = rational2heaviside(P, Q)
cfs = cfs[: len(pls)].T
if self.POD != 1:
projMat = self.samplingEngine.projectionMatrix
else:
projMat = None
foci = self.sampler.normalFoci()
plsA = self.mapParameterList(self.mapParameterList(self.mu0)(0, 0)
+ self.scaleFactor * pls, "B")(0)
idxBad = self.HFEngine.flagBadPolesResiduesAbsolute(plsA, cfs,
projMat)
if not self.HFEngine._ignoreResidues: cfs[:, idxBad] = 0.
idxBad += self.HFEngine.flagBadPolesResiduesRelative(pls, cfs,
projMat, foci)
idxBad = idxBad > 0
if not np.any(idxBad): break
vbMng(self, "MAIN",
"Removing {} spurious pole{} out of {}.".format(
np.sum(idxBad), "s" * (np.sum(idxBad) > 1), self.N), 10)
if isinstance(Q, PIN):
Q.nodes = Q.nodes[idxBad == False]
else:
Q = PI()
Q.npar, Q.polybasis = 1, self.polybasis
Q.coeffs = np.ones(1, dtype = np.complex)
for pl in pls[idxBad == False]:
Q.coeffs = polyTimes(Q.coeffs, [- pl, 1.],
Pbasis = Q.polybasis, Rbasis = Q.polybasis)
Q.coeffs /= np.linalg.norm(Q.coeffs)
self.trainedModel.data.Q = Q
self.N = Q.deg[0]
P = self._setupNumerator()
self.trainedModel.data.P = P
vbMng(self, "DEL", "Terminated approximant finalization.", 10)
def _computeInterpolantInverseBlocks(self) -> Tuple[List[Np2D], Np2D]:
"""
Compute inverse factors for rational interpolant target functional.
"""
self._setupInterpolationIndices()
pvPPar = [self.polybasis, self._derIdxs, self._reorder,
self.scaleFactorRel, self.polydegreetypeN]
lagrange = self.npar == 1 and self.S == len(self._musUniqueCN)
if self.rationalMode[:7] == "MINIMAL" and not lagrange:
E = reduceDegreeN(self.S - 1, self.S, self.npar,
self.polydegreetypeN, self.NAuxiliary)
TN = None
deg = [self.NAuxiliary] * self.npar
while self.N >= 0:
if TN is None:
deg[0] = self.N
TE = TN = pvP(self._musUniqueCN, deg, *pvPPar)
if self.rationalMode[:7] == "MINIMAL" and not lagrange:
deg[0] = E
if self.N == E:
TE = TN
else:
TE = pvP(self._musUniqueCN, deg, *pvPPar)
fitOut = pseudoInverse(TE, rcond = self.interpTol,
full = True)
degE = deg[0] if self.npar == 1 else deg
vbMng(self, "MAIN",
("Fitting {} samples with degree {} through {}... "
"Conditioning of pseudoinverse system:"
"{:.4e}.").format(TE.shape[0], degE,
polyfitname(self.polybasis),
fitOut[1][1][0] / fitOut[1][1][-1]), 5)
if fitOut[1][0] == TE.shape[1]:
if self.rationalMode[:7] == "MINIMAL" and not lagrange:
deg0 = [E - 1] + [self.NAuxiliary] * (self.npar - 1)
S0 = degreeToDOFs(deg0, self.polydegreetypeN)
fitinv = fitOut[0][S0 :]
break
if self.rationalMode[:7] == "MINIMAL" and not lagrange:
E -= 1
if E >= self.N: continue
vbMng(self, "MAIN", ("Polyfit is poorly conditioned. Reducing "
"N by 1."), 10)
self.N, TN = self.N - 1, None
if self.N < 0:
raise RROMPyException(("Instability in computation of "
"denominator. Aborting."))
if self.rationalMode[:7] == "MINIMAL":
if lagrange:
mus = self._musUniqueCN[self._reorder]
dist = baseDistanceMatrix(mus, magnitude = False)[..., 0]
dist[np.arange(self.S), np.arange(self.S)] = 1.
fitinvE = np.prod(dist, axis = 1) ** -1
vbMng(self, "MAIN",
("Evaluating quasi-Lagrangian basis of degree {} at {} "
"sample points.").format(self.S - 1, self.S), 5)
invD = [np.diag(fitinvE)]
else:
invD = vanderInvTable(fitinv, self._reorder, self._derIdxs)
else: #if self.rationalMode[: 8] == "STANDARD":
if self.rationalMode[-6:] == "OUTPUT":
dim = self.HFEngine.outputdim
else:
dim = self.HFEngine.spacedim
degN = degreeToDOFs([self.N] + [self.NAuxiliary] * (self.npar - 1),
self.polydegreetypeN) - 1
SEff = int(np.floor(self.S - degN / dim))
MEff = reduceDegreeN(SEff - 1, SEff, self.npar,
self.polydegreetypeM, self.MAuxiliary)
deg = [self.MAuxiliary] * self.npar
while MEff >= 0:
deg[0] = MEff
if MEff == self.N and self.MAuxiliary == self.NAuxiliary:
TM = TN
else:
TM = pvP(self._musUniqueCN, deg, *pvPPar)
Q, R = np.linalg.qr(TM, 'complete')
fitOut = pseudoInverse(R[: R.shape[1]], rcond = self.interpTol,
full = True)
degM = deg[0] if self.npar == 1 else deg
vbMng(self, "MAIN",
("Fitting {} samples with degree {} through {}... "
"Conditioning of pseudoinverse system:"
"{:.4e}.").format(TM.shape[0], degM,
polyfitname(self.polybasis),
fitOut[1][1][0] / fitOut[1][1][-1]),
15)
if fitOut[1][0] == R.shape[1]: break
MEff -= 1
Qo = Q[:, TM.shape[1] :].T.conj()
vbMng(self, "MAIN",
("Loss of orthogonality for adjoint vandermonde: "
"{:.4e}.").format(np.linalg.norm(Qo.dot(TM))), 5)
invD = [Qo]
return invD, TN
def findeveVGMinimal(self, sampleE:Np2D, invD:List[Np2D], TN:Np2D,
Rscaling : Np1D = None) -> Tuple[Np1D, Np2D]:
"""
Compute eigenvalues and eigenvectors of rational denominator matrix, or
of its right chol factor if POD.
"""
RROMPyAssert(self._mode, message = "Cannot solve spectral problem.")
if self.POD == 1 or self.rationalMode[-6:] == "OUTPUT":
vbMng(self, "INIT", "Building generalized half-gramian.", 10)
S, eWidth = sampleE.shape[0], len(invD)
Rstack = np.zeros((S * eWidth, TN.shape[1]),
dtype = np.complex)
for k in range(eWidth):
Rstack[k * S : (k + 1) * S, :] = dot(sampleE,
dot(invD[k], TN))
vbMng(self, "DEL", "Done building half-gramian.", 10)
if self.forceQReal:
Rstack = np.vstack((np.real(Rstack), np.imag(Rstack)))
_, s, Vh = np.linalg.svd(Rstack,
full_matrices = Rstack.shape[0] < Rstack.shape[1])
evG, eVG = s[::-1], Vh[::-1].T.conj()
if len(evG) < eVG.shape[1]:
evG = np.append(np.zeros(eVG.shape[1] - len(evG)), evG)
evExp, probKind = -2., "svd "
else:
if not hasattr(self, "_gram"):
vbMng(self, "INIT", "Building gramian matrix.", 10)
self._gram = self.HFEngine.innerProduct(sampleE, sampleE,
is_state = True)
if Rscaling is not None:
self._gram = (self._gram.T * Rscaling.conj()).T * Rscaling
vbMng(self, "DEL", "Done building gramian.", 10)
vbMng(self, "INIT", "Building generalized gramian.", 10)
G = np.zeros((TN.shape[1],) * 2, dtype = np.complex)
for k in range(len(invD)):
iDkN = dot(invD[k], TN)
G += dot(dot(self._gram, iDkN).T, iDkN.conj()).T
vbMng(self, "DEL", "Done building gramian.", 10)
if self.forceQReal: G = np.real(G)
evG, eVG = np.linalg.eigh(G)
evExp, probKind = -1., "eigen"
if (self.functionalSolve == "NORM" or evG[0] == 0.
or np.sum(np.abs(evG) < np.finfo(float).eps * np.abs(evG[-1])
* len(evG)) == 1):
eV = eVG[:, 0]
else:
eV = eVG.dot(evG ** evExp * eVG[0].conj())
evG[1 :] -= evG[0]
eV /= np.linalg.norm(eV)
if self.npar == 1:
degE = self.N
else:
degE = [self.N] + [self.NAuxiliary] * (self.npar - 1)
cond = evG[-1] / evG[1] if evG[0] > 0 and evG[1] > 0 else np.inf
vbMng(self, "MAIN",
("Solved {}problem of size {} (degree {}) with condition number "
"{:.4e}.").format(probKind, len(evG) - 1, degE, cond), 5)
return evG[1 :], eV
def findeveVGBarycentric(self, sampleE:Np2D,
Rscaling : Np1D = None) -> Tuple[Np1D, Np2D]:
"""
Compute eigenvalues and eigenvectors of rational denominator matrix, or
of its right chol factor if POD.
"""
RROMPyAssert(self._mode, message = "Cannot solve spectral problem.")
if self.POD == 1 or self.rationalMode[-6:] == "OUTPUT":
Rstack = sampleE
if self.forceQReal:
Rstack = np.vstack((np.real(Rstack), np.imag(Rstack)))
_, s, Vh = np.linalg.svd(Rstack,
full_matrices = Rstack.shape[0] < Rstack.shape[1])
evG, eVG = s[::-1], Vh[::-1].T.conj()
if len(evG) < eVG.shape[1]:
evG = np.append(np.zeros(eVG.shape[1] - len(evG)), evG)
evExp, probKind = -2., "svd "
else:
if not hasattr(self, "_gram"):
vbMng(self, "INIT", "Building gramian matrix.", 10)
self._gram = self.HFEngine.innerProduct(sampleE, sampleE,
is_state = True)
if Rscaling is not None:
self._gram = (self._gram.T * Rscaling.conj()).T * Rscaling
vbMng(self, "DEL", "Done building gramian.", 10)
G = self._gram
if self.forceQReal: G = np.real(G)
evG, eVG = np.linalg.eigh(G)
evExp, probKind = -1., "eigen"
if (self.functionalSolve == "NORM" or evG[0] == 0.
or np.sum(np.abs(evG) < np.finfo(float).eps * np.abs(evG[-1])
* len(evG)) == 1):
eV = eVG[:, 0]
else:
eV = eVG.dot(evG ** evExp * np.sum(eVG, axis = 0).conj())
_evG = np.array(evG)
evG[1 :] -= evG[0]
eV /= np.linalg.norm(eV)
if self.npar == 1:
degE = self.N
else:
degE = [self.N] + [self.NAuxiliary] * (self.npar - 1)
cond = evG[-1] / evG[1] if evG[0] > 0 and evG[1] > 0 else np.inf
vbMng(self, "MAIN",
("Solved {}problem of size {} (degree {}) with condition number "
"{:.4e}.").format(probKind, len(evG) - 1, degE, cond), 5)
S, mus = len(eV), self._musUniqueCN[self._reorder].flatten()
poles, qTm1 = barycentricRoots(eV, mus, True)
self.N = len(poles)
if self.QTol > 0:
# compare optimal score with N poles to those obtained by
# removing one of the poles
qTm1 = qTm1[1 :].conj() ** -1.
dists = mus.reshape(-1, 1) - mus
dists[np.arange(S), np.arange(S)] = 1.
dists = np.prod(dists, axis = 1).conj() ** -1.
qComp = np.empty((self.N + 1, S), dtype = np.complex)
qComp[0] = dists * np.prod(qTm1, axis = 1)
for j in range(self.N):
qTmj = np.prod(qTm1[:, np.arange(self.N) != j], axis = 1)
qComp[j + 1] = dists * qTmj
Lqs = qComp.dot(eVG)
scores = np.real(np.sum(Lqs * _evG ** -evExp * Lqs.conj(),
axis = 1))
evBad = scores[1 :] < self.QTol * scores[0]
nevBad = np.sum(evBad)
if nevBad:
vbMng(self, "MAIN",
("Suboptimal pole{} detected. Reducing N by "
"{}.").format("s" * (nevBad > 1), nevBad), 10)
self.N = self.N - nevBad
poles = poles[evBad == False]
return evG[1 :], poles
def findeveVGStandard(self, sampleE:Np2D, invD:List[Np2D], TN:Np2D,
Rscaling : Np1D = None) -> Tuple[Np1D, Np2D]:
"""
Compute eigenvalues and eigenvectors of rational denominator matrix, or
of its right chol factor if POD.
"""
RROMPyAssert(self._mode, message = "Cannot solve spectral problem.")
if self.POD == 1 or self.rationalMode[-6:] == "OUTPUT":
vbMng(self, "INIT", "Building generalized half-gramian.", 10)
Rstack = khatri_rao(sampleE, invD[0]).dot(TN)
vbMng(self, "DEL", "Done building half-gramian.", 10)
if self.forceQReal:
Rstack = np.vstack((np.real(Rstack), np.imag(Rstack)))
_, s, Vh = np.linalg.svd(Rstack,
full_matrices = Rstack.shape[0] < Rstack.shape[1])
evG, eVG = s[::-1], Vh[::-1].T.conj()
if len(evG) < eVG.shape[1]:
evG = np.append(np.zeros(eVG.shape[1] - len(evG)), evG)
evExp, probKind = -2., "svd "
else:
if not hasattr(self, "_gram"):
vbMng(self, "INIT", "Building gramian matrix.", 10)
self._gram = self.HFEngine.innerProduct(sampleE, sampleE,
is_state = True)
if Rscaling is not None:
self._gram = (self._gram.T * Rscaling.conj()).T * Rscaling
vbMng(self, "DEL", "Done building gramian.", 10)
vbMng(self, "INIT", "Building generalized gramian.", 10)
G = self._gram * invD[0].T.conj().dot(invD[0])
G = TN.T.conj().dot(G.dot(TN))
vbMng(self, "DEL", "Done building gramian.", 10)
if self.forceQReal: G = np.real(G)
evG, eVG = np.linalg.eigh(G)
evExp, probKind = -1., "eigen"
if (self.functionalSolve == "NORM" or evG[0] == 0.
or np.sum(np.abs(evG) < np.finfo(float).eps * np.abs(evG[-1])
* len(evG)) == 1):
eV = eVG[:, 0]
else:
eV = eVG.dot(evG ** evExp * eVG[-1].conj())
evG[1 :] -= evG[0]
eV /= np.linalg.norm(eV)
if self.npar == 1:
degE = self.N
else:
degE = [self.N] + [self.NAuxiliary] * (self.npar - 1)
cond = evG[-1] / evG[1] if evG[0] > 0 and evG[1] > 0 else np.inf
vbMng(self, "MAIN",
("Solved {}problem of size {} (degree {}) with condition number "
"{:.4e}.").format(probKind, len(evG) - 1, degE, cond), 5)
return evG[1 :], eV
def getResidues(self, *args, **kwargs) -> Tuple[paramList, Np2D]:
"""
Obtain approximant residues.
Returns:
Matrix with residues as columns.
"""
return self.trainedModel.getResidues(*args, **kwargs)

Event Timeline

c4science · Help