This notebook explains how to perfrom azimuthal integration only with numpy and scipy and how to get neverthless nice performances

naive_integration.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Azimuthal integration in Python\n",
    "\n",
    "## Common initialization step:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy\n",
    "npt = 1024\n",
    "y,x  = numpy.ogrid[-512:512,-512:512]\n",
    "radius = (x*x+y*y)**0.5\n",
    "rmax = radius.max()+0.1\n",
    "data = numpy.random.random((1024,1024))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Naive approach integration using corona extraction with masks:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "CPU times: user 1.05 s, sys: 0 ns, total: 1.05 s\n",
      "Wall time: 1.05 s\n"
     ]
    }
   ],
   "source": [
    "%%time\n",
    "res_loop = numpy.zeros(npt)\n",
    "for i in range(npt):\n",
    "    rinf = rmax * i / npt\n",
    "    rsup = rinf + rmax / npt\n",
    "    mask = numpy.logical_and((rinf <= radius),(radius < rsup))\n",
    "    res_loop[i] = data[mask].mean()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Vectorized version using histograms:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "CPU times: user 20.5 ms, sys: 218 µs, total: 20.7 ms\n",
      "Wall time: 19.1 ms\n"
     ]
    }
   ],
   "source": [
    "%%time\n",
    "count_of_pixels = numpy.histogram(radius, npt, range=[0,rmax] )[0]\n",
    "sum_of_intensities = numpy.histogram(radius, npt, weights=data, range=[0,rmax])[0]\n",
    "res_vec = sum_of_intensities / count_of_pixels"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "True"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Speed-up: 50x, validation:\n",
    "numpy.allclose(res_loop, res_vec)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Using a Sparse matrix multiplication\n",
    "\n",
    "Those multiplication can take advantage of parallel hardware unlike histogram which require costly *atomic* operations.\n",
    "This trick is called \"scatter to gather\" transformation in parallel programming. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(1024, 1048576)\n",
      "CPU times: user 60.5 ms, sys: 6.21 ms, total: 66.7 ms\n",
      "Wall time: 69.7 ms\n"
     ]
    }
   ],
   "source": [
    "%%time \n",
    "from scipy.sparse import csc_matrix\n",
    "positions = numpy.histogram(radius, npt, range=[0,rmax] )[1]\n",
    "row = numpy.digitize(radius.ravel(), positions) - 1\n",
    "size = row.size\n",
    "col = numpy.arange(size)\n",
    "dat = numpy.ones(size, dtype=float)\n",
    "csr = csc_matrix((dat, (row, col)), shape = (npt, data.size))\n",
    "print(csr.shape)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "CPU times: user 3.11 ms, sys: 3.1 ms, total: 6.21 ms\n",
      "Wall time: 4.69 ms\n"
     ]
    }
   ],
   "source": [
    "%%time\n",
    "count_csr = csr.dot(numpy.ones(data.size))\n",
    "sum_csr = csr.dot(data.ravel())\n",
    "res_csr = sum_csr / count_csr"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "True"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Speed-up: 5x vs histograms, validation:\n",
    "numpy.allclose(res_csr, res_vec)"
   ]
  }
 ],
 "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.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}