{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "fe8e1d06-f91e-430d-855b-b4f618232057",
   "metadata": {},
   "source": [
    "# Module 10 -- Active Information Storage in Cellular Automata\n",
    "\n",
    "In this activity we will examine Active Information Storage (AIS) in Elementary Cellular Automata rule 54. CAs are wonderful canonical systems to use to explore dynamics of self-organisation in complex systems, in particular regarding emergent structures such as blinkers, gliders and collisions. If you have not played around with CAs before, you can explore dynamics of ECAs using NetLogo via `Sample models | Computer Science | Cellular Automata | CA 1D Elementary` in the Models Library (or access via [NetLogo web](http://www.netlogoweb.org/launch#http://www.netlogoweb.org/assets/modelslib/Sample%20Models/Computer%20Science/Cellular%20Automata/CA%201D%20Elementary.nlogo) if you don't have an installation): try rules 2, 8, 18, 30, 54, 90, 110, 146, and 254 to see a variety of dynamics.\n",
    "\n",
    "In this example we will compute the appropriate embedding length manually, and then go on to examine local AIS values and how they relate to emergent structure.\n",
    "\n",
    "Make sure that you have **downloaded** the sample time-series as directed from the tutorial page. This data set was generated by running CA rule 54 from a random initial condition. The data set contains 10000 cells and 600 time steps (with periodic boundary conditions) -- we use more cells than time steps such that we are sampling the transient dynamics before a steady state is reached.\n",
    "\n",
    "1. Open the AutoAnalyser and select Active Information Storage. Select a Discrete estimator, select your CA data file, tick the checkbox for `All variables?` and click `Generate Code and Compute`.\n",
    "2. Copy and paste the generated code into the code cells below."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9c419f90-69df-439b-b2d2-986335b07d65",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Paste the import and JVM startup lines here:\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "31d1ff27-5047-4da2-816b-5c770176789b",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Paste the code loading and analysing the data here:\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7e928d5b-6bcb-4222-962a-395c79aec001",
   "metadata": {},
   "outputs": [],
   "source": [
    "# This code block plots a sample of 100 time steps for 100 cells of the CA:\n",
    "import matplotlib.pyplot as plt\n",
    "plt.imshow(data[:99,:99]) # Plotting the first 100 time steps for first 100 cells\n",
    "plt.xlabel('cells');\n",
    "plt.ylabel('timestep');\n",
    "cbar = plt.colorbar()\n",
    "cbar.set_label('CA values');"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "85cb1f19-8315-4e38-b38c-183ac26a0c27",
   "metadata": {},
   "source": [
    "3. Examine how the code is computing the AIS for each cell in the CA. Now run the code cells above. Notice the calculation for each CA cell printed above, and notice the variance in these. We would get a much more accurate result if we could make an AIS estimation by pooling all of the samples for the different CA cells together. Reflect on why pooling the samples for technically different time-series processes is valid here, and think of a data set where it would not be a valid approach.\n",
    "4. Change the code -- perhaps copy a new version into the cell below and change that -- so that, instead of making a separate calculation for each cell, we make one calculation where the samples from all cells are used together. _Hint_: think of which calculation steps should remain inside the loop over variables, and which should now go outside (before or after) -- you can go back to the flow chart for how to use the JIDT estimators in the *Introduction to JIDT lecture*. When you print the final result of the calculation using samples from all variables, also check the number of samples that were used via `calc.getNumObservations()`, which should report a number of observations consistent with having used all of the samples from all cells. If it still reports 599, obviously you've only used samples from one cell. Run the code -- your final answer (if using the sample data set) should be 0.1426 bits."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3d3925e8-d81f-4fb2-b2a6-c16b4825c4f8",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Paste the code to start to pool the samples together here:\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ed0610bb-cb0e-4b1a-b94b-888c5b3f7e3f",
   "metadata": {},
   "source": [
    "5. Now, we want to work out the most appropriate value of the history embedding `k` to use via the maximum bias-corrected AIS. JIDT doesn't yet offer an automated way of doing this for AIS Discrete, so we'll do it manually (in any case, this is instructive for us in seeing how an automated approach works).\n",
    "    First take a copy of your code pooling the samples above, and paste it into the code cell below. Then:\n",
    "    1. Add a loop (after loading in the data) over the rest of the calculation, looping over values of `k` from 1 to 20.\n",
    "    2. Supply your new value of k to the JIDT estimator. To see how if you're not sure, go back to the AutoAnalyser, change the `k` parameter, generate the code and see where the new value is inserted.\n",
    "    3. Change the print statement for your results to include `k` as well.\n",
    "    4. Now, the AIS discrete estimator does not include bias correction, so we'll do this manually. Go back to the AutoAnalyser, click the checkbox for `Add stat signif.?` and `analytically?` and generate new code. Go to the Python tab and copy the code for step 6 into the cell below. It should go after where you compute the result. Then, after the analytic surrogate distribution is returned in the object `measDist`, add a line to retrieve the bias via: `bias = measDist.getMeanOfDistribution();`\n",
    "    5. Now, store all of the computed results and bias values in an array for each value of `k`, and generate the bias-corrected AIS values (`= results - bias`).\n",
    "    6. Plot the raw computed result, the bias and bias corrected values versus `k`. Identify the optimal value of `k` to use from the maximum bias-corrected AIS. Notice how the bias rises quickly after this, indicating the onset of _undersampling_ for large `k`, with respect to the given number of samples. My result is **k=15**, but statistical fluctuations in your data set could change that slightly."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e2865a32-79dc-46d8-b643-164639927464",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Paste the code from above here:\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "328cb9c7-00de-48d3-ac41-0a5a4972eacd",
   "metadata": {},
   "source": [
    "6. Now we move on to compute the **local AIS** values using the value of `k` we have determined. Copy the code from the previous cell into a new code cell below, so that we can change it but you can come back later to what you did previously.\n",
    "    1. Remove the loop over `k`, fixing it simply to the value that you determined. You'll also need to remove the arrays where you stored values for each `k`, and where you plot the values. The code should still compute one value of AIS over all of the cells.\n",
    "    2. After the value is computed, we'll insert new code to compute the local values.  At the end of your code block after the calculation is complete, add another loop over all cells (copying the loop structure over variables that is already earlier in the code). Inside the loop, select the column data for the given cell (as is done in the earlier loop), and then compute the local AIS values along that variable's column using: `localAISForThisCell = calc.computeLocalFromPreviousObservations(variable);`.\n",
    "    3. Store all of these columns of local AIS values together in one 2D matrix for all of the cells called locals, i.e. `locals[:,v] = localAISForThisCell;`. (You should initialise locals before the loop)\n",
    "    4. After the loop, let's plot a small segment of these local values. A simple way to do this would be:\n",
    "       ```python\n",
    "       plt.imshow(locals[:99,:99]) # Plotting the first 100 time steps for first 100 cells\n",
    "       plt.xlabel('cells');\n",
    "       plt.ylabel('timestep');\n",
    "       cbar = plt.colorbar()\n",
    "       cbar.set_label('AIS (bits)');\n",
    "       ```\n",
    "       However the colours could be arranged a little better (so that we use blue for positive AIS and red for negative). So instead first import some extra libraries from matplotlib: `import matplotlib.cm as cm; import matplotlib.colors as mplcolors` and then add the following argument to the `plt.imshow()` call: `cmap=cm.coolwarm_r, norm=mplcolors.CenteredNorm()`\n",
    "    5. Run the code!\n",
    "    6. (You can ignore the first k time steps here being 0 -- it's simply not possible to compute a local value until we've accumulated k time steps.)\n",
    "    7. Observe which emergent structures execute the strongest information storage (you may find it useful to replot the raw values here again for comparison). Does this align with our intuition from the lecture? What do the areas with negative local AIS values mean? (Recall the meaning of negative values of local mutual information, and think about what the past embedded state is telling us when the local AIS goes negative)\n",
    "    8. For further details, you can compare to Fig. 4b and section IV.E (and references therein) of the JIDT paper.\n",
    "\n",
    "See sample results on the tutorial page / solution notebook. Note that your run, from a different initial condition, will likely be different in the details, but the comparison of local AIS values between blinkers/domains and gliders should qualitatively be the same."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e89b5ca3-86d0-4014-8645-f0bfea4ae102",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Paste the completed code above (which finds the optimum value of k) here:\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e367dd0c-8d8e-4ece-8147-591c317d9c81",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.10.12"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}
