backprop_ninja.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import torch \n",
    "from torch.nn import functional as F\n",
    "import matplotlib.pyplot as plt \n",
    "from icecream import ic \n",
    "from sys import exit as e\n",
    "\n",
    "plt.rcParams['figure.figsize'] = (8, 6)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "32033\n"
     ]
    }
   ],
   "source": [
    "words = open(\"./data/names.txt\", 'r').read().splitlines()\n",
    "print(len(words))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "{1: 'a',\n",
       " 2: 'b',\n",
       " 3: 'c',\n",
       " 4: 'd',\n",
       " 5: 'e',\n",
       " 6: 'f',\n",
       " 7: 'g',\n",
       " 8: 'h',\n",
       " 9: 'i',\n",
       " 10: 'j',\n",
       " 11: 'k',\n",
       " 12: 'l',\n",
       " 13: 'm',\n",
       " 14: 'n',\n",
       " 15: 'o',\n",
       " 16: 'p',\n",
       " 17: 'q',\n",
       " 18: 'r',\n",
       " 19: 's',\n",
       " 20: 't',\n",
       " 21: 'u',\n",
       " 22: 'v',\n",
       " 23: 'w',\n",
       " 24: 'x',\n",
       " 25: 'y',\n",
       " 26: 'z',\n",
       " 0: '.'}"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "unique_chars = sorted(list(set(''.join(words))))\n",
    "stoi = {c: i+1 for i, c in enumerate(unique_chars)}\n",
    "stoi['.'] = 0\n",
    "\n",
    "itos = {i:c for c,i in stoi.items()}\n",
    "\n",
    "itos"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(torch.Size([182625, 3]), torch.Size([22655, 3]), torch.Size([22866, 3]))"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "block_size = 3\n",
    "\n",
    "def build_dataset(words):\n",
    "  X, y = [], []\n",
    "  for word in words:\n",
    "    context = [0] * block_size\n",
    "    for ch in word + '.':\n",
    "      X.append(context)\n",
    "      y.append(stoi[ch])\n",
    "\n",
    "      context = context[1:] + [stoi[ch]]\n",
    "\n",
    "  X = torch.tensor(X)\n",
    "  y = torch.tensor(y)\n",
    "  return X, y\n",
    "\n",
    "import random \n",
    "random.seed(42)\n",
    "random.shuffle(words)\n",
    "t1 = int(0.8 * len(words))\n",
    "t2 = int(0.9 * len(words))\n",
    "\n",
    "Xtr, ytr = build_dataset(words[:t1])\n",
    "Xva, yva = build_dataset(words[t1:t2])\n",
    "Xte, yte = build_dataset(words[t2:])\n",
    "\n",
    "Xtr.shape, Xva.shape, Xte.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "# utility function we will use later when comparing manual gradients to PyTorch gradients\n",
    "def cmp(s, dt, t):\n",
    "  ex = torch.all(dt == t.grad).item()\n",
    "  app = torch.allclose(dt, t.grad)\n",
    "  maxdiff = (dt - t.grad).abs().max().item()\n",
    "  print(f'{s:15s} | exact: {str(ex):5s} | approximate: {str(app):5s} | maxdiff: {maxdiff}')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "4137\n"
     ]
    }
   ],
   "source": [
    "n_embed = 10 \n",
    "n_hidden = 64\n",
    "vocab_size = len(unique_chars)+1\n",
    "\n",
    "g = torch.Generator().manual_seed(2147483647)\n",
    "C = torch.randn((vocab_size, n_embed), generator=g)\n",
    "\n",
    "#Layer 1\n",
    "W1 = torch.randn((n_embed * block_size, n_hidden), generator=g) * (5/3)/(n_embed * block_size)**0.5\n",
    "b1 = torch.randn((n_hidden), generator=g) * 0.1\n",
    "#Layer 2  \n",
    "W2 = torch.randn((n_hidden, vocab_size), generator=g) * 0.1\n",
    "b2 = torch.randn((vocab_size), generator=g) * 0.1\n",
    "#BatchNorm \n",
    "bngain = torch.randn((1, n_hidden), generator=g) * 0.1 + 1.\n",
    "bnbias = torch.randn((1, n_hidden), generator=g) * 0.1\n",
    "\n",
    "parameters = [C, W1, b1, W2, b2, bngain, bnbias]\n",
    "for p in parameters:\n",
    "  p.requires_grad=True\n",
    "\n",
    "print(sum(p.nelement() for p in parameters))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "batch_size = 32 \n",
    "n = 32 \n",
    "idx = torch.randint(0, Xtr.shape[0], (n, ), generator=g)\n",
    "Xb, yb = Xtr[idx], ytr[idx]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "tensor(3.5571, grad_fn=<NegBackward0>)"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "emb = C[Xb]\n",
    "embcat = emb.view(emb.shape[0], -1)\n",
    "hprebn = embcat @ W1 + b1\n",
    "#BatchNorm\n",
    "bnmeani = (1/n) * hprebn.sum(0, keepdim=True)\n",
    "bndiff = hprebn - bnmeani\n",
    "bndiff2 = bndiff**2\n",
    "bnvar = 1/(n-1)  * (bndiff2).sum(0, keepdim=True)\n",
    "bnvar_inv = (bnvar+1e-5) ** -0.5\n",
    "bnraw = bndiff * bnvar_inv\n",
    "hpreact = bngain * bnraw + bnbias\n",
    "#Non-linearity\n",
    "h = torch.tanh(hpreact)\n",
    "# Layer 2\n",
    "logits = h @ W2 + b2 \n",
    "# Loss \n",
    "logit_maxes = logits.max(1, keepdim=True).values\n",
    "norm_logits = logits - logit_maxes\n",
    "counts = norm_logits.exp()\n",
    "counts_sum = counts.sum(1, keepdim=True)\n",
    "counts_sum_inv = counts_sum ** -1\n",
    "probs = counts * counts_sum_inv\n",
    "logprobs = probs.log()\n",
    "loss = -logprobs[range(n), yb].mean()\n",
    "\n",
    "# PyTorch backward pass\n",
    "for p in parameters:\n",
    "  p.grad = None\n",
    "for t in [logprobs, probs, counts, counts_sum, counts_sum_inv, # afaik there is no cleaner way\n",
    "          norm_logits, logit_maxes, logits, h, hpreact, bnraw,\n",
    "         bnvar_inv, bnvar, bndiff2, bndiff, hprebn, bnmeani,\n",
    "         embcat, emb]:\n",
    "  t.retain_grad()\n",
    "loss.backward()\n",
    "loss"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$\\frac{\\delta p}{\\delta csi} = \\frac{\\delta p1}{csi} + \\frac{\\delta p2}{csi} + \\frac{\\delta p3}{csi}$\n",
    "\n",
    "Because\n",
    "\n",
    "$p = [counts[:0] \\times csi , counts[:, 1] \\times csi, counts[:,2] \\times csi]$\n",
    "\n",
    "Thats why, the gradients are added over the column."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "torch.Size([32, 3, 10]) torch.Size([32, 30]) torch.Size([27, 10]) torch.Size([32, 3])\n"
     ]
    }
   ],
   "source": [
    "# embcat = emb.view(emb.shape[0], -1)\n",
    "#emb = C[Xb]\n",
    "print(emb.shape, embcat.shape, C.shape, Xb.shape)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "logprobs        | exact: True  | approximate: True  | maxdiff: 0.0\n",
      "probs           | exact: True  | approximate: True  | maxdiff: 0.0\n",
      "counts_sum_inv  | exact: True  | approximate: True  | maxdiff: 0.0\n",
      "counts_sum      | exact: True  | approximate: True  | maxdiff: 0.0\n",
      "counts          | exact: True  | approximate: True  | maxdiff: 0.0\n",
      "norm_logits     | exact: True  | approximate: True  | maxdiff: 0.0\n",
      "logit_maxes     | exact: True  | approximate: True  | maxdiff: 0.0\n",
      "logits          | exact: True  | approximate: True  | maxdiff: 0.0\n",
      "h               | exact: True  | approximate: True  | maxdiff: 0.0\n",
      "W2              | exact: True  | approximate: True  | maxdiff: 0.0\n",
      "b2              | exact: True  | approximate: True  | maxdiff: 0.0\n",
      "hpreact         | exact: False | approximate: True  | maxdiff: 4.656612873077393e-10\n",
      "bngain          | exact: False | approximate: True  | maxdiff: 1.862645149230957e-09\n",
      "bnbias          | exact: False | approximate: True  | maxdiff: 3.725290298461914e-09\n",
      "bnraw           | exact: False | approximate: True  | maxdiff: 9.313225746154785e-10\n",
      "bnvar_inv       | exact: False | approximate: True  | maxdiff: 1.862645149230957e-09\n",
      "bnvar           | exact: False | approximate: True  | maxdiff: 9.313225746154785e-10\n",
      "bndiff2         | exact: False | approximate: True  | maxdiff: 2.9103830456733704e-11\n",
      "bndiff          | exact: False | approximate: True  | maxdiff: 9.313225746154785e-10\n",
      "bnmeani         | exact: False | approximate: True  | maxdiff: 3.725290298461914e-09\n",
      "hprebn          | exact: False | approximate: True  | maxdiff: 9.313225746154785e-10\n",
      "embcat          | exact: False | approximate: True  | maxdiff: 1.3969838619232178e-09\n",
      "W1              | exact: False | approximate: True  | maxdiff: 5.587935447692871e-09\n",
      "b1              | exact: False | approximate: True  | maxdiff: 4.190951585769653e-09\n",
      "emb             | exact: False | approximate: True  | maxdiff: 1.3969838619232178e-09\n",
      "C               | exact: False | approximate: True  | maxdiff: 7.450580596923828e-09\n"
     ]
    }
   ],
   "source": [
    "# Exercise 1: backprop through the whole thing manually, \n",
    "# backpropagating through exactly all of the variables \n",
    "# as they are defined in the forward pass above, one by one\n",
    "\n",
    "# -----------------\n",
    "# YOUR CODE HERE :)\n",
    "\n",
    "# grads have same shape as original tensor. Consider gradients of only indexed tensor. Others are 0\n",
    "dlogprobs = torch.zeros_like(logprobs) # dL/dlogprobs \n",
    "dlogprobs[range(n), yb] = -1.0/n\n",
    "\n",
    "dprobs = 1.0/probs * dlogprobs # dL/dprobs = dL/dlogprobs * dlogprobs/dprobs\n",
    "\n",
    "# We sum across columns because internally, PyTorch has broadcasted the single col of counts_sum_inv\n",
    "# for multiplication (counts * counts_sum_inv)\n",
    "dcounts_sum_inv = (counts * dprobs).sum(1, keepdim=True) # dL/dcounts_sum_inv = dL/dlogprobs * dlogprobs/dprobs * dprobs/dcounts_sum_inv\n",
    "dcounts = counts_sum_inv * dprobs\n",
    "\n",
    "dcounts_sum = -1 * counts_sum**-2 * dcounts_sum_inv\n",
    "\n",
    "# Note: https://www.perplexity.ai/search/plt-draw-bar-graph-from-a-list-yH4Zu62pT3qj_qj9cC3qyQ#7\n",
    "# Second dcounts branch. So add and accumulate the gradients\n",
    "dcounts += torch.ones_like(counts) * dcounts_sum\n",
    "\n",
    "# SImple\n",
    "dnorm_logits = norm_logits.exp() * dcounts\n",
    "\n",
    "dlogits = dnorm_logits.clone()\n",
    "dlogit_maxes = (torch.ones_like(logits) * -1. * dnorm_logits).sum(1, keepdim=True)\n",
    "\n",
    "dlogits += F.one_hot(logits.max(dim=1).indices, num_classes=logits.shape[1]) * dlogit_maxes\n",
    "\n",
    "dh = dlogits@ W2.T\n",
    "dW2 = h.T @ dlogits\n",
    "db2 = dlogits.sum(0)\n",
    "\n",
    "dhpreact = (1. - h**2) * dh\n",
    "\n",
    "dbngain = (bnraw * dhpreact).sum(0, keepdim=True)\n",
    "dbnraw = bngain * dhpreact\n",
    "dbnbias = dhpreact.sum(0, keepdim=True)\n",
    "\n",
    "dbndiff = bnvar_inv * dbnraw\n",
    "dbnvar_inv = (bndiff * dbnraw).sum(0, keepdim=True)\n",
    "\n",
    "dbnvar = (-0.5 * (bnvar + 1e-5) ** -1.5) * dbnvar_inv\n",
    "dbndiff2 = (1/(n-1)) * torch.ones_like(bndiff2) * dbnvar\n",
    "dbndiff += 2*bndiff * dbndiff2\n",
    "\n",
    "dhprebn = 1. * dbndiff\n",
    "dbnmeani = (-torch.ones_like(bndiff) * dbndiff).sum(0, keepdim=True) \n",
    "\n",
    "dhprebn += torch.ones_like(hprebn) * (1/n) * dbnmeani\n",
    "\n",
    "dembcat = dhprebn @ W1.T \n",
    "dW1 = embcat.T @ dhprebn\n",
    "db1 = dhprebn.sum(0)\n",
    "\n",
    "demb = dembcat.view(emb.shape)\n",
    "\n",
    "# INDEXING\n",
    "dC = torch.zeros_like(C)\n",
    "for i in range(Xb.shape[0]):\n",
    "  for j in range(Xb.shape[1]):\n",
    "    ix = Xb[i, j]\n",
    "    dC[ix] += demb[i, j]\n",
    "# -----------------\n",
    "\n",
    "cmp('logprobs', dlogprobs, logprobs)\n",
    "cmp('probs', dprobs, probs)\n",
    "cmp('counts_sum_inv', dcounts_sum_inv, counts_sum_inv)\n",
    "cmp('counts_sum', dcounts_sum, counts_sum)\n",
    "cmp('counts', dcounts, counts)\n",
    "cmp('norm_logits', dnorm_logits, norm_logits)\n",
    "cmp('logit_maxes', dlogit_maxes, logit_maxes)\n",
    "cmp('logits', dlogits, logits)\n",
    "cmp('h', dh, h)\n",
    "cmp('W2', dW2, W2)\n",
    "cmp('b2', db2, b2)\n",
    "cmp('hpreact', dhpreact, hpreact)\n",
    "cmp('bngain', dbngain, bngain)\n",
    "cmp('bnbias', dbnbias, bnbias)\n",
    "cmp('bnraw', dbnraw, bnraw)\n",
    "cmp('bnvar_inv', dbnvar_inv, bnvar_inv)\n",
    "cmp('bnvar', dbnvar, bnvar)\n",
    "cmp('bndiff2', dbndiff2, bndiff2)\n",
    "cmp('bndiff', dbndiff, bndiff)\n",
    "cmp('bnmeani', dbnmeani, bnmeani)\n",
    "cmp('hprebn', dhprebn, hprebn)\n",
    "cmp('embcat', dembcat, embcat)\n",
    "cmp('W1', dW1, W1)\n",
    "cmp('b1', db1, b1)\n",
    "cmp('emb', demb, emb)\n",
    "cmp('C', dC, C)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "3.5571467876434326 diff: -2.384185791015625e-07\n"
     ]
    }
   ],
   "source": [
    "# Exercise 2: backprop through cross_entropy but all in one go\n",
    "# to complete this challenge look at the mathematical expression of the loss,\n",
    "# take the derivative, simplify the expression, and just write it out\n",
    "\n",
    "# forward pass\n",
    "\n",
    "# before:\n",
    "# logit_maxes = logits.max(1, keepdim=True).values\n",
    "# norm_logits = logits - logit_maxes # subtract max for numerical stability\n",
    "# counts = norm_logits.exp()\n",
    "# counts_sum = counts.sum(1, keepdims=True)\n",
    "# counts_sum_inv = counts_sum**-1 # if I use (1.0 / counts_sum) instead then I can't get backprop to be bit exact...\n",
    "# probs = counts * counts_sum_inv\n",
    "# logprobs = probs.log()\n",
    "# loss = -logprobs[range(n), Yb].mean()\n",
    "\n",
    "# now:\n",
    "loss_fast = F.cross_entropy(logits, yb)\n",
    "print(loss_fast.item(), 'diff:', (loss_fast - loss).item())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "logits          | exact: False | approximate: True  | maxdiff: 5.820766091346741e-09\n"
     ]
    }
   ],
   "source": [
    "# Backward Pass \n",
    "dlogits = F.softmax(logits, dim=-1)\n",
    "dlogits[range(n), yb] -= 1\n",
    "dlogits /= n\n",
    "\n",
    "cmp('logits', dlogits, logits) # I can only get approximate to be true, my maxdiff is 6e-9"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAbUAAAH5CAYAAAAGMKDKAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjguMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8g+/7EAAAACXBIWXMAAA9hAAAPYQGoP6dpAAAseElEQVR4nO3dfXDU5bn/8c/uZrMhTxvDQ0IkUEAEKw+dUqU5thQlh4cz488H5oy2zhzsODp6gnMk9djGn4ra9sTS32k9PUPxHw+czhRtnSk6OufgKEqobcCCIrVqCkglFBIwbbJhwyab3fv3h2NqlIfstQkbbt6vmZ2BZK9cd7773f1ks5v7CjjnnAAA8EAw1wsAAGC4EGoAAG8QagAAbxBqAABvEGoAAG8QagAAbxBqAABv5OV6AZ+WTqd15MgRlZSUKBAI5Ho5AIAcc86pu7tbVVVVCgbP/Fxs1IXakSNHVF1dnetlAABGmdbWVk2aNOmM1xl1oVZSUiJJeq75YhUVZ/bb0e/Nu8Lct+6Nd0x1xcGEuWdIts1cLgrZeyZcyFR3Ip1v7lkc7DunPa3HVZIigX5TXTKL3+QHjevtdfa7bzyL29PKGY9RJJA01aVk/02P9RzKpmeR8fsMBtLmntZzr8x4n5akD/pLM67pOZHWLVe9P5APZzLqQu3jXzkWFQdVVJLZnSAvEDb3LSyxPdgXneWp8JlY7zjFIXvPPGerdWl7z2LjMbL2zCbUCgK2nn1ZhJp1vdbbUpKUtp3v2Ugbj1GB8UH7/Au1lKkumMWrNObHoCwe94r67efeUF6S4o0iAABvjFiorVu3Tp/73OdUUFCgBQsW6PXXXx+pVgAASBqhUPvFL36h+vp6rVmzRm+88YbmzZunpUuX6tixYyPRDgAASSMUaj/60Y90++2365vf/KY+//nP64knnlBhYaH+67/+6zPX7e3tVSwWG3QBAMBi2EOtr69Pu3fvVm1t7d+aBIOqra1Vc3PzZ67f2NioaDQ6cOHt/AAAq2EPtQ8//FCpVEoVFRWDPl5RUaG2trbPXL+hoUFdXV0Dl9bW1uFeEgDgApHzt/RHIhFFIpFcLwMA4IFhf6Y2btw4hUIhtbe3D/p4e3u7Kisrh7sdAAADhj3U8vPzNX/+fG3dunXgY+l0Wlu3blVNTc1wtwMAYMCI/Pqxvr5eK1eu1Je+9CVdeeWVevzxxxWPx/XNb35zJNoBACBphELtpptu0vHjx/XQQw+pra1NX/jCF7Rly5bPvHkEAIDhFHDO2TfKGwGxWEzRaFSLL1qpvEBmm65+Z/er5r4h2faXSzj7fpPjQ3FTXU/a3jPubBvZFgXsG5h2pgtNdWXBHnNPK+vtmc2efyXGTbHTWez9aP0+redsLhxPFZlrw8aNreNp+5veqvO6THWt/VFzz6Jgr6kulcUrV0nDpurx7pRumLdfXV1dKi0984bI7P0IAPAGoQYA8AahBgDwBqEGAPAGoQYA8AahBgDwBqEGAPAGoQYA8AahBgDwBqEGAPAGoQYA8AahBgDwBqEGAPAGoQYA8MaIzFMbDne92qyiksxGFPRkMfZhfKjbVNedHmPumTCMYJDso1wkaUa4w1R3PGX/PitDMVOddZxLRxYjR0qNY2CyYRnFIWU35mRsDkbIdBrXW2AcA5MNc88sniaUBW3jr96X7fyRpKSzRcBjs75k7nl/y+sZ14QDQz82PFMDAHiDUAMAeINQAwB4g1ADAHiDUAMAeINQAwB4g1ADAHiDUAMAeINQAwB4g1ADAHiDUAMAeINQAwB4g1ADAHhj1O7SXxY6qaJQZplblMVu3q39ZaY66+7+ktSWKjXVzcn/MIuetp3Sgxnskv1phUHb7fKnZJmprizUY6qTpIQLm+r+7VL7ruX1771lrrX6s/F8D2VxHhQF+kx11p3kK0InTHWS1J4qNtVV5dkfD97vt03fKAvaz/eiQNJUt3bfdnPPtv6SjGsS6dSQr8szNQCANwg1AIA3CDUAgDcINQCANwg1AIA3CDUAgDcINQCANwg1AIA3CDUAgDcINQCANwg1AIA3CDUAgDcINQCANwg1AIA3Ru3ombQLKu0yy9xgwJn7FQZ7TXUpBcw9C4xjH1qNIyokKWX8OaYsmDD3/EuqwFT375fONdU9/v6vTXWS1JO2jeb5Tssb5p59LmSqKwjazh9JKpCttjttuy0l+/iipPH4/CWLtVrvmylnfzyw3jdDso8Dskpm+Nj8SRMMI4FOhIb+PfJMDQDgDUINAOANQg0A4A1CDQDgDUINAOANQg0A4A1CDQDgDUINAOANQg0A4A1CDQDgDUINAOANQg0A4A1CDQDgjVG7S39psFfFwcwyt7U/au5XmddtqjueKjL3tO4EXhDoN/dMONtNns2u3DHjbumBsG2t2ezObmXdgV6SZNyF3rp7vSSFZJtoUZ3Xae6ZzX3Fwvo9SlKfbMe2N4vbJJEOm+q+GOk09zySsq23I2WfFGK5XeLp1JCvyzM1AIA3CDUAgDcINQCAN4Y91B5++GEFAoFBl1mzZg13GwAAPmNE3ihy+eWX6+WXX/5bk7xR+34UAIBHRiRt8vLyVFlZORJfGgCA0xqR19T27dunqqoqTZs2TbfccosOHTp02uv29vYqFosNugAAYDHsobZgwQJt3LhRW7Zs0fr163Xw4EF99atfVXf3qf8OrLGxUdFodOBSXV093EsCAFwghj3Uli9frn/8x3/U3LlztXTpUv3P//yPOjs79ctf/vKU129oaFBXV9fApbW1dbiXBAC4QIz4OzjKysp06aWXav/+/af8fCQSUSQSGellAAAuACP+d2onTpzQgQMHNHHixJFuBQC4wA17qN17771qamrSn/70J/32t7/VDTfcoFAopK9//evD3QoAgEGG/dePhw8f1te//nV1dHRo/Pjx+spXvqIdO3Zo/Pjxw90KAIBBhj3Unn766eH+kgAADMmo3erj/gVXKy+Qn1HN//vDy2e/0ml0pzPr9bFpeV3mnvuSF5nqQkH7SA2rbMZ4hIxjWb7zzu9MddmM5kkaR/OEZB89kx8Y+liNTwpn8X2mjK88tPWXmHsWBG2jlqzCxuMqSYXGsVDneryOJB3utz+Md6bHmOqyObYWwcDQH3/Y0BgA4A1CDQDgDUINAOANQg0A4A1CDQDgDUINAOANQg0A4A1CDQDgDUINAOANQg0A4A1CDQDgDUINAOANQg0A4I1Ru0t/uuek0hnuQp5wIXO/aLDXVPcX4+7+kpQ2/kxRHkqYex4x7rLel8Wu3Nns8G9xPGXfST5o3G0/aJxEIEnjQ3FTXVuq2NyzwLgL/fTwX809rfdP6wSN8aE+U50kvZ8sNdVZpzxI9qkLMRcx9yw0Pu71pO09Ey6ccc3J9NAff3imBgDwBqEGAPAGoQYA8AahBgDwBqEGAPAGoQYA8AahBgDwBqEGAPAGoQYA8AahBgDwBqEGAPAGoQYA8AahBgDwxqjdpX/Vm++oqCSzXb2z2Z19fKjbVJdNz/LQCVPd+8lyc88Z4Q5T3fHUGHPPsuBJU11KAVNdX8o+raE0aJ+AYPWtGV8z1dW/+6a5Z0HQtiN8NpMwOo07uxcYd6+37rQvSSXG8yDsbGuVpGl5Paa6N/rGmXumA7bnNdlMpSgMZD4ZwAWH3o9nagAAbxBqAABvEGoAAG8QagAAbxBqAABvEGoAAG8QagAAbxBqAABvEGoAAG8QagAAbxBqAABvEGoAAG8QagAAbxBqAABvjNrRMwXBfhVkMG5AkhLpsLmfdcxJZShm7lkeSprqeowjPCSpM51vrrXqM/7s1JkqNNVlMz4mli4w1RUEbbelJD36x9+Y6koC9p5W7alic631PhY3nu9jQ3FTnSR1pIpMdaEsRrJ80G8b71SQg/Ogcd5XzbVr3tqWcU1SjJ4BAFyACDUAgDcINQCANwg1AIA3CDUAgDcINQCANwg1AIA3CDUAgDcINQCANwg1AIA3CDUAgDcINQCANwg1AIA3Ru0u/Yl0nkLpUEY16SwyOiRnqmtLlZp79umEqS7h7NMIqoO2qQLHU7YdxCUpP4Mdtj/Jusu6dYd1Kbsd/q0euvQqU139u2+ae1qP7fgsdr7vNO62XxDsN9V1GycuSPbjk3D2h9QpeSdNdW/02c/3okCfqa7hrV+be/Yps8f1j2qGPuGBZ2oAAG8QagAAbxBqAABvZBxq27dv17XXXquqqioFAgE9++yzgz7vnNNDDz2kiRMnasyYMaqtrdW+ffuGa70AAJxWxqEWj8c1b948rVu37pSfX7t2rX7yk5/oiSee0M6dO1VUVKSlS5cqkTj3L8ADAC4sGb9VZ/ny5Vq+fPkpP+ec0+OPP64HHnhA1113nSTpZz/7mSoqKvTss8/q5ptvzm61AACcwbC+pnbw4EG1tbWptrZ24GPRaFQLFixQc3PzKWt6e3sVi8UGXQAAsBjWUGtra5MkVVRUDPp4RUXFwOc+rbGxUdFodOBSXV09nEsCAFxAcv7ux4aGBnV1dQ1cWltbc70kAMB5alhDrbKyUpLU3t4+6OPt7e0Dn/u0SCSi0tLSQRcAACyGNdSmTp2qyspKbd26deBjsVhMO3fuVE1NzXC2AgDgMzJ+9+OJEye0f//+gf8fPHhQe/bsUXl5uSZPnqx77rlH3/ve9zRjxgxNnTpVDz74oKqqqnT99dcP57oBAPiMjENt165duvrqqwf+X19fL0lauXKlNm7cqPvuu0/xeFx33HGHOjs79ZWvfEVbtmxRQYF9c1EAAIYi41BbtGiRnDv9jvaBQECPPvqoHn300awWBgBApkbt6JmiYJ+Kgpm95Nedto9HSWUw2uCTyoI95p7WcTeFwV5zzwPJi0x1ZSH79xlztpEjSZf5iApJqszrNtVJUlt/ialuQsg2RkiSvvPeLlNdNmNgpuTZzvfXEmXmnuGAbYTMo9O+aKp7/E+/NdVJUnuq2FQXDNjGLEnSO8moqc46PkaS+oz3sbFZjGgqCKQyrokEh35cc/6WfgAAhguhBgDwBqEGAPAGoQYA8AahBgDwBqEGAPAGoQYA8AahBgDwBqEGAPAGoQYA8AahBgDwBqEGAPAGoQYA8Mao3aW/ONin4gx36Y+nbbvBS/Yd4RMubO5ZGbTt7B63be4vSSoIJk11PVkcW+vO5QUB61rtt8nFeTFTXdzZ70pVxqkCx1P2qRTxZOY7pUtSSfCkuaf1vvLQ+2+Y6j7ot02kkOzTN2Jp+9zIWfl/NdW90zfW3DMo233zgdlXn/1Kp6v9/a8zromn2aUfAHABItQAAN4g1AAA3iDUAADeINQAAN4g1AAA3iDUAADeINQAAN4g1AAA3iDUAADeINQAAN4g1AAA3iDUAADeINQAAN4YtaNnZoRDKg1nlrlhdZn7dRlHq4w1jqiQpNb+MlNdecg2skaS0s72c0w42Gfu2ZEqMtXlB2zjUVLGUTeS1JG2j3Oxau23jStZO+ML5p737dtjqrOOA5KkkGwzk6wjpazjYySp0DiiqU+2EVaStC8ZNdXNCHeYe7b2l5rqVr/1urlnyvBcKpNHAp6pAQC8QagBALxBqAEAvEGoAQC8QagBALxBqAEAvEGoAQC8QagBALxBqAEAvEGoAQC8QagBALxBqAEAvEGoAQC8MWp36W9ORFQUzmzH67KgbRdwSepMF5rqstkxP6WAqW5+vn0n8D+nbDuXv9M31tyzMhQz1SWc7fTMZqd06xSDH3zhKnPPe/f81lT30P7fmXsmXNhU153FFIPCYK+prsxY15my3acl+30zkbYdV0kKGadLWHfal6SSYMJUF85iEkZn2jaVYqh4pgYA8AahBgDwBqEGAPAGoQYA8AahBgDwBqEGAPAGoQYA8AahBgDwBqEGAPAGoQYA8AahBgDwBqEGAPAGoQYA8AahBgDwxqgdPZNyQaUyHANyPFVi7jc+1G2qy0XPl04Wm3vOCNtGz5QFbXXZeGTOQlPd6rdeN/csNY7iWPPWNnNP65iTeDpi7jk2FDfVFQRS5p6dxvWGZBsplZ/FWgsC/bbCLJ4mTMuz3cfe6Btnb2qcIBPMYvTMSOOZGgDAG4QaAMAbGYfa9u3bde2116qqqkqBQEDPPvvsoM/feuutCgQCgy7Lli0brvUCAHBaGYdaPB7XvHnztG7dutNeZ9myZTp69OjA5amnnspqkQAADEXGbxRZvny5li9ffsbrRCIRVVZWmhcFAIDFiLymtm3bNk2YMEEzZ87UXXfdpY6OjtNet7e3V7FYbNAFAACLYQ+1ZcuW6Wc/+5m2bt2qH/zgB2pqatLy5cuVSp367bWNjY2KRqMDl+rq6uFeEgDgAjHsf6d28803D/x7zpw5mjt3rqZPn65t27Zp8eLFn7l+Q0OD6uvrB/4fi8UINgCAyYi/pX/atGkaN26c9u/ff8rPRyIRlZaWDroAAGAx4qF2+PBhdXR0aOLEiSPdCgBwgcv4148nTpwY9Kzr4MGD2rNnj8rLy1VeXq5HHnlEK1asUGVlpQ4cOKD77rtPl1xyiZYuXTqsCwcA4NMyDrVdu3bp6quvHvj/x6+HrVy5UuvXr9fevXv13//93+rs7FRVVZWWLFmi7373u4pE7PvUAQAwFBmH2qJFi+Tc6TcYffHFF7NaEAAAVqN2l37ge2+/aqrrSBUN80rO7uHL/s5c++C7zcO4EuDCxobGAABvEGoAAG8QagAAbxBqAABvEGoAAG8QagAAbxBqAABvEGoAAG8QagAAbxBqAABvEGoAAG8QagAAbxBqAABvjNpd+seH4ioOZZa5nekx5n7Wnd3Hh7rNPa3rrQzFzD1b+0tNdeNDcXPP48Zj2+dCprqiYK+pTpL6ZOv56HuvmXsmjd9nNu6feqWp7lv7/2DumR9ImepKgn2munn5J0x1ktRmW2pW3u8vNNVV53Wae3an8011CRc290y5zJ9LJdJDv0F4pgYA8AahBgDwBqEGAPAGoQYA8AahBgDwBqEGAPAGoQYA8AahBgDwBqEGAPAGoQYA8AahBgDwBqEGAPAGoQYA8AahBgDwxqgdPdPr8pSX4YiClALmfqFA2lR3PFVi7jkz3GGqu3Pq18w9nzjYZKprT9nH+lhNCNlGh3SkbSM8JKkseNJca5Uw/mxZEEyae1pHyIzNYgSRdbyTdTxKS9qZ6iQp7mw9Q7L3vDKSMNU199qOqySVBG09U1k8H0q6zGPHZXBceaYGAPAGoQYA8AahBgDwBqEGAPAGoQYA8AahBgDwBqEGAPAGoQYA8AahBgDwBqEGAPAGoQYA8AahBgDwBqEGAPDGqN2lPy+QUjiQ2Y7XIWffITsc6DfVpTKcJPBJLcmxprp/O9Bs7nkkZdvBvty4m7ckxVO2Hc+tu+0XBOy71xcEUqa6znTE3DNtPIcSLmzu+fisOaa6+vfeMvdMG3+GLjLenpeG7VM73k3apjVkMylkV5/tfjI+i8kJ1gkIibT93CsxTMIIBoc+RYVnagAAbxBqAABvEGoAAG8QagAAbxBqAABvEGoAAG8QagAAbxBqAABvEGoAAG8QagAAbxBqAABvEGoAAG8QagAAbxBqAABvjNrRMxbW8TGSFM9idIiVZQTDR3X20SqhtG08T1cWx8c6IqUo0GeqKwna6iSpLVVsqstm3I11XMnjl3/B3PP+P+401Y01nrOSdDxVZK612NsXMtd25mDskXWETDZjj6pCPaa6fef48SDhhj4SimdqAABvEGoAAG9kFGqNjY264oorVFJSogkTJuj6669XS0vLoOskEgnV1dVp7NixKi4u1ooVK9Te3j6siwYA4FQyCrWmpibV1dVpx44deumll5RMJrVkyRLF43/7XfDq1av1/PPP65lnnlFTU5OOHDmiG2+8cdgXDgDAp2X0RpEtW7YM+v/GjRs1YcIE7d69WwsXLlRXV5eefPJJbdq0Sddcc40kacOGDbrsssu0Y8cOffnLXx6+lQMA8ClZvabW1dUlSSovL5ck7d69W8lkUrW1tQPXmTVrliZPnqzm5uZTfo3e3l7FYrFBFwAALMyhlk6ndc899+iqq67S7NmzJUltbW3Kz89XWVnZoOtWVFSora3tlF+nsbFR0Wh04FJdXW1dEgDgAmcOtbq6Or399tt6+umns1pAQ0ODurq6Bi6tra1ZfT0AwIXL9MfXq1at0gsvvKDt27dr0qRJAx+vrKxUX1+fOjs7Bz1ba29vV2Vl5Sm/ViQSUSRy7v/wGQDgn4yeqTnntGrVKm3evFmvvPKKpk6dOujz8+fPVzgc1tatWwc+1tLSokOHDqmmpmZ4VgwAwGlk9Eytrq5OmzZt0nPPPaeSkpKB18mi0ajGjBmjaDSq2267TfX19SovL1dpaanuvvtu1dTU8M5HAMCIyyjU1q9fL0latGjRoI9v2LBBt956qyTpxz/+sYLBoFasWKHe3l4tXbpUP/3pT4dlsQAAnElGoebc2TfDLSgo0Lp167Ru3TrzogAAsBi1u/R/7+++prxAfmY1e7ee/UqnkQra3gg6Pmjb5VqSPui/yFRXlMVO4Alnu8kLspiAENbQd9j+pD5n22U9YayTpHDAttZsJgPMDNl2vn/w3VP/7edQpIxvfP5zf6m5Z4FxukRneoytXxb3kzLj/TqWLjD3tE4xsE55kKR9SeNjULDX3NNy7oUD6SFflw2NAQDeINQAAN4g1AAA3iDUAADeINQAAN4g1AAA3iDUAADeINQAAN4g1AAA3iDUAADeINQAAN4g1AAA3iDUAADeINQAAN4YtaNn6n/TrKKSzDK3wzi6QZIuzouZ6o6nC809reMtvjNrobnnd975nbnWakLohKnOOnIk4cKmOknKN47JyWbcTY+zjfWxjsmRpCrjufe+i5p7pt358zO09dims3ie0Hf2cZWnZB3pkyuWkUApRs8AAC5EhBoAwBuEGgDAG4QaAMAbhBoAwBuEGgDAG4QaAMAbhBoAwBuEGgDAG4QaAMAbhBoAwBuEGgDAG4QaAMAbo3aX/pDSynTf83+/5HJzv//7/h5TXTY7j8eMu8lv2LfV3LMlWWqqi6ULzD1LjBMQrEIa+o7ew6U61Guu/X3fRaa6oqC9p/UIJZ39IaPEOK2hJx0x1fVl/AjyN/F0vqnOOnlDkuLO1jObx6CSYMJUVxSwTZaQpE7D7Zl0gSFfl2dqAABvEGoAAG8QagAAbxBqAABvEGoAAG8QagAAbxBqAABvEGoAAG8QagAAbxBqAABvEGoAAG8QagAAbxBqAABvjNpd+iOBfhUEMsvcR97fbe5n3SHbusu1JAWNu/S3pmy7lmcjP5Ay1x5PFZnqCrPYhd6qM1VoqtvXb98Rvixk29m9JJA096ybvshUd897vzf3TBjP92DANlOgLIv7ZkjOVNeWKjb3TBl32x+f123uWWC8Xx9PjTH3tEx6SLqhr5NnagAAbxBqAABvEGoAAG8QagAAbxBqAABvEGoAAG8QagAAbxBqAABvEGoAAG8QagAAbxBqAABvEGoAAG8QagAAbxBqAABvjNrRM+FAWuFAZjXHUyXmflPy/mqqO5JFz7Bx7EPS2cecWMd4jA/axqNI0l/SBaa6nrRtxE5KGZ44n1CUg3E3Vu1ZjDn5/r7fDuNKhqbbeB5YR9ZUh06Y6iSpx3gfs46PyUba2c936+NX4yXzzD1/8P6OjGvygkN/3OKZGgDAG4QaAMAbhBoAwBsZhVpjY6OuuOIKlZSUaMKECbr++uvV0tIy6DqLFi1SIBAYdLnzzjuHddEAAJxKRqHW1NSkuro67dixQy+99JKSyaSWLFmieDw+6Hq33367jh49OnBZu3btsC4aAIBTyejdj1u2bBn0/40bN2rChAnavXu3Fi5cOPDxwsJCVVZWDulr9vb2qrf3b+84i8VimSwJAIABWb2m1tXVJUkqLy8f9PGf//znGjdunGbPnq2Ghgb19Jz+7eCNjY2KRqMDl+rq6myWBAC4gJn/Ti2dTuuee+7RVVddpdmzZw98/Bvf+IamTJmiqqoq7d27V9/+9rfV0tKiX/3qV6f8Og0NDaqvrx/4fywWI9gAACbmUKurq9Pbb7+t1157bdDH77jjjoF/z5kzRxMnTtTixYt14MABTZ8+/TNfJxKJKBKx/ZEtAACfZPr146pVq/TCCy/o1Vdf1aRJk8543QULFkiS9u/fb2kFAMCQZfRMzTmnu+++W5s3b9a2bds0derUs9bs2bNHkjRx4kTTAgEAGKqMQq2urk6bNm3Sc889p5KSErW1tUmSotGoxowZowMHDmjTpk36h3/4B40dO1Z79+7V6tWrtXDhQs2dO3dEvgEAAD6WUaitX79e0kd/YP1JGzZs0K233qr8/Hy9/PLLevzxxxWPx1VdXa0VK1bogQceGLYFAwBwOhn/+vFMqqur1dTUlNWCPpZ0QSUz3PE6HOg39zueLjTVVefZ/67uSL9th+yxoZPmnn/uLzXVdQdsO6VLUtLZ3o9k3Z19fKjbVCfZJz0EZZt+IEkh4+QE65SHbLSlbOePlMXtYjy040P2ISQrqmtMdT882GzueaQ/aqqLOfsb7QqNUykeObDL3DMYOHOOZFvD3o8AAG8QagAAbxBqAABvEGoAAG8QagAAbxBqAABvEGoAAG8QagAAbxBqAABvEGoAAG8QagAAbxBqAABvEGoAAG8QagAAb9hnM4ywwkC/CgOZZW5JFmNgOtP5prqUC5h7pmSrjRtHuUhSkXHURE/aPt7isVlfMtV9/4+vmeoKshjJUhK0jfUpCiTNPfuMP1uGlPkIj491pIrOec/udIGpzjq66LWEbZSLJD10YLeprs04Tkqyj866OHTC3PNPxnE3BVmc752G8yCeHvr8IZ6pAQC8QagBALxBqAEAvEGoAQC8QagBALxBqAEAvEGoAQC8QagBALxBqAEAvEGoAQC8QagBALxBqAEAvEGoAQC8MWp36T/cX6qi/lBGNf9+yeXmfv928HVT3Qf9F5l7lht31/5Lqtjcc0a4w1R3PDXG3POxP/7aVGedYvDn/lJTnSSVBhOmukQWkxOs32dnFpMTxobiprpsJiBY11tmnJxgnQogSQXGHfOzeZowLa/HVPdG3zhzz6JAn6kuNYqfD43elQEAkCFCDQDgDUINAOANQg0A4A1CDQDgDUINAOANQg0A4A1CDQDgDUINAOANQg0A4A1CDQDgDUINAOANQg0A4A1CDQDgjVE7euZCkHBhU13YOhZDUqFtyklWwoG0rdCd+5+5gsa1Jl1mY5I+KWwc5xKyHlfZR+WYb0vZx5VY11oeso1y+ain7fa03paS1J6yPR6E5Mw9rWOPsrlnpg3367Qb+vfIMzUAgDcINQCANwg1AIA3CDUAgDcINQCANwg1AIA3CDUAgDcINQCANwg1AIA3CDUAgDcINQCANwg1AIA3CDUAgDdG7S79oUBaoUBmO0iv3v+uuV9neoy51sq6s/v4UNzc84N+2/cZSxeYe1p3vrfs5i1lt3t9gXGX9fHBXnPP4+mIudbquzO/bKp7sGWHuaf19vxiJGGq29lbZKqTpJKgrWdnqtDcs8h4Do0PdZt7lgX7THV/ThWbe1qmk/RmcO7wTA0A4A1CDQDgDUINAOCNjEJt/fr1mjt3rkpLS1VaWqqamhr97//+78DnE4mE6urqNHbsWBUXF2vFihVqb28f9kUDAHAqGYXapEmT9Nhjj2n37t3atWuXrrnmGl133XX6wx/+IElavXq1nn/+eT3zzDNqamrSkSNHdOONN47IwgEA+LSM3v147bXXDvr/97//fa1fv147duzQpEmT9OSTT2rTpk265pprJEkbNmzQZZddph07dujLX7a92woAgKEyv6aWSqX09NNPKx6Pq6amRrt371YymVRtbe3AdWbNmqXJkyerubn5tF+nt7dXsVhs0AUAAIuMQ+33v/+9iouLFYlEdOedd2rz5s36/Oc/r7a2NuXn56usrGzQ9SsqKtTW1nbar9fY2KhoNDpwqa6uzvibAABAMoTazJkztWfPHu3cuVN33XWXVq5cqXfeece8gIaGBnV1dQ1cWltbzV8LAHBhy3hHkfz8fF1yySWSpPnz5+t3v/ud/uM//kM33XST+vr61NnZOejZWnt7uyorK0/79SKRiCKRc7+jAgDAP1n/nVo6nVZvb6/mz5+vcDisrVu3DnyupaVFhw4dUk1NTbZtAAA4q4yeqTU0NGj58uWaPHmyuru7tWnTJm3btk0vvviiotGobrvtNtXX16u8vFylpaW6++67VVNTwzsfAQDnREahduzYMf3TP/2Tjh49qmg0qrlz5+rFF1/U3//930uSfvzjHysYDGrFihXq7e3V0qVL9dOf/nREFg4AwKdlFGpPPvnkGT9fUFCgdevWad26dVktCgAAi4BzzuV6EZ8Ui8UUjUa1KHij8gKZjSh44mCTue8R48iIqwrsL0tut023UFHANi5CkrqNI2QKsxitYhk1IdlH84SN42Mkqcc4BiYc6Df3DMl2F0wps9FMg2qNY2DSWbwMbz2HrLfJlLy/muok6Xja9nhgPWc/qrVNAvtS5C/mnvuT9pFSVpbHg3h3Sv9n3gF1dXWptLT0jNdlQ2MAgDcINQCANwg1AIA3CDUAgDcINQCANwg1AIA3CDUAgDcINQCANwg1AIA3CDUAgDcINQCANwg1AIA3bDtojqCP91fud8mMa7u70+a+8ZStNpb5Mv/W07ihsQtk8X2mbZv9uqC9Z8LZapPGrbbDWRyfk8bjk5fFJsq52dDY1tMZ1yrZzyHrbXIiL5v7ybk9ZyWp39luz+6+LL7PpL3WKuEyvz17Tny0zqHsvz/qduk/fPiwqqurc70MAMAo09raqkmTJp3xOqMu1NLptI4cOaKSkhIFAp/9ySUWi6m6ulqtra1nHUFwIeL4nBnH58w4PmfG8Tm7kThGzjl1d3erqqpKweCZXzUbdb9+DAaDZ01iSSotLeWkOgOOz5lxfM6M43NmHJ+zG+5jFI1Gh3Q93igCAPAGoQYA8MZ5F2qRSERr1qxRJGIb8e47js+ZcXzOjONzZhyfs8v1MRp1bxQBAMDqvHumBgDA6RBqAABvEGoAAG8QagAAbxBqAABvnFehtm7dOn3uc59TQUGBFixYoNdffz3XSxoVHn74YQUCgUGXWbNm5XpZObV9+3Zde+21qqqqUiAQ0LPPPjvo8845PfTQQ5o4caLGjBmj2tpa7du3LzeLzYGzHZ9bb731M+fUsmXLcrPYHGhsbNQVV1yhkpISTZgwQddff71aWloGXSeRSKiurk5jx45VcXGxVqxYofb29hyt+NwayvFZtGjRZ86hO++8c8TXdt6E2i9+8QvV19drzZo1euONNzRv3jwtXbpUx44dy/XSRoXLL79cR48eHbi89tpruV5STsXjcc2bN0/r1q075efXrl2rn/zkJ3riiSe0c+dOFRUVaenSpUokjKMTzjNnOz6StGzZskHn1FNPPXUOV5hbTU1Nqqur044dO/TSSy8pmUxqyZIlisfjA9dZvXq1nn/+eT3zzDNqamrSkSNHdOONN+Zw1efOUI6PJN1+++2DzqG1a9eO/OLceeLKK690dXV1A/9PpVKuqqrKNTY25nBVo8OaNWvcvHnzcr2MUUuS27x588D/0+m0q6ysdD/84Q8HPtbZ2ekikYh76qmncrDC3Pr08XHOuZUrV7rrrrsuJ+sZjY4dO+YkuaamJufcR+dLOBx2zzzzzMB13n33XSfJNTc352qZOfPp4+Occ1/72tfcv/zLv5zztZwXz9T6+vq0e/du1dbWDnwsGAyqtrZWzc3NOVzZ6LFv3z5VVVVp2rRpuuWWW3To0KFcL2nUOnjwoNra2gadT9FoVAsWLOB8+oRt27ZpwoQJmjlzpu666y51dHTkekk509XVJUkqLy+XJO3evVvJZHLQOTRr1ixNnjz5gjyHPn18Pvbzn/9c48aN0+zZs9XQ0KCenp4RX8uo26X/VD788EOlUilVVFQM+nhFRYXee++9HK1q9FiwYIE2btyomTNn6ujRo3rkkUf01a9+VW+//bZKSkpyvbxRp62tTZJOeT59/LkL3bJly3TjjTdq6tSpOnDggO6//34tX75czc3NCoVCuV7eOZVOp3XPPffoqquu0uzZsyV9dA7l5+errKxs0HUvxHPoVMdHkr7xjW9oypQpqqqq0t69e/Xtb39bLS0t+tWvfjWi6zkvQg1ntnz58oF/z507VwsWLNCUKVP0y1/+UrfddlsOV4bz1c033zzw7zlz5mju3LmaPn26tm3bpsWLF+dwZedeXV2d3n777Qv+derTOd3xueOOOwb+PWfOHE2cOFGLFy/WgQMHNH369BFbz3nx68dx48YpFAp95p1F7e3tqqyszNGqRq+ysjJdeuml2r9/f66XMip9fM5wPg3dtGnTNG7cuAvunFq1apVeeOEFvfrqq4PmPFZWVqqvr0+dnZ2Drn+hnUOnOz6nsmDBAkka8XPovAi1/Px8zZ8/X1u3bh34WDqd1tatW1VTU5PDlY1OJ06c0IEDBzRx4sRcL2VUmjp1qiorKwedT7FYTDt37uR8Oo3Dhw+ro6PjgjmnnHNatWqVNm/erFdeeUVTp04d9Pn58+crHA4POodaWlp06NChC+IcOtvxOZU9e/ZI0sifQ+f8rSlGTz/9tItEIm7jxo3unXfecXfccYcrKytzbW1tuV5azn3rW99y27ZtcwcPHnS/+c1vXG1trRs3bpw7duxYrpeWM93d3e7NN990b775ppPkfvSjH7k333zTffDBB8455x577DFXVlbmnnvuObd371533XXXualTp7qTJ0/meOXnxpmOT3d3t7v33ntdc3OzO3jwoHv55ZfdF7/4RTdjxgyXSCRyvfRz4q677nLRaNRt27bNHT16dODS09MzcJ0777zTTZ482b3yyitu165drqamxtXU1ORw1efO2Y7P/v373aOPPup27drlDh486J577jk3bdo0t3DhwhFf23kTas4595//+Z9u8uTJLj8/31155ZVux44duV7SqHDTTTe5iRMnuvz8fHfxxRe7m266ye3fvz/Xy8qpV1991Un6zGXlypXOuY/e1v/ggw+6iooKF4lE3OLFi11LS0tuF30Onen49PT0uCVLlrjx48e7cDjspkyZ4m6//fYL6gfIUx0bSW7Dhg0D1zl58qT753/+Z3fRRRe5wsJCd8MNN7ijR4/mbtHn0NmOz6FDh9zChQtdeXm5i0Qi7pJLLnH/+q//6rq6ukZ8bcxTAwB447x4TQ0AgKEg1AAA3iDUAADeINQAAN4g1AAA3iDUAADeINQAAN4g1AAA3iDUAADeINQAAN4g1AAA3vj/YiKFJc8epbUAAAAASUVORK5CYII=",
      "text/plain": [
       "<Figure size 800x600 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.imshow(dlogits.detach())\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "max diff: tensor(4.7684e-07, grad_fn=<MaxBackward1>)\n"
     ]
    }
   ],
   "source": [
    "# Exercise 3: backprop through batchnorm but all in one go\n",
    "# to complete this challenge look at the mathematical expression of the output of batchnorm,\n",
    "# take the derivative w.r.t. its input, simplify the expression, and just write it out\n",
    "# BatchNorm paper: https://arxiv.org/abs/1502.03167\n",
    "\n",
    "# forward pass\n",
    "\n",
    "# before:\n",
    "# bnmeani = 1/n*hprebn.sum(0, keepdim=True)\n",
    "# bndiff = hprebn - bnmeani\n",
    "# bndiff2 = bndiff**2\n",
    "# bnvar = 1/(n-1)*(bndiff2).sum(0, keepdim=True) # note: Bessel's correction (dividing by n-1, not n)\n",
    "# bnvar_inv = (bnvar + 1e-5)**-0.5\n",
    "# bnraw = bndiff * bnvar_inv\n",
    "# hpreact = bngain * bnraw + bnbias\n",
    "\n",
    "# now:\n",
    "hpreact_fast = bngain * (hprebn - hprebn.mean(0, keepdim=True)) / torch.sqrt(hprebn.var(0, keepdim=True, unbiased=True) + 1e-5) + bnbias\n",
    "print('max diff:', (hpreact_fast - hpreact).abs().max())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# backward pass\n",
    "\n",
    "# before we had:\n",
    "# dbnraw = bngain * dhpreact\n",
    "# dbndiff = bnvar_inv * dbnraw\n",
    "# dbnvar_inv = (bndiff * dbnraw).sum(0, keepdim=True)\n",
    "# dbnvar = (-0.5*(bnvar + 1e-5)**-1.5) * dbnvar_inv\n",
    "# dbndiff2 = (1.0/(n-1))*torch.ones_like(bndiff2) * dbnvar\n",
    "# dbndiff += (2*bndiff) * dbndiff2\n",
    "# dhprebn = dbndiff.clone()\n",
    "# dbnmeani = (-dbndiff).sum(0)\n",
    "# dhprebn += 1.0/n * (torch.ones_like(hprebn) * dbnmeani)\n",
    "\n",
    "# calculate dhprebn given dhpreact (i.e. backprop through the batchnorm)\n",
    "# (you'll also need to use some of the variables from the forward pass up above)\n",
    "\n",
    "# -----------------\n",
    "# YOUR CODE HERE :)\n",
    "dhprebn = None # TODO. my solution is 1 (long) line\n",
    "# -----------------\n",
    "\n",
    "cmp('hprebn', dhprebn, hprebn) # I can only get approximate to be true, my maxdiff is 9e-10"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 35,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "12297\n",
      "      0/ 200000: 3.8354\n",
      "  10000/ 200000: 2.1844\n",
      "  20000/ 200000: 2.4262\n",
      "  30000/ 200000: 2.4899\n",
      "  40000/ 200000: 2.0123\n",
      "  50000/ 200000: 2.3421\n",
      "  60000/ 200000: 2.4065\n",
      "  70000/ 200000: 2.1004\n",
      "  80000/ 200000: 2.3460\n",
      "  90000/ 200000: 2.1176\n",
      " 100000/ 200000: 1.9322\n",
      " 110000/ 200000: 2.3346\n",
      " 120000/ 200000: 2.0428\n",
      " 130000/ 200000: 2.4166\n",
      " 140000/ 200000: 2.2850\n",
      " 150000/ 200000: 2.1924\n",
      " 160000/ 200000: 1.9430\n",
      " 170000/ 200000: 1.8824\n",
      " 180000/ 200000: 2.0226\n",
      " 190000/ 200000: 1.8984\n"
     ]
    }
   ],
   "source": [
    "# Exercise 4: putting it all together!\n",
    "# Train the MLP neural net with your own backward pass\n",
    "\n",
    "# init\n",
    "n_embd = 10 # the dimensionality of the character embedding vectors\n",
    "n_hidden = 200 # the number of neurons in the hidden layer of the MLP\n",
    "\n",
    "g = torch.Generator().manual_seed(2147483647) # for reproducibility\n",
    "C  = torch.randn((vocab_size, n_embd),            generator=g)\n",
    "# Layer 1\n",
    "W1 = torch.randn((n_embd * block_size, n_hidden), generator=g) * (5/3)/((n_embd * block_size)**0.5)\n",
    "b1 = torch.randn(n_hidden,                        generator=g) * 0.1\n",
    "# Layer 2\n",
    "W2 = torch.randn((n_hidden, vocab_size),          generator=g) * 0.1\n",
    "b2 = torch.randn(vocab_size,                      generator=g) * 0.1\n",
    "# BatchNorm parameters\n",
    "bngain = torch.randn((1, n_hidden))*0.1 + 1.0\n",
    "bnbias = torch.randn((1, n_hidden))*0.1\n",
    "\n",
    "parameters = [C, W1, b1, W2, b2, bngain, bnbias]\n",
    "print(sum(p.nelement() for p in parameters)) # number of parameters in total\n",
    "for p in parameters:\n",
    "  p.requires_grad = True\n",
    "\n",
    "# same optimization as last time\n",
    "max_steps = 200000\n",
    "batch_size = 32\n",
    "n = batch_size # convenience\n",
    "lossi = []\n",
    "\n",
    "# use this context manager for efficiency once your backward pass is written (TODO)\n",
    "with torch.no_grad():\n",
    "\n",
    "  # kick off optimization\n",
    "  for i in range(max_steps):\n",
    "    # minibatch construct\n",
    "    ix = torch.randint(0, Xtr.shape[0], (batch_size,), generator=g)\n",
    "    Xb, Yb = Xtr[ix], ytr[ix] # batch X,Y\n",
    "\n",
    "    # forward pass\n",
    "    emb = C[Xb] # embed the characters into vectors\n",
    "    embcat = emb.view(emb.shape[0], -1) # concatenate the vectors\n",
    "    # Linear layer\n",
    "    hprebn = embcat @ W1 + b1 # hidden layer pre-activation\n",
    "    # BatchNorm layer\n",
    "    # -------------------------------------------------------------\n",
    "    bnmean = hprebn.mean(0, keepdim=True)\n",
    "    bnvar = hprebn.var(0, keepdim=True, unbiased=True)\n",
    "    bnvar_inv = (bnvar + 1e-5)**-0.5\n",
    "    bnraw = (hprebn - bnmean) * bnvar_inv\n",
    "    hpreact = bngain * bnraw + bnbias\n",
    "    # -------------------------------------------------------------\n",
    "    # Non-linearity\n",
    "    h = torch.tanh(hpreact) # hidden layer\n",
    "    logits = h @ W2 + b2 # output layer\n",
    "    loss = F.cross_entropy(logits, Yb) # loss function\n",
    "\n",
    "    # backward pass\n",
    "    for p in parameters:\n",
    "      p.grad = None\n",
    "    # loss.backward() # use  this for correctness comparisons, delete it later!\n",
    "\n",
    "    # manual backprop! #swole_doge_meme\n",
    "    # -----------------\n",
    "    # YOUR CODE HERE :) \n",
    "    # dC, dW1, db1, dW2, db2, dbngain, dbnbias = None, None, None, None, None, None, None\n",
    "    # grads have same shape as original tensor. Consider gradients of only indexed tensor. Others are 0\n",
    "    dlogits = F.softmax(logits, dim=-1)\n",
    "    dlogits[range(n), Yb] -= 1\n",
    "    dlogits /= n\n",
    "\n",
    "    dh = dlogits@ W2.T\n",
    "    dW2 = h.T @ dlogits\n",
    "    db2 = dlogits.sum(0)\n",
    "\n",
    "    dhpreact = (1. - h**2) * dh\n",
    "\n",
    "    dbngain = (bnraw * dhpreact).sum(0, keepdim=True)\n",
    "    dbnbias = dhpreact.sum(0, keepdim=True)\n",
    "    dhprebn = bngain*bnvar_inv/n * (n*dhpreact - dhpreact.sum(0) - n/(n-1)*bnraw*(dhpreact*bnraw).sum(0))\n",
    "\n",
    "    dembcat = dhprebn @ W1.T \n",
    "    dW1 = embcat.T @ dhprebn\n",
    "    db1 = dhprebn.sum(0)\n",
    "\n",
    "    demb = dembcat.view(emb.shape)\n",
    "    # INDEXING\n",
    "    dC = torch.zeros_like(C)\n",
    "    for k in range(Xb.shape[0]):\n",
    "      for j in range(Xb.shape[1]):\n",
    "        ix = Xb[k, j]\n",
    "        dC[ix] += demb[k, j]\n",
    "    grads = [dC, dW1, db1, dW2, db2, dbngain, dbnbias]\n",
    "\n",
    "    # -----------------\n",
    "    # -----------------\n",
    "  # YOUR CODE HERE :)\n",
    "\n",
    "\n",
    "    # -----------------\n",
    "\n",
    "    # update\n",
    "    lr = 0.1 if i < 100000 else 0.01 # step learning rate decay\n",
    "\n",
    "    for p, grad in zip(parameters, grads):\n",
    "      # p.data += -lr * p.grad # old way of cheems doge (using PyTorch grad from .backward())\n",
    "      p.data += -lr * grad # new way of swole doge TODO: enable\n",
    "\n",
    "    # track stats\n",
    "    if i % 10000 == 0: # print every once in a while\n",
    "      print(f'{i:7d}/{max_steps:7d}: {loss.item():.4f}')\n",
    "    lossi.append(loss.log10().item())\n",
    "\n",
    "    # if i >= 100: # TODO: delete early breaking when you're ready to train the full net\n",
    "    #   break"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(27, 10)        | exact: False | approximate: True  | maxdiff: 1.1175870895385742e-08\n",
      "(30, 200)       | exact: False | approximate: True  | maxdiff: 7.450580596923828e-09\n",
      "(200,)          | exact: False | approximate: True  | maxdiff: 4.6566128730773926e-09\n",
      "(200, 27)       | exact: False | approximate: True  | maxdiff: 1.1175870895385742e-08\n",
      "(27,)           | exact: False | approximate: True  | maxdiff: 7.450580596923828e-09\n",
      "(1, 200)        | exact: False | approximate: True  | maxdiff: 1.862645149230957e-09\n",
      "(1, 200)        | exact: False | approximate: True  | maxdiff: 3.725290298461914e-09\n"
     ]
    }
   ],
   "source": [
    "# useful for checking your gradients\n",
    "for p,g in zip(parameters, grads):\n",
    "  cmp(str(tuple(p.shape)), g, p)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 36,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "train 2.101194143295288\n",
      "val 2.1394879817962646\n"
     ]
    }
   ],
   "source": [
    "# evaluate train and val loss\n",
    "\n",
    "@torch.no_grad() # this decorator disables gradient tracking\n",
    "def split_loss(split):\n",
    "  x,y = {\n",
    "    'train': (Xtr, ytr),\n",
    "    'val': (Xva, yva),\n",
    "    'test': (Xte, yte),\n",
    "  }[split]\n",
    "  emb = C[x] # (N, block_size, n_embd)\n",
    "  embcat = emb.view(emb.shape[0], -1) # concat into (N, block_size * n_embd)\n",
    "  hpreact = embcat @ W1 + b1\n",
    "  hpreact = bngain * (hpreact - bnmean) * (bnvar + 1e-5)**-0.5 + bnbias\n",
    "  h = torch.tanh(hpreact) # (N, n_hidden)\n",
    "  logits = h @ W2 + b2 # (N, vocab_size)\n",
    "  loss = F.cross_entropy(logits, y)\n",
    "  print(split, loss.item())\n",
    "\n",
    "split_loss('train')\n",
    "split_loss('val')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# I achieved:\n",
    "# train 2.0718822479248047\n",
    "# val 2.1162495613098145"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# sample from the model\n",
    "g = torch.Generator().manual_seed(2147483647 + 10)\n",
    "\n",
    "for _ in range(20):\n",
    "\n",
    "    out = []\n",
    "    context = [0] * block_size # initialize with all ...\n",
    "    while True:\n",
    "      # forward pass\n",
    "      emb = C[torch.tensor([context])] # (1,block_size,d)\n",
    "      embcat = emb.view(emb.shape[0], -1) # concat into (N, block_size * n_embd)\n",
    "      hpreact = embcat @ W1 + b1\n",
    "      hpreact = bngain * (hpreact - bnmean) * (bnvar + 1e-5)**-0.5 + bnbias\n",
    "      h = torch.tanh(hpreact) # (N, n_hidden)\n",
    "      logits = h @ W2 + b2 # (N, vocab_size)\n",
    "      # sample\n",
    "      probs = F.softmax(logits, dim=1)\n",
    "      ix = torch.multinomial(probs, num_samples=1, generator=g).item()\n",
    "      context = context[1:] + [ix]\n",
    "      out.append(ix)\n",
    "      if ix == 0:\n",
    "        break\n",
    "\n",
    "    print(''.join(itos[i] for i in out))"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "torch",
   "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.13"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}