{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "eefdf8a5-2b7f-45ec-93dd-44687396bb4f",
   "metadata": {},
   "source": [
    "# Ants model analyses\n",
    "\n",
    "_Copyright (C) 2023-, Joseph T. Lizier._\n",
    "_Distributed under GNU General Public License v3_\n",
    "\n",
    "This notebook template sets up reading in the data files for the Ants model simulations.\n",
    "\n",
    "First code block is just to start up the JVM. You may need to change your paths below (use the AutoAnalyser to show them for you):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "0eb38cde-2a28-4740-acf7-a76df9babf4d",
   "metadata": {},
   "outputs": [],
   "source": [
    "from jpype import *\n",
    "import numpy\n",
    "import sys\n",
    "# Our python data file readers are a bit of a hack, python users will do better on this:\n",
    "sys.path.append(\"../../../demos/python/\")\n",
    "import readFloatsFile\n",
    "\n",
    "if (not isJVMStarted()):\n",
    "    # Add JIDT jar library to the path\n",
    "    jarLocation = \"../../../infodynamics.jar\"\n",
    "    # Start the JVM (add the \"-Xmx\" option with say 1024M if you get crashes due to not enough memory space)\n",
    "    startJVM(getDefaultJVMPath(), \"-ea\", \"-Djava.class.path=\" + jarLocation, convertStrings=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2efe2949-d8ec-48a7-b6c7-a0569944b844",
   "metadata": {},
   "source": [
    "Next load in the sample data file and check that we can plot the ants positions ok.\n",
    "\n",
    "Make sure that you have either downloaded the sample data ([positionsx.txt](https://github.com/jlizier/jidt-cc-extras/raw/main/course/Module08-InfoTheoryAndSelfOrganisation/Ants/positionsx.txt), [positionsy.txt](https://github.com/jlizier/jidt-cc-extras/raw/main/course/Module08-InfoTheoryAndSelfOrganisation/Ants/positionsy.txt) and [headings.txt](https://raw.githubusercontent.com/jlizier/jidt-cc-extras/main/course/Module08-InfoTheoryAndSelfOrganisation/Ants/headings.txt)) to the current folder or generated sample data.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "cca3e625-d2e7-4fbd-bb65-d10bcf63c794",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "0b79a04099434d55a78405256cc0ec4e",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "IntProgress(value=0, max=300)"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/html": [
       "<video src=\"Ants-sample-video.mp4\" controls  >\n",
       "      Your browser does not support the <code>video</code> element.\n",
       "    </video>"
      ],
      "text/plain": [
       "<IPython.core.display.Video object>"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "[encoded data, 26428 characters]",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# 0. Load/prepare the data:\n",
    "positionsx = readFloatsFile.readFloatsFile(\"./positionsx.txt\")\n",
    "positionsy = readFloatsFile.readFloatsFile(\"./positionsy.txt\")\n",
    "headings = readFloatsFile.readFloatsFile(\"./headings.txt\")\n",
    "# As numpy arrays:\n",
    "positionsx = numpy.array(positionsx)\n",
    "positionsy = numpy.array(positionsy)\n",
    "headings = numpy.array(headings)\n",
    "\n",
    "# Plot a sample of the ants moving:\n",
    "import matplotlib.pyplot as plt\n",
    "import matplotlib.animation as animation\n",
    "from ipywidgets import IntProgress\n",
    "from IPython.display import display\n",
    "fig, ax = plt.subplots()\n",
    "plots = []\n",
    "timeStepsForPlot = min(positionsx.shape[0], 300)\n",
    "f = IntProgress(min=0, max=timeStepsForPlot) # instantiate the progress bar\n",
    "display(f) # display progress bar\n",
    "for t in range(timeStepsForPlot):\n",
    "    im = ax.scatter(positionsx[t,:], positionsy[t,:], marker='x', color='r', animated=True)\n",
    "    if (t == 0):\n",
    "         ax.scatter(positionsx[t,:], positionsy[t,:], marker='x', color='r')  # show an initial one first\n",
    "    # If you want to keep the plots to make a video:\n",
    "    plots.append([im])\n",
    "    # To save a png of each plot as you go, uncomment the following:\n",
    "    # plt.savefig('image-%03d.png' % t);\n",
    "    f.value += 1 # update progress bar\n",
    "plt.title('Ants positions');\n",
    "plt.xlabel('x');\n",
    "plt.ylabel('y');\n",
    "# And put the animations into a short video:\n",
    "ani = animation.ArtistAnimation(fig, plots, interval=100, blit=True,\n",
    "                                repeat_delay=1000)\n",
    "videoFilename = 'Ants-sample-video.mp4';\n",
    "ani.save(videoFilename)\n",
    "from IPython.display import Video\n",
    "Video(videoFilename)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "26bd6e3c-3dfa-4e42-a384-158783eec2e2",
   "metadata": {},
   "source": [
    "# Analysis 1 -- diversity of heading versus time\n",
    "\n",
    "This sample analysis looks at how the diversity (measured by entropy) in the positions of each ant changes across time.\n",
    "\n",
    "You can see that the diversity drops when the ants are making trails, with the largest 3 drops corresponding to when the trails to each food source are being sustained."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "39fb6614-ed61-4ada-89ce-fcd5a4a7f23e",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 47212 characters]",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# 0. Load/prepare the data:\n",
    "dataRaw = readFloatsFile.readFloatsFile(\"./positionsy.txt\")\n",
    "# As numpy array:\n",
    "data = numpy.array(dataRaw)\n",
    "# Transpose so that we have a column for each time step.\n",
    "# This is different to our usual approach of each ant being a variable in a column.\n",
    "data = data.T;\n",
    "\n",
    "mUtils = JPackage('infodynamics.utils').MatrixUtils\n",
    "\n",
    "# 1. Construct the calculator:\n",
    "calcClass = JPackage(\"infodynamics.measures.discrete\").EntropyCalculatorDiscrete\n",
    "calc = calcClass(8) # Using 8 discrete bins\n",
    "# 2. No other properties to set for discrete calculators.\n",
    "results = numpy.zeros((data.shape[1])) # There will be one entropy for each time step\n",
    "for t in range(data.shape[1]):\n",
    "    # For each time step:\n",
    "    variable = mUtils.discretise(JArray(JDouble, 1)(data[:,t].tolist()), 8)\n",
    "    \n",
    "    # 3. Initialise the calculator for (re-)use:\n",
    "    calc.initialise()\n",
    "    # 4. Supply the sample data:\n",
    "    calc.addObservations(variable)\n",
    "    # 5. Compute the estimate:\n",
    "    results[t] = calc.computeAverageLocalOfObservations()\n",
    "    \n",
    "    # print(\"H_Binned(col_0) = %.4f bits\" %\\\n",
    "    #    (results[t]))\n",
    "\n",
    "import matplotlib.pyplot as plt\n",
    "plt.plot(results)\n",
    "plt.xlabel('time')\n",
    "plt.ylabel('Entropy of y-position (bits)');"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "30bb2983-6308-494d-b115-7a783fc3b1e3",
   "metadata": {},
   "source": [
    "# Analysis 2 -- relationship of position and heading\n",
    "\n",
    "This sample analysis looks at how the information in x position about the heading of each ant changes across time.\n",
    "\n",
    "You can see that the MI increases when the ants are making trails, with the largest 3 increases corresponding to when the trails to each food source are being sustained.\n",
    "\n",
    "Subseuqent analysis included in the Matlab code investigates this in more subtle ways, but does not give us additional insights."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "de24081b-04e3-4460-bf6e-6ed97a3404eb",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 53216 characters]",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# 1. Construct the calculator:\n",
    "calcClass = JPackage(\"infodynamics.measures.continuous.kraskov\").MutualInfoCalculatorMultiVariateKraskov2\n",
    "calc = calcClass()\n",
    "# 2. Set any properties to non-default values:\n",
    "# No properties were set to non-default values\n",
    "\n",
    "misVersusTime = numpy.zeros((positionsx.shape[0])); # There will be one MI for each time step\n",
    "\n",
    "# Compute for samples taken across all pairs at each time step:\n",
    "for t in range(positionsx.shape[0]):\n",
    "    source = JArray(JDouble, 1)(positionsx[t, :].tolist())\n",
    "    destination = JArray(JDouble, 1)(headings[t, :].tolist())\n",
    "\n",
    "    # 3. Initialise the calculator for (re-)use:\n",
    "    calc.initialise()\n",
    "    # 4. Supply the sample data:\n",
    "    calc.setObservations(source, destination)\n",
    "    # 5. Compute the estimate:\n",
    "    misVersusTime[t] = calc.computeAverageLocalOfObservations()\n",
    "\n",
    "    # print(\"MI_Kraskov (KSG) alg. 2(time=%d) = %.4f nats\" %\\\n",
    "    #    (t, misVersusTime[t]))\n",
    "\n",
    "plt.plot(misVersusTime);\n",
    "plt.xlabel('time');\n",
    "plt.ylabel('MI from x position to heading (nats)');\n",
    "# Compare the times of the peaks in MI to when the trails are formed and collapse.\n",
    "# What is this detecting?"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9a84fc87-db74-4d67-81ae-5e47bcb0b093",
   "metadata": {},
   "source": [
    "# Analysis 3 -- mutual information averaged over all ant pairs\n",
    "\n",
    "This sample analysis looks at how the information in x position between all pairs of ants. We use pointwise values to rearrange this to look at how MI between the ants varies with time.\n",
    "\n",
    "You can see that the MI increases when the ants are making trails, with the largest 3 increases corresponding to when the trails to each food source are being sustained."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "7abb5b6d-5fc9-48e5-8751-d81c8a75a028",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 57772 characters]",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# using positionsx as loaded above\n",
    "\n",
    "avLocals = numpy.zeros((positionsx.shape[0])); # There will be one average MI across pairs for each time step\n",
    "\n",
    "# 1. Construct the calculator:\n",
    "calcClass = JPackage(\"infodynamics.measures.discrete\").MutualInformationCalculatorDiscrete\n",
    "calc = calcClass(8, 8, 0)\n",
    "# 2. No other properties to set for discrete calculators.\n",
    "mUtils = JPackage('infodynamics.utils').MatrixUtils\n",
    "\n",
    "# Compute for all pairs:\n",
    "for s in range(positionsx.shape[1]):\n",
    "    for d in range(positionsx.shape[1]):\n",
    "        # For each source-dest pair:\n",
    "        if (s == d):\n",
    "            continue\n",
    "        source = mUtils.discretise(JArray(JDouble, 1)(positionsx[:,s].tolist()), 8)\n",
    "        destination = mUtils.discretise(JArray(JDouble, 1)(positionsx[:,d].tolist()), 8)\n",
    "\n",
    "        # 3. Initialise the calculator for (re-)use:\n",
    "        calc.initialise()\n",
    "        # 4. Supply the sample data:\n",
    "        calc.addObservations(source, destination)\n",
    "        # 5. Compute the estimate:\n",
    "        result = calc.computeAverageLocalOfObservations()\n",
    "\n",
    "        # print(\"MI_Binned(col_%d -> col_%d) = %.4f nats\" %\\\n",
    "        #    (s, d, result))\n",
    "        locals = calc.computeLocalFromPreviousObservations(source, destination);\n",
    "        avLocals = avLocals + locals;\n",
    "\n",
    "avLocals = avLocals / (positionsx.shape[1]*(positionsx.shape[1]-1)); # Dividing by the number of pairs\n",
    "plt.plot(avLocals);\n",
    "plt.xlabel('time');\n",
    "plt.ylabel('Average local MI between pairs (bits)');"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "749dbd50-4fc1-447b-a12e-9a5b9bd4e526",
   "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
}
