Page MenuHomec4science
Files F62208829

heaviside_to_from_rational.py
No OneTemporary
Actions

File Metadata

Created
Sat, May 11, 15:41

heaviside_to_from_rational.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 rrompy.utilities.poly_fitting.polynomial.polynomial_interpolator import (
PolynomialInterpolator, PolynomialInterpolatorNodal)
from rrompy.utilities.poly_fitting.polynomial.vander import polyvander
from rrompy.utilities.poly_fitting.custom_fit import customFit
from rrompy.utilities.base.types import Np1D, Np2D, Tuple, DictAny, interpEng
from rrompy.parameter.parameter_sampling import (RandomSampler as RS,
QuadratureSampler as QS)
from .val import polyval
from rrompy.utilities.exception_manager import RROMPyException, RROMPyAssert
__all__ = ['heaviside2rational', 'rational2heaviside']
def heaviside2rational(c:Np2D, poles:Np1D, murange : Np1D = None,
basis : str = "MONOMIAL_HEAVISIDE", basisD : str = None,
parameterMap : DictAny = 1.) \
-> Tuple[interpEng, interpEng]:
basisN = basis.split("_")[0]
if basisD is None: basisD = basisN
num = PolynomialInterpolator()
den = PolynomialInterpolatorNodal()
den.polybasis, den.nodes = basisD, poles
if murange is None:
multiplier = [1., 1.j]
avgs = [np.mean(np.real(poles)), np.mean(np.imag(poles))]
ranges = np.array([[np.min(np.real(poles)), np.max(np.real(poles))],
[np.min(np.imag(poles)), np.max(np.imag(poles))]])
domIdx = np.argmax(ranges[:, 1] - ranges[:, 0])
murange = [multiplier[domIdx] * x
+ multiplier[1 - domIdx] * avgs[1 - domIdx]
for x in ranges[domIdx, :]]
if basis == "CHEBYSHEV":
sampler = QS(murange, "CHEBYSHEV", parameterMap)
elif basis == "LEGENDRE":
sampler = QS(murange, "GAUSSLEGENDRE", parameterMap)
else:
sampler = RS(murange, "HALTON", parameterMap)
xAux = sampler.generatePoints(len(c))
valsAux = den(xAux) * polyval(xAux, c, poles, basis)
num.setupByInterpolation(xAux, valsAux, len(c) - 1, basisN)
return num, den
def rational2heaviside(num:interpEng, den:interpEng,
murange : Np1D = np.array([-1., 1.]), scl : Np1D = None,
parameterMap : DictAny = 1.) \
-> Tuple[Np2D, Np1D, str]:
if (not isinstance(num, PolynomialInterpolator)
or not isinstance(den, PolynomialInterpolator)):
raise RROMPyException(("Rational numerator and denominator must be "
"polynomial interpolators."))
RROMPyAssert(num.npar, 1, "Number of parameters")
RROMPyAssert(den.npar, 1, "Number of parameters")
basis = num.polybasis + "_HEAVISIDE"
c = np.zeros_like(num.coeffs)
poles = den.roots()
if len(poles) > 0:
Psp = num(poles)
Qsp = den(poles)
Qspder = den(poles, 1, scl)
polesBad = np.abs(Qsp) >= 1e-5
Psp[..., polesBad] = 0.
Qspder[polesBad] = 1.
c[: len(poles)] = (Psp / Qspder).T
if len(c) > len(poles):
from rrompy.parameter.parameter_sampling import (RandomSampler as RS,
QuadratureSampler as QS)
if num.polybasis == "CHEBYSHEV":
sampler = QS(murange, "CHEBYSHEV", parameterMap)
elif num.polybasis == "LEGENDRE":
sampler = QS(murange, "GAUSSLEGENDRE", parameterMap)
else:
sampler = RS(murange, "HALTON", parameterMap)
xAux = sampler.generatePoints(len(c))
valsAux = num(xAux) / den(xAux)
if len(poles) > 0:
valsAux -= polyval(xAux, c, poles, basis)
VanAux = polyvander(xAux, [len(c) - len(poles) - 1], num.polybasis)
c[len(poles) :] = customFit(VanAux, valsAux.T)
if len(poles) > 0: poles[polesBad] = np.inf
return c, poles, basis

Event Timeline

c4science · Help