{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "26854b66-f08d-4126-bb71-eb97744e92b0",
   "metadata": {},
   "source": [
    "# Estimators and parameters for applying TE to empirical data -- Sleep Apnoea\n",
    "\n",
    "## A. Active information storage in breath rate\n",
    "\n",
    "We'll begin our modelling process by accounting for the past of the breath rate time series, using active information storage. Recall that the AIS represents the first step in an information regression of modelling the dynamic updates of the target.\n",
    "\n",
    "1. Open the AIS AutoAnalyser, selecting the `Kraskov (KSG)` estimator, and the `SFI-heartRate_breathVol_bloodOx-extract.txt` data set.\n",
    "2. Change the variable column to 1 (breath rate).\n",
    "3. Recall from our analysis of statistical significance in Module 7 that controlling for autocorrelation in samples here is quite important. So set the property `DYN_CORR_EXCL` to 15 (which was the maximum of the autocorrelation lengths we identified for the heart and breath rates).\n",
    "4. Let's ask JIDT to set embedding parameters for the breath rate automatically for us, allowing long dimension and lag since we know that these are complex time series:\n",
    "    1. Set `AUTO_EMBED_METHOD` to `MAX_CORR_AIS` to embed both source and target.\n",
    "    2. Set `AUTO_EMBED_METHOD_K_SEARCH_MAX` to `10` to search for an embedding length up to 10.\n",
    "    3. Set `AUTO_EMBED_METHOD_TAU_SEARCH_MAX` to `10` to search for an embedding delay up to 10 (allowing a longer lag than usual since the time series are complex)\n",
    "5. Press `Generate code and compute`:\n",
    "    * Inspect the resulting embedding parameters as displayed in the `k_HISTORY` and `TAU` values in the properties section.\n",
    "    * You should see a result around 0.59 nats, which is reasonably large (as is usually the case for modelling information storage for continuous valued processes sampled in time, since they are usually slowly varying).\n",
    "6. Copy the code into the following code cells and run them."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "c980c0ae-ecbb-4c97-9ec3-eefc17b9929c",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Paste import and JVM startup lines here:\n",
    "\n",
    "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(\"/home/joseph/JIDT/infodynamics-dist-1.6.1/demos/python\")\n",
    "import readFloatsFile\n",
    "\n",
    "if (not isJVMStarted()):\n",
    "    # Add JIDT jar library to the path\n",
    "    jarLocation = \"/home/joseph/JIDT/infodynamics-dist-1.6.1/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": "code",
   "execution_count": 2,
   "id": "6641436c-ca3c-467e-b5e1-e9d939469c0b",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "AIS_Kraskov (KSG)(col_1) = 0.5848 nats\n"
     ]
    }
   ],
   "source": [
    "# Paste analysis code here:\n",
    "\n",
    "# 0. Load/prepare the data:\n",
    "dataRaw = readFloatsFile.readFloatsFile(\"/home/joseph/JIDT/infodynamics-dist-1.6.1/demos/data/SFI-heartRate_breathVol_bloodOx-extract.txt\")\n",
    "# As numpy array:\n",
    "data = numpy.array(dataRaw)\n",
    "variable = JArray(JDouble, 1)(data[:,1].tolist())\n",
    "\n",
    "# 1. Construct the calculator:\n",
    "calcClass = JPackage(\"infodynamics.measures.continuous.kraskov\").ActiveInfoStorageCalculatorKraskov\n",
    "calc = calcClass()\n",
    "# 2. Set any properties to non-default values:\n",
    "calc.setProperty(\"AUTO_EMBED_METHOD\", \"MAX_CORR_AIS\")\n",
    "calc.setProperty(\"AUTO_EMBED_K_SEARCH_MAX\", \"10\")\n",
    "calc.setProperty(\"AUTO_EMBED_TAU_SEARCH_MAX\", \"10\")\n",
    "calc.setProperty(\"DYN_CORR_EXCL\", \"15\")\n",
    "# 3. Initialise the calculator for (re-)use:\n",
    "calc.initialise()\n",
    "# 4. Supply the sample data:\n",
    "calc.setObservations(variable)\n",
    "# 5. Compute the estimate:\n",
    "result = calc.computeAverageLocalOfObservations()\n",
    "\n",
    "print(\"AIS_Kraskov (KSG)(col_1) = %.4f nats\" %\\\n",
    "    (result))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d437cf2c-58fe-4b90-b32b-4df7f38efc54",
   "metadata": {},
   "source": [
    "  The following code cell plots the breath rate versus time (like the example on the tutorial page) -- run this code cell:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "f2866af4-450b-495a-9029-ac650d858c29",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 58292 characters]",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Plot the breath rate versus sample index:\n",
    "import matplotlib.pyplot as plt\n",
    "plt.figure()\n",
    "plt.plot(variable)\n",
    "plt.title('Breath rate versus sample index')\n",
    "plt.xlabel('n')\n",
    "plt.ylabel('Breath rate')\n",
    "plt.show();"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "808795d4-cb8e-45fd-b604-7c7bb33e96e8",
   "metadata": {},
   "source": [
    "The next code cell makes an additional plot -- inspect the code in the cell below:\n",
    "  * The first two lines of second code cell show an example of how to read out an automatically determined embedding parameter (or any other property for that matter) in code. This may be useful e.g. if you want to set the determined embedding for another calculation.\n",
    "  * The subsequent lines compute and then plot the local AIS values against the corresponding next value of the breath rate, so that we can see where in the dynamics the past is most strongly informative. (See sample plot on the tutorial webpage).\n",
    "\n",
    "7. Run this next code cell, and examine the plot -- what can you see about how the AIS relates to the dynamics? I.e. for which values of the breath rate do we see extreme values of the AIS?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "19a233d5-1ad1-4983-a37c-f0b54ede5199",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 89964 characters]",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Now make an additional plot (explained in the text below)\n",
    "k_str = calc.getProperty('k_HISTORY'); # Returns a string for the property name: useful where setting the property again\n",
    "k = int(k_str); # Converts the string to an integer, useful if we need to use it in later code\n",
    "localAIS = calc.computeLocalOfPreviousObservations();\n",
    "plt.figure();\n",
    "plt.scatter(variable, localAIS, marker='o', color='blue', facecolors='none');\n",
    "plt.xlabel('breath rate');\n",
    "plt.ylabel('local AIS (nats)');\n",
    "plt.title('Local AIS versus breath rate');"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1d182f3e-7bba-4364-8d99-262dca9d3e34",
   "metadata": {},
   "source": [
    "8. There is not a one-to-one mapping between breath rate and AIS, since the AIS depends on the past values of breath rate as well. To try to get further insights, run the following code cell:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "df3bc659-8ea0-45a2-88d7-9b05e8dd9311",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 110792 characters]",
      "text/plain": [
       "<Figure size 640x480 with 2 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "[encoded data, 86092 characters]",
      "text/plain": [
       "<Figure size 640x480 with 2 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.figure();\n",
    "plt.scatter(range(1,variable.length+1), variable, c=localAIS, marker='o', facecolors='none', s=8);\n",
    "plt.xlabel('n'); plt.ylabel('breath rate'); plt.colorbar(label='Local AIS (nats)');\n",
    "plt.title('Breath rate versus sample index');\n",
    "\n",
    "plt.figure();\n",
    "plt.scatter(variable[:-1], variable[1:], c=localAIS[1:], marker='o', facecolors='none', s=8);\n",
    "plt.xlabel('breath_rate(n)'); plt.ylabel('breath_rate(n+1)'); plt.colorbar(label='Local AIS (nats)');\n",
    "plt.title('Breath rate next versus previous values');"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "30d08acb-e60a-431a-842d-cb5299f40fc8",
   "metadata": {},
   "source": [
    "  The code in the cell above plots the dynamics in two different ways, highlighting each point in the time series according to their local AIS value. On the first, you could zoom in around an apnoea event to inspect where the higher and lower AIS events occur. The second plot may make this more clear for you: certainly the larger AIS values occur at more extreme breath rates, but compare the AIS values for these extreme breath rates where the previous breath rates were or were not already extreme:\n",
    "  * For the breath rates being in the stable range, what sorts of AIS values do we have? how could we interpret that?\n",
    "  * Which parts of the breath rate dynamics have most potential for information to be added from other source variables?\n",
    "\n",
    "Sample solution plots are available on the tutorial webpage."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "41b0e272-7310-4239-9172-b058ab2a68b6",
   "metadata": {},
   "source": [
    "## B. Transfer entropy from heart rate to breath rate\n",
    "\n",
    "The next step in our information regression of modelling the dynamic updates of the target is to include contributions from an information source: we will examine the information flow from heart rate to breath here:\n",
    "1. Open the TE AutoAnalyser, selecting the `Kraskov (KSG)` estimator, and the `SFI-heartRate_breathVol_bloodOx-extract.txt` data set.\n",
    "2. Leave the source column at 0 (heart rate) and destination column at 1 (breath rate).\n",
    "3. Check the `Add stat. signif.?` button, as we wish to check the statistical significance of our measurements.\n",
    "4. Set `DYN_CORR_EXCL` to 15 as above.\n",
    "5. Leave all of the other parameters as defaults for now, to test that we can run a calculation - press `Generate code and compute`. Hmmm, the result of around 0.0112 bits is much smaller than we measured with MI before and the p-value doesn't look significant. But we haven't properly set the embedding parameters, so let's not draw a hasty conclusion.\n",
    "6. Let's ask JIDT to set embedding parameters for both the heart rate and breath rate automatically for us. We already computed these for the breath rate using AIS, and in code we could read these out as above and include them here; but let's calculate again to confirm we get the same result. We know that both are complex time series, and there are potentially multiple past values of the source contributing to the target so we will embed the source as well as target here:\n",
    "    1. Set `AUTO_EMBED_METHOD` to `MAX_CORR_AIS` to embed both source and target.\n",
    "    2. Set `AUTO_EMBED_METHOD_K_SEARCH_MAX` to 10 to search for an embedding length up to 10.\n",
    "    3. Set `AUTO_EMBED_METHOD_TAU_SEARCH_MAX` to 10 to search for an embedding delay up to 10 (allowing a longer lag than usual since the time series are complex)\n",
    "7. Press `Generate code and compute`:\n",
    "    * Check that the embedding parameters for the breath rate (the target) are the same as determined via the AIS. Note the embedding parameters determined for the heart rate as source.\n",
    "    * Now inspect the value of TE. Ok... that value looks very small and is certainly not significant. Perhaps the past history of the breath rate via AIS is accounting for everything that the heart rate could tell us about the breath?\n",
    "8. Recall however from our previous investigation that we observed a significant MI relationship over lags of 6-12 time steps from the heart rate to breath. Let's compute TE again allowing for such lags. This time - since JIDT has already worked out the embedding lengths - we'll supply the embedding parameters ourselves, and expand the code template to check the TE over various lags. Copy the code generated by the AutoAnalyser into the next code cell, and adjust as follows:\n",
    "    * Change the `DELAY` property to a non-default value, e.g. 2, so that it gets included in our template code.\n",
    "    * Uncheck `Compute result?` then press `Generate code`\n",
    "    * Copy the code in the Python tab and paste into the code cell below.\n",
    "    * Notice that the code has set the optimised embedding parameters to their non-default values (from our calculation before changing the delay), but also still sets up the parameters for auto-embedding. Let's delete the auto-embedding parameters so that it just uses the values we already determined: delete the lines which set `AUTO_EMBED_METHOD`, `AUTO_EMBED_METHOD_K_SEARCH_MAX` and `AUTO_EMBED_METHOD_TAU_SEARCH_MAX`\n",
    "    * Run the code in the cell below once to make sure that it computes a TE with the new source-target lag value that you set.\n",
    "    * Next, add a for-loop, `for delay in range(1,21):` around the lines beginning with setting the `DELAY` property and through to the end of the code. (We could leave the statistical significance calculation out, we're only really interested in it's value for the TE at the largest delay, but let's just leave it there).\n",
    "    * Change the property value supplied for `DELAY` to `str(delay)`.\n",
    "    * Change the `print` statement to include the delay in what is printed out\n",
    "    * Run the code, and observe how the TE changes with the source-target delay."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "ddbdef76-4361-492a-8831-e2ffbea746c6",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "TE_Kraskov (KSG)(delay 1) = -0.0037 nats (null: -0.0013 +/- 0.0108 std dev.; p(surrogate > measured)=0.59000 from 100 surrogates)\n",
      "TE_Kraskov (KSG)(delay 2) = -0.0118 nats (null: -0.0019 +/- 0.0108 std dev.; p(surrogate > measured)=0.84000 from 100 surrogates)\n",
      "TE_Kraskov (KSG)(delay 3) = 0.0089 nats (null: -0.0033 +/- 0.0119 std dev.; p(surrogate > measured)=0.15000 from 100 surrogates)\n",
      "TE_Kraskov (KSG)(delay 4) = 0.0253 nats (null: -0.0025 +/- 0.0117 std dev.; p(surrogate > measured)=0.01000 from 100 surrogates)\n",
      "TE_Kraskov (KSG)(delay 5) = 0.0254 nats (null: 0.0001 +/- 0.0119 std dev.; p(surrogate > measured)=0.01000 from 100 surrogates)\n",
      "TE_Kraskov (KSG)(delay 6) = 0.0546 nats (null: -0.0011 +/- 0.0125 std dev.; p(surrogate > measured)=0.00000 from 100 surrogates)\n",
      "TE_Kraskov (KSG)(delay 7) = 0.0380 nats (null: -0.0015 +/- 0.0122 std dev.; p(surrogate > measured)=0.00000 from 100 surrogates)\n",
      "TE_Kraskov (KSG)(delay 8) = 0.0400 nats (null: -0.0038 +/- 0.0113 std dev.; p(surrogate > measured)=0.00000 from 100 surrogates)\n",
      "TE_Kraskov (KSG)(delay 9) = 0.0368 nats (null: -0.0028 +/- 0.0111 std dev.; p(surrogate > measured)=0.00000 from 100 surrogates)\n",
      "TE_Kraskov (KSG)(delay 10) = 0.0304 nats (null: -0.0009 +/- 0.0113 std dev.; p(surrogate > measured)=0.00000 from 100 surrogates)\n",
      "TE_Kraskov (KSG)(delay 11) = 0.0039 nats (null: -0.0026 +/- 0.0098 std dev.; p(surrogate > measured)=0.26000 from 100 surrogates)\n",
      "TE_Kraskov (KSG)(delay 12) = 0.0207 nats (null: -0.0007 +/- 0.0129 std dev.; p(surrogate > measured)=0.09000 from 100 surrogates)\n",
      "TE_Kraskov (KSG)(delay 13) = 0.0086 nats (null: -0.0020 +/- 0.0123 std dev.; p(surrogate > measured)=0.19000 from 100 surrogates)\n",
      "TE_Kraskov (KSG)(delay 14) = 0.0102 nats (null: -0.0014 +/- 0.0113 std dev.; p(surrogate > measured)=0.13000 from 100 surrogates)\n",
      "TE_Kraskov (KSG)(delay 15) = 0.0154 nats (null: 0.0010 +/- 0.0114 std dev.; p(surrogate > measured)=0.08000 from 100 surrogates)\n",
      "TE_Kraskov (KSG)(delay 16) = 0.0252 nats (null: -0.0017 +/- 0.0114 std dev.; p(surrogate > measured)=0.02000 from 100 surrogates)\n",
      "TE_Kraskov (KSG)(delay 17) = 0.0067 nats (null: -0.0044 +/- 0.0124 std dev.; p(surrogate > measured)=0.16000 from 100 surrogates)\n",
      "TE_Kraskov (KSG)(delay 18) = -0.0005 nats (null: -0.0026 +/- 0.0124 std dev.; p(surrogate > measured)=0.45000 from 100 surrogates)\n",
      "TE_Kraskov (KSG)(delay 19) = 0.0256 nats (null: -0.0010 +/- 0.0117 std dev.; p(surrogate > measured)=0.00000 from 100 surrogates)\n",
      "TE_Kraskov (KSG)(delay 20) = 0.0146 nats (null: -0.0024 +/- 0.0096 std dev.; p(surrogate > measured)=0.06000 from 100 surrogates)\n"
     ]
    }
   ],
   "source": [
    "# Paste the generated code here -- you can skip the import and JVM startup lines and just paste the analysis code:\n",
    "\n",
    "# 0. Load/prepare the data:\n",
    "dataRaw = readFloatsFile.readFloatsFile(\"/home/joseph/JIDT/infodynamics-dist-1.6.1/demos/data/SFI-heartRate_breathVol_bloodOx-extract.txt\")\n",
    "# As numpy array:\n",
    "data = numpy.array(dataRaw)\n",
    "source = JArray(JDouble, 1)(data[:,0].tolist())\n",
    "destination = JArray(JDouble, 1)(data[:,1].tolist())\n",
    "\n",
    "maxDelay = 20;\n",
    "results = numpy.zeros(maxDelay)\n",
    "\n",
    "# 1. Construct the calculator:\n",
    "calcClass = JPackage(\"infodynamics.measures.continuous.kraskov\").TransferEntropyCalculatorKraskov\n",
    "calc = calcClass()\n",
    "# 2. Set any properties to non-default values:\n",
    "calc.setProperty(\"k_HISTORY\", \"2\")\n",
    "calc.setProperty(\"l_HISTORY\", \"2\")\n",
    "for delay in range(1,maxDelay+1):\n",
    "    calc.setProperty(\"DELAY\", str(delay))\n",
    "    calc.setProperty(\"DYN_CORR_EXCL\", \"15\")\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",
    "    result = calc.computeAverageLocalOfObservations()\n",
    "    results[delay-1] = result;\n",
    "    # 6. Compute the (statistical significance via) null distribution empirically (e.g. with 100 permutations):\n",
    "    measDist = calc.computeSignificance(100)\n",
    "    \n",
    "    print(\"TE_Kraskov (KSG)(delay %d) = %.4f nats (null: %.4f +/- %.4f std dev.; p(surrogate > measured)=%.5f from %d surrogates)\" %\\\n",
    "        (delay, result, measDist.getMeanOfDistribution(), measDist.getStdOfDistribution(), measDist.pValue, 100))\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3ce1c881-7736-40ff-b171-b98daf835494",
   "metadata": {},
   "source": [
    "9. Plot the TE versus delay for good measure. You will need to store each result in the for loop above in an array `results` and then plot them against the delay in the code block below (sample plot shown on the tutorial website). Consider:\n",
    "    * Which source-target delay do you think we should use in our best model of the dynamics of the breath rate? So which delay, and therefore what value of TE would we attribute to the heart rate - breath rate interaction?\n",
    "    * Is the p-value for this TE statistically significant against $\\alpha=0.05$? Given that we evaluated TEs for 20 potential source-target delays, how many such TEs might we expect to appear statistically significant just by chance if there was actually no relationship? We will discuss family-wise error correction to handle this in our next module."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "140dc0a9-648b-4599-af3b-237dd7568dcc",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 28912 characters]",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Add code to plot the TE versus delay here\n",
    "\n",
    "plt.figure();\n",
    "plt.scatter(range(1,maxDelay+1), results, marker='x', color='red');\n",
    "plt.xlabel('delay (time units)');\n",
    "plt.ylabel('TE (nats)');\n",
    "plt.title('Heart rate to breath rate TE versus delay');"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2bb754cd-d910-4ed4-82b5-c2199d95c5db",
   "metadata": {},
   "source": [
    "10. Next, we want to check which parts of the breath rate dynamics the heart rate is adding information to:\n",
    "    * We will need to re-run the TE calculation for the optimal delay. Copy the lines of code within the for loop above (except for the stat significance check and `print` line) and paste them in again in the code cell below. Set the variable `delay` to the optimal value you determined before the relevant property is set in these new code lines. And don't overwrite the stored values in `results`.\n",
    "    * Add the following line to compute the local TE values:\n",
    "      ```python\n",
    "      localTEs = calc.computeLocalOfPreviousObservations();\n",
    "      ```\n",
    "    * Run the code cell so that the local TE values are ready to plot and examine"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "1f0d20ec-49bf-4e7a-a52d-bb872f5ed530",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Optimal delay was 6\n"
     ]
    }
   ],
   "source": [
    "# Paste and adapt the lines to re-run the TE calculation for the optimal delay here:\n",
    "\n",
    "# Pull out the optimal delay (we know this to be 6, but this is how to do it automatically in code):\n",
    "delay = numpy.argmax(results) + 1; # Need to add one to the index here\n",
    "print('Optimal delay was %d' % delay)\n",
    "\n",
    "calc.setProperty(\"DELAY\", str(delay))\n",
    "calc.setProperty(\"DYN_CORR_EXCL\", \"15\")\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",
    "result = calc.computeAverageLocalOfObservations()\n",
    "\n",
    "# Compute the local TE values:\n",
    "localTEs = calc.computeLocalOfPreviousObservations();"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6ee40dff-291c-4aa4-9f6e-9082d5697ca6",
   "metadata": {},
   "source": [
    "  * Then we will plot the local TE values against the heart and breath rates in multiple ways to examine where the strong information flows occur. First, let's plot the breath rate transitions (which we examined for AIS) coloured now by TE -- examine the plot code in the following cell (in particular how we plot the local TE for the breath rate at time `n+1` for the breath rate transition from time `n` to time `n+1`) and then run it:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "5e7312de-1ef4-4de1-9086-777f72a4cb45",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 89084 characters]",
      "text/plain": [
       "<Figure size 640x480 with 2 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "[encoded data, 88988 characters]",
      "text/plain": [
       "<Figure size 640x480 with 2 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# First let's colour the transitions in breath rate with the TE values\n",
    "plt.figure();\n",
    "plt.scatter(destination[delay-1:-1], destination[delay:], c=localTEs[delay:], marker='o', facecolors='none', s=8);\n",
    "plt.xlabel('breath_rate(n)'); plt.ylabel('breath_rate(n+1)'); plt.colorbar(label='Local TE (nats)');\n",
    "plt.title('Breath samples coloured by local TE from heart rate');\n",
    "\n",
    "# Make a second plot where we highlight the values outside the central range:\n",
    "plt.figure();\n",
    "plt.scatter(destination[delay-1:-1], destination[delay:], c=localTEs[delay:], marker='o', facecolors='none', s=8);\n",
    "plt.xlabel('breath_rate(n)'); plt.ylabel('breath_rate(n+1)'); plt.colorbar(label='Local TE (nats)');\n",
    "plt.title('Breath samples coloured by local TE from heart rate');\n",
    "plt.clim(-0.5,0.5); # Can use this to highlight values out of the central range"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "af55ecf7-4e4d-4259-896b-000317a705c4",
   "metadata": {},
   "source": [
    "  This plot is more difficult to interpret than the similar plot for AIS we looked at before. You should be able to identify though certain parts of the dynamics where the AIS was not particularly strong where the heart rate is able to add significant information. To gain further insights, instead of plotting the raw TE values, alter the colorbar so that values outside of the central range get extreme colours (e.g. `plt.clim(-0.5,0.5);`). Can you identify general aspects of the dynamics where the source usefully adds information here?\n",
    "  * Next, we'll plot the relationship of previous heart rate to (lagged) breath rate coloured by TE -- inspect and run the following code:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "dcf5b1d4-8530-47c5-938c-3addc171b95e",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 113380 characters]",
      "text/plain": [
       "<Figure size 640x480 with 2 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "[encoded data, 115060 characters]",
      "text/plain": [
       "<Figure size 640x480 with 2 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Next let's make a scatter plot showing the heart rate values:\n",
    "plt.figure();\n",
    "plt.scatter(source[:-delay], destination[delay:], c=localTEs[delay:], marker='o', facecolors='none', s=8);\n",
    "plt.xlabel('heart_rate(n-delay+1)'); plt.ylabel('breath_rate(n+1)'); plt.colorbar(label='Local TE (nats)');\n",
    "plt.title('Heart-breath samples (delay %d) coloured by local TE' % delay);\n",
    "\n",
    "# Make a second plot where we highlight the values outside the central range:\n",
    "plt.figure();\n",
    "plt.scatter(source[:-delay], destination[delay:], c=localTEs[delay:], marker='o', facecolors='none', s=8);\n",
    "plt.xlabel('heart_rate(n-delay+1)'); plt.ylabel('breath_rate(n+1)'); plt.colorbar(label='Local TE (nats)');\n",
    "plt.title('Heart-breath samples (delay %d) coloured by local TE' % delay);\n",
    "plt.clim(-0.5,0.5); # Can use this to highlight values out of the central range"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "18eaa0e3-2463-4af0-97e1-e6d0c86816d2",
   "metadata": {},
   "source": [
    "  Again, it's little more tricky to interpret here, but it can be helped e.g. by altering the colorbar as before (e.g. `plt.clim(-0.5,0.5);`). Can you identify general aspects of the dynamics where the source usefully adds information here?\n",
    "\n",
    "Sample solution plots are available on the tutorial webpage."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7a41e6fe-e804-4c00-bb75-5be92f85eab1",
   "metadata": {},
   "source": [
    "## C. Transfer entropy from blood oxygen concentration to breath rate\n",
    "\n",
    "1. The data set also holds a blood oxygen concentration in the third column. Let's plot it, check its autocorrelation, and also plot it next to the heart and breath rates:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "866ab7c7-4967-4545-844d-c7b4efbaef77",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 75204 characters]",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "[encoded data, 38416 characters]",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "[encoded data, 89368 characters]",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Plot the blood oxygen concentration time series\n",
    "plt.figure(); plt.plot(data[:,2]); plt.xlabel('n'); plt.ylabel('Blood ox. conc.');\n",
    "plt.title('Blood oxygen concentration versus sample index');\n",
    "plt.show()\n",
    "# For the autocorrelation we will use the pandas library:\n",
    "import pandas\n",
    "pandas.plotting.autocorrelation_plot(data[:,2])\n",
    "plt.title('ACF for blood ox. conc. time series');\n",
    "plt.xlim([0,25])\n",
    "# And finally we plot it with the heart and breath rates\n",
    "normedData = (data - data.mean(axis=0)) / data.std(axis=0); # Normalise to make the time series easier to compare\n",
    "plt.figure();\n",
    "plt.plot(normedData[:,0], color='red', label='Heart rate');\n",
    "plt.plot(normedData[:,1], color='blue', label='Breath rate');\n",
    "plt.plot(normedData[:,2], color='green', label='Blood ox. conc.');\n",
    "plt.xlabel('n'); plt.ylabel('Normalised time series values');\n",
    "plt.legend()\n",
    "plt.title('Normalised series versus sample index');"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f7d27944-f3e0-49f4-9415-e70af479ac18",
   "metadata": {},
   "source": [
    "  We can see in the plots that this variable also shows peaks and troughs, with more smooth dynamics, and the troughs coinciding roughly with the apnoea events. It seems that the blood oxygen concentration falls off before the heart rate rises and then the breath rate rises. But can the blood oxygen concentration explain all of what the heart rate adds about the breath rate, does the heart rate add more information, or does the breath rate add other information? We will investigate with the **conditional transfer entropy**.\n",
    "\n",
    "2. First, let's check the embedding parameters for the blood oxygen concentration. Go back to the TE AutoAnalyser as we were using in the previous activity and:\n",
    "    * Make sure we're still set up for the `Kraskov (KSG) estimator`, and the `SFI-heartRate_breathVol_bloodOx-extract.txt` data set.\n",
    "    * Change the source column to 2 (3rd column, but we use Java zero-based indexing in the AutoAnalyser)\n",
    "    * Change `DYN_CORR_EXCL` to 21 due to the longer autocorrelation we have seen on the blood oxygen concentration variable.\n",
    "    * Leave the autoembedding set up as before (e.g. with `AUTO_EMBED_METHOD` to `MAX_CORR_AIS` etc.)\n",
    "    * Run a calculation here (you may need to check `Compute Result?` again) -- we just want to see the embedding that is calculated for the blood oxygen time series. Observe the embedding parameters computed in `l_HISTORY` and `l_TAU`. We're not interested in the result, because we have not optimised the source-target lag here.\n",
    "    * Press `Generate code and Compute` a second time so that the new embedding is added to the properties in the generated code.\n",
    "3. Now we need to compute the optimal lag. Make a copy of your code from the TE investigation of heart$\\rightarrow$breath rate interaction in the activity above (copy _all_ the code cells in activity B above, in order, into the code cell below). We are going to alter it to now show the TE from the blood oxygenation:\n",
    "    * Change the source variable to load column 2 of the data instead of column 0.\n",
    "    * From your AutoAnalyser, copy the code to set the parameters `l_HISTORY` and `l_TAU` and paste these over the top of where the properties for `l_HISTORY` and potentially `l_TAU` are first set early in the code below.\n",
    "    * Change the value set for `DYN_CORR_EXCL` to `21` (this will be in two places in your code potentially)\n",
    "    * Change any references to heart rate in your comments or figure captions to blood oxygen.\n",
    "    * Note that after the for loop, the delay may have been hard-coded to the value we found to be optimal for the heart rate. You can either run once to determine the optimal value for our new source, then go back and set it properly; or use `delay = numpy.argmax(results) + 1;` to automatically set it (if you didn't do that already)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "16924f9e-500e-4343-9ce5-01ac8292ced9",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "TE_Kraskov (KSG)(delay 1) = 0.0360 nats (null: -0.0037 +/- 0.0139 std dev.; p(surrogate > measured)=0.00000 from 100 surrogates)\n",
      "TE_Kraskov (KSG)(delay 2) = 0.0211 nats (null: -0.0040 +/- 0.0115 std dev.; p(surrogate > measured)=0.03000 from 100 surrogates)\n",
      "TE_Kraskov (KSG)(delay 3) = 0.0306 nats (null: -0.0006 +/- 0.0119 std dev.; p(surrogate > measured)=0.00000 from 100 surrogates)\n",
      "TE_Kraskov (KSG)(delay 4) = 0.0320 nats (null: -0.0015 +/- 0.0130 std dev.; p(surrogate > measured)=0.01000 from 100 surrogates)\n",
      "TE_Kraskov (KSG)(delay 5) = 0.0346 nats (null: -0.0026 +/- 0.0108 std dev.; p(surrogate > measured)=0.00000 from 100 surrogates)\n",
      "TE_Kraskov (KSG)(delay 6) = 0.0256 nats (null: -0.0023 +/- 0.0108 std dev.; p(surrogate > measured)=0.00000 from 100 surrogates)\n",
      "TE_Kraskov (KSG)(delay 7) = 0.0265 nats (null: -0.0029 +/- 0.0123 std dev.; p(surrogate > measured)=0.00000 from 100 surrogates)\n",
      "TE_Kraskov (KSG)(delay 8) = 0.0265 nats (null: -0.0021 +/- 0.0119 std dev.; p(surrogate > measured)=0.00000 from 100 surrogates)\n",
      "TE_Kraskov (KSG)(delay 9) = 0.0236 nats (null: -0.0021 +/- 0.0115 std dev.; p(surrogate > measured)=0.02000 from 100 surrogates)\n",
      "TE_Kraskov (KSG)(delay 10) = 0.0246 nats (null: -0.0022 +/- 0.0112 std dev.; p(surrogate > measured)=0.02000 from 100 surrogates)\n",
      "TE_Kraskov (KSG)(delay 11) = 0.0159 nats (null: -0.0008 +/- 0.0114 std dev.; p(surrogate > measured)=0.05000 from 100 surrogates)\n",
      "TE_Kraskov (KSG)(delay 12) = 0.0283 nats (null: -0.0021 +/- 0.0114 std dev.; p(surrogate > measured)=0.00000 from 100 surrogates)\n",
      "TE_Kraskov (KSG)(delay 13) = 0.0196 nats (null: -0.0008 +/- 0.0116 std dev.; p(surrogate > measured)=0.03000 from 100 surrogates)\n",
      "TE_Kraskov (KSG)(delay 14) = 0.0125 nats (null: -0.0020 +/- 0.0128 std dev.; p(surrogate > measured)=0.11000 from 100 surrogates)\n",
      "TE_Kraskov (KSG)(delay 15) = 0.0108 nats (null: -0.0029 +/- 0.0133 std dev.; p(surrogate > measured)=0.16000 from 100 surrogates)\n",
      "TE_Kraskov (KSG)(delay 16) = 0.0025 nats (null: -0.0005 +/- 0.0111 std dev.; p(surrogate > measured)=0.40000 from 100 surrogates)\n",
      "TE_Kraskov (KSG)(delay 17) = 0.0126 nats (null: -0.0013 +/- 0.0129 std dev.; p(surrogate > measured)=0.15000 from 100 surrogates)\n",
      "TE_Kraskov (KSG)(delay 18) = 0.0196 nats (null: -0.0019 +/- 0.0130 std dev.; p(surrogate > measured)=0.05000 from 100 surrogates)\n",
      "TE_Kraskov (KSG)(delay 19) = 0.0132 nats (null: -0.0013 +/- 0.0116 std dev.; p(surrogate > measured)=0.11000 from 100 surrogates)\n",
      "TE_Kraskov (KSG)(delay 20) = 0.0200 nats (null: -0.0020 +/- 0.0118 std dev.; p(surrogate > measured)=0.04000 from 100 surrogates)\n",
      "Optimal delay was 1\n"
     ]
    },
    {
     "data": {
      "image/png": "[encoded data, 30512 characters]",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "[encoded data, 89020 characters]",
      "text/plain": [
       "<Figure size 640x480 with 2 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "[encoded data, 88120 characters]",
      "text/plain": [
       "<Figure size 640x480 with 2 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "[encoded data, 113572 characters]",
      "text/plain": [
       "<Figure size 640x480 with 2 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "[encoded data, 114352 characters]",
      "text/plain": [
       "<Figure size 640x480 with 2 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Copy your heart->breath TE code and adapt it below:\n",
    "\n",
    "# 0. Load/prepare the data:\n",
    "dataRaw = readFloatsFile.readFloatsFile(\"/home/joseph/JIDT/infodynamics-dist-1.6.1/demos/data/SFI-heartRate_breathVol_bloodOx-extract.txt\")\n",
    "# As numpy array:\n",
    "data = numpy.array(dataRaw)\n",
    "source = JArray(JDouble, 1)(data[:,2].tolist())\n",
    "destination = JArray(JDouble, 1)(data[:,1].tolist())\n",
    "\n",
    "maxDelay = 20;\n",
    "results = numpy.zeros(maxDelay)\n",
    "\n",
    "# 1. Construct the calculator:\n",
    "calcClass = JPackage(\"infodynamics.measures.continuous.kraskov\").TransferEntropyCalculatorKraskov\n",
    "calc = calcClass()\n",
    "# 2. Set any properties to non-default values:\n",
    "calc.setProperty(\"k_HISTORY\", \"2\")\n",
    "calc.setProperty(\"l_HISTORY\", \"2\")\n",
    "calc.setProperty(\"l_TAU\", \"3\")\n",
    "for delay in range(1,maxDelay+1):\n",
    "    calc.setProperty(\"DELAY\", str(delay))\n",
    "    calc.setProperty(\"DYN_CORR_EXCL\", \"21\")\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",
    "    result = calc.computeAverageLocalOfObservations()\n",
    "    results[delay-1] = result;\n",
    "    # 6. Compute the (statistical significance via) null distribution empirically (e.g. with 100 permutations):\n",
    "    measDist = calc.computeSignificance(100)\n",
    "    \n",
    "    print(\"TE_Kraskov (KSG)(delay %d) = %.4f nats (null: %.4f +/- %.4f std dev.; p(surrogate > measured)=%.5f from %d surrogates)\" %\\\n",
    "        (delay, result, measDist.getMeanOfDistribution(), measDist.getStdOfDistribution(), measDist.pValue, 100))\n",
    "\n",
    "plt.figure();\n",
    "plt.scatter(range(1,maxDelay+1), results, marker='x', color='red');\n",
    "plt.xlabel('delay (time units)');\n",
    "plt.ylabel('TE (nats)');\n",
    "plt.title('Blood ox. conc. to breath rate TE versus delay');\n",
    "\n",
    "# Pull out the optimal delay (we know this to be 5, but this is how to do it automatically in code):\n",
    "delay = numpy.argmax(results) + 1; # Need to add one to the index here\n",
    "print('Optimal delay was %d' % delay)\n",
    "\n",
    "calc.setProperty(\"DELAY\", str(delay))\n",
    "calc.setProperty(\"DYN_CORR_EXCL\", \"21\")\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",
    "result = calc.computeAverageLocalOfObservations()\n",
    "\n",
    "# Compute the local TE values:\n",
    "localTEs = calc.computeLocalOfPreviousObservations();\n",
    "\n",
    "# First let's colour the transitions in breath rate with the TE values\n",
    "plt.figure();\n",
    "plt.scatter(destination[delay-1:-1], destination[delay:], c=localTEs[delay:], marker='o', facecolors='none', s=8);\n",
    "plt.xlabel('breath_rate(n)'); plt.ylabel('breath_rate(n+1)'); plt.colorbar(label='Local TE (nats)');\n",
    "plt.title('Breath samples coloured by local TE from blood ox.');\n",
    "\n",
    "# Make a second plot where we highlight the values outside the central range:\n",
    "plt.figure();\n",
    "plt.scatter(destination[delay-1:-1], destination[delay:], c=localTEs[delay:], marker='o', facecolors='none', s=8);\n",
    "plt.xlabel('breath_rate(n)'); plt.ylabel('breath_rate(n+1)'); plt.colorbar(label='Local TE (nats)');\n",
    "plt.title('Breath samples coloured by local TE from blood ox.');\n",
    "plt.clim(-0.5,0.5); # Can use this to highlight values out of the central range\n",
    "\n",
    "# Next let's make a scatter plot showing the blood ox. values:\n",
    "plt.figure();\n",
    "plt.scatter(source[:-delay], destination[delay:], c=localTEs[delay:], marker='o', facecolors='none', s=8);\n",
    "plt.xlabel('blood_ox(n-delay+1)'); plt.ylabel('breath_rate(n+1)'); plt.colorbar(label='Local TE (nats)');\n",
    "plt.title('BloodOx-breath samples (delay %d) coloured by local TE' % delay);\n",
    "\n",
    "# Make a second plot where we highlight the values outside the central range:\n",
    "plt.figure();\n",
    "plt.scatter(source[:-delay], destination[delay:], c=localTEs[delay:], marker='o', facecolors='none', s=8);\n",
    "plt.xlabel('blood_ox(n-delay+1)'); plt.ylabel('breath_rate(n+1)'); plt.colorbar(label='Local TE (nats)');\n",
    "plt.title('BloodOx-breath samples (delay %d) coloured by local TE' % delay);\n",
    "plt.clim(-0.5,0.5); # Can use this to highlight values out of the central range"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "de2213ee-32d5-49b3-a890-e4061bff24a4",
   "metadata": {},
   "source": [
    "5. If you didn't automate setting the `delay` variable to the optimal value, then go back to the code and change the `delay` variable after the loop to the optimal value; and then run the code again. Now, examine the scatter plots coloured by local TE values to see which parts of the dynamics are well explained by the blood oxygen concentration. (You will certainly find it useful to alter the colorbar limits as before, `plt.clim(-0.5,0.5);`). Are these the same or different than we saw for the heart rate? What further insights does this give you about how conditioning on blood oxygenation may change the TE from heart rate to breath rate?\n",
    "\n",
    "Sample solution plots are available on the tutorial webpage."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b4a4b239-834c-4d91-b071-6b8d89ad49ee",
   "metadata": {},
   "source": [
    "## D. Conditional transfer entropies\n",
    "\n",
    "Now that we have determined two significant relationships from heart rate and blood oxygen concentration to breath rate, and all of their embedding parameters, let's examine conditional transfer entropies from each of them to the breath rate.\n",
    "1. Open a new instance of the AutoAnalyser and select the Conditional Transfer Entropy.\n",
    "2. Select the `Kraskov (KSG)` estimator, and the `SFI-heartRate_breathVol_bloodOx-extract.txt` data set as before.\n",
    "3. Leave the source column at 0 (heart rate), destination column at 1 (breath rate) and conditional column at 2 (blood oxygen concentration).\n",
    "4. Check the `Add stat. signif.?` button, as we wish to check the statistical significance of our measurements.\n",
    "5. Set `DYN_CORR_EXCL` to 21 as above.\n",
    "6. Set the embedding variables for the breath rate (`k_HISTORY` and `k_TAU`) and heart rate (`l_HISTORY` and `l_TAU`) to any non-default values that were determined for their pairwise TE above.\n",
    "7. Set the `DELAY` to the optimised value determined for the pairwise TE between heart rate and breath rate above.\n",
    "8. Now set the embedding values for the blood oxygen concentration as the conditional:\n",
    "    * Set `COND_EMBED_LENGTHS` to the value determined for `l_HISTORY` in the blood oxygen to breath rate analysis.\n",
    "    * Set `COND_TAUS` to the value determined for `l_TAU` in the blood oxygen to breath rate analysis.\n",
    "    * Set `COND_DELAYS` to the value determined for `delay` in the blood oxygen to breath rate analysis.\n",
    "9. Now press `Generate code and Compute` and inspect the result. Is there a significant conditional TE here? What can you determine from that:\n",
    "    * Would the heart rate be usefully added to a model of the dynamics of breath rate, if blood oxygenation was already included in the model? Is this what you expected from above? (In other words, does the heart rate add information that blood oxygen concentration does not alone -- and could we see clues about this from our investigations above?)\n",
    "    * Does it mean there is or is not a directed relationship from the heart rate to breath rate alone?\n",
    "    * What does it mean for how we would consider both sources? If you had to choose one source only in a model of the breath rate, which would it be and why? Would you use both sources? (Can you conclude fully on this without running another calculation?)\n",
    "    * If the result was not (or had not been) significant, do you think the result might change if we had more data?\n",
    "10. (_Optional challenge_) Can you investigate the local values of TE from heart rate to breath rate, conditioned on blood oxygenation, to see which effects persisted and which did not? You would take the code generated from the previous step, paste this into the code cell below, then copy and paste the code from the pairwise TE investigation earlier which generated the local TE values and then made the plots with them. You would want to compare the plots made here for conditional TE to those for pairwise TE before. You would find the `plt.clim(-0.5,0.5);` command useful as above."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "id": "7ded30e6-df55-42a5-b607-63a9d241a5d5",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "CTE_Kraskov (KSG)(heart -> breath | bloodOx) = 0.0236 nats (null: -0.0020 +/- 0.0103 std dev.; p(surrogate > measured)=0.00000 from 100 surrogates)\n"
     ]
    },
    {
     "data": {
      "image/png": "[encoded data, 88528 characters]",
      "text/plain": [
       "<Figure size 640x480 with 2 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "[encoded data, 90012 characters]",
      "text/plain": [
       "<Figure size 640x480 with 2 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "[encoded data, 112804 characters]",
      "text/plain": [
       "<Figure size 640x480 with 2 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "[encoded data, 116364 characters]",
      "text/plain": [
       "<Figure size 640x480 with 2 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Add code for the optional challenge task of plotting local conditional TE values:\n",
    "\n",
    "# 0. Load/prepare the data:\n",
    "dataRaw = readFloatsFile.readFloatsFile(\"/home/joseph/JIDT/infodynamics-dist-1.6.1/demos/data/SFI-heartRate_breathVol_bloodOx-extract.txt\")\n",
    "# As numpy array:\n",
    "data = numpy.array(dataRaw)\n",
    "source = JArray(JDouble, 1)(data[:,0].tolist())\n",
    "destination = JArray(JDouble, 1)(data[:,1].tolist())\n",
    "conditional = JArray(JDouble, 1)(data[:,2].tolist())\n",
    "\n",
    "# 1. Construct the calculator:\n",
    "calcClass = JPackage(\"infodynamics.measures.continuous.kraskov\").ConditionalTransferEntropyCalculatorKraskov\n",
    "calc = calcClass()\n",
    "# 2. Set any properties to non-default values:\n",
    "calc.setProperty(\"k_HISTORY\", \"2\")\n",
    "calc.setProperty(\"l_HISTORY\", \"2\")\n",
    "calc.setProperty(\"DELAY\", \"6\")\n",
    "calc.setProperty(\"COND_EMBED_LENGTHS\", \"2\")\n",
    "calc.setProperty(\"COND_TAUS\", \"3\")\n",
    "calc.setProperty(\"DYN_CORR_EXCL\", \"21\")\n",
    "# 3. Initialise the calculator for (re-)use:\n",
    "calc.initialise()\n",
    "# 4. Supply the sample data:\n",
    "calc.setObservations(source, destination, conditional)\n",
    "# 5. Compute the estimate:\n",
    "result = calc.computeAverageLocalOfObservations()\n",
    "# 6. Compute the (statistical significance via) null distribution empirically (e.g. with 100 permutations):\n",
    "measDist = calc.computeSignificance(100)\n",
    "\n",
    "print(\"CTE_Kraskov (KSG)(heart -> breath | bloodOx) = %.4f nats (null: %.4f +/- %.4f std dev.; p(surrogate > measured)=%.5f from %d surrogates)\" %\\\n",
    "    (result, measDist.getMeanOfDistribution(), measDist.getStdOfDistribution(), measDist.pValue, 100))\n",
    "\n",
    "# Compute the local TE values:\n",
    "localCondTEs = calc.computeLocalOfPreviousObservations();\n",
    "\n",
    "# First let's colour the transitions in breath rate with the local conditional TE values\n",
    "delay = int(calc.getProperty('DELAY')); # Instead of hard-coding this in case of changes, we'll work it out\n",
    "plt.figure();\n",
    "plt.scatter(destination[delay-1:-1], destination[delay:], c=localCondTEs[delay:], marker='o', facecolors='none', s=8);\n",
    "plt.xlabel('breath_rate(n)'); plt.ylabel('breath_rate(n+1)'); plt.colorbar(label='Local CTE (nats)');\n",
    "plt.title('Breath samples coloured by local cond TE from heart rate');\n",
    "\n",
    "# Make a second plot where we highlight the values outside the central range:\n",
    "plt.figure();\n",
    "plt.scatter(destination[delay-1:-1], destination[delay:], c=localCondTEs[delay:], marker='o', facecolors='none', s=8);\n",
    "plt.xlabel('breath_rate(n)'); plt.ylabel('breath_rate(n+1)'); plt.colorbar(label='Local CTE (nats)');\n",
    "plt.title('Breath samples coloured by local cond TE from heart rate');\n",
    "plt.clim(-0.5,0.5); # Can use this to highlight values out of the central range\n",
    "\n",
    "# Next let's make a scatter plot showing the heart values:\n",
    "plt.figure();\n",
    "plt.scatter(source[:-delay], destination[delay:], c=localCondTEs[delay:], marker='o', facecolors='none', s=8);\n",
    "plt.xlabel('heart_rate(n-delay+1)'); plt.ylabel('breath_rate(n+1)'); plt.colorbar(label='Local TE (nats)');\n",
    "plt.title('Heart-breath samples (delay %d) coloured by local cond TE' % delay);\n",
    "\n",
    "# Make a second plot where we highlight the values outside the central range:\n",
    "plt.figure();\n",
    "plt.scatter(source[:-delay], destination[delay:], c=localCondTEs[delay:], marker='o', facecolors='none', s=8);\n",
    "plt.xlabel('heart_rate(n-delay+1)'); plt.ylabel('breath_rate(n+1)'); plt.colorbar(label='Local TE (nats)');\n",
    "plt.title('Heart-breath samples (delay %d) coloured by local cond TE' % delay);\n",
    "plt.clim(-0.5,0.5); # Can use this to highlight values out of the central range"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f6d98e67-4b12-4597-afb8-003967342987",
   "metadata": {},
   "source": [
    "11. Let's also check whether the blood oxygenation adds any information that the heart rate does not by computing TE from blood oxygenation conditioned on heart rate. Continuing in the AutoAnalyser with the values as currently set, now:\n",
    "    * Swap the source and conditional column values\n",
    "    * Swap the values of `l_HISTORY` and `COND_EMBED_LENGTHS`\n",
    "    * Swap the values of `l_TAU` and `COND_TAUS`\n",
    "    * Swap the values of `DELAY` and `COND_DELAYS`\n",
    "12. Press `Generate code and Compute` - observe the answer - do we have a statistically significant result? What can we conclude about whether the blood oxygenation would be usefully added to a model of the dynamics of breath rate, given that heart rate was already included in the model?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b054ed20-5944-4768-a704-107d2e50fcd9",
   "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
}
