diff --git a/examples/2d/pod/fracture_pod.py b/examples/2d/pod/fracture_pod.py
index 40239c9..35463e0 100644
--- a/examples/2d/pod/fracture_pod.py
+++ b/examples/2d/pod/fracture_pod.py
@@ -1,261 +1,261 @@
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.parameter.parameter_sampling import (QuadratureSampler as QS,
QuadratureSamplerTotal as QST,
ManualSampler as MS,
RandomSampler as RS)
-verb = 70
+verb = 10
size = 2
show_sample = True
show_norm = True
Delta = 0
MN = 5
R = (MN + 2) * (MN + 1) // 2
STensorized = (MN + 1) ** 2
PODTol = 1e-6
SPivot = MN + 1
-SMarginal = 3
+SMarginal = 4
MMarginal = SMarginal - 1
matchingWeight = 1.
cutOffTolerance = 1. * np.inf
cutOffType = "MAGNITUDE"
samples = "centered"
samples = "standard"
-samples = "pivoted"
+#samples = "pivoted"
algo = "rational"
-algo = "RB"
+#algo = "RB"
sampling = "quadrature"
-sampling = "quadrature_total"
+#sampling = "quadrature_total"
#sampling = "random"
samplingM = "quadrature"
#samplingM = "quadrature_total"
#samplingM = "random"
polydegreetype = "TOTAL"
-#polydegreetype = "FULL"
+polydegreetype = "FULL"
if samples in ["standard", "pivoted"]:
radial = ""
# radial = "_gaussian"
# radial = "_thinplate"
# radial = "_multiquadric"
rW0 = 5.
radialWeight = [rW0]
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':R, 'POD':True,
'polydegreetype':polydegreetype}
if samples == "standard":
params['polybasis'] = "CHEBYSHEV" + radial
# params['polybasis'] = "LEGENDRE" + radial
# params['polybasis'] = "MONOMIAL" + radial
- params['radialDirectionalWeights'] = radialWeight
+ params['radialDirectionalWeights'] = radialWeight * 2
method = RI
elif samples == "centered":
params['polybasis'] = "MONOMIAL"
method = RI
elif samples == "pivoted":
params['S'] = SPivot
params['SMarginal'] = SMarginal
params['MMarginal'] = MMarginal
params['polybasisPivot'] = "CHEBYSHEV" + radial
params['polybasisMarginal'] = "MONOMIAL" + radialM
params['radialDirectionalWeightsPivot'] = radialWeight
params['radialDirectionalWeightsMarginal'] = radialWeightM
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':R,
'POD':True, 'PODTolerance':PODTol}
if samples == "standard":
method = RB
elif samples == "centered":
method = RB
elif samples == "pivoted":
params['S'] = SPivot
params['R'] = SPivot
params['SMarginal'] = SMarginal
params['MMarginal'] = MMarginal
params['polybasisMarginal'] = "MONOMIAL" + radialM
params['radialDirectionalWeightsMarginal'] = radialWeightM
params['matchingWeight'] = matchingWeight
params['cutOffTolerance'] = cutOffTolerance
params["cutOffType"] = cutOffType
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'] = STensorized
elif sampling == "quadrature_total":
params['sampler'] = QST(murange, "CHEBYSHEV", scalingExp = rescalingExp)
else: # if sampling == "random":
params['sampler'] = RS(murange, "HALTON", scalingExp = rescalingExp)
elif samples == "centered":
params['sampler'] = MS(murange, points = [mu0], scalingExp = rescalingExp)
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")
params['S'] = STensorized
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 != "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,
+ murange, murangeEff, approx, mu0, 25,
[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 662b2c8..8c3b8cb 100644
--- a/examples/2d/pod/matrix_passive_pod.py
+++ b/examples/2d/pod/matrix_passive_pod.py
@@ -1,192 +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.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 = 5
R = (MN + 2) * (MN + 1) // 2
STensorized = (MN + 1) ** 2
PODTol = 1e-6
SPivot = MN + 1
SMarginal = 3
MMarginal = SMarginal - 1
samples = "centered"
samples = "standard"
samples = "pivoted"
algo = "rational"
algo = "RB"
sampling = "quadrature"
#sampling = "quadrature_total"
#sampling = "random"
samplingM = "quadrature"
#samplingM = "quadrature_total"
#samplingM = "random"
if samples in ["standard", "pivoted"]:
radial = ""
# radial = "_gaussian"
# radial = "_thinplate"
# radial = "_multiquadric"
rW0 = 10.
radialWeight = [rW0]
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':R, 'POD':True}
if samples == "standard":
params['polybasis'] = "CHEBYSHEV" + radial
# params['polybasis'] = "LEGENDRE" + radial
# params['polybasis'] = "MONOMIAL" + radial
- params['radialDirectionalWeights'] = radialWeight
+ params['radialDirectionalWeights'] = radialWeight * 2
method = RI
elif samples == "centered":
params['polybasis'] = "MONOMIAL"
method = RI
elif samples == "pivoted":
params['S'] = SPivot
params['SMarginal'] = SMarginal
params['MMarginal'] = MMarginal
params['polybasisPivot'] = "CHEBYSHEV" + radial
params['polybasisMarginal'] = "MONOMIAL" + radialM
params['radialDirectionalWeightsPivot'] = radialWeight
params['radialDirectionalWeightsMarginal'] = radialWeightM
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':R,
'POD':True, 'PODTolerance':PODTol}
if samples == "standard":
method = RB
elif samples == "centered":
method = RB
elif samples == "pivoted":
params['R'] = SPivot
params['S'] = SPivot
params['SMarginal'] = SMarginal
params['MMarginal'] = MMarginal
params['polybasisMarginal'] = "MONOMIAL" + radialM
params['radialDirectionalWeightsMarginal'] = radialWeightM
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")
params["S"] = STensorized
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 == "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 != "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/plot_inf_set.py b/examples/2d/pod/plot_inf_set.py
index 3f4ea96..00aa698 100644
--- a/examples/2d/pod/plot_inf_set.py
+++ b/examples/2d/pod/plot_inf_set.py
@@ -1,336 +1,362 @@
import warnings
import numpy as np
from matplotlib import pyplot as plt
def plotInfSet1FromData(mus, Z, T, R, E, beta, murange, approx, mu0,
exp = 2., normalizeDen = False):
if hasattr(approx, "mus"):
mu2x = approx.mus(0) ** exp
else:
mu2x = mu0[0] ** exp
murangeExp = [[murange[0][0] ** exp], [murange[1][0] ** exp]]
mu1 = np.real(np.power(mus, exp))
ZTmin, ZTmax = min(np.min(Z), np.min(T)), max(np.max(Z), np.max(T))
Rmin, Rmax = np.min(R), np.max(R)
Emin, Emax = np.min(E), np.max(E)
if not np.isnan(beta[0]):
eta = R / beta / E
betamin, betamax = np.min(beta), np.max(beta)
etamin, etamax = np.min(eta), np.max(eta)
plt.figure(figsize = (15, 7))
plt.jet()
plt.semilogy(mu1, Z)
plt.semilogy(mu1, T, '--')
for l_ in approx.trainedModel.getPoles():
plt.plot([np.real(l_ ** exp)] * 2, [ZTmin, ZTmax], 'b:')
plt.plot(np.real(mu2x), [ZTmin] * len(mu2x), 'kx')
+ for j, x in enumerate(np.real(mu2x)):
+ plt.annotate("{}".format(j), (x, ZTmin))
plt.plot([murangeExp[0][0]] * 2, [ZTmin, ZTmax], 'm:')
plt.plot([murangeExp[1][0]] * 2, [ZTmin, ZTmax], 'm:')
plt.xlim(mu1[0], mu1[-1])
plt.title("|u(mu)|, |u_app(mu)|")
plt.grid()
plt.show()
plt.figure(figsize = (15, 7))
plt.jet()
plt.semilogy(mu1, R)
for l_ in approx.trainedModel.getPoles():
plt.plot([np.real(l_ ** exp)] * 2, [Rmin, Rmax], 'b:')
plt.plot(np.real(mu2x), [Rmax] * len(mu2x), 'kx')
+ for j, x in enumerate(np.real(mu2x)):
+ plt.annotate("{}".format(j), (x, Rmax))
plt.plot([murangeExp[0][0]] * 2, [Rmin, Rmax], 'm:')
plt.plot([murangeExp[1][0]] * 2, [Rmin, Rmax], 'm:')
plt.xlim(mu1[0], mu1[-1])
if normalizeDen:
plt.title("|Q(mu)res(mu)|")
else:
plt.title("|res(mu)|")
plt.grid()
plt.show()
plt.figure(figsize = (15, 7))
plt.jet()
plt.semilogy(mu1, E)
for l_ in approx.trainedModel.getPoles():
plt.plot([np.real(l_ ** exp)] * 2, [Emin, Emax], 'b:')
plt.plot(np.real(mu2x), [Emax] * len(mu2x), 'kx')
+ for j, x in enumerate(np.real(mu2x)):
+ plt.annotate("{}".format(j), (x, Emax))
plt.plot([murangeExp[0][0]] * 2, [Emin, Emax], 'm:')
plt.plot([murangeExp[1][0]] * 2, [Emin, Emax], 'm:')
plt.xlim(mu1[0], mu1[-1])
if normalizeDen:
plt.title("|Q(mu)err(mu)|")
else:
plt.title("|err(mu)|")
plt.grid()
plt.show()
if not np.isnan(beta[0]):
plt.figure(figsize = (15, 7))
plt.jet()
plt.plot(mu1, beta)
for l_ in approx.trainedModel.getPoles():
plt.plot([np.real(l_ ** exp)] * 2, [betamin, betamax], 'b:')
plt.plot(np.real(mu2x), [betamax] * len(mu2x), 'kx')
+ for j, x in enumerate(np.real(mu2x)):
+ plt.annotate("{}".format(j), (x, betamax))
plt.plot([murangeExp[0][0]] * 2, [betamin, betamax], 'm:')
plt.plot([murangeExp[1][0]] * 2, [betamin, betamax], 'm:')
plt.xlim(mu1[0], mu1[-1])
plt.title("beta(mu)")
plt.grid()
plt.show()
plt.figure(figsize = (15, 7))
plt.jet()
plt.semilogy(mu1, R / beta)
plt.semilogy(mu1, E, '--')
for l_ in approx.trainedModel.getPoles():
plt.plot([np.real(l_ ** exp)] * 2, [Emin, Emax], 'b:')
plt.plot(np.real(mu2x), [Emax] * len(mu2x), 'kx')
+ for j, x in enumerate(np.real(mu2x)):
+ plt.annotate("{}".format(j), (x, Emax))
plt.plot([murangeExp[0][0]] * 2, [Emin, Emax], 'm:')
plt.plot([murangeExp[1][0]] * 2, [Emin, Emax], 'm:')
plt.xlim(mu1[0], mu1[-1])
if normalizeDen:
plt.title("|Q(mu)res(mu)/beta(mu)|")
else:
plt.title("|res(mu)/beta(mu)|")
plt.grid()
plt.show()
plt.figure(figsize = (15, 7))
plt.jet()
plt.semilogy(mu1, eta)
for l_ in approx.trainedModel.getPoles():
plt.plot([np.real(l_ ** exp)] * 2, [etamin, etamax], 'b:')
plt.plot(np.real(mu2x), [etamax] * len(mu2x), 'kx')
+ for j, x in enumerate(np.real(mu2x)):
+ plt.annotate("{}".format(j), (x, etamax))
plt.plot([murangeExp[0][0]] * 2, [etamin, etamax], 'm:')
plt.plot([murangeExp[1][0]] * 2, [etamin, etamax], 'm:')
plt.xlim(mu1[0], mu1[-1])
plt.title("eta(mu)")
plt.grid()
plt.show()
def plotInfSet1(murange, murangeEff, approx, mu0, nSamples = 20, exp = 2.,
relative = True, normalizeDen = False, nobeta = False):
mu1 = np.linspace(murangeEff[0][0] ** exp, murangeEff[1][0] ** exp,
nSamples)
mus = np.power(mu1, 1. / exp)
Z = approx.normHF(mus)
T = approx.normApprox(mus)
R = approx.normRes(mus)
E = approx.normErr(mus)
if relative:
F = approx.normRHS(mus)
R /= F
E /= Z
if normalizeDen:
Qvals = np.abs(approx.trainedModel.getQVal(mus))
R *= Qvals
E *= Qvals
if nobeta:
beta = np.empty(len(mus))
beta[:] = np.nan
else:
beta = approx.HFEngine.stabilityFactor(approx.getHF(mus), mus)
plotInfSet1FromData(mus, Z, T, R, E, beta, murange, approx, mu0,
exp, normalizeDen)
return mus, Z, T, R, E, beta
def plotInfSet2FromData(mus, Ze, Te, Re, Ee, beta, murange, approx, mu0,
exps = [2., 2.], clip = -1, normalizeDen = False):
if hasattr(approx, "mus"):
mu2x, mu2y = approx.mus(0) ** exps[0], approx.mus(1) ** exps[1]
else:
mu2x, mu2y = mu0[0] ** exps[0], mu0[1] ** exps[1]
murangeExp = [[murange[0][0] ** exps[0], murange[0][1] ** exps[1]],
[murange[1][0] ** exps[0], murange[1][1] ** exps[1]]]
mu1s = np.unique([m[0] for m in mus])
mu2s = np.unique([m[1] for m in mus])
mu1 = np.power(mu1s, exps[0])
mu2 = np.power(mu2s, exps[1])
Mu1, Mu2 = np.meshgrid(np.real(mu1), np.real(mu2))
Z = np.log10(Ze)
T = np.log10(Te)
R = np.log10(Re)
E = np.log10(Ee)
ZTmin, ZTmax = min(np.min(Z), np.min(T)), max(np.max(Z), np.max(T))
Rmin, Rmax = np.min(R), np.max(R)
Emin, Emax = np.min(E), np.max(E)
if not np.isnan(beta[0, 0]):
betamin, betamax = np.min(beta), np.max(beta)
if clip > 0:
ZTmax -= clip * (ZTmax - ZTmin)
cmap = plt.cm.bone
else:
cmap = plt.cm.jet
warnings.simplefilter("ignore", category = (UserWarning,
np.ComplexWarning))
plt.figure(figsize = (15, 7))
plt.jet()
p = plt.contourf(Mu1, Mu2, Z, cmap = cmap,
levels = np.linspace(ZTmin, ZTmax, 50))
if clip > 0:
plt.contour(Mu1, Mu2, Z, [ZTmin])
plt.plot(np.real(mu2x), np.real(mu2y), 'kx')
+ for j, xy in enumerate(zip(np.real(mu2x), np.real(mu2y))):
+ plt.annotate("{}".format(j), xy)
plt.plot([murangeExp[0][0]] * 2,
[murangeExp[0][1], murangeExp[1][1]], 'm:')
plt.plot([murangeExp[0][0], murangeExp[1][0]],
[murangeExp[1][1]] * 2, 'm:')
plt.plot([murangeExp[1][0]] * 2,
[murangeExp[1][1], murangeExp[0][1]], 'm:')
plt.plot([murangeExp[1][0], murangeExp[0][0]],
[murangeExp[0][1]] * 2, 'm:')
plt.colorbar(p)
plt.title("log10|u(mu)|")
plt.grid()
plt.show()
plt.figure(figsize = (15, 7))
plt.jet()
p = plt.contourf(Mu1, Mu2, T, cmap = cmap,
levels = np.linspace(ZTmin, ZTmax, 50))
if clip > 0:
plt.contour(Mu1, Mu2, T, [ZTmin])
plt.plot(np.real(mu2x), np.real(mu2y), 'kx')
+ for j, xy in enumerate(zip(np.real(mu2x), np.real(mu2y))):
+ plt.annotate("{}".format(j), xy)
plt.plot([murangeExp[0][0]] * 2,
[murangeExp[0][1], murangeExp[1][1]], 'm:')
plt.plot([murangeExp[0][0], murangeExp[1][0]],
[murangeExp[1][1]] * 2, 'm:')
plt.plot([murangeExp[1][0]] * 2,
[murangeExp[1][1], murangeExp[0][1]], 'm:')
plt.plot([murangeExp[1][0], murangeExp[0][0]],
[murangeExp[0][1]] * 2, 'm:')
plt.title("log10|u_app(mu)|")
plt.colorbar(p)
plt.grid()
plt.show()
plt.figure(figsize = (15, 7))
plt.jet()
p = plt.contourf(Mu1, Mu2, R, levels = np.linspace(Rmin, Rmax, 50))
plt.plot(np.real(mu2x), np.real(mu2y), 'kx')
+ for j, xy in enumerate(zip(np.real(mu2x), np.real(mu2y))):
+ plt.annotate("{}".format(j), xy)
plt.plot([murangeExp[0][0]] * 2,
[murangeExp[0][1], murangeExp[1][1]], 'm:')
plt.plot([murangeExp[0][0], murangeExp[1][0]],
[murangeExp[1][1]] * 2, 'm:')
plt.plot([murangeExp[1][0]] * 2,
[murangeExp[1][1], murangeExp[0][1]], 'm:')
plt.plot([murangeExp[1][0], murangeExp[0][0]],
[murangeExp[0][1]] * 2, 'm:')
if normalizeDen:
plt.title("log10|Q(mu)res(mu)|")
else:
plt.title("log10|res(mu)|")
plt.colorbar(p)
plt.grid()
plt.show()
plt.figure(figsize = (15, 7))
plt.jet()
p = plt.contourf(Mu1, Mu2, E, levels = np.linspace(Emin, Emax, 50))
plt.plot(np.real(mu2x), np.real(mu2y), 'kx')
+ for j, xy in enumerate(zip(np.real(mu2x), np.real(mu2y))):
+ plt.annotate("{}".format(j), xy)
plt.plot([murangeExp[0][0]] * 2,
[murangeExp[0][1], murangeExp[1][1]], 'm:')
plt.plot([murangeExp[0][0], murangeExp[1][0]],
[murangeExp[1][1]] * 2, 'm:')
plt.plot([murangeExp[1][0]] * 2,
[murangeExp[1][1], murangeExp[0][1]], 'm:')
plt.plot([murangeExp[1][0], murangeExp[0][0]],
[murangeExp[0][1]] * 2, 'm:')
if normalizeDen:
plt.title("log10|Q(mu)err(mu)|")
else:
plt.title("log10|err(mu)|")
plt.colorbar(p)
plt.grid()
plt.show()
if not np.isnan(beta[0, 0]):
plt.figure(figsize = (15, 7))
plt.jet()
p = plt.contourf(Mu1, Mu2, beta,
levels = np.linspace(betamin, betamax, 50))
plt.plot(np.real(mu2x), np.real(mu2y), 'kx')
+ for j, xy in enumerate(zip(np.real(mu2x), np.real(mu2y))):
+ plt.annotate("{}".format(j), xy)
plt.plot([murangeExp[0][0]] * 2,
[murangeExp[0][1], murangeExp[1][1]], 'm:')
plt.plot([murangeExp[0][0], murangeExp[1][0]],
[murangeExp[1][1]] * 2, 'm:')
plt.plot([murangeExp[1][0]] * 2,
[murangeExp[1][1], murangeExp[0][1]], 'm:')
plt.plot([murangeExp[1][0], murangeExp[0][0]],
[murangeExp[0][1]] * 2, 'm:')
plt.title("beta(mu)")
plt.colorbar(p)
plt.grid()
plt.show()
plt.figure(figsize = (15, 7))
plt.jet()
p = plt.contourf(Mu1, Mu2, R - np.log10(beta),
levels = np.linspace(Emin, Emax, 50))
plt.plot(np.real(mu2x), np.real(mu2y), 'kx')
+ for j, xy in enumerate(zip(np.real(mu2x), np.real(mu2y))):
+ plt.annotate("{}".format(j), xy)
plt.plot([murangeExp[0][0]] * 2,
[murangeExp[0][1], murangeExp[1][1]], 'm:')
plt.plot([murangeExp[0][0], murangeExp[1][0]],
[murangeExp[1][1]] * 2, 'm:')
plt.plot([murangeExp[1][0]] * 2,
[murangeExp[1][1], murangeExp[0][1]], 'm:')
plt.plot([murangeExp[1][0], murangeExp[0][0]],
[murangeExp[0][1]] * 2, 'm:')
if normalizeDen:
plt.title("log10|Q(mu)res(mu)/beta(mu)|")
else:
plt.title("log10|res(mu)/beta(mu)|")
plt.colorbar(p)
plt.grid()
plt.show()
plt.figure(figsize = (15, 7))
plt.jet()
p = plt.contourf(Mu1, Mu2, R - np.log10(beta) - E, 50)
plt.plot(np.real(mu2x), np.real(mu2y), 'kx')
+ for j, xy in enumerate(zip(np.real(mu2x), np.real(mu2y))):
+ plt.annotate("{}".format(j), xy)
plt.plot([murangeExp[0][0]] * 2,
[murangeExp[0][1], murangeExp[1][1]], 'm:')
plt.plot([murangeExp[0][0], murangeExp[1][0]],
[murangeExp[1][1]] * 2, 'm:')
plt.plot([murangeExp[1][0]] * 2,
[murangeExp[1][1], murangeExp[0][1]], 'm:')
plt.plot([murangeExp[1][0], murangeExp[0][0]],
[murangeExp[0][1]] * 2, 'm:')
plt.title("log10|eta(mu)|")
plt.colorbar(p)
plt.grid()
plt.show()
def plotInfSet2(murange, murangeEff, approx, mu0, nSamples = 200,
exps = [2., 2.], clip = -1, relative = True,
normalizeDen = False, nobeta = False):
mu1 = np.linspace(murangeEff[0][0] ** exps[0], murangeEff[1][0] ** exps[0],
nSamples)
mu2 = np.linspace(murangeEff[0][1] ** exps[1], murangeEff[1][1] ** exps[1],
nSamples)
mu1s = np.power(mu1, 1. / exps[0])
mu2s = np.power(mu2, 1. / exps[1])
mus = [(m1, m2) for m2 in mu2s for m1 in mu1s]
Ze = approx.normHF(mus).reshape((nSamples, nSamples))
Te = approx.normApprox(mus).reshape((nSamples, nSamples))
Re = approx.normRes(mus).reshape((nSamples, nSamples))
Ee = approx.normErr(mus).reshape((nSamples, nSamples))
if relative:
Fe = approx.normRHS(mus).reshape((nSamples, nSamples))
Re /= Fe
Ee /= Ze
if normalizeDen:
Qvals = np.abs(approx.trainedModel.getQVal(mus).reshape(
(nSamples, nSamples)))
Re *= Qvals
Ee *= Qvals
if nobeta:
betae = np.empty((nSamples, nSamples))
betae[:, :] = np.nan
else:
betae = approx.HFEngine.stabilityFactor(approx.getHF(mus), mus)\
.reshape((nSamples, nSamples))
plotInfSet2FromData(mus, Ze, Te, Re, Ee, betae, murange,
approx, mu0, exps, clip, normalizeDen)
return mus, Ze, Te, Re, Ee, betae
diff --git a/examples/2d/pod/plot_zero_set.py b/examples/2d/pod/plot_zero_set.py
index b452853..8e82cf0 100644
--- a/examples/2d/pod/plot_zero_set.py
+++ b/examples/2d/pod/plot_zero_set.py
@@ -1,95 +1,99 @@
import warnings
import numpy as np
from matplotlib import pyplot as plt
def plotZeroSet1(murange, murangeEff, approx, mu0, nSamples = 200, exp = 2.):
if hasattr(approx, "mus"):
mu2x = approx.mus(0) ** exp
else:
mu2x = mu0[0] ** exp
murangeExp = [[murange[0][0] ** exp], [murange[1][0] ** exp]]
mu1 = np.linspace(murangeEff[0][0] ** exp, murangeEff[1][0] ** exp,
nSamples)
mus = np.power(mu1, 1. / exp)
mu1 = np.real(mu1)
Z = approx.trainedModel.getQVal(mus)
Zabs = np.abs(Z)
Zmin, Zmax = np.min(Zabs), np.max(Zabs)
plt.figure(figsize = (15, 7))
plt.jet()
plt.semilogy(mu1, Zabs)
for l_ in approx.trainedModel.getPoles():
plt.plot([np.real(l_ ** exp)] * 2, [Zmin, Zmax], 'b--')
plt.plot(mu2x, [Zmax] * len(mu2x), 'kx')
+ for j, x in enumerate(mu2x):
+ plt.annotate("{}".format(j), (x, Zmax))
plt.plot([murangeExp[0][0]] * 2, [Zmin, Zmax], 'm:')
plt.plot([murangeExp[1][0]] * 2, [Zmin, Zmax], 'm:')
plt.xlim(mu1[0], mu1[-1])
plt.title("|Q(mu)|")
plt.grid()
plt.show()
return mus, Z
def plotZeroSet2(murange, murangeEff, approx, mu0, nSamples = 200,
exps = [2., 2.], clip = -1, polesImTol : float = 1e-5):
if hasattr(approx, "mus"):
mu2x, mu2y = approx.mus(0) ** exps[0], approx.mus(1) ** exps[1]
else:
mu2x, mu2y = mu0[0] ** exps[0], mu0[1] ** exps[1]
murangeExp = [[murange[0][0] ** exps[0], murange[0][1] ** exps[1]],
[murange[1][0] ** exps[0], murange[1][1] ** exps[1]]]
mu1 = np.linspace(murangeEff[0][0] ** exps[0], murangeEff[1][0] ** exps[0],
nSamples)
mu2 = np.linspace(murangeEff[0][1] ** exps[1], murangeEff[1][1] ** exps[1],
nSamples)
mu1s = np.power(mu1, 1. / exps[0])
mu2s = np.power(mu2, 1. / exps[1])
Mu1, Mu2 = np.meshgrid(np.real(mu1), np.real(mu2))
mus = [(m1, m2) for m2 in mu2s for m1 in mu1s]
poles1, poles2 = [], []
for m2 in mu2s:
pls = approx.getPoles([None, m2]) ** exps[0]
pls[np.imag(pls) > polesImTol] = np.nan
pls = np.real(pls)
poles1 += list(pls)
poles2 += [m2] * len(pls)
Z = approx.trainedModel.getQVal(mus).reshape(Mu1.shape)
Zabs = np.log10(np.abs(Z))
Zabsmin, Zabsmax = np.min(Zabs), np.max(Zabs)
if clip > 0:
Zabsmin += clip * (Zabsmax - Zabsmin)
cmap = plt.cm.bone_r
else:
cmap = plt.cm.jet
try:
pass
except:
Z = None
cmap = plt.cm.jet
warnings.simplefilter("ignore", category = (UserWarning,
np.ComplexWarning))
plt.figure(figsize = (15, 7))
plt.jet()
if Z is not None:
p = plt.contourf(Mu1, Mu2, Zabs, cmap = cmap,
levels = np.linspace(Zabsmin, Zabsmax, 50))
if clip > 0:
plt.contour(Mu1, Mu2, Zabs, [Zabsmin])
plt.plot(poles1, poles2, 'k.')
plt.plot(np.real(mu2x), np.real(mu2y), 'kx')
+ for j, xy in enumerate(zip(np.real(mu2x), np.real(mu2y))):
+ plt.annotate("{}".format(j), xy)
plt.plot([murangeExp[0][0]] * 2,
[murangeExp[0][1], murangeExp[1][1]], 'm:')
plt.plot([murangeExp[0][0], murangeExp[1][0]],
[murangeExp[1][1]] * 2, 'm:')
plt.plot([murangeExp[1][0]] * 2,
[murangeExp[1][1], murangeExp[0][1]], 'm:')
plt.plot([murangeExp[1][0], murangeExp[0][0]],
[murangeExp[0][1]] * 2, 'm:')
if Z is not None:
plt.colorbar(p)
plt.xlim(murangeExp[0][0], murangeExp[1][0])
plt.ylim(murangeExp[0][1], murangeExp[1][1])
plt.title("log10|Q(mu)|")
plt.grid()
plt.show()
return mus, Z
diff --git a/examples/from_papers/greedy_internalBox.py b/examples/from_papers/greedy_internalBox.py
index 2e92707..336f8d5 100644
--- a/examples/from_papers/greedy_internalBox.py
+++ b/examples/from_papers/greedy_internalBox.py
@@ -1,116 +1,117 @@
import numpy as np
import fenics as fen
from rrompy.hfengines.linear_problem import HelmholtzProblemEngine as HPE
from rrompy.reduction_methods.greedy import RationalInterpolantGreedy as RI
from rrompy.reduction_methods.greedy import ReducedBasisGreedy as RB
from rrompy.solver.fenics import L2NormMatrix
from rrompy.parameter.parameter_sampling import QuadratureSampler as QS
-dim = 3
+dim = 2
verb = 5
timed = False
algo = "RI"
#algo = "RB"
polyBasis = "LEGENDRE"
-#polyBasis = "CHEBYSHEV"
+polyBasis = "CHEBYSHEV"
#polyBasis = "MONOMIAL"
k0s = np.power(np.linspace(500 ** 2., 2250 ** 2., 200), .5)
k0 = np.mean(np.power(k0s, 2.)) ** .5
kl, kr = min(k0s), max(k0s)
params = {'sampler':QS([kl, kr], "UNIFORM", 2.), 'nTestPoints':500,
- 'greedyTol':1e-2, 'S':2, 'polybasis':polyBasis,
+ 'greedyTol':1e2, 'S':2, 'polybasis':polyBasis,
# 'errorEstimatorKind':'AFFINE'}
+ 'errorEstimatorKind':'DISCREPANCY'}
# 'errorEstimatorKind':'INTERPOLATORY'}
- 'errorEstimatorKind':'LOCALL2'}
+# 'errorEstimatorKind':'LOCALL2'}
if dim == 2:
mesh = fen.RectangleMesh(fen.Point(0., 0.), fen.Point(.1, .15), 10, 15)
x, y = fen.SpatialCoordinate(mesh)[:]
f = fen.exp(- 1e2 * (x + y))
else:#if dim == 3:
mesh = fen.BoxMesh(fen.Point(0., 0., 0.), fen.Point(.1, .15, .25), 4, 6,10)
x, y, z = fen.SpatialCoordinate(mesh)[:]
f = fen.exp(- 1e2 * (x + y + z))
solver = HPE(np.real(k0), verbosity = verb)
solver.V = fen.FunctionSpace(mesh, "P", 3)
solver.refractionIndex = fen.Constant(1. / 54.6)
solver.forcingTerm = f
solver.NeumannBoundary = "ALL"
#########
if algo == "RI":
approx = RI(solver, mu0 = k0, approxParameters = params, verbosity = verb)
else:
params.pop("polybasis")
params.pop("errorEstimatorKind")
approx = RB(solver, mu0 = k0, approxParameters = params, verbosity = verb)
approx.initEstimatorNormEngine(L2NormMatrix(solver.V))
if timed:
from time import clock
start_time = clock()
approx.greedy()
print("Time: ", clock() - start_time)
else:
approx.greedy(True)
approx.samplingEngine.verbosity = 0
approx.verbosity = 0
kl, kr = np.real(kl), np.real(kr)
from matplotlib import pyplot as plt
normApp = np.zeros(len(k0s))
norm = np.zeros_like(normApp)
res = np.zeros_like(normApp)
err = np.zeros_like(normApp)
for j in range(len(k0s)):
normApp[j] = approx.normApprox(k0s[j])
norm[j] = approx.normHF(k0s[j])
res[j] = (approx.estimatorNormEngine.norm(
approx.getRes(k0s[j], duality = False))
/ approx.estimatorNormEngine.norm(
approx.getRHS(k0s[j], duality = False)))
err[j] = approx.normErr(k0s[j]) / approx.normHF(k0s[j])
resApp = approx.errorEstimator(k0s)
plt.figure()
plt.plot(k0s, norm)
plt.plot(k0s, normApp, '--')
plt.plot(np.real(approx.mus(0)),
1.05*np.max(norm)*np.ones(approx.mus.shape, dtype = float), 'rx')
plt.xlim([kl, kr])
plt.grid()
plt.show()
plt.close()
plt.figure()
plt.semilogy(k0s, res)
plt.semilogy(k0s, resApp, '--')
plt.semilogy(np.real(approx.mus(0)),
4.*np.max(resApp)*np.ones(approx.mus.shape, dtype = float), 'rx')
plt.xlim([kl, kr])
plt.grid()
plt.show()
plt.close()
plt.figure()
plt.semilogy(k0s, err)
plt.xlim([kl, kr])
plt.grid()
plt.show()
plt.close()
polesApp = approx.getPoles()
mask = (np.real(polesApp) < kl) | (np.real(polesApp) > kr)
print("Outliers:", polesApp[mask])
polesApp = polesApp[~mask]
plt.figure()
plt.plot(np.real(polesApp), np.imag(polesApp), 'kx')
plt.axis('equal')
plt.grid()
plt.show()
plt.close()
diff --git a/rrompy/reduction_methods/greedy/generic_greedy_approximant.py b/rrompy/reduction_methods/greedy/generic_greedy_approximant.py
index bb81d6f..5e76d32 100644
--- a/rrompy/reduction_methods/greedy/generic_greedy_approximant.py
+++ b/rrompy/reduction_methods/greedy/generic_greedy_approximant.py
@@ -1,701 +1,702 @@
# 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.standard.generic_standard_approximant \
import GenericStandardApproximant
from rrompy.utilities.base.types import (Np1D, Np2D, DictAny, HFEng, Tuple,
List, normEng, paramVal, paramList,
sampList)
from rrompy.utilities.poly_fitting.polynomial import (polyvanderTotal as pvT,
hashIdxToDerivative as hashI,
totalDegreeN)
from rrompy.utilities.base import verbosityManager as vbMng
from rrompy.solver import normEngine
from rrompy.utilities.exception_manager import (RROMPyException, RROMPyAssert,
RROMPyWarning)
from rrompy.parameter import checkParameterList, emptyParameterList
__all__ = ['GenericGreedyApproximant']
def localL2Distance(mus:Np2D, badmus:Np2D) -> Np2D:
return np.linalg.norm(np.tile(mus[..., np.newaxis], [1, 1, len(badmus)])
- badmus[..., np.newaxis].T, axis = 1)
def pruneSamples(mus:paramList, badmus:paramList,
tol : float = 1e-8) -> Np1D:
"""Remove from mus all the elements which are too close to badmus."""
if len(badmus) == 0: return mus
musA = mus.data
badmusA = badmus.data
proximity = np.min(localL2Distance(musA, badmusA), axis = 1)
return np.arange(len(mus))[proximity <= tol]
class GenericGreedyApproximant(GenericStandardApproximant):
"""
ROM greedy interpolant computation for parametric problems
(ABSTRACT).
Args:
HFEngine: HF problem solver.
mu0(optional): Default parameter. Defaults to 0.
approxParameters(optional): Dictionary containing values for main
parameters of approximant. Recognized keys are:
- 'POD': whether to compute POD of snapshots; defaults to True;
- 'S': number of starting training points;
- 'sampler': sample point generator;
- 'greedyTol': uniform error tolerance for greedy algorithm;
defaults to 1e-2;
- 'collinearityTol': collinearity tolerance for greedy algorithm;
defaults to 0.;
- 'interactive': whether to interactively terminate greedy
algorithm; defaults to False;
- 'maxIter': maximum number of greedy steps; defaults to 1e2;
- 'refinementRatio': ratio of test points to be exhausted before
test set refinement; defaults to 0.2;
- 'nTestPoints': number of test points; defaults to 5e2;
- 'trainSetGenerator': training sample points generator; defaults
to sampler.
Defaults to empty dict.
verbosity(optional): Verbosity level. Defaults to 10.
Attributes:
HFEngine: HF problem solver.
mu0: Default parameter.
mus: Array of snapshot parameters.
approxParameters: Dictionary containing values for main parameters of
approximant. Recognized keys are in parameterList.
parameterListSoft: Recognized keys of soft approximant parameters:
- 'POD': whether to compute POD of snapshots.
- 'greedyTol': uniform error tolerance for greedy algorithm;
- 'collinearityTol': collinearity tolerance for greedy algorithm;
- 'interactive': whether to interactively terminate greedy
algorithm;
- 'maxIter': maximum number of greedy steps;
- 'refinementRatio': ratio of test points to be exhausted before
test set refinement;
- 'nTestPoints': number of test points;
- 'trainSetGenerator': training sample points generator.
parameterListCritical: Recognized keys of critical approximant
parameters:
- 'S': total number of samples current approximant relies upon;
- 'sampler': sample point generator.
POD: whether to compute POD of snapshots.
S: number of test points.
sampler: Sample point generator.
greedyTol: Uniform error tolerance for greedy algorithm.
collinearityTol: Collinearity tolerance for greedy algorithm.
interactive: whether to interactively terminate greedy algorithm.
maxIter: maximum number of greedy steps.
refinementRatio: ratio of training points to be exhausted before
training set refinement.
nTestPoints: number of starting training points.
trainSetGenerator: training sample points generator.
muBounds: list of bounds for parameter values.
samplingEngine: Sampling engine.
estimatorNormEngine: Engine for estimator norm computation.
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.
"""
TOL_INSTABILITY = 1e-6
def __init__(self, HFEngine:HFEng, mu0 : paramVal = None,
approxParameters : DictAny = {}, verbosity : int = 10,
timestamp : bool = True):
self._preInit()
self._addParametersToList(["greedyTol", "collinearityTol",
"interactive", "maxIter", "refinementRatio",
"nTestPoints"],
[1e-2, 0., False, 1e2, .2, 5e2],
["trainSetGenerator"], ["AUTO"])
super().__init__(HFEngine = HFEngine, mu0 = mu0,
approxParameters = approxParameters,
verbosity = verbosity, timestamp = timestamp)
self._postInit()
@property
def greedyTol(self):
"""Value of greedyTol."""
return self._greedyTol
@greedyTol.setter
def greedyTol(self, greedyTol):
if greedyTol < 0:
raise RROMPyException("greedyTol must be non-negative.")
if hasattr(self, "_greedyTol") and self.greedyTol is not None:
greedyTolold = self.greedyTol
else:
greedyTolold = -1
self._greedyTol = greedyTol
self._approxParameters["greedyTol"] = self.greedyTol
if greedyTolold != self.greedyTol:
self.resetSamples()
@property
def collinearityTol(self):
"""Value of collinearityTol."""
return self._collinearityTol
@collinearityTol.setter
def collinearityTol(self, collinearityTol):
if collinearityTol < 0:
raise RROMPyException("collinearityTol must be non-negative.")
if (hasattr(self, "_collinearityTol")
and self.collinearityTol is not None):
collinearityTolold = self.collinearityTol
else:
collinearityTolold = -1
self._collinearityTol = collinearityTol
self._approxParameters["collinearityTol"] = self.collinearityTol
if collinearityTolold != self.collinearityTol:
self.resetSamples()
@property
def interactive(self):
"""Value of interactive."""
return self._interactive
@interactive.setter
def interactive(self, interactive):
self._interactive = interactive
@property
def maxIter(self):
"""Value of maxIter."""
return self._maxIter
@maxIter.setter
def maxIter(self, maxIter):
if maxIter <= 0: raise RROMPyException("maxIter must be positive.")
if hasattr(self, "_maxIter") and self.maxIter is not None:
maxIterold = self.maxIter
else:
maxIterold = -1
self._maxIter = maxIter
self._approxParameters["maxIter"] = self.maxIter
if maxIterold != self.maxIter:
self.resetSamples()
@property
def refinementRatio(self):
"""Value of refinementRatio."""
return self._refinementRatio
@refinementRatio.setter
def refinementRatio(self, refinementRatio):
if refinementRatio <= 0. or refinementRatio > 1.:
raise RROMPyException(("refinementRatio must be between 0 "
"(excluded) and 1."))
if (hasattr(self, "_refinementRatio")
and self.refinementRatio is not None):
refinementRatioold = self.refinementRatio
else:
refinementRatioold = -1
self._refinementRatio = refinementRatio
self._approxParameters["refinementRatio"] = self.refinementRatio
if refinementRatioold != self.refinementRatio:
self.resetSamples()
@property
def nTestPoints(self):
"""Value of nTestPoints."""
return self._nTestPoints
@nTestPoints.setter
def nTestPoints(self, nTestPoints):
if nTestPoints <= 0:
raise RROMPyException("nTestPoints must be positive.")
if not np.isclose(nTestPoints, np.int(nTestPoints)):
raise RROMPyException("nTestPoints must be an integer.")
nTestPoints = np.int(nTestPoints)
if hasattr(self, "_nTestPoints") and self.nTestPoints is not None:
nTestPointsold = self.nTestPoints
else:
nTestPointsold = -1
self._nTestPoints = nTestPoints
self._approxParameters["nTestPoints"] = self.nTestPoints
if nTestPointsold != self.nTestPoints:
self.resetSamples()
@property
def trainSetGenerator(self):
"""Value of trainSetGenerator."""
return self._trainSetGenerator
@trainSetGenerator.setter
def trainSetGenerator(self, trainSetGenerator):
if (isinstance(trainSetGenerator, (str,))
and trainSetGenerator.upper() == "AUTO"):
trainSetGenerator = self.sampler
if 'generatePoints' not in dir(trainSetGenerator):
raise RROMPyException("trainSetGenerator type not recognized.")
if (hasattr(self, '_trainSetGenerator')
and self.trainSetGenerator not in [None, "AUTO"]):
trainSetGeneratorOld = self.trainSetGenerator
self._trainSetGenerator = trainSetGenerator
self._approxParameters["trainSetGenerator"] = self.trainSetGenerator
if (not 'trainSetGeneratorOld' in locals()
or trainSetGeneratorOld != self.trainSetGenerator):
self.resetSamples()
def resetSamples(self):
"""Reset samples."""
super().resetSamples()
self._mus = emptyParameterList()
def initEstimatorNormEngine(self, normEngn : normEng = None):
"""Initialize estimator norm engine."""
if (normEngn is not None or not hasattr(self, "estimatorNormEngine")
or self.estimatorNormEngine is None):
if normEngn is None:
if not hasattr(self.HFEngine, "energyNormPartialDualMatrix"):
self.HFEngine.buildEnergyNormPartialDualForm()
estimatorEnergyMatrix = (
self.HFEngine.energyNormPartialDualMatrix)
else:
if hasattr(normEngn, "buildEnergyNormPartialDualForm"):
if not hasattr(normEngn, "energyNormPartialDualMatrix"):
normEngn.buildEnergyNormPartialDualForm()
estimatorEnergyMatrix = (
normEngn.energyNormPartialDualMatrix)
else:
estimatorEnergyMatrix = normEngn
self.estimatorNormEngine = normEngine(estimatorEnergyMatrix)
def errorEstimator(self, mus:Np1D) -> Np1D:
"""Standard residual-based error estimator."""
self.setupApprox()
- self.assembleReducedResidualBlocks(full = True)
- nmus = len(mus)
- nAs = self.trainedModel.data.resAA.shape[1]
- nbs = self.trainedModel.data.resbb.shape[0]
mustr = mus
+ nmus = len(mus)
if nmus > 2:
mustr = "[{} ..({}).. {}]".format(mus[0], nmus - 2, mus[-1])
vbMng(self.trainedModel, "INIT",
- "Computing reduced solution at mu = {}.".format(mustr), 5)
- verb = self.trainedModel.verbosity
- self.trainedModel.verbosity = 0
+ "Evaluating error estimator at mu = {}.".format(mustr), 10)
+ self.assembleReducedResidualBlocks(full = True)
+ nAs = self.trainedModel.data.resAA.shape[1]
+ nbs = self.trainedModel.data.resbb.shape[0]
mus = checkParameterList(mus, self.npar)[0]
uApproxRs = self.getApproxReduced(mus)
- muRTest = ((mus ** self.HFEngine.rescalingExp
- - self.mu0 ** self.HFEngine.rescalingExp)
- / self.scaleFactor)
+ muRTest = (mus ** self.HFEngine.rescalingExp
+ - self.mu0 ** self.HFEngine.rescalingExp) / self.scaleFactor
vanderBase, _, vanderBaseIdxs = pvT(muRTest,
np.sum(hashI(max(nAs, nbs) - 1, self.npar)),
"MONOMIAL")
vanderBase = vanderBase[:, vanderBaseIdxs].T
radiusA = np.expand_dims(uApproxRs.data, 1) * vanderBase[: nAs, :]
radiusb = vanderBase[: nbs, :]
- self.trainedModel.verbosity = verb
- vbMng(self.trainedModel, "DEL", "Done computing reduced solution.", 5)
# 'ij,jk,ik->k', resbb, radiusb, radiusb.conj()
ff = np.sum(self.trainedModel.data.resbb.dot(radiusb) * radiusb.conj(),
axis = 0)
# 'ijk,jkl,il->l', resAb, radiusA, radiusb.conj()
Lf = np.sum(np.tensordot(self.trainedModel.data.resAb, radiusA, 2)
* radiusb.conj(), axis = 0)
# 'ijkl,klt,ijt->t', resAA, radiusA, radiusA.conj()
LL = np.sum(np.tensordot(self.trainedModel.data.resAA, radiusA, 2)
* radiusA.conj(), axis = (0, 1))
- return np.abs((LL - 2. * np.real(Lf) + ff) / ff) ** .5
+ err = np.abs((LL - 2. * np.real(Lf) + ff) / ff) ** .5
+ vbMng(self.trainedModel, "DEL", "Done evaluating error estimator", 10)
+ return err
def getMaxErrorEstimator(self, mus:paramList, n : int = 1,
plot : bool = False) -> Tuple[Np1D, int, float]:
"""
Compute maximum of (and index of maximum of) error estimator over given
parameters.
"""
errorEstTest = self.errorEstimator(mus)
errCP = copy(errorEstTest)
idxMaxEst = np.empty(n, dtype = int)
for j in range(n):
k = np.argmax(errCP)
idxMaxEst[j] = k
if j + 1 < n:
errCP *= np.linalg.norm((mus.data ** self.HFEngine.rescalingExp
- mus[k] ** self.HFEngine.rescalingExp)
/ self.scaleFactor, axis = 1)
maxEst = errorEstTest[idxMaxEst]
if plot and not np.any(np.isinf(errorEstTest)):
musre = copy(mus.re.data)
from matplotlib import pyplot as plt
plt.figure()
- plt.semilogy([musre[0, 0], musre[-1, 0]], [self.greedyTol] * 2,
- 'r--')
- plt.semilogy(self.mus.re(0),
- 2. * self.greedyTol * np.ones(len(self.mus)), '*m')
- plt.semilogy(musre[idxMaxEst, 0], maxEst, 'xr')
+ errCP = copy(errorEstTest)
while len(musre) > 0:
if self.npar == 1:
currIdx = np.arange(len(musre))
else:
currIdx = np.where(np.isclose(np.sum(
np.abs(musre[:, 1 :] - musre[0, 1 :]), 1), 0.))[0]
- plt.semilogy(musre[currIdx, 0], errorEstTest[currIdx], 'k')
+ plt.semilogy(musre[currIdx, 0], errCP[currIdx], 'k',
+ linewidth = 1)
musre = np.delete(musre, currIdx, 0)
+ errCP = np.delete(errCP, currIdx)
+ plt.semilogy([mus.re(0, 0), mus.re(-1, 0)], [self.greedyTol] * 2,
+ 'r--')
+ plt.semilogy(self.mus.re(0),
+ 2. * self.greedyTol * np.ones(len(self.mus)), '*m')
+ plt.semilogy(mus.re(idxMaxEst, 0), maxEst, 'xr')
plt.grid()
plt.show()
plt.close()
return errorEstTest, idxMaxEst, maxEst
def _isLastSampleCollinear(self) -> bool:
"""Check collinearity of last sample."""
if self.collinearityTol <= 0.: return False
if self.POD:
reff = self.samplingEngine.RPOD[:, -1]
else:
RROMPyWarning(("Repeated orthogonalization of the samples for "
"collinearity check. Consider setting POD to "
"True."))
if not hasattr(self, "_PODEngine"):
from rrompy.sampling.base.pod_engine import PODEngine
self._PODEngine = PODEngine(self.HFEngine)
reff = self._PODEngine.generalizedQR(self.samplingEngine.samples,
only_R = True)[:, -1]
return np.abs(reff[-1]) < self.collinearityTol * np.linalg.norm(reff)
- def greedyNextSample(self, muidx:int, n : int = 1, plotEst : bool = False)\
+ def greedyNextSample(self, muidx:int, plotEst : bool = False)\
-> Tuple[Np1D, int, float, paramVal]:
"""Compute next greedy snapshot of solution map."""
RROMPyAssert(self._mode, message = "Cannot add greedy sample.")
mus = copy(self.muTest[muidx])
self.muTest.pop(muidx)
for mu in mus:
vbMng(self, "MAIN",
("Adding {}-th sample point at {} to training "
"set.").format(self.samplingEngine.nsamples + 1, mu), 2)
self.mus.append(mu)
self.samplingEngine.nextSample(mu)
if self._isLastSampleCollinear():
RROMPyWarning("Collinearity above tolerance detected.")
errorEstTest = np.empty(len(self.muTest))
errorEstTest[:] = np.nan
return errorEstTest, [-1], np.nan, np.nan
- errorEstTest, muidx, maxErrorEst = self.getMaxErrorEstimator(
- self.muTest, n, plotEst)
- self.derivativeShell += 1
+ self.derivativeShellIdx += 1
self.derivativeShellSize = totalDegreeN(self.npar - 1,
- self.derivativeShell)
+ self.derivativeShellIdx)
+ errorEstTest, muidx, maxErrorEst = self.getMaxErrorEstimator(
+ self.muTest, self.derivativeShellSize, plotEst)
return errorEstTest, muidx, maxErrorEst, self.muTest[muidx]
def _preliminaryTraining(self):
"""Initialize starting snapshots of solution map."""
RROMPyAssert(self._mode, message = "Cannot start greedy algorithm.")
self.computeScaleFactor()
if self.samplingEngine.nsamples > 0:
return
self.resetSamples()
self.mus = self.trainSetGenerator.generatePoints(self.S)
- self.derivativeShell, self.derivativeShellSize, shellTot = -1, 0, 0
- while shellTot <= len(self.mus):
- self.derivativeShellSize = totalDegreeN(self.npar - 1,
- self.derivativeShell + 1)
- shellTot += self.derivativeShellSize
- self.derivativeShell += 1
- self._S = shellTot - self.derivativeShellSize
+ self.derivativeShellIdx, self.derivativeShellSize, self._S = -1, 0, 0
+ nextShellSize = 1
+ while self._S + nextShellSize <= len(self.mus):
+ self.derivativeShellIdx += 1
+ self.derivativeShellSize = nextShellSize
+ self._S += self.derivativeShellSize
+ nextShellSize = totalDegreeN(self.npar - 1,
+ self.derivativeShellIdx + 1)
self.mus = self.mus[list(range(self.S))]
muTestBase = self.sampler.generatePoints(self.nTestPoints)
idxPop = pruneSamples(muTestBase ** self.HFEngine.rescalingExp,
self.mus ** self.HFEngine.rescalingExp,
1e-10 * self.scaleFactor[0])
muTestBase.pop(idxPop)
muTestBase = muTestBase.sort()
muLast = copy(self.mus[-1])
self.mus.pop()
if len(self.mus) > 0:
vbMng(self, "MAIN",
("Adding first {} sample point{} at {} to training "
"set.").format(self.S - 1, "" + "s" * (self.S > 2),
self.mus), 2)
self.samplingEngine.iterSample(self.mus)
self.muTest = emptyParameterList()
self.muTest.reset((len(muTestBase) + 1, self.mus.shape[1]))
- self.muTest[: -1] = muTestBase
- self.muTest[-1] = muLast
+ self.muTest[: -1] = muTestBase.data
+ self.muTest[-1] = muLast.data
def _enrichTestSet(self, nTest:int):
"""Add extra elements to test set."""
RROMPyAssert(self._mode, message = "Cannot enrich test set.")
muTestExtra = self.sampler.generatePoints(2 * nTest)
muTotal = copy(self.mus)
muTotal.append(self.muTest)
idxPop = pruneSamples(muTestExtra ** self.HFEngine.rescalingExp,
muTotal ** self.HFEngine.rescalingExp,
1e-10 * self.scaleFactor[0])
muTestExtra.pop(idxPop)
- muTestNew = np.empty(len(self.muTest) + len(muTestExtra),
- dtype = np.complex)
- muTestNew[: len(self.muTest)] = self.muTest(0)
- muTestNew[len(self.muTest) :] = muTestExtra(0)
- self.muTest = checkParameterList(np.sort(muTestNew), self.npar)[0]
+ muTestNew = np.empty((len(self.muTest) + len(muTestExtra),
+ self.muTest.shape[1]), dtype = np.complex)
+ muTestNew[: len(self.muTest)] = self.muTest.data
+ muTestNew[len(self.muTest) :] = muTestExtra.data
+ self.muTest = checkParameterList(muTestNew, self.npar)[0].sort()
vbMng(self, "MAIN",
"Enriching test set by {} elements.".format(len(muTestExtra)), 5)
def greedy(self, plotEst : bool = False):
"""Compute greedy snapshots of solution map."""
RROMPyAssert(self._mode, message = "Cannot start greedy algorithm.")
if self.samplingEngine.nsamples > 0:
return
vbMng(self, "INIT", "Starting computation of snapshots.", 2)
self._preliminaryTraining()
nTest = self.nTestPoints
muT0 = copy(self.muTest[-1])
errorEstTest, muidx, maxErrorEst, mu = self.greedyNextSample(
- [len(self.muTest) - 1], self.derivativeShellSize, plotEst)
+ [len(self.muTest) - 1], plotEst)
if np.any(np.isnan(maxErrorEst)):
RROMPyWarning(("Instability in a posteriori estimator. "
"Starting preemptive greedy loop termination."))
self.muTest.append(muT0)
self.mus.pop(-1)
self.samplingEngine.popSample()
self.setupApprox()
else:
vbMng(self, "MAIN", ("Uniform testing error estimate "
"{:.4e}.").format(np.max(maxErrorEst)), 2)
trainedModelOld = copy(self.trainedModel)
while (self.samplingEngine.nsamples < self.maxIter
and np.max(maxErrorEst) > self.greedyTol):
if (1. - self.refinementRatio) * nTest > len(self.muTest):
self._enrichTestSet(nTest)
nTest = len(self.muTest)
muTestOld, maxErrorEstOld = self.muTest, np.max(maxErrorEst)
errorEstTest, muidx, maxErrorEst, mu = self.greedyNextSample(
- muidx, self.derivativeShellSize, plotEst)
+ muidx, plotEst)
vbMng(self, "MAIN", ("Uniform testing error estimate "
"{:.4e}.").format(np.max(maxErrorEst)), 2)
if (np.any(np.isnan(maxErrorEst))
or np.any(np.isinf(maxErrorEst))
or maxErrorEstOld < (np.max(maxErrorEst)
* self.TOL_INSTABILITY)):
RROMPyWarning(("Instability in a posteriori estimator. "
"Starting preemptive greedy loop "
"termination."))
self.muTest = muTestOld
self.mus.pop(-1)
self.samplingEngine.popSample()
self.trainedModel.data = copy(trainedModelOld.data)
break
trainedModelOld.data = copy(self.trainedModel.data)
if (self.interactive and np.max(maxErrorEst) <= self.greedyTol):
vbMng(self, "MAIN",
("Required tolerance {} achieved. Want to decrease "
"greedyTol and continue? "
"Y/N").format(self.greedyTol), 0, end = "")
increasemaxIter = input()
if increasemaxIter.upper() == "Y":
vbMng(self, "MAIN", "Reducing value of greedyTol...",
0)
while np.max(maxErrorEst) <= self._greedyTol:
self._greedyTol *= .5
if (self.interactive
and self.samplingEngine.nsamples >= self.maxIter):
vbMng(self, "MAIN",
("Maximum number of iterations {} reached. Want to "
"increase maxIter and continue? "
"Y/N").format(self.maxIter), 0, end = "")
increasemaxIter = input()
if increasemaxIter.upper() == "Y":
vbMng(self, "MAIN", "Doubling value of maxIter...", 0)
self._maxIter *= 2
vbMng(self, "DEL",
("Done computing snapshots (final snapshot count: "
"{}).").format(self.samplingEngine.nsamples), 2)
def checkComputedApprox(self) -> bool:
"""
Check if setup of new approximant is not needed.
Returns:
True if new setup is not needed. False otherwise.
"""
return (super().checkComputedApprox()
and len(self.mus) == self.trainedModel.data.projMat.shape[1])
def assembleReducedResidualGramian(self, pMat:sampList):
"""
Build residual gramian of reduced linear system through projections.
"""
self.initEstimatorNormEngine()
if (not hasattr(self.trainedModel.data, "gramian")
or self.trainedModel.data.gramian is None):
gramian = self.estimatorNormEngine.innerProduct(pMat, pMat)
else:
Sold = self.trainedModel.data.gramian.shape[0]
S = len(self.mus)
if Sold > S:
gramian = self.trainedModel.data.gramian[: S, : S]
else:
idxOld = list(range(Sold))
idxNew = list(range(Sold, S))
gramian = np.empty((S, S), dtype = np.complex)
gramian[: Sold, : Sold] = self.trainedModel.data.gramian
gramian[: Sold, Sold :] = (
self.estimatorNormEngine.innerProduct(pMat(idxNew),
pMat(idxOld)))
gramian[Sold :, : Sold] = gramian[: Sold, Sold :].T.conj()
gramian[Sold :, Sold :] = (
self.estimatorNormEngine.innerProduct(pMat(idxNew),
pMat(idxNew)))
self.trainedModel.data.gramian = gramian
def assembleReducedResidualBlocksbb(self, bs:List[Np1D]):
"""
Build blocks (of type bb) of reduced linear system through projections.
"""
self.initEstimatorNormEngine()
nbs = len(bs)
if (not hasattr(self.trainedModel.data, "resbb")
or self.trainedModel.data.resbb is None):
resbb = np.empty((nbs, nbs), dtype = np.complex)
for i in range(nbs):
Mbi = bs[i]
resbb[i, i] = self.estimatorNormEngine.innerProduct(Mbi, Mbi)
for j in range(i):
Mbj = bs[j]
resbb[i, j] = self.estimatorNormEngine.innerProduct(Mbj,
Mbi)
for i in range(nbs):
for j in range(i + 1, nbs):
resbb[i, j] = resbb[j, i].conj()
self.trainedModel.data.resbb = resbb
def assembleReducedResidualBlocksAb(self, As:List[Np2D], bs:List[Np1D],
pMat:sampList):
"""
Build blocks (of type Ab) of reduced linear system through projections.
"""
self.initEstimatorNormEngine()
nAs = len(As)
nbs = len(bs)
S = len(self.mus)
if (not hasattr(self.trainedModel.data, "resAb")
or self.trainedModel.data.resAb is None):
if not isinstance(pMat, (np.ndarray,)): pMat = pMat.data
resAb = np.empty((nbs, S, nAs), dtype = np.complex)
for j in range(nAs):
MAj = As[j].dot(pMat)
for i in range(nbs):
Mbi = bs[i]
resAb[i, :, j] = self.estimatorNormEngine.innerProduct(MAj,
Mbi)
else:
Sold = self.trainedModel.data.resAb.shape[1]
if Sold == S: return
if Sold > S:
resAb = self.trainedModel.data.resAb[:, : S, :]
else:
if not isinstance(pMat, (np.ndarray,)): pMat = pMat.data
resAb = np.empty((nbs, S, nAs), dtype = np.complex)
resAb[:, : Sold, :] = self.trainedModel.data.resAb
for j in range(nAs):
MAj = As[j].dot(pMat[:, Sold :])
for i in range(nbs):
Mbi = bs[i]
resAb[i, Sold :, j] = (
self.estimatorNormEngine.innerProduct(MAj, Mbi))
self.trainedModel.data.resAb = resAb
def assembleReducedResidualBlocksAA(self, As:List[Np2D], pMat:sampList):
"""
Build blocks (of type AA) of reduced linear system through projections.
"""
self.initEstimatorNormEngine()
nAs = len(As)
S = len(self.mus)
if (not hasattr(self.trainedModel.data, "resAA")
or self.trainedModel.data.resAA is None):
if not isinstance(pMat, (np.ndarray,)): pMat = pMat.data
resAA = np.empty((S, nAs, S, nAs), dtype = np.complex)
for i in range(nAs):
MAi = As[i].dot(pMat)
resAA[:, i, :, i] = (
self.estimatorNormEngine.innerProduct(MAi, MAi))
for j in range(i):
MAj = As[j].dot(pMat)
resAA[:, i, :, j] = (
self.estimatorNormEngine.innerProduct(MAj, MAi))
for i in range(nAs):
for j in range(i + 1, nAs):
resAA[:, i, :, j] = resAA[:, j, :, i].T.conj()
else:
Sold = self.trainedModel.data.resAA.shape[0]
if Sold == S: return
if Sold > S:
resAA = self.trainedModel.data.resAA[: S, :, : S, :]
else:
if not isinstance(pMat, (np.ndarray,)): pMat = pMat.data
resAA = np.empty((S, nAs, S, nAs), dtype = np.complex)
resAA[: Sold, :, : Sold, :] = self.trainedModel.data.resAA
for i in range(nAs):
MAi = As[i].dot(pMat)
resAA[: Sold, i, Sold :, i] = (
self.estimatorNormEngine.innerProduct(MAi[:, Sold :],
MAi[:, : Sold]))
resAA[Sold :, i, : Sold, i] = resAA[: Sold, i,
Sold :, i].T.conj()
resAA[Sold :, i, Sold :, i] = (
self.estimatorNormEngine.innerProduct(MAi[:, Sold :],
MAi[:, Sold :]))
for j in range(i):
MAj = As[j].dot(pMat)
resAA[: Sold, i, Sold :, j] = (
self.estimatorNormEngine.innerProduct(MAj[:, Sold :],
MAi[:, : Sold]))
resAA[Sold :, i, : Sold, j] = (
self.estimatorNormEngine.innerProduct(MAj[:, : Sold],
MAi[:, Sold :]))
resAA[Sold :, i, Sold :, j] = (
self.estimatorNormEngine.innerProduct(MAj[:, Sold :],
MAi[:, Sold :]))
for i in range(nAs):
for j in range(i + 1, nAs):
resAA[: Sold, i, Sold :, j] = (
resAA[Sold :, j, : Sold, i].T.conj())
resAA[Sold :, i, : Sold, j] = (
resAA[: Sold, j, Sold :, i].T.conj())
resAA[Sold :, i, Sold :, j] = (
resAA[Sold :, j, Sold :, i].T.conj())
self.trainedModel.data.resAA = resAA
def _setupAffineBlocks(self, full : bool = True):
bmsg = "" if full else " of RHS"
if (full and not hasattr(self, "As")) or not hasattr(self, "bs"):
vbMng(self, "INIT",
"Computing affine blocks{} of system.".format(bmsg), 10)
if full:
self.As = self.HFEngine.affineLinearSystemA(self.mu0,
self.scaleFactor)
self.bs = self.HFEngine.affineLinearSystemb(self.mu0,
self.scaleFactor)
vbMng(self, "DEL",
"Done computing affine blocks{}.".format(bmsg), 10)
def assembleReducedResidualBlocks(self, full : bool = False):
"""Build affine blocks of affine decomposition of residual."""
self._setupAffineBlocks(full)
self.assembleReducedResidualBlocksbb(self.bs)
if full:
pMat = self.trainedModel.data.projMat
self.assembleReducedResidualBlocksAb(self.As, self.bs, pMat)
self.assembleReducedResidualBlocksAA(self.As, pMat)
diff --git a/rrompy/reduction_methods/greedy/rational_interpolant_greedy.py b/rrompy/reduction_methods/greedy/rational_interpolant_greedy.py
index b9d1cad..42bf8eb 100644
--- a/rrompy/reduction_methods/greedy/rational_interpolant_greedy.py
+++ b/rrompy/reduction_methods/greedy/rational_interpolant_greedy.py
@@ -1,311 +1,385 @@
# 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_greedy_approximant import (GenericGreedyApproximant,
localL2Distance)
from rrompy.utilities.poly_fitting.polynomial import (polybases,
PolynomialInterpolator as PI,
polyvanderTotal as pvT,
hashIdxToDerivative as hashI)
from rrompy.reduction_methods.standard import RationalInterpolant
from rrompy.reduction_methods.trained_model import (
TrainedModelRational as tModel)
from rrompy.reduction_methods.trained_model import TrainedModelData
from rrompy.utilities.base.types import (Np1D, Tuple, DictAny, HFEng, paramVal,
paramList)
from rrompy.utilities.base import verbosityManager as vbMng
-from rrompy.utilities.exception_manager import RROMPyWarning
-from rrompy.utilities.exception_manager import RROMPyException, RROMPyAssert
+from rrompy.utilities.poly_fitting import customFit
+from rrompy.utilities.exception_manager import (RROMPyWarning, RROMPyException,
+ RROMPyAssert)
+from rrompy.parameter import checkParameterList
__all__ = ['RationalInterpolantGreedy']
class RationalInterpolantGreedy(GenericGreedyApproximant, RationalInterpolant):
"""
ROM greedy rational interpolant computation for parametric problems.
Args:
HFEngine: HF problem solver.
mu0(optional): Default parameter. Defaults to 0.
approxParameters(optional): Dictionary containing values for main
parameters of approximant. Recognized keys are:
- 'POD': whether to compute POD of snapshots; defaults to True;
- 'S': number of starting training points;
- 'sampler': sample point generator;
- 'radialDirectionalWeights': radial basis weights for interpolant
numerator; defaults to 0, i.e. identity;
- 'greedyTol': uniform error tolerance for greedy algorithm;
defaults to 1e-2;
+ - 'collinearityTol': collinearity tolerance for greedy algorithm;
+ defaults to 0.;
- 'interactive': whether to interactively terminate greedy
algorithm; defaults to False;
- 'maxIter': maximum number of greedy steps; defaults to 1e2;
- 'refinementRatio': ratio of training points to be exhausted
before training set refinement; defaults to 0.2;
- 'nTestPoints': number of test points; defaults to 5e2;
- 'trainSetGenerator': training sample points generator; defaults
to sampler;
- 'polybasis': type of basis for interpolation; defaults to
'MONOMIAL';
- 'errorEstimatorKind': kind of error estimator; available values
- include 'AFFINE', 'INTERPOLATORY', 'LOCALL2', and 'DIAGONAL';
- defaults to 'AFFINE';
+ include 'AFFINE', 'DISCREPANCY', 'INTERPOLATORY', 'LOCALL2',
+ and 'DIAGONAL'; defaults to 'AFFINE';
- 'interpRcond': tolerance for 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.
mus: Array of snapshot parameters.
approxParameters: Dictionary containing values for main parameters of
approximant. Recognized keys are in parameterList.
parameterListSoft: Recognized keys of soft approximant parameters:
- 'POD': whether to compute POD of snapshots.
- 'radialDirectionalWeights': radial basis weights for interpolant
numerator; defaults to 0, i.e. identity;
- 'greedyTol': uniform error tolerance for greedy algorithm;
+ - 'collinearityTol': collinearity tolerance for greedy algorithm;
- 'interactive': whether to interactively terminate greedy
algorithm;
- 'maxIter': maximum number of greedy steps;
- 'refinementRatio': ratio of test points to be exhausted before
test set refinement;
- 'nTestPoints': number of test points;
- 'trainSetGenerator': training sample points generator;
- 'errorEstimatorKind': kind of error estimator;
- 'interpRcond': tolerance for interpolation;
- 'robustTol': tolerance for robust rational denominator
management.
parameterListCritical: Recognized keys of critical approximant
parameters:
- 'S': total number of samples current approximant relies upon;
- 'sampler': sample point generator.
POD: whether to compute POD of snapshots.
S: number of test points.
sampler: Sample point generator.
radialDirectionalWeights: Radial basis weights for interpolant
numerator.
greedyTol: uniform error tolerance for greedy algorithm.
+ collinearityTol: Collinearity tolerance for greedy algorithm.
interactive: whether to interactively terminate greedy algorithm.
maxIter: maximum number of greedy steps.
refinementRatio: ratio of training points to be exhausted before
training set refinement.
nTestPoints: number of starting training points.
trainSetGenerator: training sample points generator.
robustTol: tolerance for robust rational denominator management.
errorEstimatorKind: kind of error estimator.
interpRcond: tolerance for interpolation.
robustTol: tolerance for robust rational denominator management.
muBounds: list of bounds for parameter values.
samplingEngine: Sampling engine.
estimatorNormEngine: Engine for estimator norm computation.
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.
"""
- _allowedEstimatorKinds = ["AFFINE", "INTERPOLATORY", "LOCALL2", "DIAGONAL"]
+ _allowedEstimatorKinds = ["AFFINE", "DISCREPANCY", "INTERPOLATORY",
+ "LOCALL2", "DIAGONAL"]
def __init__(self, HFEngine:HFEng, mu0 : paramVal = None,
approxParameters : DictAny = {}, verbosity : int = 10,
timestamp : bool = True):
self._preInit()
self._addParametersToList(["polybasis", "errorEstimatorKind",
"interpRcond", "robustTol"],
["MONOMIAL", "AFFINE", -1, 0])
super().__init__(HFEngine = HFEngine, mu0 = mu0,
approxParameters = approxParameters,
verbosity = verbosity, timestamp = timestamp)
self._postInit()
@property
def E(self):
"""Value of E."""
- self._E = self._S - 1
+ self._E = self.derivativeShellIdx - 1
return self._E
@E.setter
def E(self, E):
RROMPyWarning(("E is used just to simplify inheritance, and its value "
- "cannot be changed from that of S - 1."))
+ "cannot be changed from that of derivativeShellIdx "
+ "- 1."))
@property
def polybasis(self):
"""Value of polybasis."""
return self._polybasis
@polybasis.setter
def polybasis(self, polybasis):
try:
polybasis = polybasis.upper().strip().replace(" ","")
if polybasis not in polybases:
raise RROMPyException("Sample type not recognized.")
self._polybasis = polybasis
except:
RROMPyWarning(("Prescribed polybasis not recognized. Overriding "
"to 'MONOMIAL'."))
self._polybasis = "MONOMIAL"
self._approxParameters["polybasis"] = self.polybasis
@property
def errorEstimatorKind(self):
"""Value of errorEstimatorKind."""
return self._errorEstimatorKind
@errorEstimatorKind.setter
def errorEstimatorKind(self, errorEstimatorKind):
errorEstimatorKind = errorEstimatorKind.upper()
if errorEstimatorKind not in self._allowedEstimatorKinds:
RROMPyWarning(("Error estimator kind not recognized. Overriding "
"to 'AFFINE'."))
errorEstimatorKind = "AFFINE"
self._errorEstimatorKind = errorEstimatorKind
self._approxParameters["errorEstimatorKind"] = self.errorEstimatorKind
def _errorSamplingRatio(self, mus:Np1D,
muRTest:paramList) -> Tuple[Np1D, Np1D]:
"""Scalar ratio in explicit error estimator."""
- self.setupApprox()
muRTrain = self.trainedModel.centerNormalize(self.mus)
nodalVals = np.prod(localL2Distance(muRTest.data, muRTrain.data),
axis = 1)
denVals = self.trainedModel.getQVal(mus)
- return np.abs(nodalVals), np.abs(denVals)
+ return np.abs(nodalVals), np.abs(denVals)
+
+ def _errorEstimatorDiscrepancy(self, muRTest:Np1D) -> Np1D:
+ """Interpolation discrepancy-based error estimator."""
+ self.assembleReducedResidualBlocks(full = True)
+ nAs = self.trainedModel.data.resAA.shape[1]
+ nbs = self.trainedModel.data.resbb.shape[0]
+ maxDer = np.sum(hashI(max(nAs, nbs) - 2 + self.S, self.npar))
+ evalMus, _, evalMusIdxs = pvT(muRTest, maxDer, self.polybasis0)
+ evalMus = evalMus[:, evalMusIdxs[: max(nAs, nbs) + self.S - 1]].T
+ vSupp, vSuppDer, vSuppIdxs = pvT(self._musUniqueCN, maxDer,
+ self.polybasis0, self._derIdxs,
+ self._reorder)
+ vSuppBase = vSupp[:, vSuppIdxs[: self.S]]
+ vSupp = vSupp[:, vSuppIdxs[self.S : max(nAs, nbs) + self.S - 1]]
+ vSuppDer = vSuppDer[self.S : max(nAs, nbs) + self.S - 1]
+ interCoeffs = customFit(vSuppBase, vSupp,
+ rcond = self.interpRcond).T
+ discr = evalMus[self.S :] - interCoeffs.dot(evalMus[: self.S])
+ scalingA = np.zeros((self.S, nAs, nAs - 1),
+ dtype = self.trainedModel.data.P.coeffs.dtype)
+ scalingb = np.zeros((nbs, nbs - 1),
+ dtype = self.trainedModel.data.Q.coeffs.dtype)
+ for j, derIdxs in enumerate(vSuppDer):
+ for k in range(max(nAs, nbs)):
+ kIdx = hashI(k, self.npar)
+ diffIdx = [x - y for (x, y) in zip(derIdxs, kIdx)]
+ if all([x >= 0 for x in diffIdx]) and max(diffIdx) <= self.E:
+ if k < nAs and j < nAs - 1:
+ scalingA[:, k, j] = self.trainedModel.data.P.coeffs[
+ tuple(diffIdx)]
+ if k < nbs and j < nbs - 1:
+ scalingb[k, j] = self.trainedModel.data.Q.coeffs[
+ tuple(diffIdx)]
+ radiusA = np.tensordot(scalingA, discr[: nAs - 1], 1)
+ radiusb = scalingb.dot(discr[: nbs - 1])
+ # 'ij,jk,ik->k', resbb, radiusb, radiusb.conj()
+ RHSNorm2 = np.sum(self.trainedModel.data.resbb.dot(evalMus[: nbs])
+ * evalMus[: nbs].conj(), axis = 0)
+ # 'ij,jk,ik->k', resbb, radiusb, radiusb.conj()
+ ff = np.sum(self.trainedModel.data.resbb.dot(radiusb) * radiusb.conj(),
+ axis = 0)
+ # 'ijk,jkl,il->l', resAb, radiusA, radiusb.conj()
+ Lf = np.sum(np.tensordot(self.trainedModel.data.resAb, radiusA, 2)
+ * radiusb.conj(), axis = 0)
+ # 'ijkl,klt,ijt->t', resAA, radiusA, radiusA.conj()
+ LL = np.sum(np.tensordot(self.trainedModel.data.resAA, radiusA, 2)
+ * radiusA.conj(), axis = (0, 1))
+ return np.abs((LL - 2. * np.real(Lf) + ff) / RHSNorm2) ** .5
def errorEstimator(self, mus:Np1D) -> Np1D:
"""Standard residual-based error estimator."""
if self.errorEstimatorKind == "AFFINE":
return super().errorEstimator(mus)
- RROMPyAssert(self.HFEngine.npar, 1, "Parameter dimension")
self.setupApprox()
+ mustr = mus
+ if len(mus) > 2:
+ mustr = "[{} ..({}).. {}]".format(mus[0], len(mus) - 2, mus[-1])
+ vbMng(self.trainedModel, "INIT",
+ "Evaluating error estimator at mu = {}.".format(mustr), 10)
+ mus = checkParameterList(mus, self.npar)[0].data
muRTest = self.trainedModel.centerNormalize(mus)
samplingNodal, samplingQ = self._errorSamplingRatio(mus, muRTest)
+ if self.errorEstimatorKind == "DISCREPANCY":
+ return self._errorEstimatorDiscrepancy(muRTest) / samplingQ
samplingRatio = samplingNodal / samplingQ
if self.errorEstimatorKind == "DIAGONAL":
self.assembleReducedResidualBlocks(full = False)
nbs = self.HFEngine.nbs
vanderBase, _, vanderBaseIdxs = pvT(muRTest,
np.sum(hashI(nbs - 1, self.npar)),
"MONOMIAL")
radiusb = vanderBase[:, vanderBaseIdxs][:, : nbs].T
bresb = np.sum(self.trainedModel.data.resbb.dot(radiusb)
* radiusb.conj(), axis = 0)
self.assembleReducedResidualGramian(self.trainedModel.data.projMat)
pDom = self.trainedModel.data.P.domCoeff
LL = pDom.conj().dot(self.trainedModel.data.gramian.dot(pDom))
if not hasattr(self, "Anorm2Approx"):
if self.HFEngine.nAs > 1:
- Adiag = self.scaleFactor[0] * np.diagonal(
- self.HFEngine.A(self.mu0, [1]))
+ Ader = self.HFEngine.A(self.mu0,
+ [1] + [0] * (self.npar - 1))
+ try:
+ Adiag = self.scaleFactor[0] * Ader.diagonal()
+ except:
+ Adiag = self.scaleFactor[0] * np.diagonal(Ader)
self.Anorm2Approx = np.mean(np.abs(Adiag) ** 2.)
if np.isclose(self.Anorm2Approx, 0.) or self.HFEngine.nAs <= 1:
self.Anorm2Approx = 1
jOpt = np.abs(self.Anorm2Approx * LL / bresb) ** .5
err = jOpt * samplingRatio
elif self.errorEstimatorKind == "INTERPOLATORY":
self.initEstimatorNormEngine()
sampRCP = copy(samplingRatio)
idx_muTestSample = np.empty(self.derivativeShellSize, dtype = int)
for j in range(self.derivativeShellSize):
k = np.argmax(sampRCP)
idx_muTestSample[j] = k
if j + 1 < self.derivativeShellSize:
sampRCP *= np.linalg.norm(
- (mus.data ** self.HFEngine.rescalingExp
- - mus[k] ** self.HFEngine.rescalingExp)
- / self.scaleFactor, axis = 1)
+ (mus ** self.HFEngine.rescalingExp
+ - mus[k] ** self.HFEngine.rescalingExp)
+ / self.scaleFactor, axis = 1)
mu_muTestSample = mus[idx_muTestSample]
app_muTestSample = self.getApprox(mu_muTestSample)
resmus = self.HFEngine.residual(app_muTestSample, mu_muTestSample,
duality = False)
RHSmus = self.HFEngine.residual(None, mu_muTestSample,
duality = False)
ressamples = (self.estimatorNormEngine.norm(resmus)
/ self.estimatorNormEngine.norm(RHSmus))
musT = copy(self.mus)
musT.append(mu_muTestSample)
musT = self.trainedModel.centerNormalize(musT)
+ musC = self.trainedModel.centerNormalize(mus)
resT = np.zeros(len(musT), dtype = np.complex)
- resT[len(self.mus) :] = ressamples * samplingQ[idx_muTestSample]
- p = PI()
- wellCond, msg = p.setupByInterpolation(musT, resT, self.M + 1,
- self.polybasis, self.verbosity >= 25,
+ err = np.zeros(len(mus))
+ for l in range(len(mu_muTestSample)):
+ resT[len(self.mus) + l] = (ressamples[l]
+ * samplingQ[idx_muTestSample[l]])
+ p = PI()
+ wellCond, msg = p.setupByInterpolation(musT, resT, self.M + 1,
+ self.polybasis, self.verbosity >= 15,
True, {}, {"rcond": self.interpRcond})
- err = np.abs(p(self.trainedModel.centerNormalize(mus))) / samplingQ
+ err += np.abs(p(musC))
+ resT[len(self.mus) + l] = 0.
+ err /= samplingQ
+ vbMng(self, "MAIN", msg, 15)
else: #if self.errorEstimatorKind == "LOCALL2":
self.initEstimatorNormEngine()
idx_muTestSample = np.argmax(samplingRatio)
mu_muTestSample = mus[idx_muTestSample]
app_muTestSample = self.getApprox(mu_muTestSample)
resmu = self.HFEngine.residual(app_muTestSample, mu_muTestSample,
duality = False)
RHSmu = self.HFEngine.residual(None, mu_muTestSample,
duality = False)
jOpt = np.abs(self.estimatorNormEngine.norm(resmu)[0]
/ samplingRatio[idx_muTestSample]
/ self.estimatorNormEngine.norm(RHSmu)[0])
err = jOpt * samplingRatio
+ vbMng(self.trainedModel, "DEL", "Done evaluating error estimator", 10)
return err
def setupApprox(self, plotEst : bool = False):
"""
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.greedy(plotEst)
self._S = len(self.mus)
- self._N, self._M = [self._S - 1] * 2
+ self._N, self._M = [self.E] * 2
if self.trainedModel is None:
self.trainedModel = tModel()
self.trainedModel.verbosity = self.verbosity
self.trainedModel.timestamp = self.timestamp
data = TrainedModelData(self.trainedModel.name(), self.mu0,
self.samplingEngine.samples,
self.scaleFactor,
self.HFEngine.rescalingExp)
data.mus = copy(self.mus)
self.trainedModel.data = data
else:
self.trainedModel = self.trainedModel
self.trainedModel.data.projMat = copy(self.samplingEngine.samples)
self.trainedModel.data.mus = copy(self.mus)
self.catchInstability = True
if self.N > 0:
Q = self._setupDenominator()[0]
else:
Q = PI()
Q.coeffs = np.ones(1, dtype = np.complex)
Q.npar = 1
Q.polybasis = self.polybasis
self.trainedModel.data.Q = copy(Q)
self.trainedModel.data.P = copy(self._setupNumerator())
self.trainedModel.data.approxParameters = copy(self.approxParameters)
vbMng(self, "DEL", "Done setting up approximant.", 5)
diff --git a/rrompy/reduction_methods/greedy/reduced_basis_greedy.py b/rrompy/reduction_methods/greedy/reduced_basis_greedy.py
index 91fe185..61d1447 100644
--- a/rrompy/reduction_methods/greedy/reduced_basis_greedy.py
+++ b/rrompy/reduction_methods/greedy/reduced_basis_greedy.py
@@ -1,172 +1,176 @@
# 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
from .generic_greedy_approximant import GenericGreedyApproximant
from rrompy.reduction_methods.standard import ReducedBasis
from rrompy.reduction_methods.trained_model import \
TrainedModelReducedBasis as tModel
from rrompy.reduction_methods.trained_model import TrainedModelData
from rrompy.utilities.base.types import DictAny, HFEng, paramVal
from rrompy.utilities.base import verbosityManager as vbMng
from rrompy.utilities.exception_manager import RROMPyWarning, RROMPyAssert
__all__ = ['ReducedBasisGreedy']
class ReducedBasisGreedy(GenericGreedyApproximant, ReducedBasis):
"""
ROM greedy RB approximant computation for parametric problems.
Args:
HFEngine: HF problem solver.
mu0(optional): Default parameter. Defaults to 0.
approxParameters(optional): Dictionary containing values for main
parameters of approximant. Recognized keys are:
- 'POD': whether to compute POD of snapshots; defaults to True;
- 'S': number of starting training points;
- 'sampler': sample point generator;
- 'greedyTol': uniform error tolerance for greedy algorithm;
defaults to 1e-2;
+ - 'collinearityTol': collinearity tolerance for greedy algorithm;
+ defaults to 0.;
- 'interactive': whether to interactively terminate greedy
algorithm; defaults to False;
- 'maxIter': maximum number of greedy steps; defaults to 1e2;
- 'refinementRatio': ratio of test points to be exhausted before
test set refinement; defaults to 0.2;
- 'nTestPoints': number of test points; defaults to 5e2;
- 'trainSetGenerator': training sample points generator; defaults
to sampler.
Defaults to empty dict.
verbosity(optional): Verbosity level. Defaults to 10.
Attributes:
HFEngine: HF problem solver.
mu0: Default parameter.
mus: Array of snapshot parameters.
approxParameters: Dictionary containing values for main parameters of
approximant. Recognized keys are in parameterList.
parameterListSoft: Recognized keys of soft approximant parameters:
- 'POD': whether to compute POD of snapshots.
- 'greedyTol': uniform error tolerance for greedy algorithm;
+ - 'collinearityTol': collinearity tolerance for greedy algorithm;
- 'interactive': whether to interactively terminate greedy
algorithm;
- 'maxIter': maximum number of greedy steps;
- 'refinementRatio': ratio of test points to be exhausted before
test set refinement;
- 'nTestPoints': number of test points;
- 'trainSetGenerator': training sample points generator.
parameterListCritical: Recognized keys of critical approximant
parameters:
- 'S': total number of samples current approximant relies upon;
- 'sampler': sample point generator.
POD: whether to compute POD of snapshots.
S: number of test points.
sampler: Sample point generator.
greedyTol: uniform error tolerance for greedy algorithm.
+ collinearityTol: Collinearity tolerance for greedy algorithm.
interactive: whether to interactively terminate greedy algorithm.
maxIter: maximum number of greedy steps.
refinementRatio: ratio of training points to be exhausted before
training set refinement.
nTestPoints: number of starting training points.
trainSetGenerator: training sample points generator.
muBounds: list of bounds for parameter values.
samplingEngine: Sampling engine.
estimatorNormEngine: Engine for estimator norm computation.
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.
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,
approxParameters : DictAny = {}, verbosity : int = 10,
timestamp : bool = True):
self._preInit()
self._addParametersToList([], [], toBeExcluded = ["R", "PODTolerance"])
super().__init__(HFEngine = HFEngine, mu0 = mu0,
approxParameters = approxParameters,
verbosity = verbosity, timestamp = timestamp)
self._postInit()
@property
def R(self):
"""Value of R."""
self._R = self._S
return self._R
@R.setter
def R(self, R):
RROMPyWarning(("R is used just to simplify inheritance, and its value "
"cannot be changed from that of S."))
@property
def PODTolerance(self):
"""Value of PODTolerance."""
self._PODTolerance = -1
return self._PODTolerance
@PODTolerance.setter
def PODTolerance(self, PODTolerance):
RROMPyWarning(("PODTolerance is used just to simplify inheritance, "
"and its value cannot be changed from -1."))
def setupApprox(self, plotEst : bool = False):
"""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.greedy(plotEst)
vbMng(self, "INIT", "Computing projection matrix.", 7)
pMat = self.samplingEngine.samples
if self.trainedModel is None:
self.trainedModel = tModel()
self.trainedModel.verbosity = self.verbosity
self.trainedModel.timestamp = self.timestamp
data = TrainedModelData(self.trainedModel.name(), self.mu0,
pMat, self.scaleFactor,
self.HFEngine.rescalingExp)
ARBs, bRBs = self.assembleReducedSystem(pMat)
self.trainedModel.data = data
else:
self.trainedModel = self.trainedModel
pMatOld = self.trainedModel.data.projMat
Sold = pMatOld.shape[1]
idxNew = list(range(Sold, pMat.shape[1]))
ARBs, bRBs = self.assembleReducedSystem(pMat(idxNew), pMatOld)
self.trainedModel.data.projMat = copy(pMat)
self.trainedModel.data.mus = copy(self.mus)
self.trainedModel.data.ARBs = ARBs
self.trainedModel.data.bRBs = bRBs
vbMng(self, "DEL", "Done computing projection matrix.", 7)
self.trainedModel.data.approxParameters = copy(self.approxParameters)
vbMng(self, "DEL", "Done setting up approximant.", 5)
diff --git a/rrompy/utilities/poly_fitting/polynomial/polynomial_interpolator.py b/rrompy/utilities/poly_fitting/polynomial/polynomial_interpolator.py
index 21f6de7..356269f 100644
--- a/rrompy/utilities/poly_fitting/polynomial/polynomial_interpolator.py
+++ b/rrompy/utilities/poly_fitting/polynomial/polynomial_interpolator.py
@@ -1,136 +1,136 @@
# 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.utilities.base.types import (List, ListAny, DictAny, Np1D, Np2D,
paramList)
from rrompy.utilities.base import freepar as fp
from rrompy.utilities.poly_fitting.interpolator import GenericInterpolator
from rrompy.utilities.poly_fitting.custom_fit import customFit
from .base import polyfitname, polydomcoeff
from .val import polyval
from .roots import polyroots
from .degree import degreeTotalToFull
from .vander import polyvander as pv, polyvanderTotal as pvT
from rrompy.utilities.exception_manager import RROMPyAssert, RROMPyException
from rrompy.parameter import checkParameterList
__all__ = ['PolynomialInterpolator']
class PolynomialInterpolator(GenericInterpolator):
"""HERE"""
def __init__(self, other = None):
if other is None: return
self.coeffs = other.coeffs
self.npar = other.npar
self.polybasis = other.polybasis
@property
def shape(self):
if self.coeffs.ndim > self.npar:
sh = self.coeffs.shape[self.npar :]
else: sh = tuple([1])
return sh
@property
def deg(self):
return [x - 1 for x in self.coeffs.shape[: self.npar]]
def __getitem__(self, key):
return self.coeffs[key]
def __call__(self, mu:paramList, der : List[int] = None,
scl : Np1D = None):
return polyval(mu, self.coeffs, self.polybasis, der, scl)
def __copy__(self):
return PolynomialInterpolator(self)
def __deepcopy__(self, memo):
other = PolynomialInterpolator()
other.coeffs, other.npar, other.polybasis = copy(
(self.coeffs, self.npar, self.polybasis), memo)
return other
@property
def domCoeff(self):
- RROMPyAssert(self.npar, 1, "Number of parameters")
- return polydomcoeff(self.deg, self.polybasis) * self[-1]
+ return (polydomcoeff(self.deg[0], self.polybasis)
+ * self[(-1,) + (0,) * (self.npar - 1)])
def pad(self, nleft : List[int] = None, nright : List[int] = None):
if nleft is None: nleft = [0] * len(self.shape)
if nright is None: nright = [0] * len(self.shape)
if not hasattr(nleft, "__len__"): nleft = [nleft]
if not hasattr(nright, "__len__"): nright = [nright]
RROMPyAssert(len(self.shape), len(nleft), "Shape of output")
RROMPyAssert(len(self.shape), len(nright), "Shape of output")
padwidth = [(0, 0)] * self.npar
padwidth = padwidth + [(l, r) for l, r in zip(nleft, nright)]
self.coeffs = np.pad(self.coeffs, padwidth, "constant",
constant_values = (0., 0.))
def postmultiplyTensorize(self, A:Np2D):
RROMPyAssert(A.shape[0], self.shape[-1], "Shape of output")
self.coeffs = np.tensordot(self.coeffs, A, axes = (-1, 0))
def setupByInterpolation(self, support:paramList, values:ListAny,
deg:int, polybasis : str = "MONOMIAL",
verbose : bool = True, totalDegree : bool = True,
vanderCoeffs : DictAny = {},
fitCoeffs : DictAny = {}):
support = checkParameterList(support)[0]
self.npar = support.shape[1]
self.polybasis = polybasis
if totalDegree:
vander, _, reorder = pvT(support, deg, basis = polybasis,
**vanderCoeffs)
vander = vander[:, reorder]
else:
if not hasattr(deg, "__len__"): deg = [deg] * self.npar
vander = pv(support, deg, basis = polybasis, **vanderCoeffs)
outDim = values.shape[1:]
values = values.reshape(values.shape[0], -1)
fitOut = customFit(vander, values, full = True, **fitCoeffs)
P = fitOut[0]
if verbose:
msg = ("Fitting {} samples with degree {} through {}... "
"Conditioning of LS system: {:.4e}.").format(
len(vander), deg,
polyfitname(self.polybasis),
fitOut[1][2][0] / fitOut[1][2][-1])
else: msg = None
if totalDegree:
self.coeffs = degreeTotalToFull(tuple([deg + 1] * self.npar)
+ outDim, self.npar, P)
else:
self.coeffs = P.reshape(tuple([d + 1 for d in deg]) + outDim)
return fitOut[1][1] == vander.shape[1], msg
def roots(self, marginalVals : ListAny = [fp]):
RROMPyAssert(self.shape, (1,), "Shape of output")
RROMPyAssert(len(marginalVals), self.npar, "Number of parameters")
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."))
return polyroots(self.coeffs, self.polybasis, marginalVals)
diff --git a/rrompy/utilities/poly_fitting/radial_basis/radial_basis_interpolator.py b/rrompy/utilities/poly_fitting/radial_basis/radial_basis_interpolator.py
index 9f16561..f9731cb 100644
--- a/rrompy/utilities/poly_fitting/radial_basis/radial_basis_interpolator.py
+++ b/rrompy/utilities/poly_fitting/radial_basis/radial_basis_interpolator.py
@@ -1,137 +1,139 @@
# 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.utilities.base.types import (List, ListAny, DictAny, Np1D, Np2D,
paramList)
from rrompy.utilities.poly_fitting.interpolator import GenericInterpolator
from rrompy.utilities.poly_fitting.custom_fit import customFit
from .base import polyfitname
from .val import polyval
from .vander import polyvander as pv, polyvanderTotal as pvT
from rrompy.utilities.poly_fitting.polynomial.degree import degreeTotalToFull
from rrompy.utilities.exception_manager import RROMPyException, RROMPyAssert
from rrompy.parameter import checkParameterList
__all__ = ['RadialBasisInterpolator']
class RadialBasisInterpolator(GenericInterpolator):
"""HERE"""
def __init__(self, other = None):
if other is None: return
self.support = other.support
self.coeffsGlobal = other.coeffsGlobal
self.coeffsLocal = other.coeffsLocal
self.directionalWeights = other.directionalWeights
self.npar = other.npar
self.polybasis = other.polybasis
@property
def shape(self):
sh = self.coeffsLocal.shape[1 :] if self.coeffsLocal.ndim > 1 else 1
return sh
@property
def deg(self):
return [x - 1 for x in self.coeffsGlobal.shape[: self.npar]]
def __call__(self, mu:paramList, der : List[int] = None,
scl : Np1D = None):
if der is not None and np.sum(der) > 0:
raise RROMPyException(("Cannot take derivatives of radial basis "
"function."))
return polyval(mu, self.coeffsGlobal, self.coeffsLocal, self.support,
self.directionalWeights, self.polybasis)
def __copy__(self):
return RadialBasisInterpolator(self)
def __deepcopy__(self, memo):
other = RadialBasisInterpolator()
(other.support, other.coeffsGlobal, other.coeffsLocal,
other.directionalWeights, other.npar, other.polybasis) = copy(
(self.support, self.coeffsGlobal, self.coeffsLocal,
self.directionalWeights, self.npar, self.polybasis), memo)
return other
def postmultiplyTensorize(self, A:Np2D):
RROMPyAssert(A.shape[0], self.shape[-1], "Shape of output")
self.coeffsLocal = np.tensordot(self.coeffsLocal, A, axes = (-1, 0))
self.coeffsGlobal = np.tensordot(self.coeffsGlobal, A, axes = (-1, 0))
def pad(self, nleft : List[int] = None, nright : List[int] = None):
if nleft is None: nleft = [0] * len(self.shape)
if nright is None: nright = [0] * len(self.shape)
if not hasattr(nleft, "__len__"): nleft = [nleft]
if not hasattr(nright, "__len__"): nright = [nright]
RROMPyAssert(len(self.shape), len(nleft), "Shape of output")
RROMPyAssert(len(self.shape), len(nright), "Shape of output")
padwidth = [(0, 0)] + [(l, r) for l, r in zip(nleft, nright)]
self.coeffsLocal = np.pad(self.coeffsLocal, padwidth, "constant",
constant_values = (0., 0.))
padwidth = [(0, 0)] * (self.npar - 1) + padwidth
self.coeffsGlobal = np.pad(self.coeffsGlobal, padwidth, "constant",
constant_values = (0., 0.))
def setupByInterpolation(self, support:paramList, values:ListAny,
deg:int, polybasis : str = "MONOMIAL_GAUSSIAN",
directionalWeights : Np1D = None,
verbose : bool = True, totalDegree : bool = True,
vanderCoeffs : DictAny = {},
fitCoeffs : DictAny = {}):
support = checkParameterList(support)[0]
self.support = copy(support)
if "reorder" in vanderCoeffs.keys():
self.support = self.support[vanderCoeffs["reorder"]]
self.npar = support.shape[1]
if directionalWeights is None:
directionalWeights = np.ones(self.npar)
self.directionalWeights = directionalWeights
self.polybasis = polybasis
if totalDegree:
vander, _, reorder = pvT(support, deg, basis = polybasis,
directionalWeights = self.directionalWeights,
**vanderCoeffs)
vander = vander[reorder]
vander = vander[:, reorder]
else:
if not hasattr(deg, "__len__"): deg = [deg] * self.npar
- vander = pv(support, deg, basis = polybasis, **vanderCoeffs)
+ vander = pv(support, deg, basis = polybasis,
+ directionalWeights = self.directionalWeights,
+ **vanderCoeffs)
outDim = values.shape[1:]
values = values.reshape(values.shape[0], -1)
values = np.pad(values, ((0, len(vander) - len(values)), (0, 0)),
"constant")
fitOut = customFit(vander, values, full = True, **fitCoeffs)
P = fitOut[0][len(support) :]
if verbose:
msg = ("Fitting {}+{} samples with degree {} through {}... "
"Conditioning of LS system: {:.4e}.").format(
len(support), len(vander) - len(support),
deg, polyfitname(self.polybasis),
fitOut[1][2][0] / fitOut[1][2][-1])
else: msg = None
self.coeffsLocal = fitOut[0][: len(support)]
if totalDegree:
self.coeffsGlobal = degreeTotalToFull(tuple([deg + 1] * self.npar)
+ outDim, self.npar, P)
else:
self.coeffsGlobal = P.reshape(tuple([d + 1 for d in deg]) + outDim)
return fitOut[1][1] == vander.shape[1], msg