{
 "metadata": {
  "name": ""
 },
 "nbformat": 3,
 "nbformat_minor": 0,
 "worksheets": [
  {
   "cells": [
    {
     "cell_type": "heading",
     "level": 2,
     "metadata": {},
     "source": [
      "Numerical differentiation : Classes"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# Imports\n",
      "import math\n",
      "import cmath"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 1
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "### Classes \n",
      "\n",
      "Classes are a tool to make the programs modular and hence reusable. Also they lend us a bit of abstraction which makes coding a bit easier. To give a simple example, one can define a class `Quadratic` which codes quadratic polynomials. Calling something like `q = Quadratic(a, b, c)` would define a function `q` such that `q(x) = ax^x + bx + c`. Classes can be used for more complicated things like coding and drawing geometric objects and so on.\n",
      "\n",
      "Classes have their own set of variables, which we call *attribute*s and functions called *method*s. We have already encountered *method*s in our course, figure out which.\n",
      "\n",
      "**All the material in this lecture are in the 7th chapter of Langtangen's *A Primer on Scientific Programming with Python***.\n",
      "\n",
      "First let us code the `Quadratic` class to gain some experience. Note that, though it is not technically necessary, all the classes are conventionally named using a capitalized word."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "class Quadratic :\n",
      "    \"\"\"This class creates quadratic equations.\n",
      "    Attributes :\n",
      "        a, b, c : Coefficients of ax^2 + bx + c\n",
      "    Methods :\n",
      "        valueat(x) : Value of the quadrating at x.\n",
      "        roots() : Returns a tuple containing the two roots\n",
      "            of ax^2 + bx + c = 0\n",
      "        display() : Print the quadratic equation.\n",
      "    \"\"\"\n",
      "    # The doc string is similar to that of functions.\n",
      "    # In a class, we need a function to initialize the class.\n",
      "    def __init__(self, a, b, c) :\n",
      "        \"\"\"Constructor function for the class Quadratic.\"\"\"\n",
      "        # I'll explain the variable self below.\n",
      "        self.a = a\n",
      "        self.b = b\n",
      "        self.c = c\n",
      "        \n",
      "    def valueat(self, x) :\n",
      "        \"\"\"This computes the value of the quadratic at x.\"\"\"\n",
      "        # In the following line note how we refer to the elements\n",
      "        # a, b,c of the quadratic as self.a, self.b etc. self\n",
      "        # is the instance of the current class.\n",
      "        return (self.a)*x*x + (self.b)*x + (self.c)\n",
      "    \n",
      "    def roots(self) :\n",
      "        \"\"\"Returns the roots of the quadratic.\"\"\"\n",
      "        a = self.a\n",
      "        b = self.b\n",
      "        c = self.c\n",
      "        if a !=  0 :\n",
      "            disc = b*b - 4*a*c\n",
      "            if disc >= 0 :\n",
      "                r1 = float(-b + math.sqrt(disc))/2*a\n",
      "                r2 = float(-b - math.sqrt(disc))/2*a\n",
      "            else :\n",
      "                r1 = (-b + cmath.sqrt(disc))/2*a\n",
      "                r2 = (-b - cmath.sqrt(disc))/2*a\n",
      "            retval = (r1, r2)\n",
      "        else :\n",
      "            if b != 0 :\n",
      "                retval = - float(c)/b\n",
      "            else :\n",
      "                print \"No solutions.\"\n",
      "                retval = None\n",
      "        return retval\n",
      "    \n",
      "    def display(self) :\n",
      "        \"\"\"Returns a string printing the quadratic.\"\"\"\n",
      "        str = \"%g x^2 + %g x + %g\" % (self.a, self.b, self.c)\n",
      "        return str"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 2
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Now we use the class to do some computations. For fun we define another function to do the computation."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "def solve_quads(a, b, c) :\n",
      "    q = Quadratic(a, b, c)\n",
      "    # We call q an instance of the class Quadratic.\n",
      "    print \"q = \", q.display(), \"has\", q.roots(), \"as roots. q(\",\n",
      "    print math.pi,\") = \", q.valueat(math.pi)\n",
      "    return None\n",
      "\n",
      "solve_quads(1, 2, 1)\n",
      "solve_quads(1, 0, 1)\n",
      "solve_quads(1, 1, 0)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "q =  1 x^2 + 2 x + 1 has (-1.0, -1.0) as roots. q( 3.14159265359 ) =  17.1527897083\n",
        "q =  1 x^2 + 0 x + 1 has (1j, -1j) as roots. q( 3.14159265359 ) =  10.8696044011\n",
        "q =  1 x^2 + 1 x + 0 has (0.0, -1.0) as roots. q( 3.14159265359 ) =  13.0111970547\n"
       ]
      }
     ],
     "prompt_number": 3
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "### Making classes callable.\n",
      "We eventually want to diffentiate and integrate functions. We'll implement the as classes. What we want to do is to make the instances behave like functions. We do that using the `__call__` method. As an example, we just copy the `Quadratic` class to `Quadratic2` renaming the `valueat` method to be `__call__`."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "class Quadratic2 :\n",
      "    \"\"\"This class creates quadratic equations.\n",
      "    Attributes :\n",
      "        a, b, c : Coefficients of ax^2 + bx + c\n",
      "    Methods :\n",
      "        __call__(x) : Value of the quadrating at x.\n",
      "        roots() : Returns a tuple containing the two roots\n",
      "            of ax^2 + bx + c = 0\n",
      "        display() : Print the quadratic equation.\n",
      "    \"\"\"\n",
      "    # The doc string is similar to that of functions.\n",
      "    # In a class, we need a function to initialize the class.\n",
      "    def __init__(self, a, b, c) :\n",
      "        \"\"\"Constructor function for the class Quadratic.\"\"\"\n",
      "        # I'll explain the variable self below.\n",
      "        self.a = a\n",
      "        self.b = b\n",
      "        self.c = c\n",
      "        \n",
      "    def __call__(self, x) :\n",
      "        \"\"\"This computes the value of the quadratic at x.\"\"\"\n",
      "        # In the following line note how we refer to the elements\n",
      "        # a, b,c of the quadratic as self.a, self.b etc. self\n",
      "        # is the instance of the current class.\n",
      "        return (self.a)*x*x + (self.b)*x + (self.c)\n",
      "    \n",
      "    def roots(self) :\n",
      "        \"\"\"Returns the roots of the quadratic.\"\"\"\n",
      "        a = self.a\n",
      "        b = self.b\n",
      "        c = self.c\n",
      "        if a !=  0 :\n",
      "            disc = b*b - 4*a*c\n",
      "            if disc >= 0 :\n",
      "                r1 = float(-b + math.sqrt(disc))/2*a\n",
      "                r2 = float(-b - math.sqrt(disc))/2*a\n",
      "            else :\n",
      "                r1 = (-b + cmath.sqrt(disc))/2*a\n",
      "                r2 = (-b - cmath.sqrt(disc))/2*a\n",
      "            retval = (r1, r2)\n",
      "        else :\n",
      "            if b != 0 :\n",
      "                retval = - float(c)/b\n",
      "            else :\n",
      "                print \"No solutions.\"\n",
      "                retval = None\n",
      "        return retval\n",
      "    \n",
      "    def display(self) :\n",
      "        \"\"\"Returns a string printing the quadratic.\"\"\"\n",
      "        str = \"%g x^2 + %g x + %g\" % (self.a, self.b, self.c)\n",
      "        return str"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 4
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "q = Quadratic2(1, -3 , 2)\n",
      "print q(10)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "72\n"
       ]
      }
     ],
     "prompt_number": 5
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "### Numeric Differentiation\n",
      "\n",
      "We know that the mathematical definition of the derivative $f'$ of a function $f$ at $x$ is\n",
      "\n",
      "\\begin{equation}\n",
      "  \\lim_{h \\to 0} \\frac{f(x + h) - f(x)}{h}\n",
      "\\end{equation}\n",
      "\n",
      "Instead of finding limit, we can code the derivative as just the fraction $(f(x+h) - f(x))/h$ for small $h$. *Assuming* that as we reduce $h$, the value of the function will be closer and closer to $f'$, we can also write a check if reducing $h$, say by half, doesn't change the value of the fraction beyond some acceptible error. Then we can use that value as an approximate value for the derivative. "
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "class SimpleDerivative :\n",
      "    \"\"\"Implements a callable class for a naive numerical\n",
      "    derivative.\"\"\"\n",
      "    \n",
      "    def __init__(self, f, h = 1.0E-3) :\n",
      "        self.f = f\n",
      "        self.h = float(h)\n",
      "        \n",
      "    def __call__(self, x) :\n",
      "        \"\"\"Return the value of the derivative at x.\"\"\"\n",
      "        f = self.f\n",
      "        h = self.h\n",
      "        der = (f(x + h) - f(x)) / h\n",
      "        return der\n",
      "    \n",
      "    def set_deviation(self, val) :\n",
      "        \"\"\"Set the deviation.\"\"\"\n",
      "        self.h = val"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 6
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "def f(x) : \n",
      "    return x*x\n",
      "\n",
      "df = SimpleDerivative(f,.01)\n",
      "print df(2)\n",
      "df.set_deviation(df.h/100)\n",
      "print df(2)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "4.01\n",
        "4.00010000001\n"
       ]
      }
     ],
     "prompt_number": 7
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Now to end, we modify the class so that when we call, it automatically keeps on halving the deviation till two consecutive derivatives differ by a very small amount. Again as in Newton's method we have to keep track of how many iterations we are doing."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "class NaiveDerivative :\n",
      "    \"\"\"This uses a loop to compute derivatives with smaller\n",
      "    smaller deviations till the difference between successive\n",
      "    computations is less than a predetermined error\n",
      "    \n",
      "    Attributes :\n",
      "        f : function whose derivative we seek,\n",
      "        h : deviaion\n",
      "        err : error\n",
      "        N : max number of interations allowed\n",
      "        \n",
      "    Methods :\n",
      "        single_der(x) : returns (f(x+h) - f(x)) / h\n",
      "        set_dev(d) : sets the deviation to d\n",
      "        set_err(e) : sets the error allowed to e\n",
      "        allow_iter(M) : allow M iterations\n",
      "        __init__ : Constructor\n",
      "        __call__ : the actual code\n",
      "    \"\"\"\n",
      "    \n",
      "    def __init__(self, f, h=0.01, err=1E-5, N=10000) :\n",
      "        self.f = f\n",
      "        self.h = float(h)\n",
      "        self.err = err\n",
      "        self.N = N\n",
      "        \n",
      "    def set_dev(self, d) :\n",
      "        self.h = d\n",
      "        return None\n",
      "    \n",
      "    def set_err(self, e) :\n",
      "        self.err = e\n",
      "        return None\n",
      "    \n",
      "    def allow_iter(self, M) :\n",
      "        self.N = M\n",
      "        return None\n",
      "    \n",
      "    def single_der(self, x) :\n",
      "        f = self.f\n",
      "        h = self.h\n",
      "        return (f(x + h) - f(x))/h\n",
      "    \n",
      "    def __call__(self, x) :\n",
      "        N = self.N\n",
      "        err = self.err\n",
      "        \n",
      "        dfold = self.single_der(x)\n",
      "        self.set_dev(self.h / 2.0)\n",
      "        dfnew = self.single_der(x)\n",
      "        \n",
      "        no_iter = 0\n",
      "        while no_iter < N and abs(dfold - dfnew) >= err :\n",
      "            dfold = dfnew\n",
      "            self.set_dev(self.h / 2.0)\n",
      "            dfnew = self.single_der(x)\n",
      "            no_iter += 1\n",
      "            \n",
      "        if abs(dfold - dfnew) >= err :\n",
      "            print \"Could not converge after %d iterations.\" % no_iter\n",
      "            retval = None\n",
      "        else :\n",
      "            retval = dfnew\n",
      "        return retval\n"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 8
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "def myfn(x) : \n",
      "    return x*x*x + x\n",
      "\n",
      "dmyfn = NaiveDerivative(myfn)\n",
      "dmyfn.set_dev(1)\n",
      "dmyfn.set_err(.1)\n",
      "print dmyfn(1)\n",
      "dmyfn.set_err(1E-5)\n",
      "print dmyfn(1)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "4.0947265625\n",
        "4.00000572205\n"
       ]
      }
     ],
     "prompt_number": 9
    }
   ],
   "metadata": {}
  }
 ]
}