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