{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Fun with Floating Point Precision in numpy"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "I recently had a bug in my code that obviously was caused by an issue with floating point precision but had me scratching my head how it came about. Here it is:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "from astropy.table import Table\n",
    "from astropy import cosmology\n",
    "\n",
    "cosmo = cosmology.WMAP9"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "`Table` lets me read a FITS table, the standard data format in Astronomy. I draw a random sample from a table column into a `numpy.ndarray`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Downloading https://github.com/joergdietrich/joergdietrich.github.io/raw/master/data/float.fits [Done]\n"
     ]
    }
   ],
   "source": [
    "tab = Table.read(\"https://github.com/joergdietrich/joergdietrich.github.io/raw/master/data/float.fits\")\n",
    "z = np.random.choice(tab['x'], 500)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "I needed to call some `astropy` code (`angular_diameter_distance_z1z2(z1, z2)`), which takes two arrays are argument and requires that all values in `z1` are less or equal than the values in `z2`. In my case `z1` is a constant and I enforce this constraint by setting by setting all smaller values equal to this constant."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "False"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "zl = 0.32\n",
    "idx = z < zl\n",
    "z[idx] = zl\n",
    "np.any(z < zl)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "So far, so good. The code actually refused to compare a scalar with an array, so I called it like this:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "ename": "ValueError",
     "evalue": "z2 must greater than z1",
     "output_type": "error",
     "traceback": [
      "\u001b[1;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[1;31mValueError\u001b[0m                                Traceback (most recent call last)",
      "\u001b[1;32m<ipython-input-4-f71b9851644c>\u001b[0m in \u001b[0;36m<module>\u001b[1;34m()\u001b[0m\n\u001b[1;32m----> 1\u001b[1;33m \u001b[0mdls\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mcosmo\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mangular_diameter_distance_z1z2\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mnp\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mrepeat\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mzl\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mz\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0msize\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mz\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m",
      "\u001b[1;32m/home/joerg/applications/anaconda3/lib/python3.5/site-packages/astropy-1.1.dev14570-py3.5-linux-x86_64.egg/astropy/cosmology/core.py\u001b[0m in \u001b[0;36mangular_diameter_distance_z1z2\u001b[1;34m(self, z1, z2)\u001b[0m\n\u001b[0;32m   1333\u001b[0m \u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m   1334\u001b[0m         \u001b[1;32mif\u001b[0m \u001b[1;33m(\u001b[0m\u001b[0mz1\u001b[0m \u001b[1;33m>\u001b[0m \u001b[0mz2\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0many\u001b[0m\u001b[1;33m(\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m-> 1335\u001b[1;33m             \u001b[1;32mraise\u001b[0m \u001b[0mValueError\u001b[0m\u001b[1;33m(\u001b[0m\u001b[1;34m'z2 must greater than z1'\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[0;32m   1336\u001b[0m \u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m   1337\u001b[0m         \u001b[0mdm1\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mself\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mcomoving_transverse_distance\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mz1\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mvalue\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n",
      "\u001b[1;31mValueError\u001b[0m: z2 must greater than z1"
     ]
    }
   ],
   "source": [
    "dls = cosmo.angular_diameter_distance_z1z2(np.repeat(zl, z.size), z)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Oops! What happended? I had just set everything smaller than `zl` to `zl` and our check above confirmed that there are no values smaller than `zl` in `z`. Let's check again:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(True, False)"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.any(np.repeat(zl, z.size) > z), np.any(zl > z)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Weird. What ared these values?\n",
    "\n",
    "The check works for scalar/array but not for array/array. Actually the check of course works, but there's a subtlety going on, as we'll see in a moment."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "True"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.all(np.equal(np.where(np.repeat(zl, z.size) > z), np.where(idx)))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "So culprits are the elements we just set to `zl`. Lets see what they are."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(array([ 0.31999999,  0.31999999,  0.31999999,  0.31999999,  0.31999999,\n",
       "         0.31999999,  0.31999999,  0.31999999,  0.31999999,  0.31999999,\n",
       "         0.31999999,  0.31999999,  0.31999999,  0.31999999,  0.31999999,\n",
       "         0.31999999,  0.31999999,  0.31999999,  0.31999999,  0.31999999,\n",
       "         0.31999999,  0.31999999,  0.31999999,  0.31999999,  0.31999999,\n",
       "         0.31999999,  0.31999999,  0.31999999,  0.31999999,  0.31999999,\n",
       "         0.31999999,  0.31999999,  0.31999999,  0.31999999,  0.31999999,\n",
       "         0.31999999,  0.31999999,  0.31999999,  0.31999999,  0.31999999,\n",
       "         0.31999999,  0.31999999,  0.31999999,  0.31999999,  0.31999999,\n",
       "         0.31999999,  0.31999999,  0.31999999,  0.31999999], dtype=float32),\n",
       " array([ 0.32,  0.32,  0.32,  0.32,  0.32,  0.32,  0.32,  0.32,  0.32,\n",
       "         0.32,  0.32,  0.32,  0.32,  0.32,  0.32,  0.32,  0.32,  0.32,\n",
       "         0.32,  0.32,  0.32,  0.32,  0.32,  0.32,  0.32,  0.32,  0.32,\n",
       "         0.32,  0.32,  0.32,  0.32,  0.32,  0.32,  0.32,  0.32,  0.32,\n",
       "         0.32,  0.32,  0.32,  0.32,  0.32,  0.32,  0.32,  0.32,  0.32,\n",
       "         0.32,  0.32,  0.32,  0.32]))"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "z[idx], np.repeat(zl, z.size)[idx]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Clearly a round-off error. But why? We just set the array elements to exactly this value. Wait. Notice that the first array has a `dtype=float32` but the second array does not have a `dtype` specified, i.e., it is a standard `float64` array. This is where the problem is. `z` came from a FITS table, which has a 32 bit wide float column and that's just what astropy gave us. Comparing to a 64 bit float array led to inevitable round-off errors, while apparently the correct cast is made when comparing a 32 bit float array to a 64 bit scalar."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "False"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "z = np.float64(z)\n",
    "z[idx] = zl\n",
    "np.any(np.repeat(zl, z.size) > np.float64(z))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Converting to float64 before the assignment fixes this problem. This took me a while to figure out.\n",
    "\n",
    "<small>\n",
    "<a rel=\"license\" href=\"http://creativecommons.org/licenses/by-sa/4.0/\">This blog post is licensed under a <a rel=\"license\" href=\"http://creativecommons.org/licenses/by-sa/4.0/\">Creative Commons Attribution-ShareAlike 4.0 International License</a>. This post was written as a <a href=\"http://www.jupyter.org/\">Jupyter Notebook</a> in Python. You can <a href=\"https://raw.githubusercontent.com/joergdietrich/joergdietrich.github.io/master/notebooks/float_precision.ipynb\">download</a> the original notebook.</small>"
   ]
  }
 ],
 "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.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}