diff --git a/examples/2d/pod/cookie_single_pod.py b/examples/2d/pod/cookie_single_pod.py
index 6cc9f64..acd234f 100644
--- a/examples/2d/pod/cookie_single_pod.py
+++ b/examples/2d/pod/cookie_single_pod.py
@@ -1,132 +1,131 @@
import numpy as np
from rrompy.hfengines.linear_problem.bidimensional import \
CookieEngineSingle as CES
from rrompy.reduction_methods.standard import RationalInterpolant as RI
from rrompy.reduction_methods.standard import ReducedBasis as RB
from rrompy.parameter.parameter_sampling import (QuadratureSampler as QS,
QuadratureSamplerTotal as QST,
ManualSampler as MS,
RandomSampler as RS)
verb = 5
size = 1
show_sample = True
show_norm = True
clip = -1
#clip = .4
#clip = .6
Delta = -10
MN = 15
R = (MN + 2) * (MN + 1) // 2
S = [int(np.ceil(R ** .5))] * 2
PODTol = 1e-6
samples = "centered"
samples = "standard"
algo = "rational"
#algo = "RB"
sampling = "quadrature"
sampling = "quadrature_total"
sampling = "random"
if samples == "standard":
radial = ""
# radial = "_gaussian"
# radial = "_thinplate"
# radial = "_multiquadric"
rW0 = 1.
radialWeight = [rW0] * 2
assert Delta <= 0
if size == 1: # below
mu0 = [20 ** .5, 1. ** .5]
mutar = [20.5 ** .5, 1.05 ** .5]
murange = [[18.5 ** .5, .85 ** .5], [21.5 ** .5, 1.15 ** .5]]
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]]]
kappa = 20. ** .5
theta = - np.pi / 6.
n = 30
Rad = 1.
L = np.pi
nX = 2
nY = 1
solver = CES(kappa = kappa, theta = theta, n = n, R = Rad, L = L, nX = nX,
nY = nY, mu0 = mu0, verbosity = verb)
rescalingExp = [2.] * 2
if algo == "rational":
params = {'N':MN, 'M':MN + Delta, 'S':S, 'POD':True}
if samples == "standard":
params['polybasis'] = "CHEBYSHEV" + radial
# params['polybasis'] = "LEGENDRE" + radial
# params['polybasis'] = "MONOMIAL" + radial
- params['E'] = MN
params['radialDirectionalWeights'] = radialWeight
elif samples == "centered":
params['polybasis'] = "MONOMIAL"
params['S'] = R
method = RI
else: #if algo == "RB":
params = {'R':(MN + 2 + Delta) * (MN + 1 + Delta) // 2, 'S':S,
'POD':True, 'PODTolerance':PODTol}
if samples == "centered":
params['S'] = R
method = RB
if samples == "standard":
if sampling == "quadrature":
params['sampler'] = QS(murange, "CHEBYSHEV", scalingExp = rescalingExp)
# params['sampler'] = QS(murange, "GAUSSLEGENDRE", scalingExp = rescalingExp)
# params['sampler'] = QS(murange, "UNIFORM", scalingExp = rescalingExp)
params['S'] = [max(j, MN + 1) for j in params['S']]
elif sampling == "quadrature_total":
params['sampler'] = QST(murange, "CHEBYSHEV", scalingExp = rescalingExp)
params['S'] = R
else: # if sampling == "random":
params['sampler'] = RS(murange, "HALTON", scalingExp = rescalingExp)
params['S'] = R
elif samples == "centered":
params['sampler'] = MS(murange, points = [mu0], scalingExp = rescalingExp)
approx = method(solver, mu0 = mu0, approxParameters = params, verbosity = verb)
if samples == "standard": 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)))
if algo == "rational" and approx.N > 0:
from plot_zero_set import plotZeroSet2
muZeroVals, Qvals = plotZeroSet2(murange, murangeEff, approx, mu0,
200, [2., 2.], clip = clip)
if show_norm:
solver._solveBatchSize = 100
from plot_inf_set import plotInfSet2
muInfVals, normEx, normApp, normRes, normErr, beta = plotInfSet2(
murange, murangeEff, approx, mu0, 25,
[2., 2.], clip = clip, relative = False)
diff --git a/examples/2d/pod/fracture_pod.py b/examples/2d/pod/fracture_pod.py
index 7be470b..99c6baf 100644
--- a/examples/2d/pod/fracture_pod.py
+++ b/examples/2d/pod/fracture_pod.py
@@ -1,274 +1,275 @@
import numpy as np
from rrompy.hfengines.linear_problem.bidimensional import \
MembraneFractureEngine as MFE
from rrompy.reduction_methods.standard import RationalInterpolant as RI
from rrompy.reduction_methods.standard import ReducedBasis as RB
from rrompy.reduction_methods.pivoted import RationalInterpolantPivoted as RIP
from rrompy.reduction_methods.pivoted import ReducedBasisPivoted as RBP
from rrompy.reduction_methods.pole_matching import \
RationalInterpolantPoleMatching as RIPM
from rrompy.reduction_methods.pole_matching import \
ReducedBasisPoleMatching as RBPM
from rrompy.parameter.parameter_sampling import (QuadratureSampler as QS,
QuadratureSamplerTotal as QST,
ManualSampler as MS,
RandomSampler as RS)
verb = 70
size = 2
-show_sample = False
+show_sample = True
show_norm = True
Delta = 0
-MN = 8
+MN = 5
R = (MN + 2) * (MN + 1) // 2
S = [int(np.ceil(R ** .5))] * 2
-PODTol = 1e-4
+PODTol = 1e-6
SPivot = [MN + 1, 4]
MMarginal = SPivot[1] - 1
matchingWeight = 1.
cutOffTolerance = 1. * np.inf
cutOffType = "MAGNITUDE"
samples = "centered"
samples = "standard"
samples = "pivoted"
samples = "pole matching"
algo = "rational"
-#algo = "RB"
+algo = "RB"
sampling = "quadrature"
sampling = "quadrature_total"
#sampling = "random"
samplingM = "quadrature"
#samplingM = "quadrature_total"
#samplingM = "random"
+polydegreetype = "TOTAL"
+polydegreetype = "FULL"
if samples in ["standard", "pivoted", "pole matching"]:
radial = ""
# radial = "_gaussian"
# radial = "_thinplate"
# radial = "_multiquadric"
rW0 = 5.
radialWeight = [rW0]
if samples in ["pivoted", "pole matching"]:
radialM = ""
# radialM = "_gaussian"
# radialM = "_thinplate"
# radialM = "_multiquadric"
rMW0 = 2.
radialWeightM = [rMW0]
assert Delta <= 0
if size == 1: # below
mu0 = [40 ** .5, .4]
mutar = [45 ** .5, .4]
murange = [[30 ** .5, .3], [50 ** .5, .5]]
elif size == 2: # top
mu0 = [40 ** .5, .6]
mutar = [45 ** .5, .6]
murange = [[30 ** .5, .5], [50 ** .5, .7]]
elif size == 3: # interesting
mu0 = [40 ** .5, .5]
mutar = [45 ** .5, .5]
murange = [[30 ** .5, .3], [50 ** .5, .7]]
elif size == 4: # wide_low
mu0 = [40 ** .5, .2]
mutar = [45 ** .5, .2]
murange = [[10 ** .5, .1], [70 ** .5, .3]]
elif size == 5: # wide_hi
mu0 = [40 ** .5, .8]
mutar = [45 ** .5, .8]
murange = [[10 ** .5, .7], [70 ** .5, .9]]
elif size == 6: # top_zoom
mu0 = [50 ** .5, .8]
mutar = [55 ** .5, .8]
murange = [[40 ** .5, .7], [60 ** .5, .9]]
elif size == 7: # huge
mu0 = [50 ** .5, .5]
mutar = [55 ** .5, .5]
murange = [[10 ** .5, .2], [90 ** .5, .8]]
elif size == 11: #L below
mu0 = [110 ** .5, .4]
mutar = [115 ** .5, .4]
murange = [[90 ** .5, .3], [130 ** .5, .5]]
elif size == 12: #L top
mu0 = [110 ** .5, .6]
mutar = [115 ** .5, .6]
murange = [[90 ** .5, .5], [130 ** .5, .7]]
elif size == 13: #L interesting
mu0 = [110 ** .5, .5]
mutar = [115 ** .5, .5]
murange = [[90 ** .5, .3], [130 ** .5, .7]]
elif size == 14: #L belowbelow
mu0 = [110 ** .5, .2]
mutar = [115 ** .5, .2]
murange = [[90 ** .5, .1], [130 ** .5, .3]]
elif size == 15: #L toptop
mu0 = [110 ** .5, .8]
mutar = [115 ** .5, .8]
murange = [[90 ** .5, .7], [130 ** .5, .9]]
elif size == 16: #L interestinginteresting
mu0 = [110 ** .5, .5]
mutar = [115 ** .5, .6]
murange = [[90 ** .5, .1], [130 ** .5, .9]]
elif size == 17: #L interestingtop
mu0 = [110 ** .5, .7]
mutar = [115 ** .5, .6]
murange = [[90 ** .5, .5], [130 ** .5, .9]]
elif size == 18: #L interestingbelow
mu0 = [110 ** .5, .3]
mutar = [115 ** .5, .4]
murange = [[90 ** .5, .1], [130 ** .5, .5]]
elif size == 100: # tiny
mu0 = [32.5 ** .5, .5]
mutar = [34 ** .5, .5]
murange = [[30 ** .5, .3], [35 ** .5, .7]]
aEff = 1.#25
bEff = 1. - aEff
murangeEff = [[(aEff*murange[0][0]**2. + bEff*murange[1][0]**2.) ** .5,
aEff*murange[0][1] + bEff*murange[1][1]],
[(aEff*murange[1][0]**2. + bEff*murange[0][0]**2.) ** .5,
aEff*murange[1][1] + bEff*murange[0][1]]]
H = 1.
L = .75
delta = .05
n = 20
solver = MFE(mu0 = mu0, H = H, L = L, delta = delta, n = n, verbosity = verb)
rescalingExp = [2., 1.]
if algo == "rational":
- params = {'N':MN, 'M':MN + Delta, 'S':S, 'POD':True}
+ params = {'N':MN, 'M':MN + Delta, 'S':S, 'POD':True,
+ 'polydegreetype':polydegreetype}
if samples == "standard":
params['polybasis'] = "CHEBYSHEV" + radial
# params['polybasis'] = "LEGENDRE" + radial
# params['polybasis'] = "MONOMIAL" + radial
- params['E'] = MN
params['radialDirectionalWeights'] = radialWeight
method = RI
elif samples == "centered":
params['polybasis'] = "MONOMIAL"
params['S'] = R
method = RI
elif samples in ["pivoted", "pole matching"]:
params['S'] = [SPivot[0]]
- params['R'] = SPivot[0]
params['SMarginal'] = [SPivot[1]]
params['MMarginal'] = MMarginal
params['polybasisPivot'] = "CHEBYSHEV" + radial
params['polybasisMarginal'] = "MONOMIAL" + radialM
params['radialDirectionalWeightsPivot'] = radialWeight
params['radialDirectionalWeightsMarginal'] = radialWeightM
if samples == "pivoted":
method = RIP
else:
params['matchingWeight'] = matchingWeight
params['cutOffTolerance'] = cutOffTolerance
params["cutOffType"] = cutOffType
method = RIPM
else: #if algo == "RB":
params = {'R':(MN + 2 + Delta) * (MN + 1 + Delta) // 2, 'S':S,
'POD':True, 'PODTolerance':PODTol}
if samples == "standard":
method = RB
elif samples == "centered":
params['S'] = R
method = RB
elif samples in ["pivoted", "pole matching"]:
params['S'] = [SPivot[0]]
params['R'] = SPivot[0]
params['SMarginal'] = [SPivot[1]]
params['MMarginal'] = MMarginal
params['polybasisMarginal'] = "MONOMIAL" + radialM
params['radialDirectionalWeightsMarginal'] = radialWeightM
method = RBP
if samples == "pivoted":
method = RBP
else:
params['matchingWeight'] = matchingWeight
params['cutOffTolerance'] = cutOffTolerance
params["cutOffType"] = cutOffType
method = RBPM
if samples == "standard":
if sampling == "quadrature":
params['sampler'] = QS(murange, "CHEBYSHEV", scalingExp = rescalingExp)
# params['sampler'] = QS(murange, "GAUSSLEGENDRE", scalingExp = rescalingExp)
# params['sampler'] = QS(murange, "UNIFORM", scalingExp = rescalingExp)
params['S'] = [max(j, MN + 1) for j in params['S']]
elif sampling == "quadrature_total":
params['sampler'] = QST(murange, "CHEBYSHEV", scalingExp = rescalingExp)
params['S'] = R
else: # if sampling == "random":
params['sampler'] = RS(murange, "HALTON", scalingExp = rescalingExp)
params['S'] = R
elif samples == "centered":
params['sampler'] = MS(murange, points = [mu0], scalingExp = rescalingExp)
elif samples in ["pivoted", "pole matching"]:
if sampling == "quadrature":
params['samplerPivot'] = QS([murange[0][0], murange[1][0]], "CHEBYSHEV")
# params['samplerPivot'] = QS([murange[0][0], murange[1][0]], "GAUSSLEGENDRE")
# params['samplerPivot'] = QS([murange[0][0], murange[1][0]], "UNIFORM")
elif sampling == "quadrature_total":
params['samplerPivot'] = QST([murange[0][0], murange[1][0]], "CHEBYSHEV")
else: # if sampling == "random":
params['samplerPivot'] = RS([murange[0][0], murange[1][0]], "HALTON")
if samplingM == "quadrature":
params['samplerMarginal'] = QS([murange[0][1], murange[1][1]], "UNIFORM")
elif samplingM == "quadrature_total":
params['samplerMarginal'] = QST([murange[0][1], murange[1][1]], "CHEBYSHEV")
else: # if samplingM == "random":
params['samplerMarginal'] = RS([murange[0][1], murange[1][1]], "HALTON")
if samples not in ["pivoted", "pole matching"]:
approx = method(solver, mu0 = mu0, approxParameters = params,
verbosity = verb)
else:
approx = method(solver, mu0 = mu0, directionPivot = [0],
approxParameters = params, verbosity = verb)
if samples != "centered": approx.samplingEngine.allowRepeatedSamples = False
approx.setupApprox()
if show_sample:
import ufl
import fenics as fen
from rrompy.solver.fenics import affine_warping
L = mutar[1]
y = fen.SpatialCoordinate(solver.V.mesh())[1]
warp1, warpI1 = affine_warping(solver.V.mesh(),
np.array([[1, 0], [0, 2. * L]]))
warp2, warpI2 = affine_warping(solver.V.mesh(),
np.array([[1, 0], [0, 2. - 2. * L]]))
warp = ufl.conditional(ufl.ge(y, 0.), warp1, warp2)
warpI = ufl.conditional(ufl.ge(y, 0.), warpI1, warpI2)
approx.plotApprox(mutar, [warp, warpI], name = 'u_app', what = "REAL")
approx.plotHF(mutar, [warp, warpI], name = 'u_HF', what = "REAL")
approx.plotErr(mutar, [warp, warpI], name = 'err', what = "REAL")
# approx.plotRes(mutar, [warp, warpI], name = 'res', what = "REAL")
appErr = approx.normErr(mutar)
solNorm = approx.normHF(mutar)
resNorm = approx.normRes(mutar)
RHSNorm = approx.normRHS(mutar)
print(('SolNorm:\t{}\nErr:\t{}\nErrRel:\t{}').format(solNorm, appErr,
np.divide(appErr, solNorm)))
print(('RHSNorm:\t{}\nRes:\t{}\nResRel:\t{}').format(RHSNorm, resNorm,
np.divide(resNorm, RHSNorm)))
approx.verbosity = 5
from plot_zero_set import plotZeroSet2
muZeroVals, Qvals = plotZeroSet2(murange, murangeEff, approx, mu0,
200, [2., 1.])
if show_norm:
solver._solveBatchSize = 100
from plot_inf_set import plotInfSet2
muInfVals, normEx, normApp, normRes, normErr, beta = plotInfSet2(
murange, murangeEff, approx, mu0, 50,
[2., 1.], relative = True,
nobeta = True)
print(np.sort(approx.getPoles([None, .5]) ** 2.))
diff --git a/examples/2d/pod/fracture_pod_nodomain.py b/examples/2d/pod/fracture_pod_nodomain.py
index 4117c2d..d0dfd6c 100644
--- a/examples/2d/pod/fracture_pod_nodomain.py
+++ b/examples/2d/pod/fracture_pod_nodomain.py
@@ -1,167 +1,166 @@
import numpy as np
from rrompy.hfengines.linear_problem import MembraneFractureEngineNoDomain \
as MFEND
from rrompy.reduction_methods.standard import RationalInterpolant as RI
from rrompy.reduction_methods.standard import ReducedBasis as RB
from rrompy.parameter.parameter_sampling import (QuadratureSampler as QS,
QuadratureSamplerTotal as QST,
ManualSampler as MS,
RandomSampler as RS)
verb = 5
size = 1
show_sample = True
show_norm = True
ignore_forcing = True
ignore_forcing = False
clip = -1
#clip = .4
#clip = .6
homogeneize = False
#homogeneize = True
Delta = 0
MN = 6
R = MN + 1
S = R
samples = "centered"
samples = "standard"
algo = "rational"
#algo = "RB"
sampling = "quadrature"
sampling = "quadrature_total"
sampling = "random"
if samples == "standard":
radial = ""
# radial = "_gaussian"
# radial = "_thinplate"
# radial = "_multiquadric"
radialWeight = [2.]
assert Delta <= 0
if size == 1: # below
mu0Aug = [40 ** .5, .4]
mu0 = mu0Aug[0]
mutar = 45 ** .5
murange = [[30 ** .5], [50 ** .5]]
elif size == 2: # top
mu0Aug = [40 ** .5, .6]
mu0 = mu0Aug[0]
mutar = 45 ** .5
murange = [[30 ** .5], [50 ** .5]]
elif size == 3: # interesting
mu0Aug = [40 ** .5, .5]
mu0 = mu0Aug[0]
mutar = 45 ** .5
murange = [[30 ** .5], [50 ** .5]]
elif size == 4: # wide_low
mu0Aug = [40 ** .5, .2]
mu0 = mu0Aug[0]
mutar = 45 ** .5
murange = [[10 ** .5], [70 ** .5]]
elif size == 5: # wide_hi
mu0Aug = [40 ** .5, .8]
mu0 = mu0Aug[0]
mutar = 45 ** .5
murange = [[10 ** .5], [70 ** .5]]
elif size == 6: # top_zoom
mu0Aug = [50 ** .5, .8]
mu0 = mu0Aug[0]
mutar = 55 ** .5
murange = [[40 ** .5], [60 ** .5]]
aEff = 1.#25
bEff = 1. - aEff
murangeEff = [[(aEff*murange[0][0]**2. + bEff*murange[1][0]**2.) ** .5],
[(aEff*murange[1][0]**2. + bEff*murange[0][0]**2.) ** .5]]
H = 1.
L = .75
delta = .05
n = 20
solver = MFEND(mu0 = mu0Aug, H = H, L = L, delta = delta, n = n,
verbosity = verb, homogeneized = homogeneize)
rescalingExp = 2.
if algo == "rational":
params = {'N':MN, 'M':MN + Delta, 'S':S, 'POD':True}
if samples == "standard":
params['polybasis'] = "CHEBYSHEV" + radial
# params['polybasis'] = "LEGENDRE" + radial
# params['polybasis'] = "MONOMIAL" + radial
- params['E'] = MN
params['radialDirectionalWeights'] = radialWeight
elif samples == "centered":
params['polybasis'] = "MONOMIAL"
params['S'] = R
method = RI
else: #if algo == "RB":
params = {'R':R, 'S':S, 'POD':True}
if samples == "centered":
params['S'] = R
method = RB
if samples == "standard":
if sampling == "quadrature":
params['sampler'] = QS(murange, "CHEBYSHEV", scalingExp = rescalingExp)
# params['sampler'] = QS(murange, "GAUSSLEGENDRE", scalingExp = rescalingExp)
# params['sampler'] = QS(murange, "UNIFORM", scalingExp = rescalingExp)
params['S'] = [max(j, MN + 1) for j in params['S']]
elif sampling == "quadrature_total":
params['sampler'] = QST(murange, "CHEBYSHEV", scalingExp = rescalingExp)
params['S'] = R
else: # if sampling == "random":
params['sampler'] = RS(murange, "HALTON", scalingExp = rescalingExp)
params['S'] = R
elif samples == "centered":
params['sampler'] = MS(murange, points = [mu0], scalingExp = rescalingExp)
approx = method(solver, mu0 = mu0, approxParameters = params,
verbosity = verb)
if samples == "standard": approx.samplingEngine.allowRepeatedSamples = False
approx.setupApprox()
if show_sample:
import ufl
import fenics as fen
from rrompy.solver.fenics import affine_warping
L = solver.lFrac
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)))
if algo == "rational":
from plot_zero_set import plotZeroSet1
muZeroVals, Qvals = plotZeroSet1(murange, murangeEff, approx, mu0,
1000, 2.)
if show_norm:
solver._solveBatchSize = 10
from plot_inf_set import plotInfSet1
muInfVals, normEx, normApp, normRes, normErr, beta = plotInfSet1(
murange, murangeEff, approx, mu0,
200, 2., relative = False,
normalizeDen = True)
print(approx.getPoles() ** 2.)
diff --git a/examples/2d/pod/matrix_passive_pod.py b/examples/2d/pod/matrix_passive_pod.py
index d705333..8a5fa82 100644
--- a/examples/2d/pod/matrix_passive_pod.py
+++ b/examples/2d/pod/matrix_passive_pod.py
@@ -1,204 +1,203 @@
import numpy as np
from rrompy.hfengines.linear_problem.tridimensional import \
MatrixDynamicalPassive as MDP
from rrompy.reduction_methods.standard import RationalInterpolant as RI
from rrompy.reduction_methods.standard import ReducedBasis as RB
from rrompy.reduction_methods.pivoted import RationalInterpolantPivoted as RIP
from rrompy.reduction_methods.pivoted import ReducedBasisPivoted as RBP
from rrompy.reduction_methods.pole_matching import \
RationalInterpolantPoleMatching as RIPM
from rrompy.reduction_methods.pole_matching import \
ReducedBasisPoleMatching as RBPM
from rrompy.parameter.parameter_sampling import (QuadratureSampler as QS,
QuadratureSamplerTotal as QST,
ManualSampler as MS,
RandomSampler as RS)
verb = 10
size = 3
show_sample = True
show_norm = True
Delta = 0
MN = 7
R = (MN + 2) * (MN + 1) // 2
S = [int(np.ceil(R ** .5))] * 2
PODTol = 1e-6
SPivot = [MN + 1, 3]
MMarginal = SPivot[1] - 1
samples = "centered"
samples = "standard"
samples = "pivoted"
samples = "pole matching"
algo = "rational"
#algo = "RB"
sampling = "quadrature"
#sampling = "quadrature_total"
#sampling = "random"
samplingM = "quadrature"
#samplingM = "quadrature_total"
#samplingM = "random"
if samples in ["standard", "pivoted", "pole matching"]:
radial = ""
# radial = "_gaussian"
# radial = "_thinplate"
# radial = "_multiquadric"
rW0 = 10.
radialWeight = [rW0]
if samples in ["pivoted", "pole matching"]:
radialM = ""
# radialM = "_gaussian"
# radialM = "_thinplate"
# radialM = "_multiquadric"
rMW0 = 2.
radialWeightM = [rMW0]
matchingWeight = 1.
cutOffTolerance = 5.
cutOffType = "POTENTIAL"
if size == 1:
mu0 = [2.7e2, 20]
mutar = [3e2, 25]
murange = [[20., 10], [5.2e2, 30]]
elif size == 2:
mu0 = [2.7e2, 60]
mutar = [3e2, 75]
murange = [[20., 10], [5.2e2, 110]]
elif size == 3:
mu0 = [2.7e2, 160]
mutar = [3e2, 105]
murange = [[20., 10], [5.2e2, 310]]
assert Delta <= 0
aEff = 1.#25
bEff = 1. - aEff
murangeEff = [[aEff*murange[0][0] + bEff*murange[1][0],
aEff*murange[0][1] + bEff*murange[1][1]],
[aEff*murange[1][0] + bEff*murange[0][0],
aEff*murange[1][1] + bEff*murange[0][1]]]
n = 100
b = 10
solver = MDP(mu0 = mu0, n = n, b = b, verbosity = verb)
if algo == "rational":
params = {'N':MN, 'M':MN + Delta, 'S':S, 'POD':True}
if samples == "standard":
params['polybasis'] = "CHEBYSHEV" + radial
# params['polybasis'] = "LEGENDRE" + radial
# params['polybasis'] = "MONOMIAL" + radial
- params['E'] = MN
params['radialDirectionalWeights'] = radialWeight
method = RI
elif samples == "centered":
params['polybasis'] = "MONOMIAL"
params['S'] = R
method = RI
elif samples in ["pivoted", "pole matching"]:
params['S'] = [SPivot[0]]
params['SMarginal'] = [SPivot[1]]
params['MMarginal'] = MMarginal
params['polybasisPivot'] = "CHEBYSHEV" + radial
params['polybasisMarginal'] = "MONOMIAL" + radialM
params['radialDirectionalWeightsPivot'] = radialWeight
params['radialDirectionalWeightsMarginal'] = radialWeightM
if samples == "pivoted":
method = RIP
else:
params['matchingWeight'] = matchingWeight
params['cutOffTolerance'] = cutOffTolerance
params["cutOffType"] = cutOffType
method = RIPM
else: #if algo == "RB":
params = {'R':(MN + 2 + Delta) * (MN + 1 + Delta) // 2, 'S':S,
'POD':True, 'PODTolerance':PODTol}
if samples == "standard":
method = RB
elif samples == "centered":
params['S'] = R
method = RB
elif samples in ["pivoted", "pole matching"]:
params['R'] = SPivot[0]
params['S'] = [SPivot[0]]
params['SMarginal'] = [SPivot[1]]
params['MMarginal'] = MMarginal
params['polybasisMarginal'] = "MONOMIAL" + radialM
params['radialDirectionalWeightsMarginal'] = radialWeightM
if samples == "pivoted":
method = RBP
else:
params['matchingWeight'] = matchingWeight
params['cutOffTolerance'] = cutOffTolerance
params["cutOffType"] = cutOffType
method = RBPM
if samples == "standard":
if sampling == "quadrature":
params['sampler'] = QS(murange, "CHEBYSHEV")
# params['sampler'] = QS(murange, "GAUSSLEGENDRE")
# params['sampler'] = QS(murange, "UNIFORM")
elif sampling == "quadrature_total":
params['sampler'] = QST(murange, "CHEBYSHEV")
else: # if sampling == "random":
params['sampler'] = RS(murange, "HALTON")
elif samples == "centered":
params['sampler'] = MS(murange, points = [mu0])
elif samples in ["pivoted", "pole matching"]:
if sampling == "quadrature":
params['samplerPivot'] = QS([murange[0][0], murange[1][0]], "CHEBYSHEV")
# params['samplerPivot'] = QS([murange[0][0], murange[1][0]], "GAUSSLEGENDRE")
# params['samplerPivot'] = QS([murange[0][0], murange[1][0]], "UNIFORM")
elif sampling == "quadrature_total":
params['samplerPivot'] = QST([murange[0][0], murange[1][0]], "CHEBYSHEV")
else: # if sampling == "random":
params['samplerPivot'] = RS([murange[0][0], murange[1][0]], "HALTON")
if samplingM == "quadrature":
params['samplerMarginal'] = QS([murange[0][1], murange[1][1]], "UNIFORM")
elif samplingM == "quadrature_total":
params['samplerMarginal'] = QST([murange[0][1], murange[1][1]], "CHEBYSHEV")
else: # if samplingM == "random":
params['samplerMarginal'] = RS([murange[0][1], murange[1][1]], "HALTON")
if samples not in ["pivoted", "pole matching"]:
approx = method(solver, mu0 = mu0, approxParameters = params,
verbosity = verb)
else:
approx = method(solver, mu0 = mu0, directionPivot = [0],
approxParameters = params, verbosity = verb)
if samples != "centered": approx.samplingEngine.allowRepeatedSamples = False
approx.setupApprox()
if show_sample:
L = mutar[1]
approx.plotApprox(mutar, name = 'u_app', what = "REAL")
approx.plotHF(mutar, name = 'u_HF', what = "REAL")
approx.plotErr(mutar, name = 'err', what = "REAL")
# approx.plotRes(mutar, name = 'res', what = "REAL")
appErr = approx.normErr(mutar)
solNorm = approx.normHF(mutar)
resNorm = approx.normRes(mutar)
RHSNorm = approx.normRHS(mutar)
print(('SolNorm:\t{}\nErr:\t{}\nErrRel:\t{}').format(solNorm, appErr,
np.divide(appErr, solNorm)))
print(('RHSNorm:\t{}\nRes:\t{}\nResRel:\t{}').format(RHSNorm, resNorm,
np.divide(resNorm, RHSNorm)))
try:
from plot_zero_set import plotZeroSet2
muZeroVals, Qvals = plotZeroSet2(murange, murangeEff, approx, mu0,
200, [1., 1], polesImTol = 2.)
except: pass
if show_norm:
solver._solveBatchSize = 100
from plot_inf_set import plotInfSet2
muInfVals, normEx, normApp, normRes, normErr, beta = plotInfSet2(
murange, murangeEff, approx, mu0, 50,
[1., 1.], relative = False,
nobeta = True)
print(1.j * approx.getPoles([None, 50.]))
diff --git a/examples/2d/pod/square_pod.py b/examples/2d/pod/square_pod.py
index 43fc0a9..990ddf5 100644
--- a/examples/2d/pod/square_pod.py
+++ b/examples/2d/pod/square_pod.py
@@ -1,192 +1,191 @@
import numpy as np
from rrompy.hfengines.linear_problem.bidimensional import \
HelmholtzSquareDomainProblemEngine as HSDPE
from rrompy.reduction_methods.standard import RationalInterpolant as RI
from rrompy.reduction_methods.standard import ReducedBasis as RB
from rrompy.reduction_methods.pole_matching import \
RationalInterpolantPoleMatching as RIPM
from rrompy.reduction_methods.pole_matching import \
ReducedBasisPoleMatching as RBPM
from rrompy.parameter.parameter_sampling import (QuadratureSampler as QS,
QuadratureSamplerTotal as QST,
ManualSampler as MS,
RandomSampler as RS)
verb = 5
size = 7
show_sample = False
show_norm = True
ignore_forcing = True
ignore_forcing = False
clip = -1
#clip = .4
#clip = .6
MN = 6
R = (MN + 2) * (MN + 1) // 2
S = [int(np.ceil(R ** .5))] * 2
SPivot = [MN + 1, 3]
MMarginal = SPivot[1] - 2
matchingWeight = 10.
samples = "centered"
samples = "standard"
samples = "pole matching"
algo = "rational"
#algo = "RB"
sampling = "quadrature"
sampling = "quadrature_total"
#sampling = "random"
samplingM = "quadrature"
#samplingM = "quadrature_total"
#samplingM = "random"
if samples == "pole matching":
radialM = ""
radialM = "_gaussian"
# radialM = "_thinplate"
# radialM = "_multiquadric"
rMW0 = 2.
radialWeightM = [rMW0]
if size == 1: # small
mu0 = [4 ** .5, 1.5 ** .5]
mutar = [5 ** .5, 1.75 ** .5]
murange = [[2 ** .5, 1. ** .5], [6 ** .5, 2. ** .5]]
elif size == 2: # medium
mu0 = [4 ** .5, 1.75 ** .5]
mutar = [5 ** .5, 1.25 ** .5]
murange = [[1 ** .5, 1. ** .5], [7 ** .5, 2.5 ** .5]]
elif size == 3: # fat
mu0 = [6 ** .5, 4 ** .5]
mutar = [2 ** .5, 2.5 ** .5]
murange = [[0 ** .5, 2 ** .5], [12 ** .5, 6 ** .5]]
elif size == 4: # crowded
mu0 = [10 ** .5, 2 ** .5]
mutar = [9 ** .5, 2.25 ** .5]
murange = [[8 ** .5, 1.5 ** .5], [12 ** .5, 2.5 ** .5]]
elif size == 5: # tall
mu0 = [11 ** .5, 2.25 ** .5]
mutar = [10.5 ** .5, 2.5 ** .5]
murange = [[10 ** .5, 1.5 ** .5], [12 ** .5, 3 ** .5]]
elif size == 6: # taller
mu0 = [11 ** .5, 2.25 ** .5]
mutar = [10.5 ** .5, 2.5 ** .5]
murange = [[10 ** .5, 1.25 ** .5], [12 ** .5, 3.25 ** .5]]
elif size == 7: # low
mu0 = [7 ** .5, .75 ** .5]
mutar = [6.5 ** .5, .9 ** .5]
murange = [[6 ** .5, .5 ** .5], [8 ** .5, 1. ** .5]]
aEff = 1.1
bEff = 1. - aEff
murangeEff = [[(aEff*murange[0][0]**2. + bEff*murange[1][0]**2.) ** .5,
(aEff*murange[0][1]**2. + bEff*murange[1][1]**2.) ** .5],
[(aEff*murange[1][0]**2. + bEff*murange[0][0]**2.) ** .5,
(aEff*murange[1][1]**2. + bEff*murange[0][1]**2.) ** .5]]
solver = HSDPE(kappa = 2.5, theta = np.pi / 3, mu0 = mu0, n = 20,
verbosity = verb)
if ignore_forcing: solver.nbs = 1
rescalingExp = [2., 1.]
if algo == "rational":
params = {'N':MN, 'M':MN, 'S':S, 'POD':True}
if samples == "standard":
params['polybasis'] = "CHEBYSHEV"
# params['polybasis'] = "LEGENDRE"
# params['polybasis'] = "MONOMIAL"
- params['E'] = MN
method = RI
elif samples == "centered":
params['polybasis'] = "MONOMIAL"
params['S'] = R
method = RI
elif samples == "pole matching":
params['S'] = [SPivot[0]]
params['SMarginal'] = [SPivot[1]]
params['MMarginal'] = MMarginal
params['polybasisPivot'] = "CHEBYSHEV"
params['polybasisMarginal'] = "MONOMIAL" + radialM
params['matchingWeight'] = matchingWeight
params['radialDirectionalWeightsMarginal'] = radialWeightM
method = RIPM
else: #if algo == "RB":
params = {'R':R, 'S':S, 'POD':True}
if samples == "standard":
method = RB
elif samples == "centered":
params['S'] = R
method = RB
elif samples == "pole matching":
params['S'] = [SPivot[0]]
params['SMarginal'] = [SPivot[1]]
params['MMarginal'] = MMarginal
params['polybasisMarginal'] = "MONOMIAL" + radialM
params['matchingWeight'] = matchingWeight
params['radialDirectionalWeightsMarginal'] = radialWeightM
method = RBPM
if samples == "standard":
if sampling == "quadrature":
params['sampler'] = QS(murange, "CHEBYSHEV", scalingExp = rescalingExp)
# params['sampler'] = QS(murange, "GAUSSLEGENDRE", scalingExp = rescalingExp)
# params['sampler'] = QS(murange, "UNIFORM", scalingExp = rescalingExp)
params['S'] = [max(j, MN + 1) for j in params['S']]
elif sampling == "quadrature_total":
params['sampler'] = QST(murange, "CHEBYSHEV", scalingExp = rescalingExp)
params['S'] = R
else: # if sampling == "random":
params['sampler'] = RS(murange, "HALTON", scalingExp = rescalingExp)
params['S'] = R
elif samples == "centered":
params['sampler'] = MS(murange, points = [mu0], scalingExp = rescalingExp)
elif samples == "pole matching":
if sampling == "quadrature":
params['samplerPivot'] = QS([murange[0][0], murange[1][0]], "CHEBYSHEV")
# params['samplerPivot'] = QS([murange[0][0], murange[1][0]], "GAUSSLEGENDRE")
# params['samplerPivot'] = QS([murange[0][0], murange[1][0]], "UNIFORM")
elif sampling == "quadrature_total":
params['samplerPivot'] = QST([murange[0][0], murange[1][0]], "CHEBYSHEV")
else: # if sampling == "random":
params['samplerPivot'] = RS([murange[0][0], murange[1][0]], "HALTON")
if samplingM == "quadrature":
params['samplerMarginal'] = QS([murange[0][1], murange[1][1]], "UNIFORM")
elif samplingM == "quadrature_total":
params['samplerMarginal'] = QST([murange[0][1], murange[1][1]], "CHEBYSHEV")
else: # if samplingM == "random":
params['samplerMarginal'] = RS([murange[0][1], murange[1][1]], "HALTON")
if samples != "pole matching":
approx = method(solver, mu0 = mu0, approxParameters = params,
verbosity = verb)
else:
approx = method(solver, mu0 = mu0, directionPivot = [0],
approxParameters = params, verbosity = verb)
if samples != "centered": approx.samplingEngine.allowRepeatedSamples = False
approx.setupApprox()
if show_sample:
approx.plotApprox(mutar, name = 'u_app')
approx.plotHF(mutar, name = 'u_HF')
approx.plotErr(mutar, name = 'err')
approx.plotRes(mutar, name = 'res')
appErr, solNorm = approx.normErr(mutar), approx.normHF(mutar)
resNorm, RHSNorm = approx.normRes(mutar), approx.normRHS(mutar)
print(('SolNorm:\t{}\nErr:\t{}\nErrRel:\t{}').format(solNorm, appErr,
np.divide(appErr, solNorm)))
print(('RHSNorm:\t{}\nRes:\t{}\nResRel:\t{}').format(RHSNorm, resNorm,
np.divide(resNorm, RHSNorm)))
if algo == "rational":
from plot_zero_set import plotZeroSet2
muZeroVals, Qvals = plotZeroSet2(murange, murangeEff, approx, mu0,
200, [2., 1.])
if show_norm:
solver._solveBatchSize = 100
from plot_inf_set import plotInfSet2
muInfVals, normEx, normApp, normRes, normErr, beta = plotInfSet2(
murange, murangeEff, approx, mu0, 20, [2., 1.], clip = clip,
relative = False, nobeta = True)
diff --git a/examples/2d/pod/square_simplified_pod.py b/examples/2d/pod/square_simplified_pod.py
index 8273c12..9d3dbef 100644
--- a/examples/2d/pod/square_simplified_pod.py
+++ b/examples/2d/pod/square_simplified_pod.py
@@ -1,127 +1,126 @@
import numpy as np
from rrompy.hfengines.linear_problem.bidimensional import \
HelmholtzSquareSimplifiedDomainProblemEngine as HSSDPE
from rrompy.reduction_methods.standard import RationalInterpolant as RI
from rrompy.reduction_methods.standard import ReducedBasis as RB
from rrompy.parameter.parameter_sampling import (QuadratureSampler as QS,
QuadratureSamplerTotal as QST,
ManualSampler as MS,
RandomSampler as RS)
verb = 5
size = 1
show_sample = False
show_norm = True
MN = 5
R = (MN + 2) * (MN + 1) // 2
S = [int(np.ceil(R ** .5))] * 2
samples = "centered"
samples = "standard"
algo = "rational"
#algo = "RB"
sampling = "quadrature"
sampling = "quadrature_total"
sampling = "random"
if size == 1: # small
mu0 = [4 ** .5, 1.5 ** .5]
mutar = [5 ** .5, 1.75 ** .5]
murange = [[2 ** .5, 1. ** .5], [6 ** .5, 2. ** .5]]
elif size == 2: # medium
mu0 = [4 ** .5, 1.75 ** .5]
mutar = [5 ** .5, 1.25 ** .5]
murange = [[1 ** .5, 1. ** .5], [7 ** .5, 2.5 ** .5]]
elif size == 3: # fat
mu0 = [6 ** .5, 4 ** .5]
mutar = [2 ** .5, 2.5 ** .5]
murange = [[0 ** .5, 2 ** .5], [12 ** .5, 6 ** .5]]
elif size == 4: # crowded
mu0 = [10 ** .5, 2 ** .5]
mutar = [9 ** .5, 2.25 ** .5]
murange = [[8 ** .5, 1.5 ** .5], [12 ** .5, 2.5 ** .5]]
elif size == 5: # tall
mu0 = [11 ** .5, 2.25 ** .5]
mutar = [10.5 ** .5, 2.5 ** .5]
murange = [[10 ** .5, 1.5 ** .5], [12 ** .5, 3 ** .5]]
elif size == 6: # taller
mu0 = [11 ** .5, 2.25 ** .5]
mutar = [10.5 ** .5, 2.5 ** .5]
murange = [[10 ** .5, 1.25 ** .5], [12 ** .5, 3.25 ** .5]]
elif size == 7: # low
mu0 = [7 ** .5, .75 ** .5]
mutar = [8 ** .5, 1 ** .5]
murange = [[6 ** .5, .25 ** .5], [8 ** .5, 1.25 ** .5]]
aEff = 1.25
bEff = 1. - aEff
murangeEff = [[(aEff*murange[0][0]**2. + bEff*murange[1][0]**2.) ** .5,
(aEff*murange[0][1]**2. + bEff*murange[1][1]**2.) ** .5],
[(aEff*murange[1][0]**2. + bEff*murange[0][0]**2.) ** .5,
(aEff*murange[1][1]**2. + bEff*murange[0][1]**2.) ** .5]]
solver = HSSDPE(kappa = 2.5, theta = np.pi / 3, mu0 = mu0, n = 20,
verbosity = verb)
rescalingExp = [2.] * 2
if algo == "rational":
params = {'N':MN, 'M':MN, 'S':S, 'POD':True}
if samples == "standard":
params['polybasis'] = "CHEBYSHEV"
params['polybasis'] = "LEGENDRE"
params['polybasis'] = "MONOMIAL"
- params['E'] = MN
elif samples == "centered":
params['polybasis'] = "MONOMIAL"
params['S'] = R
method = RI
else: #if algo == "RB":
params = {'R':R, 'S':S, 'POD':True}
if samples == "centered":
params['S'] = R
method = RB
if samples == "standard":
if sampling == "quadrature":
params['sampler'] = QS(murange, "CHEBYSHEV", scalingExp = rescalingExp)
# params['sampler'] = QS(murange, "GAUSSLEGENDRE", scalingExp = rescalingExp)
# params['sampler'] = QS(murange, "UNIFORM", scalingExp = rescalingExp)
params['S'] = [max(j, MN + 1) for j in params['S']]
elif sampling == "quadrature_total":
params['sampler'] = QST(murange, "CHEBYSHEV", scalingExp = rescalingExp)
params['S'] = R
else: # if sampling == "random":
params['sampler'] = RS(murange, "HALTON", scalingExp = rescalingExp)
params['S'] = R
elif samples == "centered":
params['sampler'] = MS(murange, points = [mu0], scalingExp = rescalingExp)
approx = method(solver, mu0 = mu0, approxParameters = params,
verbosity = verb)
approx.setupApprox()
if show_sample:
approx.plotApprox(mutar, name = 'u_app')
approx.plotHF(mutar, name = 'u_HF')
approx.plotErr(mutar, name = 'err')
approx.plotRes(mutar, name = 'res')
appErr, solNorm = approx.normErr(mutar), approx.normHF(mutar)
resNorm, RHSNorm = approx.normRes(mutar), approx.normRHS(mutar)
print(('SolNorm:\t{}\nErr:\t{}\nErrRel:\t{}').format(solNorm, appErr,
np.divide(appErr, solNorm)))
print(('RHSNorm:\t{}\nRes:\t{}\nResRel:\t{}').format(RHSNorm, resNorm,
np.divide(resNorm, RHSNorm)))
if algo == "rational":
from plot_zero_set import plotZeroSet2
muZeroVals, Qvals = plotZeroSet2(murange, murangeEff, approx, mu0,
200, [2., 2.])
if show_norm:
solver._solveBatchSize = 100
from plot_inf_set import plotInfSet2
muInfVals, normEx, normApp, normRes, normErr, beta = plotInfSet2(
murange, murangeEff, approx, mu0, 25, [2., 2.])
diff --git a/examples/2d/pod/synthetic_pod.py b/examples/2d/pod/synthetic_pod.py
index 061550d..274ffed 100644
--- a/examples/2d/pod/synthetic_pod.py
+++ b/examples/2d/pod/synthetic_pod.py
@@ -1,138 +1,137 @@
import numpy as np
from rrompy.hfengines.linear_problem.bidimensional import \
SyntheticBivariateEngine as SBE
from rrompy.reduction_methods.standard import RationalInterpolant as RI
from rrompy.reduction_methods.standard import ReducedBasis as RB
from rrompy.parameter.parameter_sampling import (QuadratureSampler as QS,
QuadratureSamplerTotal as QST,
ManualSampler as MS,
RandomSampler as RS)
verb = 5
size = 1
show_sample = True
show_norm = True
clip = -1
#clip = .4
#clip = .6
Delta = 0
MN = 10
R = (MN + 2) * (MN + 1) // 2
S = [int(np.ceil(R ** .5))] * 2
PODTol = 1e-6
samples = "centered"
samples = "standard"
algo = "rational"
#algo = "RB"
sampling = "quadrature"
#sampling = "quadrature_total"
#sampling = "random"
if samples == "standard":
radial = ""
radial = "_gaussian"
# radial = "_thinplate"
# radial = "_multiquadric"
rW0 = 10.
radialWeight = [rW0] * 2
if size == 1: # small
mu0 = [10. ** .5, 15. ** .5]
mutar = [12. ** .5, 14. ** .5]
murange = [[5. ** .5, 10. ** .5], [15 ** .5, 20 ** .5]]
if size == 2: # large
mu0 = [15. ** .5, 17.5 ** .5]
mutar = [18. ** .5, 22. ** .5]
murange = [[5. ** .5, 10. ** .5], [25 ** .5, 25 ** .5]]
if size == 3: # medium
mu0 = [17.5 ** .5, 15 ** .5]
mutar = [20. ** .5, 18. ** .5]
murange = [[10. ** .5, 10. ** .5], [25 ** .5, 20 ** .5]]
assert Delta <= 0
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]]]
kappa = 20. ** .5
theta = - np.pi / 6.
n = 20
L = np.pi
solver = SBE(kappa = kappa, theta = theta, n = n, L = L,
mu0 = mu0, verbosity = verb)
rescalingExp = [2.] * 2
if algo == "rational":
params = {'N':MN, 'M':MN + Delta, 'S':S, 'POD':True}
if samples == "standard":
params['polybasis'] = "CHEBYSHEV" + radial
# params['polybasis'] = "LEGENDRE" + radial
# params['polybasis'] = "MONOMIAL" + radial
- params['E'] = MN
params['radialDirectionalWeights'] = radialWeight
method = RI
elif samples == "centered":
params['polybasis'] = "MONOMIAL"
params['S'] = R
method = RI
else: #if algo == "RB":
params = {'R':(MN + 2 + Delta) * (MN + 1 + Delta) // 2, 'S':S,
'POD':True, 'PODTolerance':PODTol}
if samples == "centered":
params['S'] = R
method = RB
if samples == "standard":
if sampling == "quadrature":
params['sampler'] = QS(murange, "CHEBYSHEV", scalingExp = rescalingExp)
# params['sampler'] = QS(murange, "GAUSSLEGENDRE", scalingExp = rescalingExp)
# params['sampler'] = QS(murange, "UNIFORM", scalingExp = rescalingExp)
params['S'] = [max(j, MN + 1) for j in params['S']]
elif sampling == "quadrature_total":
params['sampler'] = QST(murange, "CHEBYSHEV", scalingExp = rescalingExp)
params['S'] = R
else: # if sampling == "random":
params['sampler'] = RS(murange, "HALTON", scalingExp = rescalingExp)
params['S'] = R
elif samples == "centered":
params['sampler'] = MS(murange, points = [mu0], scalingExp = rescalingExp)
approx = method(solver, mu0 = mu0, approxParameters = params, verbosity = verb)
if samples == "standard": 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)))
if algo == "rational" and approx.N > 0:
from plot_zero_set import plotZeroSet2
muZeroVals, Qvals = plotZeroSet2(murange, murangeEff, approx, mu0,
200, [2., 2.], clip = clip)
if show_norm:
solver._solveBatchSize = 100
from plot_inf_set import plotInfSet2
muInfVals, normEx, normApp, normRes, normErr, beta = plotInfSet2(
murange, murangeEff, approx, mu0, 50,
[2., 2.], clip = clip, relative = False)
diff --git a/examples/3d/pod/fracture3_pod.py b/examples/3d/pod/fracture3_pod.py
index d1a689f..1bc1751 100644
--- a/examples/3d/pod/fracture3_pod.py
+++ b/examples/3d/pod/fracture3_pod.py
@@ -1,245 +1,251 @@
import numpy as np
from mpl_toolkits.mplot3d import Axes3D
from matplotlib import pyplot as plt
from rrompy.hfengines.linear_problem.tridimensional import \
MembraneFractureEngine3 as MFE
from rrompy.reduction_methods.standard import RationalInterpolant as RI
from rrompy.reduction_methods.standard import ReducedBasis as RB
from rrompy.reduction_methods.pole_matching import \
RationalInterpolantPoleMatching as RIPM
from rrompy.reduction_methods.pole_matching import \
ReducedBasisPoleMatching as RBPM
from rrompy.parameter.parameter_sampling import (QuadratureSampler as QS,
QuadratureSamplerTotal as QST,
ManualSampler as MS,
RandomSampler as RS)
verb = 70
size = 4
show_sample = False
show_norm = False
clip = -1
#clip = .4
#clip = .6
Delta = 0
MN = 5
R = (MN + 3) * (MN + 2) * (MN + 1) // 6
S = [int(np.ceil(R ** (1. / 3.)))] * 3
PODTol = 1e-8
-SPivot = [MN + 1, 2, 2]
+SPivot = [MN + 1, 5, 5]
MMarginal = SPivot[1] - 1
matchingWeight = 10.
samples = "centered"
samples = "standard"
-samples = "pole matching"
+#samples = "pole matching"
algo = "rational"
#algo = "RB"
sampling = "quadrature"
-sampling = "quadrature_total"
+#sampling = "quadrature_total"
#sampling = "random"
samplingM = "quadrature"
-#samplingM = "quadrature_total"
+samplingM = "quadrature_total"
#samplingM = "random"
+polydegreetype = "TOTAL"
+#polydegreetype = "FULL"
if samples == "standard":
radial = ""
- radial = "_gaussian"
+# radial = "_gaussian"
# radial = "_thinplate"
# radial = "_multiquadric"
rW0 = 5.
radialWeight = [rW0] * 3
if samples == "pole matching":
radial = ""
# radial = "_gaussian"
# radial = "_thinplate"
# radial = "_multiquadric"
rW0 = 5.
radialWeight = [rW0]
radialM = ""
- radialM = "_gaussian"
+# radialM = "_gaussian"
# radialM = "_thinplate"
# radialM = "_multiquadric"
rMW0 = 2.
radialWeightM = [rMW0] * 2
assert Delta <= 0
if size == 1:
mu0 = [47.5 ** .5, .4, .05]
mutar = [50 ** .5, .45, .07]
murange = [[40 ** .5, .3, .025], [55 ** .5, .5, .075]]
if size == 2:
mu0 = [50 ** .5, .3, .02]
mutar = [55 ** .5, .4, .03]
murange = [[40 ** .5, .1, .005], [60 ** .5, .5, .035]]
if size == 3:
mu0 = [45 ** .5, .5, .05]
mutar = [47 ** .5, .4, .045]
murange = [[40 ** .5, .3, .04], [50 ** .5, .7, .06]]
if size == 4:
mu0 = [45 ** .5, .4, .05]
mutar = [47 ** .5, .45, .045]
murange = [[40 ** .5, .3, .04], [50 ** .5, .5, .06]]
+if size == 5:
+ mu0 = [45 ** .5, .5, .05]
+ mutar = [47 ** .5, .45, .045]
+ murange = [[40 ** .5, .3, .04], [50 ** .5, .7, .06]]
aEff = 1.#25
bEff = 1. - aEff
murangeEff = [[(aEff*murange[0][0]**2. + bEff*murange[1][0]**2.) ** .5,
aEff*murange[0][1] + bEff*murange[1][1],
aEff*murange[0][2] + bEff*murange[1][2]],
[(aEff*murange[1][0]**2. + bEff*murange[0][0]**2.) ** .5,
aEff*murange[1][1] + bEff*murange[0][1],
aEff*murange[1][2] + bEff*murange[0][2]]]
H = 1.
L = .75
delta = .05
n = 50
solver = MFE(mu0 = mu0, H = H, L = L, delta = delta, n = n, verbosity = verb)
rescalingExp = [2., 1., 1.]
if algo == "rational":
- params = {'N':MN, 'M':MN + Delta, 'S':S, 'POD':True}
+ params = {'N':MN, 'M':MN + Delta, 'S':S, 'POD':True,
+ 'polydegreetype':polydegreetype}
if samples == "standard":
params['polybasis'] = "CHEBYSHEV" + radial
# params['polybasis'] = "LEGENDRE" + radial
# params['polybasis'] = "MONOMIAL" + radial
- params['E'] = MN
params['radialDirectionalWeights'] = radialWeight
method = RI
elif samples == "centered":
params['polybasis'] = "MONOMIAL"
params['S'] = R
method = RI
elif samples == "pole matching":
params['S'] = [SPivot[0]]
params['SMarginal'] = SPivot[1:]
params['MMarginal'] = MMarginal
params['polybasisPivot'] = "CHEBYSHEV" + radial
params['polybasisMarginal'] = "MONOMIAL" + radialM
params['matchingWeight'] = matchingWeight
params['radialDirectionalWeightsPivot'] = radialWeight
params['radialDirectionalWeightsMarginal'] = radialWeightM
method = RIPM
else: #if algo == "RB":
params = {'R':(MN + 3 + Delta) * (MN + 2 + Delta) * (MN + 1 + Delta) // 6,
'S':S, 'POD':True, 'PODTolerance':PODTol}
if samples == "standard":
method = RB
elif samples == "centered":
params['S'] = R
method = RB
elif samples == "pole matching":
params['S'] = [SPivot[0]]
params['SMarginal'] = SPivot[1:]
params['MMarginal'] = MMarginal
params['polybasisMarginal'] = "MONOMIAL" + radialM
params['matchingWeight'] = matchingWeight
params['radialDirectionalWeightsMarginal'] = radialWeightM
method = RBPM
if samples == "standard":
if sampling == "quadrature":
params['sampler'] = QS(murange, "CHEBYSHEV", scalingExp = rescalingExp)
# params['sampler'] = QS(murange, "GAUSSLEGENDRE",scalingExp = rescalingExp)
# params['sampler'] = QS(murange, "UNIFORM", scalingExp = rescalingExp)
params['S'] = [max(j, MN + 1) for j in params['S']]
elif sampling == "quadrature_total":
params['sampler'] = QST(murange, "CHEBYSHEV", scalingExp = rescalingExp)
params['S'] = R
else: # if sampling == "random":
params['sampler'] = RS(murange, "HALTON", scalingExp = rescalingExp)
params['S'] = R
elif samples == "centered":
params['sampler'] = MS(murange, points = [mu0], scalingExp = rescalingExp)
elif samples == "pole matching":
if sampling == "quadrature":
params['samplerPivot'] = QS([murange[0][0], murange[1][0]], "CHEBYSHEV")
# params['samplerPivot'] = QS([murange[0][0], murange[1][0]], "GAUSSLEGENDRE")
# params['samplerPivot'] = QS([murange[0][0], murange[1][0]], "UNIFORM")
elif sampling == "quadrature_total":
params['samplerPivot'] = QST([murange[0][0], murange[1][0]], "CHEBYSHEV")
params['S'] = MN + 1
else: # if sampling == "random":
params['samplerPivot'] = RS([murange[0][0], murange[1][0]], "HALTON")
params['S'] = MN + 1
if samplingM == "quadrature":
params['samplerMarginal'] = QS([murange[0][1:], murange[1][1:]], "UNIFORM")
elif samplingM == "quadrature_total":
params['samplerMarginal'] = QST([murange[0][1:], murange[1][1:]], "CHEBYSHEV")
params['SMarginal'] = (MMarginal + 2) * (MMarginal + 1) // 2
params['polybasisMarginal'] = "CHEBYSHEV" + radialM
else: # if samplingM == "random":
params['samplerMarginal'] = RS([murange[0][1:], murange[1][1:]], "HALTON")
params['SMarginal'] = (MMarginal + 2) * (MMarginal + 1) // 2
if samples != "pole matching":
approx = method(solver, mu0 = mu0, approxParameters = params,
verbosity = verb)
else:
approx = method(solver, mu0 = mu0, directionPivot = [0],
approxParameters = params, verbosity = verb)
if samples != "centered": approx.samplingEngine.allowRepeatedSamples = False
approx.setupApprox()
if show_sample:
from fracture3_warping import fracture3_warping
warps = fracture3_warping(solver.V.mesh(), L, mutar[1], delta, mutar[2])
approx.plotApprox(mutar, warps, name = 'u_app', what = "REAL")
approx.plotHF(mutar, warps, name = 'u_HF', what = "REAL")
approx.plotErr(mutar, warps, name = 'err', what = "REAL")
# approx.plotRes(mutar, warps, name = 'res', what = "REAL")
appErr = approx.normErr(mutar)
solNorm = approx.normHF(mutar)
resNorm = approx.normRes(mutar)
RHSNorm = approx.normRHS(mutar)
print(('SolNorm:\t{}\nErr:\t{}\nErrRel:\t{}').format(solNorm, appErr,
np.divide(appErr, solNorm)))
print(('RHSNorm:\t{}\nRes:\t{}\nResRel:\t{}').format(RHSNorm, resNorm,
np.divide(resNorm, RHSNorm)))
fig = plt.figure(figsize = (8, 6))
ax = Axes3D(fig)
ax.scatter(np.real(approx.mus(0) ** 2.), np.real(approx.mus(1)),
np.real(approx.mus(2)), '.')
plt.show()
plt.close()
approx.verbosity = 0
approx.trainedModel.verbosity = 0
from plot_zero_set_3 import plotZeroSet3
#muZeroVals, Qvals = plotZeroSet3(murange, murangeEff, approx, mu0, 200,
# [None, mu0[1], mu0[2]], [2., 1., 1.],
# clip = clip)
#muZeroVals, Qvals = plotZeroSet3(murange, murangeEff, approx, mu0, 200,
# [None, None, mu0[2]], [2., 1., 1.],
# clip = clip)
#muZeroVals, Qvals = plotZeroSet3(murange, murangeEff, approx, mu0, 200,
# [None, mu0[1], None], [2., 1., 1.],
# clip = clip)
muZeroScatter = plotZeroSet3(murange, murangeEff, approx, mu0, 50,
[None, None, None], [2., 1., 1.], clip = clip)
if show_norm:
solver._solveBatchSize = 25
from plot_inf_set_3 import plotInfSet3
muInfVals, normEx, normApp, normRes, normErr, beta = plotInfSet3(
murange, murangeEff, approx, mu0, 25,
[None, mu0[1], mu0[2]], [2., 1., 1.],
clip = clip, relative = False,
normalizeDen = True)
muInfVals, normEx, normApp, normRes, normErr, beta = plotInfSet3(
murange, murangeEff, approx, mu0, 25,
[None, None, mu0[2]], [2., 1., 1.],
clip = clip, relative = False,
normalizeDen = True)
muInfVals, normEx, normApp, normRes, normErr, beta = plotInfSet3(
murange, murangeEff, approx, mu0, 25,
[None, mu0[1], None], [2., 1., 1.],
clip = clip, relative = False,
normalizeDen = True)
print(approx.getPoles([None, .5, .05]) ** 2.)
diff --git a/examples/3d/pod/matrix_3_pod.py b/examples/3d/pod/matrix_3_pod.py
index f901cba..6561f79 100644
--- a/examples/3d/pod/matrix_3_pod.py
+++ b/examples/3d/pod/matrix_3_pod.py
@@ -1,186 +1,185 @@
import numpy as np
import scipy.sparse as sp
from mpl_toolkits.mplot3d import Axes3D
from matplotlib import pyplot as plt
from rrompy.hfengines.base import MatrixEngineBase as MEB
from rrompy.reduction_methods.standard import RationalInterpolant as RI
from rrompy.reduction_methods.standard import ReducedBasis as RB
from rrompy.parameter.parameter_sampling import (QuadratureSampler as QS,
QuadratureSamplerTotal as QST,
ManualSampler as MS,
RandomSampler as RS)
verb = 50
deg = 2
size = 1
show_sample = True
show_norm = False
Delta = 0
MN = 15
R = (MN + 3) * (MN + 2) * (MN + 1) // 6
S = [int(np.ceil(R ** (1. / 3.)))] * 3
PODTol = 1e-8
samples = "centered"
samples = "standard"
#samples, nDer = "standard_MMM", 2
algo = "rational"
#algo = "RB"
sampling = "quadrature"
sampling = "quadrature_total"
sampling = "random"
if samples == "standard":
radial = ""
# radial = "_gaussian"
# radial = "_thinplate"
# radial = "_multiquadric"
rW0 = 5.
radialWeight = [rW0] * 2
assert Delta <= 0
if size == 1:
mu0 = [0.] * 3
mutar = [.8, .8, .8]
murange = [[-1.] * 3, [1.] * 3]
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[0][2] + bEff*murange[1][2]],
[aEff*murange[1][0] + bEff*murange[0][0],
aEff*murange[1][1] + bEff*murange[0][1],
aEff*murange[1][2] + bEff*murange[0][2]]]
N = 150
exp = 1.05
assert exp > 1.
empty = sp.csr_matrix((np.zeros(0), np.zeros(0), np.zeros(N + 1)),
shape = (N, N))
solver = MEB(verbosity = verb)
solver.npar = 3
A0 = sp.spdiags([np.arange(1, 1 + 2 * N, 2) ** exp - (2 * (N // 4)) ** exp],
[0], N, N)
if deg == 1:
solver.nAs = 4
solver.As = [A0, sp.eye(N), - sp.eye(N), - sp.eye(N)]
elif deg == 2:
solver.nAs = 7
solver.As = [A0] + [empty] + [- sp.eye(N)] + [empty] * 3 + [sp.eye(N)]
elif deg == 3:
solver.nAs = 13
solver.As = [A0] + [empty] + [- sp.eye(N)] + [empty] * 9 + [sp.eye(N)]
np.random.seed(420)
solver.nbs = 1
solver.bs = [np.random.randn(N) + 1.j * np.random.randn(N)]
if algo == "rational":
params = {'N':MN, 'M':MN + Delta, 'S':S, 'POD':True}
if samples == "standard":
params['polybasis'] = "CHEBYSHEV" + radial
# params['polybasis'] = "LEGENDRE" + radial
# params['polybasis'] = "MONOMIAL" + radial
- params['E'] = MN
params['radialDirectionalWeights'] = radialWeight
elif samples == "centered":
params['polybasis'] = "MONOMIAL"
params['S'] = R
else: #MMM
params['S'] = nDer * int(np.ceil(R / nDer))
method = RI
else: #if algo == "RB":
params = {'R':(MN + 3 + Delta) * (MN + 2 + Delta) * (MN + 1 + Delta) // 6,
'S':S, 'POD':True, 'PODTolerance':PODTol}
if samples == "standard":
pass
elif samples == "centered":
params['S'] = R
else: #MMM
params['S'] = nDer * int(np.ceil(R / nDer))
method = RB
if samples == "standard":
if sampling == "quadrature":
params['sampler'] = QS(murange, "CHEBYSHEV")
# params['sampler'] = QS(murange, "GAUSSLEGENDRE")
# params['sampler'] = QS(murange, "UNIFORM")
params['S'] = [max(j, MN + 1) for j in params['S']]
elif sampling == "quadrature_total":
params['sampler'] = QST(murange, "CHEBYSHEV")
params['S'] = R
else: # if sampling == "random":
params['sampler'] = RS(murange, "HALTON")
params['S'] = R
elif samples == "centered":
params['sampler'] = MS(murange, points = [mu0])
elif samples == "standard_MMM":
if sampling == "quadrature":
SdirEff = int(np.ceil(int(np.ceil(R / nDer)) ** (1. / 3)))
pts = QS(murange, "CHEBYSHEV").generatePoints([SdirEff] * 3)
elif sampling == "quadrature_total":
pts = QST(murange, "CHEBYSHEV").generatePoints(int(np.ceil(R / nDer)))
else: # if sampling == "random":
pts = RS(murange, "HALTON").generatePoints(int(np.ceil(R / nDer)))
params['sampler'] = MS(murange, points = pts)
approx = method(solver, mu0 = mu0, approxParameters = params, verbosity = verb)
if samples == "standard": approx.samplingEngine.allowRepeatedSamples = False
approx.setupApprox()
if show_sample:
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)))
fig = plt.figure(figsize = (8, 6))
ax = Axes3D(fig)
ax.scatter(approx.trainedModel.data.mus(0), approx.trainedModel.data.mus(1),
approx.trainedModel.data.mus(2), '.')
ax.set_xlim3d(murangeEff[0][0], murangeEff[1][0])
ax.set_ylim3d(murangeEff[0][1], murangeEff[1][1])
ax.set_zlim3d(murangeEff[0][2], murangeEff[1][2])
plt.show()
plt.close()
approx.verbosity = 0
approx.trainedModel.verbosity = 0
if algo == "rational" and approx.N > 0:
from plot_zero_set_3 import plotZeroSet3
# muZeroVals, Qvals = plotZeroSet3(murange, murangeEff, approx, mu0, 200,
# [None, mu0[1], mu0[2]], [1., 1., 1.])
# muZeroVals, Qvals = plotZeroSet3(murange, murangeEff, approx, mu0, 200,
# [None, None, mu0[2]], [1., 1., 1.])
# muZeroVals, Qvals = plotZeroSet3(murange, murangeEff, approx, mu0, 200,
# [None, mu0[1], None], [1., 1., 1.])
plotZeroSet3(murange, murangeEff, approx, mu0, 25, [None, None, None],
[1., 1., 1.], imagTol = 1e-2)
if show_norm:
solver._solveBatchSize = 25
from plot_inf_set_3 import plotInfSet3
muInfVals, normEx, normApp, normRes, normErr, beta = plotInfSet3(
murange, murangeEff, approx, mu0, 25,
[None, mu0[1], mu0[2]], [1., 1., 1.],
relative = False, normalizeDen = True)
muInfVals, normEx, normApp, normRes, normErr, beta = plotInfSet3(
murange, murangeEff, approx, mu0, 25,
[None, None, mu0[2]], [1., 1., 1.],
relative = False, normalizeDen = True)
muInfVals, normEx, normApp, normRes, normErr, beta = plotInfSet3(
murange, murangeEff, approx, mu0, 25,
[None, mu0[1], None], [1., 1., 1.],
relative = False, normalizeDen = True)
diff --git a/examples/3d/pod/scattering1_pod.py b/examples/3d/pod/scattering1_pod.py
index 17f6574..4500cd8 100644
--- a/examples/3d/pod/scattering1_pod.py
+++ b/examples/3d/pod/scattering1_pod.py
@@ -1,185 +1,184 @@
import numpy as np
from mpl_toolkits.mplot3d import Axes3D
from matplotlib import pyplot as plt
from rrompy.hfengines.linear_problem.tridimensional import Scattering1d as S1D
from rrompy.reduction_methods.standard import RationalInterpolant as RI
from rrompy.reduction_methods.standard import ReducedBasis as RB
from rrompy.parameter.parameter_sampling import (QuadratureSampler as QS,
QuadratureSamplerTotal as QST,
ManualSampler as MS,
RandomSampler as RS)
verb = 50
size = 2
show_sample = True
show_norm = False
clip = -1
#clip = .4
#clip = .6
homogeneize = False
#homogeneize = True
Delta = 0
MN = 15
R = (MN + 3) * (MN + 2) * (MN + 1) // 6
S = [int(np.ceil(R ** (1. / 3.)))] * 3
PODTol = 1e-8
samples = "centered"
samples = "standard"
samples, nDer = "standard_MMM", 15
algo = "rational"
#algo = "RB"
sampling = "quadrature"
sampling = "quadrature_total"
sampling = "random"
if samples == "standard":
radial = ""
# radial = "_gaussian"
# radial = "_thinplate"
# radial = "_multiquadric"
rW0 = 5.
radialWeight = [rW0] * 2
assert Delta <= 0
if size == 1:
mu0 = [4., np.pi, 0.]
mutar = [4.05, .95 * np.pi, .2]
murange = [[2., .9 * np.pi, -.5], [6., 1.1 * np.pi, .5]]
if size == 2:
mu0 = [4., np.pi, .375]
mutar = [4.05, .95 * np.pi, .2]
murange = [[2., .9 * np.pi, 0.], [6., 1.1 * np.pi, .75]]
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[0][2] + bEff*murange[1][2]],
[aEff*murange[1][0] + bEff*murange[0][0],
aEff*murange[1][1] + bEff*murange[0][1],
aEff*murange[1][2] + bEff*murange[0][2]]]
n = 100
solver = S1D(mu0 = mu0, n = n, verbosity = verb, homogeneized = homogeneize)
if algo == "rational":
params = {'N':MN, 'M':MN + Delta, 'S':S, 'POD':True}
if samples == "standard":
params['polybasis'] = "CHEBYSHEV" + radial
# params['polybasis'] = "LEGENDRE" + radial
# params['polybasis'] = "MONOMIAL" + radial
- params['E'] = MN
params['radialDirectionalWeights'] = radialWeight
elif samples == "centered":
params['polybasis'] = "MONOMIAL"
params['S'] = R
else: #MMM
params['S'] = nDer * int(np.ceil(R / nDer))
method = RI
else: #if algo == "RB":
params = {'R':(MN + 3 + Delta) * (MN + 2 + Delta) * (MN + 1 + Delta) // 6,
'S':S, 'POD':True, 'PODTolerance':PODTol}
if samples == "standard":
pass
elif samples == "centered":
params['S'] = R
else: #MMM
params['S'] = nDer * int(np.ceil(R / nDer))
method = RB
if samples == "standard":
if sampling == "quadrature":
params['sampler'] = QS(murange, "CHEBYSHEV")
# params['sampler'] = QS(murange, "GAUSSLEGENDRE")
# params['sampler'] = QS(murange, "UNIFORM")
params['S'] = [max(j, MN + 1) for j in params['S']]
elif sampling == "quadrature_total":
params['sampler'] = QST(murange, "CHEBYSHEV")
params['S'] = R
else: # if sampling == "random":
params['sampler'] = RS(murange, "HALTON")
params['S'] = R
elif samples == "centered":
params['sampler'] = MS(murange, points = [mu0])
elif samples == "standard_MMM":
if sampling == "quadrature":
SdirEff = int(np.ceil(int(np.ceil(R / nDer)) ** (1. / 3)))
pts = QS(murange, "CHEBYSHEV").generatePoints([SdirEff] * 3)
elif sampling == "quadrature_total":
pts = QST(murange, "CHEBYSHEV").generatePoints(int(np.ceil(R / nDer)))
else: # if sampling == "random":
pts = RS(murange, "HALTON").generatePoints(int(np.ceil(R / nDer)))
params['sampler'] = MS(murange, points = pts)
approx = method(solver, mu0 = mu0, approxParameters = params,
verbosity = verb)
if samples == "standard": approx.samplingEngine.allowRepeatedSamples = False
approx.setupApprox()
if show_sample:
import fenics as fen
x = fen.SpatialCoordinate(solver.V.mesh())
warps = [x * mutar[1] / solver._L - x, x * solver._L / mutar[1] - x]
approx.plotApprox(mutar, warps, name = 'u_app', what = "REAL")
approx.plotHF(mutar, warps, name = 'u_HF', what = "REAL")
approx.plotErr(mutar, warps, name = 'err', what = "REAL")
# approx.plotRes(mutar, warps, name = 'res', what = "REAL")
appErr = approx.normErr(mutar)
solNorm = approx.normHF(mutar)
resNorm = approx.normRes(mutar)
RHSNorm = approx.normRHS(mutar)
print(('SolNorm:\t{}\nErr:\t{}\nErrRel:\t{}').format(solNorm, appErr,
np.divide(appErr, solNorm)))
print(('RHSNorm:\t{}\nRes:\t{}\nResRel:\t{}').format(RHSNorm, resNorm,
np.divide(resNorm, RHSNorm)))
fig = plt.figure(figsize = (8, 6))
ax = Axes3D(fig)
ax.scatter(approx.trainedModel.data.mus(0), approx.trainedModel.data.mus(1),
approx.trainedModel.data.mus(2), '.')
ax.set_xlim3d(murangeEff[0][0], murangeEff[1][0])
ax.set_ylim3d(murangeEff[0][1], murangeEff[1][1])
ax.set_zlim3d(murangeEff[0][2], murangeEff[1][2])
plt.show()
plt.close()
approx.verbosity = 0
approx.trainedModel.verbosity = 0
if algo == "rational" and approx.N > 0:
from plot_zero_set_3 import plotZeroSet3
muZeroVals, Qvals = plotZeroSet3(murange, murangeEff, approx, mu0, 200,
[None, mu0[1], mu0[2]], [1., 1., 1.],
clip = clip)
muZeroVals, Qvals = plotZeroSet3(murange, murangeEff, approx, mu0, 200,
[None, None, mu0[2]], [1., 1., 1.],
clip = clip)
muZeroVals, Qvals = plotZeroSet3(murange, murangeEff, approx, mu0, 200,
[None, mu0[1], None], [1., 1., 1.],
clip = clip)
plotZeroSet3(murange, murangeEff, approx, mu0, 100, [None, None, None],
[1., 1., 1.], clip = clip, imagTol = 1e-2)
if show_norm:
solver._solveBatchSize = 25
from plot_inf_set_3 import plotInfSet3
muInfVals, normEx, normApp, normRes, normErr, beta = plotInfSet3(
murange, murangeEff, approx, mu0, 25,
[None, mu0[1], mu0[2]], [1., 1., 1.],
clip = clip, relative = False,
normalizeDen = True)
muInfVals, normEx, normApp, normRes, normErr, beta = plotInfSet3(
murange, murangeEff, approx, mu0, 25,
[None, None, mu0[2]], [1., 1., 1.],
clip = clip, relative = False,
normalizeDen = True)
muInfVals, normEx, normApp, normRes, normErr, beta = plotInfSet3(
murange, murangeEff, approx, mu0, 25,
[None, mu0[1], None], [1., 1., 1.],
clip = clip, relative = False,
normalizeDen = True)
print(approx.getPoles([None, np.pi, 0.]))
diff --git a/rrompy/parameter/parameter_sampling/quadrature_sampler.py b/rrompy/parameter/parameter_sampling/quadrature_sampler.py
index 1ea9070..0efb5f6 100644
--- a/rrompy/parameter/parameter_sampling/quadrature_sampler.py
+++ b/rrompy/parameter/parameter_sampling/quadrature_sampler.py
@@ -1,81 +1,82 @@
# 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 .generic_sampler import GenericSampler
from rrompy.utilities.base.types import List, paramList
from rrompy.utilities.exception_manager import RROMPyException
from rrompy.utilities.numerical import lowDiscrepancy, kroneckerer
from rrompy.parameter import checkParameterList
__all__ = ['QuadratureSampler']
class QuadratureSampler(GenericSampler):
"""Generator of quadrature sample points."""
allowedKinds = ["UNIFORM", "CHEBYSHEV", "GAUSSLEGENDRE"]
def __init__(self, lims:paramList, kind : str = "UNIFORM",
scalingExp : List[float] = None):
super().__init__(lims = lims, scalingExp = scalingExp)
self.kind = kind
def __str__(self) -> str:
return "{}_{}".format(super().__str__(), self.kind)
def __repr__(self) -> str:
return self.__str__() + " at " + hex(id(self))
@property
def kind(self):
"""Value of kind."""
return self._kind
@kind.setter
def kind(self, kind):
if kind.upper() not in self.allowedKinds:
raise RROMPyException("Generator kind not recognized.")
self._kind = kind.upper()
def generatePoints(self, n:List[int], reorder : bool = True) -> paramList:
"""Array of sample points."""
if not hasattr(n, "__len__"): n = [n]
super().generatePoints(n)
nleft, nright = 1, np.prod(n)
xmat = np.empty((nright, self.npar), dtype = self.lims.dtype)
for d in range(self.npar):
nright //= n[d]
a = self.lims(0, d) ** self.scalingExp[d]
b = self.lims(1, d) ** self.scalingExp[d]
c, r = (a + b) / 2., (a - b) / 2.
if self.kind == "UNIFORM":
xd = np.linspace(a, b, n[d])
elif self.kind == "CHEBYSHEV":
nodes, _ = np.polynomial.chebyshev.chebgauss(n[d])
xd = c + r * nodes
elif self.kind == "GAUSSLEGENDRE":
nodes, _ = np.polynomial.legendre.leggauss(n[d])
xd = c + r * nodes[::-1]
xd **= 1. / self.scalingExp[d]
- if n[d] > 1 and reorder:
- fejerOrdering = [n[d] - 1] + lowDiscrepancy(n[d] - 1)
- xd = xd[fejerOrdering]
xmat[:, d] = kroneckerer(xd, nleft, nright)
nleft *= n[d]
x = checkParameterList(xmat, self.npar)[0]
+ nright = np.prod(n)
+ if nright > 1 and reorder:
+ fejerOrdering = [nright - 1] + lowDiscrepancy(nright - 1)
+ xmat = xmat[fejerOrdering, :]
return x
diff --git a/rrompy/reduction_methods/greedy/rational_interpolant_greedy.py b/rrompy/reduction_methods/greedy/rational_interpolant_greedy.py
index 2303677..3c27f9b 100644
--- a/rrompy/reduction_methods/greedy/rational_interpolant_greedy.py
+++ b/rrompy/reduction_methods/greedy/rational_interpolant_greedy.py
@@ -1,342 +1,341 @@
# 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
from rrompy.utilities.poly_fitting.polynomial import (polybases, polydomcoeff,
PolynomialInterpolator as PI)
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, Np2D, DictAny, HFEng, paramVal
from rrompy.utilities.base import verbosityManager as vbMng
from rrompy.utilities.exception_manager import RROMPyWarning
from rrompy.utilities.exception_manager import RROMPyException, RROMPyAssert
__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;
- '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 'EXACT', 'BASIC', and 'BARE'; defaults to 'EXACT';
- '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;
- '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.
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 = ["EXACT", "BASIC", "BARE"]
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", "EXACT", -1, 0],
- toBeExcluded = ["E"])
+ ["MONOMIAL", "EXACT", -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 = np.prod(self._S) - 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 prod(S) - 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 'EXACT'."))
errorEstimatorKind = "EXACT"
self._errorEstimatorKind = errorEstimatorKind
self._approxParameters["errorEstimatorKind"] = self.errorEstimatorKind
@property
def nTestPoints(self):
"""Value of nTestPoints."""
return self._nTestPoints
@nTestPoints.setter
def nTestPoints(self, nTestPoints):
if nTestPoints <= 0:
RROMPyWarning("nTestPoints must be at least 1. Overriding to 1.")
nTestPoints = 1
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()
def _errorSamplingRatio(self, mus:Np1D, muRTest:Np1D,
muRTrain:Np1D) -> Np1D:
"""Scalar ratio in explicit error estimator."""
self.setupApprox()
testTile = np.tile(np.reshape(muRTest, (-1, 1)), [1, len(muRTrain)])
nodalVals = np.prod(testTile - np.reshape(muRTrain, (1, -1)), axis = 1)
denVals = self.trainedModel.getQVal(mus)
return np.abs(nodalVals / denVals)
def _RHSNorms(self, radiusb0:Np2D) -> Np1D:
"""High fidelity system RHS norms."""
self.assembleReducedResidualBlocks(full = False)
# 'ij,jk,ik->k', resbb, radiusb0, radiusb0.conj()
b0resb0 = np.sum(self.trainedModel.data.resbb.dot(radiusb0)
* radiusb0.conj(), axis = 0)
RHSnorms = np.power(np.abs(b0resb0), .5)
return RHSnorms
def _errorEstimatorBare(self) -> Np1D:
"""Bare residual-based error estimator."""
self.setupApprox()
self.assembleReducedResidualGramian(self.trainedModel.data.projMat)
pDom = self.trainedModel.data.P.domCoeff
LL = pDom.conj().dot(self.trainedModel.data.gramian.dot(pDom))
Adiag = self.As[0].diagonal()
AnormApprox = np.linalg.norm(Adiag) ** 2. / np.size(Adiag)
return np.abs(AnormApprox * LL) ** .5
def _errorEstimatorBasic(self, muTest:paramVal, ratioTest:complex) -> Np1D:
"""Basic residual-based error estimator."""
resmu = self.HFEngine.residual(self.getApprox(muTest), muTest,
duality = False)[0]
return np.abs(self.estimatorNormEngine.norm(resmu) / ratioTest)
def _errorEstimatorExact(self, muRTrain:Np1D, vanderBase:Np2D) -> Np1D:
"""Exact residual-based error estimator."""
self.setupApprox()
self.assembleReducedResidualBlocks(full = True)
nAsM, nbsM = self.HFEngine.nAs - 1, self.HFEngine.nbs - 1
radiusA = np.zeros((len(self.mus), nAsM, vanderBase.shape[1]),
dtype = np.complex)
radiusb = np.zeros((nbsM, vanderBase.shape[1]), dtype = np.complex)
domcoeff = polydomcoeff(self.trainedModel.data.Q.deg[0],
self.trainedModel.data.Q.polybasis)
Qvals = self.trainedModel.getQVal(self.mus)
for k in range(max(nAsM, nbsM)):
if k < nAsM:
radiusA[:, k :, :] += np.multiply.outer(Qvals * self._fitinv,
vanderBase[: nAsM - k, :])
if k < nbsM:
radiusb[k :, :] += (self._fitinv.dot(Qvals)
* vanderBase[: nbsM - k, :])
Qvals = Qvals * muRTrain
if self.POD:
radiusA = np.tensordot(self.samplingEngine.RPOD, radiusA, 1)
# 'ijkl,klt,ijt->t', resAA, radiusA, radiusA.conj()
LL = np.sum(np.tensordot(self.trainedModel.data.resAA, radiusA, 2)
* radiusA.conj(), axis = (0, 1))
# 'ijk,jkl,il->l', resAb, radiusA, radiusb.conj()
Lf = np.sum(np.tensordot(self.trainedModel.data.resAb[1 :, :, :],
radiusA, 2) * radiusb.conj(), axis = 0)
# 'ij,jk,ik->k', resbb, radiusb, radiusb.conj()
ff = np.sum(self.trainedModel.data.resbb[1 :, 1 :].dot(radiusb)
* radiusb.conj(), axis = 0)
return domcoeff * np.abs(ff - 2. * np.real(Lf) + LL) ** .5
def errorEstimator(self, mus:Np1D) -> Np1D:
"""Standard residual-based error estimator."""
self.setupApprox()
nAsM, nbsM = self.HFEngine.nAs - 1, self.HFEngine.nbs - 1
muRTest = self.centerNormalize(mus)(0)
muRTrain = self.centerNormalize(self.mus)(0)
samplingRatio = self._errorSamplingRatio(mus, muRTest, muRTrain)
if self.errorEstimatorKind == "EXACT":
vanSize = max(nAsM, nbsM)
else:
vanSize = nbsM
vanderBase = np.polynomial.polynomial.polyvander(muRTest, vanSize).T
RHSnorms = self._RHSNorms(vanderBase[: nbsM + 1, :])
if self.errorEstimatorKind == "BARE":
jOpt = self._errorEstimatorBare()
elif self.errorEstimatorKind == "BASIC":
idx_muTestSample = np.argmax(samplingRatio)
jOpt = self._errorEstimatorBasic(mus[idx_muTestSample],
samplingRatio[idx_muTestSample])
else: #if self.errorEstimatorKind == "EXACT":
jOpt = self._errorEstimatorExact(muRTrain, vanderBase[: -1, :])
return jOpt * samplingRatio / RHSnorms
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()
if not hasattr(self, "As") or not hasattr(self, "bs"):
vbMng(self, "INIT", "Computing affine blocks of system.", 10)
self.As = self.HFEngine.affineLinearSystemA(self.mu0,
self.scaleFactor)[1 :]
self.bs = self.HFEngine.affineLinearSystemb(self.mu0,
self.scaleFactor)
vbMng(self, "DEL", "Done computing affine blocks.", 10)
self.greedy(plotEst)
self._S = len(self.mus)
self._N, self._M = [self._S - 1] * 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:
Qf = self._setupDenominator()
Q = Qf[0]
if self.errorEstimatorKind == "EXACT":
self._fitinv = Qf[1].flatten()
else:
Q = PI()
Q.coeffs = np.ones(1, dtype = np.complex)
Q.npar = 1
Q.polybasis = self.polybasis
if self.errorEstimatorKind == "EXACT": self._fitinv = np.ones(1)
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/pivoted/rational_interpolant_pivoted.py b/rrompy/reduction_methods/pivoted/rational_interpolant_pivoted.py
index 65e50ad..aae415a 100644
--- a/rrompy/reduction_methods/pivoted/rational_interpolant_pivoted.py
+++ b/rrompy/reduction_methods/pivoted/rational_interpolant_pivoted.py
@@ -1,594 +1,602 @@
# Copyright (C) 2018 by the RROMPy authors
#
# This file is part of RROMPy.
#
# RROMPy is free software: you can redistribute it and/or modify
# it under the terms of the GNU Lesser General Public License as published by
# the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
#
# RROMPy is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU Lesser General Public License for more details.
#
# You should have received a copy of the GNU Lesser General Public License
# along with RROMPy. If not, see .
#
from copy import deepcopy as copy
import numpy as np
from rrompy.reduction_methods.base import checkRobustTolerance
from .generic_pivoted_approximant import GenericPivotedApproximant
from rrompy.reduction_methods.standard.rational_interpolant import (
RationalInterpolant as RI)
from rrompy.utilities.poly_fitting.polynomial import (polybases as ppb,
polyfitname,
nextDerivativeIndices,
hashDerivativeToIdx as hashD,
hashIdxToDerivative as hashI,
+ polyvander as pvP,
polyvanderTotal as pvTP,
degreeTotalToFull,
+ fullDegreeMaxMask,
+ totalDegreeMaxMask,
PolynomialInterpolator as PI)
from rrompy.utilities.poly_fitting.radial_basis import (polybases as rbpb,
RadialBasisInterpolator as RBI)
from rrompy.utilities.poly_fitting.moving_least_squares import (
polybases as mlspb,
MovingLeastSquaresInterpolator as MLSI)
from rrompy.reduction_methods.trained_model import (TrainedModelPivotedData,
TrainedModelPivotedRational as tModel)
from rrompy.utilities.base.types import (Np1D, Np2D, HFEng, DictAny, Tuple,
List, ListAny, paramVal, paramList)
from rrompy.utilities.base import verbosityManager as vbMng, freepar as fp
from rrompy.utilities.numerical import multifactorial, customPInv
from rrompy.utilities.exception_manager import (RROMPyException, RROMPyAssert,
RROMPyWarning)
from rrompy.parameter import checkParameter
__all__ = ['RationalInterpolantPivoted']
class RationalInterpolantPivoted(GenericPivotedApproximant):
"""
ROM pivoted rational interpolant computation for parametric problems.
Args:
HFEngine: HF problem solver.
mu0(optional): Default parameter. Defaults to 0.
directionPivot(optional): Pivot components. Defaults to [0].
approxParameters(optional): Dictionary containing values for main
parameters of approximant. Recognized keys are:
- 'POD': whether to compute POD of snapshots; defaults to True;
- 'S': total number of pivot samples current approximant relies
upon;
- 'samplerPivot': pivot sample point generator;
- 'SMarginal': total number of marginal samples current approximant
relies upon;
- 'samplerMarginal': marginal sample point generator;
- 'polybasisPivot': type of polynomial basis for pivot
interpolation; defaults to 'MONOMIAL';
- 'polybasisMarginal': type of polynomial basis for marginal
interpolation; defaults to 'MONOMIAL';
- - 'E': number of derivatives used to compute interpolant; defaults
- to 0;
- 'M': degree of rational interpolant numerator; defaults to 0;
- 'N': degree of rational interpolant denominator; defaults to 0;
+ - 'polydegreetype': type of polynomial degree; defaults to 'TOTAL';
- 'MMarginal': degree of marginal interpolant; defaults to 0;
- 'radialDirectionalWeightsPivot': radial basis weights for pivot
numerator; defaults to 0, i.e. identity;
- 'radialDirectionalWeightsMarginal': radial basis weights for
marginal interpolant; defaults to 0, i.e. identity;
- 'interpRcondPivot': tolerance for pivot interpolation; defaults
to None;
- 'interpRcondMarginal': tolerance for marginal interpolation;
defaults to None;
- 'robustTol': tolerance for robust rational denominator
management; defaults to 0.
Defaults to empty dict.
verbosity(optional): Verbosity level. Defaults to 10.
Attributes:
HFEngine: HF problem solver.
mu0: Default parameter.
directionPivot: Pivot components.
mus: Array of snapshot parameters.
musPivot: Array of pivot snapshot parameters.
musMarginal: Array of marginal snapshot parameters.
approxParameters: Dictionary containing values for main parameters of
approximant. Recognized keys are in parameterList.
parameterListSoft: Recognized keys of soft approximant parameters:
- 'POD': whether to compute POD of snapshots;
- 'polybasisPivot': type of polynomial basis for pivot
interpolation;
- 'polybasisMarginal': type of polynomial basis for marginal
interpolation;
- - 'E': number of derivatives used to compute interpolant;
- 'M': degree of rational interpolant numerator;
- 'N': degree of rational interpolant denominator;
+ - 'polydegreetype': type of polynomial degree;
- 'MMarginal': degree of marginal interpolant;
- 'radialDirectionalWeightsPivot': radial basis weights for pivot
numerator;
- 'radialDirectionalWeightsMarginal': radial basis weights for
marginal interpolant;
- 'interpRcondPivot': tolerance for pivot interpolation;
- 'interpRcondMarginal': tolerance for marginal interpolation;
- 'robustTol': tolerance for robust rational denominator
management.
parameterListCritical: Recognized keys of critical approximant
parameters:
- 'S': total number of pivot samples current approximant relies
upon;
- 'samplerPivot': pivot sample point generator;
- 'SMarginal': total number of marginal samples current approximant
relies upon;
- 'samplerMarginal': marginal sample point generator.
POD: Whether to compute POD of snapshots.
S: Total number of pivot samples current approximant relies upon.
sampler: Pivot sample point generator.
polybasisPivot: Type of polynomial basis for pivot interpolation.
polybasisMarginal: Type of polynomial basis for marginal interpolation.
M: Numerator degree of approximant.
N: Denominator degree of approximant.
+ polydegreetype: Type of polynomial degree.
MMarginal: Degree of marginal interpolant.
radialDirectionalWeightsPivot: Radial basis weights for pivot
numerator.
radialDirectionalWeightsMarginal: Radial basis weights for marginal
interpolant.
interpRcondPivot: Tolerance for pivot interpolation.
interpRcondMarginal: Tolerance for marginal interpolation.
robustTol: Tolerance for robust rational denominator management.
- E: Complete derivative depth level of S.
muBoundsPivot: list of bounds for pivot parameter values.
muBoundsMarginal: list of bounds for marginal parameter values.
samplingEngine: Sampling engine.
uHF: High fidelity solution(s) with parameter(s) lastSolvedHF as
sampleList.
lastSolvedHF: Parameter(s) corresponding to last computed high fidelity
solution(s) as parameterList.
uApproxReduced: Reduced approximate solution(s) with parameter(s)
lastSolvedApprox as sampleList.
lastSolvedApproxReduced: Parameter(s) corresponding to last computed
reduced approximate solution(s) as parameterList.
uApprox: Approximate solution(s) with parameter(s) lastSolvedApprox as
sampleList.
lastSolvedApprox: Parameter(s) corresponding to last computed
approximate solution(s) as parameterList.
Q: Numpy 1D vector containing complex coefficients of approximant
denominator.
P: Numpy 2D vector whose columns are FE dofs of coefficients of
approximant numerator.
"""
def __init__(self, HFEngine:HFEng, mu0 : paramVal = None,
directionPivot : ListAny = [0],
approxParameters : DictAny = {}, verbosity : int = 10,
timestamp : bool = True):
self._preInit()
- self._addParametersToList(["polybasisPivot", "E", "M", "N",
+ self._addParametersToList(["polybasisPivot", "M", "N",
+ "polydegreetype",
"radialDirectionalWeightsPivot",
"interpRcondPivot", "robustTol"],
- ["MONOMIAL", -1, 0, 0, 1, -1, 0])
+ ["MONOMIAL", 0, 0, "TOTAL", 1, -1, 0])
super().__init__(HFEngine = HFEngine, mu0 = mu0,
directionPivot = directionPivot,
approxParameters = approxParameters,
verbosity = verbosity, timestamp = timestamp)
self._postInit()
@property
def polybasisPivot(self):
"""Value of polybasisPivot."""
return self._polybasisPivot
@polybasisPivot.setter
def polybasisPivot(self, polybasisPivot):
try:
polybasisPivot = polybasisPivot.upper().strip().replace(" ","")
if polybasisPivot not in ppb + rbpb + mlspb:
raise RROMPyException(
"Prescribed pivot polybasis not recognized.")
self._polybasisPivot = polybasisPivot
except:
RROMPyWarning(("Prescribed pivot polybasis not recognized. "
"Overriding to 'MONOMIAL'."))
self._polybasisPivot = "MONOMIAL"
self._approxParameters["polybasisPivot"] = self.polybasisPivot
@property
def polybasisPivot0(self):
if "_" in self.polybasisPivot:
return self.polybasisPivot.split("_")[0]
return self.polybasisPivot
@property
def radialDirectionalWeightsPivot(self):
"""Value of radialDirectionalWeightsPivot."""
return self._radialDirectionalWeightsPivot
@radialDirectionalWeightsPivot.setter
def radialDirectionalWeightsPivot(self, radialDirectionalWeightsPivot):
self._radialDirectionalWeightsPivot = radialDirectionalWeightsPivot
self._approxParameters["radialDirectionalWeightsPivot"] = (
self.radialDirectionalWeightsPivot)
@property
def interpRcondPivot(self):
"""Value of interpRcondPivot."""
return self._interpRcondPivot
@interpRcondPivot.setter
def interpRcondPivot(self, interpRcondPivot):
self._interpRcondPivot = interpRcondPivot
self._approxParameters["interpRcondPivot"] = self.interpRcondPivot
@property
def M(self):
- """Value of M. Its assignment may change S."""
+ """Value of M."""
return self._M
@M.setter
def M(self, M):
if M < 0: raise RROMPyException("M must be non-negative.")
self._M = M
self._approxParameters["M"] = self.M
- if hasattr(self, "_E") and self.E >= 0 and self.E < self.M:
- RROMPyWarning("Prescribed S is too small. Decreasing M.")
- self.M = self.E
@property
def N(self):
- """Value of N. Its assignment may change S."""
+ """Value of N."""
return self._N
@N.setter
def N(self, N):
if N < 0: raise RROMPyException("N must be non-negative.")
self._N = N
self._approxParameters["N"] = self.N
- if hasattr(self, "_E") and self.E >= 0 and self.E < self.N:
- RROMPyWarning("Prescribed S is too small. Decreasing N.")
- self.N = self.E
@property
- def E(self):
- """Value of E."""
- return self._E
- @E.setter
- def E(self, E):
- if E < 0:
- if not hasattr(self, "_S"):
- raise RROMPyException(("Value of E must be positive if S is "
- "not yet assigned."))
- E = np.sum(hashI(np.prod(self.S), len(self.directionPivot))) - 1
- self._E = E
- self._approxParameters["E"] = self.E
- if (hasattr(self, "_S")
- and self.E >= np.sum(hashI(np.prod(self.S),len(self.directionPivot)))):
- RROMPyWarning("Prescribed S is too small. Decreasing E.")
- self.E = -1
- if hasattr(self, "_M"): self.M = self.M
- if hasattr(self, "_N"): self.N = self.N
+ def polydegreetype(self):
+ """Value of polydegreetype."""
+ return self._polydegreetype
+ @polydegreetype.setter
+ def polydegreetype(self, polydegreetype):
+ try:
+ polydegreetype = polydegreetype.upper().strip().replace(" ","")
+ if polydegreetype not in ["TOTAL", "FULL"]:
+ raise RROMPyException(("Prescribed polydegreetype not "
+ "recognized."))
+ self._polydegreetype = polydegreetype
+ except:
+ RROMPyWarning(("Prescribed polydegreetype not recognized. "
+ "Overriding to 'TOTAL'."))
+ self._polydegreetype = "TOTAL"
+ self._approxParameters["polydegreetype"] = self.polydegreetype
@property
def robustTol(self):
"""Value of tolerance for robust rational denominator management."""
return self._robustTol
@robustTol.setter
def robustTol(self, robustTol):
if robustTol < 0.:
RROMPyWarning(("Overriding prescribed negative robustness "
"tolerance to 0."))
robustTol = 0.
self._robustTol = robustTol
self._approxParameters["robustTol"] = self.robustTol
@property
def S(self):
"""Value of S."""
return self._S
@S.setter
def S(self, S):
GenericPivotedApproximant.S.fset(self, S)
if hasattr(self, "_M"): self.M = self.M
if hasattr(self, "_N"): self.N = self.N
- if hasattr(self, "_E"): self.E = self.E
def resetSamples(self):
"""Reset samples."""
super().resetSamples()
self._musPUniqueCN = None
self._derPIdxs = None
self._reorderP = None
def _setupPivotInterpolationIndices(self):
"""Setup parameters for polyvander."""
RROMPyAssert(self._mode,
message = "Cannot setup interpolation indices.")
if (self._musPUniqueCN is None
or len(self._reorderP) != len(self.musPivot)):
self._musPUniqueCN, musPIdxsTo, musPIdxs, musPCount = (
self.centerNormalizePivot(self.musPivot).unique(
return_index = True,
return_inverse = True,
return_counts = True))
self._musPUnique = self.mus[musPIdxsTo]
self._derPIdxs = [None] * len(self._musPUniqueCN)
self._reorderP = np.empty(len(musPIdxs), dtype = int)
filled = 0
for j, cnt in enumerate(musPCount):
self._derPIdxs[j] = nextDerivativeIndices([],
self.musPivot.shape[1], cnt)
jIdx = np.nonzero(musPIdxs == j)[0]
self._reorderP[jIdx] = np.arange(filled, filled + cnt)
filled += cnt
def _setupDenominator(self):
"""Compute rational denominator."""
RROMPyAssert(self._mode, message = "Cannot setup denominator.")
vbMng(self, "INIT", "Starting computation of denominator.", 7)
NinvD = None
N0 = copy(self.N)
qs = []
self.verbosity -= 10
for j in range(len(self.musMarginal)):
self._N = N0
while self.N > 0:
if NinvD != self.N:
invD, fitinvP = self._computeInterpolantInverseBlocks()
NinvD = self.N
if self.POD:
ev, eV = RI.findeveVGQR(self, self.samplingEngine.RPOD[j],
invD)
else:
ev, eV = RI.findeveVGExplicit(self,
self.samplingEngine.samples[j], invD)
nevBad = checkRobustTolerance(ev, self.robustTol)
if nevBad <= 1: break
if self.catchInstability:
raise RROMPyException(("Instability in denominator "
"computation: eigenproblem is "
"poorly conditioned."))
RROMPyWarning(("Smallest {} eigenvalues below tolerance. "
"Reducing N by 1.").format(nevBad))
self.N = self.N - 1
if self.N <= 0:
self._N = 0
eV = np.ones((1, 1))
q = PI()
q.npar = self.musPivot.shape[1]
q.polybasis = self.polybasisPivot0
- q.coeffs = degreeTotalToFull(tuple([self.N + 1] * q.npar), q.npar,
- eV[:, 0])
+ if self.polydegreetype == "TOTAL":
+ q.coeffs = degreeTotalToFull(tuple([self.N + 1] * q.npar),
+ q.npar, eV[:, 0])
+ else:
+ q.coeffs = eV[:, 0].reshape([self.N + 1] * q.npar)
qs = qs + [copy(q)]
self.verbosity += 10
vbMng(self, "DEL", "Done computing denominator.", 7)
return qs, fitinvP
def _setupNumerator(self):
"""Compute rational numerator."""
RROMPyAssert(self._mode, message = "Cannot setup numerator.")
vbMng(self, "INIT", "Starting computation of numerator.", 7)
Qevaldiag = np.zeros((len(self.musPivot), len(self.musPivot)),
dtype = np.complex)
verb = self.trainedModel.verbosity
self.trainedModel.verbosity = 0
self._setupPivotInterpolationIndices()
M0 = copy(self.M)
tensor_idx = 0
ps = []
for k, muM in enumerate(self.musMarginal):
self._M = M0
idxGlob = 0
for j, derIdxs in enumerate(self._derPIdxs):
mujEff = [fp] * self.npar
for jj, kk in enumerate(self.directionPivot):
mujEff[kk] = self._musPUnique[j, jj]
for jj, kk in enumerate(self.directionMarginal):
mujEff[kk] = muM(0, jj)
mujEff = checkParameter(mujEff, self.npar)
nder = len(derIdxs)
idxLoc = np.arange(len(self.musPivot))[
(self._reorderP >= idxGlob)
* (self._reorderP < idxGlob + nder)]
idxGlob += nder
Qval = [0] * nder
for der in range(nder):
derIdx = hashI(der, self.musPivot.shape[1])
derIdxEff = [0] * self.npar
sclEff = [0] * self.npar
for jj, kk in enumerate(self.directionPivot):
derIdxEff[kk] = derIdx[jj]
sclEff[kk] = self.scaleFactorPivot[jj] ** -1.
Qval[der] = (self.trainedModel.getQVal(mujEff, derIdxEff,
scl = sclEff)
/ multifactorial(derIdx))
for derU, derUIdx in enumerate(derIdxs):
for derQ, derQIdx in enumerate(derIdxs):
diffIdx = [x - y for (x, y) in zip(derUIdx, derQIdx)]
if all([x >= 0 for x in diffIdx]):
diffj = hashD(diffIdx)
Qevaldiag[idxLoc[derU], idxLoc[derQ]] = Qval[diffj]
while self.M >= 0:
if self.polybasisPivot in ppb:
p = PI()
wellCond, msg = p.setupByInterpolation(
- self._musPUniqueCN, Qevaldiag, self.M,
- self.polybasisPivot, self.verbosity >= 5, True,
- {"derIdxs": self._derPIdxs,
- "reorder": self._reorderP,
- "scl": np.power(self.scaleFactorPivot, -1.)},
- {"rcond": self.interpRcondPivot})
+ self._musPUniqueCN, Qevaldiag, self.M,
+ self.polybasisPivot, self.verbosity >= 5,
+ self.polydegreetype == "TOTAL",
+ {"derIdxs": self._derPIdxs,
+ "reorder": self._reorderP,
+ "scl": np.power(self.scaleFactorPivot, -1.)},
+ {"rcond": self.interpRcondPivot})
elif self.polybasisPivot in rbpb:
p = RBI()
wellCond, msg = p.setupByInterpolation(
self._musPUniqueCN, Qevaldiag, self.M,
self.polybasisPivot,
self.radialDirectionalWeightsPivot,
- self.verbosity >= 5, True,
+ self.verbosity >= 5,
+ self.polydegreetype == "TOTAL",
{"derIdxs": self._derPIdxs,
"reorder": self._reorderP,
"scl": np.power(self.scaleFactorPivot, -1.)},
{"rcond": self.interpRcondPivot})
else:# if self.polybasisPivot in mlspb:
p = MLSI()
wellCond, msg = p.setupByInterpolation(
self._musPUniqueCN, Qevaldiag, self.M,
self.polybasisPivot,
self.radialDirectionalWeightsPivot,
- self.verbosity >= 5, True,
+ self.verbosity >= 5,
+ self.polydegreetype == "TOTAL",
{"derIdxs": self._derPIdxs,
"reorder": self._reorderP,
"scl": np.power(self.scaleFactorPivot, -1.)})
vbMng(self, "MAIN", msg, 5)
if wellCond: break
if self.catchInstability:
raise RROMPyException(("Instability in numerator "
"computation: polyfit is "
"poorly conditioned."))
RROMPyWarning(("Polyfit is poorly conditioned. "
"Reducing M by 1."))
self.M -= 1
if self.M < 0:
raise RROMPyException(("Instability in computation of "
"numerator. Aborting."))
tensor_idx_new = tensor_idx + Qevaldiag.shape[1]
if self.POD:
p.postmultiplyTensorize(self.samplingEngine.RPODCoalesced.T[
tensor_idx : tensor_idx_new, :])
else:
p.pad(tensor_idx, len(self.mus) - tensor_idx_new)
tensor_idx = tensor_idx_new
ps = ps + [copy(p)]
self.trainedModel.verbosity = verb
vbMng(self, "DEL", "Done computing numerator.", 7)
return ps
def setupApprox(self):
"""
Compute rational interpolant.
SVD-based robust eigenvalue management.
"""
if self.checkComputedApprox():
return
RROMPyAssert(self._mode, message = "Cannot setup approximant.")
vbMng(self, "INIT", "Setting up {}.". format(self.name()), 5)
self.computeScaleFactor()
self.computeSnapshots()
if self.trainedModel is None:
self.trainedModel = tModel()
self.trainedModel.verbosity = self.verbosity
self.trainedModel.timestamp = self.timestamp
data = TrainedModelPivotedData(self.trainedModel.name(), self.mu0,
self.samplingEngine.samplesCoalesced,
self.scaleFactor,
self.HFEngine.rescalingExp,
self.directionPivot)
data.musPivot = copy(self.musPivot)
data.musMarginal = copy(self.musMarginal)
self.trainedModel.data = data
else:
self.trainedModel = self.trainedModel
- self.trainedModel.data.projMat = copy(self.samplingEngine.samples)
+ self.trainedModel.data.projMat = copy(
+ self.samplingEngine.samplesCoalesced)
self.trainedModel.data.marginalInterp = self._setupMarginalInterp()
if self.N > 0:
Qs = self._setupDenominator()[0]
else:
Q = PI()
Q.npar = self.musPivot.shape[1]
Q.coeffs = np.ones(tuple([1] * Q.npar), dtype = np.complex)
Q.polybasis = self.polybasisPivot0
Qs = [Q for _ in range(len(self.musMarginal))]
self.trainedModel.data.Qs = Qs
self.trainedModel.data.Ps = self._setupNumerator()
self.trainedModel.data.approxParameters = copy(self.approxParameters)
vbMng(self, "DEL", "Done setting up approximant.", 5)
def _computeInterpolantInverseBlocks(self) -> Tuple[List[Np2D], Np2D]:
"""
Compute inverse factors for minimal interpolant target functional.
"""
RROMPyAssert(self._mode, message = "Cannot solve eigenvalue problem.")
self._setupPivotInterpolationIndices()
- while self.E >= 0:
- eWidth = (hashD([self.E + 1] + [0] * (self.musPivot.shape[1] - 1))
- - hashD([self.E] + [0] * (self.musPivot.shape[1] - 1)))
- TE, _, argIdxs = pvTP(self._musPUniqueCN, self.E,
- self.polybasisPivot0, self._derPIdxs,
- self._reorderP,
- scl = np.power(self.scaleFactorPivot, -1.))
- fitOut = customPInv(TE[:, argIdxs], rcond = self.interpRcondPivot,
+ nparPivot = self.musPivot.shape[1]
+ while self.N >= 0:
+ if self.polydegreetype == "TOTAL":
+ TE, _, argIdxs = pvTP(self._musPUniqueCN, self.N,
+ self.polybasisPivot0, self._derPIdxs,
+ self._reorderP,
+ scl = np.power(self.scaleFactorPivot, -1.))
+ TE = TE[:, argIdxs]
+ idxsB = totalDegreeMaxMask(self.N, nparPivot)
+ else: #if self.polydegreetype == "FULL":
+ TE = pvP(self._musPUniqueCN, [self.N] * nparPivot,
+ self.polybasisPivot0, self._derPIdxs, self._reorderP,
+ scl = np.power(self.scaleFactorPivot, -1.))
+ idxsB = fullDegreeMaxMask(self.N, nparPivot)
+ fitOut = customPInv(TE, rcond = self.interpRcondPivot,
full = True)
vbMng(self, "MAIN",
("Fitting {} samples with degree {} through {}... "
"Conditioning of pseudoinverse system: {:.4e}.").format(
- TE.shape[0], self.E,
+ TE.shape[0], self.N,
polyfitname(self.polybasisPivot0),
fitOut[1][1][0] / fitOut[1][1][-1]),
5)
- if fitOut[1][0] == len(argIdxs):
- fitinvP = fitOut[0][- eWidth : , :]
+ if fitOut[1][0] == TE.shape[1]:
+ fitinvP = fitOut[0][idxsB, :]
break
- RROMPyWarning("Polyfit is poorly conditioned. Reducing E by 1.")
- self.E -= 1
- if self.E < 0:
+ RROMPyWarning("Polyfit is poorly conditioned. Reducing N by 1.")
+ self.N -= 1
+ if self.N < 0:
raise RROMPyException(("Instability in computation of "
"denominator. Aborting."))
TN, _, argIdxs = pvTP(self._musPUniqueCN, self.N, self.polybasisPivot0,
self._derPIdxs, self._reorderP,
scl = np.power(self.scaleFactorPivot, -1.))
TN = TN[:, argIdxs]
- invD = [None] * (eWidth)
- for k in range(eWidth):
+ invD = [None] * (len(idxsB))
+ for k in range(len(idxsB)):
pseudoInv = np.diag(fitinvP[k, :])
idxGlob = 0
for j, derIdxs in enumerate(self._derPIdxs):
nder = len(derIdxs)
idxGlob += nder
if nder > 1:
idxLoc = np.arange(len(self.musPivot))[
(self._reorderP >= idxGlob - nder)
* (self._reorderP < idxGlob)]
invLoc = fitinvP[k, idxLoc]
pseudoInv[np.ix_(idxLoc, idxLoc)] = 0.
for diffj, diffjIdx in enumerate(derIdxs):
for derQ, derQIdx in enumerate(derIdxs):
derUIdx = [x - y for (x, y) in
zip(diffjIdx, derQIdx)]
if all([x >= 0 for x in derUIdx]):
derU = hashD(derUIdx)
pseudoInv[idxLoc[derU], idxLoc[derQ]] = (
invLoc[diffj])
invD[k] = pseudoInv.dot(TN)
return invD, fitinvP
def centerNormalize(self, mu : paramList = [],
mu0 : paramVal = None) -> paramList:
"""
Compute normalized parameter to be plugged into approximant.
Args:
mu: Parameter(s) 1.
mu0: Parameter(s) 2. If None, set to self.mu0.
Returns:
Normalized parameter.
"""
return self.trainedModel.centerNormalize(mu, mu0)
def centerNormalizePivot(self, mu : paramList = [],
mu0 : paramVal = None) -> paramList:
"""
Compute normalized parameter to be plugged into approximant.
Args:
mu: Parameter(s) 1.
mu0: Parameter(s) 2. If None, set to self.mu0.
Returns:
Normalized parameter.
"""
return self.trainedModel.centerNormalizePivot(mu, mu0)
def getResidues(self, *args, **kwargs) -> Np1D:
"""
Obtain approximant residues.
Returns:
Matrix with residues as columns.
"""
return self.trainedModel.getResidues(*args, **kwargs)
diff --git a/rrompy/reduction_methods/pivoted/reduced_basis_pivoted.py b/rrompy/reduction_methods/pivoted/reduced_basis_pivoted.py
index 996f942..33427ba 100644
--- a/rrompy/reduction_methods/pivoted/reduced_basis_pivoted.py
+++ b/rrompy/reduction_methods/pivoted/reduced_basis_pivoted.py
@@ -1,275 +1,273 @@
# Copyright (C) 2018 by the RROMPy authors
#
# This file is part of RROMPy.
#
# RROMPy is free software: you can redistribute it and/or modify
# it under the terms of the GNU Lesser General Public License as published by
# the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
#
# RROMPy is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU Lesser General Public License for more details.
#
# You should have received a copy of the GNU Lesser General Public License
# along with RROMPy. If not, see .
#
from copy import deepcopy as copy
import numpy as np
from .generic_pivoted_approximant import GenericPivotedApproximant
from rrompy.hfengines.base import MarginalProxyEngine
from rrompy.reduction_methods.trained_model import (TrainedModelPivotedData,
TrainedModelPivotedReducedBasis as tModel)
from rrompy.reduction_methods.base.reduced_basis_utils import \
projectAffineDecomposition
from rrompy.utilities.base.types import (Np1D, Np2D, List, ListAny, Tuple,
DictAny, HFEng, paramVal, sampList)
from rrompy.utilities.base import verbosityManager as vbMng, freepar as fp
from rrompy.utilities.numerical import customPInv
from rrompy.utilities.exception_manager import (RROMPyWarning, RROMPyException,
RROMPyAssert)
__all__ = ['ReducedBasisPivoted']
class ReducedBasisPivoted(GenericPivotedApproximant):
"""
ROM pivoted RB approximant computation for parametric problems.
Args:
HFEngine: HF problem solver.
mu0(optional): Default parameter. Defaults to 0.
directionPivot(optional): Pivot components. Defaults to [0].
approxParameters(optional): Dictionary containing values for main
parameters of approximant. Recognized keys are:
- 'POD': whether to compute POD of snapshots; defaults to True;
- 'S': total number of pivot samples current approximant relies
upon;
- 'samplerPivot': pivot sample point generator;
- 'SMarginal': total number of marginal samples current approximant
relies upon;
- 'samplerMarginal': marginal sample point generator;
- 'R': rank for pivot Galerkin projection; defaults to prod(S);
- 'PODTolerance': tolerance for pivot snapshots POD; defaults to
-1;
- 'polybasisMarginal': type of polynomial basis for marginal
interpolation; allowed values include 'MONOMIAL', 'CHEBYSHEV'
and 'LEGENDRE'; defaults to 'MONOMIAL';
- 'MMarginal': degree of marginal interpolant; defaults to 0;
- 'radialDirectionalWeightsMarginal': radial basis weights for
marginal interpolant; defaults to 0, i.e. identity;
- 'interpRcondMarginal': tolerance for marginal interpolation;
defaults to None.
Defaults to empty dict.
verbosity(optional): Verbosity level. Defaults to 10.
Attributes:
HFEngine: HF problem solver.
mu0: Default parameter.
directionPivot: Pivot components.
mus: Array of snapshot parameters.
approxRadius: Dummy radius of approximant (i.e. distance from mu0 to
farthest sample point).
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;
- 'R': rank for Galerkin projection;
- 'PODTolerance': tolerance for snapshots POD;
- 'polybasisMarginal': type of polynomial basis for marginal
interpolation;
- 'MMarginal': degree of marginal interpolant;
- 'radialDirectionalWeightsMarginal': radial basis weights for
marginal interpolant;
- 'interpRcondMarginal': tolerance for marginal interpolation.
parameterListCritical: Recognized keys of critical approximant
parameters:
- 'S': total number of pivot samples current approximant relies
upon;
- 'samplerPivot': pivot sample point generator;
- 'SMarginal': total number of marginal samples current approximant
relies upon;
- 'samplerMarginal': marginal sample point generator.
POD: Whether to compute POD of snapshots.
S: Total number of pivot samples current approximant relies upon.
samplerPivot: Pivot sample point generator.
SMarginal: Total number of marginal samples current approximant relies
upon.
samplerMarginal: Marginal sample point generator.
R: Rank for Galerkin projection.
PODTolerance: Tolerance for pivot snapshots POD.
polybasisMarginal: Type of polynomial basis for marginal interpolation.
MMarginal: Degree of marginal interpolant.
radialDirectionalWeightsMarginal: Radial basis weights for marginal
interpolant.
interpRcondMarginal: Tolerance for marginal interpolation.
muBoundsPivot: list of bounds for pivot parameter values.
muBoundsMarginal: list of bounds for marginal parameter values.
samplingEngine: Sampling engine.
uHF: High fidelity solution(s) with parameter(s) lastSolvedHF as
sampleList.
lastSolvedHF: Parameter(s) corresponding to last computed high fidelity
solution(s) as parameterList.
uApproxReduced: Reduced approximate solution(s) with parameter(s)
lastSolvedApprox as sampleList.
lastSolvedApproxReduced: Parameter(s) corresponding to last computed
reduced approximate solution(s) as parameterList.
uApprox: Approximate solution(s) with parameter(s) lastSolvedApprox as
sampleList.
lastSolvedApprox: Parameter(s) corresponding to last computed
approximate solution(s) as parameterList.
As: List of sparse matrices (in CSC format) representing coefficients
of linear system matrix.
bs: List of numpy vectors representing coefficients of linear system
RHS.
ARBs: List of sparse matrices (in CSC format) representing coefficients
of compressed linear system matrix.
bRBs: List of numpy vectors representing coefficients of compressed
linear system RHS.
"""
def __init__(self, HFEngine:HFEng, mu0 : paramVal = None,
directionPivot : ListAny = [0],
approxParameters : DictAny = {}, verbosity : int = 10,
timestamp : bool = True):
self._preInit()
self._addParametersToList(["R", "PODTolerance"], [1, -1])
super().__init__(HFEngine = HFEngine, mu0 = mu0,
directionPivot = directionPivot,
approxParameters = approxParameters,
verbosity = verbosity, timestamp = timestamp)
self._postInit()
@property
def S(self):
"""Value of S."""
return self._S
@S.setter
def S(self, S):
GenericPivotedApproximant.S.fset(self, S)
if not hasattr(self, "_R"):
self._R = np.prod(self.S) * np.prod(self.SMarginal)
@property
def R(self):
"""Value of R. Its assignment may change S."""
return self._R
@R.setter
def R(self, R):
if R < 0: raise RROMPyException("R must be non-negative.")
self._R = R
self._approxParameters["R"] = self.R
if (hasattr(self, "_S") and hasattr(self, "_SMarginal")
and np.prod(self.S) * np.prod(self.SMarginal) < self.R):
RROMPyWarning(("Prescribed S and SMarginal are too small. "
"Decreasing R."))
self.R = np.prod(self.S) * np.prod(self.SMarginal)
@property
def PODTolerance(self):
"""Value of PODTolerance."""
return self._PODTolerance
@PODTolerance.setter
def PODTolerance(self, PODTolerance):
self._PODTolerance = PODTolerance
self._approxParameters["PODTolerance"] = self.PODTolerance
def _setupProjectionMatrix(self):
"""Compute projection matrix."""
RROMPyAssert(self._mode, message = "Cannot setup numerator.")
vbMng(self, "INIT", "Starting computation of projection matrix.", 7)
if self.POD:
U, s, _ = np.linalg.svd(self.samplingEngine.RPODCoalesced)
s = s ** 2.
else:
Gramian = self.HFEngine.innerProduct(
self.samplingEngine.samplesCoalesced,
self.samplingEngine.samplesCoalesced)
U, s, _ = np.linalg.svd(Gramian)
nsamples = self.samplingEngine.samplesCoalesced.shape[1]
snorm = np.cumsum(s[::-1]) / np.sum(s)
nPODTrunc = min(nsamples - np.argmax(snorm > self.PODTolerance),
self.R)
pMat = self.samplingEngine.samplesCoalesced.dot(U[:, : nPODTrunc])
vbMng(self, "MAIN",
"Assembling {}x{} projection matrix from {} samples.".format(
*(pMat.shape), nsamples), 5)
vbMng(self, "DEL", "Done computing projection matrix.", 7)
return pMat
def _setupAffineBlocks(self):
"""Compute list of marginalized affine blocks of system."""
hasAs = hasattr(self, "AsList") and self.AsList is not None
hasbs = hasattr(self, "bsList") and self.bsList is not None
if hasAs and hasbs: return
vbMng(self, "INIT", "Computing affine blocks of system.", 10)
mu0Eff = self.mu0.data[0, self.directionPivot]
if not hasAs: self.AsList = [None] * len(self.musMarginal)
if not hasbs: self.bsList = [None] * len(self.musMarginal)
for k, muM in enumerate(self.musMarginal):
muEff = [fp] * self.npar
for jj, kk in enumerate(self.directionMarginal):
muEff[kk] = muM(0, jj)
MEnginek = MarginalProxyEngine(self.HFEngine, muEff)
if not hasAs:
self.AsList[k] = MEnginek.affineLinearSystemA(mu0Eff,
self.scaleFactorPivot)
if not hasbs:
self.bsList[k] = MEnginek.affineLinearSystemb(mu0Eff,
self.scaleFactorPivot)
vbMng(self, "DEL", "Done computing affine blocks.", 10)
def setupApprox(self):
"""Compute RB projection matrix."""
if self.checkComputedApprox():
return
RROMPyAssert(self._mode, message = "Cannot setup approximant.")
vbMng(self, "INIT", "Setting up {}.".format(self.name()), 5)
self.computeScaleFactor()
self.computeSnapshots()
self._setupAffineBlocks()
pMat = self._setupProjectionMatrix()
- ARBsList, bRBsList, pList = self.assembleReducedSystemMarginal(pMat)
+ ARBsList, bRBsList = self.assembleReducedSystemMarginal(pMat)
if self.trainedModel is None:
self.trainedModel = tModel()
self.trainedModel.verbosity = self.verbosity
self.trainedModel.timestamp = self.timestamp
data = TrainedModelPivotedData(self.trainedModel.name(), self.mu0,
- pList, self.scaleFactor,
+ pMat, self.scaleFactor,
self.HFEngine.rescalingExp,
self.directionPivot)
data.musPivot = copy(self.musPivot)
data.musMarginal = copy(self.musMarginal)
self.trainedModel.data = data
else:
self.trainedModel = self.trainedModel
- self.trainedModel.data.projMat = copy(pList)
+ self.trainedModel.data.projMat = copy(pMat)
self.trainedModel.data.marginalInterp = self._setupMarginalInterp()
self.trainedModel.data.ARBsList = ARBsList
self.trainedModel.data.bRBsList = bRBsList
self.trainedModel.data.approxParameters = copy(self.approxParameters)
vbMng(self, "DEL", "Done setting up approximant.", 5)
def assembleReducedSystemMarginal(self, pMat : sampList = None)\
- -> Tuple[List[List[Np2D]],
- List[List[Np1D]], List[Np2D]]:
+ -> Tuple[List[List[Np2D]], List[List[Np1D]]]:
"""Build affine blocks of RB linear system through projections."""
if pMat is None:
self.setupApprox()
ARBsList = self.trainedModel.data.ARBsList
bRBsList = self.trainedModel.data.bRBsList
else:
vbMng(self, "INIT", "Projecting affine terms of HF model.", 10)
ARBsList = [None] * len(self.musMarginal)
bRBsList = [None] * len(self.musMarginal)
- projList = [None] * len(self.musMarginal)
for k, (As, bs, sample) in enumerate(zip(self.AsList, self.bsList,
self.samplingEngine.samples)):
compLocal = self.HFEngine.innerProduct(sample, pMat)
- projList[k] = sample.dot(customPInv(compLocal))
+ projList = sample.dot(customPInv(compLocal))
ARBsList[k], bRBsList[k] = projectAffineDecomposition(As, bs,
- projList[k])
+ projList)
vbMng(self, "DEL", "Done projecting affine terms.", 10)
- return ARBsList, bRBsList, projList
+ return ARBsList, bRBsList
diff --git a/rrompy/reduction_methods/pole_matching/rational_interpolant_pole_matching.py b/rrompy/reduction_methods/pole_matching/rational_interpolant_pole_matching.py
index de6cc58..0f28b76 100644
--- a/rrompy/reduction_methods/pole_matching/rational_interpolant_pole_matching.py
+++ b/rrompy/reduction_methods/pole_matching/rational_interpolant_pole_matching.py
@@ -1,210 +1,209 @@
# Copyright (C) 2018 by the RROMPy authors
#
# This file is part of RROMPy.
#
# RROMPy is free software: you can redistribute it and/or modify
# it under the terms of the GNU Lesser General Public License as published by
# the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
#
# RROMPy is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU Lesser General Public License for more details.
#
# You should have received a copy of the GNU Lesser General Public License
# along with RROMPy. If not, see .
#
from copy import deepcopy as copy
import numpy as np
from rrompy.reduction_methods.pivoted.rational_interpolant_pivoted import \
RationalInterpolantPivoted
from .generic_pole_matching_approximant import GenericPoleMatchingApproximant
from rrompy.utilities.poly_fitting.polynomial import (
PolynomialInterpolator as PI)
from rrompy.reduction_methods.trained_model import (TrainedModelPivotedData,
TrainedModelPoleMatchingRational as tModel)
from rrompy.utilities.base import verbosityManager as vbMng
from rrompy.utilities.exception_manager import RROMPyAssert
__all__ = ['RationalInterpolantPoleMatching']
class RationalInterpolantPoleMatching(GenericPoleMatchingApproximant,
RationalInterpolantPivoted):
"""
ROM pivoted rational interpolant computation for parametric problems with
pole matching.
Args:
HFEngine: HF problem solver.
mu0(optional): Default parameter. Defaults to 0.
directionPivot(optional): Pivot components. Defaults to [0].
approxParameters(optional): Dictionary containing values for main
parameters of approximant. Recognized keys are:
- 'POD': whether to compute POD of snapshots; defaults to True;
- 'matchingWeight': weight for pole matching optimization; defaults
to 1;
- 'cutOffTolerance': tolerance for ignoring parasitic poles;
defaults to np.inf;
- 'cutOffType': rule for tolerance computation for parasitic poles;
defaults to 'MAGNITUDE';
- 'S': total number of pivot samples current approximant relies
upon;
- 'samplerPivot': pivot sample point generator;
- 'SMarginal': total number of marginal samples current approximant
relies upon;
- 'samplerMarginal': marginal sample point generator;
- 'polybasisPivot': type of polynomial basis for pivot
interpolation; defaults to 'MONOMIAL';
- 'polybasisMarginal': type of polynomial basis for marginal
interpolation; allowed values include 'MONOMIAL', 'CHEBYSHEV'
and 'LEGENDRE'; defaults to 'MONOMIAL';
- - 'E': number of derivatives used to compute interpolant; defaults
- to 0;
- 'M': degree of rational interpolant numerator; defaults to 0;
- 'N': degree of rational interpolant denominator; defaults to 0;
+ - 'polydegreetype': type of polynomial degree; defaults to 'TOTAL';
- 'MMarginal': degree of marginal interpolant; defaults to 0;
- 'radialDirectionalWeightsPivot': radial basis weights for pivot
numerator; defaults to 0, i.e. identity;
- 'radialDirectionalWeightsMarginal': radial basis weights for
marginal interpolant; defaults to 0, i.e. identity;
- 'interpRcondPivot': tolerance for pivot interpolation; defaults
to None;
- 'interpRcondMarginal': tolerance for marginal interpolation;
defaults to None;
- 'robustTol': tolerance for robust rational denominator
management; defaults to 0.
Defaults to empty dict.
verbosity(optional): Verbosity level. Defaults to 10.
Attributes:
HFEngine: HF problem solver.
mu0: Default parameter.
directionPivot: Pivot components.
mus: Array of snapshot parameters.
musPivot: Array of pivot snapshot parameters.
musMarginal: Array of marginal snapshot parameters.
approxParameters: Dictionary containing values for main parameters of
approximant. Recognized keys are in parameterList.
parameterListSoft: Recognized keys of soft approximant parameters:
- 'POD': whether to compute POD of snapshots;
- 'matchingWeight': weight for pole matching optimization;
- 'cutOffTolerance': tolerance for ignoring parasitic poles;
- 'cutOffType': rule for tolerance computation for parasitic poles;
- 'polybasisPivot': type of polynomial basis for pivot
interpolation;
- 'polybasisMarginal': type of polynomial basis for marginal
interpolation;
- - 'E': number of derivatives used to compute interpolant;
- 'M': degree of rational interpolant numerator;
- 'N': degree of rational interpolant denominator;
+ - 'polydegreetype': type of polynomial degree;
- 'MMarginal': degree of marginal interpolant;
- 'radialDirectionalWeightsPivot': radial basis weights for pivot
numerator;
- 'radialDirectionalWeightsMarginal': radial basis weights for
marginal interpolant;
- 'interpRcondPivot': tolerance for pivot interpolation;
- 'interpRcondMarginal': tolerance for marginal interpolation;
- 'robustTol': tolerance for robust rational denominator
management.
parameterListCritical: Recognized keys of critical approximant
parameters:
- 'S': total number of pivot samples current approximant relies
upon;
- 'samplerPivot': pivot sample point generator;
- 'SMarginal': total number of marginal samples current approximant
relies upon;
- 'samplerMarginal': marginal sample point generator.
POD: Whether to compute POD of snapshots.
matchingWeight: Weight for pole matching optimization.
cutOffTolerance: Tolerance for ignoring parasitic poles.
cutOffType: Rule for tolerance computation for parasitic poles.
S: Total number of pivot samples current approximant relies upon.
sampler: Pivot sample point generator.
polybasisPivot: Type of polynomial basis for pivot interpolation.
polybasisMarginal: Type of polynomial basis for marginal interpolation.
M: Numerator degree of approximant.
N: Denominator degree of approximant.
+ polydegreetype: Type of polynomial degree.
MMarginal: Degree of marginal interpolant.
radialDirectionalWeightsPivot: Radial basis weights for pivot
numerator.
radialDirectionalWeightsMarginal: Radial basis weights for marginal
interpolant.
interpRcondPivot: Tolerance for pivot interpolation.
interpRcondMarginal: Tolerance for marginal interpolation.
robustTol: Tolerance for robust rational denominator management.
- E: Complete derivative depth level of S.
muBoundsPivot: list of bounds for pivot parameter values.
muBoundsMarginal: list of bounds for marginal parameter values.
samplingEngine: Sampling engine.
uHF: High fidelity solution(s) with parameter(s) lastSolvedHF as
sampleList.
lastSolvedHF: Parameter(s) corresponding to last computed high fidelity
solution(s) as parameterList.
uApproxReduced: Reduced approximate solution(s) with parameter(s)
lastSolvedApprox as sampleList.
lastSolvedApproxReduced: Parameter(s) corresponding to last computed
reduced approximate solution(s) as parameterList.
uApprox: Approximate solution(s) with parameter(s) lastSolvedApprox as
sampleList.
lastSolvedApprox: Parameter(s) corresponding to last computed
approximate solution(s) as parameterList.
Q: Numpy 1D vector containing complex coefficients of approximant
denominator.
P: Numpy 2D vector whose columns are FE dofs of coefficients of
approximant numerator.
"""
def setupApprox(self):
"""
Compute rational interpolant.
SVD-based robust eigenvalue management.
"""
if self.checkComputedApprox():
return
RROMPyAssert(self._mode, message = "Cannot setup approximant.")
vbMng(self, "INIT", "Setting up {}.". format(self.name()), 5)
self.computeScaleFactor()
self.computeSnapshots()
if self.trainedModel is None:
self.trainedModel = tModel()
self.trainedModel.verbosity = self.verbosity
self.trainedModel.timestamp = self.timestamp
data = TrainedModelPivotedData(self.trainedModel.name(), self.mu0,
self.samplingEngine.samplesCoalesced,
self.scaleFactor,
self.HFEngine.rescalingExp,
self.directionPivot)
data.musPivot = copy(self.musPivot)
data.musMarginal = copy(self.musMarginal)
self.trainedModel.data = data
else:
self.trainedModel = self.trainedModel
self.trainedModel.data.projMat = copy(
self.samplingEngine.samplesCoalesced)
self.trainedModel.data.marginalInterp = self._setupMarginalInterp()
if self.N > 0:
Qs = self._setupDenominator()[0]
else:
Q = PI()
Q.npar = self.musPivot.shape[1]
Q.coeffs = np.ones(tuple([1] * Q.npar), dtype = np.complex)
Q.polybasis = self.polybasisPivot0
Qs = [Q for _ in range(len(self.musMarginal))]
self.trainedModel.data._temporary = True
self.trainedModel.data.Qs = Qs
self.trainedModel.data.Ps = self._setupNumerator()
vbMng(self, "INIT", "Matching poles.", 10)
self.trainedModel.initializeFromRational(self.HFEngine,
self.matchingWeight, self.POD)
vbMng(self, "DEL", "Done matching poles.", 10)
if not np.isinf(self.cutOffTolerance):
vbMng(self, "INIT", "Recompressing by cut-off.", 10)
msg = self.trainedModel.recompressByCutOff([-1., 1.],
self.cutOffTolerance,
self.cutOffType)
vbMng(self, "DEL", "Done recompressing." + msg, 10)
self.trainedModel.data.approxParameters = copy(self.approxParameters)
vbMng(self, "DEL", "Done setting up approximant.", 5)
diff --git a/rrompy/reduction_methods/pole_matching/reduced_basis_pole_matching.py b/rrompy/reduction_methods/pole_matching/reduced_basis_pole_matching.py
index bcc5620..4b83def 100644
--- a/rrompy/reduction_methods/pole_matching/reduced_basis_pole_matching.py
+++ b/rrompy/reduction_methods/pole_matching/reduced_basis_pole_matching.py
@@ -1,182 +1,182 @@
# Copyright (C) 2018 by the RROMPy authors
#
# This file is part of RROMPy.
#
# RROMPy is free software: you can redistribute it and/or modify
# it under the terms of the GNU Lesser General Public License as published by
# the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
#
# RROMPy is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU Lesser General Public License for more details.
#
# You should have received a copy of the GNU Lesser General Public License
# along with RROMPy. If not, see .
#
from numpy import isinf
from copy import deepcopy as copy
from rrompy.reduction_methods.pivoted.reduced_basis_pivoted import \
ReducedBasisPivoted
from .generic_pole_matching_approximant import GenericPoleMatchingApproximant
from rrompy.reduction_methods.trained_model import (TrainedModelPivotedData,
TrainedModelPoleMatchingReducedBasis as tModel)
from rrompy.utilities.base import verbosityManager as vbMng
from rrompy.utilities.exception_manager import RROMPyAssert
__all__ = ['ReducedBasisPoleMatching']
class ReducedBasisPoleMatching(GenericPoleMatchingApproximant,
ReducedBasisPivoted):
"""
ROM pivoted RB approximant computation for parametric problems with pole
matching.
Args:
HFEngine: HF problem solver.
mu0(optional): Default parameter. Defaults to 0.
directionPivot(optional): Pivot components. Defaults to [0].
approxParameters(optional): Dictionary containing values for main
parameters of approximant. Recognized keys are:
- 'POD': whether to compute POD of snapshots; defaults to True;
- 'matchingWeight': weight for pole matching optimization; defaults
to 1;
- 'cutOffTolerance': tolerance for ignoring parasitic poles;
defaults to np.inf;
- 'cutOffType': rule for tolerance computation for parasitic poles;
defaults to 'MAGNITUDE';
- 'S': total number of pivot samples current approximant relies
upon;
- 'samplerPivot': pivot sample point generator;
- 'SMarginal': total number of marginal samples current approximant
relies upon;
- 'samplerMarginal': marginal sample point generator;
- 'R': rank for pivot Galerkin projection; defaults to prod(S);
- 'PODTolerance': tolerance for pivot snapshots POD; defaults to
-1;
- 'polybasisMarginal': type of polynomial basis for marginal
interpolation; allowed values include 'MONOMIAL', 'CHEBYSHEV'
and 'LEGENDRE'; defaults to 'MONOMIAL';
- 'MMarginal': degree of marginal interpolant; defaults to 0;
- 'radialDirectionalWeightsMarginal': radial basis weights for
marginal interpolant; defaults to 0, i.e. identity;
- 'interpRcondMarginal': tolerance for marginal interpolation;
defaults to None.
Defaults to empty dict.
verbosity(optional): Verbosity level. Defaults to 10.
Attributes:
HFEngine: HF problem solver.
mu0: Default parameter.
directionPivot: Pivot components.
mus: Array of snapshot parameters.
musPivot: Array of pivot snapshot parameters.
musMarginal: Array of marginal snapshot parameters.
approxParameters: Dictionary containing values for main parameters of
approximant. Recognized keys are in parameterList.
parameterListSoft: Recognized keys of soft approximant parameters:
- 'POD': whether to compute POD of snapshots;
- 'matchingWeight': weight for pole matching optimization;
- 'cutOffTolerance': tolerance for ignoring parasitic poles;
- 'cutOffType': rule for tolerance computation for parasitic poles;
- 'R': rank for Galerkin projection;
- 'PODTolerance': tolerance for snapshots POD;
- 'polybasisMarginal': type of polynomial basis for marginal
interpolation;
- 'MMarginal': degree of marginal interpolant;
- 'radialDirectionalWeightsMarginal': radial basis weights for
marginal interpolant;
- 'interpRcondMarginal': tolerance for marginal interpolation.
parameterListCritical: Recognized keys of critical approximant
parameters:
- 'S': total number of pivot samples current approximant relies
upon;
- 'samplerPivot': pivot sample point generator;
- 'SMarginal': total number of marginal samples current approximant
relies upon;
- 'samplerMarginal': marginal sample point generator.
POD: Whether to compute POD of snapshots.
matchingWeight: Weight for pole matching optimization.
cutOffTolerance: Tolerance for ignoring parasitic poles.
cutOffType: Rule for tolerance computation for parasitic poles.
S: Total number of pivot samples current approximant relies upon.
samplerPivot: Pivot sample point generator.
SMarginal: Total number of marginal samples current approximant relies
upon.
samplerMarginal: Marginal sample point generator.
R: Rank for Galerkin projection.
PODTolerance: Tolerance for pivot snapshots POD.
polybasisMarginal: Type of polynomial basis for marginal interpolation.
MMarginal: Degree of marginal interpolant.
radialDirectionalWeightsMarginal: Radial basis weights for marginal
interpolant.
interpRcondMarginal: Tolerance for marginal interpolation.
muBoundsPivot: list of bounds for pivot parameter values.
muBoundsMarginal: list of bounds for marginal parameter values.
samplingEngine: Sampling engine.
uHF: High fidelity solution(s) with parameter(s) lastSolvedHF as
sampleList.
lastSolvedHF: Parameter(s) corresponding to last computed high fidelity
solution(s) as parameterList.
lastSolvedHF: Parameter(s) corresponding to last computed high fidelity
solution(s) as parameterList.
uApproxReduced: Reduced approximate solution(s) with parameter(s)
lastSolvedApprox as sampleList.
lastSolvedApproxReduced: Parameter(s) corresponding to last computed
reduced approximate solution(s) as parameterList.
uApprox: Approximate solution(s) with parameter(s) lastSolvedApprox as
sampleList.
lastSolvedApprox: Parameter(s) corresponding to last computed
approximate solution(s) as parameterList.
As: List of sparse matrices (in CSC format) representing coefficients
of linear system matrix.
bs: List of numpy vectors representing coefficients of linear system
RHS.
ARBs: List of sparse matrices (in CSC format) representing coefficients
of compressed linear system matrix.
bRBs: List of numpy vectors representing coefficients of compressed
linear system RHS.
"""
def setupApprox(self):
"""Compute RB projection matrix."""
if self.checkComputedApprox():
return
RROMPyAssert(self._mode, message = "Cannot setup approximant.")
vbMng(self, "INIT", "Setting up {}.".format(self.name()), 5)
self.computeScaleFactor()
self.computeSnapshots()
self._setupAffineBlocks()
pMat = self._setupProjectionMatrix()
- ARBsList, bRBsList, pList = self.assembleReducedSystemMarginal(pMat)
+ ARBsList, bRBsList = self.assembleReducedSystemMarginal(pMat)
if self.trainedModel is None:
self.trainedModel = tModel()
self.trainedModel.verbosity = self.verbosity
self.trainedModel.timestamp = self.timestamp
data = TrainedModelPivotedData(self.trainedModel.name(), self.mu0,
- pList, self.scaleFactor,
+ pMat, self.scaleFactor,
self.HFEngine.rescalingExp,
self.directionPivot)
data.musPivot = copy(self.musPivot)
data.musMarginal = copy(self.musMarginal)
self.trainedModel.data = data
else:
self.trainedModel = self.trainedModel
- self.trainedModel.data.projMat = copy(pList)
+ self.trainedModel.data.projMat = copy(pMat)
self.trainedModel.data.marginalInterp = self._setupMarginalInterp()
self.trainedModel.data.ARBsList = ARBsList
self.trainedModel.data.bRBsList = bRBsList
vbMng(self, "INIT", "Matching poles.", 10)
self.trainedModel.initializeFromAffine(self.HFEngine,
self.matchingWeight, self.POD)
vbMng(self, "DEL", "Done matching poles.", 10)
if not isinf(self.cutOffTolerance):
vbMng(self, "INIT", "Recompressing by cut-off.", 10)
msg = self.trainedModel.recompressByCutOff([- 1., 1.],
self.cutOffTolerance,
self.cutOffType)
vbMng(self, "DEL", "Done recompressing." + msg, 10)
self.trainedModel.data.approxParameters = copy(self.approxParameters)
vbMng(self, "DEL", "Done setting up approximant.", 5)
diff --git a/rrompy/reduction_methods/standard/rational_interpolant.py b/rrompy/reduction_methods/standard/rational_interpolant.py
index 95e74e8..940509c 100644
--- a/rrompy/reduction_methods/standard/rational_interpolant.py
+++ b/rrompy/reduction_methods/standard/rational_interpolant.py
@@ -1,544 +1,554 @@
# Copyright (C) 2018 by the RROMPy authors
#
# This file is part of RROMPy.
#
# RROMPy is free software: you can redistribute it and/or modify
# it under the terms of the GNU Lesser General Public License as published by
# the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
#
# RROMPy is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU Lesser General Public License for more details.
#
# You should have received a copy of the GNU Lesser General Public License
# along with RROMPy. If not, see .
#
from copy import deepcopy as copy
import numpy as np
from rrompy.reduction_methods.base import checkRobustTolerance
from .generic_standard_approximant import GenericStandardApproximant
from rrompy.utilities.poly_fitting.polynomial import (polybases as ppb,
polyfitname,
nextDerivativeIndices,
hashDerivativeToIdx as hashD,
hashIdxToDerivative as hashI,
+ polyvander as pvP,
polyvanderTotal as pvTP,
degreeTotalToFull,
+ fullDegreeMaxMask,
+ totalDegreeMaxMask,
PolynomialInterpolator as PI)
from rrompy.utilities.poly_fitting.radial_basis import (polybases as rbpb,
RadialBasisInterpolator as RBI)
from rrompy.utilities.poly_fitting.moving_least_squares import (
polybases as mlspb,
MovingLeastSquaresInterpolator as MLSI)
from rrompy.reduction_methods.trained_model import (
TrainedModelRational as tModel)
from rrompy.reduction_methods.trained_model import TrainedModelData
from rrompy.utilities.base.types import (Np1D, Np2D, HFEng, DictAny, Tuple,
List, paramVal, paramList, sampList)
from rrompy.utilities.base import verbosityManager as vbMng
from rrompy.utilities.numerical import multifactorial, customPInv
from rrompy.utilities.exception_manager import (RROMPyException, RROMPyAssert,
RROMPyWarning)
__all__ = ['RationalInterpolant']
class RationalInterpolant(GenericStandardApproximant):
"""
ROM rational interpolant computation for parametric problems.
Args:
HFEngine: HF problem solver.
mu0(optional): Default parameter. Defaults to 0.
approxParameters(optional): Dictionary containing values for main
parameters of approximant. Recognized keys are:
- 'POD': whether to compute POD of snapshots; defaults to True;
- 'S': total number of samples current approximant relies upon;
- 'sampler': sample point generator;
- 'polybasis': type of polynomial basis for interpolation; defaults
to 'MONOMIAL';
- - 'E': number of derivatives used to compute interpolant; defaults
- to 0;
- 'M': degree of rational interpolant numerator; defaults to 0;
- 'N': degree of rational interpolant denominator; defaults to 0;
+ - 'polydegreetype': type of polynomial degree; defaults to 'TOTAL';
- 'radialDirectionalWeights': radial basis weights for interpolant
numerator; defaults to 0, i.e. identity;
- '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;
- 'polybasis': type of polynomial basis for interpolation;
- - 'E': number of derivatives used to compute interpolant;
- 'M': degree of rational interpolant numerator;
- 'N': degree of rational interpolant denominator;
+ - 'polydegreetype': type of polynomial degree;
- 'radialDirectionalWeights': radial basis weights for interpolant
numerator;
- 'interpRcond': tolerance for interpolation via numpy.polyfit;
- 'robustTol': tolerance for robust rational denominator
management.
parameterListCritical: Recognized keys of critical approximant
parameters:
- 'S': total number of samples current approximant relies upon;
- 'sampler': sample point generator.
POD: Whether to compute POD of snapshots.
S: Number of solution snapshots over which current approximant is
based upon.
sampler: Sample point generator.
polybasis: type of polynomial basis for interpolation.
M: Numerator degree of approximant.
N: Denominator degree of approximant.
+ polydegreetype: Type of polynomial degree.
radialDirectionalWeights: Radial basis weights for interpolant
numerator.
interpRcond: Tolerance for interpolation via numpy.polyfit.
robustTol: Tolerance for robust rational denominator management.
- E: Complete derivative depth level of S.
muBounds: list of bounds for parameter values.
samplingEngine: Sampling engine.
uHF: High fidelity solution(s) with parameter(s) lastSolvedHF as
sampleList.
lastSolvedHF: Parameter(s) corresponding to last computed high fidelity
solution(s) as parameterList.
uApproxReduced: Reduced approximate solution(s) with parameter(s)
lastSolvedApprox as sampleList.
lastSolvedApproxReduced: Parameter(s) corresponding to last computed
reduced approximate solution(s) as parameterList.
uApprox: Approximate solution(s) with parameter(s) lastSolvedApprox as
sampleList.
lastSolvedApprox: Parameter(s) corresponding to last computed
approximate solution(s) as parameterList.
Q: Numpy 1D vector containing complex coefficients of approximant
denominator.
P: Numpy 2D vector whose columns are FE dofs of coefficients of
approximant numerator.
"""
def __init__(self, HFEngine:HFEng, mu0 : paramVal = None,
approxParameters : DictAny = {}, verbosity : int = 10,
timestamp : bool = True):
self._preInit()
- self._addParametersToList(["polybasis", "E", "M", "N",
+ self._addParametersToList(["polybasis", "M", "N", "polydegreetype",
"radialDirectionalWeights", "interpRcond",
"robustTol"],
- ["MONOMIAL", -1, 0, 0, 1, -1, 0])
+ ["MONOMIAL", 0, 0, "TOTAL", 1, -1, 0])
super().__init__(HFEngine = HFEngine, mu0 = mu0,
approxParameters = approxParameters,
verbosity = verbosity, timestamp = timestamp)
self.catchInstability = False
self._postInit()
@property
def polybasis(self):
"""Value of polybasis."""
return self._polybasis
@polybasis.setter
def polybasis(self, polybasis):
try:
polybasis = polybasis.upper().strip().replace(" ","")
if polybasis not in ppb + rbpb + mlspb:
raise RROMPyException("Prescribed polybasis not recognized.")
self._polybasis = polybasis
except:
RROMPyWarning(("Prescribed polybasis not recognized. Overriding "
"to 'MONOMIAL'."))
self._polybasis = "MONOMIAL"
self._approxParameters["polybasis"] = self.polybasis
@property
def polybasis0(self):
if "_" in self.polybasis:
return self.polybasis.split("_")[0]
return self.polybasis
@property
def interpRcond(self):
"""Value of interpRcond."""
return self._interpRcond
@interpRcond.setter
def interpRcond(self, interpRcond):
self._interpRcond = interpRcond
self._approxParameters["interpRcond"] = self.interpRcond
@property
def radialDirectionalWeights(self):
"""Value of radialDirectionalWeights."""
return self._radialDirectionalWeights
@radialDirectionalWeights.setter
def radialDirectionalWeights(self, radialDirectionalWeights):
self._radialDirectionalWeights = radialDirectionalWeights
self._approxParameters["radialDirectionalWeights"] = (
self.radialDirectionalWeights)
@property
def M(self):
- """Value of M. Its assignment may change S."""
+ """Value of M."""
return self._M
@M.setter
def M(self, M):
if M < 0: raise RROMPyException("M must be non-negative.")
self._M = M
self._approxParameters["M"] = self.M
- if hasattr(self, "_E") and self.E >= 0 and self.E < self.M:
- RROMPyWarning("Prescribed S is too small. Decreasing M.")
- self.M = self.E
@property
def N(self):
- """Value of N. Its assignment may change S."""
+ """Value of N."""
return self._N
@N.setter
def N(self, N):
if N < 0: raise RROMPyException("N must be non-negative.")
self._N = N
self._approxParameters["N"] = self.N
- if hasattr(self, "_E") and self.E >= 0 and self.E < self.N:
- RROMPyWarning("Prescribed S is too small. Decreasing N.")
- self.N = self.E
@property
- def E(self):
- """Value of E."""
- return self._E
- @E.setter
- def E(self, E):
- if E < 0:
- if not hasattr(self, "_S"):
- raise RROMPyException(("Value of E must be positive if S is "
- "not yet assigned."))
- E = np.sum(hashI(np.prod(self.S), self.npar)) - 1
- self._E = E
- self._approxParameters["E"] = self.E
- if (hasattr(self, "_S")
- and self.E >= np.sum(hashI(np.prod(self.S), self.npar))):
- RROMPyWarning("Prescribed S is too small. Decreasing E.")
- self.E = -1
- if hasattr(self, "_M"): self.M = self.M
- if hasattr(self, "_N"): self.N = self.N
+ def polydegreetype(self):
+ """Value of polydegreetype."""
+ return self._polydegreetype
+ @polydegreetype.setter
+ def polydegreetype(self, polydegreetype):
+ try:
+ polydegreetype = polydegreetype.upper().strip().replace(" ","")
+ if polydegreetype not in ["TOTAL", "FULL"]:
+ raise RROMPyException(("Prescribed polydegreetype not "
+ "recognized."))
+ self._polydegreetype = polydegreetype
+ except:
+ RROMPyWarning(("Prescribed polydegreetype not recognized. "
+ "Overriding to 'TOTAL'."))
+ self._polydegreetype = "TOTAL"
+ self._approxParameters["polydegreetype"] = self.polydegreetype
@property
def robustTol(self):
"""Value of tolerance for robust rational denominator management."""
return self._robustTol
@robustTol.setter
def robustTol(self, robustTol):
if robustTol < 0.:
RROMPyWarning(("Overriding prescribed negative robustness "
"tolerance to 0."))
robustTol = 0.
self._robustTol = robustTol
self._approxParameters["robustTol"] = self.robustTol
@property
def S(self):
"""Value of S."""
return self._S
@S.setter
def S(self, S):
GenericStandardApproximant.S.fset(self, S)
if hasattr(self, "_M"): self.M = self.M
if hasattr(self, "_N"): self.N = self.N
- if hasattr(self, "_E"): self.E = self.E
def resetSamples(self):
"""Reset samples."""
super().resetSamples()
self._musUniqueCN = None
self._derIdxs = None
self._reorder = None
def _setupInterpolationIndices(self):
"""Setup parameters for polyvander."""
RROMPyAssert(self._mode,
message = "Cannot setup interpolation indices.")
if self._musUniqueCN is None or len(self._reorder) != len(self.mus):
self._musUniqueCN, musIdxsTo, musIdxs, musCount = (
self.centerNormalize(self.mus).unique(return_index = True,
return_inverse = True,
return_counts = True))
self._musUnique = self.mus[musIdxsTo]
self._derIdxs = [None] * len(self._musUniqueCN)
self._reorder = np.empty(len(musIdxs), dtype = int)
filled = 0
for j, cnt in enumerate(musCount):
self._derIdxs[j] = nextDerivativeIndices([], self.mus.shape[1],
cnt)
jIdx = np.nonzero(musIdxs == j)[0]
self._reorder[jIdx] = np.arange(filled, filled + cnt)
filled += cnt
def _setupDenominator(self):
"""Compute rational denominator."""
RROMPyAssert(self._mode, message = "Cannot setup denominator.")
vbMng(self, "INIT", "Starting computation of denominator.", 7)
while self.N > 0:
invD, fitinv = self._computeInterpolantInverseBlocks()
if self.POD:
ev, eV = self.findeveVGQR(self.samplingEngine.RPOD, invD)
else:
ev, eV = self.findeveVGExplicit(self.samplingEngine.samples,
invD)
nevBad = checkRobustTolerance(ev, self.robustTol)
if nevBad <= 1: break
if self.catchInstability:
raise RROMPyException(("Instability in denominator "
"computation: eigenproblem is poorly "
"conditioned."))
RROMPyWarning(("Smallest {} eigenvalues below tolerance. Reducing "
"N by 1.").format(nevBad))
self.N = self.N - 1
if self.N <= 0:
self._N = 0
eV = np.ones((1, 1))
q = PI()
q.npar = self.npar
q.polybasis = self.polybasis0
- q.coeffs = degreeTotalToFull(tuple([self.N + 1] * self.npar),
- self.npar, eV[:, 0])
+ if self.polydegreetype == "TOTAL":
+ q.coeffs = degreeTotalToFull(tuple([self.N + 1] * self.npar),
+ self.npar, eV[:, 0])
+ else:
+ q.coeffs = eV[:, 0].reshape([self.N + 1] * self.npar)
vbMng(self, "DEL", "Done computing denominator.", 7)
return q, fitinv
def _setupNumerator(self):
"""Compute rational numerator."""
RROMPyAssert(self._mode, message = "Cannot setup numerator.")
vbMng(self, "INIT", "Starting computation of numerator.", 7)
Qevaldiag = np.zeros((len(self.mus), len(self.mus)),
dtype = np.complex)
verb = self.trainedModel.verbosity
self.trainedModel.verbosity = 0
-
self._setupInterpolationIndices()
idxGlob = 0
for j, derIdxs in enumerate(self._derIdxs):
nder = len(derIdxs)
idxLoc = np.arange(len(self.mus))[(self._reorder >= idxGlob)
* (self._reorder < idxGlob + nder)]
idxGlob += nder
Qval = [0] * nder
for der in range(nder):
derIdx = hashI(der, self.npar)
Qval[der] = (self.trainedModel.getQVal(
self._musUnique[j], derIdx,
scl = np.power(self.scaleFactor, -1.))
/ multifactorial(derIdx))
for derU, derUIdx in enumerate(derIdxs):
for derQ, derQIdx in enumerate(derIdxs):
diffIdx = [x - y for (x, y) in zip(derUIdx, derQIdx)]
if all([x >= 0 for x in diffIdx]):
diffj = hashD(diffIdx)
Qevaldiag[idxLoc[derU], idxLoc[derQ]] = Qval[diffj]
if self.POD:
Qevaldiag = Qevaldiag.dot(self.samplingEngine.RPOD.T)
self.trainedModel.verbosity = verb
while self.M >= 0:
if self.polybasis in ppb:
p = PI()
wellCond, msg = p.setupByInterpolation(
- self._musUniqueCN, Qevaldiag, self.M,
- self.polybasis, self.verbosity >= 5, True,
- {"derIdxs": self._derIdxs, "reorder": self._reorder,
- "scl": np.power(self.scaleFactor, -1.)},
- {"rcond": self.interpRcond})
+ self._musUniqueCN, Qevaldiag, self.M,
+ self.polybasis, self.verbosity >= 5,
+ self.polydegreetype == "TOTAL",
+ {"derIdxs": self._derIdxs,
+ "reorder": self._reorder,
+ "scl": np.power(self.scaleFactor, -1.)},
+ {"rcond": self.interpRcond})
elif self.polybasis in rbpb:
p = RBI()
wellCond, msg = p.setupByInterpolation(
- self._musUniqueCN, Qevaldiag, self.M, self.polybasis,
- self.radialDirectionalWeights, self.verbosity >= 5, True,
- {"derIdxs": self._derIdxs, "reorder": self._reorder,
- "scl": np.power(self.scaleFactor, -1.)},
- {"rcond": self.interpRcond})
+ self._musUniqueCN, Qevaldiag, self.M, self.polybasis,
+ self.radialDirectionalWeights, self.verbosity >= 5,
+ self.polydegreetype == "TOTAL",
+ {"derIdxs": self._derIdxs, "reorder": self._reorder,
+ "scl": np.power(self.scaleFactor, -1.)},
+ {"rcond": self.interpRcond})
else:# if self.polybasis in mlspb:
p = MLSI()
wellCond, msg = p.setupByInterpolation(
- self._musUniqueCN, Qevaldiag, self.M, self.polybasis,
- self.radialDirectionalWeights, self.verbosity >= 5, True,
- {"derIdxs": self._derIdxs, "reorder": self._reorder,
- "scl": np.power(self.scaleFactor, -1.)})
+ self._musUniqueCN, Qevaldiag, self.M, self.polybasis,
+ self.radialDirectionalWeights, self.verbosity >= 5,
+ self.polydegreetype == "TOTAL",
+ {"derIdxs": self._derIdxs, "reorder": self._reorder,
+ "scl": np.power(self.scaleFactor, -1.)})
vbMng(self, "MAIN", msg, 5)
if wellCond: break
if self.catchInstability:
raise RROMPyException(("Instability in numerator computation: "
"polyfit is poorly conditioned."))
RROMPyWarning("Polyfit is poorly conditioned. Reducing M by 1.")
self.M -= 1
if self.M < 0:
raise RROMPyException(("Instability in computation of numerator. "
"Aborting."))
vbMng(self, "DEL", "Done computing numerator.", 7)
return p
def setupApprox(self):
"""
Compute rational interpolant.
SVD-based robust eigenvalue management.
"""
if self.checkComputedApprox():
return
RROMPyAssert(self._mode, message = "Cannot setup approximant.")
vbMng(self, "INIT", "Setting up {}.". format(self.name()), 5)
self.computeScaleFactor()
self.computeSnapshots()
if self.trainedModel is None:
self.trainedModel = tModel()
self.trainedModel.verbosity = self.verbosity
self.trainedModel.timestamp = self.timestamp
data = 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)
if self.N > 0:
Q = self._setupDenominator()[0]
else:
Q = PI()
Q.coeffs = np.ones(tuple([1] * self.npar), dtype = np.complex)
Q.npar = self.npar
Q.polybasis = self.polybasis
self.trainedModel.data.Q = Q
self.trainedModel.data.P = self._setupNumerator()
self.trainedModel.data.approxParameters = copy(self.approxParameters)
vbMng(self, "DEL", "Done setting up approximant.", 5)
def _computeInterpolantInverseBlocks(self) -> Tuple[List[Np2D], Np2D]:
"""
Compute inverse factors for minimal interpolant target functional.
"""
RROMPyAssert(self._mode, message = "Cannot solve eigenvalue problem.")
self._setupInterpolationIndices()
- while self.E >= 0:
- eWidth = (hashD([self.E + 1] + [0] * (self.npar - 1))
- - hashD([self.E] + [0] * (self.npar - 1)))
- TE, _, argIdxs = pvTP(self._musUniqueCN, self.E, self.polybasis0,
- self._derIdxs, self._reorder,
- scl = np.power(self.scaleFactor, -1.))
- fitOut = customPInv(TE[:, argIdxs], rcond = self.interpRcond,
- full = True)
+ while self.N >= 0:
+ if self.polydegreetype == "TOTAL":
+ TE, _, argIdxs = pvTP(self._musUniqueCN, self.N,
+ self.polybasis0, self._derIdxs,
+ self._reorder,
+ scl = np.power(self.scaleFactor, -1.))
+ TE = TE[:, argIdxs]
+ idxsB = totalDegreeMaxMask(self.N, self.npar)
+ else: #if self.polydegreetype == "FULL":
+ TE = pvP(self._musUniqueCN, [self.N] * self.npar,
+ self.polybasis0, self._derIdxs, self._reorder,
+ scl = np.power(self.scaleFactor, -1.))
+ idxsB = fullDegreeMaxMask(self.N, self.npar)
+ fitOut = customPInv(TE, rcond = self.interpRcond, full = True)
vbMng(self, "MAIN",
("Fitting {} samples with degree {} through {}... "
"Conditioning of pseudoinverse system: {:.4e}.").format(
- TE.shape[0], self.E,
+ TE.shape[0], self.N,
polyfitname(self.polybasis0),
fitOut[1][1][0] / fitOut[1][1][-1]),
5)
- if fitOut[1][0] == len(argIdxs):
- fitinv = fitOut[0][- eWidth : , :]
+ if fitOut[1][0] == TE.shape[1]:
+ fitinv = fitOut[0][idxsB, :]
break
- RROMPyWarning("Polyfit is poorly conditioned. Reducing E by 1.")
- self.E -= 1
- if self.E < 0:
+ RROMPyWarning("Polyfit is poorly conditioned. Reducing N by 1.")
+ self.N -= 1
+ if self.N < 0:
raise RROMPyException(("Instability in computation of "
"denominator. Aborting."))
- TN, _, argIdxs = pvTP(self._musUniqueCN, self.N, self.polybasis0,
- self._derIdxs, self._reorder,
- scl = np.power(self.scaleFactor, -1.))
- TN = TN[:, argIdxs]
- invD = [None] * (eWidth)
- for k in range(eWidth):
+ if self.polydegreetype == "TOTAL":
+ TN, _, argIdxs = pvTP(self._musUniqueCN, self.N, self.polybasis0,
+ self._derIdxs, self._reorder,
+ scl = np.power(self.scaleFactor, -1.))
+ TN = TN[:, argIdxs]
+ else: #if self.polydegreetype == "FULL":
+ TN = pvP(self._musUniqueCN, [self.N] * self.npar,
+ self.polybasis0, self._derIdxs, self._reorder,
+ scl = np.power(self.scaleFactor, -1.))
+ invD = [None] * (len(idxsB))
+ for k in range(len(idxsB)):
pseudoInv = np.diag(fitinv[k, :])
idxGlob = 0
for j, derIdxs in enumerate(self._derIdxs):
nder = len(derIdxs)
idxGlob += nder
if nder > 1:
idxLoc = np.arange(len(self.mus))[
(self._reorder >= idxGlob - nder)
* (self._reorder < idxGlob)]
invLoc = fitinv[k, idxLoc]
pseudoInv[np.ix_(idxLoc, idxLoc)] = 0.
for diffj, diffjIdx in enumerate(derIdxs):
for derQ, derQIdx in enumerate(derIdxs):
derUIdx = [x - y for (x, y) in
- zip(diffjIdx, derQIdx)]
+ zip(diffjIdx, derQIdx)]
if all([x >= 0 for x in derUIdx]):
derU = hashD(derUIdx)
pseudoInv[idxLoc[derU], idxLoc[derQ]] = (
- invLoc[diffj])
+ invLoc[diffj])
invD[k] = pseudoInv.dot(TN)
return invD, fitinv
def findeveVGExplicit(self, sampleE:sampList,
invD:List[Np2D]) -> Tuple[Np1D, Np2D]:
"""
Compute explicitly eigenvalues and eigenvectors of rational denominator
matrix.
"""
RROMPyAssert(self._mode, message = "Cannot solve eigenvalue problem.")
nEnd = invD[0].shape[1]
eWidth = len(invD)
vbMng(self, "INIT", "Building gramian matrix.", 10)
gramian = self.HFEngine.innerProduct(sampleE, sampleE)
G = np.zeros((nEnd, nEnd), dtype = np.complex)
for k in range(eWidth):
G += invD[k].T.conj().dot(gramian.dot(invD[k]))
vbMng(self, "DEL", "Done building gramian.", 10)
vbMng(self, "INIT", "Solving eigenvalue problem for gramian matrix.",
7)
ev, eV = np.linalg.eigh(G)
vbMng(self, "MAIN",
("Solved eigenvalue problem of size {} with condition number "
"{:.4e}.").format(nEnd, ev[-1] / ev[0]), 5)
vbMng(self, "DEL", "Done solving eigenvalue problem.", 7)
return ev, eV
def findeveVGQR(self, RPODE:Np2D, invD:List[Np2D]) -> Tuple[Np1D, Np2D]:
"""
Compute eigenvalues and eigenvectors of rational denominator matrix
through SVD of R factor.
"""
RROMPyAssert(self._mode, message = "Cannot solve eigenvalue problem.")
nEnd = invD[0].shape[1]
S = RPODE.shape[0]
eWidth = len(invD)
vbMng(self, "INIT", "Building half-gramian matrix stack.", 10)
Rstack = np.zeros((S * eWidth, nEnd), dtype = np.complex)
for k in range(eWidth):
Rstack[k * S : (k + 1) * S, :] = RPODE.dot(invD[k])
vbMng(self, "DEL", "Done building half-gramian.", 10)
vbMng(self, "INIT", "Solving svd for square root of gramian matrix.",
7)
_, s, eV = np.linalg.svd(Rstack, full_matrices = False)
ev = s[::-1]
eV = eV[::-1, :].T.conj()
vbMng(self, "MAIN",
("Solved svd problem of size {} x {} with condition number "
"{:.4e}.").format(*Rstack.shape, s[0] / s[-1]), 5)
vbMng(self, "DEL", "Done solving svd.", 7)
return ev, eV
def centerNormalize(self, mu : paramList = [],
mu0 : paramVal = None) -> paramList:
"""
Compute normalized parameter to be plugged into approximant.
Args:
mu: Parameter(s) 1.
mu0: Parameter(s) 2. If None, set to self.mu0.
Returns:
Normalized parameter.
"""
return self.trainedModel.centerNormalize(mu, mu0)
def getResidues(self, *args, **kwargs) -> Np1D:
"""
Obtain approximant residues.
Returns:
Matrix with residues as columns.
"""
return self.trainedModel.getResidues(*args, **kwargs)
diff --git a/rrompy/reduction_methods/trained_model/trained_model_pivoted_reduced_basis.py b/rrompy/reduction_methods/trained_model/trained_model_pivoted_reduced_basis.py
index a6beccb..1d1b5fc 100644
--- a/rrompy/reduction_methods/trained_model/trained_model_pivoted_reduced_basis.py
+++ b/rrompy/reduction_methods/trained_model/trained_model_pivoted_reduced_basis.py
@@ -1,166 +1,148 @@
# Copyright (C) 2018 by the RROMPy authors
#
# This file is part of RROMPy.
#
# RROMPy is free software: you can redistribute it and/or modify
# it under the terms of the GNU Lesser General Public License as published by
# the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
#
# RROMPy is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU Lesser General Public License for more details.
#
# You should have received a copy of the GNU Lesser General Public License
# along with RROMPy. If not, see .
#
import numpy as np
from .trained_model_pivoted_general import TrainedModelPivotedGeneral
from .trained_model_reduced_basis import TrainedModelReducedBasis
from rrompy.utilities.base.types import (Np1D, ListAny, paramVal, paramList,
sampList)
from rrompy.utilities.base import verbosityManager as vbMng, freepar as fp
from rrompy.utilities.poly_fitting.polynomial import (
hashIdxToDerivative as hashI)
from rrompy.utilities.numerical import (eigvalsNonlinearDense,
marginalizePolyList)
from rrompy.utilities.exception_manager import RROMPyException
from rrompy.parameter import checkParameter, checkParameterList
from rrompy.sampling import emptySampleList, sampleList
__all__ = ['TrainedModelPivotedReducedBasis']
class TrainedModelPivotedReducedBasis(TrainedModelPivotedGeneral,
TrainedModelReducedBasis):
"""
ROM approximant evaluation for pivoted RB approximant.
Attributes:
Data: dictionary with all that can be pickled.
"""
def interpolateMarginal(self, mu : paramVal = [],
samples : ListAny = []) -> ListAny:
"""
Evaluate marginal interpolator at arbitrary marginal parameter.
Args:
mu: Target parameter.
samples: Objects to interpolate.
"""
mu = checkParameter(mu, self.data.nparMarginal)
sList = isinstance(samples[0], sampleList)
sEff = [None] * len(samples)
for j in range(len(samples)):
if sList:
sEff[j] = samples[j].data
else:
sEff[j] = samples[j]
+ try:
+ dtype = sEff[0].dtype
+ except:
+ dtype = sEff[0][0].dtype
vbMng(self, "INIT", "Interpolating marginal at mu = {}.".format(mu),
95)
muC = self.centerNormalizeMarginal(mu)
- p = np.zeros_like(sEff[0], dtype = sEff[0].dtype)
+ p = np.zeros_like(sEff[0], dtype = dtype)
for mIj, spj in zip(self.data.marginalInterp, sEff):
p += mIj(muC) * spj
vbMng(self, "DEL", "Done interpolating marginal.", 95)
-# if sList: p = sampleList(p)
return p
def getApproxReduced(self, mu : paramList = []) -> sampList:
"""
Evaluate reduced representation of approximant at arbitrary parameter.
Args:
mu: Target parameter.
"""
mu = checkParameterList(mu, self.data.npar)[0]
if (not hasattr(self, "lastSolvedApproxReduced")
or self.lastSolvedApproxReduced != mu):
vbMng(self, "INIT",
"Computing RB solution at mu = {}.".format(mu), 12)
self.uApproxReduced = emptySampleList()
self.uApproxReduced.reset((self.data.ARBsList[0][0].shape[0],
- len(mu)), self.data.projMat[0].dtype)
+ len(mu)), self.data.projMat.dtype)
for i, muPL in enumerate(mu):
muP = checkParameter([muPL(0, x) \
for x in self.data.directionPivot])
muM = [muPL(0, x) for x in self.data.directionMarginal]
vbMng(self, "INIT",
"Assembling reduced model for mu = {}.".format(muPL), 87)
ARBs = self.interpolateMarginal(muM, self.data.ARBsList)
bRBs = self.interpolateMarginal(muM, self.data.bRBsList)
ARBmu = ARBs[0]
bRBmu = bRBs[0]
for j in range(1, max(len(ARBs), len(bRBs))):
theta = np.prod(self.centerNormalizePivot(muP)
** hashI(j, self.data.nparPivot))
if j < len(ARBs):
ARBmu += theta * ARBs[j]
if j < len(bRBs):
bRBmu += theta * bRBs[j]
vbMng(self, "DEL", "Done assembling reduced model.", 87)
vbMng(self, "INIT",
"Solving reduced model for mu = {}.".format(muPL), 75)
self.uApproxReduced[i] = np.linalg.solve(ARBmu, bRBmu)
vbMng(self, "DEL", "Done solving reduced model.", 75)
vbMng(self, "DEL", "Done computing RB solution.", 12)
self.lastSolvedApproxReduced = mu
return self.uApproxReduced
- def getApprox(self, mu : paramList = []) -> sampList:
- """
- Evaluate approximant at arbitrary parameter.
-
- Args:
- mu: Target parameter.
- """
- mu = checkParameterList(mu, self.data.npar)[0]
- if (not hasattr(self, "lastSolvedApprox")
- or self.lastSolvedApprox != mu):
- uApproxR = self.getApproxReduced(mu)
- self.uApprox = emptySampleList()
- self.uApprox.reset((self.data.projMat[0].shape[0], len(mu)),
- self.data.projMat[0].dtype)
- for i in range(len(mu)):
- muM = [mu(i, x) for x in self.data.directionMarginal]
- projLoc = self.interpolateMarginal(muM, self.data.projMat)[0]
- self.uApprox[i] = projLoc.dot(uApproxR[i])
- self.lastSolvedApprox = mu
- return self.uApprox
-
def getPoles(self, marginalVals : ListAny = [fp], jSupp : int = 1,
**kwargs) -> Np1D:
"""
Obtain approximant poles.
Returns:
Numpy complex vector of poles.
"""
if not hasattr(marginalVals, "__len__"): marginalVals = [marginalVals]
marginalVals = list(marginalVals)
try:
rDim = marginalVals.index(fp)
if rDim < len(marginalVals) - 1 and fp in marginalVals[rDim + 1 :]:
raise
except:
raise RROMPyException(("Exactly 1 'freepar' entry in "
"marginalVals must be provided."))
if rDim in self.data.directionMarginal:
raise RROMPyException(("'freepar' entry in marginalVals must be "
"in pivot dimension."))
mValsP = [marginalVals[x] for x in self.data.directionPivot]
rDimP = mValsP.index(fp)
mValsM = [marginalVals[x] for x in self.data.directionMarginal]
ARBs = self.interpolateMarginal(mValsM, self.data.ARBsList)
if self.data.nparPivot > 1:
mValsP[rDimP] = self.data.mu0Pivot(rDimP)
mValsP = self.centerNormalizePivot(mValsP)
mValsP = mValsP.data.flatten()
mValsP[rDimP] = fp
ARBs = marginalizePolyList(ARBs, mValsP, "auto")
ev = eigvalsNonlinearDense(ARBs, jSupp = jSupp, **kwargs)
return np.power(self.data.mu0(rDim) ** self.data.rescalingExp[rDim]
- + ev * self.scaleFactor[rDim],
+ + ev * self.data.scaleFactor[rDim],
1. / self.data.rescalingExp[rDim])
diff --git a/rrompy/reduction_methods/trained_model/trained_model_pole_matching_general.py b/rrompy/reduction_methods/trained_model/trained_model_pole_matching_general.py
index 1131771..dadd780 100644
--- a/rrompy/reduction_methods/trained_model/trained_model_pole_matching_general.py
+++ b/rrompy/reduction_methods/trained_model/trained_model_pole_matching_general.py
@@ -1,241 +1,241 @@
# Copyright (C) 2018 by the RROMPy authors
#
# This file is part of RROMPy.
#
# RROMPy is free software: you can redistribute it and/or modify
# it under the terms of the GNU Lesser General Public License as published by
# the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
#
# RROMPy is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU Lesser General Public License for more details.
#
# You should have received a copy of the GNU Lesser General Public License
# along with RROMPy. If not, see .
#
import numpy as np
from .trained_model_pivoted_general import TrainedModelPivotedGeneral
from rrompy.utilities.base.types import (Np1D, Tuple, ListAny, paramVal,
paramList, sampList, HFEng)
from rrompy.utilities.base import verbosityManager as vbMng, freepar as fp
from rrompy.utilities.numerical import pointMatching
from rrompy.utilities.poly_fitting.heaviside import HeavisideInterpolator as HI
from rrompy.utilities.exception_manager import RROMPyException, RROMPyAssert
from rrompy.parameter import checkParameter, checkParameterList
from rrompy.sampling import emptySampleList
__all__ = ['TrainedModelPoleMatchingGeneral']
class TrainedModelPoleMatchingGeneral(TrainedModelPivotedGeneral):
"""
ROM approximant evaluation for pivoted approximants with pole matching.
Attributes:
Data: dictionary with all that can be pickled.
"""
def initializeFromLists(self, poles:ListAny, coeffs:ListAny, basis:str,
HFEngine:HFEng, matchingWeight : float = 1.,
POD : bool = True):
"""Initialize Heaviside representation."""
musM = self.data.musMarginal
margAbsDist = np.sum(np.abs(np.repeat(musM.data, len(musM), 0)
- np.tile(musM.data, [len(musM), 1])
), axis = 1).reshape(len(musM), len(musM))
N = len(poles[0])
explored = [0]
unexplored = list(range(1, len(musM)))
for _ in range(1, len(musM)):
minIdx = np.argmin(np.concatenate([margAbsDist[ex, unexplored] \
for ex in explored]))
minIex = explored[minIdx // len(unexplored)]
minIunex = unexplored[minIdx % len(unexplored)]
dist = np.abs(np.tile(poles[minIex].reshape(-1, 1), N)
- poles[minIunex].reshape(1, -1))
if matchingWeight != 0:
resex = coeffs[minIex][: N]
resunex = coeffs[minIunex][: N]
if POD:
distR = resex.dot(resunex.T.conj())
distR = (distR.T / np.linalg.norm(resex, axis = 1)).T
distR = distR / np.linalg.norm(resunex, axis = 1)
else:
resex = self.data.projMat.dot(resex.T)
resunex = self.data.projMat.dot(resunex.T)
distR = HFEngine.innerProduct(resex, resunex).T
distR = (distR.T / HFEngine.norm(resex)).T
distR = distR / HFEngine.norm(resunex)
distR = np.abs(distR)
distR[distR > 1.] = 1.
dist += 2. / np.pi * matchingWeight * np.arccos(distR)
reordering = pointMatching(dist)
poles[minIunex] = poles[minIunex][reordering]
coeffs[minIunex][: N] = coeffs[minIunex][reordering]
explored += [minIunex]
unexplored.remove(minIunex)
HIs = []
for pls, cfs in zip(poles, coeffs):
hsi = HI()
hsi.poles = pls
hsi.coeffs = cfs
hsi.npar = 1
hsi.polybasis = basis
HIs += [hsi]
self.data.HIs = HIs
def recompressByCutOff(self, murange : Tuple[float, float] = [- 1., 1.],
tol : float = np.inf, rtype : str = "MAGNITUDE"):
if np.isinf(tol): return " No poles erased."
N = len(self.data.HIs[0].poles)
mu0 = np.mean(murange)
musig = murange[0] - mu0
if np.isclose(musig, 0.):
radius = lambda x: np.abs(x - mu0)
else:
if rtype == "MAGNITUDE":
murdir = (murange[0] - mu0) / np.abs(musig)
def radius(x):
scalprod = (x - mu0) * murdir.conj() / np.abs(musig)
rescalepar = np.abs(np.real(scalprod))
rescaleort = np.abs(np.imag(scalprod))
return ((rescalepar - 1.) ** 2. * (rescalepar > 1.)
+ rescaleort ** 2.) ** .5
else:#if rtype == "POTENTIAL":
def radius(x):
rescale = (x - mu0) / musig
return np.max(np.abs(rescale * np.array([-1., 1.])
+ (rescale ** 2. - 1) ** .5)) - 1.
keepPole, removePole = [], []
for j in range(N):
for hi in self.data.HIs:
if radius(hi.poles[j]) <= tol:
keepPole += [j]
break
if len(keepPole) == 0 or keepPole[-1] != j: removePole += [j]
if len(keepPole) == N: return " No poles erased."
keepCoeff = keepPole + [N]
keepCoeff = keepCoeff + list(range(N + 1,len(self.data.HIs[0].coeffs)))
for hi in self.data.HIs:
polyCorrection = np.zeros_like(hi.coeffs[0, :])
for j in removePole:
polyCorrection += hi.coeffs[j, :] / (mu0 - hi.poles[j])
if len(hi.coeffs) == N:
hi.coeffs = np.vstack((hi.coeffs, polyCorrection))
else:
hi.coeffs[N, :] += polyCorrection
hi.poles = hi.poles[keepPole]
hi.coeffs = hi.coeffs[keepCoeff, :]
return (" Erased {} pole".format(len(removePole))
+ "s" * (len(removePole) > 1) + ".")
def interpolateMarginalInterpolator(self, mu : paramVal = []) -> Np1D:
"""Obtain interpolated approximant interpolator."""
mu = checkParameter(mu, self.data.nparMarginal)[0]
hsi = HI()
hsi.poles = self.interpolateMarginalPoles(mu)
hsi.coeffs = self.interpolateMarginalCoeffs(mu)
hsi.npar = 1
hsi.polybasis = self.data.HIs[0].polybasis
return hsi
def interpolateMarginalPoles(self, mu : paramList = []) -> Np1D:
"""Obtain interpolated approximant poles."""
mu = checkParameterList(mu, self.data.nparMarginal)[0]
return self.interpolateMarginal(mu, [hi.poles for hi in self.data.HIs])
def interpolateMarginalCoeffs(self, mu : paramList = []) -> Np1D:
"""Obtain interpolated approximant coefficients."""
mu = checkParameterList(mu, self.data.nparMarginal)[0]
cs = self.interpolateMarginal(mu, [hi.coeffs for hi in self.data.HIs])
if isinstance(cs, (list, tuple,)): cs = np.array(cs)
return cs.reshape(self.data.HIs[0].coeffs.shape)
def getApproxReduced(self, mu : paramList = []) -> sampList:
"""
Evaluate reduced representation of approximant at arbitrary parameter.
Args:
mu: Target parameter.
"""
RROMPyAssert(self.data.nparPivot, 1, "Number of pivot parameters")
mu = checkParameterList(mu, self.data.npar)[0]
if (not hasattr(self, "lastSolvedApproxReduced")
or self.lastSolvedApproxReduced != mu):
vbMng(self, "INIT",
"Evaluating approximant at mu = {}.".format(mu), 12)
self.uApproxReduced = emptySampleList()
self.uApproxReduced.reset((self.data.HIs[0].coeffs.shape[1],
- len(mu)), self.data.projMat[0].dtype)
+ len(mu)), self.data.projMat.dtype)
for i, muPL in enumerate(mu):
muL = self.centerNormalizePivot([muPL(0, x) \
for x in self.data.directionPivot])
muM = [muPL(0, x) for x in self.data.directionMarginal]
vbMng(self, "INIT",
"Assembling reduced model for mu = {}.".format(muPL), 87)
hsL = self.interpolateMarginalInterpolator(muM)
vbMng(self, "DEL", "Done assembling reduced model.", 87)
self.uApproxReduced[i] = hsL(muL)
vbMng(self, "DEL", "Done evaluating approximant.", 12)
self.lastSolvedApproxReduced = mu
return self.uApproxReduced
def getPoles(self, *args, **kwargs) -> Np1D:
"""
Obtain approximant poles.
Returns:
Numpy complex vector of poles.
"""
RROMPyAssert(self.data.nparPivot, 1, "Number of pivot parameters")
if len(args) + len(kwargs) > 1:
raise RROMPyException(("Wrong number of parameters passed. "
"Only 1 available."))
elif len(args) + len(kwargs) == 1:
if len(args) == 1:
mVals = args[0]
else:
mVals = kwargs["marginalVals"]
if not hasattr(mVals, "__len__"): mVals = [mVals]
mVals = list(mVals)
else:
mVals = [fp]
try:
rDim = mVals.index(fp)
if rDim < len(mVals) - 1 and fp in mVals[rDim + 1 :]:
raise
except:
raise RROMPyException(("Exactly 1 'freepar' entry in "
"marginalVals must be provided."))
if rDim != self.data.directionPivot[0]:
raise RROMPyException(("'freepar' entry in marginalVals must "
"coincide with pivot direction."))
mVals[rDim] = self.data.mu0(rDim)
mMarg = [mVals[j] for j in range(len(mVals)) if j != rDim]
roots = np.array(self.interpolateMarginalPoles(mMarg))
return np.power(self.data.mu0(rDim) ** self.data.rescalingExp[rDim]
+ self.data.scaleFactor[rDim] * roots,
1. / self.data.rescalingExp[rDim])
def getResidues(self, *args, **kwargs) -> Np1D:
"""
Obtain approximant residues.
Returns:
Numpy matrix with residues as columns.
"""
pls = self.getPoles(*args, **kwargs)
if len(args) == 1:
mVals = args[0]
else:
mVals = kwargs["marginalVals"]
if not hasattr(mVals, "__len__"): mVals = [mVals]
mVals = list(mVals)
rDim = mVals.index(fp)
mMarg = [mVals[j] for j in range(len(mVals)) if j != rDim]
residues = self.interpolateMarginalCoeffs(mMarg)[: len(pls)]
- res = self.data.projMat.dot(residues)
+ res = self.data.projMat.dot(residues.T)
return pls, res
diff --git a/rrompy/utilities/poly_fitting/heaviside/heaviside_interpolator.py b/rrompy/utilities/poly_fitting/heaviside/heaviside_interpolator.py
index b95859c..818219e 100644
--- a/rrompy/utilities/poly_fitting/heaviside/heaviside_interpolator.py
+++ b/rrompy/utilities/poly_fitting/heaviside/heaviside_interpolator.py
@@ -1,76 +1,82 @@
# 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, Np1D, paramList,
interpEng)
from rrompy.utilities.base import freepar as fp
from rrompy.utilities.poly_fitting.polynomial.polynomial_interpolator import (
PolynomialInterpolator)
from rrompy.utilities.poly_fitting.polynomial.roots import polyroots
from .val import polyval
+from .heaviside_to_from_affine import affine2heaviside
from .heaviside_to_from_rational import heaviside2rational, rational2heaviside
from rrompy.utilities.exception_manager import RROMPyAssert
__all__ = ['HeavisideInterpolator']
class HeavisideInterpolator(PolynomialInterpolator):
"""HERE"""
def __init__(self, other = None):
if other is None: return
self.poles = other.poles
super().__init__(other)
def __call__(self, mu:paramList, der : List[int] = None,
scl : Np1D = None):
return polyval(mu, self.coeffs, self.poles, self.polybasis)
def __copy__(self):
return HeavisideInterpolator(self)
def __deepcopy__(self, memo):
other = HeavisideInterpolator()
other.poles, other.coeffs, other.npar, other.polybasis = copy(
(self.poles, self.coeffs, self.npar, self.polybasis), memo)
return other
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(nleft[0], 0, "Padding in free direction")
super().pad(nleft, nright)
+ def setupFromAffine(self, As:ListAny, bs:ListAny,
+ jSupp : int = 1):
+ self.coeffs, self.poles, self.polybasis = affine2heaviside(As, bs,
+ jSupp)
+
def setupFromRational(self, num:interpEng, den:interpEng,
murange : Np1D = np.array([-1., 1.]),
scl : Np1D = None, scalingExp : List[float] = None):
self.coeffs, self.poles, self.polybasis = rational2heaviside(num, den,
murange, scl,
scalingExp)
def roots(self, marginalVals : ListAny = [fp], murange : Np1D = None,
scalingExp : List[float] = None):
RROMPyAssert(self.shape, (1,), "Shape of output")
RROMPyAssert(marginalVals, [fp], "Marginal values")
basisN = self.polybasis.split("_")[0]
coeffsN = heaviside2rational(self.coeffs, self.poles, murange, basisN,
scalingExp = scalingExp)[0]
return polyroots(coeffsN, basisN)
diff --git a/rrompy/utilities/poly_fitting/polynomial/__init__.py b/rrompy/utilities/poly_fitting/polynomial/__init__.py
index 00b3f3d..8390263 100644
--- a/rrompy/utilities/poly_fitting/polynomial/__init__.py
+++ b/rrompy/utilities/poly_fitting/polynomial/__init__.py
@@ -1,47 +1,52 @@
# 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 .base import (polybases, polyfitname, polydomcoeff)
from .der import polyder
from .val import polyval
from .marginalize import polymarginalize
from .vander import polyvander, polyvanderTotal
from .roots import polyroots
from .derivative import nextDerivativeIndices
from .hash_derivative import hashDerivativeToIdx, hashIdxToDerivative
-from .total_degree import degreeTotalToFull
+from .degree import (fullDegreeSet, totalDegreeSet, degreeTotalToFull,
+ fullDegreeMaxMask, totalDegreeMaxMask)
from .polynomial_interpolator import PolynomialInterpolator
__all__ = [
'polybases',
'polyfitname',
'polydomcoeff',
'polyder',
'polyval',
'polymarginalize',
'polyvander',
'polyvanderTotal',
'polyroots',
'nextDerivativeIndices',
'hashDerivativeToIdx',
'hashIdxToDerivative',
+ 'fullDegreeSet',
+ 'totalDegreeSet',
'degreeTotalToFull',
+ 'fullDegreeMaxMask',
+ 'totalDegreeMaxMask',
'PolynomialInterpolator'
]
diff --git a/rrompy/utilities/poly_fitting/polynomial/degree.py b/rrompy/utilities/poly_fitting/polynomial/degree.py
new file mode 100644
index 0000000..631489c
--- /dev/null
+++ b/rrompy/utilities/poly_fitting/polynomial/degree.py
@@ -0,0 +1,53 @@
+# 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 rrompy.utilities.base.types import Np2D, Tuple, List
+from .hash_derivative import (hashDerivativeToIdx as hashD,
+ hashIdxToDerivative as hashI)
+from rrompy.utilities.numerical.kroneckerer import kroneckerer
+
+__all__ = ['fullDegreeSet', 'totalDegreeSet', 'degreeTotalToFull',
+ 'fullDegreeMaxMask', 'totalDegreeMaxMask']
+
+def fullDegreeSet(deg:int, npar:int) -> List[List[int]]:
+ idxs = np.empty(((deg + 1) ** npar, npar), dtype = int)
+ for j in range(npar):
+ idxs[:, j] = kroneckerer(np.arange(deg + 1), (deg + 1) ** j,
+ (deg + 1) ** (npar - j - 1))
+ return idxs
+
+def totalDegreeSet(deg:int, npar:int,
+ return_mask : bool = False) -> List[List[int]]:
+ remN = fullDegreeSet(deg, npar)
+ mask = np.sum(remN, axis = 1) <= deg
+ if return_mask: return remN[mask], mask
+ return remN[mask]
+
+def degreeTotalToFull(shapeFull:Tuple[int], dim:int, coeffs:Np2D) -> Np2D:
+ full = np.zeros(shapeFull, dtype = coeffs.dtype)
+ for j in range(len(coeffs)):
+ full[tuple(hashI(j, dim))] = coeffs[j]
+ return full
+
+def fullDegreeMaxMask(deg:int, npar:int) -> List[int]:
+ return np.where(np.any(fullDegreeSet(deg, npar) == deg, axis = 1))[0]
+
+def totalDegreeMaxMask(deg:int, npar:int) -> List[int]:
+ return np.arange(hashD([deg] + [0] * (npar - 1)),
+ hashD([1 + deg] + [0] * (npar - 1)))
diff --git a/rrompy/utilities/poly_fitting/polynomial/polynomial_interpolator.py b/rrompy/utilities/poly_fitting/polynomial/polynomial_interpolator.py
index 3346101..21f6de7 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 .total_degree import degreeTotalToFull
+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]
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/polynomial/total_degree.py b/rrompy/utilities/poly_fitting/polynomial/total_degree.py
deleted file mode 100644
index 3a0c1bb..0000000
--- a/rrompy/utilities/poly_fitting/polynomial/total_degree.py
+++ /dev/null
@@ -1,40 +0,0 @@
-# Copyright (C) 2018 by the RROMPy authors
-#
-# This file is part of RROMPy.
-#
-# RROMPy is free software: you can redistribute it and/or modify
-# it under the terms of the GNU Lesser General Public License as published by
-# the Free Software Foundation, either version 3 of the License, or
-# (at your option) any later version.
-#
-# RROMPy is distributed in the hope that it will be useful,
-# but WITHOUT ANY WARRANTY; without even the implied warranty of
-# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
-# GNU Lesser General Public License for more details.
-#
-# You should have received a copy of the GNU Lesser General Public License
-# along with RROMPy. If not, see .
-#
-
-import numpy as np
-from rrompy.utilities.base.types import Np2D, Tuple, List
-from .hash_derivative import hashIdxToDerivative as hashI
-
-__all__ = ['degreeTotalToFull']
-
-def totalDegreeMask(degs:List[int]) -> Tuple[List[int], List[bool]]:
- dim = len(degs)
- N = np.prod([d + 1 for d in degs])
- maskN = np.arange(N, dtype = int)
- remN = np.zeros((N, dim), dtype = int)
- for j in range(dim - 1, -1, -1):
- remN[:, j] = maskN % (degs[j] + 1)
- maskN //= degs[j] + 1
- mask = np.sum(remN, axis = 1) <= min(degs)
- return remN[mask], mask
-
-def degreeTotalToFull(shapeFull:Tuple[int], dim:int, coeffs:Np2D) -> Np2D:
- full = np.zeros(shapeFull, dtype = coeffs.dtype)
- for j in range(len(coeffs)):
- full[tuple(hashI(j, dim))] = coeffs[j]
- return full
diff --git a/rrompy/utilities/poly_fitting/polynomial/vander.py b/rrompy/utilities/poly_fitting/polynomial/vander.py
index e335241..40f195d 100644
--- a/rrompy/utilities/poly_fitting/polynomial/vander.py
+++ b/rrompy/utilities/poly_fitting/polynomial/vander.py
@@ -1,138 +1,138 @@
# 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 .der import polyder
-from .total_degree import totalDegreeMask
+from .degree import totalDegreeSet
from rrompy.utilities.base.types import Np1D, Np2D, Tuple, List, paramList
from rrompy.parameter import checkParameterList
from rrompy.utilities.exception_manager import RROMPyException, RROMPyAssert
__all__ = ['polyvander', 'polyvanderTotal']
def firstDerTransition(dim:int, TDirac:List[Np2D], basis:str,
scl : Np1D = None) -> Np2D:
C_m = np.zeros((dim, len(TDirac), len(TDirac)), dtype = float)
for j, Tj in enumerate(TDirac):
m, om = [0] * dim, [(0, 0)] * dim
for idx in range(dim):
m[idx], om[idx] = 1, (0, 1)
J_der = polyder(Tj, basis, m, scl)
C_m[idx, :, j] = np.ravel(np.pad(J_der, mode = "constant",
pad_width = om))
m[idx], om[idx] = 0, (0, 0)
return C_m
def countDerDirections(n:int, base:int, digits:int, idx:int):
if digits == 0: return []
dig = n % base
return [(idx, dig)] * (dig > 0) + countDerDirections(
(n - dig) // base, base, digits - 1, idx + 1)
def polyvander(x:paramList, degs:List[int], basis:str,
derIdxs : List[List[List[int]]] = None,
reorder : List[int] = None, scl : Np1D = None) -> Np2D:
"""
Compute full Hermite-Vandermonde matrix with specified derivative
directions.
E.g. assume that we want to obtain the Vandermonde matrix for
(value, derx, derx2) at x = [0, 0],
(value, dery) at x = [1, 0],
(dery, derxy) at x = [0, 0],
of degree 3 in x and 1 in y, using Chebyshev polynomials.
This can be done by
polyvander([[0, 0], [1, 0]], # unique sample points
[3, 1], # polynomial degree
"chebyshev", # polynomial family
[
[[0, 0], [1, 0], [2, 0], [0, 1], [1, 1]],
# derivative directions at first point
[[0, 0], [0, 1]] # derivative directions at second point
],
[0, 1, 2, 5, 6, 3, 4] # reorder indices
)
"""
if not isinstance(degs, (list,tuple,np.ndarray,)): degs = [degs]
dim = len(degs)
x = checkParameterList(x, dim)[0]
x_un, idx_un = x.unique(return_inverse = True)
if len(x_un) < len(x):
raise RROMPyException("Sample points must be distinct.")
del x_un
try:
vanderbase = {"CHEBYSHEV" : np.polynomial.chebyshev.chebvander,
"LEGENDRE" : np.polynomial.legendre.legvander,
"MONOMIAL" : np.polynomial.polynomial.polyvander
}[basis.upper()]
except:
raise RROMPyException("Polynomial basis not recognized.")
VanBase = vanderbase(x(0), degs[0])
for j in range(1, dim):
VNext = vanderbase(x(j), degs[j])
for jj in range(j): VNext = np.expand_dims(VNext, 1)
VanBase = VanBase[..., None] * VNext
VanBase = VanBase.reshape((len(x), -1))
if derIdxs is None or VanBase.shape[-1] <= 1:
Van = VanBase
else:
derFlat, idxRep = [], []
for j, derIdx in enumerate(derIdxs):
derFlat += derIdx[:]
idxRep += [j] * len(derIdx[:])
for j in range(len(derFlat)):
if not hasattr(derFlat[j], "__len__"):
derFlat[j] = [derFlat[j]]
RROMPyAssert(len(derFlat[j]), dim, "Number of dimensions")
TDirac = [y.reshape([d + 1 for d in degs])
for y in np.eye(VanBase.shape[-1], dtype = int)]
Cs_loc = firstDerTransition(dim, TDirac, basis, scl)
Van = np.empty((len(derFlat), VanBase.shape[-1]),
dtype = VanBase.dtype)
for j in range(len(derFlat)):
Van[j, :] = VanBase[idxRep[j], :]
for k in range(dim):
for der in range(derFlat[j][k]):
Van[j, :] = Van[j, :].dot(Cs_loc[k]) / (der + 1)
if reorder is not None: Van = Van[reorder, :]
return Van
def polyvanderTotal(x:paramList, deg:int, basis:str,
derIdxs : List[List[List[int]]] = None,
reorder : List[int] = None, scl : Np1D = None)\
-> Tuple[Np2D, List[List[int]], List[int]]:
"""
Compute full total degree Hermite-Vandermonde matrix with specified
derivative directions.
"""
x = checkParameterList(x)[0]
degs = [deg] * x.shape[1]
VanBase = polyvander(x, degs, basis, derIdxs, reorder, scl)
- derIdxs, mask = totalDegreeMask(degs)
+ derIdxs, mask = totalDegreeSet(deg, x.shape[1], return_mask = True)
ordIdxs = np.empty(len(derIdxs), dtype = int)
derTotal = np.array([np.sum(y) for y in derIdxs])
idxPrev = 0
rangeAux = np.arange(len(derIdxs))
for j in range(deg + 1):
idxLocal = rangeAux[derTotal == j][::-1]
idxPrev += len(idxLocal)
ordIdxs[idxPrev - len(idxLocal) : idxPrev] = idxLocal
return VanBase[:, mask], derIdxs, ordIdxs
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 e0a0be5..9f16561 100644
--- a/rrompy/utilities/poly_fitting/radial_basis/radial_basis_interpolator.py
+++ b/rrompy/utilities/poly_fitting/radial_basis/radial_basis_interpolator.py
@@ -1,138 +1,137 @@
# 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.total_degree import (
- degreeTotalToFull)
+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)
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
diff --git a/tests/reduction_methods_1D/rational_interpolant_1d.py b/tests/reduction_methods_1D/rational_interpolant_1d.py
index b8bd6ac..940d715 100644
--- a/tests/reduction_methods_1D/rational_interpolant_1d.py
+++ b/tests/reduction_methods_1D/rational_interpolant_1d.py
@@ -1,68 +1,68 @@
# 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 matrix_fft import matrixFFT
from rrompy.reduction_methods.standard import RationalInterpolant as RI
from rrompy.parameter.parameter_sampling import (QuadratureSampler as QS,
ManualSampler as MS)
from rrompy.parameter import checkParameterList
def test_monomials(capsys):
mu = 1.5
solver = matrixFFT()
params = {"POD": False, "M": 9, "N": 9, "S": 10, "robustTol": 1e-6,
"interpRcond": 1e-3, "polybasis": "MONOMIAL",
"sampler": QS([1.5, 6.5], "UNIFORM")}
approx = RI(solver, 4., params, verbosity = 0)
approx.setupApprox()
out, err = capsys.readouterr()
- assert "poorly conditioned. Reducing E " in out
+ assert "poorly conditioned. Reducing N " in out
assert len(err) == 0
- assert np.isclose(approx.normErr(mu)[0], 2.9051e-4, atol = 1e-4)
+ assert np.isclose(approx.normErr(mu)[0], 1.4746e-05, atol = 1e-4)
def test_well_cond():
mu = 1.5
solver = matrixFFT()
params = {"POD": True, "M": 9, "N": 9, "S": 10, "robustTol": 1e-14,
"interpRcond": 1e-10, "polybasis": "CHEBYSHEV",
"sampler": QS([1., 7.], "CHEBYSHEV")}
approx = RI(solver, 4., params, verbosity = 0)
approx.setupApprox()
poles = approx.getPoles()
for lambda_ in np.arange(1, 8):
assert np.isclose(np.min(np.abs(poles - lambda_)), 0., atol = 1e-4)
for mu in approx.mus:
assert np.isclose(approx.normErr(mu)[0], 0., atol = 1e-8)
def test_hermite():
mu = 1.5
solver = matrixFFT()
sampler0 = QS([1., 7.], "CHEBYSHEV")
points, _ = checkParameterList(np.tile(sampler0.generatePoints(4)(0), 3))
params = {"POD": True, "M": 11, "N": 11, "S": 12, "polybasis": "CHEBYSHEV",
"sampler": MS([1., 7.], points = points)}
approx = RI(solver, 4., params, verbosity = 0)
approx.setupApprox()
poles = approx.getPoles()
for lambda_ in np.arange(1, 8):
assert np.isclose(np.min(np.abs(poles - lambda_)), 0., atol = 1e-4)
for mu in approx.mus:
assert np.isclose(approx.normErr(mu)[0], 0., atol = 1e-8)
diff --git a/tests/reduction_methods_multiD/rational_interpolant_2d.py b/tests/reduction_methods_multiD/rational_interpolant_2d.py
index 7487a8c..4f6d449 100644
--- a/tests/reduction_methods_multiD/rational_interpolant_2d.py
+++ b/tests/reduction_methods_multiD/rational_interpolant_2d.py
@@ -1,79 +1,79 @@
# 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 matrix_random import matrixRandom
from rrompy.reduction_methods.standard import RationalInterpolant as RI
from rrompy.parameter.parameter_sampling import (QuadratureSampler as QS,
ManualSampler as MS)
def test_monomials(capsys):
mu = [5.05, 7.1]
mu0 = [5., 7.]
solver = matrixRandom()
uh = solver.solve(mu)[0]
params = {"POD": False, "M": 4, "N": 4, "S": [4, 4], "robustTol": 1e-6,
"interpRcond": 1e-3, "polybasis": "MONOMIAL",
"sampler": QS([[4.9, 6.85], [5.1, 7.15]], "UNIFORM")}
approx = RI(solver, mu0, params, verbosity = 0)
approx.setupApprox()
uhP1 = approx.getApprox(mu)[0]
errP = approx.getErr(mu)[0]
errNP = approx.normErr(mu)[0]
myerrP = uhP1 - uh
assert np.allclose(np.abs(errP - myerrP), 0., rtol = 1e-3)
assert np.isclose(solver.norm(errP), errNP, rtol = 1e-3)
resP = approx.getRes(mu)[0]
resNP = approx.normRes(mu)
assert np.isclose(solver.norm(resP), resNP, rtol = 1e-3)
assert np.allclose(np.abs(resP - (solver.b(mu) - solver.A(mu).dot(uhP1))),
0., rtol = 1e-3)
assert np.isclose(errNP / solver.norm(uh), 5.2667e-05, rtol = 1e-1)
out, err = capsys.readouterr()
print(out)
- assert ("poorly conditioned. Reducing E" in out)
+ assert ("poorly conditioned. Reducing N " in out)
assert len(err) == 0
def test_well_cond():
mu = [5.05, 7.1]
mu0 = [5., 7.]
solver = matrixRandom()
params = {"POD": True, "M": 3, "N": 3, "S": [4, 4],
"interpRcond": 1e-10, "polybasis": "CHEBYSHEV",
"sampler": QS([[4.9, 6.85], [5.1, 7.15]], "UNIFORM")}
approx = RI(solver, mu0, params, verbosity = 0)
approx.setupApprox()
print(approx.normErr(mu)[0] / approx.normHF(mu)[0])
assert np.isclose(approx.normErr(mu)[0] / approx.normHF(mu)[0],
5.98695e-05, rtol = 1e-1)
def test_hermite():
mu = [5.05, 7.1]
mu0 = [5., 7.]
solver = matrixRandom()
sampler0 = QS([[4.9, 6.85], [5.1, 7.15]], "UNIFORM")
params = {"POD": True, "M": 3, "N": 3, "S": [25], "polybasis": "CHEBYSHEV",
"sampler": MS([[4.9, 6.85], [5.1, 7.15]],
points = sampler0.generatePoints([3, 3]))}
approx = RI(solver, mu0, params, verbosity = 0)
approx.setupApprox()
assert np.isclose(approx.normErr(mu)[0] / approx.normHF(mu)[0],
- 0.0001485, rtol = 5e-1)
+ 0.000224, rtol = 5e-1)
diff --git a/tests/utilities/parameter_sampling.py b/tests/utilities/parameter_sampling.py
index b1ea088..46311b9 100644
--- a/tests/utilities/parameter_sampling.py
+++ b/tests/utilities/parameter_sampling.py
@@ -1,59 +1,59 @@
# 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 rrompy.parameter.parameter_sampling import (ManualSampler,
QuadratureSampler, RandomSampler, FFTSampler)
from rrompy.parameter import checkParameter
def test_manual():
sampler = ManualSampler(lims = [0., 3.], points = np.linspace(0, 3, 101),
scalingExp = 2.)
assert sampler.name() == "ManualSampler"
x = sampler.generatePoints(10)
assert np.allclose(x(0), np.linspace(0, 3, 101)[:10], rtol = 1e-5)
def test_quadrature():
sampler = QuadratureSampler(lims = [0., 3.], kind = "CHEBYSHEV")
- x = sampler.generatePoints(9)
- assert np.isclose(x(0)[2], 1.5, rtol = 1e-5)
+ x = sampler.generatePoints(9, reorder = False)
+ assert np.isclose(x(0)[4], 1.5, rtol = 1e-5)
def test_random():
sampler = RandomSampler(lims = [0., 3.], kind = "SOBOL", seed = 13432)
x = sampler.generatePoints(100)
assert np.isclose(x(0)[47], 0.55609130859375, rtol = 1e-5)
def test_fft():
sampler = FFTSampler(lims = [-1., 1.])
x = sampler.generatePoints(100)
assert np.allclose(np.power(x(0), 100), 1., rtol = 1e-5)
def test_2D():
sampler = QuadratureSampler(lims = [(0., 0.), (3., 1.)],
kind = "GAUSSLEGENDRE")
x = sampler.generatePoints((9, 5))
assert sum(np.isclose(x(0), 1.5)) == 5
assert sum(np.isclose(x(1), .5)) == 9
def test_4D():
sampler = RandomSampler(lims = [tuple([0.] * 4), tuple([1.] * 4)],
kind = "UNIFORM", seed = 1234)
x = sampler.generatePoints(10)
assert x.shape[1] == 4
assert checkParameter([x[0]]) == checkParameter([(0.191519450378892,
0.622108771039832, 0.437727739007115, 0.785358583713769)])