Page MenuHomec4science
Files F61388208

polynomial_interpolator.py
No OneTemporary
Actions

File Metadata

Created
Mon, May 6, 08:46

polynomial_interpolator.py
View Options

# Copyright (C) 2018-2020 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 <http://www.gnu.org/licenses/>.
#
from copy import deepcopy as copy
import numpy as np
from scipy.special import factorial as fact
from collections.abc import Iterable
from itertools import combinations
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
from .polynomial_algebra import changePolyBasis, polyTimes
from rrompy.utilities.numerical import dot, baseDistanceMatrix
from rrompy.utilities.numerical.degree import degreeTotalToFull
from rrompy.utilities.exception_manager import RROMPyAssert, RROMPyException
from rrompy.parameter import checkParameterList
__all__ = ['PolynomialInterpolator', 'PolynomialInterpolatorNodal']
class PolynomialInterpolator(GenericInterpolator):
"""Function class with setup by polynomial interpolation."""
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):
if hasattr(self, "_dirPivot"):
mu = checkParameterList(mu)(self._dirPivot)
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 isinstance(nleft, Iterable): nleft = [nleft]
if not isinstance(nright, Iterable): 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)
self.npar = support.shape[1]
self.polybasis = polybasis
if not totalDegree and not isinstance(deg, Iterable):
deg = [deg] * self.npar
vander = pv(support, deg, basis = polybasis, **vanderCoeffs)
RROMPyAssert(len(vander), len(values), "Number of support values")
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")
rDim = marginalVals.index(fp)
if rDim < len(marginalVals) - 1 and fp in marginalVals[rDim + 1 :]:
raise RROMPyException(("Exactly 1 'freepar' entry in "
"marginalVals must be provided."))
return polyroots(self.coeffs, self.polybasis, marginalVals)
class PolynomialInterpolatorNodal(PolynomialInterpolator):
"""
Function class with setup by polynomial interpolation. Stores roots of
monomial polynomial instead of coefficients. Only for 1D.
"""
def __init__(self, other = None):
self.npar = 1
if other is None: return
self.nodes = other.nodes
self.polybasis = other.polybasis
@property
def nodes(self):
return self._nodes
@nodes.setter
def nodes(self, nodes):
self.coeffs = None
self._nodes = nodes
@property
def coeffs(self):
if self._coeffs is None: self.buildCoeffs()
return self._coeffs
@coeffs.setter
def coeffs(self, coeffs):
self._coeffs = coeffs
@property
def shape(self):
return (1,)
@property
def deg(self):
return [len(self.nodes)]
def __getitem__(self, key):
return self.coeffs[key]
def __call__(self, mu:paramList, der : List[int] = None,
scl : Np1D = None):
dirPivot = self._dirPivot if hasattr(self, "_dirPivot") else 0
if der is None: der = 0
elif isinstance(der, (list,tuple,np.ndarray,)): der = der[dirPivot]
if scl is None: scl = 1.
elif isinstance(scl, (list,tuple,np.ndarray,)): scl = scl[dirPivot]
mu = checkParameterList(mu)(dirPivot)
val = np.zeros(len(mu), dtype = np.complex)
if der == self.deg[0]:
val[:] = 1.
elif der >= 0 and der < self.deg[0]:
plDist = baseDistanceMatrix(mu, self.nodes, magnitude = False)[...,
0]
for terms in combinations(np.arange(self.deg[0]),
self.deg[0] - der):
val += np.prod(plDist[:, list(terms)], axis = 1)
return scl ** der * fact(der) * val
def __copy__(self):
return PolynomialInterpolatorNodal(self)
def __deepcopy__(self, memo):
other = PolynomialInterpolatorNodal()
other.nodes, other.polybasis = copy((self.nodes, self.polybasis), memo)
return other
def buildCoeffs(self):
local = [np.array([- pl, 1.], dtype = np.complex) for pl in self.nodes]
N = len(local)
while N > 1:
for j in range(N // 2):
local[j] = polyTimes(local[j], local[- 1 - j])
local = local[(N - 1) // 2 :: -1]
N = len(local)
self._coeffs = changePolyBasis(local[0], None, "MONOMIAL",
self.polybasis)
def pad(self, *args, **kwargs):
raise RROMPyException(("Padding not allowed for polynomials in nodal "
"form"))
def postmultiplyTensorize(self, *args, **kwargs):
raise RROMPyException(("Post-multiply not allowed for polynomials in "
"nodal form"))
def setupByInterpolation(self, support:paramList, *args, **kwargs):
support = checkParameterList(support)
self.npar = support.shape[1]
if self.npar > 1:
raise RROMPyException(("Polynomial in nodal form must have "
"scalar output"))
output = super().setupByInterpolation(support, *args, **kwargs)
self._nodes = super().roots()
return output
def roots(self, marginalVals : ListAny = [fp]):
return np.array(self.nodes)

Event Timeline

c4science · Help