diff --git a/rrompy/utilities/base/decorators.py b/rrompy/utilities/base/decorators.py
index 6199200..c0cb92a 100644
--- a/rrompy/utilities/base/decorators.py
+++ b/rrompy/utilities/base/decorators.py
@@ -1,27 +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 .
#
-__all__ = ['affine_construct', 'nonaffine_construct']
+from numpy import inf
+
+__all__ = ['affine_construct', 'nonaffine_construct', 'threshold']
def affine_construct(method):
method.is_affine = True
return method
def nonaffine_construct(method):
method.is_affine = False
return method
+
+def threshold(thresh = inf):
+ def method_with_thresh(method):
+ method.threshold = thresh
+ return method
+ return method_with_thresh
diff --git a/rrompy/utilities/poly_fitting/polynomial/base.py b/rrompy/utilities/poly_fitting/polynomial/base.py
index 02562a6..ba549c9 100644
--- a/rrompy/utilities/poly_fitting/polynomial/base.py
+++ b/rrompy/utilities/poly_fitting/polynomial/base.py
@@ -1,61 +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 scipy.special import binom
from rrompy.utilities.exception_manager import RROMPyException
__all__ = ['polybases', 'polyfitname', 'polydomcoeff']
polybases = ["CHEBYSHEV", "LEGENDRE", "MONOMIAL"]
def polyfitname(basis:str) -> str:
- fitnames = {"CHEBYSHEV" : "chebfit", "LEGENDRE" : "legfit",
- "MONOMIAL" : "polyfit"}
- try:
- return fitnames[basis.upper()]
- except:
+ if basis.upper() not in polybases:
raise RROMPyException("Polynomial basis not recognized.")
+ return {"CHEBYSHEV" : "chebfit", "LEGENDRE" : "legfit",
+ "MONOMIAL" : "polyfit"}[basis.upper()]
def polydomcoeff(n:int, basis:str) -> float:
basis = basis.upper()
if isinstance(n, (list, tuple, np.ndarray,)):
nv = np.array(n)
else:
nv = np.array([n])
if basis == "CHEBYSHEV":
x = np.ones_like(nv, dtype = float)
x[nv > 0] = np.power(2., nv[nv > 0] - 1)
elif basis == "LEGENDRE":
x = np.ones_like(nv, dtype = float)
x[nv > 10] = (np.power(2., nv[nv > 10])
* np.power(np.pi * nv[nv > 10], -.5))
x[nv <= 10] = (np.power(.5, nv[nv <= 10])
* binom(2 * nv[nv <= 10], nv[nv <= 10]))
elif basis == "MONOMIAL":
x = np.ones_like(nv, dtype = float)
else:
raise RROMPyException("Polynomial basis not recognized.")
if isinstance(n, (list,)):
return list(x)
if isinstance(n, (tuple,)):
return tuple(x)
if isinstance(n, (np.ndarray,)):
return x
return x[0]
diff --git a/rrompy/utilities/poly_fitting/polynomial/der.py b/rrompy/utilities/poly_fitting/polynomial/der.py
index 092d0ff..7f9a4de 100644
--- a/rrompy/utilities/poly_fitting/polynomial/der.py
+++ b/rrompy/utilities/poly_fitting/polynomial/der.py
@@ -1,43 +1,43 @@
# 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 .base import polybases
from rrompy.utilities.base.types import Np1D, Np2D, List
from rrompy.utilities.exception_manager import RROMPyException, RROMPyAssert
__all__ = ['polyder']
def polyder(c:Np2D, basis:str, m : List[int] = None,
scl : Np1D = None) -> Np2D:
c = np.array(c, ndmin = 1, copy = 1)
if m is not None:
if scl is None: scl = [1.] * c.ndim
if not isinstance(m, (list,tuple,np.ndarray,)): m = [m]
if not isinstance(scl, (list,tuple,np.ndarray,)): scl = [scl]
RROMPyAssert(c.ndim, len(m), "Number of parameters")
RROMPyAssert(c.ndim, len(scl), "Number of parameters")
- try:
- derbase = {"CHEBYSHEV" : np.polynomial.chebyshev.chebder,
- "LEGENDRE" : np.polynomial.legendre.legder,
- "MONOMIAL" : np.polynomial.polynomial.polyder
- }[basis.upper()]
- except:
+ if basis.upper() not in polybases:
raise RROMPyException("Polynomial basis not recognized.")
+ derbase = {"CHEBYSHEV" : np.polynomial.chebyshev.chebder,
+ "LEGENDRE" : np.polynomial.legendre.legder,
+ "MONOMIAL" : np.polynomial.polynomial.polyder
+ }[basis.upper()]
for j, (mj, sj) in enumerate(zip(m, scl)):
c = derbase(c, m = mj, scl = sj, axis = j)
return c
diff --git a/rrompy/utilities/poly_fitting/polynomial/marginalize.py b/rrompy/utilities/poly_fitting/polynomial/marginalize.py
index 67ce917..66da859 100644
--- a/rrompy/utilities/poly_fitting/polynomial/marginalize.py
+++ b/rrompy/utilities/poly_fitting/polynomial/marginalize.py
@@ -1,58 +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 numpy import array, polynomial as po
from copy import deepcopy as copy
+from .base import polybases
from rrompy.utilities.base.types import Np1D, Np2D
from rrompy.utilities.base import freepar as fp
from rrompy.utilities.exception_manager import RROMPyAssert, RROMPyException
__all__ = ['polymarginalize']
def polymarginalize(c:Np2D, basis:str, marginalVals : Np1D = [fp],
nMarginal : int = None) -> Np1D:
if not hasattr(c, "ndim"): c = array(c)
ndim = c.ndim
if not hasattr(marginalVals, "__len__"): marginalVals = [marginalVals]
marginalVals = list(marginalVals)
- try:
- polyvalbase = {"CHEBYSHEV" : po.chebyshev.chebval,
- "LEGENDRE" : po.legendre.legval,
- "MONOMIAL" : po.polynomial.polyval}[basis.upper()]
- except:
+ if basis.upper() not in polybases:
raise RROMPyException("Polynomial basis not recognized.")
+ polyvalbase = {"CHEBYSHEV" : po.chebyshev.chebval,
+ "LEGENDRE" : po.legendre.legval,
+ "MONOMIAL" : po.polynomial.polyval}[basis.upper()]
RROMPyAssert(ndim, len(marginalVals), "Marginalized variables")
marginalDims = []
for j in range(len(marginalVals)):
if marginalVals[j] == fp:
marginalDims += [c.shape[j]]
if nMarginal is not None and len(marginalDims) != nMarginal:
raise RROMPyException(("Exactly {} 'freepar' entries in marginalVals "
"must be provided.").format(nMarginal))
cEff = [copy(c)]
for d in range(ndim):
if marginalVals[d] != fp:
for dj in range(len(cEff)):
cEff[dj] = polyvalbase(marginalVals[d], cEff[dj],
tensor = False)
else:
cEff2 = []
for dj in range(len(cEff)):
cEff2 += list(cEff[dj])
cEff = copy(cEff2)
return array(cEff).reshape(tuple(marginalDims))
diff --git a/rrompy/utilities/poly_fitting/polynomial/roots.py b/rrompy/utilities/poly_fitting/polynomial/roots.py
index f392eb7..dbd1e9a 100644
--- a/rrompy/utilities/poly_fitting/polynomial/roots.py
+++ b/rrompy/utilities/poly_fitting/polynomial/roots.py
@@ -1,34 +1,34 @@
# 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 numpy import polynomial as po
+from .base import polybases
+from .marginalize import polymarginalize
from rrompy.utilities.base.types import Np1D, Np2D
from rrompy.utilities.base import freepar as fp
-from .marginalize import polymarginalize
from rrompy.utilities.exception_manager import RROMPyException
__all__ = ['polyroots']
def polyroots(c:Np2D, basis:str, marginalVals : Np1D = [fp]) -> Np1D:
- try:
- rootsbase = {"CHEBYSHEV" : po.chebyshev.chebroots,
- "LEGENDRE" : po.legendre.legroots,
- "MONOMIAL" : po.polynomial.polyroots}[basis.upper()]
- except:
+ if basis.upper() not in polybases:
raise RROMPyException("Polynomial basis not recognized.")
+ rootsbase = {"CHEBYSHEV" : po.chebyshev.chebroots,
+ "LEGENDRE" : po.legendre.legroots,
+ "MONOMIAL" : po.polynomial.polyroots}[basis.upper()]
return rootsbase(polymarginalize(c, basis, marginalVals, 1))
diff --git a/rrompy/utilities/poly_fitting/polynomial/val.py b/rrompy/utilities/poly_fitting/polynomial/val.py
index 71d651c..6bd85c6 100644
--- a/rrompy/utilities/poly_fitting/polynomial/val.py
+++ b/rrompy/utilities/poly_fitting/polynomial/val.py
@@ -1,43 +1,43 @@
# 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 .base import polybases
from .der import polyder
from rrompy.utilities.base.types import Np1D, Np2D, List, paramList
from rrompy.parameter import checkParameterList
from rrompy.utilities.exception_manager import RROMPyException
__all__ = ['polyval']
def polyval(x:paramList, c:Np2D, basis:str, m : List[int] = None,
scl : Np1D = None) -> Np2D:
c = polyder(c, basis, m = m, scl = scl)
x = checkParameterList(x)
if x.shape[1] > c.ndim:
raise RROMPyException("Incompatible parameter number.")
- try:
- polyvalbase = {"CHEBYSHEV" : np.polynomial.chebyshev.chebval,
- "LEGENDRE" : np.polynomial.legendre.legval,
- "MONOMIAL" : np.polynomial.polynomial.polyval
- }[basis.upper()]
- except:
+ if basis.upper() not in polybases:
raise RROMPyException("Polynomial basis not recognized.")
+ polyvalbase = {"CHEBYSHEV" : np.polynomial.chebyshev.chebval,
+ "LEGENDRE" : np.polynomial.legendre.legval,
+ "MONOMIAL" : np.polynomial.polynomial.polyval
+ }[basis.upper()]
c = polyvalbase(x(0), c, tensor = True)
for d in range(1, x.shape[1]):
c = polyvalbase(x(d), c, tensor = False)
return c
diff --git a/rrompy/utilities/poly_fitting/polynomial/vander.py b/rrompy/utilities/poly_fitting/polynomial/vander.py
index 0d3dff7..d369b11 100644
--- a/rrompy/utilities/poly_fitting/polynomial/vander.py
+++ b/rrompy/utilities/poly_fitting/polynomial/vander.py
@@ -1,138 +1,132 @@
# 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 .base import polybases
from .der import polyder
from rrompy.utilities.base.types import Np1D, Np2D, List, paramList
from rrompy.utilities.numerical.degree import totalDegreeSet
from rrompy.parameter import checkParameterList
from rrompy.utilities.exception_manager import RROMPyException, RROMPyAssert
__all__ = ['polyvander', 'polyvanderTotal']
def firstDerTransition(dim:int, TDirac:List[Np2D], basis:str,
scl : Np1D = None) -> Np2D:
C_m = np.zeros((dim, len(TDirac), len(TDirac)), dtype = float)
for j, Tj in enumerate(TDirac):
m, om = [0] * dim, [(0, 0)] * dim
for idx in range(dim):
m[idx], om[idx] = 1, (0, 1)
J_der = polyder(Tj, basis, m, scl)
if J_der.size != len(TDirac):
J_der = np.pad(J_der, mode = "constant", pad_width = om)
C_m[idx, :, j] = np.ravel(J_der)
m[idx], om[idx] = 0, (0, 0)
return C_m
-
-def countDerDirections(n:int, base:int, digits:int, idx:int):
- if digits == 0: return []
- dig = n % base
- return [(idx, dig)] * (dig > 0) + countDerDirections(
- (n - dig) // base, base, digits - 1, idx + 1)
def polyvander(x:paramList, degs:List[int], basis:str,
derIdxs : List[List[List[int]]] = None,
reorder : List[int] = None, scl : Np1D = None) -> Np2D:
"""
Compute full Hermite-Vandermonde matrix with specified derivative
directions.
E.g. assume that we want to obtain the Vandermonde matrix for
(value, derx, derx2) at x = [0, 0],
(value, dery) at x = [1, 0],
(dery, derxy) at x = [0, 0],
of degree 3 in x and 1 in y, using Chebyshev polynomials.
This can be done by
polyvander([[0, 0], [1, 0]], # unique sample points
[3, 1], # polynomial degree
"chebyshev", # polynomial family
[
[[0, 0], [1, 0], [2, 0], [0, 1], [1, 1]],
# derivative directions at first point
[[0, 0], [0, 1]] # derivative directions at second point
],
[0, 1, 2, 5, 6, 3, 4] # reorder indices
)
"""
if not isinstance(degs, (list,tuple,np.ndarray,)): degs = [degs]
dim = len(degs)
x = checkParameterList(x, dim)
x_un, idx_un = x.unique(return_inverse = True)
if len(x_un) < len(x):
raise RROMPyException("Sample points must be distinct.")
del x_un
- try:
- vanderbase = {"CHEBYSHEV" : np.polynomial.chebyshev.chebvander,
- "LEGENDRE" : np.polynomial.legendre.legvander,
- "MONOMIAL" : np.polynomial.polynomial.polyvander
- }[basis.upper()]
- except:
+ if basis.upper() not in polybases:
raise RROMPyException("Polynomial basis not recognized.")
+ vanderbase = {"CHEBYSHEV" : np.polynomial.chebyshev.chebvander,
+ "LEGENDRE" : np.polynomial.legendre.legvander,
+ "MONOMIAL" : np.polynomial.polynomial.polyvander
+ }[basis.upper()]
VanBase = vanderbase(x(0), degs[0])
for j in range(1, dim):
VNext = vanderbase(x(j), degs[j])
for jj in range(j): VNext = np.expand_dims(VNext, 1)
VanBase = VanBase[..., None] * VNext
VanBase = VanBase.reshape((len(x), -1))
if derIdxs is None or VanBase.shape[-1] <= 1:
Van = VanBase
else:
derFlat, idxRep = [], []
for j, derIdx in enumerate(derIdxs):
derFlat += derIdx[:]
idxRep += [j] * len(derIdx[:])
for j in range(len(derFlat)):
if not hasattr(derFlat[j], "__len__"):
derFlat[j] = [derFlat[j]]
RROMPyAssert(len(derFlat[j]), dim, "Number of dimensions")
TDirac = [y.reshape([d + 1 for d in degs])
for y in np.eye(VanBase.shape[-1], dtype = int)]
Cs_loc = firstDerTransition(dim, TDirac, basis, scl)
Van = np.empty((len(derFlat), VanBase.shape[-1]),
dtype = VanBase.dtype)
for j in range(len(derFlat)):
Van[j, :] = VanBase[idxRep[j], :]
for k in range(dim):
for der in range(derFlat[j][k]):
Van[j, :] = Van[j, :].dot(Cs_loc[k]) / (der + 1)
if reorder is not None: Van = Van[reorder, :]
return Van
def polyvanderTotal(x:paramList, deg:int, basis:str,
derIdxs : List[List[List[int]]] = None,
reorder : List[int] = None, scl : Np1D = None) -> Np2D:
"""
Compute full total degree Hermite-Vandermonde matrix with specified
derivative directions.
"""
x = checkParameterList(x)
degs = [deg] * x.shape[1]
VanBase = polyvander(x, degs, basis, derIdxs, reorder, scl)
derIdxs, mask = totalDegreeSet(deg, x.shape[1], return_mask = True)
ordIdxs = np.empty(len(derIdxs), dtype = int)
derTotal = np.array([np.sum(y) for y in derIdxs])
idxPrev = 0
rangeAux = np.arange(len(derIdxs))
for j in range(deg + 1):
idxLocal = rangeAux[derTotal == j][::-1]
idxPrev += len(idxLocal)
ordIdxs[idxPrev - len(idxLocal) : idxPrev] = idxLocal
return VanBase[:, mask][:, ordIdxs]
diff --git a/rrompy/utilities/poly_fitting/radial_basis/__init__.py b/rrompy/utilities/poly_fitting/radial_basis/__init__.py
index 98a445b..5b03c20 100644
--- a/rrompy/utilities/poly_fitting/radial_basis/__init__.py
+++ b/rrompy/utilities/poly_fitting/radial_basis/__init__.py
@@ -1,42 +1,38 @@
# 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,
- localWendland)
+from .kernel import kernels
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',
- 'localWendland',
+ 'kernels',
'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 2423bf0..d43f306 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 .kernel import kernels
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", "WENDLAND"]
+rbbases = list(kernels.keys())
polybases = [x + "_" + y for x, y in product(pbP, rbbases)]
def polyfitname(basis:str) -> str:
- fitrnames = {"GAUSSIAN" : "gaussian", "THINPLATE" : "thinplate",
- "MULTIQUADRIC" : "multiquadric", "WENDLAND" : "wendland"}
- basisp, basisr = basis.split("_")
try:
- return pfnP(basisp) + "_" + fitrnames[basisr]
+ basisp, basisr = basis.split("_")
+ if basisr.upper() in rbbases: basisr = basisr.lower()
+ else: raise
+ return pfnP(basisp) + "_" + 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 e349e58..2ebd45c 100644
--- a/rrompy/utilities/poly_fitting/radial_basis/kernel.py
+++ b/rrompy/utilities/poly_fitting/radial_basis/kernel.py
@@ -1,42 +1,51 @@
# 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.decorators import threshold
from rrompy.utilities.base.types import Np1D
from rrompy.utilities.exception_manager import RROMPyAssert
-__all__ = ['radialGaussian', 'thinPlateSpline', 'multiQuadric',
- 'localWendland']
+__all__ = ['kernels']
-def radialGaussian(r2:Np1D, der : int = 0) -> Np1D:
+thresholdGaussian = -2. * np.log(np.finfo(float).eps)
+@threshold(thresholdGaussian)
+def radialGaussian(r2:Np1D, der : int = 0,
+ apply_threshold : bool = True) -> Np1D:
RROMPyAssert(der, 0, "Radial basis derivative")
- return np.exp(- .5 * r2)
+ if apply_threshold: r2[r2 > thresholdGaussian] = thresholdGaussian
+ return np.exp(-.5 * r2)
-def thinPlateSpline(r2:Np1D, der : int = 0) -> Np1D:
+thresholdMultiQuadric = np.finfo(float).eps ** -2.
+@threshold(thresholdMultiQuadric)
+def multiQuadric(r2:Np1D, der : int = 0,
+ apply_threshold : bool = True) -> Np1D:
RROMPyAssert(der, 0, "Radial basis derivative")
- return .5 * r2 * np.log(np.finfo(float).eps + r2)
+ if apply_threshold: r2[r2 > thresholdMultiQuadric] = thresholdMultiQuadric
+ return (r2 + 1.) ** -.5
-def multiQuadric(r2:Np1D, der : int = 0) -> Np1D:
- RROMPyAssert(der, 0, "Radial basis derivative")
- return np.power(r2 + 1., -.5)
-
-def localWendland(r2:Np1D, der : int = 0) -> Np1D:
+@threshold(1.)
+def localWendland(r2:Np1D, der : int = 0,
+ apply_threshold : bool = True) -> Np1D:
RROMPyAssert(der, 0, "Radial basis derivative")
+ if apply_threshold: r2[r2 > 1.] = 1.
rm1 = 1. - r2 ** .5
- rm1[rm1 <= 0.] = 0.
return rm1 ** 4. * (5. - 4. * rm1)
+
+kernels = {"GAUSSIAN" : radialGaussian, "MULTIQUADRIC" : multiQuadric,
+ "WENDLAND" : localWendland}
diff --git a/rrompy/utilities/poly_fitting/radial_basis/val.py b/rrompy/utilities/poly_fitting/radial_basis/val.py
index c1b62ac..48da0d8 100644
--- a/rrompy/utilities/poly_fitting/radial_basis/val.py
+++ b/rrompy/utilities/poly_fitting/radial_basis/val.py
@@ -1,54 +1,54 @@
# 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,
- localWendland)
+from .kernel import kernels
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:
x = checkParameterList(x)
basisp, basisr = basis.split("_")
c = pvP(x, cG, basisp)
+ x = x.data
try:
- radialvalbase = {"GAUSSIAN" : radialGaussian,
- "THINPLATE" : thinPlateSpline,
- "MULTIQUADRIC" : multiQuadric,
- "WENDLAND" : localWendland}[basisr.upper()]
+ radialker = kernels[basisr.upper()]
except:
raise RROMPyException("Radial basis not recognized.")
supportPoints = checkParameterList(supportPoints)
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)
- val = np.tensordot(radialvalbase(r2j), cL, 1)
- try:
- c[..., j] += val
- except:
- c[..., j] += val.flatten()
+ radialVal = np.zeros((len(supportPoints), len(x)))
+ xDiff2V, xDiff2I, xDiff2J = [], [], []
+ for i in range(len(supportPoints)):
+ xiD2Loc = np.sum(np.abs((x - supportPoints[i])
+ * directionalWeights) ** 2., axis = 1)
+ xiD2Good = np.where(xiD2Loc <= radialker.threshold)[0]
+ xDiff2V = np.append(xDiff2V, xiD2Loc[xiD2Good])
+ xDiff2I += [i] * len(xiD2Good)
+ xDiff2J += list(xiD2Good)
+ radialVal[xDiff2I, xDiff2J] = radialker(xDiff2V, apply_threshold = False)
+ c += np.tensordot(cL, radialVal, (0, 0))
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 f5b0d8f..e76ab75 100644
--- a/rrompy/utilities/poly_fitting/radial_basis/vander.py
+++ b/rrompy/utilities/poly_fitting/radial_basis/vander.py
@@ -1,97 +1,101 @@
# 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,
- localWendland)
+from .kernel import kernels
from rrompy.utilities.poly_fitting.polynomial.vander import (polyvander as pvP,
polyvanderTotal as pvTP)
from rrompy.utilities.base.types import Np1D, Np2D, 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:
"""Compute radial-basis-Vandermonde matrix."""
x = checkParameterList(x)
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,
- "WENDLAND" : localWendland}[basis.upper()]
+ radialker = kernels[basis.upper()]
except:
raise RROMPyException("Radial basis not recognized.")
Van = np.zeros((nx, nx))
if reorder is not None: x = x[reorder]
- for j in range(nx):
- muDiff = (x - x[j]) * directionalWeights
- r2j = np.sum(np.abs(muDiff) ** 2., axis = 1).reshape(1, -1)
- Van[j] = radialkernel(r2j)
+ xDiff2V, xDiff2I, xDiff2J = [], [], []
+ for i in range(nx - 1):
+ xiD2Loc = np.sum(np.abs((x[i + 1 :] - x[i])
+ * directionalWeights) ** 2., axis = 1)
+ xiD2Good = np.where(xiD2Loc <= radialker.threshold)[0]
+ xDiff2V = np.append(xDiff2V, xiD2Loc[xiD2Good])
+ xDiff2I += [i] * len(xiD2Good)
+ xDiff2J += list(i + 1 + xiD2Good)
+ kernelV = radialker(xDiff2V, apply_threshold = False)
+ Van = np.eye(nx)
+ Van[xDiff2I, xDiff2J] = kernelV
+ Van[xDiff2J, xDiff2I] = kernelV
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 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)
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) -> Np2D:
"""
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)
VanP = pvTP(x, deg, basisp, derIdxs = derIdxs, reorder = reorder, scl = scl)
return np.block([[VanR, VanP],
[VanP.T.conj(), np.zeros(tuple([VanP.shape[1]] * 2))]])
diff --git a/tests/1_utilities/radial_fitting.py b/tests/1_utilities/radial_fitting.py
index c6b4fa3..dd00c5c 100644
--- a/tests/1_utilities/radial_fitting.py
+++ b/tests/1_utilities/radial_fitting.py
@@ -1,165 +1,130 @@
# 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,
+from rrompy.utilities.poly_fitting.radial_basis import (kernels,
polybases, polyfitname,
polydomcoeff,
polyval, polyvander,
polyvanderTotal)
from rrompy.utilities.numerical.degree import degreeTotalToFull
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)
+ radialGaussian = kernels["GAUSSIAN"]
directionalWeights = np.array([5.])
xSupp = checkParameterList(np.arange(-1, 3), 1)
cRBCoeffs = np.array([-1., 3., -3., 1., 1., 2., -.5])
globalCoeffs = cRBCoeffs[4 :]
localCoeffs = cRBCoeffs[: 4]
ySupp = 1 + 2. * xSupp.data - .5 * xSupp.data ** 2.
xx = np.linspace(-2., 3., 100)
yy = polyval(checkParameterList(xx, 1), globalCoeffs, localCoeffs,
xSupp, directionalWeights, 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 += localCoeffs[j] * rbj
ySupp += 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)
- cRBCoeffs = np.array([-1., 3., -3., 1., 1., 2., -.5])
-
- 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), globalCoeffs, localCoeffs,
- xSupp, directionalWeights, 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 += localCoeffs[j] * rbj
- ySupp += 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)
+ multiQuadric = kernels["MULTIQUADRIC"]
directionalWeights = np.array([1.])
xSupp = checkParameterList(np.arange(-1, 3), 1)
cRBCoeffs = np.array([-1., 3., -3., 1., 1., 2., -.5])
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), globalCoeffs, localCoeffs,
xSupp, directionalWeights, 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 += localCoeffs[j] * rbj
ySupp += 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)
def test_total_degree_2d():
values = lambda x, y: (x - 3.) ** 2. * y - (x + 1.) * y ** 2.
polyrbname = "CHEBYSHEV_GAUSSIAN"
xs, ys = np.meshgrid(np.linspace(0., 4., 5), np.linspace(0., 4., 4))
xySupp = np.concatenate((xs.reshape(-1, 1), ys.reshape(-1, 1)), axis = 1)
zs = values(xs, ys)
zSupp = zs.flatten()
deg = 3
directionalWeights = [2., 1.]
VanT = polyvanderTotal(xySupp, deg, polyrbname,
directionalWeights = directionalWeights)
cFit = np.linalg.solve(VanT, np.pad(zSupp, (0, len(VanT) - len(zSupp)),
'constant'))
globCoeff = degreeTotalToFull([deg + 1] * 2, 2, cFit[len(zSupp) :])
localCoeffs = cFit[: len(zSupp)]
globalCoeffs = globCoeff
xx, yy = np.meshgrid(np.linspace(0., 4., 100), np.linspace(0., 4., 100))
xxyy = np.concatenate((xx.reshape(-1, 1), yy.reshape(-1, 1)), axis = 1)
zz = polyval(xxyy, globalCoeffs, localCoeffs, xySupp, directionalWeights,
polyrbname).reshape(xx.shape)
zzex = values(xx, yy)
error = np.abs(zz - zzex)
print(np.max(error))
assert np.max(error) < 1e-10