{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "Parallel Programming\n", "====\n", "\n", "The goal is to design parallel programs that are flexible, efficient and simple.\n", "\n", "**Step 0**: Start by profiling a serial program to identify bottlenecks\n", "\n", "**Step 1**: Are there for opportunities for parallelism?\n", "\n", "- Can tasks be performed in parallel?\n", " - Function calls\n", " - Loops\n", "- Can data be split and operated on in parallel?\n", " - Decomposition of arrays along rows, columns, blocks\n", " - Decomposition of trees into sub-trees\n", "- Is there a pipeline with a sequence of stages?\n", " - Data preprocessing and analysis\n", " - Graphics rendering\n", "\n", "**Step 2**: What is the nature of the parallelism?\n", "\n", "- Linear\n", " - Embarrassingly parallel programs\n", "- Recursive\n", " - Adaptive partitioning methods\n", "\n", "**Step 3**: What is the granularity?\n", "\n", "- 10s of jobs\n", "- 1000s of jobs\n", "\n", "**Step 4**: Choose an algorithm\n", "\n", "- Organize by tasks\n", " - Task parallelism\n", " - Divide and conquer\n", "\n", "- Organize by data\n", " - Geometric decomposition\n", " - Recursive decomposition\n", "\n", "- Organize by flow\n", " - Pipeline\n", " - Event-based processing\n", "\n", "**Step 5**: Map to program and data structures\n", "\n", "- Program structures\n", " - Single program multiple data (SPMD)\n", " - Master/worker\n", " - Loop parallelism\n", " - Fork/join\n", "- Data structures \n", " - Shared data\n", " - Shared queue\n", " - Distributed array\n", "\n", "**Step 6**: Map to parallel environment\n", "\n", "- Multi-core shared memory\n", " - Cython with OpenMP\n", " - multiprocessing\n", " - IPython.cluster\n", "- Multi-computer\n", " - IPython.cluster\n", " - MPI\n", " - Hadoop / Spark\n", "- GPU\n", " - CUDA\n", " - OpenCL\n", "\n", "**Step 7**: Execute, debug, tune in parallel environment" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Embarrassingly parallel programs\n", "----\n", "\n", "Many statistical problems are embarrassingly parallel and can be easily decomposed into independent tasks or data sets. Here are several examples:\n", "\n", "- Monte Carlo integration\n", "- Multiple chains of MCMC\n", "- Bootstrap for confidence intervals\n", "- Power calculations by simulation\n", "- Permutation-resampling tests \n", "- Fitting same model on multiple data sets\n", "\n", "Other problems are serial at small scale, but can be parallelized at large scales. For example, EM and MCMC iterations are inherently serial since there is a dependence on the previous state, but within a single iteration, there can be many thousands of density calculations (one for each data point to calculate the likelihood), and this is an embarrassingly parallel problem within a single iteration. \n", "\n", "These \"low hanging fruits\" are great because they offer a path to easy parallelism with minimal complexity." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Executing parallel code\n", "----\n", "\n", "**The bigger the problem, the more scope there is for parallelism**\n", "\n", "**Amhdahls' law** says that the speedup from parallelization is bounded by the ratio of parallelizable to irreducibly serial code in the algorithm. However, for big data analysis, **Gustafson's Law** is more relevant. This says that we are nearly always interested in increasing the size of the parallelizable bits, and the ratio of parallelizable to irreducibly serial code is not a static quantity but depends on data size. For example, Gibbs sampling has an irreducibly serial nature, but for large samples, each iteration may be able perform PDF evaluations in parallel for zillions of data points." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Threading\n", "----\n", "\n", "The Python interpreter only allows one thread to execute at any one time. This is known as the Global Interpreter Lock (GIL) and generally makes the use of threads an ineffective strategy for CPU-bound tasks. However, threads are useful in I/O bound tasks such as reading/writing/downloading files since threads give up their lock once an I//O task has been initiated. " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### CPU-bound task\n", "\n", "**Estimating $\\pi$ using Monte Carlo integration**\n", "\n", "This is clearly a toy example, but the template can be used for most embarrassingly parallel problems and is a good example of a CPU-bound task." ] }, { "cell_type": "code", "execution_count": 1, "metadata": { "collapsed": true }, "outputs": [], "source": [ "def mc_pi(n):\n", " s = 0\n", " for i in range(n):\n", " x = np.random.uniform(-1, 1)\n", " y = np.random.uniform(-1, 1)\n", " if (x**2 + y**2) < 1:\n", " s += 1\n", " return 4*s/n" ] }, { "cell_type": "code", "execution_count": 2, "metadata": { "collapsed": false }, "outputs": [], "source": [ "n = int(1e7)" ] }, { "cell_type": "code", "execution_count": 3, "metadata": { "collapsed": false }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "CPU times: user 15.8 s, sys: 21.1 ms, total: 15.8 s\n", "Wall time: 15.8 s\n" ] }, { "data": { "text/plain": [ "3.1417504" ] }, "execution_count": 3, "metadata": {}, "output_type": "execute_result" } ], "source": [ "%%time\n", "mc_pi(n)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Optimizations without parallelism\n", "\n", "Review of techniques to accelerate Python code." ] }, { "cell_type": "code", "execution_count": 4, "metadata": { "collapsed": false }, "outputs": [], "source": [ "def mc_pi1(n):\n", " x = np.random.uniform(-1, 1, (n,2))\n", " return 4*np.sum((x**2).sum(1) < 1)/n" ] }, { "cell_type": "code", "execution_count": 5, "metadata": { "collapsed": false }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "CPU times: user 786 ms, sys: 300 ms, total: 1.09 s\n", "Wall time: 1.08 s\n" ] }, { "data": { "text/plain": [ "3.1411228000000002" ] }, "execution_count": 5, "metadata": {}, "output_type": "execute_result" } ], "source": [ "%%time \n", "mc_pi1(n)" ] }, { "cell_type": "code", "execution_count": 6, "metadata": { "collapsed": true }, "outputs": [], "source": [ "from numba import jit, njit" ] }, { "cell_type": "code", "execution_count": 7, "metadata": { "collapsed": true }, "outputs": [], "source": [ "@njit()\n", "def mc_pi2(n):\n", " s = 0\n", " for i in range(n):\n", " x = np.random.uniform(-1, 1)\n", " y = np.random.uniform(-1, 1)\n", " if (x**2 + y**2) < 1:\n", " s += 1\n", " return 4*s/n" ] }, { "cell_type": "code", "execution_count": 8, "metadata": { "collapsed": false }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "CPU times: user 639 ms, sys: 9.36 ms, total: 648 ms\n", "Wall time: 645 ms\n" ] }, { "data": { "text/plain": [ "3.1420776" ] }, "execution_count": 8, "metadata": {}, "output_type": "execute_result" } ], "source": [ "%%time\n", "mc_pi2(n)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Cython\n", "\n", "It is a little more work to optimize the Cython version. We do not use numpy.random as that is optimized for generating arrays rather than single numbers. Instead we use the GNU Scientific Library (wrapped for us in the package cython_gsl) to draw uniform random variates.\n", "\n", "```\n", "pip install cythongsl\n", "```" ] }, { "cell_type": "code", "execution_count": 9, "metadata": { "collapsed": true }, "outputs": [], "source": [ "%load_ext cython" ] }, { "cell_type": "code", "execution_count": 10, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/html": [ "\n", "\n", "\n", "\n", " \n", " Cython: _cython_magic_73ea126a291c7fb73e3c8fcbef013dcf.pyx\n", " \n", " \n", "\n", "\n", "

Generated by Cython 0.23.4

\n", "

\n", " Yellow lines hint at Python interaction.
\n", " Click on a line that starts with a \"+\" to see the C code that Cython generated for it.\n", "

\n", "
 01:
\n", "
 02: import cython
\n", "
 03: from cython_gsl cimport *
\n", "
 04:
\n", "
 05: @cython.cdivision(True)
\n", "
+06: def mc_pi3(int n):
\n", "
/* Python wrapper */\n",
       "static PyObject *__pyx_pw_46_cython_magic_73ea126a291c7fb73e3c8fcbef013dcf_1mc_pi3(PyObject *__pyx_self, PyObject *__pyx_arg_n); /*proto*/\n",
       "static PyMethodDef __pyx_mdef_46_cython_magic_73ea126a291c7fb73e3c8fcbef013dcf_1mc_pi3 = {\"mc_pi3\", (PyCFunction)__pyx_pw_46_cython_magic_73ea126a291c7fb73e3c8fcbef013dcf_1mc_pi3, METH_O, 0};\n",
       "static PyObject *__pyx_pw_46_cython_magic_73ea126a291c7fb73e3c8fcbef013dcf_1mc_pi3(PyObject *__pyx_self, PyObject *__pyx_arg_n) {\n",
       "  int __pyx_v_n;\n",
       "  PyObject *__pyx_r = 0;\n",
       "  __Pyx_RefNannyDeclarations\n",
       "  __Pyx_RefNannySetupContext(\"mc_pi3 (wrapper)\", 0);\n",
       "  assert(__pyx_arg_n); {\n",
       "    __pyx_v_n = __Pyx_PyInt_As_int(__pyx_arg_n); if (unlikely((__pyx_v_n == (int)-1) && PyErr_Occurred())) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 6; __pyx_clineno = __LINE__; goto __pyx_L3_error;}\n",
       "  }\n",
       "  goto __pyx_L4_argument_unpacking_done;\n",
       "  __pyx_L3_error:;\n",
       "  __Pyx_AddTraceback(\"_cython_magic_73ea126a291c7fb73e3c8fcbef013dcf.mc_pi3\", __pyx_clineno, __pyx_lineno, __pyx_filename);\n",
       "  __Pyx_RefNannyFinishContext();\n",
       "  return NULL;\n",
       "  __pyx_L4_argument_unpacking_done:;\n",
       "  __pyx_r = __pyx_pf_46_cython_magic_73ea126a291c7fb73e3c8fcbef013dcf_mc_pi3(__pyx_self, ((int)__pyx_v_n));\n",
       "  int __pyx_lineno = 0;\n",
       "  const char *__pyx_filename = NULL;\n",
       "  int __pyx_clineno = 0;\n",
       "\n",
       "  /* function exit code */\n",
       "  __Pyx_RefNannyFinishContext();\n",
       "  return __pyx_r;\n",
       "}\n",
       "\n",
       "static PyObject *__pyx_pf_46_cython_magic_73ea126a291c7fb73e3c8fcbef013dcf_mc_pi3(CYTHON_UNUSED PyObject *__pyx_self, int __pyx_v_n) {\n",
       "  gsl_rng_type *__pyx_v_T;\n",
       "  gsl_rng *__pyx_v_r;\n",
       "  double __pyx_v_s;\n",
       "  double __pyx_v_x;\n",
       "  double __pyx_v_y;\n",
       "  CYTHON_UNUSED int __pyx_v_i;\n",
       "  PyObject *__pyx_r = NULL;\n",
       "  __Pyx_RefNannyDeclarations\n",
       "  __Pyx_RefNannySetupContext(\"mc_pi3\", 0);\n",
       "/* … */\n",
       "  /* function exit code */\n",
       "  __pyx_L1_error:;\n",
       "  __Pyx_XDECREF(__pyx_t_4);\n",
       "  __Pyx_AddTraceback(\"_cython_magic_73ea126a291c7fb73e3c8fcbef013dcf.mc_pi3\", __pyx_clineno, __pyx_lineno, __pyx_filename);\n",
       "  __pyx_r = NULL;\n",
       "  __pyx_L0:;\n",
       "  __Pyx_XGIVEREF(__pyx_r);\n",
       "  __Pyx_RefNannyFinishContext();\n",
       "  return __pyx_r;\n",
       "}\n",
       "/* … */\n",
       "  __pyx_tuple_ = PyTuple_Pack(8, __pyx_n_s_n, __pyx_n_s_n, __pyx_n_s_T, __pyx_n_s_r, __pyx_n_s_s, __pyx_n_s_x, __pyx_n_s_y, __pyx_n_s_i); if (unlikely(!__pyx_tuple_)) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 6; __pyx_clineno = __LINE__; goto __pyx_L1_error;}\n",
       "  __Pyx_GOTREF(__pyx_tuple_);\n",
       "  __Pyx_GIVEREF(__pyx_tuple_);\n",
       "/* … */\n",
       "  __pyx_t_1 = PyCFunction_NewEx(&__pyx_mdef_46_cython_magic_73ea126a291c7fb73e3c8fcbef013dcf_1mc_pi3, NULL, __pyx_n_s_cython_magic_73ea126a291c7fb73e); if (unlikely(!__pyx_t_1)) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 6; __pyx_clineno = __LINE__; goto __pyx_L1_error;}\n",
       "  __Pyx_GOTREF(__pyx_t_1);\n",
       "  if (PyDict_SetItem(__pyx_d, __pyx_n_s_mc_pi3, __pyx_t_1) < 0) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 6; __pyx_clineno = __LINE__; goto __pyx_L1_error;}\n",
       "  __Pyx_DECREF(__pyx_t_1); __pyx_t_1 = 0;\n",
       "
 07:     cdef gsl_rng_type * T
\n", "
 08:     cdef gsl_rng * r
\n", "
+09:     cdef double s = 0.0
\n", "
  __pyx_v_s = 0.0;\n",
       "
 10:     cdef double x, y
\n", "
 11:     cdef int i
\n", "
 12:
\n", "
+13:     gsl_rng_env_setup()
\n", "
  gsl_rng_env_setup();\n",
       "
 14:
\n", "
+15:     T = gsl_rng_default
\n", "
  __pyx_v_T = gsl_rng_default;\n",
       "
+16:     r = gsl_rng_alloc (T)
\n", "
  __pyx_v_r = gsl_rng_alloc(__pyx_v_T);\n",
       "
 17:
\n", "
+18:     for i in range(n):
\n", "
  __pyx_t_1 = __pyx_v_n;\n",
       "  for (__pyx_t_2 = 0; __pyx_t_2 < __pyx_t_1; __pyx_t_2+=1) {\n",
       "    __pyx_v_i = __pyx_t_2;\n",
       "
+19:         x = 2*gsl_rng_uniform(r) - 1
\n", "
    __pyx_v_x = ((2.0 * gsl_rng_uniform(__pyx_v_r)) - 1.0);\n",
       "
+20:         y = 2*gsl_rng_uniform(r)- 1
\n", "
    __pyx_v_y = ((2.0 * gsl_rng_uniform(__pyx_v_r)) - 1.0);\n",
       "
+21:         if (x**2 + y**2) < 1:
\n", "
    __pyx_t_3 = (((pow(__pyx_v_x, 2.0) + pow(__pyx_v_y, 2.0)) < 1.0) != 0);\n",
       "    if (__pyx_t_3) {\n",
       "/* … */\n",
       "    }\n",
       "  }\n",
       "
+22:             s += 1
\n", "
      __pyx_v_s = (__pyx_v_s + 1.0);\n",
       "
+23:     return 4*s/n
\n", "
  __Pyx_XDECREF(__pyx_r);\n",
       "  __pyx_t_4 = PyFloat_FromDouble(((4.0 * __pyx_v_s) / __pyx_v_n)); if (unlikely(!__pyx_t_4)) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 23; __pyx_clineno = __LINE__; goto __pyx_L1_error;}\n",
       "  __Pyx_GOTREF(__pyx_t_4);\n",
       "  __pyx_r = __pyx_t_4;\n",
       "  __pyx_t_4 = 0;\n",
       "  goto __pyx_L0;\n",
       "
" ], "text/plain": [ "" ] }, "execution_count": 10, "metadata": {}, "output_type": "execute_result" } ], "source": [ "%%cython -a -lgsl\n", "\n", "import cython\n", "from cython_gsl cimport *\n", "\n", "@cython.cdivision(True) \n", "def mc_pi3(int n):\n", " cdef gsl_rng_type * T\n", " cdef gsl_rng * r\n", " cdef double s = 0.0\n", " cdef double x, y\n", " cdef int i\n", "\n", " gsl_rng_env_setup()\n", "\n", " T = gsl_rng_default\n", " r = gsl_rng_alloc (T)\n", "\n", " for i in range(n):\n", " x = 2*gsl_rng_uniform(r) - 1\n", " y = 2*gsl_rng_uniform(r)- 1\n", " if (x**2 + y**2) < 1:\n", " s += 1\n", " return 4*s/n" ] }, { "cell_type": "code", "execution_count": 11, "metadata": { "collapsed": false }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "CPU times: user 308 ms, sys: 1.39 ms, total: 309 ms\n", "Wall time: 308 ms\n" ] }, { "data": { "text/plain": [ "3.1414584" ] }, "execution_count": 11, "metadata": {}, "output_type": "execute_result" } ], "source": [ "%%time\n", "mc_pi3(n)" ] }, { "cell_type": "markdown", "metadata": { "collapsed": true }, "source": [ "Parallelization with threads\n", "----\n", "\n", "We reduce the number of Monte Carlo samples by a factor of 10, but run 10 such jobs independent. The result returned will be an generator expression with 10 values.. We then average to get the results. Notice that threading does not help much, if at all." ] }, { "cell_type": "code", "execution_count": 12, "metadata": { "collapsed": true }, "outputs": [], "source": [ "from concurrent.futures import ThreadPoolExecutor" ] }, { "cell_type": "code", "execution_count": 13, "metadata": { "collapsed": false }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "3.1427279999999995\n", "CPU times: user 311 ms, sys: 7.19 ms, total: 319 ms\n", "Wall time: 313 ms\n" ] } ], "source": [ "%%time\n", "njobs = 10\n", "with ThreadPoolExecutor(max_workers=4) as pool:\n", " res = pool.map(mc_pi3, [int(1e6)]*njobs)\n", "print(sum(res)/njobs)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Parallelization with processes\n", "\n", "We switch a multi-process model. Unlike threads, processes do run in parallel. Now you should see some speedup from parallelization." ] }, { "cell_type": "code", "execution_count": 14, "metadata": { "collapsed": true }, "outputs": [], "source": [ "from concurrent.futures import ProcessPoolExecutor" ] }, { "cell_type": "code", "execution_count": 15, "metadata": { "collapsed": false }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "3.1427279999999995\n", "CPU times: user 15.2 ms, sys: 31.5 ms, total: 46.6 ms\n", "Wall time: 138 ms\n" ] } ], "source": [ "%%time\n", "njobs = 10\n", "with ProcessPoolExecutor(max_workers=4) as pool:\n", " res = pool.map(mc_pi3, [int(1e6)]*njobs)\n", "print(sum(res)/njobs)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Using `multiprocessing`\n", "\n", "There is another module for parallel processes in the standard library that provides more data structures and control - for example, if you need the processes to share some data or communicate in some way. See [official documentation](https://docs.python.org/3/library/multiprocessing.html) for details." ] }, { "cell_type": "code", "execution_count": 16, "metadata": { "collapsed": true }, "outputs": [], "source": [ "from multiprocessing import cpu_count, Pool" ] }, { "cell_type": "code", "execution_count": 17, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "8" ] }, "execution_count": 17, "metadata": {}, "output_type": "execute_result" } ], "source": [ "cpu_count()" ] }, { "cell_type": "code", "execution_count": 18, "metadata": { "collapsed": false }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "3.1427279999999995\n", "CPU times: user 13.7 ms, sys: 43.6 ms, total: 57.3 ms\n", "Wall time: 116 ms\n" ] } ], "source": [ "%%time\n", "pool = Pool()\n", "res = pool.map(mc_pi3, [int(1e6)]*njobs)\n", "pool.close()\n", "print(sum(res)/njobs)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Running functions with multiple arguments\n", "\n", "If a function f takes multiple arguments, we can use `startmap`. The `starmap` method acts as though it applies f(*args) to the iterable of args passed to the worker pool." ] }, { "cell_type": "code", "execution_count": 19, "metadata": { "collapsed": true }, "outputs": [], "source": [ "def add(a, b):\n", " return a + b" ] }, { "cell_type": "code", "execution_count": 20, "metadata": { "collapsed": true }, "outputs": [], "source": [ "x = np.arange(10)\n", "y = x[::-1]" ] }, { "cell_type": "code", "execution_count": 21, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "[(0, 9),\n", " (1, 8),\n", " (2, 7),\n", " (3, 6),\n", " (4, 5),\n", " (5, 4),\n", " (6, 3),\n", " (7, 2),\n", " (8, 1),\n", " (9, 0)]" ] }, "execution_count": 21, "metadata": {}, "output_type": "execute_result" } ], "source": [ "list(zip(x, y))" ] }, { "cell_type": "code", "execution_count": 22, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "[9, 9, 9, 9, 9, 9, 9, 9, 9, 9]" ] }, "execution_count": 22, "metadata": {}, "output_type": "execute_result" } ], "source": [ "pool = Pool()\n", "res = pool.starmap(add, list(zip(x, y)))\n", "pool.close()\n", "res" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Using `ipyparallel`\n", "\n", "Parallel execution is tightly integrated with Jupyter in the `ipyparallel` package. Install with\n", "\n", "```bash\n", "conda install ipyparallel\n", "ipcluster nbextension enable\n", "```\n", "\n", "[Official documentation](https://ipyparallel.readthedocs.org/en/latest/)\n", "\n", "You can start a cluster on the `IPython Clusters` tab." ] }, { "cell_type": "code", "execution_count": 23, "metadata": { "collapsed": false }, "outputs": [], "source": [ "from ipyparallel import Client" ] }, { "cell_type": "code", "execution_count": 24, "metadata": { "collapsed": false }, "outputs": [], "source": [ "rc = Client()" ] }, { "cell_type": "code", "execution_count": 25, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "[0, 1, 2, 3, 4, 5, 6, 7]" ] }, "execution_count": 25, "metadata": {}, "output_type": "execute_result" } ], "source": [ "rc.ids" ] }, { "cell_type": "code", "execution_count": 26, "metadata": { "collapsed": false }, "outputs": [], "source": [ "dv = rc[:]" ] }, { "cell_type": "code", "execution_count": 27, "metadata": { "collapsed": false }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "3.1417084\n" ] } ], "source": [ "res = dv.map_sync(mc_pi1, [int(1e6)] * 10)\n", "print(sum(res)/njobs)" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": true }, "outputs": [], "source": [] } ], "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.5.1" } }, "nbformat": 4, "nbformat_minor": 0 }