{ "metadata": { "name": "Time Series Classification and Clustering" }, "nbformat": 3, "nbformat_minor": 0, "worksheets": [ { "cells": [ { "cell_type": "heading", "level": 1, "metadata": {}, "source": [ "Time Series Classification and Clustering" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "In a typical classification problem you are given a set of input features and a set of discrete output classes and you want to model the relationship between the two. There is a myriad of classification algorithms that you could use for this problem - SVMs, Naive Bayes, k-NN, etc. But what if the input features are not independent such as with time series data? In this case SVMs and Naive Bayes would not be a good choice since they assume that the input features are independent. The k-NN algorithm could still work however it relies on the notion of a similarity measure between input examples. Now the question becomes _how do we measure the similarity between two time series_?" ] }, { "cell_type": "heading", "level": 2, "metadata": {}, "source": [ "How about Euclidean distance?" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The Euclidean distance between two time series $Q$ and $C$ of length $n$ is defined as\n", "\n", "$$d(Q,C) = \\sqrt{\\sum^n_{i=1}[Q(i)-C(i)]^2}$$\n", "\n", "At first glance, it seems like simply calculating the Euclidean distance between two time series would give us a good idea of the similarity between them. After all, the Euclidean distance between identical time series is zero and the Euclidean distance between very different time series is large. However, before we settle on Euclidean distance as a similarity measure we should clearly state our desired criteria for determining the similarity between two time series " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "With a good similarity measure, small changes in two time series should result in small changes in their similarity. With respect to Euclidean distance this is true for changes in the y-axis, but it is not true for changes in the time axis (i.e. compression and stretching). Consider the following example." ] }, { "cell_type": "code", "collapsed": false, "input": [ "import pandas as pd\n", "import numpy as np\n", "import matplotlib.pylab as plt\n", "\n", "x=np.linspace(0,50,100)\n", "ts1=pd.Series(3.1*np.sin(x/1.5)+3.5)\n", "ts2=pd.Series(2.2*np.sin(x/3.5+2.4)+3.2)\n", "ts3=pd.Series(0.04*x+3.0)\n", "\n", "ts1.plot()\n", "ts2.plot()\n", "ts3.plot()\n", "\n", "plt.ylim(-2,10)\n", "plt.legend(['ts1','ts2','ts3'])\n", "plt.show()" ], "language": "python", "metadata": {}, "outputs": [ { "output_type": "display_data", "png": "iVBORw0KGgoAAAANSUhEUgAAAXYAAAD9CAYAAACoXlzKAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzsnXdYVEcXh39LEVEQrFjAXiiCgL0jimANliR2UWPsn91Y\nEkvURI29xhrsiRpj7J21oqCCrFhjBRV7R5Ey3x/jspQtd2/bXZ33eXh09947c3b27rkzZ05REEII\nGAwGg/HZYGVqARgMBoMhLkyxMxgMxmcGU+wMBoPxmcEUO4PBYHxmMMXOYDAYnxlMsTMYDMZnBifF\n3rt3b7i4uMDb2zvzvefPnyMoKAiVK1dG8+bN8fLlS8mEZDAYDAZ3OCn2Xr16Yf/+/dnemzFjBoKC\ngnD9+nU0bdoUM2bMkERABoPBYBiHgmuA0p07d9CmTRuoVCoAgLu7O44dOwYXFxckJSUhICAAV69e\nlVRYBoPBYBiGt4390aNHcHFxAQC4uLjg0aNHognFYDAYDP7YiNGIQqGAQqHQeYzBYDAYxsM34wvv\nGbvaBAMADx8+RLFixXSeSwhhf4Rg0qRJJpfBXP7YWLCxYGOh/08IvBV727ZtsXbtWgDA2rVrERoa\nKkiQL4E7d+6YWgSzgY2FBjYWGthYiAMnxd65c2fUq1cP165dg5ubG/744w+MHTsWhw4dQuXKlXH0\n6FGMHTtWalkZDAaDwQHOXjG8O1AoBC8rPheUSiUCAgJMLYZZwMZCAxsLDWwsNAjRnUyxMxgMhhki\nRHeylAIyolQqTS2C2cDGQgMbCw1OTk6ZXnZfyl+hQoVEH0dR3B0ZDAZDDF6/fv3FrfClcAlnphgG\ng2E2fIn6QtdnZqYYBoPBYGTCFLuMMFuqBjYWGthYMMSGKXYGg8H4zGCKXUaYf64GNhYa2FhYBmXL\nlsXRo0cNnrdlyxbUq1cP+fPnR5MmTWSQLDfMK4bBYDA4wHUzs3DhwhgxYgSuXLnC6UEgBWzGLiPM\nlqqBjYUGNhbmT/fu3XHv3j20adMGjo6OmDp1Krp164YiRYqgYMGCqFWrFh4/fgwAaNq0KTp27IgS\nJUqYTF6m2BkMBsMA69evR+nSpbF79268efMGxYoVw+vXr5GYmIjnz59j+fLlsLe3N7WYmTDFLiPM\nlqqBjYUGNhbcUSjE+RNKnjx58OzZM9y4cQMKhQJ+fn5wdHQU3rBIMMXOYDAsBkLE+RNK9+7dERwc\njE6dOqFUqVL44YcfkJaWJrxhkWCKXUaYLVUDGwsNbCwsg6yh/zY2Npg4cSLi4+Nx+vRp7N69G+vW\nrdN5vtwwxc5gMBgccHFxwc2bNwEAERERUKlUSE9Ph6OjI2xtbWFtbQ0AyMjIwIcPH5CamoqMjAyk\npKQgNTVVVllZrhgGg2E2mLO+2LlzJ4YMGYLXr19j5MiR2LBhAxITE+Hg4IBOnTph7ty5sLKyQnh4\nOHr37p3t2rCwMKxZs0Zru1LkimGKncFgmA1for5gScAsHGZL1cDGQgMbC4bYMMXOYDAYnxnMFMNg\nMMyGL1FfMFMMg8FgMAzCFLuMMFuqBjYWGthYMMSGKXYGg8H4zGA2dgaDYTZ8ifqC2dgZDAaDYRCm\n2GWE2VI1sLHQwMaCITZMsTMYDAYHuJbGGzVqFCpXrowCBQrAw8MD69evl0G67LDSeDLC8m5rYGOh\ngY2FZcDV5u3g4IDdu3ejcuXKiIqKQkhICCpWrIi6devKICVF8Iz9119/hZeXF7y9vdGlSxekpKSI\nIReDwWCYDcaUxps8eTIqV64MAKhVqxYaNmyIyMhIWeUVpNjv3LmDlStX4sKFC5kpLP/880+xZPvs\nYLZUDWwsNLCxMH/4lsZ7//49oqOjUbVqVVnlFWSKKVCgAGxtbZGcnAxra2skJyejVKlSYsnGYDAY\n2VBMEad4BZkkzKUya2k8b29v+Pn5aT2vf//+8PX1RfPmzQX1ZyyCFHuhQoUwcuRIlC5dGvb29ggO\nDkazZs3Eku2zg9lSNbCx0MDGgjtCFbJYdO/eHQkJCejUqRNevnyJbt26Yfr06bCx0ajU0aNH4/Ll\ny4iIiJBdPkGK/ebNm5g/fz7u3LkDJycnfP3119i4cSO6du2a7bywsDCULVsWAODs7AxfX9/Mm1m9\nDGWv2Wv2mr02Z7SVxps4cSLu3r2Lli1bokqVKpkFNiZNmoQDBw7g2LFjcHBw4NS+UqlEeHg4AGTq\nS94QAfz555+kT58+ma/XrVtHBg4cmO0cgV18VkRERJhaBLOBjYUGNhYazFlf1KlTh6xYsYIQQsjR\no0dJXFwcSUtLI8+ePSPVqlUj4eHhhBBCpk+fTipVqkSSkpI4tavrMwsZC0Gbp+7u7jhz5gzev38P\nQggOHz4MT09PYU8aBoPBMEPGjRuHadOmoWDBgjh16hS+/vprODk5wdPTEwEBAejevTsA4Mcff0RC\nQgIqVqwIR0dHODo6YsaMGbLKKjhXzKxZs7B27VpYWVnB398fq1atgq2traaDLzD3A4PB4MeXqC9Y\nzVMGg/FZ8yXqC5YEzMJRbxQx2FhkhY0FQ2yYYmcwGIzPDGaKYTAYZsOXqC+YKYbBYDAYBmGKXUaY\nLVUDGwsNbCwYYsMUO4PBYHxmMBs7g8EwG75EfcFs7AwGg8EwCFPsMsJsqRrYWGhgY2EZcC2NN2bM\nGJQuXRoFChSAq6srRowYgbS0NBkk1MAUO4PBYHCAq2mkT58+uHz5Ml6/fo2oqCgcPHgQq1atkkFC\nDUyxy4glpCaVCzYWGthYmD/GlMarUqVKZqpeQgisrKxQokQJWeVlip3BYDAMYGxpvBkzZsDR0RFu\nbm5o3bo1vvrqK1nlZYpdRpgtVQMbCw1sLIxAoRDnTyBZS+MpFAr4+fnB0dEx8/jYsWPx5s0bnD9/\nHhs3bsT27dsF92kMTLEzGAzLgRBx/gTSvXt3BAcHo1OnTihVqhR++OEHrRukfn5+GDhwINavXy+4\nT2Ngil1GmC1VAxsLDWwsLANtpfHi4+Nx+vRp7N69G+vWrdN6XWpqKvLnzy+XmACYYmcwGAxOuLi4\n4ObNmwCAiIgIqFQqpKenw9HREba2trC2tgYhBMuXL8fLly9BCEFUVBSWLl2K9u3byyorU+wywmyp\nGthYaGBjYRlwLY23Y8cOVKhQAU5OTujTpw+mTZsmu2K3kbU3BoPBsFDatm2Ltm3bZr7+8ccftZ63\nb98+uUTSCcsVw2AwzIYvUV+wXDEMBoPBMIhZK/YPH4CzZ0XxTjILuNhSr18H7t2TXhZTo20sHjwA\nRo8G3ryRXx5Toh6Lu3eBT8GLXwQfPwJHjtB/GeJidoo9PR346y/g22+B4sWB5s2BrVtNLZU8vHwJ\nNGkC+PsDdesC8+Z9OT/027eBhg2BEyeAkBDg9WtTSyQPMTHAsmWAlxfg4wO0bk1/A18C06cDXboA\npUoBQ4YA584BBQoUgEKh+KL+HB0Lij62ZqfYFy+mX3hQEJ297tsHDB9OlZ6lY8hfeexYoE0b4OFD\nYPJkIDqaKrnPZcWSlaxjcfky0KgRMHIkcPo0VXBfgnJ/84be5+7uAQgPB54/B/LkAVasMLVk0nPp\nErB0KXDhAhAVBRQtCrRqBezY8QqEkM/2LzWVICyMoHZtgqdPCc6eJcib9zkOHBB5gInEGNNFaioh\nZcoQcuZM9vf79SNkwABx5TI3TpwgpGRJQl680LyXnk6IhwchR46YTi6piY8npHhxQtat07yXkUHI\nwIGE1KlDyOvXppNNaubOJeSbb7K/FxdHSJEihCQlmUYmOUhLI6R2bUKWL8/+/sqVhLRqZRqZ5KJD\nB0JatCDk7VvNe6dO0e/88OHs5wpRz2al2P/8k5CGDXO///w5ISVKEBIZKaJgJiAiIkLr+ykphHh6\nErJlS+5jK1cS0rKltHKZAvVY9OxJyC+/5D6ekUFIu3ZU+X2OfPxIiJsbIdHRue+LUaMI6dbNNHLJ\nwbx5hDRqRCcuWUlOJsTZOYJcuWIauaTm4kX6nX/8mPvYsWO5v3Mhit1sTDGEAL/9BowalftYwYLA\n3LlAv35Aaqr8sknNrFlA+fJAx465j3XrBpw/D1y5Ir9cUpOcDPz7LxAWlvuYQgEMGABs3Ci7WLLw\n119AxYpAjRq5j02aBBw7BkREyC+X1Ny+DUybBqxcCVjl0D729kDbtsCCBaaRTWrWrwe6dgVsbXMf\na9SIHhcN3o8EjnDtQqkkpHLl3E9xNRkZhAQFEbJihYjCmQEvXhDi7EzI3bu6z5kyhZDvvpNPJrnY\ntImQ4GDdx9PSqJnm6lX5ZJKDjAxCvL0J2bdP9znbt9NzMjLkk0sO+vUjZNIk3ceTkujv4ckT2USS\nhbQ0anWIj+d+jRD1bDYz9tmz6eZZzqe4GoUC6NsX+OcfeeWSmv37gfr1gdKldZ8zYACwbRvw6JF8\ncsnBhg10RaILa2ugU6fPb9au3igLDtZ9Tmgo8OIFdSD4XMjIAHbupLNWXbi4AO3aAcuXyyeXHBw5\nApQsCXh6ytOfWSj2K1fozvinVAs6ad4cOHkSePdOHrnERpvv9q5d1BNGH0WLUvfPJUukkcsU/POP\nEqdOUQWmj65dqWL/nDyDfvuN+uurkwVquy8UCur6uHu3vLJJyYULgJMTUKmS7nOUSiWGD6f3ekqK\nfLJJzfr1hvWbmAhW7C9fvkTHjh3h4eEBT09PnDlzxug2Fi+ms9IsBUi04uQE1KoFHDrEU1gzIy2N\nzthbtzZ87vDhwO+/fz7BHEeP0gfapwpiOqleHbCxoYFqnwPXrgFXr9KViCE+N8W+cye1oRvC25v6\n9W/ZIr1McvD2LZ3Ade4sX5+CFfvQoUPRsmVLXLlyBXFxcfDw8DC6jX37tG8caqNNG8u92XP6sZ86\nBZQtSwM0DFGlCuDmRn3bPwfOng3Qa4ZRo1BoZu2fAwcPAi1aZN9A0xXfEBhIN85fvJBHNqnhsjpV\nj0XPnnRj/XNg+3agQQOgWDH5+hSk2F+9eoUTJ06gd+/eAGjyeScnJ6PauH2bekd4eXE7v3VrYM8e\naq+zdLjc6FkJDKS2Okvn2jUgIQFo2pTb+V260Nnb5+ARdeQI989tbw80bgzxg1dMwL17QGIijajm\nQpMmgFL5efzO5TbDAALT9t6+fRtFixZFr169cPHiRVSvXh0LFixAvnz5sp0XFhaGsmXLAgCcnZ3h\n6+ub+WRetkwJLy9AoaCv1fZG9fGcrxMSlMiTBzh3LgC1ahk+35xeZ7WlBgQEYNcuYORIJZRKbtc3\nbQqMHatEo0bm8Xn4vl6zhs5gbGwCOF9frlwADh8G7O1NLz/f12lpwJEjSvTsCQCa47GxsRg2bJjW\n6ytXVmL1aqBTJ9PLL+R1fHwAWrYETpzQf/78+fMz9UORIsCaNUpUrGh6+fm+3rZNiTNngJ07DZ+v\nVCoRHh4OAJn6kje8/WkIIdHR0cTGxoZERUURQggZOnQo+emnn4xy2enShQbhGMPo0YTk6MYiyBqI\ncu0ajTQ1xp3tzRtC8ucn5N078WWTk6pVCVmyJMKoaxYsICQsTBp55OLsWfrZc6IrcI0QQhISCClU\niEZlWzLBwYRs22b4vKxjMWAAIXPmSCeTHKxeTUinTvyuFaKeBZliXF1d4erqipo1awIAOnbsiAsX\nLhjxUKGbaFyXpmratKFmDEtD/ZQGqPytWxtXMN3BAfD1pbZ5S+XFC+DOHaBv3wCjrmvaFDh+XBKR\nZEOXGSbrfZETV1fqCsvDJ8FsePOG5gBq3tzwuVnHIjCQ6gdL5uRJGnwkN4IUe/HixeHm5obrn5xt\nDx8+DC+uxnJQN0d7e6BcOeP6rVtXY7OzVIy1r6tp2tSyb/bISOrZpC36Th8eHvShkJQkjVxyYIx9\nPSuW7h1z8CBQrx7g6GjcdQEBNNunJe+tnDxJzY5yI9grZtGiRejatSuqVauGuLg4jB8/nvO1R47Q\np7Kx2NhQzwJLu9nV9rQXL6hPL58fuaVvoJ48SQOysu43cMHKiioHS12tqGsLNG6c+5ihsbB0xb5z\nJ/dJTNaxKFKETvrOn5dGLql59Ah48oS7Y4iYCFbs1apVQ3R0NC5evIjt27cb5RXDdwYD0Btlzx5+\n15qagwfpD9yQ37426tShKx1LTWMsZAZTvz693hI5fRqoWhUoUMD4a2vWpAri9m3x5ZIaQqhXT6tW\n/K63ZHPMqVN0MmJlgjBQk0WepqfTREdNmvC7vlEj+mOxpIhEtf0wMpIWleCDnR01RVmivTklhc6+\n6tQxnJteG/XrW+6MXd8kxtBYWFnR38mJE+LLJTV371L5y5Thdn7OsbBkxW4qMwxgQsV+4QINzCle\nnN/1JUoA+fMDN2+KK5ccnD0L1K7N/3pLNcdcuEADrfjMWgE6c42Pp3EPloaQ1SlA9yWiosSTRy7U\n97oxTgJZadSItvHhg7hyycEXqdj5eMPkpHZtywo1VyqV+PgRiIujofJ8sdQNVLV9HTDexg5Q05WP\nj+UpuFev6ANJV3AOl7GoXdvyPjdAf5+1anE/P+dYFChAbdSW5hX07h39zj85DMqOyRQ7343TrFji\nLCYuDqhQwXCOFH34+1OPIEvL9ijGDMYSzTHHjlHzU968/Nvw86OKwtISY0VFCVudApZpjjl7lrom\nC/nOhWASxZ6WRu3M2jwEjMHSFHtAQACiooybwWjD2pouUXlMeiUlOTUZN5/fRNT9KOy7sQ8RtyPw\nLPkZALoXcuqUZsbOx8YOWOYG6okT1HVPF1zGIl8+oHJl4OJF0cSSnNRUIDZWezERXQQEBCA1PRUX\nHl7AwZsHcereKZSuFYN9pyzLt9mUZhhAYEoBvly+TAMvnJ2FtVO9Op0Bf/xIiwBbAmfPapSbEOrU\noVXdv/1WeFtCeP7+Of69+i+2XdmG43ePo2i+oiicrzAK2xfG249voXqsgkMeB3gUqAkr9+9QslRL\nCJlP1K8P9O5Nc4iYwtuAD+fO0ULlQlFPZIRODORCpaJJ7rjsqSS8SsCyc8tw7O4xXEy6iLLOZVHC\nsQSSU5Px5sM7qKo/gPviQmhVuRVaV2qNxmUbw0phvjfAyZPAkCGm698kiv38eWE2ZjUODtSsERdn\n3KzAVCiVSkRFBWDECOFt+fsDM2cKb4cvL96/wJRjUxAeG45m5Zuhu093/NnhTzjaZY9CIYTg3qt7\nmLw2AjENJ6PyomEYXGsw3N+6I6RZiNH9FitG89PHx9P0ruZORgYQE6P/flcqlZxm7bVqUbOOpcDF\nvq56pMJvp3/Dnht70LNaT3S074j9I/fnuo/cPTIw6fcLuE72YOj+obCxssHMZjMRVCFIwk/Aj7Q0\nuiewebPpZDDJI+/8eaqYxMCSzDFv39KshmIELPj7Uy8Tud090zLSsCRqCdyXuCMlPQU3htzAtm+2\noVPVTrl+jACgUChQxrkMSEwYppWOxvp263H87nF8t/M7nE44zUsGS7Kz37pF6wgUKSK8rVq1LMtZ\nQJ99/X3qewzZNwTNNzSHZ1FP3PzfTcwNngu/En5a76Ma1a2QcrsGJgVMwsX+FzG+4XgM2jsIQeuD\noHqkkviTGEdcHE2xXbiw6WQwiWK/cEGcGTtgWTd73rwB8PenkbNCKVqULnFv3RLeFleS3iYhIDwA\nf1/5G4e6H8KyVstQNH9RTteeOgU0bKhAXbe62P7tdizovwAdtnTAqIOj8D71vVFyWJKdncvqlOt+\ng6cn8OCB5eRn1+XWe+HhBVRfUR3Pkp/h8sDLGNtgLJzzUrusrrHw99dEoCoUCnT07Ij4gfEIrRKK\nwHWBWBu7VqJPYTymtq8DJlDsaWn0iebnJ057luQGZqzrlyGqV5cv3DrqfhRqrqyJoPJBONzjMHxc\nfDhf+/gx8PRp9nqPHTw7IK5/HO6+uosGfzTAo7fcXXwaNLCcGbuYq1Nra9rWuXPitCclr17RfE5V\nq2Z/f0nUEoRsCMGPjX7Epg6bUNC+IKf2qlenE8Ks2FrbYlCtQVD2VOKXk7+g3+5++JBmeof3yEga\ncWpKZFfsV6/SwCS+QSo58fKi5g1LCLHfu1cp2PUrK3Ip9vUX16P1ptZY3GIxJgVMMnrTKiaGPsiz\nbnYqlUoUzV8UWzpuwVdVvkK9NfVw49kNTu1VqUJnrU+fGiWGSeAyYzfGp99STI/nztHvXL06JYRg\nYsRELIxaiLPfnUUX7y5ar9M1Fr6+1CMoPT33Ma9iXojuG41nyc/QOLwxnr9/LtKn4EdMjHgPc77I\nrtjFnMEA9Mbx8zP/WQwhNMeLmDN2tZ1dSlaeX4kfI36EMkyJr9y/4tVGbCz9YWpDoVBgYuOJGFt/\nLBqFN0L0fcO1/xQKoFo183f9I0RcsyNgOSvUrGaYDJKBIfuGYPf13TjR6wTKFTQynSvoPkWJErT6\nljYK2BXA1q+3olGZRmi2rpnJlPu7d3Sl4u5uku4zkV2xi32jA5ZxsyckALa2AShdWrw21TN2qTZQ\nw2PD8fPxn3GkxxF4FvU0fIEOtCn2nLbUvtX7Ynnr5Wi5qSXOPTD8lPb1pe2aM7duUc8tQ7UujfHp\nV+8pmXuOJLVizyAZ6LmjJ1SPVYjoGYFi+fUPhr6xMLRCVSgUmNVsFoIqBKHpuqaZMRRyolLRFNNa\n01K/fw/cuEGjrdatA6ZPB/r3pxnSfHz4Z0rTguzujufPA6Gh4rZZq5b5FzsWmjNDGy4uNMz+7l3q\nLywmG+M2YsLRCTjS4wgqFqooqK3YWGDcOMPnta3SFivbrETbzW1xLOwYKhWupPNcX1/zj0YUy603\nK25u9N/ERM3/zQ1C6P2+cCEw8uBIJLxKwP6u+2FvyyOdaRbUK1R99UMVCgVmNJ0BBRRouq4pjvY8\nikL2hQT1y4mUFCAxEY+3JKC/QyLwawKdzSUm0n8TEqhbXKlS9Itzc6PBPGqFrq6oIhKyKvb0dLp8\nFmvjVE2tWsD//kdvKDEVp5hERQFFiyqhrnUpFupZjJiKfe+NvRh1aBQOdz8M9yLC1pTv3tEHj4dH\n9vd1+W6HuofiybsnCNkYglO9T6G4g/Yscb6+wNy5gkSTHK6KnasfO0Dvb7Wd3VwVe0IC/ffvB/Nw\n6OYhnOx9krNS1zcW1asDP/9suA2FQoFfm/6KlPQUdNzSEQe6HYCttZGVXbLy8SN1R8qpqLP+vXwJ\nlCwJjxQ3uJdyBV640Zu+eXONIi9aVDYFJativ3aNzjKFRpzmpEwZGghizrOY6GhRV1qZqBV7hw7i\ntHflyRWE7QjDjk474FVMuMO93qWpDvpW74ukt0losbEFjoUdQwG73Dvtnp40s+eHD6bLx2GICxeA\n4cPFb1dtjhHrOxebCxeAkkFbMO/MXJzqfSrTlVEofn50Y5JL1LFCocDsoNkI/SsUg/YOwvLWy6HQ\nplTT0rQr7az/f/aMpqF1ddUoaVphW/PaxQWwskL3OsDs2UBlE7s7yqrYpbCvA/Qh6OurCQwwNwih\nsv31V4Dobfv7A0uWiNPWi/cv0PbPtpjZbCbquYnjr6Vr49TQDPXHRj8i8U0iev3bC9u+3pbrR2ln\nB1SqRCNQpbinhEIId0cBY/Pm+PlRM4e5sufiGVyrMBinOh9EaSfjzAv6xqJQIRro9d9/NG+OIawJ\nsLn2bPT7vRX2qnqjlb1PbqX9+DGdSatNI25udKbYoIHmvRIlqK+pAdLT6UTGh7snsGTIqtjF9ojJ\nio8PVZ5SzIqFcv8+nbG6uIjfdtYNVCGrvLSMNHyz7Ru0qdwGvfx6iSafPo8YfSgUCiwMWYhG4Y0w\nJ3IORtUblesc9QaqOSr2O3fo/gffegP6qFaN3uvmyLPkZ9j48VsMKrcSvsV5fPEGUAcqVa6YQdOb\n5jSNZH2dlASHwoWxpqQLjnz8E3f8n6Osd0O65FHPtEuUML4Arw6uX6fNieXKLQTZFfvEidK07eMD\n7NsnTdtCiYuj8hljS+WKejIh1Aw17sg4KKDArKBZ4gkHqni7ds39PpexsLOxw9avt6LWylqoWbIm\nGpfNng60WjXz9YwxZuPU2PuiVCm6V/f4sWGPGznJIBno/k932N34Br278XONVSqVCGjcmNYC1GIa\nmXs+AQ6HEoCwB9QHMutM282NLmfU/y9VCsiTB3YAnO6dQp0tHXDu+yVwLeAq7gf/RGys+PuHfJFF\nsZ9OOI06peohNlbaGbspk2LpQ63YpUCh0Mza+Sr2gzcPYrNqM2L7x8LGSrxbQoylaWmn0ljXbh26\nbO+Cc33PoYRjicxjvr7AP/+IIKgESOERo0ah0KxQmzWTpg8+zDg5Ay+SX+PD3l9QaZOOkwihNmtd\n9uwbN4Dnz2l5tKwK29UVCAlBkq8bFv7tig1KV6M2V+qXro//1f4fum7viqM9jsLayrBpxVj4rk6l\nQEGItB6xCoUCZeeXxZbAGHzT1lmygrwpKXRT9sUL89tM69IFCAkBevSQpv2JE+nvZepU4699/O4x\n/Jb7YV3oOjQtL7CkVQ6uXqWmMTHKF/587Gco7yhxuMfhzMjXZ89oFfuXL80vhW/z5tRTq3Vradr/\n3/+oJ5QYmULFQHlHic7bOuEP98P4YwzBX7O1KG31azu73DPtrK9dXWkCei08eUL3Vl68MN70mJ6R\njuYbmqNR6UaYFDBJhE+dneBg+r2IZQ5WKBTgq55lmbG3rNQSQw/1h5//ZgDSuPvY2dGN6suXTR/O\nm5OLF4ExY6Rr398fWLnS+OsIIej9b2909+kuulIHxJ3BTGg4AQdvHsTCswsxrM4wADR7nrMzcPs2\nTd9sTki9LPfxMUEitFevtM60U+/eRqlLJ5Hw1gYZGfXhndcVWJhFWTdqlF1xCygfljX5nbHfubWV\nNda3W4/qK6qjSbkmaFSmEW85ckII9dgxlxm7LIp9dtBslDlbC8V8/gDQW7J+1MtTc1LsHz7Qm9DD\nQxobO0A/t4pH5tIl0Uvw6N0jbG+yXXSZAP2K3dixsLayxtrQtaizug6aV2ieGQmrziFiTor90SPq\nRVeyJLfz+dwXPj7A0qXGy6aTd+90b0Kq/9LTc8+w69bFQrdEvOvQHhO7rcTg8QVQsSIwbBg/MbiM\nhTqdBJ+b/z9bAAAgAElEQVTvvKRjSaxuuxrdtneDaoAKTnmd+Amag4cPqXLn+p1LjSyK3d7WHh5X\nNiHCJxAJr4Lg5iSNT6JasZsTV67QG9DOTro+ypalS9OXL7nHCFx/dh2TlZMR2ScSeaylKT8VGwsM\nGiReexUKVcD0wOno8U8PRPaJhK21baZnTPv24vUjFJWKFgGRMhbFy4uautLSOKSBTk7WKGpdtu2U\nlOymEDc3Wr0mNFSjyJ2dc32ondd2YumBm7j43UUgjwPi4qT/Lry96Rjz7adlpZZoWaklRh4ciVVt\nV4kik3qFZi4BkrJ5xdyN8kavzv9Dv939sKfLHu3BAgLx8QHmzBG9WUHExdEZBsC/zqchrKzoD/3S\nJW55oDNIBvrs7IOfGv2kN2xfKPpm7HzHoq9/X+y4ugPTTkzDlIAp8PUFwsN5iygJasXOFT5jod5b\nvHEpBR6Oelz+EhLobNzVNbvirlaNbgColXahQkZrpafJT9F/d3/81fEvOORxyIzXELJZzmUsvL2B\n7QIXmbOCZsF7mTcO3jyI5hWaC2sM5rVxCsik2F+/ppseM1qNRZ01NbE+bj16VBN/J1G9RDOn1AJS\nesRkRT2L4aLYl0UvQ3pGOgbXGiyZPElJtJixq8ieZQqFAqvbrka136shtEoofH39zM7lUaXSXTnI\nKD5+pEEQOny1o+8lwKHmK8C1ZPaZtkyh7IP2DkIX7y5oWKYhAJrVMH9+capF6cPbG5gyRVgbBewK\nYEXrFfh+1/dQDVBprdpkDLGx4ufAEoIsiv3SJRoCnjePLdZ8tQYhG0IQVD4om+uaGJQoQZX6o0fS\nBIbwIS5O47kglY0d0Ch2Q9x9eReTlJNwotcJSVy+1KhnMLr0iZCxKOFYAjOazUDfXX0R2fssXr2y\nxrNnpi1FlhWVCvjuOwMnqUPZExOhPHAAAY6OukPZsyrtSpWAJk0ANzes3uqGZzYumP6r/C5BO6/t\nRMzDGIR/FZ75nhiTGC73RZUq9CHy/j0NAuNLcMVgNC3fFGMOj8GyVsv4NwR6v0+eLKgJUZFFsWdd\nmvqX8Eff6n0xcO9AbP9mu6gmmaz+veag2AmhKwi5ZuxbtxqSh+D73d9jZN2R8Cjqof9kgcTGakxQ\nUtDLtxfWx63HknOLUK3aMMTGAk3Fd+wxmvR04Gp8OrwLPgTO6LBnJyZmD2XPm5fu+GcNZXdzozex\nnlD28g+AIzy8oYTyJuUNBu8djLWha7Ml95JrdWprSz3grlwR7igxp/kceC/zxrde3yKgbACvNt69\no18plzQHciGKYk9PT0eNGjXg6uqKXbt25Tqe0+b4U6Of4Pu7L3Ze28m7eIMu1Iq9uXCzmWAePaI/\ndPVOuVSzdUAzY9dnhtqo2ogn755oDc8Xm0uX9AfPCB0LhUKB5a2Xo97qemjl2w4qVRl5FHuGllD2\nLP/PuJ2AZ++TkCewcHbvEVdXnaHsATxFMVVqgYnKiQgsF4gm5Zpkez8uDmjbVljbXO8L9f0uVLE7\n53XGwpCFGLBnAC72v8jLkeDyZVpYQ4xaxmIhiigLFiyAp6cn3rx5o/W4SpXd/pTXJi+WtlqKXv/2\nQrPyzZA/T34xxABAFXtEhGjNCUK9cSqHvb9IEbos1ZVa4OWHlxhzaAx2dNohLIUpR1Qq/i5vXKlc\nuDKG1xmOP08Ngs2lXRAcI0FI7lD2nLbtB1lC2bP++fkBrq44Eu+G5btL4Z890ngaZaVMGepa/vw5\n3fuUg/MPzmOzajMuDbyU69jFi8CPP8ojB1fTIxdC3UOxJnYNZp+ejfENxxt9vUqVu7arqRGs2BMT\nE7F3715MmDABc3UkyNbmJRBYLhANSjfA1ONTMaPZDKFiZOLjAyxYIFpzgsi5NJXSxg5obnZtiv3H\noz/iK/evUKuUiLX5dJCaShMieeopuiTWWIyuPxqro/xx/Mk2AF/rPtFQKHtiIt2ozBnK7uZGf7VZ\n84/oCW0+GwF4GGmC4jsWVlaa77xxY8PnCyUtIw19d/XFrKBZKJIv+w5pcjLNuy+0JBzXsfD2Fi/D\npUKhwKIWi1BjRQ10rtrZ6NJ9ly4Z5wUlB4IV+/Dhw/Hbb7/h9evXOs9JTg7DsmVlAQDOzs7w9fVF\nQEAA5jSfA/dR7nB/646w0DAAmmK26i/X2NfPnilx9Srw8WMA8uQR3p6Q19TWr4RSKU9/3t7Av/8q\nkS9f9uPXnl7D33f/xuWBl2X5/HfuAK6uAciXT/f5aoT2d/rEaXxfuC9muQ7H22g3nDuiBB4/RkC+\nfEBCApRxccCTJwh4/hyws4OyYEGgWDEE+PoCrq5QuroCfn4IaN2avv5UY5GvPEePKtGoEaA2sHC5\nPjY2lnd/RYoo8fffQOPG/K435vXS6KUgdwjcKmtmDurjDg4BqFIFOHVKWH+xn1ycDJ3v7R0AlUrc\nzzey7kh0ntMZvzb9FU2aNOF8/fHjwNSpwvtXKpUI/+S7W1Zo5RwigF27dpGBAwcSQgiJiIggrVu3\nznUOABIYqLuNRWcXkcZ/NCYZGRlCRMmGuzshcXGiNcebatUIiY6Wr78//iCka9fs76Wlp5EaK2qQ\n8Jhw2eT46y9C2rUTscFXrwi5dImQffsIWbGCkIkTCenVi5BmzeiXnT8/eWlrS+67FiEkOJiQPn0I\nmTyZkNWrCTl4kJArVwh580ZEgXRTuTIhKpUsXRFCCFm6lJDvvpO+n0dvH5Eis4qQ+MfxWo+vXk1I\n9+7Sy6EmI4MQJydCnj4Vr82UtBTisdiD/H35b6OuK16ckIQE8eRQI0Q9C5qxnz59Gjt37sTevXvx\n4cMHvH79Gj169MC6deuynadvmTKgxgCEx4Zjo2ojuvl0EyJOJuoNVFMuj9TmCC/hRYg44+0NzJuX\n/b0V51fA3sZekrgBXRhlc8wZyq7NTJKRkTuUvV69bG6A7bq/w7ma3jg3YCEqFzaNe8L799QNr0oV\n+fr08ZEnQGvckXHoUa2HzqLmcnnEqFEo6D2mUtFCRmKQxzoPlrZaip47eiKkYgjy2WpPRJaVp09p\n2pBSpcSRQTTEeroolUqdM/ZVq/RfG5kQSUrOKUlef3gtiixTpxIyZowoTfFGpaKzt6xERERI2mdy\nMiH29oR8/EhfP333lBSdVZRcTLooab85CQ2ls3by7h0h164RcuQIIeHh9Ivp14+Qli1JRLlyhDg7\nU4ErVSIkMJCQnj0JmTCBkN9/J2TPHkIuXiTk+XM6PTPA6NGEtPh5NgleHyzq6s8Yzp0jxMfH+OuE\n3BevXhGSLx8h6em8mzDImYQzpMTsEuTVh1c6zwkMJGT/fuF9GTMW/foRsnCh8D5z8s3Wb8ikiEmc\nzo2IIKRBA/FlIMSEM/ac6PJJNzRzruNaB03LNcUvJ3/Br01/FSyHtzewYoXgZgRhig0Ve3s6iVWv\nFCYqJ+Jrr6/h4yLBVOpTVXZtIezTDyagckQC8OGd7qrsSUlAu3Y0qkgEtyFvb+Du3v/hUuE12HF1\nB9p5tBPhQxqHKVaJBQpQj6hbt6hvt9ikZ6Rj0N5BmNlsptbas2ouXZLfM8Tbm3riiM1vQb/Bb7kf\nwnzDUNa5rN5zzdEjBhAxQKlx48ZorGNrnos5YkazGfBZ5oM+fn1QsZCwO7RqVXqjmRJtX3iAWGtG\nPai9JNKLxGFr/FZcGXTF+Ea0hbLndAF89Yo66Gf1IPH0xIdGzfHdQTcci3cDShTRmSg9QNjHzEXV\nqsCsWbZY9Msi9NnZByEVQ7IFz8iBsTli1Ai9L7y96f0uhWJfE7MGeW3y6jWTPn5MTY9iZDY0Ziy8\nvYENG4T3mZPSTqUxrPYwjDo4Ctu+2ab33EuXpA3E44ssLvX5Obipl3QsidH1RmPEgRHY2XmnoP7K\nlaOeba9eUZdjU3DpknSFNfTh7Q3EqQiWpwzF5IDJKJwvR5x9llB2rbbthATqGK0OZc9alf1TKDtc\nXTOrsudEFQ0kewC2Mtsc3d1pkeMGpQLhX8IfcyPnYkKjCbLKoFKZpqKReiIjdq6SVx9e4aeIn7C3\n6169EeLq2brc+Zm8vWkx84wM8QutjKo3Cp5LPXHk1hG9tQpUKu2lH02NGcVKAcPqDMOqmFXY/99+\nhFQM4d2OlRX1oY6Pp3tspkDb0lQpsR870tNRs+RDbN2+Fv6Pb6Hf02Rg0/DsJhN1ocys+UfKljUq\nlF0fXJfkYo+FvT0N2Ll+nS6la66siTDfMJQqIN8Thu+MXehYeHsDO4XNhbQy9fhUtK7cGv4l9Id3\n8v3c2jBmLAoWpKaou3fpZE5M7G3tMaf5HAzdP1RnyUhCTGOC4oJZKXY7GzvMC56H4QeGo2m5poIi\nJNU75qZQ7G/f0sT7oi6NtYWy5/QiSUpCkFNhFLN5jjKkFqzzJVBFXadO9lB2CWOfTRmsoZ65dqpa\nHt9X/x7jjozDunbrDF8oAk+eUO8IsbNZcqFqVeCXX8Rt88azGwiPDUf8wHiD5166ZLriNmrTo9iK\nHQDaubfDkuglWHF+BQbWHJjr+L17gKOjfFG/xmBWih0AWlVqhUVRi/D7ud8xpPYQ3u2o7Y6mQJ07\nIuekV+dMJCNDZ1V2g6Hs/v6a2XepUpgR+RumLI/B89+3Abr3uiRDpQKGDjV8nhQrF/WPvFMnYHyD\n8XBf4o6ziWdR21WMHLr6EVJcQ+hYuLvTzdOUFPEKuow8OBJj6o+Bi4OLwXNVKvHMjsaOhfo7F5qj\nRhsKhQLzguchaH0QOlftjIL2BbMdN9fZOmCGil2hUGBu87losrYJuvp0RSF7fo/DqlWBf/8VWTiO\nZFuaGgplT0igG5UODrlD2b29Ne8ZCGUHgAdvHmD+2bnwvB+N+Higbl3pP2tOTD1jV4dQONo5Ynrg\ndAzdPxSn+5zOLIAtFWKaI4zFzo7OWK9dE8eX/NDNQ7j85DK2fm0gXSjonCQ+3nQKztsb2LNHuvZ9\nXHwQ6h6Kn4//jHnB2YNETPmdG8LsFDsAeBXzQkfPjphybAoWhPBL/KKesUtadIMQWpMuh9uf9/YE\nNPuYAFROpO/lzUvD1fPlywxlR7Nm2TP/6ajKbgwTjk5AX/++eHy1POLi5Ffsz55RM5S2XDU5kWK/\nQW1+U9OjWg8siV6CzarN6Ooj7Q6XSgVUr87vWjHGQv3ZhSr2tIw0DD8wHLObz4adjeHp/717dCFZ\nsKDBUzlh7Fh4ewMzxEs1pZWpTabCc4kn+lfvjypFNNFnhjKYmhKzVOwAMCVgCjyX0sHkkzvc5dMK\nUlDRjdevdXuOqJW5tXWuqMjjaIxGvdxQpsMnxa2uyk6TxvAURj/nH5zH/v/249rga1itEi/znTGY\nyjtCTcWK1D3+7Vs65FYKK8xtPhddt3dFO492nCIJ+XLpEhAWJlnzBhHL9Lji/AoUy18MX1Xhlk7b\n1LNWDw/xzVA5KZa/GMY2GItRh0ZhV2dNWnKVChg+XJo+haL4FOEkXQcKBfh2MTdyLo7cPoI9Xfit\ntZo0AcaPB4KCtBx8+1a3PVtbKHtWM0nW/xfIbcguXhw4d06+jTRCCBqHN0Z3n+7oW70vDh8Gpk4F\njh2Tp381ixfTm335cnn7zYq/P/D77zT1uZpvtn4D72Le+KnxT5L0mZFBZ6337ok3czWWf/4BVq8G\ndu/m38bLDy9RZXEVHOh2AL7FuRXw/PVX6h3722/8+xWKpyewebO0/uQpaSnwWuqFpa2WonmF5khN\npT/9Z89EWWxrRYjuNNsZOwAMrjUYy84tw4H/DiC4YjD3Cz9VZW/vlICMNQnAmYTcUZJZq7Kr/2rV\nAjp00Lzv5GT09FPtHSFn7oi/r/yNVymv0NuvNwBuRTekwBw2k9SfPatin9lsJmqsrIHefr0lcX+8\ncwdwdjadUgfEmbFPOz4Nbau05azUATrWwUb8NKVA/Z1LqdjtbOzwW9BvGHFgBGL7x+LGDRu4uUmn\n1IVi1oo9j3UezA6ajREHR+Bi+YvUl/TDB6qctbn9qd/7VJW9vbUb7qS7AeVcaQHONm00SptHVXYu\nqDcPtTUthV35Q9oHjD40Gqvbrs6sYeriQj0aHzyQ9wGjUgHffsvtXKl8+rVFHZcrWA59/ftiwtEJ\nCA8NF71PoeYIMcaiXDk6qXj9Wusi0iBq90ZtBTT0cekSMErEglx8xkLMohv6CHUPxcKohVh1YRUK\n3uxv8kmMPsxLsatD2bMo6bYJCShy+jFeLCyPos/e0zu3ZMnsM20vLyAkRPO6SBFAocDd08DwYUCU\nyD6++pB71rrgzAJUc6mGwHKB2d5X3+xyKXZ1sIacGf604e0N7N+f+/3xDcejyuIqOPfgHGqUrCFq\nn6a2MwN0q0cdlMdn03zM4TEYWXckijtw35BKTQVu3KB2blPi7S2P+U/t/hiyIQTdXnaGj4+Jwto5\nIJ9iT03VHsqe9fXz5zSAJouJRFGpEkr5j0O3i9OwZdhpOJWpzDl+2MuL+pRLEXKsC31LQrFnqI/e\nPsJvp39DZJ/IXMfUij2EfwCvUdy5Q4M1Chc2eCoA6fLmqFM25zRDFbArgJ8DfsbwA8NxPOy4qEXU\nVSq6GOSLWGOhXq0Yq9gjbkcg5mEMNnfYbNR1168DpUvTqF+x4DMWcs3YAcC3uC9aVW6FHfumYVaQ\nCTcWDCCPYi9Viq4Tc4aylysHNGqkmWm7uGgNZS8LoMyuq/j5xkrMKTeHc7dOTlTR3L4NVKgg3sfR\nx6VLQDdx0sob5KeIn9DTtycqFa6U65i3t7ybp3Ln49ZFiRJUqScl0f9npbdfbyyJXoJtl7fhay89\nZfSMRKWim/Smho+CS89Ix/ADwzEraBby2uiPk8iJOeypADQjxsuX1PNYjn2OaU2mwTXSG87l+gOQ\nSbEYiTzz2MhIavdOTKT/37oVmDuX+gp17AjUrk3NK3ryk0xtMhVrY9fixrMbRnUtZ6ZHQ7kjcpaF\nE0JsUiz+vfYvfmqk3dNDzlkMYLxiF3MssqJQ0BVTXFzuY9ZW1pgXPA9jDo/Bh7QPovSXkkInDkJq\nfYo1Fnzu9TUxa1DArgC+9jT+QSdFylo+Y2FlRVfncv3O82WUgHXUCCy+PlqeDnkgj2IvXRqw5Z/3\nBQBcHFwwpv4YjDpk3E5NzqAVKZErdwQhBCMOjMDERhPhnNdZ6zleXjQSMS1NWlnUmMuMHdCYY7TR\npFwT+BX3w9xI7YXXjeXKFaB8eel8qI1Bfa9z9ZB79eEVJionYn7IfF6mKXMq4uzjI9/vXKUCfN+P\nQExSDI7ePipPp0Yik+VZHIbWHor4x/E4ePMg52vkzBljaGkqli11+5XteJL8BP1q9NN5Tv78dBF0\nw7gFDm+MVexSZrnUp9gBmv1xTuQcPHzzUHBfYmycijUWJUrQ/aTHj7mdP+3ENLSs1NJg9kZdSDFj\n5zsWcq5Q4+IA36p5MTtoNobtH4a0DJlmT0ZgUYrdzsYOc4PnYviB4UhNT+V0jbe3/h+5mMhRQed9\n6nuMOjQK84Pna00lmhW5bvbkZPlrferDx0d/ZZ0KhSrgO//vMP6ocMO4OXjEqMlaB9QQ/z3/D3/E\n/IHpgdN59fXmDd3HkKK4Bx/kVuw+PkB7j/YonK8wVp5fKU/HRmBRih0A2lRug1KOpbDs3DJO53t4\nUI+N9++llQsAYmOpu7wuxLClzo2cC7/ifnqT/6uR62aPj6dK3Rhrm1Q2doCaoW7coN6zupjQcAIO\n3jyIqPtRgvoSQ7GLORbVqnErFzfy4EiMqjfKKPfGrMTF0XEWOwM037HImhtKatSKXaFQYH7wfEw+\nNhkv3r+QvmMjsDjFrvYlnXp8Kp4mPzV4fp48VOnIYY4xpNiFcv/1fcw9Mxezm8/mdL5cij0uzrzK\ng+XNSx2url7VfU4BuwL4temvGLx3MDJIBu++zGnGDgB+fvQ+1Mf+//bj8pPLGF6Hf6ITqe91Yylc\nmEaBJiRI209GRvbvvFrxamjv0R5Tjk2RtmMjsTjFDtDsj52rdsZPEdxyf/j6Gr7ZhaJOPaPPHCHU\nlvrD4R/Qv0Z/lC9YntP5cip2YzdOpa7/asjODgDdfLrB2soaa2PX8urjxQsaL1emDK/LMxFzLAzd\n6x/TP2Lo/qGYHzyfU/ZGXUil2IVWkpL6fr9zh7pUZnWr/DngZ2xUbUT8Y8NFSeTCIhU7QLM//nPl\nH5x/cN7guXIodpWKRv4JdP7RyYm7J6C8o8S4BuM4X1OxIq3k9PatNDKpMSePGDWG7OwAzf64qMUi\njD86Hq8+vDK6D5WKmiPkCn7jgqcncPMmzbyhjfln5qNSoUpoVbmVoH7MbcYOyKPYtd3rRfMXxeTG\nkzF432DeSbvExoxuSeMoaF8QvzT9BYP2DjK4lJZDsXO50fnaD9My0jBo7yDMDZ4LhzwOnK+zsaH+\n1fESTiQI4afYpbSxA7p92XNSo2QNtK7cmtdSWiwzjJhjYWcHVKqk3fR4//V9zDo1C/ND5gvqIy2N\n3lNSmKCEjIWpFDsA9K/RHy8/vMSfl/6UVgCOWKxiB4Aw3zAoFAr8EfOH3vPUP/IM/qZUg0g5g1kc\ntRguDi68gkik9gp68IA+QFwMV1CTFS6mGDXTA6djfdx6XHps3EbMxYvmt1IBdE9kxhweg++rf4+K\nhYS5sly/ToPHHR0FNSM6plTs1lbWWNJyCUYfGo03KW+kFYIDFq3YrRRWWNJyCcYfHY/n75/rPM/Z\nmeYFu3lTOlm4KHY+9sMHbx5g2vFpWNxiMa8gEj8/ICbG6Ms4w1e5SW1jd3Wl5gguPt3F8hfDzwE/\nY8CeAUZtpF64IE4RZ7HHQptij7gdgRN3T2B8Q+EunlJOYoSMhacn8N9/us1QYqDPUaCeWz0EVQjC\nz8d/lk4Ajli0YgcA/xL++Nrza0w4OkHveVKaY9LSpMtsOPrQaHxf/ftsJbmMoXp14LzhbQjemKN9\nHaA+3cbM2vvV6IfU9FSDqz81Hz/SBHPm5A2kJqdnTEpaCvrv6Y9FLRYZZcrTRWyseX7uvHmBypWl\nW6Gqs6JUyp2aKZOZzWZibexak2+kWrxiB2gemR1Xd+BM4hmd50ip2G/coFF/hvJgG2s/PHr7KE7e\nO4kJDfU/tPTh60sfOqnc4rmMhq9il9rGDnC3swN09fd7698x/uh4PHn3xOD58fHUpTJ/foFCQvyx\nyGl6nHFyBjyLeuIrd27l7gwh5Yxd6FhIOZGJj6dxMfp894vlL4YpAVPQb3c/QW60QvksFHtB+4KY\nFzwPfXf1xcd07VEpUip2KW705NRkfL/reyxusRj58/DXHo6ONHHmlSsiCpcFc/Nhz4oxM3aApmTt\n6t0VYw6PMXju+fP8i1dLTcGCNF/RrVvA9WfXsShqERaGLBSlbULM0yNGjZSKneskpl+NfiAgWH7O\ndDUiPwvFDgDfen2LMk5lMOvULK3HzUGxG2M/nKycjBola6BNFQGJvj8h1c2ekkL3LfgUWpDaxg4Y\nr9gB6kZ7+NZhKO8o9Z4nln0dkGYsfH2BmBiC/rv748dGP8LNyU2Udh8+pMq9ZElRmsuF0LEwB8Vu\npbDCyjYrMVE5Efdf35dGGEMyCLk4ISEBTZo0gZeXF6pWrYqFC8WZFfBBoVBgaaulWHB2Aa48yT09\nLV2a5jThmiDJGMSewZx/cB5rL67FwhbijKdUN/ulS9RX3hwyG2rDy4tGnxpjhnK0c8SSlkvw3c7v\n8O7jO53nmfOMHaD347q4tXj54SUG1xosWrvqe13OWrrG4ONDs5pKsYF6/jz337lnUU8MrDkQg/eJ\nN/bGIEix29raYt68eYiPj8eZM2ewZMkSXJFqzc+B0k6lMbnxZPTd1TeXfUuhoF8KlzwaxkAI9Trh\n8oVzsR+mpqeiz84++C3oNxTLX0y4gKAzywsXRGkqG1FR2YtGG4McNvb8+WlKXWNd4NpWaYs6rnV0\nJglLTaUPNbEe5lKMhatnIg6SMfjjqz8MJoszBqnNMELHwt5etx+/EFJTqe6oYURVxfENxuPq06vY\nfmW7uMJwQJBiL168OHw/fcsODg7w8PDAgwcPRBGMLwNqDgABwdLopbmOSWGOSUqim1Ri1RadfXo2\nijsUR3ef7uI0COolERcnfm52IYpdLmrVAs6eNf66hS0WYmv8Vhy/ezzXsStX6L6FuflxqyGEYMPr\n72AXOwTViou7AWLO9nU1UqxQL12ilZqM+c7tbOywss1KDN47mNOGvJiI9ii/c+cOYmJiULt27VzH\nwsLCULZsWQCAs7MzfH19M21p6ie0WK+PHzuO/kX6Y8SxEWhWvhmSLiVlHvf1BdavV6JmTfH6W7dO\niTJlAIXC8PkBAQF6j8c8jMHMjTOxvPXyTJ91scanZMkAXL0KPH0qTnsBAQGIigLq1VNCqRRvPMV+\nXbCgEjt3AgMGGH/9slbL0Hl2Z6xuuxohQSGZx/fvB/z9xZVXjRjt7b6+G2/xFFaRY/HPP0oULCje\neEZGKtGyJQCI+/mz2taVSqWg9goUAM6fF1e+a9cCUKuW8den3UpDYzTG97u/x/ZvtuPYp3qV2s5X\nKpUIDw8HgEx9yRsiAm/evCHVq1cn//zzT65jInVhNMuilxH/5f4kJS0l873YWEI8PcXt55dfCBk5\nUng77z6+Ix6LPciGixuEN6aFTp0ICQ8Xr71XrwjJn5+Qjx/Fa1MKLlwgxMOD//Vd/+5Khuwdku29\nwYMJmT1boGAScfvFbVJkVhFy6dElEhBAyMGD4rX95g0h+fIRkpoqXptSEBlJiL+/uG327k3IsmX8\nrv2Q+oH4LPMhay6sMeo6IbpTsFdMamoqOnTogG7duiE0NFRoc6LRr3o/lHAokS0HiIcHrU8pZm52\nY5amOWdnWfnh8A/U3c6nqziC5cDfX9zl6blz9HPzTXqmbyzEpGpVWgTklfE5vgBQk8y/1/7F7uu7\nM98T0yMGEG8s0jPSEbYjDKPrjYZXMS/RTY9xcYb9uIUixlj4+FBzWUqKcHnUCDE72tnYYUO7DRhz\neC8BojoAABeJSURBVAxuv7gtnlB6EKTYCSHo06cPPD09MWzYMLFkEgWFQoHVbVdjTcwanLh7AgDN\nze7pKW6IfVSUcO+IfTf2Yee1nVjaKve+gFhUry7uBqol2NcB+uDx9aUPIj4Usi+Eje03os/OPkh8\nnYj0dLqJ5ucnrpxiMOXYFNha22Jk3ZEAqIxiPszPnLGM7zxfPqBCBfE2UN+8oTEBQpKeebt4Y2z9\nseixo4c8pfR4z/UJISdOnCAKhYJUq1aN+Pr6El9fX7Jv3z7RlhNisPPqTlJmXhny5N0TQggh//sf\nIbNmidN2YiIhhQsTkpHBv40Hrx+QknNKkojbEeIIpYMXLwhxcCAkLU2c9tq1I2TzZnHakpoRIwiZ\nPl1YG9OPTycN1zQkF1WppEIFceQSkwP/HSAl55QkSW+SMt+7cYMQNzfx+mjfnpCNG8VrT0p69iRk\n+XJx2lIqCalbV3g76RnppNm6ZmTs4bGczheiOwXN2Bs0aICMjAzExsYiJiYGMTExCAkJEeeJIxJt\nqrTBt1W/RadtnZCWkYb69YGTJ8Vp+9QpoF49/j69H9M/ouPWjuhXvR8CygaII5QOnJ1pBsZr18Rp\nLyoK0LJPbpbUrk3lFcLYBmNhZ2OHiUd+Njv/9QdvHqDnjp7Y0G4DXBw0aTYrVKA5be7dE94HIfR3\nU7++8LbkQEzPGLFWp1YKK2xqvwmbVJvwz5V/hDeory9JWzcTfgn8BVYKK4w7Mg7161OFLEY+/FOn\ngAYNuJ+f0374v33/Q7H8xfBjox+FC8MBsezs9+9T+6WQjXu5bOyAxuVRyHdupbDC+nbrcfj5auSr\ntkc84SBsLFLTU9H5784YWGMgmpRrku2YQgHRJjI3b1KzVunSwtvSh1j3hTkqdoAW5dj29Tb0290P\n156KNMvSwheh2K2trLG5w2Zsu7wNJ1/8BUdHcWaup07xn8GsPL8Sx+4ew9rQtbBSyPM1iGVnV9/o\n5hp9mJMyZWisQWKisHaKOxRHxfN/Ywd64WKSyJFuPCCEYODegXDI46AzHW+DBuIodvW9binfua8v\nzb6pr6A5V8TeT6pZqiamB05H+y3t8fajROXNeBtxOCJDF5y58OACKTKrCGnRJ5qsXCmsLbXr1/v3\nxl974u4JUnRWUXL1yVVhQhjJ4cOE1K8vvJ2xYwmZMkV4O3LSujUhW7cKayMlhRBHR0JWn/mTlJ5X\nmjx4/UAc4Xgy/fh04ve7H3mT8kbnOVFRhHh7C++rb19CFiwQ3o6ceHkRcu6csDYePiSkYEFh+2i6\n6LuzL2mxoQX5mKbdZ1iI7vwiZuxq/Er4YVWbVThVpg32RAlLfXD2LPU6yJvXuOtiHsag/V/tsb7d\net451vlSpw51f3unOwUKJyzFIyYrtWvzi0DNytmzNN9379rfoq9/X7T9sy2SU5PFEdBINqk2Yfn5\n5djdZbfeHOu+vtTF9+VLYf0JWZ2aigYNgE/xQLyRcnW6pOUSWFtZo9e/vURP8ftFKXYA+Mr9K/zg\nPxO7nYNx5+Ud3u3wudHX/bsOrTa1wrJWyxBcMZh333zJn5/a2U+c4N9Gejp1HaxZU5gsctrYAfrj\nFLqBevQoEBhI/z+h4QR4FfVC+7/a432qsMAIY8fiyK0jGLZ/GHZ33o2SjvrTLNra0s8eGclfvufP\n6QasHOmZxbwvmjYFjhwR1oaUkxhba1ts6bgFd1/dxfADw0UthP3FKXYAGNuiB2yjRyMwPAhJb5N4\ntXHypHEbp3de3sGog6Pwa9Nf0cGzA68+xUDozX7tGlC0KFC4sHgyyUHNmnR/QUi+nKyKXaFQYFXb\nVSicrzBab26tNxOkmOy5vged/u6ErV9vhbcLN8dqoXb2yEiq3KQMTJKCJk3o5xZSZEZq7y97W3vs\n6rwLyjtKTD8xXbR2v0jFbmUFNMk3BHXz9kKDNQ2M3p1OT6fL8nr1uJ2veqRCQHgAfurxE3r69uQh\nsXgIVexi3ehZc4PIQcGCNId4PM+KZcnJ1Msi68PcxsoG60LXobRTabTY2IJ3EWOuY7Elfgt67+yN\n3Z13o3HZxpzbF6rYjfX+EoKY90WRIjS7J9+VWkYGEB0tfHVqCOe8ztjfdT+c7JxEa/OLVOwAvVFd\nro/HhIYT0Di8MY7d4W6MU6mokuAya917Yy8C1wXil6a/YEjtIQIkFodatWjB32fP+F1/9Kh8P3Kx\nqVuXv4I7dYraqx1ymLOtrayxuu1quBdxR9N1TZHwKkG4oFpYeX4lhu0fhkPdD6G2q3FP1jp16EOJ\nb4i9Jfmv56RZM/4TmXPnaMnLYuJkz9ZLCccSouqHL1axq/3Ze/n1wqYOm/DNtm/wR8wfnOxcXMww\nhBAsOrsIfXb2wb+d/kUX7y6y25W1kScPlZ2PKOnpwL59QKtWwuUwxVi0aAHs4emCfvQoXe1oQ10v\ntYNHB9RcWRP7buwzqm19Y/E65TW6/9Mdc8/MRUTPCPi4GF9g1tGRbvrycXX9+JFeV6eO8dfyQez7\nomlT4PBhftfu2QO0EV7AzCR8sYq9Zk2aSyI5GQgsFwhlTyXmRM5Bm81tcPflXb3XGto4vfHsBoI3\nBGN1zGqc7n0a9dw42mxkIjCQ3yzm7Fm6UpE6SEUqQkLoxjEfr6Cs9nVtWCms8EODH7D16634fvf3\nGHdkHFLShGWhir4fDf/l/rC3sce5vucEeVHxNcdcuECrZBkq1G6uNGxIPwOf73z3bnEmMSaBt6Mk\nR2Togjd16xISEaF5nZKWQqYfn04KzyxMZp+aTd6nandSd3OjeThykvwxmUxRTiGFZxYmc07PIanp\n5pnf9MIFQipXNv66ceMIGT9efHnkJDCQEC3ZpfXy8iXNs/PhA7fzH799TL7a/BUpPa80WXV+lU4/\nZV1cf3qd9NrRixSdVZRsjRfofP+Jv/4ipG1b46+bPZuQgQNFEcFkNGxISI4UVga5f5/6r5syRbEQ\n3fnFztgBOouJiNC8zmOdB+MbjseZ785AeVcJ17muGLx3MGIeatJB/vcf3WWvUIG+JoTgdMJp9Nvd\nD67zXBH3KA4x/WIwou4IUUuSiUm1atTGbmwk5q5dQOvW0sgkF23b0s9hDMeOUVME19quRfMXxY5O\nO7C5w2ZsVG2E51JPLDhDa/ESHaa+dx/f4ejto+i6vSvqramHMk5lcG3wNXT07GicsDpo0ICuNNPT\njbvu4EGgMfd9WrOEj519714gONjyPIHUKIiuO02sDhQKUf0zxeTMGSAsjOZu1haAcPflXay9uBZ/\nxP6Bj+kfUaFgBby5Vx6K5GLwrP0Ad1/dxc3nN1HArgDCfMPQzacbXAu46uwva2UYU/P119R+2KMH\nt/Pv3KEbrw8fAtbWwvs31VjcukW9mR48oN5RXBg2jCZQGzeOX59Hbx/FJtUmHLp1COkZ6ahfuj7y\n2eaDjZUNFFDgxPETuFfwHqq5VEOoeyj61+iPAnbi2z6qVwdmzqSKjguPHgFVqtDcQPnziy6OVqS4\nL06dAgYPNi5dd2go/Y10laY8AieE6E4LfR6JQ+3a1K/5/HntRWrLOJfBxMYT8WOjH5H4OhE3n99C\nl0G30K7bI9Qu74UyzmVQxqkMyjqXzSxlZymo3R65KvY9e4CWLcVR6qakfHnqzRQdzd1t8+hRYNUq\n/n0GlgtEYLlAEEJw4/kNRN+PRkp6CtIy0pCekQ6v6l74vsP3sLe1598JB7p3B9av567Yt2yhD3+5\nlLpU1KpFH+hPn1IXSEN8+EBX8qtXSy+bZIhhC9KHDF0IYtIkQoYO5XZudDQhFSpIkzdCbq5dI6RU\nKe6fJThYeK4Vc+GHH7jvFTx6RIiTk/mXg+NCUhL9LG/fcju/dm1C9u6VVia5aNmSkC1buJ27b584\nOZWEIkR3ftE2doAutf78k1tE4saNQLdulpPhTh+VKlHXRy6Vhd6+pcvZ5s2ll0sOjLGzb95MvWks\n1daaFRcX6sv/77+Gz715k85yuc7uzZ2gIOrlwoU9eyx/L+mLV+yVKtG84oZ8XdPS6I9ciM3NHPzY\n1SgUwKBBwJw5hs89fJiaLcR0eTPlWNSuDSQlAXf1e7UiLQ2YPx8YPlxaeeQcC7U5xhCbNgHffMO/\npi1fpBqL7t3pw/zBA/3nEWLhbo6f+OIVO0Bn4Rs26D/nyBHqv12pkjwyyUHfvsChQzT7nz527bLc\nQA1tWFvT/QJDs/YdO6jfvqVUiuJCaCh1GkjSkyKJELo6NeXGodgULkyV+4IF+s+7fJmmEqhaVR65\npOKL9opR8+QJVdiJiblDxtX06EE3WP/3P3llk5offgDevwcWLtR+/NEjepNHRQHlyskrm5Ts2AH8\n+itNcKXLO6ZuXWD0aKB9e3llk5qwMOryqmslcuEC0LEjNcd8DmZHNXfuUM+gW7cAJx1pWfr0oRus\nM2fKKppWhOhONmMHzVbYoAH9sWvj3Ttg507g22/llUsOhg6lqxVduWNGjQJ69fq8lDpA7ezW1sCK\nFdqPR0bSh9pXX8krlxwYMsds2gR06fJ5KXWAmlxDQoDly7UfP30a2L8fmDBBVrGkQaQNXJ3I0IUo\nbN5MSOPGhKSn5z62fDkhISHC+4jIGuZqRoSFETJtWu73Dx8mpHRp7l4UxmAOY3HpEiFFihCSmJj7\nWIcO8lUMknss0tKoR1RMTO5jT54QUrw4IfHxsoqUidRjERtLSMmSuaOIU1MJ8fEhZNMmSbs3CiG6\nk83YPxEaSv/t2TO7h8z27cBPPwHTpplGLjkYNQpYvJj676r58AEYMABYtMjy/Zh14eVFP+OQHEn1\nbt2iSdJ69zaJWJJjbQ38/DPdN7l8WfP+s2fUC6ZXL8DT03TySUm1aoC3N91DyMqSJdQE06mTaeQS\nHREfMFqRoQvRePeOkObN6WwtJYWQP/8kxMWFkPPnTS2Z9LRsST/3zp10HKZMISQ01NRSSc/794RU\nqULI338T8uIFIatWEeLvT33dP3fWr6f399mzhDx7RoivL/3cn0Ochj6OHqUr0WXLCPnvP0IePKAr\ntytXTC1ZdoToTrZ5moOUFGpLv3+fbqYeOAD4GJ8p1eJ48oTa2nfupJG4Nja0PqqlZnI0hhMn6OyV\nEDpj7dKF2uDldvUzBbt20Q3DYsWop9DMmZ+fbT0nhAB//00/+8GDNE5jyBDgl19MLVl2hOhOpti1\nkJoKTJ1Kl2ViLknNKVeMPp4/Bx4/BtzdpevD3MYiMhLw8ACcneXv29RjceIETck8cqTplbrcY0EI\nLfdYvjwN2DMnWK4YkbG1pTbIL5VChejfl0TduqaWwHQ0bEj/vkQUCmknMKaCzdgZDAbDDGF+7AwG\ng8HIRLBi379/P9zd3VGpUiXMNIdwLTPGnHLFmBo2FhrYWGhgYyEOghR7eno6Bg8ejP379+Py5cvY\nvHkzrly5IpZsDAaDweCBIBt7ZGQkpkyZgv379wMAZsyYAQAYO3aspgNmY2cwGAyjMZmN/f79+3Bz\nc8t87erqivv37wtpksFgMBgCEeTuyLUcXFhYGMqWLQsAcHZ2hq+vb6avqtqm9iW8zmo/NAd5TPla\n/Z65yGPK17GxsRg2bJjZyGPK1/Pnz/+i9UN4eDgAZOpLvggyxZw5cwaTJ0/ONMX8+uuvsLKywg8/\n/KDpgJliMlGaWVCOKWFjoYGNhQY2FhpMFnmalpaGKlWq4MiRIyhZsiRq1aqFzZs3w8PDQxThGAwG\n40vFZJGnNjY2WLx4MYKDg5Geno4+ffpkU+oMBoPBkB8WeSojbJmpgY2FBjYWGthYaGCRpwwGg8HI\nhM3YGQwGwwxhM3YG4//t3d9LU30cB/C3DwlBBGHkLGcgsh9uMzUsqbsaR0HcKvWiAgOLbiSoiP4E\nt4kXZVA3giQFzYtuIupgJkWgBrFZ4QaO2GL+vNAWGNFsfZ+L8JzHp7TnWWdz7bxfdztu53zPG/ae\nfM4ZIyIFiz2L/nkPt94xCxWzUDELbbDYiYjyDGfsREQ5iDN2IiJSsNiziPNDFbNQMQsVs9AGi52I\nKM9wxk5ElIM4YyciIgWLPYs4P1QxCxWzUDELbbDYiYjyDGfsREQ5iDN2IiJSsNiziPNDFbNQMQsV\ns9AGi52IKM9wxk5ElIM4YyciIgWLPYs4P1QxCxWzUDELbbDYiYjyDGfsREQ5iDN2IiJSsNiziPND\nFbNQMQsVs9AGi52IKM9wxk5ElIM4YyciIgWLPYs4P1QxCxWzUDELbaRd7FevXkVlZSWqq6vR0tKC\njx8/armuvDQxMbHZS8gZzELFLFTMQhtpF3tDQwMmJyfx+vVrmM1meL1eLdeVlxKJxGYvIWcwCxWz\nUDELbaRd7JIk4a+/vr+8vr4e09PTmi2KiIjSp8mMvb+/H01NTVrsKq/FYrHNXkLOYBYqZqFiFtrY\n8HZHSZIwPz//w3aPxwOXywUA6OrqQiAQwP37939+gIICjZZKRKQv6d7u+Fv3sd++fRt9fX14+vQp\ntm7dmu5uiIhIQ1vSfaEsy+jp6cHz589Z6kREOSTt/9hNJhOSySSKiooAAIcOHcKtW7c0XRwREf1/\naV88jUQieP/+PYLBIILB4A+lLssyrFYrTCYTuru7f3uhf5J4PI4jR47AbrfD4XDgxo0bAIClpSVI\nkgSz2YyGhgZd3dqVSqVQW1urXJvRaxaJRAJtbW2orKyEzWbDy5cvdZuF1+uF3W5HVVUVTp8+jS9f\nvugmi7Nnz8JgMKCqqkrZttG5e71emEwmWK1WDA0N/XL/GfnmaSqVwoULFyDLMkKhEO7du4dwOJyJ\nQ+WkwsJCXLt2DZOTkxgfH8fNmzcRDofh8/kgSRKmpqbgdDrh8/k2e6lZ09vbC5vNplxM12sWFy9e\nRFNTE8LhMN68eQOr1arLLGKxGPr6+hAIBPD27VukUin4/X7dZNHR0QFZltdsW+/cQ6EQBgcHEQqF\nIMsyOjs78e3bt40PIDJgdHRUNDY2Ko+9Xq/wer2ZONQf4dixY+LJkyfCYrGI+fl5IYQQc3NzwmKx\nbPLKsiMejwun0ylGRkZEc3OzEELoMotEIiHKy8t/2K7HLBYXF4XZbBZLS0tiZWVFNDc3i6GhIV1l\nEY1GhcPhUB6vd+4ej0f4fD7leY2NjWJsbGzDfWfkP/aZmRmUlZUpj41GI2ZmZjJxqJwXi8UQDAZR\nX1+PhYUFGAwGAIDBYMDCwsImry47Ll++jJ6eHuULbQB0mUU0GsWuXbvQ0dGB/fv34/z58/j06ZMu\nsygqKsKVK1ewd+9e7NmzBzt27IAkSbrMYtV65z47Owuj0ag877/0aUaKnfeuf7e8vIzW1lb09vZi\n+/bta/5WUFCgi5wePnyI4uJi1NbWrntPrl6y+Pr1KwKBADo7OxEIBLBt27YfRg16yeLdu3e4fv06\nYrEYZmdnsby8jLt37655jl6y+JlfnfuvcslIsZeWliIejyuP4/H4mk8cPVhZWUFrayva29tx/Phx\nAN8/hVe/8DU3N4fi4uLNXGJWjI6O4sGDBygvL8epU6cwMjKC9vZ2XWZhNBphNBpx4MABAEBbWxsC\ngQBKSkp0l8WrV69w+PBh7Ny5E1u2bEFLSwvGxsZ0mcWq9d4T/+7T6elplJaWbrivjBR7XV0dIpEI\nYrEYkskkBgcH4Xa7M3GonCSEwLlz52Cz2XDp0iVlu9vtxsDAAABgYGBAKfx85vF4EI/HEY1G4ff7\ncfToUdy5c0eXWZSUlKCsrAxTU1MAgOHhYdjtdrhcLt1lYbVaMT4+js+fP0MIgeHhYdhsNl1msWq9\n94Tb7Ybf70cymUQ0GkUkEsHBgwc33pnWFwRWPXr0SJjNZlFRUSE8Hk+mDpOTXrx4IQoKCkR1dbWo\nqakRNTU14vHjx2JxcVE4nU5hMpmEJEniw4cPm73UrHr27JlwuVxCCKHbLCYmJkRdXZ3Yt2+fOHHi\nhEgkErrNoru7W9hsNuFwOMSZM2dEMpnUTRYnT54Uu3fvFoWFhcJoNIr+/v4Nz72rq0tUVFQIi8Ui\nZFn+5f4z/tN4RESUXfwFJSKiPMNiJyLKMyx2IqI8w2InIsozLHYiojzDYiciyjN/A49+N1pMAiJT\nAAAAAElFTkSuQmCC\n" } ], "prompt_number": 34 }, { "cell_type": "markdown", "metadata": {}, "source": [ "In the above example, it is clear that $ts1$ and $ts2$ are most similar (they are both $sin$ functions under different transformations). $ts3$ is clearly the most different. Let's compute the Euclidean distance $d(ts1,ts2)$ and $d(ts1,ts3)$ to see if the Euclidean distance measure agrees with what our intuition tells us. Let's first create a function that computes the Euclidean distance between two time series." ] }, { "cell_type": "code", "collapsed": false, "input": [ "def euclid_dist(t1,t2):\n", " return sqrt(sum((t1-t2)**2))" ], "language": "python", "metadata": {}, "outputs": [], "prompt_number": 35 }, { "cell_type": "markdown", "metadata": {}, "source": [ "Let's now find the Euclidean distance between $ts1$ and $ts2$" ] }, { "cell_type": "code", "collapsed": false, "input": [ "print euclid_dist(ts1,ts2)" ], "language": "python", "metadata": {}, "outputs": [ { "output_type": "stream", "stream": "stdout", "text": [ "26.959216038\n" ] } ], "prompt_number": 36 }, { "cell_type": "markdown", "metadata": {}, "source": [ "and the Euclidean distance between $ts1$ and $ts3$" ] }, { "cell_type": "code", "collapsed": false, "input": [ "print euclid_dist(ts1,ts3)" ], "language": "python", "metadata": {}, "outputs": [ { "output_type": "stream", "stream": "stdout", "text": [ "23.1892491903\n" ] } ], "prompt_number": 37 }, { "cell_type": "markdown", "metadata": {}, "source": [ "This is not good because according to the Euclidean distance measure, $ts1$ is more similar to $ts3$ than to $ts2$ which contradicts our intuition. This is the problem with using the Euclidean distance measure. It often produced pessimistic similarity measures when it encounters distortion in the time axis. The way to deal with this is to use dynamic time warping." ] }, { "cell_type": "heading", "level": 2, "metadata": {}, "source": [ "Dynamic Time Warping" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Dynamic time warping finds the optimal non-linear alignment between two time series. The Euclidean distances between alignments are then much less susceptable to pessimistic similarity measurements due to distortion in the time axis. There is a price to pay for this, however, because dynamic time warping is quadratic in the length of the time series used." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Dynamic time warping works in the following way. Consider two time series $Q$ and $C$ of the same length $n$ where $$Q=q_1,q_2,...,q_n$$ and $$C=c_1,c_2,...,c_n$$ The first thing we do is construct an $n\\times n$ matrix whose $i,j^{th}$ element is the Euclidean distance between $q_i$ and $c_j$. We want to find a path through this matrix that minimizes the cumulative distance. This path then determines the optimal alignment between the two time series. It should be noted that it is possible for one point in a time series to be mapped to multiple points in the other time series." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Let's call the path $W$ where $$W=w_1,w_2,...,w_K$$ where each element of $W$ represents the distance between a point $i$ in $Q$ and a point $j$ in $C$ i.e. $w_k=(q_i-c_j)^2$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "So we want to find the path with the minimum Euclidean distance $$W^*=argmin_W(\\sqrt{\\sum_{k=1}^Kw_k})$$ The optimal path is found via dynamic programming, specifically the following recursive function. $$\\gamma(i,j)=d(q_i,c_j)+min ( \\gamma(i-1,j-1),\\gamma(i-1,j),\\gamma(i,j-1))$$ " ] }, { "cell_type": "code", "collapsed": false, "input": [ "def DTWDistance(s1, s2):\n", " DTW={}\n", " \n", " for i in range(len(s1)):\n", " DTW[(i, -1)] = float('inf')\n", " for i in range(len(s2)):\n", " DTW[(-1, i)] = float('inf')\n", " DTW[(-1, -1)] = 0\n", "\n", " for i in range(len(s1)):\n", " for j in range(len(s2)):\n", " dist= (s1[i]-s2[j])**2\n", " DTW[(i, j)] = dist + min(DTW[(i-1, j)],DTW[(i, j-1)], DTW[(i-1, j-1)])\n", "\t\t\n", " return sqrt(DTW[len(s1)-1, len(s2)-1])" ], "language": "python", "metadata": {}, "outputs": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now let's compute the Euclidean distance between $ts1$ and $ts2$ using dynamic time warping." ] }, { "cell_type": "code", "collapsed": false, "input": [ "print DTWDistance(ts1,ts2)" ], "language": "python", "metadata": {}, "outputs": [ { "output_type": "stream", "stream": "stdout", "text": [ "17.9297184686\n" ] } ], "prompt_number": 227 }, { "cell_type": "markdown", "metadata": {}, "source": [ "and now the dynamic time warping distance between $ts1$ and $ts3$" ] }, { "cell_type": "code", "collapsed": false, "input": [ "print DTWDistance(ts1,ts3)" ], "language": "python", "metadata": {}, "outputs": [ { "output_type": "stream", "stream": "stdout", "text": [ "21.5494948244\n" ] } ], "prompt_number": 228 }, { "cell_type": "markdown", "metadata": {}, "source": [ "As you can see, our results have changed from when we only used the Euclidean distance measure. Now, in agreement with our intuition, $ts2$ is shown to be more similar to $ts1$ than $ts3$ is." ] }, { "cell_type": "heading", "level": 3, "metadata": {}, "source": [ "Speeding Up Dynamic Time Warping" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Dynamic time warping has a complexity of $O(nm)$ where $n$ is the length of the first time series and $m$ is the length of the second time series. If you are performing dynamic time warping multiple times on long time series data, this can be prohibitively expensive. However, there are a couple of ways to speed things up. The first is to enforce a locality constraint. This works under the assumption that it is unlikely for $q_i$ and $c_j$ to be matched if $i$ and $j$ are too far apart. The threshold is determined by a window size $w$. This way, only mappings within this window are considered which speeds up the inner loop. The following is the modified code which includes the window size $w$." ] }, { "cell_type": "code", "collapsed": false, "input": [ "def DTWDistance(s1, s2,w):\n", " DTW={}\n", " \n", " w = max(w, abs(len(s1)-len(s2)))\n", " \n", " for i in range(-1,len(s1)):\n", " for j in range(-1,len(s2)):\n", " DTW[(i, j)] = float('inf')\n", " DTW[(-1, -1)] = 0\n", " \n", " for i in range(len(s1)):\n", " for j in range(max(0, i-w), min(len(s2), i+w)):\n", " dist= (s1[i]-s2[j])**2\n", " DTW[(i, j)] = dist + min(DTW[(i-1, j)],DTW[(i, j-1)], DTW[(i-1, j-1)])\n", "\t\t\n", " return sqrt(DTW[len(s1)-1, len(s2)-1])" ], "language": "python", "metadata": {}, "outputs": [], "prompt_number": 1 }, { "cell_type": "markdown", "metadata": {}, "source": [ "Let's test this faster version." ] }, { "cell_type": "code", "collapsed": false, "input": [ "print DTWDistance(ts1,ts2,10)" ], "language": "python", "metadata": {}, "outputs": [ { "output_type": "stream", "stream": "stdout", "text": [ "18.5965518384\n" ] } ], "prompt_number": 196 }, { "cell_type": "code", "collapsed": false, "input": [ "print DTWDistance(ts1,ts3,10)" ], "language": "python", "metadata": {}, "outputs": [ { "output_type": "stream", "stream": "stdout", "text": [ "22.4724828468\n" ] } ], "prompt_number": 197 }, { "cell_type": "markdown", "metadata": {}, "source": [ "Another way to speed things up is to use the _LB Keogh_ lower bound of dynamic time warping. It is defined as $$LBKeogh(Q,C)=\\sum_{i=1}^n (c_i-U_i)^2I(c_i > U_i)+(c_i-L_i)^2I(c_i < L_i)$$\n", "where $U_i$ and $L_i$ are upper and lower bounds for time series $Q$ which are defined as $U_i=max(q_{i-r}:q_{i+r})$ and $L_i=min(q_{i-r}:q_{i+r})$ for a reach $r$ and $I(\\cdot)$ is the indicator function. It can be implemented with the following function." ] }, { "cell_type": "code", "collapsed": false, "input": [ "def LB_Keogh(s1,s2,r):\n", " LB_sum=0\n", " for ind,i in enumerate(s1):\n", " \n", " lower_bound=min(s2[(ind-r if ind-r>=0 else 0):(ind+r)])\n", " upper_bound=max(s2[(ind-r if ind-r>=0 else 0):(ind+r)])\n", " \n", " if i>upper_bound:\n", " LB_sum=LB_sum+(i-upper_bound)**2\n", " elif i