Do some profiling of image combination

profile-image-combination.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import os\n",
    "import bottleneck as bn\n",
    "from ccdproc import combine, ImageFileCollection, CCDData\n",
    "from astropy.stats import median_absolute_deviation\n",
    "from astropy.io import fits\n",
    "from  pathlib import Path\n",
    "import numpy as np"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "current = os.getcwd()\n",
    "\n",
    "file_dir = 'foo4'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "os.mkdir(file_dir)\n",
    "os.chdir(file_dir)\n",
    "size = [2024, 2031]\n",
    "n_images = 25\n",
    "\n",
    "base_name = 'test-combine-{num:03d}.fits'\n",
    "\n",
    "for num in range(n_images):\n",
    "    data = np.random.normal(size=size)\n",
    "    # Now add some outlying pixels so there is something to clip\n",
    "    n_bad = 50000\n",
    "    bad_x = np.random.randint(0, high=size[0] - 1, size=n_bad)\n",
    "    bad_y = np.random.randint(0, high=size[1] - 1, size=n_bad)\n",
    "    data[bad_x, bad_y] = np.random.choice([-1, 1], size=n_bad) * (10 + np.random.rand(n_bad))\n",
    "    hdu = fits.PrimaryHDU(data=np.asarray(data, dtype='float32'))\n",
    "    hdu.header['for_prof'] = 'yes'\n",
    "    hdu.header['bunit'] = 'adu'\n",
    "    hdu.writeto(base_name.format(num=num), overwrite=True)\n",
    "\n",
    "os.chdir(current)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "ic = ImageFileCollection(file_dir, '*')\n",
    "fils = ic.files_filtered(for_prof='yes', include_path=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## clip with: numpy masked median, astropy median_absolute_deviation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      " \n",
      "*** Profile stats marshalled to file 'profile_fun_mamedian_dev_mask_med_abs_dev.prof'. \n"
     ]
    }
   ],
   "source": [
    "%%prun -D profile_fun_mamedian_dev_mask_med_abs_dev.prof\n",
    "combine(fils, sigma_clip=True, \n",
    "        sigma_clip_low_thresh=5, sigma_clip_high_thresh=5, \n",
    "        sigma_clip_func=np.ma.median, sigma_clip_dev_func=median_absolute_deviation)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## clip with: numpy unmasked median, astropy median_absolute_deviation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/Users/mcraig/conda-main/envs/ccdp-dev/lib/python3.6/site-packages/numpy/core/fromnumeric.py:664: UserWarning: Warning: 'partition' will ignore the 'mask' of the MaskedArray.\n",
      "  a.partition(kth, axis=axis, kind=kind, order=order)\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      " \n",
      "*** Profile stats marshalled to file 'profile_fun_median_dev_med_abs_dev.prof'. \n"
     ]
    }
   ],
   "source": [
    "%%prun -D profile_fun_median_dev_med_abs_dev.prof\n",
    "combine(fils, sigma_clip=True, \n",
    "        sigma_clip_low_thresh=5, sigma_clip_high_thresh=5, \n",
    "        sigma_clip_func=np.median, sigma_clip_dev_func=median_absolute_deviation)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## clip with: numpy unmasked median, numpy masked std"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/Users/mcraig/conda-main/envs/ccdp-dev/lib/python3.6/site-packages/numpy/core/fromnumeric.py:664: UserWarning: Warning: 'partition' will ignore the 'mask' of the MaskedArray.\n",
      "  a.partition(kth, axis=axis, kind=kind, order=order)\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      " \n",
      "*** Profile stats marshalled to file 'profile_fun_median_dev_mask_std.prof'. \n"
     ]
    }
   ],
   "source": [
    "%%prun -D profile_fun_median_dev_mask_std.prof\n",
    "combine(fils, sigma_clip=True, \n",
    "        sigma_clip_low_thresh=5, sigma_clip_high_thresh=5, \n",
    "        sigma_clip_func=np.median, sigma_clip_dev_func=np.ma.std)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## clip with: numpy unmasked median, numpy unmasked std"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/Users/mcraig/conda-main/envs/ccdp-dev/lib/python3.6/site-packages/numpy/core/fromnumeric.py:664: UserWarning: Warning: 'partition' will ignore the 'mask' of the MaskedArray.\n",
      "  a.partition(kth, axis=axis, kind=kind, order=order)\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      " \n",
      "*** Profile stats marshalled to file 'profile_fun_median_dev_std.prof'. \n"
     ]
    }
   ],
   "source": [
    "%%prun -D profile_fun_median_dev_std.prof\n",
    "combine(fils, sigma_clip=True, \n",
    "        sigma_clip_low_thresh=5, sigma_clip_high_thresh=5, \n",
    "        sigma_clip_func=np.median, sigma_clip_dev_func=np.std)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## clip with: bottleneck median, bottleneck nanstd"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      " \n",
      "*** Profile stats marshalled to file 'profile_bottleneck_dev_nanstd.prof'. \n"
     ]
    }
   ],
   "source": [
    "%%prun -D profile_bottleneck_dev_nanstd.prof\n",
    "combine(fils, sigma_clip=True, \n",
    "        sigma_clip_low_thresh=5, sigma_clip_high_thresh=5, \n",
    "        sigma_clip_func=bn.median, sigma_clip_dev_func=bn.nanstd)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## clip with: bottleneck median, custom mad_std"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "def bottle_mad_std(data, axis=0):\n",
    "    data_median = bn.median(data, axis=0)\n",
    "    data_median = np.expand_dims(data_median, axis=axis)\n",
    "    return bn.median(np.abs(data - data_median), axis=0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      " \n",
      "*** Profile stats marshalled to file 'profile_bottleneck_dev_custom.prof'. \n"
     ]
    }
   ],
   "source": [
    "%%prun -D profile_bottleneck_dev_custom.prof\n",
    "combine(fils, sigma_clip=True, \n",
    "        sigma_clip_low_thresh=5, sigma_clip_high_thresh=5, \n",
    "        sigma_clip_func=bn.median, sigma_clip_dev_func=bottle_mad_std)"
   ]
  }
 ],
 "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.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}