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