douglase/cudaFFTbenchmarking.ipynb

## cudaFFTbenchmarking.ipynb
{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "model name\t: Intel(R) Xeon(R) CPU @ 2.30GHz\n",
      "CPU(s):                8\n",
      "Thread(s) per core:    2\n",
      "Core(s) per socket:    4\n",
      "Socket(s):             1\n",
      "       product: GK210GL [Tesla K80]\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "WARNING: you should run this program as super-user.\n",
      "WARNING: output may be incomplete or inaccurate, you should run this program as super-user.\n"
     ]
    }
   ],
   "source": [
    "%%bash\n",
    "cat /proc/cpuinfo | grep 'model name' | uniq\n",
    " lscpu | egrep '^Thread|^Core|^Socket|^CPU\\('\n",
    " lshw -C display | grep 'product:'\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "#import accelerate\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "#import pyfftw\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.30.1\n",
      "blas_mkl_info:\n",
      "    libraries = ['mkl_intel_lp64', 'mkl_intel_thread', 'mkl_core', 'iomp5', 'pthread']\n",
      "    library_dirs = ['/home/edouglas/miniconda3/envs/accelerate/lib']\n",
      "    define_macros = [('SCIPY_MKL_H', None), ('HAVE_CBLAS', None)]\n",
      "    include_dirs = ['/home/edouglas/miniconda3/envs/accelerate/include']\n",
      "blas_opt_info:\n",
      "    libraries = ['mkl_intel_lp64', 'mkl_intel_thread', 'mkl_core', 'iomp5', 'pthread']\n",
      "    library_dirs = ['/home/edouglas/miniconda3/envs/accelerate/lib']\n",
      "    define_macros = [('SCIPY_MKL_H', None), ('HAVE_CBLAS', None)]\n",
      "    include_dirs = ['/home/edouglas/miniconda3/envs/accelerate/include']\n",
      "openblas_lapack_info:\n",
      "  NOT AVAILABLE\n",
      "lapack_mkl_info:\n",
      "    libraries = ['mkl_intel_lp64', 'mkl_intel_thread', 'mkl_core', 'iomp5', 'pthread']\n",
      "    library_dirs = ['/home/edouglas/miniconda3/envs/accelerate/lib']\n",
      "    define_macros = [('SCIPY_MKL_H', None), ('HAVE_CBLAS', None)]\n",
      "    include_dirs = ['/home/edouglas/miniconda3/envs/accelerate/include']\n",
      "lapack_opt_info:\n",
      "    libraries = ['mkl_intel_lp64', 'mkl_intel_thread', 'mkl_core', 'iomp5', 'pthread']\n",
      "    library_dirs = ['/home/edouglas/miniconda3/envs/accelerate/lib']\n",
      "    define_macros = [('SCIPY_MKL_H', None), ('HAVE_CBLAS', None)]\n",
      "    include_dirs = ['/home/edouglas/miniconda3/envs/accelerate/include']\n"
     ]
    }
   ],
   "source": [
    "import numba #http://jakevdp.github.io/blog/2013/06/15/numba-vs-cython-take-2/\n",
    "import numexpr as ne\n",
    "from  numba import vectorize\n",
    "\n",
    "#numba.vectorize?\n",
    "import math\n",
    "from numba import cuda\n",
    "\n",
    "from numba import jit\n",
    "\n",
    "#import mkl #for adjusting number of threads\n",
    "print( numba.__version__,)\n",
    "np.show_config()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true,
    "scrolled": true
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "array1024=np.ones([1024,1024],dtype=np.complex128)\n",
    "array1024[512,:]=0\n",
    "array4096=np.ones([4096,4096],dtype=np.complex128)\n",
    "#array8192=np.ones([8192,8192],dtype=np.complex128)\n",
    "out4096=array4096.copy()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.image.AxesImage at 0x7f88df18eda0>"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYYAAADxCAYAAADP0BAoAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAHcdJREFUeJzt3WuwZWV95/Hvz+ZiBBEJilxjk2nNtNQEpQuYeCkdIDQ9\nKVunapxmKkI0NS01QCRlKgPhhc4LqiiviRmEapURaxRCULTLakVgEq2pGu52tQIytIgRbJsAKgKx\noc/5zYu1tq5z2Puctc9eZ+21+vw+Vat673XZ+9m7+pz/Wc/z/P+PbBMRETHwomk3ICIiuiWBISIi\n5khgiIiIORIYIiJijgSGiIiYI4EhIiLmSGCIFUfSekkPSNop6eIhx39P0v+VtEfSX9S5VtJhkm6W\n9GD578vb+CwRyyGBIVYUSauAK4CzgLXA2ZLWzjvtSeDPgI+Oce3FwK221wC3ls8jeimBIVaak4Gd\nth+y/RxwHbCxeoLtx2zfCTw/xrUbgWvKx9cA71iuDxCx3PabdgMiFnPm2w7yE0/O1Dr37h177gV+\nVdm1xfaWyvOjgR9Xnj8CnFKzKQtde4TtXeXjnwJH1HzNiM5JYIjOe+LJGe646bha56468sFf2V63\nzE1akG1LSq2Z6K0Ehug8A7PMNvVyjwLHVp4fU+6b9Nrdko60vUvSkcBjE7c0YkoyxhCdZ8zznqm1\n1XAnsEbSakkHAJuArTWbstC1W4Fzy8fnAl+t/QEjOiZ3DNELTd0x2N4r6QLgJmAVcLXteyWdVx6/\nStKrgLuAQ4BZSRcBa20/Neza8qUvB66X9KfAj4B3NdLgiClQym5H173+9w/wt77+qlrnvuzoH989\n7TGGiL7LHUP0wiz5AyaiLQkM0XkGZhIYIlqTwBC9kDuGiPYkMETnGXg+Y2ERrUlgiM4zTldSRIsS\nGKL7DDOJCxGtSWCIzisynyOiLQkM0QNiBk27ERErRgJDdF4x+JzAENGWBIbovCKPIYEhoi0JDNEL\ns7ljiGhNAkN0Xu4YItqVwBCdZ8RMKsRHtCaBIXohXUkR7UlgiM4z4jmvmnYzIlaMBIbovCLBLV1J\nEW1JYIheyOBzRHsSGKLzbDHj3DFEtCU/bdELs6jWVoek9ZIekLRT0sVDjkvSJ8vjOyS9odz/Wknb\nK9tT5XrQSPqQpEcrxzY0+gVEtCh3DNF5xeBzM/9VJa0CrgDOAB4B7pS01fZ9ldPOAtaU2ynAlcAp\nth8ATqy8zqPAjZXrPmH7o400NGKKcscQnTcYfK6z1XAysNP2Q7afA64DNs47ZyPweRduAw6VdOS8\nc04DfmD7RxN+vIjOSWCIXpixam01HA38uPL8kXLfuOdsAq6dt+/Csuvpakkvr9OYiC5KYIjOG2Q+\n19mAwyXdVdk2N90eSQcAbwf+vrL7SuB4iq6mXcDHmn7fiLZkjCF6Ybb+rKTHba9b4PijwLGV58eU\n+8Y55yzgHtu7BzuqjyV9Gvha3QZHdE3uGKLziiJ6te8YFnMnsEbS6vIv/03A1nnnbAXOKWcnnQr8\nwvauyvGzmdeNJOnIwWyn8j1mlvhxI6YudwzReUY831BJDNt7JV0A3ASsAq62fa+k88rjVwHbgA3A\nTuBZ4D2D6yUdRDGj6X3zXvojwH8EfgjcBhwvae282U4RvZDAEJ1n02iCm+1tFL/8q/uuqjw2cP6I\na58BfnvIoSuAV9g+E0DSJRSzmxIYoncSGKIH6ievTdGwmUynVE8oB8I3A6xi1Ukv4ZD2Whcryq94\nhue8Z8k/NAkM0Xmm2TuGabG9BdgCcIgO8yk6bcotin3V7b51ousTGKIXerBQT53ZThG9kMAQnWfU\nh4V6fj3biSIgbAL+83SbFLE0CQzReQaeb6hW0nIZNdtpys2KWJJu/7RFAKBerMcwbLZTRB8lMETn\nmbEynyNiQgkM0Qt9uGOI2FckMETn2codQ0SLEhii84rB52ZKYkTE4hIYogey5nNEmxIYovOKweeM\nMUS0pfU/wxZbiD1imAbLbkfEIlq9Y6i5EHvEHD3JfI7YZ7T9J1adhdgjXmCWF9XaImJybY8xLFqa\nGOaWJ9YBB5y0/xGvbKd1MRV7n3ySmaefGXlLYMPzs/mlH9GWTg4+V8sTH3jcsT7qAxdNuUWxnH7y\nsb9e8HjRlZTAENGWtgNDShPHkiTzOaI9bQeGlCaOsWW6akS7Wr0/t70XGJQmvh+4PqWJY3FFV1Kd\nrdarLTJlWoVPlsd3SHpD5djDkr4rabukuyr7D5N0s6QHy39f3shHj5iC1jtubW+z/Rrbv2v7srbf\nP/pptlz3ebFtMZUp02cBa4GzJa2dd9pZwJpy2wxcOe/422yfaHtdZd/FwK221wC3ls8jeikjetF5\nxaykVbW2GupMmd4IfN6F24BDJR25yOtuBK4pH18DvKP+J4zolgSG6LxBgludDThc0l2VbfO8lxs2\nZfroMc4xcIuku+e99hG2d5WPfwocMcFHjpiqTk5XjZivTjdR6fF5XTxNe5PtRyW9ErhZ0vdtf7t6\ngm1L8jK2IWJZ5Y4hOm8wK6nmHcNi6kyZHnmO7cG/jwE3UnRNAewedDeV/z423qeM6I4EhuiFBmcl\n/XrKtKQDKKZMb513zlbgnHJ20qnAL2zvknSQpJcCSDoI+EPge5Vrzi0fnwt8dbJPHDE96UqKzrPF\n3oYyn23vlTSYMr0KuNr2vZLOK49fBWwDNgA7gWeB95SXHwHcKAmKn50v2v5Geexy4HpJfwr8CHhX\nIw2OmIIEhuiFJhPcbG+j+OVf3XdV5ekxFFNZ91DcVb8RuAv4OUUX0auBh4FPVa7ZDPwOMAN8xPaT\njTU4omXpSorOa3iMoY69wAdsrwVOBc4vcx2G5iqUxzYBrwPWA58q8yUieimBIXqhzcBge5fte8rH\nv6TI0j+a0bkKG4HrbO+x/UOKLqiTieipdCVF501zoR5JrwZeD9zO6FyFo4HbKpcNy42YU07+xbxk\neRoc0YAEhuiFMfIYGiPpYOBLwEW2nyoHnYGl5SpUy8kfosOS5xCdlcAQnWfD3pYX6pG0P0VQ+ILt\nL5e7d0s6spy6Ws1VSDn52KdkjCF6oc0xBhW3Bp8F7rf98cqhUbkKW4FNkg4sS8qvAe5opDERU5A7\nhui8KYwxvBF4N/BdSdvLfX/FiFyFMg/ieuA+ihlN59ueabPBEU1KYIhecIuBwfb/gZGDGqeNuOYy\nIGXkY5+w5K4kScdK+gdJ90m6V9L7y/0jFyyRdEm5+MkDks5s4gPEytDUegwRsbhJxhiSBBStsFtP\ncItY0ZYcGJIEFO0RM7MvqrVFxOQa+UkaIwlosQVSBq+3ebDQyszTzzTRxOg5W7W2iJjcxIFhfhJQ\n9ZhtU5S6GYvtLbbX2V636uCDJm1i9NwUaiVFrGgTBYaFkoDK40kCism5GGeos0XE5CaZlZQkoGhN\nZiVFtGeSPIYkAUUrXA4+R0Q7lhwYkgQUbUo3UUR7kvkcvZAZRxHtyf15dF4xsNzcdFVJ68vs+52S\nLh5yXJI+WR7fIekN5f6h2f7lsQ9JelTS9nLb0NgXENGy3DFELzRYOXUVcAVwBkUuzZ2Sttq+r3La\nWRSTI9YApwBXlv8Osv3vkfRS4G5JN1eu/YTtjzbS0Igpyh1D9EKD01VPBnbafsj2c8B1FFn5VRuB\nz7twG3DoYB2GEdn+EfuUBIboPCNmZ19UawMOH2TNl9vmeS9XJwN/0XPmZfsPXFh2PV1dLR4Z0TcJ\nDNELrrkBjw+y5sttS9NtGZHtfyVwPHAisAv4WNPvG9GWjDFE97nRWUl1MvBHnjMi2x/buwePJX0a\n+FpTDY5oW+4Yoh/GuGVYxJ3AGkmrJR1AUQp+67xztgLnlLOTTgV+Ua7zPCrbf1D+ZeCdwPfG+4AR\n3ZE7huiFpu4YbO+VdAFwE7AKuLrMyj+vPH4VsA3YQFEa/lngPeXlQ7P9bW8DPizpRIrw9DDwvkYa\nHDEFCQzReQZmZ5tLcCt/kW+bt++qymMD5w+5bmS2v+13N9bAiClLYIjuM5DM54jWZIwhemEaZbcl\nrZL0HUlfK59nPfNYERIYoh+aG3wex/spktgGsp55rAgJDNED9eokNVloT9IxwL8HPlPZnfXMY0VI\nYIh+aP+O4a+BvwRmK/smXs88og+aWPM5/bCxvAyeVa2tCZL+CHjM9t0jm7SE9cwlbR6U6niePZM2\nM2LZNHHHkH7YaIFqbo14I/B2SQ9TFNn7d5L+FxOuZ257y6BUx/4c2FRbIxo3UWBIP2y0psWuJNuX\n2D7G9qsp/pj537b/mKxnHivEpHkMg37Yl1b2LdQPe1vlvJH9sGVFzM0Aq16eIpXBcsw4WoqsZx4r\nwpIDQ7UfVtJbh51j25LG/pEuK2JuATjwuGO78SshpmeKCW62/xH4x/LxE2Q981gBJrljGPTDbgBe\nDBxS7Ycti46N3Q8bMUzTyWsRMdqSxxjSDxutmlW9LSImthy1ktIPG40bv0MyIpaqkcCQfthYVstT\n7iIiRkh11egBpbpqRIsSGKIfcscQ0ZoEhuiH2cVPiYhmJDBE92WhnohWpbpq9IJcb6v1WtL6spDj\nTkkXDzkuSZ8sj++Q9IbFrl2oeGRE3yQwRD80VCupLNx4BXAWsBY4uyzwWHUWRZ7NGorSLFfWuHZo\n8ciIPkpgiJXmZGCn7YdsP0dRPXXjvHM2Ap934Tbg0DKLf6FrRxWPjOidBIbohTG6kg4frHlQbpvn\nvVSdRXVGnbPQtaOKR0b0Tgafo/vMOOUuHre9bhlbs6ilFo+M6IoEhuiH5n7N1inmOOqc/Re4dlTx\nyIjeSVdS9EKDs5LuBNZIWi3pAIoCkFvnnbMVOKecnXQq8Iuym2iha0cVj4zondwxRD80tzrbXkkX\nADcBq4CrywKP55XHrwK2ARsoVhl8FnjPQteWLz20eGREHyUwRD802GNvexvFL//qvqsqjw2cX/fa\ncv/I4pERfZPAEJ03TvJaRExuojEGSYdKukHS9yXdL+nfLpQBKumSMmP0AUlnTt78WDGyUE9EayYd\nfP4b4Bu2fw/4feB+RmSAlhmim4DXAeuBT5WZpBGLarIkRkQsbMmBQdLLgLcAnwWw/ZztnzM6A3Qj\ncJ3tPbZ/SDGwd/JS3z9WmIZKYkTE4ia5Y1gN/DPwPyV9R9JnJB3E6AzQOhmnES9U826hyTuGdJPG\nSjZJYNgPeANwpe3XA88wr3BYObtj7B9XSZsHJQ1mnn5mgibGPqP9O4Z0k8aKNUlgeAR4xPbt5fMb\nKALF7jLzk3kZoHUyTgGwvcX2OtvrVh180ARNjH2FZuttjbxXukljhVtyYLD9U+DHkl5b7joNuI/R\nGaBbgU2SDpS0mqKk8R1Lff+IZbQs3aTVO+Hn2bOMzY+YzKR5DBcCXyjLAzxEkSH6IoZkgJbZpddT\nBI+9wPm2ZyZ8/1gp2h1YHnSTXmj7dkl/w5Bu0nEL5dneAmwBOESHZag8OmuiwGB7OzCskuXQDFDb\nlwGXTfKesQK1PxV1WDfpxYwulFe7mzSiD1JEL/qhxcHndJPGSpeSGNEP7Xe8pJs0VqwEhug80dyM\no7rSTRorWQJDdF/KXUS0KoEh+iGBIaI1CQzRDwkMEa1JYIheSFdSRHsSGKIfEhgiWpM8hug+t1Mr\naaHqqfPOW19WUd0p6eLK/o+U1Vh3SLpR0qHl/ldL+hdJ28vtqmGvG9EVCQzRD+0kuA2tnlpVVk29\nAjgLWAucXVZXBbgZOMH2vwH+H3BJ5dIf2D6x3M6buKURyyiBIXqhpfUYRlVPrToZ2Gn7IdvPAdeV\n12H7m7b3lufdRlEaI6J3EhiiH+rfMRw+qGBabpvHeJdR1VOr6i449V7g65Xnq8tupG9JevMYbYpo\nXQafo/vG6yZ63PawjGUAJN0CvGrIoUvnvOUSqqdW3uNSitIYXyh37QKOs/2EpJOAryzldSPaksAQ\nnSeam65q+/SR7yONqp5atWAlVUl/AvwRcFq5giG290CxAIPtuyX9gHQzRYelKyl6oaUxhlHVU6vu\nBNZIWl0W2NtUXoek9cBfAm+3/eyv2y69YrDUp6TjKaqvRnTWRIFB0p9LulfS9yRdK+nFWTA9lkU7\ns5IuB86Q9CBwevkcSUdJ2gZQDi5fANxEsQ709bbvLa//H8BLgZvnTUt9C7BD0naKtR0yKyk6bcld\nSZKOBv4MWGv7X8qyw5sopvDdavvyco73xcB/m7dg+lHALZJek/LEUUsLCW62n2BI9VTbPwE2VJ5v\nA7YNOe9fjXjdLwFfqu47RIdN2tyIZTNpV9J+wG9J2g94CfATsmB6NK1mN1LKZkQ0Y8mBwfajwEeB\nf6KYdfEL299kwgXTYe6i6TNPP7PUJsa+pMUV3CJWuiUHhnLsYCOwmqJr6CBJf1w9p5yVMfaPq+0t\nttfZXrfq4IOW2sTYh7RREiMiCpN0JZ0O/ND2P9t+Hvgy8AeUC6YDZMH0aEq6kiLaM0lg+CfgVEkv\nkSSKQbv7yYLp0bS63UgJDBGNWPKsJNu3S7oBuIciy/M7wBbgYLJgejQtv/QjWjNR5rPtDwIfnLd7\nD1kwPRrUZOZzRCwumc/RC5p1ra2x90vyZqxgCQzRfS2PMVSSN9fZPgFYRZGcOXS9hnnJm+uBTw1K\nYET0UQJD9MIUZiUleTNWrASG6IcW7xiWK3mzmrj5fFFsNaKTEhiiF9q8Y1iu5M1q4ub+HNhMYyOW\nQQJD9EO7eQxJ3owVLYEhus+tl8RI8masaFnBLTqv7TyGJG/GSpfAEP3gdjPckrwZK1kCQ/RCMp8j\n2pPAEN2XAnkRrcrgc/RCG4PPC5W8mHfe+rL0xc5y+drB/g9JerRc73m7pA2VYymZEb2RwBC90NKs\npKElL+a0oyh1cQVwFsX65meXJTEGPmH7xHLbVl7zgpIZE7c0YhklMET3mWLwuc42mVElL6pOBnba\nfsj2c8B15XWLve78khkRnZXAEL0wRubz4YOyE+W2eYy3GVXyomqx8hcXStoh6epKV9SwayI6a9HA\nUP4Hf0zS9yr7xi4/LOkkSd8tj32yTByKqKd+5vPjg7IT5bal+jKSbilLac/f5vzVv8T1yq8EjgdO\npKix9LFxP2ZEF9S5Y/gcRb9o1VLKD18J/BeKrNA1Q14zYqhBglsTtZJsn277hCHbVxld8qJqZPkL\n27ttz9ieBT7NbyqsDrsmorMWDQy2vw08OW/3WOWHyx+yQ2zfVv4l9nmG999GvJDrLdLTwEI9o0pe\nVN0JrJG0WtIBFH8IbYVfB5OBdwKDu+xhJTMiOmupYwzjlh8+mrn9qkPLEg9UyxPPPP3MEpsY+5R2\niuhdDpwh6UGKQnqXA0g6StI2ANt7gQuAmyjqJ11v+97y+g+X3aU7gLcBf15ecy8wKJnxDeD8iVsa\nsYwmTnCzbanZvNSyX3gLwIHHHZvUpmgl89n2EwwpeWH7J8CGyvNtwLYh5717gdeeUzLjEB02aXMj\nls1S7xjGLT/8KHP7VVOWOOozMOt6W0RMbKl3DIO+2Mt5YfnhL0r6OMUCJ2uAO2zPSHpK0qnA7cA5\nwN/WeaPnfvzI0w9f9BcPLLGdfXM48Pi0G9Giwef9nUXPzO/8iNYsGhgkXQu8lWJ++CMUFScvZ/zy\nw/+VYobTbwFfL7c6HrC9ru4H6jNJd62Uzwrjfd4U0Ytoz6KBwfbZIw6NVX7Y9l3ACWO1LqLUwIyj\niKgp1VWj+1JdNaJVfQgMWxY/ZZ+xkj4r1Py8RYJbIkNEWzofGOaXNNiXraTPCmN+3ubWc46IRXQ+\nMERA7hgi2pTAEN2XMYaIVnW27PaoVbL6TNLDZcmE7ZLuKveNXam2q5arEi+0VispIuhoYKixSlaf\nva1c3Wswf38plWq76nMsTyXeZVmoJyXlI4brZGBgaatk9dVYlWqn0L7alq0Sr5dtac/PkZLyES/Q\n1cCw2CpZfWXgFkl3V1YWG7dSbd80UYl3We4YUlI+YrgMPrfrTbYflfRK4GZJ368eXI5KtV0y0edr\n71tZKJDdVjlvEMiep2ZJ+fKPgc0AL+YlDTY5olldvWMYuUpWn9kerPT1GHAjRdfQuJVq+6aJSrxo\ndrbW1qQlLu+50OttGSw5uj8HNvWyEY3ramAYuUpWX0k6SNJLB4+BP6RY4WvUqmHDVv26o91WN2Ks\nz1f+tf6UpFPLQdxzirLbNbfJpaR8rHid7EqyvVfSYJWsVcDVlVWy+uoI4MZywsp+wBdtf0PSnYxf\nqbaTlqsSr3CbCW6tlJT/JT97+hbf0MVy8l0s/Z421VNt0+Kl7BcgJ6M0Ou5lBx3lU//15sVPBL55\n93+/u3Yp70ogA3ZTBLKvUCzDeRxlILP9ZHn+pcB7KQLZRba/Xu5fx9yS8hd6kR+srpZY72K70qZ6\nmmxTJ+8YIl5gGf6AmV9SXtJhwN9RBIWHKYLCzyrnXwZcJmk98DeS/hb4jO3LgRMk/R3wJuA7kg4F\nfm77REmvplgfenCHUB3Ejuicro4xRPxGe2MMQ3MYqhZKvrT9n8rkxROBLwFfrlz6g8Ex2+dN3NKI\nZZTAEL3Q0qykUTkMVYsmX5aD5u8Crl3gvbpaSbeL7Uqb6mmsTQkM0QM1k9sm724alcNQVSfx8M3A\nbtsPVvatLmtkfUvSm7taYr2L7Uqb6mmyTRljiO4z4/zSP3xQoLC0pfoDI+kW4FVDrrt0zltOlmx4\nNnPvFnYBx9l+QtJJwFckvc72U0t8/YhllcAQ/VC/l+jxhWZm2D591DFJuyUdaXvXvByGqgUTDyXt\nB/wH4KTKe+4B9pSP75b0A+A1QDWARXRGupKiF2TX2iY0KhmvarHky9OB79v+dZkMSa8YFNyTdDxw\nAnCtWi4pL+lYSf8g6T5J90p6f7n/Q5IeLbu6tkvaULlm2Uu/q2Pl6CW9tvJdbC/zVC5q+3vSFKv/\nJjBEP7QzxnA5cIakByl+wV8OIOkoSduKZngvMEi+vB+4fl7y5SZeOOj8FmCHpO3ADRT1lc6k/ZLy\ne4EP2F4LnAqcX3nvT1RmTW2D1ku/d6Ycve0HKrPLTgKepShhA+1+T59jStV/Exii+2yYma23TfQ2\nfsL2abbX2D59kNhm+ye2N1TO22b7NbZ/t8xtqL7Gn9i+at6+L9l+XfmL5nxgxzRKytveZfue8vEv\nKQLbQhV7p1n6vSvl6E+jmGr8owXOWZY2TbP6bwJD9EM7dwxt6EQ59TLp7vUUZTwALpS0o+y+GHRP\ntNXWLpejn38HOM3vCZopY79o2xIYoh/2ncAwdZIOpkjAu6icGXUlcDxwIsUMqo+13KQ3lXdTZ1F0\nb72lerDpKrd1lWNIbwf+vtw17e9pjuX8XhIYovsMzLre1n1TLacuaX+KoPAF218GsL3b9oztWeDT\n/KYbpJW2drgc/VnAPbZ3l+2b6vdUaqX6bwJD9IDBs/W27ptaSflyNspngfttf7yy/8jKae+kKAcP\nLZR+V7fL0c/JR5nm91QxeRn74bPt5kgeQ3SfmXhguSumXFL+jcC7ge+WM6QA/opiZtSJFN/0w8D7\nyra2Ufq9k+XoyyB1BuV3Ufpwm9+TlqmMfbkt/N4pux1d97IDjvAfHLGp1rnfeOSTtctuR8RwuWOI\nfsgfMBGtSWCIHsiMo4g2JTBE9xmYvKR2RNSUwBD9kDuGiNYkMEQPeJ+ZlRTRBwkM0X0G9yNHIWKf\nkMAQ/dCPrOaIfUICQ/RDxhgiWpPAEN1nZ1ZSRIsSGKIfcscQ0ZoEhugB45nGy+FExAgJDNF9g7Lb\nEdGKBIboh0xXjWhN1mOIzjPgWdfaJiHpMEk3S3qw/PflI867WtJjkr5X93pJl0jaKekBSWdO1NCI\nZZbAEN3n1hbquRi41fYa4Nby+TCfA9bXvV7SWooFeV5XXvcpSasmbWzEcklgiF7wzEytbUIbgWvK\nx9cA7xjaFvvbwJNjXL8RuM72Hts/BHbym2UhIzonYwzReb/kZzfd4hsOr3n6iyXdVXm+xfaWmtce\nUS6FCPBTitXFxjHq+qOB2yrnPVLui+ikBIboPNvDum2WRNItwKuGHLp03nta0pIHLSa9PmKaEhhi\nRbF9+qhjknZLOtL2rnLh98fGfPlR1z8KHFs575hyX0QnZYwh4je2AueWj88FvtrQ9VuBTZIOlLQa\nWAPcMWFbI5aNnFIDEQBI+m3geuA44EfAu2w/Keko4DO2N5TnXQu8FTgc2A180PZnR11fXnMp8F5g\nL3CR7a+3+uEixpDAEBERc6QrKSIi5khgiIiIORIYIiJijgSGiIiYI4EhIiLmSGCIiIg5EhgiImKO\n/w+e6fOqJqosWwAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7f890c815a90>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "array1024=np.ones([1024,1024],dtype=np.complex128)\n",
    "array1024[512,:]=0\n",
    "out1024=np.zeros(array1024.shape,dtype=np.complex128)\n",
    "import accelerate.cuda\n",
    "accelerate.cuda.fft.fft(array1024,out1024)\n",
    "plt.subplot(121)\n",
    "plt.imshow(((out1024.real-np.fft.fft2(array1024).real)))\n",
    "plt.colorbar()\n",
    "plt.subplot(122)\n",
    "plt.imshow(np.fft.ifft2(array1024).real)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "10 ms ± 192 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)\n",
      "77.7 ms ± 1.36 ms per loop (mean ± std. dev. of 7 runs, 10 loops each)\n"
     ]
    }
   ],
   "source": [
    "\n",
    "\n",
    "%timeit accelerate.cuda.fft.ifft_inplace(array1024)\n",
    "%timeit np.fft.ifft2(array1024)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(array([ 16777216.+0.j,         0.+0.j,         0.+0.j, ...,         0.+0.j,\n",
       "                0.+0.j,         0.+0.j]),\n",
       " array([ 16777216.+0.j,         0.+0.j,         0.+0.j, ...,         0.+0.j,\n",
       "                0.+0.j,         0.+0.j]),\n",
       " array([ 16777216.+0.j,         0.+0.j,         0.+0.j, ...,         0.+0.j,\n",
       "                0.+0.j,         0.+0.j]),\n",
       " array([ 1.+0.j,  0.+0.j,  0.+0.j, ...,  0.+0.j,  0.+0.j,  0.+0.j]))"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "plan=accelerate.cuda.fft.FFTPlan(array4096.shape,np.complex128,np.complex128)\n",
    "plan.forward(array4096,out=out4096)[0],plan.inverse(array4096,out=out4096)[0],np.fft.fft2(array4096)[0],np.fft.ifft2(array4096)[0]\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Compare cuda iFFT to numpy and pyFFTW\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "136 ms ± 7.64 ms per loop (mean ± std. dev. of 7 runs, 10 loops each)\n",
      "162 ms ± 8.57 ms per loop (mean ± std. dev. of 7 runs, 10 loops each)\n",
      "1.11 s ± 31.8 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)\n"
     ]
    }
   ],
   "source": [
    "%timeit plan.forward(array4096,out=out4096)\n",
    "%timeit accelerate.cuda.fft.ifft(array4096,out=out4096)\n",
    "#%timeit pyfftw.interfaces.numpy_fft.ifft2(array4096,threads=16)\n",
    "#%timeit pyfftw.interfaces.numpy_fft.ifft2(array4096,threads=32)\n",
    "\n",
    "\n",
    "%timeit np.fft.fft2(array4096)\n",
    "#%timeit pyfftw.interfaces.numpy_fft.ifft2(array4096)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "#%timeit pyfftw.interfaces.numpy_fft.ifft2(array4096,planner_effort=\"FFTW_EXHAUSTIVE\")\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "mkl.set_num_threads(1)\n",
    "%timeit np.fft.fft2(array4096)\n",
    "mkl.set_num_threads(16)\n",
    "%timeit np.fft.fft2(array4096)\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "#only run if memory available\n",
    "#%timeit accelerate.cuda.fft.ifft_inplace(array8192)\n",
    "#%timeit np.fft.fft2(array8192)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "## the challenge:\n",
    "\n",
    "        rsqd = (x ** 2 + y ** 2) * u.m ** 2\t\t\n",
    "        k = 2 * np.pi / wave.wavelength\n",
    "        lens_phasor = np.exp(1.j * k * rsqd / (2.0 * self.z))\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "n=2048\n",
    "y=np.ones([n,n])\n",
    "x=2*np.ones([n,n])\n",
    "z=9#2*np.ones([n,n])#*1.j\n",
    "#x.dtype,y.dtype,z.dtype"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[ 0.55555556,  0.55555556,  0.55555556, ...,  0.55555556,\n",
       "         0.55555556,  0.55555556],\n",
       "       [ 0.55555556,  0.55555556,  0.55555556, ...,  0.55555556,\n",
       "         0.55555556,  0.55555556],\n",
       "       [ 0.55555556,  0.55555556,  0.55555556, ...,  0.55555556,\n",
       "         0.55555556,  0.55555556],\n",
       "       ..., \n",
       "       [ 0.55555556,  0.55555556,  0.55555556, ...,  0.55555556,\n",
       "         0.55555556,  0.55555556],\n",
       "       [ 0.55555556,  0.55555556,  0.55555556, ...,  0.55555556,\n",
       "         0.55555556,  0.55555556],\n",
       "       [ 0.55555556,  0.55555556,  0.55555556, ...,  0.55555556,\n",
       "         0.55555556,  0.55555556]])"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "#from numba import vectorize, float64\n",
    "import math\n",
    "import cmath\n",
    "import numba\n",
    "\n",
    "dtype=numba.complex128\n",
    "@numba.vectorize([dtype(dtype, dtype,dtype)])\n",
    "def f_numba(x, y,z):\n",
    "    return (x**2 + y**2)/z\n",
    "\n",
    "@numba.vectorize([numba.float64(numba.float64, numba.float64,numba.float64)],\n",
    "                 target='cuda'\n",
    "                )\n",
    "def f_cuda(x, y,z):\n",
    "    return (x**2 + y**2)/z\n",
    "f_cuda(x,y,z) "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([ 1.+2.j])"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "\n",
    "@numba.vectorize(['complex128(float64, float64)'],target='cuda')\n",
    "def foo(real, imag):\n",
    "    a = complex(real, imag)\n",
    "    return a * a + complex(1, 0)         # complex power is not implemented yet\n",
    "foo(1,1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Time sum of squares"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "9.36 ms ± 1.18 ms per loop (mean ± std. dev. of 7 runs, 100 loops each)\n",
      "53.3 ms ± 325 µs per loop (mean ± std. dev. of 7 runs, 10 loops each)\n",
      "54.9 ms ± 286 µs per loop (mean ± std. dev. of 7 runs, 10 loops each)\n",
      "55.9 ms ± 1.88 ms per loop (mean ± std. dev. of 7 runs, 10 loops each)\n"
     ]
    }
   ],
   "source": [
    "%timeit ne.evaluate(\"(x**2 + y**2)/z\")\n",
    "%timeit   (x**2 + y**2)/z\n",
    "%timeit  f_numba(x,y,z)\n",
    "%timeit  f_cuda(x,y,z) \n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Time exponential sum of squares"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "11 ms ± 1.26 ms per loop (mean ± std. dev. of 7 runs, 100 loops each)\n",
      "180 ms ± 1.35 ms per loop (mean ± std. dev. of 7 runs, 10 loops each)\n",
      "142 ms ± 1.97 ms per loop (mean ± std. dev. of 7 runs, 10 loops each)\n"
     ]
    }
   ],
   "source": [
    "@numba.vectorize([dtype(dtype,dtype,dtype)],)\n",
    "def fexp(x, y,z):\n",
    "    return cmath.exp((x**2 + y**2)/z)\n",
    "\n",
    "%timeit ne.evaluate(\"exp((x**2 + y**2)/z)\")\n",
    "%timeit fexp(x,y,z) #numba optimized\n",
    "%timeit np.exp((x**2 + y**2)/z)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "#make complex \n",
    "x=x+1j    \n",
    "\n",
    "y=y+1j    \n",
    "   "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "21.1 ms ± 226 µs per loop (mean ± std. dev. of 7 runs, 10 loops each)\n",
      "252 ms ± 1.74 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)\n",
      "352 ms ± 4.16 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)\n"
     ]
    }
   ],
   "source": [
    "%timeit ne.evaluate(\"exp((x**2 + y**2)/z)\")\n",
    "%timeit fexp(x,y,z) #numba optimized\n",
    "%timeit np.exp((x**2 + y**2)/z)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "23.1 ms ± 473 µs per loop (mean ± std. dev. of 7 runs, 10 loops each)\n",
      "89.9 ms ± 2.59 ms per loop (mean ± std. dev. of 7 runs, 10 loops each)\n",
      "230 ms ± 1.32 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)\n",
      "222 ms ± 847 µs per loop (mean ± std. dev. of 7 runs, 1 loop each)\n",
      "764 ms ± 5.4 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)\n"
     ]
    }
   ],
   "source": [
    "#simple function to make complex values CUDA friendly:\n",
    "@vectorize([numba.complex128(numba.float64,numba.float64)],\n",
    "          target='cuda')\n",
    "def f(real,imag):\n",
    "    #work around to get vectorize to work on complex arrays, \n",
    "    #see: https://groups.google.com/a/continuum.io/forum/#!topic/numba-users/F3fYUiEXmn0\n",
    "        c=complex(real,imag)\n",
    "        return cmath.exp(c)\n",
    "\n",
    "\n",
    "a=ne.evaluate(\"(x**2 + y**2)/z\")\n",
    "\n",
    "%timeit ne.evaluate(\"exp(a)\")\n",
    "%timeit f(np.ascontiguousarray(a.real),np.ascontiguousarray(a.imag)) #sometimes get non-contigous array error\n",
    "%timeit ne.evaluate(\"cos(a)+1j*sin(a)\")\n",
    "\n",
    "%timeit np.exp(a)\n",
    "%timeit np.cos(a)+1j*np.sin(a)\n",
    "\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "163 ms ± 137 µs per loop (mean ± std. dev. of 7 runs, 10 loops each)\n"
     ]
    }
   ],
   "source": [
    "#bring it all together:\n",
    "@vectorize([numba.complex128(numba.float64,numba.float64,\n",
    "                            numba.float64,numba.float64,\n",
    "                            numba.float64,numba.float64)],\n",
    "          target=\"parallel\")\n",
    "\n",
    "def f_xyz_def(xr,xi,\n",
    "         yr,yi,\n",
    "         zr,zi):\n",
    "        x=complex(xr,xi)\n",
    "\n",
    "        y=complex(yr,yi)\n",
    "\n",
    "        z=complex(zr,zi)\n",
    "        return cmath.cos((x**2+y**2)/z) + 1j*cmath.sin((x**2+y**2)/z)\n",
    "        #return cmath.exp((x**2+y**2)/z)\n",
    "%timeit f_xyz_def(x.real,x.imag,y.real,y.imag,z.real,z.imag)\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "177 ms ± 24.2 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)\n"
     ]
    }
   ],
   "source": [
    "#bring it all together:\n",
    "@vectorize([numba.complex128(numba.float64,numba.float64,\n",
    "                            numba.float64,numba.float64,\n",
    "                            numba.float64,numba.float64)],\n",
    "          target='cuda')\n",
    "\n",
    "def f_xyz_cuda(xr,xi,\n",
    "         yr,yi,\n",
    "         zr,zi):\n",
    "        x=complex(xr,xi)\n",
    "\n",
    "        y=complex(yr,yi)\n",
    "\n",
    "        z=complex(zr,zi)\n",
    "        return cmath.cos((x**2+y**2)/z) + 1j*cmath.sin((x**2+y**2)/z)\n",
    "        #return cmath.exp((x**2+y**2)/z)\n",
    "%timeit f_xyz_cuda(x.real,x.imag,y.real,y.imag,z.real,z.imag)\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "32.7 ms ± 108 µs per loop (mean ± std. dev. of 7 runs, 10 loops each)\n",
      "12.9 ms ± 39.8 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)\n",
      "27.4 ms ± 258 µs per loop (mean ± std. dev. of 7 runs, 10 loops each)\n"
     ]
    }
   ],
   "source": [
    "#ne is still the fastest, and it's marginally better to do it all in one go:\n",
    "\n",
    "%timeit ne.evaluate(\"exp((x**2 + y**2)/z)\")\n",
    "%timeit a= ne.evaluate(\"(x**2 + y**2)\")\n",
    "%timeit ne.evaluate(\"exp(a/z)\")\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "import math\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "import poppy\n",
    "import astropy.units as u\n",
    "wf=poppy.FresnelWavefront(2.4*u.m,npix=512,oversample=2)\n",
    "wf *= poppy.CircularAperture(radius=2.4/2)\n",
    "array=wf.intensity\n",
    "tmp=array.copy()\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "@cuda.jit\n",
    "def increment_by_one(an_array):\n",
    "    pos = cuda.grid(1)\n",
    "    if pos < an_array.size:\n",
    "        an_array[pos] += 1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "#import accelerate\n",
    "import numpy as np\n",
    "#import matplotlib.pyplot as plt\n",
    "#import pyfftw\n",
    "#%matplotlib inline\n",
    "import numba #http://jakevdp.github.io/blog/2013/06/15/numba-vs-cython-take-2/\n",
    "import numexpr as ne\n",
    "from  numba import vectorize\n",
    "\n",
    "#numba.vectorize?\n",
    "import math\n",
    "\n",
    "#%timeit -n 5 cufftShift_2D_kernel[4,1](flatarr,n)\n",
    "#cufftShift_2D_kernel(array4096,4096)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from numba import cuda\n",
    "\n",
    "@cuda.jit()\n",
    "def cufftShift_2D_kernel(data, N):\n",
    "    '''\n",
    "    adopted CUDA FFT shift code from:\n",
    "    https://github.com/marwan-abdellah/cufftShift\n",
    "    (GNU Lesser Public License)\n",
    "    \n",
    "    '''\n",
    "    \n",
    "    #// 2D Slice & 1D Line\n",
    "    sLine = N\n",
    "    sSlice = N * N\n",
    "    #// Transformations Equations\n",
    "    sEq1 = int((sSlice + sLine) / 2)\n",
    "    sEq2 = int((sSlice - sLine) / 2)\n",
    "    x, y = cuda.grid(2)\n",
    "    #// Thread Index Converted into 1D Index\n",
    "    index = (y * N) + x\n",
    "    \n",
    "    if (x < N / 2):\n",
    "        if (y < N / 2):\n",
    "            #// First Quad\n",
    "            temp = data[index]\n",
    "            data[index] = data[index + sEq1]\n",
    "            #// Third Quad\n",
    "            data[index + sEq1] = temp\n",
    "    else:\n",
    "        if (y < N / 2):\n",
    "            #// Second Quad\n",
    "            temp=data[index]\n",
    "            data[index] = data[index + sEq2]\n",
    "            data[index + sEq2] = temp\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "CPU times: user 432 ms, sys: 8 ms, total: 440 ms\n",
      "Wall time: 441 ms\n"
     ]
    }
   ],
   "source": [
    "#from cudashift import cufftShift_2D_kernel\n",
    "\n",
    "    \n",
    "n=1024*8\n",
    "array=np.ones([n,n])+1j#,dtype=np.complex128)\n",
    "array[0,0]=0\n",
    "array[-1,-1]=2\n",
    "array[-1,0]=3\n",
    "array[0,-1]=4\n",
    "array=array.ravel()\n",
    "\n",
    "#calculate correct number of blocks\n",
    "#https://stackoverflow.com/questions/2392250/understanding-cuda-grid-dimensions-block-dimensions-and-threads-organization-s\n",
    "blockdim = (32, 32) # threads per block\n",
    "numBlocks = (int(n/blockdim[0]),int(n/blockdim[1]))\n",
    "%time cufftShift_2D_kernel[numBlocks, blockdim](array,n)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/home/edouglas/miniconda3/envs/accelerate/lib/python3.6/site-packages/ipykernel_launcher.py:2: VisibleDeprecationWarning: using a non-integer number instead of an integer will result in an error in the future\n",
      "  \n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "<matplotlib.colorbar.Colorbar at 0x7f88ddb01b70>"
      ]
     },
     "execution_count": 27,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAUAAAAD8CAYAAAAG730QAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAGSBJREFUeJzt3X+sHeWd3/H3J8aELaElWS9g/EOwklUJophQy6EEbcxu\nkjUuW2cllBptASGiW1Cokirb1s1KpNu/0t02rVJI3FtiBdQAQg0GixqoSVlBGkFsXMdgfmzueomw\n14nXJjFxICG+/vSPmZueHM6PuZw5955z5/OSRndmnmee8xyO/OV55plnHtkmIqKJ3jXfFYiImC8J\ngBHRWAmAEdFYCYAR0VgJgBHRWAmAEdFYAwVASe+TtFPS98u/7+2S7xVJz0naK2n3bK+PiJC0SNL/\nlfRwhzRJ+rKkKUn7JF1apcxBW4CbgW/ZXgV8qzzu5krbl9he8w6vj4hm+wzwYpe0q4BV5TYBfLVK\ngYMGwI3AXeX+XcAn5vj6iGgAScuBfwTc2SXLRuBuF54Gzpa0tF+5pw1Yr3NtHy73fwic2yWfgccl\nTQP/1fbkLK9H0gRFZEenn/4PFp97zoBVj7l01plvzncVYhZOHD7Bz3/ycw1Sxu9feaaPvTZdKe+z\n+36xH/h5y6nJljgB8J+BfwWc1aWIZcCrLccHy3OHO2cv9A2Akh4HzuuQ9CetB7Ytqdu8uitsH5J0\nDrBT0ku2n5zF9ZT/MSYB3r1yhc//3Gf7VT1GyO/8w/3zXYWYhf95w/aByzj22jTffWxlpbyLln7/\n5223x35F0tXAEdvPSlo3cMVa9A2Atj/aLU3SjyQttX24bG4e6VLGofLvEUnbgLXAk0Cl6yNi/Bg4\nxak6ivow8I8lbQDOAP6upP9u+5+25DkErGg5Xl6e62nQe4DbgRvK/RuAh9ozSDpT0lkz+8DHgeer\nXh8R48mYX3q60tazHPvf2F5u+wJgE/C/24IfFLHk+nI0+DLgeMvtta4GvQf4ReB+STcBPwA+CSDp\nfOBO2xso7uttkzTzeffYfrTX9RGxMNTUAuxI0s0AtrcAO4ANwBTwBnBjlTIGCoC2jwG/1+H835SV\nwfYBYPVsro+I8WfMdM2v27P9F8BflPtbWs4b+PRsyxu0BRgR0dUpRvt9owmAETEUBqYTACOiqdIC\njIhGMvDLEV9yIwEwIobCOF3giGgow/Rox78EwIgYjmImyGhLAIyIIRHTDPQ+haFLAIyIoSgGQRIA\nI6KBiucAEwAjoqFOpQUYEU2UFmBENJYR0yO+8GQCYEQMTbrAEdFIRrzlRfNdjZ4SACNiKIoHodMF\njoiGyiBIRDSSLaY92i3AgWon6X2Sdkr6fvn3vR3yrJD0hKQXJO2X9JmWtH8r6ZCkveW2YZD6RMRo\nOYUqbf1IOkPSdyV9r4wjf9ohzzpJx1viyW39yh20BbgZ+JbtL0raXB7/67Y8J4HP2d5Trg73rKSd\ntl8o0/+T7f8wYD0iYsQUgyC1dTJ/Afyu7ROSFgPflvSI7afb8j1l++qqhQ7aPt0I3FXu3wV8oj2D\n7cO295T7PwVepFixPSIWsJlBkCpb37IKJ8rDxeU28Mu2Bg2A57asvflDiiUwu5J0AfBB4JmW0/9c\n0j5JWzt1oSNifE1blbYqJC2StBc4Auy0/UyHbJeX8eQRSRf3K7NvAJT0uKTnO2wbW/OVy9J1jciS\n3gN8E/is7dfL018Ffhu4BDgM/Mce109I2i1p9/SJn/WrdkTMs5mZIFU2YMnMv+9ym3hbefa07UuA\n5cBaSe9vy7IHWGn7A8B/AR7sV8e+HXTbH+2WJulHkpbaPixpKUVk7pRvMUXw+4btB1rK/lFLnv8G\nPNyjHpPAJMC7V64Y8ffMRgTAqeqjwEdtr6mS0fZPJD0BrAeebzn/esv+DklfkbTE9tFuZQ3aBd4O\n3FDu3wA81J5BkoCvAS/a/lJb2tKWwz+k5ctExHgrXoZQuQXYk6TfknR2uf8bwMeAl9rynFfGGySt\npYhvx3qVO+gQzReB+yXdBPwA+GT54ecDd9reAHwYuA54ruy/A3ze9g7gzyRdQvHf6hXgnw1Yn4gY\nEUb8sr6pcEuBuyQtoghs99t+WNLNALa3ANcAt0g6CbwJbCpvzXU1UAC0fQz4vQ7n/wbYUO5/Gzo/\n6GP7ukE+PyJGl01tD0Lb3kcxgNp+fkvL/u3A7bMpNzNBImJIqj3kPJ8SACNiKEx9LcBhSQCMiKHJ\nC1EjopGM8kLUiGimYlnM0Q4xo127iBhjWRg9IhrKzGomyLxIAIyIoUkLMCIayVZagBHRTMUgSFaF\ni4hGGv01QRIAI2IoikGQ3AOMiIbKTJCIaKTMBImIRquy4NF8SgCMiKGw4ZenEgAjooGKLnACYEQ0\nVGaCREQjjcNjMLW0TyWtl/SypClJmzukS9KXy/R9ki6tem1EjKuiC1xl61uSdIak70r6nqT9kv60\nQ56ucaabgQNguUrTHcBVwEXAtZIuast2FbCq3CYoFkSvem1EjKlT5bog/bYKfgH8ru3VwCXAekmX\nteXpGGd6qaMLvBaYsn0AQNJ9wEbghZY8G4G7yyXqnpZ0drkm8AUVro2IMVSMAtczF7iMHSfKw8Xl\n1r7kZcc4Y/twt3Lr6AIvA15tOT5YnquSp8q1AEiakLRb0u7pEz8buNIRMVwzD0JX2YAlM/++y22i\nvTxJi8q1xY8AO20/05alcjyZMTaDILYngUmAd69c0XOx44gYDbNYFvOo7TW9MtieBi6RdDawTdL7\nbT8/SP3qCICHgBUtx8vLc1XyLK5wbUSMoWGNAtv+iaQngPVAawCsEot+TR1d4F3AKkkXSjod2ARs\nb8uzHbi+HKW5DDhe9surXBsRY6rGUeDfKlt+SPoN4GPAS23ZusWZrgZuAdo+KelW4DFgEbDV9n5J\nN5fpW4AdwAZgCngDuLHXtYPWKSLmny1O1jcTZClwV/nkyLuA+20/XCXO9FLLPUDbO8oPbz23pWXf\nwKerXhsRC0NdXWDb+4APdjhfKc50MzaDIBExXsZhJkgCYEQMTQJgRDRSXogaEY02i+cA50UCYEQM\nhQ0n80LUiGiqdIEjopFyDzAiGs0JgBHRVBkEiYhGsnMPMCIaS0xnFDgimir3ACOikTIXOCKay8V9\nwFGWABgRQ5NR4IhoJGcQJCKaLF3giGisUR8FrqV9Kmm9pJclTUna3CH9jyTtk/ScpO9IWt2S9kp5\nfq+k3XXUJyLmn10EwCpbP5JWSHpC0guS9kv6TIc86yQdL2PJXkm39St34BZguUjJHRSrNB0Edkna\nbvuFlmx/DXzE9o8lXUWxvu+HWtKvtH100LpExGip8TGYk8DnbO+RdBbwrKSdbXEG4CnbV1cttI4W\n4FpgyvYB228B9wEbWzPY/o7tH5eHT1Os1xkRC5xdbetfjg/b3lPu/xR4EVg2aP3qCIDLgFdbjg/S\nu2I3AY+0HBt4XNKzkia6XSRpQtJuSbunT/xsoApHxPAZcerUuyptwJKZf9/l1isWXECxQtwzHZIv\nL2+3PSLp4n51nNNBEElXUgTAK1pOX2H7kKRzgJ2SXrL9ZPu1ticpus68e+WKER9biggoWjcVHbW9\npl8mSe8Bvgl81vbrbcl7gJW2T0jaADwIrOpVXh0twEPAipbj5eW5XyPpA8CdwEbbx2bO2z5U/j0C\nbKPoUkfEuKtxEARA0mKK4PcN2w+87ePs122fKPd3AIslLelVZh0BcBewStKFkk4HNgHb2yq+EngA\nuM72X7acP7O8oYmkM4GPA8/XUKeIGAWuuPUhScDXgBdtf6lLnvPKfEhaSxHfjnXKO2PgLrDtk5Ju\nBR4DFgFbbe+XdHOZvgW4DfhN4Ctl/U6Wzd1zgW3ludOAe2w/OmidImI01Pgc4IeB64DnJO0tz30e\nWFl8jrcA1wC3SDoJvAlssnsPsdRyD7Bsbu5oO7elZf9TwKc6XHcAWN1+PiLGn4FTp+oJgLa/Db0n\nFtu+Hbh9NuVmJkhEDIeBEZ8JkgAYEUOTucAR0VwJgBHRTNUfcZkvCYARMTxpAUZEIxlc0yjwsCQA\nRsQQJQBGRFOlCxwRjZUAGBGNlAehI6LJ8iB0RDRXRoEjoqmUFmBENFLFd/3NpwTAiBgSZRAkIhos\nLcCIaKxT812B3hIAI2I4xuA5wDoWRULSekkvS5qStLlD+jpJxyXtLbfbql4bEeNLrrb1LUdaIekJ\nSS9I2i/pMx3ySNKXy1iyT9Kl/coduAUoaRFwB/AxikXRd0nabvuFtqxP2b76HV4bEeOovnuAJ4HP\n2d5TriT5rKSdbbHiKop1gFcBHwK+Wv7tqo4W4FpgyvYB228B9wEb5+DaiGgI24dt7yn3fwq8CCxr\ny7YRuNuFp4GzJS3tVW4d9wCXAa+2HB+kc9S9XNI+ikXT/9j2/llci6QJYAJg5bLT+Kt/sqVTthhR\nv3/+JfNdhZiFt1zP6MUsHoReIml3y/Gk7cmOZUoXAB8EnmlL6hRPlgGHu33oXA2C7AFW2j4haQPw\nIEUztbLyP8YkwJrVZ4z44HpEFOtiVh4EOVquFd6TpPcA3wQ+a/v1AWoH1NMFPgSsaDleXp77Fduv\n2z5R7u8AFktaUuXaiBhjrrhVIGkxRfD7hu0HOmSZdTypIwDuAlZJulDS6cAmYHtrBknnSVK5v7b8\n3GNVro2I8VXjKLCArwEv2v5Sl2zbgevL0eDLgOO2u3Z/oYYusO2Tkm4FHgMWAVtt75d0c5m+BbgG\nuEXSSeBNYJNtAx2vHbROETEi6rtZ9WHgOuA5SXvLc58HVsKv4swOYAMwBbwB3Niv0FruAZbd2h1t\n57a07N8O3F712ohYIGoKgLa/TZ8FRspG1adnU25mgkTEUFTt3s6nBMCIGJ68EDUimiotwIhorgTA\niGik3AOMiEZLAIyIptKIvxC1lvcBRkSMo7QAI2J40gWOiEbKIEhENFoCYEQ0VgJgRDSRGP1R4ATA\niBiO3AOMiEZLAIyIxkoAjIimShc4IpprxANgLVPhJK2X9LKkKUmbO6T/S0l7y+15SdOS3lemvSLp\nuTJt99tLj4ix5GIUuMrWj6Stko5Ier5L+jpJx1vizG1VqjhwC1DSIuAO4GMUCxHvkrTd9gszeWz/\nOfDnZf4/AP6F7ddairnS9tFB6xIRI6a+FuDXKdYVurtHnqdsXz2bQutoAa4FpmwfsP0WcB+wsUf+\na4F7a/jciBhxdS2LaftJ4LW+GWepjgC4DHi15fhgee5tJP0dYD3F4sYzDDwu6VlJE90+RNKEpN2S\ndv/tsekaqh0RQ1d9YfQlM/++y61rLOjhckn7JD0i6eIqF8z1IMgfAP+nrft7he1Dks4Bdkp6qYz2\nv8b2JDAJsGb1GSN+azUiWoJbFUdtrxng0/YAK22fkLQBeBBY1e+iOlqAh4AVLcfLy3OdbKKt+2v7\nUPn3CLCNoksdEWNO1NcF7sf267ZPlPs7gMWSlvS7ro4AuAtYJelCSadTBLnt7Zkk/T3gI8BDLefO\nlHTWzD7wcaDjKE9EjJ+5CoCSzpOkcn8tRWw71u+6gbvAtk9KuhV4DFgEbLW9X9LNZfqWMusfAv/L\n9s9aLj8X2FbW+zTgHtuPDlqniBgRNd2sknQvsI7iXuFB4AvAYvhVjLkGuEXSSeBNYJPtvp9eyz3A\nssm5o+3clrbjr1MMZbeeOwCsrqMOETGCagqAtq/tk347xWMys5KZIBExHHkbTEQ0WgJgRDRVXoga\nEY2VLnBENNPsHoSeFwmAETE8CYAR0UQzM0FGWQJgRAyNTo12BEwAjIjhyD3AiGiydIEjorkSACOi\nqdICjIjmSgCMiEZypsJFREPlOcCIaLb+7ySdVwmAETE0aQFGRDONwYPQdSyKhKStko5I6rigkQpf\nljRVrtt5aUvaekkvl2mb66hPRIwGnaq29S1ngBjTSy0BkGKtj/U90q+iWKNzFTABfBVA0iLgjjL9\nIuBaSRfVVKeImGd1BUDeYYzpp5YAWC5k/lqPLBuBu114Gjhb0lKKNYCnbB+w/RZwX5k3IsadKQZB\nqmz9inrnMaanulqA/SwDXm05Plie63b+bSRNSNotafffHpseWkUjoj6zWBd4ycy/73KbmOVHVY4l\nrcZmEMT2JDAJsGb1GSN+azUigNkMghy1vWaINelorgLgIWBFy/Hy8tziLucjYszN8YPQ3WJMT3PV\nBd4OXF+O1FwGHLd9GNgFrJJ0oaTTgU1l3ogYdzY6VW2rQbcY01MtLUBJ9wLrKPrxB4EvULTusL0F\n2AFsAKaAN4Aby7STkm4FHgMWAVtt76+jThExAmpqAb7TGNNPLQHQ9rV90g18ukvaDorKR8QCU1cX\neJAY08vYDIJExJgxkDVBIqKxRjv+JQBGxPDkZQgR0VhZFjMimmkM3gaTABgRQ1E8CD3aETABMCKG\nJ2uCRERTpQUYEc2Ue4AR0Vy1zfMdmgTAiBiedIEjopGyMHpENFpagBHRWKMd/xIAI2J4dGq0+8AJ\ngBExHCYPQkdEMwnnQeiIaLARD4C1LIokaaukI5Ke75L+R5L2SXpO0nckrW5Je6U8v1fS7jrqExEj\noqaF0QEkrZf0sqQpSZs7pK+TdLyMJXsl3davzLpagF8Hbgfu7pL+18BHbP9Y0lUU6/t+qCX9SttH\na6pLRIyCGu8BSloE3AF8jGLR812Sttt+oS3rU7avrlpuXYsiPSnpgh7p32k5fJpizc6IWOBqHAVe\nC0zZPgAg6T5gI9AeAGdlrtYFbnUT8EjLsYHHJT0raWIe6hMRQ1Gx+1utC7wMeLXl+GB5rt3l5e22\nRyRd3K/QOR0EkXQlRQC8ouX0FbYPSToH2CnpJdtPdrh2ApgAWLksYzcRI8/MZhBkSdsYwKTtyVl+\n4h5gpe0TkjYADwKrel0wZ5FE0geAO4GrbB+bOW/7UPn3iKRtFE3dtwXA8j/GJMCa1WeM9tBSRBSq\n94CP2l7TI/0QsKLleHl57ldsv96yv0PSVyQt6TW+MCddYEkrgQeA62z/Zcv5MyWdNbMPfBzoOJIc\nEeNHdqWtgl3AKkkXSjod2ARs/7XPks6TpHJ/LUV8O/a2klrU0gKUdC+wjqIZexD4ArAYwPYW4Dbg\nN4GvlPU7WUb7c4Ft5bnTgHtsP1pHnSJiBNT0HKDtk5JuBR4DFgFbbe+XdHOZvgW4BrhF0kngTWCT\n3bsCdY0CX9sn/VPApzqcPwCsfvsVETH2bJiuby6c7R3AjrZzW1r2b6d4HK+yjCZExPCM+EyQBMCI\nGJ4EwIhoJANZEyQimsng0X4fVgJgRAyHqXUQZBgSACNieHIPMCIaKwEwIpqp+rv+5ksCYEQMh4Es\nihQRjZUWYEQ0U71T4YYhATAihsPgPAcYEY2VmSAR0Vi5BxgRjWRnFDgiGiwtwIhoJuPp6fmuRE8J\ngBExHHkdVkQ02og/BlPLqnCStko6Iqnjim6S1kk6Lmlvud3WkrZe0suSpiRtrqM+ETH/DPiUK21V\n9IsVKny5TN8n6dJ+Zda1LObXgfV98jxl+5Jy+3cAkhYBdwBXARcB10q6qKY6RcR8cvlC1CpbHxVj\nxVUUC6GvAiaAr/Yrt5YAaPtJ4LV3cOlaYMr2AdtvAfcBG+uoU0TMP09PV9oqqBIrNgJ3u/A0cLak\npb0Knct7gJdL2kexmvsf294PLANebclzEPhQp4slTVBEdYBfLFo6tRAXUF8CdF3FfrxNLdTvtlC/\n198ftICf8uPHHvf/WFIx+xmSdrccT9qebDmuEis65VkGHO72oXMVAPcAK22fkLQBeJCimVpZ+R9j\nEkDS7nJh9QVloX4vWLjfbSF/r0HLsN3vtti8q+seYE+2X7d9otzfASyWtISiNbiiJevy8lxERKsq\nsWLW8WROAqCk8ySp3F9bfu4xYBewStKFkk4HNgHb56JOETFWqsSK7cD15WjwZcBx2127v1BTF1jS\nvcA6YImkg8AXgMUAtrcA1wC3SDoJvAlssm3gpKRbgceARcDW8t5gP5P9s4ylhfq9YOF+t3yvOWC7\nY6yQdHOZvgXYAWwApoA3gBv7lSuP+Fy9iIhhmZMucETEKEoAjIjGGosAKOl9knZK+n75971d8r0i\n6blyut3Aw/jDMowpPaOgwvfqOiVylFWY6jmWvxcMNo11QbA98hvwZ8Dmcn8z8O+75HsFWDLf9e3z\nXRYBfwX8NnA68D3gorY8G4BHAAGXAc/Md71r+l7rgIfnu67v4Lv9DnAp8HyX9LH7vWbx3cbyN6u6\njUULkGKKy13l/l3AJ+axLoMaypSeEbBgpzW6/1TPcfy9gIGmsS4I4xIAz/X/f57nh8C5XfIZeFzS\ns+XUuVHUbbrObPOMmqp1vrzsJj4i6eK5qdrQjePvNRsL8TcDRuh9gJIeB87rkPQnrQe2LanbsztX\n2D4k6Rxgp6SXyv/DxWgYeEpkzLkF/ZuNTAvQ9kdtv7/D9hDwo5kuRfn3SJcyDpV/jwDbKLplo2Yo\nU3pGQN86u/uUyHE3jr9XJQv4NwNGKAD2sR24ody/AXioPYOkMyWdNbMPfBwYxTfGDGVKzwjo+716\nTIkcd+P4e1WygH8zYIS6wH18Ebhf0k3AD4BPAkg6H7jT9gaK+4Lbyt/qNOAe24/OU3278pCm9My3\nit+r25TIkVZhqufY/V4zBpjGuiBkKlxENNa4dIEjImqXABgRjZUAGBGNlQAYEY2VABgRjZUAGBGN\nlQAYEY31/wDgym3FCz7j2QAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7f88c92406a0>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "#plt.plot(increment_by_one(array.ravel()))\n",
    "plt.imshow(array.reshape([n,n]).real[n/2-1:n/2+1,n/2-1:n/2+1])#alpha=.5)\n",
    "\n",
    "plt.colorbar()\n",
    "#plt.plot(unshifted, alpha=.5)\n",
    "\n",
    "#print(flatarr.shape)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "877 ms ± 3.98 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)\n",
      "442 ms ± 74.3 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)\n"
     ]
    }
   ],
   "source": [
    "%timeit np.fft.fftshift(array)\n",
    "%timeit cufftShift_2D_kernel[numBlocks, blockdim](array,n)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "from numba import vectorize, float64\n",
    "\n",
    "@vectorize([float64(float64, float64)])\n",
    "def f(x, y):\n",
    "    return x + y"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### check slowdown of units\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "4.51 µs ± 14.8 ns per loop (mean ± std. dev. of 7 runs, 100000 loops each)\n",
      "6.7 µs ± 63.1 ns per loop (mean ± std. dev. of 7 runs, 100000 loops each)\n",
      "44 µs ± 1.02 µs per loop (mean ± std. dev. of 7 runs, 10000 loops each)\n",
      "292 µs ± 2.62 µs per loop (mean ± std. dev. of 7 runs, 1000 loops each)\n",
      "361 µs ± 7.14 µs per loop (mean ± std. dev. of 7 runs, 1000 loops each)\n"
     ]
    }
   ],
   "source": [
    "import astropy.units as u\n",
    "x1=np.array([1])\n",
    "y1=np.array([2])\n",
    "z1=np.array([3])\n",
    "xu=x1*u.m\n",
    "yu=y1*u.m\n",
    "zu=z1*u.m\n",
    "def fun(x,y,z):\n",
    "    (x**2+y**2)/np.sqrt(z)\n",
    "    \n",
    "%timeit   fun(x1,y1,z1)\n",
    "%timeit   fun(xu.value, yu.value,zu.value)\n",
    "%timeit   fun(xu.to(u.m).value, yu.to(u.m).value, zu.to(u.m).value)\n",
    "%timeit   fun(xu,yu,zu)\n",
    "%timeit   fun(x1*u.m,y1*u.m,z1*u.m)\n",
    "\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Tentative Conclusions:\n",
    "- Units can slow things down 100x\n",
    "- fftshift calls np.ndarray.take which is slow\n",
    "- cuda FFT on K80 is >5x faster than MKL fft or pyfftw on Xeon\n",
    "- numexpr is generally 10x faster numpy and 2x faster than numba (this was also found by http://nbviewer.jupyter.org/github/rasbt/One-Python-benchmark-per-day/blob/master/ipython_nbs/day7_2_jit_numpy.ipynb) "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  },
  {
   "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.6.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}