diff --git a/rrompy/reduction_methods/distributed/rational_interpolant.py b/rrompy/reduction_methods/distributed/rational_interpolant.py
index bc8be42..3a517aa 100644
--- a/rrompy/reduction_methods/distributed/rational_interpolant.py
+++ b/rrompy/reduction_methods/distributed/rational_interpolant.py
@@ -1,534 +1,532 @@
# 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 scipy.special import factorial as fact
from rrompy.reduction_methods.base import checkRobustTolerance
from .generic_distributed_approximant import GenericDistributedApproximant
-from rrompy.utilities.poly_fitting import (polybases, polyvander, polyfitname,
- customFit)
+from rrompy.utilities.poly_fitting import customFit
+from rrompy.utilities.poly_fitting.polynomial import (polybases, polyfitname,
+ polyvander)
from rrompy.reduction_methods.trained_model import TrainedModelPade as tModel
from rrompy.reduction_methods.trained_model import TrainedModelData
from rrompy.utilities.base.types import Np1D, Np2D, HFEng, DictAny, Tuple
from rrompy.utilities.base import verbosityDepth, purgeDict
from rrompy.utilities.exception_manager import (RROMPyException, RROMPyAssert,
RROMPyWarning)
__all__ = ['RationalInterpolant']
class RationalInterpolant(GenericDistributedApproximant):
"""
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;
- 'muBounds': list of bounds for parameter values; defaults to
[0, 1];
- 'S': total number of samples current approximant relies upon;
defaults to 2;
- 'sampler': sample point generator; defaults to uniform sampler on
muBounds;
- 'polybasis': type of polynomial basis for interpolation; allowed
values include 'MONOMIAL', 'CHEBYSHEV' and 'LEGENDRE'; defaults
to 'MONOMIAL';
- 'E': coefficient of interpolant to be minimized; defaults to
min(S, M + 1);
- 'M': degree of Pade' interpolant numerator; defaults to 0;
- 'N': degree of Pade' interpolant denominator; defaults to 0;
- 'interpRcond': tolerance for interpolation; defaults to None;
- 'robustTol': tolerance for robust Pade' denominator management;
defaults to 0.
Defaults to empty dict.
homogeneized(optional): Whether to homogeneize Dirichlet BCs. Defaults
to False.
verbosity(optional): Verbosity level. Defaults to 10.
Attributes:
HFEngine: HF problem solver.
mu0: Default parameter.
mus: Array of snapshot parameters.
ws: Array of snapshot weigths.
homogeneized: Whether to homogeneize Dirichlet BCs.
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;
- 'S': total number of samples current approximant relies upon;
- 'sampler': sample point generator;
- 'polybasis': type of polynomial basis for interpolation;
- 'E': coefficient of interpolant to be minimized;
- 'M': degree of Pade' interpolant numerator;
- 'N': degree of Pade' interpolant denominator;
- '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.
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.
interpRcond: Tolerance for interpolation via numpy.polyfit.
robustTol: Tolerance for robust Pade' denominator management.
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.
uAppReduced: Reduced approximate solution(s) with parameter(s)
lastSolvedApp as sampleList.
lastSolvedAppReduced: Parameter(s) corresponding to last computed
reduced approximate solution(s) as parameterList.
uApp: Approximate solution(s) with parameter(s) lastSolvedApp as
sampleList.
lastSolvedApp: Parameter(s) corresponding to last computed approximate
solution(s) as parameterList.
Q: Numpy 1D vector containing complex coefficients of approximant
denominator.
P: Numpy 2D vector whose columns are FE dofs of coefficients of
approximant numerator.
"""
def __init__(self, HFEngine:HFEng, mu0 : complex = 0.,
approxParameters : DictAny = {}, homogeneized : bool = False,
verbosity : int = 10, timestamp : bool = True):
self._preInit()
self._addParametersToList(["polybasis", "E", "M", "N",
"interpRcond", "robustTol"])
super().__init__(HFEngine = HFEngine, mu0 = mu0,
approxParameters = approxParameters,
homogeneized = homogeneized,
verbosity = verbosity, timestamp = timestamp)
self._postInit()
@property
def approxParameters(self):
"""
Value of approximant parameters.
"""
return self._approxParameters
@approxParameters.setter
def approxParameters(self, approxParams):
approxParameters = purgeDict(approxParams, self.parameterList,
dictname = self.name() + ".approxParameters",
baselevel = 1)
approxParametersCopy = purgeDict(approxParameters, ["polybasis",
"E", "M", "N",
"interpRcond",
"robustTol"],
True, True, baselevel = 1)
if hasattr(self, "_M") and self.M is not None:
Mold = self.M
self._M = 0
if hasattr(self, "_N") and self.N is not None:
Nold = self.N
self._N = 0
if hasattr(self, "_E") and self.E is not None:
self._E = 0
GenericDistributedApproximant.approxParameters.fset(self,
approxParametersCopy)
keyList = list(approxParameters.keys())
if "polybasis" in keyList:
self.polybasis = approxParameters["polybasis"]
elif not hasattr(self, "_polybasis") or self._polybasis is None:
self.polybasis = "MONOMIAL"
if "interpRcond" in keyList:
self.interpRcond = approxParameters["interpRcond"]
elif not hasattr(self, "interpRcond") or self.interpRcond is None:
self.interpRcond = None
if "robustTol" in keyList:
self.robustTol = approxParameters["robustTol"]
elif not hasattr(self, "_robustTol") or self._robustTol is None:
self.robustTol = 0
if "M" in keyList:
self.M = approxParameters["M"]
elif hasattr(self, "_M") and self.M is not None:
self.M = Mold
else:
self.M = 0
if "N" in keyList:
self.N = approxParameters["N"]
elif hasattr(self, "_N") and self.N is not None:
self.N = Nold
else:
self.N = 0
if "E" in keyList:
self.E = approxParameters["E"]
else:
self.E = min(self.S - 1, self.M + 1)
@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("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 interpRcond(self):
"""Value of interpRcond."""
return self._interpRcond
@interpRcond.setter
def interpRcond(self, interpRcond):
self._interpRcond = interpRcond
self._approxParameters["interpRcond"] = self.interpRcond
@property
def M(self):
"""Value of M. Its assignment may change S."""
return self._M
@M.setter
def M(self, M):
if M < 0: raise RROMPyException("M must be non-negative.")
self._M = M
self._approxParameters["M"] = self.M
if hasattr(self, "_S") and self.S < self.M + 1:
RROMPyWarning("Prescribed S is too small. Updating S to M + 1.")
self.S = self.M + 1
@property
def N(self):
"""Value of N. Its assignment may change S."""
return self._N
@N.setter
def N(self, N):
if N < 0: raise RROMPyException("N must be non-negative.")
self._N = N
self._approxParameters["N"] = self.N
if hasattr(self, "_S") and self.S < self.N + 1:
RROMPyWarning("Prescribed S is too small. Updating S to N + 1.")
self.S = self.N + 1
@property
def E(self):
"""Value of E. Its assignment may change S."""
return self._E
@E.setter
def E(self, E):
if E < 0: raise RROMPyException("E must be non-negative.")
self._E = E
self._approxParameters["E"] = self.E
if hasattr(self, "_S") and self.S < self.E + 1:
RROMPyWarning("Prescribed S is too small. Updating S to E + 1.")
self.S = self.E + 1
@property
def robustTol(self):
"""Value of tolerance for robust Pade' 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 S(self):
"""Value of S."""
return self._S
@S.setter
def S(self, S):
if S <= 0: raise RROMPyException("S must be positive.")
if hasattr(self, "_S"): Sold = self.S
else: Sold = -1
vals, label = [0] * 3, {0:"M", 1:"N", 2:"E"}
if hasattr(self, "_M") and self._M is not None: vals[0] = self.M
if hasattr(self, "_N") and self._N is not None: vals[1] = self.N
if hasattr(self, "_E") and self._E is not None: vals[2] = self.E
idxmax = np.argmax(vals)
if vals[idxmax] + 1 > S:
RROMPyWarning(("Prescribed S is too small. Updating S to {} + "
"1.").format(label[idxmax]))
self.S = vals[idxmax] + 1
else:
self._S = S
self._approxParameters["S"] = self.S
if Sold != self.S:
self.resetSamples()
def _setupDenominator(self):
"""Compute Pade' denominator."""
if self.verbosity >= 7:
verbosityDepth("INIT", "Starting computation of denominator.",
timestamp = self.timestamp)
while self.N > 0:
- TE = polyvander[self.polybasis](self.centerNormalize(self.mus),
- self.E,
- scl = 1. / self.scaleFactor)
+ TE = polyvander(self.centerNormalize(self.mus), self.E,
+ self.polybasis)#, scl = 1. / self.scaleFactor)
TE = (TE.T * self.ws).T
RHS = np.zeros(self.E + 1)
RHS[-1] = 1.
fitOut = customFit(TE.T, RHS, full = True,
rcond = self.interpRcond)
if self.verbosity >= 5:
condfit = fitOut[1][2][0] / fitOut[1][2][-1]
verbosityDepth("MAIN", ("Fitting {} samples with degree {} "
"through {}... Conditioning of LS "
"system: {:.4e}.").format(
self.S, self.E,
- polyfitname[self.polybasis],
+ polyfitname(self.polybasis),
condfit),
timestamp = self.timestamp)
if fitOut[1][1] < self.E + 1:
Enew = fitOut[1][1] - 1
Nnew = min(self.N, Enew)
Mnew = min(self.M, Enew)
if Nnew == self.N:
strN = ""
else:
strN = "N from {} to {} and ".format(self.N, Nnew)
if Mnew == self.M:
strM = ""
else:
strM = "M from {} to {} and ".format(self.M, Mnew)
RROMPyWarning(("Polyfit is poorly conditioned.\nReducing {}{}"
"E from {} to {}.").format(strN, strM,
self.E, Enew))
newParams = {"N" : Nnew, "M" : Mnew, "E" : Enew}
self.approxParameters = newParams
continue
mus_un, idx_un, cnt_un = np.unique(self.mus, return_inverse = True,
return_counts = True)
- TE = polyvander[self.polybasis](self.centerNormalize(self.mus),
- self.N,
- scl = 1. / self.scaleFactor)
+ TE = polyvander(self.centerNormalize(self.mus), self.N,
+ self.polybasis)#, scl = 1. / self.scaleFactor)
TE = (TE.T * self.ws).T
if len(mus_un) == len(self.mus):
Ghalf = (TE.T * fitOut[0]).T
else:
pseudoInv = np.zeros((len(self.mus), len(self.mus)),
dtype = np.complex)
for j in range(len(mus_un)):
pseudoInv_loc = np.zeros((cnt_un[j], cnt_un[j]),
dtype = np.complex)
mask = np.arange(len(self.mus))[idx_un == j]
for der in range(cnt_un[j]):
fitderj = fitOut[0][mask[der]]
pseudoInv_loc = (pseudoInv_loc + fitderj
* np.diag(np.ones(1 + der),
k = der - cnt_un[j] + 1))
I = np.ix_(mask, mask)
pseudoInv[I] = np.flipud(pseudoInv_loc)
Ghalf = pseudoInv.dot(TE)
if self.POD:
self.Ghalf = self.samplingEngine.RPOD.dot(Ghalf)
ev, eV = self.findeveVGQR()
else:
self.Ghalf = self.samplingEngine.samples.dot(Ghalf)
ev, eV = self.findeveVGExplicit()
newParams = checkRobustTolerance(ev, self.E, self.robustTol)
if not newParams:
break
self.approxParameters = newParams
if self.N <= 0:
self._N = 0
eV = np.ones((1, 1))
if self.verbosity >= 7:
verbosityDepth("DEL", "Done computing denominator.",
timestamp = self.timestamp)
return eV[:, 0]
def _setupNumerator(self):
"""Compute Pade' numerator."""
if self.verbosity >= 7:
verbosityDepth("INIT", "Starting computation of numerator.",
timestamp = self.timestamp)
Qevaldiag = np.diag(self.trainedModel.getQVal(self.mus))
verb = self.trainedModel.verbosity
self.trainedModel.verbosity = 0
mus_un, idx_un, cnt_un = np.unique(self.mus, return_inverse = True,
return_counts = True)
for j in range(len(mus_un)):
if cnt_un[j] > 1:
Q_loc = np.zeros((cnt_un[j], cnt_un[j]), dtype = np.complex)
for der in range(1, cnt_un[j]):
Qderj = (self.trainedModel.getQVal(mus_un[j], der,
scl = 1. / self.scaleFactor)
/ fact(der))
Q_loc = Q_loc + Qderj * np.diag(np.ones(cnt_un[j] - der),
k = - der)
I = np.ix_(idx_un == j, idx_un == j)
Qevaldiag[I] = Qevaldiag[I] + Q_loc
self.trainedModel.verbosity = verb
while self.M >= 0:
- fitVander = polyvander[self.polybasis](
- self.centerNormalize(self.mus),
- self.M, scl = 1. / self.scaleFactor)
+ fitVander = polyvander(self.centerNormalize(self.mus), self.M,
+ self.polybasis)#, scl = 1. / self.scaleFactor)
fitOut = customFit(fitVander, Qevaldiag, w = self.ws, full = True,
rcond = self.interpRcond)
if self.verbosity >= 5:
condfit = fitOut[1][2][0] / fitOut[1][2][-1]
verbosityDepth("MAIN", ("Fitting {} samples with degree {} "
"through {}... Conditioning of LS "
"system: {:.4e}.").format(
self.S, self.M,
- polyfitname[self.polybasis],
+ polyfitname(self.polybasis),
condfit),
timestamp = self.timestamp)
if fitOut[1][1] == self.M + 1:
P = fitOut[0].T
break
RROMPyWarning(("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
if self.M <= 0:
raise RROMPyException(("Instability in computation of numerator. "
"Aborting."))
if self.verbosity >= 7:
verbosityDepth("DEL", "Done computing numerator.",
timestamp = self.timestamp)
return np.atleast_2d(P)
def setupApprox(self):
"""
Compute Pade' interpolant.
SVD-based robust eigenvalue management.
"""
if self.checkComputedApprox():
return
RROMPyAssert(self._mode, message = "Cannot setup approximant.")
if self.verbosity >= 5:
verbosityDepth("INIT", "Setting up {}.". format(self.name()),
timestamp = self.timestamp)
self.computeScaleFactor()
self.computeSnapshots()
if self.trainedModel is None:
self.trainedModel = tModel()
self.trainedModel.verbosity = self.verbosity
self.trainedModel.timestamp = self.timestamp
data = TrainedModelData(self.trainedModel.name(), self.mu0,
self.samplingEngine.samples,
self.HFEngine.rescalingExp)
data.polytype = self.polybasis
data.scaleFactor = self.scaleFactor
data.mus = copy(self.mus)
self.trainedModel.data = data
else:
self.trainedModel = self.trainedModel
self.trainedModel.data.projMat = copy(self.samplingEngine.samples)
if self.N > 0:
Q = self._setupDenominator()
else:
Q = np.ones(1, dtype = np.complex)
self.trainedModel.data.Q = copy(Q)
P = self._setupNumerator()
if self.POD:
P = self.samplingEngine.RPOD.dot(P)
self.trainedModel.data.P = copy(P)
self.trainedModel.data.approxParameters = copy(self.approxParameters)
if self.verbosity >= 5:
verbosityDepth("DEL", "Done setting up approximant.",
timestamp = self.timestamp)
def findeveVGExplicit(self, verbOutput : int = 5) -> Tuple[Np1D, Np2D]:
"""
Compute explicitly eigenvalues and eigenvectors of Pade' denominator
matrix.
"""
RROMPyAssert(self._mode, message = "Cannot solve eigenvalue problem.")
if self.verbosity >= 10:
verbosityDepth("INIT", "Building gramian matrix.",
timestamp = self.timestamp)
self.G = self.HFEngine.innerProduct(self.Ghalf, self.Ghalf)
if self.verbosity >= 10:
verbosityDepth("DEL", "Done building gramian.",
timestamp = self.timestamp)
if self.verbosity >= 7:
verbosityDepth("INIT", ("Solving eigenvalue problem for gramian "
"matrix."), timestamp = self.timestamp)
ev, eV = np.linalg.eigh(self.G)
if self.verbosity >= verbOutput:
try: condev = ev[-1] / ev[0]
except: condev = np.inf
verbosityDepth("MAIN", ("Solved eigenvalue problem of size {} "
"with condition number {:.4e}.").format(
self.N + 1, condev),
timestamp = self.timestamp)
if self.verbosity >= 7:
verbosityDepth("DEL", "Done solving eigenvalue problem.",
timestamp = self.timestamp)
return ev, eV
def findeveVGQR(self, verbOutput : int = 5) -> Tuple[Np1D, Np2D]:
"""
Compute eigenvalues and eigenvectors of Pade' denominator matrix
through SVD of R factor.
"""
if self.verbosity >= 7:
verbosityDepth("INIT", ("Solving svd for square root of gramian "
"matrix."), timestamp = self.timestamp)
_, s, eV = np.linalg.svd(self.Ghalf, full_matrices = False)
ev = s[::-1]
eV = eV[::-1, :].T.conj()
if self.verbosity >= verbOutput:
try: condev = s[0] / s[-1]
except: condev = np.inf
verbosityDepth("MAIN", ("Solved svd problem of size {} x {} with "
"condition number {:.4e}.").format(
self.S, self.N + 1, condev),
timestamp = self.timestamp)
if self.verbosity >= 7:
verbosityDepth("DEL", "Done solving eigenvalue problem.",
timestamp = self.timestamp)
return ev, eV
def centerNormalize(self, mu:Np1D, mu0 : float = None) -> float:
"""
Compute normalized parameter to be plugged into approximant.
Args:
mu: Parameter(s) 1.
mu0: Parameter(s) 2. If None, set to self.mu0.
Returns:
Normalized parameter.
"""
return self.trainedModel.centerNormalize(mu, mu0)
def getResidues(self) -> Np1D:
"""
Obtain approximant residues.
Returns:
Matrix with residues as columns.
"""
return self.trainedModel.getResidues()
diff --git a/rrompy/reduction_methods/distributed_greedy/rational_interpolant_greedy.py b/rrompy/reduction_methods/distributed_greedy/rational_interpolant_greedy.py
index de0426e..0522492 100644
--- a/rrompy/reduction_methods/distributed_greedy/rational_interpolant_greedy.py
+++ b/rrompy/reduction_methods/distributed_greedy/rational_interpolant_greedy.py
@@ -1,580 +1,581 @@
# 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 scipy.special import factorial as fact
from .generic_distributed_greedy_approximant import \
GenericDistributedGreedyApproximant
-from rrompy.utilities.poly_fitting import (polybases, polyvander, polydomcoeff,
- polyfitname, customFit)
+from rrompy.utilities.poly_fitting import customFit
+from rrompy.utilities.poly_fitting.polynomial import (polybases, polyfitname,
+ polydomcoeff, polyvander)
from rrompy.reduction_methods.distributed import RationalInterpolant
from rrompy.reduction_methods.trained_model import TrainedModelPade as tModel
from rrompy.reduction_methods.trained_model import TrainedModelData
from rrompy.utilities.base.types import Np1D, Np2D, DictAny, HFEng
from rrompy.utilities.base import purgeDict, verbosityDepth
from rrompy.utilities.exception_manager import RROMPyWarning
from rrompy.utilities.exception_manager import RROMPyException
__all__ = ['RationalInterpolantGreedy']
class RationalInterpolantGreedy(GenericDistributedGreedyApproximant,
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;
- 'muBounds': list of bounds for parameter values; defaults to
[0, 1];
- 'S': number of starting training points; defaults to 2;
- 'sampler': sample point generator; defaults to uniform sampler on
muBounds;
- '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;
- 'errorEstimatorKind': kind of error estimator; available values
include 'EXACT', 'BASIC', and 'BARE'; defaults to 'EXACT';
- '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;
- 'nTestPoints': number of test points; defaults to maxIter /
refinementRatio;
- 'trainSetGenerator': training sample points generator; defaults
to Chebyshev 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.
homogeneized(optional): Whether to homogeneize Dirichlet BCs. Defaults
to False.
verbosity(optional): Verbosity level. Defaults to 10.
Attributes:
HFEngine: HF problem solver.
mu0: Default parameter.
mus: Array of snapshot parameters.
homogeneized: Whether to homogeneize Dirichlet BCs.
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;
- 'S': number of starting training points;
- 'sampler': sample point generator;
- 'basis': type of basis for interpolation;
- 'Delta': difference between M and N in rational approximant;
- 'greedyTol': uniform error tolerance for greedy algorithm;
- 'errorEstimatorKind': kind of error estimator;
- 'maxIter': maximum number of greedy steps;
- 'refinementRatio': ratio of training points to be exhausted
before training set refinement;
- 'nTestPoints': number of test points;
- 'trainSetGenerator': 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.
S: number of starting training points.
sampler: Sample point generator.
greedyTol: uniform error tolerance for greedy algorithm.
errorEstimatorKind: kind of error estimator.
maxIter: maximum number of greedy steps.
refinementRatio: ratio of training points to be exhausted before
training set refinement.
nTestPoints: number of starting training points.
trainSetGenerator: training sample points generator.
interpRcond: Tolerance for interpolation via numpy.polyfit.
robustTol: Tolerance for robust Pade' denominator management.
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.
uAppReduced: Reduced approximate solution(s) with parameter(s)
lastSolvedApp as sampleList.
lastSolvedAppReduced: Parameter(s) corresponding to last computed
reduced approximate solution(s) as parameterList.
uApp: Approximate solution(s) with parameter(s) lastSolvedApp as
sampleList.
lastSolvedApp: Parameter(s) corresponding to last computed approximate
solution(s) as parameterList.
"""
_allowedEstimatorKinds = ["EXACT", "BASIC", "BARE"]
def __init__(self, HFEngine:HFEng, mu0 : complex = 0.,
approxParameters : DictAny = {}, homogeneized : bool = False,
verbosity : int = 10, timestamp : bool = True):
self._preInit()
self._addParametersToList(["polybasis", "Delta", "errorEstimatorKind",
"interpRcond", "robustTol"])
super().__init__(HFEngine = HFEngine, mu0 = mu0,
approxParameters = approxParameters,
homogeneized = homogeneized,
verbosity = verbosity, timestamp = timestamp)
if self.verbosity >= 7:
verbosityDepth("INIT", "Computing Taylor blocks of system.",
timestamp = self.timestamp)
nAs = self.HFEngine.nAs - 1
nbs = max(self.HFEngine.nbs, (nAs + 1) * self.homogeneized)
self.As = [self.HFEngine.A(self.mu0, j + 1) for j in range(nAs)]
self.bs = [self.HFEngine.b(self.mu0, j, self.homogeneized)
for j in range(nbs)]
if self.verbosity >= 7:
verbosityDepth("DEL", "Done computing Taylor blocks.",
timestamp = self.timestamp)
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, ["polybasis",
"Delta",
"errorEstimatorKind",
"interpRcond",
"robustTol"],
True, True, baselevel = 1)
if "Delta" in list(approxParameters.keys()):
self._Delta = approxParameters["Delta"]
elif not hasattr(self, "_Delta") or self._Delta is None:
self._Delta = 0
GenericDistributedGreedyApproximant.approxParameters.fset(self,
approxParametersCopy)
keyList = list(approxParameters.keys())
self.Delta = self.Delta
if "polybasis" in keyList:
self.polybasis = approxParameters["polybasis"]
elif not hasattr(self, "_polybasis") or self._polybasis is None:
self.polybasis = "MONOMIAL"
if "errorEstimatorKind" in keyList:
self.errorEstimatorKind = approxParameters["errorEstimatorKind"]
elif (not hasattr(self, "_errorEstimatorKind")
or self.errorEstimatorKind is None):
self.errorEstimatorKind = "EXACT"
if "interpRcond" in keyList:
self.interpRcond = approxParameters["interpRcond"]
elif not hasattr(self, "interpRcond") or self.interpRcond is None:
self.interpRcond = None
if "robustTol" in keyList:
self.robustTol = approxParameters["robustTol"]
elif not hasattr(self, "_robustTol") or self._robustTol is None:
self.robustTol = 0
@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 Delta(self):
"""Value of Delta."""
return self._Delta
@Delta.setter
def Delta(self, Delta):
if not np.isclose(Delta, np.floor(Delta)):
raise RROMPyException("Delta must be an integer.")
if Delta < 0:
RROMPyWarning(("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:
RROMPyWarning(("Method may be unreliable for selected Delta. "
"Suggested minimal value of Delta: {}.").format(
Deltamin))
self._Delta = Delta
self._approxParameters["Delta"] = self.Delta
@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 'EXACT'."))
errorEstimatorKind = "EXACT"
self._errorEstimatorKind = errorEstimatorKind
self._approxParameters["errorEstimatorKind"] = self.errorEstimatorKind
@property
def nTestPoints(self):
"""Value of nTestPoints."""
return self._nTestPoints
@nTestPoints.setter
def nTestPoints(self, nTestPoints):
if nTestPoints <= np.abs(self.Delta):
RROMPyWarning(("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 RROMPyException("nTestPoints must be an integer.")
nTestPoints = np.int(nTestPoints)
if hasattr(self, "_nTestPoints") and self.nTestPoints is not None:
nTestPointsold = self.nTestPoints
else: nTestPointsold = -1
self._nTestPoints = nTestPoints
self._approxParameters["nTestPoints"] = self.nTestPoints
if nTestPointsold != self.nTestPoints:
self.resetSamples()
def _errorSamplingRatio(self, mus:Np1D, muRTest:Np1D,
muRTrain:Np1D) -> Np1D:
"""Scalar ratio in explicit error estimator."""
testTile = np.tile(np.reshape(muRTest, (-1, 1)), [1, len(muRTrain)])
nodalVals = np.prod(testTile - np.reshape(muRTrain, (1, -1)), axis = 1)
denVals = self.trainedModel.getQVal(mus)
return np.abs(nodalVals / denVals)
def _RHSNorms(self, radiusb0:Np2D) -> Np1D:
"""High fidelity system RHS norms."""
# 'ij,jk,ik->k', resbb, radiusb0, radiusb0.conj()
b0resb0 = np.sum(self.trainedModel.data.resbb.dot(radiusb0)
* radiusb0.conj(), axis = 0)
RHSnorms = np.power(np.abs(b0resb0), .5)
return RHSnorms
def _errorEstimatorBare(self) -> Np1D:
"""Bare residual-based error estimator."""
self.assembleReducedResidualGramian(self.trainedModel.data.projMat)
pDom = self.trainedModel.data.P[:, -1]
LL = pDom.conj().dot(self.trainedModel.data.gramian.dot(pDom))
Adiag = self.As[0].diagonal()
LL = ((self.scaleFactor * np.linalg.norm(Adiag)) ** 2. * LL
/ np.size(Adiag))
- scalingDom = polydomcoeff[self.polybasis](len(self.mus) - 1)
+ scalingDom = polydomcoeff(len(self.mus) - 1, self.polybasis)
return scalingDom * np.power(np.abs(LL), .5)
def _errorEstimatorBasic(self, muTest:complex, ratioTest:complex) -> Np1D:
"""Basic residual-based error estimator."""
resmu = self.HFEngine.residual(self.trainedModel.getApprox(muTest),
muTest, self.homogeneized)
return np.abs(self.estimatorNormEngine.norm(resmu) / ratioTest)
def _errorEstimatorExact(self, muRTrain:Np1D, vanderBase:Np2D) -> Np1D:
"""Exact residual-based error estimator."""
nAs = self.HFEngine.nAs - 1
nbs = max(self.HFEngine.nbs - 1, nAs * self.homogeneized)
delta = len(self.mus) - len(self.trainedModel.data.Q)
nbsEff = max(0, nbs - delta)
momentQ = np.zeros(nbsEff, dtype = np.complex)
momentQu = np.zeros((len(self.mus), nAs), dtype = np.complex)
radiusbTen = np.zeros((nbsEff, nbsEff, vanderBase.shape[1]),
dtype = np.complex)
radiusATen = np.zeros((nAs, nAs, vanderBase.shape[1]),
dtype = np.complex)
if nbsEff > 0:
momentQ[0] = self.trainedModel.data.Q[-1]
radiusbTen[0, :, :] = vanderBase[: nbsEff, :]
momentQu[:, 0] = self.trainedModel.data.P[:, -1]
radiusATen[0, :, :] = vanderBase[: nAs, :]
Qvals = self.trainedModel.getQVal(self.mus)
for k in range(1, max(nAs, nbs * (nbsEff > 0))):
Qvals = Qvals * muRTrain
if k > delta and k < nbs:
momentQ[k - delta] = self._fitinv.dot(Qvals)
radiusbTen[k - delta, k :, :] = (
radiusbTen[0, : delta - k, :])
if k < nAs:
momentQu[:, k] = Qvals * self._fitinv
radiusATen[k, k :, :] = radiusATen[0, : - k, :]
if self.POD and nAs > 1:
momentQu[:, 1 :] = self.samplingEngine.RPOD.dot(
momentQu[:, 1 :])
radiusA = np.tensordot(momentQu, radiusATen, 1)
if nbsEff > 0:
radiusb = np.tensordot(momentQ, radiusbTen, 1)
# 'ij,jk,ik->k', resbb, radiusb, radiusb.conj()
ff = np.sum(self.trainedModel.data.resbb[delta + 1 :, delta + 1 :]\
.dot(radiusb) * radiusb.conj(), axis = 0)
# 'ijk,jkl,il->l', resAb, radiusA, radiusb.conj()
Lf = np.sum(np.tensordot(
self.trainedModel.data.resAb[delta :, :, :], radiusA, 2)
* radiusb.conj(), axis = 0)
else:
ff, Lf = 0., 0.
# 'ijkl,klt,ijt->t', resAA, radiusA, radiusA.conj()
LL = np.sum(np.tensordot(self.trainedModel.data.resAA, radiusA, 2)
* radiusA.conj(), axis = (0, 1))
- scalingDom = polydomcoeff[self.polybasis](len(self.mus) - 1)
+ scalingDom = polydomcoeff(len(self.mus) - 1, self.polybasis)
return scalingDom * np.power(np.abs(ff - 2. * np.real(Lf) + LL), .5)
def errorEstimator(self, mus:Np1D) -> Np1D:
"""Standard residual-based error estimator."""
self.setupApprox()
if (np.any(np.isnan(self.trainedModel.data.P[:, -1]))
or np.any(np.isinf(self.trainedModel.data.P[:, -1]))):
err = np.empty(len(mus))
err[:] = np.inf
return err
nAs = self.HFEngine.nAs - 1
nbs = max(self.HFEngine.nbs - 1, nAs * self.homogeneized)
muRTest = [x(0) for x in self.centerNormalize(mus)]
muRTrain = [x(0) for x in self.centerNormalize(self.mus)]
self.assembleReducedResidualBlocks(kind = self.errorEstimatorKind)
samplingRatio = self._errorSamplingRatio(mus, muRTest, muRTrain)
vanderBase = np.polynomial.polynomial.polyvander(muRTest,
max(nAs, nbs)).T
RHSnorms = self._RHSNorms(vanderBase[: nbs + 1, :])
if self.errorEstimatorKind == "BARE":
jOpt = self._errorEstimatorBare()
elif self.errorEstimatorKind == "BASIC":
idx_muTestSample = np.argmax(samplingRatio)
muTestSample = mus[idx_muTestSample]
samplingRatioTestSample = samplingRatio[idx_muTestSample]
jOpt = self._errorEstimatorBasic(muTestSample,
samplingRatioTestSample)
else: #if self.errorEstimatorKind == "EXACT":
jOpt = self._errorEstimatorExact(muRTrain, vanderBase[: -1, :])
return jOpt * samplingRatio / RHSnorms
def _setupDenominator(self):
"""Compute Pade' denominator."""
if self.verbosity >= 7:
verbosityDepth("INIT", "Starting computation of denominator.",
timestamp = self.timestamp)
S = len(self.mus)
- TS = polyvander[self.polybasis](self.centerNormalize(self.mus),
- S - 1).T
+ TS = polyvander(self.centerNormalize(self.mus), S - 1,
+ self.polybasis).T
RHS = np.zeros(S)
RHS[-1] = 1.
fitOut = customFit(TS, RHS, full = True, rcond = self.interpRcond)
if self.verbosity >= 2:
condfit = fitOut[1][2][0] / fitOut[1][2][-1]
verbosityDepth("MAIN", ("Fitting {} samples with degree {} "
"through {}... Conditioning of system: "
"{:.4e}.").format(S, S - 1,
- polyfitname[self.polybasis],
+ polyfitname(self.polybasis),
condfit),
timestamp = self.timestamp)
if fitOut[1][1] < S:
RROMPyWarning(("Polyfit is poorly conditioned. Starting "
"preemptive termination of computation of "
"approximant."))
Q = np.empty(max(self.N, 0) + 1, dtype = np.complex)
P = np.empty((len(self.mus), max(self.M, 0) + 1),
dtype = np.complex)
Q[:] = np.nan
P[:] = np.nan
self.trainedModel.data.Q = copy(Q)
self.trainedModel.data.P = copy(P)
self.trainedModel.data.approxParameters = copy(
self.approxParameters)
if self.verbosity >= 7:
verbosityDepth("DEL", "Aborting computation of denominator.",
timestamp = self.timestamp)
return
self._fitinv = fitOut[0]
while self.N > 0:
Ghalf = (TS[: self.N + 1, :] * self._fitinv).T
if self.POD:
self.Ghalf = self.samplingEngine.RPOD.dot(Ghalf)
ev, eV = self.findeveVGQR(2)
else:
self.Ghalf = self.samplingEngine.samples.dot(Ghalf)
ev, eV = self.findeveVGQR(2)
Nstable = np.sum(np.abs(ev) >= self.robustTol * np.linalg.norm(ev))
if self.N <= Nstable: break
if self.verbosity >= 2:
verbosityDepth("MAIN", ("Smallest {} eigenvalues below "
"tolerance. Reducing N to {}.")\
.format(self.N - Nstable + 1, Nstable),
timestamp = self.timestamp)
self._N = Nstable
if self.N <= 0:
self._N = 0
eV = np.ones((1, 1))
if self.verbosity >= 7:
verbosityDepth("DEL", "Done computing denominator.",
timestamp = self.timestamp)
return eV[:, 0]
def _setupNumerator(self):
"""Compute Pade' numerator."""
if self.verbosity >= 7:
verbosityDepth("INIT", "Starting computation of numerator.",
timestamp = self.timestamp)
Qevaldiag = np.diag(self.trainedModel.getQVal(self.mus))
verb = self.trainedModel.verbosity
self.trainedModel.verbosity = 0
mus_un, idx_un, cnt_un = np.unique(self.mus, return_inverse = True,
return_counts = True)
for j in range(len(mus_un)):
if cnt_un[j] > 1:
Q_loc = np.zeros((cnt_un[j], cnt_un[j]), dtype = np.complex)
for der in range(1, cnt_un[j]):
Qderj = (self.trainedModel.getQVal(mus_un[j], der)
/ fact(der))
Q_loc = Q_loc + Qderj * np.diag(np.ones(cnt_un[j] - der),
k = - der)
I = idx_un == j
I = np.arange(len(self.mus))[I]
I = np.ix_(I, I)
Qevaldiag[I] = Qevaldiag[I] + Q_loc
self.trainedModel.verbosity = verb
while self.M >= 0:
- fitVander = polyvander[self.polybasis](
- self.centerNormalize(self.mus), self.M)
+ fitVander = polyvander(self.centerNormalize(self.mus), self.M,
+ self.polybasis)
w = None
S = len(self.mus)
if self.M == S - 1: w = "AUTO"
fitOut = customFit(fitVander, Qevaldiag, full = True, w = w,
rcond = self.interpRcond)
if self.verbosity >= 2:
condfit = fitOut[1][2][0] / fitOut[1][2][-1]
verbosityDepth("MAIN", ("Fitting {} samples with degree {} "
"through {}... Conditioning of "
"system: {:.4e}.").format(
S, self.M,
- polyfitname[self.polybasis],
+ polyfitname(self.polybasis),
condfit),
timestamp = self.timestamp)
if fitOut[1][1] == self.M + 1:
P = fitOut[0].T
break
RROMPyWarning(("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
if self.M < 0:
raise RROMPyException(("Instability in computation of numerator. "
"Aborting."))
return np.atleast_2d(P)
def setupApprox(self, plotEst : bool = False):
"""
Compute Pade' interpolant.
SVD-based robust eigenvalue management.
"""
if self.checkComputedApprox():
return
if self.verbosity >= 5:
verbosityDepth("INIT", "Setting up {}.". format(self.name()),
timestamp = self.timestamp)
self.computeScaleFactor()
self.greedy(plotEst)
S = len(self.mus)
self._M = S - 1
self._N = S - 1
if self.Delta < 0:
self._M += self.Delta
else:
self._N -= self.Delta
if self.trainedModel is None:
self.trainedModel = tModel()
self.trainedModel.verbosity = self.verbosity
self.trainedModel.timestamp = self.timestamp
data = TrainedModelData(self.trainedModel.name(), self.mu0,
self.samplingEngine.samples,
self.HFEngine.rescalingExp)
data.polytype = self.polybasis
data.scaleFactor = self.scaleFactor
data.mus = copy(self.mus)
self.trainedModel.data = data
else:
self.trainedModel = self.trainedModel
self.trainedModel.data.projMat = copy(self.samplingEngine.samples)
self.trainedModel.data.mus = copy(self.mus)
if min(self.M, self.N) < 0:
if self.verbosity >= 5:
verbosityDepth("MAIN", "Minimal sample size not achieved.",
timestamp = self.timestamp)
Q = np.empty(max(self.N, 0) + 1, dtype = np.complex)
P = np.empty((len(self.mus), max(self.M, 0) + 1),
dtype = np.complex)
Q[:] = np.nan
P[:] = np.nan
self.trainedModel.data.Q = copy(Q)
self.trainedModel.data.P = copy(P)
self.trainedModel.data.approxParameters = copy(
self.approxParameters)
if self.verbosity >= 5:
verbosityDepth("DEL", "Aborting computation of approximant.",
timestamp = self.timestamp)
return
if self.N > 0:
Q = self._setupDenominator()
if Q is None:
if self.verbosity >= 5:
verbosityDepth("DEL",
"Aborting computation of approximant.",
timestamp = self.timestamp)
return
else:
Q = np.ones((1,), dtype = np.complex)
self.trainedModel.data.Q = copy(Q)
P = self._setupNumerator()
if self.POD:
P = self.samplingEngine.RPOD.dot(P)
self.trainedModel.data.P = copy(P)
if self.verbosity >= 7:
verbosityDepth("DEL", "Done computing numerator.",
timestamp = self.timestamp)
self.trainedModel.data.approxParameters = copy(self.approxParameters)
if self.verbosity >= 5:
verbosityDepth("DEL", "Done setting up approximant.",
timestamp = self.timestamp)
def assembleReducedResidualBlocks(self, kind : str = "EXACT"):
"""Build affine blocks of reduced linear system through projections."""
scaling = self.trainedModel.data.scaleFactor
self.assembleReducedResidualBlocksbb(self.bs, scaling)
if kind == "EXACT":
pMat = self.trainedModel.data.projMat
self.assembleReducedResidualBlocksAb(self.As, self.bs[1 :],
pMat, scaling)
self.assembleReducedResidualBlocksAA(self.As, pMat, scaling,
basic = (kind == "BASIC"))
diff --git a/rrompy/reduction_methods/trained_model/trained_model_pade.py b/rrompy/reduction_methods/trained_model/trained_model_pade.py
index 71a8738..0c8d60d 100644
--- a/rrompy/reduction_methods/trained_model/trained_model_pade.py
+++ b/rrompy/reduction_methods/trained_model/trained_model_pade.py
@@ -1,150 +1,149 @@
# 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 . import TrainedModel
-from rrompy.utilities.base.types import Np1D, paramList, sampList
+from rrompy.utilities.base.types import Np1D, List, paramList, sampList
from rrompy.utilities.base import verbosityDepth
-from rrompy.utilities.poly_fitting import polyvalder, polyroots
+from rrompy.utilities.poly_fitting.polynomial import polyval, polyroots
from rrompy.utilities.exception_manager import RROMPyAssert
from rrompy.parameter import checkParameterList
from rrompy.sampling import sampleList
__all__ = ['TrainedModelPade']
class TrainedModelPade(TrainedModel):
"""
ROM approximant evaluation for Pade' approximant.
Attributes:
Data: dictionary with all that can be pickled.
"""
def centerNormalize(self, mu:paramList, mu0 : float = None) -> float:
"""
Compute normalized parameter to be plugged into approximant.
Args:
mu: Parameter(s) 1.
mu0: Parameter(s) 2. If None, set to self.data.mu0.
Returns:
Normalized parameter.
"""
mu, wasPar = checkParameterList(mu, self.data.npar)
if mu0 is None: mu0 = self.data.mu0
rad = ((mu ** self.data.rescalingExp - mu0 ** self.data.rescalingExp)
/ self.data.scaleFactor)
if wasPar: rad = rad[0]
return rad
- def getPVal(self, mu:paramList, der : int = 0) -> sampList:
+ def getPVal(self, mu:paramList, der : List[int] = None) -> sampList:
"""
Evaluate Pade' numerator at arbitrary parameter.
Args:
mu: Target parameter.
der(optional): Derivatives to take before evaluation.
"""
mu, wasPar = checkParameterList(mu, self.data.npar)
if self.verbosity >= 10:
verbosityDepth("INIT", ("Evaluating numerator at mu = "
"{}.").format(mu),
timestamp = self.timestamp)
muCenter = self.centerNormalize(mu)
- p = sampleList([polyvalder[self.data.polytype](mC, self.data.P.T, der)
- for mC in muCenter])
+ p = sampleList(polyval(muCenter, self.data.P.T,
+ self.data.polytype, der))
if self.verbosity >= 10:
verbosityDepth("DEL", "Done evaluating numerator.",
timestamp = self.timestamp)
if wasPar: p = p[0]
return p
- def getQVal(self, mu:Np1D, der : int = 0, scl : float = 1.) -> Np1D:
+ def getQVal(self, mu:Np1D, der : List[int] = None,
+ scl : Np1D = None) -> Np1D:
"""
Evaluate Pade' denominator at arbitrary parameter.
Args:
mu: Target parameter.
der(optional): Derivatives to take before evaluation.
"""
mu, wasPar = checkParameterList(mu, self.data.npar)
if self.verbosity >= 10:
verbosityDepth("INIT", ("Evaluating denominator at mu = "
"{}.").format(mu),
timestamp = self.timestamp)
muCenter = self.centerNormalize(mu)
- q = np.array([polyvalder[self.data.polytype](mC, self.data.Q, der, scl)
- for mC in muCenter])
+ q = polyval(muCenter, self.data.Q, self.data.polytype, der, scl)
if self.verbosity >= 10:
verbosityDepth("DEL", "Done evaluating denominator.",
timestamp = self.timestamp)
if wasPar: q = q[0]
return q
def getApproxReduced(self, mu:paramList) -> sampList:
"""
Evaluate reduced representation of approximant at arbitrary parameter.
Args:
mu: Target parameter.
"""
mu, wasPar = checkParameterList(mu, self.data.npar)
if (not hasattr(self, "lastSolvedAppReduced")
or self.lastSolvedAppReduced != mu):
if self.verbosity >= 5:
verbosityDepth("INIT", ("Evaluating approximant at mu = "
"{}.").format(mu),
timestamp = self.timestamp)
self.uAppReduced = self.getPVal(mu) / self.getQVal(mu)
if self.verbosity >= 5:
verbosityDepth("DEL", "Done evaluating approximant.",
timestamp = self.timestamp)
self.lastSolvedAppReduced = mu
if wasPar: return self.uAppReduced[0]
return self.uAppReduced
def getPoles(self) -> Np1D:
"""
Obtain approximant poles.
Returns:
Numpy complex vector of poles.
"""
- RROMPyAssert(self.data.npar, 1, "Number of parameters")
return np.power(self.data.mu0(0) ** self.data.rescalingExp
+ self.data.scaleFactor
- * polyroots[self.data.polytype](self.data.Q),
+ * polyroots(self.data.Q, self.data.polytype),
1. / self.data.rescalingExp)
def getResidues(self) -> Np1D:
"""
Obtain approximant residues.
Returns:
Numpy matrix with residues as columns.
"""
RROMPyAssert(self.data.npar, 1, "Number of parameters")
pls = self.getPoles()
poles, _ = checkParameterList(pls)
print(self.data.projMat.dot(self.getPVal(poles).data).shape)
print(self.getQVal(poles, 1).shape)
res = (self.data.projMat.dot(self.getPVal(poles).data)
/ self.getQVal(poles, 1))
return pls, res
diff --git a/rrompy/utilities/poly_fitting/__init__.py b/rrompy/utilities/poly_fitting/__init__.py
index 31700b4..dc4d3ac 100644
--- a/rrompy/utilities/poly_fitting/__init__.py
+++ b/rrompy/utilities/poly_fitting/__init__.py
@@ -1,35 +1,25 @@
# 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 .custom_fit import customFit
-from .fit_utils import (polybases, polyval, polyder, polyvalder, polyvander,
- polyfitname, polyroots, polydomcoeff)
__all__ = [
- 'customFit',
- 'polybases',
- 'polyval',
- 'polyder',
- 'polyvalder',
- 'polyvander',
- 'polyfitname',
- 'polyroots',
- 'polydomcoeff'
+ 'customFit'
]
diff --git a/rrompy/utilities/poly_fitting/fit_utils.py b/rrompy/utilities/poly_fitting/fit_utils.py
deleted file mode 100644
index 0703b17..0000000
--- a/rrompy/utilities/poly_fitting/fit_utils.py
+++ /dev/null
@@ -1,90 +0,0 @@
-# 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 numpy import pi, polynomial as po
-from scipy.special import binom
-from rrompy.utilities.base.types import Np1D, Np2D
-from rrompy.parameter import parameter, parameterList
-
-__all__ = ['polybases', 'polyval', 'polyder', 'polyvalder', 'polyvander',
- 'polyfitname', 'polyroots', 'polydomcoeff']
-
-polybases = ["CHEBYSHEV", "LEGENDRE", "MONOMIAL"]
-
-polyval = {"CHEBYSHEV" : lambda x, c: po.chebyshev.chebval(flatten(x), c),
- "LEGENDRE" : lambda x, c: po.legendre.legval(flatten(x), c),
- "MONOMIAL" : lambda x, c: po.polynomial.polyval(flatten(x), c)}
-polyder = {"CHEBYSHEV" : po.chebyshev.chebder, "LEGENDRE" : po.legendre.legder,
- "MONOMIAL" : po.polynomial.polyder}
-polyvalder = {
- "CHEBYSHEV" : lambda x, c, m = 1, scl = 1.:
- po.chebyshev.chebval(flatten(x), polyder["CHEBYSHEV"](c, m, scl)),
- "LEGENDRE" : lambda x, c, m = 1, scl = 1.:
- po.legendre.legval(flatten(x), polyder["LEGENDRE"](c, m, scl)),
- "MONOMIAL" : lambda x, c, m = 1, scl = 1.:
- po.polynomial.polyval(flatten(x), polyder["MONOMIAL"](c, m, scl))}
-polyvander = {
- "CHEBYSHEV" : lambda x, deg, scl = 1.:
- polyvanderConfluence(po.chebyshev.chebvander, polyder["CHEBYSHEV"],
- flatten(x), deg, scl),
- "LEGENDRE" : lambda x, deg, scl = 1:
- polyvanderConfluence(po.legendre.legvander, polyder["LEGENDRE"],
- flatten(x), deg, scl),
- "MONOMIAL" : lambda x, deg, scl = 1:
- polyvanderConfluence(po.polynomial.polyvander, polyder["MONOMIAL"],
- flatten(x), deg, scl)}
-
-polyfitname = {"CHEBYSHEV" : "chebfit", "LEGENDRE" : "legfit",
- "MONOMIAL" : "polyfit"}
-polyroots = {"CHEBYSHEV" : po.chebyshev.chebroots,
- "LEGENDRE" : po.legendre.legroots,
- "MONOMIAL" : po.polynomial.polyroots}
-polydomcoeff = {"CHEBYSHEV" : lambda n: 2. ** (n - 1) if n > 0 else 1.,
- "LEGENDRE" : lambda n: (2. ** n * (pi * n) ** -.5 if n > 10
- else .5 ** n * binom(2 * n, n)),
- "MONOMIAL" : lambda n: 1.}
-
-def flatten(x):
- if hasattr(x, "flatten"):
- return x.flatten()
- return x
-
-def polyvanderConfluence(vander:callable, derivative:callable, x:Np1D, deg:int,
- scl : float = 1.) -> Np2D:
- """Compute Vandermonde matrix even in case of confluence."""
- x_un, idx_un, cnt_un = np.unique(x, return_inverse = True,
- return_counts = True)
- Van = vander(x, deg)
- der_max = np.max(cnt_un) - 1
- if der_max > 0:
- C_der = np.zeros((deg + 1, deg + 1), dtype = float)
- for j in range(deg + 1):
- ej = np.zeros(deg + 1)
- ej[j] = 1.
- j_der = derivative(ej, 1, scl)
- C_der[: len(j_der), j] = j_der
-
- for der in range(1, der_max + 1):
- # remove first occurrence of each node
- for i_un in np.nonzero(cnt_un > der - 1)[0]:
- idx_un[np.nonzero(idx_un == i_un)[0][0]] = -1
- idx_loc = np.nonzero(idx_un > -1)[0]
- Van[idx_loc, :] = Van[idx_loc, :].dot(C_der[:, :]) / der
- return Van
-
diff --git a/rrompy/utilities/poly_fitting/__init__.py b/rrompy/utilities/poly_fitting/polynomial/__init__.py
similarity index 76%
copy from rrompy/utilities/poly_fitting/__init__.py
copy to rrompy/utilities/poly_fitting/polynomial/__init__.py
index 31700b4..c2a210d 100644
--- a/rrompy/utilities/poly_fitting/__init__.py
+++ b/rrompy/utilities/poly_fitting/polynomial/__init__.py
@@ -1,35 +1,35 @@
# 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 .custom_fit import customFit
-from .fit_utils import (polybases, polyval, polyder, polyvalder, polyvander,
- polyfitname, polyroots, polydomcoeff)
+from .base import (polybases, polyfitname, polydomcoeff)
+from .der import polyder
+from .val import polyval
+from .vander import polyvander
+from .roots import polyroots
__all__ = [
- 'customFit',
'polybases',
- 'polyval',
+ 'polyfitname',
+ 'polydomcoeff',
'polyder',
- 'polyvalder',
+ 'polyval',
'polyvander',
- 'polyfitname',
- 'polyroots',
- 'polydomcoeff'
+ 'polyroots'
]
diff --git a/rrompy/utilities/poly_fitting/polynomial/base.py b/rrompy/utilities/poly_fitting/polynomial/base.py
new file mode 100644
index 0000000..02562a6
--- /dev/null
+++ b/rrompy/utilities/poly_fitting/polynomial/base.py
@@ -0,0 +1,61 @@
+# Copyright (C) 2018 by the RROMPy authors
+#
+# This file is part of RROMPy.
+#
+# RROMPy is free software: you can redistribute it and/or modify
+# it under the terms of the GNU Lesser General Public License as published by
+# the Free Software Foundation, either version 3 of the License, or
+# (at your option) any later version.
+#
+# RROMPy is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+# GNU Lesser General Public License for more details.
+#
+# You should have received a copy of the GNU Lesser General Public License
+# along with RROMPy. If not, see .
+#
+
+import numpy as np
+from scipy.special import binom
+from rrompy.utilities.exception_manager import RROMPyException
+
+__all__ = ['polybases', 'polyfitname', 'polydomcoeff']
+
+polybases = ["CHEBYSHEV", "LEGENDRE", "MONOMIAL"]
+
+def polyfitname(basis:str) -> str:
+ fitnames = {"CHEBYSHEV" : "chebfit", "LEGENDRE" : "legfit",
+ "MONOMIAL" : "polyfit"}
+ try:
+ return fitnames[basis.upper()]
+ except:
+ raise RROMPyException("Polynomial basis not recognized.")
+
+def polydomcoeff(n:int, basis:str) -> float:
+ basis = basis.upper()
+ if isinstance(n, (list, tuple, np.ndarray,)):
+ nv = np.array(n)
+ else:
+ nv = np.array([n])
+ if basis == "CHEBYSHEV":
+ x = np.ones_like(nv, dtype = float)
+ x[nv > 0] = np.power(2., nv[nv > 0] - 1)
+ elif basis == "LEGENDRE":
+ x = np.ones_like(nv, dtype = float)
+ x[nv > 10] = (np.power(2., nv[nv > 10])
+ * np.power(np.pi * nv[nv > 10], -.5))
+ x[nv <= 10] = (np.power(.5, nv[nv <= 10])
+ * binom(2 * nv[nv <= 10], nv[nv <= 10]))
+ elif basis == "MONOMIAL":
+ x = np.ones_like(nv, dtype = float)
+ else:
+ raise RROMPyException("Polynomial basis not recognized.")
+ if isinstance(n, (list,)):
+ return list(x)
+ if isinstance(n, (tuple,)):
+ return tuple(x)
+ if isinstance(n, (np.ndarray,)):
+ return x
+ return x[0]
+
diff --git a/rrompy/utilities/poly_fitting/polynomial/der.py b/rrompy/utilities/poly_fitting/polynomial/der.py
new file mode 100644
index 0000000..092d0ff
--- /dev/null
+++ b/rrompy/utilities/poly_fitting/polynomial/der.py
@@ -0,0 +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 .
+#
+
+import numpy as np
+from rrompy.utilities.base.types import Np1D, Np2D, List
+from rrompy.utilities.exception_manager import RROMPyException, RROMPyAssert
+
+__all__ = ['polyder']
+
+def polyder(c:Np2D, basis:str, m : List[int] = None,
+ scl : Np1D = None) -> Np2D:
+ c = np.array(c, ndmin = 1, copy = 1)
+ if m is not None:
+ if scl is None: scl = [1.] * c.ndim
+ if not isinstance(m, (list,tuple,np.ndarray,)): m = [m]
+ if not isinstance(scl, (list,tuple,np.ndarray,)): scl = [scl]
+ RROMPyAssert(c.ndim, len(m), "Number of parameters")
+ RROMPyAssert(c.ndim, len(scl), "Number of parameters")
+ try:
+ derbase = {"CHEBYSHEV" : np.polynomial.chebyshev.chebder,
+ "LEGENDRE" : np.polynomial.legendre.legder,
+ "MONOMIAL" : np.polynomial.polynomial.polyder
+ }[basis.upper()]
+ except:
+ raise RROMPyException("Polynomial basis not recognized.")
+ for j, (mj, sj) in enumerate(zip(m, scl)):
+ c = derbase(c, m = mj, scl = sj, axis = j)
+ return c
diff --git a/rrompy/utilities/poly_fitting/__init__.py b/rrompy/utilities/poly_fitting/polynomial/roots.py
similarity index 58%
copy from rrompy/utilities/poly_fitting/__init__.py
copy to rrompy/utilities/poly_fitting/polynomial/roots.py
index 31700b4..ef48d32 100644
--- a/rrompy/utilities/poly_fitting/__init__.py
+++ b/rrompy/utilities/poly_fitting/polynomial/roots.py
@@ -1,35 +1,32 @@
# 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 .custom_fit import customFit
-from .fit_utils import (polybases, polyval, polyder, polyvalder, polyvander,
- polyfitname, polyroots, polydomcoeff)
-
-__all__ = [
- 'customFit',
- 'polybases',
- 'polyval',
- 'polyder',
- 'polyvalder',
- 'polyvander',
- 'polyfitname',
- 'polyroots',
- 'polydomcoeff'
- ]
+from numpy import polynomial as po
+from rrompy.utilities.base.types import Np1D
+from rrompy.utilities.exception_manager import RROMPyException
+__all__ = ['polyroots']
+def polyroots(c:Np1D, basis:str) -> Np1D:
+ try:
+ rootsbase = {"CHEBYSHEV" : po.chebyshev.chebroots,
+ "LEGENDRE" : po.legendre.legroots,
+ "MONOMIAL" : po.polynomial.polyroots}[basis.upper()]
+ except:
+ raise RROMPyException("Polynomial basis not recognized.")
+ return rootsbase(c)
diff --git a/rrompy/utilities/poly_fitting/polynomial/val.py b/rrompy/utilities/poly_fitting/polynomial/val.py
new file mode 100644
index 0000000..f8cb52f
--- /dev/null
+++ b/rrompy/utilities/poly_fitting/polynomial/val.py
@@ -0,0 +1,44 @@
+# 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.poly_fitting.polynomial import polyder
+from rrompy.utilities.base.types import Np1D, Np2D, List, paramList
+from rrompy.parameter import checkParameterList
+from rrompy.utilities.exception_manager import RROMPyException
+
+__all__ = ['polyval']
+
+def polyval(x:paramList, c:Np2D, basis:str, m : List[int] = None,
+ scl : Np1D = None) -> Np2D:
+ c = polyder(c, basis, m = m, scl = scl)
+ x, wasPar = checkParameterList(x)
+ if x.shape[1] > c.ndim:
+ raise RROMPyException("Incompatible parameter number")
+ try:
+ polyvalbase = {"CHEBYSHEV" : np.polynomial.chebyshev.chebval,
+ "LEGENDRE" : np.polynomial.legendre.legval,
+ "MONOMIAL" : np.polynomial.polynomial.polyval
+ }[basis.upper()]
+ except:
+ raise RROMPyException("Polynomial basis not recognized.")
+ c = polyvalbase(x(0), c, tensor = True)
+ for d in range(1, x.shape[1]):
+ c = polyvalbase(x(d), c, tensor = False)
+ if wasPar: return c[..., 0]
+ return c
diff --git a/rrompy/utilities/poly_fitting/polynomial/vander.py b/rrompy/utilities/poly_fitting/polynomial/vander.py
new file mode 100644
index 0000000..3f14f55
--- /dev/null
+++ b/rrompy/utilities/poly_fitting/polynomial/vander.py
@@ -0,0 +1,69 @@
+# 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.poly_fitting.polynomial import polyder
+from rrompy.utilities.base.types import Np2D, List, paramList
+from rrompy.parameter import checkParameterList
+from rrompy.utilities.exception_manager import RROMPyException
+
+__all__ = ['polyvander']
+
+#def polyvanderconfluence(x:Np1D, deg:int, basis:str,
+# scl : float = None) -> Np2D:
+# """Compute Vandermonde matrix even in case of confluence."""
+## does not work with parameterList
+# x_un, idx_un, cnt_un = np.unique(x, return_inverse = True,
+# return_counts = True)
+# Van = polyvander(x, deg, basis)
+# der_max = np.max(cnt_un) - 1
+# if der_max > 0: # must have square-like structure
+# C_der = np.zeros((deg + 1, deg + 1), dtype = float)
+# for j in range(deg + 1):
+# ej = np.zeros(deg + 1)
+# ej[j] = 1.
+# j_der = polyder(ej, basis, 1, scl)
+# C_der[: len(j_der), j] = j_der
+# for der in range(1, der_max + 1):
+# # remove first occurrence of each node
+# for i_un in np.nonzero(cnt_un > der - 1)[0]:
+# idx_un[np.nonzero(idx_un == i_un)[0][0]] = -1
+# idx_loc = np.nonzero(idx_un > -1)[0]
+# Van[idx_loc, :] = Van[idx_loc, :].dot(C_der[:, :]) / der
+# return Van
+
+def polyvander(x:paramList, degs:List[int], basis:str) -> Np2D:
+ if not isinstance(degs, (list,tuple,np.ndarray,)): degs = [degs]
+ ideg = [int(d) for d in degs]
+ is_valid = [id == d and id >= 0 for id, d in zip(ideg, degs)]
+ dim = len(ideg)
+ if is_valid != [1] * dim:
+ raise RROMPyException("Degrees must be non-negative integers")
+ x, wasPar = checkParameterList(x, dim)
+ try:
+ vanderbase = {"CHEBYSHEV" : np.polynomial.chebyshev.chebvander,
+ "LEGENDRE" : np.polynomial.legendre.legvander,
+ "MONOMIAL" : np.polynomial.polynomial.polyvander
+ }[basis.upper()]
+ except:
+ raise RROMPyException("Polynomial basis not recognized.")
+
+ v = vanderbase(x(0), ideg[0])
+ for j, dj in enumerate(ideg[1:]):
+ v = v[..., None] * vanderbase(x(j + 1), dj)[..., None, :]
+ return v.reshape(v.shape[:-dim] + (-1,))
diff --git a/rrompy/utilities/poly_fitting/radial_basis/radial_basis_fitter.py b/rrompy/utilities/poly_fitting/radial_basis/radial_basis_fitter.py
index e6d3ed7..2a8a4a7 100644
--- a/rrompy/utilities/poly_fitting/radial_basis/radial_basis_fitter.py
+++ b/rrompy/utilities/poly_fitting/radial_basis/radial_basis_fitter.py
@@ -1,202 +1,230 @@
# 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, Np2D, List, ListAny, paramList
from rrompy.solver import Np2DLikeEye, normEngine
from rrompy.parameter import checkParameterList
from rrompy.utilities.exception_manager import RROMPyException, RROMPyAssert
__all__ = ['radialBasisFitter', 'radialGaussian', 'thinPlateSpline',
'multiQuadric']
def radialGaussian(r2):
return np.exp(- r2)
def thinPlateSpline(r2):
return .5 * r2 * np.log(np.finfo(float).eps + r2)
def multiQuadric(r2):
return np.power(r2 + 1, .5)
class radialBasisFitter:
"""HERE"""
allowedModes = ["PARAMETERS", "VALUES"]
def __init__(self, mus:paramList, basisFun : callable = radialGaussian,
- massMatrix : Np2D = None, mode : str = "PARAMETERS"):
- self.mode = mode
+ massMatrix : Np2D = None, mode : str = "PARAMETERS",
+ scl : float = 1.):
self.mus = mus
self.basisFun = basisFun
if massMatrix is None: massMatrix = normEngine(Np2DLikeEye())
self.massMatrix = massMatrix
+ self.mode = mode
+ self.scl = scl
@property
def d(self):
"""Number of parameters."""
return self.mus.shape[1]
@property
def n(self):
"""Number of parameter points."""
return len(self.mus)
@property
- def mode(self):
- """Value of mode. Its assignment may reset snapshots."""
- return self._mode
- @mode.setter
- def mode(self, mode):
- if mode.upper() not in self.allowedModes:
- raise RROMPyException("Fitting mode not recognized.")
- self._mode = mode.upper()
+ def basisFun(self):
+ """Value of basisFun. Its assignment resets all."""
+ return self._basisFun
+ @basisFun.setter
+ def basisFun(self, basisFun):
+ self.reset()
+ self._basisFun = basisFun
@property
def mus(self):
- """Value of mus. Its assignment may reset snapshots."""
+ """Value of mus. Its assignment resets all."""
return self._mus
@mus.setter
def mus(self, mus):
mus, _ = checkParameterList(mus)
self.reset()
self._mus = mus
+ @property
+ def massMatrix(self):
+ """Value of massMatrix. Its assignment resets all."""
+ return self._massMatrix
+ @massMatrix.setter
+ def massMatrix(self, massMatrix):
+ self.reset()
+ self._massMatrix = massMatrix
+
+ @property
+ def mode(self):
+ """Value of mode. Its assignment resets all."""
+ return self._mode
+ @mode.setter
+ def mode(self, mode):
+ self.reset()
+ self._mode = mode.upper()
+
+ @property
+ def scl(self):
+ """Value of scl. Its assignment resets all."""
+ return self._scl
+ @scl.setter
+ def scl(self, scl):
+ self.reset()
+ self._scl = scl
+
def reset(self):
self.vander = None
self.offDiag = None
self.offDiagT = None
self.matrixInv = None
self.probeParameters = None
self.probeValues = None
def buildMatrixBlocks(self):
if self.offDiag is None:
self.reset()
self.offDiagT = np.array([[1] + list(x.data)
for x in self.mus.data])
self.offDiag = self.offDiagT.T
muDiff = np.empty((self.d, self.n * (self.n - 1) // 2 + 1),
dtype = self.mus.dtype)
muDiff[:, 0] = 0.
idxInv = np.zeros(self.n ** 2, dtype = int)
for j in range(self.n):
idx = j * (self.n - 1) - j * (j + 1) // 2
for i in range(j + 1, self.n):
muDiff[:, idx + i] = (self.offDiag[1:, j]
- self.offDiag[1:, i])
idxInv[j * self.n + i] = idx + i
idxInv[i * self.n + j] = idx + i
- self.vander = self.basisFun(np.power(self.massMatrix.norm(muDiff),
- 2.))[idxInv]
+ self.vander = self.basisFun(np.power(self.scl *
+ self.massMatrix.norm(muDiff), 2.))[idxInv]
self.vander = self.vander.reshape((self.n, -1))
self.vanderProj = self.offDiag.dot(self.vander.dot(self.offDiag.T))
def buildMatrixInvBlocks(self):
if self.matrixInv is None:
self.buildMatrixBlocks()
vanderInv = np.linalg.inv(self.vander)
vanderProjInv = np.linalg.solve(self.vanderProj,
self.offDiag.dot(vanderInv))
self.matrixInv = np.empty((self.n + self.d + 1, self.n),
dtype = vanderProjInv.dtype)
self.matrixInv[self.n :, :] = vanderProjInv
self.matrixInv[: self.n, :] = vanderInv - vanderInv.dot(
self.offDiagT.dot(vanderProjInv))
def setupProbeParameters(self, mu:paramList) -> List[Np1D]:
mu, _ = checkParameterList(mu, self.d)
self.buildMatrixInvBlocks()
self.probeParameters = [None] * len(mu)
for j, m in enumerate(mu):
probe = np.ones(self.n + self.d + 1, dtype = m.dtype)
probe[self.n + 1 :] = m.data
mDiff = (probe[self.n + 1:] - self.offDiagT[:, 1:]).T
- probe[: self.n] = self.basisFun(np.power(
+ probe[: self.n] = self.basisFun(np.power(self.scl *
self.massMatrix.norm(mDiff), 2.))
self.probeParameters[j] = probe.dot(self.matrixInv)
def setupProbeValues(self, vals:ListAny) -> ListAny:
RROMPyAssert(len(vals), self.n, "Number of samples")
self.buildMatrixInvBlocks()
if isinstance(vals, (np.ndarray,)):
self.probeValues = np.tensordot(self.matrixInv, vals,
axes = ([-1], [0]))
else:
self.probeValues = [None] * (self.n + self.d + 1)
for j in range(self.n + self.d + 1):
self.probeValues[j] = self.matrixInv[j, 0] * vals[0]
for i in range(1, self.n):
self.probeValues[j] += self.matrixInv[j, i] * vals[i]
def interpolateParameters(self, vals:ListAny) -> ListAny:
if self.probeParameters is None:
raise RROMPyException(("Parameter probe must be set up before "
"interpolation."))
RROMPyAssert(len(vals), self.n, "Number of samples")
interpolated = [None] * len(self.probeParameters)
if isinstance(vals, (np.ndarray,)):
if vals.ndim == 1:
for j, pUp in enumerate(self.probeParameters):
interpolated[j] = pUp.dot(vals)
else:
for j, pUp in enumerate(self.probeParameters):
interpolated[j] = np.tensordot(pUp, vals,
axes = ([0], [0]))
else:
for j, pUp in enumerate(self.probeParameters):
interpolated[j] = self.probeParameters[j][0] * vals[0]
for i in range(1, self.n):
interpolated[j] += self.probeParameters[j][i] * vals[i]
return interpolated
def interpolateValues(self, mu:paramList) -> ListAny:
if self.probeValues is None:
raise RROMPyException(("Value probe must be set up before "
"interpolation."))
mu, _ = checkParameterList(mu, self.d)
probeLs = [None] * len(mu)
for j, m in enumerate(mu):
probeLs[j] = np.ones(self.n + self.d + 1, dtype = m.dtype)
probeLs[j][self.n + 1 :] = m.data
mDiff = (probeLs[j][self.n + 1:] - self.offDiagT[:, 1:]).T
- probeLs[j][: self.n] = self.basisFun(np.power(
+ probeLs[j][: self.n] = self.basisFun(np.power(self.scl *
self.massMatrix.norm(mDiff), 2.))
interpolated = [None] * len(mu)
if isinstance(self.probeValues, (np.ndarray,)):
if self.probeValues.ndim == 1:
for j, pL in enumerate(probeLs):
interpolated[j] = pL.dot(self.probeValues)
else:
for j, pL in enumerate(probeLs):
interpolated[j] = np.tensordot(pL, self.probeValues,
axes = ([-1], [0]))
else:
for j, pL in enumerate(probeLs):
interpolated[j] = pL[0] * self.probeValues[0]
for i in range(1, self.n + self.d + 1):
interpolated[j] += pL[i] * self.probeValues[i]
return interpolated
def interpolate(self, x) -> ListAny:
if self.mode == "PARAMETERS":
return self.interpolateParameters(x)
if self.mode == "VALUES":
return self.interpolateValues(x)
raise RROMPyException("Not implemented")
\ No newline at end of file
diff --git a/tests/test_1_utilities/fitting.py b/tests/test_1_utilities/fitting.py
index cb9a1b7..cc188fc 100644
--- a/tests/test_1_utilities/fitting.py
+++ b/tests/test_1_utilities/fitting.py
@@ -1,98 +1,91 @@
# 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 pytest
import numpy as np
-from rrompy.utilities.poly_fitting import (customFit, polybases, polyval,
- polyvalder, polyvander, polyfitname, polyroots, polydomcoeff)
+from rrompy.utilities.poly_fitting import customFit
+from rrompy.utilities.poly_fitting.polynomial import (polybases, polyfitname,
+ polydomcoeff, polyval,
+ polyroots, polyvander)
def test_chebyshev():
assert "CHEBYSHEV" in polybases
- val = polyval["CHEBYSHEV"]
- valder = polyvalder["CHEBYSHEV"]
- vander = polyvander["CHEBYSHEV"]
- fitname = polyfitname["CHEBYSHEV"]
- roots = polyroots["CHEBYSHEV"]
- domcoeff = polydomcoeff["CHEBYSHEV"]
+ fitname = polyfitname("CHEBYSHEV")
+ domcoeff = polydomcoeff(5, "CHEBYSHEV")
assert fitname == "chebfit"
- assert np.isclose(domcoeff(5), 16, rtol = 1e-5)
- assert np.allclose(roots((-1, 1, -1, 1)), np.array([-.5, 0., 1.]),
- rtol = 1e-5)
- Phi = vander(np.arange(5), 4)
+ assert np.isclose(domcoeff, 16, rtol = 1e-5)
+ assert np.allclose(polyroots((-1, 1, -1, 1), "CHEBYSHEV"),
+ np.array([-.5, 0., 1.]), rtol = 1e-5)
+ Phi = polyvander(np.arange(5), 4, "CHEBYSHEV")
y = 2. * np.arange(5)
cFit = customFit(Phi, y)
assert np.allclose(cFit, [0, 2, 0, 0, 0], rtol = 1e-5)
- assert np.allclose(val(np.arange(5), cFit), y, rtol = 1e-5)
- assert np.allclose(valder(np.arange(5), cFit), 2. * np.ones(5),
+ assert np.allclose(polyval(np.arange(5), cFit, "CHEBYSHEV"), y,
rtol = 1e-5)
+ assert np.allclose(polyval(np.arange(5), cFit, "CHEBYSHEV", m = 1),
+ 2. * np.ones(5), rtol = 1e-5)
def test_legendre():
assert "LEGENDRE" in polybases
- val = polyval["LEGENDRE"]
- valder = polyvalder["LEGENDRE"]
- vander = polyvander["LEGENDRE"]
- fitname = polyfitname["LEGENDRE"]
- roots = polyroots["LEGENDRE"]
- domcoeff = polydomcoeff["LEGENDRE"]
+ fitname = polyfitname("LEGENDRE")
+ domcoeff = polydomcoeff([0, 5], "LEGENDRE")
assert fitname == "legfit"
- assert np.isclose(domcoeff(5), 63. / 8, rtol = 1e-5)
- assert np.allclose(roots((1, 2, 3, 4)),
+ assert np.allclose(domcoeff, [1., 63. / 8], rtol = 1e-5)
+ assert np.allclose(polyroots((1, 2, 3, 4), "LEGENDRE"),
np.array([-0.85099543, -0.11407192, 0.51506735]),
rtol = 1e-5)
- Phi = vander(np.arange(5), 4)
+ Phi = polyvander(np.arange(5), 4, "LEGENDRE")
y = 2. * np.arange(5)
cFit = customFit(Phi, y)
assert np.allclose(cFit, [0, 2, 0, 0, 0], rtol = 1e-5)
- assert np.allclose(val(np.arange(5), cFit), y, rtol = 1e-5)
- assert np.allclose(valder(np.arange(5), cFit), 2. * np.ones(5),
- rtol = 1e-5)
+ assert np.allclose(polyval(np.arange(5), cFit, "LEGENDRE"), y, rtol = 1e-5)
+ assert np.allclose(polyval(np.arange(5), cFit, "LEGENDRE", m = 1),
+ 2. * np.ones(5), rtol = 1e-5)
+
def test_monomial():
assert "MONOMIAL" in polybases
- val = polyval["MONOMIAL"]
- valder = polyvalder["MONOMIAL"]
- vander = polyvander["MONOMIAL"]
- fitname = polyfitname["MONOMIAL"]
- roots = polyroots["MONOMIAL"]
- domcoeff = polydomcoeff["MONOMIAL"]
+ fitname = polyfitname("MONOMIAL")
+ domcoeff = polydomcoeff(5, "MONOMIAL")
assert fitname == "polyfit"
- assert np.isclose(domcoeff(5), 1., rtol = 1e-5)
- assert np.allclose(roots([0.+0.j, 1.+0.j, 0.+0.j, 1.+0.j]),
+ assert np.isclose(domcoeff, 1., rtol = 1e-5)
+ assert np.allclose(polyroots([0.+0.j, 1.+0.j, 0.+0.j, 1.+0.j], "MONOMIAL"),
np.array([0., 1.j, -1.j]), rtol = 1e-5)
- Phi = vander(np.arange(5), 4)
+ Phi = polyvander(np.arange(5), 4, "MONOMIAL")
y = 2. * np.arange(5)
cFit = customFit(Phi, y)
assert np.allclose(cFit, [0, 2, 0, 0, 0], rtol = 1e-5)
- assert np.allclose(val(np.arange(5), cFit), y, rtol = 1e-5)
- assert np.allclose(valder(np.arange(5), cFit), 2. * np.ones(5),
- rtol = 1e-5)
+ assert np.allclose(polyval(np.arange(5), cFit, "MONOMIAL"), y, rtol = 1e-5)
+ assert np.allclose(polyval(np.arange(5), cFit, "MONOMIAL", m = 1),
+ 2. * np.ones(5), rtol = 1e-5)
+@pytest.mark.xfail
def test_cheb_confluence():
- val = polyval["CHEBYSHEV"]
- vander = polyvander["CHEBYSHEV"]
- valder = polyvalder["CHEBYSHEV"]
x = np.arange(5)
x[3] = 0
- Phi = vander(x, 4)
+ Phi = polyvander(x, 4, "CHEBYSHEV")
y = 2. * x
y[3] = 2.
cFit = customFit(Phi, y)
mask = np.arange(len(x)) != 3
assert np.allclose(cFit, [0, 2, 0, 0, 0], rtol = 1e-5)
- assert np.allclose(val(x[mask], cFit), y[mask], rtol = 1e-5)
- assert np.isclose(valder(x[~mask], cFit), y[~mask], rtol = 1e-5)
+ assert np.allclose(polyval(x[mask], cFit, "CHEBYSHEV"), y[mask],
+ rtol = 1e-5)
+ assert np.allclose(polyval(x[~mask], cFit, "CHEBYSHEV", m = 1), y[~mask],
+ rtol = 1e-5)
diff --git a/tests/test_3_reduction_methods/rational_interpolant.py b/tests/test_3_reduction_methods/rational_interpolant.py
index 6deef18..f62d6ec 100644
--- a/tests/test_3_reduction_methods/rational_interpolant.py
+++ b/tests/test_3_reduction_methods/rational_interpolant.py
@@ -1,69 +1,70 @@
# 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 pytest
import numpy as np
from matrix_fft import matrixFFT
from rrompy.reduction_methods.distributed import RationalInterpolant as RI
from rrompy.parameter.parameter_sampling import (QuadratureSampler as QS,
ManualSampler as MS)
def test_monomials(capsys):
mu = 1.5
solver = matrixFFT()
params = {"POD": False, "M": 9, "N": 9, "S": 10, "robustTol": 1e-6,
"interpRcond": 1e-3, "polybasis": "MONOMIAL",
"sampler": QS([1.5, 6.5], "UNIFORM")}
approx = RI(solver, 4., params, verbosity = 0)
approx.setupApprox()
out, err = capsys.readouterr()
assert (("poorly conditioned.\nReducing N from 9 to" in out)
and ("eigenvalues below tolerance. Reducing N from" in out))
assert len(err) == 0
assert np.isclose(approx.normErr(mu), .00773727, rtol = 1e-3)
def test_well_cond():
mu = 1.5
solver = matrixFFT()
params = {"POD": True, "M": 9, "N": 9, "S": 10, "robustTol": 1e-14,
"interpRcond": 1e-10, "polybasis": "CHEBYSHEV",
"sampler": QS([1., 7.], "CHEBYSHEV")}
approx = RI(solver, 4., params, verbosity = 0)
approx.setupApprox()
for mu in approx.mus:
assert np.isclose(approx.normErr(mu), 0., atol = 1e-8)
poles = approx.getPoles()
for lambda_ in np.arange(1, 8):
assert np.isclose(np.min(np.abs(poles - lambda_)), 0., atol = 1e-4)
+@pytest.mark.xfail
def test_hermite():
mu = 1.5
solver = matrixFFT()
sampler0 = QS([1., 7.], "CHEBYSHEV")
points = np.tile(sampler0.generatePoints(4)[0], 3)
params = {"POD": True, "M": 11, "N": 11, "S": 12, "polybasis": "CHEBYSHEV",
"sampler": MS([1., 7.], points = points)}
approx = RI(solver, 4., params, verbosity = 0)
approx.setupApprox()
for mu in approx.mus:
assert np.isclose(approx.normErr(mu), 0., atol = 1e-8)
poles = approx.getPoles()
for lambda_ in np.arange(1, 8):
assert np.isclose(np.min(np.abs(poles - lambda_)), 0., atol = 1e-4)
-#_multistorage