diff --git a/rrompy/reduction_methods/greedy/rational_interpolant_greedy.py b/rrompy/reduction_methods/greedy/rational_interpolant_greedy.py index 4ff7253..c4ba5ba 100644 --- a/rrompy/reduction_methods/greedy/rational_interpolant_greedy.py +++ b/rrompy/reduction_methods/greedy/rational_interpolant_greedy.py @@ -1,353 +1,354 @@ # 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 from rrompy.utilities.poly_fitting.polynomial import (polybases, polydomcoeff, PolynomialInterpolator as PI) 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, Np2D, DictAny, HFEng, paramVal from rrompy.utilities.base import verbosityManager as vbMng from rrompy.utilities.exception_manager import RROMPyWarning from rrompy.utilities.exception_manager import RROMPyException, RROMPyAssert __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; - 'radialBasis': radial basis family for interpolant numerator; defaults to 0, i.e. no radial basis; - 'radialBasisWeights': radial basis weights for interpolant numerator; defaults to 0, i.e. identity; - 'greedyTol': uniform error tolerance for greedy algorithm; defaults to 1e-2; - '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; - 'nTestPoints': number of test points; defaults to 5e2; - 'trainSetGenerator': training sample points generator; defaults to sampler; - 'polybasis': type of basis for interpolation; allowed values include 'MONOMIAL', 'CHEBYSHEV' and 'LEGENDRE'; defaults to 'MONOMIAL'; - 'errorEstimatorKind': kind of error estimator; available values include 'EXACT', 'BASIC', and 'BARE'; defaults to 'EXACT'; - 'interpRcond': tolerance for interpolation; defaults to None; - 'robustTol': tolerance for robust rational denominator management; defaults to 0. Defaults to empty dict. homogeneized(optional): Whether to homogeneize Dirichlet BCs. Defaults to False. verbosity(optional): Verbosity level. Defaults to 10. Attributes: HFEngine: HF problem solver. mu0: Default parameter. mus: Array of snapshot parameters. homogeneized: Whether to homogeneize Dirichlet BCs. 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. - 'radialBasis': radial basis family for interpolant numerator; defaults to 0, i.e. no radial basis; - 'radialBasisWeights': radial basis weights for interpolant numerator; defaults to 0, i.e. identity; - 'greedyTol': uniform error 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. radialBasis: Radial basis family for interpolant numerator. radialBasisWeights: Radial basis weights for interpolant numerator. greedyTol: uniform error 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 = ["EXACT", "BASIC", "BARE"] def __init__(self, HFEngine:HFEng, mu0 : paramVal = None, approxParameters : DictAny = {}, homogeneized : bool = False, verbosity : int = 10, timestamp : bool = True): self._preInit() self._addParametersToList(["polybasis", "errorEstimatorKind", "interpRcond", "robustTol"], ["MONOMIAL", "EXACT", -1, 0], toBeExcluded = ["E"]) super().__init__(HFEngine = HFEngine, mu0 = mu0, approxParameters = approxParameters, homogeneized = homogeneized, verbosity = verbosity, timestamp = timestamp) self._postInit() @property def E(self): """Value of E.""" self._E = np.prod(self._S) - 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 prod(S) - 1.")) @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 'EXACT'.")) errorEstimatorKind = "EXACT" self._errorEstimatorKind = errorEstimatorKind self._approxParameters["errorEstimatorKind"] = self.errorEstimatorKind @property def nTestPoints(self): """Value of nTestPoints.""" return self._nTestPoints @nTestPoints.setter def nTestPoints(self, nTestPoints): if nTestPoints <= 0: RROMPyWarning("nTestPoints must be at least 1. Overriding to 1.") nTestPoints = 1 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() def _errorSamplingRatio(self, mus:Np1D, muRTest:Np1D, muRTrain:Np1D) -> Np1D: """Scalar ratio in explicit error estimator.""" self.setupApprox() testTile = np.tile(np.reshape(muRTest, (-1, 1)), [1, len(muRTrain)]) nodalVals = np.prod(testTile - np.reshape(muRTrain, (1, -1)), axis = 1) denVals = self.trainedModel.getQVal(mus) return np.abs(nodalVals / denVals) def _RHSNorms(self, radiusb0:Np2D) -> Np1D: """High fidelity system RHS norms.""" self.assembleReducedResidualBlocks(full = False) # 'ij,jk,ik->k', resbb, radiusb0, radiusb0.conj() b0resb0 = np.sum(self.trainedModel.data.resbb.dot(radiusb0) * radiusb0.conj(), axis = 0) RHSnorms = np.power(np.abs(b0resb0), .5) return RHSnorms def _errorEstimatorBare(self) -> Np1D: """Bare residual-based error estimator.""" self.setupApprox() self.assembleReducedResidualGramian(self.trainedModel.data.projMat) pDom = self.trainedModel.data.P.domCoeff LL = pDom.conj().dot(self.trainedModel.data.gramian.dot(pDom)) Adiag = self.As[0].diagonal() AnormApprox = np.linalg.norm(Adiag) ** 2. / np.size(Adiag) return np.abs(AnormApprox * LL) ** .5 def _errorEstimatorBasic(self, muTest:paramVal, ratioTest:complex) -> Np1D: """Basic residual-based error estimator.""" resmu = self.HFEngine.residual(self.getApprox(muTest), muTest, self.homogeneized, duality = False)[0] return np.abs(self.estimatorNormEngine.norm(resmu) / ratioTest) def _errorEstimatorExact(self, muRTrain:Np1D, vanderBase:Np2D) -> Np1D: """Exact residual-based error estimator.""" self.setupApprox() self.assembleReducedResidualBlocks(full = True) nAsM = self.HFEngine.nAs - 1 nbsM = max(self.HFEngine.nbs - 1, nAsM * self.homogeneized) radiusA = np.zeros((len(self.mus), nAsM, vanderBase.shape[1]), dtype = np.complex) radiusb = np.zeros((nbsM, vanderBase.shape[1]), dtype = np.complex) domcoeff = polydomcoeff(self.trainedModel.data.Q.deg[0], self.trainedModel.data.Q.polybasis) Qvals = self.trainedModel.getQVal(self.mus) for k in range(max(nAsM, nbsM)): if k < nAsM: radiusA[:, k :, :] += np.multiply.outer(Qvals * self._fitinv, vanderBase[: nAsM - k, :]) if k < nbsM: radiusb[k :, :] += (self._fitinv.dot(Qvals) * vanderBase[: nbsM - k, :]) Qvals = Qvals * muRTrain if self.POD: radiusA = np.tensordot(self.samplingEngine.RPOD, radiusA, 1) # 'ijkl,klt,ijt->t', resAA, radiusA, radiusA.conj() LL = np.sum(np.tensordot(self.trainedModel.data.resAA, radiusA, 2) * radiusA.conj(), axis = (0, 1)) # 'ijk,jkl,il->l', resAb, radiusA, radiusb.conj() Lf = np.sum(np.tensordot(self.trainedModel.data.resAb[1 :, :, :], radiusA, 2) * radiusb.conj(), axis = 0) # 'ij,jk,ik->k', resbb, radiusb, radiusb.conj() ff = np.sum(self.trainedModel.data.resbb[1 :, 1 :].dot(radiusb) * radiusb.conj(), axis = 0) return domcoeff * np.abs(ff - 2. * np.real(Lf) + LL) ** .5 def errorEstimator(self, mus:Np1D) -> Np1D: """Standard residual-based error estimator.""" self.setupApprox() nAsM = self.HFEngine.nAs - 1 nbsM = max(self.HFEngine.nbs - 1, nAsM * self.homogeneized) muRTest = self.centerNormalize(mus)(0) muRTrain = self.centerNormalize(self.mus)(0) samplingRatio = self._errorSamplingRatio(mus, muRTest, muRTrain) if self.errorEstimatorKind == "EXACT": vanSize = max(nAsM, nbsM) else: vanSize = nbsM vanderBase = np.polynomial.polynomial.polyvander(muRTest, vanSize).T RHSnorms = self._RHSNorms(vanderBase[: nbsM + 1, :]) if self.errorEstimatorKind == "BARE": jOpt = self._errorEstimatorBare() elif self.errorEstimatorKind == "BASIC": idx_muTestSample = np.argmax(samplingRatio) jOpt = self._errorEstimatorBasic(mus[idx_muTestSample], samplingRatio[idx_muTestSample]) else: #if self.errorEstimatorKind == "EXACT": jOpt = self._errorEstimatorExact(muRTrain, vanderBase[: -1, :]) return jOpt * samplingRatio / RHSnorms def setupApprox(self, plotEst : bool = False): """ 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.computeScaleFactor() if not hasattr(self, "As") or not hasattr(self, "bs"): vbMng(self, "INIT", "Computing affine blocks of system.", 10) self.As = self.HFEngine.affineLinearSystemA(self.mu0, self.scaleFactor)[1 :] self.bs = self.HFEngine.affineLinearSystemb(self.mu0, self.scaleFactor, self.homogeneized) vbMng(self, "DEL", "Done computing affine blocks.", 10) self.greedy(plotEst) self._S = len(self.mus) self._N, self._M = [self._S - 1] * 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: Qf = self._setupDenominator() Q = Qf[0] - if self.errorEstimatorKind == "EXACT": self._fitinv = Qf[1] + if self.errorEstimatorKind == "EXACT": + self._fitinv = Qf[1].flatten() else: Q = PI() Q.coeffs = np.ones(1, dtype = np.complex) Q.npar = 1 Q.polybasis = self.polybasis if self.errorEstimatorKind == "EXACT": self._fitinv = np.ones(1) self.trainedModel.data.Q = copy(Q) self.trainedModel.data.P = copy(self._setupNumerator()) self.trainedModel.data.approxParameters = copy(self.approxParameters) vbMng(self, "DEL", "Done setting up approximant.", 5) diff --git a/rrompy/utilities/numerical/custom_pinv.py b/rrompy/utilities/numerical/custom_pinv.py index 45b5fbd..3cda020 100644 --- a/rrompy/utilities/numerical/custom_pinv.py +++ b/rrompy/utilities/numerical/custom_pinv.py @@ -1,46 +1,48 @@ # 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 __all__ = ["customPInv"] -def customPInv(A, rcond=1e-15, full=False): +def customPInv(A, rcond=-1, full=False): """ Compute the (Moore-Penrose) pseudo-inverse of a matrix. Calculate the generalized inverse of a matrix using its singular-value decomposition (SVD) and including all *large* singular values. """ A = A.conjugate() u, s, vt = np.linalg.svd(A, full_matrices=False) + if rcond < 0: rcond = len(A) * np.finfo(A.dtype).eps + cutoff = rcond * np.amax(s) large = s > cutoff sinv = copy(s) sinv = np.divide(1, s, where = large, out = sinv) sinv[~large] = 0 res = (vt.T * sinv) @ u.T if full: return res, [np.sum(large), s, rcond] else: return res diff --git a/rrompy/utilities/poly_fitting/polynomial/polynomial_interpolator.py b/rrompy/utilities/poly_fitting/polynomial/polynomial_interpolator.py index 90f8472..7c95b03 100644 --- a/rrompy/utilities/poly_fitting/polynomial/polynomial_interpolator.py +++ b/rrompy/utilities/poly_fitting/polynomial/polynomial_interpolator.py @@ -1,141 +1,140 @@ # 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.base import freepar as fp from rrompy.utilities.poly_fitting.interpolator import GenericInterpolator from rrompy.utilities.poly_fitting.custom_fit import customFit from rrompy.utilities.poly_fitting.polynomial.base import (polyfitname, polydomcoeff) from rrompy.utilities.poly_fitting.polynomial.val import polyval from rrompy.utilities.poly_fitting.polynomial.roots import polyroots from rrompy.utilities.poly_fitting.polynomial.homogeneization import ( homogeneizedpolyvander as hpv, homogeneizedToFull) from rrompy.utilities.poly_fitting.polynomial.vander import polyvander as pv from rrompy.utilities.exception_manager import RROMPyAssert, RROMPyException from rrompy.parameter import checkParameterList __all__ = ['PolynomialInterpolator'] class PolynomialInterpolator(GenericInterpolator): """HERE""" def __init__(self, other = None): if other is None: return self.coeffs = other.coeffs self.npar = other.npar self.polybasis = other.polybasis @property def shape(self): if self.coeffs.ndim > self.npar: sh = self.coeffs.shape[self.npar :] else: sh = tuple([1]) return sh @property def deg(self): return [x - 1 for x in self.coeffs.shape[: self.npar]] def __getitem__(self, key): return self.coeffs[key] def __call__(self, mu:paramList, der : List[int] = None, scl : Np1D = None): return polyval(mu, self.coeffs, self.polybasis, der, scl) def __copy__(self): return PolynomialInterpolator(self) def __deepcopy__(self, memo): other = PolynomialInterpolator() other.coeffs, other.npar, other.polybasis = copy( (self.coeffs, self.npar, self.polybasis), memo) return other @property def domCoeff(self): RROMPyAssert(self.npar, 1, "Number of parameters") return polydomcoeff(self.deg, self.polybasis) * self[-1] 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)] * self.npar padwidth = padwidth + [(l, r) for l, r in zip(nleft, nright)] self.coeffs = np.pad(self.coeffs, padwidth, "constant", constant_values = (0., 0.)) def postmultiplyTensorize(self, A:Np2D): RROMPyAssert(A.shape[0], self.shape[-1], "Shape of output") self.coeffs = np.tensordot(self.coeffs, A, axes = (-1, 0)) def setupByInterpolation(self, support:paramList, values:ListAny, deg:int, polybasis : str = "MONOMIAL", verbose : bool = True, homogeneized : bool = True, vanderCoeffs : DictAny = {}, fitCoeffs : DictAny = {}): support = checkParameterList(support)[0] self.npar = support.shape[1] self.polybasis = polybasis if homogeneized: - vander, _, reorder = hpv(support, deg = deg, basis = polybasis, + vander, _, reorder = hpv(support, deg, basis = polybasis, **vanderCoeffs) vander = vander[:, reorder] else: if not hasattr(deg, "__len__"): deg = [deg] * self.npar - vander = pv(support, deg = deg, basis = polybasis, - **vanderCoeffs) + vander = pv(support, deg, basis = polybasis, **vanderCoeffs) outDim = values.shape[1:] values = values.reshape(values.shape[0], -1) fitOut = customFit(vander, values, full = True, **fitCoeffs) P = fitOut[0] if verbose: msg = ("Fitting {} samples with degree {} through {}... " "Conditioning of LS system: {:.4e}.").format( len(vander), deg, polyfitname(self.polybasis), fitOut[1][2][0] / fitOut[1][2][-1]) else: msg = None if homogeneized: self.coeffs = homogeneizedToFull( tuple([deg + 1] * self.npar) + outDim, self.npar, P) else: self.coeffs = P.reshape(tuple([d + 1 for d in deg]) + outDim) return fitOut[1][1] == vander.shape[1], msg def roots(self, marginalVals : ListAny = [fp]): RROMPyAssert(self.shape, (1,), "Shape of output") RROMPyAssert(len(marginalVals), self.npar, "Number of parameters") try: rDim = marginalVals.index(fp) if rDim < len(marginalVals) - 1 and fp in marginalVals[rDim + 1 :]: raise except: raise RROMPyException(("Exactly 1 'freepar' entry in " "marginalVals must be provided.")) return polyroots(self.coeffs, self.polybasis, marginalVals) 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 cd2584c..7d8593f 100644 --- a/rrompy/utilities/poly_fitting/radial_basis/radial_basis_interpolator.py +++ b/rrompy/utilities/poly_fitting/radial_basis/radial_basis_interpolator.py @@ -1,142 +1,141 @@ # 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 rrompy.utilities.poly_fitting.radial_basis.base import polyfitname from rrompy.utilities.poly_fitting.radial_basis.val import polyval from rrompy.utilities.poly_fitting.radial_basis.homogeneization import ( homogeneizedpolyvander as hpv) from rrompy.utilities.poly_fitting.radial_basis.vander import polyvander as pv from rrompy.utilities.poly_fitting.polynomial.homogeneization import ( homogeneizedToFull) 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 @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) 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) 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, homogeneized : 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"]] self.npar = support.shape[1] if directionalWeights is None: directionalWeights = np.ones(self.npar) self.directionalWeights = directionalWeights self.polybasis = polybasis if homogeneized: - vander, _, reorder = hpv(support, deg = deg, basis = polybasis, + vander, _, reorder = hpv(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 = deg, basis = polybasis, - **vanderCoeffs) + vander = pv(support, deg, basis = polybasis, **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 homogeneized: self.coeffsGlobal = homogeneizedToFull( 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/vander.py b/rrompy/utilities/poly_fitting/radial_basis/vander.py index 086d74f..d1cf051 100644 --- a/rrompy/utilities/poly_fitting/radial_basis/vander.py +++ b/rrompy/utilities/poly_fitting/radial_basis/vander.py @@ -1,79 +1,79 @@ # 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 from rrompy.utilities.base.types import Np1D, Np2D, List, paramList from rrompy.parameter import checkParameterList from rrompy.utilities.exception_manager import RROMPyException, RROMPyAssert from .kernel import radialGaussian, thinPlateSpline, multiQuadric __all__ = ['rbvander', 'polyvander'] def rbvander(x:paramList, basis:str, reorder : List[int] = None, directionalWeights : Np1D = None) -> Np2D: """Compute radial-basis-Vandermonde matrix.""" x = checkParameterList(x)[0] - x_un, idx_un = x.unique(return_inverse = True) + 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]) RROMPyAssert(len(directionalWeights), x.shape[1], "Number of directional weights") try: radialkernel = {"GAUSSIAN" : radialGaussian, "THINPLATE" : thinPlateSpline, "MULTIQUADRIC" : multiQuadric }[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] 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)) 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: """ Compute radial-basis-inclusive 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) 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))]])