{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "[Oregon Curriculum Network](http://www.4dsolutions.net/ocn)
\n", "[Discovering Math with Python](Introduction.ipynb)\n", "\n", "# The Mandelbrot Set\n", "\n", "Welcome to our Jupyter Notebook about the Mandelbrot Set. We will be taking the complex number c at each point, and adding it to z \\* z, where z starts out 0. z = z * z + c. Do that over and over, up to 80 times in this case. It's a survivor's game in that most points subjected to this algorithm spiral out of the grid after a few iterations. Others hang in there. At the end, we image how all the points did i.e. we colorize by how many iterations each survived, before having too big a magnitude.\n", "\n", "You may not remember how complex numbers multiply: imagined as arrows pointing from the origin, we multiply the arrow lengths and add their angles. Adding complex numbers is like adding two vectors, tip to tail.\n", "\n", "What are \"vectorized\" operations? That's where entire data structures act like individual objects. Every element is affected, potentially, but without our needing to write looping code to visit every cell individually. Lets start with a simple example..." ] }, { "cell_type": "code", "execution_count": 1, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "array([[1, 2, 3],\n", " [4, 5, 6]])" ] }, "execution_count": 1, "metadata": {}, "output_type": "execute_result" } ], "source": [ "import numpy as np\n", "\n", "r1 = np.array([[1,2,3],[4,5,6]])\n", "r2 = np.ones((2,3)) # the (2,3) is the \"shape\" of the array we want, 2 rows, 3 columns\n", "r1" ] }, { "cell_type": "code", "execution_count": 2, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "array([[ 1., 1., 1.],\n", " [ 1., 1., 1.]])" ] }, "execution_count": 2, "metadata": {}, "output_type": "execute_result" } ], "source": [ "r2" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Both these arrays have the same shape. We may add them...." ] }, { "cell_type": "code", "execution_count": 3, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "array([[ 2., 3., 4.],\n", " [ 5., 6., 7.]])" ] }, "execution_count": 3, "metadata": {}, "output_type": "execute_result" } ], "source": [ "r1 + r2" ] }, { "cell_type": "code", "execution_count": 4, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "array([[ 2., 5., 10.],\n", " [ 17., 26., 37.]])" ] }, "execution_count": 4, "metadata": {}, "output_type": "execute_result" } ], "source": [ "r1 ** 2 + r2 # we may also multipy r1 by itself (not matrix multiplication) and add r2" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The function below is all about vectorizing. Both z and output are 1000 x 1000 arrays! " ] }, { "cell_type": "code", "execution_count": 5, "metadata": { "collapsed": true }, "outputs": [], "source": [ "def mandelbrot_numpy(c, maxiter):\n", " output = np.zeros(c.shape)\n", " z = np.zeros(c.shape, np.complex64)\n", " for it in range(maxiter):\n", " # not done are the numbers that haven't spiralled away by this many iterations\n", " notdone = np.less(z.real*z.real + z.imag*z.imag, 4.0) # magnitude test!\n", " output[notdone] = it # lets keep track of how many iterations in output\n", " z[notdone] = z[notdone]**2 + c[notdone] # on to the next iteration!\n", " output[output == maxiter-1] = 0 # these made it to the end\n", " return output" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Lets look at a smaller slice of the complex plane to being with. This rectangle of complex numbers starts from (-1, -1j) and ends at (1, 1j), dividing both the real and imaginary axes into three intervals marked by four tick marks." ] }, { "cell_type": "code", "execution_count": 6, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "array([[-0.-1.j ],\n", " [-0.-0.33333334j],\n", " [ 0.+0.33333334j],\n", " [ 0.+1.j ]], dtype=complex64)" ] }, "execution_count": 6, "metadata": {}, "output_type": "execute_result" } ], "source": [ "r1 = np.linspace(-1, 1, 4, dtype=np.float32)\n", "r2 = np.linspace(-1, 1, 4, dtype=np.float32)\n", "r2_T = (r2 * 1j).reshape(4,1)\n", "r2_T" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "By reshaping r2 into r2_T, we make it a vertical column. Then we add. The result is a 4x4 matrix, from (-1, -1) upper left, to (1, 1) lower right." ] }, { "cell_type": "code", "execution_count": 7, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "array([[-1.00000000-1.j , -0.33333334-1.j ,\n", " 0.33333334-1.j , 1.00000000-1.j ],\n", " [-1.00000000-0.33333334j, -0.33333334-0.33333334j,\n", " 0.33333334-0.33333334j, 1.00000000-0.33333334j],\n", " [-1.00000000+0.33333334j, -0.33333334+0.33333334j,\n", " 0.33333334+0.33333334j, 1.00000000+0.33333334j],\n", " [-1.00000000+1.j , -0.33333334+1.j ,\n", " 0.33333334+1.j , 1.00000000+1.j ]], dtype=complex64)" ] }, "execution_count": 7, "metadata": {}, "output_type": "execute_result" } ], "source": [ "r1 + r2_T" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Notice how we were able to add rectangular arrays r1 and r2_T in one line, with no loops." ] }, { "cell_type": "code", "execution_count": 8, "metadata": { "collapsed": true }, "outputs": [], "source": [ "def mandelbrot_set(xmin,xmax,ymin,ymax,width,height,maxiter):\n", " r1 = np.linspace(xmin, xmax, width, dtype=np.float32)\n", " r2 = np.linspace(ymin, ymax, height, dtype=np.float32)\n", " c = r1 + (r2 * 1j).reshape(width,1) # x row + each y column\n", " n3 = mandelbrot_numpy(c, maxiter)\n", " return (r1,r2,n3) " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The Python standard library has a timeit module. Jupyter Notebooks add a layer of %magic by giving us the %timeit directive. We're able to see how long the operation takes." ] }, { "cell_type": "code", "execution_count": 9, "metadata": { "collapsed": false }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "1 loop, best of 3: 1.32 s per loop\n" ] } ], "source": [ "%timeit mandelbrot_set(-2.0,0.5,-1.25,1.25,1000,1000,80)" ] }, { "cell_type": "code", "execution_count": 10, "metadata": { "collapsed": false }, "outputs": [], "source": [ "r1, r2, result = mandelbrot_set(-2.0,0.5,-1.25,1.25,1000,1000,80)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We now have the array we want to plot in result." ] }, { "cell_type": "code", "execution_count": 11, "metadata": { "collapsed": false }, "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", " this.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 = $('