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