Page Menu
Home
c4science
Search
Configure Global Search
Log In
Files
F75819038
heaviside_to_from_rational.py
No One
Temporary
Actions
Download File
Edit File
Delete File
View Transforms
Subscribe
Mute Notifications
Award Token
Subscribers
None
File Metadata
Details
File Info
Storage
Attached
Created
Sun, Aug 4, 11:50
Size
4 KB
Mime Type
text/x-python
Expires
Tue, Aug 6, 11:50 (2 d)
Engine
blob
Format
Raw Data
Handle
19615257
Attached To
R6746 RationalROMPy
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
Event Timeline
Log In to Comment