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))]])