Page MenuHomec4science
Files F60866976

generic_pivoted_greedy_approximant.py
No OneTemporary
Actions

File Metadata

Created
Fri, May 3, 01:06

generic_pivoted_greedy_approximant.py
View Options

# 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 <http://www.gnu.org/licenses/>.
#
from copy import deepcopy as copy
import numpy as np
from matplotlib import pyplot as plt
from rrompy.reduction_methods.pivoted import (GenericPivotedApproximant,
PODGlobal)
from rrompy.utilities.base.types import Np1D, Tuple, List, paramVal, paramList
from rrompy.utilities.base import verbosityManager as vbMng
from rrompy.utilities.numerical import (pointMatching, chordalMetricAdjusted,
checkCutOff)
from rrompy.utilities.exception_manager import (RROMPyException, RROMPyWarning,
RROMPyAssert)
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;
- 'matchingWeight': weight for pole matching optimization; defaults
to 1;
- 'cutOffTolerance': tolerance for ignoring parasitic poles;
defaults to np.inf;
- 'cutOffType': rule for tolerance computation for parasitic poles;
defaults to 'MAGNITUDE';
- 'matchingWeightError': weight for pole matching optimization in
error estimation; defaults to 0;
- 'cutOffToleranceError': tolerance for ignoring parasitic poles
in error estimation; defaults to np.inf;
- 'cutOffTypeError': rule for tolerance computation for parasitic
poles in error estimation; defaults to 'MAGNITUDE';
- '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;
- '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 0;
- '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;
- 'matchingWeight': weight for pole matching optimization;
- 'cutOffTolerance': tolerance for ignoring parasitic poles;
- 'matchingWeightError': weight for pole matching optimization in
error estimation;
- 'cutOffType': rule for tolerance computation for parasitic poles;
- 'cutOffToleranceError': tolerance for ignoring parasitic poles
in error estimation;
- 'cutOffTypeError': rule for tolerance computation for parasitic
poles in error estimation;
- '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.
matchingWeight: Weight for pole matching optimization.
cutOffTolerance: Tolerance for ignoring parasitic poles.
cutOffType: Rule for tolerance computation for parasitic poles.
matchingWeightError: Weight for pole matching optimization in error
estimation.
cutOffToleranceError: Tolerance for ignoring parasitic poles in error
estimation.
cutOffTypeError: Rule for tolerance computation for 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.
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.
muBoundsPivot: list of bounds for pivot parameter values.
muBoundsMarginal: list of bounds for marginal parameter values.
samplingEngine: Sampling engine.
uHF: High fidelity solution(s) with parameter(s) lastSolvedHF as
sampleList.
lastSolvedHF: Parameter(s) corresponding to last computed high fidelity
solution(s) as parameterList.
uApproxReduced: Reduced approximate solution(s) with parameter(s)
lastSolvedApprox as sampleList.
lastSolvedApproxReduced: Parameter(s) corresponding to last computed
reduced approximate solution(s) as parameterList.
uApprox: Approximate solution(s) with parameter(s) lastSolvedApprox as
sampleList.
lastSolvedApprox: Parameter(s) corresponding to last computed
approximate solution(s) as parameterList.
"""
def __init__(self, *args, **kwargs):
self._preInit()
from rrompy.parameter import localSparseGrid as SG
SGBase = SG([[0.], [1.]], "UNIFORM")
self._addParametersToList(["matchingWeightError",
"cutOffToleranceError", "cutOffTypeError",
"greedyTolMarginal", "maxIterMarginal"],
[0., np.inf, "MAGNITUDE", 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 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 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):
self._cutOffToleranceError = cutOffToleranceError
self._approxParameters["cutOffToleranceError"] = (
self.cutOffToleranceError)
@property
def cutOffTypeError(self):
"""Value of cutOffTypeError."""
return self._cutOffTypeError
@cutOffTypeError.setter
def cutOffTypeError(self, cutOffTypeError):
try:
cutOffTypeError = cutOffTypeError.upper().strip().replace(" ","")
if cutOffTypeError not in ["MAGNITUDE", "POTENTIAL"]:
raise RROMPyException(("Prescribed cutOffTypeError not "
"recognized."))
self._cutOffTypeError = cutOffTypeError
except:
RROMPyWarning(("Prescribed cutOffTypeError not recognized. "
"Overriding to 'MAGNITUDE'."))
self._cutOffTypeError = "MAGNITUDE"
self._approxParameters["cutOffTypeError"] = self.cutOffTypeError
@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()
self._figHandle = None
if hasattr(self, "samplingEngine") and self.samplingEngine is not None:
self.samplingEngine.resetHistory()
def errorEstimator(self, return_max : bool = False) -> Np1D:
self._finalizeMarginalization()
vbMng(self.trainedModel, "INIT",
"Evaluating error estimator at mu = {}.".format(
self.trainedModel.data.musMarginal), 10)
err = np.zeros(len(self.trainedModel.data.musMarginal))
if len(err) <= 1: err[:] = np.inf
else:
_musMExcl = None
self.verbosity -= 20
for j in range(len(err)):
jEff = j - (j > 0)
muTest = self.trainedModel.data.musMarginal[jEff]
polesEx = self.trainedModel.data.HIs[jEff].poles
idxExEff = checkCutOff(polesEx, rtype = self.cutOffTypeError,
tol = self.cutOffToleranceError)
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, self.approx_state)
self._musMUniqueCN = None
self.trainedModel.data.marginalInterp = (
self._setupMarginalInterp())
polesAp = self.trainedModel.interpolateMarginalPoles(muTest)
if self.matchingWeightError != 0:
resAp = self.trainedModel.interpolateMarginalCoeffs(
muTest)[: len(polesAp)]
if self.POD != PODGlobal:
resEx = self.trainedModel.data.projMat.dot(resEx.T)
resAp = self.trainedModel.data.projMat.dot(resAp.T)
else:
resAp = None
dist = chordalMetricAdjusted(
polesEx, polesAp, self.matchingWeightError,
resEx, resAp, self.HFEngine, self.approx_state)
err[j] = np.mean(dist[np.arange(len(polesEx)),
pointMatching(dist)])
self.trainedModel.updateEffectiveSamples(self.HFEngine, None,
self.matchingWeight, self.POD == PODGlobal, self.approx_state)
self.musMarginal.append(_musMExcl)
self.verbosity += 20
self.trainedModel.data.marginalInterp = self._setupMarginalInterp()
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 plotEstimator(self, est:Np1D, idxMax:List[int], estMax:List[float]):
if not (np.any(np.isnan(est)) or np.any(np.isinf(est))):
musre = copy(self.trainedModel.data.musMarginal.re.data)
if self._figHandle is None: self._figHandle = plt.figure()
self._figHandle.clf()
ax = self._figHandle.add_subplot(1, 1, 1)
errCP = copy(est)
while len(musre) > 0:
if self.nparMarginal == 1:
currIdx = np.arange(len(musre))
else:
currIdx = np.where(np.isclose(np.sum(
np.abs(musre[:, 1 :] - musre[0, 1 :]), 1), 0.))[0]
currIdxSorted = np.argsort(musre[currIdx, 0])
ax.semilogy(musre[currIdxSorted, 0], errCP[currIdxSorted],
'k.-', linewidth = 1)
musre = np.delete(musre, currIdx, 0)
errCP = np.delete(errCP, currIdx)
ax.semilogy(self.musMarginal.re(0),
(self.greedyTolMarginal,) * len(self.musMarginal),'*m')
if len(idxMax) > 0 and estMax is not None:
ax.semilogy(self.trainedModel.data.musMarginal.re(idxMax, 0),
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), 2)
self.musMarginal.append(mus)
self.setupApproxPivoted(mus)
self._SMarginal = len(self.musMarginal)
self._approxParameters["SMarginal"] = self.SMarginal
def greedyNextSample(self, muidx:int, plotEst : str = "NONE") \
-> Tuple[Np1D, int, float, paramVal]:
RROMPyAssert(self._mode, message = "Cannot add greedy sample.")
idxAdded = self.samplerMarginalGrid.refine(muidx)
self._addMarginalSample(self.samplerMarginalGrid.points[idxAdded])
errorEstTest, muidx, maxErrorEst = self.errorEstimator(True)
if plotEst == "ALL":
self.plotEstimator(errorEstTest, muidx, maxErrorEst)
return (errorEstTest, muidx, maxErrorEst,
self.samplerMarginalGrid.points[muidx])
def _preliminaryTraining(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 setupApproxPivoted(self, mu:paramVal):
RROMPyAssert(self._mode, message = "Cannot setup approximant.")
raise RROMPyException("Must override.")
def setupApprox(self, plotEst : str = "NONE"):
"""Compute greedy snapshots of solution map."""
if self.checkComputedApprox(): return
RROMPyAssert(self._mode, message = "Cannot start greedy algorithm.")
vbMng(self, "INIT", "Starting computation of snapshots.", 2)
self._preliminaryTraining()
muidx, max2ErrorEst = [], np.inf
while (self.samplerMarginalGrid.npoints < self.maxIterMarginal
and max2ErrorEst > self.greedyTolMarginal):
errorEstTest, muidx, maxErrorEst, mu = self.greedyNextSample(muidx,
plotEst)
if len(maxErrorEst) > 0:
max2ErrorEst = np.max(maxErrorEst)
vbMng(self, "MAIN", ("Uniform testing error estimate "
"{:.4e}.").format(max2ErrorEst), 2)
else:
max2ErrorEst = 0.
if plotEst == "LAST":
self.plotEstimator(errorEstTest, muidx, maxErrorEst)
vbMng(self, "DEL",
("Done computing snapshots (final snapshot count: "
"{}).").format(np.sum(self.samplingEngine.nsamples)), 2)
def checkComputedApprox(self) -> bool:
return (super().checkComputedApprox()
and len(self.mus) == self.trainedModel.data.projMat.shape[1])

Event Timeline

c4science · Help