diff --git a/rrompy/reduction_methods/pivoted/greedy/generic_pivoted_greedy_approximant.py b/rrompy/reduction_methods/pivoted/greedy/generic_pivoted_greedy_approximant.py
index cf15d14..aad2f7f 100644
--- a/rrompy/reduction_methods/pivoted/greedy/generic_pivoted_greedy_approximant.py
+++ b/rrompy/reduction_methods/pivoted/greedy/generic_pivoted_greedy_approximant.py
@@ -1,820 +1,822 @@
# Copyright (C) 2018 by the RROMPy authors
#
# This file is part of RROMPy.
#
# RROMPy is free software: you can redistribute it and/or modify
# it under the terms of the GNU Lesser General Public License as published by
# the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
#
# RROMPy is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU Lesser General Public License for more details.
#
# You should have received a copy of the GNU Lesser General Public License
# along with RROMPy. If not, see .
#
from abc import abstractmethod
from copy import deepcopy as copy
import numpy as np
from matplotlib import pyplot as plt
from rrompy.reduction_methods.pivoted.generic_pivoted_approximant import (
GenericPivotedApproximantBase,
GenericPivotedApproximantNoMatch,
GenericPivotedApproximant)
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 import dot
from rrompy.utilities.numerical.point_matching import (pointMatching,
chordalMetricAdjusted, potential)
from rrompy.utilities.exception_manager import (RROMPyException, RROMPyAssert,
RROMPyWarning)
from rrompy.parameter import checkParameterList, emptyParameterList
__all__ = ['GenericPivotedGreedyApproximantNoMatch',
'GenericPivotedGreedyApproximant']
class GenericPivotedGreedyApproximantBase(GenericPivotedApproximantBase):
_allowedEstimatorKindsMarginal = ["LEAVE_ONE_OUT", "LOOK_AHEAD",
"LOOK_AHEAD_RECOVER", "NONE"]
def __init__(self, *args, **kwargs):
self._preInit()
self._addParametersToList(["matchingWeightError",
"cutOffToleranceError",
"errorEstimatorKindMarginal",
"greedyTolMarginal", "maxIterMarginal"],
[0., "AUTO", "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 hasattr(self, "_scaleFactorDer"):
scaleFactorDerold = self.scaleFactorDer
else: scaleFactorDerold = -1
if isinstance(scaleFactorDer, (str,)):
scaleFactorDer = scaleFactorDer.upper()
self._scaleFactorDer = scaleFactorDer
self._approxParameters["scaleFactorDer"] = self._scaleFactorDer
if scaleFactorDerold != self._scaleFactorDer: self.resetSamples()
@property
def 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 cutOffToleranceError(self):
"""Value of cutOffToleranceError."""
return self._cutOffToleranceError
@cutOffToleranceError.setter
def cutOffToleranceError(self, cutOffToleranceError):
if isinstance(cutOffToleranceError, (str,)):
cutOffToleranceError = cutOffToleranceError.upper()\
.strip().replace(" ","")
if cutOffToleranceError != "AUTO":
RROMPyWarning(("String value of cutOffToleranceError not "
"recognized. Overriding to 'AUTO'."))
cutOffToleranceError == "AUTO"
self._cutOffToleranceError = cutOffToleranceError
self._approxParameters["cutOffToleranceError"] = (
self.cutOffToleranceError)
@property
def greedyTolMarginal(self):
"""Value of greedyTolMarginal."""
return self._greedyTolMarginal
@greedyTolMarginal.setter
def greedyTolMarginal(self, greedyTolMarginal):
if greedyTolMarginal < 0:
raise RROMPyException("greedyTolMarginal must be non-negative.")
if (hasattr(self, "_greedyTolMarginal")
and self.greedyTolMarginal is not None):
greedyTolMarginalold = self.greedyTolMarginal
else:
greedyTolMarginalold = -1
self._greedyTolMarginal = greedyTolMarginal
self._approxParameters["greedyTolMarginal"] = self.greedyTolMarginal
if greedyTolMarginalold != self.greedyTolMarginal:
self.resetSamples()
@property
def maxIterMarginal(self):
"""Value of maxIterMarginal."""
return self._maxIterMarginal
@maxIterMarginal.setter
def maxIterMarginal(self, maxIterMarginal):
if maxIterMarginal <= 0:
raise RROMPyException("maxIterMarginal must be positive.")
if (hasattr(self, "_maxIterMarginal")
and self.maxIterMarginal is not None):
maxIterMarginalold = self.maxIterMarginal
else:
maxIterMarginalold = -1
self._maxIterMarginal = maxIterMarginal
self._approxParameters["maxIterMarginal"] = self.maxIterMarginal
if maxIterMarginalold != self.maxIterMarginal:
self.resetSamples()
def resetSamples(self):
"""Reset samples."""
super().resetSamples()
if not hasattr(self, "_temporaryPivot"):
self._mus = emptyParameterList()
self._musMarginal = emptyParameterList()
if hasattr(self, "samplerMarginal"): self.samplerMarginal.reset()
if hasattr(self, "samplingEngine") and self.samplingEngine is not None:
self.samplingEngine.resetHistory()
def _getPolesResExact(self, HITest, foci:Tuple[float, float],
ground:float) -> Tuple[Np1D, Np2D]:
if self.cutOffToleranceError == "AUTO":
cutOffTolErr = self.cutOffTolerance
else:
cutOffTolErr = self.cutOffToleranceError
polesEx = copy(HITest.poles)
idxExEff = np.where(potential(polesEx, foci) - ground
<= cutOffTolErr * ground)[0]
if self.matchingWeightError != 0:
resEx = HITest.coeffs[idxExEff]
else:
resEx = None
return polesEx[idxExEff], resEx
def _getDistanceApp(self, polesEx:Np1D, resEx:Np2D, muTest:paramVal,
foci:Tuple[float, float], ground:float) -> float:
if self.cutOffToleranceError == "AUTO":
cutOffTolErr = self.cutOffTolerance
else:
cutOffTolErr = self.cutOffToleranceError
polesAp = self.trainedModel.interpolateMarginalPoles(muTest)[..., 0]
idxApEff = np.where(potential(polesAp, foci) - ground
<= cutOffTolErr * ground)[0]
polesAp = polesAp[idxApEff]
if self.matchingWeightError != 0:
resAp = self.trainedModel.interpolateMarginalCoeffs(muTest)[
idxApEff, :, 0]
resEx = self.trainedModel.data.projMat[:,
: resEx.shape[1]].dot(resEx.T)
resAp = self.trainedModel.data.projMat[:,
: resAp.shape[1]].dot(resAp.T)
else:
resAp = None
dist = chordalMetricAdjusted(polesEx, polesAp,
self.matchingWeightError, resEx, resAp,
self.HFEngine, False)
pmR, pmC = pointMatching(dist)
return np.mean(dist[pmR, pmC])
def getErrorEstimatorMarginalLeaveOneOut(self) -> Np1D:
err = np.zeros(len(self.trainedModel.data.musMarginal))
self._musMarginalTestIdxs = np.arange(len(err))
if len(err) <= 1:
err[:] = np.inf
return err
_tMdataFull = copy(self.trainedModel.data)
_musMExcl = None
self.verbosity -= 35
self.trainedModel.verbosity -= 35
foci = self.samplerPivot.normalFoci()
ground = self.samplerPivot.groundPotential()
for j in range(len(err)):
jEff = j - (j > 0)
muTest = self.trainedModel.data.musMarginal[jEff]
polesEx, resEx = self._getPolesResExact(
self.trainedModel.data.HIs[jEff],
foci, ground)
if j > 0: self.musMarginal.insert(_musMExcl, j - 1)
_musMExcl = self.musMarginal[j]
self.musMarginal.pop(j)
if len(polesEx) == 0: continue
self._updateTrainedModelMarginalSamples([j])
self._finalizeMarginalization()
err[j] = self._getDistanceApp(polesEx, resEx, muTest, foci, ground)
self._updateTrainedModelMarginalSamples()
self.musMarginal.append(_musMExcl)
self.verbosity += 35
self.trainedModel.verbosity += 35
self.trainedModel.data = _tMdataFull
return err
def getErrorEstimatorMarginalLookAhead(self) -> Np1D:
if not hasattr(self.trainedModel, "_musMExcl"):
err = np.zeros(0)
err[:] = np.inf
self._musMarginalTestIdxs = np.zeros(0, dtype = int)
return err
err = np.zeros(len(self.trainedModel._musMExcl))
self._musMarginalTestIdxs = np.array(self.trainedModel._idxExcl,
dtype = int)
self.verbosity -= 35
self.trainedModel.verbosity -= 35
foci = self.samplerPivot.normalFoci()
ground = self.samplerPivot.groundPotential()
for j, (muTest, HITest) in enumerate(zip(self.trainedModel._musMExcl,
self.trainedModel._HIsExcl)):
polesEx, resEx = self._getPolesResExact(HITest, foci, ground)
if len(polesEx) == 0: continue
err[j] = self._getDistanceApp(polesEx, resEx, muTest, foci, ground)
self.verbosity += 35
self.trainedModel.verbosity += 35
return err
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 == "LEAVE_ONE_OUT":
err = self.getErrorEstimatorMarginalLeaveOneOut()
elif self.errorEstimatorKindMarginal[: 10] == "LOOK_AHEAD":
err = self.getErrorEstimatorMarginalLookAhead()
else:#if self.errorEstimatorKindMarginal == "NONE":
err = self.getErrorEstimatorMarginalNone()
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))):
fig = plt.figure(figsize = plt.figaspect(1. / self.nparMarginal))
for jpar in range(self.nparMarginal):
ax = fig.add_subplot(1, self.nparMarginal, 1 + jpar)
if self.errorEstimatorKindMarginal == "LEAVE_ONE_OUT":
musre = copy(self.trainedModel.data.musMarginal.re.data)
else:#if self.errorEstimatorKindMarginal[: 10] == "LOOK_AHEAD":
if not hasattr(self.trainedModel, "_musMExcl"): return
musre = np.real(self.trainedModel._musMExcl)
if len(idxMax) > 0 and estMax is not None:
maxrej = musre[idxMax, jpar]
errCP = copy(est)
idx = np.delete(np.arange(self.nparMarginal), jpar)
while len(musre) > 0:
if self.nparMarginal == 1:
currIdx = np.arange(len(musre))
else:
currIdx = np.where(np.isclose(np.sum(
np.abs(musre[:, idx] - musre[0, idx]), 1), 0.))[0]
currIdxSorted = currIdx[np.argsort(musre[currIdx, jpar])]
ax.semilogy(musre[currIdxSorted, jpar],
errCP[currIdxSorted], 'k.-', linewidth = 1)
musre = np.delete(musre, currIdx, 0)
errCP = np.delete(errCP, currIdx)
ax.semilogy(self.musMarginal.re(jpar),
(self.greedyTolMarginal,) * len(self.musMarginal),
'*m')
if len(idxMax) > 0 and estMax is not None:
ax.semilogy(maxrej, estMax, 'xr')
ax.grid()
plt.tight_layout()
plt.show()
def _addMarginalSample(self, mus:paramList):
mus = checkParameterList(mus, self.nparMarginal)[0]
if len(mus) == 0: return
- nmus = len(mus)
+ nmusOld, nmus = len(self.musMarginal), len(mus)
+ if (hasattr(self, "trainedModel") and self.trainedModel is not None
+ and hasattr(self.trainedModel, "_musMExcl")):
+ nmusOld += len(self.trainedModel._musMExcl)
vbMng(self, "MAIN",
("Adding marginal sample point{} no. {}{} at {} to training "
- "set.").format("s" * (nmus > 1), len(self.musMarginal) + 1,
- "--{}".format(len(self.musMarginal) + nmus) * (nmus > 1),
- mus), 3)
+ "set.").format("s" * (nmus > 1), nmusOld + 1,
+ "--{}".format(nmusOld + nmus) * (nmus > 1), mus),
+ 3)
self.musMarginal.append(mus)
self.setupApproxPivoted(mus)
self._poleMatching()
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.")
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)
self._addMarginalSample(self.samplerMarginal.points[idxAdded])
errorEstTest, muidx, maxErrorEst = self.errorEstimatorMarginal(True)
if plotEst == "ALL":
self.plotEstimatorMarginal(errorEstTest, muidx, maxErrorEst)
return (errorEstTest, self._musMarginalTestIdxs[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 _finalizeSnapshots(self):
self.samplingEngine = self._samplingEngineOld
for muM, sEN in zip(self.musMargLoc, self.samplingEngs):
self.samplingEngine.samples += [sEN.samples]
self.samplingEngine.nsamples += [sEN.nsamples]
self.samplingEngine.mus += [sEN.mus]
self.samplingEngine.musMarginal.append(muM)
self.samplingEngine._derIdxs += [[(0,) * self.npar]
for _ in range(sEN.nsamples)]
if self.POD:
self.samplingEngine.RPOD += [sEN.RPOD]
self.samplingEngine.samples_full += [copy(sEN.samples_full)]
- self.samplingEngine.coalesceSamples()
def _preSetupApproxPivoted(self, mus:paramList) -> Tuple[ListAny, ListAny]:
self.computeScaleFactor()
if self.trainedModel is None:
self.trainedModel = self.tModelType()
self.trainedModel.verbosity = self.verbosity
self.trainedModel.timestamp = self.timestamp
datadict = {"mu0": self.mu0, "mus": None,
"projMat": np.zeros((0, 0)),
"scaleFactor": self.scaleFactor,
"rescalingExp": self.HFEngine.rescalingExp,
"directionPivot": self.directionPivot}
self.trainedModel.data = self.initializeModelData(datadict)[0]
self.trainedModel.data.Qs, self.trainedModel.data.Ps = [], []
self._trainedModelOld = copy(self.trainedModel)
self._scaleFactorOldPivot = copy(self.scaleFactor)
self.scaleFactor = self.scaleFactorPivot
self._temporaryPivot = 1
self._samplingEngineOld = copy(self.samplingEngine)
self.musMargLoc, self.samplingEngs = [], [None] * len(mus)
Qs, Ps = [None] * len(mus), [None] * len(mus)
self.verbosity -= 15
return Qs, Ps
def _postSetupApproxPivoted(self, mus:paramList, Qs:ListAny, Ps:ListAny):
self.scaleFactor = self._scaleFactorOldPivot
del self._scaleFactorOldPivot, self._temporaryPivot
self._finalizeSnapshots()
del self._samplingEngineOld, self.musMargLoc, self.samplingEngs
self._mus = self.samplingEngine.musCoalesced
self.trainedModel = self._trainedModelOld
del self._trainedModelOld
padLeft = self.trainedModel.data.projMat.shape[1]
self.trainedModel.data.mus = copy(self.mus)
self.trainedModel.data.musMarginal = copy(self.musMarginal)
- padRight = self.samplingEngine.nsamplesTot - padLeft
+ padRight = self.samplingEngine.nsamplesCoalesced - padLeft
nmusMOld = len(self.trainedModel.data.Ps)
for j in range(nmusMOld):
self.trainedModel.data.Ps[j].pad(0, padRight)
self.trainedModel.data.HIs[j].pad(0, padRight)
if hasattr(self.trainedModel, "_PsExcl"):
nmusMOldExcl = len(self.trainedModel._PsExcl)
for j in range(nmusMOldExcl):
self.trainedModel._PsExcl[j].pad(0, padRight)
self.trainedModel._HIsExcl[j].pad(0, padRight)
nmusMOld += nmusMOldExcl
for j in range(len(mus)):
nsj = self.samplingEngine.nsamples[nmusMOld + j]
padRight -= nsj
Ps[j].pad(padLeft, padRight)
padLeft += nsj
pMat = self.samplingEngine.samplesCoalesced.data
pMatEff = dot(self.HFEngine.C, pMat) if self.approx_state else pMat
self.trainedModel.data.projMat = pMatEff
self.trainedModel.data.Qs += Qs
self.trainedModel.data.Ps += Ps
self.trainedModel.data.approxParameters = copy(self.approxParameters)
self.verbosity += 15
@abstractmethod
def setupApproxPivoted(self, mus:paramList) -> int:
if self.checkComputedApproxPivoted(): return -1
RROMPyAssert(self._mode, message = "Cannot setup approximant.")
vbMng(self, "INIT", "Setting up pivoted approximant.", 10)
Qs, Ps = self._preSetupApproxPivoted()
pass
self._postSetupApproxPivoted(mus, Qs, Ps)
vbMng(self, "DEL", "Done setting up pivoted approximant.", 10)
return 0
def setupApprox(self, plotEst : str = "NONE") -> int:
"""Compute greedy snapshots of solution map."""
if self.checkComputedApprox(): return -1
RROMPyAssert(self._mode, message = "Cannot start greedy algorithm.")
vbMng(self, "INIT", "Starting computation of snapshots.", 3)
max2ErrorEst, self.firstGreedyIterM = np.inf, True
self._preliminaryTrainingMarginal()
if self.errorEstimatorKindMarginal[: 10] == "LOOK_AHEAD":
muidx = np.arange(len(self.trainedModel.data.musMarginal))
else:#if self.errorEstimatorKindMarginal in ["LEAVE_ONE_OUT", "NONE"]:
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)
vbMng(self, "MAIN", ("Uniform testing error estimate "
"{:.4e}.").format(max2ErrorEst), 3)
else:
max2ErrorEst = 0.
if plotEst == "LAST":
self.plotEstimatorMarginal(errorEstTest, muidx, maxErrorEst)
vbMng(self, "DEL",
("Done computing snapshots (final snapshot count: "
"{}).").format(np.sum(self.samplingEngine.nsamples)), 3)
if (self.errorEstimatorKindMarginal == "LOOK_AHEAD_RECOVER"
and 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': whether to compute POD of snapshots; defaults to True;
- 'scaleFactorDer': scaling factors for derivative computation;
defaults to 'AUTO';
- 'cutOffTolerance': tolerance for ignoring parasitic poles;
defaults to np.inf;
- 'matchingWeightError': weight for pole matching optimization in
error estimation; defaults to 0;
- 'cutOffToleranceError': tolerance for ignoring parasitic poles
in error estimation; defaults to 'AUTO', i.e. cutOffTolerance;
- 'S': total number of pivot samples current approximant relies
upon;
- 'samplerPivot': pivot sample point generator;
- 'SMarginal': number of starting marginal samples;
- 'samplerMarginal': marginal sample point generator via sparse
grid;
- 'errorEstimatorKindMarginal': kind of marginal error estimator;
available values include 'LEAVE_ONE_OUT', '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': whether to compute POD of snapshots;
- 'scaleFactorDer': scaling factors for derivative computation;
- 'cutOffTolerance': tolerance for ignoring parasitic poles;
- 'matchingWeightError': weight for pole matching optimization in
error estimation;
- 'cutOffToleranceError': tolerance for ignoring parasitic poles
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: Whether to compute POD of snapshots.
scaleFactorDer: Scaling factors for derivative computation.
cutOffTolerance: Tolerance for ignoring parasitic poles.
matchingWeightError: Weight for pole matching optimization in error
estimation.
cutOffToleranceError: Tolerance for ignoring parasitic poles in error
estimation.
S: Total number of pivot samples current approximant relies upon.
samplerPivot: Pivot sample point generator.
SMarginal: Total number of marginal samples current approximant relies
upon.
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': whether to compute POD of snapshots; defaults to True;
- 'scaleFactorDer': scaling factors for derivative computation;
defaults to 'AUTO';
- 'matchingWeight': weight for pole matching optimization; defaults
to 1;
- 'cutOffTolerance': tolerance for ignoring parasitic poles;
defaults to np.inf;
- 'cutOffSharedRatio': required ratio of marginal points to share
resonance in cut off strategy; defaults to 1.;
- 'matchingWeightError': weight for pole matching optimization in
error estimation; defaults to 0;
- 'cutOffToleranceError': tolerance for ignoring parasitic poles
in error estimation; defaults to 'AUTO', i.e. cutOffTolerance;
- 'S': total number of pivot samples current approximant relies
upon;
- 'samplerPivot': pivot sample point generator;
- 'SMarginal': number of starting marginal samples;
- 'samplerMarginal': marginal sample point generator via sparse
grid;
- 'errorEstimatorKindMarginal': kind of marginal error estimator;
available values include 'LEAVE_ONE_OUT', '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' or
'PIECEWISE_LINEAR_*';
. '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'.
- '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': whether to compute POD of snapshots;
- 'scaleFactorDer': scaling factors for derivative computation;
- 'matchingWeight': weight for pole matching optimization;
- 'cutOffTolerance': tolerance for ignoring parasitic poles;
- 'cutOffSharedRatio': required ratio of marginal points to share
resonance in cut off strategy;
- 'matchingWeightError': weight for pole matching optimization in
error estimation;
- 'cutOffToleranceError': tolerance for ignoring parasitic poles
in error estimation;
- 'errorEstimatorKindMarginal': kind of marginal error estimator;
- 'polybasisMarginal': type of polynomial basis for marginal
interpolation;
- '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.
- '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: Whether to compute POD of snapshots.
scaleFactorDer: Scaling factors for derivative computation.
matchingWeight: Weight for pole matching optimization.
cutOffTolerance: Tolerance for ignoring parasitic poles.
cutOffSharedRatio: Required ratio of marginal points to share resonance
in cut off strategy.
matchingWeightError: Weight for pole matching optimization in error
estimation.
cutOffToleranceError: Tolerance for ignoring parasitic poles in error
estimation.
S: Total number of pivot samples current approximant relies upon.
samplerPivot: Pivot sample point generator.
SMarginal: Total number of marginal samples current approximant relies
upon.
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.HFEngine,
self.matchingWeight, False)
vbMng(self, "DEL", "Done compressing and matching poles.", 10)
def _updateTrainedModelMarginalSamples(self, idx : ListAny = []):
self.trainedModel.updateEffectiveSamples(idx, self.HFEngine,
self.matchingWeight, False)
def getErrorEstimatorMarginalLeaveOneOut(self) -> Np1D:
if self.polybasisMarginal != "NEARESTNEIGHBOR":
if not hasattr(self, "_MMarginal_isauto"):
if not hasattr(self, "_MMarginalOriginal"):
self._MMarginalOriginal = self.paramsMarginal["MMarginal"]
self.paramsMarginal["MMarginal"] = self._MMarginalOriginal
self._reduceDegreeNNoWarn = 1
err = super().getErrorEstimatorMarginalLeaveOneOut()
if self.polybasisMarginal != "NEARESTNEIGHBOR":
del self._reduceDegreeNNoWarn
return err
def setupApprox(self, *args, **kwargs) -> int:
self.purgeparamsMarginal()
return super().setupApprox(*args, **kwargs)
diff --git a/rrompy/reduction_methods/pivoted/rational_interpolant_greedy_pivoted.py b/rrompy/reduction_methods/pivoted/rational_interpolant_greedy_pivoted.py
index ce15602..a9f21c4 100644
--- a/rrompy/reduction_methods/pivoted/rational_interpolant_greedy_pivoted.py
+++ b/rrompy/reduction_methods/pivoted/rational_interpolant_greedy_pivoted.py
@@ -1,624 +1,623 @@
# 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 rrompy.reduction_methods.standard.greedy.rational_interpolant_greedy \
import RationalInterpolantGreedy
from rrompy.reduction_methods.standard.greedy.generic_greedy_approximant \
import pruneSamples
from rrompy.utilities.base.types import Np1D
from rrompy.utilities.base import verbosityManager as vbMng
from rrompy.utilities.numerical import dot
from rrompy.utilities.numerical.degree import totalDegreeN
from rrompy.utilities.poly_fitting.polynomial import polyvander as pv
from rrompy.utilities.exception_manager import RROMPyAssert, RROMPyWarning
from rrompy.parameter import emptyParameterList, checkParameterList
__all__ = ['RationalInterpolantGreedyPivotedNoMatch',
'RationalInterpolantGreedyPivoted']
class RationalInterpolantGreedyPivotedBase(GenericPivotedApproximantBase,
RationalInterpolantGreedy):
def __init__(self, *args, **kwargs):
self._preInit()
self._addParametersToList(toBeExcluded = ["sampler"])
super().__init__(*args, **kwargs)
self._postInit()
@property
def tModelType(self):
if hasattr(self, "_temporaryPivot"):
return RationalInterpolantGreedy.tModelType.fget(self)
return super().tModelType
@property
def polybasis0(self):
if "_" in self.polybasis:
return self.polybasis.split("_")[0]
return self.polybasis
@property
def correctorTol(self):
"""Value of correctorTol."""
return self._correctorTol
@correctorTol.setter
def correctorTol(self, correctorTol):
if correctorTol < 0. or (correctorTol > 0. and self.nparPivot > 1):
RROMPyWarning(("Overriding prescribed corrector tolerance "
"to 0."))
correctorTol = 0.
self._correctorTol = correctorTol
self._approxParameters["correctorTol"] = self.correctorTol
@property
def correctorMaxIter(self):
"""Value of correctorMaxIter."""
return self._correctorMaxIter
@correctorMaxIter.setter
def correctorMaxIter(self, correctorMaxIter):
if correctorMaxIter < 1 or (correctorMaxIter > 1
and self.nparPivot > 1):
RROMPyWarning(("Overriding prescribed max number of corrector "
"iterations to 1."))
correctorMaxIter = 1
self._correctorMaxIter = correctorMaxIter
self._approxParameters["correctorMaxIter"] = self.correctorMaxIter
def _polyvanderAuxiliary(self, mus, deg, *args):
degEff = [0] * self.npar
degEff[self.directionPivot[0]] = deg
return pv(mus, degEff, *args)
def _marginalizeMiscellanea(self, forward:bool):
if forward:
self._m_mu0 = copy(self.mu0)
self._m_selfmus = copy(self.mus)
self._m_HFErescalingExp = copy(self.HFEngine.rescalingExp)
self._mu0 = checkParameterList(self.mu0(self.directionPivot), 1)[0]
self._mus = checkParameterList(self.mus(self.directionPivot), 1)[0]
self.HFEngine.rescalingExp = [self.HFEngine.rescalingExp[
self.directionPivot[0]]]
else:
self._mu0 = self._m_mu0
self._mus = self._m_selfmus
self.HFEngine.rescalingExp = self._m_HFErescalingExp
del self._m_mu0, self._m_selfmus, self._m_HFErescalingExp
def _marginalizeTrainedModel(self, forward:bool):
if forward:
del self._temporaryPivot
self.trainedModel.data.mu0 = self.mu0
self.trainedModel.data.scaleFactor = [1.] * self.npar
self.trainedModel.data.scaleFactor[self.directionPivot[0]] = (
self.scaleFactor[0])
self.trainedModel.data.rescalingExp = self.HFEngine.rescalingExp
Qc = np.zeros((1,) * self.directionPivot[0]
+ (len(self.trainedModel.data.Q.coeffs),)
+ (1,) * (self.npar - self.directionPivot[0] - 1),
dtype = self.trainedModel.data.Q.coeffs.dtype)
Pc = np.zeros((1,) * self.directionPivot[0]
+ (len(self.trainedModel.data.P.coeffs),)
+ (1,) * (self.npar - self.directionPivot[0] - 1)
+ (self.trainedModel.data.P.coeffs.shape[1],),
dtype = self.trainedModel.data.P.coeffs.dtype)
for j in range(len(self.trainedModel.data.Q.coeffs)):
Qc[(0,) * self.directionPivot[0] + (j,)
+ (0,) * (self.npar - self.directionPivot[0] - 1)] = (
self.trainedModel.data.Q.coeffs[j])
for j in range(len(self.trainedModel.data.P.coeffs)):
for k in range(self.trainedModel.data.P.coeffs.shape[1]):
Pc[(0,) * self.directionPivot[0] + (j,)
+ (0,) * (self.npar - self.directionPivot[0] - 1)
+ (k,)] = self.trainedModel.data.P.coeffs[j, k]
self.trainedModel.data.Q.coeffs = Qc
self.trainedModel.data.P.coeffs = Pc
self._m_musUniqueCN = copy(self._musUniqueCN)
musUniqueCNAux = np.zeros((self.S, self.npar),
dtype = self._musUniqueCN.dtype)
musUniqueCNAux[:, self.directionPivot[0]] = self._musUniqueCN(0)
self._musUniqueCN = checkParameterList(musUniqueCNAux,
self.npar)[0]
self._m_derIdxs = copy(self._derIdxs)
for j in range(len(self._derIdxs)):
for l in range(len(self._derIdxs[j])):
derjl = self._derIdxs[j][l][0]
self._derIdxs[j][l] = [0] * self.npar
self._derIdxs[j][l][self.directionPivot[0]] = derjl
else:
self._temporaryPivot = 1
self.trainedModel.data.mu0 = checkParameterList(
self.mu0(self.directionPivot), 1)[0]
self.trainedModel.data.scaleFactor = self.scaleFactor
self.trainedModel.data.rescalingExp = self.HFEngine.rescalingExp[
self.directionPivot[0]]
self.trainedModel.data.Q.coeffs = self.trainedModel.data.Q.coeffs[
(0,) * self.directionPivot[0]
+ (slice(None),)
+ (0,) * (self.HFEngine.npar - 1
- self.directionPivot[0])]
self.trainedModel.data.P.coeffs = self.trainedModel.data.P.coeffs[
(0,) * self.directionPivot[0]
+ (slice(None),)
+ (0,) * (self.HFEngine.npar - 1
- self.directionPivot[0])]
self._musUniqueCN = copy(self._m_musUniqueCN)
self._derIdxs = copy(self._m_derIdxs)
del self._m_musUniqueCN, self._m_derIdxs
self.trainedModel.data.npar = self.npar
self.trainedModel.data.Q.npar = self.npar
self.trainedModel.data.P.npar = self.npar
def errorEstimator(self, mus:Np1D, return_max : bool = False) -> Np1D:
"""Standard residual-based error estimator."""
self._marginalizeMiscellanea(True)
setupOK = self.setupApproxLocal()
self._marginalizeMiscellanea(False)
if setupOK > 0:
err = np.empty(len(mus))
err[:] = np.nan
if not return_max: return err
return err, [- setupOK], np.nan
self._marginalizeTrainedModel(True)
errRes = super().errorEstimator(mus, return_max)
self._marginalizeTrainedModel(False)
return errRes
def _preliminaryTraining(self):
"""Initialize starting snapshots of solution map."""
RROMPyAssert(self._mode, message = "Cannot start greedy algorithm.")
S = self.S
self._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(muTestPivot ** self.HFEngine.rescalingExp[
self.directionPivot[0]],
musPivot ** self.HFEngine.rescalingExp[
self.directionPivot[0]],
1e-10 * self.scaleFactor[0])
self.mus = emptyParameterList()
self.mus.reset((self.S, self.npar + len(self.musMargLoc)))
muTestBase = emptyParameterList()
muTestBase.reset((len(muTestPivot), self.npar + len(self.musMargLoc)))
for k in range(self.S):
self.mus.data[k, self.directionPivot] = musPivot[k].data
self.mus.data[k, self.directionMarginal] = self.musMargLoc.data
for k in range(len(muTestPivot)):
muTestBase.data[k, self.directionPivot] = muTestPivot[k].data
muTestBase.data[k, self.directionMarginal] = self.musMargLoc.data
muTestBase.pop(idxPop)
muLast = copy(self.mus[-1])
self.mus.pop()
if len(self.mus) > 0:
vbMng(self, "MAIN",
("Adding first {} sample point{} at {} to training "
"set.").format(self.S - 1, "" + "s" * (self.S > 2),
self.mus), 3)
self.samplingEngine.iterSample(self.mus)
self._S = len(self.mus)
self._approxParameters["S"] = self.S
self.muTest = emptyParameterList()
self.muTest.reset((len(muTestBase) + 1, self.mus.shape[1]))
self.muTest.data[: -1] = muTestBase.data
self.muTest.data[-1] = muLast.data
self.M, self.N = ("AUTO",) * 2
def _finalizeSnapshots(self):
self.setupSampling()
self.samplingEngine.resetHistory(len(self.musMarginal))
for j in range(len(self.musMarginal)):
self.samplingEngine.setsample(self.samplingEngs[j].samples,
j, False)
self.samplingEngine.mus[j] = copy(self.samplingEngs[j].mus)
self.samplingEngine.musMarginal[j] = copy(self.musMarginal[j])
self.samplingEngine.nsamples[j] = self.samplingEngs[j].nsamples
if self.POD:
self.samplingEngine.RPOD[j] = self.samplingEngs[j].RPOD
self.samplingEngine.samples_full[j].data = (
self.samplingEngs[j].samples_full.data)
- self.samplingEngine.coalesceSamples()
def setupApprox(self, *args, **kwargs) -> int:
"""Compute rational interpolant."""
if self.checkComputedApprox(): return -1
RROMPyAssert(self._mode, message = "Cannot setup approximant.")
vbMng(self, "INIT", "Setting up {}.". format(self.name()), 5)
self._musMarginal = self.samplerMarginal.generatePoints(self.SMarginal)
while len(self.musMarginal) > self.SMarginal: self.musMarginal.pop()
S0 = copy(self.S)
Qs, Ps = [None] * len(self.musMarginal), [None] * len(self.musMarginal)
self.samplingEngs = [None] * len(self.musMarginal)
self.computeScaleFactor()
self._scaleFactorOldPivot = copy(self.scaleFactor)
self.scaleFactor = self.scaleFactorPivot
self._temporaryPivot = 1
for j in range(len(self.musMarginal)):
vbMng(self, "MAIN",
"Building marginal model no. {} at {}.".format(j + 1,
self.musMarginal[j]), 5)
self._S = S0
self.musMargLoc = self.musMarginal[j]
RationalInterpolantGreedy.setupSampling(self)
self.trainedModel = None
self.verbosity -= 5
self.samplingEngine.verbosity -= 5
super().setupApprox(*args, **kwargs)
self.verbosity += 5
self.samplingEngine.verbosity += 5
self.samplingEngs[j] = copy(self.samplingEngine)
Qs[j] = copy(self.trainedModel.data.Q)
Ps[j] = copy(self.trainedModel.data.P)
self.scaleFactor = self._scaleFactorOldPivot
del self._scaleFactorOldPivot, self._temporaryPivot
self._finalizeSnapshots()
del self.musMargLoc, self.samplingEngs
self._mus = self.samplingEngine.musCoalesced
- padLeft, padRight = 0, self.samplingEngine.nsamplesTot
+ padLeft, padRight = 0, self.samplingEngine.nsamplesCoalesced
for j in range(len(self.musMarginal)):
nsj = self.samplingEngine.nsamples[j]
padRight -= nsj
Ps[j].pad(padLeft, padRight)
padLeft += nsj
pMat = self.samplingEngine.samplesCoalesced.data
pMatEff = dot(self.HFEngine.C, pMat) if self.approx_state else pMat
self.trainedModel = self.tModelType()
self.trainedModel.verbosity = self.verbosity
self.trainedModel.timestamp = self.timestamp
datadict = {"mu0": self.mu0, "mus": copy(self.mus), "projMat": pMatEff,
"scaleFactor": self.scaleFactor,
"rescalingExp": self.HFEngine.rescalingExp,
"directionPivot": self.directionPivot}
self.trainedModel.data = self.initializeModelData(datadict)[0]
self.trainedModel.data.musMarginal = copy(self.musMarginal)
self.trainedModel.data.Qs, self.trainedModel.data.Ps = Qs, Ps
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': whether to compute POD of snapshots; defaults to True;
- 'scaleFactorDer': scaling factors for derivative computation;
defaults to 'AUTO';
- 'cutOffTolerance': tolerance for ignoring parasitic poles;
defaults to np.inf;
- '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;
- 'interpRcond': tolerance for pivot interpolation; defaults to
None;
- 'robustTol': tolerance for robust rational denominator
management; defaults to 0;
- 'correctorForce': whether corrector should forcefully delete bad
poles; defaults to False;
- 'correctorTol': tolerance for corrector step; defaults to 0.,
i.e. no bad poles;
- 'correctorMaxIter': maximum number of corrector iterations;
defaults to 1.
Defaults to empty dict.
approx_state(optional): Whether to approximate state. Defaults to
False.
verbosity(optional): Verbosity level. Defaults to 10.
Attributes:
HFEngine: HF problem solver.
mu0: Default parameter.
directionPivot: Pivot components.
mus: Array of snapshot parameters.
musMarginal: Array of marginal snapshot parameters.
approxParameters: Dictionary containing values for main parameters of
approximant. Recognized keys are in parameterList.
parameterListSoft: Recognized keys of soft approximant parameters:
- 'POD': whether to compute POD of snapshots;
- 'scaleFactorDer': scaling factors for derivative computation;
- 'cutOffTolerance': tolerance for ignoring parasitic poles;
- '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;
- 'interpRcond': tolerance for pivot interpolation;
- 'robustTol': tolerance for robust rational denominator
management;
- 'correctorForce': whether corrector should forcefully delete bad
poles;
- 'correctorTol': tolerance for corrector step;
- 'correctorMaxIter': maximum number of corrector iterations.
parameterListCritical: Recognized keys of critical approximant
parameters:
- 'S': total number of pivot samples current approximant relies
upon;
- 'samplerPivot': pivot sample point generator;
- 'SMarginal': total number of marginal samples current approximant
relies upon;
- 'samplerMarginal': marginal sample point generator.
approx_state: Whether to approximate state.
verbosity: Verbosity level.
POD: Whether to compute POD of snapshots.
scaleFactorDer: Scaling factors for derivative computation.
cutOffTolerance: Tolerance for ignoring parasitic poles.
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.
interpRcond: Tolerance for pivot interpolation.
robustTol: Tolerance for robust rational denominator management.
correctorForce: Whether corrector should forcefully delete bad poles.
correctorTol: Tolerance for corrector step.
correctorMaxIter: Maximum number of corrector iterations.
muBounds: list of bounds for pivot parameter values.
muBoundsMarginal: list of bounds for marginal parameter values.
samplingEngine: Sampling engine.
uHF: High fidelity solution(s) with parameter(s) lastSolvedHF as
sampleList.
lastSolvedHF: Parameter(s) corresponding to last computed high fidelity
solution(s) as parameterList.
uApproxReduced: Reduced approximate solution(s) with parameter(s)
lastSolvedApprox as sampleList.
lastSolvedApproxReduced: Parameter(s) corresponding to last computed
reduced approximate solution(s) as parameterList.
uApprox: Approximate solution(s) with parameter(s) lastSolvedApprox as
sampleList.
lastSolvedApprox: Parameter(s) corresponding to last computed
approximate solution(s) as parameterList.
Q: Numpy 1D vector containing complex coefficients of approximant
denominator.
P: Numpy 2D vector whose columns are FE dofs of coefficients of
approximant numerator.
"""
def _poleMatching(self):
vbMng(self, "INIT", "Compressing poles.", 10)
self.trainedModel.initializeFromRational()
vbMng(self, "DEL", "Done compressing poles.", 10)
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': whether to compute POD of snapshots; defaults to True;
- 'scaleFactorDer': scaling factors for derivative computation;
defaults to 'AUTO';
- 'matchingWeight': weight for pole matching optimization; defaults
to 1;
- 'cutOffTolerance': tolerance for ignoring parasitic poles;
defaults to np.inf;
- 'cutOffSharedRatio': required ratio of marginal points to share
resonance in cut off strategy; defaults to 1.;
- 'S': total number of pivot samples current approximant relies
upon;
- 'samplerPivot': pivot sample point generator;
- 'SMarginal': total number of marginal samples current approximant
relies upon;
- 'samplerMarginal': marginal sample point generator;
- 'polybasis': type of polynomial basis for pivot
interpolation; defaults to 'MONOMIAL';
- 'polybasisMarginal': type of polynomial basis for marginal
interpolation; allowed values include 'MONOMIAL_*',
'CHEBYSHEV_*', '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' or
'PIECEWISE_LINEAR_*';
. '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'.
- '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;
- 'interpRcond': tolerance for pivot interpolation; defaults to
None;
- 'robustTol': tolerance for robust rational denominator
management; defaults to 0;
- 'correctorForce': whether corrector should forcefully delete bad
poles; defaults to False;
- 'correctorTol': tolerance for corrector step; defaults to 0.,
i.e. no bad poles;
- 'correctorMaxIter': maximum number of corrector iterations;
defaults to 1.
Defaults to empty dict.
approx_state(optional): Whether to approximate state. Defaults to
False.
verbosity(optional): Verbosity level. Defaults to 10.
Attributes:
HFEngine: HF problem solver.
mu0: Default parameter.
directionPivot: Pivot components.
mus: Array of snapshot parameters.
musMarginal: Array of marginal snapshot parameters.
approxParameters: Dictionary containing values for main parameters of
approximant. Recognized keys are in parameterList.
parameterListSoft: Recognized keys of soft approximant parameters:
- 'POD': whether to compute POD of snapshots;
- 'scaleFactorDer': scaling factors for derivative computation;
- 'matchingWeight': weight for pole matching optimization;
- 'cutOffTolerance': tolerance for ignoring parasitic poles;
- 'cutOffSharedRatio': required ratio of marginal points to share
resonance in cut off strategy;
- 'polybasis': type of polynomial basis for pivot
interpolation;
- 'polybasisMarginal': type of polynomial basis for marginal
interpolation;
- '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.
- '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;
- 'interpRcond': tolerance for pivot interpolation;
- 'robustTol': tolerance for robust rational denominator
management;
- 'correctorForce': whether corrector should forcefully delete bad
poles;
- 'correctorTol': tolerance for corrector step;
- 'correctorMaxIter': maximum number of corrector iterations.
parameterListCritical: Recognized keys of critical approximant
parameters:
- 'S': total number of pivot samples current approximant relies
upon;
- 'samplerPivot': pivot sample point generator;
- 'SMarginal': total number of marginal samples current approximant
relies upon;
- 'samplerMarginal': marginal sample point generator.
approx_state: Whether to approximate state.
verbosity: Verbosity level.
POD: Whether to compute POD of snapshots.
scaleFactorDer: Scaling factors for derivative computation.
matchingWeight: Weight for pole matching optimization.
cutOffTolerance: Tolerance for ignoring parasitic poles.
cutOffSharedRatio: Required ratio of marginal points to share resonance
in cut off strategy.
S: Total number of pivot samples current approximant relies upon.
samplerPivot: Pivot sample point generator.
SMarginal: Total number of marginal samples current approximant relies
upon.
samplerMarginal: Marginal sample point generator.
polybasis: Type of polynomial basis for pivot interpolation.
polybasisMarginal: Type of polynomial basis for marginal interpolation.
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.
interpRcond: Tolerance for pivot interpolation.
robustTol: Tolerance for robust rational denominator management.
correctorForce: Whether corrector should forcefully delete bad poles.
correctorTol: Tolerance for corrector step.
correctorMaxIter: Maximum number of corrector iterations.
muBounds: list of bounds for pivot parameter values.
muBoundsMarginal: list of bounds for marginal parameter values.
samplingEngine: Sampling engine.
uHF: High fidelity solution(s) with parameter(s) lastSolvedHF as
sampleList.
lastSolvedHF: Parameter(s) corresponding to last computed high fidelity
solution(s) as parameterList.
uApproxReduced: Reduced approximate solution(s) with parameter(s)
lastSolvedApprox as sampleList.
lastSolvedApproxReduced: Parameter(s) corresponding to last computed
reduced approximate solution(s) as parameterList.
uApprox: Approximate solution(s) with parameter(s) lastSolvedApprox as
sampleList.
lastSolvedApprox: Parameter(s) corresponding to last computed
approximate solution(s) as parameterList.
Q: Numpy 1D vector containing complex coefficients of approximant
denominator.
P: Numpy 2D vector whose columns are FE dofs of coefficients of
approximant numerator.
"""
def _poleMatching(self):
vbMng(self, "INIT", "Compressing and matching poles.", 10)
self.trainedModel.initializeFromRational(self.HFEngine,
self.matchingWeight, False)
vbMng(self, "DEL", "Done compressing and matching poles.", 10)
def setupApprox(self, *args, **kwargs) -> int:
self.purgeparamsMarginal()
return super().setupApprox(*args, **kwargs)
diff --git a/rrompy/reduction_methods/pivoted/rational_interpolant_pivoted.py b/rrompy/reduction_methods/pivoted/rational_interpolant_pivoted.py
index 6245c33..8eb1c73 100644
--- a/rrompy/reduction_methods/pivoted/rational_interpolant_pivoted.py
+++ b/rrompy/reduction_methods/pivoted/rational_interpolant_pivoted.py
@@ -1,525 +1,519 @@
# 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 rrompy.reduction_methods.standard.rational_interpolant import (
RationalInterpolant)
from rrompy.utilities.base import verbosityManager as vbMng
from rrompy.utilities.numerical import dot
from rrompy.utilities.numerical.hash_derivative import nextDerivativeIndices
from rrompy.utilities.exception_manager import RROMPyAssert, RROMPyWarning
from rrompy.parameter import emptyParameterList
__all__ = ['RationalInterpolantPivotedNoMatch', 'RationalInterpolantPivoted']
class RationalInterpolantPivotedBase(GenericPivotedApproximantBase,
RationalInterpolant):
def __init__(self, *args, **kwargs):
self._preInit()
self._addParametersToList(toBeExcluded = ["polydegreetype", "sampler"])
super().__init__(*args, **kwargs)
self._postInit()
@property
def scaleFactorDer(self):
"""Value of scaleFactorDer."""
if self._scaleFactorDer == "NONE": return 1.
if self._scaleFactorDer == "AUTO": return self.scaleFactorPivot
return self._scaleFactorDer
@scaleFactorDer.setter
def scaleFactorDer(self, scaleFactorDer):
if hasattr(self, "_scaleFactorDer"):
scaleFactorDerold = self.scaleFactorDer
else: scaleFactorDerold = -1
if isinstance(scaleFactorDer, (str,)):
scaleFactorDer = scaleFactorDer.upper()
self._scaleFactorDer = scaleFactorDer
self._approxParameters["scaleFactorDer"] = self._scaleFactorDer
if scaleFactorDerold != self._scaleFactorDer: self.resetSamples()
@property
def polydegreetype(self):
"""Value of polydegreetype."""
return "TOTAL"
@polydegreetype.setter
def polydegreetype(self, polydegreetype):
RROMPyWarning(("polydegreetype is used just to simplify inheritance, "
"and its value cannot be changed from 'TOTAL'."))
@property
def polybasis0(self):
if "_" in self.polybasis:
return self.polybasis.split("_")[0]
return self.polybasis
@property
def correctorTol(self):
"""Value of correctorTol."""
return self._correctorTol
@correctorTol.setter
def correctorTol(self, correctorTol):
if correctorTol < 0. or (correctorTol > 0. and self.nparPivot > 1):
RROMPyWarning(("Overriding prescribed corrector tolerance "
"to 0."))
correctorTol = 0.
self._correctorTol = correctorTol
self._approxParameters["correctorTol"] = self.correctorTol
@property
def correctorMaxIter(self):
"""Value of correctorMaxIter."""
return self._correctorMaxIter
@correctorMaxIter.setter
def correctorMaxIter(self, correctorMaxIter):
if correctorMaxIter < 1 or (correctorMaxIter > 1
and self.nparPivot > 1):
RROMPyWarning(("Overriding prescribed max number of corrector "
"iterations to 1."))
correctorMaxIter = 1
self._correctorMaxIter = correctorMaxIter
self._approxParameters["correctorMaxIter"] = self.correctorMaxIter
def _setupInterpolationIndices(self):
"""Setup parameters for polyvander."""
RROMPyAssert(self._mode,
message = "Cannot setup interpolation indices.")
if (self._musUniqueCN is None
or len(self._reorder) != len(self.musPivot)):
try:
muPC = self.trainedModel.centerNormalizePivot(self.musPivot)
except:
muPC = self.trainedModel.centerNormalize(self.musPivot)
self._musUniqueCN, musIdxsTo, musIdxs, musCount = (muPC.unique(
return_index = True, return_inverse = True,
return_counts = True))
self._musUnique = self.musPivot[musIdxsTo]
self._derIdxs = [None] * len(self._musUniqueCN)
self._reorder = np.empty(len(musIdxs), dtype = int)
filled = 0
for j, cnt in enumerate(musCount):
self._derIdxs[j] = nextDerivativeIndices([], self.nparPivot,
cnt)
jIdx = np.nonzero(musIdxs == j)[0]
self._reorder[jIdx] = np.arange(filled, filled + cnt)
filled += cnt
def computeSnapshots(self):
"""Compute snapshots of solution map."""
RROMPyAssert(self._mode,
message = "Cannot start snapshot computation.")
self.computeScaleFactor()
- if self.samplingEngine.nsamplesTot != self.S * self.SMarginal:
+ if self.samplingEngine.nsamplesCoalesced != self.S * self.SMarginal:
self.resetSamples()
self.samplingEngine.scaleFactor = self.scaleFactorDer
vbMng(self, "INIT", "Starting computation of snapshots.", 5)
self.musPivot = self.samplerPivot.generatePoints(self.S)
while len(self.musPivot) > self.S: self.musPivot.pop()
self._musMarginal = self.samplerMarginal.generatePoints(
self.SMarginal)
while len(self.musMarginal) > self.SMarginal:
self.musMarginal.pop()
self.mus = emptyParameterList()
self.mus.reset((self.S * self.SMarginal, self.HFEngine.npar))
self.samplingEngine.resetHistory(self.SMarginal)
for j, muMarg in enumerate(self.musMarginal):
for k in range(j * self.S, (j + 1) * self.S):
self.mus.data[k, self.directionPivot] = (
self.musPivot[k - j * self.S].data)
self.mus.data[k, self.directionMarginal] = muMarg.data
self.samplingEngine.iterSample(self.musPivot, self.musMarginal)
- self._finalizeSnapshots()
vbMng(self, "DEL", "Done computing snapshots.", 5)
- def _finalizeSnapshots(self):
- self.samplingEngine.coalesceSamples()
-
def setupApprox(self) -> int:
"""Compute rational interpolant."""
if self.checkComputedApprox(): return -1
RROMPyAssert(self._mode, message = "Cannot setup approximant.")
vbMng(self, "INIT", "Setting up {}.". format(self.name()), 5)
self.computeSnapshots()
pMat = self.samplingEngine.samplesCoalesced.data
pMatEff = dot(self.HFEngine.C, pMat) if self.approx_state else pMat
if self.trainedModel is None:
self.trainedModel = self.tModelType()
self.trainedModel.verbosity = self.verbosity
self.trainedModel.timestamp = self.timestamp
datadict = {"mu0": self.mu0, "mus": copy(self.mus),
"projMat": pMatEff, "scaleFactor": self.scaleFactor,
"rescalingExp": self.HFEngine.rescalingExp,
"directionPivot": self.directionPivot}
self.trainedModel.data = self.initializeModelData(datadict)[0]
else:
self.trainedModel = self.trainedModel
self.trainedModel.data.projMat = copy(pMatEff)
self.trainedModel.data.mus = copy(self.mus)
N0 = copy(self.N)
Qs, Ps = [None] * len(self.musMarginal), [None] * len(self.musMarginal)
self._temporaryPivot = 1
- padLeft = 0
+ padLeft, padRight = 0, self.samplingEngine.nsamplesCoalesced
if self.POD:
- self._RPODOldPivot = copy(self.samplingEngine.RPODCoalesced)
+ self._RPODOldPivot = copy(self.samplingEngine.RPOD)
else:
self._samplesOldPivot = copy(self.samplingEngine.samples)
- padRight = self.samplingEngine.nsamplesTot
self._scaleFactorOldPivot = copy(self.scaleFactor)
self.scaleFactor = self.scaleFactorPivot
for j in range(len(self.musMarginal)):
vbMng(self, "MAIN",
"Building marginal model no. {} at {}.".format(j + 1,
self.musMarginal[j]), 5)
self.N = N0
if self.POD:
- self.samplingEngine.RPOD = (
- self._RPODOldPivot[:, padLeft : padLeft + self.S])
+ self.samplingEngine.RPOD = self._RPODOldPivot[j]
else:
self.samplingEngine.samples = self._samplesOldPivot[j]
- padRight -= self.S
self.verbosity -= 5
self._iterCorrector()
self.verbosity += 5
Qs[j] = copy(self.trainedModel.data.Q)
Ps[j] = copy(self.trainedModel.data.P)
del self.trainedModel.data.Q, self.trainedModel.data.P
- if not self.POD: Ps[j].pad(padLeft, padRight)
+ padRight -= self.S
+ Ps[j].pad(padLeft, padRight)
padLeft += self.S
if self.POD:
- self.samplingEngine.RPODCoalesced = copy(self._RPODOldPivot)
+ self.samplingEngine.RPOD = copy(self._RPODOldPivot)
del self._RPODOldPivot
else:
self.samplingEngine.samples = copy(self._samplesOldPivot)
del self._samplesOldPivot
self.scaleFactor = self._scaleFactorOldPivot
del self._temporaryPivot, self._scaleFactorOldPivot
self.trainedModel.data.musPivot = copy(self.musPivot)
self.trainedModel.data.musMarginal = copy(self.musMarginal)
self.trainedModel.data.Qs, self.trainedModel.data.Ps = Qs, Ps
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': whether to compute POD of snapshots; defaults to True;
- 'scaleFactorDer': scaling factors for derivative computation;
defaults to 'AUTO';
- 'cutOffTolerance': tolerance for ignoring parasitic poles;
defaults to np.inf;
- '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;
- 'radialDirectionalWeightsMarginal': radial basis weights for
marginal interpolant; defaults to 1;
- 'interpRcond': tolerance for pivot interpolation; defaults to
None;
- 'robustTol': tolerance for robust rational denominator
management; defaults to 0;
- 'correctorForce': whether corrector should forcefully delete bad
poles; defaults to False;
- 'correctorTol': tolerance for corrector step; defaults to 0.,
i.e. no bad poles;
- 'correctorMaxIter': maximum number of corrector iterations;
defaults to 1.
Defaults to empty dict.
approx_state(optional): Whether to approximate state. Defaults to
False.
verbosity(optional): Verbosity level. Defaults to 10.
Attributes:
HFEngine: HF problem solver.
mu0: Default parameter.
directionPivot: Pivot components.
mus: Array of snapshot parameters.
musPivot: Array of pivot snapshot parameters.
musMarginal: Array of marginal snapshot parameters.
approxParameters: Dictionary containing values for main parameters of
approximant. Recognized keys are in parameterList.
parameterListSoft: Recognized keys of soft approximant parameters:
- 'POD': whether to compute POD of snapshots;
- 'scaleFactorDer': scaling factors for derivative computation;
- 'cutOffTolerance': tolerance for ignoring parasitic poles;
- '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;
- 'radialDirectionalWeightsMarginal': radial basis weights for
marginal interpolant;
- 'interpRcond': tolerance for pivot interpolation;
- 'robustTol': tolerance for robust rational denominator
management;
- 'correctorForce': whether corrector should forcefully delete bad
poles;
- 'correctorTol': tolerance for corrector step;
- 'correctorMaxIter': maximum number of corrector iterations.
parameterListCritical: Recognized keys of critical approximant
parameters:
- 'S': total number of pivot samples current approximant relies
upon;
- 'samplerPivot': pivot sample point generator;
- 'SMarginal': total number of marginal samples current approximant
relies upon;
- 'samplerMarginal': marginal sample point generator.
approx_state: Whether to approximate state.
verbosity: Verbosity level.
POD: Whether to compute POD of snapshots.
scaleFactorDer: Scaling factors for derivative computation.
cutOffTolerance: Tolerance for ignoring parasitic poles.
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.
radialDirectionalWeightsMarginal: Radial basis weights for marginal
interpolant.
interpRcond: Tolerance for pivot interpolation.
robustTol: Tolerance for robust rational denominator management.
correctorForce: Whether corrector should forcefully delete bad poles.
correctorTol: Tolerance for corrector step.
correctorMaxIter: Maximum number of corrector iterations.
muBounds: list of bounds for pivot parameter values.
muBoundsMarginal: list of bounds for marginal parameter values.
samplingEngine: Sampling engine.
uHF: High fidelity solution(s) with parameter(s) lastSolvedHF as
sampleList.
lastSolvedHF: Parameter(s) corresponding to last computed high fidelity
solution(s) as parameterList.
uApproxReduced: Reduced approximate solution(s) with parameter(s)
lastSolvedApprox as sampleList.
lastSolvedApproxReduced: Parameter(s) corresponding to last computed
reduced approximate solution(s) as parameterList.
uApprox: Approximate solution(s) with parameter(s) lastSolvedApprox as
sampleList.
lastSolvedApprox: Parameter(s) corresponding to last computed
approximate solution(s) as parameterList.
Q: Numpy 1D vector containing complex coefficients of approximant
denominator.
P: Numpy 2D vector whose columns are FE dofs of coefficients of
approximant numerator.
"""
def _poleMatching(self):
vbMng(self, "INIT", "Compressing poles.", 10)
self.trainedModel.initializeFromRational()
vbMng(self, "DEL", "Done compressing poles.", 10)
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': whether to compute POD of snapshots; defaults to True;
- 'scaleFactorDer': scaling factors for derivative computation;
defaults to 'AUTO';
- 'matchingWeight': weight for pole matching optimization; defaults
to 1;
- 'cutOffTolerance': tolerance for ignoring parasitic poles;
defaults to np.inf;
- 'cutOffSharedRatio': required ratio of marginal points to share
resonance in cut off strategy; defaults to 1.;
- 'S': total number of pivot samples current approximant relies
upon;
- 'samplerPivot': pivot sample point generator;
- 'SMarginal': total number of marginal samples current approximant
relies upon;
- 'samplerMarginal': marginal sample point generator;
- 'polybasis': type of polynomial basis for pivot
interpolation; defaults to 'MONOMIAL';
- 'polybasisMarginal': type of polynomial basis for marginal
interpolation; allowed values include 'MONOMIAL_*',
'CHEBYSHEV_*', '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' or
'PIECEWISE_LINEAR_*';
. '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'.
- '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;
- 'radialDirectionalWeightsMarginal': radial basis weights for
marginal interpolant; defaults to 1;
- 'interpRcond': tolerance for pivot interpolation; defaults to
None;
- 'robustTol': tolerance for robust rational denominator
management; defaults to 0;
- 'correctorForce': whether corrector should forcefully delete bad
poles; defaults to False;
- 'correctorTol': tolerance for corrector step; defaults to 0.,
i.e. no bad poles;
- 'correctorMaxIter': maximum number of corrector iterations;
defaults to 1.
Defaults to empty dict.
approx_state(optional): Whether to approximate state. Defaults to
False.
verbosity(optional): Verbosity level. Defaults to 10.
Attributes:
HFEngine: HF problem solver.
mu0: Default parameter.
directionPivot: Pivot components.
mus: Array of snapshot parameters.
musPivot: Array of pivot snapshot parameters.
musMarginal: Array of marginal snapshot parameters.
approxParameters: Dictionary containing values for main parameters of
approximant. Recognized keys are in parameterList.
parameterListSoft: Recognized keys of soft approximant parameters:
- 'POD': whether to compute POD of snapshots;
- 'scaleFactorDer': scaling factors for derivative computation;
- 'matchingWeight': weight for pole matching optimization;
- 'cutOffTolerance': tolerance for ignoring parasitic poles;
- 'cutOffSharedRatio': required ratio of marginal points to share
resonance in cut off strategy;
- 'polybasis': type of polynomial basis for pivot
interpolation;
- 'polybasisMarginal': type of polynomial basis for marginal
interpolation;
- '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.
- 'M': degree of rational interpolant numerator;
- 'N': degree of rational interpolant denominator;
- 'radialDirectionalWeights': radial basis weights for pivot
numerator;
- 'radialDirectionalWeightsMarginal': radial basis weights for
marginal interpolant;
- 'interpRcond': tolerance for pivot interpolation;
- 'robustTol': tolerance for robust rational denominator
management;
- 'correctorForce': whether corrector should forcefully delete bad
poles;
- 'correctorTol': tolerance for corrector step;
- 'correctorMaxIter': maximum number of corrector iterations.
parameterListCritical: Recognized keys of critical approximant
parameters:
- 'S': total number of pivot samples current approximant relies
upon;
- 'samplerPivot': pivot sample point generator;
- 'SMarginal': total number of marginal samples current approximant
relies upon;
- 'samplerMarginal': marginal sample point generator.
approx_state: Whether to approximate state.
verbosity: Verbosity level.
POD: Whether to compute POD of snapshots.
scaleFactorDer: Scaling factors for derivative computation.
matchingWeight: Weight for pole matching optimization.
cutOffTolerance: Tolerance for ignoring parasitic poles.
cutOffSharedRatio: Required ratio of marginal points to share resonance
in cut off strategy.
S: Total number of pivot samples current approximant relies upon.
samplerPivot: Pivot sample point generator.
SMarginal: Total number of marginal samples current approximant relies
upon.
samplerMarginal: Marginal sample point generator.
polybasis: Type of polynomial basis for pivot interpolation.
polybasisMarginal: Type of polynomial basis for marginal interpolation.
paramsMarginal: Dictionary of parameters for marginal interpolation.
M: Numerator degree of approximant.
N: Denominator degree of approximant.
radialDirectionalWeights: Radial basis weights for pivot numerator.
radialDirectionalWeightsMarginal: Radial basis weights for marginal
interpolant.
interpRcond: Tolerance for pivot interpolation.
robustTol: Tolerance for robust rational denominator management.
correctorForce: Whether corrector should forcefully delete bad poles.
correctorTol: Tolerance for corrector step.
correctorMaxIter: Maximum number of corrector iterations.
muBounds: list of bounds for pivot parameter values.
muBoundsMarginal: list of bounds for marginal parameter values.
samplingEngine: Sampling engine.
uHF: High fidelity solution(s) with parameter(s) lastSolvedHF as
sampleList.
lastSolvedHF: Parameter(s) corresponding to last computed high fidelity
solution(s) as parameterList.
uApproxReduced: Reduced approximate solution(s) with parameter(s)
lastSolvedApprox as sampleList.
lastSolvedApproxReduced: Parameter(s) corresponding to last computed
reduced approximate solution(s) as parameterList.
uApprox: Approximate solution(s) with parameter(s) lastSolvedApprox as
sampleList.
lastSolvedApprox: Parameter(s) corresponding to last computed
approximate solution(s) as parameterList.
Q: Numpy 1D vector containing complex coefficients of approximant
denominator.
P: Numpy 2D vector whose columns are FE dofs of coefficients of
approximant numerator.
"""
def _poleMatching(self):
vbMng(self, "INIT", "Compressing and matching poles.", 10)
self.trainedModel.initializeFromRational(self.HFEngine,
self.matchingWeight, False)
vbMng(self, "DEL", "Done compressing and matching poles.", 10)
def setupApprox(self, *args, **kwargs) -> int:
self.purgeparamsMarginal()
return super().setupApprox(*args, **kwargs)
diff --git a/rrompy/sampling/base/sampling_engine_base_pivoted.py b/rrompy/sampling/base/sampling_engine_base_pivoted.py
index d6a4631..9c1a432 100644
--- a/rrompy/sampling/base/sampling_engine_base_pivoted.py
+++ b/rrompy/sampling/base/sampling_engine_base_pivoted.py
@@ -1,261 +1,259 @@
# Copyright (C) 2018 by the RROMPy authors
#
# This file is part of RROMPy.
#
# RROMPy is free software: you can redistribute it and/or modify
# it under the terms of the GNU Lesser General Public License as published by
# the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
#
# RROMPy is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU Lesser General Public License for more details.
#
# You should have received a copy of the GNU Lesser General Public License
# along with RROMPy. If not, see .
#
from abc import abstractmethod
import numpy as np
from rrompy.utilities.base.types import (Np1D, HFEng, List, ListAny, paramVal,
paramList, sampList, Tuple, FigHandle)
from rrompy.utilities.base import verbosityManager as vbMng, freepar as fp
from rrompy.utilities.exception_manager import RROMPyWarning
from rrompy.parameter import (emptyParameterList, checkParameter,
checkParameterList)
from rrompy.sampling import emptySampleList
from .sampling_engine_base import sampleList, SamplingEngineBase
__all__ = ['SamplingEngineBasePivoted']
class SamplingEngineBasePivoted(SamplingEngineBase):
def __init__(self, HFEngine:HFEng, directionPivot:ListAny,
*args, **kwargs):
self.directionPivot = directionPivot
self.HFEngineMarginalized = None
super().__init__(HFEngine, *args, **kwargs)
@property
def directionMarginal(self):
return tuple([x for x in range(self.HFEngine.npar) \
if x not in self.directionPivot])
@property
def nPivot(self):
return len(self.directionPivot)
@property
def nMarginal(self):
return len(self.directionMarginal)
@property
- def nsamplesTot(self):
+ def nsamplesCoalesced(self):
return np.sum(self.nsamples)
+ @property
+ def musCoalesced(self):
+ musC = emptyParameterList()
+ for j, mPs in enumerate(self.mus):
+ muEff = [fp] * self.HFEngine.npar
+ for k, x in enumerate(self.directionMarginal):
+ muEff[x] = self.musMarginal(j, k)
+ for l in range(len(mPs)):
+ for k, x in enumerate(self.directionPivot):
+ muEff[x] = mPs(l, k)
+ musC.append(muEff)
+ return musC
+
+ @property
+ def samplesCoalesced(self):
+ samplesC = emptySampleList()
+ samplesC.reset((self.samples[0].shape[0], self.nsamplesCoalesced),
+ self.samples[0].dtype)
+ run_idx = 0
+ for samp in self.samples:
+ slen = samp.shape[1]
+ samplesC.data[:, run_idx : run_idx + slen] = samp.data
+ run_idx += slen
+ return samplesC
+
def resetHistory(self, j : int = 0):
self.samples = [emptySampleList() for _ in range(j)]
self.nsamples = [0] * j
self.mus = [emptyParameterList() for _ in range(j)]
self.musMarginal = emptyParameterList()
self.musMarginal.reset((j, self.nMarginal))
self._derIdxs = [[] for _ in range(j)]
- self.resetHistoryCoalesced()
-
- def resetHistoryCoalesced(self):
- self.samplesCoalesced = emptySampleList()
def setsample(self, u:sampList, j:int, overwrite : bool = False) -> Np1D:
if overwrite:
self.samples[j][self.nsamples[j]] = u
else:
if self.nsamples[j] == 0:
self.samples[j] = sampleList(u)
else:
self.samples[j].append(u)
def popSample(self, j:int):
if hasattr(self, "nsamples") and self.nsamples[j] > 1:
if self.samples[j].shape[1] > self.nsamples[j]:
RROMPyWarning(("More than 'nsamples' memory allocated for "
"samples. Popping empty sample column."))
self.nsamples[j] += 1
self.nsamples[j] -= 1
self.samples[j].pop()
self.mus[j].pop()
self.resetHistoryCoalesced()
else:
self.resetHistory()
def preallocateSamples(self, u:sampList, mu:paramVal, n:int, j:int):
self.samples[j].reset((u.shape[0], n), u.dtype)
self.samples[j][0] = u
mu = checkParameter(mu, self.nPivot)
self.mus[j].reset((n, self.nPivot))
self.mus[j][0] = mu[0]
- def coalesceSamples(self):
- vbMng(self, "INIT", "Coalescing samples.", 7)
- self.musCoalesced = emptyParameterList()
- for j, mPs in enumerate(self.mus):
- muEff = [fp] * self.HFEngine.npar
- for k, x in enumerate(self.directionMarginal):
- muEff[x] = self.musMarginal(j, k)
- for l in range(len(mPs)):
- for k, x in enumerate(self.directionPivot):
- muEff[x] = mPs(l, k)
- self.musCoalesced.append(muEff)
- self.nsamplesCoalesced = np.sum(self.nsamples)
- self.samplesCoalesced = emptySampleList()
- self.samplesCoalesced.reset((self.samples[0].shape[0],
- self.nsamplesCoalesced),
- self.samples[0].dtype)
- run_idx = 0
- for samp in self.samples:
- slen = samp.shape[1]
- self.samplesCoalesced.data[:, run_idx : run_idx + slen] = samp.data
- run_idx += slen
- vbMng(self, "DEL", "Done coalescing samples.", 7)
-
def solveLS(self, mu : paramList = [], RHS : sampList = None) -> sampList:
"""
Solve linear system.
Args:
mu: Parameter value.
Returns:
Solution of system.
"""
mu = checkParameterList(mu, self.nPivot)[0]
vbMng(self, "INIT",
("Solving HF model for muPivot = {} and muMarginal = "
"{}.").format(mu, self.HFEngineMarginalized.muFixed), 15)
u = self.HFEngineMarginalized.solve(mu, RHS,
return_state = self.sample_state,
verbose = (self.verbosity >= 20))
vbMng(self, "DEL", "Done solving HF model.", 15)
return u
@abstractmethod
def nextSample(self, mu:paramVal, j:int, overwrite : bool = False,
postprocess : bool = True) -> Np1D:
pass
@abstractmethod
def iterSample(self, mus:paramList, musM:paramList) -> sampList:
pass
def plotSamples(self, warpings : List[List[List[callable]]] = None,
name : str = "u",
**kwargs) -> Tuple[List[List[FigHandle]], List[List[str]]]:
"""
Do some nice plots of the samples.
Args:
warpings(optional): Domain warping functions.
name(optional): Name to be shown as title of the plots. Defaults to
'u'.
Returns:
Output filenames and figure handles.
"""
if warpings is None: warpings = [[None] * self.nsamples[i] \
for i in range(len(self.nsamples))]
figs = []
filesOut = []
for i in range(len(self.nsamples)):
figsi = [None] * self.nsamples[i]
filesOuti = [None] * self.nsamples[i]
for j in range(self.nsamples[i]):
pltOut = self.HFEngine.plot(self.samples[i][j], warpings[i][j],
self.sample_state,
"{}_{}_{}".format(name, i, j),
**kwargs)
if isinstance(pltOut, (tuple,)):
figsi[j], filesOuti[j] = pltOut
else:
figsi[j] = pltOut
figs += [figsi]
filesOut += [filesOuti]
if filesOut[0][0] is None: return figs
return figs, filesOut
def outParaviewSamples(self, warpings : List[List[List[callable]]] = None,
name : str = "u", filename : str = "out",
times : List[Np1D] = None,
**kwargs) -> List[List[str]]:
"""
Output samples to ParaView file.
Args:
warpings(optional): Domain warping functions.
name(optional): Base name to be used for data output.
filename(optional): Name of output file.
times(optional): Timestamps.
Returns:
Output filenames.
"""
if warpings is None: warpings = [[None] * self.nsamples[i] \
for i in range(len(self.nsamples))]
if times is None: times = [[0.] * self.nsamples[i] \
for i in range(len(self.nsamples))]
filesOut = []
for i in range(len(self.nsamples)):
filesOuti = [None] * self.nsamples[i]
for j in range(self.nsamples[i]):
filesOuti[j] = self.HFEngine.outParaview(self.samples[i][j],
warpings[i][j], self.sample_state,
"{}_{}_{}".format(name, i, j),
"{}_{}_{}".format(filename, i, j),
times[i][j], **kwargs)
filesOut += [filesOuti]
if filesOut[0][0] is None: return None
return filesOut
def outParaviewTimeDomainSamples(self, omegas : Np1D = None,
warpings : List[List[List[callable]]] = None,
timeFinal : Np1D = None,
periodResolution : List[List[int]] = 20,
name : str = "u", filename : str = "out",
**kwargs) -> List[List[str]]:
"""
Output samples to ParaView file, converted to time domain.
Args:
omegas(optional): frequencies.
warpings(optional): Domain warping functions.
timeFinal(optional): final time of simulation.
periodResolution(optional): number of time steps per period.
name(optional): Base name to be used for data output.
filename(optional): Name of output file.
Returns:
Output filenames.
"""
if omegas is None: omegas = [np.real(self.mus[i]) \
for i in range(len(self.nsamples))]
if warpings is None: warpings = [[None] * self.nsamples[i] \
for i in range(len(self.nsamples))]
if not isinstance(timeFinal, (list, tuple,)):
timeFinal = [[timeFinal] * self.nsamples[i] \
for i in range(len(self.nsamples))]
if not isinstance(periodResolution, (list, tuple,)):
periodResolution = [[periodResolution] * self.nsamples[i] \
for i in range(len(self.nsamples))]
filesOut = []
for i in range(len(self.nsamples)):
filesOuti = [None] * self.nsamples[i]
for j in range(self.nsamples[i]):
filesOuti[j] = self.HFEngine.outParaviewTimeDomain(
self.samples[i][j], omegas[i][j],
warpings[i][j], self.sample_state,
timeFinal[i][j],
periodResolution[i][j],
"{}_{}_{}".format(name, i, j),
"{}_{}_{}".format(filename, i, j),
**kwargs)
filesOut += [filesOuti]
if filesOut[0][0] is None: return None
return filesOut
diff --git a/rrompy/sampling/pivoted/sampling_engine_pivoted_pod.py b/rrompy/sampling/pivoted/sampling_engine_pivoted_pod.py
index 4f9961e..7daa927 100644
--- a/rrompy/sampling/pivoted/sampling_engine_pivoted_pod.py
+++ b/rrompy/sampling/pivoted/sampling_engine_pivoted_pod.py
@@ -1,116 +1,104 @@
# 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.linalg import block_diag
from rrompy.sampling.base.pod_engine import PODEngine
from .sampling_engine_pivoted import SamplingEnginePivoted
from rrompy.utilities.base.types import Np1D, paramVal, sampList
from rrompy.utilities.base import verbosityManager as vbMng
from rrompy.sampling import sampleList, emptySampleList
__all__ = ['SamplingEnginePivotedPOD']
class SamplingEnginePivotedPOD(SamplingEnginePivoted):
@property
def HFEngine(self):
"""Value of HFEngine. Its assignment resets history."""
return self._HFEngine
@HFEngine.setter
def HFEngine(self, HFEngine):
self._HFEngine = HFEngine
self.resetHistory()
self.PODEngine = PODEngine(self._HFEngine)
+ @property
+ def samples_fullCoalesced(self):
+ samplesC = emptySampleList()
+ samplesC.reset((self.samples_full[0].shape[0], self.nsamplesCoalesced),
+ self.samples_full[0].dtype)
+ run_idx = 0
+ for samp in self.samples:
+ slen = samp.shape[1]
+ samplesC.data[:, run_idx : run_idx + slen] = samp.data
+ run_idx += slen
+ return samplesC
+
def resetHistory(self, j : int = 0):
super().resetHistory(j)
self.samples_full = [emptySampleList() for _ in range(j)]
self.RPOD = [np.zeros((0, 0), dtype = np.complex) for _ in range(j)]
- def resetHistoryCoalesced(self):
- super().resetHistoryCoalesced()
- self.RPODCoalesced = np.zeros((0, 0), dtype = np.complex)
- self.samples_fullCoalesced = emptySampleList()
-
def popSample(self, j:int):
if hasattr(self, "nsamples") and self.nsamples[j] > 1:
self.RPOD[j] = self.RPOD[j][: -1, : -1]
self.samples_full[j].pop()
super().popSample(j)
- def coalesceSamples(self):
- vbMng(self, "INIT", "Coalescing samples.", 7)
- verb = self.verbosity
- self.verbosity = 0
- super().coalesceSamples()
- self.verbosity = verb
- self.RPODCoalesced = block_diag(*(self.RPOD))
- self.samples_fullCoalesced = emptySampleList()
- self.samples_fullCoalesced.reset((self.samples_full[0].shape[0],
- self.nsamplesCoalesced),
- self.samples_full[0].dtype)
- ci = 0
- for j, samp_full in enumerate(self.samples_full):
- Rheg = samp_full.shape[1]
- self.samples_fullCoalesced.data[:, ci : ci + Rheg] = samp_full.data
- ci += Rheg
- vbMng(self, "DEL", "Done coalescing samples.", 7)
-
def preprocesssamples(self, idxs:Np1D, j:int) -> sampList:
if self.samples_full[j] is None or len(self.samples_full[j]) == 0:
return
return self.samples_full[j](idxs)
def setsample(self, u:sampList, j:int, overwrite : bool = False):
super().setsample(u, j, overwrite)
if overwrite:
self.samples_full[j][self.nsamples[j]] = u
else:
if self.nsamples[j] == 0:
self.samples_full[j] = sampleList(u)
else:
self.samples_full[j].append(u)
def postprocessu(self, u:sampList, j:int, overwrite : bool = False):
if overwrite:
self.samples_full[j][self.nsamples[j]] = u
else:
if self.nsamples[j] == 0:
self.samples_full[j] = sampleList(u)
else:
self.samples_full[j].append(u)
vbMng(self, "INIT", "Starting orthogonalization.", 20)
u, r, _ = self.PODEngine.GS(u, self.samples[j],
is_state = self.sample_state)
self.RPOD[j] = np.pad(self.RPOD[j], ((0, 1), (0, 1)), 'constant')
self.RPOD[j][:, -1] = r
vbMng(self, "DEL", "Done orthogonalizing.", 20)
super().setsample(u, j, overwrite)
def postprocessuBulk(self, j:int):
vbMng(self, "INIT",
"Starting orthogonalization for marginal no {}.".format(j), 40)
u, self.RPOD[j] = self.PODEngine.generalizedQR(self.samples_full[j],
is_state = self.sample_state)
vbMng(self, "DEL", "Done orthogonalizing.", 40)
self.samples[j] = sampleList(u)
def preallocateSamples(self, u:Np1D, mu:paramVal, n:int, j:int):
super().preallocateSamples(u, mu, n, j)
self.samples_full[j].reset((u.shape[0], n), u.dtype)
self.samples_full[j][0] = u