diff --git a/rrompy/reduction_methods/greedy/rational_interpolant_greedy.py b/rrompy/reduction_methods/greedy/rational_interpolant_greedy.py
index 813100e..b72494e 100644
--- a/rrompy/reduction_methods/greedy/rational_interpolant_greedy.py
+++ b/rrompy/reduction_methods/greedy/rational_interpolant_greedy.py
@@ -1,522 +1,521 @@
# 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,
PolynomialInterpolator as PI,
polyvanderTotal as pvT)
from rrompy.utilities.numerical import totalDegreeN, dot
from rrompy.utilities.expression import expressionEvaluator
from rrompy.reduction_methods.standard import RationalInterpolant
from rrompy.utilities.base.types import (Np1D, Tuple, DictAny, HFEng, paramVal,
List)
from rrompy.utilities.base import verbosityManager as vbMng
from rrompy.utilities.poly_fitting import customFit
from rrompy.utilities.exception_manager import (RROMPyWarning, RROMPyException,
RROMPyAssert, RROMPy_FRAGILE)
from rrompy.parameter import checkParameterList
from rrompy.sampling import 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;
- '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;
- 'refinementRatio': ratio of test points to be exhausted before
test set refinement; defaults to 0.2;
- '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', 'INTERPOLATORY',
'LOOK_AHEAD', and 'NONE'; defaults to 'NONE';
- 'interpRcond': tolerance for interpolation; defaults to None;
- 'robustTol': tolerance for robust rational denominator
management; defaults to 0.
Defaults to empty dict.
approx_state(optional): Whether to approximate state. Defaults 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.
- 'greedyTol': uniform error tolerance for greedy algorithm;
- 'collinearityTol': collinearity tolerance for greedy algorithm;
- 'maxIter': maximum number of greedy steps;
- 'refinementRatio': ratio of test points to be exhausted before
test set refinement;
- 'nTestPoints': number of test points;
- 'trainSetGenerator': training sample points generator;
- 'errorEstimatorKind': kind of error estimator;
- 'interpRcond': tolerance for interpolation;
- 'robustTol': tolerance for robust rational denominator
management.
parameterListCritical: Recognized keys of critical approximant
parameters:
- 'S': total number of samples current approximant relies upon;
- 'sampler': sample point generator.
approx_state: Whether to approximate state.
verbosity: Verbosity level.
POD: whether to compute POD of snapshots.
S: number of 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.
refinementRatio: ratio of training points to be exhausted before
training set refinement.
nTestPoints: number of starting training points.
trainSetGenerator: training sample points generator.
robustTol: tolerance for robust rational denominator management.
errorEstimatorKind: kind of error estimator.
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", "INTERPOLATORY",
"LOOK_AHEAD", "NONE"]
def __init__(self, HFEngine:HFEng, mu0 : paramVal = None,
approxParameters : DictAny = {}, approx_state : bool = True,
verbosity : int = 10, timestamp : bool = True):
if not approx_state: RROMPyWarning("Overriding approx_state to True.")
self._preInit()
self._addParametersToList(["errorEstimatorKind"], ["DISCREPANCY"],
toBeExcluded = ["M", "N", "polydegreetype",
"radialDirectionalWeights",
"nNearestNeighbor"])
super().__init__(HFEngine = HFEngine, mu0 = mu0,
approxParameters = approxParameters,
approx_state = True, verbosity = verbosity,
timestamp = timestamp)
self._postInit()
@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."))
@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
def _polyvanderAuxiliary(self, mus, deg, *args):
return pvT(mus, deg, *args)
def getErrorEstimatorDiscrepancy(self, mus:Np1D) -> Np1D:
"""Discrepancy-based residual estimator."""
mus = checkParameterList(mus, self.npar)[0]
muCTest = self.trainedModel.centerNormalize(mus)
verb = self.trainedModel.verbosity
self.trainedModel.verbosity = 0
QTest = self.trainedModel.getQVal(mus)
checkIfAffine(self.HFEngine, "apply discrepancy-based error estimator")
self.HFEngine.buildA()
self.HFEngine.buildb()
nAs, nbs = self.HFEngine.nAs, self.HFEngine.nbs
muTrainEff = self.mus ** self.HFEngine.rescalingExp
muTestEff = mus ** self.HFEngine.rescalingExp
PTrain = self.trainedModel.getPVal(self.mus).data.T
QTrain = self.trainedModel.getQVal(self.mus)
PTest = self.trainedModel.getPVal(mus).data
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
self.trainedModel.verbosity = verb
return err
def getErrorEstimatorInterpolatory(self, mus:Np1D) -> Np1D:
"""Interpolation-based residual estimator."""
errTest, QTest, idxMaxEst = self._EIMStep(mus)
self.initEstimatorNormEngine()
mu_muTestSample = mus[idxMaxEst]
app_muTestSample = self.getApproxReduced(mu_muTestSample)
if self._mode == RROMPy_FRAGILE:
if not self.HFEngine.isCEye:
raise RROMPyException(("Cannot compute INTERPOLATORY residual "
"estimator in fragile mode for "
"non-scalar C."))
app_muTestSample = dot(self.trainedModel.data.projMat,
app_muTestSample.data)
else:
app_muTestSample = dot(self.samplingEngine.samples,
app_muTestSample)
resmu = self.HFEngine.residual(mu_muTestSample, app_muTestSample,
post_c = False)
RHSmu = self.HFEngine.residual(mu_muTestSample, None, post_c = False)
ressamples = (self.estimatorNormEngine.norm(resmu)
/ self.estimatorNormEngine.norm(RHSmu))
musT = copy(self.mus)
musT.append(mu_muTestSample)
musT = self.trainedModel.centerNormalize(musT)
musC = self.trainedModel.centerNormalize(mus)
resT = np.zeros(len(musT), dtype = np.complex)
err = np.zeros(len(mus))
for l in range(len(mu_muTestSample)):
resT[len(self.mus) + l] = ressamples[l] * QTest[idxMaxEst[l]]
p = PI()
wellCond, msg = p.setupByInterpolation(musT, resT, self.E + 1,
self.polybasis, self.verbosity >= 15,
True, {}, {"rcond": self.interpRcond})
err += np.abs(p(musC))
resT[len(self.mus) + l] = 0.
err /= QTest
vbMng(self, "MAIN", msg, 15)
return err
def getErrorEstimatorLookAhead(self, mus:Np1D) -> Tuple[Np1D, List[int]]:
"""Residual estimator based on look-ahead idea."""
errTest, QTest, idxMaxEst = self._EIMStep(mus)
self.initEstimatorNormEngine()
mu_muTestSample = mus[idxMaxEst]
app_muTestSample = self.getApproxReduced(mu_muTestSample)
if self._mode == RROMPy_FRAGILE:
app_muTestSample = dot(self.trainedModel.data.projMat,
app_muTestSample.data)
else:
app_muTestSample = dot(self.samplingEngine.samples,
app_muTestSample)
for j, mu in enumerate(mu_muTestSample):
self.samplingEngine.nextSample(mu)
if hasattr(self.samplingEngine, "samples_full"):
uEx = self.samplingEngine.samples_full[-1]
else:
uEx = self.samplingEngine.samples[-1]
if j == 0:
solmu = emptySampleList()
solmu.reset((len(uEx), len(mu_muTestSample)),
dtype = uEx.dtype)
solmu[j] = uEx
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), dtype = np.complex)
err = np.zeros(len(mus))
for l in range(len(mu_muTestSample)):
errT[len(self.mus) + l] = errsamples[l] * QTest[idxMaxEst[l]]
p = PI()
wellCond, msg = p.setupByInterpolation(musT, errT, self.E + 1,
self.polybasis, self.verbosity >= 15,
True, {}, {"rcond": self.interpRcond})
err += np.abs(p(musC))
errT[len(self.mus) + l] = 0.
err /= QTest
vbMng(self, "MAIN", msg, 15)
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 = checkParameterList(mus, self.npar)[0]
- muCTest = self.trainedModel.centerNormalize(mus)
verb = self.trainedModel.verbosity
self.trainedModel.verbosity = 0
QTest = self.trainedModel.getQVal(mus)
QTest = np.abs(QTest)
muCTest = self.trainedModel.centerNormalize(mus)
muCTrain = self.trainedModel.centerNormalize(self.mus)
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]
coeffTest = np.linalg.solve(vanTrial, valuesTrial)
errTest = (np.abs(vanTestNextEff - vanTestEff.dot(coeffTest))
/ np.expand_dims(QTest, 1))
if only_one:
self.trainedModel.verbosity = verb
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])
self.trainedModel.verbosity = verb
return errTest, QTest, idxMaxEst
def errorEstimator(self, mus:Np1D, return_max : bool = False) -> Np1D:
"""Standard residual-based error estimator."""
setupOK = self.setupApproxLocal()
if not setupOK:
err = np.empty(len(mus))
err[:] = np.nan
if not return_max: return err
return err, 0, np.nan
mus = checkParameterList(mus, self.npar)[0]
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 == "INTERPOLATORY":
err = self.getErrorEstimatorInterpolatory(mus)
elif self.errorEstimatorKind == "LOOK_AHEAD":
err, idxMaxEst = self.getErrorEstimatorLookAhead(mus)
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 != "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 greedyNextSample(self, muidx:int, plotEst : bool = False)\
-> 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(muidx, plotEst)
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 _preliminaryTraining(self):
"""Initialize starting snapshots of solution map."""
RROMPyAssert(self._mode, message = "Cannot start greedy algorithm.")
if self.samplingEngine.nsamples > 0:
return
S = self.S
self.sampleBatchIdx, self.sampleBatchSize, self._S = -1, 0, 0
nextBatchSize = 1
while self._S + nextBatchSize <= S:
self.sampleBatchIdx += 1
self.sampleBatchSize = nextBatchSize
self._S += self.sampleBatchSize
nextBatchSize = totalDegreeN(self.npar - 1,
self.sampleBatchIdx + 1)
super()._preliminaryTraining()
def setupApproxLocal(self):
"""
Compute rational interpolant.
SVD-based robust eigenvalue management.
"""
if self.checkComputedApprox():
return True
RROMPyAssert(self._mode, message = "Cannot setup approximant.")
self.verbosity -= 10
vbMng(self, "INIT", "Setting up local approximant.", 5)
self._N, self._M = [self.E] * 2
pMat = self.samplingEngine.samples.data
pMatEff = dot(self.HFEngine.C, pMat)
if self.trainedModel is None:
self.trainedModel = self.tModelType()
self.trainedModel.verbosity = self.verbosity
self.trainedModel.timestamp = self.timestamp
datadict = {"mu0": self.mu0, "projMat": pMatEff,
"scaleFactor": self.scaleFactor,
"rescalingExp": self.HFEngine.rescalingExp}
self.trainedModel.data = self.initializeModelData(datadict)[0]
else:
self.trainedModel = self.trainedModel
self.trainedModel.data.projMat = copy(pMatEff)
self.trainedModel.data.mus = copy(self.mus)
self.trainedModel.data.mus = copy(self.mus)
self.catchInstability = 2
stable = True
if self.N > 0:
try:
Q = self._setupDenominator()[0]
except RROMPyException as RE:
RROMPyWarning(RE._msg)
vbMng(self, "DEL", "", 7)
stable = False
else:
Q = PI()
Q.coeffs = np.ones((1,) * self.npar, dtype = np.complex)
Q.npar = self.npar
Q.polybasis = self.polybasis
if stable:
self.trainedModel.data.Q = copy(Q)
try:
P = copy(self._setupNumerator())
except RROMPyException as RE:
RROMPyWarning(RE._msg)
vbMng(self, "DEL", "", 7)
stable = False
if stable:
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 stable
def loadTrainedModel(self, filename:str):
"""Load trained reduced model from file."""
super().loadTrainedModel(filename)
self.sampleBatchIdx, self.sampleBatchSize, _S = -1, 0, 0
nextBatchSize = 1
while _S + nextBatchSize <= self.S + 1:
self.sampleBatchIdx += 1
self.sampleBatchSize = nextBatchSize
_S += self.sampleBatchSize
nextBatchSize = totalDegreeN(self.npar - 1,
self.sampleBatchIdx + 1)
diff --git a/rrompy/reduction_methods/pivoted/rational_interpolant_greedy_pivoted.py b/rrompy/reduction_methods/pivoted/rational_interpolant_greedy_pivoted.py
index da79089..678ff49 100644
--- a/rrompy/reduction_methods/pivoted/rational_interpolant_greedy_pivoted.py
+++ b/rrompy/reduction_methods/pivoted/rational_interpolant_greedy_pivoted.py
@@ -1,453 +1,470 @@
# 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 GenericPivotedApproximant
from rrompy.reduction_methods.greedy.rational_interpolant_greedy import (
RationalInterpolantGreedy)
from rrompy.reduction_methods.greedy.generic_greedy_approximant import (
pruneSamples)
from rrompy.utilities.base.types import Np1D, HFEng, DictAny, ListAny, paramVal
from rrompy.utilities.base import verbosityManager as vbMng
from rrompy.utilities.numerical import totalDegreeN, dot
from rrompy.utilities.poly_fitting.polynomial import polyvander as pv
from rrompy.utilities.exception_manager import RROMPyAssert
from rrompy.parameter import emptyParameterList, checkParameterList
__all__ = ['RationalInterpolantGreedyPivoted']
class RationalInterpolantGreedyPivoted(GenericPivotedApproximant,
RationalInterpolantGreedy):
"""
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;
- 'matchingWeight': weight for pole matching optimization; defaults
to 1;
- 'cutOffTolerance': tolerance for ignoring parasitic poles;
defaults to np.inf;
- 'cutOffType': rule for tolerance computation for parasitic poles;
defaults to 'MAGNITUDE';
- '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'
and 'LEGENDRE'; 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;
- 'refinementRatio': ratio of test points to be exhausted before
test set refinement; defaults to 0.2;
- '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', 'INTERPOLATORY',
'LOOK_AHEAD', and 'NONE'; defaults to 'NONE';
- 'MMarginal': degree of marginal interpolant; defaults to 0;
- 'polydegreetypeMarginal': type of polynomial degree for marginal;
defaults to 'TOTAL';
- 'radialDirectionalWeightsMarginal': radial basis weights for
marginal interpolant; defaults to 0, i.e. identity;
- 'nNearestNeighborMarginal': number of marginal nearest neighbors
considered if polybasisMarginal allows; defaults to -1;
- 'interpRcond': tolerance for pivot interpolation; defaults to
None;
- 'interpRcondMarginal': tolerance for marginal interpolation;
defaults to None;
- 'robustTol': tolerance for robust rational denominator
management; defaults to 0.
Defaults to empty dict.
approx_state(optional): Whether to approximate state. Defaults to
False.
verbosity(optional): Verbosity level. Defaults to 10.
Attributes:
HFEngine: HF problem solver.
mu0: Default parameter.
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;
- 'matchingWeight': weight for pole matching optimization;
- 'cutOffTolerance': tolerance for ignoring parasitic poles;
- 'cutOffType': rule for tolerance computation for parasitic poles;
- 'polybasis': type of polynomial basis for pivot
interpolation;
- 'polybasisMarginal': type of polynomial basis for marginal
interpolation;
- 'greedyTol': uniform error tolerance for greedy algorithm;
- 'collinearityTol': collinearity tolerance for greedy algorithm;
- 'maxIter': maximum number of greedy steps;
- 'refinementRatio': ratio of test points to be exhausted before
test set refinement;
- 'nTestPoints': number of test points;
- 'trainSetGenerator': training sample points generator;
- 'errorEstimatorKind': kind of error estimator;
- 'MMarginal': degree of marginal interpolant;
- 'polydegreetypeMarginal': type of polynomial degree for marginal;
- 'radialDirectionalWeights': radial basis weights for pivot
numerator;
- 'radialDirectionalWeightsMarginal': radial basis weights for
marginal interpolant;
- 'nNearestNeighbor': number of pivot nearest neighbors considered
if polybasis allows;
- 'nNearestNeighborMarginal': number of marginal nearest neighbors
considered if polybasisMarginal allows;
- 'interpRcond': tolerance for pivot interpolation;
- 'interpRcondMarginal': tolerance for marginal interpolation;
- 'robustTol': tolerance for robust rational denominator
management;
- 'correctorForce': whether corrector should forcefully delete bad
poles;
- 'correctorTol': tolerance for corrector step;
- 'correctorMaxIter': maximum number of corrector iterations.
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.
matchingWeight: Weight for pole matching optimization.
cutOffTolerance: Tolerance for ignoring parasitic poles.
cutOffType: Rule for tolerance computation for parasitic poles.
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.
greedyTol: uniform error tolerance for greedy algorithm.
collinearityTol: Collinearity tolerance for greedy algorithm.
maxIter: maximum number of greedy steps.
refinementRatio: ratio of training points to be exhausted before
training set refinement.
nTestPoints: number of starting training points.
trainSetGenerator: training sample points generator.
errorEstimatorKind: kind of error estimator.
MMarginal: Degree of marginal interpolant.
polydegreetypeMarginal: Type of polynomial degree for marginal.
radialDirectionalWeightsMarginal: Radial basis weights for marginal
interpolant.
nNearestNeighborMarginal: Number of marginal nearest neighbors
considered if polybasisMarginal allows.
interpRcond: Tolerance for pivot interpolation.
interpRcondMarginal: Tolerance for marginal interpolation.
robustTol: Tolerance for robust rational denominator management.
muBoundsPivot: 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.
"""
def __init__(self, HFEngine:HFEng, mu0 : paramVal = None,
directionPivot : ListAny = [0],
approxParameters : DictAny = {}, approx_state : bool = False,
verbosity : int = 10, timestamp : bool = True):
self._preInit()
self._addParametersToList(toBeExcluded = ["sampler"])
super().__init__(HFEngine = HFEngine, mu0 = mu0,
directionPivot = directionPivot,
approxParameters = approxParameters,
approx_state = approx_state, verbosity = verbosity,
timestamp = timestamp)
self._postInit()
@property
def tModelType(self):
if hasattr(self, "_temporaryPivot"): return super().tModelType
from rrompy.reduction_methods.trained_model import \
TrainedModelPivotedRational
return TrainedModelPivotedRational
@property
def polybasis0(self):
if "_" in self.polybasis:
return self.polybasis.split("_")[0]
return self.polybasis
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_HFErescalingExp = copy(self.HFEngine.rescalingExp)
self._mu0 = checkParameterList(self.mu0(self.directionPivot), 1)[0]
self._mus = checkParameterList(self.mus(self.directionPivot), 1)[0]
self.HFEngine.rescalingExp = [self.HFEngine.rescalingExp[
self.directionPivot[0]]]
else:
self._mu0 = self._m_mu0
self._mus = self._m_selfmus
self.HFEngine.rescalingExp = self._m_HFErescalingExp
del self._m_mu0, self._m_selfmus, self._m_HFErescalingExp
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.rescalingExp = self.HFEngine.rescalingExp
Qc = np.zeros((len(self.trainedModel.data.Q.coeffs),) * self.npar,
dtype = self.trainedModel.data.Q.coeffs.dtype)
Pc = np.zeros((len(self.trainedModel.data.P.coeffs),) * self.npar
+ (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 = checkParameterList(musUniqueCNAux,
+ self.npar)[0]
+ 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 = checkParameterList(
self.mu0(self.directionPivot), 1)[0]
self.trainedModel.data.scaleFactor = self.scaleFactor
self.trainedModel.data.rescalingExp = self.HFEngine.rescalingExp[
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 not setupOK:
err = np.empty(len(mus))
err[:] = np.nan
if not return_max: return err
return err, 0, 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.")
S = self.S
self.sampleBatchIdx, self.sampleBatchSize, self._S = -1, 0, 0
nextBatchSize = 1
while self._S + nextBatchSize <= S:
self.sampleBatchIdx += 1
self.sampleBatchSize = nextBatchSize
self._S += self.sampleBatchSize
nextBatchSize = totalDegreeN(self.npar - 1,
self.sampleBatchIdx + 1)
self.resetSamples()
musPivot = self.trainSetGenerator.generatePoints(self.S)[
list(range(self.S))]
musPivot = checkParameterList(musPivot, 1)[0]
muTestPivot = self.samplerPivot.generatePoints(self.nTestPoints)
idxPop = pruneSamples(muTestPivot ** self.HFEngine.rescalingExp[
self.directionPivot[0]],
musPivot ** self.HFEngine.rescalingExp[
self.directionPivot[0]],
1e-10 * self.scaleFactor[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):
self.mus.data[k, self.directionPivot] = musPivot[k].data
self.mus.data[k, self.directionMarginal] = self.musMargLoc.data
for k in range(len(muTestPivot)):
muTestBase.data[k, self.directionPivot] = muTestPivot[k].data
muTestBase.data[k, self.directionMarginal] = self.musMargLoc.data
muTestBase.pop(idxPop)
muTestBase = muTestBase.sort()
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), 2)
self.samplingEngine.iterSample(self.mus)
self._S = len(self.mus)
self._approxParameters["S"] = self.S
self.muTest = emptyParameterList()
self.muTest.reset((len(muTestBase) + 1, self.mus.shape[1]))
self.muTest.data[: -1] = muTestBase.data
self.muTest.data[-1] = muLast.data
def _enrichTestSet(self, nTest:int):
"""Add extra elements to test set."""
RROMPyAssert(self._mode, message = "Cannot enrich test set.")
muTestExtra = self.samplerPivot.generatePoints(2 * nTest)
muTotal = copy(self.mus)
muTotal.append(self.muTest)
idxPop = pruneSamples(muTestExtra ** self.HFEngine.rescalingExp[
self.directionPivot[0]],
muTotal ** self.HFEngine.rescalingExp[
self.directionPivot[0]],
1e-10 * self.scaleFactor[0])
muTestExtra.pop(idxPop)
muTestNew = np.empty((len(self.muTest) + len(muTestExtra),
self.muTest.shape[1]), dtype = np.complex)
muTestNew[: len(self.muTest)] = self.muTest.data
muTestNew[len(self.muTest) :] = muTestExtra.data
self.muTest = checkParameterList(muTestNew, self.npar)[0].sort()
vbMng(self, "MAIN",
"Enriching test set by {} elements.".format(len(muTestExtra)), 5)
def setupApprox(self, plotEst : bool = False):
"""Compute rational interpolant."""
if self.checkComputedApprox():
return True
RROMPyAssert(self._mode, message = "Cannot setup approximant.")
vbMng(self, "INIT", "Setting up {}.". format(self.name()), 5)
self.musMarginal = self.samplerMarginal.generatePoints(self.SMarginal)
S0 = copy(self.S)
Qs, Ps = [], []
self.computeScaleFactor()
self._scaleFactorOldPivot = copy(self.scaleFactor)
self.scaleFactor = self.scaleFactorPivot
nparEff = self.npar
self._temporaryPivot = 1
samplingEngs = []
for j in range(len(self.musMarginal)):
self._S = S0
self.musMargLoc = self.musMarginal[j]
RationalInterpolantGreedy.setupSampling(self)
self.trainedModel = None
self.verbosity -= 5
super().setupApprox(plotEst)
self.verbosity += 5
samplingEngs += [copy(self.samplingEngine)]
Qs = Qs + [copy(self.trainedModel.data.Q)]
Ps = Ps + [copy(self.trainedModel.data.P)]
del self.musMargLoc
# other finalization instructions
self.setupSampling()
self.samplingEngine.resetHistory(len(self.musMarginal))
for j in range(len(self.musMarginal)):
self.samplingEngine.setsample(samplingEngs[j].samples, j, False)
self.samplingEngine.mus[j] = copy(samplingEngs[j].mus)
self.samplingEngine.nsamples[j] = samplingEngs[j].nsamples
self.samplingEngine.postprocessuBulk(j)
if j == 0:
self._mus = checkParameterList(samplingEngs[j].mus, nparEff)[0]
else:
self._mus.append(samplingEngs[j].mus)
if self.POD:
self.samplingEngine.coalesceSamples(self.interpRcondMarginal)
postR = self.samplingEngine.RPODCoalesced
else:
self.samplingEngine.coalesceSamples()
postR = np.eye(self.samplingEngine.samplesCoalesced.shape[1])
# update Ps by post-multiplication with suitable matrix
idxCurr = 0
for j in range(len(self.musMarginal)):
nsj = Ps[j].coeffs.shape[1]
Ps[j].coeffs = dot(Ps[j].coeffs,
postR[:, idxCurr : idxCurr + nsj].T)
idxCurr = idxCurr + nsj
self.scaleFactor = self._scaleFactorOldPivot
del self._scaleFactorOldPivot, self._temporaryPivot
pMat = self.samplingEngine.samplesCoalesced.data
pMatEff = dot(self.HFEngine.C, pMat)
self.trainedModel = self.tModelType()
self.trainedModel.verbosity = self.verbosity
self.trainedModel.timestamp = self.timestamp
datadict = {"mu0": self.mu0, "projMat": pMatEff,
"scaleFactor": self.scaleFactor,
"rescalingExp": self.HFEngine.rescalingExp,
"directionPivot": self.directionPivot}
self.trainedModel.data = self.initializeModelData(datadict)[0]
self.trainedModel.data.mus = copy(self.mus)
self.trainedModel.data.musMarginal = copy(self.musMarginal)
self.trainedModel.data.marginalInterp = self._setupMarginalInterp()
self.trainedModel.data.Qs, self.trainedModel.data.Ps = Qs, Ps
vbMng(self, "INIT", "Matching poles.", 10)
self.trainedModel.initializeFromRational(self.HFEngine,
self.matchingWeight, self.POD,
self.approx_state)
vbMng(self, "DEL", "Done matching poles.", 10)
vbMng(self, "INIT", "Recompressing by cut-off.", 10)
msg = self.trainedModel.recompressByCutOff([-1., 1.],
self.cutOffTolerance,
self.cutOffType)
vbMng(self, "DEL", "Done recompressing." + msg, 10)
self.trainedModel.data.approxParameters = copy(self.approxParameters)
vbMng(self, "DEL", "Done setting up approximant.", 5)
diff --git a/rrompy/utilities/poly_fitting/polynomial/vander.py b/rrompy/utilities/poly_fitting/polynomial/vander.py
index 60bb47d..3f218d9 100644
--- a/rrompy/utilities/poly_fitting/polynomial/vander.py
+++ b/rrompy/utilities/poly_fitting/polynomial/vander.py
@@ -1,137 +1,138 @@
# 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 .der import polyder
from rrompy.utilities.base.types import Np1D, Np2D, List, paramList
from rrompy.utilities.numerical import totalDegreeSet
from rrompy.parameter import checkParameterList
from rrompy.utilities.exception_manager import RROMPyException, RROMPyAssert
__all__ = ['polyvander', 'polyvanderTotal']
def firstDerTransition(dim:int, TDirac:List[Np2D], basis:str,
scl : Np1D = None) -> Np2D:
C_m = np.zeros((dim, len(TDirac), len(TDirac)), dtype = float)
for j, Tj in enumerate(TDirac):
m, om = [0] * dim, [(0, 0)] * dim
for idx in range(dim):
m[idx], om[idx] = 1, (0, 1)
J_der = polyder(Tj, basis, m, scl)
- C_m[idx, :, j] = np.ravel(np.pad(J_der, mode = "constant",
- pad_width = om))
+ if J_der.size != len(TDirac):
+ J_der = np.pad(J_der, mode = "constant", pad_width = om)
+ C_m[idx, :, j] = np.ravel(J_der)
m[idx], om[idx] = 0, (0, 0)
return C_m
def countDerDirections(n:int, base:int, digits:int, idx:int):
if digits == 0: return []
dig = n % base
return [(idx, dig)] * (dig > 0) + countDerDirections(
(n - dig) // base, base, digits - 1, idx + 1)
def polyvander(x:paramList, degs:List[int], basis:str,
derIdxs : List[List[List[int]]] = None,
reorder : List[int] = None, scl : Np1D = None) -> Np2D:
"""
Compute full Hermite-Vandermonde matrix with specified derivative
directions.
E.g. assume that we want to obtain the Vandermonde matrix for
(value, derx, derx2) at x = [0, 0],
(value, dery) at x = [1, 0],
(dery, derxy) at x = [0, 0],
of degree 3 in x and 1 in y, using Chebyshev polynomials.
This can be done by
polyvander([[0, 0], [1, 0]], # unique sample points
[3, 1], # polynomial degree
"chebyshev", # polynomial family
[
[[0, 0], [1, 0], [2, 0], [0, 1], [1, 1]],
# derivative directions at first point
[[0, 0], [0, 1]] # derivative directions at second point
],
[0, 1, 2, 5, 6, 3, 4] # reorder indices
)
"""
if not isinstance(degs, (list,tuple,np.ndarray,)): degs = [degs]
dim = len(degs)
x = checkParameterList(x, dim)[0]
x_un, idx_un = x.unique(return_inverse = True)
if len(x_un) < len(x):
raise RROMPyException("Sample points must be distinct.")
del x_un
try:
vanderbase = {"CHEBYSHEV" : np.polynomial.chebyshev.chebvander,
"LEGENDRE" : np.polynomial.legendre.legvander,
"MONOMIAL" : np.polynomial.polynomial.polyvander
}[basis.upper()]
except:
raise RROMPyException("Polynomial basis not recognized.")
VanBase = vanderbase(x(0), degs[0])
for j in range(1, dim):
VNext = vanderbase(x(j), degs[j])
for jj in range(j): VNext = np.expand_dims(VNext, 1)
VanBase = VanBase[..., None] * VNext
VanBase = VanBase.reshape((len(x), -1))
if derIdxs is None or VanBase.shape[-1] <= 1:
Van = VanBase
else:
derFlat, idxRep = [], []
for j, derIdx in enumerate(derIdxs):
derFlat += derIdx[:]
idxRep += [j] * len(derIdx[:])
for j in range(len(derFlat)):
if not hasattr(derFlat[j], "__len__"):
derFlat[j] = [derFlat[j]]
RROMPyAssert(len(derFlat[j]), dim, "Number of dimensions")
TDirac = [y.reshape([d + 1 for d in degs])
for y in np.eye(VanBase.shape[-1], dtype = int)]
Cs_loc = firstDerTransition(dim, TDirac, basis, scl)
Van = np.empty((len(derFlat), VanBase.shape[-1]),
dtype = VanBase.dtype)
for j in range(len(derFlat)):
Van[j, :] = VanBase[idxRep[j], :]
for k in range(dim):
for der in range(derFlat[j][k]):
Van[j, :] = Van[j, :].dot(Cs_loc[k]) / (der + 1)
if reorder is not None: Van = Van[reorder, :]
return Van
def polyvanderTotal(x:paramList, deg:int, basis:str,
derIdxs : List[List[List[int]]] = None,
reorder : List[int] = None, scl : Np1D = None) -> Np2D:
"""
Compute full total degree Hermite-Vandermonde matrix with specified
derivative directions.
"""
x = checkParameterList(x)[0]
degs = [deg] * x.shape[1]
VanBase = polyvander(x, degs, basis, derIdxs, reorder, scl)
derIdxs, mask = totalDegreeSet(deg, x.shape[1], return_mask = True)
ordIdxs = np.empty(len(derIdxs), dtype = int)
derTotal = np.array([np.sum(y) for y in derIdxs])
idxPrev = 0
rangeAux = np.arange(len(derIdxs))
for j in range(deg + 1):
idxLocal = rangeAux[derTotal == j][::-1]
idxPrev += len(idxLocal)
ordIdxs[idxPrev - len(idxLocal) : idxPrev] = idxLocal
return VanBase[:, mask][:, ordIdxs]
diff --git a/tests/reduction_methods_multiD/pivoted_rational_2d.py b/tests/reduction_methods_multiD/pivoted_rational_2d.py
new file mode 100644
index 0000000..647c876
--- /dev/null
+++ b/tests/reduction_methods_multiD/pivoted_rational_2d.py
@@ -0,0 +1,116 @@
+# 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 matrix_random import matrixRandom
+from rrompy.reduction_methods.pivoted import (
+ RationalInterpolantPivoted as RIP,
+ RationalInterpolantGreedyPivoted as RIGP)
+from rrompy.parameter.parameter_sampling import (QuadratureSampler as QS,
+ ManualSampler as MS)
+
+def test_pivoted_uniform():
+ mu = [5.05, 7.1]
+ mu0 = [5., 7.]
+ solver = matrixRandom()
+ uh = solver.solve(mu)[0]
+ params = {"POD": True, "M": 4, "N": 4, "S": 5, "polybasis": "CHEBYSHEV",
+ "samplerPivot": QS([4.75, 5.25], "CHEBYSHEV"),
+ "MMarginal": 4, "SMarginal": 5, "polybasisMarginal": "MONOMIAL",
+ "samplerMarginal": QS([6.75, 7.25], "UNIFORM"),
+ "matchingWeight": 1.}
+ approx = RIP(solver, mu0, [0], params, verbosity = 0)
+ approx.setupApprox()
+
+ uhP1 = approx.getApprox(mu)[0]
+ errP = approx.getErr(mu)[0]
+ errNP = approx.normErr(mu)[0]
+ myerrP = uhP1 - uh
+ assert np.allclose(np.abs(errP - myerrP), 0., rtol = 1e-3)
+ assert np.isclose(solver.norm(errP), errNP, rtol = 1e-3)
+ resP = approx.getRes(mu)[0]
+ resNP = approx.normRes(mu)
+ assert np.isclose(solver.norm(resP), resNP, rtol = 1e-3)
+ assert np.allclose(np.abs(resP - (solver.b(mu) - solver.A(mu).dot(uhP1))),
+ 0., rtol = 1e-3)
+ assert np.isclose(errNP / solver.norm(uh), 6.0631706e-04, rtol = 1)
+
+def test_pivoted_manual_grid(capsys):
+ mu = [5.05, 7.1]
+ mu0 = [5., 7.]
+ solver = matrixRandom()
+ uh = solver.solve(mu)[0]
+ params = {"POD": False, "M": 4, "N": 4, "S": 5, "polybasis": "MONOMIAL",
+ "samplerPivot": MS([4.75, 5.25], np.array([5.])),
+ "MMarginal": 4, "SMarginal": 5, "polybasisMarginal": "MONOMIAL",
+ "samplerMarginal": MS([6.75, 7.25], np.linspace(6.75, 7.25, 5)),
+ "matchingWeight": 1., "robustTol": 1e-6, "interpRcond": 1e-3,
+ "cutOffTolerance": 1., "cutOffType": "MAGNITUDE"}
+ approx = RIP(solver, mu0, [0], params, verbosity = 0)
+ approx.setupApprox()
+
+ uhP1 = approx.getApprox(mu)[0]
+ errP = approx.getErr(mu)[0]
+ errNP = approx.normErr(mu)[0]
+ myerrP = uhP1 - uh
+ assert np.allclose(np.abs(errP - myerrP), 0., rtol = 1e-3)
+ assert np.isclose(solver.norm(errP), errNP, rtol = 1e-3)
+ resP = approx.getRes(mu)[0]
+ resNP = approx.normRes(mu)
+ assert np.isclose(solver.norm(resP), resNP, rtol = 1e-3)
+ assert np.allclose(np.abs(resP - (solver.b(mu) - solver.A(mu).dot(uhP1))),
+ 0., rtol = 1e-3)
+ assert np.isclose(errNP / solver.norm(uh), .4763489, rtol = 1)
+
+ out, err = capsys.readouterr()
+ assert ("poorly conditioned" not in out)
+ assert len(err) == 0
+
+def test_pivoted_greedy(capsys):
+ mu = [5.05, 7.1]
+ mu0 = [5., 7.]
+ solver = matrixRandom()
+ uh = solver.solve(mu)[0]
+ params = {"POD": True, "nTestPoints": 100, "greedyTol": 1e-4,
+ "collinearityTol": 1e8, "errorEstimatorKind": "DISCREPANCY",
+ "S": 5, "polybasis": "CHEBYSHEV",
+ "samplerPivot": QS([4.75, 5.25], "UNIFORM"),
+ "trainSetGenerator": QS([4.75, 5.25], "CHEBYSHEV"),
+ "MMarginal": 4, "SMarginal": 5, "polybasisMarginal": "MONOMIAL",
+ "samplerMarginal": QS([6.75, 7.25], "UNIFORM"),
+ "matchingWeight": 1.,
+ "cutOffTolerance": 1., "cutOffType": "MAGNITUDE"}
+ approx = RIGP(solver, mu0, [0], params, verbosity = 100)
+ approx.setupApprox()
+
+ uhP1 = approx.getApprox(mu)[0]
+ errP = approx.getErr(mu)[0]
+ errNP = approx.normErr(mu)[0]
+ myerrP = uhP1 - uh
+ assert np.allclose(np.abs(errP - myerrP), 0., rtol = 1e-3)
+ assert np.isclose(solver.norm(errP), errNP, rtol = 1e-3)
+ resP = approx.getRes(mu)[0]
+ resNP = approx.normRes(mu)
+ assert np.isclose(solver.norm(resP), resNP, rtol = 1e-3)
+ assert np.allclose(np.abs(resP - (solver.b(mu) - solver.A(mu).dot(uhP1))),
+ 0., rtol = 1e-3)
+ assert np.isclose(errNP / solver.norm(uh), 1.181958e-02, rtol = 1)
+
+ out, err = capsys.readouterr()
+ assert ("Overriding approx_state to True." in out)
+ assert len(err) == 0
diff --git a/tests/reduction_methods_multiD/rational_interpolant_2d.py b/tests/reduction_methods_multiD/rational_interpolant_2d.py
index 9bd84a3..067880c 100644
--- a/tests/reduction_methods_multiD/rational_interpolant_2d.py
+++ b/tests/reduction_methods_multiD/rational_interpolant_2d.py
@@ -1,79 +1,77 @@
# 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 matrix_random import matrixRandom
from rrompy.reduction_methods.standard import RationalInterpolant as RI
from rrompy.parameter.parameter_sampling import (QuadratureSampler as QS,
ManualSampler as MS)
def test_monomials(capsys):
mu = [5.05, 7.1]
mu0 = [5., 7.]
solver = matrixRandom()
uh = solver.solve(mu)[0]
params = {"POD": False, "M": 4, "N": 4, "S": 16, "robustTol": 1e-6,
"interpRcond": 1e-3, "polybasis": "MONOMIAL",
"sampler": QS([[4.9, 6.85], [5.1, 7.15]], "UNIFORM")}
approx = RI(solver, mu0, params, verbosity = 0)
approx.setupApprox()
uhP1 = approx.getApprox(mu)[0]
errP = approx.getErr(mu)[0]
errNP = approx.normErr(mu)[0]
myerrP = uhP1 - uh
assert np.allclose(np.abs(errP - myerrP), 0., rtol = 1e-3)
assert np.isclose(solver.norm(errP), errNP, rtol = 1e-3)
resP = approx.getRes(mu)[0]
resNP = approx.normRes(mu)
assert np.isclose(solver.norm(resP), resNP, rtol = 1e-3)
assert np.allclose(np.abs(resP - (solver.b(mu) - solver.A(mu).dot(uhP1))),
0., rtol = 1e-3)
assert np.isclose(errNP / solver.norm(uh), 5.2667e-05, rtol = 1)
out, err = capsys.readouterr()
- print(out)
assert ("poorly conditioned. Reducing N " in out)
assert len(err) == 0
def test_well_cond():
mu = [5.05, 7.1]
mu0 = [5., 7.]
solver = matrixRandom()
params = {"POD": True, "M": 3, "N": 3, "S": 16,
"interpRcond": 1e-10, "polybasis": "CHEBYSHEV",
"sampler": QS([[4.9, 6.85], [5.1, 7.15]], "UNIFORM")}
approx = RI(solver, mu0, params, verbosity = 0)
approx.setupApprox()
- print(approx.normErr(mu)[0] / approx.normHF(mu)[0])
assert np.isclose(approx.normErr(mu)[0] / approx.normHF(mu)[0],
5.98695e-05, rtol = 1e-1)
def test_hermite():
mu = [5.05, 7.1]
mu0 = [5., 7.]
solver = matrixRandom()
sampler0 = QS([[4.9, 6.85], [5.1, 7.15]], "UNIFORM")
params = {"POD": True, "M": 3, "N": 3, "S": 25, "polybasis": "CHEBYSHEV",
"sampler": MS([[4.9, 6.85], [5.1, 7.15]],
points = sampler0.generatePoints(9))}
approx = RI(solver, mu0, params, verbosity = 0)
approx.setupApprox()
assert np.isclose(approx.normErr(mu)[0] / approx.normHF(mu)[0],
0.000224, rtol = 5e-1)