Homec4science
Diffusion RationalROMPy (master) Browse

RROMPy/examples/3d/pod/plot_zero_set_3.py
6643b38fd8aemaster

Last Change
Raw File

plot_zero_set_3.py
View Options

import warnings
from copy import deepcopy as copy
import numpy as np
from matplotlib import pyplot as plt
def plotZeroSet3(murange, murangeEff, approx, mu0, nSamples = 200,
marginal = [None, None, .05], exps = [2., 1., 1.],
clip = -1, polesImTol : float = 1e-5):
mNone = [i for i, m in enumerate(marginal) if m is None]
nNone = len(mNone)
if nNone not in [1, 2, 3]:
raise
if nNone == 1:
exp = exps[mNone[0]]
if hasattr(approx, "mus"):
mu2x = approx.mus(mNone[0]) ** exp
else:
mu2x = mu0[0] ** exp
murangeExp = [[murange[0][mNone[0]] ** exp],
[murange[1][mNone[0]] ** exp]]
mu1 = np.linspace(murangeEff[0][mNone[0]] ** exp,
murangeEff[1][mNone[0]] ** exp, nSamples)
mu1s = np.power(mu1, 1. / exp)
mus = []
mBase = copy(marginal)
for m1 in mu1s:
mBase[mNone[0]] = m1
mus += [copy(mBase)]
mu1 = np.real(mu1)
Z = approx.trainedModel.getQVal(mus)
Zabs = np.abs(Z)
Zmin, Zmax = np.min(Zabs), np.max(Zabs)
plt.figure(figsize = (15, 7))
plt.jet()
plt.semilogy(mu1, Zabs)
for l_ in approx.trainedModel.getPoles(marginal):
plt.plot([np.real(l_ ** exp)] * 2, [Zmin, Zmax], 'b--')
plt.plot(mu2x, [Zmax] * len(mu2x), 'kx')
plt.plot([murangeExp[0][0]] * 2, [Zmin, Zmax], 'm:')
plt.plot([murangeExp[1][0]] * 2, [Zmin, Zmax], 'm:')
plt.xlim(mu1[0], mu1[-1])
plt.title("|Q(mu)|")
plt.grid()
plt.show()
return mus, Z
if nNone == 2:
exps = [exps[m] for m in mNone]
if hasattr(approx, "mus"):
mu2x = approx.mus(mNone[0]) ** exps[0]
mu2y = approx.mus(mNone[1]) ** exps[1]
else:
mu2x, mu2y = mu0[mNone[0]] ** exps[0], mu0[mNone[1]] ** exps[1]
murangeExp = [[murange[0][mNone[0]] ** exps[0],
murange[0][mNone[1]] ** exps[1]],
[murange[1][mNone[0]] ** exps[0],
murange[1][mNone[1]] ** exps[1]]]
mu1 = np.linspace(murangeEff[0][mNone[0]] ** exps[0],
murangeEff[1][mNone[0]] ** exps[0], nSamples)
mu2 = np.linspace(murangeEff[0][mNone[1]] ** exps[1],
murangeEff[1][mNone[1]] ** exps[1], nSamples)
mu1s = np.power(mu1, 1. / exps[0])
mu2s = np.power(mu2, 1. / exps[1])
Mu1, Mu2 = np.meshgrid(np.real(mu1), np.real(mu2))
mus = []
mBase = copy(marginal)
for m2 in mu2s:
mBase[mNone[1]] = m2
for m1 in mu1s:
mBase[mNone[0]] = m1
mus += [copy(mBase)]
poles1, poles2 = [], []
for m in mus:
m[mNone[0]] = None
pls = approx.getPoles(m) ** exps[0]
pls[np.imag(pls) > polesImTol] = np.nan
pls = np.real(pls)
poles1 += list(pls)
poles2 += [m[mNone[1]] ** exps[1]] * len(pls)
try:
Z = approx.trainedModel.getQVal(mus).reshape(Mu1.shape)
Zabs = np.log10(np.abs(Z))
Zabsmin, Zabsmax = np.min(Zabs), np.max(Zabs)
if clip > 0:
Zabsmin += clip * (Zabsmax - Zabsmin)
cmap = plt.cm.bone_r
else:
cmap = plt.cm.jet
except:
Z = None
cmap = plt.cm.jet
warnings.simplefilter("ignore", category = (UserWarning,
np.ComplexWarning))
plt.figure(figsize = (15, 7))
plt.jet()
if Z is not None:
p = plt.contourf(Mu1, Mu2, Zabs, cmap = cmap,
levels = np.linspace(Zabsmin, Zabsmax, 50))
if clip > 0:
plt.contour(Mu1, Mu2, Zabs, [Zabsmin])
plt.plot(poles1, poles2, 'k.')
plt.plot(mu2x, mu2y, 'kx')
plt.plot([murangeExp[0][0]] * 2,
[murangeExp[0][1], murangeExp[1][1]], 'm:')
plt.plot([murangeExp[0][0], murangeExp[1][0]],
[murangeExp[1][1]] * 2, 'm:')
plt.plot([murangeExp[1][0]] * 2,
[murangeExp[1][1], murangeExp[0][1]], 'm:')
plt.plot([murangeExp[1][0], murangeExp[0][0]],
[murangeExp[0][1]] * 2, 'm:')
if Z is not None:
plt.colorbar(p)
plt.xlim(murangeExp[0][0], murangeExp[1][0])
plt.ylim(murangeExp[0][1], murangeExp[1][1])
plt.title("log10|Q(mu)|")
plt.grid()
plt.show()
return mus, Z
if nNone == 3:
exps = [exps[m] for m in mNone]
murangeExp = [[murange[0][mNone[0]] ** exps[0],
murange[0][mNone[1]] ** exps[1],
murange[0][mNone[2]] ** exps[2]],
[murange[1][mNone[0]] ** exps[0],
murange[1][mNone[1]] ** exps[1],
murange[1][mNone[2]] ** exps[2]]]
mu1 = np.linspace(murangeEff[0][mNone[0]] ** exps[0],
murangeEff[1][mNone[0]] ** exps[0], nSamples)
mu2 = np.linspace(murangeEff[0][mNone[1]] ** exps[1],
murangeEff[1][mNone[1]] ** exps[1], nSamples)
mu3 = np.linspace(murangeEff[0][mNone[2]] ** exps[2],
murangeEff[1][mNone[2]] ** exps[2], nSamples)
mu1s = np.power(mu1, 1. / exps[0])
mu2s = np.power(mu2, 1. / exps[1])
mu3s = np.power(mu3, 1. / exps[2])
Mu1, Mu2, Mu3 = np.meshgrid(np.real(mu1), np.real(mu2), np.real(mu3),
indexing = 'ij')
mus = []
mBase = copy(marginal)
for m2 in mu2s:
mBase[mNone[1]] = m2
for m3 in mu3s:
mBase[mNone[2]] = m3
mus += [copy(mBase)]
poles1, poles2, poles3 = [], [], []
for m in mus:
m[mNone[0]] = None
pls = approx.getPoles(m) ** exps[0]
pls[np.imag(pls) > polesImTol] = np.nan
pls = np.real(pls)
pls[pls < murangeExp[0][0]] = np.nan
pls[pls > murangeExp[1][0]] = np.nan
poles1 += list(pls)
poles2 += [m[mNone[1]] ** exps[1]] * len(pls)
poles3 += [m[mNone[2]] ** exps[2]] * len(pls)
pts = np.empty((len(poles1), 3))
pts[:, 0] = poles1
pts[:, 1] = poles2
pts[:, 2] = poles3
if np.size(pts) > 0:
from mpl_toolkits.mplot3d import Axes3D
fig = plt.figure(figsize = (8, 6))
ax = Axes3D(fig)
ax.scatter(pts[:, 0], pts[:, 1], pts[:, 2], '.', s = 1)
ax.set_xlim(murangeExp[0][0], murangeExp[1][0])
ax.set_ylim(murangeExp[0][1], murangeExp[1][1])
ax.set_zlim(murangeExp[0][2], murangeExp[1][2])
plt.show()
plt.close()
return pts
c4science · Help