Page MenuHomec4science

approximant_lagrange_greedy_pade_orthogonal.py
No OneTemporary

File Metadata

Created
Thu, Mar 28, 14:28

approximant_lagrange_greedy_pade_orthogonal.py

# 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 copy
import numpy as np
from rrompy.reduction_methods.base import checkRobustTolerance
from .generic_approximant_lagrange_greedy import (
GenericApproximantLagrangeGreedy)
from rrompy.reduction_methods.lagrange import ApproximantLagrangePadeOrthogonal
from rrompy.utilities.base.types import DictAny, List, HFEng
from rrompy.utilities.base import purgeDict, verbosityDepth
from rrompy.utilities.warning_manager import warn
__all__ = ['ApproximantLagrangePadeOrthogonalGreedy']
class ApproximantLagrangePadeOrthogonalGreedy(GenericApproximantLagrangeGreedy,
ApproximantLagrangePadeOrthogonal):
"""
ROM greedy Lagrange Pade' 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;
- 'muBounds': list of bounds for parameter values; defaults to
[[0], [1]];
- 'basis': type of basis for interpolation; allowed values include
'MONOMIAL', 'CHEBYSHEV' and 'LEGENDRE'; defaults to 'MONOMIAL';
- 'Delta': difference between M and N in rational approximant;
defaults to 0;
- 'greedyTol': uniform error tolerance for greedy algorithm;
defaults to 1e-2;
- 'maxIter': maximum number of greedy steps; defaults to 1e2;
- 'refinementRatio': ratio of training points to be exhausted
before training set refinement; defaults to 0.2;
- 'nTrainingPoints': number of training points; defaults to
maxIter / refinementRatio;
- 'nTestPoints': number of starting test points; defaults to 1;
- 'trainingSetGenerator': training sample points generator;
defaults to uniform sampler within muBounds;
- 'interpRcond': tolerance for interpolation via numpy.polyfit;
defaults to None;
- 'robustTol': tolerance for robust Pade' denominator management;
defaults to 0.
Defaults to empty dict.
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.
parameterList: Recognized keys of approximant parameters:
- 'POD': whether to compute POD of snapshots;
- 'muBounds': list of bounds for parameter values;
- 'basis': type of basis for interpolation;
- 'Delta': difference between M and N in rational approximant;
- 'greedyTol': uniform error tolerance for greedy algorithm;
- 'maxIter': maximum number of greedy steps;
- 'refinementRatio': ratio of training points to be exhausted
before training set refinement;
- 'nTrainingPoints': number of training points;
- 'nTestPoints': number of starting test points;
- 'trainingSetGenerator': training sample points generator;
- 'interpRcond': tolerance for interpolation via numpy.polyfit;
- 'robustTol': tolerance for robust Pade' denominator management.
extraApproxParameters: List of approxParameters keys in addition to
mother class's.
POD: whether to compute POD of snapshots.
muBounds: list of bounds for parameter values.
greedyTol: uniform error tolerance for greedy algorithm.
maxIter: maximum number of greedy steps.
refinementRatio: ratio of training points to be exhausted before
training set refinement.
nTrainingPoints: number of training points.
nTestPoints: number of starting test points.
trainingSetGenerator: training sample points generator.
interpRcond: Tolerance for interpolation via numpy.polyfit.
robustTol: Tolerance for robust Pade' denominator management.
samplingEngine: Sampling engine.
uHF: High fidelity solution with wavenumber lastSolvedHF as numpy
complex vector.
lastSolvedHF: Wavenumber corresponding to last computed high fidelity
solution.
Q: Numpy 1D vector containing complex coefficients of approximant
denominator.
P: Numpy 2D vector whose columns are FE dofs of coefficients of
approximant numerator.
uApp: Last evaluated approximant as numpy complex vector.
lastApproxParameters: List of parameters corresponding to last
computed approximant.
"""
def __init__(self, HFEngine:HFEng, mu0 : complex = 0.,
approxParameters : DictAny = {}, homogeneized : bool = False,
verbosity : int = 10):
self._preInit()
self._addParametersToList(["basis", "Delta", "interpRcond",
"robustTol"])
super().__init__(HFEngine = HFEngine, mu0 = mu0,
approxParameters = approxParameters,
homogeneized = homogeneized,
verbosity = verbosity)
self._postInit()
@property
def approxParameters(self):
"""
Value of approximant parameters. Its assignment may change robustTol.
"""
return self._approxParameters
@approxParameters.setter
def approxParameters(self, approxParams):
approxParameters = purgeDict(approxParams, self.parameterList,
dictname = self.name() + ".approxParameters",
baselevel = 1)
approxParametersCopy = purgeDict(approxParameters, ["basis",
"Delta",
"interpRcond",
"robustTol"],
True, True, baselevel = 1)
if "Delta" in list(approxParameters.keys()):
self._Delta = approxParameters["Delta"]
elif hasattr(self, "Delta"):
self._Delta = self.Delta
else:
self._Delta = 0
GenericApproximantLagrangeGreedy.approxParameters.fset(self,
approxParametersCopy)
keyList = list(approxParameters.keys())
self.Delta = self.Delta
if "basis" in keyList or not hasattr(self, "_val"):
if ("basis" in keyList
and approxParameters["basis"].upper() == "CHEBYSHEV"):
self.basis = "CHEBYSHEV"
self._val = np.polynomial.chebyshev.chebval
self._vander = np.polynomial.chebyshev.chebvander
self._fit = np.polynomial.chebyshev.chebfit
self._roots = np.polynomial.chebyshev.chebroots
self._domcoeff = lambda n: 2. ** (n - 1) if n > 0 else 1.
elif ("basis" in keyList
and approxParameters["basis"].upper() == "LEGENDRE"):
self.basis = "LEGENDRE"
self._val = np.polynomial.legendre.legval
self._vander = np.polynomial.legendre.legvander
self._fit = np.polynomial.legendre.legfit
self._roots = np.polynomial.legendre.legroots
from scipy.special import binom
self._domcoeff = lambda n: (2. ** n * (np.pi * n) ** -.5
if n > 10 else
.5 ** n * binom(2 * n, n))
else:
self.basis = "UNIFORM"
self._val = np.polynomial.polynomial.polyval
self._vander = np.polynomial.polynomial.polyvander
self._fit = np.polynomial.polynomial.polyfit
self._roots = np.polynomial.polynomial.polyroots
self._domcoeff = lambda n: 1.
if "interpRcond" in keyList:
self.interpRcond = approxParameters["interpRcond"]
elif hasattr(self, "interpRcond"):
self.interpRcond = self.interpRcond
else:
self.interpRcond = None
if "robustTol" in keyList:
self.robustTol = approxParameters["robustTol"]
elif hasattr(self, "robustTol"):
self.robustTol = self.robustTol
else:
self.robustTol = 0
@property
def Delta(self):
"""Value of Delta."""
return self._Delta
@Delta.setter
def Delta(self, Delta):
if not np.isclose(Delta, np.floor(Delta)):
raise ArithmeticError("Delta must be an integer.")
if Delta < 0:
warn(("Error estimator unreliable for Delta < 0. Overloading of "
"errorEstimator is suggested."))
else:
Deltamin = (max(self.HFEngine.nbs,
self.HFEngine.nAs * self.homogeneized)
- 1 - 1 * (self.HFEngine.nAs > 1))
if Delta < Deltamin:
warn(("Method unreliable for selected Delta. Suggested "
"minimal value of Delta: {}.").format(Deltamin))
self._Delta = Delta
self._approxParameters["Delta"] = self.Delta
@property
def nTestPoints(self):
"""Value of nTestPoints."""
return self._nTestPoints
@nTestPoints.setter
def nTestPoints(self, nTestPoints):
if nTestPoints <= np.abs(self.Delta):
warn(("nTestPoints must be at least abs(Delta) + 1. Increasing "
"value to abs(Delta) + 1."))
nTestPoints = np.abs(self.Delta) + 1
if not np.isclose(nTestPoints, np.int(nTestPoints)):
raise ArithmeticError("nTestPoints must be an integer.")
nTestPoints = np.int(nTestPoints)
if hasattr(self, "nTestPoints"):
nTestPointsold = self.nTestPoints
else: nTestPointsold = -1
self._nTestPoints = nTestPoints
self._approxParameters["nTestPoints"] = self.nTestPoints
if nTestPointsold != self.nTestPoints:
self.resetSamples()
def resetSamples(self):
"""Reset samples."""
super().resetSamples()
self.qs = np.empty(0, dtype = np.complex)
def errorEstimator(self, mus:List[np.complex]) -> List[np.complex]:
"""
Standard residual-based error estimator. Unreliable for unstable
problems.
"""
self.setupApprox()
self.initEstNormer()
nmus = len(mus)
if self.N + 1 == self.S:
QSf = self._domcoeff(self.N) * self.Q[-1] * self.HFEngine.b(
self.mu0, 1, homogeneized = self.homogeneized)
else:
QSf = 0
L1P = self._domcoeff(self.M) * self.HFEngine.A(self.mu0, 1).dot(
self.samplingEngine.samples.dot(self.P[:, -1]))
jOpt = self.estNormer.norm(QSf - L1P)
if np.isnan(jOpt) or np.isinf(jOpt):
err = np.empty(nmus)
err[:] = np.inf
return err
musTile = np.tile(self.HFEngine.rescaling(mus).reshape(-1, 1),
[1, len(self.mus)])
smusCol = self.HFEngine.rescaling(self.mus).reshape(1, -1)
mussmus = np.abs(musTile - smusCol)
num = np.prod(mussmus, axis = 1)
den = np.abs(self.getQVal(mus))
RHSnorms = np.empty(nmus)
if self.HFEngine.nbs == 1:
RHS = self.getRHS(mus[0], homogeneized = self.homogeneized)
RHSnorms[:] = self.estNormer.norm(RHS)
else:
for j in range(nmus):
RHS = self.getRHS(mus[j], homogeneized = self.homogeneized)
RHSnorms[j] = self.estNormer.norm(RHS)
return jOpt * num / den / RHSnorms / self.scaleFactor ** (self.S - 1)
def setupApprox(self):
"""
Compute Pade' interpolant.
SVD-based robust eigenvalue management.
"""
if not self.checkComputedApprox():
if self.verbosity >= 5:
verbosityDepth("INIT", "Setting up {}.". format(self.name()))
self.computeRescaleParameter()
self.S = len(self.mus)
self._M = self.S - 1
self._N = self.S - 1
if self.Delta < 0:
self._M += self.Delta
else:
self._N -= self.Delta
if min(self.M, self.N) < 0:
if self.verbosity >= 5:
verbosityDepth("MAIN", "Minimal sample size not achieved.")
self.Q = np.ones(1, dtype = np.complex)
self.P = np.diag([np.nan] * len(self.mus))
self.lastApproxParameters = copy(self.approxParameters)
if self.verbosity >= 5:
verbosityDepth("DEL", ("Aborting computation of "
"approximant.\n"))
return
self.greedy()
if self.N > 0:
if self.verbosity >= 7:
verbosityDepth("INIT", ("Starting computation of "
"denominator."))
TN = self._vander(self.radiusPade(self.mus), self.N)
while self.N > 0:
TN = TN[:, : self.N + 1]
if self.POD:
data = self.samplingEngine.RPOD
else:
data = self.samplingEngine.samples
RHSFull = np.empty((self.S, data.shape[0] * (self.N + 1)),
dtype = np.complex)
for j in range(self.S):
RHSFull[j, :] = np.kron(data[:, j], TN[j, :])
fitOut = self._fit(self.radiusPade(self.mus), RHSFull,
self.S - 1, full = True,
rcond = self.interpRcond)
if self.verbosity >= 2:
verbosityDepth("MAIN", ("Fitting {} samples with "
"degree {} through {}... "
"Conditioning of system: "
"{:.4e}.").format(self.S,
self.S - 1, self._fit.__name__,
fitOut[1][2][0] / fitOut[1][2][-1]))
if fitOut[1][1] < self.S:
warn(("Polyfit is poorly conditioned. Starting "
"preemptive termination of computation of "
"approximant."))
self.Q = np.ones(1, dtype = np.complex)
self.P = np.diag([np.nan] * len(self.mus))
self.lastApproxParameters = copy(self.approxParameters)
if hasattr(self, "lastSolvedApp"):
del self.lastSolvedApp
if self.verbosity >= 7:
verbosityDepth("DEL", ("Aborting computation of "
"denominator."))
if self.verbosity >= 5:
verbosityDepth("DEL", ("Aborting computation of "
"approximant.\n"))
return
G = fitOut[0][-1, :].reshape((data.shape[0], self.N + 1))
if self.POD:
if self.verbosity >= 7:
verbosityDepth("INIT", ("Solving svd for square "
"root of gramian matrix."),
end = "")
_, ev, eV = np.linalg.svd(G, full_matrices = False)
ev = ev[::-1]
eV = eV[::-1, :].conj().T
else:
if self.verbosity >= 10:
verbosityDepth("INIT", "Building gramian matrix.",
end = "")
G2 = self.HFEngine.innerProduct(G, G)
if self.verbosity >= 10:
verbosityDepth("DEL", "Done building gramian.",
inline = True)
if self.verbosity >= 7:
verbosityDepth("INIT", ("Solving eigenvalue "
"problem for gramian "
"matrix."), end = "")
ev, eV = np.linalg.eigh(G2)
if self.verbosity >= 7:
verbosityDepth("DEL", " Done.", inline = True)
newParameters = checkRobustTolerance(ev, self.M,
self.robustTol)
if not newParameters:
break
self._N = newParameters["N"]
self._M = newParameters["E"]
if self.N <= 0:
eV = np.ones((1, 1))
self.Q = eV[:, 0]
if self.verbosity >= 7:
verbosityDepth("DEL", "Done computing denominator.")
else:
self.Q = np.ones(1, dtype = np.complex)
if self.verbosity >= 7:
verbosityDepth("INIT", "Starting computation of numerator.")
self.lastApproxParameters = copy(self.approxParameters)
Qevaldiag = np.diag(self.getQVal(self.mus))
while self.M >= 0:
fitOut = self._fit(self.radiusPade(self.mus), Qevaldiag,
self.M, full = True,
rcond = self.interpRcond)
if fitOut[1][1] == self.M + 1:
P = fitOut[0].T
break
warn(("Polyfit is poorly conditioned. Reducing M from {} to "
"{}. Exact snapshot interpolation not guaranteed.")\
.format(self.M, fitOut[1][1] - 1))
self._M = fitOut[1][1] - 1
self.P = np.atleast_2d(P)
if self.POD:
self.P = self.samplingEngine.RPOD.dot(self.P)
if self.verbosity >= 7:
verbosityDepth("DEL", "Done computing numerator.")
self.lastApproxParameters = copy(self.approxParameters)
if hasattr(self, "lastSolvedApp"): del self.lastSolvedApp
if self.verbosity >= 5:
verbosityDepth("DEL", "Done setting up approximant.\n")

Event Timeline