diff --git a/rational_interpolation_method.pdf b/rational_interpolation_method.pdf
new file mode 100644
index 0000000..b6ea8b0
Binary files /dev/null and b/rational_interpolation_method.pdf differ
diff --git a/rational_interpolation_method.tex b/rational_interpolation_method.tex
new file mode 100644
index 0000000..7a97969
--- /dev/null
+++ b/rational_interpolation_method.tex
@@ -0,0 +1,219 @@
+\documentclass[10pt,a4paper]{article}
+\usepackage[left=1in,right=1in,top=1in,bottom=1in]{geometry}
+\usepackage{amsmath}
+\usepackage{amsfonts}
+\usepackage{amssymb}
+\usepackage{hyperref}
+\usepackage{xcolor}
+
+\setlength{\parindent}{0pt}
+\newcommand{\code}[1]{{\color{blue}\texttt{#1}}}
+\newcommand{\footpath}[1]{\footnote{\path{#1}}}
+\newcommand{\N}{\mathbb{N}}
+\newcommand{\R}{\mathbb{R}}
+\newcommand{\C}{\mathbb{C}}
+\newcommand{\I}{\mathcal{I}}
+\DeclareMathOperator*{\argmin}{arg\,min}
+\newcommand{\inner}[2]{\left\langle#1,#2\right\rangle_V}
+\newcommand{\norm}[1]{\left\|#1\right\|_V}
+
+\title{\bf The RROMPy rational interpolation method}
+\author{D. Pradovera, CSQI, EPF Lausanne -- \texttt{davide.pradovera@epfl.ch}}
+\date{}
+\begin{document}
+\maketitle
+
+\section*{Introduction}
+This document provides an explanation for the numerical method provided by the class \code{Rational Interpolant}\footpath{./rrompy/reduction_methods/standard/rational_interpolant.py} and daughters, e.g. \code{Rational Interpolant Greedy}\footpath{./rrompy/reduction_methods/standard/greedy/rational_interpolant_greedy.py}, as well as most of the pivoted approximants\footpath{./rrompy/reduction_methods/pivoted/{,greedy/}rational_interpolant_*.py}.
+
+We restrict the discussion to the single-parameter case, and most of the focus will be dedicated to the impact of the \code{functionalSolve} parameter, whose allowed values are
+\begin{itemize}
+\item \code{NORM} (default): see \ref{sec:norm}; allows for derivative information, i.e. repeated sample points.
+\item \code{DOMINANT}: see \ref{sec:dominant}; allows for derivative information, i.e. repeated sample points.
+\item \code{BARYCENTRIC\_NORM}: see \ref{sec:barycentricnorm}; does not allow for a Least Squares (LS) approach.
+\item \code{BARYCENTRIC\_AVERAGE}: see \ref{sec:barycentricaverage}; does not allow for a Least Squares (LS) approach.
+\item \code{NODAL}: see \ref{sec:nodal}; iterative method.
+\end{itemize}
+
+The main reference throughout the present document is \cite{Pradovera}.
+
+\section{Aim of approximation}
+We seek an approximation of $u:\C\to V$, with $(V,\inner{\cdot}{\cdot})$ a complex\footnote{The inner product is linear (resp. conjugate linear) in the first (resp. second) argument: $\inner{\alpha v}{\beta w}=\alpha\overline{\beta}\inner{v}{w}$.} Hilbert space (with endowed norm $\norm{\cdot}$), of the form $\widehat{p}/\widehat{q}$, where $\widehat{p}:\C\to V$ and $\widehat{q}:\C\to\C$. For a given denominator $\widehat{q}$, the numerator $\widehat{p}$ is found by interpolation (possibly, LS or based on radial basis functions) of $\widehat{q}u$. Hence, here we focus on the computation of the denominator $\widehat{q}$.
+
+Other than the choice of target function $u$, the parameters which affect the computation of $\widehat{q}$ are:
+\begin{itemize}
+\item $\code{mus}\subset\C$ ($\{\mu_j\}_{j=1}^S$ below); for all \code{functionalSolve} values but \code{NORM} and \code{DOMINANT}, the $S$ points must be distinct.
+\item $\code{N}\in\N$ ($N$ below); for \code{BARYCENTRIC}, $N$ must equal $S-1$.
+\item $\code{polybasis0}\in\{\code{"CHEBYSHEV"}, \code{"LEGENDRE"}, \code{"MONOMIAL"}\}$; only for \code{NORM} and \code{DOMINANT}.
+\end{itemize}
+For simplicity, we will consider only the case of $S$ distinct sample points. One can deal with the case of confluent sample points by extending the standard (Lagrange) interpolation steps to Hermite-Lagrange ones.
+
+The main motivation behind the method involves the modified approximation problem
+\[u\approx\I^N\left(\Big(\big(\mu_j,\widehat{q}(\mu_j)u(\mu_j)\big)\Big)_{j=1}^S\right)\Big/\widehat{q},\]
+where $\widehat{q}:\C\to\C$ is a polynomial of degree $\leq N$, and $\I^N:(\C\times V)^S\to\mathbb{P}^N(\C;V)$ is a (LS) polynomial interpolation operator, which maps $S$ samples of a function (which lie in $V$) to a polynomial of degree $N.
#
from abc import abstractmethod
from copy import deepcopy as copy
import numpy as np
from collections.abc import Iterable
from matplotlib import pyplot as plt
from rrompy.reduction_methods.pivoted.generic_pivoted_approximant import (
GenericPivotedApproximantBase,
GenericPivotedApproximantNoMatch,
GenericPivotedApproximant)
from rrompy.reduction_methods.pivoted.gather_pivoted_approximant import (
gatherPivotedApproximant)
from rrompy.utilities.base.types import (Np1D, Np2D, Tuple, List, paramVal,
paramList, ListAny)
from rrompy.utilities.base import verbosityManager as vbMng
-from rrompy.utilities.numerical.point_matching import (pointMatching,
- chordalMetricAdjusted)
+from rrompy.utilities.numerical import pointMatching
+from rrompy.utilities.numerical.point_distances import chordalMetricAngleMatrix
from rrompy.utilities.exception_manager import (RROMPyException, RROMPyAssert,
RROMPyWarning)
from rrompy.parameter import emptyParameterList
from rrompy.utilities.parallel import (masterCore, indicesScatter,
arrayGatherv, isend)
__all__ = ['GenericPivotedGreedyApproximantNoMatch',
'GenericPivotedGreedyApproximant']
class GenericPivotedGreedyApproximantBase(GenericPivotedApproximantBase):
_allowedEstimatorKindsMarginal = ["LOOK_AHEAD", "LOOK_AHEAD_RECOVER",
"NONE"]
def __init__(self, *args, **kwargs):
self._preInit()
self._addParametersToList(["matchingWeightError",
"errorEstimatorKindMarginal",
"greedyTolMarginal", "maxIterMarginal"],
[0., "NONE", 1e-1, 1e2])
super().__init__(*args, **kwargs)
self._postInit()
@property
def scaleFactorDer(self):
"""Value of scaleFactorDer."""
if self._scaleFactorDer == "NONE": return 1.
if self._scaleFactorDer == "AUTO": return self._scaleFactorOldPivot
return self._scaleFactorDer
@scaleFactorDer.setter
def scaleFactorDer(self, scaleFactorDer):
if isinstance(scaleFactorDer, (str,)):
scaleFactorDer = scaleFactorDer.upper()
elif isinstance(scaleFactorDer, Iterable):
scaleFactorDer = list(scaleFactorDer)
self._scaleFactorDer = scaleFactorDer
self._approxParameters["scaleFactorDer"] = self._scaleFactorDer
@property
def samplerMarginal(self):
"""Value of samplerMarginal."""
return self._samplerMarginal
@samplerMarginal.setter
def samplerMarginal(self, samplerMarginal):
if 'refine' not in dir(samplerMarginal):
raise RROMPyException("Marginal sampler type not recognized.")
GenericPivotedApproximantBase.samplerMarginal.fset(self,
samplerMarginal)
@property
def errorEstimatorKindMarginal(self):
"""Value of errorEstimatorKindMarginal."""
return self._errorEstimatorKindMarginal
@errorEstimatorKindMarginal.setter
def errorEstimatorKindMarginal(self, errorEstimatorKindMarginal):
errorEstimatorKindMarginal = errorEstimatorKindMarginal.upper()
if errorEstimatorKindMarginal not in (
self._allowedEstimatorKindsMarginal):
RROMPyWarning(("Marginal error estimator kind not recognized. "
"Overriding to 'NONE'."))
errorEstimatorKindMarginal = "NONE"
self._errorEstimatorKindMarginal = errorEstimatorKindMarginal
self._approxParameters["errorEstimatorKindMarginal"] = (
self.errorEstimatorKindMarginal)
@property
def matchingWeightError(self):
"""Value of matchingWeightError."""
return self._matchingWeightError
@matchingWeightError.setter
def matchingWeightError(self, matchingWeightError):
self._matchingWeightError = matchingWeightError
self._approxParameters["matchingWeightError"] = (
self.matchingWeightError)
@property
def greedyTolMarginal(self):
"""Value of greedyTolMarginal."""
return self._greedyTolMarginal
@greedyTolMarginal.setter
def greedyTolMarginal(self, greedyTolMarginal):
if greedyTolMarginal < 0:
raise RROMPyException("greedyTolMarginal must be non-negative.")
if (hasattr(self, "_greedyTolMarginal")
and self.greedyTolMarginal is not None):
greedyTolMarginalold = self.greedyTolMarginal
else:
greedyTolMarginalold = -1
self._greedyTolMarginal = greedyTolMarginal
self._approxParameters["greedyTolMarginal"] = self.greedyTolMarginal
if greedyTolMarginalold != self.greedyTolMarginal:
self.resetSamples()
@property
def maxIterMarginal(self):
"""Value of maxIterMarginal."""
return self._maxIterMarginal
@maxIterMarginal.setter
def maxIterMarginal(self, maxIterMarginal):
if maxIterMarginal <= 0:
raise RROMPyException("maxIterMarginal must be positive.")
if (hasattr(self, "_maxIterMarginal")
and self.maxIterMarginal is not None):
maxIterMarginalold = self.maxIterMarginal
else:
maxIterMarginalold = -1
self._maxIterMarginal = maxIterMarginal
self._approxParameters["maxIterMarginal"] = self.maxIterMarginal
if maxIterMarginalold != self.maxIterMarginal:
self.resetSamples()
def resetSamples(self):
"""Reset samples."""
super().resetSamples()
if not hasattr(self, "_temporaryPivot"):
self._mus = emptyParameterList()
self._musMarginal = emptyParameterList()
if hasattr(self, "samplerMarginal"): self.samplerMarginal.reset()
if hasattr(self, "samplingEngine") and self.samplingEngine is not None:
self.samplingEngine.resetHistory()
def _getDistanceApp(self, polesEx:Np1D, resEx:Np2D, muTest:paramVal,
foci:Tuple[float, float], ground:float) -> float:
polesAp = self.trainedModel.interpolateMarginalPoles(muTest)[0]
if self.matchingWeightError != 0:
resAp = self.trainedModel.interpolateMarginalCoeffs(muTest)[0][
: len(polesAp), :]
resEx = self.trainedModel.data.projMat[:,
: resEx.shape[1]].dot(resEx.T)
resAp = self.trainedModel.data.projMat[:,
: resAp.shape[1]].dot(resAp.T)
else:
resAp = None
- dist = chordalMetricAdjusted(polesEx, polesAp,
- self.matchingWeightError, resEx, resAp,
- self.HFEngine, False)
+ dist = chordalMetricAngleMatrix(polesEx, polesAp,
+ self.matchingWeightError, resEx, resAp,
+ self.HFEngine, False)
pmR, pmC = pointMatching(dist)
return np.mean(dist[pmR, pmC])
def getErrorEstimatorMarginalLookAhead(self) -> Np1D:
if not hasattr(self.trainedModel, "_musMExcl"):
err = np.zeros(0)
err[:] = np.inf
self._musMarginalTestIdxs = np.zeros(0, dtype = int)
return err
self._musMarginalTestIdxs = np.array(self.trainedModel._idxExcl,
dtype = int)
idx, sizes = indicesScatter(len(self.trainedModel._musMExcl),
return_sizes = True)
err = []
if len(idx) > 0:
self.verbosity -= 35
self.trainedModel.verbosity -= 35
foci = self.samplerPivot.normalFoci()
ground = self.samplerPivot.groundPotential()
for j in idx:
muTest = self.trainedModel._musMExcl[j]
HITest = self.trainedModel._HIsExcl[j]
polesEx = HITest.poles
idxGood = np.logical_not(np.logical_or(np.isinf(polesEx),
np.isnan(polesEx)))
polesEx = polesEx[idxGood]
if self.matchingWeightError != 0:
resEx = HITest.coeffs[np.where(idxGood)[0]]
else:
resEx = None
if len(polesEx) == 0:
err += [0.]
continue
err += [self._getDistanceApp(polesEx, resEx, muTest,
foci, ground)]
self.verbosity += 35
self.trainedModel.verbosity += 35
return arrayGatherv(np.array(err), sizes)
def getErrorEstimatorMarginalNone(self) -> Np1D:
nErr = len(self.trainedModel.data.musMarginal)
self._musMarginalTestIdxs = np.arange(nErr)
return (1. + self.greedyTolMarginal) * np.ones(nErr)
def errorEstimatorMarginal(self, return_max : bool = False) -> Np1D:
vbMng(self.trainedModel, "INIT",
"Evaluating error estimator at mu = {}.".format(
self.trainedModel.data.musMarginal), 10)
if self.errorEstimatorKindMarginal == "NONE":
nErr = len(self.trainedModel.data.musMarginal)
self._musMarginalTestIdxs = np.arange(nErr)
err = (1. + self.greedyTolMarginal) * np.ones(nErr)
else:#if self.errorEstimatorKindMarginal[: 10] == "LOOK_AHEAD":
err = self.getErrorEstimatorMarginalLookAhead()
vbMng(self.trainedModel, "DEL", "Done evaluating error estimator", 10)
if not return_max: return err
idxMaxEst = np.where(err > self.greedyTolMarginal)[0]
maxErr = err[idxMaxEst]
if self.errorEstimatorKindMarginal == "NONE": maxErr = None
return err, idxMaxEst, maxErr
def plotEstimatorMarginal(self, est:Np1D, idxMax:List[int],
estMax:List[float]):
if self.errorEstimatorKindMarginal == "NONE": return
if (not (np.any(np.isnan(est)) or np.any(np.isinf(est)))
and masterCore() and hasattr(self.trainedModel, "_musMExcl")):
fig = plt.figure(figsize = plt.figaspect(1. / self.nparMarginal))
for jpar in range(self.nparMarginal):
ax = fig.add_subplot(1, self.nparMarginal, 1 + jpar)
musre = np.real(self.trainedModel._musMExcl)
if len(idxMax) > 0 and estMax is not None:
maxrej = musre[idxMax, jpar]
errCP = copy(est)
idx = np.delete(np.arange(self.nparMarginal), jpar)
while len(musre) > 0:
if self.nparMarginal == 1:
currIdx = np.arange(len(musre))
else:
currIdx = np.where(np.isclose(np.sum(
np.abs(musre[:, idx] - musre[0, idx]), 1), 0.))[0]
currIdxSorted = currIdx[np.argsort(musre[currIdx, jpar])]
ax.semilogy(musre[currIdxSorted, jpar],
errCP[currIdxSorted], 'k.-', linewidth = 1)
musre = np.delete(musre, currIdx, 0)
errCP = np.delete(errCP, currIdx)
ax.semilogy(self.musMarginal.re(jpar),
(self.greedyTolMarginal,) * len(self.musMarginal),
'*m')
if len(idxMax) > 0 and estMax is not None:
ax.semilogy(maxrej, estMax, 'xr')
ax.set_xlim(*list(self.samplerMarginal.lims.re(jpar)))
ax.grid()
plt.tight_layout()
plt.show()
def _addMarginalSample(self, mus:paramList):
mus = self.checkParameterListMarginal(mus)
if len(mus) == 0: return
self._nmusOld, nmus = len(self.musMarginal), len(mus)
if (hasattr(self, "trainedModel") and self.trainedModel is not None
and hasattr(self.trainedModel, "_musMExcl")):
self._nmusOld += len(self.trainedModel._musMExcl)
vbMng(self, "MAIN",
("Adding marginal sample point{} no. {}{} at {} to training "
"set.").format("s" * (nmus > 1), self._nmusOld + 1,
"--{}".format(self._nmusOld + nmus) * (nmus > 1),
mus), 3)
self.musMarginal.append(mus)
self.setupApproxPivoted(mus)
self._poleMatching()
del self._nmusOld
if (self.errorEstimatorKindMarginal[: 10] == "LOOK_AHEAD"
and not self.firstGreedyIterM):
ubRange = len(self.trainedModel.data.musMarginal)
if hasattr(self.trainedModel, "_idxExcl"):
shRange = len(self.trainedModel._musMExcl)
else:
shRange = 0
testIdxs = list(range(ubRange + shRange - len(mus),
ubRange + shRange))
for j in testIdxs[::-1]:
self.musMarginal.pop(j - shRange)
if hasattr(self.trainedModel, "_idxExcl"):
testIdxs = self.trainedModel._idxExcl + testIdxs
self._updateTrainedModelMarginalSamples(testIdxs)
self._finalizeMarginalization()
self._SMarginal = len(self.musMarginal)
self._approxParameters["SMarginal"] = self.SMarginal
self.trainedModel.data.approxParameters["SMarginal"] = self.SMarginal
def greedyNextSampleMarginal(self, muidx:List[int],
plotEst : str = "NONE") \
-> Tuple[Np1D, List[int], float, paramVal]:
RROMPyAssert(self._mode, message = "Cannot add greedy sample.")
muidx = self._musMarginalTestIdxs[muidx]
if (self.errorEstimatorKindMarginal[: 10] == "LOOK_AHEAD"
and not self.firstGreedyIterM):
if not hasattr(self.trainedModel, "_idxExcl"):
raise RROMPyException(("Sample index to be added not present "
"in trained model."))
testIdxs = copy(self.trainedModel._idxExcl)
skippedIdx = 0
for cj, j in enumerate(self.trainedModel._idxExcl):
if j in muidx:
testIdxs.pop(skippedIdx)
self.musMarginal.insert(self.trainedModel._musMExcl[cj],
j - skippedIdx)
else:
skippedIdx += 1
if len(self.trainedModel._idxExcl) < (len(muidx)
+ len(testIdxs)):
raise RROMPyException(("Sample index to be added not present "
"in trained model."))
self._updateTrainedModelMarginalSamples(testIdxs)
self._SMarginal = len(self.musMarginal)
self._approxParameters["SMarginal"] = self.SMarginal
self.trainedModel.data.approxParameters["SMarginal"] = (
self.SMarginal)
self.firstGreedyIterM = False
idxAdded = self.samplerMarginal.refine(muidx)[0]
self._addMarginalSample(self.samplerMarginal.points[idxAdded])
errorEstTest, muidx, maxErrorEst = self.errorEstimatorMarginal(True)
if plotEst == "ALL":
self.plotEstimatorMarginal(errorEstTest, muidx, maxErrorEst)
return (errorEstTest, muidx, maxErrorEst,
self.samplerMarginal.points[muidx])
def _preliminaryTrainingMarginal(self):
"""Initialize starting snapshots of solution map."""
RROMPyAssert(self._mode, message = "Cannot start greedy algorithm.")
if np.sum(self.samplingEngine.nsamples) > 0: return
self.resetSamples()
self._addMarginalSample(self.samplerMarginal.generatePoints(
self.SMarginal))
def _preSetupApproxPivoted(self, mus:paramList) \
-> Tuple[ListAny, ListAny, ListAny]:
self.computeScaleFactor()
if self.trainedModel is None:
self._setupTrainedModel(np.zeros((0, 0)))
self.trainedModel.data.Qs, self.trainedModel.data.Ps = [], []
self.trainedModel.data.Psupp = []
self._trainedModelOld = copy(self.trainedModel)
self._scaleFactorOldPivot = copy(self.scaleFactor)
self.scaleFactor = self.scaleFactorPivot
self._temporaryPivot = 1
self._musLoc = copy(self.mus)
idx, sizes = indicesScatter(len(mus), return_sizes = True)
emptyCores = np.where(np.logical_not(sizes))[0]
self.verbosity -= 15
return idx, sizes, emptyCores
def _postSetupApproxPivoted(self, mus:Np2D, pMat:Np2D, Ps:ListAny,
Qs:ListAny, sizes:ListAny):
self.scaleFactor = self._scaleFactorOldPivot
del self._scaleFactorOldPivot, self._temporaryPivot
pMat, Ps, Qs, mus, nsamples = gatherPivotedApproximant(pMat, Ps, Qs,
mus, sizes,
self.polybasis)
if len(self._musLoc) > 0:
self._mus = self.checkParameterList(self._musLoc)
self._mus.append(mus)
else:
self._mus = self.checkParameterList(mus)
self.trainedModel = self._trainedModelOld
del self._trainedModelOld
padLeft = self.trainedModel.data.projMat.shape[1]
suppNew = np.append(0, np.cumsum(nsamples))
self._setupTrainedModel(pMat, padLeft > 0)
self.trainedModel.data.Qs += Qs
self.trainedModel.data.Ps += Ps
self.trainedModel.data.Psupp += list(padLeft + suppNew[: -1])
self.trainedModel.data.approxParameters = copy(self.approxParameters)
self.verbosity += 15
def _localPivotedResult(self, pMat:Np2D, req:ListAny, emptyCores:ListAny,
mus:Np2D) -> Tuple[Np2D, ListAny, Np2D]:
if pMat is None:
mus = copy(self.samplingEngine.mus.data)
pMat = copy(self.samplingEngine.projectionMatrix)
if masterCore():
for dest in emptyCores:
req += [isend((len(pMat), pMat.dtype, mus.dtype),
dest = dest, tag = dest)]
else:
mus = np.vstack((mus, self.samplingEngine.mus.data))
pMat = np.hstack((pMat,
self.samplingEngine.projectionMatrix))
return pMat, req, mus
@abstractmethod
def setupApproxPivoted(self, mus:paramList) -> int:
if self.checkComputedApproxPivoted(): return -1
RROMPyAssert(self._mode, message = "Cannot setup approximant.")
vbMng(self, "INIT", "Setting up pivoted approximant.", 10)
self._preSetupApproxPivoted()
data = []
pass
self._postSetupApproxPivoted(mus, data)
vbMng(self, "DEL", "Done setting up pivoted approximant.", 10)
return 0
def setupApprox(self, plotEst : str = "NONE") -> int:
"""Compute greedy snapshots of solution map."""
if self.checkComputedApprox(): return -1
RROMPyAssert(self._mode, message = "Cannot start greedy algorithm.")
vbMng(self, "INIT", "Starting computation of snapshots.", 3)
max2ErrorEst, self.firstGreedyIterM = np.inf, True
self._preliminaryTrainingMarginal()
if self.errorEstimatorKindMarginal == "NONE":
muidx = []
else:#if self.errorEstimatorKindMarginal[: 10] == "LOOK_AHEAD":
muidx = np.arange(len(self.trainedModel.data.musMarginal))
self._musMarginalTestIdxs = np.array(muidx)
while self.firstGreedyIterM or (max2ErrorEst > self.greedyTolMarginal
and self.samplerMarginal.npoints < self.maxIterMarginal):
errorEstTest, muidx, maxErrorEst, mu = \
self.greedyNextSampleMarginal(muidx, plotEst)
if maxErrorEst is None:
max2ErrorEst = 1. + self.greedyTolMarginal
else:
if len(maxErrorEst) > 0:
max2ErrorEst = np.max(maxErrorEst)
else:
max2ErrorEst = np.max(errorEstTest)
vbMng(self, "MAIN", ("Uniform testing error estimate "
"{:.4e}.").format(max2ErrorEst), 3)
if plotEst == "LAST":
self.plotEstimatorMarginal(errorEstTest, muidx, maxErrorEst)
vbMng(self, "DEL", ("Done computing snapshots (final snapshot count: "
"{}).").format(len(self.mus)), 3)
if (self.errorEstimatorKindMarginal == "LOOK_AHEAD_RECOVER"
and hasattr(self.trainedModel, "_idxExcl")
and len(self.trainedModel._idxExcl) > 0):
vbMng(self, "INIT", "Recovering {} test models.".format(
len(self.trainedModel._idxExcl)), 7)
for j, mu in zip(self.trainedModel._idxExcl,
self.trainedModel._musMExcl):
self.musMarginal.insert(mu, j)
self._updateTrainedModelMarginalSamples()
self._finalizeMarginalization()
self._SMarginal = len(self.musMarginal)
self._approxParameters["SMarginal"] = self.SMarginal
self.trainedModel.data.approxParameters["SMarginal"] = (
self.SMarginal)
vbMng(self, "DEL", "Done recovering test models.", 7)
return 0
def checkComputedApproxPivoted(self) -> bool:
return (super().checkComputedApprox()
and len(self.musMarginal) == len(self.trainedModel.data.musMarginal))
class GenericPivotedGreedyApproximantNoMatch(
GenericPivotedGreedyApproximantBase,
GenericPivotedApproximantNoMatch):
"""
ROM pivoted greedy interpolant computation for parametric problems (without
pole matching) (ABSTRACT).
Args:
HFEngine: HF problem solver.
mu0(optional): Default parameter. Defaults to 0.
directionPivot(optional): Pivot components. Defaults to [0].
approxParameters(optional): Dictionary containing values for main
parameters of approximant. Recognized keys are:
- 'POD': kind of snapshots orthogonalization; allowed values
include 0, 1/2, and 1; defaults to 1, i.e. POD;
- 'scaleFactorDer': scaling factors for derivative computation;
defaults to 'AUTO';
- 'matchingWeightError': weight for pole matching optimization in
error estimation; defaults to 0;
- 'S': total number of pivot samples current approximant relies
upon;
- 'samplerPivot': pivot sample point generator;
- 'SMarginal': number of starting marginal samples;
- 'samplerMarginal': marginal sample point generator via sparse
grid;
- 'errorEstimatorKindMarginal': kind of marginal error estimator;
available values include 'LOOK_AHEAD', 'LOOK_AHEAD_RECOVER',
and 'NONE'; defaults to 'NONE';
- 'greedyTolMarginal': uniform error tolerance for marginal greedy
algorithm; defaults to 1e-1;
- 'maxIterMarginal': maximum number of marginal greedy steps;
defaults to 1e2;
- 'radialDirectionalWeightsMarginal': radial basis weights for
marginal interpolant; defaults to 1.
Defaults to empty dict.
approx_state(optional): Whether to approximate state. Defaults to
False.
verbosity(optional): Verbosity level. Defaults to 10.
Attributes:
HFEngine: HF problem solver.
mu0: Default parameter.
directionPivot: Pivot components.
mus: Array of snapshot parameters.
musMarginal: Array of marginal snapshot parameters.
approxParameters: Dictionary containing values for main parameters of
approximant. Recognized keys are in parameterList.
parameterListSoft: Recognized keys of soft approximant parameters:
- 'POD': kind of snapshots orthogonalization;
- 'scaleFactorDer': scaling factors for derivative computation;
- 'matchingWeightError': weight for pole matching optimization in
error estimation;
- 'errorEstimatorKindMarginal': kind of marginal error estimator;
- 'greedyTolMarginal': uniform error tolerance for marginal greedy
algorithm;
- 'maxIterMarginal': maximum number of marginal greedy steps;
- 'radialDirectionalWeightsMarginal': radial basis weights for
marginal interpolant.
parameterListCritical: Recognized keys of critical approximant
parameters:
- 'S': total number of pivot samples current approximant relies
upon;
- 'samplerPivot': pivot sample point generator;
- 'SMarginal': total number of marginal samples current approximant
relies upon;
- 'samplerMarginal': marginal sample point generator via sparse
grid.
approx_state: Whether to approximate state.
verbosity: Verbosity level.
POD: Kind of snapshots orthogonalization.
scaleFactorDer: Scaling factors for derivative computation.
matchingWeightError: Weight for pole matching optimization in error
estimation.
S: Total number of pivot samples current approximant relies upon.
samplerPivot: Pivot sample point generator.
SMarginal: Total number of marginal samples current approximant relies
upon.
samplerMarginal: Marginal sample point generator via sparse grid.
errorEstimatorKindMarginal: Kind of marginal error estimator.
greedyTolMarginal: Uniform error tolerance for marginal greedy
algorithm.
maxIterMarginal: Maximum number of marginal greedy steps.
radialDirectionalWeightsMarginal: Radial basis weights for marginal
interpolant.
muBounds: list of bounds for pivot parameter values.
muBoundsMarginal: list of bounds for marginal parameter values.
samplingEngine: Sampling engine.
uHF: High fidelity solution(s) with parameter(s) lastSolvedHF as
sampleList.
lastSolvedHF: Parameter(s) corresponding to last computed high fidelity
solution(s) as parameterList.
uApproxReduced: Reduced approximate solution(s) with parameter(s)
lastSolvedApprox as sampleList.
lastSolvedApproxReduced: Parameter(s) corresponding to last computed
reduced approximate solution(s) as parameterList.
uApprox: Approximate solution(s) with parameter(s) lastSolvedApprox as
sampleList.
lastSolvedApprox: Parameter(s) corresponding to last computed
approximate solution(s) as parameterList.
"""
def _poleMatching(self):
vbMng(self, "INIT", "Compressing poles.", 10)
self.trainedModel.initializeFromRational()
vbMng(self, "DEL", "Done compressing poles.", 10)
def _updateTrainedModelMarginalSamples(self, idx : ListAny = []):
self.trainedModel.updateEffectiveSamples(idx)
class GenericPivotedGreedyApproximant(GenericPivotedGreedyApproximantBase,
GenericPivotedApproximant):
"""
ROM pivoted greedy interpolant computation for parametric problems (with
pole matching) (ABSTRACT).
Args:
HFEngine: HF problem solver.
mu0(optional): Default parameter. Defaults to 0.
directionPivot(optional): Pivot components. Defaults to [0].
approxParameters(optional): Dictionary containing values for main
parameters of approximant. Recognized keys are:
- 'POD': kind of snapshots orthogonalization; allowed values
include 0, 1/2, and 1; defaults to 1, i.e. POD;
- 'scaleFactorDer': scaling factors for derivative computation;
defaults to 'AUTO';
- 'matchingWeight': weight for pole matching optimization; defaults
to 1;
- 'sharedRatio': required ratio of marginal points to share
resonance; defaults to 1.;
- 'matchingWeightError': weight for pole matching optimization in
error estimation; defaults to 0;
- 'S': total number of pivot samples current approximant relies
upon;
- 'samplerPivot': pivot sample point generator;
- 'SMarginal': number of starting marginal samples;
- 'samplerMarginal': marginal sample point generator via sparse
grid;
- 'errorEstimatorKindMarginal': kind of marginal error estimator;
available values include 'LOOK_AHEAD', 'LOOK_AHEAD_RECOVER',
and 'NONE'; defaults to 'NONE';
- 'polybasisMarginal': type of polynomial basis for marginal
interpolation; allowed values include 'MONOMIAL_*',
'CHEBYSHEV_*', 'LEGENDRE_*', 'NEARESTNEIGHBOR', and
'PIECEWISE_LINEAR_*'; defaults to 'MONOMIAL';
- 'paramsMarginal': dictionary of parameters for marginal
interpolation; include:
. 'MMarginal': degree of marginal interpolant; defaults to
'AUTO', i.e. maximum allowed; not for 'NEARESTNEIGHBOR' or
'PIECEWISE_LINEAR_*';
. 'nNeighborsMarginal': number of marginal nearest neighbors;
defaults to 1; only for 'NEARESTNEIGHBOR';
. 'polydegreetypeMarginal': type of polynomial degree for
marginal; defaults to 'TOTAL'; not for 'NEARESTNEIGHBOR' or
'PIECEWISE_LINEAR_*';
. 'interpRcondMarginal': tolerance for marginal interpolation;
defaults to None; not for 'NEARESTNEIGHBOR' or
'PIECEWISE_LINEAR_*';
. 'radialDirectionalWeightsMarginalAdapt': bounds for adaptive
rescaling of marginal radial basis weights; only for
radial basis.
- 'greedyTolMarginal': uniform error tolerance for marginal greedy
algorithm; defaults to 1e-1;
- 'maxIterMarginal': maximum number of marginal greedy steps;
defaults to 1e2;
- 'radialDirectionalWeightsMarginal': radial basis weights for
marginal interpolant; defaults to 1.
Defaults to empty dict.
approx_state(optional): Whether to approximate state. Defaults to
False.
verbosity(optional): Verbosity level. Defaults to 10.
Attributes:
HFEngine: HF problem solver.
mu0: Default parameter.
directionPivot: Pivot components.
mus: Array of snapshot parameters.
musMarginal: Array of marginal snapshot parameters.
approxParameters: Dictionary containing values for main parameters of
approximant. Recognized keys are in parameterList.
parameterListSoft: Recognized keys of soft approximant parameters:
- 'POD': kind of snapshots orthogonalization;
- 'scaleFactorDer': scaling factors for derivative computation;
- 'matchingWeight': weight for pole matching optimization;
- 'sharedRatio': required ratio of marginal points to share
resonance;
- 'matchingWeightError': weight for pole matching optimization in
error estimation;
- 'errorEstimatorKindMarginal': kind of marginal error estimator;
- 'polybasisMarginal': type of polynomial basis for marginal
interpolation;
- 'paramsMarginal': dictionary of parameters for marginal
interpolation; include:
. 'MMarginal': degree of marginal interpolant;
. 'nNeighborsMarginal': number of marginal nearest neighbors;
. 'polydegreetypeMarginal': type of polynomial degree for
marginal;
. 'interpRcondMarginal': tolerance for marginal interpolation;
. 'radialDirectionalWeightsMarginalAdapt': bounds for adaptive
rescaling of marginal radial basis weights.
- 'greedyTolMarginal': uniform error tolerance for marginal greedy
algorithm;
- 'maxIterMarginal': maximum number of marginal greedy steps;
- 'radialDirectionalWeightsMarginal': radial basis weights for
marginal interpolant.
parameterListCritical: Recognized keys of critical approximant
parameters:
- 'S': total number of pivot samples current approximant relies
upon;
- 'samplerPivot': pivot sample point generator;
- 'SMarginal': total number of marginal samples current approximant
relies upon;
- 'samplerMarginal': marginal sample point generator via sparse
grid.
approx_state: Whether to approximate state.
verbosity: Verbosity level.
POD: Kind of snapshots orthogonalization.
scaleFactorDer: Scaling factors for derivative computation.
matchingWeight: Weight for pole matching optimization.
sharedRatio: Required ratio of marginal points to share resonance.
matchingWeightError: Weight for pole matching optimization in error
estimation.
S: Total number of pivot samples current approximant relies upon.
samplerPivot: Pivot sample point generator.
SMarginal: Total number of marginal samples current approximant relies
upon.
samplerMarginal: Marginal sample point generator via sparse grid.
errorEstimatorKindMarginal: Kind of marginal error estimator.
polybasisMarginal: Type of polynomial basis for marginal interpolation.
paramsMarginal: Dictionary of parameters for marginal interpolation.
greedyTolMarginal: Uniform error tolerance for marginal greedy
algorithm.
maxIterMarginal: Maximum number of marginal greedy steps.
radialDirectionalWeightsMarginal: Radial basis weights for marginal
interpolant.
muBounds: list of bounds for pivot parameter values.
muBoundsMarginal: list of bounds for marginal parameter values.
samplingEngine: Sampling engine.
uHF: High fidelity solution(s) with parameter(s) lastSolvedHF as
sampleList.
lastSolvedHF: Parameter(s) corresponding to last computed high fidelity
solution(s) as parameterList.
uApproxReduced: Reduced approximate solution(s) with parameter(s)
lastSolvedApprox as sampleList.
lastSolvedApproxReduced: Parameter(s) corresponding to last computed
reduced approximate solution(s) as parameterList.
uApprox: Approximate solution(s) with parameter(s) lastSolvedApprox as
sampleList.
lastSolvedApprox: Parameter(s) corresponding to last computed
approximate solution(s) as parameterList.
"""
def _poleMatching(self):
vbMng(self, "INIT", "Compressing and matching poles.", 10)
self.trainedModel.initializeFromRational(self.matchingWeight,
self.HFEngine, False)
vbMng(self, "DEL", "Done compressing and matching poles.", 10)
def _updateTrainedModelMarginalSamples(self, idx : ListAny = []):
self.trainedModel.updateEffectiveSamples(idx, self.matchingWeight,
self.HFEngine, False)
def setupApprox(self, *args, **kwargs) -> int:
if self.checkComputedApprox(): return -1
self.purgeparamsMarginal()
_polybasisMarginal = self.polybasisMarginal
self._polybasisMarginal = ("PIECEWISE_LINEAR_"
+ self.samplerMarginal.kind)
setupOK = super().setupApprox(*args, **kwargs)
self._polybasisMarginal = _polybasisMarginal
self._finalizeMarginalization()
return setupOK
diff --git a/rrompy/reduction_methods/pivoted/greedy/rational_interpolant_greedy_pivoted_greedy.py b/rrompy/reduction_methods/pivoted/greedy/rational_interpolant_greedy_pivoted_greedy.py
index 1066d25..8649fdb 100644
--- a/rrompy/reduction_methods/pivoted/greedy/rational_interpolant_greedy_pivoted_greedy.py
+++ b/rrompy/reduction_methods/pivoted/greedy/rational_interpolant_greedy_pivoted_greedy.py
@@ -1,500 +1,502 @@
#Copyright (C) 2018 by the RROMPy authors
#
# This file is part of RROMPy.
#
# RROMPy is free software: you can redistribute it and/or modify
# it under the terms of the GNU Lesser General Public License as published by
# the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
#
# RROMPy is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU Lesser General Public License for more details.
#
# You should have received a copy of the GNU Lesser General Public License
# along with RROMPy. If not, see .
#
from copy import deepcopy as copy
import numpy as np
from .generic_pivoted_greedy_approximant import (
GenericPivotedGreedyApproximantBase,
GenericPivotedGreedyApproximantNoMatch,
GenericPivotedGreedyApproximant)
from rrompy.reduction_methods.standard.greedy import RationalInterpolantGreedy
from rrompy.reduction_methods.standard.greedy.generic_greedy_approximant \
import pruneSamples
from rrompy.reduction_methods.pivoted import (
RationalInterpolantGreedyPivotedNoMatch,
RationalInterpolantGreedyPivoted)
from rrompy.utilities.base.types import Np1D, Tuple, paramVal, paramList
from rrompy.utilities.base import verbosityManager as vbMng
from rrompy.utilities.exception_manager import RROMPyAssert
from rrompy.parameter import emptyParameterList
from rrompy.utilities.parallel import poolRank, recv
__all__ = ['RationalInterpolantGreedyPivotedGreedyNoMatch',
'RationalInterpolantGreedyPivotedGreedy']
class RationalInterpolantGreedyPivotedGreedyBase(
GenericPivotedGreedyApproximantBase):
@property
def sampleBatchSize(self):
"""Value of sampleBatchSize."""
return 1
@property
def sampleBatchIdx(self):
"""Value of sampleBatchIdx."""
return self.S
def greedyNextSample(self, muidx:int, plotEst : str = "NONE")\
-> Tuple[Np1D, int, float, paramVal]:
"""Compute next greedy snapshot of solution map."""
RROMPyAssert(self._mode, message = "Cannot add greedy sample.")
mus = copy(self.muTest[muidx])
self.muTest.pop(muidx)
for j, mu in enumerate(mus):
vbMng(self, "MAIN",
("Adding sample point no. {} at {} to training "
"set.").format(len(self.mus) + 1, mu), 3)
self.mus.append(mu)
self._S = len(self.mus)
self._approxParameters["S"] = self.S
if (self.samplingEngine.nsamples <= len(mus) - j - 1
or not np.allclose(mu, self.samplingEngine.mus[j - len(mus)])):
self.samplingEngine.nextSample(mu)
if self._isLastSampleCollinear():
vbMng(self, "MAIN",
("Collinearity above tolerance detected. Starting "
"preemptive greedy loop termination."), 3)
self._collinearityFlag = 1
errorEstTest = np.empty(len(self.muTest))
errorEstTest[:] = np.nan
return errorEstTest, [-1], np.nan, np.nan
errorEstTest, muidx, maxErrorEst = self.errorEstimator(self.muTest,
True)
if plotEst == "ALL":
self.plotEstimator(errorEstTest, muidx, maxErrorEst)
return errorEstTest, muidx, maxErrorEst, self.muTest[muidx]
def _setSampleBatch(self, maxS:int):
return self.S
def _preliminaryTraining(self):
"""Initialize starting snapshots of solution map."""
RROMPyAssert(self._mode, message = "Cannot start greedy algorithm.")
if self.samplingEngine.nsamples > 0: return
self.resetSamples()
self.samplingEngine.scaleFactor = self.scaleFactorDer
musPivot = self.trainSetGenerator.generatePoints(self.S)
while len(musPivot) > self.S: musPivot.pop()
muTestBasePivot = self.samplerPivot.generatePoints(self.nTestPoints,
False)
idxPop = pruneSamples(self.mapParameterListPivot(muTestBasePivot),
self.mapParameterListPivot(musPivot),
1e-10 * self.scaleFactorPivot[0])
muTestBasePivot.pop(idxPop)
self._mus = emptyParameterList()
self.mus.reset((self.S - 1, self.HFEngine.npar))
self.muTest = emptyParameterList()
self.muTest.reset((len(muTestBasePivot) + 1, self.HFEngine.npar))
for k in range(self.S - 1):
muk = np.empty_like(self.mus[0])
muk[self.directionPivot] = musPivot[k]
muk[self.directionMarginal] = self.muMargLoc
self.mus[k] = muk
for k in range(len(muTestBasePivot)):
muk = np.empty_like(self.muTest[0])
muk[self.directionPivot] = muTestBasePivot[k]
muk[self.directionMarginal] = self.muMargLoc
self.muTest[k] = muk
muk = np.empty_like(self.mus[0])
muk[self.directionPivot] = musPivot[-1]
muk[self.directionMarginal] = self.muMargLoc
self.muTest[-1] = muk
if len(self.mus) > 0:
vbMng(self, "MAIN",
("Adding first {} sample point{} at {} to training "
"set.").format(self.S - 1, "" + "s" * (self.S > 2),
self.mus), 3)
self.samplingEngine.iterSample(self.mus)
self._S = len(self.mus)
self._approxParameters["S"] = self.S
self.M, self.N = ("AUTO",) * 2
def setupApproxPivoted(self, mus:paramList) -> int:
if self.checkComputedApproxPivoted(): return -1
RROMPyAssert(self._mode, message = "Cannot setup approximant.")
vbMng(self, "INIT", "Setting up pivoted approximant.", 10)
if not hasattr(self, "_plotEstPivot"): self._plotEstPivot = "NONE"
idx, sizes, emptyCores = self._preSetupApproxPivoted(mus)
S0 = copy(self.S)
pMat, Ps, Qs, req, musA = None, [], [], [], None
if len(idx) == 0:
vbMng(self, "MAIN", "Idling.", 45)
if self.storeAllSamples: self.storeSamples()
pL, pT, mT = recv(source = 0, tag = poolRank())
pMat = np.empty((pL, 0), dtype = pT)
musA = np.empty((0, self.mu0.shape[1]), dtype = mT)
else:
for i in idx:
self.muMargLoc = mus[i]
vbMng(self, "MAIN", "Building marginal model no. {} at "
"{}.".format(i + 1, self.muMargLoc), 25)
self.samplingEngine.resetHistory()
self.trainedModel = None
self.verbosity -= 5
self.samplingEngine.verbosity -= 5
RationalInterpolantGreedy.setupApprox(self, self._plotEstPivot)
self.verbosity += 5
self.samplingEngine.verbosity += 5
if self.storeAllSamples: self.storeSamples(i + self._nmusOld)
pMat, req, musA = self._localPivotedResult(pMat, req,
emptyCores, musA)
Ps += [copy(self.trainedModel.data.P)]
Qs += [copy(self.trainedModel.data.Q)]
self._S = S0
del self.muMargLoc
for r in req: r.wait()
self._postSetupApproxPivoted(musA, pMat, Ps, Qs, sizes)
vbMng(self, "DEL", "Done setting up pivoted approximant.", 10)
return 0
def setupApprox(self, plotEst : str = "NONE") -> int:
if self.checkComputedApprox(): return -1
if '_' not in plotEst: plotEst = plotEst + "_NONE"
plotEstM, self._plotEstPivot = plotEst.split("_")
val = super().setupApprox(plotEstM)
return val
class RationalInterpolantGreedyPivotedGreedyNoMatch(
RationalInterpolantGreedyPivotedGreedyBase,
GenericPivotedGreedyApproximantNoMatch,
RationalInterpolantGreedyPivotedNoMatch):
"""
ROM greedy pivoted greedy rational interpolant computation for parametric
problems (without pole matching).
Args:
HFEngine: HF problem solver.
mu0(optional): Default parameter. Defaults to 0.
directionPivot(optional): Pivot components. Defaults to [0].
approxParameters(optional): Dictionary containing values for main
parameters of approximant. Recognized keys are:
- 'POD': kind of snapshots orthogonalization; allowed values
include 0, 1/2, and 1; defaults to 1, i.e. POD;
- 'scaleFactorDer': scaling factors for derivative computation;
defaults to 'AUTO';
- 'matchingWeightError': weight for pole matching optimization in
error estimation; defaults to 0;
- 'S': total number of pivot samples current approximant relies
upon;
- 'samplerPivot': pivot sample point generator;
- 'SMarginal': number of starting marginal samples;
- 'samplerMarginal': marginal sample point generator via sparse
grid;
- 'errorEstimatorKindMarginal': kind of marginal error estimator;
available values include 'LOOK_AHEAD' and 'LOOK_AHEAD_RECOVER';
defaults to 'NONE';
- 'polybasis': type of polynomial basis for pivot interpolation;
defaults to 'MONOMIAL';
- 'greedyTol': uniform error tolerance for greedy algorithm;
defaults to 1e-2;
- 'collinearityTol': collinearity tolerance for greedy algorithm;
defaults to 0.;
- 'maxIter': maximum number of greedy steps; defaults to 1e2;
- 'nTestPoints': number of test points; defaults to 5e2;
- 'trainSetGenerator': training sample points generator; defaults
to sampler;
- 'errorEstimatorKind': kind of error estimator; available values
include 'AFFINE', 'DISCREPANCY', 'LOOK_AHEAD',
'LOOK_AHEAD_RES', and 'NONE'; defaults to 'NONE';
- 'greedyTolMarginal': uniform error tolerance for marginal greedy
algorithm; defaults to 1e-1;
- 'maxIterMarginal': maximum number of marginal greedy steps;
defaults to 1e2;
- 'radialDirectionalWeightsMarginal': radial basis weights for
marginal interpolant; defaults to 1;
- 'functionalSolve': strategy for minimization of denominator
functional; allowed values include 'NORM', 'DOMINANT', 'NODAL',
- 'LOEWNER', and 'BARYCENTRIC'; defaults to 'NORM';
+ 'BARYCENTRIC_NORM', and 'BARYCENTRIC[_AVERAGE]' (check pdf in
+ main folder for meaning); defaults to 'NORM';
- 'interpRcond': tolerance for pivot interpolation; defaults to
None;
- 'robustTol': tolerance for robust rational denominator
management; defaults to 0.
Defaults to empty dict.
approx_state(optional): Whether to approximate state. Defaults to
False.
verbosity(optional): Verbosity level. Defaults to 10.
Attributes:
HFEngine: HF problem solver.
mu0: Default parameter.
directionPivot: Pivot components.
mus: Array of snapshot parameters.
musMarginal: Array of marginal snapshot parameters.
approxParameters: Dictionary containing values for main parameters of
approximant. Recognized keys are in parameterList.
parameterListSoft: Recognized keys of soft approximant parameters:
- 'POD': kind of snapshots orthogonalization;
- 'scaleFactorDer': scaling factors for derivative computation;
- 'matchingWeightError': weight for pole matching optimization in
error estimation;
- 'errorEstimatorKindMarginal': kind of marginal error estimator;
- 'polybasis': type of polynomial basis for pivot interpolation;
- 'greedyTol': uniform error tolerance for greedy algorithm;
- 'collinearityTol': collinearity tolerance for greedy algorithm;
- 'maxIter': maximum number of greedy steps;
- 'nTestPoints': number of test points;
- 'trainSetGenerator': training sample points generator;
- 'errorEstimatorKind': kind of error estimator;
- 'greedyTolMarginal': uniform error tolerance for marginal greedy
algorithm;
- 'maxIterMarginal': maximum number of marginal greedy steps;
- 'radialDirectionalWeightsMarginal': radial basis weights for
marginal interpolant;
- 'functionalSolve': strategy for minimization of denominator
functional;
- 'interpRcond': tolerance for pivot interpolation;
- 'robustTol': tolerance for robust rational denominator
management.
parameterListCritical: Recognized keys of critical approximant
parameters:
- 'S': total number of pivot samples current approximant relies
upon;
- 'samplerPivot': pivot sample point generator;
- 'SMarginal': total number of marginal samples current approximant
relies upon;
- 'samplerMarginal': marginal sample point generator via sparse
grid.
approx_state: Whether to approximate state.
verbosity: Verbosity level.
POD: Kind of snapshots orthogonalization.
scaleFactorDer: Scaling factors for derivative computation.
matchingWeightError: Weight for pole matching optimization in error
estimation.
S: Total number of pivot samples current approximant relies upon.
samplerPivot: Pivot sample point generator.
SMarginal: Total number of marginal samples current approximant relies
upon.
samplerMarginal: Marginal sample point generator via sparse grid.
errorEstimatorKindMarginal: Kind of marginal error estimator.
polybasis: Type of polynomial basis for pivot interpolation.
greedyTol: uniform error tolerance for greedy algorithm.
collinearityTol: Collinearity tolerance for greedy algorithm.
maxIter: maximum number of greedy steps.
nTestPoints: number of starting training points.
trainSetGenerator: training sample points generator.
errorEstimatorKind: kind of error estimator.
greedyTolMarginal: Uniform error tolerance for marginal greedy
algorithm.
maxIterMarginal: Maximum number of marginal greedy steps.
radialDirectionalWeightsMarginal: Radial basis weights for marginal
interpolant.
functionalSolve: Strategy for minimization of denominator functional.
interpRcond: Tolerance for pivot interpolation.
robustTol: Tolerance for robust rational denominator management.
muBounds: list of bounds for pivot parameter values.
muBoundsMarginal: list of bounds for marginal parameter values.
samplingEngine: Sampling engine.
uHF: High fidelity solution(s) with parameter(s) lastSolvedHF as
sampleList.
lastSolvedHF: Parameter(s) corresponding to last computed high fidelity
solution(s) as parameterList.
uApproxReduced: Reduced approximate solution(s) with parameter(s)
lastSolvedApprox as sampleList.
lastSolvedApproxReduced: Parameter(s) corresponding to last computed
reduced approximate solution(s) as parameterList.
uApprox: Approximate solution(s) with parameter(s) lastSolvedApprox as
sampleList.
lastSolvedApprox: Parameter(s) corresponding to last computed
approximate solution(s) as parameterList.
"""
class RationalInterpolantGreedyPivotedGreedy(
RationalInterpolantGreedyPivotedGreedyBase,
GenericPivotedGreedyApproximant,
RationalInterpolantGreedyPivoted):
"""
ROM greedy pivoted greedy rational interpolant computation for parametric
problems (with pole matching).
Args:
HFEngine: HF problem solver.
mu0(optional): Default parameter. Defaults to 0.
directionPivot(optional): Pivot components. Defaults to [0].
approxParameters(optional): Dictionary containing values for main
parameters of approximant. Recognized keys are:
- 'POD': kind of snapshots orthogonalization; allowed values
include 0, 1/2, and 1; defaults to 1, i.e. POD;
- 'scaleFactorDer': scaling factors for derivative computation;
defaults to 'AUTO';
- 'matchingWeight': weight for pole matching optimization; defaults
to 1;
- 'sharedRatio': required ratio of marginal points to share
resonance; defaults to 1.;
- 'matchingWeightError': weight for pole matching optimization in
error estimation; defaults to 0;
- 'S': total number of pivot samples current approximant relies
upon;
- 'samplerPivot': pivot sample point generator;
- 'SMarginal': number of starting marginal samples;
- 'samplerMarginal': marginal sample point generator via sparse
grid;
- 'errorEstimatorKindMarginal': kind of marginal error estimator;
available values include 'LOOK_AHEAD' and 'LOOK_AHEAD_RECOVER';
defaults to 'NONE';
- 'polybasis': type of polynomial basis for pivot interpolation;
defaults to 'MONOMIAL';
- 'polybasisMarginal': type of polynomial basis for marginal
interpolation; allowed values include 'MONOMIAL_*',
'CHEBYSHEV_*', 'LEGENDRE_*', 'NEARESTNEIGHBOR', and
'PIECEWISE_LINEAR_*'; defaults to 'MONOMIAL';
- 'paramsMarginal': dictionary of parameters for marginal
interpolation; include:
. 'MMarginal': degree of marginal interpolant; defaults to
'AUTO', i.e. maximum allowed; not for 'NEARESTNEIGHBOR' or
'PIECEWISE_LINEAR_*';
. 'nNeighborsMarginal': number of marginal nearest neighbors;
defaults to 1; only for 'NEARESTNEIGHBOR';
. 'polydegreetypeMarginal': type of polynomial degree for
marginal; defaults to 'TOTAL'; not for 'NEARESTNEIGHBOR' or
'PIECEWISE_LINEAR_*';
. 'interpRcondMarginal': tolerance for marginal interpolation;
defaults to None; not for 'NEARESTNEIGHBOR' or
'PIECEWISE_LINEAR_*';
. 'radialDirectionalWeightsMarginalAdapt': bounds for adaptive
rescaling of marginal radial basis weights; only for
radial basis.
- 'greedyTol': uniform error tolerance for greedy algorithm;
defaults to 1e-2;
- 'collinearityTol': collinearity tolerance for greedy algorithm;
defaults to 0.;
- 'maxIter': maximum number of greedy steps; defaults to 1e2;
- 'nTestPoints': number of test points; defaults to 5e2;
- 'trainSetGenerator': training sample points generator; defaults
to sampler;
- 'errorEstimatorKind': kind of error estimator; available values
include 'AFFINE', 'DISCREPANCY', 'LOOK_AHEAD',
'LOOK_AHEAD_RES', and 'NONE'; defaults to 'NONE';
- 'greedyTolMarginal': uniform error tolerance for marginal greedy
algorithm; defaults to 1e-1;
- 'maxIterMarginal': maximum number of marginal greedy steps;
defaults to 1e2;
- 'radialDirectionalWeightsMarginal': radial basis weights for
marginal interpolant; defaults to 1;
- 'functionalSolve': strategy for minimization of denominator
functional; allowed values include 'NORM', 'DOMINANT', 'NODAL',
- 'LOEWNER', and 'BARYCENTRIC'; defaults to 'NORM';
+ 'BARYCENTRIC_NORM', and 'BARYCENTRIC[_AVERAGE]' (check pdf in
+ main folder for meaning); defaults to 'NORM';
- 'interpRcond': tolerance for pivot interpolation; defaults to
None;
- 'robustTol': tolerance for robust rational denominator
management; defaults to 0.
Defaults to empty dict.
approx_state(optional): Whether to approximate state. Defaults to
False.
verbosity(optional): Verbosity level. Defaults to 10.
Attributes:
HFEngine: HF problem solver.
mu0: Default parameter.
directionPivot: Pivot components.
mus: Array of snapshot parameters.
musMarginal: Array of marginal snapshot parameters.
approxParameters: Dictionary containing values for main parameters of
approximant. Recognized keys are in parameterList.
parameterListSoft: Recognized keys of soft approximant parameters:
- 'POD': kind of snapshots orthogonalization;
- 'scaleFactorDer': scaling factors for derivative computation;
- 'matchingWeight': weight for pole matching optimization;
- 'sharedRatio': required ratio of marginal points to share
resonance;
- 'matchingWeightError': weight for pole matching optimization in
error estimation;
- 'errorEstimatorKindMarginal': kind of marginal error estimator;
- 'polybasis': type of polynomial basis for pivot interpolation;
- 'polybasisMarginal': type of polynomial basis for marginal
interpolation;
- 'paramsMarginal': dictionary of parameters for marginal
interpolation; include:
. 'MMarginal': degree of marginal interpolant;
. 'nNeighborsMarginal': number of marginal nearest neighbors;
. 'polydegreetypeMarginal': type of polynomial degree for
marginal;
. 'interpRcondMarginal': tolerance for marginal interpolation;
. 'radialDirectionalWeightsMarginalAdapt': bounds for adaptive
rescaling of marginal radial basis weights.
- 'greedyTol': uniform error tolerance for greedy algorithm;
- 'collinearityTol': collinearity tolerance for greedy algorithm;
- 'maxIter': maximum number of greedy steps;
- 'nTestPoints': number of test points;
- 'trainSetGenerator': training sample points generator;
- 'errorEstimatorKind': kind of error estimator;
- 'greedyTolMarginal': uniform error tolerance for marginal greedy
algorithm;
- 'maxIterMarginal': maximum number of marginal greedy steps;
- 'radialDirectionalWeightsMarginal': radial basis weights for
marginal interpolant;
- 'functionalSolve': strategy for minimization of denominator
functional;
- 'interpRcond': tolerance for pivot interpolation;
- 'robustTol': tolerance for robust rational denominator
management.
parameterListCritical: Recognized keys of critical approximant
parameters:
- 'S': total number of pivot samples current approximant relies
upon;
- 'samplerPivot': pivot sample point generator;
- 'SMarginal': total number of marginal samples current approximant
relies upon;
- 'samplerMarginal': marginal sample point generator via sparse
grid.
approx_state: Whether to approximate state.
verbosity: Verbosity level.
POD: Kind of snapshots orthogonalization.
scaleFactorDer: Scaling factors for derivative computation.
matchingWeight: Weight for pole matching optimization.
sharedRatio: Required ratio of marginal points to share resonance.
matchingWeightError: Weight for pole matching optimization in error
estimation.
S: Total number of pivot samples current approximant relies upon.
samplerPivot: Pivot sample point generator.
SMarginal: Total number of marginal samples current approximant relies
upon.
samplerMarginal: Marginal sample point generator via sparse grid.
errorEstimatorKindMarginal: Kind of marginal error estimator.
polybasis: Type of polynomial basis for pivot interpolation.
polybasisMarginal: Type of polynomial basis for marginal interpolation.
paramsMarginal: Dictionary of parameters for marginal interpolation.
greedyTol: uniform error tolerance for greedy algorithm.
collinearityTol: Collinearity tolerance for greedy algorithm.
maxIter: maximum number of greedy steps.
nTestPoints: number of starting training points.
trainSetGenerator: training sample points generator.
errorEstimatorKind: kind of error estimator.
greedyTolMarginal: Uniform error tolerance for marginal greedy
algorithm.
maxIterMarginal: Maximum number of marginal greedy steps.
radialDirectionalWeightsMarginal: Radial basis weights for marginal
interpolant.
functionalSolve: Strategy for minimization of denominator functional.
interpRcond: Tolerance for pivot interpolation.
robustTol: Tolerance for robust rational denominator management.
muBounds: list of bounds for pivot parameter values.
muBoundsMarginal: list of bounds for marginal parameter values.
samplingEngine: Sampling engine.
uHF: High fidelity solution(s) with parameter(s) lastSolvedHF as
sampleList.
lastSolvedHF: Parameter(s) corresponding to last computed high fidelity
solution(s) as parameterList.
uApproxReduced: Reduced approximate solution(s) with parameter(s)
lastSolvedApprox as sampleList.
lastSolvedApproxReduced: Parameter(s) corresponding to last computed
reduced approximate solution(s) as parameterList.
uApprox: Approximate solution(s) with parameter(s) lastSolvedApprox as
sampleList.
lastSolvedApprox: Parameter(s) corresponding to last computed
approximate solution(s) as parameterList.
"""
diff --git a/rrompy/reduction_methods/pivoted/greedy/rational_interpolant_pivoted_greedy.py b/rrompy/reduction_methods/pivoted/greedy/rational_interpolant_pivoted_greedy.py
index 585da53..724a683 100644
--- a/rrompy/reduction_methods/pivoted/greedy/rational_interpolant_pivoted_greedy.py
+++ b/rrompy/reduction_methods/pivoted/greedy/rational_interpolant_pivoted_greedy.py
@@ -1,426 +1,428 @@
# Copyright (C) 2018 by the RROMPy authors
#
# This file is part of RROMPy.
#
# RROMPy is free software: you can redistribute it and/or modify
# it under the terms of the GNU Lesser General Public License as published by
# the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
#
# RROMPy is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU Lesser General Public License for more details.
#
# You should have received a copy of the GNU Lesser General Public License
# along with RROMPy. If not, see .
#
from copy import deepcopy as copy
from numpy import empty, empty_like
from .generic_pivoted_greedy_approximant import (
GenericPivotedGreedyApproximantBase,
GenericPivotedGreedyApproximantNoMatch,
GenericPivotedGreedyApproximant)
from rrompy.reduction_methods.standard import RationalInterpolant
from rrompy.reduction_methods.pivoted import (
RationalInterpolantPivotedNoMatch,
RationalInterpolantPivoted)
from rrompy.utilities.base.types import paramList
from rrompy.utilities.base import verbosityManager as vbMng
from rrompy.utilities.exception_manager import RROMPyAssert
from rrompy.parameter import emptyParameterList
from rrompy.utilities.parallel import poolRank, recv
__all__ = ['RationalInterpolantPivotedGreedyNoMatch',
'RationalInterpolantPivotedGreedy']
class RationalInterpolantPivotedGreedyBase(
GenericPivotedGreedyApproximantBase):
def computeSnapshots(self):
"""Compute snapshots of solution map."""
RROMPyAssert(self._mode,
message = "Cannot start snapshot computation.")
vbMng(self, "INIT", "Starting computation of snapshots.", 5)
self.samplingEngine.scaleFactor = self.scaleFactorDer
if not hasattr(self, "musPivot") or len(self.musPivot) != self.S:
self.musPivot = self.samplerPivot.generatePoints(self.S)
while len(self.musPivot) > self.S: self.musPivot.pop()
musLoc = emptyParameterList()
musLoc.reset((self.S, self.HFEngine.npar))
self.samplingEngine.resetHistory()
for k in range(self.S):
muk = empty_like(musLoc[0])
muk[self.directionPivot] = self.musPivot[k]
muk[self.directionMarginal] = self.muMargLoc
musLoc[k] = muk
self.samplingEngine.iterSample(musLoc)
vbMng(self, "DEL", "Done computing snapshots.", 5)
self._m_selfmus = copy(musLoc)
self._mus = self.musPivot
self._m_HFEparameterMap = copy(self.HFEngine.parameterMap)
self.HFEngine.parameterMap = {
"F": [self.HFEngine.parameterMap["F"][self.directionPivot[0]]],
"B": [self.HFEngine.parameterMap["B"][self.directionPivot[0]]]}
def setupApproxPivoted(self, mus:paramList) -> int:
if self.checkComputedApproxPivoted(): return -1
RROMPyAssert(self._mode, message = "Cannot setup approximant.")
vbMng(self, "INIT", "Setting up pivoted approximant.", 10)
idx, sizes, emptyCores = self._preSetupApproxPivoted(mus)
pMat, Ps, Qs, req, musA = None, [], [], [], None
if len(idx) == 0:
vbMng(self, "MAIN", "Idling.", 45)
if self.storeAllSamples: self.storeSamples()
pL, pT, mT = recv(source = 0, tag = poolRank())
pMat = empty((pL, 0), dtype = pT)
musA = empty((0, self.mu0.shape[1]), dtype = mT)
else:
for i in idx:
self.muMargLoc = mus[i]
vbMng(self, "MAIN", "Building marginal model no. {} at "
"{}.".format(i + 1, self.muMargLoc), 25)
self.samplingEngine.resetHistory()
self.trainedModel = None
self.verbosity -= 5
self.samplingEngine.verbosity -= 5
RationalInterpolant.setupApprox(self)
self.verbosity += 5
self.samplingEngine.verbosity += 5
self._mus = self._m_selfmus
self.HFEngine.parameterMap = self._m_HFEparameterMap
del self._m_selfmus, self._m_HFEparameterMap
if self.storeAllSamples: self.storeSamples(i + self._nmusOld)
pMat, req, musA = self._localPivotedResult(pMat, req,
emptyCores, musA)
Ps += [copy(self.trainedModel.data.P)]
Qs += [copy(self.trainedModel.data.Q)]
del self.muMargLoc
for r in req: r.wait()
self._postSetupApproxPivoted(musA, pMat, Ps, Qs, sizes)
vbMng(self, "DEL", "Done setting up pivoted approximant.", 10)
return 0
class RationalInterpolantPivotedGreedyNoMatch(
RationalInterpolantPivotedGreedyBase,
GenericPivotedGreedyApproximantNoMatch,
RationalInterpolantPivotedNoMatch):
"""
ROM pivoted greedy rational interpolant computation for parametric
problems (without pole matching).
Args:
HFEngine: HF problem solver.
mu0(optional): Default parameter. Defaults to 0.
directionPivot(optional): Pivot components. Defaults to [0].
approxParameters(optional): Dictionary containing values for main
parameters of approximant. Recognized keys are:
- 'POD': kind of snapshots orthogonalization; allowed values
include 0, 1/2, and 1; defaults to 1, i.e. POD;
- 'scaleFactorDer': scaling factors for derivative computation;
defaults to 'AUTO';
- 'matchingWeightError': weight for pole matching optimization in
error estimation; defaults to 0;
- 'S': total number of pivot samples current approximant relies
upon;
- 'samplerPivot': pivot sample point generator;
- 'SMarginal': number of starting marginal samples;
- 'samplerMarginal': marginal sample point generator via sparse
grid;
- 'errorEstimatorKindMarginal': kind of marginal error estimator;
available values include 'LOOK_AHEAD' and 'LOOK_AHEAD_RECOVER';
defaults to 'NONE';
- 'polybasis': type of polynomial basis for pivot interpolation;
defaults to 'MONOMIAL';
- 'M': degree of rational interpolant numerator; defaults to
'AUTO', i.e. maximum allowed;
- 'N': degree of rational interpolant denominator; defaults to
'AUTO', i.e. maximum allowed;
- 'greedyTolMarginal': uniform error tolerance for marginal greedy
algorithm; defaults to 1e-1;
- 'maxIterMarginal': maximum number of marginal greedy steps;
defaults to 1e2;
- 'radialDirectionalWeights': radial basis weights for pivot
numerator; defaults to 1;
- 'radialDirectionalWeightsAdapt': bounds for adaptive rescaling of
radial basis weights; defaults to [-1, -1];
- 'radialDirectionalWeightsMarginal': radial basis weights for
marginal interpolant; defaults to 1;
- 'functionalSolve': strategy for minimization of denominator
functional; allowed values include 'NORM', 'DOMINANT', 'NODAL',
- 'LOEWNER', and 'BARYCENTRIC'; defaults to 'NORM';
+ 'BARYCENTRIC_NORM', and 'BARYCENTRIC[_AVERAGE]' (check pdf in
+ main folder for meaning); defaults to 'NORM';
- 'interpRcond': tolerance for pivot interpolation; defaults to
None;
- 'robustTol': tolerance for robust rational denominator
management; defaults to 0.
Defaults to empty dict.
approx_state(optional): Whether to approximate state. Defaults to
False.
verbosity(optional): Verbosity level. Defaults to 10.
Attributes:
HFEngine: HF problem solver.
mu0: Default parameter.
directionPivot: Pivot components.
mus: Array of snapshot parameters.
musPivot: Array of pivot snapshot parameters.
musMarginal: Array of marginal snapshot parameters.
approxParameters: Dictionary containing values for main parameters of
approximant. Recognized keys are in parameterList.
parameterListSoft: Recognized keys of soft approximant parameters:
- 'POD': kind of snapshots orthogonalization;
- 'scaleFactorDer': scaling factors for derivative computation;
- 'matchingWeightError': weight for pole matching optimization in
error estimation;
- 'errorEstimatorKindMarginal': kind of marginal error estimator;
- 'polybasis': type of polynomial basis for pivot interpolation;
- 'M': degree of rational interpolant numerator;
- 'N': degree of rational interpolant denominator;
- 'greedyTolMarginal': uniform error tolerance for marginal greedy
algorithm;
- 'maxIterMarginal': maximum number of marginal greedy steps;
- 'radialDirectionalWeights': radial basis weights for pivot
numerator;
- 'radialDirectionalWeightsAdapt': bounds for adaptive rescaling of
radial basis weights;
- 'radialDirectionalWeightsMarginal': radial basis weights for
marginal interpolant;
- 'functionalSolve': strategy for minimization of denominator
functional;
- 'interpRcond': tolerance for pivot interpolation;
- 'robustTol': tolerance for robust rational denominator
management.
parameterListCritical: Recognized keys of critical approximant
parameters:
- 'S': total number of pivot samples current approximant relies
upon;
- 'samplerPivot': pivot sample point generator;
- 'SMarginal': total number of marginal samples current approximant
relies upon;
- 'samplerMarginal': marginal sample point generator via sparse
grid.
approx_state: Whether to approximate state.
verbosity: Verbosity level.
POD: Kind of snapshots orthogonalization.
scaleFactorDer: Scaling factors for derivative computation.
matchingWeightError: Weight for pole matching optimization in error
estimation.
S: Total number of pivot samples current approximant relies upon.
samplerPivot: Pivot sample point generator.
SMarginal: Total number of marginal samples current approximant relies
upon.
samplerMarginal: Marginal sample point generator via sparse grid.
errorEstimatorKindMarginal: Kind of marginal error estimator.
polybasis: Type of polynomial basis for pivot interpolation.
M: Degree of rational interpolant numerator.
N: Degree of rational interpolant denominator.
greedyTolMarginal: Uniform error tolerance for marginal greedy
algorithm.
maxIterMarginal: Maximum number of marginal greedy steps.
radialDirectionalWeights: Radial basis weights for pivot numerator.
radialDirectionalWeightsAdapt: Bounds for adaptive rescaling of radial
basis weights.
radialDirectionalWeightsMarginal: Radial basis weights for marginal
interpolant.
functionalSolve: Strategy for minimization of denominator functional.
interpRcond: Tolerance for pivot interpolation.
robustTol: Tolerance for robust rational denominator management.
muBounds: list of bounds for pivot parameter values.
muBoundsMarginal: list of bounds for marginal parameter values.
samplingEngine: Sampling engine.
uHF: High fidelity solution(s) with parameter(s) lastSolvedHF as
sampleList.
lastSolvedHF: Parameter(s) corresponding to last computed high fidelity
solution(s) as parameterList.
uApproxReduced: Reduced approximate solution(s) with parameter(s)
lastSolvedApprox as sampleList.
lastSolvedApproxReduced: Parameter(s) corresponding to last computed
reduced approximate solution(s) as parameterList.
uApprox: Approximate solution(s) with parameter(s) lastSolvedApprox as
sampleList.
lastSolvedApprox: Parameter(s) corresponding to last computed
approximate solution(s) as parameterList.
"""
class RationalInterpolantPivotedGreedy(RationalInterpolantPivotedGreedyBase,
GenericPivotedGreedyApproximant,
RationalInterpolantPivoted):
"""
ROM pivoted greedy rational interpolant computation for parametric
problems (with pole matching).
Args:
HFEngine: HF problem solver.
mu0(optional): Default parameter. Defaults to 0.
directionPivot(optional): Pivot components. Defaults to [0].
approxParameters(optional): Dictionary containing values for main
parameters of approximant. Recognized keys are:
- 'POD': kind of snapshots orthogonalization; allowed values
include 0, 1/2, and 1; defaults to 1, i.e. POD;
- 'scaleFactorDer': scaling factors for derivative computation;
defaults to 'AUTO';
- 'matchingWeight': weight for pole matching optimization; defaults
to 1;
- 'sharedRatio': required ratio of marginal points to share
resonance; defaults to 1.;
- 'matchingWeightError': weight for pole matching optimization in
error estimation; defaults to 0;
- 'S': total number of pivot samples current approximant relies
upon;
- 'samplerPivot': pivot sample point generator;
- 'SMarginal': number of starting marginal samples;
- 'samplerMarginal': marginal sample point generator via sparse
grid;
- 'errorEstimatorKindMarginal': kind of marginal error estimator;
available values include 'LOOK_AHEAD' and 'LOOK_AHEAD_RECOVER';
defaults to 'NONE';
- 'polybasis': type of polynomial basis for pivot interpolation;
defaults to 'MONOMIAL';
- 'polybasisMarginal': type of polynomial basis for marginal
interpolation; allowed values include 'MONOMIAL_*',
'CHEBYSHEV_*', 'LEGENDRE_*', 'NEARESTNEIGHBOR', and
'PIECEWISE_LINEAR_*'; defaults to 'MONOMIAL';
- 'paramsMarginal': dictionary of parameters for marginal
interpolation; include:
. 'MMarginal': degree of marginal interpolant; defaults to
'AUTO', i.e. maximum allowed; not for 'NEARESTNEIGHBOR' or
'PIECEWISE_LINEAR_*';
. 'nNeighborsMarginal': number of marginal nearest neighbors;
defaults to 1; only for 'NEARESTNEIGHBOR';
. 'polydegreetypeMarginal': type of polynomial degree for
marginal; defaults to 'TOTAL'; not for 'NEARESTNEIGHBOR' or
'PIECEWISE_LINEAR_*';
. 'interpRcondMarginal': tolerance for marginal interpolation;
defaults to None; not for 'NEARESTNEIGHBOR' or
'PIECEWISE_LINEAR_*';
. 'radialDirectionalWeightsMarginalAdapt': bounds for adaptive
rescaling of marginal radial basis weights; only for
radial basis.
- 'M': degree of rational interpolant numerator; defaults to
'AUTO', i.e. maximum allowed;
- 'N': degree of rational interpolant denominator; defaults to
'AUTO', i.e. maximum allowed;
- 'greedyTolMarginal': uniform error tolerance for marginal greedy
algorithm; defaults to 1e-1;
- 'maxIterMarginal': maximum number of marginal greedy steps;
defaults to 1e2;
- 'radialDirectionalWeights': radial basis weights for pivot
numerator; defaults to 1;
- 'radialDirectionalWeightsAdapt': bounds for adaptive rescaling of
radial basis weights; defaults to [-1, -1];
- 'radialDirectionalWeightsMarginal': radial basis weights for
marginal interpolant; defaults to 1;
- 'functionalSolve': strategy for minimization of denominator
functional; allowed values include 'NORM', 'DOMINANT', 'NODAL',
- 'LOEWNER', and 'BARYCENTRIC'; defaults to 'NORM';
+ 'BARYCENTRIC_NORM', and 'BARYCENTRIC[_AVERAGE]' (check pdf in
+ main folder for meaning); defaults to 'NORM';
- 'interpRcond': tolerance for pivot interpolation; defaults to
None;
- 'robustTol': tolerance for robust rational denominator
management; defaults to 0.
Defaults to empty dict.
approx_state(optional): Whether to approximate state. Defaults to
False.
verbosity(optional): Verbosity level. Defaults to 10.
Attributes:
HFEngine: HF problem solver.
mu0: Default parameter.
directionPivot: Pivot components.
mus: Array of snapshot parameters.
musPivot: Array of pivot snapshot parameters.
musMarginal: Array of marginal snapshot parameters.
approxParameters: Dictionary containing values for main parameters of
approximant. Recognized keys are in parameterList.
parameterListSoft: Recognized keys of soft approximant parameters:
- 'POD': kind of snapshots orthogonalization;
- 'scaleFactorDer': scaling factors for derivative computation;
- 'matchingWeight': weight for pole matching optimization;
- 'sharedRatio': required ratio of marginal points to share
resonance;
- 'matchingWeightError': weight for pole matching optimization in
error estimation;
- 'errorEstimatorKindMarginal': kind of marginal error estimator;
- 'polybasis': type of polynomial basis for pivot interpolation;
- 'polybasisMarginal': type of polynomial basis for marginal
interpolation;
- 'paramsMarginal': dictionary of parameters for marginal
interpolation; include:
. 'MMarginal': degree of marginal interpolant;
. 'nNeighborsMarginal': number of marginal nearest neighbors;
. 'polydegreetypeMarginal': type of polynomial degree for
marginal;
. 'interpRcondMarginal': tolerance for marginal interpolation;
. 'radialDirectionalWeightsMarginalAdapt': bounds for adaptive
rescaling of marginal radial basis weights.
- 'M': degree of rational interpolant numerator;
- 'N': degree of rational interpolant denominator;
- 'greedyTolMarginal': uniform error tolerance for marginal greedy
algorithm;
- 'maxIterMarginal': maximum number of marginal greedy steps;
- 'radialDirectionalWeights': radial basis weights for pivot
numerator;
- 'radialDirectionalWeightsAdapt': bounds for adaptive rescaling of
radial basis weights;
- 'radialDirectionalWeightsMarginal': radial basis weights for
marginal interpolant;
- 'functionalSolve': strategy for minimization of denominator
functional;
- 'interpRcond': tolerance for pivot interpolation;
- 'robustTol': tolerance for robust rational denominator
management.
parameterListCritical: Recognized keys of critical approximant
parameters:
- 'S': total number of pivot samples current approximant relies
upon;
- 'samplerPivot': pivot sample point generator;
- 'SMarginal': total number of marginal samples current approximant
relies upon;
- 'samplerMarginal': marginal sample point generator via sparse
grid.
approx_state: Whether to approximate state.
verbosity: Verbosity level.
POD: Kind of snapshots orthogonalization.
scaleFactorDer: Scaling factors for derivative computation.
matchingWeight: Weight for pole matching optimization.
sharedRatio: Required ratio of marginal points to share resonance.
matchingWeightError: Weight for pole matching optimization in error
estimation.
S: Total number of pivot samples current approximant relies upon.
samplerPivot: Pivot sample point generator.
SMarginal: Total number of marginal samples current approximant relies
upon.
samplerMarginal: Marginal sample point generator via sparse grid.
errorEstimatorKindMarginal: Kind of marginal error estimator.
polybasis: Type of polynomial basis for pivot interpolation.
polybasisMarginal: Type of polynomial basis for marginal interpolation.
paramsMarginal: Dictionary of parameters for marginal interpolation.
M: Degree of rational interpolant numerator.
N: Degree of rational interpolant denominator.
greedyTolMarginal: Uniform error tolerance for marginal greedy
algorithm.
maxIterMarginal: Maximum number of marginal greedy steps.
radialDirectionalWeights: Radial basis weights for pivot numerator.
radialDirectionalWeightsAdapt: Bounds for adaptive rescaling of radial
basis weights.
radialDirectionalWeightsMarginal: Radial basis weights for marginal
interpolant.
functionalSolve: Strategy for minimization of denominator functional.
interpRcond: Tolerance for pivot interpolation.
robustTol: Tolerance for robust rational denominator management.
muBounds: list of bounds for pivot parameter values.
muBoundsMarginal: list of bounds for marginal parameter values.
samplingEngine: Sampling engine.
uHF: High fidelity solution(s) with parameter(s) lastSolvedHF as
sampleList.
lastSolvedHF: Parameter(s) corresponding to last computed high fidelity
solution(s) as parameterList.
uApproxReduced: Reduced approximate solution(s) with parameter(s)
lastSolvedApprox as sampleList.
lastSolvedApproxReduced: Parameter(s) corresponding to last computed
reduced approximate solution(s) as parameterList.
uApprox: Approximate solution(s) with parameter(s) lastSolvedApprox as
sampleList.
lastSolvedApprox: Parameter(s) corresponding to last computed
approximate solution(s) as parameterList.
"""
diff --git a/rrompy/reduction_methods/pivoted/rational_interpolant_greedy_pivoted.py b/rrompy/reduction_methods/pivoted/rational_interpolant_greedy_pivoted.py
index 671362f..d196135 100644
--- a/rrompy/reduction_methods/pivoted/rational_interpolant_greedy_pivoted.py
+++ b/rrompy/reduction_methods/pivoted/rational_interpolant_greedy_pivoted.py
@@ -1,525 +1,527 @@
# Copyright (C) 2018 by the RROMPy authors
#
# This file is part of RROMPy.
#
# RROMPy is free software: you can redistribute it and/or modify
# it under the terms of the GNU Lesser General Public License as published by
# the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
#
# RROMPy is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU Lesser General Public License for more details.
#
# You should have received a copy of the GNU Lesser General Public License
# along with RROMPy. If not, see .
#
from copy import deepcopy as copy
import numpy as np
from .generic_pivoted_approximant import (GenericPivotedApproximantBase,
GenericPivotedApproximantNoMatch,
GenericPivotedApproximant)
from .gather_pivoted_approximant import gatherPivotedApproximant
from rrompy.reduction_methods.standard.greedy.rational_interpolant_greedy \
import RationalInterpolantGreedy
from rrompy.reduction_methods.standard.greedy.generic_greedy_approximant \
import pruneSamples
from rrompy.utilities.base.types import Np1D
from rrompy.utilities.base import verbosityManager as vbMng
from rrompy.utilities.poly_fitting.polynomial import polyvander as pv
from rrompy.utilities.exception_manager import RROMPyAssert
from rrompy.parameter import emptyParameterList, parameterList
from rrompy.utilities.parallel import poolRank, indicesScatter, isend, recv
__all__ = ['RationalInterpolantGreedyPivotedNoMatch',
'RationalInterpolantGreedyPivoted']
class RationalInterpolantGreedyPivotedBase(GenericPivotedApproximantBase,
RationalInterpolantGreedy):
def __init__(self, *args, **kwargs):
self._preInit()
super().__init__(*args, **kwargs)
self._ignoreResidues = self.nparPivot > 1
self._postInit()
@property
def tModelType(self):
if hasattr(self, "_temporaryPivot"):
return RationalInterpolantGreedy.tModelType.fget(self)
return super().tModelType
def _polyvanderAuxiliary(self, mus, deg, *args):
degEff = [0] * self.npar
degEff[self.directionPivot[0]] = deg
return pv(mus, degEff, *args)
def _marginalizeMiscellanea(self, forward:bool):
if forward:
self._m_selfmus = copy(self.mus)
self._m_HFEparameterMap = copy(self.HFEngine.parameterMap)
self._mus = self.checkParameterListPivot(
self.mus(self.directionPivot))
self.HFEngine.parameterMap = {
"F": [self.HFEngine.parameterMap["F"][self.directionPivot[0]]],
"B": [self.HFEngine.parameterMap["B"][self.directionPivot[0]]]}
else:
self._mus = self._m_selfmus
self.HFEngine.parameterMap = self._m_HFEparameterMap
del self._m_selfmus, self._m_HFEparameterMap
def _marginalizeTrainedModel(self, forward:bool):
if forward:
del self._temporaryPivot
self.trainedModel.data.mu0 = self.mu0
self.trainedModel.data.scaleFactor = [1.] * self.npar
self.trainedModel.data.scaleFactor[self.directionPivot[0]] = (
self.scaleFactor[0])
self.trainedModel.data.parameterMap = self.HFEngine.parameterMap
self._m_musUniqueCN = copy(self._musUniqueCN)
musUniqueCNAux = np.zeros((self.S, self.npar),
dtype = self._musUniqueCN.dtype)
musUniqueCNAux[:, self.directionPivot[0]] = self._musUniqueCN(0)
self._musUniqueCN = self.checkParameterList(musUniqueCNAux)
self._m_derIdxs = copy(self._derIdxs)
for j in range(len(self._derIdxs)):
for l in range(len(self._derIdxs[j])):
derjl = self._derIdxs[j][l][0]
self._derIdxs[j][l] = [0] * self.npar
self._derIdxs[j][l][self.directionPivot[0]] = derjl
self.trainedModel.data.Q._dirPivot = self.directionPivot[0]
self.trainedModel.data.P._dirPivot = self.directionPivot[0]
else:
self._temporaryPivot = 1
self.trainedModel.data.mu0 = self.checkParameterListPivot(
self.mu0(self.directionPivot))
self.trainedModel.data.scaleFactor = self.scaleFactor
self.trainedModel.data.parameterMap = {
"F": [self.HFEngine.parameterMap["F"][self.directionPivot[0]]],
"B": [self.HFEngine.parameterMap["B"][self.directionPivot[0]]]}
self._musUniqueCN = copy(self._m_musUniqueCN)
self._derIdxs = copy(self._m_derIdxs)
del self._m_musUniqueCN, self._m_derIdxs
del self.trainedModel.data.Q._dirPivot
del self.trainedModel.data.P._dirPivot
self.trainedModel.data.npar = self.npar
def errorEstimator(self, mus:Np1D, return_max : bool = False) -> Np1D:
"""Standard residual-based error estimator."""
setupOK = self.setupApproxLocal()
if setupOK > 0:
err = np.empty(len(mus))
err[:] = np.nan
if not return_max: return err
return err, [- setupOK], np.nan
self._marginalizeTrainedModel(True)
errRes = super().errorEstimator(mus, return_max)
self._marginalizeTrainedModel(False)
return errRes
def _preliminaryTraining(self):
"""Initialize starting snapshots of solution map."""
RROMPyAssert(self._mode, message = "Cannot start greedy algorithm.")
self._S = self._setSampleBatch(self.S)
self.resetSamples()
self.samplingEngine.scaleFactor = self.scaleFactorDer
musPivot = self.trainSetGenerator.generatePoints(self.S)
while len(musPivot) > self.S: musPivot.pop()
muTestPivot = self.samplerPivot.generatePoints(self.nTestPoints, False)
idxPop = pruneSamples(self.mapParameterListPivot(muTestPivot),
self.mapParameterListPivot(musPivot),
1e-10 * self.scaleFactorPivot[0])
self._mus = emptyParameterList()
self.mus.reset((self.S, self.npar + len(self.musMargLoc)))
muTestBase = emptyParameterList()
muTestBase.reset((len(muTestPivot), self.npar + len(self.musMargLoc)))
for k in range(self.S):
muk = np.empty_like(self.mus[0])
muk[self.directionPivot] = musPivot[k]
muk[self.directionMarginal] = self.musMargLoc
self.mus[k] = muk
for k in range(len(muTestPivot)):
muk = np.empty_like(muTestBase[0])
muk[self.directionPivot] = muTestPivot[k]
muk[self.directionMarginal] = self.musMargLoc
muTestBase[k] = muk
muTestBase.pop(idxPop)
muLast = copy(self.mus[-1])
self.mus.pop()
if len(self.mus) > 0:
vbMng(self, "MAIN",
("Adding first {} sample point{} at {} to training "
"set.").format(self.S - 1, "" + "s" * (self.S > 2),
self.mus), 3)
self.samplingEngine.iterSample(self.mus)
self._S = len(self.mus)
self._approxParameters["S"] = self.S
self.muTest = parameterList(muTestBase)
self.muTest.append(muLast)
self.M, self.N = ("AUTO",) * 2
def setupApproxLocal(self) -> int:
"""Compute rational interpolant."""
self._marginalizeMiscellanea(True)
setupOK = super().setupApproxLocal()
self._marginalizeMiscellanea(False)
return setupOK
def setupApprox(self, *args, **kwargs) -> int:
"""Compute rational interpolant."""
if self.checkComputedApprox(): return -1
RROMPyAssert(self._mode, message = "Cannot setup approximant.")
vbMng(self, "INIT", "Setting up {}.". format(self.name()), 5)
self.computeScaleFactor()
self._musMarginal = self.samplerMarginal.generatePoints(self.SMarginal)
while len(self.musMarginal) > self.SMarginal: self.musMarginal.pop()
S0 = copy(self.S)
idx, sizes = indicesScatter(len(self.musMarginal), return_sizes = True)
pMat, Ps, Qs, mus = None, [], [], None
req, emptyCores = [], np.where(np.logical_not(sizes))[0]
if len(idx) == 0:
vbMng(self, "MAIN", "Idling.", 25)
if self.storeAllSamples: self.storeSamples()
pL, pT, mT = recv(source = 0, tag = poolRank())
pMat = np.empty((pL, 0), dtype = pT)
mus = np.empty((0, self.mu0.shape[1]), dtype = mT)
else:
_scaleFactorOldPivot = copy(self.scaleFactor)
self.scaleFactor = self.scaleFactorPivot
self._temporaryPivot = 1
for i in idx:
self.musMargLoc = self.musMarginal[i]
vbMng(self, "MAIN",
"Building marginal model no. {} at {}.".format(i + 1,
self.musMargLoc), 5)
self.samplingEngine.resetHistory()
self.trainedModel = None
self.verbosity -= 5
self.samplingEngine.verbosity -= 5
super().setupApprox(*args, **kwargs)
self.verbosity += 5
self.samplingEngine.verbosity += 5
if self.storeAllSamples: self.storeSamples(i)
if pMat is None:
mus = copy(self.samplingEngine.mus.data)
pMat = copy(self.samplingEngine.projectionMatrix)
if i == 0:
for dest in emptyCores:
req += [isend((len(pMat), pMat.dtype, mus.dtype),
dest = dest, tag = dest)]
else:
mus = np.vstack((mus, self.samplingEngine.mus.data))
pMat = np.hstack((pMat,
self.samplingEngine.projectionMatrix))
Ps += [copy(self.trainedModel.data.P)]
Qs += [copy(self.trainedModel.data.Q)]
self._S = S0
del self._temporaryPivot, self.musMargLoc
self.scaleFactor = _scaleFactorOldPivot
for r in req: r.wait()
pMat, Ps, Qs, mus, nsamples = gatherPivotedApproximant(pMat, Ps, Qs,
mus, sizes,
self.polybasis)
self._mus = self.checkParameterList(mus)
Psupp = np.append(0, np.cumsum(nsamples))
self._setupTrainedModel(pMat, forceNew = True)
self.trainedModel.data.Qs, self.trainedModel.data.Ps = Qs, Ps
self.trainedModel.data.Psupp = list(Psupp[: -1])
self._poleMatching()
self._finalizeMarginalization()
vbMng(self, "DEL", "Done setting up approximant.", 5)
return 0
class RationalInterpolantGreedyPivotedNoMatch(
RationalInterpolantGreedyPivotedBase,
GenericPivotedApproximantNoMatch):
"""
ROM pivoted rational interpolant (without pole matching) computation for
parametric problems.
Args:
HFEngine: HF problem solver.
mu0(optional): Default parameter. Defaults to 0.
directionPivot(optional): Pivot components. Defaults to [0].
approxParameters(optional): Dictionary containing values for main
parameters of approximant. Recognized keys are:
- 'POD': kind of snapshots orthogonalization; allowed values
include 0, 1/2, and 1; defaults to 1, i.e. POD;
- 'scaleFactorDer': scaling factors for derivative computation;
defaults to 'AUTO';
- 'S': total number of pivot samples current approximant relies
upon;
- 'samplerPivot': pivot sample point generator;
- 'SMarginal': total number of marginal samples current approximant
relies upon;
- 'samplerMarginal': marginal sample point generator;
- 'polybasis': type of polynomial basis for pivot
interpolation; defaults to 'MONOMIAL';
- 'greedyTol': uniform error tolerance for greedy algorithm;
defaults to 1e-2;
- 'collinearityTol': collinearity tolerance for greedy algorithm;
defaults to 0.;
- 'maxIter': maximum number of greedy steps; defaults to 1e2;
- 'nTestPoints': number of test points; defaults to 5e2;
- 'trainSetGenerator': training sample points generator; defaults
to sampler;
- 'errorEstimatorKind': kind of error estimator; available values
include 'AFFINE', 'DISCREPANCY', 'LOOK_AHEAD',
'LOOK_AHEAD_RES', and 'NONE'; defaults to 'NONE';
- 'radialDirectionalWeightsMarginal': radial basis weights for
marginal interpolant; defaults to 1;
- 'functionalSolve': strategy for minimization of denominator
functional; allowed values include 'NORM', 'DOMINANT', 'NODAL',
- 'LOEWNER', and 'BARYCENTRIC'; defaults to 'NORM';
+ 'BARYCENTRIC_NORM', and 'BARYCENTRIC[_AVERAGE]' (check pdf in
+ main folder for meaning); defaults to 'NORM';
- 'interpRcond': tolerance for pivot interpolation; defaults to
None;
- 'robustTol': tolerance for robust rational denominator
management; defaults to 0.
Defaults to empty dict.
approx_state(optional): Whether to approximate state. Defaults to
False.
verbosity(optional): Verbosity level. Defaults to 10.
Attributes:
HFEngine: HF problem solver.
mu0: Default parameter.
directionPivot: Pivot components.
mus: Array of snapshot parameters.
musMarginal: Array of marginal snapshot parameters.
approxParameters: Dictionary containing values for main parameters of
approximant. Recognized keys are in parameterList.
parameterListSoft: Recognized keys of soft approximant parameters:
- 'POD': kind of snapshots orthogonalization;
- 'scaleFactorDer': scaling factors for derivative computation;
- 'polybasis': type of polynomial basis for pivot
interpolation;
- 'greedyTol': uniform error tolerance for greedy algorithm;
- 'collinearityTol': collinearity tolerance for greedy algorithm;
- 'maxIter': maximum number of greedy steps;
- 'nTestPoints': number of test points;
- 'trainSetGenerator': training sample points generator;
- 'errorEstimatorKind': kind of error estimator;
- 'radialDirectionalWeightsMarginal': radial basis weights for
marginal interpolant;
- 'functionalSolve': strategy for minimization of denominator
functional;
- 'interpRcond': tolerance for pivot interpolation;
- 'robustTol': tolerance for robust rational denominator
management.
parameterListCritical: Recognized keys of critical approximant
parameters:
- 'S': total number of pivot samples current approximant relies
upon;
- 'samplerPivot': pivot sample point generator;
- 'SMarginal': total number of marginal samples current approximant
relies upon;
- 'samplerMarginal': marginal sample point generator.
approx_state: Whether to approximate state.
verbosity: Verbosity level.
POD: Kind of snapshots orthogonalization.
scaleFactorDer: Scaling factors for derivative computation.
S: Total number of pivot samples current approximant relies upon.
samplerPivot: Pivot sample point generator.
SMarginal: Total number of marginal samples current approximant relies
upon.
samplerMarginal: Marginal sample point generator.
polybasis: Type of polynomial basis for pivot interpolation.
greedyTol: uniform error tolerance for greedy algorithm.
collinearityTol: Collinearity tolerance for greedy algorithm.
maxIter: maximum number of greedy steps.
nTestPoints: number of starting training points.
trainSetGenerator: training sample points generator.
errorEstimatorKind: kind of error estimator.
radialDirectionalWeightsMarginal: Radial basis weights for marginal
interpolant.
functionalSolve: Strategy for minimization of denominator functional.
interpRcond: Tolerance for pivot interpolation.
robustTol: Tolerance for robust rational denominator management.
muBounds: list of bounds for pivot parameter values.
muBoundsMarginal: list of bounds for marginal parameter values.
samplingEngine: Sampling engine.
uHF: High fidelity solution(s) with parameter(s) lastSolvedHF as
sampleList.
lastSolvedHF: Parameter(s) corresponding to last computed high fidelity
solution(s) as parameterList.
uApproxReduced: Reduced approximate solution(s) with parameter(s)
lastSolvedApprox as sampleList.
lastSolvedApproxReduced: Parameter(s) corresponding to last computed
reduced approximate solution(s) as parameterList.
uApprox: Approximate solution(s) with parameter(s) lastSolvedApprox as
sampleList.
lastSolvedApprox: Parameter(s) corresponding to last computed
approximate solution(s) as parameterList.
Q: Numpy 1D vector containing complex coefficients of approximant
denominator.
P: Numpy 2D vector whose columns are FE dofs of coefficients of
approximant numerator.
"""
class RationalInterpolantGreedyPivoted(RationalInterpolantGreedyPivotedBase,
GenericPivotedApproximant):
"""
ROM pivoted rational interpolant (with pole matching) computation for
parametric problems.
Args:
HFEngine: HF problem solver.
mu0(optional): Default parameter. Defaults to 0.
directionPivot(optional): Pivot components. Defaults to [0].
approxParameters(optional): Dictionary containing values for main
parameters of approximant. Recognized keys are:
- 'POD': kind of snapshots orthogonalization; allowed values
include 0, 1/2, and 1; defaults to 1, i.e. POD;
- 'scaleFactorDer': scaling factors for derivative computation;
defaults to 'AUTO';
- 'matchingWeight': weight for pole matching optimization; defaults
to 1;
- 'sharedRatio': required ratio of marginal points to share
resonance; defaults to 1.;
- 'S': total number of pivot samples current approximant relies
upon;
- 'samplerPivot': pivot sample point generator;
- 'SMarginal': total number of marginal samples current approximant
relies upon;
- 'samplerMarginal': marginal sample point generator;
- 'polybasis': type of polynomial basis for pivot
interpolation; defaults to 'MONOMIAL';
- 'polybasisMarginal': type of polynomial basis for marginal
interpolation; allowed values include 'MONOMIAL_*',
'CHEBYSHEV_*', 'LEGENDRE_*', 'NEARESTNEIGHBOR', and
'PIECEWISE_LINEAR_*'; defaults to 'MONOMIAL';
- 'paramsMarginal': dictionary of parameters for marginal
interpolation; include:
. 'MMarginal': degree of marginal interpolant; defaults to
'AUTO', i.e. maximum allowed; not for 'NEARESTNEIGHBOR' or
'PIECEWISE_LINEAR_*';
. 'nNeighborsMarginal': number of marginal nearest neighbors;
defaults to 1; only for 'NEARESTNEIGHBOR';
. 'polydegreetypeMarginal': type of polynomial degree for
marginal; defaults to 'TOTAL'; not for 'NEARESTNEIGHBOR' or
'PIECEWISE_LINEAR_*';
. 'interpRcondMarginal': tolerance for marginal interpolation;
defaults to None; not for 'NEARESTNEIGHBOR' or
'PIECEWISE_LINEAR_*';
. 'radialDirectionalWeightsMarginalAdapt': bounds for adaptive
rescaling of marginal radial basis weights; only for
radial basis.
- 'greedyTol': uniform error tolerance for greedy algorithm;
defaults to 1e-2;
- 'collinearityTol': collinearity tolerance for greedy algorithm;
defaults to 0.;
- 'maxIter': maximum number of greedy steps; defaults to 1e2;
- 'nTestPoints': number of test points; defaults to 5e2;
- 'trainSetGenerator': training sample points generator; defaults
to sampler;
- 'errorEstimatorKind': kind of error estimator; available values
include 'AFFINE', 'DISCREPANCY', 'LOOK_AHEAD',
'LOOK_AHEAD_RES', and 'NONE'; defaults to 'NONE';
- 'radialDirectionalWeightsMarginal': radial basis weights for
marginal interpolant; defaults to 1;
- 'functionalSolve': strategy for minimization of denominator
functional; allowed values include 'NORM', 'DOMINANT', 'NODAL',
- 'LOEWNER', and 'BARYCENTRIC'; defaults to 'NORM';
+ 'BARYCENTRIC_NORM', and 'BARYCENTRIC[_AVERAGE]' (check pdf in
+ main folder for meaning); defaults to 'NORM';
- 'interpRcond': tolerance for pivot interpolation; defaults to
None;
- 'robustTol': tolerance for robust rational denominator
management; defaults to 0.
Defaults to empty dict.
approx_state(optional): Whether to approximate state. Defaults to
False.
verbosity(optional): Verbosity level. Defaults to 10.
Attributes:
HFEngine: HF problem solver.
mu0: Default parameter.
directionPivot: Pivot components.
mus: Array of snapshot parameters.
musMarginal: Array of marginal snapshot parameters.
approxParameters: Dictionary containing values for main parameters of
approximant. Recognized keys are in parameterList.
parameterListSoft: Recognized keys of soft approximant parameters:
- 'POD': kind of snapshots orthogonalization;
- 'scaleFactorDer': scaling factors for derivative computation;
- 'matchingWeight': weight for pole matching optimization;
- 'sharedRatio': required ratio of marginal points to share
resonance;
- 'polybasis': type of polynomial basis for pivot
interpolation;
- 'polybasisMarginal': type of polynomial basis for marginal
interpolation;
- 'paramsMarginal': dictionary of parameters for marginal
interpolation; include:
. 'MMarginal': degree of marginal interpolant;
. 'nNeighborsMarginal': number of marginal nearest neighbors;
. 'polydegreetypeMarginal': type of polynomial degree for
marginal;
. 'interpRcondMarginal': tolerance for marginal interpolation;
. 'radialDirectionalWeightsMarginalAdapt': bounds for adaptive
rescaling of marginal radial basis weights.
- 'greedyTol': uniform error tolerance for greedy algorithm;
- 'collinearityTol': collinearity tolerance for greedy algorithm;
- 'maxIter': maximum number of greedy steps;
- 'nTestPoints': number of test points;
- 'trainSetGenerator': training sample points generator;
- 'errorEstimatorKind': kind of error estimator;
- 'radialDirectionalWeightsMarginal': radial basis weights for
marginal interpolant;
- 'functionalSolve': strategy for minimization of denominator
functional;
- 'interpRcond': tolerance for pivot interpolation;
- 'robustTol': tolerance for robust rational denominator
management.
parameterListCritical: Recognized keys of critical approximant
parameters:
- 'S': total number of pivot samples current approximant relies
upon;
- 'samplerPivot': pivot sample point generator;
- 'SMarginal': total number of marginal samples current approximant
relies upon;
- 'samplerMarginal': marginal sample point generator.
approx_state: Whether to approximate state.
verbosity: Verbosity level.
POD: Kind of snapshots orthogonalization.
scaleFactorDer: Scaling factors for derivative computation.
matchingWeight: Weight for pole matching optimization.
sharedRatio: Required ratio of marginal points to share resonance.
S: Total number of pivot samples current approximant relies upon.
samplerPivot: Pivot sample point generator.
SMarginal: Total number of marginal samples current approximant relies
upon.
samplerMarginal: Marginal sample point generator.
polybasis: Type of polynomial basis for pivot interpolation.
polybasisMarginal: Type of polynomial basis for marginal interpolation.
paramsMarginal: Dictionary of parameters for marginal interpolation.
greedyTol: uniform error tolerance for greedy algorithm.
collinearityTol: Collinearity tolerance for greedy algorithm.
maxIter: maximum number of greedy steps.
nTestPoints: number of starting training points.
trainSetGenerator: training sample points generator.
errorEstimatorKind: kind of error estimator.
radialDirectionalWeightsMarginal: Radial basis weights for marginal
interpolant.
functionalSolve: Strategy for minimization of denominator functional.
interpRcond: Tolerance for pivot interpolation.
robustTol: Tolerance for robust rational denominator management.
muBounds: list of bounds for pivot parameter values.
muBoundsMarginal: list of bounds for marginal parameter values.
samplingEngine: Sampling engine.
uHF: High fidelity solution(s) with parameter(s) lastSolvedHF as
sampleList.
lastSolvedHF: Parameter(s) corresponding to last computed high fidelity
solution(s) as parameterList.
uApproxReduced: Reduced approximate solution(s) with parameter(s)
lastSolvedApprox as sampleList.
lastSolvedApproxReduced: Parameter(s) corresponding to last computed
reduced approximate solution(s) as parameterList.
uApprox: Approximate solution(s) with parameter(s) lastSolvedApprox as
sampleList.
lastSolvedApprox: Parameter(s) corresponding to last computed
approximate solution(s) as parameterList.
Q: Numpy 1D vector containing complex coefficients of approximant
denominator.
P: Numpy 2D vector whose columns are FE dofs of coefficients of
approximant numerator.
"""
diff --git a/rrompy/reduction_methods/pivoted/rational_interpolant_pivoted.py b/rrompy/reduction_methods/pivoted/rational_interpolant_pivoted.py
index 583c821..4af819a 100644
--- a/rrompy/reduction_methods/pivoted/rational_interpolant_pivoted.py
+++ b/rrompy/reduction_methods/pivoted/rational_interpolant_pivoted.py
@@ -1,454 +1,456 @@
# 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 collections.abc import Iterable
from copy import deepcopy as copy
from .generic_pivoted_approximant import (GenericPivotedApproximantBase,
GenericPivotedApproximantNoMatch,
GenericPivotedApproximant)
from .gather_pivoted_approximant import gatherPivotedApproximant
from rrompy.reduction_methods.standard.rational_interpolant import (
RationalInterpolant)
from rrompy.utilities.base import verbosityManager as vbMng
from rrompy.utilities.numerical.hash_derivative import nextDerivativeIndices
from rrompy.utilities.exception_manager import RROMPyAssert, RROMPyWarning
from rrompy.parameter import emptyParameterList
from rrompy.utilities.parallel import poolRank, indicesScatter, isend, recv
__all__ = ['RationalInterpolantPivotedNoMatch', 'RationalInterpolantPivoted']
class RationalInterpolantPivotedBase(GenericPivotedApproximantBase,
RationalInterpolant):
def __init__(self, *args, **kwargs):
self._preInit()
self._addParametersToList(toBeExcluded = ["polydegreetype"])
super().__init__(*args, **kwargs)
self._ignoreResidues = self.nparPivot > 1
self._postInit()
@property
def scaleFactorDer(self):
"""Value of scaleFactorDer."""
if self._scaleFactorDer == "NONE": return 1.
if self._scaleFactorDer == "AUTO": return self.scaleFactorPivot
return self._scaleFactorDer
@scaleFactorDer.setter
def scaleFactorDer(self, scaleFactorDer):
if isinstance(scaleFactorDer, (str,)):
scaleFactorDer = scaleFactorDer.upper()
elif isinstance(scaleFactorDer, Iterable):
scaleFactorDer = list(scaleFactorDer)
self._scaleFactorDer = scaleFactorDer
self._approxParameters["scaleFactorDer"] = self._scaleFactorDer
@property
def polydegreetype(self):
"""Value of polydegreetype."""
return "TOTAL"
@polydegreetype.setter
def polydegreetype(self, polydegreetype):
RROMPyWarning(("polydegreetype is used just to simplify inheritance, "
"and its value cannot be changed from 'TOTAL'."))
def _setupInterpolationIndices(self):
"""Setup parameters for polyvander."""
RROMPyAssert(self._mode,
message = "Cannot setup interpolation indices.")
if (self._musUniqueCN is None
or len(self._reorder) != len(self.musPivot)):
try:
muPC = self.trainedModel.centerNormalizePivot(self.musPivot)
except:
muPC = self.trainedModel.centerNormalize(self.musPivot)
self._musUniqueCN, musIdxsTo, musIdxs, musCount = (muPC.unique(
return_index = True, return_inverse = True,
return_counts = True))
self._musUnique = self.musPivot[musIdxsTo]
self._derIdxs = [None] * len(self._musUniqueCN)
self._reorder = np.empty(len(musIdxs), dtype = int)
filled = 0
for j, cnt in enumerate(musCount):
self._derIdxs[j] = nextDerivativeIndices([], self.nparPivot,
cnt)
jIdx = np.nonzero(musIdxs == j)[0]
self._reorder[jIdx] = np.arange(filled, filled + cnt)
filled += cnt
def setupApprox(self) -> int:
"""Compute rational interpolant."""
if self.checkComputedApprox(): return -1
RROMPyAssert(self._mode, message = "Cannot setup approximant.")
vbMng(self, "INIT", "Setting up {}.". format(self.name()), 5)
self.computeScaleFactor()
self.resetSamples()
self.samplingEngine.scaleFactor = self.scaleFactorDer
self.musPivot = self.samplerPivot.generatePoints(self.S)
while len(self.musPivot) > self.S: self.musPivot.pop()
self._musMarginal = self.samplerMarginal.generatePoints(self.SMarginal)
while len(self.musMarginal) > self.SMarginal: self.musMarginal.pop()
self._mus = emptyParameterList()
self.mus.reset((self.S * self.SMarginal, self.HFEngine.npar))
for j, muMarg in enumerate(self.musMarginal):
for k in range(j * self.S, (j + 1) * self.S):
muk = np.empty_like(self.mus[0])
muk[self.directionPivot] = self.musPivot[k - j * self.S]
muk[self.directionMarginal] = muMarg
self.mus[k] = muk
N0 = copy(self.N)
self._setupTrainedModel(np.zeros((0, 0)), forceNew = True)
idx, sizes = indicesScatter(len(self.musMarginal), return_sizes = True)
pMat, Ps, Qs = None, [], []
req, emptyCores = [], np.where(np.logical_not(sizes))[0]
if len(idx) == 0:
vbMng(self, "MAIN", "Idling.", 30)
if self.storeAllSamples: self.storeSamples()
pL, pT = recv(source = 0, tag = poolRank())
pMat = np.empty((pL, 0), dtype = pT)
else:
_scaleFactorOldPivot = copy(self.scaleFactor)
self.scaleFactor = self.scaleFactorPivot
self._temporaryPivot = 1
for i in idx:
vbMng(self, "MAIN",
"Building marginal model no. {} at {}.".format(i + 1,
self.musMarginal[i]), 5)
vbMng(self, "INIT", "Starting computation of snapshots.", 10)
self.samplingEngine.resetHistory()
self.samplingEngine.iterSample(
self.mus[self.S * i : self.S * (i + 1)])
vbMng(self, "DEL", "Done computing snapshots.", 10)
self.verbosity -= 5
self.samplingEngine.verbosity -= 5
self._setupRational(self._setupDenominator())
self.verbosity += 5
self.samplingEngine.verbosity += 5
if self.storeAllSamples: self.storeSamples(i)
if pMat is None:
pMat = copy(self.samplingEngine.projectionMatrix)
if i == 0:
for dest in emptyCores:
req += [isend((len(pMat), pMat.dtype), dest = dest,
tag = dest)]
else:
pMat = np.hstack((pMat,
self.samplingEngine.projectionMatrix))
Ps += [copy(self.trainedModel.data.P)]
Qs += [copy(self.trainedModel.data.Q)]
del self.trainedModel.data.Q, self.trainedModel.data.P
self.N = N0
del self._temporaryPivot
self.scaleFactor = _scaleFactorOldPivot
for r in req: r.wait()
pMat, Ps, Qs, _, _ = gatherPivotedApproximant(pMat, Ps, Qs,
self.mus.data, sizes,
self.polybasis, False)
self._setupTrainedModel(pMat)
self.trainedModel.data.Qs, self.trainedModel.data.Ps = Qs, Ps
Psupp = np.arange(0, len(self.musMarginal) * self.S, self.S)
self.trainedModel.data.Psupp = list(Psupp)
self._poleMatching()
self._finalizeMarginalization()
vbMng(self, "DEL", "Done setting up approximant.", 5)
return 0
class RationalInterpolantPivotedNoMatch(RationalInterpolantPivotedBase,
GenericPivotedApproximantNoMatch):
"""
ROM pivoted rational interpolant (without pole matching) computation for
parametric problems.
Args:
HFEngine: HF problem solver.
mu0(optional): Default parameter. Defaults to 0.
directionPivot(optional): Pivot components. Defaults to [0].
approxParameters(optional): Dictionary containing values for main
parameters of approximant. Recognized keys are:
- 'POD': kind of snapshots orthogonalization; allowed values
include 0, 1/2, and 1; defaults to 1, i.e. POD;
- 'scaleFactorDer': scaling factors for derivative computation;
defaults to 'AUTO';
- 'S': total number of pivot samples current approximant relies
upon;
- 'samplerPivot': pivot sample point generator;
- 'SMarginal': total number of marginal samples current approximant
relies upon;
- 'samplerMarginal': marginal sample point generator;
- 'polybasis': type of polynomial basis for pivot
interpolation; defaults to 'MONOMIAL';
- 'M': degree of rational interpolant numerator; defaults to
'AUTO', i.e. maximum allowed;
- 'N': degree of rational interpolant denominator; defaults to
'AUTO', i.e. maximum allowed;
- 'radialDirectionalWeights': radial basis weights for pivot
numerator; defaults to 1;
- 'radialDirectionalWeightsAdapt': bounds for adaptive rescaling of
radial basis weights; defaults to [-1, -1];
- 'radialDirectionalWeightsMarginal': radial basis weights for
marginal interpolant; defaults to 1;
- 'functionalSolve': strategy for minimization of denominator
functional; allowed values include 'NORM', 'DOMINANT', 'NODAL',
- 'LOEWNER', and 'BARYCENTRIC'; defaults to 'NORM';
+ 'BARYCENTRIC_NORM', and 'BARYCENTRIC[_AVERAGE]' (check pdf in
+ main folder for meaning); defaults to 'NORM';
- 'interpRcond': tolerance for pivot interpolation; defaults to
None;
- 'robustTol': tolerance for robust rational denominator
management; defaults to 0.
Defaults to empty dict.
approx_state(optional): Whether to approximate state. Defaults to
False.
verbosity(optional): Verbosity level. Defaults to 10.
Attributes:
HFEngine: HF problem solver.
mu0: Default parameter.
directionPivot: Pivot components.
mus: Array of snapshot parameters.
musPivot: Array of pivot snapshot parameters.
musMarginal: Array of marginal snapshot parameters.
approxParameters: Dictionary containing values for main parameters of
approximant. Recognized keys are in parameterList.
parameterListSoft: Recognized keys of soft approximant parameters:
- 'POD': kind of snapshots orthogonalization;
- 'scaleFactorDer': scaling factors for derivative computation;
- 'polybasis': type of polynomial basis for pivot
interpolation;
- 'M': degree of rational interpolant numerator;
- 'N': degree of rational interpolant denominator;
- 'radialDirectionalWeights': radial basis weights for pivot
numerator;
- 'radialDirectionalWeightsAdapt': bounds for adaptive rescaling of
radial basis weights;
- 'radialDirectionalWeightsMarginal': radial basis weights for
marginal interpolant;
- 'functionalSolve': strategy for minimization of denominator
functional;
- 'interpRcond': tolerance for pivot interpolation;
- 'robustTol': tolerance for robust rational denominator
management.
parameterListCritical: Recognized keys of critical approximant
parameters:
- 'S': total number of pivot samples current approximant relies
upon;
- 'samplerPivot': pivot sample point generator;
- 'SMarginal': total number of marginal samples current approximant
relies upon;
- 'samplerMarginal': marginal sample point generator.
approx_state: Whether to approximate state.
verbosity: Verbosity level.
POD: Kind of snapshots orthogonalization.
scaleFactorDer: Scaling factors for derivative computation.
S: Total number of pivot samples current approximant relies upon.
samplerPivot: Pivot sample point generator.
SMarginal: Total number of marginal samples current approximant relies
upon.
samplerMarginal: Marginal sample point generator.
polybasis: Type of polynomial basis for pivot interpolation.
M: Numerator degree of approximant.
N: Denominator degree of approximant.
radialDirectionalWeights: Radial basis weights for pivot numerator.
radialDirectionalWeightsAdapt: Bounds for adaptive rescaling of radial
basis weights.
radialDirectionalWeightsMarginal: Radial basis weights for marginal
interpolant.
functionalSolve: Strategy for minimization of denominator functional.
interpRcond: Tolerance for pivot interpolation.
robustTol: Tolerance for robust rational denominator management.
muBounds: list of bounds for pivot parameter values.
muBoundsMarginal: list of bounds for marginal parameter values.
samplingEngine: Sampling engine.
uHF: High fidelity solution(s) with parameter(s) lastSolvedHF as
sampleList.
lastSolvedHF: Parameter(s) corresponding to last computed high fidelity
solution(s) as parameterList.
uApproxReduced: Reduced approximate solution(s) with parameter(s)
lastSolvedApprox as sampleList.
lastSolvedApproxReduced: Parameter(s) corresponding to last computed
reduced approximate solution(s) as parameterList.
uApprox: Approximate solution(s) with parameter(s) lastSolvedApprox as
sampleList.
lastSolvedApprox: Parameter(s) corresponding to last computed
approximate solution(s) as parameterList.
Q: Numpy 1D vector containing complex coefficients of approximant
denominator.
P: Numpy 2D vector whose columns are FE dofs of coefficients of
approximant numerator.
"""
class RationalInterpolantPivoted(RationalInterpolantPivotedBase,
GenericPivotedApproximant):
"""
ROM pivoted rational interpolant (with pole matching) computation for
parametric problems.
Args:
HFEngine: HF problem solver.
mu0(optional): Default parameter. Defaults to 0.
directionPivot(optional): Pivot components. Defaults to [0].
approxParameters(optional): Dictionary containing values for main
parameters of approximant. Recognized keys are:
- 'POD': kind of snapshots orthogonalization; allowed values
include 0, 1/2, and 1; defaults to 1, i.e. POD;
- 'scaleFactorDer': scaling factors for derivative computation;
defaults to 'AUTO';
- 'matchingWeight': weight for pole matching optimization; defaults
to 1;
- 'sharedRatio': required ratio of marginal points to share
resonance; defaults to 1.;
- 'S': total number of pivot samples current approximant relies
upon;
- 'samplerPivot': pivot sample point generator;
- 'SMarginal': total number of marginal samples current approximant
relies upon;
- 'samplerMarginal': marginal sample point generator;
- 'polybasis': type of polynomial basis for pivot
interpolation; defaults to 'MONOMIAL';
- 'polybasisMarginal': type of polynomial basis for marginal
interpolation; allowed values include 'MONOMIAL_*',
'CHEBYSHEV_*', 'LEGENDRE_*', 'NEARESTNEIGHBOR', and
'PIECEWISE_LINEAR_*'; defaults to 'MONOMIAL';
- 'paramsMarginal': dictionary of parameters for marginal
interpolation; include:
. 'MMarginal': degree of marginal interpolant; defaults to
'AUTO', i.e. maximum allowed; not for 'NEARESTNEIGHBOR' or
'PIECEWISE_LINEAR_*';
. 'nNeighborsMarginal': number of marginal nearest neighbors;
defaults to 1; only for 'NEARESTNEIGHBOR';
. 'polydegreetypeMarginal': type of polynomial degree for
marginal; defaults to 'TOTAL'; not for 'NEARESTNEIGHBOR' or
'PIECEWISE_LINEAR_*';
. 'interpRcondMarginal': tolerance for marginal interpolation;
defaults to None; not for 'NEARESTNEIGHBOR' or
'PIECEWISE_LINEAR_*';
. 'radialDirectionalWeightsMarginalAdapt': bounds for adaptive
rescaling of marginal radial basis weights; only for
radial basis.
- 'M': degree of rational interpolant numerator; defaults to
'AUTO', i.e. maximum allowed;
- 'N': degree of rational interpolant denominator; defaults to
'AUTO', i.e. maximum allowed;
- 'radialDirectionalWeights': radial basis weights for pivot
numerator; defaults to 1;
- 'radialDirectionalWeightsAdapt': bounds for adaptive rescaling of
radial basis weights; defaults to [-1, -1];
- 'radialDirectionalWeightsMarginal': radial basis weights for
marginal interpolant; defaults to 1;
- 'functionalSolve': strategy for minimization of denominator
functional; allowed values include 'NORM', 'DOMINANT', 'NODAL',
- 'LOEWNER', and 'BARYCENTRIC'; defaults to 'NORM';
+ 'BARYCENTRIC_NORM', and 'BARYCENTRIC[_AVERAGE]' (check pdf in
+ main folder for meaning); defaults to 'NORM';
- 'interpRcond': tolerance for pivot interpolation; defaults to
None;
- 'robustTol': tolerance for robust rational denominator
management; defaults to 0.
Defaults to empty dict.
approx_state(optional): Whether to approximate state. Defaults to
False.
verbosity(optional): Verbosity level. Defaults to 10.
Attributes:
HFEngine: HF problem solver.
mu0: Default parameter.
directionPivot: Pivot components.
mus: Array of snapshot parameters.
musPivot: Array of pivot snapshot parameters.
musMarginal: Array of marginal snapshot parameters.
approxParameters: Dictionary containing values for main parameters of
approximant. Recognized keys are in parameterList.
parameterListSoft: Recognized keys of soft approximant parameters:
- 'POD': kind of snapshots orthogonalization;
- 'scaleFactorDer': scaling factors for derivative computation;
- 'matchingWeight': weight for pole matching optimization;
- 'sharedRatio': required ratio of marginal points to share
resonance;
- 'polybasis': type of polynomial basis for pivot
interpolation;
- 'polybasisMarginal': type of polynomial basis for marginal
interpolation;
- 'paramsMarginal': dictionary of parameters for marginal
interpolation; include:
. 'MMarginal': degree of marginal interpolant;
. 'nNeighborsMarginal': number of marginal nearest neighbors;
. 'polydegreetypeMarginal': type of polynomial degree for
marginal;
. 'interpRcondMarginal': tolerance for marginal interpolation;
. 'radialDirectionalWeightsMarginalAdapt': bounds for adaptive
rescaling of marginal radial basis weights.
- 'M': degree of rational interpolant numerator;
- 'N': degree of rational interpolant denominator;
- 'radialDirectionalWeights': radial basis weights for pivot
numerator;
- 'radialDirectionalWeightsAdapt': bounds for adaptive rescaling of
radial basis weights;
- 'radialDirectionalWeightsMarginal': radial basis weights for
marginal interpolant;
- 'functionalSolve': strategy for minimization of denominator
functional;
- 'interpRcond': tolerance for pivot interpolation;
- 'robustTol': tolerance for robust rational denominator
management.
parameterListCritical: Recognized keys of critical approximant
parameters:
- 'S': total number of pivot samples current approximant relies
upon;
- 'samplerPivot': pivot sample point generator;
- 'SMarginal': total number of marginal samples current approximant
relies upon;
- 'samplerMarginal': marginal sample point generator.
approx_state: Whether to approximate state.
verbosity: Verbosity level.
POD: Kind of snapshots orthogonalization.
scaleFactorDer: Scaling factors for derivative computation.
matchingWeight: Weight for pole matching optimization.
sharedRatio: Required ratio of marginal points to share resonance.
S: Total number of pivot samples current approximant relies upon.
samplerPivot: Pivot sample point generator.
SMarginal: Total number of marginal samples current approximant relies
upon.
samplerMarginal: Marginal sample point generator.
polybasis: Type of polynomial basis for pivot interpolation.
polybasisMarginal: Type of polynomial basis for marginal interpolation.
paramsMarginal: Dictionary of parameters for marginal interpolation.
M: Numerator degree of approximant.
N: Denominator degree of approximant.
radialDirectionalWeights: Radial basis weights for pivot numerator.
radialDirectionalWeightsAdapt: Bounds for adaptive rescaling of radial
basis weights.
radialDirectionalWeightsMarginal: Radial basis weights for marginal
interpolant.
functionalSolve: Strategy for minimization of denominator functional.
interpRcond: Tolerance for pivot interpolation.
robustTol: Tolerance for robust rational denominator management.
muBounds: list of bounds for pivot parameter values.
muBoundsMarginal: list of bounds for marginal parameter values.
samplingEngine: Sampling engine.
uHF: High fidelity solution(s) with parameter(s) lastSolvedHF as
sampleList.
lastSolvedHF: Parameter(s) corresponding to last computed high fidelity
solution(s) as parameterList.
uApproxReduced: Reduced approximate solution(s) with parameter(s)
lastSolvedApprox as sampleList.
lastSolvedApproxReduced: Parameter(s) corresponding to last computed
reduced approximate solution(s) as parameterList.
uApprox: Approximate solution(s) with parameter(s) lastSolvedApprox as
sampleList.
lastSolvedApprox: Parameter(s) corresponding to last computed
approximate solution(s) as parameterList.
Q: Numpy 1D vector containing complex coefficients of approximant
denominator.
P: Numpy 2D vector whose columns are FE dofs of coefficients of
approximant numerator.
"""
diff --git a/rrompy/reduction_methods/standard/greedy/rational_interpolant_greedy.py b/rrompy/reduction_methods/standard/greedy/rational_interpolant_greedy.py
index 42f82cd..d1a0787 100644
--- a/rrompy/reduction_methods/standard/greedy/rational_interpolant_greedy.py
+++ b/rrompy/reduction_methods/standard/greedy/rational_interpolant_greedy.py
@@ -1,536 +1,537 @@
# Copyright (C) 2018 by the RROMPy authors
#
# This file is part of RROMPy.
#
# RROMPy is free software: you can redistribute it and/or modify
# it under the terms of the GNU Lesser General Public License as published by
# the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
#
# RROMPy is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU Lesser General Public License for more details.
#
# You should have received a copy of the GNU Lesser General Public License
# along with RROMPy. If not, see .
#
from copy import deepcopy as copy
import numpy as np
from rrompy.hfengines.base.linear_affine_engine import checkIfAffine
from .generic_greedy_approximant import GenericGreedyApproximant
from rrompy.utilities.poly_fitting.polynomial import (polybases, polyfitname,
PolynomialInterpolator as PI,
polyvander)
from rrompy.utilities.numerical import dot
from rrompy.utilities.numerical.degree import totalDegreeN
from rrompy.utilities.expression import expressionEvaluator
from rrompy.reduction_methods.standard import RationalInterpolant
from rrompy.utilities.base.types import Np1D, Tuple, paramVal, List
from rrompy.utilities.base.verbosity_depth import (verbosityManager as vbMng,
getVerbosityDepth, setVerbosityDepth)
from rrompy.utilities.poly_fitting import customFit
from rrompy.utilities.exception_manager import (RROMPyWarning, RROMPyException,
RROMPyAssert, RROMPy_FRAGILE)
from rrompy.sampling import sampleList, emptySampleList
__all__ = ['RationalInterpolantGreedy']
class RationalInterpolantGreedy(GenericGreedyApproximant, RationalInterpolant):
"""
ROM greedy rational interpolant computation for parametric problems.
Args:
HFEngine: HF problem solver.
mu0(optional): Default parameter. Defaults to 0.
approxParameters(optional): Dictionary containing values for main
parameters of approximant. Recognized keys are:
- 'POD': kind of snapshots orthogonalization; allowed values
include 0, 1/2, and 1; defaults to 1, i.e. POD;
- 'scaleFactorDer': scaling factors for derivative computation;
defaults to 'AUTO';
- 'S': number of starting training points;
- 'sampler': sample point generator;
- 'greedyTol': uniform error tolerance for greedy algorithm;
defaults to 1e-2;
- 'collinearityTol': collinearity tolerance for greedy algorithm;
defaults to 0.;
- 'maxIter': maximum number of greedy steps; defaults to 1e2;
- 'nTestPoints': number of test points; defaults to 5e2;
- 'trainSetGenerator': training sample points generator; defaults
to sampler;
- 'polybasis': type of basis for interpolation; defaults to
'MONOMIAL';
- 'errorEstimatorKind': kind of error estimator; available values
include 'AFFINE', 'DISCREPANCY', 'LOOK_AHEAD',
'LOOK_AHEAD_RES', 'LOOK_AHEAD_OUTPUT', and 'NONE'; defaults to
'NONE';
- 'functionalSolve': strategy for minimization of denominator
functional; allowed values include 'NORM', 'DOMINANT', 'NODAL',
- 'LOEWNER', and 'BARYCENTRIC'; defaults to 'NORM';
+ 'BARYCENTRIC_NORM', and 'BARYCENTRIC[_AVERAGE]' (check pdf in
+ main folder for meaning); defaults to 'NORM';
- 'interpRcond': tolerance for interpolation; defaults to None;
- 'robustTol': tolerance for robust rational denominator
management; defaults to 0.
Defaults to empty dict.
approx_state(optional): Whether to approximate state. Defaults and must
be True.
verbosity(optional): Verbosity level. Defaults to 10.
Attributes:
HFEngine: HF problem solver.
mu0: Default parameter.
mus: Array of snapshot parameters.
approxParameters: Dictionary containing values for main parameters of
approximant. Recognized keys are in parameterList.
parameterListSoft: Recognized keys of soft approximant parameters:
- 'POD': kind of snapshots orthogonalization;
- 'scaleFactorDer': scaling factors for derivative computation;
- 'greedyTol': uniform error tolerance for greedy algorithm;
- 'collinearityTol': collinearity tolerance for greedy algorithm;
- 'maxIter': maximum number of greedy steps;
- 'nTestPoints': number of test points;
- 'trainSetGenerator': training sample points generator;
- 'errorEstimatorKind': kind of error estimator;
- 'functionalSolve': strategy for minimization of denominator
functional;
- 'interpRcond': tolerance for interpolation;
- 'robustTol': tolerance for robust rational denominator
management.
parameterListCritical: Recognized keys of critical approximant
parameters:
- 'S': total number of samples current approximant relies upon;
- 'sampler': sample point generator.
approx_state: Whether to approximate state.
verbosity: Verbosity level.
POD: Kind of snapshots orthogonalization.
scaleFactorDer: Scaling factors for derivative computation.
S: number of test points.
sampler: Sample point generator.
greedyTol: uniform error tolerance for greedy algorithm.
collinearityTol: Collinearity tolerance for greedy algorithm.
maxIter: maximum number of greedy steps.
nTestPoints: number of starting training points.
trainSetGenerator: training sample points generator.
robustTol: tolerance for robust rational denominator management.
errorEstimatorKind: kind of error estimator.
functionalSolve: Strategy for minimization of denominator functional.
interpRcond: tolerance for interpolation.
robustTol: tolerance for robust rational denominator management.
muBounds: list of bounds for parameter values.
samplingEngine: Sampling engine.
estimatorNormEngine: Engine for estimator norm computation.
uHF: High fidelity solution(s) with parameter(s) lastSolvedHF as
sampleList.
lastSolvedHF: Parameter(s) corresponding to last computed high fidelity
solution(s) as parameterList.
uApproxReduced: Reduced approximate solution(s) with parameter(s)
lastSolvedApprox as sampleList.
lastSolvedApproxReduced: Parameter(s) corresponding to last computed
reduced approximate solution(s) as parameterList.
uApprox: Approximate solution(s) with parameter(s) lastSolvedApprox as
sampleList.
lastSolvedApprox: Parameter(s) corresponding to last computed
approximate solution(s) as parameterList.
"""
_allowedEstimatorKinds = ["AFFINE", "DISCREPANCY", "LOOK_AHEAD",
"LOOK_AHEAD_RES", "LOOK_AHEAD_OUTPUT", "NONE"]
def __init__(self, *args, **kwargs):
self._preInit()
self._addParametersToList(["errorEstimatorKind"], ["DISCREPANCY"],
toBeExcluded = ["M", "N", "polydegreetype",
"radialDirectionalWeights"])
super().__init__(*args, **kwargs)
if not self.approx_state and self.errorEstimatorKind not in [
"LOOK_AHEAD", "LOOK_AHEAD_OUTPUT", "NONE"]:
raise RROMPyException(("Must compute greedy approximation of "
"state, unless error estimator allows "
"otherwise."))
self._postInit()
@property
def approx_state(self):
"""Value of approx_state."""
return self._approx_state
@approx_state.setter
def approx_state(self, approx_state):
RationalInterpolant.approx_state.fset(self, approx_state)
if (not self.approx_state and hasattr(self, "_errorEstimatorKind")
and self.errorEstimatorKind not in [
"LOOK_AHEAD", "LOOK_AHEAD_OUTPUT", "NONE"]):
raise RROMPyException(("Must compute greedy approximation of "
"state, unless error estimator allows "
"otherwise."))
@property
def E(self):
"""Value of E."""
self._E = self.sampleBatchIdx - 1
return self._E
@E.setter
def E(self, E):
RROMPyWarning(("E is used just to simplify inheritance, and its value "
"cannot be changed from that of sampleBatchIdx - 1."))
def _setMAuto(self):
self.M = self.E
def _setNAuto(self):
self.N = self.E
@property
def polydegreetype(self):
"""Value of polydegreetype."""
return "TOTAL"
@polydegreetype.setter
def polydegreetype(self, polydegreetype):
RROMPyWarning(("polydegreetype is used just to simplify inheritance, "
"and its value cannot be changed from 'TOTAL'."))
@property
def polybasis(self):
"""Value of polybasis."""
return self._polybasis
@polybasis.setter
def polybasis(self, polybasis):
try:
polybasis = polybasis.upper().strip().replace(" ","")
if polybasis not in polybases:
raise RROMPyException("Sample type not recognized.")
self._polybasis = polybasis
except:
RROMPyWarning(("Prescribed polybasis not recognized. Overriding "
"to 'MONOMIAL'."))
self._polybasis = "MONOMIAL"
self._approxParameters["polybasis"] = self.polybasis
@property
def errorEstimatorKind(self):
"""Value of errorEstimatorKind."""
return self._errorEstimatorKind
@errorEstimatorKind.setter
def errorEstimatorKind(self, errorEstimatorKind):
errorEstimatorKind = errorEstimatorKind.upper()
if errorEstimatorKind not in self._allowedEstimatorKinds:
RROMPyWarning(("Error estimator kind not recognized. Overriding "
"to 'NONE'."))
errorEstimatorKind = "NONE"
self._errorEstimatorKind = errorEstimatorKind
self._approxParameters["errorEstimatorKind"] = self.errorEstimatorKind
if (self.errorEstimatorKind not in [
"LOOK_AHEAD", "LOOK_AHEAD_OUTPUT", "NONE"]
and hasattr(self, "_approx_state") and not self.approx_state):
raise RROMPyException(("Must compute greedy approximation of "
"state, unless error estimator allows "
"otherwise."))
def _polyvanderAuxiliary(self, mus, deg, *args):
return polyvander(mus, deg, *args)
def getErrorEstimatorDiscrepancy(self, mus:Np1D) -> Np1D:
"""Discrepancy-based residual estimator."""
checkIfAffine(self.HFEngine, "apply discrepancy-based error estimator")
mus = self.checkParameterList(mus)
muCTest = self.trainedModel.centerNormalize(mus)
tMverb, self.trainedModel.verbosity = self.trainedModel.verbosity, 0
QTest = self.trainedModel.getQVal(mus)
QTzero = np.where(QTest == 0.)[0]
if len(QTzero) > 0:
RROMPyWarning(("Adjusting estimator to avoid division by "
"numerically zero denominator."))
QTest[QTzero] = np.finfo(np.complex).eps / (1. + self.N)
self.HFEngine.buildA()
self.HFEngine.buildb()
nAs, nbs = self.HFEngine.nAs, self.HFEngine.nbs
muTrainEff = self.mapParameterList(self.mus)
muTestEff = self.mapParameterList(mus)
PTrain = self.trainedModel.getPVal(self.mus).data.T
QTrain = self.trainedModel.getQVal(self.mus)
QTzero = np.where(QTrain == 0.)[0]
if len(QTzero) > 0:
RROMPyWarning(("Adjusting estimator to avoid division by "
"numerically zero denominator."))
QTrain[QTzero] = np.finfo(np.complex).eps / (1. + self.N)
PTest = self.trainedModel.getPVal(mus).data
self.trainedModel.verbosity = tMverb
radiusAbTrain = np.empty((self.S, nAs * self.S + nbs),
dtype = np.complex)
radiusA = np.empty((self.S, nAs, len(mus)), dtype = np.complex)
radiusb = np.empty((nbs, len(mus)), dtype = np.complex)
for j, thA in enumerate(self.HFEngine.thAs):
idxs = j * self.S + np.arange(self.S)
radiusAbTrain[:, idxs] = expressionEvaluator(thA[0], muTrainEff,
(self.S, 1)) * PTrain
radiusA[:, j] = PTest * expressionEvaluator(thA[0], muTestEff,
(len(mus),))
for j, thb in enumerate(self.HFEngine.thbs):
idx = nAs * self.S + j
radiusAbTrain[:, idx] = QTrain * expressionEvaluator(thb[0],
muTrainEff, (self.S,))
radiusb[j] = QTest * expressionEvaluator(thb[0], muTestEff,
(len(mus),))
QRHSNorm2 = self._affineResidualMatricesContraction(radiusb)
vanTrain = self._polyvanderAuxiliary(self._musUniqueCN, self.E,
self.polybasis0, self._derIdxs,
self._reorder)
interpPQ = customFit(vanTrain, radiusAbTrain,
rcond = self.interpRcond)
vanTest = self._polyvanderAuxiliary(muCTest, self.E, self.polybasis0)
DradiusAb = vanTest.dot(interpPQ)
radiusA = (radiusA
- DradiusAb[:, : - nbs].reshape(len(mus), -1, self.S).T)
radiusb = radiusb - DradiusAb[:, - nbs :].T
ff, Lf, LL = self._affineResidualMatricesContraction(radiusb, radiusA)
err = np.abs((LL - 2. * np.real(Lf) + ff) / QRHSNorm2) ** .5
return err
def getErrorEstimatorLookAhead(self, mus:Np1D,
what : str = "") -> Tuple[Np1D, List[int]]:
"""Residual estimator based on look-ahead idea."""
errTest, QTest, idxMaxEst = self._EIMStep(mus)
_approx_state_old = self.approx_state
if what == "OUTPUT" and _approx_state_old: self._approx_state = False
self.initEstimatorNormEngine()
self._approx_state = _approx_state_old
mu_muTestSample = mus[idxMaxEst]
app_muTestSample = self.getApproxReduced(mu_muTestSample)
if self._mode == RROMPy_FRAGILE:
if what == "RES" and not self.HFEngine.isCEye:
raise RROMPyException(("Cannot compute LOOK_AHEAD_RES "
"estimator in fragile mode for "
"non-scalar C."))
app_muTestSample = dot(self.trainedModel.data.projMat[:,
: app_muTestSample.shape[0]],
app_muTestSample)
else:
app_muTestSample = dot(self.samplingEngine.projectionMatrix,
app_muTestSample)
if what == "RES":
errmu = self.HFEngine.residual(mu_muTestSample, app_muTestSample,
post_c = False)
solmu = self.HFEngine.residual(mu_muTestSample, None,
post_c = False)
else:
for j, mu in enumerate(mu_muTestSample):
uEx = self.samplingEngine.nextSample(mu)
if j == 0:
solmu = emptySampleList()
solmu.reset((len(uEx), len(mu_muTestSample)),
dtype = uEx.dtype)
solmu[j] = uEx
if what == "OUTPUT" and self.approx_state:
solmu = sampleList(self.HFEngine.applyC(solmu))
app_muTestSample = sampleList(self.HFEngine.applyC(
app_muTestSample))
errmu = solmu - app_muTestSample
errsamples = (self.estimatorNormEngine.norm(errmu)
/ self.estimatorNormEngine.norm(solmu))
musT = copy(self.mus)
musT.append(mu_muTestSample)
musT = self.trainedModel.centerNormalize(musT)
musC = self.trainedModel.centerNormalize(mus)
errT = np.zeros((len(musT), len(mu_muTestSample)), dtype = np.complex)
errT[np.arange(len(self.mus), len(musT)),
np.arange(len(mu_muTestSample))] = errsamples * QTest[idxMaxEst]
vanT = self._polyvanderAuxiliary(musT, self.E + 1, self.polybasis)
fitOut = customFit(vanT, errT, full = True, rcond = self.interpRcond)
vbMng(self, "MAIN",
("Fitting {} samples with degree {} through {}... Conditioning "
"of LS system: {:.4e}.").format(len(vanT), self.E + 1,
polyfitname(self.polybasis),
fitOut[1][2][0] / fitOut[1][2][-1]), 15)
vanC = self._polyvanderAuxiliary(musC, self.E + 1, self.polybasis)
err = np.sum(np.abs(vanC.dot(fitOut[0])), axis = -1) / QTest
return err, idxMaxEst
def getErrorEstimatorNone(self, mus:Np1D) -> Np1D:
"""EIM-based residual estimator."""
err = np.max(self._EIMStep(mus, True), axis = 1)
err *= self.greedyTol / np.mean(err)
return err
def _EIMStep(self, mus:Np1D,
only_one : bool = False) -> Tuple[Np1D, Np1D, List[int]]:
"""Residual estimator based on look-ahead idea."""
mus = self.checkParameterList(mus)
tMverb, self.trainedModel.verbosity = self.trainedModel.verbosity, 0
QTest = self.trainedModel.getQVal(mus)
QTzero = np.where(QTest == 0.)[0]
if len(QTzero) > 0:
RROMPyWarning(("Adjusting estimator to avoid division by "
"numerically zero denominator."))
QTest[QTzero] = np.finfo(np.complex).eps / (1. + self.N)
QTest = np.abs(QTest)
muCTest = self.trainedModel.centerNormalize(mus)
muCTrain = self.trainedModel.centerNormalize(self.mus)
self.trainedModel.verbosity = tMverb
vanTest = self._polyvanderAuxiliary(muCTest, self.E, self.polybasis)
vanTestNext = self._polyvanderAuxiliary(muCTest, self.E + 1,
self.polybasis)[:,
vanTest.shape[1] :]
idxsTest = np.arange(vanTestNext.shape[1])
basis = np.zeros((len(idxsTest), 0), dtype = float)
idxMaxEst = []
while len(idxsTest) > 0:
vanTrial = self._polyvanderAuxiliary(muCTrain, self.E,
self.polybasis)
vanTrialNext = self._polyvanderAuxiliary(muCTrain, self.E + 1,
self.polybasis)[:,
vanTrial.shape[1] :]
vanTrial = np.hstack((vanTrial, vanTrialNext.dot(basis).reshape(
len(vanTrialNext), basis.shape[1])))
valuesTrial = vanTrialNext[:, idxsTest]
vanTestEff = np.hstack((vanTest, vanTestNext.dot(basis).reshape(
len(vanTestNext), basis.shape[1])))
vanTestNextEff = vanTestNext[:, idxsTest]
coeffTest = np.linalg.solve(vanTrial, valuesTrial)
errTest = (np.abs(vanTestNextEff - vanTestEff.dot(coeffTest))
/ np.expand_dims(QTest, 1))
if only_one: return errTest
idxMaxErr = np.unravel_index(np.argmax(errTest), errTest.shape)
idxMaxEst += [idxMaxErr[0]]
muCTrain.append(muCTest[idxMaxErr[0]])
basis = np.pad(basis, [(0, 0), (0, 1)], "constant")
basis[idxsTest[idxMaxErr[1]], -1] = 1.
idxsTest = np.delete(idxsTest, idxMaxErr[1])
return errTest, QTest, idxMaxEst
def errorEstimator(self, mus:Np1D, return_max : bool = False) -> Np1D:
"""Standard residual-based error estimator."""
setupOK = self.setupApproxLocal()
if setupOK > 0:
err = np.empty(len(mus))
err[:] = np.nan
if not return_max: return err
return err, [- setupOK], np.nan
mus = self.checkParameterList(mus)
vbMng(self.trainedModel, "INIT",
"Evaluating error estimator at mu = {}.".format(mus), 10)
if self.errorEstimatorKind == "AFFINE":
err = self.getErrorEstimatorAffine(mus)
else:
self._setupInterpolationIndices()
if self.errorEstimatorKind == "DISCREPANCY":
err = self.getErrorEstimatorDiscrepancy(mus)
elif self.errorEstimatorKind[: 10] == "LOOK_AHEAD":
err, idxMaxEst = self.getErrorEstimatorLookAhead(mus,
self.errorEstimatorKind[11 :])
else: #if self.errorEstimatorKind == "NONE":
err = self.getErrorEstimatorNone(mus)
vbMng(self.trainedModel, "DEL", "Done evaluating error estimator", 10)
if not return_max: return err
if self.errorEstimatorKind[: 10] != "LOOK_AHEAD":
idxMaxEst = np.empty(self.sampleBatchSize, dtype = int)
errCP = copy(err)
for j in range(self.sampleBatchSize):
k = np.argmax(errCP)
idxMaxEst[j] = k
if j + 1 < self.sampleBatchSize:
musZero = self.trainedModel.centerNormalize(mus, mus[k])
errCP *= np.linalg.norm(musZero.data, axis = 1)
return err, idxMaxEst, err[idxMaxEst]
def plotEstimator(self, *args, **kwargs):
super().plotEstimator(*args, **kwargs)
if self.errorEstimatorKind == "NONE":
vbMng(self, "MAIN",
("Warning! Error estimator has been arbitrarily normalized "
"before plotting."), 15)
def greedyNextSample(self, *args,
**kwargs) -> Tuple[Np1D, int, float, paramVal]:
"""Compute next greedy snapshot of solution map."""
RROMPyAssert(self._mode, message = "Cannot add greedy sample.")
self.sampleBatchIdx += 1
self.sampleBatchSize = totalDegreeN(self.npar - 1, self.sampleBatchIdx)
err, muidx, maxErr, muNext = super().greedyNextSample(*args, **kwargs)
if maxErr is not None and (np.any(np.isnan(maxErr))
or np.any(np.isinf(maxErr))):
self.sampleBatchIdx -= 1
self.sampleBatchSize = totalDegreeN(self.npar - 1,
self.sampleBatchIdx)
if (self.errorEstimatorKind == "NONE" and not np.isnan(maxErr)
and not np.isinf(maxErr)):
maxErr = None
return err, muidx, maxErr, muNext
def _setSampleBatch(self, maxS:int):
self.sampleBatchIdx, self.sampleBatchSize, S = -1, 0, 0
nextBatchSize = 1
while S + nextBatchSize <= maxS:
self.sampleBatchIdx += 1
self.sampleBatchSize = nextBatchSize
S += self.sampleBatchSize
nextBatchSize = totalDegreeN(self.npar - 1,
self.sampleBatchIdx + 1)
return S
def _preliminaryTraining(self):
"""Initialize starting snapshots of solution map."""
RROMPyAssert(self._mode, message = "Cannot start greedy algorithm.")
if self.samplingEngine.nsamples > 0: return
self._S = self._setSampleBatch(self.S)
super()._preliminaryTraining()
self.M, self.N = ("AUTO",) * 2
def setupApproxLocal(self) -> int:
"""Compute rational interpolant."""
if self.checkComputedApprox(): return -1
RROMPyAssert(self._mode, message = "Cannot setup approximant.")
self.verbosity -= 10
vbMng(self, "INIT", "Setting up local approximant.", 5)
pMat = self.samplingEngine.projectionMatrix
if self.trainedModel is not None:
pMat = pMat[:, len(self.trainedModel.data.mus) :]
self._setupTrainedModel(pMat, self.trainedModel is not None)
self.catchInstability = 2
vbDepth = getVerbosityDepth()
unstable = 0
if self.E > 0:
try:
Q = self._setupDenominator()
except RROMPyException as RE:
if RE.critical: raise RE from None
setVerbosityDepth(vbDepth)
RROMPyWarning("Downgraded {}: {}".format(RE.__class__.__name__,
RE))
unstable = 1
else:
Q = PI()
Q.coeffs = np.ones((1,) * self.npar, dtype = np.complex)
Q.npar = self.npar
Q.polybasis = self.polybasis
if not unstable:
self.trainedModel.data.Q = copy(Q)
try:
P = copy(self._setupNumerator())
except RROMPyException as RE:
if RE.critical: raise RE from None
setVerbosityDepth(vbDepth)
RROMPyWarning("Downgraded {}: {}".format(RE.__class__.__name__,
RE))
unstable = 1
if not unstable:
self.trainedModel.data.P = copy(P)
self.trainedModel.data.approxParameters = copy(
self.approxParameters)
vbMng(self, "DEL", "Done setting up local approximant.", 5)
self.catchInstability = 0
self.verbosity += 10
return unstable
def setupApprox(self, plotEst : str = "NONE") -> int:
val = super().setupApprox(plotEst)
if val == 0:
self._setupRational(self.trainedModel.data.Q,
self.trainedModel.data.P)
self.trainedModel.data.approxParameters = copy(
self.approxParameters)
return val
def loadTrainedModel(self, filename:str):
"""Load trained reduced model from file."""
super().loadTrainedModel(filename)
self._setSampleBatch(self.S + 1)
diff --git a/rrompy/reduction_methods/standard/rational_interpolant.py b/rrompy/reduction_methods/standard/rational_interpolant.py
index 6847364..e62a496 100644
--- a/rrompy/reduction_methods/standard/rational_interpolant.py
+++ b/rrompy/reduction_methods/standard/rational_interpolant.py
@@ -1,873 +1,840 @@
# 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.linalg import eigvals
from collections.abc import Iterable
from .generic_standard_approximant import GenericStandardApproximant
from rrompy.utilities.poly_fitting.polynomial import (
polybases as ppb, polyfitname,
polyvander as pvP, polyTimes,
PolynomialInterpolator as PI,
PolynomialInterpolatorNodal as PIN)
from rrompy.utilities.poly_fitting.heaviside import rational2heaviside
from rrompy.utilities.poly_fitting.radial_basis import (polybases as rbpb,
RadialBasisInterpolator as RBI)
from rrompy.utilities.base.types import (Np1D, Np2D, Tuple, List, sampList,
interpEng)
from rrompy.utilities.base import verbosityManager as vbMng
-from rrompy.utilities.numerical import pseudoInverse, dot, potential
+from rrompy.utilities.numerical import (pseudoInverse, dot, potential,
+ distanceMatrix)
from rrompy.utilities.numerical.factorials import multifactorial
from rrompy.utilities.numerical.hash_derivative import (nextDerivativeIndices,
hashDerivativeToIdx as hashD,
hashIdxToDerivative as hashI)
from rrompy.utilities.numerical.degree import (reduceDegreeN,
degreeTotalToFull, fullDegreeMaxMask,
totalDegreeMaxMask)
from rrompy.solver import Np2DLikeGramian
from rrompy.utilities.exception_manager import (RROMPyException, RROMPyAssert,
RROMPyWarning)
__all__ = ['RationalInterpolant']
def polyTimesTable(P:interpEng, mus:Np1D, reorder:List[int],
derIdxs:List[List[List[int]]], scl : Np1D = None) -> Np2D:
"""Table of polynomial products."""
if not isinstance(P, PI):
raise RROMPyException(("Polynomial to evaluate must be a polynomial "
"interpolator."))
Pvals = [[0.] * len(derIdx) for derIdx in derIdxs]
for j, derIdx in enumerate(derIdxs):
nder = len(derIdx)
for der in range(nder):
derI = hashI(der, P.npar)
Pvals[j][der] = P([mus[j]], derI, scl) / multifactorial(derI)
return blockDiagDer(Pvals, reorder, derIdxs)
def vanderInvTable(vanInv:Np2D, idxs:List[int], reorder:List[int],
derIdxs:List[List[List[int]]]) -> Np2D:
"""Table of Vandermonde pseudo-inverse."""
S = len(reorder)
Ts = [None] * len(idxs)
for k in range(len(idxs)):
invLocs = [None] * len(derIdxs)
idxGlob = 0
for j, derIdx in enumerate(derIdxs):
nder = len(derIdx)
idxGlob += nder
idxLoc = np.arange(S)[np.logical_and(reorder >= idxGlob - nder,
reorder < idxGlob)]
invLocs[j] = vanInv[k, idxLoc]
Ts[k] = blockDiagDer(invLocs, reorder, derIdxs, [2, 1, 0])
return Ts
def blockDiagDer(vals:List[Np1D], reorder:List[int],
derIdxs:List[List[List[int]]],
permute : List[int] = None) -> Np2D:
"""Table of derivative values for point confluence."""
S = len(reorder)
T = np.zeros((S, S), dtype = np.complex)
if permute is None: permute = [0, 1, 2]
idxGlob = 0
for j, derIdx in enumerate(derIdxs):
nder = len(derIdx)
idxGlob += nder
idxLoc = np.arange(S)[np.logical_and(reorder >= idxGlob - nder,
reorder < idxGlob)]
val = vals[j]
for derI, derIdxI in enumerate(derIdx):
for derJ, derIdxJ in enumerate(derIdx):
diffIdx = [x - y for (x, y) in zip(derIdxI, derIdxJ)]
if all([x >= 0 for x in diffIdx]):
diffj = hashD(diffIdx)
i1, i2, i3 = np.array([derI, derJ, diffj])[permute]
T[idxLoc[i1], idxLoc[i2]] = val[i3]
return T
class RationalInterpolant(GenericStandardApproximant):
"""
ROM rational interpolant computation for parametric problems.
Args:
HFEngine: HF problem solver.
mu0(optional): Default parameter. Defaults to 0.
approxParameters(optional): Dictionary containing values for main
parameters of approximant. Recognized keys are:
- 'POD': kind of snapshots orthogonalization; allowed values
include 0, 1/2, and 1; defaults to 1, i.e. POD;
- 'scaleFactorDer': scaling factors for derivative computation;
defaults to 'AUTO';
- 'S': total number of samples current approximant relies upon;
- 'sampler': sample point generator;
- 'polybasis': type of polynomial basis for interpolation; defaults
to 'MONOMIAL';
- 'M': degree of rational interpolant numerator; defaults to
'AUTO', i.e. maximum allowed;
- 'N': degree of rational interpolant denominator; defaults to
'AUTO', i.e. maximum allowed;
- 'polydegreetype': type of polynomial degree; defaults to 'TOTAL';
- 'radialDirectionalWeights': radial basis weights for interpolant
numerator; defaults to 1;
- 'radialDirectionalWeightsAdapt': bounds for adaptive rescaling of
radial basis weights; defaults to [-1, -1];
- 'functionalSolve': strategy for minimization of denominator
functional; allowed values include 'NORM', 'DOMINANT', 'NODAL',
- 'LOEWNER', and 'BARYCENTRIC' (check below for meaning);
- defaults to 'NORM';
+ 'BARYCENTRIC_NORM', and 'BARYCENTRIC[_AVERAGE]' (check pdf in
+ main folder for explanation); defaults to 'NORM';
- 'interpRcond': tolerance for interpolation; defaults to None;
- 'robustTol': tolerance for robust rational denominator
management; defaults to 0.
Defaults to empty dict.
approx_state(optional): Whether to approximate state. Defaults to
False.
verbosity(optional): Verbosity level. Defaults to 10.
Attributes:
HFEngine: HF problem solver.
mu0: Default parameter.
mus: Array of snapshot parameters.
approxParameters: Dictionary containing values for main parameters of
approximant. Recognized keys are in parameterList.
parameterListSoft: Recognized keys of soft approximant parameters:
- 'POD': kind of snapshots orthogonalization;
- 'scaleFactorDer': scaling factors for derivative computation;
- 'polybasis': type of polynomial basis for interpolation;
- 'M': degree of rational interpolant numerator;
- 'N': degree of rational interpolant denominator;
- 'polydegreetype': type of polynomial degree;
- 'radialDirectionalWeights': radial basis weights for interpolant
numerator;
- 'radialDirectionalWeightsAdapt': bounds for adaptive rescaling of
radial basis weights;
- 'functionalSolve': strategy for minimization of denominator
functional;
- 'interpRcond': tolerance for interpolation via numpy.polyfit;
- 'robustTol': tolerance for robust rational denominator
management.
parameterListCritical: Recognized keys of critical approximant
parameters:
- 'S': total number of samples current approximant relies upon;
- 'sampler': sample point generator.
approx_state: Whether to approximate state.
verbosity: Verbosity level.
POD: Kind of snapshots orthogonalization.
scaleFactorDer: Scaling factors for derivative computation.
S: Number of solution snapshots over which current approximant is
based upon.
sampler: Sample point generator.
polybasis: type of polynomial basis for interpolation.
M: Numerator degree of approximant.
N: Denominator degree of approximant.
polydegreetype: Type of polynomial degree.
radialDirectionalWeights: Radial basis weights for interpolant
numerator.
radialDirectionalWeightsAdapt: Bounds for adaptive rescaling of radial
basis weights.
functionalSolve: Strategy for minimization of denominator functional.
interpRcond: Tolerance for interpolation via numpy.polyfit.
robustTol: Tolerance for robust rational denominator management.
muBounds: list of bounds for parameter values.
samplingEngine: Sampling engine.
uHF: High fidelity solution(s) with parameter(s) lastSolvedHF as
sampleList.
lastSolvedHF: Parameter(s) corresponding to last computed high fidelity
solution(s) as parameterList.
uApproxReduced: Reduced approximate solution(s) with parameter(s)
lastSolvedApprox as sampleList.
lastSolvedApproxReduced: Parameter(s) corresponding to last computed
reduced approximate solution(s) as parameterList.
uApprox: Approximate solution(s) with parameter(s) lastSolvedApprox as
sampleList.
lastSolvedApprox: Parameter(s) corresponding to last computed
approximate solution(s) as parameterList.
Q: Numpy 1D vector containing complex coefficients of approximant
denominator.
P: Numpy 2D vector whose columns are FE dofs of coefficients of
approximant numerator.
"""
+
+ _allowedFunctionalSolveKinds = ["NORM", "DOMINANT", "NODAL",
+ "BARYCENTRIC_NORM", "BARYCENTRIC_AVERAGE"]
def __init__(self, *args, **kwargs):
self._preInit()
self._addParametersToList(["polybasis", "M", "N", "polydegreetype",
"radialDirectionalWeights",
"radialDirectionalWeightsAdapt",
"functionalSolve", "interpRcond",
"robustTol"],
["MONOMIAL", "AUTO", "AUTO", "TOTAL", 1.,
[-1., -1.], "NORM", -1, 0.])
super().__init__(*args, **kwargs)
self.catchInstability = 0
self._postInit()
@property
def tModelType(self):
from .trained_model.trained_model_rational import TrainedModelRational
return TrainedModelRational
@property
def polybasis(self):
"""Value of polybasis."""
return self._polybasis
@polybasis.setter
def polybasis(self, polybasis):
try:
polybasis = polybasis.upper().strip().replace(" ","")
if polybasis not in ppb + rbpb:
raise RROMPyException("Prescribed polybasis not recognized.")
self._polybasis = polybasis
except:
RROMPyWarning(("Prescribed polybasis not recognized. Overriding "
"to 'MONOMIAL'."))
self._polybasis = "MONOMIAL"
self._approxParameters["polybasis"] = self.polybasis
@property
def polybasis0(self):
if "_" in self.polybasis:
return self.polybasis.split("_")[0]
return self.polybasis
@property
def functionalSolve(self):
"""Value of functionalSolve."""
return self._functionalSolve
@functionalSolve.setter
def functionalSolve(self, functionalSolve):
try:
functionalSolve = functionalSolve.upper().strip().replace(" ","")
- if functionalSolve not in ["NORM", "DOMINANT", "NODAL", "LOEWNER",
- "BARYCENTRIC"]:
+ if functionalSolve == "BARYCENTRIC": functionalSolve += "_AVERAGE"
+ if functionalSolve not in self._allowedFunctionalSolveKinds:
raise RROMPyException(("Prescribed functionalSolve not "
"recognized."))
self._functionalSolve = functionalSolve
except:
RROMPyWarning(("Prescribed functionalSolve not recognized. "
"Overriding to 'NORM'."))
self._functionalSolve = "NORM"
self._approxParameters["functionalSolve"] = self.functionalSolve
@property
def interpRcond(self):
"""Value of interpRcond."""
return self._interpRcond
@interpRcond.setter
def interpRcond(self, interpRcond):
self._interpRcond = interpRcond
self._approxParameters["interpRcond"] = self.interpRcond
@property
def radialDirectionalWeights(self):
"""Value of radialDirectionalWeights."""
return self._radialDirectionalWeights
@radialDirectionalWeights.setter
def radialDirectionalWeights(self, radialDirectionalWeights):
if isinstance(radialDirectionalWeights, Iterable):
radialDirectionalWeights = list(radialDirectionalWeights)
else:
radialDirectionalWeights = [radialDirectionalWeights]
self._radialDirectionalWeights = radialDirectionalWeights
self._approxParameters["radialDirectionalWeights"] = (
self.radialDirectionalWeights)
@property
def radialDirectionalWeightsAdapt(self):
"""Value of radialDirectionalWeightsAdapt."""
return self._radialDirectionalWeightsAdapt
@radialDirectionalWeightsAdapt.setter
def radialDirectionalWeightsAdapt(self, radialDirectionalWeightsAdapt):
self._radialDirectionalWeightsAdapt = radialDirectionalWeightsAdapt
self._approxParameters["radialDirectionalWeightsAdapt"] = (
self.radialDirectionalWeightsAdapt)
@property
def M(self):
"""Value of M."""
return self._M
@M.setter
def M(self, M):
if isinstance(M, str):
M = M.strip().replace(" ","")
if "-" not in M: M = M + "-0"
self._M_isauto, self._M_shift = True, int(M.split("-")[-1])
M = 0
if M < 0: raise RROMPyException("M must be non-negative.")
self._M = M
self._approxParameters["M"] = self.M
def _setMAuto(self):
self.M = max(0, reduceDegreeN(self.S, self.S, self.npar,
self.polydegreetype) - self._M_shift)
vbMng(self, "MAIN", "Automatically setting M to {}.".format(self.M),
25)
@property
def N(self):
"""Value of N."""
return self._N
@N.setter
def N(self, N):
if isinstance(N, str):
N = N.strip().replace(" ","")
if "-" not in N: N = N + "-0"
self._N_isauto, self._N_shift = True, int(N.split("-")[-1])
N = 0
if N < 0: raise RROMPyException("N must be non-negative.")
self._N = N
self._approxParameters["N"] = self.N
def _setNAuto(self):
self.N = max(0, reduceDegreeN(self.S, self.S, self.npar,
self.polydegreetype) - self._N_shift)
vbMng(self, "MAIN", "Automatically setting N to {}.".format(self.N),
25)
@property
def polydegreetype(self):
"""Value of polydegreetype."""
return self._polydegreetype
@polydegreetype.setter
def polydegreetype(self, polydegreetype):
try:
polydegreetype = polydegreetype.upper().strip().replace(" ","")
if polydegreetype not in ["TOTAL", "FULL"]:
raise RROMPyException(("Prescribed polydegreetype not "
"recognized."))
self._polydegreetype = polydegreetype
except:
RROMPyWarning(("Prescribed polydegreetype not recognized. "
"Overriding to 'TOTAL'."))
self._polydegreetype = "TOTAL"
self._approxParameters["polydegreetype"] = self.polydegreetype
@property
def robustTol(self):
"""Value of tolerance for robust rational denominator management."""
return self._robustTol
@robustTol.setter
def robustTol(self, robustTol):
if robustTol < 0.:
RROMPyWarning(("Overriding prescribed negative robustness "
"tolerance to 0."))
robustTol = 0.
self._robustTol = robustTol
self._approxParameters["robustTol"] = self.robustTol
def resetSamples(self):
"""Reset samples."""
super().resetSamples()
self._musUniqueCN = None
self._derIdxs = None
self._reorder = None
def _setupInterpolationIndices(self):
"""Setup parameters for polyvander."""
if self._musUniqueCN is None or len(self._reorder) != len(self.mus):
self._musUniqueCN, musIdxsTo, musIdxs, musCount = (
self.trainedModel.centerNormalize(self.mus).unique(
return_index = True, return_inverse = True,
return_counts = True))
self._musUnique = self.mus[musIdxsTo]
self._derIdxs = [None] * len(self._musUniqueCN)
self._reorder = np.empty(len(musIdxs), dtype = int)
filled = 0
for j, cnt in enumerate(musCount):
self._derIdxs[j] = nextDerivativeIndices([], self.mus.shape[1],
cnt)
jIdx = np.nonzero(musIdxs == j)[0]
self._reorder[jIdx] = np.arange(filled, filled + cnt)
filled += cnt
def _setupDenominator(self):
"""Compute rational denominator."""
RROMPyAssert(self._mode, message = "Cannot setup denominator.")
vbMng(self, "INIT", "Starting computation of denominator.", 7)
if hasattr(self, "_N_isauto"):
self._setNAuto()
else:
N = reduceDegreeN(self.N, self.S, self.npar, self.polydegreetype)
if N < self.N:
RROMPyWarning(("N too large compared to S. Reducing N by "
"{}").format(self.N - N))
self.N = N
while self.N > 0:
if self.functionalSolve != "NORM" and self.npar > 1:
RROMPyWarning(("Strategy for functional optimization must be "
"'NORM' for more than one parameter. "
"Overriding to 'NORM'."))
self.functionalSolve = "NORM"
- if self.functionalSolve == "BARYCENTRIC" and self.N + 1 < self.S:
+ if (self.functionalSolve[:11] == "BARYCENTRIC"
+ and self.N + 1 < self.S):
RROMPyWarning(("Barycentric strategy cannot be applied with "
"Least Squares. Overriding to 'NORM'."))
self.functionalSolve = "NORM"
- if self.functionalSolve == "BARYCENTRIC":
+ if self.functionalSolve[:11] == "BARYCENTRIC":
invD, TN = None, None
self._setupInterpolationIndices()
else:
invD, TN = self._computeInterpolantInverseBlocks()
- if (self.functionalSolve in ["NODAL", "LOEWNER", "BARYCENTRIC"]
+ if (self.functionalSolve in ["NODAL", "BARYCENTRIC_NORM",
+ "BARYCENTRIC_AVERAGE"]
and len(self._musUnique) != len(self.mus)):
- if self.functionalSolve == "BARYCENTRIC":
+ if self.functionalSolve[:11] == "BARYCENTRIC":
warnflag = "Barycentric"
+ invD, TN = self._computeInterpolantInverseBlocks()
else:
warnflag = "Iterative"
RROMPyWarning(("{} functional optimization cannot be applied "
"to repeated samples. Overriding to "
"'NORM'.").format(warnflag))
self.functionalSolve = "NORM"
idxSamplesEff = list(range(self.S))
if self.POD == 1:
ev, eV = self.findeveVGQR(
self.samplingEngine.Rscale[:, idxSamplesEff], invD, TN)
elif self.POD == 1/2:
ev, eV = self.findeveVGExplicit(
self.samplingEngine.samples_normal(idxSamplesEff),
invD, TN, self.samplingEngine.Rscale)
else:
ev, eV = self.findeveVGExplicit(
self.samplingEngine.samples(idxSamplesEff), invD, TN)
- if self.functionalSolve in ["NODAL", "LOEWNER"]: break
+ if self.functionalSolve == "NODAL": break
evR = ev / np.max(ev)
ts = self.robustTol * np.linalg.norm(evR)
nevBad = len(ev) - np.sum(np.abs(evR) >= ts)
if nevBad <= (self.functionalSolve == "NORM"): break
if self.catchInstability:
raise RROMPyException(("Instability in denominator "
"computation: eigenproblem is poorly "
"conditioned."),
self.catchInstability == 1)
vbMng(self, "MAIN", ("Smallest {} eigenvalues below tolerance. "
"Reducing N by 1.").format(nevBad), 10)
self.N = self.N - 1
if self.N <= 0:
self.N = 0
eV = np.ones((1, 1))
- if self.N > 0 and self.functionalSolve in ["NODAL", "LOEWNER",
- "BARYCENTRIC"]:
+ if self.N > 0 and self.functionalSolve in ["NODAL", "BARYCENTRIC_NORM",
+ "BARYCENTRIC_AVERAGE"]:
q = PIN()
q.polybasis, q.nodes = self.polybasis0, eV
else:
q = PI()
q.npar = self.npar
q.polybasis = self.polybasis0
if self.polydegreetype == "TOTAL":
q.coeffs = degreeTotalToFull(tuple([self.N + 1] * self.npar),
self.npar, eV)
else:
q.coeffs = eV.reshape([self.N + 1] * self.npar)
vbMng(self, "DEL", "Done computing denominator.", 7)
return q
def _setupNumerator(self):
"""Compute rational numerator."""
RROMPyAssert(self._mode, message = "Cannot setup numerator.")
vbMng(self, "INIT", "Starting computation of numerator.", 7)
self._setupInterpolationIndices()
Qevaldiag = polyTimesTable(self.trainedModel.data.Q, self._musUniqueCN,
self._reorder, self._derIdxs,
self.scaleFactorRel)
if self.POD == 1:
Qevaldiag = Qevaldiag.dot(self.samplingEngine.Rscale.T)
elif self.POD == 1/2:
Qevaldiag = Qevaldiag * self.samplingEngine.Rscale
if hasattr(self, "_M_isauto"):
self._setMAuto()
M = self.M
else:
M = reduceDegreeN(self.M, self.S, self.npar, self.polydegreetype)
if M < self.M:
RROMPyWarning(("M too large compared to S. Reducing M by "
"{}").format(self.M - M))
self.M = M
while self.M >= 0:
pParRest = [self.M, self.polybasis, self.verbosity >= 5,
self.polydegreetype == "TOTAL",
{"derIdxs": self._derIdxs, "reorder": self._reorder,
"scl": self.scaleFactorRel}]
if self.polybasis in ppb:
p = PI()
else:
self.computeScaleFactor()
rDWEff = np.array([w * f for w, f in zip(
self.radialDirectionalWeights,
self.scaleFactor)])
pParRest = pParRest[: 2] + [rDWEff] + pParRest[2 :]
pParRest[-1]["optimizeScalingBounds"] = (
self.radialDirectionalWeightsAdapt)
p = RBI()
if self.polybasis in ppb + rbpb:
pParRest += [{"rcond": self.interpRcond}]
wellCond, msg = p.setupByInterpolation(self._musUniqueCN,
Qevaldiag, *pParRest)
vbMng(self, "MAIN", msg, 5)
if wellCond: break
if self.catchInstability:
raise RROMPyException(("Instability in numerator computation: "
"polyfit is poorly conditioned."),
self.catchInstability == 1)
vbMng(self, "MAIN", ("Polyfit is poorly conditioned. Reducing M "
"by 1."), 10)
self.M = self.M - 1
if self.M < 0:
raise RROMPyException(("Instability in computation of numerator. "
"Aborting."))
self.M = M
vbMng(self, "DEL", "Done computing numerator.", 7)
return p
def setupApprox(self) -> int:
"""Compute rational interpolant."""
if self.checkComputedApprox(): return -1
RROMPyAssert(self._mode, message = "Cannot setup approximant.")
vbMng(self, "INIT", "Setting up {}.". format(self.name()), 5)
self.computeSnapshots()
self._setupTrainedModel(self.samplingEngine.projectionMatrix)
self._setupRational(self._setupDenominator())
self.trainedModel.data.approxParameters = copy(self.approxParameters)
vbMng(self, "DEL", "Done setting up approximant.", 5)
return 0
def _setupRational(self, Q:interpEng, P : interpEng = None):
vbMng(self, "INIT", "Starting approximant finalization.", 5)
self.trainedModel.data.Q = Q
if P is None: P = self._setupNumerator()
if self.N > 0 and self.npar == 1:
#check for bad poles
pls = Q.roots()
idxBad = self.HFEngine.flagBadPolesResidues(pls, relative = True)
plsN = self.mapParameterList(self.mapParameterList(self.mu0)(0, 0)
+ self.scaleFactor * pls, "B")(0)
idxBad = np.logical_or(self.HFEngine.flagBadPolesResidues(pls,
relative = True),
self.HFEngine.flagBadPolesResidues(plsN))
if np.any(idxBad):
vbMng(self, "MAIN",
"Removing {} spurious poles out of {} due to poles."\
.format(np.sum(idxBad), self.N), 10)
if isinstance(Q, PIN):
Q.nodes = Q.nodes[np.logical_not(idxBad)]
else:
Q = PI()
Q.npar = self.npar
Q.polybasis = self.polybasis0
Q.coeffs = np.ones(1, dtype = np.complex)
for pl in pls[np.logical_not(idxBad)]:
Q.coeffs = polyTimes(Q.coeffs, [- pl, 1.],
Pbasis = Q.polybasis, Rbasis = Q.polybasis)
Q.coeffs /= np.linalg.norm(Q.coeffs)
self.trainedModel.data.Q = Q
self.N = Q.deg[0]
P = self._setupNumerator()
if (not hasattr(self.HFEngine, "_ignoreResidues")
or not self.HFEngine._ignoreResidues):
#check for bad residues
cfs, pls, _ = rational2heaviside(P, Q)
cfs = cfs[: self.N].T
if self.POD != 1:
cfs = self.samplingEngine.projectionMatrix.dot(cfs)
foci = self.sampler.normalFoci()
ground = self.sampler.groundPotential()
potEff = potential(pls, foci) / ground
potEff[np.logical_or(potEff < 1., np.isinf(pls))] = 1.
cfs[:, np.isinf(pls)] = 0.
cfs /= potEff # rescale by potential
idxBad = self.HFEngine.flagBadPolesResidues(pls, cfs)
if np.any(idxBad):
vbMng(self, "MAIN",
("Removing {} spurious poles out of {} due to "
"residues.").format(np.sum(idxBad), self.N), 10)
if isinstance(Q, PIN):
Q.nodes = Q.nodes[np.logical_not(idxBad)]
else:
Q = PI()
Q.npar = self.npar
Q.polybasis = self.polybasis0
Q.coeffs = np.ones(1, dtype = np.complex)
for pl in pls[np.logical_not(idxBad)]:
Q.coeffs = polyTimes(Q.coeffs, [- pl, 1.],
Pbasis = Q.polybasis, Rbasis = Q.polybasis)
Q.coeffs /= np.linalg.norm(Q.coeffs)
self.trainedModel.data.Q = Q
self.N = Q.deg[0]
P = self._setupNumerator()
self.trainedModel.data.P = P
vbMng(self, "DEL", "Terminated approximant finalization.", 5)
def _computeInterpolantInverseBlocks(self) -> Tuple[List[Np2D], Np2D]:
"""
Compute inverse factors for minimal interpolant target functional.
"""
RROMPyAssert(self._mode, message = "Cannot solve eigenvalue problem.")
self._setupInterpolationIndices()
pvPPar = [self.polybasis0, self._derIdxs, self._reorder,
self.scaleFactorRel]
if hasattr(self, "_M_isauto"): self._setMAuto()
E = max(self.M, self.N)
- while E >= 0:
- if self.polydegreetype == "TOTAL":
- Eeff = E
- idxsB = totalDegreeMaxMask(E, self.npar)
- else: #if self.polydegreetype == "FULL":
- Eeff = [E] * self.npar
- idxsB = fullDegreeMaxMask(E, self.npar)
- TE = pvP(self._musUniqueCN, Eeff, *pvPPar)
- fitOut = pseudoInverse(TE, rcond = self.interpRcond, full = True)
+ fullE = E + 1 == len(self._reorder) == len(self._musUniqueCN)
+ if fullE:
+ mus = self._musUniqueCN[self._reorder]
+ dist = distanceMatrix(mus, magnitude = False)[..., 0]
+ dist[np.arange(E + 1), np.arange(E + 1)] = multifactorial([E])
+ fitinvE = np.prod(dist, axis = 1) ** -1
vbMng(self, "MAIN",
- ("Fitting {} samples with degree {} through {}... "
- "Conditioning of pseudoinverse system: {:.4e}.").format(
- TE.shape[0], E,
- polyfitname(self.polybasis0),
- fitOut[1][1][0] / fitOut[1][1][-1]),
- 5)
- if fitOut[1][0] == TE.shape[1]:
- fitinv = fitOut[0][idxsB, :]
- break
- if self.catchInstability:
- raise RROMPyException(("Instability in denominator "
- "computation: polyfit is poorly "
- "conditioned."),
- self.catchInstability == 1)
- EeqN = E == self.N
- vbMng(self, "MAIN", ("Polyfit is poorly conditioned. Reducing E {}"
- "by 1.").format("and N " * EeqN), 10)
- if EeqN: self.N = self.N - 1
- E -= 1
- if self.N < 0:
- raise RROMPyException(("Instability in computation of "
- "denominator. Aborting."))
- invD = vanderInvTable(fitinv, idxsB, self._reorder, self._derIdxs)
- if self.N == E:
- TN = TE
+ ("Evaluating quasi-Lagrangian basis of degree {} at {} "
+ "sample points.").format(E, E + 1), 5)
+ invD = [np.diag(fitinvE)]
else:
+ while E >= 0:
+ if self.polydegreetype == "TOTAL":
+ Eeff = E
+ idxsB = totalDegreeMaxMask(E, self.npar)
+ else: #if self.polydegreetype == "FULL":
+ Eeff = [E] * self.npar
+ idxsB = fullDegreeMaxMask(E, self.npar)
+ TE = pvP(self._musUniqueCN, Eeff, *pvPPar)
+ fitOut = pseudoInverse(TE, rcond = self.interpRcond,
+ full = True)
+ vbMng(self, "MAIN",
+ ("Fitting {} samples with degree {} through {}... "
+ "Conditioning of pseudoinverse system: {:.4e}.").format(
+ TE.shape[0], E,
+ polyfitname(self.polybasis0),
+ fitOut[1][1][0] / fitOut[1][1][-1]), 5)
+ if fitOut[1][0] == TE.shape[1]:
+ fitinv = fitOut[0][idxsB, :]
+ break
+ if self.catchInstability:
+ raise RROMPyException(("Instability in denominator "
+ "computation: polyfit is poorly "
+ "conditioned."),
+ self.catchInstability == 1)
+ EeqN = E == self.N
+ vbMng(self, "MAIN", ("Polyfit is poorly conditioned. Reducing "
+ "E {} by 1.").format("and N " * EeqN), 10)
+ if EeqN: self.N = self.N - 1
+ E -= 1
+ if self.N < 0:
+ raise RROMPyException(("Instability in computation of "
+ "denominator. Aborting."))
+ invD = vanderInvTable(fitinv, idxsB, self._reorder, self._derIdxs)
+ if self.N == E and not fullE:
+ TN = TE
+ else: #build TN from scratch
if self.polydegreetype == "TOTAL":
Neff = self.N
idxsB = totalDegreeMaxMask(self.N, self.npar)
else: #if self.polydegreetype == "FULL":
Neff = [self.N] * self.npar
idxsB = fullDegreeMaxMask(self.N, self.npar)
TN = pvP(self._musUniqueCN, Neff, *pvPPar)
return invD, TN
def findeveVGExplicit(self, sampleE:sampList, invD:List[Np2D], TN:Np2D,
Rscaling : Np1D = None) -> Tuple[Np1D, Np2D]:
"""
Compute explicitly eigenvalues and eigenvectors of rational denominator
matrix.
"""
RROMPyAssert(self._mode, message = "Cannot solve eigenvalue problem.")
vbMng(self, "INIT", "Building gramian matrix.", 10)
gramian = self.HFEngine.innerProduct(sampleE, sampleE,
is_state = self.approx_state)
if Rscaling is not None:
gramian = (gramian.T * Rscaling.conj()).T * Rscaling
if self.functionalSolve == "NODAL":
SEnd = invD[0].shape[1]
G0 = np.zeros((SEnd,) * 2, dtype = np.complex)
- elif self.functionalSolve == "LOEWNER":
- G0 = gramian
- if self.functionalSolve == "BARYCENTRIC":
+ if self.functionalSolve[:11] == "BARYCENTRIC":
+ G = gramian
nEnd = len(gramian) - 1
else:
nEnd = TN.shape[1]
G = np.zeros((nEnd, nEnd), dtype = np.complex)
for k in range(len(invD)):
iDkN = dot(invD[k], TN)
G += dot(dot(gramian, iDkN).T, iDkN.conj()).T
if self.functionalSolve == "NODAL":
G0 += dot(dot(gramian, invD[k]).T, invD[k].conj()).T
vbMng(self, "DEL", "Done building gramian.", 10)
- if self.functionalSolve == "NORM":
+ if self.functionalSolve in ["NORM", "BARYCENTRIC_NORM"]:
ev, eV = np.linalg.eigh(G)
eV = eV[:, 0]
+ if self.functionalSolve == "BARYCENTRIC_NORM":
+ eV = self.findeVBarycentric(eV)
problem = "eigenproblem"
else:
- if self.functionalSolve == "BARYCENTRIC":
- fitOut = pseudoInverse(gramian, rcond = self.interpRcond,
+ if self.functionalSolve == "BARYCENTRIC_AVERAGE":
+ fitOut = pseudoInverse(G, rcond = self.interpRcond,
full = True)
- barWeigths = fitOut[0].dot(np.ones(nEnd + 1))
- eV = self.findeVBarycentric(barWeigths / np.sum(barWeigths))
+ eV = self.findeVBarycentric(np.sum(fitOut[0], axis = 1))
else:
fitOut = pseudoInverse(G[:-1, :-1], rcond = self.interpRcond,
full = True)
eV = np.append(fitOut[0].dot(G[:-1, -1]), -1.)
ev = fitOut[1][1][::-1]
problem = "linear problem"
vbMng(self, "MAIN",
("Solved {} of size {} with condition number "
"{:.4e}.").format(problem, nEnd, ev[-1] / ev[0]), 5)
- if self.functionalSolve in ["NODAL", "LOEWNER"]:
- q = PI()
- q.npar, q.polybasis, q.coeffs = self.npar, self.polybasis0, eV
- eV, tol, niter, passed = self.findeVNewton(q.roots(), G0)
- if passed:
- vbMng(self, "MAIN",
- ("Newton algorithm for problem of size {} converged "
- "(tol = {:.4e}) in {} iterations.").format(nEnd, tol,
- niter), 5)
- else:
- RROMPyWarning(("Newton algorithm for problem of size {} did "
- "not converge (tol = {:.4e}) after {} "
- "iterations.").format(nEnd, tol, niter))
+ if self.functionalSolve == "NODAL": eV = self.findeVNodal(eV, G0)
return ev, eV
def findeveVGQR(self, RPODE:Np2D, invD:List[Np2D],
TN:Np2D) -> Tuple[Np1D, Np2D]:
"""
Compute eigenvalues and eigenvectors of rational denominator matrix
through SVD of R factor.
"""
RROMPyAssert(self._mode, message = "Cannot solve eigenvalue problem.")
vbMng(self, "INIT", "Building half-gramian matrix stack.", 10)
if self.functionalSolve == "NODAL":
gramian = Np2DLikeGramian(None, dot(RPODE, invD[0]))
- elif self.functionalSolve == "LOEWNER":
- gramian = Np2DLikeGramian(None, RPODE)
- if self.functionalSolve == "BARYCENTRIC":
+ if self.functionalSolve[:11] == "BARYCENTRIC":
+ Rstack = RPODE
nEnd = RPODE.shape[1] - 1
else:
S, nEnd, eWidth = RPODE.shape[0], TN.shape[1], len(invD)
Rstack = np.zeros((S * eWidth, nEnd), dtype = np.complex)
for k in range(eWidth):
Rstack[k * S : (k + 1) * S, :] = dot(RPODE, dot(invD[k], TN))
vbMng(self, "DEL", "Done building half-gramian.", 10)
- if self.functionalSolve in ["NORM", "BARYCENTRIC"]:
- if self.functionalSolve == "NORM":
- _, s, Vh = np.linalg.svd(Rstack, full_matrices = False)
+ if self.functionalSolve in ["NORM", "BARYCENTRIC_NORM",
+ "BARYCENTRIC_AVERAGE"]:
+ _, ev, Vh = np.linalg.svd(Rstack, full_matrices = False)
+ if self.functionalSolve in ["NORM", "BARYCENTRIC_NORM"]:
eV = Vh[-1, :].conj()
- else: #if self.functionalSolve == "BARYCENTRIC":
- _, s, Vh = np.linalg.svd(RPODE, full_matrices = False)
- s[np.logical_not(np.isclose(s, 0.))] **= -2.
- barWeigths = (Vh.T.conj() * s).dot(Vh.dot(np.ones(nEnd + 1)))
- eV = self.findeVBarycentric(barWeigths / np.sum(barWeigths))
- ev = s[::-1]
+ ev = ev[::-1]
+ else: #if self.functionalSolve == "BARYCENTRIC_AVERAGE":
+ ev[np.logical_not(np.isclose(ev, 0.))] **= -2.
+ eV = Vh.T.conj().dot(ev * np.sum(Vh, axis = 1))
+ if self.functionalSolve[:11] == "BARYCENTRIC":
+ eV = self.findeVBarycentric(eV)
problem = "svd problem"
else:
fitOut = pseudoInverse(Rstack[:, :-1], rcond = self.interpRcond,
full = True)
ev = fitOut[1][1][::-1]
eV = np.append(fitOut[0].dot(Rstack[:, -1]), -1.)
problem = "linear problem"
vbMng(self, "MAIN",
("Solved {} of size {} with condition number "
"{:.4e}.").format(problem, nEnd, ev[-1] / ev[0]), 5)
- if self.functionalSolve in ["NODAL", "LOEWNER"]:
- q = PI()
- q.npar, q.polybasis, q.coeffs = self.npar, self.polybasis0, eV
- eV, tol, niter, passed = self.findeVNewton(q.roots(), gramian)
- if passed:
- vbMng(self, "MAIN",
- ("Newton algorithm for problem of size {} converged "
- "(tol = {:.4e}) in {} iterations.").format(nEnd, tol,
- niter), 5)
- else:
- RROMPyWarning(("Newton algorithm for problem of size {} did "
- "not converge (tol = {:.2e}) after {} "
- "iterations.").format(nEnd, tol, niter))
+ if self.functionalSolve == "NODAL": eV = self.findeVNodal(eV, gramian)
return ev, eV
def findeVBarycentric(self, baryWeights:Np1D) -> Np1D:
RROMPyAssert(self._mode,
message = "Cannot solve optimization problem.")
- arrow = np.pad(np.diag(self._musUniqueCN.data[
- self._reorder].flatten()),
+ arrow = np.pad(np.diag(self._musUniqueCN[self._reorder].flatten()),
(1, 0), "constant", constant_values = 1.) + 0.j
arrow[0, 0] = 0.
arrow[0, 1:] = baryWeights
active = np.pad(np.eye(len(baryWeights)), (1, 0), "constant")
eV = eigvals(arrow, active)
return eV[np.logical_not(np.isinf(eV))]
- def findeVNewton(self, nodes0:Np1D, gram:Np2D, maxiter : int = 25,
- tolNewton : float = 1e-10) \
- -> Tuple[Np1D, float, int, bool]:
+ def findeVNodal(self, q0coeffs:Np1D, gram:Np2D, maxiter : int = 25,
+ tolNewton : float = 1e-10) -> Np1D:
RROMPyAssert(self._mode,
message = "Cannot solve optimization problem.")
- algo = self.functionalSolve
- N = len(nodes0)
- nodes = nodes0
+ q = PI()
+ q.npar, q.polybasis, q.coeffs = self.npar, self.polybasis0, q0coeffs
+ nodes = q.roots()
+ N = len(nodes)
grad = np.zeros(N, dtype = np.complex)
hess = np.zeros((N, N), dtype = np.complex)
- mu = np.repeat(self._musUniqueCN.data[self._reorder], N, axis = 1)
for niter in range(maxiter):
- if algo == "NODAL":
- plDist = mu - nodes.reshape(1, -1)
- q0, q = np.prod(plDist, axis = 1), []
- elif algo == "LOEWNER":
- loew = np.pad(np.power(mu - nodes.reshape(1, -1), -1.),
- [(0, 0), (1, 0)], 'constant',
- constant_values = 1.)
- loewI = pseudoInverse(loew)
- Ids = []
+ plDist = distanceMatrix(self._musUniqueCN[self._reorder], nodes,
+ magnitude = False)[:, :, 0]
+ q0, q = np.prod(plDist, axis = 1), []
for jS in range(N):
- if algo == "NODAL":
- mask = np.arange(N) != jS
- q += [np.prod(plDist[:, mask], axis = 1)]
- grad[jS] = q[-1].conj().dot(gram.dot(q0))
- elif algo == "LOEWNER":
- Ids += [loewI.dot(np.power(loew[:, 1 + jS], 2.))]
- zIj, jI = Ids[-1][0], loewI[1 + jS]
- grad[jS] = (zIj * jI).conj().dot(gram.dot(loewI[0]))
+ mask = np.arange(N) != jS
+ q += [np.prod(plDist[:, mask], axis = 1)]
+ grad[jS] = q[-1].conj().dot(gram.dot(q0))
for iS in range(jS + 1):
- if algo == "NODAL":
- if iS == jS:
- hij = 0.
- sij = q[-1].conj().dot(gram.dot(q[-1]))
- else:
- mask = np.logical_and(np.arange(N) != iS,
- np.arange(N) != jS)
- qij = np.prod(plDist[:, mask], axis = 1)
- hij = qij.conj().dot(gram.dot(q0))
- sij = q[-1].conj().dot(gram.dot(q[iS]))
- elif algo == "LOEWNER":
- zIi, iIj = Ids[iS][0], Ids[-1][1 + iS]
- hij = (zIi * iIj * jI).conj().dot(gram.dot(loewI[0]))
- if iS == jS:
- iI = jI
- zIdd = loewI[0].dot(np.power(loew[:, 1 + jS], 3.))
- hij += (zIdd * jI).conj().dot(gram.dot(loewI[0]))
- hij *= 2.
- else:
- jIi, iI = Ids[iS][1 + jS], loewI[1 + iS]
- hij += (zIj * jIi * iI).conj().dot(
- gram.dot(loewI[0]))
- sij = (zIj * jI).conj().dot(gram.dot(zIi * iI))
+ if iS == jS:
+ hij = 0.
+ sij = q[-1].conj().dot(gram.dot(q[-1]))
+ else:
+ mask = np.logical_and(np.arange(N) != iS,
+ np.arange(N) != jS)
+ qij = np.prod(plDist[:, mask], axis = 1)
+ hij = qij.conj().dot(gram.dot(q0))
+ sij = q[-1].conj().dot(gram.dot(q[iS]))
hess[jS, iS] = hij + sij
- if iS < jS: hess[iS, jS] = hij + sij.conj()
+ if iS < jS: hess[iS, jS] = hij + np.conj(sij)
dnodes = np.linalg.solve(hess, grad)
nodes += dnodes
tol = np.linalg.norm(dnodes) / np.linalg.norm(nodes)
if tol < tolNewton: break
- return nodes, tol, niter, niter < maxiter
+ if niter < maxiter:
+ vbMng(self, "MAIN",
+ ("Newton's method for problem of size {} converged "
+ "(err = {:.4e}) in {} iteration{}.").format(
+ N + 1, tol, niter + 1, "s" * (niter > 0)), 5)
+ else:
+ RROMPyWarning(("Newton's method for problem of size {} did "
+ "not converge (err = {:.4e}) after {} "
+ "iterations.").format(N + 1, tol, niter + 1))
+ return nodes
def getResidues(self, *args, **kwargs) -> Np2D:
"""
Obtain approximant residues.
Returns:
Matrix with residues as columns.
"""
return self.trainedModel.getResidues(*args, **kwargs)
-
-#functionalSolve flags:
-#- NORM: linearized target, denominator normalized by norm. Also for multi-D.
-#- DOMINANT: linearized target, denominator normalized by dominant coefficient.
-#- NODAL: nonlinear target, denominator normalized by dominant coefficient, in
-# nodal form.
-#- LOEWNER: nonlinear target, denominator in nodal form.
-#- BARYCENTRIC: by barycentric interpolation via eigensolve of arrowhead
-# matrix, denominator in nodal form.
diff --git a/rrompy/utilities/numerical/__init__.py b/rrompy/utilities/numerical/__init__.py
index b5ea8ad..e362722 100644
--- a/rrompy/utilities/numerical/__init__.py
+++ b/rrompy/utilities/numerical/__init__.py
@@ -1,43 +1,45 @@
# 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 .compress_matrix import compressMatrix
from .halton import haltonGenerate
from .kroneckerer import kroneckerer
from .low_discrepancy import lowDiscrepancy
+from .point_distances import distanceMatrix
from .point_matching import pointMatching, rationalFunctionMatching
from .potential import potential
from .pseudo_inverse import pseudoInverse
from .quadrature_points import quadraturePointsGenerate
from .sobol import sobolGenerate
from .tensor_la import dot, solve
__all__ = [
'compressMatrix',
'haltonGenerate',
'kroneckerer',
'lowDiscrepancy',
+ 'distanceMatrix',
'pointMatching',
'rationalFunctionMatching',
'potential',
'pseudoInverse',
'quadraturePointsGenerate',
'sobolGenerate',
'dot',
'solve'
]
diff --git a/rrompy/utilities/numerical/point_distances.py b/rrompy/utilities/numerical/point_distances.py
new file mode 100644
index 0000000..ac31169
--- /dev/null
+++ b/rrompy/utilities/numerical/point_distances.py
@@ -0,0 +1,77 @@
+# Copyright (C) 2018 by the RROMPy authors
+#
+# This file is part of RROMPy.
+#
+# RROMPy is free software: you can redistribute it and/or modify
+# it under the terms of the GNU Lesser General Public License as published by
+# the Free Software Foundation, either version 3 of the License, or
+# (at your option) any later version.
+#
+# RROMPy is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+# GNU Lesser General Public License for more details.
+#
+# You should have received a copy of the GNU Lesser General Public License
+# along with RROMPy. If not, see .
+#
+
+import numpy as np
+from rrompy.utilities.base.types import Np1D, Np2D, HFEng
+
+__all__ = ['distanceMatrix', 'vectorAngleMatrix', 'chordalMetricMatrix',
+ 'chordalMetricAngleMatrix']
+
+def distanceMatrix(x:Np2D, y : Np2D = None, npar : int = None,
+ magnitude : bool = True, weights : Np1D = None) -> Np2D:
+ if npar is None: npar = x.shape[1] if x.ndim > 1 else 1
+ if y is None: y = x
+ if x.ndim != 3 or x.shape[1] != npar: x = x.reshape(-1, 1, npar)
+ if y.ndim != 2 or y.shape[1] != npar: y = y.reshape(-1, npar)
+ dist = np.repeat(x, len(y), axis = 1) - y
+ if weights is not None: dist *= np.array(weights).flatten()
+ if magnitude:
+ if dist.shape[2] == 1:
+ dist = np.abs(dist)[..., 0]
+ else:
+ dist = np.sum(np.abs(dist) ** 2., axis = 2) ** .5
+ return dist
+
+def vectorAngleMatrix(X:Np2D, Y:Np2D, HFEngine : HFEng = None,
+ is_state : bool = True, radius : float = None) -> Np2D:
+ if HFEngine is None:
+ innerT = np.real(Y.T.conj().dot(X))
+ norm2X = np.sum(np.abs(X) ** 2., axis = 0)
+ norm2Y = np.sum(np.abs(Y) ** 2., axis = 0)
+ else:
+ innerT = np.real(HFEngine.innerProduct(X, Y, is_state = is_state))
+ norm2X = HFEngine.norm(X, is_state = is_state) ** 2.
+ norm2Y = HFEngine.norm(Y, is_state = is_state) ** 2.
+ xInf = np.where(np.isclose(norm2X, 0.))[0]
+ yInf = np.where(np.isclose(norm2Y, 0.))[0]
+ if radius is None: radius = np.mean(norm2Y) ** .5
+ dist2T = (np.tile(norm2Y.reshape(-1, 1), len(norm2X))
+ + norm2X.reshape(1, -1) - 2 * innerT)
+ dist2T[:, xInf], dist2T[yInf, :] = 1., 1.
+ dist2T[np.ix_(yInf, xInf)] = 0.
+ dist2T[dist2T < 0.] = 0.
+ return radius * ((dist2T / (norm2X + radius ** 2.)).T
+ / (norm2Y + radius ** 2.)) ** .5
+
+def chordalMetricMatrix(x:Np1D, y:Np1D, radius : float = 1.) -> Np2D:
+ x, y = np.array(x), np.array(y)
+ xInf, yInf = np.where(np.isinf(x))[0], np.where(np.isinf(y))[0]
+ x[xInf], y[yInf] = 0., 0.
+ distT = distanceMatrix(x, y)
+ distT[:, xInf], distT[yInf, :] = 1., 1.
+ distT[np.ix_(yInf, xInf)] = 0.
+ return radius * ((distT / (np.abs(x) ** 2. + radius ** 2.) ** .5).T
+ / (np.abs(y) ** 2. + radius ** 2.) ** .5)
+
+def chordalMetricAngleMatrix(x:Np1D, y:Np1D, w : float = 0, X : Np2D = None,
+ Y : Np2D = None, HFEngine : HFEng = None,
+ is_state : bool = True) -> Np2D:
+ dist = chordalMetricMatrix(x, y)
+ if w == 0: return dist
+ distAdj = vectorAngleMatrix(X, Y, HFEngine, is_state)
+ return (dist + w * distAdj) / (1. + w)
diff --git a/rrompy/utilities/numerical/point_matching.py b/rrompy/utilities/numerical/point_matching.py
index 741e6be..74082c0 100644
--- a/rrompy/utilities/numerical/point_matching.py
+++ b/rrompy/utilities/numerical/point_matching.py
@@ -1,127 +1,87 @@
# 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.optimize import linear_sum_assignment as LSA
+from .point_distances import distanceMatrix, chordalMetricAngleMatrix
from rrompy.utilities.base.types import Tuple, List, ListAny, Np1D, Np2D, HFEng
from rrompy.utilities.exception_manager import RROMPyAssert
-__all__ = ['pointMatching', 'rationalFunctionMatching', 'angleTable',
- 'chordalMetricTable', 'chordalMetricAdjusted']
+__all__ = ['pointMatching', 'rationalFunctionMatching']
-def pointMatching(distanceMatrix:Np2D) -> Tuple[Np1D, Np1D]:
- return LSA(distanceMatrix)
+def pointMatching(distMatrix:Np2D) -> Tuple[Np1D, Np1D]:
+ return LSA(distMatrix)
def rationalFunctionMatching(poles:List[Np1D], coeffs:List[Np2D],
featPts:Np2D, matchingWeight:float, supps:ListAny,
projMat:Np2D, HFEngine : HFEng = None,
is_state : bool = True, root : int = None) \
-> Tuple[List[Np1D], List[Np2D]]:
"""
Match poles and residues of a set of rational functions.
Args:
poles: List of (lists of) poles.
coeffs: List of (lists of) residues.
featPts: Marginal parameters corresponding to rational models.
matchingWeight: Matching weight in distance computation.
supps: Support indices for projection matrix.
projMat: Projection matrix for residues.
HFEngine(optional): Engine for distance evaluation. Defaults to None,
i.e. Euclidean metric.
is_state(optional): Whether residues are of system state. Defaults to
True.
root(optional): Root of search tree. Defaults to None, i.e.
automatically chosen.
Returns:
Matched list of (lists of) poles and list of (lists of) residues.
"""
M, N = len(featPts), len(poles[0])
RROMPyAssert(len(poles), M, "Number of rational functions to be matched")
RROMPyAssert(len(coeffs), M, "Number of rational functions to be matched")
if M <= 1: return poles, coeffs
- featDist = np.sum(np.abs(np.repeat(featPts, M, 0)
- - np.tile(featPts, [M, 1])), axis = 1)
+ featDist = distanceMatrix(featPts).flatten()
free = list(range(M))
if root is None: #start from sample points closest to each other
root = np.argpartition(featDist, M)[M] % M
fixed = [free.pop(root)]
featDist = featDist.reshape(M, M)
for j in range(M - 1, 0, -1):
#find closest point
idx = np.argmin(featDist[np.ix_(fixed, free)].flatten())
Ifix = fixed[idx // j]
fixed += [free.pop(idx % j)]
Ifree = fixed[-1]
plsfix, plsfree = poles[Ifix], poles[Ifree]
resfix, resfree = None, None
if matchingWeight != 0:
resfix, resfree = coeffs[Ifix][: N].T, coeffs[Ifree][: N].T
if isinstance(projMat, (np.ndarray,)):
suppfix, suppfree = supps[Ifix], supps[Ifree]
resfix = projMat[:, suppfix : suppfix + len(resfix)].dot(
resfix)
resfree = projMat[:, suppfree : suppfree + len(resfree)].dot(
resfree)
#build assignment distance matrix
- distj = chordalMetricAdjusted(plsfix, plsfree, matchingWeight, resfix,
- resfree, HFEngine, is_state)
+ distj = chordalMetricAngleMatrix(plsfix, plsfree, matchingWeight,
+ resfix, resfree, HFEngine, is_state)
reordering = pointMatching(distj)[1]
poles[Ifree] = poles[Ifree][reordering]
coeffs[Ifree][: N] = coeffs[Ifree][reordering]
return poles, coeffs
-
-def angleTable(X:Np2D, Y:Np2D, HFEngine : HFEng = None,
- is_state : bool = True, radius : float = None) -> Np2D:
- if HFEngine is None:
- innerT = np.real(Y.T.conj().dot(X))
- norm2X = np.sum(np.abs(X) ** 2., axis = 0)
- norm2Y = np.sum(np.abs(Y) ** 2., axis = 0)
- else:
- innerT = np.real(HFEngine.innerProduct(X, Y, is_state = is_state))
- norm2X = HFEngine.norm(X, is_state = is_state) ** 2.
- norm2Y = HFEngine.norm(Y, is_state = is_state) ** 2.
- xInf = np.where(np.isclose(norm2X, 0.))[0]
- yInf = np.where(np.isclose(norm2Y, 0.))[0]
- if radius is None: radius = np.mean(norm2Y) ** .5
- dist2T = (np.tile(norm2Y.reshape(-1, 1), len(norm2X))
- + norm2X.reshape(1, -1) - 2 * innerT)
- dist2T[:, xInf], dist2T[yInf, :] = 1., 1.
- dist2T[np.ix_(yInf, xInf)] = 0.
- dist2T[dist2T < 0.] = 0.
- return radius * ((dist2T / (norm2X + radius ** 2.)).T
- / (norm2Y + radius ** 2.)) ** .5
-
-def chordalMetricTable(x:Np1D, y:Np1D, radius : float = 1.) -> Np2D:
- x, y = np.array(x), np.array(y)
- xInf, yInf = np.where(np.isinf(x))[0], np.where(np.isinf(y))[0]
- x[xInf], y[yInf] = 0., 0.
- distT = np.abs(np.tile(y.reshape(-1, 1), len(x)) - x.reshape(1, -1))
- distT[:, xInf], distT[yInf, :] = 1., 1.
- distT[np.ix_(yInf, xInf)] = 0.
- return radius * ((distT / (np.abs(x) ** 2. + radius ** 2.) ** .5).T
- / (np.abs(y) ** 2. + radius ** 2.) ** .5)
-
-def chordalMetricAdjusted(x:Np1D, y:Np1D, w : float = 0, X : Np2D = None,
- Y : Np2D = None, HFEngine : HFEng = None,
- is_state : bool = True) -> Np2D:
- dist = chordalMetricTable(x, y)
- if w == 0: return dist
- distAdj = angleTable(X, Y, HFEngine, is_state)
- return (dist + w * distAdj) / (1. + w)
diff --git a/rrompy/utilities/poly_fitting/nearest_neighbor/val.py b/rrompy/utilities/poly_fitting/nearest_neighbor/val.py
index d23f3ff..4bd19ec 100644
--- a/rrompy/utilities/poly_fitting/nearest_neighbor/val.py
+++ b/rrompy/utilities/poly_fitting/nearest_neighbor/val.py
@@ -1,41 +1,39 @@
# 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.numerical import distanceMatrix
from rrompy.utilities.base.types import Np1D, Np2D, paramList
from rrompy.parameter import checkParameterList
__all__ = ['polyval']
def polyval(x:paramList, cL:Np2D, supportPoints:paramList,
nNeighbors : int = 1, directionalWeights : Np1D = None) -> Np2D:
supportPoints = checkParameterList(supportPoints, return_data = True)
if directionalWeights is None:
directionalWeights = np.ones(supportPoints.shape[1])
npar = supportPoints.shape[1]
x = checkParameterList(x, npar, return_data = True)
- muDiff = (np.repeat(supportPoints.reshape((len(supportPoints), 1,
- npar)), len(x), axis = 1) - x
- ) * directionalWeights
- dist = (np.sum(np.abs(muDiff) ** 2., axis = 2)
- + np.finfo(float).eps ** 2.) ** -.5
+ muDiff = distanceMatrix(supportPoints, x, weights = directionalWeights)
+ dist = (muDiff ** 2. + np.finfo(float).eps ** 2.) ** -.5
if len(dist) > nNeighbors:
iOut = np.argpartition(dist, - nNeighbors, axis = 0)[: - nNeighbors]
np.put_along_axis(dist, iOut, 0., 0)
dist /= np.linalg.norm(dist, axis = 0, ord = 1)
return np.moveaxis(np.tensordot(dist.T, cL, 1), 0, -1)
diff --git a/rrompy/utilities/poly_fitting/polynomial/polynomial_interpolator.py b/rrompy/utilities/poly_fitting/polynomial/polynomial_interpolator.py
index c5079f9..d434817 100644
--- a/rrompy/utilities/poly_fitting/polynomial/polynomial_interpolator.py
+++ b/rrompy/utilities/poly_fitting/polynomial/polynomial_interpolator.py
@@ -1,227 +1,226 @@
# 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 collections.abc import Iterable
from itertools import combinations
from rrompy.utilities.base.types import (List, ListAny, DictAny, Np1D, Np2D,
paramList)
from rrompy.utilities.base import freepar as fp
from rrompy.utilities.poly_fitting.interpolator import GenericInterpolator
from rrompy.utilities.poly_fitting.custom_fit import customFit
from .base import polyfitname
from .val import polyval
from .roots import polyroots
from .vander import polyvander as pv
from .polynomial_algebra import changePolyBasis, polyTimes
-from rrompy.utilities.numerical import dot
+from rrompy.utilities.numerical import dot, distanceMatrix
from rrompy.utilities.numerical.degree import degreeTotalToFull
from rrompy.utilities.exception_manager import RROMPyAssert, RROMPyException
from rrompy.parameter import checkParameterList
__all__ = ['PolynomialInterpolator', 'PolynomialInterpolatorNodal']
class PolynomialInterpolator(GenericInterpolator):
"""Function class with setup by polynomial interpolation."""
def __init__(self, other = None):
if other is None: return
self.coeffs = other.coeffs
self.npar = other.npar
self.polybasis = other.polybasis
@property
def shape(self):
if self.coeffs.ndim > self.npar:
sh = self.coeffs.shape[self.npar :]
else: sh = tuple([1])
return sh
@property
def deg(self):
return [x - 1 for x in self.coeffs.shape[: self.npar]]
def __getitem__(self, key):
return self.coeffs[key]
def __call__(self, mu:paramList, der : List[int] = None,
scl : Np1D = None):
if hasattr(self, "_dirPivot"):
mu = checkParameterList(mu)(self._dirPivot)
return polyval(mu, self.coeffs, self.polybasis, der, scl)
def __copy__(self):
return PolynomialInterpolator(self)
def __deepcopy__(self, memo):
other = PolynomialInterpolator()
other.coeffs, other.npar, other.polybasis = copy(
(self.coeffs, self.npar, self.polybasis), memo)
return other
def pad(self, nleft : List[int] = None, nright : List[int] = None):
if nleft is None: nleft = [0] * len(self.shape)
if nright is None: nright = [0] * len(self.shape)
if not isinstance(nleft, Iterable): nleft = [nleft]
if not isinstance(nright, Iterable): nright = [nright]
RROMPyAssert(len(self.shape), len(nleft), "Shape of output")
RROMPyAssert(len(self.shape), len(nright), "Shape of output")
padwidth = [(0, 0)] * self.npar
padwidth = padwidth + [(l, r) for l, r in zip(nleft, nright)]
self.coeffs = np.pad(self.coeffs, padwidth, "constant",
constant_values = (0., 0.))
def postmultiplyTensorize(self, A:Np2D):
RROMPyAssert(A.shape[0], self.shape[-1], "Shape of output")
self.coeffs = dot(self.coeffs, A)
def setupByInterpolation(self, support:paramList, values:ListAny,
deg:int, polybasis : str = "MONOMIAL",
verbose : bool = True, totalDegree : bool = True,
vanderCoeffs : DictAny = {},
fitCoeffs : DictAny = {}):
support = checkParameterList(support)
self.npar = support.shape[1]
self.polybasis = polybasis
if not totalDegree and not isinstance(deg, Iterable):
deg = [deg] * self.npar
vander = pv(support, deg, basis = polybasis, **vanderCoeffs)
RROMPyAssert(len(vander), len(values), "Number of support values")
outDim = values.shape[1:]
values = values.reshape(values.shape[0], -1)
fitOut = customFit(vander, values, full = True, **fitCoeffs)
P = fitOut[0]
if verbose:
msg = ("Fitting {} samples with degree {} through {}... "
"Conditioning of LS system: {:.4e}.").format(
len(vander), deg, polyfitname(self.polybasis),
fitOut[1][2][0] / fitOut[1][2][-1])
else: msg = None
if totalDegree:
self.coeffs = degreeTotalToFull(tuple([deg + 1] * self.npar)
+ outDim, self.npar, P)
else:
self.coeffs = P.reshape(tuple([d + 1 for d in deg]) + outDim)
return fitOut[1][1] == vander.shape[1], msg
def roots(self, marginalVals : ListAny = [fp]):
RROMPyAssert(self.shape, (1,), "Shape of output")
RROMPyAssert(len(marginalVals), self.npar, "Number of parameters")
rDim = marginalVals.index(fp)
if rDim < len(marginalVals) - 1 and fp in marginalVals[rDim + 1 :]:
raise RROMPyException(("Exactly 1 'freepar' entry in "
"marginalVals must be provided."))
return polyroots(self.coeffs, self.polybasis, marginalVals)
class PolynomialInterpolatorNodal(PolynomialInterpolator):
"""
Function class with setup by polynomial interpolation. Stores roots of
monomial polynomial instead of coefficients. Only for 1D.
"""
def __init__(self, other = None):
self.npar = 1
if other is None: return
self.nodes = other.nodes
self.polybasis = other.polybasis
@property
def nodes(self):
return self._nodes
@nodes.setter
def nodes(self, nodes):
self.coeffs = None
self._nodes = nodes
@property
def coeffs(self):
if self._coeffs is None: self.buildCoeffs()
return self._coeffs
@coeffs.setter
def coeffs(self, coeffs):
self._coeffs = coeffs
@property
def shape(self):
return (1,)
@property
def deg(self):
return [len(self.nodes)]
def __getitem__(self, key):
return self.coeffs[key]
def __call__(self, mu:paramList, der : List[int] = None,
scl : Np1D = None):
dirPivot = self._dirPivot if hasattr(self, "_dirPivot") else 0
if der is None: der = 0
elif isinstance(der, (list,tuple,np.ndarray,)): der = der[dirPivot]
if scl is None: scl = 1.
elif isinstance(scl, (list,tuple,np.ndarray,)): scl = scl[dirPivot]
mu = checkParameterList(mu)(dirPivot)
val = np.zeros(len(mu), dtype = np.complex)
if der == self.deg[0]:
val[:] = 1.
elif der >= 0 and der < self.deg[0]:
- plDist = (np.repeat(np.expand_dims(mu, 1), self.deg[0], axis = 1)
- - self.nodes.reshape(1, -1))
+ plDist = distanceMatrix(mu, self.nodes, magnitude = False)[:, :, 0]
for terms in combinations(np.arange(self.deg[0]),
self.deg[0] - der):
val += np.prod(plDist[:, list(terms)], axis = 1)
return scl ** der * fact(der) * val
def __copy__(self):
return PolynomialInterpolatorNodal(self)
def __deepcopy__(self, memo):
other = PolynomialInterpolatorNodal()
other.nodes, other.polybasis = copy((self.nodes, self.polybasis), memo)
return other
def buildCoeffs(self):
local = [np.array([- pl, 1.], dtype = np.complex) for pl in self.nodes]
N = len(local)
while N > 1:
for j in range(N // 2):
local[j] = polyTimes(local[j], local[- 1 - j])
local = local[(N - 1) // 2 :: -1]
N = len(local)
self._coeffs = changePolyBasis(local[0], None, "MONOMIAL",
self.polybasis)
def pad(self, *args, **kwargs):
raise RROMPyException(("Padding not allowed for polynomials in nodal "
"form"))
def postmultiplyTensorize(self, *args, **kwargs):
raise RROMPyException(("Post-multiply not allowed for polynomials in "
"nodal form"))
def setupByInterpolation(self, support:paramList, *args, **kwargs):
support = checkParameterList(support)
self.npar = support.shape[1]
if self.npar > 1:
raise RROMPyException(("Polynomial in nodal form must have "
"scalar output"))
output = super().setupByInterpolation(support, *args, **kwargs)
self._nodes = super().roots()
return output
def roots(self, marginalVals : ListAny = [fp]):
return self.nodes
diff --git a/tests/3_reduction_methods_1D/rational_interpolant_1d.py b/tests/3_reduction_methods_1D/rational_interpolant_1d.py
index 9f73335..50c9889 100644
--- a/tests/3_reduction_methods_1D/rational_interpolant_1d.py
+++ b/tests/3_reduction_methods_1D/rational_interpolant_1d.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 numpy as np
from matrix_fft import matrixFFT
from rrompy.reduction_methods import RationalInterpolant as RI
from rrompy.parameter.parameter_sampling import (QuadratureSampler as QS,
ManualSampler as MS)
from rrompy.parameter import checkParameterList
def test_monomials(capsys):
mu = 1.5
solver = matrixFFT()
params = {"POD": False, "S": 10, "robustTol": 1e-6, "interpRcond": 1e-3,
"polybasis": "MONOMIAL", "sampler": QS([1.5, 6.5], "UNIFORM")}
approx = RI(solver, 4., approx_state = True, approxParameters = params,
verbosity = 10)
approx.setupApprox()
out, err = capsys.readouterr()
- assert "poorly conditioned. Reducing E " in out
+ assert "below tolerance. Reducing N " in out
+ assert "poorly conditioned. Reducing M " in out
assert len(err) == 0
assert np.isclose(approx.normErr(mu)[0], 1.4746e-05, atol = 1e-4)
def test_well_cond():
mu = 1.5
solver = matrixFFT()
params = {"POD": True, "S": 10, "robustTol": 1e-14, "interpRcond": 1e-10,
"polybasis": "CHEBYSHEV", "sampler": QS([1., 7.], "CHEBYSHEV")}
approx = RI(solver, 4., approx_state = True, approxParameters = params,
verbosity = 0)
approx.setupApprox()
poles = approx.getPoles()
for lambda_ in np.arange(1, 8):
assert np.isclose(np.min(np.abs(poles - lambda_)), 0., atol = 1e-4)
for mu in approx.mus:
assert np.isclose(approx.normErr(mu)[0], 0., atol = 1e-8)
def test_hermite():
mu = 1.5
solver = matrixFFT()
sampler0 = QS([1., 7.], "CHEBYSHEV")
points = checkParameterList(np.tile(sampler0.generatePoints(4)(0), 3))
params = {"POD": True, "S": 12, "polybasis": "CHEBYSHEV",
"sampler": MS([1., 7.], points = points)}
approx = RI(solver, 4., approx_state = True, approxParameters = params,
verbosity = 0)
approx.setupApprox()
poles = approx.getPoles()
for lambda_ in np.arange(1, 8):
assert np.isclose(np.min(np.abs(poles - lambda_)), 0., atol = 1e-4)
for mu in approx.mus:
assert np.isclose(approx.normErr(mu)[0], 0., atol = 1e-8)