Page Menu
Home
c4science
Search
Configure Global Search
Log In
Files
F62610591
rational_interpolant_pivoted.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
Tue, May 14, 08:13
Size
31 KB
Mime Type
text/x-python
Expires
Thu, May 16, 08:13 (2 d)
Engine
blob
Format
Raw Data
Handle
17668022
Attached To
R6746 RationalROMPy
rational_interpolant_pivoted.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
.generic_pivoted_approximant
import
GenericPivotedApproximant
from
rrompy.reduction_methods.standard.rational_interpolant
import
(
RationalInterpolant
as
RI
)
from
rrompy.utilities.poly_fitting.polynomial
import
(
polybases
as
ppb
,
polyfitname
,
polyvander
as
pvP
,
polyvanderTotal
as
pvTP
,
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
(
polybases
as
mlspb
,
MovingLeastSquaresInterpolator
as
MLSI
)
from
rrompy.utilities.base.types
import
(
Np1D
,
Np2D
,
HFEng
,
DictAny
,
Tuple
,
List
,
ListAny
,
paramVal
)
from
rrompy.utilities.base
import
verbosityManager
as
vbMng
,
freepar
as
fp
from
rrompy.utilities.numerical
import
(
multifactorial
,
customPInv
,
dot
,
fullDegreeN
,
totalDegreeN
,
degreeTotalToFull
,
fullDegreeMaxMask
,
totalDegreeMaxMask
,
nextDerivativeIndices
,
hashDerivativeToIdx
as
hashD
,
hashIdxToDerivative
as
hashI
)
from
rrompy.utilities.exception_manager
import
(
RROMPyException
,
RROMPyAssert
,
RROMPyWarning
)
from
rrompy.parameter
import
checkParameter
__all__
=
[
'RationalInterpolantPivoted'
]
class
RationalInterpolantPivoted
(
GenericPivotedApproximant
):
"""
ROM pivoted rational interpolant (with pole matching) computation for
parametric problems.
Args:
HFEngine: HF problem solver.
mu0(optional): Default parameter. Defaults to 0.
directionPivot(optional): Pivot components. 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;
- 'matchingWeight': weight for pole matching optimization; defaults
to 1;
- 'cutOffTolerance': tolerance for ignoring parasitic poles;
defaults to np.inf;
- 'cutOffType': rule for tolerance computation for parasitic poles;
defaults to 'MAGNITUDE';
- 'S': total number of pivot samples current approximant relies
upon;
- 'samplerPivot': pivot sample point generator;
- 'SMarginal': total number of marginal samples current approximant
relies upon;
- 'samplerMarginal': marginal sample point generator;
- 'polybasisPivot': type of polynomial basis for pivot
interpolation; defaults to 'MONOMIAL';
- 'polybasisMarginal': type of polynomial basis for marginal
interpolation; allowed values include 'MONOMIAL', 'CHEBYSHEV'
and 'LEGENDRE'; defaults to 'MONOMIAL';
- 'M': degree of rational interpolant numerator; defaults to 0;
- 'N': degree of rational interpolant denominator; defaults to 0;
- 'polydegreetype': type of polynomial degree; defaults to 'TOTAL';
- 'MMarginal': degree of marginal interpolant; defaults to 0;
- 'polydegreetypeMarginal': type of polynomial degree for marginal;
defaults to 'TOTAL';
- 'radialDirectionalWeightsPivot': radial basis weights for pivot
numerator; defaults to 0, i.e. identity;
- 'radialDirectionalWeightsMarginal': radial basis weights for
marginal interpolant; defaults to 0, i.e. identity;
- 'nNearestNeighborPivot': number of pivot nearest neighbors
considered if polybasisPivot allows; defaults to -1;
- 'nNearestNeighborMarginal': number of marginal nearest neighbors
considered if polybasisMarginal allows; defaults to -1;
- 'interpRcondPivot': tolerance for pivot interpolation; defaults
to None;
- 'interpRcondMarginal': tolerance for marginal interpolation;
defaults to None;
- 'robustTol': tolerance for robust rational denominator
management; defaults to 0.
Defaults to empty dict.
force_state(optional): Whether to approximate state. Defaults to False.
verbosity(optional): Verbosity level. Defaults to 10.
Attributes:
HFEngine: HF problem solver.
mu0: Default parameter.
directionPivot: Pivot components.
mus: Array of snapshot parameters.
musPivot: Array of pivot snapshot parameters.
musMarginal: Array of marginal 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;
- 'matchingWeight': weight for pole matching optimization;
- 'cutOffTolerance': tolerance for ignoring parasitic poles;
- 'cutOffType': rule for tolerance computation for parasitic poles;
- 'polybasisPivot': type of polynomial basis for pivot
interpolation;
- 'polybasisMarginal': type of polynomial basis for marginal
interpolation;
- 'M': degree of rational interpolant numerator;
- 'N': degree of rational interpolant denominator;
- 'polydegreetype': type of polynomial degree;
- 'MMarginal': degree of marginal interpolant;
- 'polydegreetypeMarginal': type of polynomial degree for marginal;
- 'radialDirectionalWeightsPivot': radial basis weights for pivot
numerator;
- 'radialDirectionalWeightsMarginal': radial basis weights for
marginal interpolant;
- 'nNearestNeighborPivot': number of pivot nearest neighbors
considered if polybasisPivot allows;
- 'nNearestNeighborMarginal': number of marginal nearest neighbors
considered if polybasisMarginal allows;
- 'interpRcondPivot': tolerance for pivot interpolation;
- 'interpRcondMarginal': tolerance for marginal interpolation;
- 'robustTol': tolerance for robust rational denominator
management.
parameterListCritical: Recognized keys of critical approximant
parameters:
- 'S': total number of pivot samples current approximant relies
upon;
- 'samplerPivot': pivot sample point generator;
- 'SMarginal': total number of marginal samples current approximant
relies upon;
- 'samplerMarginal': marginal sample point generator.
force_state: Whether to approximate state.
verbosity: Verbosity level.
POD: Whether to compute POD of snapshots.
matchingWeight: Weight for pole matching optimization.
cutOffTolerance: Tolerance for ignoring parasitic poles.
cutOffType: Rule for tolerance computation for parasitic poles.
S: Total number of pivot samples current approximant relies upon.
sampler: Pivot sample point generator.
polybasisPivot: Type of polynomial basis for pivot interpolation.
polybasisMarginal: Type of polynomial basis for marginal interpolation.
M: Numerator degree of approximant.
N: Denominator degree of approximant.
polydegreetype: Type of polynomial degree.
MMarginal: Degree of marginal interpolant.
polydegreetypeMarginal: Type of polynomial degree for marginal.
radialDirectionalWeightsPivot: Radial basis weights for pivot
numerator.
radialDirectionalWeightsMarginal: Radial basis weights for marginal
interpolant.
nNearestNeighborPivot: Number of pivot nearest neighbors considered if
polybasisPivot allows.
nNearestNeighborMarginal: Number of marginal nearest neighbors
considered if polybasisMarginal allows.
interpRcondPivot: Tolerance for pivot interpolation.
interpRcondMarginal: Tolerance for marginal interpolation.
robustTol: Tolerance for robust rational denominator management.
muBoundsPivot: list of bounds for pivot parameter values.
muBoundsMarginal: list of bounds for marginal 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
__init__
(
self
,
HFEngine
:
HFEng
,
mu0
:
paramVal
=
None
,
directionPivot
:
ListAny
=
[
0
],
approxParameters
:
DictAny
=
{},
force_state
:
bool
=
False
,
verbosity
:
int
=
10
,
timestamp
:
bool
=
True
):
self
.
_preInit
()
self
.
_addParametersToList
([
"polybasisPivot"
,
"M"
,
"N"
,
"polydegreetype"
,
"radialDirectionalWeightsPivot"
,
"nNearestNeighborPivot"
,
"interpRcondPivot"
,
"robustTol"
],
[
"MONOMIAL"
,
0
,
0
,
"TOTAL"
,
1
,
-
1
,
-
1
,
0
])
super
()
.
__init__
(
HFEngine
=
HFEngine
,
mu0
=
mu0
,
directionPivot
=
directionPivot
,
approxParameters
=
approxParameters
,
force_state
=
force_state
,
verbosity
=
verbosity
,
timestamp
=
timestamp
)
self
.
_postInit
()
@property
def
tModelType
(
self
):
from
rrompy.reduction_methods.trained_model
import
\
TrainedModelPivotedRational
return
TrainedModelPivotedRational
def
initializeModelData
(
self
,
datadict
):
from
rrompy.reduction_methods.trained_model
import
\
TrainedModelPivotedData
return
(
TrainedModelPivotedData
(
datadict
[
"mu0"
],
datadict
.
pop
(
"projMat"
),
datadict
[
"scaleFactor"
],
datadict
.
pop
(
"rescalingExp"
),
datadict
[
"directionPivot"
]),
[
"mu0"
,
"scaleFactor"
,
"directionPivot"
,
"mus"
])
@property
def
polybasisPivot
(
self
):
"""Value of polybasisPivot."""
return
self
.
_polybasisPivot
@polybasisPivot.setter
def
polybasisPivot
(
self
,
polybasisPivot
):
try
:
polybasisPivot
=
polybasisPivot
.
upper
()
.
strip
()
.
replace
(
" "
,
""
)
if
polybasisPivot
not
in
ppb
+
rbpb
+
mlspb
:
raise
RROMPyException
(
"Prescribed pivot polybasis not recognized."
)
self
.
_polybasisPivot
=
polybasisPivot
except
:
RROMPyWarning
((
"Prescribed pivot polybasis not recognized. "
"Overriding to 'MONOMIAL'."
))
self
.
_polybasisPivot
=
"MONOMIAL"
self
.
_approxParameters
[
"polybasisPivot"
]
=
self
.
polybasisPivot
@property
def
polybasisPivot0
(
self
):
if
"_"
in
self
.
polybasisPivot
:
return
self
.
polybasisPivot
.
split
(
"_"
)[
0
]
return
self
.
polybasisPivot
@property
def
radialDirectionalWeightsPivot
(
self
):
"""Value of radialDirectionalWeightsPivot."""
return
self
.
_radialDirectionalWeightsPivot
@radialDirectionalWeightsPivot.setter
def
radialDirectionalWeightsPivot
(
self
,
radialDirectionalWeightsPivot
):
self
.
_radialDirectionalWeightsPivot
=
radialDirectionalWeightsPivot
self
.
_approxParameters
[
"radialDirectionalWeightsPivot"
]
=
(
self
.
radialDirectionalWeightsPivot
)
@property
def
nNearestNeighborPivot
(
self
):
"""Value of nNearestNeighborPivot."""
return
self
.
_nNearestNeighborPivot
@nNearestNeighborPivot.setter
def
nNearestNeighborPivot
(
self
,
nNearestNeighborPivot
):
self
.
_nNearestNeighborPivot
=
nNearestNeighborPivot
self
.
_approxParameters
[
"nNearestNeighborPivot"
]
=
(
self
.
nNearestNeighborPivot
)
@property
def
interpRcondPivot
(
self
):
"""Value of interpRcondPivot."""
return
self
.
_interpRcondPivot
@interpRcondPivot.setter
def
interpRcondPivot
(
self
,
interpRcondPivot
):
self
.
_interpRcondPivot
=
interpRcondPivot
self
.
_approxParameters
[
"interpRcondPivot"
]
=
self
.
interpRcondPivot
@property
def
M
(
self
):
"""Value of M."""
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
@property
def
N
(
self
):
"""Value of N."""
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
@property
def
polydegreetype
(
self
):
"""Value of polydegreetype."""
return
self
.
_polydegreetype
@polydegreetype.setter
def
polydegreetype
(
self
,
polydegreetype
):
try
:
polydegreetype
=
polydegreetype
.
upper
()
.
strip
()
.
replace
(
" "
,
""
)
if
polydegreetype
not
in
[
"TOTAL"
,
"FULL"
]:
raise
RROMPyException
((
"Prescribed polydegreetype not "
"recognized."
))
self
.
_polydegreetype
=
polydegreetype
except
:
RROMPyWarning
((
"Prescribed polydegreetype not recognized. "
"Overriding to 'TOTAL'."
))
self
.
_polydegreetype
=
"TOTAL"
self
.
_approxParameters
[
"polydegreetype"
]
=
self
.
polydegreetype
@property
def
robustTol
(
self
):
"""Value of tolerance for robust rational 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
def
resetSamples
(
self
):
"""Reset samples."""
super
()
.
resetSamples
()
self
.
_musPUniqueCN
=
None
self
.
_derPIdxs
=
None
self
.
_reorderP
=
None
def
_setupPivotInterpolationIndices
(
self
):
"""Setup parameters for polyvander."""
RROMPyAssert
(
self
.
_mode
,
message
=
"Cannot setup interpolation indices."
)
if
(
self
.
_musPUniqueCN
is
None
or
len
(
self
.
_reorderP
)
!=
len
(
self
.
musPivot
)):
self
.
_musPUniqueCN
,
musPIdxsTo
,
musPIdxs
,
musPCount
=
(
self
.
trainedModel
.
centerNormalizePivot
(
self
.
musPivot
)
.
unique
(
return_index
=
True
,
return_inverse
=
True
,
return_counts
=
True
))
self
.
_musPUnique
=
self
.
mus
[
musPIdxsTo
]
self
.
_derPIdxs
=
[
None
]
*
len
(
self
.
_musPUniqueCN
)
self
.
_reorderP
=
np
.
empty
(
len
(
musPIdxs
),
dtype
=
int
)
filled
=
0
for
j
,
cnt
in
enumerate
(
musPCount
):
self
.
_derPIdxs
[
j
]
=
nextDerivativeIndices
([],
self
.
nparPivot
,
cnt
)
jIdx
=
np
.
nonzero
(
musPIdxs
==
j
)[
0
]
self
.
_reorderP
[
jIdx
]
=
np
.
arange
(
filled
,
filled
+
cnt
)
filled
+=
cnt
def
_setupDenominator
(
self
):
"""Compute rational denominator."""
RROMPyAssert
(
self
.
_mode
,
message
=
"Cannot setup denominator."
)
vbMng
(
self
,
"INIT"
,
"Starting computation of denominator."
,
7
)
NinvD
=
None
N0
=
copy
(
self
.
N
)
qs
=
[]
self
.
verbosity
-=
10
for
j
in
range
(
len
(
self
.
musMarginal
)):
self
.
_N
=
N0
while
self
.
N
>
0
:
if
NinvD
!=
self
.
N
:
invD
,
fitinvP
=
self
.
_computeInterpolantInverseBlocks
()
NinvD
=
self
.
N
if
self
.
POD
:
ev
,
eV
=
RI
.
findeveVGQR
(
self
,
self
.
samplingEngine
.
RPOD
[
j
],
invD
)
else
:
ev
,
eV
=
RI
.
findeveVGExplicit
(
self
,
self
.
samplingEngine
.
samples
[
j
],
invD
)
nevBad
=
checkRobustTolerance
(
ev
,
self
.
robustTol
)
if
nevBad
<=
1
:
break
if
self
.
catchInstability
:
raise
RROMPyException
((
"Instability in denominator "
"computation: eigenproblem is "
"poorly conditioned."
))
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
.
nparPivot
q
.
polybasis
=
self
.
polybasisPivot0
if
self
.
polydegreetype
==
"TOTAL"
:
q
.
coeffs
=
degreeTotalToFull
(
tuple
([
self
.
N
+
1
]
*
q
.
npar
),
q
.
npar
,
eV
[:,
0
])
else
:
q
.
coeffs
=
eV
[:,
0
]
.
reshape
([
self
.
N
+
1
]
*
q
.
npar
)
qs
=
qs
+
[
copy
(
q
)]
self
.
verbosity
+=
10
vbMng
(
self
,
"DEL"
,
"Done computing denominator."
,
7
)
return
qs
,
fitinvP
def
_setupNumerator
(
self
):
"""Compute rational numerator."""
RROMPyAssert
(
self
.
_mode
,
message
=
"Cannot setup numerator."
)
vbMng
(
self
,
"INIT"
,
"Starting computation of numerator."
,
7
)
Qevaldiag
=
np
.
zeros
((
len
(
self
.
musPivot
),
len
(
self
.
musPivot
)),
dtype
=
np
.
complex
)
verb
=
self
.
trainedModel
.
verbosity
self
.
trainedModel
.
verbosity
=
0
self
.
_setupPivotInterpolationIndices
()
cfun
=
totalDegreeN
if
self
.
polydegreetype
==
"TOTAL"
else
fullDegreeN
M
=
copy
(
self
.
M
)
while
len
(
self
.
musPivot
)
<
cfun
(
M
,
self
.
nparPivot
):
M
-=
1
if
M
<
self
.
M
:
RROMPyWarning
((
"M too large compared to S. Reducing M by "
"{}"
)
.
format
(
self
.
M
-
M
))
self
.
M
=
M
tensor_idx
=
0
ps
=
[]
for
k
,
muM
in
enumerate
(
self
.
musMarginal
):
self
.
_M
=
M
idxGlob
=
0
for
j
,
derIdxs
in
enumerate
(
self
.
_derPIdxs
):
mujEff
=
[
fp
]
*
self
.
npar
for
jj
,
kk
in
enumerate
(
self
.
directionPivot
):
mujEff
[
kk
]
=
self
.
_musPUnique
[
j
,
jj
]
for
jj
,
kk
in
enumerate
(
self
.
directionMarginal
):
mujEff
[
kk
]
=
muM
(
0
,
jj
)
mujEff
=
checkParameter
(
mujEff
,
self
.
npar
)
nder
=
len
(
derIdxs
)
idxLoc
=
np
.
arange
(
len
(
self
.
musPivot
))[
(
self
.
_reorderP
>=
idxGlob
)
*
(
self
.
_reorderP
<
idxGlob
+
nder
)]
idxGlob
+=
nder
Qval
=
[
0
]
*
nder
for
der
in
range
(
nder
):
derIdx
=
hashI
(
der
,
self
.
nparPivot
)
derIdxEff
=
[
0
]
*
self
.
npar
sclEff
=
[
0
]
*
self
.
npar
for
jj
,
kk
in
enumerate
(
self
.
directionPivot
):
derIdxEff
[
kk
]
=
derIdx
[
jj
]
sclEff
[
kk
]
=
self
.
scaleFactorPivot
[
jj
]
**
-
1.
Qval
[
der
]
=
(
self
.
trainedModel
.
getQVal
(
mujEff
,
derIdxEff
,
scl
=
sclEff
)
/
multifactorial
(
derIdx
))
for
derU
,
derUIdx
in
enumerate
(
derIdxs
):
for
derQ
,
derQIdx
in
enumerate
(
derIdxs
):
diffIdx
=
[
x
-
y
for
(
x
,
y
)
in
zip
(
derUIdx
,
derQIdx
)]
if
all
([
x
>=
0
for
x
in
diffIdx
]):
diffj
=
hashD
(
diffIdx
)
Qevaldiag
[
idxLoc
[
derU
],
idxLoc
[
derQ
]]
=
Qval
[
diffj
]
while
self
.
M
>=
0
:
if
self
.
polybasisPivot
in
ppb
:
p
=
PI
()
wellCond
,
msg
=
p
.
setupByInterpolation
(
self
.
_musPUniqueCN
,
Qevaldiag
,
self
.
M
,
self
.
polybasisPivot
,
self
.
verbosity
>=
5
,
self
.
polydegreetype
==
"TOTAL"
,
{
"derIdxs"
:
self
.
_derPIdxs
,
"reorder"
:
self
.
_reorderP
,
"scl"
:
np
.
power
(
self
.
scaleFactorPivot
,
-
1.
)},
{
"rcond"
:
self
.
interpRcondPivot
})
elif
self
.
polybasisPivot
in
rbpb
:
p
=
RBI
()
wellCond
,
msg
=
p
.
setupByInterpolation
(
self
.
_musPUniqueCN
,
Qevaldiag
,
self
.
M
,
self
.
polybasisPivot
,
self
.
radialDirectionalWeightsPivot
,
self
.
verbosity
>=
5
,
self
.
polydegreetype
==
"TOTAL"
,
{
"derIdxs"
:
self
.
_derPIdxs
,
"reorder"
:
self
.
_reorderP
,
"scl"
:
np
.
power
(
self
.
scaleFactorPivot
,
-
1.
),
"nNearestNeighbor"
:
self
.
nNearestNeighborPivot
},
{
"rcond"
:
self
.
interpRcondPivot
})
else
:
# if self.polybasisPivot in mlspb:
p
=
MLSI
()
wellCond
,
msg
=
p
.
setupByInterpolation
(
self
.
_musPUniqueCN
,
Qevaldiag
,
self
.
M
,
self
.
polybasisPivot
,
self
.
radialDirectionalWeightsPivot
,
self
.
verbosity
>=
5
,
self
.
polydegreetype
==
"TOTAL"
,
{
"derIdxs"
:
self
.
_derPIdxs
,
"reorder"
:
self
.
_reorderP
,
"scl"
:
np
.
power
(
self
.
scaleFactorPivot
,
-
1.
),
"nNearestNeighbor"
:
self
.
nNearestNeighborPivot
})
vbMng
(
self
,
"MAIN"
,
msg
,
5
)
if
wellCond
:
break
if
self
.
catchInstability
:
raise
RROMPyException
((
"Instability in numerator "
"computation: polyfit is "
"poorly conditioned."
))
RROMPyWarning
((
"Polyfit is poorly conditioned. "
"Reducing M by 1."
))
self
.
M
=
self
.
M
-
1
tensor_idx_new
=
tensor_idx
+
Qevaldiag
.
shape
[
1
]
if
self
.
POD
:
p
.
postmultiplyTensorize
(
self
.
samplingEngine
.
RPODCoalesced
.
T
[
tensor_idx
:
tensor_idx_new
,
:])
else
:
p
.
pad
(
tensor_idx
,
len
(
self
.
mus
)
-
tensor_idx_new
)
tensor_idx
=
tensor_idx_new
ps
=
ps
+
[
copy
(
p
)]
self
.
trainedModel
.
verbosity
=
verb
vbMng
(
self
,
"DEL"
,
"Done computing numerator."
,
7
)
return
ps
def
setupApprox
(
self
):
"""
Compute rational interpolant.
SVD-based robust eigenvalue management.
"""
if
self
.
checkComputedApprox
():
return
RROMPyAssert
(
self
.
_mode
,
message
=
"Cannot setup approximant."
)
vbMng
(
self
,
"INIT"
,
"Setting up {}."
.
format
(
self
.
name
()),
5
)
self
.
computeSnapshots
()
pMat
=
self
.
samplingEngine
.
samplesCoalesced
.
data
pMatEff
=
dot
(
self
.
HFEngine
.
C
,
pMat
)
if
self
.
force_state
else
pMat
if
self
.
trainedModel
is
None
:
self
.
trainedModel
=
self
.
tModelType
()
self
.
trainedModel
.
verbosity
=
self
.
verbosity
self
.
trainedModel
.
timestamp
=
self
.
timestamp
datadict
=
{
"mu0"
:
self
.
mu0
,
"projMat"
:
pMatEff
,
"scaleFactor"
:
self
.
scaleFactor
,
"rescalingExp"
:
self
.
HFEngine
.
rescalingExp
,
"directionPivot"
:
self
.
directionPivot
}
self
.
trainedModel
.
data
=
self
.
initializeModelData
(
datadict
)[
0
]
else
:
self
.
trainedModel
=
self
.
trainedModel
self
.
trainedModel
.
data
.
projMat
=
copy
(
pMatEff
)
self
.
trainedModel
.
data
.
musPivot
=
copy
(
self
.
musPivot
)
self
.
trainedModel
.
data
.
musMarginal
=
copy
(
self
.
musMarginal
)
self
.
trainedModel
.
data
.
marginalInterp
=
self
.
_setupMarginalInterp
()
if
self
.
N
>
0
:
Qs
=
self
.
_setupDenominator
()[
0
]
else
:
Q
=
PI
()
Q
.
npar
=
self
.
nparPivot
Q
.
coeffs
=
np
.
ones
(
tuple
([
1
]
*
Q
.
npar
),
dtype
=
self
.
musMarginal
.
dtype
)
Q
.
polybasis
=
self
.
polybasisPivot0
Qs
=
[
Q
for
_
in
range
(
len
(
self
.
musMarginal
))]
self
.
trainedModel
.
data
.
_temporary
=
1
self
.
trainedModel
.
data
.
Qs
=
Qs
self
.
trainedModel
.
data
.
Ps
=
self
.
_setupNumerator
()
vbMng
(
self
,
"INIT"
,
"Matching poles."
,
10
)
self
.
trainedModel
.
initializeFromRational
(
self
.
HFEngine
,
self
.
matchingWeight
,
self
.
POD
,
self
.
force_state
)
del
self
.
trainedModel
.
data
.
_temporary
vbMng
(
self
,
"DEL"
,
"Done matching poles."
,
10
)
if
not
np
.
isinf
(
self
.
cutOffTolerance
):
vbMng
(
self
,
"INIT"
,
"Recompressing by cut-off."
,
10
)
msg
=
self
.
trainedModel
.
recompressByCutOff
([
-
1.
,
1.
],
self
.
cutOffTolerance
,
self
.
cutOffType
)
vbMng
(
self
,
"DEL"
,
"Done recompressing."
+
msg
,
10
)
self
.
trainedModel
.
data
.
approxParameters
=
copy
(
self
.
approxParameters
)
vbMng
(
self
,
"DEL"
,
"Done setting up approximant."
,
5
)
def
_computeInterpolantInverseBlocks
(
self
)
->
Tuple
[
List
[
Np2D
],
Np2D
]:
"""
Compute inverse factors for minimal interpolant target functional.
"""
RROMPyAssert
(
self
.
_mode
,
message
=
"Cannot solve eigenvalue problem."
)
self
.
_setupPivotInterpolationIndices
()
cfun
=
totalDegreeN
if
self
.
polydegreetype
==
"TOTAL"
else
fullDegreeN
N
=
copy
(
self
.
N
)
while
len
(
self
.
musPivot
)
<
cfun
(
N
,
self
.
nparPivot
):
N
-=
1
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
:
if
self
.
polydegreetype
==
"TOTAL"
:
TE
,
_
,
argIdxs
=
pvTP
(
self
.
_musPUniqueCN
,
self
.
N
,
self
.
polybasisPivot0
,
self
.
_derPIdxs
,
self
.
_reorderP
,
scl
=
np
.
power
(
self
.
scaleFactorPivot
,
-
1.
))
TE
=
TE
[:,
argIdxs
]
idxsB
=
totalDegreeMaxMask
(
self
.
N
,
self
.
nparPivot
)
else
:
#if self.polydegreetype == "FULL":
TE
=
pvP
(
self
.
_musPUniqueCN
,
[
self
.
N
]
*
self
.
nparPivot
,
self
.
polybasisPivot0
,
self
.
_derPIdxs
,
self
.
_reorderP
,
scl
=
np
.
power
(
self
.
scaleFactorPivot
,
-
1.
))
idxsB
=
fullDegreeMaxMask
(
self
.
N
,
self
.
nparPivot
)
fitOut
=
customPInv
(
TE
,
rcond
=
self
.
interpRcondPivot
,
full
=
True
)
vbMng
(
self
,
"MAIN"
,
(
"Fitting {} samples with degree {} through {}... "
"Conditioning of pseudoinverse system: {:.4e}."
)
.
format
(
TE
.
shape
[
0
],
self
.
N
,
polyfitname
(
self
.
polybasisPivot0
),
fitOut
[
1
][
1
][
0
]
/
fitOut
[
1
][
1
][
-
1
]),
5
)
if
fitOut
[
1
][
0
]
==
TE
.
shape
[
1
]:
fitinvP
=
fitOut
[
0
][
idxsB
,
:]
break
RROMPyWarning
(
"Polyfit is poorly conditioned. Reducing N by 1."
)
self
.
N
-=
1
if
self
.
N
<
0
:
raise
RROMPyException
((
"Instability in computation of "
"denominator. Aborting."
))
TN
,
_
,
argIdxs
=
pvTP
(
self
.
_musPUniqueCN
,
self
.
N
,
self
.
polybasisPivot0
,
self
.
_derPIdxs
,
self
.
_reorderP
,
scl
=
np
.
power
(
self
.
scaleFactorPivot
,
-
1.
))
TN
=
TN
[:,
argIdxs
]
invD
=
[
None
]
*
(
len
(
idxsB
))
for
k
in
range
(
len
(
idxsB
)):
pseudoInv
=
np
.
diag
(
fitinvP
[
k
,
:])
idxGlob
=
0
for
j
,
derIdxs
in
enumerate
(
self
.
_derPIdxs
):
nder
=
len
(
derIdxs
)
idxGlob
+=
nder
if
nder
>
1
:
idxLoc
=
np
.
arange
(
len
(
self
.
musPivot
))[
(
self
.
_reorderP
>=
idxGlob
-
nder
)
*
(
self
.
_reorderP
<
idxGlob
)]
invLoc
=
fitinvP
[
k
,
idxLoc
]
pseudoInv
[
np
.
ix_
(
idxLoc
,
idxLoc
)]
=
0.
for
diffj
,
diffjIdx
in
enumerate
(
derIdxs
):
for
derQ
,
derQIdx
in
enumerate
(
derIdxs
):
derUIdx
=
[
x
-
y
for
(
x
,
y
)
in
zip
(
diffjIdx
,
derQIdx
)]
if
all
([
x
>=
0
for
x
in
derUIdx
]):
derU
=
hashD
(
derUIdx
)
pseudoInv
[
idxLoc
[
derU
],
idxLoc
[
derQ
]]
=
(
invLoc
[
diffj
])
invD
[
k
]
=
dot
(
pseudoInv
,
TN
)
return
invD
,
fitinvP
def
getResidues
(
self
,
*
args
,
**
kwargs
)
->
Np1D
:
"""
Obtain approximant residues.
Returns:
Matrix with residues as columns.
"""
return
self
.
trainedModel
.
getResidues
(
*
args
,
**
kwargs
)
Event Timeline
Log In to Comment