diff --git a/rrompy/reduction_methods/pivoted/generic_pivoted_approximant.py b/rrompy/reduction_methods/pivoted/generic_pivoted_approximant.py index 7cccc39..4a9ce36 100644 --- a/rrompy/reduction_methods/pivoted/generic_pivoted_approximant.py +++ b/rrompy/reduction_methods/pivoted/generic_pivoted_approximant.py @@ -1,563 +1,484 @@ # 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, - PolynomialInterpolator as PI) -from rrompy.utilities.poly_fitting.radial_basis import (polybases as rbpb, - RadialBasisInterpolator as RBI) +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, - MovingLeastSquaresInterpolator as MLSI) + polybases as mlspb) from rrompy.sampling.pivoted import (SamplingEnginePivoted, SamplingEnginePivotedPOD, SamplingEnginePivotedPODGlobal) from rrompy.utilities.base.types import paramList, ListAny from rrompy.utilities.base import verbosityManager as vbMng -from rrompy.utilities.numerical import reduceDegreeN, nextDerivativeIndices +from rrompy.utilities.numerical 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. 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", - "polybasisMarginal", "MMarginal", - "polydegreetypeMarginal", + "cutOffKind", "polybasisMarginal", + "MMarginal", "polydegreetypeMarginal", "radialDirectionalWeightsMarginal", "nNearestNeighborMarginal", "interpRcondMarginal"], - [1, np.inf, "MONOMIAL", "AUTO", "TOTAL", [1], - -1, -1], ["samplerPivot", "SMarginal", - "samplerMarginal"], [QSBase, [1], QSBase]) + [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 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 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.""" 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 resetSamples(self): - """Reset samples.""" - super().resetSamples() - self._musMUniqueCN = None - self._derMIdxs = None - self._reorderM = None - 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 _setupMarginalInterpolationIndices(self): - """Setup parameters for polyvander.""" - RROMPyAssert(self._mode, - message = "Cannot setup interpolation indices.") - if (self._musMUniqueCN is None - or len(self._reorderM) != len(self.musMarginal)): - self._musMUniqueCN, musMIdxsTo, musMIdxs, musMCount = ( - self.trainedModel.centerNormalizeMarginal(self.musMarginal)\ - .unique(return_index = True, return_inverse = True, - return_counts = True)) - self._musMUnique = self.musMarginal[musMIdxsTo] - self._derMIdxs = [None] * len(self._musMUniqueCN) - self._reorderM = np.empty(len(musMIdxs), dtype = int) - filled = 0 - for j, cnt in enumerate(musMCount): - self._derMIdxs[j] = nextDerivativeIndices([], - self.nparMarginal, cnt) - jIdx = np.nonzero(musMIdxs == j)[0] - self._reorderM[jIdx] = np.arange(filled, filled + cnt) - filled += cnt - - def _setupMarginalInterp(self): - """Compute marginal interpolator.""" - RROMPyAssert(self._mode, message = "Cannot setup numerator.") - vbMng(self, "INIT", "Starting computation of marginal interpolator.", - 7) - self._setupMarginalInterpolationIndices() - if hasattr(self, "radialDirectionalWeightsMarginal"): - rDWM = copy(self.radialDirectionalWeightsMarginal) - if hasattr(self, "_MMarginal_isauto"): - self._setMMarginalAuto() - MM = self.MMarginal - else: - MM = reduceDegreeN(self.MMarginal, len(self.musMarginal), - self.nparMarginal, self.polydegreetypeMarginal) - if MM < self.MMarginal: - if not hasattr(self, "_reduceDegreeNNoWarn"): - RROMPyWarning(("MMarginal too large compared to SMarginal. " - "Reducing MMarginal by " - "{}").format(self.MMarginal - MM)) - self.MMarginal = MM - mI = [] - for j in range(len(self.musMarginal)): - canonicalj = 1. * (np.arange(len(self.musMarginal)) == j) - while (self.MMarginal >= 0 - and (not hasattr(self, "radialDirectionalWeightsMarginal") - or self.radialDirectionalWeightsMarginal[0] - <= rDWM[0] * 2 ** 6)): - pParRest = [self.verbosity >= 5, - self.polydegreetypeMarginal == "TOTAL", - {"derIdxs": self._derMIdxs, - "reorder": self._reorderM, - "scl": np.power(self.scaleFactorMarginal, -1.)}] - if self.polybasisMarginal in ppb: - p = PI() - else: - pParRest = ([self.radialDirectionalWeightsMarginal] - + pParRest) - pParRest[-1]["nNearestNeighbor"] = ( - self.nNearestNeighborMarginal) - p = RBI() if self.polybasisMarginal in rbpb else MLSI() - if self.polybasisMarginal in ppb + rbpb: - pParRest += [{"rcond": self.interpRcondMarginal}] - wellCond, msg = p.setupByInterpolation(self._musMUniqueCN, - canonicalj, self.MMarginal, - self.polybasisMarginal, - *pParRest) - vbMng(self, "MAIN", msg, 5) - if wellCond: break - if self.polybasisMarginal in ppb: - vbMng(self, "MAIN", ("Polyfit is poorly conditioned. " - "Reducing MMarginal by 1."), 10) - self.MMarginal = self.MMarginal - 1 - else: - vbMng(self, "MAIN", ("Polyfit is poorly conditioned. " - "Multiplying " - "radialDirectionalWeightsMarginal by " - "2."), 10) - for j in range(self.nparMarginal): - self._radialDirectionalWeightsMarginal[j] *= 2. - if (self.MMarginal < 0 - or (hasattr(self, "radialDirectionalWeightsMarginal") - and self.radialDirectionalWeightsMarginal[0] > rDWM[0] * 2 ** 6)): - raise RROMPyException(("Instability in computation of " - "interpolant. Aborting.")) - mI = mI + [copy(p)] - if self.polybasisMarginal in ppb: - self.MMarginal = MM - else: - self.radialDirectionalWeightsMarginal = rDWM - vbMng(self, "DEL", "Done computing marginal interpolator.", 7) - return mI - def _finalizeMarginalization(self): - 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) - vbMng(self, "INIT", "Recompressing by cut-off.", 10) - msg = self.trainedModel.recompressByCutOff([-1., 1.], - self.cutOffTolerance) + vbMng(self, "INIT", "Recompressing by cut off.", 10) + msg = self.trainedModel.recompressByCutOff(self.cutOffTolerance, + self.cutOffKind, [-1., 1.]) 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.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 \ No newline at end of file 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 6b3937f..31607a6 100644 --- a/rrompy/reduction_methods/pivoted/greedy/generic_pivoted_greedy_approximant.py +++ b/rrompy/reduction_methods/pivoted/greedy/generic_pivoted_greedy_approximant.py @@ -1,425 +1,440 @@ # 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 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 import (pointMatching, chordalMetricAdjusted, cutOffDistance) from rrompy.utilities.exception_manager import RROMPyException, RROMPyAssert 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 np.inf; - 'S': total number of pivot samples current approximant relies upon; - 'samplerPivot': pivot sample point generator; - 'SMarginal': number of starting marginal samples; - 'samplerMarginalGrid': marginal sample point generator via sparse grid; - '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. 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., np.inf, 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): 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() 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 not hasattr(self, "_MMarginal_isauto"): if not hasattr(self, "_MMarginalOriginal"): self._MMarginalOriginal = self.MMarginal self.MMarginal = self._MMarginalOriginal _musMExcl = None self.verbosity -= 20 + self.trainedModel.verbosity -= 20 for j in range(len(err)): - jEff = j - (j > 0) + jEff = j - (j > 0) muTest = self.trainedModel.data.musMarginal[jEff] polesEx = self.trainedModel.data.HIs[jEff].poles idxExEff = np.where(cutOffDistance(polesEx) <= self.cutOffToleranceError)[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._musMUniqueCN = None - self.trainedModel.data.marginalInterp = ( - self._setupMarginalInterp()) - polesAp = self.trainedModel.interpolateMarginalPoles(muTest) + self._finalizeMarginalization() + self._reduceDegreeNNoWarn = 1 + polesAp = self.trainedModel.interpolateMarginalPoles(muTest)[ + ..., 0] + idxApEff = np.where(cutOffDistance(polesAp) <= + self.cutOffToleranceError)[0] + polesAp = polesAp[idxApEff] if self.matchingWeightError != 0: resAp = self.trainedModel.interpolateMarginalCoeffs( - muTest)[: len(polesAp)] + 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) - err[j] = np.mean(dist[np.arange(len(polesEx)), - pointMatching(dist)]) + 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) - self.musMarginal.append(_musMExcl) if not hasattr(self, "_MMarginal_isauto"): self.MMarginal = self._MMarginalOriginal + self.musMarginal.append(_musMExcl) self.verbosity += 20 - self._reduceDegreeNNoWarn = 1 - self.trainedModel.data.marginalInterp = self._setupMarginalInterp() + self.trainedModel.verbosity += 20 + self._finalizeMarginalization() 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 = [], np.inf + max2ErrorEst = np.inf + errorEstTest, muidx, maxErrorEst, mu = \ + self.greedyNextSampleMarginal([], plotEst) while (self.samplerMarginalGrid.npoints < self.maxIterMarginal and max2ErrorEst > self.greedyTolMarginal): 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. 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 82f735b..3951874 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,369 +1,373 @@ #Copyright (C) 2018 by the RROMPy authors # # This file is part of RROMPy. # # RROMPy is free software: you can redistribute it and/or modify # it under the terms of the GNU Lesser General Public License as published by # the Free Software Foundation, either version 3 of the License, or # (at your option) any later version. # # RROMPy is distributed in the hope that it will be useful, # but WITHOUT ANY WARRANTY; without even the implied warranty of # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the # GNU Lesser General Public License for more details. # # You should have received a copy of the GNU Lesser General Public License # along with RROMPy. If not, see . # from copy import deepcopy as copy import numpy as np from .generic_pivoted_greedy_approximant import GenericPivotedGreedyApproximant from rrompy.utilities.numerical import dot from rrompy.reduction_methods.standard.greedy import RationalInterpolantGreedy from rrompy.reduction_methods.standard.greedy.generic_greedy_approximant \ import pruneSamples from rrompy.reduction_methods.pivoted import (RationalInterpolantGreedyPivoted, PODGlobal) from rrompy.utilities.base.types import Np1D, Tuple, paramVal, paramList from rrompy.utilities.base import verbosityManager as vbMng from rrompy.utilities.exception_manager import RROMPyAssert from rrompy.parameter import emptyParameterList __all__ = ['RationalInterpolantGreedyPivotedGreedy'] class RationalInterpolantGreedyPivotedGreedy(GenericPivotedGreedyApproximant, RationalInterpolantGreedyPivoted): """ ROM greedy pivoted greedy rational interpolant computation for parametric problems. Args: HFEngine: HF problem solver. mu0(optional): Default parameter. Defaults to 0. directionPivot(optional): Pivot components. Defaults to [0]. approxParameters(optional): Dictionary containing values for main parameters of approximant. Recognized keys are: - 'POD': whether to compute POD of snapshots; defaults to True; - 'matchingWeight': weight for pole matching optimization; defaults to 1; - 'cutOffTolerance': tolerance for ignoring parasitic poles; defaults to np.inf; + - '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 np.inf; - 'S': total number of pivot samples current approximant relies upon; - 'samplerPivot': pivot sample point generator; - 'SMarginal': number of starting marginal samples; - 'samplerMarginalGrid': marginal sample point generator via sparse grid; - 'polybasis': type of polynomial basis for pivot interpolation; defaults to 'MONOMIAL'; - 'polybasisMarginal': type of polynomial basis for marginal interpolation; allowed values include 'MONOMIAL', 'CHEBYSHEV' and 'LEGENDRE'; defaults to 'MONOMIAL'; - 'greedyTol': uniform error tolerance for greedy algorithm; defaults to 1e-2; - 'collinearityTol': collinearity tolerance for greedy algorithm; defaults to 0.; - 'maxIter': maximum number of greedy steps; defaults to 1e2; - '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. muBoundsMarginal: list of bounds for marginal parameter values. samplingEngine: Sampling engine. uHF: High fidelity solution(s) with parameter(s) lastSolvedHF as sampleList. lastSolvedHF: Parameter(s) corresponding to last computed high fidelity solution(s) as parameterList. uApproxReduced: Reduced approximate solution(s) with parameter(s) lastSolvedApprox as sampleList. lastSolvedApproxReduced: Parameter(s) corresponding to last computed reduced approximate solution(s) as parameterList. uApprox: Approximate solution(s) with parameter(s) lastSolvedApprox as sampleList. lastSolvedApprox: Parameter(s) corresponding to last computed approximate solution(s) as parameterList. """ @property def sampleBatchSize(self): """Value of sampleBatchSize.""" return 1 @property def sampleBatchIdx(self): """Value of sampleBatchIdx.""" return self.S def _finalizeSnapshots(self): self.samplingEngine = self._samplingEngineOld for muM, sEN in zip(self.musMargLoc, self.samplingEngs): self.samplingEngine.samples += [sEN.samples] self.samplingEngine.nsamples += [sEN.nsamples] self.samplingEngine.mus += [sEN.mus] self.samplingEngine.musMarginal.append(muM) self.samplingEngine._derIdxs += [[(0,) * self.npar] for _ in range(sEN.nsamples)] if self.POD: self.samplingEngine.RPOD += [sEN.RPOD] self.samplingEngine.samples_full += [copy(sEN.samples_full)] if self.POD == PODGlobal: self.samplingEngine.coalesceSamples(self.interpRcondMarginal) else: self.samplingEngine.coalesceSamples() def greedyNextSample(self, muidx:int, plotEst : str = "NONE")\ -> Tuple[Np1D, int, float, paramVal]: """Compute next greedy snapshot of solution map.""" RROMPyAssert(self._mode, message = "Cannot add greedy sample.") mus = copy(self.muTest[muidx]) self.muTest.pop(muidx) for j, mu in enumerate(mus): vbMng(self, "MAIN", ("Adding sample point no. {} at {} to training " "set.").format(len(self.mus) + 1, mu), 2) self.mus.append(mu) self._S = len(self.mus) self._approxParameters["S"] = self.S if (self.samplingEngine.nsamples <= len(mus) - j - 1 or not np.allclose(mu, self.samplingEngine.mus.data[j - len(mus)])): self.samplingEngine.nextSample(mu) if self._isLastSampleCollinear(): vbMng(self, "MAIN", ("Collinearity above tolerance detected. Starting " "preemptive greedy loop termination."), 2) self._collinearityFlag = 1 errorEstTest = np.empty(len(self.muTest)) errorEstTest[:] = np.nan return errorEstTest, [-1], np.nan, np.nan errorEstTest, muidx, maxErrorEst = self.errorEstimator(self.muTest, True) if plotEst == "ALL": self.plotEstimator(errorEstTest, muidx, maxErrorEst) return errorEstTest, muidx, maxErrorEst, self.muTest[muidx] def _preliminaryTraining(self): """Initialize starting snapshots of solution map.""" RROMPyAssert(self._mode, message = "Cannot start greedy algorithm.") if self.samplingEngine.nsamples > 0: return self.resetSamples() musPivot = self.trainSetGenerator.generatePoints(self.S) musPivot.data = musPivot.data[: self.S] muTestBasePivot = self.samplerPivot.generatePoints(self.nTestPoints) idxPop = pruneSamples( muTestBasePivot ** self.HFEngine.rescalingExp[self.directionPivot[0]], musPivot ** self.HFEngine.rescalingExp[self.directionPivot[0]], 1e-10 * self.scaleFactor[0]) muTestBasePivot.pop(idxPop) muTestBasePivot = muTestBasePivot.sort() self.mus = emptyParameterList() self.mus.reset((self.S - 1, self.HFEngine.npar)) self.muTest = emptyParameterList() self.muTest.reset((len(muTestBasePivot) + 1, self.HFEngine.npar)) for k in range(self.S - 1): self.mus.data[k, self.directionPivot] = musPivot[k].data self.mus.data[k, self.directionMarginal] = self.musMargLoc[-1].data for k in range(len(muTestBasePivot)): self.muTest.data[k, self.directionPivot] = muTestBasePivot[k].data self.muTest.data[k, self.directionMarginal] = ( self.musMargLoc[-1].data) self.muTest.data[-1, self.directionPivot] = musPivot[-1].data self.muTest.data[-1, self.directionMarginal] = self.musMargLoc[-1].data if len(self.mus) > 0: vbMng(self, "MAIN", ("Adding first {} sample point{} at {} to training " "set.").format(self.S - 1, "" + "s" * (self.S > 2), self.mus), 2) self.samplingEngine.iterSample(self.mus) self._S = len(self.mus) self._approxParameters["S"] = self.S def setupApproxPivoted(self, mus:paramList) -> int: if self.checkComputedApproxPivoted(): return -1 if not hasattr(self, "_plotEstPivot"): self._plotEstPivot = "NONE" RROMPyAssert(self._mode, message = "Cannot setup approximant.") vbMng(self, "INIT", "Setting up {}.". format(self.name()), 10) self.computeScaleFactor() if self.trainedModel is None: self.trainedModel = self.tModelType() self.trainedModel.verbosity = self.verbosity self.trainedModel.timestamp = self.timestamp datadict = {"mu0": self.mu0, "projMat": np.zeros((0, 0)), "scaleFactor": self.scaleFactor, "rescalingExp": self.HFEngine.rescalingExp, "directionPivot": self.directionPivot} self.trainedModel.data = self.initializeModelData(datadict)[0] self.trainedModel.data.Qs, self.trainedModel.data.Ps = [], [] _trainedModelOld = copy(self.trainedModel) self._scaleFactorOldPivot = copy(self.scaleFactor) self.scaleFactor = self.scaleFactorPivot self._temporaryPivot = 1 self._samplingEngineOld = copy(self.samplingEngine) self.musMargLoc, self.samplingEngs = [], [None] * len(mus) Qs, Ps = [None] * len(mus), [None] * len(mus) self.verbosity -= 15 S0 = copy(self.S) for j, mu in enumerate(mus): RationalInterpolantGreedy.setupSampling(self) self.trainedModel = None self.musMargLoc += [mu] RationalInterpolantGreedy.setupApprox(self, self._plotEstPivot) self.samplingEngs[j] = copy(self.samplingEngine) Qs[j] = copy(self.trainedModel.data.Q) Ps[j] = copy(self.trainedModel.data.P) self._S = S0 self.scaleFactor = self._scaleFactorOldPivot del self._scaleFactorOldPivot, self._temporaryPivot self._finalizeSnapshots() del self._samplingEngineOld, self.musMargLoc, self.samplingEngs self._mus = self.samplingEngine.musCoalesced self.trainedModel = _trainedModelOld self.trainedModel.data.mus = copy(self.mus) self.trainedModel.data.musMarginal = copy(self.musMarginal) padRight = (self.samplingEngine.nsamplesTot - self.trainedModel.data.projMat.shape[1]) nmusOld = len(self.trainedModel.data.Ps) for j in range(nmusOld): nsj = self.samplingEngine.nsamples[j] self.trainedModel.data.Ps[j].pad(0, padRight) self.trainedModel.data.HIs[j].pad(0, padRight) padLeft = self.trainedModel.data.projMat.shape[1] for j in range(len(mus)): nsj = self.samplingEngine.nsamples[nmusOld + j] if self.POD == PODGlobal: rRightj = self.samplingEngine.RPODCPart[:, padLeft : padLeft + nsj] Ps[j].postmultiplyTensorize(rRightj.T) else: padRight -= nsj Ps[j].pad(padLeft, padRight) padLeft += nsj pMat = self.samplingEngine.samplesCoalesced.data pMatEff = dot(self.HFEngine.C, pMat) if self.approx_state else pMat self.trainedModel.data.projMat = pMatEff self.trainedModel.data.Qs += Qs self.trainedModel.data.Ps += Ps self.trainedModel.data.approxParameters = copy(self.approxParameters) self.verbosity += 15 vbMng(self, "DEL", "Done setting up approximant.", 10) return 0 def setupApprox(self, plotEst : str = "NONE") -> int: if self.checkComputedApprox(): return -1 if '_' not in plotEst: plotEst = plotEst + "_NONE" plotEstM, self._plotEstPivot = plotEst.split("_") val = super().setupApprox(plotEstM) return val \ No newline at end of file diff --git a/rrompy/reduction_methods/pivoted/greedy/rational_interpolant_pivoted_greedy.py b/rrompy/reduction_methods/pivoted/greedy/rational_interpolant_pivoted_greedy.py index 7920507..e9333e9 100644 --- a/rrompy/reduction_methods/pivoted/greedy/rational_interpolant_pivoted_greedy.py +++ b/rrompy/reduction_methods/pivoted/greedy/rational_interpolant_pivoted_greedy.py @@ -1,314 +1,318 @@ # 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, 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 np.inf; - 'S': total number of pivot samples current approximant relies upon; - 'samplerPivot': pivot sample point generator; - 'SMarginal': number of starting marginal samples; - 'samplerMarginalGrid': marginal sample point generator via sparse grid; - 'polybasis': type of polynomial basis for pivot interpolation; defaults to 'MONOMIAL'; - 'polybasisMarginal': type of polynomial basis for marginal interpolation; allowed values include 'MONOMIAL', 'CHEBYSHEV' and 'LEGENDRE'; defaults to 'MONOMIAL'; - '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. 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) 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 ad7c1cf..f32922d 100644 --- a/rrompy/reduction_methods/pivoted/rational_interpolant_greedy_pivoted.py +++ b/rrompy/reduction_methods/pivoted/rational_interpolant_greedy_pivoted.py @@ -1,427 +1,435 @@ # 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 totalDegreeN, dot from rrompy.utilities.poly_fitting.polynomial import polyvander as pv from rrompy.utilities.exception_manager import RROMPyAssert from rrompy.parameter import emptyParameterList, checkParameterList __all__ = ['RationalInterpolantGreedyPivoted'] class RationalInterpolantGreedyPivoted(GenericPivotedApproximant, RationalInterpolantGreedy): """ ROM pivoted rational interpolant (with pole matching) computation for parametric problems. Args: HFEngine: HF problem solver. mu0(optional): Default parameter. Defaults to 0. directionPivot(optional): Pivot components. Defaults to [0]. approxParameters(optional): Dictionary containing values for main parameters of approximant. Recognized keys are: - 'POD': whether to compute POD of snapshots; defaults to True; - 'matchingWeight': weight for pole matching optimization; defaults to 1; - 'cutOffTolerance': tolerance for ignoring parasitic poles; defaults to np.inf; + - '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. 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, 0, np.nan self._marginalizeTrainedModel(True) errRes = super().errorEstimator(mus, return_max) self._marginalizeTrainedModel(False) return errRes def _preliminaryTraining(self): """Initialize starting snapshots of solution map.""" RROMPyAssert(self._mode, message = "Cannot start greedy algorithm.") S = self.S self.sampleBatchIdx, self.sampleBatchSize, self._S = -1, 0, 0 nextBatchSize = 1 while self._S + nextBatchSize <= S: self.sampleBatchIdx += 1 self.sampleBatchSize = nextBatchSize self._S += self.sampleBatchSize nextBatchSize = totalDegreeN(self.npar - 1, self.sampleBatchIdx + 1) self.resetSamples() musPivot = self.trainSetGenerator.generatePoints(self.S) musPivot.data = musPivot.data[: self.S] muTestPivot = self.samplerPivot.generatePoints(self.nTestPoints) idxPop = pruneSamples(muTestPivot ** self.HFEngine.rescalingExp[ self.directionPivot[0]], musPivot ** self.HFEngine.rescalingExp[ self.directionPivot[0]], 1e-10 * self.scaleFactor[0]) self.mus = emptyParameterList() self.mus.reset((self.S, self.npar + len(self.musMargLoc))) muTestBase = emptyParameterList() muTestBase.reset((len(muTestPivot), self.npar + len(self.musMargLoc))) for k in range(self.S): self.mus.data[k, self.directionPivot] = musPivot[k].data self.mus.data[k, self.directionMarginal] = self.musMargLoc.data for k in range(len(muTestPivot)): muTestBase.data[k, self.directionPivot] = muTestPivot[k].data muTestBase.data[k, self.directionMarginal] = self.musMargLoc.data muTestBase.pop(idxPop) muTestBase = muTestBase.sort() muLast = copy(self.mus[-1]) self.mus.pop() if len(self.mus) > 0: vbMng(self, "MAIN", ("Adding first {} sample point{} at {} to training " "set.").format(self.S - 1, "" + "s" * (self.S > 2), self.mus), 2) self.samplingEngine.iterSample(self.mus) self._S = len(self.mus) self._approxParameters["S"] = self.S self.muTest = emptyParameterList() self.muTest.reset((len(muTestBase) + 1, self.mus.shape[1])) self.muTest.data[: -1] = muTestBase.data self.muTest.data[-1] = muLast.data def _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) 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.marginalInterp = self._setupMarginalInterp() self.trainedModel.data.Qs, self.trainedModel.data.Ps = Qs, Ps + vbMng(self, "INIT", "Matching poles.", 10) + self.trainedModel.initializeFromRational( + self.HFEngine, self.matchingWeight, + self.POD == 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 45ed3ef..19454b0 100644 --- a/rrompy/reduction_methods/pivoted/rational_interpolant_pivoted.py +++ b/rrompy/reduction_methods/pivoted/rational_interpolant_pivoted.py @@ -1,350 +1,358 @@ # 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, 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. 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) self.musMarginal = self.samplerMarginal.generatePoints( self.SMarginal) 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.marginalInterp = self._setupMarginalInterp() self.trainedModel.data.Qs, self.trainedModel.data.Ps = Qs, Ps + vbMng(self, "INIT", "Matching poles.", 10) + self.trainedModel.initializeFromRational( + self.HFEngine, self.matchingWeight, + self.POD == 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/trained_model/trained_model_pivoted.py b/rrompy/reduction_methods/pivoted/trained_model/trained_model_pivoted.py index f485d7f..d0bd3ce 100644 --- a/rrompy/reduction_methods/pivoted/trained_model/trained_model_pivoted.py +++ b/rrompy/reduction_methods/pivoted/trained_model/trained_model_pivoted.py @@ -1,385 +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 scipy.special import factorial as fact +from copy import deepcopy as copy from itertools import combinations from rrompy.reduction_methods.standard.trained_model import ( TrainedModelRational) -from rrompy.utilities.base.types import (Np1D, Tuple, List, ListAny, paramVal, - paramList, sampList, HFEng) +from rrompy.utilities.base.types import (Np1D, Np2D, Tuple, List, ListAny, + paramVal, paramList, sampList, HFEng) from rrompy.utilities.base import verbosityManager as vbMng, freepar as fp from rrompy.utilities.numerical import (pointMatching, chordalMetricAdjusted, - cutOffDistance) + cutOffDistance, reduceDegreeN) +from rrompy.utilities.poly_fitting.polynomial import (polybases as ppb, + PolynomialInterpolator as PI) +from rrompy.utilities.poly_fitting.radial_basis import (polybases as rbpb, + RadialBasisInterpolator as RBI) +from rrompy.utilities.poly_fitting.moving_least_squares import ( + MovingLeastSquaresInterpolator as MLSI) from rrompy.utilities.poly_fitting.heaviside import (heavisideUniformShape, rational2heaviside, HeavisideInterpolator as HI) -from rrompy.utilities.exception_manager import RROMPyException, RROMPyAssert +from rrompy.utilities.exception_manager import (RROMPyException, RROMPyAssert, + RROMPyWarning) from rrompy.parameter import checkParameter, checkParameterList -from rrompy.sampling import emptySampleList, sampleList +from rrompy.sampling import emptySampleList __all__ = ['TrainedModelPivoted'] class TrainedModelPivoted(TrainedModelRational): """ ROM approximant evaluation for pivoted approximants (with pole matching). Attributes: Data: dictionary with all that can be pickled. """ def centerNormalizePivot(self, mu : paramList = [], mu0 : paramVal = None) -> paramList: """ Compute normalized parameter to be plugged into approximant. Args: mu: Parameter(s) 1. mu0: Parameter(s) 2. If None, set to self.data.mu0Pivot. Returns: Normalized parameter. """ mu = checkParameterList(mu, self.data.nparPivot)[0] if mu0 is None: mu0 = self.data.mu0Pivot rad = ((mu ** self.data.rescalingExpPivot - mu0 ** self.data.rescalingExpPivot) / self.data.scaleFactorPivot) return rad def centerNormalizeMarginal(self, mu : paramList = [], mu0 : paramVal = None) -> paramList: """ Compute normalized parameter to be plugged into approximant. Args: mu: Parameter(s) 1. mu0: Parameter(s) 2. If None, set to self.data.mu0Marginal. Returns: Normalized parameter. """ mu = checkParameterList(mu, self.data.nparMarginal)[0] if mu0 is None: mu0 = self.data.mu0Marginal rad = ((mu ** self.data.rescalingExpMarginal - mu0 ** self.data.rescalingExpMarginal) / self.data.scaleFactorMarginal) return rad def updateEffectiveSamples(self, HFEngine:HFEng, exclude : List[int] = None, matchingWeight : float = 1., POD : bool = True, is_state : bool = True): if hasattr(self, "_idxExcl"): for j, excl in enumerate(self._idxExcl): self.data.musMarginal.insert(self._musMExcl[j], excl) self.data.HIs.insert(excl, self._HIsExcl[j]) self.data.Ps.insert(excl, self._PsExcl[j]) self.data.Qs.insert(excl, self._QsExcl[j]) if exclude is None: exclude = [] self._idxExcl, self._musMExcl = list(np.sort(exclude)), [] self._HIsExcl, self._PsExcl, self._QsExcl = [], [], [] for excl in self._idxExcl[::-1]: self._musMExcl = [self.data.musMarginal[excl]] + self._musMExcl self.data.musMarginal.pop(excl) self._HIsExcl = [self.data.HIs.pop(excl)] + self._HIsExcl self._PsExcl = [self.data.Ps.pop(excl)] + self._PsExcl self._QsExcl = [self.data.Qs.pop(excl)] + self._QsExcl poles = [hi.poles for hi in self.data.HIs] coeffs = [hi.coeffs for hi in self.data.HIs] self.initializeFromLists(poles, coeffs, self.data.HIs[0].polybasis, HFEngine, matchingWeight, POD, is_state) + def setupMarginalInterp(self, approx, interpPars:ListAny, + MMAuto:bool, rDWM : Np1D = None, + noWarnReduceAuto : bool = True): + vbMng(self, "INIT", "Starting computation of marginal interpolator.", + 7) + musMCN = self.centerNormalizeMarginal(self.data.musMarginal) + pbM = approx.polybasisMarginal + if pbM not in ppb: + rDWMEf = np.array(rDWM) + self.data.marginalInterp = [] + for ipts, pts in enumerate(self.data.suppEffPts): + mI = [None] * len(musMCN) + if len(pts) > 0: + musMCNEff = musMCN[pts] + if MMAuto: + if ipts > 0: + verb = approx.verbosity + approx.verbosity = 0 + _musM = approx.musMarginal + approx.musMarginal = musMCNEff + approx._setMMarginalAuto() + if ipts > 0: + approx.musMarginal = _musM + approx.verbosity = verb + else: + MM = reduceDegreeN(approx.MMarginal, len(musMCNEff), + self.data.nparMarginal, + approx.polydegreetypeMarginal) + if MM < approx.MMarginal: + if ipts == 0 and not noWarnReduceAuto: + RROMPyWarning(("MMarginal too large compared to " + "SMarginal. Reducing MMarginal by " + "{}").format(approx.MMarginal - MM)) + approx.MMarginal = MM + MMEff = approx.MMarginal + for j in range(len(musMCNEff)): + canonicalj = 1. * (np.arange(len(musMCNEff)) == j) + while MMEff >= 0 and (pbM in ppb + or rDWMEf[0] <= rDWM[0] * 2 ** 6): + pParRest = copy(interpPars) + if pbM in ppb: + p = PI() + else: + pParRest = [rDWMEf] + pParRest + p = RBI() if pbM in rbpb else MLSI() + wellCond, msg = p.setupByInterpolation(musMCNEff, + canonicalj, MMEff, + pbM, *pParRest) + vbMng(self, "MAIN", msg, 5) + if wellCond: break + if pbM in ppb: + vbMng(self, "MAIN", + ("Polyfit is poorly conditioned. Reducing " + "MMarginal by 1."), 10) + MMEff -= 1 + else: + vbMng(self, "MAIN", + ("Polyfit is poorly conditioned. " + "Multiplying radialDirectionalWeightsMarginal " + "by 2."), 10) + rDWMEf *= 2. + if MMEff < 0 or (pbM not in ppb + and rDWMEf[0] > rDWM[0] * 2 ** 6): + raise RROMPyException(("Instability in computation of " + "interpolant. Aborting.")) + if pbM in ppb: MMEff = approx.MMarginal + else: rDWMEf = np.array(rDWM) + mI[pts[j]] = copy(p) + self.data.marginalInterp += [mI] + _MMarginalEffective = approx.MMarginal + approx.MMarginal = _MMarginalEffective + vbMng(self, "DEL", "Done computing marginal interpolator.", 7) + def initializeFromLists(self, poles:ListAny, coeffs:ListAny, basis:str, HFEngine:HFEng, matchingWeight : float = 1., POD : bool = True, is_state : bool = True): """Initialize Heaviside representation.""" musM = self.data.musMarginal margAbsDist = np.sum(np.abs(np.repeat(musM.data, len(musM), 0) - np.tile(musM.data, [len(musM), 1]) ), axis = 1).reshape(len(musM), len(musM)) explored = [0] unexplored = list(range(1, len(musM))) poles, coeffs = heavisideUniformShape(poles, coeffs) N = len(poles[0]) for _ in range(1, len(musM)): minIdx = np.argmin(np.concatenate([margAbsDist[ex, unexplored] \ for ex in explored])) minIex = explored[minIdx // len(unexplored)] minIunex = unexplored[minIdx % len(unexplored)] resex = coeffs[minIex][: N] resunex = coeffs[minIunex][: N] if matchingWeight != 0 and not POD: resex = self.data.projMat.dot(resex.T) resunex = self.data.projMat.dot(resunex.T) dist = chordalMetricAdjusted(poles[minIex], poles[minIunex], matchingWeight, resex, resunex, HFEngine, is_state) - reordering = pointMatching(dist) + reordering = pointMatching(dist)[1] poles[minIunex] = poles[minIunex][reordering] coeffs[minIunex][: N] = coeffs[minIunex][reordering] explored += [minIunex] unexplored.remove(minIunex) HIs = [] for pls, cfs in zip(poles, coeffs): hsi = HI() hsi.poles = pls hsi.coeffs = cfs hsi.npar = 1 hsi.polybasis = basis HIs += [hsi] self.data.HIs = HIs + self.data.suppEffPts = [np.arange(len(self.data.HIs))] + self.data.suppEffIdx = np.zeros(N, dtype = int) def initializeFromRational(self, HFEngine:HFEng, matchingWeight : float = 1., POD : bool = True, is_state : bool = True): """Initialize Heaviside representation.""" RROMPyAssert(self.data.nparPivot, 1, "Number of pivot parameters") poles, coeffs = [], [] for Q, P in zip(self.data.Qs, self.data.Ps): cfs, pls, basis = rational2heaviside(P, Q) poles += [pls] coeffs += [cfs] self.initializeFromLists(poles, coeffs, basis, HFEngine, matchingWeight, POD, is_state) - def recompressByCutOff(self, murange : Tuple[float, float] = [- 1., 1.], - tol : float = np.inf): + def recompressByCutOff(self, tol : float = np.inf, kind : str = "STANDARD", + murange : Tuple[float, float] = [- 1., 1.]) -> str: N = len(self.data.HIs[0].poles) mu0 = np.mean(murange) - krPoles = np.all([cutOffDistance(hi.poles, murange) <= tol - for hi in self.data.HIs], axis = 0) - keepPole = np.where(krPoles)[0] - removePole = np.where(np.logical_not(krPoles))[0] - if len(keepPole) == N: return " No poles erased." - keepCoeff = np.append(keepPole, np.append([N], - np.arange(N + 1, len(self.data.HIs[0].coeffs)))) - for hi in self.data.HIs: - polyCorrection = np.zeros_like(hi.coeffs[0, :]) - for j in removePole: - if not np.isinf(hi.poles[j]): - polyCorrection += hi.coeffs[j, :] / (mu0 - hi.poles[j]) - if len(hi.coeffs) == N: - hi.coeffs = np.vstack((hi.coeffs, polyCorrection)) - else: - hi.coeffs[N, :] += polyCorrection - hi.poles = hi.poles[keepPole] - hi.coeffs = hi.coeffs[keepCoeff, :] - return (" Erased {} pole".format(len(removePole)) - + "s" * (len(removePole) > 1) + ".") + goodLocPoles = np.array([cutOffDistance(hi.poles, murange) <= tol + for hi in self.data.HIs]) + self.data.suppEffPts = [np.arange(len(self.data.HIs))] + self.data.suppEffIdx = np.zeros(N, dtype = int) + if np.all(goodLocPoles): + return " No poles erased." + kind = kind.upper().strip().replace(" ","") + goodAllPoles = np.all(goodLocPoles, axis = 0) + badPoles = np.logical_not(goodAllPoles) + if kind == "HARD": + keepPole = np.where(goodAllPoles)[0] + halfPole = np.empty(0, dtype = int) + removePole = np.where(badPoles)[0] + elif kind == "SOFT": + goodSomePoles = np.any(goodLocPoles, axis = 0) + keepPole = np.where(goodSomePoles)[0] + halfPole = np.where(np.logical_and(badPoles, goodSomePoles))[0] + removePole = np.where(np.logical_not(goodSomePoles))[0] + else: + raise RROMPyException("Cutoff kind not recognized.") + if len(removePole) > 0: + keepCoeff = np.append(keepPole, np.append([N], + np.arange(N + 1, len(self.data.HIs[0].coeffs)))) + for hi in self.data.HIs: + polyCorrection = np.zeros_like(hi.coeffs[0, :]) + for j in removePole: + if not np.isinf(hi.poles[j]): + polyCorrection += hi.coeffs[j, :] / (mu0 - hi.poles[j]) + if len(hi.coeffs) == N: + hi.coeffs = np.vstack((hi.coeffs, polyCorrection)) + else: + hi.coeffs[N, :] += polyCorrection + hi.poles = hi.poles[keepPole] + hi.coeffs = hi.coeffs[keepCoeff, :] + for idxR in halfPole: + pts = np.where(goodLocPoles[:, idxR])[0] + idxEff = len(self.data.suppEffPts) + for idEff, prevPts in enumerate(self.data.suppEffPts): + if len(prevPts) == len(pts): + if np.allclose(prevPts, pts): + idxEff = idEff + break + if idxEff == len(self.data.suppEffPts): + self.data.suppEffPts += [pts] + self.data.suppEffIdx[idxR] = idxEff + self.data.suppEffIdx = self.data.suppEffIdx[keepPole] + return (" Hard-erased {} pole".format(len(removePole)) + + "s" * (len(removePole) != 1) + + " and soft-erased {} pole".format(len(halfPole)) + + "s" * (len(halfPole) != 1) + ".") - def interpolateMarginal(self, mu : paramList = [], - samples : ListAny = [], der : List[int] = None, - scl : Np1D = None) -> sampList: - mu = checkParameterList(mu, self.data.nparMarginal)[0] - sList = isinstance(samples[0], sampleList) - sEff = [None] * len(samples) - for j in range(len(samples)): - if sList: - sEff[j] = samples[j].data - else: - sEff[j] = samples[j] - try: - dtype = sEff[0].dtype - except: - dtype = sEff[0][0].dtype - vbMng(self, "INIT", "Interpolating marginal at mu = {}.".format(mu), - 95) - muC = self.centerNormalizeMarginal(mu) - p = emptySampleList() - p.reset((len(sEff[0]), len(muC)), dtype = dtype) - p.data[:, :] = 0. - if len(sEff[0]) > 0: - for mIj, spj in zip(self.data.marginalInterp, sEff): - p = p + spj.reshape(len(sEff[0]), - 1) * mIj(muC, der, scl) - vbMng(self, "DEL", "Done interpolating marginal.", 95) - if not sList: - p = p.data.flatten() - return p + def _interpolateMarginal(self, muC : paramList, objs : ListAny) -> Np2D: + res = np.zeros(objs[0].shape + (len(muC),), dtype = objs[0].dtype) + for suppIdx in range(len(self.data.suppEffPts)): + i = np.where(self.data.suppEffIdx == suppIdx)[0] + if suppIdx == 0: + i = np.append(i, np.arange(len(self.data.suppEffIdx), + len(res))) + if len(i) > 0: + for mIj, obj in zip(self.data.marginalInterp[suppIdx], objs): + if mIj is not None: + res[i] += np.expand_dims(obj[i], - 1) * mIj(muC) + return res def interpolateMarginalInterpolator(self, mu : paramVal = []) -> Np1D: """Obtain interpolated approximant interpolator.""" mu = checkParameter(mu, self.data.nparMarginal)[0] hsi = HI() - hsi.poles = self.interpolateMarginalPoles(mu) - hsi.coeffs = self.interpolateMarginalCoeffs(mu) + hsi.poles = self.interpolateMarginalPoles(mu)[..., 0] + hsi.coeffs = self.interpolateMarginalCoeffs(mu)[..., 0] hsi.npar = 1 hsi.polybasis = self.data.HIs[0].polybasis return hsi def interpolateMarginalPoles(self, mu : paramList = []) -> Np1D: """Obtain interpolated approximant poles.""" mu = checkParameterList(mu, self.data.nparMarginal)[0] - polesInf = np.any([np.isinf(hi.poles) for hi in self.data.HIs], - axis = 0) - intMPoles = self.interpolateMarginal(mu, [hi.poles - for hi in self.data.HIs]) - intMPoles[polesInf] = np.inf + muC = self.centerNormalizeMarginal(mu) + vbMng(self, "INIT", + "Interpolating marginal poles at mu = {}.".format(mu), 95) + intMPoles = self._interpolateMarginal(muC, + [hi.poles for hi in self.data.HIs]) + vbMng(self, "DEL", "Done interpolating marginal poles.", 95) return intMPoles def interpolateMarginalCoeffs(self, mu : paramList = []) -> Np1D: """Obtain interpolated approximant coefficients.""" mu = checkParameterList(mu, self.data.nparMarginal)[0] - cs = self.interpolateMarginal(mu, [hi.coeffs for hi in self.data.HIs]) - if isinstance(cs, (list, tuple,)): cs = np.array(cs) - return cs.reshape(self.data.HIs[0].coeffs.shape) + muC = self.centerNormalizeMarginal(mu) + vbMng(self, "INIT", + "Interpolating marginal coefficients at mu = {}.".format(mu), 95) + intMCoeffs = self._interpolateMarginal(muC, + [hi.coeffs for hi in self.data.HIs]) + vbMng(self, "DEL", "Done interpolating marginal coefficients.", 95) + return intMCoeffs def getApproxReduced(self, mu : paramList = []) -> sampList: """ Evaluate reduced representation of approximant at arbitrary parameter. Args: mu: Target parameter. """ RROMPyAssert(self.data.nparPivot, 1, "Number of pivot parameters") mu = checkParameterList(mu, self.data.npar)[0] if (not hasattr(self, "lastSolvedApproxReduced") or self.lastSolvedApproxReduced != mu): vbMng(self, "INIT", "Evaluating approximant at mu = {}.".format(mu), 12) self.uApproxReduced = emptySampleList() for i, muPL in enumerate(mu): muL = self.centerNormalizePivot([muPL(0, x) \ for x in self.data.directionPivot]) muM = [muPL(0, x) for x in self.data.directionMarginal] vbMng(self, "INIT", "Assembling reduced model for mu = {}.".format(muPL), 87) hsL = self.interpolateMarginalInterpolator(muM) vbMng(self, "DEL", "Done assembling reduced model.", 87) uAppR = hsL(muL) if i == 0: - #self.data.HIs[0].coeffs.shape[1], len(mu) self.uApproxReduced.reset((len(uAppR), len(mu)), dtype = uAppR.dtype) self.uApproxReduced[i] = uAppR vbMng(self, "DEL", "Done evaluating approximant.", 12) self.lastSolvedApproxReduced = mu return self.uApproxReduced def getPVal(self, mu : paramList = []) -> sampList: """ Evaluate rational numerator at arbitrary parameter. Args: mu: Target parameter. """ RROMPyAssert(self.data.nparPivot, 1, "Number of pivot parameters") mu = checkParameterList(mu, self.data.npar)[0] p = emptySampleList() p.reset((len(self.data.HIs[0].coeffs.shape[1]), len(mu))) for i, muPL in enumerate(mu): muL = self.centerNormalizePivot([muPL(0, x) \ for x in self.data.directionPivot]) muM = [muPL(0, x) for x in self.data.directionMarginal] hsL = self.interpolateMarginalInterpolator(muM) p[i] = hsL(muL) * np.prod(muL(0, 0) - hsL.poles) return p def getQVal(self, mu:Np1D, der : List[int] = None, scl : Np1D = None) -> Np1D: """ Evaluate rational denominator at arbitrary parameter. Args: mu: Target parameter. der(optional): Derivatives to take before evaluation. """ RROMPyAssert(self.data.nparPivot, 1, "Number of pivot parameters") mu = checkParameterList(mu, self.data.npar)[0] muP = self.centerNormalizePivot(checkParameterList( mu.data[:, self.data.directionPivot], self.data.nparPivot)[0]) muM = checkParameterList(mu.data[:, self.data.directionMarginal], self.data.nparMarginal)[0] if der is None: derP, derM = 0, [0] else: derP = der[self.data.directionPivot[0]] derM = [der[x] for x in self.data.directionMarginal] if np.any(np.array(derM) != 0): raise RROMPyException(("Derivatives of Q with respect to marginal " "parameters not allowed.")) sclP = 1 if scl is None else scl[self.data.directionPivot[0]] derVal = np.zeros(len(mu), dtype = np.complex) N = len(self.data.HIs[0].poles) if derP == N: derVal[:] = 1. elif derP >= 0 and derP < N: - pls = self.interpolateMarginalPoles(muM).reshape(-1, len(mu)).T + pls = self.interpolateMarginalPoles(muM).T plsDist = muP.data.reshape(-1, 1) - pls for terms in combinations(np.arange(N), N - derP): derVal += np.prod(plsDist[:, list(terms)], axis = 1) return sclP ** derP * fact(derP) * derVal def getPoles(self, *args, **kwargs) -> Np1D: """ Obtain approximant poles. Returns: Numpy complex vector of poles. """ RROMPyAssert(self.data.nparPivot, 1, "Number of pivot parameters") if len(args) + len(kwargs) > 1: raise RROMPyException(("Wrong number of parameters passed. " "Only 1 available.")) elif len(args) + len(kwargs) == 1: if len(args) == 1: mVals = args[0] else: mVals = kwargs["marginalVals"] if not hasattr(mVals, "__len__"): mVals = [mVals] mVals = list(mVals) else: mVals = [fp] try: rDim = mVals.index(fp) if rDim < len(mVals) - 1 and fp in mVals[rDim + 1 :]: raise except: raise RROMPyException(("Exactly 1 'freepar' entry in " "marginalVals must be provided.")) if rDim != self.data.directionPivot[0]: raise RROMPyException(("'freepar' entry in marginalVals must " "coincide with pivot direction.")) mVals[rDim] = self.data.mu0(rDim) mMarg = [mVals[j] for j in range(len(mVals)) if j != rDim] - roots = np.array(self.interpolateMarginalPoles(mMarg)) + roots = self.interpolateMarginalPoles(mMarg)[..., 0] return np.power(self.data.mu0(rDim) ** self.data.rescalingExp[rDim] + self.data.scaleFactor[rDim] * roots, 1. / self.data.rescalingExp[rDim]) def getResidues(self, *args, **kwargs) -> Np1D: """ Obtain approximant residues. Returns: Numpy matrix with residues as columns. """ pls = self.getPoles(*args, **kwargs) if len(args) == 1: mVals = args[0] elif len(args) == 0: mVals = [None] else: mVals = kwargs["marginalVals"] if not hasattr(mVals, "__len__"): mVals = [mVals] mVals = list(mVals) rDim = mVals.index(fp) mMarg = [mVals[j] for j in range(len(mVals)) if j != rDim] - residues = self.interpolateMarginalCoeffs(mMarg)[: len(pls)] + residues = self.interpolateMarginalCoeffs(mMarg)[: len(pls), :, 0] res = self.data.projMat.dot(residues.T) return pls, res diff --git a/rrompy/utilities/numerical/point_matching.py b/rrompy/utilities/numerical/point_matching.py index ba6dc9f..af2db06 100644 --- a/rrompy/utilities/numerical/point_matching.py +++ b/rrompy/utilities/numerical/point_matching.py @@ -1,76 +1,76 @@ # Copyright (C) 2018 by the RROMPy authors # # This file is part of RROMPy. # # RROMPy is free software: you can redistribute it and/or modify # it under the terms of the GNU Lesser General Public License as published by # the Free Software Foundation, either version 3 of the License, or # (at your option) any later version. # # RROMPy is distributed in the hope that it will be useful, # but WITHOUT ANY WARRANTY; without even the implied warranty of # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the # GNU Lesser General Public License for more details. # # You should have received a copy of the GNU Lesser General Public License # along with RROMPy. If not, see . # import numpy as np from scipy.optimize import linear_sum_assignment as LSA from rrompy.utilities.base.types import Tuple, Np1D, Np2D, HFEng __all__ = ['pointMatching', 'cutOffDistance', 'angleTable', 'chordalMetricTable', 'chordalMetricAdjusted'] -def pointMatching(distanceMatrix:Np2D) -> Np1D: - return LSA(distanceMatrix)[1] +def pointMatching(distanceMatrix:Np2D) -> Tuple[Np1D, Np1D]: + return LSA(distanceMatrix) def cutOffDistance(x:Np1D, murange : Tuple[float, float] = [- 1., 1.]) -> Np1D: mu0 = np.mean(murange) musig = murange[0] - mu0 isInf = np.isinf(x) dist = np.empty(len(x)) dist[isInf] = np.inf xEffR = x[np.logical_not(isInf)] - mu0 if np.isclose(musig, 0.): dist[np.logical_not(isInf)] = np.abs(xEffR) else: xEffR /= musig bernEff = (xEffR ** 2. - 1) ** .5 dist[np.logical_not(isInf)] = np.max(np.vstack(( np.abs(xEffR + bernEff), np.abs(xEffR - bernEff) )), axis = 0) - 1. return dist def angleTable(X:Np2D, Y:Np2D, HFEngine : HFEng = None, is_state : bool = True) -> Np2D: if HFEngine is None: innerT = Y.dot(X.T.conj()) normX = np.linalg.norm(X, axis = 1) normY = np.linalg.norm(Y, axis = 1) else: innerT = HFEngine.innerProduct(X, Y, is_state = is_state) normX = HFEngine.norm(X, is_state = is_state) normY = HFEngine.norm(Y, is_state = is_state) xInf = np.where(np.isclose(normX, 0.))[0] normX[xInf] = 1. return - (np.abs(innerT / normX).T - normY) def chordalMetricTable(x:Np1D, y:Np1D) -> Np2D: x, y = np.array(x), np.array(y) xInf, yInf = np.where(np.isinf(x))[0], np.where(np.isinf(y))[0] x[xInf], y[yInf] = 0., 0. T = np.abs(np.tile(x.reshape(-1, 1), len(y)) - y.reshape(1, -1)) T[xInf, :], T[:, yInf] = 1., 1. T[np.ix_(xInf, yInf)] = 0. T = (T.T * (np.abs(x) ** 2. + 1.) ** -.5).T * (np.abs(y) ** 2. + 1.) ** -.5 return T def chordalMetricAdjusted(x:Np1D, y:Np1D, w : float = 0, X : Np2D = None, Y : Np2D = None, HFEngine : HFEng = None, is_state : bool = True) -> Np2D: dist = chordalMetricTable(x, y) if w == 0: return dist distAdj = angleTable(X, Y, HFEngine, is_state) return dist + w * distAdj