Page Menu
Home
c4science
Search
Configure Global Search
Log In
Files
F84071541
radial_fitting.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
Fri, Sep 20, 14:03
Size
6 KB
Mime Type
text/x-python
Expires
Sun, Sep 22, 14:03 (2 d)
Engine
blob
Format
Raw Data
Handle
20819777
Attached To
R6746 RationalROMPy
radial_fitting.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
import
customFit
from
rrompy.utilities.poly_fitting.radial_basis
import
(
radialGaussian
,
thinPlateSpline
,
multiQuadric
,
polybases
,
polyfitname
,
polydomcoeff
,
polyval
,
polyvander
,
polyvanderTotal
)
from
rrompy.utilities.numerical.degree
import
degreeTotalToFull
from
rrompy.parameter
import
checkParameterList
def
test_monomial_gaussian
():
polyrbname
=
"MONOMIAL_GAUSSIAN"
assert
polyrbname
in
polybases
fitname
=
polyfitname
(
polyrbname
)
domcoeff
=
polydomcoeff
(
5
,
polyrbname
)
assert
fitname
==
"polyfit_gaussian"
assert
np
.
isclose
(
domcoeff
,
1.
,
rtol
=
1e-5
)
directionalWeights
=
np
.
array
([
5.
])
xSupp
=
checkParameterList
(
np
.
arange
(
-
1
,
3
),
1
)[
0
]
cRBCoeffs
=
np
.
array
([
-
1.
,
3.
,
-
3.
,
1.
,
1.
,
2.
,
-.
5
])
globalCoeffs
=
cRBCoeffs
[
4
:]
localCoeffs
=
cRBCoeffs
[:
4
]
ySupp
=
1
+
2.
*
xSupp
.
data
-
.
5
*
xSupp
.
data
**
2.
xx
=
np
.
linspace
(
-
2.
,
3.
,
100
)
yy
=
polyval
(
checkParameterList
(
xx
,
1
)[
0
],
globalCoeffs
,
localCoeffs
,
xSupp
,
directionalWeights
,
polyrbname
)
yyman
=
1.
+
2.
*
xx
-
.
5
*
xx
**
2.
for
j
,
xc
in
enumerate
(
np
.
arange
(
-
1
,
3
)):
r2j
=
(
5.
*
(
xx
-
xc
))
**
2.
rbj
=
radialGaussian
(
r2j
)
assert
np
.
allclose
(
rbj
,
np
.
exp
(
-.
5
*
r2j
))
yyman
+=
localCoeffs
[
j
]
*
rbj
ySupp
+=
localCoeffs
[
j
]
*
radialGaussian
((
directionalWeights
[
0
]
*
(
xSupp
.
data
-
xc
))
**
2.
)
assert
np
.
allclose
(
yy
,
yyman
,
atol
=
1e-5
)
VanT
=
polyvander
(
xSupp
,
[
2
],
polyrbname
,
directionalWeights
=
directionalWeights
)
ySupp
=
np
.
pad
(
ySupp
.
flatten
(),
(
0
,
len
(
VanT
)
-
len
(
xSupp
)),
"constant"
)
out
=
customFit
(
VanT
,
ySupp
)
assert
np
.
allclose
(
out
,
cRBCoeffs
,
atol
=
1e-8
)
def
test_legendre_thinplate
():
polyrbname
=
"LEGENDRE_THINPLATE"
assert
polyrbname
in
polybases
fitname
=
polyfitname
(
polyrbname
)
domcoeff
=
polydomcoeff
(
5
,
polyrbname
)
assert
fitname
==
"legfit_thinplate"
assert
np
.
isclose
(
domcoeff
,
63.
/
8
,
rtol
=
1e-5
)
directionalWeights
=
np
.
array
([
.
5
])
xSupp
=
checkParameterList
(
np
.
arange
(
-
1
,
3
),
1
)[
0
]
cRBCoeffs
=
np
.
array
([
-
1.
,
3.
,
-
3.
,
1.
,
1.
,
2.
,
-.
5
])
localCoeffs
=
cRBCoeffs
[:
4
]
globalCoeffs
=
cRBCoeffs
[
4
:]
ySupp
=
1
+
2.
*
xSupp
.
data
-
.
5
*
(
.
5
*
(
3.
*
xSupp
.
data
**
2.
-
1.
))
xx
=
np
.
linspace
(
-
2.
,
3.
,
100
)
yy
=
polyval
(
checkParameterList
(
xx
,
1
)[
0
],
globalCoeffs
,
localCoeffs
,
xSupp
,
directionalWeights
,
polyrbname
)
yyman
=
1.
+
2.
*
xx
-
.
5
*
(
.
5
*
(
3.
*
xx
**
2.
-
1.
))
for
j
,
xc
in
enumerate
(
np
.
arange
(
-
1
,
3
)):
r2j
=
(
directionalWeights
[
0
]
*
(
xx
-
xc
))
**
2.
rbj
=
thinPlateSpline
(
r2j
)
assert
np
.
allclose
(
rbj
,
.
5
*
r2j
*
np
.
log
(
np
.
finfo
(
float
)
.
eps
+
r2j
))
yyman
+=
localCoeffs
[
j
]
*
rbj
ySupp
+=
localCoeffs
[
j
]
*
thinPlateSpline
((
directionalWeights
[
0
]
*
(
xSupp
.
data
-
xc
))
**
2.
)
assert
np
.
allclose
(
yy
,
yyman
,
atol
=
1e-5
)
VanT
=
polyvander
(
xSupp
,
[
2
],
polyrbname
,
directionalWeights
=
directionalWeights
)
ySupp
=
np
.
pad
(
ySupp
.
flatten
(),
(
0
,
len
(
VanT
)
-
len
(
xSupp
)),
"constant"
)
out
=
customFit
(
VanT
,
ySupp
)
assert
np
.
allclose
(
out
,
cRBCoeffs
,
atol
=
1e-8
)
def
test_chebyshev_multiquadric
():
polyrbname
=
"CHEBYSHEV_MULTIQUADRIC"
assert
polyrbname
in
polybases
fitname
=
polyfitname
(
polyrbname
)
domcoeff
=
polydomcoeff
(
5
,
polyrbname
)
assert
fitname
==
"chebfit_multiquadric"
assert
np
.
isclose
(
domcoeff
,
16
,
rtol
=
1e-5
)
directionalWeights
=
np
.
array
([
1.
])
xSupp
=
checkParameterList
(
np
.
arange
(
-
1
,
3
),
1
)[
0
]
cRBCoeffs
=
np
.
array
([
-
1.
,
3.
,
-
3.
,
1.
,
1.
,
2.
,
-.
5
])
localCoeffs
=
cRBCoeffs
[:
4
]
globalCoeffs
=
cRBCoeffs
[
4
:]
ySupp
=
1
+
2.
*
xSupp
.
data
-
.
5
*
(
2.
*
xSupp
.
data
**
2.
-
1.
)
xx
=
np
.
linspace
(
-
2.
,
3.
,
100
)
yy
=
polyval
(
checkParameterList
(
xx
,
1
)[
0
],
globalCoeffs
,
localCoeffs
,
xSupp
,
directionalWeights
,
polyrbname
)
yyman
=
1.
+
2.
*
xx
-
.
5
*
(
2.
*
xx
**
2.
-
1.
)
for
j
,
xc
in
enumerate
(
np
.
arange
(
-
1
,
3
)):
r2j
=
(
directionalWeights
[
0
]
*
(
xx
-
xc
))
**
2.
rbj
=
multiQuadric
(
r2j
)
assert
np
.
allclose
(
rbj
,
np
.
power
(
r2j
+
1
,
-.
5
))
yyman
+=
localCoeffs
[
j
]
*
rbj
ySupp
+=
localCoeffs
[
j
]
*
multiQuadric
((
directionalWeights
[
0
]
*
(
xSupp
.
data
-
xc
))
**
2.
)
assert
np
.
allclose
(
yy
,
yyman
,
atol
=
1e-5
)
VanT
=
polyvander
(
xSupp
,
[
2
],
polyrbname
,
directionalWeights
=
directionalWeights
)
ySupp
=
np
.
pad
(
ySupp
.
flatten
(),
(
0
,
len
(
VanT
)
-
len
(
xSupp
)),
"constant"
)
out
=
customFit
(
VanT
,
ySupp
)
assert
np
.
allclose
(
out
,
cRBCoeffs
,
atol
=
1e-8
)
def
test_total_degree_2d
():
values
=
lambda
x
,
y
:
(
x
-
3.
)
**
2.
*
y
-
(
x
+
1.
)
*
y
**
2.
polyrbname
=
"CHEBYSHEV_GAUSSIAN"
xs
,
ys
=
np
.
meshgrid
(
np
.
linspace
(
0.
,
4.
,
5
),
np
.
linspace
(
0.
,
4.
,
4
))
xySupp
=
np
.
concatenate
((
xs
.
reshape
(
-
1
,
1
),
ys
.
reshape
(
-
1
,
1
)),
axis
=
1
)
zs
=
values
(
xs
,
ys
)
zSupp
=
zs
.
flatten
()
deg
=
3
directionalWeights
=
[
2.
,
1.
]
VanT
=
polyvanderTotal
(
xySupp
,
deg
,
polyrbname
,
directionalWeights
=
directionalWeights
)
cFit
=
np
.
linalg
.
solve
(
VanT
,
np
.
pad
(
zSupp
,
(
0
,
len
(
VanT
)
-
len
(
zSupp
)),
'constant'
))
globCoeff
=
degreeTotalToFull
([
deg
+
1
]
*
2
,
2
,
cFit
[
len
(
zSupp
)
:])
localCoeffs
=
cFit
[:
len
(
zSupp
)]
globalCoeffs
=
globCoeff
xx
,
yy
=
np
.
meshgrid
(
np
.
linspace
(
0.
,
4.
,
100
),
np
.
linspace
(
0.
,
4.
,
100
))
xxyy
=
np
.
concatenate
((
xx
.
reshape
(
-
1
,
1
),
yy
.
reshape
(
-
1
,
1
)),
axis
=
1
)
zz
=
polyval
(
xxyy
,
globalCoeffs
,
localCoeffs
,
xySupp
,
directionalWeights
,
polyrbname
)
.
reshape
(
xx
.
shape
)
zzex
=
values
(
xx
,
yy
)
error
=
np
.
abs
(
zz
-
zzex
)
print
(
np
.
max
(
error
))
assert
np
.
max
(
error
)
<
1e-10
Event Timeline
Log In to Comment