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