Transaction Duration Anomaly Detection.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# INITIAL WHATEVER"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "import math\n",
    "import numpy\n",
    "import numpy.random\n",
    "import matplotlib\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "from scipy.stats import norm\n",
    "from pt.plotter import *\n",
    "\n",
    "set_default_plotter_options({\n",
    "    'num_anomalies': [1, 2, 3, 4, 5],\n",
    "    'transaction_events': [\n",
    "        numpy.floor(10**x) for x in \n",
    "        [2, 2.25, 2.5, 2.75, 3, 3.25, 3.5, 3.75, 4, 4.25, 4.5]\n",
    "    ],\n",
    "})"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Generate Data with Anomalies"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAX4AAAD8CAYAAABw1c+bAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvNQv5yAAAEClJREFUeJzt3X2MZXV9x/H3pywPVm1ZYDArK1001IpGl2bcamgaiw9FbBSjpmBrNy3N2qqNpFZZ1KTQtIlalf7TqGtBtq0VKGqwQrWIUGtjsbO4IOtKeZBaZMMORapoQl349o97sMN0Hu7OPXeZu7/3K7m555z7O+d8v5nZzz17zrlzU1VIktrxE493AZKkA8vgl6TGGPyS1BiDX5IaY/BLUmMMfklqjMEvSY0x+CWpMQa/JDVmzYHc2THHHFMbNmw4kLuUpIm3Y8eO+6pqqq/tHdDg37BhAzMzMwdyl5I08ZL8R5/b81SPJDXG4Jekxhj8ktQYg1+SGmPwS1JjDH5JaozBL0mNMfglqTEGvyQ1xuCXpMYMHfxJDknytSSf7eZPSHJDktuSXJbksPGVKUnqy/4c8b8V2D1n/r3AhVV1IvBd4Ow+C5MkjcdQwZ9kPfAK4C+7+QCnAld0Q7YDZ4yjQElSv4Y94v9z4B3AI9380cADVbWvm78bOK7n2iRJY7Bs8Cf5VWBvVe2Yu3iBobXI+luSzCSZmZ2dXWGZkqS+DHPEfwrwyiR3AZcyOMXz58CRSR79e/7rgXsWWrmqtlXVdFVNT0319j0CkqQVWjb4q+q8qlpfVRuAM4EvVtWvA9cBr+2GbQauHFuVkqTejHIf/7nAHyS5ncE5/4v6KUmSNE779dWLVXU9cH03fSewqf+SJEnj5Cd3JakxBr8kNcbgl6TGGPyS1BiDX5IaY/BLUmMMfklqjMEvSY0x+CWpMQa/JDXG4Jekxhj8ktQYg1+SGmPwS1JjDH5JaozBL0mNGebL1o9I8tUkNyXZleSCbvklSb6VZGf32Dj+ciVJoxrmG7geAk6tqgeTHAp8Ock/dK+9vaquGF95kqS+LRv8VVXAg93sod2jxlmUJGl8hjrHn+SQJDuBvcA1VXVD99KfJrk5yYVJDh9blZKk3gwV/FX1cFVtBNYDm5I8BzgP+Dng+cBRwLkLrZtkS5KZJDOzs7M9lS1JWqn9uqunqh4ArgdOq6o9NfAQ8DFg0yLrbKuq6aqanpqaGrlgSdJohrmrZyrJkd30E4CXAN9Msq5bFuAM4JZxFipJ6scwd/WsA7YnOYTBG8XlVfXZJF9MMgUE2An87hjrlCT1ZJi7em4GTl5g+aljqUiSNFZ+cleSGmPwS1JjDH5JaozBL0mNMfglqTEGvyQ1xuCXpMYY/JLUGINfkhpj8EtSYwx+SWqMwS9JjTH4JakxBr8kNcbgl6TGTFTwb9h61eNdgiRNvIkKfknS6Ib5zt0jknw1yU1JdiW5oFt+QpIbktyW5LIkh42/XEnSqIY54n8IOLWqngdsBE5L8gLgvcCFVXUi8F3g7PGVKUnqy7LBXwMPdrOHdo8CTgWu6JZvB84YS4WSpF4NdY4/ySFJdgJ7gWuAO4AHqmpfN+Ru4LhF1t2SZCbJzOzsbB81S5JGMFTwV9XDVbURWA9sAp610LBF1t1WVdNVNT01NbXySiVJvdivu3qq6gHgeuAFwJFJ1nQvrQfu6bc0SdI4DHNXz1SSI7vpJwAvAXYD1wGv7YZtBq4cV5GSpP6sWX4I64DtSQ5h8EZxeVV9Nsk3gEuT/AnwNeCiMdYpSerJssFfVTcDJy+w/E4G5/slSRPET+5KUmMMfklqjMEvSY0x+CWpMQa/JDXG4Jekxhj8ktQYg1+SGmPwS1JjDH5JaozBL0mNMfglqTEGvyQ1xuCXpMYY/JLUGINfkhozzFcvPi3JdUl2J9mV5K3d8vOTfCfJzu5x+vjLlSSNapivXtwHvK2qbkzyZGBHkmu61y6sqvePrzxJUt+G+erFPcCebvr7SXYDx427MEnSeOzXOf4kGxh8/+4N3aK3JLk5ycVJ1vZcmyRpDIYO/iRPAj4JnFNV3wM+BDwD2MjgfwQfWGS9LUlmkszMzs72ULIkaRRDBX+SQxmE/ser6lMAVXVvVT1cVY8AHwU2LbRuVW2rqumqmp6amuqrbknSCg1zV0+Ai4DdVfXBOcvXzRn2auCW/suTJPVtmLt6TgHeAHw9yc5u2TuBs5JsBAq4C3jjWCqUJPVqmLt6vgxkgZeu7r8cSdK4+cldSWqMwS9JjTH4JakxBr8kNcbgl6TGGPyS1BiDX5IaY/BLUmMMfklqjMEvSY0x+CWpMQa/JDXG4Jekxhj8ktQYg1+SGmPwS1JjhvnqxacluS7J7iS7kry1W35UkmuS3NY9rx1/uZKkUQ1zxL8PeFtVPQt4AfDmJCcBW4Frq+pE4NpuXpK0yi0b/FW1p6pu7Ka/D+wGjgNeBWzvhm0HzhhXkZKk/uzXOf4kG4CTgRuAp1TVHhi8OQDH9l2cJKl/Qwd/kicBnwTOqarv7cd6W5LMJJmZnZ1dSY2SpB4NFfxJDmUQ+h+vqk91i+9Nsq57fR2wd6F1q2pbVU1X1fTU1FQfNUuSRjDMXT0BLgJ2V9UH57z0GWBzN70ZuLL/8iRJfVszxJhTgDcAX0+ys1v2TuA9wOVJzga+DbxuPCVKkvq0bPBX1ZeBLPLyi/stR5I0bn5yV5IaY/BLUmMMfklqjMEvSY0x+CWpMQa/JDXG4Jekxhj8ktQYg1+SGmPwS1JjDH5JaozBL0mNMfglqTEGvyQ1xuCXpMZMTPBv2HrV412CJB0UhvnqxYuT7E1yy5xl5yf5TpKd3eP08ZYpSerLMEf8lwCnLbD8wqra2D2u7rcsSdK4LBv8VfUl4P4DUIsk6QAY5Rz/W5Lc3J0KWttbRZKksVpp8H8IeAawEdgDfGCxgUm2JJlJMjM7O7vC3UmS+rKi4K+qe6vq4ap6BPgosGmJsduqarqqpqemplZapySpJysK/iTr5sy+GrhlsbGSpNVlzXIDknwCeBFwTJK7gT8CXpRkI1DAXcAbx1ijJKlHywZ/VZ21wOKLxlCLJOkAmJhP7kqS+mHwS1JjDH5JaozBL0mNMfglqTEGvyQ1xuCXpMYY/JLUGINfkhpj8EtSYwx+SWqMwS9JjTH4JakxBr8kNcbgl6TGGPyS1BiDX5Ias2zwJ7k4yd4kt8xZdlSSa5Lc1j2vHW+ZkqS+DHPEfwlw2rxlW4Frq+pE4NpuXpI0AZYN/qr6EnD/vMWvArZ309uBM3quS5I0Jis9x/+UqtoD0D0fu9jAJFuSzCSZmZ2dXeHuJEl9GfvF3araVlXTVTU9NTU17t1Jkpax0uC/N8k6gO55b38lSZLGaaXB/xlgcze9Gbiyn3IkSeM2zO2cnwC+Ajwzyd1JzgbeA7w0yW3AS7t5SdIEWLPcgKo6a5GXXtxzLZKkA8BP7kpSYwx+SWqMwS9JjTH4JakxBr8kNcbgl6TGGPyS1BiDX5IaY/BLUmMmLvg3bL3q8S5BkibaxAW/JGk0Br8kNcbgl6TGGPyS1BiDX5IaY/BLUmOW/SKWpSS5C/g+8DCwr6qm+yhKkjQ+IwV/55er6r4etiNJOgA81SNJjRk1+Av4xyQ7kmzpoyBJ0niNeqrnlKq6J8mxwDVJvllVX5o7oHtD2AJw/PHHj7g7SdKoRjrir6p7uue9wKeBTQuM2VZV01U1PTU1NcruJEk9WHHwJ3likic/Og28DLilr8IkSeMxyqmepwCfTvLodv62qj7XS1WSpLFZcfBX1Z3A83qsRZJ0AHg7pyQ1xuCXpMYY/JLUGINfkhpj8EtSYwx+SWqMwS9JjTH4JakxBr8kNcbgl6TGGPyS1BiDX5IaY/BLUmMMfklqjMEvSY0x+CWN3YatVz3u2x1XDat93wsZKfiTnJbk1iS3J9naV1GSpPEZ5Tt3DwH+Ang5cBJwVpKT+ipMkjQeoxzxbwJur6o7q+p/gEuBV/VTliRpXEYJ/uOA/5wzf3e3TJK0iqWqVrZi8jrgV6rqd7r5NwCbqur3543bAmzpZp8J3LqC3R0D3LeiQleng6mfg6kXsJ/VrtV+fqaqpvra6ZoR1r0beNqc+fXAPfMHVdU2YNsI+yHJTFVNj7KN1eRg6udg6gXsZ7Wzn36Mcqrn34ATk5yQ5DDgTOAz/ZQlSRqXFR/xV9W+JG8BPg8cAlxcVbt6q0ySNBajnOqhqq4Gru6plqWMdKpoFTqY+jmYegH7We3spwcrvrgrSZpM/skGSWpNVfX2AE5jcLvm7cDWBV4/HLise/0GYMOc187rlt/K4DbRJbcJnNBt47Zum4etdB+T2A/wUmAH8PXu+dRJ7mfOescDDwJ/OMm9AM8FvgLs6n5GR0xqP8ChwPauj93AeRPyu/ZLwI3APuC18/a/uRt/G7B5kvsBNs75XbsZ+LVl+1luwLAPBhd47wCeDhwG3AScNG/Mm4APd9NnApd10yd14w/vmr6j296i2wQuB87spj8M/N5K9jHB/ZwMPLWbfg7wnUn++cyp4ZPA37FE8K/2XhhcO7sZeF43fzST/bv2euDSbvongbuY94a9SvvZwOAN+K94bFAeBdzZPa/tptdOcD8/C5zYTT8V2AMcuWQeLBfowz6AFwKfnzN/HvOODBjcAfTCOf847gMyf+yj4xbbZrfOfcCa+fve331Maj/z6gjwX8Dhk9wPcAbwZ8D5LB38q7oX4HTgbw6ifztnAX/fLTsa+HfgqNXez5yxl/DYoDwL+Mic+Y8AZ01qPwvUexPdG8Fijz7P8Q/zJxx+PKaq9gH/zeAXabF1F1t+NPBAt435+9rffUxqP3O9BvhaVT00qf0keSJwLnDBEj1MRC8MjsAqyeeT3JjkHRPezxXADxgcSX4beH9V3T8B/YxS3yT182NJNjH4H8QdS40b6XbO+ftcYFkNOWax5Qu9MS01fiX7WMxq72fwYvJs4L3AyxYYN9dq7+cC4MKqejBZaMhQdQ4z5kD0sgb4ReD5wA+Ba5PsqKprFxi/1HaGGXMg+tkEPMzgNMJa4J+TfKGq7lxg/FLbGWZMn/0sZlKzYElJ1gF/zeCaxSNLje3ziH+YP+Hw4zFJ1gA/Ddy/xLqLLb8POLLbxvx97e8+JrUfkqwHPg38ZlUt+Q4/Af38AvC+JHcB5wDv7D4gOIm93A38U1XdV1U/ZPBZl59fpJdJ6Of1wOeq6kdVtRf4F2CpPzOwWvoZpb5J6ockPwVcBby7qv51ufF9nuNfw+AiyQn838WKZ88b82YeewHk8m762Tz2AsidDC5+LLpNBhcA514AedNK9jHB/RzZrf+ag+HnM6+O81n6HP+q7oXBUfGNDC6ErgG+ALxigvs5F/gYg6PRJwLfAJ672vuZs69L+P8Xd7/V/ZzWdtNLXbNY7f0cBlwLnDN0Xg87cMhwOZ3BhZ87gHd1y/4YeGU3fUTX1O3AV4Gnz1n3Xd16twIvX2qb3fKnd9u4vdvm4SvdxyT2A7ybwXnXnXMex05qP/PqPJ/lb+dc1b0Av8Hg9rpbgPdN+O/ak7rluxiE/tsnpJ/nMziy/gGDmx92zVnnt7vxtwO/Ncn9dL9rP+KxWbBxqX785K4kNcZP7kpSYwx+SWqMwS9JjTH4JakxBr8kNcbgl6TGGPyS1BiDX5Ia878uEEGzGqLXfAAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x10c84add8>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "%matplotlib inline\n",
    "rand = numpy.random.RandomState(seed=1234567)\n",
    "\n",
    "def gen_data(events, *times):\n",
    "    x = rand.normal(times[0], times[0]/5, events)\n",
    "    \n",
    "    d = 0\n",
    "    for i in range(1, len(times)):\n",
    "        t = times[i]\n",
    "        num_anomalous_events = int(rand.normal(5, 1))\n",
    "        for e in rand.normal(t, t/5, num_anomalous_events):\n",
    "            x[d] = e\n",
    "            d += 1\n",
    "\n",
    "    rand.shuffle(x)\n",
    "    return x\n",
    "\n",
    "x = gen_data(100, 1e-7, 1e-5)\n",
    "plot_data = matplotlib.pyplot.hist(x, 500)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Variables:\n",
    "1. Number of points\n",
    "2. Number of anomalies\n",
    "3. Number of necessary starting data points\n",
    "4. Frequency/Distribution of anomalies (e.g. what if a bunch happen close together vs far apart?)\n",
    "5. Frequency/Distribution of data points in time\n",
    "\n",
    "Experiment:\n",
    "- Run algorithm on a set of data N times\n",
    "- Measure the success of each run according to criteria below\n",
    "- Aggregate that into a single number we can use to compare algorithms\n",
    "\n",
    "Analysis Criteria:\n",
    "- some sort of weighted sum of True Positive, False Positive, True Negative, and False Negative\n",
    "- some sort of overhead calculation\n",
    "\n",
    "Equipment Needed:\n",
    "1. Function to generate data with anomalies in it that looks kind of realistic\n",
    "    - Scrape durations from existing data\n",
    "    - Add anomalies into it to form a composite data set\n",
    "2. Function to calculate success criteria from a given run\n",
    "3. Some algorithms"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Sanity Check\n",
    "\n",
    "Check that method makes sense; for constant-size changes to the data over time, differences made when the baseline is smaller should show up more prominently than when the baseline is larger"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "\n",
    "p1 = numpy.array([1, 2, 3, 4, 5, 6, 7, 8, 9, 10])\n",
    "p2 = numpy.array([1, 2, 4, 4, 5, 6, 7, 8, 9, 10])\n",
    "p3 = numpy.array([1, 2, 3, 4, 5, 6, 8, 8, 9, 10])\n",
    "\n",
    "matplotlib.pyplot.plot(numpy.cumsum(numpy.abs(numpy.diff(p2/p1))))\n",
    "matplotlib.pyplot.plot(numpy.cumsum(numpy.abs(numpy.diff(p2)/numpy.diff(p1))))\n",
    "#matplotlib.pyplot.plot(numpy.cumsum(numpy.abs(numpy.diff(p3/p1))))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Come up with a way to calculate bin width and kernel width"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 185,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "63\n",
      "63\n",
      "7.937253933193772\n",
      "\n",
      "1\n",
      "1\n",
      "1\n",
      "1\n",
      "1\n",
      "1\n",
      "1\n",
      "1\n",
      "1\n",
      "2\n",
      "2\n",
      "2\n",
      "2\n",
      "2\n",
      "2\n",
      "2\n",
      "2\n",
      "2\n",
      "2\n",
      "2\n",
      "2\n",
      "2\n",
      "2\n",
      "2\n",
      "2\n",
      "3\n",
      "3\n",
      "3\n",
      "3\n",
      "4\n",
      "6\n",
      "9\n",
      "\n",
      "\n",
      "IQR delta 3.49875746223e-09\n",
      "mean of deltas 9.29215313882e-09\n",
      "median of deltas 7.69404377411e-09\n",
      "sample std 2.69101820146e-08\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZYAAAD8CAYAAABU4IIeAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvNQv5yAAAIABJREFUeJzt3XucXHV9//HXZ2bvl9wTLiEhgYRLgogYQBEvBStRqqEtarwVEIsolFZLC7RVAcE29qGx/gSVCphCMSCoBAWpBarYQm4kBJKQZEnAhNxJNsneZnZmPr8/5pvNZDMzOwmzOzO77+fjsY85853v+Xw/50wynznne2bG3B0REZFiiZQ6ARERGVxUWEREpKhUWEREpKhUWEREpKhUWEREpKhUWEREpKhUWEREpKhUWEREpKhUWEREpKiqSp1AKYwZM8YnTZpU6jRERCrK0qVLd7r72L76DcnCMmnSJJYsWVLqNEREKoqZvVZIP50KExGRolJhERGRolJhERGRolJhERGRolJhERGRolJhERGRolJhERGRolJhkf73+A3pvwxzFs1hzqI5Pfe3fuMbbP3GNw7q88yDa3nmwbUHtT394zt5+sd3HtTW+ugrtD76ysFDPv44jz/++EFta9d+nbVrv95z/yvrNvGVdZsO6nPzoyu5+dGVh51/odtQzPwL3oYi5l/oNhQ1/8PYBsnv7t9v4FcrtvT7OEPyA5IywLa+eEjTy7tePuh+bPXLh/TZubHtkLbtr60/pC2+uf3QIbduPaRtX9vqg+6/1NZ5SJ9Vm/ce0lZI/lDYNhQzfyhwG4qYPxS2DUXNHwreBsnv3ude47Txw7no9GP6dRwdsYiIDBHJlBO1/h9HhUVEZIhIppxIpP8riwqLiMgQkXInaiosIiJSJMmUE9URi4iIFEvKdSpMRESKKD15r8IiIiJFolNhIiJSVCmHiI5YRESkWNJHLP0/TkFDmNlMM1tjZi1mdkOWx2vN7IHw+EIzm5Tx2I2hfY2ZXdhXTDObHGKsCzFrChjjdDN71sxWmtmLZlZ3JDtDRGQwS5bL5L2ZRYHbgQ8C04BPmNm0Xt2uAHa7+xRgLjAnrDsNmA1MB2YCd5hZtI+Yc4C57j4V2B1i5xujCrgPuMrdpwPvA7oPcz+IiAx6qTKavD8baHH39e4eB+YDs3r1mQXMC8sPAReYmYX2+e4ec/cNQEuIlzVmWOf8EIMQ8+I+xvgAsMLdXwBw9zfcPVn4LhARGRoSZTR5Px7YmHF/U2jL2sfdE8AeYHSedXO1jwZaQ4zeY+Ua4yTAzewJM3vezP6+gG0SERlSurrT77frqqP9PlYh326crbx5gX1ytWcraPn65xujCjgPOAvoAJ40s6Xu/uRBCZpdCVwJMHHixCyhREQGrz2d6RmC4fXV/T5WIUcsm4AJGfePAzbn6hPmPIYDu/Ksm6t9JzAixOg9Vr4xfuvuO929A3gMOLP3Rrj7ne4+w91njB07toDNFhEZPHZ3xAEY0VAehWUxMDVcrVVDejJ+Qa8+C4BLw/IlwFPu7qF9driiazIwFViUK2ZY5+kQgxDzkT7GeAI43cwaQsF5L7Cq8F0gIjL4vdGWLixjmmr7faw+T4W5e8LMriH9Ah4F7nb3lWZ2C7DE3RcAdwH3mlkL6aOI2WHdlWb2IOkX+gRw9f6J9Wwxw5DXA/PN7FZgWYhNnjF2m9m3SRcrBx5z91+9qb0iIjLI7GyLATCmqabfxyroFyTd/THSp5gy276asdwFfDTHurcBtxUSM7SvJ33VWO/2fGPcR/qSYxERyWLnAB6x6JP3IiJDwBttMaoixrC68phjERGRCvdGW5xRjTXl8cl7ERGpfDvbYgNyGgxUWEREhoSd7XFGD8DEPaiwiIgMCTv36YhFRESKpDOeZMueTiaMahiQ8VRYREQGudVb95JymH7ssAEZT4VFRGSQW/n6HkCFRUREimTl5r0Mr69m/Ij6ARlPhUVEZJBbs20fpxzdjA3Aj3yBCouIyKDm7rRsa+Oko5oHbEwVFhGRQWzb3hj7YglOOqppwMZUYRERGcTWbd8HwJRxOmIREZEiWLetDYCpOmIREZFiWL6xlTFNNQP2qXtQYRERGbQ640meXL2NP5521ICOq8IiIjJI/XTpRtrjST58+rEDOq4Ki4jIIPXLFVsY3VjDOSeMHtBxVVhERAahza2dLNqwi0vPnUR0AH7cK5MKi4jIIPS1BSsB+MhbB/Y0GKiwiIgMOqmUs+wPrZxydDOTxjQO+PgqLCIig0zLjjZ2tsW4/F2TSjJ+QYXFzGaa2RozazGzG7I8XmtmD4THF5rZpIzHbgzta8zswr5imtnkEGNdiFmTbwwzm2RmnWa2PPz94Eh3hojIYPCLZa9TFTHed/K4kozfZ2ExsyhwO/BBYBrwCTOb1qvbFcBud58CzAXmhHWnAbOB6cBM4A4zi/YRcw4w192nArtD7JxjBK+4+xnh76rD2gMiIoNIPJHikeWbmXbsMI4aVleSHAo5YjkbaHH39e4eB+YDs3r1mQXMC8sPARdY+vuZZwHz3T3m7huAlhAva8ywzvkhBiHmxX2MISIiwWMvbuH11k6+8N4TS5ZDIYVlPLAx4/6m0Ja1j7sngD3A6Dzr5mofDbSGGL3HyjUGwGQzW2ZmvzWzdxewTSIig04q5dz1+w001kS5cPrRJcujqoA+2Y4KvMA+udqzFbR8/fONsQWY6O5vmNnbgV+Y2XR333tQgmZXAlcCTJw4MUsoEZHKlUw5//bkOl58fQ83fXgakQH+7EqmQo5YNgETMu4fB2zO1cfMqoDhwK486+Zq3wmMCDF6j5V1jHCa7Q0Ad18KvAKc1Hsj3P1Od5/h7jPGjh1bwGaLiFSOO3+3nu8+uY6L3nIMn37H8SXNpZDCshiYGq7WqiE9Gb+gV58FwKVh+RLgKXf30D47XNE1GZgKLMoVM6zzdIhBiPlIvjHMbGy4GAAzOyGMsb7wXSAiUvkefWEzb50wgu998m1URUv7SZI+T4W5e8LMrgGeAKLA3e6+0sxuAZa4+wLgLuBeM2shfaQyO6y70sweBFYBCeBqd08CZIsZhrwemG9mtwLLQmxyjQG8B7jFzBJAErjK3Xcd+S4REaksKzfvYdWWvfz9zJMH7Hft8ylkjgV3fwx4rFfbVzOWu4CP5lj3NuC2QmKG9vWkrxrr3Z51DHd/GHi4z40QERmEdrbF+PSPFjKuuZaPzZjQ9woDoKDCIiIi5cfd+dIDy9nd0c39nztnQH/MKx99pYuISIV6/KWtPLNuJ//4oVM5d8qYUqfTQ4VFRKQCPbR0E19+cDmnHjOsZN8JlotOhYmIVJDOeJKbH13J/MUbOffE0Xxn9hklvwqsNxUWEZEKcssvV/LAko184X0n8qX3n0RNVXkVFVBhERGpGL9du4OfLNrIX757MtfPPKXU6eSkwiIiUgEeXrqJmx9dyYRR9fztB04udTp5qbCIiJSxN9pi/PPjL/PQ0k2cNWkk//xnp1NXHS11WnmpsIiIlKktezr5zF2LWL+jjc+/5wSuu/Bkqstsoj4bFRYRkTK0Zus+Pn3XQjpiCe6+7KyS/RrkkVBhEREpM8+tf4Mv3LeU6miEX1z9LqYe1VzqlA5L+R9TiYgMId//n1f45L8/hwP3fe6ciisqoCMWEZGy8cTKrcz59cuce+JofvCZtzOsrrrUKR0RFRYRkTLwfy07ufYnyzjpqCbu/IsZNNVW7suzToWJiJTY79bu4Or7n+eoYXXc97lzKrqogAqLiEhJPf3ydi7/8WKG1Vdz16UzGNdcV+qU3rTKLosiIhUqlXJu/dVq5j37Kqce08x/fPYcRjXWlDqtolBhEREZYNv3dXHdT1fwu7U7+NiM4/inP5lWsRP12aiwiIgMoMWv7uJLDyxn+94YN314Gpe9a3KpUyo6FRYRkQHw6s52vvPfa3nkhc2MbqzhwaveyRkTRpQ6rX6hwiIi0s9Wbt7DJ/99IfFEis+dN5lrzp/K8PrBc+qrNxUWEZF+tOCFzdzw8Aqaaqt4+AvvZMq4yvsk/eEq6HJjM5tpZmvMrMXMbsjyeK2ZPRAeX2hmkzIeuzG0rzGzC/uKaWaTQ4x1IWZNX2OExyeaWZuZXXe4O0FEpNjaYwm+9shLXPuTZZx6zDAWXHPekCgqUEBhMbMocDvwQWAa8Akzm9ar2xXAbnefAswF5oR1pwGzgenATOAOM4v2EXMOMNfdpwK7Q+ycY2SYCzxe6IaLiPSX1o44n5u3hHnPvsZn3zWZ+Ve+g6OHV/7nUwpVyBHL2UCLu6939zgwH5jVq88sYF5Yfgi4wMwstM9395i7bwBaQrysMcM654cYhJgX9zEGZnYxsB5YWfimi4gU3/+17OTdc57m2fVv8HcXnsxXPzytIn5DpZgK2drxwMaM+5tCW9Y+7p4A9gCj86ybq3000Bpi9B4r6xhm1ghcD9ycbyPM7EozW2JmS3bs2NHHJouIHL5fv7SVz85bzNhhtfzq2vO4+o+mlDqlkiiksFiWNi+wT7Ha841xM+lTZ21ZHj/Q0f1Od5/h7jPGjh2br6uIyGFJpZx/+sWLXHXfUqaOa+a+K85h+rHDS51WyRRyVdgmYELG/eOAzTn6bDKzKmA4sKuPdbO17wRGmFlVOCrJ7J9rjHOAS8zsm8AIIGVmXe7+vQK2TUTkTWnZ3sbNj67kmXU7uezcSfzDh06lpmponfrqrZDCshiYamaTgddJT8Z/slefBcClwLPAJcBT7u5mtgC438y+DRwLTAUWkT76OCRmWOfpEGN+iPlIvjGAd+9PwsxuAtpUVERkIPx+3U6uum8pZnDzR6bzF+88njD1O6T1WVjcPWFm1wBPAFHgbndfaWa3AEvcfQFwF3CvmbWQPoqYHdZdaWYPAquABHC1uycBssUMQ14PzDezW4FlITa5xhARGWjuzpOrt3PVfUs5engdd192FidV4C899peCPiDp7o8Bj/Vq+2rGchfw0Rzr3gbcVkjM0L6e9FVjvdtzjpHR56Z8j4uIvFkt2/dxzf3LeHnrPk45upn7//Idg+ZbiYtFn7wXESnQ4y9u4e8fWkFV1PjGn76FPztzPHXV0VKnVXZUWERE+vBGW4yv/3IVv1i+menHDuOOT53J8aMbS51W2VJhERHJY8WmVq6Yt4Td7XGu+aMpXP1HU6iv0VFKPiosIiI5vPZGO9f+ZBkdsQQ/++K5nH7c4Pya+2JTYRERyWLpa7v4xJ0LqY4aP7r0LBWVw6DCIiLSy+bWTi67ezHHjkhfSnzC2KZSp1RRhvbHQ0VEetndHufyexYD8KNLVVSOhI5YRESCddv2cdk9i9mxL8Y9l5/FlHEqKkdChUVEBNi0u4NL715EPOmD+vfoB4IKi4gMeet3tPGxHz5LLJHiJ3/5Dk4bP3S/mbgYVFhEZEiLJ1JcO38ZiZTzsy+cy1R959ebpsIiIkOWu/OPP3+Rl17fyw8+/XYVlSLRVWEiMmR998kWfrp0E19834nMPO3oUqczaKiwiMiQdOfvXmHuf6/lglPGcd0HTi51OoOKToWJyJDz7f9aw3efauHcE0fz/U+/nUhEP85VTDpiEZEh5Xdrd/Ddp1p4/6lHcdelZw35nxHuD9qjIjJkrN6yl7/6yTKmjGviWx97q76luJ+osIjIkLBpdwd//v3/ozoa4fufOpPh9dWlTmnQ0hyLiAx67s5tv1pNdzLFT696J5PH6Ee6+pOOWERk0Jv73+t4/KWtfOF9U1RUBoAKi4gMao8sf53vPrmOWWccy19fMLXU6QwJBRUWM5tpZmvMrMXMbsjyeK2ZPRAeX2hmkzIeuzG0rzGzC/uKaWaTQ4x1IWZNvjHM7GwzWx7+XjCzPz3SnSEig0sskeR7T7VwwthGvvXRtxLVZcUDos/CYmZR4Hbgg8A04BNmNq1XtyuA3e4+BZgLzAnrTgNmA9OBmcAdZhbtI+YcYK67TwV2h9g5xwBeAma4+xlhjB+ameaORIRHX9jCuu1t/N0HTqYqqhM0A6WQPX020OLu6909DswHZvXqMwuYF5YfAi4wMwvt89095u4bgJYQL2vMsM75IQYh5sX5xnD3DndPhPY6wAvdeBEZvPZ0dnPLoys5elgdF5x6VKnTGVIKKSzjgY0Z9zeFtqx9wov8HmB0nnVztY8GWjMKReZYucbAzM4xs5XAi8BVGeuLyBD18NJN7O1K8L1Pvk0fghxgheztbCclex8V5OpTrPa8ebj7QnefDpwF3Ghmdb07mtmVZrbEzJbs2LEjSygRGSw27urg279Zy1snjGDGpFGlTmfIKaSwbAImZNw/Dticq0+Y3xgO7Mqzbq72ncCIjDmSzLFyjdHD3VcD7cBpvTfC3e909xnuPmPs2LF9brSIVCZ35/P3LiXlzr9ecnqp0xmSCiksi4Gp4WqtGtKT8Qt69VkAXBqWLwGecncP7bPDFV2TganAolwxwzpPhxiEmI/kGyPEqAIws+OBk4FXC94DIjKo3LfwD6zasper3nsiJ+n3VUqiz6un3D1hZtcATwBR4G53X2lmtwBL3H0BcBdwr5m1kD6KmB3WXWlmDwKrgARwtbsnAbLFDENeD8w3s1uBZSE2ucYAzgNuMLNuIAV80d13HvkuEZFK9uTqbRwzvI6r/2hKqVMZsgq6LNfdHwMe69X21YzlLuCjOda9DbitkJihfT3pq8Z6t2cdw93vBe7tcyNEZNB7dWc7v127g7989wn6zEoJ6VIJERk07nvuNaoixufOm1zqVIY0FRYRGTSefHk7b5swknHDDrkwVAaQCouIDAp/eKODDTvb+dBb9Nv1pabCIiKDwuJX058+OOeE0SXORFRYRGRQeGbdDkY31nCyLjEuORUWERkUXti0h7MmjSKiq8FKToVFRCre9n1dbNjZzklH62ilHKiwiEjFe/zFrQCauC8TKiwiUvGW/WE3Rw+r45Sjh5U6FUGFRUQGgef/0Mpbjhte6jQkUGERkYq2p6ObP+zq4MyJI0udigQqLCJS0V7Y1ArAaeN1GqxcqLCISEV7afMeAE4/bkSJM5H9VFhEpKI9/1orE0bVM7y+utSpSKDCIiIVbXdHnImjGkqdhmRQYRGRitYZT1JfHS11GpJBhUVEKlpXIkmtCktZUWERkYrWFU9SV6XCUk5UWESkou3u6GZkgybuy4kKi4hUrGTK6exOMra5ttSpSAYVFhGpWN3JFABjmlRYyokKi4hUrO6kAzBGRyxlpaDCYmYzzWyNmbWY2Q1ZHq81swfC4wvNbFLGYzeG9jVmdmFfMc1scoixLsSsyTeGmf2xmS01sxfD7flHujNEpLLsP2IZqyOWstJnYTGzKHA78EFgGvAJM5vWq9sVwG53nwLMBeaEdacBs4HpwEzgDjOL9hFzDjDX3acCu0PsnGMAO4EPu/tbgEuBew9vF4hIpTpwKqymxJlIpkKOWM4GWtx9vbvHgfnArF59ZgHzwvJDwAVmZqF9vrvH3H0D0BLiZY0Z1jk/xCDEvDjfGO6+zN03h/aVQJ2Z6e2LyBCQ8vSpsIbaqhJnIpkKKSzjgY0Z9zeFtqx93D0B7AFG51k3V/tooDXE6D1WrjEy/TmwzN1jBWyXiFS4VPqAhboqTReXk0LKvGVp8wL75GrP9q8gX/8+8zCz6aRPj30gSz/M7ErgSoCJEydm6yIiFSblTjRiVEVVWMpJIc/GJmBCxv3jgM25+phZFTAc2JVn3VztO4ERIUbvsXKNgZkdB/wc+At3fyXbRrj7ne4+w91njB07toDNFpFyl3IdrZSjQp6RxcDUcLVWDenJ+AW9+iwgPXEOcAnwlLt7aJ8druiaDEwFFuWKGdZ5OsQgxHwk3xhmNgL4FXCju//v4Wy8iFS2lDt1+p6wstNnYQnzGdcATwCrgQfdfaWZ3WJmHwnd7gJGm1kL8GXghrDuSuBBYBXwa+Bqd0/mihliXQ98OcQaHWLnHCPEmQJ8xcyWh79xR7g/RKSCpNyp1RFL2SnoUgp3fwx4rFfbVzOWu4CP5lj3NuC2QmKG9vWkrxrr3Z51DHe/Fbi1z40QkUEn5eiIpQyp1ItIxUq5U6MjlrKjZ0REKpbmWMqTCouIVCx3NMdShvSMiEjFSqV0xFKOVFhEpGI5UB3N9tlpKSUVFhGpWOmv3lBhKTcqLCJSuRxMdaXsqLCISEVTXSk/KiwiUtF0xFJ+VFhEpII5pmOWsqPCIiIVy9ERSzlSYRGRiqbCUn5UWESkoulUWPlRYRGRyqa6UnZUWESkcrnqSjlSYRGRipWevFdpKTcqLCJS0VRWyo8Ki4hUNB2wlB8VFhGpaKor5UeFRUQqluZYypMKi4hULn27cVlSYRGRiqYPSJafggqLmc00szVm1mJmN2R5vNbMHgiPLzSzSRmP3Rja15jZhX3FNLPJIca6ELMm3xhmNtrMnjazNjP73pHuCBGpRK4jljLUZ2ExsyhwO/BBYBrwCTOb1qvbFcBud58CzAXmhHWnAbOB6cBM4A4zi/YRcw4w192nArtD7JxjAF3AV4DrDnPbRaTCOZq8L0eFHLGcDbS4+3p3jwPzgVm9+swC5oXlh4ALLD2jNguY7+4xd98AtIR4WWOGdc4PMQgxL843hru3u/vvSRcYERlidMRSfgopLOOBjRn3N4W2rH3cPQHsAUbnWTdX+2igNcToPVauMQpiZlea2RIzW7Jjx45CVxORMqc5lvJTSGHJ9qx5gX2K1V5oHjm5+53uPsPdZ4wdO7bQ1USkzOmIpfwUUlg2ARMy7h8HbM7Vx8yqgOHArjzr5mrfCYwIMXqPlWsMERmi9ENf5amQwrIYmBqu1qohPRm/oFefBcClYfkS4Cl399A+O1zRNRmYCizKFTOs83SIQYj5SB9jiMhQ5aDp+/JT1VcHd0+Y2TXAE0AUuNvdV5rZLcASd18A3AXca2YtpI8iZod1V5rZg8AqIAFc7e5JgGwxw5DXA/PN7FZgWYhNrjFCrFeBYUCNmV0MfMDdVx3pThGRyqEjlvLTZ2EBcPfHgMd6tX01Y7kL+GiOdW8DbiskZmhfT/qqsd7t+caYlHcDRGTQUl0pP/rkvYhUNB2xlB8VFhGpWOlLSVVZyo0Ki4hUNB2xlB8VFhGpaKor5UeFRUQql7t+j6UMqbCIiEhRqbCISMVyIKIjlrKjwiIiFU11pfyosIhIRVNdKT8FffJeRKQcDbUvoXR3vLOTVHs7qY6Og2+ztYXbZHs73tFBsr2dxnPewbgvf6lf81RhEZGK1B5LkEo5E0Y1lDqVnDyZzF8A2vMUhRyFggK/e9dqaog0NhJpaEjfNjYSbWomOqy5n7dahUVEKoTjbN/bxQsbW/l9y05e3rqPiBnvPal4v6+UiscPfsHvyFzOcZunzbsK/2HbSEMD1thAtKGx5zY6ZjQ1x088uEBkFIqcbQ0NWHV10fbL4VJhEZGykUo53ckU8WSKBS9s5rWd7by2q4PLNu+hI57kY994EoD66igjT4hyQnMd41MdxF/dUdBpoP23h/brSB8NdHcXlmg0mvUFvXrkyNDWQKShMX3bZ1FoJNJQj0UGz5S3CouI9Ju2WIJdbXFaO+PsaYuT6kzy8yeep2PPXjpa22iIt1Ld2cnX5z1KZ+s+utvbmBaPUZ+I0fLL9O2ZJIgk9jAsCT+LzqW2u4torItURzvmsI739JlHz2mhjBf0aFMz1Ucd3fcLf+/i0NiI1dTog5l5qLCIyCHcnVQ8jnXuxWMdrHtuOXt2ttK5tw1b+wesu5tffO07JNraSLV3MHXLHjwe4+kFXURjnRzf9QZV3d0s+n93UZ/ooj4R45hUIkT/VM845+UYv6u6lrrmJqqam4g2NhDdV0WkJkLzyaf3vMj/csuTxOuizH7bZflPEZX4tNBQpMIiMgi4O96RPp0zcvc2qmOdvPZ0kq49+xixYg/JWIoN275Hoq2NZHs75+5YR3UsyeM//AyRzk6isU7G7tiEpVIsOf3t1HZ3UeUpTg7xEz+FRtJ/Y/YPuuIlAJIWwarq8Opajq2vJ1lXj1dXEW9soOqo40k2NxFvbiLuT2H1UUYd/1kaRwyjYcQwbtq8m1h9Pd9+28npo4jGRj5573LcIjzw+Xce2MB7LkrfXv7tnqbf/3otAF+Y+Yn+3LVyBFRYRAaYO7S1d9G+Zx+drftIbX2NVLyT9Q9tIb5vH4m2dhpeXEw0HufXL23EOzqgs4OR8RasK8mT7X9GNNZJVayL65Jd1MW6WBWLEUn/Ti9/Hcbp+En6dle4P4rbAYhFqmiocbpqInRHNxOvqaertp4R1TUQifLqae/Awjv9yOspUjV1NJ87geaRw2gY3szqX94PNVVcdP3XqB8xDKupYcedLwJw2udPB+Cee+4B4PLLL+/Z7qXPrwHgxDM/3dO2Ydk6AGonTz6wf2zwzDUMVSosMuS5O91JpyuRpDuZIuWwbts+urpTdCWSTOmIkUrAc8+spnvfPpLt7dS++gY1sQT/ueKeMAHcySnrXoVUimc+dS3RrvRRwJjWDqoSMZ74aRe13TGmxvZSnUyycf4jh+Qx/OEDy8f0LP0vHVW1dFXVUN3YSaouQluym+6aZpLNY9k2qpGuunqOpT5dDBobWLI9Rnd1HReeNZnq5mbetuoWorURtsyaR+OIYYxtrOO6Z67CDO6ZeU/PSK995i8AmDHvOz1tP//W8wC8/6oze9pWLUz/8GvjUcW7GksGFxUWKXvxRIr2WIK2WIKOeJJkyvnt2h10xhN0dieZuq+LVMp58Kl1dHYn6ejqxrbsJhKLc8O//QrvaMc7Ozi3aQuRriQ/u+NWol2dnDwC6mJdfGvTbuq60xPFlyTj1Ce6eOXHceoT6bYtiU7M4QT+rCenaT1L3zwoVwdOfel/idfUkaipI5qsIVlVi40bRVddA1071pOoqiJ62ruINjYSaWpk/PoYkdp62t43mdphTdQOa+aZ3/2GVG0tsy//LHU1VVRFIyx9/pNEgFln3t8z3p+Gd/w/f9vUnrZf/PBZAC66LJxK2vstAEadPKGnj+adpT+psEi/6upOEkmmSKactRtb2dvVzb6uBNv3xUimnG/++mVa93VywdY9kEz9YriHAAALaUlEQVQy57aHSLSlPz8wY1uMSHcXf7PoYeoTMRoSMU7bs5rqZJJVi67taetOdBNJxDjj/vQkcV0i3nNaKJup/GfPcry6Gq+pJ1VXT6q+ga0pI9bYzOjxo7GG9ORv87bfEK2JsHfG56huaqKmuZF/33AXibooXzn/NqJNTUQaGtj0V9eCGdPuu7cn/v53/H/6t+l3/A/cfAMAH//aP/X02f7DFQCcPOv0nrbfr1hIFGiqr3nzT4LIAFNhGURSKac7lSKRTH8WoDvpJFIpuhO921MkUhl9wm26Pb0cT6TSf8kUse4U8WQyfZtIkuiMkepoh84OvKMD6+rEOjqxzg7o6iDS2UE01kW0q5PaeBcz/XUs4ay477Keo4Dr2EF93KmLf5TanquF4K/u/coh2/XujOVENEqyqoqTx1RjwxqINI4g2upYTR1j3nIsVU1NRBobWLpqFcmaGs57//t7rg5au/Ff8LoIbznrDiKNjXxszeukqqoOerf/8fBu/+CJ4+UAHHv553qaNv76ZwDUn3JKT9tg+hyCyJtRUGExs5nAvwFR4Efu/i+9Hq8F/gN4O/AG8HF3fzU8diNwBZAErnX3J/LFNLPJwHxgFPA88Bl3jx/JGAOpM55kd0ecXe1xWju62dvVTVd3Mn2evjtJPJl+oY4lkj0v4t3JAy/6mS/wiV5FoLvXciKZojvVuyA4ydSh79LNU9Qmu3ve3df1ut3/Dr8h0RVu04/t7z8yub9fPPSNUeWpgvZJKhIhWVtPTTRFpNo4fUyESOMYqhobWNEZJ1EX5f2n/AmRpkb2/PwXWDTKmL+6pufqoP/52R9IVddz4TUziDY2YvX1PPj1fwDg41878E9w/zv+cZ8/8I5/S5g8Hv6RjxzIp7YRgOpj0jMYqapth/MUi0iB+iwsZhYFbgf+GNgELDazBe6+KqPbFcBud59iZrOBOcDHzWwaMBuYDhwL/LeZnRTWyRVzDjDX3eeb2Q9C7O8f7hjunnwzOyaXJ1dv47EXt7K5tZPdHekisrsjTixR2IttxKCmKkJ1dP+fURWJUFMVoSpi1JjTlIzRkIwzPBmjMXngXH/6xb2L2kSc2ngXtd1d1MRj1MQ7qY7HqI53Uh3rJBrroircRmJdWIHfLURNLTQ0pCeBmxqINDYRbTyKqqb0C33BHyQLbT0fIuu5VPRnPUN989fpq4U+NfPvAOh49jkAhl90UU+f2O+iAFSPG1dY/iJSFgo5YjkbaHH39QBmNh+YBWQWllnATWH5IeB7lv5Y6ixgvrvHgA1m1hLikS2mma0Gzgc+GfrMC3G/fwRjPFvgPijYc+vf4Ip5SxjZUM2JY5s4bkQ9Z4yrZ2xVklGWZKR1M5xumr2bxmSMmngXVeGUUKQrfarIOzpI7cn9nUMejxeczyEv8iMbiTSMOfAi33Doi33OolBfrw+RiUhRFFJYxgMbM+5vAs7J1cfdE2a2Bxgd2p/rte74sJwt5mig1d0TWfofyRhFNb1zOw8v/DeaU/Ger64mmfvAKB7+gAPfLdTrRb961Kjc3ye0/6snGhuxhoaDb+sH13cLicjgUUhhyXZhYu9zK7n65GrP9oqYr/+RjHFwgmZXAlcCTJw4Mcsqfasd1sRRp08r/FtGMwrGkP5uoaPfckjTKaNOOeh+7amnHNJnzISmQ9rGHX/CIW01xzYeOuTRRx/S1tx06kH3T2uqP6TPtGOHHdJWSP5Q2DYUM38ocBuKmD8Utg1FzR8K3gYpD+Z9nH83s3cCN7n7heH+jQDu/s8ZfZ4IfZ41sypgKzAWuCGz7/5+YbVDYgL/AuwAjg5HJT1jH+4Y7p7zVNiMGTN8yZIlhewfEREJzGypu8/oq18h51IWA1PNbLKZ1ZCeKF/Qq88C4NKwfAnwlKcr1gJgtpnVhqu9pgKLcsUM6zwdYhBiPnKEY4iISAn0eSosHDlcAzxB+tLgu919pZndAixx9wXAXcC9YeJ8F+lCQej3IOmJ/gRw9f6rtbLFDENeD8w3s1uBZSE2RzKGiIgMvD5PhQ1GOhUmInL4inkqTEREpGAqLCIiUlQqLCIiUlQqLCIiUlQqLCIiUlRD8qowM9sBvPYmQowBdhYpnWJSXodHeR0e5XV4BmNex7t7nz8dOiQLy5tlZksKueRuoCmvw6O8Do/yOjxDOS+dChMRkaJSYRERkaJSYTkyd5Y6gRyU1+FRXodHeR2eIZuX5lhERKSodMQiIiJFpcJyGMxsppmtMbMWM7uhBOO/amYvmtlyM1sS2kaZ2W/MbF24HRnazcy+G3JdYWZnFjGPu81su5m9lNF22HmY2aWh/zozuzTbWEXI6yYzez3ss+Vm9qGMx24Mea0xswsz2ov6PJvZBDN72sxWm9lKM/vr0F7SfZYnr5LuMzOrM7NFZvZCyOvm0D7ZzBaGbX/A0j+5gaV/MuOBMPZCM5vUV75FzuvHZrYhY3+dEdoH7N9+iBk1s2Vm9stwv3T7y931V8Af6a/3fwU4AagBXgCmDXAOrwJjerV9E7ghLN8AzAnLHwIeJ/0Lm+8AFhYxj/cAZwIvHWkewChgfbgdGZZH9kNeNwHXZek7LTyHtcDk8NxG++N5Bo4BzgzLzcDaMH5J91mevEq6z8J2N4XlamBh2A8PArND+w+AL4TlLwI/CMuzgQfy5dsPef0YuCRL/wH7tx/ifhm4H/hluF+y/aUjlsKdDbS4+3p3jwPzgVklzgnSOcwLy/OAizPa/8PTngNGmNkxxRjQ3X9H+jdx3kweFwK/cfdd7r4b+A0wsx/yymUWMN/dY+6+AWgh/RwX/Xl29y3u/nxY3gesBsZT4n2WJ69cBmSfhe1uC3erw58D5wMPhfbe+2v/fnwIuMDMLE++xc4rlwH7t29mxwEXAT8K940S7i8VlsKNBzZm3N9E/v+E/cGB/zKzpWZ2ZWg7yt23QPqFAhgX2gc638PNYyDzuyacirh7/+mmUuUVTju8jfS73bLZZ73yghLvs3BaZzmwnfQL7ytAq7snsozRM354fA8weiDycvf9++u2sL/mmllt77x6jd8fz+N3gL8HUuH+aEq4v1RYCmdZ2gb6krp3ufuZwAeBq83sPXn6lkO+kDuPgcrv+8CJwBnAFuBbpcrLzJqAh4G/cfe9+boOZG5Z8ir5PnP3pLufARxH+l3zqXnGKFleZnYacCNwCnAW6dNb1w9kXmb2J8B2d1+a2ZxnjH7PS4WlcJuACRn3jwM2D2QC7r453G4Hfk76P9y2/ae4wu320H2g8z3cPAYkP3ffFl4MUsC/c+DQfkDzMrNq0i/e/+nuPwvNJd9n2fIql30WcmkF/of0HMUIM9v/c+qZY/SMHx4fTvqU6EDkNTOcUnR3jwH3MPD7613AR8zsVdKnIc8nfQRTuv31ZiaLhtIfUEV6km0yByYopw/g+I1Ac8by/5E+L/uvHDwB/M2wfBEHTxwuKnI+kzh4kvyw8iD9zm4D6cnLkWF5VD/kdUzG8pdIn0MGmM7BE5XrSU9CF/15Dtv+H8B3erWXdJ/lyauk+wwYC4wIy/XAM8CfAD/l4MnoL4blqzl4MvrBfPn2Q17HZOzP7wD/Uop/+yH2+zgweV+y/VW0F5qh8Ef6Ko+1pM/3/uMAj31CeNJfAFbuH5/0udEngXXhdlRoN+D2kOuLwIwi5vIT0qdIukm/y7niSPIAPkt6grAFuLyf8ro3jLsCWMDBL5r/GPJaA3ywv55n4DzSpxRWAMvD34dKvc/y5FXSfQacDiwL478EfDXj/8CisO0/BWpDe1243xIeP6GvfIuc11Nhf70E3MeBK8cG7N9+Rtz3caCwlGx/6ZP3IiJSVJpjERGRolJhERGRolJhERGRolJhERGRolJhERGRolJhERGRolJhERGRolJhERGRovr/UI3ZsuGb+LMAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x1a1ec8b668>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "seed_points = 1000\n",
    "\n",
    "seeds = numpy.array(\n",
    "        list(gen_data(seed_points, 1e-7))\n",
    "      + list(gen_data(seed_points, 2e-7))\n",
    "      + list(gen_data(seed_points, 3e-7))\n",
    "      + list(gen_data(seed_points, 4e-6)))\n",
    "\n",
    "ls = len(seeds)\n",
    "\n",
    "sorted_durations = sorted(seeds)\n",
    "matplotlib.pyplot.plot(sorted_durations)\n",
    "\n",
    "lsd = len(sorted_durations)\n",
    "num_deltas = int(lsd ** 0.5)\n",
    "points_per_bin = int(len(sorted_durations) / num_deltas)\n",
    "\n",
    "print(num_deltas)\n",
    "print(points_per_bin)\n",
    "print(points_per_bin ** 0.5)\n",
    "\n",
    "start = int(num_deltas / 4)\n",
    "finish = int(3 * num_deltas / 4) + 1\n",
    "\n",
    "mp = sorted_durations[-1]\n",
    "last_marker = sorted_durations[int(start * lsd / num_deltas)]\n",
    "\n",
    "deltas = []\n",
    "for i in range(start+1, finish):\n",
    "    x = int(lsd * i / num_deltas)\n",
    "    d = sorted_durations[x]\n",
    "    \n",
    "    deltas.append(d - last_marker)\n",
    "    last_marker = d\n",
    "    \n",
    "    matplotlib.pyplot.plot([x, x], [0, mp])\n",
    "    \n",
    "sorted_deltas = sorted(deltas)\n",
    "print()\n",
    "print(\"\\n\".join([str(int(round(x/sorted_deltas[0], 0))) for x in sorted_deltas]))\n",
    "print()\n",
    "\n",
    "midpoint = int(len(sorted_deltas) / 2)\n",
    "if len(sorted_deltas) % 2 == 1:\n",
    "    median_of_medians = sorted_deltas[midpoint]\n",
    "else:\n",
    "    median_of_medians = (sorted_deltas[midpoint - 1] + sorted_deltas[midpoint]) / 2\n",
    "\n",
    "mean_of_medians = numpy.mean(sorted_deltas)\n",
    "    \n",
    "xdata = numpy.array(range(ls))\n",
    "slope = median_of_medians / points_per_bin\n",
    "intercept = sorted_durations[int(ls/2)] - slope * xdata[int(ls/2)]\n",
    "ydata = slope * xdata + intercept\n",
    "matplotlib.pyplot.plot(ydata)\n",
    "\n",
    "p25 = sorted_durations[int(lsd * start / num_deltas)]\n",
    "p75 = sorted_durations[int(lsd * (finish - 1) / num_deltas)]\n",
    "IQR = (p75 - p25) / 1.349\n",
    "\n",
    "print()\n",
    "print(\"IQR delta\", IQR / num_deltas)\n",
    "print(\"mean of deltas\", mean_of_medians)\n",
    "print(\"median of deltas\", median_of_medians)\n",
    "print(\"sample std\", numpy.std(sorted_durations) / num_deltas)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Kernel Estimation Classes"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": 390,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "class Kernel(object):\n",
    "    def __init__(self, width, offsets=8):\n",
    "        self.width = width + 1\n",
    "        self.base_kernel = list(\n",
    "            norm.pdf(\n",
    "                numpy.linspace(\n",
    "                    norm.ppf(0.01),\n",
    "                    norm.ppf(0.99),\n",
    "                    width)))\n",
    "        \n",
    "        self.offsets = offsets\n",
    "        self.kernels = []\n",
    "        \n",
    "        start = 0.0\n",
    "        for x in range(offsets):\n",
    "            end = x/offsets\n",
    "            self.kernels.append((start, end, self._generate_shift(x/offsets)))\n",
    "        \n",
    "    def shift(self, offset):\n",
    "        low = 0\n",
    "        high = self.offsets - 1\n",
    "                                \n",
    "        while True:\n",
    "            midpoint = (low + high) // 2\n",
    "            start = self.kernels[midpoint][0]\n",
    "            end = self.kernels[midpoint][1]\n",
    "            if (offset >= start) and (offset < end):\n",
    "                return self.kernels[midpoint][2]\n",
    "                                \n",
    "            elif offset <= end:\n",
    "                if midpoint == 0:\n",
    "                    return self.kernels[0][2]\n",
    "                high = midpoint\n",
    "\n",
    "            elif offset > start:\n",
    "                if midpoint == self.offsets - 1:\n",
    "                    return self.kernels[0][2]\n",
    "                low = midpoint\n",
    "            else:\n",
    "                raise Exception(\"What?\")\n",
    "\n",
    "    def _generate_shift(self, fractional_offset):\n",
    "        leftovers = 0\n",
    "        shifted = []\n",
    "\n",
    "        for kslice in self.base_kernel:\n",
    "            shifted_part = fractional_offset * kslice\n",
    "            shifted.append(leftovers + fractional_offset * kslice)\n",
    "            leftovers = kslice - shifted_part\n",
    "\n",
    "        shifted.append(leftovers)\n",
    "        return shifted"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 515,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "class PrecisionRegion(object):\n",
    "    def __init__(self, index, precision, start, end, ratio_to_previous, total_ratio):\n",
    "        self.index = index\n",
    "        self.precision = precision\n",
    "        self.start = start\n",
    "        self.end = end\n",
    "        self.ratio_to_previous = int(ratio_to_previous)\n",
    "        self.total_ratio = int(total_ratio)\n",
    "\n",
    "class PrecisionIndex(object):\n",
    "    def __init__(self, base_precision, kernel_width):\n",
    "        self.base_precision = base_precision\n",
    "        self.kernel_width = kernel_width\n",
    "        self.hkw = self.kernel_width // 2\n",
    "        self.regions = []\n",
    "\n",
    "    def add_region(self, start, end, ratio_to_previous, total_ratio):\n",
    "        precision = self.base_precision * total_ratio\n",
    "        pr = PrecisionRegion(len(self.regions), precision, start, end, ratio_to_previous, total_ratio)\n",
    "        self.regions.append(pr)\n",
    "\n",
    "    def construct(self, points_per_tile, sorted_durations, deltas):\n",
    "        self.minimum_allowed = self.base_precision * self.hkw\n",
    "        \n",
    "        # first pass: split the deltas into chunks\n",
    "\n",
    "        chunks = [(0, 1)]\n",
    "        current_ratio = 1\n",
    "        primed = True\n",
    "        for i in range(len(deltas)):\n",
    "            cp = deltas[i]\n",
    "            new_ratio = max(1, cp // self.base_precision)\n",
    "            if primed:\n",
    "                if (new_ratio > 2*current_ratio) or (2*new_ratio < current_ratio):\n",
    "                    chunks.append([i, new_ratio])\n",
    "                    primed = False\n",
    "            else:\n",
    "                chunks[-1][-1] = new_ratio\n",
    "                primed = True\n",
    "            \n",
    "            current_ratio = new_ratio\n",
    "            \n",
    "        # second pass: calculate regions\n",
    "        \n",
    "        regions = []\n",
    "        ratio_to_previous = 1\n",
    "        \n",
    "        for c in range(len(chunks)):\n",
    "            i, ratio = chunks[c]\n",
    "            \n",
    "            if c > 0:\n",
    "                start = sorted_durations[(i+2) * points_per_tile]\n",
    "                regions[-1][1] = start\n",
    "            else:\n",
    "                start = self.minimum_allowed\n",
    "            \n",
    "            aligned_ratio = 2 ** int(math.log2(ratio))\n",
    "            ratio_to_previous = int(aligned_ratio / ratio_to_previous)\n",
    "            precision = self.base_precision * aligned_ratio\n",
    "            \n",
    "            regions.append([start, end, ratio_to_previous, aligned_ratio])\n",
    "            start = end\n",
    "            \n",
    "        regions[-1][1] = (sorted_durations[points_per_tile * (len(deltas) - 1)]\n",
    "                          + 1000 * self.base_precision * regions[-1][-1])\n",
    "            \n",
    "        # third pass: decrease the ratio of each region so that each region\n",
    "        #             never overlaps more than one on each side of it\n",
    "        \n",
    "        lr = len(regions)\n",
    "        for r in range(1, lr):\n",
    "            current = regions[r]\n",
    "            \n",
    "            self._reduce_ratio(current, regions[r-1][0], self._reduce_ratio_prev)\n",
    "            \n",
    "            if r < lr - 1:\n",
    "                self._reduce_ratio(current, regions[r+1][0], self._reduce_ratio_next)\n",
    "             \n",
    "        # fourth pass: merge any adjacent chunks with the same ratio\n",
    "        \n",
    "        current_ratio = None\n",
    "        for region in regions:\n",
    "            if region[-1] == current_ratio:\n",
    "                continue\n",
    "\n",
    "            self.add_region(*region)\n",
    "            current_ratio = region[-1]\n",
    "                    \n",
    "    def _reduce_ratio(self, region, threshold, calc_edge):\n",
    "        while True:\n",
    "            ce = calc_edge(region)\n",
    "            \n",
    "            if calc_edge(region) > threshold:\n",
    "                break\n",
    "\n",
    "            region[-1] /= 2\n",
    "            region[-2] /= 2\n",
    "\n",
    "            if region[-2] <= 1:\n",
    "                break\n",
    "                \n",
    "    def _reduce_ratio_prev(self, region):\n",
    "        return region[0] - self.base_precision * self.hkw * region[-1]\n",
    "    \n",
    "    def _reduce_ratio_next(self, region):\n",
    "        return region[1] + self.base_precision * self.hkw * region[-1]\n",
    "    \n",
    "    def lookup(self, duration):\n",
    "        if duration <= self.minimum_allowed:\n",
    "            return self.regions[0]\n",
    "        if duration <= self.maximum_allowed:\n",
    "            return self.regions[-1]\n",
    "        \n",
    "        last_node = len(self.regions) - 1\n",
    "        low = 0\n",
    "        high = last_node\n",
    "        \n",
    "        while True:\n",
    "            midpoint = (low + high) // 2\n",
    "            node = self.regions[midpoint]\n",
    "            \n",
    "            if (duration >= node.start) and (duration < node.end):\n",
    "                return node.precision\n",
    "            \n",
    "            if duration < node.end:\n",
    "                if midpoint == 0:\n",
    "                    return self.regions[0]\n",
    "                high = midpoint\n",
    "                \n",
    "            elif duration >= node.start:\n",
    "                if midpoint == last_node:\n",
    "                    return self.regions[last_node]\n",
    "                low = midpoint\n",
    "\n",
    "            else:\n",
    "                raise Exception(\"What?\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 517,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "15 8.333333333333334 8\n",
      "\n",
      "1.171875e-07\n",
      "\n",
      "[1.1718749999999988e-07, 1.171875000000001e-07, 1.171875000000001e-07, 6.7007415254237281e-06, 2.5423728813559366e-07]\n",
      "\n",
      "['8.203125e-07/8.16949152542e-06/1', '8.16949152542e-06/0.00023584375/2']\n"
     ]
    }
   ],
   "source": [
    "a = list(numpy.linspace(1e-6, 1.5e-6, 65))\n",
    "b = list(numpy.linspace(8e-6, 9e-6, 60))\n",
    "\n",
    "sd = a + b\n",
    "lsd = len(sd)\n",
    "\n",
    "points_per_tile = max(15, int(numpy.sqrt(lsd)))\n",
    "tile = lsd / points_per_tile\n",
    "num_tiles = int(tile)\n",
    "\n",
    "print(points_per_bin, tile, num_tiles)\n",
    "print()\n",
    "\n",
    "tile_deltas = []\n",
    "last = sd[points_per_tile]\n",
    "for i in range(2, num_tiles - 1):\n",
    "    d = sd[int(i * points_per_bin)]\n",
    "    tile_deltas.append(d - last)\n",
    "    last = d\n",
    "\n",
    "sorted_deltas = sorted(tile_deltas)\n",
    "base_precision = sorted_deltas[int(len(tile_deltas)/3)]\n",
    "\n",
    "pi = PrecisionIndex(base_precision, 15)\n",
    "pi.construct(points_per_tile, sd, tile_deltas)\n",
    "\n",
    "print(base_precision)\n",
    "print()\n",
    "print(tile_deltas)\n",
    "print()\n",
    "print([\"%s/%s/%s\" % (r.start, r.end, r.total_ratio) for r in pi.regions])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 528,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "class Histogram(object):\n",
    "    def __init__(self, kernel_width, precision_region):\n",
    "        self.kernel_width = kernel_width\n",
    "        self.hkw = kernel_width // 2\n",
    "        self.kernel_range = range(kernel_width)\n",
    "        \n",
    "        self.precision_region = precision_region\n",
    "        \n",
    "        p = precision_region.precision\n",
    "        start = precision_region.start\n",
    "        end = precision_region.end\n",
    "        \n",
    "        self.starting_index_offset = int(round(start / p, 0))\n",
    "        \n",
    "        self.current_maximum = (0, 0)\n",
    "        self.size = (end - start) / p\n",
    "        self.size = int(round(self.size, 0)) + self.kernel_width\n",
    "        self.trail_head = self.hkw\n",
    "        self.trail_tail = self.size - self.kernel_width - self.hkw\n",
    "        \n",
    "        self.histogram = [0] * self.size\n",
    "        self.lower_neighbor = None\n",
    "        self.upper_neighbor = None\n",
    "\n",
    "    def align(self, value):\n",
    "        p = self.precision_region.precision\n",
    "        return (value // p) * p\n",
    "    \n",
    "    def hslice(self, index):\n",
    "        return self.histogram[index : index+self.kernel_width].copy()\n",
    "    \n",
    "    # note: starting_index is the point on the left where the kernel starts;\n",
    "    #       that is to say, aligned_duration corresponds with the histogram\n",
    "    #       index of starting_index + hkw\n",
    "    def calc_starting_index(self, aligned_duration):\n",
    "        p = self.precision_region.precision\n",
    "        starting_index = int(round(aligned_duration / p, 0))\n",
    "        starting_index -= self.starting_index_offset\n",
    "        return starting_index\n",
    "        \n",
    "    def collate_values(self, aligned_duration):\n",
    "        starting_index = self.calc_starting_index(aligned_duration)\n",
    "        values = self.hslice(starting_index)\n",
    "        \n",
    "        if self.lower_neighbor and (starting_index < self.trail_head):\n",
    "            self._include_lower_overlap(starting_index, values)\n",
    "        elif self.upper_neighbor and (starting_index > self.trail_tail):\n",
    "            self._include_upper_overlap(starting_index, values)\n",
    "        \n",
    "        return values\n",
    "        \n",
    "    def _include_lower_overlap(self, starting_index, values):\n",
    "        oh = self.lower_neighbor.histogram\n",
    "        ohs = self.lower_neighbor.size - self.hkw - 1\n",
    "\n",
    "        ratio = self.precision_region.ratio_to_previous\n",
    "        factor = ratio / 2\n",
    "        overlap = self.trail_head - starting_index\n",
    "\n",
    "        for i in range(0, overlap):\n",
    "            values[overlap - i] += oh[ohs - i*ratio] * factor\n",
    "                \n",
    "    def _include_upper_overlap(self, starting_index, values):\n",
    "        oh = self.upper_neighbor.histogram\n",
    "\n",
    "        ratio = self.upper_neighbor.precision_region.ratio_to_previous\n",
    "        factor = 1 / ratio\n",
    "        overlap = starting_index - self.trail_tail\n",
    "        overlap_bins = overlap / ratio\n",
    "        whole_overlap_bins = int(math.ceil(overlap_bins))\n",
    "\n",
    "        for i in range(0, whole_overlap_bins):\n",
    "            if (i == whole_overlap_bins - 1) and (i - overlap_bins >= 0.1):\n",
    "                factor *= overlap_bins - i + 1\n",
    "            values[self.kernel_width - whole_overlap_bins + i] += oh[i] * factor\n",
    "\n",
    "    def add_kernel(self, aligned_duration, shifted_kernel):\n",
    "        starting_index = self.calc_starting_index(aligned_duration)\n",
    "        \n",
    "        for i in self.kernel_range:\n",
    "            self.histogram[starting_index + i] += shifted_kernel[i]\n",
    "            \n",
    "        return self.hslice(starting_index)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 529,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "ename": "IndexError",
     "evalue": "list index out of range",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mIndexError\u001b[0m                                Traceback (most recent call last)",
      "\u001b[0;32m<ipython-input-529-2e05c118c583>\u001b[0m in \u001b[0;36m<module>\u001b[0;34m()\u001b[0m\n\u001b[1;32m      1\u001b[0m \u001b[0mh\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mHistogram\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mpi\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mkernel_width\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mpi\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mregions\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;36m0\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m----> 2\u001b[0;31m \u001b[0mh\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0madd_kernel\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;36m5e-5\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mk\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mshift\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;36m0.5\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m",
      "\u001b[0;32m<ipython-input-528-1df0226ab5cc>\u001b[0m in \u001b[0;36madd_kernel\u001b[0;34m(self, aligned_duration, shifted_kernel)\u001b[0m\n\u001b[1;32m     79\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     80\u001b[0m         \u001b[0;32mfor\u001b[0m \u001b[0mi\u001b[0m \u001b[0;32min\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mkernel_range\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m---> 81\u001b[0;31m             \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mhistogram\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0mstarting_index\u001b[0m \u001b[0;34m+\u001b[0m \u001b[0mi\u001b[0m\u001b[0;34m]\u001b[0m \u001b[0;34m+=\u001b[0m \u001b[0mshifted_kernel\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0mi\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m     82\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     83\u001b[0m         \u001b[0;32mreturn\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mhslice\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mstarting_index\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;31mIndexError\u001b[0m: list index out of range"
     ]
    }
   ],
   "source": [
    "h = Histogram(pi.kernel_width, pi.regions[0])\n",
    "h.add_kernel(5e-5, k.shift(0.5))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "class HistogramIndex(object):\n",
    "    def __init__(self, kernel, precision_index):\n",
    "        self.kernel = kernel\n",
    "        self.precision_index = precision_index\n",
    "        self.histograms = []\n",
    "    \n",
    "        for pr in self.precision_index.regions:\n",
    "            self.histograms[pr.precision] = Histogram(self.kernel.width, pr)\n",
    "            \n",
    "        for i in range(len(self.histograms)):\n",
    "            if i > 0:\n",
    "                self.histograms[i].lower_neighbor = self.histograms[i-1]\n",
    "            if i < len(self.histograms) - 1:\n",
    "                self.histograms[i].upper_neighbor = self.histograms[i+1]\n",
    "                \n",
    "    def add_kernel(self, duration, calc_diff):\n",
    "        precision_region = self.precision_index.lookup(duration)\n",
    "        p = precision_region.precision\n",
    "        h = self.histograms[p]\n",
    "        \n",
    "        aligned_duration = h.align(duration)\n",
    "        shifted_kernel = self.kernel.shift((duration - aligned_duration) / p)\n",
    "        \n",
    "        old_values = h.collate_values(aligned_duration) if calc_diff else []\n",
    "        new_values = h.add_kernel(aligned_duration, shifted_kernel)\n",
    "        \n",
    "        return (old_values, new_values)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 361,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "class KernelEstimator(object):\n",
    "    def __init__(self):\n",
    "        self.histograms = {}\n",
    "        self.precisions = None\n",
    "        \n",
    "        self.base_precision = -1\n",
    "        self.minimum_start = -1\n",
    "        self.maximim_end = -1\n",
    "\n",
    "        self.seed_mode = True\n",
    "        self.seed_data = []\n",
    "        self.seed_size = 0\n",
    "        self.total_size = 0\n",
    "        self.last_diff_total = 0\n",
    "        \n",
    "    def finish_seed_mode(self):\n",
    "        self.seed_mode = False\n",
    "        self._calculate_precision()\n",
    "\n",
    "        self.seed_size = len(self.seed_data)\n",
    "        #for duration in self.seed_data:\n",
    "        #    self.add_data_point(duration, calc_diff=False)\n",
    "            \n",
    "    def plot_histogram(self):\n",
    "        for precision in self.histograms:\n",
    "            xdata = sorted([float(x) for x in self.histograms[precision].keys()])\n",
    "            ydata = [self.histograms[precision]['%5.3E' % x] for x in xdata]\n",
    "\n",
    "            matplotlib.pyplot.title('precision = %s' % precision)\n",
    "            matplotlib.pyplot.xlim(min(xdata)/1.2, max(xdata)*1.05)\n",
    "            matplotlib.pyplot.scatter(xdata, ydata)\n",
    "        \n",
    "    #############################################################\n",
    "        \n",
    "    def _gen_kernels(self, sizes):\n",
    "        for size in sizes:\n",
    "            n = self.kernel_width * size\n",
    "            if n%2 == 0:\n",
    "                n += 1\n",
    "                \n",
    "            self.kernels[size] = self._gen_kernel(n)\n",
    "            self.kernels[size] /= sum(self.kernels[size])\n",
    "            \n",
    "    def _gen_kernel(self, size):\n",
    "        return numpy.ones(size)\n",
    "        \n",
    "    def _align(self, value, precision):\n",
    "        dividend = value // precision\n",
    "        return dividend * precision\n",
    "    \n",
    "    #############################################################\n",
    "\n",
    "    # TODO: is there a way we can estimate the percentiles on-the-fly,\n",
    "    #       without storing every seed data point?\n",
    "    #       maybe https://en.wikipedia.org/wiki/Count%E2%80%93min_sketch ?\n",
    "    def _calculate_precision(self):\n",
    "        # TODO: this is bunk, as kernel_width should probably *SHRINK* with increasing number of points\n",
    "        #\n",
    "        # self.kernel_width = max(11, int(2 * points_per_bin ** 0.5))\n",
    "        # if self.kernel_width % 2 == 0:\n",
    "        #     self.kernel_width += 1\n",
    "        \n",
    "        self.kernel_width = 15\n",
    "        \n",
    "        sorted_durations = sorted(self.seed_data)\n",
    "        last = sorted_durations[points_per_bin]\n",
    "        \n",
    "        lsd = len(sorted_durations)\n",
    "        points_per_bin = max(self.kernel_width, int(numpy.sqrt(lsd)))\n",
    "        tile = lsd / self.points_per_bin\n",
    "        self.num_tiles = int(tile)\n",
    "        \n",
    "        tile_deltas = []\n",
    "        for i in range(2, self.num_tiles - 2):\n",
    "            d = sorted_durations[int(i * tile)]\n",
    "            tile_deltas.append(d - last)\n",
    "            last = d\n",
    "\n",
    "        sorted_deltas = sorted(tile_deltas)\n",
    "        self.base_precision = sorted_deltas[int(len(tile_deltas)/3)]\n",
    "        \n",
    "        self.precision_index = PrecisionIndex(self.base_precision, self.kernel_width)\n",
    "        self.precision_index.construct_regions(sorted_durations, tile_deltas)\n",
    "        self.histograms = HistogramIndex(self.kernel_width, self.precision_index)\n",
    "            \n",
    "    #############################################################\n",
    "\n",
    "    def add_data_point(self, duration, calc_diff=True):\n",
    "        if self.seed_mode:\n",
    "            self.seed_data.append(duration)\n",
    "            return\n",
    "\n",
    "        self.total_size += 1\n",
    "        (old_values, new_values) = self.histograms.add(duration, calc_diff)\n",
    "        \n",
    "        if not calc_diff:\n",
    "            return\n",
    "\n",
    "        abs_diff_total = 0\n",
    "        prev_ratio = 0.0\n",
    "        scaling_factor = 1 # self.total_size / self.seed_size\n",
    "        min_delta = self.min_delta\n",
    "\n",
    "        for old, new in zip(old_values, new_values):\n",
    "            new_ratio = -1 + (new + min_delta) / (old + min_delta)\n",
    "            abs_diff_total += abs(new_ratio - prev_ratio)\n",
    "            prev_ratio = new_ratio\n",
    "                \n",
    "        # TODO: need to compare this to histogram max\n",
    "        self.last_diff_total = abs_diff_total\n",
    "        return abs_diff_total > scaling_factor / 2"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Quick Test for the Estimator"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 363,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "range(-8, 4)\n",
      "\n",
      "3.17540242123e-09\n"
     ]
    }
   ],
   "source": [
    "%matplotlib inline\n",
    "\n",
    "seed_points = 1000\n",
    "data_points = 1000\n",
    "g = GaussianKernelEstimator()\n",
    "\n",
    "seeds = numpy.array(\n",
    "        list(gen_data(seed_points, 1e-7))\n",
    "      + list(gen_data(seed_points, 2e-7)))\n",
    "\n",
    "for duration in seeds:\n",
    "    g.add_data_point(duration)\n",
    "\n",
    "g.finish_seed_mode()\n",
    "\n",
    "#ax = matplotlib.pyplot.subplot(2, 1, 1)\n",
    "# ax.set_xscale(\"log\")\n",
    "#g.plot_histogram()\n",
    "#before_buckets = len(g.histogram.values())\n",
    "\n",
    "\"\"\"\n",
    "data = gen_data(data_points, 1e-7)\n",
    "for duration in data:\n",
    "    is_anomaly = g.add_data_point(duration)\n",
    "    \n",
    "    if is_anomaly:\n",
    "        print(g.total_size, duration, g.last_diff_total, g.total_size / g.seed_size)\n",
    "\n",
    "ax = matplotlib.pyplot.subplot(2, 1, 2)\n",
    "\"\"\"\n",
    "# ax.set_xscale(\"log\")\n",
    "#g.plot_histogram()\n",
    "#print()\n",
    "#print(before_buckets)\n",
    "#print(len(g.histogram.values()))\n",
    "#print()\n",
    "#print(g.kernel_width)\n",
    "#print(g.precision)\n",
    "print(g.base_precision)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 364,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[(2.2227816948610395e-08,\n",
       "  1.3971770653412249e-07,\n",
       "  '3.175E-09',\n",
       "  3.1754024212300565e-09,\n",
       "  1,\n",
       "  1),\n",
       " (1.3971770653412249e-07,\n",
       "  6.5349781828914564e-06,\n",
       "  '6.351E-09',\n",
       "  6.350804842460113e-09,\n",
       "  2,\n",
       "  2)]"
      ]
     },
     "execution_count": 364,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "g.precisions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 366,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "{'1.016E-07': {'3.175E-09': [('1.016E-07', 1.0161287747936181e-07, 0.5)],\n",
       "  '6.351E-09': [('1.016E-07', 1.0161287747936181e-07, 2),\n",
       "   ('1.048E-07', 1.0478827990059186e-07, 2)]},\n",
       " '1.048E-07': {'3.175E-09': [('1.016E-07', 1.0161287747936181e-07, 0.5)]},\n",
       " '1.080E-07': {'3.175E-09': [('1.080E-07', 1.0796368232182192e-07, 0.5)],\n",
       "  '6.351E-09': [('1.080E-07', 1.0796368232182192e-07, 2),\n",
       "   ('1.111E-07', 1.1113908474305198e-07, 2)]},\n",
       " '1.111E-07': {'3.175E-09': [('1.080E-07', 1.0796368232182192e-07, 0.5)]},\n",
       " '1.143E-07': {'3.175E-09': [('1.143E-07', 1.1431448716428203e-07, 0.5)],\n",
       "  '6.351E-09': [('1.143E-07', 1.1431448716428203e-07, 2),\n",
       "   ('1.175E-07', 1.1748988958551209e-07, 2)]},\n",
       " '1.175E-07': {'3.175E-09': [('1.143E-07', 1.1431448716428203e-07, 0.5)]},\n",
       " '1.207E-07': {'3.175E-09': [('1.207E-07', 1.2066529200674215e-07, 0.5)],\n",
       "  '6.351E-09': [('1.207E-07', 1.2066529200674215e-07, 2),\n",
       "   ('1.238E-07', 1.2384069442797219e-07, 2)]},\n",
       " '1.238E-07': {'3.175E-09': [('1.207E-07', 1.2066529200674215e-07, 0.5)]},\n",
       " '1.270E-07': {'3.175E-09': [('1.270E-07', 1.2701609684920226e-07, 0.5)],\n",
       "  '6.351E-09': [('1.270E-07', 1.2701609684920226e-07, 2),\n",
       "   ('1.302E-07', 1.3019149927043233e-07, 2)]},\n",
       " '1.302E-07': {'3.175E-09': [('1.270E-07', 1.2701609684920226e-07, 0.5)]},\n",
       " '1.334E-07': {'3.175E-09': [('1.334E-07', 1.3336690169166237e-07, 0.5)],\n",
       "  '6.351E-09': [('1.334E-07', 1.3336690169166237e-07, 2),\n",
       "   ('1.365E-07', 1.3654230411289242e-07, 2)]},\n",
       " '1.365E-07': {'3.175E-09': [('1.334E-07', 1.3336690169166237e-07, 0.5)]},\n",
       " '1.397E-07': {'3.175E-09': [('1.397E-07', 1.3971770653412249e-07, 0.5)],\n",
       "  '6.351E-09': [('1.397E-07', 1.3971770653412249e-07, 2),\n",
       "   ('1.429E-07', 1.4289310895535256e-07, 2)]},\n",
       " '1.429E-07': {'3.175E-09': [('1.397E-07', 1.3971770653412249e-07, 0.5)]},\n",
       " '1.461E-07': {'3.175E-09': [('1.461E-07', 1.460685113765826e-07, 0.5)],\n",
       "  '6.351E-09': [('1.461E-07', 1.460685113765826e-07, 2),\n",
       "   ('1.492E-07', 1.4924391379781264e-07, 2)]},\n",
       " '1.492E-07': {'3.175E-09': [('1.461E-07', 1.460685113765826e-07, 0.5)]},\n",
       " '1.524E-07': {'3.175E-09': [('1.524E-07', 1.5241931621904271e-07, 0.5)],\n",
       "  '6.351E-09': [('1.524E-07', 1.5241931621904271e-07, 2),\n",
       "   ('1.556E-07', 1.5559471864027278e-07, 2)]},\n",
       " '1.556E-07': {'3.175E-09': [('1.524E-07', 1.5241931621904271e-07, 0.5)]},\n",
       " '1.588E-07': {'3.175E-09': [('1.588E-07', 1.5877012106150282e-07, 0.5)],\n",
       "  '6.351E-09': [('1.588E-07', 1.5877012106150282e-07, 2),\n",
       "   ('1.619E-07', 1.6194552348273287e-07, 2)]},\n",
       " '1.619E-07': {'3.175E-09': [('1.588E-07', 1.5877012106150282e-07, 0.5)]},\n",
       " '8.891E-08': {'3.175E-09': [('8.891E-08', 8.8911267794441582e-08, 0.5)],\n",
       "  '6.351E-09': [('8.891E-08', 8.8911267794441582e-08, 2),\n",
       "   ('9.209E-08', 9.2086670215671638e-08, 2)]},\n",
       " '9.209E-08': {'3.175E-09': [('8.891E-08', 8.8911267794441582e-08, 0.5)]},\n",
       " '9.526E-08': {'3.175E-09': [('9.526E-08', 9.5262072636901695e-08, 0.5)],\n",
       "  '6.351E-09': [('9.526E-08', 9.5262072636901695e-08, 2),\n",
       "   ('9.844E-08', 9.8437475058131751e-08, 2)]},\n",
       " '9.844E-08': {'3.175E-09': [('9.526E-08', 9.5262072636901695e-08, 0.5)]}}"
      ]
     },
     "execution_count": 366,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "g.overlaps"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 51,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x1a21d66cf8>]"
      ]
     },
     "execution_count": 51,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAMcAAAD8CAYAAADDuLCoAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvNQv5yAAADq1JREFUeJzt3X+MXNV5xvHv4x9ExdBg88Mxxg5p40JJWgisDJHbCJLgGAvFtEKt3SpxWiSnEUhBSqXQRgoRUSWqKkRNHYEIWJiKOqglFFchgEWiENRA2HXNrxpig6BsbHkBU4NF1LDw9o85Ww3j++6M5+7uzGyfjzSaueeee8+5oMdz587sexURmNmR5vR6Amb9yuEwSzgcZgmHwyzhcJglHA6zRNtwSFom6UeSdkt6WtIXS/siSTsk7SnPC5PtN5Y+eyRtnOoDMJsuavc9h6QlwJKI2CnpeGAEuAz4HHAwIq6XdA2wMCK+3LLtImAYGAKibHteRLw25UdiNsXavnNExP6I2FlevwHsBpYC64CtpdtWGoFp9SlgR0QcLIHYAayZiombTbd5R9NZ0unAR4BHgcURsR8aAZJ0SsUmS4GXmpZHS1vVvjcBmwAWLFhw3plnnnk0UzPr2MjIyCsRcXK7fh2HQ9JxwF3A1RHxuqSONqtoqzyPi4ibgZsBhoaGYnh4uNOpmR0VSS920q+jq1WS5tMIxh0R8b3SfKB8Hpn4XDJWsekosKxp+TRgXydjmvVaJ1erBNwK7I6IG5pWbQcmrj5tBO6p2Px+YLWkheVq1urSZtb3OnnnWAV8Bvi4pF3lsRa4HrhY0h7g4rKMpCFJtwBExEHg68Bj5XFdaTPre20v5faCP3PYdJI0EhFD7fr5G3KzhMNhlnA4zBIOh1nC4TBLOBxmCYfDLOFwmCUcDrOEw2GWcDjMEg6HWcLhMEs4HGYJh8Ms4XCYJRwOs4TDYZZoW5pH0hbgUmAsIj5c2u4EzihdTgD+OyLOqdj2BeAN4G1gvJM/TTTrF53UrboN2AzcPtEQEX888VrSN4BDk2x/UUS80u0EzXqlbTgi4qFS6fAIpWzPHwEfn9ppmfVe3c8cvw8ciIg9yfoAHpA0Usp9mg2Mo6qVW2EDsG2S9asiYl+po7tD0jMR8VBVx+ZaucuXL685LbP6un7nkDQP+EPgzqxPROwrz2PA3cDKSfreHBFDETF08slta/yaTbs6p1WfBJ6JiNGqlZIWlPt5IGkBjVKgT9UYz2xGdVIrdxvwU+AMSaOSriir1tNySiXpVEn3lsXFwMOSHgd+Bnw/Iu6buqmbTa9OrlZtSNo/V9G2D1hbXj8PnF1zfmY942/IzRIOh1nC4TBLOBxmCYfDLOFwmCUcDrOEw2GWcDjMEg6HWcLhMEs4HGYJh8Ms4XCYJRwOs4TDYZZwOMwSDodZwuEwS3RSYGGLpDFJTzW1fU3SLyTtKo+1ybZrJD0raa+ka6Zy4mbTrZN3jtuANRXt34yIc8rj3taVkuYC3wYuAc4CNkg6q85kzWZS23CUCoUHu9j3SmBvRDwfEb8Cvgus62I/Zj1R5zPHVZKeKKddCyvWLwVealoeLW2VJG2SNCxp+OWXX64xLbOp0W04bgR+EzgH2A98o6KPKtoi26HLgVq/6SocEXEgIt6OiHeA71BdA3cUWNa0fBqwr5vxzHqhq3BIWtK0+AdU18B9DFgh6QOSjqFRPnR7N+OZ9UIntz3bBlwInCRpFLgWuFDSOTROk14APl/6ngrcEhFrI2Jc0lXA/cBcYEtEPD0tR2E2DRSRfgzomaGhoRgeHu71NGyWkjTSyf0p/Q25WcLhMEs4HGYJh8Ms4XCYJRwOs4TDYZZwOMwSDodZwuEwSzgcZgmHwyzhcJglHA6zhMNhlnA4zBIOh1nC4TBLdFsO9O8kPVPqVt0t6YRk2xckPVlKhvrvXm2gdFsOdAfw4Yj4XeDnwF9Nsv1FpWRo27/ZNesnXZUDjYgHImK8LD5CoyaV2awyFZ85/hz4QbIugAckjUjaNNlOXA7U+k2tcEj6CjAO3JF0WRUR59KotH6lpI9l+3I5UOs3XYdD0kbgUuBPIyl+FRH7yvMYcDfVZUPN+lK35UDXAF8GPh0RbyZ9Fkg6fuI1sJrqsqFmfamTS7nbgJ8CZ0galXQFsBk4HthRLtPeVPqeKmniRjaLgYclPQ78DPh+RNw3LUdhNg3a1sqNiA0VzbcmffcBa8vr54Gza83OrIf8DblZwuEwSzgcZgmHwyzhcJglHA6zhMNhlnA4zBIOh1nC4TBLOBxmCYfDLOFwmCUcDrOEw2GWcDjMEg6HWcLhMEt0FI6kJOgiSTsk7SnPC5NtN5Y+e0rFErOB0Ok7x20cWRL0GuDBiFgBPFiW30XSIuBa4HwaZXmuzUJk1m86CkdVSVBgHbC1vN4KXFax6aeAHRFxMCJeo1FjtzVkZn2pzmeOxRGxH6A8n1LRZynwUtPyaGk7gsuBWr+Z7g/kqmjLqiO6HKj1lTrhOCBpCUB5HqvoMwosa1o+DdhXY0yzGVMnHNuBiatPG4F7KvrcD6yWtLB8EF9d2sz6XqeXcqtKgl4PXCxpD3BxWUbSkKRbACLiIPB14LHyuK60mfU9JQXSe2poaCiGh32XNJsekkY6udOYvyE3SzgcZgmHwyzhcJglHA6zhMNhlnA4zBIOh1nC4TBLOBxmCYfDLOFwmCUcDrOEw2GWcDjMEg6HWcLhMEs4HGaJrsMh6QxJu5oer0u6uqXPhZIONfX5av0pm82Med1uGBHPAucASJoL/AK4u6LrTyLi0m7HMeuVqTqt+gTwXES8OEX7M+u5qQrHemBbsu6jkh6X9ANJH8p24HKg1m9qh0PSMcCngX+uWL0TeH9EnA38A/Cv2X5cDtT6zVS8c1wC7IyIA60rIuL1iDhcXt8LzJd00hSMaTbtpiIcG0hOqSS9T5LK65VlvFenYEyzadf11SoAScfSKAX6+aa2vwCIiJuAy4EvSBoHfgmsj34ssWhWoVY4IuJN4MSWtpuaXm8GNtcZw6xX/A25WcLhMEs4HGYJh8Ms4XCYJRwOs4TDYZZwOMwSDodZwuEwSzgcZgmHwyzhcJglHA6zhMNhlnA4zBIOh1nC4TBLTEVpnhckPVnKfQ5XrJekb0naK+kJSefWHdNsJtT6G/ImF0XEK8m6S4AV5XE+cGN5NutrM3FatQ64PRoeAU6QtGQGxjWrZSrCEcADkkYkbapYvxR4qWl5tLS9i8uBWr+ZinCsiohzaZw+XSnpYy3rVbHNEbWrXA7U+k3tcETEvvI8RuMWBCtbuowCy5qWTwP21R3XbLrVCoekBZKOn3gNrAaeaum2HfhsuWp1AXAoIvbXGddsJtS9WrUYuLuUw50H/FNE3NdSEvReYC2wF3gT+LOaY5rNiLrlQJ8Hzq5oby4JGsCVdcYx6wV/Q26WcDjMEg6HWcLhMEs4HGYJh8Ms4XCYJRwOs4TDYZZwOMwSDodZwuEwSzgcZgmHwyzhcJglHA6zhMNhlnA4zBJdh0PSMkk/krRb0tOSvljR50JJh0qp0F2SvtrteMMvHOSTN/yYvWOHu92F/T/xJ995hM0/3FN7P3X+hnwc+FJE7CwVSEYk7YiI/2zp95OIuLTGOACceNx72Dt2mIf3vAwEh//nbd56+x3eGn+Ht+PdZbBaFo8okhWt/TucQ2sBrlJYok2fqv1o0j5Vhb5aG9vto3ouFfNtM3bVflt7dTN21W6PnEv7cZrNkTj2mLn8+3OvsuqDJ03euQNdh6OU19lfXr8haTeNSoat4ZgSp594LHPniK/927Ts3maZC35jUe19TEkhaUmnAx8BHq1Y/VFJj9Mo5PaXEfF0l2PwpdW/xZOjh1j9ocW899fmM3/uHObPncO8Oe3/RWz3r12bf5Qq3n3a96rq024/re9qHW1T9d7X5t2zk/10Nv8O3oU7mO+R/x0mH7dq7PF3gv/4r9fYO3aY31l6QtVMjoqq/mcc1Q6k44AfA38TEd9rWffrwDsRcVjSWuDvI2JFsp9NwCaA5cuXn/fiiy/WmpdZRtJIRAy161e34uF84C7gjtZgAETE6xFxuLy+F5gvqfJk0LVyrd/UuVol4FZgd0TckPR5X+mHpJVlvFe7HdNsJtX5zLEK+AzwpKRdpe2vgeXwf1UPLwe+IGkc+CWwPuqex5nNkDpXqx6mzefYiNgMbO52DLNe8jfkZgmHwyzhcJglHA6zhMNhlnA4zBIOh1nC4TBLOBxmCYfDLOFwmCUcDrOEw2GWcDjMEg6HWcLhMEs4HGYJh8Ms4XCYJeqW5lkj6VlJeyVdU7H+PZLuLOsfLcXfzAZCndI8c4FvA5cAZwEbJJ3V0u0K4LWI+CDwTeBvux3PbKbVeedYCeyNiOcj4lfAd4F1LX3WAVvL638BPqHWisJmfapO3aqlwEtNy6PA+VmfiBiXdAg4EXildWfN5UCBw5KerRjzpKptZ4nZemz9eFzv76RTnXBUvQO0FmzrpE+jMeJm4OZJB5SGO6lxOohm67EN8nHVOa0aBZY1LZ9Go5J6ZR9J84D3AgdrjGk2Y+qE4zFghaQPSDoGWA9sb+mzHdhYXl8O/NDlQG1Q1CkHOi7pKuB+YC6wJSKelnQdMBwR22kUmv5HSXtpvGOsrznfSU+7BtxsPbaBPa7a9+cwm638DblZwuEwSwxMONr9VKXfSdoiaUzSU01tiyTtkLSnPC8s7ZL0rXKsT0g6t3czn1x2y+3ZcGwDEY4Of6rS724D1rS0XQM8WO6T+GBZhsZxriiPTcCNMzTHbkzccvu3gQuAK8v/m4E/toEIB539VKWvRcRDHPkdT/PPa7YClzW13x4NjwAnSFoyMzM9OhGxPyJ2ltdvABO33B74YxuUcFT9VGVpj+YylRaX+7lP3Nf9lNI+kMfbcsvtgT+2QQlHxz9DmSUG7njLLbfvAq6OiNcn61rR1pfHNijh6OSnKoPowMQpRXkeK+0DdbzJLbcH/tgGJRyd/FRlEDX/vGYjcE9T+2fLlZ0LgEMTpyj9ZpJbbg/8sRERA/EA1gI/B54DvtLr+XQx/23AfuAtGv96XkHj5/sPAnvK86LSVzSuzj0HPAkM9Xr+kxzX79E4LXoC2FUea2fDsfnnI2aJQTmtMptxDodZwuEwSzgcZgmHwyzhcJglHA6zxP8Cu4qriXcj5JoAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x1a20d9da20>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "from scipy.stats import truncnorm\n",
    "\n",
    "factor = 256\n",
    "nx = list(range(factor-1))\n",
    "ny = []\n",
    "\n",
    "for i in nx:\n",
    "    u1, v1 = truncnorm.stats(\n",
    "                 norm.ppf(max(0.5/factor, i/factor)),\n",
    "                 norm.ppf((i+1)/factor),\n",
    "                 moments='mv')\n",
    "    \n",
    "    u2, v2 = truncnorm.stats(\n",
    "                 norm.ppf((i+1)/factor),\n",
    "                 norm.ppf(min((factor-0.5)/factor, (i+2)/factor)),\n",
    "                 moments='mv')\n",
    "    \n",
    "    d = abs(u1 - u2) / ((16 * v1 * v2) ** (1/4))\n",
    "    dd = numpy.sqrt(d*d - 1)\n",
    "    t = 2 * numpy.log(d - dd) + 2 * d * dd\n",
    "    ny.append(t)\n",
    "\n",
    "ax = matplotlib.pyplot.subplot(1, 2, 1)\n",
    "ax.set_ylim(0, 20)\n",
    "matplotlib.pyplot.plot(nx, ny)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 75,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/usr/local/anaconda3/lib/python3.6/site-packages/ipykernel_launcher.py:32: RuntimeWarning: divide by zero encountered in log\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "(array([  73.,  533.,  389.,   57.,  114.,  273.,  319.,  210.,   67.,   13.]),\n",
       " array([ 0.00038404,  0.00071873,  0.00105342,  0.0013881 ,  0.00172279,\n",
       "         0.00205747,  0.00239216,  0.00272685,  0.00306153,  0.00339622,\n",
       "         0.00373091]),\n",
       " <a list of 10 Patch objects>)"
      ]
     },
     "execution_count": 75,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXwAAAD8CAYAAAB0IB+mAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvNQv5yAAAIABJREFUeJzt3Xl4XHd97/H3V5Il2VpsrbZsy5K3xLETx04UZ4VAFsji4FACN1BK0pt7XSjthcvtvYTmobS9D72BttALtwUCSQmFkkDC4jgBEpykwWRxZMdrbMe2vCmStXjT4mj/3j/m2JFtKR5JM5rlfF7PM8+cOefMOd8Z/fSZM7+zjLk7IiKS/jISXYCIiIwPBb6ISEgo8EVEQkKBLyISEgp8EZGQUOCLiITEqAPfzCrN7Dkz225m28zsM8H4YjN7xsx2BfdFsStXJP7UtiVd2WiPwzezCqDC3TeYWQGwHrgduBs44u73m9m9QJG7fz5WBYvEm9q2pKtRb+G7e6O7bwiG24HtwAxgBfBwMNvDRP5RRFKG2rakq1Fv4Z+2ELNq4AXgQuCAu08ZNO2ou5/11dfMVgIrAfLy8i5dsGDBmOsQGcr69etb3b1sNM9V25ZkNtK2nTXWFZpZPvA48Fl3bzOzqJ7n7g8ADwDU1NR4bW3tWEsRGZKZ7R/l89S2JamNtG2P6SgdM5tA5B/iR+7+s2B0U9AHerIvtHks6xBJBLVtSUdjOUrHgAeB7e7+tUGTVgF3BcN3Ab8cfXki409tW9LVWLp0rgb+CNhiZhuDcX8J3A/8xMzuAQ4AHx5biSLjTm1b0tKoA9/d1wLDdWpeP9rliiSa2rakK51pKyISEgp8EZGQUOCLiISEAl9EJCQU+CIiIaHAFxEJCQW+iEhIKPBFREJCgS8iEhIKfBGRkFDgi4iEhAJfRCQkFPgiIiGhwBcRCQkFvohISCjwRURCYiw/cfiQmTWb2dZB4/7azN40s43B7ZbYlCkyftS2JV2NZQv/+8BNQ4z/ursvCW5PjWH5IonyfdS2JQ2NOvDd/QXgSAxrEUkKatuSruLRh/9nZrY5+FpcNNxMZrbSzGrNrLalpSUOZYjEnNq2pLRYB/63gLnAEqAR+MfhZnT3B9y9xt1rysrKYlyGSMypbUvKi2ngu3uTu/e7+wDwXWBZLJcvkihq25IOYhr4ZlYx6OEHga3DzSuSStS2JR1kjfaJZvZj4D1AqZnVA18C3mNmSwAH9gF/EoMaRcaV2rakq1EHvrt/dIjRD46hFpGkoLYt6Upn2oqIhIQCX0QkJBT4IiIhocAXEQkJBb6ISEgo8EVEQkKBLyISEgp8EZGQUOCLiISEAl9EJCQU+CIiIaHAFxEJCQW+iEhIKPBFREJCgS8iEhIKfBGRkBhT4JvZQ2bWbGZbB40rNrNnzGxXcF809jJFxo/ataSrsW7hfx+46Yxx9wJr3H0+sCZ4LJJKvo/ataShMQW+u78AHDlj9Arg4WD4YeD2saxDZLypXUu6ikcf/lR3bwQI7suHmsnMVppZrZnVtrS0xKEMkZiKql2D2rYkr4TttHX3B9y9xt1rysrKElWGSMypbUuyikfgN5lZBUBw3xyHdYiMN7VrSXnxCPxVwF3B8F3AL+OwDpHxpnYtKW+sh2X+GHgJON/M6s3sHuB+4EYz2wXcGDwWSRlq15KussbyZHf/6DCTrh/LckUSSe1a0pXOtBURCQkFvohISCjwRURCQoGfYr65Zhe3fXMt7p7oUkQkxSjwU8jAgPOjVw6w5c3jvHbwWKLLEZEUo8BPIesPHOVQWxcAqzc1JrgaEUk1CvwUsnpTAzlZGVw9r4SntjQyMKBuHRGJngI/RfQPOE9uOcR1C8r5SE0lh9q6qN1/NNFliUgKUeCniFfqDtPa0c1tF0/nhgumkjshg9WbGxJdloikEAV+inhicyOTsjN57/nl5OVkcd2Ccp7acoh+deuISJQU+Cmgt3+AX29t5IYLpjIxOxOAWy+aTmtHN6/UHU5wdSKSKhT4KeDFPYc5eqKX5YsrTo27bkE5k7IzWb1FR+uISHQU+CngiU0NFORkce35b/+YxsTsTK6/YCq/3nqIvv6BBFYnIqlCgZ/kuvv6+c22Q7xv0TRysjJPm7Z8cQVHOnt4cY+6dUTk3BT4Se53b7TS3tXH8osrzpp27XllFORk6WgdEYlK3ALfzPaZ2RYz22hmtfFaT7pbvbmBKZMmcM280rOm5U7I5MaFkW6dnj5164wHtWtJZfHewn+vuy9x95o4ryctdfX288zrTdy0aBoTMof+U926uIK2rj7W7m4Z5+pCTe1aUpK6dJLY8zub6ezp57aLpw87z7vml1GYm8XqzTpaR0TeWTwD34GnzWy9ma08c6KZrTSzWjOrbWnR1ulQntjUSGl+NpfPLh52nuysDN6/aBrPbGuiq7d/HKsLrXds16C2LckrnoF/tbtfAtwMfNrM3j14ors/4O417l5TVlY29BJCrLO7jzU7mrj5wgqyhunOOWn5xdNp7+7jhTcULuPgHds1nLttuzu1+47Q0d03DuVKKttSf5zDHd0xW17cAt/dG4L7ZuDnwLJ4rSsdrdnRTFfvwGknWw3nqrklFE2aoG6dcRCLdr2toY07vv0Sv3jtzViXJ2nE3fnso6/xyR+uj9ky4xL4ZpZnZgUnh4H3AVvjsa50tXpTA1MLc7isevjunJMmZGZw04XT+O32Jt7qUbdOvMSqXS+aXsjCikJ++PJ+/XKZDOulusPsaenkIzWVMVtmvLbwpwJrzWwTsA540t1/Had1pZ32rl6ef6OFWy+aTkaGRfWc5Yunc6Knn+d2Nse5ulCLSbs2Mz5+RRU7DrWz4YAucS1D+9HLB5g8ccI7HrQxUnEJfHevc/eLg9sid/9yPNaTrp55vYmevoEhT7YazuWziynNz9ZJWHEUy3a9Ysl08nOy+Nff74thhZIuGo69xW+2HeKOS2eSOyHz3E+Ikg7LTEJPbGpgxpSJLK2cEvVzsjIzuPnCCp7d0UyndgYmvbycLO66qorVmxt5boe+lcnb3J3PP76ZrEzj7quqY7psBX6SOXaih9/tamX54grMouvOOWn54gq6egdYk4AAcXea27vYePAYx070jPv6U9F/u34+C6YV8D8f20xT8FvFIg+u3cvvdrVy3y0XUFk8KabLzorp0pLYwICzsf4YVcWTKMnPSXQ5w/rNtkP0DTjLF4+83+6y6mLKC3JYvamBD8Sw3++krt5+6o+e4MCRExw4fIIDR97iwJHOyOMjJ+jqjVzewQwWz5jM1fNKuWZ+KZdWFZ114TeBnKxM/unOJXzoX17kY999mUdWXklZQfK2TYm/x9fX8+WntnPjwql8/IqqmC8/7QO/r3+AJzY38O3n69jZ1E5hbhb33nwBd15WGfUO0fG0enMj1SWTuHBG4Yifm5Fh3HJRBf++7gDtXb0U5E4Ycz3uzm+3N/OPT+9kx6H206ZNnJBJVckkqkryeNf8MmYVT2JqYS47DrWxdlcr33mhjn95fg+5EzJYNruEd80r5ep5pSyYVpCU730iLJhWyL/+8TLuemgdH/rWi3z3EzWcP60g0WXJOHN3/uX5PfzD0zu5em4p3/zo0hF/w4+GJcNhYTU1NV5bG9vrUHX19vPT2oN854U66o++xXlT8/nEldWs2tTAur1HuLSqiC9/8EIWTBt5sAIcP9HLC7taKCvIYW5ZPqX52WP+A7V2dHP5363hU9fO5S/ef/6olrF+/xE+9K2X+NpHLuYPLpk5pnq2vnmcLz+5nZfqDjOnLI8VF8+gqmQSlcWTmFU86Zyvub2rl1fqjrB2dytrd7eyu7kDgNL8bK6aW8pl1UVcUlXE+VMLznly2ViY2fpEXfcm2ra94cBR/uTf1tPR1ccXblnAxy+v0odiSBw63sVf/nwLz+5o5raLp/PVDy0+9ct25zLStp12W/htXb388OX9PLR2H60d3SydNYUv3baI6xeUk5Fh/OHls3hsfT1/99R2ln9jLfe8azafuX4+k7Kjeyu2N7bxg5f28fPX3jzVhQFQkJvF3LJ85pblM6csj7ll+cwrz2NWcR7ZWdGF2a+2Rn6jdiRH55xpaWUR0yfnsnpz46gD/9DxLv7h6Z08vqGeKRMn8LcrFvHRZbOGvYDbcApyJ3DDwqncsHDqqeWu3d3K74Pbqk2RI4rysjO5uHIKl1ZFPgAuqSxi8qSxfztJJZfMKmL1n1/DX/x0E3/1y238/LU3+eLyhVwyqyjRpUmcdPX28/CL+/jms7vpGxjgS7ct5O6rquOyZX9S2mzht3Z089DavfzbS/tp7+7jXfNL+dP3zOOKOcVDvoFHOnu4/1fb+UltPTOmTOR/376I6xZMHXLZvf0DPL2tiYdf3Me6fUfInZDB7UtmcMelMznR009dSwd7WjrZ09JBXUsnhwbtgMvMMKqKJ3FBRSEXVBSwcHohCysmM7Uw56y6/tN3XuJwZw/P/Pd3j+mP/uUnX+f7L+6j9r4bRxScJ3r6+M5/1PHAC3X0Dzh3X13Np987j8kTYx++7k790bfYcOAo6/dHbtsb2zj5m+zzyvO5dFYRNdVFvG/RtDHVkApb+Ce5O49veJOv/HoHLe3dvPf8Mj71nnlcVl0U1yCQ8dPR3cdPXj3Id17YQ1NbN9cvKOeLyxdSXZo34mWNtG2nVOD39g/Q0t7NobYumo53cagtcqs/+ha/fb2Jnv4Bbr5wGp+6dh4XzZwc1bpfqTvMfb/Yyu7mDm5aNI0vfWAhFZMnAtDc3sUj6w7yo1f209TWTWXxRD5xRTUfrpnJlEnZwy6zo7uPuiD897R0sKupg+2H2th/+MSpeYomTQjCv5ALKgqpmDyRj33vZT5z/Xw+e8N5UdU+nE0Hj7Hin3/PV+9YHNVZev0DzuMb6vnHp3fS1NbNrRdV8PmbFjCrJLZHCJxLZ3cfm+qPsSH4ANhw4BjH3+old0IGt140nTuXVVJTNfLgS6XAP6mju4+HX9zH935Xx9ETvVw4o5CPLati+cUVFMZg34yMv+2NbTz66kEe31BPe1cfy2YX89kb5nPV3LN/6yJaaRX431yzi031xyIB39ZNa0c3Z5Y7IdMoL8jlmnmlrLx2DnPL8ke8/p6+Ab77uzq+sWYXWRnGp94zl13NHTy1pZHefufa88q466oqrj2vnMwx9Ku2d/Wy81A7rze2sb2xjdcb2thxqJ3uQT9e8tvPXcu88pG/hsHcnXf//XMc6+ylYkou+TlZ5OVkkR/c8gbd507I4Ke19bze2MaSyil8cfkFXFp17ss5jIeBAWfLm8d55NWDrNr4Jp09/cwrz+fOyyr5g0tmUpw3/IfuYKkY+Ce91dPP4xvq+eHL+9lxqJ2crAxuWDiV2xZX8O7zyqLuipTE2Nfaya+2HmLVpga2N7aRnZnBzRdN4+6rqlkag+66tAr8T//7BvY0dzC1MJdphblMnRy5nzY5h6mFuUwtzKV4UnbMdm4dOHyCL/5yK//xRgsFOVncUTOTP7qiijmj+BCJVl//APsOd7KtoQ2AFUtmxGS5z+1o5onNDXR299HZ3U9Hdx+d3X10BLfO7r5T3Sczpkzk8zcv4LZRHPs/Xjq7+1i9uYEfrzvIxoPHyM7M4H2LpvLRZbO4ck7JO7aBVA78k9ydjQeP8fPX3uTJzY0c7uwhd0IG18wr5b0Lyrn2vDJmFo3vNzI5W2//ABv2H+X5N1pYs72JN5oiByosnTWFDy6dwfLF06PeUIlGWgV+Irg7O5vaqSyaRF5O+m49uTtdvQO0d/dSNCl7xDtkE2nHoTYeWXeQn22op62rj6qSSXykppKPX1E1ZF9/OgT+YH39A6zbe4TfbDvEmh3N1B99C4A5pXlcObeEK+aUcPnsYsoLc2O6XjlbX/8A2xvbeWXvYV7ac5iX6w7T2dNPVoZxWXUx119Qzk0XTovbh7ECX0Kjq7efX21tjGz1HzjGS1+4bsiT6tIt8Adzd/a0dPIfb7Tw+92trNv79nX2q0omcemsIpZWFbG0cgrnTytIqQ/2ZHSks4dN9cd4LdjH9NqBo3QGV6idU5rHVfNKuGZeKVfNKx2XfS0KfAml5vYuyguG3qJN58A/U1//AFsb2nh17xFq9x9h/f6jtHZELnWRk5XBwumFLJ4xmUUzJrNoeiHzywuiPmw4bI509rCt4TjbGtrYUn+czW8e4+CRyLepDIPzpxVSU1XEZbOLWVZdzLTJ4/+NKvTH4Us4DRf2YZOVmcGSyiksqZzCf2XOqcNfXzt4jE0Hj7G5/hg/XV/Pwy/tByA7M4N55fksqChgwbQCzp9WyIJpBZQXnH3YcLrq6RugrrWDnYfa2XGonR2NbWxvbD/t8OqZRRO5aMZkPrasiiWVU7ho5mTyU7DLN/UqFpGomRmVxZGzo09eX6l/wE8dKPB6Q+SIsbW7WvnZhrd/gWvyxAmcNzWfeeUFzC/PZ/7UfOaV5zOtMDdlPwi6evvZ29rJ7uYOdjV3sLu5nTeaOtjX2klfcATDhExjblk+V84tYWFF4alDp4tiuKM1kRT4IiGTmWGnzgoffJG9I5097DzUzhtN7exsamdXUztPbWnk+Fu9p+bJy85kbnk+c0rzmBOcVT67NI85pflRXw4gniJXbe2mrqWTutbIuTAnT4ysP3ri1JFpGQaVxZOYX17A+xZO5fxpBZw/rYA5pflp3cUVl8A3s5uA/wtkAt9z9/vjsR6R8ZbObbs4L5sr55Zw5dySU+PcndaOHnY1t0fOJm/uYE9LB+v2HuEXG0//sZ2KybnMLs2jujSPOaV5VJXkMbs08u0illdLdXeOnuhlb2sn+w93sq+1k7rWTvYd7mRf64nTfhw+JyuDOWX5XDRzMrcvncG88nzml+czuzQvpj8skipiHvhmlgn8M3AjUA+8amar3P31WK9LZDyFsW2bGWUFOZQV5Jx1RuhbPZEukrrWSLdIZKu6kyc3n/6tIMNg+pSJVJfkUV06iariyIdCdXAhvqGC19053NnD/sMn2Hcy2A+fYP/hTva2dtLW1Xfa8mcWTWJ2aR41VcWRbxzBN4/pkyfqInSDxGMLfxmw293rAMzsEWAFkLb/FBIaatuDTMzOjPRxTz/7irNHO3siW9yHO9nb+nZoP7Hp9A8DM6gozKWqJI+S/Gz6+p2DR0+w//DpW+onPzRml+bxgSXTqS6JBHpVSR6VxRP1ewtRikfgzwAODnpcD1x+5kxmthJYGTzsMLOdwyyvFGiNaYXjLx1eA6Tu64jVL0nEqm0n0/uY8Fr2AS+9/XDYevYCvx+PgqKoJQGGq2VEbTsegT/U96ezDvZ39weAB865MLPaRB1DHSvp8BogfV7HGMSkbSfT+5hMtUBy1ZOOtcRjd3Q9MPgSjTOBhmHmFUklatuS0uIR+K8C881stpllA3cCq+KwHpHxprYtKS3mXTru3mdmfwb8hsihaw+5+7YxLPKc3T4pIB1eA6TP6xiVGLbtZHofk6kWSK560q6WpLiWjoiIxF/6nlImIiKnUeCLiIRE0ga+md1kZjvNbLeZ3ZvoekbLzPaZ2RYz22hmKXMNaDN7yMyazWzroHHFZvaMme0K7sf+G21p5lzt1sxyzOzRYPorZlY9aNoXgvE7zez941DL58zsdTPbbGZrzKxq0LT+oM1uNLMx75iOopa7zaxl0Dr/y6BpdwVtbpeZ3TXWWqKs5+uDannDzI4Nmhaz92ao/7MzppuZfSOoc7OZXTJo2sjfF3dPuhuRHWJ7gDlANrAJWJjoukb5WvYBpYmuYxR1vxu4BNg6aNxXgXuD4XuBryS6zmS6RdNugT8Fvh0M3wk8GgwvDObPAWYHy8mMcy3vBSYFw586WUvwuGOc35e7gf83xHOLgbrgvigYLop3PWfM/+dEdtDH47056//sjOm3AL8icg7IFcArY3lfknUL/9Qp7O7eA5w8hV3Gibu/ABw5Y/QK4OFg+GHg9nEtKvlF024Hv4ePAddb5HrDK4BH3L3b3fcCu4Plxa0Wd3/O3U8ED18mcl5BPIzl//n9wDPufsTdjwLPADeNcz0fBX48xnUOaZj/s8FWAD/wiJeBKWZWwSjfl2QN/KFOYY/Nr3uPPweeNrP1wSn3qWyquzcCBPflCa4n2UTTbk/N4+59wHGgJMrnxrqWwe4hsiV5Uq6Z1ZrZy2Y21g/2aGv5UNBt8ZiZnTzBLR5ZEPUyg26u2cCzg0bH8r05l+FqHdX7kqzXw4/qFPYUcbW7N5hZOfCMme0IPtUl/UTTboebJ9ZtPurlmdnHgRrg2kGjZwXtdg7wrJltcfc9cazlCeDH7t5tZp8k8i3ouiifG496TroTeMzd+weNi+V7cy4xbS/JuoWfNqewu3tDcN8M/JyxfU1PtKbg6yTBfXOC60k20bTbU/OYWRYwmchX+li3+aiWZ2Y3APcBH3D37pPjB7XbOuB5YGk8a3H3w4PW/13g0pG8jljXM8idnNGdE+P35lyGq3V070usdj7E8kbkm0cdka9SJ3eqLEp0XaN4HXlAwaDhF4GbEl3XCOqv5vSdtn/P6Tttv5roGpPpFk27BT7N6TttfxIML+L0nbZ1jG2nbTS1LCWy83L+GeOLgJxguBTYxRgOmoiylopBwx8EXg6Gi4lcKLMouO0FiuP9dwrmO5/IQRcWr/cmWM5p/2dnTLuV03farhvL+5Lwf5J3eBNuAd4IGuR9ia5nlK9hTtCYNgHbUul1ENmqaQR6iWxN3EOkr3lN0MjXjPUfLx1vQ7Vb4G+JbEED5AI/JbJTdh0wZ9Bz7wuetxO4eRxq+S3QBGwMbquC8VcBW4J2uwW4Zxxq+T/B/8gm4DlgwaDn/ufg/doN/PF4/J2Cx38N3H/G82L63gzzf/ZJ4JPBdCPyozt7gvXVjOV90aUVRERCIln78EVEJMYU+CIiIaHAFxEJiaQ4Dr+0tNSrq6sTXYakqfXr17e6e1mi6xBJtKQI/OrqamprU+a6YpJizGx/omsQSQbq0hERCQkFvohISCjwRURCIin68OOh+t4nR/W8ffffGuNKRESSg7bwRURCQoEvIhISCnwRkZBQ4IuIhIQCX0QkJBT4IiIhocAXEQkJBb6ISEgo8EVEQkKBLyISEgp8EZGQUOCLiISEAl9EJCQU+CIiIaHAFxEJCQW+iEhIKPBFREJCgS8iEhIKfBGRkIg68M0s08xeM7PVwePZZvaKme0ys0fNLDsYnxM83h1Mr45P6SIiMhIj2cL/DLB90OOvAF939/nAUeCeYPw9wFF3nwd8PZhPREQSLKrAN7OZwK3A94LHBlwHPBbM8jBwezC8InhMMP36YH4REUmgaLfw/wn4X8BA8LgEOObufcHjemBGMDwDOAgQTD8ezH8aM1tpZrVmVtvS0jLK8kVEJFrnDHwzWw40u/v6waOHmNWjmPb2CPcH3L3G3WvKysqiKlZEREYvK4p5rgY+YGa3ALlAIZEt/ilmlhVsxc8EGoL564FKoN7MsoDJwJGYVx4n1fc+OeLn7Lv/1jhUIiISW+fcwnf3L7j7THevBu4EnnX3PwSeA+4IZrsL+GUwvCp4TDD9WXc/awtfRETG11iOw/888Dkz202kj/7BYPyDQEkw/nPAvWMrUUREYiGaLp1T3P154PlguA5YNsQ8XcCHY1CbiIjEkM60FREJCQW+iEhIKPBFREJCgS8iEhIKfBGRkFDgi4iEhAJfRCQkFPgiIiExohOvJLWN5jpBoGsFiaQLbeGLiISEAl9EJCQU+CIiIaHAFxEJCQW+iEhIKPBFREJCgS8iEhIKfBGRkNCJVylotCdQiUi4aQtfRCQkFPgiIiGhwBcRCQkFvohISCjwRURCQoEvIhISCnwRkZBQ4IuIhMQ5A9/MKs3sOTPbbmbbzOwzwfhiM3vGzHYF90XBeDOzb5jZbjPbbGaXxPtFiIjIuUWzhd8H/A93vwC4Avi0mS0E7gXWuPt8YE3wGOBmYH5wWwl8K+ZVi4jIiJ3z0gru3gg0BsPtZrYdmAGsAN4TzPYw8Dzw+WD8D9zdgZfNbIqZVQTLkRSk38IVSQ8j6sM3s2pgKfAKMPVkiAf35cFsM4CDg55WH4w7c1krzazWzGpbWlpGXrmIiIxI1IFvZvnA48Bn3b3tnWYdYpyfNcL9AXevcfeasrKyaMsQEZFRiirwzWwCkbD/kbv/LBjdZGYVwfQKoDkYXw9UDnr6TKAhNuWKiMhoRXOUjgEPAtvd/WuDJq0C7gqG7wJ+OWj8J4Kjda4Ajqv/XkQk8aK5Hv7VwB8BW8xsYzDuL4H7gZ+Y2T3AAeDDwbSngFuA3cAJ4I9jWrGIiIxKNEfprGXofnmA64eY34FPj7EuERGJMZ1pKyISEgp8EZGQUOCLiISEAl9EJCQU+CIiIaHAFxEJCQW+iEhIKPBFREJCgS8iEhIKfBGRkIjmWjoSJ6P9YRERkdHQFr6ISEgo8EVEQkKBLyISEgp8EZGQUOCLiISEAl9EJCSS/rBMHbooIhIbSR/4krpG+2G97/5bY1yJiIC6dEREQkNb+DGgbicRSQXawhcRCQkFvohISCjwRURCQoEvIhISCnwRkZBQ4IuIhERcAt/MbjKznWa228zujcc6RERkZGJ+HL6ZZQL/DNwI1AOvmtkqd3891uuS9DSa8xp0dq7IucVjC38ZsNvd69y9B3gEWBGH9YiIyAjE40zbGcDBQY/rgcvPnMnMVgIrg4cdZrYzDrVEqxRoTeD6o5UKdSakRvvKO06uGqcyRJJaPALfhhjnZ41wfwB4IA7rHzEzq3X3mkTXcS6pUGcq1CgSVvHo0qkHKgc9ngk0xGE9IiIyAvEI/FeB+WY228yygTuBVXFYj4iIjEDMu3Tcvc/M/gz4DZAJPOTu22K9nhhLiq6lKKRCnalQo0gomftZ3esiIpKGdKatiEhIKPBFREIiLQP/XJd2MLMcM3s0mP6KmVUPmvaFYPxOM3v/oPEPmVmzmW1NxhrNrNLMnjOz7Wa2zcw+k6R15prZOjPbFNT5N7GoU0Si4O5pdSOyo3gPMAfIBjYBC8+Y50+BbwfDdwKPBsMLg/lzgNnBcjKDae8GLgG2JmONQAVwSTBxt7oxAAABuUlEQVRPAfDGmctMkjoNyA/mmQC8AlyR6Hajm25huKXjFn40l3ZYATwcDD8GXG9mFox/xN273X0vsDtYHu7+AnAkWWt090Z33xDU2g5sJ3LWc7LV6e7eEcw/IbjpyAGRcZCOgT/UpR3ODL5T87h7H3AcKInyuUlfY9CtspTI1nPS1WlmmWa2EWgGnnH3sdYpIlFIx8CP5tIOw80T1WUhYiBuNZpZPvA48Fl3bxt1he9cQzTzDPtcd+939yVEzsJeZmYXjqlKEYlKOgZ+NJd2ODWPmWUBk4l014zXZSHiUqOZTSAS9j9y958la50nufsx4HngphjUKiLnkI6BH82lHVYBdwXDdwDPursH4+8MjjyZDcwH1qVCjUG/+YPAdnf/WhLXWWZmUwDMbCJwA7AjRvWKyDtJ9F7jeNyAW4gcpbIHuC8Y97fAB4LhXOCnRHYkrgPmDHrufcHzdgI3Dxr/Y6AR6CWy9XpPMtUIXEOky2QzsDG43ZJs7yWwGHgtqHMr8FeJbi+66RaWmy6tICISEunYpSMiIkNQ4IuIhIQCX0QkJBT4IiIhocAXEQkJBb6ISEgo8EVEQuL/A9YR9W0k4T1oAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x1a235560f0>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "%matplotlib inline\n",
    "import numpy\n",
    "\n",
    "rand = numpy.random.RandomState(seed=1234567)\n",
    "\n",
    "size = 1024\n",
    "factor = 128\n",
    "xval = sorted(\n",
    "    list(rand.normal(1e-3, 2e-4, size))\n",
    "  + list(rand.normal(2.5e-3, 4e-4, size)))\n",
    "\n",
    "ty = []\n",
    "tx = range(int(len(xval)/factor) - 1)\n",
    "for i in tx:\n",
    "    r1 = xval[i*factor : (i+1)*factor]\n",
    "    r2 = xval[(i+1)*factor : (i+2)*factor]\n",
    "\n",
    "    u1 = numpy.mean(r1)\n",
    "    u2 = numpy.mean(r2)\n",
    "    o1 = numpy.std(r1)\n",
    "    o2 = numpy.std(r2)\n",
    "\n",
    "    d = abs(u1 - u2) / (numpy.sqrt(4 * o1 * o2))\n",
    "    dd = numpy.sqrt(d*d - 1)\n",
    "    t = 2 * numpy.log(d - dd) + 2 * d * dd\n",
    "    ty.append(t)\n",
    "\n",
    "ly = []\n",
    "lx = numpy.linspace(0.0, 1.0, 100)\n",
    "for p in lx:\n",
    "    l = numpy.log(1 - p) - numpy.log(p)\n",
    "    ly.append(numpy.abs(l))\n",
    "    \n",
    "ax = matplotlib.pyplot.subplot(2, 2, 1)\n",
    "ax.set_ylim(0, 20)\n",
    "matplotlib.pyplot.plot(tx, ty)\n",
    "\n",
    "ax = matplotlib.pyplot.subplot(2, 2, 2)\n",
    "ax.set_ylim(0, 20)\n",
    "matplotlib.pyplot.plot(lx, ly)\n",
    "\n",
    "ax = matplotlib.pyplot.subplot(2, 2, 3)\n",
    "matplotlib.pyplot.hist(xval)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 120,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAX0AAAD8CAYAAACb4nSYAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvNQv5yAAAIABJREFUeJzt3XmQHNd94Pnvr+6qvu9u3DfAEyTQAEVRkklRPKWxPGs5hl5Zkr3ape2xdseznt3wOHbHu5rxhh0xY48teczh2l7JM7bssSRKtMxT1G1LBLpB8D5wNYA+0N1o9F135m//yOzq6gtodPWFrt8nIiOPyqp8lZX1ey9fvnwpqooxxpjyEFjrBBhjjFk9FvSNMaaMWNA3xpgyYkHfGGPKiAV9Y4wpIxb0jTGmjFwz6IvIn4vIgIi8UbSsXkReFJFT/rhugfd+xl/nlIh8ZjkTbkyp7Ng25WgxJf0vAQ/PWvabwEuquhd4yZ+fQUTqgd8G7gKOAr+90B/ImDXyJezYNmXmmkFfVX8AXJm1+OPAl/3pLwM/M89bHwJeVNUrqjoMvMjcP5gxa8aObVOOQkt8X4uq9gGoap+INM+zzmbgYtF8t79sDhF5HHgcoKKi4vCBAweWmCxjrq2zs/OyqjYt8LId2+aGdI3jumCpQX8xZJ5l8/b5oKpPAk8CtLe3a0dHxwomy5Q7ETlf6kfMs8yObbOmFntcL7X1Tr+ItPkbagMG5lmnG9haNL8F6F3i9oxZLXZsmw1tqUH/aWCqxcJngG/Os87zwIMiUudf5HrQX2bMembHttnQFtNk8yvAj4H9ItItIp8Ffhd4QEROAQ/484hIu4j8KYCqXgH+LXDcHz7vLzNmXbBj25QjWW9dK1u9p1lpItKpqu2rvV07ts1KWuxxbXfkGmNMGbGgb4wxZcSCvjHGlBEL+sYYU0Ys6BtjTBmxoG+MMWXEgr4xxpQRC/rGGFNGLOgbY0wZsaBvjDFlxIK+McaUEQv6xhhTRizoG2NMGbGgb4wxZcSCvjHGlBEL+sYYU0aWHPRFZL+InCwaxkTk12etc6+IjBat82+Wur1MKk96IrfUt5syMjaUQt2lPxxotY9tY1ZTaKlvVNV3gTsARCQI9ABPzbPqD1X1Y0vdzpST377AyW9f5LYPbeaOB7aRqI6U+pFmgxnpT9L5XBfvvtzPg5+9hT2Hm5f0Oat9bBuzmpYc9Ge5HzijqueX6fPm2HO4mdGBFCe/fYHXv9fNzR/cxJ0PbKeyLrpSmzQ3iKHeCTqfPc/pjn4CoQC3/dRmWnfVLNfHr/ixbcxqWq6g/xjwlQVeu1tEXgV6gX+lqm/OXkFEHgceB9i2bdu8H9KwqZIHP3sLRz+2k85nu3j9ez288YMebnr/Jg49tI3qhvjyfBNzwxi8ME7Hs12cfWWQUDTIHR/ZxsGPbKWiZlkLAit+bBuzmkp+MLqIRPAO+ltUtX/Wa9WAq6oTIvIo8Iequvdqn7fYh0ePDqY48cJ53vnHPlDY975WDj+8ndrmRAnfxtwILp0bpfOZLrpeHyISC3LbfVs4eP9W4pWLq/Jb7AOk1+rYNmYpFntcL0dJ/xHgxOw/BYCqjhVNPyMi/0lEGlX1cqkbrWmKc98nD9D+yA5eefECb/2ol3d/3MfeIy0cfngH9ZsqSt2EWWd6Tw3T8UwXF98eJloR4ug/2cnt920hmgiv1CbX5Ng2ZiUtR9D/eRY4/RWRVqBfVVVEjuK1Fhpahm0WVNXH+NA/28fhh7dz8tsXeeP73bx3vJ/ddzbR/ugOGrdULefmzCpTVbrf8YJ976kR4lVh7v6nu7n1pzYTiS1X7eSC1vTYNmYllPSvEZEE8ADwy0XLfgVAVZ8APgH8qojkgRTwmJZan7SAipoo9/zsHg49tI1Xv32R177XzZkTg+y4vZH2R3fQsqN6JTZrVoiqcv6NITqe6aL/3BiJmggf+Lm93PzBTYQjwRXf/no6to1ZTiXX6S+35ar3TE/meP173bz60kUyyTzbbq7n8KM72LSndhlSaVaKusrZVwfpeKaLyxcnqKyPcvih7Rx4fxuh8PIE+8XWfS43q9M3K2k16/TXpVhFmCMf3cnBD2/l9e93c/LbF3nq359g875a2h/dweb9dYjIWifT+FxXOdM5QMezXVzpnaS6Kc59nzrA/rtaCYbsxnFjlsuGDfpTIvEQhx/ewe33beXNH/bwyosX+OZ/PEnrrhraH93BtlvqLfivIcdxOXWsn87nzjPSn6SuNcFHfulm9rY3EwhasDdmuW34oD8l7LfjvvWnNvP2P/Rx4vnzfOuLr9K0rYr2R3ew8/ZGJGDBf7U4OZd3fuL9DmOX0zRsruSh/+lWdt/ZZL+DMSuobIL+lFA4yG33buHmD2zi3Z9covO5Lp594nUaNldw+JEd7D7UTMCCzorJZx3e+oc+XnnhPBPDGZq3V/GBn9vLjtss0zVmNZRd0J8SDAW4+QObOHB3K6eO99Px7Hle+NM3qWs9x+GHt7P3SItVLyyjXMbhjR/0cPLFCyTHsrTtruG+XzjA1putes2Y1VS2QX9KIBhg//va2Hu0lTMnBuh8totvf+ltjn3rHIcf3sH+99mFxFJkU/nChfT0RI7N++t44LO3sHlfrQV7Y9ZA2Qf9KYGAsLe9hT2Hmjn32mU6nuniu//1HY4/c45DD27npnuWr8lgOUhP5njtOxd57bvdXpPZW+ppf3QnbbuXrSM0Y8wSWNCfRQLCrjua2HmwkQtvXqHjmXP84K/fo+PZLu58YBu3fHAz4agF/4Ukx7K8+tJFXv9+N7m0w86D3s1xzdvt5jhj1gML+gsQEbbf2sC2W+rpedfrBuAfvnqaE8+f5+D9W7nt3i2r0Q3ADWNyNMMrL1zgzR/0kM+77DnUzOFHdtC4pXKtk2aMKWJR6xpEhC0H6tlyoJ7e0yN0PtPFT75xlldeuFAI/rGKFevwa90bv5LmlefP89Y/9OE6LnuP+h3etVmHd8asRxb0r8OmPbVs+l/uoP/cGB3PdnHs785x8sUL3HbvFg5+ZPFd+24Eo4MpTjzXxTs/uQQK++9u5dBD1rW1MeudBf0laNlZzUf/+e0MXhyn89kuOp8/z6vfucit/qMcl/khHuvK8KVJOp89z3vH+5EA3PyBTdz5oD3ExpgbxQ0T9NNvvUW2q4vEkSOEmprWOjkANG2t4uHHb+NK7ySdz3V5FzC/11MIhFX1sbVO4rIZ6pmg45kuTp8YIBQKcPt9W7jzgW1U1K6fDE5VyV28SPL4cSruuYdwa+taJ8mYdeeGCfqjT/8dV770JQAiO3eSOHKExNGjJI4cIdyytAdgL5f6TRU88D/cwpGP7uTE8+d58wc9vPnDHg68v43DD22nuvHGLQUPnB+j45kuzr16mXA0yKEHt3Hw/vXxYHpVJXfhApPHjpE8fpzksePkL10CoO13/h21P/uza5xCY9afG6ZrZc3nSb/9DsmpP3hHB+7EBACR7du9DODoURJHjxBuaVntZM8wdjnFiRcu8PY/9qIu7D/awuFHdlDbcuPUd186O8rxv+/iwptDROIhbv/wFg7et5VY5dpdtFZVsl1dhQCfPHaM/MAAAMHGRiqOHikUBiK7di1485d1rWw2okU/BvRGCfqzqeN4mcBx78+f7OjAHR8HILx9G4kjR6iYOhNoa1vpZM9rYjjDKy+e580f9uLmXfYc9poxNmxen80YVZXe90Y4/kwXPe8OE6sIc/AjXgulaHz1TwpVley5Lu/39TP7/OAgAMGmRiqOTGf0kZ07F32HrwV9sxGtStAXkS5gHHCA/OwNivcv/EPgUSAJ/KKqnrjaZy71j6GOQ+bdd/1T/Q4vExgdBSC8dSsJvxRYcfQo4U2brvvzS5Ecy3Ly2xd4/fs95DMOu+7wHuXYtG19PMpRVbn41hU6nu2i7/Qo8eoId35kG7d8aNOq3ougqmTPni0E+Mljx3Eue4+cDTU3F6rzEkePENmxY8ndOCzmz7Gejm1jFmM1H6Jy31UeBv0IsNcf7gL+xB8vOwkGid18M7Gbb6bhF38RdV0y771H8tgxJo8dY+LbLzH6ta8DEN6yZcY1gciWzSuRpIJEdYT3/3d7OPTgdl79zkVe+85Fzp4cZPttDbQ/uoPWnWvTNYGq0vW690jCga4xKmqjfPCf7eXmezYRWoVHEqoq2TNnvIz62HGSx4/jDHmPmQ21tFBx990kjvoZ9bZta9FXz7o4to1ZTitdjPs48Bf+s0N/IiK1ItKmqn0rvF0kECB24ACxAweo//SnvUzg1CmSL3ulyInvfpfRp54CILxp03Qp8q6jhDdvXpEAE6sMc9dP7+KOj2zl9e91c/Kli3zt9zrZcqCOIx/dwaa9dcu+zfmoq5x5ZZCOZ7sY6p6gqiHGvZ/cz4H3tREMr1zncuq6ZE6fLgT45PHjOFeuABBqa6PyA/cUMuPw1q3rvUO2NTu2jSlFqUFfgRdERIH/rKpPznp9M3CxaL7bXzbjjyEijwOPA2zbtq3EJM1PAgFi+/cT27+f+k9/amYAOnaMie9/n9FvfAOA0KY2KqbOBI4eJbxly7IGoGgiTPujO7n9w1sL3Q0/9R9eoW1PDUce3cmWm1bmUY6u43Kqw+tJdPhSkprmOB/+9E3su6uF4Ap0I13IaP19nDx+HGdkBPD2ceUHP+jt4xXMaEuwasd253Nd1LdV2DMFzIKmqmD7zoxy10/vKumzSg3696hqr4g0Ay+KyDuq+oOi1+c7gudcRPD/UE+CV+8534b+9r2/5VtnvsXu2t3srt3Nnto97K7dTUOsYUnBQgIBYvv2Edu3j/pf+KRX1XD6dOGawMQPf8ToN58GINTaOvOawDJVNURiIQ49uJ3b7t3CWz/s5ZUXzvP0H52kZWc17Y/uYPutS/tuszmOy3svX6Lz2fOMDqaoa6vggc/ezJ7DLcv6wJjiKrWpFjbO1HWVzZupvO++6dY1y1SlpqpcmrzE6ZHTnBk5Uxj/y8P/kqNtR0v56FU5tp2cy9v/2MfogPe7HHpoG3uPrEwmbG48rqucOTHAiefPc/niBJV1Ue54YFtJDStKCvqq2uuPB0TkKeAoUPzH6Aa2Fs1vAXqXsq1oMIqiPNf1HOPZ8cLy2mgtu2p2FTKBQmYQb7iuzxcRonv3Et27l/pPfnLGRcXJY8eY/McfM/b03wFFFxX9jKCUi4oA4UiQg/dv5dYPbebtH/dx4rnz/P0fv0bj1kraH93BroNLe4Sgk3MLnzd+JU3j1koefvxWdt2xPI8knLp4PnXRdfbF88r77/fq5I8cIby5tCA/FdzPjJ6ZEdzPjJwhmU8W1muMN7K7ZndJ2/K3tyrHdjAc4L//7bs43en9sV/60tu8/M2z3PHANm6+Z5P16Fqm8jmHd358iVdevMDYYIralgT3feoA+4+2llwFu+TWOyJSAQRUddyffhH4vKo+V7TOR4HP4bVwuAv4I1W9avHrWi0cVJXB1GDhD1/85x/PTWcGddE6dtXOzQzqY/VL+r5e88FzfvPB40weP4Yz6LcsaWoqallylMjO0jKB2SXz+k0VHH5k+6JL5rmsUzhzmBzN0ryjmiOP7mD7baWdOajjkH7nnenqms5O3LExAMLbthUCfOLIkSW3kFJV+pP9M37XMyNnODN6hsncZGG9hljDnLO+3TW7qY3VXnMb12rlsJbH9vk3hjjx/Hn6To8Sqwhz+4e3cNtPbVnT+yPM6smk8rz5gx5OvnSR1FiW5u1VHHp4OzsPNl3zv7/iTTZFZBfwlD8bAv5KVX9HRH4FQFWf8Ju1fRF4GK9Z2y+p6lXbrC25yaafGRQHi6npidxEYb36WD27anbNCBZ7avdQF7u+i6iFG4WK6qsLNwo1NU5fEzhy5Ko3Cl2N67ic7hyg49nzDPdNUtMc5/DDOxasg8+m84VrBKnxXMnXCObcENfZWbgXwrshruiu6Ovs8mAquBcH9dMjpzk7cnbO77Wnds+cs7nFBPeFLCLor/mx3Xd6hBMvXKDrtcuEokFu+cAm7vjIVirrNk7XHmba5GiG177TzRvf7yabdth6Ux2HHtrO5v2L/+9u+JuzFktVGUgOTGcCVwkuU6XFpQQXVSV3/jyTxXeL9vcDEGxoKLQvrzh6lMju3dcVhGe3tqlujHHooe2F1jaZVJ7Xv3uRky9dJDOZX3JrIC/Iv12o0kp1npi+63nHjum7nq+j64vZ+//s6NnC/i8+Myve/8Ul+OvNjBfjRro5a6hngldeuOB1cCew765WDj24jbpW67p6IxgdTPLKixd55x/7cByX3Xc2c+ihbUt66JAF/WsoLmkWzg78+uLiaoSpkub1ViMUOv+autHo5WOFfmGC9fVeJuBnBNE9e5DAtevpCu3q//4cA+fHqayLsv3WBk51DJBN5b12/4/soHXX4tr9ay5H+q23Cn3XpDpP4E563z2ya1chfYkjRwg3Xz3Iz3emtVC125xqmRKq3ZbiRgr6U8Yupzj50kXe/lEv+bzLroNNHHpoOy077YlkN6LL3eOceP4Cpzv6kaBw4H1t3PnAtpK6atl4Qf/l/wydX4bqTVDdBlWb/OlNUNXmjeN1UGJrl+upUy4OWlNBrCY6f8BVVXLd3YWzgMnjx8j3eq37gnV1JNrbCyXp6N6rZwKqysW3r9DxjHcH7WLv8NVcjtQbb3h3LB87RurECdykdxE0snt34Uwk0d6+YE+mqsrl1OU5rWXOjJ6Zc4F9dhXarppd132BfV6uC8nLMNbrDeO9MNY3c/rh/wf2fGTet9+IQX9KajzLa9/t5vXvec8e3ry/lkMPbWfrTfXrrcmrmUVV6Ts9QudzF7jw5hDhaJBbP7SZg/dvXZbeajde0H/ja/D6V2Gsx/tTTw4yp4VcKOZnAJv9jKF42s8gKlsgeP2NloqbBk5VUSzYeqQoyE0Fvfkyg2x3z4x+ZXI9PQAEa2oKJezE0aNE9+2bNxNQVfI5l/ACd89qNkvqjTen6+RfeQX1g3x0757pu5Lb2wk1Ns757KH00NzgPnKGsexYYb2aaM2cKrGpkvuSglA+A+N9fhDvmX96vA/c3Mz3SRCqWv3fvA3u/hxse9+8m7iRg/6UbDrPWz/q5eSLF5gczdK4tZJDD21n96HmZW2Ga0qnrtL1+mVOPH+eS2fHiFeFuf3DXmu95Xzq3sYL+rPlszBxyQ8CvUWlvr6Z00525vskABXNc88SCtN+JhFZXJ2pqtI32Veopy4u+abyqcJ6TfGm6cygqFVRdWT69DzX0zPjmkCuuxuAQE0Nifb2Qi+S0f37keDcQO9ms6Rff326E7pXTqIpLw3RvXuL6uTbCdXXF9I/lB6aE9jPjJ5hNDNa+OzqSPWMM5vrvk9CFdKjc3+f2aX15Dy9HoQT8/xW/hnf1HRlMwQW17xxIwT9KU7O5d1jl3jlhQuM9CepaYpz54Pb2P++VkJha+65lhzH5dTxfk48f4HhvkmqGmLc+cA2bnp/24p0c7Lxg/5iqEJyqCjA9MzKJPzp9Ojc90Zr5qlK8jOFqeCTaFiwOslVl77Jvjktic6Onp2RGTTHm+dUEe2u3U1VpIpcb693PcDPCHIXLgAQqK72qoOOHCG6ZzcpP9CnXjmJptNe8vfvLwT4xJEjBGtrGUoPzciYpi5sFwf3qkjVnFL7nto9Vw/urgMTA/NXsxSX0IuqxwoSDbOC+Obp0vrUdKym5Gq7Yust6Pf8q/+NQEUFicOHSBw+TGjTpus+S3Jd5dzJQU48f56B8+MkqiOFez8ia9BDajnLZRze+gfvLGxiOEPD5goOPbSdPYebCSzhprtcfz/Jjg5SnSfIXbrE1v/0x/Out+GC/tc6u3n61V6q42Fq4iGqY2Fq4mF/Plw0771WHQ8TXOxpbnby2mcME/2g7sz3BaNelcK8Zw1+EKtshdD0A0dcdemd6J1RRXR65DTnRs/NzAwSzYWWLFPBd3u6ksCrbxfuFcieP++tLEL0wAESR9qpOHqU9C27OMflOTcyjWRGCp8/FdxnN4VsjDfODDi5VNH+KA7iRftm/BKoM3PfBELT+2PO2VTRdGhxdZk5x2UslWM0lWMsnffGhfkcY6l80XSOz923h7t2zX/9YD0FfXVdun/tcySPHy+0lAq1tpI4fJj44UMkDrdf8xrPjM9TpfvdYU48d57ud4aJxEPsOdzM5n21bNpbR2Xd+nnS2UaSHMvS894wvadGON0xQHrSazJ96KHt13VnfeGm0I5OUic6SXZ0Fqp9JZEgcccdbHniTwhE5j7EaMMF/b96+QJ/ffxC4c8+ls7juFdPe2U0RE08TFUsNCeDqI6HZmUW05lGTTxMPByc+UM5eS/wX+2MYawX8ulZqRCoaFq4FOsHQDdSQe9E74wgPZUZpJ3pz2xJtBTOCg44TbQMOlxoDfKu01M4qxjODBfWrwpXzTiTmAruTbFGJD2ycDXL1HRqmDkiVTOrVuZMb4ZEIxQFKlVlIjMVrGcG6KlAPjYrkE+tO5bOkcw6c9NRJByUot82zG88uI8P7p3/YvR6CvpT1HG8Liw6TxT+7FP3fQSqq0nceSfxw4dJtB8mduut8/7pZxs4P8YrL17gwhtDZNPe/qtpirNpXy2b99ayaV/dhnqk52qaHMnQc2qY3vdG6D01wvAl71pZKBpk64E67nhgG5v2XLu5t2azXjPpjk6SJ06Q6uws9E8VbGggcfgwicOHiB9uJ3ZgPxJa+KxtwwX92VSVyawzHSRmBY85gWNWkJm8RhAJBaSQGVTHQt54nrOK4kBTEwtRwwRVuUHCk5fmqVbyp68jkLqVrfSEw5xxU5xOD3Jm9CxnRs7MyQwqw5XTZwXVO9kdrWO3xGjOJJHxS/NfFL3ODCpb0cJoqIlRN160b6cz4WuVvq+RR1MVC827b4t/h5rEPBl1LEwsHLhhH6LyfNfzVEWq2FWzi5ZECyLitfbq6Smc1ic7O8mePeulPxIhdvttJA55mUD8jjsIVi/cdNN1lcsXx+k9NULPeyP0nR4hk8wDUN0YY9Ne7yxg875aqhpi1gpoHuNX0vSeGqH3vWF6To0wOuCdlYdjQdp213pnUvtqadpWddV+k5yJSVKvniTV2ell8K++WqiSjWzf7mXsfjVfePv2wrEwmBrk7OhZLqcu87FdH5v3szd80C9V3nHnZBDFQWp2RjKayjFe9FrOufp+S0SCC1Y71UccWgMjNOllGtwhanKXqcwOkEj3E0n2E5jo8wL1nCqTcKE6yalspbeilkvhMFszGVomryBTJfbJgXmqoiJoVRtOZRuZRCupaDPjkSZGwk0MSQOD1NPn1jKSZcFqlEx+1mfOEg0FigLxdZxdxcJUxkKLr44r0XoL+h/66w8Vzs4SoQQ7a3ays2Ynu2p2satmFztrd7K1aisyMk7qxAmSnSdInugk/eZbkM971Xv795M4dMjLBA4fvuojQ11XGeqZKJRSe04Nk5n0MoHKuiib99WxaV8tm/bWUtMUL8tMYOxyyt83XqAfu+wH5njIzyS9QN+4pfKq9fT5y5e936vTy7zT77wDjgOBALGbbipU4SUO3QkNdXSPd3Nu9BxnR89ydvQs50bPcW70XOFG0nAgzLFPHiMUmFvit6C/glSVdM6dk0GMpXOMJq9d6h1P56/6+QGB2liA7bEkOyOjbA2N0BYYpoUhGt0r1DmXqcoOUJEZIOSkyIarmYw0MRZuYjjUyGVp4JLW0evWcyFXS1eumvPpOBMZh6v93AHBL1FfLUBPn/XMztBiN0hrkfUW9P/Fv/9denWC8USedHSSfHiQjPSR1iuFdQIEaYptZkvldnbW7GJv3S5uqdzK1u4kzsk3vLOBkycLTXLDW7b41QKHid10E+EtWwjW1s4bwNVVrvRN0vPeCL2nvHrp1LjXJLaiJsKmfXX+NYFaalsSGy4TUFXGLqe87/+elwlOXMkAEK0IsWlPrZcR7q2lYUvlgk1indFRst3dZN55l2RnJ6nOzsJ1N4nFiB88SOLwIQIHb6VvRzVvTvbw3vBpzo120T3RxUC6B0enmyJHpY6othLKNxNLV1KVCtKQjfHEv/4/5v0NLOivY46rTKTnq9fOXbXOeyojmS5xK2EcckWdpSYiwTnBuFA9MqsEPrsqpSISKos23ust6E/8zm4qc15TVRdhRGoYoJ5zWs2bwQrOhsP0RpShSI6JyAT5yBjI9FmXOLVE3FaqtJUDQ3FuvZRhT88QLV2niY5Pt8xyExVoaxuyaTPBzVsIb9lCbNtWKndspWLbNkJxr35fVRnuS9J7atgv6Y6QHPOaPieqI0Ul3TrqWhM33DMAVJXRgRQ97w37Gd0IkyNekI9Vhv3rHV6VV8OmisL3c9IZJi52M9l1gfSFC+S6u8l396C9PcilXgKT0926ZCuqGNixhzObm3izLcabDSnGAv1k5RIaulKUGCGYq6YyU0l9LsKmrLArm+NmJ8VuGaWZYep0hADe752TCKH/s3/eC/sW9DewdM4pnD0ks/miC9ZhIiHrh/1a1lvQp++16Yvp45e8C+jj3j0oOt6HzLp3IQd0ReK8Ha/nvUiC0+EwF4Iu/YEM2cB0laDmo7T017N5MELzCLSMOrSMZWkeT9I6MU7UmVl9eDlew+XKBoaqGxmpaWK0tpnJ+hZS9S0E4rXUJaF6PE9sOE8w42c6AkQCSCyIRIMEYt50KBYiEAsSigcJJUKEYiHCiRChRIhILEgoGCAQEEIBISBCMCAEAxAMBAiKEAiAIDiu4qh6Y39wi+bzjksu45BL5skl8+RTefJJh3w6j5tycNIObjqPm3HRtINmHMi4hfs63UiAVF2IiaogVyqEfGaU+FA/lVf6qRoZpG5kgIbxyzSOD9GQGiVQdENoNhCkv7KK/qo4/dVR+muCDNRCT1OOSy3DSHi6NV5IA7Q4MbY5AXbl8xzITHJTepidmSSzL8lrvB6p3uTfbNjqXeObaiVY1Qptd8zbhHk1n5FrVlksHCQWDtK8Pp6rbkrVdrs3zEPAu0t5ot/PCHoJj19i73gve8cvFWUUfWh2gv5gkLORMOfCIc6Gw5xtnOBSa5h3ApAqjhMKtZNBmkdg60iIzSNBWkaKSTiUAAAapUlEQVQdGof7uWWwh+ozGYqLD7lgiMuVDQxWNtCXaKC3YhPp+BYIVxPMhAllw4QJESFAjAAxFYLzPmcGXJSUQErUGwJF0wKpgDetQEKFmAtxFW9ahbg/7w1cczsZXDI45MiT1zyOZiE3QSzZTXOyh61vDtE8OUTjxBCR/Mw7vccqI1ypC3O+RTheE6e3Ns+FWof+WhipBJUJwCvhxxQaXNicz/OBVJpdYzl25byhNe8QCCemW+y1LBDQK1uR8Mq2qLKgb8x6F4pC7TZvuArJjNM6fonWsV7eP3XGMNYHqSuQHiOZGWEoM8pQboIhJ8WQm2GoIsBQdZDenUFeDwYYCgYZCgbJuEGaRqF5VL2zhBGH5pEB2kYG2HtBiWeunmQnKOTilWTi1eTiVWRjNeSi1eQiVeTDleRDleSCCfKBBHmNk9coeQ2jLHSmqoQkS1gyhDTlDflJwvkJQrlJwtkxwtlxwukxIqlRoslxwqkJgs7VGx+kIzBYK3Q3Q+c+GKgNMFAD/XXCYA2EAnkanAwNjkuj49DgOOx2XBoyDg3ZCA2hOA2hShqjNSRitRCt9voAm7qpszigR6uX9SbDpbKgb8xGEa3yhsa9876c8IfC475cBzLjkBmD9Jh3Z7o/nUoNMTR5iaHUEEOZYYayXmbRkU9xxc0wkc4hI0ogDcG8EMhDKAeRvBDLQjSnRHMT/gCxHETHIJJTEjmI5iCahWjeGwfwalycYIxsuJJcuBIEwrlJwrkJQvk0guIC2TCkw9PjTBjGI0ImDJkEZGqml2fCAX8MTkhx/LEbBWqERDxEQzBGQyjOjlAFhyM1NMTqaIg10lDRQiLR4AXrWLV3Z/jUdLR60d1+rDdLDvoishX4C6AVcIEnVfUPZ61zL/BN4Jy/6Ouq+vmlbtOY1VA2x3YgCPFab5gljvf8xy3X+gxVL/Nwc2g+S95Jk80myeQnyeZSZPMpsrkUGSdFLpcm46TI5jNk8ylGnCwZJ002nyGXTuGkUuRTaTSTxkllUZRgtJpAPEYwFiMUTxCOxYiEYkRCUaLBKPFQjNpQjEgwTiQcIxKMEQ0nvHXCCSKhBNFIBaFgDAlFvGbPgeC6KHGvlVJK+nngN1T1hIhUAZ0i8qKqvjVrvR+q6vx3ExizPtmxvVgiXq+1wRASjhOmhjBgj3hZv5bc1ENV+1T1hD89DrwNlPb0a2PWATu2zUa2LO37RGQHcCfw8jwv3y0ir4rIsyJyywLvf1xEOkSkY3BwcDmSZMyysGPbbDQlB30RqQS+Bvy6qo7NevkEsF1VDwJfAL4x32eo6pOq2q6q7U0LPLHJmNVmx7bZiEoK+iISxvtT/KWqfn3266o6pqoT/vQzQFhEGmevZ8x6Y8e22aiWHPTF6/zhz4C3VfX3F1in1V8PETnqb29oqds0ZjXYsW02slJa79wDfAp4XURO+st+C9gGoKpPAJ8AflVE8kAKeEzXW78Pxsxlx7bZsJYc9FX1R7DA/c/T63wR+OJSt2HMWrBj22xk1juXMcaUEQv6xhhTRizoG2NMGbGgb4wxZcSCvjHGlBEL+sYYU0Ys6BtjTBmxoG+MMWXEgr4xxpQRC/rGGFNGLOgbY0wZsaBvjDFlxIK+McaUEQv6xhhTRizoG2NMGbGgb4wxZaTUZ+Q+LCLvishpEfnNeV6Pisjf+K+/LCI7StmeMavFjm2zUZXyjNwg8MfAI8DNwM+LyM2zVvssMKyqe4A/AH5vqdszZrXYsW02slJK+keB06p6VlWzwF8DH5+1zseBL/vTXwXun3qYtDHrmB3bZsMq5cHom4GLRfPdwF0LraOqeREZBRqAy8UricjjwOP+7ISIvLvANhtnv3cNrZe0rJd0wI2Tlu3XeO9aHNsraT39LtfD0n19rnVcA6UF/flKNbqEdVDVJ4Enr7lBkQ5VbV9c8lbWeknLekkHbKi0rPqxvZLW0+9yPSzdK6OU6p1uYGvR/Bagd6F1RCQE1ABXStimMavBjm2zYZUS9I8De0Vkp4hEgMeAp2et8zTwGX/6E8B3VHVOaciYdcaObbNhLbl6x6/H/BzwPBAE/lxV3xSRzwMdqvo08GfAfxGR03iloMdKTO+anibPsl7Ssl7SARskLWt0bK+k9fS7XA9L9woQK5wYY0z5sDtyjTGmjFjQN8aYMrLugr6I/LmIDIjIGwu8LiLyR/7t76+JyKE1TMu9IjIqIif94d+sUDq2ish3ReRtEXlTRP7FPOusyn5ZZFpWfL+ISExEjonIq346/u951inbrhIW8zutZyISFJFXRORba52W6yEitSLyVRF5x9/3d691muZQ1XU1AB8CDgFvLPD6o8CzeO2k3we8vIZpuRf41irskzbgkD9dBbwH3LwW+2WRaVnx/eJ/z0p/Ogy8DLxv1jr/HHjCn34M+JuV/q3Wy7CY32k9D8D/CvzVavy/ljndXwb+R386AtSudZpmD+uupK+qP+Dq7Z0/DvyFen4C1IpI2xqlZVWoap+qnvCnx4G38e4ILbYq+2WRaVlx/vec8GfD/jC7VULZdpWwXn6npRCRLcBHgT9d67RcDxGpxiso/hmAqmZVdWRtUzXXugv6izDfLfJreTDf7VcxPCsit6z0xvwqijvxSrbFVn2/XCUtsAr7xa8COAkMAC+q6oL7RFXzwFRXCWXlGr/TevQfgf8dcNc6IddpFzAI/H9+1dSfikjFWidqthsx6C/q9vdVcgLYrqoHgS8A31jJjYlIJfA14NdVdWz2y/O8ZcX2yzXSsir7RVUdVb0D747ZoyJy6+xkzve2lUjLenWN32ndEZGPAQOq2rnWaVmCEF518J+o6p3AJDCnW+61diMG/cXcIr8qVHVsqopBVZ8BwiLSuBLbEpEw3p/3L1X16/Ossmr75VppWc394m9jBPge8PCsl8q6q4RFHDPr0T3AT4tIF17vph8Wkf+6tklatG6gu+iM86t4mcC6ciMG/aeBT/utVd4HjKpq31okRERap+qIReQo3v4cWoHtCF494duq+vsLrLYq+2UxaVmN/SIiTSJS60/HgY8A78xarWy7SljkMbPuqOq/VtUtqroD7+L7d1T1F9Y4WYuiqpeAiyKy3190P/DWGiZpXqX0srkiROQreK0/GkWkG/htvIt0qOoTwDN4LVVOA0ngl9YwLZ8AflVE8kAKeGyFgso9wKeA1/06bIDfArYVpWW19sti0rIa+6UN+LJ4DzwJAP9NVb8lN25XCctt3t/JP/MyK+d/Bv5SvD6bzrKC8WmprBsGY4wpI9es3pF5blASkXoReVFETvnjugXe+xl/nVMi8pn51jHGGLN6FlOn/yXmXiD7TeAlVd0LvMQ8V6hFpB6vOuQuvMfP/fZCmYMxxpjVcc2gv8ANSsU3vXwZ+Jl53voQXtvpK6o6DLzI3MzDGGPMKlrqhdyWqZYhqtonIs3zrLPom4Wk6DmiFRUVhw8cOLDEZBlzbZ2dnZdVtWmt02HMWljJ1juLvjFGi54j2t7erh0dHSuYLFPuROT8WqfBmLWy1Hb6/VP9uvjjgXnWWTc3URljjPEsNegX3/TyGeCb86zzPPCgiNT5F3Af9JcZY4xZI4tpsvkV4MfAfhHpFpHPAr8LPCAip4AH/HlEpF1E/hRAVa8A/xbvIdPHgc/7y4wxxqyRdXdzltXpm5UmIp2q2r7W6TBmLdyIfe8YY4xZIgv6xhhTRizoG2NMGbGgb4wxZcSCvjHGlBEL+sYYU0Ys6BtjTBmxoG+MMWXEgr4xxpQRC/rGGFNGLOgbY0wZsaBvjDFlxIK+McaUEQv6xhhTRizoG2NMGbGgb4wxZWTJQV9E9ovIyaJhTER+fdY694rIaNE6/6b0JBtjjFmq0FLfqKrvAncAiEgQ6AGemmfVH6rqx5a6HWOMMctnuap37gfOqOr5Zfo8Y4wxK2C5gv5jwFcWeO1uEXlVRJ4VkVvmW0FEHheRDhHpGBwcXKYkGWOMma3koC8iEeCngb+d5+UTwHZVPQh8AfjGfJ+hqk+qaruqtjc1NZWaJGOMMQtYjpL+I8AJVe2f/YKqjqnqhD/9DBAWkcZl2KYxxpglWI6g//MsULUjIq0iIv70UX97Q8uwTWOMMUuw5NY7ACKSAB4Afrlo2a8AqOoTwCeAXxWRPJACHlNVLWWbxhhjlq6koK+qSaBh1rIniqa/CHyxlG0YY4xZPnZHrjHGlBEL+sYYU0Ys6BtjTBmxoG+MMWXEgr4xxpQRC/rGGFNGLOgbY0wZsaBvjDFlxIK+McaUEQv6xhhTRizoG2NMGbGgb4wxZcSCvjHGlBEL+sYYU0Ys6BtjTBkpKeiLSJeIvC4iJ0WkY57XRUT+SEROi8hrInKolO0ZY4wpTUkPUfHdp6qXF3jtEWCvP9wF/Ik/NsYYswZWunrn48BfqOcnQK2ItK3wNo0xxiyg1KCvwAsi0ikij8/z+mbgYtF8t79sBhF5XEQ6RKRjcHCwxCQZY4xZSKlB/x5VPYRXjfNrIvKhWa/LPO+Z82B0VX1SVdtVtb2pqanEJBljjFlISUFfVXv98QDwFHB01irdwNai+S1AbynbNMYYs3RLDvoiUiEiVVPTwIPAG7NWexr4tN+K533AqKr2LTm1xhhjSlJK650W4CkRmfqcv1LV50TkVwBU9QngGeBR4DSQBH5pqRs7PTBO1+Uk8UiQWDhILBwgHg4SjwSJh71l0VAAPz1mA3NdJZN3SeUcb8g6pHPeMDV/+5ZaWmtia51UY9adJQd9VT0LHJxn+RNF0wr82lK3UeybJ3v5wndOX3UdEYiFvAxhKhOIhoJEwwFi/jgaChAtvBYgEgwQCXlDeGo66L/mL/MGIRQIEA4FCAeEUDBAKCiEA/44GCAUEIIBIRQUguJPBwIEZ8wLIqzLzElVcRXyrovremPHVfKu4vpjpzB2yTlKzvHGeccl73rz+anlrpLLu2Qdl5zjks27ZPLeeGo+OzXOu6TzDpmct04655DJu2TyDumcN87kXVJZb3wtX/j5O/knBzetwl4z5sayHO30V8Wn7t7OAze3eKU6/88/VbIrjLNOofSXzbszgkUm5zKRyTM0MR1IigNOzg9aq0UEgiIExMsEAn6mMDUd8DMGKVp/6rq4CDOWq59sZXoaFNXpq+ZTAd11FXdqWr11HNXC9GoJ+xllIZMN+xl0aDrDroqF5mTeU2d1sXCQeDhQOPObWj515re1PrF6X8aYG8gNE/Sbq2I0V63s6brjl1SLM4Ns3iW/QKk265dq8/57XFXyznRpeL75nOPi6lQQng6+XjAuCsSuon7Inhm8geLl6mcChRMHKUwLUxmEt2AqUwmKEJiVwQRFED8TCgWnz0qCRUMoUPx6gIh/9jN1phOec/YjRIJBwiEpnFGFg16QDwTW35mOMeXghgn6q8ELbl6J0RhjNiLrcM0YY8qIBX1jjCkjFvSNMaaMWNA3xpgyYkHfGGPKiAV9Y4wpIxb0jTGmjNww7fQHL44zOpCirjVBbXOCYNjyKzNTNp1npD/J8KUkm/bWUlVvfe8YM9sNE/RPvdzLK9/uAby7TKsb49S1JqhrraC2NUF9WwW1LQliFeE1TqlZSapKcizLyKUkw5cmGS4aTwxnCuvd/+l9HHj/ljVMqTHr0w0T9I/UPMXehq8w7Gxj2N3BcGobI++1cfGNJhydDvTxSIraiklqKtNUVeaornKorlGqawNUVIeQSBxCMQj741AUgpHp6VAUglEI+cuCkeI+Dsz1cl1wspBPT4/zWXAy09P5NOQzkE9BLo2bTTEx4jA26jI2GmBsPMjYRISxiSgjk5Vk8tHCx4cDGeoi/WwK9VJXc4E6OUddqJuayn+H9/gGY0yxGyboXzybYfJUPaGKDI0VZ9gS76KiNkAwoqQCCcadakbzNYyk6xlJ1dM9Vs+Es4niyxZBslQFB6gOdlMd7KcqOEhF8AoVgStUBIeoCFwhEkjP3XgwAoEwBMPedLBounh5IATBkDcuDMGi6bA3L4HpsUyNAxAIzFomFLpWm5qeMcZ/XYt6WpuanmeZOqAuuI4/785a5s+7Drh5b3ByM+fdvD+f84K4k/MHf3r2cjc3Y1eqQkYrmHQamHTr/HE9404zo04L404z485WlOmuMASHqtAVqsKX2Vv1LrWxYapDo1QHRok4aXJpIZl2SY+7ZJIu3RNNZHonabu91KPOmI3nhgn6r470UtF7ieokVCUhC4zNWicmUFMRJl4ZobEyRq6ykmxVC9l4C/loPflAHXmquZzfR2/mdvKpuV8/HHapTOSpiGepiGWoiKdJRNLEwinioRTxUJJYcJJ4cIKQpIsCXtYLhk4OcqmZAdLJzQyW6s4KtO48y5zp4F2SokxDgkWZzVTmInMzoHkzrvDM+akzpGAYDUTIESftVJDKV5LKx0nn4qTyMZKZBJPpGJPpCBPJCMlkkLwz93pMOJIjEkkTDk3QqKcJOVcIZYaIJAeIjvcTmkgSnkgTHc8SSeVJA7Oz53wAxuMwWgFXxgdpK3HPGbMRia5mf7qL0N7erh0dHXOW9030cX78PBPZCcZTo6SHB8levkzuymV0eASGxwiOThAeTRIZTxNJ5ogn88RTLpUpSGTnbisfjJKJ1JCJ1pKK1TBRUUsqVksmVksuUks+VIMbqPaC4Xwkj4SzSCSHxBwCUSUQhUBcCMSChCqCBBMhwvEgkViISDxEJBokHA0SDAeIRMKEIgFCwRBBCRKQAEEJEgwECUrQ6x3T7yVTRAjozK6WRQvlfFRA1ZtzUX/e76lTwcXFUQfX9caOOrjqksvlyOdcctk8+bxDLu2QTTtk03lyKYdcKo+TdHCSeZyUi5sGNwNuWtBMEM2G0XwEdIF9pDmCzhih/AjhzAjRzAiJ1CgVyRES6REimVGi2VGCs84IACZiMBmDZCJIOhEimwiTqYqRr63AramEuloCdbWEGhqINjQRr22gKlZDVaSK3bW7qY/Vz/+ziXSqavv8CTZmY1tySV9EtgJ/AbQCLvCkqv7hrHXuBb4JnPMXfV1VP7+U7bVVttFWef1lN1dd0vk0k6lRJq8MkLzST/rKZbLDQ+RGR3AmJnAnJmBiksrJS1RNnkWSaYKXMwRTWUKpHOTCQAIlQS5SSS5cSTbsjXPhCnLhKrLhSvKhGE4oTiYYww1G5kmN4w+zaB5RrypEyIHmQPN4u9XFC+suWjQNWjQvCAG8LCAABPzulGcuQ8IgYXTGeKEeRf33MH29JOhkCOZThPNpQvkkkdwE4dwk4dxEYYj4Y0iCpNBIHicWJl8RxY1HcevjUBGHigpSlTvJVlaSrKwkXFtHtK6ReH0TibpmEvVN7IhVEQvG1uUDZ4y5UZVSvZMHfkNVT/jPyu0UkRdV9a1Z6/1QVT9WwnZKEpAAiXCCRDhBU3Ub7Fj6Z6nj4KbSaCqJm0rhJpO4yRRuKokzOUk+nSSfSZFPXyGXTJFO5cim8l7JOeOQy0I+B47jDa4ruG4AVwM4bgDVAA5BVAO4BFECqF+296ZDIOIvm3rNK90LLqgiXrkf8ZeJnzEISgCHgDgExPXHSiDgEhCXYFAJBCAUVMIRIRQRItEA0ViQSDxMpCJCJBYnGIsTjlUSijcTTFQQiMeReIJARYJAPO7Nx2JIwJrUGrMelfK4xD6gz58eF5G3gc3A7KC/YUgwSLCyAior1jopxhizJMtSHBORHcCdwMvzvHy3iLwqIs+KyC0LvP9xEekQkY7BwcHlSJIxxph5lBz0RaQS+Brw66o6u0HNCWC7qh4EvgB8Y77PUNUnVbVdVdubmppKTZIxxpgFlBT0RSSMF/D/UlW/Pvt1VR1T1Ql/+hkgLCKNpWzTGGPM0i056IvXpOLPgLdV9fcXWKfVXw8ROepvb2ip2zTGGFOaUlrv3AN8CnhdRE76y34L2Aagqk8AnwB+VUTyQAp4TNfbjQHGGFNGSmm98yMKt3suuM4XgS8udRvGGGOWlzWmNsaYMmJB3xhjyogFfWOMKSMW9I0xpoxY0DfGmDJiQd8YY8qIBX1jjCkjFvSNMaaMWNA3xpgyYkHfGGPKiAV9Y4wpIxb0jTGmjFjQN8aYMmJB3xhjyogFfWOMKSMW9I0xpoyU+ozch0XkXRE5LSK/Oc/rURH5G//1l0VkRynbM8YYU5pSnpEbBP4YeAS4Gfh5Ebl51mqfBYZVdQ/wB8DvLXV7xhhjSldKSf8ocFpVz6pqFvhr4OOz1vk48GV/+qvA/VMPSjfGGLP6Snkw+mbgYtF8N3DXQuuoal5ERoEG4HLxSiLyOPC4PzshIu+WkK611Mis73YD2gjfAa7+PbavZkKMWU9KCfrzldh1Ceugqk8CT5aQlnVBRDpUtX2t01GKjfAdYON8D2OWWynVO93A1qL5LUDvQuuISAioAa6UsE1jjDElKCXoHwf2ishOEYkAjwFPz1rnaeAz/vQngO+o6pySvjHGmNWx5Oodv47+c8DzQBD4c1V9U0Q+D3So6tPAnwH/RURO45XwH1uORK9jN3wVFRvjO8DG+R7GLCuxgrcxxpQPuyPXGGPKiAV9Y4wpIxb0l4GIdInI6yJyUkQ61jo9iyUify4iAyLyRtGyehF5UURO+eO6tUzjtSzwHf4vEenxf4+TIvLoWqbRmPXEgv7yuU9V77jB2oZ/CXh41rLfBF5S1b3AS/78evYl5n4HgD/wf487VPWZVU6TMeuWBf0ypqo/YO59E8VdZ3wZ+JlVTdR1WuA7GGMWYEF/eSjwgoh0+l1K3MhaVLUPwB83r3F6lupzIvKaX/2zrquojFlNFvSXxz2qegivx9FfE5EPrXWCytyfALuBO4A+4D+sbXKMWT8s6C8DVe31xwPAU3g9kN6o+kWkDcAfD6xxeq6bqvarqqOqLvD/cmP/HsYsKwv6JRKRChGpmpoGHgTeuPq71rXirjM+A3xzDdOyJFOZlu+fcmP/HsYsK7sjt0QisguvdA9etxZ/paq/s4ZJWjQR+QpwL143xP3AbwPfAP4bsA24APycqq7bC6ULfId78ap2FOgCfnnqOoUx5c6CvjHGlBGr3jHGmDJiQd8YY8qIBX1jjCkjFvSNMaaMWNA3xpgyYkHfGGPKiAV9Y4wpI/8/30JycrltNokAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x1a23f27320>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "def simulate_norm(depth):\n",
    "    for width in [32, 64, 128, 256, 1024]:\n",
    "        factor = width // depth\n",
    "        xval = norm.cdf(\n",
    "            numpy.linspace(\n",
    "                norm.ppf(1/width),\n",
    "                norm.ppf(1 - 1/width),\n",
    "                width))\n",
    "        \n",
    "        ty = []\n",
    "        tx = range(depth - 1)\n",
    "        for i in tx:\n",
    "            r1 = xval[i*factor : (i+1)*factor]\n",
    "            r2 = xval[(i+1)*factor : (i+2)*factor]\n",
    "\n",
    "            u1 = numpy.mean(r1)\n",
    "            u2 = numpy.mean(r2)\n",
    "            o1 = numpy.std(r1)\n",
    "            o2 = numpy.std(r2)\n",
    "\n",
    "            d = abs(u1 - u2) / (numpy.sqrt(4 * o1 * o2))\n",
    "            dd = numpy.sqrt(d*d - 1)\n",
    "            t = 2 * numpy.log(d - dd) + 2 * d * dd\n",
    "            ty.append(t)\n",
    "\n",
    "        matplotlib.pyplot.plot([x+1 for x in tx], ty)\n",
    "\n",
    "ax = matplotlib.pyplot.subplot(2, 2, 1)\n",
    "ax.set_ylim(0, 10)\n",
    "simulate_norm(4)\n",
    "ax = matplotlib.pyplot.subplot(2, 2, 2)\n",
    "ax.set_ylim(0, 10)\n",
    "simulate_norm(8)\n",
    "ax = matplotlib.pyplot.subplot(2, 2, 3)\n",
    "ax.set_ylim(0, 10)\n",
    "simulate_norm(16)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 64,
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAVsAAADuCAYAAACXv6SfAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvNQv5yAAAIABJREFUeJztnXl4VOXZ/7/PzJmZbBDCEpZEdgthjYTVUkT9RSJavNBUQahU4FJ74XtpbWuX69W+tn3Vtm68lldbtWq1lf70pwVFUzYlFAREdoEgSgSSsGRfZj1nnt8fwzlMklnOJDPnOZO5P9d1LkjynPM8Z5L5nnvu514Y5xwEQRBEYrGIXgBBEEQqQGJLEARhACS2BEEQBkBiSxAEYQAktgRBEAZAYksQBGEAJLYEQRAGQGJLEARhACS2BEEQBiDFOJ7SzQiC0AsTvQAzQZYtQRCEAZDYEgRBGACJLUEQhAGQ2BIEQRgAiS1BEIQBkNgSBEEYAIktQRCEAZDYEgRBGACJLUEQhAGQ2BIEQRgAiS1BEIQBkNgSBEEYAIktQRCEAZDYEgRBGACJLUEQhAGQ2BIEQRgAiS1BEIQBkNgSBEEYAIktQRCEAZDYEgRBGACJLUEQhAGQ2BIEQRhArK3MUxLOOTweT9RxFy5cQL9+/WC1Wg1YVWeam5thtVqRmZkpZH4zcO7cOQwaNEjY/PX19cjMzITD4Qg7xu/3w2KxIC0tzcCVEaIhy1YHnHNwzqOOO3PmDGRZNmBFoamvr0dTU5Ow+YHAayCSU6dOCZ2/uro66oOZcw6v12vQigizwPSISBAxDe4p+P1+VFRUoLm5OeK4hoYG9O7dW5hl63Q6wRhDenq6kPkBoLa2Fv3790/Z+RsbG5GVlQVJCv+hkXOOfv36YfTo0WCMGbg6w+nRNxcr5EbQyaBBgzBw4MCIY44dO4Zhw4YJE7vq6mowxjB48GAh8wNAS0sLrrzyypSd/+jRoxgxYkREF0FrayvOnDmDkSNHCnswE8ZDYqsDxhjS0tJgsUT2ujgcDjgcDmE+07S0NHDOhfpsRfuMRc9vsViQlZUV0WcrSRIOHToERVFIbFMI8tnGEYvFAkVRhM3PGNPlWyYSh7r5FQmLxUK/pxSExDaOWCwW+P1+ofPTm1gsTdeuwLacIl1jRW6mEsZDYhtHRIstABLbJIBzDsYYnE6n6KUQBkJiqwN1xziakIkWW7JsxbI5a5LusRaLhcQ2xSCx1YmeEB2r1SpUbFPdZ5tM905im3qQ2MYRM2yQiXZjiET9eC6CWKxadZ0ulyupHhBE9yCxjSOi3Qhk2fKokQBmQI1YcLvdKf37SjXM/5eZRJDYikWUZRuLVQsExFaSJBLbFIOSGnQiy3JUH5vb7YbP5xNWn8DpdMLtdgutj6AoirD5ZVkWOn8wkdbgdDphsVggyzL8fj8lNqQIJLY6cbvdOH36dETLyeVywefzCfPbejweuFwuocVgvF6vsPkVRYHb7TZ0/nMz7wz5/UhrqKurg8PhAGMMsizDZrMlanmEiSCx1UmvXr0wYcKEiGJ78eJFNDQ04Fvf+paBK7tMQ0MDampqMG7cOCHzA4EyjxMmTBAyt9vtxhdffGHo/OfCfD/SGg4dOoS2tjbY7XY0NTUJLRxEGAf5bHWixxco2meb6nG2RvtsY/XVqthsNu3v5JtvvonnkggTQ2IbA2ZPagCSK9Y03ogM/YoFSZLg9/ths9koZTeFILHVSbJYtqLFXiRGim1XrVrgstharVYoipLSD8hUgsQ2jpghgyyVSRbLVk1+kSSJxDaFILGNI5RBJhajxLY7Vm0wFosFPp+PxDZFILGNI6I/xlNSg/kzyDjnUBRF28zknKf0AzKVMPdfpolIBp8tiW3iLdt4WLXq5pjX69VibYmeD4ltHCGxFYvf708Kn63f74fdbofX64XNZoPb7Ra9JMIASGx1Qpat+Um0ZRsvX60qtrIsw+FwoK2tLS7XJcwNtTLXiaIo2LVrF7xeb8Rxra2tyMrKMmhV7eGcw+l0Cm14KPL+ZVmGLMsRO9t2B+cN9+kem7HxxZDfV322ubm5uHjxIvr27YvMzEyMGTMmXss0E+b/mGEglK4bA1OmTIm6AbNz505cffXVBq2oPbIsY9++fZg+fbqQ+QGx919bW4u6urqECFesVm2418Dv9+Pjjz+Gw+HQHgxqXdtkcIEQXYfcCD0IciMkj2ClpaVBURRkZGTA4/GIXg5hACS2OkmGNzHF2SZGbOPlqwXah6dxzpGeno76+nr4fL64zUGYExLbHkQyPBASSbJYtowxTVwlSYIsyyn9iSRVILHtYaTymzYRYhtPq1ZFzRwLjrFN5d9bqkAbZAlAlIWVDFZdIkmWDDLVsmWMwePxgDEGt9udsCiKZKXIksmbefzT30/C8y/OeUncLxwFElud6BUyNdaWWp0YT7wfcomwalU8Ho/W9NFisaCtrQ19+vRJ2HzJSDMUrE4fHvfr3uSq6B/3i+qAxDbOkNiKI1kyyABoqbputxuMMUpsCAGzMth6JUCiXPG/pB5IbOOMWqOU+koZTzwt20RZtapvNlhsg9uaJ8vDwgiYhcGabm63UCyQ2MYZ0Sm7qUwyiJW6RtW/rIqtyyXI3DIzFpDYpirnz5/HhQsXIo5pamrCsWPHIEliXlqn04mDBw8KmVv0/G1tbWCMob6+vlvXufjt73d7LeFeA7/fD0VRYLfboSgKvF4vJEmCx+Mh91MHGGOQ0nrO60FiqxPGGHJycpCVlRXRejpx4gSGDBkirD5Ac3MzRo0aJWRu0fNXVVVBkiQMHDiwW9e5GIe1hHsNzp07h4aGBjgcDi1NVw0FI7FtD2OA1d5zXg8S2xhwOBxwOBwRxdbhcCAtLU2Y2FqtVmFzi57fbrfD4XB0a/54+WrDraF3797gnMNut8PlckGSJFitVsiyTO6njjAGq43cCEQYRLfGSWWSwWdrt9u1Eouq6KpuBbfbDYfDIXqJpoFZAKvDeIlijKUBKAfgQEAj3+Gc/6rDmDkAngMwCcAizvk70a5LYhtnRDd9TGW6K7bxsmq/84ebwv5MkiRwztu1M/f5fMjIyEBrayuys7PjsoaeABNn2XoAXMc5b2WM2QD8mzH2Eed8V9CY0wB+AOAnei9KYhtnKBpBHGbIIIsktCqMMdhsNnDOtX+tViu++uor5OXlGbDKJIExSAIsWx6Iz2u99KXt0sE7jKkEAMaY7jc7iW0M6ClhSGIrju5YtonMFgtGURQwxjR3k81m0/6lXmTtYQywSAnZIOvPGNsb9PWfOed/bj83swL4HMBoAGs457u7OymJbZwxg9gmg+8yEYjOINNj1SqKAqvVqlX9Ut1OkiRBUZSU/d2FhDFYpIR8UqnlnE+NNIBzrgAoZIz1AfAeY2wC5/xIdyYlsY0zojfIVOs7Fd+wXb1vo6xa4LLYqgXD1b8X1YdL4V+XYQxC3AjBcM4bGWOfACgBQGJrFHrcCMFWiwhSuVuDyIeMHqsWgGbFejweSJKk9bRTPxGl6u8uJIwlyo0QZVo2AIDvktCmA/g/AH7X3euS2MZItDe0aDcCiW1sYhsPq1av0ALQ/LNer1drZ64mNag/F5V9aDYCPlshG56DAbx+yW9rAfB/OecfMMZ+DWAv53w9Y2wagPcA5AD4LmPsMc75+EgXpd9qDMiyjBMnTkR8QzudTni93qhdeBOFy+XCiRMnhO3Ku91uHD9+XMjcDQ0N8Pl8OH/+vJD5gwn3GjidTthsNjidTjgcDng8Hk10JUmiWNtgmAVWu/EFnTjnhwBcFeL7jwb9/zMA+bFcl8Q2BiwWCwYMGBBRbJuamtDS0oIBAwYYuLLL1NfXo1+/fsKsowsXLgi799bWVuTk5Ohu5b5/2LXdnjOcVRvuNTh8+DD69+8Pn8+HXr16wel0wm63w+PxwGazoa2tjWJtL8EYwKyUQZaSWK1W9O3bN6rV6PF40K9fP4NW1R6Hw4Hs7GxhVf8lSRJ279XV1ejTpw969+5tyHyR3AfhXgM1c0xRFDgcDjQ1NSEzM1OzaKn6VxCCfLaJgsQ2zlitVqE+01SMQlCJxWdrZASCipq8wDkH5xwOhwOyLCMtLQ3Nzc3Izs6G0+k0fF1mJVCIpudIVM+5EwPQ80Y2Q+iX6DhfURgZjRDLplgwalytxWLRYmuDRbetrS1lQ/c6w8B6UBgciW0M6BVbikYQg9503e5atV0VWuByJw/1ocw5R3p6uia2dXV1FGurQm4EIhIktuIwu0WoPgyCw70YY9qnEYfDAa/Xm7K/v44wxmDpQe2lSGzjDImtOIxI1+2OVQugXRtzWZa1f4HLLqimpiZhm4ymggGsB1m2PSeuwiSILrGYymKrx7LtjgshFqH1vPCLkN9XIxEYY1rTR5/PB4vFohUQb25u7vIaexQs4LON9yEKsmxjIFk2yEhszYkahaD+jtSsMY/Ho4mu6HRvMxFwI/Qcieo5d2ISRLsRLBYLiW0YjLJqI6HWRvD7/fB4PJoPV7V01a+JAD3JjUBiGyNnzpxBW1tbxDFutxtHjx41aEXtaWhogMfjEZbyKfLe29raoqZTG0nH14FzDpfLpYV4AQHrze12A7hs4b722mt46qmnDF+v6WAMTKINspQlJycnaoZSbW0tBg8ebNCK2uPxeJCTkyMs5VPkvdfX12PQoEEhw78+v+KaLl+3q1Ztx9fB5XLh3LlzyM7ORmtrqxaJ4Ha7YbFY4Ha7tXAwApfEVkjVLz09yBwA/gqgCEAdgDvU7g3hILGNkaysrKixnJIkIScnx6AVtefChQvo1auXsPlF3rvVakVOTk5ci/B0x33Q8XXIyMiA3+9HWlqa1vRRFViLxQKPx4OzZ8/igw8+wNNPP93dpSc/DICYDS09PchWAGjgnI9mjC1CoATjHZEuSmIbA2b5eBoJyiDr/DsSkZobCjVVNz09HZxzpKWlweVyaX5+n8+HzMxMU1QtMwWC3Ah6epABuAXAf136/zsA/sgYYzzCxxIS2x5GKkcjAPF9IMZrU6wjwe3M/X4/MjIytAI0Ho8HRUVFCZk32WCJS9eNRw+yPABnAIBzLjPGmgD0A1AbblIS2x5GqottR7pq1SZCaFXLW40YsdvtAAKV2tra2pCRkQGv14srrrgi7nMnJYwlyo0Qjx5koZ7qEd94JLYJQi02YjQktuYlOF1XFVu1+hcApKWlQVEUZGRkCF6pSWAMEByNEKEH2VkAVwA4yxiTAGQDqI90Lcog6wJ6+pCJ8puS2F7GTFatCmNMC/FSfbjBFi7nnMQ2GKs1/kcUGGMDLlm0COpB1rH1xnoAyy79vxTA1kj+WoAs25hQQ3WiITKxgcS2eyRSaAG0C/UCAp+AbDabtmEGgMRWJXFuhGhE7UEG4BUAbzDGTiJg0S6KdlES2wRAYises0QgBKMWynG73e3amdtsNi0kzGKxkNiqCHIj6OxB5gbwvViuS26EBCBSbFM5Xbe7JNqqVTs1eL1eTWxVC1d1J0iSFFZs3W43pk+fjsmTJ2P8+PH41a9+FXIcALzzzjtgjGHv3sub7ocOHcKsWbMwfvx4TJw4UctcMzUWa/wPQZBlGyP19fU4depUxDEtLS04fPiwkALQ6hvo3Llzhs8NBLrHBr/BRczdOHd5zOcmQmg7vg5qtS+1uaPH49HEFwgkhNhstrBi63A4sHXrVmRlZcHn82H27Nm48cYbMXPmzHbjWlpa8D//8z+YMWOG9j1ZlrF06VK88cYbmDx5Murq6mAze61YcW6EhEBiGyPZ2dmYMGFCRN/t8ePHMWTIEMMaDwZTU1MDv9+PvLw8w+cGgM8++wwTJ04UOvd2IbN3puPrUFFRAZ/PpxWj8Xq9sNlsWmKDz+cD5zxsd2DGGLKysgAAPp9Pq4vbkUceeQQPP/xwu/oKGzduxKRJkzB58mQA4RtSmgkOBm41+QMhBkhsY0SSJFit1ohia7PZIEmSkGIwNpsNPp9PWCEai8UidO7t/abFfF6i3AcdXwc12gCA1rbcZrPB6/VCkiS4XC54vd6ID2lFUVBUVISTJ09i1apV7axXANi/fz/OnDmDm2++uZ3YqgV65s2bh4sXL2LRokV4+OGH43i3CYAxcIEf++MNiW0CoA2y1IbNviHk99WoA/WBrabnqm4Fl8uFt956C6WlpWGvbbVaceDAATQ2NmLhwoU4cuQIJkyYACCwAfejH/0Ir732WqfzZFnGv//9b3z22WfIyMjA9ddfj6KiIlx//fVxuefEwMClniNRtEEWIxT6ZV6cN9wX8znxtmrDCS1wuTaCdElA1JRdr9cLu90Oj8eD8+fPo7GxMeo8ffr0wdy5c1FWVqZ9r6WlBUeOHMHcuXMxfPhw7Nq1CwsWLMDevXuRn5+Pa665Bv3790dGRgbmz5+Pffv2df+GEwkDuEWK+yEKEtsEQEkNxvP2HnvM5yQ6+qAj6t+F2mGXc665EdQKYH379sXIkSNDnn/x4kVNiF0uFzZv3oyxY8dqP8/OzkZtbS0qKytRWVmJmTNnYv369Zg6dSrmzZuHQ4cOwel0QpZlbNu2DePGjTPkvrtOwI0Q70MUPcdGNwhVzCJZuCJb46Sa2HZFZBNFJKtWRRVb9V9VeB0OB1pbW5GRkYHCwsKQ59bU1GDZsmVQFAV+vx+33347br75Zjz66KOYOnUqFixYEHbenJwcPPTQQ5g2bRoYY5g/fz5uusnYh02scMbArT1HonrOnZgIciMYQ3eE1kj3gYr6kA5u8uj3+7X6CF6vFx6PJ2zo16RJk7B///5O3//1r38dcvwnn3zS7uulS5di6dKl0W/GRHBGG2REBEhsE0t3rVkRQgtAKxIOBDasVKtWrZOgKEpEsU09mFAfa7zpOXdiItSYSVFz92SxNZPbANAvtJxzrQ6C3++HLMuQJEnrqCvLMgCQ2AbBGYOf3AhEJILz3o2mJ3dqiIfQxtOq1Su0KoqiaGKr1kmQZVn712q1Ii0tTVicsvlg8JMbIXVhjGH37t0RrUefzwdFUYSkzMqyDJ/Ph4aGBsPnBoDW1lbs3LkzrteskubG5TpGC23H10FtxulyubQNVrV9uc/ngyRJyM7ODrn56na7MWfOHHg8HsiyjNLSUjz22GPtxrz44otYs2YNrFYrsrKy8Oc//xnjxo3D3/72N/zhD3/Qxh06dAj79u0LuxFnGhiERA8wxq5AoJnjIAB+BDo5rO4wJgfAXwCMAuAGsLxDcfFOkNjGCGMMM2bMiFgYvLa2FnV1dRgzZoyBKwvQ2NiI6upqYWE9O3fuxNVXXx2365nNbRALwa8D5xwff/wxbDYbWlpatAQHNeVWlmVYLBb06dMn5LX01EW48847cd99gVjj9evX46GHHkJZWRmWLFmCJUuWAAAOHz6MW265xfxCi0C6rmIRkq4rA/gx53wfY6wXgM8ZY5s458G96X8J4ADnfCFjbCyANQAiZoiQ2CYA0RtkPcGNEG+RFek+AC5HIkiSBEVRkJmZCafTCQBaZ91IYqunLkJwmm9bW1tIC/mtt97C4sWLY16/GMS4ETjnNQBqLv2/hTF2DIGeY8FiOw7AE5fGHGeMDWeMDeSch+3WSWKbAESLbbJvkPU0oVWxWCxaDHZaWhqam5shSRI451p9hOzs7LDnR6uLAABr1qzBM888A6/Xi61bt3b6+T/+8Q+sW7euy/dgJByAnyUk7ypqw0cVxthwBGrbdmz4eBDArQi0OZ8OYBiAfABhxZYyyBIAZZB1jbf32E3tNuiO0AJo151BLUrjcDigKAoURUFDQwO2bdsW9ny1LsLZs2exZ88eHDnS2UW4atUqfPXVV/jd736H3/72t+1+tnv3bmRkZGi1FEwPC7gR4n3gUsPHoCOc0GYB+H8AHuScN3f48ZMAchhjBwD8B4D9CLgfwkKWbYzoqY3AGKMMshhJlMjGy6rtrtCqzR7VNN309HRNbNU25larFYcPH456reC6COGEc9GiRfjhD3/Y7ntr165NIhdCwGcrKhqBMWZDQGj/xjl/t+PPL4nv3ZfGMgCnLh1hIcs2AZBlGxtmtmaB7gutdp1LG2EAYLfbtVY4akEan8+HQYMGhTw3Wl0EAPjyyy+1/2/YsAFXXnml9rXf78fbb7+NRYuitsoyFX5mifsRjUvi+QqAY5zzZ8KM6cMYU/9wVwIoD2H9toMs2wRAPlt9JFpk42HVxktoVdQ0XQCaZas2e3S73Rg/fnzI8/TURfjjH/+IzZs3w2azIScnB6+//rp2fnl5OfLz88MWuTEjHAwKFyJR3wbwfQCHL7kJgED0wVAA4Jy/CKAAwF8ZYwoCG2crol2UxDYBUA+y6CRCaBu2HkfOdQFrz+iKXtFQfydqFIHqTlAtW4fDAZ/P1626CKtXr+70c5W5c+di165d3bwL4/EL+PDNOf83gIj+Qs75pwCujDSmIyS2CYDcCOExu8sgmHhbtWpcrcVi0XrFBbcx9/v9YVvipCIBy5YyyFIWxhiam5sjdiblnMPj8QjLIHO5XMIaPvp8vrBzbz891JA1mMV9EPw6qB//1apfbrcbFotFeziqKbrp6endnrcnIcKyTRQktl3A6/Wira0t4hi/3x91TCJQ8+5FzK3OH2rufXUFhsxvFqEF0O51aG1thdfrRVpaWjuxVTvrqt8ny/YyAcuWxDal6d+/P3JzcyOOOX/+PEaNGmXQii4jyzLq6+uFzA10vu9kchsA8XUdXOE8DfvEawEA9fX1uHjxohaF4PF4IElSO3eC1WoN6bPVUxfhmWeewcsvvwxJkjBgwAD85S9/wbBhwwAADz/8MDZs2AC/34/i4mKsXr1aVwijcDh6lBuh5zw2DMasflEz+WyNFNoL9YGjO8TbRxuM3W7X2uAwxuD1emG1WuFyudqVWAwltmpdhIMHD+LAgQMoKyvrtNl11VVXYe/evTh06BBKS0u1zrk7d+7Ejh07cOjQIRw5cgSfffZZxMQJM8HB4OeWuB+iIMs2RsxuEZhBbI0W2XiQSKEFLm+E2e12LRrB4XBoFq7L5YKiKCHFVk9dhGuvvVb7/8yZM/Hmm29q57rdbni9Xm2DbuDAgQm80/iicHO/32KBLNsuYGbBFS228SqHGI14WLJGov7NqLUQgIAAB4ut2+3WRLUjiqKgsLAQubm5KC4uDlkXQeWVV17BjTfeCACYNWsWrr32WgwePBiDBw/GvHnzUFBgjP+8u3AwKH5r3A9RkGXbwzD6QbD4p2fafX3r9xLvK463yCbaqgWghQKq7cytVqu2QaZ21v3nP/+Jb3/72yHPV+siNDY2YuHChThy5EjIVN0333wTe/fu1VwFJ0+exLFjx3D27FkAQHFxMcrLyzFnzpwE3Wkc4T3LsiWxTRBqqcNIdW+TkY7iaiSJsGSNEFoA2t+C2h7H4XBo6btZWVlwu93tUnLDEakuwubNm/Hf//3f2LZtmxZK9t5772HmzJmaxXzjjTdi165dSSG2HBDqY403JLYJQs0iS2axFSmswSTKXWCU0AIBN4Da1DG4nTkQcC20traib9++GDFiRKdzL168CJvNhj59+mh1EX72s5+1G7N//37ce++9KCsraxcpM3ToULz00kv4xS9+Ac45tm3bhgcffDCxNxsnOBhkf/K+fzpCYtsF9PhFRabsdhWziKtKIn2yRghtY/ZQqLKniq0sy5obQe22a7PZ4PV6kZmZiWnTpnW6jp66CD/96U/R2tqK733vewACIrt+/XqUlpZi69atmDhxIhhjKCkpwXe/+92E33u88JMbgYiGyJRdPZhNWINJ9MaXUUIbjGrNql2XVX+tJElaIkq4zrp66iJs3rw55DqsViv+9Kc/dedWhME5oPhJbFOa5ubmkIWbg3E6ndi7d68QN0KopovPv3eF4evoCskutB1FVv09qJ111fha9XvqQ5kxFrEQTaoiCxBbnQ0fswG8iUAlMAnAU5zzVyNdl8S2C2RnZ2PWrFkRd/4PHTqEESNGoFevXgauLMDOnTuTRlxVkimMKxwdhRa43PTxs88+g6Io8Pl8WsytLMtIT0+HoiiQJAk2m43ENgjOmSjLVk/Dx1UAjnLOv8sYGwCggjH2N865N9xFSWy7wOOPP465c+dG7FDq9/vhcrkgSd1/iW/+/sGYxueOSB6hNVJkE2nVhhJaAFoXBjUSwe12IzMzEz6fT7NoVQvX4XCQ2AbBIcay1dnwkQPodanQeBaAelBbnPjzxRdfYMGCBRGFVFEUpKentxtTsvhzI5aXFBhtyYoQWgDa71/tqqt2ZVBrIqjuA0mSkJ6eHrLql57aCOXl5XjwwQdx6NAhrF27FqWlpdrPkrY2AoAEbXvEo+HjHwGsB1ANoBeAOzjnEVdLYtsFqqurMWTIkIj+2Ief8AH4wrhFJRGpIrTA5SaPavQBEMgca2lpgc1mAxBIwU1PT0evXr00f24wam2ErKws+Hw+zJ49GzfeeCNmzpypjRk6dChee+01PPXUU+3ODa6NAACzZ8/Gtm3bMHfu3C7fs1EkcIOslnM+NdqgKA0f5wE4AOA6AKMAbGKMbY/UGofENkbUttNpaWnazjKhDxF+2UQJbTSRVVH/RtQkF8YYbDab1s5clmUtBKxPnz4hr6GnNsLw4cMBoJMBkMy1ETgAWRFjgUdr+IhAs8cneSAG9CRj7BSAsQD2hLsmiW0MlJSUoLa2FpWVlZg5cyays7M7WRJEZ0RtfokWWgCaRalW9bJYLFqXXYfDgZaWFgDQKoKFQ1EUFBUV4eTJk1i1alXE2gjBBNdG4Jzj/vvvT5raCOCAIiB6Uk/DRwCnAVwPYDtjbCCAMQC+jnRdEtsYKCsrQ0NDA2655RaUlZWFHRf4uLjPuIWZFJERBmYQWgAoKioCAHg8Huzbt0+zPFXLtrGxUfPn/uUvf2kXOxuM3toIHUnm2ggcCfPZRkNPw8ffAHiNMXYYgX5lP+Oc10a6KIltjFRVVYVtN62iVt9PVXpCGFcoYhXaYNRqX6o7AYDmRsjMzISiKKirq4t6nUi1EUKR1LUROOCLuL+fqHl1NXysBhDTE73nJB4bRHV1dVSxveUHkRMeEknuiDxhc/s8PlMIbSJhBmsAAAAa90lEQVSs2u4IrYrqZ/X5fFrfMb/fj7S0NLhcLs3v2pHgAjVqbYSxY8fqmnPo0KHYtm0bZFmGz+fDtm3bkseNgIBlG+9DFCS2MVJVVYXBgweLXoZwfB5fp0M0TU2yqYVW7dagdthV62ekpaXB4/GE/buqqanBtddei0mTJmHatGkoLi7WaiOsX78eQCBpIj8/H2+//TbuvfdejB8/HgBQWlqKUaNGYeLEiZg8eTImT56cNLURAtEIPO6HKMiNECM1NTUpJ7ZmENJINDUl5rNmPES2rqoS/fKGAwi0xnG5XFo0gd/v1yxctRBNKPTURpg2bZrmlw0m2WsjyALcCImCxDZGqqurMWXKFNHLSAhmF9WOJEpkgfgIbTAWi0Vr9uj1erU+ZEAgjjZcS5xUR0Q0QqIgsY2RmpqaqD7bZCDZhDWYRIosEH+hVbHbA73Z1DrHHo9Hi7vlnFMb8w5wDiiKOZqXxgMS2xhJNjdCMotqMIkWWJVECS0AbUNM3ShTxVb9Hlm27eEAZLnniC1tkMUA5xxOpzNsUz4zYLZNq+7S1CT3CKEFAmIb3M7c7XZrdW39fn9IsXW73Zg+fTomT56M8ePH41e/+lWnMR6PB3fccQdGjx6NGTNmoLKyEkAgBPHuu+/WNsc++eSThN5f3OHx3xwTuUGW0mJbUVGBwsJC7ejduzeee+65sOPV7gzJUsQjmTFSZIHECy1wWWztdrvmRrBYLHC73VovslDnbN26FQcPHsSBAwdQVlaGXbt2tRvzyiuvICcnBydPnsSPfvQjrWXOSy+9BAA4fPgwNm3ahB//+MemLmjfEQ7Ar/C4H6JIaTfCmDFjcOBAIEFEURTk5eVh4cKFYcc3NDSEzV8nuo+R4qpihMiqqKFfalquLMuQJAkejwcejwe9e/fudI6eugjr1q3Df/3XfwEIhHrdf//94Jzj6NGjuP766wEAubm56NOnD/bu3Yvp06cn8C7jRyCpgdwIPY4tW7Zg1KhRGDZsWNgxVVVVUYt4ROtNRnTGaCtWxUihBS5/IlLdCGovMrfbjWPHjqGtrS3keYqioLCwELm5uSguLu5UF6GqqgpXXBGoYSxJErKzs1FXV4fJkydj3bp1kGUZp06dwueff44zZ8zbDqkTnCzbHsnatWuxePHiiGNqampQXFzc6WNcMIGPaeELihCXESGwKkYJ7ZdnzuHLM+cABB7EFosFVqtVE1rVsj18+DDy8/NDXiNaXYRQD3jGGJYvX45jx45h6tSpGDZsGK6++uq4FLM3Cs45ZF/yuD2ikTyvfALxer1Yv349nnjiibBjSkpKcPz4cciyjIKCAqxbty7kuBvu+CxRy+wRJFpgnzp6PX4ybkvEMUYIrZrIoP4LBB7E+/fvB+e800YZ5xytra0RrxmuLkJ+fj7OnDmD/Px8yLKMpqYm9O3bF4wxPPvss9q4q6++GldeeWVc7zPRKAJ8zDp7kP0UwJJLX0oACgAM4JyHTVgnNwKAjz76CFOmTInoIigrK8PKlSvx+OOPhxVaIjyiXAUdMVJoQ2G327WODarY+v1+5OXlYenSpZ3G66mLsGDBArz++usAgHfeeQfXXXcdGGNwOp2aa2LTpk2QJAnjxo2L010mHjXONt6HDtQeZAUAZgJYxRhr98Jxzv/AOS/knBcC+AWAbZGEFiDLFgDw1ltvRXUhAIHssWnTphmwop6DGQQWMEZkvzxzrl33hI5YLBYtW0xtbS7LcsTOujU1NVi2bJl2zu23367VRZg6dSoWLFiAFStW4Pvf/z5Gjx6Nvn37Yu3atQCACxcuYN68ebBYLMjLy8Mbb7yRsHtPBAE3giJiXj09yIJZDOCtaNdNebF1Op3YtGmTrvzxc+fOJVVCg0jMIrJA4oVWtWRV32wk7Ha7VjBckiStKajVag0ptnrqIqSlpeHtt9/uNGb48OGoqKjQexvmI3EZZPHoQab+PANACYD7o02a8mKbkZGhq44okHzZY0aiFsU2k8gCxgktcLn1TaTedA6HA/X19WCMaZ11bTYbJEmiDLIOcHD4E1McIR49yFS+C2BHNBcCQGKrG845XC5Xyr8hIgmJ2Uik0Ibyy9psNvh8PjgcjrDnORwO+Hw+rR6CWoDGbren/N9WRziHEDcCoKsHmcoi6HAhACS2ukml7LFkEtRwGC20QMBF4PV6I4qt3W7XBFlRFCiKAqvVCofDQWLbEc6hCCj7pbMHGRhj2QCuAdB5ZzMEyf+uMoi6ujrk5OSIXkbcUItXhzqSmYwMa8KEtl/e8KiRBtFaIqnRCKroAoE42vT09JBie+bMGVx77bUoKCjA+PHjsXr16k5jGhoasHDhQkyaNAnTp0/HkSOXO4WsXr0aEyZMwPjx4yOmopuRQLquP+6HDtQeZNcxxg5cOuYzxu5jjN0XNG4hgI2c89DZKB0gy1YnetrhmI1kF069ZGRYE3r9SAIbjB6xVRMa7HY7nE6n5rvNyMgIKbaSJOHpp5/GlClT0NLSgqKiIhQXF7cL4Xr88cdRWFiI9957D8ePH8eqVauwZcsWHDlyBC+99BL27NkDu92OkpIS3HTTTckTaysuGiFqD7JL414D8Jre66bGuzEEjY2NKC0txdixY1FQUIBPP/004vjq6uqom2OJStW12my6j55ioUYjI8OqHYlEr9AC+sRWFVfVv6sWFO/Vq1fIVuaDBw/WitX36tULBQUFqKqqajcmuAbC2LFjUVlZifPnz+PYsWOYOXMmMjIyIEkSrrnmGrz33nu670c0nAOKrMT9EEXPfkdG4IEHHtCywg4ePBi1CV51dXXUugiKosQkjHoPIoBRAgtEdxmEQm1vEw31YSjLMux2O2RZxgsvvBB1P6CyshL79+/vVBth8uTJePfdwB7Onj178M033+Ds2bOYMGECysvLUVdXB6fTiQ8//DDJaiNwUW6EhJCSboTm5maUl5fjtddeAxCwSNQq+uGoqqrCiBEjIo75rsCuuj0VI4Q1mFgFNhg9li0QEFs17MvhcMDtduP06dNaq/NQtLa24rbbbsNzzz3XqTrYz3/+czzwwAMoLCzExIkTcdVVV0GSJBQUFOBnP/sZiouLkZWVhcmTJydZbQRAFtHLPEEkzysfR77++msMGDAAd999Nw4ePIiioiKsXr06YluSc+cC2UEej8fAlcZGn4H9RC8hLhgtsCpZ/Qd36/fLOYfb7Y56Dc45/H4//H6/luTQt2/fsJEMPp8Pt912G5YsWYJbb72108979+6NV199Vbv2iBEjNMNgxYoVWLFiBQDgl7/8ZdhiN2YkgXG2QkhJsZVlGfv27cPzzz+PGTNm4IEHHsCTTz6J3/zmN53GlpSUoLa2FpWVldi8eTNycnLw9NNPh7lyynpluo0ogQWAb85fikc/HzUuPSKcczQ1NeGLL76IOM5isWjFaNLT0yHLMvr16xdSaDnnWLFiBQoKCvDQQw+FvF5jY6MWq/vyyy9jzpw5mvV74cIF5Obm4vTp03j33Xej7k2YiUC6bvJ3G1FJSbHNz89Hfn6+5vsqLS3Fk08+GXJsWVkZAGDGjBn48MMPkZ6eHuHKB+K91B6NSIEFQlfm6i67du2K2n353LlzaGxsBGMMkiTB5/OF/Xi/Y8cOvPHGG5g4cSIKCwsBBKIPTp8+DQC47777cOzYMdx1112wWq0YN24cXnnlFe382267DXV1dbDZbFizZk1yhS9ywC+TZZvUDBo0CFdccQUqKiowZswYbNmyJWI1JM45PB4P0tLSDFxlz0S0wKrEU2A7Esn3CgR8u6qrQd0oC+fCmj17dtQol1mzZuHLL78M+bPt27frXLX5UDPsegopKbYA8Pzzz2PJkiXwer0YOXKk5vMKRSpljyUCswgskFiRBQJxsbIshwzjUklLS4PX69VqKagpu0QHOIdCboTkp7CwEHv37o0+EEBtbS369esZm09GYSaBBRIvsipqREIksVVdB8GddSNtzqYqapxtT4F2dHRQVVWVdNljRtJ/SF/tMCoONhaMElpAfxYZEPik5PV6oSgKsrOzO43rbqpuWVkZxowZg9GjR4fdkzA3HH5FifshipS1bGMhGVN140n/IX1FL6FLGCmyKnqzyKxWK/x+PzweD2RZRt++nV/j7qTqKoqCVatWYdOmTcjPz8e0adOwYMGCpOrU0NpY8a/t//xO/wRcujYB14wKia0OampqeqzYJquQhmNG/3LkDJ4qLF1Zb2KDOk5tYx7Ksh08eLCWIh6cqhssmEePHsUvfvELAO1Tdb/++muMHj0aI0eOBAAsWrQI69atSyqx5ZyXiF5DPCGx1UFVVRVGjx4dccyNS80X9tXThDQU/zHf3e7rzz+3QlEUoWKrdmIIh8Vi0Xy6Ho8HFRUVcLlcEc+Jlqo7e/bsdqm6we3NgUC44+7dIZsNEAZBYquDmpoazJ4927D5Bo4Y0qXzevXp+ZssHcW1I1arVet+IAK9lm1ws8fTp0+jvj58QkVXUnXDtTcnxJFyYjt8+HD06tULVqsVkiTpikjQ47PtqkASkYkmrh1RxVYUsbgRVKxWa9hzupqq63Q62xWdOXv2LIYMob9RkaSc2ALAxx9/jP799fvdqdGjccQqrh0RLbZ6K3+p1qfNZsOgQYNCfnLqTqrutGnT8OWXX+LUqVPIy8vD2rVr8fe//73b90d0nZQU21jgnEdtdUJ0ne6Ka0fU9uCi0Cv2wWJrsViQlZXVaUx3UnUlScIf//hHzJs3D4qiYPny5Rg/fnwc75SIlZQTW8YYbrjhBjDGcO+99+Kee+6JON7v92sFn4nuEW9hDYVoy1YvwR0bwrUx726q7vz58zF//vy4rJfoPikntjt27MCQIUNw4cIFFBcXY+zYsZgzZ07Y8RcvXsQtt9yCXbt2hR3z4oaRiVhqTJhxc2xa30/afR3hJYwbHo8HNTU1+OqrrxI/WRja2trw6aefRnxA+3w+cM4hSRK1MU8RUk5s1U2C3NxcLFy4EHv27AkrtiUlJTh9+jTOnz+P9evXo1+/fli3bl2ncS9uuJDQNZud5398Weh37foEM2fOvPTVzNAnJJDTp09DkiShm0EHDhzA2LFjIxYuUhQFu3fvhiRJsNlsJLYpQEqJbVtbm9bvqa2tDRs3bsSjjz4adnxZWRk++OADbNu2DY899piBKxVPsIAmE6J9tsDliIRIYqtmkTHGYLfbSWxTgJQS2/Pnz2PhwoUAAgXE77zzTpSURE5S6Ul1EZJVQGMhUhiVUQSXUIw2zu/3Iy0tLWwb87vuugvnzp2DxWLBPffcgwceeKDdmKamJixduhSnT5+GLMv4yU9+grvvvhsA8Prrr+O3v/0tAOA///M/sWzZsjjcHdFVUkpsR44ciYMHD8Z0TnV1NcaMGZOgFenjtd/kRh2jKAo+//xzTJ8+3YAVmRczbJDZ7Xb4opQGtFgsWrPHcGKrpzbCmjVrMG7cOLz//vu4ePEixowZgyVLlqC1tRWPPfYY9u7dC8YYioqKsGDBguQqHt7DSCmx7Qo1NTW45pprIo559dcDsHv37iBfpfEoiqJVk0plzCK20dJvgUBMrsfjwddffx1SbPXURmCMoaWlBZxztLa2om/fvpAkCf/6179QXFysFbgpLi5GWVkZFi9eHKe7JGKFSixGoaamJupmi8/nE5YeqkJiG8AsYhtLMZrXX389atfbcLUR7r//fhw7dgxDhgzBxIkTsXr1algslpC1Eaqqqrp2Q0RcILGNwrlz5zBw4MCIY7xeb9RW6IlGUZSkalOdKJJJbB0OB3w+H9ra2tDU1BR2XKTaCP/6179QWFiI6upqHDhwAPfffz+am5upNoIJIbGNAOccPp8vqpCaRWxFVboyE8kmtrIsY9CgQWEf6NFqI7z66qu49dZbwRjD6NGjMWLECBw/fhz5+flUG8Fk0LszAoqi6MoeM4vYkhvBHGJrs9mibpABl8XW4XCEDBPTUxth6NCh2LJlC4BAtE1FRQVGjhyJefPmYePGjWhoaEBDQwM2btyIefPmde/GiG5BnzsjcOHCBQwYMCDqOBJb82AGsdX7cV3t1hBuvJ7aCI888gh+8IMfYOLEieCc43e/+51WZOmRRx7BtGnTAACPPvpoyG4QhHGklNgqioKpU6ciLy8PH3zwQdTxetvheL3ekJX2jYR8tgEsFgv8fr/oZQCI3tJcrY8Qboye2ghDhgzBxo0bQ/5s+fLlWL58uf4FEwklpdwIq1evRkFBge7xsYgtWbbmwCybQHpcCepaI2WaET2HlBHbs2fPYsOGDVi5cqXuc/Rmj5mhBCOJrbnQW9cWQEy1lYnkJWXE9sEHH8Tvf//7mHbsq6urdRUN93q9FGdLtMNms0UVW4vFAlmWQ4qtnjbmf/jDH1BYWIjCwkJMmDABVqtVa6+zfPly5ObmYsKECfG5IaLbsGg+oQ7ENNgsfPDBB/jwww/xv//7v/jkk0/w1FNPRfTZlpSUoLa2FmfOnEHv3r2Rk5ODJ598Muz4trY2ZGaKrTvg8Xi0Vj8iMcNr0dbWhoyMDKEuBY/H066xYzhqa2vx6aef4tlnn233/ZqaGtTU1LRL1f3nP/8Ztjvu+++/j2effRZbt24FAJSXlyMrKwt33XUXjhw5Ep+bih1z+HRMQkrsqOzYsQPr16/Hhx9+CLfbjebmZixduhRvvvlmyPFlZWUAgJtuugnPPvss8vPzI15/165dQlN1AaCiogK5ubnCc9/N8Frs2bMHRUVFQi396upqyLKMoUOHRhxXXl4eUpD1pOoG89Zbb7VLxZ0zZw4qKyu7fgNE3EkJN8ITTzyBs2fPorKyEmvXrsV1110XVmiDOX/+PHJzIxeBMcvHd7OswwyYIfxLb2LDF198gX//+98Rx4RL1VVxOp0oKyvDbbfd1qW1EsaQEmLbFTjnkGU56sfAtrY24W9sAGhubha9BNPg8Xh0lThMJJIkoaGhIeq43r1748KF8MXnI6Xqqrz//vv49re/TXG0JiflxHbu3Lm6YmzV9FeLxaJlkYU69u/fjxdffDHiGCOOxx57DBcvXhS+jp///OfC1/DMM8/gyy+/FLqGuro6PProo1HH9e7dG1OmTAn5NxgtVVdl7dq1VM0rCUiJDbKuIMsyZs+eHbXqf2NjI5xOp/C88xMnTmDkyJHCN8i++eYbDBs2TPga+vfvL3Sjzu/348SJExg7dmzEcT6fD71798b27dvbfZ9zjmXLlqFv37547rnnwp7f1NSEESNG4MyZM53ut7KyEjfffDNtkJkFznksB9GBEydO8H379oleBn/33Xe5z+cTvQw+b9480Uvg5eXlvLq6WvQy+D/+8Y8un7t9+3YOgE+cOJFPnjyZT548mW/YsIG/8MIL/IUXXtDGvfrqq/yOO+7odP6iRYv4oEGDuCRJPC8vj7/88stdXks3iFVfevRBli1BEImCLNsgUs5nSxAEIQISW4IgCAMgsSUIE6InXRcAPvnkExQWFmL8+PHteuVRuq756PE+W7fbjTlz5sDj8UCWZZSWluKxxx4TvSyCiIiedN3GxkZcffXVKCsrw9ChQ3HhwgUtCYfSdc1Hj7dsHQ4Htm7dioMHD+LAgQMoKyvDrl27RC+LICIyePBgLf42OF03mL///e+49dZbtZTg4GzHOXPmUJKDyejxYssYQ1ZWFoBATKPP5wNj9MAlkodw6bonTpxAQ0MD5s6di6KiIvz1r38VtEJCDylRiEZRFBQVFeHkyZNYtWpV2BxzgjAbkdJ1ZVnG559/ji1btsDlcmHWrFmYOXMmvvWtbwlaLRGJHm/ZAoHCJAcOHMDZs2exZ88ekT4sgtBNtHTd/Px8lJSUIDMzE/3798ecOXNw8OBBASsl9JASYqvSp08fzJ07VyuhSBBmhevorHvLLbdg+/btkGUZTqcTu3fvjqntE2EsPV5sL168iMbGRgCAy+XC5s2bo+arE4Ro1M66W7du1boxfPjhh3jxxRfx4osvAgAKCgpQUlKCSZMmYfr06Vi5cqUW6rV48WLMmjULFRUVyM/PxyuvvCLydgikQOjXoUOHsGzZMiiKAr/fj9tvvx2PPvqo6GURRCpAO9FB9HixJQhCGCS2QfR4NwJBEIQZILElCIIwABLbbvDss89i/PjxmDBhAhYvXgy32234GioqKrQNlMLCQvTu3TtiselE0tjYiNLSUowdOxYFBQX49NNPhaxj+PDhmDhxIgoLCzF16lQhawAC8d1XXXUVbr75ZmFrIMwD+Wy7SFVVFWbPno2jR48iPT0dt99+O+bPn48f/OAHwtakKAry8vKwe/duId0Sli1bhu985ztYuXIlvF4vnE4n+vTpY/g6hg8fjr1796J///6Gzx3MM888g71796K5uVlXK6YeCPlsgyDLthvIsgyXy6XFOYpujbNlyxaMGjVKiNA2NzejvLwcK1asABDoLitCaM3C2bNnsWHDBqxcuVL0UgiTQGLbRfLy8vCTn/wEQ4cOxeDBg5GdnY0bbrhB6JpENv77+uuvMWDAANx999246qqrsHLlSrS1tQlZC2MMN9xwA4qKivDnP/9ZyBoefPBB/P73v4fFQm8xIgD9JXSRhoYGrFu3DqdOnUJ1dTXa2trw5ptvCluP1+vF+vXr8b3vfU/I/LIsY9++ffjhD3+I/fv3IzMzE08++aSQtezYsQP79u3DRx99hDVr1qC8vNzQ+T/44APk5uaiqKjI0HkJc0Ni20U2b96MESNGYMCAAbDZbLj11luxc+dOYev56KOPMGXKFAwcOFDI/Pn5+cjPz9eK/JSWlmLfvn1C1qK6c3Jzc7Fw4ULs2bPH0Pl37NiB9evXY/jw4Vi0aBG2bt2KpUuXGroGwnyQ2HaRoUOHYteuXXA6neCcY8uWLULz0t966y1hLgQAGDRoEK644gpUVFQACPiPgwtdG0VbWxtaWlq0/2/cuNHwbgVPPPEEzp49i8rKSqxduxbXXXed0E89hDlIiRKLiWDGjBkoLS3FlClTIEkSrrrqKtxzzz1C1uJ0OrFp0yb86U9/EjK/yvPPP48lS5bA6/Vi5MiRePXVVw1fw/nz57Fw4UIAAdfGnXfeiZKSEsPXQRAdodAvgiASBYV+BRGrZUsvHkEQRBcgny1BEIQBkNgSBEEYAIktQRCEAZDYEgRBGACJLUEQhAGQ2BIEQRgAiS1BEIQBkNgSBEEYAIktQRCEAZDYEgRBGMD/B9b4S/hfaCsxAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x11b325da0>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "from mpl_toolkits.mplot3d import Axes3D\n",
    "import matplotlib.pyplot as plt\n",
    "from matplotlib import cm\n",
    "from matplotlib.ticker import LinearLocator, FormatStrFormatter\n",
    "import numpy\n",
    "from scipy.stats import norm\n",
    "\n",
    "fig = plt.figure()\n",
    "ax = fig.gca(projection='3d')\n",
    "\n",
    "slice_n = 0\n",
    "Xarray = numpy.arange(3.0, 8.0, 1.0) # 2.0 is OK\n",
    "Yarray = numpy.arange(4.0, 9.0, 1.0) # 4.0 is OK\n",
    "X, Y = numpy.meshgrid(Xarray, Yarray)\n",
    "Z = X + 0.0\n",
    "W = X + 0.0\n",
    "\n",
    "for i,num_slices_pow in enumerate(list(Yarray)):\n",
    "    depth = 2 ** int(num_slices_pow)\n",
    "\n",
    "    for j,slice_width_pow in enumerate(list(Xarray)):\n",
    "        factor = 2 ** int(slice_width_pow)\n",
    "        width = depth * factor\n",
    "        W[i][j] = width\n",
    "\n",
    "        xval = norm.cdf(\n",
    "            numpy.linspace(\n",
    "                norm.ppf(1/width),\n",
    "                norm.ppf(1 - 1/width),\n",
    "                width))\n",
    "\n",
    "        r1 = xval[slice_n*factor : (slice_n+1)*factor]\n",
    "        r2 = xval[(slice_n+1)*factor : (slice_n+2)*factor]\n",
    "\n",
    "        u1 = numpy.mean(r1)\n",
    "        u2 = numpy.mean(r2)\n",
    "        o1 = numpy.std(r1)\n",
    "        o2 = numpy.std(r2)\n",
    "\n",
    "        d = abs(u1 - u2) / numpy.sqrt(4 * o1 * o2)\n",
    "        dd = numpy.sqrt(d*d - 1)\n",
    "        t = 2 * numpy.log(d - dd) + 2 * d * dd\n",
    "        \n",
    "        Z[i][j] = t\n",
    "\n",
    "# Plot the surface.\n",
    "surf = ax.plot_surface(X, Y, Z, cmap=cm.coolwarm,\n",
    "                       linewidth=0, antialiased=False)\n",
    "\n",
    "# Customize the z axis.\n",
    "ax.zaxis.set_major_locator(LinearLocator(10))\n",
    "ax.zaxis.set_major_formatter(FormatStrFormatter('%.02f'))\n",
    "ax.view_init(30, 180)\n",
    "\n",
    "# Add a color bar which maps values to colors.\n",
    "fig.colorbar(surf, shrink=0.5, aspect=5)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "ename": "SyntaxError",
     "evalue": "invalid syntax (<ipython-input-11-d862d2fd1795>, line 61)",
     "output_type": "error",
     "traceback": [
      "\u001b[0;36m  File \u001b[0;32m\"<ipython-input-11-d862d2fd1795>\"\u001b[0;36m, line \u001b[0;32m61\u001b[0m\n\u001b[0;31m    slice1, slice2 =\u001b[0m\n\u001b[0m                     ^\u001b[0m\n\u001b[0;31mSyntaxError\u001b[0m\u001b[0;31m:\u001b[0m invalid syntax\n"
     ]
    }
   ],
   "source": [
    "class HistogramSlice(object):\n",
    "    test_threshold = 16\n",
    "    test_half_threshold = 8\n",
    "    \n",
    "    def __init__(self, count=0, full_sum=0.0, squared_sum=0.0):\n",
    "        self.count = count\n",
    "        self.sum = full_sum\n",
    "        self.is_bimodal = False\n",
    "        \n",
    "        half_count = count // 2\n",
    "        divisor = half_count / count\n",
    "        half_sum = full_sum * divisor\n",
    "        half_squared_sum = squared_sum * divisor\n",
    "        \n",
    "        self.left_count = half_count\n",
    "        self.left_sum = half_sum\n",
    "        self.left_squaredsum = half_squared_sum\n",
    "        \n",
    "        self.right_count = half_count\n",
    "        self.right_sum = half_sum\n",
    "        self.right_squaredsum = half_squared_sum\n",
    "        \n",
    "    def insert(self, val):\n",
    "        assert self.is_bimodal is False\n",
    "        \n",
    "        self.count += 1\n",
    "        self.sum += val\n",
    "        mean = self.sum / self.count\n",
    "\n",
    "        if val < mean:\n",
    "            self.left_count += 1\n",
    "            self.left_sum += val\n",
    "            self.left_squaredsum += val*val\n",
    "        else:\n",
    "            self.right_count += 1\n",
    "            self.right_sum += val\n",
    "            self.right_squaredsum += val*val\n",
    "            \n",
    "        # TODO:\n",
    "        #    add count-min-sketch for median?\n",
    "        \n",
    "    @property\n",
    "    def is_bimodal(self):\n",
    "        if self.is_bimodal:\n",
    "            return True\n",
    "        \n",
    "        if ((self.count < self.test_threshold)\n",
    "                or (self.test_left_count < self.test_half_threshold)\n",
    "                or (self.test_right_count < self.test_half_threshold)):\n",
    "            return False\n",
    "        \n",
    "        # TODO: all this jazz can be algebraically optimized\n",
    "        u1 = self.left_sum / self.left_count \n",
    "        u2 = self.right_sum / self.right_count \n",
    "        v1 = (self.left_squaredsum - u1 * self.left_sum) / (self.left_count - 1)\n",
    "        v2 = (self.right_squaredsum - u2 * self.right_sum) / (self.right_count - 1)\n",
    "\n",
    "        d = 0.5 * abs(u1 - u2) / ((v1 * v2) ** 0.25)\n",
    "        dd = numpy.sqrt(d*d - 1)\n",
    "        t = 2 * numpy.log(d - dd) + 2 * d * dd\n",
    "        \n",
    "        # correct t for standard normal of slice\n",
    "        # compare against threshold\n",
    "        \n",
    "        self.is_bimodal = ...\n",
    "        return self.is_bimodal\n",
    "    \n",
    "    def split(self):\n",
    "        return (Slice(self.left_count, self.left_sum, self.left_squaredsum),\n",
    "                Slice(self.right_count, self.right_sum, self.right_squaredsum))\n",
    "\n",
    "class HistogramTree(object):\n",
    "    def __init__(self):\n",
    "        pass\n",
    "        \n",
    "    def insert_point(self, val):\n",
    "        pass\n",
    "    \n",
    "        # 1.) find non-bimodal slice\n",
    "        # 2.) insert into it\n",
    "        # 3.) check it for bimodality\n",
    "        # 4.) split if necessary\n",
    "        # 5.) get bin in slice\n",
    "        # 6.) add to bin, checking for anomaly\n",
    "    \n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 100,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "nan"
      ]
     },
     "execution_count": 100,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "%matplotlib inline\n",
    "import numpy\n",
    "\n",
    "rand = numpy.random.RandomState(seed=1234567)\n",
    "\n",
    "size = 1024\n",
    "factor = 128\n",
    "xval = sorted(\n",
    "    list(rand.normal(1e-3, 2e-4, size))\n",
    "  + list(rand.normal(2.5e-3, 4e-4, size)))\n",
    "\n",
    "ty = []\n",
    "tx = range(int(len(xval)/factor) - 1)\n",
    "for i in tx:\n",
    "    r1 = xval[i*factor : (i+1)*factor]\n",
    "    r2 = xval[(i+1)*factor : (i+2)*factor]\n",
    "\n",
    "    u1 = numpy.mean(r1)\n",
    "    u2 = numpy.mean(r2)\n",
    "    o1 = numpy.std(r1)\n",
    "    o2 = numpy.std(r2)\n",
    "\n",
    "    d = abs(u1 - u2) / (numpy.sqrt(4 * o1 * o2))\n",
    "    dd = numpy.sqrt(d*d - 1)\n",
    "    t = 2 * numpy.log(d - dd) + 2 * d * dd\n",
    "    ty.append(t)\n",
    "\n",
    "ly = []\n",
    "lx = numpy.linspace(0.0, 1.0, 100)\n",
    "for p in lx:\n",
    "    l = numpy.log(1 - p) - numpy.log(p)\n",
    "    ly.append(numpy.abs(l))\n",
    "    \n",
    "ax = matplotlib.pyplot.subplot(2, 2, 1)\n",
    "ax.set_ylim(0, 20)\n",
    "matplotlib.pyplot.plot(tx, ty)\n",
    "\n",
    "ax = matplotlib.pyplot.subplot(2, 2, 2)\n",
    "ax.set_ylim(0, 20)\n",
    "matplotlib.pyplot.plot(lx, ly)\n",
    "\n",
    "ax = matplotlib.pyplot.subplot(2, 2, 3)\n",
    "matplotlib.pyplot.hist(xval)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "%matplotlib inline\n",
    "import numpy\n",
    "import numpy.random\n",
    "\n",
    "rand = numpy.random.RandomState(seed=1234567)\n",
    "\n",
    "size = 1024\n",
    "factor = 128\n",
    "xval = numpy.random.shuffle(\n",
    "           numpy.concatenate((\n",
    "               rand.normal(1e-3, 2e-4, size),\n",
    "               rand.normal(2.5e-3, 4e-4, size)\n",
    "           )))\n",
    "           "
   ]
  }
 ],
 "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.4"
  },
  "toc": {
   "base_numbering": 1,
   "nav_menu": {},
   "number_sections": true,
   "sideBar": true,
   "skip_h1_title": false,
   "title_cell": "Table of Contents",
   "title_sidebar": "Contents",
   "toc_cell": false,
   "toc_position": {},
   "toc_section_display": true,
   "toc_window_display": false
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}