{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Module 3 - What is Information? I\n",
    "\n",
    "Author: Julio Correa, 2020; based on the original Matlab tutorials.<br/>\n",
    "Adaptations by: J. Lizier, 2023-\n",
    "\n",
    "The following block aims to import all the relevant libraries to analyse data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import math"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Preparing your environment\n",
    "\n",
    "As per the previous module, we want to use functions we have defined in our previous work in other notebooks.\n",
    "\n",
    "You have several options on handling this as per last week. I suggest you update `simpleinfotheory.py` script to **add the new functions you wrote in the last module**, and import the required functions from this."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Option 1: your notebook from Module 1 is complete:\n",
    "# from ipynb.fs.full.Module_1_notebook import entropy\n",
    "# Option 2: you use the Module 1 note book solutions: (if so, ignore the out\n",
    "# from ipynb.fs.full.Module_1_notebook_solutions import entropy\n",
    "# Option 3: edit simpleinfotheory.py and past your functions into that as you write them\n",
    "from simpleinfotheory import entropy, entropyempirical, jointentropy, jointentropyempirical, conditionalentropy, conditionalentropyempirical"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 4. Coding mutual information\n",
    "\n",
    "In this exercise we continue to alter the Python code to measure the mutual information between variables x and y for a distribution p(x,y):\n",
    "\n",
    "$I\\left(X;Y\\right)=H\\left(X\\right)+H\\left(Y\\right)-H\\left(X,Y\\right)$\n",
    "\n",
    "Your task is to edit the Python function `mutualinformation(p)` in the next cell to return the mutual information for the given distribution $p(x,y)$ over joint outcomes $\\{x,y\\}$ of variables $X,Y$.\n",
    "\n",
    "As before, the input argument to the function is a matrix `p`, representing the probability mass for each joint outcome of $\\{x,y\\}$. That is, `p` is a matrix with the $(i,j)$th entry in the matrix giving the probability for the joint outcome of the $i$th value that $x$ may take along with the $j$th value that $y$ may take. The sum of the items in the matrix `p` must be 1.\n",
    "For example, for the probability table:\n",
    "<table style=\"margin: 20px;\">\n",
    "    <tbody>\n",
    "        <tr style=\"border-bottom: solid;\">\n",
    "            <td style=\"margin: 20px;\">p(x,y)</td>\n",
    "            <td style=\"border-left: solid; margin: 20px;\">y=0</td>\n",
    "            <td style=\"margin: 20px;\">y=1</td>\n",
    "        </tr>\n",
    "        <tr>\n",
    "            <td>x=0</td>\n",
    "            <td style=\"border-left: solid;\">0.2</td>\n",
    "            <td>0.3</td>\n",
    "        </tr>\n",
    "        <tr>\n",
    "            <td>x=1</td>\n",
    "            <td style=\"border-left: solid;\">0.1</td>\n",
    "            <td>0.4</td>\n",
    "        </tr>\n",
    "    </tbody>\n",
    "</table>\n",
    "where we have a binary $x$ and $y$ we have <code>p=np.array([[0.2, 0.3],[0.1,0.4]])</code> where $p(x=0,y=0) = 0.2$, $p(x=0,y=1) = 0.3$, $p(x=1,y=0) = 0.1$, and $p(x=1,y=1) = 0.4$. If the variable $x$ can take more than two values for example, then we will have more than two rows in <code>p</code> (e.g. <code>p = np.array([[0.15, 0.1], [0.1, 0.3], [0.15, 0.2]])</code>).\n",
    "\n",
    "1. To fill in the template, you will need to call your existing functions `jointentropy(p)` for $H(X,Y)$ and `entropy(p)` for $H(X)$ and $H(Y)$ to provide the calculations needed. Note that to compute $H(Y)$ you will need to extract $p(y)$ from the $p(x,y)$ matrix by summing over all $x$ rows (as per the activity for conditional entropy in the previous module), whilst for $H(X)$ you will need to extract $p(x)$ from the $p(x,y)$ matrix by summing over all $y$ columns."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "\"\"\"function mutualinformation(p)\n",
    "Computes the mutual information over all outcomes x of a random\n",
    "variable X with outcomes y of a random variable Y.\n",
    "Probability matrix p(x,y) is given for each candidate outcome\n",
    "(x,y).\n",
    "\n",
    "Inputs:\n",
    "- p - 2D probability distribution function over all outcomes (x,y).\n",
    "   p is a numpy matrix (or list of lists) over all combinations of x and y,\n",
    "where p[0,2] gives the probability of the first symbol of variable\n",
    "x co-occuring with the third symbol of variable y.\n",
    "   E.g. p = np.array([[0.2, 0.3], [0.1, 0.4]]).\n",
    "   The sum over p must be 1.\n",
    "\n",
    "Outputs:\n",
    "- result - mutual information of X with Y\n",
    "\n",
    "Copyright (C) 2020-, Julio Correa, Joseph T. Lizier\n",
    "Distributed under GNU General Public License v3\n",
    "\"\"\"\n",
    "def mutualinformation(p):\n",
    "    \n",
    "    # First make sure the array is now a numpy array\n",
    "    if type(p) != np.array:\n",
    "        p = np.array(p)\n",
    "\n",
    "    # Should we check any potential error conditions on the input?\n",
    "    # a. Should we check p is a matrix, not a vector?\n",
    "    # Actually we won't since a vector would be valid if one variable only ever took one value.\n",
    "    # b. Check that the probabilities normalise to 1:\n",
    "    if (abs(np.sum(p) - 1) > 0.00001):\n",
    "        raise Exception(\"Probability distribution must sum to 1: sum is %.4f\" % np.sum(p))\n",
    "\n",
    "    # We need to compute H(X) + H(Y) - H(X,Y):\n",
    "    # 1. joint entropy:\n",
    "    H_XY = jointentropy(p)\n",
    "\n",
    "    # 2. marginal entropy of X:\n",
    "    # But how to get p_x???\n",
    "    p_x = p.sum(axis=1); # Since x changes along the rows, summing over the y's (dimension 1 argument in the sum) will just return p(x)\n",
    "    H_X = entropy(p_x);\n",
    "\n",
    "    # 2. marginal entropy of Y:\n",
    "    # But how to get p_y???\n",
    "    p_y = p.sum(axis=0); # Since y changes along the columns, summing over the x's (dimension 0 argument in the sum) will just return p(y)\n",
    "    H_Y = entropy(p_y);\n",
    "\n",
    "    result = H_X + H_Y - H_XY\n",
    "    \n",
    "    return result"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "2. Test that your code works, e.g. by running:\n",
    "    1. `mutualinformation(np.array([[0.2, 0.3],[ 0.1, 0.4]]))` and validating that you get the result 0.0349 bits. Recall that the conditional entropy computed for this $p(x,y)$ probability table in the previous module was 0.965 bits - is the result for MI sensible with respect to that result?<br/>\n",
    "    Confirm that MI is symmetric with respect to the input variables by computing `mutualinformation(np.array([[0.2, 0.1], [0.3, 0.4]]))`\n",
    "    1. `mutualinformation(np.array([[0.5, 0],[ 0, 0.5]]))` and validating that you get the result 1 bit.\n",
    "    1. `mutualinformation(np.array([[0.25, 0.25],[ 0.25, 0.25]]))` and validating that you get the result 0 bits. Can you explain this and the previous result?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.034851554559677256\n",
      "1.0\n",
      "0.0\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/home/joseph/JIDT/course/course/Module01-Entropy/PythonSimpleFunctions/completed/simpleinfotheory.py:23: RuntimeWarning: divide by zero encountered in log2\n",
      "  return -np.log2(p)\n",
      "/home/joseph/JIDT/course/course/Module01-Entropy/PythonSimpleFunctions/completed/simpleinfotheory.py:52: RuntimeWarning: invalid value encountered in multiply\n",
      "  weightedShannonInfos = p*(infocontent(p))\n"
     ]
    }
   ],
   "source": [
    "# Test the code here:\n",
    "print( mutualinformation(np.array([[0.2, 0.3],[ 0.1, 0.4]])) )\n",
    "print( mutualinformation(np.array([[0.5, 0],[ 0, 0.5]])) )\n",
    "print( mutualinformation(np.array([[0.25, 0.25],[ 0.25, 0.25]])) )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "3. Coming back to the Guess Who? example using the [Kooky character sheet](https://web.archive.org/web/20170215034006/http://www.hasbro.com/upload/guesswho/GWc_Kooky-en_GB.pdf):\n",
    "    1. Compute the mutual information between whether the character has horns and whether they have eyebrows, i.e. $I(horns ; eyebrows)$? As per the exercise in the previous module, construct first the table $p(horns, eyebrows)$ for all 4 combinations of these two binary variables, then pass this to your function.\n",
    "    1. This next question requires you to have watched the next video lecture on Pointwise Mutual Information first.<br/>\n",
    "    From your constructed table $p(horns, eyebrows)$, first construct the distribution $p(eyebrows)$ and evaluate $p(eyebrows = true)$. Next, construct the conditional probability distribution $p(eyebrows | horns)$, and then evaluate the conditional probabilities given that the character has horns, i.e. $p(eyebrows | horns = true)$. Finally, evaluate $p(eyebrows = true | horns = true)$, and then use $p(eyebrows = true | horns = true)$ and $p(eyebrows = true)$ to compute $i(eyebrows = true ; horns = true)$. Use your result to explain how helpful or unhelpful knowing that the character has horns, i.e. $horns = true$, was in determining whether the character had eyebrows, i.e. $eyebrows = true$.\n",
    "    1. Is $I(eyebrows ; horns)$ the same as $I(horns ; eyebrows)$? Recall that $H(horns | eyebrows) != H(eyebrows | horns)$, so explain why the result for MI is the same or different? \n",
    "    1. _Challenge_: Can you find a pair of traits (or pairs of sets of traits) that appear to have high mutual information? What does it mean for these traits to have high mutual information in the individuals in this sheet?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "I(horns ; eyebrows) = 0.0153 bits\n",
      "p(eyebrows=true)=0.4583\n",
      "Pointwise info from horns = true to eyebrows = true is 0.3886 bits\n",
      "I(eyebrows ; horns) = 0.0153 bits\n"
     ]
    }
   ],
   "source": [
    "# Construct the table p(horns,eyebrows) -- or grab it from last week\n",
    "#                  [[h_0_e_0, h_0_e_1], [h_1_e_0, h_1_e_1]]\n",
    "p_horns_eyebrows = np.array([[11, 8], [2, 3]]) / 24\n",
    "# Compute I(horns ; eyebrows)\n",
    "print(\"I(horns ; eyebrows) = %.4f bits\" % mutualinformation(p_horns_eyebrows))\n",
    "\n",
    "# Compute the probabilities for the pointwise mutual information calculations\n",
    "p_eyebrows = p_horns_eyebrows.sum(axis=0) # Since eyebrows changes along the columns, summing over the horns's (dimension 0 argument in the sum) will just return p(eyebrows)\n",
    "p_eyebrows_true = p_eyebrows[1]\n",
    "print(\"p(eyebrows=true)=%.4f\" % p_eyebrows_true)\n",
    "p_horns = p_horns_eyebrows.sum(axis=1) # Since horns changes along the rows, summing over the eyebrows's (dimension 1 argument in the sum) will just return p(horns)\n",
    "p_eyebrows_given_horns = [p_horns_eyebrows[0,:]/p_horns[0], p_horns_eyebrows[1,:]/p_horns[1]]\n",
    "p_eyebrows_given_horns_true = p_horns_eyebrows[1,:]/p_horns[1]\n",
    "p_eyebrows_true_given_horns_true = p_eyebrows_given_horns_true[1]\n",
    "i_e_1_h_1 = np.log2( p_eyebrows_true_given_horns_true / p_eyebrows_true )\n",
    "print('Pointwise info from horns = true to eyebrows = true is %.4f bits' % i_e_1_h_1)\n",
    "\n",
    "# Compute I(eyebrows ; horns) should be exactly the same\n",
    "print(\"I(eyebrows ; horns) = %.4f bits\" % mutualinformation(p_horns_eyebrows.T))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "4. Finally, let's code mutual information $I(X;Y)$ for empirical samples `xn` and `yn` in the cell below.<br/>\n",
    "_Hint_: You can call your existing code `jointentropyempirical()` to compute $H(X,Y)$, $H(X)$ and $H(Y)$ respectively, by passing in `[xn,yn]`, `xn` and `yn` as arguments to these functions respectively."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "\"\"\"function mutualinformationempirical(xn,yn)\n",
    "Computes the mutual information over all samples xn of a random\n",
    "variable X with samples yn of a random variable Y.\n",
    "\n",
    "Inputs:\n",
    "- xn - numpy matrix of samples of outcomes x. May be a 1D vector of samples, or\n",
    "    a 2D matrix, where each row is a vector sample for a multivariate X.\n",
    "- yn - numpy matrix of samples of outcomes x. May be a 1D vector of samples, or\n",
    "    a 2D matrix, where each row is a vector sample for a multivariate Y.\n",
    "    Must have the same number of rows as X.\n",
    "\n",
    "Outputs:\n",
    "- result - mutual information of X with Y\n",
    "\n",
    "Copyright (C) 2020-, Julio Correa, Joseph T. Lizier\n",
    "Distributed under GNU General Public License v3\n",
    "\"\"\"\n",
    "def mutualinformationempirical(xn,yn):\n",
    "    \n",
    "    # First, error checking, and converting argument into standard form:    \n",
    "    xn = np.array(xn)\n",
    "    # Convert to column vectors if not already:\n",
    "    if xn.ndim == 1:\n",
    "        xn = np.reshape(xn,(len(xn),1))\n",
    "    yn = np.array(yn)\n",
    "    if yn.ndim == 1:\n",
    "        yn = np.reshape(yn,(len(yn),1))\n",
    "    [rx,cx] = xn.shape\n",
    "    [ry,cy] = yn.shape\n",
    "\n",
    "    # Should we check any potential error conditions on the input?\n",
    "    # Check that their number of rows are the same:\n",
    "    assert(rx == ry)\n",
    "\n",
    "    # We need to compute H(X) + H(Y) - H(X,Y):\n",
    "    # 1. joint entropy:\n",
    "    (H_XY, xySymbols, xyProbs) = jointentropyempirical(xn, yn); # How to compute this empirically ...?\n",
    "    # 2. marginal entropy of Y: (call 'joint' in case yn is multivariate)\n",
    "    (H_Y, ySymbols, yProbs) = jointentropyempirical(yn)\n",
    "    # 3. marginal entropy of X: (call 'joint' in case xn is multivariate)\n",
    "    (H_X, xSymbols, xProbs) = jointentropyempirical(xn);\n",
    "\t\n",
    "    result = H_X + H_Y - H_XY;\n",
    "    return result"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "5. Test that your code works by running, e.g.:\n",
    "    1. `mutualinformationempirical([0,0,1,1],[0,1,0,1])` and validating that you get the result 0 bits.\n",
    "    1. `mutualinformationempirical([0,0,1,1],[0,0,1,1])` and validating that you get the result 1 bit.\n",
    "    1. Can you explain the expected results for these boundary cases?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.0\n",
      "1.0\n"
     ]
    }
   ],
   "source": [
    "# Test the code here:\n",
    "print( mutualinformationempirical([0,0,1,1],[0,1,0,1]) )\n",
    "print( mutualinformationempirical([0,0,1,1],[0,0,1,1]) )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "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": 4
}
