<jupyter_start><jupyter_text>Simple Parameter Study: 2-D Model Grid Often we want to study how the 21-cm signal changes over a range of parameters. We can do so using the ModelGrid class, and use numpy arrays to represent the range of values we’re interested in.Before we start, the few usual imports:<jupyter_code>%pylab inline
import ares
import numpy as np
import matplotlib.pyplot as pl<jupyter_output>Populating the interactive namespace from numpy and matplotlib<jupyter_text>Efficient Example: *tanh* model for the global 21-cm signal Before we run a set of models, we need to decide what quantities we’d like to save. For a detailed description of how to do this in general cases, check out [Inline Analysis](../example_inline_analyis).For now, let’s save the redshift and brightness temperature of the global 21-cm emission maximum, which we dub "Turning Point D", and the CMB optical depth,<jupyter_code>blobs_scalar = ['z_D', 'dTb_D', 'tau_e']<jupyter_output><empty_output><jupyter_text>in addition to the ionization, thermal, and global 21-cm histories at redshifts between 5 and 20 (at $\Delta z = 1$ increments),<jupyter_code>blobs_1d = ['cgm_h_2', 'igm_Tk', 'dTb']
blobs_1d_z = np.arange(5, 21)<jupyter_output><empty_output><jupyter_text>**NOTE:** For a complete listing of ideas for 1-D blobs see the [Fields Listing](../fields.html). Now, we’ll make a dictionary full of parameters that will get passed to every global 21-cm signal calculation. In addition to the blobs, we’ll set ``tanh_model=True`` to speed things up (see next section regarding physical models), and ``problem_type=101``:<jupyter_code>base_pars = \
{
 'problem_type': 101,
 'tanh_model': True,
 'blob_names': [blobs_scalar, blobs_1d],
 'blob_ivars': [None, [('z', blobs_1d_z)]],
 'blob_funcs': None,
}<jupyter_output><empty_output><jupyter_text>and create the ``ModelGrid`` instance,<jupyter_code>mg = ares.inference.ModelGrid(**base_pars)<jupyter_output><empty_output><jupyter_text>At this point we have yet to specify which parameters will define the axes of the model grid. Since we set ``tanh_model=True``, we have 9 parameters to choose from: a step height, step width, and step redshift for the Ly-$\alpha$, thermal, and ionization histories:* Ly-$\alpha$ history parameters: ``tanh_J0``, ``tanh_Jz0``, ``tanh_Jdz``.* Thermal history parameters: ``tanh_T0``, ``tanh_Tz0``, ``tanh_Tdz``.* Ionization history parameters: ``tanh_x0``, ``tanh_xz0``, ``tanh_xdz``.The Ly-$\alpha$ step height, ``tanh_J0``, must be provided in units of $J_{21} = 10^{-21} \mathrm{erg} \ \mathrm{s}^{-1} \ \mathrm{cm}^{-2} \ \mathrm{Hz}^{-1} \ \mathrm{sr}^{-1}$, while the temperature step height is assumed to be in Kelvin. The ionization step height should not exceed unity -- in fact it's safe to assume ``tanh_x0=1`` (we know that reionization *ends*!). See [Harker et al. (2016)]( ) for more information about the *tanh* model.Let’s take the reionization redshift, ``tanh_xz0``, and duration, ``tanh_xdz``, and sample them over a reasonable redshift interval with a spacing of $\Delta z = 0.2$,<jupyter_code>z0 = np.arange(6, 12.2, 0.2)
dz = np.arange(0.2, 8.2, 0.2)<jupyter_output><empty_output><jupyter_text>Now, we just set the ``axes`` attribute to a dictionary containing the array of values for each parameter:<jupyter_code>mg.axes = {'tanh_xz0': z0, 'tanh_xdz': dz}<jupyter_output><empty_output><jupyter_text>To run,<jupyter_code>mg.run('test_2d_grid', clobber=True, save_freq=100)<jupyter_output>grid:   0% |                                                                                  | ETA:   0:02:12<jupyter_text>To speed things up, you could increase the grid spacing. Or, execute the above in parallel as a Python script (assuming you have MPI and mpi4py installed).**NOTE:** If the model grid doesn’t finish running, that’s OK! Simply re-execute the above command with ``restart=True`` as an additional keyword argument and it will pick up where it left off.To analyze the results, create an analysis instance,<jupyter_code>anl = ares.analysis.ModelSet('test_2d_grid')<jupyter_output>##############################################################################################################
####                                         Analysis: Model Set                                          ####
##############################################################################################################
#### ---------------------------------------------------------------------------------------------------- ####
#### Basic Information                                                                                    ####
#### ---------------------------------------------------------------------------------------------------- ####
#### path        : ./                                                                                     ####
#### prefix      : test_2d_grid                                                                           ####
#### N-d         : 2                                                                                      ####
[...]<jupyter_text>and, for example, plot the 2-d parameter space with points color-coded by ``tau_e``,<jupyter_code>ax1 = anl.Scatter(anl.parameters, c='tau_e', fig=1)<jupyter_output># Loaded test_2d_grid.000.blob_0d.tau_e.pkl<jupyter_text>or instead, the position of the emission maximum with the same color coding:<jupyter_code>ax2 = anl.Scatter(['z_D', 'dTb_D'], c='tau_e', fig=2)<jupyter_output># Loaded test_2d_grid.000.blob_0d.z_D.pkl
# Loaded test_2d_grid.000.blob_0d.dTb_D.pkl<jupyter_text>**NOTE:** In general, you may want to save ``'z_C'`` and ``'dTb_C'``, i.e., the location of the global 21-cm absorption feature. But, in this case since we've only varied parameters of the ionization history, that point will not change and so saving it as a blob is unnecessary.     See :doc:`example_grid_analysis` for more information.  Accessing the Data Directly If you'd like to access the data directly for further manipulation, you'll be looking at the following attributes of the ``ModelSet`` class:* ``chain``, which is a 2-D array with dimensions (number of models, number dimensions).* ``get_blob``, which is a function that can be used to read-in blobs from disk.**NOTE:** The ``chain`` attribute is referred to as such because is analogous to an MCMC chain, but rather than random samples of the posterior distribution, it represents "samples" on a structured mesh.For example, to retrieve the samples of the ``test_2d_grid`` dataset above, you could do:<jupyter_code># Just the names of the axes
x, y = anl.parameters 
    
xdata, ydata = anl.chain[:,0], anl.chain[:,1]<jupyter_output><empty_output><jupyter_text>or equivalently,<jupyter_code>xdata, ydata = anl.chain.T<jupyter_output><empty_output><jupyter_text>And to plot the samples,<jupyter_code>pl.scatter(xdata, ydata)
pl.xlabel(x)
pl.ylabel(y)<jupyter_output><empty_output><jupyter_text>To extract blobs, you could do<jupyter_code>QHII = anl.get_blob('cgm_h_2')
print(QHII.shape)<jupyter_output># Loaded test_2d_grid.000.blob_1d.cgm_h_2.pkl
(1240, 16)<jupyter_text>Notice that the first dimension of ``QHII`` is the same as the first dimension of ``chain`` -- just the number of samples in the ModelGrid. The second dimension, however, is different. Now, rather than representing the dimensionality of the parameter space, it represents the dimensionality of this particular blob. Why 16 elements? Because our blobs were setup such that the quantities ``cgm_h_2``, ``igm_Tk``, and ``dTb`` were recorded at all redshifts in ``np.arange(5, 21)``, which has 16 elements.So, we could for example color-code the points in our previous plot by the volume-averaged ionization fraction at :math:`z=10` by doing:<jupyter_code>pl.scatter(xdata, ydata, c=QHII[:,5], edgecolors='none')<jupyter_output><empty_output><jupyter_text>If you forget the properties of a blob, you can type<jupyter_code>group, element, nd, shape = anl.blob_info('cgm_h_2')
print(group, element, nd, shape)<jupyter_output>1 0 1 (16,)<jupyter_text>which returns the index of blob group, index of the element within that group, dimensionality of the blob, and the shape of blob. This can be useful, for example, to automatically figure out the independent variables for a blob:<jupyter_code># Should be 10 (redshift of interest above)
anl.blob_ivars[group][element][5]<jupyter_output><empty_output><jupyter_text>All of the built-in analysis routines are structured so that you don't have to think about these things on a regular basis if you don't want to!      More Expensive Models Setting ``tanh_model=True`` sped things up considerably in the previous example. In general, you can run grids varying any *ARES* parameters you like, just know that physical models (i.e., those with ``tanh_model=False``) generally take a few seconds each, whereas the $tanh$ model takes much less than a second for one model.For example, to repeat the previous exercise for a physical model, you could replace this commands:<jupyter_code>z0 = np.arange(6, 12, 0.1)
dz = np.arange(0.1, 8.1, 0.1)
mg.axes = {'tanh_xz0': z0, 'tanh_xdz': dz}<jupyter_output><empty_output><jupyter_text>with (for example)<jupyter_code>fX = np.logspace(-1, 1, 21)
Tmin = np.logspace(3, 5, 21)
mg.axes = {'fX': z0, 'Tmin': dz}<jupyter_output><empty_output><jupyter_text>In some cases -- e.g., when ``Tmin`` or ``pop_Tmin`` is an axis of the model grid -- load-balancing can be very advantageous. Just execute the following command before running the grid:<jupyter_code>mg.LoadBalance(method=1, par='Tmin') # or 'pop_Tmin'<jupyter_output><empty_output><jupyter_text>The ``method=1`` setting assigns all models with common ``Tmin`` values to the same processor. This helps because *ARES* knows that it need only generate lookup tables for $df_{\mathrm{coll}} / dz$ (which determines the star formation rate density in the simplest models; see [this example](example_gs_standard)) once per value of ``Tmin``, which means you save a little bit of runtime at the outset of each calculation.    There is also a ``method=2`` option for load balancing, which is advantageous if the runtime of individual models is strongly correlated with a given parameter. In this case, the models will be sorted such that each processor gets a (roughly) equal share of the models for each value of the input ``par``. It helps to imagine the grid points of our 2-D parameter space color-coded by processor ID number: the resulting image for ``method=2`` is simply the transpose of the image you'd get for ``method=1``.If the edges of your parameter space correspond to rather extreme    models you might find that the calculations grind to a halt. This can be a big problem because you'll end up with one or more processors spinning their wheels while the rest of the processors continue. One way of dealing with this is to set an "alarm" that will be tripped if the runtime of a particular model exceeds some user-defined value. For example, before running a model grid, you might set:<jupyter_code>mg.timeout = 60<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.1384254666786906, 'fraction_numerical': 0.015871584453061593, 'mean_word_length': 3.9865107913669067, 'pattern_counts': {'":': 0, '<': 82, '<?xml version=': 0, '>': 82, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 22, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '22', 'dir': 'jupyter-structured-clean-dedup', 'id': '274671', 'license': 'MIT', 'n_tokens_mistral': 3365, 'n_tokens_neox': 3189, 'n_words': 1249, 'path': 'docs/examples/example_grid.ipynb', 'repo_name': 'mirochaj/ares', 'size': '11089'}