_notebooks/2020-08-04-dask-array-rounding.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Dask array rounding\n",
    "> Dask array rounding float32\n",
    "- categories: [python, dask]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This notebook is an extract from the [Dask array Tutorial notebook](https://github.com/dask/dask-tutorial/blob/master/03_array.ipynb), see also the Youtube SciPy 2020 class at <https://www.youtube.com/watch?v=mqdglv9GnM8>.\n",
    "\n",
    "We notice that `dask` is automatically rounding `float32` numbers to machine precision, which I think is the most sensible choice, but surprising difference compared to `numpy`."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Create the input data\n",
    "\n",
    "Needs only to be done once, defaults to ~4GB of data, you can reduce it by setting `blocksize` to a smaller number, e.g. 1000"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "import os\n",
    "import h5py\n",
    "import numpy as np"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%time\n",
    "\n",
    "blocksize = 1000000\n",
    "nblocks = 1000\n",
    "shape = nblocks * blocksize\n",
    "\n",
    "if not os.path.exists('random.hdf5'):\n",
    "\n",
    "    with h5py.File('random.hdf5', mode='w') as f:\n",
    "        dset = f.create_dataset('/x', shape=(shape,), dtype='f4')\n",
    "        for i in range(0, shape, blocksize):\n",
    "            dset[i: i + blocksize] = np.random.exponential(size=blocksize)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Setup"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "from dask.distributed import Client\n",
    "\n",
    "client = Client(n_workers=24, processes=False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Load data with h5py\n",
    "# this creates a pointer to the data, but does not actually load\n",
    "import h5py\n",
    "import os\n",
    "f = h5py.File('random.hdf5', mode='r')\n",
    "dset = f['/x']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "dtype('<f4')"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "dset.dtype"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "-rw-r--r-- 1 zonca csb148 3.8G Jul 24 22:51 data/random.hdf5\n"
     ]
    }
   ],
   "source": [
    "!ls -lah data/random.hdf5"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Compute sum using blocked algorithm**"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Before using dask, lets consider the concept of blocked algorithms. We can compute the sum of a large number of elements by loading them chunk-by-chunk, and keeping a running total.\n",
    "\n",
    "Here we compute the sum of this large array on disk by \n",
    "\n",
    "1.  Computing the sum of each 1,000,000 sized chunk of the array\n",
    "2.  Computing the sum of the 1,000 intermediate sums\n",
    "\n",
    "Note that this is a sequential process in the notebook kernel, both the loading and summing."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "1000000000"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "len(dset)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "999976587.6875\n"
     ]
    }
   ],
   "source": [
    "# Compute sum of large array, one million numbers at a time\n",
    "sums = []\n",
    "for i in range(0, 1000000000, 1000000):\n",
    "    chunk = dset[i: i + 1000000]  # pull out numpy array\n",
    "    sums.append(chunk.sum())\n",
    "\n",
    "total = sum(sums)\n",
    "print(total)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Create `dask.array` object**"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You can create a `dask.array` `Array` object with the `da.from_array` function.  This function accepts\n",
    "\n",
    "1.  `data`: Any object that supports NumPy slicing, like `dset`\n",
    "2.  `chunks`: A chunk size to tell us how to block up our array, like `(1000000,)`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<table>\n",
       "<tr>\n",
       "<td>\n",
       "<table>\n",
       "  <thead>\n",
       "    <tr><td> </td><th> Array </th><th> Chunk </th></tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr><th> Bytes </th><td> 4.00 GB </td> <td> 40.00 MB </td></tr>\n",
       "    <tr><th> Shape </th><td> (1000000000,) </td> <td> (10000000,) </td></tr>\n",
       "    <tr><th> Count </th><td> 101 Tasks </td><td> 100 Chunks </td></tr>\n",
       "    <tr><th> Type </th><td> float32 </td><td> numpy.ndarray </td></tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</td>\n",
       "<td>\n",
       "<svg width=\"170\" height=\"75\" style=\"stroke:rgb(0,0,0);stroke-width:1\" >\n",
       "\n",
       "  <!-- Horizontal lines -->\n",
       "  <line x1=\"0\" y1=\"0\" x2=\"120\" y2=\"0\" style=\"stroke-width:2\" />\n",
       "  <line x1=\"0\" y1=\"25\" x2=\"120\" y2=\"25\" style=\"stroke-width:2\" />\n",
       "\n",
       "  <!-- Vertical lines -->\n",
       "  <line x1=\"0\" y1=\"0\" x2=\"0\" y2=\"25\" style=\"stroke-width:2\" />\n",
       "  <line x1=\"1\" y1=\"0\" x2=\"1\" y2=\"25\" />\n",
       "  <line x1=\"2\" y1=\"0\" x2=\"2\" y2=\"25\" />\n",
       "  <line x1=\"3\" y1=\"0\" x2=\"3\" y2=\"25\" />\n",
       "  <line x1=\"4\" y1=\"0\" x2=\"4\" y2=\"25\" />\n",
       "  <line x1=\"6\" y1=\"0\" x2=\"6\" y2=\"25\" />\n",
       "  <line x1=\"7\" y1=\"0\" x2=\"7\" y2=\"25\" />\n",
       "  <line x1=\"8\" y1=\"0\" x2=\"8\" y2=\"25\" />\n",
       "  <line x1=\"9\" y1=\"0\" x2=\"9\" y2=\"25\" />\n",
       "  <line x1=\"10\" y1=\"0\" x2=\"10\" y2=\"25\" />\n",
       "  <line x1=\"12\" y1=\"0\" x2=\"12\" y2=\"25\" />\n",
       "  <line x1=\"13\" y1=\"0\" x2=\"13\" y2=\"25\" />\n",
       "  <line x1=\"14\" y1=\"0\" x2=\"14\" y2=\"25\" />\n",
       "  <line x1=\"15\" y1=\"0\" x2=\"15\" y2=\"25\" />\n",
       "  <line x1=\"16\" y1=\"0\" x2=\"16\" y2=\"25\" />\n",
       "  <line x1=\"18\" y1=\"0\" x2=\"18\" y2=\"25\" />\n",
       "  <line x1=\"19\" y1=\"0\" x2=\"19\" y2=\"25\" />\n",
       "  <line x1=\"20\" y1=\"0\" x2=\"20\" y2=\"25\" />\n",
       "  <line x1=\"21\" y1=\"0\" x2=\"21\" y2=\"25\" />\n",
       "  <line x1=\"22\" y1=\"0\" x2=\"22\" y2=\"25\" />\n",
       "  <line x1=\"24\" y1=\"0\" x2=\"24\" y2=\"25\" />\n",
       "  <line x1=\"25\" y1=\"0\" x2=\"25\" y2=\"25\" />\n",
       "  <line x1=\"26\" y1=\"0\" x2=\"26\" y2=\"25\" />\n",
       "  <line x1=\"27\" y1=\"0\" x2=\"27\" y2=\"25\" />\n",
       "  <line x1=\"28\" y1=\"0\" x2=\"28\" y2=\"25\" />\n",
       "  <line x1=\"30\" y1=\"0\" x2=\"30\" y2=\"25\" />\n",
       "  <line x1=\"31\" y1=\"0\" x2=\"31\" y2=\"25\" />\n",
       "  <line x1=\"32\" y1=\"0\" x2=\"32\" y2=\"25\" />\n",
       "  <line x1=\"33\" y1=\"0\" x2=\"33\" y2=\"25\" />\n",
       "  <line x1=\"34\" y1=\"0\" x2=\"34\" y2=\"25\" />\n",
       "  <line x1=\"36\" y1=\"0\" x2=\"36\" y2=\"25\" />\n",
       "  <line x1=\"37\" y1=\"0\" x2=\"37\" y2=\"25\" />\n",
       "  <line x1=\"38\" y1=\"0\" x2=\"38\" y2=\"25\" />\n",
       "  <line x1=\"39\" y1=\"0\" x2=\"39\" y2=\"25\" />\n",
       "  <line x1=\"40\" y1=\"0\" x2=\"40\" y2=\"25\" />\n",
       "  <line x1=\"42\" y1=\"0\" x2=\"42\" y2=\"25\" />\n",
       "  <line x1=\"43\" y1=\"0\" x2=\"43\" y2=\"25\" />\n",
       "  <line x1=\"44\" y1=\"0\" x2=\"44\" y2=\"25\" />\n",
       "  <line x1=\"45\" y1=\"0\" x2=\"45\" y2=\"25\" />\n",
       "  <line x1=\"46\" y1=\"0\" x2=\"46\" y2=\"25\" />\n",
       "  <line x1=\"48\" y1=\"0\" x2=\"48\" y2=\"25\" />\n",
       "  <line x1=\"49\" y1=\"0\" x2=\"49\" y2=\"25\" />\n",
       "  <line x1=\"50\" y1=\"0\" x2=\"50\" y2=\"25\" />\n",
       "  <line x1=\"51\" y1=\"0\" x2=\"51\" y2=\"25\" />\n",
       "  <line x1=\"52\" y1=\"0\" x2=\"52\" y2=\"25\" />\n",
       "  <line x1=\"54\" y1=\"0\" x2=\"54\" y2=\"25\" />\n",
       "  <line x1=\"55\" y1=\"0\" x2=\"55\" y2=\"25\" />\n",
       "  <line x1=\"56\" y1=\"0\" x2=\"56\" y2=\"25\" />\n",
       "  <line x1=\"57\" y1=\"0\" x2=\"57\" y2=\"25\" />\n",
       "  <line x1=\"58\" y1=\"0\" x2=\"58\" y2=\"25\" />\n",
       "  <line x1=\"60\" y1=\"0\" x2=\"60\" y2=\"25\" />\n",
       "  <line x1=\"61\" y1=\"0\" x2=\"61\" y2=\"25\" />\n",
       "  <line x1=\"62\" y1=\"0\" x2=\"62\" y2=\"25\" />\n",
       "  <line x1=\"63\" y1=\"0\" x2=\"63\" y2=\"25\" />\n",
       "  <line x1=\"64\" y1=\"0\" x2=\"64\" y2=\"25\" />\n",
       "  <line x1=\"66\" y1=\"0\" x2=\"66\" y2=\"25\" />\n",
       "  <line x1=\"67\" y1=\"0\" x2=\"67\" y2=\"25\" />\n",
       "  <line x1=\"68\" y1=\"0\" x2=\"68\" y2=\"25\" />\n",
       "  <line x1=\"69\" y1=\"0\" x2=\"69\" y2=\"25\" />\n",
       "  <line x1=\"70\" y1=\"0\" x2=\"70\" y2=\"25\" />\n",
       "  <line x1=\"72\" y1=\"0\" x2=\"72\" y2=\"25\" />\n",
       "  <line x1=\"73\" y1=\"0\" x2=\"73\" y2=\"25\" />\n",
       "  <line x1=\"74\" y1=\"0\" x2=\"74\" y2=\"25\" />\n",
       "  <line x1=\"75\" y1=\"0\" x2=\"75\" y2=\"25\" />\n",
       "  <line x1=\"76\" y1=\"0\" x2=\"76\" y2=\"25\" />\n",
       "  <line x1=\"78\" y1=\"0\" x2=\"78\" y2=\"25\" />\n",
       "  <line x1=\"79\" y1=\"0\" x2=\"79\" y2=\"25\" />\n",
       "  <line x1=\"80\" y1=\"0\" x2=\"80\" y2=\"25\" />\n",
       "  <line x1=\"81\" y1=\"0\" x2=\"81\" y2=\"25\" />\n",
       "  <line x1=\"82\" y1=\"0\" x2=\"82\" y2=\"25\" />\n",
       "  <line x1=\"84\" y1=\"0\" x2=\"84\" y2=\"25\" />\n",
       "  <line x1=\"85\" y1=\"0\" x2=\"85\" y2=\"25\" />\n",
       "  <line x1=\"86\" y1=\"0\" x2=\"86\" y2=\"25\" />\n",
       "  <line x1=\"87\" y1=\"0\" x2=\"87\" y2=\"25\" />\n",
       "  <line x1=\"88\" y1=\"0\" x2=\"88\" y2=\"25\" />\n",
       "  <line x1=\"90\" y1=\"0\" x2=\"90\" y2=\"25\" />\n",
       "  <line x1=\"91\" y1=\"0\" x2=\"91\" y2=\"25\" />\n",
       "  <line x1=\"92\" y1=\"0\" x2=\"92\" y2=\"25\" />\n",
       "  <line x1=\"93\" y1=\"0\" x2=\"93\" y2=\"25\" />\n",
       "  <line x1=\"94\" y1=\"0\" x2=\"94\" y2=\"25\" />\n",
       "  <line x1=\"96\" y1=\"0\" x2=\"96\" y2=\"25\" />\n",
       "  <line x1=\"97\" y1=\"0\" x2=\"97\" y2=\"25\" />\n",
       "  <line x1=\"98\" y1=\"0\" x2=\"98\" y2=\"25\" />\n",
       "  <line x1=\"99\" y1=\"0\" x2=\"99\" y2=\"25\" />\n",
       "  <line x1=\"100\" y1=\"0\" x2=\"100\" y2=\"25\" />\n",
       "  <line x1=\"102\" y1=\"0\" x2=\"102\" y2=\"25\" />\n",
       "  <line x1=\"103\" y1=\"0\" x2=\"103\" y2=\"25\" />\n",
       "  <line x1=\"104\" y1=\"0\" x2=\"104\" y2=\"25\" />\n",
       "  <line x1=\"105\" y1=\"0\" x2=\"105\" y2=\"25\" />\n",
       "  <line x1=\"106\" y1=\"0\" x2=\"106\" y2=\"25\" />\n",
       "  <line x1=\"108\" y1=\"0\" x2=\"108\" y2=\"25\" />\n",
       "  <line x1=\"109\" y1=\"0\" x2=\"109\" y2=\"25\" />\n",
       "  <line x1=\"110\" y1=\"0\" x2=\"110\" y2=\"25\" />\n",
       "  <line x1=\"111\" y1=\"0\" x2=\"111\" y2=\"25\" />\n",
       "  <line x1=\"112\" y1=\"0\" x2=\"112\" y2=\"25\" />\n",
       "  <line x1=\"114\" y1=\"0\" x2=\"114\" y2=\"25\" />\n",
       "  <line x1=\"115\" y1=\"0\" x2=\"115\" y2=\"25\" />\n",
       "  <line x1=\"116\" y1=\"0\" x2=\"116\" y2=\"25\" />\n",
       "  <line x1=\"117\" y1=\"0\" x2=\"117\" y2=\"25\" />\n",
       "  <line x1=\"118\" y1=\"0\" x2=\"118\" y2=\"25\" />\n",
       "  <line x1=\"120\" y1=\"0\" x2=\"120\" y2=\"25\" style=\"stroke-width:2\" />\n",
       "\n",
       "  <!-- Colored Rectangle -->\n",
       "  <polygon points=\"0.000000,0.000000 120.000000,0.000000 120.000000,25.412617 0.000000,25.412617\" style=\"fill:#ECB172A0;stroke-width:0\"/>\n",
       "\n",
       "  <!-- Text -->\n",
       "  <text x=\"60.000000\" y=\"45.412617\" font-size=\"1.0rem\" font-weight=\"100\" text-anchor=\"middle\" >1000000000</text>\n",
       "  <text x=\"140.000000\" y=\"12.706308\" font-size=\"1.0rem\" font-weight=\"100\" text-anchor=\"middle\" transform=\"rotate(0,140.000000,12.706308)\">1</text>\n",
       "</svg>\n",
       "</td>\n",
       "</tr>\n",
       "</table>"
      ],
      "text/plain": [
       "dask.array<array, shape=(1000000000,), dtype=float32, chunksize=(10000000,), chunktype=numpy.ndarray>"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "import dask.array as da\n",
    "import numpy as np\n",
    "x = da.from_array(dset, chunks=(10000000,))\n",
    "x"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [],
   "source": [
    "x_float64 = x.astype(np.float64)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "999976700.0"
      ]
     },
     "execution_count": 18,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "x.sum().compute()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The machine resolution of `float32` is `1e-6`, therefore everything after the 7th digit is garbage, so it is reasonable to remove it, otherwise it gives you the impression that the computation is more precise than it actually it.\n",
    "Still I am surprised `dask` does it, `numpy` above instead doesn't care about that and prints all the digits.\n",
    "\n",
    "If we need more precision, we need to increase the precision of the calculation, see below, but we are going to use a lot more memory, also, the input data were `float32`, so it is not very useful anyway, we should generate again the input with higher precision."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "finfo(resolution=1e-06, min=-3.4028235e+38, max=3.4028235e+38, dtype=float32)"
      ]
     },
     "execution_count": 19,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.finfo(np.float32)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<table>\n",
       "<tr>\n",
       "<td>\n",
       "<table>\n",
       "  <thead>\n",
       "    <tr><td> </td><th> Array </th><th> Chunk </th></tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr><th> Bytes </th><td> 8 B </td> <td> 8 B </td></tr>\n",
       "    <tr><th> Shape </th><td> () </td> <td> () </td></tr>\n",
       "    <tr><th> Count </th><td> 336 Tasks </td><td> 1 Chunks </td></tr>\n",
       "    <tr><th> Type </th><td> float64 </td><td> numpy.ndarray </td></tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</td>\n",
       "<td>\n",
       "\n",
       "</td>\n",
       "</tr>\n",
       "</table>"
      ],
      "text/plain": [
       "dask.array<sum-aggregate, shape=(), dtype=float64, chunksize=(), chunktype=numpy.ndarray>"
      ]
     },
     "execution_count": 20,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "x_float64.sum()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "999976584.1788422"
      ]
     },
     "execution_count": 21,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "x_float64.sum().compute()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "client.shutdown()"
   ]
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "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.7.8"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}