Page Menu
Home
c4science
Search
Configure Global Search
Log In
Files
F61807183
rational_pade.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
Thu, May 9, 02:24
Size
17 KB
Mime Type
text/x-python
Expires
Sat, May 11, 02:24 (2 d)
Engine
blob
Format
Raw Data
Handle
17561036
Attached To
R6746 RationalROMPy
rational_pade.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/>.
#
from
copy
import
deepcopy
as
copy
import
numpy
as
np
from
rrompy.reduction_methods.base
import
checkRobustTolerance
from
rrompy.reduction_methods.trained_model
import
(
TrainedModelData
,
TrainedModelPade
as
tModel
)
from
.generic_centered_approximant
import
GenericCenteredApproximant
from
rrompy.utilities.base.types
import
(
Np1D
,
Np2D
,
Tuple
,
DictAny
,
HFEng
,
paramVal
,
paramList
,
sampList
)
from
rrompy.utilities.base
import
verbosityDepth
,
purgeDict
from
rrompy.utilities.exception_manager
import
(
RROMPyException
,
RROMPyAssert
,
RROMPyWarning
)
__all__
=
[
'RationalPade'
]
class
RationalPade
(
GenericCenteredApproximant
):
"""
ROM single-point fast Pade' approximant computation for parametric
problems.
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': whether to compute POD of snapshots; defaults to True;
- 'rho': weight for computation of original Pade' approximant;
defaults to np.inf, i.e. fast approximant;
- 'M': degree of Pade' approximant numerator; defaults to 0;
- 'N': degree of Pade' approximant denominator; defaults to 0;
- 'E': total number of derivatives current approximant relies upon;
defaults to 1;
- 'robustTol': tolerance for robust Pade' denominator management;
defaults to 0.
Defaults to empty dict.
homogeneized(optional): Whether to homogeneize Dirichlet BCs. Defaults
to False.
verbosity(optional): Verbosity level. Defaults to 10.
Attributes:
HFEngine: HF problem solver.
mu0: Default parameter.
homogeneized: Whether to homogeneize Dirichlet BCs.
approxParameters: Dictionary containing values for main parameters of
approximant. Recognized keys are in parameterList.
parameterList: Recognized keys of approximant parameters:
- 'POD': whether to compute POD of snapshots;
- 'rho': weight for computation of original Pade' approximant;
- 'M': degree of Pade' approximant numerator;
- 'N': degree of Pade' approximant denominator;
- 'E': total number of derivatives current approximant relies upon;
- 'robustTol': tolerance for robust Pade' denominator management.
POD: Whether to compute QR factorization of derivatives.
rho: Weight of approximant.
M: Numerator degree of approximant.
N: Denominator degree of approximant.
E: Number of solution derivatives over which current approximant is
based upon.
robustTol: Tolerance for robust Pade' denominator management.
initialHFData: HF problem initial data.
uHF: High fidelity solution with wavenumber lastSolvedHF as numpy
complex vector.
lastSolvedHF: Wavenumber corresponding to last computed high fidelity
solution.
G: Square Numpy 2D vector of size (N+1) corresponding to Pade'
denominator matrix (see paper).
uApp: Last evaluated approximant as numpy complex vector.
"""
def
__init__
(
self
,
HFEngine
:
HFEng
,
mu0
:
paramVal
=
0
,
approxParameters
:
DictAny
=
{},
homogeneized
:
bool
=
False
,
verbosity
:
int
=
10
,
timestamp
:
bool
=
True
):
self
.
_preInit
()
self
.
_addParametersToList
([
"M"
,
"N"
,
"robustTol"
,
"rho"
])
super
()
.
__init__
(
HFEngine
=
HFEngine
,
mu0
=
mu0
,
approxParameters
=
approxParameters
,
homogeneized
=
homogeneized
,
verbosity
=
verbosity
,
timestamp
=
timestamp
)
self
.
_postInit
()
@property
def
approxParameters
(
self
):
"""Value of approximant parameters."""
return
self
.
_approxParameters
@approxParameters.setter
def
approxParameters
(
self
,
approxParams
):
approxParameters
=
purgeDict
(
approxParams
,
self
.
parameterList
,
dictname
=
self
.
name
()
+
".approxParameters"
,
baselevel
=
1
)
approxParametersCopy
=
purgeDict
(
approxParameters
,
[
"M"
,
"N"
,
"robustTol"
,
"rho"
],
True
,
True
,
baselevel
=
1
)
keyList
=
list
(
approxParameters
.
keys
())
if
"rho"
in
keyList
:
self
.
_rho
=
approxParameters
[
"rho"
]
elif
not
hasattr
(
self
,
"_rho"
)
or
self
.
rho
is
None
:
self
.
_rho
=
np
.
inf
GenericCenteredApproximant
.
approxParameters
.
fset
(
self
,
approxParametersCopy
)
self
.
rho
=
self
.
_rho
if
"robustTol"
in
keyList
:
self
.
robustTol
=
approxParameters
[
"robustTol"
]
elif
not
hasattr
(
self
,
"_robustTol"
)
or
self
.
_robustTol
is
None
:
self
.
robustTol
=
0
self
.
_ignoreParWarnings
=
True
if
"M"
in
keyList
:
self
.
M
=
approxParameters
[
"M"
]
elif
hasattr
(
self
,
"_M"
)
and
self
.
_M
is
not
None
:
self
.
M
=
self
.
M
else
:
self
.
M
=
0
del
self
.
_ignoreParWarnings
if
"N"
in
keyList
:
self
.
N
=
approxParameters
[
"N"
]
elif
hasattr
(
self
,
"_N"
)
and
self
.
_N
is
not
None
:
self
.
N
=
self
.
N
else
:
self
.
N
=
0
@property
def
rho
(
self
):
"""Value of rho."""
return
self
.
_rho
@rho.setter
def
rho
(
self
,
rho
):
self
.
_rho
=
np
.
abs
(
rho
)
self
.
_approxParameters
[
"rho"
]
=
self
.
rho
@property
def
M
(
self
):
"""Value of M. Its assignment may change E."""
return
self
.
_M
@M.setter
def
M
(
self
,
M
):
if
M
<
0
:
raise
RROMPyException
(
"M must be non-negative."
)
self
.
_M
=
M
self
.
_approxParameters
[
"M"
]
=
self
.
M
if
not
hasattr
(
self
,
"_ignoreParWarnings"
):
self
.
checkMNE
()
@property
def
N
(
self
):
"""Value of N. Its assignment may change E."""
return
self
.
_N
@N.setter
def
N
(
self
,
N
):
if
N
<
0
:
raise
RROMPyException
(
"N must be non-negative."
)
self
.
_N
=
N
self
.
_approxParameters
[
"N"
]
=
self
.
N
if
not
hasattr
(
self
,
"_ignoreParWarnings"
):
self
.
checkMNE
()
def
checkMNE
(
self
):
"""Check consistency of M, N, and E."""
if
not
hasattr
(
self
,
"_E"
)
or
self
.
E
is
None
:
return
M
=
self
.
M
if
(
hasattr
(
self
,
"_M"
)
and
self
.
M
is
not
None
)
else
0
N
=
self
.
N
if
(
hasattr
(
self
,
"_N"
)
and
self
.
N
is
not
None
)
else
0
msg
=
"max(M, N)"
if
self
.
rho
==
np
.
inf
else
"M + N"
bound
=
eval
(
msg
)
if
self
.
E
<
bound
:
RROMPyWarning
((
"Prescribed E is too small. Updating E to "
"{}."
)
.
format
(
msg
))
self
.
E
=
bound
del
M
,
N
@property
def
robustTol
(
self
):
"""Value of tolerance for robust Pade' denominator management."""
return
self
.
_robustTol
@robustTol.setter
def
robustTol
(
self
,
robustTol
):
if
robustTol
<
0.
:
RROMPyWarning
((
"Overriding prescribed negative robustness "
"tolerance to 0."
))
robustTol
=
0.
self
.
_robustTol
=
robustTol
self
.
_approxParameters
[
"robustTol"
]
=
self
.
robustTol
@property
def
E
(
self
):
"""Value of E."""
return
self
.
_E
@E.setter
def
E
(
self
,
E
):
GenericCenteredApproximant
.
E
.
fset
(
self
,
E
)
self
.
checkMNE
()
def
_setupDenominator
(
self
):
"""Compute Pade' denominator."""
if
self
.
verbosity
>=
7
:
verbosityDepth
(
"INIT"
,
"Starting computation of denominator."
,
timestamp
=
self
.
timestamp
)
while
self
.
N
>
0
:
if
self
.
POD
:
ev
,
eV
=
self
.
findeveVGQR
()
else
:
ev
,
eV
=
self
.
findeveVGExplicit
()
newParameters
=
checkRobustTolerance
(
ev
,
self
.
E
,
self
.
robustTol
)
if
not
newParameters
:
break
self
.
approxParameters
=
newParameters
if
self
.
N
<=
0
:
eV
=
np
.
ones
((
1
,
1
))
if
self
.
verbosity
>=
7
:
verbosityDepth
(
"DEL"
,
"Done computing denominator."
,
timestamp
=
self
.
timestamp
)
return
eV
[::
-
1
,
0
]
def
_setupNumerator
(
self
):
"""Compute Pade' numerator."""
if
self
.
verbosity
>=
7
:
verbosityDepth
(
"INIT"
,
"Starting computation of numerator."
,
timestamp
=
self
.
timestamp
)
P
=
np
.
zeros
((
self
.
E
+
1
,
self
.
M
+
1
),
dtype
=
np
.
complex
)
for
i
in
range
(
self
.
E
+
1
):
l
=
min
(
self
.
M
+
1
,
i
+
self
.
N
+
1
)
if
i
<
l
:
P
[
i
,
i
:
l
]
=
self
.
trainedModel
.
data
.
Q
[:
l
-
i
]
if
self
.
verbosity
>=
7
:
verbosityDepth
(
"DEL"
,
"Done computing numerator."
,
timestamp
=
self
.
timestamp
)
return
self
.
rescaleParameter
(
P
.
T
)
.
T
def
setupApprox
(
self
):
"""
Compute Pade' approximant. SVD-based robust eigenvalue management.
"""
if
self
.
checkComputedApprox
():
return
if
self
.
verbosity
>=
5
:
verbosityDepth
(
"INIT"
,
"Setting up {}."
.
format
(
self
.
name
()),
timestamp
=
self
.
timestamp
)
self
.
computeDerivatives
()
if
self
.
trainedModel
is
None
:
self
.
trainedModel
=
tModel
()
self
.
trainedModel
.
verbosity
=
self
.
verbosity
self
.
trainedModel
.
timestamp
=
self
.
timestamp
data
=
TrainedModelData
(
self
.
trainedModel
.
name
(),
self
.
mu0
,
None
,
self
.
HFEngine
.
rescalingExp
)
data
.
polytype
=
"MONOMIAL"
self
.
trainedModel
.
data
=
data
if
self
.
N
>
0
:
Q
=
self
.
_setupDenominator
()
else
:
Q
=
np
.
ones
(
1
,
dtype
=
np
.
complex
)
self
.
trainedModel
.
data
.
Q
=
copy
(
Q
)
self
.
trainedModel
.
data
.
scaleFactor
=
self
.
scaleFactor
self
.
trainedModel
.
data
.
projMat
=
copy
(
self
.
samplingEngine
.
samples
(
list
(
range
(
self
.
E
+
1
))))
P
=
self
.
_setupNumerator
()
if
self
.
POD
:
P
=
self
.
samplingEngine
.
RPOD
.
dot
(
P
)
self
.
trainedModel
.
data
.
P
=
copy
(
P
)
self
.
trainedModel
.
data
.
approxParameters
=
copy
(
self
.
approxParameters
)
if
self
.
verbosity
>=
5
:
verbosityDepth
(
"DEL"
,
"Done setting up approximant."
,
timestamp
=
self
.
timestamp
)
def
rescaleParameter
(
self
,
R
:
Np2D
,
A
:
Np2D
=
None
,
exponent
:
float
=
1.
)
->
Np2D
:
"""
Prepare parameter rescaling.
Args:
R: Matrix whose columns need rescaling.
A(optional): Matrix whose diagonal defines scaling factor. If None,
previous value of scaleFactor is used. Defaults to None.
exponent(optional): Exponent of scaling factor in matrix diagonal.
Defaults to 1.
Returns:
Rescaled matrix.
"""
RROMPyAssert
(
self
.
_mode
,
message
=
"Cannot compute rescaling factor."
)
if
A
is
not
None
:
aDiag
=
np
.
diag
(
A
)
scaleCoeffs
=
np
.
polyfit
(
np
.
arange
(
A
.
shape
[
1
]),
np
.
log
(
aDiag
),
1
)
self
.
scaleFactor
=
np
.
exp
(
-
scaleCoeffs
[
0
]
/
exponent
)
return
R
*
np
.
power
(
self
.
scaleFactor
,
np
.
arange
(
R
.
shape
[
1
]))
def
buildG
(
self
):
"""Assemble Pade' denominator matrix."""
RROMPyAssert
(
self
.
_mode
,
message
=
"Cannot compute G matrix."
)
self
.
computeDerivatives
()
if
self
.
verbosity
>=
10
:
verbosityDepth
(
"INIT"
,
"Building gramian matrix."
,
timestamp
=
self
.
timestamp
)
if
self
.
rho
==
np
.
inf
:
Nmin
=
self
.
E
-
self
.
N
else
:
Nmin
=
self
.
M
-
self
.
N
+
1
if
self
.
POD
:
RPODE
=
self
.
samplingEngine
.
RPOD
[:
self
.
E
+
1
,
Nmin
:
self
.
E
+
1
]
RPODE
=
self
.
rescaleParameter
(
RPODE
,
RPODE
[
Nmin
:,
:])
G
=
RPODE
.
T
.
conj
()
.
dot
(
RPODE
)
else
:
DerE
=
self
.
samplingEngine
.
samples
(
list
(
range
(
Nmin
,
self
.
E
+
1
)))
G
=
self
.
HFEngine
.
innerProduct
(
DerE
,
DerE
)
DerE
=
self
.
rescaleParameter
(
DerE
,
G
,
2.
)
G
=
self
.
HFEngine
.
innerProduct
(
DerE
,
DerE
)
if
self
.
rho
==
np
.
inf
:
self
.
G
=
G
else
:
Gbig
=
G
self
.
G
=
np
.
zeros
((
self
.
N
+
1
,
self
.
N
+
1
),
dtype
=
np
.
complex
)
for
k
in
range
(
self
.
E
-
self
.
M
):
self
.
G
+=
self
.
rho
**
(
2
*
k
)
*
Gbig
[
k
:
k
+
self
.
N
+
1
,
k
:
k
+
self
.
N
+
1
]
if
self
.
verbosity
>=
10
:
verbosityDepth
(
"DEL"
,
"Done building gramian."
,
timestamp
=
self
.
timestamp
)
def
findeveVGExplicit
(
self
)
->
Tuple
[
Np1D
,
Np2D
]:
"""
Compute explicitly eigenvalues and eigenvectors of Pade' denominator
matrix.
"""
RROMPyAssert
(
self
.
_mode
,
message
=
"Cannot solve eigenvalue problem."
)
self
.
buildG
()
if
self
.
verbosity
>=
7
:
verbosityDepth
(
"INIT"
,
"Solving eigenvalue problem for gramian matrix."
,
timestamp
=
self
.
timestamp
)
ev
,
eV
=
np
.
linalg
.
eigh
(
self
.
G
)
if
self
.
verbosity
>=
5
:
try
:
condev
=
ev
[
-
1
]
/
ev
[
0
]
except
:
condev
=
np
.
inf
verbosityDepth
(
"MAIN"
,
(
"Solved eigenvalue problem of size {} "
"with condition number {:.4e}."
)
.
format
(
self
.
N
+
1
,
condev
),
timestamp
=
self
.
timestamp
)
if
self
.
verbosity
>=
7
:
verbosityDepth
(
"DEL"
,
"Done solving eigenvalue problem."
,
timestamp
=
self
.
timestamp
)
return
ev
,
eV
def
_buildRStack
(
self
,
R
:
Np2D
)
->
Np2D
:
if
self
.
verbosity
>=
10
:
verbosityDepth
(
"INIT"
,
(
"Building matrix stack for square "
"root of gramian."
),
timestamp
=
self
.
timestamp
)
REff
=
np
.
zeros
((
R
.
shape
[
0
]
*
(
self
.
E
-
self
.
M
),
self
.
N
+
1
),
dtype
=
np
.
complex
)
for
k
in
range
(
self
.
E
-
self
.
M
):
RTleft
=
max
(
0
,
self
.
N
-
self
.
M
-
k
)
Rleft
=
max
(
0
,
self
.
M
-
self
.
N
+
k
)
REff
[
k
*
R
.
shape
[
0
]
:
(
k
+
1
)
*
R
.
shape
[
0
],
RTleft
:]
=
(
self
.
rho
**
k
*
R
[:,
Rleft
:
self
.
M
+
1
+
k
])
if
self
.
verbosity
>=
10
:
verbosityDepth
(
"DEL"
,
"Done building matrix stack."
,
timestamp
=
self
.
timestamp
)
return
REff
def
findeveVGQR
(
self
)
->
Tuple
[
Np1D
,
Np2D
]:
"""
Compute eigenvalues and eigenvectors of Pade' denominator matrix
through SVD of R factor.
Returns:
Eigenvalues in ascending order and corresponding eigenvector
matrix.
"""
RROMPyAssert
(
self
.
_mode
,
message
=
"Cannot solve eigenvalue problem."
)
RROMPyAssert
(
self
.
POD
,
True
,
"POD value"
)
self
.
computeDerivatives
()
if
self
.
rho
==
np
.
inf
:
Nmin
=
self
.
E
-
self
.
N
else
:
Nmin
=
self
.
M
-
self
.
N
+
1
R
=
self
.
samplingEngine
.
RPOD
[:
self
.
E
+
1
,
Nmin
:
self
.
E
+
1
]
R
=
self
.
rescaleParameter
(
R
,
R
[
R
.
shape
[
0
]
-
R
.
shape
[
1
]
:,
:])
REff
=
R
if
self
.
rho
==
np
.
inf
else
self
.
_buildRStack
(
R
)
if
self
.
verbosity
>=
7
:
verbosityDepth
(
"INIT"
,
(
"Solving svd for square root of "
"gramian matrix."
),
timestamp
=
self
.
timestamp
)
sizeI
=
REff
.
shape
[
0
]
_
,
s
,
V
=
np
.
linalg
.
svd
(
REff
,
full_matrices
=
False
)
eV
=
V
[::
-
1
,
:]
.
T
.
conj
()
if
self
.
verbosity
>=
5
:
try
:
condev
=
s
[
0
]
/
s
[
-
1
]
except
:
condev
=
np
.
inf
verbosityDepth
(
"MAIN"
,
(
"Solved svd problem of size {} x {} with "
"condition number {:.4e}."
)
.
format
(
sizeI
,
self
.
N
+
1
,
condev
),
timestamp
=
self
.
timestamp
)
if
self
.
verbosity
>=
7
:
verbosityDepth
(
"DEL"
,
"Done solving eigenvalue problem."
,
timestamp
=
self
.
timestamp
)
return
s
[::
-
1
],
eV
def
centerNormalize
(
self
,
mu
:
paramList
,
mu0
:
paramVal
=
None
)
->
paramList
:
"""
Compute normalized parameter to be plugged into approximant.
Args:
mu: Parameter(s) 1.
mu0: Parameter(s) 2. If None, set to self.mu0.
Returns:
Normalized parameter.
"""
return
self
.
trainedModel
.
centerNormalize
(
mu
,
mu0
)
def
getResidues
(
self
)
->
sampList
:
"""
Obtain approximant residues.
Returns:
Matrix with residues as columns.
"""
return
self
.
trainedModel
.
getResidues
()
Event Timeline
Log In to Comment