diff --git a/examples/all_forcing/all_forcing_engine.py b/examples/all_forcing/all_forcing_engine.py index e26119f..0074346 100644 --- a/examples/all_forcing/all_forcing_engine.py +++ b/examples/all_forcing/all_forcing_engine.py @@ -1,42 +1,43 @@ import numpy as np import fenics as fen from rrompy.hfengines.linear_problem import LaplaceBaseProblemEngine as LBPE from rrompy.utilities.base import verbosityManager as vbMng from rrompy.solver.fenics import fenics2Vector class AllForcingEngine(LBPE): def __init__(self, mu0:float, n : int = 30, degree_threshold : int = np.inf, verbosity : int = 10, timestamp : bool = True): super().__init__(mu0 = mu0, degree_threshold = degree_threshold, homogeneized = False, verbosity = verbosity, timestamp = timestamp) mesh = fen.RectangleMesh(fen.Point(-5., -5.), fen.Point(5., 5.), n, n) self.V = fen.FunctionSpace(mesh, "P", 1) self.nAs, self.nbs = 1, 4 x, y = fen.SpatialCoordinate(mesh)[:] scaling = (2. * np.pi) ** -1. r2 = x ** 2. + y ** 2. self.forcingCoeffs = [ scaling * fen.exp(- (r2 + 1. - 2. * x + 1. - 2. * y) / 2. / 4.) / 2., scaling * fen.exp(- (r2 + 1. + 2. * x + 1. + 2. * y) / 2. / 16.) / 4., - scaling * fen.exp(- (r2 + 1. + 2. * x + 1. - 2. * y) / 2. / 9.) / 30., scaling * fen.exp(- (r2 + 1. - 2. * x + 1. + 2. * y) / 2. / 25.) / 120.] def buildb(self): """Build terms of operator of linear system.""" if self.thbs[0] is None: - self.thbs = [None] * (self.nbs + self.homogeneized * self.nAs) + self.thbs = (self.getMonomialWeights(self.nbs) + + [None] * (self.homogeneized * self.nAs)) for j in range(self.nbs): if self.bs[j] is None: vbMng(self, "INIT", "Assembling forcing term b{}.".format(j), 20) parsRe = self.iterReduceQuadratureDegree([( self.forcingCoeffs[j], "forcingCoefficient")]) u0Re = self.DirichletDatum[0] L0Re = fen.dot(self.forcingCoeffs[j], self.v) * fen.dx DBCR = fen.DirichletBC(self.V, u0Re, self.DirichletBoundary) self.bs[j] = fenics2Vector(L0Re, parsRe, DBCR, 1) vbMng(self, "DEL", "Done assembling forcing term.", 20) self.setbHomogeneized() diff --git a/examples/pod/with_error_plot.py b/examples/pod/with_error_plot.py new file mode 100644 index 0000000..c436d68 --- /dev/null +++ b/examples/pod/with_error_plot.py @@ -0,0 +1,111 @@ +import numpy as np +from rrompy.hfengines.linear_problem import \ + HelmholtzSquareBubbleProblemEngine as HSBPE +from rrompy.reduction_methods.standard import RationalInterpolant as RI +from rrompy.reduction_methods.standard import RationalMovingLeastSquares as RIM +from rrompy.reduction_methods.standard import ReducedBasis as RB +from rrompy.parameter.parameter_sampling import QuadratureSampler as QS + +verb = 100 +sol = "single" +sol = "sweep" +algo = "RI" +algo = "RIM" +#algo = "RB" +polyBasis = "LEGENDRE" +polyBasis = "CHEBYSHEV" +#polyBasis = "MONOMIAL" +radialBasis = "GAUSSIAN" +#radialBasis = "THINPLATE" +#radialBasis = "MULTIQUADRIC" +#radialBasis = "NEARESTNEIGHBOR" +radialBasisDen = "GAUSSIAN" +radialBasisDen = "THINPLATE" +radialBasisDen = "MULTIQUADRIC" +radialBasisDen = "NEARESTNEIGHBOR" + +k0sC = np.power([7 + 0.j, 55 + 0.j], .5) +k0 = np.mean(k0sC ** 2.) ** .5 +ktar = 14 ** .5 + +n = 20 +solver = HSBPE(kappa = 12 ** .5, theta = np.pi / 3, n = 40, + verbosity = verb) + +params = {'N':1, 'M':1, 'S':30, 'POD':True, 'polybasis':polyBasis, + 'sampler':QS(k0sC, "CHEBYSHEV", 2.)} + +if algo == "RI": + approx = RI(solver, mu0 = k0, approxParameters = params, verbosity = verb) +elif algo == "RIM": + params["radialBasis"] = radialBasis + params["radialDirectionalWeights"] = [.75 * params["S"]] + params["radialBasisDen"] = radialBasisDen + params["nNearestNeighborDen"] = params["N"] + 1 + approx = RIM(solver, mu0 = k0, approxParameters = params, verbosity = verb) +else: + params.pop("N") + params.pop("M") + params.pop("polybasis") + approx = RB(solver, mu0 = k0, approxParameters = params, verbosity = verb) + +approx.setupApprox() +if sol == "single": + approx.plotSamples(what = "REAL") + approx.plotApprox(ktar, what = "REAL", name = "uApp") + approx.plotHF(ktar, what = "REAL", name = "uHF") + approx.plotErr(ktar, what = "REAL", name = "err") + approx.plotRes(ktar, what = "REAL", name = "res") + + appErr, solNorm = approx.normErr(ktar), approx.normHF(ktar) + resNorm, RHSNorm = approx.normRes(ktar), approx.normRHS(ktar) + print(('SolNorm:\t{}\nErr:\t{}\nErrRel:\t{}').format(solNorm, appErr, + np.divide(appErr, solNorm))) + print(('RHSNorm:\t{}\nRes:\t{}\nResRel:\t{}').format(RHSNorm, resNorm, + np.divide(resNorm, RHSNorm))) +poles = approx.getPoles() +try: + print('Poles:', poles ** 2.) +except: pass + +if sol == "sweep": + z0s = np.real(np.linspace(k0sC[0] ** 2., k0sC[1] ** 2., 100)) + k0s = z0s ** .5 + zl, zr = min(z0s), max(z0s) + approx.samplingEngine.verbosity = 0 + approx.trainedModel.verbosity = 0 + approx.verbosity = 0 + from matplotlib import pyplot as plt +# normRHS = approx.normRHS(k0s) + norm = approx.normHF(k0s) + normApp = approx.normApprox(k0s) +# res = approx.normRes(k0s) / normRHS +# err = approx.normErr(k0s) / norm + + plt.figure() + plt.semilogy(z0s, norm) + plt.semilogy(z0s, normApp, '--') + plt.semilogy(np.real(approx.mus.data) ** 2., + 1.05*np.max(norm)*np.ones_like(approx.mus.data, dtype = float), + 'rx') + plt.xlim([zl, zr]) + plt.grid() + plt.show() + plt.close() + +# plt.figure() +# plt.semilogy(z0s, res) +# plt.xlim([zl, zr]) +# plt.grid() +# plt.show() +# plt.close() +# +# plt.figure() +# plt.semilogy(z0s, err) +# plt.xlim([zl, zr]) +# plt.grid() +# plt.show() +# plt.close() + +#for j, k in enumerate(k0s): +# print(k ** 2., approx.getPoles(mu = k) ** 2., norm[j], normApp[j]) diff --git a/rrompy/parameter/parameter_sampling/quadrature_sampler.py b/rrompy/parameter/parameter_sampling/quadrature_sampler.py index 8bc4bb7..fe65111 100644 --- a/rrompy/parameter/parameter_sampling/quadrature_sampler.py +++ b/rrompy/parameter/parameter_sampling/quadrature_sampler.py @@ -1,86 +1,86 @@ # Copyright (C) 2018 by the RROMPy authors # # This file is part of RROMPy. # # RROMPy is free software: you can redistribute it and/or modify # it under the terms of the GNU Lesser General Public License as published by # the Free Software Foundation, either version 3 of the License, or # (at your option) any later version. # # RROMPy is distributed in the hope that it will be useful, # but WITHOUT ANY WARRANTY; without even the implied warranty of # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the # GNU Lesser General Public License for more details. # # You should have received a copy of the GNU Lesser General Public License # along with RROMPy. If not, see . # import numpy as np from .generic_sampler import GenericSampler from rrompy.utilities.base.types import List, paramList from rrompy.utilities.exception_manager import RROMPyException from rrompy.utilities.numerical import lowDiscrepancy, kroneckerer from rrompy.parameter import checkParameterList __all__ = ['QuadratureSampler'] class QuadratureSampler(GenericSampler): """Generator of quadrature sample points.""" _allowedKinds = ["UNIFORM", "CHEBYSHEV", "EXTENDEDCHEBYSHEV", "GAUSSLEGENDRE", "EXTENDEDGAUSSLEGENDRE"] def __init__(self, lims:paramList, kind : str = "UNIFORM", scalingExp : List[float] = None): super().__init__(lims = lims, scalingExp = scalingExp) self.kind = kind def __str__(self) -> str: return "{}_{}".format(super().__str__(), self.kind) def __repr__(self) -> str: return self.__str__() + " at " + hex(id(self)) @property def kind(self): """Value of kind.""" return self._kind @kind.setter def kind(self, kind): if kind.upper() not in self._allowedKinds: raise RROMPyException("Generator kind not recognized.") self._kind = kind.upper() def generatePoints(self, n:int, reorder : bool = True) -> paramList: """Array of sample points.""" n1d = int(np.ceil(n ** (1. / self.npar))) nleft, nright = 1, n1d ** self.npar xmat = np.empty((nright, self.npar), dtype = self.lims.dtype) for d in range(self.npar): nright //= n1d a = self.lims(0, d) ** self.scalingExp[d] b = self.lims(1, d) ** self.scalingExp[d] c, r = (a + b) / 2., (a - b) / 2. if self.kind == "UNIFORM": xd = np.linspace(a, b, n1d) elif self.kind in ["CHEBYSHEV", "EXTENDEDCHEBYSHEV"]: nodes = np.polynomial.chebyshev.chebgauss(n1d)[0] if n1d > 1 and self.kind == "EXTENDEDCHEBYSHEV": nodes /= nodes[0] xd = c + r * nodes elif self.kind in ["GAUSSLEGENDRE", "EXTENDEDGAUSSLEGENDRE"]: nodes = np.polynomial.legendre.leggauss(n1d)[0][::-1] if n1d > 1 and self.kind == "EXTENDEDCHEBYSHEV": nodes /= nodes[0] xd = c + r * nodes xd **= 1. / self.scalingExp[d] xmat[:, d] = kroneckerer(xd, nleft, nright) nleft *= n1d - x = checkParameterList(xmat, self.npar)[0] nright = n1d ** self.npar if nright > 1 and reorder: fejerOrdering = [nright - 1] + lowDiscrepancy(nright - 1) xmat = xmat[fejerOrdering, :] + x = checkParameterList(xmat, self.npar)[0] return x diff --git a/rrompy/reduction_methods/base/reduced_basis_utils.py b/rrompy/reduction_methods/base/reduced_basis_utils.py index c9740d1..37e8665 100644 --- a/rrompy/reduction_methods/base/reduced_basis_utils.py +++ b/rrompy/reduction_methods/base/reduced_basis_utils.py @@ -1,67 +1,67 @@ # Copyright (C) 2018 by the RROMPy authors # # This file is part of RROMPy. # # RROMPy is free software: you can redistribute it and/or modify # it under the terms of the GNU Lesser General Public License as published by # the Free Software Foundation, either version 3 of the License, or # (at your option) any later version. # # RROMPy is distributed in the hope that it will be useful, # but WITHOUT ANY WARRANTY; without even the implied warranty of # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the # GNU Lesser General Public License for more details. # # You should have received a copy of the GNU Lesser General Public License # along with RROMPy. If not, see . # import numpy as np from rrompy.utilities.base.types import (Np1D, Np2D, Tuple, List, ListAny, sampList) from rrompy.utilities.numerical import hashIdxToDerivative as hashI from rrompy.utilities.exception_manager import RROMPyAssert from rrompy.sampling import sampleList __all__ = ['projectAffineDecomposition'] def projectAffineDecomposition(As:List[Np2D], bs:List[Np1D], pMat:sampList, ARBsOld : List[Np2D] = None, bRBsOld : List[Np1D] = None, pMatOld : sampList = None)\ -> Tuple[List[Np2D], List[Np1D]]: """Project affine decomposition of linear system onto basis.""" RROMPyAssert((ARBsOld is None, bRBsOld is None), (pMatOld is None, pMatOld is None), "Old affine projected terms") if isinstance(pMat, (sampleList,)): pMat = pMat.data pMatH = pMat.T.conj() ARBs = [None] * len(As) bRBs = [None] * len(bs) if pMatOld is None: for j in range(max(len(As), len(bs))): if j < len(As): ARBs[j] = pMatH.dot(As[j].dot(pMat)) if j < len(bs): bRBs[j] = pMatH.dot(bs[j]) else: RROMPyAssert((len(ARBsOld), len(bRBsOld)), (len(As), len(bs)), "Old affine projected terms") if isinstance(pMatOld, (sampleList,)): pMatOld = pMatOld.data pMatOldH = pMatOld.T.conj() Sold = pMatOld.shape[1] Snew = pMat.shape[1] for j in range(max(len(As), len(bs))): if j < len(As): ARBs[j] = np.empty((Sold + Snew, Sold + Snew), - dtype = np.complex) + dtype = ARBsOld[j].dtype) ARBs[j][: Sold, : Sold] = ARBsOld[j] ARBs[j][: Sold, Sold :] = pMatOldH.dot(As[j].dot(pMat)) ARBs[j][Sold :, : Sold] = pMatH.dot(As[j].dot(pMatOld)) ARBs[j][Sold :, Sold :] = pMatH.dot(As[j].dot(pMat)) if j < len(bs): - bRBs[j] = np.empty((Sold + Snew), dtype = np.complex) + bRBs[j] = np.empty((Sold + Snew), dtype = bRBsOld[j].dtype) bRBs[j][: Sold] = bRBsOld[j] bRBs[j][Sold :] = pMatH.dot(bs[j]) return ARBs, bRBs diff --git a/rrompy/reduction_methods/greedy/generic_greedy_approximant.py b/rrompy/reduction_methods/greedy/generic_greedy_approximant.py index 93df37b..8ff5904 100644 --- a/rrompy/reduction_methods/greedy/generic_greedy_approximant.py +++ b/rrompy/reduction_methods/greedy/generic_greedy_approximant.py @@ -1,668 +1,671 @@ # Copyright (C) 2018 by the RROMPy authors # # This file is part of RROMPy. # # RROMPy is free software: you can redistribute it and/or modify # it under the terms of the GNU Lesser General Public License as published by # the Free Software Foundation, either version 3 of the License, or # (at your option) any later version. # # RROMPy is distributed in the hope that it will be useful, # but WITHOUT ANY WARRANTY; without even the implied warranty of # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the # GNU Lesser General Public License for more details. # # You should have received a copy of the GNU Lesser General Public License # along with RROMPy. If not, see . # from copy import deepcopy as copy import numpy as np from rrompy.reduction_methods.standard.generic_standard_approximant \ import GenericStandardApproximant from rrompy.utilities.base.types import (Np1D, Np2D, DictAny, HFEng, Tuple, List, normEng, paramVal, paramList, sampList) from rrompy.utilities.base import verbosityManager as vbMng from rrompy.utilities.expression import expressionEvaluator from rrompy.solver import normEngine from rrompy.utilities.exception_manager import (RROMPyException, RROMPyAssert, RROMPyWarning) from rrompy.parameter import checkParameterList, emptyParameterList __all__ = ['GenericGreedyApproximant'] def localL2Distance(mus:Np2D, badmus:Np2D) -> Np2D: return np.linalg.norm(np.tile(mus[..., np.newaxis], [1, 1, len(badmus)]) - badmus[..., np.newaxis].T, axis = 1) def pruneSamples(mus:paramList, badmus:paramList, tol : float = 1e-8) -> Np1D: """Remove from mus all the elements which are too close to badmus.""" if len(badmus) == 0: return mus proximity = np.min(localL2Distance(mus.data, badmus.data), axis = 1) return np.arange(len(mus))[proximity <= tol] class GenericGreedyApproximant(GenericStandardApproximant): """ ROM greedy interpolant computation for parametric problems (ABSTRACT). Args: HFEngine: HF problem solver. mu0(optional): Default parameter. Defaults to 0. approxParameters(optional): Dictionary containing values for main parameters of approximant. Recognized keys are: - 'POD': whether to compute POD of snapshots; defaults to True; - 'S': number of starting training points; - 'sampler': sample point generator; - 'greedyTol': uniform error tolerance for greedy algorithm; defaults to 1e-2; - 'collinearityTol': collinearity tolerance for greedy algorithm; defaults to 0.; - 'interactive': whether to interactively terminate greedy algorithm; defaults to False; - 'maxIter': maximum number of greedy steps; defaults to 1e2; - 'refinementRatio': ratio of test points to be exhausted before test set refinement; defaults to 0.2; - 'nTestPoints': number of test points; defaults to 5e2; - 'trainSetGenerator': training sample points generator; defaults to sampler. Defaults to empty dict. verbosity(optional): Verbosity level. Defaults to 10. Attributes: HFEngine: HF problem solver. mu0: Default parameter. mus: Array of snapshot parameters. approxParameters: Dictionary containing values for main parameters of approximant. Recognized keys are in parameterList. parameterListSoft: Recognized keys of soft approximant parameters: - 'POD': whether to compute POD of snapshots. - 'greedyTol': uniform error tolerance for greedy algorithm; - 'collinearityTol': collinearity tolerance for greedy algorithm; - 'interactive': whether to interactively terminate greedy algorithm; - 'maxIter': maximum number of greedy steps; - 'refinementRatio': ratio of test points to be exhausted before test set refinement; - 'nTestPoints': number of test points; - 'trainSetGenerator': training sample points generator. parameterListCritical: Recognized keys of critical approximant parameters: - 'S': total number of samples current approximant relies upon; - 'sampler': sample point generator. POD: whether to compute POD of snapshots. S: number of test points. sampler: Sample point generator. greedyTol: Uniform error tolerance for greedy algorithm. collinearityTol: Collinearity tolerance for greedy algorithm. interactive: whether to interactively terminate greedy algorithm. maxIter: maximum number of greedy steps. refinementRatio: ratio of training points to be exhausted before training set refinement. nTestPoints: number of starting training points. trainSetGenerator: training sample points generator. muBounds: list of bounds for parameter values. samplingEngine: Sampling engine. estimatorNormEngine: Engine for estimator norm computation. uHF: High fidelity solution(s) with parameter(s) lastSolvedHF as sampleList. lastSolvedHF: Parameter(s) corresponding to last computed high fidelity solution(s) as parameterList. uApproxReduced: Reduced approximate solution(s) with parameter(s) lastSolvedApprox as sampleList. lastSolvedApproxReduced: Parameter(s) corresponding to last computed reduced approximate solution(s) as parameterList. uApprox: Approximate solution(s) with parameter(s) lastSolvedApprox as sampleList. lastSolvedApprox: Parameter(s) corresponding to last computed approximate solution(s) as parameterList. """ TOL_INSTABILITY = 1e-6 def __init__(self, HFEngine:HFEng, mu0 : paramVal = None, approxParameters : DictAny = {}, verbosity : int = 10, timestamp : bool = True): self._preInit() self._addParametersToList(["greedyTol", "collinearityTol", "interactive", "maxIter", "refinementRatio", "nTestPoints"], [1e-2, 0., False, 1e2, .2, 5e2], ["trainSetGenerator"], ["AUTO"]) super().__init__(HFEngine = HFEngine, mu0 = mu0, approxParameters = approxParameters, verbosity = verbosity, timestamp = timestamp) self._postInit() @property def greedyTol(self): """Value of greedyTol.""" return self._greedyTol @greedyTol.setter def greedyTol(self, greedyTol): if greedyTol < 0: raise RROMPyException("greedyTol must be non-negative.") if hasattr(self, "_greedyTol") and self.greedyTol is not None: greedyTolold = self.greedyTol else: greedyTolold = -1 self._greedyTol = greedyTol self._approxParameters["greedyTol"] = self.greedyTol if greedyTolold != self.greedyTol: self.resetSamples() @property def collinearityTol(self): """Value of collinearityTol.""" return self._collinearityTol @collinearityTol.setter def collinearityTol(self, collinearityTol): if collinearityTol < 0: raise RROMPyException("collinearityTol must be non-negative.") if (hasattr(self, "_collinearityTol") and self.collinearityTol is not None): collinearityTolold = self.collinearityTol else: collinearityTolold = -1 self._collinearityTol = collinearityTol self._approxParameters["collinearityTol"] = self.collinearityTol if collinearityTolold != self.collinearityTol: self.resetSamples() @property def interactive(self): """Value of interactive.""" return self._interactive @interactive.setter def interactive(self, interactive): self._interactive = interactive @property def maxIter(self): """Value of maxIter.""" return self._maxIter @maxIter.setter def maxIter(self, maxIter): if maxIter <= 0: raise RROMPyException("maxIter must be positive.") if hasattr(self, "_maxIter") and self.maxIter is not None: maxIterold = self.maxIter else: maxIterold = -1 self._maxIter = maxIter self._approxParameters["maxIter"] = self.maxIter if maxIterold != self.maxIter: self.resetSamples() @property def refinementRatio(self): """Value of refinementRatio.""" return self._refinementRatio @refinementRatio.setter def refinementRatio(self, refinementRatio): if refinementRatio <= 0. or refinementRatio > 1.: raise RROMPyException(("refinementRatio must be between 0 " "(excluded) and 1.")) if (hasattr(self, "_refinementRatio") and self.refinementRatio is not None): refinementRatioold = self.refinementRatio else: refinementRatioold = -1 self._refinementRatio = refinementRatio self._approxParameters["refinementRatio"] = self.refinementRatio if refinementRatioold != self.refinementRatio: self.resetSamples() @property def nTestPoints(self): """Value of nTestPoints.""" return self._nTestPoints @nTestPoints.setter def nTestPoints(self, nTestPoints): if nTestPoints <= 0: raise RROMPyException("nTestPoints must be positive.") if not np.isclose(nTestPoints, np.int(nTestPoints)): raise RROMPyException("nTestPoints must be an integer.") nTestPoints = np.int(nTestPoints) if hasattr(self, "_nTestPoints") and self.nTestPoints is not None: nTestPointsold = self.nTestPoints else: nTestPointsold = -1 self._nTestPoints = nTestPoints self._approxParameters["nTestPoints"] = self.nTestPoints if nTestPointsold != self.nTestPoints: self.resetSamples() @property def trainSetGenerator(self): """Value of trainSetGenerator.""" return self._trainSetGenerator @trainSetGenerator.setter def trainSetGenerator(self, trainSetGenerator): if (isinstance(trainSetGenerator, (str,)) and trainSetGenerator.upper() == "AUTO"): trainSetGenerator = self.sampler if 'generatePoints' not in dir(trainSetGenerator): raise RROMPyException("trainSetGenerator type not recognized.") if (hasattr(self, '_trainSetGenerator') and self.trainSetGenerator not in [None, "AUTO"]): trainSetGeneratorOld = self.trainSetGenerator self._trainSetGenerator = trainSetGenerator self._approxParameters["trainSetGenerator"] = self.trainSetGenerator if (not 'trainSetGeneratorOld' in locals() or trainSetGeneratorOld != self.trainSetGenerator): self.resetSamples() def resetSamples(self): """Reset samples.""" super().resetSamples() self._mus = emptyParameterList() def initEstimatorNormEngine(self, normEngn : normEng = None): """Initialize estimator norm engine.""" if (normEngn is not None or not hasattr(self, "estimatorNormEngine") or self.estimatorNormEngine is None): if normEngn is None: if not hasattr(self.HFEngine, "energyNormPartialDualMatrix"): self.HFEngine.buildEnergyNormPartialDualForm() estimatorEnergyMatrix = ( self.HFEngine.energyNormPartialDualMatrix) else: if hasattr(normEngn, "buildEnergyNormPartialDualForm"): if not hasattr(normEngn, "energyNormPartialDualMatrix"): normEngn.buildEnergyNormPartialDualForm() estimatorEnergyMatrix = ( normEngn.energyNormPartialDualMatrix) else: estimatorEnergyMatrix = normEngn self.estimatorNormEngine = normEngine(estimatorEnergyMatrix) def _affineResidualMatricesContraction(self, rb:Np2D, rA : Np2D = None) \ -> Tuple[Np1D, Np1D, Np1D]: self.assembleReducedResidualBlocks(full = True) # 'ij,jk,ik->k', resbb, radiusb, radiusb.conj() ff = np.sum(self.trainedModel.data.resbb.dot(rb) * rb.conj(), axis = 0) if rA is None: return ff # 'ijk,jkl,il->l', resAb, radiusA, radiusb.conj() Lf = np.sum(np.tensordot(self.trainedModel.data.resAb, rA, 2) * rb.conj(), axis = 0) # 'ijkl,klt,ijt->t', resAA, radiusA, radiusA.conj() LL = np.sum(np.tensordot(self.trainedModel.data.resAA, rA, 2) * rA.conj(), axis = (0, 1)) return ff, Lf, LL def errorEstimator(self, mus:Np1D) -> Np1D: """Standard residual-based error estimator.""" self.setupApprox() mus = checkParameterList(mus, self.npar)[0] vbMng(self.trainedModel, "INIT", "Evaluating error estimator at mu = {}.".format(mus), 10) verb = self.trainedModel.verbosity self.trainedModel.verbosity = 0 uApproxRs = self.getApproxReduced(mus) muTestEff = mus ** self.HFEngine.rescalingExp radiusA = np.empty((len(self.HFEngine.thAs), len(mus)), - dtype = mus.dtype) + dtype = np.complex) radiusb = np.empty((len(self.HFEngine.thbs), len(mus)), - dtype = mus.dtype) + dtype = np.complex) for j, thA in enumerate(self.HFEngine.thAs): radiusA[j] = expressionEvaluator(thA[0], muTestEff) for j, thb in enumerate(self.HFEngine.thbs): radiusb[j] = expressionEvaluator(thb[0], muTestEff) radiusA = np.expand_dims(uApproxRs.data, 1) * radiusA ff, Lf, LL = self._affineResidualMatricesContraction(radiusb, radiusA) err = np.abs((LL - 2. * np.real(Lf) + ff) / ff) ** .5 self.trainedModel.verbosity = verb vbMng(self.trainedModel, "DEL", "Done evaluating error estimator", 10) return err def getMaxErrorEstimator(self, mus:paramList) -> Tuple[Np1D, int, float]: """ Compute maximum of (and index of maximum of) error estimator over given parameters. """ errorEstTest = self.errorEstimator(mus) idxMaxEst = [np.argmax(errorEstTest)] return errorEstTest, idxMaxEst, errorEstTest[idxMaxEst] def _isLastSampleCollinear(self) -> bool: """Check collinearity of last sample.""" if self.collinearityTol <= 0.: return False if self.POD: reff = self.samplingEngine.RPOD[:, -1] else: RROMPyWarning(("Repeated orthogonalization of the samples for " "collinearity check. Consider setting POD to " "True.")) if not hasattr(self, "_PODEngine"): from rrompy.sampling.base.pod_engine import PODEngine self._PODEngine = PODEngine(self.HFEngine) reff = self._PODEngine.generalizedQR(self.samplingEngine.samples, only_R = True)[:, -1] - return np.abs(reff[-1]) < self.collinearityTol * np.linalg.norm(reff) + cLevel = np.abs(reff[-1]) / np.linalg.norm(reff) + vbMng(self, "MAIN", "Collinearity indicator {:.4e}.".format(cLevel), 5) + return cLevel < self.collinearityTol def greedyNextSample(self, muidx:int, plotEst : bool = False)\ -> Tuple[Np1D, int, float, paramVal]: """Compute next greedy snapshot of solution map.""" RROMPyAssert(self._mode, message = "Cannot add greedy sample.") mus = copy(self.muTest[muidx]) self.muTest.pop(muidx) for mu in mus: vbMng(self, "MAIN", - ("Adding {}-th sample point at {} to training " + ("Adding sample point no. {} at {} to training " "set.").format(self.samplingEngine.nsamples + 1, mu), 2) self.mus.append(mu) self.samplingEngine.nextSample(mu) if self._isLastSampleCollinear(): RROMPyWarning("Collinearity above tolerance detected.") errorEstTest = np.empty(len(self.muTest)) errorEstTest[:] = np.nan return errorEstTest, [-1], np.nan, np.nan errorEstTest, muidx, maxErrorEst = self.getMaxErrorEstimator( self.muTest) - if plotEst and not np.any(np.isinf(errorEstTest)): + if (plotEst and not np.any(np.isnan(errorEstTest)) + and not np.any(np.isinf(errorEstTest))): musre = copy(self.muTest.re.data) from matplotlib import pyplot as plt plt.figure() errCP = copy(errorEstTest) while len(musre) > 0: if self.npar == 1: currIdx = np.arange(len(musre)) else: currIdx = np.where(np.isclose(np.sum( np.abs(musre[:, 1 :] - musre[0, 1 :]), 1), 0.))[0] plt.semilogy(musre[currIdx, 0], errCP[currIdx], 'k', linewidth = 1) musre = np.delete(musre, currIdx, 0) errCP = np.delete(errCP, currIdx) plt.semilogy([self.muTest.re(0, 0), self.muTest.re(-1, 0)], [self.greedyTol] * 2, 'r--') plt.semilogy(self.mus.re(0), 2. * self.greedyTol * np.ones(len(self.mus)), '*m') plt.semilogy(self.muTest.re(muidx, 0), maxErrorEst, 'xr') plt.grid() plt.show() plt.close() return errorEstTest, muidx, maxErrorEst, self.muTest[muidx] def _preliminaryTraining(self): """Initialize starting snapshots of solution map.""" RROMPyAssert(self._mode, message = "Cannot start greedy algorithm.") if self.samplingEngine.nsamples > 0: return self.computeScaleFactor() self.resetSamples() self.mus = self.trainSetGenerator.generatePoints(self.S)[ list(range(self.S))] muTestBase = self.sampler.generatePoints(self.nTestPoints) idxPop = pruneSamples(muTestBase ** self.HFEngine.rescalingExp, self.mus ** self.HFEngine.rescalingExp, 1e-10 * self.scaleFactor[0]) muTestBase.pop(idxPop) muTestBase = muTestBase.sort() muLast = copy(self.mus[-1]) self.mus.pop() if len(self.mus) > 0: vbMng(self, "MAIN", ("Adding first {} sample point{} at {} to training " "set.").format(self.S - 1, "" + "s" * (self.S > 2), self.mus), 2) self.samplingEngine.iterSample(self.mus) self.muTest = emptyParameterList() self.muTest.reset((len(muTestBase) + 1, self.mus.shape[1])) self.muTest[: -1] = muTestBase.data self.muTest[-1] = muLast.data def _enrichTestSet(self, nTest:int): """Add extra elements to test set.""" RROMPyAssert(self._mode, message = "Cannot enrich test set.") muTestExtra = self.sampler.generatePoints(2 * nTest) muTotal = copy(self.mus) muTotal.append(self.muTest) idxPop = pruneSamples(muTestExtra ** self.HFEngine.rescalingExp, muTotal ** self.HFEngine.rescalingExp, 1e-10 * self.scaleFactor[0]) muTestExtra.pop(idxPop) muTestNew = np.empty((len(self.muTest) + len(muTestExtra), self.muTest.shape[1]), dtype = np.complex) muTestNew[: len(self.muTest)] = self.muTest.data muTestNew[len(self.muTest) :] = muTestExtra.data self.muTest = checkParameterList(muTestNew, self.npar)[0].sort() vbMng(self, "MAIN", "Enriching test set by {} elements.".format(len(muTestExtra)), 5) def greedy(self, plotEst : bool = False): """Compute greedy snapshots of solution map.""" RROMPyAssert(self._mode, message = "Cannot start greedy algorithm.") if self.samplingEngine.nsamples > 0: return vbMng(self, "INIT", "Starting computation of snapshots.", 2) self._preliminaryTraining() nTest = self.nTestPoints muT0 = copy(self.muTest[-1]) errorEstTest, muidx, maxErrorEst, mu = self.greedyNextSample( [len(self.muTest) - 1], plotEst) if np.any(np.isnan(maxErrorEst)): RROMPyWarning(("Instability in a posteriori estimator. " "Starting preemptive greedy loop termination.")) self.muTest.append(muT0) self.mus.pop(-1) self.samplingEngine.popSample() self.setupApprox() else: vbMng(self, "MAIN", ("Uniform testing error estimate " "{:.4e}.").format(np.max(maxErrorEst)), 2) trainedModelOld = copy(self.trainedModel) while (self.samplingEngine.nsamples < self.maxIter and np.max(maxErrorEst) > self.greedyTol): if (1. - self.refinementRatio) * nTest > len(self.muTest): self._enrichTestSet(nTest) nTest = len(self.muTest) muTestOld, maxErrorEstOld = self.muTest, np.max(maxErrorEst) errorEstTest, muidx, maxErrorEst, mu = self.greedyNextSample( muidx, plotEst) vbMng(self, "MAIN", ("Uniform testing error estimate " "{:.4e}.").format(np.max(maxErrorEst)), 2) if (np.any(np.isnan(maxErrorEst)) or np.any(np.isinf(maxErrorEst)) or maxErrorEstOld < (np.max(maxErrorEst) * self.TOL_INSTABILITY)): RROMPyWarning(("Instability in a posteriori estimator. " "Starting preemptive greedy loop " "termination.")) self.muTest = muTestOld self.mus.pop(-1) self.samplingEngine.popSample() self.trainedModel.data = copy(trainedModelOld.data) break trainedModelOld.data = copy(self.trainedModel.data) if (self.interactive and np.max(maxErrorEst) <= self.greedyTol): vbMng(self, "MAIN", ("Required tolerance {} achieved. Want to decrease " "greedyTol and continue? " "Y/N").format(self.greedyTol), 0, end = "") increasemaxIter = input() if increasemaxIter.upper() == "Y": vbMng(self, "MAIN", "Reducing value of greedyTol...", 0) while np.max(maxErrorEst) <= self._greedyTol: self._greedyTol *= .5 if (self.interactive and self.samplingEngine.nsamples >= self.maxIter): vbMng(self, "MAIN", ("Maximum number of iterations {} reached. Want to " "increase maxIter and continue? " "Y/N").format(self.maxIter), 0, end = "") increasemaxIter = input() if increasemaxIter.upper() == "Y": vbMng(self, "MAIN", "Doubling value of maxIter...", 0) self._maxIter *= 2 vbMng(self, "DEL", ("Done computing snapshots (final snapshot count: " "{}).").format(self.samplingEngine.nsamples), 2) def checkComputedApprox(self) -> bool: """ Check if setup of new approximant is not needed. Returns: True if new setup is not needed. False otherwise. """ return (super().checkComputedApprox() and len(self.mus) == self.trainedModel.data.projMat.shape[1]) def assembleReducedResidualGramian(self, pMat:sampList): """ Build residual gramian of reduced linear system through projections. """ self.initEstimatorNormEngine() if (not hasattr(self.trainedModel.data, "gramian") or self.trainedModel.data.gramian is None): gramian = self.estimatorNormEngine.innerProduct(pMat, pMat) else: Sold = self.trainedModel.data.gramian.shape[0] S = len(self.mus) if Sold > S: gramian = self.trainedModel.data.gramian[: S, : S] else: idxOld = list(range(Sold)) idxNew = list(range(Sold, S)) gramian = np.empty((S, S), dtype = np.complex) gramian[: Sold, : Sold] = self.trainedModel.data.gramian gramian[: Sold, Sold :] = ( self.estimatorNormEngine.innerProduct(pMat(idxNew), pMat(idxOld))) gramian[Sold :, : Sold] = gramian[: Sold, Sold :].T.conj() gramian[Sold :, Sold :] = ( self.estimatorNormEngine.innerProduct(pMat(idxNew), pMat(idxNew))) self.trainedModel.data.gramian = gramian def assembleReducedResidualBlocksbb(self, bs:List[Np1D]): """ Build blocks (of type bb) of reduced linear system through projections. """ self.initEstimatorNormEngine() nbs = len(bs) if (not hasattr(self.trainedModel.data, "resbb") or self.trainedModel.data.resbb is None): resbb = np.empty((nbs, nbs), dtype = np.complex) for i in range(nbs): Mbi = bs[i] resbb[i, i] = self.estimatorNormEngine.innerProduct(Mbi, Mbi) for j in range(i): Mbj = bs[j] resbb[i, j] = self.estimatorNormEngine.innerProduct(Mbj, Mbi) for i in range(nbs): for j in range(i + 1, nbs): resbb[i, j] = resbb[j, i].conj() self.trainedModel.data.resbb = resbb def assembleReducedResidualBlocksAb(self, As:List[Np2D], bs:List[Np1D], pMat:sampList): """ Build blocks (of type Ab) of reduced linear system through projections. """ self.initEstimatorNormEngine() nAs = len(As) nbs = len(bs) S = len(self.mus) if (not hasattr(self.trainedModel.data, "resAb") or self.trainedModel.data.resAb is None): if not isinstance(pMat, (np.ndarray,)): pMat = pMat.data resAb = np.empty((nbs, S, nAs), dtype = np.complex) for j in range(nAs): MAj = As[j].dot(pMat) for i in range(nbs): Mbi = bs[i] resAb[i, :, j] = self.estimatorNormEngine.innerProduct(MAj, Mbi) else: Sold = self.trainedModel.data.resAb.shape[1] if Sold == S: return if Sold > S: resAb = self.trainedModel.data.resAb[:, : S, :] else: if not isinstance(pMat, (np.ndarray,)): pMat = pMat.data resAb = np.empty((nbs, S, nAs), dtype = np.complex) resAb[:, : Sold, :] = self.trainedModel.data.resAb for j in range(nAs): MAj = As[j].dot(pMat[:, Sold :]) for i in range(nbs): Mbi = bs[i] resAb[i, Sold :, j] = ( self.estimatorNormEngine.innerProduct(MAj, Mbi)) self.trainedModel.data.resAb = resAb def assembleReducedResidualBlocksAA(self, As:List[Np2D], pMat:sampList): """ Build blocks (of type AA) of reduced linear system through projections. """ self.initEstimatorNormEngine() nAs = len(As) S = len(self.mus) if (not hasattr(self.trainedModel.data, "resAA") or self.trainedModel.data.resAA is None): if not isinstance(pMat, (np.ndarray,)): pMat = pMat.data resAA = np.empty((S, nAs, S, nAs), dtype = np.complex) for i in range(nAs): MAi = As[i].dot(pMat) resAA[:, i, :, i] = ( self.estimatorNormEngine.innerProduct(MAi, MAi)) for j in range(i): MAj = As[j].dot(pMat) resAA[:, i, :, j] = ( self.estimatorNormEngine.innerProduct(MAj, MAi)) for i in range(nAs): for j in range(i + 1, nAs): resAA[:, i, :, j] = resAA[:, j, :, i].T.conj() else: Sold = self.trainedModel.data.resAA.shape[0] if Sold == S: return if Sold > S: resAA = self.trainedModel.data.resAA[: S, :, : S, :] else: if not isinstance(pMat, (np.ndarray,)): pMat = pMat.data resAA = np.empty((S, nAs, S, nAs), dtype = np.complex) resAA[: Sold, :, : Sold, :] = self.trainedModel.data.resAA for i in range(nAs): MAi = As[i].dot(pMat) resAA[: Sold, i, Sold :, i] = ( self.estimatorNormEngine.innerProduct(MAi[:, Sold :], MAi[:, : Sold])) resAA[Sold :, i, : Sold, i] = resAA[: Sold, i, Sold :, i].T.conj() resAA[Sold :, i, Sold :, i] = ( self.estimatorNormEngine.innerProduct(MAi[:, Sold :], MAi[:, Sold :])) for j in range(i): MAj = As[j].dot(pMat) resAA[: Sold, i, Sold :, j] = ( self.estimatorNormEngine.innerProduct(MAj[:, Sold :], MAi[:, : Sold])) resAA[Sold :, i, : Sold, j] = ( self.estimatorNormEngine.innerProduct(MAj[:, : Sold], MAi[:, Sold :])) resAA[Sold :, i, Sold :, j] = ( self.estimatorNormEngine.innerProduct(MAj[:, Sold :], MAi[:, Sold :])) for i in range(nAs): for j in range(i + 1, nAs): resAA[: Sold, i, Sold :, j] = ( resAA[Sold :, j, : Sold, i].T.conj()) resAA[Sold :, i, : Sold, j] = ( resAA[: Sold, j, Sold :, i].T.conj()) resAA[Sold :, i, Sold :, j] = ( resAA[Sold :, j, Sold :, i].T.conj()) self.trainedModel.data.resAA = resAA def assembleReducedResidualBlocks(self, full : bool = False): """Build affine blocks of affine decomposition of residual.""" self.assembleReducedResidualBlocksbb(self.HFEngine.bs) if full: pMat = self.trainedModel.data.projMat self.assembleReducedResidualBlocksAb(self.HFEngine.As, self.HFEngine.bs, pMat) self.assembleReducedResidualBlocksAA(self.HFEngine.As, pMat) diff --git a/rrompy/reduction_methods/greedy/rational_interpolant_greedy.py b/rrompy/reduction_methods/greedy/rational_interpolant_greedy.py index bcdd485..01514cb 100644 --- a/rrompy/reduction_methods/greedy/rational_interpolant_greedy.py +++ b/rrompy/reduction_methods/greedy/rational_interpolant_greedy.py @@ -1,466 +1,490 @@ # Copyright (C) 2018 by the RROMPy authors # # This file is part of RROMPy. # # RROMPy is free software: you can redistribute it and/or modify # it under the terms of the GNU Lesser General Public License as published by # the Free Software Foundation, either version 3 of the License, or # (at your option) any later version. # # RROMPy is distributed in the hope that it will be useful, # but WITHOUT ANY WARRANTY; without even the implied warranty of # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the # GNU Lesser General Public License for more details. # # You should have received a copy of the GNU Lesser General Public License # along with RROMPy. If not, see . # from copy import deepcopy as copy import numpy as np from .generic_greedy_approximant import (GenericGreedyApproximant, localL2Distance as lL2D) from rrompy.utilities.poly_fitting.polynomial import (polybases, polydomcoeff, PolynomialInterpolator as PI, polyvanderTotal as pvT) from rrompy.utilities.numerical import totalDegreeN from rrompy.utilities.expression import expressionEvaluator from rrompy.reduction_methods.standard import RationalInterpolant from rrompy.reduction_methods.trained_model import ( TrainedModelRational as tModel) from rrompy.reduction_methods.trained_model import TrainedModelData from rrompy.utilities.base.types import (Np1D, Tuple, DictAny, HFEng, paramVal, paramList) from rrompy.utilities.base import verbosityManager as vbMng from rrompy.utilities.poly_fitting import customFit from rrompy.utilities.exception_manager import (RROMPyWarning, RROMPyException, RROMPyAssert) from rrompy.parameter import checkParameterList __all__ = ['RationalInterpolantGreedy'] class RationalInterpolantGreedy(GenericGreedyApproximant, RationalInterpolant): """ ROM greedy rational interpolant computation for parametric problems. Args: HFEngine: HF problem solver. mu0(optional): Default parameter. Defaults to 0. approxParameters(optional): Dictionary containing values for main parameters of approximant. Recognized keys are: - 'POD': whether to compute POD of snapshots; defaults to True; - 'S': number of starting training points; - 'sampler': sample point generator; - - 'radialDirectionalWeights': radial basis weights for interpolant - numerator; defaults to 0, i.e. identity; - 'greedyTol': uniform error tolerance for greedy algorithm; defaults to 1e-2; - 'collinearityTol': collinearity tolerance for greedy algorithm; defaults to 0.; - 'interactive': whether to interactively terminate greedy algorithm; defaults to False; - 'maxIter': maximum number of greedy steps; defaults to 1e2; - - 'refinementRatio': ratio of training points to be exhausted - before training set refinement; defaults to 0.2; + - 'refinementRatio': ratio of test points to be exhausted before + test set refinement; defaults to 0.2; - 'nTestPoints': number of test points; defaults to 5e2; - 'trainSetGenerator': training sample points generator; defaults to sampler; - 'polybasis': type of basis for interpolation; defaults to 'MONOMIAL'; - 'errorEstimatorKind': kind of error estimator; available values include 'AFFINE', 'DISCREPANCY', 'INTERPOLATORY', 'EIM_INTERPOLATORY', and 'EIM_DIAGONAL'; defaults to 'AFFINE'; - 'interpRcond': tolerance for interpolation; defaults to None; - 'robustTol': 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. mus: Array of snapshot parameters. approxParameters: Dictionary containing values for main parameters of approximant. Recognized keys are in parameterList. parameterListSoft: Recognized keys of soft approximant parameters: - 'POD': whether to compute POD of snapshots. - - 'radialDirectionalWeights': radial basis weights for interpolant - numerator; defaults to 0, i.e. identity; - 'greedyTol': uniform error tolerance for greedy algorithm; - 'collinearityTol': collinearity tolerance for greedy algorithm; - 'interactive': whether to interactively terminate greedy algorithm; - 'maxIter': maximum number of greedy steps; - 'refinementRatio': ratio of test points to be exhausted before test set refinement; - 'nTestPoints': number of test points; - 'trainSetGenerator': training sample points generator; - 'errorEstimatorKind': kind of error estimator; - 'interpRcond': tolerance for interpolation; - 'robustTol': tolerance for robust rational denominator management. parameterListCritical: Recognized keys of critical approximant parameters: - 'S': total number of samples current approximant relies upon; - 'sampler': sample point generator. POD: whether to compute POD of snapshots. S: number of test points. sampler: Sample point generator. - radialDirectionalWeights: Radial basis weights for interpolant - numerator. greedyTol: uniform error tolerance for greedy algorithm. collinearityTol: Collinearity tolerance for greedy algorithm. interactive: whether to interactively terminate greedy algorithm. maxIter: maximum number of greedy steps. refinementRatio: ratio of training points to be exhausted before training set refinement. nTestPoints: number of starting training points. trainSetGenerator: training sample points generator. robustTol: tolerance for robust rational denominator management. errorEstimatorKind: kind of error estimator. interpRcond: tolerance for interpolation. robustTol: tolerance for robust rational denominator management. muBounds: list of bounds for parameter values. samplingEngine: Sampling engine. estimatorNormEngine: Engine for estimator norm computation. uHF: High fidelity solution(s) with parameter(s) lastSolvedHF as sampleList. lastSolvedHF: Parameter(s) corresponding to last computed high fidelity solution(s) as parameterList. uApproxReduced: Reduced approximate solution(s) with parameter(s) lastSolvedApprox as sampleList. lastSolvedApproxReduced: Parameter(s) corresponding to last computed reduced approximate solution(s) as parameterList. uApprox: Approximate solution(s) with parameter(s) lastSolvedApprox as sampleList. lastSolvedApprox: Parameter(s) corresponding to last computed approximate solution(s) as parameterList. """ _allowedEstimatorKinds = ["AFFINE", "DISCREPANCY", "INTERPOLATORY", "EIM_INTERPOLATORY", "EIM_DIAGONAL"] def __init__(self, HFEngine:HFEng, mu0 : paramVal = None, approxParameters : DictAny = {}, verbosity : int = 10, timestamp : bool = True): self._preInit() - self._addParametersToList(["polybasis", "errorEstimatorKind", - "interpRcond", "robustTol"], - ["MONOMIAL", "AFFINE", -1, 0]) + self._addParametersToList(["errorEstimatorKind"], ["AFFINE"], + toBeExcluded = ["M", "N", "polydegreetype", + "radialDirectionalWeights", + "nNearestNeighbor"]) super().__init__(HFEngine = HFEngine, mu0 = mu0, approxParameters = approxParameters, verbosity = verbosity, timestamp = timestamp) self._postInit() @property def E(self): """Value of E.""" self._E = self.sampleBatchIdx - 1 return self._E @E.setter def E(self, E): RROMPyWarning(("E is used just to simplify inheritance, and its value " "cannot be changed from that of sampleBatchIdx - 1.")) + @property + def polydegreetype(self): + """Value of polydegreetype.""" + return "TOTAL" + @polydegreetype.setter + def polydegreetype(self, polydegreetype): + RROMPyWarning(("polydegreetype is used just to simplify inheritance, " + "and its value cannot be changed from 'TOTAL'.")) + @property def polybasis(self): """Value of polybasis.""" return self._polybasis @polybasis.setter def polybasis(self, polybasis): try: polybasis = polybasis.upper().strip().replace(" ","") if polybasis not in polybases: raise RROMPyException("Sample type not recognized.") self._polybasis = polybasis except: RROMPyWarning(("Prescribed polybasis not recognized. Overriding " "to 'MONOMIAL'.")) self._polybasis = "MONOMIAL" self._approxParameters["polybasis"] = self.polybasis @property def errorEstimatorKind(self): """Value of errorEstimatorKind.""" return self._errorEstimatorKind @errorEstimatorKind.setter def errorEstimatorKind(self, errorEstimatorKind): errorEstimatorKind = errorEstimatorKind.upper() if errorEstimatorKind not in self._allowedEstimatorKinds: RROMPyWarning(("Error estimator kind not recognized. Overriding " "to 'AFFINE'.")) errorEstimatorKind = "AFFINE" self._errorEstimatorKind = errorEstimatorKind self._approxParameters["errorEstimatorKind"] = self.errorEstimatorKind def errorEstimator(self, mus:Np1D) -> Np1D: """Standard residual-based error estimator.""" if self.errorEstimatorKind == "AFFINE": return super().errorEstimator(mus) - self.setupApprox() + setupOK = self.setupApprox() + if not setupOK: + err = np.empty(len(mus)) + err[:] = np.nan + return err if self.errorEstimatorKind == "DIAGONAL": return self.errorEstimatorEIM(mus) mus = checkParameterList(mus, self.npar)[0] muCTest = self.trainedModel.centerNormalize(mus) vbMng(self.trainedModel, "INIT", "Evaluating error estimator at mu = {}.".format(mus), 10) verb = self.trainedModel.verbosity self.trainedModel.verbosity = 0 QTest = self.trainedModel.getQVal(mus) if self.errorEstimatorKind == "DISCREPANCY": nAs, nbs = len(self.HFEngine.thAs), len(self.HFEngine.thbs) muTrainEff = self.mus ** self.HFEngine.rescalingExp muTestEff = mus ** self.HFEngine.rescalingExp PTrain = self.trainedModel.getPVal(self.mus).data.T QTrain = self.trainedModel.getQVal(self.mus) PTest = self.trainedModel.getPVal(mus).data radiusAbTrain = np.empty((self.S, nAs * self.S + nbs), - dtype = self.mus.dtype) - radiusA = np.empty((self.S, nAs, len(mus)), dtype = mus.dtype) - radiusb = np.empty((nbs, len(mus)), dtype = mus.dtype) + dtype = np.complex) + radiusA = np.empty((self.S, nAs, len(mus)), dtype = np.complex) + radiusb = np.empty((nbs, len(mus)), dtype = np.complex) for j, thA in enumerate(self.HFEngine.thAs): idxs = j * self.S + np.arange(self.S) radiusAbTrain[:, idxs] = expressionEvaluator(thA[0], muTrainEff, (self.S, 1)) * PTrain radiusA[:, j] = PTest * expressionEvaluator(thA[0], muTestEff, (len(mus),)) for j, thb in enumerate(self.HFEngine.thbs): idx = nAs * self.S + j radiusAbTrain[:, idx] = QTrain * expressionEvaluator(thb[0], muTrainEff, (self.S,)) radiusb[j] = QTest * expressionEvaluator(thb[0], muTestEff, (len(mus),)) QRHSNorm2 = self._affineResidualMatricesContraction(radiusb) vanTrain, _, vanTrainIdxs = pvT(self._musUniqueCN, self.N, self.polybasis0, self._derIdxs, self._reorder) interpPQ = customFit(vanTrain[:, vanTrainIdxs], radiusAbTrain, rcond = self.interpRcond) vanTest, _, vanTestIdxs = pvT(muCTest, self.N, self.polybasis0) DradiusAb = vanTest[:, vanTestIdxs].dot(interpPQ) - radiusA -= DradiusAb[:, : - nbs].reshape(len(mus), -1, self.S).T - radiusb -= DradiusAb[:, - nbs :].T + radiusA = (radiusA + - DradiusAb[:, : - nbs].reshape(len(mus), -1, self.S).T) + radiusb = radiusb - DradiusAb[:, - nbs :].T ff, Lf, LL = self._affineResidualMatricesContraction(radiusb, radiusA) err = np.abs((LL - 2. * np.real(Lf) + ff) / QRHSNorm2) ** .5 else: #if self.errorEstimatorKind == "INTERPOLATORY": muCTrain = self.trainedModel.centerNormalize(self.mus) samplingRatio = np.prod(lL2D(muCTest.data, muCTrain.data), axis = 1) / np.abs(QTest) self.initEstimatorNormEngine() QTest = np.abs(QTest) sampRCP = copy(samplingRatio) idx_muTestSample = np.empty(self.sampleBatchSize, dtype = int) for j in range(self.sampleBatchSize): k = np.argmax(sampRCP) idx_muTestSample[j] = k if j + 1 < self.sampleBatchSize: musZero = self.trainedModel.centerNormalize(mus, mus[k]) sampRCP *= np.linalg.norm(musZero.data, axis = 1) mu_muTestSample = mus[idx_muTestSample] app_muTestSample = self.getApprox(mu_muTestSample) resmus = self.HFEngine.residual(mu_muTestSample, app_muTestSample, duality = False) RHSmus = self.HFEngine.residual(mu_muTestSample, None, duality = False) ressamples = (self.estimatorNormEngine.norm(resmus) / self.estimatorNormEngine.norm(RHSmus)) musT = copy(self.mus) musT.append(mu_muTestSample) musT = self.trainedModel.centerNormalize(musT) musC = self.trainedModel.centerNormalize(mus) resT = np.zeros(len(musT), dtype = np.complex) err = np.zeros(len(mus)) for l in range(len(mu_muTestSample)): resT[len(self.mus) + l] = (ressamples[l] * QTest[idx_muTestSample[l]]) p = PI() wellCond, msg = p.setupByInterpolation(musT, resT, self.M + 1, self.polybasis, self.verbosity >= 15, True, {}, {"rcond": self.interpRcond}) err += np.abs(p(musC)) resT[len(self.mus) + l] = 0. err /= QTest vbMng(self, "MAIN", msg, 15) self.trainedModel.verbosity = verb vbMng(self.trainedModel, "DEL", "Done evaluating error estimator", 10) return err def errorEstimatorEIM(self, mus:Np1D, return_max_idxs : bool = False) -> Np1D: """EIM-like interpolation error estimator.""" - self.setupApprox() + setupOK = self.setupApprox() + if not setupOK: + err = np.empty(len(mus)) + err[:] = np.nan + return err mus = checkParameterList(mus, self.npar)[0] vbMng(self.trainedModel, "INIT", "Evaluating error estimator at mu = {}.".format(mus), 10) verb = self.trainedModel.verbosity self.trainedModel.verbosity = 0 QTest = self.trainedModel.getQVal(mus) muCTest = self.trainedModel.centerNormalize(mus) muCTrain = self.trainedModel.centerNormalize(self.mus) vanderTest, _, vanderTestR = pvT(muCTest, self.E, self.polybasis) vanderTest = vanderTest[:, vanderTestR] vanderTestNext, _, vanderTestNextR = pvT(muCTest, self.E + 1, self.polybasis) vanderTestNext = vanderTestNext[:, vanderTestNextR[ vanderTest.shape[1] :]] idxsTest = np.arange(vanderTestNext.shape[1]) basis = np.zeros((len(idxsTest), 0), dtype = float) idxMaxEst = [] err = None while len(idxsTest) > 0: vanderTrial, _, vanderTrialR = pvT(muCTrain, self.E, self.polybasis) vanderTrial = vanderTrial[:, vanderTrialR] vanderTrialNext, _, vanderTrialNextR = pvT(muCTrain, self.E + 1, self.polybasis) vanderTrialNext = vanderTrialNext[:, vanderTrialNextR[ vanderTrial.shape[1] :]] vanderTrial = np.hstack((vanderTrial, vanderTrialNext.dot(basis).reshape( len(vanderTrialNext), basis.shape[1]))) valuesTrial = vanderTrialNext[:, idxsTest] vanderTestEff = np.hstack((vanderTest, vanderTestNext.dot(basis).reshape( len(vanderTestNext), basis.shape[1]))) vanderTestNextEff = vanderTestNext[:, idxsTest] coeffTest = np.linalg.solve(vanderTrial, valuesTrial) errTest = np.abs((vanderTestNextEff - vanderTestEff.dot(coeffTest)) / np.expand_dims(QTest, 1)) idxMaxErr = np.unravel_index(np.argmax(errTest), errTest.shape) idxMaxEst += [idxMaxErr[0]] if err is None: err = np.max(errTest, axis = 1) if not return_max_idxs: break muCTrain.append(muCTest[idxMaxErr[0]]) basis = np.pad(basis, [(0, 0), (0, 1)], "constant") basis[idxsTest[idxMaxErr[1]], -1] = 1. idxsTest = np.delete(idxsTest, idxMaxErr[1]) if self.errorEstimatorKind == "EIM_DIAGONAL": self.assembleReducedResidualBlocks(full = False) muTestEff = mus ** self.HFEngine.rescalingExp radiusb = np.empty((len(self.HFEngine.thbs), len(mus)), - dtype = mus.dtype) + dtype = np.complex) for j, thb in enumerate(self.HFEngine.thbs): radiusb[j] = expressionEvaluator(thb[0], muTestEff) bresb = self._affineResidualMatricesContraction(radiusb) self.assembleReducedResidualGramian(self.trainedModel.data.projMat) pDom = (polydomcoeff(self.E, self.polybasis) * self.trainedModel.data.P[(-1,) + (0,) * (self.npar - 1)]) LL = pDom.conj().dot(self.trainedModel.data.gramian.dot(pDom)) if not hasattr(self, "Anorm2Approx"): if self.HFEngine.nAs > 1: Ader = self.HFEngine.A(self.mu0, [1] + [0] * (self.npar - 1)) try: Adiag = self.scaleFactor[0] * Ader.diagonal() except: Adiag = self.scaleFactor[0] * np.diagonal(Ader) self.Anorm2Approx = np.mean(np.abs(Adiag) ** 2.) if (np.isclose(self.Anorm2Approx, 0.) or self.HFEngine.nAs <= 1): self.Anorm2Approx = 1 jOpt = np.abs(self.Anorm2Approx * LL / bresb) ** .5 err = jOpt * err else: #if self.errorEstimatorKind == "EIM_INTERPOLATORY": self.initEstimatorNormEngine() mu_muTestSample = mus[idxMaxEst[0]] app_muTestSample = self.getApprox(mu_muTestSample) resmu = self.HFEngine.residual(mu_muTestSample, app_muTestSample, duality = False) RHSmu = self.HFEngine.residual(mu_muTestSample, None, duality = False) jOpt = np.abs(self.estimatorNormEngine.norm(resmu)[0] / err[idxMaxEst[0]] / self.estimatorNormEngine.norm(RHSmu)[0]) err = jOpt * err self.trainedModel.verbosity = verb vbMng(self.trainedModel, "DEL", "Done evaluating error estimator", 10) if return_max_idxs: return err, idxMaxEst return err def getMaxErrorEstimator(self, mus:paramList) -> Tuple[Np1D, int, float]: """ Compute maximum of (and index of maximum of) error estimator over given parameters. """ if self.errorEstimatorKind[: 4] == "EIM_": errorEstTest, idxMaxEst = self.errorEstimatorEIM(mus, True) else: errorEstTest = self.errorEstimator(mus) idxMaxEst = np.empty(self.sampleBatchSize, dtype = int) errCP = copy(errorEstTest) for j in range(self.sampleBatchSize): k = np.argmax(errCP) idxMaxEst[j] = k if j + 1 < self.sampleBatchSize: musZero = self.trainedModel.centerNormalize(mus, mus[k]) errCP *= np.linalg.norm(musZero.data, axis = 1) maxEst = errorEstTest[idxMaxEst] return errorEstTest, idxMaxEst, maxEst def greedyNextSample(self, muidx:int, plotEst : bool = False)\ -> Tuple[Np1D, int, float, paramVal]: """Compute next greedy snapshot of solution map.""" RROMPyAssert(self._mode, message = "Cannot add greedy sample.") self.sampleBatchIdx += 1 self.sampleBatchSize = totalDegreeN(self.npar - 1, self.sampleBatchIdx) return super().greedyNextSample(muidx, plotEst) def _preliminaryTraining(self): """Initialize starting snapshots of solution map.""" RROMPyAssert(self._mode, message = "Cannot start greedy algorithm.") if self.samplingEngine.nsamples > 0: return S = self.S self.sampleBatchIdx, self.sampleBatchSize, self._S = -1, 0, 0 nextBatchSize = 1 while self._S + nextBatchSize <= S: self.sampleBatchIdx += 1 self.sampleBatchSize = nextBatchSize self._S += self.sampleBatchSize nextBatchSize = totalDegreeN(self.npar - 1, self.sampleBatchIdx + 1) super()._preliminaryTraining() def setupApprox(self, plotEst : bool = False): """ Compute rational interpolant. SVD-based robust eigenvalue management. """ if self.checkComputedApprox(): - return + return True RROMPyAssert(self._mode, message = "Cannot setup approximant.") vbMng(self, "INIT", "Setting up {}.".format(self.name()), 5) self.greedy(plotEst) self._S = len(self.mus) self._N, self._M = [self.E] * 2 if self.trainedModel is None: self.trainedModel = tModel() self.trainedModel.verbosity = self.verbosity self.trainedModel.timestamp = self.timestamp data = TrainedModelData(self.trainedModel.name(), self.mu0, self.samplingEngine.samples, self.scaleFactor, self.HFEngine.rescalingExp) data.mus = copy(self.mus) self.trainedModel.data = data else: self.trainedModel = self.trainedModel self.trainedModel.data.projMat = copy(self.samplingEngine.samples) self.trainedModel.data.mus = copy(self.mus) self.catchInstability = True if self.N > 0: - Q = self._setupDenominator()[0] + try: + Q = self._setupDenominator()[0] + except RROMPyException as RE: + RROMPyWarning(RE) + vbMng(self, "DEL", "Done setting up approximant.", 5) + return False else: Q = PI() Q.coeffs = np.ones(1, dtype = np.complex) Q.npar = 1 Q.polybasis = self.polybasis self.trainedModel.data.Q = copy(Q) - self.trainedModel.data.P = copy(self._setupNumerator()) + try: + self.trainedModel.data.P = copy(self._setupNumerator()) + except RROMPyException as RE: + RROMPyWarning(RE) + vbMng(self, "DEL", "Done setting up approximant.", 5) + return False self.trainedModel.data.approxParameters = copy(self.approxParameters) vbMng(self, "DEL", "Done setting up approximant.", 5) + return True diff --git a/rrompy/reduction_methods/pivoted/generic_pivoted_approximant.py b/rrompy/reduction_methods/pivoted/generic_pivoted_approximant.py index c579d10..13ff50d 100644 --- a/rrompy/reduction_methods/pivoted/generic_pivoted_approximant.py +++ b/rrompy/reduction_methods/pivoted/generic_pivoted_approximant.py @@ -1,529 +1,548 @@ # Copyright (C) 2018 by the RROMPy authors # # This file is part of RROMPy. # # RROMPy is free software: you can redistribute it and/or modify # it under the terms of the GNU Lesser General Public License as published by # the Free Software Foundation, either version 3 of the License, or # (at your option) any later version. # # RROMPy is distributed in the hope that it will be useful, # but WITHOUT ANY WARRANTY; without even the implied warranty of # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the # GNU Lesser General Public License for more details. # # You should have received a copy of the GNU Lesser General Public License # along with RROMPy. If not, see . # import numpy as np from copy import deepcopy as copy from rrompy.reduction_methods.base.generic_approximant import ( GenericApproximant) from rrompy.utilities.poly_fitting.polynomial import (polybases as ppb, PolynomialInterpolator as PI) from rrompy.utilities.poly_fitting.radial_basis import (polybases as rbpb, RadialBasisInterpolator as RBI) from rrompy.utilities.poly_fitting.moving_least_squares import ( polybases as mlspb, MovingLeastSquaresInterpolator as MLSI) from rrompy.sampling.pivoted import (SamplingEnginePivoted, SamplingEnginePivotedPOD) from rrompy.utilities.base.types import (ListAny, DictAny, HFEng, paramVal, paramList) from rrompy.utilities.base import verbosityManager as vbMng from rrompy.utilities.numerical import (fullDegreeN, totalDegreeN, nextDerivativeIndices) from rrompy.utilities.exception_manager import (RROMPyException, RROMPyAssert, RROMPyWarning) from rrompy.parameter import emptyParameterList __all__ = ['GenericPivotedApproximant'] class GenericPivotedApproximant(GenericApproximant): """ ROM pivoted approximant (with pole matching) computation for parametric problems (ABSTRACT). Args: HFEngine: HF problem solver. mu0(optional): Default parameter. Defaults to 0. directionPivot(optional): Pivot components. Defaults to [0]. approxParameters(optional): Dictionary containing values for main parameters of approximant. Recognized keys are: - 'POD': whether to compute POD of snapshots; defaults to True; - 'matchingWeight': weight for pole matching optimization; defaults to 1; - 'cutOffTolerance': tolerance for ignoring parasitic poles; defaults to np.inf; - 'cutOffType': rule for tolerance computation for parasitic poles; defaults to 'MAGNITUDE'; - 'S': total number of pivot samples current approximant relies upon; - 'samplerPivot': pivot sample point generator; - 'SMarginal': total number of marginal samples current approximant relies upon; - 'samplerMarginal': marginal sample point generator; - 'polybasisMarginal': type of polynomial basis for marginal interpolation; allowed values include 'MONOMIAL', 'CHEBYSHEV' and 'LEGENDRE'; defaults to 'MONOMIAL'; - 'MMarginal': degree of marginal interpolant; defaults to 0; - 'polydegreetypeMarginal': type of polynomial degree for marginal; defaults to 'TOTAL'; - 'radialDirectionalWeightsMarginal': radial basis weights for marginal interpolant; defaults to 0, i.e. identity; + - 'nNearestNeighborMarginal': number of marginal nearest neighbors + considered if polybasisMarginal allows; defaults to -1; - 'interpRcondMarginal': tolerance for marginal interpolation; defaults to None. 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': whether to compute POD of snapshots; - 'matchingWeight': weight for pole matching optimization; - 'cutOffTolerance': tolerance for ignoring parasitic poles; - 'cutOffType': rule for tolerance computation for parasitic poles; - 'polybasisMarginal': type of polynomial basis for marginal interpolation; - 'MMarginal': degree of marginal interpolant; - 'polydegreetypeMarginal': type of polynomial degree for marginal; - 'radialDirectionalWeightsMarginal': radial basis weights for marginal interpolant; + - 'nNearestNeighborMarginal': number of marginal nearest neighbors + considered if polybasisMarginal allows; - 'interpRcondMarginal': tolerance for marginal interpolation. parameterListCritical: Recognized keys of critical approximant parameters: - 'S': total number of pivot samples current approximant relies upon; - 'samplerPivot': pivot sample point generator; - 'SMarginal': total number of marginal samples current approximant relies upon; - 'samplerMarginal': marginal sample point generator. POD: Whether to compute POD of snapshots. matchingWeight: Weight for pole matching optimization. cutOffTolerance: Tolerance for ignoring parasitic poles. cutOffType: Rule for tolerance computation for parasitic poles. S: Total number of pivot samples current approximant relies upon. samplerPivot: Pivot sample point generator. SMarginal: Total number of marginal samples current approximant relies upon. samplerMarginal: Marginal sample point generator. polybasisMarginal: Type of polynomial basis for marginal interpolation. MMarginal: Degree of marginal interpolant. polydegreetypeMarginal: Type of polynomial degree for marginal. radialDirectionalWeightsMarginal: Radial basis weights for marginal interpolant. + nNearestNeighborMarginal: Number of marginal nearest neighbors + considered if polybasisMarginal allows. interpRcondMarginal: Tolerance for marginal interpolation. muBoundsPivot: list of bounds for pivot parameter values. muBoundsMarginal: list of bounds for marginal parameter values. samplingEngine: Sampling engine. uHF: High fidelity solution(s) with parameter(s) lastSolvedHF as sampleList. lastSolvedHF: Parameter(s) corresponding to last computed high fidelity solution(s) as parameterList. uApproxReduced: Reduced approximate solution(s) with parameter(s) lastSolvedApprox as sampleList. lastSolvedApproxReduced: Parameter(s) corresponding to last computed reduced approximate solution(s) as parameterList. uApprox: Approximate solution(s) with parameter(s) lastSolvedApprox as sampleList. lastSolvedApprox: Parameter(s) corresponding to last computed approximate solution(s) as parameterList. """ def __init__(self, HFEngine:HFEng, mu0 : paramVal = None, directionPivot : ListAny = [0], approxParameters : DictAny = {}, verbosity : int = 10, timestamp : bool = True): self._preInit() if len(directionPivot) > 1: raise RROMPyException(("Exactly 1 pivot parameter allowed in pole " "matching.")) from rrompy.parameter.parameter_sampling import QuadratureSampler as QS QSBase = QS([[0], [1]], "UNIFORM") self._addParametersToList(["matchingWeight", "cutOffTolerance", "cutOffType", "polybasisMarginal", "MMarginal", "polydegreetypeMarginal", "radialDirectionalWeightsMarginal", + "nNearestNeighborMarginal", "interpRcondMarginal"], [1, np.inf, "MAGNITUDE", "MONOMIAL", 0, - "TOTAL", 1, -1], + "TOTAL", 1, -1, -1], ["samplerPivot", "SMarginal", "samplerMarginal"], [QSBase, [1], QSBase]) del QS self._directionPivot = directionPivot super().__init__(HFEngine = HFEngine, mu0 = mu0, approxParameters = approxParameters, verbosity = verbosity, timestamp = timestamp) self._postInit() def setupSampling(self): """Setup sampling engine.""" RROMPyAssert(self._mode, message = "Cannot setup sampling engine.") if not hasattr(self, "_POD") or self._POD is None: return if self.POD: SamplingEngine = SamplingEnginePivotedPOD else: SamplingEngine = SamplingEnginePivoted self.samplingEngine = SamplingEngine(self.HFEngine, self.directionPivot, verbosity = self.verbosity, allowRepeatedSamples = True) @property def matchingWeight(self): """Value of matchingWeight.""" return self._matchingWeight @matchingWeight.setter def matchingWeight(self, matchingWeight): self._matchingWeight = matchingWeight self._approxParameters["matchingWeight"] = self.matchingWeight @property def cutOffTolerance(self): """Value of cutOffTolerance.""" return self._cutOffTolerance @cutOffTolerance.setter def cutOffTolerance(self, cutOffTolerance): self._cutOffTolerance = cutOffTolerance self._approxParameters["cutOffTolerance"] = self.cutOffTolerance @property def cutOffType(self): """Value of cutOffType.""" return self._cutOffType @cutOffType.setter def cutOffType(self, cutOffType): try: cutOffType = cutOffType.upper().strip().replace(" ","") if cutOffType not in ["MAGNITUDE", "POTENTIAL"]: raise RROMPyException("Prescribed cutOffType not recognized.") self._cutOffType = cutOffType except: RROMPyWarning(("Prescribed cutOffType not recognized. Overriding " "to 'MAGNITUDE'.")) self._cutOffType = "MAGNITUDE" self._approxParameters["cutOffType"] = self.cutOffType @property def SMarginal(self): """Value of SMarginal.""" return self._SMarginal @SMarginal.setter def SMarginal(self, SMarginal): if SMarginal <= 0: raise RROMPyException("SMarginal must be positive.") if hasattr(self, "_SMarginal") and self._SMarginal is not None: Sold = self.SMarginal else: Sold = -1 self._SMarginal = SMarginal self._approxParameters["SMarginal"] = self.SMarginal if Sold != self.SMarginal: self.resetSamples() @property def polybasisMarginal(self): """Value of polybasisMarginal.""" return self._polybasisMarginal @polybasisMarginal.setter def polybasisMarginal(self, polybasisMarginal): try: polybasisMarginal = polybasisMarginal.upper().strip().replace(" ", "") if polybasisMarginal not in ppb + rbpb + mlspb: raise RROMPyException( "Prescribed marginal polybasis not recognized.") self._polybasisMarginal = polybasisMarginal except: RROMPyWarning(("Prescribed marginal polybasis not recognized. " "Overriding to 'MONOMIAL'.")) self._polybasisMarginal = "MONOMIAL" self._approxParameters["polybasisMarginal"] = self.polybasisMarginal @property def MMarginal(self): """Value of MMarginal.""" return self._MMarginal @MMarginal.setter def MMarginal(self, MMarginal): if MMarginal < 0: raise RROMPyException("MMarginal must be non-negative.") self._MMarginal = MMarginal self._approxParameters["MMarginal"] = self.MMarginal @property def polydegreetypeMarginal(self): """Value of polydegreetypeMarginal.""" return self._polydegreetypeMarginal @polydegreetypeMarginal.setter def polydegreetypeMarginal(self, polydegreetypeM): try: polydegreetypeM = polydegreetypeM.upper().strip().replace(" ","") if polydegreetypeM not in ["TOTAL", "FULL"]: raise RROMPyException(("Prescribed polydegreetypeMarginal not " "recognized.")) self._polydegreetypeMarginal = polydegreetypeM except: RROMPyWarning(("Prescribed polydegreetypeMarginal not recognized. " "Overriding to 'TOTAL'.")) self._polydegreetypeMarginal = "TOTAL" self._approxParameters["polydegreetypeMarginal"] = ( self.polydegreetypeMarginal) @property def radialDirectionalWeightsMarginal(self): """Value of radialDirectionalWeightsMarginal.""" return self._radialDirectionalWeightsMarginal @radialDirectionalWeightsMarginal.setter def radialDirectionalWeightsMarginal(self, radialDirWeightsMarginal): self._radialDirectionalWeightsMarginal = radialDirWeightsMarginal self._approxParameters["radialDirectionalWeightsMarginal"] = ( self.radialDirectionalWeightsMarginal) + @property + def nNearestNeighborMarginal(self): + """Value of nNearestNeighborMarginal.""" + return self._nNearestNeighborMarginal + @nNearestNeighborMarginal.setter + def nNearestNeighborMarginal(self, nNearestNeighborMarginal): + self._nNearestNeighborMarginal = nNearestNeighborMarginal + self._approxParameters["nNearestNeighborMarginal"] = ( + self.nNearestNeighborMarginal) + @property def interpRcondMarginal(self): """Value of interpRcondMarginal.""" return self._interpRcondMarginal @interpRcondMarginal.setter def interpRcondMarginal(self, interpRcondMarginal): self._interpRcondMarginal = interpRcondMarginal self._approxParameters["interpRcondMarginal"] = ( self.interpRcondMarginal) @property def directionPivot(self): """Value of directionPivot. Its assignment may reset snapshots.""" return self._directionPivot @directionPivot.setter def directionPivot(self, directionPivot): if hasattr(self, '_directionPivot'): directionPivotOld = copy(self.directionPivot) else: directionPivotOld = None if (directionPivotOld is None or len(directionPivot) != len(directionPivotOld) or not directionPivot == directionPivotOld): self.resetSamples() self._directionPivot = directionPivot @property def directionMarginal(self): return [x for x in range(self.HFEngine.npar) \ if x not in self.directionPivot] @property def nparPivot(self): return len(self.directionPivot) @property def nparMarginal(self): return self.npar - self.nparPivot @property def rescalingExpPivot(self): return [self.HFEngine.rescalingExp[x] for x in self.directionPivot] @property def rescalingExpMarginal(self): return [self.HFEngine.rescalingExp[x] for x in self.directionMarginal] @property def muBoundsPivot(self): """Value of muBoundsPivot.""" return self.samplerPivot.lims @property def muBoundsMarginal(self): """Value of muBoundsMarginal.""" return self.samplerMarginal.lims @property def samplerPivot(self): """Value of samplerPivot.""" return self._samplerPivot @samplerPivot.setter def samplerPivot(self, samplerPivot): if 'generatePoints' not in dir(samplerPivot): raise RROMPyException("Pivot sampler type not recognized.") if hasattr(self, '_samplerPivot') and self._samplerPivot is not None: samplerOld = self.samplerPivot self._samplerPivot = samplerPivot self._approxParameters["samplerPivot"] = self.samplerPivot.__str__() if not 'samplerOld' in locals() or samplerOld != self.samplerPivot: self.resetSamples() @property def samplerMarginal(self): """Value of samplerMarginal.""" return self._samplerMarginal @samplerMarginal.setter def samplerMarginal(self, samplerMarginal): if 'generatePoints' not in dir(samplerMarginal): raise RROMPyException("Marginal sampler type not recognized.") if (hasattr(self, '_samplerMarginal') and self._samplerMarginal is not None): samplerOld = self.samplerMarginal self._samplerMarginal = samplerMarginal self._approxParameters["samplerMarginal"] = ( self.samplerMarginal.__str__()) if not 'samplerOld' in locals() or samplerOld != self.samplerMarginal: self.resetSamples() def resetSamples(self): """Reset samples.""" super().resetSamples() self._musMUniqueCN = None self._derMIdxs = None self._reorderM = None def setSamples(self, samplingEngine): """Copy samplingEngine and samples.""" super().setSamples(samplingEngine) self.mus = copy(self.samplingEngine[0].mus) for sEj in self.samplingEngine[1:]: self.mus.append(sEj.mus) def computeSnapshots(self): """Compute snapshots of solution map.""" RROMPyAssert(self._mode, message = "Cannot start snapshot computation.") if self.samplingEngine.nsamplesTot != self.S * self.SMarginal: self.computeScaleFactor() self.resetSamples() vbMng(self, "INIT", "Starting computation of snapshots.", 5) if self.HFEngine.As[0] is None: self.HFEngine.A(self.mu0) if self.HFEngine.bs[0] is None: self.HFEngine.b(self.mu0) self.musPivot = self.samplerPivot.generatePoints(self.S) self.musMarginal = self.samplerMarginal.generatePoints( self.SMarginal) self.mus = emptyParameterList() self.mus.reset((self.S * self.SMarginal, self.HFEngine.npar)) self.samplingEngine.resetHistory(self.SMarginal) for j, muMarg in enumerate(self.musMarginal): for k in range(j * self.S, (j + 1) * self.S): self.mus.data[k, self.directionPivot] = ( self.musPivot[k - j * self.S].data) self.mus.data[k, self.directionMarginal] = muMarg.data self.samplingEngine.iterSample(self.musPivot, self.musMarginal) if self.POD: self.samplingEngine.coalesceSamples(self.interpRcondMarginal) else: self.samplingEngine.coalesceSamples() vbMng(self, "DEL", "Done computing snapshots.", 5) def _setupMarginalInterpolationIndices(self): """Setup parameters for polyvander.""" RROMPyAssert(self._mode, message = "Cannot setup interpolation indices.") if (self._musMUniqueCN is None or len(self._reorderM) != len(self.musMarginal)): self._musMUniqueCN, musMIdxsTo, musMIdxs, musMCount = ( self.trainedModel.centerNormalizeMarginal(self.musMarginal)\ .unique(return_index = True, return_inverse = True, return_counts = True)) self._musMUnique = self.musMarginal[musMIdxsTo] self._derMIdxs = [None] * len(self._musMUniqueCN) self._reorderM = np.empty(len(musMIdxs), dtype = int) filled = 0 for j, cnt in enumerate(musMCount): self._derMIdxs[j] = nextDerivativeIndices([], self.nparMarginal, cnt) jIdx = np.nonzero(musMIdxs == j)[0] self._reorderM[jIdx] = np.arange(filled, filled + cnt) filled += cnt def _setupMarginalInterp(self): """Compute marginal interpolator.""" RROMPyAssert(self._mode, message = "Cannot setup numerator.") vbMng(self, "INIT", "Starting computation of marginal interpolator.", 7) self._setupMarginalInterpolationIndices() if self.polydegreetypeMarginal == "TOTAL": cfun = totalDegreeN else: cfun = fullDegreeN MM = copy(self.MMarginal) while len(self.musMarginal) < cfun(MM, self.nparMarginal): MM -= 1 if MM < self.MMarginal: RROMPyWarning( ("MMarginal too large compared to SMarginal. " "Reducing MMarginal by {}").format(self.MMarginal - MM)) self.MMarginal = MM mI = [] for j in range(len(self.musMarginal)): canonicalj = 1. * (np.arange(len(self.musMarginal)) == j) self._MMarginal = MM while self.MMarginal >= 0: if self.polybasisMarginal in ppb: p = PI() wellCond, msg = p.setupByInterpolation( self._musMUniqueCN, canonicalj, self.MMarginal, self.polybasisMarginal, self.verbosity >= 5, self.polydegreetypeMarginal == "TOTAL", {"derIdxs": self._derMIdxs, "reorder": self._reorderM, "scl": np.power(self.scaleFactorMarginal, -1.)}, {"rcond": self.interpRcondMarginal}) elif self.polybasisMarginal in rbpb: p = RBI() wellCond, msg = p.setupByInterpolation( self._musMUniqueCN, canonicalj, self.MMarginal, self.polybasisMarginal, self.radialDirectionalWeightsMarginal, self.verbosity >= 5, self.polydegreetypeMarginal == "TOTAL", {"derIdxs": self._derMIdxs, - "reorder": self._reorderM, - "scl": np.power(self.scaleFactorMarginal, -1.)}, + "reorder": self._reorderM, + "scl": np.power(self.scaleFactorMarginal, -1.), + "nNearestNeighbor" : self.nNearestNeighborMarginal}, {"rcond": self.interpRcondMarginal}) else:# if self.polybasisMarginal in mlspb: p = MLSI() wellCond, msg = p.setupByInterpolation( self._musMUniqueCN, canonicalj, self.MMarginal, self.polybasisMarginal, self.radialDirectionalWeightsMarginal, self.verbosity >= 5, self.polydegreetypeMarginal == "TOTAL", {"derIdxs": self._derMIdxs, - "reorder": self._reorderM, - "scl": np.power(self.scaleFactorMarginal, -1.)}) + "reorder": self._reorderM, + "scl": np.power(self.scaleFactorMarginal, -1.), + "nNearestNeighbor" : self.nNearestNeighborMarginal}) vbMng(self, "MAIN", msg, 5) if wellCond: break RROMPyWarning(("Polyfit is poorly conditioned. Reducing " "MMarginal by 1.")) self.MMarginal = self.MMarginal - 1 mI = mI + [copy(p)] vbMng(self, "DEL", "Done computing marginal interpolator.", 7) return mI def normApprox(self, mu:paramList) -> float: """ Compute norm of approximant at arbitrary parameter. Args: mu: Target parameter. Returns: Target norm of approximant. """ if not self.POD: return super().normApprox(mu) return np.linalg.norm(self.getApproxReduced(mu).data, axis = 0) def computeScaleFactor(self): """Compute parameter rescaling factor.""" RROMPyAssert(self._mode, message = "Cannot compute rescaling factor.") self.scaleFactorPivot = .5 * np.abs( self.muBoundsPivot[0] ** self.rescalingExpPivot - self.muBoundsPivot[1] ** self.rescalingExpPivot) self.scaleFactorMarginal = .5 * np.abs( self.muBoundsMarginal[0] ** self.rescalingExpMarginal - self.muBoundsMarginal[1] ** self.rescalingExpMarginal) self.scaleFactor = np.empty(self.npar) self.scaleFactor[self.directionPivot] = self.scaleFactorPivot self.scaleFactor[self.directionMarginal] = self.scaleFactorMarginal diff --git a/rrompy/reduction_methods/pivoted/rational_interpolant_pivoted.py b/rrompy/reduction_methods/pivoted/rational_interpolant_pivoted.py index 90a73fa..95da4a7 100644 --- a/rrompy/reduction_methods/pivoted/rational_interpolant_pivoted.py +++ b/rrompy/reduction_methods/pivoted/rational_interpolant_pivoted.py @@ -1,598 +1,624 @@ # Copyright (C) 2018 by the RROMPy authors # # This file is part of RROMPy. # # RROMPy is free software: you can redistribute it and/or modify # it under the terms of the GNU Lesser General Public License as published by # the Free Software Foundation, either version 3 of the License, or # (at your option) any later version. # # RROMPy is distributed in the hope that it will be useful, # but WITHOUT ANY WARRANTY; without even the implied warranty of # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the # GNU Lesser General Public License for more details. # # You should have received a copy of the GNU Lesser General Public License # along with RROMPy. If not, see . # from copy import deepcopy as copy import numpy as np from rrompy.reduction_methods.base import checkRobustTolerance from .generic_pivoted_approximant import GenericPivotedApproximant from rrompy.reduction_methods.standard.rational_interpolant import ( RationalInterpolant as RI) from rrompy.utilities.poly_fitting.polynomial import ( polybases as ppb, polyfitname, polyvander as pvP, polyvanderTotal as pvTP, PolynomialInterpolator as PI) from rrompy.utilities.poly_fitting.radial_basis import (polybases as rbpb, RadialBasisInterpolator as RBI) from rrompy.utilities.poly_fitting.moving_least_squares import ( polybases as mlspb, MovingLeastSquaresInterpolator as MLSI) from rrompy.reduction_methods.trained_model import (TrainedModelPivotedData, TrainedModelPivotedRational as tModel) from rrompy.utilities.base.types import (Np1D, Np2D, HFEng, DictAny, Tuple, List, ListAny, paramVal) from rrompy.utilities.base import verbosityManager as vbMng, freepar as fp from rrompy.utilities.numerical import (multifactorial, customPInv, fullDegreeN, totalDegreeN, degreeTotalToFull, fullDegreeMaxMask, totalDegreeMaxMask, nextDerivativeIndices, hashDerivativeToIdx as hashD, hashIdxToDerivative as hashI) from rrompy.utilities.exception_manager import (RROMPyException, RROMPyAssert, RROMPyWarning) from rrompy.parameter import checkParameter __all__ = ['RationalInterpolantPivoted'] class RationalInterpolantPivoted(GenericPivotedApproximant): """ ROM pivoted rational interpolant (with pole matching) computation for parametric problems. Args: HFEngine: HF problem solver. mu0(optional): Default parameter. Defaults to 0. directionPivot(optional): Pivot components. Defaults to [0]. approxParameters(optional): Dictionary containing values for main parameters of approximant. Recognized keys are: - 'POD': whether to compute POD of snapshots; defaults to True; - 'matchingWeight': weight for pole matching optimization; defaults to 1; - 'cutOffTolerance': tolerance for ignoring parasitic poles; defaults to np.inf; - 'cutOffType': rule for tolerance computation for parasitic poles; defaults to 'MAGNITUDE'; - 'S': total number of pivot samples current approximant relies upon; - 'samplerPivot': pivot sample point generator; - 'SMarginal': total number of marginal samples current approximant relies upon; - 'samplerMarginal': marginal sample point generator; - 'polybasisPivot': type of polynomial basis for pivot interpolation; defaults to 'MONOMIAL'; - 'polybasisMarginal': type of polynomial basis for marginal interpolation; allowed values include 'MONOMIAL', 'CHEBYSHEV' and 'LEGENDRE'; defaults to 'MONOMIAL'; - 'M': degree of rational interpolant numerator; defaults to 0; - 'N': degree of rational interpolant denominator; defaults to 0; - 'polydegreetype': type of polynomial degree; defaults to 'TOTAL'; - 'MMarginal': degree of marginal interpolant; defaults to 0; - 'polydegreetypeMarginal': type of polynomial degree for marginal; defaults to 'TOTAL'; - 'radialDirectionalWeightsPivot': radial basis weights for pivot numerator; defaults to 0, i.e. identity; - 'radialDirectionalWeightsMarginal': radial basis weights for marginal interpolant; defaults to 0, i.e. identity; + - 'nNearestNeighborPivot': number of pivot nearest neighbors + considered if polybasisPivot allows; defaults to -1; + - 'nNearestNeighborMarginal': number of marginal nearest neighbors + considered if polybasisMarginal allows; defaults to -1; - 'interpRcondPivot': tolerance for pivot interpolation; defaults to None; - 'interpRcondMarginal': tolerance for marginal interpolation; defaults to None; - 'robustTol': tolerance for robust rational denominator management; defaults to 0. Defaults to empty dict. verbosity(optional): Verbosity level. Defaults to 10. Attributes: HFEngine: HF problem solver. mu0: Default parameter. directionPivot: Pivot components. mus: Array of snapshot parameters. musPivot: Array of pivot snapshot parameters. musMarginal: Array of marginal snapshot parameters. approxParameters: Dictionary containing values for main parameters of approximant. Recognized keys are in parameterList. parameterListSoft: Recognized keys of soft approximant parameters: - 'POD': whether to compute POD of snapshots; - 'matchingWeight': weight for pole matching optimization; - 'cutOffTolerance': tolerance for ignoring parasitic poles; - 'cutOffType': rule for tolerance computation for parasitic poles; - 'polybasisPivot': type of polynomial basis for pivot interpolation; - 'polybasisMarginal': type of polynomial basis for marginal interpolation; - 'M': degree of rational interpolant numerator; - 'N': degree of rational interpolant denominator; - 'polydegreetype': type of polynomial degree; - 'MMarginal': degree of marginal interpolant; - 'polydegreetypeMarginal': type of polynomial degree for marginal; - 'radialDirectionalWeightsPivot': radial basis weights for pivot numerator; - 'radialDirectionalWeightsMarginal': radial basis weights for marginal interpolant; + - 'nNearestNeighborPivot': number of pivot nearest neighbors + considered if polybasisPivot allows; + - 'nNearestNeighborMarginal': number of marginal nearest neighbors + considered if polybasisMarginal allows; - 'interpRcondPivot': tolerance for pivot interpolation; - 'interpRcondMarginal': tolerance for marginal interpolation; - 'robustTol': tolerance for robust rational denominator management. parameterListCritical: Recognized keys of critical approximant parameters: - 'S': total number of pivot samples current approximant relies upon; - 'samplerPivot': pivot sample point generator; - 'SMarginal': total number of marginal samples current approximant relies upon; - 'samplerMarginal': marginal sample point generator. POD: Whether to compute POD of snapshots. matchingWeight: Weight for pole matching optimization. cutOffTolerance: Tolerance for ignoring parasitic poles. cutOffType: Rule for tolerance computation for parasitic poles. S: Total number of pivot samples current approximant relies upon. sampler: Pivot sample point generator. polybasisPivot: Type of polynomial basis for pivot interpolation. polybasisMarginal: Type of polynomial basis for marginal interpolation. M: Numerator degree of approximant. N: Denominator degree of approximant. polydegreetype: Type of polynomial degree. MMarginal: Degree of marginal interpolant. polydegreetypeMarginal: Type of polynomial degree for marginal. radialDirectionalWeightsPivot: Radial basis weights for pivot numerator. radialDirectionalWeightsMarginal: Radial basis weights for marginal interpolant. + nNearestNeighborPivot: Number of pivot nearest neighbors considered if + polybasisPivot allows. + nNearestNeighborMarginal: Number of marginal nearest neighbors + considered if polybasisMarginal allows. interpRcondPivot: Tolerance for pivot interpolation. interpRcondMarginal: Tolerance for marginal interpolation. robustTol: Tolerance for robust rational denominator management. muBoundsPivot: list of bounds for pivot parameter values. muBoundsMarginal: list of bounds for marginal parameter values. samplingEngine: Sampling engine. uHF: High fidelity solution(s) with parameter(s) lastSolvedHF as sampleList. lastSolvedHF: Parameter(s) corresponding to last computed high fidelity solution(s) as parameterList. uApproxReduced: Reduced approximate solution(s) with parameter(s) lastSolvedApprox as sampleList. lastSolvedApproxReduced: Parameter(s) corresponding to last computed reduced approximate solution(s) as parameterList. uApprox: Approximate solution(s) with parameter(s) lastSolvedApprox as sampleList. lastSolvedApprox: Parameter(s) corresponding to last computed approximate solution(s) as parameterList. Q: Numpy 1D vector containing complex coefficients of approximant denominator. P: Numpy 2D vector whose columns are FE dofs of coefficients of approximant numerator. """ def __init__(self, HFEngine:HFEng, mu0 : paramVal = None, directionPivot : ListAny = [0], approxParameters : DictAny = {}, verbosity : int = 10, timestamp : bool = True): self._preInit() self._addParametersToList(["polybasisPivot", "M", "N", "polydegreetype", "radialDirectionalWeightsPivot", + "nNearestNeighborPivot", "interpRcondPivot", "robustTol"], - ["MONOMIAL", 0, 0, "TOTAL", 1, -1, 0]) + ["MONOMIAL", 0, 0, "TOTAL", 1, -1, -1, 0]) super().__init__(HFEngine = HFEngine, mu0 = mu0, directionPivot = directionPivot, approxParameters = approxParameters, verbosity = verbosity, timestamp = timestamp) self._postInit() @property def polybasisPivot(self): """Value of polybasisPivot.""" return self._polybasisPivot @polybasisPivot.setter def polybasisPivot(self, polybasisPivot): try: polybasisPivot = polybasisPivot.upper().strip().replace(" ","") if polybasisPivot not in ppb + rbpb + mlspb: raise RROMPyException( "Prescribed pivot polybasis not recognized.") self._polybasisPivot = polybasisPivot except: RROMPyWarning(("Prescribed pivot polybasis not recognized. " "Overriding to 'MONOMIAL'.")) self._polybasisPivot = "MONOMIAL" self._approxParameters["polybasisPivot"] = self.polybasisPivot @property def polybasisPivot0(self): if "_" in self.polybasisPivot: return self.polybasisPivot.split("_")[0] return self.polybasisPivot @property def radialDirectionalWeightsPivot(self): """Value of radialDirectionalWeightsPivot.""" return self._radialDirectionalWeightsPivot @radialDirectionalWeightsPivot.setter def radialDirectionalWeightsPivot(self, radialDirectionalWeightsPivot): self._radialDirectionalWeightsPivot = radialDirectionalWeightsPivot self._approxParameters["radialDirectionalWeightsPivot"] = ( self.radialDirectionalWeightsPivot) + @property + def nNearestNeighborPivot(self): + """Value of nNearestNeighborPivot.""" + return self._nNearestNeighborPivot + @nNearestNeighborPivot.setter + def nNearestNeighborPivot(self, nNearestNeighborPivot): + self._nNearestNeighborPivot = nNearestNeighborPivot + self._approxParameters["nNearestNeighborPivot"] = ( + self.nNearestNeighborPivot) + @property def interpRcondPivot(self): """Value of interpRcondPivot.""" return self._interpRcondPivot @interpRcondPivot.setter def interpRcondPivot(self, interpRcondPivot): self._interpRcondPivot = interpRcondPivot self._approxParameters["interpRcondPivot"] = self.interpRcondPivot @property def M(self): """Value of M.""" return self._M @M.setter def M(self, M): if M < 0: raise RROMPyException("M must be non-negative.") self._M = M self._approxParameters["M"] = self.M @property def N(self): """Value of N.""" return self._N @N.setter def N(self, N): if N < 0: raise RROMPyException("N must be non-negative.") self._N = N self._approxParameters["N"] = self.N @property def polydegreetype(self): """Value of polydegreetype.""" return self._polydegreetype @polydegreetype.setter def polydegreetype(self, polydegreetype): try: polydegreetype = polydegreetype.upper().strip().replace(" ","") if polydegreetype not in ["TOTAL", "FULL"]: raise RROMPyException(("Prescribed polydegreetype not " "recognized.")) self._polydegreetype = polydegreetype except: RROMPyWarning(("Prescribed polydegreetype not recognized. " "Overriding to 'TOTAL'.")) self._polydegreetype = "TOTAL" self._approxParameters["polydegreetype"] = self.polydegreetype @property def robustTol(self): """Value of tolerance for robust rational denominator management.""" return self._robustTol @robustTol.setter def robustTol(self, robustTol): if robustTol < 0.: RROMPyWarning(("Overriding prescribed negative robustness " "tolerance to 0.")) robustTol = 0. self._robustTol = robustTol self._approxParameters["robustTol"] = self.robustTol def resetSamples(self): """Reset samples.""" super().resetSamples() self._musPUniqueCN = None self._derPIdxs = None self._reorderP = None def _setupPivotInterpolationIndices(self): """Setup parameters for polyvander.""" RROMPyAssert(self._mode, message = "Cannot setup interpolation indices.") if (self._musPUniqueCN is None or len(self._reorderP) != len(self.musPivot)): self._musPUniqueCN, musPIdxsTo, musPIdxs, musPCount = ( self.trainedModel.centerNormalizePivot(self.musPivot).unique( return_index = True, return_inverse = True, return_counts = True)) self._musPUnique = self.mus[musPIdxsTo] self._derPIdxs = [None] * len(self._musPUniqueCN) self._reorderP = np.empty(len(musPIdxs), dtype = int) filled = 0 for j, cnt in enumerate(musPCount): self._derPIdxs[j] = nextDerivativeIndices([], self.nparPivot, cnt) jIdx = np.nonzero(musPIdxs == j)[0] self._reorderP[jIdx] = np.arange(filled, filled + cnt) filled += cnt def _setupDenominator(self): """Compute rational denominator.""" RROMPyAssert(self._mode, message = "Cannot setup denominator.") vbMng(self, "INIT", "Starting computation of denominator.", 7) NinvD = None N0 = copy(self.N) qs = [] self.verbosity -= 10 for j in range(len(self.musMarginal)): self._N = N0 while self.N > 0: if NinvD != self.N: invD, fitinvP = self._computeInterpolantInverseBlocks() NinvD = self.N if self.POD: ev, eV = RI.findeveVGQR(self, self.samplingEngine.RPOD[j], invD) else: ev, eV = RI.findeveVGExplicit(self, self.samplingEngine.samples[j], invD) nevBad = checkRobustTolerance(ev, self.robustTol) if nevBad <= 1: break if self.catchInstability: raise RROMPyException(("Instability in denominator " "computation: eigenproblem is " "poorly conditioned.")) RROMPyWarning(("Smallest {} eigenvalues below tolerance. " "Reducing N by 1.").format(nevBad)) self.N = self.N - 1 if self.N <= 0: self._N = 0 eV = np.ones((1, 1)) q = PI() q.npar = self.nparPivot q.polybasis = self.polybasisPivot0 if self.polydegreetype == "TOTAL": q.coeffs = degreeTotalToFull(tuple([self.N + 1] * q.npar), q.npar, eV[:, 0]) else: q.coeffs = eV[:, 0].reshape([self.N + 1] * q.npar) qs = qs + [copy(q)] self.verbosity += 10 vbMng(self, "DEL", "Done computing denominator.", 7) return qs, fitinvP def _setupNumerator(self): """Compute rational numerator.""" RROMPyAssert(self._mode, message = "Cannot setup numerator.") vbMng(self, "INIT", "Starting computation of numerator.", 7) Qevaldiag = np.zeros((len(self.musPivot), len(self.musPivot)), dtype = np.complex) verb = self.trainedModel.verbosity self.trainedModel.verbosity = 0 self._setupPivotInterpolationIndices() cfun = totalDegreeN if self.polydegreetype == "TOTAL" else fullDegreeN M = copy(self.M) while len(self.musPivot) < cfun(M, self.nparPivot): M -= 1 if M < self.M: RROMPyWarning(("M too large compared to S. Reducing M by " "{}").format(self.M - M)) self.M = M tensor_idx = 0 ps = [] for k, muM in enumerate(self.musMarginal): self._M = M idxGlob = 0 for j, derIdxs in enumerate(self._derPIdxs): mujEff = [fp] * self.npar for jj, kk in enumerate(self.directionPivot): mujEff[kk] = self._musPUnique[j, jj] for jj, kk in enumerate(self.directionMarginal): mujEff[kk] = muM(0, jj) mujEff = checkParameter(mujEff, self.npar) nder = len(derIdxs) idxLoc = np.arange(len(self.musPivot))[ (self._reorderP >= idxGlob) * (self._reorderP < idxGlob + nder)] idxGlob += nder Qval = [0] * nder for der in range(nder): derIdx = hashI(der, self.nparPivot) derIdxEff = [0] * self.npar sclEff = [0] * self.npar for jj, kk in enumerate(self.directionPivot): derIdxEff[kk] = derIdx[jj] sclEff[kk] = self.scaleFactorPivot[jj] ** -1. Qval[der] = (self.trainedModel.getQVal(mujEff, derIdxEff, scl = sclEff) / multifactorial(derIdx)) for derU, derUIdx in enumerate(derIdxs): for derQ, derQIdx in enumerate(derIdxs): diffIdx = [x - y for (x, y) in zip(derUIdx, derQIdx)] if all([x >= 0 for x in diffIdx]): diffj = hashD(diffIdx) Qevaldiag[idxLoc[derU], idxLoc[derQ]] = Qval[diffj] while self.M >= 0: if self.polybasisPivot in ppb: p = PI() wellCond, msg = p.setupByInterpolation( self._musPUniqueCN, Qevaldiag, self.M, self.polybasisPivot, self.verbosity >= 5, self.polydegreetype == "TOTAL", {"derIdxs": self._derPIdxs, "reorder": self._reorderP, "scl": np.power(self.scaleFactorPivot, -1.)}, {"rcond": self.interpRcondPivot}) elif self.polybasisPivot in rbpb: p = RBI() wellCond, msg = p.setupByInterpolation( self._musPUniqueCN, Qevaldiag, self.M, self.polybasisPivot, self.radialDirectionalWeightsPivot, self.verbosity >= 5, self.polydegreetype == "TOTAL", {"derIdxs": self._derPIdxs, - "reorder": self._reorderP, - "scl": np.power(self.scaleFactorPivot, -1.)}, + "reorder": self._reorderP, + "scl": np.power(self.scaleFactorPivot, -1.), + "nNearestNeighbor" : self.nNearestNeighborPivot}, {"rcond": self.interpRcondPivot}) else:# if self.polybasisPivot in mlspb: p = MLSI() wellCond, msg = p.setupByInterpolation( self._musPUniqueCN, Qevaldiag, self.M, self.polybasisPivot, self.radialDirectionalWeightsPivot, self.verbosity >= 5, self.polydegreetype == "TOTAL", {"derIdxs": self._derPIdxs, - "reorder": self._reorderP, - "scl": np.power(self.scaleFactorPivot, -1.)}) + "reorder": self._reorderP, + "scl": np.power(self.scaleFactorPivot, -1.), + "nNearestNeighbor" : self.nNearestNeighborPivot}) vbMng(self, "MAIN", msg, 5) if wellCond: break if self.catchInstability: raise RROMPyException(("Instability in numerator " "computation: polyfit is " "poorly conditioned.")) RROMPyWarning(("Polyfit is poorly conditioned. " "Reducing M by 1.")) self.M = self.M - 1 tensor_idx_new = tensor_idx + Qevaldiag.shape[1] if self.POD: p.postmultiplyTensorize(self.samplingEngine.RPODCoalesced.T[ tensor_idx : tensor_idx_new, :]) else: p.pad(tensor_idx, len(self.mus) - tensor_idx_new) tensor_idx = tensor_idx_new ps = ps + [copy(p)] self.trainedModel.verbosity = verb vbMng(self, "DEL", "Done computing numerator.", 7) return ps def setupApprox(self): """ Compute rational interpolant. SVD-based robust eigenvalue management. """ if self.checkComputedApprox(): return RROMPyAssert(self._mode, message = "Cannot setup approximant.") vbMng(self, "INIT", "Setting up {}.". format(self.name()), 5) self.computeSnapshots() if self.trainedModel is None: self.trainedModel = tModel() self.trainedModel.verbosity = self.verbosity self.trainedModel.timestamp = self.timestamp data = TrainedModelPivotedData(self.trainedModel.name(), self.mu0, self.samplingEngine.samplesCoalesced, self.scaleFactor, self.HFEngine.rescalingExp, self.directionPivot) data.musPivot = copy(self.musPivot) data.musMarginal = copy(self.musMarginal) self.trainedModel.data = data else: self.trainedModel = self.trainedModel self.trainedModel.data.projMat = copy( self.samplingEngine.samplesCoalesced) self.trainedModel.data.marginalInterp = self._setupMarginalInterp() if self.N > 0: Qs = self._setupDenominator()[0] else: Q = PI() Q.npar = self.nparPivot - Q.coeffs = np.ones(tuple([1] * Q.npar), dtype = np.complex) + Q.coeffs = np.ones(tuple([1] * Q.npar), + dtype = self.musMarginal.dtype) Q.polybasis = self.polybasisPivot0 Qs = [Q for _ in range(len(self.musMarginal))] self.trainedModel.data._temporary = True self.trainedModel.data.Qs = Qs self.trainedModel.data.Ps = self._setupNumerator() vbMng(self, "INIT", "Matching poles.", 10) self.trainedModel.initializeFromRational(self.HFEngine, self.matchingWeight, self.POD) vbMng(self, "DEL", "Done matching poles.", 10) if not np.isinf(self.cutOffTolerance): vbMng(self, "INIT", "Recompressing by cut-off.", 10) msg = self.trainedModel.recompressByCutOff([-1., 1.], self.cutOffTolerance, self.cutOffType) vbMng(self, "DEL", "Done recompressing." + msg, 10) self.trainedModel.data.approxParameters = copy(self.approxParameters) vbMng(self, "DEL", "Done setting up approximant.", 5) def _computeInterpolantInverseBlocks(self) -> Tuple[List[Np2D], Np2D]: """ Compute inverse factors for minimal interpolant target functional. """ RROMPyAssert(self._mode, message = "Cannot solve eigenvalue problem.") self._setupPivotInterpolationIndices() cfun = totalDegreeN if self.polydegreetype == "TOTAL" else fullDegreeN N = copy(self.N) while len(self.musPivot) < cfun(N, self.nparPivot): N -= 1 if N < self.N: RROMPyWarning(("N too large compared to S. Reducing N by " "{}").format(self.N - N)) self.N = N while self.N >= 0: if self.polydegreetype == "TOTAL": TE, _, argIdxs = pvTP(self._musPUniqueCN, self.N, self.polybasisPivot0, self._derPIdxs, self._reorderP, scl = np.power(self.scaleFactorPivot, -1.)) TE = TE[:, argIdxs] idxsB = totalDegreeMaxMask(self.N, self.nparPivot) else: #if self.polydegreetype == "FULL": TE = pvP(self._musPUniqueCN, [self.N] * self.nparPivot, self.polybasisPivot0, self._derPIdxs, self._reorderP, scl = np.power(self.scaleFactorPivot, -1.)) idxsB = fullDegreeMaxMask(self.N, self.nparPivot) fitOut = customPInv(TE, rcond = self.interpRcondPivot, full = True) vbMng(self, "MAIN", ("Fitting {} samples with degree {} through {}... " "Conditioning of pseudoinverse system: {:.4e}.").format( TE.shape[0], self.N, polyfitname(self.polybasisPivot0), fitOut[1][1][0] / fitOut[1][1][-1]), 5) if fitOut[1][0] == TE.shape[1]: fitinvP = fitOut[0][idxsB, :] break RROMPyWarning("Polyfit is poorly conditioned. Reducing N by 1.") self.N -= 1 if self.N < 0: raise RROMPyException(("Instability in computation of " "denominator. Aborting.")) TN, _, argIdxs = pvTP(self._musPUniqueCN, self.N, self.polybasisPivot0, self._derPIdxs, self._reorderP, scl = np.power(self.scaleFactorPivot, -1.)) TN = TN[:, argIdxs] invD = [None] * (len(idxsB)) for k in range(len(idxsB)): pseudoInv = np.diag(fitinvP[k, :]) idxGlob = 0 for j, derIdxs in enumerate(self._derPIdxs): nder = len(derIdxs) idxGlob += nder if nder > 1: idxLoc = np.arange(len(self.musPivot))[ (self._reorderP >= idxGlob - nder) * (self._reorderP < idxGlob)] invLoc = fitinvP[k, idxLoc] pseudoInv[np.ix_(idxLoc, idxLoc)] = 0. for diffj, diffjIdx in enumerate(derIdxs): for derQ, derQIdx in enumerate(derIdxs): derUIdx = [x - y for (x, y) in zip(diffjIdx, derQIdx)] if all([x >= 0 for x in derUIdx]): derU = hashD(derUIdx) pseudoInv[idxLoc[derU], idxLoc[derQ]] = ( invLoc[diffj]) invD[k] = pseudoInv.dot(TN) return invD, fitinvP def getResidues(self, *args, **kwargs) -> Np1D: """ Obtain approximant residues. Returns: Matrix with residues as columns. """ return self.trainedModel.getResidues(*args, **kwargs) diff --git a/rrompy/reduction_methods/standard/__init__.py b/rrompy/reduction_methods/standard/__init__.py index 0ddf98c..a9a5e6a 100644 --- a/rrompy/reduction_methods/standard/__init__.py +++ b/rrompy/reduction_methods/standard/__init__.py @@ -1,29 +1,31 @@ # Copyright (C) 2018 by the RROMPy authors # # This file is part of RROMPy. # # RROMPy is free software: you can redistribute it and/or modify # it under the terms of the GNU Lesser General Public License as published by # the Free Software Foundation, either version 3 of the License, or # (at your option) any later version. # # RROMPy is distributed in the hope that it will be useful, # but WITHOUT ANY WARRANTY; without even the implied warranty of # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the # GNU Lesser General Public License for more details. # # You should have received a copy of the GNU Lesser General Public License # along with RROMPy. If not, see . # from .generic_standard_approximant import GenericStandardApproximant from .rational_interpolant import RationalInterpolant +from .rational_moving_least_squares import RationalMovingLeastSquares from .reduced_basis import ReducedBasis __all__ = [ 'GenericStandardApproximant', 'RationalInterpolant', + 'RationalMovingLeastSquares', 'ReducedBasis' ] diff --git a/rrompy/reduction_methods/standard/rational_interpolant.py b/rrompy/reduction_methods/standard/rational_interpolant.py index 55cb25c..b0ae423 100644 --- a/rrompy/reduction_methods/standard/rational_interpolant.py +++ b/rrompy/reduction_methods/standard/rational_interpolant.py @@ -1,539 +1,561 @@ # Copyright (C) 2018 by the RROMPy authors # # This file is part of RROMPy. # # RROMPy is free software: you can redistribute it and/or modify # it under the terms of the GNU Lesser General Public License as published by # the Free Software Foundation, either version 3 of the License, or # (at your option) any later version. # # RROMPy is distributed in the hope that it will be useful, # but WITHOUT ANY WARRANTY; without even the implied warranty of # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the # GNU Lesser General Public License for more details. # # You should have received a copy of the GNU Lesser General Public License # along with RROMPy. If not, see . # from copy import deepcopy as copy import numpy as np from rrompy.reduction_methods.base import checkRobustTolerance from .generic_standard_approximant import GenericStandardApproximant from rrompy.utilities.poly_fitting.polynomial import ( polybases as ppb, polyfitname, polyvander as pvP, polyvanderTotal as pvTP, PolynomialInterpolator as PI) from rrompy.utilities.poly_fitting.radial_basis import (polybases as rbpb, RadialBasisInterpolator as RBI) from rrompy.utilities.poly_fitting.moving_least_squares import ( polybases as mlspb, MovingLeastSquaresInterpolator as MLSI) from rrompy.reduction_methods.trained_model import ( TrainedModelRational as tModel) from rrompy.reduction_methods.trained_model import TrainedModelData from rrompy.utilities.base.types import (Np1D, Np2D, HFEng, DictAny, Tuple, List, paramVal, sampList) from rrompy.utilities.base import verbosityManager as vbMng from rrompy.utilities.numerical import (multifactorial, customPInv, fullDegreeN, totalDegreeN, degreeTotalToFull, fullDegreeMaxMask, totalDegreeMaxMask, nextDerivativeIndices, hashDerivativeToIdx as hashD, hashIdxToDerivative as hashI) from rrompy.utilities.exception_manager import (RROMPyException, RROMPyAssert, RROMPyWarning) __all__ = ['RationalInterpolant'] class RationalInterpolant(GenericStandardApproximant): """ ROM rational interpolant computation for parametric problems. Args: HFEngine: HF problem solver. mu0(optional): Default parameter. Defaults to 0. approxParameters(optional): Dictionary containing values for main parameters of approximant. Recognized keys are: - 'POD': whether to compute POD of snapshots; defaults to True; - 'S': total number of samples current approximant relies upon; - 'sampler': sample point generator; - 'polybasis': type of polynomial basis for interpolation; defaults to 'MONOMIAL'; - 'M': degree of rational interpolant numerator; defaults to 0; - 'N': degree of rational interpolant denominator; defaults to 0; - 'polydegreetype': type of polynomial degree; defaults to 'TOTAL'; - 'radialDirectionalWeights': radial basis weights for interpolant numerator; defaults to 0, i.e. identity; + - 'nNearestNeighbor': mumber of nearest neighbors considered in + numerator if polybasis allows; defaults to -1; - 'interpRcond': tolerance for interpolation; defaults to None; - 'robustTol': 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. mus: Array of snapshot parameters. approxParameters: Dictionary containing values for main parameters of approximant. Recognized keys are in parameterList. parameterListSoft: Recognized keys of soft approximant parameters: - 'POD': whether to compute POD of snapshots; - 'polybasis': type of polynomial basis for interpolation; - 'M': degree of rational interpolant numerator; - 'N': degree of rational interpolant denominator; - 'polydegreetype': type of polynomial degree; - 'radialDirectionalWeights': radial basis weights for interpolant numerator; + - 'nNearestNeighbor': mumber of nearest neighbors considered in + numerator if polybasis allows; - 'interpRcond': tolerance for interpolation via numpy.polyfit; - 'robustTol': tolerance for robust rational denominator management. parameterListCritical: Recognized keys of critical approximant parameters: - 'S': total number of samples current approximant relies upon; - 'sampler': sample point generator. POD: Whether to compute POD of snapshots. S: Number of solution snapshots over which current approximant is based upon. sampler: Sample point generator. polybasis: type of polynomial basis for interpolation. M: Numerator degree of approximant. N: Denominator degree of approximant. polydegreetype: Type of polynomial degree. radialDirectionalWeights: Radial basis weights for interpolant numerator. + nNearestNeighbor: Number of nearest neighbors considered in numerator + if polybasis allows. interpRcond: Tolerance for interpolation via numpy.polyfit. robustTol: Tolerance for robust rational denominator management. muBounds: list of bounds for parameter values. samplingEngine: Sampling engine. uHF: High fidelity solution(s) with parameter(s) lastSolvedHF as sampleList. lastSolvedHF: Parameter(s) corresponding to last computed high fidelity solution(s) as parameterList. uApproxReduced: Reduced approximate solution(s) with parameter(s) lastSolvedApprox as sampleList. lastSolvedApproxReduced: Parameter(s) corresponding to last computed reduced approximate solution(s) as parameterList. uApprox: Approximate solution(s) with parameter(s) lastSolvedApprox as sampleList. lastSolvedApprox: Parameter(s) corresponding to last computed approximate solution(s) as parameterList. Q: Numpy 1D vector containing complex coefficients of approximant denominator. P: Numpy 2D vector whose columns are FE dofs of coefficients of approximant numerator. """ def __init__(self, HFEngine:HFEng, mu0 : paramVal = None, approxParameters : DictAny = {}, verbosity : int = 10, timestamp : bool = True): self._preInit() self._addParametersToList(["polybasis", "M", "N", "polydegreetype", - "radialDirectionalWeights", "interpRcond", + "radialDirectionalWeights", + "nNearestNeighbor", "interpRcond", "robustTol"], - ["MONOMIAL", 0, 0, "TOTAL", 1, -1, 0]) + ["MONOMIAL", 0, 0, "TOTAL", 1, -1, -1, 0]) super().__init__(HFEngine = HFEngine, mu0 = mu0, approxParameters = approxParameters, verbosity = verbosity, timestamp = timestamp) self.catchInstability = False self._postInit() @property def polybasis(self): """Value of polybasis.""" return self._polybasis @polybasis.setter def polybasis(self, polybasis): try: polybasis = polybasis.upper().strip().replace(" ","") if polybasis not in ppb + rbpb + mlspb: raise RROMPyException("Prescribed polybasis not recognized.") self._polybasis = polybasis except: RROMPyWarning(("Prescribed polybasis not recognized. Overriding " "to 'MONOMIAL'.")) self._polybasis = "MONOMIAL" self._approxParameters["polybasis"] = self.polybasis @property def polybasis0(self): if "_" in self.polybasis: return self.polybasis.split("_")[0] return self.polybasis @property def interpRcond(self): """Value of interpRcond.""" return self._interpRcond @interpRcond.setter def interpRcond(self, interpRcond): self._interpRcond = interpRcond self._approxParameters["interpRcond"] = self.interpRcond @property def radialDirectionalWeights(self): """Value of radialDirectionalWeights.""" return self._radialDirectionalWeights @radialDirectionalWeights.setter def radialDirectionalWeights(self, radialDirectionalWeights): self._radialDirectionalWeights = radialDirectionalWeights self._approxParameters["radialDirectionalWeights"] = ( self.radialDirectionalWeights) + @property + def nNearestNeighbor(self): + """Value of nNearestNeighbor.""" + return self._nNearestNeighbor + @nNearestNeighbor.setter + def nNearestNeighbor(self, nNearestNeighbor): + self._nNearestNeighbor = nNearestNeighbor + self._approxParameters["nNearestNeighbor"] = self.nNearestNeighbor + @property def M(self): """Value of M.""" return self._M @M.setter def M(self, M): if M < 0: raise RROMPyException("M must be non-negative.") self._M = M self._approxParameters["M"] = self.M @property def N(self): """Value of N.""" return self._N @N.setter def N(self, N): if N < 0: raise RROMPyException("N must be non-negative.") self._N = N self._approxParameters["N"] = self.N @property def polydegreetype(self): """Value of polydegreetype.""" return self._polydegreetype @polydegreetype.setter def polydegreetype(self, polydegreetype): try: polydegreetype = polydegreetype.upper().strip().replace(" ","") if polydegreetype not in ["TOTAL", "FULL"]: raise RROMPyException(("Prescribed polydegreetype not " "recognized.")) self._polydegreetype = polydegreetype except: RROMPyWarning(("Prescribed polydegreetype not recognized. " "Overriding to 'TOTAL'.")) self._polydegreetype = "TOTAL" self._approxParameters["polydegreetype"] = self.polydegreetype @property def robustTol(self): """Value of tolerance for robust rational denominator management.""" return self._robustTol @robustTol.setter def robustTol(self, robustTol): if robustTol < 0.: RROMPyWarning(("Overriding prescribed negative robustness " "tolerance to 0.")) robustTol = 0. self._robustTol = robustTol self._approxParameters["robustTol"] = self.robustTol def resetSamples(self): """Reset samples.""" super().resetSamples() self._musUniqueCN = None self._derIdxs = None self._reorder = None 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.mus): self._musUniqueCN, musIdxsTo, musIdxs, musCount = ( self.trainedModel.centerNormalize(self.mus).unique( return_index = True, return_inverse = True, return_counts = True)) self._musUnique = self.mus[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.mus.shape[1], cnt) jIdx = np.nonzero(musIdxs == j)[0] self._reorder[jIdx] = np.arange(filled, filled + cnt) filled += cnt def _setupDenominator(self): """Compute rational denominator.""" RROMPyAssert(self._mode, message = "Cannot setup denominator.") vbMng(self, "INIT", "Starting computation of denominator.", 7) while self.N > 0: invD, fitinv = self._computeInterpolantInverseBlocks() if self.POD: ev, eV = self.findeveVGQR(self.samplingEngine.RPOD, invD) else: ev, eV = self.findeveVGExplicit(self.samplingEngine.samples, invD) nevBad = checkRobustTolerance(ev, self.robustTol) if nevBad <= 1: break if self.catchInstability: raise RROMPyException(("Instability in denominator " "computation: eigenproblem is poorly " "conditioned.")) RROMPyWarning(("Smallest {} eigenvalues below tolerance. Reducing " "N by 1.").format(nevBad)) self.N = self.N - 1 if self.N <= 0: self._N = 0 eV = np.ones((1, 1)) q = PI() q.npar = self.npar q.polybasis = self.polybasis0 if self.polydegreetype == "TOTAL": q.coeffs = degreeTotalToFull(tuple([self.N + 1] * self.npar), self.npar, eV[:, 0]) else: q.coeffs = eV[:, 0].reshape([self.N + 1] * self.npar) vbMng(self, "DEL", "Done computing denominator.", 7) return q, fitinv def _setupNumerator(self): """Compute rational numerator.""" RROMPyAssert(self._mode, message = "Cannot setup numerator.") vbMng(self, "INIT", "Starting computation of numerator.", 7) Qevaldiag = np.zeros((len(self.mus), len(self.mus)), dtype = np.complex) verb = self.trainedModel.verbosity self.trainedModel.verbosity = 0 self._setupInterpolationIndices() idxGlob = 0 for j, derIdxs in enumerate(self._derIdxs): nder = len(derIdxs) idxLoc = np.arange(len(self.mus))[(self._reorder >= idxGlob) * (self._reorder < idxGlob + nder)] idxGlob += nder Qval = [0] * nder for der in range(nder): derIdx = hashI(der, self.npar) Qval[der] = (self.trainedModel.getQVal( self._musUnique[j], derIdx, scl = np.power(self.scaleFactor, -1.)) / multifactorial(derIdx)) for derU, derUIdx in enumerate(derIdxs): for derQ, derQIdx in enumerate(derIdxs): diffIdx = [x - y for (x, y) in zip(derUIdx, derQIdx)] if all([x >= 0 for x in diffIdx]): diffj = hashD(diffIdx) Qevaldiag[idxLoc[derU], idxLoc[derQ]] = Qval[diffj] if self.POD: Qevaldiag = Qevaldiag.dot(self.samplingEngine.RPOD.T) self.trainedModel.verbosity = verb cfun = totalDegreeN if self.polydegreetype == "TOTAL" else fullDegreeN M = copy(self.M) while len(self.mus) < cfun(M, self.npar): M -= 1 if M < self.M: RROMPyWarning(("M too large compared to S. Reducing M by " "{}").format(self.M - M)) self.M = M while self.M >= 0: if self.polybasis in ppb: p = PI() wellCond, msg = p.setupByInterpolation( self._musUniqueCN, Qevaldiag, self.M, self.polybasis, self.verbosity >= 5, self.polydegreetype == "TOTAL", {"derIdxs": self._derIdxs, "reorder": self._reorder, "scl": np.power(self.scaleFactor, -1.)}, {"rcond": self.interpRcond}) elif self.polybasis in rbpb: p = RBI() wellCond, msg = p.setupByInterpolation( self._musUniqueCN, Qevaldiag, self.M, self.polybasis, self.radialDirectionalWeights, self.verbosity >= 5, self.polydegreetype == "TOTAL", {"derIdxs": self._derIdxs, "reorder": self._reorder, - "scl": np.power(self.scaleFactor, -1.)}, + "scl": np.power(self.scaleFactor, -1.), + "nNearestNeighbor": self.nNearestNeighbor}, {"rcond": self.interpRcond}) else:# if self.polybasis in mlspb: p = MLSI() wellCond, msg = p.setupByInterpolation( self._musUniqueCN, Qevaldiag, self.M, self.polybasis, self.radialDirectionalWeights, self.verbosity >= 5, self.polydegreetype == "TOTAL", {"derIdxs": self._derIdxs, "reorder": self._reorder, - "scl": np.power(self.scaleFactor, -1.)}) + "scl": np.power(self.scaleFactor, -1.), + "nNearestNeighbor": self.nNearestNeighbor}) vbMng(self, "MAIN", msg, 5) if wellCond: break if self.catchInstability: raise RROMPyException(("Instability in numerator computation: " "polyfit is poorly conditioned.")) RROMPyWarning("Polyfit is poorly conditioned. Reducing M by 1.") self.M = self.M - 1 if self.M < 0: raise RROMPyException(("Instability in computation of numerator. " "Aborting.")) vbMng(self, "DEL", "Done computing numerator.", 7) return p def setupApprox(self): """ Compute rational interpolant. SVD-based robust eigenvalue management. """ if self.checkComputedApprox(): return RROMPyAssert(self._mode, message = "Cannot setup approximant.") vbMng(self, "INIT", "Setting up {}.". format(self.name()), 5) self.computeSnapshots() if self.trainedModel is None: self.trainedModel = tModel() self.trainedModel.verbosity = self.verbosity self.trainedModel.timestamp = self.timestamp data = TrainedModelData(self.trainedModel.name(), self.mu0, self.samplingEngine.samples, self.scaleFactor, self.HFEngine.rescalingExp) data.mus = copy(self.mus) self.trainedModel.data = data else: self.trainedModel = self.trainedModel self.trainedModel.data.projMat = copy(self.samplingEngine.samples) if self.N > 0: Q = self._setupDenominator()[0] else: Q = PI() Q.coeffs = np.ones(tuple([1] * self.npar), dtype = np.complex) Q.npar = self.npar Q.polybasis = self.polybasis self.trainedModel.data.Q = Q self.trainedModel.data.P = self._setupNumerator() self.trainedModel.data.approxParameters = copy(self.approxParameters) vbMng(self, "DEL", "Done setting up approximant.", 5) def _computeInterpolantInverseBlocks(self) -> Tuple[List[Np2D], Np2D]: """ Compute inverse factors for minimal interpolant target functional. """ RROMPyAssert(self._mode, message = "Cannot solve eigenvalue problem.") self._setupInterpolationIndices() cfun = totalDegreeN if self.polydegreetype == "TOTAL" else fullDegreeN N = copy(self.N) while len(self.mus) < cfun(N, self.npar): N -= 1 if N < self.N: RROMPyWarning(("N too large compared to S. Reducing N by " "{}").format(self.N - N)) self.N = N while self.N >= 0: if self.polydegreetype == "TOTAL": TE, _, argIdxs = pvTP(self._musUniqueCN, self.N, self.polybasis0, self._derIdxs, self._reorder, scl = np.power(self.scaleFactor, -1.)) TE = TE[:, argIdxs] idxsB = totalDegreeMaxMask(self.N, self.npar) else: #if self.polydegreetype == "FULL": TE = pvP(self._musUniqueCN, [self.N] * self.npar, self.polybasis0, self._derIdxs, self._reorder, scl = np.power(self.scaleFactor, -1.)) idxsB = fullDegreeMaxMask(self.N, self.npar) fitOut = customPInv(TE, rcond = self.interpRcond, full = True) vbMng(self, "MAIN", ("Fitting {} samples with degree {} through {}... " "Conditioning of pseudoinverse system: {:.4e}.").format( TE.shape[0], self.N, polyfitname(self.polybasis0), fitOut[1][1][0] / fitOut[1][1][-1]), 5) if fitOut[1][0] == TE.shape[1]: fitinv = fitOut[0][idxsB, :] break + if self.catchInstability: + raise RROMPyException(("Instability in denominator " + "computation: polyfit is poorly " + "conditioned.")) RROMPyWarning("Polyfit is poorly conditioned. Reducing N by 1.") self.N = self.N - 1 if self.polydegreetype == "TOTAL": TN, _, argIdxs = pvTP(self._musUniqueCN, self.N, self.polybasis0, self._derIdxs, self._reorder, scl = np.power(self.scaleFactor, -1.)) TN = TN[:, argIdxs] else: #if self.polydegreetype == "FULL": TN = pvP(self._musUniqueCN, [self.N] * self.npar, self.polybasis0, self._derIdxs, self._reorder, scl = np.power(self.scaleFactor, -1.)) invD = [None] * (len(idxsB)) for k in range(len(idxsB)): pseudoInv = np.diag(fitinv[k, :]) idxGlob = 0 for j, derIdxs in enumerate(self._derIdxs): nder = len(derIdxs) idxGlob += nder if nder > 1: idxLoc = np.arange(len(self.mus))[ (self._reorder >= idxGlob - nder) * (self._reorder < idxGlob)] invLoc = fitinv[k, idxLoc] pseudoInv[np.ix_(idxLoc, idxLoc)] = 0. for diffj, diffjIdx in enumerate(derIdxs): for derQ, derQIdx in enumerate(derIdxs): derUIdx = [x - y for (x, y) in zip(diffjIdx, derQIdx)] if all([x >= 0 for x in derUIdx]): derU = hashD(derUIdx) pseudoInv[idxLoc[derU], idxLoc[derQ]] = ( invLoc[diffj]) invD[k] = pseudoInv.dot(TN) return invD, fitinv def findeveVGExplicit(self, sampleE:sampList, invD:List[Np2D]) -> Tuple[Np1D, Np2D]: """ Compute explicitly eigenvalues and eigenvectors of rational denominator matrix. """ RROMPyAssert(self._mode, message = "Cannot solve eigenvalue problem.") nEnd = invD[0].shape[1] eWidth = len(invD) vbMng(self, "INIT", "Building gramian matrix.", 10) gramian = self.HFEngine.innerProduct(sampleE, sampleE) G = np.zeros((nEnd, nEnd), dtype = np.complex) for k in range(eWidth): G += invD[k].T.conj().dot(gramian.dot(invD[k])) vbMng(self, "DEL", "Done building gramian.", 10) vbMng(self, "INIT", "Solving eigenvalue problem for gramian matrix.", 7) ev, eV = np.linalg.eigh(G) vbMng(self, "MAIN", ("Solved eigenvalue problem of size {} with condition number " "{:.4e}.").format(nEnd, ev[-1] / ev[0]), 5) vbMng(self, "DEL", "Done solving eigenvalue problem.", 7) return ev, eV def findeveVGQR(self, RPODE:Np2D, invD:List[Np2D]) -> Tuple[Np1D, Np2D]: """ Compute eigenvalues and eigenvectors of rational denominator matrix through SVD of R factor. """ RROMPyAssert(self._mode, message = "Cannot solve eigenvalue problem.") nEnd = invD[0].shape[1] S = RPODE.shape[0] eWidth = len(invD) vbMng(self, "INIT", "Building half-gramian matrix stack.", 10) Rstack = np.zeros((S * eWidth, nEnd), dtype = np.complex) for k in range(eWidth): Rstack[k * S : (k + 1) * S, :] = RPODE.dot(invD[k]) vbMng(self, "DEL", "Done building half-gramian.", 10) vbMng(self, "INIT", "Solving svd for square root of gramian matrix.", 7) _, s, eV = np.linalg.svd(Rstack, full_matrices = False) ev = s[::-1] eV = eV[::-1, :].T.conj() vbMng(self, "MAIN", ("Solved svd problem of size {} x {} with condition number " "{:.4e}.").format(*Rstack.shape, s[0] / s[-1]), 5) vbMng(self, "DEL", "Done solving svd.", 7) return ev, eV def getResidues(self, *args, **kwargs) -> Np1D: """ Obtain approximant residues. Returns: Matrix with residues as columns. """ return self.trainedModel.getResidues(*args, **kwargs) diff --git a/rrompy/reduction_methods/standard/rational_moving_least_squares.py b/rrompy/reduction_methods/standard/rational_moving_least_squares.py new file mode 100644 index 0000000..63b24ef --- /dev/null +++ b/rrompy/reduction_methods/standard/rational_moving_least_squares.py @@ -0,0 +1,300 @@ +# Copyright (C) 2018 by the RROMPy authors +# +# This file is part of RROMPy. +# +# RROMPy is free software: you can redistribute it and/or modify +# it under the terms of the GNU Lesser General Public License as published by +# the Free Software Foundation, either version 3 of the License, or +# (at your option) any later version. +# +# RROMPy is distributed in the hope that it will be useful, +# but WITHOUT ANY WARRANTY; without even the implied warranty of +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +# GNU Lesser General Public License for more details. +# +# You should have received a copy of the GNU Lesser General Public License +# along with RROMPy. If not, see . +# + +from copy import deepcopy as copy +import numpy as np +from .rational_interpolant import RationalInterpolant +from rrompy.utilities.poly_fitting.polynomial import (polybases as ppb, + polyvander as pvP, + polyvanderTotal as pvTP) +from rrompy.reduction_methods.trained_model import ( + TrainedModelRationalMLS as tModel) +from rrompy.reduction_methods.trained_model import TrainedModelData +from rrompy.utilities.base.types import Np2D, HFEng, DictAny, paramVal +from rrompy.utilities.base import verbosityManager as vbMng +from rrompy.utilities.numerical import fullDegreeMaxMask, totalDegreeMaxMask +from rrompy.utilities.exception_manager import (RROMPyException, RROMPyAssert, + RROMPyWarning) + +__all__ = ['RationalMovingLeastSquares'] + +class RationalMovingLeastSquares(RationalInterpolant): + """ + ROM rational moving LS interpolant computation for parametric problems. + + Args: + HFEngine: HF problem solver. + mu0(optional): Default parameter. Defaults to 0. + approxParameters(optional): Dictionary containing values for main + parameters of approximant. Recognized keys are: + - 'POD': whether to compute POD of snapshots; defaults to True; + - 'S': total number of samples current approximant relies upon; + - 'sampler': sample point generator; + - 'polybasis': type of polynomial basis for interpolation; defaults + to 'MONOMIAL'; + - 'M': degree of rational interpolant numerator; defaults to 0; + - 'N': degree of rational interpolant denominator; defaults to 0; + - 'polydegreetype': type of polynomial degree; defaults to 'TOTAL'; + - 'radialBasis': numerator radial basis type; defaults to + 'GAUSSIAN'; + - 'radialDirectionalWeights': radial basis weights for interpolant + numerator; defaults to 0, i.e. identity; + - 'nNearestNeighbor': number of nearest neighbors considered in + numerator if radialBasis allows; defaults to -1; + - 'radialBasisDen': denominator radial basis type; defaults to + 'GAUSSIAN'; + - 'radialDirectionalWeightsDen': radial basis weights for + interpolant denominator; defaults to 0, i.e. identity; + - 'nNearestNeighborDen': number of nearest neighbors considered in + denominator if radialBasisDen allows; defaults to -1; + - 'interpRcond': tolerance for interpolation; defaults to None; + - 'robustTol': 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. + mus: Array of snapshot parameters. + approxParameters: Dictionary containing values for main parameters of + approximant. Recognized keys are in parameterList. + parameterListSoft: Recognized keys of soft approximant parameters: + - 'POD': whether to compute POD of snapshots; + - 'polybasis': type of polynomial basis for interpolation; + - 'M': degree of rational interpolant numerator; + - 'N': degree of rational interpolant denominator; + - 'polydegreetype': type of polynomial degree; + - 'radialBasis': numerator radial basis type; + - 'radialDirectionalWeights': radial basis weights for interpolant + numerator; + - 'nNearestNeighbor': number of nearest neighbors considered in + numerator if radialBasis allows; + - 'radialBasisDen': denominator radial basis type; + - 'radialDirectionalWeightsDen': radial basis weights for + interpolant denominator; + - 'nNearestNeighborDen': number of nearest neighbors considered in + denominator if radialBasisDen allows; + - 'interpRcond': tolerance for interpolation via numpy.polyfit; + - 'robustTol': tolerance for robust rational denominator + management. + parameterListCritical: Recognized keys of critical approximant + parameters: + - 'S': total number of samples current approximant relies upon; + - 'sampler': sample point generator. + POD: Whether to compute POD of snapshots. + S: Number of solution snapshots over which current approximant is + based upon. + sampler: Sample point generator. + polybasis: type of polynomial basis for interpolation. + M: Numerator degree of approximant. + N: Denominator degree of approximant. + polydegreetype: Type of polynomial degree. + radialBasis: Numerator radial basis type. + radialDirectionalWeights: Radial basis weights for interpolant + numerator. + nNearestNeighbor: Number of nearest neighbors considered in numerator + if radialBasis allows. + radialBasisDen: Denominator radial basis type. + radialDirectionalWeightsDen: Radial basis weights for interpolant + denominator. + nNearestNeighborDen: Number of nearest neighbors considered in + denominator if radialBasisDen allows. + interpRcond: Tolerance for interpolation via numpy.polyfit. + robustTol: Tolerance for robust rational denominator management. + muBounds: list of bounds for parameter values. + samplingEngine: Sampling engine. + uHF: High fidelity solution(s) with parameter(s) lastSolvedHF as + sampleList. + lastSolvedHF: Parameter(s) corresponding to last computed high fidelity + solution(s) as parameterList. + uApproxReduced: Reduced approximate solution(s) with parameter(s) + lastSolvedApprox as sampleList. + lastSolvedApproxReduced: Parameter(s) corresponding to last computed + reduced approximate solution(s) as parameterList. + uApprox: Approximate solution(s) with parameter(s) lastSolvedApprox as + sampleList. + lastSolvedApprox: Parameter(s) corresponding to last computed + approximate solution(s) as parameterList. + Q: Numpy 1D vector containing complex coefficients of approximant + denominator. + P: Numpy 2D vector whose columns are FE dofs of coefficients of + approximant numerator. + """ + + def __init__(self, HFEngine:HFEng, mu0 : paramVal = None, + approxParameters : DictAny = {}, verbosity : int = 10, + timestamp : bool = True): + self._preInit() + self._addParametersToList(["radialBasis", "radialBasisDen", + "radialDirectionalWeightsDen", + "nNearestNeighborDen"], + ["GAUSSIAN", "GAUSSIAN", 1, -1]) + super().__init__(HFEngine = HFEngine, mu0 = mu0, + approxParameters = approxParameters, + verbosity = verbosity, timestamp = timestamp) + self.catchInstability = False + self._postInit() + + @property + def polybasis(self): + """Value of polybasis.""" + return self._polybasis + @polybasis.setter + def polybasis(self, polybasis): + try: + polybasis = polybasis.upper().strip().replace(" ","") + if polybasis not in ppb: + raise RROMPyException("Prescribed polybasis not recognized.") + self._polybasis = polybasis + except: + RROMPyWarning(("Prescribed polybasis not recognized. Overriding " + "to 'MONOMIAL'.")) + self._polybasis = "MONOMIAL" + self._approxParameters["polybasis"] = self.polybasis + + @property + def radialBasis(self): + """Value of radialBasis.""" + return self._radialBasis + @radialBasis.setter + def radialBasis(self, radialBasis): + self._radialBasis = radialBasis + self._approxParameters["radialBasis"] = self.radialBasis + + @property + def radialBasisDen(self): + """Value of radialBasisDen.""" + return self._radialBasisDen + @radialBasisDen.setter + def radialBasisDen(self, radialBasisDen): + self._radialBasisDen = radialBasisDen + self._approxParameters["radialBasisDen"] = self.radialBasisDen + + @property + def radialDirectionalWeightsDen(self): + """Value of radialDirectionalWeightsDen.""" + return self._radialDirectionalWeightsDen + @radialDirectionalWeightsDen.setter + def radialDirectionalWeightsDen(self, radialDirectionalWeightsDen): + self._radialDirectionalWeightsDen = radialDirectionalWeightsDen + self._approxParameters["radialDirectionalWeightsDen"] = ( + self.radialDirectionalWeightsDen) + + @property + def nNearestNeighborDen(self): + """Value of nNearestNeighborDen.""" + return self._nNearestNeighborDen + @nNearestNeighborDen.setter + def nNearestNeighborDen(self, nNearestNeighborDen): + self._nNearestNeighborDen = nNearestNeighborDen + self._approxParameters["nNearestNeighborDen"] = ( + self.nNearestNeighborDen) + + def _setupDenominator(self) -> Np2D: + """Compute rational denominator.""" + RROMPyAssert(self._mode, message = "Cannot setup denominator.") + vbMng(self, "INIT", + "Starting computation of denominator-related blocks.", 7) + self._setupInterpolationIndices() + if self.polydegreetype == "TOTAL": + TN, _, argIdxs = pvTP(self._musUniqueCN, self.N, self.polybasis0, + self._derIdxs, self._reorder, + scl = np.power(self.scaleFactor, -1.)) + TN = TN[:, argIdxs] + else: #if self.polydegreetype == "FULL": + TN = pvP(self._musUniqueCN, [self.N] * self.npar, + self.polybasis0, self._derIdxs, self._reorder, + scl = np.power(self.scaleFactor, -1.)) + TNTen = np.zeros((self.S, self.S, TN.shape[1]), dtype = TN.dtype) + TNTen[np.arange(self.S), np.arange(self.S)] = TN + if self.POD: TNTen = np.tensordot(self.samplingEngine.RPOD, TNTen, 1) + vbMng(self, "DEL", "Done computing denominator-related blocks.", 7) + return TN, TNTen + + def _setupNumerator(self) -> Np2D: + """Compute rational numerator.""" + RROMPyAssert(self._mode, message = "Cannot setup numerator.") + vbMng(self, "INIT", + "Starting computation of denominator-related blocks.", 7) + self._setupInterpolationIndices() + if self.polydegreetype == "TOTAL": + TM, _, argIdxs = pvTP(self._musUniqueCN, self.M, self.polybasis0, + self._derIdxs, self._reorder, + scl = np.power(self.scaleFactor, -1.)) + TM = TM[:, argIdxs] + else: #if self.polydegreetype == "FULL": + TM = pvP(self._musUniqueCN, [self.M] * self.npar, + self.polybasis0, self._derIdxs, self._reorder, + scl = np.power(self.scaleFactor, -1.)) + vbMng(self, "DEL", "Done computing denominator-related blocks.", 7) + return TM + + def setupApprox(self): + """ + Compute rational interpolant. + SVD-based robust eigenvalue management. + """ + if self.checkComputedApprox(): + return + RROMPyAssert(self._mode, message = "Cannot setup approximant.") + vbMng(self, "INIT", "Setting up {}.". format(self.name()), 5) + self.computeSnapshots() + if self.trainedModel is None: + self.trainedModel = tModel() + self.trainedModel.verbosity = self.verbosity + self.trainedModel.timestamp = self.timestamp + data = TrainedModelData(self.trainedModel.name(), self.mu0, + self.samplingEngine.samples, + self.scaleFactor, + self.HFEngine.rescalingExp) + data.POD = self.POD + data.polybasis = self.polybasis + data.polydegreetype = self.polydegreetype + data.radialBasis = self.radialBasis + data.radialWeights = self.radialDirectionalWeights + data.nNearestNeighbor = self.nNearestNeighbor + data.radialBasisDen = self.radialBasisDen + data.radialWeightsDen = self.radialDirectionalWeightsDen + data.nNearestNeighborDen = self.nNearestNeighborDen + data.interpRcond = self.interpRcond + self.trainedModel.data = data + else: + self.trainedModel = self.trainedModel + self.trainedModel.data.projMat = copy(self.samplingEngine.samples) + if not self.POD: + self.trainedModel.data.gramian = self.HFEngine.innerProduct( + self.samplingEngine.samples, + self.samplingEngine.samples) + self.trainedModel.data.mus = copy(self.mus) + self.trainedModel.data.M = self.M + self.trainedModel.data.N = self.N + QVan, self.trainedModel.data.QBlocks = self._setupDenominator() + self.trainedModel.data.PVan = self._setupNumerator() + if self.polydegreetype == "TOTAL": + degreeMaxMask = totalDegreeMaxMask + else: #if self.polydegreetype == "FULL": + degreeMaxMask = fullDegreeMaxMask + if self.N > self.M: + self.trainedModel.data.QVan = QVan + self.trainedModel.data.domQIdxs = degreeMaxMask(self.N, self.npar) + else: + self.trainedModel.data.domQIdxs = degreeMaxMask(self.M, self.npar) + self.trainedModel.data.approxParameters = copy(self.approxParameters) + vbMng(self, "DEL", "Done setting up approximant.", 5) + diff --git a/rrompy/reduction_methods/trained_model/__init__.py b/rrompy/reduction_methods/trained_model/__init__.py index 1e8d8b4..ca0c25b 100644 --- a/rrompy/reduction_methods/trained_model/__init__.py +++ b/rrompy/reduction_methods/trained_model/__init__.py @@ -1,33 +1,35 @@ # Copyright (C) 2018 by the RROMPy authors # # This file is part of RROMPy. # # RROMPy is free software: you can redistribute it and/or modify # it under the terms of the GNU Lesser General Public License as published by # the Free Software Foundation, either version 3 of the License, or # (at your option) any later version. # # RROMPy is distributed in the hope that it will be useful, # but WITHOUT ANY WARRANTY; without even the implied warranty of # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the # GNU Lesser General Public License for more details. # # You should have received a copy of the GNU Lesser General Public License # along with RROMPy. If not, see . # from .trained_model_data import TrainedModelData from .trained_model_rational import TrainedModelRational +from .trained_model_rational_mls import TrainedModelRationalMLS from .trained_model_reduced_basis import TrainedModelReducedBasis from .trained_model_pivoted_data import TrainedModelPivotedData from .trained_model_pivoted_rational import TrainedModelPivotedRational __all__ = [ 'TrainedModelData', 'TrainedModelRational', + 'TrainedModelRationalMLS', 'TrainedModelReducedBasis', 'TrainedModelPivotedData', 'TrainedModelPivotedRational' ] diff --git a/rrompy/reduction_methods/trained_model/trained_model.py b/rrompy/reduction_methods/trained_model/trained_model.py index dfc80ee..0c6a4bb 100644 --- a/rrompy/reduction_methods/trained_model/trained_model.py +++ b/rrompy/reduction_methods/trained_model/trained_model.py @@ -1,89 +1,96 @@ # Copyright (C) 2018 by the RROMPy authors # # This file is part of RROMPy. # # RROMPy is free software: you can redistribute it and/or modify # it under the terms of the GNU Lesser General Public License as published by # the Free Software Foundation, either version 3 of the License, or # (at your option) any later version. # # RROMPy is distributed in the hope that it will be useful, # but WITHOUT ANY WARRANTY; without even the implied warranty of # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the # GNU Lesser General Public License for more details. # # You should have received a copy of the GNU Lesser General Public License # along with RROMPy. If not, see . # from abc import abstractmethod from rrompy.utilities.base.types import Np1D, paramList, sampList from rrompy.parameter import checkParameterList from rrompy.sampling import sampleList, emptySampleList __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)) + @property + def npar(self): + """Number of parameters.""" + return self.data.mu0.shape[1] + @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 = checkParameterList(mu, self.data.npar)[0] if (not hasattr(self, "lastSolvedApprox") or self.lastSolvedApprox != mu): uApproxR = self.getApproxReduced(mu) self.uApprox = emptySampleList() - self.uApprox.reset((self.data.projMat.shape[0], len(mu)), - self.data.projMat.dtype) for i in range(len(mu)): if isinstance(self.data.projMat, (list, sampleList,)): - self.uApprox[i] = uApproxR[i][0] * self.data.projMat[0] + uApp = uApproxR[i][0] * self.data.projMat[0] for j in range(1, uApproxR.shape[0]): - self.uApprox[i] += (uApproxR[i][j] - * self.data.projMat[j]) + uApp += uApproxR[i][j] * self.data.projMat[j] else: - self.uApprox[i] = self.data.projMat.dot(uApproxR[i]) + uApp = self.data.projMat.dot(uApproxR[i]) + if i == 0: + #self.data.projMat.shape[0], len(mu) + self.uApprox.reset((len(uApp), len(mu)), + dtype = uApp.dtype) + self.uApprox[i] = uApp self.lastSolvedApprox = mu return self.uApprox @abstractmethod def getPoles(self) -> Np1D: """ Obtain approximant poles. Returns: Numpy complex vector of poles. """ pass diff --git a/rrompy/reduction_methods/trained_model/trained_model_pivoted_general.py b/rrompy/reduction_methods/trained_model/trained_model_pivoted_general.py index af04fb4..96b219a 100644 --- a/rrompy/reduction_methods/trained_model/trained_model_pivoted_general.py +++ b/rrompy/reduction_methods/trained_model/trained_model_pivoted_general.py @@ -1,373 +1,375 @@ # Copyright (C) 2018 by the RROMPy authors # # This file is part of RROMPy. # # RROMPy is free software: you can redistribute it and/or modify # it under the terms of the GNU Lesser General Public License as published by # the Free Software Foundation, either version 3 of the License, or # (at your option) any later version. # # RROMPy is distributed in the hope that it will be useful, # but WITHOUT ANY WARRANTY; without even the implied warranty of # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the # GNU Lesser General Public License for more details. # # You should have received a copy of the GNU Lesser General Public License # along with RROMPy. If not, see . # import numpy as np from scipy.special import factorial as fact from itertools import combinations from .trained_model import TrainedModel from rrompy.utilities.base.types import (Np1D, Tuple, List, ListAny, paramVal, paramList, sampList, HFEng) from rrompy.utilities.base import verbosityManager as vbMng, freepar as fp from rrompy.utilities.numerical import pointMatching from rrompy.utilities.poly_fitting.heaviside import HeavisideInterpolator as HI from rrompy.utilities.exception_manager import RROMPyException, RROMPyAssert from rrompy.parameter import checkParameter, checkParameterList from rrompy.sampling import emptySampleList, sampleList __all__ = ['TrainedModelPivotedGeneral'] class TrainedModelPivotedGeneral(TrainedModel): """ ROM approximant evaluation for pivoted approximants (with pole matching). Attributes: Data: dictionary with all that can be pickled. """ def centerNormalizePivot(self, mu : paramList = [], mu0 : paramVal = None) -> paramList: """ Compute normalized parameter to be plugged into approximant. Args: mu: Parameter(s) 1. mu0: Parameter(s) 2. If None, set to self.data.mu0Pivot. Returns: Normalized parameter. """ mu = checkParameterList(mu, self.data.nparPivot)[0] if mu0 is None: mu0 = self.data.mu0Pivot rad = ((mu ** self.data.rescalingExpPivot - mu0 ** self.data.rescalingExpPivot) / self.data.scaleFactorPivot) return rad def centerNormalizeMarginal(self, mu : paramList = [], mu0 : paramVal = None) -> paramList: """ Compute normalized parameter to be plugged into approximant. Args: mu: Parameter(s) 1. mu0: Parameter(s) 2. If None, set to self.data.mu0Marginal. Returns: Normalized parameter. """ mu = checkParameterList(mu, self.data.nparMarginal)[0] if mu0 is None: mu0 = self.data.mu0Marginal rad = ((mu ** self.data.rescalingExpMarginal - mu0 ** self.data.rescalingExpMarginal) / self.data.scaleFactorMarginal) return rad def initializeFromLists(self, poles:ListAny, coeffs:ListAny, basis:str, HFEngine:HFEng, matchingWeight : float = 1., POD : bool = True): """Initialize Heaviside representation.""" musM = self.data.musMarginal margAbsDist = np.sum(np.abs(np.repeat(musM.data, len(musM), 0) - np.tile(musM.data, [len(musM), 1]) ), axis = 1).reshape(len(musM), len(musM)) N = len(poles[0]) explored = [0] unexplored = list(range(1, len(musM))) for _ in range(1, len(musM)): minIdx = np.argmin(np.concatenate([margAbsDist[ex, unexplored] \ for ex in explored])) minIex = explored[minIdx // len(unexplored)] minIunex = unexplored[minIdx % len(unexplored)] dist = np.abs(np.tile(poles[minIex].reshape(-1, 1), N) - poles[minIunex].reshape(1, -1)) if matchingWeight != 0: resex = coeffs[minIex][: N] resunex = coeffs[minIunex][: N] if POD: distR = resex.dot(resunex.T.conj()) distR = (distR.T / np.linalg.norm(resex, axis = 1)).T distR = distR / np.linalg.norm(resunex, axis = 1) else: resex = self.data.projMat.dot(resex.T) resunex = self.data.projMat.dot(resunex.T) distR = HFEngine.innerProduct(resex, resunex).T distR = (distR.T / HFEngine.norm(resex)).T distR = distR / HFEngine.norm(resunex) distR = np.abs(distR) distR[distR > 1.] = 1. dist += 2. / np.pi * matchingWeight * np.arccos(distR) reordering = pointMatching(dist) poles[minIunex] = poles[minIunex][reordering] coeffs[minIunex][: N] = coeffs[minIunex][reordering] explored += [minIunex] unexplored.remove(minIunex) HIs = [] for pls, cfs in zip(poles, coeffs): hsi = HI() hsi.poles = pls hsi.coeffs = cfs hsi.npar = 1 hsi.polybasis = basis HIs += [hsi] self.data.HIs = HIs def recompressByCutOff(self, murange : Tuple[float, float] = [- 1., 1.], tol : float = np.inf, rtype : str = "MAGNITUDE"): if np.isinf(tol): return " No poles erased." N = len(self.data.HIs[0].poles) mu0 = np.mean(murange) musig = murange[0] - mu0 if np.isclose(musig, 0.): radius = lambda x: np.abs(x - mu0) else: if rtype == "MAGNITUDE": murdir = (murange[0] - mu0) / np.abs(musig) def radius(x): scalprod = (x - mu0) * murdir.conj() / np.abs(musig) rescalepar = np.abs(np.real(scalprod)) rescaleort = np.abs(np.imag(scalprod)) return ((rescalepar - 1.) ** 2. * (rescalepar > 1.) + rescaleort ** 2.) ** .5 else:#if rtype == "POTENTIAL": def radius(x): rescale = (x - mu0) / musig return np.max(np.abs(rescale * np.array([-1., 1.]) + (rescale ** 2. - 1) ** .5)) - 1. keepPole, removePole = [], [] for j in range(N): for hi in self.data.HIs: if radius(hi.poles[j]) <= tol: keepPole += [j] break if len(keepPole) == 0 or keepPole[-1] != j: removePole += [j] if len(keepPole) == N: return " No poles erased." keepCoeff = keepPole + [N] keepCoeff = keepCoeff + list(range(N + 1,len(self.data.HIs[0].coeffs))) for hi in self.data.HIs: polyCorrection = np.zeros_like(hi.coeffs[0, :]) for j in removePole: polyCorrection += hi.coeffs[j, :] / (mu0 - hi.poles[j]) if len(hi.coeffs) == N: hi.coeffs = np.vstack((hi.coeffs, polyCorrection)) else: hi.coeffs[N, :] += polyCorrection hi.poles = hi.poles[keepPole] hi.coeffs = hi.coeffs[keepCoeff, :] return (" Erased {} pole".format(len(removePole)) + "s" * (len(removePole) > 1) + ".") def interpolateMarginal(self, mu : paramList = [], samples : ListAny = [], der : List[int] = None, scl : Np1D = None) -> sampList: """ Evaluate marginal interpolator at arbitrary marginal parameter. Args: mu: Target parameter. samples: Objects to interpolate. der(optional): Derivatives to take before evaluation. """ mu = checkParameterList(mu, self.data.nparMarginal)[0] sList = isinstance(samples[0], sampleList) sEff = [None] * len(samples) for j in range(len(samples)): if sList: sEff[j] = samples[j].data else: sEff[j] = samples[j] try: dtype = sEff[0].dtype except: dtype = sEff[0][0].dtype vbMng(self, "INIT", "Interpolating marginal at mu = {}.".format(mu), 95) muC = self.centerNormalizeMarginal(mu) p = emptySampleList() p.reset((len(sEff[0]), len(muC)), dtype = dtype) p.data[:] = 0. if len(sEff[0]) > 0: for mIj, spj in zip(self.data.marginalInterp, sEff): p = p + spj.reshape(len(sEff[0]), - 1) * mIj(muC, der, scl) vbMng(self, "DEL", "Done interpolating marginal.", 95) if not sList: p = p.data.flatten() return p def interpolateMarginalInterpolator(self, mu : paramVal = []) -> Np1D: """Obtain interpolated approximant interpolator.""" mu = checkParameter(mu, self.data.nparMarginal)[0] hsi = HI() hsi.poles = self.interpolateMarginalPoles(mu) hsi.coeffs = self.interpolateMarginalCoeffs(mu) hsi.npar = 1 hsi.polybasis = self.data.HIs[0].polybasis return hsi def interpolateMarginalPoles(self, mu : paramList = []) -> Np1D: """Obtain interpolated approximant poles.""" mu = checkParameterList(mu, self.data.nparMarginal)[0] return self.interpolateMarginal(mu, [hi.poles for hi in self.data.HIs]) def interpolateMarginalCoeffs(self, mu : paramList = []) -> Np1D: """Obtain interpolated approximant coefficients.""" mu = checkParameterList(mu, self.data.nparMarginal)[0] cs = self.interpolateMarginal(mu, [hi.coeffs for hi in self.data.HIs]) if isinstance(cs, (list, tuple,)): cs = np.array(cs) return cs.reshape(self.data.HIs[0].coeffs.shape) def getApproxReduced(self, mu : paramList = []) -> sampList: """ Evaluate reduced representation of approximant at arbitrary parameter. Args: mu: Target parameter. """ RROMPyAssert(self.data.nparPivot, 1, "Number of pivot parameters") mu = checkParameterList(mu, self.data.npar)[0] if (not hasattr(self, "lastSolvedApproxReduced") or self.lastSolvedApproxReduced != mu): vbMng(self, "INIT", "Evaluating approximant at mu = {}.".format(mu), 12) self.uApproxReduced = emptySampleList() - self.uApproxReduced.reset((self.data.HIs[0].coeffs.shape[1], - len(mu)), self.data.projMat.dtype) - for i, muPL in enumerate(mu): muL = self.centerNormalizePivot([muPL(0, x) \ for x in self.data.directionPivot]) muM = [muPL(0, x) for x in self.data.directionMarginal] vbMng(self, "INIT", "Assembling reduced model for mu = {}.".format(muPL), 87) hsL = self.interpolateMarginalInterpolator(muM) vbMng(self, "DEL", "Done assembling reduced model.", 87) - self.uApproxReduced[i] = hsL(muL) + uAppR = hsL(muL) + if i == 0: + #self.data.HIs[0].coeffs.shape[1], len(mu) + self.uApproxReduced.reset((len(uAppR), len(mu)), + dtype = uAppR.dtype) + self.uApproxReduced[i] = uAppR vbMng(self, "DEL", "Done evaluating approximant.", 12) self.lastSolvedApproxReduced = mu return self.uApproxReduced def getPVal(self, mu : paramList = []) -> sampList: """ Evaluate rational numerator at arbitrary parameter. Args: mu: Target parameter. """ RROMPyAssert(self.data.nparPivot, 1, "Number of pivot parameters") mu = checkParameterList(mu, self.data.npar)[0] p = emptySampleList() p.reset((len(self.data.HIs[0].coeffs.shape[1]), len(mu))) for i, muPL in enumerate(mu): muL = self.centerNormalizePivot([muPL(0, x) \ for x in self.data.directionPivot]) muM = [muPL(0, x) for x in self.data.directionMarginal] hsL = self.interpolateMarginalInterpolator(muM) p[i] = hsL(muL) * np.prod(muL(0, 0) - hsL.poles) return p def getQVal(self, mu:Np1D, der : List[int] = None, scl : Np1D = None) -> Np1D: """ Evaluate rational denominator at arbitrary parameter. Args: mu: Target parameter. der(optional): Derivatives to take before evaluation. """ RROMPyAssert(self.data.nparPivot, 1, "Number of pivot parameters") mu = checkParameterList(mu, self.data.npar)[0] muP = self.centerNormalizePivot(checkParameterList( mu.data[:, self.data.directionPivot], self.data.nparPivot)[0]) muM = checkParameterList(mu.data[:, self.data.directionMarginal], self.data.nparMarginal)[0] if der is None: derP, derM = 0, [0] else: derP = der[self.data.directionPivot[0]] derM = [der[x] for x in self.data.directionMarginal] if np.any(np.array(derM) != 0): raise RROMPyException(("Derivatives of Q with respect to marginal " "parameters not allowed.")) sclP = 1 if scl is None else scl[self.data.directionPivot[0]] - derVal = np.zeros(len(mu), dtype = self.data.HIs[0].poles.dtype) + derVal = np.zeros(len(mu), dtype = np.complex) N = len(self.data.HIs[0].poles) if derP == N: derVal[:] = 1. elif derP >= 0 and derP < N: pls = self.interpolateMarginalPoles(muM).reshape(-1, len(mu)).T plsDist = muP.data.reshape(-1, 1) - pls for terms in combinations(np.arange(N), N - derP): derVal += np.prod(plsDist[:, list(terms)], axis = 1) return sclP ** derP * fact(derP) * derVal def getPoles(self, *args, **kwargs) -> Np1D: """ Obtain approximant poles. Returns: Numpy complex vector of poles. """ RROMPyAssert(self.data.nparPivot, 1, "Number of pivot parameters") if len(args) + len(kwargs) > 1: raise RROMPyException(("Wrong number of parameters passed. " "Only 1 available.")) elif len(args) + len(kwargs) == 1: if len(args) == 1: mVals = args[0] else: mVals = kwargs["marginalVals"] if not hasattr(mVals, "__len__"): mVals = [mVals] mVals = list(mVals) else: mVals = [fp] try: rDim = mVals.index(fp) if rDim < len(mVals) - 1 and fp in mVals[rDim + 1 :]: raise except: raise RROMPyException(("Exactly 1 'freepar' entry in " "marginalVals must be provided.")) if rDim != self.data.directionPivot[0]: raise RROMPyException(("'freepar' entry in marginalVals must " "coincide with pivot direction.")) mVals[rDim] = self.data.mu0(rDim) mMarg = [mVals[j] for j in range(len(mVals)) if j != rDim] roots = np.array(self.interpolateMarginalPoles(mMarg)) return np.power(self.data.mu0(rDim) ** self.data.rescalingExp[rDim] + self.data.scaleFactor[rDim] * roots, 1. / self.data.rescalingExp[rDim]) def getResidues(self, *args, **kwargs) -> Np1D: """ Obtain approximant residues. Returns: Numpy matrix with residues as columns. """ pls = self.getPoles(*args, **kwargs) if len(args) == 1: mVals = args[0] else: mVals = kwargs["marginalVals"] if not hasattr(mVals, "__len__"): mVals = [mVals] mVals = list(mVals) rDim = mVals.index(fp) mMarg = [mVals[j] for j in range(len(mVals)) if j != rDim] residues = self.interpolateMarginalCoeffs(mMarg)[: len(pls)] res = self.data.projMat.dot(residues.T) return pls, res diff --git a/rrompy/reduction_methods/trained_model/trained_model_rational_mls.py b/rrompy/reduction_methods/trained_model/trained_model_rational_mls.py new file mode 100644 index 0000000..f012a55 --- /dev/null +++ b/rrompy/reduction_methods/trained_model/trained_model_rational_mls.py @@ -0,0 +1,174 @@ +# Copyright (C) 2018 by the RROMPy authors +# +# This file is part of RROMPy. +# +# RROMPy is free software: you can redistribute it and/or modify +# it under the terms of the GNU Lesser General Public License as published by +# the Free Software Foundation, either version 3 of the License, or +# (at your option) any later version. +# +# RROMPy is distributed in the hope that it will be useful, +# but WITHOUT ANY WARRANTY; without even the implied warranty of +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +# GNU Lesser General Public License for more details. +# +# You should have received a copy of the GNU Lesser General Public License +# along with RROMPy. If not, see . +# + +import numpy as np +from .trained_model_rational import TrainedModelRational +from rrompy.utilities.base.types import Np1D, paramVal, paramList, sampList +from rrompy.utilities.base import verbosityManager as vbMng +from rrompy.utilities.poly_fitting.moving_least_squares import mlsweights +from rrompy.utilities.poly_fitting.polynomial import ( + PolynomialInterpolator as PI) +from rrompy.utilities.numerical import customPInv, degreeTotalToFull +from rrompy.parameter import checkParameterList +from rrompy.sampling import emptySampleList + +__all__ = ['TrainedModelRationalMLS'] + +class TrainedModelRationalMLS(TrainedModelRational): + """ + ROM approximant evaluation for rational moving least squares approximant. + + Attributes: + Data: dictionary with all that can be pickled. + """ + + def assembleReducedModel(self, mu:paramVal): + if not hasattr(self, "lastSetupMu") or self.lastSetupMu != mu: + vbMng(self, "INIT", "Assembling reduced model for mu = {}."\ + .format(mu), 17) + vbMng(self, "INIT", "Starting computation of denominator.", 35) + muC = self.centerNormalize(mu) + muSC = self.centerNormalize(self.data.mus) + wQ = mlsweights(muC, muSC, self.data.radialBasisDen, + directionalWeights = self.data.radialWeightsDen, + nNearestNeighbor = self.data.nNearestNeighborDen) + if self.data.N > self.data.M: + PQVan = self.data.QVan + else: + PQVan = self.data.PVan + VQAdjW = PQVan.conj().T * wQ + VQAdjWVQ = VQAdjW.dot(PQVan) + interpPseudoInverse, info = customPInv(VQAdjWVQ, full = True, + rcond = self.data.interpRcond) + interpPseudoInverse = interpPseudoInverse.dot(VQAdjW).dot( + self.data.QBlocks) + if info[0] < interpPseudoInverse.shape[-1]: + q = np.zeros(interpPseudoInverse.shape[-1], dtype = np.complex) + q[0] = 1. + else: + halfGram = interpPseudoInverse[self.data.domQIdxs] + if self.data.POD: + Rstack = halfGram.reshape(-1, halfGram.shape[-1]) + vbMng(self, "INIT", + "Solving svd for square root of gramian matrix.", 67) + _, s, eV = np.linalg.svd(Rstack, full_matrices = False) + condN = s[0] / s[-1] + q = eV[-1, :].T.conj() + vbMng(self, "MAIN", + ("Solved svd problem of size {} x {} with condition " + "number {:.4e}.").format(*Rstack.shape, condN), 55) + vbMng(self, "DEL", "Done solving svd.", 67) + else: + RRstack = np.tensordot(self.trainedModel.gramian, halfGram, + 1).reshape(-1, halfGram.shape[-1]) + RLstack = halfGram.reshape(-1, halfGram.shape[-1]) + gram = RLstack.T.conj().dot(RRstack) + vbMng(self, "INIT", + "Solving eigenvalue problem for gramian matrix.", 67) + ev, eV = np.linalg.eigh(gram) + condN = ev[-1] / ev[0] + q = eV[:, 0] + vbMng(self, "MAIN", + ("Solved eigenvalue problem of size {} with " + "condition number {:.4e}.").format(gram.shape[0], + condN), 55) + vbMng(self, "DEL", "Done solving eigenvalue problem.", 67) + self.data.Q = PI() + self.data.Q.npar = self.npar + self.data.Q.polybasis = self.data.polybasis + if self.data.polydegreetype == "TOTAL": + self.data.Q.coeffs = degreeTotalToFull( + (self.data.N + 1,) * self.npar, + self.npar, q) + else: + self.data.Q.coeffs = q.reshape((self.data.N + 1,) * self.npar) + vbMng(self, "DEL", "Done computing denominator.", 35) + vbMng(self, "INIT", "Starting computation of numerator.", 35) + self.data.P = PI() + self.data.P.npar = self.npar + self.data.P.polybasis = self.data.polybasis + wP = mlsweights(muC, muSC, self.data.radialBasis, + directionalWeights = self.data.radialWeights, + nNearestNeighbor = self.data.nNearestNeighbor) + VAdjW = self.data.PVan.conj().T * wP + VAdjWV = VAdjW.dot(self.data.PVan) + interpPPseudoInverse = customPInv(VAdjWV, self.data.interpRcond) + Pcoeffs = np.tensordot(interpPPseudoInverse.dot(VAdjW), + self.data.QBlocks.dot(q), ([1], [1])) + if self.data.polydegreetype == "TOTAL": + self.data.P.coeffs = degreeTotalToFull( + (self.data.M + 1,) * self.npar + + (self.data.QBlocks.shape[0],), + self.npar, Pcoeffs) + else: + self.data.P.coeffs = Pcoeffs.reshape( + (self.data.M + 1,) * self.npar + + (self.data.QBlocks.shape[0],)) + vbMng(self, "DEL", "Done computing numerator.", 35) + vbMng(self, "DEL", "Done assembling reduced model.", 17) + self.lastSetupMu = mu + + def getApproxReduced(self, mu : paramList = []) -> sampList: + """ + Evaluate reduced representation of approximant at arbitrary parameter. + + Args: + mu: Target parameter. + """ + mu = checkParameterList(mu, self.data.npar)[0] + if (not hasattr(self, "lastSolvedApproxReduced") + or self.lastSolvedApproxReduced != mu): + vbMng(self, "INIT", + "Evaluating approximant at mu = {}.".format(mu), 12) + self.uApproxReduced = emptySampleList() + for i in range(len(mu)): + self.assembleReducedModel(mu[i]) + vbMng(self, "INIT", + "Solving reduced model for mu = {}.".format(mu[i]), 15) + uAppR = self.getPVal(mu[i]) / self.getQVal(mu[i]) + if i == 0: + #self.data.P.shape[-1], len(mu) + self.uApproxReduced.reset((len(uAppR), len(mu)), + dtype = uAppR.dtype) + self.uApproxReduced[i] = uAppR + vbMng(self, "DEL", "Done solving reduced model.", 15) + vbMng(self, "DEL", "Done evaluating approximant.", 12) + self.lastSolvedApproxReduced = mu + return self.uApproxReduced + + def getPoles(self, *args, mu : paramVal = None, **kwargs) -> Np1D: + """ + Obtain approximant poles. + + Returns: + Numpy complex vector of poles. + """ + if mu is None: mu = self.data.mu0 + self.assembleReducedModel(mu) + return super().getPoles(*args, **kwargs) + + def getResidues(self, *args, mu : paramVal = None, **kwargs) -> Np1D: + """ + Obtain approximant residues. + + Returns: + Numpy matrix with residues as columns. + """ + if mu is None: mu = self.data.mu0 + self.assembleReducedModel(mu) + return super().getResidues(*args, **kwargs) diff --git a/rrompy/reduction_methods/trained_model/trained_model_reduced_basis.py b/rrompy/reduction_methods/trained_model/trained_model_reduced_basis.py index 53a20ff..f3cf511 100644 --- a/rrompy/reduction_methods/trained_model/trained_model_reduced_basis.py +++ b/rrompy/reduction_methods/trained_model/trained_model_reduced_basis.py @@ -1,108 +1,115 @@ # Copyright (C) 2018 by the RROMPy authors # # This file is part of RROMPy. # # RROMPy is free software: you can redistribute it and/or modify # it under the terms of the GNU Lesser General Public License as published by # the Free Software Foundation, either version 3 of the License, or # (at your option) any later version. # # RROMPy is distributed in the hope that it will be useful, # but WITHOUT ANY WARRANTY; without even the implied warranty of # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the # GNU Lesser General Public License for more details. # # You should have received a copy of the GNU Lesser General Public License # along with RROMPy. If not, see . # import numpy as np from .trained_model import TrainedModel -from rrompy.utilities.base.types import Np1D, ListAny, paramList, sampList +from rrompy.utilities.base.types import (Np1D, ListAny, paramVal, paramList, + sampList) from rrompy.utilities.base import verbosityManager as vbMng, freepar as fp from rrompy.utilities.numerical import (eigvalsNonlinearDense, marginalizePolyList) from rrompy.utilities.expression import expressionEvaluator from rrompy.utilities.exception_manager import RROMPyException, RROMPyWarning from rrompy.parameter import checkParameter, checkParameterList from rrompy.sampling import emptySampleList __all__ = ['TrainedModelReducedBasis'] class TrainedModelReducedBasis(TrainedModel): """ ROM approximant evaluation for RB approximant. Attributes: Data: dictionary with all that can be pickled. """ + def assembleReducedModel(self, mu:paramVal): + mu = checkParameter(mu, self.data.npar) + if not hasattr(self, "lastSetupMu") or self.lastSetupMu != mu: + vbMng(self, "INIT", "Assembling reduced model for mu = {}."\ + .format(mu), 17) + muEff = mu ** self.data.rescalingExp + self.data.ARBmu, self.data.bRBmu = 0., 0. + for thA, ARB in zip(self.data.thAs, self.data.ARBs): + self.data.ARBmu = (expressionEvaluator(thA[0], muEff) * ARB + + self.data.ARBmu) + for thb, bRB in zip(self.data.thbs, self.data.bRBs): + self.data.bRBmu = (expressionEvaluator(thb[0], muEff) * bRB + + self.data.bRBmu) + vbMng(self, "DEL", "Done assembling reduced model.", 17) + self.lastSetupMu = mu + def getApproxReduced(self, mu : paramList = []) -> sampList: """ Evaluate reduced representation of approximant at arbitrary parameter. Args: mu: Target parameter. """ mu = checkParameterList(mu, self.data.npar)[0] if (not hasattr(self, "lastSolvedApproxReduced") or self.lastSolvedApproxReduced != mu): vbMng(self, "INIT", "Computing RB solution at mu = {}.".format(mu), 12) - ARBs, bRBs = self.data.ARBs, self.data.bRBs self.uApproxReduced = emptySampleList() - self.uApproxReduced.reset((ARBs[0].shape[0], len(mu)), - self.data.projMat.dtype) - muEff = mu ** self.data.rescalingExp - thAsVals = [expressionEvaluator(thA[0], muEff, (len(mu),)) \ - for thA in self.data.thAs] - thbsVals = [expressionEvaluator(thb[0], muEff, (len(mu),)) \ - for thb in self.data.thbs] for i in range(len(mu)): - vbMng(self, "INIT", "Assembling reduced model for mu = {}."\ - .format(mu[i]), 17) - ARBmu, bRBmu = 0., 0. - for thAv, ARB in zip(thAsVals, ARBs): - ARBmu = thAv[i] * ARB + ARBmu - for thbv, bRB in zip(thbsVals, bRBs): - bRBmu = thbv[i] * bRB + bRBmu - vbMng(self, "DEL", "Done assembling reduced model.", 17) + self.assembleReducedModel(mu[i]) vbMng(self, "INIT", "Solving reduced model for mu = {}.".format(mu[i]), 15) - self.uApproxReduced[i] = np.linalg.solve(ARBmu, bRBmu) + uAppR = np.linalg.solve(self.data.ARBmu, self.data.bRBmu) + if i == 0: + #self.data.ARBs[0].shape[-1], len(mu) + self.uApproxReduced.reset((len(uAppR), len(mu)), + dtype = uAppR.dtype) + self.uApproxReduced[i] = uAppR vbMng(self, "DEL", "Done solving reduced model.", 15) vbMng(self, "DEL", "Done computing RB solution.", 12) self.lastSolvedApproxReduced = mu return self.uApproxReduced def getPoles(self, marginalVals : ListAny = [fp], jSupp : int = 1, **kwargs) -> Np1D: """ Obtain approximant poles. Returns: Numpy complex vector of poles. """ if not self.data.affinePoly: RROMPyWarning(("Unable to compute approximate poles due " "to parametric dependence (detected non-" "polynomial). Change HFEngine.affinePoly to True " "if necessary.")) return if not hasattr(marginalVals, "__len__"): marginalVals = [marginalVals] mVals = list(marginalVals) try: rDim = mVals.index(fp) if rDim < len(mVals) - 1 and fp in mVals[rDim + 1 :]: raise except: raise RROMPyException(("Exactly 1 'freepar' entry in " "marginalVals must be provided.")) ARBs = self.data.ARBs if self.data.npar > 1: mVals[rDim] = self.data.mu0(rDim) mVals = checkParameter(mVals).data.flatten() mVals[rDim] = fp ARBs = marginalizePolyList(ARBs, mVals, "auto") ev = eigvalsNonlinearDense(ARBs, jSupp = jSupp, **kwargs) return np.power(ev, 1. / self.data.rescalingExp[rDim]) diff --git a/rrompy/sampling/sample_list.py b/rrompy/sampling/sample_list.py index 79dc269..f050c88 100644 --- a/rrompy/sampling/sample_list.py +++ b/rrompy/sampling/sample_list.py @@ -1,222 +1,222 @@ # Copyright (C) 2018 by the RROMPy authors # # This file is part of RROMPy. # # RROMPy is free software: you can redistribute it and/or modify # it under the terms of the GNU Lesser General Public License as published by # the Free Software Foundation, either version 3 of the License, or # (at your option) any later version. # # RROMPy is distributed in the hope that it will be useful, # but WITHOUT ANY WARRANTY; without even the implied warranty of # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the # GNU Lesser General Public License for more details. # # You should have received a copy of the GNU Lesser General Public License # along with RROMPy. If not, see . # from copy import deepcopy as copy import numpy as np from rrompy.utilities.exception_manager import RROMPyAssert from rrompy.utilities.base.types import Np1D, List __all__ = ['emptySampleList', 'sampleList'] def emptySampleList(): return sampleList(np.empty((0, 0))) class sampleList: """HERE""" def __init__(self, data:List[Np1D], lengthCheck : int = None, deep : bool = True): if isinstance(data, (self.__class__,)): data = data.data if isinstance(data, (np.ndarray,)): self.data = copy(data) if deep else data if self.data.ndim <= 1: self.data.shape = (self.data.shape[0], 1) else: if not isinstance(data, (list,)): data = [data] self.data = np.empty((len(data[0]), len(data)), dtype = data[0].dtype) for j, par in enumerate(data): self[j] = copy(data[j]) if deep else data[j] if j == 0 and lengthCheck is None: lengthCheck = self.shape[0] RROMPyAssert(len(data[j]), lengthCheck, "Number of parameters") def __len__(self): return self.shape[1] def __str__(self): return str(self.data) def __repr__(self): return repr(self.data) @property def shape(self): return self.data.shape @property def re(self): return sampleList(np.real(self.data)) @property def im(self): return sampleList(np.imag(self.data)) @property def abs(self): return sampleList(np.abs(self.data)) @property def angle(self): return sampleList(np.angle(self.data)) def conj(self): return sampleList(np.conj(self.data)) @property def T(self): return sampleList(self.data.T) @property def H(self): return sampleList(self.data.T.conj()) @property def dtype(self): return self.data.dtype @dtype.setter def dtype(self, dtype): self.data.dtype = dtype def __getitem__(self, key): return self.data[:, key] def __call__(self, key): return sampleList(self.data[:, key]) def __setitem__(self, key, value): if isinstance(value, self.__class__): value = value.data if isinstance(key, (tuple, list,)): RROMPyAssert(len(key), len(value), "Slice length") for k, val in zip(key, value): self[k] = val else: self.data[:, key] = value.flatten() def __iter__(self): return self.data.T.__iter__() def __eq__(self, other): if not hasattr(other, "shape") or self.shape != other.shape: return False if isinstance(other, self.__class__): fac = other.data else: fac = other return np.allclose(self.data, fac) def __ne__(self, other): return not self == other def __copy__(self): return sampleList(self.data) def __deepcopy__(self, memo): return sampleList(copy(self.data, memo)) def __add__(self, other): if isinstance(other, self.__class__): RROMPyAssert(self.shape, other.shape, "Sample shape") fac = other.data else: fac = other return sampleList(self.data + fac) def __iadd__(self, other): self.data = (self + other).data return self def __sub__(self, other): if isinstance(other, self.__class__): RROMPyAssert(self.shape, other.shape, "Sample shape") fac = other.data else: fac = other return sampleList(self.data - fac) def __isub__(self, other): self.data = (self - other).data return self def __mul__(self, other): if isinstance(other, self.__class__): RROMPyAssert(self.shape, other.shape, "Sample shape") fac = other.data else: fac = other return sampleList(self.data * fac) def __imul__(self, other): self.data = (self * other).data return self def __truediv__(self, other): if isinstance(other, self.__class__): RROMPyAssert(self.shape, other.shape, "Sample shape") fac = other.data else: fac = other return sampleList(self.data / fac) def __idiv__(self, other): self.data = (self / other).data return self def __pow__(self, other): if isinstance(other, self.__class__): RROMPyAssert(self.shape, other.shape, "Sample shape") fac = other.data else: fac = other return sampleList(np.power(self.data, fac)) def __ipow__(self, other): self.data = (self ** other).data return self def __neg__(self): return sampleList(- self.data) def __pos__(self): return sampleList(self.data) - def reset(self, size, dtype = np.float): + def reset(self, size, dtype = np.complex): self.data = np.empty(size, dtype = dtype) self.data[:] = np.nan def append(self, items): if isinstance(items, self.__class__): fac = items.data else: fac = np.array(items, ndmin = 2) self.data = np.append(self.data, fac, axis = 1) def pop(self, idx = -1): self.data = np.delete(self.data, idx, axis = 1) def dot(self, other, sampleListOut : bool = True): if isinstance(other, self.__class__): other = other.data prod = self.data.dot(other) if sampleListOut: prod = sampleList(prod) return prod diff --git a/rrompy/solver/norm_utilities.py b/rrompy/solver/norm_utilities.py index 96b151f..8e579fb 100644 --- a/rrompy/solver/norm_utilities.py +++ b/rrompy/solver/norm_utilities.py @@ -1,89 +1,94 @@ # Copyright (C) 2018 by the RROMPy authors # # This file is part of RROMPy. # # RROMPy is free software: you can redistribute it and/or modify # it under the terms of the GNU Lesser General Public License as published by # the Free Software Foundation, either version 3 of the License, or # (at your option) any later version. # # RROMPy is distributed in the hope that it will be useful, # but WITHOUT ANY WARRANTY; without even the implied warranty of # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the # GNU Lesser General Public License for more details. # # You should have received a copy of the GNU Lesser General Public License # along with RROMPy. If not, see . # from abc import abstractmethod import numpy as np from copy import deepcopy as copy from rrompy.utilities.base.types import Np1D, Np2D, DictAny from rrompy.solver.linear_solver import setupSolver from rrompy.utilities.exception_manager import RROMPyException __all__ = ['Np2DLike', 'Np2DLikeEye', 'Np2DLikeInv', 'Np2DLikeInvLowRank', 'normEngine'] @abstractmethod class Np2DLike: def dot(self, u:Np2D) -> Np2D: pass class Np2DLikeEye(Np2DLike): @property def T(self): return self - def conj(self): return self - def dot(self, u:Np2D) -> Np2D: return u class Np2DLikeInv(Np2DLike): def __init__(self, K:Np2D, M:Np2D, solverType:str, solverArgs:DictAny): self.K, self.M, self.MH = K, M, M.T.conj() try: self.solver, self.solverArgs = setupSolver(solverType, solverArgs) except: self.solver, self.solverArgs = solverType, solverArgs def dot(self, u:Np2D) -> Np2D: return self.MH.dot(self.solver(self.K, self.M.dot(u), self.solverArgs)).reshape(u.shape) + @property + def shape(self): + return (self.MH.shape[0], self.M.shape[1]) class Np2DLikeInvLowRank(Np2DLike): def __init__(self, K:Np2D, M:Np2D, solverType:str, solverArgs:DictAny, rank:int, oversampling : int = 10, seed : int = 420): - if rank > M.shape[1]: + sizeO = K.shape[1] if hasattr(K, "shape") else M.shape[1] + if rank > sizeO: raise RROMPyException(("Cannot select compressed rank larger than " "original size.")) if oversampling < 0: raise RROMPyException("Oversampling parameter must be positive.") HF = Np2DLikeInv(K, M, solverType, solverArgs) np.random.seed(seed) - xs = np.random.randn(M.shape[1], rank + oversampling) + xs = np.random.randn(sizeO, rank + oversampling) samples = HF.dot(xs) Q, _ = np.linalg.qr(samples, mode = "reduced") R = HF.dot(Q).T.conj() # assuming HF (i.e. K) hermitian... U, s, Vh = np.linalg.svd(R) self.L = Q.dot(U[:, : rank]) * s[: rank] self.R = Vh[: rank, :] def dot(self, u:Np2D) -> Np2D: return self.L.dot(self.R.dot(u)).reshape(u.shape) + @property + def shape(self): + return (self.L.shape[0], self.R.shape[1]) class normEngine: def __init__(self, energyNormMatrix:Np2D): self.energyNormMatrix = copy(energyNormMatrix) def innerProduct(self, u:Np2D, v:Np2D, onlyDiag : bool = False) -> Np2D: if not isinstance(u, (np.ndarray,)): u = u.data if not isinstance(v, (np.ndarray,)): v = v.data if onlyDiag: return np.sum(self.energyNormMatrix.dot(u) * v.conj(), axis = 0) return v.T.conj().dot(self.energyNormMatrix.dot(u)) def norm(self, u:Np2D) -> Np1D: return np.power(np.abs(self.innerProduct(u, u, onlyDiag = True)), .5) diff --git a/rrompy/utilities/poly_fitting/heaviside/heaviside_to_from_affine.py b/rrompy/utilities/poly_fitting/heaviside/heaviside_to_from_affine.py index 11cbc3b..2057f10 100644 --- a/rrompy/utilities/poly_fitting/heaviside/heaviside_to_from_affine.py +++ b/rrompy/utilities/poly_fitting/heaviside/heaviside_to_from_affine.py @@ -1,94 +1,94 @@ # Copyright (C) 2018 by the RROMPy authors # # This file is part of RROMPy. # # RROMPy is free software: you can redistribute it and/or modify # it under the terms of the GNU Lesser General Public License as published by # the Free Software Foundation, either version 3 of the License, or # (at your option) any later version. # # RROMPy is distributed in the hope that it will be useful, # but WITHOUT ANY WARRANTY; without even the implied warranty of # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the # GNU Lesser General Public License for more details. # # You should have received a copy of the GNU Lesser General Public License # along with RROMPy. If not, see . # import numpy as np from scipy.special import binom import scipy.sparse as sp from rrompy.utilities.base.types import (Np1D, Np2D, List, ListAny, Tuple, paramVal) from rrompy.utilities.numerical import eigNonlinearDense from rrompy.utilities.exception_manager import RROMPyException from rrompy.parameter import checkParameter __all__ = ['heaviside2affine', 'affine2heaviside'] def heaviside2affine(c:Np2D, poles:Np1D, mu : paramVal = [], basis : str = "MONOMIAL_HEAVISIDE", sparse : bool = False) \ -> Tuple[Np2D, List[Np2D], List[Np1D]]: mu = checkParameter(mu, 1)(0, 0) n, d = len(poles), len(c) - len(poles) basisN = basis.split("_")[0] if basisN not in ["MONOMIAL", "CHEBYSHEV", "LEGENDRE"]: raise RROMPyException("Polynomial basis not recognized.") if sparse: A0 = sp.spdiags([np.concatenate((- mu - poles, np.ones(d)))], [0], n + d, n + d) A1 = sp.spdiags([np.concatenate((np.ones(n), np.zeros(d)))], [0], n + d, n + d) else: A0 = np.diag(np.concatenate((mu - poles, np.ones(d)))) A1 = np.diag(np.concatenate((np.ones(n), np.zeros(d)))) As = [A0, A1] bs = np.zeros((d, n + d), dtype = poles.dtype) bs[0, :] = 1. if d > 0: bs[0, n + 1 :] = 0. if d > 1: bs[1, n + 1] = 1. for j in range(2, d): if basisN == "MONOMIAL": bs[j, n + j] = 1. else: alpha = - 1. if basisN == "CHEBYSHEV" else 1. / j - 1. bs[:, n + j] = alpha * bs[:, n + j - 2] bs[1 :, n + j] += (1. - alpha) * bs[: -1, n + j - 1] bs = list(bs) return c.reshape(c.shape[0], -1).T, As, bs def affine2heaviside(As:ListAny, bs:ListAny, jSupp : int = 1) -> Tuple[Np2D, Np1D, str]: if jSupp != 1 and not (isinstance(jSupp, (int,)) and jSupp.upper() == "COMPANION"): raise RROMPyException(("Affine to heaviside conversion does not allow " "nonlinear eigenproblem support outside first " "block row.")) N = len(As) M = len(bs) n = As[0].shape[0] if N == 1: poles = np.empty(0, dtype = np.complex) Q = np.eye(n) else: basis = "MONOMIAL_HEAVISIDE" poles, P, Q = eigNonlinearDense(As, jSupp = jSupp, return_inverse = True) P = P[- n :, :] Q = Q[:, : n] bEffs = np.array([Q.dot(np.linalg.solve(As[-1], b)) for b in bs]) if N == 1: c = bEffs else: - c = np.zeros((len(poles) + M - 1, As[0].shape[1]), dtype = As[0].dtype) + c = np.zeros((len(poles) + M - 1, As[0].shape[1]), dtype = np.complex) for l, pl in enumerate(poles): for i in range(M): c[l, :] = pl ** i * bEffs[i, l] * P[:, l] for l in range(M - 1): for i in range(l + 1, M): c[len(poles) + l, :] = P.dot(poles ** (i- 1 - l) * bEffs[i, :]) return c, poles, basis diff --git a/rrompy/utilities/poly_fitting/heaviside/vander.py b/rrompy/utilities/poly_fitting/heaviside/vander.py index 1480a8e..ce816d3 100644 --- a/rrompy/utilities/poly_fitting/heaviside/vander.py +++ b/rrompy/utilities/poly_fitting/heaviside/vander.py @@ -1,89 +1,89 @@ # Copyright (C) 2018 by the RROMPy authors # # This file is part of RROMPy. # # RROMPy is free software: you can redistribute it and/or modify # it under the terms of the GNU Lesser General Public License as published by # the Free Software Foundation, either version 3 of the License, or # (at your option) any later version. # # RROMPy is distributed in the hope that it will be useful, # but WITHOUT ANY WARRANTY; without even the implied warranty of # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the # GNU Lesser General Public License for more details. # # You should have received a copy of the GNU Lesser General Public License # along with RROMPy. If not, see . # import numpy as np from rrompy.utilities.poly_fitting.polynomial.vander import (polyvander as pvP, polyvanderTotal as pvTP) from rrompy.utilities.base.types import Np1D, Np2D, Tuple, List, paramList from rrompy.parameter import checkParameterList from rrompy.utilities.exception_manager import RROMPyException __all__ = ['heavisidevander', 'polyvander', 'polyvanderTotal'] def heavisidevander(x:paramList, poles:Np1D, reorder : List[int] = None) -> Np2D: """Compute Heaviside-Vandermonde matrix.""" x = checkParameterList(x, 1)[0] x_un, idx_un = x.unique(return_inverse = True) nx = len(x) if len(x_un) < nx: raise RROMPyException("Sample points must be distinct.") del x_un x = x.data.flatten() if reorder is not None: x = x[reorder] poles = poles.flatten() - Van = np.empty((len(x), len(poles)), dtype = poles.dtype) + Van = np.empty((len(x), len(poles)), dtype = np.complex) for j in range(len(x)): Van[j, :] = 1. / (x[j] - poles) return Van def polyvander(x:paramList, poles:Np1D, degs : List[int] = None, basis : str = "MONOMIAL_HEAVISIDE", derIdxs : List[List[List[int]]] = None, reorder : List[int] = None, scl : Np1D = None) -> Np2D: """ Compute full Hermite-Vandermonde matrix with specified derivative directions. """ if derIdxs is not None and np.sum(np.sum(derIdxs)) > 0: raise RROMPyException(("Cannot take derivatives of heaviside " "function.")) basisp = basis.split("_")[0] VanH = heavisidevander(x, poles, reorder = reorder) if degs is None or np.sum(degs) < 0: VanP = np.empty((len(x), 0)) else: VanP = pvP(x, degs, basisp, derIdxs = derIdxs, reorder = reorder, scl = scl) return np.block([[VanH, VanP]]) def polyvanderTotal(x:paramList, poles:Np1D, deg : int = None, basis : str = "MONOMIAL_HEAVISIDE", derIdxs : List[List[List[int]]] = None, reorder : List[int] = None, scl : Np1D = None)\ -> Tuple[Np2D, List[List[int]], List[int]]: """ Compute full total degree Hermite-Vandermonde matrix with specified derivative directions. """ if derIdxs is not None and np.sum(np.sum(derIdxs)) > 0: raise RROMPyException(("Cannot take derivatives of radial basis " "function.")) basisp = basis.split("_")[0] VanR = heavisidevander(x, poles, reorder = reorder) if deg is None or deg < 0: VanP = np.empty((len(x), 0)) derIdxs, ordIdxs = np.zeros(0, dtype = int), np.zeros(0, dtype = int) else: VanP, derIdxs, ordIdxs = pvTP(x, deg, basisp, derIdxs = derIdxs, reorder = reorder, scl = scl) ordIdxsEff = np.concatenate((np.arange(len(VanR)), ordIdxs + len(VanR))) return (np.block([[VanR, VanP], [VanP.T.conj(), np.zeros(tuple([VanP.shape[1]] * 2))]]), derIdxs, ordIdxsEff) diff --git a/rrompy/utilities/poly_fitting/moving_least_squares/__init__.py b/rrompy/utilities/poly_fitting/moving_least_squares/__init__.py index 8b8fd3f..d41094e 100644 --- a/rrompy/utilities/poly_fitting/moving_least_squares/__init__.py +++ b/rrompy/utilities/poly_fitting/moving_least_squares/__init__.py @@ -1,34 +1,40 @@ # Copyright (C) 2018 by the RROMPy authors # # This file is part of RROMPy. # # RROMPy is free software: you can redistribute it and/or modify # it under the terms of the GNU Lesser General Public License as published by # the Free Software Foundation, either version 3 of the License, or # (at your option) any later version. # # RROMPy is distributed in the hope that it will be useful, # but WITHOUT ANY WARRANTY; without even the implied warranty of # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the # GNU Lesser General Public License for more details. # # You should have received a copy of the GNU Lesser General Public License # along with RROMPy. If not, see . # +from .kernel import (radialGaussian, thinPlateSpline, multiQuadric, + nearestNeighbor) from .base import mlsbases, polybases, polyfitname, polydomcoeff from .vander import mlsweights, polyvander, polyvanderTotal from .moving_least_squares_interpolator import MovingLeastSquaresInterpolator __all__ = [ + 'radialGaussian', + 'thinPlateSpline', + 'multiQuadric', + 'nearestNeighbor', 'mlsbases', 'polybases', 'polyfitname', 'polydomcoeff', 'mlsweights', 'polyvander', 'polyvanderTotal', 'MovingLeastSquaresInterpolator' ] diff --git a/rrompy/reduction_methods/standard/__init__.py b/rrompy/utilities/poly_fitting/moving_least_squares/kernel.py similarity index 71% copy from rrompy/reduction_methods/standard/__init__.py copy to rrompy/utilities/poly_fitting/moving_least_squares/kernel.py index 0ddf98c..626f30b 100644 --- a/rrompy/reduction_methods/standard/__init__.py +++ b/rrompy/utilities/poly_fitting/moving_least_squares/kernel.py @@ -1,29 +1,24 @@ # Copyright (C) 2018 by the RROMPy authors # # This file is part of RROMPy. # # RROMPy is free software: you can redistribute it and/or modify # it under the terms of the GNU Lesser General Public License as published by # the Free Software Foundation, either version 3 of the License, or # (at your option) any later version. # # RROMPy is distributed in the hope that it will be useful, # but WITHOUT ANY WARRANTY; without even the implied warranty of # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the # GNU Lesser General Public License for more details. # # You should have received a copy of the GNU Lesser General Public License # along with RROMPy. If not, see . # -from .generic_standard_approximant import GenericStandardApproximant -from .rational_interpolant import RationalInterpolant -from .reduced_basis import ReducedBasis - -__all__ = [ - 'GenericStandardApproximant', - 'RationalInterpolant', - 'ReducedBasis' - ] +from rrompy.utilities.poly_fitting.radial_basis.kernel import (radialGaussian, + thinPlateSpline, multiQuadric, nearestNeighbor) +__all__ = ['radialGaussian', 'thinPlateSpline', 'multiQuadric', + 'nearestNeighbor'] diff --git a/rrompy/utilities/poly_fitting/moving_least_squares/moving_least_squares_interpolator.py b/rrompy/utilities/poly_fitting/moving_least_squares/moving_least_squares_interpolator.py index 2b27805..5753d91 100644 --- a/rrompy/utilities/poly_fitting/moving_least_squares/moving_least_squares_interpolator.py +++ b/rrompy/utilities/poly_fitting/moving_least_squares/moving_least_squares_interpolator.py @@ -1,142 +1,145 @@ # Copyright (C) 2018 by the RROMPy authors # # This file is part of RROMPy. # # RROMPy is free software: you can redistribute it and/or modify # it under the terms of the GNU Lesser General Public License as published by # the Free Software Foundation, either version 3 of the License, or # (at your option) any later version. # # RROMPy is distributed in the hope that it will be useful, # but WITHOUT ANY WARRANTY; without even the implied warranty of # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the # GNU Lesser General Public License for more details. # # You should have received a copy of the GNU Lesser General Public License # along with RROMPy. If not, see . # import numpy as np from copy import deepcopy as copy from rrompy.utilities.base.types import (List, ListAny, DictAny, Np1D, Np2D, paramList) from rrompy.utilities.numerical import customPInv from .vander import mlsweights from rrompy.utilities.poly_fitting.interpolator import GenericInterpolator from rrompy.utilities.poly_fitting.polynomial.vander import (polyvander as pv, polyvanderTotal as pvT) from rrompy.utilities.exception_manager import RROMPyException, RROMPyAssert from rrompy.parameter import checkParameterList __all__ = ['MovingLeastSquaresInterpolator'] class MovingLeastSquaresInterpolator(GenericInterpolator): """HERE""" def __init__(self, other = None): if other is None: return self.support = other.support self.localProjector = other.localProjector self.localVanders = other.localVanders self.suppValues = other.suppValues self.directionalWeights = other.directionalWeights self.degree = other.degree self.npar = other.npar self.radialbasis = other.radialbasis self.polybasis = other.polybasis self.evalParams = other.evalParams self.totalDegree = other.totalDegree @property def shape(self): sh = self.suppValues.shape[1 :] if self.suppValues.ndim > 1 else 1 return sh @property def deg(self): return self.degree def __call__(self, mu:paramList, der : List[int] = None, scl : Np1D = None): if der is not None and np.sum(der) > 0: raise RROMPyException(("Cannot take derivatives of moving least " "squares function.")) mu = checkParameterList(mu, self.npar)[0] sh = self.shape if sh == 1: sh = tuple([]) - values = np.empty((len(mu),) + sh, dtype = self.suppValues.dtype) + values = np.empty((len(mu),) + sh, dtype = np.complex) for i, m in enumerate(mu): weights = mlsweights(m, self.support, self.radialbasis, - directionalWeights = self.directionalWeights) + directionalWeights = self.directionalWeights, + nNearestNeighbor = self.evalParams["nNearestNeighbor"]) weights /= np.linalg.norm(weights) vanderLS = np.sum(self.localVanders * weights, axis = 2) RHSLS = np.tensordot(self.localProjector * weights, self.suppValues, 1) if self.totalDegree: vanderEval, _, _ = pvT(m, self.deg[0], self.polybasis, **self.evalParams) else: vanderEval = pv(m, self.deg, self.polybasis, **self.evalParams) vanderEval = vanderEval.flatten() values[i] = vanderEval.dot(customPInv(vanderLS).dot(RHSLS)) return values def __copy__(self): return MovingLeastSquaresInterpolator(self) def __deepcopy__(self, memo): other = MovingLeastSquaresInterpolator() (other.support, other.localProjector, other.localVanders, - other.suppValues, other.directionalWeights, other.degree, - other.npar, other.radialbasis, other.polybasis, other.evalParams, + other.suppValues, other.directionalWeights, other.degree, other.npar, + other.radialbasis, other.polybasis, other.evalParams, other.totalDegree) = copy( - self.support, self.localProjector, self.localVanders, + (self.support, self.localProjector, self.localVanders, self.suppValues, self.directionalWeights, self.degree, self.npar, self.radialbasis, self.polybasis, - self.evalParams, self.totalDegree, memo) + self.evalParams, self.totalDegree), memo) return other def postmultiplyTensorize(self, A:Np2D): RROMPyAssert(A.shape[0], self.shape[-1], "Shape of output") self.suppValues = np.tensordot(self.suppValues, A, axes = (-1, 0)) def pad(self, nleft : List[int] = None, nright : List[int] = None): if nleft is None: nleft = [0] * len(self.shape) if nright is None: nright = [0] * len(self.shape) if not hasattr(nleft, "__len__"): nleft = [nleft] if not hasattr(nright, "__len__"): nright = [nright] RROMPyAssert(len(self.shape), len(nleft), "Shape of output") RROMPyAssert(len(self.shape), len(nright), "Shape of output") padwidth = [(0, 0)] + [(l, r) for l, r in zip(nleft, nright)] self.suppValues = np.pad(self.suppValues, padwidth, "constant", constant_values = (0., 0.)) def setupByInterpolation(self, support:paramList, values:ListAny, deg:int, polybasis : str = "MONOMIAL_GAUSSIAN", directionalWeights : Np1D = None, totalDegree : bool = True, vanderCoeffs : DictAny = {}): support = checkParameterList(support)[0] self.support = copy(support) if "reorder" in vanderCoeffs.keys(): self.support = self.support[vanderCoeffs["reorder"]] + if "nNearestNeighbor" not in vanderCoeffs.keys(): + vanderCoeffs["nNearestNeighbor"] = -1 self.npar = support.shape[1] if directionalWeights is None: directionalWeights = np.ones(self.npar) self.directionalWeights = directionalWeights self.polybasis, self.radialbasis, _ = polybasis.split("_") self.totalDegree = totalDegree self.evalParams = vanderCoeffs if totalDegree: vander, _, _ = pvT(support, deg, self.polybasis, **vanderCoeffs) if not hasattr(deg, "__len__"): deg = [deg] * self.npar else: if not hasattr(deg, "__len__"): deg = [deg] * self.npar vander = pv(support, deg, self.polybasis, **vanderCoeffs) self.degree = deg self.localProjector = vander.T.conj() self.localVanders = np.array([np.outer(van, van.conj()) \ for van in vander]) self.localVanders = np.swapaxes(self.localVanders, 0, 2) self.suppValues = np.array(values) diff --git a/rrompy/utilities/poly_fitting/moving_least_squares/vander.py b/rrompy/utilities/poly_fitting/moving_least_squares/vander.py index 66f6c3d..04094e1 100644 --- a/rrompy/utilities/poly_fitting/moving_least_squares/vander.py +++ b/rrompy/utilities/poly_fitting/moving_least_squares/vander.py @@ -1,83 +1,97 @@ # Copyright (C) 2018 by the RROMPy authors # # This file is part of RROMPy. # # RROMPy is free software: you can redistribute it and/or modify # it under the terms of the GNU Lesser General Public License as published by # the Free Software Foundation, either version 3 of the License, or # (at your option) any later version. # # RROMPy is distributed in the hope that it will be useful, # but WITHOUT ANY WARRANTY; without even the implied warranty of # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the # GNU Lesser General Public License for more details. # # You should have received a copy of the GNU Lesser General Public License # along with RROMPy. If not, see . # import numpy as np -from rrompy.utilities.poly_fitting.radial_basis.kernel import (radialGaussian, - thinPlateSpline, multiQuadric) +from .kernel import (radialGaussian, thinPlateSpline, multiQuadric, + nearestNeighbor) from rrompy.utilities.poly_fitting.polynomial.vander import (polyvander as pvP, polyvanderTotal as pvTP) from rrompy.utilities.base.types import (Np1D, Np2D, Tuple, List, paramVal, paramList) from rrompy.parameter import checkParameter, checkParameterList from rrompy.utilities.exception_manager import RROMPyException, RROMPyAssert __all__ = ['mlsweights', 'polyvander', 'polyvanderTotal'] def mlsweights(x:paramVal, xSupp:paramList, basis:str, - reorder : List[int] = None, - directionalWeights : Np1D = None) -> Np2D: + reorder : List[int] = None, directionalWeights : Np1D = None, + nNearestNeighbor : int = -1) -> Np2D: """Compute moving least squares weight vector.""" x = checkParameter(x) xSupp = checkParameterList(xSupp)[0] x = x.data xSupp = xSupp.data if directionalWeights is None: directionalWeights = np.ones(x.shape[1]) + elif not hasattr(directionalWeights, "__len__"): + directionalWeights = directionalWeights * np.ones(x.shape[1]) RROMPyAssert(len(directionalWeights), x.shape[1], "Number of directional weights") try: radialkernel = {"GAUSSIAN" : radialGaussian, "THINPLATE" : thinPlateSpline, - "MULTIQUADRIC" : multiQuadric - }[basis.upper()] + "MULTIQUADRIC" : multiQuadric, + "NEARESTNEIGHBOR" : nearestNeighbor}[basis.upper()] except: raise RROMPyException("Radial basis not recognized.") if reorder is not None: xSupp = xSupp[reorder] - r2 = np.sum(np.abs((xSupp - x) * directionalWeights) ** 2., axis = 1) - return radialkernel(r2) + muDiff = (xSupp.data - x) * directionalWeights + r2 = np.sum(np.abs(muDiff) ** 2., axis = 1).reshape(1, -1) + if basis.upper() == "NEARESTNEIGHBOR": + if nNearestNeighbor > 0 and nNearestNeighbor < len(xSupp): + cutoffValue = np.partition(r2, nNearestNeighbor - 1)[0, + nNearestNeighbor - 1] + r2 /= cutoffValue + else: + r2[0, :] = 1. * (nNearestNeighbor == 0) + return radialkernel(r2)[0] def polyvander(x:paramVal, xSupp:paramList, degs:List[int], basis:str, derIdxs : List[List[List[int]]] = None, reorder : List[int] = None, directionalWeights : Np1D = None, - scl : Np1D = None) -> Tuple[Np2D, Np2D]: + scl : Np1D = None, + nNearestNeighbor : int = -1) -> Tuple[Np2D, Np2D]: """ Compute full Hermite-Vandermonde matrix with specified derivative directions. """ basisp, basisr, _ = basis.split("_") - Weights = mlsweights(x, xSupp, basisr, reorder, directionalWeights) + Weights = mlsweights(x, xSupp, basisr, reorder, directionalWeights, + nNearestNeighbor) VanP = pvP(xSupp, degs, basisp, derIdxs = derIdxs, reorder = reorder, scl = scl) RHP = VanP.T.conj() * Weights return RHP.dot(VanP), RHP def polyvanderTotal(x:paramList, xSupp:paramList, deg:int, basis:str, derIdxs : List[List[List[int]]] = None, reorder : List[int] = None, - directionalWeights : Np1D = None, scl : Np1D = None)\ + directionalWeights : Np1D = None, scl : Np1D = None, + nNearestNeighbor : int = -1)\ -> Tuple[Np2D, Np2D, List[List[int]], List[int]]: """ Compute full total degree Hermite-Vandermonde matrix with specified derivative directions. """ basisp, basisr, _ = basis.split("_") - Weights = mlsweights(x, xSupp, basisr, reorder, directionalWeights) + Weights = mlsweights(x, xSupp, basisr, reorder, directionalWeights, + nNearestNeighbor) VanP, derIdxs, ordIdxs = pvTP(x, deg, basisp, derIdxs = derIdxs, reorder = reorder, scl = scl) RHP = VanP.T.conj() * Weights return RHP.dot(VanP), RHP, derIdxs, ordIdxs diff --git a/rrompy/utilities/poly_fitting/radial_basis/__init__.py b/rrompy/utilities/poly_fitting/radial_basis/__init__.py index c94c22a..c660671 100644 --- a/rrompy/utilities/poly_fitting/radial_basis/__init__.py +++ b/rrompy/utilities/poly_fitting/radial_basis/__init__.py @@ -1,40 +1,42 @@ # Copyright (C) 2018 by the RROMPy authors # # This file is part of RROMPy. # # RROMPy is free software: you can redistribute it and/or modify # it under the terms of the GNU Lesser General Public License as published by # the Free Software Foundation, either version 3 of the License, or # (at your option) any later version. # # RROMPy is distributed in the hope that it will be useful, # but WITHOUT ANY WARRANTY; without even the implied warranty of # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the # GNU Lesser General Public License for more details. # # You should have received a copy of the GNU Lesser General Public License # along with RROMPy. If not, see . # -from .kernel import radialGaussian, thinPlateSpline, multiQuadric +from .kernel import (radialGaussian, thinPlateSpline, multiQuadric, + nearestNeighbor) from .base import rbbases, polybases, polyfitname, polydomcoeff from .val import polyval from .vander import rbvander, polyvander, polyvanderTotal from .radial_basis_interpolator import RadialBasisInterpolator __all__ = [ 'radialGaussian', 'thinPlateSpline', 'multiQuadric', + 'nearestNeighbor', 'rbbases', 'polybases', 'polyfitname', 'polydomcoeff', 'polyval', 'rbvander', 'polyvander', 'polyvanderTotal', 'RadialBasisInterpolator' ] diff --git a/rrompy/utilities/poly_fitting/radial_basis/base.py b/rrompy/utilities/poly_fitting/radial_basis/base.py index 2638b7e..29c4ba6 100644 --- a/rrompy/utilities/poly_fitting/radial_basis/base.py +++ b/rrompy/utilities/poly_fitting/radial_basis/base.py @@ -1,43 +1,44 @@ # Copyright (C) 2018 by the RROMPy authors # # This file is part of RROMPy. # # RROMPy is free software: you can redistribute it and/or modify # it under the terms of the GNU Lesser General Public License as published by # the Free Software Foundation, either version 3 of the License, or # (at your option) any later version. # # RROMPy is distributed in the hope that it will be useful, # but WITHOUT ANY WARRANTY; without even the implied warranty of # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the # GNU Lesser General Public License for more details. # # You should have received a copy of the GNU Lesser General Public License # along with RROMPy. If not, see . # from itertools import product from rrompy.utilities.exception_manager import RROMPyException from rrompy.utilities.poly_fitting.polynomial.base import (polybases as pbP, polyfitname as pfnP, polydomcoeff as polydomcoeffB) __all__ = ['rbbases', 'polybases', 'polyfitname', 'polydomcoeff'] -rbbases = ["GAUSSIAN", "THINPLATE", "MULTIQUADRIC"] +rbbases = ["GAUSSIAN", "THINPLATE", "MULTIQUADRIC", "NEARESTNEIGHBOR"] polybases = [x + "_" + y for x, y in product(pbP, rbbases)] def polyfitname(basis:str) -> str: fitrnames = {"GAUSSIAN" : "gaussian", "THINPLATE" : "thinplate", - "MULTIQUADRIC" : "multiquadric"} + "MULTIQUADRIC" : "multiquadric", + "NEARESTNEIGHBOR" : "nearestneighbor"} basisp, basisr = basis.split("_") try: return pfnP(basisp) + "_" + fitrnames[basisr] except: raise RROMPyException("Polynomial-radial basis combination not " "recognized.") def polydomcoeff(n:int, basis:str) -> float: return polydomcoeffB(n, basis.split("_")[0]) diff --git a/rrompy/utilities/poly_fitting/radial_basis/kernel.py b/rrompy/utilities/poly_fitting/radial_basis/kernel.py index d0e1db2..d51bddb 100644 --- a/rrompy/utilities/poly_fitting/radial_basis/kernel.py +++ b/rrompy/utilities/poly_fitting/radial_basis/kernel.py @@ -1,35 +1,40 @@ # Copyright (C) 2018 by the RROMPy authors # # This file is part of RROMPy. # # RROMPy is free software: you can redistribute it and/or modify # it under the terms of the GNU Lesser General Public License as published by # the Free Software Foundation, either version 3 of the License, or # (at your option) any later version. # # RROMPy is distributed in the hope that it will be useful, # but WITHOUT ANY WARRANTY; without even the implied warranty of # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the # GNU Lesser General Public License for more details. # # You should have received a copy of the GNU Lesser General Public License # along with RROMPy. If not, see . # import numpy as np from rrompy.utilities.base.types import Np1D from rrompy.utilities.exception_manager import RROMPyAssert -__all__ = ['radialGaussian', 'thinPlateSpline', 'multiQuadric'] +__all__ = ['radialGaussian', 'thinPlateSpline', 'multiQuadric', + 'nearestNeighbor'] def radialGaussian(r2:Np1D, der : int = 0) -> Np1D: RROMPyAssert(der, 0, "Radial basis derivative") return np.exp(- .5 * r2) def thinPlateSpline(r2:Np1D, der : int = 0) -> Np1D: RROMPyAssert(der, 0, "Radial basis derivative") return .5 * r2 * np.log(np.finfo(float).eps + r2) def multiQuadric(r2:Np1D, der : int = 0) -> Np1D: RROMPyAssert(der, 0, "Radial basis derivative") return np.power(r2 + 1., -.5) + +def nearestNeighbor(r2:Np1D, der : int = 0) -> Np1D: + RROMPyAssert(der, 0, "Radial basis derivative") + return 1. * (r2 <= 1.) diff --git a/rrompy/utilities/poly_fitting/radial_basis/radial_basis_interpolator.py b/rrompy/utilities/poly_fitting/radial_basis/radial_basis_interpolator.py index 325de85..1a1b74d 100644 --- a/rrompy/utilities/poly_fitting/radial_basis/radial_basis_interpolator.py +++ b/rrompy/utilities/poly_fitting/radial_basis/radial_basis_interpolator.py @@ -1,139 +1,147 @@ # Copyright (C) 2018 by the RROMPy authors # # This file is part of RROMPy. # # RROMPy is free software: you can redistribute it and/or modify # it under the terms of the GNU Lesser General Public License as published by # the Free Software Foundation, either version 3 of the License, or # (at your option) any later version. # # RROMPy is distributed in the hope that it will be useful, # but WITHOUT ANY WARRANTY; without even the implied warranty of # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the # GNU Lesser General Public License for more details. # # You should have received a copy of the GNU Lesser General Public License # along with RROMPy. If not, see . # import numpy as np from copy import deepcopy as copy from rrompy.utilities.base.types import (List, ListAny, DictAny, Np1D, Np2D, paramList) from rrompy.utilities.poly_fitting.interpolator import GenericInterpolator from rrompy.utilities.poly_fitting.custom_fit import customFit from .base import polyfitname from .val import polyval from .vander import polyvander as pv, polyvanderTotal as pvT from rrompy.utilities.numerical import degreeTotalToFull from rrompy.utilities.exception_manager import RROMPyException, RROMPyAssert from rrompy.parameter import checkParameterList __all__ = ['RadialBasisInterpolator'] class RadialBasisInterpolator(GenericInterpolator): """HERE""" def __init__(self, other = None): if other is None: return self.support = other.support self.coeffsGlobal = other.coeffsGlobal self.coeffsLocal = other.coeffsLocal self.directionalWeights = other.directionalWeights self.npar = other.npar self.polybasis = other.polybasis + self.nNearestNeighbor = other.nNearestNeighbor @property def shape(self): sh = self.coeffsLocal.shape[1 :] if self.coeffsLocal.ndim > 1 else 1 return sh @property def deg(self): return [x - 1 for x in self.coeffsGlobal.shape[: self.npar]] def __call__(self, mu:paramList, der : List[int] = None, scl : Np1D = None): if der is not None and np.sum(der) > 0: raise RROMPyException(("Cannot take derivatives of radial basis " "function.")) return polyval(mu, self.coeffsGlobal, self.coeffsLocal, self.support, - self.directionalWeights, self.polybasis) + self.directionalWeights, self.polybasis, + self.nNearestNeighbor) def __copy__(self): return RadialBasisInterpolator(self) def __deepcopy__(self, memo): other = RadialBasisInterpolator() (other.support, other.coeffsGlobal, other.coeffsLocal, - other.directionalWeights, other.npar, other.polybasis) = copy( - (self.support, self.coeffsGlobal, self.coeffsLocal, - self.directionalWeights, self.npar, self.polybasis), memo) + other.directionalWeights, other.npar, other.polybasis, + other.nNearestNeighbor) = copy( + (self.support, self.coeffsGlobal, self.coeffsLocal, + self.directionalWeights, self.npar, self.polybasis, + self.nNearestNeighbor), memo) return other def postmultiplyTensorize(self, A:Np2D): RROMPyAssert(A.shape[0], self.shape[-1], "Shape of output") self.coeffsLocal = np.tensordot(self.coeffsLocal, A, axes = (-1, 0)) self.coeffsGlobal = np.tensordot(self.coeffsGlobal, A, axes = (-1, 0)) def pad(self, nleft : List[int] = None, nright : List[int] = None): if nleft is None: nleft = [0] * len(self.shape) if nright is None: nright = [0] * len(self.shape) if not hasattr(nleft, "__len__"): nleft = [nleft] if not hasattr(nright, "__len__"): nright = [nright] RROMPyAssert(len(self.shape), len(nleft), "Shape of output") RROMPyAssert(len(self.shape), len(nright), "Shape of output") padwidth = [(0, 0)] + [(l, r) for l, r in zip(nleft, nright)] self.coeffsLocal = np.pad(self.coeffsLocal, padwidth, "constant", constant_values = (0., 0.)) padwidth = [(0, 0)] * (self.npar - 1) + padwidth self.coeffsGlobal = np.pad(self.coeffsGlobal, padwidth, "constant", constant_values = (0., 0.)) def setupByInterpolation(self, support:paramList, values:ListAny, deg:int, polybasis : str = "MONOMIAL_GAUSSIAN", directionalWeights : Np1D = None, verbose : bool = True, totalDegree : bool = True, vanderCoeffs : DictAny = {}, fitCoeffs : DictAny = {}): support = checkParameterList(support)[0] self.support = copy(support) if "reorder" in vanderCoeffs.keys(): self.support = self.support[vanderCoeffs["reorder"]] + if "nNearestNeighbor" in vanderCoeffs.keys(): + self.nNearestNeighbor = vanderCoeffs["nNearestNeighbor"] + else: + self.nNearestNeighbor = -1 self.npar = support.shape[1] if directionalWeights is None: directionalWeights = np.ones(self.npar) self.directionalWeights = directionalWeights self.polybasis = polybasis if totalDegree: vander, _, reorder = pvT(support, deg, basis = polybasis, directionalWeights = self.directionalWeights, **vanderCoeffs) vander = vander[reorder] vander = vander[:, reorder] else: if not hasattr(deg, "__len__"): deg = [deg] * self.npar vander = pv(support, deg, basis = polybasis, directionalWeights = self.directionalWeights, **vanderCoeffs) outDim = values.shape[1:] values = values.reshape(values.shape[0], -1) values = np.pad(values, ((0, len(vander) - len(values)), (0, 0)), "constant") fitOut = customFit(vander, values, full = True, **fitCoeffs) P = fitOut[0][len(support) :] if verbose: msg = ("Fitting {}+{} samples with degree {} through {}... " "Conditioning of LS system: {:.4e}.").format( len(support), len(vander) - len(support), deg, polyfitname(self.polybasis), fitOut[1][2][0] / fitOut[1][2][-1]) else: msg = None self.coeffsLocal = fitOut[0][: len(support)] if totalDegree: self.coeffsGlobal = degreeTotalToFull(tuple([deg + 1] * self.npar) + outDim, self.npar, P) else: self.coeffsGlobal = P.reshape(tuple([d + 1 for d in deg]) + outDim) return fitOut[1][1] == vander.shape[1], msg diff --git a/rrompy/utilities/poly_fitting/radial_basis/val.py b/rrompy/utilities/poly_fitting/radial_basis/val.py index 46847a4..170e821 100644 --- a/rrompy/utilities/poly_fitting/radial_basis/val.py +++ b/rrompy/utilities/poly_fitting/radial_basis/val.py @@ -1,52 +1,62 @@ # Copyright (C) 2018 by the RROMPy authors # # This file is part of RROMPy. # # RROMPy is free software: you can redistribute it and/or modify # it under the terms of the GNU Lesser General Public License as published by # the Free Software Foundation, either version 3 of the License, or # (at your option) any later version. # # RROMPy is distributed in the hope that it will be useful, # but WITHOUT ANY WARRANTY; without even the implied warranty of # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the # GNU Lesser General Public License for more details. # # You should have received a copy of the GNU Lesser General Public License # along with RROMPy. If not, see . # from copy import deepcopy as copy import numpy as np -from .kernel import radialGaussian, thinPlateSpline, multiQuadric +from .kernel import (radialGaussian, thinPlateSpline, multiQuadric, + nearestNeighbor) from rrompy.utilities.poly_fitting.polynomial.val import polyval as pvP from rrompy.utilities.base.types import Np1D, Np2D, paramList from rrompy.parameter import checkParameterList from rrompy.utilities.exception_manager import RROMPyException __all__ = ['polyval'] def polyval(x:paramList, cG:Np2D, cL:Np2D, supportPoints:paramList, - directionalWeights:Np1D, basis:str) -> Np2D: + directionalWeights:Np1D, basis:str, + nNearestNeighbor : int = -1) -> Np2D: x = checkParameterList(x)[0] basisp, basisr = basis.split("_") c = pvP(x, cG, basisp) try: radialvalbase = {"GAUSSIAN" : radialGaussian, "THINPLATE" : thinPlateSpline, - "MULTIQUADRIC" : multiQuadric - }[basisr.upper()] + "MULTIQUADRIC" : multiQuadric, + "NEARESTNEIGHBOR" : nearestNeighbor}[basisr.upper()] except: raise RROMPyException("Radial basis not recognized.") + isnearestneighbor = basisr.upper() == "NEARESTNEIGHBOR" csh = copy(c.shape) if len(csh) == 1: c = c.reshape(1, -1) for j in range(len(x)): muDiff = (supportPoints.data - x[j]) * directionalWeights r2j = np.sum(np.abs(muDiff) ** 2., axis = 1).reshape(1, -1) + if isnearestneighbor: + if nNearestNeighbor > 0 and nNearestNeighbor < len(supportPoints): + cutoffValue = np.partition(r2j, nNearestNeighbor - 1)[0, + nNearestNeighbor - 1] + r2j /= cutoffValue + else: + r2j[0, :] = 1. * (nNearestNeighbor == 0) val = radialvalbase(r2j).dot(cL) try: c[..., j] += val except: c[..., j] += val.flatten() if len(csh) == 1: c = c.flatten() return c diff --git a/rrompy/utilities/poly_fitting/radial_basis/vander.py b/rrompy/utilities/poly_fitting/radial_basis/vander.py index c5b9d62..9738162 100644 --- a/rrompy/utilities/poly_fitting/radial_basis/vander.py +++ b/rrompy/utilities/poly_fitting/radial_basis/vander.py @@ -1,102 +1,111 @@ # Copyright (C) 2018 by the RROMPy authors # # This file is part of RROMPy. # # RROMPy is free software: you can redistribute it and/or modify # it under the terms of the GNU Lesser General Public License as published by # the Free Software Foundation, either version 3 of the License, or # (at your option) any later version. # # RROMPy is distributed in the hope that it will be useful, # but WITHOUT ANY WARRANTY; without even the implied warranty of # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the # GNU Lesser General Public License for more details. # # You should have received a copy of the GNU Lesser General Public License # along with RROMPy. If not, see . # import numpy as np -from .kernel import radialGaussian, thinPlateSpline, multiQuadric +from .kernel import (radialGaussian, thinPlateSpline, multiQuadric, + nearestNeighbor) from rrompy.utilities.poly_fitting.polynomial.vander import (polyvander as pvP, polyvanderTotal as pvTP) from rrompy.utilities.base.types import Np1D, Np2D, Tuple, List, paramList from rrompy.parameter import checkParameterList from rrompy.utilities.exception_manager import RROMPyException, RROMPyAssert __all__ = ['rbvander', 'polyvander', 'polyvanderTotal'] def rbvander(x:paramList, basis:str, reorder : List[int] = None, - directionalWeights : Np1D = None) -> Np2D: + directionalWeights : Np1D = None, + nNearestNeighbor : int = -1) -> Np2D: """Compute radial-basis-Vandermonde matrix.""" x = checkParameterList(x)[0] x_un = x.unique() nx = len(x) if len(x_un) < nx: raise RROMPyException("Sample points must be distinct.") del x_un x = x.data if directionalWeights is None: directionalWeights = np.ones(x.shape[1]) + elif not hasattr(directionalWeights, "__len__"): + directionalWeights = directionalWeights * np.ones(x.shape[1]) RROMPyAssert(len(directionalWeights), x.shape[1], "Number of directional weights") try: radialkernel = {"GAUSSIAN" : radialGaussian, "THINPLATE" : thinPlateSpline, - "MULTIQUADRIC" : multiQuadric - }[basis.upper()] + "MULTIQUADRIC" : multiQuadric, + "NEARESTNEIGHBOR" : nearestNeighbor}[basis.upper()] except: raise RROMPyException("Radial basis not recognized.") - r2 = np.zeros((nx * (nx - 1) // 2 + 1), dtype = float) - idxInv = np.zeros(nx ** 2, dtype = int) - if reorder is not None: x = x[reorder] + isnearestneighbor = basis.upper() == "NEARESTNEIGHBOR" + Van = np.zeros((nx, nx)) for j in range(nx): - idx = j * (nx - 1) - j * (j + 1) // 2 - II = np.arange(j + 1, nx) - r2[idx + II] = np.sum(np.abs((x[II] - x[j]) - * directionalWeights) ** 2., axis = 1) - idxInv[j * nx + II] = idx + II - idxInv[II * nx + j] = idx + II - Van = radialkernel(r2)[idxInv].reshape((nx, nx)) + muDiff = (x.data - x[j]) * directionalWeights + r2j = np.sum(np.abs(muDiff) ** 2., axis = 1).reshape(1, -1) + if isnearestneighbor: + if nNearestNeighbor > 0 and nNearestNeighbor < len(x): + cutoffValue = np.partition(r2j, nNearestNeighbor - 1)[0, + nNearestNeighbor - 1] + r2j /= cutoffValue + else: + r2j[0, :] = 1. * (nNearestNeighbor == 0) + Van[j] = radialkernel(r2j) return Van def polyvander(x:paramList, degs:List[int], basis:str, derIdxs : List[List[List[int]]] = None, reorder : List[int] = None, directionalWeights : Np1D = None, - scl : Np1D = None) -> Np2D: + scl : Np1D = None, nNearestNeighbor : int = -1) -> Np2D: """ Compute full Hermite-Vandermonde matrix with specified derivative directions. """ if derIdxs is not None and np.sum(np.sum(derIdxs)) > 0: raise RROMPyException(("Cannot take derivatives of radial basis " "function.")) basisp, basisr = basis.split("_") VanR = rbvander(x, basisr, reorder = reorder, - directionalWeights = directionalWeights) + directionalWeights = directionalWeights, + nNearestNeighbor = nNearestNeighbor) VanP = pvP(x, degs, basisp, derIdxs = derIdxs, reorder = reorder, scl = scl) return np.block([[VanR, VanP], [VanP.T.conj(), np.zeros(tuple([VanP.shape[1]] * 2))]]) def polyvanderTotal(x:paramList, deg:int, basis:str, derIdxs : List[List[List[int]]] = None, reorder : List[int] = None, - directionalWeights : Np1D = None, scl : Np1D = None)\ + directionalWeights : Np1D = None, scl : Np1D = None, + nNearestNeighbor : int = -1)\ -> Tuple[Np2D, List[List[int]], List[int]]: """ Compute full total degree Hermite-Vandermonde matrix with specified derivative directions. """ if derIdxs is not None and np.sum(np.sum(derIdxs)) > 0: raise RROMPyException(("Cannot take derivatives of radial basis " "function.")) basisp, basisr = basis.split("_") VanR = rbvander(x, basisr, reorder = reorder, - directionalWeights = directionalWeights) + directionalWeights = directionalWeights, + nNearestNeighbor = nNearestNeighbor) VanP, derIdxs, ordIdxs = pvTP(x, deg, basisp, derIdxs = derIdxs, reorder = reorder, scl = scl) ordIdxsEff = np.concatenate((np.arange(len(VanR)), ordIdxs + len(VanR))) return (np.block([[VanR, VanP], [VanP.T.conj(), np.zeros(tuple([VanP.shape[1]] * 2))]]), derIdxs, ordIdxsEff) diff --git a/tests/utilities/parameter_sampling.py b/tests/utilities/parameter_sampling.py index d7d1766..cdc3009 100644 --- a/tests/utilities/parameter_sampling.py +++ b/tests/utilities/parameter_sampling.py @@ -1,59 +1,59 @@ # Copyright (C) 2018 by the RROMPy authors # # This file is part of RROMPy. # # RROMPy is free software: you can redistribute it and/or modify # it under the terms of the GNU Lesser General Public License as published by # the Free Software Foundation, either version 3 of the License, or # (at your option) any later version. # # RROMPy is distributed in the hope that it will be useful, # but WITHOUT ANY WARRANTY; without even the implied warranty of # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the # GNU Lesser General Public License for more details. # # You should have received a copy of the GNU Lesser General Public License # along with RROMPy. If not, see . # import numpy as np from rrompy.parameter.parameter_sampling import (ManualSampler, QuadratureSampler, RandomSampler, FFTSampler) from rrompy.parameter import checkParameter def test_manual(): sampler = ManualSampler(lims = [0., 3.], points = np.linspace(0, 3, 101), scalingExp = 2.) assert sampler.name() == "ManualSampler" x = sampler.generatePoints(10) assert np.allclose(x(0), np.linspace(0, 3, 101)[:10], rtol = 1e-5) def test_quadrature(): sampler = QuadratureSampler(lims = [0., 3.], kind = "CHEBYSHEV") x = sampler.generatePoints(9, reorder = False) assert np.isclose(x(0)[4], 1.5, rtol = 1e-5) def test_random(): sampler = RandomSampler(lims = [0., 3.], kind = "SOBOL", seed = 13432) x = sampler.generatePoints(100) assert np.isclose(x(0)[47], 0.55609130859375, rtol = 1e-5) def test_fft(): sampler = FFTSampler(lims = [-1., 1.]) x = sampler.generatePoints(100) assert np.allclose(np.power(x(0), 100), 1., rtol = 1e-5) def test_2D(): sampler = QuadratureSampler(lims = [(0., 0.), (3., 1.)], kind = "GAUSSLEGENDRE") x = sampler.generatePoints(81) assert sum(np.isclose(x(0), 1.5)) == 9 assert sum(np.isclose(x(1), .5)) == 9 def test_4D(): - sampler = RandomSampler(lims = [tuple([0.] * 4), tuple([1.] * 4)], + sampler = RandomSampler(lims = [(0.,) * 4, (1.,) * 4], kind = "UNIFORM", seed = 1234) x = sampler.generatePoints(10) assert x.shape[1] == 4 assert checkParameter([x[0]]) == checkParameter([(0.191519450378892, 0.622108771039832, 0.437727739007115, 0.785358583713769)])