Etude1.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import torch\n",
    "import numpy as np"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "my_device = torch.device(\"cuda:0\" if torch.cuda.is_available() else \"cpu\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Etude 1\n",
    "Let $u$ be a vector of $2D$ points. We need to get a matrix $A$ whose element $(i,j)$ depends on euclidean distance between $u_i$ and $u_j$.\n",
    "\n",
    "$$u \\in \\mathbb{R}^{N \\times 2}$$\n",
    "$$A \\in \\mathbb{R}^{N \\times N}, \\ A_{ij} = f(\\|u_i - u_j\\|) $$\n",
    "\n",
    "Particularly, let's take $f(x) = exp(-\\frac{x}{4})$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Auxiliary stuff for performance measurement"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "import time\n",
    "import matplotlib.pyplot as plt\n",
    "%matplotlib inline\n",
    "\n",
    "\n",
    "def bench(func, n_range, repeat=3, device=None, warmup=0):\n",
    "    \"\"\"\n",
    "    Run 'func' with vector sizes in 'n_range'.\n",
    "    For each vector size average time of 'repeat' runs.\n",
    "    \"\"\"\n",
    "    if warmup:\n",
    "        for i in range(warmup):\n",
    "            size = np.random.choice(n_range)\n",
    "            func(torch.randn(size, 2, device=device))\n",
    "    t = []\n",
    "    for n in n_range:\n",
    "        t0 = []\n",
    "        for i in range(repeat):\n",
    "            u = torch.randn(n, 2, device=device)\n",
    "            start = time.time()\n",
    "            func(u)\n",
    "            end = time.time()\n",
    "            t0.append(end-start)\n",
    "        t.append(sum(t0)/len(t0))\n",
    "    return t\n",
    "\n",
    "\n",
    "def print_interval(ts):\n",
    "    \"\"\"\n",
    "    Print (mean ± std. dev.) interval \n",
    "    for the input set of execution times.\n",
    "    \"\"\"\n",
    "    m, dev = np.mean(ts), np.std(ts)\n",
    "    unit, factor = 's', 1\n",
    "    if m < 1e-3:\n",
    "        unit, factor = 'µs', 10**6\n",
    "    elif m < 1:\n",
    "        unit, factor = 'ms', 10**3\n",
    "    m *= factor\n",
    "    dev *= factor\n",
    "    print(f\"{m:.2f} ± {dev:.2f} {unit} per run (mean ± std. dev. of {len(ts)} runs)\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Naive approach"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "def naive(u):\n",
    "    N = u.shape[0]\n",
    "    result = torch.zeros(N, N)\n",
    "    for i in range(N):\n",
    "        for j in range(N):\n",
    "            result[i,j] = torch.exp(-0.25 * torch.dist(u[i], u[j]))\n",
    "    return result"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Vectorized approaches"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "def fully_vectorized(u):\n",
    "    u_ = u.reshape(1, *u.shape) # add one dimension for transposing\n",
    "    A_ = u_.permute((1, 0, 2)) - u_ # row - column: broadcasting subtraction\n",
    "    return A_.pow_(2).sum(dim=2).sqrt_().mul_(-0.25).exp_()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "def one_dim_vectorized(u):\n",
    "    N = u.shape[0]\n",
    "    result = torch.zeros(N, N)\n",
    "    for i in range(N):\n",
    "        du = u[i] - u\n",
    "        result[i,:] = du.pow_(2).sum(dim=1).sqrt_().mul_(-0.25).exp_()\n",
    "    return result"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Comparison"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### fully vectorized vs. naive"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "N = 10**3\n",
    "u = torch.randn(N, 2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "CPU times: user 12.1 s, sys: 392 ms, total: 12.5 s\n",
      "Wall time: 12.5 s\n"
     ]
    }
   ],
   "source": [
    "%%time\n",
    "A = naive(u)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "CPU times: user 64 ms, sys: 0 ns, total: 64 ms\n",
      "Wall time: 36.9 ms\n"
     ]
    }
   ],
   "source": [
    "%%time\n",
    "B = fully_vectorized(u)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "assert torch.equal(A, B)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 1-dim vectorized vs. fully vectorized"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "N = 10**4\n",
    "u = torch.randn(N, 2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "2.34 s ± 184 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)\n"
     ]
    }
   ],
   "source": [
    "%%timeit\n",
    "A = fully_vectorized(u)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "2.97 s ± 189 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)\n"
     ]
    }
   ],
   "source": [
    "%%timeit\n",
    "B = one_dim_vectorized(u)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [],
   "source": [
    "assert torch.equal(A, B)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [],
   "source": [
    "ns = list(range(10**3, 6*10**3, 10**3))\n",
    "t_fv = bench(func=fully_vectorized, n_range=ns)\n",
    "t_1dv = bench(func=one_dim_vectorized, n_range=ns)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.legend.Legend at 0x7f089017cda0>"
      ]
     },
     "execution_count": 16,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYUAAAEKCAYAAAD9xUlFAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4wLCBo\ndHRwOi8vbWF0cGxvdGxpYi5vcmcvpW3flQAAIABJREFUeJzt3Xd8VfX9x/HXh0AIOxDCymDJ3pAw\n1IoLxQValaVWQEvdba22+nOitVU7tLa4FbGCAa0DcaAWrRNIwpIwJAbIhIRgAgESMj6/P84lxpiQ\nADk59yaf5+ORB/eee3LvO4fkfu73fL/n+xVVxRhjjAFo4nUAY4wx/sOKgjHGmHJWFIwxxpSzomCM\nMaacFQVjjDHlrCgYY4wpZ0XBGGNMOSsKxhhjyllRMMYYU66p1wGOVceOHbVHjx5exzDGmICSmJi4\nR1XDa9ov4IpCjx49SEhI8DqGMcYEFBHZWZv97PSRMcaYclYUjDHGlLOiYIwxppyrfQoiMhH4BxAE\nPK+qD1d6PBpYAIT69rlDVd871tcpLi4mPT2dwsLCOkjtP0JCQoiMjKRZs2ZeRzHGNBKuFQURCQLm\nAROAdCBeRJaq6qYKu90NLFHVp0RkIPAe0ONYXys9PZ02bdrQo0cPRKQO0ntPVcnNzSU9PZ2ePXt6\nHccY00i4efpoNJCsqimqehiIAyZX2keBtr7b7YDM43mhwsJCwsLCGkxBABARwsLCGlzrxxjj39w8\nfRQBpFW4nw6MqbTP/cCHInIz0Ao4+3hfrCEVhCMa4s9kjPFvbrYUqnpHq7z253TgJVWNBM4H/i0i\nP8kkInNEJEFEEnJyclyIaowx/qu0THno3U1k5B1y/bXcLArpQFSF+5H89PTQNcASAFX9GggBOlZ+\nIlV9VlVjVDUmPLzGC/I8MXv2bDp16sTgwYOrfHzHjh3VPmaMMdVRVea+k8Rzn2/nf1vd/1DsZlGI\nB/qISE8RCQamAUsr7ZMKnAUgIgNwikJANgVmzpzJBx984HUMY0wD88xnKbz89U5++bOezBgT7frr\nuVYUVLUEuAlYDmzGGWWUJCIPiMgk326/A34pIuuBV4GZqlr5FFNAOO200+jQocOPtiUmJjJs2DDG\njRvHvHnzyrePGTOGpKSk8vunn346iYmJ9ZbVGBMY3l6XwcPvb+GiYd2487wB9fKarl6n4Lvm4L1K\n2+6tcHsTcEpdvubcd5LYlLmvLp+Sgd3act9Fg475+2bNmsU///lPxo8fz+23316+fdq0aSxZsoS5\nc+eSlZVFZmYmo0aNqsvIxpgA91XyHm57bT1jenbgr5cPpUmT+hl4Ylc0uyQ/P5+8vDzGjx8PwFVX\nXVX+2JQpU3jttdcAWLJkCZdffrknGY0x/mlz1j5+9e9EenZsxbO/iKF506B6e+2AmyW1Jsfzid4N\nqlrtkNKIiAjCwsLYsGEDixcv5plnnqnndMYYf5WZd4hZ8+Np1bwpL80aTbsW9TujgbUUXBIaGkq7\ndu344osvAFi4cOGPHp82bRqPPvoo+fn5DBkyxIuIxhg/k3+omJnzV3OgqISXZsfSLbRFvWewolBH\npk+fzrhx49i6dSuRkZG88MILzJ8/nxtvvJFx48bRosWP/3Mvu+wy4uLimDJlikeJjTH+pKiklDkv\nJ7B9zwGeuWoU/bu0rfmbXNDgTh955dVXX61y+/r168tv33///eW3O3fuTElJiduxjDEBoKxM+d2S\n9azavpd/TBvOySf95HKtemMtBWOM8djDH2xh2YYs7jivP5OHR3iaxYqCMcZ4aP6X23n2sxR+Ma47\nvzqtl9dxrCgYY4xX3v8miweWbeLcQZ2576JBfjEJphUFY4zxQPyOvfx68TpGRIXyj2kjCKqni9Nq\nYkXBGGPqWXL2fq5dkEBkaAteuDqWkGb1d3FaTawoGGNMPcreV8jVL8bTLKgJC2aPpn2rYK8j/YgV\nhTryxBNPMGDAAK644gqvoxhj/FRBUQkz58fz/cHDzJ8ZS1SHll5H+gm7TqGOPPnkk7z//vu2nrIx\npkrFpWVc/0oiW3fv5/mrYxgS2c7rSFWylkIduO6660hJSWHSpEm0a9eOvLy88sdOOukkdu/e7WE6\nY4zXVJU//GcDn2/bw59/PoQz+nXyOlK1Gl5L4f07YNc3dfucXYbAeQ9X+/DTTz/NBx98wCeffMKD\nDz7Im2++yaxZs1i1ahU9evSgc+fOdZvHGBNQ/v7Rt7yxJoPfnt2XKTFRNX+Dh6ylUMemTp3K4sWL\nAYiLi2Pq1KkeJzLGeGnhqp38c0Uy02KjuOWsk7yOU6OG11I4yif6+jBu3DiSk5PJycnhrbfe4u67\n7/Y0jzHGOx9v2s09b23kjH7h/PHiwX5xcVpNXG0piMhEEdkqIskickcVjz8mIut8X9+KSF5VzxNI\nRIRLLrmEW2+9lQEDBhAWFuZ1JGOMB9al5XHTq2sYHNGOf80YSdOgwDgx41pLQUSCgHnABCAdiBeR\npb4lOAFQ1d9W2P9mYIRbeerT1KlTiY2N5aWXXvI6ijHGAzv2HGD2S/F0ahPCC1fH0qp54JyUcTPp\naCBZVVMARCQOmAxsqmb/6cB9LuZx1Y4dO8pvx8TEoKrehTHGeGZPQRFXz1+NqvLSrFjC2zT3OtIx\ncbM9EwGkVbif7tv2EyLSHegJrKjm8TkikiAiCTk5OXUe1Bhj6sLBwyVc81I8u/cV8sLMWHqFt/Y6\n0jFzsyhU1aNS3cfnacDrqlpa1YOq+qyqxqhqTHh4eJ0FNMaYulJSWsbNi9byTUY+/5w+kpHR7b2O\ndFzcLArpQMUBuZFAZjX7TgOqXrqslhri6ZqG+DMZ0xCpKve8ncR/t2Qzd/JgJgwM3GuT3CwK8UAf\nEekpIsE4b/xLK+8kIv2A9sDXx/tCISEh5ObmNqg3UVUlNzeXkJAQr6MYY2ow75NkXl2dyg2n9+aq\nsd29jnNCXOtoVtUSEbkJWA4EAS+qapKIPAAkqOqRAjEdiNMTeEePjIwkPT2dhtbfEBISQmRkpNcx\njDFH8XpiOn/98Ft+PiKC28/t53WcEyaB9uk6JiZGExISvI5hjDF89m0Os1+KZ2yvMF6cGUtwU/+9\nFkFEElU1pqb9/PcnMMYYP7YxI5/rX0mkT+c2PHXlSL8uCMeiYfwUxhhTj9L2HmTWS/G0a9GMl2bF\n0iakmdeR6owVBWOMOQZ5Bw8zc/5qiopLWTB7NJ3bNqzBIIFz7bUxxnissLiUaxckkLb3EP++ZjR9\nOrfxOlKds6JgjDG1UFqm/HbxOhJ2fs+8GSMZ06thTnZpp4+MMaYGqsqDyzbx/sZd3H3BAC4Y2tXr\nSK6xomCMMTV47vMUXvpqB9ec2pNrf9bL6ziusqJgjDFHsXR9Jn96bwsXDO3KXecP8DqO66woGGNM\nNb7+LpfblqxndM8O/O3yYTRp4v8rp50oKwrGGFOFLbv2MeffCXQPa8lzV8UQ0izI60j1woqCMcZU\nkpV/iFnz42nRLIiXZo+mXcuGc3FaTWxIqjHGVLCvsJhZ8+PZX1jCkl+NIyK0hdeR6pW1FIwxxqeo\npJRfvZxIcnYBT185ioHd2nodqd5ZS8EYY4CyMuX21zbwdUouf58yjFP7dPQ6kiespWCMMcAjy7ew\ndH0mt5/bj5+PbLzrmFhRMMY0egu+2sEz/0vhyrHR3HB6b6/jeMrVoiAiE0Vkq4gki8gd1ewzRUQ2\niUiSiCxyM48xxlT2wcZd3P9OEhMGdmbupMGINPxrEY7GtT4FEQkC5gETgHQgXkSWquqmCvv0Ae4E\nTlHV70Wkk1t5jDGmssSde/l13FqGR4XyxLQRBDWCi9Nq4mZLYTSQrKopqnoYiAMmV9rnl8A8Vf0e\nQFWzXcxjjDHlvssp4JoFCXQLbcELV8fSIrhxXJxWEzeLQgSQVuF+um9bRX2BviLypYisFJGJLuYx\nxhgAsvcXcvWLq2naRFgwazQdWgV7HclvuDkktap2mFbx+n2A04FI4HMRGayqeT96IpE5wByA6Ojo\nuk9qjGk0CopKmP1SPLkFh4mbM5bosJZeR/IrbrYU0oGoCvcjgcwq9nlbVYtVdTuwFadI/IiqPquq\nMaoaEx4e7lpgY0zDVlxaxg0L17A5az9PXjGSYVGhXkfyO24WhXigj4j0FJFgYBqwtNI+bwFnAIhI\nR5zTSSkuZjLGNFKqyv+98Q2ffZvDQxcP5oz+Nq6lKq4VBVUtAW4ClgObgSWqmiQiD4jIJN9uy4Fc\nEdkEfALcrqq5bmUyxjRej328jdcS07nlrD5MG22noasjqpVP8/u3mJgYTUhI8DqGMSaAvLo6lTvf\n+IYpMZE8cunQRnktgogkqmpMTfvZFc3GmAZtxZbd3P3WRsb3DeehS4Y0yoJwLKwoGGMarPVpedy4\ncC0DurbhyStG0izI3vJqYkfIGNMg7cw9wOyX4glrHcyLM2Np1TyAJ4U+uBfeugFyvnX9pQL4KBlj\nTNVyC4q4+sXVlKqyYPZoOrUJ8TrS8dvyHiz7DRzMhe6nQHhfV1/OioIxpkE5dLiUaxYkkJVfyKJf\njqF3eGuvIx2fQ9/D+3fAhjjoPBiueB26DnX9Za0oGGMajJLSMm5+dS3r0/N46opRjOrewetIx+fb\n5fDOr6EgG077PZx2OzStn6k4rCgYYxoEVeX+d5L4ePNuHpg8iImDu3gd6dgdyoPld8G6V6DTQJj+\nKnQbUa8RrCgYYxqEJz/9jldWpvKr8b34xbgeXsc5dts+hqU3Q8Eu+NnvYPwfoGnzeo9hRcEYE/De\nWJPOX5Zv5eLh3fjDuf29jnNsCvfBh3fBmpchvD9MewUiRnkWx4qCMSagfb4th9+/voGTe4fx6GXD\naBJIC+V8twLevhn2Z8Ipv4HT74Rm3o6UsqJgjAlYSZn5XP/KGk7q1JqnrxpFcNMAufSqaD98eA8k\nzoewPnDNRxBZ4wwU9cKKgjEmIKV/f5BZ8+NpE9KU+bNiaRvSzOtItZPyP3j7JshPg5NvhjPugmYt\nvE5VzoqCMSbg5B8sZub8eA4Vl/L6dSfTtZ3/vKlWq6gAPr4P4p+HDr1h9nKIHuN1qp+womCMCSiF\nxaX88uUEUnMPsmD2aPp1aeN1pJrt+MKZpiIvFcbeCGfeDcH+ueKbFQVjTMAoK1N+t2Q9q3fs5Ynp\nIxjXO8zrSEd3+AD89wFY9TS07wmz3oPuJ3ud6qisKBhjAsZD723m3W+yuOv8AUwa1s3rOEe38yun\ndfD9dhhzHZx1LwS38jpVjawoGGMCwvOfp/DCF9uZeXIPrv1ZT6/jVO/wQVjxR1j5JIRGw8x3ocep\nXqeqNVfHb4nIRBHZKiLJInJHFY/PFJEcEVnn+7rWzTzGmMC0bEMmf3x3M+cN7sI9Fw7034VyUlfB\n06fCynkQey1c/1VAFQRwsaUgIkHAPGACkA7Ei8hSVd1UadfFqnqTWzmMMYFtVUouty5eT2yP9jw2\ndThB/nhxWvEh+OQh+Opf0C4KfrEUeo33OtVxcfP00WggWVVTAEQkDpgMVC4KxhhTpW937+eXLycQ\n1aEFz/0ihpBmQV5H+qn0BHjzOsjdBjGzYcID0DwARkRVw83TRxFAWoX76b5tlV0qIhtE5HURiXIx\njzEmgOzKL2Tmi6tp3iyIl2aNJrRl/UwdXWvFhfDRffDCBKelcNWbcOFjAV0QwN2iUFUbTyvdfwfo\noapDgY+BBVU+kcgcEUkQkYScnJw6jmmM8Tf7CouZOX81+YeKmT8zlqgOfjamPyMRnh0PXz4OI66E\nG76G3md6napOuFkU0oGKn/wjgcyKO6hqrqoW+e4+B1Q5NaCqPquqMaoaEx4e7kpYY4x/OFxSxvWv\nJJKcXcBTV45icEQ7ryP9oKTIue7g+QnO7KZX/gcm/RNC2nqdrM642acQD/QRkZ5ABjANmFFxBxHp\nqqpZvruTgM0u5jHG+DlV5Q//2cCXybn89fJhnNbXjz4EZq6Dt66H7E0w/Eo49yFoEep1qjrnWlFQ\n1RIRuQlYDgQBL6pqkog8ACSo6lLgFhGZBJQAe4GZbuUxxvi/R5dv5c21Gdx2Tl8uGxXpdRxHyWH4\n7C/w+d+gdSeY8Rr0PcfrVK4R1cqn+f1bTEyMJiQkeB3DGFPH/v31Du55O4npo6P50yWD/eNahKwN\nzlXJu7+BYdNh4p+hRXuvUx0XEUlU1Rrn57Yrmo0xnvswaRf3LU3i7AGdeHDyIO8LQmmx0zL47C/Q\nMgymx0G/87zNVE+sKBhjPJW483tufnUtQyJDeWL6CJoGebxQzq6NTt/Brg0wZAqc9wi07OBtpnpk\nRcEY45mUnAKuXRBPl3YhvHB1DC2DPXxLKi2BLx+DTx9xOpCnLoQBF3qXxyNWFIwxnsjZX8TV81fT\nRIQFs0bTsXVz78Ls3uS0DrLWweBL4by/QCs/n5bbJVYUjDH1Liv/EHNeTiRnfxFxc8bRo6NHU0qX\nlsBXT8Cnf4bmbWHKyzBwsjdZ/IQVBWNMvSkqKeX5z7fzrxXJlKny5BUjGR7l0Vj/nK1O6yAjEQZe\nDBf8DVp19CaLH6mxKIjIWCBJVff77rcBBqrqKrfDGWMajk+3ZjP3nU1s33OAcwd15u4LBnozfUVZ\nKXz9L1jxkLPozWXzYfDP6z+Hn6pNS+EpYGSF+weq2GaMMVVK23uQB5Zt4qNNu+nVsRULZo9mvFdX\nKu/Z5rQO0uOh/4XOBHatO3mTxU/VpiiIVrjCTVXLRMROOxljjurQ4VKe+t93PP2/72jaRLjjvP7M\nPqUnwU09GHJaVgorn4IVD0KzFnDpC06HstfXQ/ih2ry5p4jILTitA4AbgBT3IhljApmqsjxpNw8u\n20RG3iEmDevG/50/gC7tQrwJlPudc1Vy2krodz5c+Di06exNlgBQm6JwHfAEcDfO1Nf/Bea4GcoY\nE5i+yyng/qVJfL5tD/06t+HVX45lXG+PhnaWlcHqZ+DjudA0GC55FoZOsdZBDWosCqqajTPDqTHG\nVOlAUQlPrNjGi19sJ6RpEPddNJCrxnb37urkvSnw9k2w80voO9FpHbTt6k2WAFOb0Ud9cU4ddVbV\nwSIyFJikqn90PZ0xxq+pKkvXZ/Kn9zaze18Rl4+K5PcT+xPexqML0crKIOEF+OheaNIMLn7KmcjO\nWge1VpvTR88BtwPPAKjqBhFZBFhRMKYR27JrH/e+ncTq7XsZHNGWp64cxchoD2cQ/X6H0zrY8Tmc\ndDZc9AS0q2oFYHM0tSkKLVV1daVZC0tcymOM8XP5h4p57KNv+ffKnbQJacqfLhnC1Ngogpp49Gm8\nrAwSX4QP74UmQTDpX84SmdY6OC61KQp7RKQ3vvWVReQyIOvo32KMaWjKypTX16TzyPtb2HvwMFeM\nieZ3E/rRvlWwd6HyUp3Wwfb/Qa8znKUxQ6Nq/j5TrdoUhRuBZ4H+IpIBbAeurM2Ti8hE4B84K689\nr6oPV7PfZcBrQKyq2go6xviZDel53Pt2EuvS8hgZHcqC2aO9XTtZFdYsgOW+QZEXPg6jZlrroA7U\nZvRRCnC2iLQCmhyZ7qImIhIEzAMmAOlAvIgsVdVNlfZrA9wC2LQZxviZvQcO85flW4mLTyWsVXP+\ndvkwLhkRQROvThUB5KfD0pvhuxXQczxM/heERnuXp4GpzeijXwPzgf3AcyIyErhDVT+s4VtHA8m+\nooKIxAGTgU2V9nsQeBS47RizG2NcUlqmLFqdyl+Xb6WgqITZp/Tk12f3oW1IM+9CqcLaV2D5/zlX\nKF/wd4iZba2DOlab00ezVfUfInIu0AmYhVMkaioKEUBahfvpwJiKO4jICCBKVZeJiBUFY/xA4s69\n3PNWEpuy9jGuVxhzJw+ib+c23obalwlLb4Hkj6DHz5zWQfse3mZqoGo195Hv3/OB+aq6Xmq3gGpV\n+5TPoSQiTYDHgJk1PpHIHHxXUUdHWzPRGDdk7y/k4fe38MaaDLq2C+FfM0ZwwZCu3q6XrArrX4X3\n74CyYmfxm9hroYnHS3Y2YLUpCoki8iHQE7jT1wdQVovvSwcqDgOIBDIr3G8DDAY+9f3SdQGWisik\nyp3NqvosTmc3MTExijGmzhSXlrHgqx08/vE2ikpKueH03tx4xkm0au7xvJf7smDZb+DbDyD6ZLh4\nHnTo5W2mRqA2/+vXAMOBFFU9KCJhOKeQahIP9BGRnkAGzlQZM448qKr5QPmKFiLyKXCbjT4ypv58\nlbyH+5YmsS27gPF9w7nvooH0Cm/tbShV2LAE3v89lBTBxIdh9K+sdVBPajP6qAxYU+F+LpBbi+8r\nEZGbgOU4Q1JfVNUkEXkASFDVpccf2xhzIjLzDvHQe5t5d0MWUR1a8NwvYjh7QCdvTxUB7N8Ny34L\nW9+FqLFw8ZMQ1tvbTI2Mq+1DVX0PeK/Stnur2fd0N7MYY366HOZvz+7Lr8b3IqRZkLfBVGHjf+C9\n26D4EJzzEIy93rlC2dQrWyzHmEbiky3ZzH0niR25B71dDrOyghx497ew+R2IjHUmsevYx+tUjZYV\nBWMauNTcgzywLImPN2d7vxxmZUlvwru/g6ICmPAAjLvJWgceO+aiICKbfTfnqeq/6jiPMaaO+NVy\nmJXlZzgXoW16C7qNhEuehvB+XqcyHEdRUNUBvhFIY13IY4w5Qc5ymLt4cNlm/1gOs6Ki/fDF4/D1\nPNAyOOs+OPkWCLKTFv6iNtNcPKKqf6i0+fdVbDPGeCw5u4C57/ywHGbcnLGM7eXRcpgVlZbA2pfh\nkz/BgRwYcjmcda/NWeSHalOeJwCVC8B5VWwzxnikoKiEf/53Gy98sZ0WwUHcf9FArvRyOcwjVGHb\nh/DhPbBnK3Q/BWYshohR3uYy1aq2KIjI9cANQC8R2VDhoTbAl24HM8bU7MhymA+9u5ns/c5ymH84\nrz8dW3u0HGZFWRvgw7udtQ7CToJpi6Df+TaBnZ87WkthEfA+8Gfgjgrb96vqXldTGWNqtDlrH/ct\ndZbDHBLRjqev8ng5zCPyM2DFH505i1q0d+YripkFQR7OsGpqrdqi4JuGIh+YXn9xjDE1ObIc5stf\n76Bti2beL4d5ROVO5FN+DT+7FUI8XIzHHDPr8jcmQJSVKa8npvPIB1v4/uBhZvjDcphgncgNjBUF\nYwLAhvQ87nk7ifVpeYzq3p4FkzxeDhOsE7mBsqJgjB9zlsPcQlx8GmGtmvP3Kc5ymJ5PXJe13teJ\n/Jl1IjcwVhSM8UOlZcqiVTv564ffUlBUwjWn9OQWr5fDhB93IrfsAOf/FUbNtE7kBsSKgjF+JmHH\nXu5921kO8+TeYdw/yQ+Ww7RO5EbDioIxfiJ7n285zLXOcpjzZozk/CFdvD1VZJ3IjY4VBWM8VnE5\nzMMlZdx4hrMcZstgD/88rRO50bKiYIyHvvQth5mcXcDp/cK576JB9OzYyttQ1oncqLlaFERkIvAP\nnOU4n1fVhys9fh1wI1AKFABzVHWTm5mM8QcZeYf407ubefcbZznM538Rw1leL4eZnwErHoT1cdaJ\n3Ii5VhREJAiYhzOhXjoQLyJLK73pL1LVp337TwL+Dkx0K5MxXqu8HOatE/oy5zSPl8Ms70T+l3Pa\nyDqRGzU3WwqjgWRVTQEQkThgMlBeFFR1X4X9WwHqYh5jPLViy24eeGcTO3IPMnFQF+66YIC3y2Fa\nJ7KpgptFIQJIq3A/HRhTeScRuRG4FQgGznQxjzGe2Jl7gAfe2cR/t2TTK7wVL88ezWleLodpncjm\nKNwsClWdHP1JS0BV5wHzRGQGcDdw9U+eSGQOMAcgOto+xZjAcOhwKU9+mswzn6XQtIlw53n9meX1\ncpjWiWxq4GZRSAeiKtyPBDKPsn8c8FRVD6jqs8CzADExMXaKyfg1VeWDjbv447vOcpiTh3fjzvM8\nXg7TOpFNLblZFOKBPiLSE8gApgEzKu4gIn1UdZvv7gXANowJYMnZ+7l/6Sa+SN5D/y5tWDxnLGO8\nXA7TOpHNMXKtKKhqiYjcBCzHGZL6oqomicgDQIKqLgVuEpGzgWLge6o4dWRMIMjeV8hzn6cw/8sd\n/rEcZmkJrFkAn/7Z14k8Bc66xzqRTY1cvU5BVd8D3qu07d4Kt3/t5usb46aS0jI+3ZpDXHwan2zN\nprRMmRITye8nergcpnUimxNkVzQbc4x27DnAkoQ0Xk9MJ3t/ER1bN+eXP+vFlJhIeoW39i6YdSKb\nOmBFwZhaKCwu5YONu4iLT2Vlyl6aCJzRrxNTY6M4o38nmnl1mgisE9nUKSsKxhxFUmY+i+PTeGtt\nBvsKS4ju0JLbz+3HpSMjvR1NBNaJbFxhRcGYSvYVFvP2ukyWxKfxTUY+wU2bcN7gLkyNjWJszzCa\nNPH4dIx1IhsXWVEwBufagtXb97I4IY33vsmisLiMAV3bMnfSICYP70Zoy2CvIzqtgW+Xw0f3VuhE\nXgIRI71OZhoQKwqmUcveX8gbazJYEp9Gyp4DtGnelEtHRjItNprBEW29Xwv5COtENvXEioJpdEpK\ny/hsWw5xq9P47xZnKOnoHh244YyTOH9IF28Xt6nMOpFNPfOj335j3JWae5AlCWm8lpjG7n1FdGwd\nzLWn9mRKbBS9vRxKWhXrRDYesaJgGrTC4lKWJ+1icXwaX32XSxOB8X3DmTspmrMGeDyUtCrWiWw8\nZkXBNEibs/axOD6NN9dmkH+omMj2LfjdhL5cFhNJ13YtvI73U9aJbPyEFQXTYOwvLGbpemco6fr0\nfIKDmnDu4C5Mi41iXC8/GEpaHetENn7EioIJaKpKws7vWRyfxrsbsjhUXEr/Lm2476KBXDw8gvat\n/GAoaXXy02HFH60T2fgVKwomIO0pKOKNNenExaeRknOAVsFBXDwigmmxUQyNbOc/Q0mrUrQfvngM\nvp5nncjG71hRMAGjtEz5bFsOi1en8fHm3ZSUKTHd23PdZb25YEhXWjX3819n60Q2AcDP/4qMgbS9\nB3ktIY3XEtPJyi8krFUws07pwdTYKE7q1MbreDWzTmQTQKwoGL9UVFLKh0m7WRyfxpff7QHgtD7h\n3HvhQM4a0NnbdY6PhXUimwBjRcH4la279hMXn8qbazPIO1hMRGgLfnOWM5Q0ItQPh5JWxzqRTYBy\ntSiIyETgHzjLcT6vqg9XevyG/UE6AAAUPElEQVRW4FqgBMgBZqvqTjczGf9TUFTCO+szWRyfxrq0\nPJoFCecMcoaSntK7o/8OJa1K4T748nHrRDYBy7WiICJBwDxgApAOxIvIUlXdVGG3tUCMqh4UkeuB\nR4GpbmUy/kNVWZPqDCVdtiGLg4dL6du5NfdcOJBLRkTQwZ+Hklal+BCsWwifPmydyCagudlSGA0k\nq2oKgIjEAZOB8qKgqp9U2H8lcKWLeYwfyC0o4s21GcTFp5GcXUDL4CAmDevG1NgohkeF+vdQ0spU\nISMR1r4CG9+AonzrRDYBz82iEAGkVbifDow5yv7XAO9X9YCIzAHmAERH2yevQFNapny+LYclCWl8\ntGk3xaXKyOhQHrl0CBcM7UZrfx9KWtm+LNiwGNYtckYTNW0BAyfDiCugx8+sE9kENDf/Gqv6y9Aq\ndxS5EogBxlf1uKo+CzwLEBMTU+VzGP+T/v1BXktI57WENDLzC2nfshlXj3OGkvbpHABDSSsqKYKt\n7zmFIPlj0DKIHgeT/gkDL4aQtl4nNKZOuFkU0oGoCvcjgczKO4nI2cBdwHhVLXIxj6kHRSWlfLwp\nm7j4VL5IdoaSnnpSR+66YCBnD+xE86ZBHic8BqqQudYpBN+8BoV50DYCTr0Vhs+AsN5eJzSmzrlZ\nFOKBPiLSE8gApgEzKu4gIiOAZ4CJqprtYhbjsm937y+flXTvgcN0axfCLWf24fKYSCLbt/Q63rEp\nyP7h9FD2JmgaAv0vdE4P9RwPTQKosBlzjFwrCqpaIiI3ActxhqS+qKpJIvIAkKCqS4G/AK2B13wd\njKmqOsmtTKZuHSgqYdkGZyjpmlRnKOmEgZ2ZGhvNqSd1JCiQhpKWHIZty2HtQtj2IWgpRMbChY/B\noJ9Di1CvExpTL0Q1sE7Rx8TEaEJCgtcxGi1VZW1aHotXp7FsQyYHDpdyUqfWTIuN4pIREYS1bu51\nxGOTtcF3emgJHMyF1l1g2DTn9FB4P6/TGVNnRCRRVWNq2i/Ahn0Yr+w9cJg312awOD6Vb3cX0KJZ\nEBcN68rU2GhGRgfYUNIDe5w+grULYfc3EBQM/S+A4VdArzMgyP4sTONlv/2mWmVlyhfJe1ickMZH\nSbs5XFrG8KhQHv75EC4cFmBDSUuLYdtHzgVm3y6HsmLoNsKZfmLwpc5UFMYYKwrmp7L3FbI4Po24\n+DQy8g4R2rIZV47tztTYKPp1CbChpLs3OYVgw2LnSuNW4TDmV06roPNAr9MZ43esKBjAaRV8+d0e\nFq5MLV+r4JSTwrjjvP6cM6hzYA0lPbgXvnndKQZZ66BJM+g30SkEJ51tk9IZcxRWFBq5PQVFvJ6Y\nzqurU9mZe5AOrYK55tSeTB8dTY+OrbyOV3ulJfDdClj3Cmx9H0oPQ5ehMPERGHI5tArzOqExAcGK\nQiOkqqxM2cui1al8sDGL4lJldM8O3DqhLxMHdwmsVkHOVqdFsH4xFOyClmEQe60zeqjLEK/TGRNw\nrCg0InkHD/N6YjqLVqeSknOAtiFNuXJsd64YEx0YK5gdcSgPNv7HKQYZiSBB0Pdc5/RQn3OgaYDN\nsGqMH7Gi0MAdmaJ64cpUln2TxeGSMkZGh/LXy4dx4dCuhDQLkFZBWSmkfOJcU7B5GZQWQaeBcM5D\nMHQKtO7kdUJjGgQrCg3UvsJi3lqbwcKVqWzdvZ/WzZsyNSaKGWOiGdA1gCZv25PsOz0UB/szoUV7\nGHW1c3qo63CbkdSYOmZFoYHZkJ7HwpWpLF2fyaHiUoZEtOPhnw/homHdaBUo1xUU7oOkN51ikLYK\npIkzamjin6HfedA0wK6aNiaABMi7hDmaA0UlLF2fycJVO9mYsY8WzYKYPLwbM8ZEMzQyQObsKSuD\nHZ85p4c2LYWSQ9CxH5w915l2ok0XrxMa0yhYUQhgmzL3sWj1Tt5am0lBUQn9u7ThwcmDmDwigrYh\nATIWf28KrHsV1r8K+WnQvB0Mnw7Dr3RWL7PTQ8bUKysKAebQ4VKWbchk0epU1qbm0bxpEy4Y2pUr\nxnQPnDmIigpg01tOq2Dnl4BA7zPh7PudOYiatfA4oDGNlxWFALFt934WrkrljTXp7CssoXd4K+65\ncCCXjowgtGUADMEsK4PUr5xJ6Da9DcUHoENvOOteGDoN2kV4ndAYgxUFv1ZUUsoHG3excGUqq3fs\npVmQcN7grswYE82Ynh0Co1Xw/U5n5NC6hZC3E4LbwJBLndNDUaPt9JAxfsaKgh/avucAr65O5bWE\nNL4/WEz3sJbceV5/LhsVGRjrFRw+AJvfgbWvwI7PAYGep8EZd8GAiyA4wFZiM6YRcbUoiMhE4B84\nK689r6oPV3r8NOBxYCgwTVVfdzOPPysuLeOjTbtZuGonXybnEtREOGdgZ64Y052Te4fRxN9XMVOF\n1JVOiyDpLTi8H9r3cArBsGkQGu11QmNMLbhWFEQkCJgHTADSgXgRWaqqmyrslgrMBG5zK4e/S9t7\nkLj4VBbHp7OnoIiI0Bbcdk5fpsRE0altiNfxapaf7owcWrfIGUnUrBUMusS5uKz7yXZ6yJgA42ZL\nYTSQrKopACISB0wGyouCqu7wPVbmYg6/U1Jaxoot2Sxancr/vs1BgDP7d+KKMd05rW+4/69tXHzI\nmWpi3UJI+RRQ6H4qnHY7DJgEzVt7ndAYc5zcLAoRQFqF++nAGBdfz+9l5R8ibnUai+PT2LWvkM5t\nm3PzmX2YFhtFt1A/H4apCukJztTUG9+Aon3QLhrG/x6GTYcOPb1OaIypA24Whao+7upxPZHIHGAO\nQHR0YJ2bLi1TPtuWw8KVqazYshsFTusTztzJgzirfyeaBjXxOuLR7cuCDXHO6aE930LTFjBwsnN6\nqMfPoImf5zfGHBM3i0I6EFXhfiSQeTxPpKrPAs8CxMTEHFdhqW/Z+wt5LSGdRatSycg7RMfWwVw3\nvjfTR0cT1cGPR9+UlsDujU6ncfJHzsI1WgbR42DSP2HgxRASQBPqGWOOiZtFIR7oIyI9gQxgGjDD\nxdfzXFmZ8nVKLgtX7eTDJGdJy5N7h/F/5w9gwsDOBDf1w0/VRQWQkeAUgdSvnVNEhwucx0Kj4dTf\nOusUhPX2Nqcxpl64VhRUtUREbgKW4wxJfVFVk0TkASBBVZeKSCzwJtAeuEhE5qrqILcyuWXvgcO8\nnpjGolWp7Mg9SGjLZsw6pQfTR0fTK9zPOl337/IVgJWQthKyNoCWAgKdBzv9A9Fjna92kV6nNcbU\nM1ENiLMx5WJiYjQhIcHrGKgqq7c7S1q+/80uDpeWEdujPVeM6c7EwV38Y/EaVacf4EgRSP0avt/u\nPNa0BUTGOG/+UWMhKhZC2nmb1xjjGhFJVNWYmvazK5qPUf7BYv6zxlnSMjm7gDYhTZkxJpoZY6Lp\n29njJS1LiiBrvfPmf6QQHNrrPNayo1MAYq9x+ge6DLVlK40xP2FFoRacJS3zWLQqlWUbMikqKWN4\nVCiPXjaUi4Z2o0WwR62CQ99DWvwPRSAj0VmmEpzJ5vqd7zsVNM7pE7ALyYwxNbCicBT7C4t5a10m\nC1fuZMuu/bQKDuKyUZHMGBPNoG71fKpF1Vlv4MhpoNSVkL0ZUGjSFLoOg9G/9J0OGmNrFhtjjosV\nhSp8k57PotU7eXtdJgcPlzKoW1v+dMkQJg3vRuv6WtKyrBR2J/3QIZy6EvZlOI8Ft3FmGB30c4ge\nAxGjILhV/eQyxjRoVhR8Dh4uYek6Z/GaDen5hDRrwqRh3ZgxpjvDItu5P0314YO+oaGrnJZA2mpn\nUjmANt2g+zinQzh6LHQeBE38oCPbGNPgNPqisDlrH4tWpfLW2gz2F5XQt3Nr5k4axMUjImjXwsUl\nLQtyfmgBpH7tdBCXlQACnQbC0MudvoDosdAuyvoDjDH1olEWhcLiUt7dkMXCVTtZk5pHcNMmXDCk\nK1eMiWZU9/Z13ypQhdzvKowK+hr2fuc8FtTcOf1z8i1OEYiKhRbt6/b1jTGmlhpVUUjOLmDRqlT+\nsyad/EPF9OrYirsvGMClIyNp36oOh2eWHIZdG348NPTgHuexFu2dN/9RVzv/dh0GTQNg4RxjTKPQ\naIrCk58m8+gHW2kWJJw7qAszxkQzrldY3bQKCvN/GBqatsqZKqLkkPNY+57Q5xynQzh6HIT1sUnk\njDF+q9EUhdP6hCMIl8dE0vFEl7TMT69wlfBKZwI5FCQIug6FmFnOsNDosdCmS53kN8aY+tBoisLg\niHYMjjiOawvKyiBn849PBeX7lokIbu1MFXH6HU4BiIixBWaMMQGt0RSFWis+BBlrfigCaauhKN95\nrHUX581/3E2+oaGDIcgOoTGm4bB3tAO5vqGhviKQuQ7Kip3HwvvD4Et+GBoa2t2GhhpjGrTGVRRU\nncXlK14lvOdb57GgYOg2Esbd6BsaOhpadvA2rzHG1LPGUxQSF8CKP8KBbOd+SKjz6X/YdKcIdBsB\nzUK8zWiMMR5rPEWhTVfofcYPs4Z27GdDQ40xppLGUxT6nuN8GWOMqZarH5VFZKKIbBWRZBG5o4rH\nm4vIYt/jq0Skh5t5jDHGHJ1rRUFEgoB5wHnAQGC6iAystNs1wPeqehLwGPCIW3mMMcbUzM2Wwmgg\nWVVTVPUwEAdMrrTPZGCB7/brwFni+hzVxhhjquNmUYgA0ircT/dtq3IfVS0B8oGwyk8kInNEJEFE\nEnJyclyKa4wxxs2iUNUnfj2OfVDVZ1U1RlVjwsPD6yScMcaYn3KzKKQDURXuRwKZ1e0jIk2BdsBe\nFzMZY4w5CjeLQjzQR0R6ikgwMA1YWmmfpcDVvtuXAStU9SctBWOMMfXDtesUVLVERG4ClgNBwIuq\nmiQiDwAJqroUeAH4t4gk47QQprmVxxhjTM0k0D6Yi0gOsPM4v70jsKcO49QVy3VsLNex89dsluvY\nnEiu7qpaY6dswBWFEyEiCaoa43WOyizXsbFcx85fs1muY1MfuWzyH2OMMeWsKBhjjCnX2IrCs14H\nqIblOjaW69j5azbLdWxcz9Wo+hSMMcYcXWNrKRhjjDmKgC8KIvKiiGSLyMYK2zqIyEciss33b3vf\ndhGRJ3xTdW8QkZEVvudq3/7bROTqql6rDnLdLyIZIrLO93V+hcfu9OXaKiLnVth+1OnHjyNXlIh8\nIiKbRSRJRH7t2+7pMTtKLk+PmYiEiMhqEVnvyzXXt72nb7r3bb7p34N926udDr66vHWc6yUR2V7h\neA33ba+3333fcwaJyFoRWea77+nxqiaTvxyrHSLyjS9Dgm+bd3+PqhrQX8BpwEhgY4VtjwJ3+G7f\nATziu30+8D7OnEtjgVW+7R2AFN+/7X2327uQ637gtir2HQisB5oDPYHvcC74C/Ld7gUE+/YZeIK5\nugIjfbfbAN/6Xt/TY3aUXJ4eM9/P3dp3uxmwyncclgDTfNufBq733b4BeNp3exqw+Gh5Xcj1EnBZ\nFfvX2+++73lvBRYBy3z3PT1e1WTyl2O1A+hYaZtnf48B31JQ1c/46XxJFafkXgBcXGH7y+pYCYSK\nSFfgXOAjVd2rqt8DHwETXchVnclAnKoWqep2IBln6vHaTD9+rLmyVHWN7/Z+YDPObLWeHrOj5KpO\nvRwz389d4LvbzPelwJk4073DT49XVdPBV5e3rnNVp95+90UkErgAeN53X/D4eFXOVIN6O1Y1ZPDk\n7zHgi0I1OqtqFjhvNkAn3/bqpvOuzTTfdeUmX7PvxSNNQq9y+ZrqI3A+ZfrNMauUCzw+Zr7TDuuA\nbJw/tu+APHWme6/8GtVNB+96LlU9crwe8h2vx0SkeeVclV7fjf/Hx4HfA2W++2F4f7wqZzrC62MF\nTjH/UEQSRWSOb5tnf48NtShUp7qpums1hXcdeAroDQwHsoC/eZVLRFoD/wF+o6r7jrZrfWarIpfn\nx0xVS1V1OM5Mv6OBAUd5Dc9yichg4E6gPxCLcyrhD/WZS0QuBLJVNbHi5qO8huu5qskEHh+rCk5R\n1ZE4q1TeKCKnHWVf17M11KKw29ekwvdvtm97ddN512aa7xOmqrt9f8hlwHP80Byu11wi0gznjXeh\nqr7h2+z5Masql78cM1+WPOBTnHO5oeJM9175NaqbDr4+ck30nYZTVS0C5lP/x+sUYJKI7MA5dXcm\nzqd0L4/XTzKJyCt+cKwAUNVM37/ZwJu+HN79PR5PR4S/fQE9+HGH7l/4cSfNo77bF/DjTprV+kMn\nzXacDpr2vtsdXMjVtcLt3+KcMwUYxI871VJwOkyb+m735IdO00EnmEmAl4HHK2339JgdJZenxwwI\nB0J9t1sAnwMXAq/x447TG3y3b+THHadLjpbXhVxdKxzPx4GHvfjd9z336fzQqevp8aomk+fHCmgF\ntKlw+yucvgDP/h5P+D/e6y/gVZzTCsU41fIanHOS/wW2+f7tUOE/fx7OOeFvgJgKzzMbpzMrGZjl\nUq5/+153A85aEhXf8O7y5doKnFdh+/k4I3G+A+6qg1yn4jQrNwDrfF/ne33MjpLL02MGDAXW+l5/\nI3Cvb3svYLXvZ38NaO7bHuK7n+x7vFdNees41wrf8doIvMIPI5Tq7Xe/wvOezg9vwJ4er2oyeX6s\nfMdlve8r6cjvKx7+PdoVzcYYY8o11D4FY4wxx8GKgjHGmHJWFIwxxpSzomCMMaacFQVjjDHlrCgY\nc4JEREXkbxXu3yYi93sYyZjjZkXBmBNXBPxcRDp6HcSYE2VFwZgTV4KzTOJvvQ5izImyomBM3ZgH\nXCEi7bwOYsyJsKJgTB1QZ0bXl4FbvM5izImwomBM3XkcZ46rVl4HMeZ4WVEwpo6o6l6cZSev8TqL\nMcfLioIxdetvgI1CMgHLZkk1xhhTzloKxhhjyllRMMYYU86KgjHGmHJWFIwxxpSzomCMMaacFQVj\njDHlrCgYY4wpZ0XBGGNMuf8HrErHoEgowvEAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7f08e053f4a8>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.plot(ns, t_1dv, label='1dv')\n",
    "plt.plot(ns, t_fv, label='fv')\n",
    "plt.ylabel('t, sec')\n",
    "plt.xlabel('N')\n",
    "plt.legend()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### cpu vs. gpu\n",
    "_Note_: don't use `timeit` for GPU because it disables gc and causes memory error"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [],
   "source": [
    "N = 10**4\n",
    "ns = [N] * 20"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Fully vectorized on CPU"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "2.37 ± 0.23 s per run (mean ± std. dev. of 20 runs)\n"
     ]
    }
   ],
   "source": [
    "ts = bench(fully_vectorized, n_range=ns, repeat=1, warmup=3)\n",
    "print_interval(ts)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 1-dim vectorized on GPU"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1.49 ± 0.11 s per run (mean ± std. dev. of 20 runs)\n"
     ]
    }
   ],
   "source": [
    "ts = bench(one_dim_vectorized, n_range=ns, repeat=1, device=my_device, warmup=3)\n",
    "print_interval(ts)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Fully vectorized on GPU"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "99.11 ± 3.32 µs per run (mean ± std. dev. of 20 runs)\n"
     ]
    }
   ],
   "source": [
    "ts = bench(fully_vectorized, n_range=ns, repeat=1, device=my_device, warmup=3)\n",
    "print_interval(ts)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}