diff --git a/rrompy/reduction_methods/pivoted/generic_pivoted_approximant.py b/rrompy/reduction_methods/pivoted/generic_pivoted_approximant.py
index 5d71dde..7fa339b 100644
--- a/rrompy/reduction_methods/pivoted/generic_pivoted_approximant.py
+++ b/rrompy/reduction_methods/pivoted/generic_pivoted_approximant.py
@@ -1,493 +1,473 @@
# 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.reduction_methods.base.generic_approximant import (
GenericApproximant)
from rrompy.utilities.poly_fitting.polynomial import polybases as ppb
from rrompy.utilities.poly_fitting.radial_basis import polybases as rbpb
from rrompy.utilities.poly_fitting.moving_least_squares import (
polybases as mlspb)
from rrompy.sampling import (SamplingEnginePivoted, SamplingEnginePivotedPOD,
SamplingEnginePivotedPODGlobal)
from rrompy.utilities.base.types import paramList, ListAny
from rrompy.utilities.base import verbosityManager as vbMng
from rrompy.utilities.numerical.degree import reduceDegreeN
from rrompy.utilities.exception_manager import (RROMPyException, RROMPyAssert,
RROMPyWarning)
__all__ = ['GenericPivotedApproximant', 'PODGlobal']
PODGlobal = 2
class GenericPivotedApproximant(GenericApproximant):
"""
ROM pivoted approximant (with pole matching) computation for parametric
problems (ABSTRACT).
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;
- 'cutOffTolerance': tolerance for ignoring parasitic poles;
defaults to np.inf;
- 'cutOffSharedRatio': required ratio of marginal points to share
resonance in cut off strategy; 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;
- 'polybasisMarginal': type of polynomial basis for marginal
interpolation; allowed values include 'MONOMIAL', 'CHEBYSHEV'
and 'LEGENDRE'; defaults to 'MONOMIAL';
- 'MMarginal': degree of marginal interpolant; defaults to 'AUTO',
i.e. maximum allowed;
- 'polydegreetypeMarginal': type of polynomial degree for marginal;
defaults to 'TOTAL';
- 'radialDirectionalWeightsMarginal': radial basis weights for
marginal interpolant; defaults to 1;
- - 'nNearestNeighborMarginal': number of marginal nearest neighbors
- considered if polybasisMarginal allows; defaults to -1;
- 'interpRcondMarginal': tolerance for marginal interpolation;
defaults to None.
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;
- 'cutOffTolerance': tolerance for ignoring parasitic poles;
- 'cutOffSharedRatio': required ratio of marginal points to share
resonance in cut off strategy;
- 'polybasisMarginal': type of polynomial basis for marginal
interpolation;
- '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;
- 'interpRcondMarginal': tolerance for marginal interpolation.
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.
cutOffTolerance: Tolerance for ignoring parasitic poles.
cutOffSharedRatio: Required ratio of marginal points to share resonance
in cut off strategy.
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.
polybasisMarginal: Type of polynomial basis for marginal interpolation.
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.
interpRcondMarginal: Tolerance for marginal interpolation.
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.
"""
def __init__(self, directionPivot:ListAny, *args, **kwargs):
self._preInit()
if len(directionPivot) > 1:
raise RROMPyException(("Exactly 1 pivot parameter allowed in pole "
"matching."))
from rrompy.parameter.parameter_sampling import QuadratureSampler as QS
QSBase = QS([[0], [1]], "UNIFORM")
self._addParametersToList(["matchingWeight", "cutOffTolerance",
"cutOffSharedRatio", "polybasisMarginal",
"MMarginal", "polydegreetypeMarginal",
"radialDirectionalWeightsMarginal",
- "nNearestNeighborMarginal",
- "interpRcondMarginal"],
- [1., np.inf, 1., "MONOMIAL", "AUTO", "TOTAL",
- [1.], -1, -1], ["samplerPivot", "SMarginal",
- "samplerMarginal"],
- [QSBase, [1], QSBase])
+ "interpRcondMarginal"], [1., np.inf, 1.,
+ "MONOMIAL", "AUTO", "TOTAL", [1.], -1],
+ ["samplerPivot", "SMarginal",
+ "samplerMarginal"], [QSBase, [1], QSBase])
del QS
self._directionPivot = directionPivot
super().__init__(*args, **kwargs)
self._postInit()
@property
def tModelType(self):
from .trained_model.trained_model_pivoted import TrainedModelPivoted
return TrainedModelPivoted
def setupSampling(self):
"""Setup sampling engine."""
RROMPyAssert(self._mode, message = "Cannot setup sampling engine.")
if not hasattr(self, "_POD") or self._POD is None: return
if self.POD:
if self.POD == PODGlobal:
SamplingEngine = SamplingEnginePivotedPODGlobal
else:
SamplingEngine = SamplingEnginePivotedPOD
else:
SamplingEngine = SamplingEnginePivoted
self.samplingEngine = SamplingEngine(self.HFEngine,
self.directionPivot,
sample_state = self.approx_state,
verbosity = self.verbosity)
def initializeModelData(self, datadict):
if "directionPivot" in datadict.keys():
from .trained_model.trained_model_pivoted_data import (
TrainedModelPivotedData)
return (TrainedModelPivotedData(datadict["mu0"],
datadict.pop("projMat"),
datadict["scaleFactor"],
datadict.pop("rescalingExp"),
datadict["directionPivot"]),
["mu0", "scaleFactor", "directionPivot", "mus"])
else:
return super().initializeModelData(datadict)
@property
def npar(self):
"""Number of parameters."""
if hasattr(self, "_temporaryPivot"): return self.nparPivot
return super().npar
@property
def mus(self):
"""Value of mus. Its assignment may reset snapshots."""
return self._mus
@mus.setter
def mus(self, mus):
musOld = copy(self.mus) if hasattr(self, '_mus') else None
if (musOld is None or len(mus) != len(musOld) or not mus == musOld):
self.resetSamples()
self._mus = mus
@property
def matchingWeight(self):
"""Value of matchingWeight."""
return self._matchingWeight
@matchingWeight.setter
def matchingWeight(self, matchingWeight):
self._matchingWeight = matchingWeight
self._approxParameters["matchingWeight"] = self.matchingWeight
@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 cutOffSharedRatio(self):
"""Value of cutOffSharedRatio."""
return self._cutOffSharedRatio
@cutOffSharedRatio.setter
def cutOffSharedRatio(self, cutOffSharedRatio):
if cutOffSharedRatio > 1.:
RROMPyWarning("Cut off shared ratio too large. Clipping to 1.")
cutOffSharedRatio = 1.
elif cutOffSharedRatio < 0.:
RROMPyWarning("Cut off shared ratio too small. Clipping to 0.")
cutOffSharedRatio = 0.
self._cutOffSharedRatio = cutOffSharedRatio
self._approxParameters["cutOffSharedRatio"] = self.cutOffSharedRatio
@property
def SMarginal(self):
"""Value of SMarginal."""
return self._SMarginal
@SMarginal.setter
def SMarginal(self, SMarginal):
if SMarginal <= 0:
raise RROMPyException("SMarginal must be positive.")
if hasattr(self, "_SMarginal") and self._SMarginal is not None:
Sold = self.SMarginal
else: Sold = -1
self._SMarginal = SMarginal
self._approxParameters["SMarginal"] = self.SMarginal
if Sold != self.SMarginal: self.resetSamples()
@property
def polybasisMarginal(self):
"""Value of polybasisMarginal."""
return self._polybasisMarginal
@polybasisMarginal.setter
def polybasisMarginal(self, polybasisMarginal):
try:
polybasisMarginal = polybasisMarginal.upper().strip().replace(" ",
"")
if polybasisMarginal not in ppb + rbpb + mlspb:
raise RROMPyException(
"Prescribed marginal polybasis not recognized.")
self._polybasisMarginal = polybasisMarginal
except:
RROMPyWarning(("Prescribed marginal polybasis not recognized. "
"Overriding to 'MONOMIAL'."))
self._polybasisMarginal = "MONOMIAL"
self._approxParameters["polybasisMarginal"] = self.polybasisMarginal
@property
def MMarginal(self):
"""Value of MMarginal."""
return self._MMarginal
@MMarginal.setter
def MMarginal(self, MMarginal):
if isinstance(MMarginal, str):
MMarginal = MMarginal.strip().replace(" ","")
if "-" not in MMarginal: MMarginal = MMarginal + "-0"
self._MMarginal_isauto = True
self._MMarginal_shift = int(MMarginal.split("-")[-1])
MMarginal = 0
if MMarginal < 0:
raise RROMPyException("MMarginal must be non-negative.")
self._MMarginal = MMarginal
self._approxParameters["MMarginal"] = self.MMarginal
def _setMMarginalAuto(self):
self.MMarginal = max(0, reduceDegreeN(
len(self.musMarginal), len(self.musMarginal),
self.nparMarginal, self.polydegreetypeMarginal
) - self._MMarginal_shift)
vbMng(self, "MAIN", ("Automatically setting MMarginal to "
"{}.").format(self.MMarginal), 25)
@property
def polydegreetypeMarginal(self):
"""Value of polydegreetypeMarginal."""
return self._polydegreetypeMarginal
@polydegreetypeMarginal.setter
def polydegreetypeMarginal(self, polydegreetypeM):
try:
polydegreetypeM = polydegreetypeM.upper().strip().replace(" ","")
if polydegreetypeM not in ["TOTAL", "FULL"]:
raise RROMPyException(("Prescribed polydegreetypeMarginal not "
"recognized."))
self._polydegreetypeMarginal = polydegreetypeM
except:
RROMPyWarning(("Prescribed polydegreetypeMarginal not recognized. "
"Overriding to 'TOTAL'."))
self._polydegreetypeMarginal = "TOTAL"
self._approxParameters["polydegreetypeMarginal"] = (
self.polydegreetypeMarginal)
@property
def radialDirectionalWeightsMarginal(self):
"""Value of radialDirectionalWeightsMarginal."""
return self._radialDirectionalWeightsMarginal
@radialDirectionalWeightsMarginal.setter
def radialDirectionalWeightsMarginal(self, radialDirWeightsMarginal):
if not hasattr(radialDirWeightsMarginal, "__len__"):
radialDirWeightsMarginal = [radialDirWeightsMarginal]
self._radialDirectionalWeightsMarginal = radialDirWeightsMarginal
self._approxParameters["radialDirectionalWeightsMarginal"] = (
self.radialDirectionalWeightsMarginal)
- @property
- def nNearestNeighborMarginal(self):
- """Value of nNearestNeighborMarginal."""
- return self._nNearestNeighborMarginal
- @nNearestNeighborMarginal.setter
- def nNearestNeighborMarginal(self, nNearestNeighborMarginal):
- self._nNearestNeighborMarginal = nNearestNeighborMarginal
- self._approxParameters["nNearestNeighborMarginal"] = (
- self.nNearestNeighborMarginal)
-
@property
def interpRcondMarginal(self):
"""Value of interpRcondMarginal."""
return self._interpRcondMarginal
@interpRcondMarginal.setter
def interpRcondMarginal(self, interpRcondMarginal):
self._interpRcondMarginal = interpRcondMarginal
self._approxParameters["interpRcondMarginal"] = (
self.interpRcondMarginal)
@property
def directionPivot(self):
"""Value of directionPivot. Its assignment may reset snapshots."""
return self._directionPivot
@directionPivot.setter
def directionPivot(self, directionPivot):
if hasattr(self, '_directionPivot'):
directionPivotOld = copy(self.directionPivot)
else:
directionPivotOld = None
if (directionPivotOld is None
or len(directionPivot) != len(directionPivotOld)
or not directionPivot == directionPivotOld):
self.resetSamples()
self._directionPivot = directionPivot
@property
def directionMarginal(self):
return [x for x in range(self.HFEngine.npar) \
if x not in self.directionPivot]
@property
def nparPivot(self):
return len(self.directionPivot)
@property
def nparMarginal(self):
return self.npar - self.nparPivot
@property
def rescalingExpPivot(self):
return [self.HFEngine.rescalingExp[x] for x in self.directionPivot]
@property
def rescalingExpMarginal(self):
return [self.HFEngine.rescalingExp[x] for x in self.directionMarginal]
@property
def muBounds(self):
"""Value of muBounds."""
return self.samplerPivot.lims
@property
def muBoundsMarginal(self):
"""Value of muBoundsMarginal."""
return self.samplerMarginal.lims
@property
def samplerPivot(self):
"""Value of samplerPivot."""
return self._samplerPivot
@samplerPivot.setter
def samplerPivot(self, samplerPivot):
if 'generatePoints' not in dir(samplerPivot):
raise RROMPyException("Pivot sampler type not recognized.")
if hasattr(self, '_samplerPivot') and self._samplerPivot is not None:
samplerOld = self.samplerPivot
self._samplerPivot = samplerPivot
self._approxParameters["samplerPivot"] = self.samplerPivot.__str__()
if not 'samplerOld' in locals() or samplerOld != self.samplerPivot:
self.resetSamples()
@property
def samplerMarginal(self):
"""Value of samplerMarginal."""
return self._samplerMarginal
@samplerMarginal.setter
def samplerMarginal(self, samplerMarginal):
if 'generatePoints' not in dir(samplerMarginal):
raise RROMPyException("Marginal sampler type not recognized.")
if (hasattr(self, '_samplerMarginal')
and self._samplerMarginal is not None):
samplerOld = self.samplerMarginal
self._samplerMarginal = samplerMarginal
self._approxParameters["samplerMarginal"] = (
self.samplerMarginal.__str__())
if not 'samplerOld' in locals() or samplerOld != self.samplerMarginal:
self.resetSamples()
def setSamples(self, samplingEngine):
"""Copy samplingEngine and samples."""
self.mus = copy(samplingEngine.mus[0])
for sEj in samplingEngine.mus[1:]:
self.mus.append(sEj)
super().setSamples(samplingEngine)
def _finalizeMarginalization(self):
vbMng(self, "INIT", "Recompressing by cut off.", 10)
msg = self.trainedModel.recompressByCutOff(self.cutOffTolerance,
self.cutOffSharedRatio,
self.samplerPivot.normalFoci(),
self.samplerPivot.groundPotential())
vbMng(self, "DEL", "Done recompressing." + msg, 10)
interpPars = [self.verbosity >= 5,
self.polydegreetypeMarginal == "TOTAL", {}]
- if self.polybasisMarginal not in ppb:
- interpPars[-1]["nNearestNeighbor"] = self.nNearestNeighborMarginal
if self.polybasisMarginal in ppb + rbpb:
interpPars += [{"rcond": self.interpRcondMarginal}]
self.trainedModel.setupMarginalInterp(self, interpPars,
hasattr(self, "_MMarginal_isauto"),
self.radialDirectionalWeightsMarginal,
hasattr(self, "_reduceDegreeNNoWarn"))
self.trainedModel.data.approxParameters = copy(self.approxParameters)
def computeScaleFactor(self):
"""Compute parameter rescaling factor."""
RROMPyAssert(self._mode, message = "Cannot compute rescaling factor.")
self.scaleFactorPivot = .5 * np.abs(
self.muBounds[0] ** self.rescalingExpPivot
- self.muBounds[1] ** self.rescalingExpPivot)
self.scaleFactorMarginal = .5 * np.abs(
self.muBoundsMarginal[0] ** self.rescalingExpMarginal
- self.muBoundsMarginal[1] ** self.rescalingExpMarginal)
self.scaleFactor = np.empty(self.npar)
self.scaleFactor[self.directionPivot] = self.scaleFactorPivot
self.scaleFactor[self.directionMarginal] = self.scaleFactorMarginal
def normApprox(self, mu:paramList) -> float:
_PODOld = self.POD
self._POD = self.POD == PODGlobal
result = super().normApprox(mu)
self._POD = _PODOld
return result
diff --git a/rrompy/reduction_methods/pivoted/greedy/generic_pivoted_greedy_approximant.py b/rrompy/reduction_methods/pivoted/greedy/generic_pivoted_greedy_approximant.py
index debfe17..17f417f 100644
--- a/rrompy/reduction_methods/pivoted/greedy/generic_pivoted_greedy_approximant.py
+++ b/rrompy/reduction_methods/pivoted/greedy/generic_pivoted_greedy_approximant.py
@@ -1,724 +1,718 @@
# 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 abc import abstractmethod
from copy import deepcopy as copy
import numpy as np
from matplotlib import pyplot as plt
from rrompy.reduction_methods.pivoted.generic_pivoted_approximant import (
GenericPivotedApproximant, PODGlobal)
from rrompy.utilities.base.types import (Np1D, Tuple, List, paramVal,
paramList, ListAny)
from rrompy.utilities.base import verbosityManager as vbMng
from rrompy.utilities.numerical import dot
from rrompy.utilities.numerical.point_matching import (pointMatching,
chordalMetricAdjusted, potential)
from rrompy.utilities.exception_manager import (RROMPyException, RROMPyAssert,
RROMPyWarning)
from rrompy.parameter import checkParameterList, emptyParameterList
__all__ = ['GenericPivotedGreedyApproximant']
class GenericPivotedGreedyApproximant(GenericPivotedApproximant):
"""
ROM pivoted greedy interpolant computation for parametric problems
(ABSTRACT).
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;
- 'cutOffTolerance': tolerance for ignoring parasitic poles;
defaults to np.inf;
- 'cutOffSharedRatio': required ratio of marginal points to share
resonance in cut off strategy; 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;
- 'samplerMarginalGrid': 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';
- 'polybasisMarginal': type of polynomial basis for marginal
interpolation; allowed values include 'MONOMIAL', 'CHEBYSHEV'
and 'LEGENDRE'; defaults to 'MONOMIAL';
- 'MMarginal': degree of marginal interpolant; 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;
- 'polydegreetypeMarginal': type of polynomial degree for marginal;
defaults to 'TOTAL';
- 'radialDirectionalWeightsMarginal': radial basis weights for
marginal interpolant; defaults to 1;
- - 'nNearestNeighborMarginal': number of marginal nearest neighbors
- considered if polybasisMarginal allows; defaults to -1;
- 'interpRcondMarginal': tolerance for marginal interpolation;
defaults to None.
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;
- 'cutOffTolerance': tolerance for ignoring parasitic poles;
- 'cutOffSharedRatio': required ratio of marginal points to share
resonance in cut off strategy;
- 'matchingWeightError': weight for pole matching optimization in
error estimation;
- 'cutOffToleranceError': tolerance for ignoring parasitic poles
in error estimation;
- 'errorEstimatorKindMarginal': kind of marginal error estimator;
- 'polybasisMarginal': type of polynomial basis for marginal
interpolation;
- 'MMarginal': degree of marginal interpolant;
- 'greedyTolMarginal': uniform error tolerance for marginal greedy
algorithm;
- 'maxIterMarginal': maximum number of marginal greedy steps;
- 'polydegreetypeMarginal': type of polynomial degree for marginal;
- 'radialDirectionalWeightsMarginal': radial basis weights for
marginal interpolant;
- - 'nNearestNeighborMarginal': number of marginal nearest neighbors
- considered if polybasisMarginal allows;
- 'interpRcondMarginal': tolerance for marginal interpolation.
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;
- 'samplerMarginalGrid': 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.
cutOffTolerance: Tolerance for ignoring parasitic poles.
cutOffSharedRatio: Required ratio of marginal points to share resonance
in cut off strategy.
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.
samplerMarginalGrid: Marginal sample point generator via sparse grid.
errorEstimatorKindMarginal: Kind of marginal error estimator.
polybasisMarginal: Type of polynomial basis for marginal interpolation.
MMarginal: Degree of marginal interpolant.
greedyTolMarginal: Uniform error tolerance for marginal greedy
algorithm.
maxIterMarginal: Maximum number of marginal greedy steps.
polydegreetypeMarginal: Type of polynomial degree for marginal.
radialDirectionalWeightsMarginal: Radial basis weights for marginal
interpolant.
- nNearestNeighborMarginal: Number of marginal nearest neighbors
- considered if polybasisMarginal allows.
interpRcondMarginal: Tolerance for marginal interpolation.
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.
"""
_allowedEstimatorKindsMarginal = ["LEAVE_ONE_OUT", "LOOK_AHEAD",
"LOOK_AHEAD_RECOVER"]
def __init__(self, *args, **kwargs):
self._preInit()
from rrompy.parameter import localSparseGrid as SG
SGBase = SG([[0.], [1.]], "UNIFORM")
self._addParametersToList(["matchingWeightError",
"cutOffToleranceError",
"errorEstimatorKindMarginal",
"greedyTolMarginal",
"maxIterMarginal"],
[0., "AUTO", "LEAVE_ONE_OUT", 1e-1, 1e2],
["samplerMarginalGrid"], [SGBase],
toBeExcluded = ["samplerMarginal"])
super().__init__(*args, **kwargs)
self._postInit()
@property
def muBoundsMarginal(self):
"""Value of muBoundsMarginal."""
return self.samplerMarginalGrid.lims
@property
def scaleFactorDer(self):
"""Value of scaleFactorDer."""
if self._scaleFactorDer == "NONE": return 1.
if self._scaleFactorDer == "AUTO": return self._scaleFactorOldPivot
return self._scaleFactorDer
@scaleFactorDer.setter
def scaleFactorDer(self, scaleFactorDer):
if hasattr(self, "_scaleFactorDer"):
scaleFactorDerold = self.scaleFactorDer
else: scaleFactorDerold = -1
if isinstance(scaleFactorDer, (str,)):
scaleFactorDer = scaleFactorDer.upper()
self._scaleFactorDer = scaleFactorDer
self._approxParameters["scaleFactorDer"] = self._scaleFactorDer
if scaleFactorDerold != self._scaleFactorDer: self.resetSamples()
@property
def samplerMarginalGrid(self):
"""Value of samplerMarginalGrid."""
return self._samplerMarginalGrid
@samplerMarginalGrid.setter
def samplerMarginalGrid(self, samplerMarginalGrid):
if 'refine' not in dir(samplerMarginalGrid):
raise RROMPyException("Marginal sampler type not recognized.")
if (hasattr(self, '_samplerMarginalGrid')
and self._samplerMarginalGrid is not None):
samplerOld = self.samplerMarginalGrid
self._samplerMarginalGrid = samplerMarginalGrid
self._approxParameters["samplerMarginalGrid"] = (
self.samplerMarginalGrid.__str__())
if (not 'samplerOld' in locals()
or samplerOld != self.samplerMarginalGrid):
self.resetSamples()
@property
def errorEstimatorKindMarginal(self):
"""Value of errorEstimatorKindMarginal."""
return self._errorEstimatorKindMarginal
@errorEstimatorKindMarginal.setter
def errorEstimatorKindMarginal(self, errorEstimatorKindMarginal):
errorEstimatorKindMarginal = errorEstimatorKindMarginal.upper()
if errorEstimatorKindMarginal not in (
self._allowedEstimatorKindsMarginal):
RROMPyWarning(("Marginal error estimator kind not recognized. "
"Overriding to 'LEAVE_ONE_OUT'."))
errorEstimatorKindMarginal = "LEAVE_ONE_OUT"
self._errorEstimatorKindMarginal = errorEstimatorKindMarginal
self._approxParameters["errorEstimatorKindMarginal"] = (
self.errorEstimatorKindMarginal)
@property
def matchingWeightError(self):
"""Value of matchingWeightError."""
return self._matchingWeightError
@matchingWeightError.setter
def matchingWeightError(self, matchingWeightError):
self._matchingWeightError = matchingWeightError
self._approxParameters["matchingWeightError"] = (
self.matchingWeightError)
@property
def cutOffToleranceError(self):
"""Value of cutOffToleranceError."""
return self._cutOffToleranceError
@cutOffToleranceError.setter
def cutOffToleranceError(self, cutOffToleranceError):
if isinstance(cutOffToleranceError, (str,)):
cutOffToleranceError = cutOffToleranceError.upper()\
.strip().replace(" ","")
if cutOffToleranceError != "AUTO":
RROMPyWarning(("String value of cutOffToleranceError not "
"recognized. Overriding to 'AUTO'."))
cutOffToleranceError == "AUTO"
self._cutOffToleranceError = cutOffToleranceError
self._approxParameters["cutOffToleranceError"] = (
self.cutOffToleranceError)
@property
def greedyTolMarginal(self):
"""Value of greedyTolMarginal."""
return self._greedyTolMarginal
@greedyTolMarginal.setter
def greedyTolMarginal(self, greedyTolMarginal):
if greedyTolMarginal < 0:
raise RROMPyException("greedyTolMarginal must be non-negative.")
if (hasattr(self, "_greedyTolMarginal")
and self.greedyTolMarginal is not None):
greedyTolMarginalold = self.greedyTolMarginal
else:
greedyTolMarginalold = -1
self._greedyTolMarginal = greedyTolMarginal
self._approxParameters["greedyTolMarginal"] = self.greedyTolMarginal
if greedyTolMarginalold != self.greedyTolMarginal:
self.resetSamples()
@property
def maxIterMarginal(self):
"""Value of maxIterMarginal."""
return self._maxIterMarginal
@maxIterMarginal.setter
def maxIterMarginal(self, maxIterMarginal):
if maxIterMarginal <= 0:
raise RROMPyException("maxIterMarginal must be positive.")
if (hasattr(self, "_maxIterMarginal")
and self.maxIterMarginal is not None):
maxIterMarginalold = self.maxIterMarginal
else:
maxIterMarginalold = -1
self._maxIterMarginal = maxIterMarginal
self._approxParameters["maxIterMarginal"] = self.maxIterMarginal
if maxIterMarginalold != self.maxIterMarginal:
self.resetSamples()
def resetSamples(self):
"""Reset samples."""
super().resetSamples()
if not hasattr(self, "_temporaryPivot"):
self._mus = emptyParameterList()
self.musMarginal = emptyParameterList()
if hasattr(self, "samplerMarginalGrid"):
self.samplerMarginalGrid.reset()
if hasattr(self, "samplingEngine") and self.samplingEngine is not None:
self.samplingEngine.resetHistory()
def getErrorEstimatorMarginalLookAhead(self) -> Np1D:
if not hasattr(self.trainedModel, "_musMExcl"):
err = np.zeros(0)
err[:] = np.inf
self._musMarginalTestIdxs = np.zeros(0, dtype = int)
return err
err = np.zeros(len(self.trainedModel._musMExcl))
self._musMarginalTestIdxs = np.array(self.trainedModel._idxExcl,
dtype = int)
self.verbosity -= 35
self.trainedModel.verbosity -= 35
if self.cutOffToleranceError == "AUTO":
cutOffTolErr = self.cutOffTolerance
else:
cutOffTolErr = self.cutOffToleranceError
foci = self.samplerPivot.normalFoci()
ground = self.samplerPivot.groundPotential()
for j, (muTest, HITest) in enumerate(zip(self.trainedModel._musMExcl,
self.trainedModel._HIsExcl)):
polesEx = HITest.poles
idxExEff = np.where(potential(polesEx, foci) - ground
<= cutOffTolErr * ground)[0]
polesEx = polesEx[idxExEff]
if self.matchingWeightError != 0:
resEx = HITest.coeffs[idxExEff]
else:
resEx = None
if len(polesEx) == 0: continue
polesAp = self.trainedModel.interpolateMarginalPoles(muTest)[...,
0]
idxApEff = np.where(potential(polesAp, foci) - ground
<= cutOffTolErr * ground)[0]
polesAp = polesAp[idxApEff]
if self.matchingWeightError != 0:
resAp = self.trainedModel.interpolateMarginalCoeffs(
muTest)[idxApEff, :, 0]
if self.POD != PODGlobal:
resEx = self.trainedModel.data.projMat[:,
: resEx.shape[1]].dot(resEx.T)
resAp = self.trainedModel.data.projMat[:,
: resAp.shape[1]].dot(resAp.T)
else:
resAp = None
dist = chordalMetricAdjusted(polesEx, polesAp,
self.matchingWeightError, resEx,
resAp, self.HFEngine, False)
pmR, pmC = pointMatching(dist)
err[j] = np.mean(dist[pmR, pmC])
self.verbosity += 35
self.trainedModel.verbosity += 35
return err
def getErrorEstimatorMarginalLeaveOneOut(self) -> Np1D:
err = np.zeros(len(self.trainedModel.data.musMarginal))
self._musMarginalTestIdxs = np.arange(len(err))
if len(err) <= 1:
err[:] = np.inf
return err
_tMdataFull = copy(self.trainedModel.data)
if self.cutOffToleranceError == "AUTO":
cutOffTolErr = self.cutOffTolerance
else:
cutOffTolErr = self.cutOffToleranceError
if not hasattr(self, "_MMarginal_isauto"):
if not hasattr(self, "_MMarginalOriginal"):
self._MMarginalOriginal = self.MMarginal
self.MMarginal = self._MMarginalOriginal
_musMExcl = None
self.verbosity -= 35
self.trainedModel.verbosity -= 35
foci = self.samplerPivot.normalFoci()
ground = self.samplerPivot.groundPotential()
for j in range(len(err)):
jEff = j - (j > 0)
muTest = self.trainedModel.data.musMarginal[jEff]
polesEx = self.trainedModel.data.HIs[jEff].poles
idxExEff = np.where(potential(polesEx, foci) - ground
<= cutOffTolErr * ground)[0]
polesEx = polesEx[idxExEff]
if self.matchingWeightError != 0:
resEx = self.trainedModel.data.HIs[jEff].coeffs[idxExEff]
else:
resEx = None
if j > 0: self.musMarginal.insert(_musMExcl, j - 1)
_musMExcl = self.musMarginal[j]
self.musMarginal.pop(j)
if len(polesEx) == 0: continue
self.trainedModel.updateEffectiveSamples(self.HFEngine, [j],
self.matchingWeight,
self.POD == PODGlobal, False)
self._reduceDegreeNNoWarn = 1
self._finalizeMarginalization()
polesAp = self.trainedModel.interpolateMarginalPoles(muTest)[...,
0]
idxApEff = np.where(potential(polesAp, foci) - ground
<= cutOffTolErr * ground)[0]
polesAp = polesAp[idxApEff]
if self.matchingWeightError != 0:
resAp = self.trainedModel.interpolateMarginalCoeffs(
muTest)[idxApEff, :, 0]
if self.POD != PODGlobal:
resEx = self.trainedModel.data.projMat[:,
: resEx.shape[1]].dot(resEx.T)
resAp = self.trainedModel.data.projMat[:,
: resAp.shape[1]].dot(resAp.T)
else:
resAp = None
dist = chordalMetricAdjusted(polesEx, polesAp,
self.matchingWeightError, resEx,
resAp, self.HFEngine, False)
pmR, pmC = pointMatching(dist)
err[j] = np.mean(dist[pmR, pmC])
self.trainedModel.updateEffectiveSamples(self.HFEngine, [],
self.matchingWeight,
self.POD == PODGlobal, False)
if not hasattr(self, "_MMarginal_isauto"):
self.MMarginal = self._MMarginalOriginal
self.musMarginal.append(_musMExcl)
self.verbosity += 35
self.trainedModel.verbosity += 35
self.trainedModel.data = _tMdataFull
if hasattr(self, "_reduceDegreeNNoWarn"): del self._reduceDegreeNNoWarn
return err
def errorEstimatorMarginal(self, return_max : bool = False) -> Np1D:
vbMng(self.trainedModel, "INIT",
"Evaluating error estimator at mu = {}.".format(
self.trainedModel.data.musMarginal), 10)
if self.errorEstimatorKindMarginal == "LEAVE_ONE_OUT":
err = self.getErrorEstimatorMarginalLeaveOneOut()
else:#if self.errorEstimatorKindMarginal[: 10] == "LOOK_AHEAD":
err = self.getErrorEstimatorMarginalLookAhead()
vbMng(self.trainedModel, "DEL", "Done evaluating error estimator", 10)
if not return_max: return err
idxMaxEst = np.where(err > self.greedyTolMarginal)[0]
return err, idxMaxEst, err[idxMaxEst]
def plotEstimatorMarginal(self, est:Np1D, idxMax:List[int],
estMax:List[float]):
if not (np.any(np.isnan(est)) or np.any(np.isinf(est))):
fig = plt.figure(figsize = plt.figaspect(1. / self.nparMarginal))
for jpar in range(self.nparMarginal):
ax = fig.add_subplot(1, self.nparMarginal, 1 + jpar)
if self.errorEstimatorKindMarginal == "LEAVE_ONE_OUT":
musre = copy(self.trainedModel.data.musMarginal.re.data)
else:#if self.errorEstimatorKindMarginal[: 10] == "LOOK_AHEAD":
musre = np.real(self.trainedModel._musMExcl)
if len(idxMax) > 0 and estMax is not None:
maxrej = musre[idxMax, jpar]
errCP = copy(est)
idx = np.delete(np.arange(self.nparMarginal), jpar)
while len(musre) > 0:
if self.nparMarginal == 1:
currIdx = np.arange(len(musre))
else:
currIdx = np.where(np.isclose(np.sum(
np.abs(musre[:, idx] - musre[0, idx]), 1), 0.))[0]
currIdxSorted = currIdx[np.argsort(musre[currIdx, jpar])]
ax.semilogy(musre[currIdxSorted, jpar],
errCP[currIdxSorted], 'k.-', linewidth = 1)
musre = np.delete(musre, currIdx, 0)
errCP = np.delete(errCP, currIdx)
ax.semilogy(self.musMarginal.re(jpar),
(self.greedyTolMarginal,) * len(self.musMarginal),
'*m')
if len(idxMax) > 0 and estMax is not None:
ax.semilogy(maxrej, estMax, 'xr')
ax.grid()
plt.tight_layout()
plt.show()
def _addMarginalSample(self, mus:paramList):
mus = checkParameterList(mus, self.nparMarginal)[0]
if len(mus) == 0: return
nmus = len(mus)
vbMng(self, "MAIN",
("Adding marginal sample point{} no. {}{} at {} to training "
"set.").format("s" * (nmus > 1), len(self.musMarginal) + 1,
"--{}".format(len(self.musMarginal) + nmus) * (nmus > 1),
mus), 3)
self.musMarginal.append(mus)
self.setupApproxPivoted(mus)
vbMng(self, "INIT", "Matching poles.", 10)
self.trainedModel.initializeFromRational(self.HFEngine,
self.matchingWeight,
self.POD == PODGlobal, False)
vbMng(self, "DEL", "Done matching poles.", 10)
if (self.errorEstimatorKindMarginal[: 10] == "LOOK_AHEAD"
and not self.firstGreedyIterM):
ubRange = len(self.trainedModel.data.musMarginal)
if hasattr(self.trainedModel, "_idxExcl"):
shRange = len(self.trainedModel._musMExcl)
else:
shRange = 0
testIdxs = list(range(ubRange + shRange - len(mus),
ubRange + shRange))
for j in testIdxs[::-1]:
self.musMarginal.pop(j - shRange)
if hasattr(self.trainedModel, "_idxExcl"):
testIdxs = self.trainedModel._idxExcl + testIdxs
self.trainedModel.updateEffectiveSamples(self.HFEngine, testIdxs,
self.matchingWeight,
self.POD == PODGlobal, False)
self._finalizeMarginalization()
self._SMarginal = len(self.musMarginal)
self._approxParameters["SMarginal"] = self.SMarginal
self.trainedModel.data.approxParameters["SMarginal"] = self.SMarginal
def greedyNextSampleMarginal(self, muidx:List[int],
plotEst : str = "NONE") \
-> Tuple[Np1D, List[int], float, paramVal]:
RROMPyAssert(self._mode, message = "Cannot add greedy sample.")
if (self.errorEstimatorKindMarginal[: 10] == "LOOK_AHEAD"
and not self.firstGreedyIterM):
if not hasattr(self.trainedModel, "_idxExcl"):
raise RROMPyException(("Sample index to be added not present "
"in trained model."))
testIdxs = copy(self.trainedModel._idxExcl)
skippedIdx = 0
for cj, j in enumerate(self.trainedModel._idxExcl):
if j in muidx:
testIdxs.pop(skippedIdx)
self.musMarginal.insert(self.trainedModel._musMExcl[cj],
j - skippedIdx)
else:
skippedIdx += 1
if len(self.trainedModel._idxExcl) < (len(muidx)
+ len(testIdxs)):
raise RROMPyException(("Sample index to be added not present "
"in trained model."))
self.trainedModel.updateEffectiveSamples(self.HFEngine, testIdxs,
self.matchingWeight,
self.POD == PODGlobal, False)
self._SMarginal = len(self.musMarginal)
self._approxParameters["SMarginal"] = self.SMarginal
self.trainedModel.data.approxParameters["SMarginal"] = (
self.SMarginal)
self.firstGreedyIterM = False
idxAdded = self.samplerMarginalGrid.refine(muidx)
self._addMarginalSample(self.samplerMarginalGrid.points[idxAdded])
errorEstTest, muidx, maxErrorEst = self.errorEstimatorMarginal(True)
if plotEst == "ALL":
self.plotEstimatorMarginal(errorEstTest, muidx, maxErrorEst)
return (errorEstTest, self._musMarginalTestIdxs[muidx], maxErrorEst,
self.samplerMarginalGrid.points[muidx])
def _preliminaryTrainingMarginal(self):
"""Initialize starting snapshots of solution map."""
RROMPyAssert(self._mode, message = "Cannot start greedy algorithm.")
if np.sum(self.samplingEngine.nsamples) > 0: return
self.resetSamples()
idx = [0]
while self.samplerMarginalGrid.npoints < self.SMarginal:
idx = self.samplerMarginalGrid.refine(idx)
self._addMarginalSample(self.samplerMarginalGrid.points)
def _finalizeSnapshots(self):
self.samplingEngine = self._samplingEngineOld
for muM, sEN in zip(self.musMargLoc, self.samplingEngs):
self.samplingEngine.samples += [sEN.samples]
self.samplingEngine.nsamples += [sEN.nsamples]
self.samplingEngine.mus += [sEN.mus]
self.samplingEngine.musMarginal.append(muM)
self.samplingEngine._derIdxs += [[(0,) * self.npar]
for _ in range(sEN.nsamples)]
if self.POD:
self.samplingEngine.RPOD += [sEN.RPOD]
self.samplingEngine.samples_full += [copy(sEN.samples_full)]
if self.POD == PODGlobal:
self.samplingEngine.coalesceSamples(self.interpRcondMarginal)
else:
self.samplingEngine.coalesceSamples()
def _preSetupApproxPivoted(self, mus:paramList) -> Tuple[ListAny, ListAny]:
self.computeScaleFactor()
if self.trainedModel is None:
self.trainedModel = self.tModelType()
self.trainedModel.verbosity = self.verbosity
self.trainedModel.timestamp = self.timestamp
datadict = {"mu0": self.mu0, "projMat": np.zeros((0, 0)),
"scaleFactor": self.scaleFactor,
"rescalingExp": self.HFEngine.rescalingExp,
"directionPivot": self.directionPivot}
self.trainedModel.data = self.initializeModelData(datadict)[0]
self.trainedModel.data.Qs, self.trainedModel.data.Ps = [], []
self._trainedModelOld = copy(self.trainedModel)
self._scaleFactorOldPivot = copy(self.scaleFactor)
self.scaleFactor = self.scaleFactorPivot
self._temporaryPivot = 1
self._samplingEngineOld = copy(self.samplingEngine)
self.musMargLoc, self.samplingEngs = [], [None] * len(mus)
Qs, Ps = [None] * len(mus), [None] * len(mus)
self.verbosity -= 15
return Qs, Ps
def _postSetupApproxPivoted(self, mus:paramList, Qs:ListAny, Ps:ListAny):
self.scaleFactor = self._scaleFactorOldPivot
del self._scaleFactorOldPivot, self._temporaryPivot
self._finalizeSnapshots()
del self._samplingEngineOld, self.musMargLoc, self.samplingEngs
self._mus = self.samplingEngine.musCoalesced
self.trainedModel = self._trainedModelOld
del self._trainedModelOld
self.trainedModel.data.mus = copy(self.mus)
self.trainedModel.data.musMarginal = copy(self.musMarginal)
padRight = (self.samplingEngine.nsamplesTot
- self.trainedModel.data.projMat.shape[1])
nmusOld = len(self.trainedModel.data.Ps)
for j in range(nmusOld):
self.trainedModel.data.Ps[j].pad(0, padRight)
self.trainedModel.data.HIs[j].pad(0, padRight)
if hasattr(self.trainedModel, "_PsExcl"):
nmusOldExcl = len(self.trainedModel._PsExcl)
for j in range(nmusOldExcl):
self.trainedModel._PsExcl[j].pad(0, padRight)
self.trainedModel._HIsExcl[j].pad(0, padRight)
nmusOld += nmusOldExcl
padLeft = self.trainedModel.data.projMat.shape[1]
for j in range(len(mus)):
nsj = self.samplingEngine.nsamples[nmusOld + j]
if self.POD == PODGlobal:
rRightj = self.samplingEngine.RPODCPart[:,
padLeft : padLeft + nsj]
Ps[j].postmultiplyTensorize(rRightj.T)
else:
padRight -= nsj
Ps[j].pad(padLeft, padRight)
padLeft += nsj
pMat = self.samplingEngine.samplesCoalesced.data
pMatEff = dot(self.HFEngine.C, pMat) if self.approx_state else pMat
self.trainedModel.data.projMat = pMatEff
self.trainedModel.data.Qs += Qs
self.trainedModel.data.Ps += Ps
self.trainedModel.data.approxParameters = copy(self.approxParameters)
self.verbosity += 15
@abstractmethod
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)
Qs, Ps = self._preSetupApproxPivoted()
pass
self._postSetupApproxPivoted(mus, Qs, Ps)
vbMng(self, "DEL", "Done setting up pivoted approximant.", 10)
return 0
def setupApprox(self, plotEst : str = "NONE") -> int:
"""Compute greedy snapshots of solution map."""
if self.checkComputedApprox(): return -1
RROMPyAssert(self._mode, message = "Cannot start greedy algorithm.")
vbMng(self, "INIT", "Starting computation of snapshots.", 3)
max2ErrorEst, self.firstGreedyIterM = np.inf, True
self._preliminaryTrainingMarginal()
if self.errorEstimatorKindMarginal == "LEAVE_ONE_OUT":
muidx = []
else:#if self.errorEstimatorKindMarginal[: 10] == "LOOK_AHEAD":
muidx = np.arange(len(self.trainedModel.data.musMarginal))
while self.firstGreedyIterM or (max2ErrorEst > self.greedyTolMarginal
and self.samplerMarginalGrid.npoints < self.maxIterMarginal):
errorEstTest, muidx, maxErrorEst, mu = \
self.greedyNextSampleMarginal(muidx, plotEst)
if len(maxErrorEst) > 0:
max2ErrorEst = np.max(maxErrorEst)
vbMng(self, "MAIN", ("Uniform testing error estimate "
"{:.4e}.").format(max2ErrorEst), 3)
else:
max2ErrorEst = 0.
if plotEst == "LAST":
self.plotEstimatorMarginal(errorEstTest, muidx, maxErrorEst)
vbMng(self, "DEL",
("Done computing snapshots (final snapshot count: "
"{}).").format(np.sum(self.samplingEngine.nsamples)), 3)
if (self.errorEstimatorKindMarginal == "LOOK_AHEAD_RECOVER"
and len(self.trainedModel._idxExcl) > 0):
vbMng(self, "INIT", "Recovering {} test models.".format(
len(self.trainedModel._idxExcl)), 7)
for j, mu in zip(self.trainedModel._idxExcl,
self.trainedModel._musMExcl):
self.musMarginal.insert(mu, j)
self.trainedModel.updateEffectiveSamples(self.HFEngine, [],
self.matchingWeight,
self.POD == PODGlobal, False)
self._finalizeMarginalization()
self._SMarginal = len(self.musMarginal)
self._approxParameters["SMarginal"] = self.SMarginal
self.trainedModel.data.approxParameters["SMarginal"] = (
self.SMarginal)
vbMng(self, "DEL", "Done recovering test models.", 7)
return 0
def checkComputedApprox(self) -> bool:
return (super().checkComputedApprox()
and len(self.mus) == len(self.trainedModel.data.mus))
def checkComputedApproxPivoted(self) -> bool:
return (super().checkComputedApprox()
and len(self.musMarginal) == len(self.trainedModel.data.musMarginal))
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 9305ec1..c277e2d 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,316 +1,310 @@
#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 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 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
__all__ = ['RationalInterpolantGreedyPivotedGreedy']
class RationalInterpolantGreedyPivotedGreedy(GenericPivotedGreedyApproximant,
RationalInterpolantGreedyPivoted):
"""
ROM greedy pivoted greedy rational interpolant 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;
- 'cutOffTolerance': tolerance for ignoring parasitic poles;
defaults to np.inf;
- 'cutOffSharedRatio': required ratio of marginal points to share
resonance in cut off strategy; 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;
- 'samplerMarginalGrid': 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'
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;
- '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';
- 'MMarginal': degree of marginal interpolant; 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;
- 'polydegreetypeMarginal': type of polynomial degree for marginal;
defaults to 'TOTAL';
- 'radialDirectionalWeightsMarginal': radial basis weights for
marginal interpolant; defaults to 1;
- - '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.
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;
- 'cutOffTolerance': tolerance for ignoring parasitic poles;
- 'cutOffSharedRatio': required ratio of marginal points to share
resonance in cut off strategy;
- '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;
- '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;
- 'MMarginal': degree of marginal interpolant;
- 'greedyTolMarginal': uniform error tolerance for marginal greedy
algorithm;
- 'maxIterMarginal': maximum number of marginal greedy steps;
- '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.
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;
- 'samplerMarginalGrid': 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.
cutOffTolerance: Tolerance for ignoring parasitic poles.
cutOffSharedRatio: Required ratio of marginal points to share resonance
in cut off strategy.
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.
samplerMarginalGrid: 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.
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.
MMarginal: Degree of marginal interpolant.
greedyTolMarginal: Uniform error tolerance for marginal greedy
algorithm.
maxIterMarginal: Maximum number of marginal greedy steps.
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.
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.
"""
@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.data[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 _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(
muTestBasePivot ** self.HFEngine.rescalingExp[self.directionPivot[0]],
musPivot ** self.HFEngine.rescalingExp[self.directionPivot[0]],
1e-10 * self.scaleFactor[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):
self.mus.data[k, self.directionPivot] = musPivot[k].data
self.mus.data[k, self.directionMarginal] = self.musMargLoc[-1].data
for k in range(len(muTestBasePivot)):
self.muTest.data[k, self.directionPivot] = muTestBasePivot[k].data
self.muTest.data[k, self.directionMarginal] = (
self.musMargLoc[-1].data)
self.muTest.data[-1, self.directionPivot] = musPivot[-1].data
self.muTest.data[-1, self.directionMarginal] = self.musMargLoc[-1].data
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"
Qs, Ps = self._preSetupApproxPivoted(mus)
S0 = copy(self.S)
for j, mu in enumerate(mus):
RationalInterpolantGreedy.setupSampling(self)
self.trainedModel = None
self.musMargLoc += [mu]
RationalInterpolantGreedy.setupApprox(self, self._plotEstPivot)
self.samplingEngs[j] = copy(self.samplingEngine)
Qs[j] = copy(self.trainedModel.data.Q)
Ps[j] = copy(self.trainedModel.data.P)
self._S = S0
self._postSetupApproxPivoted(mus, Qs, Ps)
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
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 bec9970..e80d12a 100644
--- a/rrompy/reduction_methods/pivoted/greedy/rational_interpolant_pivoted_greedy.py
+++ b/rrompy/reduction_methods/pivoted/greedy/rational_interpolant_pivoted_greedy.py
@@ -1,262 +1,250 @@
# 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 .generic_pivoted_greedy_approximant import GenericPivotedGreedyApproximant
from rrompy.reduction_methods.standard import RationalInterpolant
from rrompy.reduction_methods.pivoted import 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 checkParameterList, emptyParameterList
__all__ = ['RationalInterpolantPivotedGreedy']
class RationalInterpolantPivotedGreedy(GenericPivotedGreedyApproximant,
RationalInterpolantPivoted):
"""
ROM pivoted greedy rational interpolant 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;
- 'cutOffTolerance': tolerance for ignoring parasitic poles;
defaults to np.inf;
- 'cutOffSharedRatio': required ratio of marginal points to share
resonance in cut off strategy; 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;
- 'samplerMarginalGrid': 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'
and 'LEGENDRE'; 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;
- 'MMarginal': degree of marginal interpolant; 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;
- 'polydegreetypeMarginal': type of polynomial degree for marginal;
defaults to 'TOTAL';
- 'radialDirectionalWeights': radial basis weights for pivot
numerator; defaults to 1;
- 'radialDirectionalWeightsMarginal': radial basis weights for
marginal interpolant; defaults to 1;
- - 'nNearestNeighbor': number of pivot nearest neighbors considered
- if polybasis allows; defaults to -1;
- - '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;
- 'correctorForce': whether corrector should forcefully delete bad
poles; defaults to False;
- 'correctorTol': tolerance for corrector step; defaults to 0.,
i.e. no bad poles;
- 'correctorMaxIter': maximum number of corrector iterations;
defaults to 1.
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;
- 'cutOffTolerance': tolerance for ignoring parasitic poles;
- 'cutOffSharedRatio': required ratio of marginal points to share
resonance in cut off strategy;
- '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;
- 'M': degree of rational interpolant numerator;
- 'N': degree of rational interpolant denominator;
- 'MMarginal': degree of marginal interpolant;
- 'greedyTolMarginal': uniform error tolerance for marginal greedy
algorithm;
- 'maxIterMarginal': maximum number of marginal greedy steps;
- '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;
- 'samplerMarginalGrid': 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.
cutOffTolerance: Tolerance for ignoring parasitic poles.
cutOffSharedRatio: Required ratio of marginal points to share resonance
in cut off strategy.
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.
samplerMarginalGrid: 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.
M: Degree of rational interpolant numerator.
N: Degree of rational interpolant denominator.
MMarginal: Degree of marginal interpolant.
greedyTolMarginal: Uniform error tolerance for marginal greedy
algorithm.
maxIterMarginal: Maximum number of marginal greedy steps.
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.
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.
"""
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
self.musPivot = self.samplerPivot.generatePoints(self.S)
while len(self.musPivot) > self.S: self.musPivot.pop()
self.mus = emptyParameterList()
self.mus.reset((self.S, self.HFEngine.npar))
self.samplingEngine.resetHistory()
for k in range(self.S):
self.mus.data[k, self.directionPivot] = self.musPivot[k].data
self.mus.data[k, self.directionMarginal] = self.musMargLoc[-1].data
self.samplingEngine.iterSample(self.mus)
vbMng(self, "DEL", "Done computing snapshots.", 5)
self._m_selfmus = copy(self.mus)
self._mus = self.musPivot
self._m_mu0 = copy(self.mu0)
self._m_HFErescalingExp = copy(self.HFEngine.rescalingExp)
self._mu0 = checkParameterList(self.mu0(self.directionPivot), 1)[0]
self.HFEngine.rescalingExp = [self.HFEngine.rescalingExp[
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)
Qs, Ps = self._preSetupApproxPivoted(mus)
for j, mu in enumerate(mus):
RationalInterpolant.setupSampling(self)
self.trainedModel = None
self.musMargLoc += [mu]
RationalInterpolant.setupApprox(self)
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
self.samplingEngs[j] = copy(self.samplingEngine)
Qs[j] = copy(self.trainedModel.data.Q)
Ps[j] = copy(self.trainedModel.data.P)
self._postSetupApproxPivoted(mus, Qs, Ps)
vbMng(self, "DEL", "Done setting up pivoted approximant.", 10)
return 0
diff --git a/rrompy/reduction_methods/pivoted/rational_interpolant_greedy_pivoted.py b/rrompy/reduction_methods/pivoted/rational_interpolant_greedy_pivoted.py
index 6d3dc94..8b28a11 100644
--- a/rrompy/reduction_methods/pivoted/rational_interpolant_greedy_pivoted.py
+++ b/rrompy/reduction_methods/pivoted/rational_interpolant_greedy_pivoted.py
@@ -1,447 +1,439 @@
# 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, PODGlobal
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.numerical import dot
from rrompy.utilities.numerical.degree import totalDegreeN
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;
- 'scaleFactorDer': scaling factors for derivative computation;
defaults to 'AUTO';
- 'matchingWeight': weight for pole matching optimization; defaults
to 1;
- 'cutOffTolerance': tolerance for ignoring parasitic poles;
defaults to np.inf;
- 'cutOffSharedRatio': required ratio of marginal points to share
resonance in cut off strategy; 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'
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;
- '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';
- 'MMarginal': degree of marginal interpolant; defaults to 'AUTO',
i.e. maximum allowed;
- 'polydegreetypeMarginal': type of polynomial degree for marginal;
defaults to 'TOTAL';
- 'radialDirectionalWeightsMarginal': radial basis weights for
marginal interpolant; defaults to 1;
- - '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.
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;
- 'cutOffTolerance': tolerance for ignoring parasitic poles;
- 'cutOffSharedRatio': required ratio of marginal points to share
resonance in cut off strategy;
- '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;
- '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;
- '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.
scaleFactorDer: Scaling factors for derivative computation.
matchingWeight: Weight for pole matching optimization.
cutOffTolerance: Tolerance for ignoring parasitic poles.
cutOffSharedRatio: Required ratio of marginal points to share resonance
in cut off strategy.
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.
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.
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.
"""
def __init__(self, *args, **kwargs):
self._preInit()
self._addParametersToList(toBeExcluded = ["sampler"])
super().__init__(*args, **kwargs)
self._postInit()
@property
def tModelType(self):
if hasattr(self, "_temporaryPivot"):
return RationalInterpolantGreedy.tModelType.fget(self)
return super().tModelType
@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 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.")
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()
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(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)
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 = emptyParameterList()
self.muTest.reset((len(muTestBase) + 1, self.mus.shape[1]))
self.muTest.data[: -1] = muTestBase.data
self.muTest.data[-1] = muLast.data
self.M, self.N = ("AUTO",) * 2
def _finalizeSnapshots(self):
self.setupSampling()
self.samplingEngine.resetHistory(len(self.musMarginal))
for j in range(len(self.musMarginal)):
self.samplingEngine.setsample(self.samplingEngs[j].samples,
j, False)
self.samplingEngine.mus[j] = copy(self.samplingEngs[j].mus)
self.samplingEngine.musMarginal[j] = copy(self.musMarginal[j])
self.samplingEngine.nsamples[j] = self.samplingEngs[j].nsamples
if self.POD:
self.samplingEngine.RPOD[j] = self.samplingEngs[j].RPOD
self.samplingEngine.samples_full[j].data = (
self.samplingEngs[j].samples_full.data)
if self.POD == PODGlobal:
self.samplingEngine.coalesceSamples(self.interpRcondMarginal)
else:
self.samplingEngine.coalesceSamples()
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.musMarginal = self.samplerMarginal.generatePoints(self.SMarginal)
while len(self.musMarginal) > self.SMarginal: self.musMarginal.pop()
S0 = copy(self.S)
Qs, Ps = [None] * len(self.musMarginal), [None] * len(self.musMarginal)
self.samplingEngs = [None] * len(self.musMarginal)
self.computeScaleFactor()
self._scaleFactorOldPivot = copy(self.scaleFactor)
self.scaleFactor = self.scaleFactorPivot
self._temporaryPivot = 1
for j in range(len(self.musMarginal)):
vbMng(self, "MAIN",
"Building marginal model no. {} at {}.".format(j + 1,
self.musMarginal[j]), 5)
self._S = S0
self.musMargLoc = self.musMarginal[j]
RationalInterpolantGreedy.setupSampling(self)
self.trainedModel = None
self.verbosity -= 5
self.samplingEngine.verbosity -= 5
super().setupApprox(*args, **kwargs)
self.verbosity += 5
self.samplingEngine.verbosity += 5
self.samplingEngs[j] = copy(self.samplingEngine)
Qs[j] = copy(self.trainedModel.data.Q)
Ps[j] = copy(self.trainedModel.data.P)
self.scaleFactor = self._scaleFactorOldPivot
del self._scaleFactorOldPivot, self._temporaryPivot
self._finalizeSnapshots()
del self.musMargLoc, self.samplingEngs
self._mus = self.samplingEngine.musCoalesced
padLeft = 0
if self.POD != PODGlobal: padRight = self.samplingEngine.nsamplesTot
for j in range(len(self.musMarginal)):
nsj = self.samplingEngine.nsamples[j]
if self.POD == PODGlobal:
rRightj = self.samplingEngine.RPODCPart[:,
padLeft : padLeft + nsj]
Ps[j].postmultiplyTensorize(rRightj.T)
else:
padRight -= nsj
Ps[j].pad(padLeft, padRight)
padLeft += nsj
pMat = self.samplingEngine.samplesCoalesced.data
pMatEff = dot(self.HFEngine.C, pMat) if self.approx_state else 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.Qs, self.trainedModel.data.Ps = Qs, Ps
vbMng(self, "INIT", "Matching poles.", 10)
self.trainedModel.initializeFromRational(self.HFEngine,
self.matchingWeight,
self.POD == PODGlobal, False)
vbMng(self, "DEL", "Done matching poles.", 10)
self._finalizeMarginalization()
vbMng(self, "DEL", "Done setting up approximant.", 5)
return 0
diff --git a/rrompy/reduction_methods/pivoted/rational_interpolant_pivoted.py b/rrompy/reduction_methods/pivoted/rational_interpolant_pivoted.py
index a8ef072..5b03be1 100644
--- a/rrompy/reduction_methods/pivoted/rational_interpolant_pivoted.py
+++ b/rrompy/reduction_methods/pivoted/rational_interpolant_pivoted.py
@@ -1,389 +1,377 @@
# 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, PODGlobal
from rrompy.reduction_methods.standard.rational_interpolant import (
RationalInterpolant)
from rrompy.utilities.base import verbosityManager as vbMng
from rrompy.utilities.numerical import dot
from rrompy.utilities.numerical.hash_derivative import nextDerivativeIndices
from rrompy.utilities.exception_manager import RROMPyAssert, RROMPyWarning
from rrompy.parameter import emptyParameterList
__all__ = ['RationalInterpolantPivoted']
class RationalInterpolantPivoted(GenericPivotedApproximant,
RationalInterpolant):
"""
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;
- 'cutOffTolerance': tolerance for ignoring parasitic poles;
defaults to np.inf;
- 'cutOffSharedRatio': required ratio of marginal points to share
resonance in cut off strategy; 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'
and 'LEGENDRE'; 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;
- 'MMarginal': degree of marginal interpolant; defaults to 'AUTO',
i.e. maximum allowed;
- 'polydegreetypeMarginal': type of polynomial degree for marginal;
defaults to 'TOTAL';
- 'radialDirectionalWeights': radial basis weights for pivot
numerator; defaults to 1;
- 'radialDirectionalWeightsMarginal': radial basis weights for
marginal interpolant; defaults to 1;
- - 'nNearestNeighbor': number of pivot nearest neighbors considered
- if polybasis allows; defaults to -1;
- - '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;
- 'correctorForce': whether corrector should forcefully delete bad
poles; defaults to False;
- 'correctorTol': tolerance for corrector step; defaults to 0.,
i.e. no bad poles;
- 'correctorMaxIter': maximum number of corrector iterations;
defaults to 1.
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;
- 'cutOffTolerance': tolerance for ignoring parasitic poles;
- 'cutOffSharedRatio': required ratio of marginal points to share
resonance in cut off strategy;
- 'polybasis': type of polynomial basis for pivot
interpolation;
- 'polybasisMarginal': type of polynomial basis for marginal
interpolation;
- 'M': degree of rational interpolant numerator;
- 'N': degree of rational interpolant denominator;
- '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.
scaleFactorDer: Scaling factors for derivative computation.
matchingWeight: Weight for pole matching optimization.
cutOffTolerance: Tolerance for ignoring parasitic poles.
cutOffSharedRatio: Required ratio of marginal points to share resonance
in cut off strategy.
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.
M: Numerator degree of approximant.
N: Denominator degree of approximant.
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.
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.
"""
def __init__(self, *args, **kwargs):
self._preInit()
self._addParametersToList(toBeExcluded = ["polydegreetype", "sampler"])
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 hasattr(self, "_scaleFactorDer"):
scaleFactorDerold = self.scaleFactorDer
else: scaleFactorDerold = -1
if isinstance(scaleFactorDer, (str,)):
scaleFactorDer = scaleFactorDer.upper()
self._scaleFactorDer = scaleFactorDer
self._approxParameters["scaleFactorDer"] = self._scaleFactorDer
if scaleFactorDerold != self._scaleFactorDer: self.resetSamples()
@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 polybasis0(self):
if "_" in self.polybasis:
return self.polybasis.split("_")[0]
return self.polybasis
@property
def correctorTol(self):
"""Value of correctorTol."""
return self._correctorTol
@correctorTol.setter
def correctorTol(self, correctorTol):
if correctorTol < 0. or (correctorTol > 0. and self.nparPivot > 1):
RROMPyWarning(("Overriding prescribed corrector tolerance "
"to 0."))
correctorTol = 0.
self._correctorTol = correctorTol
self._approxParameters["correctorTol"] = self.correctorTol
@property
def correctorMaxIter(self):
"""Value of correctorMaxIter."""
return self._correctorMaxIter
@correctorMaxIter.setter
def correctorMaxIter(self, correctorMaxIter):
if correctorMaxIter < 1 or (correctorMaxIter > 1
and self.nparPivot > 1):
RROMPyWarning(("Overriding prescribed max number of corrector "
"iterations to 1."))
correctorMaxIter = 1
self._correctorMaxIter = correctorMaxIter
self._approxParameters["correctorMaxIter"] = self.correctorMaxIter
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 computeSnapshots(self):
"""Compute snapshots of solution map."""
RROMPyAssert(self._mode,
message = "Cannot start snapshot computation.")
self.computeScaleFactor()
if self.samplingEngine.nsamplesTot != self.S * self.SMarginal:
self.resetSamples()
self.samplingEngine.scaleFactor = self.scaleFactorDer
vbMng(self, "INIT", "Starting computation of snapshots.", 5)
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))
self.samplingEngine.resetHistory(self.SMarginal)
for j, muMarg in enumerate(self.musMarginal):
for k in range(j * self.S, (j + 1) * self.S):
self.mus.data[k, self.directionPivot] = (
self.musPivot[k - j * self.S].data)
self.mus.data[k, self.directionMarginal] = muMarg.data
self.samplingEngine.iterSample(self.musPivot, self.musMarginal)
self._finalizeSnapshots()
vbMng(self, "DEL", "Done computing snapshots.", 5)
def _finalizeSnapshots(self):
if self.POD == PODGlobal:
self.samplingEngine.coalesceSamples(self.interpRcondMarginal)
else:
self.samplingEngine.coalesceSamples()
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()
pMat = self.samplingEngine.samplesCoalesced.data
pMatEff = dot(self.HFEngine.C, pMat) if self.approx_state else pMat
if self.trainedModel is None:
self.trainedModel = self.tModelType()
self.trainedModel.verbosity = self.verbosity
self.trainedModel.timestamp = self.timestamp
datadict = {"mu0": self.mu0, "projMat": pMatEff,
"scaleFactor": self.scaleFactor,
"rescalingExp": self.HFEngine.rescalingExp,
"directionPivot": self.directionPivot}
self.trainedModel.data = self.initializeModelData(datadict)[0]
else:
self.trainedModel = self.trainedModel
self.trainedModel.data.projMat = copy(pMatEff)
N0 = copy(self.N)
Qs, Ps = [None] * len(self.musMarginal), [None] * len(self.musMarginal)
self._temporaryPivot = 1
padLeft = 0
if self.POD:
self._RPODOldPivot = copy(self.samplingEngine.RPODCoalesced)
else:
self._samplesOldPivot = copy(self.samplingEngine.samples)
padRight = self.samplingEngine.nsamplesTot
self._scaleFactorOldPivot = copy(self.scaleFactor)
self.scaleFactor = self.scaleFactorPivot
for j in range(len(self.musMarginal)):
vbMng(self, "MAIN",
"Building marginal model no. {} at {}.".format(j + 1,
self.musMarginal[j]), 5)
self.N = N0
if self.POD:
self.samplingEngine.RPOD = (
self._RPODOldPivot[:, padLeft : padLeft + self.S])
else:
self.samplingEngine.samples = self._samplesOldPivot[j]
padRight -= self.S
self.verbosity -= 5
self._iterCorrector()
self.verbosity += 5
Qs[j] = copy(self.trainedModel.data.Q)
Ps[j] = copy(self.trainedModel.data.P)
del self.trainedModel.data.Q, self.trainedModel.data.P
if not self.POD: Ps[j].pad(padLeft, padRight)
padLeft += self.S
if self.POD:
self.samplingEngine.RPODCoalesced = copy(self._RPODOldPivot)
del self._RPODOldPivot
else:
self.samplingEngine.samples = copy(self._samplesOldPivot)
del self._samplesOldPivot
self.scaleFactor = self._scaleFactorOldPivot
del self._temporaryPivot, self._scaleFactorOldPivot
self.trainedModel.data.mus = copy(self.mus)
self.trainedModel.data.musPivot = copy(self.musPivot)
self.trainedModel.data.musMarginal = copy(self.musMarginal)
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 == PODGlobal, False)
vbMng(self, "DEL", "Done matching poles.", 10)
self._finalizeMarginalization()
vbMng(self, "DEL", "Done setting up approximant.", 5)
return 0
diff --git a/rrompy/reduction_methods/standard/greedy/rational_interpolant_greedy.py b/rrompy/reduction_methods/standard/greedy/rational_interpolant_greedy.py
index 068b5e1..3952777 100644
--- a/rrompy/reduction_methods/standard/greedy/rational_interpolant_greedy.py
+++ b/rrompy/reduction_methods/standard/greedy/rational_interpolant_greedy.py
@@ -1,552 +1,551 @@
# 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 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 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 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';
- '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;
- '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;
- '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.
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.
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", "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",
- "nNearestNeighbor"])
+ "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 pvT(mus, deg, *args)
def getErrorEstimatorDiscrepancy(self, mus:Np1D) -> Np1D:
"""Discrepancy-based residual estimator."""
checkIfAffine(self.HFEngine, "apply discrepancy-based error 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)
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.mus ** self.HFEngine.rescalingExp
muTestEff = mus ** self.HFEngine.rescalingExp
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
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 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.data)
else:
app_muTestSample = dot(self.samplingEngine.samples,
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 hasattr(self.samplingEngine, "samples_full"):
uEx = self.samplingEngine.samples_full[-1]
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.data))
app_muTestSample = sampleList(self.HFEngine.applyC(
app_muTestSample.data))
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]
verb = 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)
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:
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 setupOK > 0:
err = np.empty(len(mus))
err[:] = np.nan
if not return_max: return err
return err, [- setupOK], 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[: 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 _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()
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.samples.data
pMatEff = dot(self.HFEngine.C, pMat) if self.approx_state else pMat
if self.trainedModel is None:
self.trainedModel = self.tModelType()
self.trainedModel.verbosity = self.verbosity
self.trainedModel.timestamp = self.timestamp
datadict = {"mu0": self.mu0, "projMat": pMatEff,
"scaleFactor": self.scaleFactor,
"rescalingExp": self.HFEngine.rescalingExp}
self.trainedModel.data = self.initializeModelData(datadict)[0]
else:
self.trainedModel = self.trainedModel
self.trainedModel.data.projMat = copy(pMatEff)
self.trainedModel.data.mus = copy(self.mus)
self.trainedModel.data.mus = copy(self.mus)
self.catchInstability = 2
unstable = False
if self.E > 0:
try:
Q = self._setupDenominator()[0]
except RROMPyException as RE:
RROMPyWarning("Downgraded {}: {}".format(RE.__class__.__name__,
RE))
vbMng(self, "DEL", "", 7)
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:
RROMPyWarning("Downgraded {}: {}".format(RE.__class__.__name__,
RE))
vbMng(self, "DEL", "", 7)
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._iterCorrector()
self.trainedModel.data.approxParameters = copy(
self.approxParameters)
return val
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/standard/rational_interpolant.py b/rrompy/reduction_methods/standard/rational_interpolant.py
index 4f189ae..9969d38 100644
--- a/rrompy/reduction_methods/standard/rational_interpolant.py
+++ b/rrompy/reduction_methods/standard/rational_interpolant.py
@@ -1,694 +1,677 @@
# Copyright (C) 2018 by the RROMPy authors
#
# This file is part of RROMPy.
#
# RROMPy is free software: you can redistribute it and/or modify
# it under the terms of the GNU Lesser General Public License as published by
# the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
#
# RROMPy is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU Lesser General Public License for more details.
#
# You should have received a copy of the GNU Lesser General Public License
# along with RROMPy. If not, see .
#
from copy import deepcopy as copy
import numpy as np
from rrompy.reduction_methods.base import checkRobustTolerance
from .generic_standard_approximant import GenericStandardApproximant
from rrompy.utilities.poly_fitting.polynomial import (
polybases as ppb, polyfitname,
polyvander as pvP, polyvanderTotal as pvTP,
polyTimes, polyTimesTable, vanderInvTable,
PolynomialInterpolator as PI)
from rrompy.utilities.poly_fitting.heaviside import rational2heaviside
from rrompy.utilities.poly_fitting.radial_basis import (polybases as rbpb,
RadialBasisInterpolator as RBI)
from rrompy.utilities.poly_fitting.moving_least_squares import (
polybases as mlspb,
MovingLeastSquaresInterpolator as MLSI)
from rrompy.utilities.base.types import Np1D, Np2D, Tuple, List, sampList
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;
- - 'nNearestNeighbor': mumber of nearest neighbors considered in
- numerator if polybasis allows; defaults to -1;
- 'interpRcond': tolerance for interpolation; defaults to None;
- 'robustTol': tolerance for robust rational denominator
management; defaults to 0;
- 'correctorForce': whether corrector should forcefully delete bad
poles; defaults to False;
- 'correctorTol': tolerance for corrector step; defaults to 0.,
i.e. no bad poles;
- 'correctorMaxIter': maximum number of corrector iterations;
defaults to 1.
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;
- - 'nNearestNeighbor': mumber of nearest neighbors considered in
- numerator if polybasis allows;
- 'interpRcond': tolerance for interpolation via numpy.polyfit;
- 'robustTol': tolerance for robust rational denominator
management;
- '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 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.
- nNearestNeighbor: Number of nearest neighbors considered in numerator
- if polybasis allows.
interpRcond: Tolerance for interpolation via numpy.polyfit.
robustTol: Tolerance for robust rational denominator management.
correctorForce: Whether corrector should forcefully delete bad poles.
correctorTol: Tolerance for corrector step.
correctorMaxIter: Maximum number of corrector iterations.
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",
- "nNearestNeighbor", "interpRcond",
+ "radialDirectionalWeights", "interpRcond",
"robustTol", "correctorForce",
"correctorTol", "correctorMaxIter"],
["MONOMIAL", "AUTO", "AUTO", "TOTAL", [1.],
- -1, -1, 0., False, 0., 1])
+ -1, 0., False, 0., 1])
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 + mlspb:
raise RROMPyException("Prescribed polybasis not recognized.")
self._polybasis = polybasis
except:
RROMPyWarning(("Prescribed polybasis not recognized. Overriding "
"to 'MONOMIAL'."))
self._polybasis = "MONOMIAL"
self._approxParameters["polybasis"] = self.polybasis
@property
def polybasis0(self):
if "_" in self.polybasis:
return self.polybasis.split("_")[0]
return self.polybasis
@property
def interpRcond(self):
"""Value of interpRcond."""
return self._interpRcond
@interpRcond.setter
def interpRcond(self, interpRcond):
self._interpRcond = interpRcond
self._approxParameters["interpRcond"] = self.interpRcond
@property
def radialDirectionalWeights(self):
"""Value of radialDirectionalWeights."""
return self._radialDirectionalWeights
@radialDirectionalWeights.setter
def radialDirectionalWeights(self, radialDirectionalWeights):
if not hasattr(radialDirectionalWeights, "__len__"):
radialDirectionalWeights = [radialDirectionalWeights]
self._radialDirectionalWeights = radialDirectionalWeights
self._approxParameters["radialDirectionalWeights"] = (
self.radialDirectionalWeights)
- @property
- def nNearestNeighbor(self):
- """Value of nNearestNeighbor."""
- return self._nNearestNeighbor
- @nNearestNeighbor.setter
- def nNearestNeighbor(self, nNearestNeighbor):
- self._nNearestNeighbor = nNearestNeighbor
- self._approxParameters["nNearestNeighbor"] = self.nNearestNeighbor
-
@property
def M(self):
"""Value of M."""
return self._M
@M.setter
def M(self, M):
if 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 correctorForce(self):
"""Value of correctorForce."""
return self._correctorForce
@correctorForce.setter
def correctorForce(self, correctorForce):
self._correctorForce = correctorForce
self._approxParameters["correctorForce"] = self.correctorForce
@property
def correctorTol(self):
"""Value of correctorTol."""
return self._correctorTol
@correctorTol.setter
def correctorTol(self, correctorTol):
if correctorTol < 0. or (correctorTol > 0. and self.npar > 1):
RROMPyWarning(("Overriding prescribed corrector tolerance "
"to 0."))
correctorTol = 0.
self._correctorTol = correctorTol
self._approxParameters["correctorTol"] = self.correctorTol
@property
def correctorMaxIter(self):
"""Value of correctorMaxIter."""
return self._correctorMaxIter
@correctorMaxIter.setter
def correctorMaxIter(self, correctorMaxIter):
if correctorMaxIter < 1 or (correctorMaxIter > 1 and self.npar > 1):
RROMPyWarning(("Overriding prescribed max number of corrector "
"iterations to 1."))
correctorMaxIter = 1
self._correctorMaxIter = correctorMaxIter
self._approxParameters["correctorMaxIter"] = self.correctorMaxIter
def resetSamples(self):
"""Reset samples."""
super().resetSamples()
self._musUniqueCN = None
self._derIdxs = None
self._reorder = None
def _setupInterpolationIndices(self):
"""Setup parameters for polyvander."""
if self._musUniqueCN is None or len(self._reorder) != len(self.mus):
self._musUniqueCN, musIdxsTo, musIdxs, musCount = (
self.trainedModel.centerNormalize(self.mus).unique(
return_index = True, return_inverse = True,
return_counts = True))
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()
idxSamplesEff = list(range(self.S))
if self.POD:
ev, eV = self.findeveVGQR(
self.samplingEngine.RPOD[:, idxSamplesEff], invD)
else:
ev, eV = self.findeveVGExplicit(
self.samplingEngine.samples(idxSamplesEff), invD)
nevBad = checkRobustTolerance(ev, self.robustTol)
if nevBad <= 1: break
if self.catchInstability > 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])
else:
q.coeffs = eV[:, 0].reshape([self.N + 1] * self.npar)
vbMng(self, "DEL", "Done computing denominator.", 7)
return q, fitinv
def _setupNumerator(self):
"""Compute rational numerator."""
RROMPyAssert(self._mode, message = "Cannot setup numerator.")
vbMng(self, "INIT", "Starting computation of numerator.", 7)
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.verbosity >= 5, self.polydegreetype == "TOTAL",
{"derIdxs": self._derIdxs, "reorder": self._reorder,
"scl": self.scaleFactorRel}]
if self.polybasis in ppb:
p = PI()
else:
pParRest = [self.radialDirectionalWeights] + pParRest
- pParRest[-1]["nNearestNeighbor"] = self.nNearestNeighbor
p = RBI() if self.polybasis in rbpb else MLSI()
if self.polybasis in ppb + rbpb:
pParRest += [{"rcond": self.interpRcond}]
wellCond, msg = p.setupByInterpolation(self._musUniqueCN,
Qevaldiag, self.M,
self.polybasis, *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()
pMat = self.samplingEngine.samples.data
pMatEff = dot(self.HFEngine.C, pMat) if self.approx_state else pMat
if self.trainedModel is None:
self.trainedModel = self.tModelType()
self.trainedModel.verbosity = self.verbosity
self.trainedModel.timestamp = self.timestamp
datadict = {"mu0": self.mu0, "projMat": pMatEff,
"scaleFactor": self.scaleFactor,
"rescalingExp": self.HFEngine.rescalingExp}
self.trainedModel.data = self.initializeModelData(datadict)[0]
else:
self.trainedModel = self.trainedModel
self.trainedModel.data.projMat = copy(pMatEff)
self.trainedModel.data.mus = copy(self.mus)
self._iterCorrector()
self.trainedModel.data.approxParameters = copy(self.approxParameters)
vbMng(self, "DEL", "Done setting up approximant.", 5)
return 0
def _iterCorrector(self):
if self.correctorTol > 0. and (self.correctorMaxIter > 1
or self.correctorForce):
vbMng(self, "INIT", "Starting correction iterations.", 5)
self._Qhat = PI()
self._Qhat.npar = self.npar
self._Qhat.polybasis = "MONOMIAL"
self._Qhat.coeffs = np.ones(1, dtype = np.complex)
if self.POD:
self._RPODOld = copy(self.samplingEngine.RPOD)
else:
self._samplesOld = copy(self.samplingEngine.samples)
else: vbMng(self, "INIT", "Starting approximant finalization.", 5)
for j in range(self.correctorMaxIter):
if self.N > 0 or (hasattr(self, "_N_isauto")
and self.S > self.npar):
Q = self._setupDenominator()[0]
else:
Q = PI()
Q.coeffs = np.ones((1,) * self.npar, dtype = np.complex)
Q.npar = self.npar
Q.polybasis = self.polybasis
self.N = 0
if j == 0: _N0 = self.N
self.trainedModel.data.Q = Q
self.trainedModel.data.P = self._setupNumerator()
self._applyCorrector(j)
if self.N <= 0: break
self.N = _N0
if self.correctorTol <= 0. or (self.correctorMaxIter <= 1
and not self.correctorForce):
vbMng(self, "DEL", "Terminated approximant finalization.", 5)
return
if self.POD:
self.samplingEngine.RPOD = self._RPODOld
del self._RPODOld
else:
self.samplingEngine.samples = self._samplesOld
del self._samplesOld
if self.correctorForce:
QOld, QOldBasis = [1.], "MONOMIAL"
else:
QOld, QOldBasis = Q.coeffs, self.polybasis
Q = polyTimes(self._Qhat.coeffs, QOld, Pbasis = self._Qhat.polybasis,
Qbasis = QOldBasis, Rbasis = self.polybasis)
del self._Qhat
gamma = np.linalg.norm(Q)
self.trainedModel.data.Q.coeffs = np.pad(Q, (0, self.N - len(Q) + 1),
"constant") / gamma
if self.correctorForce:
self.trainedModel.data.P = self._setupNumerator()
else:
self.trainedModel.data.P.coeffs /= gamma
vbMng(self, "DEL", "Terminated correction iterations.", 5)
def _applyCorrector(self, j:int):
if self.correctorTol <= 0. or (j >= self.correctorMaxIter - 1
and not self.correctorForce):
self.N = 0
return
cfs, pls, _ = rational2heaviside(self.trainedModel.data.P,
self.trainedModel.data.Q)
cfs = cfs[: self.N]
if self.POD:
resEff = np.linalg.norm(cfs, axis = 1)
else:
resEff = self.HFEngine.norm(self.samplingEngine.samples.dot(cfs.T),
is_state = self.approx_state)
goodPole = np.logical_not(np.isinf(pls))
potentialGood = (potential(pls[goodPole], self.sampler.normalFoci())
/ self.sampler.groundPotential())
potentialGood[potentialGood < 1.] = 1.
resEff[goodPole] /= potentialGood
resEff /= np.max(resEff)
idxKeep = np.where(resEff >= self.correctorTol)[0]
vbMng(self, "MAIN",
("Correction iteration no. {}: {} out of {} residuals satisfy "
"tolerance.").format(j + 1, len(idxKeep), self.N), 10)
self.N -= len(idxKeep)
if self.N <= 0 and not self.correctorForce: return
for i in idxKeep:
self._Qhat.coeffs = polyTimes(self._Qhat.coeffs, [- pls[i], 1.],
Pbasis = self._Qhat.polybasis,
Rbasis = self._Qhat.polybasis)
self._Qhat.coeffs /= np.linalg.norm(self._Qhat.coeffs)
if self.N <= 0: return
vbMng(self, "MAIN",
("Removing poles at (normalized positions): "
"{}.").format(pls[resEff < self.correctorTol]), 65)
That = polyTimesTable(self._Qhat, self._musUniqueCN,
self._reorder, self._derIdxs,
self.scaleFactorRel).T
if self.POD:
self.samplingEngine.RPOD = self._RPODOld.dot(That)
else:
self.samplingEngine.samples = self._samplesOld.dot(That)
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()
TEGen = pvTP if self.polydegreetype == "TOTAL" else pvP
TEGenPar = [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 = TEGen(self._musUniqueCN, Eeff, *TEGenPar)
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 = TEGen(self._musUniqueCN, Neff, *TEGenPar)
for k in range(len(invD)): invD[k] = dot(invD[k], TN)
return invD, fitinv
def findeveVGExplicit(self, sampleE:sampList,
invD:List[Np2D]) -> Tuple[Np1D, Np2D]:
"""
Compute explicitly eigenvalues and eigenvectors of rational denominator
matrix.
"""
RROMPyAssert(self._mode, message = "Cannot solve eigenvalue problem.")
nEnd = invD[0].shape[1]
eWidth = len(invD)
vbMng(self, "INIT", "Building gramian matrix.", 10)
gramian = self.HFEngine.innerProduct(sampleE, sampleE,
is_state = self.approx_state)
G = np.zeros((nEnd, nEnd), dtype = np.complex)
for k in range(eWidth):
G += dot(dot(gramian, invD[k]).T, invD[k].conj()).T
vbMng(self, "DEL", "Done building gramian.", 10)
vbMng(self, "INIT", "Solving eigenvalue problem for gramian matrix.",
7)
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)
return ev, eV
def findeveVGQR(self, RPODE:Np2D, invD:List[Np2D]) -> Tuple[Np1D, Np2D]:
"""
Compute eigenvalues and eigenvectors of rational denominator matrix
through SVD of R factor.
"""
RROMPyAssert(self._mode, message = "Cannot solve eigenvalue problem.")
nEnd = invD[0].shape[1]
S = RPODE.shape[0]
eWidth = len(invD)
vbMng(self, "INIT", "Building half-gramian matrix stack.", 10)
Rstack = np.zeros((S * eWidth, nEnd), dtype = np.complex)
for k in range(eWidth):
Rstack[k * S : (k + 1) * S, :] = dot(RPODE, invD[k])
vbMng(self, "DEL", "Done building half-gramian.", 10)
vbMng(self, "INIT", "Solving svd for square root of gramian matrix.",
7)
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)
return ev, eV
def getResidues(self, *args, **kwargs) -> Np1D:
"""
Obtain approximant residues.
Returns:
Matrix with residues as columns.
"""
return self.trainedModel.getResidues(*args, **kwargs)
diff --git a/rrompy/reduction_methods/standard/rational_moving_least_squares.py b/rrompy/reduction_methods/standard/rational_moving_least_squares.py
index d965c4a..fecdbec 100644
--- a/rrompy/reduction_methods/standard/rational_moving_least_squares.py
+++ b/rrompy/reduction_methods/standard/rational_moving_least_squares.py
@@ -1,337 +1,312 @@
# Copyright (C) 2018 by the RROMPy authors
#
# This file is part of RROMPy.
#
# RROMPy is free software: you can redistribute it and/or modify
# it under the terms of the GNU Lesser General Public License as published by
# the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
#
# RROMPy is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU Lesser General Public License for more details.
#
# You should have received a copy of the GNU Lesser General Public License
# along with RROMPy. If not, see .
#
from copy import deepcopy as copy
import numpy as np
from .rational_interpolant import RationalInterpolant
from rrompy.utilities.poly_fitting.polynomial import (polybases as ppb,
polyvander as pvP,
polyvanderTotal as pvTP)
from rrompy.utilities.base.types import Np2D
from rrompy.utilities.base import verbosityManager as vbMng
from rrompy.utilities.numerical import dot
from rrompy.utilities.numerical.degree import (fullDegreeMaxMask,
totalDegreeMaxMask)
from rrompy.utilities.exception_manager import (RROMPyException, RROMPyAssert,
RROMPyWarning)
__all__ = ['RationalMovingLeastSquares']
class RationalMovingLeastSquares(RationalInterpolant):
"""
ROM rational moving LS interpolant computation for parametric problems.
Args:
HFEngine: HF problem solver.
mu0(optional): Default parameter. Defaults to 0.
approxParameters(optional): Dictionary containing values for main
parameters of approximant. Recognized keys are:
- 'POD': whether to compute POD of snapshots; defaults to True;
- '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';
- 'radialBasis': numerator radial basis type; defaults to
'GAUSSIAN';
- 'radialDirectionalWeights': radial basis weights for interpolant
numerator; defaults to 1;
- - 'nNearestNeighbor': number of nearest neighbors considered in
- numerator if radialBasis allows; defaults to -1;
- 'radialBasisDen': denominator radial basis type; defaults to
'GAUSSIAN';
- 'radialDirectionalWeightsDen': radial basis weights for
interpolant denominator; defaults to 1;
- - 'nNearestNeighborDen': number of nearest neighbors considered in
- denominator if radialBasisDen allows; defaults to -1;
- 'interpRcond': tolerance for interpolation; defaults to None;
- 'robustTol': tolerance for robust rational denominator
management; defaults to 0.
Defaults to empty dict.
approx_state(optional): Whether to approximate state. Defaults to
False.
verbosity(optional): Verbosity level. Defaults to 10.
Attributes:
HFEngine: HF problem solver.
mu0: Default parameter.
mus: Array of snapshot parameters.
approxParameters: Dictionary containing values for main parameters of
approximant. Recognized keys are in parameterList.
parameterListSoft: Recognized keys of soft approximant parameters:
- 'POD': whether to compute POD of snapshots;
- '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;
- 'radialBasis': numerator radial basis type;
- 'radialDirectionalWeights': radial basis weights for interpolant
numerator;
- - 'nNearestNeighbor': number of nearest neighbors considered in
- numerator if radialBasis allows;
- 'radialBasisDen': denominator radial basis type;
- 'radialDirectionalWeightsDen': radial basis weights for
interpolant denominator;
- - 'nNearestNeighborDen': number of nearest neighbors considered in
- denominator if radialBasisDen allows;
- 'interpRcond': tolerance for interpolation via numpy.polyfit;
- 'robustTol': tolerance for robust rational denominator
management.
parameterListCritical: Recognized keys of critical approximant
parameters:
- 'S': total number of samples current approximant relies upon;
- 'sampler': sample point generator.
approx_state: Whether to approximate state.
verbosity: Verbosity level.
POD: Whether to compute POD of snapshots.
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.
radialBasis: Numerator radial basis type.
radialDirectionalWeights: Radial basis weights for interpolant
numerator.
- nNearestNeighbor: Number of nearest neighbors considered in numerator
- if radialBasis allows.
radialBasisDen: Denominator radial basis type.
radialDirectionalWeightsDen: Radial basis weights for interpolant
denominator.
- nNearestNeighborDen: Number of nearest neighbors considered in
- denominator if radialBasisDen allows.
interpRcond: Tolerance for interpolation via numpy.polyfit.
robustTol: Tolerance for robust rational denominator management.
muBounds: list of bounds for parameter values.
samplingEngine: Sampling engine.
uHF: High fidelity solution(s) with parameter(s) lastSolvedHF as
sampleList.
lastSolvedHF: Parameter(s) corresponding to last computed high fidelity
solution(s) as parameterList.
uApproxReduced: Reduced approximate solution(s) with parameter(s)
lastSolvedApprox as sampleList.
lastSolvedApproxReduced: Parameter(s) corresponding to last computed
reduced approximate solution(s) as parameterList.
uApprox: Approximate solution(s) with parameter(s) lastSolvedApprox as
sampleList.
lastSolvedApprox: Parameter(s) corresponding to last computed
approximate solution(s) as parameterList.
Q: Numpy 1D vector containing complex coefficients of approximant
denominator.
P: Numpy 2D vector whose columns are FE dofs of coefficients of
approximant numerator.
"""
def __init__(self, *args, **kwargs):
self._preInit()
self._addParametersToList(["radialBasis", "radialBasisDen",
- "radialDirectionalWeightsDen",
- "nNearestNeighborDen"],
- ["GAUSSIAN", "GAUSSIAN", 1., -1],
+ "radialDirectionalWeightsDen"],
+ ["GAUSSIAN", "GAUSSIAN", 1.],
toBeExcluded = ["correctorForce",
"correctorTol",
"correctorMaxIter"])
super().__init__(*args, **kwargs)
self.catchInstability = 0
self._postInit()
@property
def correctorForce(self):
"""Value of correctorForce."""
return False
@correctorForce.setter
def correctorForce(self, correctorForce):
RROMPyWarning(("correctorForce is used just to simplify inheritance, "
"and its value cannot be changed from False."))
@property
def correctorTol(self):
"""Value of correctorTol."""
return 0.
@correctorTol.setter
def correctorTol(self, correctorTol):
RROMPyWarning(("correctorTol is used just to simplify inheritance, "
"and its value cannot be changed from 0."))
@property
def correctorMaxIter(self):
"""Value of correctorMaxIter."""
return 1
@correctorMaxIter.setter
def correctorMaxIter(self, correctorMaxIter):
RROMPyWarning(("correctorMaxIter is used just to simplify "
"inheritance, and its value cannot be changed from 1."))
@property
def tModelType(self):
from .trained_model.trained_model_rational_mls import (
TrainedModelRationalMLS)
return TrainedModelRationalMLS
@property
def polybasis(self):
"""Value of polybasis."""
return self._polybasis
@polybasis.setter
def polybasis(self, polybasis):
try:
polybasis = polybasis.upper().strip().replace(" ","")
if polybasis not in ppb:
raise RROMPyException("Prescribed polybasis not recognized.")
self._polybasis = polybasis
except:
RROMPyWarning(("Prescribed polybasis not recognized. Overriding "
"to 'MONOMIAL'."))
self._polybasis = "MONOMIAL"
self._approxParameters["polybasis"] = self.polybasis
@property
def radialBasis(self):
"""Value of radialBasis."""
return self._radialBasis
@radialBasis.setter
def radialBasis(self, radialBasis):
self._radialBasis = radialBasis
self._approxParameters["radialBasis"] = self.radialBasis
@property
def radialBasisDen(self):
"""Value of radialBasisDen."""
return self._radialBasisDen
@radialBasisDen.setter
def radialBasisDen(self, radialBasisDen):
self._radialBasisDen = radialBasisDen
self._approxParameters["radialBasisDen"] = self.radialBasisDen
@property
def radialDirectionalWeightsDen(self):
"""Value of radialDirectionalWeightsDen."""
return self._radialDirectionalWeightsDen
@radialDirectionalWeightsDen.setter
def radialDirectionalWeightsDen(self, radialDirectionalWeightsDen):
self._radialDirectionalWeightsDen = radialDirectionalWeightsDen
self._approxParameters["radialDirectionalWeightsDen"] = (
self.radialDirectionalWeightsDen)
- @property
- def nNearestNeighborDen(self):
- """Value of nNearestNeighborDen."""
- return self._nNearestNeighborDen
- @nNearestNeighborDen.setter
- def nNearestNeighborDen(self, nNearestNeighborDen):
- self._nNearestNeighborDen = nNearestNeighborDen
- self._approxParameters["nNearestNeighborDen"] = (
- self.nNearestNeighborDen)
-
def _setupDenominator(self) -> Np2D:
"""Compute rational denominator."""
RROMPyAssert(self._mode, message = "Cannot setup denominator.")
vbMng(self, "INIT",
"Starting computation of denominator-related blocks.", 7)
self._setupInterpolationIndices()
pPar = [self._musUniqueCN, self.N, self.polybasis0, self._derIdxs,
self._reorder, self.scaleFactorRel]
if self.polydegreetype == "TOTAL":
TN = pvTP(*pPar)
else: #if self.polydegreetype == "FULL":
pPar[1] = [pPar[1]] * self.npar
TN = pvP(*pPar)
TNTen = np.zeros((self.S, self.S, TN.shape[1]), dtype = TN.dtype)
TNTen[np.arange(self.S), np.arange(self.S)] = TN
if self.POD: TNTen = dot(self.samplingEngine.RPOD, TNTen)
vbMng(self, "DEL", "Done computing denominator-related blocks.", 7)
return TN, TNTen
def _setupNumerator(self) -> Np2D:
"""Compute rational numerator."""
RROMPyAssert(self._mode, message = "Cannot setup numerator.")
vbMng(self, "INIT",
"Starting computation of denominator-related blocks.", 7)
self._setupInterpolationIndices()
pPar = [self._musUniqueCN, self.M, self.polybasis0, self._derIdxs,
self._reorder, self.scaleFactorRel]
if self.polydegreetype == "TOTAL":
TM = pvTP(*pPar)
else: #if self.polydegreetype == "FULL":
pPar[1] = [pPar[1]] * self.npar
TM = pvP(*pPar)
vbMng(self, "DEL", "Done computing denominator-related blocks.", 7)
return TM
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()
pMat = self.samplingEngine.samples.data
pMatEff = dot(self.HFEngine.C, pMat) if self.approx_state else pMat
if self.trainedModel is None:
self.trainedModel = self.tModelType()
self.trainedModel.verbosity = self.verbosity
self.trainedModel.timestamp = self.timestamp
datadict = {"mu0": self.mu0, "projMat": pMatEff,
"scaleFactor": self.scaleFactor,
"rescalingExp": self.HFEngine.rescalingExp}
data = self.initializeModelData(datadict)[0]
data.POD = self.POD
data.polybasis = self.polybasis
data.polydegreetype = self.polydegreetype
data.radialBasis = self.radialBasis
data.radialWeights = self.radialDirectionalWeights
- data.nNearestNeighbor = self.nNearestNeighbor
data.radialBasisDen = self.radialBasisDen
data.radialWeightsDen = self.radialDirectionalWeightsDen
- data.nNearestNeighborDen = self.nNearestNeighborDen
data.interpRcond = self.interpRcond
self.trainedModel.data = data
else:
self.trainedModel = self.trainedModel
self.trainedModel.data.projMat = copy(pMatEff)
if not self.POD:
self.trainedModel.data.gramian = self.HFEngine.innerProduct(
self.samplingEngine.samples,
self.samplingEngine.samples,
is_state = self.approx_state)
self.trainedModel.data.mus = copy(self.mus)
self.trainedModel.data.M = self.M
self.trainedModel.data.N = self.N
QVan, self.trainedModel.data.QBlocks = self._setupDenominator()
self.trainedModel.data.PVan = self._setupNumerator()
if self.polydegreetype == "TOTAL":
degreeMaxMask = totalDegreeMaxMask
else: #if self.polydegreetype == "FULL":
degreeMaxMask = fullDegreeMaxMask
if self.N > self.M:
self.trainedModel.data.QVan = QVan
self.trainedModel.data.domQIdxs = degreeMaxMask(self.N, self.npar)
else:
self.trainedModel.data.domQIdxs = degreeMaxMask(self.M, self.npar)
self.trainedModel.data.approxParameters = copy(self.approxParameters)
vbMng(self, "DEL", "Done setting up approximant.", 5)
return 0
diff --git a/rrompy/reduction_methods/standard/rational_pade.py b/rrompy/reduction_methods/standard/rational_pade.py
index 3260f7f..d75461a 100644
--- a/rrompy/reduction_methods/standard/rational_pade.py
+++ b/rrompy/reduction_methods/standard/rational_pade.py
@@ -1,317 +1,309 @@
# 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 .rational_interpolant import RationalInterpolant
from rrompy.utilities.poly_fitting.polynomial import (
polybases as ppb, polyfitname,
polyvander as pvP, polyvanderTotal as pvTP,
polyTimesTable, vanderInvTable,
PolynomialInterpolator as PI)
from rrompy.utilities.poly_fitting.radial_basis import (polybases as rbpb,
RadialBasisInterpolator as RBI)
from rrompy.utilities.poly_fitting.moving_least_squares import (
MovingLeastSquaresInterpolator as MLSI)
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;
- - 'nNearestNeighbor': mumber of nearest neighbors considered in
- numerator if polybasis allows; defaults to -1;
- 'interpRcond': tolerance for interpolation; defaults to None;
- 'robustTol': tolerance for robust rational denominator
management; defaults to 0;
- 'correctorForce': whether corrector should forcefully delete bad
poles; defaults to False;
- 'correctorTol': tolerance for corrector step; defaults to 0.,
i.e. no bad poles;
- 'correctorMaxIter': maximum number of corrector iterations;
defaults to 1.
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;
- - 'nNearestNeighbor': mumber of nearest neighbors considered in
- numerator if polybasis allows;
- 'interpRcond': tolerance for interpolation via numpy.polyfit;
- 'robustTol': tolerance for robust rational denominator
management;
- '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 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.
- nNearestNeighbor: Number of nearest neighbors considered in numerator
- if polybasis allows.
interpRcond: Tolerance for interpolation via numpy.polyfit.
robustTol: Tolerance for robust rational denominator management.
correctorForce: Whether corrector should forcefully delete bad poles.
correctorTol: Tolerance for corrector step.
correctorMaxIter: Maximum number of corrector iterations.
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()
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)
else:
ev, eV = self.findeveVGExplicit(
self.samplingEngine.samples(idxSamplesEff), invD)
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])
else:
q.coeffs = eV[:, 0].reshape([self.N + 1] * self.npar)
vbMng(self, "DEL", "Done computing denominator.", 7)
return q, fitinv
def _setupNumerator(self):
"""Compute rational numerator."""
RROMPyAssert(self._mode, message = "Cannot setup numerator.")
vbMng(self, "INIT", "Starting computation of numerator.", 7)
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.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:
pParRest = [self.radialDirectionalWeights] + pParRest
- pParRest[-1]["nNearestNeighbor"] = self.nNearestNeighbor
p = RBI() if self.polybasis in rbpb else MLSI()
if self.polybasis in ppb + rbpb:
pParRest += [{"rcond": self.interpRcond}]
wellCond, msg = p.setupByInterpolation(self._musUniqueCN,
Qevaldiag[: Seff, : Seff],
self.M, self.polybasis,
*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, TEGen = totalDegreeN, pvTP
else:
cfun, TEGen = fullDegreeN, pvP
E = max(self.M, self.N)
while E >= 0:
Seff = cfun(E, self.npar)
TEGenPar = [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 = TEGen(self._musUniqueCN, Eeff, *TEGenPar)
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 = TEGen(self._musUniqueCN, Neff, *TEGenPar)
for k in range(len(invD)): invD[k] = dot(invD[k], TN)
return invD, fitinv
diff --git a/rrompy/reduction_methods/standard/trained_model/trained_model_rational_mls.py b/rrompy/reduction_methods/standard/trained_model/trained_model_rational_mls.py
index 329c5b5..d2fcaa3 100644
--- a/rrompy/reduction_methods/standard/trained_model/trained_model_rational_mls.py
+++ b/rrompy/reduction_methods/standard/trained_model/trained_model_rational_mls.py
@@ -1,191 +1,189 @@
# 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 .trained_model_rational import TrainedModelRational
from rrompy.utilities.base.types import Np1D, paramVal, paramList, sampList
from rrompy.utilities.base import verbosityManager as vbMng
from rrompy.utilities.poly_fitting.moving_least_squares import mlsweights
from rrompy.utilities.poly_fitting.polynomial import (
PolynomialInterpolator as PI)
from rrompy.utilities.numerical import customPInv
from rrompy.utilities.numerical.degree import degreeTotalToFull
from rrompy.parameter import checkParameterList
from rrompy.sampling import emptySampleList
from rrompy.utilities.exception_manager import RROMPyException, RROMPyWarning
__all__ = ['TrainedModelRationalMLS']
class TrainedModelRationalMLS(TrainedModelRational):
"""
ROM approximant evaluation for rational moving least squares approximant.
Attributes:
Data: dictionary with all that can be pickled.
"""
def reset(self):
super().reset()
self.lastSetupMu = None
def assembleReducedModel(self, mu:paramVal):
if not (hasattr(self.data, "lastSetupMu")
and self.data.lastSetupMu == mu):
vbMng(self, "INIT", "Assembling reduced model for mu = {}."\
.format(mu), 17)
vbMng(self, "INIT", "Starting computation of denominator.", 35)
muC = self.centerNormalize(mu)
muSC = self.centerNormalize(self.data.mus)
wQ = mlsweights(muC, muSC, self.data.radialBasisDen,
- directionalWeights = self.data.radialWeightsDen,
- nNearestNeighbor = self.data.nNearestNeighborDen)
+ directionalWeights = self.data.radialWeightsDen)
if self.data.N > self.data.M:
PQVan = self.data.QVan
else:
PQVan = self.data.PVan
VQAdjW = PQVan.conj().T * wQ
VQAdjWVQ = VQAdjW.dot(PQVan)
interpPseudoInverse, info = customPInv(VQAdjWVQ, full = True,
rcond = self.data.interpRcond)
interpPseudoInverse = interpPseudoInverse.dot(VQAdjW).dot(
self.data.QBlocks)
if info[0] < interpPseudoInverse.shape[-1]:
q = np.zeros(interpPseudoInverse.shape[-1], dtype = np.complex)
q[0] = 1.
else:
halfGram = interpPseudoInverse[self.data.domQIdxs]
if self.data.POD:
Rstack = halfGram.reshape(-1, halfGram.shape[-1])
vbMng(self, "INIT",
"Solving svd for square root of gramian matrix.", 67)
try:
_, s, eV = np.linalg.svd(Rstack, full_matrices = False)
except np.linalg.LinAlgError as e:
raise RROMPyException(e)
condN = s[0] / s[-1]
q = eV[-1, :].T.conj()
vbMng(self, "MAIN",
("Solved svd problem of size {} x {} with condition "
"number {:.4e}.").format(*Rstack.shape, condN), 55)
vbMng(self, "DEL", "Done solving svd.", 67)
else:
RRstack = np.tensordot(self.trainedModel.gramian, halfGram,
1).reshape(-1, halfGram.shape[-1])
RLstack = halfGram.reshape(-1, halfGram.shape[-1])
gram = RLstack.T.conj().dot(RRstack)
vbMng(self, "INIT",
"Solving eigenvalue problem for gramian matrix.", 67)
try:
ev, eV = np.linalg.eigh(gram)
except np.linalg.LinAlgError as e:
raise RROMPyException(e)
condN = ev[-1] / ev[0]
q = eV[:, 0]
vbMng(self, "MAIN",
("Solved eigenvalue problem of size {} with "
"condition number {:.4e}.").format(gram.shape[0],
condN), 55)
vbMng(self, "DEL", "Done solving eigenvalue problem.", 67)
self.data.Q = PI()
self.data.Q.npar = self.npar
self.data.Q.polybasis = self.data.polybasis
if self.data.polydegreetype == "TOTAL":
self.data.Q.coeffs = degreeTotalToFull(
(self.data.N + 1,) * self.npar,
self.npar, q)
else:
self.data.Q.coeffs = q.reshape((self.data.N + 1,) * self.npar)
vbMng(self, "DEL", "Done computing denominator.", 35)
vbMng(self, "INIT", "Starting computation of numerator.", 35)
self.data.P = PI()
self.data.P.npar = self.npar
self.data.P.polybasis = self.data.polybasis
wP = mlsweights(muC, muSC, self.data.radialBasis,
- directionalWeights = self.data.radialWeights,
- nNearestNeighbor = self.data.nNearestNeighbor)
+ directionalWeights = self.data.radialWeights)
VAdjW = self.data.PVan.conj().T * wP
VAdjWV = VAdjW.dot(self.data.PVan)
interpPPseudoInverse = customPInv(VAdjWV, self.data.interpRcond)
Pcoeffs = np.tensordot(interpPPseudoInverse.dot(VAdjW),
self.data.QBlocks.dot(q), ([1], [1]))
if self.data.polydegreetype == "TOTAL":
self.data.P.coeffs = degreeTotalToFull(
(self.data.M + 1,) * self.npar
+ (self.data.QBlocks.shape[0],),
self.npar, Pcoeffs)
else:
self.data.P.coeffs = Pcoeffs.reshape(
(self.data.M + 1,) * self.npar
+ (self.data.QBlocks.shape[0],))
vbMng(self, "DEL", "Done computing numerator.", 35)
vbMng(self, "DEL", "Done assembling reduced model.", 17)
self.data.lastSetupMu = mu
def getApproxReduced(self, mu : paramList = []) -> sampList:
"""
Evaluate reduced representation of approximant at arbitrary parameter.
Args:
mu: Target parameter.
"""
mu = checkParameterList(mu, self.data.npar)[0]
if (not hasattr(self, "lastSolvedApproxReduced")
or self.lastSolvedApproxReduced != mu):
vbMng(self, "INIT",
"Evaluating approximant at mu = {}.".format(mu), 12)
self.uApproxReduced = emptySampleList()
for i in range(len(mu)):
self.assembleReducedModel(mu[i])
vbMng(self, "INIT",
"Solving reduced model for mu = {}.".format(mu[i]), 15)
Qv = self.getQVal(mu[i])
if Qv == 0.:
RROMPyWarning(("Adjusting approximation to avoid division "
"by numerically zero denominator."))
Qv = np.finfo(np.complex).eps / (1. + self.data.Q.deg[0])
uAppR = self.getPVal(mu[i]) / Qv
if i == 0:
self.uApproxReduced.reset((uAppR.shape[0], len(mu)),
dtype = uAppR.dtype)
self.uApproxReduced[i] = uAppR
vbMng(self, "DEL", "Done solving reduced model.", 15)
vbMng(self, "DEL", "Done evaluating approximant.", 12)
self.lastSolvedApproxReduced = mu
return self.uApproxReduced
def getPoles(self, *args, mu : paramVal = None, **kwargs) -> Np1D:
"""
Obtain approximant poles.
Returns:
Numpy complex vector of poles.
"""
if mu is None: mu = self.data.mu0
self.assembleReducedModel(mu)
return super().getPoles(*args, **kwargs)
def getResidues(self, *args, mu : paramVal = None, **kwargs) -> Np1D:
"""
Obtain approximant residues.
Returns:
Numpy matrix with residues as columns.
"""
if mu is None: mu = self.data.mu0
self.assembleReducedModel(mu)
return super().getResidues(*args, **kwargs)
diff --git a/rrompy/utilities/poly_fitting/moving_least_squares/__init__.py b/rrompy/utilities/poly_fitting/moving_least_squares/__init__.py
index e230f15..c126e42 100644
--- a/rrompy/utilities/poly_fitting/moving_least_squares/__init__.py
+++ b/rrompy/utilities/poly_fitting/moving_least_squares/__init__.py
@@ -1,41 +1,40 @@
# 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 .kernel import (radialGaussian, thinPlateSpline, multiQuadric,
- localWendland, nearestNeighbor)
+ localWendland)
from .base import mlsbases, polybases, polyfitname, polydomcoeff
from .vander import mlsweights, polyvander, polyvanderTotal
from .moving_least_squares_interpolator import MovingLeastSquaresInterpolator
__all__ = [
'radialGaussian',
'thinPlateSpline',
'multiQuadric',
'localWendland',
- 'nearestNeighbor',
'mlsbases',
'polybases',
'polyfitname',
'polydomcoeff',
'mlsweights',
'polyvander',
'polyvanderTotal',
'MovingLeastSquaresInterpolator'
]
diff --git a/rrompy/utilities/poly_fitting/moving_least_squares/kernel.py b/rrompy/utilities/poly_fitting/moving_least_squares/kernel.py
index 1708271..207a8f1 100644
--- a/rrompy/utilities/poly_fitting/moving_least_squares/kernel.py
+++ b/rrompy/utilities/poly_fitting/moving_least_squares/kernel.py
@@ -1,24 +1,24 @@
# 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 rrompy.utilities.poly_fitting.radial_basis.kernel import (radialGaussian,
- thinPlateSpline, multiQuadric, localWendland, nearestNeighbor)
+ thinPlateSpline, multiQuadric, localWendland)
__all__ = ['radialGaussian', 'thinPlateSpline', 'multiQuadric',
- 'localWendland', 'nearestNeighbor']
+ 'localWendland']
diff --git a/rrompy/utilities/poly_fitting/moving_least_squares/moving_least_squares_interpolator.py b/rrompy/utilities/poly_fitting/moving_least_squares/moving_least_squares_interpolator.py
index a333338..7bd4ef8 100644
--- a/rrompy/utilities/poly_fitting/moving_least_squares/moving_least_squares_interpolator.py
+++ b/rrompy/utilities/poly_fitting/moving_least_squares/moving_least_squares_interpolator.py
@@ -1,142 +1,139 @@
# 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 (List, ListAny, DictAny, Np1D, Np2D,
paramList)
from rrompy.utilities.numerical import customPInv, dot
from .vander import mlsweights
from rrompy.utilities.poly_fitting.interpolator import GenericInterpolator
from rrompy.utilities.poly_fitting.polynomial.vander import (polyvander as pv,
polyvanderTotal as pvT)
from rrompy.utilities.exception_manager import RROMPyException, RROMPyAssert
from rrompy.parameter import checkParameterList
__all__ = ['MovingLeastSquaresInterpolator']
class MovingLeastSquaresInterpolator(GenericInterpolator):
def __init__(self, other = None):
if other is None: return
self.support = other.support
self.localProjector = other.localProjector
self.localVanders = other.localVanders
self.suppValues = other.suppValues
self.directionalWeights = other.directionalWeights
self.degree = other.degree
self.npar = other.npar
self.radialbasis = other.radialbasis
self.polybasis = other.polybasis
self.evalParams = other.evalParams
self.totalDegree = other.totalDegree
@property
def shape(self):
sh = self.suppValues.shape[1 :] if self.suppValues.ndim > 1 else 1
return sh
@property
def deg(self):
return self.degree
def __call__(self, mu:paramList, der : List[int] = None,
scl : Np1D = None):
if der is not None and np.sum(der) > 0:
raise RROMPyException(("Cannot take derivatives of moving least "
"squares function."))
mu = checkParameterList(mu, self.npar)[0]
sh = self.shape
if sh == 1: sh = tuple([])
values = np.empty((len(mu),) + sh, dtype = np.complex)
for i, m in enumerate(mu):
weights = mlsweights(m, self.support, self.radialbasis,
- directionalWeights = self.directionalWeights,
- nNearestNeighbor = self.evalParams["nNearestNeighbor"])
+ directionalWeights = self.directionalWeights)
weights /= np.linalg.norm(weights)
vanderLS = np.sum(self.localVanders * weights, axis = 2)
RHSLS = dot(self.localProjector * weights, self.suppValues)
if self.totalDegree:
vanderEval = pvT(m, self.deg[0], self.polybasis,
**self.evalParams)
else:
vanderEval = pv(m, self.deg, self.polybasis, **self.evalParams)
vanderEval = vanderEval.flatten()
values[i] = dot(vanderEval, dot(customPInv(vanderLS), RHSLS))
return values
def __copy__(self):
return MovingLeastSquaresInterpolator(self)
def __deepcopy__(self, memo):
other = MovingLeastSquaresInterpolator()
(other.support, other.localProjector, other.localVanders,
other.suppValues, other.directionalWeights, other.degree, other.npar,
other.radialbasis, other.polybasis, other.evalParams,
other.totalDegree) = copy(
(self.support, self.localProjector, self.localVanders,
self.suppValues, self.directionalWeights, self.degree,
self.npar, self.radialbasis, self.polybasis,
self.evalParams, self.totalDegree), memo)
return other
def postmultiplyTensorize(self, A:Np2D):
RROMPyAssert(A.shape[0], self.shape[-1], "Shape of output")
self.suppValues = self.suppValues.dot(A)
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)] + [(l, r) for l, r in zip(nleft, nright)]
self.suppValues = np.pad(self.suppValues, padwidth, "constant",
constant_values = (0., 0.))
def setupByInterpolation(self, support:paramList, values:ListAny,
deg:int, polybasis : str = "MONOMIAL_GAUSSIAN",
directionalWeights : Np1D = None,
totalDegree : bool = True,
vanderCoeffs : DictAny = {}):
support = checkParameterList(support)[0]
self.support = copy(support)
if "reorder" in vanderCoeffs.keys():
self.support = self.support[vanderCoeffs["reorder"]]
- if "nNearestNeighbor" not in vanderCoeffs.keys():
- vanderCoeffs["nNearestNeighbor"] = -1
self.npar = support.shape[1]
if directionalWeights is None:
directionalWeights = np.ones(self.npar)
self.directionalWeights = directionalWeights
self.polybasis, self.radialbasis, _ = polybasis.split("_")
self.totalDegree = totalDegree
self.evalParams = vanderCoeffs
if totalDegree:
vander = pvT(support, deg, self.polybasis, **vanderCoeffs)
if not hasattr(deg, "__len__"): deg = [deg] * self.npar
else:
if not hasattr(deg, "__len__"): deg = [deg] * self.npar
vander = pv(support, deg, self.polybasis, **vanderCoeffs)
self.degree = deg
self.localProjector = vander.T.conj()
self.localVanders = np.array([np.outer(van, van.conj()) \
for van in vander])
self.localVanders = np.swapaxes(self.localVanders, 0, 2)
self.suppValues = np.array(values)
diff --git a/rrompy/utilities/poly_fitting/moving_least_squares/vander.py b/rrompy/utilities/poly_fitting/moving_least_squares/vander.py
index 141023b..f378f57 100644
--- a/rrompy/utilities/poly_fitting/moving_least_squares/vander.py
+++ b/rrompy/utilities/poly_fitting/moving_least_squares/vander.py
@@ -1,96 +1,86 @@
# 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 .kernel import (radialGaussian, thinPlateSpline, multiQuadric,
- localWendland, nearestNeighbor)
+ localWendland)
from rrompy.utilities.poly_fitting.polynomial.vander import (polyvander as pvP,
polyvanderTotal as pvTP)
from rrompy.utilities.base.types import (Np1D, Np2D, Tuple, List, paramVal,
paramList)
from rrompy.parameter import checkParameter, checkParameterList
from rrompy.utilities.exception_manager import RROMPyException, RROMPyAssert
__all__ = ['mlsweights', 'polyvander', 'polyvanderTotal']
def mlsweights(x:paramVal, xSupp:paramList, basis:str,
- reorder : List[int] = None, directionalWeights : Np1D = None,
- nNearestNeighbor : int = -1) -> Np2D:
+ reorder : List[int] = None,
+ directionalWeights : Np1D = None) -> Np2D:
"""Compute moving least squares weight vector."""
x = checkParameter(x)
xSupp = checkParameterList(xSupp)[0]
x = x.data
xSupp = xSupp.data
if directionalWeights is None:
directionalWeights = np.ones(x.shape[1])
elif not hasattr(directionalWeights, "__len__"):
directionalWeights = directionalWeights * np.ones(x.shape[1])
RROMPyAssert(len(directionalWeights), x.shape[1],
"Number of directional weights")
try:
- radialkernel = {
- "GAUSSIAN" : radialGaussian, "THINPLATE" : thinPlateSpline,
- "MULTIQUADRIC" : multiQuadric, "WENDLAND" : localWendland,
- "NEARESTNEIGHBOR" : nearestNeighbor}[basis.upper()]
+ radialkernel = {"GAUSSIAN" : radialGaussian,
+ "THINPLATE" : thinPlateSpline,
+ "MULTIQUADRIC" : multiQuadric,
+ "WENDLAND" : localWendland}[basis.upper()]
except:
raise RROMPyException("Radial basis not recognized.")
if reorder is not None: xSupp = xSupp[reorder]
muDiff = (xSupp.data - x) * directionalWeights
r2 = np.sum(np.abs(muDiff) ** 2., axis = 1).reshape(1, -1)
- if basis.upper() == "NEARESTNEIGHBOR":
- if nNearestNeighbor > 0 and nNearestNeighbor < len(xSupp):
- cutoffValue = np.partition(r2, nNearestNeighbor - 1)[0,
- nNearestNeighbor - 1]
- r2 /= cutoffValue
- else:
- r2[0, :] = 1. * (nNearestNeighbor == 0)
return radialkernel(r2)[0]
def polyvander(x:paramVal, xSupp:paramList, degs:List[int], basis:str,
derIdxs : List[List[List[int]]] = None,
reorder : List[int] = None, directionalWeights : Np1D = None,
- scl : Np1D = None,
- nNearestNeighbor : int = -1) -> Tuple[Np2D, Np2D]:
+ scl : Np1D = None) -> Tuple[Np2D, Np2D]:
"""
Compute full Hermite-Vandermonde matrix with specified derivative
directions.
"""
basisp, basisr, _ = basis.split("_")
- Weights = mlsweights(x, xSupp, basisr, reorder, directionalWeights,
- nNearestNeighbor)
+ Weights = mlsweights(x, xSupp, basisr, reorder, directionalWeights)
VanP = pvP(xSupp, degs, basisp, derIdxs = derIdxs, reorder = reorder,
scl = scl)
RHP = VanP.T.conj() * Weights
return RHP.dot(VanP), RHP
def polyvanderTotal(x:paramList, xSupp:paramList, deg:int, basis:str,
derIdxs : List[List[List[int]]] = None,
reorder : List[int] = None,
- directionalWeights : Np1D = None, scl : Np1D = None,
- nNearestNeighbor : int = -1) -> Tuple[Np2D, Np2D]:
+ directionalWeights : Np1D = None,
+ scl : Np1D = None) -> Tuple[Np2D, Np2D]:
"""
Compute full total degree Hermite-Vandermonde matrix with specified
derivative directions.
"""
basisp, basisr, _ = basis.split("_")
- Weights = mlsweights(x, xSupp, basisr, reorder, directionalWeights,
- nNearestNeighbor)
+ Weights = mlsweights(x, xSupp, basisr, reorder, directionalWeights)
VanP = pvTP(x, deg, basisp, derIdxs = derIdxs,
reorder = reorder, scl = scl)
RHP = VanP.T.conj() * Weights
return RHP.dot(VanP), RHP
diff --git a/rrompy/utilities/poly_fitting/radial_basis/__init__.py b/rrompy/utilities/poly_fitting/radial_basis/__init__.py
index 8830a91..98a445b 100644
--- a/rrompy/utilities/poly_fitting/radial_basis/__init__.py
+++ b/rrompy/utilities/poly_fitting/radial_basis/__init__.py
@@ -1,43 +1,42 @@
# 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 .kernel import (radialGaussian, thinPlateSpline, multiQuadric,
- localWendland, nearestNeighbor)
+ localWendland)
from .base import rbbases, polybases, polyfitname, polydomcoeff
from .val import polyval
from .vander import rbvander, polyvander, polyvanderTotal
from .radial_basis_interpolator import RadialBasisInterpolator
__all__ = [
'radialGaussian',
'thinPlateSpline',
'multiQuadric',
'localWendland',
- 'nearestNeighbor',
'rbbases',
'polybases',
'polyfitname',
'polydomcoeff',
'polyval',
'rbvander',
'polyvander',
'polyvanderTotal',
'RadialBasisInterpolator'
]
diff --git a/rrompy/utilities/poly_fitting/radial_basis/base.py b/rrompy/utilities/poly_fitting/radial_basis/base.py
index 4fc5818..2423bf0 100644
--- a/rrompy/utilities/poly_fitting/radial_basis/base.py
+++ b/rrompy/utilities/poly_fitting/radial_basis/base.py
@@ -1,45 +1,43 @@
# 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 itertools import product
from rrompy.utilities.exception_manager import RROMPyException
from rrompy.utilities.poly_fitting.polynomial.base import (polybases as pbP,
polyfitname as pfnP,
polydomcoeff as polydomcoeffB)
__all__ = ['rbbases', 'polybases', 'polyfitname', 'polydomcoeff']
-rbbases = ["GAUSSIAN", "THINPLATE", "MULTIQUADRIC", "WENDLAND",
- "NEARESTNEIGHBOR"]
+rbbases = ["GAUSSIAN", "THINPLATE", "MULTIQUADRIC", "WENDLAND"]
polybases = [x + "_" + y for x, y in product(pbP, rbbases)]
def polyfitname(basis:str) -> str:
fitrnames = {"GAUSSIAN" : "gaussian", "THINPLATE" : "thinplate",
- "MULTIQUADRIC" : "multiquadric", "WENDLAND" : "wendland",
- "NEARESTNEIGHBOR" : "nearestneighbor"}
+ "MULTIQUADRIC" : "multiquadric", "WENDLAND" : "wendland"}
basisp, basisr = basis.split("_")
try:
return pfnP(basisp) + "_" + fitrnames[basisr]
except:
raise RROMPyException("Polynomial-radial basis combination not "
"recognized.")
def polydomcoeff(n:int, basis:str) -> float:
return polydomcoeffB(n, basis.split("_")[0])
diff --git a/rrompy/utilities/poly_fitting/radial_basis/kernel.py b/rrompy/utilities/poly_fitting/radial_basis/kernel.py
index e2af2b6..e349e58 100644
--- a/rrompy/utilities/poly_fitting/radial_basis/kernel.py
+++ b/rrompy/utilities/poly_fitting/radial_basis/kernel.py
@@ -1,46 +1,42 @@
# 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.utilities.base.types import Np1D
from rrompy.utilities.exception_manager import RROMPyAssert
__all__ = ['radialGaussian', 'thinPlateSpline', 'multiQuadric',
- 'localWendland', 'nearestNeighbor']
+ 'localWendland']
def radialGaussian(r2:Np1D, der : int = 0) -> Np1D:
RROMPyAssert(der, 0, "Radial basis derivative")
return np.exp(- .5 * r2)
def thinPlateSpline(r2:Np1D, der : int = 0) -> Np1D:
RROMPyAssert(der, 0, "Radial basis derivative")
return .5 * r2 * np.log(np.finfo(float).eps + r2)
def multiQuadric(r2:Np1D, der : int = 0) -> Np1D:
RROMPyAssert(der, 0, "Radial basis derivative")
return np.power(r2 + 1., -.5)
def localWendland(r2:Np1D, der : int = 0) -> Np1D:
RROMPyAssert(der, 0, "Radial basis derivative")
rm1 = 1. - r2 ** .5
rm1[rm1 <= 0.] = 0.
return rm1 ** 4. * (5. - 4. * rm1)
-
-def nearestNeighbor(r2:Np1D, der : int = 0) -> Np1D:
- RROMPyAssert(der, 0, "Radial basis derivative")
- return 1. * (r2 <= 1.)
diff --git a/rrompy/utilities/poly_fitting/radial_basis/radial_basis_interpolator.py b/rrompy/utilities/poly_fitting/radial_basis/radial_basis_interpolator.py
index 9b2ba80..bd84a36 100644
--- a/rrompy/utilities/poly_fitting/radial_basis/radial_basis_interpolator.py
+++ b/rrompy/utilities/poly_fitting/radial_basis/radial_basis_interpolator.py
@@ -1,144 +1,137 @@
# 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 (List, ListAny, DictAny, Np1D, Np2D,
paramList)
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 .vander import polyvander as pv, polyvanderTotal as pvT
from rrompy.utilities.numerical import dot
from rrompy.utilities.numerical.degree import degreeTotalToFull
from rrompy.utilities.exception_manager import RROMPyException, RROMPyAssert
from rrompy.parameter import checkParameterList
__all__ = ['RadialBasisInterpolator']
class RadialBasisInterpolator(GenericInterpolator):
def __init__(self, other = None):
if other is None: return
self.support = other.support
self.coeffsGlobal = other.coeffsGlobal
self.coeffsLocal = other.coeffsLocal
self.directionalWeights = other.directionalWeights
self.npar = other.npar
self.polybasis = other.polybasis
- self.nNearestNeighbor = other.nNearestNeighbor
@property
def shape(self):
sh = self.coeffsLocal.shape[1 :] if self.coeffsLocal.ndim > 1 else 1
return sh
@property
def deg(self):
return [x - 1 for x in self.coeffsGlobal.shape[: self.npar]]
def __call__(self, mu:paramList, der : List[int] = None,
scl : Np1D = None):
if der is not None and np.sum(der) > 0:
raise RROMPyException(("Cannot take derivatives of radial basis "
"function."))
return polyval(mu, self.coeffsGlobal, self.coeffsLocal, self.support,
- self.directionalWeights, self.polybasis,
- self.nNearestNeighbor)
+ self.directionalWeights, self.polybasis)
def __copy__(self):
return RadialBasisInterpolator(self)
def __deepcopy__(self, memo):
other = RadialBasisInterpolator()
(other.support, other.coeffsGlobal, other.coeffsLocal,
- other.directionalWeights, other.npar, other.polybasis,
- other.nNearestNeighbor) = copy(
- (self.support, self.coeffsGlobal, self.coeffsLocal,
- self.directionalWeights, self.npar, self.polybasis,
- self.nNearestNeighbor), memo)
+ other.directionalWeights, other.npar, other.polybasis) = copy(
+ (self.support, self.coeffsGlobal,
+ self.coeffsLocal, self.directionalWeights,
+ self.npar, self.polybasis), memo)
return other
def postmultiplyTensorize(self, A:Np2D):
RROMPyAssert(A.shape[0], self.shape[-1], "Shape of output")
self.coeffsLocal = dot(self.coeffsLocal, A)
self.coeffsGlobal = dot(self.coeffsGlobal, A)
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)] + [(l, r) for l, r in zip(nleft, nright)]
self.coeffsLocal = np.pad(self.coeffsLocal, padwidth, "constant",
constant_values = (0., 0.))
padwidth = [(0, 0)] * (self.npar - 1) + padwidth
self.coeffsGlobal = np.pad(self.coeffsGlobal, padwidth, "constant",
constant_values = (0., 0.))
def setupByInterpolation(self, support:paramList, values:ListAny,
deg:int, polybasis : str = "MONOMIAL_GAUSSIAN",
directionalWeights : Np1D = None,
verbose : bool = True, totalDegree : bool = True,
vanderCoeffs : DictAny = {},
fitCoeffs : DictAny = {}):
support = checkParameterList(support)[0]
self.support = copy(support)
if "reorder" in vanderCoeffs.keys():
self.support = self.support[vanderCoeffs["reorder"]]
- if "nNearestNeighbor" in vanderCoeffs.keys():
- self.nNearestNeighbor = vanderCoeffs["nNearestNeighbor"]
- else:
- self.nNearestNeighbor = -1
self.npar = support.shape[1]
if directionalWeights is None:
directionalWeights = np.ones(self.npar)
self.directionalWeights = directionalWeights
self.polybasis = polybasis
if totalDegree:
vander = pvT(support, deg, basis = polybasis,
directionalWeights = self.directionalWeights,
**vanderCoeffs)
else:
if not hasattr(deg, "__len__"): deg = [deg] * self.npar
vander = pv(support, deg, basis = polybasis,
directionalWeights = self.directionalWeights,
**vanderCoeffs)
outDim = values.shape[1:]
values = values.reshape(values.shape[0], -1)
values = np.pad(values, ((0, len(vander) - len(values)), (0, 0)),
"constant")
fitOut = customFit(vander, values, full = True, **fitCoeffs)
P = fitOut[0][len(support) :]
if verbose:
msg = ("Fitting {}+{} samples with degree {} through {}... "
"Conditioning of LS system: {:.4e}.").format(
len(support), len(vander) - len(support),
deg, polyfitname(self.polybasis),
fitOut[1][2][0] / fitOut[1][2][-1])
else: msg = None
self.coeffsLocal = fitOut[0][: len(support)]
if totalDegree:
self.coeffsGlobal = degreeTotalToFull(tuple([deg + 1] * self.npar)
+ outDim, self.npar, P)
else:
self.coeffsGlobal = P.reshape(tuple([d + 1 for d in deg]) + outDim)
return fitOut[1][1] == vander.shape[1], msg
diff --git a/rrompy/utilities/poly_fitting/radial_basis/val.py b/rrompy/utilities/poly_fitting/radial_basis/val.py
index 25d9ac5..a2eb877 100644
--- a/rrompy/utilities/poly_fitting/radial_basis/val.py
+++ b/rrompy/utilities/poly_fitting/radial_basis/val.py
@@ -1,63 +1,54 @@
# 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 .kernel import (radialGaussian, thinPlateSpline, multiQuadric,
- localWendland, nearestNeighbor)
+ localWendland)
from rrompy.utilities.poly_fitting.polynomial.val import polyval as pvP
from rrompy.utilities.base.types import Np1D, Np2D, paramList
from rrompy.parameter import checkParameterList
from rrompy.utilities.exception_manager import RROMPyException
__all__ = ['polyval']
def polyval(x:paramList, cG:Np2D, cL:Np2D, supportPoints:paramList,
- directionalWeights:Np1D, basis:str,
- nNearestNeighbor : int = -1) -> Np2D:
+ directionalWeights:Np1D, basis:str) -> Np2D:
x = checkParameterList(x)[0]
basisp, basisr = basis.split("_")
c = pvP(x, cG, basisp)
try:
- radialvalbase = {
- "GAUSSIAN" : radialGaussian, "THINPLATE" : thinPlateSpline,
- "MULTIQUADRIC" : multiQuadric, "WENDLAND" : localWendland,
- "NEARESTNEIGHBOR" : nearestNeighbor}[basisr.upper()]
+ radialvalbase = {"GAUSSIAN" : radialGaussian,
+ "THINPLATE" : thinPlateSpline,
+ "MULTIQUADRIC" : multiQuadric,
+ "WENDLAND" : localWendland}[basisr.upper()]
except:
raise RROMPyException("Radial basis not recognized.")
supportPoints = checkParameterList(supportPoints)[0]
- isnearestneighbor = basisr.upper() == "NEARESTNEIGHBOR"
csh = copy(c.shape)
if len(csh) == 1: c = c.reshape(1, -1)
for j in range(len(x)):
muDiff = (supportPoints.data - x[j]) * directionalWeights
r2j = np.sum(np.abs(muDiff) ** 2., axis = 1).reshape(1, -1)
- if isnearestneighbor:
- if nNearestNeighbor > 0 and nNearestNeighbor < len(supportPoints):
- cutoffValue = np.partition(r2j, nNearestNeighbor - 1)[0,
- nNearestNeighbor - 1]
- r2j /= cutoffValue
- else:
- r2j[0, :] = 1. * (nNearestNeighbor == 0)
val = radialvalbase(r2j).dot(cL)
try:
c[..., j] += val
except:
c[..., j] += val.flatten()
if len(csh) == 1: c = c.flatten()
return c
diff --git a/rrompy/utilities/poly_fitting/radial_basis/vander.py b/rrompy/utilities/poly_fitting/radial_basis/vander.py
index ba246de..773a5b3 100644
--- a/rrompy/utilities/poly_fitting/radial_basis/vander.py
+++ b/rrompy/utilities/poly_fitting/radial_basis/vander.py
@@ -1,108 +1,97 @@
# 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 .kernel import (radialGaussian, thinPlateSpline, multiQuadric,
- localWendland, nearestNeighbor)
+ localWendland)
from rrompy.utilities.poly_fitting.polynomial.vander import (polyvander as pvP,
polyvanderTotal as pvTP)
from rrompy.utilities.base.types import Np1D, Np2D, List, paramList
from rrompy.parameter import checkParameterList
from rrompy.utilities.exception_manager import RROMPyException, RROMPyAssert
__all__ = ['rbvander', 'polyvander', 'polyvanderTotal']
def rbvander(x:paramList, basis:str, reorder : List[int] = None,
- directionalWeights : Np1D = None,
- nNearestNeighbor : int = -1) -> Np2D:
+ directionalWeights : Np1D = None) -> Np2D:
"""Compute radial-basis-Vandermonde matrix."""
x = checkParameterList(x)[0]
x_un = x.unique()
nx = len(x)
if len(x_un) < nx:
raise RROMPyException("Sample points must be distinct.")
del x_un
x = x.data
if directionalWeights is None:
directionalWeights = np.ones(x.shape[1])
elif not hasattr(directionalWeights, "__len__"):
directionalWeights = directionalWeights * np.ones(x.shape[1])
RROMPyAssert(len(directionalWeights), x.shape[1],
"Number of directional weights")
try:
- radialkernel = {
- "GAUSSIAN" : radialGaussian, "THINPLATE" : thinPlateSpline,
- "MULTIQUADRIC" : multiQuadric, "WENDLAND" : localWendland,
- "NEARESTNEIGHBOR" : nearestNeighbor}[basis.upper()]
+ radialkernel = {"GAUSSIAN" : radialGaussian,
+ "THINPLATE" : thinPlateSpline,
+ "MULTIQUADRIC" : multiQuadric,
+ "WENDLAND" : localWendland}[basis.upper()]
except:
raise RROMPyException("Radial basis not recognized.")
- isnearestneighbor = basis.upper() == "NEARESTNEIGHBOR"
Van = np.zeros((nx, nx))
if reorder is not None: x = x[reorder]
for j in range(nx):
muDiff = (x - x[j]) * directionalWeights
r2j = np.sum(np.abs(muDiff) ** 2., axis = 1).reshape(1, -1)
- if isnearestneighbor:
- if nNearestNeighbor > 0 and nNearestNeighbor < len(x):
- cutoffValue = np.partition(r2j, nNearestNeighbor - 1)[0,
- nNearestNeighbor - 1]
- r2j /= cutoffValue
- else:
- r2j[0, :] = 1. * (nNearestNeighbor == 0)
Van[j] = radialkernel(r2j)
return Van
def polyvander(x:paramList, degs:List[int], basis:str,
derIdxs : List[List[List[int]]] = None,
reorder : List[int] = None, directionalWeights : Np1D = None,
- scl : Np1D = None, nNearestNeighbor : int = -1) -> Np2D:
+ scl : Np1D = None) -> Np2D:
"""
Compute full Hermite-Vandermonde matrix with specified derivative
directions.
"""
if derIdxs is not None and np.sum(np.sum(derIdxs)) > 0:
raise RROMPyException(("Cannot take derivatives of radial basis "
"function."))
basisp, basisr = basis.split("_")
VanR = rbvander(x, basisr, reorder = reorder,
- directionalWeights = directionalWeights,
- nNearestNeighbor = nNearestNeighbor)
+ directionalWeights = directionalWeights)
VanP = pvP(x, degs, basisp, derIdxs = derIdxs, reorder = reorder,
scl = scl)
return np.block([[VanR, VanP],
[VanP.T.conj(), np.zeros(tuple([VanP.shape[1]] * 2))]])
def polyvanderTotal(x:paramList, deg:int, basis:str,
derIdxs : List[List[List[int]]] = None,
reorder : List[int] = None,
- directionalWeights : Np1D = None, scl : Np1D = None,
- nNearestNeighbor : int = -1) -> Np2D:
+ directionalWeights : Np1D = None,
+ scl : Np1D = None) -> Np2D:
"""
Compute full total degree Hermite-Vandermonde matrix with specified
derivative directions.
"""
if derIdxs is not None and np.sum(np.sum(derIdxs)) > 0:
raise RROMPyException(("Cannot take derivatives of radial basis "
"function."))
basisp, basisr = basis.split("_")
VanR = rbvander(x, basisr, reorder = reorder,
- directionalWeights = directionalWeights,
- nNearestNeighbor = nNearestNeighbor)
+ directionalWeights = directionalWeights)
VanP = pvTP(x, deg, basisp, derIdxs = derIdxs, reorder = reorder, scl = scl)
return np.block([[VanR, VanP],
[VanP.T.conj(), np.zeros(tuple([VanP.shape[1]] * 2))]])