Page Menu
Home
c4science
Search
Configure Global Search
Log In
Files
F80270701
sampling_engine.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, Aug 30, 20:24
Size
15 KB
Mime Type
text/x-python
Expires
Sun, Sep 1, 20:24 (1 d, 23 h)
Engine
blob
Format
Raw Data
Handle
20332011
Attached To
R6746 RationalROMPy
sampling_engine.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
numbers
import
Number
from
copy
import
deepcopy
as
copy
import
numpy
as
np
from
collections.abc
import
Iterable
from
warnings
import
catch_warnings
from
rrompy.utilities.base.types
import
(
Np1D
,
Np2D
,
HFEng
,
List
,
paramVal
,
Any
,
paramList
,
sampList
,
Tuple
,
TupleAny
,
DictAny
,
FigHandle
)
from
rrompy.utilities.base.data_structures
import
getNewFilename
from
rrompy.utilities.base
import
verbosityManager
as
vbMng
from
rrompy.utilities.exception_manager
import
(
RROMPyWarning
,
RROMPyException
,
RROMPyAssert
,
RROMPy_READY
,
RROMPy_FRAGILE
)
from
rrompy.utilities.numerical.hash_derivative
import
nextDerivativeIndices
from
rrompy.utilities.base.pickle_utilities
import
pickleDump
,
pickleLoad
from
rrompy.parameter
import
(
emptyParameterList
,
checkParameter
,
checkParameterList
)
from
rrompy.sampling
import
sampleList
,
emptySampleList
from
rrompy.utilities.parallel
import
bcast
,
masterCore
__all__
=
[
'SamplingEngine'
]
class
SamplingEngine
:
def
__init__
(
self
,
HFEngine
:
HFEng
,
verbosity
:
int
=
10
,
timestamp
:
bool
=
True
,
scaleFactor
:
Np1D
=
None
):
self
.
verbosity
=
verbosity
self
.
timestamp
=
timestamp
vbMng
(
self
,
"INIT"
,
"Initializing sampling engine of type {}."
.
format
(
self
.
name
()),
10
)
self
.
HFEngine
=
HFEngine
vbMng
(
self
,
"DEL"
,
"Done initializing sampling engine."
,
10
)
self
.
scaleFactor
=
scaleFactor
def
name
(
self
)
->
str
:
return
self
.
__class__
.
__name__
def
__str__
(
self
)
->
str
:
return
self
.
name
()
def
__repr__
(
self
)
->
str
:
return
self
.
__str__
()
+
" at "
+
hex
(
id
(
self
))
@property
def
HFEngine
(
self
):
"""Value of HFEngine. Its assignment resets history."""
return
self
.
_HFEngine
@HFEngine.setter
def
HFEngine
(
self
,
HFEngine
):
self
.
_HFEngine
=
HFEngine
self
.
resetHistory
()
@property
def
scaleFactor
(
self
):
"""Value of scaleFactor."""
return
self
.
_scaleFactor
@scaleFactor.setter
def
scaleFactor
(
self
,
scaleFactor
):
if
scaleFactor
is
None
:
scaleFactor
=
1.
if
not
isinstance
(
scaleFactor
,
Iterable
):
scaleFactor
=
[
scaleFactor
]
self
.
_scaleFactor
=
scaleFactor
def
scaleDer
(
self
,
derIdx
:
Np1D
):
if
not
isinstance
(
self
.
scaleFactor
,
Number
):
RROMPyAssert
(
len
(
derIdx
),
len
(
self
.
scaleFactor
),
"Number of derivative indices"
)
with
catch_warnings
(
record
=
True
)
as
w
:
res
=
np
.
prod
(
np
.
power
(
self
.
scaleFactor
,
derIdx
))
if
len
(
w
)
==
0
:
return
res
raise
RROMPyException
((
"Error in computing derivative scaling "
"factor: {}"
.
format
(
str
(
w
[
-
1
]
.
message
))))
@property
def
feature_keys
(
self
)
->
TupleAny
:
return
[
"mus"
,
"samples"
,
"nsamples"
,
"_derIdxs"
]
@property
def
feature_vals
(
self
)
->
DictAny
:
return
{
"mus"
:
self
.
mus
,
"samples"
:
self
.
samples
,
"nsamples"
:
self
.
nsamples
,
"_derIdxs"
:
self
.
_derIdxs
,
"_scaleFactor"
:
self
.
scaleFactor
}
def
_mergeFeature
(
self
,
name
:
str
,
value
:
Any
):
if
name
in
[
"mus"
,
"samples"
]:
getattr
(
self
,
name
)
.
append
(
value
)
elif
name
==
"nsamples"
:
self
.
nsamples
+=
value
elif
name
==
"_derIdxs"
:
self
.
_derIdxs
+=
value
else
:
raise
RROMPyException
((
"Invalid key {} in sampling engine "
"merge."
.
format
(
name
)))
def
store
(
self
,
filenameBase
:
str
=
"sampling_engine"
,
forceNewFile
:
bool
=
True
,
local
:
bool
=
False
)
->
str
:
"""Store sampling engine to file."""
filename
=
None
if
masterCore
():
vbMng
(
self
,
"INIT"
,
"Storing sampling engine to file."
,
20
)
if
forceNewFile
:
filename
=
getNewFilename
(
filenameBase
,
"pkl"
)
else
:
filename
=
"{}.pkl"
.
format
(
filenameBase
)
pickleDump
(
self
.
feature_vals
,
filename
)
vbMng
(
self
,
"DEL"
,
"Done storing engine."
,
20
)
if
local
:
return
filename
=
bcast
(
filename
)
return
filename
def
load
(
self
,
filename
:
str
,
merge
:
bool
=
False
):
"""Load sampling engine from file."""
if
isinstance
(
filename
,
(
list
,
tuple
,)):
self
.
load
(
filename
[
0
],
merge
)
for
filen
in
filename
[
1
:]:
self
.
load
(
filen
,
True
)
return
vbMng
(
self
,
"INIT"
,
"Loading stored sampling engine from file."
,
20
)
datadict
=
pickleLoad
(
filename
)
for
key
in
datadict
:
if
key
in
self
.
feature_keys
:
if
merge
and
key
!=
"_scaleFactor"
:
self
.
_mergeFeature
(
key
,
datadict
[
key
])
else
:
setattr
(
self
,
key
,
datadict
[
key
])
self
.
_mode
=
RROMPy_FRAGILE
vbMng
(
self
,
"DEL"
,
"Done loading stored engine."
,
20
)
@property
def
projectionMatrix
(
self
)
->
Np2D
:
return
self
.
samples
.
data
def
resetHistory
(
self
):
self
.
_mode
=
RROMPy_READY
self
.
samples
=
emptySampleList
()
self
.
nsamples
=
0
self
.
mus
=
emptyParameterList
()
self
.
_derIdxs
=
[]
def
setsample
(
self
,
u
:
sampList
,
overwrite
:
bool
=
False
):
if
overwrite
:
self
.
samples
[
self
.
nsamples
]
=
u
else
:
if
self
.
nsamples
==
0
:
self
.
samples
=
sampleList
(
u
)
else
:
self
.
samples
.
append
(
u
)
def
popSample
(
self
):
if
hasattr
(
self
,
"nsamples"
)
and
self
.
nsamples
>
1
:
if
self
.
samples
.
shape
[
1
]
>
self
.
nsamples
:
RROMPyWarning
((
"More than 'nsamples' memory allocated for "
"samples. Popping empty sample column."
))
self
.
nsamples
+=
1
self
.
nsamples
-=
1
self
.
samples
.
pop
()
self
.
mus
.
pop
()
else
:
self
.
resetHistory
()
def
preallocateSamples
(
self
,
u
:
sampList
,
mu
:
paramVal
,
n
:
int
):
self
.
_mode
=
RROMPy_READY
self
.
samples
.
reset
((
u
.
shape
[
0
],
n
),
u
.
dtype
)
self
.
samples
[
0
]
=
u
mu
=
checkParameter
(
mu
,
self
.
HFEngine
.
npar
)
self
.
mus
.
reset
((
n
,
self
.
HFEngine
.
npar
))
self
.
mus
[
0
]
=
mu
[
0
]
def
postprocessu
(
self
,
u
:
sampList
,
overwrite
:
bool
=
False
):
self
.
setsample
(
u
,
overwrite
)
def
postprocessuBulk
(
self
):
pass
def
solveLS
(
self
,
mu
:
paramList
=
[],
RHS
:
sampList
=
None
)
->
sampList
:
"""
Solve linear system.
Args:
mu: Parameter value.
Returns:
Solution of system.
"""
mu
=
checkParameterList
(
mu
,
self
.
HFEngine
.
npar
)
vbMng
(
self
,
"INIT"
,
"Solving HF model for mu = {}."
.
format
(
mu
),
15
)
u
=
self
.
HFEngine
.
solve
(
mu
,
RHS
,
return_state
=
True
)
vbMng
(
self
,
"DEL"
,
"Done solving HF model."
,
15
)
return
u
def
_getSampleConcurrence
(
self
,
mu
:
paramVal
,
previous
:
Np1D
)
->
sampList
:
"""
Compute sample after checking if it is a derivative.
Args:
mu: Parameter value.
previous: Indices of previous unique samples.
Returns:
Snapshot.
"""
RROMPyAssert
(
self
.
_mode
,
message
=
"Cannot add samples."
)
from
rrompy.utilities.numerical
import
dot
if
len
(
previous
)
>=
len
(
self
.
_derIdxs
):
self
.
_derIdxs
+=
nextDerivativeIndices
(
self
.
_derIdxs
,
len
(
self
.
scaleFactor
),
len
(
previous
)
+
1
-
len
(
self
.
_derIdxs
))
derIdx
=
self
.
_derIdxs
[
len
(
previous
)]
mu
=
checkParameter
(
mu
,
self
.
HFEngine
.
npar
)
samplesOld
=
self
.
samples
(
previous
)
RHS
=
self
.
scaleDer
(
derIdx
)
*
self
.
HFEngine
.
b
(
mu
,
derIdx
)
for
j
,
derP
in
enumerate
(
self
.
_derIdxs
[:
len
(
previous
)]):
diffP
=
[
x
-
y
for
(
x
,
y
)
in
zip
(
derIdx
,
derP
)]
if
np
.
all
([
x
>=
0
for
x
in
diffP
]):
RHS
-=
self
.
scaleDer
(
diffP
)
*
dot
(
self
.
HFEngine
.
A
(
mu
,
diffP
),
samplesOld
[
j
])
return
self
.
solveLS
(
mu
,
RHS
=
RHS
)
def
nextSample
(
self
,
mu
:
paramVal
,
overwrite
:
bool
=
False
,
postprocess
:
bool
=
True
)
->
Np1D
:
"""
Compute one sample.
Args:
mu: Parameter value.
overwrite(optional): Whether to overwrite sample in self.samples.
Defaults to False.
postprocess(optional): Whether to perform post-processing step.
Defaults to True.
Returns:
Snapshot.
"""
RROMPyAssert
(
self
.
_mode
,
message
=
"Cannot add samples."
)
mu
=
checkParameter
(
mu
,
self
.
HFEngine
.
npar
)
muidxs
=
self
.
mus
.
findall
(
mu
[
0
])
if
len
(
muidxs
)
>
0
:
u
=
self
.
_getSampleConcurrence
(
mu
,
np
.
sort
(
muidxs
))
else
:
u
=
self
.
solveLS
(
mu
)
if
postprocess
:
self
.
postprocessu
(
u
,
overwrite
=
overwrite
)
else
:
self
.
setsample
(
u
,
overwrite
)
if
overwrite
:
self
.
mus
[
self
.
nsamples
]
=
mu
[
0
]
else
:
self
.
mus
.
append
(
mu
)
self
.
nsamples
+=
1
return
self
.
samples
[
self
.
nsamples
-
1
]
def
iterSample
(
self
,
mus
:
paramList
)
->
sampList
:
"""
Compute set of samples.
Args:
mus: Parameter values.
Returns:
Snapshots.
"""
mus
=
checkParameterList
(
mus
,
self
.
HFEngine
.
npar
)
vbMng
(
self
,
"INIT"
,
"Starting sampling iterations."
,
5
)
n
=
len
(
mus
)
if
n
<=
0
:
raise
RROMPyException
((
"Number of samples must be positive."
))
self
.
resetHistory
()
if
len
(
mus
.
unique
())
!=
n
:
for
j
in
range
(
n
):
vbMng
(
self
,
"MAIN"
,
"Computing sample {} / {}."
.
format
(
j
+
1
,
n
),
7
)
self
.
nextSample
(
mus
[
j
],
overwrite
=
(
j
>
0
),
postprocess
=
False
)
if
n
>
1
and
j
==
0
:
self
.
preallocateSamples
(
self
.
samples
[
0
],
mus
[
0
],
n
)
else
:
self
.
setsample
(
self
.
solveLS
(
mus
),
overwrite
=
False
)
self
.
mus
=
copy
(
mus
)
self
.
nsamples
=
n
self
.
postprocessuBulk
()
vbMng
(
self
,
"DEL"
,
"Finished sampling iterations."
,
5
)
return
self
.
samples
def
plotSamples
(
self
,
warpings
:
List
[
List
[
callable
]]
=
None
,
name
:
str
=
"u"
,
**
kwargs
)
->
Tuple
[
List
[
FigHandle
],
List
[
str
]]:
"""
Do some nice plots of the samples.
Args:
warpings(optional): Domain warping functions.
name(optional): Name to be shown as title of the plots. Defaults to
'u'.
Returns:
Output filenames and figure handles.
"""
if
warpings
is
None
:
warpings
=
[
None
]
*
self
.
nsamples
figs
=
[
None
]
*
self
.
nsamples
filesOut
=
[
None
]
*
self
.
nsamples
for
j
in
range
(
self
.
nsamples
):
pltOut
=
self
.
HFEngine
.
plot
(
self
.
samples
[
j
],
warpings
[
j
],
True
,
"{}_{}"
.
format
(
name
,
j
),
**
kwargs
)
if
isinstance
(
pltOut
,
(
tuple
,)):
figs
[
j
],
filesOut
[
j
]
=
pltOut
else
:
figs
[
j
]
=
pltOut
if
filesOut
[
0
]
is
None
:
return
figs
return
figs
,
filesOut
def
outParaviewSamples
(
self
,
warpings
:
List
[
List
[
callable
]]
=
None
,
name
:
str
=
"u"
,
filename
:
str
=
"out"
,
times
:
Np1D
=
None
,
**
kwargs
)
->
List
[
str
]:
"""
Output samples to ParaView file.
Args:
warpings(optional): Domain warping functions.
name(optional): Base name to be used for data output.
filename(optional): Name of output file.
times(optional): Timestamps.
Returns:
Output filenames.
"""
if
warpings
is
None
:
warpings
=
[
None
]
*
self
.
nsamples
if
times
is
None
:
times
=
[
0.
]
*
self
.
nsamples
filesOut
=
[
None
]
*
self
.
nsamples
for
j
in
range
(
self
.
nsamples
):
filesOut
[
j
]
=
self
.
HFEngine
.
outParaview
(
self
.
samples
[
j
],
warpings
[
j
],
True
,
"{}_{}"
.
format
(
name
,
j
),
"{}_{}"
.
format
(
filename
,
j
),
times
[
j
],
**
kwargs
)
if
filesOut
[
0
]
is
None
:
return
None
return
filesOut
def
outParaviewTimeDomainSamples
(
self
,
omegas
:
Np1D
=
None
,
warpings
:
List
[
List
[
callable
]]
=
None
,
timeFinal
:
Np1D
=
None
,
periodResolution
:
List
[
int
]
=
20
,
name
:
str
=
"u"
,
filename
:
str
=
"out"
,
**
kwargs
)
->
List
[
str
]:
"""
Output samples to ParaView file, converted to time domain.
Args:
omegas(optional): frequencies.
timeFinal(optional): final time of simulation.
periodResolution(optional): number of time steps per period.
name(optional): Base name to be used for data output.
filename(optional): Name of output file.
Returns:
Output filename.
"""
if
omegas
is
None
:
omegas
=
np
.
real
(
self
.
mus
)
if
warpings
is
None
:
warpings
=
[
None
]
*
self
.
nsamples
if
not
isinstance
(
timeFinal
,
(
list
,
tuple
,)):
timeFinal
=
[
timeFinal
]
*
self
.
nsamples
if
not
isinstance
(
periodResolution
,
(
list
,
tuple
,)):
periodResolution
=
[
periodResolution
]
*
self
.
nsamples
filesOut
=
[
None
]
*
self
.
nsamples
for
j
in
range
(
self
.
nsamples
):
filesOut
[
j
]
=
self
.
HFEngine
.
outParaviewTimeDomain
(
self
.
samples
[
j
],
omegas
[
j
],
warpings
[
j
],
True
,
timeFinal
[
j
],
periodResolution
[
j
],
"{}_{}"
.
format
(
name
,
j
),
"{}_{}"
.
format
(
filename
,
j
),
**
kwargs
)
if
filesOut
[
0
]
is
None
:
return
None
return
filesOut
Event Timeline
Log In to Comment