Bitshuffle/LZ4 and precision reduction

precision.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Reducing the precision of floats to improve the compression rate ...\n",
    "\n",
    "This notebook tried to evaluate compression rate of bitshuffle and LZ4 scheme, when coupled with precision reduction, i.e. reduction of the number of bits of the mantissa of the floating point.\n",
    "\n",
    "In this document we will work with IEEE754 floating point values, mainly:\n",
    "* Single precision float, stored in 32 bits (4 bytes) and containing 23+1 bits of mantissa for a precision of 1.2e-7\n",
    "* Double precision float, stored in 64 bits (8 bytes) and containing 52+1 bits of mantissa for a precision of 2.2e-16\n",
    "\n",
    "## Rational:\n",
    "* Our detector work with ADC of precision 12-20 bits, often. \n",
    "* For sake on convieniance, most scientists perform their calculation in double precision floating points which offer 53 bits of mantissa (64 bits)\n",
    "* This represents a large increase in the size of the data. Moreover floating point data are hardly compressible.\n",
    "\n",
    "We focuses only on LZ4 compression because it is already used in Eiger detector on the one hand and it provides a high speed compression and decompression on the other making it suitable for on-the-fly compression.\n",
    "\n",
    "This work is an extension of https://arxiv.org/pdf/1503.00638.pdf to floating point data to characterize different preprocessor..\n",
    "\n",
    "## Double precision, single precision and 32 bits data in 64 bits containers\n",
    "\n",
    "The easiest way to evaluate the effect of the mantissa-size on the compression is to compare double precision data and single precision data stored in 64 bits containers.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy\n",
    "import bitshuffle\n",
    "import lz4\n",
    "from lz4 import block, frame"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "shape = 1024,1024\n",
    "random_64 = numpy.random.random(shape)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Total size in float64: 8388608\n",
      "Total size of buffer: 8388608 100.0\n",
      "Total size in LZ4 : 8389139 100.01%\n",
      "Total size in bitshuffle/LZ4 : 7287095 86.87%\n"
     ]
    }
   ],
   "source": [
    "print(\"Total size in float64:\", random_64.nbytes)\n",
    "ref = len(random_64.tobytes())\n",
    "print(\"Total size of buffer:\", ref, 100.*ref/ref)\n",
    "lz4_64 = len(frame.compress(random_64.tobytes()))\n",
    "print(\"Total size in LZ4 : %i %.2f%%\"%(lz4_64, 100.*lz4_64/ref))\n",
    "blz4_64 = len(frame.compress(bitshuffle.bitshuffle(random_64).tobytes()))\n",
    "print(\"Total size in bitshuffle/LZ4 : %i %.2f%%\"%(blz4_64, 100.*blz4_64/ref))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Maximum error when considering 32bits floats 2.9802308953996715e-08\n"
     ]
    }
   ],
   "source": [
    "#Now work with 32bit data, and store them in a 64bit containers ...\n",
    "random_32 = random_64.astype(\"float32\")\n",
    "random_64_32 = random_32.astype(\"float64\")\n",
    "maximum_error = abs(random_64_32-random_64).max()\n",
    "print(\"Maximum error when considering 32bits floats\", maximum_error)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Total size in float32: 4194304\n",
      "Total size in float64: 8388608\n",
      "Total size of LZ4_32: 8388608 100.01%\n",
      "Total size in LZ4 64 truncated to 32: 6099766 72.71% 145.43%\n",
      "Total size in bitshuffle/LZ4_32 : 3460057 82.49% \n",
      "Total size in bitshuffle/LZ4_64_32 : 3482568 41.52% 83.03%\n"
     ]
    }
   ],
   "source": [
    "print(\"Total size in float32:\", random_32.nbytes)\n",
    "ref32 = len(random_32.tobytes())\n",
    "print(\"Total size in float64:\", random_64_32.nbytes)\n",
    "lz4_32 = len(frame.compress(random_32.tobytes()))\n",
    "print(\"Total size of LZ4_32: %i %.2f%%\"%(ref, 100.*lz4_32/ref32))\n",
    "lz4_64_32 = len(frame.compress(random_64_32.tobytes()))\n",
    "print(\"Total size in LZ4 64 truncated to 32: %i %.2f%% %.2f%%\"%(lz4_64_32, 100.*lz4_64_32/ref, 100.*lz4_64_32/ref32))\n",
    "blz4_32 = len(frame.compress(bitshuffle.bitshuffle(random_32).tobytes()))\n",
    "print(\"Total size in bitshuffle/LZ4_32 : %i %.2f%% \"%\n",
    "      (blz4_32, 100.*blz4_32/ref32))\n",
    "blz4_64_32 = len(frame.compress(bitshuffle.bitshuffle(random_64_32).tobytes()))\n",
    "print(\"Total size in bitshuffle/LZ4_64_32 : %i %.2f%% %.2f%%\"%\n",
    "      ( blz4_64_32, 100.*blz4_64_32/ref, 100.*blz4_64_32/ref32))\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "## Analysis\n",
    "\n",
    "* Floating point data compresses very bady, if at all (here we used white noise which is the worse case !)\n",
    "* Bitshuffle preprocessor helps in reformating the floating point data to compress togeather exponant and mantissa. All data being in the range 0-1 the exponant are the same but the mantissa contains white noise. So the theoretical maximum compression rate (output size/inputsize) should be 81% for double and 71% for single. There is no compression without preprocessing.\n",
    "* Storing limited precision float in larger containers compresses to the same size as when working with smaller containers (when using bitshuffle/LZ4). \n",
    "\n",
    "So it looks interesting to perform all calculation in double precision and to limit the size of the mantissa just prior to saving as suggested in the publication, even when working with floating point data.\n",
    "\n",
    "## Reducing the precision of the mantissa\n",
    "\n",
    "In computer, floating point are represented according to IEEE754 specification (https://en.wikipedia.org/wiki/IEEE_754), so the mantissa is stored in the least significant bits of the structure, and the heading \"1\" is omitted. The idea of the publication was to round the value so that the n-least significant bits are set to zeros to improve compression. \n",
    "\n",
    "We will treat this mantissa as an integrer number and use the same strategy, Kiyo Masui published in:\n",
    "https://gist.github.com/kiyo-masui/b61c7fa4f11fca453bdd"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0b100000 0b11111\n",
      "12345 0b11000000111001\n",
      "False 0b110010\n",
      "12352 0b11000001000000\n"
     ]
    }
   ],
   "source": [
    "# let n be the minimum number of bits set to 0\n",
    "n = 5\n",
    "gran = (1<<n)\n",
    "bitmask = gran - 1\n",
    "print(bin(gran), bin(bitmask))\n",
    "val = 12345\n",
    "print(val, bin(val))\n",
    "tie = ((val & bitmask) << 1) == gran\n",
    "print(tie, bin(((val & bitmask) << 1)))\n",
    "val_t = (val - (gran >> 1)) | bitmask\n",
    "val_t += 1\n",
    "val_t -= (gran >> 1) == 0\n",
    "val_t -= val_t & (tie * gran)\n",
    "print(val_t, bin(val_t))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "def round_int(val, nbits=5):\n",
    "    \"Round an integer so that the last nbits are set to zero\"\n",
    "    gran = (1<<nbits)\n",
    "    bitmask = gran - 1\n",
    "    tie = ((val & bitmask) << 1) == gran\n",
    "    val_t = (val - (gran >> 1)) | bitmask\n",
    "    val_t += 1\n",
    "    val_t -= (gran >> 1) == 0\n",
    "    val_t -= val_t & (tie * gran)\n",
    "    return val_t\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "#test round_int\n",
    "n = 8\n",
    "gran = (1<<n)\n",
    "bitmask = gran - 1\n",
    "for i in range(-65000, 65000):\n",
    "    j = round_int(i, n)\n",
    "    delta = abs(j-i)\n",
    "    if (delta>gran/2) or (j&bitmask):\n",
    "        print(i,j, delta, j&bitmask == 0)\n",
    "    \n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "def round_float(val, nbits=5):\n",
    "    mlen = None\n",
    "    try:\n",
    "        if val.dtype==numpy.float32:\n",
    "            mlen = 23\n",
    "        elif val.dtype == numpy.float64:\n",
    "            mlen = 52\n",
    "    except:\n",
    "        pass\n",
    "    if mlen is None:\n",
    "        mlen = 52\n",
    "    val = numpy.float64(val) \n",
    "    intval = int(val.view(int))\n",
    "    mask = 1<<(mlen)\n",
    "    bitmask = mask - 1\n",
    "    bigint = intval & bitmask\n",
    "    \n",
    "    bigint |= mask #set the bit at mantissa-length\n",
    "    #print(bin(bigint))\n",
    "    rndint = round_int(bigint, mlen-nbits)\n",
    "    #print(bin(rndint))\n",
    "    if rndint & (1<<(mlen+1)):\n",
    "        #print(val, bin(bigint), bin(rndint))\n",
    "        #Handle the rare case where offset of exponent is needed\n",
    "        bitmask = (1<<63) - mask\n",
    "        expo = (intval & bitmask)\n",
    "        expo += 1<<mlen\n",
    "        rndint = rndint>>1 #drop the right-handside bit\n",
    "        expo|= intval & (1<<63) #copy the sign which is left-hand bit\n",
    "        intval = expo\n",
    "        #print(val, bin(intval))\n",
    "    rndint &= ~mask #clear the the bit at mantissa-length\n",
    "    bitmask = (1<<64) - (1<<mlen)\n",
    "    rndint |= intval & bitmask\n",
    "    return numpy.uint64(rndint).view(\"float64\")\n",
    "    #return rndint"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "3.141592653589793\n",
      "Remaining Mantissa size | Value | Mesured error\n",
      "52 3.141592653589793 0.0\n",
      "51 3.141592653589793 0.0\n",
      "50 3.141592653589793 0.0\n",
      "49 3.141592653589793 0.0\n",
      "48 3.1415926535897967 -3.552713678800501e-15\n",
      "47 3.1415926535897967 -3.552713678800501e-15\n",
      "46 3.1415926535897825 1.0658141036401503e-14\n",
      "45 3.1415926535897825 1.0658141036401503e-14\n",
      "44 3.1415926535897825 1.0658141036401503e-14\n",
      "43 3.141592653589896 -1.0302869668521453e-13\n",
      "42 3.1415926535896688 1.2434497875801753e-13\n",
      "41 3.1415926535901235 -3.304023721284466e-13\n",
      "40 3.1415926535901235 -3.304023721284466e-13\n",
      "39 3.1415926535883045 1.4885870314174099e-12\n",
      "38 3.1415926535919425 -2.149391775674303e-12\n",
      "37 3.1415926535846666 5.126565838509123e-12\n",
      "36 3.1415926535846666 5.126565838509123e-12\n",
      "35 3.1415926535846666 5.126565838509123e-12\n",
      "34 3.1415926535846666 5.126565838509123e-12\n",
      "33 3.141592653701082 -1.1128875598842569e-10\n",
      "32 3.1415926534682512 1.2154188766544394e-10\n",
      "31 3.1415926534682512 1.2154188766544394e-10\n",
      "30 3.1415926534682512 1.2154188766544394e-10\n",
      "29 3.1415926553308964 -1.741103261565513e-09\n",
      "28 3.141592651605606 1.984187036896401e-09\n",
      "27 3.141592651605606 1.984187036896401e-09\n",
      "26 3.141592651605606 1.984187036896401e-09\n",
      "25 3.1415926814079285 -2.781813535079891e-08\n",
      "24 3.1415926218032837 3.1786509424591713e-08\n",
      "23 3.1415927410125732 -8.742278012618954e-08\n",
      "22 3.141592502593994 1.5099579897537296e-07\n",
      "21 3.1415929794311523 -3.2584135922775204e-07\n",
      "20 3.141592025756836 6.27832957178498e-07\n",
      "19 3.1415939331054688 -1.279515675634002e-06\n",
      "18 3.1415939331054688 -1.279515675634002e-06\n",
      "17 3.1415863037109375 6.349878855615998e-06\n",
      "16 3.1416015625 -8.908910206884002e-06\n",
      "15 3.1416015625 -8.908910206884002e-06\n",
      "14 3.1416015625 -8.908910206884002e-06\n",
      "13 3.1416015625 -8.908910206884002e-06\n",
      "12 3.1416015625 -8.908910206884002e-06\n",
      "11 3.1416015625 -8.908910206884002e-06\n",
      "10 3.140625 0.000967653589793116\n",
      "9 3.140625 0.000967653589793116\n",
      "8 3.140625 0.000967653589793116\n",
      "7 3.140625 0.000967653589793116\n",
      "6 3.15625 -0.014657346410206884\n",
      "5 3.125 0.016592653589793116\n",
      "4 3.125 0.016592653589793116\n",
      "3 3.25 -0.10840734641020688\n",
      "2 3.0 0.14159265358979312\n",
      "1 3.0 0.14159265358979312\n"
     ]
    }
   ],
   "source": [
    "pi = numpy.pi\n",
    "print(pi)\n",
    "print(\"Remaining Mantissa size | Value | Mesured error\")\n",
    "for i in range(52, 0, -1):\n",
    "    print( i, round_float(pi, i),pi - round_float(pi, i))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "3.75 µs ± 90.3 ns per loop (mean ± std. dev. of 7 runs, 100000 loops each)\n"
     ]
    }
   ],
   "source": [
    "%timeit round_float(pi, 20)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Analysis:\n",
    "We are now able to round floating point values to a given number of digits in the mantissa. Let's see how it behaves when compressing the data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "def round_array(ary, precision):\n",
    "    res = numpy.empty_like(ary)\n",
    "    flat_ary = ary.ravel()\n",
    "    flat_res = res.ravel()\n",
    "    for i in range(ary.size):\n",
    "        flat_res[i] = round_float(flat_ary[i], precision)\n",
    "    return res    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "CPU times: user 3.94 s, sys: 0 ns, total: 3.94 s\n",
      "Wall time: 3.94 s\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "7.957555701794154e-08"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "%time abs(round_array(random_64, 20)- random_64).mean()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Make it useable\n",
    "\n",
    "Let's port this to Cython to make it a bit faster"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [],
   "source": [
    "%load_ext Cython"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%cython\n",
    "\n",
    "import numpy\n",
    "cdef union float_and_int:\n",
    "    long long int   integer\n",
    "    double floating\n",
    "    \n",
    "cdef inline _round_int(val, nbits=5):\n",
    "    \"Round an integer so that the last nbits are set to zero\"\n",
    "#     cdef: \n",
    "#         working_t gran, bitmask, tie, val_t\n",
    "    gran = (1<<nbits)\n",
    "    bitmask = gran - 1\n",
    "    tie = ((val & bitmask) << 1) == gran\n",
    "    val_t = (val - (gran >> 1)) | bitmask\n",
    "    val_t += 1\n",
    "    val_t -= (gran >> 1) == 0\n",
    "    val_t -= val_t & (tie * gran)\n",
    "    return val_t\n",
    "\n",
    "def cround_int(val, nbits=5):\n",
    "    \"Round an integer so that the last nbits are set to zero\"\n",
    "    return _round_int(val, nbits)\n",
    "\n",
    "cdef inline double _round_float(double val, int nbits=5):\n",
    "    cdef:\n",
    "        float_and_int inp, out\n",
    "        int mlen = 52\n",
    "    inp.floating = val\n",
    "    intval = int(inp.integer)\n",
    "    mask = 1<<(mlen)\n",
    "    bitmask = mask - 1\n",
    "    bigint = intval & bitmask\n",
    "    \n",
    "    bigint |= mask #set the bit at mantissa-length\n",
    "    rndint = _round_int(bigint, mlen-nbits)\n",
    "    if rndint & (1<<(mlen+1)):\n",
    "        #Handle the rare case where offset of exponent is needed\n",
    "        bitmask = (1<<63) - mask\n",
    "        expo = (intval & bitmask)\n",
    "        expo += 1<<mlen\n",
    "        rndint = rndint>>1 #drop the right-handside bit\n",
    "        expo|= intval & (1<<63) #copy the sign which is left-hand bit\n",
    "        intval = expo\n",
    "    rndint &= ~mask #clear the the bit at mantissa-length\n",
    "    bitmask = (1<<64) - (1<<mlen)\n",
    "    rndint |= intval & bitmask\n",
    "    out.integer = rndint\n",
    "    return out.floating\n",
    "        \n",
    "    \n",
    "def cround_float(val, nbits=5):\n",
    "    \"Round the floating point value with nbits remaining in the mantissa\"\n",
    "    return _round_float(val, nbits)\n",
    "\n",
    "\n",
    "def cround_array(ary, precision):\n",
    "    cdef:\n",
    "        int i\n",
    "        double[:] flat_ary, flat_res\n",
    "    dary = numpy.ascontiguousarray(ary, dtype=numpy.float64)\n",
    "    res = numpy.zeros(ary.shape, dtype=numpy.float64)\n",
    "    flat_ary = dary.ravel()\n",
    "    flat_res = res.ravel()\n",
    "    for i in range(ary.size):\n",
    "        flat_res[i] = _round_float(flat_ary[i], precision)\n",
    "    return res    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "265 ns ± 8.39 ns per loop (mean ± std. dev. of 7 runs, 1000000 loops each)\n"
     ]
    }
   ],
   "source": [
    "#test round_int\n",
    "n = 8\n",
    "gran = (1<<n)\n",
    "bitmask = gran - 1\n",
    "for i in range(0, 65000):\n",
    "    j = cround_int(i, n)\n",
    "    delta = abs(j-i)\n",
    "    if (delta>gran/2) or (j&bitmask):\n",
    "        print(i,j, delta, j&bitmask == 0)\n",
    "\n",
    "%timeit cround_int(i, n)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "3.141592653589793\n",
      "Remaining Mantissa size | Value | Mesured error\n",
      "52 3.141592653589793 0.0\n",
      "51 3.141592653589793 0.0\n",
      "50 3.141592653589793 0.0\n",
      "49 3.141592653589793 0.0\n",
      "48 3.1415926535897967 -3.552713678800501e-15\n",
      "47 3.1415926535897967 -3.552713678800501e-15\n",
      "46 3.1415926535897825 1.0658141036401503e-14\n",
      "45 3.1415926535897825 1.0658141036401503e-14\n",
      "44 3.1415926535897825 1.0658141036401503e-14\n",
      "43 3.141592653589896 -1.0302869668521453e-13\n",
      "42 3.1415926535896688 1.2434497875801753e-13\n",
      "41 3.1415926535901235 -3.304023721284466e-13\n",
      "40 3.1415926535901235 -3.304023721284466e-13\n",
      "39 3.1415926535883045 1.4885870314174099e-12\n",
      "38 3.1415926535919425 -2.149391775674303e-12\n",
      "37 3.1415926535846666 5.126565838509123e-12\n",
      "36 3.1415926535846666 5.126565838509123e-12\n",
      "35 3.1415926535846666 5.126565838509123e-12\n",
      "34 3.1415926535846666 5.126565838509123e-12\n",
      "33 3.141592653701082 -1.1128875598842569e-10\n",
      "32 3.1415926534682512 1.2154188766544394e-10\n",
      "31 3.1415926534682512 1.2154188766544394e-10\n",
      "30 3.1415926534682512 1.2154188766544394e-10\n",
      "29 3.1415926553308964 -1.741103261565513e-09\n",
      "28 3.141592651605606 1.984187036896401e-09\n",
      "27 3.141592651605606 1.984187036896401e-09\n",
      "26 3.141592651605606 1.984187036896401e-09\n",
      "25 3.1415926814079285 -2.781813535079891e-08\n",
      "24 3.1415926218032837 3.1786509424591713e-08\n",
      "23 3.1415927410125732 -8.742278012618954e-08\n",
      "22 3.141592502593994 1.5099579897537296e-07\n",
      "21 3.1415929794311523 -3.2584135922775204e-07\n",
      "20 3.141592025756836 6.27832957178498e-07\n",
      "19 3.1415939331054688 -1.279515675634002e-06\n",
      "18 3.1415939331054688 -1.279515675634002e-06\n",
      "17 3.1415863037109375 6.349878855615998e-06\n",
      "16 3.1416015625 -8.908910206884002e-06\n",
      "15 3.1416015625 -8.908910206884002e-06\n",
      "14 3.1416015625 -8.908910206884002e-06\n",
      "13 3.1416015625 -8.908910206884002e-06\n",
      "12 3.1416015625 -8.908910206884002e-06\n",
      "11 3.1416015625 -8.908910206884002e-06\n",
      "10 3.140625 0.000967653589793116\n",
      "9 3.140625 0.000967653589793116\n",
      "8 3.140625 0.000967653589793116\n",
      "7 3.140625 0.000967653589793116\n",
      "6 3.15625 -0.014657346410206884\n",
      "5 3.125 0.016592653589793116\n",
      "4 3.125 0.016592653589793116\n",
      "3 3.25 -0.10840734641020688\n",
      "2 3.0 0.14159265358979312\n",
      "1 3.0 0.14159265358979312\n",
      "661 ns ± 20.4 ns per loop (mean ± std. dev. of 7 runs, 1000000 loops each)\n"
     ]
    }
   ],
   "source": [
    "pi = numpy.pi\n",
    "print(pi)\n",
    "print(\"Remaining Mantissa size | Value | Mesured error\")\n",
    "for i in range(52, 0, -1):\n",
    "    print( i, cround_float(pi, i),pi - cround_float(pi, i))\n",
    "%timeit cround_float(pi, 20)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "CPU times: user 576 ms, sys: 12 ms, total: 588 ms\n",
      "Wall time: 582 ms\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "2.9802308953996715e-08"
      ]
     },
     "execution_count": 18,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "%time abs(cround_array(random_64, 23)- random_64).max()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Compressed size as a function on the mantissa length\n",
    "\n",
    "We have now the tools to benchmark the compressed size as function of the mantissa size of the float64."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [],
   "source": [
    "sizes = [len(frame.compress(bitshuffle.bitshuffle(cround_array(random_64, i)).tobytes())) for i in range(53)]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[5.3472161293029785, 7.01749324798584, 8.63410234451294, 10.21111011505127, 11.792588233947754, 13.352382183074951, 14.932501316070557, 16.48489236831665, 18.059325218200684, 19.612765312194824, 21.192705631256104, 22.74419069290161, 24.32166337966919, 25.873327255249023, 27.442800998687744, 29.002833366394043, 30.56471347808838, 32.130515575408936, 33.69598388671875, 35.25642156600952, 36.82820796966553, 38.3844256401062, 39.947545528411865, 41.515445709228516, 43.084120750427246, 44.64414119720459, 46.21344804763794, 47.77116775512695, 49.33271408081055, 50.906574726104736, 52.45786905288696, 54.038870334625244, 55.59808015823364, 57.15765953063965, 58.72153043746948, 60.29253005981445, 61.835408210754395, 63.41524124145508, 64.97175693511963, 66.55751466751099, 68.1148886680603, 69.67103481292725, 71.21807336807251, 72.80173301696777, 74.35799837112427, 75.9243369102478, 77.46745347976685, 79.04058694839478, 80.62218427658081, 82.20208883285522, 83.7560772895813, 85.3075385093689, 86.86894178390503]\n"
     ]
    }
   ],
   "source": [
    "compression_rate = [100*i/ref for i in sizes]\n",
    "print(compression_rate)\n",
    "theoretical = [(64-(52-i))*100/64 for i in range(53)]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "from matplotlib.pyplot import subplots"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/mntdirect/_scisoft/users/jupyter/jupy35/lib/python3.5/site-packages/matplotlib/figure.py:459: UserWarning: matplotlib is currently using a non-GUI backend, so cannot show the figure\n",
      "  \"matplotlib is currently using a non-GUI backend, \"\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYgAAAEKCAYAAAAIO8L1AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzs3Xd4FOX2wPHvSSGhhBoIHUKvoYMCoRdFFBWwXERRlOu1YMMrFgQVEAQVwZ8FywV7r/QSOiKEXhJIgAChE0ggvez7+2MWjBhgA9lMNjmf58mzO7MzO2eGZc/OnHnfV4wxKKWUUhfzsjsApZRSBZMmCKWUUjnSBKGUUipHmiCUUkrlSBOEUkqpHGmCUEoplSNNEEoppXKkCUIppVSONEEopZTKkY/dAVyLwMBAU7t2bbvDUEopj7Jx48ZTxpiKV1rOoxNE7dq1CQ8PtzsMpZTyKCJywJXl9BKTUkqpHGmCUEoplSNNEEoppXLk0TWInGRkZBAbG0tqaqrdoRR6/v7+VK9eHV9fX7tDUUq5QaFLELGxsQQEBFC7dm1ExO5wCi1jDHFxccTGxhIcHGx3OEopN3DbJSYR+VRETojIjmzzyovIYhGJcj6Wc84XEZkuItEisk1EWl/tdlNTU6lQoYImBzcTESpUqKBnakoVYu6sQcwCbrho3mhgqTGmPrDUOQ1wI1Df+TcCeP9aNqzJIX/ocVaqcHNbgjDGrAROXzR7ADDb+Xw2cGu2+Z8ZyzqgrIhUcVdsSinlsVITYMk4OBPj9k3l911MQcaYo87nx4Ag5/NqwKFsy8U65/2DiIwQkXARCT958qT7Ii2iunXrpo0PlSqIsjJg/UcwvRWsfhuil7h9k7YVqY0xRkTMVaw3E5gJ0LZt21yvXxRlZmbi41Po7kdQqmgwBnbPh8UvQ1wU1A6FPuOhaku3bzq/zyCOn7905Hw84Zx/GKiRbbnqznkeKSYmhkaNGjFs2DAaNGjAkCFDWLJkCZ06daJ+/fqsX7+epKQkHnjgAdq3b0+rVq349ddfL6wbGhpK69atad26NWvXrgXg6NGjdOnShZYtW9KsWTNWrVoFQKlSpS5s94cffmDYsGEADBs2jIcffpgOHTrw3//+95LbS0lJ4a677qJx48bcdtttpKSk5OORUkpd1pEtMPtm+OZuwMBdX8N9v+dLcoD8P4P4DbgPmOR8/DXb/MdE5BugA5CQ7VLUVXvl953sOnL2Wt/mb5pULc3Ym5tecbno6Gi+//57Pv30U9q1a8dXX33F6tWr+e2335g4cSJNmjShR48efPrpp8THx9O+fXt69epFpUqVWLx4Mf7+/kRFRXH33XcTHh7OV199Rd++fXnxxRfJysoiOTn5ijHExsaydu1avL29eeGFF3Lc3ocffkiJEiWIiIhg27ZttG591TeQKaXySkIsLH0Ntn0DxcvDjVOg7f3gnb9tjtyWIETka6AbECgiscBYrMTwnYgMBw4AdzgXnwf0A6KBZOB+d8WVX4KDg2nevDkATZs2pWfPnogIzZs3JyYmhtjYWH777TemTp0KWLfnHjx4kKpVq/LYY4+xZcsWvL292bNnDwDt2rXjgQceICMjg1tvvZWWLa/8C2Lw4MF4e3sDsGjRohy3t3LlSkaOHAlASEgIISEheX4slFIuSjsHq6fBH++CcUDHkRD6DBQva0s4bksQxpi7L/FSzxyWNcCjeR2DK7/03cXPz+/Ccy8vrwvTXl5eZGZm4u3tzY8//kjDhg3/tt64ceMICgpi69atOBwO/P39AejSpQsrV65k7ty5DBs2jKeffpp77733b7eaXtwmoWTJkheeG2Ny3J5SqgDIyoTNn8OyiZB0ApoNgp4vQ7latoalfTHZpG/fvsyYMQMrN8LmzZsBSEhIoEqVKnh5efH555+TlZUFwIEDBwgKCuKhhx7iwQcfZNOmTQAEBQURERGBw+Hg559/zvX2unTpwldffQXAjh072LZtm3t2WCn1T8ZA1GL4oDPMeRIq1IUHl8KgT2xPDqAJwjZjxowhIyODkJAQmjZtypgxYwB45JFHmD17Ni1atCAyMvLCWcDy5ctp0aIFrVq14ttvv+WJJ54AYNKkSfTv35+OHTtSpcqlm45canv/+c9/SExMpHHjxrz88su0adPGzXuulALg2A74/Db4chBkpsIdn8P986F6W7sju0DO/6L0RG3btjUX37MfERFB48aNbYqo6NHjrVQunT0Ky8bD5i/Bvwx0Gw1th4NPsXwLQUQ2GmOumIn05nillMoP6Umwdgasecdq9Hb9o9BlFBQvZ3dkl6QJQiml3MmRBVu+grDxkHgMmgyAXuOgfB27I7siTRBKKeUue5fBopfg+A6o1hbumA01r7M7KpdpglBKqbx2IhIWj4GoRVC2Jgz6FJreDh7WA7ImCKWUyiuJJ6y2DJtmQ7EA6P0atB8Bvv52R3ZVNEEopdS1Sk+Gdf9ntYLOTIV2D0HX56BkBbsjuyaaIJRS6mo5HLD9O1j6Kpw9DI36Q69XILCe3ZHlCU0QBVRWVtaFfpSUUgXQ/lWw6EU4uhWqtITbZ0LtznZHlae0JXUemzJlCtOnTwfgqaeeokePHgCEhYUxZMgQAL7++muaN29Os2bNeO655y6sW6pUKZ555hlatGjBH3/8wejRo2nSpAkhISGMGjUKgJMnTzJw4EDatWtHu3btWLNmzT9iyMrKYtSoUTRr1oyQkBBmzJgBwKuvvkq7du1o1qwZI0aMuNDtxvTp0y9s56677gK4ZPfgShV5p6Lg67thdn9IioPbZsJDywpdcoDCfgYxfzQc256371m5Odw46ZIvh4aG8uabbzJy5EjCw8NJS0sjIyODVatW0aVLF44cOcJzzz3Hxo0bKVeuHH369OGXX37h1ltvJSkpiQ4dOvDmm28SFxfH8OHDiYyMRESIj48H4IknnuCpp56ic+fOHDx4kL59+xIREfG3GGbOnElMTAxbtmzBx8eH06etkV8fe+wxXn75ZQCGDh3KnDlzuPnmm5k0aRL79+/Hz8/vwnYmTJiQY/fg2TsAVKpISYqD5a9D+KfgW8LqTO+6R8C3uN2RuY2eQeSxNm3asHHjRs6ePYufnx/XX3894eHhrFq1itDQUDZs2EC3bt2oWLEiPj4+DBkyhJUrVwLg7e3NwIEDAShTpgz+/v4MHz6cn376iRIlSgCwZMkSHnvsMVq2bMktt9zC2bNnSUxM/FsMS5Ys4d///veFUeTKly8PwLJly+jQoQPNmzcnLCyMnTt3AlY330OGDOGLL764sM6iRYuYNGkSLVu2pFu3bhe6B1eqyMlItYrP01tayaHNMBi52eqGuxAnByjsZxCX+aXvLr6+vgQHBzNr1iw6duxISEgIy5YtIzo6msaNGxMVFXXJdf39/S/UHXx8fFi/fj1Lly7lhx9+4N133yUsLAyHw8G6desudAPuqtTUVB555BHCw8OpUaMG48aNu9A9+Ny5c1m5ciW///47EyZMYPv27do9uFLGwI4fYckrkHAQ6veFPq9BxaLzf0LPINwgNDSUqVOn0qVLF0JDQ/nggw9o1aoVIkL79u1ZsWIFp06dIisri6+//pquXbv+4z0SExNJSEigX79+vP3222zduhWAPn36XKgpAGzZsuUf6/bu3ZsPP/yQzMxMAE6fPn0hGQQGBpKYmMgPP/wAgMPh4NChQ3Tv3p3JkyeTkJBAYmLiJbsHV6pIOPAHfNwTfhwOxcvAvb/BkO+KVHKAwn4GYZPQ0FAmTJjA9ddfT8mSJfH39yc0NBSAKlWqMGnSJLp3744xhptuuokBAwb84z3OnTvHgAEDSE1NxRjDW2+9BVgF5UcffZSQkBAyMzPp0qULH3zwwd/WffDBB9mzZw8hISH4+vry0EMP8dhjj/HQQw/RrFkzKleuTLt27QCroH3PPfeQkJCAMYaRI0dStmxZxowZw5NPPklISAgOh4Pg4GDmzJnj5iOnlM3i9sKSsRDxOwRUhVvfh5C7wKto/pbW7r7VNdHjrQqF5NOwcgqs/wi8i0Hnp6zeVouVsDsyt9DuvpVS6koy06yksPINazzoVkOh+4sQEGR3ZAWCJgilVNFjDOz6FZaMgzP7oV4vq9+koCZ2R1agFMoEYYxBPKzXRE/kyZcnVREWGw4LX4RD66BSU7jnJ6jX0+6ocuVgXDJVyvrj6+3e2kihq7z4+/sTFxenX15uZowhLi4u17fbKmWbMzHw/f3W3Uln9sPN0+HhVR6VHBKSM3htzi56vrWcr/50f7ukQncGUb16dWJjYzl58qTdoRR6/v7+VK9e3e4wlLq8lHhY9Sb8+QGIt9XLaseR4FfK7shclp7p4PN1B5i+NIqzqRnc0aYGNzSr7PbtFroEcb6hmlKqiMvKsFo+L58EKWeg5RDo8SKUrmp3ZC4zxrBw5zEmzY8kJi6Z0PqBvNCvMY2rlM6X7Re6BKGUKuKMgci5sPhlOL0XgrtAnwlQJcTuyHJl88EzTJwXwYaYMzQIKsWs+9vRtUHFfK2vaoJQShUehzdZY0AfWAOBDeFf30H9Ph411OfBuGTeWBjJnG1HCSzlx8TbmnNH2+r4uLkgnRNNEEopzxd/yBq0Z/t3UCIQbnoLWt8H3p7zFRefnM67YdHM/iMGHy8vRvasz7+71KGkn3374DlHTymlLpZ6Fla/BX+8Z013ftpqBe2fP9fo80JaZhaf/3GAGWHRFwrQT/dpQFBp++8Q1AShlPI8WZmwaRYsex2ST0HIndBjDJStYXdkLjPGMHf7USYviOTQ6ZR8L0C7QhOEUspzGANRi2DRGDi1G2p1gj7fQ7XWdkeWK+Exp5kwL4LNB+NpVDmAzx5oT5cGFe0O6x80QSilPMPRbdYY0PtXQvm6cOeX0OgmjypAx5xKYvKCSObvOEalAD/eGBjCwDbV8fYqmPugCUIpVbCdPQJh42HLV1C8HNwwGdo+AD7F7I7MZQnJGcwIi2L2HzH4envxdO8GPBgaTIliBfsr2JboROQp4EHAANuB+4EqwDdABWAjMNQYk25HfEqpAiAtEda8A2tngMmCjo9B6CgoXtbuyFyWkeXgy3UHmLY0ioQUqwD9TJ8GVCoABWhX5HuCEJFqwEigiTEmRUS+A+4C+gFvG2O+EZEPgOHA+/kdn1LKZo4s2Pw5hE2ApBPQ9HboNRbK1bY7MpcZYwiLPMGEeRHsO5lEx7oVeOmmJjSpWnAK0K6w6/zGByguIhlACeAo0AP4l/P12cA4NEEoVbREL7EK0Cd2QY0OcNeXUKO93VHlSsTRs0yYG8Hq6FPUCSzJx/e2pWfjSh7Zw3S+JwhjzGERmQocBFKARViXlOKNMZnOxWKBavkdm1LKJsd3Wi2g94ZZZwqDZ0OTAR5VgD5xLpW3Fu3hu/BDBPj7MvbmJtxzXS23d8ntTnZcYioHDACCgXjge+CGXKw/AhgBULNmTXeEqJTKL+eOwbIJsPkL8Auw+kxq/xD4+NkdmctSM7L4ZPV+3lsWTVqmg/s7BfN4j3qULeE5RfRLseMSUy9gvzHmJICI/AR0AsqKiI/zLKI6cDinlY0xM4GZYI1JnT8hK6XyVHoSrH3XKkJnpUOHh6HLs1CivN2RucwYw29bj/DGgt0cjk+hT5Mgnu/XmODAknaHlmfsSBAHgetEpATWJaaeQDiwDBiEdSfTfcCvNsSmlHInhwO2fg1hr8G5o9D4Fug1DirUtTuyXNl44DSvzYlgy6F4mlYtzZTBIXSsG2h3WHnOjhrEnyLyA7AJyAQ2Y50RzAW+EZHxznmf5HdsSik32rfCauh2bDtUawOD/ge1rrc7qlw5GJfM5AWRzN1+lKDSfkwd3ILbW1XDq4A2dLtWttzFZIwZC4y9aPY+wLNuV1BKXdnJ3dbYDHsWQJmaMPAT69ZVL88p3iakZPB/y6KZtSYGby/hyV71GdGlToFv6HatCvfeKaXsk3gSlr8OG2dBsZLQ6xWr1uDrGY3E4K+eVt9dFk1CSgYDW1dnVJ+GVC7jOftwLTRBKKXyVkYKrHsPVr0NGcnQbjh0HQ0lK9gdmcty6ml19I2NaFq1jN2h5StNEEqpvOFwwI4frIF7Eg5Bw37Q+1UIrG93ZLmyIeY0E+ZaBehGlQOY/UB7uhbAnlbzgyYIpdS1i1ljFaCPbIYqLeDW9yE41O6ocmX/qSQmzY9g4c7jBJX2441BIQxsXXB7Ws0PmiCUUlfvVDQsGQuRc6B0NbjtQ2h+h0cVoOOT05m+NJrP11k9rT7TuwEPhtaheDFvu0OznSYIpVTuJcXBiskQ/gn4+EOPl+C6R6FYCbsjc1l6poPP/ohhRlg051IzuLNdDZ7q3YBKAUWjAO0KTRBKKddlpsGfH8LKqZB+DlrfB91fgFKV7I7MZcYY5u84xuQFkRyIS6ZLg4q82K8xDSsH2B1agaMJQil1ZcbAzp9gyTiIPwj1+0Dv16BSI7sjy5WNB84wYe4uNh2Mp2FQALPub0e3hp6T3PKbSwlCRIKAds7J9caYE+4LSSlVoBz80ypAx26AoGYw9Beo293uqHLlYFwykxdGMnfbUSoG+DHp9uYMblujSBegXXHFBCEidwBTgOWAADNE5FljzA9ujk0pZafT+2DJK7DrFyhVGW55F1r+C7w8p3ibkJLBu2FRzF57AG8v4YmeVgvokn568cQVrhylF4F2588aRKQisATQBKFUYZRyxqox/PkhePtCt+eh4+NWa2gPkZHl4Ks/DzJtyR7iUzIY3KY6z/RpSJCHDPVZULiSILwuuqQUB3jOPWxKKddkpsOGj627k1IToNUQ6P4SlK5id2QuM8awbPcJJsyNYO/JJK6vU4GX+jcuci2g84orCWKBiCwEvnZO3wnMc19ISql8ZQxE/AaLx8KZ/VCnmzVwT+VmdkeWK5HHrKE+V0WdIjiwJB/d25ZeHjrUZ0FxxQRhjHlWRAZiDeoDMNMY87N7w1JK5YvYjVYB+uAfULERDPkB6vXyqKE+jyWk8tbi3Xy/MZYAPx/G9G/C0OtqUcxHL3RcK5cqNcaYH4Ef3RyLUiq/nDkAS1+BHT9CyYrQfxq0GgrenlO8TUzLZOaKvcxctY8sh+GBQjTUZ0FxyU+DiKw2xnQWkXNA9qE9BTDGmNJuj04plbdSE2DVm7DuA+ssIXQUdH7SGg/aQ2RmOfg2/BBvL47iVGIa/UOq8N++jahZwXNacXuKSyYIY0xn56PnfHKUUjnLyrDGZVj+OiTHQYu7occYKFPN7shcZowhLPIEr8+PJPpEIu1ql+Oje9vQqmY5u0MrtFxpB/G5MWboleYppQogY2D3fGtEt7goqB0KfcZD1ZZ2R5YrOw4nMHFeBGv3xhEcWJIP7mlD36ZBWoB2M1cuODbNPiEiPkAb94SjlMozR7bAopcgZhVUqA93fwMNbvCoAvSR+BSmLtrNz5sPU7a4L6/c0pR/daiJr7cWoPPD5WoQzwMvAMVF5Oz52UA6MDMfYlNKXY2EWFj6Gmz7BkpUgH5Toc0wq9GbhziXmsEHK/by8ar9GGBElzo82r0epf09Zx8Kg8vVIF4HXheR140xz+djTEqpq5F2Dla/DX/8n3VpqdOTEPo0+HtOI7HMLAffbDjEtCV7OJWYzi0tqvJs34bUKK8FaDu40g7ieREpB9QH/LPNX+nOwJRSLsrKhM2fwbKJkHQSmg2Cni9DuVp2R+ay8wXoifOsFtDta5fnk/sa06JGWbtDK9JcKVI/CDwBVAe2ANcBfwA93BuaUuqyjIHoJVad4WQk1Lwe7v4WqntWiXDnkQQmzP2rAP3h0Db0aaIF6ILAlSL1E1hdfa8zxnQXkUbARPeGpZS6rGPbrcSwbzmUrwN3fA6Nb/aoAvSJs6lMXWS1gNYCdMHkSoJINcakiggi4meMiRSRhm6PTCn1T2ePwrLxsPlLKF4WbpgEbYeDj+e0Hk7NyOLjVft4b/leMrIcPNg5mMd61KdMcS1AFzSuJIhYESkL/AIsFpEzwAH3hqWU+pv0JFgzHdZOtxq9Xf8odBkFxT2nkZgxht+2HuGNBbs5HJ9C36ZBPH9jY2oHek434kWNK0Xq25xPx4nIMqAMsMCtUSmlLI4s2PIVhI2HxGPQ9DboORbKB9sdWa6s33+aCfMi2HoonqZVSzN1cAuur1vB7rDUFVw2QYiIN7DTGNMIwBizIl+iUkrB3jBYNAaO74Dq7eDOz6FGe7ujypV9JxOZvCCShTuPU7m0P1MHt+C2VtV0qE8PcdkEYYzJEpHdIlLTGHMwv4JSqkg7EWEVoKOXQNlaMOh/1pmDBxWgTyelM31pFF+sO4Cfjxej+jRgeOc6FC/mOcOVKtdqEOWAnSKyHkg6P9MYc4vbolKqKEo8AcsmwKbPrN5V+4yH9iPAx8/uyFyWlpnF7LUxzAiLJiktk7vb1+TJXg2oGOA5+6D+4kqCGOP2KJQqytKTrdbPa6ZBZqqVFLo+ByXK2x2Zy4wxLNhxjNfnR3LwdDLdG1bkhX6NqR+knUF7MleK1Fp3UModHA7Y9i2EvQZnD0Oj/tDrFQisZ3dkubItNp7xcyJYH3OaBkGl+OyB9nRpUNHusFQe8Jzho5QqTPavsob6PLoVqraC2z+C2p2uvF4Bcjg+hakLrZ5WA0sVY+JtzbmjbXV8tKFboWFLgnC2q/gYaIY1Wt0DwG7gW6A2EAPcYYw5Y0d8SrnNqSjrzqQ986FMDbj9Y2g2ELw850s1ISWD95ZH8781MQA83LUuj3avS4D2tFro2HUG8Q6wwBgzSESKASWwuhZfaoyZJCKjgdHAczbFp1TeSjoFyydB+KfgW8Jqy3Ddf8C3uN2RuSw908EX6w4wIyyK+JQMbmtVjWf6NKRaWc/ZB5U7rnTW1wkYB9RyLn9+TOo6V7NBESkDdAGGYb1ROpAuIgOAbs7FZgPL0QShPF1GKvz5Pqx6y2oN3fZ+6DoaSnnONXpjDPN3HGPygkgOxCXTuV4go29sRLNqntONuLo6rpxBfAI8BWwEsvJgm8HASeB/ItLC+b5PAEHGmKPOZY4BQXmwLaXs4XDAzp9gySuQcBAa3Ai9X4GKntWN2eaDZ5gwN4LwA2doGBTArPvb0bVBRe1ptYhwJUEkGGPm5/E2WwOPG2P+FJF3sC4nXWCMMSJiclpZREYAIwBq1qyZh2EplUcO/GEVoA9vhMohMOBdqNPV7qhyJfZMMm8s2M1vW48QWMqPSbc3Z3DbGtoCuohxJUEsE5EpwE9A2vmZxphNV7nNWCDWGPOnc/oHrARxXESqGGOOikgV4EROKxtjZuIc8rRt27Y5JhGlbBG3F5aMhYjfIaAq3Po+hNzlUQXoc6kZvL98Lx+v3o8Aj3Wvx8Pd6lLKT294LIpc+Vfv4Hxsm22e4SoHDDLGHBORQyLS0BizG+gJ7HL+3QdMcj7+ejXvr1S+Sz4NK96ADR+Btx90f8nqbbWY5wyTmZnl4OsNh5i2eA9xSenc1qoaz/ZtSFUtQBdprjSU6+6G7T4OfOm8g2kfcD/gBXwnIsOxuhO/ww3bVSrvZKbB+pmwcoo1HnSrodD9RQjwnPKZMYZlu08wcV4k0ScSaV+7PJ8O06E+lcWVu5jKAGOx7jwCWAG8aoxJuNqNGmO28PczkvN6Xu17KpVvjIFdv8CScXAmBur1gt6vQVATuyPLlZ1HEpg4L4I10TrUp8qZK5eYPgV28Ncv+qHA/4Db3RWUUgXWofWw8EWIXQ+VmsI9P1oJwoMcS7CG+vxxUyxlivsy9uYmDOlQi2I+nlMrUfnDlQRR1xgzMNv0KyKyxV0BKVUgnYmxzhh2/gylguCWGdByCHh5TvfVSWmZfLhiLzNX7cPhwBrqs3t9ypTQFtAqZ64kiBQR6WyMWQ0XGs6luDcspQqIlDOw6k3480MQb6uX1Y4jwa+U3ZG5LMth+C78EG8u2sOpxDT6h1ThuRsaUaO85xTRlT1cSRD/AWY7axECnMbZClqpQisz3eoWY8UkSIm3zhZ6vAilq9odWa4s332C1+dFsvv4OdrUKsfMe9vQuqbnjGOt7OXKXUxbgBYiUto5fdbtUSllF2Mgcg4sHgun90JwV2vgniohdkeWKxFHzzJxXgSrok5Rs3wJ3h/SmhuaVdYCtMqVSyYIEbnHGPOFiDx90XwAjDFvuTk2pfLX4Y1WT6sH1kBgQ/jX91C/t0cN9Xn8bCpvLtrN9xtjKe3vy5j+TRh6nRag1dW53BlESeejDgmlCrf4g7D0Vdj+PZQIhJvegtb3gbfntB5OTs/kwxX7mLlyH5kOB8M7BfN4Dy1Aq2tzyf8BxpgPnY+v5F84SuWj1ARY/Tb88Z51lhD6DHR6EvxL2x2ZyxwOw0+bDzNlYSTHz6ZxU3OrAF2zghag1bVzpaHcG8B4rDuXFgAhwFPGmC/cHJtS7pGVARtnwfLXITnO6i+px0tQtobdkeXKun1xjJ+7ix2Hz9KiRlneG9KaNrU8ZxxrVfC5cg7dxxjzXxG5DWukt9uBlYAmCOVZjIE9C2Dxy3BqD9TqDH3HW0N+epCYU0m8Pj+ChTuPU7WMP+/c1ZKbQ6ripT2tqjzmSoI4v8xNwPfGmAS9E0J5nKNbrRbQMaugQj246yto2M+jCtBnktKZHhbFF+sO4Ovtxag+DXgwtA7+vp7TWE95FlcSxBwRicS6xPQfEakIpLo3LKXyyNkjsPQ12Po1FC8HN06xRnXz9pzibWpGFrPXxvDusmiS0jK5s11Nnupdn0oB/naHpgo5V9pBjHbWIRKMMVkikgQMcH9oSl2DtERY8w6snQEmCzo+bhWhi3tOL6UOh+H3bUd4Y8FuDsen0L1hRZ7v15gGQXpjocofrhSpBwMLnMnhJazR4MZjDQuqVMHiyILNn0PYBEg6Ac0GQs+xUK6W3ZHlyoaY04yfs4utsQk0qVKaNwaF0KleoN1hqSLGlUtMY4wx34tIZ6AXMAV4n78GElKqYIhaAovHwIldUKMD3P01VM+pV/mC60BcEpPmRzJ/xzEql/Zn6uAW3N6qmhaglS1cSRBZzsebgJnGmLkiMt6NMSmVO8d3wqKXYG8YlKsNg2dDkwEeVYBOSMng3bAoZq2NwcfLi6d7N+Ch0DpLnAo+AAAepUlEQVQUL6YFaGUfVxLEYRH5EOgNTBYRP6zR35Sy17ljsGwCbP4C/EpDnwnQ/iHw8bM7MpdlZDn4ct0B3lkaRXxKBoPbVOeZPg0JKq0FaGU/VxLEHcANwFRjTLyIVAGedW9YSl1GehKsfdcqQmelQ4eHocuzUMJzGokZY1i06ziT5key/1QSHetW4MWbGtO0ahm7Q1PqAlfuYkoWkRNAZyAKyHQ+KpW/HFmw9RsIew3OHYXGt0CvcVChrt2R5cq22HgmzI3gz/2nqVepFJ8Oa0v3hpW0p1VV4LhyF9NYrPGjG2INNeqL1Yq6k3tDUyqbfcutOsOx7VCtDQz6H9S63u6ocuVwfApTFkTyy5YjVChZjPG3NuOudjXw8dYrtqpgcuUS021AK2ATgDHmiIjojdgqf5zcbXXBHbUQytSEgZ9Yt6560K/tc6kZvL98L5+s3g/Ao93r8nDXugT4e05jPVU0uZIg0o0xRkQMgIiUvNIKSl2zxJOwfCJsnA3FSkHvV6H9v8HXc4q3GVkOvll/kGlLoohLSuf2VtV4pm9DqpUtbndoSrnElQTxnfMuprIi8hDwAPCRe8NSRVZGCqx7D1a9DZkp0G44dB0NJSvYHZnLjDEsjTjBxPkR7DuZxHV1yjOrXxOaV9cCtPIsrhSpp4pIb+AsVh3iZWPMYrdHpooWh8MasGfpq3A2FhreBL1fgcD6dkeWKzsOJzB+7i7W7TtNnYol+fjetvRsrAVo5ZkumyBExBtYYozpDmhSUO4Rs9rqafXoFqjSAm77AIJD7Y4qV44lpDJl4W5+2hxLuRLFeHVAU+5uXxNfLUArD3bZBOHsf8khImWMMQn5FZQqIk5FW2Mz7J4LpavBbTOh+WDw8pwv1aS0TD5cuY+ZK/ficMCILnV4tHs9SmsBWhUCrtQgEoHtIrIYSDo/0xgz0m1RqcItKQ5WTIbwT8DHH3qMgesfBV/PKd5mOQw/bDzEm4v2cOJcGv1DrKE+a5TXoT5V4eFKgvjJ+afUtclIhfUfwso3If0ctBkG3Z6HUpXsjixXVu45ycR5EUQeO0frmmV5/542tKlVzu6wlMpzrhSpZ4tIMaARYIDdxph0t0emCg9jYOdPsGQcxB+E+n2g92tQqZHdkeXK7mPnmDAvgpV7TlKzfAneG9KaG5tV1gK0KrRcaUndD/gQ2AsIECwi/zbGzHd3cKoQOPgnLHwBDodDUDMY+gvU7W53VLly4lwqby3aw3fhhwjw9+Wlmxoz9Ppa+PloT6uqcHPlEtNbQHdjTDSAiNQF5gKaINSlnd5nnTHs+hUCqsCA/4MWd4OX53yppqRn8dGqfXywYi8ZWQ6GdQxmZM96lC1RzO7QlMoXriSIc+eTg9M+4Jyb4lGeLvk0rJwK62eCdzHo9gJ0fAyKeU4DfIfD8PPmw0xZuJtjZ1O5oWllRt/YiNqBnrMPSuUFVxJEuIjMA77DqkEMBjaIyO0AxhgtYCvITIcNH8GKNyA1AVrdAz1egoDKdkeWK3/sjWPCvF3sOHyWFtXLMP3uVrQP9pxuxJXKS64kCH/gONDVOX0SKA7cjJUwripBOBvhhQOHjTH9RSQY+AaoAGwEhmox3AMYAxG/weKxcGY/1O1hFaArN7M7slzZezKR1+dFsiTiOFXL+DPtzpbc0qKqDvWpijRX7mK6303bfgKIAEo7pycDbxtjvhGRD4DhWGNfq4IqNtxqAX1oHVRsDEN+hPq97I4qV84kpfPO0ii+WHcAPx8vnu3bkOGdg/H39ZxaiVLu4spdTMHA40Dt7MsbY2652o2KSHWsMa4nAE+LdZ9gD+BfzkVmA+PQBFEwnTkAS1+BHT9CyUrQfxq0GgrerpyQFgxpmVnMXhvDjLBoktIyuat9TZ7q1YCKAZ4zXKlS7ubK/+hfgE+A3wFHHm13GvBf4Py4EhWAeGNMpnM6FqiWR9tSeSUlHla9CX9+AOJtDfPZ6Qnw85zhQYwxzNt+jMkLIjl4OpluDSvyQr/GNAjynH1QKr+4kiBSjTHT82qDItIfOGGM2Sgi3a5i/RHACICaNWvmVVjqcrIyIPx/sPx1SDlj3a7a4yUo41k5fOOBM0yYu4tNB+NpVDmAzx5oT5cGFe0OS6kCy5UE8Y5z2NFFQNr5mcaYTVe5zU7ALc4GeP5YNYh3sMab8HGeRVQHDue0sjFmJjAToG3btuYqY1CuMAZ2z7c61IuLgtqh0HeC1eOqBzkYl8zkhZHM3XaUigF+TB7YnEFtauCtBWilLsuVBNEcGIpVIzh/ick4p3PNGPM88DyA8wxilDFmiIh8DwzCupPpPuDXq3l/lUeObLaG+oxZBYEN4O5voMENHjXUZ0JyBu8ui2L22gN4ecHInvX5d5c6lPTznFqJUnZy5X/KYKBOPtxy+hzwjYiMBzZj1T1UfkuItQbt2fYtlAiEm96E1veBt+d0X52e6eDzdQeYvjSKs6kZDGxdnVF9GlK5jOcMV6pUQeBKgtgBlAVO5PXGjTHLgeXO5/uA9nm9DeWi1LOwZhr88X/WpaXOT0PnJ8Hfc4bJNMawYMcxJi2I5EBcMp3rBfJCv8Y0qVr6yisrpf7BlQRRFogUkQ38vQZx1be5qgIkKxM2zbYK0Eknofkd0PNlKFvD7shyZfPBM4yfG8HGA2doEFSKWfe3o2uDitrTqlLXwJUEMdbtUaj8ZwxELbLqDKd2Q82O8K9voVobuyPLlUOnk3lj4W5+33qEwFJ+vH57cwa3qY6PDvWp1DVzpSX1ChEJAto5Z603xuT55SaVj45th0Uvwb7lUL4O3PkFNOrvUQXoc6kZvLd8L5+s3o+XwOM96vHvrnUppQVopfKMKy2p7wCmYNUKBJghIs8aY35wc2wqr509CmHjYcuXULws3DAJ2g4HH8/pvjozy8E3Gw7x9uI9xCWlc3uraozq25CqZT1nuFKlPIUrP7deBNqdP2sQkYrAEkAThKdIS4S102HtDHBkWuM/dxkFxT1nmExjDMt2n2DivEiiTyTSIbg8s25qQvPqnlNEV8rTuJIgvC66pBQH6AVeT+DIss4WwiZA4jFoehv0HAvlg+2OLFd2Hklg4rwI1kTHERxYkg/uaUPfpkFagFbKzVxJEAtEZCHwtXP6TnQ0uYJvb5hVgD6+A6q3hzs/hxqedRfxsYRUpi7azY+bYilb3JdxNzdhyHW18NUCtFL5wpUi9bPOwYE6O2fNNMb87N6w1FU7vgsWj4HoJVC2FgyeBU1u9agCdFJaJh+u2MvMVftwOGBEaB0e6V6PMsU9p7GeUoXBJROEiNQDgowxa5yjxv3knN9ZROoaY/bmV5DKBeeOw7IJsPlzq3fVPuOh/Qjw8Zzuq7Mchu/DD/Hm4j2cPJdG/5AqPHdDI2qUL2F3aEoVSZc7g5iGs8+kiyQ4X7vZLRGp3ElPtlo/r5kGmanQ/t/Q9b9QwrOGyVy55yQT50UQeewcbWqV48OhbWhd03OK6EoVRpdLEEHGmO0XzzTGbBeR2m6LSLnG4YBt38DS1+DcEasdQ+9XoUJduyPLlT3HzzFhbgQr9pykZvkSvDekNTc2q6wFaKUKgMsliLKXeU1vOrfTvhWw6EWrwVvVVjDwY6jdye6ocuXkuTTeWryHbzccpJSfDy/2a8y9HWvh56NDfSpVUFwuQYSLyEPGmI+yzxSRB4GN7g1L5ejkHqsAvWcBlKkBt38MzQaCl+fc1ZOSnsUnq/fx/vK9pGU6uK9jbUb2qE+5kp7TWE+pouJyCeJJ4GcRGcJfCaEtUAy4zd2BqWySTlmd6YX/D4qVtNoyXPcf8PWcEzmHw/DLlsNMWbibowmp9G0axOgbGxMcWNLu0JRSl3DJBGGMOQ50FJHuQDPn7LnGmLB8iUxBRiqsew9WvQUZydD2Aeg2GkoG2h1ZrvyxN44J83ax4/BZmlcrw7Q7W9KhTgW7w1JKXYEr7SCWAcvyIRZ1nsMBO36Epa9AwiFocKNVgK7YwO7IcmXvyURenxfJkojjVC3jz7Q7W3JLi6p46VCfSnkE7fqyoDmwFha+CEc2QeUQuPU9CO5id1S5cjopnWlL9vDlnwcp7uvNf29oyAOdgvH31QK0Up5EE0RBEbcXFr8MkXMgoCrc+gGE3OlRBejUjCxmrY3h/8KiSc7I4u72NXiyVwMCS3lOYz2l1F80Qdgt+TSsmAwbPgZvP+j+ktXbajHPaT3scBh+33aENxbs5nB8Cj0bVeL5fo2oVynA7tCUUtdAE4RdMtNg/UxYOQXSzkHre6HbCxAQZHdkubJ+/2kmzN3F1tgEmlYtzZRBIXSs51lFdKVUzjRB5DdjYNcvsHgsxB+Aer2g92sQ1MTuyHJl/6kkJs2PYOHO41Qu7c+bg1twW6tqWoBWqhDRBJGfDm2wWkAf+hMqNYV7foJ6Pe2OKlfik9OZvjSaz9fF4OvtxTO9G/BgaB2KF9MCtFKFjSaI/HAmBpaMg50/Q6kguGUGtBwCXp7zpZqe6eDzdQeYvjSKc6kZ3NmuBk/1bkClAH+7Q1NKuYkmCHdKOQMrp1q1Bi8f6PocdBwJfqXsjsxlxhgW7jzGpPmRxMQlE1o/kBdvakyjyqXtDk0p5WaaINwhMx3CP4UVkyAlHloNge4vQumqdkeWK5sOnmHi3AjCD5yhfqVSzLq/Hd0aVrI7LKVUPtEEkZeMsdoxLB4Lp/dCnW7WwD2Vm9sdWa4cjEtm8sJI5m47SsUAPybd3pxBbarjo0N9KlWkaILIK4c3wsKX4OBaqNgI/vU91O/tUUN9JiRnMCMsitl/xODj5cUTPeszoksdSvrpx0Spokj/51+r+IOw9FXY/j2UrAj934ZW94K35xzajCwHX6w7wDtLo0hIyeCONjV4uk8DgkprAVqposxzvsUKmtQEq5fVde9bZwmhz0CnJ8Hfc4q3xhiWRJzg9XkR7DuVROd6gbzQrzFNqnrOPiil3EcTRG5lZcDGWdb4DMlxVn9JPV+GMtXtjixXdh5JYPycCP7YF0fdiiX537B2dGtYUYf6VEpdoAnCVcZYI7ktfhlO7YFanaHveGvITw9yNCGFqQv38NPmWMoW9+W1AU25q31NfLUArZS6iCYIVxzdanXBHbMKKtSDu76Ghjd6VAH6XGoGH6zYy8er9mMMjAitwyPd61GmuK/doSmlCihNEJeTcBjCxsPWr6F4ObhxCrS9H7w950s1I8vBN+sPMm1JFHFJ6QxoWZVRfRpSo7zn9BarlLJHvicIEakBfAYEAQaYaYx5R0TKA98CtYEY4A5jzJn8jg+weldd8w6sfReMAzqNhM5PQ/GytoRzNYwxhEWeYMK8CPadTKJ9cHk+7deYFjU8Zx+UUvay4wwiE3jGGLNJRAKAjSKyGBgGLDXGTBKR0cBo4Ll8jSwrEzZ/DssmQtIJaDbIKkCXq5WvYVyriKNnmTA3gtXRpwgOLMnMoW3o3SRIC9BKqVzJ9wRhjDkKHHU+PyciEUA1YADQzbnYbGA5+ZkgopbAopfgZATUuA7u/hqqt823zeeFE+dSeWvRHr4LP0SAvy9jb27CkA61KOajBWilVO7ZWoMQkdpAK+BPIMiZPACOYV2CymmdEcAIgJo1a157EMd2wOIxsDcMygXDHZ9B41s8qgCdkp7Fp2v2896yaNIyHdzfKZjHe9SjbIlidoemlPJgtiUIESkF/Ag8aYw5m/3yhzHGiIjJaT1jzExgJkDbtm1zXMYl545ZBegtX4Jfaej7OrR7EHw850vV4TD8vPkwUxft5mhCKr2bBPFCv8YEB5a0OzSlVCFgS4IQEV+s5PClMeYn5+zjIlLFGHNURKoAJ9wWwOYvYN5/ISsdrnvEagVdorzbNucOa6NPMWFeBDuPnKVF9TJMu7MlHepUsDsspVQhYsddTAJ8AkQYY97K9tJvwH3AJOfjr24LomwtqN8Leo2D8nXcthl3iDp+jtfnRxIWeYJqZYvzzl0tuTmkqg71qZTKc3acQXQChgLbRWSLc94LWInhOxEZDhwA7nBbBMGh1p8HOZWYxrQle/h6/SFKFPPm+RsbcV/H2vj7es6odEopz2LHXUyrgUv93PWsAZrzQWpGFv9bE8N7y6JJzsjing41eaJXA8qX9JxaiVLKM2lL6gLKGMOcbUeZND+Sw/Ep9GpcidE3NqZeJc8ZrlQp5dk0QRRA4TGnmTAvgs0H42lcpTRvDAqhU71Au8NSShUxmiAKkP2nkpg8P5IFO49RKcCPNwaGMLBNdby1AK2UsoEmiALgdFI605dG8cW6AxTz8eLp3g14MDSYEsX0n0cpZR/9BrJRWmYWs9fGMCMsmqS0TO5sV5OnetenUoAO9amUsp8mCBsYY5i/4xivz4/g0OkUujWsyAv9GtMgKMDu0JRS6gJNEPls66F4xs/dxYaYMzQIKsVnD7SnS4OKdoellFL/oAkinxyJT2Hqwt38tPkwgaWKMfG25tzRtjo+OtSnUqqA0gThZolpmby/PNoa6hP4T7e6PNKtLgH+njMqnVKqaNIE4SaZWQ6+2XCIaUv2cCoxnVtbVmVU34ZUL6dDfSqlPIMmiDxmjGHZ7hNMnBdJ9IlE2tcuzyf36VCfSinPowkiD0UeO8v4OX8N9fnh0Db00aE+lVIeShNEHjh5Lo23Fu/h2w0HdahPpVShoQniGqRmZPHJ6r+G+hzWMZiRPXWoT6VU4aAJ4ioYY/h921EmO3ta7d0kiOdvbESditrTqlKq8NAEkUsbD5xh/NxdbD4YT5MqpZkyOISOdbWnVaVU4aMJwkWHTiczeUEkc7YdpVKAH1MGhXB7a+1pVSlVeGmCuILEtEzeDYvm0zX78RJ4omd9RnSpQ0k/PXRKqcJNv+UuweEw/LgpljcW7ubkuTRub1WNZ29oSJUyxe0OTSml8oUmiBxsPHCaV37fxbbYBFrVLMtH97alpTZ0U0oVMZogsjmakMKk+ZH8uuUIQaX9mHZnS25pURUvrTMopYogTRBYA/d8sno/M5ZG4zCGx3vU4+GudbXOoJQq0or8N+CKPSd55bed7DuVRJ8mQYzp34Qa5bVDPaWUKrIJ4tDpZMbP3cXCnccJDizJrPvb0a1hJbvDUkqpAqNIJojvNhxizK878BLh2b4NeTA0GD8fb7vDUkqpAqVIJojagSXp2bgSL97UhGpl9bZVpZTKSZFMEO2Dy9M+uLzdYSilVIGm/VErpZTKkSYIpZRSOdIEoZRSKkeaIJRSSuVIE4RSSqkcaYJQSimVI00QSimlcqQJQimlVI7EGGN3DFdNRE4CB65y9UDgVB6GU5AVlX0tKvsJRWdfi8p+Qv7uay1jTMUrLeTRCeJaiEi4Maat3XHkh6Kyr0VlP6Ho7GtR2U8omPuql5iUUkrlSBOEUkqpHBXlBDHT7gDyUVHZ16Kyn1B09rWo7CcUwH0tsjUIpZRSl1eUzyCUUkpdRpFMECJyg4jsFpFoERltdzx5SUQ+FZETIrIj27zyIrJYRKKcj+XsjDEviEgNEVkmIrtEZKeIPOGcX6j2VUT8RWS9iGx17ucrzvnBIvKn8zP8rYgUszvWvCIi3iKyWUTmOKcL3b6KSIyIbBeRLSIS7pxX4D67RS5BiIg38H/AjUAT4G4RaWJvVHlqFnDDRfNGA0uNMfWBpc5pT5cJPGOMaQJcBzzq/HcsbPuaBvQwxrQAWgI3iMh1wGTgbWNMPeAMMNzGGPPaE0BEtunCuq/djTEts93aWuA+u0UuQQDtgWhjzD5jTDrwDTDA5pjyjDFmJXD6otkDgNnO57OBW/M1KDcwxhw1xmxyPj+H9YVSjUK2r8aS6Jz0df4ZoAfwg3O+x+/neSJSHbgJ+Ng5LRTSfc1BgfvsFsUEUQ04lG061jmvMAsyxhx1Pj8GBNkZTF4TkdpAK+BPCuG+Oi+5bAFOAIuBvUC8MSbTuUhh+gxPA/4LOJzTFSic+2qARSKyUURGOOcVuM9ukRyTuigzxhgRKTS3rolIKeBH4EljzFnrB6elsOyrMSYLaCkiZYGfgUY2h+QWItIfOGGM2Sgi3eyOx806G2MOi0glYLGIRGZ/saB8doviGcRhoEa26erOeYXZcRGpAuB8PGFzPHlCRHyxksOXxpifnLML5b4CGGPigWXA9UBZETn/A6+wfIY7AbeISAzWpd8ewDsUwn01xhx2Pp7ASvrtKYCf3aKYIDYA9Z13RhQD7gJ+szkmd/sNuM/5/D7gVxtjyRPOa9OfABHGmLeyvVSo9lVEKjrPHBCR4kBvrHrLMmCQczGP308AY8zzxpjqxpjaWP8vw4wxQyhk+yoiJUUk4PxzoA+wgwL42S2SDeVEpB/WtU5v4FNjzASbQ8ozIvI10A2rZ8jjwFjgF+A7oCZW77d3GGMuLmR7FBHpDKwCtvPX9eoXsOoQhWZfRSQEq2DpjfWD7jtjzKsiUgfrV3Z5YDNwjzEmzb5I85bzEtMoY0z/wravzv352TnpA3xljJkgIhUoYJ/dIpkglFJKXVlRvMSklFLKBZoglFJK5UgThFJKqRxpglBKKZUjTRBKKaVypAlCuUREjIi8mW16lIiMy6P3niUig6685DVvZ7CIRIjIsovmdzvfc2gO68wTkbLOv0dyub2Pc9MRpDOOjtmm8+W45JaIPCkiJbJNzzvfVkMVLpoglKvSgNtFJNDuQLLL1sLWFcOBh4wx3V1dwRjTz9mCuSyQqwRhjHnQGLMrF6t0AzpeaaEC4EngQoLIdoxUIaMJQrkqE2tIxKcufuHiX7oikuh87CYiK0TkVxHZJyKTRGSIc3yD7SJSN9vb9BKRcBHZ4+yT53wndVNEZIOIbBORf2d731Ui8hvwjy9gEbnb+f47RGSyc97LQGfgExGZksP+lRaRuWKNE/KBiHg514txJsVJQF1n//1TRKSKiKx0Tu8QkdAc4lguIm3PHxMRmSDWuA7rRCToomVrAw8DTznf8/z7dRGRtc7jl/0YP5vtuLySw/6c3+YUscaRWCIi7Z0x7RORW85v13ksNzn/OmY7xstF5AcRiRSRL8UyEqgKLDt/Jnb+GDlbCM917uMOEbnT+fokscbt2CYiU53zbhZrjIfNzths75hO5cAYo3/6d8U/IBEoDcQAZYBRwDjna7OAQdmXdT52A+KBKoAfVh86rzhfewKYlm39BVg/WOpj9djpD4wAXnIu4weEA8HO900CgnOIsypwEKiI1Uo1DLjV+dpyoG0O63QDUoE6WC2WF5/fH+f+BgK1gR3Z1nkGeNH53BsIyOF9L2wPq/fOm53P3zi/XxctPw6r9TDZjsv3zuPSBKuberC6ZpgJiPO1OUCXHN7PADc6n/8MLMLqLrwFsMU5vwTg73xeHwjPdkwSsPo+8gL+wOpg7sIxybad88doIPBRtvllsHpj3c1fjXLLOh/LZZv3IPCm3Z9x/fvnn55BKJcZY84CnwEjc7HaBmON3ZCG1U31Iuf87Vhfuud9Z4xxGGOigH1YPZb2Ae4Vq6vrP7G+bOo7l19vjNmfw/baAcuNMSeN1UX0l0AXF+Jcb6wxQrKAr7HONi67X8D9zjpMc2ONSXE56Vhf5AAb+fu+X84vzuOyi7+6f+7j/NsMbMI6VvVzWDcdK/GCdbxXGGMy+Pux9wU+EpHtWMkoe81kvTEm1hjjALa4EPN2oLeITBaRUGNMAlaSScU6c7sdSHYuWx1Y6Nzus0DTK7y3soEmCJVb07Cu5ZfMNi8T52fJeWkm+5CQ2fvMcWSbdvD37uYv7vPFYP1CftxYo261NMYEG2POJ5ika9qLf8pp+5de2BqYqQvWWdEsEbn3Cu+fYZw/l4EsXO9qP/vxk2yPr2c7LvWMMZ9cYZsXjr3zC//89p/C6rOrBdCWS//bXTFmY8weoDVWohgvIi87k3R7rAF/+vNXwpoBvGuMaQ78G+uMURUwmiBUrhir87Dv+PuwjzFAG+fzW7B+lebWYBHxctYl6mBdllgI/Eesbr0RkQZi9X55OeuBrs5r4t7A3cAKF7bfXqwefr2AO4HVF71+Dgg4PyEitYDjxpiPsEY/a+3CNq7kb9u4jIXAA2KNhYGIVBNrXIGrUQY46kwaQ7Eul11VnCJSFUg2xnwBTAFaO2MsY4yZh5WMWmTb7vluu++7+L1UwaADBqmr8SbwWLbpj4BfRWQr1i/Eq/l1fxDry7008LAxJlVEPsa6rLFJRAQ4yRWGYTTGHBWR0VhdRAsw1xjjSrfJG4B3gXrOdX/O/qIxJk5E1ojIDmA+VvfMz4pIBlZ95kpnEK74HfhBRAYAj19qIWPMIpH/b+8ObRAGwjAMf79jDSSCjVgAywR4PIIVCALHBmyCxeBIANEakl+QJrjn0RWXnnhz18u1Fkmuw2vJI8kq0/4fsE9yGldAv87dIcmlqm7v7xNhyyS7qnoleSZZZwjJuapmGeZjMz67TXKsqnuG70TzCWPnz9zmCkDLFhMALYEAoCUQALQEAoCWQADQEggAWgIBQEsgAGh9AL0/qg6m34MNAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, ax = subplots()\n",
    "ax.plot(compression_rate, label=\"measured\")\n",
    "ax.plot(theoretical, label=\"worse case\")\n",
    "ax.set_xlabel(\"Number of bits in the mantissa\")\n",
    "ax.set_ylabel(\"Compression ratio\")\n",
    "ax.legend()\n",
    "fig.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/mntdirect/_scisoft/users/jupyter/jupy35/lib/python3.5/site-packages/matplotlib/figure.py:459: UserWarning: matplotlib is currently using a non-GUI backend, so cannot show the figure\n",
      "  \"matplotlib is currently using a non-GUI backend, \"\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYsAAAEKCAYAAADjDHn2AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzt3XucXHV9//HXey/JbshudpOAkhubQKzGCoFuAt4gKmKsFbRCgUpNLS0g0kq1Vnr5gaD+0CI/8WFpBblWUQSsNbWxSLkVuWaBlEAoEEKADWpi7iHZzV4+vz/O2WQy2c3MbvbsbGbez0enmXOdzxnWfe8533O+X0UEZmZm+1JV6gLMzGz0c1iYmVlBDgszMyvIYWFmZgU5LMzMrCCHhZmZFeSwMDOzghwWZmZWkMPCzMwKqil1AcNl8uTJ0dLSUuoyzMwOKI8//vhvIuLgQuuVTVi0tLTQ1tZW6jLMzA4okl4uZj1fhjIzs4IcFmZmVpDDwszMCiqbNgszs8Hq6uqivb2djo6OUpeSubq6OqZNm0Ztbe2QtndYmFnFam9vp6GhgZaWFiSVupzMRATr16+nvb2dmTNnDmkfmV6GkrRQ0nOSVkq6qJ/ln5W0QtJTku6WdFjOsh5Jy9LX4izrNLPK1NHRwaRJk8o6KAAkMWnSpP06g8rszEJSNXA18H6gHVgqaXFErMhZ7UmgNSK2S/oU8A/A6emyHRExN6v6zMyAsg+KPvt7nFmeWcwHVkbEqojYCdwKnJK7QkTcGxHb08lHgGkZ1tOvrR1dfOOu51n26qaR/mgzswNGlmExFXg1Z7o9nTeQs4Gf5UzXSWqT9Iikj/S3gaRz0nXa1q1bN6Qie3qDb979Ak+8vHFI25uZVYJR0cAt6SygFTghZ/ZhEbFG0izgHknLI+LF3O0i4lrgWoDW1tYYymePH5t8BVs7uoeyuZlZpnp6eqiurh5weiDd3d3U1Azfr/gszyzWANNzpqel8/Yg6UTg74CTI6Kzb35ErEn/XQXcBxydRZE11VWMG1PN1o6uLHZvZrZP3/ve95g/fz5z587l3HPPpaenh/Hjx/O5z32Oo446iocffpiWlha+8IUvcMwxx3D77bezbNkyjjvuOI488kg++tGPsnFjcmVkwYIFXHjhhbS2tvLNb35zWOvM8sxiKTBb0kySkDgD+MPcFSQdDVwDLIyItTnzm4HtEdEpaTLwTpLG70w01NX4zMKswl3678+w4rUtw7rPOVMaueTDbx1w+bPPPssPf/hDHnzwQWprazn//PO55ZZbeP311zn22GO58sord607adIknnjiCQCOPPJIvvWtb3HCCSdw8cUXc+mll3LVVVcBsHPnzkz6ycssLCKiW9IFwJ1ANXBDRDwj6TKgLSIWA1cA44Hb05b6VyLiZOAtwDWSeknOfr6adxfVsGqoq2WLzyzMbITdfffdPP7448ybNw+AHTt2cMghh1BdXc3HPvaxPdY9/fTkRtHNmzezadMmTjghuWq/aNEiTjvttL3WG26ZtllExBJgSd68i3PenzjAdg8Bb8uytlw+szCzfZ0BZCUiWLRoEZdffvke87/+9a/v1S5x0EEHFbXPYtcbLPcNBTTW1brNwsxG3Pve9z7uuOMO1q5NrsJv2LCBl1/ed4/hEyZMoLm5mQceeACA7373u7vOMrI0Ku6GKrWGuhpe3bC98IpmZsNozpw5fPnLX+akk06it7eX2tparr766oLb3XzzzZx33nls376dWbNmceONN2Zeq8OCvjYLX4Yys5F3+umn79XOsG3btj2mV69evcf03LlzeeSRR/ba13333Tfc5e3iy1BAY12NL0OZme2Dw4LkMlRndy+d3T2lLsXMbFRyWJBchgI/xW1WiSKG1PnDAWd/j9NhATTWu8sPs0pUV1fH+vXryz4w+sazqKurG/I+3MANNIztO7Nwu4VZJZk2bRrt7e0MtSPSA0nfSHlD5bAgabMAn1mYVZra2tohjxxXaXwZitw2C59ZmJn1x2HB7jOLLTt8ZmFm1h+HBUl3H4A7EzQzG4DDAhjvNgszs31yWADVVWL8WPc8a2Y2EIdFqsFdfpiZDchhkfKYFmZmA3NYpDxanpnZwBwWKZ9ZmJkNzGGRavBoeWZmA3JYpBp9ZmFmNiCHRSo5s3BYmJn1x2GRaqirYWdPLx1dHgDJzCyfwyLV2Nc/lNstzMz24rBIebQ8M7OBOSxSHtPCzGxgDotUY73HtDAzG4jDIuUzCzOzgTksUh4tz8xsYA6LlEfLMzMbmMMiNX5MDZLPLMzM+pNpWEhaKOk5SSslXdTP8s9KWiHpKUl3SzosZ9kiSS+kr0VZ1glQVSXGj6lhi9sszMz2kllYSKoGrgY+CMwBzpQ0J2+1J4HWiDgSuAP4h3TbicAlwLHAfOASSc1Z1dqnsd5dfpiZ9SfLM4v5wMqIWBURO4FbgVNyV4iIeyNiezr5CDAtff8B4K6I2BARG4G7gIUZ1gp4tDwzs4FkGRZTgVdzptvTeQM5G/jZELcdFh7TwsysfzWlLgBA0llAK3DCILc7BzgHYMaMGftdR0NdLb/e0rHf+zEzKzdZnlmsAabnTE9L5+1B0onA3wEnR0TnYLaNiGsjojUiWg8++OD9LthnFmZm/csyLJYCsyXNlDQGOANYnLuCpKOBa0iCYm3OojuBkyQ1pw3bJ6XzMuU2CzOz/mV2GSoiuiVdQPJLvhq4ISKekXQZ0BYRi4ErgPHA7ZIAXomIkyNig6QvkQQOwGURsSGrWvs0pgMgRQRpPWZmRsZtFhGxBFiSN+/inPcn7mPbG4Absqtubw11tXT3Bh1dvdSPqR7JjzYzG9X8BHeO3Z0J+lKUmVkuh0WOBo+WZ2bWL4dFjsa051l3+WFmtieHRQ6PaWFm1j+HRQ6PaWFm1j+HRY7Gep9ZmJn1x2GRw2cWZmb9c1jkOGhMNVXyaHlmZvkcFjkkMX6su/wwM8vnsMjTUOcBkMzM8jks8jTUeWhVM7N8Dos8ydCqvgxlZpbLYZGn0WNamJntxWGRp6Gu1n1DmZnlcVjk8Wh5ZmZ7c1jkaairYVtnMgCSmZklHBZ5Gupq6ekNtu/sKXUpZmajhsMiT+OuLj98KcrMrI/DIo9HyzMz25vDIs/u0fJ8ZmFm1sdhkadh12h5PrMwM+vjsMjT6NHyzMz24rDI4zEtzMz25rDI49HyzMz25rDIU19bTXWVfGZhZpbDYZFHkrv8MDPL47DoR0NdDVt2+MzCzKyPw6IfDWM9Wp6ZWS6HRT98GcrMbE8Oi354TAszsz3tMywkVUl6x0gVM1o01vvMwsws1z7DIiJ6gauHunNJCyU9J2mlpIv6WX68pCckdUs6NW9Zj6Rl6WvxUGsYisY6j8NtZparmMtQd0v6mCQNZseSqkmC5oPAHOBMSXPyVnsF+GPg+/3sYkdEzE1fJw/ms/dXQ10NWzu76e31AEhmZlBcWJwL3A7slLRF0lZJW4rYbj6wMiJWRcRO4FbglNwVImJ1RDwF9A628Cw11NUQAa/v9KUoMzMoIiwioiEiqiKiNiIa0+nGIvY9FXg1Z7o9nVesOkltkh6R9JH+VpB0TrpO27p16wax631r8ABIZmZ7qClmJUknA8enk/dFxE+zK2mXwyJijaRZwD2SlkfEi7krRMS1wLUAra2tw3bNqME9z5qZ7aHgmYWkrwKfAVakr89IuryIfa8BpudMT0vnFSUi1qT/rgLuA44udtv91eieZ83M9lBMm8XvAu+PiBsi4gZgIfChIrZbCsyWNFPSGOAMoKi7miQ1Sxqbvp8MvJMkqEaEzyzMzPZU7EN5TTnvJxSzQUR0AxcAdwLPArdFxDOSLksvayFpnqR24DTgGknPpJu/BWiT9D/AvcBXI2IEw8Kj5ZmZ5SqmzeJy4ElJ9wIiabvY65mJ/kTEEmBJ3ryLc94vJbk8lb/dQ8DbivmMLDR6HG4zsz3sMyzSZyt+ARwHzEtnfyEifpV1YaXk0fLMzPa0z7CIiJC0JCLeRpHtDeWgrraKmiq5zcLMLFVMm8UTkuYVXq18SKKx3l1+mJn1KabN4ljg45JeBl4nabeIiDgy08pKzN2Um5ntVkxYfCDzKkYhj5ZnZrZboQbuauDOiHjzCNUzani0PDOz3Qp1Ud4DPCdpxgjVM2r4MpSZ2W7FXIZqBp6R9BhJmwUAI91t+Ehr8JgWZma7FBMW/yfzKkYhj5ZnZrZbwbCIiPslHQbMjoj/kjQOqM6+tNJqqKtl285kAKSqqkGN+2RmVnaK6XX2z4A7gGvSWVOBf8uyqNGgMR0AaWunzy7MzIp5KO/TJL2+bgGIiBeAQ7IsajTY3fOs2y3MzIoJi850WFQAJNUAZT84tUfLMzPbrZiwuF/S3wL1kt5PMh73v2dbVul5TAszs92KCYuLgHXAcuBcki7H/z7LokaDiQeNAeA32zpLXImZWekVczdUL/Cd9FUxpjbVA/Daph0lrsTMrPSKHSmv4kyor+WgMdWscViYmTksBiKJKU31PrMwM6O45yzq+pk3OZtyRpcpTfU+szAzo7gzi6WSjuubkPQx4KHsSho9pjbX89qmjlKXYWZWcsX0DfWHwA2S7gOmAJOA92ZZ1GgxtameDa/vZMfOHurHlH0PJ2ZmAyrmbqjlkr4CfBfYChwfEe2ZVzYKTGlKrsCt2bSDIw4ZX+JqzMxKp5g2i+uBC4EjgU8CP5X06awLGw2mNo0DfPusmVkxbRbLgfdExEsRcSfJmNzHZFvW6NB3ZuGwMLNKV8xlqKvypjcDZ2dW0SjyhsY6qoTviDKzilcwLCTNBi4H5gC7bqONiFkZ1jUq1FZX8cbGOoeFmVW8Yi5D3Qj8M9ANvAf4F+B7WRY1mvjBPDOz4sKiPiLuBhQRL0fEF4EPZVvW6JGEhZ+1MLPKVtR4FpKqgBckXSDpo0DF3Ec6tbmeX27eQW9v2Q/hYWY2oGLC4jPAOOAvgN8B/ghYlGVRo8mUpnq6eoJ17qrczCpYwbCIiKURsS0i2iPikxHx+xHxSDE7l7RQ0nOSVkq6qJ/lx0t6QlK3pFPzli2S9EL6Klk4Tc15MM/MrFIV81Beq6Qfp7/Un+p7FbFdNXA18EGSO6nOlDQnb7VXgD8Gvp+37UTgEpJnOuYDl0hqLuaAhlvfg3lrNjoszKxyFdM31C3A50kezusdxL7nAysjYhWApFuBU4AVfStExOp0Wf5+PwDcFREb0uV3AQuBHwzi84eFH8wzMysuLNZFxOIh7Hsq8GrOdDvJmcJQt506hBr2W0NdLQ11NQ4LM6toxYTFJZKuA+4GdrXyRsS/ZlZVkSSdA5wDMGPGjMw+Z6rHtTCzCldMWHwSeDNQy+7LUAEUCos1wPSc6WnpvGKsARbkbXtf/koRcS1wLUBra2tm97YmYeFnLcyschUTFvMi4reGsO+lwGxJM0l++Z9BMjZGMe4E/m9Oo/ZJwN8MoYZhMaWpnraXN5bq483MSq6Y5ywe6ucupoIiohu4gOQX/7PAbRHxjKTLJJ0MIGmepHbgNOAaSc+k224AvkQSOEuBy/oau0thSlM9m3d0sa2zu1QlmJmVVDFnFscByyS9RNJmISAi4shCG0bEEmBJ3ryLc94vJbnE1N+2NwA3FFFf5qY21wPJHVFvekNDiasxMxt5xYTFwsyrGOVyH8xzWJhZJSpmPIuXR6KQ0WxKU3Jm4QfzzKxSFdNmUfEOaaijpkp+1sLMKpbDogjVVeKNE+ocFmZWsRwWRfK4FmZWyRwWRZrmp7jNrII5LIo0pameX23poLtnMH0pmpmVB4dFkaY01dPTG6zd6kGQzKzyOCyK1Pdgni9FmVklclgUaarHtTCzCuawKNKuB/McFmZWgRwWRRo3pobmcbV+itvMKpLDYhCSZy0cFmZWeRwWg+AH88ysUjksBqFveNWIzAblMzMblRwWgzC1qZ5tnd1s6fAgSGZWWRwWg9B3R5TbLcys0jgsBmHXg3m+I8rMKozDYhCm9D2Yt9lhYWaVxWExCJMPGsuY6io/mGdmFcdhMQhVVeLQpjrfPmtmFcdhMUhTm+pZs3F7qcswMxtRDotB8oN5ZlaJHBaDNKWpnl9v7aDLgyCZWQVxWAzS9OZ6IuCVDb4UZWaVw2ExSHOnNwHw+OqNJa7EzGzkOCwG6YhDxtM8rpalqzeUuhQzsxHjsBgkSbS2THRYmFlFcVgMwbyWZlav387arb4ryswqg8NiCOa1TASgze0WZlYhHBZD8NtTJ1BXW8VjL/lSlJlVhkzDQtJCSc9JWinpon6Wj5X0w3T5o5Ja0vktknZIWpa+vp1lnYNVW13F0dObaXvZYWFmlSGzsJBUDVwNfBCYA5wpaU7eamcDGyPiCOAbwNdylr0YEXPT13lZ1TlU82ZOZMVrW9ja0VXqUszMMpflmcV8YGVErIqIncCtwCl565wC3Jy+vwN4nyRlWNOwmdfSTG/AE69sKnUpZmaZyzIspgKv5ky3p/P6XSciuoHNwKR02UxJT0q6X9K7+/sASedIapPUtm7duuGtvoBjZjRTXSXafAutmVWA0drA/UtgRkQcDXwW+L6kxvyVIuLaiGiNiNaDDz54RAs8aGwNb53S6EZuM6sIWYbFGmB6zvS0dF6/60iqASYA6yOiMyLWA0TE48CLwJsyrHVIWg+byLJXN7Gz250Kmll5yzIslgKzJc2UNAY4A1ict85iYFH6/lTgnogISQenDeRImgXMBlZlWOuQzJ/ZTGd3L8vXbC51KWZmmcosLNI2iAuAO4Fngdsi4hlJl0k6OV3temCSpJUkl5v6bq89HnhK0jKShu/zImLUXe9pTR/Oc9cfZlbuarLceUQsAZbkzbs4530HcFo/2/0I+FGWtQ2HyePHMmvyQUkj9wmHl7ocM7PMjNYG7gPGvJaJLF29kd7eKHUpZmaZcVjsp9aWZjbv6OKFtdtKXYqZWWYcFvtp/ky3W5hZ+XNY7KcZE8dxSMNYh4WZlTWHxX6SlLRb+OE8MytjDothMK+lmdc2d7Bm045Sl2JmlgmHxTCY19du4bMLMytTDoth8OY3NtIwtsbtFmZWthwWw6C6ShxzWLPDwszKlsNimLx79mSe//U2lre7nygzKz8Oi2Fy+rzpNNTV8M/3ryx1KWZmw85hMUwa6mr5xNsP42dP/4oX1/lpbjMrLw6LYfTJd85kbE0V377vxVKXYmY2rBwWw2jy+LGcMW8GP35yDa/5mQszKyMOi2H2Z8fPAuA7D4y6sZrMzIbMYTHMpjbVc8rcqfzgsVdYv62z1OWYmQ0Lh0UGPrVgFp3dvdz00OpSl2JmNiwcFhk44pAGPjDnjdz80Gq2dnSVuhwzs/3msMjI+e85nC0d3Xz/0VdKXYqZ2X5zWGTkyGlNvOuIyVz3i5fo6OopdTlmZvvFYZGh8xcczrqtndzxeHupSzEz2y8Oiwy9/fBJzJ3exDfuep5VfqrbzA5gDosMSeLKPzgKgLOue5T2jdtLXJGZ2dA4LDJ2+MHj+Zez57Ots5uPX/coa7d0lLokM7NBc1iMgLdOmcBNfzKfdVs7Oev6R9n4+s5Sl2RmNigOixFyzIxmrvtEK6vXb2fRjY/5+QszO6A4LEbQO46YzD9//BhWvLaFs29qY8dO31JrZgcGh8UIe99b3sA3Tp9L28sb+Og/Pci9z60lIkpdlpnZPjksSuDDR03h22f9Dq/v7OaTNy7l9Gseoc3jd5vZKOawKJGT3vpG7v7sAr50ylt5af3rnPrth/mTm5ay4rUtpS7NzGwvmYaFpIWSnpO0UtJF/SwfK+mH6fJHJbXkLPubdP5zkj6QZZ2lMqamij96ewv3f34Bf73wt2hbvYEPfesBTvv2Q3zjrudZunoDXT29pS7TzAxldb1cUjXwPPB+oB1YCpwZESty1jkfODIizpN0BvDRiDhd0hzgB8B8YArwX8CbImLAFuHW1tZoa2vL5FhGyubtXdz40Evc879rWb5mMxEwbkw1x86cyDsOn8xvvbGBGRPHMaWpnjE1Pik0s/0n6fGIaC20Xk2GNcwHVkbEqrSgW4FTgBU565wCfDF9fwfwj5KUzr81IjqBlyStTPf3cIb1ltyEcbVceOKbuPDEN7Fp+04eWbWeB1eu58GVv+He557dtV6V4NAJ9UyfWM+05nFMPGgME+praayvpbGuZtf7uppqxtZWMbamirraasbWVDGmporaqiqqqlTCIzWzA02WYTEVeDVnuh04dqB1IqJb0mZgUjr/kbxtp2ZX6ujTNG4MC3/7UBb+9qEArN3SwUu/eZ1XNmzn1Y07eHXDdl7ZsJ1fvPAbNu3YSUfX4C5XSVBTJWqqqqipElVVokpQpT3fi6Tbkr5tJBAinbXn8l3/r5/PG7CO7ELLcWiV4s2HNvKtM4/O9DOyDIvMSToHOAdgxowZJa4mW4c01nFIYx3HzprU7/KOrh62dHSxZUcXm3d0saWjm86uHjq7e+ns6qWju4fOrl46u3vo6YXu3l66e4Oe3qC7J+jp7aU3oDeC3oCIZFkAERAE6f/tutW3b9nu9/1f0hzwQmeGdwxHljs3G2WmN9dn/hlZhsUaYHrO9LR0Xn/rtEuqASYA64vcloi4FrgWkjaLYav8AFRXW01dbTWHNNSVuhQzK0NZtpIuBWZLmilpDHAGsDhvncXAovT9qcA9kfx5uhg4I71baiYwG3gsw1rNzGwfMjuzSNsgLgDuBKqBGyLiGUmXAW0RsRi4Hvhu2oC9gSRQSNe7jaQxvBv49L7uhDIzs2xlduvsSCuHW2fNzEZasbfO+mZ9MzMryGFhZmYFOSzMzKwgh4WZmRXksDAzs4LK5m4oSeuAl/djF5OB3wxTOaNZpRwnVM6xVspxQuUc60ge52ERcXChlcomLPaXpLZibh870FXKcULlHGulHCdUzrGOxuP0ZSgzMyvIYWFmZgU5LHa7ttQFjJBKOU6onGOtlOOEyjnWUXecbrMwM7OCfGZhZmYFVXxYSFoo6TlJKyVdVOp6hpOkGyStlfR0zryJku6S9EL6b3MpaxwOkqZLulfSCknPSPpMOr8cj7VO0mOS/ic91kvT+TMlPZr+HP8wHRbggCepWtKTkn6aTpfrca6WtFzSMklt6bxR9fNb0WEhqRq4GvggMAc4U9Kc0lY1rG4CFubNuwi4OyJmA3en0we6buBzETEHOA74dPrfsRyPtRN4b0QcBcwFFko6Dvga8I2IOALYCJxdwhqH02eAZ3Omy/U4Ad4TEXNzbpkdVT+/FR0WwHxgZUSsioidwK3AKSWuadhExH+TjBOS6xTg5vT9zcBHRrSoDETELyPiifT9VpJfLlMpz2ONiNiWTtamrwDeC9yRzi+LY5U0DfgQcF06LcrwOPdhVP38VnpYTAVezZluT+eVszdExC/T978C3lDKYoabpBbgaOBRyvRY00szy4C1wF3Ai8CmiOhOVymXn+OrgL8GetPpSZTncUIS+D+X9Likc9J5o+rnN8sxuG2Ui4iQVDa3w0kaD/wIuDAitiR/iCbK6VjTUSPnSmoCfgy8ucQlDTtJvwesjYjHJS0odT0j4F0RsUbSIcBdkv43d+Fo+Pmt9DOLNcD0nOlp6bxy9mtJhwKk/64tcT3DQlItSVDcEhH/ms4uy2PtExGbgHuBtwNNkvr++CuHn+N3AidLWk1yefi9wDcpv+MEICLWpP+uJfkDYD6j7Oe30sNiKTA7vcNiDMkY4ItLXFPWFgOL0veLgJ+UsJZhkV7Lvh54NiL+X86icjzWg9MzCiTVA+8naaO5Fzg1Xe2AP9aI+JuImBYRLST/u7wnIj5OmR0ngKSDJDX0vQdOAp5mlP38VvxDeZJ+l+TaaDVwQ0R8pcQlDRtJPwAWkPRg+WvgEuDfgNuAGSS99P5BROQ3gh9QJL0LeABYzu7r239L0m5Rbsd6JEljZzXJH3u3RcRlkmaR/AU+EXgSOCsiOktX6fBJL0P9VUT8XjkeZ3pMP04na4DvR8RXJE1iFP38VnxYmJlZYZV+GcrMzIrgsDAzs4IcFmZmVpDDwszMCnJYmJlZQQ4LGxJJIenKnOm/kvTFYdr3TZJOLbzmfn/OaZKelXRv3vwFfb2c9rPNEklN6ev8QX7edYPpqDKt4x050yPyvQyWpAsljcuZXtL3LIiVD4eFDVUn8PuSJpe6kFw5T/cW42zgzyLiPcVuEBG/mz453QQMKiwi4k8jYsUgNlkAvKPQSqPAhcCusMj5jqyMOCxsqLpJhn78y/wF+X8BS9qW/rtA0v2SfiJplaSvSvp4Oj7DckmH5+zmREltkp5P+wnq60DvCklLJT0l6dyc/T4gaTGw1y9jSWem+39a0tfSeRcD7wKul3RFP8fXKOk/lIx18m1JVel2q9OA/CpweDr+wBWSDpX03+n005Le3U8d90lq7ftOJH1FybgUj0h6Q966LcB5wF+m++zb3/GSHkq/v9zv+PM538ul/RxP32deoWQcjP+SND+taZWkk/s+N/0un0hf78j5ju+TdIek/5V0ixJ/AUwB7u07Q+v7jtInk/8jPcanJZ2eLv+qkrFHnpL09XTeh5WMU/FkWltZdPpYViLCL78G/QK2AY3AamAC8FfAF9NlNwGn5q6b/rsA2AQcCowl6dfn0nTZZ4Crcrb/T5I/ZmaT9C5aB5wD/H26zligDZiZ7vd1YGY/dU4BXgEOJnk69h7gI+my+4DWfrZZAHQAs0ielL6r73jS450MtABP52zzOeDv0vfVQEM/+931eSS9jH44ff8PfceVt/4XSZ5cJud7uT39XuaQdK8PSfcQ1wJKl/0UOL6f/QXwwfT9j4Gfk3RxfhSwLJ0/DqhL388G2nK+k80k/TFVAQ+TdH636zvJ+Zy+7+hjwHdy5k8g6Tn2OXY/ENyU/tucM+9PgStL/TPu154vn1nYkEXEFuBfgL8YxGZLIxl/opOka+2fp/OXk/wC7nNbRPRGxAvAKpKeVU8CPqGke+5HSX7xzE7XfywiXurn8+YB90XEuki6tr4FOL6IOh+LZJyTHuAHJGch+zwu4JNpu83bIhlXY192kvxSB3icPY99X/4t/V5WsLvL6pPS15PAEyTf1ex+tt0+Rg3CAAACdElEQVRJEsKQfN/3R0QXe373tcB3JC0nCabcNpbHIqI9InqBZUXUvBx4v6SvSXp3RGwmCZwOkjO63we2p+tOA+5MP/fzwFsL7NtGmMPC9tdVJNf+D8qZ1036s5Vevskd+jK3H5/enOle9uwyP78fmiD5y/nPIxlNbG5EzIyIvrB5fb+OYm/9ff7AKycDTR1PcrZ0k6RPFNh/V6R/RgM9FD9cQO73p5x/L8/5Xo6IiOsLfOau7z795d/3+X9J0o/YUUArA/+3K1hzRDwPHEMSGl+WdHEa2PNJBjD6PXaH17eAf4yItwHnkpxJ2ijisLD9EknHZrex5/CWq4HfSd+fTPLX6mCdJqkqbceYRXLp4k7gU0q6I0fSm5T00rkvjwEnpNfQq4EzgfuL+Pz5SnojrgJOB36Rt3wr0NA3Iekw4NcR8R2Skd2OKeIzCtnjM/bhTuBPlIzngaSpSsZFGIoJwC/TAPkjkktqQ6pT0hRge0R8D7gCOCatcUJELCEJpqNyPrevu/FF+fuy0vPgRzYcrgQuyJn+DvATSf9D8pfjUP7qf4XkF30jcF5EdEi6juTSxxOSBKyjwFCTEfFLSReRdG0t4D8iopiunpcC/wgckW7749yFEbFe0oOSngZ+RtKl9OcldZG05xQ6syjGvwN3SDoF+POBVoqIn0t6C/Bw8rWwDTiLoY1/8E/Aj9Izo2L/210L/Kek12LPO8veBlwhqRfoAj5FEio/kVRH8t/js+m6XwRul7SRpF1p5hBqtwy511kzMyvIl6HMzKwgh4WZmRXksDAzs4IcFmZmVpDDwszMCnJYmJlZQQ4LMzMryGFhZmYF/X/vp8aFxXEtRQAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "error = [abs(cround_array(random_64, i) - random_64 ).max() for i in range(53)]\n",
    "\n",
    "fig, ax = subplots()\n",
    "ax.plot(error, label=\"error\")\n",
    "ax.set_xlabel(\"Number of bits in the mantissa\")\n",
    "ax.set_ylabel(\"max error\")\n",
    "ax.legend()\n",
    "fig.show()\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Conclusions:\n",
    "\n",
    "Bitshuffle is needed to make loseless algorithms like LZ4 efficient on floating point data.\n",
    "Zeroing out the last bits of the mantissa helps the bitshuffle/LZ4 algorithm to better compress data. \n",
    "A demonstrator is proposed and works on float64. \n",
    "The compression ratio varies linearly with the size of mantissa used. The precision of the data drops accordingly"
   ]
  }
 ],
 "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.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}