Page Menu
Home
c4science
Search
Configure Global Search
Log In
Files
F85133684
generic_greedy_approximant.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 27, 00:24
Size
30 KB
Mime Type
text/x-python
Expires
Sun, Sep 29, 00:24 (1 d, 23 h)
Engine
blob
Format
Raw Data
Handle
21133177
Attached To
R6746 RationalROMPy
generic_greedy_approximant.py
View Options
# Copyright (C) 2018-2020 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/>.
#
from
abc
import
abstractmethod
from
copy
import
deepcopy
as
copy
import
numpy
as
np
from
matplotlib
import
pyplot
as
plt
from
rrompy.hfengines.base.linear_affine_engine
import
checkIfAffine
from
rrompy.reduction_methods.standard.generic_standard_approximant
import
(
GenericStandardApproximant
)
from
rrompy.utilities.base.types
import
(
Np1D
,
Np2D
,
Tuple
,
List
,
paramVal
,
paramList
,
sampList
)
from
rrompy.utilities.base
import
verbosityManager
as
vbMng
from
rrompy.utilities.numerical
import
dot
from
rrompy.utilities.expression
import
expressionEvaluator
from
rrompy.utilities.exception_manager
import
(
RROMPyException
,
RROMPyAssert
,
RROMPyWarning
)
from
rrompy.sampling.sample_list
import
sampleList
from
rrompy.parameter
import
emptyParameterList
,
parameterList
from
rrompy.utilities.parallel
import
masterCore
__all__
=
[
'GenericGreedyApproximant'
]
def
localL2Distance
(
mus
:
Np2D
,
badmus
:
Np2D
)
->
Np2D
:
return
np
.
linalg
.
norm
(
np
.
tile
(
mus
[
...
,
np
.
newaxis
],
[
1
,
1
,
len
(
badmus
)])
-
badmus
[
...
,
np
.
newaxis
]
.
T
,
axis
=
1
)
def
pruneSamples
(
mus
:
paramList
,
badmus
:
paramList
,
tol
:
float
=
1e-8
)
->
Np1D
:
"""Remove from mus all the elements which are too close to badmus."""
if
isinstance
(
mus
,
(
parameterList
,
sampleList
)):
mus
=
mus
.
data
if
isinstance
(
badmus
,
(
parameterList
,
sampleList
)):
badmus
=
badmus
.
data
if
len
(
badmus
)
==
0
:
return
np
.
arange
(
len
(
mus
))
proximity
=
np
.
min
(
localL2Distance
(
mus
,
badmus
),
axis
=
1
)
return
np
.
where
(
proximity
<=
tol
)[
0
]
class
GenericGreedyApproximant
(
GenericStandardApproximant
):
"""
ROM greedy interpolant computation for parametric problems
(ABSTRACT).
Args:
HFEngine: HF problem solver.
mu0(optional): Default parameter. Defaults to 0.
approxParameters(optional): Dictionary containing values for main
parameters of approximant. Recognized keys are:
- 'POD': kind of snapshots orthogonalization; allowed values
include 0, 1/2, and 1; defaults to 1, i.e. POD;
- 'scaleFactorDer': scaling factors for derivative computation;
defaults to 'AUTO';
- 'S': number of starting training points;
- 'sampler': sample point generator;
- 'greedyTol': uniform error tolerance for greedy algorithm;
defaults to 1e-2;
- 'collinearityTol': collinearity tolerance for greedy algorithm;
defaults to 0.;
- 'maxIter': maximum number of greedy steps; defaults to 1e2;
- 'nTestPoints': number of test points; defaults to 5e2;
- 'samplerTrainSet': training sample points generator; defaults to
sampler.
Defaults to empty dict.
verbosity(optional): Verbosity level. Defaults to 10.
Attributes:
HFEngine: HF problem solver.
mu0: Default parameter.
mus: Array of snapshot parameters.
approxParameters: Dictionary containing values for main parameters of
approximant. Recognized keys are in parameterList.
parameterListSoft: Recognized keys of soft approximant parameters:
- 'POD': kind of snapshots orthogonalization;
- 'scaleFactorDer': scaling factors for derivative computation;
- 'greedyTol': uniform error tolerance for greedy algorithm;
- 'collinearityTol': collinearity tolerance for greedy algorithm;
- 'maxIter': maximum number of greedy steps;
- 'nTestPoints': number of test points;
- 'samplerTrainSet': training sample points generator.
parameterListCritical: Recognized keys of critical approximant
parameters:
- 'S': total number of samples current approximant relies upon;
- 'sampler': sample point generator.
verbosity: Verbosity level.
POD: Kind of snapshots orthogonalization.
scaleFactorDer: Scaling factors for derivative computation.
S: number of test points.
sampler: Sample point generator.
greedyTol: Uniform error tolerance for greedy algorithm.
collinearityTol: Collinearity tolerance for greedy algorithm.
maxIter: maximum number of greedy steps.
nTestPoints: number of starting training points.
samplerTrainSet: training sample points generator.
muBounds: list of bounds for parameter values.
samplingEngine: Sampling engine.
uHF: High fidelity solution(s) with parameter(s) lastSolvedHF as
sampleList.
lastSolvedHF: Parameter(s) corresponding to last computed high fidelity
solution(s) as parameterList.
uApproxReduced: Reduced approximate solution(s) with parameter(s)
lastSolvedApprox as sampleList.
lastSolvedApproxReduced: Parameter(s) corresponding to last computed
reduced approximate solution(s) as parameterList.
uApprox: Approximate solution(s) with parameter(s) lastSolvedApprox as
sampleList.
lastSolvedApprox: Parameter(s) corresponding to last computed
approximate solution(s) as parameterList.
"""
def
__init__
(
self
,
*
args
,
**
kwargs
):
self
.
_preInit
()
if
not
hasattr
(
self
,
"_affine_lvl"
):
self
.
_affine_lvl
=
[]
self
.
_affine_lvl
+=
[
1
]
self
.
_addParametersToList
([
"greedyTol"
,
"collinearityTol"
,
"maxIter"
,
"nTestPoints"
,
"samplerTrainSet"
],
[
1e-2
,
0.
,
1e2
,
5e2
,
"AUTO"
])
super
()
.
__init__
(
*
args
,
**
kwargs
)
self
.
_postInit
()
@property
def
greedyTol
(
self
):
"""Value of greedyTol."""
return
self
.
_greedyTol
@greedyTol.setter
def
greedyTol
(
self
,
greedyTol
):
if
greedyTol
<
0
:
raise
RROMPyException
(
"greedyTol must be non-negative."
)
if
hasattr
(
self
,
"_greedyTol"
)
and
self
.
greedyTol
is
not
None
:
greedyTolold
=
self
.
greedyTol
else
:
greedyTolold
=
-
1
self
.
_greedyTol
=
greedyTol
self
.
_approxParameters
[
"greedyTol"
]
=
self
.
greedyTol
if
greedyTolold
!=
self
.
greedyTol
:
self
.
resetSamples
()
@property
def
collinearityTol
(
self
):
"""Value of collinearityTol."""
return
self
.
_collinearityTol
@collinearityTol.setter
def
collinearityTol
(
self
,
collinearityTol
):
if
collinearityTol
<
0
:
raise
RROMPyException
(
"collinearityTol must be non-negative."
)
if
(
hasattr
(
self
,
"_collinearityTol"
)
and
self
.
collinearityTol
is
not
None
):
collinearityTolold
=
self
.
collinearityTol
else
:
collinearityTolold
=
-
1
self
.
_collinearityTol
=
collinearityTol
self
.
_approxParameters
[
"collinearityTol"
]
=
self
.
collinearityTol
if
collinearityTolold
!=
self
.
collinearityTol
:
self
.
resetSamples
()
@property
def
maxIter
(
self
):
"""Value of maxIter."""
return
self
.
_maxIter
@maxIter.setter
def
maxIter
(
self
,
maxIter
):
if
maxIter
<=
0
:
raise
RROMPyException
(
"maxIter must be positive."
)
if
hasattr
(
self
,
"_maxIter"
)
and
self
.
maxIter
is
not
None
:
maxIterold
=
self
.
maxIter
else
:
maxIterold
=
-
1
self
.
_maxIter
=
maxIter
self
.
_approxParameters
[
"maxIter"
]
=
self
.
maxIter
if
maxIterold
!=
self
.
maxIter
:
self
.
resetSamples
()
@property
def
nTestPoints
(
self
):
"""Value of nTestPoints."""
return
self
.
_nTestPoints
@nTestPoints.setter
def
nTestPoints
(
self
,
nTestPoints
):
if
nTestPoints
<=
0
:
raise
RROMPyException
(
"nTestPoints must be positive."
)
if
not
np
.
isclose
(
nTestPoints
,
np
.
int
(
nTestPoints
)):
raise
RROMPyException
(
"nTestPoints must be an integer."
)
nTestPoints
=
np
.
int
(
nTestPoints
)
if
hasattr
(
self
,
"_nTestPoints"
)
and
self
.
nTestPoints
is
not
None
:
nTestPointsold
=
self
.
nTestPoints
else
:
nTestPointsold
=
-
1
self
.
_nTestPoints
=
nTestPoints
self
.
_approxParameters
[
"nTestPoints"
]
=
self
.
nTestPoints
if
nTestPointsold
!=
self
.
nTestPoints
:
self
.
resetSamples
()
@property
def
samplerTrainSet
(
self
):
"""Value of samplerTrainSet."""
return
self
.
_samplerTrainSet
@samplerTrainSet.setter
def
samplerTrainSet
(
self
,
samplerTrainSet
):
if
(
isinstance
(
samplerTrainSet
,
(
str
,))
and
samplerTrainSet
.
upper
()
==
"AUTO"
):
samplerTrainSet
=
self
.
sampler
if
'generatePoints'
not
in
dir
(
samplerTrainSet
):
raise
RROMPyException
(
"samplerTrainSet type not recognized."
)
if
(
hasattr
(
self
,
'_samplerTrainSet'
)
and
self
.
samplerTrainSet
not
in
[
None
,
"AUTO"
]):
samplerTrainSetOld
=
self
.
samplerTrainSet
self
.
_samplerTrainSet
=
samplerTrainSet
self
.
_approxParameters
[
"samplerTrainSet"
]
=
self
.
samplerTrainSet
if
(
not
'samplerTrainSetOld'
in
locals
()
or
samplerTrainSetOld
!=
self
.
samplerTrainSet
):
self
.
resetSamples
()
def
resetSamples
(
self
):
"""Reset samples."""
super
()
.
resetSamples
()
self
.
_mus
=
emptyParameterList
()
def
_affineResidualMatricesContraction
(
self
,
rb
:
Np2D
,
rA
:
Np2D
=
None
)
\
->
Tuple
[
Np1D
,
Np1D
,
Np1D
]:
self
.
assembleReducedResidualBlocks
(
full
=
rA
is
not
None
)
# 'ij,jk,ik->k', resbb, radiusb, radiusb.conj()
ff
=
np
.
sum
(
self
.
trainedModel
.
data
.
resbb
.
dot
(
rb
)
*
rb
.
conj
(),
axis
=
0
)
if
rA
is
None
:
return
ff
# 'ijk,jkl,il->l', resAb, radiusA, radiusb.conj()
Lf
=
np
.
sum
(
np
.
tensordot
(
self
.
trainedModel
.
data
.
resAb
,
rA
,
2
)
*
rb
.
conj
(),
axis
=
0
)
# 'ijkl,klt,ijt->t', resAA, radiusA, radiusA.conj()
LL
=
np
.
sum
(
np
.
tensordot
(
self
.
trainedModel
.
data
.
resAA
,
rA
,
2
)
*
rA
.
conj
(),
axis
=
(
0
,
1
))
return
ff
,
Lf
,
LL
def
getErrorEstimatorAffine
(
self
,
mus
:
Np1D
)
->
Np1D
:
"""Standard residual estimator."""
checkIfAffine
(
self
.
HFEngine
,
"apply affinity-based error estimator"
,
False
,
self
.
_affine_lvl
)
self
.
HFEngine
.
buildA
()
self
.
HFEngine
.
buildb
()
mus
=
self
.
checkParameterList
(
mus
)
tMverb
,
self
.
trainedModel
.
verbosity
=
self
.
trainedModel
.
verbosity
,
0
uApproxRs
=
self
.
getApproxReduced
(
mus
)
.
data
self
.
trainedModel
.
verbosity
=
tMverb
muTestEff
=
self
.
mapParameterList
(
mus
)
radiusA
=
np
.
empty
((
len
(
self
.
HFEngine
.
thAs
),
len
(
mus
)),
dtype
=
np
.
complex
)
radiusb
=
np
.
empty
((
len
(
self
.
HFEngine
.
thbs
),
len
(
mus
)),
dtype
=
np
.
complex
)
for
j
,
thA
in
enumerate
(
self
.
HFEngine
.
thAs
):
radiusA
[
j
]
=
expressionEvaluator
(
thA
[
0
],
muTestEff
)
for
j
,
thb
in
enumerate
(
self
.
HFEngine
.
thbs
):
radiusb
[
j
]
=
expressionEvaluator
(
thb
[
0
],
muTestEff
)
radiusA
=
np
.
expand_dims
(
uApproxRs
,
1
)
*
radiusA
ff
,
Lf
,
LL
=
self
.
_affineResidualMatricesContraction
(
radiusb
,
radiusA
)
err
=
np
.
abs
((
LL
-
2.
*
np
.
real
(
Lf
)
+
ff
)
/
ff
)
**
.
5
return
err
def
errorEstimator
(
self
,
mus
:
Np1D
,
return_max
:
bool
=
False
)
->
Np1D
:
setupOK
=
self
.
setupApproxLocal
()
if
setupOK
>
0
:
err
=
np
.
empty
(
len
(
mus
))
err
[:]
=
np
.
nan
if
not
return_max
:
return
err
return
err
,
[
-
setupOK
],
np
.
nan
mus
=
self
.
checkParameterList
(
mus
)
vbMng
(
self
.
trainedModel
,
"INIT"
,
"Evaluating error estimator at mu = {}."
.
format
(
mus
),
10
)
err
=
self
.
getErrorEstimatorAffine
(
mus
)
vbMng
(
self
.
trainedModel
,
"DEL"
,
"Done evaluating error estimator."
,
10
)
if
not
return_max
:
return
err
idxMaxEst
=
[
np
.
argmax
(
err
)]
return
err
,
idxMaxEst
,
err
[
idxMaxEst
]
def
_isLastSampleCollinear
(
self
)
->
bool
:
"""Check collinearity of last sample."""
if
self
.
collinearityTol
<=
0.
:
return
False
if
self
.
POD
==
1
:
reff
=
self
.
samplingEngine
.
Rscale
[:,
-
1
]
else
:
RROMPyWarning
((
"Repeated orthogonalization of the samples for "
"collinearity check. Consider setting POD to "
"True."
))
if
not
hasattr
(
self
,
"_PODEngine"
):
from
rrompy.sampling
import
PODEngine
self
.
_PODEngine
=
PODEngine
(
self
.
HFEngine
)
reff
=
self
.
_PODEngine
.
generalizedQR
(
self
.
samplingEngine
.
samples
,
only_R
=
True
,
is_state
=
True
)[:,
-
1
]
cLevel
=
np
.
abs
(
reff
[
-
1
])
/
np
.
linalg
.
norm
(
reff
)
cLevel
=
np
.
inf
if
np
.
isclose
(
cLevel
,
0.
,
atol
=
1e-15
)
else
1
/
cLevel
vbMng
(
self
,
"MAIN"
,
"Collinearity indicator {:.4e}."
.
format
(
cLevel
),
3
)
return
cLevel
>
self
.
collinearityTol
def
plotEstimator
(
self
,
est
:
Np1D
,
idxMax
:
List
[
int
],
estMax
:
List
[
float
]):
if
(
not
(
np
.
any
(
np
.
isnan
(
est
))
or
np
.
any
(
np
.
isinf
(
est
)))
and
masterCore
()):
fig
=
plt
.
figure
(
figsize
=
plt
.
figaspect
(
1.
/
self
.
npar
))
for
jpar
in
range
(
self
.
npar
):
ax
=
fig
.
add_subplot
(
1
,
self
.
npar
,
1
+
jpar
)
musre
=
np
.
array
(
self
.
muTest
.
re
.
data
)
errCP
=
copy
(
est
)
idx
=
np
.
delete
(
np
.
arange
(
self
.
npar
),
jpar
)
while
len
(
musre
)
>
0
:
if
self
.
npar
==
1
:
currIdx
=
np
.
arange
(
len
(
musre
))
else
:
currIdx
=
np
.
where
(
np
.
isclose
(
np
.
sum
(
np
.
abs
(
musre
[:,
idx
]
-
musre
[
0
,
idx
]),
1
),
0.
,
atol
=
1e-15
))[
0
]
ax
.
semilogy
(
musre
[
currIdx
,
jpar
],
errCP
[
currIdx
],
'k'
,
linewidth
=
1
)
musre
=
np
.
delete
(
musre
,
currIdx
,
0
)
errCP
=
np
.
delete
(
errCP
,
currIdx
)
ax
.
semilogy
([
self
.
muBounds
.
re
(
0
,
jpar
),
self
.
muBounds
.
re
(
-
1
,
jpar
)],
[
self
.
greedyTol
]
*
2
,
'r--'
)
ax
.
semilogy
(
self
.
mus
.
re
(
jpar
),
2.
*
self
.
greedyTol
*
np
.
ones
(
len
(
self
.
mus
)),
'*m'
)
if
len
(
idxMax
)
>
0
and
estMax
is
not
None
:
ax
.
semilogy
(
self
.
muTest
.
re
(
idxMax
,
jpar
),
estMax
,
'xr'
)
ax
.
set_xlim
(
*
list
(
self
.
sampler
.
lims
.
re
(
jpar
)))
ax
.
grid
()
plt
.
tight_layout
()
plt
.
show
()
def
greedyNextSample
(
self
,
muidx
:
int
,
plotEst
:
str
=
"NONE"
)
\
->
Tuple
[
Np1D
,
int
,
float
,
paramVal
]:
"""Compute next greedy snapshot of solution map."""
RROMPyAssert
(
self
.
_mode
,
message
=
"Cannot add greedy sample."
)
mus
=
copy
(
self
.
muTest
[
muidx
])
self
.
muTest
.
pop
(
muidx
)
for
j
,
mu
in
enumerate
(
mus
):
vbMng
(
self
,
"MAIN"
,
(
"Adding sample point no. {} at {} to training "
"set."
)
.
format
(
len
(
self
.
mus
)
+
1
,
mu
),
3
)
self
.
mus
.
append
(
mu
)
self
.
_S
=
len
(
self
.
mus
)
self
.
_approxParameters
[
"S"
]
=
self
.
S
if
(
self
.
samplingEngine
.
nsamples
<=
len
(
mus
)
-
j
-
1
or
not
np
.
allclose
(
mu
,
self
.
samplingEngine
.
mus
[
j
-
len
(
mus
)])):
self
.
samplingEngine
.
nextSample
(
mu
)
if
self
.
_isLastSampleCollinear
():
vbMng
(
self
,
"MAIN"
,
(
"Collinearity above tolerance detected. Starting "
"preemptive greedy loop termination."
),
3
)
self
.
_collinearityFlag
=
1
errorEstTest
=
np
.
empty
(
len
(
self
.
muTest
))
errorEstTest
[:]
=
np
.
nan
return
errorEstTest
,
[
-
1
],
np
.
nan
,
np
.
nan
errorEstTest
,
muidx
,
maxErrorEst
=
self
.
errorEstimator
(
self
.
muTest
,
True
)
if
plotEst
==
"ALL"
:
self
.
plotEstimator
(
errorEstTest
,
muidx
,
maxErrorEst
)
return
errorEstTest
,
muidx
,
maxErrorEst
,
self
.
muTest
[
muidx
]
def
_preliminaryTraining
(
self
):
"""Initialize starting snapshots of solution map."""
RROMPyAssert
(
self
.
_mode
,
message
=
"Cannot start greedy algorithm."
)
if
self
.
samplingEngine
.
nsamples
>
0
:
return
self
.
resetSamples
()
self
.
computeScaleFactor
()
self
.
samplingEngine
.
scaleFactor
=
self
.
scaleFactorDer
self
.
mus
=
self
.
samplerTrainSet
.
generatePoints
(
self
.
S
)
while
len
(
self
.
mus
)
>
self
.
S
:
self
.
mus
.
pop
()
muTestBase
=
self
.
sampler
.
generatePoints
(
self
.
nTestPoints
,
False
)
idxPop
=
pruneSamples
(
self
.
mapParameterList
(
muTestBase
),
self
.
mapParameterList
(
self
.
mus
),
1e-10
*
self
.
scaleFactor
[
0
])
muTestBase
.
pop
(
idxPop
)
muLast
=
copy
(
self
.
mus
[
-
1
])
self
.
mus
.
pop
()
if
len
(
self
.
mus
)
>
0
:
vbMng
(
self
,
"MAIN"
,
(
"Adding first {} sample point{} at {} to training "
"set."
)
.
format
(
self
.
S
-
1
,
""
+
"s"
*
(
self
.
S
>
2
),
self
.
mus
),
3
)
self
.
samplingEngine
.
iterSample
(
self
.
mus
)
self
.
_S
=
len
(
self
.
mus
)
self
.
_approxParameters
[
"S"
]
=
self
.
S
self
.
muTest
=
emptyParameterList
()
self
.
muTest
.
reset
((
len
(
muTestBase
)
+
1
,
self
.
mus
.
shape
[
1
]))
self
.
muTest
.
data
[:
-
1
]
=
muTestBase
.
data
self
.
muTest
.
data
[
-
1
]
=
muLast
.
data
@abstractmethod
def
setupApproxLocal
(
self
)
->
int
:
if
self
.
checkComputedApprox
():
return
-
1
RROMPyAssert
(
self
.
_mode
,
message
=
"Cannot setup approximant."
)
vbMng
(
self
,
"INIT"
,
"Setting up local approximant."
,
5
)
pass
vbMng
(
self
,
"DEL"
,
"Done setting up local approximant."
,
5
)
return
0
def
setupApprox
(
self
,
plotEst
:
str
=
"NONE"
)
->
int
:
"""Compute greedy snapshots of solution map."""
if
self
.
checkComputedApprox
():
return
-
1
RROMPyAssert
(
self
.
_mode
,
message
=
"Cannot start greedy algorithm."
)
vbMng
(
self
,
"INIT"
,
"Setting up {}."
.
format
(
self
.
name
()),
5
)
vbMng
(
self
,
"INIT"
,
"Starting computation of snapshots."
,
5
)
self
.
_collinearityFlag
=
0
self
.
_preliminaryTraining
()
muidx
,
self
.
firstGreedyIter
=
[
len
(
self
.
muTest
)
-
1
],
True
errorEstTest
,
maxErrorEst
=
[
np
.
inf
],
np
.
inf
max2ErrorEst
,
trainedModelOld
=
np
.
inf
,
None
while
self
.
firstGreedyIter
or
(
len
(
self
.
muTest
)
>
0
and
(
maxErrorEst
is
None
or
max2ErrorEst
>
self
.
greedyTol
)
and
self
.
samplingEngine
.
nsamples
<
self
.
maxIter
):
muTestOld
,
errorEstTestOld
=
self
.
muTest
,
errorEstTest
muidxOld
,
maxErrorEstOld
=
muidx
,
maxErrorEst
errorEstTest
,
muidx
,
maxErrorEst
,
mu
=
self
.
greedyNextSample
(
muidx
,
plotEst
)
if
maxErrorEst
is
not
None
and
(
np
.
any
(
np
.
isnan
(
maxErrorEst
))
or
np
.
any
(
np
.
isinf
(
maxErrorEst
))):
if
self
.
_collinearityFlag
==
0
and
not
self
.
firstGreedyIter
:
RROMPyWarning
((
"Instability in a posteriori "
"estimator. Starting preemptive greedy "
"loop termination."
))
self
.
muTest
,
errorEstTest
=
muTestOld
,
errorEstTestOld
if
self
.
firstGreedyIter
and
muidx
[
0
]
<
0
:
self
.
trainedModel
=
None
raise
RROMPyException
((
"Instability in approximant "
"computation. Aborting greedy "
"iterations."
))
self
.
_S
=
trainedModelOld
.
data
.
approxParameters
[
"S"
]
self
.
_approxParameters
[
"S"
]
=
self
.
S
while
self
.
samplingEngine
.
nsamples
>
self
.
S
:
self
.
samplingEngine
.
popSample
()
while
len
(
self
.
mus
)
>
self
.
S
:
self
.
mus
.
pop
(
-
1
)
muidx
,
maxErrorEst
=
muidxOld
,
maxErrorEstOld
break
if
maxErrorEst
is
not
None
:
max2ErrorEst
=
np
.
max
(
maxErrorEst
)
vbMng
(
self
,
"MAIN"
,
(
"Uniform testing error estimate "
"{:.4e}."
)
.
format
(
max2ErrorEst
),
5
)
if
self
.
firstGreedyIter
:
trainedModelOld
=
copy
(
self
.
trainedModel
)
else
:
trainedModelOld
.
data
=
copy
(
self
.
trainedModel
.
data
)
self
.
firstGreedyIter
=
False
vbMng
(
self
,
"DEL"
,
(
"Done computing snapshots (final snapshot count: "
"{})."
)
.
format
(
self
.
samplingEngine
.
nsamples
),
5
)
if
(
maxErrorEst
is
None
or
max2ErrorEst
<=
self
.
greedyTol
or
np
.
any
(
np
.
isnan
(
maxErrorEst
))
or
np
.
any
(
np
.
isinf
(
maxErrorEst
))):
while
self
.
samplingEngine
.
nsamples
>
self
.
S
:
self
.
samplingEngine
.
popSample
()
while
len
(
self
.
mus
)
>
self
.
S
:
self
.
mus
.
pop
(
-
1
)
else
:
self
.
_S
=
self
.
samplingEngine
.
nsamples
while
len
(
self
.
mus
)
<
self
.
S
:
self
.
mus
.
append
(
self
.
samplingEngine
.
mus
[
len
(
self
.
mus
)])
self
.
trainedModel
=
None
self
.
setupApproxLocal
()
if
plotEst
==
"LAST"
:
self
.
plotEstimator
(
errorEstTest
,
muidx
,
maxErrorEst
)
vbMng
(
self
,
"DEL"
,
"Done setting up approximant."
,
5
)
return
0
def
assembleReducedResidualGramian
(
self
,
pMat
:
sampList
):
"""
Build residual gramian of reduced linear system through projections.
"""
if
(
not
hasattr
(
self
.
trainedModel
.
data
,
"gramian"
)
or
self
.
trainedModel
.
data
.
gramian
is
None
):
gramian
=
self
.
HFEngine
.
innerProduct
(
pMat
,
pMat
,
dual
=
True
)
else
:
Sold
=
self
.
trainedModel
.
data
.
gramian
.
shape
[
0
]
S
=
len
(
self
.
mus
)
if
Sold
>
S
:
gramian
=
self
.
trainedModel
.
data
.
gramian
[:
S
,
:
S
]
else
:
idxOld
=
list
(
range
(
Sold
))
idxNew
=
list
(
range
(
Sold
,
S
))
gramian
=
np
.
empty
((
S
,
S
),
dtype
=
np
.
complex
)
gramian
[:
Sold
,
:
Sold
]
=
self
.
trainedModel
.
data
.
gramian
gramian
[:
Sold
,
Sold
:]
=
self
.
HFEngine
.
innerProduct
(
pMat
(
idxNew
),
pMat
(
idxOld
),
dual
=
True
)
gramian
[
Sold
:,
:
Sold
]
=
gramian
[:
Sold
,
Sold
:]
.
T
.
conj
()
gramian
[
Sold
:,
Sold
:]
=
self
.
HFEngine
.
innerProduct
(
pMat
(
idxNew
),
pMat
(
idxNew
),
dual
=
True
)
self
.
trainedModel
.
data
.
gramian
=
gramian
def
assembleReducedResidualBlocksbb
(
self
,
bs
:
List
[
Np1D
]):
"""
Build blocks (of type bb) of reduced linear system through projections.
"""
nbs
=
len
(
bs
)
if
(
not
hasattr
(
self
.
trainedModel
.
data
,
"resbb"
)
or
self
.
trainedModel
.
data
.
resbb
is
None
):
resbb
=
np
.
empty
((
nbs
,
nbs
),
dtype
=
np
.
complex
)
for
i
in
range
(
nbs
):
Mbi
=
bs
[
i
]
resbb
[
i
,
i
]
=
self
.
HFEngine
.
innerProduct
(
Mbi
,
Mbi
,
dual
=
True
)
for
j
in
range
(
i
):
Mbj
=
bs
[
j
]
resbb
[
i
,
j
]
=
self
.
HFEngine
.
innerProduct
(
Mbj
,
Mbi
,
dual
=
True
)
for
i
in
range
(
nbs
):
for
j
in
range
(
i
+
1
,
nbs
):
resbb
[
i
,
j
]
=
resbb
[
j
,
i
]
.
conj
()
self
.
trainedModel
.
data
.
resbb
=
resbb
def
assembleReducedResidualBlocksAb
(
self
,
As
:
List
[
Np2D
],
bs
:
List
[
Np1D
],
pMat
:
sampList
):
"""
Build blocks (of type Ab) of reduced linear system through projections.
"""
nAs
=
len
(
As
)
nbs
=
len
(
bs
)
S
=
len
(
self
.
mus
)
if
(
not
hasattr
(
self
.
trainedModel
.
data
,
"resAb"
)
or
self
.
trainedModel
.
data
.
resAb
is
None
):
if
isinstance
(
pMat
,
(
parameterList
,
sampleList
)):
pMat
=
pMat
.
data
resAb
=
np
.
empty
((
nbs
,
S
,
nAs
),
dtype
=
np
.
complex
)
for
j
in
range
(
nAs
):
MAj
=
dot
(
As
[
j
],
pMat
)
for
i
in
range
(
nbs
):
Mbi
=
bs
[
i
]
resAb
[
i
,
:,
j
]
=
self
.
HFEngine
.
innerProduct
(
MAj
,
Mbi
,
dual
=
True
)
else
:
Sold
=
self
.
trainedModel
.
data
.
resAb
.
shape
[
1
]
if
Sold
==
S
:
return
if
Sold
>
S
:
resAb
=
self
.
trainedModel
.
data
.
resAb
[:,
:
S
,
:]
else
:
if
isinstance
(
pMat
,
(
parameterList
,
sampleList
)):
pMat
=
pMat
.
data
resAb
=
np
.
empty
((
nbs
,
S
,
nAs
),
dtype
=
np
.
complex
)
resAb
[:,
:
Sold
,
:]
=
self
.
trainedModel
.
data
.
resAb
for
j
in
range
(
nAs
):
MAj
=
dot
(
As
[
j
],
pMat
[:,
Sold
:])
for
i
in
range
(
nbs
):
Mbi
=
bs
[
i
]
resAb
[
i
,
Sold
:,
j
]
=
self
.
HFEngine
.
innerProduct
(
MAj
,
Mbi
,
dual
=
True
)
self
.
trainedModel
.
data
.
resAb
=
resAb
def
assembleReducedResidualBlocksAA
(
self
,
As
:
List
[
Np2D
],
pMat
:
sampList
):
"""
Build blocks (of type AA) of reduced linear system through projections.
"""
nAs
=
len
(
As
)
S
=
len
(
self
.
mus
)
if
(
not
hasattr
(
self
.
trainedModel
.
data
,
"resAA"
)
or
self
.
trainedModel
.
data
.
resAA
is
None
):
if
isinstance
(
pMat
,
(
parameterList
,
sampleList
)):
pMat
=
pMat
.
data
resAA
=
np
.
empty
((
S
,
nAs
,
S
,
nAs
),
dtype
=
np
.
complex
)
for
i
in
range
(
nAs
):
MAi
=
dot
(
As
[
i
],
pMat
)
resAA
[:,
i
,
:,
i
]
=
self
.
HFEngine
.
innerProduct
(
MAi
,
MAi
,
dual
=
True
)
for
j
in
range
(
i
):
MAj
=
dot
(
As
[
j
],
pMat
)
resAA
[:,
i
,
:,
j
]
=
self
.
HFEngine
.
innerProduct
(
MAj
,
MAi
,
dual
=
True
)
for
i
in
range
(
nAs
):
for
j
in
range
(
i
+
1
,
nAs
):
resAA
[:,
i
,
:,
j
]
=
resAA
[:,
j
,
:,
i
]
.
T
.
conj
()
else
:
Sold
=
self
.
trainedModel
.
data
.
resAA
.
shape
[
0
]
if
Sold
==
S
:
return
if
Sold
>
S
:
resAA
=
self
.
trainedModel
.
data
.
resAA
[:
S
,
:,
:
S
,
:]
else
:
if
isinstance
(
pMat
,
(
parameterList
,
sampleList
)):
pMat
=
pMat
.
data
resAA
=
np
.
empty
((
S
,
nAs
,
S
,
nAs
),
dtype
=
np
.
complex
)
resAA
[:
Sold
,
:,
:
Sold
,
:]
=
self
.
trainedModel
.
data
.
resAA
for
i
in
range
(
nAs
):
MAi
=
dot
(
As
[
i
],
pMat
)
resAA
[:
Sold
,
i
,
Sold
:,
i
]
=
self
.
HFEngine
.
innerProduct
(
MAi
[:,
Sold
:],
MAi
[:,
:
Sold
],
dual
=
True
)
resAA
[
Sold
:,
i
,
:
Sold
,
i
]
=
resAA
[:
Sold
,
i
,
Sold
:,
i
]
.
T
.
conj
()
resAA
[
Sold
:,
i
,
Sold
:,
i
]
=
self
.
HFEngine
.
innerProduct
(
MAi
[:,
Sold
:],
MAi
[:,
Sold
:],
dual
=
True
)
for
j
in
range
(
i
):
MAj
=
dot
(
As
[
j
],
pMat
)
resAA
[:
Sold
,
i
,
Sold
:,
j
]
=
(
self
.
HFEngine
.
innerProduct
(
MAj
[:,
Sold
:],
MAi
[:,
:
Sold
],
dual
=
True
))
resAA
[
Sold
:,
i
,
:
Sold
,
j
]
=
(
self
.
HFEngine
.
innerProduct
(
MAj
[:,
:
Sold
],
MAi
[:,
Sold
:],
dual
=
True
))
resAA
[
Sold
:,
i
,
Sold
:,
j
]
=
(
self
.
HFEngine
.
innerProduct
(
MAj
[:,
Sold
:],
MAi
[:,
Sold
:],
dual
=
True
))
for
i
in
range
(
nAs
):
for
j
in
range
(
i
+
1
,
nAs
):
resAA
[:
Sold
,
i
,
Sold
:,
j
]
=
(
resAA
[
Sold
:,
j
,
:
Sold
,
i
]
.
T
.
conj
())
resAA
[
Sold
:,
i
,
:
Sold
,
j
]
=
(
resAA
[:
Sold
,
j
,
Sold
:,
i
]
.
T
.
conj
())
resAA
[
Sold
:,
i
,
Sold
:,
j
]
=
(
resAA
[
Sold
:,
j
,
Sold
:,
i
]
.
T
.
conj
())
self
.
trainedModel
.
data
.
resAA
=
resAA
def
assembleReducedResidualBlocks
(
self
,
full
:
bool
=
False
):
"""Build affine blocks of affine decomposition of residual."""
if
full
:
checkIfAffine
(
self
.
HFEngine
,
"assemble reduced residual blocks"
,
False
,
self
.
_affine_lvl
)
else
:
checkIfAffine
(
self
.
HFEngine
,
"assemble reduced RHS blocks"
,
True
,
self
.
_affine_lvl
)
self
.
HFEngine
.
buildb
()
self
.
assembleReducedResidualBlocksbb
(
self
.
HFEngine
.
bs
)
if
full
:
pMat
=
self
.
samplingEngine
.
projectionMatrix
self
.
HFEngine
.
buildA
()
self
.
assembleReducedResidualBlocksAb
(
self
.
HFEngine
.
As
,
self
.
HFEngine
.
bs
,
pMat
)
self
.
assembleReducedResidualBlocksAA
(
self
.
HFEngine
.
As
,
pMat
)
Event Timeline
Log In to Comment