Page Menu
Home
c4science
Search
Configure Global Search
Log In
Files
F68344876
matrix_greedy.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, Jun 26, 23:53
Size
3 KB
Mime Type
text/x-python
Expires
Fri, Jun 28, 23:53 (2 d)
Engine
blob
Format
Raw Data
Handle
18567752
Attached To
R6746 RationalROMPy
matrix_greedy.py
View Options
import
numpy
as
np
import
scipy.sparse
as
sp
from
matplotlib
import
pyplot
as
plt
from
rrompy.hfengines.base
import
MatrixEngineBase
as
MEB
from
rrompy.reduction_methods.distributed_greedy
import
\
RationalInterpolantGreedy
as
Pade
from
rrompy.reduction_methods.distributed_greedy
import
\
RBDistributedGreedy
as
RB
test
=
1
timed
=
False
method
=
"Pade"
method
=
"RB"
verb
=
200
errorEstimatorKind
=
"BARE"
#errorEstimatorKind = "BASIC"
#errorEstimatorKind = "EXACT"
N
=
100
solver
=
MEB
(
verbosity
=
verb
)
solver
.
npar
=
1
solver
.
nAs
=
2
if
test
==
1
:
solver
.
As
=
[
sp
.
spdiags
([
np
.
arange
(
1
,
1
+
N
)],
[
0
],
N
,
N
),
-
sp
.
eye
(
N
)]
elif
test
==
2
:
solver
.
setSolver
(
"SOLVE"
)
fftB
=
np
.
fft
.
fft
(
np
.
eye
(
N
))
*
N
**-.
5
solver
.
As
=
[
fftB
.
dot
(
np
.
multiply
(
np
.
arange
(
1
,
1
+
N
),
fftB
.
conj
())
.
T
),
-
np
.
eye
(
N
)]
np
.
random
.
seed
(
420
)
solver
.
nbs
=
1
solver
.
bs
=
[
np
.
random
.
randn
(
N
)
+
1.j
*
np
.
random
.
randn
(
N
)]
mu0
=
10.25
murange
=
[
1.25
,
19.25
]
mutars
=
np
.
linspace
(
murange
[
0
],
murange
[
1
],
500
)
if
method
==
"Pade"
:
params
=
{
'muBounds'
:
murange
,
'nTestPoints'
:
200
,
'Delta'
:
0
,
'S'
:
5
,
'greedyTol'
:
1e-2
,
'polybasis'
:
"CHEBYSHEV"
,
'errorEstimatorKind'
:
errorEstimatorKind
}
approx
=
Pade
(
solver
,
mu0
=
mu0
,
approxParameters
=
params
,
verbosity
=
verb
)
elif
method
==
"RB"
:
params
=
{
'muBounds'
:
murange
,
'nTestPoints'
:
500
,
'greedyTol'
:
1e-2
,
'S'
:
5
}
approx
=
RB
(
solver
,
mu0
=
mu0
,
approxParameters
=
params
,
verbosity
=
verb
)
if
timed
:
from
time
import
clock
start_time
=
clock
()
approx
.
greedy
()
print
(
"Time: "
,
clock
()
-
start_time
)
else
:
approx
.
greedy
(
True
)
approx
.
samplingEngine
.
verbosity
=
0
approx
.
trainedModel
.
verbosity
=
0
approx
.
verbosity
=
0
normApp
=
np
.
zeros
(
len
(
mutars
))
norm
=
np
.
zeros_like
(
normApp
)
res
=
np
.
zeros_like
(
normApp
)
err
=
np
.
zeros_like
(
normApp
)
for
j
in
range
(
len
(
mutars
)):
normApp
[
j
]
=
approx
.
normApprox
(
mutars
[
j
])
norm
[
j
]
=
approx
.
normHF
(
mutars
[
j
])
res
[
j
]
=
(
approx
.
estimatorNormEngine
.
norm
(
approx
.
getRes
(
mutars
[
j
]))
/
approx
.
estimatorNormEngine
.
norm
(
approx
.
getRHS
(
mutars
[
j
])))
err
[
j
]
=
approx
.
normErr
(
mutars
[
j
])
/
approx
.
normHF
(
mutars
[
j
])
resApp
=
approx
.
errorEstimator
(
mutars
)
plt
.
figure
()
plt
.
semilogy
(
mutars
,
norm
)
plt
.
semilogy
(
mutars
,
normApp
,
'--'
)
plt
.
semilogy
(
np
.
real
(
approx
.
mus
(
0
)),
1.25
*
np
.
max
(
norm
)
*
np
.
ones
(
approx
.
mus
.
shape
,
dtype
=
float
),
'rx'
)
plt
.
xlim
(
murange
)
plt
.
grid
()
plt
.
show
()
plt
.
close
()
plt
.
figure
()
plt
.
semilogy
(
mutars
,
res
)
plt
.
semilogy
(
mutars
,
resApp
,
'--'
)
plt
.
semilogy
(
np
.
real
(
approx
.
mus
(
0
)),
4.
*
np
.
max
(
resApp
)
*
np
.
ones
(
approx
.
mus
.
shape
,
dtype
=
float
),
'rx'
)
plt
.
xlim
(
murange
)
plt
.
grid
()
plt
.
show
()
plt
.
close
()
plt
.
figure
()
plt
.
semilogy
(
mutars
,
err
)
plt
.
xlim
(
murange
)
plt
.
grid
()
plt
.
show
()
plt
.
close
()
polesTrue
=
np
.
arange
(
1
,
1
+
N
)
polesTrue
=
polesTrue
[
polesTrue
>=
murange
[
0
]]
polesTrue
=
polesTrue
[
polesTrue
<=
murange
[
1
]]
polesApp
=
approx
.
getPoles
()
mask
=
(
np
.
real
(
polesApp
)
<
murange
[
0
])
|
(
np
.
real
(
polesApp
)
>
murange
[
1
])
print
(
"Outliers:"
,
polesApp
[
mask
])
polesApp
=
polesApp
[
~
mask
]
plt
.
figure
()
plt
.
plot
(
np
.
real
(
polesApp
),
np
.
imag
(
polesApp
),
'kx'
)
plt
.
plot
(
polesTrue
,
np
.
zeros_like
(
polesTrue
),
'm.'
)
plt
.
axis
(
'equal'
)
plt
.
grid
()
plt
.
show
()
plt
.
close
()
Event Timeline
Log In to Comment