Page Menu
Home
c4science
Search
Configure Global Search
Log In
Files
F60192563
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
Sun, Apr 28, 04:52
Size
16 KB
Mime Type
text/x-python
Expires
Tue, Apr 30, 04:52 (2 d)
Engine
blob
Format
Raw Data
Handle
17315572
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
.rational_interpolant
import
RationalInterpolant
from
rrompy.utilities.poly_fitting.polynomial
import
(
polybases
as
ppb
,
polyfitname
,
polyvander
as
pvP
,
polyvanderTotal
as
pvTP
,
polyTimesTable
,
vanderInvTable
,
PolynomialInterpolator
as
PI
)
from
rrompy.utilities.poly_fitting.radial_basis
import
(
polybases
as
rbpb
,
RadialBasisInterpolator
as
RBI
)
from
rrompy.utilities.poly_fitting.moving_least_squares
import
(
MovingLeastSquaresInterpolator
as
MLSI
)
from
rrompy.utilities.base.types
import
Np2D
,
Tuple
,
List
from
rrompy.utilities.base
import
verbosityManager
as
vbMng
from
rrompy.utilities.numerical
import
customPInv
,
dot
from
rrompy.utilities.numerical.degree
import
(
fullDegreeN
,
totalDegreeN
,
reduceDegreeN
,
degreeTotalToFull
,
fullDegreeMaxMask
,
totalDegreeMaxMask
)
from
rrompy.utilities.exception_manager
import
(
RROMPyException
,
RROMPyAssert
,
RROMPyWarning
)
__all__
=
[
'RationalPade'
]
class
RationalPade
(
RationalInterpolant
):
"""
ROM rational Pade' 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;
- 'scaleFactorDer': scaling factors for derivative computation;
defaults to 'AUTO';
- 'S': total number of samples current approximant relies upon;
- 'sampler': sample point generator;
- 'polybasis': type of polynomial basis for interpolation; defaults
to 'MONOMIAL';
- 'M': degree of rational interpolant numerator; defaults to
'AUTO', i.e. maximum allowed;
- 'N': degree of rational interpolant denominator; defaults to
'AUTO', i.e. maximum allowed;
- 'polydegreetype': type of polynomial degree; defaults to 'TOTAL';
- 'radialDirectionalWeights': radial basis weights for interpolant
numerator; defaults to 1;
- 'nNearestNeighbor': mumber of nearest neighbors considered in
numerator if polybasis allows; defaults to -1;
- 'interpRcond': tolerance for interpolation; defaults to None;
- 'robustTol': tolerance for robust rational denominator
management; defaults to 0;
- 'correctorForce': whether corrector should forcefully delete bad
poles; defaults to False;
- 'correctorTol': tolerance for corrector step; defaults to 0.,
i.e. no bad poles;
- 'correctorMaxIter': maximum number of corrector iterations;
defaults to 1.
Defaults to empty dict.
approx_state(optional): Whether to approximate state. Defaults to
False.
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': whether to compute POD of snapshots;
- 'scaleFactorDer': scaling factors for derivative computation;
- 'polybasis': type of polynomial basis for interpolation;
- 'M': degree of rational interpolant numerator;
- 'N': degree of rational interpolant denominator;
- 'polydegreetype': type of polynomial degree;
- 'radialDirectionalWeights': radial basis weights for interpolant
numerator;
- 'nNearestNeighbor': mumber of nearest neighbors considered in
numerator if polybasis allows;
- 'interpRcond': tolerance for interpolation via numpy.polyfit;
- 'robustTol': tolerance for robust rational denominator
management;
- 'correctorForce': whether corrector should forcefully delete bad
poles;
- 'correctorTol': tolerance for corrector step;
- 'correctorMaxIter': maximum number of corrector iterations.
parameterListCritical: Recognized keys of critical approximant
parameters:
- 'S': total number of samples current approximant relies upon;
- 'sampler': sample point generator.
approx_state: Whether to approximate state.
verbosity: Verbosity level.
POD: Whether to compute POD of snapshots.
scaleFactorDer: Scaling factors for derivative computation.
S: Number of solution snapshots over which current approximant is
based upon.
sampler: Sample point generator.
polybasis: type of polynomial basis for interpolation.
M: Numerator degree of approximant.
N: Denominator degree of approximant.
polydegreetype: Type of polynomial degree.
radialDirectionalWeights: Radial basis weights for interpolant
numerator.
nNearestNeighbor: Number of nearest neighbors considered in numerator
if polybasis allows.
interpRcond: Tolerance for interpolation via numpy.polyfit.
robustTol: Tolerance for robust rational denominator management.
correctorForce: Whether corrector should forcefully delete bad poles.
correctorTol: Tolerance for corrector step.
correctorMaxIter: Maximum number of corrector iterations.
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.
Q: Numpy 1D vector containing complex coefficients of approximant
denominator.
P: Numpy 2D vector whose columns are FE dofs of coefficients of
approximant numerator.
"""
def
_setupInterpolationIndices
(
self
):
"""Setup parameters for polyvander."""
super
()
.
_setupInterpolationIndices
()
if
len
(
self
.
_musUniqueCN
)
>
1
:
raise
RROMPyException
((
"Cannot apply centered-like method with "
"more than one distinct sample point."
))
def
_setupDenominator
(
self
):
"""Compute rational denominator."""
RROMPyAssert
(
self
.
_mode
,
message
=
"Cannot setup denominator."
)
vbMng
(
self
,
"INIT"
,
"Starting computation of denominator."
,
7
)
cfun
=
totalDegreeN
if
self
.
polydegreetype
==
"TOTAL"
else
fullDegreeN
if
hasattr
(
self
,
"_N_isauto"
):
self
.
_setNAuto
()
else
:
N
=
reduceDegreeN
(
self
.
N
,
self
.
S
,
self
.
npar
,
self
.
polydegreetype
)
if
N
<
self
.
N
:
RROMPyWarning
((
"N too large compared to S. Reducing N by "
"{}"
)
.
format
(
self
.
N
-
N
))
self
.
N
=
N
while
self
.
N
>
0
:
invD
,
fitinv
=
self
.
_computeInterpolantInverseBlocks
()
Seff
=
cfun
(
self
.
N
,
self
.
npar
)
idxSamplesEff
=
list
(
range
(
self
.
S
-
Seff
,
self
.
S
))
if
self
.
POD
:
ev
,
eV
=
self
.
findeveVGQR
(
self
.
samplingEngine
.
RPOD
[:,
idxSamplesEff
],
invD
)
else
:
ev
,
eV
=
self
.
findeveVGExplicit
(
self
.
samplingEngine
.
samples
(
idxSamplesEff
),
invD
)
nevBad
=
checkRobustTolerance
(
ev
,
self
.
robustTol
)
if
nevBad
<=
1
:
break
if
self
.
catchInstability
>
0
:
raise
RROMPyException
((
"Instability in denominator "
"computation: eigenproblem is poorly "
"conditioned."
),
self
.
catchInstability
==
1
)
RROMPyWarning
((
"Smallest {} eigenvalues below tolerance. Reducing "
"N by 1."
)
.
format
(
nevBad
))
self
.
N
=
self
.
N
-
1
if
self
.
N
<=
0
:
self
.
N
=
0
eV
=
np
.
ones
((
1
,
1
))
q
=
PI
()
q
.
npar
=
self
.
npar
q
.
polybasis
=
self
.
polybasis0
if
self
.
polydegreetype
==
"TOTAL"
:
q
.
coeffs
=
degreeTotalToFull
(
tuple
([
self
.
N
+
1
]
*
self
.
npar
),
self
.
npar
,
eV
[:,
0
])
else
:
q
.
coeffs
=
eV
[:,
0
]
.
reshape
([
self
.
N
+
1
]
*
self
.
npar
)
vbMng
(
self
,
"DEL"
,
"Done computing denominator."
,
7
)
return
q
,
fitinv
def
_setupNumerator
(
self
):
"""Compute rational numerator."""
RROMPyAssert
(
self
.
_mode
,
message
=
"Cannot setup numerator."
)
vbMng
(
self
,
"INIT"
,
"Starting computation of numerator."
,
7
)
self
.
_setupInterpolationIndices
()
Qevaldiag
=
polyTimesTable
(
self
.
trainedModel
.
data
.
Q
,
self
.
_musUniqueCN
,
self
.
_reorder
,
self
.
_derIdxs
,
self
.
scaleFactorRel
)
if
self
.
POD
:
Qevaldiag
=
Qevaldiag
.
dot
(
self
.
samplingEngine
.
RPOD
.
T
)
if
hasattr
(
self
,
"radialDirectionalWeights"
):
rDW
=
copy
(
self
.
radialDirectionalWeights
)
cfun
=
totalDegreeN
if
self
.
polydegreetype
==
"TOTAL"
else
fullDegreeN
if
hasattr
(
self
,
"_M_isauto"
):
self
.
_setMAuto
()
M
=
self
.
M
else
:
M
=
reduceDegreeN
(
self
.
M
,
self
.
S
,
self
.
npar
,
self
.
polydegreetype
)
if
M
<
self
.
M
:
RROMPyWarning
((
"M too large compared to S. Reducing M by "
"{}"
)
.
format
(
self
.
M
-
M
))
self
.
M
=
M
while
(
self
.
M
>=
0
and
(
not
hasattr
(
self
,
"radialDirectionalWeights"
)
or
self
.
radialDirectionalWeights
[
0
]
<=
rDW
[
0
]
*
2
**
6
)):
Seff
=
cfun
(
self
.
M
,
self
.
npar
)
pParRest
=
[
self
.
verbosity
>=
5
,
self
.
polydegreetype
==
"TOTAL"
,
{
"derIdxs"
:
[
self
.
_derIdxs
[
0
][:
Seff
]],
"reorder"
:
self
.
_reorder
[:
Seff
],
"scl"
:
self
.
scaleFactorRel
}]
if
self
.
polybasis
in
ppb
:
p
=
PI
()
else
:
pParRest
=
[
self
.
radialDirectionalWeights
]
+
pParRest
pParRest
[
-
1
][
"nNearestNeighbor"
]
=
self
.
nNearestNeighbor
p
=
RBI
()
if
self
.
polybasis
in
rbpb
else
MLSI
()
if
self
.
polybasis
in
ppb
+
rbpb
:
pParRest
+=
[{
"rcond"
:
self
.
interpRcond
}]
wellCond
,
msg
=
p
.
setupByInterpolation
(
self
.
_musUniqueCN
,
Qevaldiag
[:
Seff
,
:
Seff
],
self
.
M
,
self
.
polybasis
,
*
pParRest
)
vbMng
(
self
,
"MAIN"
,
msg
,
5
)
if
wellCond
:
break
if
self
.
catchInstability
>
0
:
raise
RROMPyException
((
"Instability in numerator computation: "
"polyfit is poorly conditioned."
),
self
.
catchInstability
==
1
)
if
self
.
polybasis
in
ppb
:
vbMng
(
self
,
"MAIN"
,
(
"Polyfit is poorly conditioned. Reducing "
"M by 1."
),
10
)
self
.
M
=
self
.
M
-
1
else
:
vbMng
(
self
,
"MAIN"
,
(
"Polyfit is poorly conditioned. "
"Multiplying radialDirectionalWeights "
"by 2."
),
10
)
for
j
in
range
(
self
.
npar
):
self
.
_radialDirectionalWeights
[
j
]
*=
2.
if
self
.
M
<
0
or
(
hasattr
(
self
,
"radialDirectionalWeights"
)
and
self
.
radialDirectionalWeights
[
0
]
>
rDW
[
0
]
*
2
**
6
):
raise
RROMPyException
((
"Instability in computation of numerator. "
"Aborting."
))
if
self
.
polybasis
in
ppb
:
self
.
M
=
M
else
:
self
.
radialDirectionalWeights
=
rDW
vbMng
(
self
,
"DEL"
,
"Done computing numerator."
,
7
)
return
p
def
_computeInterpolantInverseBlocks
(
self
)
->
Tuple
[
List
[
Np2D
],
Np2D
]:
"""
Compute inverse factors for minimal interpolant target functional.
"""
RROMPyAssert
(
self
.
_mode
,
message
=
"Cannot solve eigenvalue problem."
)
self
.
_setupInterpolationIndices
()
if
self
.
polydegreetype
==
"TOTAL"
:
cfun
,
TEGen
=
totalDegreeN
,
pvTP
else
:
cfun
,
TEGen
=
fullDegreeN
,
pvP
E
=
max
(
self
.
M
,
self
.
N
)
while
E
>=
0
:
Seff
=
cfun
(
E
,
self
.
npar
)
TEGenPar
=
[
self
.
polybasis0
,
[
self
.
_derIdxs
[
0
][:
Seff
]],
self
.
_reorder
[:
Seff
],
self
.
scaleFactorRel
]
if
self
.
polydegreetype
==
"TOTAL"
:
Eeff
=
E
idxsB
=
totalDegreeMaxMask
(
E
,
self
.
npar
)
else
:
#if self.polydegreetype == "FULL":
Eeff
=
[
E
]
*
self
.
npar
idxsB
=
fullDegreeMaxMask
(
E
,
self
.
npar
)
TE
=
TEGen
(
self
.
_musUniqueCN
,
Eeff
,
*
TEGenPar
)
fitOut
=
customPInv
(
TE
,
rcond
=
self
.
interpRcond
,
full
=
True
)
vbMng
(
self
,
"MAIN"
,
(
"Fitting {} samples with degree {} through {}... "
"Conditioning of pseudoinverse system: {:.4e}."
)
.
format
(
TE
.
shape
[
0
],
E
,
polyfitname
(
self
.
polybasis0
),
fitOut
[
1
][
1
][
0
]
/
fitOut
[
1
][
1
][
-
1
]),
5
)
if
fitOut
[
1
][
0
]
==
TE
.
shape
[
1
]:
fitinv
=
fitOut
[
0
][
idxsB
,
:]
break
if
self
.
catchInstability
>
0
:
raise
RROMPyException
((
"Instability in denominator "
"computation: polyfit is poorly "
"conditioned."
),
self
.
catchInstability
==
1
)
EeqN
=
E
==
self
.
N
vbMng
(
self
,
"MAIN"
,
(
"Polyfit is poorly conditioned. Reducing E {}"
"by 1."
)
.
format
(
"and N "
*
EeqN
),
10
)
if
EeqN
:
self
.
N
=
self
.
N
-
1
E
-=
1
if
self
.
N
<
0
:
raise
RROMPyException
((
"Instability in computation of "
"denominator. Aborting."
))
invD
=
vanderInvTable
(
fitinv
,
idxsB
,
self
.
_reorder
[:
Seff
],
[
self
.
_derIdxs
[
0
][:
Seff
]])
if
self
.
N
==
E
:
TN
=
TE
else
:
if
self
.
polydegreetype
==
"TOTAL"
:
Neff
=
self
.
N
idxsB
=
totalDegreeMaxMask
(
self
.
N
,
self
.
npar
)
else
:
#if self.polydegreetype == "FULL":
Neff
=
[
self
.
N
]
*
self
.
npar
idxsB
=
fullDegreeMaxMask
(
self
.
N
,
self
.
npar
)
TN
=
TEGen
(
self
.
_musUniqueCN
,
Neff
,
*
TEGenPar
)
for
k
in
range
(
len
(
invD
)):
invD
[
k
]
=
dot
(
invD
[
k
],
TN
)
return
invD
,
fitinv
Event Timeline
Log In to Comment