diff --git a/rrompy/reduction_methods/pivoted/greedy/rational_interpolant_greedy_pivoted_greedy.py b/rrompy/reduction_methods/pivoted/greedy/rational_interpolant_greedy_pivoted_greedy.py index b6a27cc..a33214d 100644 --- a/rrompy/reduction_methods/pivoted/greedy/rational_interpolant_greedy_pivoted_greedy.py +++ b/rrompy/reduction_methods/pivoted/greedy/rational_interpolant_greedy_pivoted_greedy.py @@ -1,519 +1,531 @@ #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 .generic_pivoted_greedy_approximant import ( GenericPivotedGreedyApproximantBase, GenericPivotedGreedyApproximantNoMatch, GenericPivotedGreedyApproximant) from rrompy.reduction_methods.standard.greedy import RationalInterpolantGreedy from rrompy.reduction_methods.standard.greedy.generic_greedy_approximant \ import pruneSamples from rrompy.reduction_methods.pivoted import ( RationalInterpolantGreedyPivotedNoMatch, RationalInterpolantGreedyPivoted) from rrompy.utilities.base.types import Np1D, Tuple, paramVal, paramList from rrompy.utilities.base import verbosityManager as vbMng from rrompy.utilities.exception_manager import RROMPyAssert from rrompy.parameter import emptyParameterList from rrompy.utilities.parallel import poolRank, recv __all__ = ['RationalInterpolantGreedyPivotedGreedyNoMatch', 'RationalInterpolantGreedyPivotedGreedy'] class RationalInterpolantGreedyPivotedGreedyBase( GenericPivotedGreedyApproximantBase): @property def sampleBatchSize(self): """Value of sampleBatchSize.""" return 1 @property def sampleBatchIdx(self): """Value of sampleBatchIdx.""" return self.S def greedyNextSample(self, muidx:int, plotEst : str = "NONE")\ -> Tuple[Np1D, int, float, paramVal]: """Compute next greedy snapshot of solution map.""" RROMPyAssert(self._mode, message = "Cannot add greedy sample.") mus = copy(self.muTest[muidx]) self.muTest.pop(muidx) for j, mu in enumerate(mus): vbMng(self, "MAIN", ("Adding sample point no. {} at {} to training " "set.").format(len(self.mus) + 1, mu), 3) self.mus.append(mu) self._S = len(self.mus) self._approxParameters["S"] = self.S if (self.samplingEngine.nsamples <= len(mus) - j - 1 or not np.allclose(mu, self.samplingEngine.mus[j - len(mus)])): self.samplingEngine.nextSample(mu) if self._isLastSampleCollinear(): vbMng(self, "MAIN", ("Collinearity above tolerance detected. Starting " "preemptive greedy loop termination."), 3) self._collinearityFlag = 1 errorEstTest = np.empty(len(self.muTest)) errorEstTest[:] = np.nan return errorEstTest, [-1], np.nan, np.nan errorEstTest, muidx, maxErrorEst = self.errorEstimator(self.muTest, True) if plotEst == "ALL": self.plotEstimator(errorEstTest, muidx, maxErrorEst) return errorEstTest, muidx, maxErrorEst, self.muTest[muidx] def _setSampleBatch(self, maxS:int): return self.S def _preliminaryTraining(self): """Initialize starting snapshots of solution map.""" RROMPyAssert(self._mode, message = "Cannot start greedy algorithm.") if self.samplingEngine.nsamples > 0: return self.resetSamples() self.samplingEngine.scaleFactor = self.scaleFactorDer musPivot = self.trainSetGenerator.generatePoints(self.S) while len(musPivot) > self.S: musPivot.pop() muTestBasePivot = self.samplerPivot.generatePoints(self.nTestPoints, False) idxPop = pruneSamples(self.HFEngine.mapParameterList(muTestBasePivot, idx = self.directionPivot), self.HFEngine.mapParameterList(musPivot, idx = self.directionPivot), 1e-10 * self.scaleFactorPivot[0]) muTestBasePivot.pop(idxPop) self._mus = emptyParameterList() self.mus.reset((self.S - 1, self.HFEngine.npar)) self.muTest = emptyParameterList() self.muTest.reset((len(muTestBasePivot) + 1, self.HFEngine.npar)) for k in range(self.S - 1): muk = np.empty_like(self.mus[0]) muk[self.directionPivot] = musPivot[k] muk[self.directionMarginal] = self.muMargLoc self.mus[k] = muk for k in range(len(muTestBasePivot)): muk = np.empty_like(self.muTest[0]) muk[self.directionPivot] = muTestBasePivot[k] muk[self.directionMarginal] = self.muMargLoc self.muTest[k] = muk muk = np.empty_like(self.mus[0]) muk[self.directionPivot] = musPivot[-1] muk[self.directionMarginal] = self.muMargLoc self.muTest[-1] = muk if len(self.mus) > 0: vbMng(self, "MAIN", ("Adding first {} sample point{} at {} to training " "set.").format(self.S - 1, "" + "s" * (self.S > 2), self.mus), 3) self.samplingEngine.iterSample(self.mus) self._S = len(self.mus) self._approxParameters["S"] = self.S self.M, self.N = ("AUTO",) * 2 def setupApproxPivoted(self, mus:paramList) -> int: if self.checkComputedApproxPivoted(): return -1 RROMPyAssert(self._mode, message = "Cannot setup approximant.") vbMng(self, "INIT", "Setting up pivoted approximant.", 10) if not hasattr(self, "_plotEstPivot"): self._plotEstPivot = "NONE" idx, sizes, emptyCores = self._preSetupApproxPivoted(mus) S0 = copy(self.S) pMat, Ps, Qs, req, musA = None, [], [], [], None if len(idx) == 0: vbMng(self, "MAIN", "Idling.", 45) if self.storeAllSamples: self.storeSamples() pL, pT, mT = recv(source = 0, tag = poolRank()) pMat = np.empty((pL, 0), dtype = pT) musA = np.empty((0, self.mu0.shape[1]), dtype = mT) else: for i in idx: self.muMargLoc = mus[i] vbMng(self, "MAIN", "Building marginal model no. {} at " "{}.".format(i + 1, self.muMargLoc), 25) self.samplingEngine.resetHistory() self.trainedModel = None self.verbosity -= 5 self.samplingEngine.verbosity -= 5 RationalInterpolantGreedy.setupApprox(self, self._plotEstPivot) self.verbosity += 5 self.samplingEngine.verbosity += 5 if self.storeAllSamples: self.storeSamples(i + self._nmusOld) pMat, req, musA = self._localPivotedResult(pMat, req, emptyCores, musA) Ps += [copy(self.trainedModel.data.P)] Qs += [copy(self.trainedModel.data.Q)] self._S = S0 del self.muMargLoc for r in req: r.wait() self._postSetupApproxPivoted(musA, pMat, Ps, Qs, sizes) vbMng(self, "DEL", "Done setting up pivoted approximant.", 10) return 0 def setupApprox(self, plotEst : str = "NONE") -> int: if self.checkComputedApprox(): return -1 if '_' not in plotEst: plotEst = plotEst + "_NONE" plotEstM, self._plotEstPivot = plotEst.split("_") val = super().setupApprox(plotEstM) return val class RationalInterpolantGreedyPivotedGreedyNoMatch( RationalInterpolantGreedyPivotedGreedyBase, GenericPivotedGreedyApproximantNoMatch, RationalInterpolantGreedyPivotedNoMatch): """ ROM greedy pivoted greedy rational interpolant computation for parametric problems (without pole matching). Args: HFEngine: HF problem solver. mu0(optional): Default parameter. Defaults to 0. directionPivot(optional): Pivot components. 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; - 'scaleFactorDer': scaling factors for derivative computation; defaults to 'AUTO'; - 'matchingWeightError': weight for pole matching optimization in error estimation; defaults to 0; - 'cutOffToleranceError': tolerance for ignoring parasitic poles in error estimation; defaults to 'AUTO', i.e. cutOffTolerance; - 'S': total number of pivot samples current approximant relies upon; - 'samplerPivot': pivot sample point generator; - 'SMarginal': number of starting marginal samples; - 'samplerMarginal': marginal sample point generator via sparse grid; - 'errorEstimatorKindMarginal': kind of marginal error estimator; available values include 'LEAVE_ONE_OUT', 'LOOK_AHEAD', and 'LOOK_AHEAD_RECOVER'; defaults to 'LEAVE_ONE_OUT'; - 'polybasis': type of polynomial basis for pivot interpolation; defaults to 'MONOMIAL'; - 'greedyTol': uniform error tolerance for greedy algorithm; defaults to 1e-2; - 'collinearityTol': collinearity tolerance for greedy algorithm; defaults to 0.; - 'maxIter': maximum number of greedy steps; defaults to 1e2; - 'nTestPoints': number of test points; defaults to 5e2; - 'trainSetGenerator': training sample points generator; defaults to sampler; - 'errorEstimatorKind': kind of error estimator; available values include 'AFFINE', 'DISCREPANCY', 'LOOK_AHEAD', 'LOOK_AHEAD_RES', and 'NONE'; defaults to 'NONE'; - 'greedyTolMarginal': uniform error tolerance for marginal greedy algorithm; defaults to 1e-1; - 'maxIterMarginal': maximum number of marginal greedy steps; defaults to 1e2; - 'radialDirectionalWeightsMarginal': radial basis weights for marginal interpolant; defaults to 1; + - 'functionalSolve': strategy for minimization of denominator + functional; allowed values include 'NORM' and 'DOMINANT'; + defaults to 'NORM'; - 'interpRcond': tolerance for pivot interpolation; defaults to None; - 'robustTol': tolerance for robust rational denominator management; defaults to 0; - 'cutOffTolerance': tolerance for ignoring parasitic poles; defaults to np.inf; - 'residueTol': tolerance for residue elimination; defaults to 0., i.e. no bad residues. 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. directionPivot: Pivot components. mus: Array of snapshot parameters. musMarginal: Array of marginal 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; - 'scaleFactorDer': scaling factors for derivative computation; - 'matchingWeightError': weight for pole matching optimization in error estimation; - 'cutOffToleranceError': tolerance for ignoring parasitic poles in error estimation; - 'errorEstimatorKindMarginal': kind of marginal error estimator; - 'polybasis': type of polynomial basis for pivot interpolation; - 'greedyTol': uniform error tolerance for greedy algorithm; - 'collinearityTol': collinearity tolerance for greedy algorithm; - 'maxIter': maximum number of greedy steps; - 'nTestPoints': number of test points; - 'trainSetGenerator': training sample points generator; - 'errorEstimatorKind': kind of error estimator; - 'greedyTolMarginal': uniform error tolerance for marginal greedy algorithm; - 'maxIterMarginal': maximum number of marginal greedy steps; - 'radialDirectionalWeightsMarginal': radial basis weights for marginal interpolant; + - 'functionalSolve': strategy for minimization of denominator + functional; - 'interpRcond': tolerance for pivot interpolation; - 'robustTol': tolerance for robust rational denominator management; - 'cutOffTolerance': tolerance for ignoring parasitic poles; - 'residueTol': tolerance for residue elimination. parameterListCritical: Recognized keys of critical approximant parameters: - 'S': total number of pivot samples current approximant relies upon; - 'samplerPivot': pivot sample point generator; - 'SMarginal': total number of marginal samples current approximant relies upon; - 'samplerMarginal': marginal sample point generator via sparse grid. approx_state: Whether to approximate state. verbosity: Verbosity level. POD: Whether to compute POD of snapshots. scaleFactorDer: Scaling factors for derivative computation. matchingWeightError: Weight for pole matching optimization in error estimation. cutOffToleranceError: Tolerance for ignoring parasitic poles in error estimation. S: Total number of pivot samples current approximant relies upon. samplerPivot: Pivot sample point generator. SMarginal: Total number of marginal samples current approximant relies upon. samplerMarginal: Marginal sample point generator via sparse grid. errorEstimatorKindMarginal: Kind of marginal error estimator. polybasis: Type of polynomial basis for pivot interpolation. greedyTol: uniform error tolerance for greedy algorithm. collinearityTol: Collinearity tolerance for greedy algorithm. maxIter: maximum number of greedy steps. nTestPoints: number of starting training points. trainSetGenerator: training sample points generator. errorEstimatorKind: kind of error estimator. greedyTolMarginal: Uniform error tolerance for marginal greedy algorithm. maxIterMarginal: Maximum number of marginal greedy steps. radialDirectionalWeightsMarginal: Radial basis weights for marginal interpolant. + functionalSolve: Strategy for minimization of denominator functional. interpRcond: Tolerance for pivot interpolation. robustTol: Tolerance for robust rational denominator management. cutOffTolerance: Tolerance for ignoring parasitic poles. residueTol: Tolerance for residue elimination. muBounds: list of bounds for pivot parameter values. muBoundsMarginal: list of bounds for marginal 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. """ class RationalInterpolantGreedyPivotedGreedy( RationalInterpolantGreedyPivotedGreedyBase, GenericPivotedGreedyApproximant, RationalInterpolantGreedyPivoted): """ ROM greedy pivoted greedy rational interpolant computation for parametric problems (with pole matching). Args: HFEngine: HF problem solver. mu0(optional): Default parameter. Defaults to 0. directionPivot(optional): Pivot components. 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; - 'scaleFactorDer': scaling factors for derivative computation; defaults to 'AUTO'; - 'matchingWeight': weight for pole matching optimization; defaults to 1; - 'matchingMode': mode for pole matching optimization; allowed values include 'NONE' and 'SHIFT'; defaults to 'NONE'; - 'sharedRatio': required ratio of marginal points to share resonance; defaults to 1.; - 'matchingWeightError': weight for pole matching optimization in error estimation; defaults to 0; - 'cutOffToleranceError': tolerance for ignoring parasitic poles in error estimation; defaults to 'AUTO', i.e. cutOffTolerance; - 'S': total number of pivot samples current approximant relies upon; - 'samplerPivot': pivot sample point generator; - 'SMarginal': number of starting marginal samples; - 'samplerMarginal': marginal sample point generator via sparse grid; - 'errorEstimatorKindMarginal': kind of marginal error estimator; available values include 'LEAVE_ONE_OUT', 'LOOK_AHEAD', and 'LOOK_AHEAD_RECOVER'; defaults to 'LEAVE_ONE_OUT'; - 'polybasis': type of polynomial basis for pivot interpolation; defaults to 'MONOMIAL'; - 'polybasisMarginal': type of polynomial basis for marginal interpolation; allowed values include 'MONOMIAL_*', 'CHEBYSHEV_*', 'LEGENDRE_*', 'NEARESTNEIGHBOR', and 'PIECEWISE_LINEAR_*'; defaults to 'MONOMIAL'; - 'paramsMarginal': dictionary of parameters for marginal interpolation; include: . 'MMarginal': degree of marginal interpolant; defaults to 'AUTO', i.e. maximum allowed; not for 'NEARESTNEIGHBOR' or 'PIECEWISE_LINEAR_*'; . 'nNeighborsMarginal': number of marginal nearest neighbors; defaults to 1; only for 'NEARESTNEIGHBOR'; . 'polydegreetypeMarginal': type of polynomial degree for marginal; defaults to 'TOTAL'; not for 'NEARESTNEIGHBOR' or 'PIECEWISE_LINEAR_*'; . 'interpRcondMarginal': tolerance for marginal interpolation; defaults to None; not for 'NEARESTNEIGHBOR'; . 'radialDirectionalWeightsMarginalAdapt': bounds for adaptive rescaling of marginal radial basis weights; only for radial basis. - 'greedyTol': uniform error tolerance for greedy algorithm; defaults to 1e-2; - 'collinearityTol': collinearity tolerance for greedy algorithm; defaults to 0.; - 'maxIter': maximum number of greedy steps; defaults to 1e2; - 'nTestPoints': number of test points; defaults to 5e2; - 'trainSetGenerator': training sample points generator; defaults to sampler; - 'errorEstimatorKind': kind of error estimator; available values include 'AFFINE', 'DISCREPANCY', 'LOOK_AHEAD', 'LOOK_AHEAD_RES', and 'NONE'; defaults to 'NONE'; - 'greedyTolMarginal': uniform error tolerance for marginal greedy algorithm; defaults to 1e-1; - 'maxIterMarginal': maximum number of marginal greedy steps; defaults to 1e2; - 'radialDirectionalWeightsMarginal': radial basis weights for marginal interpolant; defaults to 1; + - 'functionalSolve': strategy for minimization of denominator + functional; allowed values include 'NORM' and 'DOMINANT'; + defaults to 'NORM'; - 'interpRcond': tolerance for pivot interpolation; defaults to None; - 'robustTol': tolerance for robust rational denominator management; defaults to 0; - 'cutOffTolerance': tolerance for ignoring parasitic poles; defaults to np.inf; - 'residueTol': tolerance for residue elimination; defaults to 0., i.e. no bad residues. 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. directionPivot: Pivot components. mus: Array of snapshot parameters. musMarginal: Array of marginal 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; - 'scaleFactorDer': scaling factors for derivative computation; - 'matchingWeight': weight for pole matching optimization; - 'matchingMode': mode for pole matching optimization; - 'sharedRatio': required ratio of marginal points to share resonance; - 'matchingWeightError': weight for pole matching optimization in error estimation; - 'cutOffToleranceError': tolerance for ignoring parasitic poles in error estimation; - 'errorEstimatorKindMarginal': kind of marginal error estimator; - 'polybasis': type of polynomial basis for pivot interpolation; - 'polybasisMarginal': type of polynomial basis for marginal interpolation; - 'paramsMarginal': dictionary of parameters for marginal interpolation; include: . 'MMarginal': degree of marginal interpolant; . 'nNeighborsMarginal': number of marginal nearest neighbors; . 'polydegreetypeMarginal': type of polynomial degree for marginal; . 'interpRcondMarginal': tolerance for marginal interpolation; . 'radialDirectionalWeightsMarginalAdapt': bounds for adaptive rescaling of marginal radial basis weights. - 'greedyTol': uniform error tolerance for greedy algorithm; - 'collinearityTol': collinearity tolerance for greedy algorithm; - 'maxIter': maximum number of greedy steps; - 'nTestPoints': number of test points; - 'trainSetGenerator': training sample points generator; - 'errorEstimatorKind': kind of error estimator; - 'greedyTolMarginal': uniform error tolerance for marginal greedy algorithm; - 'maxIterMarginal': maximum number of marginal greedy steps; - 'radialDirectionalWeightsMarginal': radial basis weights for marginal interpolant; + - 'functionalSolve': strategy for minimization of denominator + functional; - 'interpRcond': tolerance for pivot interpolation; - 'robustTol': tolerance for robust rational denominator management; - 'cutOffTolerance': tolerance for ignoring parasitic poles; - 'residueTol': tolerance for residue elimination. parameterListCritical: Recognized keys of critical approximant parameters: - 'S': total number of pivot samples current approximant relies upon; - 'samplerPivot': pivot sample point generator; - 'SMarginal': total number of marginal samples current approximant relies upon; - 'samplerMarginal': marginal sample point generator via sparse grid. approx_state: Whether to approximate state. verbosity: Verbosity level. POD: Whether to compute POD of snapshots. scaleFactorDer: Scaling factors for derivative computation. matchingWeight: Weight for pole matching optimization. matchingMode: Mode for pole matching optimization. sharedRatio: Required ratio of marginal points to share resonance. matchingWeightError: Weight for pole matching optimization in error estimation. cutOffToleranceError: Tolerance for ignoring parasitic poles in error estimation. S: Total number of pivot samples current approximant relies upon. samplerPivot: Pivot sample point generator. SMarginal: Total number of marginal samples current approximant relies upon. samplerMarginal: Marginal sample point generator via sparse grid. errorEstimatorKindMarginal: Kind of marginal error estimator. polybasis: Type of polynomial basis for pivot interpolation. polybasisMarginal: Type of polynomial basis for marginal interpolation. paramsMarginal: Dictionary of parameters for marginal interpolation. greedyTol: uniform error tolerance for greedy algorithm. collinearityTol: Collinearity tolerance for greedy algorithm. maxIter: maximum number of greedy steps. nTestPoints: number of starting training points. trainSetGenerator: training sample points generator. errorEstimatorKind: kind of error estimator. greedyTolMarginal: Uniform error tolerance for marginal greedy algorithm. maxIterMarginal: Maximum number of marginal greedy steps. radialDirectionalWeightsMarginal: Radial basis weights for marginal interpolant. + functionalSolve: Strategy for minimization of denominator functional. interpRcond: Tolerance for pivot interpolation. robustTol: Tolerance for robust rational denominator management. cutOffTolerance: Tolerance for ignoring parasitic poles. residueTol: Tolerance for residue elimination. muBounds: list of bounds for pivot parameter values. muBoundsMarginal: list of bounds for marginal 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. """ diff --git a/rrompy/reduction_methods/pivoted/greedy/rational_interpolant_pivoted_greedy.py b/rrompy/reduction_methods/pivoted/greedy/rational_interpolant_pivoted_greedy.py index 26ad0ea..90c65c6 100644 --- a/rrompy/reduction_methods/pivoted/greedy/rational_interpolant_pivoted_greedy.py +++ b/rrompy/reduction_methods/pivoted/greedy/rational_interpolant_pivoted_greedy.py @@ -1,446 +1,458 @@ # 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 from numpy import empty, empty_like from .generic_pivoted_greedy_approximant import ( GenericPivotedGreedyApproximantBase, GenericPivotedGreedyApproximantNoMatch, GenericPivotedGreedyApproximant) from rrompy.reduction_methods.standard import RationalInterpolant from rrompy.reduction_methods.pivoted import ( RationalInterpolantPivotedNoMatch, RationalInterpolantPivoted) from rrompy.utilities.base.types import paramList from rrompy.utilities.base import verbosityManager as vbMng from rrompy.utilities.exception_manager import RROMPyAssert from rrompy.parameter import emptyParameterList from rrompy.utilities.parallel import poolRank, recv __all__ = ['RationalInterpolantPivotedGreedyNoMatch', 'RationalInterpolantPivotedGreedy'] class RationalInterpolantPivotedGreedyBase( GenericPivotedGreedyApproximantBase): def computeSnapshots(self): """Compute snapshots of solution map.""" RROMPyAssert(self._mode, message = "Cannot start snapshot computation.") vbMng(self, "INIT", "Starting computation of snapshots.", 5) self.samplingEngine.scaleFactor = self.scaleFactorDer if not hasattr(self, "musPivot") or len(self.musPivot) != self.S: self.musPivot = self.samplerPivot.generatePoints(self.S) while len(self.musPivot) > self.S: self.musPivot.pop() musLoc = emptyParameterList() musLoc.reset((self.S, self.HFEngine.npar)) self.samplingEngine.resetHistory() for k in range(self.S): muk = empty_like(musLoc[0]) muk[self.directionPivot] = self.musPivot[k] muk[self.directionMarginal] = self.muMargLoc musLoc[k] = muk self.samplingEngine.iterSample(musLoc) vbMng(self, "DEL", "Done computing snapshots.", 5) self._m_selfmus = copy(musLoc) self._mus = self.musPivot self._m_mu0 = copy(self.mu0) self._m_HFEparameterMap = copy(self.HFEngine.parameterMap) self._mu0 = self.checkParameterListPivot(self.mu0(self.directionPivot)) self.HFEngine.parameterMap = { "F": [self.HFEngine.parameterMap["F"][self.directionPivot[0]]], "B": [self.HFEngine.parameterMap["B"][self.directionPivot[0]]]} def setupApproxPivoted(self, mus:paramList) -> int: if self.checkComputedApproxPivoted(): return -1 RROMPyAssert(self._mode, message = "Cannot setup approximant.") vbMng(self, "INIT", "Setting up pivoted approximant.", 10) idx, sizes, emptyCores = self._preSetupApproxPivoted(mus) pMat, Ps, Qs, req, musA = None, [], [], [], None if len(idx) == 0: vbMng(self, "MAIN", "Idling.", 45) if self.storeAllSamples: self.storeSamples() pL, pT, mT = recv(source = 0, tag = poolRank()) pMat = empty((pL, 0), dtype = pT) musA = empty((0, self.mu0.shape[1]), dtype = mT) else: for i in idx: self.muMargLoc = mus[i] vbMng(self, "MAIN", "Building marginal model no. {} at " "{}.".format(i + 1, self.muMargLoc), 25) self.samplingEngine.resetHistory() self.trainedModel = None self.verbosity -= 5 self.samplingEngine.verbosity -= 5 RationalInterpolant.setupApprox(self) self.verbosity += 5 self.samplingEngine.verbosity += 5 self._mu0 = self._m_mu0 self._mus = self._m_selfmus self.HFEngine.parameterMap = self._m_HFEparameterMap del self._m_mu0, self._m_selfmus, self._m_HFEparameterMap if self.storeAllSamples: self.storeSamples(i + self._nmusOld) pMat, req, musA = self._localPivotedResult(pMat, req, emptyCores, musA) Ps += [copy(self.trainedModel.data.P)] Qs += [copy(self.trainedModel.data.Q)] del self.muMargLoc for r in req: r.wait() self._postSetupApproxPivoted(musA, pMat, Ps, Qs, sizes) vbMng(self, "DEL", "Done setting up pivoted approximant.", 10) return 0 class RationalInterpolantPivotedGreedyNoMatch( RationalInterpolantPivotedGreedyBase, GenericPivotedGreedyApproximantNoMatch, RationalInterpolantPivotedNoMatch): """ ROM pivoted greedy rational interpolant computation for parametric problems (without pole matching). Args: HFEngine: HF problem solver. mu0(optional): Default parameter. Defaults to 0. directionPivot(optional): Pivot components. 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; - 'scaleFactorDer': scaling factors for derivative computation; defaults to 'AUTO'; - 'matchingWeightError': weight for pole matching optimization in error estimation; defaults to 0; - 'cutOffToleranceError': tolerance for ignoring parasitic poles in error estimation; defaults to 'AUTO', i.e. cutOffTolerance; - 'S': total number of pivot samples current approximant relies upon; - 'samplerPivot': pivot sample point generator; - 'SMarginal': number of starting marginal samples; - 'samplerMarginal': marginal sample point generator via sparse grid; - 'errorEstimatorKindMarginal': kind of marginal error estimator; available values include 'LEAVE_ONE_OUT', 'LOOK_AHEAD', and 'LOOK_AHEAD_RECOVER'; defaults to 'LEAVE_ONE_OUT'; - 'polybasis': type of polynomial basis for pivot interpolation; defaults to 'MONOMIAL'; - '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; - 'greedyTolMarginal': uniform error tolerance for marginal greedy algorithm; defaults to 1e-1; - 'maxIterMarginal': maximum number of marginal greedy steps; defaults to 1e2; - 'radialDirectionalWeights': radial basis weights for pivot numerator; defaults to 1; - 'radialDirectionalWeightsAdapt': bounds for adaptive rescaling of radial basis weights; defaults to [-1, -1]; - 'radialDirectionalWeightsMarginal': radial basis weights for marginal interpolant; defaults to 1; + - 'functionalSolve': strategy for minimization of denominator + functional; allowed values include 'NORM' and 'DOMINANT'; + defaults to 'NORM'; - 'interpRcond': tolerance for pivot interpolation; defaults to None; - 'robustTol': tolerance for robust rational denominator management; defaults to 0; - 'cutOffTolerance': tolerance for ignoring parasitic poles; defaults to np.inf; - 'residueTol': tolerance for residue elimination; defaults to 0., i.e. no bad residues. 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. directionPivot: Pivot components. mus: Array of snapshot parameters. musPivot: Array of pivot snapshot parameters. musMarginal: Array of marginal 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; - 'scaleFactorDer': scaling factors for derivative computation; - 'matchingWeightError': weight for pole matching optimization in error estimation; - 'cutOffToleranceError': tolerance for ignoring parasitic poles in error estimation; - 'errorEstimatorKindMarginal': kind of marginal error estimator; - 'polybasis': type of polynomial basis for pivot interpolation; - 'M': degree of rational interpolant numerator; - 'N': degree of rational interpolant denominator; - 'greedyTolMarginal': uniform error tolerance for marginal greedy algorithm; - 'maxIterMarginal': maximum number of marginal greedy steps; - 'radialDirectionalWeights': radial basis weights for pivot numerator; - 'radialDirectionalWeightsAdapt': bounds for adaptive rescaling of radial basis weights; - 'radialDirectionalWeightsMarginal': radial basis weights for marginal interpolant; + - 'functionalSolve': strategy for minimization of denominator + functional; - 'interpRcond': tolerance for pivot interpolation; - 'robustTol': tolerance for robust rational denominator management; - 'cutOffTolerance': tolerance for ignoring parasitic poles; - 'residueTol': tolerance for residue elimination. parameterListCritical: Recognized keys of critical approximant parameters: - 'S': total number of pivot samples current approximant relies upon; - 'samplerPivot': pivot sample point generator; - 'SMarginal': total number of marginal samples current approximant relies upon; - 'samplerMarginal': marginal sample point generator via sparse grid. approx_state: Whether to approximate state. verbosity: Verbosity level. POD: Whether to compute POD of snapshots. scaleFactorDer: Scaling factors for derivative computation. matchingWeightError: Weight for pole matching optimization in error estimation. cutOffToleranceError: Tolerance for ignoring parasitic poles in error estimation. S: Total number of pivot samples current approximant relies upon. samplerPivot: Pivot sample point generator. SMarginal: Total number of marginal samples current approximant relies upon. samplerMarginal: Marginal sample point generator via sparse grid. errorEstimatorKindMarginal: Kind of marginal error estimator. polybasis: Type of polynomial basis for pivot interpolation. M: Degree of rational interpolant numerator. N: Degree of rational interpolant denominator. greedyTolMarginal: Uniform error tolerance for marginal greedy algorithm. maxIterMarginal: Maximum number of marginal greedy steps. radialDirectionalWeights: Radial basis weights for pivot numerator. radialDirectionalWeightsAdapt: Bounds for adaptive rescaling of radial basis weights. radialDirectionalWeightsMarginal: Radial basis weights for marginal interpolant. + functionalSolve: Strategy for minimization of denominator functional. interpRcond: Tolerance for pivot interpolation. robustTol: Tolerance for robust rational denominator management. cutOffTolerance: Tolerance for ignoring parasitic poles. residueTol: Tolerance for residue elimination. muBounds: list of bounds for pivot parameter values. muBoundsMarginal: list of bounds for marginal 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. """ class RationalInterpolantPivotedGreedy(RationalInterpolantPivotedGreedyBase, GenericPivotedGreedyApproximant, RationalInterpolantPivoted): """ ROM pivoted greedy rational interpolant computation for parametric problems (with pole matching). Args: HFEngine: HF problem solver. mu0(optional): Default parameter. Defaults to 0. directionPivot(optional): Pivot components. 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; - 'scaleFactorDer': scaling factors for derivative computation; defaults to 'AUTO'; - 'matchingWeight': weight for pole matching optimization; defaults to 1; - 'matchingMode': mode for pole matching optimization; allowed values include 'NONE' and 'SHIFT'; defaults to 'NONE'; - 'sharedRatio': required ratio of marginal points to share resonance; defaults to 1.; - 'matchingWeightError': weight for pole matching optimization in error estimation; defaults to 0; - 'cutOffToleranceError': tolerance for ignoring parasitic poles in error estimation; defaults to 'AUTO', i.e. cutOffTolerance; - 'S': total number of pivot samples current approximant relies upon; - 'samplerPivot': pivot sample point generator; - 'SMarginal': number of starting marginal samples; - 'samplerMarginal': marginal sample point generator via sparse grid; - 'errorEstimatorKindMarginal': kind of marginal error estimator; available values include 'LEAVE_ONE_OUT', 'LOOK_AHEAD', and 'LOOK_AHEAD_RECOVER'; defaults to 'LEAVE_ONE_OUT'; - 'polybasis': type of polynomial basis for pivot interpolation; defaults to 'MONOMIAL'; - 'polybasisMarginal': type of polynomial basis for marginal interpolation; allowed values include 'MONOMIAL_*', 'CHEBYSHEV_*', 'LEGENDRE_*', 'NEARESTNEIGHBOR', and 'PIECEWISE_LINEAR_*'; defaults to 'MONOMIAL'; - 'paramsMarginal': dictionary of parameters for marginal interpolation; include: . 'MMarginal': degree of marginal interpolant; defaults to 'AUTO', i.e. maximum allowed; not for 'NEARESTNEIGHBOR' or 'PIECEWISE_LINEAR_*'; . 'nNeighborsMarginal': number of marginal nearest neighbors; defaults to 1; only for 'NEARESTNEIGHBOR'; . 'polydegreetypeMarginal': type of polynomial degree for marginal; defaults to 'TOTAL'; not for 'NEARESTNEIGHBOR' or 'PIECEWISE_LINEAR_*'; . 'interpRcondMarginal': tolerance for marginal interpolation; defaults to None; not for 'NEARESTNEIGHBOR'; . 'radialDirectionalWeightsMarginalAdapt': bounds for adaptive rescaling of marginal radial basis weights; only for radial basis. - '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; - 'greedyTolMarginal': uniform error tolerance for marginal greedy algorithm; defaults to 1e-1; - 'maxIterMarginal': maximum number of marginal greedy steps; defaults to 1e2; - 'radialDirectionalWeights': radial basis weights for pivot numerator; defaults to 1; - 'radialDirectionalWeightsAdapt': bounds for adaptive rescaling of radial basis weights; defaults to [-1, -1]; - 'radialDirectionalWeightsMarginal': radial basis weights for marginal interpolant; defaults to 1; + - 'functionalSolve': strategy for minimization of denominator + functional; allowed values include 'NORM' and 'DOMINANT'; + defaults to 'NORM'; - 'interpRcond': tolerance for pivot interpolation; defaults to None; - 'robustTol': tolerance for robust rational denominator management; defaults to 0; - 'cutOffTolerance': tolerance for ignoring parasitic poles; defaults to np.inf; - 'residueTol': tolerance for residue elimination; defaults to 0., i.e. no bad residues. 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. directionPivot: Pivot components. mus: Array of snapshot parameters. musPivot: Array of pivot snapshot parameters. musMarginal: Array of marginal 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; - 'scaleFactorDer': scaling factors for derivative computation; - 'matchingWeight': weight for pole matching optimization; - 'matchingMode': mode for pole matching optimization; - 'sharedRatio': required ratio of marginal points to share resonance; - 'matchingWeightError': weight for pole matching optimization in error estimation; - 'cutOffToleranceError': tolerance for ignoring parasitic poles in error estimation; - 'errorEstimatorKindMarginal': kind of marginal error estimator; - 'polybasis': type of polynomial basis for pivot interpolation; - 'polybasisMarginal': type of polynomial basis for marginal interpolation; - 'paramsMarginal': dictionary of parameters for marginal interpolation; include: . 'MMarginal': degree of marginal interpolant; . 'nNeighborsMarginal': number of marginal nearest neighbors; . 'polydegreetypeMarginal': type of polynomial degree for marginal; . 'interpRcondMarginal': tolerance for marginal interpolation; . 'radialDirectionalWeightsMarginalAdapt': bounds for adaptive rescaling of marginal radial basis weights. - 'M': degree of rational interpolant numerator; - 'N': degree of rational interpolant denominator; - 'greedyTolMarginal': uniform error tolerance for marginal greedy algorithm; - 'maxIterMarginal': maximum number of marginal greedy steps; - 'radialDirectionalWeights': radial basis weights for pivot numerator; - 'radialDirectionalWeightsAdapt': bounds for adaptive rescaling of radial basis weights; - 'radialDirectionalWeightsMarginal': radial basis weights for marginal interpolant; + - 'functionalSolve': strategy for minimization of denominator + functional; - 'interpRcond': tolerance for pivot interpolation; - 'robustTol': tolerance for robust rational denominator management; - 'cutOffTolerance': tolerance for ignoring parasitic poles; - 'residueTol': tolerance for residue elimination. parameterListCritical: Recognized keys of critical approximant parameters: - 'S': total number of pivot samples current approximant relies upon; - 'samplerPivot': pivot sample point generator; - 'SMarginal': total number of marginal samples current approximant relies upon; - 'samplerMarginal': marginal sample point generator via sparse grid. approx_state: Whether to approximate state. verbosity: Verbosity level. POD: Whether to compute POD of snapshots. scaleFactorDer: Scaling factors for derivative computation. matchingWeight: Weight for pole matching optimization. matchingMode: Mode for pole matching optimization. sharedRatio: Required ratio of marginal points to share resonance. matchingWeightError: Weight for pole matching optimization in error estimation. cutOffToleranceError: Tolerance for ignoring parasitic poles in error estimation. S: Total number of pivot samples current approximant relies upon. samplerPivot: Pivot sample point generator. SMarginal: Total number of marginal samples current approximant relies upon. samplerMarginal: Marginal sample point generator via sparse grid. errorEstimatorKindMarginal: Kind of marginal error estimator. polybasis: Type of polynomial basis for pivot interpolation. polybasisMarginal: Type of polynomial basis for marginal interpolation. paramsMarginal: Dictionary of parameters for marginal interpolation. M: Degree of rational interpolant numerator. N: Degree of rational interpolant denominator. greedyTolMarginal: Uniform error tolerance for marginal greedy algorithm. maxIterMarginal: Maximum number of marginal greedy steps. radialDirectionalWeights: Radial basis weights for pivot numerator. radialDirectionalWeightsAdapt: Bounds for adaptive rescaling of radial basis weights. radialDirectionalWeightsMarginal: Radial basis weights for marginal interpolant. + functionalSolve: Strategy for minimization of denominator functional. interpRcond: Tolerance for pivot interpolation. robustTol: Tolerance for robust rational denominator management. cutOffTolerance: Tolerance for ignoring parasitic poles. residueTol: Tolerance for residue elimination. muBounds: list of bounds for pivot parameter values. muBoundsMarginal: list of bounds for marginal 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. """ diff --git a/rrompy/reduction_methods/pivoted/rational_interpolant_greedy_pivoted.py b/rrompy/reduction_methods/pivoted/rational_interpolant_greedy_pivoted.py index fb9f77e..1edf004 100644 --- a/rrompy/reduction_methods/pivoted/rational_interpolant_greedy_pivoted.py +++ b/rrompy/reduction_methods/pivoted/rational_interpolant_greedy_pivoted.py @@ -1,567 +1,579 @@ # 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 .generic_pivoted_approximant import (GenericPivotedApproximantBase, GenericPivotedApproximantNoMatch, GenericPivotedApproximant) from .gather_pivoted_approximant import gatherPivotedApproximant from rrompy.reduction_methods.standard.greedy.rational_interpolant_greedy \ import RationalInterpolantGreedy from rrompy.reduction_methods.standard.greedy.generic_greedy_approximant \ import pruneSamples from rrompy.utilities.base.types import Np1D from rrompy.utilities.base import verbosityManager as vbMng from rrompy.utilities.poly_fitting.polynomial import polyvander as pv from rrompy.utilities.exception_manager import RROMPyAssert, RROMPyWarning from rrompy.parameter import emptyParameterList, parameterList from rrompy.utilities.parallel import poolRank, indicesScatter, isend, recv __all__ = ['RationalInterpolantGreedyPivotedNoMatch', 'RationalInterpolantGreedyPivoted'] class RationalInterpolantGreedyPivotedBase(GenericPivotedApproximantBase, RationalInterpolantGreedy): @property def tModelType(self): if hasattr(self, "_temporaryPivot"): return RationalInterpolantGreedy.tModelType.fget(self) return super().tModelType @property def residueTol(self): """Value of residueTol.""" return self._residueTol @residueTol.setter def residueTol(self, residueTol): if residueTol < 0. or (residueTol > 0. and self.nparPivot > 1): RROMPyWarning("Overriding prescribed residue tolerance to 0.") residueTol = 0. self._residueTol = residueTol self._approxParameters["residueTol"] = self.residueTol def _polyvanderAuxiliary(self, mus, deg, *args): degEff = [0] * self.npar degEff[self.directionPivot[0]] = deg return pv(mus, degEff, *args) def _marginalizeMiscellanea(self, forward:bool): if forward: self._m_mu0 = copy(self.mu0) self._m_selfmus = copy(self.mus) self._m_HFEparameterMap = copy(self.HFEngine.parameterMap) self._mu0 = self.checkParameterListPivot( self.mu0(self.directionPivot)) self._mus = self.checkParameterListPivot( self.mus(self.directionPivot)) self.HFEngine.parameterMap = { "F": [self.HFEngine.parameterMap["F"][self.directionPivot[0]]], "B": [self.HFEngine.parameterMap["B"][self.directionPivot[0]]]} else: self._mu0 = self._m_mu0 self._mus = self._m_selfmus self.HFEngine.parameterMap = self._m_HFEparameterMap del self._m_mu0, self._m_selfmus, self._m_HFEparameterMap def _marginalizeTrainedModel(self, forward:bool): if forward: del self._temporaryPivot self.trainedModel.data.mu0 = self.mu0 self.trainedModel.data.scaleFactor = [1.] * self.npar self.trainedModel.data.scaleFactor[self.directionPivot[0]] = ( self.scaleFactor[0]) self.trainedModel.data.parameterMap = self.HFEngine.parameterMap Qc = np.zeros((1,) * self.directionPivot[0] + (len(self.trainedModel.data.Q.coeffs),) + (1,) * (self.npar - self.directionPivot[0] - 1), dtype = self.trainedModel.data.Q.coeffs.dtype) Pc = np.zeros((1,) * self.directionPivot[0] + (len(self.trainedModel.data.P.coeffs),) + (1,) * (self.npar - self.directionPivot[0] - 1) + (self.trainedModel.data.P.coeffs.shape[1],), dtype = self.trainedModel.data.P.coeffs.dtype) for j in range(len(self.trainedModel.data.Q.coeffs)): Qc[(0,) * self.directionPivot[0] + (j,) + (0,) * (self.npar - self.directionPivot[0] - 1)] = ( self.trainedModel.data.Q.coeffs[j]) for j in range(len(self.trainedModel.data.P.coeffs)): for k in range(self.trainedModel.data.P.coeffs.shape[1]): Pc[(0,) * self.directionPivot[0] + (j,) + (0,) * (self.npar - self.directionPivot[0] - 1) + (k,)] = self.trainedModel.data.P.coeffs[j, k] self.trainedModel.data.Q.coeffs = Qc self.trainedModel.data.P.coeffs = Pc self._m_musUniqueCN = copy(self._musUniqueCN) musUniqueCNAux = np.zeros((self.S, self.npar), dtype = self._musUniqueCN.dtype) musUniqueCNAux[:, self.directionPivot[0]] = self._musUniqueCN(0) self._musUniqueCN = self.checkParameterList(musUniqueCNAux) self._m_derIdxs = copy(self._derIdxs) for j in range(len(self._derIdxs)): for l in range(len(self._derIdxs[j])): derjl = self._derIdxs[j][l][0] self._derIdxs[j][l] = [0] * self.npar self._derIdxs[j][l][self.directionPivot[0]] = derjl else: self._temporaryPivot = 1 self.trainedModel.data.mu0 = self.checkParameterListPivot( self.mu0(self.directionPivot)) self.trainedModel.data.scaleFactor = self.scaleFactor self.trainedModel.data.parameterMap = { "F": [self.HFEngine.parameterMap["F"][self.directionPivot[0]]], "B": [self.HFEngine.parameterMap["B"][self.directionPivot[0]]]} self.trainedModel.data.Q.coeffs = self.trainedModel.data.Q.coeffs[ (0,) * self.directionPivot[0] + (slice(None),) + (0,) * (self.HFEngine.npar - 1 - self.directionPivot[0])] self.trainedModel.data.P.coeffs = self.trainedModel.data.P.coeffs[ (0,) * self.directionPivot[0] + (slice(None),) + (0,) * (self.HFEngine.npar - 1 - self.directionPivot[0])] self._musUniqueCN = copy(self._m_musUniqueCN) self._derIdxs = copy(self._m_derIdxs) del self._m_musUniqueCN, self._m_derIdxs self.trainedModel.data.npar = self.npar self.trainedModel.data.Q.npar = self.npar self.trainedModel.data.P.npar = self.npar def errorEstimator(self, mus:Np1D, return_max : bool = False) -> Np1D: """Standard residual-based error estimator.""" self._marginalizeMiscellanea(True) setupOK = self.setupApproxLocal() self._marginalizeMiscellanea(False) if setupOK > 0: err = np.empty(len(mus)) err[:] = np.nan if not return_max: return err return err, [- setupOK], np.nan self._marginalizeTrainedModel(True) errRes = super().errorEstimator(mus, return_max) self._marginalizeTrainedModel(False) return errRes def _preliminaryTraining(self): """Initialize starting snapshots of solution map.""" RROMPyAssert(self._mode, message = "Cannot start greedy algorithm.") self._S = self._setSampleBatch(self.S) self.resetSamples() self.samplingEngine.scaleFactor = self.scaleFactorDer musPivot = self.trainSetGenerator.generatePoints(self.S) while len(musPivot) > self.S: musPivot.pop() muTestPivot = self.samplerPivot.generatePoints(self.nTestPoints, False) idxPop = pruneSamples(self.HFEngine.mapParameterList(muTestPivot, idx = self.directionPivot), self.HFEngine.mapParameterList(musPivot, idx = self.directionPivot), 1e-10 * self.scaleFactorPivot[0]) self._mus = emptyParameterList() self.mus.reset((self.S, self.npar + len(self.musMargLoc))) muTestBase = emptyParameterList() muTestBase.reset((len(muTestPivot), self.npar + len(self.musMargLoc))) for k in range(self.S): muk = np.empty_like(self.mus[0]) muk[self.directionPivot] = musPivot[k] muk[self.directionMarginal] = self.musMargLoc self.mus[k] = muk for k in range(len(muTestPivot)): muk = np.empty_like(muTestBase[0]) muk[self.directionPivot] = muTestPivot[k] muk[self.directionMarginal] = self.musMargLoc muTestBase[k] = muk muTestBase.pop(idxPop) muLast = copy(self.mus[-1]) self.mus.pop() if len(self.mus) > 0: vbMng(self, "MAIN", ("Adding first {} sample point{} at {} to training " "set.").format(self.S - 1, "" + "s" * (self.S > 2), self.mus), 3) self.samplingEngine.iterSample(self.mus) self._S = len(self.mus) self._approxParameters["S"] = self.S self.muTest = parameterList(muTestBase) self.muTest.append(muLast) self.M, self.N = ("AUTO",) * 2 def setupApprox(self, *args, **kwargs) -> 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.computeScaleFactor() self._musMarginal = self.samplerMarginal.generatePoints(self.SMarginal) while len(self.musMarginal) > self.SMarginal: self.musMarginal.pop() S0 = copy(self.S) idx, sizes = indicesScatter(len(self.musMarginal), return_sizes = True) pMat, Ps, Qs, mus = None, [], [], None req, emptyCores = [], np.where(np.logical_not(sizes))[0] if len(idx) == 0: vbMng(self, "MAIN", "Idling.", 25) if self.storeAllSamples: self.storeSamples() pL, pT, mT = recv(source = 0, tag = poolRank()) pMat = np.empty((pL, 0), dtype = pT) mus = np.empty((0, self.mu0.shape[1]), dtype = mT) else: _scaleFactorOldPivot = copy(self.scaleFactor) self.scaleFactor = self.scaleFactorPivot self._temporaryPivot = 1 for i in idx: self.musMargLoc = self.musMarginal[i] vbMng(self, "MAIN", "Building marginal model no. {} at {}.".format(i + 1, self.musMargLoc), 5) self.samplingEngine.resetHistory() self.trainedModel = None self.verbosity -= 5 self.samplingEngine.verbosity -= 5 super().setupApprox(*args, **kwargs) self.verbosity += 5 self.samplingEngine.verbosity += 5 if self.storeAllSamples: self.storeSamples(i) if pMat is None: mus = copy(self.samplingEngine.mus.data) pMat = copy(self.samplingEngine.projectionMatrix) if i == 0: for dest in emptyCores: req += [isend((len(pMat), pMat.dtype, mus.dtype), dest = dest, tag = dest)] else: mus = np.vstack((mus, self.samplingEngine.mus.data)) pMat = np.hstack((pMat, self.samplingEngine.projectionMatrix)) Ps += [copy(self.trainedModel.data.P)] Qs += [copy(self.trainedModel.data.Q)] self._S = S0 del self._temporaryPivot, self.musMargLoc self.scaleFactor = _scaleFactorOldPivot for r in req: r.wait() pMat, Ps, Qs, mus, nsamples = gatherPivotedApproximant(pMat, Ps, Qs, mus, sizes, self.polybasis) self._mus = self.checkParameterList(mus) Psupp = np.append(0, np.cumsum(nsamples)) self._setupTrainedModel(pMat, forceNew = True) self.trainedModel.data.Qs, self.trainedModel.data.Ps = Qs, Ps self.trainedModel.data.Psupp = list(Psupp[: -1]) self._poleMatching() self._finalizeMarginalization() vbMng(self, "DEL", "Done setting up approximant.", 5) return 0 class RationalInterpolantGreedyPivotedNoMatch( RationalInterpolantGreedyPivotedBase, GenericPivotedApproximantNoMatch): """ ROM pivoted rational interpolant (without pole matching) computation for parametric problems. Args: HFEngine: HF problem solver. mu0(optional): Default parameter. Defaults to 0. directionPivot(optional): Pivot components. 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; - 'scaleFactorDer': scaling factors for derivative computation; defaults to 'AUTO'; - 'S': total number of pivot samples current approximant relies upon; - 'samplerPivot': pivot sample point generator; - 'SMarginal': total number of marginal samples current approximant relies upon; - 'samplerMarginal': marginal sample point generator; - 'polybasis': type of polynomial basis for pivot interpolation; defaults to 'MONOMIAL'; - 'greedyTol': uniform error tolerance for greedy algorithm; defaults to 1e-2; - 'collinearityTol': collinearity tolerance for greedy algorithm; defaults to 0.; - 'maxIter': maximum number of greedy steps; defaults to 1e2; - 'nTestPoints': number of test points; defaults to 5e2; - 'trainSetGenerator': training sample points generator; defaults to sampler; - 'errorEstimatorKind': kind of error estimator; available values include 'AFFINE', 'DISCREPANCY', 'LOOK_AHEAD', 'LOOK_AHEAD_RES', and 'NONE'; defaults to 'NONE'; - 'radialDirectionalWeightsMarginal': radial basis weights for marginal interpolant; defaults to 1; + - 'functionalSolve': strategy for minimization of denominator + functional; allowed values include 'NORM' and 'DOMINANT'; + defaults to 'NORM'; - 'interpRcond': tolerance for pivot interpolation; defaults to None; - 'robustTol': tolerance for robust rational denominator management; defaults to 0; - 'cutOffTolerance': tolerance for ignoring parasitic poles; defaults to np.inf; - 'residueTol': tolerance for residue elimination; defaults to 0., i.e. no bad residues. 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. directionPivot: Pivot components. mus: Array of snapshot parameters. musMarginal: Array of marginal 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; - 'scaleFactorDer': scaling factors for derivative computation; - 'polybasis': type of polynomial basis for pivot interpolation; - 'greedyTol': uniform error tolerance for greedy algorithm; - 'collinearityTol': collinearity tolerance for greedy algorithm; - 'maxIter': maximum number of greedy steps; - 'nTestPoints': number of test points; - 'trainSetGenerator': training sample points generator; - 'errorEstimatorKind': kind of error estimator; - 'radialDirectionalWeightsMarginal': radial basis weights for marginal interpolant; + - 'functionalSolve': strategy for minimization of denominator + functional; - 'interpRcond': tolerance for pivot interpolation; - 'robustTol': tolerance for robust rational denominator management; - 'cutOffTolerance': tolerance for ignoring parasitic poles; - 'residueTol': tolerance for residue elimination. parameterListCritical: Recognized keys of critical approximant parameters: - 'S': total number of pivot samples current approximant relies upon; - 'samplerPivot': pivot sample point generator; - 'SMarginal': total number of marginal samples current approximant relies upon; - 'samplerMarginal': marginal sample point generator. approx_state: Whether to approximate state. verbosity: Verbosity level. POD: Whether to compute POD of snapshots. scaleFactorDer: Scaling factors for derivative computation. S: Total number of pivot samples current approximant relies upon. samplerPivot: Pivot sample point generator. SMarginal: Total number of marginal samples current approximant relies upon. samplerMarginal: Marginal sample point generator. polybasis: Type of polynomial basis for pivot interpolation. greedyTol: uniform error tolerance for greedy algorithm. collinearityTol: Collinearity tolerance for greedy algorithm. maxIter: maximum number of greedy steps. nTestPoints: number of starting training points. trainSetGenerator: training sample points generator. errorEstimatorKind: kind of error estimator. radialDirectionalWeightsMarginal: Radial basis weights for marginal interpolant. + functionalSolve: Strategy for minimization of denominator functional. interpRcond: Tolerance for pivot interpolation. robustTol: Tolerance for robust rational denominator management. cutOffTolerance: Tolerance for ignoring parasitic poles. residueTol: Tolerance for residue elimination. muBounds: list of bounds for pivot parameter values. muBoundsMarginal: list of bounds for marginal 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. """ class RationalInterpolantGreedyPivoted(RationalInterpolantGreedyPivotedBase, GenericPivotedApproximant): """ ROM pivoted rational interpolant (with pole matching) computation for parametric problems. Args: HFEngine: HF problem solver. mu0(optional): Default parameter. Defaults to 0. directionPivot(optional): Pivot components. 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; - 'scaleFactorDer': scaling factors for derivative computation; defaults to 'AUTO'; - 'matchingWeight': weight for pole matching optimization; defaults to 1; - 'matchingMode': mode for pole matching optimization; allowed values include 'NONE' and 'SHIFT'; defaults to 'NONE'; - 'sharedRatio': required ratio of marginal points to share resonance; defaults to 1.; - 'S': total number of pivot samples current approximant relies upon; - 'samplerPivot': pivot sample point generator; - 'SMarginal': total number of marginal samples current approximant relies upon; - 'samplerMarginal': marginal sample point generator; - 'polybasis': type of polynomial basis for pivot interpolation; defaults to 'MONOMIAL'; - 'polybasisMarginal': type of polynomial basis for marginal interpolation; allowed values include 'MONOMIAL_*', 'CHEBYSHEV_*', 'LEGENDRE_*', 'NEARESTNEIGHBOR', and 'PIECEWISE_LINEAR_*'; defaults to 'MONOMIAL'; - 'paramsMarginal': dictionary of parameters for marginal interpolation; include: . 'MMarginal': degree of marginal interpolant; defaults to 'AUTO', i.e. maximum allowed; not for 'NEARESTNEIGHBOR' or 'PIECEWISE_LINEAR_*'; . 'nNeighborsMarginal': number of marginal nearest neighbors; defaults to 1; only for 'NEARESTNEIGHBOR'; . 'polydegreetypeMarginal': type of polynomial degree for marginal; defaults to 'TOTAL'; not for 'NEARESTNEIGHBOR' or 'PIECEWISE_LINEAR_*'; . 'interpRcondMarginal': tolerance for marginal interpolation; defaults to None; not for 'NEARESTNEIGHBOR'; . 'radialDirectionalWeightsMarginalAdapt': bounds for adaptive rescaling of marginal radial basis weights; only for radial basis. - 'greedyTol': uniform error tolerance for greedy algorithm; defaults to 1e-2; - 'collinearityTol': collinearity tolerance for greedy algorithm; defaults to 0.; - 'maxIter': maximum number of greedy steps; defaults to 1e2; - 'nTestPoints': number of test points; defaults to 5e2; - 'trainSetGenerator': training sample points generator; defaults to sampler; - 'errorEstimatorKind': kind of error estimator; available values include 'AFFINE', 'DISCREPANCY', 'LOOK_AHEAD', 'LOOK_AHEAD_RES', and 'NONE'; defaults to 'NONE'; - 'radialDirectionalWeightsMarginal': radial basis weights for marginal interpolant; defaults to 1; + - 'functionalSolve': strategy for minimization of denominator + functional; allowed values include 'NORM' and 'DOMINANT'; + defaults to 'NORM'; - 'interpRcond': tolerance for pivot interpolation; defaults to None; - 'robustTol': tolerance for robust rational denominator management; defaults to 0; - 'cutOffTolerance': tolerance for ignoring parasitic poles; defaults to np.inf; - 'residueTol': tolerance for residue elimination; defaults to 0., i.e. no bad residues. 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. directionPivot: Pivot components. mus: Array of snapshot parameters. musMarginal: Array of marginal 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; - 'scaleFactorDer': scaling factors for derivative computation; - 'matchingWeight': weight for pole matching optimization; - 'matchingMode': mode for pole matching optimization; - 'sharedRatio': required ratio of marginal points to share resonance; - 'polybasis': type of polynomial basis for pivot interpolation; - 'polybasisMarginal': type of polynomial basis for marginal interpolation; - 'paramsMarginal': dictionary of parameters for marginal interpolation; include: . 'MMarginal': degree of marginal interpolant; . 'nNeighborsMarginal': number of marginal nearest neighbors; . 'polydegreetypeMarginal': type of polynomial degree for marginal; . 'interpRcondMarginal': tolerance for marginal interpolation; . 'radialDirectionalWeightsMarginalAdapt': bounds for adaptive rescaling of marginal radial basis weights. - 'greedyTol': uniform error tolerance for greedy algorithm; - 'collinearityTol': collinearity tolerance for greedy algorithm; - 'maxIter': maximum number of greedy steps; - 'nTestPoints': number of test points; - 'trainSetGenerator': training sample points generator; - 'errorEstimatorKind': kind of error estimator; - 'radialDirectionalWeightsMarginal': radial basis weights for marginal interpolant; + - 'functionalSolve': strategy for minimization of denominator + functional; - 'interpRcond': tolerance for pivot interpolation; - 'robustTol': tolerance for robust rational denominator management; - 'cutOffTolerance': tolerance for ignoring parasitic poles; - 'residueTol': tolerance for residue elimination. parameterListCritical: Recognized keys of critical approximant parameters: - 'S': total number of pivot samples current approximant relies upon; - 'samplerPivot': pivot sample point generator; - 'SMarginal': total number of marginal samples current approximant relies upon; - 'samplerMarginal': marginal sample point generator. approx_state: Whether to approximate state. verbosity: Verbosity level. POD: Whether to compute POD of snapshots. scaleFactorDer: Scaling factors for derivative computation. matchingWeight: Weight for pole matching optimization. matchingMode: Mode for pole matching optimization. sharedRatio: Required ratio of marginal points to share resonance. S: Total number of pivot samples current approximant relies upon. samplerPivot: Pivot sample point generator. SMarginal: Total number of marginal samples current approximant relies upon. samplerMarginal: Marginal sample point generator. polybasis: Type of polynomial basis for pivot interpolation. polybasisMarginal: Type of polynomial basis for marginal interpolation. paramsMarginal: Dictionary of parameters for marginal interpolation. greedyTol: uniform error tolerance for greedy algorithm. collinearityTol: Collinearity tolerance for greedy algorithm. maxIter: maximum number of greedy steps. nTestPoints: number of starting training points. trainSetGenerator: training sample points generator. errorEstimatorKind: kind of error estimator. radialDirectionalWeightsMarginal: Radial basis weights for marginal interpolant. + functionalSolve: Strategy for minimization of denominator functional. interpRcond: Tolerance for pivot interpolation. robustTol: Tolerance for robust rational denominator management. cutOffTolerance: Tolerance for ignoring parasitic poles. residueTol: Tolerance for residue elimination. muBounds: list of bounds for pivot parameter values. muBoundsMarginal: list of bounds for marginal 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. """ diff --git a/rrompy/reduction_methods/pivoted/rational_interpolant_pivoted.py b/rrompy/reduction_methods/pivoted/rational_interpolant_pivoted.py index 71ba112..9af24b2 100644 --- a/rrompy/reduction_methods/pivoted/rational_interpolant_pivoted.py +++ b/rrompy/reduction_methods/pivoted/rational_interpolant_pivoted.py @@ -1,469 +1,481 @@ # 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 .generic_pivoted_approximant import (GenericPivotedApproximantBase, GenericPivotedApproximantNoMatch, GenericPivotedApproximant) from .gather_pivoted_approximant import gatherPivotedApproximant from rrompy.reduction_methods.standard.rational_interpolant import ( RationalInterpolant) from rrompy.utilities.base import verbosityManager as vbMng from rrompy.utilities.numerical.hash_derivative import nextDerivativeIndices from rrompy.utilities.exception_manager import RROMPyAssert, RROMPyWarning from rrompy.parameter import emptyParameterList from rrompy.utilities.parallel import poolRank, indicesScatter, isend, recv __all__ = ['RationalInterpolantPivotedNoMatch', 'RationalInterpolantPivoted'] class RationalInterpolantPivotedBase(GenericPivotedApproximantBase, RationalInterpolant): def __init__(self, *args, **kwargs): self._preInit() self._addParametersToList(toBeExcluded = ["polydegreetype"]) super().__init__(*args, **kwargs) self._postInit() @property def scaleFactorDer(self): """Value of scaleFactorDer.""" if self._scaleFactorDer == "NONE": return 1. if self._scaleFactorDer == "AUTO": return self.scaleFactorPivot return self._scaleFactorDer @scaleFactorDer.setter def scaleFactorDer(self, scaleFactorDer): if isinstance(scaleFactorDer, (str,)): scaleFactorDer = scaleFactorDer.upper() elif hasattr(scaleFactorDer, "__len__"): scaleFactorDer = list(scaleFactorDer) self._scaleFactorDer = scaleFactorDer self._approxParameters["scaleFactorDer"] = self._scaleFactorDer @property def polydegreetype(self): """Value of polydegreetype.""" return "TOTAL" @polydegreetype.setter def polydegreetype(self, polydegreetype): RROMPyWarning(("polydegreetype is used just to simplify inheritance, " "and its value cannot be changed from 'TOTAL'.")) @property def residueTol(self): """Value of residueTol.""" return self._residueTol @residueTol.setter def residueTol(self, residueTol): if residueTol < 0. or (residueTol > 0. and self.nparPivot > 1): RROMPyWarning("Overriding prescribed residue tolerance to 0.") residueTol = 0. self._residueTol = residueTol self._approxParameters["residueTol"] = self.residueTol def _setupInterpolationIndices(self): """Setup parameters for polyvander.""" RROMPyAssert(self._mode, message = "Cannot setup interpolation indices.") if (self._musUniqueCN is None or len(self._reorder) != len(self.musPivot)): try: muPC = self.trainedModel.centerNormalizePivot(self.musPivot) except: muPC = self.trainedModel.centerNormalize(self.musPivot) self._musUniqueCN, musIdxsTo, musIdxs, musCount = (muPC.unique( return_index = True, return_inverse = True, return_counts = True)) self._musUnique = self.musPivot[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.nparPivot, cnt) jIdx = np.nonzero(musIdxs == j)[0] self._reorder[jIdx] = np.arange(filled, filled + cnt) filled += cnt 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.computeScaleFactor() self.resetSamples() self.samplingEngine.scaleFactor = self.scaleFactorDer self.musPivot = self.samplerPivot.generatePoints(self.S) while len(self.musPivot) > self.S: self.musPivot.pop() self._musMarginal = self.samplerMarginal.generatePoints(self.SMarginal) while len(self.musMarginal) > self.SMarginal: self.musMarginal.pop() self._mus = emptyParameterList() self.mus.reset((self.S * self.SMarginal, self.HFEngine.npar)) for j, muMarg in enumerate(self.musMarginal): for k in range(j * self.S, (j + 1) * self.S): muk = np.empty_like(self.mus[0]) muk[self.directionPivot] = self.musPivot[k - j * self.S] muk[self.directionMarginal] = muMarg self.mus[k] = muk N0 = copy(self.N) self._setupTrainedModel(np.zeros((0, 0)), forceNew = True) idx, sizes = indicesScatter(len(self.musMarginal), return_sizes = True) pMat, Ps, Qs = None, [], [] req, emptyCores = [], np.where(np.logical_not(sizes))[0] if len(idx) == 0: vbMng(self, "MAIN", "Idling.", 30) if self.storeAllSamples: self.storeSamples() pL, pT = recv(source = 0, tag = poolRank()) pMat = np.empty((pL, 0), dtype = pT) else: _scaleFactorOldPivot = copy(self.scaleFactor) self.scaleFactor = self.scaleFactorPivot self._temporaryPivot = 1 for i in idx: vbMng(self, "MAIN", "Building marginal model no. {} at {}.".format(i + 1, self.musMarginal[i]), 5) vbMng(self, "INIT", "Starting computation of snapshots.", 10) self.samplingEngine.resetHistory() self.samplingEngine.iterSample( self.mus[self.S * i : self.S * (i + 1)]) vbMng(self, "DEL", "Done computing snapshots.", 10) self.verbosity -= 5 self.samplingEngine.verbosity -= 5 - self._setupRational(self._setupDenominator()[0]) + self._setupRational(self._setupDenominator()) self.verbosity += 5 self.samplingEngine.verbosity += 5 if self.storeAllSamples: self.storeSamples(i) if pMat is None: pMat = copy(self.samplingEngine.projectionMatrix) if i == 0: for dest in emptyCores: req += [isend((len(pMat), pMat.dtype), dest = dest, tag = dest)] else: pMat = np.hstack((pMat, self.samplingEngine.projectionMatrix)) Ps += [copy(self.trainedModel.data.P)] Qs += [copy(self.trainedModel.data.Q)] del self.trainedModel.data.Q, self.trainedModel.data.P self.N = N0 del self._temporaryPivot self.scaleFactor = _scaleFactorOldPivot for r in req: r.wait() pMat, Ps, Qs, _, _ = gatherPivotedApproximant(pMat, Ps, Qs, self.mus.data, sizes, self.polybasis, False) self._setupTrainedModel(pMat) self.trainedModel.data.Qs, self.trainedModel.data.Ps = Qs, Ps Psupp = np.arange(0, len(self.musMarginal) * self.S, self.S) self.trainedModel.data.Psupp = list(Psupp) self._poleMatching() self._finalizeMarginalization() vbMng(self, "DEL", "Done setting up approximant.", 5) return 0 class RationalInterpolantPivotedNoMatch(RationalInterpolantPivotedBase, GenericPivotedApproximantNoMatch): """ ROM pivoted rational interpolant (without pole matching) computation for parametric problems. Args: HFEngine: HF problem solver. mu0(optional): Default parameter. Defaults to 0. directionPivot(optional): Pivot components. 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; - 'scaleFactorDer': scaling factors for derivative computation; defaults to 'AUTO'; - 'S': total number of pivot samples current approximant relies upon; - 'samplerPivot': pivot sample point generator; - 'SMarginal': total number of marginal samples current approximant relies upon; - 'samplerMarginal': marginal sample point generator; - 'polybasis': type of polynomial basis for pivot interpolation; defaults to 'MONOMIAL'; - '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; - 'radialDirectionalWeights': radial basis weights for pivot numerator; defaults to 1; - 'radialDirectionalWeightsAdapt': bounds for adaptive rescaling of radial basis weights; defaults to [-1, -1]; - 'radialDirectionalWeightsMarginal': radial basis weights for marginal interpolant; defaults to 1; + - 'functionalSolve': strategy for minimization of denominator + functional; allowed values include 'NORM' and 'DOMINANT'; + defaults to 'NORM'; - 'interpRcond': tolerance for pivot interpolation; defaults to None; - 'robustTol': tolerance for robust rational denominator management; defaults to 0; - 'cutOffTolerance': tolerance for ignoring parasitic poles; defaults to np.inf; - 'residueTol': tolerance for residue elimination; defaults to 0., i.e. no bad residues. 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. directionPivot: Pivot components. mus: Array of snapshot parameters. musPivot: Array of pivot snapshot parameters. musMarginal: Array of marginal 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; - 'scaleFactorDer': scaling factors for derivative computation; - 'polybasis': type of polynomial basis for pivot interpolation; - 'M': degree of rational interpolant numerator; - 'N': degree of rational interpolant denominator; - 'radialDirectionalWeights': radial basis weights for pivot numerator; - 'radialDirectionalWeightsAdapt': bounds for adaptive rescaling of radial basis weights; - 'radialDirectionalWeightsMarginal': radial basis weights for marginal interpolant; + - 'functionalSolve': strategy for minimization of denominator + functional; - 'interpRcond': tolerance for pivot interpolation; - 'robustTol': tolerance for robust rational denominator management; - 'cutOffTolerance': tolerance for ignoring parasitic poles; - 'residueTol': tolerance for residue elimination. parameterListCritical: Recognized keys of critical approximant parameters: - 'S': total number of pivot samples current approximant relies upon; - 'samplerPivot': pivot sample point generator; - 'SMarginal': total number of marginal samples current approximant relies upon; - 'samplerMarginal': marginal sample point generator. approx_state: Whether to approximate state. verbosity: Verbosity level. POD: Whether to compute POD of snapshots. scaleFactorDer: Scaling factors for derivative computation. S: Total number of pivot samples current approximant relies upon. samplerPivot: Pivot sample point generator. SMarginal: Total number of marginal samples current approximant relies upon. samplerMarginal: Marginal sample point generator. polybasis: Type of polynomial basis for pivot interpolation. M: Numerator degree of approximant. N: Denominator degree of approximant. radialDirectionalWeights: Radial basis weights for pivot numerator. radialDirectionalWeightsAdapt: Bounds for adaptive rescaling of radial basis weights. radialDirectionalWeightsMarginal: Radial basis weights for marginal interpolant. + functionalSolve: Strategy for minimization of denominator functional. interpRcond: Tolerance for pivot interpolation. robustTol: Tolerance for robust rational denominator management. cutOffTolerance: Tolerance for ignoring parasitic poles. residueTol: Tolerance for residue elimination. muBounds: list of bounds for pivot parameter values. muBoundsMarginal: list of bounds for marginal 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. """ class RationalInterpolantPivoted(RationalInterpolantPivotedBase, GenericPivotedApproximant): """ ROM pivoted rational interpolant (with pole matching) computation for parametric problems. Args: HFEngine: HF problem solver. mu0(optional): Default parameter. Defaults to 0. directionPivot(optional): Pivot components. 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; - 'scaleFactorDer': scaling factors for derivative computation; defaults to 'AUTO'; - 'matchingWeight': weight for pole matching optimization; defaults to 1; - 'matchingMode': mode for pole matching optimization; allowed values include 'NONE' and 'SHIFT'; defaults to 'NONE'; - 'sharedRatio': required ratio of marginal points to share resonance; defaults to 1.; - 'S': total number of pivot samples current approximant relies upon; - 'samplerPivot': pivot sample point generator; - 'SMarginal': total number of marginal samples current approximant relies upon; - 'samplerMarginal': marginal sample point generator; - 'polybasis': type of polynomial basis for pivot interpolation; defaults to 'MONOMIAL'; - 'polybasisMarginal': type of polynomial basis for marginal interpolation; allowed values include 'MONOMIAL_*', 'CHEBYSHEV_*', 'LEGENDRE_*', 'NEARESTNEIGHBOR', and 'PIECEWISE_LINEAR_*'; defaults to 'MONOMIAL'; - 'paramsMarginal': dictionary of parameters for marginal interpolation; include: . 'MMarginal': degree of marginal interpolant; defaults to 'AUTO', i.e. maximum allowed; not for 'NEARESTNEIGHBOR' or 'PIECEWISE_LINEAR_*'; . 'nNeighborsMarginal': number of marginal nearest neighbors; defaults to 1; only for 'NEARESTNEIGHBOR'; . 'polydegreetypeMarginal': type of polynomial degree for marginal; defaults to 'TOTAL'; not for 'NEARESTNEIGHBOR' or 'PIECEWISE_LINEAR_*'; . 'interpRcondMarginal': tolerance for marginal interpolation; defaults to None; not for 'NEARESTNEIGHBOR'; . 'radialDirectionalWeightsMarginalAdapt': bounds for adaptive rescaling of marginal radial basis weights; only for radial basis. - '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; - 'radialDirectionalWeights': radial basis weights for pivot numerator; defaults to 1; - 'radialDirectionalWeightsAdapt': bounds for adaptive rescaling of radial basis weights; defaults to [-1, -1]; - 'radialDirectionalWeightsMarginal': radial basis weights for marginal interpolant; defaults to 1; + - 'functionalSolve': strategy for minimization of denominator + functional; allowed values include 'NORM' and 'DOMINANT'; + defaults to 'NORM'; - 'interpRcond': tolerance for pivot interpolation; defaults to None; - 'robustTol': tolerance for robust rational denominator management; defaults to 0; - 'cutOffTolerance': tolerance for ignoring parasitic poles; defaults to np.inf; - 'residueTol': tolerance for residue elimination; defaults to 0., i.e. no bad residues. 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. directionPivot: Pivot components. mus: Array of snapshot parameters. musPivot: Array of pivot snapshot parameters. musMarginal: Array of marginal 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; - 'scaleFactorDer': scaling factors for derivative computation; - 'matchingWeight': weight for pole matching optimization; - 'matchingMode': mode for pole matching optimization; - 'sharedRatio': required ratio of marginal points to share resonance; - 'polybasis': type of polynomial basis for pivot interpolation; - 'polybasisMarginal': type of polynomial basis for marginal interpolation; - 'paramsMarginal': dictionary of parameters for marginal interpolation; include: . 'MMarginal': degree of marginal interpolant; . 'nNeighborsMarginal': number of marginal nearest neighbors; . 'polydegreetypeMarginal': type of polynomial degree for marginal; . 'interpRcondMarginal': tolerance for marginal interpolation; . 'radialDirectionalWeightsMarginalAdapt': bounds for adaptive rescaling of marginal radial basis weights. - 'M': degree of rational interpolant numerator; - 'N': degree of rational interpolant denominator; - 'radialDirectionalWeights': radial basis weights for pivot numerator; - 'radialDirectionalWeightsAdapt': bounds for adaptive rescaling of radial basis weights; - 'radialDirectionalWeightsMarginal': radial basis weights for marginal interpolant; + - 'functionalSolve': strategy for minimization of denominator + functional; - 'interpRcond': tolerance for pivot interpolation; - 'robustTol': tolerance for robust rational denominator management; - 'cutOffTolerance': tolerance for ignoring parasitic poles; - 'residueTol': tolerance for residue elimination. parameterListCritical: Recognized keys of critical approximant parameters: - 'S': total number of pivot samples current approximant relies upon; - 'samplerPivot': pivot sample point generator; - 'SMarginal': total number of marginal samples current approximant relies upon; - 'samplerMarginal': marginal sample point generator. approx_state: Whether to approximate state. verbosity: Verbosity level. POD: Whether to compute POD of snapshots. scaleFactorDer: Scaling factors for derivative computation. matchingWeight: Weight for pole matching optimization. matchingMode: Mode for pole matching optimization. sharedRatio: Required ratio of marginal points to share resonance. S: Total number of pivot samples current approximant relies upon. samplerPivot: Pivot sample point generator. SMarginal: Total number of marginal samples current approximant relies upon. samplerMarginal: Marginal sample point generator. polybasis: Type of polynomial basis for pivot interpolation. polybasisMarginal: Type of polynomial basis for marginal interpolation. paramsMarginal: Dictionary of parameters for marginal interpolation. M: Numerator degree of approximant. N: Denominator degree of approximant. radialDirectionalWeights: Radial basis weights for pivot numerator. radialDirectionalWeightsAdapt: Bounds for adaptive rescaling of radial basis weights. radialDirectionalWeightsMarginal: Radial basis weights for marginal interpolant. + functionalSolve: Strategy for minimization of denominator functional. interpRcond: Tolerance for pivot interpolation. robustTol: Tolerance for robust rational denominator management. cutOffTolerance: Tolerance for ignoring parasitic poles. residueTol: Tolerance for residue elimination. muBounds: list of bounds for pivot parameter values. muBoundsMarginal: list of bounds for marginal 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. """ diff --git a/rrompy/reduction_methods/standard/greedy/rational_interpolant_greedy.py b/rrompy/reduction_methods/standard/greedy/rational_interpolant_greedy.py index 232335e..1171232 100644 --- a/rrompy/reduction_methods/standard/greedy/rational_interpolant_greedy.py +++ b/rrompy/reduction_methods/standard/greedy/rational_interpolant_greedy.py @@ -1,538 +1,544 @@ # 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.hfengines.base.linear_affine_engine import checkIfAffine from .generic_greedy_approximant import GenericGreedyApproximant from rrompy.utilities.poly_fitting.polynomial import (polybases, polyfitname, PolynomialInterpolator as PI, polyvander) from rrompy.utilities.numerical import dot from rrompy.utilities.numerical.degree import totalDegreeN from rrompy.utilities.expression import expressionEvaluator from rrompy.reduction_methods.standard import RationalInterpolant from rrompy.utilities.base.types import Np1D, Tuple, paramVal, List from rrompy.utilities.base.verbosity_depth import (verbosityManager as vbMng, getVerbosityDepth, setVerbosityDepth) from rrompy.utilities.poly_fitting import customFit from rrompy.utilities.exception_manager import (RROMPyWarning, RROMPyException, RROMPyAssert, RROMPy_FRAGILE) from rrompy.sampling import sampleList, emptySampleList __all__ = ['RationalInterpolantGreedy'] class RationalInterpolantGreedy(GenericGreedyApproximant, RationalInterpolant): """ ROM greedy 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; - 'scaleFactorDer': scaling factors for derivative computation; defaults to 'AUTO'; - 'S': number of starting training points; - 'sampler': sample point generator; - 'greedyTol': uniform error tolerance for greedy algorithm; defaults to 1e-2; - 'collinearityTol': collinearity tolerance for greedy algorithm; defaults to 0.; - 'maxIter': maximum number of greedy steps; defaults to 1e2; - 'nTestPoints': number of test points; defaults to 5e2; - 'trainSetGenerator': training sample points generator; defaults to sampler; - 'polybasis': type of basis for interpolation; defaults to 'MONOMIAL'; - 'errorEstimatorKind': kind of error estimator; available values include 'AFFINE', 'DISCREPANCY', 'LOOK_AHEAD', 'LOOK_AHEAD_RES', 'LOOK_AHEAD_OUTPUT', and 'NONE'; defaults to 'NONE'; + - 'functionalSolve': strategy for minimization of denominator + functional; allowed values include 'NORM' and 'DOMINANT'; + defaults to 'NORM'; - 'interpRcond': tolerance for interpolation; defaults to None; - 'robustTol': tolerance for robust rational denominator management; defaults to 0; - 'cutOffTolerance': tolerance for ignoring parasitic poles; defaults to np.inf; - 'residueTol': tolerance for residue elimination; defaults to 0., i.e. no bad residues. Defaults to empty dict. approx_state(optional): Whether to approximate state. Defaults and must be True. 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; - 'scaleFactorDer': scaling factors for derivative computation; - 'greedyTol': uniform error tolerance for greedy algorithm; - 'collinearityTol': collinearity tolerance for greedy algorithm; - 'maxIter': maximum number of greedy steps; - 'nTestPoints': number of test points; - 'trainSetGenerator': training sample points generator; - 'errorEstimatorKind': kind of error estimator; + - 'functionalSolve': strategy for minimization of denominator + functional; - 'interpRcond': tolerance for interpolation; - 'robustTol': tolerance for robust rational denominator management; - 'cutOffTolerance': tolerance for ignoring parasitic poles; - 'residueTol': tolerance for residue elimination. 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. scaleFactorDer: Scaling factors for derivative computation. S: number of test points. sampler: Sample point generator. greedyTol: uniform error tolerance for greedy algorithm. collinearityTol: Collinearity tolerance for greedy algorithm. maxIter: maximum number of greedy steps. nTestPoints: number of starting training points. trainSetGenerator: training sample points generator. robustTol: tolerance for robust rational denominator management. cutOffTolerance: Tolerance for ignoring parasitic poles. residueTol: Tolerance for residue elimination. errorEstimatorKind: kind of error estimator. + functionalSolve: Strategy for minimization of denominator functional. interpRcond: tolerance for interpolation. robustTol: tolerance for robust rational denominator management. muBounds: list of bounds for parameter values. samplingEngine: Sampling engine. estimatorNormEngine: Engine for estimator norm computation. 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. """ _allowedEstimatorKinds = ["AFFINE", "DISCREPANCY", "LOOK_AHEAD", "LOOK_AHEAD_RES", "LOOK_AHEAD_OUTPUT", "NONE"] def __init__(self, *args, **kwargs): self._preInit() self._addParametersToList(["errorEstimatorKind"], ["DISCREPANCY"], toBeExcluded = ["M", "N", "polydegreetype", "radialDirectionalWeights"]) super().__init__(*args, **kwargs) if not self.approx_state and self.errorEstimatorKind not in [ "LOOK_AHEAD", "LOOK_AHEAD_OUTPUT", "NONE"]: raise RROMPyException(("Must compute greedy approximation of " "state, unless error estimator allows " "otherwise.")) self._postInit() @property def approx_state(self): """Value of approx_state.""" return self._approx_state @approx_state.setter def approx_state(self, approx_state): RationalInterpolant.approx_state.fset(self, approx_state) if (not self.approx_state and hasattr(self, "_errorEstimatorKind") and self.errorEstimatorKind not in [ "LOOK_AHEAD", "LOOK_AHEAD_OUTPUT", "NONE"]): raise RROMPyException(("Must compute greedy approximation of " "state, unless error estimator allows " "otherwise.")) @property def E(self): """Value of E.""" self._E = self.sampleBatchIdx - 1 return self._E @E.setter def E(self, E): RROMPyWarning(("E is used just to simplify inheritance, and its value " "cannot be changed from that of sampleBatchIdx - 1.")) def _setMAuto(self): self.M = self.E def _setNAuto(self): self.N = self.E @property def polydegreetype(self): """Value of polydegreetype.""" return "TOTAL" @polydegreetype.setter def polydegreetype(self, polydegreetype): RROMPyWarning(("polydegreetype is used just to simplify inheritance, " "and its value cannot be changed from 'TOTAL'.")) @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("Sample type not recognized.") self._polybasis = polybasis except: RROMPyWarning(("Prescribed polybasis not recognized. Overriding " "to 'MONOMIAL'.")) self._polybasis = "MONOMIAL" self._approxParameters["polybasis"] = self.polybasis @property def errorEstimatorKind(self): """Value of errorEstimatorKind.""" return self._errorEstimatorKind @errorEstimatorKind.setter def errorEstimatorKind(self, errorEstimatorKind): errorEstimatorKind = errorEstimatorKind.upper() if errorEstimatorKind not in self._allowedEstimatorKinds: RROMPyWarning(("Error estimator kind not recognized. Overriding " "to 'NONE'.")) errorEstimatorKind = "NONE" self._errorEstimatorKind = errorEstimatorKind self._approxParameters["errorEstimatorKind"] = self.errorEstimatorKind if (self.errorEstimatorKind not in [ "LOOK_AHEAD", "LOOK_AHEAD_OUTPUT", "NONE"] and hasattr(self, "_approx_state") and not self.approx_state): raise RROMPyException(("Must compute greedy approximation of " "state, unless error estimator allows " "otherwise.")) def _polyvanderAuxiliary(self, mus, deg, *args): return polyvander(mus, deg, *args) def getErrorEstimatorDiscrepancy(self, mus:Np1D) -> Np1D: """Discrepancy-based residual estimator.""" checkIfAffine(self.HFEngine, "apply discrepancy-based error estimator") mus = self.checkParameterList(mus) muCTest = self.trainedModel.centerNormalize(mus) tMverb, self.trainedModel.verbosity = self.trainedModel.verbosity, 0 QTest = self.trainedModel.getQVal(mus) QTzero = np.where(QTest == 0.)[0] if len(QTzero) > 0: RROMPyWarning(("Adjusting estimator to avoid division by " "numerically zero denominator.")) QTest[QTzero] = np.finfo(np.complex).eps / (1. + self.N) self.HFEngine.buildA() self.HFEngine.buildb() nAs, nbs = self.HFEngine.nAs, self.HFEngine.nbs muTrainEff = self.HFEngine.mapParameterList(self.mus) muTestEff = self.HFEngine.mapParameterList(mus) PTrain = self.trainedModel.getPVal(self.mus).data.T QTrain = self.trainedModel.getQVal(self.mus) QTzero = np.where(QTrain == 0.)[0] if len(QTzero) > 0: RROMPyWarning(("Adjusting estimator to avoid division by " "numerically zero denominator.")) QTrain[QTzero] = np.finfo(np.complex).eps / (1. + self.N) PTest = self.trainedModel.getPVal(mus).data self.trainedModel.verbosity = tMverb radiusAbTrain = np.empty((self.S, nAs * self.S + nbs), dtype = np.complex) radiusA = np.empty((self.S, nAs, len(mus)), dtype = np.complex) radiusb = np.empty((nbs, len(mus)), dtype = np.complex) for j, thA in enumerate(self.HFEngine.thAs): idxs = j * self.S + np.arange(self.S) radiusAbTrain[:, idxs] = expressionEvaluator(thA[0], muTrainEff, (self.S, 1)) * PTrain radiusA[:, j] = PTest * expressionEvaluator(thA[0], muTestEff, (len(mus),)) for j, thb in enumerate(self.HFEngine.thbs): idx = nAs * self.S + j radiusAbTrain[:, idx] = QTrain * expressionEvaluator(thb[0], muTrainEff, (self.S,)) radiusb[j] = QTest * expressionEvaluator(thb[0], muTestEff, (len(mus),)) QRHSNorm2 = self._affineResidualMatricesContraction(radiusb) vanTrain = self._polyvanderAuxiliary(self._musUniqueCN, self.E, self.polybasis0, self._derIdxs, self._reorder) interpPQ = customFit(vanTrain, radiusAbTrain, rcond = self.interpRcond) vanTest = self._polyvanderAuxiliary(muCTest, self.E, self.polybasis0) DradiusAb = vanTest.dot(interpPQ) radiusA = (radiusA - DradiusAb[:, : - nbs].reshape(len(mus), -1, self.S).T) radiusb = radiusb - DradiusAb[:, - nbs :].T ff, Lf, LL = self._affineResidualMatricesContraction(radiusb, radiusA) err = np.abs((LL - 2. * np.real(Lf) + ff) / QRHSNorm2) ** .5 return err def getErrorEstimatorLookAhead(self, mus:Np1D, what : str = "") -> Tuple[Np1D, List[int]]: """Residual estimator based on look-ahead idea.""" errTest, QTest, idxMaxEst = self._EIMStep(mus) _approx_state_old = self.approx_state if what == "OUTPUT" and _approx_state_old: self._approx_state = False self.initEstimatorNormEngine() self._approx_state = _approx_state_old mu_muTestSample = mus[idxMaxEst] app_muTestSample = self.getApproxReduced(mu_muTestSample) if self._mode == RROMPy_FRAGILE: if what == "RES" and not self.HFEngine.isCEye: raise RROMPyException(("Cannot compute LOOK_AHEAD_RES " "estimator in fragile mode for " "non-scalar C.")) app_muTestSample = dot(self.trainedModel.data.projMat[:, : app_muTestSample.shape[0]], app_muTestSample) else: app_muTestSample = dot(self.samplingEngine.projectionMatrix, app_muTestSample) if what == "RES": errmu = self.HFEngine.residual(mu_muTestSample, app_muTestSample, post_c = False) solmu = self.HFEngine.residual(mu_muTestSample, None, post_c = False) else: for j, mu in enumerate(mu_muTestSample): uEx = self.samplingEngine.nextSample(mu) if j == 0: solmu = emptySampleList() solmu.reset((len(uEx), len(mu_muTestSample)), dtype = uEx.dtype) solmu[j] = uEx if what == "OUTPUT" and self.approx_state: solmu = sampleList(self.HFEngine.applyC(solmu)) app_muTestSample = sampleList(self.HFEngine.applyC( app_muTestSample)) errmu = solmu - app_muTestSample errsamples = (self.estimatorNormEngine.norm(errmu) / self.estimatorNormEngine.norm(solmu)) musT = copy(self.mus) musT.append(mu_muTestSample) musT = self.trainedModel.centerNormalize(musT) musC = self.trainedModel.centerNormalize(mus) errT = np.zeros((len(musT), len(mu_muTestSample)), dtype = np.complex) errT[np.arange(len(self.mus), len(musT)), np.arange(len(mu_muTestSample))] = errsamples * QTest[idxMaxEst] vanT = self._polyvanderAuxiliary(musT, self.E + 1, self.polybasis) fitOut = customFit(vanT, errT, full = True, rcond = self.interpRcond) vbMng(self, "MAIN", ("Fitting {} samples with degree {} through {}... Conditioning " "of LS system: {:.4e}.").format(len(vanT), self.E + 1, polyfitname(self.polybasis), fitOut[1][2][0] / fitOut[1][2][-1]), 15) vanC = self._polyvanderAuxiliary(musC, self.E + 1, self.polybasis) err = np.sum(np.abs(vanC.dot(fitOut[0])), axis = -1) / QTest return err, idxMaxEst def getErrorEstimatorNone(self, mus:Np1D) -> Np1D: """EIM-based residual estimator.""" err = np.max(self._EIMStep(mus, True), axis = 1) err *= self.greedyTol / np.mean(err) return err def _EIMStep(self, mus:Np1D, only_one : bool = False) -> Tuple[Np1D, Np1D, List[int]]: """Residual estimator based on look-ahead idea.""" mus = self.checkParameterList(mus) tMverb, self.trainedModel.verbosity = self.trainedModel.verbosity, 0 QTest = self.trainedModel.getQVal(mus) QTzero = np.where(QTest == 0.)[0] if len(QTzero) > 0: RROMPyWarning(("Adjusting estimator to avoid division by " "numerically zero denominator.")) QTest[QTzero] = np.finfo(np.complex).eps / (1. + self.N) QTest = np.abs(QTest) muCTest = self.trainedModel.centerNormalize(mus) muCTrain = self.trainedModel.centerNormalize(self.mus) self.trainedModel.verbosity = tMverb vanTest = self._polyvanderAuxiliary(muCTest, self.E, self.polybasis) vanTestNext = self._polyvanderAuxiliary(muCTest, self.E + 1, self.polybasis)[:, vanTest.shape[1] :] idxsTest = np.arange(vanTestNext.shape[1]) basis = np.zeros((len(idxsTest), 0), dtype = float) idxMaxEst = [] while len(idxsTest) > 0: vanTrial = self._polyvanderAuxiliary(muCTrain, self.E, self.polybasis) vanTrialNext = self._polyvanderAuxiliary(muCTrain, self.E + 1, self.polybasis)[:, vanTrial.shape[1] :] vanTrial = np.hstack((vanTrial, vanTrialNext.dot(basis).reshape( len(vanTrialNext), basis.shape[1]))) valuesTrial = vanTrialNext[:, idxsTest] vanTestEff = np.hstack((vanTest, vanTestNext.dot(basis).reshape( len(vanTestNext), basis.shape[1]))) vanTestNextEff = vanTestNext[:, idxsTest] try: coeffTest = np.linalg.solve(vanTrial, valuesTrial) except np.linalg.LinAlgError as e: raise RROMPyException(e) errTest = (np.abs(vanTestNextEff - vanTestEff.dot(coeffTest)) / np.expand_dims(QTest, 1)) if only_one: return errTest idxMaxErr = np.unravel_index(np.argmax(errTest), errTest.shape) idxMaxEst += [idxMaxErr[0]] muCTrain.append(muCTest[idxMaxErr[0]]) basis = np.pad(basis, [(0, 0), (0, 1)], "constant") basis[idxsTest[idxMaxErr[1]], -1] = 1. idxsTest = np.delete(idxsTest, idxMaxErr[1]) return errTest, QTest, idxMaxEst def errorEstimator(self, mus:Np1D, return_max : bool = False) -> Np1D: """Standard residual-based error estimator.""" setupOK = self.setupApproxLocal() if setupOK > 0: err = np.empty(len(mus)) err[:] = np.nan if not return_max: return err return err, [- setupOK], np.nan mus = self.checkParameterList(mus) vbMng(self.trainedModel, "INIT", "Evaluating error estimator at mu = {}.".format(mus), 10) if self.errorEstimatorKind == "AFFINE": err = self.getErrorEstimatorAffine(mus) else: self._setupInterpolationIndices() if self.errorEstimatorKind == "DISCREPANCY": err = self.getErrorEstimatorDiscrepancy(mus) elif self.errorEstimatorKind[: 10] == "LOOK_AHEAD": err, idxMaxEst = self.getErrorEstimatorLookAhead(mus, self.errorEstimatorKind[11 :]) else: #if self.errorEstimatorKind == "NONE": err = self.getErrorEstimatorNone(mus) vbMng(self.trainedModel, "DEL", "Done evaluating error estimator", 10) if not return_max: return err if self.errorEstimatorKind[: 10] != "LOOK_AHEAD": idxMaxEst = np.empty(self.sampleBatchSize, dtype = int) errCP = copy(err) for j in range(self.sampleBatchSize): k = np.argmax(errCP) idxMaxEst[j] = k if j + 1 < self.sampleBatchSize: musZero = self.trainedModel.centerNormalize(mus, mus[k]) errCP *= np.linalg.norm(musZero.data, axis = 1) return err, idxMaxEst, err[idxMaxEst] def plotEstimator(self, *args, **kwargs): super().plotEstimator(*args, **kwargs) if self.errorEstimatorKind == "NONE": vbMng(self, "MAIN", ("Warning! Error estimator has been arbitrarily normalized " "before plotting."), 15) def greedyNextSample(self, *args, **kwargs) -> Tuple[Np1D, int, float, paramVal]: """Compute next greedy snapshot of solution map.""" RROMPyAssert(self._mode, message = "Cannot add greedy sample.") self.sampleBatchIdx += 1 self.sampleBatchSize = totalDegreeN(self.npar - 1, self.sampleBatchIdx) err, muidx, maxErr, muNext = super().greedyNextSample(*args, **kwargs) if maxErr is not None and (np.any(np.isnan(maxErr)) or np.any(np.isinf(maxErr))): self.sampleBatchIdx -= 1 self.sampleBatchSize = totalDegreeN(self.npar - 1, self.sampleBatchIdx) if (self.errorEstimatorKind == "NONE" and not np.isnan(maxErr) and not np.isinf(maxErr)): maxErr = None return err, muidx, maxErr, muNext def _setSampleBatch(self, maxS:int): self.sampleBatchIdx, self.sampleBatchSize, S = -1, 0, 0 nextBatchSize = 1 while S + nextBatchSize <= maxS: self.sampleBatchIdx += 1 self.sampleBatchSize = nextBatchSize S += self.sampleBatchSize nextBatchSize = totalDegreeN(self.npar - 1, self.sampleBatchIdx + 1) return S def _preliminaryTraining(self): """Initialize starting snapshots of solution map.""" RROMPyAssert(self._mode, message = "Cannot start greedy algorithm.") if self.samplingEngine.nsamples > 0: return self._S = self._setSampleBatch(self.S) super()._preliminaryTraining() self.M, self.N = ("AUTO",) * 2 def setupApproxLocal(self) -> int: """Compute rational interpolant.""" if self.checkComputedApprox(): return -1 RROMPyAssert(self._mode, message = "Cannot setup approximant.") self.verbosity -= 10 vbMng(self, "INIT", "Setting up local approximant.", 5) pMat = self.samplingEngine.projectionMatrix if self.trainedModel is not None: pMat = pMat[:, len(self.trainedModel.data.mus) :] self._setupTrainedModel(pMat, self.trainedModel is not None) self.catchInstability = 2 vbDepth = getVerbosityDepth() unstable = False if self.E > 0: try: - Q = self._setupDenominator()[0] + Q = self._setupDenominator() except RROMPyException as RE: setVerbosityDepth(vbDepth) RROMPyWarning("Downgraded {}: {}".format(RE.__class__.__name__, RE)) unstable = True else: Q = PI() Q.coeffs = np.ones((1,) * self.npar, dtype = np.complex) Q.npar = self.npar Q.polybasis = self.polybasis if not unstable: self.trainedModel.data.Q = copy(Q) try: P = copy(self._setupNumerator()) except RROMPyException as RE: setVerbosityDepth(vbDepth) RROMPyWarning("Downgraded {}: {}".format(RE.__class__.__name__, RE)) unstable = True if not unstable: self.trainedModel.data.P = copy(P) self.trainedModel.data.approxParameters = copy( self.approxParameters) vbMng(self, "DEL", "Done setting up local approximant.", 5) self.catchInstability = 0 self.verbosity += 10 return 1 * unstable def setupApprox(self, plotEst : str = "NONE") -> int: val = super().setupApprox(plotEst) if val == 0: self._setupRational(self.trainedModel.data.Q, self.trainedModel.data.P) self.trainedModel.data.approxParameters = copy( self.approxParameters) return val def loadTrainedModel(self, filename:str): """Load trained reduced model from file.""" super().loadTrainedModel(filename) self._setSampleBatch(self.S + 1) diff --git a/rrompy/reduction_methods/standard/rational_interpolant.py b/rrompy/reduction_methods/standard/rational_interpolant.py index 6cfedf6..1a53788 100644 --- a/rrompy/reduction_methods/standard/rational_interpolant.py +++ b/rrompy/reduction_methods/standard/rational_interpolant.py @@ -1,631 +1,751 @@ # 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, polyTimes, - polyTimesTable, vanderInvTable, - PolynomialInterpolator as PI) + polybases as ppb, polyfitname, + polyvander as pvP, polyTimes, + polyTimesTable, vanderInvTable, + PolynomialInterpolator as PI, + PolynomialInterpolatorNodal as PIN) 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.base.types import (Np1D, Np2D, Tuple, List, sampList, interpEng) from rrompy.utilities.base import verbosityManager as vbMng from rrompy.utilities.numerical import customPInv, dot, potential from rrompy.utilities.numerical.hash_derivative import nextDerivativeIndices from rrompy.utilities.numerical.degree import (reduceDegreeN, degreeTotalToFull, fullDegreeMaxMask, totalDegreeMaxMask) 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; - 'scaleFactorDer': scaling factors for derivative computation; defaults to 'AUTO'; - '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 'AUTO', i.e. maximum allowed; - 'N': degree of rational interpolant denominator; defaults to 'AUTO', i.e. maximum allowed; - 'polydegreetype': type of polynomial degree; defaults to 'TOTAL'; - 'radialDirectionalWeights': radial basis weights for interpolant numerator; defaults to 1; - 'radialDirectionalWeightsAdapt': bounds for adaptive rescaling of radial basis weights; defaults to [-1, -1]; + - 'functionalSolve': strategy for minimization of denominator + functional; allowed values include 'NORM' and 'DOMINANT'; + defaults to 'NORM'; - 'interpRcond': tolerance for interpolation; defaults to None; - 'robustTol': tolerance for robust rational denominator management; defaults to 0; - 'cutOffTolerance': tolerance for ignoring parasitic poles; defaults to np.inf; - 'residueTol': tolerance for residue elimination; defaults to 0., i.e. no bad residues. 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; - 'scaleFactorDer': scaling factors for derivative computation; - '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; - 'radialDirectionalWeightsAdapt': bounds for adaptive rescaling of radial basis weights; + - 'functionalSolve': strategy for minimization of denominator + functional; - 'interpRcond': tolerance for interpolation via numpy.polyfit; - 'robustTol': tolerance for robust rational denominator management; - 'cutOffTolerance': tolerance for ignoring parasitic poles; - 'residueTol': tolerance for residue elimination. 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. scaleFactorDer: Scaling factors for derivative computation. 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. radialDirectionalWeightsAdapt: Bounds for adaptive rescaling of radial basis weights. + functionalSolve: Strategy for minimization of denominator functional. interpRcond: Tolerance for interpolation via numpy.polyfit. robustTol: Tolerance for robust rational denominator management. cutOffTolerance: Tolerance for ignoring parasitic poles. residueTol: Tolerance for residue elimination. 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(["polybasis", "M", "N", "polydegreetype", "radialDirectionalWeights", "radialDirectionalWeightsAdapt", - "interpRcond", "robustTol", - "cutOffTolerance", "residueTol"], + "functionalSolve", "interpRcond", + "robustTol", "cutOffTolerance", + "residueTol"], ["MONOMIAL", "AUTO", "AUTO", "TOTAL", 1., - [-1., -1.], -1, 0., np.inf, 0.]) + [-1., -1.], "NORM", -1, 0., np.inf, 0.]) super().__init__(*args, **kwargs) self.catchInstability = 0 self._postInit() @property def tModelType(self): from .trained_model.trained_model_rational 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: 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 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", "NODAL"]: + raise RROMPyException(("Prescribed functionalSolve not " + "recognized.")) + self._functionalSolve = functionalSolve + except: + RROMPyWarning(("Prescribed functionalSolve not recognized. " + "Overriding to 'NORM'.")) + self._functionalSolve = "NORM" + self._approxParameters["functionalSolve"] = self.functionalSolve + @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): if hasattr(radialDirectionalWeights, "__len__"): radialDirectionalWeights = list(radialDirectionalWeights) else: radialDirectionalWeights = [radialDirectionalWeights] self._radialDirectionalWeights = radialDirectionalWeights self._approxParameters["radialDirectionalWeights"] = ( self.radialDirectionalWeights) @property def radialDirectionalWeightsAdapt(self): """Value of radialDirectionalWeightsAdapt.""" return self._radialDirectionalWeightsAdapt @radialDirectionalWeightsAdapt.setter def radialDirectionalWeightsAdapt(self, radialDirectionalWeightsAdapt): self._radialDirectionalWeightsAdapt = radialDirectionalWeightsAdapt self._approxParameters["radialDirectionalWeightsAdapt"] = ( self.radialDirectionalWeightsAdapt) @property def M(self): """Value of M.""" return self._M @M.setter def M(self, M): if isinstance(M, str): M = M.strip().replace(" ","") if "-" not in M: M = M + "-0" self._M_isauto, self._M_shift = True, int(M.split("-")[-1]) M = 0 if M < 0: raise RROMPyException("M must be non-negative.") self._M = M self._approxParameters["M"] = self.M def _setMAuto(self): self.M = max(0, reduceDegreeN(self.S, self.S, self.npar, self.polydegreetype) - self._M_shift) vbMng(self, "MAIN", "Automatically setting M to {}.".format(self.M), 25) @property def N(self): """Value of N.""" return self._N @N.setter def N(self, N): if isinstance(N, str): N = N.strip().replace(" ","") if "-" not in N: N = N + "-0" self._N_isauto, self._N_shift = True, int(N.split("-")[-1]) N = 0 if N < 0: raise RROMPyException("N must be non-negative.") self._N = N self._approxParameters["N"] = self.N def _setNAuto(self): self.N = max(0, reduceDegreeN(self.S, self.S, self.npar, self.polydegreetype) - self._N_shift) vbMng(self, "MAIN", "Automatically setting N to {}.".format(self.N), 25) @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 cutOffTolerance(self): """Value of cutOffTolerance.""" return self._cutOffTolerance @cutOffTolerance.setter def cutOffTolerance(self, cutOffTolerance): self._cutOffTolerance = cutOffTolerance self._approxParameters["cutOffTolerance"] = self.cutOffTolerance @property def residueTol(self): """Value of residueTol.""" return self._residueTol @residueTol.setter def residueTol(self, residueTol): if residueTol < 0. or (residueTol > 0. and self.npar > 1): RROMPyWarning("Overriding prescribed residue tolerance to 0.") residueTol = 0. self._residueTol = residueTol self._approxParameters["residueTol"] = self.residueTol 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 hasattr(self, "_N_isauto"): self._setNAuto() else: N = reduceDegreeN(self.N, self.S, self.npar, self.polydegreetype) 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: - invD, fitinv = self._computeInterpolantInverseBlocks() + invD, TN = self._computeInterpolantInverseBlocks() + if self.functionalSolve != "NORM" and len(invD) > 1: + RROMPyWarning(("Strategy for functional optimization must be " + "'NORM' for more than one parameter. " + "Overriding to 'NORM'.")) + self.functionalSolve = "NORM" idxSamplesEff = list(range(self.S)) if self.POD: ev, eV = self.findeveVGQR( - self.samplingEngine.RPOD[:, idxSamplesEff], invD) + self.samplingEngine.RPOD[:, idxSamplesEff], invD, TN) else: ev, eV = self.findeveVGExplicit( - self.samplingEngine.samples(idxSamplesEff), invD) + self.samplingEngine.samples(idxSamplesEff), invD, TN) + if self.functionalSolve == "NODAL": break nevBad = checkRobustTolerance(ev, self.robustTol) if nevBad <= 1: break if self.catchInstability > 0: raise RROMPyException(("Instability in denominator " "computation: eigenproblem is poorly " "conditioned."), self.catchInstability == 1) vbMng(self, "MAIN", ("Smallest {} eigenvalues below tolerance. " "Reducing N by 1.").format(nevBad), 10) 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]) + if self.N > 0 and self.functionalSolve == "NODAL": + q = PIN() + q.polybasis, q.nodes = self.polybasis0, eV else: - q.coeffs = eV[:, 0].reshape([self.N + 1] * self.npar) + 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) + else: + q.coeffs = eV.reshape([self.N + 1] * self.npar) vbMng(self, "DEL", "Done computing denominator.", 7) - return q, fitinv + 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.POD: Qevaldiag = Qevaldiag.dot(self.samplingEngine.RPOD.T) if hasattr(self, "_M_isauto"): self._setMAuto() M = self.M else: M = reduceDegreeN(self.M, self.S, self.npar, self.polydegreetype) 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: pParRest = [self.M, self.polybasis, self.verbosity >= 5, self.polydegreetype == "TOTAL", {"derIdxs": self._derIdxs, "reorder": self._reorder, "scl": self.scaleFactorRel}] if self.polybasis in ppb: p = PI() else: self.computeScaleFactor() rDWEff = np.array([w * f for w, f in zip( self.radialDirectionalWeights, self.scaleFactor)]) pParRest = pParRest[: 2] + [rDWEff] + pParRest[2 :] pParRest[-1]["optimizeScalingBounds"] = ( self.radialDirectionalWeightsAdapt) p = RBI() if self.polybasis in ppb + rbpb: pParRest += [{"rcond": self.interpRcond}] wellCond, msg = p.setupByInterpolation(self._musUniqueCN, Qevaldiag, *pParRest) vbMng(self, "MAIN", msg, 5) if wellCond: break if self.catchInstability > 0: raise RROMPyException(("Instability in numerator computation: " "polyfit is poorly conditioned."), self.catchInstability == 1) 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.")) self.M = M 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() self._setupTrainedModel(self.samplingEngine.projectionMatrix) - self._setupRational(self._setupDenominator()[0]) + self._setupRational(self._setupDenominator()) self.trainedModel.data.approxParameters = copy(self.approxParameters) vbMng(self, "DEL", "Done setting up approximant.", 5) return 0 def _setupRational(self, Q:interpEng, P : interpEng = None): vbMng(self, "INIT", "Starting approximant finalization.", 5) computeNum = P is None if self.N > 0 and self.npar == 1: foci = self.sampler.normalFoci() ground = self.sampler.groundPotential() if not np.isinf(self.cutOffTolerance): pls = Q.roots() ground = self.sampler.groundPotential() idKeep = np.logical_and(np.logical_not(np.isinf(pls)), potential(pls, foci) / ground - 1. <= self.cutOffTolerance) if np.sum(idKeep) < self.N: vbMng(self, "MAIN", ("Removing {} poles out of {} due to cut " "off.").format(np.sum(idKeep), self.N), 10) - Q = PI() - Q.npar = self.npar - Q.polybasis = self.polybasis0 - Q.coeffs = np.ones(1, dtype = np.complex) - for pl in pls[idKeep]: - Q.coeffs = polyTimes(Q.coeffs, [- pl, 1.], + if isinstance(Q, PIN): + Q.nodes = Q.nodes[idKeep] + else: + Q = PI() + Q.npar = self.npar + Q.polybasis = self.polybasis0 + Q.coeffs = np.ones(1, dtype = np.complex) + for pl in pls[idKeep]: + Q.coeffs = polyTimes(Q.coeffs, [- pl, 1.], Pbasis = Q.polybasis, Rbasis = Q.polybasis) - Q.coeffs /= np.linalg.norm(Q.coeffs) + Q.coeffs /= np.linalg.norm(Q.coeffs) self.N = np.sum(idKeep) computeNum = True if self.residueTol > 0: if computeNum: P = self._setupNumerator() cfs, pls, _ = rational2heaviside(P, Q) cfs = cfs[: self.N] if self.POD: resEff = np.linalg.norm(cfs, axis = 1) else: resEff = self.HFEngine.norm( self.samplingEngine.projectionMatrix.dot(cfs.T), is_state = self.approx_state) potEff = potential(pls, foci) / ground potEff[np.logical_or(potEff < 1., np.isinf(pls))] = 1. resEff[np.isinf(pls)] = 0. resEff /= potEff idKeep = resEff >= self.residueTol * np.max(resEff) if np.sum(idKeep) < self.N: vbMng(self, "MAIN", ("Removing {} poles out of {} due to residue " "magnitude.").format(np.sum(idKeep), self.N), 10) - Q = PI() - Q.npar = self.npar - Q.polybasis = self.polybasis0 - Q.coeffs = np.ones(1, dtype = np.complex) - for pl in pls[idKeep]: - Q.coeffs = polyTimes(Q.coeffs, [- pl, 1.], + if isinstance(Q, PIN): + Q.nodes = Q.nodes[idKeep] + else: + Q = PI() + Q.npar = self.npar + Q.polybasis = self.polybasis0 + Q.coeffs = np.ones(1, dtype = np.complex) + for pl in pls[idKeep]: + Q.coeffs = polyTimes(Q.coeffs, [- pl, 1.], Pbasis = Q.polybasis, Rbasis = Q.polybasis) - Q.coeffs /= np.linalg.norm(Q.coeffs) - self.N = np.sum(idKeep) + Q.coeffs /= np.linalg.norm(Q.coeffs) + self.N = np.sum(idKeep) else: computeNum = False self.trainedModel.data.Q = Q if computeNum: P = self._setupNumerator() self.trainedModel.data.P = P vbMng(self, "DEL", "Terminated approximant finalization.", 5) 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() pvPPar = [self.polybasis0, self._derIdxs, self._reorder, self.scaleFactorRel] if hasattr(self, "_M_isauto"): self._setMAuto() E = max(self.M, self.N) while E >= 0: if self.polydegreetype == "TOTAL": Eeff = E idxsB = totalDegreeMaxMask(E, self.npar) else: #if self.polydegreetype == "FULL": Eeff = [E] * self.npar idxsB = fullDegreeMaxMask(E, self.npar) TE = pvP(self._musUniqueCN, Eeff, *pvPPar) 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], E, 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 > 0: raise RROMPyException(("Instability in denominator " "computation: polyfit is poorly " "conditioned."), self.catchInstability == 1) EeqN = E == self.N vbMng(self, "MAIN", ("Polyfit is poorly conditioned. Reducing E {}" "by 1.").format("and N " * EeqN), 10) if EeqN: self.N = self.N - 1 E -= 1 if self.N < 0: raise RROMPyException(("Instability in computation of " "denominator. Aborting.")) invD = vanderInvTable(fitinv, idxsB, self._reorder, self._derIdxs) if self.N == E: TN = TE else: if self.polydegreetype == "TOTAL": Neff = self.N idxsB = totalDegreeMaxMask(self.N, self.npar) else: #if self.polydegreetype == "FULL": Neff = [self.N] * self.npar idxsB = fullDegreeMaxMask(self.N, self.npar) TN = pvP(self._musUniqueCN, Neff, *pvPPar) - for k in range(len(invD)): invD[k] = dot(invD[k], TN) - return invD, fitinv + return invD, TN def findeveVGExplicit(self, sampleE:sampList, - invD:List[Np2D]) -> Tuple[Np1D, Np2D]: + invD:List[Np2D], TN: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] + nEnd = TN.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) + if self.functionalSolve == "NODAL": + SEnd = invD[0].shape[1] + G0 = np.zeros((SEnd,) * 2, dtype = np.complex) for k in range(eWidth): - G += dot(dot(gramian, invD[k]).T, invD[k].conj()).T + iDkN = dot(invD[k], TN) + G += dot(dot(gramian, iDkN).T, iDkN.conj()).T + if self.functionalSolve == "NODAL": + G0 += 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) - try: - ev, eV = np.linalg.eigh(G) - except np.linalg.LinAlgError as e: - raise RROMPyException(e) - 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) + if self.functionalSolve == "NORM": + try: + ev, eV = np.linalg.eigh(G) + eV = eV[:, 0] + except np.linalg.LinAlgError as e: + raise RROMPyException(e) + problem = "eigenproblem" + else: # if self.functionalSolve in ["DOMINANT", "NODAL"]: + fitOut = customPInv(G[:-1, :-1], rcond = self.interpRcond, + full = True) + ev = fitOut[1][1][::-1] + eV = np.append(fitOut[0].dot(G[:-1, -1]), -1.) + problem = "linear problem" + if self.functionalSolve == "NODAL": + q = PI() + q.npar, q.polybasis, q.coeffs = self.npar, self.polybasis0, eV + eV, tol, niter, passed = self.findeVNewton(q.roots(), G0) + if passed: + vbMng(self, "MAIN", + ("Newton algorithm for problem of size {} converged " + "(tol = {:.4e}) in {} iterations.").format(nEnd, tol, + niter), 5) + else: + RROMPyWarning(("Newton algorithm for problem of size {} did " + "not converge (tol = {:.4e}) after {} " + "iterations.").format(nEnd, tol, niter)) + else: + vbMng(self, "MAIN", + ("Solved {} of size {} with condition number " + "{:.4e}.").format(problem, nEnd, ev[-1] / ev[0]), 5) return ev, eV - def findeveVGQR(self, RPODE:Np2D, invD:List[Np2D]) -> Tuple[Np1D, Np2D]: + def findeveVGQR(self, RPODE:Np2D, invD:List[Np2D], + TN: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] + nEnd = TN.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) + if self.functionalSolve == "NODAL": + SEnd = invD[0].shape[1] + G0 = np.zeros((SEnd,) * 2, dtype = np.complex) for k in range(eWidth): - Rstack[k * S : (k + 1) * S, :] = dot(RPODE, invD[k]) + Rstack[k * S : (k + 1) * S, :] = dot(RPODE, dot(invD[k], TN)) + if self.functionalSolve == "NODAL": + rDk = dot(RPODE, invD[k]) + G0 += dot(rDk.T.conj(), rDk) vbMng(self, "DEL", "Done building half-gramian.", 10) - vbMng(self, "INIT", "Solving svd for square root of gramian matrix.", - 7) - try: - _, s, eV = np.linalg.svd(Rstack, full_matrices = False) - except np.linalg.LinAlgError as e: - raise RROMPyException(e) - 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) + if self.functionalSolve == "NORM": + try: + _, s, Vh = np.linalg.svd(Rstack, full_matrices = False) + except np.linalg.LinAlgError as e: + raise RROMPyException(e) + ev, eV = s[::-1], Vh[-1, :].conj() + problem = "svd problem" + else: # if self.functionalSolve in ["DOMINANT", "NODAL"]: + fitOut = customPInv(Rstack[:, :-1], rcond = self.interpRcond, + full = True) + ev = fitOut[1][1][::-1] + eV = np.append(fitOut[0].dot(Rstack[:, -1]), -1.) + problem = "linear problem" + if self.functionalSolve == "NODAL": + q = PI() + q.npar, q.polybasis, q.coeffs = self.npar, self.polybasis0, eV + eV, tol, niter, passed = self.findeVNewton(q.roots(), G0) + if passed: + vbMng(self, "MAIN", + ("Newton algorithm for problem of size {} converged " + "(tol = {:.4e}) in {} iterations.").format(nEnd, tol, + niter), 5) + else: + RROMPyWarning(("Newton algorithm for problem of size {} did " + "not converge (tol = {:.4e}) after {} " + "iterations.").format(nEnd, tol, niter)) + else: + vbMng(self, "MAIN", + ("Solved {} of size {} with condition number " + "{:.4e}.").format(problem, nEnd, ev[-1] / ev[0]), 5) return ev, eV + def findeVNewton(self, nodes0:Np1D, gram:Np2D, maxiter : int = 10, + tolNewton : float = 1e-10) -> Tuple[Np1D, float, int, bool]: + N = len(nodes0) + nodes = nodes0 + grad = np.zeros(N, dtype = np.complex) + hess = np.zeros((N, N), dtype = np.complex) + mu = np.repeat(self._musUniqueCN.data[self._reorder], N, axis = 1) + for niter in range(maxiter): + plDist = mu - nodes.reshape(1, -1) + q0, q = np.prod(plDist, axis = 1), [] + for jS in range(N): + mask = np.arange(N) != jS + q += [np.prod(plDist[:, mask], axis = 1)] + grad[jS] = q[-1].conj().dot(gram.dot(q0)) + hess[jS, jS] = q[-1].conj().dot(gram.dot(q[-1])) + for iS in range(jS - 1): + mask = np.logical_and(np.arange(N) != iS, + np.arange(N) != jS) + qij = np.prod(plDist[:, mask], axis = 1) + hij = qij.conj().dot(gram.dot(q0)) + hess[iS, jS] = hij + q[-1].conj().dot(gram.dot(q[iS])) + hess[jS, iS] = hij + q[iS].conj().dot(gram.dot(q[-1])) + dnodes = np.linalg.solve(hess, grad) + nodes += dnodes + tol = np.linalg.norm(dnodes) / np.linalg.norm(nodes) + if tol < tolNewton: break + return nodes, tol, niter, niter < maxiter + def getResidues(self, *args, **kwargs) -> Np2D: """ Obtain approximant residues. Returns: Matrix with residues as columns. """ return self.trainedModel.getResidues(*args, **kwargs) diff --git a/rrompy/reduction_methods/standard/rational_pade.py b/rrompy/reduction_methods/standard/rational_pade.py index 6852778..e97d3ef 100644 --- a/rrompy/reduction_methods/standard/rational_pade.py +++ b/rrompy/reduction_methods/standard/rational_pade.py @@ -1,313 +1,326 @@ # 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 rrompy.reduction_methods.base import checkRobustTolerance from .rational_interpolant import RationalInterpolant from rrompy.utilities.poly_fitting.polynomial import (polybases as ppb, polyfitname, polyvander as pvP, polyTimesTable, vanderInvTable, PolynomialInterpolator as PI) from rrompy.utilities.poly_fitting.radial_basis import (polybases as rbpb, RadialBasisInterpolator as RBI) from rrompy.utilities.base.types import Np2D, Tuple, List from rrompy.utilities.base import verbosityManager as vbMng from rrompy.utilities.numerical import customPInv, dot from rrompy.utilities.numerical.degree import (fullDegreeN, totalDegreeN, reduceDegreeN, degreeTotalToFull, fullDegreeMaxMask, totalDegreeMaxMask) from rrompy.utilities.exception_manager import (RROMPyException, RROMPyAssert, RROMPyWarning) __all__ = ['RationalPade'] class RationalPade(RationalInterpolant): """ ROM rational Pade' 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; - 'scaleFactorDer': scaling factors for derivative computation; defaults to 'AUTO'; - '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 'AUTO', i.e. maximum allowed; - 'N': degree of rational interpolant denominator; defaults to 'AUTO', i.e. maximum allowed; - 'polydegreetype': type of polynomial degree; defaults to 'TOTAL'; - 'radialDirectionalWeights': radial basis weights for interpolant numerator; defaults to 1; - 'radialDirectionalWeightsAdapt': bounds for adaptive rescaling of radial basis weights; defaults to [-1, -1]; + - 'functionalSolve': strategy for minimization of denominator + functional; allowed values include 'NORM' and 'DOMINANT'; + defaults to 'NORM'; - 'interpRcond': tolerance for interpolation; defaults to None; - 'robustTol': tolerance for robust rational denominator management; defaults to 0; - 'cutOffTolerance': tolerance for ignoring parasitic poles; defaults to np.inf; - 'residueTol': tolerance for residue elimination; defaults to 0., i.e. no bad residues. 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; - 'scaleFactorDer': scaling factors for derivative computation; - '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; - 'radialDirectionalWeightsAdapt': bounds for adaptive rescaling of radial basis weights; + - 'functionalSolve': strategy for minimization of denominator + functional; - 'interpRcond': tolerance for interpolation via numpy.polyfit; - 'robustTol': tolerance for robust rational denominator management; - 'cutOffTolerance': tolerance for ignoring parasitic poles; - 'residueTol': tolerance for residue elimination. 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. scaleFactorDer: Scaling factors for derivative computation. 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. radialDirectionalWeightsAdapt: Bounds for adaptive rescaling of radial basis weights. + functionalSolve: Strategy for minimization of denominator functional. interpRcond: Tolerance for interpolation via numpy.polyfit. robustTol: Tolerance for robust rational denominator management. cutOffTolerance: Tolerance for ignoring parasitic poles. residueTol: Tolerance for residue elimination. 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 _setupInterpolationIndices(self): """Setup parameters for polyvander.""" super()._setupInterpolationIndices() if len(self._musUniqueCN) > 1: raise RROMPyException(("Cannot apply centered-like method with " "more than one distinct sample point.")) def _setupDenominator(self): """Compute rational denominator.""" RROMPyAssert(self._mode, message = "Cannot setup denominator.") vbMng(self, "INIT", "Starting computation of denominator.", 7) cfun = totalDegreeN if self.polydegreetype == "TOTAL" else fullDegreeN if hasattr(self, "_N_isauto"): self._setNAuto() else: N = reduceDegreeN(self.N, self.S, self.npar, self.polydegreetype) 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: - invD, fitinv = self._computeInterpolantInverseBlocks() + invD, TN = self._computeInterpolantInverseBlocks() + if self.functionalSolve != "NORM" and len(invD) > 1: + RROMPyWarning(("Strategy for functional optimization must be " + "'NORM' for more than one parameter. " + "Overriding to 'NORM'.")) + self.functionalSolve = "NORM" Seff = cfun(self.N, self.npar) idxSamplesEff = list(range(self.S - Seff, self.S)) if self.POD: ev, eV = self.findeveVGQR( - self.samplingEngine.RPOD[:, idxSamplesEff], invD) + self.samplingEngine.RPOD[:, idxSamplesEff], invD, TN) else: ev, eV = self.findeveVGExplicit( - self.samplingEngine.samples(idxSamplesEff), invD) + self.samplingEngine.samples(idxSamplesEff), invD, TN) nevBad = checkRobustTolerance(ev, self.robustTol) if nevBad <= 1: break if self.catchInstability > 0: raise RROMPyException(("Instability in denominator " "computation: eigenproblem is poorly " "conditioned."), self.catchInstability == 1) 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]) + eV = np.ones(1) + if self.N > 0 and self.functionalSolve == "NODAL": + q = eV else: - q.coeffs = eV[:, 0].reshape([self.N + 1] * self.npar) + 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) + else: + q.coeffs = eV.reshape([self.N + 1] * self.npar) vbMng(self, "DEL", "Done computing denominator.", 7) - return q, fitinv + 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.POD: Qevaldiag = Qevaldiag.dot(self.samplingEngine.RPOD.T) cfun = totalDegreeN if self.polydegreetype == "TOTAL" else fullDegreeN if hasattr(self, "_M_isauto"): self._setMAuto() M = self.M else: M = reduceDegreeN(self.M, self.S, self.npar, self.polydegreetype) 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: Seff = cfun(self.M, self.npar) pParRest = [self.M, self.polybasis, self.verbosity >= 5, self.polydegreetype == "TOTAL", {"derIdxs": [self._derIdxs[0][: Seff]], "reorder": self._reorder[: Seff], "scl": self.scaleFactorRel}] if self.polybasis in ppb: p = PI() else: self.computeScaleFactor() rDWEff = np.array([w * f for w, f in zip( self.radialDirectionalWeights, self.scaleFactor)]) pParRest = pParRest[: 2] + [rDWEff] + pParRest[2 :] pParRest[-1]["optimizeScalingBounds"] = ( self.radialDirectionalWeightsAdapt) p = RBI() if self.polybasis in ppb + rbpb: pParRest += [{"rcond": self.interpRcond}] wellCond, msg = p.setupByInterpolation(self._musUniqueCN, Qevaldiag[: Seff, : Seff], *pParRest) vbMng(self, "MAIN", msg, 5) if wellCond: break if self.catchInstability > 0: raise RROMPyException(("Instability in numerator computation: " "polyfit is poorly conditioned."), self.catchInstability == 1) 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.")) self.M = M vbMng(self, "DEL", "Done computing numerator.", 7) return p 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() if self.polydegreetype == "TOTAL": cfun = totalDegreeN else: cfun = fullDegreeN E = max(self.M, self.N) while E >= 0: Seff = cfun(E, self.npar) pvPPar = [self.polybasis0, [self._derIdxs[0][: Seff]], self._reorder[: Seff], self.scaleFactorRel] if self.polydegreetype == "TOTAL": Eeff = E idxsB = totalDegreeMaxMask(E, self.npar) else: #if self.polydegreetype == "FULL": Eeff = [E] * self.npar idxsB = fullDegreeMaxMask(E, self.npar) TE = pvP(self._musUniqueCN, Eeff, *pvPPar) 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], E, 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 > 0: raise RROMPyException(("Instability in denominator " "computation: polyfit is poorly " "conditioned."), self.catchInstability == 1) EeqN = E == self.N vbMng(self, "MAIN", ("Polyfit is poorly conditioned. Reducing E {}" "by 1.").format("and N " * EeqN), 10) if EeqN: self.N = self.N - 1 E -= 1 if self.N < 0: raise RROMPyException(("Instability in computation of " "denominator. Aborting.")) invD = vanderInvTable(fitinv, idxsB, self._reorder[: Seff], [self._derIdxs[0][: Seff]]) if self.N == E: TN = TE else: if self.polydegreetype == "TOTAL": Neff = self.N idxsB = totalDegreeMaxMask(self.N, self.npar) else: #if self.polydegreetype == "FULL": Neff = [self.N] * self.npar idxsB = fullDegreeMaxMask(self.N, self.npar) TN = pvP(self._musUniqueCN, Neff, *pvPPar) - for k in range(len(invD)): invD[k] = dot(invD[k], TN) - return invD, fitinv + return invD, TN diff --git a/rrompy/utilities/poly_fitting/polynomial/__init__.py b/rrompy/utilities/poly_fitting/polynomial/__init__.py index 8bd3678..b1fab51 100644 --- a/rrompy/utilities/poly_fitting/polynomial/__init__.py +++ b/rrompy/utilities/poly_fitting/polynomial/__init__.py @@ -1,47 +1,49 @@ # 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 from .roots import polyroots from .polynomial_algebra import (changePolyBasis, polyTimes, polyDivide, polyTimesTable, vanderInvTable, blockDiagDer) -from .polynomial_interpolator import PolynomialInterpolator +from .polynomial_interpolator import (PolynomialInterpolator, + PolynomialInterpolatorNodal) __all__ = [ 'polybases', 'polyfitname', 'polydomcoeff', 'polyder', 'polyval', 'polymarginalize', 'polyvander', 'polyroots', 'changePolyBasis', 'polyTimes', 'polyDivide', 'polyTimesTable', 'vanderInvTable', 'blockDiagDer', - 'PolynomialInterpolator' + 'PolynomialInterpolator', + 'PolynomialInterpolatorNodal' ] diff --git a/rrompy/utilities/poly_fitting/polynomial/polynomial_algebra.py b/rrompy/utilities/poly_fitting/polynomial/polynomial_algebra.py index d50df84..303733e 100644 --- a/rrompy/utilities/poly_fitting/polynomial/polynomial_algebra.py +++ b/rrompy/utilities/poly_fitting/polynomial/polynomial_algebra.py @@ -1,146 +1,146 @@ # 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, Tuple, List, interpEng from .vander import polyvander -from .polynomial_interpolator import PolynomialInterpolator from rrompy.utilities.numerical import customPInv from rrompy.utilities.numerical.factorials import multifactorial from rrompy.utilities.numerical.hash_derivative import ( hashDerivativeToIdx as hashD, hashIdxToDerivative as hashI) from rrompy.utilities.exception_manager import RROMPyException __all__ = ['changePolyBasis', 'polyTimes', 'polyDivide', 'polyTimesTable', 'vanderInvTable', 'blockDiagDer'] 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") -> 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: + from .polynomial_interpolator import PolynomialInterpolator if not isinstance(P, PolynomialInterpolator): 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) def vanderInvTable(vanInv:Np2D, idxs:List[int], reorder:List[int], derIdxs:List[List[List[int]]]) -> Np2D: S = len(reorder) Ts = [None] * len(idxs) for k in range(len(idxs)): invLocs = [None] * len(derIdxs) idxGlob = 0 for j, derIdx in enumerate(derIdxs): nder = len(derIdx) idxGlob += nder idxLoc = np.arange(S)[np.logical_and(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) -> Np2D: S = len(reorder) 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)[np.logical_and(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] T[idxLoc[i1], idxLoc[i2]] = val[i3] return T diff --git a/rrompy/utilities/poly_fitting/polynomial/polynomial_interpolator.py b/rrompy/utilities/poly_fitting/polynomial/polynomial_interpolator.py index 5b80325..96d3ffb 100644 --- a/rrompy/utilities/poly_fitting/polynomial/polynomial_interpolator.py +++ b/rrompy/utilities/poly_fitting/polynomial/polynomial_interpolator.py @@ -1,126 +1,221 @@ # 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 scipy.special import factorial as fact +from itertools import combinations from copy import deepcopy as copy from rrompy.utilities.base.types import (List, ListAny, DictAny, Np1D, Np2D, paramList) from rrompy.utilities.base import freepar as fp from rrompy.utilities.poly_fitting.interpolator import GenericInterpolator from rrompy.utilities.poly_fitting.custom_fit import customFit from .base import polyfitname from .val import polyval from .roots import polyroots from .vander import polyvander as pv +from .polynomial_algebra import changePolyBasis, polyTimes from rrompy.utilities.numerical import dot from rrompy.utilities.numerical.degree import degreeTotalToFull from rrompy.utilities.exception_manager import RROMPyAssert, RROMPyException from rrompy.parameter import checkParameterList -__all__ = ['PolynomialInterpolator'] +__all__ = ['PolynomialInterpolator', 'PolynomialInterpolatorNodal'] class PolynomialInterpolator(GenericInterpolator): def __init__(self, other = None): if other is None: return self.coeffs = other.coeffs self.npar = other.npar self.polybasis = other.polybasis @property def shape(self): if self.coeffs.ndim > self.npar: sh = self.coeffs.shape[self.npar :] else: sh = tuple([1]) return sh @property def deg(self): return [x - 1 for x in self.coeffs.shape[: self.npar]] def __getitem__(self, key): return self.coeffs[key] def __call__(self, mu:paramList, der : List[int] = None, scl : Np1D = None): return polyval(mu, self.coeffs, self.polybasis, der, scl) def __copy__(self): return PolynomialInterpolator(self) def __deepcopy__(self, memo): other = PolynomialInterpolator() other.coeffs, other.npar, other.polybasis = copy( (self.coeffs, self.npar, self.polybasis), memo) return other def pad(self, nleft : List[int] = None, nright : List[int] = None): if nleft is None: nleft = [0] * len(self.shape) if nright is None: nright = [0] * len(self.shape) if not hasattr(nleft, "__len__"): nleft = [nleft] if not hasattr(nright, "__len__"): nright = [nright] RROMPyAssert(len(self.shape), len(nleft), "Shape of output") RROMPyAssert(len(self.shape), len(nright), "Shape of output") padwidth = [(0, 0)] * self.npar padwidth = padwidth + [(l, r) for l, r in zip(nleft, nright)] self.coeffs = np.pad(self.coeffs, padwidth, "constant", constant_values = (0., 0.)) def postmultiplyTensorize(self, A:Np2D): RROMPyAssert(A.shape[0], self.shape[-1], "Shape of output") self.coeffs = dot(self.coeffs, A) def setupByInterpolation(self, support:paramList, values:ListAny, deg:int, polybasis : str = "MONOMIAL", verbose : bool = True, totalDegree : bool = True, vanderCoeffs : DictAny = {}, fitCoeffs : DictAny = {}): support = checkParameterList(support) self.npar = support.shape[1] self.polybasis = polybasis if not totalDegree and not hasattr(deg, "__len__"): deg = [deg] * self.npar vander = pv(support, deg, basis = polybasis, **vanderCoeffs) RROMPyAssert(len(vander), len(values), "Number of support values") outDim = values.shape[1:] values = values.reshape(values.shape[0], -1) fitOut = customFit(vander, values, full = True, **fitCoeffs) P = fitOut[0] if verbose: msg = ("Fitting {} samples with degree {} through {}... " "Conditioning of LS system: {:.4e}.").format( len(vander), deg, polyfitname(self.polybasis), fitOut[1][2][0] / fitOut[1][2][-1]) else: msg = None if totalDegree: self.coeffs = degreeTotalToFull(tuple([deg + 1] * self.npar) + outDim, self.npar, P) else: self.coeffs = P.reshape(tuple([d + 1 for d in deg]) + outDim) return fitOut[1][1] == vander.shape[1], msg def roots(self, marginalVals : ListAny = [fp]): RROMPyAssert(self.shape, (1,), "Shape of output") RROMPyAssert(len(marginalVals), self.npar, "Number of parameters") try: rDim = marginalVals.index(fp) if rDim < len(marginalVals) - 1 and fp in marginalVals[rDim + 1 :]: raise except: raise RROMPyException(("Exactly 1 'freepar' entry in " "marginalVals must be provided.")) return polyroots(self.coeffs, self.polybasis, marginalVals) + +class PolynomialInterpolatorNodal(PolynomialInterpolator): + def __init__(self, other = None): + self.npar = 1 + if other is None: return + self.nodes = other.nodes + self.polybasis = other.polybasis + + @property + def nodes(self): + return self._nodes + @nodes.setter + def nodes(self, nodes): + self.coeffs = None + self._nodes = nodes + + @property + def coeffs(self): + if self._coeffs is None: self.buildCoeffs() + return self._coeffs + @coeffs.setter + def coeffs(self, coeffs): + self._coeffs = coeffs + + @property + def shape(self): + return (1,) + + @property + def deg(self): + return [len(self.nodes)] + + def __getitem__(self, key): + return self.coeffs[key] + + def __call__(self, mu:paramList, der : List[int] = None, + scl : Np1D = None): + if der is None: der = 0 + elif isinstance(der, (list,tuple,np.ndarray,)): der = der[0] + if scl is None: scl = 1. + elif isinstance(scl, (list,tuple,np.ndarray,)): scl = scl[0] + mu = checkParameterList(mu, 1) + val = np.zeros(len(mu), dtype = np.complex) + if der == self.deg[0]: + val[:] = 1. + elif der >= 0 and der < self.deg[0]: + plDist = (np.repeat(mu.data, self.deg[0], axis = 1) + - self.nodes.reshape(1, -1)) + for terms in combinations(np.arange(self.deg[0]), + self.deg[0] - der): + val += np.prod(plDist[:, list(terms)], axis = 1) + return scl ** der * fact(der) * val + + def __copy__(self): + return PolynomialInterpolatorNodal(self) + + def __deepcopy__(self, memo): + other = PolynomialInterpolatorNodal() + other.nodes, other.polybasis = copy((self.nodes, self.polybasis), memo) + return other + + def buildCoeffs(self): + local = [np.array([- pl, 1.], dtype = np.complex) for pl in self.nodes] + N = len(local) + while N > 1: + for j in range(N // 2): + local[j] = polyTimes(local[j], local[- 1 - j]) + local = local[(N - 1) // 2 :: -1] + N = len(local) + self._coeffs = changePolyBasis(local[0], None, "MONOMIAL", + self.polybasis) + + def pad(self, nleft : List[int] = None, nright : List[int] = None): + raise RROMPyException(("Padding not allowed for polynomials in nodal " + "form")) + + def postmultiplyTensorize(self, A:Np2D): + raise RROMPyException(("Post-multiply not allowed for polynomials in " + "nodal form")) + + def setupByInterpolation(self, support:paramList, *args, **kwargs): + support = checkParameterList(support) + self.npar = support.shape[1] + if self.npar > 1: + raise RROMPyException(("Polynomial in nodal form must have " + "scalar output")) + output = super().setupByInterpolation(support, *args, **kwargs) + self._nodes = super().roots() + return output + + def roots(self, marginalVals : ListAny = [fp]): + return self.nodes