10-generators.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Iterators and Generators"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In the section on loops we introduced the `range` function, and said that you should think about it as creating a list of numbers. In Python `2.X` this is exactly what it does. In Python `3.X` this is *not* what it does. Instead it creates the numbers one at a time. The difference in speed and memory usage is enormous for very large lists - examples are given [here](http://justindailey.blogspot.se/2011/09/python-range-vs-xrange.html) and [here](https://asmeurer.github.io/python3-presentation/slides.html#42).\n",
    "\n",
    "We can recreate one of the examples from [Meuer's slides](https://asmeurer.github.io/python3-presentation/slides.html#44) in detail:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def naivesum_list(N):\n",
    "  \"\"\"\n",
    "  Naively sum the first N integers\n",
    "  \"\"\"\n",
    "  A = 0\n",
    "  for i in list(range(N + 1)):\n",
    "      A += i\n",
    "  return A"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We will now see how much memory this uses:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "%load_ext memory_profiler"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "peak memory: 32.38 MiB, increment: 0.50 MiB\n"
     ]
    }
   ],
   "source": [
    "%memit naivesum_list(10**4)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "peak memory: 35.95 MiB, increment: 3.57 MiB\n"
     ]
    }
   ],
   "source": [
    "%memit naivesum_list(10**5)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "peak memory: 70.75 MiB, increment: 34.79 MiB\n"
     ]
    }
   ],
   "source": [
    "%memit naivesum_list(10**6)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "peak memory: 426.25 MiB, increment: 382.75 MiB\n"
     ]
    }
   ],
   "source": [
    "%memit naivesum_list(10**7)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "peak memory: 3856.96 MiB, increment: 3744.59 MiB\n"
     ]
    }
   ],
   "source": [
    "%memit naivesum_list(10**8)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We see that the memory usage is growing very rapidly - as the list gets large it's growing as $N$.\n",
    "\n",
    "Instead we can use the `range` function that yields one integer at a time:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "def naivesum(N):\n",
    "  \"\"\"\n",
    "  Naively sum the first N integers\n",
    "  \"\"\"\n",
    "  A = 0\n",
    "  for i in range(N + 1):\n",
    "      A += i\n",
    "  return A"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "peak memory: 34.35 MiB, increment: 0.13 MiB\n"
     ]
    }
   ],
   "source": [
    "%memit naivesum(10**4)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "peak memory: 34.38 MiB, increment: 0.01 MiB\n"
     ]
    }
   ],
   "source": [
    "%memit naivesum(10**5)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "peak memory: 34.40 MiB, increment: 0.01 MiB\n"
     ]
    }
   ],
   "source": [
    "%memit naivesum(10**6)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "peak memory: 34.41 MiB, increment: 0.00 MiB\n"
     ]
    }
   ],
   "source": [
    "%memit naivesum(10**7)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "peak memory: 34.41 MiB, increment: 0.00 MiB\n"
     ]
    }
   ],
   "source": [
    "%memit naivesum(10**8)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We see that the *memory* usage is unchanged with $N$, making it practical to run much larger calculations. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Iterators"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The `range` function is returning an [*iterator*](https://docs.python.org/3/glossary.html#term-iterator) here. This is an object - a general thing - that represents a stream, or a sequence, of data. The iterator knows how to create the first element of the stream, and it knows how to get the next element. It does not, in general, need to know all of the elements at once.\n",
    "\n",
    "As we've seen above this can save a lot of memory. It can also save time: the code does not need to construct all of the members of the sequence before starting, and it's quite possible you don't need all of them (think about the \"Shortest published mathematical paper\" exercise).\n",
    "\n",
    "An iterator such as `range` is very useful, and there's a lot more useful ways to work with iterators in the `itertools` module. These functions that return iterators, such as `range`, are called [*generators*](https://docs.python.org/3/glossary.html#term-generator), and it's useful to be able to make your own."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Making your own generators"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's look at an example: finding all primes less than $N$ that can be written in the form $4 k - 1$, where $k$ is an integer.\n",
    "\n",
    "We're going to need to calculate all prime numbers less than or equal to $N$. We could write a function that returns all these numbers as a list. However, if $N$ gets large then this will be expensive, both in time and memory. As we only need one number at a time, we can use a generator."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def all_primes(N):\n",
    "    \"\"\"\n",
    "    Return all primes less than or equal to N.\n",
    "    \n",
    "    Parameters\n",
    "    ----------\n",
    "    \n",
    "    N : int\n",
    "        Maximum number\n",
    "        \n",
    "    Returns\n",
    "    -------\n",
    "    \n",
    "    prime : generator\n",
    "        Prime numbers\n",
    "    \"\"\"\n",
    "    \n",
    "    primes = []\n",
    "    for n in range(2, N+1):\n",
    "        is_n_prime = True\n",
    "        for p in primes:\n",
    "            if n%p == 0:\n",
    "                is_n_prime = False\n",
    "                break\n",
    "        if is_n_prime:\n",
    "            primes.append(n)\n",
    "            yield n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This code needs careful examination. First it defines the list of all prime numbers that it currently knows, `primes` (which is initially empty). Then it loops through all integers $n$ from $2$ to $N$ (ignoring $1$ as we know it's not prime).\n",
    "\n",
    "Inside this loop it initially assumes that $n$ is prime. It then checks if any of the known primes exactly divides $n$ (`n%p == 0` checks if $n \\bmod p = 0$). As soon as it finds such a prime divisor it knows that $n$ is not prime it resets the assumption with this new knowledge, then `break`s out of the loop. This statement stops the `for p in primes` loop early, as we don't need to look at later primes. \n",
    "\n",
    "If no known prime ever divides $n$ then at the end of the `for p in primes` loop we will still have `is_n_prime` being `True`. In this case we must have $n$ being prime, so we add it to the list of known primes and return it.\n",
    "\n",
    "It is precisely this point which makes the code above define a generator. We return the value of the prime number found \n",
    "\n",
    "1. using the `yield` keyword, not the `return` keyword, and\n",
    "2. we return the value as soon as it is known.\n",
    "\n",
    "It is the use of the `yield` keyword that makes this function a generator.\n",
    "\n",
    "This means that only the latest prime number is stored for return."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To use the iterator within a loop, we code it in the same way as with the `range` function:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "All prime numbers less than or equal to 20:\n",
      "2\n",
      "3\n",
      "5\n",
      "7\n",
      "11\n",
      "13\n",
      "17\n",
      "19\n"
     ]
    }
   ],
   "source": [
    "print(\"All prime numbers less than or equal to 20:\")\n",
    "for p in all_primes(20):\n",
    "    print(p)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To see what the generator is actually doing, we can step through it one call at a time using the built in `next` function:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "a = all_primes(10)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "2"
      ]
     },
     "execution_count": 17,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "next(a)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "3"
      ]
     },
     "execution_count": 18,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "next(a)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "5"
      ]
     },
     "execution_count": 19,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "next(a)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "7"
      ]
     },
     "execution_count": 20,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "next(a)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "ename": "StopIteration",
     "evalue": "",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mStopIteration\u001b[0m                             Traceback (most recent call last)",
      "\u001b[0;32m<ipython-input-21-3f6e2eea332d>\u001b[0m in \u001b[0;36m<module>\u001b[0;34m()\u001b[0m\n\u001b[0;32m----> 1\u001b[0;31m \u001b[0mnext\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0ma\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m",
      "\u001b[0;31mStopIteration\u001b[0m: "
     ]
    }
   ],
   "source": [
    "next(a)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "So, when the generator gets to the end of its iteration it raises an exception. As seen in previous sections, we could surround the `next` call with a `try` block to capture the `StopIteration` so that we can continue after it finishes. This is effectively what the `for` loop is doing."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can now find all primes (less than or equal to 100, for example) that have the form $4 k - 1$ using"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The prime 3 is 4 * 1 - 1.\n",
      "The prime 7 is 4 * 2 - 1.\n",
      "The prime 11 is 4 * 3 - 1.\n",
      "The prime 19 is 4 * 5 - 1.\n",
      "The prime 23 is 4 * 6 - 1.\n",
      "The prime 31 is 4 * 8 - 1.\n",
      "The prime 43 is 4 * 11 - 1.\n",
      "The prime 47 is 4 * 12 - 1.\n",
      "The prime 59 is 4 * 15 - 1.\n",
      "The prime 67 is 4 * 17 - 1.\n",
      "The prime 71 is 4 * 18 - 1.\n",
      "The prime 79 is 4 * 20 - 1.\n",
      "The prime 83 is 4 * 21 - 1.\n"
     ]
    }
   ],
   "source": [
    "for p in all_primes(100):\n",
    "    if (1+p)%4 == 0:\n",
    "        print(\"The prime {} is 4 * {} - 1.\".format(p, int((1+p)/4)))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Exercise : twin primes"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "A *twin prime* is a pair $(p_1, p_2)$ such that both $p_1$ and $p_2$ are prime and $p_2 = p_1 + 2$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Exercise 1\n",
    "\n",
    "Write a generator that returns twin primes. You can use the generators above, and may want to look at the [itertools](https://docs.python.org/3/library/itertools.html) module together with [its recipes](https://docs.python.org/3/library/itertools.html#itertools-recipes), particularly the `pairwise` recipe."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Exercise 2\n",
    "\n",
    "Find how many twin primes there are with $p_2 < 1000$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Exercise 3\n",
    "\n",
    "Let $\\pi_N$ be the number of twin primes such that $p_2 < N$. Plot how $\\pi_N / N$ varies with $N$ for $N=2^k$ and $k = 4, 5, \\dots 16$. (You should use a logarithmic scale where appropriate!)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Exercise : a basis for the polynomials"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In the section on classes we defined a `Monomial` class to represent a polynomial with leading coefficient $1$. As the $N+1$ monomials $1, x, x^2, \\dots, x^N$ form a basis for the vector space of polynomials of order $N$, $\\mathbb{P}^N$, we can use the `Monomial` class to return this basis.\n",
    "\n",
    "## Exercise 1\n",
    "\n",
    "Define a generator that will iterate through this basis of $\\mathbb{P}^N$ and test it on $\\mathbb{P}^3$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Exercise 2\n",
    "\n",
    "An alternative basis is given by the monomials\n",
    "\n",
    "$$ \\begin{aligned}  p_0(x) &= 1, \\\\ p_1(x) &= 1-x, \\\\ p_2(x) &= (1-x)(2-x), \\\\ \\dots & \\quad \\dots, \\\\ p_N(x) &= \\prod_{n=1}^N (n-x). \\end{aligned} $$\n",
    "\n",
    "Define a generator that will iterate through this basis of $\\mathbb{P}^N$ and test it on $\\mathbb{P}^4$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Exercise 3\n",
    "\n",
    "Use these generators to write another generator that produces a basis of $\\mathbb{P}^3 \\times \\mathbb{P}^4$."
   ]
  }
 ],
 "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"
  },
  "nbconvert": {
   "title": "Iterators and Generators"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}