{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Chapter 10. Big Entropy and the Generalized Linear Model"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "!pip install -q numpyro arviz"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import os\n",
    "\n",
    "import arviz as az\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "import jax.numpy as jnp\n",
    "from jax import random, tree_map, vmap\n",
    "\n",
    "import numpyro\n",
    "import numpyro.distributions as dist\n",
    "\n",
    "if \"SVG\" in os.environ:\n",
    "    %config InlineBackend.figure_formats = [\"svg\"]\n",
    "az.style.use(\"arviz-darkgrid\")\n",
    "numpyro.set_platform(\"cpu\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Code 10.1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "p = {}\n",
    "p[\"A\"] = jnp.array([0, 0, 10, 0, 0])\n",
    "p[\"B\"] = jnp.array([0, 1, 8, 1, 0])\n",
    "p[\"C\"] = jnp.array([0, 2, 6, 2, 0])\n",
    "p[\"D\"] = jnp.array([1, 2, 4, 2, 1])\n",
    "p[\"E\"] = jnp.array([2, 2, 2, 2, 2])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Code 10.2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "p_norm = tree_map(lambda q: q / jnp.sum(q), p)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Code 10.3"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "{'A': DeviceArray(-0., dtype=float32),\n",
       " 'B': DeviceArray(0.6390318, dtype=float32),\n",
       " 'C': DeviceArray(0.95027053, dtype=float32),\n",
       " 'D': DeviceArray(1.4708084, dtype=float32),\n",
       " 'E': DeviceArray(1.609438, dtype=float32)}"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "H = tree_map(lambda q: -jnp.sum(jnp.where(q == 0, 0, q * jnp.log(q))), p_norm)\n",
    "H"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Code 10.4"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "ways = jnp.array([1, 90, 1260, 37800, 113400])\n",
    "logwayspp = jnp.log(ways) / 10"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Code 10.5"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "{1: DeviceArray(1., dtype=float32),\n",
       " 2: DeviceArray(1., dtype=float32),\n",
       " 3: DeviceArray(1., dtype=float32),\n",
       " 4: DeviceArray(1., dtype=float32)}"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# build list of the candidate distributions\n",
    "p = {}\n",
    "p[1] = jnp.array([1 / 4, 1 / 4, 1 / 4, 1 / 4])\n",
    "p[2] = jnp.array([2 / 6, 1 / 6, 1 / 6, 2 / 6])\n",
    "p[3] = jnp.array([1 / 6, 2 / 6, 2 / 6, 1 / 6])\n",
    "p[4] = jnp.array([1 / 8, 4 / 8, 2 / 8, 1 / 8])\n",
    "\n",
    "# compute expected value of each\n",
    "tree_map(lambda p: jnp.sum(p * jnp.array([0, 1, 1, 2])), p)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Code 10.6"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "{1: DeviceArray(1.3862944, dtype=float32),\n",
       " 2: DeviceArray(1.3296614, dtype=float32),\n",
       " 3: DeviceArray(1.3296614, dtype=float32),\n",
       " 4: DeviceArray(1.2130076, dtype=float32)}"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# compute entropy of each distribution\n",
    "tree_map(lambda p: -jnp.sum(p * jnp.log(p)), p)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Code 10.7"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "DeviceArray([0.09, 0.21, 0.21, 0.49], dtype=float32)"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "p = 0.7\n",
    "A = jnp.array([(1 - p) ** 2, p * (1 - p), (1 - p) * p, p**2])\n",
    "A"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Code 10.8"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "DeviceArray(1.2217286, dtype=float32)"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "-jnp.sum(A * jnp.log(A))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Code 10.9"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "def sim_p(i, G=1.4):\n",
    "    x123 = dist.Uniform().sample(random.PRNGKey(i), (3,))\n",
    "    x4 = (G * jnp.sum(x123, keepdims=True) - x123[1] - x123[2]) / (2 - G)\n",
    "    z = jnp.sum(jnp.concatenate([x123, x4]))\n",
    "    p = jnp.concatenate([x123, x4]) / z\n",
    "    return {\"H\": -jnp.sum(p * jnp.log(p)), \"p\": p}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Code 10.10"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAtsAAAHrCAYAAAAe4lGYAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAAPYQAAD2EBqD+naQAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAgAElEQVR4nOzdd5RV1dnH8e+Zoc4w4NCEIIhUARFEFMUaRWKvMRo1xtgTNW9MjDWamBjRmBgbsWuM0USxd7FiVDoiKiBKFZQiZWgCU/b7x2YYEVDQuXPLfD9rse6555w7dzt7udaPh+fsnYQQApIkSZKqXV66ByBJkiTlKsO2JEmSlCKGbUmSJClFDNuSJElSihi2JUmSpBQxbEuSJEkpYtiWJEmSUsSwLUmSJKVInXQP4JssXrx4g3NNmjShpKQkDaPR5nB+Mpvzk9mcn8zl3GQ25yez5er8FBcXf+M9WVnZzsvLymHXGs5PZnN+Mpvzk7mcm8zm/GS22jw/tfe/XJIkSUoxw7YkSZKUIoZtSZIkKUUM25IkSVKKGLYlSZKkFDFsS5IkSSli2JYkSZJSxLAtSZIkpYhhW5IkSUoRw7YkSZKUIoZtSZIkKUUM25IkSVKKGLYlSZKkFDFsS5IkSSli2JYkSZJSxLAtSZIkpYhhW5IkSUoRw7YkSZKyTkVF4J3xgZUrQ7qH8rUM25IkSco6I0fBeb8K3DzYsC1JkiRVq9lz4uvKlekdxzcxbEuSJCnrLF8eXxs1Su84volhW5IkSVln+fLYPlJo2JYkSZKq1/IV8bWoUZLegXwDw7YkSZKyjm0kkiRJUopUhu3CwvSO45sYtiVJkpR1KsN2kZVtSZIkqXpZ2ZYkSZJSxJ5tSZIkKQVCCF9ajSS9Y/kmhm1JkiRllVWroLw8HttGIkmSJFWjFWur2vl50LBhesfyTQzbkiRJyirLvtSvnSRuaiNJkiRVm3UrkWR4vzYYtiVJkpRlKh+OzPSVSMCwLUmSpCyTLRvagGFbkiRJWSZbNrQBw7YkSZKyTLZsaAOGbUmSJGWZ5csDYBuJJEmSVO2q2kgye9k/MGxLkiQpy7gaiSRJkpQirkYiSZIkpYirkUiSJEkpssI2EkmSJCk1lrn0nyRJkpQaVrYlSZKkFCgtDaxaFY8N25IkSVI1qlz2D6CwIH3j2FyGbUmSJGWNFWv7tQsKID/fTW0kSZKkarM8ix6OBMO2JEmSssi6lUiyYI1tMGxLkiQpi2TTSiRg2JYkSVIWsY1EkiRJShHbSCRJkqQUWbEiAFa2JUmSpGqXbW0kdbb0A08++SRjx47l/fffZ8qUKZSWljJo0CCOPvrojd6/fPlybr75ZoYOHcqCBQto0aIFAwcO5LzzzqNRtvyWJEmSlBHWtZE0yvw1tuFbhO0bb7yROXPmUFxcTMuWLZkzZ84m7125ciUnnXQSkyZNYo899uCQQw5h8uTJ/POf/2TkyJE8+OCDFBRkwdY/kiRJygg5vxrJVVddxauvvsqIESM4/vjjv/beu+66i0mTJnH66adzzz33cMEFF3DXXXdxzjnnMGnSJO66665vPXBJkiTVPpVtJIW5+oBk//79adOmzTfeF0JgyJAhFBQUcM4556x37ayzzqJJkyY88sgjhBC2dAiSJEmqpSrbSIpytbK9uWbMmMH8+fPp06fPBq0i9evXp2/fvsybN4+ZM2emagiSJEnKMSuy7AHJlIXtyhDdvn37jV7fdttt17tPkiRJ+ibLs2yd7S1+QHJzLVu2DGCTK45Unq+8b1OaNGlCXt6GfycoLi7+jiNUKjk/mc35yWzOT+ZybjKb85PZqmN+KioCK1YuAqDNNltRXJz5q1inLGxXl5KSkg3OFRcXs3jx4jSMRpvD+clszk9mc34yl3OT2ZyfzFZd87N8eaDycb+y0iUsXpze5f825y8QKfvrQFFRERDX2d6YyvOV90mSJElfpzJW1qsL9etnxzrbKQvblT3ZM2bM2Oj1yl7tyvskSZKkr7M8y9bYhhSG7fbt29OyZUvGjRvHypUr17u2evVqxowZQ8uWLQ3bkiRJ2izZtlU7pDBsJ0nCsccey8qVKxk8ePB6126//XZKSko49thjSZLs+CcASZIkpde6DW2yKGxv8QOSQ4YMYezYsQBMmTJl3blRo0YBMGDAAAYMGADA6aefzquvvrpuJ8kePXowefJk3njjDbp168bpp59eXf8dkiRJynGVbSTZsqENfIuwPXbsWB5//PH1zo0bN45x48YB0KZNm3Vhu6CggPvvv59bbrmFF198kVGjRtG8eXNOOeUUzj333A02u5EkSZI2JRvbSLY4bF9zzTVcc801m31/UVERl1xyCZdccsmWfpUkSZK0zmdz47p/2RS2M38lcEmSJNV6K1YEnn0uHvfbJXue+TNsS5IkKeM9+XRsI2nXFvbaM92j2XyGbUmSJGW0NWsCDw2JLSQnnpCQl2dlW5IkSaoWz78ICxdCyxYwcEC6R7NlDNuSJEnKWGVlgQf/E6vax/8ooW7d7Klqg2FbkiRJGWzYGzDnU2jcGA47NN2j2XKGbUmSJGWsl1+JVe2jj4SGDbOrqg2GbUmSJGWoiorA+AnxuP/u2Re0wbAtSZKkDDV1KixbBg0bQpfO6R7Nt2PYliRJUkZ659342mtHqFPHyrYkSZJUbd4ZH/u1e/fKzqANhm1JkiRloIqKwLtr+7V36p3esXwXhm1JkiRlnGnTYelSaNgAunZJ92i+PcO2JEmSMs474+Nrz57Z268Nhm1JkiRloPHvxn7tnXpnb9AGw7YkSZIyTEVFYPzaynY292uDYVuSJEkZZvoMKFkKDRrA9l3TPZrvxrAtSZKkjFJZ1e65Q3b3a4NhW5IkSRnmnRzp1wbDtiRJkjLMxInxtecO6R1HdTBsS5IkKWN8vjAwfwHk5WX3+tqVDNuSJEnKGJMnx9dt20FBgW0kkiRJUrWZNDn2a3frluaBVBPDtiRJkjLGpLWV7W7bZ39VGwzbkiRJyhAhBCZ/GI+7bZ/esVQXw7YkSZIywqefwtKlULcudOyQ7tFUD8O2JEmSMsLEtS0knTtB3bq2kUiSJEnf2mefBaZNC+veT658ODJHWkjAsC1JkqQ0qKgInPN/gZ+dHtaF7ImT4rVu3XKjqg2GbUmSJKXBZ3Nh/nwor4Drrg+sWROY8lG81t3KtiRJkvTtfTy16vjDKfD3mwKrV0NhIWyzTfrGVd0M25IkSapxH38cW0e2ahLfP/1MfN2+K+Tl2UYiSZIkfWsffxxfTzoxYfuuVedz6eFIMGxLkiQpDT5aG7a7doELfp2QrC1m58rOkZXqpHsAkiRJql2WLQvMnRePO3aExkUJv/olTJgQ2K1fesdW3QzbkiRJqlGVD0e22joGbYBjjko45qjcqmqDbSSSJEmqYZVhu1On9I6jJhi2JUmSVKM+WrsSSWfDtiRJklS91lW2O+Ze28hXGbYlSZJUY8rKAtOnx2PbSCRJkqRvobQ0UFERNjg/axaUlkJBAbRulYaB1TDDtiRJkqrVrE8CPzgkcP2NG4btyvW1O3XMrZ0iN8WwLUmSpGo1YiSsWQMvvQzl5esH7o+nxvedOqZjZDXPsC1JkqRqVRmoV6yAadPWv1ZZ2e7cKfer2mDYliRJUjWbOrXqePyEquMQQq1aYxsM25IkSapGX15tBODdCVVtJLPnwJIlkJ8PHbZLw+DSwLAtSZKkavPJbFhTWvV+woRY0QZ4fVg8t3MfqF/fNhJJkiRpi1S2kHTpDHXrwqLFMGtWBQCvD4uhe999akfQBsO2JEmSqlHlw5Hdto9/AMaOK+XTzwIfToG8PNhrzzQOsIbVSfcAJEmSlDumrl19pGPHhMaNAxPeg7Hjyvhsbjy/U28o3qr2VLYN25IkSao2H39p05rWrRN4IDD2nVKmfFT7WkjAsC1JkqRqUlISWPB5PO7YAUKAJIFPPok920kCe9eiFhKwZ1uSJEnVpLKFpHVrKCxMaNQoWW+nyF47QrNmtauybdiWJElStfhyC0mlXjtWHde2FhIwbEuSJKmafDwt9mWvH7arAvY+e9X0iNLPnm1JkiRVi3VbsXesCti77gLdu+XTqVM5LVrUvsq2YVuSJEnfWVlZYMbabdq/XNkuLEwY8t+tWLx4cXoGlma2kUiSJGmzlZcHPv88bHC+cpv2hg3jA5KKaqSyHULgpZde4v7772f69OksW7aMVq1a0a9fP8444wzatm1bE8OQJEnSd3TbHYH/PAT/dx4ce0xVW8j4d+Nrxw6Ql1f72kU2pUYq29deey3nnXce06dPZ//99+ekk05im2224eGHH+aII45gypQpNTEMSZIkfY13JwRGjd6wal1pxYrAE0/G41sGB8aOi/d+9FFg8K3xeLd+Bu0vS3lle8GCBdx33320adOGp556ikaNGq279s9//pNBgwZx7733MmjQoFQPRZIkSZuwalXgNxcGVq2CO2+DbttvGJqHvgxfrIrH5RXw+ysD1w6Cy/8QP9d3ZzjphBoeeIZLeWV7zpw5VFRU0KdPn/WCNsC+++4LwKJFi1I9DEmSJH2NiZNg1dogfevtgRDWr3CHEHjiyXju7DMTunaBJSVw9jmB+fOhXVv44x8S6tSxsv1lKQ/b2267LXXr1mXcuHEsX758vWvDhg0DYLfddkv1MCRJkvQ1JrxXdTzuHRg5av3r770fd4isXx+OOAz+/KeErZrELdmLiuDaQQmNiwzaX5XyNpLi4mLOP/98/vKXv3DwwQez3377UVhYyJQpUxg+fDjHHXccJ510UqqHIUmSpK8x4b1YtW7RHBZ8DrfeEdh1l6qHHZ98Kl4fsD8UFSVrAzY8+N/Aj49LaLuNQXtjkvDVfyNIkaeffporrriClStXrju30047ccEFF9C3b99Nfq6iooK8PFcolCRJSpXy8kD/vRezfHng7jsa86vfLGPZssCgqxpx+GH1Wby4gu8fsJjSUvjvA03ouYNbtWyuGvlN/eMf/+Af//gH5557LkceeSSNGzdm0qRJXHPNNZx88snccMMNDBw4cKOfLSkp2eBccXFxrV0YPRs4P5nN+clszk/mcm4ym/Pz3Xz0cWD58kBBAXTssIwTfwy33QHX/nU5Q19ezpIlUFoKXbvANm2WsaW/6lydn+Li4m+8J+Ul4+HDh3PjjTdy4okncvbZZ9OqVSsKCgrYeeeduf3226lfv74rkUiSJKXRe+/H1x7doU6dhGOPgZYtYckSeO11eGd8vH7UEbaKbKmUV7YrH4Ls16/fBteaNm1K165deeedd1i0aBFNmzZN9XAkSZL0FZX92jv2jGG6fv2Ef9wcQ/bKFbBiJRQUwEEHpnOU2SnlYbu0tBTY9PJ+lefr1auX6qFIkiRpIypXIum5Q9W5VlsnHPSD9Iwnl6S8jaRPnz5A3MBm2bJl6117/PHHmTlzJj169NhgDW5JkiSl3tx5cZ3s/Dzo3i3do8k9Ka9sH3jggfz3v/9l1KhRDBw4kP3224/GjRvz4Ycf8tZbb1GvXj0uvfTSVA9DkiRJG1HZr925MxQU2JNd3VIetvPz87n77ru57777eP7553n22WcpLS2lWbNmHHrooZx11ll06dIl1cOQJEnSRry3tl/7yy0kqj41svRfvXr1OOOMMzjjjDNq4uskSZL0FZMmB2bMgPkLYMHngaIi2G7bhLHj4vXKhyNVvVyRXJIkKceNHhM4/4KN7WNYda5nz5obT21i2JYkScpxb7wZQ3XbtrFdpHkzKCmBGTNh5izouzM0b2ZlOxUM25IkSTluzNj4+ouzEvba01Bdk1K+9J8kSZLSZ+68wCefQF4e7NQ73aOpfQzbkiRJOayyqt1te2jUyKp2TTNsS5Ik5bAxY2O/9i590zyQWsqwLUmSlKMqKsK6ynbfna1qp4NhW5IkKUdNmw5LlkCDBtCje7pHUzsZtiVJknJUZVW7dy+oW9fKdjoYtiVJknLU6DGxX9sWkvQxbEuSJOWgNWsC706Ix313Tu9YajPDtiRJUg76YCKsWgXFxdCxQ7pHU3sZtiVJknLQ8y9UtpBAkthGki6GbUmSpBwzfUbghaHx+NhjDNrpZNiWJEnKMXffE6iogL32hO7dDNvpZNiWJEnKIZMnB15/A5IEzjjNoJ1uhm1JkqQccsfdsVd74AHQYTvDdroZtiVJknLEuHcCo0ZDfj6ceopBOxMYtiVJknLAihWBa6+LVe3DD4M23zNsZwLDtiRJUpYLIXDd9YE5n0Krre3VziSGbUmSpCz33Avw8iuQnwd/uCKhcZFhO1MYtiVJkrLYjJmBv98Y20dOPy1hhx4G7Uxi2JYkScpSy5YFLv1dYNWquFPkiT9O94j0VYZtSZKkLFRWFrjiysCsT6BlC7j80oS8PKvamcawLUmSlIVuHhwYPQYaNIBrr05o1sygnYkM25IkSVnmqWcCjz4ej6+4LKFzZ4N2pjJsS5IkZZl/3R8fiDzjtIS99zJoZzLDtiRJUhaZOzcwd15c5u/YY9I9Gn0Tw7YkSVIWeXdCfO3SFQoKrGpnOsO2JElSFhk/IbaQ9N4xzQPRZjFsS5IkZZF3342vvXpZ1c4Ghm1JkqQssWhRXFc7SWDHnukejTaHYVuSJClLVPZrd+gAjYusbGcDw7YkSVKWeNd+7axj2JYkScoS4+3XzjqGbUmSpCywdFlg6rR43Mt+7axh2JYkScoC770HIUDbttCsmZXtbGHYliRJyiCfzA7Mnx82OG+/dnaqk+4BSJIkKfpkduCnpwbKy+GIwwI/+2lCkyYwdhwMeyPeY792djFsS5IkZYgnnwqsWROPH3sCnn8x0LgxzJsXz+XnQ5/e6RuftpxtJJIkSRlg9erA8y/E45/9FLbvCl98EYN2o0I48gi4+46Eli2tbGcTK9uSJEkZ4I03oWQptGwBp5yc8LOfwqjRsHo17NYP6tc3ZGcjw7YkSVIGeOrp+ADkoYck5OfHYL1bv3SOSNXBNhJJkqQ0mzUr8M54yMuDQw5O92hUnQzbkiRJafbUM7GqvVs/2Nqe7JxiG4kkSVKalJcHJn/IugcjjzjMoJ1rDNuSJEk1rKQk8PebAiNHwbJl8VzLFtBv1/SOS9XPsC1JklTDnnoGXn4lHjcqhJ13jiuQ1KljZTvXGLYlSZJq2NvDY4/26acmnHQChuwc5gOSkiRJNaikJPDBxHh84A8M2rnOsC1JklSDRo2BigrosB202tqgnesM25IkSTVo+NoWkt13S/NAVCMM25IkSTWkvDyuQAKw+25WtWsDw7YkSVINmTgJSpZCo0awQ490j0Y1wbAtSZJUQ4aPiC0k/XbxwcjawrAtSZJUQ4aPiK+7727Qri0M25IkSTVgwYLARx9DkrhTZG3ipjaSJEkpMnxk4PY7AkuXVW3L3m17KN7KynZtUaNh+6WXXuLBBx9k4sSJfPHFFzRv3pzevXvz29/+ltatW9fkUCRJklJq1arAX64LLPh8/fOHHGzQrk1qJGyHEPj973/PQw89RLt27Tj44IMpLCxk/vz5jB49mjlz5hi2JUlSTvnvw7Dgc2i1NfzpyoTCQmjSGJo0MWzXJjUStu+//34eeughTjzxRC677DLy8/PXu15WVlYTw5AkSaoRCxcGHngwrjxy9pkJ3bY3YNdWKX9ActWqVQwePJi2bdty6aWXbhC0AerUsXVckiTljrvuDXyxCrp3g/33S/dolE4pT7lvvfUWS5Ys4aijjqKiooKhQ4cyY8YMioqK6N+/P9tuu22qhyBJklRjpk0LPPtcPD73FwlJYlW7Nkt52H7//fcByM/P5/DDD2f69OnrruXl5XHKKadw0UUXpXoYkiRJNeKOuwMVFbDv3rBjT4N2bZfysL1w4UIA7r33Xrp3786QIUPo2LEjkyZN4vLLL+eee+6hbdu2nHDCCRv9fJMmTcjL27Dbpbi4OKXj1nfj/GQ25yezOT+Zy7nJbJkwP9Oml/PmW0tIErjgN1tRXLxh+2xtlQnzkw4pD9shxIcD6taty+DBg9l6660B6Nu3LzfddBOHH34499577ybDdklJyQbniouLWbx4ceoGre/E+clszk9mc34yl3OT2TJlfu68uwKAPfrDVk2WkgFDygiZMj/VbXP+ApHyByQbNWoEwA477LAuaFfq3Lkzbdu2ZdasWSxdujTVQ5EkSUqZxUsCL7wYj4//ke0jilIetjt06ABAUVHRRq9Xnl+1alWqhyJJkpQyTzwJa9bA9l2h147pHo0yRcrbSPr16wfAtGnTNrhWWlrKrFmzKCgooGnTpqkeiiRJUrW5/4HA1GmBY45K6NIZHnsits4e9yNXIFGVlFe227Vrx5577snMmTMZMmTIetfuuOMOli5dyoABA1xrW5IkZY2SksDtdwZefgV+fm7glNMDixdDy5bw/X3SPTplkhpJuL///e85/vjj+d3vfsfLL79Mhw4dmDhxIiNGjKBNmzZceOGFNTEMSZKkajF+QnwtKIDVq+GTT+L7Y49JqFPHqraq1EjYbteuHY8++ig33XQT//vf/3jrrbdo3rw5J554Iueccw7NmjWriWFIkiRVi3fGx5aRHwyEE49P+O/DgWXL4IjD0jwwZZwa691o3bo1gwYNqqmvkyRJSpnx4+Nr714JrVol/OqXVrO1cSnv2ZYkScolS5cGpq5d92GnXukdizKfYVuSJGkTVq8OXHxZBffeF9ade3cChADbtoOmTa1o6+sZtiVJkjZhzFh48y2455+BmTNj4K7s196pdzpHpmxh2JYkSdqEiZNisA4BHvhPZdiO13bqbVVb38ywLUmStAkTJ1Udv/gSfDw18PHU+L63/draDIZtSZKkjaioCEyaHI+33hrKy+EPfwqEAO3aQrNmVrb1zQzbkiRJGzF7NixfDvXqwQXnx2A9Y0a81tt+bW0mw7YkSdJGfLC2haRrF9itH3TpUnXNfm1tLsO2JEnSRlQ+HNm9OyRJwk9OqArYrq+tzVVjO0hKkiRlk4kT42v3bjFk770XHHkENGkMzZtb2dbmMWxLkiR9xerVVauOdN8+vubnJ+t6t6XNZRuJJEnSV3z0cVx9pLgYWrVK92iUzQzbkiSpVgohsGBBoLQ0bHCtqoUk9mtL35ZtJJIkqVZ66hm47m+BenWha9dAj+5wxOEJbbdJ+KDy4chuBm19N4ZtSZJUKw17IwbqNaXw3vvxzxNPBc47p2rnyO7d0jhA5QTDtiRJymlDXwq88WbgkgsTCgtjpbq8PPDB2laRP/8xYeUX8NzzgXfGx2p3pe27pmPEyiWGbUmSlLOWLw/89e+BlSuhz05w9JHx/IyZsGIFNGwIe/SHOnUSfnAAPDQE7rgrUFoK27aDoiLbSPTd+ICkJEnKWU8/CytXxuPRo6sq1u+9H1+7d4tBGyAvL+HHxyXccWvCXnvCaacatPXdWdmWJEk5qawsMOTRqoA99p14rk6dhPc/iOd77rDh5zp3Shh0lUFb1cPKtiRJykmvDYP58+Na2U0axwp3ZZ/2+2sr2zv0MFQrtQzbkiQp54QQ+O9DsXp9zFEJffvG86NGBxYtCsyeE9/36J6mAarWMGxLkqSc8+4E+HAK1KsHRxwOu+4SK9ijxsD7H8R7tmvvA5BKPcO2JEnKOQ89HKvaB/0AirdK2HVtZXvyZHh7eLy2w0b6taXqZtiWJEk5pawsMHJ0PD7qyFi5btEiYbv2EAI8/2K81tN+bdUAw7YkScopsz6BNWugYQPosF3V+V13ia/l5fG1Z8+aH5tqH8O2JEnKKR9/HF87dYprZ1eq7NsG2KoJbNOmpkem2siwLUmScsqUj2JPdudO65/vtSPUqxuPd9gBksQ2EqWeYVuSJOWUj9ZWtjt3Wj9MN2iQ0Lt3PN6xp0FbNcMdJCVJUs4IIVSF7c4bXv/VeQkvDA0cdUTNjku1l5VtSZKUlcrKApddUcGNN1esOzdvPixdCvl5cR3tr2rXLuHM0/No2NDKtmqGlW1JkpSVPpwCw96Ix6edWk6jwqqHI9u3h/r1DdRKPyvbkiQpK1W2iwC89PIaAKZ8FN9vrIVESgfDtiRJykoffxzWHQ99KYbtjz6uXInEqrYyg2FbkiRlpS9Xtt97v4y5cwMfVVa2O238M1JNM2xLkqSsU14emDotHm+9dXx95rnA3HnxuJNhWxnCsC1JkrLOnDmwahU0aADHHRtbRh5+JF5r3QoaF9lGosxg2JYkSVmnsoWkYwfYb994vHJlfLWqrUxi2JYkSVmn6kFIaN48YafeVasZd+lsVVuZw7AtSZKyTmVlu9PaVUcGHlBv3TUfjlQmMWxLkqSsU7l5TWWwPmBAVdju4hrbyiDuIClJkrLKwoWBhYsgLy/2bAO0bpXPpRclrCmFli1tI1HmMGxLkqSMtmZN4PEnod+u0H7bhI+nxvNtt4EGDaqC9cEHGbKVeWwjkSRJGe1f/w7cPDjwq98EFi8JX+rXTu+4pM1hZVuSJGWs+fMD/3koHn/+OVx9TaCgYXzvluzKBoZtSZKUse68O7B6NXTYDmbPhuEjIH/tv8u76oiygW0kkiQpI035KPDC0Hh88YUJ550bK9nlFfGcYVvZwLAtSZIyTgiBwbcGQoAB+0P3bglHHg777h2vNy2Gpk1tI1Hms41EkiRllBACQx6FseOgXl046/QYqpMk4aLfQv0Ggb59DNrKDoZtSZKUMVauDFx3feCll+P7Hx8PrVtXBeuiooTLLzVoK3sYtiVJUkb47LPAby8OzJgZH4I868yEHx+X7lFJ341hW5IkZYS77olBu3lzuPKKhF47WsFW9vMBSUmSlHbl5YERI+PxHy43aCt3GLYlSVLaTZoMJUuhUSHs0CPdo5Gqj2FbkiSl3chRAYBddoE6daxqK3cYtiVJUtoNHxFfd+9n0FZuMWxLkqS0WrQoMPnDeNxv1/SORapuhm1JkpRWI0fH1y5doFkzK9vKLYZtSZKUViNGxn7t3axqKwfVeNi+88476dq1K127dmX8+PE1/fWSJCmDlJUFRq2tbO++m1Vt5Z4aDdtTp07lpptuoqCgoCa/VpIkZZjS0ljNnjQZli2DoiLo3i3Ng5JSoMZ2kFcWfkkAACAASURBVCwvL+eiiy5i++23p3379jz11FM19dWSJCmDXH1tBc89D0VFgbprk8iuu0B+vpVt5Z4aq2zfeeedTJ48mauvvpr8/Pya+lpJkpRBFiwIvPBiPF62DBYtjsf77G3QVm6qkcr2lClTuOWWW/j5z39O586da+IrJUlSmq1cGVtFCgqqgvSLL0FFBfTcAS68IGHBAgghVralXJTysF1WVsbFF19Mx44dOfPMM1P9dZIkKc3mzQ88+J/A08/AVsXwz7uhcVFCCIHnX4gB/JCDE7Zrn7Bd+7QOVUq5lIft2267jQ8//JCHH36YunXrbvHnmzRpQl7eht0uxcXF1TE8pYjzk9mcn8zm/GQu5+brVVQE/vb3lfz7wVWUlcVz8+fDfx6qz8W/LWTChFJmzlpKwwZw9JHFFBZWb+uI85PZauv8pDRsT548mdtuu41TTz2VHj16fKufUVJSssG54uJiFi9e/F2HpxRxfjKb85PZnJ/M5dx8s2H/C/zzX7FyvVNv6Ldrwm13BB78zyoOPGA1Qx6N1/beG9asWcKaNdX33c5PZsvV+dmcv0CkNGxfdNFFtG3blvPOOy+VXyNJkjLAkEdimD7hePjF2fFfpSe8F3h7ONxwc2DSpHjfwQf6MKRqj5RXtgF69uy50evHHXccAIMHD2bAgAGpHIokSUqhD6cExr8L+flw7DFVYfrcXySMGh0YPSa+33rrWPWWaouUhu0f/vCHGz0/ZswYZsyYwX777UfTpk1p06ZNKochSZJSrLJFZL/vQ4sWVWG7XduEY44OPPRwfH/QDyAvz8q2ao+Uhu0///nPGz1/8cUXM2PGDM466yx69/avt5IkZbPPFwZefiUef7mqXemUnyS8/HJg2TI4yBYS1TI1toOkJEnKDbNnB/7450D79jFcD3sjUFYW187u3m3DMF1UlHDX7bBiJbT5nmFbtYthW5IkbZGHHwlMnAQTJ8FzzwcqN4Y+9oebDtItWiS0qKHxSZmkxrZr/7JrrrmGDz/80BYSSZKyTHl54PVh8XjHnpCfB+Xl8cHHvfdM79ikTGRlW5IkbbYJ78GixdCoEdx4fcKixfD6MOizE9SpY4uI9FWGbUmStNlefX3txjR7Qd26CVu3hOOOTfOgpAyWljYSSZKUfcrLA8PWtpB8f1+r2NLmMGxLkqTN8u6E2EJSVAR9+6R7NFJ2MGxLkqTN8tqw9VtIJH0zw7YkSfpG67WQ7GPQljaXD0hKkqSNKisLfDAR5s6FKR+HqhaSndM9Mil7GLYlSdJG3XJr4JFH1z+3794u8SdtCcO2JEnawNJlgaefice9e0HrVtCqFRx9pEFb2hKGbUmStIFnnoXVq6FjB7j5hoQkMWRL34YPSEqSpPWUlwcefyKuPPLDow3a0ndh2JYkqZYYPSbw3vvhG+8bPgI+mwuNG8MBA2pgYFIOM2xLklQLfPZZ4DcXBn55fmDuvPUD91+vr+Ckn1Yw7I14/pHH4uuhh0CDBla1pe/CsC1JUi3w5ttQUQGlpfCvf1eF7dFjAk88BTNmwmVXBH5zYQVjxkJeHhx9hEFb+q4M25Ik1QJvvV0VsJ99DuZ8GigvD9xyazzfpTPk58HIUfGePfeAVq0M29J3ZdiWJCnHrVgRGP9uPO6wHZSXw33/CrwwFKZOhUaFcP11CXfentClM9SvDyedYNCWqoNL/0mSlONGjoayMmjXFi76bcJZv4hB++0Rsap98k8SttoqYaut4O47YM0aqF/fsC1VByvbkiTluLfeiqF6j/7Qo3tC/91i//aSJXGzmmOOqro3SRKDtlSNDNuSJOWwsrLA8JHxeI/+MUSfdmpVmD77TMO1lEq2kUiSlMM+mAhLl0JREezQI57r2iXh17+CkhLY7/vpHZ+U6wzbkiTlsMpVSHbvB3XqVFWwjz7SarZUE2wjkSQph731dnzdYw/DtZQOhm1JknLUiJGBmbMgPx/67ZLu0Ui1k2FbkqQcNO6dwGVXxBaSHwyERo2sbEvpYNiWJClLPPpY4MSTK3jjf+Fr73vv/cBFlwRWr4b+u8Fvf23QltLFsC1JUhYYPSZww82xLeSyKwIPP7Jh4J4+IzD41gp+/dvAF6tgl77wpysT6tY1bEvp4mokkiRluHnzA1f+KRACbNMGZs+Bm24JfDI70H7bhJkzAxMnweQPqz7TZycYdJVraEvpZtiWJCmDlZYGfn9lYEkJdOkMt96SMORRuO2OwONPAFRVuPPzoP/ucPBBCbvvtv5Sf5LSw7AtSVKGqqgIXH9j4P0PoFEh/OkPsVJ90gnwve/BkEcCjRtD+22hffuEfrtA06YGbCmTGLYlScpAZWWBQdcGXnwpvr/skoQ2baqC9H77Juy3r8FaynSGbUmSMsCsTwKzZ0OLFlC8Ffzlb4G3h8fWkMsuSdhrT4O1lI0M25IkpVFZWeC++wP/uh/KK9a/Vq8eXHVlQv/dDdpStjJsS5KUJnPmBP7458AHE+P7dm1h6TJYsgSaFsdl+3rtaNCWsplhW5KkNFiyJHDmzwMlS+PDj7/5dcIB+8dgvWZNID8f8vMN2lK2M2xLkpQGr7wKJUthm23ghr8mtGpVFazr1TNkS7nCHSQlSUqDF1+K62MffcT6QVtSbjFsS5JUwz6ZHXd8zM+DAfunezSSUsmwLUlSDRu6tqrdt6+b0Ei5zrAtSVINCiEwdO1GNT84wKAt5TrDtiRJNeiDiTDnU2jYAPbaM92jkZRqhm1JkmpQ5YORe+8NDRta2ZZynWFbkqQaUloaePXVeGwLiVQ7GLYlSfqW3no7cPFlFYx7J2zW/ffdHzexadYUdu6T4sFJyghuaiNJ0rdQURG4/sbAvHnw5luBIw4L/OLshMLCjVesh/0v8M9/xeNfnJ24O6RUS1jZliTpW3h3AsybB3XWlq2efBp+ckpg9JgNq9wzZgauujqeP/YY+MFAg7ZUWxi2JUn6FirXyj5wINz094Q234P5C+D8CwI3D65g9epAWVkM35f8LvDFF9C7F5zzc4O2VJvYRiJJ0hZavTrw2uvxeOABCX12SrjvHhh8W+DxJ+ChIbG1ZPlyKFka72vZEv70h4Q6dQzbUm1i2JYkaQsNHwHLV8QA3btXPNegQcJvfpWwe7/AoL8E5nwaz2/VJK6nfdIJCcXFBm2ptjFsS5K0hV4cGltIBg6AvLz1A3T/3RPuuxvefBu2aQM79sRqtlSLGbYlSdoCJSWB4SPj8cBNrJXdtGnC4YfW4KAkZSwfkJQkaQu89jqUlUHnTtBhOyvWkr6eYVuSpK+oqAg8/EjgqWfWX8avvDzw5NNrW0jcAVLSZrCNRJKkLwkhcP0NgSeeiu/r16taF/vhR+Cjj6FRYVzyT5K+iWFbkqS1QgjcdEtV0Ab4698D3btDSUk5d94dq9rnnuPKIpI2j2FbkiRi0L7tjsCQR+P7Cy9IGPpSYPy78Ps/BhoVLmfNGtilLxxyUHrHKil7GLYlSQLuvQ8e+E88/s35CYcfmrDbrvCz0wNTpgCU0bABXPibhCSxqi1p8/iApCSp1vv3g4F7/hlbRM47J+GoI2KYbtky4dKLq4L12WcltG5t0Ja0+axsS5JqrRACQx6B2+6IQfusMxKOO3b9ML1H/4SLL4SVKxtw1BGr0jFMSVks5WF73rx5PP/887zxxhtMmzaNzz//nCZNmtCnTx9OP/10evXqleohSJJqidLSwPQZUFEe3xcUQttt2KDtY/GSwNCX4NnnAtOmx3M/+yn85MSNV60PPTihuLiAxYtXp3D0knJRysP2/fffz5133km7du3o378/zZo1Y+bMmbz88su8/PLL/O1vf+Pggw9O9TAkSbXAVVcHXnlt/XN77Qm/OR+aN0tYuTLwz/sDDw+JG9MA1KsLPzkp4ZSTa368knJfysP2jjvuyAMPPEDfvn3XOz9mzBhOOeUUrrzySgYMGEC9evVSPRRJUg5bvCTw+rB43LIFkMDChfC/N+Gd8YGjjww8/wIs+Dzes31XOOTghP33g8ZF9mFLSo2Uh+2BAze+6n/fvn3p168fb775Jh9++CE9e/ZM9VAkSTnstdehvAK6dIF77ojP/0+dFhh0bWDyh/Cvf8f7vvc9+L9zE/bob8CWlHppfUCyTp06671KkvRtvfTy2m3UB1SF6I4dEm4bHHd+fObZwID9E078MdSvb9CWVDPSlnI//fRT3n77bVq0aEGXLl3SNQxJUg747LPAe+9DksD+31//Wp06CSccDyccb8CWVPPSErZLS0u58MILWbNmDRdccAH5+fmbvLdJkybk5W24HHhxcXEqh6jvyPnJbM5PZnN+ttwjj30BrGTXXerQpUuTlH2Pc5PZnJ/MVlvnp8bDdkVFBZdeeimjR4/mRz/6EUceeeTX3l9SUrLBueLiYhYvXpyqIeo7cn4ym/OT2ZyfLRdC4MmnYgvJ9/cpT9nvz7nJbM5PZsvV+dmcv0DU6A6SIQR+97vf8dRTT3H44Ydz5ZVX1uTXS5Jy0MdTYfoMqFsX9tk73aORpPXVWNiurGg/+uijHHrooVxzzTUbbQ+RJGlLvPRKrGrvvhsUuYSfpAxTI2m3oqKCyy67jMcee4yDDz6Yv/zlL1/bpy1J0uZYvTrw4tB4/OVVSCQpU6Q8bH85aB944IFcd911Bm1JUrV45tm4cU3LltB/93SPRpI2lPIHJAcPHsxjjz1GQUEB7du359Zbb93gngEDBtCtW7dUD0WSlEPWrAn8+8HYQvKTExPq1bOyLSnzpDxsz5kzB4CVK1dy2223bfSeNm3aGLYlSZsUQuDfD0KDBnDMUZCXl/DMc3Hr9ZYt4JCD0j1CSdq4lIfta665hmuuuSbVXyNJymET3oPb74xV7HfGw0UXwL8fiO9PPMGqtqTM5T7pkqSM99TTYd3xG/+Dd8YHli2D5s3h0IPTODBJ+gauvSdJymglJYHXXo/H5/8yoVlTWLYsvj/phIT69a1qS8pcVrYlSRntxaGwphQ6d4Kjj4K990q45rpARQUcdki6RydJX8+wLUnKWCEEnnomtpAcflhCkiS0aAF/+4vVbEnZwTYSSVLGmvAezJgZVyEZOCDdo5GkLWfYliRlrMoHIwfsB4WFVrMlZR/DtiQpI02dVvVg5BGHG7QlZSfDtiQprSa8Fzj0yAouvqyCyZNjJfv5FwJn/jywphR6dIftu6Z5kJL0LfmApCQprf7z38CSJfDmW/DmW4GOHQJTp8Vru+4Cl18WH4yUpGxkZVuSlDZLlwaGj4zHe+0J+XkwdRokCZx+asJ11yQUb2XQlpS9rGxLktLm9TegrAw6doRBV+UxZ07guRcCu/RN6N3LkC0p+xm2JUlpM/Sl2KM9cEAM1m3aJJxxmiFbUu6wjUSSlBbz5gfenRCP998vvWORpFQxbEuS0uKVVyEE6N0LWm1tNVtSbjJsS5JqRFlZYMbMQEVFbB156eW1G9bsb9CWlLvs2ZYkpdzq1YELLwmMHQfNmsGufQMffQx16sD390n36CQpdaxsS5JSqqwscOVVMWgDLFwIz78Yj/vtCk2aWNmWlLusbEuSUiaEwF+vD7zxP6hbFwZdlRACvDYsMHUqnHKyQVtSbjNsS5JS5t774JnnIC8P/nB5wm79YrjefTdDtqTawTYSSVJKzJwZuO/++BDkBb9O2GdvA7ak2sewLUlKiVtuDZSXQ//d4fBDDdqSaifDtiSp2o0YGRg+Iq42cu4vDNqSai/DtiTpW3n/g7hu9leVlQVuHhzPH3M0tGtr2JZUexm2JUlb7P0PAj8/N3DKaYHnX1w/cD/6OMycBVs1gVN+YtCWVLu5GokkaYvde18gBCgrgz8PCsydC/t/H26/KzDsjXjP6aclFBUZtiXVbla2JUlb5IOJgZGjID8PDjs0nrv73sAJJ8egnZcHRx4Bhx2S3nFKUiawsi1J2qSyssCYsbB9V9hqq1ilvve+2Dbyg4Fw0QV5dO0S+PsNgfKKuPLI2WckdOhgRVuSwLAtSdqIEAKvDYM77gzMngNbbQWXXhRfR4yMVe2TT4qB+sjDE3p0h/Iy2H57Q7YkfZlhW5IExIA9ew6MGg0vDA1MmhTP5+XBkiVw4SWBZs3iuYEHwDbbVAXrzp0M2ZK0MYZtSRJvDw/8/abAZ59VnWvYAI4/Do45KuH+BwIPDYGFC2P4PtlVRiRpsxi2JamWmzw5cPkfAqtXQ926sGNP6LtzwsEHQrNmMVSfd07Cbv0Ct98Z2HOPhLbbGLYlaXMYtiWpFluwIHDx72LQ3n03+OPvExo23HiQ3qVvwi59DdmStCUM25JUC4UQWLYcLvld4PPPoX17+MPlmw7akqRvx7AtSTlo0aJARYDmzarC89x5gX/cFhg1GlauhIqKeL5JY7j26oTCQoO2JFU3w7Yk5Zh58wM/PTWwYgXs3Cdw8EEJc+fCfffHdpEva9Ysto60+Z5BW5JSwbAtSTlm8K2B5cvj8ZixMGZsWHet145w9pkJrVtDYQE0aABJYtCWpFQxbEtSDhk7LvDqa3F5vquvSvjww8CLL0FeAqeeknDAAMO1JNUkw7YkZanhIwL3PxDotj2ccnJCw4bw95tiFfvIw2HP/gl79k847WdpHqgk1WKGbUnKMp99FrjplsD/3orvJ7wHQ18O9NkJZsyArZrA6adavZakTGDYlqQM9/4HgcG3xiX6VqyE5cvjSiL5+XDYITDuHZj1Cbzyarz/rDMSGjc2bEtSJjBsS1IGW7QocOnvAosWr39+5z7wq18mbNc+obQ0bqV+3/2BHt3hkIPTM1ZJ0oYM25JUg8rKAqPHQL160Genr39YsaIicNWgGLQ7bAcX/TahsBCKGlVtow5Qt27CSSfAj4+L7/PyrGpLUqYwbEtSDZgxM/D0s4EXh8KSJfFcn53g/P+D7drHcFxeHlixomqZvgf/C6NGQ/36cOXvk3X3bUp+viFbkjKNYVuSUmjuvMBd98SQHdbm6KbFsHxF7LU+5bRAv10D8+fDrFmwpnQRzZpCu3YwYUK8//xffnPQliRlJsO2JKVAaWkM2UMegTWl8dyee8ChhyTstivMXwA3r11R5O3h63924aL4B2D//ezBlqRsZtiWpO9g5szA/Q8G+u+e8P19Yg/20mWByy4PvDM+3rNTb/jF2Qndtq+qTn+vNQz6c8LYcYGPPoa220D7baFdu6147/0lzJgRVx057FA3oZGkbGbYllTrLV8eCAGKirYs1H4yO/DL8wMLF8ELLwZ27gMnnQA33hyYMRMKCuCySxL23nPTgXnnPgk796l637hxHt27JXTv9l3+iyRJmcKwLSlnrVoVtyrv2wfatNl42F24MHDqGYGly2D//QLHHpPQtcuG906fERg+AvruDJ07wWdz4f/WBu3WrWLbx9hxcbt0gBbN4S/XJHTuZFVakmozw7aknFRaGrj08sCo0VBUBNddAzv0WD/4hhD42w1hXX/0Cy/GCnXvXoGzz0zW3f/s84G//T2wZk2873vfg7LS2Hfdflu4+YaElV/ATbcE3nobOnaE6wYltGxp0Jak2s6wLSnnlJcHrro6Bm2AZcvgV78J/PmP0G/XqgD82jB4439xJ8bLLkkYMSLwymsw/l04+5zA3nsFCgvg+Rfj/e3bw6efxj8A27SBG/6WUFycUFwM116dMOuTQOtWce1rSZIM25JySgiBG26OoblOnbg+9VNPB0aOgosuDZx5OhxyUFyG7/obYsvHySfBwAEJAwcknH1m4N77As8+H4M4QF4enPazhJ+cCKtWwYhR8NFHgaOPTGjefP1Q3a6tIVuSVCUJIYRvvi19Fi9evMG54uLijZ5XZnB+Mls2zU9ZWeCOuwILF8JJJybrbf7y+jB4Z3ygfgMoLEgIITD5Q5g4EUqWQpLAHy5P2H+/uJ35VYMCr7waf27durB1S5g9J+7MePcdyQaV6Okz4nd/8gn8+lcJfXaqmRCdTfNT2zg3mc35yWy5Oj/FxcXfeI+VbUk1IoTAokUwcxZMnxF3VJw9Oz5I2L1bQo/u0KFD1S6Iq1cHrrgy9kADvPRK4OADAzv0SHjwv4FZn6z309f7rvr14f/Oi0EbYkvHFZdB717w1NNxqb3Zc2LF+tKLNgzaEHd1HHSVVWpJ0ndj2JZULca/Gxg5KtCpUwzOzZrC+x/AqDGBd9+NAXvZso1/9rkXYlhu1gwOHBj4/r4Jt/wjMP5dqFcvhuRRo+GZ5+CZ5+K9RUVw8IGQXwdWroCyMtZ9d6eOG/ZM5+cnHHUEHHVEwkcfBV55LdC5c8L22xuoJUmpY9iWBMCKFYHHnoAe3dlky8RnnwUeeiRAgFNPSWjcON738iuBP10dKC+Hyipzfj5r31fJy4ububRvD9u2g7bbJHw2NzBxEnwwERYuhAf+Aw/8J/6MwsL40GHvXgkT3gvcfmdg9hz44dEJRx8JhYXfLih37pzQubMhW5KUeoZtSSxeErjgwsCHU+L7XXcJ/PyshE4dYcUKmPMpPPZE4IUXqwL0q68Ffv0rWLwkPmgYAuzYE1athqkfx/uKi2GXnaHvzgmdO0O7tlC//ldDbnxfWhpbRp57PjBiFGzVBP76l4Qua0Pxjj0TBt9kQJYkZRfDtpRlKh8EbNCAdQ8Mfhdz5wXOvyA+CNioUVxtY9RoGD0m0KABfPHF+vf33Rk+/xxmzITf/b6qV/qoI+H8Xybk5SWsWlW12Ute3uaNsW7dhH33gX33SSgpCdStCwUFhmtJUnYzbEvVbOmyQP16G6vgRmVlganTYOnSWAne1H0b8/4HgX/cFpjwXny/U++442GvHeHzhXGTlbp1oOcOVT83hMCsWfDxNAgVUFi4miUlcYWPBQsCb/wPFnwOW28Nf/9rQn4e3HFXXDqvMmg3bhzHetIJcaOX1asD990feOBBKK+AU06OS+NVbkneoEFCm+99+99hkyaGbElSbqixpf8mTJjAzTffzPjx4yktLaVTp0789Kc/5bDDDvvaz7n0X/bJ5vmZOy8w7A349NPA/vsl7Nhz80LfvPnxc68PC7z3PhQUxGB67DExeM6bH3j+BRg1OlalK3ciLCqCA/aHvfdKWLQYZswIzJsf12ru0X3ttuCfwQeTYOSowNvD4+fq1YPyshh0N6ZeXejVC5o3g7HvwPz5Xz/+9tvC9detv+Ph3LmB0rK4WkiDBhv/PcyaFVi0GHr3MhxXl2z+/yfXOTeZzfnJbLk6P5uz9F+NhO2RI0dy2mmnUbduXQ455BCKiooYOnQos2fP5vzzz+fss8/e5GcN29kn2+ZnzZrAcy/As88FJk1e/9rOfWJobt1qw8+tXgMjR8WAPXHSxn92i+aw7bYwdlzcRKVSo0ZQvx7rtgnfXHl5cQWO036WUBHg8ScCTz0Tq+RNGkOLFlBSEivVX1avLnTpEkN6nTp1qKgoo3nzeH/rVgn77QuNGhmYM0G2/f9Tmzg3mc35yWy5Oj8ZEbbLyso46KCDmDt3Lg899BDdu3cHYPny5Rx//PFMnz6dZ599lvbt22/084bt7PPV+Vm9OlZzy8tjxbewMG4oUq/e14e7srLAO+Nh2BuBKR9Bw4bx87GlIWGXnaFly7iZybz5MHt2DKOFhfG+srL4cN/KlfF1xcq4RFydujEEt2gBH3wA/34wrAunSQK9doSWLeCV1zZcTWNTkiS2buy7d8Jee8KE9+DOuwNz51Xd07sX/OCA2PKxzTYxfI8ZGx8IfO8DaLV1XKWjZYuE6dMDH0yEufNi9bt7t/hnv+8nG/Rpl5cHysrWbxuZOQtGj4HFiwO9dozfWVmd9v+fzOb8ZC7nJrM5P5ktV+cnIza1GTFiBLNmzeLoo49eF7QBGjVqxC9+8QvOP/98HnvsMX7961+neig5r7w88MUXUFHBuiXZvqvVq2O/7/SZsOhLVdi8vLXBuaAq3BYWxsrpzFmlTJ0WmDMn7vA3/t2qtolK9etD716BXXdJqFsXPpgYw+WSJfHnFBbEZeBKlm58XM+uXWu5davAkpINH+LbUs2bw/E/Sjhgf2jWLP7uzjw98O//BF57DUrLNvxMkkDXLvGBvr33gubNqn7nrVvDvvvA8y/GtaX32xfatNlwTvrtCv123dhcxXMrVgQaNvz6hwzz8xPy8788roT228bWkMqfI0mS0iPlYXvUqP9v796Doqr7P4C/l7vrcveOl6x1EXSs8AFRaLyP+XgpNKl0kBll1LykzFhjk009PE7aHxmDxQxWOkahJthYXkrNKwiuYj4F8tNqENhFAUWRBWV34fz+WHd13SVu57C4+37NOML3e/bs5/Bhl89+z/d8jxoAEBsba9MXExNjtU1PYTAIOHUGaGoyFX1yOeDx8CclCKZ28ygpZI+KQ5mbaRS1scF0sdr16wKul5lO8Y9QAuHhMiiVptHSxocjrm2dVhAE4N49ATU1povf6upMz32/0TSNwazZCNx/8Oj7fv0EjAoHQlUy+Pi0/RwPHgCNjQIaGoDaO0BNjWkqQk2N9fSH9rGtkIODAV+F6Zjr602xnleb5iE/Sad79HWAP/BSLPCvf8nQ0mIaoa6qFlBYCPzfVeDGTdN27u5AyCDThwDzKLanByDv/egDQW850EsOGAyw/DzlciD+NRn+/bLthYoDBsiwPlmG9ckdPX4Tb28ZXp3buceadXYdaSIiIuoZJC+2r1+/DgAYNmyYTZ+/vz8CAwNRVlYmdRgdci4f+M9/xZ1dU3MLOFfQLdeiAjBdEFddDZw81fXn9PMzjZL26weYB1ibm4HG+49Nz3j4IeNBExAc7IbgoBb06wuMCpchMtL0ePNKFYIgoLQUUF8ELhaa1mcODzNtO2CAaZS6ocE0Sh4eBnh42FmXOcm06se1a6YpISEh9rYjIiIicizJi23dw2FKX19fu/0KhQI3b95s9fH+/v5wc3OzaW/PHJnOmjSxBa/HN+LGjRY0969wJQAADh9JREFUNAjQNQhoeWzurrc3oFC4obfc9L3usW0UChl695YhMMANzz7rjueedYevrwxXrhjxvz+MKC1tho+PDAqFDHK5DO1ZgtjXV4b+/d3Qv78bggLdLM/h7S2zTBJw9wAUvU3tBiNQXGzE/3434s+/jFaxt8bHx/TY3r1lCAyUYUB/N/Tv54ZBIe4IDnq0pJtYgoKAsWO7to/AQNNdCKnjpHz9UNcxPz0Xc9OzMT89m6vmp8evs11XV2fT1h2T7NesbGuL1ipY4bF+g6V1+DPArH8/uY34GhpM/49Qmv6J4e7djm3vrBdBOAvmp2djfnou5qZnY356NmfNT3s+QNgOGYtMoVAAAOrr6+3263S6Vke9iYiIiIieZpIX2+Yl/ezNy66rq8OdO3fszucmIiIiInraSV5sR0ZGAgByc3Nt+vLy8gAAUVFRUodBRERERNTtJC+2x48fjyFDhuDgwYMoKXl0mz2dTof09HR4eHggLi5O6jCIiIiIiLqd5BdIenh4YNOmTUhKSsLChQsxe/ZsKBQKy+3a161bh+HDh0sdBhERERFRt+uW1Uiio6ORlZWFtLQ0HDlyBAaDAUqlEmvXrsXcuV286wcRERERUQ/VbUv/jRkzBl999VV3PR0RERERkcNJPmebiIiIiMhVsdgmIiIiIpIIi20iIiIiIomw2CYiIiIikgiLbSIiIiIiibDYJiIiIiKSCIttIiIiIiKJsNgmIiIiIpIIi20iIiIiIomw2CYiIiIikgiLbSIiIiIiicgEQRAcHQQRERERkTPiyDYRERERkURYbBMRERERSYTFNhERERGRRFhsExERERFJhMU2EREREZFEWGwTEREREUnEw9EBAMDvv/+Obdu24fLlyzAYDFAqlUhMTMScOXM6tB+dTocdO3bg6NGjqKiogKenJ4YMGYKpU6di9erVEkXv/Lqan4SEBKjV6n/c5pNPPsGrr74qRrguRYzXzr1797Bz504cP34cGo0GXl5eGDx4MOLi4rBgwQJ4e3tLeATOTYz83Lx5E+np6Thz5gxu3bqFgIAAvPTSS3j77bcxcOBACaN3bgcOHEBhYSGKiopw7do1GAwGbN68GfPmzevQflpaWpCVlYW9e/eirKwMcrkc48aNQ3JyMp555hlpgncBYuTn9u3byM7ORnFxMYqKiqDVagEAV69elSpslyFGfi5evIjjx49DrVZDq9WisbERISEhmDp1KpYvXw4/Pz8Jj6B7OXyd7fPnz2Pp0qXw9PTErFmz4Ovri6NHj0Kj0SA5ORkrVqxo134qKyuRmJiIiooKTJgwAWFhYdDr9SgvL0dlZSV++ukniY/EOYmRn/3791ve5B5nNBqRkZEBNzc3nDx5Ev3795fiEJyWGLm5d+8e5s2bh4qKCowdOxbPP/889Ho9zpw5g/LyckRHR2Pnzp1wc+NJsI4SIz/l5eV44403cPv2bcTExCA0NBRlZWU4ceIEgoKCsGfPHgwdOrQbjsb5TJkyBVqtFoGBgZDL5dBqtZ0qtj/44AN8//33UCqVmDhxIm7fvo3Dhw/D29sbe/bsgVKplOgInJsY+Tl//jwWL14MmUyGYcOGoaqqCvfv32exLQIx8hMTE4M7d+5g7NixCAsLg0wmg1qtxpUrVzB06FDs2bMHwcHBEh5FNxIcyGAwCNOmTRNGjx4tFBcXW9rr6+uFWbNmCeHh4UJpaWmb+zEajcL8+fOFMWPGCPn5+XafhzpOrPy05ueffxZUKpWwfPlyEaJ1LWLlZvv27YJKpRI+/vhjq/ampiZh/vz5gkqlEtRqtdjhOz2x8rNs2TJBpVIJu3btsmo/fPiwoFKphCVLlogdusvIy8sTNBqNIAiCkJGRIahUKiEnJ6dD+8jPzxdUKpWwcOFCoampydJ+7tw5ITQ0VFi0aJGoMbsSMfJTU1MjqNVqob6+XhAEQZgxY4agUqlEj9UViZGfjIwMoaqqyqqtpaVF+PDDDwWVSiV89NFHosXraA4driooKEB5eTlmz56N8PBwS7tCocDKlSthNBqxf//+Nvfzyy+/4I8//sCSJUsQHR1t0+/h0SNmyzx1xMpPa/bt2wcAeO2117ocq6sRKzcVFRUAgIkTJ1q1e3l5ISYmBoDpVCx1jBj5aWpqQm5uLvr06YOEhASrvpkzZyIsLAy5ubmWHFLHTJgwASEhIV3ah/k9bN26dfDy8rK0jx8/HrGxsbhw4QJKS0u79ByuSoz89OnTB5GRkVAoFCJFRWZi5GfZsmXo16+fVZtMJsPKlSsBABcuXOjS/nsShxbb5nm8sbGxNn3mP/RtzfUFgMOHDwMAXn75Zdy4cQO7d+/G9u3bceTIETQ0NIgYsWsRKz/23Lx5E3l5eejbty8mTZrU6RhdlVi5GTFiBADg7NmzVu0GgwHnzp2Dj48PXnzxxa6G63LEyM/du3dhNBoxaNAgyGQym/7BgwcDMBX25Bjnz5+HXC5HRESETZ85985UMBB1B/MAqbu7u4MjEY9Dh3yvX78OABg2bJhNn7+/PwIDA1FWVtbmfoqKigAAhYWF2Lx5M/R6vaUvKCgIqampGDdunDhBuxCx8mNPTk4OWlpaEBcXxzMPnSBWbhYsWIADBw5gx44dKCoqwujRo2EwGHD27FnU1dXh008/5Vz6ThAjP35+fnB3d0dlZSUEQbApuDUajdVzUfdqbGxETU0NVCqV3aLAfHEk80PUMTk5OQAeDUw4A4eObOt0OgCAr6+v3X6FQoH6+vo292M+zb1p0yYkJibi9OnTyM/Px8aNG1FfX49Vq1ahurpavMBdhFj5eZIgCJZT6JxC0jli5cbHxweZmZmYO3cu1Go1duzYgczMTMsUCHsjdtQ2MfLTq1cvREZG4tatW8jKyrLqO3r0KEpKSgCgU69B6jrzz721KQrmdvPvAhG1raSkBF988QWCg4ORlJTk6HBE4xRDisLDBVUmTZqE9evXW9oTEhJQVVWFL7/8EtnZ2ZZ5QORYBQUF0Gg0iIqKsjvyR92ntrYWK1euRG1tLbZv346IiAg0NTXhxIkT2LJlC06dOoWcnBz4+/s7OlSX9N577+HNN99ESkoKTpw4gdDQUJSXl+PXX39FaGgorl69ypViiMgpVFRUYPny5WhubsbWrVsRFBTk6JBE49B3afMn/9ZGZnQ6XasjQ/b2M2XKFJu+yZMnA3g01YTaT6z8PIkXRnadWLnZsmULfvvtN6SlpWHixInw9fVFnz59EB8fj3feeQcVFRXYtWuXqLG7ArHyM3LkSGRnZ2PmzJm4cuUKvvnmG5SWliIlJQWvvPIKADjVH6SniTl/rY1cm9t5cR5R27RaLRITE1FbW4u0tDS7i108zRxabJvntNmbu1hXV4c7d+60a+Rz+PDhAGB3AXRzW1NTUxcidU1i5efJxx07dgx+fn6YMWOGGGG6JLFyc/r0aQQEBGDkyJE2feY3u+Li4q4F64LEfO0899xzSE1NRX5+PoqKinDo0CEsWLAAf/75JwBg9OjRosVN7SeXy9G3b19oNBo0Nzfb9JvnavPGNkT/TKPRICEhAdXV1UhNTbUMkjoThxbbkZGRAIDc3Fybvry8PABAVFRUm/sxFwV//fWXTZ+5ratL1LgisfLzuB9//BF6vR5z5syBj49P14N0UWLlRq/XQ6fTWV1UbFZbWwsAVkuaUftI8dp5nE6nw8mTJxEQEOBUFxE9baKiotDY2IhLly7Z9Jlzb/5dICJbGo0GixcvRnV1NT777DNMmzbN0SFJwqHF9vjx4zFkyBAcPHjQcrEPYPpDkp6eDg8PD8TFxVnaa2tr8ffff1uKALN58+bBy8sL3377Laqqqqz2k5GRAcC0Li11jFj5eVx2djYATiHpKrFyExERAaPRiPT0dKt2vV5vaeNKPh0nVn4ePHgAo9Fo1abX6/H+++/j7t27WLVqFby9vaU9GGo1P/Hx8QCA1NRUqw+s+fn5yM3NRWRkpOXMK0mnPX97yHFay4+50K6qqsLWrVsxffp0B0UoPYffrr2goABJSUnw9PTE7NmzoVAoLLc0XrduHd566y3Lttu2bcPnn3+O1atXY82aNVb7yczMxKZNmxAQEIDp06fDy8sLp06dglarxeuvv46UlJTuPjSnIFZ+ANO8+fnz52PUqFFduhkOmYiRm5KSEixatAgNDQ0YM2aM5QJJ881SRo0ahd27d7Og6wQx8nPx4kWsWbMGEyZMwMCBA6HT6XD69GlUVlYiPj4eKSkpdtfgprbt27cPhYWFAIBr166huLgYERERluk906ZNs4yy/dN728aNG7Fv3z7erl1kYuVnw4YNlq+PHTsGnU5n9UH33Xff5XUPnSBGfsy3fH/hhRfs3pMAgN1a4mnk8NVIoqOjkZWVhbS0NBw5cgQGgwFKpRJr167F3Llz272fhIQEhISE4Ouvv8ahQ4fQ3NwMpVKJFStWWEYfqOPEyg/AUW2xiZGbsLAw7N+/HxkZGSgoKMB3330Hd3d3DB06FGvWrMHSpUtZaHeSGPkZNGgQoqKiUFhYiFu3bqFXr14IDw/Hhg0beM1DFxUWFuKHH36wart06ZJlSkhISEi7TmmnpKQgNDQUe/fuRWZmJuRyOSZPnozk5GSOaneBWPl5ch9Ptq1evZrFdieIkR+tVgsAuHz5Mi5fvmx3G2cpth0+sk1ERERE5Ky4QCsRERERkURYbBMRERERSYTFNhERERGRRFhsExERERFJhMU2EREREZFEWGwTEREREUmExTYRERERkURYbBMRERERSYTFNhERERGRRFhsExERERFJhMU2EREREZFEWGwTEREREUnk/wGIJ0mHpBUAhQAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 720x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "H = vmap(lambda i: sim_p(i, G=1.4))(jnp.arange(int(1e5)))\n",
    "az.plot_kde(H[\"H\"], bw=0.25)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Code 10.11"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "entropies = H[\"H\"]\n",
    "distributions = H[\"p\"]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Code 10.12"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "DeviceArray(1.2217282, dtype=float32)"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "jnp.max(entropies)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Code 10.13"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "DeviceArray([0.09018064, 0.20994425, 0.20969447, 0.49018064], dtype=float32)"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "distributions[jnp.argmax(entropies)]"
   ]
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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.11.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}