diff --git a/rrompy/utilities/poly_fitting/heaviside/vander.py b/rrompy/utilities/poly_fitting/heaviside/vander.py
index ce816d3..3d0a2f7 100644
--- a/rrompy/utilities/poly_fitting/heaviside/vander.py
+++ b/rrompy/utilities/poly_fitting/heaviside/vander.py
@@ -1,89 +1,87 @@
# Copyright (C) 2018 by the RROMPy authors
#
# This file is part of RROMPy.
#
# RROMPy is free software: you can redistribute it and/or modify
# it under the terms of the GNU Lesser General Public License as published by
# the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
#
# RROMPy is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU Lesser General Public License for more details.
#
# You should have received a copy of the GNU Lesser General Public License
# along with RROMPy. If not, see .
#
import numpy as np
from rrompy.utilities.poly_fitting.polynomial.vander import (polyvander as pvP,
polyvanderTotal as pvTP)
from rrompy.utilities.base.types import Np1D, Np2D, Tuple, List, paramList
from rrompy.parameter import checkParameterList
from rrompy.utilities.exception_manager import RROMPyException
__all__ = ['heavisidevander', 'polyvander', 'polyvanderTotal']
def heavisidevander(x:paramList, poles:Np1D,
reorder : List[int] = None) -> Np2D:
"""Compute Heaviside-Vandermonde matrix."""
x = checkParameterList(x, 1)[0]
x_un, idx_un = x.unique(return_inverse = True)
nx = len(x)
if len(x_un) < nx:
raise RROMPyException("Sample points must be distinct.")
del x_un
x = x.data.flatten()
if reorder is not None: x = x[reorder]
poles = poles.flatten()
Van = np.empty((len(x), len(poles)), dtype = np.complex)
for j in range(len(x)):
Van[j, :] = 1. / (x[j] - poles)
return Van
def polyvander(x:paramList, poles:Np1D, degs : List[int] = None,
basis : str = "MONOMIAL_HEAVISIDE",
derIdxs : List[List[List[int]]] = None,
reorder : List[int] = None, scl : Np1D = None) -> Np2D:
"""
Compute full Hermite-Vandermonde matrix with specified derivative
directions.
"""
if derIdxs is not None and np.sum(np.sum(derIdxs)) > 0:
raise RROMPyException(("Cannot take derivatives of heaviside "
"function."))
basisp = basis.split("_")[0]
VanH = heavisidevander(x, poles, reorder = reorder)
if degs is None or np.sum(degs) < 0:
VanP = np.empty((len(x), 0))
else:
VanP = pvP(x, degs, basisp, derIdxs = derIdxs, reorder = reorder,
scl = scl)
return np.block([[VanH, VanP]])
def polyvanderTotal(x:paramList, poles:Np1D, deg : int = None,
basis : str = "MONOMIAL_HEAVISIDE",
derIdxs : List[List[List[int]]] = None,
- reorder : List[int] = None, scl : Np1D = None)\
- -> Tuple[Np2D, List[List[int]], List[int]]:
+ reorder : List[int] = 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 = basis.split("_")[0]
VanR = heavisidevander(x, poles, reorder = reorder)
if deg is None or deg < 0:
VanP = np.empty((len(x), 0))
derIdxs, ordIdxs = np.zeros(0, dtype = int), np.zeros(0, dtype = int)
else:
- VanP, derIdxs, ordIdxs = pvTP(x, deg, basisp, derIdxs = derIdxs,
- reorder = reorder, scl = scl)
+ VanP = pvTP(x, deg, basisp, derIdxs = derIdxs,
+ reorder = reorder, scl = scl)
ordIdxsEff = np.concatenate((np.arange(len(VanR)), ordIdxs + len(VanR)))
- return (np.block([[VanR, VanP],
- [VanP.T.conj(), np.zeros(tuple([VanP.shape[1]] * 2))]]),
- derIdxs, ordIdxsEff)
+ return np.block([[VanR, VanP],
+ [VanP.T.conj(), np.zeros(tuple([VanP.shape[1]] * 2))]])
diff --git a/rrompy/utilities/poly_fitting/moving_least_squares/moving_least_squares_interpolator.py b/rrompy/utilities/poly_fitting/moving_least_squares/moving_least_squares_interpolator.py
index bab415d..5dd78d1 100644
--- a/rrompy/utilities/poly_fitting/moving_least_squares/moving_least_squares_interpolator.py
+++ b/rrompy/utilities/poly_fitting/moving_least_squares/moving_least_squares_interpolator.py
@@ -1,144 +1,144 @@
# Copyright (C) 2018 by the RROMPy authors
#
# This file is part of RROMPy.
#
# RROMPy is free software: you can redistribute it and/or modify
# it under the terms of the GNU Lesser General Public License as published by
# the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
#
# RROMPy is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU Lesser General Public License for more details.
#
# You should have received a copy of the GNU Lesser General Public License
# along with RROMPy. If not, see .
#
import numpy as np
from copy import deepcopy as copy
from rrompy.utilities.base.types import (List, ListAny, DictAny, Np1D, Np2D,
paramList)
from rrompy.utilities.numerical import customPInv, dot
from .vander import mlsweights
from rrompy.utilities.poly_fitting.interpolator import GenericInterpolator
from rrompy.utilities.poly_fitting.polynomial.vander import (polyvander as pv,
polyvanderTotal as pvT)
from rrompy.utilities.exception_manager import RROMPyException, RROMPyAssert
from rrompy.parameter import checkParameterList
__all__ = ['MovingLeastSquaresInterpolator']
class MovingLeastSquaresInterpolator(GenericInterpolator):
"""HERE"""
def __init__(self, other = None):
if other is None: return
self.support = other.support
self.localProjector = other.localProjector
self.localVanders = other.localVanders
self.suppValues = other.suppValues
self.directionalWeights = other.directionalWeights
self.degree = other.degree
self.npar = other.npar
self.radialbasis = other.radialbasis
self.polybasis = other.polybasis
self.evalParams = other.evalParams
self.totalDegree = other.totalDegree
@property
def shape(self):
sh = self.suppValues.shape[1 :] if self.suppValues.ndim > 1 else 1
return sh
@property
def deg(self):
return self.degree
def __call__(self, mu:paramList, der : List[int] = None,
scl : Np1D = None):
if der is not None and np.sum(der) > 0:
raise RROMPyException(("Cannot take derivatives of moving least "
"squares function."))
mu = checkParameterList(mu, self.npar)[0]
sh = self.shape
if sh == 1: sh = tuple([])
values = np.empty((len(mu),) + sh, dtype = np.complex)
for i, m in enumerate(mu):
weights = mlsweights(m, self.support, self.radialbasis,
directionalWeights = self.directionalWeights,
nNearestNeighbor = self.evalParams["nNearestNeighbor"])
weights /= np.linalg.norm(weights)
vanderLS = np.sum(self.localVanders * weights, axis = 2)
RHSLS = dot(self.localProjector * weights, self.suppValues)
if self.totalDegree:
- vanderEval, _, _ = pvT(m, self.deg[0], self.polybasis,
- **self.evalParams)
+ vanderEval = pvT(m, self.deg[0], self.polybasis,
+ **self.evalParams)
else:
vanderEval = pv(m, self.deg, self.polybasis, **self.evalParams)
vanderEval = vanderEval.flatten()
values[i] = dot(vanderEval, dot(customPInv(vanderLS), RHSLS))
return values
def __copy__(self):
return MovingLeastSquaresInterpolator(self)
def __deepcopy__(self, memo):
other = MovingLeastSquaresInterpolator()
(other.support, other.localProjector, other.localVanders,
other.suppValues, other.directionalWeights, other.degree, other.npar,
other.radialbasis, other.polybasis, other.evalParams,
other.totalDegree) = copy(
(self.support, self.localProjector, self.localVanders,
self.suppValues, self.directionalWeights, self.degree,
self.npar, self.radialbasis, self.polybasis,
self.evalParams, self.totalDegree), memo)
return other
def postmultiplyTensorize(self, A:Np2D):
RROMPyAssert(A.shape[0], self.shape[-1], "Shape of output")
self.suppValues = self.suppValues.dot(A)
def pad(self, nleft : List[int] = None, nright : List[int] = None):
if nleft is None: nleft = [0] * len(self.shape)
if nright is None: nright = [0] * len(self.shape)
if not hasattr(nleft, "__len__"): nleft = [nleft]
if not hasattr(nright, "__len__"): nright = [nright]
RROMPyAssert(len(self.shape), len(nleft), "Shape of output")
RROMPyAssert(len(self.shape), len(nright), "Shape of output")
padwidth = [(0, 0)] + [(l, r) for l, r in zip(nleft, nright)]
self.suppValues = np.pad(self.suppValues, padwidth, "constant",
constant_values = (0., 0.))
def setupByInterpolation(self, support:paramList, values:ListAny,
deg:int, polybasis : str = "MONOMIAL_GAUSSIAN",
directionalWeights : Np1D = None,
totalDegree : bool = True,
vanderCoeffs : DictAny = {}):
support = checkParameterList(support)[0]
self.support = copy(support)
if "reorder" in vanderCoeffs.keys():
self.support = self.support[vanderCoeffs["reorder"]]
if "nNearestNeighbor" not in vanderCoeffs.keys():
vanderCoeffs["nNearestNeighbor"] = -1
self.npar = support.shape[1]
if directionalWeights is None:
directionalWeights = np.ones(self.npar)
self.directionalWeights = directionalWeights
self.polybasis, self.radialbasis, _ = polybasis.split("_")
self.totalDegree = totalDegree
self.evalParams = vanderCoeffs
if totalDegree:
- vander, _, _ = pvT(support, deg, self.polybasis, **vanderCoeffs)
+ vander = pvT(support, deg, self.polybasis, **vanderCoeffs)
if not hasattr(deg, "__len__"): deg = [deg] * self.npar
else:
if not hasattr(deg, "__len__"): deg = [deg] * self.npar
vander = pv(support, deg, self.polybasis, **vanderCoeffs)
self.degree = deg
self.localProjector = vander.T.conj()
self.localVanders = np.array([np.outer(van, van.conj()) \
for van in vander])
self.localVanders = np.swapaxes(self.localVanders, 0, 2)
self.suppValues = np.array(values)
diff --git a/rrompy/utilities/poly_fitting/moving_least_squares/vander.py b/rrompy/utilities/poly_fitting/moving_least_squares/vander.py
index 04094e1..1ab35c8 100644
--- a/rrompy/utilities/poly_fitting/moving_least_squares/vander.py
+++ b/rrompy/utilities/poly_fitting/moving_least_squares/vander.py
@@ -1,97 +1,96 @@
# Copyright (C) 2018 by the RROMPy authors
#
# This file is part of RROMPy.
#
# RROMPy is free software: you can redistribute it and/or modify
# it under the terms of the GNU Lesser General Public License as published by
# the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
#
# RROMPy is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU Lesser General Public License for more details.
#
# You should have received a copy of the GNU Lesser General Public License
# along with RROMPy. If not, see .
#
import numpy as np
from .kernel import (radialGaussian, thinPlateSpline, multiQuadric,
nearestNeighbor)
from rrompy.utilities.poly_fitting.polynomial.vander import (polyvander as pvP,
polyvanderTotal as pvTP)
from rrompy.utilities.base.types import (Np1D, Np2D, Tuple, List, paramVal,
paramList)
from rrompy.parameter import checkParameter, checkParameterList
from rrompy.utilities.exception_manager import RROMPyException, RROMPyAssert
__all__ = ['mlsweights', 'polyvander', 'polyvanderTotal']
def mlsweights(x:paramVal, xSupp:paramList, basis:str,
reorder : List[int] = None, directionalWeights : Np1D = None,
nNearestNeighbor : int = -1) -> Np2D:
"""Compute moving least squares weight vector."""
x = checkParameter(x)
xSupp = checkParameterList(xSupp)[0]
x = x.data
xSupp = xSupp.data
if directionalWeights is None:
directionalWeights = np.ones(x.shape[1])
elif not hasattr(directionalWeights, "__len__"):
directionalWeights = directionalWeights * np.ones(x.shape[1])
RROMPyAssert(len(directionalWeights), x.shape[1],
"Number of directional weights")
try:
radialkernel = {"GAUSSIAN" : radialGaussian,
"THINPLATE" : thinPlateSpline,
"MULTIQUADRIC" : multiQuadric,
"NEARESTNEIGHBOR" : nearestNeighbor}[basis.upper()]
except:
raise RROMPyException("Radial basis not recognized.")
if reorder is not None: xSupp = xSupp[reorder]
muDiff = (xSupp.data - x) * directionalWeights
r2 = np.sum(np.abs(muDiff) ** 2., axis = 1).reshape(1, -1)
if basis.upper() == "NEARESTNEIGHBOR":
if nNearestNeighbor > 0 and nNearestNeighbor < len(xSupp):
cutoffValue = np.partition(r2, nNearestNeighbor - 1)[0,
nNearestNeighbor - 1]
r2 /= cutoffValue
else:
r2[0, :] = 1. * (nNearestNeighbor == 0)
return radialkernel(r2)[0]
def polyvander(x:paramVal, xSupp:paramList, degs:List[int], basis:str,
derIdxs : List[List[List[int]]] = None,
reorder : List[int] = None, directionalWeights : Np1D = None,
scl : Np1D = None,
nNearestNeighbor : int = -1) -> Tuple[Np2D, Np2D]:
"""
Compute full Hermite-Vandermonde matrix with specified derivative
directions.
"""
basisp, basisr, _ = basis.split("_")
Weights = mlsweights(x, xSupp, basisr, reorder, directionalWeights,
nNearestNeighbor)
VanP = pvP(xSupp, degs, basisp, derIdxs = derIdxs, reorder = reorder,
scl = scl)
RHP = VanP.T.conj() * Weights
return RHP.dot(VanP), RHP
def polyvanderTotal(x:paramList, xSupp:paramList, deg:int, basis:str,
derIdxs : List[List[List[int]]] = None,
reorder : List[int] = None,
directionalWeights : Np1D = None, scl : Np1D = None,
- nNearestNeighbor : int = -1)\
- -> Tuple[Np2D, Np2D, List[List[int]], List[int]]:
+ nNearestNeighbor : int = -1) -> Tuple[Np2D, Np2D]:
"""
Compute full total degree Hermite-Vandermonde matrix with specified
derivative directions.
"""
basisp, basisr, _ = basis.split("_")
Weights = mlsweights(x, xSupp, basisr, reorder, directionalWeights,
nNearestNeighbor)
- VanP, derIdxs, ordIdxs = pvTP(x, deg, basisp, derIdxs = derIdxs,
- reorder = reorder, scl = scl)
+ VanP = pvTP(x, deg, basisp, derIdxs = derIdxs,
+ reorder = reorder, scl = scl)
RHP = VanP.T.conj() * Weights
- return RHP.dot(VanP), RHP, derIdxs, ordIdxs
+ return RHP.dot(VanP), RHP
diff --git a/rrompy/utilities/poly_fitting/polynomial/polynomial_interpolator.py b/rrompy/utilities/poly_fitting/polynomial/polynomial_interpolator.py
index 3b46e64..0eb05ff 100644
--- a/rrompy/utilities/poly_fitting/polynomial/polynomial_interpolator.py
+++ b/rrompy/utilities/poly_fitting/polynomial/polynomial_interpolator.py
@@ -1,131 +1,129 @@
# Copyright (C) 2018 by the RROMPy authors
#
# This file is part of RROMPy.
#
# RROMPy is free software: you can redistribute it and/or modify
# it under the terms of the GNU Lesser General Public License as published by
# the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
#
# RROMPy is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU Lesser General Public License for more details.
#
# You should have received a copy of the GNU Lesser General Public License
# along with RROMPy. If not, see .
#
import numpy as np
from copy import deepcopy as copy
from rrompy.utilities.base.types import (List, ListAny, DictAny, Np1D, Np2D,
paramList)
from rrompy.utilities.base import freepar as fp
from rrompy.utilities.poly_fitting.interpolator import GenericInterpolator
from rrompy.utilities.poly_fitting.custom_fit import customFit
from .base import polyfitname
from .val import polyval
from .roots import polyroots
from .vander import polyvander as pv, polyvanderTotal as pvT
from rrompy.utilities.numerical import degreeTotalToFull, dot
from rrompy.utilities.exception_manager import RROMPyAssert, RROMPyException
from rrompy.parameter import checkParameterList
__all__ = ['PolynomialInterpolator']
class PolynomialInterpolator(GenericInterpolator):
"""HERE"""
def __init__(self, other = None):
if other is None: return
self.coeffs = other.coeffs
self.npar = other.npar
self.polybasis = other.polybasis
@property
def shape(self):
if self.coeffs.ndim > self.npar:
sh = self.coeffs.shape[self.npar :]
else: sh = tuple([1])
return sh
@property
def deg(self):
return [x - 1 for x in self.coeffs.shape[: self.npar]]
def __getitem__(self, key):
return self.coeffs[key]
def __call__(self, mu:paramList, der : List[int] = None,
scl : Np1D = None):
return polyval(mu, self.coeffs, self.polybasis, der, scl)
def __copy__(self):
return PolynomialInterpolator(self)
def __deepcopy__(self, memo):
other = PolynomialInterpolator()
other.coeffs, other.npar, other.polybasis = copy(
(self.coeffs, self.npar, self.polybasis), memo)
return other
def pad(self, nleft : List[int] = None, nright : List[int] = None):
if nleft is None: nleft = [0] * len(self.shape)
if nright is None: nright = [0] * len(self.shape)
if not hasattr(nleft, "__len__"): nleft = [nleft]
if not hasattr(nright, "__len__"): nright = [nright]
RROMPyAssert(len(self.shape), len(nleft), "Shape of output")
RROMPyAssert(len(self.shape), len(nright), "Shape of output")
padwidth = [(0, 0)] * self.npar
padwidth = padwidth + [(l, r) for l, r in zip(nleft, nright)]
self.coeffs = np.pad(self.coeffs, padwidth, "constant",
constant_values = (0., 0.))
def postmultiplyTensorize(self, A:Np2D):
RROMPyAssert(A.shape[0], self.shape[-1], "Shape of output")
self.coeffs = dot(self.coeffs, A)
def setupByInterpolation(self, support:paramList, values:ListAny,
deg:int, polybasis : str = "MONOMIAL",
verbose : bool = True, totalDegree : bool = True,
vanderCoeffs : DictAny = {},
fitCoeffs : DictAny = {}):
support = checkParameterList(support)[0]
self.npar = support.shape[1]
self.polybasis = polybasis
if totalDegree:
- vander, _, reorder = pvT(support, deg, basis = polybasis,
- **vanderCoeffs)
- vander = vander[:, reorder]
+ vander = pvT(support, deg, basis = polybasis, **vanderCoeffs)
else:
if not hasattr(deg, "__len__"): deg = [deg] * self.npar
vander = pv(support, deg, basis = polybasis, **vanderCoeffs)
outDim = values.shape[1:]
values = values.reshape(values.shape[0], -1)
fitOut = customFit(vander, values, full = True, **fitCoeffs)
P = fitOut[0]
if verbose:
msg = ("Fitting {} samples with degree {} through {}... "
"Conditioning of LS system: {:.4e}.").format(
len(vander), deg,
polyfitname(self.polybasis),
fitOut[1][2][0] / fitOut[1][2][-1])
else: msg = None
if totalDegree:
self.coeffs = degreeTotalToFull(tuple([deg + 1] * self.npar)
+ outDim, self.npar, P)
else:
self.coeffs = P.reshape(tuple([d + 1 for d in deg]) + outDim)
return fitOut[1][1] == vander.shape[1], msg
def roots(self, marginalVals : ListAny = [fp]):
RROMPyAssert(self.shape, (1,), "Shape of output")
RROMPyAssert(len(marginalVals), self.npar, "Number of parameters")
try:
rDim = marginalVals.index(fp)
if rDim < len(marginalVals) - 1 and fp in marginalVals[rDim + 1 :]:
raise
except:
raise RROMPyException(("Exactly 1 'freepar' entry in "
"marginalVals must be provided."))
return polyroots(self.coeffs, self.polybasis, marginalVals)
diff --git a/rrompy/utilities/poly_fitting/polynomial/vander.py b/rrompy/utilities/poly_fitting/polynomial/vander.py
index 4f20015..e3e1e2b 100644
--- a/rrompy/utilities/poly_fitting/polynomial/vander.py
+++ b/rrompy/utilities/poly_fitting/polynomial/vander.py
@@ -1,138 +1,137 @@
# 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 .der import polyder
from rrompy.utilities.base.types import Np1D, Np2D, Tuple, List, paramList
from rrompy.utilities.numerical 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)
C_m[idx, :, j] = np.ravel(np.pad(J_der, mode = "constant",
pad_width = om))
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)[0]
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:
raise RROMPyException("Polynomial basis not recognized.")
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)\
- -> Tuple[Np2D, List[List[int]], List[int]]:
+ reorder : List[int] = None, scl : Np1D = None) -> Np2D:
"""
Compute full total degree Hermite-Vandermonde matrix with specified
derivative directions.
"""
x = checkParameterList(x)[0]
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], derIdxs, ordIdxs
+ return VanBase[:, mask][:, ordIdxs]
diff --git a/rrompy/utilities/poly_fitting/radial_basis/radial_basis_interpolator.py b/rrompy/utilities/poly_fitting/radial_basis/radial_basis_interpolator.py
index e1e0023..23134ec 100644
--- a/rrompy/utilities/poly_fitting/radial_basis/radial_basis_interpolator.py
+++ b/rrompy/utilities/poly_fitting/radial_basis/radial_basis_interpolator.py
@@ -1,147 +1,145 @@
# Copyright (C) 2018 by the RROMPy authors
#
# This file is part of RROMPy.
#
# RROMPy is free software: you can redistribute it and/or modify
# it under the terms of the GNU Lesser General Public License as published by
# the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
#
# RROMPy is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU Lesser General Public License for more details.
#
# You should have received a copy of the GNU Lesser General Public License
# along with RROMPy. If not, see .
#
import numpy as np
from copy import deepcopy as copy
from rrompy.utilities.base.types import (List, ListAny, DictAny, Np1D, Np2D,
paramList)
from rrompy.utilities.poly_fitting.interpolator import GenericInterpolator
from rrompy.utilities.poly_fitting.custom_fit import customFit
from .base import polyfitname
from .val import polyval
from .vander import polyvander as pv, polyvanderTotal as pvT
from rrompy.utilities.numerical import degreeTotalToFull, dot
from rrompy.utilities.exception_manager import RROMPyException, RROMPyAssert
from rrompy.parameter import checkParameterList
__all__ = ['RadialBasisInterpolator']
class RadialBasisInterpolator(GenericInterpolator):
"""HERE"""
def __init__(self, other = None):
if other is None: return
self.support = other.support
self.coeffsGlobal = other.coeffsGlobal
self.coeffsLocal = other.coeffsLocal
self.directionalWeights = other.directionalWeights
self.npar = other.npar
self.polybasis = other.polybasis
self.nNearestNeighbor = other.nNearestNeighbor
@property
def shape(self):
sh = self.coeffsLocal.shape[1 :] if self.coeffsLocal.ndim > 1 else 1
return sh
@property
def deg(self):
return [x - 1 for x in self.coeffsGlobal.shape[: self.npar]]
def __call__(self, mu:paramList, der : List[int] = None,
scl : Np1D = None):
if der is not None and np.sum(der) > 0:
raise RROMPyException(("Cannot take derivatives of radial basis "
"function."))
return polyval(mu, self.coeffsGlobal, self.coeffsLocal, self.support,
self.directionalWeights, self.polybasis,
self.nNearestNeighbor)
def __copy__(self):
return RadialBasisInterpolator(self)
def __deepcopy__(self, memo):
other = RadialBasisInterpolator()
(other.support, other.coeffsGlobal, other.coeffsLocal,
other.directionalWeights, other.npar, other.polybasis,
other.nNearestNeighbor) = copy(
(self.support, self.coeffsGlobal, self.coeffsLocal,
self.directionalWeights, self.npar, self.polybasis,
self.nNearestNeighbor), memo)
return other
def postmultiplyTensorize(self, A:Np2D):
RROMPyAssert(A.shape[0], self.shape[-1], "Shape of output")
self.coeffsLocal = dot(self.coeffsLocal, A)
self.coeffsGlobal = dot(self.coeffsGlobal, A)
def pad(self, nleft : List[int] = None, nright : List[int] = None):
if nleft is None: nleft = [0] * len(self.shape)
if nright is None: nright = [0] * len(self.shape)
if not hasattr(nleft, "__len__"): nleft = [nleft]
if not hasattr(nright, "__len__"): nright = [nright]
RROMPyAssert(len(self.shape), len(nleft), "Shape of output")
RROMPyAssert(len(self.shape), len(nright), "Shape of output")
padwidth = [(0, 0)] + [(l, r) for l, r in zip(nleft, nright)]
self.coeffsLocal = np.pad(self.coeffsLocal, padwidth, "constant",
constant_values = (0., 0.))
padwidth = [(0, 0)] * (self.npar - 1) + padwidth
self.coeffsGlobal = np.pad(self.coeffsGlobal, padwidth, "constant",
constant_values = (0., 0.))
def setupByInterpolation(self, support:paramList, values:ListAny,
deg:int, polybasis : str = "MONOMIAL_GAUSSIAN",
directionalWeights : Np1D = None,
verbose : bool = True, totalDegree : bool = True,
vanderCoeffs : DictAny = {},
fitCoeffs : DictAny = {}):
support = checkParameterList(support)[0]
self.support = copy(support)
if "reorder" in vanderCoeffs.keys():
self.support = self.support[vanderCoeffs["reorder"]]
if "nNearestNeighbor" in vanderCoeffs.keys():
self.nNearestNeighbor = vanderCoeffs["nNearestNeighbor"]
else:
self.nNearestNeighbor = -1
self.npar = support.shape[1]
if directionalWeights is None:
directionalWeights = np.ones(self.npar)
self.directionalWeights = directionalWeights
self.polybasis = polybasis
if totalDegree:
- vander, _, reorder = pvT(support, deg, basis = polybasis,
- directionalWeights = self.directionalWeights,
- **vanderCoeffs)
- vander = vander[reorder]
- vander = vander[:, reorder]
+ vander = pvT(support, deg, basis = polybasis,
+ directionalWeights = self.directionalWeights,
+ **vanderCoeffs)
else:
if not hasattr(deg, "__len__"): deg = [deg] * self.npar
vander = pv(support, deg, basis = polybasis,
directionalWeights = self.directionalWeights,
**vanderCoeffs)
outDim = values.shape[1:]
values = values.reshape(values.shape[0], -1)
values = np.pad(values, ((0, len(vander) - len(values)), (0, 0)),
"constant")
fitOut = customFit(vander, values, full = True, **fitCoeffs)
P = fitOut[0][len(support) :]
if verbose:
msg = ("Fitting {}+{} samples with degree {} through {}... "
"Conditioning of LS system: {:.4e}.").format(
len(support), len(vander) - len(support),
deg, polyfitname(self.polybasis),
fitOut[1][2][0] / fitOut[1][2][-1])
else: msg = None
self.coeffsLocal = fitOut[0][: len(support)]
if totalDegree:
self.coeffsGlobal = degreeTotalToFull(tuple([deg + 1] * self.npar)
+ outDim, self.npar, P)
else:
self.coeffsGlobal = P.reshape(tuple([d + 1 for d in deg]) + outDim)
return fitOut[1][1] == vander.shape[1], msg
diff --git a/rrompy/utilities/poly_fitting/radial_basis/vander.py b/rrompy/utilities/poly_fitting/radial_basis/vander.py
index 9738162..a504ac5 100644
--- a/rrompy/utilities/poly_fitting/radial_basis/vander.py
+++ b/rrompy/utilities/poly_fitting/radial_basis/vander.py
@@ -1,111 +1,107 @@
# 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,
nearestNeighbor)
from rrompy.utilities.poly_fitting.polynomial.vander import (polyvander as pvP,
polyvanderTotal as pvTP)
from rrompy.utilities.base.types import Np1D, Np2D, Tuple, List, paramList
from rrompy.parameter import checkParameterList
from rrompy.utilities.exception_manager import RROMPyException, RROMPyAssert
__all__ = ['rbvander', 'polyvander', 'polyvanderTotal']
def rbvander(x:paramList, basis:str, reorder : List[int] = None,
directionalWeights : Np1D = None,
nNearestNeighbor : int = -1) -> Np2D:
"""Compute radial-basis-Vandermonde matrix."""
x = checkParameterList(x)[0]
x_un = x.unique()
nx = len(x)
if len(x_un) < nx:
raise RROMPyException("Sample points must be distinct.")
del x_un
x = x.data
if directionalWeights is None:
directionalWeights = np.ones(x.shape[1])
elif not hasattr(directionalWeights, "__len__"):
directionalWeights = directionalWeights * np.ones(x.shape[1])
RROMPyAssert(len(directionalWeights), x.shape[1],
"Number of directional weights")
try:
radialkernel = {"GAUSSIAN" : radialGaussian,
"THINPLATE" : thinPlateSpline,
"MULTIQUADRIC" : multiQuadric,
"NEARESTNEIGHBOR" : nearestNeighbor}[basis.upper()]
except:
raise RROMPyException("Radial basis not recognized.")
isnearestneighbor = basis.upper() == "NEARESTNEIGHBOR"
Van = np.zeros((nx, nx))
for j in range(nx):
muDiff = (x.data - x[j]) * directionalWeights
r2j = np.sum(np.abs(muDiff) ** 2., axis = 1).reshape(1, -1)
if isnearestneighbor:
if nNearestNeighbor > 0 and nNearestNeighbor < len(x):
cutoffValue = np.partition(r2j, nNearestNeighbor - 1)[0,
nNearestNeighbor - 1]
r2j /= cutoffValue
else:
r2j[0, :] = 1. * (nNearestNeighbor == 0)
Van[j] = radialkernel(r2j)
return Van
def polyvander(x:paramList, degs:List[int], basis:str,
derIdxs : List[List[List[int]]] = None,
reorder : List[int] = None, directionalWeights : Np1D = None,
scl : Np1D = None, nNearestNeighbor : int = -1) -> Np2D:
"""
Compute full Hermite-Vandermonde matrix with specified derivative
directions.
"""
if derIdxs is not None and np.sum(np.sum(derIdxs)) > 0:
raise RROMPyException(("Cannot take derivatives of radial basis "
"function."))
basisp, basisr = basis.split("_")
VanR = rbvander(x, basisr, reorder = reorder,
directionalWeights = directionalWeights,
nNearestNeighbor = nNearestNeighbor)
VanP = pvP(x, degs, basisp, derIdxs = derIdxs, reorder = reorder,
scl = scl)
return np.block([[VanR, VanP],
[VanP.T.conj(), np.zeros(tuple([VanP.shape[1]] * 2))]])
def polyvanderTotal(x:paramList, deg:int, basis:str,
derIdxs : List[List[List[int]]] = None,
reorder : List[int] = None,
directionalWeights : Np1D = None, scl : Np1D = None,
- nNearestNeighbor : int = -1)\
- -> Tuple[Np2D, List[List[int]], List[int]]:
+ nNearestNeighbor : int = -1) -> 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,
nNearestNeighbor = nNearestNeighbor)
- VanP, derIdxs, ordIdxs = pvTP(x, deg, basisp, derIdxs = derIdxs,
- reorder = reorder, scl = scl)
- ordIdxsEff = np.concatenate((np.arange(len(VanR)), ordIdxs + len(VanR)))
- return (np.block([[VanR, VanP],
- [VanP.T.conj(), np.zeros(tuple([VanP.shape[1]] * 2))]]),
- derIdxs, ordIdxsEff)
+ 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))]])