Homec4science
Diffusion RationalROMPy (develop) Browse

RROMPy/rrompy/utilities/poly_fitting/polynomial/polynomial_algebra.py
9253ed2b70a5develop

Last Change
Raw File

polynomial_algebra.py
View Options

# 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 <http://www.gnu.org/licenses/>.
#
import numpy as np
from copy import deepcopy as copy
from rrompy.utilities.base.types import Np1D, Np2D, Tuple, List, interpEng
from .vander import polyvander
from .polynomial_interpolator import PolynomialInterpolator
from rrompy.utilities.numerical import (multifactorial, customPInv,
hashDerivativeToIdx as hashD,
hashIdxToDerivative as hashI)
from rrompy.utilities.exception_manager import RROMPyException
__all__ = ['changePolyBasis', 'polyTimes', 'polyDivide', 'polyTimesTable']
def changePolyBasis(P:Np2D, dim : int = None, basis0 : str = "MONOMIAL",
basisF : str = "MONOMIAL") -> Np2D:
if basis0 == basisF: return P
if dim is None: dim = P.ndim
if basis0 != "MONOMIAL" and basisF != "MONOMIAL":
return changePolyBasis(changePolyBasis(P, dim, basis0, "MONOMIAL"),
dim, "MONOMIAL", basisF)
basisD = basisF if basis0 == "MONOMIAL" else basis0
R = copy(P)
N = np.max(P.shape[: dim]) - 1
vander = polyvander([0], N, basisD, [list(range(N + 1))])
if basis0 == "MONOMIAL": vander = customPInv(vander)
for j in range(dim):
R = np.tensordot(vander, R, (-1, j))
return R
def polyTimes(P:Np2D, Q:Np2D, dim : int = None, Pbasis : str = "MONOMIAL",
Qbasis : str = "MONOMIAL", Rbasis : str = "MONOMIAL") -> Np2D:
if not isinstance(P, (np.ndarray,)): P = np.array(P)
if not isinstance(Q, (np.ndarray,)): Q = np.array(Q)
P = changePolyBasis(P, dim, Pbasis, "MONOMIAL")
Q = changePolyBasis(Q, dim, Qbasis, "MONOMIAL")
if dim is None: dim = P.ndim
if dim <= 0: return
R = np.zeros([x + y - 1 for (x, y) in zip(P.shape[: dim], Q.shape[: dim])],
dtype = P.dtype)
if dim == 1:
for j, Qj in enumerate(Q):
R[j : j + len(P)] = R[j : j + len(P)] + Qj * P
else:
for j, Qj in enumerate(Q):
for l, Pl in enumerate(P):
R[j + l] = R[j + l] + polyTimes(Pl, Qj, dim - 1)
return changePolyBasis(R, dim, "MONOMIAL", Rbasis)
def polyDivide(P:Np2D, Q:Np2D, dim : int = None, Pbasis : str = "MONOMIAL",
Qbasis : str = "MONOMIAL",
Rbasis : str = "MONOMIAL") -> Tuple[Np2D, Np2D]:
if not isinstance(P, (np.ndarray,)): P = np.array(P)
if not isinstance(Q, (np.ndarray,)): Q = np.array(Q)
P = changePolyBasis(P, dim, Pbasis, "MONOMIAL")
Pc = copy(P)
Q = changePolyBasis(Q, dim, Qbasis, "MONOMIAL")
if dim is None: dim = P.ndim
if dim <= 0: return
R = np.zeros([x - y + 1 for (x, y) in zip(P.shape[: dim], Q.shape[: dim])],
dtype = P.dtype)
if dim == 1:
for i in range(len(R) - 1, -1, -1):
try:
R[i] = Pc[-1] / Q[-1]
except:
raise RROMPyException(("Numerical instability in polynomial "
"quotient."))
Pc = Pc[: -1]
for j, Qj in enumerate(Q[::-1]):
if j > 0: Pc[-j] = Pc[-j] - R[i] * Qj
else:
raise RROMPyException(("Quotient of multivariate polynomials not "
"supported."))
return (changePolyBasis(R, dim, "MONOMIAL", Rbasis),
changePolyBasis(Pc, dim, "MONOMIAL", Rbasis))
def polyTimesTable(P:interpEng, mus:Np1D, reorder:List[int],
derIdxs:List[List[List[int]]], scl : Np1D = None) -> Np2D:
if not isinstance(P, PolynomialInterpolator):
raise RROMPyException(("Polynomial to evaluate must be a polynomial "
"interpolator."))
S = len(reorder)
T = np.zeros((S, S), dtype = np.complex)
idxGlob = 0
for j, derIdx in enumerate(derIdxs):
nder = len(derIdx)
idxLoc = np.arange(S)[np.logical_and(reorder >= idxGlob,
reorder < idxGlob + nder)]
idxGlob += nder
Pval = []
for der in range(nder):
derI = hashI(der, P.npar)
Pval += [P(mus[j], derI, scl) / multifactorial(derI)]
for derI, derIdxI in enumerate(derIdx):
for derJ, derIdxJ in enumerate(derIdx):
diffIdx = [x - y for (x, y) in zip(derIdxI, derIdxJ)]
if all([x >= 0 for x in diffIdx]):
diffj = hashD(diffIdx)
T[idxLoc[derI], idxLoc[derJ]] = Pval[diffj]
return T
c4science · Help