diff --git a/examples/5_anisotropic_square/anisotropic_square.py b/examples/5_anisotropic_square/anisotropic_square.py index 038a654..949d536 100644 --- a/examples/5_anisotropic_square/anisotropic_square.py +++ b/examples/5_anisotropic_square/anisotropic_square.py @@ -1,81 +1,82 @@ ### example from Smetana, Zahm, Patera. Randomized residual-based error ### estimators for parametrized equations. import numpy as np import matplotlib.pyplot as plt from itertools import product from anisotropic_square_engine import (AnisotropicSquareEngine as engine, AnisotropicSquareEnginePoles as plsEx) from rrompy.reduction_methods import ( RationalInterpolantGreedyPivotedGreedyPoleMatch as RIGPG) from rrompy.parameter.parameter_sampling import (QuadratureSampler as QS, SparseGridSampler as SGS) zs, Ls = [10., 50.], [.2, 1.2] z0, L0, n = np.mean(zs), np.mean(Ls), 50 murange = [[zs[0], Ls[0]], [zs[-1], Ls[-1]]] np.random.seed(4020) mu = [zs[0] + np.random.rand() * (zs[-1] - zs[0]), Ls[0] + np.random.rand() * (Ls[-1] - Ls[0])] solver = engine(z0, L0, n) fighandles = [] params = {"POD": True, "nTestPoints": 100, "greedyTol": 1e-4, "S": 3, "polybasisMarginal": "PIECEWISE_LINEAR_UNIFORM", "polybasis": "LEGENDRE", "samplerPivot":QS(zs, "UNIFORM"), "samplerTrainSet":QS(zs, "UNIFORM"), "errorEstimatorKind":"LOOK_AHEAD_RES", "errorEstimatorKindMarginal":"LOOK_AHEAD_RECOVER", "matchingChordalRadius": [1., "AUTO"], "SMarginal": 3, "paramsMarginal": {"MMarginal": 2, "radialDirectionalWeightsMarginalAdapt": [1e9, 1e12]}, "greedyTolMarginal": 1e-2, "samplerMarginal":SGS(Ls), - "radialDirectionalWeightsMarginal": [4.], "matchingWeight": 1.} + "radialDirectionalWeightsMarginal": [4.], "matchingWeight": 1., + "badPoleCorrection": "POLYNOMIAL"} for shared, tol in product([1., 0.], [1., 3.]): print("Testing cutoff tolerance {} with shared ratio {}.".format(tol, shared)) solver.cutOffPolesRMinRel = - 1. - tol solver.cutOffPolesRMaxRel = 1. + tol params["matchingShared"] = shared approx = RIGPG([0], solver, mu0 = [z0, L0], approxParameters = params, verbosity = 5) approx.setupApprox("ALL") verb = approx.verbosity approx.verbosity = 0 tspace = np.linspace(Ls[0], Ls[-1], 100) for j, t in enumerate(tspace): plsE = plsEx(t, 0., zs[-1]) pls = approx.getPoles([None, t]) pls[np.abs(np.imag(pls)) > 1e-5] = np.nan if j == 0: polesE = np.empty((len(tspace), len(plsE))) poles = np.empty((len(tspace), len(pls))) polesE[:] = np.nan if len(plsE) > polesE.shape[1]: nanR = np.empty((len(tspace), len(plsE) - polesE.shape[1])) nanR[:] = np.nan polesE = np.hstack((polesE, nanR)) polesE[j, : len(plsE)] = np.real(plsE) poles[j] = np.real(pls) approx.verbosity = verb fighandles += [plt.figure(figsize = (17, 5))] ax1 = fighandles[-1].add_subplot(1, 2, 1) ax2 = fighandles[-1].add_subplot(1, 2, 2) ax1.plot(poles, tspace) ax1.set_ylim(Ls) ax1.set_xlabel("mu_1") ax1.set_ylabel("mu_2") ax1.grid() ax2.plot(polesE, tspace, "k-.", linewidth = 1) ax2.plot(poles, tspace) for mm in approx.musMarginal: ax2.plot(zs, [mm[0, 0]] * 2, "k--", linewidth = 1) ax2.set_xlim(zs) ax2.set_ylim(Ls) ax2.set_xlabel("mu_1") ax2.set_ylabel("mu_2") ax2.grid() plt.show() print("\n") diff --git a/rrompy/parameter/parameter_sampling/sparse_grid/sparse_grid_sampler.py b/rrompy/parameter/parameter_sampling/sparse_grid/sparse_grid_sampler.py index 9580348..72861d4 100644 --- a/rrompy/parameter/parameter_sampling/sparse_grid/sparse_grid_sampler.py +++ b/rrompy/parameter/parameter_sampling/sparse_grid/sparse_grid_sampler.py @@ -1,108 +1,108 @@ # Copyright (C) 2018-2020 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 from itertools import product import numpy as np from rrompy.parameter.parameter_sampling.generic_sampler import GenericSampler from rrompy.utilities.poly_fitting.piecewise_linear import (sparsekinds, sparseMap) from rrompy.utilities.base.types import Tuple, List, Np1D, DictAny, paramList from rrompy.utilities.exception_manager import RROMPyException __all__ = ['SparseGridSampler'] class SparseGridSampler(GenericSampler): """Generator of sparse grid sample points.""" def __init__(self, lims:paramList, kind : str = "UNIFORM", parameterMap : DictAny = 1.): super().__init__(lims = lims, parameterMap = parameterMap) self.kind = kind self.reset() def __str__(self) -> str: return "{}[{}_{}]_{}".format(self.name(), self.lims[0], self.lims[1], self.kind) @property def npoints(self): """Number of points.""" return len(self.points) @property def kind(self): """Value of kind.""" return self._kind @kind.setter def kind(self, kind): if kind.upper() not in [sk.split("_")[2] + extra for sk, extra in product(sparsekinds, ["", "-HAAR"])]: raise RROMPyException("Generator kind not recognized.") self._kind = kind.upper() self._noBoundary = "HAAR" in self._kind def reset(self): limsE = self.mapParameterList(self.lims) centerEff = .5 * (limsE[0] + limsE[1]) self.points = self.mapParameterList(centerEff, "B") self.depth = np.array([[self._noBoundary] * self.npar], dtype = int) def refine(self, active : List[int] = None) -> Tuple[List[int], List[int]]: if active is None: active = np.arange(self.npoints) active = np.array(active) if np.any(active < 0) or np.any(active >= self.npoints): raise RROMPyException(("Active indices must be between 0 " "(included) and npoints (excluded).")) newIdxs, oldIdxs = [], [] for act in active: point, dpt = self.points[act], self.depth[act] for jdelta, signdelta in product(range(self.npar), [-1., 1.]): idx = self.addForwardPoint(point, dpt, jdelta, signdelta) if idx is not None: if idx > 0: newIdxs += [idx] - else: oldIdxs += [- idx] + elif - idx not in newIdxs: oldIdxs += [- idx] return newIdxs, oldIdxs def addForwardPoint(self, basepoint:Np1D, basedepth:Np1D, index:int, sign:float) -> int: if basedepth[index] < self._noBoundary: return None #makeshift skip for wrong boundary points at lvl 1 limd = self.mapParameterList(self.lims(index), idx = [index])(0) xd0 = sparseMap(self.mapParameterList(basepoint[index], idx = [index])(0, 0), limd, self.kind, False) + .5 ** basedepth[index] * sign if np.abs(xd0) >= 1. + 1e-15 * (1 - 2 * self._noBoundary): return None #point out of bounds pt = copy(basepoint) pt[index] = self.mapParameterList(sparseMap(xd0, limd, self.kind), "B", [index])(0, 0) dist = np.sum(np.abs(self.points.data - pt.reshape(1, -1)), axis = 1) samePt = np.where(np.isclose(dist, 0., atol = 1e-15))[0] if len(samePt) > 0: #point already exists return - samePt[0] self.points.append(pt) self.depth = np.append(self.depth, [basedepth], 0) self.depth[-1, index] += 1 return self.npoints - 1 def generatePoints(self, n:int, reorder = None) -> paramList: if self.npoints > n: self.reset() idx = np.arange(self.npoints) while self.npoints < n: idx = self.refine(idx)[0] return self.points diff --git a/rrompy/parameter/parameter_sampling/sparse_grid/sparse_grid_sampler_two_way.py b/rrompy/parameter/parameter_sampling/sparse_grid/sparse_grid_sampler_two_way.py index a18b359..7c0ec01 100644 --- a/rrompy/parameter/parameter_sampling/sparse_grid/sparse_grid_sampler_two_way.py +++ b/rrompy/parameter/parameter_sampling/sparse_grid/sparse_grid_sampler_two_way.py @@ -1,47 +1,47 @@ # Copyright (C) 2018-2020 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 itertools import product import numpy as np from .sparse_grid_sampler import SparseGridSampler from rrompy.utilities.base.types import Tuple, List from rrompy.utilities.exception_manager import RROMPyException __all__ = ['SparseGridSamplerTwoWay'] class SparseGridSamplerTwoWay(SparseGridSampler): """Generator of sparse grid sample points with two-way refinement.""" def refine(self, active : List[int] = None) -> Tuple[List[int], List[int]]: if active is None: active = np.arange(self.npoints) active = np.array(active) if np.any(active < 0) or np.any(active >= self.npoints): raise RROMPyException(("Active indices must be between 0 " "(included) and npoints (excluded).")) newIdxs, oldIdxs = [], [] for act in active: point, dpt = self.points[act], self.depth[act] for jdelta, signdelta, backwards in product(range(self.npar), [-1., 1.], [0, 1]): if backwards: dpt[jdelta] -= 1 idx = self.addForwardPoint(point, dpt, jdelta, signdelta) if idx is not None: if idx > 0: newIdxs += [idx] - else: oldIdxs += [- idx] + elif - idx not in newIdxs: oldIdxs += [- idx] if backwards: dpt[jdelta] += 1 return newIdxs, oldIdxs diff --git a/rrompy/reduction_methods/base/trained_model/trained_model.py b/rrompy/reduction_methods/base/trained_model/trained_model.py index eb5fbf9..3102e22 100644 --- a/rrompy/reduction_methods/base/trained_model/trained_model.py +++ b/rrompy/reduction_methods/base/trained_model/trained_model.py @@ -1,114 +1,115 @@ # Copyright (C) 2018-2020 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 abc import abstractmethod import numpy as np from rrompy.utilities.base.types import Np1D, List, paramList, sampList from rrompy.utilities.expression import expressionEvaluator +from rrompy.utilities.numerical import dot from rrompy.parameter import checkParameterList, emptyParameterList from rrompy.sampling import sampleList __all__ = ['TrainedModel'] class TrainedModel: """ ABSTRACT ROM approximant evaluation. Attributes: Data: dictionary with all that can be pickled. """ def name(self) -> str: return self.__class__.__name__ def __str__(self) -> str: return self.name() def __repr__(self) -> str: return self.__str__() + " at " + hex(id(self)) def reset(self): self.lastSolvedApproxReduced = None self.lastSolvedApprox = None def compress(self, collapse : bool = False, tol : float = 0.): if collapse: self.data.projMat = 1. self.data._collapsed = True if tol > 0.: self.data._compressTol = tol @property def npar(self): """Number of parameters.""" return self.data.mu0.shape[1] def checkParameterList(self, mu:paramList, check_if_single : bool = False) -> paramList: return checkParameterList(mu, self.data.npar, check_if_single) def mapParameterList(self, mu:paramList, direct : str = "F", idx : List[int] = None) -> paramList: if idx is None: idx = np.arange(self.npar) muMapped = checkParameterList(mu, len(idx)) for j, d in enumerate(idx): muMapped.data[:, j] = expressionEvaluator( self.data.parameterMap[direct][d], muMapped(j)).flatten() return muMapped @abstractmethod def getApproxReduced(self, mu : paramList = []) -> sampList: """ Evaluate reduced representation of approximant at arbitrary parameter. (ABSTRACT) Args: mu: Target parameter. """ pass def getApprox(self, mu : paramList = []) -> sampList: """ Evaluate approximant at arbitrary parameter. Args: mu: Target parameter. """ mu = self.checkParameterList(mu) if (not hasattr(self, "lastSolvedApprox") or self.lastSolvedApprox != mu): uApproxR = self.getApproxReduced(mu) if self.data._collapsed: self.uApprox = uApproxR else: for i, uApR in enumerate(uApproxR): uApREff = uApR - uAp = self.data.projMat.dot(uApREff) + uAp = dot(self.data.projMat, uApREff) if i == 0: uApprox = np.empty((len(uAp), len(uApproxR)), dtype = uAp.dtype) uApprox[:, i] = uAp self.uApprox = sampleList(uApprox) self.lastSolvedApprox = mu return self.uApprox def getPoles(self, *args, **kwargs) -> Np1D: """Obtain approximant poles.""" return emptyParameterList() 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 66760b6..59d16ad 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,361 +1,361 @@ #Copyright (C) 2018-2020 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 ( GenericPivotedGreedyApproximantBase, GenericPivotedGreedyApproximantPoleMatch) 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 ( RationalInterpolantGreedyPivotedPoleMatch) 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 from rrompy.utilities.parallel import poolRank, recv __all__ = ['RationalInterpolantGreedyPivotedGreedyPoleMatch'] class RationalInterpolantGreedyPivotedGreedyBase( GenericPivotedGreedyApproximantBase): @property def sampleBatchSize(self): """Value of sampleBatchSize.""" return 1 @property def sampleBatchIdx(self): """Value of sampleBatchIdx.""" return self.S 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), 3) 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[j - len(mus)])): self.samplingEngine.nextSample(mu) if self._isLastSampleCollinear(): vbMng(self, "MAIN", ("Collinearity above tolerance detected. Starting " "preemptive greedy loop termination."), 3) 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 _setSampleBatch(self, maxS:int): return self.S 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() self.samplingEngine.scaleFactor = self.scaleFactorDer musPivot = self.samplerTrainSet.generatePoints(self.S) while len(musPivot) > self.S: musPivot.pop() muTestBasePivot = self.samplerPivot.generatePoints(self.nTestPoints, False) idxPop = pruneSamples(self.mapParameterListPivot(muTestBasePivot), self.mapParameterListPivot(musPivot), 1e-10 * self.scaleFactorPivot[0]) muTestBasePivot.pop(idxPop) self._mus = emptyParameterList() self.mus.reset((self.S - 1, self.HFEngine.npar)) self.muTest = emptyParameterList() self.muTest.reset((len(muTestBasePivot) + 1, self.HFEngine.npar)) self.mus.data[:, self.directionPivot] = musPivot[: -1] self.mus.data[:, self.directionMarginal] = np.repeat(self.muMargLoc, self.S - 1, axis = 0) self.muTest.data[: -1, self.directionPivot] = muTestBasePivot.data self.muTest.data[-1, self.directionPivot] = musPivot[-1] self.muTest.data[:, self.directionMarginal] = np.repeat(self.muMargLoc, len(muTestBasePivot) + 1, axis = 0) 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), 3) self.samplingEngine.iterSample(self.mus) self._S = len(self.mus) self._approxParameters["S"] = self.S self.M, self.N = ("AUTO",) * 2 def setupApproxPivoted(self, mus:paramList) -> int: if self.checkComputedApproxPivoted(): return -1 RROMPyAssert(self._mode, message = "Cannot setup approximant.") vbMng(self, "INIT", "Setting up pivoted approximant.", 10) if not hasattr(self, "_plotEstPivot"): self._plotEstPivot = "NONE" idx, sizes, emptyCores = self._preSetupApproxPivoted(mus) S0 = copy(self.S) pMat, Ps, Qs, req, musA = None, [], [], [], None if len(idx) == 0: vbMng(self, "MAIN", "Idling.", 45) if self.storeAllSamples: self.storeSamples() pL, pT, mT = recv(source = 0, tag = poolRank()) pMat = np.empty((pL, 0), dtype = pT) musA = np.empty((0, self.mu0.shape[1]), dtype = mT) else: for i in idx: self.muMargLoc = mus[[i]] vbMng(self, "MAIN", "Building marginal model no. {} at " - "{}.".format(i + 1, self.muMargLoc), 25) + "{}.".format(i + 1, self.muMargLoc[0]), 25) self.samplingEngine.resetHistory() self.trainedModel = None self.verbosity -= 5 - self.samplingEngine.verbosity -= 10 + self.samplingEngine.verbosity -= 5 RationalInterpolantGreedy.setupApprox(self, self._plotEstPivot) self.verbosity += 5 - self.samplingEngine.verbosity += 10 + self.samplingEngine.verbosity += 5 if self.storeAllSamples: self.storeSamples(i + self._nmusOld) pMat, req, musA = self._localPivotedResult(pMat, req, emptyCores, musA) Ps += [copy(self.trainedModel.data.P)] Qs += [copy(self.trainedModel.data.Q)] self._S = S0 del self.muMargLoc for r in req: r.wait() self._postSetupApproxPivoted(musA, pMat, Ps, Qs, sizes) vbMng(self, "DEL", "Done setting up pivoted 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 class RationalInterpolantGreedyPivotedGreedyPoleMatch( RationalInterpolantGreedyPivotedGreedyBase, GenericPivotedGreedyApproximantPoleMatch, RationalInterpolantGreedyPivotedPoleMatch): """ ROM greedy pivoted greedy rational interpolant computation for parametric problems (with pole matching). 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': kind of snapshots orthogonalization; allowed values include 0, 1/2, and 1; defaults to 1, i.e. POD; - 'scaleFactorDer': scaling factors for derivative computation; defaults to 'AUTO'; - 'matchState': whether to match the system state rather than the system output; defaults to False; - 'matchingWeight': weight for pole matching optimization; defaults to 1; - 'matchingChordalRadius': radius to be used in chordal metric for poles and residues; if <= 0, Euclidean metric is used; if 'AUTO', automatically selected; defaults to -1; - 'matchingShared': required ratio of marginal points to share resonance; defaults to 1.; - 'badPoleCorrection': strategy for correction of bad poles; available values include 'ERASE', 'RATIONAL', and 'POLYNOMIAL'; defaults to 'ERASE'; - 'matchingWeightError': weight for pole matching optimization in error estimation; defaults to 0; - 'S': total number of pivot samples current approximant relies upon; - 'samplerPivot': pivot sample point generator; - 'SMarginal': number of starting marginal samples; - 'samplerMarginal': marginal sample point generator via sparse grid; - 'errorEstimatorKindMarginal': kind of marginal error estimator; available values include 'LOOK_AHEAD' and 'LOOK_AHEAD_RECOVER'; defaults to 'NONE'; - 'polybasis': type of polynomial basis for pivot interpolation; defaults to 'MONOMIAL'; - 'polybasisMarginal': type of polynomial basis for marginal interpolation; allowed values include 'MONOMIAL_*', 'CHEBYSHEV_*', 'LEGENDRE_*', 'NEARESTNEIGHBOR', and 'PIECEWISE_LINEAR_*'; defaults to 'MONOMIAL'; - 'paramsMarginal': dictionary of parameters for marginal interpolation; include: . 'MMarginal': degree of marginal interpolant; defaults to 'AUTO', i.e. maximum allowed; not for 'NEARESTNEIGHBOR' or 'PIECEWISE_LINEAR_*'; . 'nNeighborsMarginal': number of marginal nearest neighbors; defaults to 1; only for 'NEARESTNEIGHBOR'; . 'polydegreetypeMarginal': type of polynomial degree for marginal; defaults to 'TOTAL'; not for 'NEARESTNEIGHBOR' or 'PIECEWISE_LINEAR_*'; . 'interpTolMarginal': tolerance for marginal interpolation; defaults to None; not for 'NEARESTNEIGHBOR' or 'PIECEWISE_LINEAR_*'; . 'radialDirectionalWeightsMarginalAdapt': bounds for adaptive rescaling of marginal radial basis weights; only for radial basis. - '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; - 'samplerTrainSet': training sample points generator; defaults to samplerPivot; - 'errorEstimatorKind': kind of error estimator; available values include 'AFFINE', 'DISCREPANCY', 'LOOK_AHEAD', 'LOOK_AHEAD_RES', and 'NONE'; defaults to 'NONE'; - 'greedyTolMarginal': uniform error tolerance for marginal greedy algorithm; defaults to 1e-1; - 'maxIterMarginal': maximum number of marginal greedy steps; defaults to 1e2; - 'radialDirectionalWeightsMarginal': radial basis weights for marginal interpolant; defaults to 1; - 'functionalSolve': strategy for minimization of denominator functional; allowed values include 'NORM', 'DOMINANT', 'BARYCENTRIC[_NORM]', and 'BARYCENTRIC_AVERAGE' (check pdf in main folder for explanation); defaults to 'NORM'; - 'interpTol': tolerance for pivot interpolation; defaults to None; - 'QTol': tolerance for robust rational denominator management; defaults to 0. Defaults to empty dict. 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': kind of snapshots orthogonalization; - 'scaleFactorDer': scaling factors for derivative computation; - 'matchState': whether to match the system state rather than the system output; - 'matchingWeight': weight for pole matching optimization; - 'matchingChordalRadius': radius to be used in chordal metric for poles and residues; - 'matchingShared': required ratio of marginal points to share resonance; - 'badPoleCorrection': strategy for correction of bad poles; - 'matchingWeightError': weight for pole matching optimization in error estimation; - 'errorEstimatorKindMarginal': kind of marginal error estimator; - 'polybasis': type of polynomial basis for pivot interpolation; - 'polybasisMarginal': type of polynomial basis for marginal interpolation; - 'paramsMarginal': dictionary of parameters for marginal interpolation; include: . 'MMarginal': degree of marginal interpolant; . 'nNeighborsMarginal': number of marginal nearest neighbors; . 'polydegreetypeMarginal': type of polynomial degree for marginal; . 'interpTolMarginal': tolerance for marginal interpolation; . 'radialDirectionalWeightsMarginalAdapt': bounds for adaptive rescaling of marginal radial basis weights. - 'greedyTol': uniform error tolerance for greedy algorithm; - 'collinearityTol': collinearity tolerance for greedy algorithm; - 'maxIter': maximum number of greedy steps; - 'nTestPoints': number of test points; - 'samplerTrainSet': training sample points generator; - 'errorEstimatorKind': kind of error estimator; - 'greedyTolMarginal': uniform error tolerance for marginal greedy algorithm; - 'maxIterMarginal': maximum number of marginal greedy steps; - 'radialDirectionalWeightsMarginal': radial basis weights for marginal interpolant; - 'functionalSolve': strategy for minimization of denominator functional; - 'interpTol': tolerance for pivot interpolation; - 'QTol': 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; - 'samplerMarginal': marginal sample point generator via sparse grid. verbosity: Verbosity level. POD: Kind of snapshots orthogonalization. scaleFactorDer: Scaling factors for derivative computation. matchState: Whether to match the system state rather than the system output. matchingWeight: Weight for pole matching optimization. matchingChordalRadius: Radius to be used in chordal metric for poles and residues. matchingShared: Required ratio of marginal points to share resonance. badPoleCorrection: Strategy for correction of bad poles. matchingWeightError: Weight for pole matching optimization 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. samplerMarginal: Marginal sample point generator via sparse grid. errorEstimatorKindMarginal: Kind of marginal error estimator. polybasis: Type of polynomial basis for pivot interpolation. polybasisMarginal: Type of polynomial basis for marginal interpolation. paramsMarginal: Dictionary of parameters 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. samplerTrainSet: training sample points generator. errorEstimatorKind: kind of error estimator. greedyTolMarginal: Uniform error tolerance for marginal greedy algorithm. maxIterMarginal: Maximum number of marginal greedy steps. radialDirectionalWeightsMarginal: Radial basis weights for marginal interpolant. functionalSolve: Strategy for minimization of denominator functional. interpTol: Tolerance for pivot interpolation. QTol: Tolerance for robust rational denominator management. muBounds: list of bounds for pivot parameter values. muBoundsMarginal: list of bounds for marginal parameter values. samplingEngine: Sampling engine. uHF: High fidelity solution(s) with parameter(s) lastSolvedHF as sampleList. lastSolvedHF: Parameter(s) corresponding to last computed high fidelity solution(s) as parameterList. uApproxReduced: Reduced approximate solution(s) with parameter(s) lastSolvedApprox as sampleList. lastSolvedApproxReduced: Parameter(s) corresponding to last computed reduced approximate solution(s) as parameterList. uApprox: Approximate solution(s) with parameter(s) lastSolvedApprox as sampleList. lastSolvedApprox: Parameter(s) corresponding to last computed approximate solution(s) as parameterList. """ 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 1aad812..6160230 100644 --- a/rrompy/reduction_methods/pivoted/greedy/rational_interpolant_pivoted_greedy.py +++ b/rrompy/reduction_methods/pivoted/greedy/rational_interpolant_pivoted_greedy.py @@ -1,295 +1,295 @@ # Copyright (C) 2018-2020 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 ( GenericPivotedGreedyApproximantBase, GenericPivotedGreedyApproximantPoleMatch) from rrompy.reduction_methods.standard import RationalInterpolant from rrompy.reduction_methods.pivoted import ( RationalInterpolantPivotedPoleMatch) 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 emptyParameterList from rrompy.utilities.parallel import poolRank, recv __all__ = ['RationalInterpolantPivotedGreedyPoleMatch'] class RationalInterpolantPivotedGreedyBase( GenericPivotedGreedyApproximantBase): 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.samplingEngine.scaleFactor = self.scaleFactorDer if not hasattr(self, "musPivot") or len(self.musPivot) != self.S: self.musPivot = self.samplerPivot.generatePoints(self.S) while len(self.musPivot) > self.S: self.musPivot.pop() musLoc = emptyParameterList() musLoc.reset((self.S, self.HFEngine.npar)) self.samplingEngine.resetHistory() musLoc.data[:, self.directionPivot] = self.musPivot.data musLoc.data[:, self.directionMarginal] = np.repeat(self.muMargLoc, self.S, axis = 0) self.samplingEngine.iterSample(musLoc) vbMng(self, "DEL", "Done computing snapshots.", 5) self._m_selfmus = copy(musLoc) self._mus = self.musPivot self._m_HFEparameterMap = copy(self.HFEngine.parameterMap) self.HFEngine.parameterMap = { "F": [self.HFEngine.parameterMap["F"][self.directionPivot[0]]], "B": [self.HFEngine.parameterMap["B"][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 pivoted approximant.", 10) idx, sizes, emptyCores = self._preSetupApproxPivoted(mus) pMat, Ps, Qs, req, musA = None, [], [], [], None if len(idx) == 0: vbMng(self, "MAIN", "Idling.", 45) if self.storeAllSamples: self.storeSamples() pL, pT, mT = recv(source = 0, tag = poolRank()) pMat = np.empty((pL, 0), dtype = pT) musA = np.empty((0, self.mu0.shape[1]), dtype = mT) else: for i in idx: self.muMargLoc = mus[[i]] vbMng(self, "MAIN", "Building marginal model no. {} at " - "{}.".format(i + 1, self.muMargLoc), 25) + "{}.".format(i + 1, self.muMargLoc[0]), 25) self.samplingEngine.resetHistory() self.trainedModel = None self.verbosity -= 5 self.samplingEngine.verbosity -= 5 RationalInterpolant.setupApprox(self) self.verbosity += 5 self.samplingEngine.verbosity += 5 self._mus = self._m_selfmus self.HFEngine.parameterMap = self._m_HFEparameterMap del self._m_selfmus, self._m_HFEparameterMap if self.storeAllSamples: self.storeSamples(i + self._nmusOld) pMat, req, musA = self._localPivotedResult(pMat, req, emptyCores, musA) Ps += [copy(self.trainedModel.data.P)] Qs += [copy(self.trainedModel.data.Q)] del self.muMargLoc for r in req: r.wait() self._postSetupApproxPivoted(musA, pMat, Ps, Qs, sizes) vbMng(self, "DEL", "Done setting up pivoted approximant.", 10) return 0 class RationalInterpolantPivotedGreedyPoleMatch( RationalInterpolantPivotedGreedyBase, GenericPivotedGreedyApproximantPoleMatch, RationalInterpolantPivotedPoleMatch): """ ROM pivoted greedy rational interpolant computation for parametric problems (with pole matching). 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': kind of snapshots orthogonalization; allowed values include 0, 1/2, and 1; defaults to 1, i.e. POD; - 'scaleFactorDer': scaling factors for derivative computation; defaults to 'AUTO'; - 'matchState': whether to match the system state rather than the system output; defaults to False; - 'matchingWeight': weight for pole matching optimization; defaults to 1; - 'matchingChordalRadius': radius to be used in chordal metric for poles and residues; if <= 0, Euclidean metric is used; if 'AUTO', automatically selected; defaults to -1; - 'matchingShared': required ratio of marginal points to share resonance; defaults to 1.; - 'badPoleCorrection': strategy for correction of bad poles; available values include 'ERASE', 'RATIONAL', and 'POLYNOMIAL'; defaults to 'ERASE'; - 'matchingWeightError': weight for pole matching optimization in error estimation; defaults to 0; - 'S': total number of pivot samples current approximant relies upon; - 'samplerPivot': pivot sample point generator; - 'SMarginal': number of starting marginal samples; - 'samplerMarginal': marginal sample point generator via sparse grid; - 'errorEstimatorKindMarginal': kind of marginal error estimator; available values include 'LOOK_AHEAD' and 'LOOK_AHEAD_RECOVER'; defaults to 'NONE'; - 'polybasis': type of polynomial basis for pivot interpolation; defaults to 'MONOMIAL'; - 'polybasisMarginal': type of polynomial basis for marginal interpolation; allowed values include 'MONOMIAL_*', 'CHEBYSHEV_*', 'LEGENDRE_*', 'NEARESTNEIGHBOR', and 'PIECEWISE_LINEAR_*'; defaults to 'MONOMIAL'; - 'paramsMarginal': dictionary of parameters for marginal interpolation; include: . 'MMarginal': degree of marginal interpolant; defaults to 'AUTO', i.e. maximum allowed; not for 'NEARESTNEIGHBOR' or 'PIECEWISE_LINEAR_*'; . 'nNeighborsMarginal': number of marginal nearest neighbors; defaults to 1; only for 'NEARESTNEIGHBOR'; . 'polydegreetypeMarginal': type of polynomial degree for marginal; defaults to 'TOTAL'; not for 'NEARESTNEIGHBOR' or 'PIECEWISE_LINEAR_*'; . 'interpTolMarginal': tolerance for marginal interpolation; defaults to None; not for 'NEARESTNEIGHBOR' or 'PIECEWISE_LINEAR_*'; . 'radialDirectionalWeightsMarginalAdapt': bounds for adaptive rescaling of marginal radial basis weights; only for radial basis. - '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; - 'greedyTolMarginal': uniform error tolerance for marginal greedy algorithm; defaults to 1e-1; - 'maxIterMarginal': maximum number of marginal greedy steps; defaults to 1e2; - 'radialDirectionalWeights': radial basis weights for pivot numerator; defaults to 1; - 'radialDirectionalWeightsAdapt': bounds for adaptive rescaling of radial basis weights; defaults to [-1, -1]; - 'radialDirectionalWeightsMarginal': radial basis weights for marginal interpolant; defaults to 1; - 'functionalSolve': strategy for minimization of denominator functional; allowed values include 'NORM', 'DOMINANT', 'BARYCENTRIC[_NORM]', and 'BARYCENTRIC_AVERAGE' (check pdf in main folder for explanation); defaults to 'NORM'; - 'interpTol': tolerance for pivot interpolation; defaults to None; - 'QTol': tolerance for robust rational denominator management; defaults to 0. Defaults to empty dict. 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': kind of snapshots orthogonalization; - 'scaleFactorDer': scaling factors for derivative computation; - 'matchState': whether to match the system state rather than the system output; - 'matchingWeight': weight for pole matching optimization; - 'matchingChordalRadius': radius to be used in chordal metric for poles and residues; - 'matchingShared': required ratio of marginal points to share resonance; - 'badPoleCorrection': strategy for correction of bad poles; - 'matchingWeightError': weight for pole matching optimization in error estimation; - 'errorEstimatorKindMarginal': kind of marginal error estimator; - 'polybasis': type of polynomial basis for pivot interpolation; - 'polybasisMarginal': type of polynomial basis for marginal interpolation; - 'paramsMarginal': dictionary of parameters for marginal interpolation; include: . 'MMarginal': degree of marginal interpolant; . 'nNeighborsMarginal': number of marginal nearest neighbors; . 'polydegreetypeMarginal': type of polynomial degree for marginal; . 'interpTolMarginal': tolerance for marginal interpolation; . 'radialDirectionalWeightsMarginalAdapt': bounds for adaptive rescaling of marginal radial basis weights. - 'M': degree of rational interpolant numerator; - 'N': degree of rational interpolant denominator; - 'greedyTolMarginal': uniform error tolerance for marginal greedy algorithm; - 'maxIterMarginal': maximum number of marginal greedy steps; - 'radialDirectionalWeights': radial basis weights for pivot numerator; - 'radialDirectionalWeightsAdapt': bounds for adaptive rescaling of radial basis weights; - 'radialDirectionalWeightsMarginal': radial basis weights for marginal interpolant; - 'functionalSolve': strategy for minimization of denominator functional; - 'interpTol': tolerance for pivot interpolation; - 'QTol': 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; - 'samplerMarginal': marginal sample point generator via sparse grid. verbosity: Verbosity level. POD: Kind of snapshots orthogonalization. scaleFactorDer: Scaling factors for derivative computation. matchState: Whether to match the system state rather than the system output. matchingWeight: Weight for pole matching optimization. matchingChordalRadius: Radius to be used in chordal metric for poles and residues. matchingShared: Required ratio of marginal points to share resonance. badPoleCorrection: Strategy for correction of bad poles. matchingWeightError: Weight for pole matching optimization 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. samplerMarginal: Marginal sample point generator via sparse grid. errorEstimatorKindMarginal: Kind of marginal error estimator. polybasis: Type of polynomial basis for pivot interpolation. polybasisMarginal: Type of polynomial basis for marginal interpolation. paramsMarginal: Dictionary of parameters for marginal interpolation. M: Degree of rational interpolant numerator. N: Degree of rational interpolant denominator. greedyTolMarginal: Uniform error tolerance for marginal greedy algorithm. maxIterMarginal: Maximum number of marginal greedy steps. radialDirectionalWeights: Radial basis weights for pivot numerator. radialDirectionalWeightsAdapt: Bounds for adaptive rescaling of radial basis weights. radialDirectionalWeightsMarginal: Radial basis weights for marginal interpolant. functionalSolve: Strategy for minimization of denominator functional. interpTol: Tolerance for pivot interpolation. QTol: Tolerance for robust rational denominator management. muBounds: list of bounds for pivot parameter values. muBoundsMarginal: list of bounds for marginal parameter values. samplingEngine: Sampling engine. uHF: High fidelity solution(s) with parameter(s) lastSolvedHF as sampleList. lastSolvedHF: Parameter(s) corresponding to last computed high fidelity solution(s) as parameterList. uApproxReduced: Reduced approximate solution(s) with parameter(s) lastSolvedApprox as sampleList. lastSolvedApproxReduced: Parameter(s) corresponding to last computed reduced approximate solution(s) as parameterList. uApprox: Approximate solution(s) with parameter(s) lastSolvedApprox as sampleList. lastSolvedApprox: Parameter(s) corresponding to last computed approximate solution(s) as parameterList. """ diff --git a/rrompy/reduction_methods/pivoted/rational_interpolant_greedy_pivoted.py b/rrompy/reduction_methods/pivoted/rational_interpolant_greedy_pivoted.py index 2e7cb2b..a05143c 100644 --- a/rrompy/reduction_methods/pivoted/rational_interpolant_greedy_pivoted.py +++ b/rrompy/reduction_methods/pivoted/rational_interpolant_greedy_pivoted.py @@ -1,571 +1,571 @@ # Copyright (C) 2018-2020 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 (GenericPivotedApproximantBase, GenericPivotedApproximantNoMatch, GenericPivotedApproximantPoleMatch) from .gather_pivoted_approximant import gatherPivotedApproximant 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.poly_fitting.polynomial import polyvander as pv from rrompy.utilities.exception_manager import RROMPyException, RROMPyAssert from rrompy.parameter import emptyParameterList from rrompy.utilities.parallel import poolRank, indicesScatter, isend, recv __all__ = ['RationalInterpolantGreedyPivotedNoMatch', 'RationalInterpolantGreedyPivotedPoleMatch'] class RationalInterpolantGreedyPivotedBase(GenericPivotedApproximantBase, RationalInterpolantGreedy): def __init__(self, *args, **kwargs): self._preInit() super().__init__(*args, **kwargs) if self.nparPivot > 1: self.HFEngine._ignoreResidues = 1 self._postInit() @property def tModelType(self): if hasattr(self, "_temporaryPivot"): return RationalInterpolantGreedy.tModelType.fget(self) return super().tModelType 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_selfmus = copy(self.mus) self._m_HFEparameterMap = copy(self.HFEngine.parameterMap) self._mus = self.checkParameterListPivot( self.mus(self.directionPivot)) self.HFEngine.parameterMap = { "F": [self.HFEngine.parameterMap["F"][self.directionPivot[0]]], "B": [self.HFEngine.parameterMap["B"][self.directionPivot[0]]]} else: self._mus = self._m_selfmus self.HFEngine.parameterMap = self._m_HFEparameterMap del self._m_selfmus, self._m_HFEparameterMap 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.parameterMap = self.HFEngine.parameterMap 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 = self.checkParameterList(musUniqueCNAux) 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 self.trainedModel.data.Q._dirPivot = self.directionPivot[0] self.trainedModel.data.P._dirPivot = self.directionPivot[0] # tell greedy error estimator that operator / RHS is pivot-affine if hasattr(self.HFEngine.A, "is_affine"): self._A_is_affine = self.HFEngine.A.is_affine else: self._A_is_affine = 0 if hasattr(self.HFEngine.b, "is_affine"): self._b_is_affine = self.HFEngine.b.is_affine else: self._b_is_affine = 0 if self._A_is_affine >= 1 / 2 and self._b_is_affine >= 1 / 2: self._affine_lvl += [1 / 2] else: self._temporaryPivot = 1 self.trainedModel.data.mu0 = self.checkParameterListPivot( self.mu0(self.directionPivot)) self.trainedModel.data.scaleFactor = self.scaleFactor self.trainedModel.data.parameterMap = { "F": [self.HFEngine.parameterMap["F"][self.directionPivot[0]]], "B": [self.HFEngine.parameterMap["B"][self.directionPivot[0]]]} self._musUniqueCN = copy(self._m_musUniqueCN) self._derIdxs = copy(self._m_derIdxs) del self._m_musUniqueCN, self._m_derIdxs del self.trainedModel.data.Q._dirPivot del self.trainedModel.data.P._dirPivot if self._A_is_affine >= 1 / 2 and self._b_is_affine >= 1 / 2: self._affine_lvl.pop() del self._A_is_affine, self._b_is_affine self.trainedModel.data.npar = self.npar def errorEstimator(self, mus:Np1D, return_max : bool = False) -> Np1D: """Standard residual-based error estimator.""" setupOK = self.setupApproxLocal() if setupOK > 0: err = np.empty(len(mus)) err[:] = np.nan if not return_max: return err return err, [- setupOK], np.nan 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.") self._S = self._setSampleBatch(self.S) self.resetSamples() self.samplingEngine.scaleFactor = self.scaleFactorDer musPivot = self.samplerTrainSet.generatePoints(self.S) while len(musPivot) > self.S: musPivot.pop() muTestPivot = self.samplerPivot.generatePoints(self.nTestPoints, False) idxPop = pruneSamples(self.mapParameterListPivot(muTestPivot), self.mapParameterListPivot(musPivot), 1e-10 * self.scaleFactorPivot[0]) muTestPivot.pop(idxPop) self._mus = emptyParameterList() self.mus.reset((self.S - 1, self.HFEngine.npar)) self.muTest = emptyParameterList() self.muTest.reset((len(muTestPivot) + 1, self.HFEngine.npar)) self.mus.data[:, self.directionPivot] = musPivot[: -1] self.mus.data[:, self.directionMarginal] = np.repeat(self.muMargLoc, self.S - 1, axis = 0) self.muTest.data[: -1, self.directionPivot] = muTestPivot.data self.muTest.data[-1, self.directionPivot] = musPivot[-1] self.muTest.data[:, self.directionMarginal] = np.repeat(self.muMargLoc, len(muTestPivot) + 1, axis = 0) 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), 3) self.samplingEngine.iterSample(self.mus) self._S = len(self.mus) self._approxParameters["S"] = self.S self.M, self.N = ("AUTO",) * 2 def setupApproxLocal(self) -> int: """Compute rational interpolant.""" self._marginalizeMiscellanea(True) setupOK = super().setupApproxLocal() self._marginalizeMiscellanea(False) return setupOK 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.computeScaleFactor() self._musMarginal = self.samplerMarginal.generatePoints(self.SMarginal) while len(self.musMarginal) > self.SMarginal: self.musMarginal.pop() S0 = copy(self.S) idx, sizes = indicesScatter(len(self.musMarginal), return_sizes = True) pMat, Ps, Qs, mus = None, [], [], None req, emptyCores = [], np.where(sizes == 0)[0] if len(idx) == 0: vbMng(self, "MAIN", "Idling.", 25) if self.storeAllSamples: self.storeSamples() pL, pT, mT = recv(source = 0, tag = poolRank()) pMat = np.empty((pL, 0), dtype = pT) mus = np.empty((0, self.mu0.shape[1]), dtype = mT) else: _scaleFactorOldPivot = copy(self.scaleFactor) self.scaleFactor = self.scaleFactorPivot self._temporaryPivot = 1 for i in idx: self.muMargLoc = self.musMarginal[[i]] vbMng(self, "MAIN", "Building marginal model no. {} at {}.".format(i + 1, - self.muMargLoc), 5) + self.musMarginal[i]), 5) self.samplingEngine.resetHistory() self.trainedModel = None self.verbosity -= 5 - self.samplingEngine.verbosity -= 10 + self.samplingEngine.verbosity -= 5 RationalInterpolantGreedy.setupApprox(self, *args, **kwargs) self.verbosity += 5 - self.samplingEngine.verbosity += 10 + self.samplingEngine.verbosity += 5 if self.storeAllSamples: self.storeSamples(i) musi = self.samplingEngine.mus pMati = self.samplingEngine.projectionMatrix if not hasattr(self, "matchState") or not self.matchState: if self.POD == 1 and not ( hasattr(self.HFEngine.C, "is_mu_independent") and self.HFEngine.C.is_mu_independent in self._output_lvl): raise RROMPyException(("Cannot apply mu-dependent C " "to orthonormalized samples.")) vbMng(self, "INIT", "Extracting system output from state.", 35) pMatiEff = None for j, mu in enumerate(musi): pMij = np.expand_dims(self.HFEngine.applyC( pMati[:, j], mu), -1) if pMatiEff is None: pMatiEff = np.array(pMij) else: pMatiEff = np.append(pMatiEff, pMij, axis = 1) pMati = pMatiEff vbMng(self, "DEL", "Done extracting system output.", 35) if pMat is None: mus = copy(musi.data) pMat = copy(pMati) if i == 0: for dest in emptyCores: req += [isend((len(pMat), pMat.dtype, mus.dtype), dest = dest, tag = dest)] else: mus = np.vstack((mus, musi.data)) pMat = np.hstack((pMat, pMati)) Ps += [copy(self.trainedModel.data.P)] Qs += [copy(self.trainedModel.data.Q)] self._S = S0 del self._temporaryPivot, self.muMargLoc self.scaleFactor = _scaleFactorOldPivot for r in req: r.wait() pMat, Ps, Qs, mus, nsamples = gatherPivotedApproximant(pMat, Ps, Qs, mus, sizes, self.polybasis) self._mus = self.checkParameterList(mus) Psupp = np.append(0, np.cumsum(nsamples)) self._setupTrainedModel(pMat, forceNew = True) self.trainedModel.data.Qs, self.trainedModel.data.Ps = Qs, Ps self.trainedModel.data.Psupp = list(Psupp[: -1]) self._preliminaryMarginalFinalization() self._finalizeMarginalization() vbMng(self, "DEL", "Done setting up approximant.", 5) return 0 class RationalInterpolantGreedyPivotedNoMatch( RationalInterpolantGreedyPivotedBase, GenericPivotedApproximantNoMatch): """ ROM pivoted rational interpolant (without 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': kind of snapshots orthogonalization; allowed values include 0, 1/2, and 1; defaults to 1, i.e. POD; - 'scaleFactorDer': scaling factors for derivative computation; defaults to 'AUTO'; - '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'; - '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; - 'samplerTrainSet': training sample points generator; defaults to samplerPivot; - 'errorEstimatorKind': kind of error estimator; available values include 'AFFINE', 'DISCREPANCY', 'LOOK_AHEAD', 'LOOK_AHEAD_RES', and 'NONE'; defaults to 'NONE'; - 'radialDirectionalWeightsMarginal': radial basis weights for marginal interpolant; defaults to 1; - 'functionalSolve': strategy for minimization of denominator functional; allowed values include 'NORM', 'DOMINANT', 'BARYCENTRIC[_NORM]', and 'BARYCENTRIC_AVERAGE' (check pdf in main folder for explanation); defaults to 'NORM'; - 'interpTol': tolerance for pivot interpolation; defaults to None; - 'QTol': tolerance for robust rational denominator management; defaults to 0. Defaults to empty dict. 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': kind of snapshots orthogonalization; - 'scaleFactorDer': scaling factors for derivative computation; - 'polybasis': type of polynomial basis for pivot 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; - 'samplerTrainSet': training sample points generator; - 'errorEstimatorKind': kind of error estimator; - 'radialDirectionalWeightsMarginal': radial basis weights for marginal interpolant; - 'functionalSolve': strategy for minimization of denominator functional; - 'interpTol': tolerance for pivot interpolation; - 'QTol': 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; - 'samplerMarginal': marginal sample point generator. verbosity: Verbosity level. POD: Kind of snapshots orthogonalization. scaleFactorDer: Scaling factors for derivative computation. 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. 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. samplerTrainSet: training sample points generator. errorEstimatorKind: kind of error estimator. radialDirectionalWeightsMarginal: Radial basis weights for marginal interpolant. functionalSolve: Strategy for minimization of denominator functional. interpTol: Tolerance for pivot interpolation. QTol: Tolerance for robust rational denominator management. muBounds: list of bounds for pivot parameter values. muBoundsMarginal: list of bounds for marginal parameter values. samplingEngine: Sampling engine. uHF: High fidelity solution(s) with parameter(s) lastSolvedHF as sampleList. lastSolvedHF: Parameter(s) corresponding to last computed high fidelity solution(s) as parameterList. uApproxReduced: Reduced approximate solution(s) with parameter(s) lastSolvedApprox as sampleList. lastSolvedApproxReduced: Parameter(s) corresponding to last computed reduced approximate solution(s) as parameterList. uApprox: Approximate solution(s) with parameter(s) lastSolvedApprox as sampleList. lastSolvedApprox: Parameter(s) corresponding to last computed approximate solution(s) as parameterList. Q: Numpy 1D vector containing complex coefficients of approximant denominator. P: Numpy 2D vector whose columns are FE dofs of coefficients of approximant numerator. """ class RationalInterpolantGreedyPivotedPoleMatch( RationalInterpolantGreedyPivotedBase, GenericPivotedApproximantPoleMatch): """ 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': kind of snapshots orthogonalization; allowed values include 0, 1/2, and 1; defaults to 1, i.e. POD; - 'scaleFactorDer': scaling factors for derivative computation; defaults to 'AUTO'; - 'matchState': whether to match the system state rather than the system output; defaults to False; - 'matchingWeight': weight for pole matching optimization; defaults to 1; - 'matchingChordalRadius': radius to be used in chordal metric for poles and residues; if <= 0, Euclidean metric is used; if 'AUTO', automatically selected; defaults to -1; - 'matchingShared': required ratio of marginal points to share resonance; defaults to 1.; - 'badPoleCorrection': strategy for correction of bad poles; available values include 'ERASE', 'RATIONAL', and 'POLYNOMIAL'; defaults to 'ERASE'; - '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_*', 'LEGENDRE_*', 'NEARESTNEIGHBOR', and 'PIECEWISE_LINEAR_*'; defaults to 'MONOMIAL'; - 'paramsMarginal': dictionary of parameters for marginal interpolation; include: . 'MMarginal': degree of marginal interpolant; defaults to 'AUTO', i.e. maximum allowed; not for 'NEARESTNEIGHBOR' or 'PIECEWISE_LINEAR_*'; . 'nNeighborsMarginal': number of marginal nearest neighbors; defaults to 1; only for 'NEARESTNEIGHBOR'; . 'polydegreetypeMarginal': type of polynomial degree for marginal; defaults to 'TOTAL'; not for 'NEARESTNEIGHBOR' or 'PIECEWISE_LINEAR_*'; . 'interpTolMarginal': tolerance for marginal interpolation; defaults to None; not for 'NEARESTNEIGHBOR' or 'PIECEWISE_LINEAR_*'; . 'radialDirectionalWeightsMarginalAdapt': bounds for adaptive rescaling of marginal radial basis weights; only for radial basis. - '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; - 'samplerTrainSet': training sample points generator; defaults to samplerPivot; - 'errorEstimatorKind': kind of error estimator; available values include 'AFFINE', 'DISCREPANCY', 'LOOK_AHEAD', 'LOOK_AHEAD_RES', and 'NONE'; defaults to 'NONE'; - 'radialDirectionalWeightsMarginal': radial basis weights for marginal interpolant; defaults to 1; - 'functionalSolve': strategy for minimization of denominator functional; allowed values include 'NORM', 'DOMINANT', 'BARYCENTRIC[_NORM]', and 'BARYCENTRIC_AVERAGE' (check pdf in main folder for explanation); defaults to 'NORM'; - 'interpTol': tolerance for pivot interpolation; defaults to None; - 'QTol': tolerance for robust rational denominator management; defaults to 0. Defaults to empty dict. 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': kind of snapshots orthogonalization; - 'scaleFactorDer': scaling factors for derivative computation; - 'matchState': whether to match the system state rather than the system output; - 'matchingWeight': weight for pole matching optimization; - 'matchingChordalRadius': radius to be used in chordal metric for poles and residues; - 'matchingShared': required ratio of marginal points to share resonance; - 'badPoleCorrection': strategy for correction of bad poles; - 'polybasis': type of polynomial basis for pivot interpolation; - 'polybasisMarginal': type of polynomial basis for marginal interpolation; - 'paramsMarginal': dictionary of parameters for marginal interpolation; include: . 'MMarginal': degree of marginal interpolant; . 'nNeighborsMarginal': number of marginal nearest neighbors; . 'polydegreetypeMarginal': type of polynomial degree for marginal; . 'interpTolMarginal': tolerance for marginal interpolation; . 'radialDirectionalWeightsMarginalAdapt': bounds for adaptive rescaling of marginal radial basis weights. - 'greedyTol': uniform error tolerance for greedy algorithm; - 'collinearityTol': collinearity tolerance for greedy algorithm; - 'maxIter': maximum number of greedy steps; - 'nTestPoints': number of test points; - 'samplerTrainSet': training sample points generator; - 'errorEstimatorKind': kind of error estimator; - 'radialDirectionalWeightsMarginal': radial basis weights for marginal interpolant; - 'functionalSolve': strategy for minimization of denominator functional; - 'interpTol': tolerance for pivot interpolation; - 'QTol': 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; - 'samplerMarginal': marginal sample point generator. verbosity: Verbosity level. POD: Kind of snapshots orthogonalization. scaleFactorDer: Scaling factors for derivative computation. matchState: Whether to match the system state rather than the system output. matchingWeight: Weight for pole matching optimization. matchingChordalRadius: Radius to be used in chordal metric for poles and residues. matchingShared: Required ratio of marginal points to share resonance. badPoleCorrection: Strategy for correction of bad poles. S: Total number of pivot samples current approximant relies upon. samplerPivot: Pivot sample point generator. SMarginal: Total number of marginal samples current approximant relies upon. samplerMarginal: Marginal sample point generator. polybasis: Type of polynomial basis for pivot interpolation. polybasisMarginal: Type of polynomial basis for marginal interpolation. paramsMarginal: Dictionary of parameters 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. samplerTrainSet: training sample points generator. errorEstimatorKind: kind of error estimator. radialDirectionalWeightsMarginal: Radial basis weights for marginal interpolant. functionalSolve: Strategy for minimization of denominator functional. interpTol: Tolerance for pivot interpolation. QTol: Tolerance for robust rational denominator management. muBounds: list of bounds for pivot parameter values. muBoundsMarginal: list of bounds for marginal parameter values. samplingEngine: Sampling engine. uHF: High fidelity solution(s) with parameter(s) lastSolvedHF as sampleList. lastSolvedHF: Parameter(s) corresponding to last computed high fidelity solution(s) as parameterList. uApproxReduced: Reduced approximate solution(s) with parameter(s) lastSolvedApprox as sampleList. lastSolvedApproxReduced: Parameter(s) corresponding to last computed reduced approximate solution(s) as parameterList. uApprox: Approximate solution(s) with parameter(s) lastSolvedApprox as sampleList. lastSolvedApprox: Parameter(s) corresponding to last computed approximate solution(s) as parameterList. Q: Numpy 1D vector containing complex coefficients of approximant denominator. P: Numpy 2D vector whose columns are FE dofs of coefficients of approximant numerator. """ def setupApprox(self, *args, **kwargs) -> int: if self.checkComputedApprox(): return -1 self.purgeparamsMarginal() setupOK = super().setupApprox(*args, **kwargs) if self.matchState: self._postApplyC() return setupOK diff --git a/rrompy/reduction_methods/pivoted/rational_interpolant_pivoted.py b/rrompy/reduction_methods/pivoted/rational_interpolant_pivoted.py index bb07351..dc66a48 100644 --- a/rrompy/reduction_methods/pivoted/rational_interpolant_pivoted.py +++ b/rrompy/reduction_methods/pivoted/rational_interpolant_pivoted.py @@ -1,491 +1,491 @@ # Copyright (C) 2018-2020 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 collections.abc import Iterable from copy import deepcopy as copy from .generic_pivoted_approximant import (GenericPivotedApproximantBase, GenericPivotedApproximantNoMatch, GenericPivotedApproximantPoleMatch) from .gather_pivoted_approximant import gatherPivotedApproximant from rrompy.reduction_methods.standard.rational_interpolant import ( RationalInterpolant) from rrompy.utilities.base import verbosityManager as vbMng from rrompy.utilities.numerical.hash_derivative import nextDerivativeIndices from rrompy.utilities.exception_manager import (RROMPyException, RROMPyAssert, RROMPyWarning) from rrompy.parameter import emptyParameterList from rrompy.utilities.parallel import poolRank, indicesScatter, isend, recv __all__ = ['RationalInterpolantPivotedNoMatch', 'RationalInterpolantPivotedPoleMatch'] class RationalInterpolantPivotedBase(GenericPivotedApproximantBase, RationalInterpolant): def __init__(self, *args, **kwargs): self._preInit() self._addParametersToList(toBeExcluded = ["polydegreetype"]) super().__init__(*args, **kwargs) if self.nparPivot > 1: self.HFEngine._ignoreResidues = 1 self._postInit() @property def scaleFactorDer(self): """Value of scaleFactorDer.""" if self._scaleFactorDer == "NONE": return 1. if self._scaleFactorDer == "AUTO": return self.scaleFactorPivot return self._scaleFactorDer @scaleFactorDer.setter def scaleFactorDer(self, scaleFactorDer): if isinstance(scaleFactorDer, (str,)): scaleFactorDer = scaleFactorDer.upper() elif isinstance(scaleFactorDer, Iterable): scaleFactorDer = list(scaleFactorDer) self._scaleFactorDer = scaleFactorDer self._approxParameters["scaleFactorDer"] = self._scaleFactorDer @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'.")) 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 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.computeScaleFactor() self.resetSamples() self.samplingEngine.scaleFactor = self.scaleFactorDer self.musPivot = self.samplerPivot.generatePoints(self.S) while len(self.musPivot) > self.S: self.musPivot.pop() self._musMarginal = self.samplerMarginal.generatePoints(self.SMarginal) while len(self.musMarginal) > self.SMarginal: self.musMarginal.pop() self._mus = emptyParameterList() self.mus.reset((self.S * self.SMarginal, self.HFEngine.npar)) self.mus.data[:, self.directionPivot] = np.tile(self.musPivot.data, (self.SMarginal, 1)) self.mus.data[:, self.directionMarginal] = np.repeat( self.musMarginal.data, self.S, axis = 0) N0 = copy(self.N) self._setupTrainedModel(np.zeros((0, 0)), forceNew = True) idx, sizes = indicesScatter(len(self.musMarginal), return_sizes = True) pMat, Ps, Qs = None, [], [] req, emptyCores = [], np.where(sizes == 0)[0] if len(idx) == 0: vbMng(self, "MAIN", "Idling.", 30) if self.storeAllSamples: self.storeSamples() pL, pT = recv(source = 0, tag = poolRank()) pMat = np.empty((pL, 0), dtype = pT) else: _scaleFactorOldPivot = copy(self.scaleFactor) self.scaleFactor = self.scaleFactorPivot self._temporaryPivot = 1 for i in idx: musi = self.mus[self.S * i : self.S * (i + 1)] vbMng(self, "MAIN", "Building marginal model no. {} at {}.".format(i + 1, self.musMarginal[i]), 5) vbMng(self, "INIT", "Starting computation of snapshots.", 10) self.samplingEngine.resetHistory() self.samplingEngine.iterSample(musi) vbMng(self, "DEL", "Done computing snapshots.", 10) self.verbosity -= 5 - self.samplingEngine.verbosity -= 10 + self.samplingEngine.verbosity -= 5 self._setupRational(self._setupDenominator()) self.verbosity += 5 - self.samplingEngine.verbosity += 10 + self.samplingEngine.verbosity += 5 if self.storeAllSamples: self.storeSamples(i) pMati = self.samplingEngine.projectionMatrix if not hasattr(self, "matchState") or not self.matchState: if self.POD == 1 and not ( hasattr(self.HFEngine.C, "is_mu_independent") and self.HFEngine.C.is_mu_independent in self._output_lvl): raise RROMPyException(("Cannot apply mu-dependent C " "to orthonormalized samples.")) vbMng(self, "INIT", "Extracting system output from state.", 35) pMatiEff = None for j, mu in enumerate(musi): pMij = np.expand_dims(self.HFEngine.applyC( pMati[:, j], mu), -1) if pMatiEff is None: pMatiEff = np.array(pMij) else: pMatiEff = np.append(pMatiEff, pMij, axis = 1) pMati = pMatiEff vbMng(self, "DEL", "Done extracting system output.", 35) if pMat is None: pMat = copy(pMati) if i == 0: for dest in emptyCores: req += [isend((len(pMat), pMat.dtype), dest = dest, tag = dest)] else: pMat = np.hstack((pMat, pMati)) Ps += [copy(self.trainedModel.data.P)] Qs += [copy(self.trainedModel.data.Q)] del self.trainedModel.data.Q, self.trainedModel.data.P self.N = N0 del self._temporaryPivot - self.scaleFactor = _scaleFactorOldPivot + self.scaleFactor = _scaleFactorOldPivot for r in req: r.wait() pMat, Ps, Qs, _, _ = gatherPivotedApproximant(pMat, Ps, Qs, self.mus.data, sizes, self.polybasis, False) self._setupTrainedModel(pMat) self.trainedModel.data.Qs, self.trainedModel.data.Ps = Qs, Ps Psupp = np.arange(0, len(self.musMarginal) * self.S, self.S) self.trainedModel.data.Psupp = list(Psupp) self._preliminaryMarginalFinalization() self._finalizeMarginalization() vbMng(self, "DEL", "Done setting up approximant.", 5) return 0 class RationalInterpolantPivotedNoMatch(RationalInterpolantPivotedBase, GenericPivotedApproximantNoMatch): """ ROM pivoted rational interpolant (without 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': kind of snapshots orthogonalization; allowed values include 0, 1/2, and 1; defaults to 1, i.e. POD; - 'scaleFactorDer': scaling factors for derivative computation; defaults to 'AUTO'; - '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'; - '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; - 'radialDirectionalWeights': radial basis weights for pivot numerator; defaults to 1; - 'radialDirectionalWeightsAdapt': bounds for adaptive rescaling of radial basis weights; defaults to [-1, -1]; - 'radialDirectionalWeightsMarginal': radial basis weights for marginal interpolant; defaults to 1; - 'functionalSolve': strategy for minimization of denominator functional; allowed values include 'NORM', 'DOMINANT', 'BARYCENTRIC[_NORM]', and 'BARYCENTRIC_AVERAGE' (check pdf in main folder for explanation); defaults to 'NORM'; - 'interpTol': tolerance for pivot interpolation; defaults to None; - 'QTol': tolerance for robust rational denominator management; defaults to 0. Defaults to empty dict. 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': kind of snapshots orthogonalization; - 'scaleFactorDer': scaling factors for derivative computation; - 'polybasis': type of polynomial basis for pivot interpolation; - 'M': degree of rational interpolant numerator; - 'N': degree of rational interpolant denominator; - 'radialDirectionalWeights': radial basis weights for pivot numerator; - 'radialDirectionalWeightsAdapt': bounds for adaptive rescaling of radial basis weights; - 'radialDirectionalWeightsMarginal': radial basis weights for marginal interpolant; - 'functionalSolve': strategy for minimization of denominator functional; - 'interpTol': tolerance for pivot interpolation; - 'QTol': 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; - 'samplerMarginal': marginal sample point generator. verbosity: Verbosity level. POD: Kind of snapshots orthogonalization. scaleFactorDer: Scaling factors for derivative computation. 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. M: Numerator degree of approximant. N: Denominator degree of approximant. radialDirectionalWeights: Radial basis weights for pivot numerator. radialDirectionalWeightsAdapt: Bounds for adaptive rescaling of radial basis weights. radialDirectionalWeightsMarginal: Radial basis weights for marginal interpolant. functionalSolve: Strategy for minimization of denominator functional. interpTol: Tolerance for pivot interpolation. QTol: Tolerance for robust rational denominator management. muBounds: list of bounds for pivot parameter values. muBoundsMarginal: list of bounds for marginal parameter values. samplingEngine: Sampling engine. uHF: High fidelity solution(s) with parameter(s) lastSolvedHF as sampleList. lastSolvedHF: Parameter(s) corresponding to last computed high fidelity solution(s) as parameterList. uApproxReduced: Reduced approximate solution(s) with parameter(s) lastSolvedApprox as sampleList. lastSolvedApproxReduced: Parameter(s) corresponding to last computed reduced approximate solution(s) as parameterList. uApprox: Approximate solution(s) with parameter(s) lastSolvedApprox as sampleList. lastSolvedApprox: Parameter(s) corresponding to last computed approximate solution(s) as parameterList. Q: Numpy 1D vector containing complex coefficients of approximant denominator. P: Numpy 2D vector whose columns are FE dofs of coefficients of approximant numerator. """ class RationalInterpolantPivotedPoleMatch(RationalInterpolantPivotedBase, GenericPivotedApproximantPoleMatch): """ 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': kind of snapshots orthogonalization; allowed values include 0, 1/2, and 1; defaults to 1, i.e. POD; - 'scaleFactorDer': scaling factors for derivative computation; defaults to 'AUTO'; - 'matchState': whether to match the system state rather than the system output; defaults to False; - 'matchingWeight': weight for pole matching optimization; defaults to 1; - 'matchingChordalRadius': radius to be used in chordal metric for poles and residues; if <= 0, Euclidean metric is used; if 'AUTO', automatically selected; defaults to -1; - 'matchingShared': required ratio of marginal points to share resonance; defaults to 1.; - 'badPoleCorrection': strategy for correction of bad poles; available values include 'ERASE', 'RATIONAL', and 'POLYNOMIAL'; defaults to 'ERASE'; - '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_*', 'LEGENDRE_*', 'NEARESTNEIGHBOR', and 'PIECEWISE_LINEAR_*'; defaults to 'MONOMIAL'; - 'paramsMarginal': dictionary of parameters for marginal interpolation; include: . 'MMarginal': degree of marginal interpolant; defaults to 'AUTO', i.e. maximum allowed; not for 'NEARESTNEIGHBOR' or 'PIECEWISE_LINEAR_*'; . 'nNeighborsMarginal': number of marginal nearest neighbors; defaults to 1; only for 'NEARESTNEIGHBOR'; . 'polydegreetypeMarginal': type of polynomial degree for marginal; defaults to 'TOTAL'; not for 'NEARESTNEIGHBOR' or 'PIECEWISE_LINEAR_*'; . 'interpTolMarginal': tolerance for marginal interpolation; defaults to None; not for 'NEARESTNEIGHBOR' or 'PIECEWISE_LINEAR_*'; . 'radialDirectionalWeightsMarginalAdapt': bounds for adaptive rescaling of marginal radial basis weights; only for radial basis. - '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; - 'radialDirectionalWeights': radial basis weights for pivot numerator; defaults to 1; - 'radialDirectionalWeightsAdapt': bounds for adaptive rescaling of radial basis weights; defaults to [-1, -1]; - 'radialDirectionalWeightsMarginal': radial basis weights for marginal interpolant; defaults to 1; - 'functionalSolve': strategy for minimization of denominator functional; allowed values include 'NORM', 'DOMINANT', 'BARYCENTRIC[_NORM]', and 'BARYCENTRIC_AVERAGE' (check pdf in main folder for explanation); defaults to 'NORM'; - 'interpTol': tolerance for pivot interpolation; defaults to None; - 'QTol': tolerance for robust rational denominator management; defaults to 0. Defaults to empty dict. 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': kind of snapshots orthogonalization; - 'scaleFactorDer': scaling factors for derivative computation; - 'matchState': whether to match the system state rather than the system output; - 'matchingWeight': weight for pole matching optimization; - 'matchingChordalRadius': radius to be used in chordal metric for poles and residues; - 'matchingShared': required ratio of marginal points to share resonance; - 'badPoleCorrection': strategy for correction of bad poles; - 'polybasis': type of polynomial basis for pivot interpolation; - 'polybasisMarginal': type of polynomial basis for marginal interpolation; - 'paramsMarginal': dictionary of parameters for marginal interpolation; include: . 'MMarginal': degree of marginal interpolant; . 'nNeighborsMarginal': number of marginal nearest neighbors; . 'polydegreetypeMarginal': type of polynomial degree for marginal; . 'interpTolMarginal': tolerance for marginal interpolation; . 'radialDirectionalWeightsMarginalAdapt': bounds for adaptive rescaling of marginal radial basis weights. - 'M': degree of rational interpolant numerator; - 'N': degree of rational interpolant denominator; - 'radialDirectionalWeights': radial basis weights for pivot numerator; - 'radialDirectionalWeightsAdapt': bounds for adaptive rescaling of radial basis weights; - 'radialDirectionalWeightsMarginal': radial basis weights for marginal interpolant; - 'functionalSolve': strategy for minimization of denominator functional; - 'interpTol': tolerance for pivot interpolation; - 'QTol': 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; - 'samplerMarginal': marginal sample point generator. verbosity: Verbosity level. POD: Kind of snapshots orthogonalization. scaleFactorDer: Scaling factors for derivative computation. matchState: Whether to match the system state rather than the system output. matchingWeight: Weight for pole matching optimization. matchingChordalRadius: Radius to be used in chordal metric for poles and residues. matchingShared: Required ratio of marginal points to share resonance. badPoleCorrection: Strategy for correction of bad poles. S: Total number of pivot samples current approximant relies upon. samplerPivot: Pivot sample point generator. SMarginal: Total number of marginal samples current approximant relies upon. samplerMarginal: Marginal sample point generator. polybasis: Type of polynomial basis for pivot interpolation. polybasisMarginal: Type of polynomial basis for marginal interpolation. paramsMarginal: Dictionary of parameters for marginal interpolation. M: Numerator degree of approximant. N: Denominator degree of approximant. radialDirectionalWeights: Radial basis weights for pivot numerator. radialDirectionalWeightsAdapt: Bounds for adaptive rescaling of radial basis weights. radialDirectionalWeightsMarginal: Radial basis weights for marginal interpolant. functionalSolve: Strategy for minimization of denominator functional. interpTol: Tolerance for pivot interpolation. QTol: Tolerance for robust rational denominator management. muBounds: list of bounds for pivot parameter values. muBoundsMarginal: list of bounds for marginal parameter values. samplingEngine: Sampling engine. uHF: High fidelity solution(s) with parameter(s) lastSolvedHF as sampleList. lastSolvedHF: Parameter(s) corresponding to last computed high fidelity solution(s) as parameterList. uApproxReduced: Reduced approximate solution(s) with parameter(s) lastSolvedApprox as sampleList. lastSolvedApproxReduced: Parameter(s) corresponding to last computed reduced approximate solution(s) as parameterList. uApprox: Approximate solution(s) with parameter(s) lastSolvedApprox as sampleList. lastSolvedApprox: Parameter(s) corresponding to last computed approximate solution(s) as parameterList. Q: Numpy 1D vector containing complex coefficients of approximant denominator. P: Numpy 2D vector whose columns are FE dofs of coefficients of approximant numerator. """ def setupApprox(self, *args, **kwargs) -> int: if self.checkComputedApprox(): return -1 self.purgeparamsMarginal() setupOK = super().setupApprox(*args, **kwargs) if self.matchState: self._postApplyC() return setupOK diff --git a/rrompy/reduction_methods/pivoted/trained_model/trained_model_pivoted_rational_nomatch.py b/rrompy/reduction_methods/pivoted/trained_model/trained_model_pivoted_rational_nomatch.py index 5ea1d21..cc18e84 100644 --- a/rrompy/reduction_methods/pivoted/trained_model/trained_model_pivoted_rational_nomatch.py +++ b/rrompy/reduction_methods/pivoted/trained_model/trained_model_pivoted_rational_nomatch.py @@ -1,229 +1,232 @@ # Copyright (C) 2018-2020 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 rrompy.reduction_methods.standard.trained_model.trained_model_rational \ import TrainedModelRational from rrompy.utilities.base.types import (Np1D, Np2D, Tuple, List, ListAny, paramVal, paramList, sampList) from rrompy.utilities.base import verbosityManager as vbMng, freepar as fp from rrompy.utilities.numerical import dot from rrompy.utilities.numerical.compress_matrix import compressMatrix from rrompy.utilities.poly_fitting.heaviside import rational2heaviside from rrompy.utilities.poly_fitting.nearest_neighbor import ( NearestNeighborInterpolator as NNI) from rrompy.utilities.exception_manager import (RROMPyException, RROMPyAssert, RROMPyWarning) from rrompy.parameter import checkParameterList from rrompy.sampling import emptySampleList __all__ = ['TrainedModelPivotedRationalNoMatch'] class TrainedModelPivotedRationalNoMatch(TrainedModelRational): """ ROM approximant evaluation for pivoted approximants based on interpolation of rational approximants (without pole matching). Attributes: Data: dictionary with all that can be pickled. """ def checkParameterListPivot(self, mu:paramList, check_if_single : bool = False) -> paramList: return checkParameterList(mu, self.data.nparPivot, check_if_single) def checkParameterListMarginal(self, mu:paramList, check_if_single : bool = False) -> paramList: return checkParameterList(mu, self.data.nparMarginal, check_if_single) def compress(self, collapse : bool = False, tol : float = 0., returnRMat : bool = False, **compressMatrixkwargs): if not collapse and tol <= 0.: return RMat = self.data.projMat if not collapse: if hasattr(self.data, "_compressTol"): RROMPyWarning(("Recompressing already compressed model is " "ineffective. Aborting.")) return self.data.projMat, RMat, _ = compressMatrix(RMat, tol, **compressMatrixkwargs) if hasattr(self.data, "Ps"): for obj, suppj in zip(self.data.Ps, self.data.Psupp): obj.postmultiplyTensorize(RMat.T[suppj : suppj + obj.shape[0]]) if hasattr(self, "_PsExcl"): for obj, suppj in zip(self._PsExcl, self._PsuppExcl): obj.postmultiplyTensorize(RMat.T[suppj : suppj + obj.shape[0]]) self._PsuppExcl = [0] * len(self._PsuppExcl) self.data.Psupp = [0] * len(self.data.Psupp) super(TrainedModelRational, self).compress(collapse, tol) if returnRMat: return RMat 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 = self.checkParameterListPivot(mu) if mu0 is None: mu0 = self.checkParameterListPivot( self.data.mu0(0, self.data.directionPivot)) return (self.mapParameterList(mu, idx = self.data.directionPivot) - self.mapParameterList(mu0, idx = self.data.directionPivot) ) / [self.data.scaleFactor[x] for x in self.data.directionPivot] def setupMarginalInterp(self, interpPars:ListAny): self.data.marginalInterp = NNI() self.data.marginalInterp.setupByInterpolation(self.data.musMarginal, np.arange(len(self.data.musMarginal)), 1, *interpPars) def updateEffectiveSamples(self, exclude:List[int]): if hasattr(self, "_idxExcl"): for j, excl in enumerate(self._idxExcl): self.data.musMarginal.insert(self._musMExcl[j], excl) self.data.Ps.insert(excl, self._PsExcl[j]) self.data.Qs.insert(excl, self._QsExcl[j]) self.data.Psupp.insert(excl, self._PsuppExcl[j]) self._idxExcl, self._musMExcl = list(np.sort(exclude)), [] self._PsExcl, self._QsExcl, self._PsuppExcl = [], [], [] for excl in self._idxExcl[::-1]: self._musMExcl = [self.data.musMarginal[excl]] + self._musMExcl self.data.musMarginal.pop(excl) self._PsExcl = [self.data.Ps.pop(excl)] + self._PsExcl self._QsExcl = [self.data.Qs.pop(excl)] + self._QsExcl self._PsuppExcl = [self.data.Psupp.pop(excl)] + self._PsuppExcl 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 = self.checkParameterList(mu) muP = self.centerNormalizePivot(mu(self.data.directionPivot)) muM = self.checkParameterListMarginal(mu(self.data.directionMarginal)) idxMUnique, idxMmap = np.unique(self.data.marginalInterp(muM), return_inverse = True) idxMUnique = np.array(idxMUnique, dtype = int) p = emptySampleList() vbMng(self, "INIT", "Evaluating numerator at mu = {}.".format(mu), 17) for i, iM in enumerate(idxMUnique): idx = np.where(idxMmap == i)[0] Pval, supp = self.data.Ps[iM](muP[idx]), self.data.Psupp[iM] if i == 0: - p.reset((self.data.projMat.shape[1], len(mu)), - dtype = Pval.dtype) + if hasattr(self.data.projMat, "shape"): + plen = self.data.projMat.shape[1] + else: + plen = len(Pval) + p.reset((plen, len(mu)), dtype = Pval.dtype) p.data[:] = 0. p.data[supp : supp + len(Pval), idx] = Pval vbMng(self, "DEL", "Done evaluating numerator.", 17) 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 = self.checkParameterList(mu) muP = self.centerNormalizePivot(mu(self.data.directionPivot)) muM = self.checkParameterListMarginal(mu(self.data.directionMarginal)) 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]] idxMUnique, idxMmap = np.unique(self.data.marginalInterp(muM), return_inverse = True) idxMUnique = np.array(idxMUnique, dtype = int) vbMng(self, "INIT", "Evaluating denominator at mu = {}.".format(mu), 17) for i, iM in enumerate(idxMUnique): idx = np.where(idxMmap == i)[0] Qval = self.data.Qs[iM](muP[idx], derP, sclP) if i == 0: q = np.empty(len(mu), dtype = Qval.dtype) q[idx] = Qval vbMng(self, "DEL", "Done evaluating denominator.", 17) return q def getPoles(self, marginalVals : ListAny = [fp]) -> paramList: """ Obtain approximant poles. Returns: Numpy complex vector of poles. """ RROMPyAssert(self.data.nparPivot, 1, "Number of pivot parameters") mVals = list(marginalVals) rDim = mVals.index(fp) if rDim < len(mVals) - 1 and fp in mVals[rDim + 1 :]: 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)[0] muM = self.checkParameterListMarginal([mVals[j] for j in range(len(mVals)) if j != rDim]) iM = int(self.data.marginalInterp(muM)) roots = self.data.scaleFactor[rDim] * self.data.Qs[iM].roots() return self.mapParameterList(self.mapParameterList(self.data.mu0(rDim), idx = [rDim])(0, 0) + roots, "B", [rDim])(0) def getResidues(self, *args, **kwargs) -> Tuple[paramList, Np2D]: """ Obtain approximant residues. Returns: Numpy matrix with residues as columns. """ pls = self.getPoles(*args, **kwargs) if len(pls) == 0: return pls, np.empty((0, 0), dtype = self.data.Ps[0].coeffs.dtype) if len(args) == 1: mVals = args[0] elif len(args) == 0: mVals = [fp] else: mVals = kwargs["marginalVals"] rDim = mVals.index(fp) mVals[rDim] = self.data.mu0(rDim)[0] muM = self.checkParameterListMarginal([mVals[j] for j in range(len(mVals)) if j != rDim]) iM = int(self.data.marginalInterp(muM)) res = rational2heaviside(self.data.Ps[iM], self.data.Qs[iM])[0] res = res[: len(pls), :].T if not self.data._collapsed: res = dot(self.data.projMat, res).T return pls, res diff --git a/rrompy/reduction_methods/pivoted/trained_model/trained_model_pivoted_rational_polematch.py b/rrompy/reduction_methods/pivoted/trained_model/trained_model_pivoted_rational_polematch.py index 8e24cbd..91dd9d0 100644 --- a/rrompy/reduction_methods/pivoted/trained_model/trained_model_pivoted_rational_polematch.py +++ b/rrompy/reduction_methods/pivoted/trained_model/trained_model_pivoted_rational_polematch.py @@ -1,594 +1,658 @@ # Copyright (C) 2018-2020 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 warnings import numpy as np from scipy.special import factorial as fact from scipy.sparse import csr_matrix, hstack, SparseEfficiencyWarning from collections.abc import Iterable from copy import deepcopy as copy from itertools import combinations from .trained_model_pivoted_rational_nomatch import ( TrainedModelPivotedRationalNoMatch) from rrompy.utilities.base.types import (Tuple, Np1D, Np2D, List, ListAny, paramVal, paramList, sampList, HFEng) from rrompy.utilities.base import verbosityManager as vbMng, freepar as fp from rrompy.utilities.numerical import dot from rrompy.utilities.numerical.point_matching import rationalFunctionMatching from rrompy.utilities.numerical.degree import 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.heaviside import (rational2heaviside, polyval as heavival, heavisideUniformShape, HeavisideInterpolator as HI) from rrompy.utilities.poly_fitting.nearest_neighbor import ( NearestNeighborInterpolator as NNI) from rrompy.utilities.poly_fitting.piecewise_linear import (sparsekinds, PiecewiseLinearInterpolator as PLI) from rrompy.utilities.exception_manager import (RROMPyException, RROMPyAssert, RROMPyWarning) from rrompy.sampling import sampleList, emptySampleList __all__ = ['TrainedModelPivotedRationalPoleMatch'] def getChordalScaling(x:Np2D, r:float, projGramian : Np2D = 1., projHalfGramian : Np2D = None) -> Tuple[Np2D, Np2D]: goodX = np.where(np.isinf(x[:, 0]) == False)[0] normX = np.empty(len(x)) if projGramian is None: normX[goodX] = np.sum(np.abs(dot(projHalfGramian, x[goodX].T)) ** 2., axis = 0) ** .5 else: normX[goodX] = np.abs(np.sum(dot(projGramian, x[goodX].T) * x[goodX].T.conj(), axis = 0)) ** .5 scale = np.ones((len(normX), 1)) scale[goodX, 0] = 1. / ((normX[goodX] / r) ** 2. + 1.) xscaled = np.zeros_like(x) for j in goodX: xscaled[j] = x[j] * scale[j] return xscaled, r * (1 - scale) def normalizeChordal(x:Np2D, r:float, projGramian : Np2D = 1., projHalfGramian : Np2D = None) -> Np2D: for j in range(x.shape[0]): x[j, -1] -= .5 * r if projGramian is None: norm2xj = np.sum(np.abs(dot(projHalfGramian, x[j, : -1])) ** 2.) else: norm2xj = np.abs(np.sum(dot(projGramian, x[j, : -1]) * x[j, : -1].conj())) normxj = (norm2xj + np.abs(x[j, -1]) ** 2.) ** .5 if normxj < 1e-15: normxj += np.finfo(float).eps x[j] *= .5 * r / normxj x[j, -1] += .5 * r return x def pullbackChordal(x:Np2D, r:float) -> Np2D: y = copy(x[:, : -1]) for j, p in enumerate(x[:, -1]): scalexj = 1. - p / r y[j] = np.inf if scalexj < 1e-15 else y[j] / scalexj return y class TrainedModelPivotedRationalPoleMatch(TrainedModelPivotedRationalNoMatch): """ ROM approximant evaluation for pivoted approximants based on interpolation of rational approximants (with pole matching). Attributes: Data: dictionary with all that can be pickled. """ def compress(self, collapse : bool = False, tol : float = 0., returnRMat : bool = False, **compressMatrixkwargs): Psupp = copy(self.data.Psupp) RMat = super().compress(collapse, tol, True, **compressMatrixkwargs) if RMat is None: return + for j in range(len(self.data.coeffsEff)): + self.data.coeffsEff[j] = dot(self.data.coeffsEff[j], RMat.T) for obj, suppj in zip(self.data.HIs, Psupp): obj.postmultiplyTensorize(RMat.T[suppj : suppj + obj.shape[0]]) if hasattr(self, "_HIsExcl"): for obj, suppj in zip(self._HIsExcl, Psupp): obj.postmultiplyTensorize(RMat.T[suppj : suppj + obj.shape[0]]) if not hasattr(self, "_PsExcl"): self._PsuppExcl = [0] * len(self._PsuppExcl) if returnRMat: return RMat 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 = self.checkParameterListMarginal(mu) if mu0 is None: mu0 = self.checkParameterListMarginal( self.data.mu0(0, self.data.directionMarginal)) return (self.mapParameterList(mu, idx = self.data.directionMarginal) - self.mapParameterList(mu0, idx = self.data.directionMarginal) ) / [self.data.scaleFactor[x] for x in self.data.directionMarginal] def setupMarginalInterp(self, approx, interpPars:ListAny, extraPar = None): vbMng(self, "INIT", "Starting computation of marginal interpolator.", 12) musMCN = self.centerNormalizeMarginal(self.data.musMarginal) nM, pbM = len(musMCN), approx.polybasisMarginal if pbM in ppb + rbpb: if extraPar: approx._setMMarginalAuto() _MMarginalEff = approx.paramsMarginal["MMarginal"] if pbM in ppb: p = PI() elif pbM in rbpb: p = RBI() else: # if pbM in sparsekinds + ["NEARESTNEIGHBOR"]: if pbM == "NEARESTNEIGHBOR": p = NNI() else: # if pbM in sparsekinds: pllims = [[-1.] * self.data.nparMarginal, [1.] * self.data.nparMarginal] p = PLI() for ipts, pts in enumerate(self.data.suppEffPts): if len(pts) == 0: raise RROMPyException("Empty list of support points.") musMCNEff, valsEff = musMCN[pts], np.eye(len(pts)) if pbM in ppb + rbpb: if extraPar: if ipts > 0: verb = approx.verbosity approx.verbosity = 0 _musM = approx.musMarginal approx.musMarginal = musMCNEff approx._setMMarginalAuto() approx.musMarginal = _musM approx.verbosity = verb else: approx.paramsMarginal["MMarginal"] = reduceDegreeN( _MMarginalEff, len(musMCNEff), self.data.nparMarginal, approx.paramsMarginal["polydegreetypeMarginal"]) MMEff = approx.paramsMarginal["MMarginal"] while MMEff >= 0: wellCond, msg = p.setupByInterpolation(musMCNEff, valsEff, MMEff, *interpPars) vbMng(self, "MAIN", msg, 30) if wellCond: break vbMng(self, "MAIN", ("Polyfit is poorly conditioned. Reducing " "MMarginal by 1."), 35) MMEff -= 1 if MMEff < 0: raise RROMPyException(("Instability in computation of " "interpolant. Aborting.")) if (pbM in rbpb and len(interpPars) > 4 and "optimizeScalingBounds" in interpPars[4].keys()): interpPars[4]["optimizeScalingBounds"] = [-1., -1.] elif pbM == "NEARESTNEIGHBOR": if ipts > 0: interpPars[0] = 1 p.setupByInterpolation(musMCNEff, valsEff, *interpPars) elif ipts == 0: # and pbM in sparsekinds: p.setupByInterpolation(musMCNEff, valsEff, pllims, extraPar[pts], *interpPars) if ipts == 0: self.data.marginalInterp = copy(p) self.data.coeffsEff, self.data.polesEff = [], [] N = len(self.data.suppEffIdx) goodIdx = np.where(self.data.suppEffIdx != -1)[0] for hi, sup in zip(self.data.HIs, self.data.Psupp): pEff, cEff = hi.poles.reshape(-1, 1), hi.coeffs if self.data.chordalRadius[0] > 0.: pEff = np.hstack(getChordalScaling(pEff, self.data.chordalRadius[0])) if self.data.chordalRadius[1] > 0.: if self.data.projGramian is None: projGramian = None projHalfGramian = self.data.projMat[:, sup : sup + cEff.shape[1]] else: projGramian = self.data.projGramian[ sup : sup + cEff.shape[1]][:, sup : sup + cEff.shape[1]] projHalfGramian = None cEff, cEffH = getChordalScaling(cEff, self.data.chordalRadius[1], projGramian, projHalfGramian) else: cEffH = np.empty((cEff.shape[0], 0)) if (self.data._collapsed or self.data.projMat.shape[1] == cEff.shape[1]): cEff = np.hstack([cEff, cEffH]) else: supC = self.data.projMat.shape[1] - sup - cEff.shape[1] cEff = hstack((csr_matrix((len(cEff), sup)), csr_matrix(cEff), csr_matrix((len(cEff), supC)), cEffH), "csr") goodIdxC = np.append(goodIdx, np.arange(N, cEff.shape[0])) self.data.coeffsEff += [cEff[goodIdxC, :]] self.data.polesEff += [pEff[goodIdx]] else: ptsBad = [i for i in range(nM) if i not in pts] idxBad = np.where(self.data.suppEffIdx[goodIdx] == ipts)[0] warnings.simplefilter('ignore', SparseEfficiencyWarning) if pbM in sparsekinds: for ij, j in enumerate(ptsBad): nearest = pts[np.argmin(np.sum(np.abs(musMCNEff.data - np.tile(musMCN[j], [len(pts), 1]) ), axis = 1).flatten())] self.data.coeffsEff[j][idxBad] = copy( self.data.coeffsEff[nearest][idxBad]) self.data.polesEff[j][idxBad] = copy( self.data.polesEff[nearest][idxBad]) else: if (self.data._collapsed or self.data.projMat.shape[1] == cEff.shape[1]): cfBase = np.zeros((len(idxBad), cEff.shape[1]), dtype = cEff.dtype) else: cfBase = csr_matrix((len(idxBad), self.data.coeffsEff[0].shape[1]), dtype = cEff.dtype) valMuMBad = p(musMCN[ptsBad]) for ijb, jb in enumerate(ptsBad): self.data.coeffsEff[jb][idxBad] = copy(cfBase) self.data.polesEff[jb][idxBad] = 0. for ij, j in enumerate(pts): val = valMuMBad[ij][ijb] if not np.isclose(val, 0., atol = 1e-15): self.data.coeffsEff[jb][idxBad] += (val * self.data.coeffsEff[j][idxBad]) self.data.polesEff[jb][idxBad] += (val * self.data.polesEff[j][idxBad]) if self.data.chordalRadius[0] > 0: self.data.polesEff[jb][idxBad] = normalizeChordal( self.data.polesEff[jb][idxBad], self.data.chordalRadius[0]) if self.data.chordalRadius[1] > 0: self.data.coeffsEff[jb][idxBad] = normalizeChordal( self.data.coeffsEff[jb][idxBad], self.data.chordalRadius[1], self.data.projGramian, self.data.projMat) warnings.filters.pop(0) if pbM in ppb + rbpb: approx.paramsMarginal["MMarginal"] = _MMarginalEff vbMng(self, "DEL", "Done computing marginal interpolator.", 12) def updateEffectiveSamples(self, exclude:List[int], *args, **kwargs): if hasattr(self, "_idxExcl"): for j, excl in enumerate(self._idxExcl): self.data.HIs.insert(excl, self._HIsExcl[j]) super().updateEffectiveSamples(exclude) self._HIsExcl = [] for excl in self._idxExcl[::-1]: self._HIsExcl = [self.data.HIs.pop(excl)] + self._HIsExcl poles = [hi.poles for hi in self.data.HIs] coeffs = [hi.coeffs for hi in self.data.HIs] self.initializeFromLists(poles, coeffs, self.data.Psupp, self.data.HIs[0].polybasis, *args, **kwargs) def initializeFromRational(self, *args, **kwargs): """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, self.data.Psupp, basis, *args, **kwargs) def initializeFromLists(self, poles:ListAny, coeffs:ListAny, supps:ListAny, basis:str, matchingWeight:float, HFEngine:HFEng, is_state:bool, chordalRadius:Tuple[float, float]): """Initialize Heaviside representation.""" poles, coeffs = heavisideUniformShape(poles, coeffs) N = len(poles[0]) if chordalRadius[0] == "AUTO": chordalRadius[0] = 1. if chordalRadius[1] == "AUTO": norm2s = 0. for c, sup in zip(coeffs, self.data.Psupp): if self.data.projGramian is None: gramEff = self.data.projMat[:, sup : sup + c.shape[1]] norm2s += np.sum(np.abs(dot(gramEff, c[: N].T)) ** 2.) else: gramEff = self.data.projGramian[sup : sup + c.shape[1]][:, sup : sup + c.shape[1]] norm2s += np.sum(np.abs(dot(gramEff, c[: N].T) * c[: N].T.conj())) chordalRadius[1] = (norm2s / N / len(coeffs)) ** .5 self.data.chordalRadius = copy(chordalRadius) if is_state and chordalRadius[1] > 0: chordalRadius[1] = "AUTO" poles, coeffs = rationalFunctionMatching(poles, coeffs, self.data.musMarginal.data, matchingWeight, supps, self.data.projMat, HFEngine, is_state, None, chordalRadius) self.data.HIs = [] for pls, cfs in zip(poles, coeffs): hsi = HI() hsi.poles = pls if len(cfs) == len(pls): cfs = np.pad(cfs, ((0, 1), (0, 0)), "constant") hsi.coeffs = cfs hsi.npar = 1 hsi.polybasis = basis self.data.HIs += [hsi] self.data.suppEffPts = [np.arange(len(self.data.HIs))] self.data.suppEffIdx = np.zeros(len(poles[0]), dtype = int) def checkShared(self, shared:float, correction : str = "ERASE") -> str: N = len(self.data.HIs[0].poles) M = len(self.data.HIs) correction = correction.upper().strip().replace(" ","") if correction not in ["ERASE", "RATIONAL", "POLYNOMIAL"]: RROMPyWarning(("Correction kind not recognized. Overriding to " "'ERASE'.")) correction = "ERASE" goodLocPoles = np.array([np.isinf(hi.poles) == False for hi in self.data.HIs]) self.data.suppEffPts = [np.arange(len(self.data.HIs))] self.data.suppEffIdx = - np.ones(N, dtype = int) goodGlobPoles = np.sum(goodLocPoles, axis = 0) goodEnoughPoles = goodGlobPoles >= max(1., 1. * shared * M) keepPole = np.where(goodEnoughPoles)[0] halfPole = np.where(goodEnoughPoles * (goodGlobPoles < M))[0] self.data.suppEffIdx[keepPole] = 0 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 degBad = len(self.data.HIs[0].coeffs) - N - 1 for pt in range(len(self.data.HIs)): idxR = np.where(goodLocPoles[pt] * (goodEnoughPoles == False))[0] - if len(idxR) == 0: continue - if correction != "ERASE": - muM, musEff = self.data.musMarginal[pt], [] - for mu in self.data.mus: - if np.allclose(mu(self.data.directionMarginal), muM): - musEff += [mu(self.data.directionPivot[0])] - musEff = self.centerNormalizePivot(musEff) - idxK = np.where(goodLocPoles[pt] * goodEnoughPoles)[0] - plsBad = self.data.HIs[pt].poles[idxR] - plsGood = self.data.HIs[pt].poles[idxK] - corrVals = heavival(musEff, self.data.HIs[pt].coeffs[idxR], - plsBad, self.data.HIs[pt].polybasis).T - if correction == "RATIONAL": - hi = HI() - hi.setupByInterpolation(musEff, plsGood, corrVals, degBad, - self.data.HIs[pt].polybasis) - self.data.HIs[pt].coeffs[idxK] += hi.coeffs[: len(idxK)] - polyCorr = hi.coeffs[len(idxK) :] - elif correction == "POLYNOMIAL": - pi = PI() - pi.setupByInterpolation(musEff, corrVals, degBad, - self.data.HIs[pt].polybasis.split("_")[0]) - polyCorr = pi.coeffs - self.data.HIs[pt].coeffs[N : N + degBad + 1] += polyCorr - self.data.HIs[pt].poles[idxR] = np.inf - self.data.HIs[pt].coeffs[idxR] = 0. + self.removePoleResLocal(idxR, pt, degBad, correction, True) return ("Hard-erased {} pole".format(N - len(keepPole)) + "s" * (N - len(keepPole) != 1) + " and soft-erased {} pole".format(len(halfPole)) + "s" * (len(halfPole) != 1) + ".") + def removePoleResLocal(self, badidx:List[int], margidx:int, + degcorr : int = None, correction : str = "ERASE", + hidden : bool = False): + if not hasattr(badidx, "__len__"): badidx = [badidx] + badidx = np.array(badidx) + if len(badidx) == 0: return + correction = correction.upper().strip().replace(" ","") + if correction not in ["ERASE", "RATIONAL", "POLYNOMIAL"]: + RROMPyWarning(("Correction kind not recognized. Overriding to " + "'ERASE'.")) + correction = "ERASE" + if hidden: + N = len(self.data.HIs[margidx].poles) + else: + N = len(self.data.polesEff[margidx]) + goodidx = [j for j in range(N) if j not in badidx] + if correction != "ERASE": + if degcorr is None: + if hidden: + degcorr = len(self.data.HIs[margidx].coeffs) - N - 1 + else: + degcorr = self.data.coeffsEff[margidx].shape[0] - N - 1 + muM, musEff = self.data.musMarginal[margidx], [] + polybasis = self.data.HIs[margidx].polybasis + for mu in self.data.mus: + if np.allclose(mu(self.data.directionMarginal), muM): + musEff += [mu(self.data.directionPivot[0])] + musEff = self.centerNormalizePivot(musEff) + if hidden: + plsBad = self.data.HIs[margidx].poles[badidx] + else: + plsBad = self.data.polesEff[margidx][badidx, 0] + plsBadEff = np.isinf(plsBad) == False + plsBad, badidx = plsBad[plsBadEff], badidx[plsBadEff] + if hidden: + plsGood = self.data.HIs[margidx].poles[goodidx] + corrVals = heavival(musEff, + self.data.HIs[margidx].coeffs[badidx], + plsBad, polybasis).T + else: + plsGood = self.data.polesEff[margidx][goodidx] + corrVals = heavival(musEff, + self.data.coeffsEff[margidx].toarray()[badidx], + plsBad, polybasis).T + if correction == "RATIONAL": + hi = HI() + hi.setupByInterpolation(musEff, plsGood, corrVals, degcorr, + polybasis) + if hidden: + self.data.HIs[margidx].coeffs[goodidx] += ( + hi.coeffs[: len(goodidx)]) + else: + self.data.coeffsEff[margidx][goodidx, :] += ( + hi.coeffs[: len(goodidx)]) + polyCorr = hi.coeffs[len(goodidx) :] + elif correction == "POLYNOMIAL": + pi = PI() + pi.setupByInterpolation(musEff, corrVals, degcorr, + polybasis.split("_")[0]) + polyCorr = pi.coeffs + if hidden: + self.data.HIs[margidx].coeffs[N : N + degcorr + 1] += polyCorr + else: + self.data.coeffsEff[margidx][N : N + degcorr + 1, :] += ( + polyCorr) + if hidden: + self.data.HIs[margidx].poles[badidx] = np.inf + self.data.HIs[margidx].coeffs[badidx] = 0. + else: + self.data.polesEff[margidx] = self.data.polesEff[margidx][goodidx] + goodidx += list(range(N, self.data.coeffsEff[margidx].shape[0])) + self.data.coeffsEff[margidx] = ( + self.data.coeffsEff[margidx][goodidx, :]) + + def removePoleResGlobal(self, badidx:List[int], degcorr : int = None, + correction : str = "ERASE", hidden : bool = False): + if not hasattr(badidx, "__len__"): badidx = [badidx] + if len(badidx) == 0: return + correction = correction.upper().strip().replace(" ","") + if correction not in ["ERASE", "RATIONAL", "POLYNOMIAL"]: + RROMPyWarning(("Correction kind not recognized. Overriding to " + "'ERASE'.")) + correction = "ERASE" + for margidx in range(len(self.data.HIs)): + self.removePoleResLocal(badidx, margidx, degcorr, correction, + hidden) + 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 = self.checkParameterList(mu) if (not hasattr(self, "lastSolvedApproxReduced") or self.lastSolvedApproxReduced != mu): vbMng(self, "INIT", "Evaluating approximant at mu = {}.".format(mu), 12) muP = self.centerNormalizePivot(mu(self.data.directionPivot)) muM = mu(self.data.directionMarginal) his = self.interpolateMarginalInterpolator(muM) for i, (mP, hi) in enumerate(zip(muP, his)): uAppR = hi(mP)[:, 0] if i == 0: uApproxR = np.empty((len(uAppR), len(mu)), dtype = uAppR.dtype) uApproxR[:, i] = uAppR self.uApproxReduced = sampleList(uApproxR) vbMng(self, "DEL", "Done evaluating approximant.", 12) self.lastSolvedApproxReduced = mu return self.uApproxReduced def interpolateMarginalInterpolator(self, mu : paramList = []) -> ListAny: """Obtain interpolated approximant interpolator.""" mu = self.checkParameterListMarginal(mu) vbMng(self, "INIT", "Interpolating marginal models at mu = {}.".format(mu), 95) his = [] muC = self.centerNormalizeMarginal(mu) mIvals = self.data.marginalInterp(muC) verb, self.verbosity = self.verbosity, 0 poless = self.interpolateMarginalPoles(mu, mIvals) coeffss = self.interpolateMarginalCoeffs(mu, mIvals) self.verbosity = verb for j in range(len(mu)): his += [HI()] his[-1].poles = poless[j] his[-1].coeffs = coeffss[j] his[-1].npar = 1 his[-1].polybasis = self.data.HIs[0].polybasis vbMng(self, "DEL", "Done interpolating marginal models.", 95) return his def interpolateMarginalPoles(self, mu : paramList = [], mIvals : Np2D = None) -> ListAny: """Obtain interpolated approximant poles.""" mu = self.checkParameterListMarginal(mu) vbMng(self, "INIT", "Interpolating marginal poles at mu = {}.".format(mu), 95) intMPoles = np.zeros((len(mu),) + self.data.polesEff[0].shape, dtype = self.data.polesEff[0].dtype) if mIvals is None: muC = self.centerNormalizeMarginal(mu) mIvals = self.data.marginalInterp(muC) for pEff, mI in zip(self.data.polesEff, mIvals): for j, m in enumerate(mI): intMPoles[j] += m * pEff rCP = self.data.chordalRadius[0] if rCP > 0: for j in range(len(mu)): intMPoles[j, ..., 0] = pullbackChordal( normalizeChordal(intMPoles[j], rCP), rCP)[..., 0] vbMng(self, "DEL", "Done interpolating marginal poles.", 95) return intMPoles[..., 0] def interpolateMarginalCoeffs(self, mu : paramList = [], mIvals : Np2D = None) -> ListAny: """Obtain interpolated approximant coefficients.""" mu = self.checkParameterListMarginal(mu) vbMng(self, "INIT", "Interpolating marginal coefficients at mu = {}.".format(mu), 95) intMCoeffs = np.zeros((len(mu),) + self.data.coeffsEff[0].shape, dtype = self.data.coeffsEff[0].dtype) if mIvals is None: muC = self.centerNormalizeMarginal(mu) mIvals = self.data.marginalInterp(muC) for cEff, mI in zip(self.data.coeffsEff, mIvals): for j, m in enumerate(mI): intMCoeffs[j] += m * cEff rCC = self.data.chordalRadius[1] if rCC > 0: for j in range(len(mu)): intMCoeffs[j, ..., : -1] = pullbackChordal( normalizeChordal(intMCoeffs[j], rCC, self.data.projGramian, self.data.projMat), rCC) intMCoeffs = intMCoeffs[..., : -1] vbMng(self, "DEL", "Done interpolating marginal coefficients.", 95) return intMCoeffs 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 = self.checkParameterList(mu) p = emptySampleList() muP = self.centerNormalizePivot(mu(self.data.directionPivot)) muM = mu(self.data.directionMarginal) his = self.interpolateMarginalInterpolator(muM) for i, (mP, hi) in enumerate(zip(muP, his)): Pval = hi(mP) * np.prod(mP[0] - hi.poles) if i == 0: p.reset((len(Pval), len(mu)), dtype = Pval.dtype) p[i] = Pval 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 = self.checkParameterList(mu) muP = self.centerNormalizePivot(mu(self.data.directionPivot)) muM = mu(self.data.directionMarginal) 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) pls = self.interpolateMarginalPoles(muM) for i, (mP, pl) in enumerate(zip(muP, pls)): N = len(pl) if derP == N: derVal[i] = 1. elif derP >= 0 and derP < N: plDist = mP[0] - pl for terms in combinations(np.arange(N), N - derP): derVal[i] += np.prod(plDist[list(terms)]) return sclP ** derP * fact(derP) * derVal def getPoles(self, marginalVals : ListAny = [fp]) -> paramList: """ Obtain approximant poles. Returns: Numpy complex vector of poles. """ RROMPyAssert(self.data.nparPivot, 1, "Number of pivot parameters") mVals = list(marginalVals) rDim = mVals.index(fp) if rDim < len(mVals) - 1 and fp in mVals[rDim + 1 :]: 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)[0] mMarg = [mVals[j] for j in range(len(mVals)) if j != rDim] roots = (self.data.scaleFactor[rDim] * self.interpolateMarginalPoles(mMarg)[0]) return self.mapParameterList(self.mapParameterList(self.data.mu0(rDim), idx = [rDim])(0, 0) + roots, "B", [rDim])(0) def getResidues(self, *args, **kwargs) -> Tuple[paramList, Np2D]: """ 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 isinstance(mVals, Iterable): mVals = [mVals] mVals = list(mVals) rDim = mVals.index(fp) mMarg = [mVals[j] for j in range(len(mVals)) if j != rDim] res = self.interpolateMarginalCoeffs(mMarg)[0][: len(pls), :].T if not self.data._collapsed: res = dot(self.data.projMat, res).T return pls, res diff --git a/rrompy/reduction_methods/standard/trained_model/trained_model_nearest_neighbor.py b/rrompy/reduction_methods/standard/trained_model/trained_model_nearest_neighbor.py index 4a894ed..4c2c141 100644 --- a/rrompy/reduction_methods/standard/trained_model/trained_model_nearest_neighbor.py +++ b/rrompy/reduction_methods/standard/trained_model/trained_model_nearest_neighbor.py @@ -1,94 +1,97 @@ # Copyright (C) 2018-2020 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 rrompy.reduction_methods.base.trained_model.trained_model import ( TrainedModel) from rrompy.utilities.numerical.compress_matrix import compressMatrix from rrompy.utilities.base.types import Np1D, paramList, sampList from rrompy.utilities.base import verbosityManager as vbMng from rrompy.utilities.exception_manager import RROMPyWarning from rrompy.sampling import emptySampleList __all__ = ['TrainedModelNearestNeighbor'] class TrainedModelNearestNeighbor(TrainedModel): """ ROM approximant evaluation for nearest neighbor approximant. Attributes: Data: dictionary with all that can be pickled. """ def compress(self, collapse : bool = False, tol : float = 0., *args, **kwargs): if not collapse and tol <= 0.: return RMat = self.data.projMat if not collapse: if hasattr(self.data, "_compressTol"): RROMPyWarning(("Recompressing already compressed model is " "ineffective. Aborting.")) return self.data.projMat, RMat, _ = compressMatrix(RMat, tol, *args, **kwargs) for j, suppj in enumerate(self.data.supp): self.data.vals[j] = RMat[:, suppj : suppj + len(self.data.vals[j]) ].dot(self.data.vals[j]) self.data.supp[j] = [0] super().compress(collapse, tol) def getNearestNeighbor(self, mu : paramList = []) -> Np1D: """ Find nearest neighbor to arbitrary parameter. Args: mu: Target parameter. """ mu = self.checkParameterList(mu) idxUnique, idxmap = np.unique(self.data.NN(mu), return_inverse = True) idxUnique = np.array(idxUnique, dtype = int) vbMng(self, "INIT", "Finding nearest neighbor to mu = {}.".format(mu), 22) nn = emptySampleList() for i, iM in enumerate(idxUnique): idx = np.where(idxmap == i)[0] val, supp = self.data.vals[iM], self.data.supp[iM] if i == 0: - nn.reset((self.data.projMat.shape[1], len(mu)), - dtype = val.dtype) + if hasattr(self.data.projMat, "shape"): + nnlen = self.data.projMat.shape[1] + else: + nnlen = len(val) + nn.reset((nnlen, len(mu)), dtype = val.dtype) nn.data[:] = 0. for i in idx: nn.data[supp : supp + len(val), i] = val vbMng(self, "DEL", "Done finding nearest neighbor.", 22) return nn def getApproxReduced(self, mu : paramList = []) -> sampList: """ Evaluate reduced representation of approximant at arbitrary parameter. Args: mu: Target parameter. """ mu = self.checkParameterList(mu) if (not hasattr(self, "lastSolvedApproxReduced") or self.lastSolvedApproxReduced != mu): vbMng(self, "INIT", "Evaluating approximant at mu = {}.".format(mu), 12) self.uApproxReduced = self.getNearestNeighbor(mu) vbMng(self, "DEL", "Done evaluating approximant.", 12) self.lastSolvedApproxReduced = mu return self.uApproxReduced