Homec4science
Diffusion RationalROMPy (master) Browse

RROMPy/rrompy/utilities/poly_fitting/heaviside/heaviside_to_from_rational.py
6c4eee04e0dfmaster

Last Change
Raw File

heaviside_to_from_rational.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 rrompy.utilities.poly_fitting.polynomial.polynomial_interpolator import \
PolynomialInterpolator
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]:
num = PolynomialInterpolator()
den = PolynomialInterpolator()
basisN = basis.split("_")[0]
if basisD is None: basisD = basisN
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, :]]
extraPt = None
while extraPt is None or np.any(np.isclose(extraPt, poles)):
extraPt = murange[0] + (murange[1] - murange[0]) * np.random.rand(1)
denAuxPts = np.concatenate((poles, extraPt))
denAuxVals = np.concatenate((np.zeros(len(poles)), [1.]))
den.setupByInterpolation(denAuxPts, denAuxVals, len(poles), basisD)
den.coeffs /= np.linalg.norm(den.coeffs)
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()
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)
- polyval(xAux, c, poles, basis)).T
VanAux = polyvander(xAux, [len(c) - len(poles) - 1], num.polybasis)
c[len(poles) :] = customFit(VanAux, valsAux)
poles[polesBad] = np.inf
return c, poles, basis
c4science · Help