diff --git a/examples/1_symmetric_disk/symmetric_disk.py b/examples/1_symmetric_disk/symmetric_disk.py index 898e2b1..7cd9be6 100644 --- a/examples/1_symmetric_disk/symmetric_disk.py +++ b/examples/1_symmetric_disk/symmetric_disk.py @@ -1,84 +1,86 @@ import numpy as np from symmetric_disk_engine import SymmetricDiskEngine as engine from rrompy.sampling import SamplingEngineStandard as SES from rrompy.reduction_methods import ( NearestNeighbor as NN, RationalInterpolant as RI, ReducedBasis as RB, RationalInterpolantGreedy as RIG, ReducedBasisGreedy as RBG) from rrompy.parameter.parameter_sampling import QuadratureSampler as QS ks = [10., 20.] k0, n = np.mean(np.power(ks, 2.)) ** .5, 150 solver = engine(k0, n) k = 12. for method in ["RI", "RB", "RI_GREEDY", "RB_GREEDY"]: print("Testing {} method".format(method)) if method == "RI": params = {'S':40, 'POD':True, 'polybasis':"CHEBYSHEV", 'sampler':QS(ks, "CHEBYSHEV", scalingExp = 2.)} algo = RI if method == "RB": params = {'S':40, 'POD':True, 'sampler':QS(ks, "CHEBYSHEV", scalingExp = 2.)} algo = RB if method == "RI_GREEDY": params = {'S':10, 'POD':True, 'polybasis':"LEGENDRE", 'greedyTol':1e-2, 'sampler':QS(ks, "UNIFORM", scalingExp = 2.), - 'errorEstimatorKind':"DISCREPANCY"} + 'errorEstimatorKind':"DISCREPANCY", + 'trainSetGenerator':QS(ks, "CHEBYSHEV", scalingExp = 2.)} algo = RIG if method == "RB_GREEDY": params = {'S':10, 'POD':True, 'greedyTol':1e-2, - 'sampler':QS(ks, "UNIFORM", scalingExp = 2.)} + 'sampler':QS(ks, "UNIFORM", scalingExp = 2.), + 'trainSetGenerator':QS(ks, "CHEBYSHEV", scalingExp = 2.)} algo = RBG approx = algo(solver, mu0 = k0, approx_state = True, approxParameters = params, verbosity = 20) if len(method) == 2: approx.setupApprox() else: approx.setupApprox("LAST") print("--- Approximant ---") approx.plotApprox(k, name = 'u_app') approx.plotHF(k, name = 'u_HF') approx.plotErr(k, name = 'err_app') approx.plotRes(k, name = 'res_app') normErr = approx.normErr(k)[0] normSol = approx.normHF(k)[0] normRes = approx.normRes(k)[0] normRHS = approx.normRHS(k)[0] print("SolNorm:\t{:.5e}\nErr_app: \t{:.5e}\nErrRel_app:\t{:.5e}".format( normSol, normErr, normErr / normSol)) print("RHSNorm:\t{:.5e}\nRes_app: \t{:.5e}\nResRel_app:\t{:.5e}".format( normRHS, normRes, normRes / normRHS)) print("--- Closest snapshot ---") eng = SES(solver, verbosity = 0) approxNN = NN(solver, mu0 = k0, approx_state = True, verbosity = 0, approxParameters = {'S':approx.S, 'POD':True}) approxNN.setSamples(approx.samplingEngine) approxNN.plotApprox(k, name = 'u_close') approxNN.plotHF(k, name = 'u_HF') approxNN.plotErr(k, name = 'err_close') approxNN.plotRes(k, name = 'res_close') normErr = approxNN.normErr(k)[0] normSol = approxNN.normHF(k)[0] normRes = approxNN.normRes(k)[0] normRHS = approxNN.normRHS(k)[0] print("SolNorm:\t{:.5e}\nErr_close:\t{:.5e}\nErrRel_close:\t{:.5e}".format( normSol, normErr, normErr / normSol)) print("RHSNorm:\t{:.5e}\nRes_close:\t{:.5e}\nResRel_close:\t{:.5e}".format( normRHS, normRes, normRes / normRHS)) if method[:2] == "RI": poles, residues = approx.getResidues() if method[:2] == "RB": poles = approx.getPoles() print("Poles:\n{}".format(poles)) if method[:2] == "RI": for pol, res in zip(poles, residues.T): solver.plot(res) print("pole = {:.5e}".format(pol)) print("\n") diff --git a/examples/7_MHD/mhd.py b/examples/7_MHD/mhd.py index c4d3232..8e58b3d 100644 --- a/examples/7_MHD/mhd.py +++ b/examples/7_MHD/mhd.py @@ -1,65 +1,73 @@ import numpy as np import matplotlib.pyplot as plt from mhd_engine import MHDEngine as engine from rrompy.reduction_methods import (RationalInterpolant as RI, RationalInterpolantGreedy as RIG) from rrompy.parameter.parameter_sampling import (FFTSampler as FFTS, QuadratureCircleSampler as QCS, QuadratureBoxSampler as QBS) ks = [-.35 + .5j, 0. + .5j] k0 = np.mean(ks) solver = engine(5) kEffDelta = .1 * (ks[1] - ks[0]) kEff = np.real([ks[0] - kEffDelta, ks[1] + kEffDelta]) iEff = kEff - .5 * np.sum(np.real(ks)) + np.imag(ks[0]) nPoles = 50 polesEx = solver.getPolesExact(nPoles, k0) -for method in ["FFT", "BOX", "GREEDY"]: - print("Testing {} method".format(method)) - - if method == "FFT": - params = {'S':64, 'POD':True, 'polybasis':"MONOMIAL", - 'sampler':FFTS(ks)} - algo = RI - if method == "BOX": - params = {'S':64, 'POD':True, 'polybasis':"MONOMIAL", - 'sampler':QBS(ks)} - algo = RI - if method == "GREEDY": - params = {'S':30, 'POD':True, 'greedyTol':1e-2, 'polybasis':"MONOMIAL", - 'sampler':QCS(ks), 'errorEstimatorKind':"LOOK_AHEAD", - 'nTestPoints':10000, 'trainSetGenerator':FFTS(ks)} - algo = RIG - - approx = algo(solver, mu0 = k0, approx_state = True, - approxParameters = params, verbosity = 10) - approx.setupApprox() +for corrector in [False, True]: + for method in ["FFT", "BOX", "GREEDY"]: + print("Testing {} method with{} corrector".format(method, + "out" * (not corrector))) - poles, residues = approx.getResidues() - inRange = np.logical_and( + if method == "FFT": + params = {'S':64, 'POD':True, 'polybasis':"MONOMIAL", + 'sampler':FFTS(ks), 'correctorTol':1e-5} + algo = RI + if method == "BOX": + params = {'S':64, 'POD':True, 'polybasis':"MONOMIAL", + 'sampler':QBS(ks), 'correctorTol':1e-5} + algo = RI + if method == "GREEDY": + params = {'S':30, 'POD':True, 'greedyTol':1e-2, + 'polybasis':"MONOMIAL", 'sampler':QCS(ks), + 'errorEstimatorKind':"LOOK_AHEAD", 'nTestPoints':10000, + 'trainSetGenerator':FFTS(ks), 'correctorTol':1e-5} + algo = RIG + params['correctorForce'] = corrector + + approx = algo(solver, mu0 = k0, approx_state = True, + approxParameters = params, verbosity = 10) + approx.setupApprox() + + poles, residues = approx.getResidues() + inRange = np.logical_and( np.logical_and(np.real(poles) >= kEff[0], np.real(poles) <= kEff[1]), np.logical_and(np.imag(poles) >= iEff[0], np.imag(poles) <= iEff[1])) - polesEff = poles[inRange] - resNormEff = np.linalg.norm(residues, axis = 0)[inRange] - rLm = np.min(np.log(resNormEff)) - rLmM = np.max(np.log(resNormEff)) - rLm - fig = plt.figure(figsize = (10, 10)) - ax = fig.add_subplot(1, 1, 1) - if method == "GREEDY": - ax.plot(approx.muTest.re.data.flatten(), - approx.muTest.im.data.flatten(), 'k,', alpha = 0.25) - for pl, rN in zip(polesEff, resNormEff): - alpha = .1 + .65 * (np.log(rN) - rLm) / rLmM - ax.annotate("{0:.0e}".format(rN), (np.real(pl), np.imag(pl)), - alpha = alpha) - ax.plot(np.real(pl), np.imag(pl), 'r+', alpha = alpha + .25) - ax.plot(approx.mus.re.data.flatten(), approx.mus.im.data.flatten(), 'k.') - ax.plot(np.real(polesEx), np.imag(polesEx), 'bx') - ax.set_xlim(kEff) - ax.set_ylim(iEff) - ax.grid() - plt.tight_layout() - plt.show() - - print("\n") + polesEff = poles[inRange] + resNormEff = np.linalg.norm(residues, axis = 0)[inRange] + rLm = np.min(np.log(resNormEff)) + rLmM = np.max(np.log(resNormEff)) - rLm + fig = plt.figure(figsize = (10, 10)) + ax = fig.add_subplot(1, 1, 1) + if method == "GREEDY": + ax.plot(approx.muTest.re.data.flatten(), + approx.muTest.im.data.flatten(), 'k,', alpha = 0.25) + for pl, rN in zip(polesEff, resNormEff): + if corrector: + alpha = .35 + .4 * (np.log(rN) - rLm) / rLmM + else: + alpha = .1 + .65 * (np.log(rN) - rLm) / rLmM + ax.annotate("{0:.0e}".format(rN), (np.real(pl), np.imag(pl)), + alpha = alpha) + ax.plot(np.real(pl), np.imag(pl), 'r+', alpha = alpha + .25) + ax.plot(approx.mus.re.data.flatten(), + approx.mus.im.data.flatten(), 'k.') + ax.plot(np.real(polesEx), np.imag(polesEx), 'bx') + ax.set_xlim(kEff) + ax.set_ylim(iEff) + ax.grid() + plt.tight_layout() + plt.show() + + print("\n") diff --git a/rrompy/reduction_methods/pivoted/generic_pivoted_approximant.py b/rrompy/reduction_methods/pivoted/generic_pivoted_approximant.py index e43b294..b086c2f 100644 --- a/rrompy/reduction_methods/pivoted/generic_pivoted_approximant.py +++ b/rrompy/reduction_methods/pivoted/generic_pivoted_approximant.py @@ -1,486 +1,486 @@ # Copyright (C) 2018 by the RROMPy authors # # This file is part of RROMPy. # # RROMPy is free software: you can redistribute it and/or modify # it under the terms of the GNU Lesser General Public License as published by # the Free Software Foundation, either version 3 of the License, or # (at your option) any later version. # # RROMPy is distributed in the hope that it will be useful, # but WITHOUT ANY WARRANTY; without even the implied warranty of # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the # GNU Lesser General Public License for more details. # # You should have received a copy of the GNU Lesser General Public License # along with RROMPy. If not, see . # import numpy as np from copy import deepcopy as copy from rrompy.reduction_methods.base.generic_approximant import ( GenericApproximant) from rrompy.utilities.poly_fitting.polynomial import polybases as ppb from rrompy.utilities.poly_fitting.radial_basis import polybases as rbpb from rrompy.utilities.poly_fitting.moving_least_squares import ( polybases as mlspb) from rrompy.sampling import (SamplingEnginePivoted, SamplingEnginePivotedPOD, SamplingEnginePivotedPODGlobal) from rrompy.utilities.base.types import paramList, ListAny from rrompy.utilities.base import verbosityManager as vbMng from rrompy.utilities.numerical.degree import reduceDegreeN from rrompy.utilities.exception_manager import (RROMPyException, RROMPyAssert, RROMPyWarning) __all__ = ['GenericPivotedApproximant', 'PODGlobal'] PODGlobal = 2 class GenericPivotedApproximant(GenericApproximant): """ ROM pivoted approximant (with pole matching) computation for parametric problems (ABSTRACT). Args: HFEngine: HF problem solver. mu0(optional): Default parameter. Defaults to 0. directionPivot(optional): Pivot components. Defaults to [0]. approxParameters(optional): Dictionary containing values for main parameters of approximant. Recognized keys are: - 'POD': whether to compute POD of snapshots; defaults to True; - 'matchingWeight': weight for pole matching optimization; defaults to 1; - 'cutOffTolerance': tolerance for ignoring parasitic poles; defaults to np.inf; - 'cutOffKind': kind of cut off strategy; available values include 'SOFT' and 'HARD'; defaults to 'HARD'; - 'S': total number of pivot samples current approximant relies upon; - 'samplerPivot': pivot sample point generator; - 'SMarginal': total number of marginal samples current approximant relies upon; - 'samplerMarginal': marginal sample point generator; - 'polybasisMarginal': type of polynomial basis for marginal interpolation; allowed values include 'MONOMIAL', 'CHEBYSHEV' and 'LEGENDRE'; defaults to 'MONOMIAL'; - 'MMarginal': degree of marginal interpolant; defaults to 'AUTO', i.e. maximum allowed; - '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; - 'interpRcondMarginal': tolerance for marginal interpolation; defaults to None. 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; - 'cutOffKind': kind of cut off strategy; - 'polybasisMarginal': type of polynomial basis for marginal interpolation; - '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; - 'interpRcondMarginal': tolerance for marginal interpolation. parameterListCritical: Recognized keys of critical approximant parameters: - 'S': total number of pivot samples current approximant relies upon; - 'samplerPivot': pivot sample point generator; - 'SMarginal': total number of marginal samples current approximant relies upon; - 'samplerMarginal': marginal sample point generator. 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. cutOffKind: Kind of cut off strategy. S: Total number of pivot samples current approximant relies upon. samplerPivot: Pivot sample point generator. SMarginal: Total number of marginal samples current approximant relies upon. samplerMarginal: Marginal sample point generator. polybasisMarginal: Type of polynomial basis for marginal interpolation. 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. interpRcondMarginal: Tolerance for marginal interpolation. - muBoundsPivot: list of bounds for pivot parameter values. + muBounds: list of bounds for pivot parameter values. muBoundsMarginal: list of bounds for marginal parameter values. samplingEngine: Sampling engine. uHF: High fidelity solution(s) with parameter(s) lastSolvedHF as sampleList. lastSolvedHF: Parameter(s) corresponding to last computed high fidelity solution(s) as parameterList. uApproxReduced: Reduced approximate solution(s) with parameter(s) lastSolvedApprox as sampleList. lastSolvedApproxReduced: Parameter(s) corresponding to last computed reduced approximate solution(s) as parameterList. uApprox: Approximate solution(s) with parameter(s) lastSolvedApprox as sampleList. lastSolvedApprox: Parameter(s) corresponding to last computed approximate solution(s) as parameterList. """ def __init__(self, directionPivot:ListAny, *args, **kwargs): self._preInit() if len(directionPivot) > 1: raise RROMPyException(("Exactly 1 pivot parameter allowed in pole " "matching.")) from rrompy.parameter.parameter_sampling import QuadratureSampler as QS QSBase = QS([[0], [1]], "UNIFORM") self._addParametersToList(["matchingWeight", "cutOffTolerance", "cutOffKind", "polybasisMarginal", "MMarginal", "polydegreetypeMarginal", "radialDirectionalWeightsMarginal", "nNearestNeighborMarginal", "interpRcondMarginal"], [1, np.inf, "HARD", "MONOMIAL", "AUTO", "TOTAL", [1], -1, -1], ["samplerPivot", "SMarginal", "samplerMarginal"], [QSBase, [1], QSBase]) del QS self._directionPivot = directionPivot super().__init__(*args, **kwargs) self._postInit() @property def tModelType(self): from .trained_model.trained_model_pivoted import TrainedModelPivoted return TrainedModelPivoted def setupSampling(self): """Setup sampling engine.""" RROMPyAssert(self._mode, message = "Cannot setup sampling engine.") if not hasattr(self, "_POD") or self._POD is None: return if self.POD: if self.POD == PODGlobal: SamplingEngine = SamplingEnginePivotedPODGlobal else: SamplingEngine = SamplingEnginePivotedPOD else: SamplingEngine = SamplingEnginePivoted self.samplingEngine = SamplingEngine(self.HFEngine, self.directionPivot, sample_state = self.approx_state, verbosity = self.verbosity) def initializeModelData(self, datadict): if "directionPivot" in datadict.keys(): from .trained_model.trained_model_pivoted_data import ( TrainedModelPivotedData) return (TrainedModelPivotedData(datadict["mu0"], datadict.pop("projMat"), datadict["scaleFactor"], datadict.pop("rescalingExp"), datadict["directionPivot"]), ["mu0", "scaleFactor", "directionPivot", "mus"]) else: return super().initializeModelData(datadict) @property def npar(self): """Number of parameters.""" if hasattr(self, "_temporaryPivot"): return self.nparPivot return super().npar @property def mus(self): """Value of mus. Its assignment may reset snapshots.""" return self._mus @mus.setter def mus(self, mus): musOld = copy(self.mus) if hasattr(self, '_mus') else None if (musOld is None or len(mus) != len(musOld) or not mus == musOld): self.resetSamples() self._mus = mus @property def matchingWeight(self): """Value of matchingWeight.""" return self._matchingWeight @matchingWeight.setter def matchingWeight(self, matchingWeight): self._matchingWeight = matchingWeight self._approxParameters["matchingWeight"] = self.matchingWeight @property def cutOffTolerance(self): """Value of cutOffTolerance.""" return self._cutOffTolerance @cutOffTolerance.setter def cutOffTolerance(self, cutOffTolerance): self._cutOffTolerance = cutOffTolerance self._approxParameters["cutOffTolerance"] = self.cutOffTolerance @property def cutOffKind(self): """Value of cutOffKind.""" return self._cutOffKind @cutOffKind.setter def cutOffKind(self, cutOffKind): cutOffKind = cutOffKind.upper() if cutOffKind not in ["SOFT", "HARD"]: RROMPyWarning(("Cut off kind not recognized. Overriding to " "'HARD'.")) cutOffKind = "HARD" self._cutOffKind = cutOffKind self._approxParameters["cutOffKind"] = self.cutOffKind @property def SMarginal(self): """Value of SMarginal.""" return self._SMarginal @SMarginal.setter def SMarginal(self, SMarginal): if SMarginal <= 0: raise RROMPyException("SMarginal must be positive.") if hasattr(self, "_SMarginal") and self._SMarginal is not None: Sold = self.SMarginal else: Sold = -1 self._SMarginal = SMarginal self._approxParameters["SMarginal"] = self.SMarginal if Sold != self.SMarginal: self.resetSamples() @property def polybasisMarginal(self): """Value of polybasisMarginal.""" return self._polybasisMarginal @polybasisMarginal.setter def polybasisMarginal(self, polybasisMarginal): try: polybasisMarginal = polybasisMarginal.upper().strip().replace(" ", "") if polybasisMarginal not in ppb + rbpb + mlspb: raise RROMPyException( "Prescribed marginal polybasis not recognized.") self._polybasisMarginal = polybasisMarginal except: RROMPyWarning(("Prescribed marginal polybasis not recognized. " "Overriding to 'MONOMIAL'.")) self._polybasisMarginal = "MONOMIAL" self._approxParameters["polybasisMarginal"] = self.polybasisMarginal @property def MMarginal(self): """Value of MMarginal.""" return self._MMarginal @MMarginal.setter def MMarginal(self, MMarginal): if isinstance(MMarginal, str): MMarginal = MMarginal.strip().replace(" ","") if "-" not in MMarginal: MMarginal = MMarginal + "-0" self._MMarginal_isauto = True self._MMarginal_shift = int(MMarginal.split("-")[-1]) MMarginal = 0 if MMarginal < 0: raise RROMPyException("MMarginal must be non-negative.") self._MMarginal = MMarginal self._approxParameters["MMarginal"] = self.MMarginal def _setMMarginalAuto(self): self.MMarginal = max(0, reduceDegreeN( len(self.musMarginal), len(self.musMarginal), self.nparMarginal, self.polydegreetypeMarginal ) - self._MMarginal_shift) vbMng(self, "MAIN", ("Automatically setting MMarginal to " "{}.").format(self.MMarginal), 25) @property def polydegreetypeMarginal(self): """Value of polydegreetypeMarginal.""" return self._polydegreetypeMarginal @polydegreetypeMarginal.setter def polydegreetypeMarginal(self, polydegreetypeM): try: polydegreetypeM = polydegreetypeM.upper().strip().replace(" ","") if polydegreetypeM not in ["TOTAL", "FULL"]: raise RROMPyException(("Prescribed polydegreetypeMarginal not " "recognized.")) self._polydegreetypeMarginal = polydegreetypeM except: RROMPyWarning(("Prescribed polydegreetypeMarginal not recognized. " "Overriding to 'TOTAL'.")) self._polydegreetypeMarginal = "TOTAL" self._approxParameters["polydegreetypeMarginal"] = ( self.polydegreetypeMarginal) @property def radialDirectionalWeightsMarginal(self): """Value of radialDirectionalWeightsMarginal.""" return self._radialDirectionalWeightsMarginal @radialDirectionalWeightsMarginal.setter def radialDirectionalWeightsMarginal(self, radialDirWeightsMarginal): if not hasattr(radialDirWeightsMarginal, "__len__"): radialDirWeightsMarginal = [radialDirWeightsMarginal] self._radialDirectionalWeightsMarginal = radialDirWeightsMarginal self._approxParameters["radialDirectionalWeightsMarginal"] = ( self.radialDirectionalWeightsMarginal) @property def nNearestNeighborMarginal(self): """Value of nNearestNeighborMarginal.""" return self._nNearestNeighborMarginal @nNearestNeighborMarginal.setter def nNearestNeighborMarginal(self, nNearestNeighborMarginal): self._nNearestNeighborMarginal = nNearestNeighborMarginal self._approxParameters["nNearestNeighborMarginal"] = ( self.nNearestNeighborMarginal) @property def interpRcondMarginal(self): """Value of interpRcondMarginal.""" return self._interpRcondMarginal @interpRcondMarginal.setter def interpRcondMarginal(self, interpRcondMarginal): self._interpRcondMarginal = interpRcondMarginal self._approxParameters["interpRcondMarginal"] = ( self.interpRcondMarginal) @property def directionPivot(self): """Value of directionPivot. Its assignment may reset snapshots.""" return self._directionPivot @directionPivot.setter def directionPivot(self, directionPivot): if hasattr(self, '_directionPivot'): directionPivotOld = copy(self.directionPivot) else: directionPivotOld = None if (directionPivotOld is None or len(directionPivot) != len(directionPivotOld) or not directionPivot == directionPivotOld): self.resetSamples() self._directionPivot = directionPivot @property def directionMarginal(self): return [x for x in range(self.HFEngine.npar) \ if x not in self.directionPivot] @property def nparPivot(self): return len(self.directionPivot) @property def nparMarginal(self): return self.npar - self.nparPivot @property def rescalingExpPivot(self): return [self.HFEngine.rescalingExp[x] for x in self.directionPivot] @property def rescalingExpMarginal(self): return [self.HFEngine.rescalingExp[x] for x in self.directionMarginal] @property - def muBoundsPivot(self): - """Value of muBoundsPivot.""" + def muBounds(self): + """Value of muBounds.""" return self.samplerPivot.lims @property def muBoundsMarginal(self): """Value of muBoundsMarginal.""" return self.samplerMarginal.lims @property def samplerPivot(self): """Value of samplerPivot.""" return self._samplerPivot @samplerPivot.setter def samplerPivot(self, samplerPivot): if 'generatePoints' not in dir(samplerPivot): raise RROMPyException("Pivot sampler type not recognized.") if hasattr(self, '_samplerPivot') and self._samplerPivot is not None: samplerOld = self.samplerPivot self._samplerPivot = samplerPivot self._approxParameters["samplerPivot"] = self.samplerPivot.__str__() if not 'samplerOld' in locals() or samplerOld != self.samplerPivot: self.resetSamples() @property def samplerMarginal(self): """Value of samplerMarginal.""" return self._samplerMarginal @samplerMarginal.setter def samplerMarginal(self, samplerMarginal): if 'generatePoints' not in dir(samplerMarginal): raise RROMPyException("Marginal sampler type not recognized.") if (hasattr(self, '_samplerMarginal') and self._samplerMarginal is not None): samplerOld = self.samplerMarginal self._samplerMarginal = samplerMarginal self._approxParameters["samplerMarginal"] = ( self.samplerMarginal.__str__()) if not 'samplerOld' in locals() or samplerOld != self.samplerMarginal: self.resetSamples() def setSamples(self, samplingEngine): """Copy samplingEngine and samples.""" self.mus = copy(samplingEngine.mus[0]) for sEj in samplingEngine.mus[1:]: self.mus.append(sEj) super().setSamples(samplingEngine) def _finalizeMarginalization(self): vbMng(self, "INIT", "Recompressing by cut off.", 10) msg = self.trainedModel.recompressByCutOff( self.cutOffTolerance, self.cutOffKind, self.samplerPivot.normalFoci(), self.samplerPivot.groundPotential()) vbMng(self, "DEL", "Done recompressing." + msg, 10) interpPars = [self.verbosity >= 5, self.polydegreetypeMarginal == "TOTAL", {}] if self.polybasisMarginal not in ppb: interpPars[-1]["nNearestNeighbor"] = self.nNearestNeighborMarginal if self.polybasisMarginal in ppb + rbpb: interpPars += [{"rcond": self.interpRcondMarginal}] self.trainedModel.setupMarginalInterp(self, interpPars, hasattr(self, "_MMarginal_isauto"), self.radialDirectionalWeightsMarginal, hasattr(self, "_reduceDegreeNNoWarn")) self.trainedModel.data.approxParameters = copy(self.approxParameters) def computeScaleFactor(self): """Compute parameter rescaling factor.""" RROMPyAssert(self._mode, message = "Cannot compute rescaling factor.") self.scaleFactorPivot = .5 * np.abs( - self.muBoundsPivot[0] ** self.rescalingExpPivot - - self.muBoundsPivot[1] ** self.rescalingExpPivot) + self.muBounds[0] ** self.rescalingExpPivot + - self.muBounds[1] ** self.rescalingExpPivot) self.scaleFactorMarginal = .5 * np.abs( self.muBoundsMarginal[0] ** self.rescalingExpMarginal - self.muBoundsMarginal[1] ** self.rescalingExpMarginal) self.scaleFactor = np.empty(self.npar) self.scaleFactor[self.directionPivot] = self.scaleFactorPivot self.scaleFactor[self.directionMarginal] = self.scaleFactorMarginal def normApprox(self, mu:paramList) -> float: _PODOld = self.POD self._POD = self.POD == PODGlobal result = super().normApprox(mu) self._POD = _PODOld return result diff --git a/rrompy/reduction_methods/pivoted/greedy/generic_pivoted_greedy_approximant.py b/rrompy/reduction_methods/pivoted/greedy/generic_pivoted_greedy_approximant.py index 6286694..5a07185 100644 --- a/rrompy/reduction_methods/pivoted/greedy/generic_pivoted_greedy_approximant.py +++ b/rrompy/reduction_methods/pivoted/greedy/generic_pivoted_greedy_approximant.py @@ -1,454 +1,454 @@ # 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 matplotlib import pyplot as plt from rrompy.reduction_methods.pivoted.generic_pivoted_approximant import ( GenericPivotedApproximant, PODGlobal) from rrompy.utilities.base.types import Np1D, Tuple, List, paramVal, paramList from rrompy.utilities.base import verbosityManager as vbMng from rrompy.utilities.numerical.point_matching import (pointMatching, chordalMetricAdjusted, potential) from rrompy.utilities.exception_manager import (RROMPyException, RROMPyAssert, RROMPyWarning) from rrompy.parameter import checkParameterList, emptyParameterList __all__ = ['GenericPivotedGreedyApproximant'] class GenericPivotedGreedyApproximant(GenericPivotedApproximant): """ ROM pivoted greedy interpolant computation for parametric problems (ABSTRACT). Args: HFEngine: HF problem solver. mu0(optional): Default parameter. Defaults to 0. directionPivot(optional): Pivot components. Defaults to [0]. approxParameters(optional): Dictionary containing values for main parameters of approximant. Recognized keys are: - 'POD': whether to compute POD of snapshots; defaults to True; - 'matchingWeight': weight for pole matching optimization; defaults to 1; - 'cutOffTolerance': tolerance for ignoring parasitic poles; defaults to np.inf; - 'cutOffKind': kind of cut off strategy; available values include 'SOFT' and 'HARD'; defaults to 'HARD'; - 'matchingWeightError': weight for pole matching optimization in error estimation; defaults to 0; - 'cutOffToleranceError': tolerance for ignoring parasitic poles in error estimation; defaults to 'AUTO', i.e. cutOffTolerance; - 'S': total number of pivot samples current approximant relies upon; - 'samplerPivot': pivot sample point generator; - 'SMarginal': number of starting marginal samples; - 'samplerMarginalGrid': marginal sample point generator via sparse grid; - 'polybasisMarginal': type of polynomial basis for marginal interpolation; allowed values include 'MONOMIAL', 'CHEBYSHEV' and 'LEGENDRE'; defaults to 'MONOMIAL'; - 'MMarginal': degree of marginal interpolant; defaults to 'AUTO', i.e. maximum allowed; - '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; - 'interpRcondMarginal': tolerance for marginal interpolation; defaults to None. 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; - 'cutOffKind': kind of cut off strategy; - 'matchingWeightError': weight for pole matching optimization in error estimation; - 'cutOffToleranceError': tolerance for ignoring parasitic poles in error estimation; - 'polybasisMarginal': type of polynomial basis for marginal interpolation; - '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; - 'interpRcondMarginal': tolerance for marginal interpolation. parameterListCritical: Recognized keys of critical approximant parameters: - 'S': total number of pivot samples current approximant relies upon; - 'samplerPivot': pivot sample point generator; - 'SMarginal': total number of marginal samples current approximant relies upon; - '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. cutOffKind: Kind of cut off strategy. matchingWeightError: Weight for pole matching optimization in error estimation. cutOffToleranceError: Tolerance for ignoring parasitic poles in error estimation. S: Total number of pivot samples current approximant relies upon. samplerPivot: Pivot sample point generator. SMarginal: Total number of marginal samples current approximant relies upon. samplerMarginalGrid: Marginal sample point generator via sparse grid. polybasisMarginal: Type of polynomial basis for marginal interpolation. 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. interpRcondMarginal: Tolerance for marginal interpolation. - muBoundsPivot: list of bounds for pivot parameter values. + muBounds: list of bounds for pivot parameter values. muBoundsMarginal: list of bounds for marginal parameter values. samplingEngine: Sampling engine. uHF: High fidelity solution(s) with parameter(s) lastSolvedHF as sampleList. lastSolvedHF: Parameter(s) corresponding to last computed high fidelity solution(s) as parameterList. uApproxReduced: Reduced approximate solution(s) with parameter(s) lastSolvedApprox as sampleList. lastSolvedApproxReduced: Parameter(s) corresponding to last computed reduced approximate solution(s) as parameterList. uApprox: Approximate solution(s) with parameter(s) lastSolvedApprox as sampleList. lastSolvedApprox: Parameter(s) corresponding to last computed approximate solution(s) as parameterList. """ def __init__(self, *args, **kwargs): self._preInit() from rrompy.parameter import localSparseGrid as SG SGBase = SG([[0.], [1.]], "UNIFORM") self._addParametersToList(["matchingWeightError", "cutOffToleranceError", "greedyTolMarginal", "maxIterMarginal"], [0., "AUTO", 1e-1, 1e2], ["samplerMarginalGrid"], [SGBase], toBeExcluded = ["samplerMarginal"]) super().__init__(*args, **kwargs) self._postInit() @property def muBoundsMarginal(self): """Value of muBoundsMarginal.""" return self.samplerMarginalGrid.lims @property def samplerMarginalGrid(self): """Value of samplerMarginalGrid.""" return self._samplerMarginalGrid @samplerMarginalGrid.setter def samplerMarginalGrid(self, samplerMarginalGrid): if 'refine' not in dir(samplerMarginalGrid): raise RROMPyException("Marginal sampler type not recognized.") if (hasattr(self, '_samplerMarginalGrid') and self._samplerMarginalGrid is not None): samplerOld = self.samplerMarginalGrid self._samplerMarginalGrid = samplerMarginalGrid self._approxParameters["samplerMarginalGrid"] = ( self.samplerMarginalGrid.__str__()) if (not 'samplerOld' in locals() or samplerOld != self.samplerMarginalGrid): self.resetSamples() @property def matchingWeightError(self): """Value of matchingWeightError.""" return self._matchingWeightError @matchingWeightError.setter def matchingWeightError(self, matchingWeightError): self._matchingWeightError = matchingWeightError self._approxParameters["matchingWeightError"] = ( self.matchingWeightError) @property def cutOffToleranceError(self): """Value of cutOffToleranceError.""" return self._cutOffToleranceError @cutOffToleranceError.setter def cutOffToleranceError(self, cutOffToleranceError): if isinstance(cutOffToleranceError, (str,)): cutOffToleranceError = cutOffToleranceError.upper()\ .strip().replace(" ","") if cutOffToleranceError != "AUTO": RROMPyWarning(("String value of cutOffToleranceError not " "recognized. Overriding to 'AUTO'.")) cutOffToleranceError == "AUTO" self._cutOffToleranceError = cutOffToleranceError self._approxParameters["cutOffToleranceError"] = ( self.cutOffToleranceError) @property def greedyTolMarginal(self): """Value of greedyTolMarginal.""" return self._greedyTolMarginal @greedyTolMarginal.setter def greedyTolMarginal(self, greedyTolMarginal): if greedyTolMarginal < 0: raise RROMPyException("greedyTolMarginal must be non-negative.") if (hasattr(self, "_greedyTolMarginal") and self.greedyTolMarginal is not None): greedyTolMarginalold = self.greedyTolMarginal else: greedyTolMarginalold = -1 self._greedyTolMarginal = greedyTolMarginal self._approxParameters["greedyTolMarginal"] = self.greedyTolMarginal if greedyTolMarginalold != self.greedyTolMarginal: self.resetSamples() @property def maxIterMarginal(self): """Value of maxIterMarginal.""" return self._maxIterMarginal @maxIterMarginal.setter def maxIterMarginal(self, maxIterMarginal): if maxIterMarginal <= 0: raise RROMPyException("maxIterMarginal must be positive.") if (hasattr(self, "_maxIterMarginal") and self.maxIterMarginal is not None): maxIterMarginalold = self.maxIterMarginal else: maxIterMarginalold = -1 self._maxIterMarginal = maxIterMarginal self._approxParameters["maxIterMarginal"] = self.maxIterMarginal if maxIterMarginalold != self.maxIterMarginal: self.resetSamples() def resetSamples(self): """Reset samples.""" super().resetSamples() if not hasattr(self, "_temporaryPivot"): self._mus = emptyParameterList() self.musMarginal = emptyParameterList() if hasattr(self, "samplerMarginalGrid"): self.samplerMarginalGrid.reset() if hasattr(self, "samplingEngine") and self.samplingEngine is not None: self.samplingEngine.resetHistory() def errorEstimatorMarginal(self, return_max : bool = False) -> Np1D: vbMng(self, "INIT", "Matching poles.", 10) self.trainedModel.initializeFromRational( self.HFEngine, self.matchingWeight, self.POD == PODGlobal, self.approx_state) vbMng(self, "DEL", "Done matching poles.", 10) self._finalizeMarginalization() _tMdataFull = copy(self.trainedModel.data) vbMng(self.trainedModel, "INIT", "Evaluating error estimator at mu = {}.".format( self.trainedModel.data.musMarginal), 10) err = np.zeros(len(self.trainedModel.data.musMarginal)) if len(err) <= 1: err[:] = np.inf else: if self.cutOffToleranceError == "AUTO": cutOffTolErr = self.cutOffTolerance else: cutOffTolErr = self.cutOffToleranceError if not hasattr(self, "_MMarginal_isauto"): if not hasattr(self, "_MMarginalOriginal"): self._MMarginalOriginal = self.MMarginal self.MMarginal = self._MMarginalOriginal _musMExcl = None self.verbosity -= 35 self.trainedModel.verbosity -= 35 foci = self.samplerPivot.normalFoci() ground = self.samplerPivot.groundPotential() for j in range(len(err)): jEff = j - (j > 0) muTest = self.trainedModel.data.musMarginal[jEff] polesEx = self.trainedModel.data.HIs[jEff].poles idxExEff = np.where(potential(polesEx, foci) - ground <= cutOffTolErr * ground)[0] polesEx = polesEx[idxExEff] if self.matchingWeightError != 0: resEx = self.trainedModel.data.HIs[jEff].coeffs[idxExEff] else: resEx = None if j > 0: self.musMarginal.insert(_musMExcl, j - 1) _musMExcl = self.musMarginal[j] self.musMarginal.pop(j) if len(polesEx) == 0: continue self.trainedModel.updateEffectiveSamples( self.HFEngine, [j], self.matchingWeight, self.POD == PODGlobal, self.approx_state) self._reduceDegreeNNoWarn = 1 self._finalizeMarginalization() polesAp = self.trainedModel.interpolateMarginalPoles(muTest)[ ..., 0] idxApEff = np.where(potential(polesAp, foci) - ground <= cutOffTolErr * ground)[0] polesAp = polesAp[idxApEff] if self.matchingWeightError != 0: resAp = self.trainedModel.interpolateMarginalCoeffs( muTest)[idxApEff, :, 0] if self.POD != PODGlobal: resEx = self.trainedModel.data.projMat.dot(resEx.T) resAp = self.trainedModel.data.projMat.dot(resAp.T) else: resAp = None dist = chordalMetricAdjusted( polesEx, polesAp, self.matchingWeightError, resEx, resAp, self.HFEngine, self.approx_state) pmR, pmC = pointMatching(dist) err[j] = np.mean(dist[pmR, pmC]) self.trainedModel.updateEffectiveSamples(self.HFEngine, None, self.matchingWeight, self.POD == PODGlobal, self.approx_state) if not hasattr(self, "_MMarginal_isauto"): self.MMarginal = self._MMarginalOriginal self.musMarginal.append(_musMExcl) self.verbosity += 35 self.trainedModel.verbosity += 35 self.trainedModel.data = _tMdataFull del self._reduceDegreeNNoWarn vbMng(self.trainedModel, "DEL", "Done evaluating error estimator", 10) if not return_max: return err idxMaxEst = np.where(err > self.greedyTolMarginal)[0] return err, idxMaxEst, err[idxMaxEst] def plotEstimatorMarginal(self, est:Np1D, idxMax:List[int], estMax:List[float]): if not (np.any(np.isnan(est)) or np.any(np.isinf(est))): fig = plt.figure(figsize = plt.figaspect(1. / self.nparMarginal)) for jpar in range(self.nparMarginal): ax = fig.add_subplot(1, self.nparMarginal, 1 + jpar) musre = copy(self.trainedModel.data.musMarginal.re.data) errCP = copy(est) idx = np.delete(np.arange(self.nparMarginal), jpar) while len(musre) > 0: if self.nparMarginal == 1: currIdx = np.arange(len(musre)) else: currIdx = np.where(np.isclose(np.sum( np.abs(musre[:, idx] - musre[0, idx]), 1), 0.))[0] currIdxSorted = currIdx[np.argsort(musre[currIdx, jpar])] ax.semilogy(musre[currIdxSorted, jpar], errCP[currIdxSorted], 'k.-', linewidth = 1) musre = np.delete(musre, currIdx, 0) errCP = np.delete(errCP, currIdx) ax.semilogy(self.musMarginal.re(jpar), (self.greedyTolMarginal,) * len(self.musMarginal), '*m') if len(idxMax) > 0 and estMax is not None: ax.semilogy(self.trainedModel.data.musMarginal.re( idxMax, jpar), estMax, 'xr') ax.grid() plt.tight_layout() plt.show() def _addMarginalSample(self, mus:paramList): mus = checkParameterList(mus, self.nparMarginal)[0] if len(mus) == 0: return nmus = len(mus) vbMng(self, "MAIN", ("Adding marginal sample point{} no. {}{} at {} to training " "set.").format("s" * (nmus > 1), len(self.musMarginal) + 1, "--{}".format(len(self.musMarginal) + nmus) * (nmus > 1), mus), 2) self.musMarginal.append(mus) self.setupApproxPivoted(mus) self._SMarginal = len(self.musMarginal) self._approxParameters["SMarginal"] = self.SMarginal def greedyNextSampleMarginal(self, muidx:int, plotEst : str = "NONE") \ -> Tuple[Np1D, int, float, paramVal]: RROMPyAssert(self._mode, message = "Cannot add greedy sample.") idxAdded = self.samplerMarginalGrid.refine(muidx) self._addMarginalSample(self.samplerMarginalGrid.points[idxAdded]) errorEstTest, muidx, maxErrorEst = self.errorEstimatorMarginal(True) if plotEst == "ALL": self.plotEstimatorMarginal(errorEstTest, muidx, maxErrorEst) return (errorEstTest, muidx, maxErrorEst, self.samplerMarginalGrid.points[muidx]) def _preliminaryTrainingMarginal(self): """Initialize starting snapshots of solution map.""" RROMPyAssert(self._mode, message = "Cannot start greedy algorithm.") if np.sum(self.samplingEngine.nsamples) > 0: return self.resetSamples() idx = [0] while self.samplerMarginalGrid.npoints < self.SMarginal: idx = self.samplerMarginalGrid.refine(idx) self._addMarginalSample(self.samplerMarginalGrid.points) def setupApproxPivoted(self, mu:paramVal) -> int: if self.checkComputedApproxPivoted(): return -1 RROMPyAssert(self._mode, message = "Cannot setup approximant.") raise RROMPyException("Must override.") def setupApprox(self, plotEst : str = "NONE") -> int: """Compute greedy snapshots of solution map.""" if self.checkComputedApprox(): return -1 RROMPyAssert(self._mode, message = "Cannot start greedy algorithm.") vbMng(self, "INIT", "Starting computation of snapshots.", 2) self._preliminaryTrainingMarginal() muidx, max2ErrorEst, firstGreedyIter = [], np.inf, True while firstGreedyIter or (max2ErrorEst > self.greedyTolMarginal and self.samplerMarginalGrid.npoints < self.maxIterMarginal): errorEstTest, muidx, maxErrorEst, mu = \ self.greedyNextSampleMarginal(muidx, plotEst) if len(maxErrorEst) > 0: max2ErrorEst = np.max(maxErrorEst) vbMng(self, "MAIN", ("Uniform testing error estimate " "{:.4e}.").format(max2ErrorEst), 2) else: max2ErrorEst = 0. firstGreedyIter = False if plotEst == "LAST": self.plotEstimatorMarginal(errorEstTest, muidx, maxErrorEst) vbMng(self, "DEL", ("Done computing snapshots (final snapshot count: " "{}).").format(np.sum(self.samplingEngine.nsamples)), 2) return 0 def checkComputedApprox(self) -> bool: return (super().checkComputedApprox() and len(self.mus) == len(self.trainedModel.data.mus)) def checkComputedApproxPivoted(self) -> bool: return (super().checkComputedApprox() and len(self.musMarginal) == len(self.trainedModel.data.musMarginal)) 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 index ee7c925..6974400 100644 --- a/rrompy/reduction_methods/pivoted/greedy/rational_interpolant_greedy_pivoted_greedy.py +++ b/rrompy/reduction_methods/pivoted/greedy/rational_interpolant_greedy_pivoted_greedy.py @@ -1,374 +1,374 @@ #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 from rrompy.reduction_methods.pivoted.generic_pivoted_approximant import ( 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; - 'cutOffKind': kind of cut off strategy; available values include 'SOFT' and 'HARD'; defaults to 'HARD'; - 'matchingWeightError': weight for pole matching optimization in error estimation; defaults to 0; - 'cutOffToleranceError': tolerance for ignoring parasitic poles in error estimation; defaults to 'AUTO', i.e. cutOffTolerance; - 'S': total number of pivot samples current approximant relies upon; - 'samplerPivot': pivot sample point generator; - 'SMarginal': number of starting marginal samples; - '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; - '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 'AUTO', i.e. maximum allowed; - '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; - 'cutOffKind': kind of cut off strategy; - '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; - '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. cutOffKind: Kind of cut off strategy. matchingWeightError: Weight for pole matching optimization in error estimation. cutOffToleranceError: Tolerance for ignoring parasitic poles in error estimation. S: Total number of pivot samples current approximant relies upon. samplerPivot: Pivot sample point generator. SMarginal: Total number of marginal samples current approximant relies upon. 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. 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. + muBounds: list of bounds for pivot parameter values. muBoundsMarginal: list of bounds for marginal parameter values. samplingEngine: Sampling engine. uHF: High fidelity solution(s) with parameter(s) lastSolvedHF as sampleList. lastSolvedHF: Parameter(s) corresponding to last computed high fidelity solution(s) as parameterList. uApproxReduced: Reduced approximate solution(s) with parameter(s) lastSolvedApprox as sampleList. lastSolvedApproxReduced: Parameter(s) corresponding to last computed reduced approximate solution(s) as parameterList. uApprox: Approximate solution(s) with parameter(s) lastSolvedApprox as sampleList. lastSolvedApprox: Parameter(s) corresponding to last computed approximate solution(s) as parameterList. """ @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) while len(musPivot) > self.S: musPivot.pop() muTestBasePivot = self.samplerPivot.generatePoints(self.nTestPoints, False) idxPop = pruneSamples( muTestBasePivot ** self.HFEngine.rescalingExp[self.directionPivot[0]], musPivot ** self.HFEngine.rescalingExp[self.directionPivot[0]], 1e-10 * self.scaleFactor[0]) muTestBasePivot.pop(idxPop) 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 self.checkComputedApprox(): return -1 if '_' not in plotEst: plotEst = plotEst + "_NONE" plotEstM, self._plotEstPivot = plotEst.split("_") val = super().setupApprox(plotEstM) return val diff --git a/rrompy/reduction_methods/pivoted/greedy/rational_interpolant_pivoted_greedy.py b/rrompy/reduction_methods/pivoted/greedy/rational_interpolant_pivoted_greedy.py index 64ac2f9..60cc759 100644 --- a/rrompy/reduction_methods/pivoted/greedy/rational_interpolant_pivoted_greedy.py +++ b/rrompy/reduction_methods/pivoted/greedy/rational_interpolant_pivoted_greedy.py @@ -1,320 +1,320 @@ # 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 import RationalInterpolant from rrompy.reduction_methods.pivoted import RationalInterpolantPivoted from rrompy.reduction_methods.pivoted.generic_pivoted_approximant import ( PODGlobal) from rrompy.utilities.base.types import paramList from rrompy.utilities.base import verbosityManager as vbMng from rrompy.utilities.exception_manager import RROMPyAssert from rrompy.parameter import checkParameterList, emptyParameterList __all__ = ['RationalInterpolantPivotedGreedy'] class RationalInterpolantPivotedGreedy(GenericPivotedGreedyApproximant, RationalInterpolantPivoted): """ ROM 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; - 'cutOffKind': kind of cut off strategy; available values include 'SOFT' and 'HARD'; defaults to 'HARD'; - 'matchingWeightError': weight for pole matching optimization in error estimation; defaults to 0; - 'cutOffToleranceError': tolerance for ignoring parasitic poles in error estimation; defaults to 'AUTO', i.e. cutOffTolerance; - 'S': total number of pivot samples current approximant relies upon; - 'samplerPivot': pivot sample point generator; - 'SMarginal': number of starting marginal samples; - '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'; - 'M': degree of rational interpolant numerator; defaults to 'AUTO', i.e. maximum allowed; - 'N': degree of rational interpolant denominator; defaults to 'AUTO', i.e. maximum allowed; - 'MMarginal': degree of marginal interpolant; defaults to 'AUTO', i.e. maximum allowed; - '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'; - 'radialDirectionalWeights': radial basis weights for pivot numerator; defaults to 1; - 'radialDirectionalWeightsMarginal': radial basis weights for marginal interpolant; defaults to 1; - 'nNearestNeighbor': number of pivot nearest neighbors considered if polybasis allows; 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; - 'correctorForce': whether corrector should forcefully delete bad poles; defaults to False; - 'correctorTol': tolerance for corrector step; defaults to 0., i.e. no bad poles; - 'correctorMaxIter': maximum number of corrector iterations; defaults to 1. Defaults to empty dict. approx_state(optional): Whether to approximate state. Defaults to False. verbosity(optional): Verbosity level. Defaults to 10. Attributes: HFEngine: HF problem solver. mu0: Default parameter. directionPivot: Pivot components. mus: Array of snapshot parameters. musPivot: Array of pivot snapshot parameters. musMarginal: Array of marginal snapshot parameters. approxParameters: Dictionary containing values for main parameters of approximant. Recognized keys are in parameterList. parameterListSoft: Recognized keys of soft approximant parameters: - 'POD': whether to compute POD of snapshots; - 'matchingWeight': weight for pole matching optimization; - 'cutOffTolerance': tolerance for ignoring parasitic poles; - 'cutOffKind': kind of cut off strategy; - '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; - 'M': degree of rational interpolant numerator; - 'N': degree of rational interpolant denominator; - '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; - '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; - '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. cutOffKind: Kind of cut off strategy. matchingWeightError: Weight for pole matching optimization in error estimation. cutOffToleranceError: Tolerance for ignoring parasitic poles in error estimation. S: Total number of pivot samples current approximant relies upon. samplerPivot: Pivot sample point generator. SMarginal: Total number of marginal samples current approximant relies upon. samplerMarginalGrid: Marginal sample point generator via sparse grid. polybasis: Type of polynomial basis for pivot interpolation. polybasisMarginal: Type of polynomial basis for marginal interpolation. M: Degree of rational interpolant numerator. N: Degree of rational interpolant denominator. 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. 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. - muBoundsPivot: list of bounds for pivot parameter values. + muBounds: list of bounds for pivot parameter values. muBoundsMarginal: list of bounds for marginal parameter values. samplingEngine: Sampling engine. uHF: High fidelity solution(s) with parameter(s) lastSolvedHF as sampleList. lastSolvedHF: Parameter(s) corresponding to last computed high fidelity solution(s) as parameterList. uApproxReduced: Reduced approximate solution(s) with parameter(s) lastSolvedApprox as sampleList. lastSolvedApproxReduced: Parameter(s) corresponding to last computed reduced approximate solution(s) as parameterList. uApprox: Approximate solution(s) with parameter(s) lastSolvedApprox as sampleList. lastSolvedApprox: Parameter(s) corresponding to last computed approximate solution(s) as parameterList. """ def _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 computeSnapshots(self): """Compute snapshots of solution map.""" RROMPyAssert(self._mode, message = "Cannot start snapshot computation.") vbMng(self, "INIT", "Starting computation of snapshots.", 5) self.musPivot = self.samplerPivot.generatePoints(self.S) while len(self.musPivot) > self.S: self.musPivot.pop() self.mus = emptyParameterList() self.mus.reset((self.S, self.HFEngine.npar)) self.samplingEngine.resetHistory() for k in range(self.S): self.mus.data[k, self.directionPivot] = self.musPivot[k].data self.mus.data[k, self.directionMarginal] = self.musMargLoc[-1].data self.samplingEngine.iterSample(self.mus) vbMng(self, "DEL", "Done computing snapshots.", 5) self._m_selfmus = copy(self.mus) self._mus = self.musPivot self._m_mu0 = copy(self.mu0) self._m_HFErescalingExp = copy(self.HFEngine.rescalingExp) self._mu0 = checkParameterList(self.mu0(self.directionPivot), 1)[0] self.HFEngine.rescalingExp = [self.HFEngine.rescalingExp[ self.directionPivot[0]]] def setupApproxPivoted(self, mus:paramList) -> int: if self.checkComputedApproxPivoted(): return -1 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 for j, mu in enumerate(mus): RationalInterpolant.setupSampling(self) self.trainedModel = None self.musMargLoc += [mu] RationalInterpolant.setupApprox(self) 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 self.samplingEngs[j] = copy(self.samplingEngine) Qs[j] = copy(self.trainedModel.data.Q) Ps[j] = copy(self.trainedModel.data.P) self.scaleFactor = self._scaleFactorOldPivot del self._scaleFactorOldPivot, self._temporaryPivot self._finalizeSnapshots() del self._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.verbosity += 15 vbMng(self, "DEL", "Done setting up approximant.", 10) return 0 diff --git a/rrompy/reduction_methods/pivoted/rational_interpolant_greedy_pivoted.py b/rrompy/reduction_methods/pivoted/rational_interpolant_greedy_pivoted.py index bb1cc6d..b5e106d 100644 --- a/rrompy/reduction_methods/pivoted/rational_interpolant_greedy_pivoted.py +++ b/rrompy/reduction_methods/pivoted/rational_interpolant_greedy_pivoted.py @@ -1,436 +1,436 @@ # 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, PODGlobal from rrompy.reduction_methods.standard.greedy.rational_interpolant_greedy \ import RationalInterpolantGreedy from rrompy.reduction_methods.standard.greedy.generic_greedy_approximant \ import pruneSamples from rrompy.utilities.base.types import Np1D from rrompy.utilities.base import verbosityManager as vbMng from rrompy.utilities.numerical import dot from rrompy.utilities.numerical.degree import totalDegreeN from rrompy.utilities.poly_fitting.polynomial import polyvander as pv from rrompy.utilities.exception_manager import RROMPyAssert 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; - 'cutOffKind': kind of cut off strategy; available values include 'SOFT' and 'HARD'; defaults to 'HARD'; - '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; - '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 'AUTO', i.e. maximum allowed; - '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; - 'cutOffKind': kind of cut off strategy; - '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; - '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; - '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. cutOffKind: Kind of cut off strategy. S: Total number of pivot samples current approximant relies upon. samplerPivot: Pivot sample point generator. SMarginal: Total number of marginal samples current approximant relies upon. samplerMarginal: Marginal sample point generator. polybasis: Type of polynomial basis for pivot interpolation. polybasisMarginal: Type of polynomial basis for marginal interpolation. greedyTol: uniform error tolerance for greedy algorithm. collinearityTol: Collinearity tolerance for greedy algorithm. maxIter: maximum number of greedy steps. nTestPoints: number of starting training points. trainSetGenerator: training sample points generator. errorEstimatorKind: kind of error estimator. 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. + muBounds: list of bounds for pivot parameter values. muBoundsMarginal: list of bounds for marginal parameter values. samplingEngine: Sampling engine. uHF: High fidelity solution(s) with parameter(s) lastSolvedHF as sampleList. lastSolvedHF: Parameter(s) corresponding to last computed high fidelity solution(s) as parameterList. uApproxReduced: Reduced approximate solution(s) with parameter(s) lastSolvedApprox as sampleList. lastSolvedApproxReduced: Parameter(s) corresponding to last computed reduced approximate solution(s) as parameterList. uApprox: Approximate solution(s) with parameter(s) lastSolvedApprox as sampleList. lastSolvedApprox: Parameter(s) corresponding to last computed approximate solution(s) as parameterList. Q: Numpy 1D vector containing complex coefficients of approximant denominator. P: Numpy 2D vector whose columns are FE dofs of coefficients of approximant numerator. """ def __init__(self, *args, **kwargs): self._preInit() self._addParametersToList(toBeExcluded = ["sampler"]) super().__init__(*args, **kwargs) self._postInit() @property def tModelType(self): if hasattr(self, "_temporaryPivot"): return RationalInterpolantGreedy.tModelType.fget(self) return super().tModelType @property def polybasis0(self): if "_" in self.polybasis: return self.polybasis.split("_")[0] return self.polybasis 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 setupOK > 0: err = np.empty(len(mus)) err[:] = np.nan if not return_max: return err return err, - setupOK, np.nan self._marginalizeTrainedModel(True) errRes = super().errorEstimator(mus, return_max) self._marginalizeTrainedModel(False) return errRes def _preliminaryTraining(self): """Initialize starting snapshots of solution map.""" RROMPyAssert(self._mode, message = "Cannot start greedy algorithm.") S = self.S self.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) while len(musPivot) > self.S: musPivot.pop() muTestPivot = self.samplerPivot.generatePoints(self.nTestPoints, False) idxPop = pruneSamples(muTestPivot ** self.HFEngine.rescalingExp[ self.directionPivot[0]], musPivot ** self.HFEngine.rescalingExp[ self.directionPivot[0]], 1e-10 * self.scaleFactor[0]) self.mus = emptyParameterList() self.mus.reset((self.S, self.npar + len(self.musMargLoc))) muTestBase = emptyParameterList() muTestBase.reset((len(muTestPivot), self.npar + len(self.musMargLoc))) for k in range(self.S): self.mus.data[k, self.directionPivot] = musPivot[k].data self.mus.data[k, self.directionMarginal] = self.musMargLoc.data for k in range(len(muTestPivot)): muTestBase.data[k, self.directionPivot] = muTestPivot[k].data muTestBase.data[k, self.directionMarginal] = self.musMargLoc.data muTestBase.pop(idxPop) muLast = copy(self.mus[-1]) self.mus.pop() if len(self.mus) > 0: vbMng(self, "MAIN", ("Adding first {} sample point{} at {} to training " "set.").format(self.S - 1, "" + "s" * (self.S > 2), self.mus), 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 _finalizeSnapshots(self): self.setupSampling() self.samplingEngine.resetHistory(len(self.musMarginal)) for j in range(len(self.musMarginal)): self.samplingEngine.setsample(self.samplingEngs[j].samples, j, False) self.samplingEngine.mus[j] = copy(self.samplingEngs[j].mus) self.samplingEngine.musMarginal[j] = copy(self.musMarginal[j]) self.samplingEngine.nsamples[j] = self.samplingEngs[j].nsamples if self.POD: self.samplingEngine.RPOD[j] = self.samplingEngs[j].RPOD self.samplingEngine.samples_full[j].data = ( self.samplingEngs[j].samples_full.data) if self.POD == PODGlobal: self.samplingEngine.coalesceSamples(self.interpRcondMarginal) else: self.samplingEngine.coalesceSamples() def setupApprox(self, *args, **kwargs) -> int: """Compute rational interpolant.""" if self.checkComputedApprox(): return -1 RROMPyAssert(self._mode, message = "Cannot setup approximant.") vbMng(self, "INIT", "Setting up {}.". format(self.name()), 5) self.musMarginal = self.samplerMarginal.generatePoints(self.SMarginal) while len(self.musMarginal) > self.SMarginal: self.musMarginal.pop() S0 = copy(self.S) Qs, Ps = [None] * len(self.musMarginal), [None] * len(self.musMarginal) self.samplingEngs = [None] * len(self.musMarginal) self.computeScaleFactor() self._scaleFactorOldPivot = copy(self.scaleFactor) self.scaleFactor = self.scaleFactorPivot self._temporaryPivot = 1 for j in range(len(self.musMarginal)): self._S = S0 self.musMargLoc = self.musMarginal[j] RationalInterpolantGreedy.setupSampling(self) self.trainedModel = None self.verbosity -= 5 self.samplingEngine.verbosity -= 5 super().setupApprox(*args, **kwargs) self.verbosity += 5 self.samplingEngine.verbosity += 5 self.samplingEngs[j] = copy(self.samplingEngine) Qs[j] = copy(self.trainedModel.data.Q) Ps[j] = copy(self.trainedModel.data.P) self.scaleFactor = self._scaleFactorOldPivot del self._scaleFactorOldPivot, self._temporaryPivot self._finalizeSnapshots() del self.musMargLoc, self.samplingEngs self._mus = self.samplingEngine.musCoalesced padLeft = 0 if self.POD != PODGlobal: padRight = self.samplingEngine.nsamplesTot for j in range(len(self.musMarginal)): nsj = self.samplingEngine.nsamples[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 = 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.Qs, self.trainedModel.data.Ps = Qs, Ps vbMng(self, "INIT", "Matching poles.", 10) self.trainedModel.initializeFromRational( self.HFEngine, self.matchingWeight, self.POD == PODGlobal, self.approx_state) vbMng(self, "DEL", "Done matching poles.", 10) self._finalizeMarginalization() vbMng(self, "DEL", "Done setting up approximant.", 5) return 0 diff --git a/rrompy/reduction_methods/pivoted/rational_interpolant_pivoted.py b/rrompy/reduction_methods/pivoted/rational_interpolant_pivoted.py index b33b33b..aff0f4b 100644 --- a/rrompy/reduction_methods/pivoted/rational_interpolant_pivoted.py +++ b/rrompy/reduction_methods/pivoted/rational_interpolant_pivoted.py @@ -1,362 +1,362 @@ # 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, PODGlobal from rrompy.reduction_methods.standard.rational_interpolant import ( RationalInterpolant) from rrompy.utilities.base import verbosityManager as vbMng from rrompy.utilities.numerical import dot from rrompy.utilities.numerical.hash_derivative import nextDerivativeIndices from rrompy.utilities.exception_manager import RROMPyAssert, RROMPyWarning from rrompy.parameter import emptyParameterList __all__ = ['RationalInterpolantPivoted'] class RationalInterpolantPivoted(GenericPivotedApproximant, RationalInterpolant): """ 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; - 'cutOffKind': kind of cut off strategy; available values include 'SOFT' and 'HARD'; defaults to 'HARD'; - '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'; - 'M': degree of rational interpolant numerator; defaults to 'AUTO', i.e. maximum allowed; - 'N': degree of rational interpolant denominator; defaults to 'AUTO', i.e. maximum allowed; - 'MMarginal': degree of marginal interpolant; defaults to 'AUTO', i.e. maximum allowed; - 'polydegreetypeMarginal': type of polynomial degree for marginal; defaults to 'TOTAL'; - 'radialDirectionalWeights': radial basis weights for pivot numerator; defaults to 1; - 'radialDirectionalWeightsMarginal': radial basis weights for marginal interpolant; defaults to 1; - 'nNearestNeighbor': number of pivot nearest neighbors considered if polybasis allows; 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; - 'correctorForce': whether corrector should forcefully delete bad poles; defaults to False; - 'correctorTol': tolerance for corrector step; defaults to 0., i.e. no bad poles; - 'correctorMaxIter': maximum number of corrector iterations; defaults to 1. Defaults to empty dict. approx_state(optional): Whether to approximate state. Defaults to False. verbosity(optional): Verbosity level. Defaults to 10. Attributes: HFEngine: HF problem solver. mu0: Default parameter. directionPivot: Pivot components. mus: Array of snapshot parameters. musPivot: Array of pivot snapshot parameters. musMarginal: Array of marginal snapshot parameters. approxParameters: Dictionary containing values for main parameters of approximant. Recognized keys are in parameterList. parameterListSoft: Recognized keys of soft approximant parameters: - 'POD': whether to compute POD of snapshots; - 'matchingWeight': weight for pole matching optimization; - 'cutOffTolerance': tolerance for ignoring parasitic poles; - 'cutOffKind': kind of cut off strategy; - 'polybasis': type of polynomial basis for pivot interpolation; - 'polybasisMarginal': type of polynomial basis for marginal interpolation; - 'M': degree of rational interpolant numerator; - 'N': degree of rational interpolant denominator; - '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. cutOffKind: Kind of cut off strategy. S: Total number of pivot samples current approximant relies upon. samplerPivot: Pivot sample point generator. SMarginal: Total number of marginal samples current approximant relies upon. samplerMarginal: Marginal sample point generator. polybasis: Type of polynomial basis for pivot interpolation. polybasisMarginal: Type of polynomial basis for marginal interpolation. M: Numerator degree of approximant. N: Denominator degree of approximant. 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. - muBoundsPivot: list of bounds for pivot parameter values. + muBounds: list of bounds for pivot parameter values. muBoundsMarginal: list of bounds for marginal parameter values. samplingEngine: Sampling engine. uHF: High fidelity solution(s) with parameter(s) lastSolvedHF as sampleList. lastSolvedHF: Parameter(s) corresponding to last computed high fidelity solution(s) as parameterList. uApproxReduced: Reduced approximate solution(s) with parameter(s) lastSolvedApprox as sampleList. lastSolvedApproxReduced: Parameter(s) corresponding to last computed reduced approximate solution(s) as parameterList. uApprox: Approximate solution(s) with parameter(s) lastSolvedApprox as sampleList. lastSolvedApprox: Parameter(s) corresponding to last computed approximate solution(s) as parameterList. Q: Numpy 1D vector containing complex coefficients of approximant denominator. P: Numpy 2D vector whose columns are FE dofs of coefficients of approximant numerator. """ def __init__(self, *args, **kwargs): self._preInit() self._addParametersToList(toBeExcluded = ["polydegreetype", "sampler"]) super().__init__(*args, **kwargs) self._postInit() @property def polydegreetype(self): """Value of polydegreetype.""" return "TOTAL" @polydegreetype.setter def polydegreetype(self, polydegreetype): RROMPyWarning(("polydegreetype is used just to simplify inheritance, " "and its value cannot be changed from 'TOTAL'.")) @property def polybasis0(self): if "_" in self.polybasis: return self.polybasis.split("_")[0] return self.polybasis @property def correctorTol(self): """Value of correctorTol.""" return self._correctorTol @correctorTol.setter def correctorTol(self, correctorTol): if correctorTol < 0. or (correctorTol > 0. and self.nparPivot > 1): RROMPyWarning(("Overriding prescribed corrector tolerance " "to 0.")) correctorTol = 0. self._correctorTol = correctorTol self._approxParameters["correctorTol"] = self.correctorTol @property def correctorMaxIter(self): """Value of correctorMaxIter.""" return self._correctorMaxIter @correctorMaxIter.setter def correctorMaxIter(self, correctorMaxIter): if correctorMaxIter < 1 or (correctorMaxIter > 1 and self.nparPivot > 1): RROMPyWarning(("Overriding prescribed max number of corrector " "iterations to 1.")) correctorMaxIter = 1 self._correctorMaxIter = correctorMaxIter self._approxParameters["correctorMaxIter"] = self.correctorMaxIter def _setupInterpolationIndices(self): """Setup parameters for polyvander.""" RROMPyAssert(self._mode, message = "Cannot setup interpolation indices.") if (self._musUniqueCN is None or len(self._reorder) != len(self.musPivot)): try: muPC = self.trainedModel.centerNormalizePivot(self.musPivot) except: muPC = self.trainedModel.centerNormalize(self.musPivot) self._musUniqueCN, musIdxsTo, musIdxs, musCount = (muPC.unique( return_index = True, return_inverse = True, return_counts = True)) self._musUnique = self.musPivot[musIdxsTo] self._derIdxs = [None] * len(self._musUniqueCN) self._reorder = np.empty(len(musIdxs), dtype = int) filled = 0 for j, cnt in enumerate(musCount): self._derIdxs[j] = nextDerivativeIndices([], self.nparPivot, cnt) jIdx = np.nonzero(musIdxs == j)[0] self._reorder[jIdx] = np.arange(filled, filled + cnt) filled += cnt def computeSnapshots(self): """Compute snapshots of solution map.""" RROMPyAssert(self._mode, message = "Cannot start snapshot computation.") if self.samplingEngine.nsamplesTot != self.S * self.SMarginal: self.computeScaleFactor() self.resetSamples() vbMng(self, "INIT", "Starting computation of snapshots.", 5) self.musPivot = self.samplerPivot.generatePoints(self.S) while len(self.musPivot) > self.S: self.musPivot.pop() self.musMarginal = self.samplerMarginal.generatePoints( self.SMarginal) while len(self.musMarginal) > self.SMarginal: self.musMarginal.pop() self.mus = emptyParameterList() self.mus.reset((self.S * self.SMarginal, self.HFEngine.npar)) self.samplingEngine.resetHistory(self.SMarginal) for j, muMarg in enumerate(self.musMarginal): for k in range(j * self.S, (j + 1) * self.S): self.mus.data[k, self.directionPivot] = ( self.musPivot[k - j * self.S].data) self.mus.data[k, self.directionMarginal] = muMarg.data self.samplingEngine.iterSample(self.musPivot, self.musMarginal) self._finalizeSnapshots() vbMng(self, "DEL", "Done computing snapshots.", 5) def _finalizeSnapshots(self): if self.POD == PODGlobal: self.samplingEngine.coalesceSamples(self.interpRcondMarginal) else: self.samplingEngine.coalesceSamples() def setupApprox(self) -> int: """Compute rational interpolant.""" if self.checkComputedApprox(): return -1 RROMPyAssert(self._mode, message = "Cannot setup approximant.") vbMng(self, "INIT", "Setting up {}.". format(self.name()), 5) self.computeSnapshots() pMat = self.samplingEngine.samplesCoalesced.data pMatEff = dot(self.HFEngine.C, pMat) if self.approx_state else pMat if self.trainedModel is None: self.trainedModel = self.tModelType() self.trainedModel.verbosity = self.verbosity self.trainedModel.timestamp = self.timestamp datadict = {"mu0": self.mu0, "projMat": pMatEff, "scaleFactor": self.scaleFactor, "rescalingExp": self.HFEngine.rescalingExp, "directionPivot": self.directionPivot} self.trainedModel.data = self.initializeModelData(datadict)[0] else: self.trainedModel = self.trainedModel self.trainedModel.data.projMat = copy(pMatEff) N0 = copy(self.N) Qs, Ps = [None] * len(self.musMarginal), [None] * len(self.musMarginal) self._temporaryPivot = 1 padLeft = 0 if self.POD: self._RPODOldPivot = copy(self.samplingEngine.RPODCoalesced) else: self._samplesOldPivot = copy(self.samplingEngine.samples) padRight = self.samplingEngine.nsamplesTot self._scaleFactorOldPivot = copy(self.scaleFactor) self.scaleFactor = self.scaleFactorPivot for j in range(len(self.musMarginal)): self.N = N0 if self.POD: self.samplingEngine.RPOD = ( self._RPODOldPivot[:, padLeft : padLeft + self.S]) else: self.samplingEngine.samples = self._samplesOldPivot[j] padRight -= self.S self.verbosity -= 5 self._iterCorrector() self.verbosity += 5 Qs[j] = copy(self.trainedModel.data.Q) Ps[j] = copy(self.trainedModel.data.P) del self.trainedModel.data.Q, self.trainedModel.data.P if not self.POD: Ps[j].pad(padLeft, padRight) padLeft += self.S if self.POD: self.samplingEngine.RPODCoalesced = copy(self._RPODOldPivot) del self._RPODOldPivot else: self.samplingEngine.samples = copy(self._samplesOldPivot) del self._samplesOldPivot self.scaleFactor = self._scaleFactorOldPivot del self._temporaryPivot, self._scaleFactorOldPivot self.trainedModel.data.mus = copy(self.mus) self.trainedModel.data.musPivot = copy(self.musPivot) self.trainedModel.data.musMarginal = copy(self.musMarginal) 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 == PODGlobal, self.approx_state) vbMng(self, "DEL", "Done matching poles.", 10) self._finalizeMarginalization() vbMng(self, "DEL", "Done setting up approximant.", 5) return 0 diff --git a/rrompy/reduction_methods/standard/greedy/rational_interpolant_greedy.py b/rrompy/reduction_methods/standard/greedy/rational_interpolant_greedy.py index f5ca226..816460c 100644 --- a/rrompy/reduction_methods/standard/greedy/rational_interpolant_greedy.py +++ b/rrompy/reduction_methods/standard/greedy/rational_interpolant_greedy.py @@ -1,572 +1,576 @@ # 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 dot from rrompy.utilities.numerical.degree import totalDegreeN from rrompy.utilities.expression import expressionEvaluator from rrompy.reduction_methods.standard import RationalInterpolant from rrompy.utilities.base.types import Np1D, Tuple, 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 sampleList, 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; - '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', 'LOOK_AHEAD_OUTPUT', 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; - '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. 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", "LOOK_AHEAD_OUTPUT", "NONE"] def __init__(self, *args, **kwargs): self._preInit() self._addParametersToList(["errorEstimatorKind"], ["DISCREPANCY"], toBeExcluded = ["M", "N", "polydegreetype", "radialDirectionalWeights", "nNearestNeighbor"]) super().__init__(*args, **kwargs) if not self.approx_state and self.errorEstimatorKind not in [ "LOOK_AHEAD", "LOOK_AHEAD_OUTPUT", "NONE"]: raise RROMPyException(("Must compute greedy approximation of " "state, unless error estimator allows " "otherwise.")) self.M, self.N = ("AUTO",) * 2 self._postInit() @property def approx_state(self): """Value of approx_state.""" return self._approx_state @approx_state.setter def approx_state(self, approx_state): RationalInterpolant.approx_state.fset(self, approx_state) if (not self.approx_state and hasattr(self, "_errorEstimatorKind") and self.errorEstimatorKind not in [ "LOOK_AHEAD", "LOOK_AHEAD_OUTPUT", "NONE"]): raise RROMPyException(("Must compute greedy approximation of " "state, unless error estimator allows " "otherwise.")) @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.")) def _setMAuto(self): self.M = self.E def _setNAuto(self): self.N = self.E @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 if (self.errorEstimatorKind not in [ "LOOK_AHEAD", "LOOK_AHEAD_OUTPUT", "NONE"] and hasattr(self, "_approx_state") and not self.approx_state): raise RROMPyException(("Must compute greedy approximation of " "state, unless error estimator allows " "otherwise.")) def _polyvanderAuxiliary(self, mus, deg, *args): return pvT(mus, deg, *args) def getErrorEstimatorDiscrepancy(self, mus:Np1D) -> Np1D: """Discrepancy-based residual estimator.""" checkIfAffine(self.HFEngine, "apply discrepancy-based error 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) QTzero = np.where(QTest == 0.)[0] if len(QTzero) > 0: RROMPyWarning(("Adjusting estimator to avoid division by " "numerically zero denominator.")) QTest[QTzero] = np.finfo(np.complex).eps / (1. + self.N) 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) QTzero = np.where(QTrain == 0.)[0] if len(QTzero) > 0: RROMPyWarning(("Adjusting estimator to avoid division by " "numerically zero denominator.")) QTrain[QTzero] = np.finfo(np.complex).eps / (1. + self.N) 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, force_output : bool = False) \ -> Tuple[Np1D, List[int]]: """Residual estimator based on look-ahead idea.""" errTest, QTest, idxMaxEst = self._EIMStep(mus) _approx_state_old = self.approx_state if force_output and _approx_state_old: self._approx_state = False self.initEstimatorNormEngine() self._approx_state = _approx_state_old 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): uEx = self.samplingEngine.nextSample(mu) if hasattr(self.samplingEngine, "samples_full"): uEx = self.samplingEngine.samples_full[-1] if j == 0: solmu = emptySampleList() solmu.reset((len(uEx), len(mu_muTestSample)), dtype = uEx.dtype) solmu[j] = uEx if force_output and self.approx_state: solmu = sampleList(self.HFEngine.applyC(solmu.data)) app_muTestSample = sampleList(self.HFEngine.applyC( app_muTestSample.data)) 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] verb = self.trainedModel.verbosity self.trainedModel.verbosity = 0 QTest = self.trainedModel.getQVal(mus) QTzero = np.where(QTest == 0.)[0] if len(QTzero) > 0: RROMPyWarning(("Adjusting estimator to avoid division by " "numerically zero denominator.")) QTest[QTzero] = np.finfo(np.complex).eps / (1. + self.N) 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] try: coeffTest = np.linalg.solve(vanTrial, valuesTrial) except np.linalg.LinAlgError as e: raise RROMPyException(e) 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 setupOK > 0: err = np.empty(len(mus)) err[:] = np.nan if not return_max: return err return err, - setupOK, 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[: 10] == "LOOK_AHEAD": err, idxMaxEst = self.getErrorEstimatorLookAhead(mus, self.errorEstimatorKind == "LOOK_AHEAD_OUTPUT") 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[: 10] != "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 plotEstimator(self, *args, **kwargs): super().plotEstimator(*args, **kwargs) if self.errorEstimatorKind == "NONE": vbMng(self, "MAIN", ("Warning! Error estimator has been arbitrarily normalized " "before plotting."), 5) def greedyNextSample(self, *args, **kwargs) -> 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(*args, **kwargs) 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) -> int: """Compute rational interpolant.""" if self.checkComputedApprox(): return -1 RROMPyAssert(self._mode, message = "Cannot setup approximant.") self.verbosity -= 10 vbMng(self, "INIT", "Setting up local approximant.", 5) 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.trainedModel.data.mus = copy(self.mus) self.catchInstability = 2 unstable = False if self.E > 0: try: Q = self._setupDenominator()[0] except RROMPyException as RE: RROMPyWarning("Downgraded {}: {}".format(RE.__class__.__name__, RE)) vbMng(self, "DEL", "", 7) unstable = True else: Q = PI() Q.coeffs = np.ones((1,) * self.npar, dtype = np.complex) Q.npar = self.npar Q.polybasis = self.polybasis if not unstable: self.trainedModel.data.Q = copy(Q) try: P = copy(self._setupNumerator()) except RROMPyException as RE: RROMPyWarning("Downgraded {}: {}".format(RE.__class__.__name__, RE)) vbMng(self, "DEL", "", 7) unstable = True if not unstable: 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 1 * unstable + def setupApprox(self, plotEst : str = "NONE") -> int: + val = super().setupApprox(plotEst) + if val == 0: self._iterCorrector() + return val + 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) -