{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "# Standard numpy for matrices, vectors, etc..\n",
    "import numpy as np\n",
    "\n",
    "# Visualisation (plotting, etc..)\n",
    "import matplotlib.pyplot as plt\n",
    "%matplotlib inline\n",
    "\n",
    "# New for today! Tensorflow from Google: \n",
    "# https://www.tensorflow.org\n",
    "import tensorflow as tf\n",
    "\n",
    "# Please see the following page for getting \n",
    "# started guide and tutorials:\n",
    "# https://www.tensorflow.org/get_started/"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "In numpy:\n",
      "a =  [1. 2. 3.]\n",
      "a =  [2. 2. 2.]\n",
      "a + b =  [3. 4. 5.]\n",
      "a ** b =  [1. 4. 9.]\n",
      "c = (a + b) * (a ** b) =  [ 3. 16. 45.]\n"
     ]
    }
   ],
   "source": [
    "# What we are used to in standard programming:\n",
    "\n",
    "a = np.array([1.0, 2.0, 3.0])\n",
    "\n",
    "b = np.array([2.0, 2.0, 2.0])\n",
    "\n",
    "a_plus_b = a + b\n",
    "\n",
    "a_power_b = a ** b\n",
    "\n",
    "c = a_plus_b * a_power_b\n",
    "\n",
    "print('In numpy:')\n",
    "print('a = ', a)\n",
    "print('a = ', b)\n",
    "print('a + b = ', a_plus_b)\n",
    "print('a ** b = ', a_power_b)\n",
    "print('c = (a + b) * (a ** b) = ', c)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "In tensorflow:\n",
      "a =  Tensor(\"a:0\", shape=(3,), dtype=float64)\n",
      "a =  Tensor(\"b:0\", shape=(3,), dtype=float64)\n",
      "a + b =  Tensor(\"add:0\", shape=(3,), dtype=float64)\n",
      "a ** b =  Tensor(\"pow:0\", shape=(3,), dtype=float64)\n",
      "c = (a + b) * (a ** b) =  Tensor(\"mul:0\", shape=(3,), dtype=float64)\n"
     ]
    }
   ],
   "source": [
    "# Let's do this with tensorflow:\n",
    "#\n",
    "# (Using t_ to denote tensorflow)\n",
    "#\n",
    "\n",
    "# Reset tensorflow!\n",
    "tf.reset_default_graph()\n",
    "\n",
    "t_a = tf.constant(a, name='a')\n",
    "\n",
    "t_b = tf.constant(b, name='b')\n",
    "\n",
    "t_a_plus_b = t_a + t_b\n",
    "\n",
    "t_a_power_b = t_a ** t_b\n",
    "\n",
    "t_c = t_a_plus_b * t_a_power_b\n",
    "\n",
    "print('In tensorflow:')\n",
    "print('a = ', t_a)\n",
    "print('a = ', t_b)\n",
    "print('a + b = ', t_a_plus_b)\n",
    "print('a ** b = ', t_a_power_b)\n",
    "print('c = (a + b) * (a ** b) = ', t_c)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Can ignore this code for now, this just creates a \n",
    "# visualisation file so we can see what is going on!\n",
    "\n",
    "with tf.Session() as session:\n",
    "    summary_file_writer = tf.summary.FileWriter(\n",
    "        'visualisation_files', session.graph)\n",
    "    summary_file_writer.flush()\n",
    "\n",
    "# This has created a folder called \"visualisation_files\" \n",
    "# with some data in it that we can view with the command:\n",
    "#\n",
    "# tensorboard --logdir=visualisation_files\n",
    "#"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### This has created a computational graph!\n",
    "\n",
    "<img src=\"graph_01.png\" width=\"75%\">\n",
    "\n",
    "### Why would we want to do this?\n",
    "\n",
    "..."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "result =  [ 3. 16. 45.]\n"
     ]
    }
   ],
   "source": [
    "# First - does it actually work?\n",
    "\n",
    "# In order to use the graph we have to \"run\" it!\n",
    "# In tensorflow, we need to run things inside of a session\n",
    "\n",
    "# This line creates a session..\n",
    "with tf.Session() as session:\n",
    "    # Inside here we can use the \"session\" object created..\n",
    "    \n",
    "    # Let's run our graph to actually compute something!\n",
    "    result = session.run(t_c)\n",
    "    \n",
    "    print(\"result = \", result)\n",
    "    \n",
    "# Outside of the \"with\" statement the session object is \n",
    "# deleted and we can no longer use it\n",
    "#\n",
    "# This line would cause an error:\n",
    "# result_again = session.run(t_c)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "In numpy:\n",
      "a =  [1. 2. 3.]\n",
      "a =  [2. 2. 2.]\n",
      "a + b =  [3. 4. 5.]\n",
      "a ** b =  [1. 4. 9.]\n",
      "c = (a + b) * (a ** b) =  [ 3. 16. 45.]\n",
      "In tensorflow session:\n",
      "a =  [1. 2. 3.]\n",
      "a =  [2. 2. 2.]\n",
      "a + b =  [3. 4. 5.]\n",
      "a ** b =  [1. 4. 9.]\n",
      "c = (a + b) * (a ** b) =  [ 3. 16. 45.]\n"
     ]
    }
   ],
   "source": [
    "# Let's check everything makes sense:\n",
    "\n",
    "print('In numpy:')\n",
    "print('a = ', a)\n",
    "print('a = ', b)\n",
    "print('a + b = ', a_plus_b)\n",
    "print('a ** b = ', a_power_b)\n",
    "print('c = (a + b) * (a ** b) = ', c)\n",
    "\n",
    "with tf.Session() as session:\n",
    "    print('In tensorflow session:')   \n",
    "    print('a = ', session.run(t_a))\n",
    "    print('a = ', session.run(t_b))\n",
    "    print('a + b = ', session.run(t_a_plus_b))\n",
    "    print('a ** b = ', session.run(t_a_power_b))\n",
    "    print('c = (a + b) * (a ** b) = ', session.run(t_c))\n",
    "\n",
    "# Everything should be the same!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### So it seems to work!\n",
    "\n",
    "### But isn't this more effort than the numpy version?\n",
    "\n",
    "..."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### What if we are doing an optimistion?\n",
    "\n",
    "We want to fit a Gaussian distribution $\\mathcal{N}(\\mu, \\sigma^2)$ to a set of numbers $X = \\{ x_0, x_1, \\dots, x_{N-1} \\}$\n",
    "\n",
    "If we assume the numbers are i.i.d. (indentically and independently distributed) samples from a Gaussian then the likelihood of $X$ is given by:\n",
    "\n",
    "\\begin{align}\n",
    "p(X) &= p(x_0) \\cdot p(x_1) \\cdot \\dots \\cdot p(x_{N-1}) \\\\\n",
    " &= \\mathcal{N}(x_0 \\,|\\, \\mu, \\sigma^2) \\cdot \\mathcal{N}(x_1 \\,|\\, \\mu, \\sigma^2) \\cdot \\dots \\cdot \\mathcal{N}(x_{N-1} \\,|\\, \\mu, \\sigma^2) \\\\\n",
    " &= \\prod_{n=0}^{N-1} \\mathcal{N}(x_{n} \\,|\\, \\mu, \\sigma^2) \\\\\n",
    " &= \\prod_{n=0}^{N-1} \\frac{1}{\\sqrt{2\\pi\\sigma^2}} \n",
    "    \\exp{\\left( - \\frac{(x_{n} - \\mu)^2}{2\\sigma^2} \\right)}\n",
    "\\end{align}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Top Tip!** When working with exponential family of distributions it often helps to work in the log domain..\n",
    "\n",
    "\\begin{align}\n",
    "\\log  p(X)  &= \\sum_{n=0}^{N-1} \n",
    "    -\\frac{1}{2} \\log{\\left( 2\\pi\\sigma^2 \\right)}\n",
    "    -\\frac{(x_{n} - \\mu)^2}{2\\sigma^2}\n",
    "\\end{align}\n",
    "\n",
    "So, we have the *maximum likelihood* fit to the parameters when we find the values of $\\mu$ and $\\sigma^2$ that maximise $p(X)$ which (since $\\log\\,(\\cdot)$ is a concave function) occurs at the same time that $\\log p(X)$ is maximised."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In our case we can find an analytic solution for\n",
    "\n",
    "\\begin{align}\n",
    "\\mu^{*} &= {\\arg\\max}_{\\mu} \\, \\log p(X) \\\\\n",
    "{\\sigma^{*}}^2 &= {\\arg\\max}_{\\sigma^2} \\, \\log p(X) \\\\\n",
    "\\end{align}\n",
    "\n",
    "But let's pretend that the problem was more complicated and we needed to use *optimisation* to solve the problem.."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In previous lectures you have seen that you need to be able to calculate gradients of the objective function ($\\log p(X)$) wrt the parameters that you are optimising ($\\mu$ and $\\sigma^2$).\n",
    "\n",
    "*Let's see how to do this in tensorflow..*"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "X = \n",
      " [2.321 5.012 3.248 2.968 3.611 3.659 0.697 1.405\n",
      " 2.598 2.57  1.359 3.093 3.052 2.728 4.011 0.96\n",
      " 4.414 4.034 4.248 3.415]\n"
     ]
    }
   ],
   "source": [
    "# First let's generate some numbers to fit the data to..\n",
    "\n",
    "# How many values of x?\n",
    "N = 20\n",
    "\n",
    "# Pick the real mean and variance..\n",
    "mu_true = 2.5\n",
    "sigma_true = 1.5\n",
    "\n",
    "x_n = np.random.normal(mu_true, sigma_true, N)\n",
    "\n",
    "np.set_printoptions(precision=3, linewidth=50)\n",
    "print('X = \\n', np.transpose(x_n))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We are now going to build our tensorflow graph but we are going to account for the fact that $\\mu$ and $\\sigma^2$ are no longer constants since we wish to vary their values to find the maximum of $\\log p(X)$. As with the optimisation lectures, we need to start with a guess for the values of $\\mu$ and $\\sigma^2$; in this case, we will start with\n",
    "\\begin{align}\n",
    "\\mu_{\\mathrm{initial}} &= 1\\\\\n",
    "\\sigma^2_{\\mathrm{initial}} &= 1\n",
    "\\end{align}\n",
    "\n",
    "**Top Tip!** Care needs to be taken with $\\sigma$ since it can only be a positive value (unlike $\\mu$ which can be any real number). In general, tensorflow variables can be positive or negative. In this example we square the value of `t_sigma` before using it to ensure that `t_sigma_2` is a positive value but we shouldn't, therefore, use the value for `t_sigma` directly in calculations..\n",
    "\n",
    "As a reminder, we want to find:\n",
    "\\begin{align}\n",
    "\\log  p(X)  &= \\sum_{n=0}^{N-1} \n",
    "    -\\frac{1}{2} \\log{\\left( 2\\pi\\sigma^2 \\right)}\n",
    "    -\\frac{(x_{n} - \\mu)^2}{2\\sigma^2} \n",
    "%    \\\\\n",
    "%    &=  -\\frac{N}{2} \\log{\\left( 2\\pi\\sigma^2 \\right)}\n",
    "%    - \\frac{1}{2\\sigma^2} \\sum_{n=0}^{N-1}\\left(x_{n} - \\mu\\right)^2\n",
    "\\end{align}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Reset tensorflow to remove our old a, b, etc..\n",
    "tf.reset_default_graph()\n",
    "\n",
    "# Our initial guesses..\n",
    "mu_initial_guess = 1.0\n",
    "sigma_initial_guess = np.sqrt(1.0)\n",
    "\n",
    "# The data to fit to\n",
    "t_x_n = tf.constant(x_n, name='X')\n",
    "\n",
    "# Note: mu and sigma are now *variables* not constants!\n",
    "# We need to specify their data type and initial value.. \n",
    "t_mu = tf.Variable(mu_initial_guess, \n",
    "                   dtype=tf.float64, \n",
    "                   name=\"mu\")\n",
    "t_sigma = tf.Variable(sigma_initial_guess, \n",
    "                      dtype=tf.float64, \n",
    "                      name=\"sigma\")\n",
    "\n",
    "# Note: this step is important - don't use t_sigma directly!! \n",
    "t_sigma_2 = t_sigma ** 2.0\n",
    "\n",
    "# Calculate log p(X) terms..\n",
    "\n",
    "t_x_minus_mu_2 = (t_x_n - t_mu) ** 2.0\n",
    "t_denom = 2.0 * t_sigma_2\n",
    "t_sigma_term = - 0.5 * tf.log(2.0 * np.pi * t_sigma_2)\n",
    "\n",
    "t_log_P_terms = t_sigma_term - (t_x_minus_mu_2 / t_denom)\n",
    "\n",
    "# The sum is performed by a reduction in tensorflow \n",
    "# (since a vector goes in and a scalar comes out)\n",
    "# but this is effectively the same as np.sum(...)\n",
    "t_log_P = tf.reduce_sum(t_log_P_terms)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Tensorflow log p(X) =  -70.35135977725858\n",
      "(using initial guesses for mu and sigma)\n",
      "\n",
      "Value from scipy stats package =  -70.35135977725858\n",
      "\n",
      "Everything working!\n"
     ]
    }
   ],
   "source": [
    "# Let's just check that we calculated things correctly:\n",
    "\n",
    "with tf.Session() as session:\n",
    "    # IMPORTANT! Need to run this at the start to\n",
    "    # initialise the values for the variables\n",
    "    # t_mu and t_sigma. You will get an error if \n",
    "    # you forget!\n",
    "    session.run(tf.global_variables_initializer())\n",
    "    \n",
    "    test_value = session.run(t_log_P)\n",
    "    print('Tensorflow log p(X) = ', test_value)\n",
    "    print('(using initial guesses for mu and sigma)\\n')\n",
    "\n",
    "# Check with scipy..\n",
    "from scipy.stats import norm\n",
    "check_value = np.sum(norm.logpdf(x_n, \n",
    "                                 mu_initial_guess, \n",
    "                                 sigma_initial_guess))\n",
    "print('Value from scipy stats package = ', check_value)\n",
    "\n",
    "assert(np.isclose(test_value, check_value))\n",
    "\n",
    "print('\\nEverything working!')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "New tensorflow log p(X) =  -61.492242187791284\n",
      "(using new values for mu and sigma)\n",
      "\n",
      "New value from scipy stats package =  -61.49224218779127\n",
      "\n",
      "Everything working!\n"
     ]
    }
   ],
   "source": [
    "# We can even go crazy and check with different \n",
    "# values of the parameters..\n",
    "\n",
    "mu_new_test_value = 3.3\n",
    "sigma_new_test_value = 0.5\n",
    "\n",
    "with tf.Session() as session:\n",
    "    # IMPORTANT! (see above..)\n",
    "    session.run(tf.global_variables_initializer())\n",
    "    \n",
    "    # Change the values of the variables while the\n",
    "    # session is running..\n",
    "    session.run(t_mu.assign(mu_new_test_value))\n",
    "    session.run(t_sigma.assign(sigma_new_test_value))\n",
    "    \n",
    "    test_value = session.run(t_log_P)\n",
    "    print('New tensorflow log p(X) = ', test_value)\n",
    "    print('(using new values for mu and sigma)\\n')\n",
    "\n",
    "# Check with scipy..\n",
    "from scipy.stats import norm\n",
    "check_value = np.sum(norm.logpdf(x_n, \n",
    "                                 mu_new_test_value, \n",
    "                                 sigma_new_test_value))\n",
    "print('New value from scipy stats package = ', check_value)\n",
    "\n",
    "assert(np.isclose(test_value, check_value))\n",
    "\n",
    "print('\\nEverything working!')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "*So finally, we get to the advantage of tensorflow!*\n",
    "\n",
    "Now, we can calcuate the objective function and we can calculate the value of the objective when changing the input parameters. \n",
    "\n",
    "This is great for optimisation (since we are going to need to change the parameters to increase the objective) but what we really need for the optimisation is to calculate the **gradient of the objective wrt to the parameters**. Let's see how to do that in tensorflow.."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Gradient wrt mu =  [39.403689273669976]\n",
      "Gradient wrt sigma =  [83.94517822633026]\n"
     ]
    }
   ],
   "source": [
    "with tf.Session() as session:\n",
    "    # IMPORTANT! (see above..)\n",
    "    session.run(tf.global_variables_initializer())\n",
    "    \n",
    "    t_gradient_wrt_mu = tf.gradients(t_log_P, \n",
    "                                     t_mu)\n",
    "    t_gradient_wrt_sigma = tf.gradients(t_log_P, \n",
    "                                        t_sigma)\n",
    "    \n",
    "    grad_mu = session.run(t_gradient_wrt_mu)\n",
    "    grad_sigma = session.run(t_gradient_wrt_sigma)\n",
    "    \n",
    "    print('Gradient wrt mu = ', grad_mu)\n",
    "    print('Gradient wrt sigma = ', grad_sigma)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Shall we check that result. Remember we have:\n",
    "\n",
    "\\begin{align}\n",
    "\\log  p(X)  &= \\sum_{n=0}^{N-1} \n",
    "    -\\frac{1}{2} \\log{\\left( 2\\pi\\sigma^2 \\right)}\n",
    "    -\\frac{(x_{n} - \\mu)^2}{2\\sigma^2} \\\\\n",
    "    &=  -\\frac{N}{2} \\log{\\left( 2\\pi\\sigma^2 \\right)}\n",
    "    - \\frac{1}{2\\sigma^2} \\sum_{n=0}^{N-1}\\left(x_{n} - \\mu\\right)^2\n",
    "\\end{align}\n",
    "\n",
    "So for $\\mu$ we have:\n",
    "\n",
    "\\begin{align}\n",
    "\\frac{\\partial \\log  p(X)}{\\partial \\mu}   \n",
    "    &= - 0\n",
    "    - \\frac{1}{2\\sigma^2} \\frac{\\partial}{\\partial \\mu}  \\sum_{n=0}^{N-1} \\left(x_{n} - \\mu\\right)^2 \\\\\n",
    "    &= - \\frac{1}{2\\sigma^2} \\sum_{n=0}^{N-1} \\frac{\\partial}{\\partial \\mu} \\left(x_{n} - \\mu\\right)^2 \\\\\n",
    "    &= - \\frac{1}{2\\sigma^2} \\sum_{n=0}^{N-1} 2 \\left(x_{n} - \\mu\\right) \\frac{\\partial}{\\partial \\mu}\\left(x_{n} - \\mu\\right)  \\\\\n",
    "    &= \\frac{1}{\\sigma^2} \\sum_{n=0}^{N-1} \\left(x_{n} - \\mu\\right)\n",
    "\\end{align}\n",
    "\n",
    "where we used the \"chain rule\" a number of times.."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Our analytic gradient wrt mu =  39.403689273669976\n",
      "Tensorflow gradient wrt mu =  [39.403689273669976]\n",
      "\n",
      "Excellent! tensorflow calculated the gradient for us :)\n"
     ]
    }
   ],
   "source": [
    "# numpy check of gradient wrt mu\n",
    "\n",
    "grad_mu_check = np.sum(x_n - mu_initial_guess) / \\\n",
    "                (sigma_initial_guess ** 2)\n",
    "    \n",
    "print('Our analytic gradient wrt mu = ', grad_mu_check)\n",
    "\n",
    "print('Tensorflow gradient wrt mu = ', grad_mu)\n",
    "\n",
    "assert(np.isclose(grad_mu, grad_mu_check))\n",
    "\n",
    "print('\\nExcellent! tensorflow calculated the gradient for us :)')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Everyone should now be in awe!\n",
    "\n",
    "This might seem like something trivial but hopefully you can see that actually quite a lot of maths and then coding went into determining the gradient. \n",
    "\n",
    "In fact, you can do the same to check the value for the gradient wrt $\\sigma^2$.\n",
    "\n",
    "When we calculated the result using the chain rule. Since tensorflow built up a graph of the operations, it is able to apply the chain rule results for us automatically.\n",
    "\n",
    "This:\n",
    "\\begin{align}\n",
    "\\log  p(X)  &= \\sum_{n=0}^{N-1} \n",
    "    -\\frac{1}{2} \\log{\\left( 2\\pi\\sigma^2 \\right)}\n",
    "    -\\frac{(x_{n} - \\mu)^2}{2\\sigma^2} \n",
    "%    \\\\\n",
    "%    &=  -\\frac{N}{2} \\log{\\left( 2\\pi\\sigma^2 \\right)}\n",
    "%    - \\frac{1}{2\\sigma^2} \\sum_{n=0}^{N-1}\\left(x_{n} - \\mu\\right)^2\n",
    "\\end{align}\n",
    "has become this:\n",
    "<img src=\"graph_02.png\" width=\"50%\">\n",
    "\n",
    "For example, the `pow` operation represents $r = a^b$ for the inputs $a,b$ and result $r$. Tensorflow then knows that $\\frac{\\partial r}{\\partial a} = b a^{b-1}$, and by chaining these operations together it can work backwards through the graph (from $\\log p(X)$ at the top to $\\mu$ at the bottom) to calculate the gradient.\n",
    "\n",
    "Therefore, the tensorflow graph has multiple uses. A forward pass can calculate the objective for the current set of parameters and a backwards pass can calcuate the gradients of an objective wrt any of the parameters."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## But the fun doesn't end here!\n",
    "\n",
    "In fact tensorflow has actually done all the work to do the optimisation part, not just calculate the derivatives. So we can now run a full optimisation with our graph and it will use the gradients internally."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "iter   10, log P(X) = -44.641\n",
      "iter   20, log P(X) = -31.123\n",
      "iter   30, log P(X) = -31.122\n",
      "iter   40, log P(X) = -31.122\n",
      "iter   50, log P(X) = -31.122\n",
      "\n",
      "After optimisation:\n",
      "Tensorflow mu =  2.970184463683499\n",
      "Tensorflow sigma =  1.1470100654949587\n",
      "\n",
      "Analytic estimates:\n",
      "Estimated mu =  2.970184463683499\n",
      "Estimated std =  1.1470100655080047\n",
      "\n",
      "Ground truth values:\n",
      "True mu =  2.5\n",
      "True sigma =  1.5\n"
     ]
    }
   ],
   "source": [
    "# Create a gradient descent optimiser that uses a\n",
    "# certain step size (learning_rate)..\n",
    "optimiser = tf.train.GradientDescentOptimizer(learning_rate=0.05)\n",
    "\n",
    "# We want to maximise log p(X) therefore we \n",
    "# need to minimise - log p(X)\n",
    "t_objective = - t_log_P\n",
    "\n",
    "# We want to optimise wrt mu and sigma\n",
    "vars_to_optimise = [t_mu, t_sigma]\n",
    "\n",
    "minimize_operation = optimiser.minimize(t_objective,\n",
    "                                       var_list=vars_to_optimise)\n",
    "\n",
    "# Number of iterations to perform\n",
    "num_iterations = 50\n",
    "\n",
    "with tf.Session() as session:\n",
    "    # IMPORTANT! (see above..)\n",
    "    session.run(tf.global_variables_initializer())\n",
    "    \n",
    "    # Run a number of iterations of gradient descent..\n",
    "    for iteration in range(num_iterations):\n",
    "        # At each iteration evaluate the minimize_operation\n",
    "        # to perform the gradient descent step and also\n",
    "        # keep track of the current value..\n",
    "        step, cost = session.run([minimize_operation, t_log_P])\n",
    "        \n",
    "        # Print out the value of log P every 10 iterations..\n",
    "        if ((iteration + 1) % 10 == 0):\n",
    "            print('iter %4d, log P(X) = %0.3f' % \n",
    "                  (iteration + 1, cost))\n",
    "    \n",
    "    # Get the final results of the optimisation..\n",
    "    mu_optimised = session.run(t_mu)\n",
    "    sigma_optimised = session.run(t_sigma)\n",
    "    \n",
    "    print('\\nAfter optimisation:')\n",
    "    print('Tensorflow mu = ', mu_optimised)\n",
    "    print('Tensorflow sigma = ', sigma_optimised)\n",
    "\n",
    "print('\\nAnalytic estimates:')\n",
    "print('Estimated mu = ', np.mean(x_n))\n",
    "print('Estimated std = ', np.std(x_n))\n",
    "\n",
    "print('\\nGround truth values:')\n",
    "print('True mu = ', mu_true)\n",
    "print('True sigma = ', sigma_true)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Excellent! So we agree with the analytic estimate :)\n",
    "\n",
    "Of course, the values don't match the true estimate since we didn't have a very large sample size.\n",
    "\n",
    "What if we want to run again with more samples?\n",
    "\n",
    "Unfortunately, we made `t_x_n` a constant at the start of our tensorflow code so now we can't change it. Instead, we could have made it a `placeholder`. This tells tensorflow \"there will be some data here but I'm going to give it to you later\".\n",
    "\n",
    "*How do we give the data later on?*\n",
    "\n",
    "We can provide the values to placeholders by specifying a \"feed dictionary\" to `session.run`. This means, \"during this session use the following values to replace all the placeholders\".\n",
    "\n",
    "Let's do our example again.."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Reset tensorflow to remove our old a, b, etc..\n",
    "tf.reset_default_graph()\n",
    "\n",
    "# THIS TIME USE A PLACEHOLDER!\n",
    "#\n",
    "# The data to fit to is provided as a placeholder.\n",
    "# We need to tell it what type of data we will provide..\n",
    "t_x_n = tf.placeholder(dtype=tf.float64, name='X')\n",
    "\n",
    "# EVERYTHING ELSE IS AS IT WAS BEFORE..\n",
    "\n",
    "# Note: mu and sigma are now *variables* not constants!\n",
    "# We need to specify their data type and initial value.. \n",
    "t_mu = tf.Variable(mu_initial_guess, \n",
    "                   dtype=tf.float64, \n",
    "                   name=\"mu\")\n",
    "t_sigma = tf.Variable(sigma_initial_guess, \n",
    "                      dtype=tf.float64, \n",
    "                      name=\"sigma\")\n",
    "\n",
    "# Note: this step is important - don't use t_sigma directly!! \n",
    "t_sigma_2 = t_sigma ** 2.0\n",
    "\n",
    "# Calculate log p(X) terms..\n",
    "\n",
    "t_x_minus_mu_2 = (t_x_n - t_mu) ** 2.0\n",
    "t_denom = 2.0 * t_sigma_2\n",
    "t_sigma_term = - 0.5 * tf.log(2.0 * np.pi * t_sigma_2)\n",
    "\n",
    "t_log_P_terms = t_sigma_term - (t_x_minus_mu_2 / t_denom)\n",
    "\n",
    "# The sum is performed by a reduction in tensorflow \n",
    "# (since a vector goes in and a scalar comes out)\n",
    "# but this is effectively the same as np.sum(...)\n",
    "t_log_P = tf.reduce_sum(t_log_P_terms)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "iter   10, log P(X) = -44.641\n",
      "iter   20, log P(X) = -31.123\n",
      "iter   30, log P(X) = -31.122\n",
      "iter   40, log P(X) = -31.122\n",
      "iter   50, log P(X) = -31.122\n",
      "\n",
      "After optimisation:\n",
      "Tensorflow mu =  2.970184463683499\n",
      "Tensorflow sigma =  1.1470100654949587\n",
      "\n",
      "Analytic estimates:\n",
      "Estimated mu =  2.970184463683499\n",
      "Estimated std =  1.1470100655080047\n",
      "\n",
      "Ground truth values:\n",
      "True mu =  2.5\n",
      "True sigma =  1.5\n"
     ]
    }
   ],
   "source": [
    "# NOW WHEN WE RUN WE NEED TO FILL IN THE PLACEHOLDER..\n",
    "\n",
    "# Create a gradient descent optimiser that uses a\n",
    "# certain step size (learning_rate)..\n",
    "optimiser = tf.train.GradientDescentOptimizer(learning_rate=0.05)\n",
    "\n",
    "# We want to maximise log p(X) therefore we \n",
    "# need to minimise - log p(X)\n",
    "t_objective = - t_log_P\n",
    "\n",
    "# We want to optimise wrt mu and sigma\n",
    "vars_to_optimise = [t_mu, t_sigma]\n",
    "\n",
    "minimize_operation = optimiser.minimize(t_objective,\n",
    "                                       var_list=vars_to_optimise)\n",
    "\n",
    "# Number of iterations to perform\n",
    "num_iterations = 50\n",
    "\n",
    "with tf.Session() as session:\n",
    "    # IMPORTANT! (see above..)\n",
    "    session.run(tf.global_variables_initializer())\n",
    "    \n",
    "    # Run a number of iterations of gradient descent..\n",
    "    for iteration in range(num_iterations):\n",
    "        # At each iteration evaluate the minimize_operation\n",
    "        # to perform the gradient descent step and also\n",
    "        # keep track of the current value..\n",
    "        #\n",
    "        # NEED TO ADD THE FEED DICTIONARY OTHERWISE WE\n",
    "        # DON'T KNOW WHAT VALUE TO USE FOR t_x_n..\n",
    "        #\n",
    "        step, cost = session.run([minimize_operation, t_log_P],\n",
    "                                 feed_dict={ t_x_n : x_n })\n",
    "        \n",
    "        # Print out the value of log P every 10 iterations..\n",
    "        if ((iteration + 1) % 10 == 0):\n",
    "            print('iter %4d, log P(X) = %0.3f' % \n",
    "                  (iteration + 1, cost))\n",
    "    \n",
    "    # Get the final results of the optimisation..\n",
    "    mu_optimised = session.run(t_mu)\n",
    "    sigma_optimised = session.run(t_sigma)\n",
    "    \n",
    "    print('\\nAfter optimisation:')\n",
    "    print('Tensorflow mu = ', mu_optimised)\n",
    "    print('Tensorflow sigma = ', sigma_optimised)\n",
    "\n",
    "print('\\nAnalytic estimates:')\n",
    "print('Estimated mu = ', np.mean(x_n))\n",
    "print('Estimated std = ', np.std(x_n))\n",
    "\n",
    "print('\\nGround truth values:')\n",
    "print('True mu = ', mu_true)\n",
    "print('True sigma = ', sigma_true)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can now even make this a function and call it with lots of different data.."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [],
   "source": [
    "def find_parameters_using_tensorflow(x_input, \n",
    "                                     learning_rate=0.05):\n",
    "    # Create a gradient descent optimiser that uses a\n",
    "    # certain step size (learning_rate)..\n",
    "    optimiser = tf.train.GradientDescentOptimizer(\n",
    "        learning_rate=learning_rate)\n",
    "\n",
    "    # We want to maximise log p(X) therefore we \n",
    "    # need to minimise - log p(X)\n",
    "    t_objective = - t_log_P\n",
    "\n",
    "    # We want to optimise wrt mu and sigma\n",
    "    vars_to_optimise = [t_mu, t_sigma]\n",
    "\n",
    "    minimize_operation = optimiser.minimize(t_objective,\n",
    "                                           var_list=vars_to_optimise)\n",
    "\n",
    "    # Number of iterations to perform\n",
    "    num_iterations = 50\n",
    "\n",
    "    with tf.Session() as session:\n",
    "        # IMPORTANT! (see above..)\n",
    "        session.run(tf.global_variables_initializer())\n",
    "\n",
    "        # Run a number of iterations of gradient descent..\n",
    "        for iteration in range(num_iterations):\n",
    "            # At each iteration evaluate the minimize_operation\n",
    "            # to perform the gradient descent step and also\n",
    "            # keep track of the current value..\n",
    "            #\n",
    "            # PASS THE ARGUMENT TO THE FUNCTION INTO THE FEED\n",
    "            # DICTIONARY..\n",
    "            #\n",
    "            step, cost = session.run([minimize_operation, t_log_P],\n",
    "                                     feed_dict={ t_x_n : x_input })\n",
    "\n",
    "            # Print out the value of log P every 10 iterations..\n",
    "            if ((iteration + 1) % 10 == 0):\n",
    "                print('iter %4d, log P(X) = %0.3f' % \n",
    "                      (iteration + 1, cost))\n",
    "\n",
    "        # Get the final results of the optimisation..\n",
    "        mu_optimised = session.run(t_mu)\n",
    "        sigma_optimised = session.run(t_sigma)\n",
    "\n",
    "    return mu_optimised, sigma_optimised"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "iter   10, log P(X) = -2070.007\n",
      "iter   20, log P(X) = -1813.668\n",
      "iter   30, log P(X) = -1813.668\n",
      "iter   40, log P(X) = -1813.668\n",
      "iter   50, log P(X) = -1813.668\n",
      "Tensorflow estimates:\n",
      "Tensorflow mu =  2.5218790682954064\n",
      "Tensorflow sigma =  1.48398279889398\n",
      "\n",
      "Analytic estimates:\n",
      "Estimated mu =  2.5218790682954064\n",
      "Estimated std =  1.48398279889398\n",
      "\n",
      "Ground truth values:\n",
      "True mu =  2.5\n",
      "True sigma =  1.5\n"
     ]
    }
   ],
   "source": [
    "# Let's try with a larger N\n",
    "\n",
    "N_bigger = 1000\n",
    "\n",
    "x_bigger = np.random.normal(mu_true, sigma_true, N_bigger)\n",
    "\n",
    "new_mu, new_sigma = find_parameters_using_tensorflow(x_bigger, \n",
    "                                                     learning_rate=0.001)\n",
    "\n",
    "print('Tensorflow estimates:')\n",
    "print('Tensorflow mu = ', new_mu)\n",
    "print('Tensorflow sigma = ', new_sigma)\n",
    "\n",
    "print('\\nAnalytic estimates:')\n",
    "print('Estimated mu = ', np.mean(x_bigger))\n",
    "print('Estimated std = ', np.std(x_bigger))\n",
    "\n",
    "print('\\nGround truth values:')\n",
    "print('True mu = ', mu_true)\n",
    "print('True sigma = ', sigma_true)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## All sorts of more advanced topics\n",
    "\n",
    "- Visualise parts of computation (e.g. Tensorboard)\n",
    "- Reusable components (e.g. modules for neural networks / classifiers / etc..)\n",
    "- Run computations on the GPU instead of the CPU (often faster)\n",
    "- Easy to scale; can distribute computations over an entire cluster!"
   ]
  }
 ],
 "metadata": {
  "hide_input": false,
  "kernelspec": {
   "display_name": "Python 3",
   "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.6.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
