{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "# Automatic Differentiation"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This is the concept: https://en.wikipedia.org/wiki/Automatic_differentiation\n",
    "It's nicely done here: http://codon.com/automatic-differentiation-in-ruby, where I was looking for help.\n",
    "The implementation is based on a concept of Dual Numbers: https://en.wikipedia.org/wiki/Dual_number.\n",
    "The basic idea is to (not like in numeric differentiation) compute change and value together (Dual Numbers!). The first thing - implement a dual numbers class!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# dual numbers class\n",
    "import math\n",
    "class Dual_number:\n",
    "    \"\"\"Dual numbers API\"\"\"\n",
    "    def __init__(self, _real, _dual=0):\n",
    "        self.real = _real\n",
    "        self.dual = _dual\n",
    "    \n",
    "    # present a number to the world:\n",
    "    def __str__(self):\n",
    "        def sign_help(x):\n",
    "            return '+' if self.dual >= 0 else '-' \n",
    "        return str(self.real) + ' ' + sign_help(self.dual) + ' ' + str(abs(self.dual)) + 'e'\n",
    "    \n",
    "    def __add__(self, other):\n",
    "        if not isinstance(other, Dual_number):\n",
    "            other = dual_number(other)\n",
    "            return Dual_number(self.real + other.real, self.dual + other.dual)\n",
    "        else:\n",
    "            return Dual_number(self.real + other.real, self.dual + other.dual)\n",
    "    \n",
    "    def __radd__(self, other):\n",
    "        return Dual_number(self.real + other, self.dual)\n",
    "    \n",
    "    def __sub__(self, other):\n",
    "        if not isinstance(other, Dual_number):\n",
    "            other = dual_number(other)\n",
    "            return Dual_number(self.real - other.real, self.dual - other.dual)\n",
    "        else:\n",
    "            return Dual_number(self.real - other.real, self.dual - other.dual)\n",
    "    \n",
    "    def __rsub__(self, other):\n",
    "        return Dual_number(-(self.real - other), -self.dual)\n",
    "    \n",
    "    def __mul__(self, other):\n",
    "        if not isinstance(other, Dual_number):\n",
    "            other = dual_number(other)\n",
    "            return Dual_number(self.real * other.real, self.real * other.dual + \n",
    "                          self.dual * other.real)\n",
    "        else:\n",
    "            return Dual_number(self.real * other.real, self.real * other.dual + \n",
    "                          self.dual * other.real)\n",
    "    \n",
    "    def __rmul__(self, other):\n",
    "        return Dual_number(self.real * other, self.dual * other)\n",
    "    \n",
    "    def __truediv__(self, other):\n",
    "        if not isinstance(other, Dual_number):\n",
    "            other = dual_number(other)\n",
    "            return Dual_number(self.real / other.real ,  (self.dual * other.real - \n",
    "                                                     self.real * other.dual) / (other.real * other.real))\n",
    "        else:\n",
    "            return Dual_number(self.real / other.real ,  (self.dual * other.real - \n",
    "                                                     self.real * other.dual) / (other.real * other.real))\n",
    "    def __rtruediv__(self, other):\n",
    "        return Dual_number((1 / self.real) * other, ((-1 * self.dual) / self.real * self.real) * other )\n",
    "    \n",
    "def dual_number(x):\n",
    "    return Dual_number(x)\n",
    "        \n",
    "        \n",
    "        \n",
    "    \n",
    "\n",
    "# dual number with dual part eq to zero is  a real number (similar to complex), constant value \n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As seen arithmetic operations are constructed in similar way like in complex class, just have in mind, that e is infinitesimal so  e*e is zero; add and substract like normal algebraic expressions, multiplication is based on a matrix representation[of dual number]. Let's check if everything's all right.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1 + 3e\n",
      "0 - 1e\n",
      "2 + 1e\n",
      "0 - 1e\n",
      "9 + 9e\n",
      "10 + 0e\n",
      "1.0 - 1.0e\n",
      "2.09 + 1e\n",
      "-2 + 2e\n",
      "0.5 + 0.5e\n",
      "2.0 - 2.0e\n"
     ]
    }
   ],
   "source": [
    "d3 = Dual_number(1, 1)\n",
    "d4 = Dual_number(1, 2)\n",
    "print(d3 * d4)\n",
    "print(d3 - d4)\n",
    "print(1 + d3)\n",
    "print(1 - d3)\n",
    "print(3 * d3 * 3)\n",
    "print((3 - d4) * (5 * d3))\n",
    "print(d3 / d4)\n",
    "print(d3 + 1.09)\n",
    "print(d4 - 3)\n",
    "print(d3 / 2)\n",
    "print(2 / d3)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "All seem to be fine; see how it works.\n",
    "The function is 2t^2 + 3 and we have:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# distance function:\n",
    "def s(t):\n",
    "    return 2 * (t * t) + 3\n",
    "    \n",
    "def ds(t):\n",
    "    \"\"\"numeric differentiation\"\"\"\n",
    "    dt = 0.001\n",
    "    s0 = s(t)\n",
    "    s1 = s(t + dt)\n",
    "    return (s1 - s0)/dt\n",
    "\n",
    "def ds_aut(x):\n",
    "    \"automatic differnetiation\"\n",
    "    return 2 * x * x + 3"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Time = 0:\n",
      "0.0019999999998354667\n",
      "3 + 0e\n",
      "Time = 4:\n",
      "16.0020000000074\n",
      "35 + 16e\n"
     ]
    }
   ],
   "source": [
    "print(\"Time = 0:\")\n",
    "t0 = Dual_number(0, 1)\n",
    "print(ds(0))\n",
    "print(ds_aut(t0))\n",
    "print(\"Time = 4:\")\n",
    "t0 = Dual_number(4, 1)\n",
    "print(ds(4))\n",
    "print(ds_aut(t0))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# Another function -  tangent \n",
    "def s(t):\n",
    "    return tan(t)\n",
    "def ds(t):\n",
    "    \"\"\"numeric differentiation\"\"\"\n",
    "    dt = 0.001\n",
    "    s0 = s(t)\n",
    "    s1 = s(t + dt)\n",
    "    return (s1 - s0)/dt\n",
    "\n",
    "def ds_aut(x):\n",
    "    \"automatic differnetiation\"\n",
    "    return tan(x)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Time = 0:\n",
      "1.0000003333334668\n",
      "0.0 + 1.0e\n",
      "Time = PI / 3:\n",
      "4.006941562015198\n",
      "1.7320508075688767 + 3.9999999999999982e\n"
     ]
    }
   ],
   "source": [
    "print(\"Time = 0:\")\n",
    "t0 = Dual_number(0, 1)\n",
    "print(ds(0))\n",
    "print(ds_aut(t0))\n",
    "print(\"Time = PI / 3:\")\n",
    "t0 = Dual_number(math.pi / 3, 1)\n",
    "print(ds(math.pi / 3))\n",
    "print(ds_aut(t0))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As seen our dual numbers do the job, we set dual part of time to one: this maybe the change for one second, the real part is the distance travelled or whatever - the point where we calculate the derivative. The cool thing is, that the result is accurate - like in symbolic calculation; but numeric is as acurate as the step (dt) is.    \n",
    "All what is left is define functions, I've done a few here:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "def sin(x):\n",
    "    if isinstance(x, Dual_number):\n",
    "        return Dual_number(math.sin(x.real), x.dual * math.cos(x.real))\n",
    "    else:\n",
    "        return math.sin(x)\n",
    "    \n",
    "def cos(x):\n",
    "    if isinstance(x, Dual_number):\n",
    "        return Dual_number(math.cos(x.real), -x.dual * math.sin(x.real))\n",
    "    else:\n",
    "        return math.sin(x)\n",
    "\n",
    "def tan(x):\n",
    "    \"\"\"tangent x for |x|  < pi / 2 \"\"\"\n",
    "    if isinstance(x, Dual_number):\n",
    "        return sin(x) / cos(x)\n",
    "    else:\n",
    "        return math.tan(x)\n",
    "    \n",
    "def cot(x):\n",
    "    if isinstance(x, Dual_number):\n",
    "        return cos(x) / sin(x)\n",
    "    else:\n",
    "        return math.cos(x) / math.sin(x)\n",
    "\n",
    "def sinh(x):\n",
    "    if isinstance(x, Dual_number):\n",
    "        return  Dual_number(math.sinh(x.real), x.dual * math.cosh(x.real))\n",
    "    else:\n",
    "        return math.sinh(x)\n",
    "\n",
    "def cosh(x):\n",
    "    if isinstance(x, Dual_number):\n",
    "        return Dual_number(math.cosh(x.real), -x.dual * math.sinh(x.real))\n",
    "    else:\n",
    "        return math.cosh(x)\n",
    "\n",
    "def tanh(x):\n",
    "    if isinstance(x, Dual_number):\n",
    "        return sinh(x) / cosh(x)\n",
    "    else:\n",
    "        return math.tanh(x)\n",
    "\n",
    "def coth(x):\n",
    "    \"\"\"coth if x or x.real not eq to 0\"\"\"\n",
    "    if isinstance(x, Dual_number):\n",
    "        return cosh(x) / sinh(x)\n",
    "    else:\n",
    "        return math.cosh(x) / math.sinh(x)\n",
    "        \n",
    "def sqrt(x):\n",
    "    if isinstance(x, Dual_number):\n",
    "        return Dual_number(math.sqrt(x.real), (x.dual) / (2 * sqrt(x.real)))\n",
    "    else:\n",
    "        return math.sqrt(x)\n",
    "def cube_root(x):\n",
    "    if isinstance(x, Dual_number):\n",
    "        return Dual_number(x.real ** (1/3), (x.dual) / (3 * (x.real ** 2)**(1 / 3)))\n",
    "    else:\n",
    "        return x ** (1 / 3)\n",
    "\n",
    "def exp(x):\n",
    "    if isinstance(x, Dual_number):\n",
    "        return Dual_number(math.exp(x.real), math.exp(x.real) * x.dual)\n",
    "    else:\n",
    "        return math.exp(x)\n",
    "    \n",
    "def ln(x):\n",
    "    if isinstance(x, Dual_number):\n",
    "        if not x.dual is 0:\n",
    "            return Dual_number(math.log(x.real) ,  (x.real / x.dual))\n",
    "        else:\n",
    "            return Dual_number(math.log(x.real))\n",
    "    else:\n",
    "        return math.log(x)\n",
    "\n",
    "\n",
    "def pow(a, x):\n",
    "    if isinstance(x, Dual_number):\n",
    "        if not x.dual is 0:\n",
    "            return Dual_number(math.pow(a, x.real), math.pow(a, x.real) * x.dual * ln(x.dual))\n",
    "        else:\n",
    "            return Dual_number(math.pow(a, x.real))\n",
    "    else:\n",
    "        return math.pow(a, x)\n",
    "\n",
    "        "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "To derrive them formulas we can used Maclaurin expansions with all factors with e^2 and higher dropped. That's all."
   ]
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "kernelspec": {
   "display_name": "Python [default]",
   "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.5.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}