Exploring the Impact of Mini-Halos¶
[17]:
import matplotlib
import numpy as np
import attrs
import matplotlib.pyplot as plt
import py21cmfast as p21c
from tempfile import mkdtemp
Basic Setup¶
Let’s fix the initial condition for this tutorial.
[3]:
output_dir = "cache"
cache = p21c.OutputCache(output_dir)
[5]:
# USE_FFTW_WISDOM makes the FFT faster
inputs = p21c.InputParameters.from_template(
'Qin20', HII_DIM=128, LOWRES_CELL_SIZE_MPC=2.0, HIRES_TO_LOWRES_FACTOR=2, USE_FFTW_WISDOM=True, N_THREADS=4, USE_RELATIVE_VELOCITIES=True, random_seed=1993
)
initial_conditions = p21c.compute_initial_conditions(
inputs=inputs, cache=cache, write=True
)output_dir
Fiducial Lightcone¶
Let’s run a ‘fiducial’ model and see its lightcones
Note that the reference model has
pow(10, "F_STAR7_MINI") = pow(10, "F_STAR10") / pow(1000,ALPHA_STAR) * 10 # 10 times enhancement
pow(10, "F_ESC7_MINI" ) = pow(10, "F_ESC10" ) / pow(1000,ALPHA_ESC ) / 10 # 0.1 times enhancement to balance the 10 times enhanced Ngamma
pow(10, "L_X_MINI" ) = pow(10, "L_X")
1 - "F_H2_SHIELD" = 1
[9]:
# the lightcones we want to plot later together with their color maps and min/max
lightcone_quantities = (
"brightness_temp",
"spin_temperature",
"neutral_fraction",
"cumulative_recombinations",
"z_reion",
"ionisation_rate_G12",
"J_21_LW",
"density",
)
global_quantities = (
"brightness_temp",
"spin_temperature",
"neutral_fraction",
"cumulative_recombinations",
"z_reion",
"ionisation_rate_G12",
"density",
)
lcn = p21c.RectilinearLightconer.between_redshifts(
min_redshift=min(inputs.node_redshifts)+0.1,
max_redshift=max(inputs.node_redshifts)-0.1,
quantities=lightcone_quantities,
resolution=inputs.simulation_options.cell_size,
)
# This takes about 40 minutes.
lightcone_fid = p21c.run_lightcone(
lightconer=lcn,
inputs=inputs,
initial_conditions=initial_conditions,
cache=cache,
write=True,
progressbar=True
)
lightcone_fid.save(cache.direc / "mini-halo-fiducial-lightcone.h5")
[10]:
plt.style.use("dark_background")
EoR_colour = matplotlib.colors.LinearSegmentedColormap.from_list(
"mycmap",
[
(0, "white"),
(0.33, "yellow"),
(0.5, "orange"),
(0.68, "red"),
(0.83333, "black"),
(0.9, "blue"),
(1, "cyan"),
],
)
cmaps = [
EoR_colour,
"Reds",
"magma",
"magma",
"magma",
"cubehelix",
"cubehelix",
"viridis",
]
vmins = [-150, 1e1, 0, 0, 5, 0, 0, -1]
vmaxs = [30, 1e3, 1, 2, 9, 1, 10, 1]
bt = lightcone_fid.lightcones['brightness_temp']
fig, axs = plt.subplots(
len(lightcone_quantities),
1,
figsize=(
12,
12 * (bt.shape[0]/bt.shape[2]) * len(lightcone_quantities),
),
sharex=True,
sharey=True,
layout='tight',
gridspec_kw={"hspace": 0, "wspace": 0}
)
for ii, lightcone_quantity in enumerate(lightcone_quantities):
axs[ii].imshow(
lightcone_fid.lightcones[lightcone_quantity][1],
vmin=vmins[ii],
vmax=vmaxs[ii],
cmap=cmaps[ii],
)
t = axs[ii].text(
0.99,
0.1,
lightcone_quantity,
horizontalalignment="right",
verticalalignment="bottom",
transform=axs[ii].transAxes,
color="w" if lightcone_quantity in ['brightness_temp', 'density'] else 'C0',
#backgroundcolor="white",
fontsize=15,
)
# t.set_bbox(dict(facecolor='white', alpha=0.5, edgecolor='none'))
axs[ii].set_xticks([])
axs[ii].set_yticks([])
Varying Parameters¶
Let’s vary paremeters that describe mini-halos and see the impact on the global signal.
We keep other parameters fixed and vary one of following by a factor of 0.1, 0.5, 2 and 10:
pow(10, "F_STAR7_MINI")
pow(10, "F_ESC7_MINI")
pow(10, "L_X_MINI")
1 - "F_H2_SHIELD"
We also have a NO-mini model where mini-halos are not included
[11]:
# defining those color, linstyle, b
linestyles = ["-", "-", ":", "-.", "-.", ":"]
colors = ["gray", "black", "#e41a1c", "#377eb8", "#e41a1c", "#377eb8"]
lws = [1, 3, 2, 2, 2, 2]
texts = [
"varying " + r"$f_{*,7}^{\rm mol}$",
"varying " + r"$f_{\rm esc}^{\rm mol}$",
"varying " + r"$L_{\rm x}^{\rm mol}$",
"varying " + r"$1-f_{\rm H_2}^{\rm shield}$",
]
var_params = ["F_STAR7_MINI","F_ESC7_MINI","L_X_MINI","F_H2_SHIELD"]
factors = [0, 1, 0.1, 0.5, 2, 10]
labels = ["NO-mini", "reference", "x0.1", "x0.5", "x2", "x10"]
Note that I’ve run these simulations in parallel before this tutorial. With these setup, each took ~6h to finish. Here, running means read the cached outputs.
[18]:
all_lightcones = {}
for pp, text in enumerate(texts):
all_lightcones[text] = {}
for kk, label in enumerate(labels):
print(text, label)
factor = factors[kk]
if label!='NO-mini':
gq = global_quantities + ("J_21_LW",)
_lcn = lcn
else:
gq = global_quantities
_lcn = attrs.evolve(lcn, quantities=[q for q in lcn.quantities if q != 'J_21_LW'])
if var_params[pp] == "F_H2_SHIELD":
if factor > 1:
continue
inputs_var = inputs.evolve_input_structs(
USE_MINI_HALOS=(label!="NO-mini"),
F_H2_SHIELD=1. - (1. - inputs.astro_params.F_H2_SHIELD)*factor,
)
else:
inputs_var = inputs.evolve_input_structs(
USE_MINI_HALOS=(label!="NO-mini"),
**{var_params[pp] : getattr(inputs.astro_params, var_params[pp]) * factor}
)
if label == "reference":
lightcone = lightcone_fid
else:
lightcone = p21c.run_lightcone(
lightconer=_lcn,
initial_conditions=initial_conditions,
inputs=inputs_var,
global_quantities=gq,
cache=cache,
write=True,
progressbar=True
)
all_lightcones[text][label] = lightcone
Evolving Astrophysics ━━━━━━━━━━╸━━━━━━━━━━━━━━━━━━━━━━━━━━ 30% • 26/88 redshifts • 0:05:41 • 0:23:17 remaining
IOPub message rate exceeded.
The Jupyter server will temporarily stop sending output
to the client in order to avoid crashing it.
To change this limit, set the config variable
`--ServerApp.iopub_msg_rate_limit`.
Current values:
ServerApp.iopub_msg_rate_limit=1000.0 (msgs/sec)
ServerApp.rate_limit_window=3.0 (secs)
IOPub message rate exceeded.
The Jupyter server will temporarily stop sending output
to the client in order to avoid crashing it.
To change this limit, set the config variable
`--ServerApp.iopub_msg_rate_limit`.
Current values:
ServerApp.iopub_msg_rate_limit=1000.0 (msgs/sec)
ServerApp.rate_limit_window=3.0 (secs)
IOPub message rate exceeded.
The Jupyter server will temporarily stop sending output
to the client in order to avoid crashing it.
To change this limit, set the config variable
`--ServerApp.iopub_msg_rate_limit`.
Current values:
ServerApp.iopub_msg_rate_limit=1000.0 (msgs/sec)
ServerApp.rate_limit_window=3.0 (secs)
Evolving Astrophysics ━━━━━━━━━━━━━━━━━━━━━━╺━━━━━━━━━━━━━━ 60% • 53/88 redshifts • 0:19:11 • 0:17:27 remaining
IOPub message rate exceeded.
The Jupyter server will temporarily stop sending output
to the client in order to avoid crashing it.
To change this limit, set the config variable
`--ServerApp.iopub_msg_rate_limit`.
Current values:
ServerApp.iopub_msg_rate_limit=1000.0 (msgs/sec)
ServerApp.rate_limit_window=3.0 (secs)
Global Averages¶
[27]:
all_lightcones.keys()
[27]:
dict_keys(['varying $f_{*,7}^{\\rm mol}$', 'varying $f_{\\rm esc}^{\\rm mol}$', 'varying $L_{\\rm x}^{\\rm mol}$'])
[28]:
lightcone_fid.global_quantities.keys()
[28]:
dict_keys(['brightness_temp', 'neutral_fraction'])
[26]:
# choose some to plot...
plot_quantities = (
"brightness_temp",
"spin_temperature",
"neutral_fraction",
"cumulative_recombinations",
"ionisation_rate_G12",
"J_21_LW",
)
ymins = [-120, 1e1, 0, 0, 0, 0]
ymaxs = [30, 1e3, 1, 1, 1, 10]
fig, axss = plt.subplots(
len(plot_quantities),
len(labels),
sharex=True,
figsize=(4 * len(labels), 2 * len(global_quantities)),
)
for pp, text in enumerate(texts):
axs = axss[:, pp]
for kk, label in enumerate(labels):
factor = factors[kk]
lightcone = all_lightcones[text][label]
freqs = 1420.4 / (np.array(lightcone.inputs.node_redshifts) + 1.0)
for jj, global_quantity in enumerate(plot_quantities):
if global_quantity in lightcone.global_quantities:
axs[jj].plot(
freqs,
lightcone.global_quantities[global_quantity],
color=colors[kk],
linestyle=linestyles[kk],
label=labels[kk],
lw=lws[kk],
)
else:
print(global_quantity)
axs[0].text(
0.01,
0.99,
texts,
horizontalalignment="left",
verticalalignment="top",
transform=axs[0].transAxes,
fontsize=15,
)
for jj, global_quantity in enumerate(plot_quantities):
axs[jj].set_ylim(ymins[jj], ymaxs[jj])
axs[-1].set_xlabel("Frequency/MHz", fontsize=15)
axs[-1].xaxis.set_tick_params(labelsize=15)
axs[0].set_xlim(1420.4 / (35 + 1.0), 1420.4 / (5.5 + 1.0))
zlabels = np.array([6, 7, 8, 10, 13, 18, 25, 35])
ax2 = axs[0].twiny()
ax2.set_xlim(axs[0].get_xlim())
ax2.set_xticks(1420.4 / (zlabels + 1.0))
ax2.set_xticklabels(zlabels.astype(str))
ax2.set_xlabel("redshift", fontsize=15)
ax2.xaxis.set_tick_params(labelsize=15)
ax2.grid(False)
if pp == 0:
axs[0].legend(
loc="lower right", ncol=2, fontsize=13, fancybox=True, frameon=True
)
for jj, global_quantity in enumerate(plot_quantities):
axs[jj].set_ylabel(
"global_{}".format(global_quantity.replace("_box", "")), fontsize=15
)
axs[jj].yaxis.set_tick_params(labelsize=15)
else:
for jj, global_quantity in enumerate(plot_quantities):
axs[jj].set_ylabel(
"global_{}".format(global_quantity.replace("_box", "")), fontsize=0
)
axs[jj].yaxis.set_tick_params(labelsize=0)
plt.tight_layout()
fig.subplots_adjust(hspace=0.0, wspace=0.0)
J_21_LW
spin_temperature
cumulative_recombinations
ionisation_rate_G12
J_21_LW
J_21_LW
spin_temperature
cumulative_recombinations
ionisation_rate_G12
J_21_LW
---------------------------------------------------------------------------
KeyError Traceback (most recent call last)
Cell In [26], line 25
22 for kk, label in enumerate(labels):
23 factor = factors[kk]
---> 25 lightcone = all_lightcones[text][label]
27 freqs = 1420.4 / (np.array(lightcone.inputs.node_redshifts) + 1.0)
28 for jj, global_quantity in enumerate(plot_quantities):
KeyError: 'NO-mini'
Error in callback <function _draw_all_if_interactive at 0x7f694b6efd80> (for post_execute):
---------------------------------------------------------------------------
ValueError Traceback (most recent call last)
File /data4/smurray/miniforge3/envs/21cmfast/lib/python3.11/site-packages/matplotlib/pyplot.py:279, in _draw_all_if_interactive()
277 def _draw_all_if_interactive() -> None:
278 if matplotlib.is_interactive():
--> 279 draw_all()
File /data4/smurray/miniforge3/envs/21cmfast/lib/python3.11/site-packages/matplotlib/_pylab_helpers.py:131, in Gcf.draw_all(cls, force)
129 for manager in cls.get_all_fig_managers():
130 if force or manager.canvas.figure.stale:
--> 131 manager.canvas.draw_idle()
File /data4/smurray/miniforge3/envs/21cmfast/lib/python3.11/site-packages/matplotlib/backend_bases.py:1891, in FigureCanvasBase.draw_idle(self, *args, **kwargs)
1889 if not self._is_idle_drawing:
1890 with self._idle_draw_cntx():
-> 1891 self.draw(*args, **kwargs)
File /data4/smurray/miniforge3/envs/21cmfast/lib/python3.11/site-packages/matplotlib/backends/backend_agg.py:382, in FigureCanvasAgg.draw(self)
379 # Acquire a lock on the shared font cache.
380 with (self.toolbar._wait_cursor_for_draw_cm() if self.toolbar
381 else nullcontext()):
--> 382 self.figure.draw(self.renderer)
383 # A GUI class may be need to update a window using this draw, so
384 # don't forget to call the superclass.
385 super().draw()
File /data4/smurray/miniforge3/envs/21cmfast/lib/python3.11/site-packages/matplotlib/artist.py:94, in _finalize_rasterization.<locals>.draw_wrapper(artist, renderer, *args, **kwargs)
92 @wraps(draw)
93 def draw_wrapper(artist, renderer, *args, **kwargs):
---> 94 result = draw(artist, renderer, *args, **kwargs)
95 if renderer._rasterizing:
96 renderer.stop_rasterizing()
File /data4/smurray/miniforge3/envs/21cmfast/lib/python3.11/site-packages/matplotlib/artist.py:71, in allow_rasterization.<locals>.draw_wrapper(artist, renderer)
68 if artist.get_agg_filter() is not None:
69 renderer.start_filter()
---> 71 return draw(artist, renderer)
72 finally:
73 if artist.get_agg_filter() is not None:
File /data4/smurray/miniforge3/envs/21cmfast/lib/python3.11/site-packages/matplotlib/figure.py:3257, in Figure.draw(self, renderer)
3254 # ValueError can occur when resizing a window.
3256 self.patch.draw(renderer)
-> 3257 mimage._draw_list_compositing_images(
3258 renderer, self, artists, self.suppressComposite)
3260 renderer.close_group('figure')
3261 finally:
File /data4/smurray/miniforge3/envs/21cmfast/lib/python3.11/site-packages/matplotlib/image.py:134, in _draw_list_compositing_images(renderer, parent, artists, suppress_composite)
132 if not_composite or not has_images:
133 for a in artists:
--> 134 a.draw(renderer)
135 else:
136 # Composite any adjacent images together
137 image_group = []
File /data4/smurray/miniforge3/envs/21cmfast/lib/python3.11/site-packages/matplotlib/artist.py:71, in allow_rasterization.<locals>.draw_wrapper(artist, renderer)
68 if artist.get_agg_filter() is not None:
69 renderer.start_filter()
---> 71 return draw(artist, renderer)
72 finally:
73 if artist.get_agg_filter() is not None:
File /data4/smurray/miniforge3/envs/21cmfast/lib/python3.11/site-packages/matplotlib/axes/_base.py:3210, in _AxesBase.draw(self, renderer)
3207 if artists_rasterized:
3208 _draw_rasterized(self.get_figure(root=True), artists_rasterized, renderer)
-> 3210 mimage._draw_list_compositing_images(
3211 renderer, self, artists, self.get_figure(root=True).suppressComposite)
3213 renderer.close_group('axes')
3214 self.stale = False
File /data4/smurray/miniforge3/envs/21cmfast/lib/python3.11/site-packages/matplotlib/image.py:134, in _draw_list_compositing_images(renderer, parent, artists, suppress_composite)
132 if not_composite or not has_images:
133 for a in artists:
--> 134 a.draw(renderer)
135 else:
136 # Composite any adjacent images together
137 image_group = []
File /data4/smurray/miniforge3/envs/21cmfast/lib/python3.11/site-packages/matplotlib/artist.py:71, in allow_rasterization.<locals>.draw_wrapper(artist, renderer)
68 if artist.get_agg_filter() is not None:
69 renderer.start_filter()
---> 71 return draw(artist, renderer)
72 finally:
73 if artist.get_agg_filter() is not None:
File /data4/smurray/miniforge3/envs/21cmfast/lib/python3.11/site-packages/matplotlib/text.py:752, in Text.draw(self, renderer)
749 renderer.open_group('text', self.get_gid())
751 with self._cm_set(text=self._get_wrapped_text()):
--> 752 bbox, info, descent = self._get_layout(renderer)
753 trans = self.get_transform()
755 # don't use self.get_position here, which refers to text
756 # position in Text:
File /data4/smurray/miniforge3/envs/21cmfast/lib/python3.11/site-packages/matplotlib/text.py:382, in Text._get_layout(self, renderer)
380 clean_line, ismath = self._preprocess_math(line)
381 if clean_line:
--> 382 w, h, d = _get_text_metrics_with_cache(
383 renderer, clean_line, self._fontproperties,
384 ismath=ismath, dpi=self.get_figure(root=True).dpi)
385 else:
386 w = h = d = 0
File /data4/smurray/miniforge3/envs/21cmfast/lib/python3.11/site-packages/matplotlib/text.py:69, in _get_text_metrics_with_cache(renderer, text, fontprop, ismath, dpi)
66 """Call ``renderer.get_text_width_height_descent``, caching the results."""
67 # Cached based on a copy of fontprop so that later in-place mutations of
68 # the passed-in argument do not mess up the cache.
---> 69 return _get_text_metrics_with_cache_impl(
70 weakref.ref(renderer), text, fontprop.copy(), ismath, dpi)
File /data4/smurray/miniforge3/envs/21cmfast/lib/python3.11/site-packages/matplotlib/text.py:77, in _get_text_metrics_with_cache_impl(renderer_ref, text, fontprop, ismath, dpi)
73 @functools.lru_cache(4096)
74 def _get_text_metrics_with_cache_impl(
75 renderer_ref, text, fontprop, ismath, dpi):
76 # dpi is unused, but participates in cache invalidation (via the renderer).
---> 77 return renderer_ref().get_text_width_height_descent(text, fontprop, ismath)
File /data4/smurray/miniforge3/envs/21cmfast/lib/python3.11/site-packages/matplotlib/backends/backend_agg.py:215, in RendererAgg.get_text_width_height_descent(self, s, prop, ismath)
211 return super().get_text_width_height_descent(s, prop, ismath)
213 if ismath:
214 ox, oy, width, height, descent, font_image = \
--> 215 self.mathtext_parser.parse(s, self.dpi, prop)
216 return width, height, descent
218 font = self._prepare_font(prop)
File /data4/smurray/miniforge3/envs/21cmfast/lib/python3.11/site-packages/matplotlib/mathtext.py:86, in MathTextParser.parse(self, s, dpi, prop, antialiased)
81 from matplotlib.backends import backend_agg
82 load_glyph_flags = {
83 "vector": LoadFlags.NO_HINTING,
84 "raster": backend_agg.get_hinting_flag(),
85 }[self._output_type]
---> 86 return self._parse_cached(s, dpi, prop, antialiased, load_glyph_flags)
File /data4/smurray/miniforge3/envs/21cmfast/lib/python3.11/site-packages/matplotlib/mathtext.py:100, in MathTextParser._parse_cached(self, s, dpi, prop, antialiased, load_glyph_flags)
97 if self._parser is None: # Cache the parser globally.
98 self.__class__._parser = _mathtext.Parser()
--> 100 box = self._parser.parse(s, fontset, fontsize, dpi)
101 output = _mathtext.ship(box)
102 if self._output_type == "vector":
File /data4/smurray/miniforge3/envs/21cmfast/lib/python3.11/site-packages/matplotlib/_mathtext.py:2173, in Parser.parse(self, s, fonts_object, fontsize, dpi)
2170 result = self._expression.parseString(s)
2171 except ParseBaseException as err:
2172 # explain becomes a plain method on pyparsing 3 (err.explain(0)).
-> 2173 raise ValueError("\n" + ParseException.explain(err, 0)) from None
2174 self._state_stack = []
2175 self._in_subscript_or_superscript = False
ValueError:
f_{*,7}^{\\rm mol}
^
ParseSyntaxException: Expected {accent | symbol | function | operatorname | group | frac | dfrac | binom | genfrac | overset | underset | sqrt | overline | text | boldsymbol | substack}, found '\' (at char 9), (line:1, col:10)
---------------------------------------------------------------------------
ValueError Traceback (most recent call last)
File /data4/smurray/miniforge3/envs/21cmfast/lib/python3.11/site-packages/IPython/core/formatters.py:338, in BaseFormatter.__call__(self, obj)
336 pass
337 else:
--> 338 return printer(obj)
339 # Finally look for special method names
340 method = get_real_method(obj, self.print_method)
File /data4/smurray/miniforge3/envs/21cmfast/lib/python3.11/site-packages/IPython/core/pylabtools.py:152, in print_figure(fig, fmt, bbox_inches, base64, **kwargs)
149 from matplotlib.backend_bases import FigureCanvasBase
150 FigureCanvasBase(fig)
--> 152 fig.canvas.print_figure(bytes_io, **kw)
153 data = bytes_io.getvalue()
154 if fmt == 'svg':
File /data4/smurray/miniforge3/envs/21cmfast/lib/python3.11/site-packages/matplotlib/backend_bases.py:2155, in FigureCanvasBase.print_figure(self, filename, dpi, facecolor, edgecolor, orientation, format, bbox_inches, pad_inches, bbox_extra_artists, backend, **kwargs)
2152 # we do this instead of `self.figure.draw_without_rendering`
2153 # so that we can inject the orientation
2154 with getattr(renderer, "_draw_disabled", nullcontext)():
-> 2155 self.figure.draw(renderer)
2156 if bbox_inches:
2157 if bbox_inches == "tight":
File /data4/smurray/miniforge3/envs/21cmfast/lib/python3.11/site-packages/matplotlib/artist.py:94, in _finalize_rasterization.<locals>.draw_wrapper(artist, renderer, *args, **kwargs)
92 @wraps(draw)
93 def draw_wrapper(artist, renderer, *args, **kwargs):
---> 94 result = draw(artist, renderer, *args, **kwargs)
95 if renderer._rasterizing:
96 renderer.stop_rasterizing()
File /data4/smurray/miniforge3/envs/21cmfast/lib/python3.11/site-packages/matplotlib/artist.py:71, in allow_rasterization.<locals>.draw_wrapper(artist, renderer)
68 if artist.get_agg_filter() is not None:
69 renderer.start_filter()
---> 71 return draw(artist, renderer)
72 finally:
73 if artist.get_agg_filter() is not None:
File /data4/smurray/miniforge3/envs/21cmfast/lib/python3.11/site-packages/matplotlib/figure.py:3257, in Figure.draw(self, renderer)
3254 # ValueError can occur when resizing a window.
3256 self.patch.draw(renderer)
-> 3257 mimage._draw_list_compositing_images(
3258 renderer, self, artists, self.suppressComposite)
3260 renderer.close_group('figure')
3261 finally:
File /data4/smurray/miniforge3/envs/21cmfast/lib/python3.11/site-packages/matplotlib/image.py:134, in _draw_list_compositing_images(renderer, parent, artists, suppress_composite)
132 if not_composite or not has_images:
133 for a in artists:
--> 134 a.draw(renderer)
135 else:
136 # Composite any adjacent images together
137 image_group = []
File /data4/smurray/miniforge3/envs/21cmfast/lib/python3.11/site-packages/matplotlib/artist.py:71, in allow_rasterization.<locals>.draw_wrapper(artist, renderer)
68 if artist.get_agg_filter() is not None:
69 renderer.start_filter()
---> 71 return draw(artist, renderer)
72 finally:
73 if artist.get_agg_filter() is not None:
File /data4/smurray/miniforge3/envs/21cmfast/lib/python3.11/site-packages/matplotlib/axes/_base.py:3210, in _AxesBase.draw(self, renderer)
3207 if artists_rasterized:
3208 _draw_rasterized(self.get_figure(root=True), artists_rasterized, renderer)
-> 3210 mimage._draw_list_compositing_images(
3211 renderer, self, artists, self.get_figure(root=True).suppressComposite)
3213 renderer.close_group('axes')
3214 self.stale = False
File /data4/smurray/miniforge3/envs/21cmfast/lib/python3.11/site-packages/matplotlib/image.py:134, in _draw_list_compositing_images(renderer, parent, artists, suppress_composite)
132 if not_composite or not has_images:
133 for a in artists:
--> 134 a.draw(renderer)
135 else:
136 # Composite any adjacent images together
137 image_group = []
File /data4/smurray/miniforge3/envs/21cmfast/lib/python3.11/site-packages/matplotlib/artist.py:71, in allow_rasterization.<locals>.draw_wrapper(artist, renderer)
68 if artist.get_agg_filter() is not None:
69 renderer.start_filter()
---> 71 return draw(artist, renderer)
72 finally:
73 if artist.get_agg_filter() is not None:
File /data4/smurray/miniforge3/envs/21cmfast/lib/python3.11/site-packages/matplotlib/text.py:752, in Text.draw(self, renderer)
749 renderer.open_group('text', self.get_gid())
751 with self._cm_set(text=self._get_wrapped_text()):
--> 752 bbox, info, descent = self._get_layout(renderer)
753 trans = self.get_transform()
755 # don't use self.get_position here, which refers to text
756 # position in Text:
File /data4/smurray/miniforge3/envs/21cmfast/lib/python3.11/site-packages/matplotlib/text.py:382, in Text._get_layout(self, renderer)
380 clean_line, ismath = self._preprocess_math(line)
381 if clean_line:
--> 382 w, h, d = _get_text_metrics_with_cache(
383 renderer, clean_line, self._fontproperties,
384 ismath=ismath, dpi=self.get_figure(root=True).dpi)
385 else:
386 w = h = d = 0
File /data4/smurray/miniforge3/envs/21cmfast/lib/python3.11/site-packages/matplotlib/text.py:69, in _get_text_metrics_with_cache(renderer, text, fontprop, ismath, dpi)
66 """Call ``renderer.get_text_width_height_descent``, caching the results."""
67 # Cached based on a copy of fontprop so that later in-place mutations of
68 # the passed-in argument do not mess up the cache.
---> 69 return _get_text_metrics_with_cache_impl(
70 weakref.ref(renderer), text, fontprop.copy(), ismath, dpi)
File /data4/smurray/miniforge3/envs/21cmfast/lib/python3.11/site-packages/matplotlib/text.py:77, in _get_text_metrics_with_cache_impl(renderer_ref, text, fontprop, ismath, dpi)
73 @functools.lru_cache(4096)
74 def _get_text_metrics_with_cache_impl(
75 renderer_ref, text, fontprop, ismath, dpi):
76 # dpi is unused, but participates in cache invalidation (via the renderer).
---> 77 return renderer_ref().get_text_width_height_descent(text, fontprop, ismath)
File /data4/smurray/miniforge3/envs/21cmfast/lib/python3.11/site-packages/matplotlib/backends/backend_agg.py:215, in RendererAgg.get_text_width_height_descent(self, s, prop, ismath)
211 return super().get_text_width_height_descent(s, prop, ismath)
213 if ismath:
214 ox, oy, width, height, descent, font_image = \
--> 215 self.mathtext_parser.parse(s, self.dpi, prop)
216 return width, height, descent
218 font = self._prepare_font(prop)
File /data4/smurray/miniforge3/envs/21cmfast/lib/python3.11/site-packages/matplotlib/mathtext.py:86, in MathTextParser.parse(self, s, dpi, prop, antialiased)
81 from matplotlib.backends import backend_agg
82 load_glyph_flags = {
83 "vector": LoadFlags.NO_HINTING,
84 "raster": backend_agg.get_hinting_flag(),
85 }[self._output_type]
---> 86 return self._parse_cached(s, dpi, prop, antialiased, load_glyph_flags)
File /data4/smurray/miniforge3/envs/21cmfast/lib/python3.11/site-packages/matplotlib/mathtext.py:100, in MathTextParser._parse_cached(self, s, dpi, prop, antialiased, load_glyph_flags)
97 if self._parser is None: # Cache the parser globally.
98 self.__class__._parser = _mathtext.Parser()
--> 100 box = self._parser.parse(s, fontset, fontsize, dpi)
101 output = _mathtext.ship(box)
102 if self._output_type == "vector":
File /data4/smurray/miniforge3/envs/21cmfast/lib/python3.11/site-packages/matplotlib/_mathtext.py:2173, in Parser.parse(self, s, fonts_object, fontsize, dpi)
2170 result = self._expression.parseString(s)
2171 except ParseBaseException as err:
2172 # explain becomes a plain method on pyparsing 3 (err.explain(0)).
-> 2173 raise ValueError("\n" + ParseException.explain(err, 0)) from None
2174 self._state_stack = []
2175 self._in_subscript_or_superscript = False
ValueError:
f_{*,7}^{\\rm mol}
^
ParseSyntaxException: Expected {accent | symbol | function | operatorname | group | frac | dfrac | binom | genfrac | overset | underset | sqrt | overline | text | boldsymbol | substack}, found '\' (at char 9), (line:1, col:10)
<Figure size 2400x1400 with 38 Axes>
21-cm power spectra¶
[87]:
# define functions to calculate PS, following py21cmmc
from powerbox.tools import get_power
def compute_power(
box,
length,
n_psbins,
log_bins=True,
ignore_kperp_zero=True,
ignore_kpar_zero=False,
ignore_k_zero=False,
):
# Determine the weighting function required from ignoring k's.
k_weights = np.ones(box.shape, dtype=int)
n0 = k_weights.shape[0]
n1 = k_weights.shape[-1]
if ignore_kperp_zero:
k_weights[n0 // 2, n0 // 2, :] = 0
if ignore_kpar_zero:
k_weights[:, :, n1 // 2] = 0
if ignore_k_zero:
k_weights[n0 // 2, n0 // 2, n1 // 2] = 0
res = get_power(
box,
boxlength=length,
bins=n_psbins,
bin_ave=False,
get_variance=False,
log_bins=log_bins,
k_weights=k_weights,
)
res = list(res)
k = res[1]
if log_bins:
k = np.exp((np.log(k[1:]) + np.log(k[:-1])) / 2)
else:
k = (k[1:] + k[:-1]) / 2
res[1] = k
return res
def powerspectra(
brightness_temp, n_psbins=50, nchunks=10, min_k=0.1, max_k=1.0, logk=True
):
data = []
chunk_indices = list(
range(
0,
brightness_temp.n_slices,
round(brightness_temp.n_slices / nchunks),
)
)
if len(chunk_indices) > nchunks:
chunk_indices = chunk_indices[:-1]
chunk_indices.append(brightness_temp.n_slices)
for i in range(nchunks):
start = chunk_indices[i]
end = chunk_indices[i + 1]
chunklen = (end - start) * brightness_temp.cell_size
power, k = compute_power(
brightness_temp.brightness_temp[:, :, start:end],
(BOX_LEN, BOX_LEN, chunklen),
n_psbins,
log_bins=logk,
)
data.append({"k": k, "delta": power * k**3 / (2 * np.pi**2)})
return data
[ ]:
# do 5 chunks but only plot 1 - 4, the 0th has no power for minihalo models where xH=0
nchunks = 4
fig, axss = plt.subplots(
nchunks,
len(labelss),
sharex=True,
sharey=True,
figsize=(4 * len(labelss), 3 * (nchunks)),
subplot_kw={"xscale": "log", "yscale": "log"},
)
for pp, text in enumerate(texts):
labels = labels[pp]
factors = factors[pp]
axs = axss[:, pp]
for kk, label in enumerate(labels):
factor = factors[kk]
if var_params[pp] == "F_H2_SHIELD":
if factor > 1:
continue
inputs_var = inputs.evolve_input_structs(
USE_MINI_HALOS=(label!="NOmini"),
F_H2_SHIELD=1. - (1. - inputs.astro_params.F_H2_SHIELD)*factor,
)
else:
inputs_var = inputs.evolve_input_structs(
USE_MINI_HALOS=(label!="NOmini"),
**{var_params[pp] : inputs.astro_params.get(var_params[pp]) * factor}
)
if label == "reference":
lightcone = lightcone_fid
else:
lightcone = p21c.run_lightcone(
initial_conditions=initial_conditions,
inputs=inputs_var,
direc=cache,
)
PS = powerspectra(lightcone)
for ii in range(nchunks):
axs[ii].plot(
PS[ii + 1]["k"],
PS[ii + 1]["delta"],
color=colors[kk],
linestyle=linestyles[kk],
label=labels[kk],
lw=lws[kk],
)
if pp == len(texts) - 1 and kk == 0:
axs[ii].text(
0.99,
0.01,
"Chunk-%02d" % (ii + 1),
horizontalalignment="right",
verticalalignment="bottom",
transform=axs[ii].transAxes,
fontsize=15,
)
axs[0].text(
0.01,
0.99,
texts,
horizontalalignment="left",
verticalalignment="top",
transform=axs[0].transAxes,
fontsize=15,
)
axs[-1].set_xlabel("$k$ [Mpc$^{-3}$]", fontsize=15)
axs[-1].xaxis.set_tick_params(labelsize=15)
if pp == 0:
for ii in range(nchunks):
axs[ii].set_ylim(2e-1, 2e2)
axs[ii].set_ylabel("$k^3 P$", fontsize=15)
axs[ii].yaxis.set_tick_params(labelsize=15)
else:
for ii in range(nchunks - 1):
axs[ii].set_ylim(2e-1, 2e2)
axs[ii].set_ylabel("$k^3 P$", fontsize=0)
axs[ii].yaxis.set_tick_params(labelsize=0)
axss[0, 0].legend(loc="lower left", ncol=2, fontsize=13, fancybox=True, frameon=True)
plt.tight_layout()
fig.subplots_adjust(hspace=0.0, wspace=0.0)
Now you know how minihalo can shape the 21-cm signal!
[ ]: