diff --git a/rrompy/reduction_methods/greedy/rational_interpolant_greedy.py b/rrompy/reduction_methods/greedy/rational_interpolant_greedy.py index 813100e..b72494e 100644 --- a/rrompy/reduction_methods/greedy/rational_interpolant_greedy.py +++ b/rrompy/reduction_methods/greedy/rational_interpolant_greedy.py @@ -1,522 +1,521 @@ # 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: 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 = 2 stable = True if self.N > 0: try: Q = self._setupDenominator()[0] except RROMPyException as RE: RROMPyWarning(RE._msg) 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 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 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/pivoted/rational_interpolant_greedy_pivoted.py b/rrompy/reduction_methods/pivoted/rational_interpolant_greedy_pivoted.py index da79089..678ff49 100644 --- a/rrompy/reduction_methods/pivoted/rational_interpolant_greedy_pivoted.py +++ b/rrompy/reduction_methods/pivoted/rational_interpolant_greedy_pivoted.py @@ -1,453 +1,470 @@ # 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 GenericPivotedApproximant from rrompy.reduction_methods.greedy.rational_interpolant_greedy import ( RationalInterpolantGreedy) from rrompy.reduction_methods.greedy.generic_greedy_approximant import ( pruneSamples) from rrompy.utilities.base.types import Np1D, HFEng, DictAny, ListAny, paramVal from rrompy.utilities.base import verbosityManager as vbMng from rrompy.utilities.numerical import totalDegreeN, dot from rrompy.utilities.poly_fitting.polynomial import polyvander as pv from rrompy.utilities.exception_manager import RROMPyAssert from rrompy.parameter import emptyParameterList, checkParameterList __all__ = ['RationalInterpolantGreedyPivoted'] class RationalInterpolantGreedyPivoted(GenericPivotedApproximant, RationalInterpolantGreedy): """ 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; - 'matchingWeight': weight for pole matching optimization; defaults to 1; - 'cutOffTolerance': tolerance for ignoring parasitic poles; defaults to np.inf; - 'cutOffType': rule for tolerance computation for parasitic poles; defaults to 'MAGNITUDE'; - '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' and 'LEGENDRE'; 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; - '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; - 'errorEstimatorKind': kind of error estimator; available values include 'AFFINE', 'DISCREPANCY', 'INTERPOLATORY', 'LOOK_AHEAD', and 'NONE'; defaults to 'NONE'; - 'MMarginal': degree of marginal interpolant; defaults to 0; - 'polydegreetypeMarginal': type of polynomial degree for marginal; defaults to 'TOTAL'; - 'radialDirectionalWeightsMarginal': radial basis weights for marginal interpolant; defaults to 0, i.e. identity; - 'nNearestNeighborMarginal': number of marginal nearest neighbors considered if polybasisMarginal allows; defaults to -1; - 'interpRcond': tolerance for pivot interpolation; defaults to None; - 'interpRcondMarginal': tolerance for marginal 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. 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; - 'matchingWeight': weight for pole matching optimization; - 'cutOffTolerance': tolerance for ignoring parasitic poles; - 'cutOffType': rule for tolerance computation for parasitic poles; - 'polybasis': type of polynomial basis for pivot interpolation; - 'polybasisMarginal': type of polynomial basis for marginal interpolation; - '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; - 'MMarginal': degree of marginal interpolant; - 'polydegreetypeMarginal': type of polynomial degree for marginal; - 'radialDirectionalWeights': radial basis weights for pivot numerator; - 'radialDirectionalWeightsMarginal': radial basis weights for marginal interpolant; - 'nNearestNeighbor': number of pivot nearest neighbors considered if polybasis allows; - 'nNearestNeighborMarginal': number of marginal nearest neighbors considered if polybasisMarginal allows; - 'interpRcond': tolerance for pivot interpolation; - 'interpRcondMarginal': tolerance for marginal 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. matchingWeight: Weight for pole matching optimization. cutOffTolerance: Tolerance for ignoring parasitic poles. cutOffType: Rule for tolerance computation for 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. polybasisMarginal: Type of polynomial basis for marginal interpolation. 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. errorEstimatorKind: kind of error estimator. MMarginal: Degree of marginal interpolant. polydegreetypeMarginal: Type of polynomial degree for marginal. radialDirectionalWeightsMarginal: Radial basis weights for marginal interpolant. nNearestNeighborMarginal: Number of marginal nearest neighbors considered if polybasisMarginal allows. interpRcond: Tolerance for pivot interpolation. interpRcondMarginal: Tolerance for marginal interpolation. robustTol: Tolerance for robust rational denominator management. muBoundsPivot: list of bounds for pivot parameter values. muBoundsMarginal: list of bounds for marginal parameter values. samplingEngine: Sampling engine. uHF: High fidelity solution(s) with parameter(s) lastSolvedHF as sampleList. lastSolvedHF: Parameter(s) corresponding to last computed high fidelity solution(s) as parameterList. uApproxReduced: Reduced approximate solution(s) with parameter(s) lastSolvedApprox as sampleList. lastSolvedApproxReduced: Parameter(s) corresponding to last computed reduced approximate solution(s) as parameterList. uApprox: Approximate solution(s) with parameter(s) lastSolvedApprox as sampleList. lastSolvedApprox: Parameter(s) corresponding to last computed approximate solution(s) as parameterList. Q: Numpy 1D vector containing complex coefficients of approximant denominator. P: Numpy 2D vector whose columns are FE dofs of coefficients of approximant numerator. """ def __init__(self, HFEngine:HFEng, mu0 : paramVal = None, directionPivot : ListAny = [0], approxParameters : DictAny = {}, approx_state : bool = False, verbosity : int = 10, timestamp : bool = True): self._preInit() self._addParametersToList(toBeExcluded = ["sampler"]) super().__init__(HFEngine = HFEngine, mu0 = mu0, directionPivot = directionPivot, approxParameters = approxParameters, approx_state = approx_state, verbosity = verbosity, timestamp = timestamp) self._postInit() @property def tModelType(self): if hasattr(self, "_temporaryPivot"): return super().tModelType from rrompy.reduction_methods.trained_model import \ TrainedModelPivotedRational return TrainedModelPivotedRational @property def polybasis0(self): if "_" in self.polybasis: return self.polybasis.split("_")[0] return self.polybasis 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((len(self.trainedModel.data.Q.coeffs),) * self.npar, dtype = self.trainedModel.data.Q.coeffs.dtype) Pc = np.zeros((len(self.trainedModel.data.P.coeffs),) * self.npar + (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 not setupOK: err = np.empty(len(mus)) err[:] = np.nan if not return_max: return err return err, 0, 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.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) self.resetSamples() musPivot = self.trainSetGenerator.generatePoints(self.S)[ list(range(self.S))] musPivot = checkParameterList(musPivot, 1)[0] muTestPivot = self.samplerPivot.generatePoints(self.nTestPoints) 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) muTestBase = muTestBase.sort() 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), 2) 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 def _enrichTestSet(self, nTest:int): """Add extra elements to test set.""" RROMPyAssert(self._mode, message = "Cannot enrich test set.") muTestExtra = self.samplerPivot.generatePoints(2 * nTest) muTotal = copy(self.mus) muTotal.append(self.muTest) idxPop = pruneSamples(muTestExtra ** self.HFEngine.rescalingExp[ self.directionPivot[0]], muTotal ** self.HFEngine.rescalingExp[ self.directionPivot[0]], 1e-10 * self.scaleFactor[0]) muTestExtra.pop(idxPop) muTestNew = np.empty((len(self.muTest) + len(muTestExtra), self.muTest.shape[1]), dtype = np.complex) muTestNew[: len(self.muTest)] = self.muTest.data muTestNew[len(self.muTest) :] = muTestExtra.data self.muTest = checkParameterList(muTestNew, self.npar)[0].sort() vbMng(self, "MAIN", "Enriching test set by {} elements.".format(len(muTestExtra)), 5) def setupApprox(self, plotEst : bool = False): """Compute rational interpolant.""" if self.checkComputedApprox(): return True RROMPyAssert(self._mode, message = "Cannot setup approximant.") vbMng(self, "INIT", "Setting up {}.". format(self.name()), 5) self.musMarginal = self.samplerMarginal.generatePoints(self.SMarginal) S0 = copy(self.S) Qs, Ps = [], [] self.computeScaleFactor() self._scaleFactorOldPivot = copy(self.scaleFactor) self.scaleFactor = self.scaleFactorPivot nparEff = self.npar self._temporaryPivot = 1 samplingEngs = [] for j in range(len(self.musMarginal)): self._S = S0 self.musMargLoc = self.musMarginal[j] RationalInterpolantGreedy.setupSampling(self) self.trainedModel = None self.verbosity -= 5 super().setupApprox(plotEst) self.verbosity += 5 samplingEngs += [copy(self.samplingEngine)] Qs = Qs + [copy(self.trainedModel.data.Q)] Ps = Ps + [copy(self.trainedModel.data.P)] del self.musMargLoc # other finalization instructions self.setupSampling() self.samplingEngine.resetHistory(len(self.musMarginal)) for j in range(len(self.musMarginal)): self.samplingEngine.setsample(samplingEngs[j].samples, j, False) self.samplingEngine.mus[j] = copy(samplingEngs[j].mus) self.samplingEngine.nsamples[j] = samplingEngs[j].nsamples self.samplingEngine.postprocessuBulk(j) if j == 0: self._mus = checkParameterList(samplingEngs[j].mus, nparEff)[0] else: self._mus.append(samplingEngs[j].mus) if self.POD: self.samplingEngine.coalesceSamples(self.interpRcondMarginal) postR = self.samplingEngine.RPODCoalesced else: self.samplingEngine.coalesceSamples() postR = np.eye(self.samplingEngine.samplesCoalesced.shape[1]) # update Ps by post-multiplication with suitable matrix idxCurr = 0 for j in range(len(self.musMarginal)): nsj = Ps[j].coeffs.shape[1] Ps[j].coeffs = dot(Ps[j].coeffs, postR[:, idxCurr : idxCurr + nsj].T) idxCurr = idxCurr + nsj self.scaleFactor = self._scaleFactorOldPivot del self._scaleFactorOldPivot, self._temporaryPivot pMat = self.samplingEngine.samplesCoalesced.data pMatEff = dot(self.HFEngine.C, pMat) 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, "directionPivot": self.directionPivot} self.trainedModel.data = self.initializeModelData(datadict)[0] self.trainedModel.data.mus = copy(self.mus) self.trainedModel.data.musMarginal = copy(self.musMarginal) self.trainedModel.data.marginalInterp = self._setupMarginalInterp() self.trainedModel.data.Qs, self.trainedModel.data.Ps = Qs, Ps vbMng(self, "INIT", "Matching poles.", 10) self.trainedModel.initializeFromRational(self.HFEngine, self.matchingWeight, self.POD, self.approx_state) vbMng(self, "DEL", "Done matching poles.", 10) vbMng(self, "INIT", "Recompressing by cut-off.", 10) msg = self.trainedModel.recompressByCutOff([-1., 1.], self.cutOffTolerance, self.cutOffType) vbMng(self, "DEL", "Done recompressing." + msg, 10) self.trainedModel.data.approxParameters = copy(self.approxParameters) vbMng(self, "DEL", "Done setting up approximant.", 5) diff --git a/rrompy/utilities/poly_fitting/polynomial/vander.py b/rrompy/utilities/poly_fitting/polynomial/vander.py index 60bb47d..3f218d9 100644 --- a/rrompy/utilities/poly_fitting/polynomial/vander.py +++ b/rrompy/utilities/poly_fitting/polynomial/vander.py @@ -1,137 +1,138 @@ # 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 .der import polyder from rrompy.utilities.base.types import Np1D, Np2D, List, paramList from rrompy.utilities.numerical import totalDegreeSet from rrompy.parameter import checkParameterList from rrompy.utilities.exception_manager import RROMPyException, RROMPyAssert __all__ = ['polyvander', 'polyvanderTotal'] def firstDerTransition(dim:int, TDirac:List[Np2D], basis:str, scl : Np1D = None) -> Np2D: C_m = np.zeros((dim, len(TDirac), len(TDirac)), dtype = float) for j, Tj in enumerate(TDirac): m, om = [0] * dim, [(0, 0)] * dim for idx in range(dim): m[idx], om[idx] = 1, (0, 1) J_der = polyder(Tj, basis, m, scl) - C_m[idx, :, j] = np.ravel(np.pad(J_der, mode = "constant", - pad_width = om)) + if J_der.size != len(TDirac): + J_der = np.pad(J_der, mode = "constant", pad_width = om) + C_m[idx, :, j] = np.ravel(J_der) m[idx], om[idx] = 0, (0, 0) return C_m def countDerDirections(n:int, base:int, digits:int, idx:int): if digits == 0: return [] dig = n % base return [(idx, dig)] * (dig > 0) + countDerDirections( (n - dig) // base, base, digits - 1, idx + 1) def polyvander(x:paramList, degs:List[int], basis:str, derIdxs : List[List[List[int]]] = None, reorder : List[int] = None, scl : Np1D = None) -> Np2D: """ Compute full Hermite-Vandermonde matrix with specified derivative directions. E.g. assume that we want to obtain the Vandermonde matrix for (value, derx, derx2) at x = [0, 0], (value, dery) at x = [1, 0], (dery, derxy) at x = [0, 0], of degree 3 in x and 1 in y, using Chebyshev polynomials. This can be done by polyvander([[0, 0], [1, 0]], # unique sample points [3, 1], # polynomial degree "chebyshev", # polynomial family [ [[0, 0], [1, 0], [2, 0], [0, 1], [1, 1]], # derivative directions at first point [[0, 0], [0, 1]] # derivative directions at second point ], [0, 1, 2, 5, 6, 3, 4] # reorder indices ) """ if not isinstance(degs, (list,tuple,np.ndarray,)): degs = [degs] dim = len(degs) x = checkParameterList(x, dim)[0] x_un, idx_un = x.unique(return_inverse = True) if len(x_un) < len(x): raise RROMPyException("Sample points must be distinct.") del x_un try: vanderbase = {"CHEBYSHEV" : np.polynomial.chebyshev.chebvander, "LEGENDRE" : np.polynomial.legendre.legvander, "MONOMIAL" : np.polynomial.polynomial.polyvander }[basis.upper()] except: raise RROMPyException("Polynomial basis not recognized.") VanBase = vanderbase(x(0), degs[0]) for j in range(1, dim): VNext = vanderbase(x(j), degs[j]) for jj in range(j): VNext = np.expand_dims(VNext, 1) VanBase = VanBase[..., None] * VNext VanBase = VanBase.reshape((len(x), -1)) if derIdxs is None or VanBase.shape[-1] <= 1: Van = VanBase else: derFlat, idxRep = [], [] for j, derIdx in enumerate(derIdxs): derFlat += derIdx[:] idxRep += [j] * len(derIdx[:]) for j in range(len(derFlat)): if not hasattr(derFlat[j], "__len__"): derFlat[j] = [derFlat[j]] RROMPyAssert(len(derFlat[j]), dim, "Number of dimensions") TDirac = [y.reshape([d + 1 for d in degs]) for y in np.eye(VanBase.shape[-1], dtype = int)] Cs_loc = firstDerTransition(dim, TDirac, basis, scl) Van = np.empty((len(derFlat), VanBase.shape[-1]), dtype = VanBase.dtype) for j in range(len(derFlat)): Van[j, :] = VanBase[idxRep[j], :] for k in range(dim): for der in range(derFlat[j][k]): Van[j, :] = Van[j, :].dot(Cs_loc[k]) / (der + 1) if reorder is not None: Van = Van[reorder, :] return Van def polyvanderTotal(x:paramList, deg:int, basis:str, derIdxs : List[List[List[int]]] = None, reorder : List[int] = None, scl : Np1D = None) -> Np2D: """ Compute full total degree Hermite-Vandermonde matrix with specified derivative directions. """ x = checkParameterList(x)[0] degs = [deg] * x.shape[1] VanBase = polyvander(x, degs, basis, derIdxs, reorder, scl) derIdxs, mask = totalDegreeSet(deg, x.shape[1], return_mask = True) ordIdxs = np.empty(len(derIdxs), dtype = int) derTotal = np.array([np.sum(y) for y in derIdxs]) idxPrev = 0 rangeAux = np.arange(len(derIdxs)) for j in range(deg + 1): idxLocal = rangeAux[derTotal == j][::-1] idxPrev += len(idxLocal) ordIdxs[idxPrev - len(idxLocal) : idxPrev] = idxLocal return VanBase[:, mask][:, ordIdxs] diff --git a/tests/reduction_methods_multiD/pivoted_rational_2d.py b/tests/reduction_methods_multiD/pivoted_rational_2d.py new file mode 100644 index 0000000..647c876 --- /dev/null +++ b/tests/reduction_methods_multiD/pivoted_rational_2d.py @@ -0,0 +1,116 @@ +# 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_random import matrixRandom +from rrompy.reduction_methods.pivoted import ( + RationalInterpolantPivoted as RIP, + RationalInterpolantGreedyPivoted as RIGP) +from rrompy.parameter.parameter_sampling import (QuadratureSampler as QS, + ManualSampler as MS) + +def test_pivoted_uniform(): + mu = [5.05, 7.1] + mu0 = [5., 7.] + solver = matrixRandom() + uh = solver.solve(mu)[0] + params = {"POD": True, "M": 4, "N": 4, "S": 5, "polybasis": "CHEBYSHEV", + "samplerPivot": QS([4.75, 5.25], "CHEBYSHEV"), + "MMarginal": 4, "SMarginal": 5, "polybasisMarginal": "MONOMIAL", + "samplerMarginal": QS([6.75, 7.25], "UNIFORM"), + "matchingWeight": 1.} + approx = RIP(solver, mu0, [0], params, verbosity = 0) + approx.setupApprox() + + uhP1 = approx.getApprox(mu)[0] + errP = approx.getErr(mu)[0] + errNP = approx.normErr(mu)[0] + myerrP = uhP1 - uh + assert np.allclose(np.abs(errP - myerrP), 0., rtol = 1e-3) + assert np.isclose(solver.norm(errP), errNP, rtol = 1e-3) + resP = approx.getRes(mu)[0] + resNP = approx.normRes(mu) + assert np.isclose(solver.norm(resP), resNP, rtol = 1e-3) + assert np.allclose(np.abs(resP - (solver.b(mu) - solver.A(mu).dot(uhP1))), + 0., rtol = 1e-3) + assert np.isclose(errNP / solver.norm(uh), 6.0631706e-04, rtol = 1) + +def test_pivoted_manual_grid(capsys): + mu = [5.05, 7.1] + mu0 = [5., 7.] + solver = matrixRandom() + uh = solver.solve(mu)[0] + params = {"POD": False, "M": 4, "N": 4, "S": 5, "polybasis": "MONOMIAL", + "samplerPivot": MS([4.75, 5.25], np.array([5.])), + "MMarginal": 4, "SMarginal": 5, "polybasisMarginal": "MONOMIAL", + "samplerMarginal": MS([6.75, 7.25], np.linspace(6.75, 7.25, 5)), + "matchingWeight": 1., "robustTol": 1e-6, "interpRcond": 1e-3, + "cutOffTolerance": 1., "cutOffType": "MAGNITUDE"} + approx = RIP(solver, mu0, [0], params, verbosity = 0) + approx.setupApprox() + + uhP1 = approx.getApprox(mu)[0] + errP = approx.getErr(mu)[0] + errNP = approx.normErr(mu)[0] + myerrP = uhP1 - uh + assert np.allclose(np.abs(errP - myerrP), 0., rtol = 1e-3) + assert np.isclose(solver.norm(errP), errNP, rtol = 1e-3) + resP = approx.getRes(mu)[0] + resNP = approx.normRes(mu) + assert np.isclose(solver.norm(resP), resNP, rtol = 1e-3) + assert np.allclose(np.abs(resP - (solver.b(mu) - solver.A(mu).dot(uhP1))), + 0., rtol = 1e-3) + assert np.isclose(errNP / solver.norm(uh), .4763489, rtol = 1) + + out, err = capsys.readouterr() + assert ("poorly conditioned" not in out) + assert len(err) == 0 + +def test_pivoted_greedy(capsys): + mu = [5.05, 7.1] + mu0 = [5., 7.] + solver = matrixRandom() + uh = solver.solve(mu)[0] + params = {"POD": True, "nTestPoints": 100, "greedyTol": 1e-4, + "collinearityTol": 1e8, "errorEstimatorKind": "DISCREPANCY", + "S": 5, "polybasis": "CHEBYSHEV", + "samplerPivot": QS([4.75, 5.25], "UNIFORM"), + "trainSetGenerator": QS([4.75, 5.25], "CHEBYSHEV"), + "MMarginal": 4, "SMarginal": 5, "polybasisMarginal": "MONOMIAL", + "samplerMarginal": QS([6.75, 7.25], "UNIFORM"), + "matchingWeight": 1., + "cutOffTolerance": 1., "cutOffType": "MAGNITUDE"} + approx = RIGP(solver, mu0, [0], params, verbosity = 100) + approx.setupApprox() + + uhP1 = approx.getApprox(mu)[0] + errP = approx.getErr(mu)[0] + errNP = approx.normErr(mu)[0] + myerrP = uhP1 - uh + assert np.allclose(np.abs(errP - myerrP), 0., rtol = 1e-3) + assert np.isclose(solver.norm(errP), errNP, rtol = 1e-3) + resP = approx.getRes(mu)[0] + resNP = approx.normRes(mu) + assert np.isclose(solver.norm(resP), resNP, rtol = 1e-3) + assert np.allclose(np.abs(resP - (solver.b(mu) - solver.A(mu).dot(uhP1))), + 0., rtol = 1e-3) + assert np.isclose(errNP / solver.norm(uh), 1.181958e-02, rtol = 1) + + out, err = capsys.readouterr() + assert ("Overriding approx_state to True." in out) + assert len(err) == 0 diff --git a/tests/reduction_methods_multiD/rational_interpolant_2d.py b/tests/reduction_methods_multiD/rational_interpolant_2d.py index 9bd84a3..067880c 100644 --- a/tests/reduction_methods_multiD/rational_interpolant_2d.py +++ b/tests/reduction_methods_multiD/rational_interpolant_2d.py @@ -1,79 +1,77 @@ # 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_random import matrixRandom from rrompy.reduction_methods.standard import RationalInterpolant as RI from rrompy.parameter.parameter_sampling import (QuadratureSampler as QS, ManualSampler as MS) def test_monomials(capsys): mu = [5.05, 7.1] mu0 = [5., 7.] solver = matrixRandom() uh = solver.solve(mu)[0] params = {"POD": False, "M": 4, "N": 4, "S": 16, "robustTol": 1e-6, "interpRcond": 1e-3, "polybasis": "MONOMIAL", "sampler": QS([[4.9, 6.85], [5.1, 7.15]], "UNIFORM")} approx = RI(solver, mu0, params, verbosity = 0) approx.setupApprox() uhP1 = approx.getApprox(mu)[0] errP = approx.getErr(mu)[0] errNP = approx.normErr(mu)[0] myerrP = uhP1 - uh assert np.allclose(np.abs(errP - myerrP), 0., rtol = 1e-3) assert np.isclose(solver.norm(errP), errNP, rtol = 1e-3) resP = approx.getRes(mu)[0] resNP = approx.normRes(mu) assert np.isclose(solver.norm(resP), resNP, rtol = 1e-3) assert np.allclose(np.abs(resP - (solver.b(mu) - solver.A(mu).dot(uhP1))), 0., rtol = 1e-3) assert np.isclose(errNP / solver.norm(uh), 5.2667e-05, rtol = 1) out, err = capsys.readouterr() - print(out) assert ("poorly conditioned. Reducing N " in out) assert len(err) == 0 def test_well_cond(): mu = [5.05, 7.1] mu0 = [5., 7.] solver = matrixRandom() params = {"POD": True, "M": 3, "N": 3, "S": 16, "interpRcond": 1e-10, "polybasis": "CHEBYSHEV", "sampler": QS([[4.9, 6.85], [5.1, 7.15]], "UNIFORM")} approx = RI(solver, mu0, params, verbosity = 0) approx.setupApprox() - print(approx.normErr(mu)[0] / approx.normHF(mu)[0]) assert np.isclose(approx.normErr(mu)[0] / approx.normHF(mu)[0], 5.98695e-05, rtol = 1e-1) def test_hermite(): mu = [5.05, 7.1] mu0 = [5., 7.] solver = matrixRandom() sampler0 = QS([[4.9, 6.85], [5.1, 7.15]], "UNIFORM") params = {"POD": True, "M": 3, "N": 3, "S": 25, "polybasis": "CHEBYSHEV", "sampler": MS([[4.9, 6.85], [5.1, 7.15]], points = sampler0.generatePoints(9))} approx = RI(solver, mu0, params, verbosity = 0) approx.setupApprox() assert np.isclose(approx.normErr(mu)[0] / approx.normHF(mu)[0], 0.000224, rtol = 5e-1)