Homec4science
Diffusion RationalROMPy (feature_greedy_coarsening) Browse

RROMPy/rrompy/utilities/poly_fitting/fit_utils.py
c9551b918a85feature_greedy_coarseni…

Last Change
Raw File

fit_utils.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 numpy import pi, polynomial as po
from scipy.special import binom
from rrompy.utilities.base.types import Np1D, Np2D
__all__ = ['polybases', 'polyval', 'polyder', 'polyvalder', 'polyvander',
'polyfitname', 'polyroots', 'polydomcoeff']
polybases = ["CHEBYSHEV", "LEGENDRE", "MONOMIAL"]
polyval = {"CHEBYSHEV" : po.chebyshev.chebval, "LEGENDRE" : po.legendre.legval,
"MONOMIAL" : po.polynomial.polyval}
polyder = {"CHEBYSHEV" : po.chebyshev.chebder, "LEGENDRE" : po.legendre.legder,
"MONOMIAL" : po.polynomial.polyder}
polyvalder = {
"CHEBYSHEV" : lambda x, c, m = 1, scl = 1.:
po.chebyshev.chebval(x, polyder["CHEBYSHEV"](c, m, scl)),
"LEGENDRE" : lambda x, c, m = 1, scl = 1.:
po.legendre.legval(x, polyder["LEGENDRE"](c, m, scl)),
"MONOMIAL" : lambda x, c, m = 1, scl = 1.:
po.polynomial.polyval(x, polyder["MONOMIAL"](c, m, scl))}
polyvander = {
"CHEBYSHEV" : lambda x, deg, scl = 1.:
polyvanderConfluence(po.chebyshev.chebvander, polyder["CHEBYSHEV"],
x, deg, scl),
"LEGENDRE" : lambda x, deg, scl = 1:
polyvanderConfluence(po.legendre.legvander, polyder["LEGENDRE"],
x, deg, scl),
"MONOMIAL" : lambda x, deg, scl = 1:
polyvanderConfluence(po.polynomial.polyvander, polyder["MONOMIAL"],
x, deg, scl)}
polyfitname = {"CHEBYSHEV" : "chebfit", "LEGENDRE" : "legfit",
"MONOMIAL" : "polyfit"}
polyroots = {"CHEBYSHEV" : po.chebyshev.chebroots,
"LEGENDRE" : po.legendre.legroots,
"MONOMIAL" : po.polynomial.polyroots}
polydomcoeff = {"CHEBYSHEV" : lambda n: 2. ** (n - 1) if n > 0 else 1.,
"LEGENDRE" : lambda n: (2. ** n * (pi * n) ** -.5 if n > 10
else .5 ** n * binom(2 * n, n)),
"MONOMIAL" : lambda n: 1.}
def polyvanderConfluence(vander:callable, derivative:callable, x:Np1D, deg:int,
scl : float = 1.) -> Np2D:
"""Compute Vandermonde matrix even in case of confluence."""
x_un, idx_un, cnt_un = np.unique(x, return_inverse = True,
return_counts = True)
Van = vander(x, deg)
der_max = np.max(cnt_un) - 1
if der_max > 0:
C_der = np.zeros((deg + 1, deg + 1), dtype = float)
for j in range(deg + 1):
ej = np.zeros(deg + 1)
ej[j] = 1.
j_der = derivative(ej, 1, scl)
C_der[: len(j_der), j] = j_der
for der in range(1, der_max + 1):
# remove first occurrence of each node
for i_un in np.nonzero(cnt_un > der - 1)[0]:
idx_un[np.nonzero(idx_un == i_un)[0][0]] = -1
idx_loc = np.nonzero(idx_un > -1)[0]
Van[idx_loc, :] = Van[idx_loc, :].dot(C_der[:, :]) / der
return Van
c4science · Help