Merge pull request JendaPlhak#1 from Marigold/mojmir

mojmir's work
quapka · Jun 28, 2014 · 0d1a569 · 0d1a569
2 parents 612d500 + cc3b808
commit 0d1a569
Show file tree

Hide file tree

Showing 27 changed files with 365,950 additions and 0 deletions.
diff --git a/MojmirsWork/chaos_game.ipynb b/MojmirsWork/chaos_game.ipynb
diff --git a/MojmirsWork/data/maze1.txt b/MojmirsWork/data/maze1.txt
@@ -0,0 +1,4 @@
+012F
+1233
+1313
+S132
diff --git a/MojmirsWork/data/maze2.txt b/MojmirsWork/data/maze2.txt
@@ -0,0 +1,5 @@
+1112F
+11313
+00122
+12122
+S1321
diff --git a/MojmirsWork/feigenbaum.ipynb b/MojmirsWork/feigenbaum.ipynb
diff --git a/MojmirsWork/fractal_optimization.ipynb b/MojmirsWork/fractal_optimization.ipynb
@@ -0,0 +1,307 @@
+{
+ "metadata": {
+  "name": "",
+  "signature": "sha256:08617e7a508b2f1420c72fc6874d7bf350c1cc031e86aa6c179d6313fd255f72"
+ },
+ "nbformat": 3,
+ "nbformat_minor": 0,
+ "worksheets": [
+  {
+   "cells": [
+    {
+     "cell_type": "markdown",
+     "metadata": {},
+     "source": [
+      "# Fractal optimization\n",
+      "\n",
+      "The purpose of this notebook is to try to speed up fractal generation of mandelbrot set as much as possible. I got inspiration from this code in [Cython](https://github.com/cython/cython/wiki/examples-mandelbrot) and tried to write a code which would run as fast using [pypy](http://pypy.org/).\n",
+      "\n",
+      "*Note: The scripts below are not supposed to be run in iPython, use pypy instead*\n",
+      "\n",
+      "## PyPy\n",
+      "\n",
+      "From all the things I tried, those are the ones that had the greatest effect on speedup:\n",
+      "\n",
+      "* Image is not being created with `putpixel` method, but rather buffered into an array.\n",
+      "* Coloring is not done continously, but a palette of colors is pre-generated and we pick colors from it.\n",
+      "* Type `complex` was replaced with tuple of numbers to avoid expensive square root when calculating complex number size.\n",
+      "* Slicing of numpy arrays is replaced by assignment (see code).\n",
+      "\n",
+      "Although `pypy` should run without any code modifications, I got it working with decent speedup only after several hours due to various bugs and unexpected problems. It took at least 10x more time than writing the python code itself. Applying `pypy` on an original code `fractals.py` was as fast as regular python.\n",
+      "\n",
+      "## Parallelization\n",
+      "I tried to parallelize it using multiprocessing module by splitting images into blocks which are then processed separately and concatenated at last. Unfortunately, PyPy didn't generate an image due to problems with array pickling.\n",
+      "\n",
+      "## Cython\n",
+      "See this [blog post](https://github.com/cython/cython/wiki/examples-mandelbrot).\n",
+      "\n",
+      "## Comparison\n",
+      "Performance comparison is done for 500px X 500px image with 500 iterations, center in $-0.75-0.75i$ and radius $0.75$. I have an old `Intel\u00ae Core\u21222 Duo CPU P8600 @ 2.40GHz \u00d7 2`.\n",
+      "\n",
+      "|     Implementation    | Time   |\n",
+      "|:---------------------:|--------|\n",
+      "|       Python 2.7      | 190.4s |\n",
+      "|        PyPy 2.0       | 0.44s  |\n",
+      "|         Cython        | 0.29s  |\n",
+      "| PyPy 2.0, 2 processes | 0.36s  |\n",
+      "| PyPy 2.0, 4 processes | 0.39s  |"
+     ]
+    },
+    {
+     "cell_type": "markdown",
+     "metadata": {},
+     "source": [
+      "## Pallete of colors"
+     ]
+    },
+    {
+     "cell_type": "code",
+     "collapsed": false,
+     "input": [
+      "from __future__ import division\n",
+      "import Image\n",
+      "import colorsys\n",
+      "\n",
+      "colscale = lambda c: (int(c[0] * 255), int(c[1] * 255), int(c[2] * 255))\n",
+      "\n",
+      "def linspace(start, stop, n):\n",
+      "    if n == 1:\n",
+      "        yield stop\n",
+      "        return\n",
+      "    h = (stop - start) / float(n - 1)\n",
+      "    for i in range(n):\n",
+      "        yield start + h * i\n",
+      "\n",
+      "\n",
+      "def palette(n=100, rng=3, system='RGB', scaled=False):\n",
+      "    h = map(lambda x: x % 1, linspace(0 + 0.6, 3 + 0.6, n))\n",
+      "    s = [1.] * n\n",
+      "    v = map(lambda x: x**0.2, linspace(0, 1, n))\n",
+      "    if system == 'RGB':\n",
+      "        cols = [colorsys.hsv_to_rgb(*x) for x in zip(h, s, v)]\n",
+      "    else:\n",
+      "        cols = zip(h, s, v)\n",
+      "\n",
+      "    if scaled:\n",
+      "        cols = map(colscale, cols)\n",
+      "\n",
+      "    return cols"
+     ],
+     "language": "python",
+     "metadata": {},
+     "outputs": []
+    },
+    {
+     "cell_type": "markdown",
+     "metadata": {},
+     "source": [
+      "## PyPy script"
+     ]
+    },
+    {
+     "cell_type": "code",
+     "collapsed": false,
+     "input": [
+      "try:\n",
+      "    import numpypy as np\n",
+      "    def linspace(start, stop, n):\n",
+      "        if n == 1:\n",
+      "            yield stop\n",
+      "            return\n",
+      "        h = (stop - start) / (n - 1)\n",
+      "        for i in range(n):\n",
+      "            yield start + h * i\n",
+      "except ImportError:\n",
+      "    import numpy as np\n",
+      "    linspace = np.linspace\n",
+      "\n",
+      "import Image\n",
+      "import time\n",
+      "\n",
+      "PALETTE = palette(100, rng=2, system='RGB', scaled=True)\n",
+      "\n",
+      "def toimage(arr):\n",
+      "    data = np.asarray(arr)\n",
+      "    shape = list(data.shape)\n",
+      "    shape = (shape[1],shape[0])\n",
+      "    image = Image.fromstring('RGB', shape, data.tostring())\n",
+      "    return image\n",
+      "\n",
+      "\n",
+      "def mandelbrot(re, im, max_iter=500, escape=4.):\n",
+      "    z_im, z_re = 0., 0.\n",
+      "    i = 0\n",
+      "    while i < max_iter:\n",
+      "        z_re, z_im = (z_re * z_re - z_im * z_im + re, 2 * z_re * z_im + im)\n",
+      "        if z_re * z_re + z_im * z_im >= escape:\n",
+      "            return i\n",
+      "        i += 1\n",
+      "    return -1\n",
+      "\n",
+      "\n",
+      "def complex_image(size, center=(0,0), d=2):\n",
+      "\n",
+      "    byte_array = (np.zeros((size, size, 3))).astype('uint8')\n",
+      "\n",
+      "    palette_len = len(PALETTE)\n",
+      "\n",
+      "    for i, a in enumerate(linspace(center[0] - d, center[0] + d, size)):\n",
+      "        for j, b in enumerate(linspace(center[1] - d, center[1] + d, size)):\n",
+      "            k = mandelbrot(a, b)\n",
+      "\n",
+      "            if k != -1:\n",
+      "                byte_array[i, j, 0] = PALETTE[k % palette_len][0]\n",
+      "                byte_array[i, j, 1] = PALETTE[k % palette_len][1]\n",
+      "                byte_array[i, j, 2] = PALETTE[k % palette_len][2]\n",
+      "\n",
+      "    im = toimage(byte_array)\n",
+      "    return im\n",
+      "\n",
+      "\n",
+      "if __name__ == '__main__':\n",
+      "    start = time.time()\n",
+      "    im = complex_image(500, center=(-0.75, -0.75), d=0.75)\n",
+      "    print \"it took %f seconds to run\"%( time.time() - start)\n",
+      "\n",
+      "    im.show()\n",
+      "    ##im.save('mandel_fast.bmp')\n",
+      "else:\n",
+      "    import cProfile, pstats\n",
+      "\n",
+      "    stats = cProfile.runctx('complex_image(500, center=(-0.75, -0.75), d=0.75)',\n",
+      "        globals(), locals(), 'stats')\n",
+      "\n",
+      "    p = pstats.Stats('stats')\n",
+      "    p.strip_dirs().sort_stats('time').print_stats()"
+     ],
+     "language": "python",
+     "metadata": {},
+     "outputs": []
+    },
+    {
+     "cell_type": "markdown",
+     "metadata": {},
+     "source": [
+      "## PyPy multiprocessing"
+     ]
+    },
+    {
+     "cell_type": "code",
+     "collapsed": false,
+     "input": [
+      "from __future__ import division\n",
+      "\n",
+      "try:\n",
+      "    import numpypy as np\n",
+      "    def linspace(start, stop, n):\n",
+      "        h = (stop - start) / float(n - 1)\n",
+      "        return [start + h * i for i in range(n)]\n",
+      "\n",
+      "except ImportError:\n",
+      "    import numpy as np\n",
+      "    linspace = np.linspace\n",
+      "\n",
+      "import Image\n",
+      "import time\n",
+      "\n",
+      "import multiprocessing\n",
+      "\n",
+      "PALETTE = palette(100, rng=2, system='RGB', scaled=True)\n",
+      "\n",
+      "def toimage(arr):\n",
+      "    data = np.asarray(arr)\n",
+      "    shape = list(data.shape)\n",
+      "    shape = (shape[1],shape[0])\n",
+      "    image = Image.fromstring('RGB', shape, data.tostring())\n",
+      "    return image\n",
+      "\n",
+      "\n",
+      "def mandelbrot(re, im, max_iter=500, escape=4.):\n",
+      "    z_im, z_re = 0., 0.\n",
+      "    i = 0\n",
+      "    while i < max_iter:\n",
+      "        z_re, z_im = (z_re * z_re - z_im * z_im + re, 2 * z_re * z_im + im)\n",
+      "        if z_re * z_re + z_im * z_im >= escape:\n",
+      "            return i\n",
+      "        i += 1\n",
+      "    return -1\n",
+      "\n",
+      "\n",
+      "def complex_image_process(queue, size=(500,250), center=(-0.75, -0.75), d=0.75):\n",
+      "    height, width = size\n",
+      "    d_height, d_width = d, d * width / height\n",
+      "    byte_array = (np.zeros((height, width, 3))).astype('uint8')\n",
+      "    palette_len = len(PALETTE)\n",
+      "\n",
+      "    for i, a in enumerate(linspace(center[1] - d_height, center[1] + d_height, height)):\n",
+      "        for j, b in enumerate(linspace(center[0] - d_width, center[0] + d_width, width)):\n",
+      "            k = mandelbrot(a, b)\n",
+      "\n",
+      "            if k != -1:\n",
+      "                byte_array[i, j, 0] = PALETTE[k % palette_len][0]\n",
+      "                byte_array[i, j, 1] = PALETTE[k % palette_len][1]\n",
+      "                byte_array[i, j, 2] = PALETTE[k % palette_len][2]\n",
+      "\n",
+      "    #queue.put(byte_array)\n",
+      "\n",
+      "\n",
+      "def complex_image(size=500, center=(-0.75, -0.75), d=0.75, processes=2):\n",
+      "    assert size % processes == 0, 'Processes must divide image size.'\n",
+      "\n",
+      "    jobs = []\n",
+      "\n",
+      "    queue = multiprocessing.Queue()\n",
+      "\n",
+      "    for i in range(processes):\n",
+      "        re = center[0] - d + d / processes * (2*i + 1)\n",
+      "        kwargs = {'size': (size, int(size / processes)),\n",
+      "                  'center': (re, center[1]),\n",
+      "                  'd': d,\n",
+      "                  'queue': queue}\n",
+      "        p = multiprocessing.Process(target=complex_image_process, kwargs=kwargs)\n",
+      "        jobs.append(p)\n",
+      "        p.start()\n",
+      "\n",
+      "    for j in jobs:\n",
+      "        j.join()\n",
+      "\n",
+      "    block_list = []\n",
+      "    while not queue.empty():\n",
+      "        block_list.append(queue.get())\n",
+      "\n",
+      "    if block_list:\n",
+      "        b = np.hstack(block_list)\n",
+      "        im = toimage(b)\n",
+      "        return im\n",
+      "    else:\n",
+      "        return None\n",
+      "\n",
+      "\n",
+      "if __name__ == '__main__':\n",
+      "    start = time.time()\n",
+      "\n",
+      "    im = complex_image(500, center=(-0.75, -0.75), d=0.75, processes=2)\n",
+      "\n",
+      "    print \"it took %f seconds to run\" % ( time.time() - start)\n",
+      "\n",
+      "    if im:\n",
+      "        im.show()\n",
+      "\n",
+      "else:\n",
+      "    import cProfile, pstats\n",
+      "\n",
+      "    stats = cProfile.runctx('complex_image(500, center=(-0.75, -0.75), d=0.75)',\n",
+      "        globals(), locals(), 'stats')\n",
+      "\n",
+      "    p = pstats.Stats('stats')\n",
+      "    p.strip_dirs().sort_stats('time').print_stats()"
+     ],
+     "language": "python",
+     "metadata": {},
+     "outputs": []
+    }
+   ],
+   "metadata": {}
+  }
+ ]
+}