<jupyter_start><jupyter_text>Analyzing Model Grids / Monte Carlo Simulations Once you have a model grid in hand, there are a slew of built-in analysis routines that you might consider using. For the rest of this example,we'll assume you have completed the [Example 2-D Grid](example_grid), and have the associated set of fileswith prefix ``test_2d_grid``. To begin, initialize an instance of the analysis class, ``ModelSet``:<jupyter_code>import ares

anl = ares.analysis.ModelSet('test_2d_grid')<jupyter_output>##############################################################################################################
####                                         Analysis: Model Set                                          ####
##############################################################################################################
#### ---------------------------------------------------------------------------------------------------- ####
#### Basic Information                                                                                    ####
#### ---------------------------------------------------------------------------------------------------- ####
#### path        : ./                                                                                     ####
#### prefix      : test_2d_grid                                                                           ####
#### N-d         : 2                                                                                      ####
[...]<jupyter_text>First, let's verify that we've surveyed the part of parameter space we intended to:<jupyter_code>ax = anl.Scatter(anl.parameters)<jupyter_output><empty_output><jupyter_text>Basic Inspection Now, the kind of analysis you can do will be limited by what quantitieswere saved for each model. Recall (from [Example 2-D Grid](example_grid)) that we have the following quantities at our disposal:* 21-cm brightness temperature, ``dTb``.* Kinetic temperature of the IGM, ``igm_Tk``.* HII region volume filling factor, ``cgm_h_2``.* CMB optical depth, ``tau_e``.* Position of 21-cm emission maximum, ``z_D`` and ``dTb_D``.Let's have a look at how the ionization history depends on our two parameters ofinterest in this example, ``tanh_xz0`` and ``tanh_xdz``,<jupyter_code>ax = anl.Scatter(['tanh_xz0', 'tanh_xdz'], c='cgm_h_2', ivar=[None,None,10], fig=2, edgecolors='none')<jupyter_output># Loaded test_2d_grid.000.blob_1d.cgm_h_2.pkl<jupyter_text>**NOTE:** All ``ModelSet`` functions accept ``fig`` as an optional keyword     argument, which you can set to any integer to open plots in a new window.    The keyword argument ``ivar`` is short for "independent variables" -- it is ``None`` by default. However, because we have chosen to plot ``cgm_h_2``, which is a 1-D blob, we must specify the redshift of interest. Recall that we have access to integer redshifts in the interval $5 \leq z \leq 20$, or check for yourself:<jupyter_code>print(anl.blob_ivars)<jupyter_output>(None, [array([ 5,  6,  7,  8,  9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20])])<jupyter_text>So our choice of $z=10$ should be OK.  If you forget what fields are available for analysis, see:<jupyter_code>print(anl.blob_names, anl.blob_ivars)<jupyter_output>(['z_D', 'dTb_D', 'tau_e'], ['cgm_h_2', 'igm_Tk', 'dTb']) (None, [array([ 5,  6,  7,  8,  9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20])])<jupyter_text>**NOTE:** Calling quantities of interest `blobs` was inspired by the arbitrary meta-data blobs in [emcee](http://dan.iel.fm/emcee/current/).    For more 21-cm-focused analyses, you may want to view how the extrema in theglobal 21-cm signal change as a function of the model parameters:<jupyter_code># Scatterplot showing where emission peak occurs
ax = anl.Scatter(['z_D', 'dTb_D'], fig=4)<jupyter_output># Loaded test_2d_grid.000.blob_0d.z_D.pkl
# Loaded test_2d_grid.000.blob_0d.dTb_D.pkl<jupyter_text>or, color-code points by CMB optical depth,<jupyter_code>ax = anl.Scatter(['z_D', 'dTb_D'], c='tau_e', fig=5, edgecolors='none')<jupyter_output># Loaded test_2d_grid.000.blob_0d.tau_e.pkl<jupyter_text>You can also create your own derived quantities. A very simple example is to convert redshifts to observed frequencies,<jupyter_code># 1420.4057 is the rest frequency of the 21-cm line in MHz
anl.DeriveBlob(expr='1420.4057 / (1. + x)', varmap={'x': 'z_D'}, name='nu_D')<jupyter_output><empty_output><jupyter_text>This will create a new blob, called ``nu_D``, that can be used for subsequent analysis. For example,<jupyter_code># Scatterplot showing where emission peak occurs
ax = anl.Scatter(['nu_D', 'dTb_D'], c='tau_e', fig=6, edgecolors='none')<jupyter_output><empty_output><jupyter_text>Problem Realizations You may have noticed that in this model grid there are three realizations whose emission maxima seem to occur at $\delta T_b \approx 0$. In general, this is possible, but given the regularity of the grid points in parameter space it seems unlikely that any individual model would stray substantially from the locus of all other models.To inspect potentially problematic realizations, it is first useful to isolate them from the rest. You can select them visually by first invoking (**NOTE:** this won't work in a Jupyter notebook):<jupyter_code>anl.SelectModels()<jupyter_output><empty_output><jupyter_text>and then clicking and dragging within the plot window to define a rectangle, starting from its upper left corner (click) and ending with its bottom right corner (release). The set of models bounded by this rectangle will be saved as a new ``ModelSet`` object that can be used just like the original one. Each successive "slice" will be saved as attributes ``slice_0``, ``slice_1``, etc. that you can assign to a new variable, as, e.g.<jupyter_code>#slc0 = anl.slice_0
#slc0.Scatter(['nu_D', 'dTb_D'], c='tau_e', fig=7, edgecolors='none')<jupyter_output><empty_output><jupyter_text>Alternatively, you can specify a rectangle by hand. For example,<jupyter_code>slc = anl.Slice([100, 120, 0, 10], pars=['nu_D', 'dTb_D'])<jupyter_output>##############################################################################################################
####                                         Analysis: Model Set                                          ####
##############################################################################################################
#### ---------------------------------------------------------------------------------------------------- ####
#### Basic Information                                                                                    ####
#### ---------------------------------------------------------------------------------------------------- ####
#### path        : ./                                                                                     ####
#### prefix      : test_2d_grid                                                                           ####
#### N-d         : 2                                                                                      ####
[...]<jupyter_text>extracts all models with ``100 <= nu_D <= 120`` and ``0 <= dTb_D <= 10``. Check:<jupyter_code>slc.Scatter(['nu_D', 'dTb_D'], c='tau_e', fig=8, edgecolors='none')<jupyter_output># Loaded test_2d_grid.000.blob_0d.dTb_D.pkl
# Loaded test_2d_grid.000.blob_0d.tau_e.pkl<jupyter_text>If you wanted to examine models in more detail, you could re-run them. Collecting the parameter dictionaries required to do so is easy<jupyter_code>kwargs_list = slc.AssembleParametersList(include_bkw=True)<jupyter_output>WARNING: Any un-pickleable kwargs will not have been saved in test_2d_grid.binfo.pkl!<jupyter_text>This routine returns a list in which each element is a dictionary of parameters for a single model. The keyword argument ``include_bkw`` controls whether the "base kwargs," i.e., those that are shared by all models in the grid, are included in each list element. If they are (as above), then any individual dictionary can be used to initialize a simulation. For example:<jupyter_code>ax = None
for kwargs in kwargs_list[0:3]:
    sim = ares.simulations.Global21cm(**kwargs)
    sim.run()
    ax, zax = sim.GlobalSignature(color='b', alpha=0.5, ax=ax)<jupyter_output><empty_output><jupyter_text>If you've got models that seem to have something wrong with them, sending me the dictionary (or a list of them as above) will help a lot. Just do something like:<jupyter_code>import pickle
    f = open('problematic_models.pkl', 'wb')
    pickle.dump(kwargs_list, f)
    f.close()<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.19951083158630328, 'fraction_numerical': 0.019916142557651992, 'mean_word_length': 2.7530594405594404, 'pattern_counts': {'":': 0, '<': 60, '<?xml version=': 0, '>': 56, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 44, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '16', 'dir': 'jupyter-structured-clean-dedup', 'id': '522653', 'license': 'MIT', 'n_tokens_mistral': 2466, 'n_tokens_neox': 2307, 'n_words': 797, 'path': 'docs/examples/example_grid_analysis.ipynb', 'repo_name': 'mirochaj/ares', 'size': '8586'}