Skip to content

Instantly share code, notes, and snippets.

@liweitianux
Created September 12, 2018 04:13
Show Gist options
  • Save liweitianux/29d90a9f13967f297a4851805e678931 to your computer and use it in GitHub Desktop.
Save liweitianux/29d90a9f13967f297a4851805e678931 to your computer and use it in GitHub Desktop.
Matplotlib example
import matplotlib as mpl
import matplotlib.pyplot as plt
# %matplotlib inline
# Matplotlib Customization
# References: https://matplotlib.org/users/customizing.html
# refresh the font cache:
#mpl.font_manager._rebuild()
# Reset to the defaults
#matplotlib.rcdefaults()
mpl.style.use("ggplot")
for k, v in [("font.family", "Inconsolata"),
("font.size", 14.0),
("pdf.fonttype", 42), # Type 42 (a.k.a. TrueType)
("figure.figsize", [8, 6]),
("image.cmap", "jet"),
("xtick.labelsize", "large"),
("xtick.major.size", 7.0),
("xtick.major.width", 2.0),
("xtick.minor.size", 4.0),
("xtick.minor.width", 1.5),
("ytick.labelsize", "large"),
("ytick.major.size", 7.0),
("ytick.major.width", 2.0),
("ytick.minor.size", 4.0),
("ytick.minor.width", 1.5)]:
mpl.rcParams[k] = v
def plot_fluxfunc(data, nbin=50, q=None, figsize=(6, 6), figfile=None):
fidx158 = data['fidx158']
fidx1400 = data['fidx1400']
halos_df = data['halos_df']
halos = data['halos']
area_asec2 = data['area_asec2']
allskyfactor = data['allskyfactor']
testdirs = data['testdirs']
ntests = data['ntests']
flux158 = [None]*ntests
flux1400 = [None]*ntests
power158 = [None]*ntests
power1400 = [None]*ntests
for i, df in enumerate(halos):
power158[i] = np.array(df["power[%d]" % fidx158]) / 1e24 # [1e24 W/Hz]
flux158[i] = 1e3 * np.array(df["flux[%d]" % fidx158]) # [mJy]
power1400[i] = np.array(df["power[%d]" % fidx1400]) / 1e24 # [1e24 W/Hz]
flux1400[i] = 1e3 * np.array(df["flux[%d]" % fidx1400]) # [mJy]
Tb158 = [Fnu_to_Tb(flux.sum(), area_asec2, 158) for flux in flux158] # [mK]
Tb1400 = [Fnu_to_Tb(flux.sum(), area_asec2, 1400) for flux in flux1400]
fig, ax = plt.subplots(figsize=figsize)
# observed halos
cobs1400, xobs1400, __ = countdist_integrated(obs2017["S1.4[mJy]"], nbin=20)
ax.loglog(xobs1400, cobs1400, color="C2", ls="--", lw=2.5, label="Observations at 1400 MHz")
try:
qlow, qhigh = q
except TypeError:
qlow, qhigh = 0, 100
vlow, vhigh = np.percentile(Tb158, q=(qlow, qhigh))
idx_keep = [idx for idx, Tb in enumerate(Tb158) if (Tb>=vlow and Tb<=vhigh)]
Tb158_ = [Tb for idx, Tb in enumerate(Tb158) if idx in idx_keep]
Tb1400_ = [Tb for idx, Tb in enumerate(Tb1400) if idx in idx_keep]
power158_ = [power for idx, power in enumerate(power158) if idx in idx_keep]
power1400_ = [power for idx, power in enumerate(power1400) if idx in idx_keep]
flux158_ = [flux for idx, flux in enumerate(flux158) if idx in idx_keep]
flux1400_ = [flux for idx, flux in enumerate(flux1400) if idx in idx_keep]
print("----------------------------------")
print("q: %d-%d; #tests: %d" % (qlow, qhigh, len(flux158_)))
print("Tb@158: %.3f +/- %.3f [mK]" % (np.mean(Tb158_), np.std(Tb158_)))
print("Tb@158: %.3f (%.3f - %.3f) [mK]" % tuple(np.percentile(Tb158_, q=(50, 25, 75))))
print("Tb@1400: %.5f +/- %.5f [mK]" % (np.mean(Tb1400_), np.std(Tb1400_)))
print("Tb@1400: %.5f (%.5f - %.5f) [mK]" % tuple(np.percentile(Tb1400_, q=(50, 25, 75))))
flux = np.concatenate(flux158_)
f158min, f158max = flux.min(), flux.max()
flux = np.concatenate(flux1400_)
f1400min, f1400max = flux.min(), flux.max()
ntests = len(flux158_)
ff158 = np.zeros(shape=(ntests, nbin))
ff1400 = np.zeros(shape=(ntests, nbin))
pp158 = np.zeros(shape=(ntests, nbin))
pp1400 = np.zeros(shape=(ntests, nbin))
for i in range(ntests):
# 1400 [MHz]
flux = flux1400_[i]
c1400, x1400, __ = countdist_integrated(flux, nbin=nbin, xmin=f1400min, xmax=f1400max)
ff1400[i, :] = c1400
#ax.loglog(x1400, c1400*allskyfactor, alpha=0.1, color="C0")
# 158 [MHz]
flux = flux158_[i]
c158, x158, __ = countdist_integrated(flux, nbin=nbin, xmin=f158min, xmax=f158max)
ff158[i, :] = c158
#ax.loglog(x158, c158*allskyfactor, alpha=0.1, color="C1")
# mean flux function
ff1400mean = ff1400.mean(axis=0)
ff158mean = ff158.mean(axis=0)
ax.loglog(x1400, ff1400mean*allskyfactor, color="C2", ls="-", lw=2, label="Simulations at 1400 MHz")
ax.loglog(x158, ff158mean*allskyfactor, color="C1", ls="-", lw=2, label="Simulations at 158 MHz")
# uncertainty ranges
ff1400std = ff1400.std(axis=0)
ff158std = ff158.std(axis=0)
ax.fill_between(x1400, (ff1400mean-ff1400std)*allskyfactor,
(ff1400mean+ff1400std)*allskyfactor, alpha=0.3, color="C2")
ax.fill_between(x158, (ff158mean-ff158std)*allskyfactor,
(ff158mean+ff158std)*allskyfactor, alpha=0.3, color="C1")
ax.set(xscale="log", yscale="log",
xlim=(1e-3, 3e3), ylim=(2, 1e5),
xlabel="Flux Density [mJy]",
ylabel="Integrated Counts (>flux)")
ax.legend()
if figfile:
fig.savefig(figfile, dpi=150, bbox_inches='tight')
print("saved figure to: %s" % figfile)
figfile = path.splitext(figfile)[0] + '.pdf'
fig.savefig(figfile, bbox_inches='tight')
print("saved figure to: %s" % figfile)
else:
plt.show()
plot_fluxfunc(data, q=(2.5, 92.5), figsize=(8, 8), figfile='fluxfunc-simucomp-1400.png')
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment