Page Menu
Home
c4science
Search
Configure Global Search
Log In
Files
F61750249
rational_interpolant.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
Wed, May 8, 17:59
Size
37 KB
Mime Type
text/x-python
Expires
Fri, May 10, 17:59 (2 d)
Engine
blob
Format
Raw Data
Handle
17553663
Attached To
R6746 RationalROMPy
rational_interpolant.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
scipy.linalg
import
eigvals
from
collections.abc
import
Iterable
from
rrompy.reduction_methods.base
import
checkRobustTolerance
from
.generic_standard_approximant
import
GenericStandardApproximant
from
rrompy.utilities.poly_fitting.polynomial
import
(
polybases
as
ppb
,
polyfitname
,
polyvander
as
pvP
,
polyTimes
,
polyTimesTable
,
vanderInvTable
,
PolynomialInterpolator
as
PI
,
PolynomialInterpolatorNodal
as
PIN
)
from
rrompy.utilities.poly_fitting.heaviside
import
rational2heaviside
from
rrompy.utilities.poly_fitting.radial_basis
import
(
polybases
as
rbpb
,
RadialBasisInterpolator
as
RBI
)
from
rrompy.utilities.base.types
import
(
Np1D
,
Np2D
,
Tuple
,
List
,
sampList
,
interpEng
)
from
rrompy.utilities.base
import
verbosityManager
as
vbMng
from
rrompy.utilities.numerical
import
pseudoInverse
,
dot
,
potential
from
rrompy.utilities.numerical.hash_derivative
import
nextDerivativeIndices
from
rrompy.utilities.numerical.degree
import
(
reduceDegreeN
,
degreeTotalToFull
,
fullDegreeMaxMask
,
totalDegreeMaxMask
)
from
rrompy.solver
import
Np2DLikeGramian
from
rrompy.utilities.exception_manager
import
(
RROMPyException
,
RROMPyAssert
,
RROMPyWarning
)
__all__
=
[
'RationalInterpolant'
]
class
RationalInterpolant
(
GenericStandardApproximant
):
"""
ROM rational interpolant 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;
- 'radialDirectionalWeightsAdapt': bounds for adaptive rescaling of
radial basis weights; defaults to [-1, -1];
- 'functionalSolve': strategy for minimization of denominator
functional; allowed values include 'NORM', 'DOMINANT', 'NODAL',
'LOEWNER', and 'BARYCENTRIC'; defaults to 'NORM';
- 'interpRcond': tolerance for interpolation; defaults to None;
- 'robustTol': tolerance for robust rational denominator
management; defaults to 0.
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;
- 'radialDirectionalWeightsAdapt': bounds for adaptive rescaling of
radial basis weights;
- 'functionalSolve': strategy for minimization of denominator
functional;
- 'interpRcond': tolerance for interpolation via numpy.polyfit;
- 'robustTol': tolerance for robust rational denominator
management.
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.
radialDirectionalWeightsAdapt: Bounds for adaptive rescaling of radial
basis weights.
functionalSolve: Strategy for minimization of denominator functional.
interpRcond: Tolerance for interpolation via numpy.polyfit.
robustTol: Tolerance for robust rational denominator management.
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
__init__
(
self
,
*
args
,
**
kwargs
):
self
.
_preInit
()
self
.
_addParametersToList
([
"polybasis"
,
"M"
,
"N"
,
"polydegreetype"
,
"radialDirectionalWeights"
,
"radialDirectionalWeightsAdapt"
,
"functionalSolve"
,
"interpRcond"
,
"robustTol"
],
[
"MONOMIAL"
,
"AUTO"
,
"AUTO"
,
"TOTAL"
,
1.
,
[
-
1.
,
-
1.
],
"NORM"
,
-
1
,
0.
])
super
()
.
__init__
(
*
args
,
**
kwargs
)
self
.
catchInstability
=
0
self
.
_postInit
()
@property
def
tModelType
(
self
):
from
.trained_model.trained_model_rational
import
TrainedModelRational
return
TrainedModelRational
@property
def
polybasis
(
self
):
"""Value of polybasis."""
return
self
.
_polybasis
@polybasis.setter
def
polybasis
(
self
,
polybasis
):
try
:
polybasis
=
polybasis
.
upper
()
.
strip
()
.
replace
(
" "
,
""
)
if
polybasis
not
in
ppb
+
rbpb
:
raise
RROMPyException
(
"Prescribed polybasis not recognized."
)
self
.
_polybasis
=
polybasis
except
:
RROMPyWarning
((
"Prescribed polybasis not recognized. Overriding "
"to 'MONOMIAL'."
))
self
.
_polybasis
=
"MONOMIAL"
self
.
_approxParameters
[
"polybasis"
]
=
self
.
polybasis
@property
def
polybasis0
(
self
):
if
"_"
in
self
.
polybasis
:
return
self
.
polybasis
.
split
(
"_"
)[
0
]
return
self
.
polybasis
@property
def
functionalSolve
(
self
):
"""Value of functionalSolve."""
return
self
.
_functionalSolve
@functionalSolve.setter
def
functionalSolve
(
self
,
functionalSolve
):
try
:
functionalSolve
=
functionalSolve
.
upper
()
.
strip
()
.
replace
(
" "
,
""
)
if
functionalSolve
not
in
[
"NORM"
,
"DOMINANT"
,
"NODAL"
,
"LOEWNER"
,
"BARYCENTRIC"
]:
raise
RROMPyException
((
"Prescribed functionalSolve not "
"recognized."
))
self
.
_functionalSolve
=
functionalSolve
except
:
RROMPyWarning
((
"Prescribed functionalSolve not recognized. "
"Overriding to 'NORM'."
))
self
.
_functionalSolve
=
"NORM"
self
.
_approxParameters
[
"functionalSolve"
]
=
self
.
functionalSolve
@property
def
interpRcond
(
self
):
"""Value of interpRcond."""
return
self
.
_interpRcond
@interpRcond.setter
def
interpRcond
(
self
,
interpRcond
):
self
.
_interpRcond
=
interpRcond
self
.
_approxParameters
[
"interpRcond"
]
=
self
.
interpRcond
@property
def
radialDirectionalWeights
(
self
):
"""Value of radialDirectionalWeights."""
return
self
.
_radialDirectionalWeights
@radialDirectionalWeights.setter
def
radialDirectionalWeights
(
self
,
radialDirectionalWeights
):
if
isinstance
(
radialDirectionalWeights
,
Iterable
):
radialDirectionalWeights
=
list
(
radialDirectionalWeights
)
else
:
radialDirectionalWeights
=
[
radialDirectionalWeights
]
self
.
_radialDirectionalWeights
=
radialDirectionalWeights
self
.
_approxParameters
[
"radialDirectionalWeights"
]
=
(
self
.
radialDirectionalWeights
)
@property
def
radialDirectionalWeightsAdapt
(
self
):
"""Value of radialDirectionalWeightsAdapt."""
return
self
.
_radialDirectionalWeightsAdapt
@radialDirectionalWeightsAdapt.setter
def
radialDirectionalWeightsAdapt
(
self
,
radialDirectionalWeightsAdapt
):
self
.
_radialDirectionalWeightsAdapt
=
radialDirectionalWeightsAdapt
self
.
_approxParameters
[
"radialDirectionalWeightsAdapt"
]
=
(
self
.
radialDirectionalWeightsAdapt
)
@property
def
M
(
self
):
"""Value of M."""
return
self
.
_M
@M.setter
def
M
(
self
,
M
):
if
isinstance
(
M
,
str
):
M
=
M
.
strip
()
.
replace
(
" "
,
""
)
if
"-"
not
in
M
:
M
=
M
+
"-0"
self
.
_M_isauto
,
self
.
_M_shift
=
True
,
int
(
M
.
split
(
"-"
)[
-
1
])
M
=
0
if
M
<
0
:
raise
RROMPyException
(
"M must be non-negative."
)
self
.
_M
=
M
self
.
_approxParameters
[
"M"
]
=
self
.
M
def
_setMAuto
(
self
):
self
.
M
=
max
(
0
,
reduceDegreeN
(
self
.
S
,
self
.
S
,
self
.
npar
,
self
.
polydegreetype
)
-
self
.
_M_shift
)
vbMng
(
self
,
"MAIN"
,
"Automatically setting M to {}."
.
format
(
self
.
M
),
25
)
@property
def
N
(
self
):
"""Value of N."""
return
self
.
_N
@N.setter
def
N
(
self
,
N
):
if
isinstance
(
N
,
str
):
N
=
N
.
strip
()
.
replace
(
" "
,
""
)
if
"-"
not
in
N
:
N
=
N
+
"-0"
self
.
_N_isauto
,
self
.
_N_shift
=
True
,
int
(
N
.
split
(
"-"
)[
-
1
])
N
=
0
if
N
<
0
:
raise
RROMPyException
(
"N must be non-negative."
)
self
.
_N
=
N
self
.
_approxParameters
[
"N"
]
=
self
.
N
def
_setNAuto
(
self
):
self
.
N
=
max
(
0
,
reduceDegreeN
(
self
.
S
,
self
.
S
,
self
.
npar
,
self
.
polydegreetype
)
-
self
.
_N_shift
)
vbMng
(
self
,
"MAIN"
,
"Automatically setting N to {}."
.
format
(
self
.
N
),
25
)
@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
.
_musUniqueCN
=
None
self
.
_derIdxs
=
None
self
.
_reorder
=
None
def
_setupInterpolationIndices
(
self
):
"""Setup parameters for polyvander."""
if
self
.
_musUniqueCN
is
None
or
len
(
self
.
_reorder
)
!=
len
(
self
.
mus
):
self
.
_musUniqueCN
,
musIdxsTo
,
musIdxs
,
musCount
=
(
self
.
trainedModel
.
centerNormalize
(
self
.
mus
)
.
unique
(
return_index
=
True
,
return_inverse
=
True
,
return_counts
=
True
))
self
.
_musUnique
=
self
.
mus
[
musIdxsTo
]
self
.
_derIdxs
=
[
None
]
*
len
(
self
.
_musUniqueCN
)
self
.
_reorder
=
np
.
empty
(
len
(
musIdxs
),
dtype
=
int
)
filled
=
0
for
j
,
cnt
in
enumerate
(
musCount
):
self
.
_derIdxs
[
j
]
=
nextDerivativeIndices
([],
self
.
mus
.
shape
[
1
],
cnt
)
jIdx
=
np
.
nonzero
(
musIdxs
==
j
)[
0
]
self
.
_reorder
[
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
)
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
:
if
self
.
functionalSolve
!=
"NORM"
and
self
.
npar
>
1
:
RROMPyWarning
((
"Strategy for functional optimization must be "
"'NORM' for more than one parameter. "
"Overriding to 'NORM'."
))
self
.
functionalSolve
=
"NORM"
if
self
.
functionalSolve
==
"BARYCENTRIC"
and
self
.
N
+
1
<
self
.
S
:
RROMPyWarning
((
"Barycentric strategy cannot be applied with "
"Least Squares. Overriding to 'NORM'."
))
self
.
functionalSolve
=
"NORM"
if
self
.
functionalSolve
==
"BARYCENTRIC"
:
invD
,
TN
=
None
,
None
self
.
_setupInterpolationIndices
()
else
:
invD
,
TN
=
self
.
_computeInterpolantInverseBlocks
()
if
(
self
.
functionalSolve
in
[
"NODAL"
,
"LOEWNER"
,
"BARYCENTRIC"
]
and
len
(
self
.
_musUnique
)
!=
len
(
self
.
mus
)):
if
self
.
functionalSolve
==
"BARYCENTRIC"
:
warnflag
=
"Barycentric"
else
:
warnflag
=
"Iterative"
RROMPyWarning
((
"{} functional optimization cannot be applied "
"to repeated samples. Overriding to "
"'NORM'."
)
.
format
(
warnflag
))
self
.
functionalSolve
=
"NORM"
idxSamplesEff
=
list
(
range
(
self
.
S
))
if
self
.
POD
:
ev
,
eV
=
self
.
findeveVGQR
(
self
.
samplingEngine
.
RPOD
[:,
idxSamplesEff
],
invD
,
TN
)
else
:
ev
,
eV
=
self
.
findeveVGExplicit
(
self
.
samplingEngine
.
samples
(
idxSamplesEff
),
invD
,
TN
)
if
self
.
functionalSolve
in
[
"NODAL"
,
"LOEWNER"
]:
break
nevBad
=
checkRobustTolerance
(
ev
,
self
.
robustTol
)
if
nevBad
<=
(
self
.
functionalSolve
==
"NORM"
):
break
if
self
.
catchInstability
:
raise
RROMPyException
((
"Instability in denominator "
"computation: eigenproblem is poorly "
"conditioned."
),
self
.
catchInstability
==
1
)
vbMng
(
self
,
"MAIN"
,
(
"Smallest {} eigenvalues below tolerance. "
"Reducing N by 1."
)
.
format
(
nevBad
),
10
)
self
.
N
=
self
.
N
-
1
if
self
.
N
<=
0
:
self
.
N
=
0
eV
=
np
.
ones
((
1
,
1
))
if
self
.
N
>
0
and
self
.
functionalSolve
in
[
"NODAL"
,
"LOEWNER"
,
"BARYCENTRIC"
]:
q
=
PIN
()
q
.
polybasis
,
q
.
nodes
=
self
.
polybasis0
,
eV
else
:
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
)
else
:
q
.
coeffs
=
eV
.
reshape
([
self
.
N
+
1
]
*
self
.
npar
)
vbMng
(
self
,
"DEL"
,
"Done computing denominator."
,
7
)
return
q
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
,
"_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
:
pParRest
=
[
self
.
M
,
self
.
polybasis
,
self
.
verbosity
>=
5
,
self
.
polydegreetype
==
"TOTAL"
,
{
"derIdxs"
:
self
.
_derIdxs
,
"reorder"
:
self
.
_reorder
,
"scl"
:
self
.
scaleFactorRel
}]
if
self
.
polybasis
in
ppb
:
p
=
PI
()
else
:
self
.
computeScaleFactor
()
rDWEff
=
np
.
array
([
w
*
f
for
w
,
f
in
zip
(
self
.
radialDirectionalWeights
,
self
.
scaleFactor
)])
pParRest
=
pParRest
[:
2
]
+
[
rDWEff
]
+
pParRest
[
2
:]
pParRest
[
-
1
][
"optimizeScalingBounds"
]
=
(
self
.
radialDirectionalWeightsAdapt
)
p
=
RBI
()
if
self
.
polybasis
in
ppb
+
rbpb
:
pParRest
+=
[{
"rcond"
:
self
.
interpRcond
}]
wellCond
,
msg
=
p
.
setupByInterpolation
(
self
.
_musUniqueCN
,
Qevaldiag
,
*
pParRest
)
vbMng
(
self
,
"MAIN"
,
msg
,
5
)
if
wellCond
:
break
if
self
.
catchInstability
:
raise
RROMPyException
((
"Instability in numerator computation: "
"polyfit is poorly conditioned."
),
self
.
catchInstability
==
1
)
vbMng
(
self
,
"MAIN"
,
(
"Polyfit is poorly conditioned. Reducing M "
"by 1."
),
10
)
self
.
M
=
self
.
M
-
1
if
self
.
M
<
0
:
raise
RROMPyException
((
"Instability in computation of numerator. "
"Aborting."
))
self
.
M
=
M
vbMng
(
self
,
"DEL"
,
"Done computing numerator."
,
7
)
return
p
def
setupApprox
(
self
)
->
int
:
"""Compute rational interpolant."""
if
self
.
checkComputedApprox
():
return
-
1
RROMPyAssert
(
self
.
_mode
,
message
=
"Cannot setup approximant."
)
vbMng
(
self
,
"INIT"
,
"Setting up {}."
.
format
(
self
.
name
()),
5
)
self
.
computeSnapshots
()
self
.
_setupTrainedModel
(
self
.
samplingEngine
.
projectionMatrix
)
self
.
_setupRational
(
self
.
_setupDenominator
())
self
.
trainedModel
.
data
.
approxParameters
=
copy
(
self
.
approxParameters
)
vbMng
(
self
,
"DEL"
,
"Done setting up approximant."
,
5
)
return
0
def
_setupRational
(
self
,
Q
:
interpEng
,
P
:
interpEng
=
None
):
vbMng
(
self
,
"INIT"
,
"Starting approximant finalization."
,
5
)
self
.
trainedModel
.
data
.
Q
=
Q
if
P
is
None
:
P
=
self
.
_setupNumerator
()
if
self
.
N
>
0
and
self
.
npar
==
1
:
pls
=
Q
.
roots
()
idxBad
=
self
.
HFEngine
.
flagBadPolesResidues
(
pls
,
relative
=
True
)
plsN
=
self
.
mapParameterList
(
self
.
mapParameterList
(
self
.
mu0
)(
0
,
0
)
+
self
.
scaleFactor
*
pls
,
"B"
)(
0
)
idxBad
=
np
.
logical_or
(
self
.
HFEngine
.
flagBadPolesResidues
(
pls
,
relative
=
True
),
self
.
HFEngine
.
flagBadPolesResidues
(
plsN
))
if
np
.
any
(
idxBad
):
vbMng
(
self
,
"MAIN"
,
"Removing {} spurious poles out of {} due to poles."
\
.
format
(
np
.
sum
(
idxBad
),
self
.
N
),
10
)
if
isinstance
(
Q
,
PIN
):
Q
.
nodes
=
Q
.
nodes
[
np
.
logical_not
(
idxBad
)]
else
:
Q
=
PI
()
Q
.
npar
=
self
.
npar
Q
.
polybasis
=
self
.
polybasis0
Q
.
coeffs
=
np
.
ones
(
1
,
dtype
=
np
.
complex
)
for
pl
in
pls
[
np
.
logical_not
(
idxBad
)]:
Q
.
coeffs
=
polyTimes
(
Q
.
coeffs
,
[
-
pl
,
1.
],
Pbasis
=
Q
.
polybasis
,
Rbasis
=
Q
.
polybasis
)
Q
.
coeffs
/=
np
.
linalg
.
norm
(
Q
.
coeffs
)
self
.
trainedModel
.
data
.
Q
=
Q
self
.
N
=
Q
.
deg
[
0
]
P
=
self
.
_setupNumerator
()
if
(
not
hasattr
(
self
.
HFEngine
,
"_ignoreResidues"
)
or
not
self
.
HFEngine
.
_ignoreResidues
):
cfs
,
pls
,
_
=
rational2heaviside
(
P
,
Q
)
cfs
=
cfs
[:
self
.
N
]
.
T
if
not
self
.
POD
:
cfs
=
self
.
samplingEngine
.
projectionMatrix
.
dot
(
cfs
)
foci
=
self
.
sampler
.
normalFoci
()
ground
=
self
.
sampler
.
groundPotential
()
potEff
=
potential
(
pls
,
foci
)
/
ground
potEff
[
np
.
logical_or
(
potEff
<
1.
,
np
.
isinf
(
pls
))]
=
1.
cfs
[:,
np
.
isinf
(
pls
)]
=
0.
cfs
/=
potEff
# rescale by potential
idxBad
=
self
.
HFEngine
.
flagBadPolesResidues
(
pls
,
cfs
)
if
np
.
any
(
idxBad
):
vbMng
(
self
,
"MAIN"
,
(
"Removing {} spurious poles out of {} due to "
"residues."
)
.
format
(
np
.
sum
(
idxBad
),
self
.
N
),
10
)
if
isinstance
(
Q
,
PIN
):
Q
.
nodes
=
Q
.
nodes
[
np
.
logical_not
(
idxBad
)]
else
:
Q
=
PI
()
Q
.
npar
=
self
.
npar
Q
.
polybasis
=
self
.
polybasis0
Q
.
coeffs
=
np
.
ones
(
1
,
dtype
=
np
.
complex
)
for
pl
in
pls
[
np
.
logical_not
(
idxBad
)]:
Q
.
coeffs
=
polyTimes
(
Q
.
coeffs
,
[
-
pl
,
1.
],
Pbasis
=
Q
.
polybasis
,
Rbasis
=
Q
.
polybasis
)
Q
.
coeffs
/=
np
.
linalg
.
norm
(
Q
.
coeffs
)
self
.
trainedModel
.
data
.
Q
=
Q
self
.
N
=
Q
.
deg
[
0
]
P
=
self
.
_setupNumerator
()
self
.
trainedModel
.
data
.
P
=
P
vbMng
(
self
,
"DEL"
,
"Terminated approximant finalization."
,
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
.
_setupInterpolationIndices
()
pvPPar
=
[
self
.
polybasis0
,
self
.
_derIdxs
,
self
.
_reorder
,
self
.
scaleFactorRel
]
if
hasattr
(
self
,
"_M_isauto"
):
self
.
_setMAuto
()
E
=
max
(
self
.
M
,
self
.
N
)
while
E
>=
0
:
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
=
pvP
(
self
.
_musUniqueCN
,
Eeff
,
*
pvPPar
)
fitOut
=
pseudoInverse
(
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
:
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
,
self
.
_derIdxs
)
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
=
pvP
(
self
.
_musUniqueCN
,
Neff
,
*
pvPPar
)
return
invD
,
TN
def
findeveVGExplicit
(
self
,
sampleE
:
sampList
,
invD
:
List
[
Np2D
],
TN
:
Np2D
)
->
Tuple
[
Np1D
,
Np2D
]:
"""
Compute explicitly eigenvalues and eigenvectors of rational denominator
matrix.
"""
RROMPyAssert
(
self
.
_mode
,
message
=
"Cannot solve eigenvalue problem."
)
vbMng
(
self
,
"INIT"
,
"Building gramian matrix."
,
10
)
gramian
=
self
.
HFEngine
.
innerProduct
(
sampleE
,
sampleE
,
is_state
=
self
.
approx_state
)
if
self
.
functionalSolve
==
"NODAL"
:
SEnd
=
invD
[
0
]
.
shape
[
1
]
G0
=
np
.
zeros
((
SEnd
,)
*
2
,
dtype
=
np
.
complex
)
elif
self
.
functionalSolve
==
"LOEWNER"
:
G0
=
gramian
if
self
.
functionalSolve
==
"BARYCENTRIC"
:
nEnd
=
len
(
gramian
)
-
1
else
:
nEnd
=
TN
.
shape
[
1
]
G
=
np
.
zeros
((
nEnd
,
nEnd
),
dtype
=
np
.
complex
)
for
k
in
range
(
len
(
invD
)):
iDkN
=
dot
(
invD
[
k
],
TN
)
G
+=
dot
(
dot
(
gramian
,
iDkN
)
.
T
,
iDkN
.
conj
())
.
T
if
self
.
functionalSolve
==
"NODAL"
:
G0
+=
dot
(
dot
(
gramian
,
invD
[
k
])
.
T
,
invD
[
k
]
.
conj
())
.
T
vbMng
(
self
,
"DEL"
,
"Done building gramian."
,
10
)
if
self
.
functionalSolve
==
"NORM"
:
ev
,
eV
=
np
.
linalg
.
eigh
(
G
)
eV
=
eV
[:,
0
]
problem
=
"eigenproblem"
else
:
if
self
.
functionalSolve
==
"BARYCENTRIC"
:
fitOut
=
pseudoInverse
(
gramian
,
rcond
=
self
.
interpRcond
,
full
=
True
)
barWeigths
=
fitOut
[
0
]
.
dot
(
np
.
ones
(
nEnd
+
1
))
eV
=
self
.
findeVBarycentric
(
barWeigths
/
np
.
sum
(
barWeigths
))
else
:
fitOut
=
pseudoInverse
(
G
[:
-
1
,
:
-
1
],
rcond
=
self
.
interpRcond
,
full
=
True
)
eV
=
np
.
append
(
fitOut
[
0
]
.
dot
(
G
[:
-
1
,
-
1
]),
-
1.
)
ev
=
fitOut
[
1
][
1
][::
-
1
]
problem
=
"linear problem"
vbMng
(
self
,
"MAIN"
,
(
"Solved {} of size {} with condition number "
"{:.4e}."
)
.
format
(
problem
,
nEnd
,
ev
[
-
1
]
/
ev
[
0
]),
5
)
if
self
.
functionalSolve
in
[
"NODAL"
,
"LOEWNER"
]:
q
=
PI
()
q
.
npar
,
q
.
polybasis
,
q
.
coeffs
=
self
.
npar
,
self
.
polybasis0
,
eV
eV
,
tol
,
niter
,
passed
=
self
.
findeVNewton
(
q
.
roots
(),
G0
)
if
passed
:
vbMng
(
self
,
"MAIN"
,
(
"Newton algorithm for problem of size {} converged "
"(tol = {:.4e}) in {} iterations."
)
.
format
(
nEnd
,
tol
,
niter
),
5
)
else
:
RROMPyWarning
((
"Newton algorithm for problem of size {} did "
"not converge (tol = {:.4e}) after {} "
"iterations."
)
.
format
(
nEnd
,
tol
,
niter
))
return
ev
,
eV
def
findeveVGQR
(
self
,
RPODE
:
Np2D
,
invD
:
List
[
Np2D
],
TN
:
Np2D
)
->
Tuple
[
Np1D
,
Np2D
]:
"""
Compute eigenvalues and eigenvectors of rational denominator matrix
through SVD of R factor.
"""
RROMPyAssert
(
self
.
_mode
,
message
=
"Cannot solve eigenvalue problem."
)
vbMng
(
self
,
"INIT"
,
"Building half-gramian matrix stack."
,
10
)
if
self
.
functionalSolve
==
"NODAL"
:
gramian
=
Np2DLikeGramian
(
None
,
dot
(
RPODE
,
invD
[
0
]))
elif
self
.
functionalSolve
==
"LOEWNER"
:
gramian
=
Np2DLikeGramian
(
None
,
RPODE
)
if
self
.
functionalSolve
==
"BARYCENTRIC"
:
nEnd
=
RPODE
.
shape
[
1
]
-
1
else
:
S
,
nEnd
,
eWidth
=
RPODE
.
shape
[
0
],
TN
.
shape
[
1
],
len
(
invD
)
Rstack
=
np
.
zeros
((
S
*
eWidth
,
nEnd
),
dtype
=
np
.
complex
)
for
k
in
range
(
eWidth
):
Rstack
[
k
*
S
:
(
k
+
1
)
*
S
,
:]
=
dot
(
RPODE
,
dot
(
invD
[
k
],
TN
))
vbMng
(
self
,
"DEL"
,
"Done building half-gramian."
,
10
)
if
self
.
functionalSolve
in
[
"NORM"
,
"BARYCENTRIC"
]:
if
self
.
functionalSolve
==
"NORM"
:
_
,
s
,
Vh
=
np
.
linalg
.
svd
(
Rstack
,
full_matrices
=
False
)
eV
=
Vh
[
-
1
,
:]
.
conj
()
else
:
#if self.functionalSolve == "BARYCENTRIC":
_
,
s
,
Vh
=
np
.
linalg
.
svd
(
RPODE
,
full_matrices
=
False
)
s
[
np
.
logical_not
(
np
.
isclose
(
s
,
0.
))]
**=
-
2.
barWeigths
=
(
Vh
.
T
.
conj
()
*
s
)
.
dot
(
Vh
.
dot
(
np
.
ones
(
nEnd
+
1
)))
eV
=
self
.
findeVBarycentric
(
barWeigths
/
np
.
sum
(
barWeigths
))
ev
=
s
[::
-
1
]
problem
=
"svd problem"
else
:
fitOut
=
pseudoInverse
(
Rstack
[:,
:
-
1
],
rcond
=
self
.
interpRcond
,
full
=
True
)
ev
=
fitOut
[
1
][
1
][::
-
1
]
eV
=
np
.
append
(
fitOut
[
0
]
.
dot
(
Rstack
[:,
-
1
]),
-
1.
)
problem
=
"linear problem"
vbMng
(
self
,
"MAIN"
,
(
"Solved {} of size {} with condition number "
"{:.4e}."
)
.
format
(
problem
,
nEnd
,
ev
[
-
1
]
/
ev
[
0
]),
5
)
if
self
.
functionalSolve
in
[
"NODAL"
,
"LOEWNER"
]:
q
=
PI
()
q
.
npar
,
q
.
polybasis
,
q
.
coeffs
=
self
.
npar
,
self
.
polybasis0
,
eV
eV
,
tol
,
niter
,
passed
=
self
.
findeVNewton
(
q
.
roots
(),
gramian
)
if
passed
:
vbMng
(
self
,
"MAIN"
,
(
"Newton algorithm for problem of size {} converged "
"(tol = {:.4e}) in {} iterations."
)
.
format
(
nEnd
,
tol
,
niter
),
5
)
else
:
RROMPyWarning
((
"Newton algorithm for problem of size {} did "
"not converge (tol = {:.2e}) after {} "
"iterations."
)
.
format
(
nEnd
,
tol
,
niter
))
return
ev
,
eV
def
findeVBarycentric
(
self
,
baryWeights
:
Np1D
)
->
Np1D
:
RROMPyAssert
(
self
.
_mode
,
message
=
"Cannot solve optimization problem."
)
arrow
=
np
.
pad
(
np
.
diag
(
self
.
_musUniqueCN
.
data
[
self
.
_reorder
]
.
flatten
()),
(
1
,
0
),
"constant"
,
constant_values
=
1.
)
+
0.j
arrow
[
0
,
0
]
=
0.
arrow
[
0
,
1
:]
=
baryWeights
active
=
np
.
pad
(
np
.
eye
(
len
(
baryWeights
)),
(
1
,
0
),
"constant"
)
eV
=
eigvals
(
arrow
,
active
)
return
eV
[
np
.
logical_not
(
np
.
isinf
(
eV
))]
def
findeVNewton
(
self
,
nodes0
:
Np1D
,
gram
:
Np2D
,
maxiter
:
int
=
25
,
tolNewton
:
float
=
1e-10
)
\
->
Tuple
[
Np1D
,
float
,
int
,
bool
]:
RROMPyAssert
(
self
.
_mode
,
message
=
"Cannot solve optimization problem."
)
algo
=
self
.
functionalSolve
N
=
len
(
nodes0
)
nodes
=
nodes0
grad
=
np
.
zeros
(
N
,
dtype
=
np
.
complex
)
hess
=
np
.
zeros
((
N
,
N
),
dtype
=
np
.
complex
)
mu
=
np
.
repeat
(
self
.
_musUniqueCN
.
data
[
self
.
_reorder
],
N
,
axis
=
1
)
for
niter
in
range
(
maxiter
):
if
algo
==
"NODAL"
:
plDist
=
mu
-
nodes
.
reshape
(
1
,
-
1
)
q0
,
q
=
np
.
prod
(
plDist
,
axis
=
1
),
[]
elif
algo
==
"LOEWNER"
:
loew
=
np
.
pad
(
np
.
power
(
mu
-
nodes
.
reshape
(
1
,
-
1
),
-
1.
),
[(
0
,
0
),
(
1
,
0
)],
'constant'
,
constant_values
=
1.
)
loewI
=
pseudoInverse
(
loew
)
Ids
=
[]
for
jS
in
range
(
N
):
if
algo
==
"NODAL"
:
mask
=
np
.
arange
(
N
)
!=
jS
q
+=
[
np
.
prod
(
plDist
[:,
mask
],
axis
=
1
)]
grad
[
jS
]
=
q
[
-
1
]
.
conj
()
.
dot
(
gram
.
dot
(
q0
))
elif
algo
==
"LOEWNER"
:
Ids
+=
[
loewI
.
dot
(
np
.
power
(
loew
[:,
1
+
jS
],
2.
))]
zIj
,
jI
=
Ids
[
-
1
][
0
],
loewI
[
1
+
jS
]
grad
[
jS
]
=
(
zIj
*
jI
)
.
conj
()
.
dot
(
gram
.
dot
(
loewI
[
0
]))
for
iS
in
range
(
jS
+
1
):
if
algo
==
"NODAL"
:
if
iS
==
jS
:
hij
=
0.
sij
=
q
[
-
1
]
.
conj
()
.
dot
(
gram
.
dot
(
q
[
-
1
]))
else
:
mask
=
np
.
logical_and
(
np
.
arange
(
N
)
!=
iS
,
np
.
arange
(
N
)
!=
jS
)
qij
=
np
.
prod
(
plDist
[:,
mask
],
axis
=
1
)
hij
=
qij
.
conj
()
.
dot
(
gram
.
dot
(
q0
))
sij
=
q
[
-
1
]
.
conj
()
.
dot
(
gram
.
dot
(
q
[
iS
]))
elif
algo
==
"LOEWNER"
:
zIi
,
iIj
=
Ids
[
iS
][
0
],
Ids
[
-
1
][
1
+
iS
]
hij
=
(
zIi
*
iIj
*
jI
)
.
conj
()
.
dot
(
gram
.
dot
(
loewI
[
0
]))
if
iS
==
jS
:
iI
=
jI
zIdd
=
loewI
[
0
]
.
dot
(
np
.
power
(
loew
[:,
1
+
jS
],
3.
))
hij
+=
(
zIdd
*
jI
)
.
conj
()
.
dot
(
gram
.
dot
(
loewI
[
0
]))
hij
*=
2.
else
:
jIi
,
iI
=
Ids
[
iS
][
1
+
jS
],
loewI
[
1
+
iS
]
hij
+=
(
zIj
*
jIi
*
iI
)
.
conj
()
.
dot
(
gram
.
dot
(
loewI
[
0
]))
sij
=
(
zIj
*
jI
)
.
conj
()
.
dot
(
gram
.
dot
(
zIi
*
iI
))
hess
[
jS
,
iS
]
=
hij
+
sij
if
iS
<
jS
:
hess
[
iS
,
jS
]
=
hij
+
sij
.
conj
()
dnodes
=
np
.
linalg
.
solve
(
hess
,
grad
)
nodes
+=
dnodes
tol
=
np
.
linalg
.
norm
(
dnodes
)
/
np
.
linalg
.
norm
(
nodes
)
if
tol
<
tolNewton
:
break
return
nodes
,
tol
,
niter
,
niter
<
maxiter
def
getResidues
(
self
,
*
args
,
**
kwargs
)
->
Np2D
:
"""
Obtain approximant residues.
Returns:
Matrix with residues as columns.
"""
return
self
.
trainedModel
.
getResidues
(
*
args
,
**
kwargs
)
Event Timeline
Log In to Comment