diff --git a/rrompy/reduction_methods/greedy/rational_interpolant_greedy.py b/rrompy/reduction_methods/greedy/rational_interpolant_greedy.py index e8b67ba..813100e 100644 --- a/rrompy/reduction_methods/greedy/rational_interpolant_greedy.py +++ b/rrompy/reduction_methods/greedy/rational_interpolant_greedy.py @@ -1,514 +1,522 @@ # Copyright (C) 2018 by the RROMPy authors # # This file is part of RROMPy. # # RROMPy is free software: you can redistribute it and/or modify # it under the terms of the GNU Lesser General Public License as published by # the Free Software Foundation, either version 3 of the License, or # (at your option) any later version. # # RROMPy is distributed in the hope that it will be useful, # but WITHOUT ANY WARRANTY; without even the implied warranty of # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the # GNU Lesser General Public License for more details. # # You should have received a copy of the GNU Lesser General Public License # along with RROMPy. If not, see . # from copy import deepcopy as copy import numpy as np from rrompy.hfengines.base.linear_affine_engine import checkIfAffine from .generic_greedy_approximant import GenericGreedyApproximant from rrompy.utilities.poly_fitting.polynomial import (polybases, PolynomialInterpolator as PI, polyvanderTotal as pvT) from rrompy.utilities.numerical import totalDegreeN, dot from rrompy.utilities.expression import expressionEvaluator from rrompy.reduction_methods.standard import RationalInterpolant from rrompy.utilities.base.types import (Np1D, Tuple, DictAny, HFEng, paramVal, List) from rrompy.utilities.base import verbosityManager as vbMng from rrompy.utilities.poly_fitting import customFit from rrompy.utilities.exception_manager import (RROMPyWarning, RROMPyException, RROMPyAssert, RROMPy_FRAGILE) from rrompy.parameter import checkParameterList from rrompy.sampling import emptySampleList __all__ = ['RationalInterpolantGreedy'] class RationalInterpolantGreedy(GenericGreedyApproximant, RationalInterpolant): """ ROM greedy rational interpolant computation for parametric problems. Args: HFEngine: HF problem solver. mu0(optional): Default parameter. Defaults to 0. approxParameters(optional): Dictionary containing values for main parameters of approximant. Recognized keys are: - 'POD': whether to compute POD of snapshots; defaults to True; - 'S': number of starting training points; - 'sampler': sample point generator; - 'greedyTol': uniform error tolerance for greedy algorithm; defaults to 1e-2; - 'collinearityTol': collinearity tolerance for greedy algorithm; defaults to 0.; - 'maxIter': maximum number of greedy steps; defaults to 1e2; - 'refinementRatio': ratio of test points to be exhausted before test set refinement; defaults to 0.2; - 'nTestPoints': number of test points; defaults to 5e2; - 'trainSetGenerator': training sample points generator; defaults to sampler; - 'polybasis': type of basis for interpolation; defaults to 'MONOMIAL'; - 'errorEstimatorKind': kind of error estimator; available values include 'AFFINE', 'DISCREPANCY', 'INTERPOLATORY', 'LOOK_AHEAD', and 'NONE'; defaults to 'NONE'; - 'interpRcond': tolerance for interpolation; defaults to None; - 'robustTol': tolerance for robust rational denominator management; defaults to 0. Defaults to empty dict. approx_state(optional): Whether to approximate state. Defaults and must be True. verbosity(optional): Verbosity level. Defaults to 10. Attributes: HFEngine: HF problem solver. mu0: Default parameter. mus: Array of snapshot parameters. approxParameters: Dictionary containing values for main parameters of approximant. Recognized keys are in parameterList. parameterListSoft: Recognized keys of soft approximant parameters: - 'POD': whether to compute POD of snapshots. - 'greedyTol': uniform error tolerance for greedy algorithm; - 'collinearityTol': collinearity tolerance for greedy algorithm; - 'maxIter': maximum number of greedy steps; - 'refinementRatio': ratio of test points to be exhausted before test set refinement; - 'nTestPoints': number of test points; - 'trainSetGenerator': training sample points generator; - 'errorEstimatorKind': kind of error estimator; - 'interpRcond': tolerance for interpolation; - 'robustTol': tolerance for robust rational denominator management. parameterListCritical: Recognized keys of critical approximant parameters: - 'S': total number of samples current approximant relies upon; - 'sampler': sample point generator. approx_state: Whether to approximate state. verbosity: Verbosity level. POD: whether to compute POD of snapshots. S: number of test points. sampler: Sample point generator. greedyTol: uniform error tolerance for greedy algorithm. collinearityTol: Collinearity tolerance for greedy algorithm. maxIter: maximum number of greedy steps. refinementRatio: ratio of training points to be exhausted before training set refinement. nTestPoints: number of starting training points. trainSetGenerator: training sample points generator. robustTol: tolerance for robust rational denominator management. errorEstimatorKind: kind of error estimator. interpRcond: tolerance for interpolation. robustTol: tolerance for robust rational denominator management. muBounds: list of bounds for parameter values. samplingEngine: Sampling engine. estimatorNormEngine: Engine for estimator norm computation. uHF: High fidelity solution(s) with parameter(s) lastSolvedHF as sampleList. lastSolvedHF: Parameter(s) corresponding to last computed high fidelity solution(s) as parameterList. uApproxReduced: Reduced approximate solution(s) with parameter(s) lastSolvedApprox as sampleList. lastSolvedApproxReduced: Parameter(s) corresponding to last computed reduced approximate solution(s) as parameterList. uApprox: Approximate solution(s) with parameter(s) lastSolvedApprox as sampleList. lastSolvedApprox: Parameter(s) corresponding to last computed approximate solution(s) as parameterList. """ _allowedEstimatorKinds = ["AFFINE", "DISCREPANCY", "INTERPOLATORY", "LOOK_AHEAD", "NONE"] def __init__(self, HFEngine:HFEng, mu0 : paramVal = None, approxParameters : DictAny = {}, approx_state : bool = True, verbosity : int = 10, timestamp : bool = True): if not approx_state: RROMPyWarning("Overriding approx_state to True.") self._preInit() self._addParametersToList(["errorEstimatorKind"], ["DISCREPANCY"], toBeExcluded = ["M", "N", "polydegreetype", "radialDirectionalWeights", "nNearestNeighbor"]) super().__init__(HFEngine = HFEngine, mu0 = mu0, approxParameters = approxParameters, approx_state = True, verbosity = verbosity, timestamp = timestamp) self._postInit() @property def E(self): """Value of E.""" self._E = self.sampleBatchIdx - 1 return self._E @E.setter def E(self, E): RROMPyWarning(("E is used just to simplify inheritance, and its value " "cannot be changed from that of sampleBatchIdx - 1.")) @property def polydegreetype(self): """Value of polydegreetype.""" return "TOTAL" @polydegreetype.setter def polydegreetype(self, polydegreetype): RROMPyWarning(("polydegreetype is used just to simplify inheritance, " "and its value cannot be changed from 'TOTAL'.")) @property def polybasis(self): """Value of polybasis.""" return self._polybasis @polybasis.setter def polybasis(self, polybasis): try: polybasis = polybasis.upper().strip().replace(" ","") if polybasis not in polybases: raise RROMPyException("Sample type not recognized.") self._polybasis = polybasis except: RROMPyWarning(("Prescribed polybasis not recognized. Overriding " "to 'MONOMIAL'.")) self._polybasis = "MONOMIAL" self._approxParameters["polybasis"] = self.polybasis @property def errorEstimatorKind(self): """Value of errorEstimatorKind.""" return self._errorEstimatorKind @errorEstimatorKind.setter def errorEstimatorKind(self, errorEstimatorKind): errorEstimatorKind = errorEstimatorKind.upper() if errorEstimatorKind not in self._allowedEstimatorKinds: RROMPyWarning(("Error estimator kind not recognized. Overriding " "to 'NONE'.")) errorEstimatorKind = "NONE" self._errorEstimatorKind = errorEstimatorKind self._approxParameters["errorEstimatorKind"] = self.errorEstimatorKind def _polyvanderAuxiliary(self, mus, deg, *args): return pvT(mus, deg, *args) def getErrorEstimatorDiscrepancy(self, mus:Np1D) -> Np1D: """Discrepancy-based residual estimator.""" mus = checkParameterList(mus, self.npar)[0] muCTest = self.trainedModel.centerNormalize(mus) verb = self.trainedModel.verbosity self.trainedModel.verbosity = 0 QTest = self.trainedModel.getQVal(mus) checkIfAffine(self.HFEngine, "apply discrepancy-based error estimator") self.HFEngine.buildA() self.HFEngine.buildb() nAs, nbs = self.HFEngine.nAs, self.HFEngine.nbs muTrainEff = self.mus ** self.HFEngine.rescalingExp muTestEff = mus ** self.HFEngine.rescalingExp PTrain = self.trainedModel.getPVal(self.mus).data.T QTrain = self.trainedModel.getQVal(self.mus) PTest = self.trainedModel.getPVal(mus).data radiusAbTrain = np.empty((self.S, nAs * self.S + nbs), dtype = np.complex) radiusA = np.empty((self.S, nAs, len(mus)), dtype = np.complex) radiusb = np.empty((nbs, len(mus)), dtype = np.complex) for j, thA in enumerate(self.HFEngine.thAs): idxs = j * self.S + np.arange(self.S) radiusAbTrain[:, idxs] = expressionEvaluator(thA[0], muTrainEff, (self.S, 1)) * PTrain radiusA[:, j] = PTest * expressionEvaluator(thA[0], muTestEff, (len(mus),)) for j, thb in enumerate(self.HFEngine.thbs): idx = nAs * self.S + j radiusAbTrain[:, idx] = QTrain * expressionEvaluator(thb[0], muTrainEff, (self.S,)) radiusb[j] = QTest * expressionEvaluator(thb[0], muTestEff, (len(mus),)) QRHSNorm2 = self._affineResidualMatricesContraction(radiusb) vanTrain = self._polyvanderAuxiliary(self._musUniqueCN, self.E, self.polybasis0, self._derIdxs, self._reorder) interpPQ = customFit(vanTrain, radiusAbTrain, rcond = self.interpRcond) vanTest = self._polyvanderAuxiliary(muCTest, self.E, self.polybasis0) DradiusAb = vanTest.dot(interpPQ) radiusA = (radiusA - DradiusAb[:, : - nbs].reshape(len(mus), -1, self.S).T) radiusb = radiusb - DradiusAb[:, - nbs :].T ff, Lf, LL = self._affineResidualMatricesContraction(radiusb, radiusA) err = np.abs((LL - 2. * np.real(Lf) + ff) / QRHSNorm2) ** .5 self.trainedModel.verbosity = verb return err def getErrorEstimatorInterpolatory(self, mus:Np1D) -> Np1D: """Interpolation-based residual estimator.""" errTest, QTest, idxMaxEst = self._EIMStep(mus) self.initEstimatorNormEngine() mu_muTestSample = mus[idxMaxEst] app_muTestSample = self.getApproxReduced(mu_muTestSample) if self._mode == RROMPy_FRAGILE: if not self.HFEngine.isCEye: raise RROMPyException(("Cannot compute INTERPOLATORY residual " "estimator in fragile mode for " "non-scalar C.")) app_muTestSample = dot(self.trainedModel.data.projMat, app_muTestSample.data) else: app_muTestSample = dot(self.samplingEngine.samples, app_muTestSample) resmu = self.HFEngine.residual(mu_muTestSample, app_muTestSample, post_c = False) RHSmu = self.HFEngine.residual(mu_muTestSample, None, post_c = False) ressamples = (self.estimatorNormEngine.norm(resmu) / self.estimatorNormEngine.norm(RHSmu)) musT = copy(self.mus) musT.append(mu_muTestSample) musT = self.trainedModel.centerNormalize(musT) musC = self.trainedModel.centerNormalize(mus) resT = np.zeros(len(musT), dtype = np.complex) err = np.zeros(len(mus)) for l in range(len(mu_muTestSample)): resT[len(self.mus) + l] = ressamples[l] * QTest[idxMaxEst[l]] p = PI() wellCond, msg = p.setupByInterpolation(musT, resT, self.E + 1, self.polybasis, self.verbosity >= 15, True, {}, {"rcond": self.interpRcond}) err += np.abs(p(musC)) resT[len(self.mus) + l] = 0. err /= QTest vbMng(self, "MAIN", msg, 15) return err def getErrorEstimatorLookAhead(self, mus:Np1D) -> Tuple[Np1D, List[int]]: """Residual estimator based on look-ahead idea.""" errTest, QTest, idxMaxEst = self._EIMStep(mus) self.initEstimatorNormEngine() mu_muTestSample = mus[idxMaxEst] app_muTestSample = self.getApproxReduced(mu_muTestSample) if self._mode == RROMPy_FRAGILE: app_muTestSample = dot(self.trainedModel.data.projMat, app_muTestSample.data) else: app_muTestSample = dot(self.samplingEngine.samples, app_muTestSample) for j, mu in enumerate(mu_muTestSample): self.samplingEngine.nextSample(mu) if hasattr(self.samplingEngine, "samples_full"): uEx = self.samplingEngine.samples_full[-1] else: uEx = self.samplingEngine.samples[-1] if j == 0: solmu = emptySampleList() solmu.reset((len(uEx), len(mu_muTestSample)), dtype = uEx.dtype) solmu[j] = uEx errmu = solmu - app_muTestSample errsamples = (self.estimatorNormEngine.norm(errmu) / self.estimatorNormEngine.norm(solmu)) musT = copy(self.mus) musT.append(mu_muTestSample) musT = self.trainedModel.centerNormalize(musT) musC = self.trainedModel.centerNormalize(mus) errT = np.zeros(len(musT), dtype = np.complex) err = np.zeros(len(mus)) for l in range(len(mu_muTestSample)): errT[len(self.mus) + l] = errsamples[l] * QTest[idxMaxEst[l]] p = PI() wellCond, msg = p.setupByInterpolation(musT, errT, self.E + 1, self.polybasis, self.verbosity >= 15, True, {}, {"rcond": self.interpRcond}) err += np.abs(p(musC)) errT[len(self.mus) + l] = 0. err /= QTest vbMng(self, "MAIN", msg, 15) return err, idxMaxEst def getErrorEstimatorNone(self, mus:Np1D) -> Np1D: """EIM-based residual estimator.""" err = np.max(self._EIMStep(mus, True), axis = 1) err *= self.greedyTol / np.mean(err) return err def _EIMStep(self, mus:Np1D, only_one : bool = False) -> Tuple[Np1D, Np1D, List[int]]: """Residual estimator based on look-ahead idea.""" mus = checkParameterList(mus, self.npar)[0] muCTest = self.trainedModel.centerNormalize(mus) verb = self.trainedModel.verbosity self.trainedModel.verbosity = 0 QTest = self.trainedModel.getQVal(mus) QTest = np.abs(QTest) muCTest = self.trainedModel.centerNormalize(mus) muCTrain = self.trainedModel.centerNormalize(self.mus) vanTest = self._polyvanderAuxiliary(muCTest, self.E, self.polybasis) vanTestNext = self._polyvanderAuxiliary(muCTest, self.E + 1, self.polybasis)[:, vanTest.shape[1] :] idxsTest = np.arange(vanTestNext.shape[1]) basis = np.zeros((len(idxsTest), 0), dtype = float) idxMaxEst = [] while len(idxsTest) > 0: vanTrial = self._polyvanderAuxiliary(muCTrain, self.E, self.polybasis) vanTrialNext = self._polyvanderAuxiliary(muCTrain, self.E + 1, self.polybasis)[:, vanTrial.shape[1] :] vanTrial = np.hstack((vanTrial, vanTrialNext.dot(basis).reshape( len(vanTrialNext), basis.shape[1]))) valuesTrial = vanTrialNext[:, idxsTest] vanTestEff = np.hstack((vanTest, vanTestNext.dot(basis).reshape( len(vanTestNext), basis.shape[1]))) vanTestNextEff = vanTestNext[:, idxsTest] coeffTest = np.linalg.solve(vanTrial, valuesTrial) errTest = (np.abs(vanTestNextEff - vanTestEff.dot(coeffTest)) / np.expand_dims(QTest, 1)) - if only_one: return errTest + if only_one: + self.trainedModel.verbosity = verb + return errTest idxMaxErr = np.unravel_index(np.argmax(errTest), errTest.shape) idxMaxEst += [idxMaxErr[0]] muCTrain.append(muCTest[idxMaxErr[0]]) basis = np.pad(basis, [(0, 0), (0, 1)], "constant") basis[idxsTest[idxMaxErr[1]], -1] = 1. idxsTest = np.delete(idxsTest, idxMaxErr[1]) self.trainedModel.verbosity = verb return errTest, QTest, idxMaxEst def errorEstimator(self, mus:Np1D, return_max : bool = False) -> Np1D: """Standard residual-based error estimator.""" setupOK = self.setupApproxLocal() if not setupOK: err = np.empty(len(mus)) err[:] = np.nan if not return_max: return err return err, 0, np.nan mus = checkParameterList(mus, self.npar)[0] vbMng(self.trainedModel, "INIT", "Evaluating error estimator at mu = {}.".format(mus), 10) if self.errorEstimatorKind == "AFFINE": err = self.getErrorEstimatorAffine(mus) else: self._setupInterpolationIndices() if self.errorEstimatorKind == "DISCREPANCY": err = self.getErrorEstimatorDiscrepancy(mus) elif self.errorEstimatorKind == "INTERPOLATORY": err = self.getErrorEstimatorInterpolatory(mus) elif self.errorEstimatorKind == "LOOK_AHEAD": err, idxMaxEst = self.getErrorEstimatorLookAhead(mus) else: #if self.errorEstimatorKind == "NONE": err = self.getErrorEstimatorNone(mus) vbMng(self.trainedModel, "DEL", "Done evaluating error estimator", 10) if not return_max: return err if self.errorEstimatorKind != "LOOK_AHEAD": idxMaxEst = np.empty(self.sampleBatchSize, dtype = int) errCP = copy(err) for j in range(self.sampleBatchSize): k = np.argmax(errCP) idxMaxEst[j] = k if j + 1 < self.sampleBatchSize: musZero = self.trainedModel.centerNormalize(mus, mus[k]) errCP *= np.linalg.norm(musZero.data, axis = 1) return err, idxMaxEst, err[idxMaxEst] def greedyNextSample(self, muidx:int, plotEst : bool = False)\ -> Tuple[Np1D, int, float, paramVal]: """Compute next greedy snapshot of solution map.""" RROMPyAssert(self._mode, message = "Cannot add greedy sample.") self.sampleBatchIdx += 1 self.sampleBatchSize = totalDegreeN(self.npar - 1, self.sampleBatchIdx) err, muidx, maxErr, muNext = super().greedyNextSample(muidx, plotEst) if maxErr is not None and (np.any(np.isnan(maxErr)) or np.any(np.isinf(maxErr))): self.sampleBatchIdx -= 1 self.sampleBatchSize = totalDegreeN(self.npar - 1, self.sampleBatchIdx) if (self.errorEstimatorKind == "NONE" and not np.isnan(maxErr) and not np.isinf(maxErr)): maxErr = None return err, muidx, maxErr, muNext def _preliminaryTraining(self): """Initialize starting snapshots of solution map.""" RROMPyAssert(self._mode, message = "Cannot start greedy algorithm.") if self.samplingEngine.nsamples > 0: return S = self.S self.sampleBatchIdx, self.sampleBatchSize, self._S = -1, 0, 0 nextBatchSize = 1 while self._S + nextBatchSize <= S: self.sampleBatchIdx += 1 self.sampleBatchSize = nextBatchSize self._S += self.sampleBatchSize nextBatchSize = totalDegreeN(self.npar - 1, self.sampleBatchIdx + 1) super()._preliminaryTraining() def setupApproxLocal(self): """ Compute rational interpolant. SVD-based robust eigenvalue management. """ if self.checkComputedApprox(): return True RROMPyAssert(self._mode, message = "Cannot setup approximant.") self.verbosity -= 10 vbMng(self, "INIT", "Setting up local approximant.", 5) self._N, self._M = [self.E] * 2 pMat = self.samplingEngine.samples.data pMatEff = dot(self.HFEngine.C, pMat) if self.trainedModel is None: self.trainedModel = self.tModelType() self.trainedModel.verbosity = self.verbosity self.trainedModel.timestamp = self.timestamp datadict = {"mu0": self.mu0, "projMat": pMatEff, "scaleFactor": self.scaleFactor, "rescalingExp": self.HFEngine.rescalingExp} 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) self.trainedModel.data.mus = copy(self.mus) - self.catchInstability = True + self.catchInstability = 2 + stable = True if self.N > 0: try: Q = self._setupDenominator()[0] except RROMPyException as RE: RROMPyWarning(RE._msg) - vbMng(self, "DEL", "Done setting up local approximant.", 5) - return False + vbMng(self, "DEL", "", 7) + stable = False else: Q = PI() Q.coeffs = np.ones((1,) * self.npar, dtype = np.complex) Q.npar = self.npar Q.polybasis = self.polybasis - self.trainedModel.data.Q = copy(Q) - try: - self.trainedModel.data.P = copy(self._setupNumerator()) - except RROMPyException as RE: - RROMPyWarning(RE._msg) - vbMng(self, "DEL", "Done setting up local approximant.", 5) - return False - self.trainedModel.data.approxParameters = copy(self.approxParameters) + if stable: + self.trainedModel.data.Q = copy(Q) + try: + P = copy(self._setupNumerator()) + except RROMPyException as RE: + RROMPyWarning(RE._msg) + vbMng(self, "DEL", "", 7) + stable = False + if stable: + self.trainedModel.data.P = copy(P) + self.trainedModel.data.approxParameters = copy( + self.approxParameters) vbMng(self, "DEL", "Done setting up local approximant.", 5) + self.catchInstability = 0 self.verbosity += 10 - return True + return stable def loadTrainedModel(self, filename:str): """Load trained reduced model from file.""" super().loadTrainedModel(filename) self.sampleBatchIdx, self.sampleBatchSize, _S = -1, 0, 0 nextBatchSize = 1 while _S + nextBatchSize <= self.S + 1: self.sampleBatchIdx += 1 self.sampleBatchSize = nextBatchSize _S += self.sampleBatchSize nextBatchSize = totalDegreeN(self.npar - 1, self.sampleBatchIdx + 1) diff --git a/rrompy/reduction_methods/standard/rational_interpolant.py b/rrompy/reduction_methods/standard/rational_interpolant.py index 36f8a84..23cd860 100644 --- a/rrompy/reduction_methods/standard/rational_interpolant.py +++ b/rrompy/reduction_methods/standard/rational_interpolant.py @@ -1,631 +1,634 @@ # Copyright (C) 2018 by the RROMPy authors # # This file is part of RROMPy. # # RROMPy is free software: you can redistribute it and/or modify # it under the terms of the GNU Lesser General Public License as published by # the Free Software Foundation, either version 3 of the License, or # (at your option) any later version. # # RROMPy is distributed in the hope that it will be useful, # but WITHOUT ANY WARRANTY; without even the implied warranty of # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the # GNU Lesser General Public License for more details. # # You should have received a copy of the GNU Lesser General Public License # along with RROMPy. If not, see . # from copy import deepcopy as copy import numpy as np from rrompy.reduction_methods.base import checkRobustTolerance from .generic_standard_approximant import GenericStandardApproximant from rrompy.utilities.poly_fitting.polynomial import ( polybases as ppb, polyfitname, polyvander as pvP, polyvanderTotal as pvTP, polyTimes, polyTimesTable, vanderInvTable, PolynomialInterpolator as PI) from rrompy.utilities.poly_fitting.heaviside import rational2heaviside from rrompy.utilities.poly_fitting.radial_basis import (polybases as rbpb, RadialBasisInterpolator as RBI) from rrompy.utilities.poly_fitting.moving_least_squares import ( polybases as mlspb, MovingLeastSquaresInterpolator as MLSI) from rrompy.utilities.base.types import (Np1D, Np2D, HFEng, DictAny, Tuple, List, paramVal, sampList) from rrompy.utilities.base import verbosityManager as vbMng from rrompy.utilities.numerical import (customPInv, dot, fullDegreeN, totalDegreeN, degreeTotalToFull, fullDegreeMaxMask, totalDegreeMaxMask, nextDerivativeIndices) from rrompy.utilities.exception_manager import (RROMPyException, RROMPyAssert, RROMPyWarning) __all__ = ['RationalInterpolant'] class RationalInterpolant(GenericStandardApproximant): """ ROM rational interpolant computation for parametric problems. Args: HFEngine: HF problem solver. mu0(optional): Default parameter. Defaults to 0. approxParameters(optional): Dictionary containing values for main parameters of approximant. Recognized keys are: - 'POD': whether to compute POD of snapshots; defaults to True; - 'S': total number of samples current approximant relies upon; - 'sampler': sample point generator; - 'polybasis': type of polynomial basis for interpolation; defaults to 'MONOMIAL'; - 'M': degree of rational interpolant numerator; defaults to 0; - 'N': degree of rational interpolant denominator; defaults to 0; - 'polydegreetype': type of polynomial degree; defaults to 'TOTAL'; - 'radialDirectionalWeights': radial basis weights for interpolant numerator; defaults to 0, i.e. identity; - 'nNearestNeighbor': mumber of nearest neighbors considered in numerator if polybasis allows; defaults to -1; - 'interpRcond': tolerance for 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. mus: Array of snapshot parameters. approxParameters: Dictionary containing values for main parameters of approximant. Recognized keys are in parameterList. parameterListSoft: Recognized keys of soft approximant parameters: - 'POD': whether to compute POD of snapshots; - 'polybasis': type of polynomial basis for interpolation; - 'M': degree of rational interpolant numerator; - 'N': degree of rational interpolant denominator; - 'polydegreetype': type of polynomial degree; - 'radialDirectionalWeights': radial basis weights for interpolant numerator; - 'nNearestNeighbor': mumber of nearest neighbors considered in numerator if polybasis allows; - 'interpRcond': tolerance for interpolation via numpy.polyfit; - '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 samples current approximant relies upon; - 'sampler': sample point generator. approx_state: Whether to approximate state. verbosity: Verbosity level. POD: Whether to compute POD of snapshots. S: Number of solution snapshots over which current approximant is based upon. sampler: Sample point generator. polybasis: type of polynomial basis for interpolation. M: Numerator degree of approximant. N: Denominator degree of approximant. polydegreetype: Type of polynomial degree. radialDirectionalWeights: Radial basis weights for interpolant numerator. nNearestNeighbor: Number of nearest neighbors considered in numerator if polybasis allows. interpRcond: Tolerance for interpolation via numpy.polyfit. 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 parameter values. samplingEngine: Sampling engine. uHF: High fidelity solution(s) with parameter(s) lastSolvedHF as sampleList. lastSolvedHF: Parameter(s) corresponding to last computed high fidelity solution(s) as parameterList. uApproxReduced: Reduced approximate solution(s) with parameter(s) lastSolvedApprox as sampleList. lastSolvedApproxReduced: Parameter(s) corresponding to last computed reduced approximate solution(s) as parameterList. uApprox: Approximate solution(s) with parameter(s) lastSolvedApprox as sampleList. lastSolvedApprox: Parameter(s) corresponding to last computed approximate solution(s) as parameterList. Q: Numpy 1D vector containing complex coefficients of approximant denominator. P: Numpy 2D vector whose columns are FE dofs of coefficients of approximant numerator. """ def __init__(self, HFEngine:HFEng, mu0 : paramVal = None, approxParameters : DictAny = {}, approx_state : bool = False, verbosity : int = 10, timestamp : bool = True): self._preInit() self._addParametersToList(["polybasis", "M", "N", "polydegreetype", "radialDirectionalWeights", "nNearestNeighbor", "interpRcond", "robustTol", "correctorForce", "correctorTol", "correctorMaxIter"], ["MONOMIAL", 0, 0, "TOTAL", 1, -1, -1, 0, False, 0., 1]) super().__init__(HFEngine = HFEngine, mu0 = mu0, approxParameters = approxParameters, approx_state = approx_state, verbosity = verbosity, timestamp = timestamp) - self.catchInstability = False + self.catchInstability = 0 self._postInit() @property def tModelType(self): from rrompy.reduction_methods.trained_model import TrainedModelRational return TrainedModelRational @property def polybasis(self): """Value of polybasis.""" return self._polybasis @polybasis.setter def polybasis(self, polybasis): try: polybasis = polybasis.upper().strip().replace(" ","") if polybasis not in ppb + rbpb + mlspb: raise RROMPyException("Prescribed polybasis not recognized.") self._polybasis = polybasis except: RROMPyWarning(("Prescribed polybasis not recognized. Overriding " "to 'MONOMIAL'.")) self._polybasis = "MONOMIAL" self._approxParameters["polybasis"] = self.polybasis @property def polybasis0(self): if "_" in self.polybasis: return self.polybasis.split("_")[0] return self.polybasis @property def interpRcond(self): """Value of interpRcond.""" return self._interpRcond @interpRcond.setter def interpRcond(self, interpRcond): self._interpRcond = interpRcond self._approxParameters["interpRcond"] = self.interpRcond @property def radialDirectionalWeights(self): """Value of radialDirectionalWeights.""" return self._radialDirectionalWeights @radialDirectionalWeights.setter def radialDirectionalWeights(self, radialDirectionalWeights): self._radialDirectionalWeights = radialDirectionalWeights self._approxParameters["radialDirectionalWeights"] = ( self.radialDirectionalWeights) @property def nNearestNeighbor(self): """Value of nNearestNeighbor.""" return self._nNearestNeighbor @nNearestNeighbor.setter def nNearestNeighbor(self, nNearestNeighbor): self._nNearestNeighbor = nNearestNeighbor self._approxParameters["nNearestNeighbor"] = self.nNearestNeighbor @property def M(self): """Value of M.""" return self._M @M.setter def M(self, M): if M < 0: raise RROMPyException("M must be non-negative.") self._M = M self._approxParameters["M"] = self.M @property def N(self): """Value of N.""" return self._N @N.setter def N(self, N): if N < 0: raise RROMPyException("N must be non-negative.") self._N = N self._approxParameters["N"] = self.N @property def polydegreetype(self): """Value of polydegreetype.""" return self._polydegreetype @polydegreetype.setter def polydegreetype(self, polydegreetype): try: polydegreetype = polydegreetype.upper().strip().replace(" ","") if polydegreetype not in ["TOTAL", "FULL"]: raise RROMPyException(("Prescribed polydegreetype not " "recognized.")) self._polydegreetype = polydegreetype except: RROMPyWarning(("Prescribed polydegreetype not recognized. " "Overriding to 'TOTAL'.")) self._polydegreetype = "TOTAL" self._approxParameters["polydegreetype"] = self.polydegreetype @property def robustTol(self): """Value of tolerance for robust rational denominator management.""" return self._robustTol @robustTol.setter def robustTol(self, robustTol): if robustTol < 0.: RROMPyWarning(("Overriding prescribed negative robustness " "tolerance to 0.")) robustTol = 0. self._robustTol = robustTol self._approxParameters["robustTol"] = self.robustTol @property def correctorForce(self): """Value of correctorForce.""" return self._correctorForce @correctorForce.setter def correctorForce(self, correctorForce): self._correctorForce = correctorForce self._approxParameters["correctorForce"] = self.correctorForce @property def correctorTol(self): """Value of correctorTol.""" return self._correctorTol @correctorTol.setter def correctorTol(self, correctorTol): if correctorTol < 0. or (correctorTol > 0. and self.npar > 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.npar > 1): RROMPyWarning(("Overriding prescribed max number of corrector " "iterations to 1.")) correctorMaxIter = 1 self._correctorMaxIter = correctorMaxIter self._approxParameters["correctorMaxIter"] = self.correctorMaxIter def resetSamples(self): """Reset samples.""" super().resetSamples() self._musUniqueCN = None self._derIdxs = None self._reorder = None def _setupInterpolationIndices(self): """Setup parameters for polyvander.""" if self._musUniqueCN is None or len(self._reorder) != len(self.mus): self._musUniqueCN, musIdxsTo, musIdxs, musCount = ( self.trainedModel.centerNormalize(self.mus).unique( return_index = True, return_inverse = True, return_counts = True)) self._musUnique = self.mus[musIdxsTo] self._derIdxs = [None] * len(self._musUniqueCN) self._reorder = np.empty(len(musIdxs), dtype = int) filled = 0 for j, cnt in enumerate(musCount): self._derIdxs[j] = nextDerivativeIndices([], self.mus.shape[1], cnt) jIdx = np.nonzero(musIdxs == j)[0] self._reorder[jIdx] = np.arange(filled, filled + cnt) filled += cnt def _setupDenominator(self): """Compute rational denominator.""" RROMPyAssert(self._mode, message = "Cannot setup denominator.") vbMng(self, "INIT", "Starting computation of denominator.", 7) while self.N > 0: invD, fitinv = self._computeInterpolantInverseBlocks() idxSamplesEff = list(range(self.S)) if self.POD: ev, eV = self.findeveVGQR( self.samplingEngine.RPOD[:, idxSamplesEff], invD) else: ev, eV = self.findeveVGExplicit( self.samplingEngine.samples(idxSamplesEff), invD) nevBad = checkRobustTolerance(ev, self.robustTol) if nevBad <= 1: break - if self.catchInstability: + if self.catchInstability > 0: raise RROMPyException(("Instability in denominator " "computation: eigenproblem is poorly " - "conditioned.")) + "conditioned."), + self.catchInstability == 1) RROMPyWarning(("Smallest {} eigenvalues below tolerance. Reducing " "N by 1.").format(nevBad)) self.N = self.N - 1 if self.N <= 0: self._N = 0 eV = np.ones((1, 1)) q = PI() q.npar = self.npar q.polybasis = self.polybasis0 if self.polydegreetype == "TOTAL": q.coeffs = degreeTotalToFull(tuple([self.N + 1] * self.npar), self.npar, eV[:, 0]) else: q.coeffs = eV[:, 0].reshape([self.N + 1] * self.npar) vbMng(self, "DEL", "Done computing denominator.", 7) return q, fitinv def _setupNumerator(self): """Compute rational numerator.""" RROMPyAssert(self._mode, message = "Cannot setup numerator.") vbMng(self, "INIT", "Starting computation of numerator.", 7) self._setupInterpolationIndices() Qevaldiag = polyTimesTable(self.trainedModel.data.Q, self._musUniqueCN, self._reorder, self._derIdxs, np.power(self.scaleFactor, -1.)) if self.POD: Qevaldiag = Qevaldiag.dot(self.samplingEngine.RPOD.T) cfun = totalDegreeN if self.polydegreetype == "TOTAL" else fullDegreeN M = copy(self.M) while self.S < cfun(M, self.npar): M -= 1 if M < self.M: RROMPyWarning(("M too large compared to S. Reducing M by " "{}").format(self.M - M)) self.M = M while self.M >= 0: pParRest = [self.verbosity >= 5, self.polydegreetype == "TOTAL", {"derIdxs": self._derIdxs, "reorder": self._reorder, "scl": np.power(self.scaleFactor, -1.)}] if self.polybasis in ppb: p = PI() else: pParRest = [self.radialDirectionalWeights] + pParRest pParRest[-1]["nNearestNeighbor"] = self.nNearestNeighbor p = RBI() if self.polybasis in rbpb else MLSI() if self.polybasis in ppb + rbpb: pParRest += [{"rcond": self.interpRcond}] wellCond, msg = p.setupByInterpolation(self._musUniqueCN, Qevaldiag, self.M, self.polybasis, *pParRest) vbMng(self, "MAIN", msg, 5) if wellCond: break - if self.catchInstability: + if self.catchInstability > 0: raise RROMPyException(("Instability in numerator computation: " - "polyfit is poorly conditioned.")) + "polyfit is poorly conditioned."), + self.catchInstability == 1) RROMPyWarning("Polyfit is poorly conditioned. Reducing M by 1.") self.M = self.M - 1 if self.M < 0: raise RROMPyException(("Instability in computation of numerator. " "Aborting.")) vbMng(self, "DEL", "Done computing numerator.", 7) return p def setupApprox(self): """Compute rational interpolant.""" if self.checkComputedApprox(): return RROMPyAssert(self._mode, message = "Cannot setup approximant.") vbMng(self, "INIT", "Setting up {}.". format(self.name()), 5) self.computeSnapshots() pMat = self.samplingEngine.samples.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, "projMat": pMatEff, "scaleFactor": self.scaleFactor, "rescalingExp": self.HFEngine.rescalingExp} 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) self._iterCorrector() self.trainedModel.data.approxParameters = copy(self.approxParameters) vbMng(self, "DEL", "Done setting up approximant.", 5) def _iterCorrector(self): self._N0 = self.N if self.correctorTol > 0. and (self.correctorMaxIter > 1 or self.correctorForce): vbMng(self, "INIT", "Starting correction iterations.", 7) self._Qhat = PI() self._Qhat.npar = self.npar self._Qhat.polybasis = "MONOMIAL" self._Qhat.coeffs = np.ones(1, dtype = np.complex) if self.POD: self._RPODOld = copy(self.samplingEngine.RPOD) else: self._samplesOld = copy(self.samplingEngine.samples) for j in range(self.correctorMaxIter): if self.N > 0: Q = self._setupDenominator()[0] else: Q = PI() Q.coeffs = np.ones((1,) * self.npar, dtype = np.complex) Q.npar = self.npar Q.polybasis = self.polybasis self.trainedModel.data.Q = Q self.trainedModel.data.P = self._setupNumerator() self._applyCorrector(j) if self.N <= 0: break self._N = self._N0 del self._N0 if self.correctorTol <= 0. or (self.correctorMaxIter <= 1 and not self.correctorForce): return if self.POD: self.samplingEngine.RPOD = self._RPODOld del self._RPODOld else: self.samplingEngine.samples = self._samplesOld del self._samplesOld if self.correctorForce: QOld, QOldBasis = [1.], "MONOMIAL" else: QOld, QOldBasis = Q.coeffs, self.polybasis Q = polyTimes(self._Qhat.coeffs, QOld, Pbasis = self._Qhat.polybasis, Qbasis = QOldBasis, Rbasis = self.polybasis) del self._Qhat gamma = np.linalg.norm(Q) self.trainedModel.data.Q.coeffs = np.pad(Q, (0, self._N - len(Q) + 1), "constant") / gamma if self.correctorForce: self.trainedModel.data.P = self._setupNumerator() else: self.trainedModel.data.P.coeffs /= gamma vbMng(self, "DEL", "Terminated correction iterations.", 7) def _applyCorrector(self, j:int): if self.correctorTol <= 0. or (j >= self.correctorMaxIter - 1 and not self.correctorForce): self._N = 0 return cfs, pls, _ = rational2heaviside(self.trainedModel.data.P, self.trainedModel.data.Q) cfs = cfs[: self.N] if self.POD: resEff = np.linalg.norm(cfs, axis = 1) else: resEff = self.HFEngine.norm(self.samplingEngine.samples.dot(cfs.T), is_state = self.approx_state) resEff /= np.max(np.hstack((np.ones((self.N, 1)), np.abs(pls).reshape(-1, 1))), axis = 1) resEff /= np.mean(resEff) idxKeep = np.where(resEff >= self.correctorTol)[0] vbMng(self, "MAIN", ("Correction iteration no. {}: {} out of {} residuals satisfy " "tolerance.").format(j + 1, len(idxKeep), self.N), 10) self._N -= len(idxKeep) if self.N <= 0 and not self.correctorForce: return for i in idxKeep: self._Qhat.coeffs = polyTimes(self._Qhat.coeffs, [- pls[i], 1.], Pbasis = self._Qhat.polybasis, Rbasis = self._Qhat.polybasis) self._Qhat.coeffs /= np.linalg.norm(self._Qhat.coeffs) if self.N <= 0: return vbMng(self, "MAIN", ("Removing poles at (normalized positions): " "{}.").format(pls[resEff < self.correctorTol]), 65) That = polyTimesTable(self._Qhat, self._musUniqueCN, self._reorder, self._derIdxs, np.power(self.scaleFactor, -1.)).T if self.POD: self.samplingEngine.RPOD = self._RPODOld.dot(That) else: self.samplingEngine.samples = self._samplesOld.dot(That) def _computeInterpolantInverseBlocks(self) -> Tuple[List[Np2D], Np2D]: """ Compute inverse factors for minimal interpolant target functional. """ RROMPyAssert(self._mode, message = "Cannot solve eigenvalue problem.") self._setupInterpolationIndices() cfun = totalDegreeN if self.polydegreetype == "TOTAL" else fullDegreeN N = copy(self.N) while self.S < cfun(N, self.npar): N -= 1 if N < self.N: RROMPyWarning(("N too large compared to S. Reducing N by " "{}").format(self.N - N)) self.N = N TEGen = pvTP if self.polydegreetype == "TOTAL" else pvP TEGenPar = [self.polybasis0, self._derIdxs, self._reorder, np.power(self.scaleFactor, -1.)] while self.N >= 0: if self.polydegreetype == "TOTAL": Neff = self.N idxsB = totalDegreeMaxMask(self.N, self.npar) else: #if self.polydegreetype == "FULL": Neff = [self.N] * self.npar idxsB = fullDegreeMaxMask(self.N, self.npar) TE = TEGen(self._musUniqueCN, Neff, *TEGenPar) fitOut = customPInv(TE, rcond = self.interpRcond, full = True) vbMng(self, "MAIN", ("Fitting {} samples with degree {} through {}... " "Conditioning of pseudoinverse system: {:.4e}.").format( TE.shape[0], self.N, polyfitname(self.polybasis0), fitOut[1][1][0] / fitOut[1][1][-1]), 5) if fitOut[1][0] == TE.shape[1]: fitinv = fitOut[0][idxsB, :] break - if self.catchInstability: + if self.catchInstability > 0: raise RROMPyException(("Instability in denominator " "computation: polyfit is poorly " - "conditioned.")) + "conditioned."), + self.catchInstability == 1) RROMPyWarning("Polyfit is poorly conditioned. Reducing N by 1.") self.N = self.N - 1 invD = vanderInvTable(fitinv, idxsB, self._reorder, self._derIdxs) for k in range(len(invD)): invD[k] = dot(invD[k], TE) return invD, fitinv def findeveVGExplicit(self, sampleE:sampList, invD:List[Np2D]) -> Tuple[Np1D, Np2D]: """ Compute explicitly eigenvalues and eigenvectors of rational denominator matrix. """ RROMPyAssert(self._mode, message = "Cannot solve eigenvalue problem.") nEnd = invD[0].shape[1] eWidth = len(invD) vbMng(self, "INIT", "Building gramian matrix.", 10) gramian = self.HFEngine.innerProduct(sampleE, sampleE, is_state = self.approx_state) G = np.zeros((nEnd, nEnd), dtype = np.complex) for k in range(eWidth): G += dot(dot(gramian, invD[k]).T, invD[k].conj()).T vbMng(self, "DEL", "Done building gramian.", 10) vbMng(self, "INIT", "Solving eigenvalue problem for gramian matrix.", 7) ev, eV = np.linalg.eigh(G) vbMng(self, "MAIN", ("Solved eigenvalue problem of size {} with condition number " "{:.4e}.").format(nEnd, ev[-1] / ev[0]), 5) vbMng(self, "DEL", "Done solving eigenvalue problem.", 7) return ev, eV def findeveVGQR(self, RPODE:Np2D, invD:List[Np2D]) -> Tuple[Np1D, Np2D]: """ Compute eigenvalues and eigenvectors of rational denominator matrix through SVD of R factor. """ RROMPyAssert(self._mode, message = "Cannot solve eigenvalue problem.") nEnd = invD[0].shape[1] S = RPODE.shape[0] eWidth = len(invD) vbMng(self, "INIT", "Building half-gramian matrix stack.", 10) Rstack = np.zeros((S * eWidth, nEnd), dtype = np.complex) for k in range(eWidth): Rstack[k * S : (k + 1) * S, :] = dot(RPODE, invD[k]) vbMng(self, "DEL", "Done building half-gramian.", 10) vbMng(self, "INIT", "Solving svd for square root of gramian matrix.", 7) _, s, eV = np.linalg.svd(Rstack, full_matrices = False) ev = s[::-1] eV = eV[::-1, :].T.conj() vbMng(self, "MAIN", ("Solved svd problem of size {} x {} with condition number " "{:.4e}.").format(*Rstack.shape, s[0] / s[-1]), 5) vbMng(self, "DEL", "Done solving svd.", 7) return ev, eV def getResidues(self, *args, **kwargs) -> Np1D: """ Obtain approximant residues. Returns: Matrix with residues as columns. """ return self.trainedModel.getResidues(*args, **kwargs) diff --git a/rrompy/reduction_methods/standard/rational_moving_least_squares.py b/rrompy/reduction_methods/standard/rational_moving_least_squares.py index cb0377c..ecf846b 100644 --- a/rrompy/reduction_methods/standard/rational_moving_least_squares.py +++ b/rrompy/reduction_methods/standard/rational_moving_least_squares.py @@ -1,340 +1,340 @@ # 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 .rational_interpolant import RationalInterpolant from rrompy.utilities.poly_fitting.polynomial import (polybases as ppb, polyvander as pvP, polyvanderTotal as pvTP) from rrompy.utilities.base.types import Np2D, HFEng, DictAny, paramVal from rrompy.utilities.base import verbosityManager as vbMng from rrompy.utilities.numerical import (fullDegreeMaxMask, totalDegreeMaxMask, dot) from rrompy.utilities.exception_manager import (RROMPyException, RROMPyAssert, RROMPyWarning) __all__ = ['RationalMovingLeastSquares'] class RationalMovingLeastSquares(RationalInterpolant): """ ROM rational moving LS interpolant computation for parametric problems. Args: HFEngine: HF problem solver. mu0(optional): Default parameter. Defaults to 0. approxParameters(optional): Dictionary containing values for main parameters of approximant. Recognized keys are: - 'POD': whether to compute POD of snapshots; defaults to True; - 'S': total number of samples current approximant relies upon; - 'sampler': sample point generator; - 'polybasis': type of polynomial basis for interpolation; defaults to 'MONOMIAL'; - 'M': degree of rational interpolant numerator; defaults to 0; - 'N': degree of rational interpolant denominator; defaults to 0; - 'polydegreetype': type of polynomial degree; defaults to 'TOTAL'; - 'radialBasis': numerator radial basis type; defaults to 'GAUSSIAN'; - 'radialDirectionalWeights': radial basis weights for interpolant numerator; defaults to 0, i.e. identity; - 'nNearestNeighbor': number of nearest neighbors considered in numerator if radialBasis allows; defaults to -1; - 'radialBasisDen': denominator radial basis type; defaults to 'GAUSSIAN'; - 'radialDirectionalWeightsDen': radial basis weights for interpolant denominator; defaults to 0, i.e. identity; - 'nNearestNeighborDen': number of nearest neighbors considered in denominator if radialBasisDen allows; defaults to -1; - 'interpRcond': tolerance for interpolation; defaults to None; - 'robustTol': tolerance for robust rational denominator management; defaults to 0. Defaults to empty dict. approx_state(optional): Whether to approximate state. Defaults to False. verbosity(optional): Verbosity level. Defaults to 10. Attributes: HFEngine: HF problem solver. mu0: Default parameter. mus: Array of snapshot parameters. approxParameters: Dictionary containing values for main parameters of approximant. Recognized keys are in parameterList. parameterListSoft: Recognized keys of soft approximant parameters: - 'POD': whether to compute POD of snapshots; - 'polybasis': type of polynomial basis for interpolation; - 'M': degree of rational interpolant numerator; - 'N': degree of rational interpolant denominator; - 'polydegreetype': type of polynomial degree; - 'radialBasis': numerator radial basis type; - 'radialDirectionalWeights': radial basis weights for interpolant numerator; - 'nNearestNeighbor': number of nearest neighbors considered in numerator if radialBasis allows; - 'radialBasisDen': denominator radial basis type; - 'radialDirectionalWeightsDen': radial basis weights for interpolant denominator; - 'nNearestNeighborDen': number of nearest neighbors considered in denominator if radialBasisDen allows; - 'interpRcond': tolerance for interpolation via numpy.polyfit; - 'robustTol': tolerance for robust rational denominator management. parameterListCritical: Recognized keys of critical approximant parameters: - 'S': total number of samples current approximant relies upon; - 'sampler': sample point generator. approx_state: Whether to approximate state. verbosity: Verbosity level. POD: Whether to compute POD of snapshots. S: Number of solution snapshots over which current approximant is based upon. sampler: Sample point generator. polybasis: type of polynomial basis for interpolation. M: Numerator degree of approximant. N: Denominator degree of approximant. polydegreetype: Type of polynomial degree. radialBasis: Numerator radial basis type. radialDirectionalWeights: Radial basis weights for interpolant numerator. nNearestNeighbor: Number of nearest neighbors considered in numerator if radialBasis allows. radialBasisDen: Denominator radial basis type. radialDirectionalWeightsDen: Radial basis weights for interpolant denominator. nNearestNeighborDen: Number of nearest neighbors considered in denominator if radialBasisDen allows. interpRcond: Tolerance for interpolation via numpy.polyfit. robustTol: Tolerance for robust rational denominator management. muBounds: list of bounds for parameter values. samplingEngine: Sampling engine. uHF: High fidelity solution(s) with parameter(s) lastSolvedHF as sampleList. lastSolvedHF: Parameter(s) corresponding to last computed high fidelity solution(s) as parameterList. uApproxReduced: Reduced approximate solution(s) with parameter(s) lastSolvedApprox as sampleList. lastSolvedApproxReduced: Parameter(s) corresponding to last computed reduced approximate solution(s) as parameterList. uApprox: Approximate solution(s) with parameter(s) lastSolvedApprox as sampleList. lastSolvedApprox: Parameter(s) corresponding to last computed approximate solution(s) as parameterList. Q: Numpy 1D vector containing complex coefficients of approximant denominator. P: Numpy 2D vector whose columns are FE dofs of coefficients of approximant numerator. """ def __init__(self, HFEngine:HFEng, mu0 : paramVal = None, approxParameters : DictAny = {}, approx_state : bool = False, verbosity : int = 10, timestamp : bool = True): self._preInit() self._addParametersToList(["radialBasis", "radialBasisDen", "radialDirectionalWeightsDen", "nNearestNeighborDen"], ["GAUSSIAN", "GAUSSIAN", 1, -1], toBeExcluded = ["correctorForce", "correctorTol", "correctorMaxIter"]) super().__init__(HFEngine = HFEngine, mu0 = mu0, approxParameters = approxParameters, approx_state = approx_state, verbosity = verbosity, timestamp = timestamp) - self.catchInstability = False + self.catchInstability = 0 self._postInit() @property def correctorForce(self): """Value of correctorForce.""" return False @correctorForce.setter def correctorForce(self, correctorForce): RROMPyWarning(("correctorForce is used just to simplify inheritance, " "and its value cannot be changed from False.")) @property def correctorTol(self): """Value of correctorTol.""" return 0. @correctorTol.setter def correctorTol(self, correctorTol): RROMPyWarning(("correctorTol is used just to simplify inheritance, " "and its value cannot be changed from 0.")) @property def correctorMaxIter(self): """Value of correctorMaxIter.""" return 1 @correctorMaxIter.setter def correctorMaxIter(self, correctorMaxIter): RROMPyWarning(("correctorMaxIter is used just to simplify " "inheritance, and its value cannot be changed from 1.")) @property def tModelType(self): from rrompy.reduction_methods.trained_model import \ TrainedModelRationalMLS return TrainedModelRationalMLS @property def polybasis(self): """Value of polybasis.""" return self._polybasis @polybasis.setter def polybasis(self, polybasis): try: polybasis = polybasis.upper().strip().replace(" ","") if polybasis not in ppb: raise RROMPyException("Prescribed polybasis not recognized.") self._polybasis = polybasis except: RROMPyWarning(("Prescribed polybasis not recognized. Overriding " "to 'MONOMIAL'.")) self._polybasis = "MONOMIAL" self._approxParameters["polybasis"] = self.polybasis @property def radialBasis(self): """Value of radialBasis.""" return self._radialBasis @radialBasis.setter def radialBasis(self, radialBasis): self._radialBasis = radialBasis self._approxParameters["radialBasis"] = self.radialBasis @property def radialBasisDen(self): """Value of radialBasisDen.""" return self._radialBasisDen @radialBasisDen.setter def radialBasisDen(self, radialBasisDen): self._radialBasisDen = radialBasisDen self._approxParameters["radialBasisDen"] = self.radialBasisDen @property def radialDirectionalWeightsDen(self): """Value of radialDirectionalWeightsDen.""" return self._radialDirectionalWeightsDen @radialDirectionalWeightsDen.setter def radialDirectionalWeightsDen(self, radialDirectionalWeightsDen): self._radialDirectionalWeightsDen = radialDirectionalWeightsDen self._approxParameters["radialDirectionalWeightsDen"] = ( self.radialDirectionalWeightsDen) @property def nNearestNeighborDen(self): """Value of nNearestNeighborDen.""" return self._nNearestNeighborDen @nNearestNeighborDen.setter def nNearestNeighborDen(self, nNearestNeighborDen): self._nNearestNeighborDen = nNearestNeighborDen self._approxParameters["nNearestNeighborDen"] = ( self.nNearestNeighborDen) def _setupDenominator(self) -> Np2D: """Compute rational denominator.""" RROMPyAssert(self._mode, message = "Cannot setup denominator.") vbMng(self, "INIT", "Starting computation of denominator-related blocks.", 7) self._setupInterpolationIndices() if self.polydegreetype == "TOTAL": TN = pvTP(self._musUniqueCN, self.N, self.polybasis0, self._derIdxs, self._reorder, scl = np.power(self.scaleFactor, -1.)) else: #if self.polydegreetype == "FULL": TN = pvP(self._musUniqueCN, [self.N] * self.npar, self.polybasis0, self._derIdxs, self._reorder, scl = np.power(self.scaleFactor, -1.)) TNTen = np.zeros((self.S, self.S, TN.shape[1]), dtype = TN.dtype) TNTen[np.arange(self.S), np.arange(self.S)] = TN if self.POD: TNTen = dot(self.samplingEngine.RPOD, TNTen) vbMng(self, "DEL", "Done computing denominator-related blocks.", 7) return TN, TNTen def _setupNumerator(self) -> Np2D: """Compute rational numerator.""" RROMPyAssert(self._mode, message = "Cannot setup numerator.") vbMng(self, "INIT", "Starting computation of denominator-related blocks.", 7) self._setupInterpolationIndices() if self.polydegreetype == "TOTAL": TM = pvTP(self._musUniqueCN, self.M, self.polybasis0, self._derIdxs, self._reorder, scl = np.power(self.scaleFactor, -1.)) else: #if self.polydegreetype == "FULL": TM = pvP(self._musUniqueCN, [self.M] * self.npar, self.polybasis0, self._derIdxs, self._reorder, scl = np.power(self.scaleFactor, -1.)) vbMng(self, "DEL", "Done computing denominator-related blocks.", 7) return TM def setupApprox(self): """ Compute rational interpolant. SVD-based robust eigenvalue management. """ if self.checkComputedApprox(): return RROMPyAssert(self._mode, message = "Cannot setup approximant.") vbMng(self, "INIT", "Setting up {}.". format(self.name()), 5) self.computeSnapshots() pMat = self.samplingEngine.samples.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, "projMat": pMatEff, "scaleFactor": self.scaleFactor, "rescalingExp": self.HFEngine.rescalingExp} data = self.initializeModelData(datadict)[0] data.POD = self.POD data.polybasis = self.polybasis data.polydegreetype = self.polydegreetype data.radialBasis = self.radialBasis data.radialWeights = self.radialDirectionalWeights data.nNearestNeighbor = self.nNearestNeighbor data.radialBasisDen = self.radialBasisDen data.radialWeightsDen = self.radialDirectionalWeightsDen data.nNearestNeighborDen = self.nNearestNeighborDen data.interpRcond = self.interpRcond self.trainedModel.data = data else: self.trainedModel = self.trainedModel self.trainedModel.data.projMat = copy(pMatEff) if not self.POD: self.trainedModel.data.gramian = self.HFEngine.innerProduct( self.samplingEngine.samples, self.samplingEngine.samples, is_state = self.approx_state) self.trainedModel.data.mus = copy(self.mus) self.trainedModel.data.M = self.M self.trainedModel.data.N = self.N QVan, self.trainedModel.data.QBlocks = self._setupDenominator() self.trainedModel.data.PVan = self._setupNumerator() if self.polydegreetype == "TOTAL": degreeMaxMask = totalDegreeMaxMask else: #if self.polydegreetype == "FULL": degreeMaxMask = fullDegreeMaxMask if self.N > self.M: self.trainedModel.data.QVan = QVan self.trainedModel.data.domQIdxs = degreeMaxMask(self.N, self.npar) else: self.trainedModel.data.domQIdxs = degreeMaxMask(self.M, self.npar) self.trainedModel.data.approxParameters = copy(self.approxParameters) vbMng(self, "DEL", "Done setting up approximant.", 5) diff --git a/rrompy/reduction_methods/standard/rational_pade.py b/rrompy/reduction_methods/standard/rational_pade.py index 9c84461..200e235 100644 --- a/rrompy/reduction_methods/standard/rational_pade.py +++ b/rrompy/reduction_methods/standard/rational_pade.py @@ -1,282 +1,285 @@ # Copyright (C) 2018 by the RROMPy authors # # This file is part of RROMPy. # # RROMPy is free software: you can redistribute it and/or modify # it under the terms of the GNU Lesser General Public License as published by # the Free Software Foundation, either version 3 of the License, or # (at your option) any later version. # # RROMPy is distributed in the hope that it will be useful, # but WITHOUT ANY WARRANTY; without even the implied warranty of # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the # GNU Lesser General Public License for more details. # # You should have received a copy of the GNU Lesser General Public License # along with RROMPy. If not, see . # from copy import deepcopy as copy import numpy as np from rrompy.reduction_methods.base import checkRobustTolerance from .rational_interpolant import RationalInterpolant from rrompy.utilities.poly_fitting.polynomial import ( polybases as ppb, polyfitname, polyvander as pvP, polyvanderTotal as pvTP, polyTimesTable, vanderInvTable, PolynomialInterpolator as PI) from rrompy.utilities.poly_fitting.radial_basis import (polybases as rbpb, RadialBasisInterpolator as RBI) from rrompy.utilities.poly_fitting.moving_least_squares import ( MovingLeastSquaresInterpolator as MLSI) from rrompy.utilities.base.types import Np2D, Tuple, List from rrompy.utilities.base import verbosityManager as vbMng from rrompy.utilities.numerical import (customPInv, dot, fullDegreeN, totalDegreeN, degreeTotalToFull, fullDegreeMaxMask, totalDegreeMaxMask) from rrompy.utilities.exception_manager import (RROMPyException, RROMPyAssert, RROMPyWarning) __all__ = ['RationalPade'] class RationalPade(RationalInterpolant): """ ROM rational Pade' computation for parametric problems. Args: HFEngine: HF problem solver. mu0(optional): Default parameter. Defaults to 0. approxParameters(optional): Dictionary containing values for main parameters of approximant. Recognized keys are: - 'POD': whether to compute POD of snapshots; defaults to True; - 'S': total number of samples current approximant relies upon; - 'sampler': sample point generator; - 'polybasis': type of polynomial basis for interpolation; defaults to 'MONOMIAL'; - 'M': degree of rational interpolant numerator; defaults to 0; - 'N': degree of rational interpolant denominator; defaults to 0; - 'polydegreetype': type of polynomial degree; defaults to 'TOTAL'; - 'radialDirectionalWeights': radial basis weights for interpolant numerator; defaults to 0, i.e. identity; - 'nNearestNeighbor': mumber of nearest neighbors considered in numerator if polybasis allows; defaults to -1; - 'interpRcond': tolerance for 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. mus: Array of snapshot parameters. approxParameters: Dictionary containing values for main parameters of approximant. Recognized keys are in parameterList. parameterListSoft: Recognized keys of soft approximant parameters: - 'POD': whether to compute POD of snapshots; - 'polybasis': type of polynomial basis for interpolation; - 'M': degree of rational interpolant numerator; - 'N': degree of rational interpolant denominator; - 'polydegreetype': type of polynomial degree; - 'radialDirectionalWeights': radial basis weights for interpolant numerator; - 'nNearestNeighbor': mumber of nearest neighbors considered in numerator if polybasis allows; - 'interpRcond': tolerance for interpolation via numpy.polyfit; - '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 samples current approximant relies upon; - 'sampler': sample point generator. approx_state: Whether to approximate state. verbosity: Verbosity level. POD: Whether to compute POD of snapshots. S: Number of solution snapshots over which current approximant is based upon. sampler: Sample point generator. polybasis: type of polynomial basis for interpolation. M: Numerator degree of approximant. N: Denominator degree of approximant. polydegreetype: Type of polynomial degree. radialDirectionalWeights: Radial basis weights for interpolant numerator. nNearestNeighbor: Number of nearest neighbors considered in numerator if polybasis allows. interpRcond: Tolerance for interpolation via numpy.polyfit. 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 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 _setupInterpolationIndices(self): """Setup parameters for polyvander.""" super()._setupInterpolationIndices() if len(self._musUniqueCN) > 1: raise RROMPyException(("Cannot apply centered-like method with " "more than one distinct sample point.")) def _setupDenominator(self): """Compute rational denominator.""" RROMPyAssert(self._mode, message = "Cannot setup denominator.") vbMng(self, "INIT", "Starting computation of denominator.", 7) while self.N > 0: invD, fitinv = self._computeInterpolantInverseBlocks() if self.polydegreetype == "TOTAL": Seff = totalDegreeN(self.N, self.npar) else: Seff = fullDegreeN(self.N, self.npar) idxSamplesEff = list(range(self.S - Seff, self.S)) if self.POD: ev, eV = self.findeveVGQR( self.samplingEngine.RPOD[:, idxSamplesEff], invD) else: ev, eV = self.findeveVGExplicit( self.samplingEngine.samples(idxSamplesEff), invD) nevBad = checkRobustTolerance(ev, self.robustTol) if nevBad <= 1: break - if self.catchInstability: + if self.catchInstability > 0: raise RROMPyException(("Instability in denominator " "computation: eigenproblem is poorly " - "conditioned.")) + "conditioned."), + self.catchInstability == 1) RROMPyWarning(("Smallest {} eigenvalues below tolerance. Reducing " "N by 1.").format(nevBad)) self.N = self.N - 1 if self.N <= 0: self._N = 0 eV = np.ones((1, 1)) q = PI() q.npar = self.npar q.polybasis = self.polybasis0 if self.polydegreetype == "TOTAL": q.coeffs = degreeTotalToFull(tuple([self.N + 1] * self.npar), self.npar, eV[:, 0]) else: q.coeffs = eV[:, 0].reshape([self.N + 1] * self.npar) vbMng(self, "DEL", "Done computing denominator.", 7) return q, fitinv def _setupNumerator(self): """Compute rational numerator.""" RROMPyAssert(self._mode, message = "Cannot setup numerator.") vbMng(self, "INIT", "Starting computation of numerator.", 7) self._setupInterpolationIndices() Qevaldiag = polyTimesTable(self.trainedModel.data.Q, self._musUniqueCN, self._reorder, self._derIdxs, np.power(self.scaleFactor, -1.)) if self.POD: Qevaldiag = Qevaldiag.dot(self.samplingEngine.RPOD.T) cfun = totalDegreeN if self.polydegreetype == "TOTAL" else fullDegreeN M = copy(self.M) while self.S < cfun(M, self.npar): M -= 1 if M < self.M: RROMPyWarning(("M too large compared to S. Reducing M by " "{}").format(self.M - M)) self.M = M while self.M >= 0: Seff = cfun(self.M, self.npar) pParRest = [self.verbosity >= 5, self.polydegreetype == "TOTAL", {"derIdxs": [self._derIdxs[0][: Seff]], "reorder": self._reorder[: Seff], "scl": np.power(self.scaleFactor, -1.)}] if self.polybasis in ppb: p = PI() else: pParRest = [self.radialDirectionalWeights] + pParRest pParRest[-1]["nNearestNeighbor"] = self.nNearestNeighbor p = RBI() if self.polybasis in rbpb else MLSI() if self.polybasis in ppb + rbpb: pParRest += [{"rcond": self.interpRcond}] wellCond, msg = p.setupByInterpolation(self._musUniqueCN, Qevaldiag[: Seff, : Seff], self.M, self.polybasis, *pParRest) vbMng(self, "MAIN", msg, 5) if wellCond: break - if self.catchInstability: + if self.catchInstability > 0: raise RROMPyException(("Instability in numerator computation: " - "polyfit is poorly conditioned.")) + "polyfit is poorly conditioned."), + self.catchInstability == 1) RROMPyWarning("Polyfit is poorly conditioned. Reducing M by 1.") self.M = self.M - 1 if self.M < 0: raise RROMPyException(("Instability in computation of numerator. " "Aborting.")) vbMng(self, "DEL", "Done computing numerator.", 7) return p def _computeInterpolantInverseBlocks(self) -> Tuple[List[Np2D], Np2D]: """ Compute inverse factors for minimal interpolant target functional. """ RROMPyAssert(self._mode, message = "Cannot solve eigenvalue problem.") self._setupInterpolationIndices() cfun = totalDegreeN if self.polydegreetype == "TOTAL" else fullDegreeN N = copy(self.N) while self.S < cfun(N, self.npar): N -= 1 if N < self.N: RROMPyWarning(("N too large compared to S. Reducing N by " "{}").format(self.N - N)) self.N = N TEGen = pvTP if self.polydegreetype == "TOTAL" else pvP while self.N >= 0: Seff = cfun(self.N, self.npar) TEGenPar = [self.polybasis0, [self._derIdxs[0][: Seff]], self._reorder[: Seff], np.power(self.scaleFactor, -1.)] if self.polydegreetype == "TOTAL": Neff = self.N idxsB = totalDegreeMaxMask(self.N, self.npar) else: #if self.polydegreetype == "FULL": Neff = [self.N] * self.npar idxsB = fullDegreeMaxMask(self.N, self.npar) TE = TEGen(self._musUniqueCN, Neff, *TEGenPar) fitOut = customPInv(TE, rcond = self.interpRcond, full = True) vbMng(self, "MAIN", ("Fitting {} samples with degree {} through {}... " "Conditioning of pseudoinverse system: {:.4e}.").format( TE.shape[0], self.N, polyfitname(self.polybasis0), fitOut[1][1][0] / fitOut[1][1][-1]), 5) if fitOut[1][0] == TE.shape[1]: fitinv = fitOut[0][idxsB, :] break - if self.catchInstability: + if self.catchInstability > 0: raise RROMPyException(("Instability in denominator " "computation: polyfit is poorly " - "conditioned.")) + "conditioned."), + self.catchInstability == 1) RROMPyWarning("Polyfit is poorly conditioned. Reducing N by 1.") self.N = self.N - 1 invD = vanderInvTable(fitinv, idxsB, self._reorder, self._derIdxs) for k in range(len(invD)): invD[k] = dot(invD[k], TE) return invD, fitinv diff --git a/rrompy/utilities/exception_manager/exception_manager.py b/rrompy/utilities/exception_manager/exception_manager.py index c9c4f1b..86eec92 100644 --- a/rrompy/utilities/exception_manager/exception_manager.py +++ b/rrompy/utilities/exception_manager/exception_manager.py @@ -1,33 +1,33 @@ # 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 . # __all__ = ["RROMPyException"] def purgeVerbosityDepth(): from rrompy.utilities.base.verbosity_depth import verbosityDepth while True: try: verbosityDepth("DEL", "", "", False) except: break class RROMPyException(Exception): - def __init__(self, message): - purgeVerbosityDepth() + def __init__(self, message, purge : bool = True): + if purge: purgeVerbosityDepth() self._msg = message super().__init__(message) diff --git a/tests/reduction_methods_1D/rational_interpolant_greedy_1d.py b/tests/reduction_methods_1D/rational_interpolant_greedy_1d.py index 7e6118e..93cf745 100644 --- a/tests/reduction_methods_1D/rational_interpolant_greedy_1d.py +++ b/tests/reduction_methods_1D/rational_interpolant_greedy_1d.py @@ -1,94 +1,94 @@ # Copyright (C) 2018 by the RROMPy authors # # This file is part of RROMPy. # # RROMPy is free software: you can redistribute it and/or modify # it under the terms of the GNU Lesser General Public License as published by # the Free Software Foundation, either version 3 of the License, or # (at your option) any later version. # # RROMPy is distributed in the hope that it will be useful, # but WITHOUT ANY WARRANTY; without even the implied warranty of # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the # GNU Lesser General Public License for more details. # # You should have received a copy of the GNU Lesser General Public License # along with RROMPy. If not, see . # import numpy as np from matrix_fft import matrixFFT from rrompy.reduction_methods.greedy import RationalInterpolantGreedy as RIG from rrompy.parameter.parameter_sampling import QuadratureSampler as QS def test_lax_tolerance(capsys): mu = 2.25 solver = matrixFFT() params = {"POD": True, "sampler": QS([1.5, 6.5], "UNIFORM"), "S": 4, "polybasis": "CHEBYSHEV", "greedyTol": 1e-2, - "errorEstimatorKind": "DIAGONAL", + "errorEstimatorKind": "LOOK_AHEAD", "trainSetGenerator": QS([1.5, 6.5], "CHEBYSHEV")} - approx = RIG(solver, 4., params, verbosity = 100) + approx = RIG(solver, 4., params, verbosity = 1000) approx.setupApprox() out, err = capsys.readouterr() assert "Done computing snapshots (final snapshot count: 10)." in out assert len(err) == 0 assert np.isclose(approx.normErr(mu)[0], 2.169678e-4, rtol = 1e-1) def test_samples_at_poles(): solver = matrixFFT() params = {"POD": True, "sampler": QS([1.5, 6.5], "UNIFORM"), "S": 4, "nTestPoints": 100, "polybasis": "CHEBYSHEV", "greedyTol": 1e-5, "errorEstimatorKind": "AFFINE", "trainSetGenerator": QS([1.5, 6.5], "CHEBYSHEV")} approx = RIG(solver, 4., params, verbosity = 0) approx.setupApprox() for mu in approx.mus: assert np.isclose(approx.normErr(mu)[0] / (1e-15+approx.normHF(mu)[0]), 0., atol = 1e-4) poles = approx.getPoles() for lambda_ in range(2, 7): assert np.isclose(np.min(np.abs(poles - lambda_)), 0., atol = 1e-3) assert np.isclose(np.min(np.abs(np.array(approx.mus(0)) - lambda_)), 0., atol = 1e-1) def test_maxIter(): solver = matrixFFT() params = {"POD": True, "sampler": QS([1.5, 6.5], "UNIFORM"), "S": 5, "nTestPoints": 500, "polybasis": "CHEBYSHEV", "greedyTol": 1e-6, "maxIter": 10, "errorEstimatorKind": "INTERPOLATORY", "trainSetGenerator": QS([1.5, 6.5], "CHEBYSHEV")} approx = RIG(solver, 4., params, verbosity = 0) approx.input = lambda: "N" approx.setupApprox() assert len(approx.mus) == 10 _, _, maxEst = approx.errorEstimator(approx.muTest, True) assert maxEst > 1e-6 def test_load_copy(capsys): mu = 3. solver = matrixFFT() params = {"POD": True, "sampler": QS([1.5, 6.5], "UNIFORM"), "S": 4, "nTestPoints": 100, "polybasis": "CHEBYSHEV", "greedyTol": 1e-5, "errorEstimatorKind": "AFFINE", "trainSetGenerator": QS([1.5, 6.5], "CHEBYSHEV")} approx1 = RIG(solver, 4., params, verbosity = 100) approx1.setupApprox() err1 = approx1.normErr(mu)[0] out, err = capsys.readouterr() assert "Solving HF model for mu =" in out assert len(err) == 0 approx2 = RIG(solver, 4., params, verbosity = 100) approx2.setTrainedModel(approx1) approx2.setHF(mu, approx1.uHF) err2 = approx2.normErr(mu)[0] out, err = capsys.readouterr() assert "Solving HF model for mu =" not in out assert len(err) == 0 assert np.isclose(err1, err2, rtol = 1e-10)