FFI.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# Interfacing between Python, C and Fortran"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Foreign Function Interface (FFI)\n",
    "There are a variety of reasons we might wish to interface between code written in C or Fortran when using Python. The most common reasons (and best reasons in my opinion) are:\n",
    "- Speed - Both Fortran and C are (significantly) faster than Fortran for numerical programming\n",
    "- Existing code - If we have code that already does the job, but we just want to incorporate it into python code (or make it easier to use).\n",
    "\n",
    "Note that I'm using Python 3 here (hence some of the strings being declared as `b'string'` - this is because Python 3 uses unicode by default whereas Python 2 uses ASCII).\n",
    "\n",
    "Python has a built-in (part of the standard library) module `ctypes`"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## ctypes"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": true,
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "import ctypes\n",
    "from sys import platform\n",
    "from datetime import datetime\n",
    "core_libraries = {\n",
    "    'darwin': 'libc.dylib',\n",
    "    'linux': 'libc.so',\n",
    "    'linux2': 'libc.so',\n",
    "    'cygwin': 'libc.so',\n",
    "    'win32': 'msvcrt' \n",
    "}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "We can use c data types with python types, and extract the value with `.value`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "c_ushort(65526) 65526 65526\n"
     ]
    }
   ],
   "source": [
    "print(ctypes.c_ushort(-10), 2**16 - 10, ctypes.c_ushort(-10).value)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "c_wchar_p(4449400792) Hello, World\n"
     ]
    }
   ],
   "source": [
    "s = \"Hello, World\"\n",
    "c_string = ctypes.c_wchar_p(s)\n",
    "print(c_string, c_string.value) # wraps around the pointer, but can still access the underlying string"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "We can call functions from dynamic c libraries:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "<CDLL 'libc.dylib', handle 10d9a9770 at 0x109347f28>\n",
      "2017-12-13 11:31:26\n"
     ]
    }
   ],
   "source": [
    "corelib = ctypes.cdll.LoadLibrary(core_libraries[platform])\n",
    "print(corelib)\n",
    "t = corelib.time(None) # Note do NOT call corelib.time() or you'll crash\n",
    "print(datetime.fromtimestamp(t)) "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## cffi"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "35"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "from cffi import FFI\n",
    "ffi = FFI()\n",
    "ffi.cdef(\"\"\"\n",
    "    int printf(const char *format, ...);   // copy-pasted from the man page\n",
    "\"\"\")\n",
    "C = ffi.dlopen(None)                     # loads the entire C namespace\n",
    "arg = ffi.new(\"char[]\", b'ffi created char array')  # equivalent to C code: char arg[] = \"ffi created char array\";\n",
    "C.printf(b\"we can print %s\", arg)         # call printf"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "We can compile our own library and wrapper"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": true,
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "ffibuilder = FFI()\n",
    "\n",
    "ffibuilder.set_source(\"_example\", # name of our library\n",
    "   r\"\"\" // passed to the real C compiler\n",
    "        #include <sys/types.h>\n",
    "        #include <pwd.h>\n",
    "    \"\"\",\n",
    "    libraries=[])   # or a list of libraries to link with\n",
    "    # (more arguments like setup.py's Extension class:\n",
    "    # include_dirs=[..], extra_objects=[..], and so on)\n",
    "\n",
    "ffibuilder.cdef(\"\"\"     // some declarations from the man page\n",
    "    struct passwd {\n",
    "        char *pw_name;\n",
    "        ...;     // literally dot-dot-dot\n",
    "    };\n",
    "    struct passwd *getpwuid(int uid);\n",
    "\"\"\")\n",
    "\n",
    "lib_location = ffibuilder.compile(verbose=False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "_example.c                     _example.o\r\n",
      "\u001b[31m_example.cpython-36m-darwin.so\u001b[m\u001b[m\r\n"
     ]
    }
   ],
   "source": [
    "# if we list the directory we will see the files here\n",
    "!ls _example*\n",
    "# you can see the contents of _example.c"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "b'root'\n",
      "b'daemon'\n"
     ]
    }
   ],
   "source": [
    "from _example import ffi, lib\n",
    "\n",
    "p = lib.getpwuid(0) # user with uid 0 should be root\n",
    "# p is now a c struct (passwd struct as defined above)\n",
    "print(ffi.string(p.pw_name))\n",
    "print(ffi.string(lib.getpwuid(1).pw_name))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "<cdata 'struct passwd *' 0x7fc61d76a050>\n"
     ]
    }
   ],
   "source": [
    "print(p)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "These modules can be used to wrap fortran code, once we have generated c interfaces using `iso_c_binding`. In order to make this more straightforward for some use cases, `f2py` can be used."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## f2py"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "\n",
    "In order to use f2py you probably need to have installed `numpy` and/or `scipy`, which is very likely already the case if you actually want to use Python with Fortran. For more complete documentation have a look [here](https://docs.scipy.org/doc/numpy-dev/f2py/python-usage.html)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Overwriting mod.f90\n"
     ]
    }
   ],
   "source": [
    "%%writefile mod.f90\n",
    "module mod\n",
    "  integer :: i\n",
    "  integer :: x(4)\n",
    "  real, dimension(2,3) :: a\n",
    "  \n",
    "contains\n",
    "  subroutine foo\n",
    "    a(1,2) = a(1,2)+3\n",
    "  end subroutine foo\n",
    "  \n",
    "  subroutine assumed(arr, n)\n",
    "  integer, intent(in) :: n\n",
    "  real(8), dimension(n), intent(inout) :: arr\n",
    "    arr = arr + 1\n",
    "  end subroutine assumed\n",
    "\n",
    "  subroutine assumed2d(arr, m, n)\n",
    "  integer, intent(in) :: m, n\n",
    "  real(8), dimension(m, n), intent(inout) :: arr\n",
    "    arr = arr * 1.5\n",
    "  end subroutine assumed2d\n",
    "\n",
    "  subroutine assumed_shape(arr)\n",
    "    real(8), intent(inout) :: arr(:,:)\n",
    "    arr = arr * 3.1415\n",
    "  end subroutine assumed_shape\n",
    "end module mod"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "collapsed": true,
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "# wrap the module using f2py \n",
    "# (note that on my system I've got multiple f2py versions so I called f2py3.6 rather than just f2py)\n",
    "!f2py3.6 -c -m mod mod.f90 -DF2PY_REPORT_ON_ARRAY_COPY=1 2>&1 > f2py.log"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "assumed(arr,[n])\n",
      "\n",
      "Wrapper for ``assumed``.\n",
      "\n",
      "Parameters\n",
      "----------\n",
      "arr : in/output rank-1 array('d') with bounds (n)\n",
      "\n",
      "Other Parameters\n",
      "----------------\n",
      "n : input int, optional\n",
      "    Default: len(arr)\n",
      "\n"
     ]
    }
   ],
   "source": [
    "from mod import mod\n",
    "import numpy as np\n",
    "print(mod.assumed.__doc__)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([1, 2, 0, 0], dtype=int32)"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "mod.x[:2] = [1,2]\n",
    "mod.a = [[1,2,3],[4,5,6]]\n",
    "mod.x"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[ 1.,  2.,  3.],\n",
       "       [ 4.,  5.,  6.]], dtype=float32)"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "mod.a "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {
    "collapsed": true,
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "mod.foo()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[ 1.,  5.,  3.],\n",
       "       [ 4.,  5.,  6.]], dtype=float32)"
      ]
     },
     "execution_count": 16,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "mod.a"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[ 2.  3.  4.  5.]\n"
     ]
    }
   ],
   "source": [
    "a = np.array([1,2,3,4], dtype=np.float64)\n",
    "mod.assumed(a, len(a))\n",
    "print(a)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[ 1.5  1.5  1.5  1.5  1.5]\n",
      " [ 1.5  1.5  1.5  1.5  1.5]\n",
      " [ 1.5  1.5  1.5  1.5  1.5]\n",
      " [ 1.5  1.5  1.5  1.5  1.5]]\n"
     ]
    }
   ],
   "source": [
    "x = np.asfortranarray(np.ones(20, dtype=np.float64).reshape((4,5)))\n",
    "mod.assumed2d(x, 4, 5)\n",
    "print(x)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[ 4.71224999  4.71224999  4.71224999  4.71224999  4.71224999]\n",
      " [ 4.71224999  4.71224999  4.71224999  4.71224999  4.71224999]\n",
      " [ 4.71224999  4.71224999  4.71224999  4.71224999  4.71224999]\n",
      " [ 4.71224999  4.71224999  4.71224999  4.71224999  4.71224999]]\n"
     ]
    }
   ],
   "source": [
    "mod.assumed_shape(x)\n",
    "print(x)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "f2py does NOT support derived types natively, but something like [f90wrap](https://github.com/jameskermode/f90wrap) might be useful here"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "For an example of a fortran wrapped module look at [python-dftd3](https://github.com/peterspackman/python-dftd3)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Further reading\n",
    "http://lucumr.pocoo.org/2013/8/18/beautiful-native-libraries/\n",
    "\n",
    "https://github.com/bast/python-cffi-demo"
   ]
  }
 ],
 "metadata": {
  "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.6.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}

Idioms.ipynb