diff --git a/rrompy/reduction_methods/standard/rational_interpolant.py b/rrompy/reduction_methods/standard/rational_interpolant.py
index d29ef58..d93444d 100644
--- a/rrompy/reduction_methods/standard/rational_interpolant.py
+++ b/rrompy/reduction_methods/standard/rational_interpolant.py
@@ -1,686 +1,720 @@
# 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_standard_approximant import GenericStandardApproximant
from rrompy.utilities.poly_fitting.polynomial import (
polybases as ppb, polyfitname,
polyvander as pvP, polyvanderTotal as pvTP,
polyTimes, polyTimesTable,
PolynomialInterpolator as PI)
from rrompy.utilities.poly_fitting.heaviside import rational2heaviside
from rrompy.utilities.poly_fitting.radial_basis import (polybases as rbpb,
RadialBasisInterpolator as RBI)
from rrompy.utilities.poly_fitting.moving_least_squares import (
polybases as mlspb,
MovingLeastSquaresInterpolator as MLSI)
from rrompy.utilities.base.types import (Np1D, Np2D, HFEng, DictAny, Tuple,
List, paramVal, sampList)
from rrompy.utilities.base import verbosityManager as vbMng
from rrompy.utilities.numerical import (customPInv, dot, fullDegreeN,
totalDegreeN, degreeTotalToFull,
fullDegreeMaxMask, totalDegreeMaxMask,
nextDerivativeIndices,
hashDerivativeToIdx as hashD)
from rrompy.utilities.exception_manager import (RROMPyException, RROMPyAssert,
RROMPyWarning)
__all__ = ['RationalInterpolant']
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': 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; defaults
to 'MONOMIAL';
- 'M': degree of rational interpolant numerator; defaults to 0;
- 'N': degree of rational interpolant denominator; defaults to 0;
- 'polydegreetype': type of polynomial degree; defaults to 'TOTAL';
- 'radialDirectionalWeights': radial basis weights for interpolant
numerator; defaults to 0, i.e. identity;
- 'nNearestNeighbor': mumber of nearest neighbors considered in
numerator if polybasis allows; defaults to -1;
- 'interpRcond': tolerance for interpolation; defaults to None;
- 'robustTol': tolerance for robust rational denominator
management; defaults to 0;
- 'centeredLike': whether samples should be managed as if centered;
involves making svd and interpolation problems square; defaults
to False.
Defaults to empty dict.
approx_state(optional): Whether to approximate state. Defaults to
False.
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': whether to compute POD of snapshots;
- 'polybasis': type of polynomial basis for interpolation;
- 'M': degree of rational interpolant numerator;
- 'N': degree of rational interpolant denominator;
- 'polydegreetype': type of polynomial degree;
- 'radialDirectionalWeights': radial basis weights for interpolant
numerator;
- 'nNearestNeighbor': mumber of nearest neighbors considered in
numerator if polybasis allows;
- 'interpRcond': tolerance for interpolation via numpy.polyfit;
- 'robustTol': tolerance for robust rational denominator
management;
- 'centeredLike': whether samples should be managed as if centered;
involves making svd and interpolation problems square.
parameterListCritical: Recognized keys of critical approximant
parameters:
- 'S': total number of samples current approximant relies upon;
- 'sampler': sample point generator.
approx_state: Whether to approximate state.
verbosity: Verbosity level.
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.
polydegreetype: Type of polynomial degree.
radialDirectionalWeights: Radial basis weights for interpolant
numerator.
nNearestNeighbor: Number of nearest neighbors considered in numerator
if polybasis allows.
interpRcond: Tolerance for interpolation via numpy.polyfit.
robustTol: Tolerance for robust rational denominator management.
centeredLike: Whether samples should be managed as if centered;
involves making svd and interpolation problems square.
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 = {}, approx_state : bool = False,
verbosity : int = 10, timestamp : bool = True):
self._preInit()
- self._addParametersToList(["polybasis", "M", "N", "polydegreetype",
+ self._addParametersToList(["polybasis", "M", "N", "polydegreetype",
"radialDirectionalWeights",
"nNearestNeighbor", "interpRcond",
- "robustTol", "centeredLike", "correctorTol",
+ "robustTol", "centeredLike",
+ "correctorKind", "correctorTol",
"correctorMaxIter"],
["MONOMIAL", 0, 0, "TOTAL", 1, -1, -1, 0,
- False, 0., 1])
+ False, "CONSERVATIVE", 0., 1])
super().__init__(HFEngine = HFEngine, mu0 = mu0,
approxParameters = approxParameters,
approx_state = approx_state, verbosity = verbosity,
timestamp = timestamp)
self.catchInstability = False
self._postInit()
@property
def tModelType(self):
from rrompy.reduction_methods.trained_model import TrainedModelRational
return TrainedModelRational
@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 ppb + rbpb + mlspb:
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 polybasis0(self):
if "_" in self.polybasis:
return self.polybasis.split("_")[0]
return self.polybasis
@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 radialDirectionalWeights(self):
"""Value of radialDirectionalWeights."""
return self._radialDirectionalWeights
@radialDirectionalWeights.setter
def radialDirectionalWeights(self, radialDirectionalWeights):
self._radialDirectionalWeights = radialDirectionalWeights
self._approxParameters["radialDirectionalWeights"] = (
self.radialDirectionalWeights)
@property
def nNearestNeighbor(self):
"""Value of nNearestNeighbor."""
return self._nNearestNeighbor
@nNearestNeighbor.setter
def nNearestNeighbor(self, nNearestNeighbor):
self._nNearestNeighbor = nNearestNeighbor
self._approxParameters["nNearestNeighbor"] = self.nNearestNeighbor
@property
def M(self):
"""Value of M."""
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
@property
def N(self):
"""Value of N."""
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
@property
def polydegreetype(self):
"""Value of polydegreetype."""
return self._polydegreetype
@polydegreetype.setter
def polydegreetype(self, polydegreetype):
try:
polydegreetype = polydegreetype.upper().strip().replace(" ","")
if polydegreetype not in ["TOTAL", "FULL"]:
raise RROMPyException(("Prescribed polydegreetype not "
"recognized."))
self._polydegreetype = polydegreetype
except:
RROMPyWarning(("Prescribed polydegreetype not recognized. "
"Overriding to 'TOTAL'."))
self._polydegreetype = "TOTAL"
self._approxParameters["polydegreetype"] = self.polydegreetype
@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 centeredLike(self):
"""Whether samples should be managed as if centered."""
return self._centeredLike
@centeredLike.setter
def centeredLike(self, centeredLike):
if centeredLike and not hasattr(self, "_centeredLike"):
RROMPyWarning(("Centered-like method is unstable for more than "
"one parameter."))
self._centeredLike = centeredLike
self._approxParameters["centeredLike"] = self.centeredLike
+ @property
+ def correctorKind(self):
+ """Value of correctorKind."""
+ return self._correctorKind
+ @correctorKind.setter
+ def correctorKind(self, correctorKind):
+ try:
+ correctorKind = correctorKind.upper().strip().replace(" ","")
+ if correctorKind not in ["CONSERVATIVE", "MINIMAL"]:
+ raise RROMPyException(("Prescribed correctorKind not "
+ "recognized."))
+ self._correctorKind = correctorKind
+ except:
+ RROMPyWarning(("Overriding prescribed corrector tolerance "
+ "to 'CONSERVATIVE'."))
+ self._correctorKind = "CONSERVATIVE"
+ self._approxParameters["correctorKind"] = self.correctorKind
+
@property
def correctorTol(self):
"""Value of correctorTol."""
return self._correctorTol
@correctorTol.setter
def correctorTol(self, correctorTol):
if correctorTol < 0. or (correctorTol > 0. and self.npar > 1):
RROMPyWarning(("Overriding prescribed corrector tolerance "
"to 0."))
correctorTol = 0.
self._correctorTol = correctorTol
self._approxParameters["correctorTol"] = self.correctorTol
@property
def correctorMaxIter(self):
"""Value of correctorMaxIter."""
return self._correctorMaxIter
@correctorMaxIter.setter
def correctorMaxIter(self, correctorMaxIter):
if correctorMaxIter < 1 or (correctorMaxIter > 1 and self.npar > 1):
RROMPyWarning(("Overriding prescribed max number of corrector "
"iterations to 1."))
correctorMaxIter = 1
self._correctorMaxIter = correctorMaxIter
self._approxParameters["correctorMaxIter"] = self.correctorMaxIter
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))
if self.centeredLike and len(self._musUniqueCN) > 1:
raise RROMPyException(("Cannot apply centered-like method "
"with more than one distinct sample "
"point."))
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, fitinv = self._computeInterpolantInverseBlocks()
if self.centeredLike:
if self.polydegreetype == "TOTAL":
Seff = totalDegreeN(self.N, self.npar)
else:
Seff = fullDegreeN(self.N, self.npar)
else:
Seff = self.S
idxSamplesEff = list(range(self.S - Seff, self.S))
if self.POD:
ev, eV = self.findeveVGQR(
self.samplingEngine.RPOD[:, idxSamplesEff], invD)
else:
ev, eV = self.findeveVGExplicit(
self.samplingEngine.samples(idxSamplesEff), invD)
nevBad = checkRobustTolerance(ev, self.robustTol)
if nevBad <= 1: break
if self.catchInstability:
raise RROMPyException(("Instability in denominator "
"computation: eigenproblem is poorly "
"conditioned."))
RROMPyWarning(("Smallest {} eigenvalues below tolerance. Reducing "
"N by 1.").format(nevBad))
self.N = self.N - 1
if self.N <= 0:
self._N = 0
eV = np.ones((1, 1))
q = PI()
q.npar = self.npar
q.polybasis = self.polybasis0
if self.polydegreetype == "TOTAL":
q.coeffs = degreeTotalToFull(tuple([self.N + 1] * self.npar),
self.npar, eV[:, 0])
else:
q.coeffs = eV[:, 0].reshape([self.N + 1] * self.npar)
vbMng(self, "DEL", "Done computing denominator.", 7)
return q, fitinv
def _setupNumerator(self):
"""Compute rational numerator."""
RROMPyAssert(self._mode, message = "Cannot setup numerator.")
vbMng(self, "INIT", "Starting computation of numerator.", 7)
verb = self.trainedModel.verbosity
self.trainedModel.verbosity = 0
self._setupInterpolationIndices()
Qevaldiag = polyTimesTable(self.trainedModel.data.Q, self._musUniqueCN,
self._reorder, self._derIdxs,
np.power(self.scaleFactor, -1.))
if self.POD:
Qevaldiag = Qevaldiag.dot(self.samplingEngine.RPOD.T)
self.trainedModel.verbosity = verb
cfun = totalDegreeN if self.polydegreetype == "TOTAL" else fullDegreeN
M = copy(self.M)
- while len(self.mus) < cfun(M, self.npar): M -= 1
+ while self.S < cfun(M, self.npar): M -= 1
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:
if self.centeredLike:
Seff = cfun(self.M, self.npar)
derIdxsEff = [self._derIdxs[0][: Seff]]
reorder = self._reorder[: Seff]
QevaldiagEff = Qevaldiag[: Seff, : Seff]
else:
derIdxsEff = self._derIdxs
reorder = self._reorder
QevaldiagEff = Qevaldiag
if self.polybasis in ppb:
p = PI()
wellCond, msg = p.setupByInterpolation(
self._musUniqueCN, QevaldiagEff, self.M,
self.polybasis, self.verbosity >= 5,
self.polydegreetype == "TOTAL",
{"derIdxs": derIdxsEff,
"reorder": reorder,
"scl": np.power(self.scaleFactor, -1.)},
{"rcond": self.interpRcond})
elif self.polybasis in rbpb:
p = RBI()
wellCond, msg = p.setupByInterpolation(
self._musUniqueCN, QevaldiagEff, self.M,
self.polybasis, self.radialDirectionalWeights,
self.verbosity >= 5, self.polydegreetype == "TOTAL",
{"derIdxs": derIdxsEff, "reorder": reorder,
"scl": np.power(self.scaleFactor, -1.),
"nNearestNeighbor": self.nNearestNeighbor},
{"rcond": self.interpRcond})
else:# if self.polybasis in mlspb:
p = MLSI()
wellCond, msg = p.setupByInterpolation(
self._musUniqueCN, QevaldiagEff, self.M,
self.polybasis, self.radialDirectionalWeights,
self.verbosity >= 5, self.polydegreetype == "TOTAL",
{"derIdxs": derIdxsEff, "reorder": reorder,
"scl": np.power(self.scaleFactor, -1.),
"nNearestNeighbor": self.nNearestNeighbor})
vbMng(self, "MAIN", msg, 5)
if wellCond: break
if self.catchInstability:
raise RROMPyException(("Instability in numerator computation: "
"polyfit is poorly conditioned."))
RROMPyWarning("Polyfit is poorly conditioned. Reducing M by 1.")
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):
"""Compute rational interpolant."""
if self.checkComputedApprox():
return
RROMPyAssert(self._mode, message = "Cannot setup approximant.")
vbMng(self, "INIT", "Setting up {}.". format(self.name()), 5)
self.computeSnapshots()
pMat = self.samplingEngine.samples.data
pMatEff = dot(self.HFEngine.C, pMat) if self.approx_state else pMat
if self.trainedModel is None:
self.trainedModel = self.tModelType()
self.trainedModel.verbosity = self.verbosity
self.trainedModel.timestamp = self.timestamp
datadict = {"mu0": self.mu0, "projMat": pMatEff,
"scaleFactor": self.scaleFactor,
"rescalingExp": self.HFEngine.rescalingExp}
self.trainedModel.data = self.initializeModelData(datadict)[0]
else:
self.trainedModel = self.trainedModel
self.trainedModel.data.projMat = copy(pMatEff)
- self._initCorrector()
+ self.trainedModel.data.mus = copy(self.mus)
+ self._iterCorrector()
+ self.trainedModel.data.approxParameters = copy(self.approxParameters)
+ vbMng(self, "DEL", "Done setting up approximant.", 5)
+
+ def _iterCorrector(self):
+ self._N0 = self.N
+ if self.correctorTol > 0. and (self.correctorMaxIter > 1
+ or self.correctorKind == "MINIMAL"):
+ vbMng(self, "INIT", "Starting correction iterations.", 15)
+ self._Qhat = PI()
+ self._Qhat.npar = self.npar
+ self._Qhat.polybasis = "MONOMIAL"
+ self._Qhat.coeffs = np.ones(1, dtype = np.complex)
+ if self.POD:
+ self._RPODOld = copy(self.samplingEngine.RPOD)
+ else:
+ self._samplesOld = copy(self.samplingEngine.samples)
for j in range(self.correctorMaxIter):
if self.N > 0:
Q = self._setupDenominator()[0]
else:
Q = PI()
- Q.coeffs = np.ones(tuple([1] * self.npar), dtype = np.complex)
+ Q.coeffs = np.ones((1,) * self.npar, dtype = np.complex)
Q.npar = self.npar
Q.polybasis = self.polybasis
- self.trainedModel.data.mus = copy(self.mus)
self.trainedModel.data.Q = Q
self.trainedModel.data.P = self._setupNumerator()
self._applyCorrector(j)
if self.N <= 0: break
- self._finalCorrector()
- self.trainedModel.data.approxParameters = copy(self.approxParameters)
- vbMng(self, "DEL", "Done setting up approximant.", 5)
-
- def _initCorrector(self):
- self._N0 = self.N
- if self.correctorTol <= 0. or self.correctorMaxIter <= 1: return
- vbMng(self, "INIT", "Starting correction iterations.", 15)
- self._Qhat = PI()
- self._Qhat.npar = self.npar
- self._Qhat.polybasis = "MONOMIAL"
- self._Qhat.coeffs = np.ones((1,) * self.npar, dtype = np.complex)
+ self._N = self._N0
+ del self._N0
+ if self.correctorTol <= 0. or (self.correctorMaxIter <= 1
+ and self.correctorKind == "CONSERVATIVE"):
+ return
if self.POD:
- self._RPODOld = copy(self.samplingEngine.RPOD)
+ self.samplingEngine.RPOD = self._RPODOld
+ del self._RPODOld
+ else:
+ self.samplingEngine.samples = self._samplesOld
+ del self._samplesOld
+ if self.correctorKind == "CONSERVATIVE":
+ QOld, QOldBasis = self.trainedModel.data.Q.coeffs, self.polybasis
+ else:
+ QOld, QOldBasis = [1.], "MONOMIAL"
+ Q = polyTimes(self._Qhat.coeffs, QOld, Pbasis = self._Qhat.polybasis,
+ Qbasis = QOldBasis, Rbasis = self.polybasis)
+ del self._Qhat
+ gamma = np.linalg.norm(Q)
+ self.trainedModel.data.Q.coeffs = np.pad(Q, (0, self._N - len(Q) + 1),
+ "constant") / gamma
+ if self.correctorKind == "CONSERVATIVE":
+ self.trainedModel.data.P.coeffs /= gamma
else:
- self._samplesOld = copy(self.samplingEngine.samples)
+ self.trainedModel.data.P = self._setupNumerator()
+ vbMng(self, "DEL", "Terminated correction iterations.", 15)
def _applyCorrector(self, j:int):
- if self.correctorTol <= 0. or j >= self.correctorMaxIter - 1:
+ if self.correctorTol <= 0. or (j >= self.correctorMaxIter - 1
+ and self.correctorKind == "CONSERVATIVE"):
self._N = 0
return
cfs, pls, _ = rational2heaviside(self.trainedModel.data.P,
self.trainedModel.data.Q)
- resEff = np.linalg.norm(cfs[: self.N], axis = 1)
+ cfs = cfs[: self.N]
+ if self.POD:
+ resEff = np.linalg.norm(cfs, axis = 1)
+ else:
+ resEff = self.HFEngine.norm(self.samplingEngine.samples.dot(cfs.T),
+ is_state = self.approx_state)
resEff /= np.max(np.hstack((np.ones((self.N, 1)),
np.abs(pls).reshape(-1, 1))), axis = 1)
resEff /= np.linalg.norm(resEff)
idxKeep = np.where(resEff >= self.correctorTol)[0]
vbMng(self, "MAIN",
("Correction iteration no. {}: {} out of {} residuals satisfy "
"tolerance.").format(j + 1, len(idxKeep), self.N), 15)
self._N -= len(idxKeep)
- if self.N <= 0: return
+ if self.N <= 0 and self.correctorKind == "CONSERVATIVE": return
for i in idxKeep:
self._Qhat.coeffs = polyTimes(self._Qhat.coeffs, [- pls[i], 1.],
Pbasis = self._Qhat.polybasis,
Rbasis = self._Qhat.polybasis)
self._Qhat.coeffs /= np.linalg.norm(self._Qhat.coeffs)
+ if self.N <= 0: return
+ vbMng(self, "MAIN",
+ ("Removing poles at (normalized positions): "
+ "{}.").format(pls[resEff < self.correctorTol]), 75)
That = polyTimesTable(self._Qhat, self._musUniqueCN,
self._reorder, self._derIdxs,
np.power(self.scaleFactor, -1.)).T
if self.POD:
self.samplingEngine.RPOD = self._RPODOld.dot(That)
else:
self.samplingEngine.samples = self._samplesOld.dot(That)
- def _finalCorrector(self):
- self._N = self._N0
- del self._N0
- if self.correctorTol <= 0. or self.correctorMaxIter <= 1: return
- if self.POD:
- self.samplingEngine.RPOD = self._RPODOld
- del self._RPODOld
- else:
- self.samplingEngine.samples = self._samplesOld
- del self._samplesOld
- Q = polyTimes(self.trainedModel.data.Q.coeffs, self._Qhat.coeffs,
- Pbasis = self.trainedModel.data.Q.polybasis,
- Qbasis = self._Qhat.polybasis,
- Rbasis = self.trainedModel.data.Q.polybasis)
- del self._Qhat
- gamma = np.linalg.norm(Q)
- self.trainedModel.data.P.coeffs /= gamma
- self.trainedModel.data.Q.coeffs = Q / gamma
- vbMng(self, "DEL", "Terminated correction iterations.", 15)
-
def _computeInterpolantInverseBlocks(self) -> Tuple[List[Np2D], Np2D]:
"""
Compute inverse factors for minimal interpolant target functional.
"""
RROMPyAssert(self._mode, message = "Cannot solve eigenvalue problem.")
self._setupInterpolationIndices()
cfun = totalDegreeN if self.polydegreetype == "TOTAL" else fullDegreeN
N = copy(self.N)
- while len(self.mus) < cfun(N, self.npar): N -= 1
+ while self.S < cfun(N, self.npar): N -= 1
if N < self.N:
RROMPyWarning(("N too large compared to S. Reducing N by "
"{}").format(self.N - N))
self.N = N
while self.N >= 0:
if self.centeredLike:
Seff = cfun(self.N, self.npar)
- #derIdxsEff = [self._derIdxs[0][- Seff :]]
derIdxsEff = [self._derIdxs[0][: Seff]]
reorder = self._reorder[: Seff]
else:
- Seff = len(self.mus)
+ Seff = self.S
derIdxsEff = self._derIdxs
reorder = self._reorder
if self.polydegreetype == "TOTAL":
TE = pvTP(self._musUniqueCN, self.N, self.polybasis0,
derIdxsEff, reorder,
scl = np.power(self.scaleFactor, -1.))
idxsB = totalDegreeMaxMask(self.N, self.npar)
else: #if self.polydegreetype == "FULL":
TE = pvP(self._musUniqueCN, [self.N] * self.npar,
self.polybasis0, derIdxsEff, reorder,
scl = np.power(self.scaleFactor, -1.))
idxsB = fullDegreeMaxMask(self.N, self.npar)
fitOut = customPInv(TE, rcond = self.interpRcond, full = True)
vbMng(self, "MAIN",
("Fitting {} samples with degree {} through {}... "
"Conditioning of pseudoinverse system: {:.4e}.").format(
TE.shape[0], self.N,
polyfitname(self.polybasis0),
fitOut[1][1][0] / fitOut[1][1][-1]),
5)
if fitOut[1][0] == TE.shape[1]:
fitinv = fitOut[0][idxsB, :]
break
if self.catchInstability:
raise RROMPyException(("Instability in denominator "
"computation: polyfit is poorly "
"conditioned."))
RROMPyWarning("Polyfit is poorly conditioned. Reducing N by 1.")
self.N = self.N - 1
if self.polydegreetype == "TOTAL":
TN = pvTP(self._musUniqueCN, self.N, self.polybasis0, derIdxsEff,
reorder, scl = np.power(self.scaleFactor, -1.))
else: #if self.polydegreetype == "FULL":
TN = pvP(self._musUniqueCN, [self.N] * self.npar, self.polybasis0,
derIdxsEff, reorder,
scl = np.power(self.scaleFactor, -1.))
invD = [None] * (len(idxsB))
for k in range(len(idxsB)):
pseudoInv = np.diag(fitinv[k, :])
idxGlob = 0
for j, derIdxs in enumerate(derIdxsEff):
nder = len(derIdxs)
idxGlob += nder
if nder > 1:
idxLoc = np.arange(Seff)[(reorder >= idxGlob - nder)
* (reorder < idxGlob)]
invLoc = 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] = dot(pseudoInv, TN)
return invD, fitinv
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,
is_state = self.approx_state)
G = np.zeros((nEnd, nEnd), dtype = np.complex)
for k in range(eWidth):
G += dot(dot(gramian, invD[k]).T, invD[k].conj()).T
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, :] = dot(RPODE, 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 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/standard/rational_moving_least_squares.py b/rrompy/reduction_methods/standard/rational_moving_least_squares.py
index 0491d84..8cd3013 100644
--- a/rrompy/reduction_methods/standard/rational_moving_least_squares.py
+++ b/rrompy/reduction_methods/standard/rational_moving_least_squares.py
@@ -1,330 +1,340 @@
# 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 import RationalInterpolant
from rrompy.utilities.poly_fitting.polynomial import (polybases as ppb,
polyvander as pvP,
polyvanderTotal as pvTP)
from rrompy.utilities.base.types import Np2D, HFEng, DictAny, paramVal
from rrompy.utilities.base import verbosityManager as vbMng
from rrompy.utilities.numerical import (fullDegreeMaxMask, totalDegreeMaxMask,
dot)
from rrompy.utilities.exception_manager import (RROMPyException, RROMPyAssert,
RROMPyWarning)
__all__ = ['RationalMovingLeastSquares']
class RationalMovingLeastSquares(RationalInterpolant):
"""
ROM rational moving LS 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; defaults
to 'MONOMIAL';
- 'M': degree of rational interpolant numerator; defaults to 0;
- 'N': degree of rational interpolant denominator; defaults to 0;
- 'polydegreetype': type of polynomial degree; defaults to 'TOTAL';
- 'radialBasis': numerator radial basis type; defaults to
'GAUSSIAN';
- 'radialDirectionalWeights': radial basis weights for interpolant
numerator; defaults to 0, i.e. identity;
- 'nNearestNeighbor': number of nearest neighbors considered in
numerator if radialBasis allows; defaults to -1;
- 'radialBasisDen': denominator radial basis type; defaults to
'GAUSSIAN';
- 'radialDirectionalWeightsDen': radial basis weights for
interpolant denominator; defaults to 0, i.e. identity;
- 'nNearestNeighborDen': number of nearest neighbors considered in
denominator if radialBasisDen allows; defaults to -1;
- 'interpRcond': tolerance for interpolation; defaults to None;
- 'robustTol': tolerance for robust rational denominator
management; defaults to 0.
Defaults to empty dict.
approx_state(optional): Whether to approximate state. Defaults to
False.
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': whether to compute POD of snapshots;
- 'polybasis': type of polynomial basis for interpolation;
- 'M': degree of rational interpolant numerator;
- 'N': degree of rational interpolant denominator;
- 'polydegreetype': type of polynomial degree;
- 'radialBasis': numerator radial basis type;
- 'radialDirectionalWeights': radial basis weights for interpolant
numerator;
- 'nNearestNeighbor': number of nearest neighbors considered in
numerator if radialBasis allows;
- 'radialBasisDen': denominator radial basis type;
- 'radialDirectionalWeightsDen': radial basis weights for
interpolant denominator;
- 'nNearestNeighborDen': number of nearest neighbors considered in
denominator if radialBasisDen allows;
- '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.
approx_state: Whether to approximate state.
verbosity: Verbosity level.
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.
polydegreetype: Type of polynomial degree.
radialBasis: Numerator radial basis type.
radialDirectionalWeights: Radial basis weights for interpolant
numerator.
nNearestNeighbor: Number of nearest neighbors considered in numerator
if radialBasis allows.
radialBasisDen: Denominator radial basis type.
radialDirectionalWeightsDen: Radial basis weights for interpolant
denominator.
nNearestNeighborDen: Number of nearest neighbors considered in
denominator if radialBasisDen allows.
interpRcond: Tolerance for interpolation via numpy.polyfit.
robustTol: 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, HFEngine:HFEng, mu0 : paramVal = None,
approxParameters : DictAny = {}, approx_state : bool = False,
verbosity : int = 10, timestamp : bool = True):
self._preInit()
self._addParametersToList(["radialBasis", "radialBasisDen",
"radialDirectionalWeightsDen",
"nNearestNeighborDen"],
["GAUSSIAN", "GAUSSIAN", 1, -1],
- toBeExcluded = ["correctorTol",
+ toBeExcluded = ["correctorKind",
+ "correctorTol",
"correctorMaxIter"])
super().__init__(HFEngine = HFEngine, mu0 = mu0,
approxParameters = approxParameters,
approx_state = approx_state, verbosity = verbosity,
timestamp = timestamp)
self.catchInstability = False
self._postInit()
+ @property
+ def correctorKind(self):
+ """Value of correctorKind."""
+ return "CONSERVATIVE"
+ @correctorKind.setter
+ def correctorKind(self, correctorKind):
+ RROMPyWarning(("correctorKind is used just to simplify inheritance, "
+ "and its value cannot be changed from 'CONSERVATIVE'."))
+
@property
def correctorTol(self):
"""Value of correctorTol."""
return 0.
@correctorTol.setter
def correctorTol(self, correctorTol):
RROMPyWarning(("correctorTol is used just to simplify inheritance, "
"and its value cannot be changed from 0."))
@property
def correctorMaxIter(self):
"""Value of correctorMaxIter."""
return 1
@correctorMaxIter.setter
def correctorMaxIter(self, correctorMaxIter):
RROMPyWarning(("correctorMaxIter is used just to simplify "
"inheritance, and its value cannot be changed from 1."))
@property
def tModelType(self):
from rrompy.reduction_methods.trained_model import \
TrainedModelRationalMLS
return TrainedModelRationalMLS
@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 ppb:
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):
self._radialBasis = radialBasis
self._approxParameters["radialBasis"] = self.radialBasis
@property
def radialBasisDen(self):
"""Value of radialBasisDen."""
return self._radialBasisDen
@radialBasisDen.setter
def radialBasisDen(self, radialBasisDen):
self._radialBasisDen = radialBasisDen
self._approxParameters["radialBasisDen"] = self.radialBasisDen
@property
def radialDirectionalWeightsDen(self):
"""Value of radialDirectionalWeightsDen."""
return self._radialDirectionalWeightsDen
@radialDirectionalWeightsDen.setter
def radialDirectionalWeightsDen(self, radialDirectionalWeightsDen):
self._radialDirectionalWeightsDen = radialDirectionalWeightsDen
self._approxParameters["radialDirectionalWeightsDen"] = (
self.radialDirectionalWeightsDen)
@property
def nNearestNeighborDen(self):
"""Value of nNearestNeighborDen."""
return self._nNearestNeighborDen
@nNearestNeighborDen.setter
def nNearestNeighborDen(self, nNearestNeighborDen):
self._nNearestNeighborDen = nNearestNeighborDen
self._approxParameters["nNearestNeighborDen"] = (
self.nNearestNeighborDen)
def _setupDenominator(self) -> Np2D:
"""Compute rational denominator."""
RROMPyAssert(self._mode, message = "Cannot setup denominator.")
vbMng(self, "INIT",
"Starting computation of denominator-related blocks.", 7)
self._setupInterpolationIndices()
if self.polydegreetype == "TOTAL":
TN = pvTP(self._musUniqueCN, self.N, self.polybasis0,
self._derIdxs, self._reorder,
scl = np.power(self.scaleFactor, -1.))
else: #if self.polydegreetype == "FULL":
TN = pvP(self._musUniqueCN, [self.N] * self.npar,
self.polybasis0, self._derIdxs, self._reorder,
scl = np.power(self.scaleFactor, -1.))
TNTen = np.zeros((self.S, self.S, TN.shape[1]), dtype = TN.dtype)
TNTen[np.arange(self.S), np.arange(self.S)] = TN
if self.POD: TNTen = dot(self.samplingEngine.RPOD, TNTen)
vbMng(self, "DEL", "Done computing denominator-related blocks.", 7)
return TN, TNTen
def _setupNumerator(self) -> Np2D:
"""Compute rational numerator."""
RROMPyAssert(self._mode, message = "Cannot setup numerator.")
vbMng(self, "INIT",
"Starting computation of denominator-related blocks.", 7)
self._setupInterpolationIndices()
if self.polydegreetype == "TOTAL":
TM = pvTP(self._musUniqueCN, self.M, self.polybasis0,
self._derIdxs, self._reorder,
scl = np.power(self.scaleFactor, -1.))
else: #if self.polydegreetype == "FULL":
TM = pvP(self._musUniqueCN, [self.M] * self.npar,
self.polybasis0, self._derIdxs, self._reorder,
scl = np.power(self.scaleFactor, -1.))
vbMng(self, "DEL", "Done computing denominator-related blocks.", 7)
return TM
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.computeSnapshots()
pMat = self.samplingEngine.samples.data
pMatEff = dot(self.HFEngine.C, pMat) if self.approx_state else pMat
if self.trainedModel is None:
self.trainedModel = self.tModelType()
self.trainedModel.verbosity = self.verbosity
self.trainedModel.timestamp = self.timestamp
datadict = {"mu0": self.mu0, "projMat": pMatEff,
"scaleFactor": self.scaleFactor,
"rescalingExp": self.HFEngine.rescalingExp}
data = self.initializeModelData(datadict)[0]
data.POD = self.POD
data.polybasis = self.polybasis
data.polydegreetype = self.polydegreetype
data.radialBasis = self.radialBasis
data.radialWeights = self.radialDirectionalWeights
data.nNearestNeighbor = self.nNearestNeighbor
data.radialBasisDen = self.radialBasisDen
data.radialWeightsDen = self.radialDirectionalWeightsDen
data.nNearestNeighborDen = self.nNearestNeighborDen
data.interpRcond = self.interpRcond
self.trainedModel.data = data
else:
self.trainedModel = self.trainedModel
self.trainedModel.data.projMat = copy(pMatEff)
if not self.POD:
self.trainedModel.data.gramian = self.HFEngine.innerProduct(
self.samplingEngine.samples,
self.samplingEngine.samples,
is_state = self.approx_state)
self.trainedModel.data.mus = copy(self.mus)
self.trainedModel.data.M = self.M
self.trainedModel.data.N = self.N
QVan, self.trainedModel.data.QBlocks = self._setupDenominator()
self.trainedModel.data.PVan = self._setupNumerator()
if self.polydegreetype == "TOTAL":
degreeMaxMask = totalDegreeMaxMask
else: #if self.polydegreetype == "FULL":
degreeMaxMask = fullDegreeMaxMask
if self.N > self.M:
self.trainedModel.data.QVan = QVan
self.trainedModel.data.domQIdxs = degreeMaxMask(self.N, self.npar)
else:
self.trainedModel.data.domQIdxs = degreeMaxMask(self.M, self.npar)
self.trainedModel.data.approxParameters = copy(self.approxParameters)
vbMng(self, "DEL", "Done setting up approximant.", 5)
diff --git a/rrompy/utilities/poly_fitting/polynomial/__init__.py b/rrompy/utilities/poly_fitting/polynomial/__init__.py
index 9710535..2d252cd 100644
--- a/rrompy/utilities/poly_fitting/polynomial/__init__.py
+++ b/rrompy/utilities/poly_fitting/polynomial/__init__.py
@@ -1,44 +1,46 @@
# Copyright (C) 2018 by the RROMPy authors
#
# This file is part of RROMPy.
#
# RROMPy is free software: you can redistribute it and/or modify
# it under the terms of the GNU Lesser General Public License as published by
# the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
#
# RROMPy is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU Lesser General Public License for more details.
#
# You should have received a copy of the GNU Lesser General Public License
# along with RROMPy. If not, see .
#
from .base import (polybases, polyfitname, polydomcoeff)
from .der import polyder
from .val import polyval
from .marginalize import polymarginalize
from .vander import polyvander, polyvanderTotal
from .roots import polyroots
-from .polynomial_algebra import changePolyBasis, polyTimes, polyTimesTable
+from .polynomial_algebra import (changePolyBasis, polyTimes, polyDivide,
+ polyTimesTable)
from .polynomial_interpolator import PolynomialInterpolator
__all__ = [
'polybases',
'polyfitname',
'polydomcoeff',
'polyder',
'polyval',
'polymarginalize',
'polyvander',
'polyvanderTotal',
'polyroots',
'changePolyBasis',
'polyTimes',
+ 'polyDivide',
'polyTimesTable',
'PolynomialInterpolator'
]
diff --git a/rrompy/utilities/poly_fitting/polynomial/polynomial_algebra.py b/rrompy/utilities/poly_fitting/polynomial/polynomial_algebra.py
index 9bd87a0..708de80 100644
--- a/rrompy/utilities/poly_fitting/polynomial/polynomial_algebra.py
+++ b/rrompy/utilities/poly_fitting/polynomial/polynomial_algebra.py
@@ -1,90 +1,117 @@
# Copyright (C) 2018 by the RROMPy authors
#
# This file is part of RROMPy.
#
# RROMPy is free software: you can redistribute it and/or modify
# it under the terms of the GNU Lesser General Public License as published by
# the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
#
# RROMPy is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU Lesser General Public License for more details.
#
# You should have received a copy of the GNU Lesser General Public License
# along with RROMPy. If not, see .
#
import numpy as np
from copy import deepcopy as copy
-from rrompy.utilities.base.types import Np1D, Np2D, List, interpEng
+from rrompy.utilities.base.types import Np1D, Np2D, Tuple, List, interpEng
from .vander import polyvander
from .polynomial_interpolator import PolynomialInterpolator
from rrompy.utilities.numerical import (multifactorial, customPInv,
hashDerivativeToIdx as hashD,
hashIdxToDerivative as hashI)
from rrompy.utilities.exception_manager import RROMPyException
-__all__ = ['changePolyBasis', 'polyTimes', 'polyTimesTable']
+__all__ = ['changePolyBasis', 'polyTimes', 'polyDivide', 'polyTimesTable']
def changePolyBasis(P:Np2D, dim : int = None, basis0 : str = "MONOMIAL",
basisF : str = "MONOMIAL") -> Np2D:
if basis0 == basisF: return P
if dim is None: dim = P.ndim
if basis0 != "MONOMIAL" and basisF != "MONOMIAL":
return changePolyBasis(changePolyBasis(P, dim, basis0, "MONOMIAL"),
dim, "MONOMIAL", basisF)
basisD = basisF if basis0 == "MONOMIAL" else basis0
R = copy(P)
N = np.max(P.shape[: dim]) - 1
vander = polyvander([0], N, basisD, [list(range(N + 1))])
if basis0 == "MONOMIAL": vander = customPInv(vander)
for j in range(dim):
R = np.tensordot(vander, R, (-1, j))
return R
def polyTimes(P:Np2D, Q:Np2D, dim : int = None, Pbasis : str = "MONOMIAL",
- Qbasis : str = "MONOMIAL", Rbasis : str = "MONOMIAL",
- allow_augmentation : bool = False) -> Np2D:
+ Qbasis : str = "MONOMIAL", Rbasis : str = "MONOMIAL") -> Np2D:
if not isinstance(P, (np.ndarray,)): P = np.array(P)
if not isinstance(Q, (np.ndarray,)): Q = np.array(Q)
P = changePolyBasis(P, dim, Pbasis, "MONOMIAL")
Q = changePolyBasis(Q, dim, Qbasis, "MONOMIAL")
if dim is None: dim = P.ndim
if dim <= 0: return
R = np.zeros([x + y - 1 for (x, y) in zip(P.shape[: dim], Q.shape[: dim])],
dtype = P.dtype)
if dim == 1:
for j, Qj in enumerate(Q):
R[j : j + len(P)] = R[j : j + len(P)] + Qj * P
else:
for j, Qj in enumerate(Q):
for l, Pl in enumerate(P):
R[j + l] = R[j + l] + polyTimes(Pl, Qj, dim - 1)
return changePolyBasis(R, dim, "MONOMIAL", Rbasis)
+def polyDivide(P:Np2D, Q:Np2D, dim : int = None, Pbasis : str = "MONOMIAL",
+ Qbasis : str = "MONOMIAL",
+ Rbasis : str = "MONOMIAL") -> Tuple[Np2D, Np2D]:
+ if not isinstance(P, (np.ndarray,)): P = np.array(P)
+ if not isinstance(Q, (np.ndarray,)): Q = np.array(Q)
+ P = changePolyBasis(P, dim, Pbasis, "MONOMIAL")
+ Pc = copy(P)
+ Q = changePolyBasis(Q, dim, Qbasis, "MONOMIAL")
+ if dim is None: dim = P.ndim
+ if dim <= 0: return
+ R = np.zeros([x - y + 1 for (x, y) in zip(P.shape[: dim], Q.shape[: dim])],
+ dtype = P.dtype)
+ if dim == 1:
+ for i in range(len(R) - 1, -1, -1):
+ try:
+ R[i] = Pc[-1] / Q[-1]
+ except:
+ raise RROMPyException(("Numerical instability in polynomial "
+ "quotient."))
+ Pc = Pc[: -1]
+ for j, Qj in enumerate(Q[::-1]):
+ if j > 0: Pc[-j] = Pc[-j] - R[i] * Qj
+ else:
+ raise RROMPyException(("Quotient of multivariate polynomials not "
+ "supported."))
+ return (changePolyBasis(R, dim, "MONOMIAL", Rbasis),
+ changePolyBasis(Pc, dim, "MONOMIAL", Rbasis))
+
def polyTimesTable(P:interpEng, mus:Np1D, reorder:List[int],
derIdxs:List[List[List[int]]], scl : Np1D = None) -> Np2D:
if not isinstance(P, PolynomialInterpolator):
raise RROMPyException(("Polynomial to evaluate must be a polynomial "
"interpolator."))
S = len(reorder)
T = np.zeros((S, S), dtype = np.complex)
idxGlob = 0
for j, derIdx in enumerate(derIdxs):
nder = len(derIdx)
idxLoc = np.arange(S)[np.logical_and(reorder >= idxGlob,
reorder < idxGlob + nder)]
idxGlob += nder
Pval = []
for der in range(nder):
derI = hashI(der, P.npar)
Pval += [P(mus[j], derI, scl) / multifactorial(derI)]
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)
T[idxLoc[derI], idxLoc[derJ]] = Pval[diffj]
return T