8-low-level.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# Low-level hackery"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "_(nothing to install; just numpy, matplotlib from before)_"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "One function I've used without much comment is `numpy.frombuffer`, which lets us wrap arbitrary regions of memory as Numpy arrays. We can \"peek\" at any memory we want; we can also \"poke\" it, changing values, byte by byte."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "Consider, for instance, a byte string. These are immutable (cannot be changed) in Python:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "hello = b\"Hello, world!\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "try:\n",
    "    hello[4:8] = b\"????\"\n",
    "except TypeError as err:\n",
    "    print(\"Nope: \" + str(err))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "import numpy\n",
    "a = numpy.frombuffer(hello, dtype=numpy.uint8)\n",
    "a"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "a.view(\"S1\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "By default, Numpy tries to protect you from doing evil things."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "try:\n",
    "    a[4:8] = [69, 86, 73, 76]\n",
    "except ValueError as err:\n",
    "    print(\"Nope: \" + str(err))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "But this is Python: we can shoot our feet if we want to."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "a.flags.writeable = True"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "a[4:8] = [69, 86, 73, 76]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "hello"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "This messes with Python's internal data model."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "hello = b\"Hello, world!\"\n",
    "a = numpy.frombuffer(hello, dtype=numpy.uint8)\n",
    "a.flags.writeable = True\n",
    "a[4:8] = [69, 86, 73, 76]\n",
    "print(hello == b\"Hello, world!\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "exec(\"\"\"\n",
    "hello = b\"Hello, world!\"\n",
    "a = numpy.frombuffer(hello, dtype=numpy.uint8)\n",
    "a.flags.writeable = True\n",
    "a[4:8] = [69, 86, 73, 76]\n",
    "print(hello == b\"Hello, world!\")\n",
    "\"\"\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "(The second example was interpreted as a `.pyc` script, in which all instances of the literal `b\"Hello, world!\"` were replaced by a single object: modifying that object in line 4 changed it in line 5!)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "With the help of ctypes, a built-in Python library, Numpy can wrap any address at all. (Some will cause segmentation faults, so be careful!)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "x = 12345"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import ctypes\n",
    "import sys\n",
    "\n",
    "ptr = ctypes.cast(id(x), ctypes.POINTER(ctypes.c_uint8))\n",
    "a = numpy.ctypeslib.as_array(ptr, (sys.getsizeof(x),))\n",
    "a"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "We're looking at a Python object header, a pointer to the `int` type (also a Python object), and then the number itself: `12345` in little endian bytes is `57, 48, 0, 0`. Do you see it?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "a[-4:].view(numpy.int32)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "Let's try a string."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "y = \"Hey there.\"\n",
    "ptr = ctypes.cast(id(y), ctypes.POINTER(ctypes.c_uint8))\n",
    "a = numpy.ctypeslib.as_array(ptr, (sys.getsizeof(y),))\n",
    "a"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "a[-11:].tostring()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "By wrapping a pointer as a Numpy array (and maybe setting `array.flags.writeable = True`), we can do anything."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "We can break anything."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "Seriously, just letting `numpy.ctypeslib.as_array(0, (8,))` repr itself on the screen would cause a segmentation fault. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "Is there anything useful we can do with this power?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "The original purpose of the ctypes library was to run code in compiled C libraries (not C++, that's much more complicated).\n",
    "\n",
    "For a busy data analyst, the advantage of that is that you don't need to write (or wait for) official bindings to use a C library."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import ctypes\n",
    "libc = ctypes.cdll.LoadLibrary(\"/lib/x86_64-linux-gnu/libc.so.6\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "We can run arbitrary functions from `libc.so.6`, but we have to tell Python what its argument types are (that's not stored in the shared object file)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "libc.malloc.argtypes = (ctypes.c_size_t,)                      # argument types (only one)\n",
    "libc.malloc.restype = ctypes.POINTER(ctypes.c_double)          # return type"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "ptr = libc.malloc(100 * numpy.dtype(numpy.float64).itemsize)   # pass number of bytes\n",
    "ptr"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "a = numpy.ctypeslib.as_array(ptr, (100,))\n",
    "a"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "a[:] = 0.0\n",
    "a"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "a.flags"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "We've just used low-level `libc.malloc` to allocate data for an array. This array doesn't \"own\" its memory, so it doesn't go away when the Python object (`a`) gets garbage collected. We'd have to call `libc.free`, like any C programmer."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "Why do this? I've used it to allocate arrays on NUMA hardware and Knight's Landing MCDRAM. Specialty memory allocations will probably get more important, not less, as architectures get more heterogeneous."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "(No reason to write an entire analysis in C just to get special allocators.)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "Another possible reason is to access special devices; device drivers are often written in C and distributed with C programmers in mind, but they may be the only thing between you and an important dataset."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "Another example: you can wrap structs. (The snake eats its own tail again.)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "class PyObject(ctypes.Structure): pass\n",
    "PyObject._fields_ = [(\"ob_refcnt\", ctypes.c_size_t),\n",
    "                     (\"ob_type\", ctypes.POINTER(PyObject))]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "which is equivalent to\n",
    "\n",
    "```c\n",
    "struct PyObject {\n",
    "    size_t ob_refcnt;\n",
    "    PyObject* ob_type;\n",
    "    // the rest depends on the type of object\n",
    "}\n",
    "```"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "hello = b\"Hello, world!\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ptr = PyObject.from_address(id(hello))\n",
    "ptr"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "ptr.ob_refcnt"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "This `ob_refcnt` is the number of Python references to a given object. There's a way to do it with a `sys` call:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sys.getrefcount(hello)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "but it's always one too high because you create a reference to pass it to that function!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "biglist = [hello] * 1000"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ptr.ob_refcnt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sys.getrefcount(hello)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "The ctypes library can wrap numbers, strings, pointers, arrays, and structs, which is just about everything you might encounter in C. The reason C++ isn't supported is because C shoehorns its much larger type system into shared object files by \"name mangling.\"\n",
    "\n",
    "Even if you reverse the name mangling with the `c++filt` program, those names cannot be uniquely identified without interpreting the C++ headers. By that point, you might as well use Cython or pybind11."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "**One last example:** the quickest, dirtiest way possible to call out to compiled code:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import os\n",
    "with open(\"tmp.cpp\", \"w\") as cfile:\n",
    "    cfile.write(\"\"\"\n",
    "#include <complex>\n",
    "extern \"C\" {    \n",
    "    void quick(int height, int width, int maxiterations, double* c, int* fractal) {\n",
    "        for (int h = 0;  h < height;  h++) {\n",
    "            for (int w = 0;  w < width;  w++) {\n",
    "                double creal = c[2 * (h + height*w)];\n",
    "                double cimag = c[2 * (h + height*w) + 1];\n",
    "                std::complex<double> ci = std::complex<double>(creal, cimag);\n",
    "                std::complex<double> z = ci;\n",
    "                for (int i = 0;  i < maxiterations;  i++) {\n",
    "                    z = z * z + ci;\n",
    "                    if (std::abs(z) > 2) {\n",
    "                        fractal[h + height*w] = i;\n",
    "                        break;\n",
    "                    }\n",
    "                }\n",
    "            }\n",
    "        }\n",
    "    }\n",
    "}\n",
    "\"\"\")\n",
    "assert os.system(\"gcc -O3 -fPIC -shared tmp.cpp -o libtmp.so\") == 0"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "libtmp = ctypes.cdll.LoadLibrary(os.path.join(os.getcwd(), \"libtmp.so\"))\n",
    "libtmp.quick.argtypes = (ctypes.c_int, ctypes.c_int, ctypes.c_int,\n",
    "                         ctypes.POINTER(ctypes.c_double),\n",
    "                         ctypes.POINTER(ctypes.c_int))\n",
    "libtmp.quick.restype = None"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def run_dirty(height, width, maxiterations=20):\n",
    "    y, x = numpy.ogrid[-1:0:height*1j, -1.5:0:width*1j]\n",
    "    c = x + y*1j\n",
    "    fractal = numpy.zeros(c.shape, dtype=numpy.int32) + maxiterations\n",
    "    libtmp.quick(height, width, maxiterations,\n",
    "                 ctypes.cast(c.ctypes.data, ctypes.POINTER(ctypes.c_double)),\n",
    "                 ctypes.cast(fractal.ctypes.data, ctypes.POINTER(ctypes.c_int)))\n",
    "    return fractal"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "import time\n",
    "starttime = time.time()\n",
    "fractal = run_dirty(8000, 12000)\n",
    "time.time() - starttime"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "import matplotlib.pyplot as plt\n",
    "fig, ax = plt.subplots(figsize=(12, 8))\n",
    "ax.imshow(fractal)\n",
    "# ax.imshow(fractal[-200:, :300])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "Not very good time, but the right answer."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "With that monstrosity, I'll end the course.\n",
    "\n",
    "Cheers!"
   ]
  }
 ],
 "metadata": {
  "celltoolbar": "Slideshow",
  "kernelspec": {
   "display_name": "Python 3",
   "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.7.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}