{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "##
Python Generators in Numerical Study of Discrete Dynamical Systems
" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "In this IPython Notebook we present generators as an excelent tool provided by Python for computational mathematics in general, and for the numerical study of \n", "discrete dynamical systems in particular. " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Generators functions form a special class of functions. In order to understand objects that are returned by these functions we recall a few details about iterable objects and iterators in Python.\n", "\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "A container in Python is an object that contains other objects (example of containers are lists, tuples, strings, files, etc). \n", "\n", "A container, `container`, is iterable if it is allowed to run over its items \n", "with a for statement:\n", "
\n",
    "for item in  container:\n",
    "    print(item)\n",
    "
\n", "The lists, strings and tuples are iterable objects. Dictionaries and files are also iterable." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Iteration is the process of accessing and returning one item at a time from a container." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Examples:" ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "2\n", "5\n", "1\n" ] } ], "source": [ "L = [2, 5, 1]#list\n", "for k in L:\n", " print(k)" ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "s\n", "t\n", "r\n", "i\n", "n\n", "g\n" ] } ], "source": [ "S = \"string\"\n", "for c in S:\n", " print(c)" ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "4\n", "9\n", "7\n" ] } ], "source": [ "T=(4, 9, 7)#tuple\n", "for t in T:\n", " print(t)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The iteration over a dictionary `d` returns its keys:" ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "lst\n", "today\n", "nr\n" ] } ], "source": [ "d={'lst': [4,7,1], 'today': 'monday', 'nr': 10}\n", "\n", "for item in d:\n", " print(item)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "An iteration over a file object returns the file's lines:" ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "this is the first line\n", "\n", "second line\n", "\n", "final line\n", "\n" ] } ], "source": [ "f = open(\"toyfile.txt\")# f is a file object i.e. a container for the file's lines \n", "for line in f:\n", " print(line)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The above examples are built-in iterable objects in Python. \n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Iterators" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ " Intuitively, an iterator is an object that produces a stream of items, according to a prescribed protocol. \n", "\n", "From the programming language point of view an iterator is an object of a class that has the methods `__iter__`, and `__next__`. By the iterator protocol these methods act as follows:\n", "\n", "The `__iter__` method returns the iterator object itself, and the `__next__` method returns the next item of the stream.\n", "\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "How does a container, that is an iterable object, relate with an iterator?\n", "\n", "The built-in function `iter` takes a container and returns an iterator." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "In all the above examples, `for item in container` launches the iterator protocol.\n", " The `for` statement calls `iter()` on the container object. \n", "\n", " - `it=iter(container)` is an iterator object and by calling `next(it)`\n", " are accessed the container's items, one at a time. \n", "\n", " - When there are no more items, `next(it)`\n", " raises a `StopIteration` exception, which forces the `for` loop to terminate.\n", " " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Let us explicitly call `iter()` on the list, L, and file object f, defined above:\n" ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "2 5 1\n" ] } ], "source": [ "it_list=iter(L)\n", "print (next(it_list), next(it_list), next(it_list))" ] }, { "cell_type": "code", "execution_count": 7, "metadata": {}, "outputs": [ { "ename": "StopIteration", "evalue": "", "output_type": "error", "traceback": [ "\u001b[1;31m---------------------------------------------------------------------------\u001b[0m", "\u001b[1;31mStopIteration\u001b[0m Traceback (most recent call last)", "\u001b[1;32m\u001b[0m in \u001b[0;36m\u001b[1;34m\u001b[0m\n\u001b[0;32m 1\u001b[0m \u001b[0mit_file\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0miter\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mf\u001b[0m\u001b[1;33m)\u001b[0m \u001b[1;31m# iterator associated to the file object f\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m----> 2\u001b[1;33m \u001b[0mprint\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mnext\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mit_file\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m", "\u001b[1;31mStopIteration\u001b[0m: " ] } ], "source": [ "it_file = iter(f) # iterator associated to the file object f\n", "print(next(it_file))\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We notice that calling three times the next() method for the iterator associated to the list, L,\n", " all three items are successively returned, while a single call `next(it_file)` for the iterator associated to the file object, raises the `StopIteration`\n", "exception. This means that the stream of lines was exhausted when we called `for line in f`, and no \n", " line can be accessed again.\n", "\n", "\n", "If we reopen the file its lines are also accessible:\n" ] }, { "cell_type": "code", "execution_count": 8, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "this is the first line\n", "\n", "second line\n", "\n", "final line\n", "\n" ] } ], "source": [ "it_f = iter(open('toyfile.txt'))\n", "for _ in range(3):\n", " print(next(it_f))\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Hence the iterators associated to built-in containers can exhibit different behaviours. While the lists can be iterated as many times as we want, a file can be iterated only once after it was open." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The reader can experiment the behavior of the iterators `iter(S)`, `iter(T)`, and `iter(d)`, associated to the string S,\n", "tuple, T, and dictionary, d, defined above." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We conclude that:" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Iterables are objects that define an iterator when they are passed to the `iter()` built-in function." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Generators" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Generators functions provide a powerful tool for creating iterators.\n", "\n", "A generator function is defined as any function in Python using `def`.\n", "\n", "A function is a generator function if it contains one or more yield statements within its body." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The generator function below will generate multiple of 3 of positive integer numbers, less than $n$:" ] }, { "cell_type": "code", "execution_count": 9, "metadata": {}, "outputs": [], "source": [ "def gen_mult3(n):\n", " k = 0\n", " while k < n:\n", " yield 3*k\n", " k += 1" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The first call of the generator function defines a generator, i.e. a particular iterator object:" ] }, { "cell_type": "code", "execution_count": 10, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "\n" ] } ], "source": [ "G = gen_mult3(5)\n", "print(G)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "After this call the numbers 0, 3, 6, 9, 12 are not yet generated as we could expect." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Iterating the generator G (which is an iterator), by calling successively `next(G)`,\n", " these numbers are produced one at atime, and at the sixth call a `StopIteration` exception is raised, because the generator exhausted:" ] }, { "cell_type": "code", "execution_count": 11, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "0\n", "3\n", "6\n", "9\n", "12\n" ] }, { "ename": "StopIteration", "evalue": "", "output_type": "error", "traceback": [ "\u001b[1;31m---------------------------------------------------------------------------\u001b[0m", "\u001b[1;31mStopIteration\u001b[0m Traceback (most recent call last)", "\u001b[1;32m\u001b[0m in \u001b[0;36m\u001b[1;34m\u001b[0m\n\u001b[0;32m 1\u001b[0m \u001b[1;32mfor\u001b[0m \u001b[0m_\u001b[0m \u001b[1;32min\u001b[0m \u001b[0mrange\u001b[0m\u001b[1;33m(\u001b[0m\u001b[1;36m6\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m----> 2\u001b[1;33m \u001b[0mprint\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mnext\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mG\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m", "\u001b[1;31mStopIteration\u001b[0m: " ] } ], "source": [ "for _ in range(6):\n", " print(next(G))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "A generator makes a lazzy evaluation, meaning that it yields a value (an object) only on demand\n", "(when it is needed)." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "In order to understand the difference between a Python generator function and a regular function, we recall how behaves\n", "a regular function that contains at least a return statement. \n", "\n", "At each call of a regular function a private namespace is associated. \n", "When a return line is reached in execution, the local variables are destroyed\n", "and the computed object (value) is returned.\n", "\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "In the case of a generator, after each yielded value (with `yield`) the function freezes, i.e. \n", " stops its execution. When it is called again, it resumes its execution at the point where it stopped (it \"remembers\" the last yielded value, and the last executed statement)." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "To illustrate this behaviour, we insert a few `print` lines in the definition of the generator function above:" ] }, { "cell_type": "code", "execution_count": 12, "metadata": {}, "outputs": [], "source": [ "def Gen_mult3(n):\n", " print('this is the first line of the generator function')\n", " k = 0\n", " while k < n:\n", " print('line before yield')\n", " yield 3*k\n", " print ('line after yield')\n", " k+=1\n", " print( 'this is the last line of the generator function') " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Define the generator:" ] }, { "cell_type": "code", "execution_count": 13, "metadata": {}, "outputs": [], "source": [ "gen = Gen_mult3(3)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Running the above cell no string was printed. Hence at the first call the generator is only created and no line of the function's body is executed." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now we iterate the generator `gen`:" ] }, { "cell_type": "code", "execution_count": 14, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "this is the first line of the generator function\n", "line before yield\n", "0\n" ] } ], "source": [ "print(next(gen))" ] }, { "cell_type": "code", "execution_count": 15, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "line after yield\n", "line before yield\n", "3\n" ] } ], "source": [ "print(next(gen))" ] }, { "cell_type": "code", "execution_count": 16, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "line after yield\n", "line before yield\n", "6\n" ] } ], "source": [ "print(next(gen))" ] }, { "cell_type": "code", "execution_count": 17, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "line after yield\n", "this is the last line of the generator function\n" ] }, { "ename": "StopIteration", "evalue": "", "output_type": "error", "traceback": [ "\u001b[1;31m---------------------------------------------------------------------------\u001b[0m", "\u001b[1;31mStopIteration\u001b[0m Traceback (most recent call last)", "\u001b[1;32m\u001b[0m in \u001b[0;36m\u001b[1;34m\u001b[0m\n\u001b[1;32m----> 1\u001b[1;33m \u001b[0mprint\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mnext\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mgen\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m", "\u001b[1;31mStopIteration\u001b[0m: " ] } ], "source": [ "print(next(gen))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Note that Python generators do not give access to previously generated items,\n", "unless the user explicitly stored them in some container (for example in a list).\n", "\n", "If we want to iterate again over data produced by a generator we have to\n", "re-call the generator function, i.e to re-define the generator." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Here is an example of storing the generated values in a list:" ] }, { "cell_type": "code", "execution_count": 18, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "[0, 3, 6, 9, 12, 15]\n" ] } ], "source": [ "ge = gen_mult3(6)\n", "L = [next(ge) for _ in range(6)]\n", "print(L)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "A natural question is whether the builtin `range` function is a generator or not:" ] }, { "cell_type": "code", "execution_count": 19, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "\n" ] } ], "source": [ "R = range(-4, 6, 2)\n", "print(type(R))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "R can be converted to a list:" ] }, { "cell_type": "code", "execution_count": 20, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "[-4, -2, 0, 2, 4]" ] }, "execution_count": 20, "metadata": {}, "output_type": "execute_result" } ], "source": [ "list(R)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "or to an iterator:" ] }, { "cell_type": "code", "execution_count": 21, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "-4 -2\n" ] } ], "source": [ "it_r = iter(R)\n", "print(next(it_r), next(it_r))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Hence `range` is not a generator. It is a class of immutable iterable objects." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Python generators in Computational Mathematics" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Generators allow to deal with infinite streams of \"mathematical objects\"." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "For example, a sequence of real or complex numbers is an infinite mathematical object.\n", "\n", "Let $a_n=\\left(1+\\displaystyle\\frac{1}{n}\\right)^n$, $n\\geq 1$, be the sequence that is known as being convergent to the number e.\n", "\n", "Its terms can be generated by the folowing generator function:" ] }, { "cell_type": "code", "execution_count": 22, "metadata": {}, "outputs": [], "source": [ "def genseq_e():\n", " n = 1\n", " while True:\n", " yield (1.0+1.0/n)**n\n", " n += 1 " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Suppose we want to see how far is $a_{1000}$ from the sequence's limit, $e$. To generate $a_{1000}$ we simply can\n", "iterate the generator `G=genseq_e()`:" ] }, { "cell_type": "code", "execution_count": 23, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "2.7169239322355936\n" ] } ], "source": [ "G = genseq_e()\n", "for _ in range(1, 1000):\n", " next(G)\n", " # the last generated value is a999\n", "a1000 = next(G)\n", "print(a1000)" ] }, { "cell_type": "code", "execution_count": 24, "metadata": {}, "outputs": [], "source": [ "import math" ] }, { "cell_type": "code", "execution_count": 25, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "2.7169239322355936 2.718281828459045 0.0013578962234515046\n" ] } ], "source": [ "print (a1000, math.exp(1), math.fabs(a1000-math.exp(1.0)) )" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "A more elegant solution is to slice the infinite generator, using `itertools.islice`\n", " [https://docs.python.org/3/library/itertools.html](https://docs.python.org/3/library/itertools.html).\n", "\n", "`itertools` is a Python module that contains functions for creating iterators. In particular it allows \n", " to work effectively with generators.\n", "\n", "\n", "When a generator `G` is iterated, the values returned by `next(G)` are counted like the items in any Python sequence.\n", "\n", "A slice:\n", "\n", " `itertools.islice(generator, start, stop, step)`, \n", " \n", " (where `start`, `stop`, and `step` have the same significance as for the \n", "`range(start, stop, step)`)\n", "is an iterator that produces the generator values (objects) counted from index equal to `start`.\n", "\n", "The iterator `itertools.islice(generator, n)` is supposed to have `start=0`, `stop=n`, and `step=1`. " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Hence to get the term $a_{1000}$ of the sequence $a_n$, we proceed as follows:" ] }, { "cell_type": "code", "execution_count": 26, "metadata": {}, "outputs": [], "source": [ "from itertools import islice" ] }, { "cell_type": "code", "execution_count": 27, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "[2.716925287534764]" ] }, "execution_count": 27, "metadata": {}, "output_type": "execute_result" } ], "source": [ "gg = genseq_e()\n", "list(islice(gg, 1000, 1001)) #start=1000, stop=1001" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Hence the slice contains just an element: a1000. To get it we could proceed as follows:" ] }, { "cell_type": "code", "execution_count": 28, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "2.716925287534764 2.718281828459045 0.0013565409242812798\n" ] } ], "source": [ "gg = genseq_e()\n", "a1000 = next(islice(gg, 1000, 1001)) # i.e. the next element of the iterator is its unique element a1000\n", "print(a1000, math.exp(1), math.fabs(a1000-math.exp(1.0)) )" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Python generators can be defined in many contexts of the Computational Math, but here we show how they can be defined and manipulated to study a discrete dynamical system." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Python generators in the study of discrete dynamical systems" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "A discrete dynamical system on the space $\\mathbb{R}^m$, is defined usually by a one-to-one\n", "differentiable or only a continuous map, $f:\\mathbb{R}^m\\to\\mathbb{R}^m$. If $x_0\\in\\mathbb{R}^m$ is the initial state, its orbit or trajectory is the sequence $(x_n)_{n\\in\\mathbb{N}}$, $x_n=f^n(x_0)$, where $f^n$, $n\\geq 1$, is the $n^{th}$ iterate of the map.\n", "\n", "The goal of the dynamical system theory is to understand the asymptotic behaviour of the systems' orbits.\n", "\n", "In the case of a two-dimensional dynamical system we can plot a set of orbits, and the visual representation we get is called the *phase portrait* of that system." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "In order to illustrate the asymptotic behaviour of an orbit or to draw the phase portrait of a two-dimensional dynamical system, it is very useful to define a generator\n", "for an orbit.\n", "\n", "The generator will produce one point at a time or if it is necessary, we can take a slice in the infinite generator and store its points in a `list` or a `numpy.array`." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Experimenting with generator functions that generate points on an orbit, we deduced that there is a trick in the process of creation of a list of points resulted by iterating a slice in that generator.\n", "\n", "To illustrate how we were lead to that trick we define a generator of an orbit for a version of the [Hénon map](http://en.wikipedia.org/wiki/H%C3%A9non_map): \n", "\n", "$$(x,y)\\mapsto(y, 1-1.4 y^2+0.3x)$$\n", "\n" ] }, { "cell_type": "code", "execution_count": 29, "metadata": {}, "outputs": [], "source": [ "def gen_Henon(pt0): #pt0 is the initial point given as 2-list\n", " if not isinstance(pt0, list):\n", " raise ValueError(\"pt0 must be a tuple\")\n", " pt = pt0[:] # a copy of the initial point\n", " while True:\n", " yield pt\n", " tmp = pt[0]\n", " pt[0] = pt[1]\n", " pt[1] = 1 - 1.4 * pt[1]**2 + 0.3*tmp\n", " " ] }, { "cell_type": "code", "execution_count": 30, "metadata": {}, "outputs": [], "source": [ "pt0 = [0.0, 0.0]\n", "G = gen_Henon(pt0)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Let us iterate a slice in this generator:" ] }, { "cell_type": "code", "execution_count": 31, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "[-0.7408864000000001, 0.554322279213056]\n", "[0.554322279213056, 0.3475516150752599]\n", "[0.3475516150752599, 0.9971877085659265]\n" ] } ], "source": [ "for pt in islice(G, 5, 8):\n", " print (pt)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Usually to plot an orbit of some length we store its points in a list, and plot it with a single call of `plt.scatter`, instead of plotting one point at a time as it is generated.\n", "\n", "Trying to store the above three points in a list created by comprehension, we notice that the list\n", "contains only the last point, displayed three times:" ] }, { "cell_type": "code", "execution_count": 32, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "[0.0, 0.0]\n", "[[0.3475516150752599, 0.9971877085659265], [0.3475516150752599, 0.9971877085659265], [0.3475516150752599, 0.9971877085659265]]\n" ] } ], "source": [ "print(pt0)\n", "G = gen_Henon(pt0)# to reiterate the slice we have to call the generator function again\n", "L1 = [pt for pt in islice(G, 5, 8)]\n", "print(L1)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Why does the iteration work well when we only iterate the generator, but it fails when trying to create a list from items resulted from the iteration?" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The same drawback is manifested when we define the generator function with the argument\n", " `pt0` given as a `numpy.array`,\n", "and replace the line `pt=pt0[:]` with `pt=pt0.copy()`." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "After trial and error experimentation we got a correct list of items in the following way:" ] }, { "cell_type": "code", "execution_count": 33, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "[[-0.7408864000000001, 0.554322279213056], [0.554322279213056, 0.3475516150752599], [0.3475516150752599, 0.9971877085659265]]\n" ] } ], "source": [ "pt0 = [0.0, 0.0]\n", "G = gen_Henon(pt0)\n", "L3 = [list(pt) for pt in islice(G, 5, 8)]\n", "print(L3)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "This behaviour is really strange, because each `pt` was a list object and this conversion appears\n", "to be redundant." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "A similar right result we get converting each `pt` to a `numpy.array`:" ] }, { "cell_type": "code", "execution_count": 34, "metadata": {}, "outputs": [], "source": [ "import numpy as np" ] }, { "cell_type": "code", "execution_count": 35, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "[array([-0.7408864 , 0.55432228]), array([0.55432228, 0.34755162]), array([0.34755162, 0.99718771])]\n" ] } ], "source": [ "G = gen_Henon(pt0)\n", "L4 = [np.array(pt) for pt in islice(G, 5, 8)]\n", "print(L4)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "A third solution is to explicitly set the type of the object to be yielded by the generator:" ] }, { "cell_type": "code", "execution_count": 36, "metadata": {}, "outputs": [], "source": [ "def gen_Henon_new(pt0):\n", " if not isinstance(pt0, list):\n", " raise ValueError(\"pt0 must be a list\")\n", " pt = pt0[:]\n", " while True:\n", " yield pt\n", " tmp = pt[0]\n", " pt = [pt[1], 1-1.4*pt[1]**2+0.3*tmp]#pt is explicitly set as list " ] }, { "cell_type": "code", "execution_count": 37, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "[[-0.7408864000000001, 0.554322279213056], [0.554322279213056, 0.3475516150752599], [0.3475516150752599, 0.9971877085659265]]\n" ] } ], "source": [ "G = gen_Henon_new(pt0)\n", "L5 = [pt for pt in islice(G, 5,8)]\n", "print (L5)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We note that when the yielded object is explicitly set as list, no conversion is needed when we create the list of \n", " items resulted from the iteration." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Let us check this observed particularity for another generator." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Namely, we define the generator of pairs of consecutive terms in [Fibonacci sequence](http://en.wikipedia.org/wiki/Fibonacci_number): " ] }, { "cell_type": "code", "execution_count": 38, "metadata": {}, "outputs": [], "source": [ "def pair_fibonacci(a=1, b=1):\n", " while True:\n", " yield (a,b)# the yielded object is explicitly set as a tuple\n", " a, b = b, a+b\n" ] }, { "cell_type": "code", "execution_count": 39, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "[(1, 2), (2, 3), (3, 5), (5, 8), (8, 13), (13, 21), (21, 34), (34, 55), (55, 89), (89, 144), (144, 233), (233, 377), (377, 610), (610, 987)]\n" ] } ], "source": [ "genp = pair_fibonacci()\n", "L = [pair for pair in islice(genp, 1,15)]\n", "print (L)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The sequence of pairs $(p_n, q_n)$, $p_1=1, q_1=2$, produced by this generator, defines the convergents, \n", "$p_n/q_n$, of the continued fraction expansion of the inverse $1/\\phi$ for the [golden ratio](http://mathworld.wolfram.com/GoldenRatio.html)\n", "$\\phi=(1+\\sqrt{5})/2$.\n", "\n", "In the study of the [dynamics of twist maps](http://amath.colorado.edu/pub/dynamics/papers/Raglan2011_pt2.pdf), these pairs represent the rotation numbers\n", "of the periodic orbits, converging to the invariant 1-dimensional torus (or cantorus) having the rotation number,\n", "$1/\\phi$." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now we illustrate the phase portrait of the [standard map](http://www.scholarpedia.org/article/Chirikov_standard_map) whose orbits are generated by a generator function." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ " The standard map, $f$, is defined on the infinite annulus $\\mathbb{S}^1\\times \\mathbb{R}$.\n", "$\\mathbb{S}^1=\\mathbb{R}/\\mathbb{Z}$ is the unit circle identified with the interval $[0,1)$ or any other real \n", " interval of length 1.\n", " \n", "$f$ maps $(x,y)$ to $(x', y')$, where:\n", "\n", "$$ \\begin{align}\n", "x'&= x+y'\\: (modulo\\: 1)\\\\\n", "y'&= y-\\displaystyle\\frac{\\epsilon}{2\\pi}\\sin(2\\pi x), \\quad \\epsilon\\in\\mathbb{R}\n", "\\end{align}$$\n", "\n" ] }, { "cell_type": "code", "execution_count": 40, "metadata": {}, "outputs": [], "source": [ "def stdmap(pt):\n", " global eps\n", " pt[1] -= eps*np.sin(2*np.pi*pt[0])/(2.0*np.pi)\n", " pt[0 ]= np.mod(pt[0]+pt[1], 1.0)\n", " return pt" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We define an orbit generator, and a function that returns a `numpy.array` containing the points of an orbit:" ] }, { "cell_type": "code", "execution_count": 41, "metadata": {}, "outputs": [], "source": [ "def gen_stdmap(pt0) :\n", " if not isinstance(pt0, list):\n", " raise ValueError(\"pt0 must be a list\")\n", " pt = pt0[:]\n", " while True:\n", " yield pt\n", " pt = stdmap(pt) # here pt is not explicitly defined as a list" ] }, { "cell_type": "code", "execution_count": 42, "metadata": {}, "outputs": [], "source": [ "def orbit_stdmap(pt0, n):\n", " G = gen_stdmap(pt0)\n", " return np.array([list(next(G)) for _ in range(n)]) #conversion list(G.next()) is needed for\n", " # a right iteration" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Let us check first that the function `orbit_stdmap` works correctly:" ] }, { "cell_type": "code", "execution_count": 43, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "[[0. 0.57 ]\n", " [0.57 0.57 ]\n", " [0.20584273 0.63584273]\n", " [0.69295879 0.48711606]\n", " [0.32488936 0.63193057]\n", " [0.81898545 0.49409609]\n", " [0.45342134 0.63443589]\n", " [0.04324307 0.58982173]\n", " [0.59156338 0.54832031]\n", " [0.22402276 0.63245939]]\n" ] } ], "source": [ "eps = 0.971635# define the global parameter for the standard map\n", "pt0 = [0, 0.57]\n", "orb= orbit_stdmap(pt0, 10)\n", "print(orb)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "In the sequel we plot interactively a few orbits of the standard map, corresponding to the fixed parameter $\\epsilon$.\n", "Each orbit is colored with a color from a list of colors that is cyclically iterated.\n", "\n", "If we set the matplotlib backend `notebook`:\n", "\n", "`%matplotlib notebook` \n", "\n", "the phase portrait (i.e. the resulted figure) is inserted inline." ] }, { "cell_type": "code", "execution_count": 44, "metadata": {}, "outputs": [], "source": [ "from itertools import cycle" ] }, { "cell_type": "code", "execution_count": 64, "metadata": {}, "outputs": [], "source": [ "cols= ['#636efa',\n", " '#EF553B',\n", " '#00cc96',\n", " '#ab63fa',\n", " '#19d3f3',\n", " '#e763fa',\n", " '#FECB52',\n", " '#FF6692',\n", " '#B6E880']\n", "colors = cycle(cols)# colors= cyclic iterator" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The initial point of each orbit is chosen with a left mouse button click within the pop-up window\n", "that opens when running the cell below. \n", "\n", "When we have drawn a sufficient number of orbits such that \n", "the structure of the phase portrait is evident, \n", " we press the right button of the mouse and run the next cell.\n" ] }, { "cell_type": "code", "execution_count": 65, "metadata": {}, "outputs": [], "source": [ "%matplotlib notebook\n", "import matplotlib.pyplot as plt" ] }, { "cell_type": "code", "execution_count": 66, "metadata": {}, "outputs": [ { "data": { "application/javascript": [ "/* Put everything inside the global mpl namespace */\n", "window.mpl = {};\n", "\n", "\n", "mpl.get_websocket_type = function() {\n", " if (typeof(WebSocket) !== 'undefined') {\n", " return WebSocket;\n", " } else if (typeof(MozWebSocket) !== 'undefined') {\n", " return MozWebSocket;\n", " } else {\n", " alert('Your browser does not have WebSocket support.' +\n", " 'Please try Chrome, Safari or Firefox ≥ 6. ' +\n", " 'Firefox 4 and 5 are also supported but you ' +\n", " 'have to enable WebSockets in about:config.');\n", " };\n", "}\n", "\n", "mpl.figure = function(figure_id, websocket, ondownload, parent_element) {\n", " this.id = figure_id;\n", "\n", " this.ws = websocket;\n", "\n", " this.supports_binary = (this.ws.binaryType != undefined);\n", "\n", " if (!this.supports_binary) {\n", " var warnings = document.getElementById(\"mpl-warnings\");\n", " if (warnings) {\n", " warnings.style.display = 'block';\n", " warnings.textContent = (\n", " \"This browser does not support binary websocket messages. \" +\n", " \"Performance may be slow.\");\n", " }\n", " }\n", "\n", " this.imageObj = new Image();\n", "\n", " this.context = undefined;\n", " this.message = undefined;\n", " this.canvas = undefined;\n", " this.rubberband_canvas = undefined;\n", " this.rubberband_context = undefined;\n", " this.format_dropdown = undefined;\n", "\n", " this.image_mode = 'full';\n", "\n", " this.root = $('
');\n", " this._root_extra_style(this.root)\n", " this.root.attr('style', 'display: inline-block');\n", "\n", " $(parent_element).append(this.root);\n", "\n", " this._init_header(this);\n", " this._init_canvas(this);\n", " this._init_toolbar(this);\n", "\n", " var fig = this;\n", "\n", " this.waiting = false;\n", "\n", " this.ws.onopen = function () {\n", " fig.send_message(\"supports_binary\", {value: fig.supports_binary});\n", " fig.send_message(\"send_image_mode\", {});\n", " if (mpl.ratio != 1) {\n", " fig.send_message(\"set_dpi_ratio\", {'dpi_ratio': mpl.ratio});\n", " }\n", " fig.send_message(\"refresh\", {});\n", " }\n", "\n", " this.imageObj.onload = function() {\n", " if (fig.image_mode == 'full') {\n", " // Full images could contain transparency (where diff images\n", " // almost always do), so we need to clear the canvas so that\n", " // there is no ghosting.\n", " fig.context.clearRect(0, 0, fig.canvas.width, fig.canvas.height);\n", " }\n", " fig.context.drawImage(fig.imageObj, 0, 0);\n", " };\n", "\n", " this.imageObj.onunload = function() {\n", " fig.ws.close();\n", " }\n", "\n", " this.ws.onmessage = this._make_on_message_function(this);\n", "\n", " this.ondownload = ondownload;\n", "}\n", "\n", "mpl.figure.prototype._init_header = function() {\n", " var titlebar = $(\n", " '
');\n", " var titletext = $(\n", " '
');\n", " titlebar.append(titletext)\n", " this.root.append(titlebar);\n", " this.header = titletext[0];\n", "}\n", "\n", "\n", "\n", "mpl.figure.prototype._canvas_extra_style = function(canvas_div) {\n", "\n", "}\n", "\n", "\n", "mpl.figure.prototype._root_extra_style = function(canvas_div) {\n", "\n", "}\n", "\n", "mpl.figure.prototype._init_canvas = function() {\n", " var fig = this;\n", "\n", " var canvas_div = $('
');\n", "\n", " canvas_div.attr('style', 'position: relative; clear: both; outline: 0');\n", "\n", " function canvas_keyboard_event(event) {\n", " return fig.key_event(event, event['data']);\n", " }\n", "\n", " canvas_div.keydown('key_press', canvas_keyboard_event);\n", " canvas_div.keyup('key_release', canvas_keyboard_event);\n", " this.canvas_div = canvas_div\n", " this._canvas_extra_style(canvas_div)\n", " this.root.append(canvas_div);\n", "\n", " var canvas = $('');\n", " canvas.addClass('mpl-canvas');\n", " canvas.attr('style', \"left: 0; top: 0; z-index: 0; outline: 0\")\n", "\n", " this.canvas = canvas[0];\n", " this.context = canvas[0].getContext(\"2d\");\n", "\n", " var backingStore = this.context.backingStorePixelRatio ||\n", "\tthis.context.webkitBackingStorePixelRatio ||\n", "\tthis.context.mozBackingStorePixelRatio ||\n", "\tthis.context.msBackingStorePixelRatio ||\n", "\tthis.context.oBackingStorePixelRatio ||\n", "\tthis.context.backingStorePixelRatio || 1;\n", "\n", " mpl.ratio = (window.devicePixelRatio || 1) / backingStore;\n", "\n", " var rubberband = $('');\n", " rubberband.attr('style', \"position: absolute; left: 0; top: 0; z-index: 1;\")\n", "\n", " var pass_mouse_events = true;\n", "\n", " canvas_div.resizable({\n", " start: function(event, ui) {\n", " pass_mouse_events = false;\n", " },\n", " resize: function(event, ui) {\n", " fig.request_resize(ui.size.width, ui.size.height);\n", " },\n", " stop: function(event, ui) {\n", " pass_mouse_events = true;\n", " fig.request_resize(ui.size.width, ui.size.height);\n", " },\n", " });\n", "\n", " function mouse_event_fn(event) {\n", " if (pass_mouse_events)\n", " return fig.mouse_event(event, event['data']);\n", " }\n", "\n", " rubberband.mousedown('button_press', mouse_event_fn);\n", " rubberband.mouseup('button_release', mouse_event_fn);\n", " // Throttle sequential mouse events to 1 every 20ms.\n", " rubberband.mousemove('motion_notify', mouse_event_fn);\n", "\n", " rubberband.mouseenter('figure_enter', mouse_event_fn);\n", " rubberband.mouseleave('figure_leave', mouse_event_fn);\n", "\n", " canvas_div.on(\"wheel\", function (event) {\n", " event = event.originalEvent;\n", " event['data'] = 'scroll'\n", " if (event.deltaY < 0) {\n", " event.step = 1;\n", " } else {\n", " event.step = -1;\n", " }\n", " mouse_event_fn(event);\n", " });\n", "\n", " canvas_div.append(canvas);\n", " canvas_div.append(rubberband);\n", "\n", " this.rubberband = rubberband;\n", " this.rubberband_canvas = rubberband[0];\n", " this.rubberband_context = rubberband[0].getContext(\"2d\");\n", " this.rubberband_context.strokeStyle = \"#000000\";\n", "\n", " this._resize_canvas = function(width, height) {\n", " // Keep the size of the canvas, canvas container, and rubber band\n", " // canvas in synch.\n", " canvas_div.css('width', width)\n", " canvas_div.css('height', height)\n", "\n", " canvas.attr('width', width * mpl.ratio);\n", " canvas.attr('height', height * mpl.ratio);\n", " canvas.attr('style', 'width: ' + width + 'px; height: ' + height + 'px;');\n", "\n", " rubberband.attr('width', width);\n", " rubberband.attr('height', height);\n", " }\n", "\n", " // Set the figure to an initial 600x600px, this will subsequently be updated\n", " // upon first draw.\n", " this._resize_canvas(600, 600);\n", "\n", " // Disable right mouse context menu.\n", " $(this.rubberband_canvas).bind(\"contextmenu\",function(e){\n", " return false;\n", " });\n", "\n", " function set_focus () {\n", " canvas.focus();\n", " canvas_div.focus();\n", " }\n", "\n", " window.setTimeout(set_focus, 100);\n", "}\n", "\n", "mpl.figure.prototype._init_toolbar = function() {\n", " var fig = this;\n", "\n", " var nav_element = $('
')\n", " nav_element.attr('style', 'width: 100%');\n", " this.root.append(nav_element);\n", "\n", " // Define a callback function for later on.\n", " function toolbar_event(event) {\n", " return fig.toolbar_button_onclick(event['data']);\n", " }\n", " function toolbar_mouse_event(event) {\n", " return fig.toolbar_button_onmouseover(event['data']);\n", " }\n", "\n", " for(var toolbar_ind in mpl.toolbar_items) {\n", " var name = mpl.toolbar_items[toolbar_ind][0];\n", " var tooltip = mpl.toolbar_items[toolbar_ind][1];\n", " var image = mpl.toolbar_items[toolbar_ind][2];\n", " var method_name = mpl.toolbar_items[toolbar_ind][3];\n", "\n", " if (!name) {\n", " // put a spacer in here.\n", " continue;\n", " }\n", " var button = $('