dtw.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 74,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "import seaborn\n",
    "import numpy, scipy, scipy.spatial, matplotlib.pyplot as plt\n",
    "plt.rcParams['figure.figsize'] = (14, 3)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "[&larr; Back to Index](index.html)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Dynamic Time Warping "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In MIR, we often want to compare two sequences of different lengths. For example, we may want to compute a similarity measure between two versions of the same song. These two signals, $x$ and $y$, may have similar sequences of chord progressions and instrumentations, but there may be timing deviations between the two. Even if we were to express the two audio signals using the same feature space (e.g. chroma or MFCCs), we cannot simply sum their pairwise distances because the signals have different lengths.\n",
    "\n",
    "As another example, you might want to align two different performances of the same musical work, e.g. so you can hop from one performance to another at any moment in the work. This problem is known as **music synchronization** (FMP, p. 115)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Dynamic time warping (DTW)** ([Wikipedia](https://en.wikipedia.org/wiki/Dynamic_time_warping); FMP, p. 131) is an algorithm used to align two sequences of similar content but possibly different lengths. \n",
    "\n",
    "Given two sequences, $x[n], n \\in \\{0, ..., N_x - 1\\}$, and $y[n], n \\in \\{0, ..., N_y - 1\\}$, DTW produces a set of index coordinate pairs $\\{ (i, j) ... \\}$ such that $x[i]$ and $y[j]$ are similar.\n",
    "\n",
    "We will use the same approach described in the notebooks [Dynamic Programming](dp.html) and [Longest Common Subsequence](lcs.html)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Example"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Create two arrays, $x$ and $y$, of lengths $N_x$ and $N_y$, respectively."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 123,
   "metadata": {},
   "outputs": [],
   "source": [
    "x = [0, 4, 4, 0, -4, -4, 0]\n",
    "y = [1, 3, 4, 3, 1, -1, -2, -1, 0]\n",
    "nx = len(x)\n",
    "ny = len(y)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 127,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.legend.Legend at 0x11ce5d910>"
      ]
     },
     "execution_count": 127,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAy4AAADBCAYAAAApSeRhAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzs3Xd8VGXa//HPtMxMKi2h9zIB0gMK2CsgRSkqooKw4uqq\nuO1ZH3d119/q7j7bFRVdERCQovRiWRVR7EB6Ahl6DYSQQtpMMuX8/piQEJQWcjJnJtf79doX6+TM\n5M53TiZzzbmv+9YpioIQQgghhBBCaJne3wMQQgghhBBCiIuRwkUIIYQQQgiheVK4CCGEEEIIITRP\nChchhBBCCCGE5knhIoQQQgghhNA8KVyEEEIIIYQQmmdsqW9UVFShmXWX27YNpbS02t/DCGqSsfok\nY/VJxuqSfNUnGatPMlafZKwureUbHR2hO9/XWuUVF6PR4O8hBD3JWH2SsfokY3VJvuqTjNUnGatP\nMlZXIOXbKgsXIYQQQgghRGCRwkUIIYQQQgiheVfU42Kz2WKANOA2u92e3zxDEkIIIYQQQojGmnzF\nxWazmYD/AI7mG44QQgghhBBC/NCVXHH5B/AG8EwzjUUEAbfHy7ovD3D0VBW1tW5/DyeohFWU0v1w\nPt0P59P12F5cfbsTcd1w3MkpuGMHgcnk7yEKIYQQQqhGpyiXv0qxzWZ7COhmt9tftNlsnwOPXmyq\nmNvtUQJp1QJx+U5X1vCXRdvJ21/s76EEPGutg76Fe+l/Yi8DTuyhf+EeOpYXnf8OFgskJ8NVV8HQ\nob5/+/UD3XlXFBRCCCGE0KLzvnlpauGyFVDq/pcE7AbG2+32E+e7j5b2cYmOjqCoqMLfwwgqx05V\nMWdVFkVlTlJt0Tw9/SrKT2tnTXBNq6nBuDMXU2Y6pswMTBlpGHbb0Z31u+lt3wFXcgqu5FRcSSkU\n9rSx6r1vidqVzbDqwySWHcSUvxOdx9Nwn6g2uJOScSWn4k5OxZ2cgrdTZ3/8hAFLXivUJfmqTzJW\nn2SsPslYXVrL90L7uDRpqpjdbr/+zP8/64rLeYsWEdxy9hfzxvpcHDUexo3oxZ3X9cZqNlIpn/b/\nkNeLYe8ejOk7MGWmY8xIw5iXi662tuGQsHBcw6/BnZyKKzkFd1IK3u49Gl09aQ/84v9ieXH+d/zx\nQAldO4Tx1Jh+dDqyF1NmGsb0NIyZ6YR8sYWQL7bU38/TqXN9EeNKTsWdlIwS1aYlExBCCCGEaJIr\nWlVMtG6KovDpjqOs+GwPBr2eR8YNYtjgTv4elnYoCvpjRzFmpGHKqCtSsjLRVzZ8qqGYTLgHx9UV\nKam4k1Lw9B8AhotPqwyzmnjq7gTe3byXT9OO8sf3dvLExHgGXHV1/TG6slKMdVdxjHVjMH+4CfOH\nm+qPcffp27iYiUsAq7V5sxBCCCGEuEJNmirWFDJVLLi4PV6WfrKbLzILiAwL4clJ8fTtElX/9daY\nsa6kGGNmOqa6qx2m9DT0pxr6UhSdDk//AbiT6gqE5BTcg+PBbG7S9zs7488zjrH0k90ATBtl47qE\nLue9n/54AcaM9IZiJjMdffnphnEajbgHDsadlII7xTc1zWOLBWPr+5yjNZ7HLUnyVZ9krD7JWH2S\n8ZV777O9bM8/+aNfMxh0eDyX/zZ9aGwM99zc77xfX736PXJysnj++T/x4ot/YNCgOCZOvPuij9vs\nU8VE61bpcDF3bQ75h8voERPO7MkJtIu0+HtYLauqClNOVt2UrDRM6ekYDh9sdIinW3dqxt6Jq64A\ncCcmoUREqjKcG5O70rGtlbnrcln4QT7HT1Uz+ca+6PU//N33du5Cbecu1N4xtu4GL4YD++qnl5ky\n0jHmZGHKyYIlCwFQQkNxxyfW/yyupBS8vXpL878QQgghftSkSfewY8f3/OlPz+NyuS6paLkYueIi\nLsvx4ipeXpnNyTIHKQOimTV2EOaQH05rCqqMXS6Mu/IavbE32Heh83rrD/G2beub7lX/xj4VJSZG\n1WH9WMaFJdW8vCqbEyXVJPZtzyPjB2M1N+HzCZcLY/7Ohp85Pe3Hf+b6q0e+n13p2PFKfyxNCarz\nWIMkX/VJxuqTjNUnGatLzXxzc3N49NEZzJ//DjZb7KWOp3lXFWsKKVwCX+6BYl5fl4ejxs2Y4T2Z\ncH0f9Of5xD1gM/Z6Mezf5+tHOdObkpuNrqam/hAlNBRXQlKjqVTenr1a/OrD+TKudrp4fV0ueQdL\n6RodxlOTEujQphl6Vs5cZcpIr7vKlIbh0MFGh3i6dms8FS4xCSUy6scfLwAE7HkcICRf9UnG6pOM\n1ScZq0utfF0uF48/PosxY8azadN65s59C9Ml7Dknhcs55Bfg8iiKwmfpx1j+6R70eh0zRscyPO7C\nTfiBkrH+eAHG9DTfCl/paRizMn7Y7zEozlek1DWvewbYNNHvcaGMPV4vKz7dy+b0o0SEmnhiYjz9\nuzX/6mH1fT11jf+mjHT0RQ1zaBWdDk+//r5iJsW3+IB7cLxv35kAECjncaCSfNUnGatPMlafZKwu\ntfKdM+efREZG8dBDDzN//n+orq7iySd/eSnjkcLlbPILcOncHi/LPt3D5xnHiAw18cSkBPp1vfgn\n6FrMWFda4lth68wyxBnpGAobr+Lt7tuv8Qpbg+M1u8LWpWT8WfpRln2yB70epo+K5Zp4lfdxURT0\nBccaisGMNIyZGT9cSW1QXEPGyamXvJJaS9PieRxMJF/1Scbqk4zVJxmrS2v5SnO+aJJKh2/K0a5D\npXSLDmf25Hg6RGnzTfwPVFdjzMn27WlyZvWsA/sbHeLp3IWaO8bV75USjHua3JzSjY7tQnl9bS7z\n399FQXEVk27oe94pfldMp8PbtRu1XbtRO+5O321n9q7JSPOtZJaZjjE3B1NWBta35/sOCQvHnZDY\nqGg8d+8aIYQQQrRuUriIH3W8uIo5q7IpLHWQ3L8Ds8YNwhKi0dPF7caQv6v+TbEpPQ3Dj+wiX3vD\nTXXTlVrXLvKDe7Xjd9NSmbMqmw+/O8yJ4uqWfT71ejwDbHgG2Ki5d6rvttpajDtzG5ZlzkzH9N03\nhHz7df3dvO3b1+9tU7/gQYcOLTNmIYQQQmiOTBUTP5B3sITX1+ZSXePmjmE9mXjD+Zvwz0e1jBUF\n/YH9DW9265bu1TkcDYdYLL6le8/0VCSn4OndN+g+vb/cjM+9gvbU5ATaR2mn10RXWYExO6t+o0xT\nZjqGw4caHePp3qNRMeNOSEQJj1BtTPJaoS7JV32SsfokY/VJxurSWr4yVUxcsrN7In4yZqD6PREX\noS88UfdGdoevSMlMR19WVv91xWDAEzvIN92rbkleT+xAuIRVK1qbcKuJX9yTyPJP97Al4xgvLNp+\nyT1LLUEJj8A14lpcI66tv01XVIQpK/2spajTsGxYCxvW+u6j0+Gxxfqa/88UM4PiICTEXz+GEEII\nIVQihYsAfKtQLf90D5+lH1N1FaoL0ZWfxpiZ0bAMcUYahuMFjY5x9+5D7c231hUpqbjjEyA0tEXH\nGciMBj0PjrTRpUMYyz7dzd+WpTNj9MCLrhLnL0p0NLW3jqT21pF1NyjojxyunxJozEzHlJmBMX8X\nlhVLfYeEhOCOi2+0x4ynX3/Q6/34kwghhBDiSknhIqhyunijbt+PbtFhzG6ufT8uxOnEmJvd6A2o\nce+eRod4YjpSM+qOhjegSckobdupO65W4pbUbnRsZ+X1dXnM27STguKqC+7Loxk6Hd4ePant0ZPa\n8RN8t3k8GPbsbmj+z0j3LcyQnoaVeQB4IyJxJyY12mPG27Vb0E0fFEIIIYKZFC6tXGFJNS+tyqaw\npJqkfr4m/CbttH4hHg8Ge37DXimZ6Rh35qJzu+sP8UZEUnvdjWctkZuCt3MXeWOporje7Xl2Wiov\nr8zm/W8Pcby4mofHDtTuIgznYzDgiR2IJ3YgNfc94LvN6axr/m+4emf6+ktCvtpafzdvdEzDFMO6\nleWUdu399EMIIYQQ4mKkOb8V23WwhLnrcqlyuhl1dQ8m39AXvf7KCwVdYSEd8tKo/uJr3xvG7Cx0\n1VX1X1fMZtxxCfVvGt3JqXj69JWpPJepuc7jSoeLuWtzyD9cRo+YcGZPTqBdpHaa9puLrvw0xqxM\n30pmdXvMGI4dbXSMp2evuvNyCK6UIbQdfTNFJdV+GnHwk9di9UnG6pOM1ScZq0tr+coGlOfQ2hPk\nD59nHGPpJ7sB38aE1yZcYRO+omD6+kusC+YR8uGm+qWIFb0ejy228UpQsYOkeboZNOd57PZ4WfrJ\nbr7ILCAyLIQnJ8XTt4s2mvbVpCssrC9izqxUpy8tbTigVy8qH5iB8/5pKO3lakxzk9di9UnG6pOM\n1ScZX7mw55/FvHHdj37NoNfh8V7+2/SacXdR9fyL5/3688//jttvH82IEddy8OABXnvtJf7+95cv\n+riyqpio5/F6eXfzXj5NO0q41deEP6B705vwdRXlmN9bgfXttzDa8wFwD47HOP1BymITcMUlQHh4\ncw1fqMRo0DOtrml/xeY9/HVpBjPHxDJskDab9puL0rEjtSNHUztydN0NCvpDBzFlpGHa+jnWtasI\nf/EPhP39z9TcORHHzFm4U4b4d9BCCCFEABg/fgJr165ixIhref/9DYwde+cVP6ZccWlFqp0u3lif\nR+6BErp2CGP25ASim9iEb8jfhXXhPMzvrUBfVYliMlEz7k4cMx7BfdXVRMdEtsqMW5Ja53HO/mLe\nWJ+Lo8bD2BG9uOu63tpv2ldJtMlD5atvYFn4Fsb9+wBwJSXjmPkINXdOBKvKi1gEudb6WtySJGP1\nScbqk4zVpVa+iqIwffoUXn75dX7+88eZP38JRuPFr5nIVLFztMZfgMLSauasyuZ4cTUJfdvz0/GD\nL78J3+Ui5MNNWBe+RcjXXwLg6dIV57QZOB54CCUmpv7Q1phxS1Mz42OnqpizKouiMieptmgeHjMI\nc4hBle+lZfUZe72YvtiCdeE8Qj7+CJ3Xi7dtW5xTp+GYPhNvr97+HmpAktcJ9UnG6pOM1ScZq0vN\nfN9552327NlNTExHHn/8qUsdjxQuZ2ttvwD5h0p5bW0OVU43I6/qzt039rusJnx94QksixdiWfI2\nhhPHAai97kYcM2f5ptj8SPXc2jL2B7UzrqiuZe7aXOxHyujZMYInJ8UHZdP+hfxYxvojh7EuXojl\nnbfRFxej6HTU3nIbzpmzqL35Nllk4jLI64T6JGP1Scbqk4zVpWa+JSXFTJw4hkWLVtCzZ69LHc95\n36TKX9gg90XmMf75bibOWg8PjY7l3pv7X1rRoiiYvv2aiFkP0S55EGF//wu6qiqqH/4pJV/v4PTq\nDdSOGfejRYsIDhGhIfxqShLXJXTmUGEFLyzewf6Ccn8Py++83XtQ9bs/UJyZT/lrb+JOGYL504+J\nmno37a5OwvraHHSlJf4ephBCCOF3Ho+HxMTkSy5aLkYKlyDl9Sos/3QPiz6yYzUb+fWUJK5P7HLx\nO1ZWYln4Fm1vHE6bO0djWb8GT/8BVPzt3xRn5VP157/j6T9A/R9AaILRoOeh0bFMubkf5VW1/HVZ\nOtt2Ffp7WNpgNlNz9xTKPtxM6adbcUx9EH3hCcL/37O0T4wl/KmfYczK8PcohRBCCL/4/PPN/OpX\nT/LTnz7ebI8pU8WCULXTzX825JGzv5gudU34MRdpwjfstvua7d9djr6yAsVopGbseJwzZuEaNuKy\nN4IM9oy1oKUzzt53ijfW5+Gs9TD+ml6Mvzb4m/YvN2NdaQmW5Uuxvv0WhoMHAHClpOKYMcvXzG9p\nXVPtLkZeJ9QnGatPMlafZKwureUrPS7n0NoT1JxOljmYsyqbglNVxPfxNeGHWs4zncvtJuSjD3wN\nx19+AYCnU2ec02bgfPAhvB2bvhRuMGesFf7I+FhRJS+vyubUaSdDY2OYOWYgZlPwNu03OWOvF9Pn\nm337Gn3yX3SKgrd9+4Zm/h49m3+wAUheJ9QnGatPMlafZKwureUr+7i0EvbDpby2NpdKh4vbhnTn\n3pt/vAlfd/Ik1nfexrJ4IYaCYwDUXnOdr9l+1BgwmVp66CJAdI0O59npQ5i7Joft+ScpKnPw5KQE\n2kaY/T00bdHrcd18G66bb0N/6CDWRQuwLFtM6Cv/xvrqS9TePgrHjFm4brxZmvmFEEKISyRXXILE\nl1kFLP6vHYAHbh/ADUldGx+gKBi3fY914ZuYN65H53LhDQun5p4pOGbMwhM7sFnHE4wZa40/M3a5\nvSz5r52vco7TJjyE2ZMT6NUp0i9jUVOzZux0Yl63GuvCeZgy0gFw9+6Dc8bDOKfcj9KmbfN8nwAi\nrxPqk4zVJxmrTzJWl9bylali59DaE3QlvF6FlZ/v5b/bjhBmMfL4hHhie571BqiqCsvq97AufAtj\nXg4Ablusb8793feiRKjzZjOYMtYqf2esKAr/3XaElVv2YjLq+cnYQQyNjbn4HQOIWhkbM9KwLnwL\n89pV6GpqUKxWnBPvxjlzFu74xGb/flrl73O4NZCM1ScZq08yVpfW8pWpYkHKUeNrws/eV0zn9qHM\nnpxAx7ahABj27cGy8C0sK5ahLz+NYjBQM+4uHDNn4Rpx7WU32wtxLp1Ox6ire9CpfSj/2ZDH6+ty\nOX5tb8Zd0wudnF8X5E5OpSI5lcrnX8Sy7B2sb8/HunQx1qWLcQ25CsfMWdSMuwvMMgVPCCGEOEOu\nuASoorom/GOnqojr3Y5H7xxMqFFHyCf/xbrgTUK+2AKAJ6YjzgcfwjltBt7Ol7AccjMJhoy1TksZ\nHz3pa9ovLndy1cAYZt4xkJAgaNpvsYw9HkI++wTLgnmEfPapr5m/QwccD9T97nbrrv4Y/EBL53Cw\nkozVJxmrTzJWl9bylSsuQWb3kTJeXZNDpcPFranduC8hirA3XvY12x89AkDt8GtwzpxFzR3jpNle\nqK5bTDjPTR/Cq2tz2LaroWm/TbhcMbgkBgO1t42i9rZR6A/sr2/mD3vpH4TO+Re1t4/2XS294Sa5\nWiqEEKLVkisuAear7OMs+igfxaswu3s113y1DvPGdehqa1FCw3DePQXHjIfxDBrs13EGcsaBQosZ\nu9xeFn2Uzze5J2gbYWb2pAR6dorw97CazK8ZOxy+Zv4F8zDVbWTp7tvP18x/71SUqDb+GVcz0uI5\nHGwkY/VJxuqTjNWltXybvTnfZrOZgAVAL8AMvGi32zdc6D5SuFwZr1dh1Rf72PLVHm7b9zVT920m\nwp4HgLv/ABwzHqbmnvtQIqP8PFKfQMw40Gg1Y0VR+Oj7w6z6fB8mk55ZYweRagvMpn1NZKwoGNN3\nYF0wD/P6NXUfUoTinHQvjpmz8AyO8+/4roAm8g1ykrH6JGP1Scbq0lq+akwVewAottvtD9pstvZA\nBnDBwkU0naPGzer5H9Nr4zIW7dxCmKMCRa+n5o5xvukj190g00eEZuh0OkYP60mndqG8uXEnr63N\nZcL1fRg7vKc07TeFToc7dSgVqUOp/H9/xrJsMdZFC7AuWYh1yUJcVw/3NfOPGQ8hIf4erRBCCKGa\nphYuK4FVZ/23uxnGIs7l8eBcv5GKf73CL3Zv993UIZqqR3+Kc9pMvF27+XmAQpxf8oBonnkghVdW\nZ7N2636OF1cxY3QsJmPgN+37i9KhA47Zv8Tx+FMNC3F8/hmm77/FGx2D48xCHF26XvSxhBBCiEBz\nRT0uNpstAt+Vlnl2u33ZhY51uz2KUd6wXJpTp2D+fGpfnUvI0cMAHI9NIubZX2O4+275VFUElNIK\nJ39euI38Q6XYerTldzOuom2kxd/DCh579sDrr8PChVBWBgYD3HknPP443CTN/EIIIQJO829AabPZ\nugNrgbl2u33BxY6XHpeLazSPvaYGp8nMF7E3wGOPkjL5Vn8P77JoNeNgEkgZu9we3v4wn2/zCmkX\n6Wva79FR+037gZQx1dVY1qzEsmAeptxsANwDbA39byptNnslAirfACUZq08yVp9krC6t5XuhHhd9\nUx7QZrN1BD4Gnr6UokVcgMOBefk7tLn9BtqOuhnLe8s53a4Tb974Ex57YhHm+W8GXNEixLlMRgMP\njx3EpBv6UFJew5/fSSN9d5G/hxVcQkNxPjCdss1fUrrpE5yT7sFwYD8Rz/wP7RJiCf/NLzDs2unv\nUQohhBBN1tRVxV4G7gXyz7p5tN1ud5zvPnLFpTH9wQNY356PZfkS9KWlKHo9jttGsdJ2G6sNPYlp\nF8bsyQl0bh/m13E2lRYyDnaBmnGavYh5m/KodXmZdEMf7him3ab9QM34DF1REdali7AsWoDh2FFA\nW3s8BXq+gUAyVp9krD7JWF1ay7fZl0NuCilcAK+3YXfszZ/U747tvH86BROm8q/vSjlyspKBPdvy\n2F1xhFsDd+NIrf0SBKNAzvjQiQrmrM6mtKKG4YM78dBomyab9gM540bcbkI+/gjrgnmEbN0CgKdj\nJ5xnmvk7dfbLsIImXw2TjNUnGatPMlaX1vJt9qli4vLoSkuwvjaHdlcnETX1bsyffow7ZQjlr71J\nccYush/6OX/YXMiRk5XcmNyVX9yTGNBFixAX07NTBL+fPoQ+XSL5Nu8Ef1uewemqWn8PK3gZjdTe\nMZbTq9ZT8k0a1bMeRVddTdg//o92KYOJeHg6pm++ghb6IEsIIYRoCrnioiJjVgaWBfOwrF2FzulE\nsVhwTrwb58xZuBOSAPgu7wQLPsjH4/Uy9dYB3JzSVbPTZi6H1qr3YBQMGde6PCz8MJ/vdxbSPtLM\n7MmJdI8J9/ew6gVDxudVWYll9XtYF8zDuKtuM9vYgThmzKLm7ntRwtVfPCGo89UIyVh9krH6JGN1\naS1fmSp2DlWfIKcT8/o1WBfOw5SeBoCnV28cM2bhnDIVpW07ALyKwtqt+3n/20NYzUYeu2swcb3b\nqzMmP9DaL0EwCpaMFUVh07eHWLt1P2aTgUfGDyK5f7S/hwUET8YXpCgYv/8O68I3MW9cj87txhse\nQc09U3DMfATPAJtq37pV5OtnkrH6JGP1Scbq0lq+FypcmroBpTiH/vAhrIsWYFm2GH1xMYpOR83t\no3w72994C+gbZuXV1Hp4a9NO0nYXEdPGylN3B24TvhBXSqfTMW5ELzq3C+WtTTt5dXUOk2/sy6ir\newTF1UfN0+lwDxtOxbDhVP6xEOs7b2NZvBDrgnlYF8yj9trrccyYRe3oMWCUPxlCCCH8R/4KXQmv\nF9Pnn2FdOI+Qjz/yNdu3a0f1Ez/HMX0m3p69fnCXknInc1Znc7iwktgebfjZhHjpZxECGBIbQ3Qb\nK3NWZ7Py830UnKpi2qhYTEZpxWspSseOVP/qaaqf+hUhH77ve237aishX23F07kLzmkzcDzwEErH\njv4eqhBCiFZIpoo1ga6sFMuKpVgWvoXxwH4AXCmpvrnhd04Ey4/vCr6/oJxXVmdzuqqW6xO78MDt\nAzAagvNNmdYuOwajYM24rLKGV1Znc+B4Bf26RfHExHgiQ0P8MpZgzfhyGHbbsS6ch/nd5egrK1BM\nJmrGjscx4xHcVw+DK7gqJvmqTzJWn2SsPslYXVrLV3pcztHUJ8iYk+Vrtl+zEp3DgWI2UzNhMo6Z\ns3AnpVzwvt/vLGTBB7twe7xMubk/tw7pFtTTYLT2SxCMgjnjWpeHBR/sYtuuk3SIsjB7cgLdolu+\naT+YM75cusoKzCvfxbpwHsb8XQC4B8XhmPEwzkn3QPjlPz+Sr/okY/VJxuqTjNWltXylcDnHZT1B\nNTWYN67DumAeph3bAPD06IXjoZ/gnPoASrsLN9R7FYUNXx1gw9cHsZoNPHpnHPF9gqcJ/3y09ksQ\njII9Y0VR2Pj1QdZ9dQBziIGfjh9MUr8OLTqGYM+4SRQF07dfY1kwD/MHG33N/BGROKdMxTljFp5+\n/S/5oSRf9UnG6pOM1ScZq0tr+UpzfhPojx7xNai+8zb6U6d8zfa33o5zxsPU3nwbGC6+WV6Ny8P8\nTTvZYS8iuo2F2ZMT6dpBmvCFuBQ6nY7x1/amc4cw5m/aySursrn7pn6MvKp7UF+t1DydDteIa3GN\nuJaqE8exLPE184fOe4PQeW9Qe/1NOGbOovb2UdLML4QQolnJX5WzKQqmL7b4dpf++EN0Xi/etm2p\n/tlsX7N97z6X/FClFTXMWZ3NoRMVDOjehscnxBHhp3n6QgSyobExdIiy8MrqbN7bsreuad8WtP1h\ngcTbqTPV//MM1T//NSEfbvK9dm7dQsjWLXi6dsM5fSaO+6ejRGtjeWshhBCBTaaKAbrTZVjeXeZr\ntt+3FwBXYjKOmbOouWsSWK2X9fgHjpczZ3U2pytruS6hMw+ObH1vsrR22TEYtbaMG30Y0C2KxyfG\nq/5hQGvLuDkYdu30NfOvfBd9VaWvmX/cXThmPoJ76FWNmvklX/VJxuqTjNUnGatLa/lKj8s5zjxB\nhrxcrAvmYVn9LrrqapSQEGrunOhrtk8Z0qTVcrbtKmT++74m/Htv6sdtQ1vntBat/RIEo9aYcY3L\nw/z3d7Ej39e0/9TkBLqq2LTfGjNuLrqKcszvLce68C2Mu+0AuOIScM6chXPi3RAaKvm2AMlYfZKx\n+iRjdWktXylczuZyEb31Y1wvzcH0/bcAeLr3wDH9JzinPojSoWnNv4qisOHrg6z/6gCWukbixBZu\nJNYSrf0SBKPWmvHZC15YQnwLXiT0VWfBi9aacbNSFExfbcW68C1CPtyEzuPBG9UG55T7Cf3VUxS1\n6eTvEQY1OYfVJxmrTzJWl9bylcLlLGEvPk/onH8BUHvTLThmPkLtrbdfUrP9+Zy7dKvanwIHAq39\nEgSj1p7x2UuMq3V1s7Vn3Nz0Bcd8i54seRt90Umg7nV4xixqbxt5Ra/D4sfJOaw+yVh9krG6tJav\nrCp2lpo7xhLaPoqSkePw9Ol3xY9XWlHDq2t8m+X1r5t376/N8oRoTa4e1JHoNlZeWZPNis/2UlBc\nHdSbugYDb5euVP/vs1T/8jeY399A5JIFhGzZTMiWzXVXvmfinDqtyVe+hRBCBLdW9xfenTIEfv/7\nZilaDp39p+/VAAAgAElEQVSo4MXFOzhwvIJr4zvz6ynJUrQI0YL6dInkuWlD6NExnK1ZBfxzRSaV\nDpe/hyUuJiSEmgmT4csvKdnyDY5pM9EXnyL8xedpnxRLxOOPYEzbDi00I0AIIURgaHWFS3PZkX+S\nv7yTRllFDffc1I8Zd8RiMkqcQrS0dpEWnrk/lVRbNPYjZby4aAcFp6r8PSxxiTyD46j8x0sUZ9up\n/NNf8fToiWXlCtqOvoU2t9+Iefk74HD4e5hCCCE0QN5pXybfbt4HmLsuF51ex5OTExh1dY9WuXKY\nEFphDjHw2F1xjB3Ri5NlDv60ZAc5+4v9PSxxGZTIKByzHqP06x2UrVxPzeixGHOyiHzqZ7RPtBH2\nh9+hP7Df38MUQgjhR1K4XIZal4c3N+5k7ZcHaB9p4XcPpJLUilcOE0JL9DodE6/vwyPjBuFyK7y0\nMotPdhyhpRYgEc1Ep8N1w02UL1pGyY4cqn7+azAaCX39FdoNSybyvkmEfPIReDz+HqkQQogWJoXL\nJSqrrOGvyzL4fmch/bpG8dz0IXSLad0rhwmhRcMGd+Lp+5OJCA1h+ad7WPJfO26P19/DEk3g7dad\n6t/+nuKMXZTPnYd7yFWYN39C1P330O7qZKyvvoyuRK6sCSFEayGFyyU4dKKCFxbt4MDxckbEdeJ/\n7ksmMkya8IXQqr5dovj99CF0jwnn88wC/v1eljTtBzKzmZrJ91L2/ieUbv4SxwPT0RcVEv7H52if\nGEvEk49izEjz9yiFEEKoTAqXi0izn+QvS31N+Hff2JefjBkoTfhCBIB2kRaeeSCF5P4d2HWolBcX\n7+B4sTTtBzp3fCKV/3qF4qx8Kv/4ZzxdumJ5dxltR95Em5E3Yl6xFJxOfw9TCCGECuQd+HkoisKm\nbw7y2tpcdOh4YmI8o4f1lCZ8IQKIJcTI4xPjGTO8JydLHby4OI28AyX+HpZoBkqbtjgefYLSb9Mp\nW7GGmlF3YMzKJHL2Y7RPiiXsj79Hf+igv4cphBCiGUnh8iNcbg/zNu1kzdb9tI80+z61HRDt72EJ\nIZpAr9Mx6Ya+PDx2IC63h3+/l8XmtKP+HpZoLno9rptvpXzxCkq2ZVE9+5eg0xH66ku0uyqRyAfu\nwfTZJ+CVPichhAh0Uric43RVLX9blsF3eYX07RrJs9OH0qNjhL+HJYS4QiPiOvObqSmEW40s/WQ3\nSz6Wpv1g4+3Rk6pnn/c187/6H9wpqZg//og2UybRblgy1rmvoCuVK25CCBGopHA5y+HCCl5YtJ19\nBeUMH9yR39yXTJQ04QsRNPp1jeLZ6UPoFh3OlvRj/Pu9LKqc0rQfdCwWau65j7IPP6P0ky9wTH0Q\n/YnjhD//O9onDST8549jzM709yiFEEJcJilc6mTsLuIv76RTUl7DpBv68PDYQZiMBn8PSwjRzDpE\nWfntgykk9TvTtJ/GiZJqfw9LqMSdmEzlS6/5mvmf/xPemI5Yly2h7a3X02b0LZhXroCaGn8PUwgh\nxCVo9YWLoih88N0hXl2Tg4LC4xPiGTO8lzThCxHELCFGnpgUz+hhPSgsqebFRTvYeVCmEAUzpW07\nHD97kpLvMzm9fBU1t43EmL6DyMcf8TXzv/g8+iOH/T1MIYQQF9CqCxeX28v893ex6vN9tIkw88z9\nqaTapAlfiNZAr9Nx9439+MmYgdS6Pfzr3Sy2pEvTftDT66m95XbKl66k5PtMqh9/CrxeQuf8i3ZD\nE4icNgXTls3SzC+EEBpkbOodbTabHpgLJAI1wMN2u31vcw1MbeVVtby6Joe9x07Tp0skT06MJyrc\n7O9hCSFa2DXxnYlpa+XVNTks+Xg3BaeqmXJrPwz6Vv25Tqvg7dWbqj+8QNVvfot5/RqsC97E/NEH\nmD/6AHefvjhnPIxzyv0oUW38PVQhhBBc2RWXuwCL3W4fDvwv8M/mGZL6DhSc5oVFO9h77DRXD6pr\nwpeiRYhWq3+3Njw3bQhdo8PYnH6Ul1ZmUy1N+62H1UrNlPsp+/gLSv+7Bee9UzEcO0r4c8/QPjGW\n8F/NxpCb4+9RCiFEq3clhcu1wEcAdrv9O2BIs4xIZfbDpTz96pcUlzuZcH0fHhk3iBCTNOEL0dp1\naGPltw+kkti3PXkHSnhxcRoFpyr9PSzRwtzJqVS88gbFmflUPvdHvB2isS55m3Y3X0ObsbdjXrMS\namv9PUyhAW6Plw+/P0TuvlP+HooQl6+yEtO3X2N9bQ789rcBs0iJTlGUJt3RZrO9Bay22+0f1v33\nYaCP3W53/9jxbrdHMWpgla5F7+9k41f7+cV9KVyT0MXfwxFCaIzHq7Do/Z2s/Xwv4VYTzzw0lIR+\n0vvWank88OGHMHeu71+AmBiYNQt++lPo3t2/4xN+UVFdy/8t2k723lOMHtGLn01K9PeQhDi/2lrI\nyYFt22D7dt+/u3Y19PKZTJCXB/37+3ecDc67QtaVFC7/Ar6z2+3v1f33Ubvd3u18xxcVVTTtGzUz\nRVGIjAqlotzh76EEtejoCIqKKvw9jKAmGavry+wClvzXjqLA/bcP4Makrv4eUtAJtHNYf2A/1rfn\nY1m+BH1ZGYpeT+2oMThmzsJ13Q2gwdUoAy3jQHC8uIqXV2VzstRByoBonnnoKnlPoTI5jy+D14th\n316MGWmYMtIwZqZjzM1Bd9YVFSU0FFdCEu7kVNzJKUSOuoUii3Z6+aKjI877Ytrk5nzga2Ac8J7N\nZhsGBMQEYJ1Oh8VsRE5/IcSFXJfQBVvvDry44HsWf2Sn4FQV994sTfutmbd3H6r+35+oevp3WNat\nxrJgHuYPNmL+YCPufv19zfz3TkWJjPL3UIVK8g6UMHddLo4aN2OG92TC9X3kPYXwH0VBX3AMY0Z6\nQ5GSmYG+orzhEKMR96C4+iLFlZSCZ4ANjGeVANERECCF4ZVccTmzqlgCvks6M+x2e/75jtfKFReQ\nyr0lSMbqk4zVFx0dQd6ek8xZlU3BqSrierfj0TvjCLVcyWc+4oyAP4cVBWPadqwL5mHesBZdbS1K\naBjOyffimDkLz6DB/h5h4GesIZvTjrL80z3o9TpmjI5leFwnQDJuCZKxj660pFGRYkpPQ190stEx\n7n79cSel4EpJxZ2UgjsuASyWCz6u1vK90BWXJhcul0sKl9ZFMlafZKy+Mxk7atz8Z0Me2fuK6dw+\nlKcmJxDTNtTfwwt4wXQO64qKsCxfgvXt+RiOHgGgdtgInDNnUXPHOAgJ8cu4giljf3F7vCzfvIct\n6ceIDDXxxKQE+nVtuKomGauvVWZcXY0xOwtTZlrdtK90DAcPNDrE06Vr4yIlMalJy7drLV8pXM6h\ntScoGEnG6pOM1Xd2xl6vwntb9vLx9iOEWYw8MTEeW4+2fh5hYAvKc9jjIeST/2Jd8CYhn3/muymm\nI84HH8I5bQbezi27KExQZtyCqpwu5q7NZdehUrpFhzN7cjwdoqyNjpGM1Rf0GbtcGPN3YsxIbyhS\n7LvQeTz1h3jbtPEVKckpuJOH4E5OwduxU7N8e63lK4XLObT2BAUjyVh9krH6fizjrVm+pn2AB0fa\nuD5RVidsqmA/hw379mB5ez6W5UvRl59GMRiovWMcjhkP47rmuhZp5g/2jNV0oqSal1dlU1hSTXL/\nDswaNwhLyA+niUrG6guqjL1eDAf2+YqUuulextxsdE5n/SGK1Yo7PrGuSEnFlZSCt3cf1V4ztJav\nWs35QgjR6lyf2IWOba28uiaHtz/Mp+BUFffc1A+9XnsrSgn/8vTtT9UL/0fV/z6HZc1KXy/MxnWY\nN67DbYvF8dDD1NwzBSUi0t9DFefYebCEuWtzqa5xc8ewnky8oQ96Da4aJ7RPf+I4xvSGnhRjVgb6\n02X1X1cMBjyxg+qne7mSU/HEDmzcPC/qyRUXoQrJWH2SsfoulPHJUt+nsceLq0no256fjh+M1Sx/\naC5HqzuHFQXj9m1YF7yJeeM6dC4X3rBwau6ZgmPGLN+blWbW6jJuBlvSj7L0kz3o9TB9VCzXxHe+\n4PGSsfoCJWNdWSnGzAxMmen1074MJ443Osbdu0/DCl/JQ3DHxUOof3smtZavTBU7h9aeoGAkGatP\nMlbfxTKudrp5Y30uuQdK6NohjCcnJxDTxnre40Vjrfkc1p08iXXpIiyLFmAoOAZA7TXX4Zg5i9pR\nY3wbwjWD1pzx5fJ4vaz4dC+b048SEWriiYnx9O928UZnyVh9mszY4cCYm+1b4auuSDHu39foEE/H\nTo2WIXYnJaO0beenAZ+f1vKVwuUcWnuCgpFkrD7JWH2XkrHH6+Xdz/by6Y6jhFt9b3YGdNfORl5a\nJucw4HYT8t8PsS58i5CtWwDwdOrc0Mx/hc23kvGlqXa6eH19HnkHSugWHcbsSQl0uMQPISRj9fk9\nY7cbgz2/oUjJTMe4Kw+d211/iDcyCndiMu4UX0+KOyW1xRfjaCq/53sOKVzOobUnKBhJxuqTjNV3\nORl/nnGMpZ/sBmDaKBvXJQTGHyx/knO4McOe3VjefgvLimXoK8pRjEZqxozHOXMWrmEjmtSYKxlf\nXGFpNS+vzOZESTWJfdvzyGVO+5SM1deiGSsK+oMHfNO90tN8/+ZkoauubjjEbMYdl9CwDHFyKp4+\nfSFANyjW2jkszflCCKGyG5O70rFdKHPX5rDwg3yOn6pm8o19pWlfXDJP/wFU/elvVD3zeyyr38O6\nYB6W9WuwrF+De+BgHDMexjn5XggP9/dQg8auQ6XMXZtDldPNqKt7MPkG+Z1tbXSFhXU9KTsw1V1N\n0ZeW1n9d0evx2AbWr/DlTk7BPXBws03nFJdHrrgIVUjG6pOM1deUjAvrllBt6qe3rYmcwxehKJi+\n/xbLgjcxb9qAzu3GGxGJ8977cD70MJ4Btos+hGR8fmdfJZ0+KpZrEy7chH8+krH6mitjXflpjFmZ\njXafNxw72ugYT89ejfZKccUnQljYFX9vLdPaOSxTxc6htScoGEnG6pOM1dfUjM+eL981OoynLmO+\nfGsi5/Cl0xeewLLkbSyLF9avUlR73Q04ZsyidtQd5106VTL+oebuS5OM1dekjJ1OjHk5DcsQZ6Zj\n2LsH3Vnve70dohuWIU5JxZ2YgtK+fTOPXvu0dg5L4XIOrT1BwUgyVp9krL4ryfjcFYoenyBN++eS\nc7gJXC5CPvoA68J5hHy1FQBPl644p83A8cBDKDExjQ6XjBs7dyXA2ZMTiL7CDxUkY/VdNGOPB8Nu\nu69IOdNAvzMXnctVf4g3PAJ3UnL9Xinu5BS8Xbu1yEawWqe1c1gKl3No7QkKRpKx+iRj9TVHxpe7\nJ0RrIufwlTHY87EunIf5vRXoKytQTCZqxt2JY8YjuK+6GnQ6yfgsau29JBmrr1HGioL+8KFGe6WY\nsjLRVVfVH6+EhOCOi/cVKUkpuFOG4OnXP2Cb59WmtXNYCpdzaO0JCkaSsfokY/U1V8Z5B0t4vW4X\n7tFX92CSNAADcg43F11lBeb3VmBdOA+jPR8A9+B4HDMeJuKnMylyaObPr9/YD5fy6hpfE/7tQ7tz\nz039mu13UM5jdemKiuhwYBdVn3/lK1Iy09EXF9d/XdHp8Ayw4U5uWIbYPXAwmM1+HHVg0do5LIXL\nObT2BAUjyVh9krH6mjPjEyXVvLwyi8JSB0n9OvDI+EFYQlp3076cw81MUTB9+zWWBfMwv78BnccD\nRiOugYN9S7bW7S/hscWetycmGG3NKmDJf+0APDjSxvWJzbtUuZzHzUdXWYExOwtjXU+KKSMNw5HD\njY7xdO/hm+p1pkhJSEQJj/DTiIOD1s5hKVzOobUnKBhJxuqTjNXX3BlXOV3MXZvLrkOldIsOZ/bk\neDpEtd6mfTmH1aM/cRzLkrcJ++pzlPR0dDU19V9TQkNxxyfWfzrtSkrB26t30M3193oV3tuyl4+3\nHyHcauLxCXHYerRt9u8j53ET1dZi3JnbsFdKRhqG3fbGzfPt2+NKSsF87QhODxiMKykVJTraj4MO\nTlo7h6VwOYfWnqBgJBmrTzJWnxoZuz1elm/ew5b0Y0SGmnhiUgL9ukY16/cIFHIOqy86OoKighKM\n+TvP+hQ7HUP+TnReb/1x3rZtz2pa9hUzSseOfhz5lXHUuPnPhjyy9xXTpa4JP0allf3kPL4EXi+G\nvXt8U73qliE25uagq62tP0QJDcOVmFS/V4orORVv9x7Sq9UCtJavbEAphBAaYTToefB2G13ah7H8\n0z38bVk6M0YPZHhcJ38PTQQrkwl3fCLu+ESYPtN3W1UVxpzsujeRaZjS0wjZspmQLZvr7+bp2q3R\nCkzupGSUiEg//RCX7mSZgzmrsik4VUVcn3Y8Oj6OUIu83WkxioL+2NFGe6UYMzPQVza8MVZMJtyD\n4hqWIU5K8e1LZDD4ceAiEMhvshBC+MEtqd3o1C6UuetymbdpJwXFVUy4vg/6IJuuIzQqLAz3sOG4\nhw2vv0lXUlx/RebM3hfm9zdgfn8DUNcE3a+/74pMcoqvx2BwPFgs/vopfmD3kTJeXZNDpcPFbUO6\nc8/NfTHISlKqanTeZKRhykhHX3Sy0THu/gOoTRpTX6Ro7bwRgUMKFyGE8JPBvdvx7LRUXl6Vzfvf\nHqLgVBWzxknTvvAPpV17XDffhuvm2+puUNAXHGvUg2DMzMCyZzmW95b7DjnzyXlywzQzT/8Bfvnk\n/MusAhbXNeFPG2XjxqSuLT6GoFdVhSknq24Z4h2+aYeHDjY6xNO1GzVjxjdcqUtMQolsndNhRfOT\nHhehCslYfZKx+loq40qHi7lrc8g/XEaPmHBmT06gXWTwfxop57D6mj1jrxfDvr0Y03c0FDPn9Cp4\nw8JxJyT+aK+CGrxehVWf7+OjbYcJsxj52YR4BvZs/ib88wna89jlatwblZ6Gwb6rcW9UmzYNV+CS\nh6jWGxW0GWuE1vKVHhchhNCwcKuJX96bxNJPdvNFZgF/XLSDJyfG07eVNu0LDdPr8fQfgKf/AGru\nneq77czqUGf1NJi++4aQb7+uv5u3fftGS9i6klJROnS44uE4aty8uSGPrH3FdG4fyuzJCXRsG3rF\nj9vqeL0YDuxrVKQY83LQOZ31hyhWK+6hVwf9anRC26RwEUIIDTAa9EwbaaNLhzBWbN7DX5dlMPOO\nWIYNlqZ9oXEhIb6CJCkF54yHgbP24zjT95CZjvnTjzF/+nH93a50P45TZQ5eXp3NsaIqBvdux2N3\nDibUYmr2Hy8Y6Y8XnDUF0NfTpC8/Xf91xWCob54/c9Wste3/I7RJzkAhhNAInU7HbUO606ldKG+s\nz+XNjTspKK7mrut6S9O+CChKeASuEdfiGnFt/W26U6cwZaY1KmYsG9bChrW+++h0eGyxvpWmzhQz\ng+IgJOQHj7/7SBmvrc2hotrFLandmHJLP2nCPw9dWSnGzAzf1bC67A2FJxod4+7Tl9rbRjb0KsUl\ngLX17jEltEsKFyGE0Jj4Pu353YNDeHlVFpu+Ocjx4ioeHjMIc4gsFSoCl9KhA7W3jqT21pF1Nyjo\njxxumJqUmY4pMwNj/i4sK5b6DgkJwR0X32iPma2OcBZ9vBuvFx4caeOmZGnCr+dw+Ja5zkyrn/Zl\n3L+v0SGeTp2pGT22oUhJTEJp03I9QUJcCSlchBBCg7p0COO56UN5bU0OafYiisrSmD2pdTTti1ZC\np8Pboye1PXpSO36C7zaPB8Oe3XXFzA7fG++cbEzpaViZB8BtIaH069yPyBtGEHW0HHe0F2/Xbq2v\n18LtxpC/q2GRhIx0jLvy0Hk89Yd4o9pQe/1NDcsQJ6fg7dzFj4MW4srIqmJCFZKx+iRj9WkhY7fH\nyzsf29madZyosBCenJRAny7a3wTwUmgh32AXFBk7nXgys9jxzgdYczMZWLSPzqeOoDvr/Ys3OqZu\nZauGPWaUdu1bZHgtkrGioD+wv74nxZSRhjEnC53D0XCIxYI7LqFRkeLp3ReCYApdUJzHGqa1fGVV\nMSGECFBGg57po2Lp0iGcdz/bw1+XpTPzjoFcPaj5lxwVQotO1SjMydNxtNsNDLr2LvreFUexy4Ex\nK7PRHjPmjz/C/PFH9ffz9OxVv0yvOzkFV3wihIX58Se5dPrCE3VN82m+aXRZGehLS+u/ruj1eGIH\n1Rdr7uQU3LGDwCSLE4jgJoWLEEJonE6n4/ah3enUzsob6/P4z4Y8jhdXMf5aadoXwW3vsdO8ujqb\n8moXN6V05b5b+mM06FEsJlzXXo/r2us5c81BV1hYX8ScWZbZsm4NrFsD1L3Ztw1suCKRkqqJN/u6\n8tMYMzMa7T5vKDjW6BhPr944b7zZd0UpKRV3fELAFGFCNCcpXIQQIkAk9O3A7x5M5eVV2Wz4+iAF\nxdX8ZMxAzCZp2hfB55vc47z9YT5eL9x/2wBuSe12weOVjh2pHTma2pGj625Q0B86WL+a1pnpVcZd\nebB0se8QiwX34PhGxYyq06ucTox5OXXFVV1vyt49jQ7xRsdQM3J0w4IEScktNu1NCK1rUuFis9mi\ngHeASCAE+KXdbv+2OQcmhBDih7pGh/Pc9CG8tiaHHfknKSpzMHtSAm0jzP4emhDNwqsorN26n/e/\nPUSo2chjd8UxuHe7y38gnQ5vr97U9OpNzYTJvtvcbgz2/Ib9SzJ807BMadsbvn9kVMP+JXXFTJMa\n2j0eDLvtjVZNM+7MRedyNXyv8Ahqr7vhrFXTUvB26dr6FhoQ4hI19YrLL4HNdrv9JZvNZgOWAynN\nNywhhBDnExEawq/vS2bxR3a+yjnOHxdtZ/akBHp3Do6mfdF6OWvdzNu4k4w9p+jY1srsyQl0bt+M\nU6KMRjyD4/AMjoP7p/luczgw5mb7ipm6AiNk6xZCtm6pv5unY6f6XpL6qyBnLyGsKOgPH2r0GKas\nTHTVVQ2HhITgTkg8a5+aIXj69guK5nkhWkpTC5d/AzVnPYazeYYjhBDiUhgNembcEUuXDmGs3LKX\nvy5NZ+aYgVw1UJr2RWAqPu1kzupsjpysZGDPtjx2Vxzh1hboP7FacQ+9GvfQq+tv0p0ua+g7qStE\nzB+9j/mj9+uPcffugzspGZzVtN+2DX1xcf3Xzmym6Sty6vppBg7+0c00hRCX7qLLIdtstp8Avzjn\n5hl2u327zWbrBHwI/Nxut39xocdxuz2K0SjzsIUQorlt23mCf7yzA0eNh6kjY5ly2wB0MtVEBJD8\nQyX8aeE2yipqGDW8Fz+dEI/RoLErEQUFsH07bNvm+3f7digr832tVy+46ioYOtT3b0oKhIf7dbhC\nBLDz/gFr8j4uNpstHlgB/Nput394seNlH5fWRTJWn2SsvkDK+OjJSuaszubUaSdXDYxh5h0DCdF4\n034g5RuoAiHj7/JOsOCDfDxeL/fd0p9bUrsFRuFd1/zfvldnipCNYdUUCOdxINNavhfax6VJH2fY\nbLZBwEpg6qUULUIIIdTVLSacZ6cNoV+3KLbtOslfl6VTVllz8TsK4SdeRWHN1n28uXEnJqOOX9yd\nyK1DugdG0QL1zf9ER/t7JEK0Gk29DvsXwAK8bLPZPrfZbOubcUxCCCGaIDIshP+Zksw1cZ04cLyC\nFxbt4NAJ7XyKJsQZNbUeXl+by6ZvDhHTxsrvHhxCXB9Z8lcIcWFNas632+13NvdAhBBCXDmTUc/M\nMQPpEh3Gqi37+Ms7aTw8dhBDYmP8PTQhACgp9zXhHy6sJLZHG342Ib5lmvCFEAFPY51vQgghrpRO\np2P01T15YlI8Op2Ouety2fj1AZra0yhEc9lfUM4Li3ZwuLCS6xO78Mt7k6RoEUJcMilchBAiSCX3\nj+a3D6bSPtLM2i8PMG/jTmpdHn8PS7RS3+8s5K/L0imvrmXKLf2ZPsqmvZXDhBCaJq8YQggRxLrH\nhPPs9KH07RrJdzsL+dvyDE5L075oQV5FYd2X+/nPhjwMeh1PTU7k9qEB1IQvhNAMKVyEECLIRYWF\n8Jv7khk+uJNvqs7iHRwulKZ9ob4al4c31uWy4euDdIiy8LsHU0noK034QoimkcJFCCFaAZPRwMNj\nBzLphj6UlNfw53fSSLMX+XtYIoiVVtTwf0vT2WEvYkC3KJ6bPoSu0bIpoxCi6aRwEUKIVkKn0zFm\neC8enxAPwGtrc3j/24PStC+a3YHj5fxx0XYOnajg2oTO/Pq+ZCJCQ/w9LCFEgGvScshCCCECV6ot\nmug2qcxZnc3qL/ZTcKqKh0bHYjIa/D00EQS27Spk/vu7cLu93HtzP+lnEUI0G7niIoQQrVCPjhE8\nN20IfbpE8m1eXdN+Va2/hyUCmKIorP/qAG+sz0Ov1zF7cgIjr+ohRYsQotlI4SKEEK1UVLiZp6cm\nM2xQR/YdK+fFRdulaV80Sa3Lw3825LH+qwP1TfiJ/Tr4e1hCiCAjhYsQQrRiJqOBWeMGMeH6PhSX\n1/CXd9LJ2CNN++LSlVbU8Ndl6WzbdZL+3aJ4dvoQukkTvhBCBVK4CCFEK6fT6Rg3ohc/uysOBYVX\nV+fw4XeHpGlfXNShExW8uHgHB45XcE18J349JZlIacIXQqhEmvOFEEIAMCQ2hug2Vuaszmbl5/s4\ndqqK6aNiMRnlMy7xQzvyT/LWpp243F7uvqkvo6SfRQihMvlrJIQQol7PThE8N30IvTtH8E3uCf6+\nIoNyadoXZ1EUhY1fH2Duulx0eh1PTIpn9NU9pWgRQqhOChchhBCNtAk38/TUFK4aGMPeo6d5YdEO\njp6s9PewhAbUujy8uXEna788QPtIM799IJXk/tH+HpYQopWQwkUIIcQPhJgM/HT8YO66rjfF5U7+\n9E4amXtP+XtYwo/KKmv467IMvt9ZSL+uUTw3fSjdY6QJXwjRcqRwEUII8aN0Oh3jr+nNY3fFoXgV\nXlmVzUffH5am/Vbo0IkKXli0gwPHyxk+uBP/c18SkWHShC+EaFnSnC+EEOKChsbG0CHKwiurs3lv\ny2SkuSsAAAisSURBVF4KTlUxbZQNo0E++2oN0uxFzNuUh8vlZfKNfRl9tTThCyH8Q/7qCCGEuKje\nnSN5bvpQenWK4Kuc4/xjeQbl1dK0H8wURWHTNwd5bW0OOnQ8PjGeO4ZJE74Qwn+kcBFCCHFJ2kaY\nefr+FIbExrD76GleXLSDo0XStB+MXG4P8zbtZM3W/bSLNPPMAymkDJAmfCGEf0nhIoQQ4pKZTQYe\nvXMw46/pxanTTv68JI3sfdK0H0xOV9Xyt2UZfJdXSJ8ukTw3bQg9Okb4e1hCCCGFixBCiMuj1+m4\n67o+PHrnYDxehZdXZfPxNmnaDwaHCyt4YdF29hWUM2xwR56emkxUuNnfwxJCCECa84UQQjTRVQM7\nEt3GypzV2az4bC8FxVU8cLs07QeqjN1FvLlxJzUuDxOv78OY4dLPIoTQFvnrIoQQosl6d/ZNJerZ\nMYKtWcf554pMKqRpP6AoisIH3x3i1TU5KCg8PiGOsSN6SdEihNAcKVyEEEJckXaRFv73/hRSbdHY\nj5Tx4uIdHDtV5e9hiUvgcnuZ//4uVn2+jzYRZp65P5VUW4y/hyWEED9KChchhBBXzBxi4LG7fJ/U\nF5U5+fOSHeTsL/b3sMQFlFfV8vcVGXyTe6Juuesh9OwkTfhCCO2SwkUIIUSz0Ot0TLy+D4+MG4TL\nrfDSyiw+2X5EmvY16MjJSl5YtIO9R09z1cAYnp6aTBtpwhdCaJw05wshhGhWwwZ3IrqtlVdW57B8\n8x4Kiqu4/7YB0rSvEZl7TvGfjXnU1HqYcF1v6WcRQgQM+SsihBCi2fXtEsXvpw+hR0w4X2QW8K93\nM6l0uPw9rFZNURQ+/P4Qr6zORvEq/OyuOMZd01uKFiFEwJDCRQghhCraRVp45oFUUgZEk3/Y17R/\nvFia9v3B5fay4INdrNyyj6jwEP73gRSGxEoTvhAisFxR4WKz2WJtNttpm81maa4BCSGECB7mEAM/\nmxDHmOE9OVnq4MXFaeQekKb9llReXcs/VmTwdc4JenWK4LnpQ+nVKdLfwxJCiMvW5MLFZrNFAv8E\nappvOEIIIYKNXqdj0g19mTV2EC63h5fey2Zz2lFp2m8Bh46X8+KiHew5epqhsTE8fX8KbSOkCV8I\nEZia1Jxvs9l0wJvAb4H1zToiIYQQQWl4nK9p/9XV2Sz9ZDdHTlXRLjzE38MKWm6Pl81px3DUuLnz\n2t6Mv0aa8IUQgU13sU+8bDbbT4BfnHPzIWCF3W5fYrPZDgKxdrvdeaHHcbs9itFouIKhCiGECAYn\nS6p5YcH3HDxe7u+hBL0Qo56f35fCdUld/T0UIYS4VOf9hOWihcuPsdlse4Gjdf85DNhmt9uvv9B9\niooqNDMnIDo6gqKiCn8PI6hJxuqTjNUnGavH5fZQXO2mpESa9dU0qF8MuN3+HkZQk9cJ9UnG6tJa\nvtHREectXJo0Vcxut/c78//rrrjc3pTHEUII0TqZjAbi+7ahKFL6LdQU3daqqTckQghxJWQ5ZCGE\nEEIIIYTmNemKy9nsdnuvZhiHEEIIIYQQQpyXXHERQgghhBBCaJ4ULkIIIYQQQgjNk8JFCCGEEEII\noXlNWg5ZCCGEEEIIIVqSXHERQgghhBBCaJ4ULkIIIYQQQgjNk8JFCCGEEEIIoXlSuAghhBBCCCE0\nTwoX8f/buZcQK+s4jOPf0bGGRCnoYpHUyieIsIVQaU2zETEoImhXpEOFVHTZRA0ZFEUEZVAhlWlO\nt0U3o1ykVEZZRCAu3PgMWbsuxHS1snKaFu8ZxpWTcOr/zv88HxiYd/dlOJxzfv/3905EREREROtl\ncImIiIiIiNbrLx3wf5I0B9gILAX+AK63/XnZqvpIugB42PZQ6ZbaSJoHbAHOBo4HHrD9VtGoykia\nC2wCBEwAa20fKFtVJ0mnAnuAlbb3l+6pjaS9wE+dyy9try3ZUxtJdwNXAMcBG21vLpxUFUlrgDWd\nywHgfGCR7R9LNdWm851ilOY7xQRwQ9vfi3vtjsuVwIDti4C7gEcL91RH0p3AszRvMtF91wDjti8B\nVgNPFu6p0eUAtlcA9wIbyubUqfOB+TTwe+mWGkkaALA91PnJ0NJFkoaA5cAK4FJgcdGgCtneOvX6\npTnguDVDS9ddBvTbXg7cDzxYuGdGvTa4XAy8A2D7U2BZ2ZwqHQCuKh1RsVeB9UdcHy4VUivbbwI3\ndi7PAr4tmFOzR4CngK9Kh1RqKXCCpJ2S3pd0YemgyqwC9gHbgLeB7WVz6iVpGXCu7WdKt1RoDOjv\nbCQtBP4q3DOjXhtcFjJ92xxgQlJPrcv912y/zix44c9Wtg/a/kXSAuA14J7STTWyfVjSKPAEzd85\nuqizAvKd7R2lWyr2G81wuApYB7yUz7uuOpnm8PNqpv++fWWTqjUC3Fc6olIHadbE9tOsSD9etOZf\n6LXB5WdgwRHXc2znxDpmFUmLgV3AC7ZfLt1TK9vXAUuATZLml+6pzDCwUtIHNHvrz0taVDapOmPA\ni7YnbY8B48DphZtqMg7ssP2nbQOHgFMKN1VH0onAObZ3lW6p1B00r+MlNHdpR6fWTNuq105fPqbZ\nX3+lc9t8X+GeiGMi6TRgJ3CL7fdK99RI0rXAmbYfojm1/pvmocXoEtuDU793hpd1tr8pV1SlYeA8\n4CZJZ9BsHHxdNqkqu4HbJG2gGQjn0wwz0V2DwLulIyr2A9NbMt8D84C55XJm1muDyzaaU75PgD4g\nDyvGbDMCnASslzT1rMtq23nAuXveAJ6T9CHNm/jttg8Vboo4VpuBrZJ2A5PAcDYMusf2dkmDwGc0\n2ys3284BR/cJ+KJ0RMUeA7ZI+ojmv+ON2P61cNNR9U1OTpZuiIiIiIiIOKpee8YlIiIiIiJmoQwu\nERERERHRehlcIiIiIiKi9TK4RERERERE62VwiYiIiIiI1svgEhERERERrZfBJSIiIiIiWi+DS0RE\nREREtN4/u6cPBCJ5QIcAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x11cd40410>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.plot(x)\n",
    "plt.plot(y, c='r')\n",
    "plt.legend(('x', 'y'))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In this simple example, there is only one value or \"feature\" at each time index. However, in practice, you can use sequences of *vectors*, e.g. spectrograms, chromagrams, or MFCC-grams."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Distance Metric"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "DTW requires the use of a distance metric between corresponding observations of `x` and `y`. One common choice is the **Euclidean distance** ([Wikipedia](https://en.wikipedia.org/wiki/Euclidean_distance); FMP, p. 454):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 128,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "5.0"
      ]
     },
     "execution_count": 128,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "scipy.spatial.distance.euclidean(0, [3, 4])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 129,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "13.0"
      ]
     },
     "execution_count": 129,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "scipy.spatial.distance.euclidean([0, 0], [5, 12])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Another choice might be the **cosine distance** ([Wikipedia](https://en.wikipedia.org/wiki/Cosine_similarity); FMP, p. 376) which can be interpreted as the (normalized) angle between two vectors:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 134,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.0"
      ]
     },
     "execution_count": 134,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "scipy.spatial.distance.cosine([1, 0], [100, 0])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 132,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "1.0"
      ]
     },
     "execution_count": 132,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "scipy.spatial.distance.cosine([1, 0, 0], [0, 0, -1])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 133,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "2.0"
      ]
     },
     "execution_count": 133,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "scipy.spatial.distance.cosine([1, 0], [-1, 0])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Step 1: Cost Table Construction"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As described in the notebooks [Dynamic Programming](dp.html) and [Longest Common Subsequence](lcs.html), we will use dynamic programming to solve this problem. First, we create a table which stores the solutions to all subproblems. Then, we will use this table to solve each larger subproblem until the problem is solved for the full original inputs."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The basic idea of DTW is to find a path of index coordinate pairs the sum of distances along the path $P$ is minimized:\n",
    "\n",
    "$$ \\min \\sum_{(i, j) \\in P} d(x[i], y[j]) $$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The path constraint is that, at $(i, j)$, the valid steps are $(i+1, j)$, $(i, j+1)$, and $(i+1, j+1)$. In other words, the alignment always moves forward in time for at least one of the signals. It never goes forward in time for one signal and backward in time for the other signal."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here is the optimal substructure. Suppose that the best alignment contains index pair `(i, j)`, i.e., `x[i]` and `y[j]` are part of the optimal DTW path. Then, we prepend to the optimal path \n",
    "\n",
    "$$ \\mathrm{argmin} \\ \\{ d(x[i-1], y[j]), d(x[i], y[j-1]), d(x[i-1], j-1]) \\} $$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We create a table where cell `(i, j)` stores the optimum cost of `dtw(x[:i], y[:j])`, i.e. the optimum cost from `(0, 0)` to `(i, j)`. First, we solve for the boundary cases, i.e. when either one of the two sequences is empty. Then we populate the table from the top left to the bottom right."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 135,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def dtw_table(x, y):\n",
    "    nx = len(x)\n",
    "    ny = len(y)\n",
    "    table = [[0 for _ in range(ny+1)] for _ in range(nx+1)]\n",
    "    \n",
    "    # Compute left column separately, i.e. j=0.\n",
    "    for i in range(1, nx+1):\n",
    "        d = scipy.spatial.distance.euclidean(x[i-1], 0)\n",
    "        table[i][0] = d + table[i-1][0]\n",
    "        \n",
    "    # Compute top row separately, i.e. i=0.\n",
    "    for j in range(1, ny+1):\n",
    "        d = scipy.spatial.distance.euclidean(0, y[j-1])\n",
    "        table[0][j] = d + table[0][j-1]\n",
    "        \n",
    "    # Fill in the rest.\n",
    "    for i in range(1, nx+1):\n",
    "        for j in range(1, ny+1):\n",
    "            d = scipy.spatial.distance.euclidean(x[i-1], y[j-1])\n",
    "            table[i][j] = d + min(table[i-1][j], table[i][j-1], table[i-1][j-1])\n",
    "    return table"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 136,
   "metadata": {},
   "outputs": [],
   "source": [
    "table = dtw_table(x, y)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's visualize this table:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 142,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "             1   3   4   3   1  -1  -2  -1   0\n",
      "     +----------------------------------------\n",
      "     |   0   1   4   8  11  12  13  15  16  16\n",
      "   0 |   0   1   4   8  11  12  13  15  16  16\n",
      "   4 |   4   3   2   2   3   6  11  17  20  20\n",
      "   4 |   8   6   3   2   3   6  11  17  22  24\n",
      "   0 |   8   7   6   6   5   4   5   7   8   8\n",
      "  -4 |  12  12  13  14  12   9   7   7  10  12\n",
      "  -4 |  16  17  19  21  19  14  10   9  10  14\n",
      "   0 |  16  17  20  23  22  15  11  11  10  10\n"
     ]
    }
   ],
   "source": [
    "print '         ', ''.join('%4d' % n for n in y)\n",
    "print '     +' + '----' * (ny+1)\n",
    "for i, row in enumerate(table):\n",
    "    if i == 0:\n",
    "        z0 = ''\n",
    "    else:\n",
    "        z0 = x[i-1]\n",
    "    print ('%4s |' % z0) + ''.join('%4d' % z for z in row)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The time complexity of this operation is $O(N_x N_y)$. The space complexity is $O(N_x N_y)$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Step 2: Backtracking"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To assemble the best path, we use **backtracking** (FMP, p. 139). We will start at the end, $(N_x - 1, N_y - 1)$, and backtrack to the beginning, $(0, 0)$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Finally, just read off the sequences of time index pairs starting at the end."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 138,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def dtw(x, y, table):\n",
    "    i = len(x)\n",
    "    j = len(y)\n",
    "    path = [(i, j)]\n",
    "    while i > 0 or j > 0:\n",
    "        minval = numpy.inf\n",
    "        if table[i-1][j] < minval:\n",
    "            minval = table[i-1][j]\n",
    "            step = (i-1, j)\n",
    "        if table[i][j-1] < minval:\n",
    "            minval = table[i][j-1]\n",
    "            step = (i, j-1)\n",
    "        if table[i-1][j-1] < minval:\n",
    "            minval = table[i-1][j-1]\n",
    "            step = (i-1, j-1)\n",
    "        path.insert(0, step)\n",
    "        i, j = step\n",
    "    return path"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 139,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[(0, 0),\n",
       " (1, 0),\n",
       " (1, 1),\n",
       " (2, 2),\n",
       " (2, 3),\n",
       " (3, 3),\n",
       " (3, 4),\n",
       " (4, 5),\n",
       " (4, 6),\n",
       " (5, 7),\n",
       " (6, 7),\n",
       " (7, 8),\n",
       " (7, 9)]"
      ]
     },
     "execution_count": 139,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "path = dtw(x, y, table)\n",
    "path"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The time complexity of this operation is $O(N_x + N_y)$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As a sanity check, compute the total distance of this alignment:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 144,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "10"
      ]
     },
     "execution_count": 144,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "sum(abs(x[i-1] - y[j-1]) for (i, j) in path if i >= 0 and j >= 0)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Indeed, that is the same as the cumulative distance of the optimal path computed earlier:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 146,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "10.0"
      ]
     },
     "execution_count": 146,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "table[-1][-1]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "[&larr; Back to Index](index.html)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 2",
   "language": "python",
   "name": "python2"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.13"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}