diff --git a/rrompy/utilities/poly_fitting/radial_basis/__init__.py b/rrompy/utilities/poly_fitting/radial_basis/__init__.py index 485b731..e4d799d 100644 --- a/rrompy/utilities/poly_fitting/radial_basis/__init__.py +++ b/rrompy/utilities/poly_fitting/radial_basis/__init__.py @@ -1,50 +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 .base import polybases, polyfitname, polydomcoeff, radialFunction +from .base import rbbases, polybases, polyfitname, polydomcoeff, radialFunction from .der import polyder from .val import polyval -#from .vander import polyvander -#from .roots import polyroots -#from .derivative import nextDerivativeIndices -#from .hash_derivative import hashDerivativeToIdx, hashIdxToDerivative -#from .homogeneization import (homogeneizationMask, homogeneizedpolyvander, -# homogeneizedToFull) +from .vander import rbvander, polyvander +from .homogeneization import homogeneizedpolyvander __all__ = [ 'radialGaussian', 'thinPlateSpline', 'multiQuadric', + 'rbbases', 'polybases', 'polyfitname', 'polydomcoeff', 'radialFunction', 'polyder', 'polyval', -# 'polyvander', -# 'polyroots', -# 'nextDerivativeIndices', -# 'hashDerivativeToIdx', -# 'hashIdxToDerivative', -# 'homogeneizationMask', -# 'homogeneizedpolyvander', -# 'homogeneizedToFull' + 'rbvander', + 'polyvander', + 'homogeneizedpolyvander' ] diff --git a/rrompy/utilities/poly_fitting/radial_basis/base.py b/rrompy/utilities/poly_fitting/radial_basis/base.py index d3ea64d..57b49e2 100644 --- a/rrompy/utilities/poly_fitting/radial_basis/base.py +++ b/rrompy/utilities/poly_fitting/radial_basis/base.py @@ -1,55 +1,58 @@ # 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.base.types import Tuple +from rrompy.utilities.base.types import Np1D, Np2D, paramList from rrompy.utilities.exception_manager import RROMPyException from rrompy.utilities.poly_fitting.polynomial.base import polydomcoeff as \ polydomcoeffB -__all__ = ['polybases', 'polyfitname', 'polydomcoeff', 'radialFunction'] +__all__ = ['rbbases', 'polybases', 'polyfitname', 'polydomcoeff', + 'radialFunction'] -polybases = [x + "_" + y for x, y in product( - ["CHEBYSHEV", "LEGENDRE", "MONOMIAL"], - ["GAUSSIAN", "THINPLATE", "MULTIQUADRIC"])] +rbbases = ["GAUSSIAN", "THINPLATE", "MULTIQUADRIC"] -def splitpolybasis(basis:str) -> Tuple[str, str]: - basisUnder = basis.find("_") - return basis[: basisUnder].upper(), basis[basisUnder + 1 :].upper() +polybases = [x + "_" + y for x, y in product( + ["CHEBYSHEV", "LEGENDRE", "MONOMIAL"], rbbases)] def polyfitname(basis:str) -> str: fitpnames = {"CHEBYSHEV" : "chebfit", "LEGENDRE" : "legfit", "MONOMIAL" : "polyfit"} fitrnames = {"GAUSSIAN" : "gaussian", "THINPLATE" : "thinplate", "MULTIQUADRIC" : "multiquadric"} - basisp, basisr = splitpolybasis(basis) + basisp, basisr = basis.split("_") try: - return fitpnames[basisp] + fitrnames[basisr] + return fitpnames[basisp] + "_" + fitrnames[basisr] except: - raise RROMPyException("Polynomial basis not recognized.") + raise RROMPyException("Polynomial-radial basis combination not " + "recognized.") def polydomcoeff(n:int, basis:str) -> float: - return polydomcoeffB(splitpolybasis(basis)[0]) + return polydomcoeffB(n, basis.split("_")[0]) class radialFunction: - supportPoints = None - localBasis = None - globalBasis = None - localCoeffs = None # 1-dimensional vector with rb coefficients - globalCoeffs = None # d-dimensional tensor with poly coefficients - + def __init__(self, supportPoints : paramList = None, + directionalWeights : Np1D = None, localBasis : str = None, + globalBasis : str = None, localCoeffs : Np1D = None, + globalCoeffs : Np2D = None): + self.supportPoints = supportPoints + self.directionalWeights = directionalWeights + self.localBasis = localBasis + self.globalBasis = globalBasis + self.localCoeffs = localCoeffs + self.globalCoeffs = globalCoeffs diff --git a/rrompy/utilities/poly_fitting/radial_basis/der.py b/rrompy/utilities/poly_fitting/radial_basis/der.py index 6826460..89f2aed 100644 --- a/rrompy/utilities/poly_fitting/radial_basis/der.py +++ b/rrompy/utilities/poly_fitting/radial_basis/der.py @@ -1,28 +1,30 @@ # 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, List, radialFun -from rrompy.utilities.exception_manager import RROMPyAssert +from rrompy.utilities.exception_manager import RROMPyException __all__ = ['polyder'] def polyder(c:radialFun, basis:str, m : List[int] = None, scl : Np1D = None) -> radialFun: - RROMPyAssert(np.sum(m), 0., "Radial basis derivative") + if m is not None and np.sum(m) > 0: + raise RROMPyException(("Cannot take derivatives of radial basis " + "function.")) return c diff --git a/rrompy/utilities/poly_fitting/radial_basis/homogeneization.py b/rrompy/utilities/poly_fitting/radial_basis/homogeneization.py new file mode 100644 index 0000000..6fa44b9 --- /dev/null +++ b/rrompy/utilities/poly_fitting/radial_basis/homogeneization.py @@ -0,0 +1,49 @@ +# 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, paramList +from rrompy.utilities.poly_fitting.polynomial.homogeneization import ( + homogeneizedpolyvander as hpvP) +from rrompy.utilities.exception_manager import RROMPyException +from .vander import rbvander + +__all__ = ['homogeneizedpolyvander'] + +def homogeneizedpolyvander(x:paramList, degs:List[int], basis:str, + derIdxs : List[List[List[int]]] = None, + reorder : List[int] = None, + directionalWeights : Np1D = None, + scl : Np1D = None)\ + -> Tuple[Np2D, List[List[int]], List[int]]: + """ + Compute radial-basis-inclusive homogeneized 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, derIdxs, ordIdxs = hpvP(x, degs, basisp, derIdxs = derIdxs, + reorder = reorder, scl = scl) + ordIdxsEff = np.concatenate((np.arange(len(VanR)), ordIdxs + len(VanR))) + return (np.bmat([[VanR, VanP], + [VanP.T.conj(), np.zeros(tuple([VanP.shape[1]] * 2))]]), + derIdxs, ordIdxsEff) diff --git a/rrompy/utilities/poly_fitting/radial_basis/kernel.py b/rrompy/utilities/poly_fitting/radial_basis/kernel.py index 3ef9941..d0e1db2 100644 --- a/rrompy/utilities/poly_fitting/radial_basis/kernel.py +++ b/rrompy/utilities/poly_fitting/radial_basis/kernel.py @@ -1,35 +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 . # import numpy as np from rrompy.utilities.base.types import Np1D from rrompy.utilities.exception_manager import RROMPyAssert __all__ = ['radialGaussian', 'thinPlateSpline', 'multiQuadric'] 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) + return np.power(r2 + 1., -.5) diff --git a/rrompy/utilities/poly_fitting/radial_basis/radial_basis_fitter.py b/rrompy/utilities/poly_fitting/radial_basis/radial_basis_fitter.py deleted file mode 100644 index dec76f5..0000000 --- a/rrompy/utilities/poly_fitting/radial_basis/radial_basis_fitter.py +++ /dev/null @@ -1,229 +0,0 @@ -# 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, List, ListAny, paramList -from rrompy.solver import Np2DLikeEye, normEngine -from rrompy.parameter import checkParameterList -from rrompy.utilities.exception_manager import RROMPyException, RROMPyAssert - -__all__ = ['radialBasisFitter', 'radialGaussian', 'thinPlateSpline', - 'multiQuadric'] - -def radialGaussian(r2): - return np.exp(- r2) - -def thinPlateSpline(r2): - return .5 * r2 * np.log(np.finfo(float).eps + r2) - -def multiQuadric(r2): - return np.power(r2 + 1, .5) - -class radialBasisFitter: - """HERE""" - - allowedModes = ["PARAMETERS", "VALUES"] - - def __init__(self, mus:paramList, basisFun : callable = radialGaussian, - massMatrix : Np2D = None, mode : str = "PARAMETERS", - scl : float = 1.): - self.mus = mus - self.basisFun = basisFun - if massMatrix is None: massMatrix = normEngine(Np2DLikeEye()) - self.massMatrix = massMatrix - self.mode = mode - self.scl = scl - - @property - def d(self): - """Number of parameters.""" - return self.mus.shape[1] - - @property - def n(self): - """Number of parameter points.""" - return len(self.mus) - - @property - def basisFun(self): - """Value of basisFun. Its assignment resets all.""" - return self._basisFun - @basisFun.setter - def basisFun(self, basisFun): - self.reset() - self._basisFun = basisFun - - @property - def mus(self): - """Value of mus. Its assignment resets all.""" - return self._mus - @mus.setter - def mus(self, mus): - mus, _ = checkParameterList(mus) - self.reset() - self._mus = mus - - @property - def massMatrix(self): - """Value of massMatrix. Its assignment resets all.""" - return self._massMatrix - @massMatrix.setter - def massMatrix(self, massMatrix): - self.reset() - self._massMatrix = massMatrix - - @property - def mode(self): - """Value of mode. Its assignment resets all.""" - return self._mode - @mode.setter - def mode(self, mode): - self.reset() - self._mode = mode.upper() - - @property - def scl(self): - """Value of scl. Its assignment resets all.""" - return self._scl - @scl.setter - def scl(self, scl): - self.reset() - self._scl = scl - - def reset(self): - self.vander = None - self.offDiag = None - self.offDiagT = None - self.matrixInv = None - self.probeParameters = None - self.probeValues = None - - def buildMatrixBlocks(self): - if self.offDiag is None: - self.reset() - self.offDiagT = np.array([[1] + list(x[0]) for x in self.mus]) - self.offDiag = self.offDiagT.T - muDiff = np.empty((self.d, self.n * (self.n - 1) // 2 + 1), - dtype = self.mus.dtype) - muDiff[:, 0] = 0. - idxInv = np.zeros(self.n ** 2, dtype = int) - for j in range(self.n): - idx = j * (self.n - 1) - j * (j + 1) // 2 - for i in range(j + 1, self.n): - muDiff[:, idx + i] = (self.offDiag[1:, j] - - self.offDiag[1:, i]) - idxInv[j * self.n + i] = idx + i - idxInv[i * self.n + j] = idx + i - self.vander = self.basisFun(np.power(self.scl * - self.massMatrix.norm(muDiff), 2.))[idxInv] - self.vander = self.vander.reshape((self.n, -1)) - self.vanderProj = self.offDiag.dot(self.vander.dot(self.offDiag.T)) - - def buildMatrixInvBlocks(self): - if self.matrixInv is None: - self.buildMatrixBlocks() - vanderInv = np.linalg.inv(self.vander) - vanderProjInv = np.linalg.solve(self.vanderProj, - self.offDiag.dot(vanderInv)) - self.matrixInv = np.empty((self.n + self.d + 1, self.n), - dtype = vanderProjInv.dtype) - self.matrixInv[self.n :, :] = vanderProjInv - self.matrixInv[: self.n, :] = vanderInv - vanderInv.dot( - self.offDiagT.dot(vanderProjInv)) - - def setupProbeParameters(self, mu : paramList = []) -> List[Np1D]: - mu, _ = checkParameterList(mu, self.d) - self.buildMatrixInvBlocks() - self.probeParameters = [None] * len(mu) - for j, m in enumerate(mu): - probe = np.ones(self.n + self.d + 1, dtype = m.dtype) - probe[self.n + 1 :] = m.data # flatten? - mDiff = (probe[self.n + 1:] - self.offDiagT[:, 1:]).T - probe[: self.n] = self.basisFun(np.power(self.scl * - self.massMatrix.norm(mDiff), 2.)) - self.probeParameters[j] = probe.dot(self.matrixInv) - - def setupProbeValues(self, vals:ListAny) -> ListAny: - RROMPyAssert(len(vals), self.n, "Number of samples") - self.buildMatrixInvBlocks() - if isinstance(vals, (np.ndarray,)): - self.probeValues = np.tensordot(self.matrixInv, vals, - axes = ([-1], [0])) - else: - self.probeValues = [None] * (self.n + self.d + 1) - for j in range(self.n + self.d + 1): - self.probeValues[j] = self.matrixInv[j, 0] * vals[0] - for i in range(1, self.n): - self.probeValues[j] += self.matrixInv[j, i] * vals[i] - - def interpolateParameters(self, vals:ListAny) -> ListAny: - if self.probeParameters is None: - raise RROMPyException(("Parameter probe must be set up before " - "interpolation.")) - RROMPyAssert(len(vals), self.n, "Number of samples") - interpolated = [None] * len(self.probeParameters) - if isinstance(vals, (np.ndarray,)): - if vals.ndim == 1: - for j, pUp in enumerate(self.probeParameters): - interpolated[j] = pUp.dot(vals) - else: - for j, pUp in enumerate(self.probeParameters): - interpolated[j] = np.tensordot(pUp, vals, - axes = ([0], [0])) - else: - for j, pUp in enumerate(self.probeParameters): - interpolated[j] = self.probeParameters[j][0] * vals[0] - for i in range(1, self.n): - interpolated[j] += self.probeParameters[j][i] * vals[i] - return interpolated - - def interpolateValues(self, mu : paramList = []) -> ListAny: - if self.probeValues is None: - raise RROMPyException(("Value probe must be set up before " - "interpolation.")) - mu, _ = checkParameterList(mu, self.d) - probeLs = [None] * len(mu) - for j, m in enumerate(mu): - probeLs[j] = np.ones(self.n + self.d + 1, dtype = m.dtype) - probeLs[j][self.n + 1 :] = m.data # flatten? - mDiff = (probeLs[j][self.n + 1:] - self.offDiagT[:, 1:]).T - probeLs[j][: self.n] = self.basisFun(np.power(self.scl * - self.massMatrix.norm(mDiff), 2.)) - interpolated = [None] * len(mu) - if isinstance(self.probeValues, (np.ndarray,)): - if self.probeValues.ndim == 1: - for j, pL in enumerate(probeLs): - interpolated[j] = pL.dot(self.probeValues) - else: - for j, pL in enumerate(probeLs): - interpolated[j] = np.tensordot(pL, self.probeValues, - axes = ([-1], [0])) - else: - for j, pL in enumerate(probeLs): - interpolated[j] = pL[0] * self.probeValues[0] - for i in range(1, self.n + self.d + 1): - interpolated[j] += pL[i] * self.probeValues[i] - return interpolated - - def interpolate(self, x) -> ListAny: - if self.mode == "PARAMETERS": - return self.interpolateParameters(x) - if self.mode == "VALUES": - return self.interpolateValues(x) - raise RROMPyException("Not implemented") - diff --git a/rrompy/utilities/poly_fitting/radial_basis/val.py b/rrompy/utilities/poly_fitting/radial_basis/val.py index ac927da..66ad07a 100644 --- a/rrompy/utilities/poly_fitting/radial_basis/val.py +++ b/rrompy/utilities/poly_fitting/radial_basis/val.py @@ -1,59 +1,57 @@ # 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 import polyder from rrompy.utilities.base.types import Np1D, Np2D, List, paramList, radialFun from rrompy.parameter import checkParameterList from rrompy.utilities.exception_manager import RROMPyException -from .base import splitpolybasis +from .der import polyder from .kernel import radialGaussian, thinPlateSpline, multiQuadric __all__ = ['polyval'] def polyval(x:paramList, c:radialFun, basis:str, m : List[int] = None, scl : Np1D = None) -> Np2D: cFun = polyder(c, basis, m = m, scl = scl) c = cFun.globalCoeffs x, _ = checkParameterList(x) if x.shape[1] > c.ndim: raise RROMPyException("Incompatible parameter number.") - basisp, basisr = splitpolybasis(basis) + basisp, basisr = basis.split("_") try: polyvalbase = {"CHEBYSHEV" : np.polynomial.chebyshev.chebval, "LEGENDRE" : np.polynomial.legendre.legval, "MONOMIAL" : np.polynomial.polynomial.polyval }[basisp.upper()] except: raise RROMPyException("Polynomial basis not recognized.") try: radialvalbase = {"GAUSSIAN" : radialGaussian, "THINPLATE" : thinPlateSpline, "MULTIQUADRIC" : multiQuadric }[basisr.upper()] except: raise RROMPyException("Radial basis not recognized.") c = polyvalbase(x(0), c, tensor = True) for d in range(1, x.shape[1]): c = polyvalbase(x(d), c, tensor = False) - print(c.shape) - for j, xp in enumerate(x): - muDiff = cFun.supportPoints - xp - c[j] += radialvalbase(np.sum(np.abs(muDiff) ** 2., axis = 1)).dot( - cFun.localCoeffs) + for j in range(len(x)): + muDiff = (cFun.supportPoints.data - x[j]) * cFun.directionalWeights + r2j = np.sum(np.abs(muDiff) ** 2., axis = 1).reshape(1, -1) + c[j] += radialvalbase(r2j).dot(cFun.localCoeffs) return c diff --git a/rrompy/utilities/poly_fitting/radial_basis/vander.py b/rrompy/utilities/poly_fitting/radial_basis/vander.py new file mode 100644 index 0000000..e798b23 --- /dev/null +++ b/rrompy/utilities/poly_fitting/radial_basis/vander.py @@ -0,0 +1,78 @@ +# 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) + 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 + 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) + for j in range(nx): + idx = j * (nx - 1) - j * (j + 1) // 2 + for i in range(j + 1, nx): + r2[idx + i] = np.sum(np.abs((x[j] - x[i]) * directionalWeights + ) ** 2.) + idxInv[j * nx + i] = idx + i + idxInv[i * nx + j] = idx + i + Van = radialkernel(r2)[idxInv].reshape((nx, nx)) + if reorder is not None: Van = Van[reorder, :] + 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.bmat([[VanR, VanP], + [VanP.T.conj(), np.zeros(tuple([VanP.shape[1]] * 2))]]) diff --git a/tests/test_1_utilities/fitting.py b/tests/test_1_utilities/poly_fitting.py similarity index 100% rename from tests/test_1_utilities/fitting.py rename to tests/test_1_utilities/poly_fitting.py diff --git a/tests/test_1_utilities/radial_fitting.py b/tests/test_1_utilities/radial_fitting.py new file mode 100644 index 0000000..1736734 --- /dev/null +++ b/tests/test_1_utilities/radial_fitting.py @@ -0,0 +1,131 @@ +# 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 import customFit +from rrompy.utilities.poly_fitting.radial_basis import (radialGaussian, + thinPlateSpline, + multiQuadric, + polybases, polyfitname, + polydomcoeff, + radialFunction, + polyval, polyvander) +from rrompy.parameter import checkParameterList + +def test_monomial_gaussian(): + polyrbname = "MONOMIAL_GAUSSIAN" + assert polyrbname in polybases + fitname = polyfitname(polyrbname) + domcoeff = polydomcoeff(5, polyrbname) + assert fitname == "polyfit_gaussian" + assert np.isclose(domcoeff, 1., rtol = 1e-5) + + directionalWeights = np.array([5.]) + xSupp = checkParameterList(np.arange(-1, 3), 1)[0] + cRBCoeffs = np.array([-1., 3., -3., 1., 1., 2., -.5]) + cRB = radialFunction(directionalWeights = directionalWeights, + localBasis = "GAUSSIAN", globalBasis = "MONOMIAL", + supportPoints = xSupp, localCoeffs = cRBCoeffs[: 4], + globalCoeffs = cRBCoeffs[4 :]) + ySupp = 1 + 2. * xSupp.data - .5 * xSupp.data ** 2. + xx = np.linspace(-2., 3., 100) + yy = polyval(checkParameterList(xx, 1)[0], cRB, polyrbname) + yyman = 1. + 2. * xx - .5 * xx ** 2. + for j, xc in enumerate(np.arange(-1, 3)): + r2j = (5. * (xx - xc)) ** 2. + rbj = radialGaussian(r2j) + assert np.allclose(rbj, np.exp(-.5 * r2j)) + yyman += cRB.localCoeffs[j] * rbj + ySupp += cRB.localCoeffs[j] * radialGaussian((directionalWeights[0] + * (xSupp.data - xc)) ** 2.) + assert np.allclose(yy, yyman, atol = 1e-5) + + VanT = polyvander(xSupp, [2], polyrbname, + directionalWeights = directionalWeights) + ySupp = np.pad(ySupp.flatten(), (0, len(VanT) - len(xSupp)), "constant") + out = customFit(VanT, ySupp) + assert np.allclose(out, cRBCoeffs, atol = 1e-8) + +def test_legendre_thinplate(): + polyrbname = "LEGENDRE_THINPLATE" + assert polyrbname in polybases + fitname = polyfitname(polyrbname) + domcoeff = polydomcoeff(5, polyrbname) + assert fitname == "legfit_thinplate" + assert np.isclose(domcoeff, 63. / 8, rtol = 1e-5) + + directionalWeights = np.array([.5]) + xSupp = checkParameterList(np.arange(-1, 3), 1)[0] + cRBCoeffs = np.array([-1., 3., -3., 1., 1., 2., -.5]) + cRB = radialFunction(directionalWeights = directionalWeights, + localBasis = "THINPLATE", globalBasis = "LEGENDRE", + supportPoints = xSupp, localCoeffs = cRBCoeffs[: 4], + globalCoeffs = cRBCoeffs[4 :]) + ySupp = 1 + 2. * xSupp.data - .5 * (.5 * (3. * xSupp.data ** 2. - 1.)) + xx = np.linspace(-2., 3., 100) + yy = polyval(checkParameterList(xx, 1)[0], cRB, polyrbname) + yyman = 1. + 2. * xx - .5 * (.5 * (3. * xx ** 2. - 1.)) + for j, xc in enumerate(np.arange(-1, 3)): + r2j = (directionalWeights[0] * (xx - xc)) ** 2. + rbj = thinPlateSpline(r2j) + assert np.allclose(rbj, .5 * r2j * np.log(np.finfo(float).eps + r2j)) + yyman += cRB.localCoeffs[j] * rbj + ySupp += cRB.localCoeffs[j] * thinPlateSpline((directionalWeights[0] + * (xSupp.data - xc)) ** 2.) + assert np.allclose(yy, yyman, atol = 1e-5) + + VanT = polyvander(xSupp, [2], polyrbname, + directionalWeights = directionalWeights) + ySupp = np.pad(ySupp.flatten(), (0, len(VanT) - len(xSupp)), "constant") + out = customFit(VanT, ySupp) + assert np.allclose(out, cRBCoeffs, atol = 1e-8) + +def test_chebyshev_multiquadric(): + polyrbname = "CHEBYSHEV_MULTIQUADRIC" + assert polyrbname in polybases + fitname = polyfitname(polyrbname) + domcoeff = polydomcoeff(5, polyrbname) + assert fitname == "chebfit_multiquadric" + assert np.isclose(domcoeff, 16, rtol = 1e-5) + + directionalWeights = np.array([1.]) + xSupp = checkParameterList(np.arange(-1, 3), 1)[0] + cRBCoeffs = np.array([-1., 3., -3., 1., 1., 2., -.5]) + cRB = radialFunction(directionalWeights = directionalWeights, + localBasis = "MULTIQUADRIC", + globalBasis = "CHEBYSHEV", supportPoints = xSupp, + localCoeffs = cRBCoeffs[: 4], + globalCoeffs = cRBCoeffs[4 :]) + ySupp = 1 + 2. * xSupp.data - .5 * (2. * xSupp.data ** 2. - 1.) + xx = np.linspace(-2., 3., 100) + yy = polyval(checkParameterList(xx, 1)[0], cRB, polyrbname) + yyman = 1. + 2. * xx - .5 * (2. * xx ** 2. - 1.) + for j, xc in enumerate(np.arange(-1, 3)): + r2j = (directionalWeights[0] * (xx - xc)) ** 2. + rbj = multiQuadric(r2j) + assert np.allclose(rbj, np.power(r2j + 1, -.5)) + yyman += cRB.localCoeffs[j] * rbj + ySupp += cRB.localCoeffs[j] * multiQuadric((directionalWeights[0] + * (xSupp.data - xc)) ** 2.) + assert np.allclose(yy, yyman, atol = 1e-5) + + VanT = polyvander(xSupp, [2], polyrbname, + directionalWeights = directionalWeights) + ySupp = np.pad(ySupp.flatten(), (0, len(VanT) - len(xSupp)), "constant") + out = customFit(VanT, ySupp) + assert np.allclose(out, cRBCoeffs, atol = 1e-8)