diff --git a/examples/3d/pod/plot_zero_set_3.py b/examples/3d/pod/plot_zero_set_3.py
index d3576f2..19b2856 100644
--- a/examples/3d/pod/plot_zero_set_3.py
+++ b/examples/3d/pod/plot_zero_set_3.py
@@ -1,167 +1,167 @@
import warnings
from copy import deepcopy as copy
import numpy as np
from matplotlib import pyplot as plt
def plotZeroSet3(murange, murangeEff, approx, mu0, nSamples = 200,
marginal = [None, None, .05], exps = [2., 1., 1.],
clip = -1, imagTol = 1e-8):
mNone = [i for i, m in enumerate(marginal) if m is None]
nNone = len(mNone)
if nNone not in [1, 2, 3]:
raise
if nNone == 1:
exp = exps[mNone[0]]
if hasattr(approx, "mus"):
mu2x = approx.mus(mNone[0]) ** exp
else:
mu2x = mu0[0] ** exp
murangeExp = [[murange[0][mNone[0]] ** exp],
[murange[1][mNone[0]] ** exp]]
mu1 = np.linspace(murangeEff[0][mNone[0]] ** exp,
murangeEff[1][mNone[0]] ** exp, nSamples)
mu1s = np.power(mu1, 1. / exp)
mus = []
mBase = copy(marginal)
for m1 in mu1s:
mBase[mNone[0]] = m1
mus += [copy(mBase)]
mu1 = np.real(mu1)
Z = approx.trainedModel.getQVal(mus)
Zabs = np.abs(Z)
Zmin, Zmax = np.min(Zabs), np.max(Zabs)
plt.figure(figsize = (15, 7))
plt.jet()
plt.semilogy(mu1, Zabs)
for l_ in approx.trainedModel.getPoles(marginal):
plt.plot([np.real(l_ ** exp)] * 2, [Zmin, Zmax], 'b--')
plt.plot(mu2x, [Zmax] * len(mu2x), 'kx')
plt.plot([murangeExp[0][0]] * 2, [Zmin, Zmax], 'm:')
plt.plot([murangeExp[1][0]] * 2, [Zmin, Zmax], 'm:')
plt.xlim(mu1[0], mu1[-1])
plt.title("|Q(mu)|")
plt.grid()
plt.show()
return mus, Z
if nNone == 2:
exps = [exps[m] for m in mNone]
if hasattr(approx, "mus"):
mu2x = approx.mus(mNone[0]) ** exps[0]
mu2y = approx.mus(mNone[1]) ** exps[1]
else:
mu2x, mu2y = mu0[mNone[0]] ** exps[0], mu0[mNone[1]] ** exps[1]
murangeExp = [[murange[0][mNone[0]] ** exps[0],
murange[0][mNone[1]] ** exps[1]],
[murange[1][mNone[0]] ** exps[0],
murange[1][mNone[1]] ** exps[1]]]
mu1 = np.linspace(murangeEff[0][mNone[0]] ** exps[0],
murangeEff[1][mNone[0]] ** exps[0], nSamples)
mu2 = np.linspace(murangeEff[0][mNone[1]] ** exps[1],
murangeEff[1][mNone[1]] ** exps[1], nSamples)
mu1s = np.power(mu1, 1. / exps[0])
mu2s = np.power(mu2, 1. / exps[1])
Mu1, Mu2 = np.meshgrid(np.real(mu1), np.real(mu2))
mus = []
mBase = copy(marginal)
for m2 in mu2s:
mBase[mNone[1]] = m2
for m1 in mu1s:
mBase[mNone[0]] = m1
mus += [copy(mBase)]
Z = approx.trainedModel.getQVal(mus).reshape(Mu1.shape)
Zabs = np.log10(np.abs(Z))
Zabsmin, Zabsmax = np.min(Zabs), np.max(Zabs)
if clip > 0:
Zabsmin += clip * (Zabsmax - Zabsmin)
cmap = plt.cm.bone_r
else:
cmap = plt.cm.jet
warnings.simplefilter("ignore", category = UserWarning)
plt.figure(figsize = (15, 7))
plt.jet()
p = plt.contourf(Mu1, Mu2, Zabs, cmap = cmap,
levels = np.linspace(Zabsmin, Zabsmax, 50))
if clip > 0:
plt.contour(Mu1, Mu2, Zabs, [Zabsmin])
plt.contour(Mu1, Mu2, np.real(Z), [0.], linestyles = 'dashed')
plt.contour(Mu1, Mu2, np.imag(Z), [0.], linewidths = 1,
linestyles = 'dotted')
plt.plot(mu2x, mu2y, 'kx')
plt.plot([murangeExp[0][0]] * 2,
[murangeExp[0][1], murangeExp[1][1]], 'm:')
plt.plot([murangeExp[0][0], murangeExp[1][0]],
[murangeExp[1][1]] * 2, 'm:')
plt.plot([murangeExp[1][0]] * 2,
[murangeExp[1][1], murangeExp[0][1]], 'm:')
plt.plot([murangeExp[1][0], murangeExp[0][0]],
[murangeExp[0][1]] * 2, 'm:')
plt.colorbar(p)
plt.title("log10|Q(mu)|")
plt.grid()
plt.show()
return mus, Z
if nNone == 3:
exps = [exps[m] for m in mNone]
if hasattr(approx, "mus"):
mu2x = approx.mus(mNone[0]) ** exps[0]
mu2y = approx.mus(mNone[1]) ** exps[1]
mu2z = approx.mus(mNone[2]) ** exps[2]
else:
mu2x = mu0[mNone[0]] ** exps[0]
mu2y = mu0[mNone[1]] ** exps[1]
mu2z = mu0[mNone[2]] ** exps[2]
murangeExp = [[murange[0][mNone[0]] ** exps[0],
murange[0][mNone[1]] ** exps[1],
murange[0][mNone[2]] ** exps[2]],
[murange[1][mNone[0]] ** exps[0],
murange[1][mNone[1]] ** exps[1],
murange[1][mNone[2]] ** exps[2]]]
mu1 = np.linspace(murangeEff[0][mNone[0]] ** exps[0],
murangeEff[1][mNone[0]] ** exps[0], nSamples)
mu2 = np.linspace(murangeEff[0][mNone[1]] ** exps[1],
murangeEff[1][mNone[1]] ** exps[1], nSamples)
mu3 = np.linspace(murangeEff[0][mNone[2]] ** exps[2],
murangeEff[1][mNone[2]] ** exps[2], nSamples)
mu1s = np.power(mu1, 1. / exps[0])
mu2s = np.power(mu2, 1. / exps[1])
mu3s = np.power(mu3, 1. / exps[2])
Mu1, Mu2, Mu3 = np.meshgrid(np.real(mu1), np.real(mu2), np.real(mu3),
indexing = 'ij')
mus = []
mBase = copy(marginal)
for m1 in mu1s:
mBase[mNone[0]] = m1
for m2 in mu2s:
mBase[mNone[1]] = m2
for m3 in mu3s:
mBase[mNone[2]] = m3
mus += [copy(mBase)]
Z = approx.trainedModel.getQVal(mus).reshape(Mu1.shape)
pts = []
for jdir in range(len(Mu3)):
z = Mu3[0, 0, jdir]
xis, yis, fis = Mu1[:, :, jdir], Mu2[:, :, jdir], Z[:, :, jdir]
plt.figure()
CS = plt.contour(xis, yis, np.real(fis), levels = [0.])
plt.close()
for i, CSc in enumerate(CS.allsegs):
if np.isclose(CS.cvalues[i], 0.):
for j, CScj in enumerate(CSc):
for k in range(len(CScj)):
pts += [[CScj[k, 0], CScj[k, 1], z]]
pts += [[np.nan] * 3]
pts = np.array(pts)
if np.size(pts) > 0:
from mpl_toolkits.mplot3d import Axes3D
fig = plt.figure(figsize = (8, 6))
ax = Axes3D(fig)
ax.plot(pts[:, 0], pts[:, 1], pts[:, 2], 'k-')
- ax.scatter(mu2x, mu2y, mu2z, '.')
+ ax.scatter(np.real(mu2x), np.real(mu2y), np.real(mu2z), '.')
ax.set_xlim(murangeExp[0][0], murangeExp[1][0])
ax.set_ylim(murangeExp[0][1], murangeExp[1][1])
ax.set_zlim(murangeExp[0][2], murangeExp[1][2])
plt.show()
plt.close()
- return pts
\ No newline at end of file
+ return pts
diff --git a/rrompy/reduction_methods/pivoting/generic_pivoted_approximant.py b/rrompy/reduction_methods/pivoting/generic_pivoted_approximant.py
index 03039fd..d9c2cc6 100644
--- a/rrompy/reduction_methods/pivoting/generic_pivoted_approximant.py
+++ b/rrompy/reduction_methods/pivoting/generic_pivoted_approximant.py
@@ -1,484 +1,481 @@
# 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 copy import deepcopy as copy
from rrompy.reduction_methods.base.generic_approximant import (
GenericApproximant)
from rrompy.hfengines.base import MarginalProxyEngine
from rrompy.utilities.poly_fitting.polynomial import (polybases,
nextDerivativeIndices,
PolynomialInterpolator as PI)
from rrompy.utilities.poly_fitting.radial_basis import (
RadialBasisInterpolator as RBI)
from rrompy.utilities.poly_fitting.radial_basis import rbbases
from rrompy.sampling.pivoted import (SamplingEnginePivoted,
SamplingEnginePivotedPOD)
from rrompy.utilities.base.types import (ListAny, DictAny, HFEng, paramVal,
paramList)
from rrompy.utilities.base import verbosityManager as vbMng, freepar as fp
from rrompy.utilities.exception_manager import (RROMPyException, RROMPyAssert,
RROMPyWarning)
from rrompy.parameter import emptyParameterList
__all__ = ['GenericPivotedApproximant']
class GenericPivotedApproximant(GenericApproximant):
"""
ROM pivoted approximant computation for parametric problems (ABSTRACT).
Args:
HFEngine: HF problem solver.
mu0(optional): Default parameter. Defaults to 0.
directionPivot(optional): Pivot components. Defaults to [0].
approxParameters(optional): Dictionary containing values for main
parameters of approximant. Recognized keys are:
- 'POD': whether to compute POD of snapshots; defaults to True;
- '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;
- 'polybasisMarginal': type of polynomial basis for marginal
interpolation; allowed values include 'MONOMIAL', 'CHEBYSHEV'
and 'LEGENDRE'; defaults to 'MONOMIAL';
- 'MMarginal': degree of marginal interpolant; defaults to 0;
- 'radialBasisMarginal': radial basis family for marginal
interpolant; defaults to 0, i.e. no radial basis;
- 'radialBasisWeightsMarginal': radial basis weights for marginal
interpolant; defaults to 0, i.e. identity;
- 'interpRcondMarginal': tolerance for marginal interpolation;
defaults to None.
Defaults to empty dict.
homogeneized(optional): Whether to homogeneize Dirichlet BCs. Defaults
to False.
verbosity(optional): Verbosity level. Defaults to 10.
Attributes:
HFEngine: HF problem solver.
mu0: Default parameter.
directionPivot: Pivot components.
mus: Array of snapshot parameters.
musPivot: Array of pivot snapshot parameters.
musMarginal: Array of marginal snapshot parameters.
homogeneized: Whether to homogeneize Dirichlet BCs.
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;
- 'polybasisMarginal': type of polynomial basis for marginal
interpolation;
- 'MMarginal': degree of marginal interpolant;
- 'radialBasisMarginal': radial basis family for marginal
interpolant;
- 'radialBasisWeightsMarginal': radial basis weights for marginal
interpolant;
- 'interpRcondMarginal': tolerance for marginal interpolation.
parameterListCritical: Recognized keys of critical approximant
parameters:
- 'S': total number of pivot samples current approximant relies
upon;
- 'samplerPivot': pivot sample point generator;
- 'SMarginal': total number of marginal samples current approximant
relies upon;
- 'samplerMarginal': marginal sample point generator.
POD: Whether to compute POD of snapshots.
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.
polybasisMarginal: Type of polynomial basis for marginal interpolation.
MMarginal: Degree of marginal interpolant.
radialBasisMarginal: Radial basis family for marginal interpolant.
radialBasisWeightsMarginal: Radial basis weights for marginal
interpolant.
interpRcondMarginal: Tolerance for marginal interpolation.
muBoundsPivot: list of bounds for pivot parameter values.
muBoundsMarginal: list of bounds for marginal parameter values.
samplingEngine: Sampling engine.
uHF: High fidelity solution(s) with parameter(s) lastSolvedHF as
sampleList.
lastSolvedHF: Parameter(s) corresponding to last computed high fidelity
solution(s) as parameterList.
uApproxReduced: Reduced approximate solution(s) with parameter(s)
lastSolvedApprox as sampleList.
lastSolvedApproxReduced: Parameter(s) corresponding to last computed
reduced approximate solution(s) as parameterList.
uApprox: Approximate solution(s) with parameter(s) lastSolvedApprox as
sampleList.
lastSolvedApprox: Parameter(s) corresponding to last computed
approximate solution(s) as parameterList.
"""
def __init__(self, HFEngine:HFEng, mu0 : paramVal = None,
directionPivot : ListAny = [0],
approxParameters : DictAny = {}, homogeneized : bool = False,
verbosity : int = 10, timestamp : bool = True):
self._preInit()
from rrompy.parameter.parameter_sampling import QuadratureSampler as QS
QSBase = QS([[0], [1]], "UNIFORM")
self._addParametersToList(
["polybasisMarginal", "MMarginal", "radialBasisMarginal",
"radialBasisWeightsMarginal", "interpRcondMarginal"],
["MONOMIAL", 0, 0, 1, -1],
["samplerPivot", "SMarginal", "samplerMarginal"],
[QSBase, [1], QSBase])
del QS
self._directionPivot = directionPivot
super().__init__(HFEngine = HFEngine, mu0 = mu0,
approxParameters = approxParameters,
homogeneized = homogeneized,
verbosity = verbosity, timestamp = timestamp)
self._postInit()
def setupSampling(self):
"""Setup sampling engine."""
RROMPyAssert(self._mode, message = "Cannot setup sampling engine.")
if not hasattr(self, "_POD") or self._POD is None: return
if self.POD:
SamplingEngine = SamplingEnginePivotedPOD
else:
SamplingEngine = SamplingEnginePivoted
self.samplingEngine = SamplingEngine(self.HFEngine,
self.directionPivot,
verbosity = self.verbosity,
allowRepeatedSamples = True)
@property
def SMarginal(self):
"""Value of SMarginal."""
return self._SMarginal
@SMarginal.setter
def SMarginal(self, SMarginal):
if not hasattr(SMarginal, "__len__"): SMarginal = [SMarginal]
if any([s <= 0 for s in SMarginal]):
raise RROMPyException("SMarginal must be positive.")
if hasattr(self, "_SMarginal") and self._SMarginal is not None:
Sold = tuple(self.SMarginal)
else: Sold = -1
self._SMarginal = SMarginal
self._approxParameters["SMarginal"] = self.SMarginal
if Sold != tuple(self.SMarginal):
self.resetSamples()
@property
def polybasisMarginal(self):
"""Value of polybasisMarginal."""
return self._polybasisMarginal
@polybasisMarginal.setter
def polybasisMarginal(self, polybasisMarginal):
try:
polybasisMarginal = polybasisMarginal.upper().strip().replace(" ",
"")
if polybasisMarginal not in polybases:
raise RROMPyException(
"Prescribed marginal polybasis not recognized.")
self._polybasisMarginal = polybasisMarginal
except:
RROMPyWarning(("Prescribed marginal polybasis not recognized. "
"Overriding to 'MONOMIAL'."))
self._polybasisMarginal = "MONOMIAL"
self._approxParameters["polybasisMarginal"] = self.polybasisMarginal
@property
def MMarginal(self):
"""Value of MMarginal."""
return self._MMarginal
@MMarginal.setter
def MMarginal(self, MMarginal):
if MMarginal < 0:
raise RROMPyException("MMarginal must be non-negative.")
self._MMarginal = MMarginal
self._approxParameters["MMarginal"] = self.MMarginal
@property
def radialBasisMarginal(self):
"""Value of radialBasisMarginal."""
return self._radialBasisMarginal
@radialBasisMarginal.setter
def radialBasisMarginal(self, radialBasisMarginal):
try:
if radialBasisMarginal != 0:
radialBasisMarginal = radialBasisMarginal.upper().strip()\
.replace(" ","")
if radialBasisMarginal not in rbbases:
raise RROMPyException(("Prescribed marginal radialBasis "
"not recognized."))
self._radialBasisMarginal = radialBasisMarginal
except:
RROMPyWarning(("Prescribed marginal radialBasis not recognized. "
"Overriding to 0."))
self._radialBasisMarginal = 0
self._approxParameters["radialBasisMarginal"] = (
self.radialBasisMarginal)
@property
def radialBasisWeightsMarginal(self):
"""Value of radialBasisWeightsMarginal."""
return self._radialBasisWeightsMarginal
@radialBasisWeightsMarginal.setter
def radialBasisWeightsMarginal(self, radialBasisWeightsMarginal):
self._radialBasisWeightsMarginal = radialBasisWeightsMarginal
self._approxParameters["radialBasisWeightsMarginal"] = (
self.radialBasisWeightsMarginal)
@property
def polybasisMarginalP(self):
if self.radialBasisMarginal == 0:
return self._polybasisMarginal
return self._polybasisMarginal + "_" + self.radialBasisMarginal
@property
def interpRcondMarginal(self):
"""Value of interpRcondMarginal."""
return self._interpRcondMarginal
@interpRcondMarginal.setter
def interpRcondMarginal(self, interpRcondMarginal):
self._interpRcondMarginal = interpRcondMarginal
self._approxParameters["interpRcondMarginal"] = (
self.interpRcondMarginal)
@property
def directionPivot(self):
"""Value of directionPivot. Its assignment may reset snapshots."""
return self._directionPivot
@directionPivot.setter
def directionPivot(self, directionPivot):
if hasattr(self, '_directionPivot'):
directionPivotOld = copy(self.directionPivot)
else:
directionPivotOld = None
if (directionPivotOld is None
or len(directionPivot) != len(directionPivotOld)
or not directionPivot == directionPivotOld):
self.resetSamples()
self._directionPivot = directionPivot
@property
def directionMarginal(self):
return [x for x in range(self.HFEngine.npar) \
if x not in self.directionPivot]
@property
def rescalingExpPivot(self):
return [self.HFEngine.rescalingExp[x] for x in self.directionPivot]
@property
def rescalingExpMarginal(self):
return [self.HFEngine.rescalingExp[x] for x in self.directionMarginal]
@property
def muBoundsPivot(self):
"""Value of muBoundsPivot."""
return self.samplerPivot.lims
@property
def muBoundsMarginal(self):
"""Value of muBoundsMarginal."""
return self.samplerMarginal.lims
@property
def samplerPivot(self):
"""Value of samplerPivot."""
return self._samplerPivot
@samplerPivot.setter
def samplerPivot(self, samplerPivot):
if 'generatePoints' not in dir(samplerPivot):
raise RROMPyException("Pivot sampler type not recognized.")
if hasattr(self, '_samplerPivot') and self._samplerPivot is not None:
samplerOld = self.samplerPivot
self._samplerPivot = samplerPivot
self._approxParameters["samplerPivot"] = self.samplerPivot.__str__()
if not 'samplerOld' in locals() or samplerOld != self.samplerPivot:
self.resetSamples()
@property
def samplerMarginal(self):
"""Value of samplerMarginal."""
return self._samplerMarginal
@samplerMarginal.setter
def samplerMarginal(self, samplerMarginal):
if 'generatePoints' not in dir(samplerMarginal):
raise RROMPyException("Marginal sampler type not recognized.")
if (hasattr(self, '_samplerMarginal')
and self._samplerMarginal is not None):
samplerOld = self.samplerMarginal
self._samplerMarginal = samplerMarginal
self._approxParameters["samplerMarginal"] = (
self.samplerMarginal.__str__())
if not 'samplerOld' in locals() or samplerOld != self.samplerMarginal:
self.resetSamples()
def resetSamples(self):
"""Reset samples."""
super().resetSamples()
self._musMUniqueCN = None
self._derMIdxs = None
self._reorderM = None
def setSamples(self, samplingEngine):
"""Copy samplingEngine and samples."""
super().setSamples(samplingEngine)
self.mus = copy(self.samplingEngine[0].mus)
for sEj in self.samplingEngine[1:]:
self.mus.append(sEj.mus)
def computeSnapshots(self):
"""Compute snapshots of solution map."""
RROMPyAssert(self._mode,
message = "Cannot start snapshot computation.")
Sp = np.prod(self.S)
Seff = Sp * np.prod(self.SMarginal)
if self.samplingEngine.nsamplesTot != Seff:
self.resetSamples()
vbMng(self, "INIT", "Starting computation of snapshots.", 5)
+ if self.HFEngine.As[0] is None: self.HFEngine.A(self.mu0)
+ if self.HFEngine.bs[0] is None: self.HFEngine.b(self.mu0)
self.musPivot = self.samplerPivot.generatePoints(self.S)
self.musMarginal = self.samplerMarginal.generatePoints(
self.SMarginal)
self.mus = emptyParameterList()
self.mus.reset((Seff, self.HFEngine.npar))
self.samplingEngine.resetHistory(Seff // Sp)
for j, muMarg in enumerate(self.musMarginal):
for k in range(j * Sp, (j + 1) * Sp):
self.mus.data[k, self.directionPivot] = (
- self.musPivot[k - j*Sp].data)
+ self.musPivot[k - j * Sp].data)
self.mus.data[k, self.directionMarginal] = muMarg.data
- muMEff = [fp] * self.npar
- for k, x in enumerate(self.directionMarginal):
- muMEff[x] = muMarg(0, k)
- self.samplingEngine.HFEngineMarginalized = MarginalProxyEngine(
- self.HFEngine, list(muMEff))
- self.samplingEngine.iterSample(self.musPivot, j,
- homogeneized = self.homogeneized)
+ self.samplingEngine.iterSample(self.musPivot, self.musMarginal,
+ homogeneized = self.homogeneized)
if self.POD:
self.samplingEngine.coalesceSamples(self.interpRcondMarginal)
else:
self.samplingEngine.coalesceSamples()
vbMng(self, "DEL", "Done computing snapshots.", 5)
def _setupMarginalInterpolationIndices(self):
"""Setup parameters for polyvander."""
RROMPyAssert(self._mode,
message = "Cannot setup interpolation indices.")
if (self._musMUniqueCN is None
or len(self._reorderM) != len(self.musMarginal)):
self._musMUniqueCN, musMIdxsTo, musMIdxs, musMCount = (
self.centerNormalizeMarginal(self.musMarginal).unique(
return_index = True,
return_inverse = True,
return_counts = True))
self._musMUnique = self.musMarginal[musMIdxsTo]
self._derMIdxs = [None] * len(self._musMUniqueCN)
self._reorderM = np.empty(len(musMIdxs), dtype = int)
filled = 0
for j, cnt in enumerate(musMCount):
self._derMIdxs[j] = nextDerivativeIndices([],
self.musMarginal.shape[1], cnt)
jIdx = np.nonzero(musMIdxs == j)[0]
self._reorderM[jIdx] = np.arange(filled, filled + cnt)
filled += cnt
def _setupMarginalInterp(self):
"""Compute marginal interpolator."""
RROMPyAssert(self._mode, message = "Cannot setup numerator.")
vbMng(self, "INIT", "Starting computation of marginal interpolator.",
7)
self._setupMarginalInterpolationIndices()
MMarginal0 = copy(self.MMarginal)
mI = []
for j in range(len(self.musMarginal)):
canonicalj = 1. * (np.arange(len(self.musMarginal)) == j)
self._MMarginal = MMarginal0
while self.MMarginal >= 0:
if self.radialBasisMarginal == 0:
p = PI()
wellCond, msg = p.setupByInterpolation(
self._musMUniqueCN, canonicalj, self.MMarginal,
self.polybasisMarginal, self.verbosity >= 5, True,
{"derIdxs": self._derMIdxs,
"reorder": self._reorderM,
"scl": np.power(self.scaleFactorMarginal, -1.)},
{"rcond": self.interpRcondMarginal})
else:
p = RBI()
wellCond, msg = p.setupByInterpolation(
self._musMUniqueCN, canonicalj, self.MMarginal,
self.polybasisMarginalP,
self.radialBasisWeightsMarginal,
self.verbosity >= 5, True,
{"derIdxs": self._derMIdxs,
"reorder": self._reorderM,
"scl": np.power(self.scaleFactorMarginal, -1.)},
{"rcond": self.interpRcondMarginal})
vbMng(self, "MAIN", msg, 5)
if wellCond: break
RROMPyWarning(("Polyfit is poorly conditioned. Reducing "
"MMarginal by 1."))
self.MMarginal -= 1
if self.MMarginal < 0:
raise RROMPyException(("Instability in computation of "
"marginal interpolator. Aborting."))
mI = mI + [copy(p)]
vbMng(self, "DEL", "Done computing marginal interpolator.", 7)
return mI
def normApprox(self, mu:paramList, homogeneized : bool = False) -> float:
"""
Compute norm of approximant at arbitrary parameter.
Args:
mu: Target parameter.
homogeneized(optional): Whether to remove Dirichlet BC. Defaults to
False.
Returns:
Target norm of approximant.
"""
if not self.POD or self.homogeneized != homogeneized:
return super().normApprox(mu, homogeneized)
return np.linalg.norm(self.getApproxReduced(mu).data, axis = 0)
def computeScaleFactor(self):
"""Compute parameter rescaling factor."""
RROMPyAssert(self._mode, message = "Cannot compute rescaling factor.")
self.scaleFactorPivot = .5 * np.abs(
self.muBoundsPivot[0] ** self.rescalingExpPivot
- self.muBoundsPivot[1] ** self.rescalingExpPivot)
self.scaleFactorMarginal = .5 * np.abs(
self.muBoundsMarginal[0] ** self.rescalingExpMarginal
- self.muBoundsMarginal[1] ** self.rescalingExpMarginal)
self.scaleFactor = np.empty(self.npar)
self.scaleFactor[self.directionPivot] = self.scaleFactorPivot
self.scaleFactor[self.directionMarginal] = self.scaleFactorMarginal
def centerNormalizeMarginal(self, mu : paramList = [],
mu0 : paramVal = None) -> paramList:
"""
Compute normalized parameter to be plugged into approximant.
Args:
mu: Parameter(s) 1.
mu0: Parameter(s) 2. If None, set to self.mu0.
Returns:
Normalized parameter.
"""
return self.trainedModel.centerNormalizeMarginal(mu, mu0)
diff --git a/rrompy/reduction_methods/pivoting/rational_interpolant_pivoted.py b/rrompy/reduction_methods/pivoting/rational_interpolant_pivoted.py
index 0916139..af5207a 100644
--- a/rrompy/reduction_methods/pivoting/rational_interpolant_pivoted.py
+++ b/rrompy/reduction_methods/pivoting/rational_interpolant_pivoted.py
@@ -1,600 +1,602 @@
# Copyright (C) 2018 by the RROMPy authors
#
# This file is part of RROMPy.
#
# RROMPy is free software: you can redistribute it and/or modify
# it under the terms of the GNU Lesser General Public License as published by
# the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
#
# RROMPy is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU Lesser General Public License for more details.
#
# You should have received a copy of the GNU Lesser General Public License
# along with RROMPy. If not, see .
#
from copy import deepcopy as copy
import numpy as np
from rrompy.reduction_methods.base import checkRobustTolerance
from .generic_pivoted_approximant import GenericPivotedApproximant
from rrompy.reduction_methods.distributed.rational_interpolant import (
RationalInterpolant as RI)
from rrompy.utilities.poly_fitting import customPInv
from rrompy.utilities.poly_fitting.polynomial import (polybases, polyfitname,
nextDerivativeIndices,
hashDerivativeToIdx as hashD,
hashIdxToDerivative as hashI,
homogeneizedpolyvander as hpvP,
homogeneizedToFull,
PolynomialInterpolator as PI)
from rrompy.utilities.poly_fitting.radial_basis import (rbbases,
RadialBasisInterpolator as RBI)
from rrompy.reduction_methods.trained_model import (TrainedModelPivotedData,
TrainedModelPivotedPade as tModel)
from rrompy.utilities.base.types import (Np1D, Np2D, HFEng, DictAny, List,
ListAny, paramVal, paramList)
from rrompy.utilities.base import (verbosityManager as vbMng, multifactorial,
freepar as fp)
from rrompy.utilities.exception_manager import (RROMPyException, RROMPyAssert,
RROMPyWarning)
from rrompy.parameter import checkParameter
__all__ = ['RationalInterpolantPivoted']
class RationalInterpolantPivoted(GenericPivotedApproximant):
"""
ROM pivoted rational interpolant 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;
- '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;
- 'polybasisPivot': type of polynomial basis for pivot
interpolation; allowed values include 'MONOMIAL', 'CHEBYSHEV'
and 'LEGENDRE'; defaults to 'MONOMIAL';
- 'polybasisMarginal': type of polynomial basis for marginal
interpolation; allowed values include 'MONOMIAL', 'CHEBYSHEV'
and 'LEGENDRE'; defaults to 'MONOMIAL';
- 'E': number of derivatives used to compute interpolant; defaults
to 0;
- 'M': degree of Pade' interpolant numerator; defaults to 0;
- 'N': degree of Pade' interpolant denominator; defaults to 0;
- 'radialBasisPivot': radial basis family for pivot numerator;
defaults to 0, i.e. no radial basis;
- 'radialBasisMarginal': radial basis family for marginal
interpolant; defaults to 0, i.e. no radial basis;
- 'radialBasisWeightsPivot': radial basis weights for pivot
numerator; defaults to 0, i.e. identity;
- 'radialBasisWeightsMarginal': radial basis weights for marginal
interpolant; defaults to 0, i.e. identity;
- 'interpRcondPivot': tolerance for pivot interpolation; defaults
to None;
- 'interpRcondMarginal': tolerance for marginal interpolation;
defaults to None;
- 'robustTol': tolerance for robust Pade' denominator management;
defaults to 0.
Defaults to empty dict.
homogeneized(optional): Whether to homogeneize Dirichlet BCs. Defaults
to False.
verbosity(optional): Verbosity level. Defaults to 10.
Attributes:
HFEngine: HF problem solver.
mu0: Default parameter.
directionPivot: Pivot components.
mus: Array of snapshot parameters.
musPivot: Array of pivot snapshot parameters.
musMarginal: Array of marginal snapshot parameters.
homogeneized: Whether to homogeneize Dirichlet BCs.
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;
- 'polybasisPivot': type of polynomial basis for pivot
interpolation;
- 'polybasisMarginal': type of polynomial basis for marginal
interpolation;
- 'E': number of derivatives used to compute interpolant;
- 'M': degree of Pade' interpolant numerator;
- 'N': degree of Pade' interpolant denominator;
- 'radialBasisPivot': radial basis family for pivot numerator;
- 'radialBasisMarginal': radial basis family for marginal
interpolant;
- 'radialBasisWeightsPivot': radial basis weights for pivot
numerator;
- 'radialBasisWeightsMarginal': radial basis weights for marginal
interpolant;
- 'interpRcondPivot': tolerance for pivot interpolation;
- 'interpRcondMarginal': tolerance for marginal interpolation;
- 'robustTol': tolerance for robust Pade' denominator management.
parameterListCritical: Recognized keys of critical approximant
parameters:
- 'S': total number of pivot samples current approximant relies
upon;
- 'samplerPivot': pivot sample point generator;
- 'SMarginal': total number of marginal samples current approximant
relies upon;
- 'samplerMarginal': marginal sample point generator.
POD: Whether to compute POD of snapshots.
S: Total number of pivot samples current approximant relies upon.
sampler: Pivot sample point generator.
polybasisPivot: Type of polynomial basis for pivot interpolation.
polybasisMarginal: Type of polynomial basis for marginal interpolation.
M: Numerator degree of approximant.
N: Denominator degree of approximant.
radialBasisPivot: Radial basis family for pivot numerator.
radialBasisMarginal: Radial basis family for marginal interpolant.
radialBasisWeightsPivot: Radial basis weights for pivot numerator.
radialBasisWeightsMarginal: Radial basis weights for marginal
interpolant.
interpRcondPivot: Tolerance for pivot interpolation.
interpRcondMarginal: Tolerance for marginal interpolation.
robustTol: Tolerance for robust Pade' denominator management.
E: Complete derivative depth level of S.
muBoundsPivot: list of bounds for pivot parameter values.
muBoundsMarginal: list of bounds for marginal parameter values.
samplingEngine: Sampling engine.
uHF: High fidelity solution(s) with parameter(s) lastSolvedHF as
sampleList.
lastSolvedHF: Parameter(s) corresponding to last computed high fidelity
solution(s) as parameterList.
uApproxReduced: Reduced approximate solution(s) with parameter(s)
lastSolvedApprox as sampleList.
lastSolvedApproxReduced: Parameter(s) corresponding to last computed
reduced approximate solution(s) as parameterList.
uApprox: Approximate solution(s) with parameter(s) lastSolvedApprox as
sampleList.
lastSolvedApprox: Parameter(s) corresponding to last computed
approximate solution(s) as parameterList.
Q: Numpy 1D vector containing complex coefficients of approximant
denominator.
P: Numpy 2D vector whose columns are FE dofs of coefficients of
approximant numerator.
"""
def __init__(self, HFEngine:HFEng, mu0 : paramVal = None,
directionPivot : ListAny = [0],
approxParameters : DictAny = {}, homogeneized : bool = False,
verbosity : int = 10, timestamp : bool = True):
self._preInit()
self._addParametersToList(
["polybasisPivot", "E", "M", "N",
"radialBasisPivot", "radialBasisWeightsPivot",
"interpRcondPivot", "robustTol"],
["MONOMIAL", -1, 0, 0, 0, 1, -1, 0])
super().__init__(HFEngine = HFEngine, mu0 = mu0,
directionPivot = directionPivot,
approxParameters = approxParameters,
homogeneized = homogeneized,
verbosity = verbosity, timestamp = timestamp)
self._postInit()
@property
def polybasisPivot(self):
"""Value of polybasisPivot."""
return self._polybasisPivot
@polybasisPivot.setter
def polybasisPivot(self, polybasisPivot):
try:
polybasisPivot = polybasisPivot.upper().strip().replace(" ","")
if polybasisPivot not in polybases:
raise RROMPyException(
"Prescribed pivot polybasis not recognized.")
self._polybasisPivot = polybasisPivot
except:
RROMPyWarning(("Prescribed pivot polybasis not recognized. "
"Overriding to 'MONOMIAL'."))
self._polybasisPivot = "MONOMIAL"
self._approxParameters["polybasisPivot"] = self.polybasisPivot
@property
def radialBasisPivot(self):
"""Value of radialBasisPivot."""
return self._radialBasisPivot
@radialBasisPivot.setter
def radialBasisPivot(self, radialBasisPivot):
try:
if radialBasisPivot != 0:
radialBasisPivot = radialBasisPivot.upper().strip().replace(
" ","")
if radialBasisPivot not in rbbases:
raise RROMPyException(("Prescribed pivot radialBasis not "
"recognized."))
self._radialBasisPivot = radialBasisPivot
except:
RROMPyWarning(("Prescribed pivot radialBasis not recognized. "
"Overriding to 0."))
self._radialBasisPivot = 0
self._approxParameters["radialBasisPivot"] = self.radialBasisPivot
@property
def radialBasisWeightsPivot(self):
"""Value of radialBasisWeightsPivot."""
return self._radialBasisWeightsPivot
@radialBasisWeightsPivot.setter
def radialBasisWeightsPivot(self, radialBasisWeightsPivot):
self._radialBasisWeightsPivot = radialBasisWeightsPivot
self._approxParameters["radialBasisWeightsPivot"] = (
self.radialBasisWeightsPivot)
@property
def polybasisPivotP(self):
if self.radialBasisPivot == 0:
return self._polybasisPivot
return self._polybasisPivot + "_" + self.radialBasisPivot
@property
def interpRcondPivot(self):
"""Value of interpRcondPivot."""
return self._interpRcondPivot
@interpRcondPivot.setter
def interpRcondPivot(self, interpRcondPivot):
self._interpRcondPivot = interpRcondPivot
self._approxParameters["interpRcondPivot"] = self.interpRcondPivot
@property
def M(self):
"""Value of M. Its assignment may change S."""
return self._M
@M.setter
def M(self, M):
if M < 0: raise RROMPyException("M must be non-negative.")
self._M = M
self._approxParameters["M"] = self.M
if hasattr(self, "_E") and self.E >= 0 and self.E < self.M:
RROMPyWarning("Prescribed S is too small. Decreasing M.")
self.M = self.E
@property
def N(self):
"""Value of N. Its assignment may change S."""
return self._N
@N.setter
def N(self, N):
if N < 0: raise RROMPyException("N must be non-negative.")
self._N = N
self._approxParameters["N"] = self.N
if hasattr(self, "_E") and self.E >= 0 and self.E < self.N:
RROMPyWarning("Prescribed S is too small. Decreasing N.")
self.N = self.E
@property
def E(self):
"""Value of E."""
return self._E
@E.setter
def E(self, E):
if E < 0:
if not hasattr(self, "_S"):
raise RROMPyException(("Value of E must be positive if S is "
"not yet assigned."))
E = np.sum(hashI(np.prod(self.S), len(self.directionPivot))) - 1
self._E = E
self._approxParameters["E"] = self.E
if (hasattr(self, "_S")
and self.E >= np.sum(hashI(np.prod(self.S),len(self.directionPivot)))):
RROMPyWarning("Prescribed S is too small. Decreasing E.")
self.E = -1
if hasattr(self, "_M"): self.M = self.M
if hasattr(self, "_N"): self.N = self.N
@property
def robustTol(self):
"""Value of tolerance for robust Pade' denominator management."""
return self._robustTol
@robustTol.setter
def robustTol(self, robustTol):
if robustTol < 0.:
RROMPyWarning(("Overriding prescribed negative robustness "
"tolerance to 0."))
robustTol = 0.
self._robustTol = robustTol
self._approxParameters["robustTol"] = self.robustTol
@property
def S(self):
"""Value of S."""
return self._S
@S.setter
def S(self, S):
GenericPivotedApproximant.S.fset(self, S)
if hasattr(self, "_M"): self.M = self.M
if hasattr(self, "_N"): self.N = self.N
if hasattr(self, "_E"): self.E = self.E
def resetSamples(self):
"""Reset samples."""
super().resetSamples()
self._musPUniqueCN = None
self._derPIdxs = None
self._reorderP = None
def _setupPivotInterpolationIndices(self):
"""Setup parameters for polyvander."""
RROMPyAssert(self._mode,
message = "Cannot setup interpolation indices.")
if (self._musPUniqueCN is None
or len(self._reorderP) != len(self.musPivot)):
self._musPUniqueCN, musPIdxsTo, musPIdxs, musPCount = (
self.centerNormalizePivot(self.musPivot).unique(
return_index = True,
return_inverse = True,
return_counts = True))
self._musPUnique = self.mus[musPIdxsTo]
self._derPIdxs = [None] * len(self._musPUniqueCN)
self._reorderP = np.empty(len(musPIdxs), dtype = int)
filled = 0
for j, cnt in enumerate(musPCount):
self._derPIdxs[j] = nextDerivativeIndices([],
self.musPivot.shape[1], cnt)
jIdx = np.nonzero(musPIdxs == j)[0]
self._reorderP[jIdx] = np.arange(filled, filled + cnt)
filled += cnt
def _setupDenominator(self):
"""Compute Pade' denominator."""
RROMPyAssert(self._mode, message = "Cannot setup denominator.")
vbMng(self, "INIT", "Starting computation of denominator.", 7)
NinvD = None
N0 = copy(self.N)
qs = []
+ self.verbosity -= 10
for j in range(len(self.musMarginal)):
self._N = N0
while self.N > 0:
if NinvD != self.N:
invD = self._computeInterpolantInverseBlocks()
NinvD = self.N
if self.POD:
ev, eV = RI.findeveVGQR(self, self.samplingEngine.RPOD[j],
invD)
else:
ev, eV = RI.findeveVGExplicit(self,
self.samplingEngine.samples[j], invD)
newParams = checkRobustTolerance(ev, self.N, self.robustTol)
if not newParams:
break
self.approxParameters = newParams
if self.N <= 0:
self._N = 0
eV = np.ones((1, 1))
q = PI()
q.npar = self.musPivot.shape[1]
q.polybasis = self.polybasisPivot
q.coeffs = homogeneizedToFull(tuple([self.N + 1] * q.npar),
q.npar, eV[:, 0])
qs = qs + [copy(q)]
+ self.verbosity += 10
vbMng(self, "DEL", "Done computing denominator.", 7)
return qs
def _setupNumerator(self):
"""Compute Pade' numerator."""
RROMPyAssert(self._mode, message = "Cannot setup numerator.")
vbMng(self, "INIT", "Starting computation of numerator.", 7)
Qevaldiag = np.zeros((len(self.musPivot), len(self.musPivot)),
dtype = np.complex)
verb = self.trainedModel.verbosity
self.trainedModel.verbosity = 0
self._setupPivotInterpolationIndices()
M0 = copy(self.M)
tensor_idx = 0
ps = []
for k, muM in enumerate(self.musMarginal):
self._M = M0
idxGlob = 0
for j, derIdxs in enumerate(self._derPIdxs):
mujEff = [fp] * self.npar
for jj, kk in enumerate(self.directionPivot):
mujEff[kk] = self._musPUnique[j, jj]
for jj, kk in enumerate(self.directionMarginal):
mujEff[kk] = muM(0, jj)
mujEff = checkParameter(mujEff, self.npar)
nder = len(derIdxs)
idxLoc = np.arange(len(self.musPivot))[
(self._reorderP >= idxGlob)
* (self._reorderP < idxGlob + nder)]
idxGlob += nder
Qval = [0] * nder
for der in range(nder):
derIdx = hashI(der, self.musPivot.shape[1])
derIdxEff = [0] * self.npar
sclEff = [0] * self.npar
for jj, kk in enumerate(self.directionPivot):
derIdxEff[kk] = derIdx[jj]
sclEff[kk] = self.scaleFactorPivot[jj] ** -1.
Qval[der] = (self.trainedModel.getQVal(mujEff, derIdxEff,
scl = sclEff)
/ multifactorial(derIdx))
for derU, derUIdx in enumerate(derIdxs):
for derQ, derQIdx in enumerate(derIdxs):
diffIdx = [x - y for (x, y) in zip(derUIdx, derQIdx)]
if all([x >= 0 for x in diffIdx]):
diffj = hashD(diffIdx)
Qevaldiag[idxLoc[derU], idxLoc[derQ]] = Qval[diffj]
- self.trainedModel.verbosity = verb
while self.M >= 0:
if self.radialBasisPivot == 0:
p = PI()
wellCond, msg = p.setupByInterpolation(
self._musPUniqueCN, Qevaldiag, self.M,
self.polybasisPivotP, self.verbosity >= 5, True,
{"derIdxs": self._derPIdxs,
"reorder": self._reorderP,
"scl": np.power(self.scaleFactorPivot, -1.)},
{"rcond": self.interpRcondPivot})
else:
p = RBI()
wellCond, msg = p.setupByInterpolation(
self._musPUniqueCN, Qevaldiag, self.M,
self.polybasisPivotP,
self.radialBasisWeightsPivot,
self.verbosity >= 5, True,
{"derIdxs": self._derPIdxs,
"reorder": self._reorderP,
"scl": np.power(self.scaleFactorPivot, -1.)},
{"rcond": self.interpRcondPivot})
vbMng(self, "MAIN", msg, 5)
if wellCond: break
RROMPyWarning(("Polyfit is poorly conditioned. "
"Reducing M by 1."))
self.M -= 1
if self.M < 0:
raise RROMPyException(("Instability in computation of "
"numerator. Aborting."))
tensor_idx_new = tensor_idx + Qevaldiag.shape[1]
if self.POD:
p.postmultiplyTensorize(self.samplingEngine.RPODCoalesced.T[
tensor_idx : tensor_idx_new, :])
else:
p.pad(tensor_idx, len(self.mus) - tensor_idx_new)
tensor_idx = tensor_idx_new
ps = ps + [copy(p)]
+ self.trainedModel.verbosity = verb
vbMng(self, "DEL", "Done computing numerator.", 7)
return ps
def setupApprox(self):
"""
Compute Pade' interpolant.
SVD-based robust eigenvalue management.
"""
if self.checkComputedApprox():
return
RROMPyAssert(self._mode, message = "Cannot setup approximant.")
vbMng(self, "INIT", "Setting up {}.". format(self.name()), 5)
self.computeScaleFactor()
self.computeSnapshots()
if self.trainedModel is None:
self.trainedModel = tModel()
self.trainedModel.verbosity = self.verbosity
self.trainedModel.timestamp = self.timestamp
data = TrainedModelPivotedData(self.trainedModel.name(), self.mu0,
self.samplingEngine.samplesCoalesced,
self.HFEngine.rescalingExp,
self.directionPivot,
self.scaleFactor)
data.musPivot = copy(self.musPivot)
data.musMarginal = copy(self.musMarginal)
self.trainedModel.data = data
else:
self.trainedModel = self.trainedModel
self.trainedModel.data.projMat = copy(self.samplingEngine.samples)
self.trainedModel.data.marginalInterp = self._setupMarginalInterp()
if self.N > 0:
Qs = self._setupDenominator()
else:
Q = PI()
Q.npar = self.musPivot.shape[1]
Q.coeffs = np.ones(tuple([1] * Q.npar), dtype = np.complex)
Q.polybasis = self.polybasisPivot
Qs = [Q for _ in range(len(self.musMarginal))]
self.trainedModel.data.Qs = Qs
self.trainedModel.data.Ps = self._setupNumerator()
self.trainedModel.data.approxParameters = copy(self.approxParameters)
vbMng(self, "DEL", "Done setting up approximant.", 5)
def _computeInterpolantInverseBlocks(self) -> List[Np2D]:
"""
Compute inverse factors for minimal interpolant target functional.
"""
RROMPyAssert(self._mode, message = "Cannot solve eigenvalue problem.")
self._setupPivotInterpolationIndices()
while self.E >= 0:
eWidth = (hashD([self.E + 1] + [0] * (self.musPivot.shape[1] - 1))
- hashD([self.E] + [0] * (self.musPivot.shape[1] - 1)))
TE, _, argIdxs = hpvP(self._musPUniqueCN, self.E,
self.polybasisPivot, self._derPIdxs,
self._reorderP,
scl = np.power(self.scaleFactorPivot, -1.))
fitOut = customPInv(TE[:, argIdxs], rcond = self.interpRcondPivot,
full = True)
vbMng(self, "MAIN",
("Fitting {} samples with degree {} through {}... "
"Conditioning of pseudoinverse system: {:.4e}.").format(
TE.shape[0], self.E,
polyfitname(self.polybasisPivot),
fitOut[1][1][0] / fitOut[1][1][-1]),
5)
if fitOut[1][0] == len(argIdxs):
self._fitinvP = fitOut[0][- eWidth : , :]
break
RROMPyWarning("Polyfit is poorly conditioned. Reducing E by 1.")
self.E -= 1
if self.E < 0:
raise RROMPyException(("Instability in computation of "
"denominator. Aborting."))
TN, _, argIdxs = hpvP(self._musPUniqueCN, self.N, self.polybasisPivot,
self._derPIdxs, self._reorderP,
scl = np.power(self.scaleFactorPivot, -1.))
TN = TN[:, argIdxs]
invD = [None] * (eWidth)
for k in range(eWidth):
pseudoInv = np.diag(self._fitinvP[k, :])
idxGlob = 0
for j, derIdxs in enumerate(self._derPIdxs):
nder = len(derIdxs)
idxGlob += nder
if nder > 1:
idxLoc = np.arange(len(self.musPivot))[
(self._reorderP >= idxGlob - nder)
* (self._reorderP < idxGlob)]
invLoc = self._fitinvP[k, idxLoc]
pseudoInv[np.ix_(idxLoc, idxLoc)] = 0.
for diffj, diffjIdx in enumerate(derIdxs):
for derQ, derQIdx in enumerate(derIdxs):
derUIdx = [x - y for (x, y) in
zip(diffjIdx, derQIdx)]
if all([x >= 0 for x in derUIdx]):
derU = hashD(derUIdx)
pseudoInv[idxLoc[derU], idxLoc[derQ]] = (
invLoc[diffj])
invD[k] = pseudoInv.dot(TN)
return invD
def centerNormalize(self, mu : paramList = [],
mu0 : paramVal = None) -> paramList:
"""
Compute normalized parameter to be plugged into approximant.
Args:
mu: Parameter(s) 1.
mu0: Parameter(s) 2. If None, set to self.mu0.
Returns:
Normalized parameter.
"""
return self.trainedModel.centerNormalize(mu, mu0)
def centerNormalizePivot(self, mu : paramList = [],
mu0 : paramVal = None) -> paramList:
"""
Compute normalized parameter to be plugged into approximant.
Args:
mu: Parameter(s) 1.
mu0: Parameter(s) 2. If None, set to self.mu0.
Returns:
Normalized parameter.
"""
return self.trainedModel.centerNormalizePivot(mu, mu0)
def getResidues(self) -> Np1D:
"""
Obtain approximant residues.
Returns:
Matrix with residues as columns.
"""
return self.trainedModel.getResidues()
diff --git a/rrompy/reduction_methods/trained_model/trained_model_pivoted_pade.py b/rrompy/reduction_methods/trained_model/trained_model_pivoted_pade.py
index 63adfa0..5da4e6a 100644
--- a/rrompy/reduction_methods/trained_model/trained_model_pivoted_pade.py
+++ b/rrompy/reduction_methods/trained_model/trained_model_pivoted_pade.py
@@ -1,238 +1,238 @@
# 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 .trained_model_pade import TrainedModelPade
from rrompy.utilities.base.types import (Np1D, List, ListAny, paramVal,
paramList, sampList)
from rrompy.utilities.base import verbosityManager as vbMng, freepar as fp
from rrompy.utilities.poly_fitting.polynomial import polyroots
from rrompy.utilities.exception_manager import RROMPyException
from rrompy.parameter import checkParameter, checkParameterList
from rrompy.sampling import sampleList, emptySampleList
__all__ = ['TrainedModelPivotedPade']
class TrainedModelPivotedPade(TrainedModelPade):
"""
ROM approximant evaluation for pivoted Pade' approximant.
Attributes:
Data: dictionary with all that can be pickled.
"""
def centerNormalizePivot(self, mu : paramList = [],
mu0 : paramVal = None) -> paramList:
"""
Compute normalized parameter to be plugged into approximant.
Args:
mu: Parameter(s) 1.
mu0: Parameter(s) 2. If None, set to self.data.mu0.
Returns:
Normalized parameter.
"""
mu = checkParameterList(mu, self.data.nparPivot)[0]
if mu0 is None: mu0 = self.data.mu0(self.data.directionPivot)
rad = ((mu ** self.data.rescalingExpPivot
- mu0 ** self.data.rescalingExpPivot)
/ self.data.scaleFactorPivot)
return rad
def centerNormalizeMarginal(self, mu : paramList = [],
mu0 : paramVal = None) -> paramList:
"""
Compute normalized parameter to be plugged into approximant.
Args:
mu: Parameter(s) 1.
mu0: Parameter(s) 2. If None, set to self.data.mu0.
Returns:
Normalized parameter.
"""
mu = checkParameterList(mu, self.data.nparMarginal)[0]
if mu0 is None: mu0 = self.data.mu0(self.data.directionMarginal)
rad = ((mu ** self.data.rescalingExpMarginal
- mu0 ** self.data.rescalingExpMarginal)
/ self.data.scaleFactorMarginal)
return rad
def interpolateMarginal(self, mu : paramList = [],
samples : ListAny = [], der : List[int] = None,
scl : Np1D = None) -> sampList:
"""
Evaluate marginal interpolator at arbitrary marginal parameter.
Args:
mu: Target parameter.
samples: Objects to interpolate.
der(optional): Derivatives to take before evaluation.
"""
mu = checkParameterList(mu, self.data.nparMarginal)[0]
sList = isinstance(samples[0], sampleList)
sEff = [None] * len(samples)
for j in range(len(samples)):
if sList:
- sEff[j] = samples[j].data.flatten()
+ sEff[j] = samples[j].data
else:
- sEff[j] = samples[j].flatten()
+ sEff[j] = samples[j]
vbMng(self, "INIT", "Interpolating marginal at mu = {}.".format(mu),
35)
muC = self.centerNormalizeMarginal(mu)
p = emptySampleList()
p.reset((len(sEff[0]), len(muC)))
p.data[:] = 0.
for mIj, spj in zip(self.data.marginalInterp, sEff):
- p = p + np.outer(spj, mIj(muC, der, scl))
+ p = p + spj * mIj(muC, der, scl)
vbMng(self, "DEL", "Done interpolating marginal.", 35)
if not sList:
p = p.data.flatten()
return p
def getPValsPivot(self, mu : paramList = [], der : List[int] = None,
scl : Np1D = None) -> sampList:
"""
Evaluate Pade' numerators at arbitrary pivot parameter.
Args:
mu: Target parameter.
der(optional): Derivatives to take before evaluation.
"""
mu = checkParameterList(mu, self.data.nparPivot)[0]
vbMng(self, "INIT", "Evaluating numerator at mu = {}.".format(mu), 17)
muC = self.centerNormalizePivot(mu)
pP = [sampleList(P(muC, der, scl)) for P in self.data.Ps]
vbMng(self, "DEL", "Done evaluating numerator.", 17)
return pP
def getPVal(self, mu : paramList = [], der : List[int] = None,
scl : Np1D = None) -> sampList:
"""
Evaluate Pade' numerator at arbitrary parameter.
Args:
mu: Target parameter.
der(optional): Derivatives to take before evaluation.
"""
mu = checkParameterList(mu, self.data.npar)[0]
muP = checkParameterList(mu.data[:, self.data.directionPivot],
self.data.nparPivot)[0]
muM = checkParameterList(mu.data[:, self.data.directionMarginal],
self.data.nparMarginal)[0]
if der is None:
derP, derM = None, None
else:
derP = [der[x] for x in self.data.directionPivot]
derM = [der[x] for x in self.data.directionMarginal]
if scl is None:
sclP, sclM = None, None
else:
sclP = [scl[x] for x in self.data.directionPivot]
sclM = [scl[x] for x in self.data.directionMarginal]
pP = self.getPValsPivot(muP, derP, sclP)
return self.interpolateMarginal(muM, pP, derM, sclM)
def getQValsPivot(self, mu:Np1D, der : List[int] = None,
scl : Np1D = None) -> Np1D:
"""
Evaluate Pade' denominators at arbitrary pivot parameter.
Args:
mu: Target parameter.
der(optional): Derivatives to take before evaluation.
"""
mu = checkParameterList(mu, self.data.nparPivot)[0]
vbMng(self, "INIT", "Evaluating denominator at mu = {}.".format(mu),
17)
muC = self.centerNormalizePivot(mu)
- qP = [Q(muC, der, scl) for Q in self.data.Qs]
+ qP = [Q(muC, der, scl).reshape(1, -1) for Q in self.data.Qs]
vbMng(self, "DEL", "Done evaluating denominator.", 17)
return qP
def getQVal(self, mu:Np1D, der : List[int] = None,
scl : Np1D = None) -> Np1D:
"""
Evaluate Pade' denominator at arbitrary parameter.
Args:
mu: Target parameter.
der(optional): Derivatives to take before evaluation.
"""
mu = checkParameterList(mu, self.data.npar)[0]
muP = checkParameterList(mu.data[:, self.data.directionPivot],
self.data.nparPivot)[0]
muM = checkParameterList(mu.data[:, self.data.directionMarginal],
self.data.nparMarginal)[0]
if der is None:
derP, derM = None, None
else:
derP = [der[x] for x in self.data.directionPivot]
derM = [der[x] for x in self.data.directionMarginal]
if scl is None:
sclP, sclM = None, None
else:
sclP = [scl[x] for x in self.data.directionPivot]
sclM = [scl[x] for x in self.data.directionMarginal]
qP = self.getQValsPivot(muP, derP, sclP)
return self.interpolateMarginal(muM, qP, derM, sclM)
def getPoles(self, *args, **kwargs) -> Np1D:
"""
Obtain approximant poles.
Returns:
Numpy complex vector of poles.
"""
if len(args) + len(kwargs) > 1:
raise RROMPyException(("Wrong number of parameters passed. "
"Only 1 available."))
elif len(args) + len(kwargs) == 1:
if len(args) == 1:
mVals = args[0]
else:
mVals = kwargs["marginalVals"]
if not hasattr(mVals, "__len__"): mVals = [mVals]
mVals = list(mVals)
else:
mVals = [fp]
mValsP = [mVals[x] for x in self.data.directionPivot]
mValsM = [mVals[x] for x in self.data.directionMarginal]
try:
rDim = mValsP.index(fp)
if rDim < len(mValsP) - 1 and fp in mValsP[rDim + 1 :]:
raise
except:
raise RROMPyException(("Exactly 1 'freepar' entry in "
"marginalVals must be provided."))
if fp in mValsM:
raise RROMPyException(("'freepar' entry in marginalVals must be "
"pivot dimension."))
rDimB = self.data.directionPivot[rDim]
- Qeff = self.interpolateMarginal(mValsM,
- [Q.coeffs for Q in self.data.Qs])
+ Qeff = self.interpolateMarginal(mValsM, [Q.coeffs[:, np.newaxis] \
+ for Q in self.data.Qs])
Qeff = Qeff.reshape(self.data.Qs[0].coeffs.shape)
mValsP[rDim] = self.data.mu0(rDimB)
mValsP = self.centerNormalizePivot(checkParameter(mValsP))
mValsP = list(mValsP.data.flatten())
mValsP[rDim] = fp
return np.power(self.data.mu0(rDimB) ** self.data.rescalingExp[rDimB]
+ self.data.scaleFactor[rDimB]
* polyroots(Qeff, self.data.Qs[0].polybasis, mValsP),
1. / self.data.rescalingExp[rDimB])
diff --git a/rrompy/sampling/base/sampling_engine_base_pivoted.py b/rrompy/sampling/base/sampling_engine_base_pivoted.py
index 418b54b..1d6c889 100644
--- a/rrompy/sampling/base/sampling_engine_base_pivoted.py
+++ b/rrompy/sampling/base/sampling_engine_base_pivoted.py
@@ -1,202 +1,208 @@
# Copyright (C) 2018 by the RROMPy authors
#
# This file is part of RROMPy.
#
# RROMPy is free software: you can redistribute it and/or modify
# it under the terms of the GNU Lesser General Public License as published by
# the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
#
# RROMPy is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU Lesser General Public License for more details.
#
# You should have received a copy of the GNU Lesser General Public License
# along with RROMPy. If not, see .
#
import numpy as np
from rrompy.utilities.base.types import (Np1D, HFEng, List, ListAny, strLst,
paramVal, paramList, sampList)
from rrompy.utilities.base import verbosityManager as vbMng
from rrompy.utilities.exception_manager import RROMPyWarning
from rrompy.parameter import (emptyParameterList, checkParameter,
checkParameterList)
from rrompy.sampling import emptySampleList
from .sampling_engine_base import SamplingEngineBase
__all__ = ['SamplingEngineBasePivoted']
class SamplingEngineBasePivoted(SamplingEngineBase):
"""HERE"""
def __init__(self, HFEngine:HFEng, directionPivot:ListAny,
verbosity : int = 10, timestamp : bool = True,
allowRepeatedSamples : bool = True):
super().__init__(HFEngine, verbosity, timestamp)
self.directionPivot = directionPivot
self.HFEngineMarginalized = None
self.resetHistory()
@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):
return np.sum(self.nsamples)
def resetHistory(self, j : int = 1):
self.samples = [emptySampleList() for _ in range(j)]
self.nsamples = [0] * j
self.mus = [emptyParameterList() for _ in range(j)]
self._derIdxs = [[] for _ in range(j)]
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()
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):
self.samplesCoalesced = emptySampleList()
self.samplesCoalesced.reset((self.samples[0].shape[0],
np.sum([samp.shape[1] \
for samp in self.samples])),
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
def solveLS(self, mu : paramList = [], RHS : sampList = None,
homogeneized : bool = False) -> 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 mu = {}.".format(mu), 15)
+ vbMng(self, "INIT",
+ ("Solving HF model for mu = {} for marginal "
+ "{}.").format(mu, self.HFEngineMarginalized.muFixed), 15)
u = self.HFEngineMarginalized.solve(mu, RHS, homogeneized)
vbMng(self, "DEL", "Done solving HF model.", 15)
return u
def plotSamples(self, warping : List[callable] = None, name : str = "u",
save : str = None, what : strLst = 'all',
saveFormat : str = "eps", saveDPI : int = 100,
show : bool = True, plotArgs : dict = {}, **figspecs):
"""
Do some nice plots of the samples.
Args:
warping(optional): Domain warping functions.
name(optional): Name to be shown as title of the plots. Defaults to
'u'.
save(optional): Where to save plot(s). Defaults to None, i.e. no
saving.
what(optional): Which plots to do. If list, can contain 'ABS',
'PHASE', 'REAL', 'IMAG'. If str, same plus wildcard 'ALL'.
Defaults to 'ALL'.
saveFormat(optional): Format for saved plot(s). Defaults to "eps".
saveDPI(optional): DPI for saved plot(s). Defaults to 100.
show(optional): Whether to show figure. Defaults to True.
plotArgs(optional): Optional arguments for fen/pyplot.
figspecs(optional key args): Optional arguments for matplotlib
figure creation.
"""
for i in range(len(self.nsamples)):
for j in range(self.nsamples[i]):
self.HFEngine.plot(self.samples[i][j], warping,
"{}_{}_{}".format(name, i, j), save, what,
saveFormat, saveDPI, show, plotArgs,
**figspecs)
def outParaviewSamples(self, name : str = "u", folders : bool = True,
filename : str = "out", times : Np1D = None,
what : strLst = 'all', forceNewFile : bool = True,
filePW = None):
"""
Output samples to ParaView file.
Args:
name(optional): Base name to be used for data output.
folders(optional): Whether to split output in folders.
filename(optional): Name of output file.
times(optional): Timestamps.
what(optional): Which plots to do. If list, can contain 'MESH',
'ABS', 'PHASE', 'REAL', 'IMAG'. If str, same plus wildcard
'ALL'. Defaults to 'ALL'.
forceNewFile(optional): Whether to create new output file.
filePW(optional): Fenics File entity (for time series).
"""
if times is None: times = [[0.] * self.nsamples[i] \
for i in range(len(self.nsamples))]
for i in range(len(self.nsamples)):
for j in range(self.nsamples[i]):
self.HFEngine.outParaview(self.samples[i][j],
name = "{}_{}_{}".format(name, i, j),
filename = "{}_{}_{}".format(filename,
i, j),
time = times[i][j], what = what,
forceNewFile = forceNewFile,
folder = folders, filePW = filePW)
def outParaviewTimeDomainSamples(self, omegas : Np1D = None,
timeFinal : Np1D = None,
periodResolution : int = 20,
name : str = "u", folders : bool = True,
filename : str = "out",
forceNewFile : bool = True):
"""
Output samples to ParaView file, converted to time domain.
Args:
omegas(optional): frequencies.
timeFinal(optional): final time of simulation.
periodResolution(optional): number of time steps per period.
name(optional): Base name to be used for data output.
folders(optional): Whether to split output in folders.
filename(optional): Name of output file.
forceNewFile(optional): Whether to create new output file.
"""
if omegas is None: omegas = [[np.real(self.mus[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))]
for i in range(len(self.nsamples)):
for j in range(self.nsamples[i]):
self.HFEngine.outParaviewTimeDomain(self.samples[i][j],
omega = omegas[i][j],
timeFinal = timeFinal[i][j],
periodResolution = periodResolution,
name = "{}_{}_{}".format(name, i, j),
filename = "{}_{}_{}".format(filename,
i, j),
forceNewFile = forceNewFile,
folder = folders)
diff --git a/rrompy/sampling/pivoted/sampling_engine_pivoted.py b/rrompy/sampling/pivoted/sampling_engine_pivoted.py
index 34275b0..d164869 100644
--- a/rrompy/sampling/pivoted/sampling_engine_pivoted.py
+++ b/rrompy/sampling/pivoted/sampling_engine_pivoted.py
@@ -1,119 +1,128 @@
# Copyright (C) 2018 by the RROMPy authors
#
# This file is part of RROMPy.
#
# RROMPy is free software: you can redistribute it and/or modify
# it under the terms of the GNU Lesser General Public License as published by
# the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
#
# RROMPy is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU Lesser General Public License for more details.
#
# You should have received a copy of the GNU Lesser General Public License
# along with RROMPy. If not, see .
#
from copy import deepcopy as copy
import numpy as np
from rrompy.sampling.base.sampling_engine_base_pivoted import (
SamplingEngineBasePivoted)
+from rrompy.hfengines.base import MarginalProxyEngine
from rrompy.utilities.base.types import Np1D, paramVal, paramList, sampList
-from rrompy.utilities.base import verbosityManager as vbMng
+from rrompy.utilities.base import verbosityManager as vbMng, freepar as fp
from rrompy.utilities.exception_manager import RROMPyException
from rrompy.utilities.poly_fitting.polynomial import nextDerivativeIndices
from rrompy.parameter import checkParameter, checkParameterList
from rrompy.sampling import sampleList
__all__ = ['SamplingEnginePivoted']
class SamplingEnginePivoted(SamplingEngineBasePivoted):
"""HERE"""
def preprocesssamples(self, idxs:Np1D, j:int) -> sampList:
if self.samples[j] is None or len(self.samples[j]) == 0: return
return self.samples[j](idxs)
def postprocessu(self, u:sampList, j:int,
overwrite : bool = False) -> Np1D:
return copy(u)
def postprocessuBulk(self, u:sampList, j:int) -> sampList:
return copy(u)
def lastSampleManagement(self, j:int):
pass
def _getSampleConcurrence(self, mu:paramVal, j:int, previous:Np1D,
homogeneized : bool = False) -> sampList:
if len(previous) >= len(self._derIdxs[j]):
self._derIdxs[j] += nextDerivativeIndices(
self._derIdxs[j], self.nPivot,
len(previous) + 1 - len(self._derIdxs[j]))
derIdx = self._derIdxs[j][len(previous)]
mu = checkParameter(mu, self.nPivot)
samplesOld = self.preprocesssamples(previous, j)
RHS = self.HFEngineMarginalized.b(mu, derIdx,
homogeneized = homogeneized)
for j, derP in enumerate(self._derIdxs[j][: len(previous)]):
diffP = [x - y for (x, y) in zip(derIdx, derP)]
if np.all([x >= 0 for x in diffP]):
RHS -= self.HFEngineMarginalized.A(mu, diffP).dot(
samplesOld[j])
return self.solveLS(mu, RHS = RHS, homogeneized = homogeneized)
def nextSample(self, mu:paramVal, j:int, overwrite : bool = False,
homogeneized : bool = False,
lastSample : bool = True) -> Np1D:
mu = checkParameter(mu, self.nPivot)
ns = self.nsamples[j]
muidxs = self.mus[j].findall(mu[0])
if len(muidxs) > 0:
u = self._getSampleConcurrence(mu, j, np.sort(muidxs),
homogeneized)
else:
u = self.solveLS(mu, homogeneized = homogeneized)
u = self.postprocessu(u, j, overwrite = overwrite)
if overwrite:
self.samples[j][ns] = u
self.mus[j][ns] = mu[0]
else:
if ns == 0:
self.samples[j] = sampleList(u)
else:
self.samples[j].append(u)
self.mus[j].append(mu)
self.nsamples[j] += 1
if lastSample: self.lastSampleManagement(j)
return u
- def iterSample(self, mus:paramList, j:int,
+ def iterSample(self, mus:paramList, musM:paramList,
homogeneized : bool = False) -> sampList:
- if self.HFEngineMarginalized is None:
- raise RROMPyException(("Marginalized HFEngine must be assigned "
- "before sample computation."))
mus = checkParameterList(mus, self.nPivot)[0]
+ musM = checkParameterList(musM, self.nMarginal)[0]
vbMng(self, "INIT", "Starting sampling iterations.", 5)
n = len(mus)
+ m = len(musM)
if n <= 0:
- raise RROMPyException(("Number of samples must be positive."))
-
- if self.allowRepeatedSamples:
- for k in range(n):
- vbMng(self, "MAIN",
- "Computing sample {} / {}.".format(k + 1, n), 7)
- self.nextSample(mus[k], j, overwrite = (k > 0),
- homogeneized = homogeneized,
- lastSample = (n == k + 1))
- if k == 0:
- self.preallocateSamples(self.samples[j][0], mus[0], n, j)
- else:
- self.samples[j] = self.postprocessuBulk(self.solveLS(mus,
+ raise RROMPyException("Number of samples must be positive.")
+ if m <= 0:
+ raise RROMPyException(("Number of marginal samples must be "
+ "positive."))
+ for j in range(m):
+ muMEff = [fp] * self.HFEngine.npar
+ for k, x in enumerate(self.directionMarginal):
+ muMEff[x] = musM(j, k)
+ self.HFEngineMarginalized = MarginalProxyEngine(self.HFEngine,
+ list(muMEff))
+ if self.allowRepeatedSamples:
+ for k in range(n):
+ vbMng(self, "MAIN",
+ "Computing sample {} / {} for marginal {} / {}."\
+ .format(k + 1, n, j, m), 10)
+ self.nextSample(mus[k], j, overwrite = (k > 0),
+ homogeneized = homogeneized,
+ lastSample = (n == k + 1))
+ if k == 0:
+ self.preallocateSamples(self.samples[j][0], mus[0],
+ n, j)
+ else:
+ self.samples[j] = self.postprocessuBulk(self.solveLS(mus,
homogeneized = homogeneized), j)
- self.mus[j] = copy(mus)
- self.nsamples[j] = n
+ self.mus[j] = copy(mus)
+ self.nsamples[j] = n
vbMng(self, "DEL", "Finished sampling iterations.", 5)
return self.samples[j]
-
diff --git a/rrompy/sampling/pivoted/sampling_engine_pivoted_pod.py b/rrompy/sampling/pivoted/sampling_engine_pivoted_pod.py
index e738e9d..af32fb3 100644
--- a/rrompy/sampling/pivoted/sampling_engine_pivoted_pod.py
+++ b/rrompy/sampling/pivoted/sampling_engine_pivoted_pod.py
@@ -1,113 +1,114 @@
# 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 copy import deepcopy as copy
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):
"""HERE"""
def resetHistory(self, j : int = 1):
super().resetHistory(j)
self.samples_full = [None] * j
self.RPOD = [None] * j
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, tol : float = 1e-12):
super().coalesceSamples()
self.samplesCoalesced, RPODC = (
self.PODEngine.generalizedQR(self.samplesCoalesced))
self.RPODCoalesced = np.zeros((self.samplesCoalesced.shape[1],
self.samplesCoalesced.shape[1]),
dtype = self.RPOD[0].dtype)
self.samples_fullCoalesced = emptySampleList()
self.samples_fullCoalesced.reset((self.samples_full[0].shape[0],
self.samplesCoalesced.shape[1]),
self.samples_full[0].dtype)
ci = 0
for j, (Rloc, samp) in enumerate(zip(self.RPOD, self.samples_full)):
ri = 0
Rheg = Rloc.shape[1]
for k, Rloc2 in enumerate(self.RPOD[: j + 1]):
Rlen = Rloc2.shape[1]
self.RPODCoalesced[ri : ri + Rlen, ci : ci + Rheg] = (
RPODC[ri : ri + Rlen, ci : ci + Rheg].dot(Rloc))
ri += Rlen
self.samples_fullCoalesced.data[:, ci : ci + Rheg] = samp.data
ci += Rheg
RCdiag = np.abs(np.diag(self.RPODCoalesced))
RCdiag /= RCdiag[0]
ntrunc = np.nonzero(RCdiag < tol)[0]
if len(ntrunc) == 0: return
self.samplesCoalesced.data = self.samplesCoalesced.data[:, : ntrunc[0]]
self.RPODCoalesced = self.RPODCoalesced[: ntrunc[0], :]
@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)
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 postprocessu(self, u:sampList, j:int,
overwrite : bool = False) -> Np1D:
ns = self.nsamples[j]
if overwrite:
self.samples_full[j][ns] = copy(u)
else:
if ns == 0:
self.samples_full[j] = sampleList(u)
else:
self.samples_full[j].append(u)
return u
def postprocessuBulk(self, u:sampList, j:int) -> sampList:
self.samples_full[j] = copy(u)
- vbMng(self, "INIT", "Starting orthogonalization.", 10)
+ vbMng(self, "INIT",
+ "Starting orthogonalization for marginal {}.".format(j), 40)
u, self.RPOD[j] = self.PODEngine.generalizedQR(self.samples_full[j])
- vbMng(self, "DEL", "Done orthogonalizing.", 10)
+ vbMng(self, "DEL", "Done orthogonalizing.", 40)
return u
def lastSampleManagement(self, j:int):
self.samples[j] = self.postprocessuBulk(self.samples_full[j], j)
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