diff --git a/rrompy/reduction_methods/pivoted/greedy/__init__.py b/rrompy/reduction_methods/pivoted/greedy/__init__.py index 5b55bed..5d473ca 100644 --- a/rrompy/reduction_methods/pivoted/greedy/__init__.py +++ b/rrompy/reduction_methods/pivoted/greedy/__init__.py @@ -1,27 +1,29 @@ # 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 .generic_pivoted_greedy_approximant import GenericPivotedGreedyApproximant from .rational_interpolant_pivoted_greedy import RationalInterpolantPivotedGreedy +from .rational_interpolant_greedy_pivoted_greedy import RationalInterpolantGreedyPivotedGreedy __all__ = [ 'GenericPivotedGreedyApproximant', - 'RationalInterpolantPivotedGreedy' + 'RationalInterpolantPivotedGreedy', + 'RationalInterpolantGreedyPivotedGreedy' ] diff --git a/rrompy/reduction_methods/pivoted/greedy/rational_interpolant_greedy_pivoted_greedy.py b/rrompy/reduction_methods/pivoted/greedy/rational_interpolant_greedy_pivoted_greedy.py new file mode 100644 index 0000000..ff8d188 --- /dev/null +++ b/rrompy/reduction_methods/pivoted/greedy/rational_interpolant_greedy_pivoted_greedy.py @@ -0,0 +1,373 @@ +#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_greedy_approximant import GenericPivotedGreedyApproximant +from rrompy.utilities.numerical import dot +from rrompy.reduction_methods.standard.greedy import RationalInterpolantGreedy +from rrompy.reduction_methods.standard.greedy.generic_greedy_approximant \ + import pruneSamples +from rrompy.reduction_methods.pivoted import (RationalInterpolantGreedyPivoted, + PODGlobal) +from rrompy.utilities.base.types import Np1D, Tuple, paramVal, paramList +from rrompy.utilities.base import verbosityManager as vbMng +from rrompy.utilities.exception_manager import RROMPyAssert +from rrompy.parameter import emptyParameterList + +__all__ = ['RationalInterpolantGreedyPivotedGreedy'] + +class RationalInterpolantGreedyPivotedGreedy(GenericPivotedGreedyApproximant, + RationalInterpolantGreedyPivoted): + """ + ROM greedy pivoted greedy rational interpolant computation for parametric + problems. + + Args: + HFEngine: HF problem solver. + mu0(optional): Default parameter. Defaults to 0. + directionPivot(optional): Pivot components. Defaults to [0]. + approxParameters(optional): Dictionary containing values for main + parameters of approximant. Recognized keys are: + - 'POD': whether to compute POD of snapshots; defaults to True; + - 'matchingWeight': weight for pole matching optimization; defaults + to 1; + - 'cutOffTolerance': tolerance for ignoring parasitic poles; + defaults to np.inf; + - 'matchingWeightError': weight for pole matching optimization in + error estimation; defaults to 0; + - 'cutOffToleranceError': tolerance for ignoring parasitic poles + in error estimation; defaults to np.inf; + - 'S': total number of pivot samples current approximant relies + upon; + - 'samplerPivot': pivot sample point generator; + - 'SMarginal': number of starting marginal samples; + - 'samplerMarginalGrid': marginal sample point generator via sparse + grid; + - '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; + - 'greedyTolMarginal': uniform error tolerance for marginal greedy + algorithm; defaults to 1e-1; + - 'maxIterMarginal': maximum number of marginal greedy steps; + defaults to 1e2; + - 'polydegreetypeMarginal': type of polynomial degree for marginal; + defaults to 'TOTAL'; + - 'radialDirectionalWeightsMarginal': radial basis weights for + marginal interpolant; defaults to 1; + - '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. + 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; + - 'matchingWeightError': weight for pole matching optimization in + error estimation; + - 'cutOffToleranceError': tolerance for ignoring parasitic poles + in error estimation; + - '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; + - 'greedyTolMarginal': uniform error tolerance for marginal greedy + algorithm; + - 'maxIterMarginal': maximum number of marginal greedy steps; + - '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. + 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; + - 'samplerMarginalGrid': marginal sample point generator via sparse + grid. + 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. + matchingWeightError: Weight for pole matching optimization in error + estimation. + cutOffToleranceError: Tolerance for ignoring parasitic poles in error + estimation. + S: Total number of pivot samples current approximant relies upon. + samplerPivot: Pivot sample point generator. + SMarginal: Total number of marginal samples current approximant relies + upon. + samplerMarginalGrid: Marginal sample point generator via sparse grid. + 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. + greedyTolMarginal: Uniform error tolerance for marginal greedy + algorithm. + maxIterMarginal: Maximum number of marginal greedy steps. + 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. + """ + + @property + def sampleBatchSize(self): + """Value of sampleBatchSize.""" + return 1 + + @property + def sampleBatchIdx(self): + """Value of sampleBatchIdx.""" + return self.S + + def _finalizeSnapshots(self): + self.samplingEngine = self._samplingEngineOld + for muM, sEN in zip(self.musMargLoc, self.samplingEngs): + self.samplingEngine.samples += [sEN.samples] + self.samplingEngine.nsamples += [sEN.nsamples] + self.samplingEngine.mus += [sEN.mus] + self.samplingEngine.musMarginal.append(muM) + self.samplingEngine._derIdxs += [[(0,) * self.npar] + for _ in range(sEN.nsamples)] + if self.POD: + self.samplingEngine.RPOD += [sEN.RPOD] + self.samplingEngine.samples_full += [copy(sEN.samples_full)] + if self.POD == PODGlobal: + self.samplingEngine.coalesceSamples(self.interpRcondMarginal) + else: + self.samplingEngine.coalesceSamples() + + def greedyNextSample(self, muidx:int, plotEst : str = "NONE")\ + -> Tuple[Np1D, int, float, paramVal]: + """Compute next greedy snapshot of solution map.""" + RROMPyAssert(self._mode, message = "Cannot add greedy sample.") + mus = copy(self.muTest[muidx]) + self.muTest.pop(muidx) + for j, mu in enumerate(mus): + vbMng(self, "MAIN", + ("Adding sample point no. {} at {} to training " + "set.").format(len(self.mus) + 1, mu), 2) + self.mus.append(mu) + self._S = len(self.mus) + self._approxParameters["S"] = self.S + if (self.samplingEngine.nsamples <= len(mus) - j - 1 + or not np.allclose(mu, + self.samplingEngine.mus.data[j - len(mus)])): + self.samplingEngine.nextSample(mu) + if self._isLastSampleCollinear(): + vbMng(self, "MAIN", + ("Collinearity above tolerance detected. Starting " + "preemptive greedy loop termination."), 2) + self._collinearityFlag = 1 + errorEstTest = np.empty(len(self.muTest)) + errorEstTest[:] = np.nan + return errorEstTest, [-1], np.nan, np.nan + errorEstTest, muidx, maxErrorEst = self.errorEstimator(self.muTest, + True) + if plotEst == "ALL": + self.plotEstimator(errorEstTest, muidx, maxErrorEst) + return errorEstTest, muidx, maxErrorEst, self.muTest[muidx] + + def _preliminaryTraining(self): + """Initialize starting snapshots of solution map.""" + RROMPyAssert(self._mode, message = "Cannot start greedy algorithm.") + if self.samplingEngine.nsamples > 0: return + self.resetSamples() + musPivot = self.trainSetGenerator.generatePoints(self.S) + musPivot.data = musPivot.data[: self.S] + muTestBasePivot = self.samplerPivot.generatePoints(self.nTestPoints) + idxPop = pruneSamples( + muTestBasePivot ** self.HFEngine.rescalingExp[self.directionPivot[0]], + musPivot ** self.HFEngine.rescalingExp[self.directionPivot[0]], + 1e-10 * self.scaleFactor[0]) + muTestBasePivot.pop(idxPop) + muTestBasePivot = muTestBasePivot.sort() + self.mus = emptyParameterList() + self.mus.reset((self.S - 1, self.HFEngine.npar)) + self.muTest = emptyParameterList() + self.muTest.reset((len(muTestBasePivot) + 1, self.HFEngine.npar)) + for k in range(self.S - 1): + self.mus.data[k, self.directionPivot] = musPivot[k].data + self.mus.data[k, self.directionMarginal] = self.musMargLoc[-1].data + for k in range(len(muTestBasePivot)): + self.muTest.data[k, self.directionPivot] = muTestBasePivot[k].data + self.muTest.data[k, self.directionMarginal] = ( + self.musMargLoc[-1].data) + self.muTest.data[-1, self.directionPivot] = musPivot[-1].data + self.muTest.data[-1, self.directionMarginal] = self.musMargLoc[-1].data + 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 + + def setupApproxPivoted(self, mus:paramList) -> int: + if self.checkComputedApproxPivoted(): return -1 + if not hasattr(self, "_plotEstPivot"): self._plotEstPivot = "NONE" + RROMPyAssert(self._mode, message = "Cannot setup approximant.") + vbMng(self, "INIT", "Setting up {}.". format(self.name()), 10) + self.computeScaleFactor() + if self.trainedModel is None: + self.trainedModel = self.tModelType() + self.trainedModel.verbosity = self.verbosity + self.trainedModel.timestamp = self.timestamp + datadict = {"mu0": self.mu0, "projMat": np.zeros((0, 0)), + "scaleFactor": self.scaleFactor, + "rescalingExp": self.HFEngine.rescalingExp, + "directionPivot": self.directionPivot} + self.trainedModel.data = self.initializeModelData(datadict)[0] + self.trainedModel.data.Qs, self.trainedModel.data.Ps = [], [] + _trainedModelOld = copy(self.trainedModel) + self._scaleFactorOldPivot = copy(self.scaleFactor) + self.scaleFactor = self.scaleFactorPivot + self._temporaryPivot = 1 + self._samplingEngineOld = copy(self.samplingEngine) + self.musMargLoc, self.samplingEngs = [], [None] * len(mus) + Qs, Ps = [None] * len(mus), [None] * len(mus) + self.verbosity -= 15 + S0 = copy(self.S) + for j, mu in enumerate(mus): + RationalInterpolantGreedy.setupSampling(self) + self.trainedModel = None + self.musMargLoc += [mu] + RationalInterpolantGreedy.setupApprox(self, self._plotEstPivot) + self.samplingEngs[j] = copy(self.samplingEngine) + Qs[j] = copy(self.trainedModel.data.Q) + Ps[j] = copy(self.trainedModel.data.P) + self._S = S0 + self.scaleFactor = self._scaleFactorOldPivot + del self._scaleFactorOldPivot, self._temporaryPivot + self._finalizeSnapshots() + del self._samplingEngineOld, self.musMargLoc, self.samplingEngs + self._mus = self.samplingEngine.musCoalesced + self.trainedModel = _trainedModelOld + self.trainedModel.data.mus = copy(self.mus) + self.trainedModel.data.musMarginal = copy(self.musMarginal) + padRight = (self.samplingEngine.nsamplesTot + - self.trainedModel.data.projMat.shape[1]) + nmusOld = len(self.trainedModel.data.Ps) + for j in range(nmusOld): + nsj = self.samplingEngine.nsamples[j] + self.trainedModel.data.Ps[j].pad(0, padRight) + self.trainedModel.data.HIs[j].pad(0, padRight) + padLeft = self.trainedModel.data.projMat.shape[1] + for j in range(len(mus)): + nsj = self.samplingEngine.nsamples[nmusOld + j] + if self.POD == PODGlobal: + rRightj = self.samplingEngine.RPODCPart[:, + padLeft : padLeft + nsj] + Ps[j].postmultiplyTensorize(rRightj.T) + else: + padRight -= nsj + Ps[j].pad(padLeft, padRight) + padLeft += nsj + pMat = self.samplingEngine.samplesCoalesced.data + pMatEff = dot(self.HFEngine.C, pMat) if self.approx_state else pMat + self.trainedModel.data.projMat = pMatEff + self.trainedModel.data.Qs += Qs + self.trainedModel.data.Ps += Ps + self.trainedModel.data.approxParameters = copy(self.approxParameters) + self.verbosity += 15 + vbMng(self, "DEL", "Done setting up approximant.", 10) + return 0 + + def setupApprox(self, plotEst : str = "NONE") -> int: + if '_' not in plotEst: plotEst = plotEst + "_NONE" + plotEstM, self._plotEstPivot = plotEst.split("_") + val = super().setupApprox(plotEstM) + return val + \ No newline at end of file diff --git a/tests/reduction_methods_multiD/greedy_pivoted_rational_2d.py b/tests/reduction_methods_multiD/greedy_pivoted_rational_2d.py new file mode 100644 index 0000000..dc9509e --- /dev/null +++ b/tests/reduction_methods_multiD/greedy_pivoted_rational_2d.py @@ -0,0 +1,83 @@ +# 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.greedy import ( + RationalInterpolantPivotedGreedy as RIPG, + RationalInterpolantGreedyPivotedGreedy as RIGPG) +from rrompy.parameter.parameter_sampling import QuadratureSampler as QS +from rrompy.parameter import localSparseGrid as LSG + +def test_greedy_pivoted(): + 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": 1, "SMarginal": 3, "greedyTolMarginal": 1e-2, + "radialDirectionalWeightsMarginal": 2., + "polybasisMarginal": "MONOMIAL_GAUSSIAN", + "matchingWeight": 1., "samplerMarginalGrid":LSG([6.75, 7.25])} + approx = RIPG([0], solver, mu0, approx_state = True, + approxParameters = 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_greedy_pivoted_greedy(): + mu = [5.05, 7.1] + mu0 = [5., 7.] + solver = matrixRandom() + uh = solver.solve(mu)[0] + params = {"POD": True, "nTestPoints": 100, "greedyTol": 1e-3, "S": 2, + "polybasis": "CHEBYSHEV", + "samplerPivot": QS([4.75, 5.25], "CHEBYSHEV"), + "trainSetGenerator": QS([4.75, 5.25], "CHEBYSHEV"), + "MMarginal": 1, "SMarginal": 3, "greedyTolMarginal": 1e-2, + "radialDirectionalWeightsMarginal": 2., + "polybasisMarginal": "MONOMIAL_GAUSSIAN", + "matchingWeight": 1., "samplerMarginalGrid":LSG([6.75, 7.25])} + approx = RIGPG([0], solver, mu0, approx_state = True, + approxParameters = 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)