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)