Page Menu
Home
c4science
Search
Configure Global Search
Log In
Files
F61662789
matrix_3_pod.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
Wed, May 8, 04:21
Size
7 KB
Mime Type
text/x-python
Expires
Fri, May 10, 04:21 (2 d)
Engine
blob
Format
Raw Data
Handle
17542057
Attached To
R6746 RationalROMPy
matrix_3_pod.py
View Options
import
numpy
as
np
import
scipy.sparse
as
sp
from
mpl_toolkits.mplot3d
import
Axes3D
from
matplotlib
import
pyplot
as
plt
from
rrompy.hfengines.base
import
MatrixEngineBase
as
MEB
from
rrompy.reduction_methods.standard
import
RationalInterpolant
as
RI
from
rrompy.reduction_methods.standard
import
ReducedBasis
as
RB
from
rrompy.parameter.parameter_sampling
import
(
QuadratureSampler
as
QS
,
QuadratureSamplerTotal
as
QST
,
ManualSampler
as
MS
,
RandomSampler
as
RS
)
verb
=
50
deg
=
2
size
=
1
show_sample
=
True
show_norm
=
False
Delta
=
0
MN
=
15
R
=
(
MN
+
3
)
*
(
MN
+
2
)
*
(
MN
+
1
)
//
6
S
=
[
int
(
np
.
ceil
(
R
**
(
1.
/
3.
)))]
*
3
PODTol
=
1e-8
samples
=
"centered"
samples
=
"standard"
#samples, nDer = "standard_MMM", 2
algo
=
"rational"
#algo = "RB"
sampling
=
"quadrature"
sampling
=
"quadrature_total"
sampling
=
"random"
if
samples
==
"standard"
:
radial
=
""
# radial = "_gaussian"
# radial = "_thinplate"
# radial = "_multiquadric"
rW0
=
5.
radialWeight
=
[
rW0
]
*
2
assert
Delta
<=
0
if
size
==
1
:
mu0
=
[
0.
]
*
3
mutar
=
[
.
8
,
.
8
,
.
8
]
murange
=
[[
-
1.
]
*
3
,
[
1.
]
*
3
]
aEff
=
1.
#25
bEff
=
1.
-
aEff
murangeEff
=
[[
aEff
*
murange
[
0
][
0
]
+
bEff
*
murange
[
1
][
0
],
aEff
*
murange
[
0
][
1
]
+
bEff
*
murange
[
1
][
1
],
aEff
*
murange
[
0
][
2
]
+
bEff
*
murange
[
1
][
2
]],
[
aEff
*
murange
[
1
][
0
]
+
bEff
*
murange
[
0
][
0
],
aEff
*
murange
[
1
][
1
]
+
bEff
*
murange
[
0
][
1
],
aEff
*
murange
[
1
][
2
]
+
bEff
*
murange
[
0
][
2
]]]
N
=
150
exp
=
1.05
assert
exp
>
1.
empty
=
sp
.
csr_matrix
((
np
.
zeros
(
0
),
np
.
zeros
(
0
),
np
.
zeros
(
N
+
1
)),
shape
=
(
N
,
N
))
solver
=
MEB
(
verbosity
=
verb
)
solver
.
npar
=
3
A0
=
sp
.
spdiags
([
np
.
arange
(
1
,
1
+
2
*
N
,
2
)
**
exp
-
(
2
*
(
N
//
4
))
**
exp
],
[
0
],
N
,
N
)
if
deg
==
1
:
solver
.
nAs
=
4
solver
.
As
=
[
A0
,
sp
.
eye
(
N
),
-
sp
.
eye
(
N
),
-
sp
.
eye
(
N
)]
elif
deg
==
2
:
solver
.
nAs
=
7
solver
.
As
=
[
A0
]
+
[
empty
]
+
[
-
sp
.
eye
(
N
)]
+
[
empty
]
*
3
+
[
sp
.
eye
(
N
)]
elif
deg
==
3
:
solver
.
nAs
=
13
solver
.
As
=
[
A0
]
+
[
empty
]
+
[
-
sp
.
eye
(
N
)]
+
[
empty
]
*
9
+
[
sp
.
eye
(
N
)]
np
.
random
.
seed
(
420
)
solver
.
nbs
=
1
solver
.
bs
=
[
np
.
random
.
randn
(
N
)
+
1.j
*
np
.
random
.
randn
(
N
)]
if
algo
==
"rational"
:
params
=
{
'N'
:
MN
,
'M'
:
MN
+
Delta
,
'S'
:
S
,
'POD'
:
True
}
if
samples
==
"standard"
:
params
[
'polybasis'
]
=
"CHEBYSHEV"
+
radial
# params['polybasis'] = "LEGENDRE" + radial
# params['polybasis'] = "MONOMIAL" + radial
params
[
'radialDirectionalWeights'
]
=
radialWeight
elif
samples
==
"centered"
:
params
[
'polybasis'
]
=
"MONOMIAL"
params
[
'S'
]
=
R
else
:
#MMM
params
[
'S'
]
=
nDer
*
int
(
np
.
ceil
(
R
/
nDer
))
method
=
RI
else
:
#if algo == "RB":
params
=
{
'R'
:(
MN
+
3
+
Delta
)
*
(
MN
+
2
+
Delta
)
*
(
MN
+
1
+
Delta
)
//
6
,
'S'
:
S
,
'POD'
:
True
,
'PODTolerance'
:
PODTol
}
if
samples
==
"standard"
:
pass
elif
samples
==
"centered"
:
params
[
'S'
]
=
R
else
:
#MMM
params
[
'S'
]
=
nDer
*
int
(
np
.
ceil
(
R
/
nDer
))
method
=
RB
if
samples
==
"standard"
:
if
sampling
==
"quadrature"
:
params
[
'sampler'
]
=
QS
(
murange
,
"CHEBYSHEV"
)
# params['sampler'] = QS(murange, "GAUSSLEGENDRE")
# params['sampler'] = QS(murange, "UNIFORM")
params
[
'S'
]
=
[
max
(
j
,
MN
+
1
)
for
j
in
params
[
'S'
]]
elif
sampling
==
"quadrature_total"
:
params
[
'sampler'
]
=
QST
(
murange
,
"CHEBYSHEV"
)
params
[
'S'
]
=
R
else
:
# if sampling == "random":
params
[
'sampler'
]
=
RS
(
murange
,
"HALTON"
)
params
[
'S'
]
=
R
elif
samples
==
"centered"
:
params
[
'sampler'
]
=
MS
(
murange
,
points
=
[
mu0
])
elif
samples
==
"standard_MMM"
:
if
sampling
==
"quadrature"
:
SdirEff
=
int
(
np
.
ceil
(
int
(
np
.
ceil
(
R
/
nDer
))
**
(
1.
/
3
)))
pts
=
QS
(
murange
,
"CHEBYSHEV"
)
.
generatePoints
([
SdirEff
]
*
3
)
elif
sampling
==
"quadrature_total"
:
pts
=
QST
(
murange
,
"CHEBYSHEV"
)
.
generatePoints
(
int
(
np
.
ceil
(
R
/
nDer
)))
else
:
# if sampling == "random":
pts
=
RS
(
murange
,
"HALTON"
)
.
generatePoints
(
int
(
np
.
ceil
(
R
/
nDer
)))
params
[
'sampler'
]
=
MS
(
murange
,
points
=
pts
)
approx
=
method
(
solver
,
mu0
=
mu0
,
approxParameters
=
params
,
verbosity
=
verb
)
if
samples
==
"standard"
:
approx
.
samplingEngine
.
allowRepeatedSamples
=
False
approx
.
setupApprox
()
if
show_sample
:
approx
.
plotApprox
(
mutar
,
name
=
'u_app'
,
what
=
"REAL"
)
approx
.
plotHF
(
mutar
,
name
=
'u_HF'
,
what
=
"REAL"
)
approx
.
plotErr
(
mutar
,
name
=
'err'
,
what
=
"REAL"
)
# approx.plotRes(mutar, name = 'res', what = "REAL")
appErr
=
approx
.
normErr
(
mutar
)
solNorm
=
approx
.
normHF
(
mutar
)
resNorm
=
approx
.
normRes
(
mutar
)
RHSNorm
=
approx
.
normRHS
(
mutar
)
print
((
'SolNorm:
\t
{}
\n
Err:
\t
{}
\n
ErrRel:
\t
{}'
)
.
format
(
solNorm
,
appErr
,
np
.
divide
(
appErr
,
solNorm
)))
print
((
'RHSNorm:
\t
{}
\n
Res:
\t
{}
\n
ResRel:
\t
{}'
)
.
format
(
RHSNorm
,
resNorm
,
np
.
divide
(
resNorm
,
RHSNorm
)))
fig
=
plt
.
figure
(
figsize
=
(
8
,
6
))
ax
=
Axes3D
(
fig
)
ax
.
scatter
(
approx
.
trainedModel
.
data
.
mus
(
0
),
approx
.
trainedModel
.
data
.
mus
(
1
),
approx
.
trainedModel
.
data
.
mus
(
2
),
'.'
)
ax
.
set_xlim3d
(
murangeEff
[
0
][
0
],
murangeEff
[
1
][
0
])
ax
.
set_ylim3d
(
murangeEff
[
0
][
1
],
murangeEff
[
1
][
1
])
ax
.
set_zlim3d
(
murangeEff
[
0
][
2
],
murangeEff
[
1
][
2
])
plt
.
show
()
plt
.
close
()
approx
.
verbosity
=
0
approx
.
trainedModel
.
verbosity
=
0
if
algo
==
"rational"
and
approx
.
N
>
0
:
from
plot_zero_set_3
import
plotZeroSet3
# muZeroVals, Qvals = plotZeroSet3(murange, murangeEff, approx, mu0, 200,
# [None, mu0[1], mu0[2]], [1., 1., 1.])
# muZeroVals, Qvals = plotZeroSet3(murange, murangeEff, approx, mu0, 200,
# [None, None, mu0[2]], [1., 1., 1.])
# muZeroVals, Qvals = plotZeroSet3(murange, murangeEff, approx, mu0, 200,
# [None, mu0[1], None], [1., 1., 1.])
plotZeroSet3
(
murange
,
murangeEff
,
approx
,
mu0
,
25
,
[
None
,
None
,
None
],
[
1.
,
1.
,
1.
],
imagTol
=
1e-2
)
if
show_norm
:
solver
.
_solveBatchSize
=
25
from
plot_inf_set_3
import
plotInfSet3
muInfVals
,
normEx
,
normApp
,
normRes
,
normErr
,
beta
=
plotInfSet3
(
murange
,
murangeEff
,
approx
,
mu0
,
25
,
[
None
,
mu0
[
1
],
mu0
[
2
]],
[
1.
,
1.
,
1.
],
relative
=
False
,
normalizeDen
=
True
)
muInfVals
,
normEx
,
normApp
,
normRes
,
normErr
,
beta
=
plotInfSet3
(
murange
,
murangeEff
,
approx
,
mu0
,
25
,
[
None
,
None
,
mu0
[
2
]],
[
1.
,
1.
,
1.
],
relative
=
False
,
normalizeDen
=
True
)
muInfVals
,
normEx
,
normApp
,
normRes
,
normErr
,
beta
=
plotInfSet3
(
murange
,
murangeEff
,
approx
,
mu0
,
25
,
[
None
,
mu0
[
1
],
None
],
[
1.
,
1.
,
1.
],
relative
=
False
,
normalizeDen
=
True
)
Event Timeline
Log In to Comment