Page Menu
Home
c4science
Search
Configure Global Search
Log In
Files
F81415923
polynomial_algebra.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
Thu, Sep 5, 23:12
Size
5 KB
Mime Type
text/x-python
Expires
Sat, Sep 7, 23:12 (2 d)
Engine
blob
Format
Raw Data
Handle
20563266
Attached To
R6746 RationalROMPy
polynomial_algebra.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
copy
import
deepcopy
as
copy
from
rrompy.utilities.base.types
import
Np1D
,
Np2D
,
Tuple
,
List
,
interpEng
from
.vander
import
polyvander
from
.polynomial_interpolator
import
PolynomialInterpolator
from
rrompy.utilities.numerical
import
(
multifactorial
,
customPInv
,
hashDerivativeToIdx
as
hashD
,
hashIdxToDerivative
as
hashI
)
from
rrompy.utilities.exception_manager
import
RROMPyException
__all__
=
[
'changePolyBasis'
,
'polyTimes'
,
'polyDivide'
,
'polyTimesTable'
]
def
changePolyBasis
(
P
:
Np2D
,
dim
:
int
=
None
,
basis0
:
str
=
"MONOMIAL"
,
basisF
:
str
=
"MONOMIAL"
)
->
Np2D
:
if
basis0
==
basisF
:
return
P
if
dim
is
None
:
dim
=
P
.
ndim
if
basis0
!=
"MONOMIAL"
and
basisF
!=
"MONOMIAL"
:
return
changePolyBasis
(
changePolyBasis
(
P
,
dim
,
basis0
,
"MONOMIAL"
),
dim
,
"MONOMIAL"
,
basisF
)
basisD
=
basisF
if
basis0
==
"MONOMIAL"
else
basis0
R
=
copy
(
P
)
N
=
np
.
max
(
P
.
shape
[:
dim
])
-
1
vander
=
polyvander
([
0
],
N
,
basisD
,
[
list
(
range
(
N
+
1
))])
if
basis0
==
"MONOMIAL"
:
vander
=
customPInv
(
vander
)
for
j
in
range
(
dim
):
R
=
np
.
tensordot
(
vander
,
R
,
(
-
1
,
j
))
return
R
def
polyTimes
(
P
:
Np2D
,
Q
:
Np2D
,
dim
:
int
=
None
,
Pbasis
:
str
=
"MONOMIAL"
,
Qbasis
:
str
=
"MONOMIAL"
,
Rbasis
:
str
=
"MONOMIAL"
)
->
Np2D
:
if
not
isinstance
(
P
,
(
np
.
ndarray
,)):
P
=
np
.
array
(
P
)
if
not
isinstance
(
Q
,
(
np
.
ndarray
,)):
Q
=
np
.
array
(
Q
)
P
=
changePolyBasis
(
P
,
dim
,
Pbasis
,
"MONOMIAL"
)
Q
=
changePolyBasis
(
Q
,
dim
,
Qbasis
,
"MONOMIAL"
)
if
dim
is
None
:
dim
=
P
.
ndim
if
dim
<=
0
:
return
R
=
np
.
zeros
([
x
+
y
-
1
for
(
x
,
y
)
in
zip
(
P
.
shape
[:
dim
],
Q
.
shape
[:
dim
])],
dtype
=
P
.
dtype
)
if
dim
==
1
:
for
j
,
Qj
in
enumerate
(
Q
):
R
[
j
:
j
+
len
(
P
)]
=
R
[
j
:
j
+
len
(
P
)]
+
Qj
*
P
else
:
for
j
,
Qj
in
enumerate
(
Q
):
for
l
,
Pl
in
enumerate
(
P
):
R
[
j
+
l
]
=
R
[
j
+
l
]
+
polyTimes
(
Pl
,
Qj
,
dim
-
1
)
return
changePolyBasis
(
R
,
dim
,
"MONOMIAL"
,
Rbasis
)
def
polyDivide
(
P
:
Np2D
,
Q
:
Np2D
,
dim
:
int
=
None
,
Pbasis
:
str
=
"MONOMIAL"
,
Qbasis
:
str
=
"MONOMIAL"
,
Rbasis
:
str
=
"MONOMIAL"
)
->
Tuple
[
Np2D
,
Np2D
]:
if
not
isinstance
(
P
,
(
np
.
ndarray
,)):
P
=
np
.
array
(
P
)
if
not
isinstance
(
Q
,
(
np
.
ndarray
,)):
Q
=
np
.
array
(
Q
)
P
=
changePolyBasis
(
P
,
dim
,
Pbasis
,
"MONOMIAL"
)
Pc
=
copy
(
P
)
Q
=
changePolyBasis
(
Q
,
dim
,
Qbasis
,
"MONOMIAL"
)
if
dim
is
None
:
dim
=
P
.
ndim
if
dim
<=
0
:
return
R
=
np
.
zeros
([
x
-
y
+
1
for
(
x
,
y
)
in
zip
(
P
.
shape
[:
dim
],
Q
.
shape
[:
dim
])],
dtype
=
P
.
dtype
)
if
dim
==
1
:
for
i
in
range
(
len
(
R
)
-
1
,
-
1
,
-
1
):
try
:
R
[
i
]
=
Pc
[
-
1
]
/
Q
[
-
1
]
except
:
raise
RROMPyException
((
"Numerical instability in polynomial "
"quotient."
))
Pc
=
Pc
[:
-
1
]
for
j
,
Qj
in
enumerate
(
Q
[::
-
1
]):
if
j
>
0
:
Pc
[
-
j
]
=
Pc
[
-
j
]
-
R
[
i
]
*
Qj
else
:
raise
RROMPyException
((
"Quotient of multivariate polynomials not "
"supported."
))
return
(
changePolyBasis
(
R
,
dim
,
"MONOMIAL"
,
Rbasis
),
changePolyBasis
(
Pc
,
dim
,
"MONOMIAL"
,
Rbasis
))
def
polyTimesTable
(
P
:
interpEng
,
mus
:
Np1D
,
reorder
:
List
[
int
],
derIdxs
:
List
[
List
[
List
[
int
]]],
scl
:
Np1D
=
None
)
->
Np2D
:
if
not
isinstance
(
P
,
PolynomialInterpolator
):
raise
RROMPyException
((
"Polynomial to evaluate must be a polynomial "
"interpolator."
))
S
=
len
(
reorder
)
T
=
np
.
zeros
((
S
,
S
),
dtype
=
np
.
complex
)
idxGlob
=
0
for
j
,
derIdx
in
enumerate
(
derIdxs
):
nder
=
len
(
derIdx
)
idxLoc
=
np
.
arange
(
S
)[
np
.
logical_and
(
reorder
>=
idxGlob
,
reorder
<
idxGlob
+
nder
)]
idxGlob
+=
nder
Pval
=
[]
for
der
in
range
(
nder
):
derI
=
hashI
(
der
,
P
.
npar
)
Pval
+=
[
P
(
mus
[
j
],
derI
,
scl
)
/
multifactorial
(
derI
)]
for
derI
,
derIdxI
in
enumerate
(
derIdx
):
for
derJ
,
derIdxJ
in
enumerate
(
derIdx
):
diffIdx
=
[
x
-
y
for
(
x
,
y
)
in
zip
(
derIdxI
,
derIdxJ
)]
if
all
([
x
>=
0
for
x
in
diffIdx
]):
diffj
=
hashD
(
diffIdx
)
T
[
idxLoc
[
derI
],
idxLoc
[
derJ
]]
=
Pval
[
diffj
]
return
T
Event Timeline
Log In to Comment