{
 "metadata": {
  "name": "",
  "signature": "sha256:1c86b687bfabaaedabefe03c147b1c71b50eafd5b35fe011e57ba8d464df715d"
 },
 "nbformat": 3,
 "nbformat_minor": 0,
 "worksheets": [
  {
   "cells": [
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "%pylab nbagg"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "Populating the interactive namespace from numpy and matplotlib\n"
       ]
      }
     ],
     "prompt_number": 1
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "import os\n",
      "#os.environ[\"PYOPENCL_CTX\"] = \"0:1\"\n",
      "os.environ[\"PYOPENCL_COMPILER_OUTPUT\"]=\"1\""
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 2
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "import pyopencl,numpy\n",
      "from pyopencl import array\n",
      "shape=1024,1024\n",
      "ctx=pyopencl.create_some_context()\n",
      "queue=pyopencl.CommandQueue(ctx, properties=pyopencl.command_queue_properties.PROFILING_ENABLE)\n",
      "data=numpy.exp(20*numpy.random.random(shape)**2)\n",
      "x = numpy.random.randint(0, shape[1], size=10*shape[1])\n",
      "y = numpy.random.randint(0, shape[0], size=10*shape[0])\n",
      "data[y,x] = numpy.NaN\n",
      "datag=array.to_device(queue, data.astype(\"float32\"))\n",
      "prg=pyopencl.Program(ctx, open(\"../pyFAI/resources/openCL/sigma_clip.cl\").read()).build();\n",
      "m=array.empty(queue,1024,\"float32\");\n",
      "d=array.empty(queue,1024,\"float32\");\n",
      "ws=shape[0]//8;\n",
      "local_mem = pyopencl.LocalMemory(ws * 20);\n",
      "print(abs(numpy.nanmean(data, dtype=\"float32\",axis=0)-numpy.nanmean(data, dtype=\"float64\",axis=0)).max(),\n",
      "      abs(numpy.nanmean(data, dtype=\"float32\",axis=-1)-numpy.nanmean(data, dtype=\"float64\",axis=-1)).max(),\n",
      "      abs(numpy.nanstd(data, dtype=\"float32\",axis=0)-numpy.nanstd(data, dtype=\"float64\",axis=0)).max(), \n",
      "      abs(numpy.nanstd(data, dtype=\"float32\",axis=-1)-numpy.nanstd(data, dtype=\"float64\",axis=-1)).max())\n"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "(28.443244431167841, 27.278686532750726, 326.02058730274439, 319.114636272192)\n"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stderr",
       "text": [
        "/usr/lib/python2.7/dist-packages/pkg_resources/__init__.py:1869: UserWarning: /usr/lib/pymodules/python2.7/rpl-1.5.5.egg-info could not be properly decoded in UTF-8\n",
        "  warnings.warn(msg)\n"
       ]
      }
     ],
     "prompt_number": 3
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "e=prg.mean_std_horizontal(queue, (1024, 1024//8), (1, 1024//8), datag.data, m.data, d.data, numpy.float32(-1), local_mem); \n",
      "print(abs(numpy.nanmean(data, axis=1, dtype=\"float64\")-m.get()).max(), abs(numpy.nanstd(data, axis=1, dtype=\"float64\")-d.get()).max())\n",
      "print(1e-6 * (e.profile.end - e.profile.start),\"ms\")\n",
      "print(m)\n",
      "print(d)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "(3.2379647828638554, 9.2297753766179085)\n",
        "(0.109792, 'ms')\n",
        "[ 13357314.  15665482.  10535890. ...,  14345746.  11944847.  11272129.]\n",
        "[ 55772528.  59189568.  44902948. ...,  58260868.  51066760.  49634728.]\n"
       ]
      }
     ],
     "prompt_number": 4
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "e = prg.mean_std_vertical(queue, (1024//8, 1024), (1024//8, 1), datag.data, m.data, d.data, numpy.float32(-1), local_mem)\n",
      "print(abs(numpy.nanmean(data, axis=0, dtype=\"float64\")-m.get()).max(),abs(numpy.nanstd(data, axis=0, dtype=\"float64\")-d.get()).max())\n",
      "print(1e-6 * (e.profile.end - e.profile.start),\"ms\")\n",
      "print(m)\n",
      "print(d)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "(2.6332576386630535, 9.4150413945317268)\n",
        "(0.233376, 'ms')\n",
        "[  9952937.  13959638.  15037384. ...,  12514131.  13894042.  13005303.]\n",
        "[ 46625852.  56019904.  60556288. ...,  54266764.  58167156.  55265980.]\n"
       ]
      }
     ],
     "prompt_number": 5
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "import warnings\n",
      "as_strided = numpy.lib.stride_tricks.as_strided\n",
      "\n",
      "def sigma_clip(image, sigma_lo=3, sigma_hi=3, max_iter=5, axis=0):\n",
      "    image=image.copy()\n",
      "    mask = numpy.logical_not(numpy.isfinite(image))\n",
      "    dummies = mask.sum()\n",
      "    image[mask] = numpy.NaN\n",
      "    mean = numpy.nanmean(image, axis=axis)\n",
      "    std = numpy.nanstd(image, axis=axis)\n",
      "    for i in range(max_iter):\n",
      "        if axis==0:\n",
      "            mean2d = as_strided(mean, image.shape, (0, mean.strides[0]))\n",
      "            std2d = as_strided(std, image.shape, (0, std.strides[0]))\n",
      "        else:\n",
      "            mean2d = as_strided(mean, image.shape, (mean.strides[0], 0))\n",
      "            std2d = as_strided(std, image.shape, (std.strides[0], 0))\n",
      "        with warnings.catch_warnings():\n",
      "           warnings.simplefilter(\"ignore\")\n",
      "           delta = (image - mean2d) / std2d\n",
      "           mask = numpy.logical_or(delta > sigma_hi,\n",
      "                                   delta < -sigma_lo)\n",
      "        dummies = mask.sum()\n",
      "        if dummies == 0:\n",
      "            break\n",
      "        image[mask] = numpy.NaN\n",
      "        mean = numpy.nanmean(image, axis=axis)\n",
      "        std = numpy.nanstd(image, axis=axis)\n",
      "    return mean, std"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 6
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "#work on the vertical direction\n",
      "%timeit sigma_clip(data, max_iter=5)\n",
      "mean_n, std_n = sigma_clip(data, max_iter=5)\n",
      "print(mean_n)\n",
      "print(std_n)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "10 loops, best of 3: 117 ms per loop\n",
        "[  92733.44112401  159155.93516633  158456.13784525 ...,  159711.39281794\n",
        "  159685.13550263  191007.00782656]"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "\n",
        "[ 343436.52753755  557663.12703671  577658.10171051 ...,  602293.83770182\n",
        "  550346.87793645  630643.26263446]\n"
       ]
      }
     ],
     "prompt_number": 7
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "e = prg.sigma_clip_vertical(queue, (1024//8, 1024), (1024//8, 1), \n",
      "                            datag.data, m.data, d.data, \n",
      "                            numpy.float32(-1), numpy.float32(3.0), numpy.float32(3.0), numpy.int32(5), local_mem)\n",
      "e.wait()\n",
      "#print(abs(data.mean(axis=0, dtype=\"float64\")-m.get()).max(),abs(data.std(axis=0, dtype=\"float64\")-d.get()).max())\n",
      "print(1e-6 * (e.profile.end - e.profile.start),\"ms\")\n",
      "print(m)\n",
      "print(d)\n",
      "print(\"error:\", abs(mean_n-m.get()).max(),abs(std_n-d.get()).max())\n"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "(0.514336, 'ms')\n",
        "[  92733.4375    159155.9375    158456.140625 ...,  159711.390625\n",
        "  159685.140625  191007.03125 ]\n",
        "[ 343436.5625  557663.125   577658.125  ...,  602293.8125  550346.875\n",
        "  630643.3125]\n",
        "('error:', 0.025501322001218796, 0.10890190454665571)\n"
       ]
      }
     ],
     "prompt_number": 8
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "print(\"Vertical speed-up\", 113/(1e-6 * (e.profile.end - e.profile.start)))"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "('Vertical speed-up', 219.7007403720525)\n"
       ]
      }
     ],
     "prompt_number": 9
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "#work on the horizontal direction\n",
      "%timeit sigma_clip(data, max_iter=5, axis=1)\n",
      "mean_n, std_n = sigma_clip(data, max_iter=5, axis=1)\n",
      "print(mean_n)\n",
      "print(std_n)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "10 loops, best of 3: 125 ms per loop\n",
        "[ 151751.73710174  136632.06422316  102139.70055282 ...,  109451.53472975\n",
        "  117544.86283701   98161.13061407]"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "\n",
        "[ 510740.43844003  471928.27273414  356969.13879386 ...,  419526.83165962\n",
        "  447362.00584674  354266.71424213]\n"
       ]
      }
     ],
     "prompt_number": 10
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "e = prg.sigma_clip_horizontal(queue, (1024, 1024//8), (1, 1024//8), \n",
      "                            datag.data, m.data, d.data, \n",
      "                            numpy.float32(-1), numpy.float32(3.0), numpy.float32(3.0), numpy.int32(5), local_mem)\n",
      "e.wait()\n",
      "print(1e-6 * (e.profile.end - e.profile.start),\"ms\")\n",
      "print(m)\n",
      "print(d)\n",
      "print(\"error:\", abs(mean_n-m.get()).max(),abs(std_n-d.get()).max())\n"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "(0.419968, 'ms')\n",
        "[ 151751.734375   136632.0625     102139.6953125 ...,  109451.53125\n",
        "  117544.875       98161.1328125]\n",
        "[ 510740.40625  471928.3125   356969.15625 ...,  419526.84375  447361.96875\n",
        "  354266.75   ]\n",
        "('error:', 0.028056760405888781, 0.10537082934752107)\n"
       ]
      }
     ],
     "prompt_number": 11
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "print(\"Horizontal speed-up\", 124/(1e-6 * (e.profile.end - e.profile.start)))"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "('Horizontal speed-up', 295.2605912831454)\n"
       ]
      }
     ],
     "prompt_number": 12
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "%%timeit \n",
      "prg.sigma_clip_horizontal(queue, (1024, 1024//8), (1, 1024//8), \n",
      "                            datag.data, m.data, d.data, \n",
      "                            numpy.float32(-1), numpy.float32(3.0), numpy.float32(3.0), numpy.int32(5), local_mem).wait()"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "1000 loops, best of 3: 453 \u00b5s per loop\n"
       ]
      }
     ],
     "prompt_number": 13
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "%%timeit \n",
      "prg.sigma_clip_vertical(queue, (1024//8, 1024), (1024//8, 1), \n",
      "                            datag.data, m.data, d.data, \n",
      "                            numpy.float32(-1), numpy.float32(3.0), numpy.float32(3.0), numpy.int32(5), local_mem).wait()"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "1000 loops, best of 3: 548 \u00b5s per loop\n"
       ]
      }
     ],
     "prompt_number": 14
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 14
    }
   ],
   "metadata": {}
  }
 ]
}