Page MenuHomec4science
Files F60694453

vander.py
No OneTemporary
Actions

File Metadata

Created
Thu, May 2, 00:35

vander.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/>.
#
import numpy as np
from collections.abc import Iterable
from .base import polybases
from .der import polyder
from rrompy.utilities.base.types import Np1D, Np2D, List, paramList
from rrompy.utilities.numerical.degree import degreeSet, degreeSetMixed
from rrompy.parameter import checkParameterList
from rrompy.utilities.exception_manager import RROMPyException, RROMPyAssert
__all__ = ['polyvander']
def firstDerTransition(npar:int, TDirac:List[Np2D], basis:str,
scl : Np1D = None) -> Np2D:
"""Manage step from function samples to function derivatives."""
C_m = np.zeros((npar, len(TDirac), len(TDirac)), dtype = float)
for j, Tj in enumerate(TDirac):
m, om = [0] * npar, [(0, 0)] * npar
for idx in range(npar):
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 polyvander(x:paramList, degs:List[int], basis:str,
derIdxs : List[List[List[int]]] = None,
reorder : List[int] = None, scl : Np1D = None,
degreetype : str = "TENSOR", mixedIdx : int = -1) -> 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
)
"""
x = checkParameterList(x)
npar = x.shape[1]
degreetype = degreetype.upper().strip().replace(" ","")
if degreetype not in ["TENSOR", "TOTAL", "MIXED"]:
raise RROMPyException("Degree type not recognized.")
if degreetype == "MIXED" and npar == 1: degreetype = "TOTAL"
if not isinstance(degs, (list, tuple, np.ndarray,)):
degreetype, degs = "TOTAL", [degs] * npar
if degreetype == "TOTAL":
if np.any(np.array(degs) != degs[0]):
raise RROMPyException(("Cannot force total degree if prescribed "
"degrees are different"))
if degreetype == "MIXED":
if mixedIdx == -1:
raise RROMPyException("Must specify mixedIdx for mixed vander.")
RROMPyAssert(len(degs), npar, "Number of parameters")
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
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, npar):
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 isinstance(derFlat[j], Iterable):
derFlat[j] = [derFlat[j]]
RROMPyAssert(len(derFlat[j]), npar, "Number of dimensions")
#manage mixed derivatives
TDirac = [y.reshape([d + 1 for d in degs])
for y in np.eye(VanBase.shape[-1], dtype = int)]
Cs_loc = firstDerTransition(npar, 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(npar):
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, :]
if degreetype != "TENSOR":
if degreetype == "TOTAL":
derIdxs, mask = degreeSet(degs[0], npar, degreetype, True)
degM = degs[0]
else: #if degreetype == "MIXED":
degS = degs[0] if mixedIdx else degs[1]
derIdxs, mask = degreeSetMixed(degs[mixedIdx], degS, npar,
mixedIdx, degreetype, True)
degM = max(degs)
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(degs[0] + 1):
idxLocal = rangeAux[derTotal == j][::-1]
idxPrev += len(idxLocal)
ordIdxs[idxPrev - len(idxLocal) : idxPrev] = idxLocal
Van = Van[:, mask][:, ordIdxs]
return Van

Event Timeline

c4science · Help