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