{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "header",
   "metadata": {},
   "source": [
    "# Categorical Sampling\n",
    "\n",
    "Lecture 9 | CMU ANLP Fall 2025 | Instructor: Sean Welleck\n",
    "\n",
    "Sampling from a categorical distribution over vocabulary tokens."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "setup",
   "metadata": {},
   "source": [
    "## Setup\n",
    "\n",
    "Given a probability distribution over tokens, compute the cumulative distribution function (CDF)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "cf48e08b",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([0.1, 0.3, 0.4, 0.8, 1. ])"
      ]
     },
     "execution_count": 1,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "import numpy as np\n",
    "\n",
    "vocab = ['a', 'b', 'c', 'd', 'e']\n",
    "probs = np.array(\n",
    "        [0.1, 0.2, 0.1, 0.4, 0.2]\n",
    ")\n",
    "cdf = probs.cumsum()\n",
    "cdf"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "ef7db492",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.19169709790148137"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.random.random()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "linear",
   "metadata": {},
   "source": [
    "## Linear search\n",
    "\n",
    "Sample a uniform random number and find the first CDF value that exceeds it."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "1d921cec",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "'d'"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "def sample_linear(cdf):\n",
    "    r = np.random.random()\n",
    "    for i in range(len(cdf)):\n",
    "        if cdf[i] >= r:\n",
    "            return i\n",
    "\n",
    "idx = sample_linear(cdf)\n",
    "next_token = vocab[idx]\n",
    "next_token"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "binary",
   "metadata": {},
   "source": [
    "## Binary search\n",
    "\n",
    "More efficient O(log n) sampling for large vocabularies."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "9c323319",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Sampled token: c\n"
     ]
    }
   ],
   "source": [
    "def sample_binary_search(cdf):\n",
    "    r = np.random.random()\n",
    "    left, right = 0, len(cdf) - 1\n",
    "    \n",
    "    while left < right:\n",
    "        mid = (left + right) // 2\n",
    "        if cdf[mid] < r:\n",
    "            left = mid + 1\n",
    "        else:\n",
    "            right = mid\n",
    "    \n",
    "    return left\n",
    "\n",
    "# Test binary search sampling\n",
    "idx = sample_binary_search(cdf)\n",
    "next_token = vocab[idx]\n",
    "print(f\"Sampled token: {next_token}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "empirical",
   "metadata": {},
   "source": [
    "## Empirical distribution\n",
    "\n",
    "Sampling approximates the true distribution."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "ec2mk261uis",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAABKUAAAGGCAYAAACqvTJ0AAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjYsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvq6yFwwAAAAlwSFlzAAAPYQAAD2EBqD+naQAAPDFJREFUeJzt3QuYVVXdP/AFKCAqoCLgHW+pZIJKIL4aXsNUzGvGXwVR6eI1UVM0QUXDvPCiaeINy0ylzDK6oIWalRQGiYp30/AOpIKigsL8n9/uPdPMMAPDXPY5zPl8nuc4c/bsc846Z8+4v/zW2mu1qqioqEgAAAAAkKPWeb4YAAAAAARFKQAAAABypygFAAAAQO4UpQAAAADInaIUAAAAALlTlAIAAAAgd4pSAAAAAOROUQoAAACA3ClKAQAAAJA7RSmgRerRo0c6/vjjm/11XnnlldSqVav0wx/+sHJbvO4666yT8hKvf9FFF+X2egCwOnn44Yezc2V8rY+99toru61uGaW2TNJc4jXiteI1q76vgw8+OJXiMW1qy5YtSzvuuGO67LLLUrmYMGFC2nzzzdPixYuL3RRaGEUpKGNxMq/PrVgn/IIIhoW2tG7dOnXs2DFtt9126bjjjku///3vm+x1fvvb35ZscaeU2wYANYsVdd3++te/FruJq42qn9saa6yR1l9//bTrrrumM844Iz399NNN9jo/+MEPcilktaS23XXXXenVV19Np556aq6vO2nSpHTsscembbfdNvu9WFHxNIpH5557btp4443TWmutlfr161dnbn700UfTHnvskTp06JC6d++eTj/99PTBBx9U2ycKqUuWLEk33nhjk78vylurioqKimI3AiiOO+64o9r922+/PTtZ/fjHP662ff/990/dunVLxRIn3JdeeimNHTs2u79o0aL04osvpnvvvTf985//TF/5yley97LmmmtWOxFHAavqtpWJYHH99denVfnfYuwbrxWv06ZNm8qT9j333LPcybwxVtS2jz/+OAurcQOAYooCwrBhw9Ill1ySttxyy+V+fsABB6QuXbrkPqol/jHdtm3bLBusTOwbYv/mECOKItusrNgSRYfIYEOGDMnO/wsWLEizZs1KP/vZz7Is9L3vfS+NGDFihZmkPmLETxyTVemEXLp0afrkk09Su3btsnYW3lc8169//et6P09D27aqx7Sp9e7dOyvy5F2gid+bGTNmpM9//vPp8ccfTzvttFOdx23w4MFZHv3Wt76VFbHi9+2xxx5LDz30UFaAKojn6d+/f9phhx3S1772tfTaa6+lq666Ku29997pd7/7XbXnjCJXFMZefvnlyuMOjeVfMFDGoqelqui9jKJUze01ffjhh1lPSp46deq0XLsuv/zyrCcnetEiCEU4K4iQ1Jw+/fTTLBBFGGrfvn0qpmK/PgDU9KUvfSn16dMnlYIoWtTnXFnIN81VjGqIz3zmM7Xmn0GDBqWzzjorbb/99unAAw/MtkeRoLkzQRTD1l577azotSqFr2Id0+bwj3/8IysOXn311bm/dnQcb7LJJtn7j4JdXaZPn57uvvvudOWVV6azzz472xbFzXjMt7/97WxkVMH555+f1ltvvay4FVcjhMjVw4cPTw888ED64he/WLlvdARfccUVWWFrn332adb3Svlw+R6w0h6ZOIFFr8wXvvCFLKzFyWtFcxnVNlfCe++9l/XUbLbZZlnBaJtttsmKSFHYaagIQ9dee23q2bNnuu6667IexLraEL15F198cdZTFCFmgw02yHqJCsOYY98YiVR4X4Vb1Tkaotdo/Pjxaeutt87eQwydX9H8DTGKa+DAgVl4i6HT0WtcdaRTXfMh1HzOFbWtruMQgSn+QRDhIua32nfffZe7ZKJwicVf/vKXrKd1ww03zNp62GGHpXnz5jXomABAfVQ9t8Y5bquttsoyRvwDOC6LivPlmDFj0qabbppdevTlL385vfPOO9WeozCHUfzDOUauxPk9MkGMpK6qtvPtivJNbXNKxajkONdGkSheZ6ONNkqHH354NpK7IN7L7rvvnmWMaHNcahcjVZpaPH8UHGKEdNU5jWrLJG+99VY2ci0+x8gu0e74LAtzQcVnOHv27PTHP/6xMl8U3nshJ8TPTj755NS1a9fseeqaU6pgZccjPsfaRtnUfM4Vta2uDBWjyOJzj88/RlhFQe/111+vtk9h7s/Yfuihh2bfRwaK4k2MAFuZX/7yl1nhMn5vantfMZo/XqNz585Zp2p8/lHwbAqRo+szMix+7yInx8ingjgeJ554Ypo2bVr2NxYWLlxY2SFdKEgVCljxufz0pz+t9rzx2cZlpPfdd1+TvB8IRkoBK/Xvf/87K3B89atfzU5aq3opX5yIBwwYkJ38v/71r2eTJEYPzciRI9Obb76ZFXoaKk64MTz5wgsvTH/+85/TQQcdVOt+ERTi8r+TTjop9e3bNzsJ//3vf08zZ87MhsZHu954441aL18suO2227JQGif4CHZxUq6rqBahJi5R2G233bIepSlTpqTRo0dnI6yiOLUq6tO2qiLA7bnnnlm4iN6wGMYfw8sjyEWwi+HmVZ122mlZD1m0L4JgHI+4XDCGZwNAQ0RH0fz586tti3+wR0Glqp/85CfZZVhxLoqiU5wzYzRGjMKIgkNcLhT/yP/+97+fFQ0mTpxY7fEvvPBCOvroo9M3vvGNNHTo0OxcfdRRR2Xn3Ti/N0W+iXN6FL+mTp2a7RtzOr3//vvZefmpp57KOqvCNddckw455JB0zDHHZO8pCkfRlricra580lCRpSJbxYiVyDRVCwpVHXHEEVkuiM83ijxz587N2j1nzpzsfpzz42dRgLjggguyx9T8HKIgFUWbUaNGZSOlVqQxx6Om+rSttktH49K2yHxvv/12dkyi8y0666JIVPWYRsdhZKIoJv7hD3/IRj7FsfzmN7+5wnZFho2CZl1TRMTvb1y6Gm2InHnLLbdkBb2qI/rj7yM6TFcmCkkNWTwn3m8UUGv+XkQGLlyyFwWuJ598MsumNUc1RtEtCovxPDXtsssu2WcKTSbmlAIIp5xySgzjqbZtwIAB2bYJEyYst39sHz169HLbt9hii4qhQ4dW3h8zZkzF2muvXfH8889X2++8886raNOmTcWcOXNW2K5ow2c/+9k6f/6LX/wia8s111xTZxt69epVcdBBB63y+w8vv/xytr1jx44Vc+fOrfVnt912W+W2eN3Ydtppp1VuW7ZsWfb6bdu2rZg3b1627aGHHsr2i68re8662lbbcTj00EOz13nppZcqt73xxhsV6667bsUXvvCFym3x/PHY/fbbL2tfwZlnnpkdl/fee2+FnxcA1FQ4t9R2a9eu3XLnug033LDa+WbkyJHZ9jhvf/LJJ5XbBw8enJ3bPv7442rn+tj35z//eeW2BQsWVGy00UYVO++8c+W22s63K8o38bO4FUycODHbd9y4ccvtW/X8+eGHH1b72ZIlSyp23HHHin322afa9poZpS7xmnH+r8sZZ5yR7TNr1qxa88O7776b3b/yyitX+DqRsaq+35rHco899qj49NNPa/1ZvOaqHo/ILLVlmtqes6621Tym8Vl37do1+7w/+uijyv1+/etfZ/uNGjVquZx2ySWXVHvOaOOuu+5asTKbbrppxRFHHLHc9sL7OuGEE6ptP+ywwyo22GCDatsKv38ru63o96Suz6bws5q/d2H27NnVfu9/9rOfZfcfeeSR5fY96qijKrp3777c9q997WsVa621Vp3tglXl8j1gpWJUUPQ8NVQMpY6ROzEaJ3pNC7f99tsv66l65JFHGtW+Qg9S9FrWJXrHoqcwevAaKnobo6ewvqquyBK9w3E/ek6jN665xOcZw+ZjOHpcClEQw/X/3//7f9losuhRrSpGflUdRh/HKp7nX//6V7O1E4CWLS7JixE5VW81J00OMYomLnEqKIzmjZFLVRfwiO1xDq15KVZcHh+XnRfEyJC49ChGeMSla02Rb37+859nl4LFqJ2aqp4/45KxgnfffTcbDRPn1Bgt0xxWln+iPTHiJUacRXsaKuYWqu/8UY05Ho0Ro99jFFiM6qo611SMUIt5t37zm98s95gYzVVVHKuYemFlYoRdZNq61Pa88Ziq+StGZdX8+6jtFiPeG+Kjjz6qdX7VwmcTP6/6ta59Cz+vKt57bG+qSxLB5XvASsWEio2Z9DMKQU888USdBZ0IEY1RWOVu3XXXrXOfuGQu5lCIocwx5DourTvuuOOyVUvqq7ZVhOoS1/tXLQqFeO1Q2/wLTSXmgoqQsN122y33s1hVJS43jHkEPvvZz1a7BKCqQtBqTIAFoLzFZUL1mei85jmoUKCKS4tq217z3BRzVNacn6jq+TaWt29svol5o+K8urJVbuMyvUsvvTS7NCpWwStorlXKVpZ/otAQl4zFhOhx2VtMKRCXIUaRaEWfS2PyT2OOR2MUOtJqyz9RlIpOuZoFl5q5NPJPfbPPilZqXlGuKlxOF3MzNacoSFb9HSyIaSgKP6/6ta59qxZaa753q+/RVBSlgJWq7YS0IjUniYxCSMwjUFdvTyGsNFTM51AIQnWJySgjVMbEjDGSKK7v/9///d80YcKEbJ6p5vgcVqauk3l9JtlsSnX1fq4ocAFAc56D8jg3NeV5/U9/+lM2n1TkjVgVOEYox5xDMafSnXfemZpD5J/4nFZUNIpFZmKlvpic+/7778/m4Iy5jh588MG08847l23+aczKgTEv2oqKV/X53Y3502LkX30++6ojCesrfv9qjioMMZdrYURbYb+q22vuW9ivqnjvsTBAU/9eUL5cvgc0WPT8xKp6VcUJtuaJLSaNjN68uFyvtlvNHqVVEQEmwl6cHGM1vRWJicljmP5dd92VjRaKUVJVV61ryh6fKMTVHAL+/PPPZ19jYtGqPWc1P8PaLpurb9ui1y8+i+eee265nz377LPZCK6avc8AsLqKSdBrFqpqnm8bK3JMnFdXNDF1XOIXo2+i8HPCCSdkE6hHxmkuMVF5LF7Sv3//FY4UL7Q/RktFp1wUsiKrxeVjzZF/6nM8miP/bLHFFtnX2vJPbCv8vCnEyKuXX365Uc8RKzdGQWhlt5hUvyFikvL43GtO2fC3v/2t8uchrh6IEYBx+WNV8TsSI/4K+1UV7z1G30NTUZQCGixCTs35oG666ablerpiFZJYfjaCWk0RSGLVj4aI1zn99NPTM888k32ta+WZENfy15yHIUZWVR2uvPbaa1e2qSlcd911ld9HQIv70Wu67777ZtsiIEVvWs3PMHpYa6pv2+L5YjntGBFW9TLBWIEmindRuFvR5wQAq5NYnfYXv/hF5f34R/jtt9+e/WO6qS4VizklYy7Mquf1gkIBJs6/UUCpmoHiPBwjlJpajLKJlYfjtQqr0tUmLucvXK5VNbtFEatm/mmq7FOf41FYrbBq/olV/X70ox8t93z1bVtcKhor3MUI+KrvLeYxi5zYlKsfRiEwinu1XfJWX809p9SRRx6Z/X5ELi+I9sbIvZifrdBBGaOwonh6xx13VJubLFZ7jg7lmPOtppgjbffdd29Qu6A2Lt8DGiwue4vJHCOsxeV5s2bNygpPMRloVeecc0761a9+lc1jcPzxx2fX0Uf4iGVo77nnniy01XxMTTFZaJwwCyEreuLuvffe7JK8WJ55zJgxK3x8z54901577ZW9doyYih6heO2qk5EXru+PAlcsExwBM567IaK3NJY/juWQ4+QfoSgm2Tz//PMr5zCIIBAn+1jmOoJshLSYj6K2ObZWpW0xn0UEmShAxYSf0QN24403ZmEkltoGgOYW570YoVtT/GO25pyLjRFTAJx44onpsccey+ZNmjhxYtYRE//4bioxB1MUVkaMGJGmT5+eTVwdOSYWLonzbMxZGUWPcePGZXNWxsIicS6Pyd6jAyzm1WyoGO0S+SeKX1HgiawVC8hEwaDweit6bHSERedg5KDIA1Ewis+naoaIjHHDDTdk+SHaG8WdffbZp9mOR3SexSj52C8yYmSa2C/yUYwAq6q+bYtOv5g/K0bEDxgwICvaxetec8012QitM888MzWVON6RO2OkWryXhmjonFJRyCsU82Ie0fg9jM8mxKWjcQuRPSNjjhw5MvtdjM8uin6RuW+99dZqz3nZZZdlf5fxucXiN6+99lpWNIv3VvP3a8aMGVlRND4DaDKrvF4f0GLFssM1/7cQS83GsrK1Wbp0acW5555b0aVLl4oOHTpUDBw4sOLFF1+sdanj999/P1vmeZtttsmWdI7H7L777hVXXXVVtozvitRcNnedddap2HbbbSuOPfbYigceeKDWx9Rsw6WXXlrRt2/fis6dO2fL2G6//fYVl112WbXXjuWOTzvttGx56latWlV+FoUllmtbUrnm8sshXnfttdeueOmllyq++MUvZp9Nt27dsqWC4zOrat68edmywrHPeuutV/H1r3+94qmnnlruOetqW4jv47mrmjlzZnY84rOK5957770rHn300VqXXn7ssceqba9t6WwAqI/CuaWuW+HcVte5tXAOiqXqV3bOinP9QQcdVHH//fdX7LTTThXt2rXLzu81H1vbeW1F+SZ+FreqPvzww4oLLrigYsstt6xYc801K7p3715x5JFHZuf6gltvvTXLJ4V2RJvj/FwzW9WWk2pT9XNr3bp1lmF23nnnijPOOKNi9uzZK80k8+fPz7JdtCVySadOnSr69etX8dOf/rTa4956663sc1x33XWzxxfee105oerP4jVX9XiEGTNmZG2JTLj55ptXjBs3rtbnrKttdWWVSZMmZZ9RvPb6669fccwxx1S89tpr1fYp5LSaajtWdYn3d+KJJ9b6+Mh2K/usGqrwGrXdambBjz76qOLss8/Oflfj8/j85z9fMWXKlFqf909/+lOWy9u3b59lzfi9Wbhw4XL7Re6P47Vs2bJGvxcoaBX/aboSFwAAQPOLETAxJ06MMoY8xeVtp5xySjayq3PnzqkcxIj7+Js777zzGjzXFdTGnFIAAABQT8ccc0x2CWJcolku4hLMuEwypu6ApmSkFAAAsNoxUgpg9WekFAAAAADlWZSKYY/R0xGrVcVKAbGqRV1++MMfZqtUVb3F4wAAyon8RLmLlcSMkgJYvRW9KDVp0qRsedXRo0enmTNnpl69emXLnde2JHpBx44d05tvvll5+9e//pVrmwEAikl+AgBagqIXpcaNG5eGDx+ehg0blnr27JkmTJiQOnTokCZOnFjnY6J3r3v37pW3bt265dpmAIBikp8AgJZgjWK++JIlS9KMGTPSyJEjK7e1bt067bfffmnatGl1Pu6DDz5IW2yxRVq2bFnaZZdd0ne/+9302c9+ts6lK+NWEI9555130gYbbJCFMwCAFYk1Yd5///208cYbZzml2PLIT0GGAgCaOz8VtSg1f/78tHTp0uV66uL+s88+W+tjtttuu6wXcKeddkoLFixIV111Vdp9993T7Nmz06abbrrc/mPHjk0XX3xxs70HAKA8vPrqq7VmjZaYn4IMBQA0d34qalGqIfr375/dCiJQ7bDDDunGG29MY8aMWW7/6EWMORcKIohtvvnm2QcTcysAAKzIwoUL02abbZbWXXfdtLpa1fwUZCgAoLnzU1GLUl26dElt2rRJb7/9drXtcT/mOqiPNddcM+28887pxRdfrPXn7dq1y241RZgSqACA+iqVS9byyE9BhgIAmjs/FXVihLZt26Zdd901TZ06tdp8BXG/am/eisTw9SeffDJttNFGzdhSAIDSID8BAC1F0S/fi2HhQ4cOTX369El9+/ZN48ePT4sWLcpWkwlDhgxJm2yySTavQbjkkkvSbrvtlrbZZpv03nvvpSuvvDJb0vikk04q8jsBAMiH/AQAtARFL0odffTRad68eWnUqFHprbfeSr17905TpkypnLxzzpw51WZqf/fdd7MlkGPf9dZbL+spfPTRR7PlkAEAyoH8BAC0BK0qYp2+Mptsq1OnTtlkneZDAABWRnb4D58DANDUuaGoc0oBAAAAUJ4UpQAAAADInaIUAAAAALlTlAIAAAAgd4pSAAAAAOROUQoAAACA3ClKAQAAAJA7RSkAAAAAcqcoBQAAAEDuFKUAAAAAyJ2iFAAAAAC5U5QCAAAAIHeKUgAAAADkTlEKAAAAgNwpSgEAAACQO0UpAAAAAHKnKAUAAABA7hSlAAAAAMidohQAAAAAuVOUAgAAACB3ilIAAAAA5E5RCgAAAIDcKUoBAAAAkDtFKQAAAABypygFAAAAQO4UpQAAAADInaIUAAAAALlTlAIAAAAgd4pSAAAAAOROUQoAAACA3ClKAQAAAJA7RSkAAAAAcqcoBQAAAEDuFKUAAAAAyJ2iFAAAAAC5U5QCAAAAIHeKUgAAAADkTlEKAAAAgNwpSgEAAACQO0UpAAAAAHKnKAUAAABA7hSlAAAAAMidohQAAAAAuVOUAgAAACB3ilIAAAAA5E5RCgAAAIDcKUoBAAAAkDtFKQAAAABypygFAAAAQO4UpQAAAAAoz6LU9ddfn3r06JHat2+f+vXrl6ZPn16vx919992pVatW6dBDD232NgIAlBL5CQBY3RW9KDVp0qQ0YsSINHr06DRz5szUq1evNHDgwDR37twVPu6VV15JZ599dtpzzz1zaysAQCmQnwCAlqDoRalx48al4cOHp2HDhqWePXumCRMmpA4dOqSJEyfW+ZilS5emY445Jl188cVpq622yrW9AADFJj8BAC1BUYtSS5YsSTNmzEj77bfffxvUunV2f9q0aXU+7pJLLkldu3ZNJ5544kpfY/HixWnhwoXVbgAAq6s88lOQoQCA5rZGKqL58+dnvXbdunWrtj3uP/vss7U+5s9//nO69dZb0+OPP16v1xg7dmzWIwhAdYMGFbsFq7fJk4vdAspVHvkpyFAAtRt0lxDVUJMHC1CU2OV7q+L9999Pxx13XLr55ptTly5d6vWYkSNHpgULFlTeXn311WZvJwDA6pyfggwFALTokVIRjNq0aZPefvvtatvjfvfu3Zfb/6WXXsom6BxUpXt/2bJl2dc11lgjPffcc2nrrbeu9ph27dplNwCAliCP/BRkKACgRY+Uatu2bdp1113T1KlTq4WkuN+/f//l9t9+++3Tk08+mQ09L9wOOeSQtPfee2ffb7bZZjm/AwCAfMlPAEBLUdSRUiGWMx46dGjq06dP6tu3bxo/fnxatGhRtppMGDJkSNpkk02yeQ3at2+fdtxxx2qP79y5c/a15nYAgJZKfgIAWoKiF6WOPvroNG/evDRq1Kj01ltvpd69e6cpU6ZUTt45Z86cbEUZAAD+Q34CAFqCVhUVFRWpjMRyxp06dcom7OzYsWOxmwNQNFbfaxyr75UP2eE/fA4A/2H1vYaz+l75WFjP3KALDQAAAIDcKUoBAAAAkDtFKQAAAABypygFAAAAQO4UpQAAAADInaIUAAAAALlTlAIAAAAgd4pSAAAAAOROUQoAAACA3ClKAQAAAJA7RSkAAAAAcqcoBQAAAEDuFKUAAAAAyJ2iFAAAAAC5U5QCAAAAIHeKUgAAAADkTlEKAAAAgNwpSgEAAACQO0UpAAAAAHKnKAUAAABA7hSlAAAAAMidohQAAAAAuVOUAgAAACB3ilIAAAAA5E5RCgAAAIDcKUoBAAAAkDtFKQAAAABypygFAAAAQO4UpQAAAADInaIUAAAAALlTlAIAAAAgd4pSAAAAAOROUQoAAACA3ClKAQAAAJA7RSkAAAAAcqcoBQAAAEDuFKUAAAAAyJ2iFAAAAAC5U5QCAAAAIHeKUgAAAADkTlEKAAAAgNwpSgEAAACQO0UpAAAAAHKnKAUAAABA7hSlAAAAAMidohQAAAAAuVOUAgAAACB3ilIAAAAA5E5RCgAAAIDcKUoBAAAAsHoUpR566KEmbcT111+fevTokdq3b5/69euXpk+fXue+9957b+rTp0/q3LlzWnvttVPv3r3Tj3/84yZtDwBAc2jKDCU/AQBlWZQ64IAD0tZbb50uvfTS9OqrrzaqAZMmTUojRoxIo0ePTjNnzky9evVKAwcOTHPnzq11//XXXz9dcMEFadq0aemJJ55Iw4YNy273339/o9oBANDcmipDyU8AQNkWpV5//fV06qmnpnvuuSdttdVWWQj66U9/mpYsWbLKzzVu3Lg0fPjwLBj17NkzTZgwIXXo0CFNnDix1v332muvdNhhh6UddtghC3VnnHFG2mmnndKf//znhrwVAIDcNFWGkp8AgLItSnXp0iWdeeaZ6fHHH09/+9vf0mc+85l08sknp4033jidfvrpadasWfV6nghgM2bMSPvtt99/G9S6dXY/evJWpqKiIk2dOjU999xz6Qtf+EKt+yxevDgtXLiw2g0AoBiaIkPlkZ+CDAUAlPxE57vssksaOXJk1uv3wQcfZD10u+66a9pzzz3T7NmzV/jY+fPnp6VLl6Zu3bpV2x7333rrrToft2DBgrTOOuuktm3bpoMOOih9//vfT/vvv3+t+44dOzZ16tSp8rbZZps18J0CADSdhmaoPPJTkKEAgJItSn3yySfZ0PMDDzwwbbHFFtmcBNddd116++2304svvphtO+qoo1JzWHfddbMexsceeyxddtll2ZwKDz/8cK37RtiLEFa4NXYOLACAxihWhlqV/BRkKACgua3RkAeddtpp6a677sqGfx933HHpiiuuSDvuuGPlz2NVl6uuuiobir6yIext2rTJQlhVcb979+51Pi6GqG+zzTbZ97F6zDPPPJP15sV8CTW1a9cuuwEAFFtTZKg88lOQoQCAkhwp9fTTT2dDvt944400fvz4amGqamBa2bLHMXw8hqnHvAYFy5Yty+7379+/3u2Jx8S8BwAApawpMpT8BACU9UipWH549913T2usUf3hn376aXr00UezSTPjZwMGDFjpc8XQ8aFDh6Y+ffqkvn37ZgFt0aJF2WoyYciQIWmTTTbJevJCfI19Y+WYCFK//e1v049//ON0ww03NOStAADkpqkylPwEAJRtUWrvvfdOb775ZuratWu17THfQPwsJt+sr6OPPjrNmzcvjRo1KpucM4aTT5kypXLyzjlz5mTDzQsicMUqNa+99lpaa6210vbbb5/uuOOO7HkAAEpZU2Uo+QkAaAlaVcSkBqsoQk7MW7DhhhtW2/78889nvXClvGRwtC1WkInw17Fjx2I3B6BoBg0qdgtWb5MnF7sFrI7ZQYYCWP0NukuIaqjJgwWocrGwnrlhlUZKHX744dnXVq1apeOPP77a5JfRs/fEE09kQ9IBAPgvGQoAoJFFqahyhRhcFcsKx/DvqpNu7rbbbmn48OGr8pQAAC2eDAUA0Mii1G233ZZ97dGjRzr77LOzZYsBAFgxGQoAoAlX3wMAYNXIUAAADShK7bLLLmnq1KlpvfXWSzvvvHM2J0JdZs6cWd+nBQBo0WQoAIBGFqW+/OUvV07Keeihh9b3YQAAZU2GAgCoXauKmHGzjFjOGOA/BlnNuFEmW9G4bMgO/+FzAPiPQXcJUQ01ebAAVS4W1jM3tM61VQAAAACwKpfvxTwIK5oDoap33nmnMW0CAGgxZCgAgEYWpcaPH1/fXQEA+D8yFABAI4tSQ4cOre+uAAD8HxkKAKCRRamYpKowOVV8vyImvwQA+A8ZCgCgCeaUevPNN1PXrl1T586da50bIRbyi+1Lly6t79MCALRoMhQAQCOLUg8++GBaf/31s+8feuih+j4MAKCsyVAAAI0sSg0YMKDW7wEAqJsMBQDQyKJUTe+++2669dZb0zPPPJPd79mzZxo2bFhlTyAAAMuToQAA/qN1aoBHHnkk9ejRI1177bVZsIpbfL/llltmPwMAYHkyFABAI0dKnXLKKenoo49ON9xwQ2rTpk22LSbmPPnkk7OfPfnkkw15WgCAFk2GAgBo5EipF198MZ111lmVYSrE9yNGjMh+BgDA8mQoAIBGFqV22WWXynkQqoptvXr1ashTAgC0eDIUAEADLt974oknKr8//fTT0xlnnJH16O22227Ztr/+9a/p+uuvT5dffnl9nxIAoMWToQAAateqoqKiItVD69atU6tWrdLKdo99Ym6EUrVw4cLUqVOntGDBgtSxY8diNwegaAYNKnYLVm+TJxe7Bawu2UGGAmhZBt0lRDXU5MECVLlYWM/cUO+RUi+//HJTtQ0AoGzIUAAAjSxKbbHFFvXdFQCA/yNDAQA0sihVm6effjrNmTMnLVmypNr2Qw45pDFPCwDQoslQAAANLEr985//TIcddlh68sknq82REN+HUp4PAQCgWGQoAID/ap0aIFaN2XLLLdPcuXNThw4d0uzZs9MjjzyS+vTpkx5++OGGPCUAQIsnQwEANHKk1LRp09KDDz6YunTpkq0oE7c99tgjjR07Nlvq+B//+EdDnhYAoEWToQAAGjlSKoaWr7vuutn3EareeOONyok8n3vuuYY8JQBAiydDAQA0cqTUjjvumGbNmpUNP+/Xr1+64oorUtu2bdNNN92Uttpqq4Y8JQBAiydDAQA0sij1ne98Jy1atCj7/pJLLkkHH3xw2nPPPdMGG2yQJk2a1JCnBABo8WQoAIBGFqUGDhxY+f0222yTnn322fTOO++k9dZbr3L1GAAAqpOhAAAaWZSq6tVXX82+brbZZo19KgCAsiFDAQDlrkETnX/66afpwgsvTJ06dUo9evTIbvF9DEn/5JNPmr6VAAAtgAwFANDIkVKnnXZauvfee7PJOfv371+5xPFFF12U/v3vf6cbbrihIU8LANCiyVAAAI0sSt15553p7rvvTl/60pcqt+20007Z8PPBgwcLVAAAtZChAAAaefleu3btsuHmNcXyxrGsMQAAy5OhAAAaWZQ69dRT05gxY9LixYsrt8X3l112WfYzAACWJ0MBADTg8r3DDz+82v0//OEPadNNN029evXK7s+aNSstWbIk7bvvvvV9SgCAFk+GAgBoZFEqVoap6ogjjqh233LGAADLk6EAABpZlLrtttvquysAAP9HhgIAaMLV9wrmzZuXnnvuuez77bbbLm244YaNeToAgLIgQwEANHCi80WLFqUTTjghbbTRRukLX/hCdtt4443TiSeemD788MOmbyUAQAsgQwEANLIoNWLEiPTHP/4xTZ48Ob333nvZ7b777su2nXXWWQ15SgCAFk+GAgBo5OV7P//5z9M999yT9tprr8ptBx54YFprrbXSV77ylXTDDTc05GkBAFo0GQoAoJEjpWJ4ebdu3Zbb3rVrV0PPAQDqIEMBADSyKNW/f/80evTo9PHHH1du++ijj9LFF1+c/QwAgOXJUAAAjbx8b/z48emAAw5Im266aerVq1e2bdasWal9+/bp/vvvb8hTAgC0eDIUAEAji1Kf+9zn0gsvvJB+8pOfpGeffTbbNnjw4HTMMcdkcyIAALA8GYrVxaC7BhW7Cau1yYMnF7sJAC3z8r1PPvkkbb311ulf//pXGj58eLr66quz20knndTgMHX99denHj16ZL2E/fr1S9OnT69z35tvvjntueeeab311stu++233wr3BwAoBU2doeQnAKDsilJrrrlmtXkQGmvSpEnZ8sgxv8LMmTOzoewDBw5Mc+fOrXX/hx9+OOtRfOihh9K0adPSZpttlr74xS+m119/vcnaBADQ1JoyQ8lPAEDZTnR+yimnpO9973vp008/bXQDxo0bl/UWDhs2LPXs2TNNmDAhdejQIU2cOLHW/WO4+8knn5x69+6dtt9++3TLLbekZcuWpalTpza6LQAAzampMpT8BACU7ZxSjz32WBZiHnjggWxuhLXXXrvaz++99956Pc+SJUvSjBkz0siRIyu3tW7dOhtSHr149RHLJ8dw+PXXX38V3wUAQL6aIkPJTwBAWRelOnfunI444ohGv/j8+fPT0qVLU7du3aptj/uFyT9X5txzz00bb7xxFsRqs3jx4uxWsHDhwka2GgCgYZoiQ+WRn4IMBQCUVFEqhnlfeeWV6fnnn8966fbZZ5900UUXFW21mMsvvzzdfffd2TwJMclnbcaOHZsuvvji3NtG6Rhk8ZhGmWzxGMiF/1e17P9XlVKGqk9+CjIUAFBSc0pddtll6fzzz0/rrLNO2mSTTdK1116bzY3QUF26dElt2rRJb7/9drXtcb979+4rfOxVV12VhaoY/r7TTjvVuV8MbV+wYEHl7dVXX21wewEAGqIpM1Qe+SnIUABASRWlbr/99vSDH/wg3X///emXv/xlmjx5cjZxZvT+NUTbtm3TrrvuWm2SzcKkm/3796/zcVdccUUaM2ZMmjJlSurTp88KX6Ndu3apY8eO1W4AAHlqygyVR34KMhQAUFJFqTlz5qQDDzyw8n7MQ9CqVav0xhtvNLgBsZzxzTffnH70ox+lZ555Jn3zm99MixYtylaTCUOGDKk2kWesWHPhhRdmq8v06NEjvfXWW9ntgw8+aHAbAACaU1NnKPkJACi7OaVi+eKacw+sueaa2eotDXX00UenefPmpVGjRmXhKJYqjh68wuSdEeJiRZmCG264IZuL4cgjj6z2PKNHj87mZgAAKDVNnaHkJwCg7IpSFRUV6fjjj8+Gcxd8/PHH6Rvf+Ea1JY3rs5xxVaeeemp2q01MwlnVK6+8skrPDQBQbM2RoeQnAKCsilJDhw5dbtuxxx7blO0BAGhxZCgAaHqD7rJ8cWNMHjx59SpK3Xbbbc3XEgCAFkqGAgBo5ETnAAAAANAUFKUAAAAAyJ2iFAAAAAC5U5QCAAAAIHeKUgAAAADkTlEKAAAAgNwpSgEAAACQO0UpAAAAAHKnKAUAAABA7hSlAAAAAMidohQAAAAAuVOUAgAAACB3ilIAAAAA5E5RCgAAAIDcKUoBAAAAkDtFKQAAAABypygFAAAAQO4UpQAAAADInaIUAAAAALlTlAIAAAAgd4pSAAAAAOROUQoAAACA3ClKAQAAAJA7RSkAAAAAcqcoBQAAAEDuFKUAAAAAyJ2iFAAAAAC5U5QCAAAAIHeKUgAAAADkTlEKAAAAgNwpSgEAAACQO0UpAAAAAHKnKAUAAABA7hSlAAAAAMidohQAAAAAuVOUAgAAACB3ilIAAAAA5E5RCgAAAIDcKUoBAAAAkDtFKQAAAABypygFAAAAQO7WyP8lAQAAqK9Bdw0qdhNWa5MHTy52E4A6GCkFAAAAQO4UpQAAAADInaIUAAAAALlTlAIAAAAgd4pSAAAAAOROUQoAAACA8itKXX/99alHjx6pffv2qV+/fmn69Ol17jt79ux0xBFHZPu3atUqjR8/Pte2AgCUChkKAFjdFbUoNWnSpDRixIg0evToNHPmzNSrV680cODANHfu3Fr3//DDD9NWW22VLr/88tS9e/fc2wsAUApkKACgJShqUWrcuHFp+PDhadiwYalnz55pwoQJqUOHDmnixIm17v/5z38+XXnllemrX/1qateuXe7tBQAoBTIUANASFK0otWTJkjRjxoy03377/bcxrVtn96dNm1asZgEAlDQZCgBoKdYo1gvPnz8/LV26NHXr1q3a9rj/7LPPNtnrLF68OLsVLFy4sMmeGwAgbzIUANBSFK0olZexY8emiy++ONfXHDQo15drcSZPLnYLaE7+PhrO3wbQ0jMUAFBeinb5XpcuXVKbNm3S22+/XW173G/KCThHjhyZFixYUHl79dVXm+y5AQDyJkMBAC1F0YpSbdu2TbvuumuaOnVq5bZly5Zl9/v3799krxOTeXbs2LHaDQBgdSVDAQAtRVEv34uljIcOHZr69OmT+vbtm8aPH58WLVqUrSQThgwZkjbZZJNs+HhhYs+nn3668vvXX389Pf7442mdddZJ22yzTTHfCgBAbmQoAKAlKGpR6uijj07z5s1Lo0aNSm+99Vbq3bt3mjJlSuXEnXPmzMlWkyl444030s4771x5/6qrrspuAwYMSA8//HBR3gMAQN5kKACgJSj6ROennnpqdqtNzZDUo0ePVFFRkVPLAABKlwwFAKzuijanFAAAAADlS1EKAAAAgNwpSgEAAACQO0UpAAAAAHKnKAUAAABA7hSlAAAAAMjdGvm/JAAALG/QXYOK3YTV1uTBk4vdBABYZUZKAQAAAJA7RSkAAAAAcqcoBQAAAEDuFKUAAAAAyJ2iFAAAAAC5U5QCAAAAIHeKUgAAAADkTlEKAAAAgNwpSgEAAACQO0UpAAAAAHKnKAUAAABA7hSlAAAAAMidohQAAAAAuVOUAgAAACB3ilIAAAAA5E5RCgAAAIDcKUoBAAAAkDtFKQAAAABypygFAAAAQO4UpQAAAADInaIUAAAAALlTlAIAAAAgd4pSAAAAAOROUQoAAACA3ClKAQAAAJA7RSkAAAAAcqcoBQAAAEDuFKUAAAAAyJ2iFAAAAAC5U5QCAAAAIHeKUgAAAADkTlEKAAAAgNwpSgEAAACQO0UpAAAAAHKnKAUAAABA7hSlAAAAAMidohQAAAAAuVOUAgAAACB3ilIAAAAA5E5RCgAAAIDcKUoBAAAAkDtFKQAAAABypygFAAAAQHkWpa6//vrUo0eP1L59+9SvX780ffr0Fe7/s5/9LG2//fbZ/p/73OfSb3/729zaCgBQCuQnAGB1V/Si1KRJk9KIESPS6NGj08yZM1OvXr3SwIED09y5c2vd/9FHH02DBw9OJ554YvrHP/6RDj300Oz21FNP5d52AIBikJ8AgJag6EWpcePGpeHDh6dhw4alnj17pgkTJqQOHTqkiRMn1rr/Nddckw444IB0zjnnpB122CGNGTMm7bLLLum6667Lve0AAMUgPwEALUFRi1JLlixJM2bMSPvtt99/G9S6dXZ/2rRptT4mtlfdP0TPYF37AwC0JPITANBSrFHMF58/f35aunRp6tatW7Xtcf/ZZ5+t9TFvvfVWrfvH9tosXrw4uxUsWLAg+7pw4cLUXD75pNmeuiw09aFxPBrH8SgdjkVpcTxKSzOe1iszQ0VFRSoFeeSnomWoD/0hNFRTHxfHonEcj9LieJQOx6K0LGzGc3p981NRi1J5GDt2bLr44ouX277ZZpsVpT2sXKdOxW4BVTkepcOxKC2OR/kdj/fffz91KqMDL0OtXjqdVD6/m6sDx6O0OB6lw7Eov+Px/kryU1GLUl26dElt2rRJb7/9drXtcb979+61Pia2r8r+I0eOzCYCLVi2bFl655130gYbbJBatWqVyk1UKyNMvvrqq6ljx47Fbk7ZczxKh2NRWhyP0lLuxyN6+CJQbbzxxqkU5JGfggxVXbn/HZQSx6K0OB6lxfEoHeV+LCrqmZ+KWpRq27Zt2nXXXdPUqVOzFWAKgSfun3rqqbU+pn///tnPv/Wtb1Vu+/3vf59tr027du2yW1WdO3dO5S7+KMrxD6NUOR6lw7EoLY5HaSnn41FKI6TyyE9BhqpdOf8dlBrHorQ4HqXF8Sgd5XwsOtUjPxX98r3ogRs6dGjq06dP6tu3bxo/fnxatGhRtppMGDJkSNpkk02yIeThjDPOSAMGDEhXX311Ouigg9Ldd9+d/v73v6ebbrqpyO8EACAf8hMA0BIUvSh19NFHp3nz5qVRo0Zlk2327t07TZkypXIyzjlz5mQryhTsvvvu6c4770zf+c530vnnn5+23Xbb9Mtf/jLtuOOORXwXAAD5kZ8AgJag6EWpEEPN6xpu/vDDDy+37aijjspurLoYhj969OjlhuNTHI5H6XAsSovjUVocj9IkP+XL30HpcCxKi+NRWhyP0uFY1E+rilJZ3xgAAACAsvHfcd0AAAAAkBNFKQAAAABypygFRbDXXntVW5YbqJu/FwCC8wGsGn8zrA4UpQAAAADInaIUAAAAALlTlCojU6ZMSXvssUfq3Llz2mCDDdLBBx+cXnrppWI3q2x9+umn2VLenTp1Sl26dEkXXnhhshhmcSxbtixdccUVaZtttsmWbN18883TZZddVuxmlaVFixalIUOGpHXWWSdttNFG6eqrry52k1K5/22MHTs2bbnllmmttdZKvXr1Svfcc0+xmwW5kp9Ki/xUWmSo0iFDlRYZqv4Upcrsf1QjRoxIf//739PUqVNT69at02GHHZb9wZC/H/3oR2mNNdZI06dPT9dcc00aN25cuuWWW4rdrLI0cuTIdPnll2fB9umnn0533nln6tatW7GbVZbOOeec9Mc//jHdd9996YEHHkgPP/xwmjlzZrGbVbYiTN1+++1pwoQJafbs2enMM89Mxx57bHaMoFzIT6VFfiotMlTpkKFKiwxVf60qdC2Urfnz56cNN9wwPfnkk2nHHXcsdnPKbtLBuXPnZv+DatWqVbbtvPPOS7/61a+yEzr5ef/997O/g+uuuy6ddNJJxW5OWfvggw+yUQh33HFHOuqoo7Jt77zzTtp0003T1772tTR+/PhiN7GsLF68OK2//vrpD3/4Q+rfv3/l9vg7+fDDD7N/eEA5kp+KR34qLTJU6ZChSosMtWqMlCojL7zwQho8eHDaaqutUseOHVOPHj2y7XPmzCl208rSbrvtVhmoQvwPK47R0qVLi9qucvPMM89kJ45999232E0pe3E5zJIlS1K/fv0qt8UJfbvttitqu8rViy++mAWn/fffP7sUoHCLXj+XLlFO5KfSIj+VDhmqdMhQpUWGWjVrrOL+rMYGDRqUtthii3TzzTenjTfeOBt2Hj188T8wKFdxjTdQe69r+M1vfpM22WSTaj+LeUOgXMhPUDsZCmonQ60aI6XKxL///e/03HPPpe985ztZb8YOO+yQ3n333WI3q6z97W9/q3b/r3/9a9p2221TmzZtitamchSfeYSqmCeE4tp6663TmmuuWe1vI/4/9fzzzxe1XeWqZ8+eWXCK0SAxgW3V22abbVbs5kEu5KfSIz+VDhmqdMhQpUWGWjVGSpWJ9dZbL7vO+KabbspWY4g/kLgGn+KJYxATp37961/PJiH8/ve/b5WMImjfvn0699xz07e//e3Utm3b9D//8z9p3rx52XwVJ554YrGbV1ZiWHN85jFRZ/z/qmvXrumCCy7IJhUmf+uuu246++yzs4k5Y2RIrD62YMGC9Je//CW7hGno0KHFbiI0O/mp9MhPpUOGKh0yVGmRoVaNolSZiP8h3X333en000/PhpzH9cXXXnttNmEkxRFLtn700Uepb9++We/eGWeckU1ESP5ixZhYyWfUqFHpjTfeyP7h8Y1vfKPYzSpLV155ZTbkOS6XiRP6WWedlZ3EKY4xY8Zkk9jGCjL//Oc/U+fOndMuu+ySzj///GI3DXIhP5Ue+am0yFClQ4YqLTJU/Vl9DwAAAIDcGc8HAAAAQO4UpQAAAADInaIUAAAAALlTlAIAAAAgd4pSAAAAAOROUQoAAACA3ClKAQAAAJA7RSkAAAAAcqcoBZStV155JbVq1So9/vjjxW4KAMBqQX4CmpKiFLBai1C0ottFF11U7CYCAJQU+QkoFWsUuwEAjfHmm29Wfj9p0qQ0atSo9Nxzz1VuW2eddYrUMgCA0iQ/AaXCSClgtda9e/fKW6dOnbLevcL9rl27pnHjxqVNN900tWvXLvXu3TtNmTKlzudaunRpOuGEE9L222+f5syZk22777770i677JLat2+fttpqq3TxxRenTz/9tPIx8Xq33HJLOuyww1KHDh3Stttum371q1/l8t4BABpCfgJKhaIU0GJdc8016eqrr05XXXVVeuKJJ9LAgQPTIYcckl544YXl9l28eHE66qijsvkR/vSnP6XNN988+zpkyJB0xhlnpKeffjrdeOON6Yc//GG67LLLqj02gtZXvvKV7DUOPPDAdMwxx6R33nknx3cKANA05CcgT4pSQIsVYercc89NX/3qV9N2222Xvve972W9fePHj6+23wcffJAOOuigNG/evPTQQw+lDTfcsDIsnXfeeWno0KFZL9/++++fxowZk4Wrqo4//vg0ePDgtM0226Tvfve72fNNnz491/cKANAU5CcgT+aUAlqkhQsXpjfeeCP9z//8T7XtcX/WrFnVtkUgiiHqDz74YFprrbUqt8d+f/nLX6r17MUQ9Y8//jh9+OGH2XDzsNNOO1X+fO21104dO3ZMc+fObcZ3BwDQ9OQnIG+KUkDZiyHjd9xxR5o2bVraZ599KrdHj1309h1++OHLPSbmSChYc801q/0s5klYtmxZM7caAKB45CegKShKAS1S9LZtvPHGWU/dgAEDKrfH/b59+1bb95vf/Gbacccds/kSfvOb31TuHxN0xko0MawcAKClk5+AvClKAS3WOeeck0aPHp223nrrbC6E2267LZuI8yc/+cly+5522mnZ0PKDDz44/e53v0t77LFHtjxy3I9JO4888sjUunXrbEj6U089lS699NKivCcAgOYkPwF5UpQCWqzTTz89LViwIJ111lnZHAU9e/bMlhuOZYdr861vfSsbNh7D0WPp41ht5te//nW65JJLskk+Y5h5LHd80kkn5f5eAADyID8BeWpVUVFRkesrAgAAAFD2Whe7AQAAAACUH0UpAAAAAHKnKAUAAABA7hSlAAAAAMidohQAAAAAuVOUAgAAACB3ilIAAAAA5E5RCgAAAIDcKUoBAAAAkDtFKQAAAABypygFAAAAQO4UpQAAAABIefv/YVPJi+AnACwAAAAASUVORK5CYII=",
      "text/plain": [
       "<Figure size 1200x400 with 2 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Token counts from 100 samples:\n",
      "  a: 8 (0.08)\n",
      "  b: 21 (0.21)\n",
      "  c: 12 (0.12)\n",
      "  d: 40 (0.40)\n",
      "  e: 19 (0.19)\n"
     ]
    }
   ],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "\n",
    "# Sample 100 tokens and plot the empirical distribution\n",
    "n_samples = 100\n",
    "samples = [sample_binary_search(cdf) for _ in range(n_samples)]\n",
    "\n",
    "# Count occurrences\n",
    "counts = np.zeros(len(vocab))\n",
    "for s in samples:\n",
    "    counts[s] += 1\n",
    "\n",
    "# Create bar plot comparing true vs empirical distribution\n",
    "fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(12, 4))\n",
    "\n",
    "# True distribution\n",
    "ax1.bar(vocab, probs, color='blue', alpha=0.7)\n",
    "ax1.set_title('True Distribution')\n",
    "ax1.set_xlabel('Token')\n",
    "ax1.set_ylabel('Probability')\n",
    "ax1.set_ylim([0, 0.5])\n",
    "\n",
    "# Empirical distribution\n",
    "empirical_probs = counts / n_samples\n",
    "ax2.bar(vocab, empirical_probs, color='green', alpha=0.7)\n",
    "ax2.set_title(f'Empirical Distribution (n={n_samples})')\n",
    "ax2.set_xlabel('Token')\n",
    "ax2.set_ylabel('Probability')\n",
    "ax2.set_ylim([0, 0.5])\n",
    "\n",
    "plt.tight_layout()\n",
    "plt.show()\n",
    "\n",
    "print(\"Token counts from 100 samples:\")\n",
    "for token, count in zip(vocab, counts):\n",
    "    print(f\"  {token}: {int(count)} ({count/n_samples:.2f})\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "pytorch",
   "metadata": {},
   "source": [
    "## PyTorch\n",
    "\n",
    "Categorical sampling is implemented in common libraries like PyTorch."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "a5nlvbb56mu",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "tensor([3, 3, 1, 3, 1, 2, 3, 2, 0, 1, 3, 3, 3, 0, 3, 3, 2, 3, 3, 1, 0, 2, 3, 1,\n",
       "        3, 3, 3, 2, 0, 4, 3, 0, 2, 3, 4, 3, 2, 2, 4, 3, 2, 1, 1, 0, 4, 3, 3, 3,\n",
       "        3, 2, 3, 4, 3, 3, 4, 4, 2, 4, 4, 4, 3, 3, 4, 1, 4, 1, 3, 3, 4, 1, 3, 3,\n",
       "        1, 1, 3, 3, 4, 3, 3, 1, 1, 3, 3, 3, 3, 1, 0, 4, 1, 0, 3, 3, 4, 2, 3, 3,\n",
       "        1, 3, 3, 4])"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "import torch\n",
    "\n",
    "# Sample 100 times using PyTorch\n",
    "torch_probs = torch.tensor(probs)\n",
    "categorical = torch.distributions.Categorical(probs=torch_probs)\n",
    "categorical.sample((100,))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6c252fa3",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "anlp",
   "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.10.18"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}