Profiling calculation algorithms in Ginga's LineProfile plugin

LineProfile_profiling.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "from astropy.io import fits"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "filename = 'NRCN815A_LIN_20160115_uncal.fits'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "CPU times: user 14.1 ms, sys: 275 µs, total: 14.4 ms\n",
      "Wall time: 14.1 ms\n"
     ]
    }
   ],
   "source": [
    "%time bigcube = fits.getdata(filename)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(1, 190, 2048, 2048)"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "bigcube.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# A useful function from HSI project (zmath.py) since we cannot use Ginga canvas object here.\n",
    "def circular_mask(arr_shape, r, x_offset=0, y_offset=0):\n",
    "    \"\"\"Generate circular mask indices for 2D image.\n",
    "\n",
    "    Parameters\n",
    "    ----------\n",
    "    arr_shape : tuple of int\n",
    "        Shape of the array to use the mask.\n",
    "\n",
    "    r : number\n",
    "        Radius of the mask in pixels.\n",
    "\n",
    "    x_offset, y_offset : number, optional\n",
    "        Mask offset relative to image center.\n",
    "\n",
    "    Returns\n",
    "    -------\n",
    "    idx : tuple of arrays\n",
    "        Numpy indices of the mask, rounded to nearest\n",
    "        integer.\n",
    "\n",
    "    \"\"\"\n",
    "    # http://mail.scipy.org/pipermail/numpy-discussion/2011-January/054470.html\n",
    "\n",
    "    assert len(arr_shape) == 2, 'Image is not 2-D'\n",
    "\n",
    "    ny, nx = arr_shape\n",
    "    assert nx > 1 and ny > 1, 'Image is too small'\n",
    "\n",
    "    r = int(round(r))\n",
    "    assert r > 0, 'Radius must be int > 0'\n",
    "\n",
    "    xcen = np.round(0.5 * nx - 0.5 + x_offset).astype('int')\n",
    "    ycen = np.round(0.5 * ny - 0.5 + y_offset).astype('int')\n",
    "\n",
    "    x1, x2 = xcen - r, xcen + r\n",
    "    y1, y2 = ycen - r, ycen + r\n",
    "\n",
    "    y, x = np.ogrid[-r:r, -r:r]\n",
    "    i = np.where(x**2 + y**2 <= r**2)\n",
    "\n",
    "    if y1 >= 0 and y2 < ny and x1 >= 0 and x2 < nx:\n",
    "        # Mask contained in image.\n",
    "        a = np.zeros(arr_shape).astype('bool')\n",
    "\n",
    "        # populate array with mask position\n",
    "        a[y1:y2, x1:x2][i] = True\n",
    "\n",
    "        idx = np.where(a)\n",
    "\n",
    "    else:\n",
    "        # Mask falls outside image bounds.\n",
    "        # compute number of pixels mask extends beyond current array size\n",
    "        xout = np.abs(min(0 - x1,nx - x2))+1\n",
    "        yout = np.abs(min(0 - y1,ny - y2))+1\n",
    "\n",
    "        # derive size of new array that will contain entire mask\n",
    "        ny_new = ny + yout*2\n",
    "        nx_new = nx + xout*2\n",
    "        # initialize larger array\n",
    "        a = np.zeros((ny_new,nx_new)).astype('bool')\n",
    "\n",
    "        # recompute positions relative to larger array size\n",
    "        xcen = np.round(0.5 * nx_new - 0.5 + x_offset).astype('int')\n",
    "        ycen = np.round(0.5 * ny_new - 0.5 + y_offset).astype('int')\n",
    "\n",
    "        x1, x2 = xcen - r, xcen + r\n",
    "        y1, y2 = ycen - r, ycen + r\n",
    "\n",
    "        # populate array with mask position\n",
    "        a[y1:y2, x1:x2][i] = True\n",
    "\n",
    "        idx = np.where(a[yout:-yout,xout:-xout])\n",
    "\n",
    "    return idx"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "mask = np.zeros((2048, 2048), dtype=bool)\n",
    "i = circular_mask(mask.shape, 5, x_offset=100, y_offset=200)\n",
    "mask[i] = True"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(190, 2048, 2048)"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "data = bigcube[0]\n",
    "data.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "CPU times: user 7.73 s, sys: 1.95 s, total: 9.69 s\n",
      "Wall time: 51.3 s\n"
     ]
    }
   ],
   "source": [
    "%%time\n",
    "mask3d = np.zeros(data.shape, dtype=np.bool)\n",
    "mask3d[:, :, :] = mask[np.newaxis, :, :]  # z, y, x\n",
    "masked_data = np.ma.masked_array(data, mask=~mask3d)\n",
    "masked_mean = masked_data.mean(axis=(1, 2))\n",
    "plot_y_axis_data = masked_mean[~masked_mean.mask].data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "175 ms ± 599 µs per loop (mean ± std. dev. of 7 runs, 10 loops each)\n"
     ]
    }
   ],
   "source": [
    "%timeit result2 = [data[i][mask].mean() for i in range(data.shape[0])]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.testing.assert_allclose(result2, plot_y_axis_data, rtol=1e-6)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "172 ms ± 416 µs per loop (mean ± std. dev. of 7 runs, 10 loops each)\n"
     ]
    }
   ],
   "source": [
    "%timeit result3 = list(map(np.mean, [data[i][mask] for i in range(data.shape[0])]))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "np.testing.assert_allclose(result3, plot_y_axis_data, rtol=1e-6)"
   ]
  }
 ],
 "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
}