{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div>\n",
    "<a href=\"http://www.music-processing.de/\"><img style=\"float:left;\" src=\"../data/FMP_Teaser_Cover.png\" width=40% alt=\"FMP\"></a>\n",
    "<a href=\"https://www.audiolabs-erlangen.de\"><img src=\"../data/Logo_AudioLabs_Long.png\" width=59% style=\"float: right;\" alt=\"AudioLabs\"></a>\n",
    "</div>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div>\n",
    "<a href=\"../C3/C3.html\"><img src=\"../data/C3_nav.png\" width=\"100\"  style=\"float:right;\" alt=\"C3\"></a>\n",
    "<div><h1>DTW Variants</h1> </div>\n",
    "</div>\n",
    "\n",
    "<p>\n",
    "Various modifications have been proposed in order to speed up DTW computations as well as to \n",
    "better control the overall course of the warping paths. Following Section 3.2.2 of <a href=\"http://www.music-processing.de\">[Müller, FMP, Springer 2015]</a>, we discuss in this notebook some of these DTW variants. Efficient and robust multiscale versions of DTW are described in the following articles.   \n",
    "    \n",
    "<ul>    \n",
    "<li><span style=\"color:black\">\n",
    "Thomas Prätzlich, Jonathan Driedger, and Meinard Müller:  <a href=\"https://www.audiolabs-erlangen.de/fau/professor/mueller/publications/2016_PraetzlichDriedgerMueller_MrMsDTW_ICASSP.pdf\"><strong>Memory-Restricted Multiscale Dynamic Time Warping.</strong></a> Proceedings of the IEEE International Conference on Acoustics, Speech, and Signal Processing (ICASSP): 569&ndash;573, 2016. \n",
    "<br>\n",
    "<a type=\"button\" class=\"btn btn-default btn-xs\" target=\"_blank\" href=\"../data/bibtex/FMP_bibtex_PraetzlichDM16_MsDTW_ICASSP.txt\"> Bibtex </a>\n",
    "</span></li>\n",
    "    \n",
    "<li><span style=\"color:black\">\n",
    "Meinard Müller, Henning Mattes, and Frank Kurth: <a href=\"https://www.audiolabs-erlangen.de/fau/professor/mueller/publications/2006_MuellerMattesKurth_MultiscaleAudioSynchronization_ISMIR.pdf\"><strong>An Efficient Multiscale Approach to Audio Synchronization.</strong></a> Proceedings of the International Conference on Music Information Retrieval (ISMIR): 192–197, 2006.\n",
    "<br>\n",
    "<a type=\"button\" class=\"btn btn-default btn-xs\" target=\"_blank\" href=\"../data/bibtex/FMP_bibtex_MuellerMK06_EfficientMultiscaleApproach_ISMIR.txt\"> Bibtex </a>\n",
    "</span></li>\n",
    "</ul>   \n",
    "</p> "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Step Size Condition\n",
    "\n",
    "In [classical DTW](../C3/C3S2_DTWbasic.html), the step size condition is expressed by the set $\\Sigma=\\{(1,0),(0,1),(1,1)\\}$. Introducing a kind of local continuity condition, this condition ensures that a warping path aligns each element of the sequence $X=(x_1,x_2,\\ldots,x_N)$ to an element of $Y=(y_1,y_2,\\ldots,y_M)$ and vice versa. One drawback of this condition is that a single element of one sequence may be assigned to many consecutive elements of the other sequence, which leads to vertical and horizontal sections in the warping path. Intuitively, in such cases the warping path is stuck at some position in one of the sequences, while moving on in the other sequence. In terms of physical time, this situation corresponds to a strong temporal deformation in the alignment of the two time series. To avoid such degenerations, one can modify the step size condition by constraining the slope of the admissible warping paths, which can be done by replacing the set $\\Sigma$. This is illustrated by the following figure, which shows three different step size conditions (top) along with some typical warping paths (bottom):\n",
    "\n",
    "<img src=\"../data/C3/FMP_C3_F17.png\" width=\"600px\" align=\"center\" alt=\"C1\">\n",
    "\n",
    "In the first case, the original set $\\Sigma=\\{(1,0),(0,1),(1,1)\\}$ is used (observe the degenerations in the warping paths). Replacing this set by \n",
    "\n",
    "\\begin{equation}\n",
    "   \\Sigma = \\{(2,1),(1,2),(1,1)\\}\n",
    "\\end{equation}\n",
    "\n",
    "leads to warping paths having a local slope within the bounds $1/2$ and $2$. For computing an optimal warping path that satisfies the new step size constraints, one only has to slightly modify the [original DTW algorithm](../C3/C3S2_DTWbasic.html). To compute the accumulated cost matrix $\\mathbf{D}$, one can use the following recursion:\n",
    "\n",
    "\\begin{equation}\n",
    "  \\mathbf{D}(n,m)= \\mathbf{C}(n,m) + \\min\\left\\{\n",
    "             \\begin{array}{l}\\mathbf{D}(n-1,m-1),\\\\ \\mathbf{D}(n-2,m-1),\\\\\\mathbf{D}(n-1,m-2) \\end{array}\\right.\n",
    "\\end{equation}\n",
    "\n",
    "for $n\\in [1:N]$ and $m\\in [1:N]$ with $(n,m)\\not=(1,1)$. For the initialization, one can use a trick of extending $\\mathbf{D}$ by two additional rows and columns (indexed by $-1$ and $0$) and set $\\mathbf{D}(1,1):=\\mathbf{C}(1,1)$, $\\mathbf{D}(n,-1):=\\mathbf{D}(n,0):=\\infty$ for $n\\in [-1:N]$, and $\\mathbf{D}(-1,m):=\\mathbf{D}(0,m):=\\infty$ for $m\\in [-1:M]$. Note that, with respect to the modified step size condition, there is a warping path of finite total cost between two sequences $X$ and $Y$. Furthermore, note that not all elements of $X$ need to be assigned to some element of $Y$ and vice versa, which is illustrated by the figure above. The third case of the above figure shows a step size condition which avoids such omissions while imposing constraints on the slope of the warping path. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## libfmp Implemetation\n",
    "\n",
    "In the next code cell, we provide an implementation the DTW algorithm with step size condition $\\Sigma=\\{(1,1),(2,1),(1,2)\\}$. Again note (as already discussed in the [implementation of the classical DTW algorithm](../C3/C3S2_DTWbasic.html)) that there is a systematic index shift of minus one in relation to the algorithmic description given above due to Python indexing conventions."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2024-02-15T08:50:49.281415Z",
     "iopub.status.busy": "2024-02-15T08:50:49.281201Z",
     "iopub.status.idle": "2024-02-15T08:50:53.298648Z",
     "shell.execute_reply": "2024-02-15T08:50:53.298023Z"
    }
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAagAAACICAYAAACyaX9CAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAAAndklEQVR4nO3deVyVVf7A8c8BQVRcEi1NXMilnBy10tQsTMSFJTXN0am0aabsV7ZpyziaU1NWP0ctbfpNamZmNZmaWoLghitqLmnlimmaZJOK+wIC9/v74yikIqLcy/MA3/fr9byEy3Of5+sV7/eec77nHCMiKKWUUm7j53QASimlVF40QSmllHIlTVBKKaVcSROUUkopV9IEpZRSypXK+OKi1apVk3r16vni0koppUqY9evXHxSR6hc+7pMEVa9ePdatW+eLSyullCphjDF78npcu/iUUkq5kiYopZRSrqQJSpVIs2fPJj4+3ukwlFKF4JMxKKWcIiK8+uqrvPLKKxhj+OSTT7j//vudDkspn8rMzCQ1NZX09HSnQ8lXUFAQoaGhBAQEFOh8TVCqxMjKyuLxxx9n4sSJ9OvXj59++ol+/fpRtmxZevbs6XR4SvlMamoqFStWpF69ehhjnA4nTyJCWloaqamphIWFFeg52sWnSoQTJ07QrVs3Jk6cyLBhw5g8eTJz5syhVatW9OnThzlz5jgdolI+k56eTkhIiGuTE4AxhpCQkCtq5WmCUsXer7/+yt13301iYiLjx4/n1VdfxRhDcHAwc+fOpXnz5tx3333Mnz/f6VCV8hk3J6dzrjRGTVCqWEtJSaFNmzZs3bqVL7/8kv79+5/388qVKzNv3jwaN25Mt27dWLJkiTOBKqWumCYoVWytWrWKO+64gxMnTrB48WJiY2PzPK9q1aosWLCAG264gdjYWJKTk4s4UqVKtr179xIWFsahQ4cAOHz4MGFhYezZk+f82wLTBKWKpdmzZxMREcE111zDypUruf322/M9v3r16ixatIhatWoRFRXFmjVriihSpUq+2rVr8/jjjzN48GAABg8eTP/+/albt26hrqsJShU7//73v+nZsyfNmjVj5cqVNGjQoEDPq1GjBosWLaJatWp07tyZjRs3+jZQpUqRgQMHsnr1asaMGcOKFSt47rnnCn3NApWZG2MGAo8AAnwPPCwi7i64VyWOx+NhyJAhjBgxgnvuuYepU6dSvnz5K7pGaGgoSUlJhIeHExkZyZIlS2jSpImPIlaq6D377LNe//DVvHlzxowZk+85AQEBjBw5ki5dujB//nwCAwMLfd/LtqCMMbWAp4EWItIE8Af6FPrOSl2BM2fO0K9fP0aMGMFjjz3GzJkzrzg5nVOvXj2SkpIIDAwkMjKS7du3ezlapUqnhIQEatasyaZNm7xyvYJO1C0DlDPGZALlgX1eubtSBXD06FF69OhBUlISb7zxBoMHDy50SW2DBg1ISkqiXbt2REREsGzZMurXr++liN3vn/+E99+HevWgUSN7NGxo/6xXD8roFP5i63ItHV/ZuHEjCxYsYPXq1dx555306dOHmjVrFu6iInLZA3gGOAEcAD69xDn9gXXAujp16ohS3pCamipNmzaVMmXKyJQpU7x+/e+++06qVq0qderUkT179nj9+m40ebIIiLRqJdKypUjlyvb7c0eZMiI33ihyzz0igwaJjBsnkpQkkpoq4vE4Hb3Ky5YtWxy9v8fjkdatW8v8+fNFROSdd96R+++/P89z84oVWCd55JXLfk4yxlwDdAPCgCPAdGPMgyLyyQWJbgIwAaBFixZSuLSpFGzevJkuXbpw9OhR5s6dS8eOHb1+j9///vcsWLCAiIiInJbU9ddf7/X7uMWyZfDooxARAYmJEBBg09KBA7BjB6Sk5B47dsCCBfDbif/ly+e2tH7b6mrUCEJCnPt7KWe9//771KlTJ+f/6BNPPMHkyZNZunQp7dq1u+rrGpu88jnBmF5AFxH5y9nv+wGtReSJSz2nRYsWohsWqsJYsmQJ3bt3p3z58jmrQfjS119/TWRkJKGhoSxZsoTrrrvOp/dzwg8/QKtWUL06rFoF11xz+ed4PJCampuwfpvAfvwRsrNzz61a9fyEdS6BNWwIwcG++3sp2Lp1K40bN3Y6jALJK1ZjzHoRaXHhuQXpaf4JaG2MKQ+cBjpgu/KU8ompU6fy0EMPUb9+fRITE6lTp47P79mqVSvmzp1Lly5diIyMZPHixVSrVs3n9y0qhw9DTIz9Oi6uYMkJwM8P6tSxR2Tk+T/LzLRJ6rctrpQUWLwYPv74/HOvvz7vVtcNN4AXir1UCXXZBCUiXxtjZgDfAFnABs525SnlTSLCW2+9xfPPP89dd93Fl19+yTUFfSf1grvuuos5c+YQExNDp06dSEpKokqVKkV2f1/JzIT77rPJZOFCKOC0scsKCMhNNBc6dcq22C7sMpw5Ew4ezD3Pz+/8Qo3fJrDatcHf3zuxquKpQLU6IvIy8LKPY1GlWHZ2Ns899xxjx46lV69eTJkyhaCgoCKPIyIiglmzZtGtW7ec+RyVKlUq8ji8RQSeeAKSkuCjjyA8vGjuW748NG1qjwsdOpTb2vptt+Hy5XDyZO55ZcvaZHpht2GjRnDttVAM1kZVhaTFpMpxp0+fpm/fvnzxxRcMHDiQUaNG4efn3CInXbp0Ydq0adx3333ExMSQmJhIhQoVHIunMEaPhokTYcgQ6NfP6WisqlXtWFirVuc/LgL//e/Fra7t2yE+3rYEz6lYMe9WV8OGUAIaveosTVDKUYcOHaJr166sXLmSt956i4EDBzodEgDdunXjP//5D3369KFr167ExcVRrlw5p8O6IrNnw4sv2u69115zOprLMwZq1rTHhYVf2dmwZ8/FhRqrV8PUqTa5nXPttXkXazRoAMXsn7DU0wSlHLN7926ioqLYtWsXn3/+Ob169XI6pPP06tWLjIwM+vXrR48ePZg9ezZly5Z1OqwC+eYbeOABaNHCdu052CD1Cn9/W1Bxww3QufP5P0tPh127Lu4yTEiADz/MPc8YO66VV7GGTk52J/0nUY7YsGED0dHRpKens2DBAsKLanDkCj344INkZGTwyCOP8Ic//IEZM2YQEBDgdFj5+vlnuOceOy/pq6/seFBJFhQEv/udPS507Nj5xRrnEtinn8LRo7nnlSljk99vW13Nm1/cDanyJiLcddddDB06lKioKACmTZvGpEmTSExMLNyFvX3cdtttVz0jWZV8iYmJEhwcLHXq1JHNmzc7HU6BvPvuuwJIr169JDMz0+lwLunECZFbbhEJDhb59luno3Evj0dk/36R5GSRDz8U+dvfRHr2FGnaVKRcudxVNR5+2L6mbuf0ShIiIt9//73cdNNNcvr0aTlx4oQ0aNBAfvjhh4vO8+pKEkp50+TJk3n00Ue5+eabmTt3brFZtWHAgAFkZGTw3HPPERgYyEcffYS/y2qgPR548EH49lvbcsqrgk5ZxtgJy9Wrwx13nP8zj8e2QsePhzfesONcn38Ov/+9M7EWF02aNOGee+5hxIgRnDx5kn79+hV6fUtNUKpIiAivv/46w4YNIzIyki+++KLYlW8PGjSI9PR0hg4dSlBQEBMmTHC02vBCgwfbwogxY3In5aor5+dnx6qGD4f27W3Sv/12+7r27+/+8vZnnwVvb3XWvLn9+1/Oyy+/zK233kpgYCDeWE1IE5TyuaysLAYMGMCECRPo27cvEydO9MpeMU4YMmQI6enpvPbaa5QtW5Z333230Cure8PEiTByJDz+ODz9tNPRlBwdOtg3+3794H/+BxYtggkTtJT9UipUqEDv3r0JDg72SkGRJijlUydPnqR3797Ex8czZMgQhg8f7oo39ML4xz/+QXp6OiNHjqRs2bKMHj3a0b9TUpJNTJ06wTvvuP8TfnFz3XW2InDkSBg6FNautV1+t9/udGR5c2i3jRx+fn5e61lwT/+EKnH2799P+/btSUhI4L333uP1118v9skJwBjDiBEjeOqpp3j77bcZOnTouS1nitz27dCzp606mzZNS6V9xc8P/vpXu9qFCLRtC6NG2fEq5Tv666x8YseOHURFRbFv3z5mzZpF165dnQ7Jq4wxjB07loyMDN58803KlSvHsGHDijSGtDQ71hQQYBeArVy5SG9fKrVpAxs2wCOPwAsv5C4hVb2605GVTJqglNd9/fXXxMbGApCUlETr1q0djsg3jDG89957ZGRk8Pe//52yZcvy4osvFsm9MzKgRw+7FUZSEoSFFcltFXYl+Bkz4L33YNAgaNbMzqtq397pyNzhlVde8dq1CtTFZ4ypYoyZYYzZZozZaoxp47UIVIny1Vdf0b59eypXrszKlStLbHI6x8/Pjw8++IA+ffrw17/+lXfeecfn9xSBxx6zmw9OmnRxmbTyPWPsIryrV9t1ATt0gJdfhqwspyMrWQo6BjUWSBSRm4BmwFbfhaSKq3HjxnHvvffSpEkTVq5cScOGDZ0OqUj4+/szZcoU7r33Xp555hnGjx/v0/v97//abqWXX4b77/fprdRlNG8O69fbKr9XX7WJKjXV6ahKkLxm7/72ACoBP3J2992CHLqSROni8XhkyJAhAkhsbKycKA5T730gIyNDYmJiBJDJkyf75B7Tp9sVDv74R7sagnKPKVNEKlQQCQkRmTOnaO+9ZcsW8RSDXwiPx3NFK0kUpAV1A3AA+NAYs8EYM9EYc9HeA8aY/saYdcaYdQcOHPBeBlWudubMGR566CHeeOMN+vfvz6xZs4rt1hSFFRgYyIwZM4iMjOTPf/4zn332mVevv3Yt9O1rB+onTdJycrfp29e2pkJD7VqIgwbBmTNFc++goCDS0tIcqyYtCBEhLS3tivZ5M5f7CxljWgCrgbZid9cdCxwTkUuWLLVo0UK8MYtYuduxY8fo2bMnCxcu5LXXXmPo0KElooy8sE6dOkVUVBTJyclMmzaNHj16FPqae/faeTdBQfD113ZLCeVO6em2wu/dd+1q8lOnQiFX/LmszMxMUlNTSU9P9+2NCikoKIjQ0NCLFlw2xqwXkRYXPSGvZpWc38VXA9j9m+/vAuLze4528ZV8P//8szRr1kzKlCkjH374odPhuM6xY8ekTZs2EhAQIHMK2d9z/LhIs2YilSqJbNrknfiU782cKVKlikjFiiJTpzodjbtxtV18IvJfYK8x5sazD3UAthQmi6ribcuWLbRu3ZqdO3cSFxfHn/70J6dDcp2KFSuSkJBAs2bN6NmzJ/Pnz7+q62Rnwx//CJs22Ym4N9/s5UCVz9x7r10mqUkT6NPHruN36pTTURUvBa3iewr41BjzHdAceMNnESlXW7ZsGW3btiUzM5Nly5bR+cLd41SOypUrM2/ePG666Sa6d+/OkiVLrvgaL7xgJ+G+887FG/Up96tbF5YutQv5vv++7abdvNnpqIqPAiUoEdkoIi1EpKmIdBeRw74OTLnP9OnT6dixIzVq1GDVqlXccsstTofkelWrVmXhwoWEhYURGxvLypUrC/zc8ePh7bft4q9PPOHDIJVPBQTAm2/CvHlw4AC0bGmTlYvrGVxD1+JTBfL222/Tu3dvWrZsSXJyMvXq1XM6pGKjevXqLFy4kOuvv56oqCjWrl172ecsWAADBkB0NLz1VhEEqXyuUye7V1fbtra7749/tDv+qkvTBKXy5fF4GDhwIIMGDaJHjx4sXLiQqlWrOh1WsVOzZk2SkpIICQmhc+fObMxnw56tW6FXL7uF+dSp4LJ9EVUh1KhhW1Kvv26XS7rlFtCC50vTBKUuKT09nT59+jBmzBiefvppPv/88yuaw6DOFxoaSlJSEsHBwXTs2JHNeQxGHDhgF4ANCrJjTxUrOhCo8ik/PxgyxI5NZWbapareflu7/PKiCUrl6fDhw3Tu3Jnp06czatQoxowZ47otzoujevXqsWjRIgICAujQoQMpKSk5P8vIsJVfv/xit2yvU8fBQJXPtW1rq/yio+2k3q5d4eBBp6NyF01Q6iJ79uyhbdu2rF69ms8++4znnntOJ+B6UcOGDVm0aBEej4eIiAh27dqFiN3CITkZpkxx72Z4yruqVoVZs2yV5vz5dm2/Zcucjso9NEGp82zcuJE2bdqwb98+5s2bR58+fZwOqURq3LgxCxcu5PTp00RERPD880f45BMYPtyOP6nSwxh46ilYtQrKlbPbdrz6qp0DV9ppglI5FixYQHh4OP7+/iQnJ3P33Xc7HVKJ1rRpU+bPn8/+/e15660q3HffKYYMcToq5ZRbb4VvvrHVfS+/DJGRsG+f01E5SxOUAuDjjz8mOjqaevXqsXr1am7WJQuKRGbmbWRnf4CfXzLffdeaAwf2Ox2SclDFivDxx/Dhh7Bmjd0MMSHB6aicowmqlBMR3njjDfr160d4eDjLly+nVq1aTodVKuzeDd26Qe3afsye7Udq6k4iIyNJS0tzOjTlIGPgT3+y5ec1a9oiihdeKLqV0d1EE1QplpWVxRNPPMHQoUN54IEHSEhIoHLlyk6HVSocO2a3ZDhzxpaT33NPG7766itSUlLo1KkTR44ccTpE5bDGje3K9Y8/DqNGwV13wY8/Oh1V0dIEVUqdOnWKHj16MG7cOAYPHsyUKVMIDAx0OqxSISsLeveGbdvsZM2bbrKPd+jQgZkzZ/L999/TpUsXjh8/7mygynHlysG//w3Tp8P27XZi74wZTkdVdDRBlUIHDhwgIiKC+Ph4/u///o8333wTPz/9VSgqAwdCYqJ94+nQ4fyfRUdHM23aNNatW0dMTAwnT550JkjlKvfdBxs22A8zvXrZVtXp005H5XsFflcyxvif3VE3zpcBKd/auXMnd9xxB99++y0zZ87kCV2FtEi9+649nnsOHn0073O6d+/Op59+SnJyMl27duV0aXgnUpcVFgbLl9vxqHHjoFUruyxWSXYlH5ufAUr4y1GyrVmzhjZt2nD48GGSkpLo1q2b0yGVKgkJ8MwzdsWAESPyP7d3795MnjyZxYsX07NnTzIyMoomSOVqAQHwz3/a36X//tfu2DtpUsldJqlACcoYEwrEABN9G47ylbi4ONq3b09wcDArV66kTZs2TodUqmzaZMedmjaFTz8t2AKwffv2Zfz48SQkJNC7d28yMzN9H6gqFrp0scsktWoFf/kLPPgglMQhy4K2oMYALwKeS51gjOlvjFlnjFl34MABb8SmvGTChAl069aNxo0bs2rVKho1auR0SKXKr79CbKyd4zJnDgQHF/y5jz76KP/617/48ssveeCBB8jKyvJdoKpYuf56uy3Lq6/aVe/PTfQtSS6boIwxscB+EVmf33kiMuHspoYtqlev7rUA1dUTEYYNG8Zjjz1Gly5dWLJkCdddd53TYZUqp0/buU4HDtgFYENDr/waTz75JKNGjWL69Ok8/PDDZOsaOOosf38YNgyWLLG/a23a2HX9SkyXn4jkewBvAqnAbuC/wCngk/yec9ttt4ly1pkzZ+Shhx4SQP7yl79IZmam0yGVOtnZIr17ixgjMnNm4a83fPhwAeSRRx6R7Ozswl9QlSgHD4rExoqASLduImlpTkdUcMA6ySv/5PXgpQ7gbiDucudpgnLWsWPHpFOnTgLIK6+8Ih6Px+mQSqVhw+z/sBEjvHfNl156SQAZMGCA/ruqi3g8Im+/LRIQIFK7tsjy5U5HVDCaoEqJffv2SfPmzcXf318++OADp8MptT7+2P7v+vOf7ZuGt3g8Hnn++ecFkEGDBmmSUnlau1akfn0Rf3+R4cNFsrKcjih/XklQBT00QTljy5YtUrduXalQoYIkJCQ4HU6ptXy5SGCgyN13i2RkeP/6Ho9HnnzySQFk6NCh3r+BKhGOHhXp08e+y3foIPLLL05HdGmXSlBlimqsS/nWihUr6Nq1K4GBgSxdupTbbrvN6ZBKpV277K64devCF1+AL1aPMsYwduxYMjIyeP311wkKCuKll17y/o1UsVapEvznP3bbjqeesiujT5kCnTs7HVnB6fo2JcAXX3xBZGQk1157LatWrdLk5JAjRyAmxm40Fx9vd0v1FT8/P8aNG0ffvn0ZNmwYI0eO9N3NVLFljJ0ntXYtVK9u508NHgzFZUqdJqhibuzYsfTq1YvbbruN5ORkwsLCnA6pVMrMtGuk7dwJM2dCw4a+v6efnx+TJk2id+/evPjii/zrX//y/U1VsXTzzXZ/qf797Som7drBnj1OR3V5mqCKKY/Hw/PPP8+zzz5L9+7dWbhwISEhIU6HVSqJ2C6UhQth/Hgoyo2Iy5Qpw8cff0z37t15+umnmTBhQtHdXBUr5cvb38+pU2HzZmje3H6YcrW8BqYKe2iRhO8cPXpUpk+fLlFRUQLIk08+KVluL9Ep4d56yw5E//WvzsWQnp4u0dHRYoyRyZMnOxeIKhZ27hRp0cL+3g4YIHL6tLPxoFV8xdf27dtl9OjREhERIWXKlBFArrnmGhk9erSWGTvsq6/sRNwePezEXCedPn1aIiMjxc/PTz777DNng1Gul5EhMmiQzQLNmols2+ZcLJdKUMb+zLtatGgh69at8/p1S4szZ86wfPly4uLiiIuL44cffgCgSZMmxMTEEBsbS+vWrSlTRoswnfTtt9C2rd2jZ9ky24XitJMnTxIdHU1ycjLTp0/n3nvvdTok5XLx8fDQQ5Cebvco69ev6GMwxqwXkRYXPa4Jyh1+/fVX5s6dS3x8PPPnz+f48eOULVuWiIgIYmNjiYmJoW7duk6Hqc765Re4/Xb79ddf24U73eL48eN06tSJ9evXM2vWLGJiYpwOSbnczz/DAw/A0qXQt69NVFeyqHFhaYJyGY/Hw4YNG4iPjycuLo61a9cCUKtWrZxWUkREBBUqVHA4UnWhU6dsFdTWrbBihR1sdpsjR44QGRnJpk2bmDNnDh07dnQ6JOVy2dnw2mv2aNAAPv+86H63NUG5wIkTJ1i4cCFxcXHMnTuXX375BWMMrVq1ymklNWvWDGOM06GqS/B44A9/sNVPX34J99zjdESXlpaWRkREBDt27CAhIYF27do5HZIqBpYssa2ptDQYPRqeeMLOp/KlSyUoHcTwsZ07dxIfH098fDxLlizhzJkzVKpUiS5duhATE0NUVBS6PUnx8dJLdoWI0aPdnZwAQkJCWLBgAXfffTcxMTHMnz+fO+64w+mwlMvdfbfdDPFPf4Inn4RFi+CDD+Caa4o+Fm1BeVlmZibJyck5XXfbtm0D4Kabbsrpumvbti0BAQEOR6qu1OTJ8PDDdrLjuHG+/1TpLb/88gvh4eHs37+fRYsW0aLFRR9UlbqIxwNvv21Xnrj+ejt/ylcbcWsXnw8dPHiQhIQE4uLimDdvHkePHiUwMJB27drldN3Vr1/f6TBVISxdCh072rGnuXOhuH2+2Lt3L+Hh4Rw9epTFixfTrFkzp0NSxcSaNdCnD/z0EwwfDi++CH5eXuLhqhOUMaY2MAWogd3yfYKIjM3vOSU9QYkI3333HXFxccTHx7N69WpEhBo1ahATE0NMTAyRkZFUrFjR6VCVF+zYAa1bw7XXwqpVUKWK0xFdnR9//JHw8HDS09NZunQpv/vd75wOSRUTR4/Co4/C9OnQqZNddNabm3MXJkHVBGqKyDfGmIrAeqC7iGy51HO8kaD277dvCG5x6tQpFi1alDOelJqaCkDLli1zuu5uueUW/Lz90UI56tAh262RlmbLyYt7QzglJSWnWGLp0qU0atTI4YhUcSEC778PzzxjP6R9/LFdKd0bvNbFZ4z5EnhXRBZc6pzCJqiff4bQUFviGBtrV4hu2RL8/a/6kldlz549OWNJixcvJj09neDgYDp16kRsbCxRUVHUqFGjaINSRebMGbv6c3KyHSi+806nI/KOLVu20K5dO44cOcINN9xAw4YNadSoUc7RsGFDatWqpR+2VJ6+/x5694Zt22DaNLjvvsJf0ysJyhhTD1gGNBGRYxf8rD/QH6BOnTq37SnEUrkHD8KkSXaGc3Kyrc+vXh2iomzC6tQJKle+6stfUlZWFqtXr87putu0aRMADRo0yBlLCg8PJ9AXm/woVxGxXRoffGA/KT74oNMReVdKSgofffQRKSkp7Nixg5SUFE6fPp3z83LlyuUkrgsTWEhIiE6FKOVOnoTXX4e//Q28MZJR6ARljAkGlgKvi0i+a+B6cwzq0CGYN88mq4QE+32ZMnDXXbZlFRsLjRpdfUXVoUOHSExMJD4+noSEBA4fPkyZMmUIDw/P6brTbpDSZ+RIOxj80kt24mJJ5/F42Ldv33kJ69yxa9cusrKycs6tUqXKRS2uc3/quKu6GoVKUMaYACAOmCcib13ufF8VSWRnw+rVEBdnE9b339vH69fP7QoMD4eyZS99DRFh8+bNOV13K1euxOPxUL16daKjo4mNjaVjx45U9kUTTRULs2dDjx52f6fPPvN+xVJxk5mZyZ49e85LWueS2E8//XTeuTVr1syz1XXDDTdQNr//mKpUK0yRhAE+Ag6JyLMFuVlRVfH99JNNVHFxkJRkFzsMDrblwLGxEB0NNWpAeno6ixcvzum6O9f9eMstt+R03bVs2VL73BXffGNb57//PSxeDOXKOR2Ru50+fZoffvjholbXjh072L9/f855fn5+1K1b96JWV6NGjahTpw7+RT3ArFylMAnqTmA58D22zBxgiIjMvdRznCgzP3XKJqlzCetskR2VK+/g5MlpZGXNoly5bXTqFElMTAzR0dHUqlWrSGNU7vbzz3YB2DJl7NwPb5bRlkZHjhw5L3H99uvjx4/nnBcYGEj9+vXz7DasUaOGjneVAle91JGIrABc/xtStmw21aqtJSQkjpCQeFJTs4FY0tN7kJX1N2AolSoJISGGkBCoVMnpiJWbnDhhly46ftwW5mhyKrwqVarQsmVLWrZsed7jIsKvv/6aZ6srMTGRjIyMnHODg4Pz7DJs2LAh1zix9o4qUsV6JYkjR44wf/584uPjmTt3LgcPHsTf35+2bdvmdN01btyYtDRDYqJtWc2bB0eO2JUA2rXLHbtq0MDn4SqXys6Gnj1hzhz7OxIV5XREpVd2djZ79+7Ns1hj9+7deDyenHOrVauWZ5dhgwYNKO+GzblUgZWIpY5EhO3bt+cUOKxYsYKsrCyqVq1KdHQ0MTExdO7cOd9PVpmZsHJlblfg1q328RtvzK0KvPPO4reUjbp6L7wAo0bBO+/AU085HY26lIyMDH788cc8izX27dt33rmhoaF5dhmGhYXpOpguVGwTVEZGBsuWLcspcNi5cycATZs2zSkDb9Wq1VUPsu7aZZNVfLwdFD9zxnb/de5sE1ZUlLtWtFDeNXGine80YAC8+67T0airdeLECXbs2HFRqyslJYXDhw/nnOfv73/JycmhoaFaKOWQYpWgjh49yowZM3J2lz158iRBQUF06NCB2NhYoqOjqVOnjhcjtk6csCsGnCtj/+UXO7/q9ttzuwKbNy8+q1ir/C1aZFeKiIy03XtldPOZEiktLS3PQo0dO3Zw6tSpnPOCgoIumpzco0cPnXJSBIpVgtq9ezdhYWHUrl07Zyypffv2RdqvLAIbNuR2Ba5dax+rVcuWr8fGQocOoBveFk/bttk19mrVskUR+h5U+ohIzuTkC7sNd+7cSVZWFqmpqVrtWwSKVYIC2LZtGzfeeKNrSkx//dWuZBEfbwstjh+3E4Lbt7ctq5gYCAtzOkpVEAcP2tXJjx+3C8DWq+d0RMptsrKy2L17N/Xr13fNe1BJVuwSlJudOQMrVuR2Baak2Md/97vcrsA77tAuIzfKyLATudessWOOvtqATSlVcJqgfGjHjtyuwGXLbKVglSp2fCM21v4ZEuJ0lErEbmM9ZYpdwqhPH6cjUkqBJqgic+wYLFiQWxm4f79dy61Nm9wy9iZNtNDCCW+8AUOHwj/+AX//u9PRKKXO0QTlAI8H1q/P7Qpcv94+Xrt2bldgRISu91YUpk+HP/wBHnjAbp+hHxCUcg9NUC6wb58ttIiLs62skydtcoqIyE1YtWs7HWXJs2aNXTXk1lttaXlQkNMRKaV+SxOUy2RkwNKluWNXu3bZx5s2ze0KbNWq6HcRLml++snOYytf3lbsVa/udERKqQsVdj+oLsBYwB+YKCL/m9/5mqCujAhs357bFbh8uV0fLiTErmQRE2NXttC1Ma/M8ePQti3s2QOrVtkqS6WU+xRmuw1/IAXoCKQCa4E/isiWSz1HE1ThHDkC8+fbhJWQYOft+PvbN9tzXYGNG+s4Sn6ys6FbN0hMtK9hx45OR6SUupTCJKg2wCsi0vns938DEJE3L/UcTVDek51tx1DOdQV++619PCxMJ5jm5/Bh2LgR3nsP/ud/nI5GKZWfq94PCqgF7P3N96lAqzxu0B/oD/hknbzSyt/flqi3aQPDh8PevTB3rm0ZpKU5HZ17VawI//ynJielirOCJKi8OpIuanaJyARgAtgWVCHjUpdQuzY89pg9lFKqJCvI2vKpwG+Ln0OBfZc4VymllPKKgiSotUBDY0yYMSYQ6AN85duwlFJKlXaX7eITkSxjzJPAPGyZ+SQR2ezzyJRSSpVqPpmoa4w5AOzxwqWqAQe9cJ2SSl+f/Onrkz99ffKnr8/lees1qisiF02j90mC8hZjzLq8Sg+Vpa9P/vT1yZ++PvnT1+fyfP0aFWQMSimllCpymqCUUkq5ktsT1ASnA3A5fX3yp69P/vT1yZ++Ppfn09fI1WNQSimlSi+3t6CUUkqVUpqglFJKuZIrE5QxposxZrsx5gdjzGCn43EbY8wkY8x+Y8wmp2NxI2NMbWPMYmPMVmPMZmPMM07H5CbGmCBjzBpjzLdnX59/OB2TGxlj/I0xG4wxcU7H4jbGmN3GmO+NMRuNMT7busJ1Y1BXs/9UaWOMCQdOAFNEpInT8biNMaYmUFNEvjHGVATWA931d8gyxhiggoicMMYEACuAZ0RktcOhuYoxZhDQAqgkIrFOx+MmxpjdQAsR8elEZje2oG4HfhCRXSJyBpgKdHM4JlcRkWXAIafjcCsR+UVEvjn79XFgK3bbGAWIdeLstwFnD3d9UnWYMSYUiAEmOh1LaebGBJXX/lP65qKuijGmHnAL8LXDobjK2e6rjcB+YIGI6OtzvjHAi4DH4TjcSoD5xpj1Z/cC9Ak3JqgC7T+l1OUYY4KBL4BnReSY0/G4iYhki0hz7PY5txtjtKv4LGNMLLBfRNY7HYuLtRWRW4EoYMDZYQevc2OC0v2nVKGdHVv5AvhURGY6HY9bicgRYAnQxdlIXKUt0PXsOMtUIMIY84mzIbmLiOw7++d+YBZ2aMbr3JigdP8pVShniwA+ALaKyFtOx+M2xpjqxpgqZ78uB0QC2xwNykVE5G8iEioi9bDvP0ki8qDDYbmGMabC2eIjjDEVgE6ATyqKXZegRCQLOLf/1FZgmu4/dT5jzGfAKuBGY0yqMeYvTsfkMm2BvthPvhvPHtFOB+UiNYHFxpjvsB8IF4iIllKrgroOWGGM+RZYA8SLSKIvbuS6MnOllFIKXNiCUkoppUATlFJKKZfSBKWUUsqVNEEppZRyJU1QSimlXEkTlFJKKVfSBKWUUsqV/h8cfzW+MxmEgQAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x144 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAkcAAADbCAYAAABneDXmAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAAAwKklEQVR4nO3de5RU5ZX38e+mAQG5thiVS4KaiBcUNLyOSpQgHRTj4DjRRFfimOiMr28mxkzISMzIksxr5g1jRjSTK4kakxg1El2jjoiCrRhHjQ3iDYwag9JglIuAKIp07/ePOoVF29V1r/OcU7/PWr3oOl3n7KeK6t377HN5zN0RERERkYxecQ9AREREJCQqjkRERERyqDgSERERyaHiSERERCSHiiMRERGRHCqORERERHKoOBIRERHJoeJIREREJIeKoxxm9qyZfbKHn682s5b6jWi32D2OrYLtxvaa8qnVa622EN87kXwaMb9F2w7u91Q5LnyJLY7MbJCZ/ZuZvWhmb5rZn83sB2a2d7nbdPfD3P2BaPtxJooPxM4dW9qF+FobOUlIfZnZMDNzM9sWfb1iZgvMbFwl2+36exVSjgvxd76WQny9ynG7S2RxZGZDgYeAg4Hp7j4IOB7oA3wkxqFJD8ysd9xjEEmACcAmdx/o7gOBI4EngcfM7OBYRyY9Uo5Lj0QWR8A8YBNwhru/AODu7e7+v929LfeJZvYlM7sz5/GLZvbbnMdrzGxC9P1qM2sxs18BHwbujPbcLsnZ5AQze8rMtpjZLWbWr7sBmtkhZvaAmW2OWqgzcn622swuNbOVZvaGmV2f3U6+2LlVffT9P0fjeMvMrjWzfcxsYdRFW2xmw3LifdPM/hT9bKWZnV7oDS7xfcu7/Wiss8zsKeAtM+vd0+vP81q/ke89N7OjzOyJKPat0c+vyPOaCsXt9nVU4/MgUoIJwIrsA3ff6O7/F1gOnN/1ycX+rnb5vcr3mS7681zNHGdduhbKccpxsXP3RH0Bo4GdwLFFPv8AYDOZQnA/4GVgbc7P3gB6RY9XAy1dv8/Z1mrgD8AIoBlYBVzYTcw+wIvAt4C+wInAm8DYnO08E72WZuBh4IoucbqLnTu2R4F9gJHA62QS55HAHsD9wOU5654ZjbkX8DngLWC/fLHKeN8KbX9F9Fr7l/r6e3rPo/f2ZeDi6D3/W2BH7ra6eQ97ilvS+1Ts50Ff+irlC/gl8B/dLP8Z8Jtulhf1u9r1M5zncVGfZ6qc4/I8Vo5TjovtK4mdoxZgvbs/UsyT3f0lMr+0E4DJwCJgrWXa05OBh9y9s4T433f3de6+Cbgz2m5XxwADge+6+w53vx+4Czg75zk/cPc10Xa+0+VnxfhPd3/N3deSOcT4mLs/4e7vAreTSSIAuPut0Zg73f0W4AXg6J42Xsr7VsT2vx+91u1lvv587/kxQO/o5++5+21kfpF7kjduOe9TD2MTKdcEcjpHOYYA67surHKOK/bzrBynHJdqSTw+ug/wSonrPAh8Evho9P1mMh/+Y6PHpfhLzvdvk6mouxoBrOmSkF4msweUtabLz7rbTk9ey/l+ezePB2YfmNnfAV8HxkSLBgLDi4hR1PtWxPZzX2t3ywq9/nzv+Qgye3peIFZRcct8n4r5PIgUxcz2AA4hc45R7vIm4Djgq3lWrVaOK/bzrBynHJdqSewcvQKMNLNSxp79BTg++v5BMr8Ak8mfODzP8mKsA0Z3GeOHgbU5j0d3+dm6KsXejZl9hEw7/ivAXu4+lEzb1YpYveD7VuT2u3s9Pb3+Yr1K5rOQG2t0vif3FLeI11G1/xORHowDOskcvsh1IZnDKXd+YI0M5TjluB7jKseVJonF0V3Rv981s8Fm1sfMDo9O2Mt3Gf+DwBQyx4PbybRoTwb2Ap7Is85rZI47l+MxMsdyL4nG90ngr4Gbc57zj2Y2ysyayRy3v6VKsbvak8yHfj1kTkIkk4CLUcz7Vu72e3r9xXoE6AC+Ep0EeRqFW8T54hZ6HdX8PxHJ50jgWXd/D8DMRkcn334bOCu7vBvKccpxheIqx5UgccWRu28lc/LfQWSOl24k8wv5mrt/4Hh8tM7zwDYyH/zsNl4CHnb3jjyh/h9wmWWuxPhGiWPcAcwApgMbgB8Bf+fuz+U87TfAvdE4XgJyrz4oO3Y3Y1kJ/AeZX7LXgMPJnKRXzLoF37cKtt/T6y9K9D7/LZkreDYDXyBTPL9batwiXkfV/k9EejABOCK6mugNYDEwDJjo7nnPNVGOU44rFFc5rjS2+6FMqQczWw38vbsvjnsscajl6zezx4CfuPv19YwrIu9r9N815bjkS1znSCSXmU02s32jlvO5wBHAPXGPS0SkGpTj4pHEq9VEco0Ffkvmqos/kbkx6KvxDklEpGqU42Kgw2oiIiIiOXRYTURERCSHiiMRERGRHEGdc9Tc3OyjRo2KexgN56WXXootdt++fWOLDRDX5629vZ1NmzYVc5O63ZhZoePgi9z95DKHJTU0fPhwHzNmTNzDiMUzzzwTW+x33+3pqvfa2mOPPWKLDTB06NBY4m7dupW333675PwG4eS4oIqjUaNGcffdd8c9jIZz5plnxhZ79OhCN3utrauuuiqWuKecckrZ6/bqlb/h29nZWcyUCRKDMWPG0NbWFvcwYnHwwQfHFvuPf/xjbLHjLoZPO+20WOLecMMNFa0fQo4LqjgSkcIKJI46jkREpPpCyHEqjkQSxMzYfZolEZH0CCXHqTgSSZie9qpERJIuhByn4kgkYUJIHCIitRJCjlNxJJIgobScRURqIZQcp+JIJGFC2KsSEamVEHJc/CMQkZJk96y6+xIRSbpKcpyZXWdmr5vZMznLrjSz58zsKTO73cyGFtqOiiORBDEzevXqlfdLRCTJqpDjfgF0vUnkfcA4dz8CeB64tNBGdFhNJGFUBIlImlWS49x9qZmN6bLs3pyHjwJnFNqOiiORhNHhMxFJsxrnuPOAWwo9qebFkZk1AW3AWnc/tdbxRGqt/+23M2TuXJrWraNjxAi2zJrF9tNPr0vsbMtZwqD8JmkzYeVKTl66lKFbt7J58GDuOeEEVhx6aN3iF5HjhptZ7jw88919fpHb/hdgJ3BjoefWo3N0MbAKGFyHWCI11f/22xk2axa9tm8HoPfatQybNQugbgWSiqOgKL9JakxYuZLP3HMPfXfuBGDY1q185p57AOpaIBXIcRvcfWKp2zSzc4FTganuXmhy29qekG1mo4BPAz+vZRyRehkyd+6uwiir1/btDJk7t25jqPRqNTP7JzN71syeMbObzKxfjYecSspvkjYnL126qzDK6rtzJycvXVrXcVT7ilwzOxmYBcxw97eLWafWu6BXA5cAeWeKM7MLzKzNzNo2bdpU4+GIVKZp3bqSlldbpVdymNlI4KvARHcfBzQBZ9V42Gl1NSXkt/Xr15ccYOHChUyZMoXZs2fT2trK5MmTmTlzZvkjTkhsicfQrVtLWl4LVchxNwGPAGPNrN3Mzgd+AAwC7jOzFWb2k0LbqVlxZGanAq+7+7Kenufu8919ortPbG5urtVwRKqiY999u18+YkTdxlCFS/l7A/3NrDcwAKhPZZci5eS3vffeu+Q4S5YsobW1lUGDBjFv3jwWL17M+PHjWblyZblDT0Rsicfmwd0fHc63vFYqyXHufra77+fufdx9lLtf6+4fdffR7j4h+rqw4Biq8kq6NwmYYWargZuBE83s1zWMJ1JzO8aNo+vB6s7+/dkSnXdUDwVazsOznYro64Lcdd19LfA94BXgVWBLl8tcpTh1yW/ZUyOam5vZtm0bHR0duDtFnDKR6NgSjxUHH/yBZTt69+aeE06o6zhCuNFtzYojd780qtrGkGnb3+/uX6hVPJFa6/XGG/T7n/9hx5FHsnPkSNyMnSNH8sbcuXW/Wq2HvaoN2U5F9DW/y/rDgNOA/YERwJ5mpt/LEtUrv7W0tDB16lRWrVrFnDlzmDZtGm1tbRx22GHVDhVUbImBOwe0t/PWHnvwxqBBOPDG4MH87uSTY7laLe4b3eo+RyJFGjh/Pvb227xx5ZXsHDs2tnFUuPfUAvzZ3ddH27oNOA5QVzdA06dPZ/r06bseL63jibFxxpb6O2j1aj6ybh23TZvGYxMmxDqWEO7lVpfiyN0fAB6oRyyRWuj1xhsMvP56tp96aqyFEVR8Kf8rwDFmNgDYDkwlc58eKZPymySeOy0PP8zmQYNoO/zwuEcTxO1K1DkSKUK2a7T14otjHUelN4F098fMbAGwnMzN0J4AirqBmoikU27XqKOpKdaxhHKjWxVHIgWE1DWCylvO7n45cHl1RiMiiRZY1wga6LCaSJKF0jXKCmGvSkTSIaSuUVYIOU7FkUgPQuwahZA4RCQFAu0ahZDj4h+BSMBC6xpBGPcAEZHky3aN7j/22GC6RhBGjlPnSCSP0LpGWSHsVYlIwgXYNcoKIcepOBLJI9SuUQiJQ0SSLcRzjSCcHKfiSKQboXaNIIwrOUQkwQLuGkEYOU7FkUg3QuwaZYWwVyUiyRVq1ygrhBwXVHHUt29fRo0aFUvsRx99NJa4IRg9enRsseP6/+4x/saN8ItfwGc/y75Tp9Ykbt++fcteN4S9KpFSnHfeebHFvuOOO2KLHbex3XW93Tl1wQK2NTez9TOfYWzv6pcB/fr1q2j9EHJcUMWRSBCuugreegtmz457JB8QyvF4EUmmkc8+y4deeomHzzmHzhoURpUKJcfFPwKRkGzcCN//Pnz2sxDo7OMhzFgtIgnkzpH/9V9sa27mhU98Iu7R5FVJjjOz68zsdTN7JmdZs5ndZ2YvRP8OKziGCl+DSLoE3DXKCuEeICKSPNmu0ZOf/nSQXaOsCnPcL4CTuyz7JrDE3T8GLIke90jFkUhWArpG2ZazOkciUpKEdI0qzXHuvhTY1GXxacAN0fc3AH9TaDvhlo4i9ZaArhGEcSWHiCRL6Oca5apBjtvH3V8FcPdXzexDhVYI+x0SqZcEdI2ydPhMREqSkK5RVoEcN9zM2nIez3f3+dUeg4ojEUhM1yiUKzlEJDmS1DUqIsdtcPeJJW72NTPbL+oa7Qe8XmgFZVmRBHWNQFeriUgJEtY1gprkuDuAc6PvzwX+q+AYyo0kkhoJ6Rpl6Wo1ESlWUq5Qy1VJjjOzm4BHgLFm1m5m5wPfBT5lZi8An4oe9ygZ75RIrSSsa6TDaiJStAR2jSrNce5+dp4flTTdgbKsNLaEdY1AnSMRKU4Su0YQRo5LzrslUmUDtm+Hn/0sMV2jLHWORKSgBHaNskLIcTUrjsysH7AU2COKs8DdL69VPJFSTV6+PJFdoxASR6NTfpPQHb15c2KuUMsVSo6r5Tv2LnCiu28zsz7A781sobs/WsOYIkUZsH07xz/xROK6RqD7HAVC+U3C5c55a9YksmsEYeS4mhVH7u7Atuhhn+jLaxVPpBSTly+n73vvJaprlBXCXlWjU36TkB29eTPj3nwzcV2jrBByXE1HYGZNZraCzA2X7nP3x7p5zgVm1mZmbevXr6/lcESA97tGTx50UCK7RrrPURjqkd8WLlzIlClTmD17Nq2trUyePJmZM2dWPvjAY0sFoq7Ra337JrZrFEKOq2kkd+9w9wnAKOBoMxvXzXPmu/tEd5+4995713I4IsD7XaN7/+qv4h5KWUK4kkPqk9+WLFlCa2srgwYNYt68eSxevJjx48ezcuXKyl9AwLGlfNmu0S9Hj05k1wjCyHF1KcPcfTPwAHByPeKJ5JPbNXpt+PC4h1OWSveqzGyomS0ws+fMbJWZHVvjIadaLfNb5ugdNDc3s23bNjo6OnD3XctrKc7YUqacrtF/f6jg3KrBSnXnyMz2NrOh0ff9gRbguVrFEylGGrpGVUgc1wD3uPvBwHhgVc0GnFL1ym8tLS1MnTqVVatWMWfOHKZNm0ZbWxuH1eFwcJyxpTy5XaOdCT3MHsphtbw9NzObCVzt7h1dlu8F/Lu7n19g2/sBN5hZE5ki7LfuflelAxYpVxq6RlDZlRxmNhg4AfgigLvvAHZUZWAJkpT8Nn36dKZPn77r8dKlS6sdIsjYUoaUdI0gjKvVeirDxgLLzGxSdoGZfRloA54utGF3f8rdj3T3I9x9nLv/a+XDFSlf0rtGUJW9qgOA9cD1ZvaEmf3czPas7aiDpPwmqZKGrhEkoHPk7heY2XHAD8zsWeBg4AXgOHd/tV4DFKmGtHSNoOBe1XAza8t5PN/d5+c87g0cBVzk7o+Z2TXAN4Hk3dOgAspvkiop6hpBGJ2jQqeyPwM8TuZEQwNmKnFIEqWha5RVYO9pg7tP7OHn7UB7zmXnC8gUR41I+U1SIds1uvLAAxPdNcoK4bYkeUdgZl8AVgAvAQcCpwP/bma/NLPkl6bSMNLWNaqk5ezufwHWmNnYaNFUoOGuzVZ+k9RIYdco6MNqwJnAFHd/OXq8LLrk90LgUTLnLogEL01dI6hKy/ki4EYz60umOPhSxYNKHuU3SYW0dY0g8MNq7n5aN8sc+LGZLajpqESqJE1do6xK957cfQXQ06G31FN+k1RIWdcoq5IcZ2b/BPw9mel8nga+5O7vlDyGcoK7u+b5kERIY9cohJZzmim/SVKk5Qq1XJXkODMbCXwVmOju44Am4KxyxpHMe4uLFCGNXSMIo+UsIjFLadcIKs5xvYH+ZvYeMABYV85G0lFqinQjbV2jLHWORCSNXaOscnOcu68Fvge8ArwKbHH3e8saQ6EnmNkAM5ttZj+LHn/MzE4tJ5hIvaS5a6TiqHqU3ySRUt41KpDjhptZW87XBTnrDgNOA/YHRgB7RlemlqyYbHo98C6QnZyyHbiinGAi9ZLWrhGEMWN1iii/SeKkuWsEBXPcBnefmPOVe5PbFuDP7r7e3d8DbgOOK2cMxZxzdKC7f87MzgZw9+1Woyz83HPPceyx8UwQ/uijj8YSNwRr1qyJLXZ7e3vVt9l7yxaO/OlP2dTSwuQvf7nH58b1eXvuufLnKFWHqKrqlt/i9IlPfCLW+A8//HBssa+99trYYm/YsKH6G3Xn8z/6EVuHDMHOO48ZvfP/GT/vvPOqH78I559//rJK1q8gx70CHGNmA4DtZO7j1tbzKt0rpjjaEc067QBmdiCZPS2RIO1300302r6dtV9K3+17si1nqRrlN0mUMc8/z8hXXmHR6afT2UNhlFSV5LhoSqQFwHJgJ/AEML/ntbpXzDt7OXAPMNrMbgQmEc3oLRKa3lu2sO+tt7Jx6lS2H5DO+/ilsLERJ+U3SQ53Ji1ezNYhQ3h6YnpvVVZJjnP3y8n8XlekYHHk7veZ2XLgGDLzD13s7jXoFYpULs1doyx1jqpH+U2SJO1do6wQclwxV6udDux09/9297uAnWb2NzUfmUiJGqFrBDohu5qU3yQxGqRrBGHkuGLKs8vdfUv2gbtvpgotK5Fqa4SukS7lrzrlN0mEbNfokRNPTHXXKJQcV8w73N1o0vs/I4nUKF0jCKPlnCLKbxK+BuoaQRg5rpgRtJnZVWZ2oJkdYGbzgIou0xOptkboGmWF0HJOEeU3CV6jdI2yQshxxRRHFwE7gFuAW4F3gH+s5aBEStFIXaNQWs4povwmYWuwrlEoOa6Yq9XeAr5Zh7GIlKWRukYQRss5LZTfJHSNcoVarhByXMF32swOAr4BjMl9vrufWLthiRSnkbpGWTp8Vj3KbxK0BusaZYWQ44opQ28FfgL8HOio7XBEStNoXSPdIbvqlN8kWI3YNQolxxXzbu909x+XumEzGw38EtgX6ATmu/s1pW5HJJ9G7BpBGC3nFFF+kzA1aNcIwshxxRRHd5rZl4HbyZlzyN03FVhvJzDT3Zeb2SBgmZnd5+4ryx+uCOy1aBEf/slP6PuXvwDw1kEHxTyi+gqh5Zwiym8SlEOWL+eERYsYvHkzBjw5cWLDdI2yQshxxbzj50b//nPOMgd63FV391eBV6Pv3zSzVcBIQMlDyrbXokUc8N3v0vTOO7uWjbruOnZ86ENsPOmkGEdWH6G0nFNE+U2Cccjy5Zx82230ee+9XcsOffJJXjnwQFYddVSMI6ufUHJcwRG4+/7dfJV0DMPMxgBHAo+VOU4RAD784x/vVhgBNL3zDh/+yU9iGlH9hXAPkLQIPb8tXLiQKVOmMHv2bFpbW5k8eTIzZ86sdhgJxAmLFu1WGAH0ee89Tli0qC7xQ/m8hZDjiplbbYCZXWZm86PHHzOzU4sNYGYDgd8BX3P3rd38/AIzazOztve6fChEdunsZNgDD9D3tde6/XG+5WkUwj1A0qKe+W39+vUlj2/JkiW0trYyaNAg5s2bx+LFixk/fjwrV6pBlTb7tLczePPmbn+Wb3m1hfJ5CyHHFXNY7Xoyd4w9LnrcTuYKj7sKrWhmfcgkjhvd/bbunuPu84H5AAMHDvQixiONpLOTYUuXMuq669jzhRfobGrCOj54UdGOffaJYXD1F0rLOUXqlt8mTpxYcn5zz6zS3NzMtm3b6OjowN13LZfk26e9nUmLF/PRVavoNMO6+b/dOnRoXcYSwuctlBxXTHF0oLt/zszOBnD37VZEbyt6zrXAKne/qsJxSqPpUhRt//CHefHyy/HOTg648srdDq119OvHKxdeGONg60uHz6oq6PzW0tLC1KlTmTBhAnPmzGHatGmMHz+ec889t/DKErTcomh7//4sPekk3tpzT1ruvHO3Q2vv9enD0jqdTxnK562SHGdmQ8ncmmMcmfMHz3P3R0rdTjHF0Q4z6x8FwcwOJOeqjh5MAs4BnjazFdGyb7n73aUOUhpInqJoQ0sLZK/YaGrKXK322mvs2GcfXrnwwoY4GTsrhL2qFAk6v02fPp3p06fverx06dJqbl5i0F1RtPy449jRrx8AO/v02XW12tahQ1l60kl1Oxk7lM9bhTnuGuAedz/DzPoCA8rZSDHF0eXAPcBoM7uRTFL4YqGV3P33gHZxpTjFFEWRjSed1FDFUK5QWs4povwmdVGoKMpaddRRDXNlWncqyXFmNhg4geh32N13kJk7sWTFzK12n5ktB44hkwwudvcN5QQT+YASiiLJqMZhNTNrAtqAte5e9AnIaaP8JrVWbFEk76sgxx0ArAeuN7PxZM4nvDiaQ7EkxcytdkL07ZvRv4eaGe6u/q6Ur7OTfosWcfiVV6ooKlGVOkcXA6uAwdXYWFIpv0mt7LV6NScsWKCiqAwFctxwM2vLeTw/uvABMjXNUcBF7v6YmV1DZmLp2aWOoZi/Qrk3R+sHHE2mGtPEjFK6qCgafPXV9F25UkVRiapxWM3MRgGfBr4DfL0a40ow5Tepqr1Wr+bIO+7gw08+qaKoDEXkuA3unm8+lXag3d2z9xxbQKY4Klkxh9X+OvdxNKfQv5cTTBpYl6LovQMOYNM11/D8UUepKCpRgZZzT3tVWVcDlwCDqjy0xFF+k2rJLYreHTCAZaefzkMTJqgoKkO5h9Xc/S9mtsbMxrr7H4GplHnX+nL+KrWTuUROpLA8RdHbM2ZkiqL29rhHmDgV7FUR3eDwdXdfZmafrPLQ0kD5TUrSXVG0cupU3uvfnx0bdPpaOSrsjl8E3BhdqfYS8KVyNlLMOUf/SXSZK5k7ak8AniwnmDSQQkWRlKUKt9CfBMwws1PIHEYabGa/dvcvVGWACaP8JuXqqSiS8lWa49x9BZB3B7FYxfyVym3R7wRucveHKw0sKaWiqOYq2aty90uBSwGiztE3GrUwiii/SUlUFNVeCLcrKeacoxvqMRBJOBVFdRNC4kgL5Tcploqi+gkhxxVzWO1p3m877/YjwN39iKqPSpJDRVFdVXNmand/AHigKhtLKOU3KURFUX1VM8dVopi/Xgujf38V/ft54G1Ae1yNTEVRbELYq0oR5Tfploqi+ISQ44r5KzbJ3SflPP6mmT3s7v9a7cEccMAB3HrrrdXebFEeeaTkeelS46qrip8309wZ96c/Me3RRxm+fj2vDxvGfSefzBNjx9K5ejV8//slxR49enSJo62uuD5vp5xyStnrhpA4UqRu+S1Ol1xySdxDiM2MGTNKW2HZMvj2t+HOO2HYMLjiCva46CI+PngwH6/NEKWLEHJcMcXRnmb2iWguIczsOGDP2g5LQpNbFI2MiqIbs0VRAB/kRhFKyzlFlN8ko5uiiIsugsENfRP5ugslxxVTHJ0PXGdmQ8gcm98CnFfTUUkwVBSFJ4S9qhRRfmt0KoqCE0KOK+ZqtWXA+Gi2W3P3LbUflsRNRVG4QkgcaaH81sBUFAUrhBxXzNVq+wD/Boxw9+lmdihwrLtfW/PRSd2pKApbKC3ntFB+a0AqioIWSo4r5rDaL4DrgX+JHj8P3AIoeaSIiqLkCGGvKkV+gfJbY1BRlBgh5LhiiqPh7v5bM7sUwN13mllHjccl9dLZyeEvvqiiKEFC2KtKEeW3tFNRlDgh5LhiiqO3zGwvohulmdkxZE5alCTr7GREWxuH3XorQ19+WUVRQphZEHtVKaL8llJDXnyRsTffDI8/rqIoQULJccUUR18H7gAONLOHgb2BM2o6KqmdLkXRm/vtp6IoYUJIHCmi/JYy2aJo38cfZ8fAgSqKEiiEHFfM1WrLzWwyMJbMLfX/6O7v1XxkUl3dFEWPfeUrrJk0iWV/+EPco5MShNByTgvlt/ToWhSt+vzn+fOpp3LKWWfFPTQpUQg5Lm9xZGb/C1jj7n+JjsN/HPgM8LKZzXH3TXUbpZSvh6LIm5riHp2UKJSWc9Ipv6VHvqJo54ABcQ9NyhBKjuupc/RToAXAzE4AvgtcBEwA5qPWc9hUFKVWCIkjBZTfEk5FUXpVmuPMrAloA9a6+6nlbKOn4qgpZ+/pc8B8d/8d8DszW1FOMKkDFUWpF0LLOQWU3xJKRVH6VSHHXQysAso+0azH4sjMerv7TmAqcEGR6wFgZtcBpwKvu/u4cgco+Y1+6CEOv+kmBmzcyNt77UX70Uezz7PPqihKsVBazilQUX4D5bhaG/nggxzyq1/Rf8MGtg8fzupPfYrmF15QUZRyleY4MxsFfBr4DpkLLsrSUxK4CXjQzDYA24GHosAfpbhLXX8B/AD4ZbmDk/xGP/QQE3/6U3rv2AHAnhs2MPbuu9k+ZIiKopRTcVQVleY3UI6rmZEPPsj4H/6Q3u++C8CA9es59De/Yecee6goagAV5rirgUuAQZVsJG9x5O7fMbMlwH7Ave7u0Y96kTk23yN3X2pmYyoZnOR3+E037SqMcnX26cMrJ5wQw4ikXnRYrXKV5rdoG8pxNXLIr361qzDKtWPQIF747GdjGJHUU4EcN9zM2nIez3f3+dF62U7uMjP7ZCVj6LE8c/dH3f12d38rZ9nz7r68kqC5zOwCM2szs7ZNm3SBSEGdnYz4wx8YsGFDtz8esHFjnQck9ZRtOef7kuLVO7+tX7++5PUXLlzIlClTmD17Nq2trUyePJmZM2dWa3hBGvLii/TP8171r1N+i/N9b9TYWUXkuA3uPjHna37O6pOAGWa2GrgZONHMfl3OOGLPpu4+P/sim5ub4x5OuKKi6FOzZjHpe9/D8/whfHuvveo8MKm37MSM3X1JWHLz2957713y+kuWLKG1tZVBgwYxb948Fi9ezPjx41m5cmUNRhuvIS++yNFXXMHkmTMhz2d5+/DhdRlLnO97o8bOVW6Oc/dL3X2Uu48BzgLud/cvlDOGok48lBjlufoMdz7+s5/tdmhtZ9++PH322TEOVupBHaLGkT3a19zczLZt2+jo6MDdef8oYPJ1d/XZO0OHcvjPf77bobWde+zBqnPOqcuY4nzfGzV2rhBynIqjUBVxSb736rXb1WpPn302a44/PuaBSy3parXG0tLSwtSpU5kwYQJz5sxh2rRpjB8/nnPPPTfuoVWs0CX5nXvssdvVaqvOOYe1kyfXZWxxvu+NGjurWjnO3R8AHih3/ZoVR2Z2E/BJMidPtQOXu/u1tYqXGiXcp2jN8cerGGpAOnwWhnrkuOnTpzN9+vRdj5cuXVrNzcei2PsUrZ08uW7FUFdxvu+NGjtXCDmuZsWRu+v4Til080YpkjpHYVCOK41u3ijFCiHH6bBa3FQUSQl0WE2SRkWRlCKUHKfiKC4qiqRMIbScRQpRUSTlCiHHqTiqNxVFUqEQ9qpE8lFRJJUKIcepOKoXFUVSBVWYd2g0meku9gU6ydxd9poqDU8amIoiqQYdVmsUKoqkyipsOe8EZrr7cjMbBCwzs/vcPX13FZS6UFEk1abDammmokhqpJK9Knd/FXg1+v5NM1sFjARUHElJVBRJrahzlEYqiqSGimg5552UsZttjQGOBB6r3ggl7VQUSS3psFraqCiSOinQct7g7hOL2MZA4HfA19x9a7XGJumlokjqRYfV0kBFkdRZpXtVZtaHTGF0o7vfVpVBSWqpKJJ6U+coyVQUSQyKmZm6wPoGXAuscverqjYwSR0VRRKHSnNctag4ihx77LHFPbGzk36LFjH46qvpu3Il7+2/P5uuuYa3Z8xgZO/ejKztMGviqqvi+xvZ3t4eW2yAM888M9b45ahwr2oScA7wtJmtiJZ9y93vrnRc0rOdO3eycePGWGLPmDGj+CcvWwbf/jbceScMGwZXXEHfiy7ikMGDOaR2QxTZRZ2jJOmhKKK33kapnwqvVvs9EP9umYSnm6KIiy6CwYPjHpk0GBVHSaCiSAISSstZUkRFkQQklBynv+75qCiSQIWwVyUpoKJIAlVujqvmDAD6K9+ViiIJnIojqYiKIglcBTmuajMA6K99looiSYBQWs6SQCqKJAEqyXHVnAFAf/VVFEnCqHMkpWhasQKuuUZFkSRGNXJcpTMANO5ffxVFklAqjqQYTStWMODKK+m7aJGKIkmUSqdIqsYMAI1XBagokgTTYTUpJLco6hw6lLe/9S0GzJqlokgSoYgc1+MUSdWaAaBxqgEVRZIS6hxJd7orit75h3/ABw1igAojSZAKrlar2gwA6a8KVBRJyqhzJLl6KopEkqiCHFe1GQDSWx2oKJIUMjN1jgRQUSTpVEmOq+YMADWtEszsZOAaoAn4ubt/t9ox+t9+O0PmzqVp3To6RoxgyyWX4P37qyiS1FJxFIZ65Le+CxYw4Ior6LV2LZ0jR/L2ZZfR8dGPqiiSVAshx9WsWjCzJuCHwKeAduBxM7ujnJsx5dP/9tsZNmsWvbZvB6D32rU0f+1rmLuKIkktHVaLXz3yW98FCxj4T/+ERfmtqb2dgV/+MtbZqaJIUi2EHFfLquFo4EV3fwnAzG4GTqOMmzHlM2Tu3F2FUZa50zFsGK/df7+KIkkdHVYLRs3z24ArrthVGGVZZyedgwez+YknVBRJKoWS42o5gpHAmpzH7dGy3ZjZBWbWZmZtmzZtKilA07p13S7vtXmzCiNJrV69euX9kropOb9t3LixpAC91q7tdrm9+WZdCqOFCxcyZcoUZs+eTWtrK5MnT2bmzJk1jxt37LjjN2rsXCHkuFpG6q4v5h9Y4D7f3Se6+8Tm5uaSAnSMGFHScpE0yN4HpLsvqZuS89tee+1VUoDOkR+otXpcXm1LliyhtbWVQYMGMW/ePBYvXsz48eNZubJqzbEgY8cdv1Fj5wohx9WyOGoHRuc8HgV03+op05ZZs+js33+3ZZ39+7Nl1qxqhhEJRrblHPdeldQ+v7192WV4l/zm/fvz9mWXVTNMXu6ZWq+5uZlt27bR0dGBu+9antbYccdv1NhZoeS4WkZ6HPiYme1vZn2Bs4A7qhlg++mn88bcuewcORI3Y+fIkbwxdy7bTz+9mmFEghJC4pDa57cdZ5zBtnnz6Bg1CjejY9Qots2bx44zzqhmmLxaWlqYOnUqq1atYs6cOUybNo22tjYOO+ywVMeOO36jxs4VQo6zWlaEZnYKcDWZS12vc/fv9PT8I444wu++u+R7NUmFzjzzzNhijx49uvCTauiqqyq6iWrZTjnlFJ566qmSe8QTJkzwJUuW5P358OHDl/V0a32pnlLzW6H/u1oq9ZCeSCXMrOw8FEqOq+lZy9FdKVXtiFRJKFdyiPKbSC2EkuN0SZdIwujEaxFJsxBynIojkYQJYa9KRKRWQshxKo5EEiSUlrOISC2EkuNUHIkkTAgtZxGRWgkhx6k4EkmYEPaqRERqJYQcp+JIJEFCaTmLiNRCKDku/hGISEkqvbW+mZ1sZn80sxfN7Js1Hq6ISEkqyXHVym/qHIkkTCV7VWbWBPwQ+BSZKTAeN7M73L2+kyeJiORRbo6rZn5T50gkQaow79DRwIvu/pK77wBuBk6r6aBFRIpUYY6rWn5TcSSSMBUeVhsJrMl53B4tExEJQgU5rmr5raZzq5XKzNYDL5e5+nBgQxWHk5TYccdv1NiVxv+Iu+9d6kpmdk8UN59+wDs5j+e7+/yc9c8ETnL3v48enwMc7e4XlToWKU2C81vc8RU7efHLym9QWY6rZn4L6pyjct9MADNri2vCzThjxx2/UWPHFd/dT65wE+1A7my/o4B1FW5TipDU/BZ3fMWOR1zxK8xxVctvOqwm0lgeBz5mZvubWV/gLOCOmMckIlINVctvQXWORKS23H2nmX0FWAQ0Ade5+7MxD0tEpGLVzG9pKo7mF35KKmPHHb9RY4cQvyzufjdwd9zjkJLE/Vlr1N/zRo0dQvyyVCu/BXVCtoiIiEjcdM6RiIiISI5UFEdxTYdgZteZ2etm9ky9YubEHm1mrWa2ysyeNbOL6xy/n5n9wcyejOJ/u87xm8zsCTO7q55xo9irzexpM1thZm31ji+NJc7pXho1x8Wd36IxxJLjlN8yEn9YLbpd+PPk3C4cOLse0yGY2QnANuCX7j6u1vG6xN4P2M/dl5vZIGAZ8Df1mgbCMnfj2tPdt5lZH+D3wMXu/mid4n8dmAgMdvdT6xEzJ/ZqYKK7x3kPEmkAcea3KH5D5ri481s0hlhynPJbRho6R7FNh+DuS4FN9YjVTexX3X159P2bwCrqeKdjz9gWPewTfdWl0jazUcCngZ/XI55IjGKd7qVRc1yc+Q2U40KQhuKo4adDMLMxwJHAY3WO22RmK4DXgfvcvV7xrwYuATrrFK8rB+41s2VmdkFMY5DG0PD5DeLJcTHmN4g3xym/kY7iqLvJVpJ9rLAEZjYQ+B3wNXffWs/Y7t7h7hPI3IX0aDOredvdzE4FXnf3ZbWO1YNJ7n4UMB34x+jQg0gtNHR+g/hyXBz5DYLIccpvpKM4atjpEKJj4b8DbnT32+Iah7tvBh4AKp3aohiTgBnRcfGbgRPN7Nd1iLuLu6+L/n0duJ3MoQ+RWmjY/AZh5Lg65zeIOccpv2WkoThqyOkQohMGrwVWuftVMcTf28yGRt/3B1qA52od190vdfdR7j6GzP/1/e7+hVrHzTKzPaOTQzGzPYFpQN2v5JGG0ZD5DeLNcXHlN4g3xym/vS/xxZG77wSytwtfBfy2XtMhmNlNwCPAWDNrN7Pz6xE3Mgk4h8xexYro65Q6xt8PaDWzp8gk8Pvcve6X1cdgH+D3ZvYk8Afgv939npjHJCkVZ36Dhs5xym8Nnt8Sfym/iIiISDUlvnMkIiIiUk0qjkRERERyqDgSERERyaHiSERERCSHiiMRERGRHCqOAmNm/xLNAv1UdOnqX8U9pnKZ2QVmdkvO48Fm9icz2z/OcYlIPJTfJClUHAXEzI4FTgWOcvcjyNx4bE3PawXtZ8AoM2uJHv8rcJ27/znGMYlIDJTfJElUHIVlP2CDu78L4O4bsrdyN7OPm9mD0WSAi8xsv5zlT5rZI2Z2pZk9Ey3/opn9ILthM7vLzD4ZfT8tev5yM7s1mrsIM1ttZt+Olj9tZgdHywea2fXRsqfM7DM9bSfLMzfR+j/A1WY2EZgKXFnD909EwqX8Jomh4igs9wKjzex5M/uRmU2GXfML/Sdwhrt/HLgO+E60zvXAV9392GICmNlw4DKgJZpcsA34es5TNkTLfwx8I1o2G9ji7odHe3z3F7EdANz9KTJ3910SjXNHsW+GiKSK8pskRu+4ByDvc/dtZvZx4HhgCnCLmX2TzC/mOOA+MwNoAl41syHAUHd/MNrEr8jMpNyTY4BDgYejbfUlMz1AVnZyx2XA30bft5CZ4yc7zjcsM3N0T9vJ9UNguru3FhibiKSU8pskiYqjwLh7B5kZoB8ws6eBc8n8Ij/bde/JMhMj5pv/ZSe7dwb7ZVcjM0/Q2XnWezf6t4P3Px/WTZxC28nVGX2JSANTfpOk0GG1gJjZWDP7WM6iCcDLwB+BvaMTGjGzPmZ2mLtvBraY2Sei538+Z93VwAQz62Vmo4Gjo+WPApPM7KPRtgaY2UEFhnYvmckvs+McVuZ2RKRBKb9Jkqg4CstA4AYzW2mZ2aAPBeZEx7HPAOZaZrbkFcBx0TpfAn5oZo8A23O29TDwZ+Bp4HvAcgB3Xw98EbgpivEocHCBcV0BDDOzZ6L4U8rcjog0LuU3SQzLnHAvaWBmY4C73H1c3GMREakm5TepJ3WORERERHKocyQiIiKSQ50jERERkRwqjkRERERyqDgSERERyaHiSERERCSHiiMRERGRHCqORERERHL8f9P5xKo/gZdCAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 648x216 with 4 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "import numpy as np\n",
    "import scipy.spatial\n",
    "import librosa\n",
    "import matplotlib.pyplot as plt\n",
    "from numba import jit\n",
    "%matplotlib inline\n",
    "\n",
    "import sys\n",
    "sys.path.append('..')\n",
    "import libfmp.c3\n",
    "\n",
    "@jit(nopython=True)\n",
    "def compute_accumulated_cost_matrix_21(C):\n",
    "    \"\"\"Compute the accumulated cost matrix given the cost matrix\n",
    "\n",
    "    Notebook: C3/C3S2_DTWvariants.ipynb\n",
    "\n",
    "    Args:\n",
    "        C (np.ndarray): Cost matrix\n",
    "\n",
    "    Returns:\n",
    "        D (np.ndarray): Accumulated cost matrix\n",
    "    \"\"\"\n",
    "    N = C.shape[0]\n",
    "    M = C.shape[1]\n",
    "    D = np.zeros((N + 2, M + 2))\n",
    "    D[:, 0:2] = np.inf\n",
    "    D[0:2, :] = np.inf\n",
    "    D[2, 2] = C[0, 0]\n",
    "\n",
    "    for n in range(N):\n",
    "        for m in range(M):\n",
    "            if n == 0 and m == 0:\n",
    "                continue\n",
    "            D[n+2, m+2] = C[n, m] + min(D[n-1+2, m-1+2], D[n-2+2, m-1+2], D[n-1+2, m-2+2])\n",
    "    D = D[2:, 2:]\n",
    "    return D\n",
    "\n",
    "@jit(nopython=True)\n",
    "def compute_optimal_warping_path_21(D):\n",
    "    \"\"\"Compute the warping path given an accumulated cost matrix\n",
    "\n",
    "    Notebook: C3/C3S2_DTWvariants.ipynb\n",
    "\n",
    "    Args:\n",
    "        D (np.ndarray): Accumulated cost matrix\n",
    "\n",
    "    Returns:\n",
    "        P (np.ndarray): Optimal warping path\n",
    "    \"\"\"\n",
    "    N = D.shape[0]\n",
    "    M = D.shape[1]\n",
    "    n = N - 1\n",
    "    m = M - 1\n",
    "    P = [(n, m)]\n",
    "    while n > 0 or m > 0:\n",
    "        if n == 0:\n",
    "            cell = (0, m - 1)\n",
    "        elif m == 0:\n",
    "            cell = (n - 1, 0)\n",
    "        else:\n",
    "            val = min(D[n-1, m-1], D[n-2, m-1], D[n-1, m-2])\n",
    "            if val == D[n-1, m-1]:\n",
    "                cell = (n-1, m-1)\n",
    "            elif val == D[n-2, m-1]:\n",
    "                cell = (n-2, m-1)\n",
    "            else:\n",
    "                cell = (n-1, m-2)\n",
    "        P.append(cell)\n",
    "        (n, m) = cell\n",
    "    P.reverse()\n",
    "    P = np.array(P)\n",
    "    return P\n",
    "\n",
    "# Sequences\n",
    "X = [1, 3, 9, 2, 1]\n",
    "Y = [2, 0, 0, 8, 7, 2]\n",
    "N, M = len(X), len(Y)\n",
    "\n",
    "C = libfmp.c3.compute_cost_matrix(X, Y, metric='euclidean')\n",
    "D = compute_accumulated_cost_matrix_21(C)\n",
    "P = compute_optimal_warping_path_21(D)  \n",
    "    \n",
    "plt.figure(figsize=(6, 2))\n",
    "plt.plot(X, c='k', label='X')\n",
    "plt.plot(Y, c='b', label='Y')\n",
    "plt.legend()\n",
    "plt.tight_layout()\n",
    "\n",
    "plt.figure(figsize=(9, 3))\n",
    "ax = plt.subplot(1, 2, 1)\n",
    "libfmp.c3.plot_matrix_with_points(C, P, linestyle='-', \n",
    "    ax=[ax], aspect='equal', clim=[0, np.max(C)],\n",
    "    title='$C$ with optimal warping path', xlabel='Sequence Y', ylabel='Sequence X');\n",
    "\n",
    "ax = plt.subplot(1, 2, 2)\n",
    "D_max = np.nanmax(D[D != np.inf])\n",
    "libfmp.c3.plot_matrix_with_points(D, P, linestyle='-', \n",
    "    ax=[ax], aspect='equal', clim=[0, D_max],\n",
    "    title='$D$ with optimal warping path', xlabel='Sequence Y', ylabel='Sequence X');\n",
    "for x, y in zip(*np.where(np.isinf(D))):\n",
    "    plt.text(y, x, '$\\infty$', horizontalalignment='center', verticalalignment='center')\n",
    "\n",
    "plt.tight_layout()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## LibROSA Implemetation\n",
    "\n",
    "The DTW function of [LibROSA](https://librosa.org/doc) allows for specifying arbitrary step size conditions. The following code cell calls the function `librosa.sequence.dtw` with different step size conditions $\\Sigma$.\n",
    "\n",
    "<font color=\"red\"><strong>Note:</strong> For general $\\Sigma$, there may not always be a valid warping path fulling both the boundary and the step-size condition. In this case, the function `librosa.sequence.dtw` outputs an **error message** (\"No valid sub-sequence warping path could be constructed with the given step sizes\").\n",
    "</font>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2024-02-15T08:50:53.301340Z",
     "iopub.status.busy": "2024-02-15T08:50:53.301091Z",
     "iopub.status.idle": "2024-02-15T08:50:53.917858Z",
     "shell.execute_reply": "2024-02-15T08:50:53.917289Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Accumulated cost matrix and optimal warping path for different step size conditions:\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAApUAAADQCAYAAABWW5pRAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAABBkElEQVR4nO3deXxU9b3/8dcnG2sQIigIVq5rFTUKuGsxQqlQ19alaFu1dblVe7XSW21/5jYu5WFbFb23tRWlaltFsct1DSgxSrFaG9xQckVRrCyyyJYoIdvn98c5E4cw+zkz58zM5/l45EHmzMz3fBLmnfnO95zz/YqqYowxxhhjjBclQRdgjDHGGGPyn3UqjTHGGGOMZ9apNMYYY4wxnlmn0hhjjDHGeGadSmOMMcYY45l1Ko0xxhhjjGcF16kUkfNF5Jkc7Wu0iKiItIrIpVlov4/bdoeI3Oxuq3Nvt4rIgCzs87tu2yoi+7rbVojINhH5g9/7c9t/TkTaRGRRNto34WaZ9bxPy6zJGcur530WdF49dSpF5HgR+buIbBGRjSLyoogc4d53YRB/cFT1QVWdnOPdDlbVWQAiUiEif3JfJCoiJyZ6ovui/p2IbBWRj0Xkmsh9qrpdVQcCD/Z62iOqOlBVP3XbEBH5uYh84n79QkQkwT5nicg7ItItIhdG36eqs9199naqqn4rqo2bRGSJiHSKSF2Sn/FgEZkvIhtEZKeJUVX1JODfE7Vh/GGZ7eEls+e4v8PPROT56PvSyOx/ishbItIiIh+IyH8m2WeuM3uBiCx2/y6tdP+mlEXt0zKbA5bXHkHn9WoRed/Nw2oRmRmdhxj7zHVev+Hub4uIrBORB0RkUNQ+c5bXjDuVbsFPAv8DVAEjgRuA7f6UltcWAd8EPk7hsXXAfsBeQA3wIxE5Oc39XQqcAVQDhwKnAJclePwbwOXAq2nuJ9p7wI+Ap1J4bAcwF/iuh/0ZjyyzCaWT2Y3AHcAtHvYnwLeBIcDJwJUi8o0Ej891ZvsDVwNDgaOAicAPPezbpMnymlCu8/oEMFZVBwEH47zX/keCx+c6ry8Cx6nqLsDeQBlws4d9Z8zLSOX+AKo6R1W7VHWbqj6jqm+KyIHAb4FjxBnm3Qw9o3K3isi/RGStiPxWRPq5953ofiL+iTuitUJEzo+3c/dT2vtRn/TPj9q+yP3+R+7+I18dInK/e98uIjJbRNaIyCoRuVlESt379hWRF9xe/wYReSTVX4qqtqvqHaq6COhK4SnfBm5S1U2q2gzcA1yY6v5cFwC3qepKVV0F3JaoDVX9tao2AG1p7ie6jQdUtR5oSeGx76jqbODtTPdnfGGZjSHdzKrqAlWdC6xOdR8x2viFqr6qqp2q+g7wGHBcgsfnOrO/UdW/ub+bVTgjOXHrM1lheY0hoLwuV9XN7k0BuoF9Ezw+13n9SFU3RG3qSlRfNnnpVC4DusQZZp0iIkMid7ido38HXnKHkAe7d/0cJyiH4fzAI4H/impzOM4n45E4HaVZInJA7x2Lc57DfwNTVLUSOBZ4vffj3D/cA92h5gOB9TgjZgAPAJ1uHYcDk4GL3ftuAp7BGUUYhfNJMbLvJ0XkuhR+P0m5v7M9cD7VRLwBjEmzqTE+tGEKn2U2hEREgBMI94euLxHu+gqR5TVEROQ8EdkKbMAZqbw74JJ2IM6pEltwOqFfxxmdzbmMO5WquhU4HlCc0bX1IvK4iOwe6/HuH85LgB+o6kZVbQFmAL0P+dS65zm8gDPse06cErqBg0Wkn6quUdW4f/DcT2r/C9ypqk+7NU4BrlbVT1V1HTAzqpYOnMPRe6hqm/uJKPJzn6KqXobRo0XOq9gStW0LUJlBO73bGOj+zo0BLLNxfzHBq8P5W3xfwHXEJCIXAeOBW4OupZhYXsNFVR9yD3/vjzNKvDbgknagqovcw9+jgF8CK4Kow9OFOqrarKoXquoonPMM9iB+73gYznk6i0VksztcP8/dHrFJ3RNjXR+6bfbe76fAuTif1NaIyFMi8sUEpc4G3lHVn7u39wLK3edGarkb2M29/0c4Q9yviMjbIvKdBG170er+Oyhq2yBSGO6O0U7vNlpVdaeLYkxxs8yGi4hciXMKzFdVNXTnyonIGTjnok3pdXjN5IDlNXxU9V2cUfu7gq4lFvd0lXnAw0Hs37cphVT1/4D7cV744Hy6irYB2AaMUdXB7tcuuuNVUENkx0v4v0Cc8yBUdb6qfhkYAfwfzie5nbjD6Aew40UiH+Gc7Dw0qpZBqjrGbftjVb1EVffAueDlLnEv/feTqm4C1uAMpUdUk/5hprd9aMMUGctssNw30uuAiaq6Muh6ehPngsF7cK5KXRJ0PcXO8hoqZcA+QReRQGD1ebn6+4siMl1ERrm39wSmAS+7D1kLjBKRCgBV7cZ5Uc4Ukd3c54wUka/0avoGcaYMOAHnKuZHY+x7dxE5zQ3HdpyRup1O2BWRKThXaJ2hqtsi21V1Dc75HLeJyCARKRGRfURkgvu8syM/F7AJJ7ypXHQT2W8fEenr3qwQkb4JDkX/HrheRIa4nwQvwfnDkY7fA9e4v889gOmJ2nB/v31xPimWu/Wl9VoQkXK3jRKgzG2jNM5jxX1shXu7r4j0SWd/xjvLbMLfTcqZFZFS97FlQIn72PJU9+W2cT7Oockvq+r7KTw+15k9CefinK+r6ivp7Mf4w/Ka8HeT67xeHPU7PQj4MdCQ4PG5zuv5IvIF9712L+BnierLJi8jlS04U038Q0Q+xXmhv4XToQF4Dme07GMRiRw2uRbnMvmXxTnhdQHOJ5yIj3FeYKtx/qD9u/vpLFbd093HbQQm4Fy+39u5OEP/zfL51Wm/de/7Nk4nZ6m7zz/hfCIDOML9uVqBx4GrVPUDABGpF5GfJPndvIPziXEkMN/9fq84j/0psBznMMQLwC9VdV6S9nu7G2fKgyU4/wdPkfgk4mfcmo4FZrnffynNfd7jPm8a8P/c778V57F7ufdHRk+34fyOTG5ZZuNLJ7Pfcu//Dc4FNtuIM4qTwM3ArsA/Y/ycseQ6s7XALsDTUfXVp7k/443lNb5c5/U4YIn7//C0+5Woxlzn9SDg7zid/xdxfj+XpLk/f6hqKL6AE4GVQdeRZs174UwZsBm4JAvt93Hb/hT4qbvtevf2ZmBAFvZ5kdt2G7C3u+0dYCvwQJZ+j8/i/AFtCPr/1L7S+n+zzO7cvmXWvkL5ZXmN2b7l1ecvcXcYOHFmxf+jOickG2NCzjJrTP6wvJpcKLi1v40xxgRLRPYUkUYRaRbn6t6r3O1VIvKsiLzr/jskWVvGmOzzK7OhGak0xhhTGERkBDBCVV8VkUpgMc5SshcCG1X1FnGuGh6iqtcGV6kxBvzLrI1UGmOM8ZU6k2W/6n7fAjTjXFRxOs5KK7j/nhFIgcaYHfiV2VCNVJaWlmpZWZnndrq6Up6ZIGdtjRrl32ksW7du9aWdwYMH+9LOJ598Qmtra9LVe0Qk0Yttvqqe7EtBJieGDh2qo0ePDrqM0Hvttdd8acevvAJs3Lhxg6oOS/SYJHkF58rj6LWNZ6nqrBjtjAYW4syv+C/9fElBRGSTqtoh8Bwo5Lx+9tlnvrXVv39/X9ppaUl3DZP4li1bljSvkDSzKeXVbWc0GWbWew/OR2VlZb50vjZu3OhDNY7Nmzf70s4PfvADX9oBmD9/vi/tnHHGGb60c8stqa+oFWcqMVR1qC/FmJwZPXo0TU1NQZcRepWV6a66GttXvtJ7usHMzZkz58NUHhcvrwDqLK83PsnzBwJ/xlmub2ui9kx2FXJeFy9e7Ftb48aN86Wd559/3pd2AGpqalLKKyR8j02aV/f5njJrh79NTolIzC9jTPjEy2sqmXUnmP4z8KCq/sXdvNY9dytyDte6rBVvTBHy8h7rR2atU2lyqqSkJOaXMSZ84uU1WWbd1U1mA82qenvUXY8DF7jfXwA8lpXCjSlSmb7H+pXZUB3+NoVNRKwDaUye8JjX43BW/1giIq+7234C3ALMFZHvAv8CzvZapzHGEYbMWqfS5JQd6jYmf2SaV1VdhLPucSwTMy7IGJNQ0Jm1TqXJKRupNCZ/WF6NyS9BZ9Y6lSamfV95haMff5yBGzfSWlXFy6edxntHHumpTTv8bUx2HLtiBee+8Qa7fvYZn/TvzyPV1fzd4/Qxlldjsme3BQvY+9576bNuHdt32433L76YdZMmeWozDJm1TqXZyb6vvMKJDz1EeXs7AJUbN3LiQw8BeO5YBv2CN6bQHLtiBZe88gp93Dl1h332GZe88gqA546l5dUY/+22YAEH3Horpdu3A9B37VoOuPVWAM8dy6Aza38xzE6Ofvzxng5lRHl7O0c//rjntm1KoeJQX19PTU0NtbW1NDY2MmHCBKZPn5737fjdlh/OfeONng5lRJ+uLs594w3PbXuZUsjkD8tZbmva+957ezqUEaXbt7P3vfd6rjHovFqn0uxkYJzJ4+NtT1VkaN6mFCp8DQ0NNDY2UllZycyZM1mwYAHV1dUsXbo0r9vxuy0vRJUjPvqIoXFWE9nV4yojifJqmS0slrPc1tRnXeypHuNtT1UY3mPt8LfZSWtVFZUxOpCtVVWe27Y3o+IQWf61qqqK1tZWurq6UFXSXRY2bO343VYmRJXxK1fytbfeYq/Nm+kUoSzGvj/xYbk5y2txsJzltqbtu+1G37VrY273KujM2l8Ms5OXTzuNjoqKHbZ1VFTw8mmneWo3DJ+iTG5MmjSJiRMn0tzcTF1dHZMnT6apqYkxY8bkdTt+t5WOyMjkjHnz+MGiRZR3dfHro4/m7qOOYntp6Q6P3V5ayiPV1d72ZyOVRcNyltuaVpx/Pr27oV19+vD+xRd7qi8M77ESZK+/tz59+mihrv192223+dIO5Gbt731feYWTfv97Srq7k179fcstt/Dhhx8mPWmjvLxchwyJvQ79+vXrF6eyLqkJj/Hjx2uhriXsJ69rf4sqp3Z1MaNvX/bavJnVlZX8dcwYXtprL7rdN4t0r/6eM2dO0rwlyitYZvNNIec139b+Hj5vHl/8+c/ZXlVFxaZNSa/+rqmpSSlrYXiPtcPfJqb3jjyScfPmsXGPPXjW46enaDbCYUxqIp3J69rbOaS7m9Xl5fz66KN36ExG/H30aM9XesdieTXGf8Oef55tw4fzj4ceAp8vogk6s9apNDkThjm0jAm73p3Jd0W4uE8f2qZO3akzmdU6LK/G+K6spYUhixez8utf971DGYbMWqfS5FTQL3hjwipeZ/JPZWV0iXBqANmxvBrjr6EvvkhJZyfrTzwxK+0HnVnrVJqcsvntjNlRss5koLVZXo3xVeTQd8sBB2Sl/aAzax9DTc54uTJNRPYUkUYRaRaRt0XkKnd7nYisEpHX3a+pWf9BjPGBqHJaZycvbtvGg21t9FXl4j59OKJ/fx4pLw9Fh9Ku/jbGP5FD3+snTPD90DeE4+rvUI1UVlZWMmHCBM/t+Hn1t59t+cWvKRUSXdkJUFpaSkV5eUqPS5WHF3cnMF1VXxWRSmCxiDzr3jdTVW/NtGFjevNz5Y7W1tYdbgtwBvBToBp4B/gm8LAqXdu3Q6+VNiLOO+8832qaM2dOSo+zzqPJpssuu8yXdu6++25f2gEYMWKEL+2sWbNm54333w+dnXxh+nS+cMQRvuynt6AzG6pOpSl8mb7gVXUNsMb9vkVEmoGRPpZmTFbF7UwCXfGfFqig36CMKShz58Lo0TA+ezP7BJ1Z+4thcibJOsJDRaQp6uvSBO2MBg4H/uFuulJE3hSR34lI4mFVY3JMgDOB14C/AH1xOpNjgAcJb4cyUV6DPm/LmLyzaRMsWABnn52VQ9+Q9D02J2yk0uRUgk9RG1KZmFVEBgJ/Bq5W1a0i8hvgJkDdf28DvuNTucZkTFQ5+L33eI38GZnsLehRD2MKxmOPQUeH06nMoqAza51Kk1NeXvAiUo7ToXxQVf8CoKpro+6/B3jSa43GeBHpTE5+6SX22LAhLzuTEUG/QRlTMHJw6BuCz6x1Kk3OeBmGF+eJs4FmVb09avsI93xLcI4yvuW5UGMy0LszuW7IEB48+WQumDcv7zqT4C2vxpgokUPfV1+dtUPfEI7MZrVTKSKDgXuBg3EOT35HVV/K5j6NP76waBGVH3/MoFWrOPX73+eNc8/lX8cf77ldD5+ijgO+BSwRkdfdbT8BponIYTivrxWAP5cTFiHLa2oOb25m6qJFDGlpYVNlJfXHHUdHeflOncnXv/hFuktK6Jo3L+iSMxb0qIdJzDKbJ3J06BuCz2y2RyrvBOap6lkiUgH0z/L+jA++sGgRR95zDyVdzvjKgA0bOPKeewA8dyw9XP29COeah96e9lSQiWZ5TeLw5mbOefZZKjo7AahqaWHavHmUwE6dyUIQ9BuUScoymw9ydOgbgs9s1jqVIjII+BJwIYCqtgPt2dqf8U/1I49Q1r7jf1VZezvVjzziqVMZhnVJTWyW19RMXbSop0MZUQK09u3LLy+4oGA6k2B5DTvLbJ7I0aFvCEdms7n3vYH1wH0i8pqI3CsiA3o/SEQujUwj09bWlsVyTKr6b9gQe/snn3huO+jpDkxcaed1/fr1cRurr6+npqaG2tpaGhsbmTBhQkYTioetnSEtLTG3D2hrK6gOZUQYpigxcSXNbL7m1U+B1xTj0Hc2awo6r9n8K1gGjAV+o6qHA58C1/V+kKrOUtXxqjq+b9++WSzHpKS7m84+fWLe9dmuu3puPuglpExcaed12LBhcRtraGigsbGRyspKZs6cyYIFC6iurmbp0qVpFRW2djZVVqa1Pd/ZMo2hljSz+ZpXPwVeU4xD39msKei8ZnNPK4GVqhqZoPpPOAEwYdXdzRH33kv59u1091p6sbOigjfOPddT82FYl9TE5WteVRWAqqoqWltb6erqQlV7tudrO08ff/xOa3K3l5XxtA8XsYWNrf0der5lNmw581OgNcWZ8DxbNYXhPTZre1LVj4GPROQAd9NEILiPKyYxt0O5T2Mjb33ta7x82WV8OnQoKsKnQ4fyyiWX+Hb1t71BhY/feZ00aRITJ06kubmZuro6Jk+eTFNTU9rr1oetnQ9HjEBUaSsvR4GNlZXM/fKXee3AA9NqJ19YpzK8/Mxs2HLmp0BrinPVdzZrCjqv2b76+/vAg+5Vae8DF2V5fyYTvTqUb511FojwrxNO8H1Xdi5WqPmW1ylTpjBlypSe2wsXLiyIdia98gpdpaX8/MIL2Vqgh7yjWV5Dz5fMhi1nfgq0pjhXfWezpqAzm9VOpaq+DmT/GnqTuTgdymwIw5VpJj7La2JVmzczfulS/n7ooUXTobS8hptlNsRyeNV3RBgya38xilkOO5QRQQ/NG5OpSa+8QrcIzx1xRNCl5IyXw98i8jsRWScib0VtqxORVSLyuvs1Nas/gDFByeGE59GCzqu9mxcpUc15hzIMJxEbk4nIKOXLhxxSFKOU4MuFOvcDJ8fYPlNVD3O/bPECU5hyOOF5hMf32PvxIa+29ncRElXOXrCAfd56K2cdyp592zlaJg8V4ygleMurqi4UkdH+VWNMftiluzvnh74jMs2sX3kNVady4MCBnODDxSGbNm3yoRp/vffee761dcYZZ2T+5O5uDrrzTkYF0KGE4JeQMiYVra2tPd8P3bqV8W+/zQsHHcRqEYi6LxXTpk3zpabVq1f70k46kuR1qIg0Rd2epaqzUmj2ShH5NtAETFfV8P3BNjlx9913+9JOXV2dL+0ArFmzxnsj998PF13E/EGD2Dhnjvf20pAgsznJa6g6lSbLIh3KefNYft55vHXKKTntUIbhJGJj0jXltdfoLilh/mGHBV1KTqWQ1w2qmu6xvd8ANwHq/nsb8J3MKjQmpObOpXXYMDbuvXdOd5sksznJq73DF4teHcrl3/52zoflwS7UMfll6NatHLNsGX/74hfZPGCnVSsLnt/zVKrqWlXtUtVu4B7gSF8LNiZo7lXfHx11VN6/x2aSV3s3LwYh6VBC8OuSGpOOYh2ljPB77W8RGRF180zgrXiPNSYvuVd9/+uoowLZfdB5tcPfhS5kHUoblTT5IjJK+cJBBxXlKKXXvIrIHOBEnHO5VgI/BU4UkcNwDqetAC7zXKgxYeJe9Z3rQ9/gLbN+5dU6lYUsRB3KCOtUmnxR7KOU4C2vqhrrCqXZmVdjTMgFMOF5bx4OdfuSV+tUFpDhzz3HfvfdR9/162kbNozPhg9n1zffDE2HEjKf7kBE9gR+DwwHunGuXLtTRKqAR4DROJ+kzrGrSY1XxT5KGWGnphiThugJz32c8SUdQWfWho0KxPDnnmPMHXfQb906RJV+69ax65tvsvaYY0LVofRwEnEnznQGBwJHA1eIyEHAdUCDqu4HNLi3jfHERil9mfzcmOISwITn0cKwwIj9ZSgQ+913H6Xbt++0fdDy5aHoUEZk+oJX1TWq+qr7fQvQDIwETgcecB/2AHBGdio3ReP994v6iu9o1qk0JkWRQ99nnx3oe27QebW/DAWi7/r1aW0PSoIX/FARaYr6ujReG+6s/4cD/wB2V9U14HQ8gd1y8GOYJOrr66mpqaG2tpbGxkYmTJjA9OnT86OdGTPoLilhXhGPUkZYp7I4hC2vfrfll4Q1BbTWd29B5zXunkRkuoiUxti+q4jYydYh0zZsWFrbg5BkepINqjo+6ivmTP8iMhD4M3C1qm7NZf1hF6bMNjQ00NjYSGVlJTNnzmTBggVUV1ezdOnScLfz/vvwwAMsPPBAthT5KGWivAZ93lYhsLzmri2/JKwp4EPf4P8UYJlI1H09AFgsIsdFNojI5ThL9SzJdmEmPe9edBFdffrssK2rTx/eveiigCqKzcunKBEpx+lQPqiqf3E3r43MpeX+uy4rheeH0GRWVQGoqqqitbWVrq4uVLVne2jbmTEDSkuZX12dVvuFykYqs8rymqO2/BK3ps2bQ3HoG0I8UqmqlwKXA/8jIn8QkX8CxwPHquodOarPpOjjk07i7auvZvvgwQBsHzyYt6++mo9POinYwqJ4OYlYnI9as4FmVb096q7HgQvc7y8AHvO98DwRpsxOmjSJiRMn0tzcTF1dHZMnT6apqYkxY8aEtx13lJLLLiv6UUqwC3WyzfKau7b8EremZctCceg7DBfqJJtS6C3gn8DJgOBcfevDausmGz4+6STahwxh/HXX8cb117P5kEOCLmknHl7cxwHfApaIyOvutp8AtwBzReS7wL+AYFMdvFBkdsqUKUyZMqXn9sKFC8PfjjtKybXXwg03ZLSfQmOdx6yzvOagLb/ErWnq1MAPfUcEndlE51R+E3gdeB/YB2eJnl+IyO9FxC6GMBnJ9HwPVV2kqqKqh6rqYe7X06r6iapOVNX93H835uDHCCXLrAdRo5TssUfQ1YRGGM7RKlSW1wIRkqu+I4LOa6KRyrOBGlX90L29WESOAf4deBnI/RpEJq/ZMo1ZZ5nNVPQopQEsrzlgeS0EIbnqG8KR2bidSlU9PcY2BX4jIn/KalWmYAX9gi9kltkMRUYpL7/cRil7sbxmj+W1QITgqu9oQWc2o2UaVTVckx+avBH0C75YWWYTsFHKuCyvwbC85okQrPXdW9CZtbW/Tc7YuVgmdGyUMi7LqzFJhOjQN4Qjs6HqVJaUlNCvXz/P7fjRRsS2bdt8aWfTpk2+tAPOHFnxVFZWAjBo0CBKEjwOYPny5b7U09XVlfJjg/4UZQpbuitunPPMM4xVZUZbG1ujnrt69Wrfaqqrq/OlnYceesiXdtJheTXZtHjxYl/a8StjAMuWLUv5sSPvu4+KkSP5YNAg6PW8adOm+VbTeeedl/Jjg85s0r2LSH8RqRWRe9zb+4nIKdkvzRSioOfQKgaW2dRUbd7M+KVLeenQQ9k6cGDQ5YSSzVOZfZbX/FSyZQv9X3qJlpNPDs2hbwj+PTaVPd0HbAeOcW+vBG7OWkWmYIVhYtYiYZlNwaRXXqFbhMYjjgi6lFCyyc9zxvKahwY2NCAdHU6nMiTC8B6byp72UdVfAB0AqroNZ5JWY9IW9BxaRcIym4SNUqbG5qnMCctrHqqsr6dj5Ei2h2yRkaDzmso5le0i0g9QABHZB+dTlTFpsxGOnLDMJmGjlKmxvOaE5TXPRA59b/r2t0N16BuCz2wqncqfAvOAPUXkQZzl8i7MZlGmMIVhYtYiYZlNIDJK+WJ1tY1SJmB5zRnLa54J46FvCEdmk3YqVfVZEXkVOBpnSP4qVd2Q6g5EpBRoAlapqp18nEVV9fXsOXMmAPtcdx0fXX01G6PWKQ0DO2yWfV4yW8h5Pby5mamLFjGkpQWA9YMHB1tQHrC8Zp/lNX9UPvEEQ2+/nbLVq9HSUipWrGD7oYcGXdYOgs5sKld/nwl0qupTqvok0CkiZ6Sxj6uA5gzrMymqqq9nrxkzKHenLirfuJG9Zsygqr4+4Mp2FPRJxMXAY2YLMq+HNzdzzrPPUtXSguC8c5/6t79xeHPB/ai+sgt1ss/ymh8qn3iC3a+/nvLVq52/IV1d7F5bS+UTTwRd2g6Czmsqe/qpqm6J3FDVzTjD9UmJyCjgq8C9GVVnUjbyrrsobWvbYVtpWxsj77oroIp2FoYr04pERpkt5LxOXbSIis7OHbZVdHYyddGigCoKP7v6O2csr3lg6G23UdLrPbakrY2ht98eUEU7C8N7bCp7ivWYVCdNvwP4EdAd7wEicqmINIlIU4t7WMqkr2Lt2rS2ByXoF3yRyDSzd5BGXtevj7+SXH19PTU1NdTW1tLY2MiECRPSnpjcz3aGxPnbEm+7cVinMicsr1nkpSbZvp0BCxYwfPp0ytasifmYeNuzVVMyQec1lT01icjtIrKPiOwtIjOBpNPgu5O3rlPVhI9V1VmqOl5Vx0dWgzHpa99997S2ByXo6Q6KRNqZzSSvw4YNi/u4hoYGGhsbqaysZObMmSxYsIDq6mqWLl2a1g/iVzub4vxtibfdOGxKoZywvGZRujVFdyT3PuYYRl5xBQMWLUL794/5+M4RI7JeUzqCzmsqncrvA+3AI8CjQBtwRQrPOw44TURWAA8DJ4nIHzOs0ySx6vLL6erbd4dtXX37suryywOqaGdhGJovEplk1te8qirgLCna2tpKV1cXqtqzPdftvDB2LL2f0V5WxtPHH59WO8XEDn/njOU1i1Kqqa0N/vd/4fzzd+hItk6dysrZs1m+aBFrb7yR7l7vsd19+7LhmmuyU1MGwvAem8rV358C16XbsKr+GPgxgIicCPxQVb+ZbjsmNZGrvPecOZPyTZvoqKoK5dXfmb64ReR3QOTT+cHutjrgEiByXOcnqvq0D2XmtUwy63deJ02axMSJEznssMOoq6tj8uTJVFdXc8EFFwTSzh6ffEKXCK0DBrBLayubKit5+vjjee3AA9Nqp9hY5zH7LK/ZFbemc891OpKPPgpPPAEtLbDrrrROnUrLySfz2VFHQXl5Tzstp54K4Fz9vWYNnSNGsOGaa3q2+1KTD7+noDObtFMpIvsDPwRGRz9eVU/KXlkmExunTKFj11054IorWH7LLbQefnjQJe3A4xxa9wO/An7fa/tMVb3VS12FJgyZnTJlClOiPtAsXLgwsHZ65qU87DAeq6nJqI5iFIY574qB5TW7dqiprY2F11zjdCR3262nI8m558LZZ0NNDWs/+CBuWy2nnppRJzJhTfj3ewpDZlM5GfhR4Lc4V5h1ZbITVX0eeD6T55rCkum5Haq6UERG+1tNwfKU2ULLq62ekzk7dzInLK/Z1NYG8+btNCIZ3ZGMHpHMd0FnNpVOZaeq/ibrlZiikIVPUVeKyLdxJgCerqqb/N5BHrLMumz1HG+85DXOKStVOOcOjgZWAOdYZi2vviuyjmQ0n08xSzuvqez9CRG5XERGiEhV5Cujqk1RS3IS8dDI1Bfu16UpNPkbYB/gMGANcFv2qs8rllmXjVJmzocLde4Heq9jdx3QoKr7AQ1kcL5+AbK8+iHqYht22w3OPBPmz3c6kvPnw5o1cM89MHlywXYoPV6ocz8+5DWVkcrImaP/GbVNgb1TeK4xO0jw4t6gquPTaUtVeybhFJF7gCc9lFZILLPYKKUfvIxUxjll5XTgRPf7B3AO216b8U4Kg+U1Q7J9O4Neegluv73oRiTjyTSzfuU1lau//y398jJTUlJCv379PLezceNGH6pxbNrkz5GZL33pS76046e+vaZHyFQ653D4eb6HiIxQ1cjMs2cCb/nWeB7LZWbDpibqQpxD//u/obSU7VddRc2uuwZWk1+LOuyxxx6+tJOOJHkdKiJNUbdnqeqsJE3uHsmsqq4Rkd281pjv8i2vy5Yt862tcePGpf+kGIe2NwB/xTk5tfGTT+i89164N7OFhvya/uiAAw7wpZ10JchsTvKaytXf/YFrgC+o6qUish9wgLtGqTEp83JlmojMwfnENFREVuIsY3aiiByG86l+BXCZL4XmOcss9P/4Y0Y99xwfTp3K9gA7lPkshbymfXTB7MzymoJY50hWVfWMSI74ylfoTN5KwUuS2ZzkNZXD3/fhzO5/rHt7Jc4HAnvBm7R5GJqfFmPzbG/VFKyiz+y+c+eiJSW89/WvB11KXsvChXVrI0cYRGQEsM7vHeShos9rTEk6ktGHtq1D+TmfM5t2XlPpVO6jqueKyDQAVd0mQV+zbvJW0HNoFYmizqyNUvonC3l9HOccwlvcfx/zewd5qKjzuoM0OpImNp8zm3ZeU+lUtotIP5xDjIjIPsB2D0WaImVrBudMUWfWRin94TWvcU5ZuQWYKyLfBf4FnO1DqfmuqPNqHUn/eMmsX3lNpVP5U2AesKeIPIiz5uiFGVVtip6NVOZE0WbWRin95fHq71inrABMzLjRwlR0eZXt2+n/t7/BDTdYR9JnPp9iBmnmNZWrv58VkVeBowEBrlLVDensxORGVX09e86cCcA+111XUGt/m9QVc2ZtlNJfltfsK7S8Vj7xRMz1sSMdycr6egY0NlL66afWkcyCoDObytXfkblwIvNiHCQiqGrwi3qaHlX19ew1YwalbW0AlG/cyF4zZgCEpmNph79zo2gz+/77NkrpI8trbhRSXiufeILdr7+eksj70OrV7P7jH7PLH/9In3ffpfTTT+kaPJjWqVNpOflkRn3rW9aR9FEYMpvK4e/oCVn7AkfiXKmWs8XuTXIj77qrp0MZUdrWxsi77gpNpxKC/xRVJIozszNm2CilzyyvOVEweR16++09HcqIko4O+r3xBlvPOouWk0/ms6OO+rwjaR1K3wWd2aR7V9VTo76+DBwMrE32PJNbFWtj/5fE2x4UD0tImRSFIbP19fXU1NRQW1tLY2MjEyZMYPr06dlr5/334YEH+NdXvmKjlD7yuEyjSUEh5bVszZq49629+WY+O/74lDuSftXkpzDW1FvQec1kTytxXvQmRNp3iz3Rffvuu+e4kvg8rktqMpfzzDY0NNDY2EhlZSUzZ85kwYIFVFdXs3Tp0uy0M2MGlJay/KyzfPwpipsPa3+bzORtXjtHjEhrey5q8lMYa4oWhvfYVM6p/B/cqQ5wOqGHAW9ksSaTJmlro2vAABTnLO+Irr59WXX55UGVFVPQ53sUgzBkNrLUWVVVFa2trXR1daGqaS+BllI77igl3/sebTZK6SvLa/YVUl43XHMNw6+7Dun8fDry7r592XDNNYHV5Kcw1tRb0JlNpfvahHN+x2LgJeBaVf1mVqsyKZO2NvabPp1+H3zA+jPPZPvw4agI24cP58Of/CRU51NC8EPzRSLwzE6aNImJEyfS3NxMXV0dkydPpqmpiTFjxvjfjjtKybXX+vxTGBupzImCyWvLKafQNXgw3RUVqAgde+zB2ptvpuXUUwOryU9hrKm3oPOaypRCD+SiEJO+SIey8p//ZEVtLZ9kENxc8rL2t0ldGDI7ZcoUpkR9oFm4MLMLWZO2EzVKyciR8NprGe3H7MzymhuFlNeKZcso27CBtXV1bJkWb9rD3NbkpzDWFC0MmU3l8PcSPh+a3+EuQFX1UN+rMknlW4cyIugXfDEoqszaKGVWWV6zr5DyWllfj5aU0Dp5ctClFK2gM5vKlEL17r9/cP89H/gMCPzTVbHK1w4lBH++R5Eojsz2HqU0vrO85kRh5FWVynnz2HbEEXTZuc2BCTqzqXQqj1PV46JuXyciL6rqjdkqysRX2t6e1x3KoD9FFYniyKyNUmaV5TVnCiKvFcuWUfHBB2y64IKgSylaYchsKnsfICLHR26IyLHAgOyVZOIpbW9nym9/m5cdyoigTyIuEoWf2cgo5aWX2ihlFtmFOjlREHm1Q9/hEHReUxmp/C7wOxHZBee8jy3Ad7JaldlJpEM56p138rZDGYZPUUWi8DNro5RZZ3nNmfzPqx36DoUwZDaVq78XA9UiMggQVd2S/bIKzwknnJDxc6Wtjb2+/30GvPMOq268kb/27w8huuqsq6sr5ccGfb5HMci3zNbV1aX1+CGbNnHlfffRNH489ffcs8N9y5Yt86UmP1fJmDVrli/t3H333b60A/DDH/4wpcdZXrMv3/K6//7777xxyRL44AMqfvSj2PfHcd555/lSk5/zRF522WW+tPPOO+/40g6kl8OgM5u0Sysiu4vIbOARVd0iIgeJyHdzUJshqkP5j3+w6sYb2XzGGUGX5EmmQ/Mi8jsRWScib0VtqxKRZ0XkXfffIVktPk8UemZP+NvfUBEWHXdc8gcbT+zwd/YVRF7nzoWSEvja14KupOgFnddU9nQ/MB/Yw729DLg6S/WYKIXWofS4hNT9wMm9tl0HNKjqfkCDe9sUcGaHbNpE9RtvsHjcOFoGDQq6nIJmyzTmzP3kc15V4dFHYcIEiLNcsMmNMCzTmMqehqrqXKAbQFU7gdSPd5qMFFqHMiLTF7yqLgQ29tp8Op9Pu/EAcIavxeavgs2sjVLmlnUqcyK/8/rWW/DOO3DOOUFXYgh+pDKVC3U+FZFdcSdnFZGjcU4kNllSqB1KSHi+x1ARaYq6PUtVk52MtruqrgFQ1TUiYh+THQWZ2cgoZdP48TZKmSNBn59VJPI7r3boO1SCzmwqncprgMeBfUTkRWAYcFZWqypihd6hTPCJaYOqjs9lPQWsIDNro5S5FYYrSYtE/ubVDn2HShgym8rV36+KyATgAJxlo95R1Y5kzxORPYHfA8NxhvVnqeqdHustSLs89RS733kn5R9/jJaXI+3trLrppoLqUEb4/IJfKyIj3FHKEcA6PxvPV5lkNsx5PWTJEr787LNUtrTQXlHB6A8/ZMkhhwRdVlEI+g2qGOT1e2zk0PfVV+d0tya+oDMbd+8icoSIDIeeczzGAT8DbhORqhTa7gSmq+qBwNHAFSJykA81F5RdnnqKkXV1VKxZg6hS0t6Olpej5eVBl5YVIhLzK0OPA5HlGy4AHvOlyDzlMbOhzOshS5Zw6hNPMKilBQH6tLdz6hNPcMiSJUGXVhTi5TXoQ2yFoCDeY+3Qd+gEnddEXdq7gXYAEfkScAvOp6ItQNKJ11R1jaq+6n7fAjQDtvRFL7vfeSclbW07bCvp6GD3O0MxSOQrL1emicgc4CXgABFZ6U65cQvwZRF5F/iye7uYZZxZv/NaX19PTU0NtbW1NDY2MmHChIzmfpzY0EBFx46DNhUdHUxsaMi0tILk1+87ml39nXWheY/N6PUT49B3Nl6HXhVTTWG/+rtUVSNX256LM7T+Z1WtBfZNZyciMho4HPhHjPsuFZEmEWnaunVrOs0WhPKPP05re77zcPX3NFUdoarlqjpKVWer6ieqOlFV93P/7X11eLHxJbOp5nX9+vVx22hoaKCxsZHKykpmzpzJggULqK6uZunSpWn9QLtsiX29Qrztxcqv33dv1qnMqqy/x2Y1rzGu+s7W69CLYqsp6Lwm7FSKSOScy4nAc1H3pXKBDwAiMhD4M3C1qu7Ua1TVWao6XlXHDyqyKzqlrS3uYe6O4cNzXE32heFTVIHznNl08jps2LC47URWuKiqqqK1tZWuri5UNe2VL1oqK2Nu37LLLmm1U+j8+n1Hs5HKrMv6e2xW8xrj0Hc2XodeFVNNYXiPTbSnOcALIvIYsA34G4CI7EuK0x2ISDnOi/1BVf2Lx1oLSuQqb2lvp7tXx7K7b1/WXnVVQJVlV9DnexQ4T5n1M6+TJk1i4sSJNDc3U1dXx+TJk2lqamLMmDFptfNJVRW9/8y2l5fTMHGil/IKjl+/797snMqsCs17bNqvnzhXfWfrdehFsdUUdF7jfhpS1Z+JSAMwAnhGP+9ClwDfT9awOD/FbKBZVW/3o9hCscO0QTffjJaV9Vz93TF8OGuvuootX/1q0GVmhY1wZI+XzPqd1ylTpjBlypSe2wszWKt+yKZNfOGjj1i+994M/eQTdtmyhS277ELDxIl29Xcvfvy+Y7G8Zk+Y3mPTfv3Eueo7W69DL4qtpqAzm3CIXVVfjrFtWYptHwd8C1giIq+7236iqk+nVWGB2aFDedNNbD79dICC7URGC8McWoXOQ2ZDl9fIvJSPnX66TXYeAMtr9uXte6xd9R1KYchsyudtpEtVF+HMuWVc8TqUxSToF7yJLWx5tdVzwsHyGl6BZdYmPA+1oDObtU6l2ZF1KB12LpZJha2eEw5e8ioiK4AWnHWsO23FrAJhE56HWtCZtU5lDpR1dFiHknAMzZvws1HKcPAprzWqusGPekxI2KHv0ApDZq1TmWVlHR2cPns2A957r6g7lBHWqTTJ2ChleFhezQ7s0HfoBZ3ZUHUqKysr+dKXvuS5nY0b/ZsDe9u2bRk/t+eQ93vv8dr3v89HBx8My5d7rmnIkCGe2wBY7kMtAJ2dnSk/NugXvAmfU089tef7ipUrOfjmm1l/1lmceP75gdU0btw439q6++67fWsr15LkdaiINEXdnqWq0SvBKPCMiChwd6/7TD5yD31/eOaZbFi82HNzDz30kA9FwZw5c3xpB/I7r5Aws8nyCj5kNlSdykLS+xzKjw4+OOiSAmfz25lkRtx3H1payscXXJD8wSarUsjrhiTnXB2nqqtFZDfgWRH5P1UNfj4Xk7lHH0VLSth80klBV2JiSJLZZHkFHzJrw0ZZYBflxBf0bP8mvCpWrmTXJ59kw5ln0mGH1kLBy4o6qrra/Xcd8FfgyCyXa7JJFebOpWXsWDqrqoKuxsTh5T3Wj8zau7nPrEOZmHUqTTw2Shk+mXYqRWSAiFRGvgcmA2/loGSTLe6h702TJgVdiUkg0/dYvzJrh799ZB3KxOzqbxNPZJRy/Vln2ShlSHjM6+7AX91DcWXAQ6o6z6/aTAAefRTs0HeohSGz1qn0iXUoU2PnVJpYbJQynDLNq6q+D1T7W40JjHvomwkT7NB3yAWdWetUerDLU0/1rNmt5eVIezurbr7ZOpQJ2Eil6c1GKcPL8moAm/A8jwSdWetUZmiXp55iZF0dJW1tAEh7O93l5WiZ/UrjscPfJhYbpQwny6vp4R765mtfg48+CroaE0cYMmt/MTK0+5139nQoI0o6Otj9zjsDqig/RKY86P2V4nNXiMgSEXm913xbJl+9/75d8R1i8fJqp7EUkahD3zbhefgFnVfrVGao/OOP09puHD5c/V2jqofZOsLhVl9fT01NDbW1tTQ2NjJhwgSmT5++8wNnzMi7UcqUf7YcteN3W9G8TClk8kfC10/k0Pc554SnpoCEMfu9BZ1X+8uQAWlrQ8vLY97XMXx4jqvJH5GheXuDKnwNDQ00NjZSWVnJzJkzWbBgAdXV1SxduvTzB73/PjzwQN6NUqb0s+WwHb/bikiUV8tsYUn4+ok+9B2WmgISxuxHC8N7rJ0AmKbIVd6RcyhLOjp67uvu25e1V10VYHXhF/QSUiY3VBWAqqoqWltb6erqQlV7tgMwYwbk2SglpPiz5bAdv9uKZp3H4hD39dPdHdih72y9psNQUzZ/tqAza38x0rDDtEE338yqm26ifcQIVIT2ESNYVVfHlq9+NegyQy3B+R4bVHV81FesDuNxqjoWmAJcISLeF4o3WTFp0iQmTpxIc3MzdXV1TJ48maamJsaMGeM8wB2l5NJL82qUElL42XLcjt9tRbNzKotD3NePaiCHvhPW5PE1HYaasvmzBZ1XCbLX39vhhx+uzz33nOd2Nm7c6EM1jm3btgHe56Fcvny5bzX59fP5VdPs2bNZs2ZN0lftwQcfrI8++mjM+w466KDF6ZwnKSJ1QKuq3ppyocZX48eP16amDK+Xuvhi+OMfYflyFofwPORx48YFXUJWiUjSvCXKK6SfWROsjPL6X/8FP/sZrFnTM1K5ePFi32ryK2dz5szxpR2AadOm+daWX1LJK/j7HpspG6lMgU1s7p+gl5AyIRA1SsnIkUFXYxKwcyqLmF31nZeCzqudU5mEdSj9E4YlpEwIuOdScu21QVdiEgjDnHcmQDbhed4JQ2atU5mAdSj9F/QSUiZgkVHK733PRinzgJ07WcQCuurbeBN0Zq1TGYe0tTH8ssvoZx1KXwX9KcoEzEYp84rltUjZoe+8FXRmrVMZQ0+H8u9/tw6lj8IwNG8CZKOUecXyWsTs0HdeCkNmQ9WpLCkpoV+/fp7bqaqqyvzJ27Yx8Lvfpezvf+ezX/2KkjPOwENrPf75z3/60Iq//LqKvKurK+XHBv2CN/7p7u6mtbU15cf3ueEGykpL+ezKK9Go5xX6ldb5zPJaONLJa8WDD1JeUsJnX/nKDlmFcOY1jFdsByXozIaqUxm4bdsYeN55lC1cyGe/+hXt06aBO6WQ8UfQ53uYYMgHH1D20EN0XHwxusceQZdjUmR5LUKqlP3lL3Qdfzw6bFjQ1Zg0BZ1Z61RGxOpQGl+FYWjeBKPi1luhtJSOH/wg6FJMiiyvxalk6VJK3n2X9iuuCLoUk6YwZNY6lWAdyhwK+gVvcs9GKfOX5bX4lP31r2hJCV2nnRZ0KSYDQWfWOpXWocypoF/wJvdslDJ/WV6LjB36zntBZzarexeRk0XkHRF5T0Suy+a+MmIdypyydYTDLRt57RmlvOgiG6XMM4nyapkNB78zGzn03WlzU+alMOQ1a51KESkFfg1MAQ4CponIQZm2V/Lww1Tsvz99+venYv/9KXn44YzaKX/0UQYdeiiDd92VwXvvTdkLL1iHMoeCXkLKxOZ3XsvmzqX/QQfR/9BDobOT7n339atUk0O2TGN4ZSOz/SZPRoGKX/yCsrlzfarU5FLQec3m4e8jgffclVAQkYeB04Gl6TZU8vDDlF9xBfLZZwDIRx9R/r3v0blqFd1Tp+78+La2mO2Uz59Pv1tuQbZvdzZs346Wl6NldhZArtibUWj5lteyuXPpc+WVSNTMCX1qa2HIEDrPOcevek0OWF5DLWuZldWr6XPllQCW2TwTdGazufeRwEdRt1e629JW9l//1dOhjJC2Nsqvv54+Y8fu9LXLscfG/Op/ww2fdygj7XR00O+mmzIpy6QpDEPzJi7f8lpRV7dDhxJAtm2joq4u7bbq6+upqamhtraWxsZGJkyYwPTp0wNrp9BrimaHv0MvdJkN42u60GuKFoa8ZrNTGeun0J0eJHKpiDSJSNP69etjN7RyZcztCrT/4Q87fbXOnh3za6edu0pWrUrxRzJeBT00b+JKO68bNmyI3VCcvMbbnkhDQwONjY1UVlYyc+ZMFixYQHV1NUuXpjcY41c7hV5Tb3b4O9SSZjaVvIJ/mQ3ja7rQa+ot6Lxm87jvSmDPqNujgNW9H6Sqs4BZAOPGjYvZ79NRo5CPPtp5+5570n3WWTtt74gzYXn3T39KaYyQdNuScTkRhjm0TFxp53Xs2LHp5XXUqLSLUnV2UVVVRWtrK11dXahqz/Zct1PoNUWzvIZe0symklfwL7NhfE0Xek3RwpDZbO79n8B+IvJvIlIBfAN4PJOGOm+8Ee3ff4dt2r8/nTfemFY722pr0V7LQGq/fmyrrc2kLJOBoD9Fmbh8y2t7XV3MnLVncPh70qRJTJw4kebmZurq6pg8eTJNTU2MGTMmkHYKvabebKQy1EKX2TC+pgu9pt6Czqt47RknbFxkKnAHUAr8TlV/lujx48aN0xdffDHmfSUPP+ycW7lyJTpqFJ033kj3N74R87HbEiytWP7oo/S76SZKVq2ie+RIttXW0nH22XEfn6itdMyfP9+Xdvzk13rkf/rTn1i3bl3SkzbGjh2rixYtinnfgAEDFqvqeF8KMhlJN69jx47VhQsXxryvbO5c5zwtN6/tdXVxT/gfOHCgt8JN2kQkad4S5RUss2GQTmYT5RVSz6zlNfdSySuE4z02q5c9q+rTwNN+tNX9jW/QHqcTmY6Os89O2Ik02ROGoXkTn5957TznHLtqNM9ZXsPPMmuihSGzNpeOyamgX/DGmNRZXo3JL0Fn1v5imJzycr5H6FdoMqbAeDmn0vJqTO4F/R5rnUqTM17m0PJ79QhjTGJe5qm0vBqTe2F4j7VOpckpD5+ielaPUNV2ILJ6hDEmSzyMVFpejQlA0O+xoTqn8tVXX93Qr1+/D5M8bCgQfxbXYIStplzXs1cqD1q8ePH8kpKSoXHu7isiTVG3Z7lzrEXEWj3iqPTKNH567bXXNlRWVlpe/RG6zCbJKyTOrOU1ZFLMK4QvH2GrB0KYVwjHe2yoOpWqOizZY0SkKWzTWIStprDVE6GqJ3t4ekorvpjcsbz6J4w1WV4LSyp5hfC9FsNWD4SzJghHZu3wt8kXKa34YowJBcurMfnFl8xap9LkC99WjzDGZJ3l1Zj84ktmQ3X4O0Wzkj8k58JWU9jq8UxVO0XkSmA+n68e8XbAZZnkwvhatJqyzPKa18L2WgxbPRDOmjzxK7NZXabRGGOMMcYUBzv8bYwxxhhjPLNOpTHGGGOM8SxvOpVhW/JLRPYUkUYRaRaRt0XkqqBrihCRUhF5TUSeDLoWU7wss6mxvJowsLymzjIbX150KkO65FcnMF1VDwSOBq4IQU0RVwHNQRdhipdlNi2WVxMoy2vaLLNx5EWnkhAu+aWqa1T1Vff7FpwX2MggawIQkVHAV4F7g67FFDXLbAosryYkLK8psswmli+dyljLBwX+4ooQkdHA4cA/Ai4F4A7gR0B3wHWY4maZTc0dWF5N8CyvqbsDy2xc+dKpDO2SXyIyEPgzcLWqbg24llOAdaq6OMg6jMEym0odllcTFpbX1GqxzCaRL53KUC75JSLlOC/2B1X1L0HXAxwHnCYiK3AOX5wkIn8MtiRTpCyzyVleTVhYXlNjmU0iLyY/F5EyYBkwEViFs5zQeUGu0CAiAjwAbFTVq4OqIx4RORH4oaqeEnAppghZZtNjeTVBsrymzzIbW16MVKpqJxBZPqgZmBuCJb+OA76F80nldfdrasA1GRMKlllj8ofl1fglL0YqjTHGGGNMuOXFSKUxxhhjjAk361QaY4wxxhjPrFNpjDHGGGM8s06lMcYYY4zxzDqVxhhjjDHGs6LtVIrI/xORt0XkTXeqgqOCrilTInKpiDwSdXuQiCwXkX8Lsi5j/GJ5NSZ/WF6LV1F2KkXkGOAUYKyqHgpMYsd1T/PNPcAoEZnk3r4R+J2qfhBgTcb4wvJqTP6wvBa3ouxUAiOADaq6HUBVN6jqagARGSciL4jIYhGZLyIjora/ISIvicgvReQtd/uFIvKrSMMi8qQ70z4iMtl9/Ksi8qi7hikiskJEbnC3LxGRL7rbB4rIfe62N0Xk64naiVBnstHvAXeIyHicVRF+mcXfnzG5ZHk1Jn9YXotYsXYqnwH2FJFlInKXiEyAnnVG/wc4S1XHAb8DfuY+5z7gP1T1mFR2ICJDgeuBSao6FmgCrol6yAZ3+2+AH7rbaoEtqnqI+wnvuRTaAUBV38RZDaHBrbM91V+GMSFneTUmf1hei1hZ0AUEQVVbRWQccAJQAzwiItfhvKAOBp4VEYBSYI2I7AIMVtUX3Cb+AExJspujgYOAF922KoCXou7/i/vvYuBr7veTgG9E1blJRE5J0k60XwNTVLUxSW3G5A3LqzH5w/Ja3IqyUwmgql3A88DzIrIEuADnBfh2709LIjIYiLeeZSc7jvj2jTwNeFZVp8V53nb33y4+/3+QGPtJ1k60bvfLmIJieTUmf1hei1dRHv4WkQNEZL+oTYcBHwLvAMPEOdEYESkXkTGquhnYIiLHu48/P+q5K4DDRKRERPYEjnS3vwwcJyL7um31F5H9k5T2DHBlVJ1DMmzHmIJheTUmf1hei1tRdiqBgcADIrJURN7EGf6uc8+TOAv4uYi8AbwOHOs+5yLg1yLyErAtqq0XgQ+AJcCtwKsAqroeuBCY4+7jZeCLSeq6GRgiIm+5+6/JsB1jConl1Zj8YXktYuJc2GTSISKjgSdV9eCgazHGJGZ5NSZ/WF7zW7GOVBpjjDHGGB/ZSKUxxhhjjPHMRiqNMcYYY4xn1qk0xhhjjDGeWafSGGOMMcZ4Zp1KY4wxxhjjmXUqjTHGGGOMZ/8fDT7m7P526pYAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 720x216 with 6 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "def compute_plot_D_P(X, Y, ax, step_sizes_sigma=np.array([[1, 1], [0, 1], [1, 0]]),\n",
    "                     weights_mul=np.array([1, 1, 1]), title='',\n",
    "                     global_constraints=False, band_rad=0.25):\n",
    "    D, P = librosa.sequence.dtw(X, Y, metric='euclidean', weights_mul=weights_mul,\n",
    "                    step_sizes_sigma=step_sizes_sigma, \n",
    "                    global_constraints=global_constraints, band_rad=band_rad)\n",
    "    D_max = np.nanmax(D[D != np.inf])\n",
    "    libfmp.c3.plot_matrix_with_points(D, P, linestyle='-', \n",
    "        ax=[ax], aspect='equal', clim=[0, D_max],\n",
    "        title= title, xlabel='Sequence Y', ylabel='Sequence X');\n",
    "    for x, y in zip(*np.where(np.isinf(D))):\n",
    "        plt.text(y, x, '$\\infty$', horizontalalignment='center', verticalalignment='center')\n",
    "    \n",
    "X = [1, 3, 9, 2, 1, 3, 9, 9]\n",
    "Y = [2, 0, 0, 9, 1, 7]  \n",
    "\n",
    "print(r'Accumulated cost matrix and optimal warping path for different step size conditions:')\n",
    "plt.figure(figsize=(10, 3))\n",
    "\n",
    "ax = plt.subplot(1, 3, 1)\n",
    "step_sizes_sigma = np.array([[1, 0], [0, 1], [1, 1]])\n",
    "title='Step sizes:'+''.join(str(s) for s in step_sizes_sigma)\n",
    "compute_plot_D_P(X, Y, ax=ax, step_sizes_sigma=step_sizes_sigma, title=title)\n",
    "\n",
    "ax = plt.subplot(1, 3, 2)\n",
    "step_sizes_sigma = np.array([[1, 1], [2, 1], [1, 2]])\n",
    "title='Step sizes:'+''.join(str(s) for s in step_sizes_sigma)\n",
    "compute_plot_D_P(X, Y, ax=ax, step_sizes_sigma=step_sizes_sigma, title=title)\n",
    "\n",
    "ax = plt.subplot(1, 3, 3)\n",
    "step_sizes_sigma = np.array([[1, 1], [3, 1], [1, 3]])\n",
    "title='Step sizes:'+''.join(str(s) for s in step_sizes_sigma)                                                       \n",
    "compute_plot_D_P(X, Y, ax=ax, step_sizes_sigma=step_sizes_sigma, title=title)\n",
    "\n",
    "plt.tight_layout()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Local Weights\n",
    "\n",
    "To favor the vertical, horizontal, or diagonal direction in the alignment, one can introduce additional **local weights** $w_{\\mathrm{d}},w_{\\mathrm{h}},w_{\\mathrm{v}}\\in\\mathbb{R}$. To compute the accumulated cost matrix $\\mathbf{D}$, one uses the following initialization and recursion:\n",
    "\n",
    "\\begin{eqnarray}\n",
    "  \\mathbf{D}(1,1)&:=&\\mathbf{C}(1,1)\\\\\n",
    "  \\mathbf{D}(n,1)&=&\\sum_{k=1}^{n} w_{\\mathrm{h}}\\cdot \\mathbf{C}(k,1) \\,\\,\\mbox{for}\\,\\, n\\in [2:N]\\\\\n",
    "  \\mathbf{D}(1,m)&=&\\sum_{k=1}^{m} w_{\\mathrm{v}}\\cdot \\mathbf{C}(1,k) \\,\\,\\mbox{for}\\,\\, m\\in [2:M]\\\\\n",
    "  \\mathbf{D}(n,m)&=&\\min \\left\\{\n",
    "   \\begin{array}{l}\n",
    "    \\mathbf{D}(n-1,m-1) + w_{\\mathrm{d}}\\cdot \\mathbf{C}(n,m)\\\\\n",
    "    \\mathbf{D}(n-1,m)   + w_{\\mathrm{v}}\\cdot \\mathbf{C}(n,m)\\\\\n",
    "    \\mathbf{D}(n,m-1)   + w_{\\mathrm{h}}\\cdot \\mathbf{C}(n,m)\n",
    "   \\end{array}\n",
    "   \\right.            \n",
    "\\end{eqnarray}\n",
    "\n",
    "for $n\\in[2:N]$ and $m\\in[2:M]$. The case $w_{\\mathrm{d}}=w_{\\mathrm{h}}=w_{\\mathrm{v}}=1$ \n",
    "reduces to classical DTW. Note that in the classical case one has a preference for the diagonal alignment direction, since one diagonal step (cost of one cell) corresponds to the combination of one horizontal and one vertical step (cost of two cells). To balance out this preference, one often chooses $w_{\\mathrm{d}}=2$ and $w_{\\mathrm{h}}=w_{\\mathrm{v}}=1$. Similarly, one can introduce weights for other step size conditions. In the following code cell, we call the function `librosa.sequence.dtw` with different weight settings."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2024-02-15T08:50:53.920511Z",
     "iopub.status.busy": "2024-02-15T08:50:53.920311Z",
     "iopub.status.idle": "2024-02-15T08:50:54.340945Z",
     "shell.execute_reply": "2024-02-15T08:50:54.340340Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Accumulated cost matrix and optimal warping path for different local weights:\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAroAAADgCAYAAADsd8WhAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAABIxElEQVR4nO3deZwcdbX//9eZfbLvYbIgiwGEyBr8sl1kCfzCooAiJCgisngv+yKySgISLyIEVAIaEiQiEkOAy76ESAQUEBLZQgJRCJCFbGbPTE9m5vz+qOqhM+meXqqq+9Pd5/l4zGPSNTWfPjPJO32quurzEVXFGGOMMcaYUlNR6AKMMcYYY4yJgjW6xhhjjDGmJFmja4wxxhhjSpI1usYYY4wxpiRZo2uMMcYYY0qSNbrGGGOMMaYkWaObIRH5rog8n+G+PxCRV6KuyX+uRSLSKCL35+P5ghCRG0Rkk4ioiFQVuh5TeiynwVlOTdQsp8FZTjNX0o2uiFwtIk932LYwxbbRnY2lqg+o6tEh1TVbRM4OYyzfN1T19E6eb5KIfCAibSLygzS1HS4iL4rIOhFZlO6JsxlbVccCe6Qb05QXyymIyC4i8piIrBSR/4jIcyKyaye1WU5NXllO258vm9fTS0TkIxFZLyJLReT2VE2piPQTkb+JyGoRWSsir4rIwanGtpxmrqQbXeAl4GARqQQQke2AamDfDtu+7O9bqt4GzgPmZrDvJuBe4IoIxjYmGcsp9AIeB3YFBgL/AB7rZH/Lqck3y6knmyw9Aeyrqj2A4cBewEUp9t0I/BDoD/QGfgE8YWdrgyv1RvcNvCDu7T8+FHgR+KDDtn+r6lIR6SkiU0RkmYgsEZGbEgK81dsnInK0f1S3TkTuEpG/djyqFJFbRWSNiHwsIsf428YD/wXcKSIbReRO8dwuIiv88d4RkeFh/RJUdaKqzgKaMtj3H6p6P/BR2GMbk0LZ59TP3RRV/Y+qbgFuB3YVkb6d7G85NflU9jmFrF9P/62qa+M/AtCGdyCQbN8mVf1AVdv8fVvxGt4+oRRexkq60VXVZuB1vPDhf34ZeKXDtvjR51SgBe8f4j7A0cA2b4mISD9gBnA10Bcv6Ad12O3/+dv7AbcAU0REVPVav4YLVLWbql7gP8+hwC54Z3ZOBVb7z3WaiLyT8y/BGMdZTpM6FPhcVVeHOKYxObOc5sZ/zvXAKrwzur9Ls/87eE3048BkVV0RfZWlraQbXd9f+SKE/4UXipc7bPuriAwEjgEuUdVN/j+u24Fk1xodC8xT1UdUtQX4NfB5h30+UdV7VLUVL/ANeG9JJrMF6A7sBoiqzlfVZQCq+idV3TPrn9qY4mI59YnIEGAicFkY4xkTIstplvzn7IHXeP8WWJ5m/z2BHsBpeAcRJqByaHRfAg4Rkd5Af1VdCPwdOMjfNtzf50t4b8ssE+9C8LV4R14Dkow5CPgs/kBVFVjcYZ/PE76+2f9jt2QFqupfgDvxXtyWi3exe49sf1BjipjlFBCR/sDzwF2q+mCYYxsTAstpjvzf1Tzgrgz2bfLzf5WI7BV5cSWuHBrdV4GewLnA3wBUdT2w1N+2VFU/xgtaDOinqr38jx6qmuyuxmXAkPgDEZHExxnQbTao/lpV98O7i3IXMr/JxJhSUPY59RuF54HHVXV8WOMaE6Kyz2lAVcDOWexfDewUUS1lo+QbXVVtBN7Eexvw5YQvveJve8nfbxnei8xtItJDRCpEZGcR+XqSYZ8CvioiJ4p3R+T5wHZZlLWchH+8IrK/iPw/EanGu5u6Ce9C9FCISI2I1OFd4F4tInUikvTv3v+56/ACJv6+NWGMbUwq5Z5T/4zTc8DfVPWqDPa3nJq8K/ec+uNn83p6togM8P+8O951yLNS7HuAiBzij18vIlfiXZ7xeli1l6ty+Y/ur3hvmSRe7/Kyvy1xGpTvAzXA+8AavAvkGzoOpqqrgO/gXRS/GtgdL/yxDOv5FXCyeHeQ/hrvepx7/Of8xB/zVmifWHtehuOm8jzQiHeB/yT/z4em2PdQ/+tPA9v7f+5sYu9sxjamM+Wc05OA/YEzxbt7PP6xfYr9LaemUMo5p5Bdlg4G3hWRTXhZfRq4JsW+tXiXW6wGluBdu3ycqi4NWG/ZE+9yGBOEfzS3GPiuqr6Y5+f+AO8/j0dV9Yx8Pne2RGQs3lF/LdDVv7HAmLywnGbGcmoKyXKaGctp5qzRzZGI/H94byk04l3/cz6wk//WjjHGAZZTY9xnOTVRKpdLF6JwIPBvvLnxvgGcaKEsD/41Wf8QkbdFZJ6I3OBvHyfexOhv+R/HFrpWYzktdyJSKSL/FJEn/cd9RGSmeEvVzvRvAjSFZTk1kbEzusZkyb8ruKuqbvRveHgFuBgYBWxU1VsLWqAxpp2IXAaMAHqo6vEicgvwH1W9WUSuAnqr6pWFrdIYExU7o2tMltSz0X9Y7X/YEaMxjhFv8Y3jgMkJm0/AW3QA//OJeS7LGJNH1ugakwP/7dC3gBXATFWNTwFzgXhrq99rb4kaU3B3AD8B2hK2DUxYKWsZyRcxMMaUiKpCF5Covr5ee/QIvoBJXV1dCNWEp5QvD2lra0u/UwbC+B2tXbuWTZs2Sbr9RCTdk83Dm3sxbpKqTkrcwb/DdW8R6QU8KiLDgbuBn+Gd3f0ZcBvww8x/guLQpUsX7dWrV+BxqqurgxdDePkKa5yWlpZQxglzLO9qm+DC+h2tWrVqlar272yfoDkVkeOBFao6R0QOy7XWYhXW62m3bkkXIMtaZWVlKOOE9W85rHEAKirCOWfo0uspwIcffhhGTgGeU9VRoRSVA6ca3R49enDqqacGHmf48OEhVBOeWCzT6QCLT1NTU/qdMhDG7+juu+/OeN/O/pPzl18ckck4qrpWRGYDoxKvzRWRe4AnMy6oiPTq1Ytzzz038DiDBg0KoRpobm4OZZywcrp69epQxglzrJqalGtJZCWs3/WkSZM+yWS/gDk9GPimf1NoHdBDRP6Ityxsg6ouE5EGvHdlSk6PHj0YM2ZM4HEOOeSQEKqBnj17hjJOWE1lWAfaAF26dAllHJdeTwFGjhwZOKcAqtovzfffC8QPTIf72/oAfwZ2ABYBp6jqGv9rVwNn4S0EcpGqPtfZ+HbpgilLFRUVKT/SEZH+/plcRKQeGAks8F80404C3ouidmPKRZCcqurVqjpEVXcARgN/UdXvAY8D8TlSzwAei6p+Y8pBZznN8MDkPrybuRNdBcxS1WF4q8ldBe0rzI3GW955FHCXiHT6doFTZ3SNyZeAb1s1AFP9cFUA01X1SRG5X0T2xrt0YRHwo6B1GlPOwnx7OcHNwHQROQv4FG9VLmNMjoLmVFVfEpEdOmw+ATjM//NUYDZwpb99mqrGgI9F5F/A14BXU41vja4pOyIS6O0vVX0H2CfJ9tOD1GWM+ULQnCZS1dl4L5So6mrgyFAGNqbMhZnTDra6aVRE4jeNDgZeS9hvsb8tpUgvXRCRS/0J9d8TkQdFxK27xEzZEpGUH+XGcmpcZTn9guXUuKqznPpZ7ScibyZ8BLnJI1n4O70hLrIzuiIyGLgI2F1VG0VkOt51FfdF9Zym9Ozx9tscPnMmPdetY13Pnrx41FHM22uvwONGdARadCynJgz7L1zISW+8QZ+NG/lPt248uv/+vDFsWOBxLacey6kJy8BZs9h5yhTqVq6kqX9//n3WWSw/MtgbHBnkdFWmN3gnSHXT6GJgaMJ+Q4ClndaX5RNnqwqoF5EqoEu6YoxJtMfbb3PcY4/Ra906BOi1bh3HPfYYe7z9dqBxMzj6LDeWU5Oz/Rcu5PSXX6bvxo0I0HfjRk5/+WX2X7gw0LiW021YTk0gA2fN4isTJlC/YgWiSv2KFXxlwgQGzpqV85jpchogq6luGn0cGC0itSKyIzAM+EdnA0XW6KrqEuBWvIv9lwHrVPX5qJ7PlJ7DZ86kZsuWrbbVbNnC4TNnBh474B2iJcNyaoI66Y03qO0w329tSwsnvfFG4LEtpx7LqQnDzlOmUNlh6rHKWIydp0wJNG7QWRdE5EG8m8l2FZHF/o2iNwNHichC4Cj/Mao6D5gOvA88C5zvz2ufur5AP13nhffGuztuR2AQ0FVEvpdkv3Pj1200NjZGVY4pQj3XrctqezbsBdSTS043b96c7zKNw/ps3JjV9mxYTj32emrCULdyZVbbMxW00VXVMaraoKrV/pSAU1R1taoeqarD/M//Sdh/vKrurKq7quozaesL9NN1biTwsaquVNUtwCPAQR13UtVJqjpCVUfU19dHWI4pNptTTMK9LuDE4/aW6FayzmlYk6Ob0rA2xb+H/wRcUctyuhV7PTWBSHMzbSkWyWjq3+niZ52PG92lC6GJcnqxT4EDRKQL0Ig3ncubET6fKSHDFiygvrGRNhEqEpYzbK6u5sWjjgo8frmdEeqE5dTkrO+GDVS2tqJsfSt0rKqKR/ffP/D4ltN2llOTM2lu5qs33khlczNtVVVUJFxq1Fpby7/POivQ+K7nNMprdF8HZgBzgXf955rU6TcZg9fknjxtGssGDeKpb3yDtT17osDanj156oQTQpl1weWjz3yynJpc9d2wgcueeIJK4P9GjGB1t24osLpbN+7/r/8KZdYFy6nHcmpyFW9y+7/2Ggsuvpj3r7iCxgEDUBEaBwxg/mWXBZ51oZzP6KKqY4GxUT6HKS3xJvfz7bbjT2ecQay+nrdDODOUSKKb4LooWU5NtuJNbv2WLdx+3HF81q8fz+67b6jPYTndmuXUZKtjk7vkG98ACNzYbvUcRZBTWxnNOCNZkxsVF44yjSlGyZrcqFhOjclNqiY3kudyPKfW6Bon5LPJBfevKTLGRflscsFyakwu8tnkgvs5tUbXFNzO77/PCXlscovhrRZjXNNn/XouzGOTazk1Jnv5bnKLIaduV2dK3s7vv88J99+ftyY3zuUL541xTZ/167no0Ufz1uTGBcmpiNSJyD9E5G0RmSciN/jbx4nIEhF5y/84NvIfxJg8kOZmvvLTn+atyW1/3nK+GS1bdXV1DB8+PPA4AwYMCKEaiHVYQaSUrFixIv1OEYs3uSsGDeJPp5+etyYX3H+rxWV1dXXsuuuugcepqakJoRpobm52ahwXBfm/rPe6dZz3f/9HXXMzvznhBJYPGEA4f3PpBcxpDDhCVTeKSDXwiojEJ5e/XVVvDVygw+rr6/nqV78aeJywXk/D+j83rP83wpwPfMOGDaGM09JhhcFsSHMze4wbR59XX2XeBRfw2ahR0GFl0ai4/nrqVKNrykdik/vQ2WcTy3NQghxlikgd8BJQi5ehGao6VkT6AH8GdgAWAaeo6prAxRpTIL3XreO8hx6iLhbjtyefzOLevfP6/EFyqqoKxJdnq/Y/NPV3GFOc4k1uv1df5cNLL+Wzo4/O7/M7cNa2M2634aYkbdPk5nkFn/g1RQGWFo2fKdoL2BsYJSIHAFcBs1R1GDDLf2xMUerY5C4ZODCvz59BTvvFl7v1P85NMkaliLwFrABm+vPRAlwgIu+IyL3+8rrGFKWOTe7Sb34zv8+fJqcunO0tfAWmrBS6yY0Lcj2RepKdKToBmOpvnwqcGEHpxkSu0E1uXJqcroovd+t/bLOAgqq2qurewBDgayIyHLgb2BnvIHUZcFvefiBjQlToJre9Dsev0bVG1+SNK00uEPjoM8WZooGqugzA/xzOxW3G5JErTS4Ez2mcqq4FZgOjVHW53wC3AfcAXwu9cGMi5kqTC53n1M7omrLhUpMbxluiKc4UGVPUXGpyg15iJCL9RaSX/+d6YCSwQEQaEnY7CXgvivqNiYpLTW4xXLpgN6OViO1feYU9p02jy+rVbO7bl3dGj+bTQw7Jepxhb7zBQU8+Sfc1a9jQuzd/P/54FuawBO9uc+dy6LPP0mPtWjZ37Ur95s0sHzKk4E1uXJq3U1ap6ohMxlHVtSIyGxgFLBeRBlVd5r+YFn5qC+OU3ebM4ZBnnqHH2rWs79WLV445hgX77Zf1OHvNm8fRs2fTa/161vbowfOHHcbbe+yR9Tj7zJ/Psa+8Qu8NG1jXtStVra1UqBa8yY0L+LZnAzBVRCrxTupMV9UnReR+Edkb73KjRcCPgtZpTNQGvPACO02eTO2KFbRVV1PZ3FzwJjfOhcsTOmONbgnY/pVX2H/SJKr86ZG6rlrF/pO8y9WyaXaHvfEGR06bRrU/JUmPNWs4cto0gKya3d3mzmXUww+3j9N10ybaRHh7//2daHIh2HQoItIf2OI3ufEzRb8AHgfOAG72Pz8WQqmmROw2Zw5Hz5jRnouea9dy9IwZAFk1u3vNm8dJTz9NjT8VUe/16znp6acBsmp295k/n1Nmzmwfp9emTSjw9EEHOdHkQrCcquo7wD5Jtp8epCZj8m3ACy+w6623UulPE1jZ3ExbVRUtIU6RFoQLZ207Y41uCdhz2rT2JjeuqrmZEZMn03/BgqTf09jYuM22Xd98s/1FOK56yxaOmD6dwR99lHSc1iTz/u3+z39uM06FKgf+5S+8e8ABnf4s+RDCBfKpzhS9CkwXkbOAT4HvBK/WlIpDnnkmab6Oevhhhnz8cdLvaUoy/+3e773X3pzG1bS0cOIzz7DDZ58lHae1tXWbbfvOn7/NOAIc+O67/KU0cmpMSdhp8uT2JjeuoqWFnSZPZsXIkQWqylMMObVGtwR0Wb066faqpiYGv/FG0q+1tbVts606xaTy1bEYO73zTvIn122npaxOMfF+j7Vrk49RABGdKVoNHBmgLFPCUv37r25u5svz5iX9WrKc1qSYBL5myxZ2//DDpF/TJDmtTTFO75Amvw+D62eKjMmH2hQLPKXanm+u59Qa3RKwpUsXajZt2mb75n79ePLOO5N+T7KV0c4YO5Yea7Zd32BD795MveGGpOM0NTVts+3cn/+cnkle1Nf36pV0jEJw/QjUlJ7N3bvTNUkTub5XLyZfd13S71md5CD2iokT6b1+/Tbb1/bowS/PPz/pOMlWRrv2nnvok6SeNd27Jx2jECynptxJc3P7NbkdxUJatS4o13Pqdhtu0trlqaeo2bSJtg5HVC01NbwzenRWY/39+OPZUl291bYt1dX8/fjjsxrnpVGjko7z0qhRWY0TlRAWjDAmK4MWLaJm8+ZtluXaUl3NK8cck9VYzx92GM1VW5+jaK6q4vnDDstqnKcPOSTpOE/ncBNrFCynptzFZ1eIX5ObqLW2lo/OPrtAlX3BZl0wkdrlqafY5/77+fSAA1iy777sOX16oFkX4jecBZ11YcG++wK0z7qwvlcvXho1qn27C1wInykPgxYt4lv33MOG3r2Ze8gh7D97dqBZF+I3nAWddeGfX/kKQPusC2u6d+fpQw5p3+4Cy6kpVx2nEGvp0qV91oXYgAF8dPbZBb8+N871nFqjW6QSm9zXLrwQrazk00MPDTzuwv33z2k6sY4W7LuvU41tR66/1WJKQ7zJ3dS9Ow/9z/+wsWdP3g7hjOnbe+yR03RiHf3zK19xqrHtyHJqylGqeXJdaWw7cj2n1ugWoWRNrslc/K0WY6KUrMk1mbOcmnLk0mIQmSiGnLpdndmGNbnhcHldblP8rMkNh+XUlJNia3LjOstpJlkVkUtFZJ6IvCciD4pInYj0EZGZIrLQ/9w71/qs0S0i1uSGx+UL501xsyY3PJZTUy6KtcmFznOaLqsiMhi4CBihqsOBSmA0cBUwS1WHAbP8xzmxSxeKhDW54SmGt1pMceq/cCFHWpMbCsupKRfS3Mwu115LnyJsckPKaRVQLyJbgC7AUuBq4DD/61OB2cCVuQ7ujPr6eoYPHx54nGRzRhazoQ89xJfvv58Fe+7JkyecgC5ZEnjMOXPmhFAZLFy4MJRx/vWvfwUeY+PGjRnva2995q6uro5ddtkl8DgbQlqYoDnFAiX5Hqf3/Pl8bcIE1nftyn2nn86GigoI+DO+++67gb4/7sMUC0kUapxMWU5zV19fz1e/+tXA44R1sFHdYcrJXHUJadnbbF4v0tmSYvGVTEhzM7tffz19Xn2V+RddxOJjjoEA40HyefKjlEFO+4nImwmPJ6nqJABVXSIit+KtJtoIPK+qz4vIQFVd5u+zTERynjTYqUbXbGvoQw/x5bvvZsXXv86To0bZmdyQ2JkiE6be8+fztbFjifXuzX2nnsqGHj0KXVJJCJJTEakDXgJq8V7rZqjqWBHpA/wZ2AFYBJyiqtuulGNMxOJNbt9XX2XhZZex2JG55rOVQU5XqeqIZF/wr709AdgRWAs8JCLfC7W+MAcz4Upsct+/7jprckNkN7mYsCQ2ua/9/OfW5IYoYE5jwBGquhewNzBKRA4gxGv/jMlVxyZ32QknFLqknAW8GW0k8LGqrlTVLcAjwEHAchFp8MdvAHI+TW2NrqOsyY2OrbhkwtKxyW3q27fQJZWMoDlVT/z96Wr/Q/HOHk31t08FToygfGNSKrUmN+DKaJ8CB4hIF/G64iOB+cDjwBn+PmcAj+Vao1264CBrcqNnZ25NUNbkRi9oTkWkEpgDfBmYqKqvh3ntnzHZKqUmNy5ITv1MzgDmAi3AP4FJQDdguoichdcMfyfX54i00RWRXsBkYDjekfQPVfXVKJ+z2FmTmx925vYLltPsWZObH2lymvIGlzhVbQX29v+NPyoiwe92LhDLafErxSYXgr+equpYYGyHzTG8s7uBRX1G91fAs6p6sojU4E0bYRIMeOEFdpo8mbqVK2np2pXqjRutyY2YTVu0DctpGoNmz2a3P/yB+lWriPXsSdWmTTT1729NboQyyGnKG1w6UtW1IjIbGIV/7Z9/NjfQtX95ZjktQv1nzmTHe+6hdsUK2qqrqWxuLqkmtxheTyOrTkR6AIcCUwBUtVlV10b1fMVowAsvsNttt1G/YgWiSvXGjWhFBasOPNCa3IgFXMVlqIi8KCLz/dVcLva3jxORJSLylv9xbOQ/SECW0/QGzZ7NnnfeSZeVKxFV6taupbKlhY+PP96a3IgFzGl//ywoIlKPd9PLAkK89i9fLKfFqf/Mmezyy19St3w5okplczNtVVW0hDRFmiuCrowWtSjb8J2AlcDvReSfIjJZRLpG+HxFZ6fJk6nsMOevtLWx4733Fqii8hHwZrQW4HJV/QpwAHC+iOzuf+12Vd3b/3g6qvpDZDlNY7c//IGqjjlVZedHHy1QReUjYE4bgBdF5B3gDWCmqj4J3AwcJSILgaP8x66znBahHe+5Z5vX+IqWFna8554CVRSNgDejRV9fhGNXAfsCd6vqPsAmkkzjIiLnisibIvLm2rVrIyzHPXUrV2a13YQj6NGnqi5T1bn+nzfg3SE6OOKyo5J1TtesKa8pR+tXrcpquwlHCDl9R1X3UdU9VXW4qt7ob1+tqkeq6jD/838i/2GCs9fTIlSbYuGGVNuLUbqclvoZ3cXAYlV93X88Ay+oW1HVSao6QlVH9OrVK8Jy3NPSNfkBeVP//nmupPyEdfQpIjsA+wDxf+cXiMg7InKveBNhuy7rnPbuXQw/Vnga+/XLarsJj8tnifLMXk+LjDQ305ZiNbjYgNKa6KNsz+iq6ufAZyKyq7/pSOD9qJ6v2AydPr39mtxErbW1fHT22QWqqnykOfrsFz8r4n+cm2KMbsDDwCWquh64G9gZb3L6ZcBt+flpcmc5TW/B979PS23tVttaamtZ8P3vF6ii8uHyWaJ8spwWl/jsCvFrchO11tby8TnnFKiyaLh+RjfqWRcuBB7w7xD9CDgz4ucrCkOnT+fLv/0tK77+dVYdeCA73nsvdStX0tS/Px+dfTYrRo4sdIklLYy7uUWkGq/JfUBVHwFQ1eUJX78HeDKEcvPBctqJpYcdBsCeEydS2dREY//+LPj+99u3m2gUw93ceWY5LQIdpxBr6dKlfdaF2IABfHzOOaw86qhClxmaYshppI2uqr4FZDT9S7lIbHLjU4gtP/roQpdVdoIEU7xD1CnAfFWdkLC9IT4RPXAS8F6gIvPEcpre0sMOo+977zHwH//gL3azaN64/gKaT5ZT96WaJ7eUGttkXM+prYyWR8maXFMYAd9OORg4HXhXRN7yt10DjBGRvfEmc18E/CjIkxhT7lx429OYTJTqYhCZcD2n1ujmiTW57gj6VouqvgIkS3YxTCdmTFEohrdEjQFrcl3PqTW6eWBNrntcPwI1xlhOjfvKucmNcz2nKdtwEblcRLbpyESkr4hMibas0mFNrptcngolG5ZTU8pKJadgWS1F0tzMjldcUdZNLrg/vVhnZ3R3BeaIyPmq+jcAETkPuAJvze3QVVRUUNthGp9chDFGGAb88Y8M/e1vWXXEEfxr7FhqqoKdQK+rqwupsvDGCut3HcY4mzdvzmg/V6Y8CUnecyoizmQMoLm5OZRxYh1WMOqouroaEaFbt26d7te9e/dQ6glzrHQ153ucjRs3pt2nxHIKec5qWK+n1Snmgs1Wl5CWva0M6WTRli1bAn2/NDcz5Mor6f7yy3x85ZWs+da3CPqqGtZiPJm+FoahGHKasvNS1XNF5CDgThGZB+wGLAQOSriz3KQw4I9/ZOgdd/CfkSP58Kc/hYBNrgmXC0eZYbCcmlJWKjkFy2opkeZmhlxyCd3/+leWXX89K487rtAlFZTrOU3Xfb2Ht0b4KLybby63QKaX2OR+fNNN0Npa6JJMB64HM0uWU1OSSiynYFkteh2b3DWnngoZvENRylzPaWfX6H4PeAtvYuqd8eYFvUVE/iAipbV+XYi2aXLtTK5zXF/FJRuWU1OqguZURIaKyIsiMl9E5onIxf72cSKyRETe8j+OjfyHwbJaCpI2uWUuXU5deE3trAv7DnC4qn7iP54jIgcC/w28BuwUdXHFxprc4uH6EWgWLKemZAXMaQveGdO5ItIdLxsz/a/drqq3Bi4wO5bVImZNbmquv552do3uNrcPqqoCd4vIjEirKkLW5BYXF44yw2A5NaUsSE79SwKW+X/eICLzgcEhlZZLPZbVImVNbudcfz3NqRtT1ZVhF1Jsej/zDIMnTqRm+XJau3WjasMGa3KLRDFMcB0Gy2l4BrzwAgNmzaKysZEDRo/mo7PPZsXIkYUuq6SFmVMR2QHYB3gdb2XDC0Tk+8CbeGd9w7ndPUeWVff0ePJJBtxxB9Wff47W1FARi1mTm0QxvJ66XZ2jej/zDF8aP57azz9HVKnasAGtqGDtIYdYk1skXL6eyLhlwAsvsOutt1LV2IgAdcuXs+uttzLghRcKXVrJS5PTfiLyZsLHuSnG6AY8DFyiquuBu/Gukd0b74zvbfn5aUyx6PHkkwwaO5aaZcsQVSpiMdqqq2nt2rXQpTnJ9Wt0rdHNweCJE6lsatpqm7S1Mfi3vy1QRSZbLk9ubdyy0+TJVHaYY7cyFmOnyZMLVFH5SJPTVao6IuFjUsfvF5FqvCb3AVV9BEBVl6tqq6q2AfcAX8vnz2TcN+COO6jo8BpfsWULA+64ozAFOc71BSPSViAiXUTkpyJyj/94mIgcH31p7qpZvjyr7cYt8bdaXA1lLiyn0aldsSKr7SYcQXMq3qmkKcB8VZ2QsL0hYbeT8Kb8yhvLqvuqP/88q+3lLF1OXXhNzaSC3wMx4ED/8WLgpsgqKgKtKVYHah44MM+VmFy5/DZLjiynEYkNSD7zU6rtJjwBc3owcDpwRIepxG4RkXdF5B3gcODSCH+EZCyrDpPmZrSmJunXtmy3XZ6rKQ4hTAXYS0RmiMgCfzrAA0Wkj4jMFJGF/ufeudaXSaO7s6reAmwBUNVGvImuy9KAP/6x/ZrcRK11dSw5//wCVWWy5fLRZ44spxH56Oyzae2wlGprbS0fnX12gSoqH0FyqqqvqKqo6p6qurf/8bSqnq6qX/W3f7MACzZYVh0Vn10hfk1uora6OlZccklhCnNcCGd0fwU8q6q7AXsB84GrgFmqOgyY5T/Orb4M9mkWkXpAAURkZ7yj0bKz1RRi119PbLvtUBFi223HJ9dey5pjjil0iSZDAY8+U01EH9oRaA4spxFZMXIkH/z4x7TU16NA08CBfPDjH9usC3lQgu+8gGXVSR2nEFt60000NzSgIjQ3NLD0hhtYf7xdYZJMkDO6ItIDOBTvMiNUtVlV1wInAFP93aYCJ+ZaXyZTBIwFngWGisgDeG8H/SDXJyxWyebJXWP/6ItSCNOhpJqI/gd4R6A3i8hVeEegVwYuODOW0witGDmSXm+/Td+//53Xpk0rdDllIYScusqy6phU8+RaY5teCDndCVgJ/F5E9gLmABcDA+PvtqjqMgmwemDaRldVZ4rIXOAAvLdXLlbVVbk+YTGyxSBKT5AzQp1MRH8CcJi/21RgNnlqdC2nphQV+ZnbpCyrbrHFIILLIKf9ROTNhMeTEmZJqQL2BS5U1ddF5FcEuEwhmUxmXTgJaFHVp1T1SaBFRE4MswiXWZNbmsK6Rle2noh+qyNQIG93K5V7Tk1pKsFr6S2rDrEmNxwZXKPb2VSAi4HFqvq6/3gGXuO7PD5Div8552luMvnfYqyqros/8K+dGJvrExYTa3JLUwbTFuU6EX0hlW1OTWkqxWkAfZZVB1iTG46g04up6ufAZyKyq7/pSOB94HHgDH/bGcBjudaYSeeWrMqS7/isyS1tad5qWaWqI9J8/zYT0eMfgfrXEwU6As1BWebUlLZSvHQBy2rhxWLW5IYohJxeCDwgIjXAR8CZeDmZLiJnAZ8C38l18EzC9aaITAAm4t0leiHexcKhExFqO0zjk4u6urpA399ryhT63XEH644+mqW/+AW1jjS5Yfxuwh4r6O86zHGyCVuQM0IiySei54sj0JsJeASag7zmtCbFPJPZCGMMgObm5lDGSVdPdXU1FSJ079690/369esXSj0ADQ0N6XfKQCzm1k39c+fOzWi/Ij9zm0pesioiofy/2i3FvPHZqqysDGWcwP8mYjF6nncetX/9Kyt/9jMaTzuNoL+lFSEtHrN582anxslU0L8TVX0LSHZy6chAA/syqe5CoBn4M/AQ0ASU7ISxvaZMod/NN7PhmGP47Be/sDO5JSjo5Naknoj+ZuAoEVkIHOU/zpeyyqkpfSHk1FWW1UKJxeh51lnUzpzJyp/9jA2nnVboiopeupy6kNVMZl3YRMh3wLkqscldPmECtLQUuiQTkSBHoKr6CqkneA/lCDRb5ZRTUz5K8YyuZbVAEprc9bfcwoZvf7vQFZUM13OattEVkV2AHwM7JO6vqkdEV1b+bdPkVlVZo1vCXDjKDFO55NSUl1LLKVhWC6JDk9t0xhmwcWOhqyoZruc0k/flHwJ+C0wGWqMtpzCSNrmmZJXoRPQln1NTXko0p2BZza9kTa4JTTHkNJOOrkVV7468kgKxJrc8uR7MHJR0Tk15KsGcgmU1f6zJzQvXc5pJdU+IyHki0iAifeIfmT6BiFSKyD9F5MkAdYam2+OP86Wvf52dd9mFHffbz5rcMuXyhfM5Kqmcdn/iCXY8/HCG7bYbOx5+ON2feCKncXo+9RS7HH00e+y5J7scfTQ9n3oqp3F6P/MMfZ59lqrVqxl+/PH0fuaZnMYx2QmSUxEZKiIvish8EZknIhf72/uIyEwRWeh/7h35D7K1nLPqWk5dU/vww/Tdbz/6b7cdfffdl97HHmtNbh4U/c1ofDFh7xUJ2xRvfeJMXAzMB3pkUVckuj3+OAOuvZaKpiYAKtevRysq2HT44dbklpFieKslByWT0+5PPMHA665rz2n10qUMvO46ADZ84xsZj9PzqacYPG5c+zg1y5YxeNw4ANYdd1zG4/R+5hm+NH48lf44tZ9/zpfGjwdgzTHHZDyOyU4IOW0BLlfVuSLSHZgjIjOBHwCzVPVmEbkK78awvCzV7QuSVWdy6prahx+mx+WXI42NAFQuWULlkiVsHj3amtwIFcPraSazLuyY6+AiMgQ4DhgPXJbrOGHpe9tt7S96cdLWRt877mDjSScVqCpTCC4cZYaplHLab8KEbXJa0dRE//Hj0SRzcbakuGm04eabk47TcPPNWY0z9Lbb2pvcuMqmJgZPnGiNbsSC5NRfhju+JPcGEZkPDAZOAA7zd5sKzCaPjW6uWXUtp67p9vOftze5iWpfeQW77Sxarr+eZjLrQhe8UG2vqueKyDBgV3+N7nTuAH4CpJxhXbzlVc+F8CZHT6Vq2bKstpvS5foRaLbymdNBgwYFL7gTKXO6Zg2DLr00+Phr17L9FVek3zGNmuXLA49hOpcmp/1E5M2Ex5NUdVKyHUVkB2Af4HVgoN8E469iOCCkcjMSIKt34FBOXVOxZElW2014XH89zeT9+t/jrdpykP94Md5do52GUkSOB1ao6hwROSzVfv5/TJMA9thjD82gnpy1de9O5fr122xvibjBNu5x/Qg0B3nL6fDhwyPNaWvfvlStWrXN9pb+/Vl8333bbN+yZUvScXY45xyqV67cdv/+/Vl0zz3bbE+1wtqw886jJkk9zQMHJt3fhCdNTtMu1e2P0Q1vue5LVHW9A9nPOqsu5tQ1bYMGUZmkqW0bPLgA1ZQXBzLVqUwa3Z1V9VQRGQOgqo2S2U91MPBNf8WoOqCHiPxRVb8XoN6c9Zoypf2aXGlra9/eVlfH6ssvL0RJpkCK4ZqiHJRETuvmzqVi/XpUBNEvXqfb6upYeeWVNH/5y9t8T6oG9fPLL9/qGt34OJ9ffjmxnXfeZv9Uy+Quvvjira7RBWitq2PJ+baYVZTCyKmIVOM1uQ+o6iP+5uUi0uCfzW0Awlm/NXO5ZNWpnDonFqOtd+9tGl2tr2fjNdcUqKjyUAyvp5lU1ywi9XgXyyMiOwNpF05X1atVdYiq7gCMBv5SyCa3fXaFm29my6BBqAhbBg1ixfjxbPzmNwtRlimgioqKlB9FquhzWjd3LoPPPpuWQYNYcd11W+V0+U03ZXUjGng3nC0ZN47mhgZUhOaGBpaMG5fVjWjg3XD2ybXXEttuO1SE2Hbb8cm119r1uXkQJKd+8zgFmK+qExK+9Dhf3BB2BvBY6IV3LuusupRT5/hTiFW/9x6bR4+mdcgQVITWIUNYf9ttxGwFtMh1llMXXlMzOaM7FngWGCoiD+AdWf4gyqLClGyeXLvxzLj+VksOijqn8Sa3tX9/Fv/hD7QMHMi67wV/HV933HFZN7bJrDnmGGtsCyBgTg8GTgfeFZG3/G3XADcD00XkLOBT4DtBniQHRZ1VpySZJ9duPMs/119PM5l1YaaIzAUOAAS4WFW3vWCt8zFm493Zmle2GIRJphjeaslWMec0WZNrTNCcquoreFlI5sicBw4oaFYLlVPn2GIQTiiG19NMZl041P/jBv/z7iKCqr4UXVnBWZNrOuP6EWi2ijWn1uSazpRaTqF4s+oUa3Kd4npOM+n+EufhqQO+hnfH6BGRVBQCa3JNOq4fgeag6HJqTa5JpwRzCkWYVadYk+sc13OayaULW90BIiJDgVsiqyigvlOn0u/WW63JNSkFXZZQRO4F4tP9DPe3jQPOAeLzWV2jqk8HLDVjxZZTa3JNOq4sHxq2YsuqU2Ixep5zjjW5DimGnObSBS4GhoddSBj6Tp1KgzW5JgMBj0DvA+4E/tBh++2qemuQgUPkbE6tyTWZcv1MUUiczapTYjF6WZPrJNdzmsk1ur/BnwYFbzqyvYG3oyimsrKSnj175vS93SdNotett7L5uONYfuutJdnkpprnMxfr1q0LZZwVK8KZgvKzzz4LZZxMBbzJ5SV/pSVn5Dun3bunXJypU9VvvEGfc86hbeBAFt93XyhNbli52LBhQ/qdMrAqyeISuVoW0qqNixYtCmWcDz/8MJRxMuX6C2gu8pXVIK+nicL6O6hMsux2xmIxup97LjUzZ7LptttoPfNMqgPWE1bewxxr48Zw5owIK++Zcj2nmXSDiUsstgAPqurfIqonJ90nTaLXTTex+bjjWP2b30Bra6FLMg6L8K2WC0Tk+3iZuVxV10TxJCk4n9PqN96gz3e/S9uAAax+6CFaunUrdEnGYcXwlmiOnM+qU2Ixup95JjXPP8/GX/6SLWeeWeiKTIJiyGkm1+hOzUchudqmya2qskbXpJXmCLSfiCS+GE3yl9bszN3Az/DO1PwMuA34YaAis+B6Tjs2uW0NDRDiGRVTmlw/U5QL17PqlA5NbuzMMzNa5crkl+s5zeTShXf54m2Wrb4EqKruGXpVGUra5BqTgTRHoKtUdUQ246nq8oSx76GTdeuj4HJOkza5xmTA9TNFuXA5q05J0uQaN7me00w6w2f8z/f7n78LbAYKelRqTa7JVRQTXItIg6rGL6g8CXgv1CdIz8mcWpNrclUME9HnyMmsOsWa3KJRDDnNpDs8WFUPTnh8lYj8TVVvjKqoZLo8+ig9b7mFyqVLaevRg8p166zJNTkLOL3Yg8BheJc4LMZb0vMwEdkb70zNIuBHgYvMjhM5rXvkEbr/7/96Oe3XD1m/nrbBg63JNTlx/UxRjpzIqktqZsygy003UbFkCW2DBqG9e1P13nvW5BaJMHIqIpV4168vUdXjRaQP8GdgB7zX1FNyve8lkw6xq4gc4i+niIgcBHTN5cly1eXRR+l91VVUNDYCULluHVpZSePIkdbkmpwEnHVhTJLNU3KvJhQFz2ndI4/Q84orvsjpypWoCJt++ENrck1OXD9TlKOCZ9UlNTNm0O3SS5H4/xtLlsCSJTSNGWNNbpEIKacXA/OBHv7jq4BZqnqziFzlP74yp/oy2OcsYKKILBKRj4G7yONNNgA9b7ml/cUzTlpb6XmrK1OWmmISf6sl1UeRKnhOu//v/26bU1W63n13PsswJSJoTkXkXhFZISLvJWwbJyJLROQt/+PYSH+I5AqeVZd0uemm9iY3UfXLLxegGpOtdDnNMKtDgOOAyQmbT+CLy3mmAifmWmMmsy7MAfYSkR6AqGo4E7BmoXLp0qy2G5NOqb0lajk1pShgTu/DwYVdXMiqSyqWLMlqu3FPCK+ndwA/ARInaB8Yv+9FVZeJyIBcB0/baovIQBGZAvxZVdeJyO4iclauT5iL1kGDstpuTDqldkbXcmpKUZCcqupLwH+irzI7LmTVJW0p/n9oGzw4z5WYXGVwRrefiLyZ8HFu/HtF5HhghX8AGE19GexzH/AcEP/X+CFwSUT1JLXuJz+hrbZ2q21t9fWs+8lP8lmGKSHxSa6TfRSp+yhwTjdcfTVt9fVbbWurr2fD1VfnswxTQiLK6QUi8o5/aUPvsGrNwn0UOKvOiMXQ3tv+FWh9PZuvu64ABZlcdJZTP6urVHVEwkfivPQHA98UkUXANOAIEfkjsFxEGvzxG4Ccl2HNpNHtp6rTgTYAVW0B8roiw+aTTmpvahVoGTyYNTffzOaTTspnGaZElOg1ugXPadO3vsW6X/4Sraxsz+m6X/6Spm99K59lmBKRQU5TniXqxN3AznjL7i7DW9gl3wqeVSf4U4hVvfceTWPG0DpkCCpC65AhbLz9dppPPrnQFZoMBL1GV1WvVtUhqroDMBr4i6p+D3gcOMPf7QzgsVxrzGTKgk0i0hd/gmsROQDI+zVFjUcfTe+f/Yz/TJjAZguACaiIz9ym4kROm771LVpvu40te+/N2okT8/30psSU2sIuPieyWlA2T25Jiej19GZgun9Zz6fAd3IdKJNG9zK8znpnEfkb0B+wTtMUrWKY4DoHllNTUkp0YRco96xak1tSwsypqs4GZvt/Xg0cGca4mcy6MFdEvg7sirdE4QequiWMJzemUEqt0bWcmlIUJKeOLuxS3lm1Jrckuf56mrLRFZH9gc9U9XNVbRGR/YBvA5+IyDhVde5uVmMyVSqXLlhOTSkLklPXFnYp+6xak1uyXH897awN/x3QDCAih+JdL/EHvGuJJnXyfcY4rcRuRrOcmpJUYjmFcs6qNbklK4wFI6LW2aULlQlHmKcCk1T1YeBhEXkr8sqMiZDrR6BZsJyaklVCOYVyzWosRs+zzqJm5kxrckuU6znttNEVkSp/6pMjgcSpWzK5iS1rFRUV1NXVJf+aP49udXV1yn1cFYvFCl2C6cCFo8yQ5D2nIkJNTU3qr+P9fjvbB0j79UzVdphjO1fNzc2hjBNWPRDe78i1cTJVQjmFPGe1oqIilH+LlZWVuX9zLEb9uedSPXMmTXfcAWedRdCK1q0LZ4KKzZs3hzIOQGOSJYxz8emnnzo1TqZcz2ln4XoQ+KuIrAIagZcBROTLlNtUKKbkuH4EmgXLqSlZJZRTKLesxmLUf//7VD/7LI0TJtB6Vtku/lbyXM9pykZXVceLyCygAXheVdX/UgVwYT6KMyYKpTS9mOXUlKpSyimUWVY7NLlbzjoro9WpTPEphpx2+naJqr6WZNuH0ZVjTH4EnLboXiC+Pvdwf1sf4M/ADnjTFp2iqmsCF5oBy6kpVa6/gGarLLKapMk1pc31nLpdnTERSbMudzr3AaM6bLsKmKWqw4BZ/mNjTAABc2ryzZrcstRZTl3IamSNrogMFZEXRWS+iMwTkYujei5jshF02iJVfQnoOOflCcBU/89TgRNDLToillPjqhKcXixnRZFTa3LLUjFMLxZlBS3A5ar6FeAA4HwR2T3XwWqefRaA7hdeSJ9996X24YfDqdKUpQiOPgfGlxb1Pw8IrdhohZvTGTOo+PRTah5+mF57703NjBmhFWrKj8tnifIs1JyGpWr6dLoNH073Xr3ovv321uSWKdfP6EYy/RC0v9jHX/g3iMh8YDDwfrZj1T78MN1+/nPAm7qocvFiul92GQCxb387rJJNGUlzlNlPRN5MeDxJVUtyQvcwc1ozYwbdLr0UaW0FvJx2u/RSNgLNJ58cYtWmXLhwNsgFYeY0LFXTp1N/0UVIfGqtpia0uhrt3r1QJZkCcT2nealORHYA9gFez+X7u44fjzQ1bT1mYyNdx48PXpwpOxkcfa5S1REJH5k0uctFpMEfvwFYEeXPEIWgOe1y001fvOjFx2xspMtNNwUvzpQd188SFUrQnIal7sYbt837li3U3XhjgSoyhZAupy5kNfJGV0S6AQ8Dl6jq+iRfP1dE3hSRN1evXp10jIolS7Labkw6EVxP9Dhwhv/nM4DHQik0TyynxkVBcioi94rIChF5L2FbHxGZKSIL/c+9I/0BQpZNTletWhVtLYsXZ7XdlK5yvkYXEanGC+UDqvpIsn1UdVL8zFnfvn2TjtM2eHBW241JJ+AL6IPAq8CuIrJYRM7CW7f+KBFZCBzlPy4KllPjqoAvnvdRQrOjZJvTfv36RVdMLAYpVl3TIUOie17jpLJtdMU7Xz0FmK+qE4KMtenaa9EOy/5qfT2brr02yLCmTAV9m0VVx6hqg6pWq+oQVZ2iqqtV9UhVHeZ/7jgrg5PCzOnm665D6+u32qb19Wy+7rogw5oyFUJOS2l2lNByGpg/u4L41+Qm0vp6mq6/vkCFmUIo90sXDgZOB44Qkbf8j2NzGSj27W+z8ZprAFCgdcgQNkyYYDeimZy5fPSZZ6HltPnkk9l4++1oZWV7TjfefrvdiGZyFkFOi3V2lNByGkiHKcQa77qLtqFDURHahg6l8de/puWUU/Jeliks18/oRjnrwit4kySEonnUKLj+ejb85jfETj01rGFNmXLhKNMFoef05JNp+8UvaNl3Xzb+7ndhDWvKVJqcltPsKKHmNCcp5sndaI1t2XP99TSyRtcYV0kRrM1tTLnLIKerVHVElsMuF5EGVV1WrLOjFIQtBmFSKIbXU7erMyYiLl9PZIzxRJDTop4dpSCsyTVpuH6Nrp3RNWXJ9SNQY0ywnPqzoxyGd4nDYmAs3mwo0/2ZUj4FvhNCmaUrFqP+zDOtyTWdcv311BpdU3aK4a0WY8pd0Jyq6pgUXzoy50HLSSxG1x/8gOrnnrMm16RUDK+nTjW6IkJtirn5pKYGgOrq6pTz94UtFovl5XlM/rnwdkqxqqioSJlTAKmooKKystN9AJqbm0OpJ6xxavz/Y1wZB0j7O8z3OGH+bJmwnOZORKjrMC1nxmIxan74Qyqfe44tv/kNFeecQ9B/QevWrQs4gmfz5s1OjQOwYMGCkhwnU0FyKiJDgT8A2wFteDeV/kpE+gB/BnYAFgGnqOqaXJ7D7TbcmIi4PBWKMcZjOS2AWIya006j8plnaP7Vr2g755xCV2QcF3B6sRbgclX9CnAAcL6I7E6Ii7s4dUbXmHyxM0XGuM9ymmcdmtzWc86xs2EmrSA59eezjs9tvUFE5gOD8RZ3OczfbSowG7gyl+ewRteUnWK4psiYcmc5zbMkTa4x6YSZUxHZAdgHeJ0Oi7uISM6Lu1ija8qSnSkyxn2W0zyxJtcEkEFO0y7uIiLdgIeBS1R1fZjZt0bXlCU7U2SM+yyneWBNrgkog5x2uriLiFTjNbkPqOoj/ubQFncpmv9Fqp5+GoC6//5vug0fTtX06QWuyBSr+FstdpNL+KqmT0c++YSq6dPpsvvullOTM8tpdCqnTaN2112p69qVuoYGa3JNztLlNF1WxTt1OwWYr6oTEr4U2uIuRXFGt2r6dOpuvBHwFvuWzz6j/qKLaARabJ1tkwN7SzR8VdOnU3vBBUhLC+DltPaCCwDLqcmN5TR8ldOmUX3++Uhjo7ehqQmtrobu3QtbmClaAXN6MHA68K6IvOVvu4YQF3cpisPiuhtvRJqattomjY3tza8x2bIzReGrGTfuixdPnzQ2UjNuXGEKMkXPchq+qrFjt83pli1UjR1boIpMsQtyRldVX1FVUdU9VXVv/+NpVV2tqkeq6jD/839yra8ozujK4sVZbTemM2Gsvy0ii4ANQCvQ0tn1R+XCcmrCFEZOzbYspyZMxZDTomh0dcgQ5LPPkm43JhchnRE6XFVXhTFQKbCcmrDZmdvwWU5N2FzPqdvV+Zquvx7tsJSh1tfTdP31BarIFLv4UWiyD5Ob5nHj0Pr6rbZpfT3NdumCyZHlNHwtN9yQNKctN9xQoIpMsesspy5ktSga3ZZTTmlvahVoGzqUxl//2m5wMTkJ6W5uBZ4XkTkicm6E5RaNllNOIXbnnWhVVXtOY3feaTk1OQkjpyKySETeFZG3OszjWbZaR49my8SJaLdu7TndMnEiraNHF7o0U4SCzrqQD0Vx6QJAy7HHwjXX0PTb37JlzJhCl2OKXJrwpZ3cGjhYVZf6q7XMFJEFqvpS6IUWmZZTTqHm5z+ndb/9iE2ZUuhyTJGzS4yi0Tp6NBWvvUblI48Q++CDQpdjipwLzWxniqbRNSZMad5O6XRyawBVXep/XiEijwJfA8q+0TUmTC687WmM6ZzrOXW7DTcmAkHfEhWRriLSPf5n4GjgvYjLNqasZJDTfiLyZsJHskuI7BIjYyJkly4Y46iAR6ADgUf9MaqAP6nqs2HUZYz5QtB3XrBLjIyJnOtndK3RNWUpyFGmqn4E7BVeNcaYZIKeDbJLjIyJngtnbTvjVKNbUVFBbW1t8i/W1ABQVVVFRap9fLFYLOzSTAmJv9ViciMi1Ph5TPp1vCx3tg+Q9uuZCmuclP/3FGgcKO2fLZ2gOfUvK6pQ1Q0JlxiVzXKanb6eAlJZCaT/O123bl0o9WzevNmpcT4I8Sa8BQsWODVOmD9bOsXweupUo2tMvrj+Vosxxi4xMqYYuP56ao2uKUuuH4EaY+wSI2OKgeuvp9bomrLk+hGoMcZyakwxcD2n1uiaslMM1xQZU+4sp8a4rxhyGml1IjJKRD4QkX+JyFVRPpcx2XB5Xe58s5waV1lOv2A5Na7qLKcuZDWyRldEKoGJwDHA7sAYEdk91/EqnnwSgKqzz6Zml12omDYtlDpNeXJ5cut8Cj2nDz4In3xCxbRpVA8b5j02JkeWU0/YOeVPf4L774dVq5Add/QeG5Ojcl4w4mvAv/wbAhCRacAJwPvZDlQxbRpV118PeFMX8emnVJ13Hi1A2+jRoRVsykMxvNWSR+Hl9MEHqTzvPKSlxdvw6adUnnceAG1jxoRVrykTltOthJZT/vQn5Ec/QuLTdH36KfzoRyjAaaeFVa8pE8WQ0yirGwx8lvB4sb8ta1XXX480NW21TTZvbm9+jcmWy2+z5FloOa28/vovXjx9snkzlZZTkyPLabvQcirXXps0p3LttblXZ8qa65cuRHlGN9lPp9vs5K0/fi7A0KFDk4/02WfZbTemE8VwBJpHWed0++23Tz6S5dSEyHK6FcupcVIx5DTK6hYDiZ3rEGBpx51UdZKqjlDVEf37908+UqoGONV2Y9Jw+egzz7LOab9+/ZKPZDk1IbOctrPXU+Ms18/oRtnovgEME5EdRaQGGA08nstALTfeiHbpstU27dKFlhvLZjVHEzKXL5zPs9By2poip62WU5Mjy2m70HKq48cnzamOHx+8SlOWXL8ZLbIKVLUFuAB4DpgPTFfVebmM1TZ6NC133YVuvz0qgm6/PS133WU3opmcuXz0mU+h5nTMGFo75LT1rrvsRjSTM8upJ8ycctpp6O9+t1VO9Xe/sxvRTM5cP6Mb6YIRqvo08HQYY7WNHk2zNbYmBMVwTVE+hZrTMWOssTWhsJxuLcycctppqDW2JgTFkFO3qzMmIkHfZrHJ242JnuXUGPcFvXQh6pxao2vKUpC3WUKfvN0Yk5Tl1Bj3Bbl0IR85jfTSBWNcFMJbLeFN3m6MScpyaoz7iiGndkbXlKWAF86HNnm7MSY1y6kx7gt4M1rkOXXqjO7cuXNX1dXVfZJmt37AqnzUkyGrp3P5rOdLmew0Z86c5yoqKlJMBgtAnYi8mfB4kqpOSnic0eTtpWru3LmramtrLafBlHs9abNqOQ1mzpw5qyorKy2nwblWk1OvqRnkFDrPauQ5darRVdUUM1x/QUTeVNUR+agnE1ZP51yrB0BVRwUcIqPJ20uV5TQ4qyc9y2kwltNwuFaTa/UUQ07t0gVjshfa5O3GmMhYTo1xX+Q5deqMrjHFQFVbRCQ+eXslcG/Ok7cbYyJhOTXGffnIaTE2upPS75JXVk/nXKsnFKFO3l6aXPt7t3o651o9obCcpuXa37tr9YB7NblWT2BR51RUy+bafGOMMcYYU0bsGl1jjDHGGFOSiqbRdWkpRxEZKiIvish8EZknIhcXsp44EakUkX+KyJOFrgVARHqJyAwRWeD/rg4sdE0mWi7l1K/HuaxaTo0LXMqqizkFt7JqOc1dUVy64C8R9yFwFN5UFG8AY1S1ICvciEgD0KCqc0WkOzAHOLFQ9STUdRkwAuihqscXsha/nqnAy6o62b+bsouqri1wWSYiruXUr8m5rFpOTaG5llUXc+rX5UxWLae5K5Yzuu1LxKlqMxBfIq4gVHWZqs71/7wBmE+BV9wRkSHAccDkQtYRJyI9gEOBKQCq2myhLHlO5RTcy6rl1DjCqay6llNwK6uW02CKpdF1dilHEdkB2Ad4vcCl3AH8BGgrcB1xOwErgd/7b/1MFpGuhS7KRMrZnIIzWb0Dy6kpPGez6khOwa2sWk4DKJZG18mlHEWkG/AwcImqri9gHccDK1R1TqFqSKIK2Be4W1X3ATYBBb9m00TKyZyCG1m1nBqHOJlVF3Lq1+FaVi2nARRLo+vcUo4iUo0XyAdU9ZFC1gIcDHxTRBbhvQV1hIj8sbAlsRhYrKrxo/IZeEE1pcu5nIJTWbWcGlc4l1WHcgruZdVyGkCxNLpOLeUoIoJ3rcx8VZ1QqDriVPVqVR2iqjvg/W7+oqrfK3BNnwOficiu/qYjgYLeWGAi51ROwa2sWk6NQ5zKqks5BfeyajkNpihWRnNwKceDgdOBd0XkLX/bNf7qHuYLFwIP+P+RfgScWeB6TIQczClYVjNhOS0zDmbVcpqe5TRHRTG9mDHGGGOMMdkqlksXjDHGGGOMyYo1usYYY4wxpiRZo2uMMcYYY0qSNbrGGGOMMaYkWaNrjDHGGGNKkjW6AYjItSIyT0TeEZG3ROT/FbqmXInIuSLy54THPUTk3yKyYyHrMiYoy6kx7rOcmqhYo5sjETkQOB7YV1X3BEay9drhxeYeYIiIjPQf34g3t+LHBazJmEAsp8a4z3JqomSNbu4agFWqGgNQ1VWquhRARPYTkb+KyBwReU5EGhK2vy0ir4rIL0XkPX/7D0TkzvjAIvKkiBzm//lof/+5IvKQvxY4IrJIRG7wt78rIrv527uJyO/9be+IyLc7GydOvQmV/we4Q0RG4K288ssIf3/G5IPl1Bj3WU5NZKzRzd3zwFAR+VBE7hKRr0P7et2/AU5W1f2Ae4Hx/vf8HrhIVQ/M5AlEpB9wHTBSVfcF3gQuS9hllb/9buDH/rafAutU9av+kfFfMhgHAFV9B2+lnFl+nc2Z/jKMcZTl1Bj3WU5NZIpiCWAXqepGEdkP+C/gcODPInIV3j/64cBMEQFvecVlItIT6KWqf/WHuB84Js3THADsDvzNH6sGeDXh64/4n+cA3/L/PBJvbe54nWtE5Pg04ySaCByjqi+mqc0Y51lOjXGf5dREyRrdAFS1FZgNzBaRd4Ez8EIyr+NRpoj0AlKtt9zC1mfX6+LfBsxU1TEpvi/mf27li79LSfI86cZJ1OZ/GFMSLKfGuM9yaqJily7kSER2FZFhCZv2Bj4BPgD6i3dxPSJSLSJ7qOpaYJ2IHOLv/92E710E7C0iFSIyFPiav/014GAR+bI/VhcR2SVNac8DFyTU2TvHcYwpepZTY9xnOTVRskY3d92AqSLyvoi8g/dWxjj/OpyTgV+IyNvAW8BB/vecCUwUkVeBxoSx/gZ8DLwL3ArMBVDVlcAPgAf953gN2C1NXTcBvUXkPf/5D89xHGNKgeXUGPdZTk1kxLs50OSbiOwAPKmqwwtdizEmOcupMe6znJrO2BldY4wxxhhTkuyMrjHGGGOMKUl2RtcYY4wxxpQka3SNMcYYY0xJskbXGGOMMcaUJGt0jTHGGGNMSbJG1xhjjDHGlCRrdI0xxhhjTEn6/wE8Nws2/gLfHwAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 720x216 with 6 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "X = [1, 1, 1, 1, 2, 1, 1, 6, 6]\n",
    "Y = [0, 3, 3, 3, 9, 9, 7, 7]\n",
    "\n",
    "print(r'Accumulated cost matrix and optimal warping path for different local weights:')\n",
    "plt.figure(figsize=(10, 3))\n",
    "\n",
    "ax = plt.subplot(1, 3, 1)\n",
    "weights_mul = np.array([1, 1, 1])\n",
    "title='Weights: '+'[ '+''.join(str(s)+' ' for s in weights_mul)+']'\n",
    "compute_plot_D_P(X, Y, ax=ax, weights_mul=weights_mul, title=title)\n",
    "\n",
    "ax = plt.subplot(1, 3, 2)\n",
    "weights_mul = np.array([2, 1, 1])\n",
    "title='Weights: '+'[ '+''.join(str(s)+' ' for s in weights_mul)+']'\n",
    "compute_plot_D_P(X, Y, ax=ax, weights_mul=weights_mul, title=title)\n",
    "\n",
    "ax = plt.subplot(1, 3, 3)\n",
    "weights_mul = np.array([1, 3, 3])\n",
    "title='Weights: '+'[ '+''.join(str(s)+' ' for s in weights_mul)+']'                                                      \n",
    "compute_plot_D_P(X, Y, ax=ax, weights_mul=weights_mul, title=title)\n",
    "\n",
    "plt.tight_layout()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Global Constraints\n",
    "\n",
    "One common DTW variant is to impose **global constraints** on the admissible warping paths. Such constraints not only speed up DTW computations but also prevent \"pathological\" alignments\n",
    "by globally controlling the overall course of a warping path. More precisely, let $R\\subseteq [1:N]\\times[1:M]$ be a subset referred to as a global **constraint region**. Then a **warping path relative to $R$** is a warping path that entirely runs within the region $R$. The **optimal warping path relative to $R$**, denoted by $P^\\ast_{R}$, is the cost-minimizing warping path among all warping paths relative to $R$. The following figure shows two global constraint regions known as the **Sakoe&ndash;Chiba band** and the **Itakura parallelogram**. Alignments of cells can be selected only from the respective shaded region.\n",
    "\n",
    "<img src=\"../data/C3/FMP_C3_F18_text.png\" width=\"500px\" align=\"center\" alt=\"C1\">\n",
    "\n",
    "For a general constraint region $R$, the path $P^\\ast_{R}$ can be computed similarly to the unconstrained case by formally setting $\\mathbf{C}(n,m):=\\infty$ for all $(n,m)\\in [1:N]\\times[1:M]\\setminus R$. Therefore, in the computation of $P^\\ast_{R}$ only the cells that lie in $R$ need to be evaluated. This may significantly speed up the DTW computation. However, the usage of global constraint regions is also problematic, since the unconstrained optimal warping path $P^\\ast$ may traverse cells outside the specified constraint region. In this case, the constrained optimal warping path $P^\\ast_{R}$  does not coincide with $P^\\ast$ (see the last case of the above figure). In the following code cell, we call the function `librosa.sequence.dtw` with different constraint regions determined by a a Sakoe-Chiba band. \n",
    "\n",
    "<!--**Note**: In the following example, the warping path does not begin at $(1, 1)$. The reason for this is a bug in the backtracking procedure in `librosa`. This was [fixed](https://github.com/librosa/librosa/pull/920) and will be solved in the next version of `librosa`. However, currently the fix is not included the the latest version.-->"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2024-02-15T08:50:54.343494Z",
     "iopub.status.busy": "2024-02-15T08:50:54.343233Z",
     "iopub.status.idle": "2024-02-15T08:50:54.985461Z",
     "shell.execute_reply": "2024-02-15T08:50:54.984909Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Accumulated cost matrix and optimal warping path for different constraint regions:\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAArgAAADgCAYAAADoghWcAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAABJ50lEQVR4nO3deZwU1b3//9dnFmYGGPZFEBLcjRhxIcbthiBLnGhcrho1uUZN1BsTzSIm12hAohiNUfH+jBcvbiFoRBG9bsEEyBjinoGAIhj1q6gosgkjA8MAM5/fH1U9Nk33dPdUdffp6s/z8egHM9Xdp07PzJv6VNWpU6KqGGOMMcYYExVlhe6AMcYYY4wxYbIC1xhjjDHGRIoVuMYYY4wxJlKswDXGGGOMMZFiBa4xxhhjjIkUK3CNMcYYY0ykFKTAFZGVIjK2EOtO6Mf5IvJcB8/PFZHzMnltmOsNUybrEpHnReSwENf5rIhcGFZ7+SQiVSLyhogMKHRfXGKZtcy6yjK7O8ur5dVV+cxrpwtcETlORF4QkUYR+cT/BX4pzM6FQUS+JiILRWSziKwTkb+JyMmZvFdV61R1Rq77WEgi8g1gs6r+s9B96YiIjBaRev/vbWUGrx/jh2ir/77Pxz0nIvIbEdngP24SEQFQ1RbgXuC/cvZhCsQyGw1FlNmUOUvy2mEioiLSFPeYmElbUc2s5TUaIprXo0Rknv93uU5EZovIoLjnJ4vIjoQ87w35zWunClwR6QE8BdwO9AH2BH4FtITXteBE5AxgNvAHYAgwEJgEfKOQ/XLM94GZqZ4UkYo89qUjW/BC8bN0LxSRfsCjwES8v88G4KG4l1wMnAqMAA4BTgL+M+75PwLniUhVGB13gWU2Uools+lylkwvVe3uP67Loq1IZdbyGilRzGtvYDowDPg8sBm4L+E1D8VlubuqvhP3XF7y2tkjuPsDqOqDqtqqqs2q+hdVfRVARPYRkb/6ewHrReQBEemVrCEROVBE3hWRs/3vZ4vIx/5e60IRGR732p4i8gd/j+E9EfmliCT9DP6ex63Adap6t6o2qmqbqv5NVS9KeO3NIrLR70dd3PLE0wAiIrf7fXtDRMbEPXGBiKzw92LfEZF0/5F3qi0R+aqIrBKRCSKyVkRWi8gFcc/3FZEnRORTEXkF2KeDDnQBjgf+Frdssog8IiL3i8inwPkicqSIvCgim/z1/c5/b+w94/zP0CgivwOS7vUFoaqvqOpM4J20L4Z/B15X1dmqug2YDIwQkQP9588DblHVVar6IXALcH7culYBG4GjQvwIhWaZtczmNbOkyVmYbUUws5ZXy6uzeVXVuf729VNV3Qr8Djg20xXlLa+qmvUD6AFsAGYAdUDvhOf3BcYBVUB/YCFwW9zzK4GxwOHA+8BJcc99F6j133sbsCTuuT8Aj/vPDwPeBL6Xoo8HAgrs1cHnOB/YAVwElAOXAB8B4j//LHBh3Gt3Aj8FKoGzgEagj//8iXh/6AKMArYCh3ew3k61BXzVf++1/nu/7j/f239+FvAw0A04GPgQeC5FP4YDWxKWTfZ/Jqfi7QDVAEfg/SFW+D/3FcBP/Nf3Az4FzvD781O/fxemWOe3gE0dPD6X5m9vLLAyzWv+G5iWsGwZcLr/dSPw5bjnRuKdQop//RPAjzqTDxcfWGYts3nOLBnkLO65Yf7v/kNgFd7RoH7ZtEWEMovl1fLqcF6TvPcnwEsJn7ER+AR4HbgkyXtyntcgAfwC8Hu8/4x2+p0dmOK1pwL/TAjfr/z3ju5gHb3wAtQTLxwtwEFxz/8n8GyK9x7rv7c6Tfjejvu+q/+ePVKErz2Y/rJXgHNTtP1/wI87WG+n2sILXzNQEff8Wj8c5XjBOTDuuV+TOnzHAh8nLJsMLMzgj/kx/+vvJPxhi/97TRq+wH+wmRW49wA3Jix7Hjjf/7o14We0n/97j/99PABMysVnKNTDMmuZ9b/OS2YzyVncc93xNqgVeKe5HwH+nE1bUcus5dXy6n/tXF4T3ncIXiH7b3HLDgIG+z+vY4DVwDkJ78t5Xjt9kZmqrlDV81V1CN5ezGC8vUFEZICIzBKRD/1D8Pfj7YXE+z7wgqrWxxaISLmI3Cgi/89/30r/qX7+owvwXlwb7+GNTUJE7pTPBjNfhbf3CzCIjn0c95m2+l92T/HaD9X/zcStf7C//joReUm8Qdeb8Pb6Ej9zWG1tUNWdcd9v9fvcH28D8UFCu6lsxNtTTxT/fkRkfxF5yj+t9SleoGP9GRz/ev8z7fL+AmjCOwISrwfeOKFkz/cAmhJ+H7V4e7uRYZltX79lNj+ZzSRnsT40qWqDqu5U1TXApcB48caiZtpWpDJreW1fv+XVsbzGiMi+wFy8nYO/x/Vxuap+pN7wmhfwzqqekfD2nOc1lGnCVPUNvD3Ng/1FN+BV/oeoag/gP9h9zMj3gc+JyNS4Zd8CTsE7StcT71A9/nvX4+05fT7u9Z/DOz2Aqn5fPxvM/GvgX3h/BKeH8BFj9hTZ5arCzwEfiTdQeg5wM94edi/gT3Q8TibMtmLW4e3pD01oN5W38MYp7ZmwPPEPehrwBrCf//u8Kq4/q+PX53+moaQgIt+WXa+sTHx01N9MvY43UD62zm54p6NeT/a8//Xr7OoLwNIQ+uIky6xlNvbCHGY2k5ylEvs8sT6XdGYtr5bX2Atdyat4MxPNxxuDnfIiOp+y+88353nt7CwKB4o3AHuI//1Q4BzgJf8ltXh7A5v8X2yyK983AycAXxGRG+Pe14K3Z9gVby8GAFVtxRv3cr2I1Po/3Mvx9lx34+91XA5MFG9AeQ8RKRNv6pXpnfncwADgRyJSKSJn4v2C/oS311uF/8cv3iD68XlsC2j/GT0KTBaRriJyEN7A8VSv34H3BzoqTdO1eGOAmsS7UOuSuOeeBoaLyL+LdzXoj4A9OljnA7rrlZWJj/eTvc//3VXjjUESEamWuEH4CR4DDhaR0/33TAJe9TcS4I0zu1xE9hSRwcAEvI1HbF174l25/BIRYZm1zMY9l5fMkiZn8UTkyyJygP/77gv8f3inxhszaStqmbW8Wl7jnnMxr3sCfwXuUNU7kzx/ioj0Fs+Rfp8fT3h/zvPa2SO4m4EvAy+LyBa8Ti7D+4GAN/bncLxBxk/j/UHsRlU34Q2UrxOR6/B+wO/h7TEuZ/cPfxnedFHvAM/hTTVxb6pOquojeIPLv4s3HmcNMIW4H3SWXsYbl7IeuB44Q1U3qOpmvF/gw3inJL6FN14qX23FuxTvVMrHeH+c96V5/f8C56Z5zRV+PzYDdxE35ZaqrgfOBG7E+09zP7zxrmH7Ct64qD/h7TE3A3+JPSkir4vIt/0+rcM7qnA93s/wy8DZcW39L/Ak8Bre3+3T/rKYbwEz1JuvLyoss5ZZIK+Z7TBn8ZkF9gae8fu7DK8IOyfTtoheZi2vllfA2bxeiJfZa+KPDse1dTbwtv95/gD8Rned7zgveY1dyWhKmHh3YblMHZ+IOh/EO321FPiKqq4tdH+MScYy+xnLrHGd5fUz+cyrFbjGGGOMMSZSQrnIzJhS4o//fUVElvqnbX7lL58s3lXNS/zH1wvdV2OMMaYU2RFcY7IkIgJ0U9UmEanEG6v2Y7wLOppU9eaCdtAYY4wpca7cA9mYouFfPRwbUF/pP2xP0RhjjHGEDVEwphPEmzB9Cd4dbuap6sv+U5eKyKsicq+I9C5cD40xxpjS5dQQhZqaGu3RI/EGVNmrrq4OoTfhcelnHLa2trZQ2gnjZ7Rp0ya2bNmSdrJuEUm3steBbXHfT1fVpPM6ikgvvHl3L8ObV3E93tHc64BBqvrdDLpelPr166fDhg0rdDectWbNmkJ3YTcu5RVg9erV61W1f0evySCv4N3W94RQOhVRltfi8957Hd0kLXNh5f6TTz5Jm1fIKLN5yatTQxR69OjBWWedFbidgw8+OP2L8qilJSpTM+5u27Zt6V+UgTB+RtOmTcv4tSKp62BV3aaqIzNpR1U3icizwAnxY29F5C7gqYw7VISGDRtGQ0NDobvhrJtvdm8otkt5BZgyZUpGW/CO8gqgqh3dstVgeS1GF110USjthJXXmTNnZlxxp9nG5iWvNkTBlKSysrKUj3REpL9/5BYRqcG77eUbIhJ/T/bT8CbLNsYE1FFeM8msMSa/XMirU0dwjcmXdEeE0hgEzBCRcrydxIdV9SkRmSkih+INUVgJ/GfQfhpjAufVGJNnLmTWClxTckQk0F6kqr4KHJZkebrbMRpjshQ0r8aY/HIlszntgYj81J8If5mIPCgibl39ZUqWiKR8lCrLq3FVR3m1zFpmjXtcyGvOjuCKyJ7Aj4CDVLVZRB4GzgZ+n6t1mugZvnQpo+fNo2djI409e1I/bhyvjxgRuF0X9i5dYnk1YTlk2TLG1tfTq7GRTT17Mn/0aF4NeOGv5XV3llkThqPfeYczlyyh75YtbOjWjdmHHsqLe+8duF0XMpvrIQoVQI2I7AC6Ah/leH0mQoYvXcqJjz9Olx07AOjV2MiJjz8OEKjILfWjPh2wvJpADlm2jFOefro9s70bGznl6acBOl3kWl47ZJk1nXb0O+/w3Zdeoqq1FYB+W7bw3ZdeAghU5LqS2ZyV2Kr6IXAz8D6wGmhU1b/kan0mekbPm9e+oYzpsmMHo+fNC9y2C1d4uiTMvM6dO5fRo0czceJE6uvrGTVqFBMmTChYO672KYrG1tcnzezY+vpA7dosCruLcmatnfw4c8mS9uI2pqq1lTOXLAnctgt5zdma/Ls4nQLsBQwGuonIfyR53cUi0iAiDc3NzbnqjilCPRsbs1qeDRfC55LO5HXdunVJ21qwYAH19fXU1tYydepU5s+fz4gRI1i+fHlWfQqrHVf7FEW9UmQz1fJMWYG7u0wym0lewb18WDv50XfLlqyWZ8OFvOZyTWOBd1V1naruAB4Fjkl8kapOV9WRqjqypqYmh90xxWZr165Jlzf27BmoXbtgJams89q/f/Ib2sTuctWnTx+amppobW1FVbO++1VY7bjap6ipaW6mNcXGa1OAzKbLq2U2dWYzyav/OsCdfFg7+bG1sjLp8g3dugVq15W85rLAfR84SkS6iveJxgArcrg+EyH7vfEGNc3NtCWEYXtlJfXjxgVu34W9S8eEltexY8cyZswYVqxYweTJkxk/fjwNDQ0MHz68IO242qcoqWlu5vwHHgBVdpaX7/Lc9spK5o8eHaj9IEdwRaRaRF4RkaX+jAO/8pf3EZF5IvKW/2/vQJ3Mv8hm1trJvW+89hrdduygNWEb21JezuxDDw3cfsCbKYWSWcnlnoPfqbOAncA/gQtVNeU94wYOHKh2q97ikotbf+73xhucMWsWH++xB4uPOIJ/+9vfMppFYdq0aXz44Ydpdw8rKiq0ZwdHlD755JNFmd6qN0qyzevIkSPVbv2ZWqncqjdW3PZft44/fvObdG1uzngWhSlTpqTNWrq8QseZ9Yu/bqraJCKVwHPAj4F/Bz5R1RtF5Eqgt6r+V4crckw2mbW8Fp9c3ar3G6+9xplLlvD8Xnvx6qBBnLl0aUazKMycOTOjbWPQbWxYmc3pLAqqeg1wTS7XYaIlvrj943nn0VJTw9IvfSnUdYgjk1C7xvJqspVY3L69zz5A52dMSCZoXtU7itPkf1vpPxRv/OpX/eUzgGeBoipwLbMmW/HF7fRjjkHLynjRz21YXMmsbeWNM5IVt7niwvggY4pZquI2F4KO6RORchFZAqwF5qnqy8BAVV0N4P87IGcfwBgHJCtucyXoNjaMzNqteo0T8lncghuTUBtTrPJZ3EJGee0nIvHn36er6vTYN6raChwqIr2Ax0TErXFsxuRYPotbSJvZDvMK4WTWClxTcPssX84peSxubYiCMZ1XvXUr38pjcZthXtdnMjZQVTeJyLPACcAaERmkqqtFZBDekSJjIiffxW0Gmc0orxAss7aVNwW1z/LlnDJzZt6K2xgbomBM9qq3buXbv/993orbmCBDFESkv38UCBGpwZte6w3gCeA8/2XnAY/n7hMYUxhfX7Ikr8VtTMAhRaFk1qkjuNXV1aHMgDBgQDhDqaI8+8HatYU/WBErbtcOHswfzz03b8Ut2BAFk9rMmTNDaWfgwIGhtAOwZs2a0NrqrPji9oE8FrcQOK+DgBkiUo53UOdhVX1KRF4EHhaR7+FNuXVm8J4aA5dddlmhuwB4xe1pDQ15L27Bjcw6VeCa0hFf3M6+8EJa8lxwBjlSKyLVwEKgCi9Dj6jqNSLSB3gIGAasBL6pqhsDd9aYAosvbh/+1rd4+3Ofy+v6g+RVVV8FDkuyfAPe3LHGRE6suH1pn32YftRReS1uwY3M2mEsk3e7Fbd5voNdbHxQgBs9tADHq+oI4FDgBBE5CrgSWKCq+wEL/O+NKWqJxe07++2X1/Wny6udjTFmV/HF7b2jRhWkuHUhr/Y/g8mrQhe3MUHGB6kn1Rx9M/zlM4BTc9B1Y/Km0MVtTNBpwowpFYUubmNcyKsNUTB540pxC8HH4PpjgxYB+wJ3qOrLIrLLHH0iYvNqmqLlSnELNmbemEy4UtyCG5m1AtfkhUvFbQZTmORljj5jXOVScWvT+hmTnkvFrSuZtQI3Ij733HMcMmsWXTdsYGvfvrx69tm8f9xxWbez3z/+wTFPPUXtxo1s7t2bF046ibc6cavcAxcv5ivPPEOPTZvY2q0bNVu3smbIkIIXtzFpTpPkZY4+U7qGvfACh82eTbcNG9jSty//PPNMVh5zTNbtHLh4Mf82dy49Nm3i0169+HtdHW8cfnin+jR86VJGz5tHz8ZGWsvKEFUeOvfcgha3MTYMwUTFEf/6Fye/9BK9N29mY20tTxx1FIsOOCDrdo58+21Oa2igb1MTzV260HX7dieK2xgXMmsFbgR87rnn+NL06VRs3w5At/Xr+dJ074BjNkXufv/4B2NmzaJyxw4AemzcyJhZswCyKnIPXLyYE+bMaW+n25YttImw9EtfcqK4hWCnT0SkP7DDL25jc/T9hs/m6LsRm1fTpDDshRc4+t572/PafcMGjr73XoCsitwDFy/ma4880p6znps28bVHHgHIusgdvnQpJz7+OF38tira2thZXk7N1q1ZtZMrLhwNMiaoI/71L86pr6dq504A+mzezDn19QBZFblHvv023/n736lqbQWg6/bttIqwbMgQJ4pbcCOzVuBGwCGzZrVvLGMqtm9n5N130/+NN5K+p7m5ebdlBzQ0tG8sYyp37OD4hx9mz3feSdpOqx/UeAf985+7tVOmytF//SuvHXVUh58lH0IY6G7zanZg7ty53HTTTRx33HEcf/zxTJ48mZEjR3LLLbcUpJ2w2wrqsNmzk+b1qPvuY+C//pX0PVuT5PWgxYuT5nX8nDkMeffdpO20+hvERAcvXdpe3Lb3qbWV0fPm8fqIESk/Sz7YhWS551pmo9rOyS+91F7cxlTt3MnZzz7Lvh99tNvrU+X1y2+/3V7cxpSrcuqiRbzsyBkXFzJrBW4EdN2wIenyim3b2PMf/0j6XFtb227LKlPc2KKypYW9X301+cpVd399wsY7psemTcnbKIAge5c2r2bHFixYQH19PTfddBNTp05l/vz5PPjggyxfvpyDDjoo7+2E3VZQ3VLltaWFoYsXJ32uLcmGLlXOKrdvZ79ly5I+t3taPV1StNWzsTHFO/LLhaNBUeZaZqPaTu/Nm5Mur9qxg0OS7JRqku0rsFuRHNO3qSnp8kJwIbNW4EbAjq5d6bJly27Lt/brx1O/+13S9yS7k9l511xDj42735dgc+/ezPjVr5K2s23btt2WXfzrX9MzSTH7aa9eSdsoBBf2LqMq9p9ynz59aGpqorW1FVVN+Z91rtsJu62gmnv1omuSfGzp25fHpk5N+p5kdzK76PrrU+bsrquvTtpOsrwCXHrzzfRKUsw29uyZ9PX5ZnnNLdcyG8l2VNlWWUlNwpkSgI21tVxz3nm7LU+V1xtmzaJfkmJ2Q/fumfcnx1zIbOFLbBPI/k8/TZctW2hL2Fva2aULr559dlZtvXDSSeyorNxl2Y7KSl446aSs2ll4wglJ21l4wglZtZMrrkxCHVVjx45lzJgxrFixgsmTJzN+/HgaGhoYPnx4QdoJu60galevprylZbcjqTu7dOGfZ2Y3ouXvdXVJc/b3urqs+1U/bhzbE9raXllJ/bhxWbcVNrvRQ+65ltnItaPKvz/3HDU7dtCaUPi1VFTwRJZD9x4bOZKW8vJd2ykv57GRGV0bnXOubGOlEEcwUvn85z+vV6c48pCNAQPCmX60JcUpe1fs//TTHDZzJu8fdRQfHn44hzz8cMazKCQ7ggvZz6KQag8zfhaFT3v1YuEJJ3R44UsYP+tp06bx4Ycfpt1t7Nq1q+6///4pn1+6dOmiTGdRKGUjR47UhoaG9C8sMjNnzsxJu7WrVzP+hhuQtjaWf+1rHPDXv2Y8i0KyI7iQ/SwKqfIKu86i0NizJ/XjxqUcfxvW/41TpkxJm7V0eQXLbCaimlcXXXbZZZ994xe3o5cupX7ECN4bMCDjWRQ6ymv8LAobunfnsZEjeWXffZO+Nqy8zpw5M6OcubKNtSEKRSq+uH3pssvQ8nLe/8pXArf71pe+1KlpwRK9cfjhnZ6uKB9cOH1iSkd8cTvvyitpHDKE5d/4RuB2w8zZ6yNGFPyCslQsr6YoJRS3jx53HIh0alqwRK/su2/KgtYFLmTWCtwilKy4NZlzZRJqUxqSFbcmc5ZXU5RSFLelwJXMWoFbZKy4DYcLe5cm+qy4DYfl1RSVEi5uY1zIrBW4RcSK2/C4sHdpos2K2/BYXk3RsOIWcCOzVuAWCStuw+PK6RMTXVbchsfyaoqGKlx+eckXt65k1qkCt6amhoMPPjhwO67PfpCtobNns+/MmbxxyCE8dcop6IcfBm5z0aJFIfQM3nrrrVDaefvttwO30ZTFJNcunD4x4Zo3b14o7eyxxx6B3t911SpG/uY37Nixg/svuID1ra3w3nuB2nzllVcCvT/mXynulJatN998M5R2MhUkryIyFPgDsAfQBkxX1f8WkcnARcA6/6VXqeqfAnbVFKkpU6YEa0CVcXPn8uUXX+SJvfbi7qFD4f33AzVZrHkFNzLrVIFrdjd09mz2nTaNtaNG8dQJJ9iR25C4sHdpoqfrqlWMvOIKaGvzituBAwvdpUgImNedwARVXSwitcAiEYntDU1V1ZsDd9CUtrji9uWjj+buvn1L8shtPBcya1t5h8UXt8t/+UsrbkMUu1d2socxnRFf3C666SYrbkPUUV7TZVZVV6vqYv/rzcAKYM88dNuUgoTidl5dXckXtxBsGxtWZq3AdZQVt7njyl1WTHQkFrdbhg0rdJciI8M7mfUTkYa4x8Up2hoGHAa87C+6VEReFZF7RaR3fj6RiQwrbpPKYBubUV79tobRycza1txBVtzmnh3BNWGx4jb3MjiCu15VR8Y9pidpozswB/iJqn4KTAP2AQ4FVgO35O8TmaJnxW2HgubVbyNQZnNa4IpILxF5RETeEJEVInJ0LtcXBVbc5ocdwd2d5TV7VtzmRwZHcDskIpV4G8oHVPVRAFVdo6qtqtoG3AUcmdMPkQOW2QKx4jatoNvYMDKb64vM/ht4RlXPEJEuQNccr6/oDJg/n73vvpvqdevY2a0blU1NVtzmmCtTmDjI8prGHgsWsO9991G9bh0tffpQvm0bbZWVVtzmUNC8infI6B5ghareGrd8kKqu9r89DVgWqKOFYZnNk+FLlzJ63jx6NjayvUsXqrZvt+I2BVcym7MCV0R6AF8BzgdQ1e3A9lytrxgNmD+fA2+5hXJ/WrPKpia0rIz1Rx9txW2OuTCFiUssr+ntsWABB912W3teqzdsQIG3vvc9K25zLODQoWOBc4HXRGSJv+wq4BwRORRQYCXwn0FWkm+W2fwZvnQpJz7+OF127ACgavt2WsvK+GjwYCtuU3Ahs7k8jLU33ob+PhH5p4jcLSLdcri+orP33Xe3byxjpK2Nve69t0A9Kh0BT5/EpjD5AnAU8EMROch/bqqqHuo/iqK49YWW17lz5zJ69GgmTpxIfX09o0aNYsKECQVrJyz73nff7nkFhj75ZGE6VEKCDFFQ1edUVVT1kPhsquq5qvpFf/nJcUeGikVkM+taO6PnzWsvbmPK29oYPX9+1m2ViiDb2LAym8sCtwI4HJimqocBW4ArE18kIhfHrqTbtGlTDrvjnup167JabsIRZMohiOy0Q1nndV2Kv9MFCxZQX19PbW0tU6dOZf78+YwYMYLly5dn1aGw2gmL5bUw0uW1hC8MTZvZTPIK7mXWtXZ6NjZmtbzUuZLXXBa4q4BVqhqb2uERvDDuQlWnx66k69WrVw67456d3ZLvbG/r3z/PPSk9YV1kFqFph7LOa/8Uf6eqCkCfPn1oamqitbUVVW1fnqmw2glLS58+SZdbXnMv6EVmEZU2s5nk1X8d4E5mnWpHle1duiR9qrFnz6z6U0pcyGvO1qSqHwMfiMgB/qIxQGEOvTho6MMPt4+5jddaVcU7F15YoF6VjjR7l5nOqRmZaYfCzOvYsWMZM2YMK1asYPLkyYwfP56GhgaGDx9ekHbC0HXVKsq3bSNxs9haVcXbF1yQ9/6UGleOCLkkypl1ph1/toTYmNt42ysrqR83Lqv+lBIX8prrWRQuAx7wr+58B7AtAV5xu++dd7J21CjWH300e917L9Xr1rGtf3/eufBC1o4dW+guRloGV3iuV9WRadpIOoVJ3PN3AU+F0N18CiWvdXV11NXVtX+/cOHCTnUmrHaCik0F1lZZyVvf+x5Dn3yyPa9vX3ABH48ZU5B+lQqb9aRDkcysE+0kTAX20eDBjJ4/n56NjTT27En9uHG8PmJEp/oVda5kNqcFrqouATosFEpNfHEbmwpszfjxhe5WyXFhChPXWF53l2ye2/fOOqvQ3So5LmwsXWSZzZEU89y+fuihhe5Z0XAhs7k+gmviJCtuTWG4MIWJcZvdxMEdpToMwRSA3cQhFC5k1grcPLHi1h1BT5+o6nN4M0QlKqZpwUwHrLh1hyunO00JsOI2FK5k1grcPLDi1j0u7F0aN1lx6x7Lq8k5K25D5UJmU5bYIjJBRHarxESkr4jck9tuRYcVt25yYQqTsFlmg7Pi1k1RnCbM8uoQK25D50JeOzqCewCwSER+qKrPA4jID4Cf4d3/OnRlZWVUVVUFbieMNsIw4P77GXrnnaw//njevuYaulQEO2BeXV0dUs/Cayusn3UY7WzdujWj10V4aqG8ZzYMr7zySijt9Aw4J2X1++9z4M9/Dqq8ettt7NxrL4L8VbqYV9dy35JwZ7hkLK+mI4FnV1Fl7zvu4Msvvsii447j7yefTHXAvzfX8hrW3OGZ5tCVzKasuFT1YhE5BvidiLwOHAi8BRxThLc0zLsB99/P0Ntu45OxY3lz4kQIWNyacBXrUZ+OWGY7r/r99znwkkuQ1laWTp3K1r32KnSXTBzLq8kJv7gdMns2i447jvqTT7YjtyFxIbPpqq5lwD+AE/AuqplgwUsvvrh9d8oUaG0tdJdMAhfClyOW2SzFF7dv/M//sHXw4EJ3ySSwvJrQxRW3q844g/ojj7TiNkQuZLajMbj/ASzBmzx6H7x5PW8SkT+IyID8dK/47Fbc2pFb50T1rkiW2ewlFrfNe+9d6C6ZBOnyWqyZtbwWUEJx+86ll1pxGyJX8tpR9XUmMFpV3/O/XyQiRwPfB14CbEuQwIrb4uHC3mUOWGazYMVt8bC8mtBYcZsXLmS2ozG4pyRZpsA0EXkkp70qQlbcFpdiPerTEcts5qy4LS6WVxMKK27zxoXMdqoKU9V1YXek2PSeO5c977iDLmvW0Nq9OxWbN1txWyRcmYQ6n0o9s32eeYah06bRZc0advTtS9m2bWhlpRW3RSBoXkVkKPAHYA+gDZiuqv8tIn2Ah4BheHce/Kaqbgzc4RCUel7D1H/ePPa66y6q1q6ltaaGiq1brbjNMVcyW1pb+ZD0njuXz19/PVUff4yoUrF5M1pWxqbjjrPitki4MD7I5EefZ55hrxtuaM9rl/XrKW9qYvW3vmXFbZEIOKZvJ97FW18AjgJ+KCIHAVcCC1R1P2CB/72JkP7z5rH/b39L9Zo13rZ661baysvZfMABVtzmWMBtbCiZtQK3E/a84w7Kt23bZZm0tbHnnXcWqEcmWy5MQh1Vc+fOZfTo0UycOJH6+npGjRrFhAkTCtafodOm7Z5XYOCcOYXpkMlakBs9qOpqVV3sf70ZWAHsCZwCzPBfNgM4NXefwG1hZda1dva66y7KE+ZaLmttZa+77866rSjK5f/VQbaxYWU27ZpEpKuITBSRu/zv9xORk9L2MMK6rFmT1XLjltjpk6gWuIXO7IIFC6ivr6e2tpapU6cyf/58RowYwfLly/PVhV1YXotburxmk1kRGQYcBrwMDIxNyeX/W5CZCwqdVwgvs661U7V2bVbLS02u/q8OcxsbJLOZrOk+oAU42v9+FTAlqx5GTGv37kmXbx84MM89MZ0V8SEKBc1s7K45ffr0oampidbWVlQ1tLvpZGtH375Jl1tei0cGQxT6iUhD3OPiJG10B+YAP1HVT/P9GTpQ8G1sWJl1qh1VWmtqkj7VMsBmYYPc/l8dNK9+G4Eym0mBu4+q3gTsAFDVZrwzfCVpwP33t4+5jddaXc2HP/xhgXplshXlI7gUOLNjx45lzJgxrFixgsmTJzN+/HgaGhoYPnx4vrrQrvr9970LyhKWt1ZX88Ell+S9P6ZzMjiCu15VR8Y9pse/X0Qq8TaUD6jqo/7iNSIyyH9+EFCow3oF38aGlVln2vFnS4iNuY3XWlXFuxddlFV/oiqX/1cHySuEk9lMrojaLiI14G0jRGQfvL3NkhM/Fdim445jzzvvpMuaNWwfOJAPf/hDNtbVFbqLJkNBjtQWwVXZBc1sXV0ddXFZCHyv+E6KTQWmlZV88MMfMnDOnPa8fnDJJXxywgkF6ZfJXsC8CnAPsEJVb4176gngPOBG/9/Hg/QxgIJvY8PKrBPtJEwFtvmAA9jr7rupWruWlgEDePeii1g3blyn+hU1ufy/2oXMZlLgXgM8AwwVkQeAY4HzO9HnopZsntuNJ5X0UOSiFcI0YbErPBeLSC3eBO3z8HKxQFVvFJEr8a7w/K/AHc5eyWc22Ty3H3/nO4XulumEEPJ6LHAu8JqILPGXXYW3kXxYRL4HvI9344VCKPm8hibFPLfrxo8vdM9KiiuZTVvgquo8EVmMN1WDAD9W1fUBOl507CYO0RNk79If3B4b6L5ZROKv8Pyq/7IZwLMUoMAt9czaTRyiJ2BenyP1Kf8xnW44JKWe19DYTRyc4kJmM5lF4TRgp6o+rapPATtF5NRMV1DsrLiNJheu8MyVUs6sFbfRFNYsCi4q5byGxopb57iQ10zWdI2qNsa+UdVNeKdUIs+K22jKYAqTvFzhmUMlmVkrbqMpzGnCHFWSeQ2NFbfOcWUqzkwqtmS9iXylZ8VttKU5fbJeVUemeX/KKzxVdXWBr8ouucxacRttEZm+L5WSy2torLh1lguZzSREDSJyK3AH3lWelwGLctEZEaGqqipwO9XV1YHe3+uee+h32200jh/PR7/5DVWOFLdh/GzCbivozzrMdrIJVMD7ZLt+VXZeMtvc3Myrr74auJ1evXoFen/lu+8y9Ac/QNvaeOfuu9F99yWcv8pgwspGmG251k5LS2aTBUTgKG1H8raNdUng/ztUGXTzzQyZPZuPzzqLj376U6oDFFUu5nXp0qWhtFMILmQ2kx5cBmzHm/5oNrANiOyEr73uuYd+N97I5ro6PvjNb+zIbQRlMGl8OrErPI8XkSX+4+t4he04EXkLGOd/Xwglk9nKd99l6LnnQlsb79x1Fy377lvoLpmQpcurC0eKAiqZvIbGL277338/H591Fu//9Kd25NYhruQ1k1kUtuBNdxR58cXtmltvhZ07C90lkyNB9i6L4KrskshsfHG7asYMWoYOLXSXTI64cDQoV0olr6GJK27XffvbvH/ZZVbcOsiFzKYtcEVkf+AKvMnr21+vqsfnrlv5t1txW1FhBW6EReCoT0qlkNnE4nb7fvvBtm2F7pbJEcurAXYrblf/7GeQ4TAXk18uZDaT8++zgTuBu4HW3HanMJIWtyayQpiE2nWRzmzS4tZEluXVAMmLWweKKLM7VzKbSSW3U1Wn5bwnBWLFbWlyIXw5FNnMWnFbmiyvJc6K26LjQmYzqeaeFJEfAI8Rd39sVf0kkxWISDnQAHyoqgW/t233J56g7y23ULF6NW21tZR/+qkVtyXIhdMnOdTpzLqW19onn6TfrbdSsXo1O/v3p6y5Ge3SxYrbEmN5Tc61vIap19NPs8ftt1P58ce0de1K+ZYtVtwWERcym0lFd57/78/ilimQ6USTPwZWAD2y6FdOdH/iCQZcfTVl/li98k8/RcvK2DJ6tBW3JcSV0yc5FCSzzuS19sknGfjLX7bntXLtWhRYP2GCFbclxPLaIWfyGqZeTz/NkGuv/WxbvWULWl5O8/DhVtwWAVcym7YHqrpXkkdGxa2IDAFOxBtbVHB9b7mlPTAx0tZG39tuK0yHTMG4MIVJrnQ2s67ltd+tt+6eV6DXgw8WpkOmYFyZdigXopLXMO1x++27Z7+1lT1uv71APTLZciGvaQtcEekqIr8Uken+9/uJSKanQm4Dfg60ddD+xbFbom7cuDHDZjunYvXqrJab6HLhNoK5EiCzt2F5NQ6K8q1685XXdevWpWxo7ty5jB49mokTJ1JfX8+oUaOYMGFClp8kvHYqP/44q+WlxrXfVzIu5DWTNd2HNwn1Mf73q4Ap6d7kB3StqnZ4RxZVna6qI1V1ZO/evTPoTue11dYmXb5z0KCcrte4x4W9yxzKOrMu5nVn//7Jl1teS06Uj+CSp7z2T5EngAULFlBfX09tbS1Tp05l/vz5jBgxguXLl2f1QUJpR5W2rl2TPrVjjz2y6k9UOfX7SsGFvGZS4O6jqjcBOwBUtZnUk9zHOxY4WURWArPw7vp0f2c7GlSve+5pH3Mbr626mg0h7bGY4hAbH1Tovcsc6kxmncpr5bvveheUJSxvq65m/eWXF6RPpjDS5TUCmS14XlW9pPXp04empiZaW1tR1fbleWvHny0hNuY2Xlt1NR9fdllW/YkqZ35fKbiyjc1kTdtFpAZv0Dsisg9xV3qmoqq/UNUhqjoMOBv4q6r+R5DOdtYuU4HdeCM7Bg9GRdgxeDBrr7+eppNPLkS3TAG5EL4cyjqzLuU1NhWYdunC+gkTdsnrmilT2PyNbxSiW6aAIl7gFjyvY8eOZcyYMaxYsYLJkyczfvx4GhoaGD58eP7aSZgK7IPrrmP7oEGoCNsHDWLVpElsOvHETn7CaHHi95WGC3nNZOqAa4BngKEi8gDenuP5uexUmJLNc9t02mmF7pYpsAic1uxI0WY22Ty3Gy++uNDdMgUWNK8ici8QO61/sL9sMnAREBucepWq/inQijqn4Hmtq6ujrq6u/fuFCxfmt50U89xaQZtcwX9fGQiS2bDymrbAVdV5IrIYOArvtMmPVXV9Np1V1WeBZ7N5TxjsJg4mGVemMMmVoJktVF7tJg4mmZDy+nvgd8AfEpZPVdWbgzYeRLHmNTR2E4fICSGzvyeEvKat+ETkK/6Xm/1/DxIRVDW8Uj8HrLg1HYnyEdxizKwVt6YjQfOqqgtFZFg4vQlXMeY1NFbcRlaQzIaV10yqvvjJp6uBI4FFwPFBV54rVtyadKJ8BJciy6wVtyadHOb1UhH5Dt7dwCaoam7nvkuuqPIaGituIy1Hmc0qr5nc6OEbcY9xwMHAmnD6Gr6+M2ZYcWs6FHTKIRG5V0TWisiyuGWTReRDEVniP76e0w/RgWLKbJeVK624NR1Kl1c/s/1i8736j0wGbk8D9gEOBVYDt+TuU6RWTHkNjRW3keZKXjtT/a3CC6Bz+s6YwaCbb7bi1qTlwvigPHIys11WrmSfCy+04taklUFe16vqyGzaVNX2IlJE7gKe6kTXcsHJvIZGFS6/3IrbiEuT2bzkNZMxuLdD+3SUZXjV89JsOpap8vJyevbs2an31k6fTq+bb2briSey5uabI1nctrSknZ0tY42NjaG0s3bt2lDa+eCDD0JpJ1NBClyXx/NB/jIbJK8V77zDwIsuQlT5IKLF7baEW40GsWnTplDaWbMmnAODYfUn03F6uTjdKSKDVDV2W7zTgGUdvT5X8rmNDcObb77Z+Ter0v+GG+g9Y4ZzxW2Yef3DHxKPfZSesDPbmbxmUgU2xH29E3hQVZ/vRP9ypnb6dHpNmcLWE09kw+23Q2trobtkHJbDu6m4MJ4PHM9sxTvvMPDss5G2Ntb88Y9s//znC90l47Aw8ioiDwJfxTs1ugpvaq6visiheMXlSuA/A62k85zOa2jiituN3/kOqy+/3Jni1oQraGbDymsm04TN6HQv82C34raiwgpck1aavct+IhK/0ZmuqtPTNDkNuA4vfNfhjQ/6bqBOdpLLmU0sbnfsvz+EeGbCRFPQo0Gqek6SxfcEajQkLuc1NAnF7bqrrrLcR1zAs6Sh5DWTIQqvwW53zARvvj5V1UOyXWlYkha3xmQgzd5lUY/nczWzSYtbYzIQ8Wn9nMxraJIVtxH+fRqPC5nNpCKc6/870//328BWoKB7nVbcms7KxY0eXBnP53Mus1bcms6K+o1ZcDCvobHitiS5ktlMqsJjVfXYuO+vFJHnVfXaXHUqma6PPUbPm26i/KOPaOvRg/LGRituTae5MD4ohwqe2a7/93/0/u1vKf/oI1r790eam6Gqyopb0ykuHA3KoYLnNUy1Tz5Jv1tvpWL1atq6dqV8yxYrbkuQC5nNpDLsJiLHqepzACJyDNAtt93aVdfHHqP3lVdS1twMQHljI1peTvPYsVbcmk5xYXxQDhU0s13/7//o+4tftOe1Yu1aFNj0gx9YcWs6xYWjQTlU8G1sWGqffJKBv/wlZf6MBOVbtqDl5Wz74hetuC0xLmQ2k+rwe8C9ItIT7+hUI3m+eKbnTTe1byxjpLWVnjffzNbTT89nV0wEuHL6JIcKmtnev/3t7nkFau+/n08vuSRf3TARYXktHv1uvbW9uI2R1lb6TZ3K5pNPLlCvTL65ktlM7mS2SFVHAIcAh6rqoaq6OPdd+0z5Rx9ltdyYdILcycx1hc6s5bV0zZ07l9GjRzNx4kTq6+sZNWoUEyZMCNxu0LsPuqzQeYXwfm8Vq1dntbzUhPVzdq2dZFzIa9oCV0QGisg9wEOq2igiB4nI9/LQt3atgwdntdyYdMrKylI+il2hM2t5LV0LFiygvr6e2tpapk6dyvz58xkxYgTLly8P1G5HeS32zBY6rxDS702Vtq5dkz61c9CgkHpa3MLKh2vtJONCXjNZ0++BPwOxrdObwE9y1J+kGn/+c9qqqnZZ1lZTQ+PPf57PbpgIcWHvMod+TwEzu/FnP6OtpmaXZW01NWz82c/y1QVTIKrebFd9+vShqamJ1tZWVLV9eWdF+QguDmxjA//e/NkSYmNu47VVV7P+8svD7nJRCisfrrWTjAt5zaTA7aeqDwNtAKq6E8jrnRS2nnZaezGrwM4992TjjTey9bTT8tkNExGx8UGF3rvMoYJmduupp7Lhhhtoq6pqz+uGG25g66mn5qsLpkDGjh3LmDFjWLFiBZMnT2b8+PE0NDQwfPjwTreZLq8RyGzBt7GBfm8JU4F9fOON7Bg8GBVhx+DBrJkyhc3f+EbuP0QRCCsfrrWTyJVtbCYXmW0Rkb74E1GLyFF4g+Dzqnn8eHpfdx2f3HorW884I9+rNxETgaM+HSl4ZreeeirdH3kEaW5mzZw5+Vy1KaC6ujrq6urav1+4cGEo7Vpec6vTv7cU89zaBWXJhZUP19pJxoXMZlLgXg48AewjIs8D/QGrME3RcuUKzxyyzJrIsLw6ym7iYFJwJbNpC1xVXSwio4AD8Gb7+Zeq7sh5z4zJIRfClyuWWRM1llfHWHFr0nAhsyl7ICJfEpE9oH1M0BHA9cAtItInT/0zJidcGAAfNsusiaooXmRWtHm14tZkwIW8dlRi/y+wHUBEvgLcCPwBb2zQ9Nx3zZjccGUAfA5YZk3kRPgis+LLqxW3JgOubGM7GqJQrqqf+F+fBUxX1TnAHBFZkvOeGZNDxXrUJw3LrIkky6sDVOHyy624NRlxIbMdFrgiUuGfOhkDXJzh+zqtrKyM6urq5M/58+BWVlamfI2rWlpaCt0Fk6CIj/p0JK+ZFZEOs1hWVoZ0kGlXbUu41WiUrFy5stBd6BTLa3A7duxgzZo1nXuzKt2vuYZu06fz6QUX8OmkSVQ5UMBAeHm94oorQmnHeFzIbEc9eBD4m4g8DjQDfwcQkX0pwDRhxoTJhfFBOWCZNZEUdAyuiNwrImtFZFncsj4iMk9E3vL/7Z3TD7G74shrXHG75aKL2Dhpkh25NWm5kNeUBa6qXg9MwLvLynH62a0tyoDL0vbQGEe5Mj4obJZZE0UhjcH9PXBCwrIrgQWquh+wwP8+b4oirwnFbdO111pxa9IKYRv7e0LIa4enQVT1pSTL3sykd8a4LEghKyL3AicBa1X1YH9ZH+AhYBiwEvimqm4M3NEsWWZNFAXd8VTVhSIyLGHxKcBX/a9nAM8C/xVoRVlyOq9W3JoAgmQ2rLwW7+EqYwIIOETh9zh2NMiYKMtgiEI/EWmIe1ycrk1goKquBvD/HZDLz1BUrLg1AbmQ15xcLAYgIkPxpjzZA+8e29NV9b9ztT5jMhU7fdJZrh4NCsLyalyVYV7Xq+rIfPTHFTnLrBW3JqAMMpuXvObyCO5OYIKqfgE4CvihiBzU2ca6PPMMALWXXUafww+nyu5vbwLIwUVmxX40KNS8Vs+ZQ5eXX6byH/+g38iRVFteTQBBLzJLYY2IDPLbHwSsDa3D+RFaZqvnzKHfyJEMGDSI/vvtZ8WtCcyFvOaswFXV1aq62P96M7AC2LMzbVXNmUP3X/8a8O5jWL5qFbWXX25Frum0NAPgO3P6pKiFmdfqOXOoveIKpKXls7xecYUVuabTcnSjhyeA8/yvzwMeD6WzeRJWZmN5LV+1ClGlrKkJrahg56GHWnFrOs2FvOZlDK5/Ovcw4OXOvL/b9dcjCXPdSXMz3a6/PnjnTMnJ4GjQelUdGffI5K5CxX40qF3QvHa/4QbKmpt3WVbW3Ez3G24I3jmTE3PnzmX06NFMnDiR+vp6Ro0axYQJEwrWTrx0ec1w2qEHgReBA0RklYh8D+/OYeNE5C1gnP99UQqS2WR5lZ07La9xXMuHa+0kciWvOS9wRaQ7MAf4iap+muT5i2NHyjZs2JC0jbIPP8xquTHpuLB36aJs8vrJJ5/s3gCW12K0YMEC6uvrqa2tZerUqcyfP58RI0awfPnygrSTKOgRXFU9R1UHqWqlqg5R1XtUdYOqjlHV/fx/k/9BO66jzFpew+FaPlxrJxkX8prTAldEKvGC94CqPprsNao6PXakrG/fvknbadsz+VmXVMuNSSdI+KJ6NCjbvPbp0ydpO5bX4hObgrVPnz40NTXR2tqKqrYvz3c7iXI0RKHopcus5TUcruXDtXaScSGvOVuTeMeh7wFWqOqtQdracvXVaMLtPrWmhi1XXx2kWVOigp4+ieLRoDDz2vSLX9BWU7PLsraaGpp+8YsgzZocGjt2LGPGjGHFihVMnjyZ8ePH09DQwPDhwwvSTrwwhihEUViZtbym51o+XGsnkSt5zdk0YcCxwLnAayKyxF92lar+KduGWk4/HVm3jtpJk1CgbcgQtlx9NS2nnx5id00pKeWjPimEltdtfi57XH45tLTQNmQITb/4Rfty4566ujrq6urav1+4cGFB20lkeU0qlMzGctl94kTKP/mE1oEDaZo0yfIax7V8uNZOMi5kNmcFrqo+hzfpQSi2n3ACTJrE5ttvp+Wss8Jq1pSoUj3qk0rYed12+ulUP/QQ0tzMxiefDKtZU6Isr7sLM7PbTj8drayk18UXs/Ghh2g98MAwmjUlzIXM5vIIrjFOymASamOMIyyvxhQXVzJrBa4pSS7sXRpjMmN5Naa4uJBZK3BNSXJh79IYkxnLqzHFxYXMWoFrSo4rp0+MMelZXo0pLq5k1qkCV0SoqqpK/lyXLgBUVlZCiteEraWlJS/rMfnnwumTYldWVkZ1wvR9ic9LmteEKcp5feKJJwrdhYKyvAaXLq8V/ja2qqqKtjxkNqy8HnnkkaG0Y8LlQmadKnCNyRcX9i6NMZmxvBpTXFzIrBW4piS5sHdpjMmM5dWY4uJCZq3ANSXHlfFBxpj0LK/GFBdXMmsFrilJLuxdGmMyY3k1pri4kFkrcE1JcmHv0hiTGcurMcXFhcwWvgcZqviTd3vt6u9/n+4HH0zFww8XuEemWMVOn6R6mOAqZ8+m4sUXKX/5ZWq/+EUqZ88udJdMkUqXV8tscJWzZ1MzYQIA3U47zfJqAnFlG1sUR3ArHn6Y6muvBbwbb8sHH1Dzox/RDOz85jcL2jdTnFw4fRJVlbNnU/OjHyH+NECxvALsOPPMQnbNFKmgeRWRlcBmoBXYqaojQ+hWJLTntbkZgLKPP7a8msBcyGxR7PpWX3stsm3bLsukubm96DUmWy7sXUZV9bXXtm8sYyyvuTF37lxGjx7NxIkTqa+vZ9SoUUzwj8QVczuJQjqCO1pVD7XidldRzqtrf9dRbSeZkLaxgTJbFFtzWbUqq+XGdEREOnxk2MZKEXlNRJaISEOOu1xULK/5s2DBAurr66mtrWXq1KnMnz+fESNGsHz58qJuJ166vNrZmGCinFfX/q6j2k4iV/JaFAWuDhmS1XJj0nFh7zKqLK/5o6oA9OnTh6amJlpbW1HV9uXF2k6iEI7gKvAXEVkkIhcH6kzERDmvrv1dR7WdZELYxgbObFEUuNsmTUITbh2oNTVsmzSpQD0yxc6Fvcuo2jZpElpTs8syy2tujB07ljFjxrBixQomT57M+PHjaWhoYPjw4UXdTqIMjgj1E5GGuEfiBvFYVT0cqAN+KCJfCdShCIlyXl37u45qO8kEzCuEkFkJo1IPy2GHHaZ/+9vfkj5Xeccd1Fx1FQro0KFsmzQp5QVmYd3jurGx0al2XnjhhVDaAXj++edDaWfZsmWhtPP666+H0o6qpq1Qv/jFL+qjjz6a8vn9999/UbqjsiLyLrARby/zf1V1erZ9LXaHHXaYPvvss0mfq5w9m5pLL4WWlva8prpgZVvC+PrOCitnmzZtCqWdhQsXhtIOwBVXXBFaWy4RkbRZS5dXyCyzceucDDSp6s0ZdzQC0uW1+sorKduwgbY99mDbddcVTV7333//UNox6WWSVwhnG5uw3sl0IrNFMYsCwM6vfx2uuoptd97JjnPOKXR3TJFLc5qkX8K42ulJCthjVfUjERkAzBORN1Q1vIqmyO0480wqH3gAaW5my5//XOjumCIX5OJPEekGlKnqZv/r8UDxX0EVoh1nnolWVtLt/PPZ8thjtH3hC4XukilyLmS2aApcY8KUZijC+nR7l6r6kf/vWhF5DDgSsALXmBwIOHRoIPCY30YF8EdVfSaMfhljknMhs1bgmpIT9D7ZdkTImPwJmldVfQcYEV6PjDEdcSWzVuCakuTC3qUxJjN28acxxcWFzFqBa0qSC3uXxpjM2A1YjCkuLmTWqQK3rKyMqqqq5E926QJARUUFZale4wtrFgUTTUFPnxhPWVkZ1QnT98UrLysDkQ5fA+Fdle2aqM58kG+W13Cky2tZZSUAVVVVu03LGS+svNrsB9HlSmadKnCNyRcXTp8YYzJjeTWmuLiQWStwTUlyYe/SGJMZy6sxxcWFzFqBa0qSC3uXxpjMWF6NKS4uZNYKXFNyXBkfZIxJz/JqTHFxJbM57YGInCAi/xKRt0Xkylyuy5hspLlPdkmyvBpXdZRXy6xl1rjHhbzmrMAVkXLgDqAOOAg4R0QO6mx7ZU89BUDFhRfSZf/9KZs1K5R+mtJUVlaW8lGKQs/rrFnI888jL71keTWBdZRXy2zwzJbNmkXFj38MQOWJJ1peTWAu5DWXQxSOBN725wxFRGYBpwDLs22obNYsKiZNAkAA3n+fih/8gJ1A29lnh9ZhUxpcOX3imHDz+oMfILHp+iyvJgDLa0qhZLY9r1u3AiCrV1teTSCuZDaXPdgT+CDu+1X+sqxVTJqEJMy9J1u3the9xmTLhdMnjgk3r/7GMsbyuqu5c+cyevRoJk6cSH19PaNGjWLChAnWTgo2RCGpUDLrYl5d+3u0drLnQl5zWeAm+xS624tELhaRBhFpWLduXfKWPvggu+XGdCC2d1no0yeOsbzm0YIFC6ivr6e2tpapU6cyf/58RowYwfLl2R0wj2o78dLl1TK7i10yW6x5de3v0drJjivb2FyuaRUwNO77IcBHiS9S1emqOlJVR/bv3z95S0OHZrfcmDRc2Lt0jOU1j1S9OqRPnz40NTXR2tqKqrYvL/V2EtkR3KTSZrZY8+ra36O1kz0X8prLAvcfwH4ispeIdAHOBp7oTEM7r70W7dp1l2XatSs7r702eC9NSXJh79Ixltc8Gjt2LGPGjGHFihVMnjyZ8ePH09DQwPDhw62dJOwIblKhZNbFvLr292jtZM+FvEoYlXrKxkW+DtwGlAP3qur1Hb3+iCOO0BdeeCHpc+0Xmn3wAQwdys5rr005AP7TTz8N1nFfY2OjU+2k+tl0xvPPPx9KO8uWLQulnddffz2UdlQ17e7h4Ycfrh19/q5duy5S1ZGhdKiIFCqvmzZtCtTvmLByFlZ/jjzyyFDaiTIRSZu1dHkFyywZZNa1vA4cODCUdkz+ZJJXcGcbm9MbPajqn4A/hdFW29lns92u6DQhcOUKT9dYXo2LLK+phZVZy6sJkyuZLXwPjCmAoKdPbIJ1Y/In6BAFy6sx+eXCNtYKXFOSggyAD/umCMaYjgW5yMzyakz+ubCNzekQBWNcFMLpk9BuimCM6Zjl1Zji4kpm7QiuKUkBpzAJ7aYIxpj0Ak4TZnk1Js9c2MY6dQR38eLF66urq99L87J+wPp89CdD1p+O5bM/n8/kRYsWLfpzWVlZvw5eUi0iDXHfT1fV6XHfZ3RThKizvIbCtf6AY5nNIK/QcWYtrxRtXsG9PpVyf4pqG+tUgauqKWai/oyINLg0HYz1p2Ou9QdAVU8I2ERGN0WIOstrcK71B9zrk+U1HMWYV3CvT9af9FzJrA1RMCZ7od0UwRiTc5ZXY4pLKJl16giuMcVAVXeKyKXAn/lsgvVw7lRhjAmV5dWY4hJWZouxwJ2e/iV5Zf3pmGv9CUWYN0WIONd+/9af9FzsUyCW14y5+Lt3rU/WnzwII7M5vVWvMcYYY4wx+WZjcI0xxhhjTKQUTYHr0q0WRWSoiNSLyAoReV1EflzI/sSISLmI/FNEnip0XwBEpJeIPCIib/g/q6ML3SeTHy7l1e+PZTZ9XyyvJcylzFpeM+qL5TWNohii4N+27U1gHN70Ef8AzlHVgtyJRkQGAYNUdbGI1AKLgFML1Z+4fl0OjAR6qOpJheyL358ZwN9V9W7/SsiuqrqpwN0yOeZaXv0+WWbT98XyWqJcy6zlNaO+WF7TKJYjuO23bVPV7UDstm0FoaqrVXWx//VmYAUFvjOOiAwBTgTuLmQ/YkSkB/AV4B4AVd1u4SsZTuUVLLMZ9MXyWtqcyqzlNW1fLK8ZKJYC19lbLYrIMOAw4OUCd+U24OdAW4H7EbM3sA64zz+lc7eIdCt0p0xeOJtXsMymYHktbc5m1vKalOU1A8VS4Dp5q0UR6Q7MAX6iqp8WsB8nAWtVdVGh+pBEBXA4ME1VDwO2AAUfi2nywsm8gmW2A5bX0uZkZi2vKVleM1AsBa5zt1oUkUq84D2gqo8Wsi/AscDJIrIS79TS8SJyf2G7xCpglarG9rofwQukiT7n8gqW2TQsr6XNucxaXjtkec1AsRS4Tt1qUUQEb+zLClW9tVD9iFHVX6jqEFUdhvez+auq/keB+/Qx8IGIHOAvGgMU9AIBkzdO5RUssxn0x/Ja2pzKrOU1bX8srxkoijuZOXirxWOBc4HXRGSJv+wq/84b5jOXAQ/4/2G+A1xQ4P6YPHAwr2CZzYTltUQ5mFnLa3qW1zSKYpowY4wxxhhjMlUsQxSMMcYYY4zJiBW4xhhjjDEmUqzANcYYY4wxkWIFrjHGGGOMiRQrcI0xxhhjTKRYgRuAiFwtIq+LyKsiskREvlzoPnWWiFwsIg/Ffd9DRP6fiOxVyH4ZExbLqzHFxTJrgrACt5NE5GjgJOBwVT0EGMuu9/IuNncBQ0RkrP/9tXhzIb5bwD4ZEwrLqzHFxTJrgrICt/MGAetVtQVAVder6kcAInKEiPxNRBaJyJ9FZFDc8qUi8qKI/FZElvnLzxeR38UaFpGnROSr/tfj/dcvFpHZ/r25EZGVIvIrf/lrInKgv7y7iNznL3tVRE7vqJ0Y9SZEvgS4TURG4t0Z5bc5/PkZk0+WV2OKi2XWBGIFbuf9BRgqIm+KyP+IyChov3/27cAZqnoEcC9wvf+e+4AfqerRmaxARPoBvwTGqurhQANwedxL1vvLpwFX+MsmAo2q+kV/r/evGbQDgKq+incnmwV+P7dn+sMwxnGWV2OKi2XWBFIUt+p1kao2icgRwL8Bo4GHRORKvD/sg4F5IgLebQ9Xi0hPoJeq/s1vYiZQl2Y1RwEHAc/7bXUBXox7/lH/30XAv/tfj8W7V3asnxtF5KQ07cS7A6hT1fo0fTOmaFhejSkullkTlBW4AahqK/As8KyIvAachxeE1xP3IEWkF5Dqvsg72fVoenXsbcA8VT0nxfta/H9b+ex3KUnWk66deG3+w5hIsbwaU1wssyYIG6LQSSJygIjsF7foUOA94F9Af/EGyCMilSIyXFU3AY0icpz/+m/HvXclcKiIlInIUOBIf/lLwLEisq/fVlcR2T9N1/4CXBrXz96dbMeYyLC8GlNcLLMmKCtwO687MENElovIq3inJyb7Y2rOAH4jIkuBJcAx/nsuAO4QkReB5ri2ngfeBV4DbgYWA6jqOuB84EF/HS8BB6bp1xSgt4gs89c/upPtGBMllldjiotl1gQi3oV9Jt9EZBjwlKoeXOi+GGM6Znk1prhYZo0dwTXGGGOMMZFiR3CNMcYYY0yk2BFcY4wxxhgTKVbgGmOMMcaYSLEC1xhjjDHGRIoVuMYYY4wxJlKswDXGGGOMMZFiBa4xxhhjjImU/x9JjZKqpxcRpAAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 720x216 with 6 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "X = [1, 1, 1, 1, 2, 1, 1, 6, 6]\n",
    "Y = [0, 3, 3, 3, 9, 9, 7, 7]\n",
    "\n",
    "print(r'Accumulated cost matrix and optimal warping path for different constraint regions:')\n",
    "plt.figure(figsize=(10, 3))\n",
    "global_constraints = True\n",
    "\n",
    "ax = plt.subplot(1, 3, 1)\n",
    "band_rad = 1\n",
    "title='Sakao-Chiba band (rad = %.2f)'%band_rad\n",
    "compute_plot_D_P(X, Y, ax=ax, global_constraints=global_constraints, band_rad=band_rad,title=title)\n",
    "\n",
    "ax = plt.subplot(1, 3, 2)\n",
    "band_rad = 0.5\n",
    "title='Sakao-Chiba band (rad = %.2f)'%band_rad\n",
    "compute_plot_D_P(X, Y, ax=ax, global_constraints=global_constraints, band_rad=band_rad,title=title)\n",
    "\n",
    "ax = plt.subplot(1, 3, 3)\n",
    "band_rad = 0.25\n",
    "title='Sakao-Chiba band (rad = %.2f)'%band_rad\n",
    "compute_plot_D_P(X, Y, ax=ax, global_constraints=global_constraints, band_rad=band_rad,title=title)\n",
    "\n",
    "plt.tight_layout()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Multiscale DTW"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "When using the concept of global constraint regions, one needs to make sure that the optimal warping path to be computed actually lies  within this region. Since this path is not known \n",
    "a priori, it is often difficult to find a good trade-off between choosing the constraint region as small as possible (to speed up computations) but large enough to contain the desired path. One possible strategy to increase the probability of finding the \"right\" path is to use data-dependent constraint regions instead of a data-independent, fixed constraint region. This idea can be realized by a [**multiscale approach**](https://www.audiolabs-erlangen.de/fau/professor/mueller/publications/2006_MuellerMattesKurth_MultiscaleAudioSynchronization_ISMIR.pdf) to DTW, where the general strategy is to recursively project an optimal warping path computed at a coarse resolution level to the next higher level and then to refine the projected path. The following figure summarizes the main steps of such an approach: An optimal warping computed on a coarse resolution level is projected onto a finer resolution level. Along with a small neighborhood (to increase robustness of the overall procedure), this defines a constraint region used to compute a warping path on the finer resolution level. For details, see Section 3.2.2.4 of <a href=\"http://www.music-processing.de\">[Müller, FMP, Springer 2015]</a>.\n",
    "\n",
    "<img src=\"../data/C3/FMP_C3_F19.png\" width=\"500px\" align=\"center\" alt=\"C1\">"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Further Notes\n",
    "\n",
    "* The article [\"An Efficient Multiscale Approach to Audio Synchronization\"](https://www.audiolabs-erlangen.de/fau/professor/mueller/publications/2006_MuellerMattesKurth_MultiscaleAudioSynchronization_ISMIR.pdf) shows how the multiscale approach can be applied for [music synchronization](../C3/C3_MusicSynchronization.html).\n",
    "* Extending this work, the article [\"Memory-Restricted Multiscale Dynamic Time Warping\"](https://www.audiolabs-erlangen.de/fau/professor/mueller/publications/2016_PraetzlichDriedgerMueller_MrMsDTW_ICASSP.pdf) introduces a DTW approach that works with a pre-defined memory quota."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div class=\"alert\" style=\"background-color:#F5F5F5; border-color:#C8C8C8\">\n",
    "<strong>Acknowledgment:</strong> This notebook was created by <a href=\"https://www.audiolabs-erlangen.de/fau/assistant/zalkow\">Frank Zalkow</a> and <a href=\"https://www.audiolabs-erlangen.de/fau/professor/mueller\">Meinard Müller</a>.\n",
    "</div> "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<table style=\"border:none\">\n",
    "<tr style=\"border:none\">\n",
    "    <td style=\"min-width:50px; border:none\" bgcolor=\"white\"><a href=\"../C0/C0.html\"><img src=\"../data/C0_nav.png\" style=\"height:50px\" alt=\"C0\"></a></td>\n",
    "    <td style=\"min-width:50px; border:none\" bgcolor=\"white\"><a href=\"../C1/C1.html\"><img src=\"../data/C1_nav.png\" style=\"height:50px\" alt=\"C1\"></a></td>\n",
    "    <td style=\"min-width:50px; border:none\" bgcolor=\"white\"><a href=\"../C2/C2.html\"><img src=\"../data/C2_nav.png\" style=\"height:50px\" alt=\"C2\"></a></td>\n",
    "    <td style=\"min-width:50px; border:none\" bgcolor=\"white\"><a href=\"../C3/C3.html\"><img src=\"../data/C3_nav.png\" style=\"height:50px\" alt=\"C3\"></a></td>\n",
    "    <td style=\"min-width:50px; border:none\" bgcolor=\"white\"><a href=\"../C4/C4.html\"><img src=\"../data/C4_nav.png\" style=\"height:50px\" alt=\"C4\"></a></td>\n",
    "    <td style=\"min-width:50px; border:none\" bgcolor=\"white\"><a href=\"../C5/C5.html\"><img src=\"../data/C5_nav.png\" style=\"height:50px\" alt=\"C5\"></a></td>\n",
    "    <td style=\"min-width:50px; border:none\" bgcolor=\"white\"><a href=\"../C6/C6.html\"><img src=\"../data/C6_nav.png\" style=\"height:50px\" alt=\"C6\"></a></td>\n",
    "    <td style=\"min-width:50px; border:none\" bgcolor=\"white\"><a href=\"../C7/C7.html\"><img src=\"../data/C7_nav.png\" style=\"height:50px\" alt=\"C7\"></a></td>\n",
    "    <td style=\"min-width:50px; border:none\" bgcolor=\"white\"><a href=\"../C8/C8.html\"><img src=\"../data/C8_nav.png\" style=\"height:50px\" alt=\"C8\"></a></td>\n",
    "</tr>\n",
    "</table>"
   ]
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "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.8.16"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}