diff --git a/examples/2d/pod/fracture_pod.py b/examples/2d/pod/fracture_pod.py
index 99c6baf..36c04e2 100644
--- a/examples/2d/pod/fracture_pod.py
+++ b/examples/2d/pod/fracture_pod.py
@@ -1,275 +1,263 @@
import numpy as np
from rrompy.hfengines.linear_problem.bidimensional import \
MembraneFractureEngine as MFE
from rrompy.reduction_methods.standard import RationalInterpolant as RI
from rrompy.reduction_methods.standard import ReducedBasis as RB
from rrompy.reduction_methods.pivoted import RationalInterpolantPivoted as RIP
from rrompy.reduction_methods.pivoted import ReducedBasisPivoted as RBP
-from rrompy.reduction_methods.pole_matching import \
- RationalInterpolantPoleMatching as RIPM
-from rrompy.reduction_methods.pole_matching import \
- ReducedBasisPoleMatching as RBPM
from rrompy.parameter.parameter_sampling import (QuadratureSampler as QS,
QuadratureSamplerTotal as QST,
ManualSampler as MS,
RandomSampler as RS)
verb = 70
size = 2
show_sample = True
show_norm = True
Delta = 0
MN = 5
R = (MN + 2) * (MN + 1) // 2
S = [int(np.ceil(R ** .5))] * 2
PODTol = 1e-6
SPivot = [MN + 1, 4]
MMarginal = SPivot[1] - 1
matchingWeight = 1.
cutOffTolerance = 1. * np.inf
cutOffType = "MAGNITUDE"
samples = "centered"
samples = "standard"
samples = "pivoted"
-samples = "pole matching"
algo = "rational"
algo = "RB"
sampling = "quadrature"
sampling = "quadrature_total"
#sampling = "random"
samplingM = "quadrature"
#samplingM = "quadrature_total"
#samplingM = "random"
polydegreetype = "TOTAL"
polydegreetype = "FULL"
-if samples in ["standard", "pivoted", "pole matching"]:
+if samples in ["standard", "pivoted"]:
radial = ""
# radial = "_gaussian"
# radial = "_thinplate"
# radial = "_multiquadric"
rW0 = 5.
radialWeight = [rW0]
-if samples in ["pivoted", "pole matching"]:
+if samples == "pivoted":
radialM = ""
# radialM = "_gaussian"
# radialM = "_thinplate"
# radialM = "_multiquadric"
rMW0 = 2.
radialWeightM = [rMW0]
assert Delta <= 0
if size == 1: # below
mu0 = [40 ** .5, .4]
mutar = [45 ** .5, .4]
murange = [[30 ** .5, .3], [50 ** .5, .5]]
elif size == 2: # top
mu0 = [40 ** .5, .6]
mutar = [45 ** .5, .6]
murange = [[30 ** .5, .5], [50 ** .5, .7]]
elif size == 3: # interesting
mu0 = [40 ** .5, .5]
mutar = [45 ** .5, .5]
murange = [[30 ** .5, .3], [50 ** .5, .7]]
elif size == 4: # wide_low
mu0 = [40 ** .5, .2]
mutar = [45 ** .5, .2]
murange = [[10 ** .5, .1], [70 ** .5, .3]]
elif size == 5: # wide_hi
mu0 = [40 ** .5, .8]
mutar = [45 ** .5, .8]
murange = [[10 ** .5, .7], [70 ** .5, .9]]
elif size == 6: # top_zoom
mu0 = [50 ** .5, .8]
mutar = [55 ** .5, .8]
murange = [[40 ** .5, .7], [60 ** .5, .9]]
elif size == 7: # huge
mu0 = [50 ** .5, .5]
mutar = [55 ** .5, .5]
murange = [[10 ** .5, .2], [90 ** .5, .8]]
elif size == 11: #L below
mu0 = [110 ** .5, .4]
mutar = [115 ** .5, .4]
murange = [[90 ** .5, .3], [130 ** .5, .5]]
elif size == 12: #L top
mu0 = [110 ** .5, .6]
mutar = [115 ** .5, .6]
murange = [[90 ** .5, .5], [130 ** .5, .7]]
elif size == 13: #L interesting
mu0 = [110 ** .5, .5]
mutar = [115 ** .5, .5]
murange = [[90 ** .5, .3], [130 ** .5, .7]]
elif size == 14: #L belowbelow
mu0 = [110 ** .5, .2]
mutar = [115 ** .5, .2]
murange = [[90 ** .5, .1], [130 ** .5, .3]]
elif size == 15: #L toptop
mu0 = [110 ** .5, .8]
mutar = [115 ** .5, .8]
murange = [[90 ** .5, .7], [130 ** .5, .9]]
elif size == 16: #L interestinginteresting
mu0 = [110 ** .5, .5]
mutar = [115 ** .5, .6]
murange = [[90 ** .5, .1], [130 ** .5, .9]]
elif size == 17: #L interestingtop
mu0 = [110 ** .5, .7]
mutar = [115 ** .5, .6]
murange = [[90 ** .5, .5], [130 ** .5, .9]]
elif size == 18: #L interestingbelow
mu0 = [110 ** .5, .3]
mutar = [115 ** .5, .4]
murange = [[90 ** .5, .1], [130 ** .5, .5]]
elif size == 100: # tiny
mu0 = [32.5 ** .5, .5]
mutar = [34 ** .5, .5]
murange = [[30 ** .5, .3], [35 ** .5, .7]]
aEff = 1.#25
bEff = 1. - aEff
murangeEff = [[(aEff*murange[0][0]**2. + bEff*murange[1][0]**2.) ** .5,
aEff*murange[0][1] + bEff*murange[1][1]],
[(aEff*murange[1][0]**2. + bEff*murange[0][0]**2.) ** .5,
aEff*murange[1][1] + bEff*murange[0][1]]]
H = 1.
L = .75
delta = .05
n = 20
solver = MFE(mu0 = mu0, H = H, L = L, delta = delta, n = n, verbosity = verb)
rescalingExp = [2., 1.]
if algo == "rational":
params = {'N':MN, 'M':MN + Delta, 'S':S, 'POD':True,
'polydegreetype':polydegreetype}
if samples == "standard":
params['polybasis'] = "CHEBYSHEV" + radial
# params['polybasis'] = "LEGENDRE" + radial
# params['polybasis'] = "MONOMIAL" + radial
params['radialDirectionalWeights'] = radialWeight
method = RI
elif samples == "centered":
params['polybasis'] = "MONOMIAL"
params['S'] = R
method = RI
- elif samples in ["pivoted", "pole matching"]:
+ elif samples == "pivoted":
params['S'] = [SPivot[0]]
params['SMarginal'] = [SPivot[1]]
params['MMarginal'] = MMarginal
params['polybasisPivot'] = "CHEBYSHEV" + radial
params['polybasisMarginal'] = "MONOMIAL" + radialM
params['radialDirectionalWeightsPivot'] = radialWeight
params['radialDirectionalWeightsMarginal'] = radialWeightM
- if samples == "pivoted":
- method = RIP
- else:
- params['matchingWeight'] = matchingWeight
- params['cutOffTolerance'] = cutOffTolerance
- params["cutOffType"] = cutOffType
- method = RIPM
+ params['matchingWeight'] = matchingWeight
+ params['cutOffTolerance'] = cutOffTolerance
+ params["cutOffType"] = cutOffType
+ method = RIP
else: #if algo == "RB":
params = {'R':(MN + 2 + Delta) * (MN + 1 + Delta) // 2, 'S':S,
'POD':True, 'PODTolerance':PODTol}
if samples == "standard":
method = RB
elif samples == "centered":
params['S'] = R
method = RB
- elif samples in ["pivoted", "pole matching"]:
+ elif samples == "pivoted":
params['S'] = [SPivot[0]]
params['R'] = SPivot[0]
params['SMarginal'] = [SPivot[1]]
params['MMarginal'] = MMarginal
params['polybasisMarginal'] = "MONOMIAL" + radialM
params['radialDirectionalWeightsMarginal'] = radialWeightM
+ params['matchingWeight'] = matchingWeight
+ params['cutOffTolerance'] = cutOffTolerance
+ params["cutOffType"] = cutOffType
method = RBP
- if samples == "pivoted":
- method = RBP
- else:
- params['matchingWeight'] = matchingWeight
- params['cutOffTolerance'] = cutOffTolerance
- params["cutOffType"] = cutOffType
- method = RBPM
if samples == "standard":
if sampling == "quadrature":
params['sampler'] = QS(murange, "CHEBYSHEV", scalingExp = rescalingExp)
# params['sampler'] = QS(murange, "GAUSSLEGENDRE", scalingExp = rescalingExp)
# params['sampler'] = QS(murange, "UNIFORM", scalingExp = rescalingExp)
params['S'] = [max(j, MN + 1) for j in params['S']]
elif sampling == "quadrature_total":
params['sampler'] = QST(murange, "CHEBYSHEV", scalingExp = rescalingExp)
params['S'] = R
else: # if sampling == "random":
params['sampler'] = RS(murange, "HALTON", scalingExp = rescalingExp)
params['S'] = R
elif samples == "centered":
params['sampler'] = MS(murange, points = [mu0], scalingExp = rescalingExp)
-elif samples in ["pivoted", "pole matching"]:
+elif samples == "pivoted":
if sampling == "quadrature":
params['samplerPivot'] = QS([murange[0][0], murange[1][0]], "CHEBYSHEV")
# params['samplerPivot'] = QS([murange[0][0], murange[1][0]], "GAUSSLEGENDRE")
# params['samplerPivot'] = QS([murange[0][0], murange[1][0]], "UNIFORM")
elif sampling == "quadrature_total":
params['samplerPivot'] = QST([murange[0][0], murange[1][0]], "CHEBYSHEV")
else: # if sampling == "random":
params['samplerPivot'] = RS([murange[0][0], murange[1][0]], "HALTON")
if samplingM == "quadrature":
params['samplerMarginal'] = QS([murange[0][1], murange[1][1]], "UNIFORM")
elif samplingM == "quadrature_total":
params['samplerMarginal'] = QST([murange[0][1], murange[1][1]], "CHEBYSHEV")
else: # if samplingM == "random":
params['samplerMarginal'] = RS([murange[0][1], murange[1][1]], "HALTON")
-if samples not in ["pivoted", "pole matching"]:
+if samples != "pivoted":
approx = method(solver, mu0 = mu0, approxParameters = params,
verbosity = verb)
else:
approx = method(solver, mu0 = mu0, directionPivot = [0],
approxParameters = params, verbosity = verb)
if samples != "centered": approx.samplingEngine.allowRepeatedSamples = False
approx.setupApprox()
if show_sample:
import ufl
import fenics as fen
from rrompy.solver.fenics import affine_warping
L = mutar[1]
y = fen.SpatialCoordinate(solver.V.mesh())[1]
warp1, warpI1 = affine_warping(solver.V.mesh(),
np.array([[1, 0], [0, 2. * L]]))
warp2, warpI2 = affine_warping(solver.V.mesh(),
np.array([[1, 0], [0, 2. - 2. * L]]))
warp = ufl.conditional(ufl.ge(y, 0.), warp1, warp2)
warpI = ufl.conditional(ufl.ge(y, 0.), warpI1, warpI2)
approx.plotApprox(mutar, [warp, warpI], name = 'u_app', what = "REAL")
approx.plotHF(mutar, [warp, warpI], name = 'u_HF', what = "REAL")
approx.plotErr(mutar, [warp, warpI], name = 'err', what = "REAL")
# approx.plotRes(mutar, [warp, warpI], name = 'res', what = "REAL")
appErr = approx.normErr(mutar)
solNorm = approx.normHF(mutar)
resNorm = approx.normRes(mutar)
RHSNorm = approx.normRHS(mutar)
print(('SolNorm:\t{}\nErr:\t{}\nErrRel:\t{}').format(solNorm, appErr,
np.divide(appErr, solNorm)))
print(('RHSNorm:\t{}\nRes:\t{}\nResRel:\t{}').format(RHSNorm, resNorm,
np.divide(resNorm, RHSNorm)))
approx.verbosity = 5
from plot_zero_set import plotZeroSet2
muZeroVals, Qvals = plotZeroSet2(murange, murangeEff, approx, mu0,
200, [2., 1.])
if show_norm:
solver._solveBatchSize = 100
from plot_inf_set import plotInfSet2
muInfVals, normEx, normApp, normRes, normErr, beta = plotInfSet2(
murange, murangeEff, approx, mu0, 50,
[2., 1.], relative = True,
nobeta = True)
print(np.sort(approx.getPoles([None, .5]) ** 2.))
diff --git a/examples/2d/pod/matrix_passive_pod.py b/examples/2d/pod/matrix_passive_pod.py
index 8a5fa82..74e68f7 100644
--- a/examples/2d/pod/matrix_passive_pod.py
+++ b/examples/2d/pod/matrix_passive_pod.py
@@ -1,203 +1,192 @@
import numpy as np
from rrompy.hfengines.linear_problem.tridimensional import \
MatrixDynamicalPassive as MDP
from rrompy.reduction_methods.standard import RationalInterpolant as RI
from rrompy.reduction_methods.standard import ReducedBasis as RB
from rrompy.reduction_methods.pivoted import RationalInterpolantPivoted as RIP
from rrompy.reduction_methods.pivoted import ReducedBasisPivoted as RBP
-from rrompy.reduction_methods.pole_matching import \
- RationalInterpolantPoleMatching as RIPM
-from rrompy.reduction_methods.pole_matching import \
- ReducedBasisPoleMatching as RBPM
from rrompy.parameter.parameter_sampling import (QuadratureSampler as QS,
QuadratureSamplerTotal as QST,
ManualSampler as MS,
RandomSampler as RS)
verb = 10
size = 3
show_sample = True
show_norm = True
Delta = 0
MN = 7
R = (MN + 2) * (MN + 1) // 2
S = [int(np.ceil(R ** .5))] * 2
PODTol = 1e-6
SPivot = [MN + 1, 3]
MMarginal = SPivot[1] - 1
samples = "centered"
samples = "standard"
samples = "pivoted"
-samples = "pole matching"
algo = "rational"
#algo = "RB"
sampling = "quadrature"
#sampling = "quadrature_total"
#sampling = "random"
samplingM = "quadrature"
#samplingM = "quadrature_total"
#samplingM = "random"
-if samples in ["standard", "pivoted", "pole matching"]:
+if samples in ["standard", "pivoted"]:
radial = ""
# radial = "_gaussian"
# radial = "_thinplate"
# radial = "_multiquadric"
rW0 = 10.
radialWeight = [rW0]
-if samples in ["pivoted", "pole matching"]:
+if samples == "pivoted":
radialM = ""
# radialM = "_gaussian"
# radialM = "_thinplate"
# radialM = "_multiquadric"
rMW0 = 2.
radialWeightM = [rMW0]
matchingWeight = 1.
cutOffTolerance = 5.
cutOffType = "POTENTIAL"
if size == 1:
mu0 = [2.7e2, 20]
mutar = [3e2, 25]
murange = [[20., 10], [5.2e2, 30]]
elif size == 2:
mu0 = [2.7e2, 60]
mutar = [3e2, 75]
murange = [[20., 10], [5.2e2, 110]]
elif size == 3:
mu0 = [2.7e2, 160]
mutar = [3e2, 105]
murange = [[20., 10], [5.2e2, 310]]
assert Delta <= 0
aEff = 1.#25
bEff = 1. - aEff
murangeEff = [[aEff*murange[0][0] + bEff*murange[1][0],
aEff*murange[0][1] + bEff*murange[1][1]],
[aEff*murange[1][0] + bEff*murange[0][0],
aEff*murange[1][1] + bEff*murange[0][1]]]
n = 100
b = 10
solver = MDP(mu0 = mu0, n = n, b = b, verbosity = verb)
if algo == "rational":
params = {'N':MN, 'M':MN + Delta, 'S':S, 'POD':True}
if samples == "standard":
params['polybasis'] = "CHEBYSHEV" + radial
# params['polybasis'] = "LEGENDRE" + radial
# params['polybasis'] = "MONOMIAL" + radial
params['radialDirectionalWeights'] = radialWeight
method = RI
elif samples == "centered":
params['polybasis'] = "MONOMIAL"
params['S'] = R
method = RI
- elif samples in ["pivoted", "pole matching"]:
+ elif samples == "pivoted":
params['S'] = [SPivot[0]]
params['SMarginal'] = [SPivot[1]]
params['MMarginal'] = MMarginal
params['polybasisPivot'] = "CHEBYSHEV" + radial
params['polybasisMarginal'] = "MONOMIAL" + radialM
params['radialDirectionalWeightsPivot'] = radialWeight
params['radialDirectionalWeightsMarginal'] = radialWeightM
- if samples == "pivoted":
- method = RIP
- else:
- params['matchingWeight'] = matchingWeight
- params['cutOffTolerance'] = cutOffTolerance
- params["cutOffType"] = cutOffType
- method = RIPM
+ params['matchingWeight'] = matchingWeight
+ params['cutOffTolerance'] = cutOffTolerance
+ params["cutOffType"] = cutOffType
+ method = RIP
else: #if algo == "RB":
params = {'R':(MN + 2 + Delta) * (MN + 1 + Delta) // 2, 'S':S,
'POD':True, 'PODTolerance':PODTol}
if samples == "standard":
method = RB
elif samples == "centered":
params['S'] = R
method = RB
- elif samples in ["pivoted", "pole matching"]:
+ elif samples == "pivoted":
params['R'] = SPivot[0]
params['S'] = [SPivot[0]]
params['SMarginal'] = [SPivot[1]]
params['MMarginal'] = MMarginal
params['polybasisMarginal'] = "MONOMIAL" + radialM
params['radialDirectionalWeightsMarginal'] = radialWeightM
- if samples == "pivoted":
- method = RBP
- else:
- params['matchingWeight'] = matchingWeight
- params['cutOffTolerance'] = cutOffTolerance
- params["cutOffType"] = cutOffType
- method = RBPM
+ params['matchingWeight'] = matchingWeight
+ params['cutOffTolerance'] = cutOffTolerance
+ params["cutOffType"] = cutOffType
+ method = RBP
if samples == "standard":
if sampling == "quadrature":
params['sampler'] = QS(murange, "CHEBYSHEV")
# params['sampler'] = QS(murange, "GAUSSLEGENDRE")
# params['sampler'] = QS(murange, "UNIFORM")
elif sampling == "quadrature_total":
params['sampler'] = QST(murange, "CHEBYSHEV")
else: # if sampling == "random":
params['sampler'] = RS(murange, "HALTON")
elif samples == "centered":
params['sampler'] = MS(murange, points = [mu0])
-elif samples in ["pivoted", "pole matching"]:
+elif samples == "pivoted":
if sampling == "quadrature":
params['samplerPivot'] = QS([murange[0][0], murange[1][0]], "CHEBYSHEV")
# params['samplerPivot'] = QS([murange[0][0], murange[1][0]], "GAUSSLEGENDRE")
# params['samplerPivot'] = QS([murange[0][0], murange[1][0]], "UNIFORM")
elif sampling == "quadrature_total":
params['samplerPivot'] = QST([murange[0][0], murange[1][0]], "CHEBYSHEV")
else: # if sampling == "random":
params['samplerPivot'] = RS([murange[0][0], murange[1][0]], "HALTON")
if samplingM == "quadrature":
params['samplerMarginal'] = QS([murange[0][1], murange[1][1]], "UNIFORM")
elif samplingM == "quadrature_total":
params['samplerMarginal'] = QST([murange[0][1], murange[1][1]], "CHEBYSHEV")
else: # if samplingM == "random":
params['samplerMarginal'] = RS([murange[0][1], murange[1][1]], "HALTON")
-if samples not in ["pivoted", "pole matching"]:
+if samples != "pivoted":
approx = method(solver, mu0 = mu0, approxParameters = params,
verbosity = verb)
else:
approx = method(solver, mu0 = mu0, directionPivot = [0],
approxParameters = params, verbosity = verb)
if samples != "centered": approx.samplingEngine.allowRepeatedSamples = False
approx.setupApprox()
if show_sample:
L = mutar[1]
approx.plotApprox(mutar, name = 'u_app', what = "REAL")
approx.plotHF(mutar, name = 'u_HF', what = "REAL")
approx.plotErr(mutar, name = 'err', what = "REAL")
# approx.plotRes(mutar, name = 'res', what = "REAL")
appErr = approx.normErr(mutar)
solNorm = approx.normHF(mutar)
resNorm = approx.normRes(mutar)
RHSNorm = approx.normRHS(mutar)
print(('SolNorm:\t{}\nErr:\t{}\nErrRel:\t{}').format(solNorm, appErr,
np.divide(appErr, solNorm)))
print(('RHSNorm:\t{}\nRes:\t{}\nResRel:\t{}').format(RHSNorm, resNorm,
np.divide(resNorm, RHSNorm)))
try:
from plot_zero_set import plotZeroSet2
muZeroVals, Qvals = plotZeroSet2(murange, murangeEff, approx, mu0,
200, [1., 1], polesImTol = 2.)
except: pass
if show_norm:
solver._solveBatchSize = 100
from plot_inf_set import plotInfSet2
muInfVals, normEx, normApp, normRes, normErr, beta = plotInfSet2(
murange, murangeEff, approx, mu0, 50,
[1., 1.], relative = False,
nobeta = True)
print(1.j * approx.getPoles([None, 50.]))
diff --git a/examples/2d/pod/square_pod.py b/examples/2d/pod/square_pod.py
index 990ddf5..77a00fe 100644
--- a/examples/2d/pod/square_pod.py
+++ b/examples/2d/pod/square_pod.py
@@ -1,191 +1,189 @@
import numpy as np
from rrompy.hfengines.linear_problem.bidimensional import \
HelmholtzSquareDomainProblemEngine as HSDPE
from rrompy.reduction_methods.standard import RationalInterpolant as RI
from rrompy.reduction_methods.standard import ReducedBasis as RB
-from rrompy.reduction_methods.pole_matching import \
- RationalInterpolantPoleMatching as RIPM
-from rrompy.reduction_methods.pole_matching import \
- ReducedBasisPoleMatching as RBPM
+from rrompy.reduction_methods.pivoted import RationalInterpolantPivoted as RIP
+from rrompy.reduction_methods.pivoted import ReducedBasisPivoted as RBP
from rrompy.parameter.parameter_sampling import (QuadratureSampler as QS,
QuadratureSamplerTotal as QST,
ManualSampler as MS,
RandomSampler as RS)
verb = 5
size = 7
show_sample = False
show_norm = True
ignore_forcing = True
ignore_forcing = False
clip = -1
#clip = .4
#clip = .6
MN = 6
R = (MN + 2) * (MN + 1) // 2
S = [int(np.ceil(R ** .5))] * 2
SPivot = [MN + 1, 3]
MMarginal = SPivot[1] - 2
matchingWeight = 10.
samples = "centered"
samples = "standard"
-samples = "pole matching"
+samples = "pivoted"
algo = "rational"
#algo = "RB"
sampling = "quadrature"
sampling = "quadrature_total"
#sampling = "random"
samplingM = "quadrature"
#samplingM = "quadrature_total"
#samplingM = "random"
-if samples == "pole matching":
+if samples == "pivoted":
radialM = ""
radialM = "_gaussian"
# radialM = "_thinplate"
# radialM = "_multiquadric"
rMW0 = 2.
radialWeightM = [rMW0]
if size == 1: # small
mu0 = [4 ** .5, 1.5 ** .5]
mutar = [5 ** .5, 1.75 ** .5]
murange = [[2 ** .5, 1. ** .5], [6 ** .5, 2. ** .5]]
elif size == 2: # medium
mu0 = [4 ** .5, 1.75 ** .5]
mutar = [5 ** .5, 1.25 ** .5]
murange = [[1 ** .5, 1. ** .5], [7 ** .5, 2.5 ** .5]]
elif size == 3: # fat
mu0 = [6 ** .5, 4 ** .5]
mutar = [2 ** .5, 2.5 ** .5]
murange = [[0 ** .5, 2 ** .5], [12 ** .5, 6 ** .5]]
elif size == 4: # crowded
mu0 = [10 ** .5, 2 ** .5]
mutar = [9 ** .5, 2.25 ** .5]
murange = [[8 ** .5, 1.5 ** .5], [12 ** .5, 2.5 ** .5]]
elif size == 5: # tall
mu0 = [11 ** .5, 2.25 ** .5]
mutar = [10.5 ** .5, 2.5 ** .5]
murange = [[10 ** .5, 1.5 ** .5], [12 ** .5, 3 ** .5]]
elif size == 6: # taller
mu0 = [11 ** .5, 2.25 ** .5]
mutar = [10.5 ** .5, 2.5 ** .5]
murange = [[10 ** .5, 1.25 ** .5], [12 ** .5, 3.25 ** .5]]
elif size == 7: # low
mu0 = [7 ** .5, .75 ** .5]
mutar = [6.5 ** .5, .9 ** .5]
murange = [[6 ** .5, .5 ** .5], [8 ** .5, 1. ** .5]]
aEff = 1.1
bEff = 1. - aEff
murangeEff = [[(aEff*murange[0][0]**2. + bEff*murange[1][0]**2.) ** .5,
(aEff*murange[0][1]**2. + bEff*murange[1][1]**2.) ** .5],
[(aEff*murange[1][0]**2. + bEff*murange[0][0]**2.) ** .5,
(aEff*murange[1][1]**2. + bEff*murange[0][1]**2.) ** .5]]
solver = HSDPE(kappa = 2.5, theta = np.pi / 3, mu0 = mu0, n = 20,
verbosity = verb)
if ignore_forcing: solver.nbs = 1
rescalingExp = [2., 1.]
if algo == "rational":
params = {'N':MN, 'M':MN, 'S':S, 'POD':True}
if samples == "standard":
params['polybasis'] = "CHEBYSHEV"
# params['polybasis'] = "LEGENDRE"
# params['polybasis'] = "MONOMIAL"
method = RI
elif samples == "centered":
params['polybasis'] = "MONOMIAL"
params['S'] = R
method = RI
- elif samples == "pole matching":
+ elif samples == "pivoted":
params['S'] = [SPivot[0]]
params['SMarginal'] = [SPivot[1]]
params['MMarginal'] = MMarginal
params['polybasisPivot'] = "CHEBYSHEV"
params['polybasisMarginal'] = "MONOMIAL" + radialM
params['matchingWeight'] = matchingWeight
params['radialDirectionalWeightsMarginal'] = radialWeightM
- method = RIPM
+ method = RIP
else: #if algo == "RB":
params = {'R':R, 'S':S, 'POD':True}
if samples == "standard":
method = RB
elif samples == "centered":
params['S'] = R
method = RB
- elif samples == "pole matching":
+ elif samples == "pivoted":
params['S'] = [SPivot[0]]
params['SMarginal'] = [SPivot[1]]
params['MMarginal'] = MMarginal
params['polybasisMarginal'] = "MONOMIAL" + radialM
params['matchingWeight'] = matchingWeight
params['radialDirectionalWeightsMarginal'] = radialWeightM
- method = RBPM
+ method = RBP
if samples == "standard":
if sampling == "quadrature":
params['sampler'] = QS(murange, "CHEBYSHEV", scalingExp = rescalingExp)
# params['sampler'] = QS(murange, "GAUSSLEGENDRE", scalingExp = rescalingExp)
# params['sampler'] = QS(murange, "UNIFORM", scalingExp = rescalingExp)
params['S'] = [max(j, MN + 1) for j in params['S']]
elif sampling == "quadrature_total":
params['sampler'] = QST(murange, "CHEBYSHEV", scalingExp = rescalingExp)
params['S'] = R
else: # if sampling == "random":
params['sampler'] = RS(murange, "HALTON", scalingExp = rescalingExp)
params['S'] = R
elif samples == "centered":
params['sampler'] = MS(murange, points = [mu0], scalingExp = rescalingExp)
-elif samples == "pole matching":
+elif samples == "pivoted":
if sampling == "quadrature":
params['samplerPivot'] = QS([murange[0][0], murange[1][0]], "CHEBYSHEV")
# params['samplerPivot'] = QS([murange[0][0], murange[1][0]], "GAUSSLEGENDRE")
# params['samplerPivot'] = QS([murange[0][0], murange[1][0]], "UNIFORM")
elif sampling == "quadrature_total":
params['samplerPivot'] = QST([murange[0][0], murange[1][0]], "CHEBYSHEV")
else: # if sampling == "random":
params['samplerPivot'] = RS([murange[0][0], murange[1][0]], "HALTON")
if samplingM == "quadrature":
params['samplerMarginal'] = QS([murange[0][1], murange[1][1]], "UNIFORM")
elif samplingM == "quadrature_total":
params['samplerMarginal'] = QST([murange[0][1], murange[1][1]], "CHEBYSHEV")
else: # if samplingM == "random":
params['samplerMarginal'] = RS([murange[0][1], murange[1][1]], "HALTON")
-if samples != "pole matching":
+if samples != "pivoted":
approx = method(solver, mu0 = mu0, approxParameters = params,
verbosity = verb)
else:
approx = method(solver, mu0 = mu0, directionPivot = [0],
approxParameters = params, verbosity = verb)
if samples != "centered": approx.samplingEngine.allowRepeatedSamples = False
approx.setupApprox()
if show_sample:
approx.plotApprox(mutar, name = 'u_app')
approx.plotHF(mutar, name = 'u_HF')
approx.plotErr(mutar, name = 'err')
approx.plotRes(mutar, name = 'res')
appErr, solNorm = approx.normErr(mutar), approx.normHF(mutar)
resNorm, RHSNorm = approx.normRes(mutar), approx.normRHS(mutar)
print(('SolNorm:\t{}\nErr:\t{}\nErrRel:\t{}').format(solNorm, appErr,
np.divide(appErr, solNorm)))
print(('RHSNorm:\t{}\nRes:\t{}\nResRel:\t{}').format(RHSNorm, resNorm,
np.divide(resNorm, RHSNorm)))
if algo == "rational":
from plot_zero_set import plotZeroSet2
muZeroVals, Qvals = plotZeroSet2(murange, murangeEff, approx, mu0,
200, [2., 1.])
if show_norm:
solver._solveBatchSize = 100
from plot_inf_set import plotInfSet2
muInfVals, normEx, normApp, normRes, normErr, beta = plotInfSet2(
murange, murangeEff, approx, mu0, 20, [2., 1.], clip = clip,
relative = False, nobeta = True)
diff --git a/examples/3d/pod/fracture3_pod.py b/examples/3d/pod/fracture3_pod.py
index 1bc1751..be59a22 100644
--- a/examples/3d/pod/fracture3_pod.py
+++ b/examples/3d/pod/fracture3_pod.py
@@ -1,251 +1,249 @@
import numpy as np
from mpl_toolkits.mplot3d import Axes3D
from matplotlib import pyplot as plt
from rrompy.hfengines.linear_problem.tridimensional import \
MembraneFractureEngine3 as MFE
from rrompy.reduction_methods.standard import RationalInterpolant as RI
from rrompy.reduction_methods.standard import ReducedBasis as RB
-from rrompy.reduction_methods.pole_matching import \
- RationalInterpolantPoleMatching as RIPM
-from rrompy.reduction_methods.pole_matching import \
- ReducedBasisPoleMatching as RBPM
+from rrompy.reduction_methods.pivoted import RationalInterpolantPivoted as RIP
+from rrompy.reduction_methods.pivoted import ReducedBasisPivoted as RBP
from rrompy.parameter.parameter_sampling import (QuadratureSampler as QS,
QuadratureSamplerTotal as QST,
ManualSampler as MS,
RandomSampler as RS)
verb = 70
size = 4
show_sample = False
show_norm = False
clip = -1
#clip = .4
#clip = .6
Delta = 0
MN = 5
R = (MN + 3) * (MN + 2) * (MN + 1) // 6
S = [int(np.ceil(R ** (1. / 3.)))] * 3
PODTol = 1e-8
SPivot = [MN + 1, 5, 5]
MMarginal = SPivot[1] - 1
matchingWeight = 10.
samples = "centered"
samples = "standard"
-#samples = "pole matching"
+#samples = "pivoted"
algo = "rational"
#algo = "RB"
sampling = "quadrature"
#sampling = "quadrature_total"
#sampling = "random"
samplingM = "quadrature"
samplingM = "quadrature_total"
#samplingM = "random"
polydegreetype = "TOTAL"
#polydegreetype = "FULL"
if samples == "standard":
radial = ""
# radial = "_gaussian"
# radial = "_thinplate"
# radial = "_multiquadric"
rW0 = 5.
radialWeight = [rW0] * 3
-if samples == "pole matching":
+if samples == "pivoted":
radial = ""
# radial = "_gaussian"
# radial = "_thinplate"
# radial = "_multiquadric"
rW0 = 5.
radialWeight = [rW0]
radialM = ""
# radialM = "_gaussian"
# radialM = "_thinplate"
# radialM = "_multiquadric"
rMW0 = 2.
radialWeightM = [rMW0] * 2
assert Delta <= 0
if size == 1:
mu0 = [47.5 ** .5, .4, .05]
mutar = [50 ** .5, .45, .07]
murange = [[40 ** .5, .3, .025], [55 ** .5, .5, .075]]
if size == 2:
mu0 = [50 ** .5, .3, .02]
mutar = [55 ** .5, .4, .03]
murange = [[40 ** .5, .1, .005], [60 ** .5, .5, .035]]
if size == 3:
mu0 = [45 ** .5, .5, .05]
mutar = [47 ** .5, .4, .045]
murange = [[40 ** .5, .3, .04], [50 ** .5, .7, .06]]
if size == 4:
mu0 = [45 ** .5, .4, .05]
mutar = [47 ** .5, .45, .045]
murange = [[40 ** .5, .3, .04], [50 ** .5, .5, .06]]
if size == 5:
mu0 = [45 ** .5, .5, .05]
mutar = [47 ** .5, .45, .045]
murange = [[40 ** .5, .3, .04], [50 ** .5, .7, .06]]
aEff = 1.#25
bEff = 1. - aEff
murangeEff = [[(aEff*murange[0][0]**2. + bEff*murange[1][0]**2.) ** .5,
aEff*murange[0][1] + bEff*murange[1][1],
aEff*murange[0][2] + bEff*murange[1][2]],
[(aEff*murange[1][0]**2. + bEff*murange[0][0]**2.) ** .5,
aEff*murange[1][1] + bEff*murange[0][1],
aEff*murange[1][2] + bEff*murange[0][2]]]
H = 1.
L = .75
delta = .05
n = 50
solver = MFE(mu0 = mu0, H = H, L = L, delta = delta, n = n, verbosity = verb)
rescalingExp = [2., 1., 1.]
if algo == "rational":
params = {'N':MN, 'M':MN + Delta, 'S':S, 'POD':True,
'polydegreetype':polydegreetype}
if samples == "standard":
params['polybasis'] = "CHEBYSHEV" + radial
# params['polybasis'] = "LEGENDRE" + radial
# params['polybasis'] = "MONOMIAL" + radial
params['radialDirectionalWeights'] = radialWeight
method = RI
elif samples == "centered":
params['polybasis'] = "MONOMIAL"
params['S'] = R
method = RI
- elif samples == "pole matching":
+ elif samples == "pivoted":
params['S'] = [SPivot[0]]
params['SMarginal'] = SPivot[1:]
params['MMarginal'] = MMarginal
params['polybasisPivot'] = "CHEBYSHEV" + radial
params['polybasisMarginal'] = "MONOMIAL" + radialM
params['matchingWeight'] = matchingWeight
params['radialDirectionalWeightsPivot'] = radialWeight
params['radialDirectionalWeightsMarginal'] = radialWeightM
- method = RIPM
+ method = RIP
else: #if algo == "RB":
params = {'R':(MN + 3 + Delta) * (MN + 2 + Delta) * (MN + 1 + Delta) // 6,
'S':S, 'POD':True, 'PODTolerance':PODTol}
if samples == "standard":
method = RB
elif samples == "centered":
params['S'] = R
method = RB
- elif samples == "pole matching":
+ elif samples == "pivoted":
params['S'] = [SPivot[0]]
params['SMarginal'] = SPivot[1:]
params['MMarginal'] = MMarginal
params['polybasisMarginal'] = "MONOMIAL" + radialM
params['matchingWeight'] = matchingWeight
params['radialDirectionalWeightsMarginal'] = radialWeightM
- method = RBPM
+ method = RBP
if samples == "standard":
if sampling == "quadrature":
params['sampler'] = QS(murange, "CHEBYSHEV", scalingExp = rescalingExp)
# params['sampler'] = QS(murange, "GAUSSLEGENDRE",scalingExp = rescalingExp)
# params['sampler'] = QS(murange, "UNIFORM", scalingExp = rescalingExp)
params['S'] = [max(j, MN + 1) for j in params['S']]
elif sampling == "quadrature_total":
params['sampler'] = QST(murange, "CHEBYSHEV", scalingExp = rescalingExp)
params['S'] = R
else: # if sampling == "random":
params['sampler'] = RS(murange, "HALTON", scalingExp = rescalingExp)
params['S'] = R
elif samples == "centered":
params['sampler'] = MS(murange, points = [mu0], scalingExp = rescalingExp)
-elif samples == "pole matching":
+elif samples == "pivoted":
if sampling == "quadrature":
params['samplerPivot'] = QS([murange[0][0], murange[1][0]], "CHEBYSHEV")
# params['samplerPivot'] = QS([murange[0][0], murange[1][0]], "GAUSSLEGENDRE")
# params['samplerPivot'] = QS([murange[0][0], murange[1][0]], "UNIFORM")
elif sampling == "quadrature_total":
params['samplerPivot'] = QST([murange[0][0], murange[1][0]], "CHEBYSHEV")
params['S'] = MN + 1
else: # if sampling == "random":
params['samplerPivot'] = RS([murange[0][0], murange[1][0]], "HALTON")
params['S'] = MN + 1
if samplingM == "quadrature":
params['samplerMarginal'] = QS([murange[0][1:], murange[1][1:]], "UNIFORM")
elif samplingM == "quadrature_total":
params['samplerMarginal'] = QST([murange[0][1:], murange[1][1:]], "CHEBYSHEV")
params['SMarginal'] = (MMarginal + 2) * (MMarginal + 1) // 2
params['polybasisMarginal'] = "CHEBYSHEV" + radialM
else: # if samplingM == "random":
params['samplerMarginal'] = RS([murange[0][1:], murange[1][1:]], "HALTON")
params['SMarginal'] = (MMarginal + 2) * (MMarginal + 1) // 2
-if samples != "pole matching":
+if samples != "pivoted":
approx = method(solver, mu0 = mu0, approxParameters = params,
verbosity = verb)
else:
approx = method(solver, mu0 = mu0, directionPivot = [0],
approxParameters = params, verbosity = verb)
if samples != "centered": approx.samplingEngine.allowRepeatedSamples = False
approx.setupApprox()
if show_sample:
from fracture3_warping import fracture3_warping
warps = fracture3_warping(solver.V.mesh(), L, mutar[1], delta, mutar[2])
approx.plotApprox(mutar, warps, name = 'u_app', what = "REAL")
approx.plotHF(mutar, warps, name = 'u_HF', what = "REAL")
approx.plotErr(mutar, warps, name = 'err', what = "REAL")
# approx.plotRes(mutar, warps, name = 'res', what = "REAL")
appErr = approx.normErr(mutar)
solNorm = approx.normHF(mutar)
resNorm = approx.normRes(mutar)
RHSNorm = approx.normRHS(mutar)
print(('SolNorm:\t{}\nErr:\t{}\nErrRel:\t{}').format(solNorm, appErr,
np.divide(appErr, solNorm)))
print(('RHSNorm:\t{}\nRes:\t{}\nResRel:\t{}').format(RHSNorm, resNorm,
np.divide(resNorm, RHSNorm)))
fig = plt.figure(figsize = (8, 6))
ax = Axes3D(fig)
ax.scatter(np.real(approx.mus(0) ** 2.), np.real(approx.mus(1)),
np.real(approx.mus(2)), '.')
plt.show()
plt.close()
approx.verbosity = 0
approx.trainedModel.verbosity = 0
from plot_zero_set_3 import plotZeroSet3
#muZeroVals, Qvals = plotZeroSet3(murange, murangeEff, approx, mu0, 200,
# [None, mu0[1], mu0[2]], [2., 1., 1.],
# clip = clip)
#muZeroVals, Qvals = plotZeroSet3(murange, murangeEff, approx, mu0, 200,
# [None, None, mu0[2]], [2., 1., 1.],
# clip = clip)
#muZeroVals, Qvals = plotZeroSet3(murange, murangeEff, approx, mu0, 200,
# [None, mu0[1], None], [2., 1., 1.],
# clip = clip)
muZeroScatter = plotZeroSet3(murange, murangeEff, approx, mu0, 50,
[None, None, None], [2., 1., 1.], clip = clip)
if show_norm:
solver._solveBatchSize = 25
from plot_inf_set_3 import plotInfSet3
muInfVals, normEx, normApp, normRes, normErr, beta = plotInfSet3(
murange, murangeEff, approx, mu0, 25,
[None, mu0[1], mu0[2]], [2., 1., 1.],
clip = clip, relative = False,
normalizeDen = True)
muInfVals, normEx, normApp, normRes, normErr, beta = plotInfSet3(
murange, murangeEff, approx, mu0, 25,
[None, None, mu0[2]], [2., 1., 1.],
clip = clip, relative = False,
normalizeDen = True)
muInfVals, normEx, normApp, normRes, normErr, beta = plotInfSet3(
murange, murangeEff, approx, mu0, 25,
[None, mu0[1], None], [2., 1., 1.],
clip = clip, relative = False,
normalizeDen = True)
print(approx.getPoles([None, .5, .05]) ** 2.)
diff --git a/rrompy/reduction_methods/pivoted/generic_pivoted_approximant.py b/rrompy/reduction_methods/pivoted/generic_pivoted_approximant.py
index ea31f75..b77b853 100644
--- a/rrompy/reduction_methods/pivoted/generic_pivoted_approximant.py
+++ b/rrompy/reduction_methods/pivoted/generic_pivoted_approximant.py
@@ -1,451 +1,532 @@
# 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 .
#
import numpy as np
from copy import deepcopy as copy
from rrompy.reduction_methods.base.generic_approximant import (
GenericApproximant)
from rrompy.utilities.poly_fitting.polynomial import (polybases as ppb,
nextDerivativeIndices,
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.sampling.pivoted import (SamplingEnginePivoted,
SamplingEnginePivotedPOD)
from rrompy.utilities.base.types import (ListAny, DictAny, HFEng, paramVal,
paramList)
from rrompy.utilities.base import verbosityManager as vbMng
from rrompy.utilities.exception_manager import (RROMPyException, RROMPyAssert,
RROMPyWarning)
from rrompy.parameter import emptyParameterList
__all__ = ['GenericPivotedApproximant']
class GenericPivotedApproximant(GenericApproximant):
"""
- ROM pivoted approximant computation for parametric problems (ABSTRACT).
+ ROM pivoted approximant (with pole matching) computation for parametric
+ problems (ABSTRACT).
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;
- 'polybasisMarginal': type of polynomial basis for marginal
- interpolation; defaults to 'MONOMIAL';
+ interpolation; allowed values include 'MONOMIAL', 'CHEBYSHEV'
+ and 'LEGENDRE'; defaults to 'MONOMIAL';
- 'MMarginal': degree of marginal interpolant; defaults to 0;
+ - 'polydegreetypeMarginal': type of polynomial degree for marginal;
+ defaults to 'TOTAL';
- 'radialDirectionalWeightsMarginal': radial basis weights for
marginal interpolant; defaults to 0, i.e. identity;
- 'interpRcondMarginal': tolerance for marginal interpolation;
defaults to None.
Defaults to empty dict.
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;
- 'polybasisMarginal': type of polynomial basis for marginal
interpolation;
- 'MMarginal': degree of marginal interpolant;
+ - 'polydegreetypeMarginal': type of polynomial degree for marginal;
- 'radialDirectionalWeightsMarginal': radial basis weights for
marginal interpolant;
- 'interpRcondMarginal': tolerance for marginal interpolation.
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.
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.
samplerPivot: Pivot sample point generator.
SMarginal: Total number of marginal samples current approximant relies
upon.
samplerMarginal: Marginal sample point generator.
polybasisMarginal: Type of polynomial basis for marginal interpolation.
MMarginal: Degree of marginal interpolant.
+ polydegreetypeMarginal: Type of polynomial degree for marginal.
radialDirectionalWeightsMarginal: Radial basis weights for marginal
interpolant.
interpRcondMarginal: Tolerance for marginal interpolation.
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.
"""
def __init__(self, HFEngine:HFEng, mu0 : paramVal = None,
directionPivot : ListAny = [0],
approxParameters : DictAny = {}, verbosity : int = 10,
timestamp : bool = True):
self._preInit()
+ if len(directionPivot) > 1:
+ raise RROMPyException(("Exactly 1 pivot parameter allowed in pole "
+ "matching."))
from rrompy.parameter.parameter_sampling import QuadratureSampler as QS
QSBase = QS([[0], [1]], "UNIFORM")
- self._addParametersToList(["polybasisMarginal", "MMarginal",
+ self._addParametersToList(["matchingWeight", "cutOffTolerance",
+ "cutOffType", "polybasisMarginal",
+ "MMarginal", "polydegreetypeMarginal",
"radialDirectionalWeightsMarginal",
"interpRcondMarginal"],
- ["MONOMIAL", 0, 1, -1],
+ [1, np.inf, "MAGNITUDE", "MONOMIAL", 0,
+ "TOTAL", 1, -1],
["samplerPivot", "SMarginal",
"samplerMarginal"], [QSBase, [1], QSBase])
del QS
self._directionPivot = directionPivot
super().__init__(HFEngine = HFEngine, mu0 = mu0,
approxParameters = approxParameters,
verbosity = verbosity, timestamp = timestamp)
self._postInit()
def setupSampling(self):
"""Setup sampling engine."""
RROMPyAssert(self._mode, message = "Cannot setup sampling engine.")
if not hasattr(self, "_POD") or self._POD is None: return
if self.POD:
SamplingEngine = SamplingEnginePivotedPOD
else:
SamplingEngine = SamplingEnginePivoted
self.samplingEngine = SamplingEngine(self.HFEngine,
self.directionPivot,
verbosity = self.verbosity,
allowRepeatedSamples = True)
+ @property
+ def matchingWeight(self):
+ """Value of matchingWeight."""
+ return self._matchingWeight
+ @matchingWeight.setter
+ def matchingWeight(self, matchingWeight):
+ self._matchingWeight = matchingWeight
+ self._approxParameters["matchingWeight"] = self.matchingWeight
+
+ @property
+ def cutOffTolerance(self):
+ """Value of cutOffTolerance."""
+ return self._cutOffTolerance
+ @cutOffTolerance.setter
+ def cutOffTolerance(self, cutOffTolerance):
+ self._cutOffTolerance = cutOffTolerance
+ self._approxParameters["cutOffTolerance"] = self.cutOffTolerance
+
+ @property
+ def cutOffType(self):
+ """Value of cutOffType."""
+ return self._cutOffType
+ @cutOffType.setter
+ def cutOffType(self, cutOffType):
+ try:
+ cutOffType = cutOffType.upper().strip().replace(" ","")
+ if cutOffType not in ["MAGNITUDE", "POTENTIAL"]:
+ raise RROMPyException("Prescribed cutOffType not recognized.")
+ self._cutOffType = cutOffType
+ except:
+ RROMPyWarning(("Prescribed cutOffType not recognized. Overriding "
+ "to 'MAGNITUDE'."))
+ self._cutOffType = "MAGNITUDE"
+ self._approxParameters["cutOffType"] = self.cutOffType
+
@property
def SMarginal(self):
"""Value of SMarginal."""
return self._SMarginal
@SMarginal.setter
def SMarginal(self, SMarginal):
if not hasattr(SMarginal, "__len__"): SMarginal = [SMarginal]
if any([s <= 0 for s in SMarginal]):
raise RROMPyException("SMarginal must be positive.")
if hasattr(self, "_SMarginal") and self._SMarginal is not None:
Sold = tuple(self.SMarginal)
else: Sold = -1
self._SMarginal = SMarginal
self._approxParameters["SMarginal"] = self.SMarginal
if Sold != tuple(self.SMarginal):
self.resetSamples()
@property
def polybasisMarginal(self):
"""Value of polybasisMarginal."""
return self._polybasisMarginal
@polybasisMarginal.setter
def polybasisMarginal(self, polybasisMarginal):
try:
polybasisMarginal = polybasisMarginal.upper().strip().replace(" ",
"")
if polybasisMarginal not in ppb + rbpb + mlspb:
raise RROMPyException(
"Prescribed marginal polybasis not recognized.")
self._polybasisMarginal = polybasisMarginal
except:
RROMPyWarning(("Prescribed marginal polybasis not recognized. "
"Overriding to 'MONOMIAL'."))
self._polybasisMarginal = "MONOMIAL"
self._approxParameters["polybasisMarginal"] = self.polybasisMarginal
@property
def MMarginal(self):
"""Value of MMarginal."""
return self._MMarginal
@MMarginal.setter
def MMarginal(self, MMarginal):
if MMarginal < 0:
raise RROMPyException("MMarginal must be non-negative.")
self._MMarginal = MMarginal
self._approxParameters["MMarginal"] = self.MMarginal
+ @property
+ def polydegreetypeMarginal(self):
+ """Value of polydegreetypeMarginal."""
+ return self._polydegreetypeMarginal
+ @polydegreetypeMarginal.setter
+ def polydegreetypeMarginal(self, polydegreetypeM):
+ try:
+ polydegreetypeM = polydegreetypeM.upper().strip().replace(" ","")
+ if polydegreetypeM not in ["TOTAL", "FULL"]:
+ raise RROMPyException(("Prescribed polydegreetypeMarginal not "
+ "recognized."))
+ self._polydegreetypeMarginal = polydegreetypeM
+ except:
+ RROMPyWarning(("Prescribed polydegreetypeMarginal not recognized. "
+ "Overriding to 'TOTAL'."))
+ self._polydegreetypeMarginal = "TOTAL"
+ self._approxParameters["polydegreetypeMarginal"] = (
+ self.polydegreetypeMarginal)
+
@property
def radialDirectionalWeightsMarginal(self):
"""Value of radialDirectionalWeightsMarginal."""
return self._radialDirectionalWeightsMarginal
@radialDirectionalWeightsMarginal.setter
def radialDirectionalWeightsMarginal(self, radialDirWeightsMarginal):
self._radialDirectionalWeightsMarginal = radialDirWeightsMarginal
self._approxParameters["radialDirectionalWeightsMarginal"] = (
self.radialDirectionalWeightsMarginal)
@property
def interpRcondMarginal(self):
"""Value of interpRcondMarginal."""
return self._interpRcondMarginal
@interpRcondMarginal.setter
def interpRcondMarginal(self, interpRcondMarginal):
self._interpRcondMarginal = interpRcondMarginal
self._approxParameters["interpRcondMarginal"] = (
self.interpRcondMarginal)
@property
def directionPivot(self):
"""Value of directionPivot. Its assignment may reset snapshots."""
return self._directionPivot
@directionPivot.setter
def directionPivot(self, directionPivot):
if hasattr(self, '_directionPivot'):
directionPivotOld = copy(self.directionPivot)
else:
directionPivotOld = None
if (directionPivotOld is None
or len(directionPivot) != len(directionPivotOld)
or not directionPivot == directionPivotOld):
self.resetSamples()
self._directionPivot = directionPivot
@property
def directionMarginal(self):
return [x for x in range(self.HFEngine.npar) \
if x not in self.directionPivot]
@property
def rescalingExpPivot(self):
return [self.HFEngine.rescalingExp[x] for x in self.directionPivot]
@property
def rescalingExpMarginal(self):
return [self.HFEngine.rescalingExp[x] for x in self.directionMarginal]
@property
def muBoundsPivot(self):
"""Value of muBoundsPivot."""
return self.samplerPivot.lims
@property
def muBoundsMarginal(self):
"""Value of muBoundsMarginal."""
return self.samplerMarginal.lims
@property
def samplerPivot(self):
"""Value of samplerPivot."""
return self._samplerPivot
@samplerPivot.setter
def samplerPivot(self, samplerPivot):
if 'generatePoints' not in dir(samplerPivot):
raise RROMPyException("Pivot sampler type not recognized.")
if hasattr(self, '_samplerPivot') and self._samplerPivot is not None:
samplerOld = self.samplerPivot
self._samplerPivot = samplerPivot
self._approxParameters["samplerPivot"] = self.samplerPivot.__str__()
if not 'samplerOld' in locals() or samplerOld != self.samplerPivot:
self.resetSamples()
@property
def samplerMarginal(self):
"""Value of samplerMarginal."""
return self._samplerMarginal
@samplerMarginal.setter
def samplerMarginal(self, samplerMarginal):
if 'generatePoints' not in dir(samplerMarginal):
raise RROMPyException("Marginal sampler type not recognized.")
if (hasattr(self, '_samplerMarginal')
and self._samplerMarginal is not None):
samplerOld = self.samplerMarginal
self._samplerMarginal = samplerMarginal
self._approxParameters["samplerMarginal"] = (
self.samplerMarginal.__str__())
if not 'samplerOld' in locals() or samplerOld != self.samplerMarginal:
self.resetSamples()
def resetSamples(self):
"""Reset samples."""
super().resetSamples()
self._musMUniqueCN = None
self._derMIdxs = None
self._reorderM = None
def setSamples(self, samplingEngine):
"""Copy samplingEngine and samples."""
super().setSamples(samplingEngine)
self.mus = copy(self.samplingEngine[0].mus)
for sEj in self.samplingEngine[1:]:
self.mus.append(sEj.mus)
def computeSnapshots(self):
"""Compute snapshots of solution map."""
RROMPyAssert(self._mode,
message = "Cannot start snapshot computation.")
Sp = np.prod(self.S)
Seff = Sp * np.prod(self.SMarginal)
if self.samplingEngine.nsamplesTot != Seff:
self.resetSamples()
vbMng(self, "INIT", "Starting computation of snapshots.", 5)
if self.HFEngine.As[0] is None: self.HFEngine.A(self.mu0)
if self.HFEngine.bs[0] is None: self.HFEngine.b(self.mu0)
self.musPivot = self.samplerPivot.generatePoints(self.S)
self.musMarginal = self.samplerMarginal.generatePoints(
self.SMarginal)
self.mus = emptyParameterList()
self.mus.reset((Seff, self.HFEngine.npar))
self.samplingEngine.resetHistory(Seff // Sp)
for j, muMarg in enumerate(self.musMarginal):
for k in range(j * Sp, (j + 1) * Sp):
self.mus.data[k, self.directionPivot] = (
self.musPivot[k - j * Sp].data)
self.mus.data[k, self.directionMarginal] = muMarg.data
self.samplingEngine.iterSample(self.musPivot, self.musMarginal)
if self.POD:
self.samplingEngine.coalesceSamples(self.interpRcondMarginal)
else:
self.samplingEngine.coalesceSamples()
vbMng(self, "DEL", "Done computing snapshots.", 5)
def _setupMarginalInterpolationIndices(self):
"""Setup parameters for polyvander."""
RROMPyAssert(self._mode,
message = "Cannot setup interpolation indices.")
if (self._musMUniqueCN is None
or len(self._reorderM) != len(self.musMarginal)):
self._musMUniqueCN, musMIdxsTo, musMIdxs, musMCount = (
self.centerNormalizeMarginal(self.musMarginal).unique(
return_index = True,
return_inverse = True,
return_counts = True))
self._musMUnique = self.musMarginal[musMIdxsTo]
self._derMIdxs = [None] * len(self._musMUniqueCN)
self._reorderM = np.empty(len(musMIdxs), dtype = int)
filled = 0
for j, cnt in enumerate(musMCount):
self._derMIdxs[j] = nextDerivativeIndices([],
self.musMarginal.shape[1], cnt)
jIdx = np.nonzero(musMIdxs == j)[0]
self._reorderM[jIdx] = np.arange(filled, filled + cnt)
filled += cnt
def _setupMarginalInterp(self):
"""Compute marginal interpolator."""
RROMPyAssert(self._mode, message = "Cannot setup numerator.")
vbMng(self, "INIT", "Starting computation of marginal interpolator.",
7)
self._setupMarginalInterpolationIndices()
MMarginal0 = copy(self.MMarginal)
mI = []
for j in range(len(self.musMarginal)):
canonicalj = 1. * (np.arange(len(self.musMarginal)) == j)
self._MMarginal = MMarginal0
while self.MMarginal >= 0:
if self.polybasisMarginal in ppb:
p = PI()
wellCond, msg = p.setupByInterpolation(
- self._musMUniqueCN, canonicalj, self.MMarginal,
- self.polybasisMarginal, self.verbosity >= 5, True,
- {"derIdxs": self._derMIdxs,
- "reorder": self._reorderM,
- "scl": np.power(self.scaleFactorMarginal, -1.)},
- {"rcond": self.interpRcondMarginal})
+ self._musMUniqueCN, canonicalj, self.MMarginal,
+ self.polybasisMarginal, self.verbosity >= 5,
+ self.polydegreetypeMarginal == "TOTAL",
+ {"derIdxs": self._derMIdxs,
+ "reorder": self._reorderM,
+ "scl": np.power(self.scaleFactorMarginal, -1.)},
+ {"rcond": self.interpRcondMarginal})
elif self.polybasisMarginal in rbpb:
p = RBI()
wellCond, msg = p.setupByInterpolation(
self._musMUniqueCN, canonicalj, self.MMarginal,
self.polybasisMarginal,
self.radialDirectionalWeightsMarginal,
- self.verbosity >= 5, True,
+ self.verbosity >= 5,
+ self.polydegreetypeMarginal == "TOTAL",
{"derIdxs": self._derMIdxs,
"reorder": self._reorderM,
"scl": np.power(self.scaleFactorMarginal, -1.)},
{"rcond": self.interpRcondMarginal})
else:# if self.polybasisMarginal in mlspb:
p = MLSI()
wellCond, msg = p.setupByInterpolation(
self._musMUniqueCN, canonicalj, self.MMarginal,
self.polybasisMarginal,
self.radialDirectionalWeightsMarginal,
- self.verbosity >= 5, True,
+ self.verbosity >= 5,
+ self.polydegreetypeMarginal == "TOTAL",
{"derIdxs": self._derMIdxs,
"reorder": self._reorderM,
"scl": np.power(self.scaleFactorMarginal, -1.)})
vbMng(self, "MAIN", msg, 5)
if wellCond: break
RROMPyWarning(("Polyfit is poorly conditioned. Reducing "
"MMarginal by 1."))
self.MMarginal -= 1
if self.MMarginal < 0:
raise RROMPyException(("Instability in computation of "
"marginal interpolator. Aborting."))
mI = mI + [copy(p)]
vbMng(self, "DEL", "Done computing marginal interpolator.", 7)
return mI
def normApprox(self, mu:paramList) -> float:
"""
Compute norm of approximant at arbitrary parameter.
Args:
mu: Target parameter.
Returns:
Target norm of approximant.
"""
if not self.POD: return super().normApprox(mu)
return np.linalg.norm(self.getApproxReduced(mu).data, axis = 0)
def computeScaleFactor(self):
"""Compute parameter rescaling factor."""
RROMPyAssert(self._mode, message = "Cannot compute rescaling factor.")
self.scaleFactorPivot = .5 * np.abs(
self.muBoundsPivot[0] ** self.rescalingExpPivot
- self.muBoundsPivot[1] ** self.rescalingExpPivot)
self.scaleFactorMarginal = .5 * np.abs(
self.muBoundsMarginal[0] ** self.rescalingExpMarginal
- self.muBoundsMarginal[1] ** self.rescalingExpMarginal)
self.scaleFactor = np.empty(self.npar)
self.scaleFactor[self.directionPivot] = self.scaleFactorPivot
self.scaleFactor[self.directionMarginal] = self.scaleFactorMarginal
def centerNormalizeMarginal(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.centerNormalizeMarginal(mu, mu0)
diff --git a/rrompy/reduction_methods/pivoted/rational_interpolant_pivoted.py b/rrompy/reduction_methods/pivoted/rational_interpolant_pivoted.py
index aae415a..dfb4219 100644
--- a/rrompy/reduction_methods/pivoted/rational_interpolant_pivoted.py
+++ b/rrompy/reduction_methods/pivoted/rational_interpolant_pivoted.py
@@ -1,602 +1,630 @@
# 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 .
#
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,
nextDerivativeIndices,
hashDerivativeToIdx as hashD,
hashIdxToDerivative as hashI,
polyvander as pvP,
polyvanderTotal as pvTP,
degreeTotalToFull,
fullDegreeMaxMask,
totalDegreeMaxMask,
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.reduction_methods.trained_model import (TrainedModelPivotedData,
TrainedModelPivotedRational as tModel)
from rrompy.utilities.base.types import (Np1D, Np2D, HFEng, DictAny, Tuple,
List, ListAny, paramVal, paramList)
from rrompy.utilities.base import verbosityManager as vbMng, freepar as fp
from rrompy.utilities.numerical import multifactorial, customPInv
from rrompy.utilities.exception_manager import (RROMPyException, RROMPyAssert,
RROMPyWarning)
from rrompy.parameter import checkParameter
__all__ = ['RationalInterpolantPivoted']
class RationalInterpolantPivoted(GenericPivotedApproximant):
"""
- ROM pivoted rational interpolant computation for parametric problems.
+ 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; defaults to 'MONOMIAL';
+ 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;
- '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.
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;
- '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.
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.
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 = {}, verbosity : int = 10,
timestamp : bool = True):
self._preInit()
self._addParametersToList(["polybasisPivot", "M", "N",
"polydegreetype",
"radialDirectionalWeightsPivot",
"interpRcondPivot", "robustTol"],
["MONOMIAL", 0, 0, "TOTAL", 1, -1, 0])
super().__init__(HFEngine = HFEngine, mu0 = mu0,
directionPivot = directionPivot,
approxParameters = approxParameters,
verbosity = verbosity, timestamp = timestamp)
self._postInit()
@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 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
@property
def S(self):
"""Value of S."""
return self._S
@S.setter
def S(self, S):
GenericPivotedApproximant.S.fset(self, S)
if hasattr(self, "_M"): self.M = self.M
if hasattr(self, "_N"): self.N = self.N
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.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.musPivot.shape[1], 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.musPivot.shape[1]
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()
M0 = copy(self.M)
tensor_idx = 0
ps = []
for k, muM in enumerate(self.musMarginal):
self._M = M0
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.musPivot.shape[1])
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.)},
{"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.)})
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 -= 1
if self.M < 0:
raise RROMPyException(("Instability in computation of "
"numerator. Aborting."))
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.computeScaleFactor()
self.computeSnapshots()
if self.trainedModel is None:
self.trainedModel = tModel()
self.trainedModel.verbosity = self.verbosity
self.trainedModel.timestamp = self.timestamp
data = TrainedModelPivotedData(self.trainedModel.name(), self.mu0,
self.samplingEngine.samplesCoalesced,
self.scaleFactor,
self.HFEngine.rescalingExp,
self.directionPivot)
data.musPivot = copy(self.musPivot)
data.musMarginal = copy(self.musMarginal)
self.trainedModel.data = data
else:
self.trainedModel = self.trainedModel
self.trainedModel.data.projMat = copy(
self.samplingEngine.samplesCoalesced)
self.trainedModel.data.marginalInterp = self._setupMarginalInterp()
if self.N > 0:
Qs = self._setupDenominator()[0]
else:
Q = PI()
Q.npar = self.musPivot.shape[1]
Q.coeffs = np.ones(tuple([1] * Q.npar), dtype = np.complex)
Q.polybasis = self.polybasisPivot0
Qs = [Q for _ in range(len(self.musMarginal))]
+ self.trainedModel.data._temporary = True
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)
+ 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()
nparPivot = self.musPivot.shape[1]
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, nparPivot)
else: #if self.polydegreetype == "FULL":
TE = pvP(self._musPUniqueCN, [self.N] * nparPivot,
self.polybasisPivot0, self._derPIdxs, self._reorderP,
scl = np.power(self.scaleFactorPivot, -1.))
idxsB = fullDegreeMaxMask(self.N, 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] = pseudoInv.dot(TN)
return invD, fitinvP
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 centerNormalizePivot(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.centerNormalizePivot(mu, mu0)
def getResidues(self, *args, **kwargs) -> Np1D:
"""
Obtain approximant residues.
Returns:
Matrix with residues as columns.
"""
return self.trainedModel.getResidues(*args, **kwargs)
diff --git a/rrompy/reduction_methods/pivoted/reduced_basis_pivoted.py b/rrompy/reduction_methods/pivoted/reduced_basis_pivoted.py
index 33427ba..4ce4b49 100644
--- a/rrompy/reduction_methods/pivoted/reduced_basis_pivoted.py
+++ b/rrompy/reduction_methods/pivoted/reduced_basis_pivoted.py
@@ -1,273 +1,302 @@
# 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 .
#
from copy import deepcopy as copy
import numpy as np
from .generic_pivoted_approximant import GenericPivotedApproximant
from rrompy.hfengines.base import MarginalProxyEngine
from rrompy.reduction_methods.trained_model import (TrainedModelPivotedData,
TrainedModelPivotedReducedBasis as tModel)
from rrompy.reduction_methods.base.reduced_basis_utils import \
projectAffineDecomposition
from rrompy.utilities.base.types import (Np1D, Np2D, List, ListAny, Tuple,
DictAny, HFEng, paramVal, sampList)
from rrompy.utilities.base import verbosityManager as vbMng, freepar as fp
from rrompy.utilities.numerical import customPInv
from rrompy.utilities.exception_manager import (RROMPyWarning, RROMPyException,
RROMPyAssert)
__all__ = ['ReducedBasisPivoted']
class ReducedBasisPivoted(GenericPivotedApproximant):
"""
- ROM pivoted RB approximant computation for parametric problems.
+ ROM pivoted RB approximant (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;
- 'R': rank for pivot Galerkin projection; defaults to prod(S);
- 'PODTolerance': tolerance for pivot snapshots POD; defaults to
-1;
- 'polybasisMarginal': type of polynomial basis for marginal
interpolation; allowed values include 'MONOMIAL', 'CHEBYSHEV'
and 'LEGENDRE'; defaults to 'MONOMIAL';
- 'MMarginal': degree of marginal interpolant; defaults to 0;
+ - 'polydegreetypeMarginal': type of polynomial degree for marginal;
+ defaults to 'TOTAL';
- 'radialDirectionalWeightsMarginal': radial basis weights for
marginal interpolant; defaults to 0, i.e. identity;
- 'interpRcondMarginal': tolerance for marginal interpolation;
defaults to None.
Defaults to empty dict.
verbosity(optional): Verbosity level. Defaults to 10.
Attributes:
HFEngine: HF problem solver.
mu0: Default parameter.
directionPivot: Pivot components.
mus: Array of snapshot parameters.
- approxRadius: Dummy radius of approximant (i.e. distance from mu0 to
- farthest sample point).
+ 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;
- 'R': rank for Galerkin projection;
- 'PODTolerance': tolerance for snapshots POD;
- 'polybasisMarginal': type of polynomial basis for marginal
interpolation;
- 'MMarginal': degree of marginal interpolant;
+ - 'polydegreetypeMarginal': type of polynomial degree for marginal;
- 'radialDirectionalWeightsMarginal': radial basis weights for
marginal interpolant;
- 'interpRcondMarginal': tolerance for marginal interpolation.
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.
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.
samplerPivot: Pivot sample point generator.
SMarginal: Total number of marginal samples current approximant relies
upon.
samplerMarginal: Marginal sample point generator.
R: Rank for Galerkin projection.
PODTolerance: Tolerance for pivot snapshots POD.
polybasisMarginal: Type of polynomial basis for marginal interpolation.
MMarginal: Degree of marginal interpolant.
+ polydegreetypeMarginal: Type of polynomial degree for marginal.
radialDirectionalWeightsMarginal: Radial basis weights for marginal
interpolant.
interpRcondMarginal: Tolerance for marginal interpolation.
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.
+ 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.
As: List of sparse matrices (in CSC format) representing coefficients
of linear system matrix.
bs: List of numpy vectors representing coefficients of linear system
RHS.
ARBs: List of sparse matrices (in CSC format) representing coefficients
of compressed linear system matrix.
bRBs: List of numpy vectors representing coefficients of compressed
linear system RHS.
"""
def __init__(self, HFEngine:HFEng, mu0 : paramVal = None,
directionPivot : ListAny = [0],
approxParameters : DictAny = {}, verbosity : int = 10,
timestamp : bool = True):
self._preInit()
self._addParametersToList(["R", "PODTolerance"], [1, -1])
super().__init__(HFEngine = HFEngine, mu0 = mu0,
directionPivot = directionPivot,
approxParameters = approxParameters,
verbosity = verbosity, timestamp = timestamp)
self._postInit()
@property
def S(self):
"""Value of S."""
return self._S
@S.setter
def S(self, S):
GenericPivotedApproximant.S.fset(self, S)
if not hasattr(self, "_R"):
self._R = np.prod(self.S) * np.prod(self.SMarginal)
@property
def R(self):
"""Value of R. Its assignment may change S."""
return self._R
@R.setter
def R(self, R):
if R < 0: raise RROMPyException("R must be non-negative.")
self._R = R
self._approxParameters["R"] = self.R
if (hasattr(self, "_S") and hasattr(self, "_SMarginal")
and np.prod(self.S) * np.prod(self.SMarginal) < self.R):
RROMPyWarning(("Prescribed S and SMarginal are too small. "
"Decreasing R."))
self.R = np.prod(self.S) * np.prod(self.SMarginal)
@property
def PODTolerance(self):
"""Value of PODTolerance."""
return self._PODTolerance
@PODTolerance.setter
def PODTolerance(self, PODTolerance):
self._PODTolerance = PODTolerance
self._approxParameters["PODTolerance"] = self.PODTolerance
def _setupProjectionMatrix(self):
"""Compute projection matrix."""
RROMPyAssert(self._mode, message = "Cannot setup numerator.")
vbMng(self, "INIT", "Starting computation of projection matrix.", 7)
if self.POD:
U, s, _ = np.linalg.svd(self.samplingEngine.RPODCoalesced)
s = s ** 2.
else:
Gramian = self.HFEngine.innerProduct(
self.samplingEngine.samplesCoalesced,
self.samplingEngine.samplesCoalesced)
U, s, _ = np.linalg.svd(Gramian)
nsamples = self.samplingEngine.samplesCoalesced.shape[1]
snorm = np.cumsum(s[::-1]) / np.sum(s)
nPODTrunc = min(nsamples - np.argmax(snorm > self.PODTolerance),
self.R)
pMat = self.samplingEngine.samplesCoalesced.dot(U[:, : nPODTrunc])
vbMng(self, "MAIN",
"Assembling {}x{} projection matrix from {} samples.".format(
*(pMat.shape), nsamples), 5)
vbMng(self, "DEL", "Done computing projection matrix.", 7)
return pMat
def _setupAffineBlocks(self):
"""Compute list of marginalized affine blocks of system."""
hasAs = hasattr(self, "AsList") and self.AsList is not None
hasbs = hasattr(self, "bsList") and self.bsList is not None
if hasAs and hasbs: return
vbMng(self, "INIT", "Computing affine blocks of system.", 10)
mu0Eff = self.mu0.data[0, self.directionPivot]
if not hasAs: self.AsList = [None] * len(self.musMarginal)
if not hasbs: self.bsList = [None] * len(self.musMarginal)
for k, muM in enumerate(self.musMarginal):
muEff = [fp] * self.npar
for jj, kk in enumerate(self.directionMarginal):
muEff[kk] = muM(0, jj)
MEnginek = MarginalProxyEngine(self.HFEngine, muEff)
if not hasAs:
self.AsList[k] = MEnginek.affineLinearSystemA(mu0Eff,
self.scaleFactorPivot)
if not hasbs:
self.bsList[k] = MEnginek.affineLinearSystemb(mu0Eff,
self.scaleFactorPivot)
vbMng(self, "DEL", "Done computing affine blocks.", 10)
def setupApprox(self):
"""Compute RB projection matrix."""
if self.checkComputedApprox():
return
RROMPyAssert(self._mode, message = "Cannot setup approximant.")
vbMng(self, "INIT", "Setting up {}.".format(self.name()), 5)
self.computeScaleFactor()
self.computeSnapshots()
self._setupAffineBlocks()
pMat = self._setupProjectionMatrix()
ARBsList, bRBsList = self.assembleReducedSystemMarginal(pMat)
if self.trainedModel is None:
self.trainedModel = tModel()
self.trainedModel.verbosity = self.verbosity
self.trainedModel.timestamp = self.timestamp
data = TrainedModelPivotedData(self.trainedModel.name(), self.mu0,
pMat, self.scaleFactor,
self.HFEngine.rescalingExp,
self.directionPivot)
data.musPivot = copy(self.musPivot)
data.musMarginal = copy(self.musMarginal)
self.trainedModel.data = data
else:
self.trainedModel = self.trainedModel
self.trainedModel.data.projMat = copy(pMat)
self.trainedModel.data.marginalInterp = self._setupMarginalInterp()
self.trainedModel.data.ARBsList = ARBsList
self.trainedModel.data.bRBsList = bRBsList
+ vbMng(self, "INIT", "Matching poles.", 10)
+ self.trainedModel.initializeFromAffine(self.HFEngine,
+ self.matchingWeight, self.POD)
+ 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 assembleReducedSystemMarginal(self, pMat : sampList = None)\
-> Tuple[List[List[Np2D]], List[List[Np1D]]]:
"""Build affine blocks of RB linear system through projections."""
if pMat is None:
self.setupApprox()
ARBsList = self.trainedModel.data.ARBsList
bRBsList = self.trainedModel.data.bRBsList
else:
vbMng(self, "INIT", "Projecting affine terms of HF model.", 10)
ARBsList = [None] * len(self.musMarginal)
bRBsList = [None] * len(self.musMarginal)
for k, (As, bs, sample) in enumerate(zip(self.AsList, self.bsList,
self.samplingEngine.samples)):
compLocal = self.HFEngine.innerProduct(sample, pMat)
projList = sample.dot(customPInv(compLocal))
ARBsList[k], bRBsList[k] = projectAffineDecomposition(As, bs,
projList)
vbMng(self, "DEL", "Done projecting affine terms.", 10)
return ARBsList, bRBsList
diff --git a/rrompy/reduction_methods/pole_matching/__init__.py b/rrompy/reduction_methods/pole_matching/__init__.py
deleted file mode 100644
index 0370cf1..0000000
--- a/rrompy/reduction_methods/pole_matching/__init__.py
+++ /dev/null
@@ -1,29 +0,0 @@
-# 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 .
-#
-
-from .generic_pole_matching_approximant import GenericPoleMatchingApproximant
-from .rational_interpolant_pole_matching import RationalInterpolantPoleMatching
-from .reduced_basis_pole_matching import ReducedBasisPoleMatching
-
-__all__ = [
- 'GenericPoleMatchingApproximant',
- 'RationalInterpolantPoleMatching',
- 'ReducedBasisPoleMatching'
- ]
-
-
diff --git a/rrompy/reduction_methods/pole_matching/generic_pole_matching_approximant.py b/rrompy/reduction_methods/pole_matching/generic_pole_matching_approximant.py
deleted file mode 100644
index 3976dbb..0000000
--- a/rrompy/reduction_methods/pole_matching/generic_pole_matching_approximant.py
+++ /dev/null
@@ -1,170 +0,0 @@
-# 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 .
-#
-
-from numpy import inf
-from rrompy.reduction_methods.pivoted.generic_pivoted_approximant import (
- GenericPivotedApproximant)
-from rrompy.utilities.base.types import ListAny, DictAny, HFEng, paramVal
-from rrompy.utilities.exception_manager import RROMPyException, RROMPyWarning
-
-__all__ = ['GenericPoleMatchingApproximant']
-
-class GenericPoleMatchingApproximant(GenericPivotedApproximant):
- """
- ROM pole matching approximant computation for parametric problems
- (ABSTRACT).
-
- 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;
- - 'polybasisMarginal': type of polynomial basis for marginal
- interpolation; allowed values include 'MONOMIAL', 'CHEBYSHEV'
- and 'LEGENDRE'; defaults to 'MONOMIAL';
- - 'MMarginal': degree of marginal interpolant; defaults to 0;
- - 'radialDirectionalWeightsMarginal': radial basis weights for
- marginal interpolant; defaults to 0, i.e. identity;
- - 'interpRcondMarginal': tolerance for marginal interpolation;
- defaults to None.
- Defaults to empty dict.
- 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;
- - 'polybasisMarginal': type of polynomial basis for marginal
- interpolation;
- - 'MMarginal': degree of marginal interpolant;
- - 'radialDirectionalWeightsMarginal': radial basis weights for
- marginal interpolant;
- - 'interpRcondMarginal': tolerance for marginal interpolation.
- 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.
- 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.
- samplerPivot: Pivot sample point generator.
- SMarginal: Total number of marginal samples current approximant relies
- upon.
- samplerMarginal: Marginal sample point generator.
- polybasisMarginal: Type of polynomial basis for marginal interpolation.
- MMarginal: Degree of marginal interpolant.
- radialDirectionalWeightsMarginal: Radial basis weights for marginal
- interpolant.
- interpRcondMarginal: Tolerance for marginal interpolation.
- 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.
- """
-
- def __init__(self, HFEngine:HFEng, mu0 : paramVal = None,
- directionPivot : ListAny = [0],
- approxParameters : DictAny = {}, verbosity : int = 10,
- timestamp : bool = True):
- self._preInit()
- if len(directionPivot) > 1:
- raise RROMPyException(("Exactly 1 pivot parameter allowed in pole "
- "matching."))
- self._addParametersToList(["matchingWeight", "cutOffTolerance",
- "cutOffType"], [1, inf, "MAGNITUDE"])
- super().__init__(HFEngine = HFEngine, mu0 = mu0,
- directionPivot = directionPivot,
- approxParameters = approxParameters,
- verbosity = verbosity, timestamp = timestamp)
- self._postInit()
-
- @property
- def matchingWeight(self):
- """Value of matchingWeight."""
- return self._matchingWeight
- @matchingWeight.setter
- def matchingWeight(self, matchingWeight):
- self._matchingWeight = matchingWeight
- self._approxParameters["matchingWeight"] = self.matchingWeight
-
- @property
- def cutOffTolerance(self):
- """Value of cutOffTolerance."""
- return self._cutOffTolerance
- @cutOffTolerance.setter
- def cutOffTolerance(self, cutOffTolerance):
- self._cutOffTolerance = cutOffTolerance
- self._approxParameters["cutOffTolerance"] = self.cutOffTolerance
-
- @property
- def cutOffType(self):
- """Value of cutOffType."""
- return self._cutOffType
- @cutOffType.setter
- def cutOffType(self, cutOffType):
- try:
- cutOffType = cutOffType.upper().strip().replace(" ","")
- if cutOffType not in ["MAGNITUDE", "POTENTIAL"]:
- raise RROMPyException("Prescribed cutOffType not recognized.")
- self._cutOffType = cutOffType
- except:
- RROMPyWarning(("Prescribed cutOffType not recognized. Overriding "
- "to 'MAGNITUDE'."))
- self._cutOffType = "MAGNITUDE"
- self._approxParameters["cutOffType"] = self.cutOffType
diff --git a/rrompy/reduction_methods/pole_matching/rational_interpolant_pole_matching.py b/rrompy/reduction_methods/pole_matching/rational_interpolant_pole_matching.py
deleted file mode 100644
index 0f28b76..0000000
--- a/rrompy/reduction_methods/pole_matching/rational_interpolant_pole_matching.py
+++ /dev/null
@@ -1,209 +0,0 @@
-# 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 .
-#
-
-from copy import deepcopy as copy
-import numpy as np
-from rrompy.reduction_methods.pivoted.rational_interpolant_pivoted import \
- RationalInterpolantPivoted
-from .generic_pole_matching_approximant import GenericPoleMatchingApproximant
-from rrompy.utilities.poly_fitting.polynomial import (
- PolynomialInterpolator as PI)
-from rrompy.reduction_methods.trained_model import (TrainedModelPivotedData,
- TrainedModelPoleMatchingRational as tModel)
-from rrompy.utilities.base import verbosityManager as vbMng
-from rrompy.utilities.exception_manager import RROMPyAssert
-
-__all__ = ['RationalInterpolantPoleMatching']
-
-class RationalInterpolantPoleMatching(GenericPoleMatchingApproximant,
- RationalInterpolantPivoted):
- """
- ROM pivoted rational interpolant computation for parametric problems with
- pole matching.
-
- 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;
- - '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;
- - '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.
- 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;
- - 'radialDirectionalWeightsPivot': radial basis weights for pivot
- numerator;
- - 'radialDirectionalWeightsMarginal': radial basis weights for
- marginal interpolant;
- - '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.
- 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.
- radialDirectionalWeightsPivot: Radial basis weights for pivot
- numerator.
- radialDirectionalWeightsMarginal: Radial basis weights for marginal
- interpolant.
- 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 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.computeScaleFactor()
- self.computeSnapshots()
- if self.trainedModel is None:
- self.trainedModel = tModel()
- self.trainedModel.verbosity = self.verbosity
- self.trainedModel.timestamp = self.timestamp
- data = TrainedModelPivotedData(self.trainedModel.name(), self.mu0,
- self.samplingEngine.samplesCoalesced,
- self.scaleFactor,
- self.HFEngine.rescalingExp,
- self.directionPivot)
- data.musPivot = copy(self.musPivot)
- data.musMarginal = copy(self.musMarginal)
- self.trainedModel.data = data
- else:
- self.trainedModel = self.trainedModel
- self.trainedModel.data.projMat = copy(
- self.samplingEngine.samplesCoalesced)
- self.trainedModel.data.marginalInterp = self._setupMarginalInterp()
- if self.N > 0:
- Qs = self._setupDenominator()[0]
- else:
- Q = PI()
- Q.npar = self.musPivot.shape[1]
- Q.coeffs = np.ones(tuple([1] * Q.npar), dtype = np.complex)
- Q.polybasis = self.polybasisPivot0
- Qs = [Q for _ in range(len(self.musMarginal))]
- self.trainedModel.data._temporary = True
- 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)
- 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)
-
diff --git a/rrompy/reduction_methods/pole_matching/reduced_basis_pole_matching.py b/rrompy/reduction_methods/pole_matching/reduced_basis_pole_matching.py
deleted file mode 100644
index 4b83def..0000000
--- a/rrompy/reduction_methods/pole_matching/reduced_basis_pole_matching.py
+++ /dev/null
@@ -1,182 +0,0 @@
-# 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 .
-#
-
-from numpy import isinf
-from copy import deepcopy as copy
-from rrompy.reduction_methods.pivoted.reduced_basis_pivoted import \
- ReducedBasisPivoted
-from .generic_pole_matching_approximant import GenericPoleMatchingApproximant
-from rrompy.reduction_methods.trained_model import (TrainedModelPivotedData,
- TrainedModelPoleMatchingReducedBasis as tModel)
-from rrompy.utilities.base import verbosityManager as vbMng
-from rrompy.utilities.exception_manager import RROMPyAssert
-
-__all__ = ['ReducedBasisPoleMatching']
-
-class ReducedBasisPoleMatching(GenericPoleMatchingApproximant,
- ReducedBasisPivoted):
- """
- ROM pivoted RB approximant computation for parametric problems with pole
- matching.
-
- 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;
- - 'R': rank for pivot Galerkin projection; defaults to prod(S);
- - 'PODTolerance': tolerance for pivot snapshots POD; defaults to
- -1;
- - 'polybasisMarginal': type of polynomial basis for marginal
- interpolation; allowed values include 'MONOMIAL', 'CHEBYSHEV'
- and 'LEGENDRE'; defaults to 'MONOMIAL';
- - 'MMarginal': degree of marginal interpolant; defaults to 0;
- - 'radialDirectionalWeightsMarginal': radial basis weights for
- marginal interpolant; defaults to 0, i.e. identity;
- - 'interpRcondMarginal': tolerance for marginal interpolation;
- defaults to None.
- Defaults to empty dict.
- 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;
- - 'R': rank for Galerkin projection;
- - 'PODTolerance': tolerance for snapshots POD;
- - 'polybasisMarginal': type of polynomial basis for marginal
- interpolation;
- - 'MMarginal': degree of marginal interpolant;
- - 'radialDirectionalWeightsMarginal': radial basis weights for
- marginal interpolant;
- - 'interpRcondMarginal': tolerance for marginal interpolation.
- 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.
- 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.
- samplerPivot: Pivot sample point generator.
- SMarginal: Total number of marginal samples current approximant relies
- upon.
- samplerMarginal: Marginal sample point generator.
- R: Rank for Galerkin projection.
- PODTolerance: Tolerance for pivot snapshots POD.
- polybasisMarginal: Type of polynomial basis for marginal interpolation.
- MMarginal: Degree of marginal interpolant.
- radialDirectionalWeightsMarginal: Radial basis weights for marginal
- interpolant.
- interpRcondMarginal: Tolerance for marginal interpolation.
- 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.
- 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.
- As: List of sparse matrices (in CSC format) representing coefficients
- of linear system matrix.
- bs: List of numpy vectors representing coefficients of linear system
- RHS.
- ARBs: List of sparse matrices (in CSC format) representing coefficients
- of compressed linear system matrix.
- bRBs: List of numpy vectors representing coefficients of compressed
- linear system RHS.
- """
-
- def setupApprox(self):
- """Compute RB projection matrix."""
- if self.checkComputedApprox():
- return
- RROMPyAssert(self._mode, message = "Cannot setup approximant.")
- vbMng(self, "INIT", "Setting up {}.".format(self.name()), 5)
- self.computeScaleFactor()
- self.computeSnapshots()
- self._setupAffineBlocks()
- pMat = self._setupProjectionMatrix()
- ARBsList, bRBsList = self.assembleReducedSystemMarginal(pMat)
- if self.trainedModel is None:
- self.trainedModel = tModel()
- self.trainedModel.verbosity = self.verbosity
- self.trainedModel.timestamp = self.timestamp
- data = TrainedModelPivotedData(self.trainedModel.name(), self.mu0,
- pMat, self.scaleFactor,
- self.HFEngine.rescalingExp,
- self.directionPivot)
- data.musPivot = copy(self.musPivot)
- data.musMarginal = copy(self.musMarginal)
- self.trainedModel.data = data
- else:
- self.trainedModel = self.trainedModel
- self.trainedModel.data.projMat = copy(pMat)
- self.trainedModel.data.marginalInterp = self._setupMarginalInterp()
- self.trainedModel.data.ARBsList = ARBsList
- self.trainedModel.data.bRBsList = bRBsList
- vbMng(self, "INIT", "Matching poles.", 10)
- self.trainedModel.initializeFromAffine(self.HFEngine,
- self.matchingWeight, self.POD)
- vbMng(self, "DEL", "Done matching poles.", 10)
- if not 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)
diff --git a/rrompy/reduction_methods/trained_model/__init__.py b/rrompy/reduction_methods/trained_model/__init__.py
index 4354bc8..562d08a 100644
--- a/rrompy/reduction_methods/trained_model/__init__.py
+++ b/rrompy/reduction_methods/trained_model/__init__.py
@@ -1,42 +1,36 @@
# 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 .
#
from .trained_model_data import TrainedModelData
from .trained_model_rational import TrainedModelRational
from .trained_model_reduced_basis import TrainedModelReducedBasis
from .trained_model_pivoted_data import TrainedModelPivotedData
from .trained_model_pivoted_rational import TrainedModelPivotedRational
from .trained_model_pivoted_reduced_basis import \
TrainedModelPivotedReducedBasis
-from .trained_model_pole_matching_rational import \
- TrainedModelPoleMatchingRational
-from .trained_model_pole_matching_reduced_basis import \
- TrainedModelPoleMatchingReducedBasis
__all__ = [
'TrainedModelData',
'TrainedModelRational',
'TrainedModelReducedBasis',
'TrainedModelPivotedData',
'TrainedModelPivotedRational',
'TrainedModelPivotedReducedBasis',
- 'TrainedModelPoleMatchingRational',
- 'TrainedModelPoleMatchingReducedBasis'
]
diff --git a/rrompy/reduction_methods/trained_model/trained_model_pivoted_general.py b/rrompy/reduction_methods/trained_model/trained_model_pivoted_general.py
index d27a318..f2e7141 100644
--- a/rrompy/reduction_methods/trained_model/trained_model_pivoted_general.py
+++ b/rrompy/reduction_methods/trained_model/trained_model_pivoted_general.py
@@ -1,76 +1,286 @@
# 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 .
#
from abc import abstractmethod
+import numpy as np
from .trained_model import TrainedModel
-from rrompy.utilities.base.types import paramVal, paramList
-from rrompy.parameter import checkParameterList
+from rrompy.utilities.base.types import (Np1D, Tuple, ListAny, paramVal,
+ paramList, sampList, HFEng)
+from rrompy.utilities.base import verbosityManager as vbMng, freepar as fp
+from rrompy.utilities.numerical import pointMatching
+from rrompy.utilities.poly_fitting.heaviside import HeavisideInterpolator as HI
+from rrompy.utilities.exception_manager import RROMPyException, RROMPyAssert
+from rrompy.parameter import checkParameter, checkParameterList
+from rrompy.sampling import emptySampleList
+
__all__ = ['TrainedModelPivotedGeneral']
class TrainedModelPivotedGeneral(TrainedModel):
"""
- ROM approximant evaluation for generic pivoted approximant.
+ ROM approximant evaluation for pivoted approximants (with pole matching).
Attributes:
Data: dictionary with all that can be pickled.
"""
def centerNormalizePivot(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.data.mu0Pivot.
Returns:
Normalized parameter.
"""
mu = checkParameterList(mu, self.data.nparPivot)[0]
if mu0 is None: mu0 = self.data.mu0Pivot
rad = ((mu ** self.data.rescalingExpPivot
- mu0 ** self.data.rescalingExpPivot)
/ self.data.scaleFactorPivot)
return rad
def centerNormalizeMarginal(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.data.mu0Marginal.
Returns:
Normalized parameter.
"""
mu = checkParameterList(mu, self.data.nparMarginal)[0]
if mu0 is None: mu0 = self.data.mu0Marginal
rad = ((mu ** self.data.rescalingExpMarginal
- mu0 ** self.data.rescalingExpMarginal)
/ self.data.scaleFactorMarginal)
return rad
+ def initializeFromLists(self, poles:ListAny, coeffs:ListAny, basis:str,
+ HFEngine:HFEng, matchingWeight : float = 1.,
+ POD : bool = True):
+ """Initialize Heaviside representation."""
+ musM = self.data.musMarginal
+ margAbsDist = np.sum(np.abs(np.repeat(musM.data, len(musM), 0)
+ - np.tile(musM.data, [len(musM), 1])
+ ), axis = 1).reshape(len(musM), len(musM))
+ N = len(poles[0])
+ explored = [0]
+ unexplored = list(range(1, len(musM)))
+ for _ in range(1, len(musM)):
+ minIdx = np.argmin(np.concatenate([margAbsDist[ex, unexplored] \
+ for ex in explored]))
+ minIex = explored[minIdx // len(unexplored)]
+ minIunex = unexplored[minIdx % len(unexplored)]
+ dist = np.abs(np.tile(poles[minIex].reshape(-1, 1), N)
+ - poles[minIunex].reshape(1, -1))
+ if matchingWeight != 0:
+ resex = coeffs[minIex][: N]
+ resunex = coeffs[minIunex][: N]
+ if POD:
+ distR = resex.dot(resunex.T.conj())
+ distR = (distR.T / np.linalg.norm(resex, axis = 1)).T
+ distR = distR / np.linalg.norm(resunex, axis = 1)
+ else:
+ resex = self.data.projMat.dot(resex.T)
+ resunex = self.data.projMat.dot(resunex.T)
+ distR = HFEngine.innerProduct(resex, resunex).T
+ distR = (distR.T / HFEngine.norm(resex)).T
+ distR = distR / HFEngine.norm(resunex)
+ distR = np.abs(distR)
+ distR[distR > 1.] = 1.
+ dist += 2. / np.pi * matchingWeight * np.arccos(distR)
+ reordering = pointMatching(dist)
+ poles[minIunex] = poles[minIunex][reordering]
+ coeffs[minIunex][: N] = coeffs[minIunex][reordering]
+ explored += [minIunex]
+ unexplored.remove(minIunex)
+ HIs = []
+ for pls, cfs in zip(poles, coeffs):
+ hsi = HI()
+ hsi.poles = pls
+ hsi.coeffs = cfs
+ hsi.npar = 1
+ hsi.polybasis = basis
+ HIs += [hsi]
+ self.data.HIs = HIs
+
+ def recompressByCutOff(self, murange : Tuple[float, float] = [- 1., 1.],
+ tol : float = np.inf, rtype : str = "MAGNITUDE"):
+ if np.isinf(tol): return " No poles erased."
+ N = len(self.data.HIs[0].poles)
+ mu0 = np.mean(murange)
+ musig = murange[0] - mu0
+ if np.isclose(musig, 0.):
+ radius = lambda x: np.abs(x - mu0)
+ else:
+ if rtype == "MAGNITUDE":
+ murdir = (murange[0] - mu0) / np.abs(musig)
+ def radius(x):
+ scalprod = (x - mu0) * murdir.conj() / np.abs(musig)
+ rescalepar = np.abs(np.real(scalprod))
+ rescaleort = np.abs(np.imag(scalprod))
+ return ((rescalepar - 1.) ** 2. * (rescalepar > 1.)
+ + rescaleort ** 2.) ** .5
+ else:#if rtype == "POTENTIAL":
+ def radius(x):
+ rescale = (x - mu0) / musig
+ return np.max(np.abs(rescale * np.array([-1., 1.])
+ + (rescale ** 2. - 1) ** .5)) - 1.
+ keepPole, removePole = [], []
+ for j in range(N):
+ for hi in self.data.HIs:
+ if radius(hi.poles[j]) <= tol:
+ keepPole += [j]
+ break
+ if len(keepPole) == 0 or keepPole[-1] != j: removePole += [j]
+ if len(keepPole) == N: return " No poles erased."
+ keepCoeff = keepPole + [N]
+ keepCoeff = keepCoeff + list(range(N + 1,len(self.data.HIs[0].coeffs)))
+ for hi in self.data.HIs:
+ polyCorrection = np.zeros_like(hi.coeffs[0, :])
+ for j in removePole:
+ polyCorrection += hi.coeffs[j, :] / (mu0 - hi.poles[j])
+ if len(hi.coeffs) == N:
+ hi.coeffs = np.vstack((hi.coeffs, polyCorrection))
+ else:
+ hi.coeffs[N, :] += polyCorrection
+ hi.poles = hi.poles[keepPole]
+ hi.coeffs = hi.coeffs[keepCoeff, :]
+ return (" Erased {} pole".format(len(removePole))
+ + "s" * (len(removePole) > 1) + ".")
+
@abstractmethod
def interpolateMarginal(self, *args, **kwargs):
""" Evaluate marginal interpolator at arbitrary marginal parameter."""
pass
+ def interpolateMarginalInterpolator(self, mu : paramVal = []) -> Np1D:
+ """Obtain interpolated approximant interpolator."""
+ mu = checkParameter(mu, self.data.nparMarginal)[0]
+ hsi = HI()
+ hsi.poles = self.interpolateMarginalPoles(mu)
+ hsi.coeffs = self.interpolateMarginalCoeffs(mu)
+ hsi.npar = 1
+ hsi.polybasis = self.data.HIs[0].polybasis
+ return hsi
+
+ def interpolateMarginalPoles(self, mu : paramList = []) -> Np1D:
+ """Obtain interpolated approximant poles."""
+ mu = checkParameterList(mu, self.data.nparMarginal)[0]
+ return self.interpolateMarginal(mu, [hi.poles for hi in self.data.HIs])
+
+ def interpolateMarginalCoeffs(self, mu : paramList = []) -> Np1D:
+ """Obtain interpolated approximant coefficients."""
+ mu = checkParameterList(mu, self.data.nparMarginal)[0]
+ cs = self.interpolateMarginal(mu, [hi.coeffs for hi in self.data.HIs])
+ if isinstance(cs, (list, tuple,)): cs = np.array(cs)
+ return cs.reshape(self.data.HIs[0].coeffs.shape)
+
+ def getApproxReduced(self, mu : paramList = []) -> sampList:
+ """
+ Evaluate reduced representation of approximant at arbitrary parameter.
+
+ Args:
+ mu: Target parameter.
+ """
+ RROMPyAssert(self.data.nparPivot, 1, "Number of pivot parameters")
+ mu = checkParameterList(mu, self.data.npar)[0]
+ if (not hasattr(self, "lastSolvedApproxReduced")
+ or self.lastSolvedApproxReduced != mu):
+ vbMng(self, "INIT",
+ "Evaluating approximant at mu = {}.".format(mu), 12)
+ self.uApproxReduced = emptySampleList()
+ self.uApproxReduced.reset((self.data.HIs[0].coeffs.shape[1],
+ len(mu)), self.data.projMat.dtype)
+
+ for i, muPL in enumerate(mu):
+ muL = self.centerNormalizePivot([muPL(0, x) \
+ for x in self.data.directionPivot])
+ muM = [muPL(0, x) for x in self.data.directionMarginal]
+ vbMng(self, "INIT",
+ "Assembling reduced model for mu = {}.".format(muPL), 87)
+ hsL = self.interpolateMarginalInterpolator(muM)
+ vbMng(self, "DEL", "Done assembling reduced model.", 87)
+ self.uApproxReduced[i] = hsL(muL)
+ vbMng(self, "DEL", "Done evaluating approximant.", 12)
+ self.lastSolvedApproxReduced = mu
+ return self.uApproxReduced
+
+ def getPoles(self, *args, **kwargs) -> Np1D:
+ """
+ Obtain approximant poles.
+
+ Returns:
+ Numpy complex vector of poles.
+ """
+ RROMPyAssert(self.data.nparPivot, 1, "Number of pivot parameters")
+ if len(args) + len(kwargs) > 1:
+ raise RROMPyException(("Wrong number of parameters passed. "
+ "Only 1 available."))
+ elif len(args) + len(kwargs) == 1:
+ if len(args) == 1:
+ mVals = args[0]
+ else:
+ mVals = kwargs["marginalVals"]
+ if not hasattr(mVals, "__len__"): mVals = [mVals]
+ mVals = list(mVals)
+ else:
+ mVals = [fp]
+ try:
+ rDim = mVals.index(fp)
+ if rDim < len(mVals) - 1 and fp in mVals[rDim + 1 :]:
+ raise
+ except:
+ raise RROMPyException(("Exactly 1 'freepar' entry in "
+ "marginalVals must be provided."))
+ if rDim != self.data.directionPivot[0]:
+ raise RROMPyException(("'freepar' entry in marginalVals must "
+ "coincide with pivot direction."))
+ mVals[rDim] = self.data.mu0(rDim)
+ mMarg = [mVals[j] for j in range(len(mVals)) if j != rDim]
+ roots = np.array(self.interpolateMarginalPoles(mMarg))
+ return np.power(self.data.mu0(rDim) ** self.data.rescalingExp[rDim]
+ + self.data.scaleFactor[rDim] * roots,
+ 1. / self.data.rescalingExp[rDim])
+
+ def getResidues(self, *args, **kwargs) -> Np1D:
+ """
+ Obtain approximant residues.
+
+ Returns:
+ Numpy matrix with residues as columns.
+ """
+ pls = self.getPoles(*args, **kwargs)
+ if len(args) == 1:
+ mVals = args[0]
+ else:
+ mVals = kwargs["marginalVals"]
+ if not hasattr(mVals, "__len__"): mVals = [mVals]
+ mVals = list(mVals)
+ rDim = mVals.index(fp)
+ mMarg = [mVals[j] for j in range(len(mVals)) if j != rDim]
+ residues = self.interpolateMarginalCoeffs(mMarg)[: len(pls)]
+ res = self.data.projMat.dot(residues.T)
+ return pls, res
diff --git a/rrompy/reduction_methods/trained_model/trained_model_pivoted_rational.py b/rrompy/reduction_methods/trained_model/trained_model_pivoted_rational.py
index f459de4..da08aea 100644
--- a/rrompy/reduction_methods/trained_model/trained_model_pivoted_rational.py
+++ b/rrompy/reduction_methods/trained_model/trained_model_pivoted_rational.py
@@ -1,189 +1,219 @@
# 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 .
#
import numpy as np
+from scipy.special import factorial as fact
+from itertools import combinations
from .trained_model_pivoted_general import TrainedModelPivotedGeneral
from .trained_model_rational import TrainedModelRational
from rrompy.utilities.base.types import (Np1D, List, ListAny, paramList,
- sampList)
+ sampList, HFEng)
from rrompy.utilities.base import verbosityManager as vbMng, freepar as fp
from rrompy.utilities.poly_fitting.polynomial import polyroots
-from rrompy.utilities.exception_manager import RROMPyException
+from rrompy.utilities.poly_fitting.heaviside import rational2heaviside
+from rrompy.utilities.exception_manager import RROMPyException, RROMPyAssert
from rrompy.parameter import checkParameter, checkParameterList
from rrompy.sampling import sampleList, emptySampleList
__all__ = ['TrainedModelPivotedRational']
class TrainedModelPivotedRational(TrainedModelPivotedGeneral,
TrainedModelRational):
"""
- ROM approximant evaluation for pivoted Rational approximant.
+ ROM approximant evaluation for pivoted Rational approximant (with pole
+ matching).
Attributes:
Data: dictionary with all that can be pickled.
"""
+ def initializeFromRational(self, HFEngine:HFEng,
+ matchingWeight : float = 1., POD : bool = True):
+ """Initialize Heaviside representation."""
+ RROMPyAssert(self.data.nparPivot, 1, "Number of pivot parameters")
+ poles, coeffs = [], []
+ for Q, P in zip(self.data.Qs, self.data.Ps):
+ cfs, pls, basis = rational2heaviside(P, Q)
+ poles += [pls]
+ coeffs += [cfs]
+ self.initializeFromLists(poles, coeffs, basis, HFEngine,
+ matchingWeight, POD)
+ self.data._temporary = False
+
def interpolateMarginal(self, mu : paramList = [],
samples : ListAny = [], der : List[int] = None,
scl : Np1D = None) -> sampList:
"""
Evaluate marginal interpolator at arbitrary marginal parameter.
Args:
mu: Target parameter.
samples: Objects to interpolate.
der(optional): Derivatives to take before evaluation.
"""
mu = checkParameterList(mu, self.data.nparMarginal)[0]
sList = isinstance(samples[0], sampleList)
sEff = [None] * len(samples)
for j in range(len(samples)):
if sList:
sEff[j] = samples[j].data
else:
sEff[j] = samples[j]
vbMng(self, "INIT", "Interpolating marginal at mu = {}.".format(mu),
95)
muC = self.centerNormalizeMarginal(mu)
p = emptySampleList()
p.reset((len(sEff[0]), len(muC)))
p.data[:] = 0.
if len(sEff[0]) > 0:
for mIj, spj in zip(self.data.marginalInterp, sEff):
p = p + spj.reshape(len(sEff[0]), - 1) * mIj(muC, der, scl)
vbMng(self, "DEL", "Done interpolating marginal.", 95)
if not sList:
p = p.data.flatten()
return p
def getPValsPivot(self, mu : paramList = []) -> sampList:
"""
Evaluate rational numerators at arbitrary pivot parameter.
Args:
mu: Target parameter.
"""
mu = checkParameterList(mu, self.data.nparPivot)[0]
vbMng(self, "INIT", "Evaluating numerator at mu = {}.".format(mu), 17)
muC = self.centerNormalizePivot(mu)
pP = [sampleList(P(muC)) for P in self.data.Ps]
vbMng(self, "DEL", "Done evaluating numerator.", 17)
return pP
def getPVal(self, mu : paramList = []) -> sampList:
"""
Evaluate rational numerator at arbitrary parameter.
Args:
mu: Target parameter.
"""
+ if self.data._temporary: return super().getPVal(mu)
+ RROMPyAssert(self.data.nparPivot, 1, "Number of pivot parameters")
mu = checkParameterList(mu, self.data.npar)[0]
- muP = checkParameterList(mu.data[:, self.data.directionPivot],
- self.data.nparPivot)[0]
- muM = checkParameterList(mu.data[:, self.data.directionMarginal],
- self.data.nparMarginal)[0]
- pP = self.getPValsPivot(muP)
- return self.interpolateMarginal(muM, pP)
+ p = emptySampleList()
+ p.reset((len(self.data.HIs[0].coeffs.shape[1]), len(mu)))
+ for i, muPL in enumerate(mu):
+ muL = self.centerNormalizePivot([muPL(0, x) \
+ for x in self.data.directionPivot])
+ muM = [muPL(0, x) for x in self.data.directionMarginal]
+ hsL = self.interpolateMarginalInterpolator(muM)
+ p[i] = hsL(muL) * np.prod(muL(0, 0) - hsL.poles)
+ return p
def getQValsPivot(self, mu:Np1D, der : List[int] = None,
scl : Np1D = None) -> Np1D:
"""
Evaluate rational denominators at arbitrary pivot parameter.
Args:
mu: Target parameter.
der(optional): Derivatives to take before evaluation.
"""
mu = checkParameterList(mu, self.data.nparPivot)[0]
vbMng(self, "INIT", "Evaluating denominator at mu = {}.".format(mu),
17)
muC = self.centerNormalizePivot(mu)
qP = [Q(muC, der, scl).reshape(1, -1) for Q in self.data.Qs]
vbMng(self, "DEL", "Done evaluating denominator.", 17)
return qP
def getQVal(self, mu:Np1D, der : List[int] = None,
scl : Np1D = None) -> Np1D:
"""
Evaluate rational denominator at arbitrary parameter.
Args:
mu: Target parameter.
der(optional): Derivatives to take before evaluation.
"""
+ if self.data._temporary: return super().getQVal(mu, der, scl)
+ RROMPyAssert(self.data.nparPivot, 1, "Number of pivot parameters")
mu = checkParameterList(mu, self.data.npar)[0]
- muP = checkParameterList(mu.data[:, self.data.directionPivot],
- self.data.nparPivot)[0]
+ muP = self.centerNormalizePivot(checkParameterList(
+ mu.data[:, self.data.directionPivot], 1)[0])
muM = checkParameterList(mu.data[:, self.data.directionMarginal],
- self.data.nparMarginal)[0]
+ self.data.npar - 1)[0]
if der is None:
- derP, derM = None, None
+ derP, derM = 0, [0]
else:
- derP = [der[x] for x in self.data.directionPivot]
+ derP = der[self.data.directionPivot[0]]
derM = [der[x] for x in self.data.directionMarginal]
- if scl is None:
- sclP, sclM = None, None
- else:
- sclP = [scl[x] for x in self.data.directionPivot]
- sclM = [scl[x] for x in self.data.directionMarginal]
- qP = self.getQValsPivot(muP, derP, sclP)
- return self.interpolateMarginal(muM, qP, derM, sclM)
+ if np.any(np.array(derM) != 0):
+ raise RROMPyException(("Derivatives of Q with respect to marginal "
+ "parameters not allowed."))
+ sclP = 1 if scl is None else scl[self.data.directionPivot[0]]
+ derVal = np.zeros(len(mu), dtype = self.data.HIs[0].poles.dtype)
+ N = len(self.data.HIs[0].poles)
+ if derP == N: derVal[:] = 1.
+ elif derP >= 0 and derP < N:
+ pls = self.interpolateMarginalPoles(muM).reshape(-1, len(mu)).T
+ plsDist = muP.data.reshape(-1, 1) - pls
+ for terms in combinations(np.arange(N), N - derP):
+ derVal += np.prod(plsDist[:, list(terms)], axis = 1)
+ return sclP ** derP * fact(derP) * derVal
def getPoles(self, *args, **kwargs) -> Np1D:
"""
Obtain approximant poles.
Returns:
Numpy complex vector of poles.
"""
if len(args) + len(kwargs) > 1:
raise RROMPyException(("Wrong number of parameters passed. "
"Only 1 available."))
elif len(args) + len(kwargs) == 1:
if len(args) == 1:
mVals = args[0]
else:
mVals = kwargs["marginalVals"]
if not hasattr(mVals, "__len__"): mVals = [mVals]
mVals = list(mVals)
else:
mVals = [fp]
mValsP = [mVals[x] for x in self.data.directionPivot]
mValsM = [mVals[x] for x in self.data.directionMarginal]
try:
rDim = mValsP.index(fp)
if rDim < len(mValsP) - 1 and fp in mValsP[rDim + 1 :]:
raise
except:
raise RROMPyException(("Exactly 1 'freepar' entry in "
"marginalVals must be provided."))
if fp in mValsM:
raise RROMPyException(("'freepar' entry in marginalVals must be "
"pivot dimension."))
rDimB = self.data.directionPivot[rDim]
Qeff = self.interpolateMarginal(mValsM, [Q.coeffs[:, np.newaxis] \
for Q in self.data.Qs])
Qeff = Qeff.reshape(self.data.Qs[0].coeffs.shape)
mValsP[rDim] = self.data.mu0(rDimB)
mValsP = self.centerNormalizePivot(checkParameter(mValsP))
mValsP = list(mValsP.data.flatten())
mValsP[rDim] = fp
return np.power(self.data.mu0(rDimB) ** self.data.rescalingExp[rDimB]
+ self.data.scaleFactor[rDimB]
* polyroots(Qeff, self.data.Qs[0].polybasis, mValsP),
1. / self.data.rescalingExp[rDimB])
diff --git a/rrompy/reduction_methods/trained_model/trained_model_pivoted_reduced_basis.py b/rrompy/reduction_methods/trained_model/trained_model_pivoted_reduced_basis.py
index 1d1b5fc..e6ac78f 100644
--- a/rrompy/reduction_methods/trained_model/trained_model_pivoted_reduced_basis.py
+++ b/rrompy/reduction_methods/trained_model/trained_model_pivoted_reduced_basis.py
@@ -1,148 +1,168 @@
# 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 .
#
import numpy as np
from .trained_model_pivoted_general import TrainedModelPivotedGeneral
from .trained_model_reduced_basis import TrainedModelReducedBasis
from rrompy.utilities.base.types import (Np1D, ListAny, paramVal, paramList,
- sampList)
+ sampList, HFEng)
from rrompy.utilities.base import verbosityManager as vbMng, freepar as fp
from rrompy.utilities.poly_fitting.polynomial import (
hashIdxToDerivative as hashI)
+from rrompy.utilities.poly_fitting.heaviside import affine2heaviside
from rrompy.utilities.numerical import (eigvalsNonlinearDense,
marginalizePolyList)
-from rrompy.utilities.exception_manager import RROMPyException
+from rrompy.utilities.exception_manager import RROMPyException, RROMPyAssert
from rrompy.parameter import checkParameter, checkParameterList
from rrompy.sampling import emptySampleList, sampleList
__all__ = ['TrainedModelPivotedReducedBasis']
class TrainedModelPivotedReducedBasis(TrainedModelPivotedGeneral,
TrainedModelReducedBasis):
"""
- ROM approximant evaluation for pivoted RB approximant.
+ ROM approximant evaluation for pivoted RB approximant (with pole matching).
Attributes:
Data: dictionary with all that can be pickled.
"""
+ def initializeFromAffine(self, HFEngine:HFEng, matchingWeight : float = 1.,
+ POD : bool = True, jSupp : int = 1):
+ """Initialize Heaviside representation."""
+ RROMPyAssert(self.data.nparPivot, 1, "Number of pivot parameters")
+ poles, coeffs = [], []
+ nbadpls = 0
+ for As, bs in zip(self.data.ARBsList, self.data.bRBsList):
+ cfs, pls, basis = affine2heaviside(As, bs, jSupp)
+ poles += [pls]
+ coeffs += [cfs]
+ nbadpls = max(nbadpls, np.sum(np.isinf(pls)))
+ if nbadpls > 0:
+ for j in range(len(poles)):
+ plsgood = np.argsort(np.abs(pls))[: - nbadpls]
+ poles[j] = poles[j][plsgood]
+ coeffs[j] = coeffs[j][plsgood, :]
+ self.initializeFromLists(poles, coeffs, basis, HFEngine,
+ matchingWeight, POD)
+
def interpolateMarginal(self, mu : paramVal = [],
samples : ListAny = []) -> ListAny:
"""
Evaluate marginal interpolator at arbitrary marginal parameter.
Args:
mu: Target parameter.
samples: Objects to interpolate.
"""
mu = checkParameter(mu, self.data.nparMarginal)
sList = isinstance(samples[0], sampleList)
sEff = [None] * len(samples)
for j in range(len(samples)):
if sList:
sEff[j] = samples[j].data
else:
sEff[j] = samples[j]
try:
dtype = sEff[0].dtype
except:
dtype = sEff[0][0].dtype
vbMng(self, "INIT", "Interpolating marginal at mu = {}.".format(mu),
95)
muC = self.centerNormalizeMarginal(mu)
p = np.zeros_like(sEff[0], dtype = dtype)
for mIj, spj in zip(self.data.marginalInterp, sEff):
p += mIj(muC) * spj
vbMng(self, "DEL", "Done interpolating marginal.", 95)
return p
def getApproxReduced(self, mu : paramList = []) -> sampList:
"""
Evaluate reduced representation of approximant at arbitrary parameter.
Args:
mu: Target parameter.
"""
mu = checkParameterList(mu, self.data.npar)[0]
if (not hasattr(self, "lastSolvedApproxReduced")
or self.lastSolvedApproxReduced != mu):
vbMng(self, "INIT",
"Computing RB solution at mu = {}.".format(mu), 12)
self.uApproxReduced = emptySampleList()
self.uApproxReduced.reset((self.data.ARBsList[0][0].shape[0],
len(mu)), self.data.projMat.dtype)
for i, muPL in enumerate(mu):
muP = checkParameter([muPL(0, x) \
for x in self.data.directionPivot])
muM = [muPL(0, x) for x in self.data.directionMarginal]
vbMng(self, "INIT",
"Assembling reduced model for mu = {}.".format(muPL), 87)
ARBs = self.interpolateMarginal(muM, self.data.ARBsList)
bRBs = self.interpolateMarginal(muM, self.data.bRBsList)
ARBmu = ARBs[0]
bRBmu = bRBs[0]
for j in range(1, max(len(ARBs), len(bRBs))):
theta = np.prod(self.centerNormalizePivot(muP)
** hashI(j, self.data.nparPivot))
if j < len(ARBs):
ARBmu += theta * ARBs[j]
if j < len(bRBs):
bRBmu += theta * bRBs[j]
vbMng(self, "DEL", "Done assembling reduced model.", 87)
vbMng(self, "INIT",
"Solving reduced model for mu = {}.".format(muPL), 75)
self.uApproxReduced[i] = np.linalg.solve(ARBmu, bRBmu)
vbMng(self, "DEL", "Done solving reduced model.", 75)
vbMng(self, "DEL", "Done computing RB solution.", 12)
self.lastSolvedApproxReduced = mu
return self.uApproxReduced
def getPoles(self, marginalVals : ListAny = [fp], jSupp : int = 1,
**kwargs) -> Np1D:
"""
Obtain approximant poles.
Returns:
Numpy complex vector of poles.
"""
if not hasattr(marginalVals, "__len__"): marginalVals = [marginalVals]
marginalVals = list(marginalVals)
try:
rDim = marginalVals.index(fp)
if rDim < len(marginalVals) - 1 and fp in marginalVals[rDim + 1 :]:
raise
except:
raise RROMPyException(("Exactly 1 'freepar' entry in "
"marginalVals must be provided."))
if rDim in self.data.directionMarginal:
raise RROMPyException(("'freepar' entry in marginalVals must be "
"in pivot dimension."))
mValsP = [marginalVals[x] for x in self.data.directionPivot]
rDimP = mValsP.index(fp)
mValsM = [marginalVals[x] for x in self.data.directionMarginal]
ARBs = self.interpolateMarginal(mValsM, self.data.ARBsList)
if self.data.nparPivot > 1:
mValsP[rDimP] = self.data.mu0Pivot(rDimP)
mValsP = self.centerNormalizePivot(mValsP)
mValsP = mValsP.data.flatten()
mValsP[rDimP] = fp
ARBs = marginalizePolyList(ARBs, mValsP, "auto")
ev = eigvalsNonlinearDense(ARBs, jSupp = jSupp, **kwargs)
return np.power(self.data.mu0(rDim) ** self.data.rescalingExp[rDim]
+ ev * self.data.scaleFactor[rDim],
1. / self.data.rescalingExp[rDim])
diff --git a/rrompy/reduction_methods/trained_model/trained_model_pole_matching_general.py b/rrompy/reduction_methods/trained_model/trained_model_pole_matching_general.py
deleted file mode 100644
index dadd780..0000000
--- a/rrompy/reduction_methods/trained_model/trained_model_pole_matching_general.py
+++ /dev/null
@@ -1,241 +0,0 @@
-# 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 .
-#
-
-import numpy as np
-from .trained_model_pivoted_general import TrainedModelPivotedGeneral
-from rrompy.utilities.base.types import (Np1D, Tuple, ListAny, paramVal,
- paramList, sampList, HFEng)
-from rrompy.utilities.base import verbosityManager as vbMng, freepar as fp
-from rrompy.utilities.numerical import pointMatching
-from rrompy.utilities.poly_fitting.heaviside import HeavisideInterpolator as HI
-from rrompy.utilities.exception_manager import RROMPyException, RROMPyAssert
-from rrompy.parameter import checkParameter, checkParameterList
-from rrompy.sampling import emptySampleList
-
-__all__ = ['TrainedModelPoleMatchingGeneral']
-
-class TrainedModelPoleMatchingGeneral(TrainedModelPivotedGeneral):
- """
- ROM approximant evaluation for pivoted approximants with pole matching.
-
- Attributes:
- Data: dictionary with all that can be pickled.
- """
-
- def initializeFromLists(self, poles:ListAny, coeffs:ListAny, basis:str,
- HFEngine:HFEng, matchingWeight : float = 1.,
- POD : bool = True):
- """Initialize Heaviside representation."""
- musM = self.data.musMarginal
- margAbsDist = np.sum(np.abs(np.repeat(musM.data, len(musM), 0)
- - np.tile(musM.data, [len(musM), 1])
- ), axis = 1).reshape(len(musM), len(musM))
- N = len(poles[0])
- explored = [0]
- unexplored = list(range(1, len(musM)))
- for _ in range(1, len(musM)):
- minIdx = np.argmin(np.concatenate([margAbsDist[ex, unexplored] \
- for ex in explored]))
- minIex = explored[minIdx // len(unexplored)]
- minIunex = unexplored[minIdx % len(unexplored)]
- dist = np.abs(np.tile(poles[minIex].reshape(-1, 1), N)
- - poles[minIunex].reshape(1, -1))
- if matchingWeight != 0:
- resex = coeffs[minIex][: N]
- resunex = coeffs[minIunex][: N]
- if POD:
- distR = resex.dot(resunex.T.conj())
- distR = (distR.T / np.linalg.norm(resex, axis = 1)).T
- distR = distR / np.linalg.norm(resunex, axis = 1)
- else:
- resex = self.data.projMat.dot(resex.T)
- resunex = self.data.projMat.dot(resunex.T)
- distR = HFEngine.innerProduct(resex, resunex).T
- distR = (distR.T / HFEngine.norm(resex)).T
- distR = distR / HFEngine.norm(resunex)
- distR = np.abs(distR)
- distR[distR > 1.] = 1.
- dist += 2. / np.pi * matchingWeight * np.arccos(distR)
- reordering = pointMatching(dist)
- poles[minIunex] = poles[minIunex][reordering]
- coeffs[minIunex][: N] = coeffs[minIunex][reordering]
- explored += [minIunex]
- unexplored.remove(minIunex)
- HIs = []
- for pls, cfs in zip(poles, coeffs):
- hsi = HI()
- hsi.poles = pls
- hsi.coeffs = cfs
- hsi.npar = 1
- hsi.polybasis = basis
- HIs += [hsi]
- self.data.HIs = HIs
-
- def recompressByCutOff(self, murange : Tuple[float, float] = [- 1., 1.],
- tol : float = np.inf, rtype : str = "MAGNITUDE"):
- if np.isinf(tol): return " No poles erased."
- N = len(self.data.HIs[0].poles)
- mu0 = np.mean(murange)
- musig = murange[0] - mu0
- if np.isclose(musig, 0.):
- radius = lambda x: np.abs(x - mu0)
- else:
- if rtype == "MAGNITUDE":
- murdir = (murange[0] - mu0) / np.abs(musig)
- def radius(x):
- scalprod = (x - mu0) * murdir.conj() / np.abs(musig)
- rescalepar = np.abs(np.real(scalprod))
- rescaleort = np.abs(np.imag(scalprod))
- return ((rescalepar - 1.) ** 2. * (rescalepar > 1.)
- + rescaleort ** 2.) ** .5
- else:#if rtype == "POTENTIAL":
- def radius(x):
- rescale = (x - mu0) / musig
- return np.max(np.abs(rescale * np.array([-1., 1.])
- + (rescale ** 2. - 1) ** .5)) - 1.
- keepPole, removePole = [], []
- for j in range(N):
- for hi in self.data.HIs:
- if radius(hi.poles[j]) <= tol:
- keepPole += [j]
- break
- if len(keepPole) == 0 or keepPole[-1] != j: removePole += [j]
- if len(keepPole) == N: return " No poles erased."
- keepCoeff = keepPole + [N]
- keepCoeff = keepCoeff + list(range(N + 1,len(self.data.HIs[0].coeffs)))
- for hi in self.data.HIs:
- polyCorrection = np.zeros_like(hi.coeffs[0, :])
- for j in removePole:
- polyCorrection += hi.coeffs[j, :] / (mu0 - hi.poles[j])
- if len(hi.coeffs) == N:
- hi.coeffs = np.vstack((hi.coeffs, polyCorrection))
- else:
- hi.coeffs[N, :] += polyCorrection
- hi.poles = hi.poles[keepPole]
- hi.coeffs = hi.coeffs[keepCoeff, :]
- return (" Erased {} pole".format(len(removePole))
- + "s" * (len(removePole) > 1) + ".")
-
- def interpolateMarginalInterpolator(self, mu : paramVal = []) -> Np1D:
- """Obtain interpolated approximant interpolator."""
- mu = checkParameter(mu, self.data.nparMarginal)[0]
- hsi = HI()
- hsi.poles = self.interpolateMarginalPoles(mu)
- hsi.coeffs = self.interpolateMarginalCoeffs(mu)
- hsi.npar = 1
- hsi.polybasis = self.data.HIs[0].polybasis
- return hsi
-
- def interpolateMarginalPoles(self, mu : paramList = []) -> Np1D:
- """Obtain interpolated approximant poles."""
- mu = checkParameterList(mu, self.data.nparMarginal)[0]
- return self.interpolateMarginal(mu, [hi.poles for hi in self.data.HIs])
-
- def interpolateMarginalCoeffs(self, mu : paramList = []) -> Np1D:
- """Obtain interpolated approximant coefficients."""
- mu = checkParameterList(mu, self.data.nparMarginal)[0]
- cs = self.interpolateMarginal(mu, [hi.coeffs for hi in self.data.HIs])
- if isinstance(cs, (list, tuple,)): cs = np.array(cs)
- return cs.reshape(self.data.HIs[0].coeffs.shape)
-
- def getApproxReduced(self, mu : paramList = []) -> sampList:
- """
- Evaluate reduced representation of approximant at arbitrary parameter.
-
- Args:
- mu: Target parameter.
- """
- RROMPyAssert(self.data.nparPivot, 1, "Number of pivot parameters")
- mu = checkParameterList(mu, self.data.npar)[0]
- if (not hasattr(self, "lastSolvedApproxReduced")
- or self.lastSolvedApproxReduced != mu):
- vbMng(self, "INIT",
- "Evaluating approximant at mu = {}.".format(mu), 12)
- self.uApproxReduced = emptySampleList()
- self.uApproxReduced.reset((self.data.HIs[0].coeffs.shape[1],
- len(mu)), self.data.projMat.dtype)
-
- for i, muPL in enumerate(mu):
- muL = self.centerNormalizePivot([muPL(0, x) \
- for x in self.data.directionPivot])
- muM = [muPL(0, x) for x in self.data.directionMarginal]
- vbMng(self, "INIT",
- "Assembling reduced model for mu = {}.".format(muPL), 87)
- hsL = self.interpolateMarginalInterpolator(muM)
- vbMng(self, "DEL", "Done assembling reduced model.", 87)
- self.uApproxReduced[i] = hsL(muL)
- vbMng(self, "DEL", "Done evaluating approximant.", 12)
- self.lastSolvedApproxReduced = mu
- return self.uApproxReduced
-
- def getPoles(self, *args, **kwargs) -> Np1D:
- """
- Obtain approximant poles.
-
- Returns:
- Numpy complex vector of poles.
- """
- RROMPyAssert(self.data.nparPivot, 1, "Number of pivot parameters")
- if len(args) + len(kwargs) > 1:
- raise RROMPyException(("Wrong number of parameters passed. "
- "Only 1 available."))
- elif len(args) + len(kwargs) == 1:
- if len(args) == 1:
- mVals = args[0]
- else:
- mVals = kwargs["marginalVals"]
- if not hasattr(mVals, "__len__"): mVals = [mVals]
- mVals = list(mVals)
- else:
- mVals = [fp]
- try:
- rDim = mVals.index(fp)
- if rDim < len(mVals) - 1 and fp in mVals[rDim + 1 :]:
- raise
- except:
- raise RROMPyException(("Exactly 1 'freepar' entry in "
- "marginalVals must be provided."))
- if rDim != self.data.directionPivot[0]:
- raise RROMPyException(("'freepar' entry in marginalVals must "
- "coincide with pivot direction."))
- mVals[rDim] = self.data.mu0(rDim)
- mMarg = [mVals[j] for j in range(len(mVals)) if j != rDim]
- roots = np.array(self.interpolateMarginalPoles(mMarg))
- return np.power(self.data.mu0(rDim) ** self.data.rescalingExp[rDim]
- + self.data.scaleFactor[rDim] * roots,
- 1. / self.data.rescalingExp[rDim])
-
- def getResidues(self, *args, **kwargs) -> Np1D:
- """
- Obtain approximant residues.
-
- Returns:
- Numpy matrix with residues as columns.
- """
- pls = self.getPoles(*args, **kwargs)
- if len(args) == 1:
- mVals = args[0]
- else:
- mVals = kwargs["marginalVals"]
- if not hasattr(mVals, "__len__"): mVals = [mVals]
- mVals = list(mVals)
- rDim = mVals.index(fp)
- mMarg = [mVals[j] for j in range(len(mVals)) if j != rDim]
- residues = self.interpolateMarginalCoeffs(mMarg)[: len(pls)]
- res = self.data.projMat.dot(residues.T)
- return pls, res
diff --git a/rrompy/reduction_methods/trained_model/trained_model_pole_matching_rational.py b/rrompy/reduction_methods/trained_model/trained_model_pole_matching_rational.py
deleted file mode 100644
index c24171f..0000000
--- a/rrompy/reduction_methods/trained_model/trained_model_pole_matching_rational.py
+++ /dev/null
@@ -1,109 +0,0 @@
-# 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 .
-#
-
-import numpy as np
-from scipy.special import factorial as fact
-from itertools import combinations
-from .trained_model_pole_matching_general import \
- TrainedModelPoleMatchingGeneral
-from .trained_model_pivoted_rational import TrainedModelPivotedRational
-from rrompy.utilities.base.types import Np1D, List, paramList, sampList, HFEng
-from rrompy.utilities.poly_fitting.heaviside import rational2heaviside
-from rrompy.utilities.exception_manager import RROMPyException, RROMPyAssert
-from rrompy.parameter import checkParameterList
-from rrompy.sampling import emptySampleList
-
-__all__ = ['TrainedModelPoleMatchingRational']
-
-class TrainedModelPoleMatchingRational(TrainedModelPoleMatchingGeneral,
- TrainedModelPivotedRational):
- """
- ROM approximant evaluation for pivoted Rational approximant with pole
- matching.
-
- Attributes:
- Data: dictionary with all that can be pickled.
- """
-
- def initializeFromRational(self, HFEngine:HFEng,
- matchingWeight : float = 1., POD : bool = True):
- """Initialize Heaviside representation."""
- RROMPyAssert(self.data.nparPivot, 1, "Number of pivot parameters")
- poles, coeffs = [], []
- for Q, P in zip(self.data.Qs, self.data.Ps):
- cfs, pls, basis = rational2heaviside(P, Q)
- poles += [pls]
- coeffs += [cfs]
- self.initializeFromLists(poles, coeffs, basis, HFEngine,
- matchingWeight, POD)
- self.data._temporary = False
-
- def getPVal(self, mu : paramList = []) -> sampList:
- """
- Evaluate rational numerator at arbitrary parameter.
-
- Args:
- mu: Target parameter.
- """
- if self.data._temporary: return super().getPVal(mu)
- RROMPyAssert(self.data.nparPivot, 1, "Number of pivot parameters")
- mu = checkParameterList(mu, self.data.npar)[0]
- p = emptySampleList()
- p.reset((len(self.data.HIs[0].coeffs.shape[1]), len(mu)))
- for i, muPL in enumerate(mu):
- muL = self.centerNormalizePivot([muPL(0, x) \
- for x in self.data.directionPivot])
- muM = [muPL(0, x) for x in self.data.directionMarginal]
- hsL = self.interpolateMarginalInterpolator(muM)
- p[i] = hsL(muL) * np.prod(muL(0, 0) - hsL.poles)
- return p
-
- def getQVal(self, mu:Np1D, der : List[int] = None,
- scl : Np1D = None) -> Np1D:
- """
- Evaluate rational denominator at arbitrary parameter.
-
- Args:
- mu: Target parameter.
- der(optional): Derivatives to take before evaluation.
- """
- if self.data._temporary: return super().getQVal(mu, der, scl)
- RROMPyAssert(self.data.nparPivot, 1, "Number of pivot parameters")
- mu = checkParameterList(mu, self.data.npar)[0]
- muP = self.centerNormalizePivot(checkParameterList(
- mu.data[:, self.data.directionPivot], 1)[0])
- muM = checkParameterList(mu.data[:, self.data.directionMarginal],
- self.data.npar - 1)[0]
- if der is None:
- derP, derM = 0, [0]
- else:
- derP = der[self.data.directionPivot[0]]
- derM = [der[x] for x in self.data.directionMarginal]
- if np.any(np.array(derM) != 0):
- raise RROMPyException(("Derivatives of Q with respect to marginal "
- "parameters not allowed."))
- sclP = 1 if scl is None else scl[self.data.directionPivot[0]]
- derVal = np.zeros(len(mu), dtype = self.data.HIs[0].poles.dtype)
- N = len(self.data.HIs[0].poles)
- if derP == N: derVal[:] = 1.
- elif derP >= 0 and derP < N:
- pls = self.interpolateMarginalPoles(muM).reshape(-1, len(mu)).T
- plsDist = muP.data.reshape(-1, 1) - pls
- for terms in combinations(np.arange(N), N - derP):
- derVal += np.prod(plsDist[:, list(terms)], axis = 1)
- return sclP ** derP * fact(derP) * derVal
diff --git a/rrompy/reduction_methods/trained_model/trained_model_pole_matching_reduced_basis.py b/rrompy/reduction_methods/trained_model/trained_model_pole_matching_reduced_basis.py
deleted file mode 100644
index 23802b1..0000000
--- a/rrompy/reduction_methods/trained_model/trained_model_pole_matching_reduced_basis.py
+++ /dev/null
@@ -1,56 +0,0 @@
-# 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 .
-#
-
-import numpy as np
-from .trained_model_pole_matching_general import \
- TrainedModelPoleMatchingGeneral
-from .trained_model_pivoted_reduced_basis import \
- TrainedModelPivotedReducedBasis
-from rrompy.utilities.base.types import HFEng
-from rrompy.utilities.poly_fitting.heaviside import affine2heaviside
-from rrompy.utilities.exception_manager import RROMPyAssert
-
-__all__ = ['TrainedModelPoleMatchingReducedBasis']
-
-class TrainedModelPoleMatchingReducedBasis(TrainedModelPoleMatchingGeneral,
- TrainedModelPivotedReducedBasis):
- """
- ROM approximant evaluation for pivoted RB approximant with pole matching.
-
- Attributes:
- Data: dictionary with all that can be pickled.
- """
-
- def initializeFromAffine(self, HFEngine:HFEng, matchingWeight : float = 1.,
- POD : bool = True, jSupp : int = 1):
- """Initialize Heaviside representation."""
- RROMPyAssert(self.data.nparPivot, 1, "Number of pivot parameters")
- poles, coeffs = [], []
- nbadpls = 0
- for As, bs in zip(self.data.ARBsList, self.data.bRBsList):
- cfs, pls, basis = affine2heaviside(As, bs, jSupp)
- poles += [pls]
- coeffs += [cfs]
- nbadpls = max(nbadpls, np.sum(np.isinf(pls)))
- if nbadpls > 0:
- for j in range(len(poles)):
- plsgood = np.argsort(np.abs(pls))[: - nbadpls]
- poles[j] = poles[j][plsgood]
- coeffs[j] = coeffs[j][plsgood, :]
- self.initializeFromLists(poles, coeffs, basis, HFEngine,
- matchingWeight, POD)