{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Signed Cumulative Distribution Transform Nearest Subspace (SCDT-NS) Classifier\n",
"\n",
"This tutorial will demonstrate how to use the SCDT-NS classifier from the *PyTransKit* package to classify 1D signals."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Class:: SCDT_NS\n",
"**Functions**:\n",
"\n",
"1. Constructor function:\n",
" scdt_ns_obj = SCDT_NS(num_classes, rm_edge)\n",
" \n",
" Inputs:\n",
" ----------------\n",
" num_classes : integer value\n",
" totale number of classes in the dataset.\n",
" rm_edge : [optional] boolean \n",
" IF TRUE the first and last points of CDTs will be removed.\n",
" \n",
" Outputs:\n",
" ----------------\n",
" scdt_ns_obj : class object\n",
" Instance of the class SCDT_NS.\n",
" \n",
"2. Fit function:\n",
" scdt_ns_obj.fit(Xtrain, Ytrain, Ttrain)\n",
" \n",
" Inputs:\n",
" ----------------\n",
" Xtrain : array-like, shape (n_samples, n_columns)\n",
" 1D data for training.\n",
" Ytrain : ndarray of shape (n_samples,)\n",
" Labels of the training samples.\n",
" Ttrain : [optional] array-like, shape (n_samples, n_columns)\n",
" domain for corresponding training signals.\n",
" \n",
"3. Predict function:\n",
" preds = scdt_ns_obj.predict(Xtest, Ttest, use_gpu)\n",
" \n",
" Inputs:\n",
" ----------------\n",
" Xtest : array-like, shape (n_samples, n_columns)\n",
" 1D data for testing.\n",
" Ttest : [optional] array-like, shape (n_samples, n_columns)\n",
" domain for corresponding test signals.\n",
" use_gpu: [optional] boolean flag; IF TRUE, use gpu for calculations\n",
" default = False.\n",
" \n",
" Outputs:\n",
" ----------------\n",
" preds : 1d array, shape (n_samples,)\n",
" Predicted labels for test samples.\n",
" "
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Example\n",
"The following example will demonstrate how to:\n",
"* create and initialize an instance of the class SCDT_NS\n",
"* train the model using 1D signals\n",
"* apply the model to predict calss labels of the test 1D samples\n",
"In this example we have used a synthetic dataset (1D). The dataset contains three classes.
\n",
"Class 0: 4-th degree polynomial applied to a Gabor wave
\n",
"Class 1: 4-th degree polynomial applied to a sawtooth wave
\n",
"Class 2: 4-th degree polynomial applied to a square wave"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Import python libraries"
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {},
"outputs": [],
"source": [
"import numpy as np\n",
"from numpy import interp\n",
"import math\n",
"import matplotlib.pyplot as plt\n",
"from scipy.linalg import lstsq\n",
"from scipy import signal\n",
"\n",
"import sys\n",
"sys.path.append('../')\n",
"from pytranskit.classification.utils import *"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Import SCDT-NS class from PyTransKit package"
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {},
"outputs": [],
"source": [
"from pytranskit.classification.scdt_ns import SCDT_NS"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Define a three template signals: Gabor wave, sawtooth wave, square wave"
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {},
"outputs": [],
"source": [
"def signal_gabor(t):\n",
" w = 30\n",
" s1 = np.real(math.pi**-0.25 * np.exp(2*math.pi*1j*w*(t-0.5)) * np.exp(-200*((t-0.5)**2)))\n",
" return s1\n",
"\n",
"def signal_sawtooth(t):\n",
" s2 = signal.sawtooth(2 * np.pi * 15 * (t-0.5))* np.exp(-250*((t-0.5)**2))\n",
" return s2\n",
"\n",
"def signal_square(t):\n",
" s3 = signal.square(2 * np.pi * 15 * (t-0.5))*np.exp(-250*((t-0.5)**2))\n",
" return s3"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Generate template signals for three classes"
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAABYQAAADSCAYAAAD6xMQPAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+j8jraAAAgAElEQVR4nOzdeZyVZf3/8fc1G9sAoggIyKaComlqKVquaWqampm7CJk/s2zxW37FvmpmllpalGYqXwTUr2BSGZqpZQZmCi4oLsiI7CAw7MPMMMuZ6/fHdW7mzJmz3PeZs5/X8/GYx32W+9z3dW4drrnf53M+t7HWCgAAAAAAAABQ/MpyPQAAAAAAAAAAQHYQCAMAAAAAAABAiSAQBgAAAAAAAIASQSAMAAAAAAAAACWCQBgAAAAAAAAASgSBMAAAAAAAAACUCAJhIIoxZoQxxhpjKnI9FgAAkH3hvwP2z/U4AAAAgEwgEEZRMsZcZIyZb4ypN8ZsDN/+ljHG5HpsAAAUM2PM540x/zHGbDfGbDHGvGKM+Wya95G2wNYY8y9jzDfSsS0AAIpBNuZyALlFIIyiY4z5gaTfSPqlpEGSBkr6pqTPSarK8lioMgYAlAxjTB9Jz0i6V9KekoZI+omkplyOCwAA+JNPc7kxpjzb+wRKBYEwiooxpq+k2yR9y1o721pbZ52F1tpLrbVN4fXONMYsNMbsMMasNsbcGmNzXzfGrDPGfBIOmb19dDPGTA4/ty58u1v4uRONMWuMMTcYY9ZLmhZjjCuNMUeGb18WrnIaG77/DWPMU+HbRxljXjXGbAuP4T5jTFX4uQeMMXdHbfcvxpj/Ct8ebIz5ozGm1hiz3Bjz3S4fXAAAkhstSdbamdbakLW20Vr7grV2kSQZY/YzxvzTGLPZGLPJGPN/xpg9ws9NNMY87W3IGLPUGPOHiPurjTGfNsbMCz/0jjFmpzHmwvDzV4Vfs8UYM8cYMzjitccaY14PVzq9bow5Nvz4zyQdJ+m+8Lbui3gvpxhjPjLGbDXG/C7Wt4yMMd2NMY3GmP7h+zcZY1rDJ9MyxtxujJkcvh33bw9jzHPGmGujtv2OMea88O0DjTF/D7+3JcaYCwL+dwEAwK9kc3m5Mebu8Dy+zBjzbRPRctEYs8IYc4q3MWPMrcaYxyLuP2mMWR+ek+cZYw6OeG66Meb3xphnjTH1kk7ye25rjBkZPncuC9//X2PMxojnHzPGfD98e6IxZrExpi78Hq6OWG+xMeasiPsV4fd6RPj+OOOqp7eF5+oTu3a4gdwgEEaxOUZSN0l/SbJevaTxkvaQdKaka4wx50atc5KkAyR9UdKkiEntfySNk/RpSYdJOkrSTRGvGyT3SepwSf8vxr7nSjoxfPt4ScsknRBxf274dkjSdZL6h9/XFyR9K/zc45Iu9E5OjTH9wuOcFZ4An5b0jtynuV+Q9H1jzGlJjgkAAF1VIylkjJlhjDkjPD9FMpLukDRY0kGS9pV0a/i5uZKOM8aUGWP2kVQp9+0eGWNGSaqWtMhae3x4/cOstdXW2ieMMSeHt3uBpH0krZQ0K/zaPSX9VdJvJe0l6VeS/mqM2cta+z+SXpZ0bXhbkaHsWZI+KzfXXyCp0zxqrd0l6XV1nMdXeuNWx3k90d8ej0u6ePdBch8UDw+Ps5ekv4fXGRBe7/7IE2gAANIo2Vx+ldwcebikz0g6P+D2/yZ3nj1A0luS/i/q+Usk/UxSb0n/kc9zW2vtckk7wuOS3Ae+O40xB4XvR87JG8PvoY+kiZJ+7QW+kmYqYk6Wm/83WWvfMsYMkfub4na5c/4fSvqjMWbvgMcAyDkCYRSb/nL/WLd6D0R8etdojDlekqy1/7LWvmutbQt/0jlT7Sdznp9Ya+utte/KVfp6k8Klkm6z1m601tbKfX3m8ojXtUn6sbW2yVrbGGOMcyP2dZzcCax3/4Tw87LWvmmtfc1a22qtXSHpwYj1XpZkw6+X3CT8qrV2ndzJ697W2tustc3W2mWSpki6KOnRAwCgC6y1OyR9Xm6OmiKpNlytOzD8/FJr7d/Dc2StXDh7Qvi5ZZLq5D5wPUHS85LWGmMODN9/2VrbFmfXl0p62Fr7VvjbQDdKOsYYM0IufP3IWvtoeE6dKelDSV9O8nbutNZus9aukvRSeFyxzJV0Qrgy6lC54PkEY0x3uTn55fD7S/S3x58lfdoYMzzi/fwp/F7OkrTCWjstPP63JP1RwU/AAQBIKtlcLvch6WRr7Wpr7Ra589kg2384/E3eJrkPhQ8z7pu+nr9Ya18Jz/mfUrBzW29OHhS+Pzt8f6Rc+PtOeAx/tdZ+HP428VxJL6j93PpxSWcbY3qG718SfkySLpP0rLX22fB8/ndJb0j6UpBjAOQDAmEUm82S+puI3r3W2mOttXuEn/O+PnK0Meal8NdOtsv1GO4fta3VEbdXylUzKbxcGec5SaoNVwzF41VADZJULukJSZ8Ln7T2lfR2eIyjjTHPhL9Os0PSz70xWmutXOWTF1JfovZPVodLGhwOwbcZY7ZJ+pFcL2UAADLKWrvYWjvBWjtU0iFyc6TXNmGAMWaWMWZteG57TB3nX+9bNF4Vz7/kQtPdH5jG0WFuttbulJv3h0Q/F7Yy/Fwi6yNuN8hVKMfijfkISe/KVfOeIPdtoqXW2k1S4r89rLV1chVH3gnuReo4rx8dNa9fKveNJAAA0i7RXB6+HX2u7Eu43cSdxpiPw38HrAg/Ffm3QOS2g57bRv4dMU8d/47Y/cFyuPL5NeNaMW2TC3S9OXmppMWSvhwOhc9WeyA8XNLXosbzeblvJwEFhUAYxeZVuWb35yRZ73FJcyTta63tK+kBua+xRto34vYwSevCt9fJTQSxnpPcJ6lxhSeYBknflTQvfBK4Xq69xL8jqp9+L1fBdIC1to/cxBc5xpmSzg9XEx0tVy0kuQl0ubV2j4if3tZaPrUEAGSVtfZDSdPlTiYlV0VkJR0antsuU8e5zTuROy582/tWTbJAuMPcHG6zsJektdHPhQ0LPyclmbd9+I+kMZK+ImmutfaD8PbPjBpzsr89Zkq62BhzjKQeclXJkpvX50bN69XW2mu6OG4AAJKKMZd/os7nypHqJfWMuB/5AeYlcufqp8gVQ40IPx45H0bOy0HPbefK/Q1xYvj2v+XaOO3+O8K46//8UdLdkgaGi8eeVYw5OTzWD8Ln8N54Ho0aTy9r7Z1xxgPkLQJhFBVr7Ta5Fg73G2PON8ZUh3sRflpSr4hVe0vaYq3dZYw5Sm5iinazMaZnuEffRLlKXslNDjcZY/Y27iIyt8hVOAUxV9K1aj9R/FfUfW+MO+T6Hh0oqcOJn7V2oaRaSf8r6fnwe5ekBZJ2GHdhux7hT2EPMcZ8NuAYAQAIxLiLn/3AGDM0fH9fuROq18Kr9Ja0U9K2cB++66M2MVeuh38Pa+0auXYLp8uFuwsj1tsgaVTE/cclTTTuonPd5L5VMz/cculZSaONMZeELwxzoaSxcldQj7WtQKy1DZLelPRttc/j/5F0tTrP64n+9nhWLri+TdITER8QPxMe/+XGmMrwz2cjeiICAJA2PubyP0j6rjFmaLi/8KSoTbwt6aLwfBXdY7i3XAHXZrnQ+OdJhhPo3NZa+5GkRrkPnOeF219skPRVtc/JVXLXHaqV1GqMOUPuejyRZoUfu0bt1cGSO+//sjHmtPBYuht3YfmhSd4HkHcIhFF0rLW/kPRfkv5brln8Brn+uzfInaBJ7uJstxlj6uQC3T/E2NRcSUslvSjpbmvtC+HHb5frE7RI7quhb4UfC2Ku3GQ4L859yTWov0Sun+IUtQfSkWbKfbq6e5Ky1obk+iJ+WtJySZvkQuO+MV4PAEA61cl9a2W+cVcHf03Se5J+EH7+J3KtFbbLtUj4U+SLrbU1coGx13d3h9zFV18Jz2+eWyXNCH9d8wJr7YuSbpar+PlE0n4Kt1+w1m6W68P7A7kT0P+WdJbXykHSb+S+cbPVGPPbFN/3XLmL4C2IuB89ryf82yPcS/FP6jyv18mdlF4kV+28XtJdciezAACkW7K5fIpcn/935M6F/xT1+pvl5uGtcvN+ZKD6iFyLibWSPlB7yBxTiue2cyVtDl8DwLtvFP5gOTyvflduHt4qd849J2q/n8h9+/hYRZyHW2tXy1UN/0guUF4t9+E22RoKjnGtSAEAAAAAAAD/wtfCWS6p0kZc3B1AfuNTDAAAAAAAAAAoEQTCAAAAAAAAAFAiaBkBAAAAAAAAACWCCmEAAAAAAAAAKBEEwgAAAAAAAABQIiqCrNy/f387YsSIDA0FAFDq3nzzzU3W2r1zPY5Cx3wNAMgk5uv0YL4GAGRavDk7UCA8YsQIvfHGG+kbFQAAEYwxK3M9hmLAfA0AyCTm6/RgvgYAZFq8OZuWEQAAAAAAAABQIgiEAQAAAAAAAKBEEAgDAAAAAAAAQIkgEAYAAAAAAACAEkEgDGTAokXSN74hff/7UmNjrkcDAABiam2VfvpTadu2XI8EAAAksXDhQt18881asGBBrocCFDwCYSADHn1UmjpV+s1vpDffzPVoAABATK++Kt1yi/SPf+R6JAAAIImHH35Yt99+u44++mhdeuml2rp1a66HBBQsAmEgA3bubL9dX5+7cQAAgARqatwyFMrtOAAAQFKtra3q16+fbr31Vj355JM6+eSTtWnTplwPCyhIBMJABkSGwATCAADkKS8QbmvL7TgAAEBSbW1tqqqq0o9//GM9/fTT+vDDD3XmmWeqqakp10MDCg6BMJAB9fVS9+7ttwEAQB4iEAYAoGC0tbWprMzFWKeddpr+7//+TwsWLND3vve9HI8MKDwEwkAG1NdLe+/dfhsAAOQhAmEAAApGW1ubysvLd98/77zzdP311+vBBx/Uc889l8ORAYWHQBjIgPp6acCA9tsAACDPhELS0qXuNoEwAAB5L7JC2HPbbbdpzJgxuvrqq1XPyTfgG4EwkAFUCAMAkOdWrZKam91tAmEAAPJerEC4e/fumjJlilatWqV77rknRyMDCg+BMJABDQ1Snz5St24EwgAA5CWvXYREIAwAQAGIFQhL0nHHHafzzz9fv/jFL/TJJ5/kYGRA4SEQBjKgvl7q1cv9EAgDAJCHIgPhUCh34wAAAL6EQqGYgbAk3XHHHWpqatJdd92V5VEBhYlAGMgAAmEAAPIcFcIAABSUeBXCkrT//vvrkksu0ZQpU7Rp06YsjwwoPATCQAYQCAMAkOdqaqR993W3CYQBAMh7bW1tKi8vj/v8DTfcoIaGBt17771ZHBVQmAiEgTRrbXXXqPEC4YaGXI8IAAB0UlMjHXigu00gDABA3ktUISxJY8eO1TnnnKN7771XO3fuzOLIgMJDIAz4tG2btH178vW8iuCgFcJr1tDCEACArNi1S1q5UjroIHefQBgAgLyXLBCWpEmTJmnr1q166KGHsjQqoDARCAM+9esnHXpo8vW8ALhnT/fjJxBuaHDfWp0woUtDBAAAfnz8sWQtgTAAAAXETyA8btw4nXjiifrVr36llpaWLI0MKDwEwoAP3nVnVq1Kvm4qFcKbN7vlY4+581MAAJBB3sROywgAAApGKBRKGghL0nXXXae1a9dqzpw5WRgVUJgIhAEfHn/cLbt3T75uKoFwZCuKt98OPj4AABDAkiVuOWaMW9KzCQCAvOenQliSzjzzTA0fPly/+93vsjAqoDARCAM+eIFwr17J1/UuIpdqIOztCwAAZEhNjbTPPlLfvu4+FcIAAOQ9v4FweXm5rr76ar300ktau3ZtFkYGFB4CYSCJpibpo4/c7e3bk7d06GqF8LvvpjZOAADgU02NNHq05J1UEggDAJD32traVF5e7mvdUaNGSZJ27NiRySEBBYtAGEiirs4tBw+WWlvdhckTiQ6EGxuTn2d6c9Tgwe37AwAAGUIgDABAwfFbISxpd3Acoi0UEBOBMJCEF9AOGeKWkdW8sUQHwlJ7G4l4vG0OHSrt3JnaOAEAgA9bt0q1tQTCAADkyBtvvKERI0Zo/fr1gV5HIAykD4EwkIQX0A4d6papBMLJ2kZEBsJUCAMAkEFeH6jRoyXva6cEwgAAZM19992nlStXatWqVYFeRyAMpA+BMJCEF9AGDYR79nQ/kY/Fs327OycdOJBAGACAjKqpccvRoyVj3G1OFgEAyIqdO3dq9uzZkoKHtaFQiEAYSBMCYSCJdLSM8BMI9+0r9e5NIAwAQEbV1LhWEeGLzaisjAphAACyZPbs2aoPnyAHDWuDVAh76xEIA7ERCANJpFIhXFEhVVUF6yHct69UXS01NUktLV0bMwAAiKOmRho50k3UEoEwAABZNG3aNJnwN3RSCYS9yt9kqBAGEiMQBpJIpYewFwSnUiEcuU8AAJBmNTWuXYSHQBgAgKz4+OOPNW/ePJ166qmSXMAbRCo9hIPuAygVBMJAEtEVwjt2JF6/oSF4ILxjh9SnT3sgTNsIAAAywFoCYQAAcuSRRx6RMUYTJ06UlNmWEVQIA4kRCANJeOHsPvu4JRXCAAAUqE8+cZNyZCBcXk4gDABAhrW1tWnGjBk69dRTNXz4cEkEwkAuEQgDSdTVuZ7APXq4wDaTgXB1dfs+AQBAmtXUuCUVwgAAZNXcuXO1cuVKTZgwIeWwNhQKEQgDaUIgDCSxc6cLgo1xoW02KoQJhAEAyIB4gTAniwAAZNT06dPVp08fnXvuuSmHtUEqhL31CISB2AiEgSTq6tqDWj+BcEOD1LOnu+0tEwXC1hIIAwCQFTU1Uvfu7RcGkKgQBgAgw+rq6jR79mxddNFF6tGjR8phLS0jgPQhEAaSqKtrb+XQp0/yQHjXLqlbN3fbWzY1xV+/ocEVJkW2jKCHMAAAGVBTIx1wgAuBPQTCAABk1OzZs9XQ0KAJEyZISj2sbWtr2/3aZLz12pjjgZgIhIEkoiuEd+xIvH5zc3sQXFbm+g83N8df39seFcIAAGTYkiUd20VIBMIAAGTYtGnTNHr0aI0bN05S6mEtFcJA+hAIA0l4PYQlfy0jmpraA2HJ3U5UIextr08fAmEAADKmpUVatqxzIFxeTiAMAECGfPzxx3r55Zc1YcIEGWMkda1CmEAYSA8CYSCJoD2EUw2E+/aVevRwhUq0jAAAIM1WrJBaW6kQBgAgi6ZNm6aysjJdfvnlux8jEAZyj0AYSCKyh3CmA2Fj3L6oEAYAIM1qatySQBgAgKwIhUKaPn26TjvtNA2NuKBrqmFtKBQiEAbShEAYSCKyQri62l00LtF5Y9BAuL7eLXv1csvevQmEAQBIu0SBMCeLAACk3QsvvKC1a9fqyiuv7PC4F+pmskI41X0ApYJAGEgisodw9+5umSjgDRoIe8/16OGWvXvTMgIAgLSrqZH69ZP22qvj41QIAwCQEVOnTlX//v315S9/ucPjtIwAco9AGEigudn9eIGwF/Tu2hV7fWuDB8LetrzXUCEMAEAG1NS46uDwBW12IxAGACDtamtrNWfOHF1++eWqqqrq8JwX1rYFnH/b2tp2vzaZVPcBlAoCYSABL5j1egh7FcLxAuGWFrdMJRD2tk0PYQAAMsALhKMRCAMAkHaPPfaYWlpaOrWLkKgQBvIBgTCQgBfMRreMiBcIe8FvVwJhKoQBAEiz+nppzZrYgXB5OYEwAABpZK3V1KlTdfTRR+vggw/u9DyBMJB7BMJAAl4vX789hL3HI78RU1Xlr4dwZCBMD2EAANJo6VK3pEIYAICMW7Bggd5//319/etfj/l8qmFtKBQiEAbShEAYSCC6QjhZD+GuVAh7ITIVwgAApFlNjVsSCAMAkHEPP/ywevbsqYsuuijm89moEPbWIxAGYiMQRknatk365jelLVsSr1df75a9erllqi0jmpvj72PXLhcGe/Nar17+KoT/9S/pZz9Lvh4AACXvww/d8oADOj9XViZxsggAQFrs3LlTM2fO1Ne+9jX16dMn5jqphrW0jADSh0AYJemRR6QHH5T+8Y/E6zU0uGWPHm6ZLBD2gt+gFcLedr19NTZK1iYe25Qp0k03ScuXJ14PAICSt3ixNHx4+ye8kagQBgAgbR599FHV1dXp6quvjrtOVyqEvdcm463XxhwPxEQgjJI0c6ZbrliReL3GRreMDoST9RAOEgg3NXUOhK1NXFUstY/9iScSrwcAQMlbvFg68MDYzxEIAwCQFtZa3XfffTryyCM1bty4uOulGtZSIQykD4EwSs7y5dJrr7nbQQPhTPUQjlzf25e373i8sc+alXg9AABKWlubtGSJdNBBsZ8vLycQBgAgDV566SV98MEHuvbaa2WMibteNnoIEwgDiREIo+TMmeOW++yTvN1CvArhdAfC0RXCkfuO95p169x7eOcd2kYAABDXqlVuUo0XCFMhDABAWtx7773aa6+94l5MzkMPYSD3CIRRclaskKqrpWOPTb1lRK4D4VWr3PL0091y5cr46wIAUNIWL3bLXATCra3SwoXSf/5D6AwAKGorV67UnDlzdNVVV6l75AluDMYYlZWVBQ5rQ6EQgTCQJgTCKDnr10uDBkkjRrhAONHF2+K1jEh3D+GgLSO8imCvLdP69fHXBQCgpPkJhDNxsvjnP0tjxkhHHCF97nMuFAYAoEg9/vjjamtr0zXXXONr/fLy8oxWCHstKwiEgdgIhFFyNmxwgfDIka46d8OG+Os2NrrzxKoqdz/VCuFQKP65ZnSFcM+e7fuOx6ts9gLhRO8BAICStnix1L+/+4klExXCv/+99NWvSr17S7fc4h7bsiW9+wAAII9s2rRJvXr10rBhw3ytn0qFcJBAWEotdAZKBYEwSs769dLAga5CWErcf7ex0VXsev3wUw2EI5+LlkrLiOXLpcpK6eCD3ZIKYQBAUVq9uuvVux9+GL86WEp/IPzSS9K3vy2deab06qvS+ee7x5ub07cPAADyTHNzs6q8Siofgoa11lpZa3e3gvC7j7aAc/yaNWu0cOHCQK8BChGBMEqO1zJi5Eh3P1EfYS8Q9njhbq4D4RUrpOHD3YXRBw6kQhgAUITeflsaNkx66KGubWfx4uwFwtu3S+PHSwccIM2a5Sb1ZH8IAABQBJqamtQt8kQ4iaBhrQ33esx0hfBPfvITff7zn9fGjRsDvQ4oNATCKClNTdLWrcErhD3l5VJFRfIewpEfjHq3E70mVg/hhob441qxon38AwdSIQwAKEJ//KNbJpqok6mtlTZvlg48MP465eXpC4Tvvltas0Z65BGpVy/3mPeHABXCAIAilukKYS88znQgXFdXp4aGBt11112BXgcUGgJhlBTvQ75Bg1yv3upqadOm+OtHB8KSq+ZNpUI43nlgKhXCtbUuCJbce6FCGABQdObNc8t99019G8kuKCelr0K4tlaaPFn62teko49uf5wKYQBACUilQjhIWOutm+lAuCk8X99///365JNPAr0WKCQEwigpXnDqhal9+7pvd8aTzkA4nS0jtm93Y5eoEAYAFKHWVumVV9ztrgSp2QyE771Xqq+XfvKTjo8TCAMASkBzc3NGA+FsVQg3NTVp0KBBamlp0R133BHotUAhIRBGSfGC00GD3DJZINzQkH+BsLUdA+FBg1zlc7ovkA4AQM68+Wb7xeS60mrhww/dV4ISVRmXlXX9wnUtLdL//q90xhmdw2daRgAASkBTU1OglhFlZWUZD4SD7kNywfaoUaM0ceJEPfjgg1q9enWg1wOFgkAYJSVoIByrQrhbt/jhrneuV1nZcX0peA/heIFwY6MrnIoMhFtbpS1b4r8PAAAKyrPPSsa4212tED7wQBf6xpOOCuFnnpE++US6+urOz1EhDAAoAZluGeEFwuXl5Rnbh9QebN90002SpEmTJgV6PVAoCIRRUryWEQMGuGUmWkZ069Z+Diulv0LYG29kywiJthEAgCLy1FPScce56tquVNYuXpy4XYSUnkB46lRpyBDpS1/q/FxFhVtSIQwAKGKpXFSuLcD8m2rLiCD7kNpbXwwfPlw33HCDHn/8cc2dOzfQNoBCQCCMkrJ+vbTHHu0BbKYC4UiJAuFQyFX3RgbCFRXux28g7FU7c2E5AEBRWLZMWrRIOvfcxF/LSWbbNmnVKulTn0q8Xnl51wLhbdukF16QLr64PfyNZEzX3gcAAAUgWxXC2egh7L2PSZMmafjw4fr2t7+tlpaWQNsB8h2BMErK+vXtFbWS1KdPeltGBA2EY/Ucltw+qRAGAJSkv/zFLc85x02QqVbWLlrklocemni9rlYIP/206yF8/vnx1+lqpTMAAHkulQrhfAyEI99Hz549NXnyZL3//vu67777Am0HyHcEwigpGze2t4uQcl8h7G0nskJYChYIe+9n48bY6wMAUFBmzpQOO0waNcoFqalW1r7zjlsedlji9boaCM+e7S5ad9RR8dehQhgAUOQyXSHsrZvNCmFJOuecc3TGGWfopptu0pIlSwJtC8hnBMIoKZs3S3vt1X6/b18XysYr2mlsdBcnj5SNQLhnT6mhIfY+duxwyz592t9DWZl7bwAAFLTFi6XXX5fGj3f3uxKkLlok9e8v7bNP4vW6Egg3NEjPPy995SsdLyAQjUAYAFDkvN67fmWjQrisrCzli8p5jDGaMmWKunfvrosvvlhNzOcoEgTCKCmxAmGpPWSNVggVwmVl0p57EggDAIrAjBmup++ll7r7XWm18M47rl1EoqBWchNpwJPF3ebNcxN8rIvJRaJlBACgyEUHqckEDWuz2TIiOtgeMmSIpk2bpoULF+rGG28MtD0gXxEIo2RYGz8QjtU2IhRyLQHT1UM41nlgOnoIS+49EQgDAApac7P0yCPS6ae3N8hPtbI2FJLeey95uwipaxXCzz/vxnj88YnXo0IYAFDksnVRufLy8kD7aAs4x8d7H2effbauvfZa/frXv9Yf//jHQNsE8hGBMAqetdI110gTJiReb+dOF/D6DYS9QDZohXD0h6Le/XRWCBsj9e7d/pifQLi2VjrySOnJJxOvBwBATvzhD9Inn0jXXtv+WKqVtUuXuonUTyBcXp56IPzCCy4Mjv5jIRoVwgCAIpfKReWChLW5uOvTjVsAACAASURBVKhctF/+8pcaN26cLr30Us2dOzfQdoF8QyCMgjdjhvTAA2759tvx1/MC08hA2OvDm85AOBstI3r3dgVNHj+B8P33S2+9JV15pbRiReJ1AQDIKmulX/9aOugg6bTT2h9PtbLWu6DcoYcmXzfVCuE1a6QPPug43nioEAYAFLlsVQhnMhC21iZ8H927d9czzzyjkSNH6pxzztGiRYt8bxvINwTCKGgtLdJ110nHHCNVV0v33BN/3ViBcCYqhLMRCEe2i5CSB8KNjdLvfieNG+fu33RT/HUBAMi6Z591n1r+13917PnbrVtqlbXvvCNVVEhjxyZfN9VA2KsMOvnk5OtSIQwAKHL5eFG5oPsIhUKy1iasdN5rr730/PPPq7q6WqeeeqoWLlzoe/tAPiEQRkGbP1/atk364Q+liROlWbNca4hYggbCDQ1uma4ewrFek2oP4aCB8DPPuJYRt9/uLoT+/POpfzsWAIC0CoWkSZOk/feXrrii43NVValXCB94YOcJNpZUA+GXX3ZfNfJThUyFMACgiIVCIYVCocAtI4KGtVJmA+Gm8FydLNgeNmyYXnzxRXXr1k0nn3yydsYLIYA8RiCMgvbCC+487uST3U9rq/v2ZizprBBubo597tjc3Pncs6zMFSllukK4sTH+a95+243huOOkU06RNm1K3F4DAICseeghdwG4n/9cqqzs+FyqQeqiRf6CWslN1AH7C0qS5s2TPvc514M4GQJhAEABaGtr0913362RI0dq2bJlvl/XHP4WTJAK4bKysoxXCAfdhxcI+wm2x4wZo3vuuUfbtm3Txx9/7HsfQL4gEEZB+/vfpc9+VtpjD+ngg91j770Xe910BsJS/IrfWHNgvPPAdAbCUvwq4ffekw44wBVanXKKe+zvf4+9LgAAWbN8uXT99W5yOv/8zs+n0mphyxZp9Wp/F5STUqsQ3rRJWrzYfdLqBy0jAAB5bt68eTrqqKN0/fXXa8WKFXrjjTd8vzZIkOpJtWVEuZ8PYlPcR9Bge+jQoZKkdevW+d4HkC8IhFGwtm2TFiyQTj3V3R81ygWr778fe30vLN1zz/bHKitd+BokEPbmhlh9hFMNhGO1jPBaVkRLJRB+/33pkEPc7X32cbcJhAEAOdXQIF1wgQtkp07t2DvYk0plrXeBl0wGwv/+t1v6DYSpEAYA5CFrrf71r3/p3HPP1QknnKANGzbo/vvvlxQs5PTbaiFSeXn57pDXj1R7CAfZR9D3MWTIEEnS2rVrfe8DyBcEwihYL7/szt++8AV3v7zcXZw8USDcp0/nb6P27Zu7CmHvsXgVwtZ2fs327e59REoUCDc0SMuWtVdQS64Q65VXODcFAORIU5N0ySXSm29Kjz0mDRsWe71UKmu9nkh+W0aUlwcPhF991f1B8ZnP+FufCmEAQJ6w1uqtt97ST3/6Ux188ME66aSTNG/ePN12221asmSJvvnNb6pbt26BQk6vsjYbFcKZ7CEc9H0MGjRIEoEwClNFrgcARNu1y10s7tOf7lwJG2nuXHd+NW5c+2MHHyy99FLs9Tdv7tguwtOnj7RjR+fHkwXC6awQjhUIt7VJLS3uPUbasSNYhfDixS5YjgyETzhBmjxZev116fOf7/waz8cfu0rsI46IXbgFAChRW7ZIK1dKn/qUa1IfxMaNLgx+8UXp3nuls8+Ov24qlbULFkhDh7qvxPiRSoXwggXuD5XoCTyeVCuEd+2S/vlPF0CvWiVt2OA+Ga6ocIH0Hnu4PwL695dGjpQmTOj8hwMAoORt3LhRDzzwgBYsWKDXX39dGzdulCQdc8wxmjZtmi688EL1iDjxHTJkSFYqhPMtEA76PqqqqjRgwIDALSMaGxv1/vvva5999tHAgQNVEfRvKSAN+L8OeWXaNOmHP3TnmT16SL/6lfTNb8Zed9486eijO56LHXKIKzTats2dI0WKFwgnqxDu2bPj49kIhL19NjZ2PK9ranI/QQJhr2LaaxkhtYfA8+bFDoRbW925+pNPuvuHHio98oj/b98CKALNzdKdd3Z8LNbXFtL5WKa3zz5Te6y11U2sW7a4iWbdOtdDV5L+3/+THnyw82tiaWpyE/0tt7iJd/p06YorEr8m1UD4qKP8rx80EA6FpDfeSD72SEErhDdskO66y7XS2LHDVTEPGSINHOj+wAmF3PY++kh67TX336OlRVqyRLrnHvcHxJo17f/NtmxxXxlqaXH/PSOXQcPwRNL16TGfQmfOKadIxxyT61EgzWbNmqWampoOj5kYv0epPlYor8uHMWT7dZWVlerRo4d69uy5+6dPnz7ae++9teeee6q8vFwtLS0655xzNH/+fI0dO1ann366TjrpJJ1xxhkaOHBgp31ILhBOpUI4k4Gwt24+VQhLwY9VY2Ojxo0bp0XhFlfGGA0cOFD9+vVTdXW1evfurerqalVXV6tbt26qqKhQRUWFysvLd9+OfizW/yvJ2Fh/82XhtfBn0KBBuuqqqzK6DwLhFKXr3CuT53W5HmMo5M4/NmxwwWffvq6lQ3V159c2N0vXXSfdf7+rXv3Wt6QpU6Rrr5XGjpWOP77j+nV10ltvSTfe2PFxrwr2gw+kY4/t+FyqgbDfHsLWJg6EY50HJuoh7I0hMvz1KpmDtIx4/313Hrr//u2P9e/vAuJ586Qf/ajza266yYXBkya53sy33uqO57Rprt1jNGvdf+elS935Zo8ebox9+ki9e7slBUtAgWluln784/RuM9Yfq9GP+Vkn3Y+xz8SPlZVJ/fq5JvwjR7qv5owe7Sbbhx5yn2rut59bx1uvvNxNinV17usm8+dLzz7rguVjj3UhcuQnlfEEDVI3b3b7C/IHdNBA+MMPpZ073afSfgUJtp94wn0aXlcnXXSRdPnl7tPbXr3iv8Za90fTr34lzZ7tLqrHyRpi6dmTQLgIPf7443r66adzPQzkGWOM9tprL/Xu3VvLly/XrFmzdOGFF/p67eDBg/Xmm2/63lc2LyoXJBAuKyvLaIWw5I6VnwrhTZs26ZVXXtGMGTO0aNEiTZ48Wd27d9e6deu0bt06bdu2TTt37tTOnTu1cuVK1dXVqaWlRa2trR1+QqHQ7ttB+iOjsBxxxBHFFwhPnSp95zsdH8u3IJO/nzPHGPc36DnnuIvBDR0qLVwo3XyzK+i5/nrp5z9334Q8/XTXmm/8eFfwEvlv8n/+4wLn6KB4zBi3XLo0diA8enTnMVVXS+vXd348aA/hlha3DNpDuLLSnYtGigyEI9XVuWXv3p2336tXe8FWpKVLXagb/S2U4493Vb+trR2fmz/fFSVdfbV0xx3usbPOchd/v/BC6fnnpe9/3xUpvfuu9Ne/SnPmuPPvRLp1c+OuqnLHqqXFneNHFydF5hDpuo3MW7bM/T+BItKrl/uHNprfUBHFr6XFVQv/9rfJ1x00yLWGGD9eOvlk///PeBOotf5e8/rrbhm0QjjAyaIWLAi+Dz+BsLXSz37m/igaN859Cnvggf62b4yrDN6xw02qo0e78L5/fxfQ77mn+4PHazURuQxwtfak48+n7SA2/r0uSk899VSH+7Gq91J9rFBelw9jyMV7bm5uVmNjoxoaGtTY2Kj6+npt375dtbW1qq2t1caNG7V69Wp985vf9B0GS67q9emnn5a11lf1aSoVwkHD2nxsGSG5Y/W69/dHBGutFi5cqJkzZ+pvf/ub3g9/dbesrEy33nqrvve97/neRzxtbW2B3l+0VCqL0/Fa5IesB8IHH+wKGKKlUsSSyWKYfFsn1/tPZYzGuMrVAQNcwLl5s7t2zF//Kt1wg/vxDB7sCmIiq0/79HEVw6eeKv3+9y6E9Pztby5YjC5wGDbM7Xf58s5ji1ch3Lu3K/SJFrSHsHeeF7RlRKz2g/ECYW+c0YGw5N5brArh5cvdOWG0k092x/ff/5ZOPNE9Zq377zJggHT33e3r7rOP6838ox+5do8PP9z+XFWVu7Dftde2V4A3Nrpz0ro6t4y83dzsXlNZ2f7jnYtG/o2TrtvIjujfExQBYwgOkFhlpZuQW1vdV222bnU/W7a4gNX7tHL4cBcIp/L/U1WV+0c9FPLXq3jBArcfvxd7k4JXCC9Y4L6+c8AB/l/jVTonCrbvvtuFwZdf7ioooq+Cm0z37tKjjwZ7TTrx7wWQM0ECMsCPwYMHq6GhQdu3b9ce0b0YY8hmhXB5gA8yy8vLA1XRptIyYvDgwdq4caOam5tVVVWl5uZmPfroo7r77rv14YcfqqKiQieddJIuvfRSHX/88Tr88MPVM7ovZYrKysr4/UfKsh4IjxvX8SJgKC3nniv99Kfum4yvveYKi4YPl047LXagdMop7udnP5O+/nUXEre1uXYGZ5zRuf1Et24uXF6xouPjTU0ujNx77877qK5ur7yN1Njozm2iA95kgXCsuaOqKlgg7M0PDQ0dH/fGGavtxoABUm1t58eXL4/9zcDTT3f7eeKJ9kD4uefcxfruvbfzPqqq3Lnq97/vwuEtW1xF9uc/H3s8AIASUVHhPpWM9alrV3mTcFOTv0D41Vddr6lYn5zGU14eLBB+4w3pyCM7f70nEe99xLparCQ984z7RPaCC1xvZU7uAAA5NGTIEEnSunXrAgXCQXsIBwlr87lCWJI++eQTvfTSS7rlllu0evVqHXnkkXrooYf01a9+VXvuuafv7QHZQg9h5MS++7ofP+680xX63H23dNttrqJ13TrXVi+WESM6Vwhv2OCWsS44Hq9CuL7eBabRBS+R56aRklUIe/1/o18Ta/14gXCiCuFBg9z1YyJt2+Z+RozovH6vXu7bu7NnuwC4rMz1ZB41yl0jKJ6hQ13xEgAAGeeFp83NiXvoSq6K+D//kS6+ONg+glQIt7a6fkmxvu6WSOT7iA6EN22SrrzSXbmVMBgAkAe8kHPt2rUaO3Zs0vVTqazNRg/hbF1UTpLOPPNMvf/++xo3bpymTJmiL37xi7RVQF7jL07kvSOPdAUz99zjgt7f/95VE591Vuz1R47sXCHs9QgeNKjz+tXVLpj1egB7Ghrag9lI2WgZkUqF8KBBnXshe8chVssIyfUE3rTJhcLTpknvvOMquLkAHAAgL8T7FDaWRYvcp6/HHRdsH0EC4Q8/dGM5/PBg+0j0Pq67zrXaeOQR+u8AAPLC4MGDJblA2I9UK4TzLRBO5X0MHTpUkrRy5Urdf//9euWVV3TaaacRBiPvEQijIPzsZ+4bnWPHSrNmuXOneG0KRoxwLSkiA95EgbBXbRtdJdzQELsYKV4g7F0EPZOBcLIK4Y0bO14XxwuEY1UIS67txuGHu2v8XHWVuxBfvMprAACyLkgg/PLLbplKIGytvwb0Cxe65ac/HWwf3vvw/ljwvPGG9Nhj7qq6n/pUsG0CAJAhXiD8wgsv+GrrkMpF5YKGtd66+RYIH3LIIZo+fbree+89XXPNNfT0RcHg/1QUhP33d9dw+dSnpEmTpNtvj7/uyJGu0CeyfUKyCmGpcyDstYyI5s0N6aoQjrW+F0TX13d8PFmFcFubq/j1eK0z4lUId+sm/eMf7tx5/HjphRf4pioAII9EtlpI5uWX3dVlhw0Ltg9v4vMbCHfv7proB+G9j+g/BiZNkvr373ilXQAAcqxHjx664YYbNHPmTE2cODFpKJzNi8oFCVzLysoy3jLCGKMrrrhCw4cP9/0aIB/QQxgF48ADXSicjFcNu3x5exDqBcIDBnRe36u2jb6wXLKWEUF7CMcKhJuaUqsQjhcIS+69Dhzobi9f7tZN1MN+zz2lF1+M/zwAADnjt0K4rc0FwqecEnwf3ollKJT8U9G335YOPdTfBe4ixXof8+e7Cfiee9xVcwEAyCN33HGHevXqpVtuuUXdu3fXAw88ELcNQiqVtUHDWi8QLi8v9/2abFQIA4WKWkAUHS8EjuwjvH69u/h5rA/64lUIB20Zka0ewhUVsfcRGQh7Vqxwx4P2RQCAguS3QnjRIncF2a4Ewsm+EmutC4SDtouQYr+PX/9a6tvX9WwCACDPGGN0880368Ybb9RDDz2kyZMnx103ny8q56flhSeV9wEUKgJhFJ2hQ12/4WXL2h9bvz52uwgpfoVwtlpGxAqE47WM2LnTBdixAt5YgfDy5fH7BwMAkPf8Vgi/8IJbfvGLwffhVRolO2Fct85d/C2VXr/R72P1andF16uuin1hAAAA8sTtt9+u8847Tz/84Q/1j3/8I+Y6qV5ULkhYm68XlQMKFYEwik5lpQtBly5tfyxRIJyoQjhWIFxR4c4dgwbCsYqb4vUQrqx0+4nVMiLeeaPXJsILhNva3DE44IDY6wMAkPf8Vgg//7wLasMXwQnEb4Xw+++75SGHBN9H9PuYMcO1qPj2t4NvCwCALCorK9OMGTM0duxYXXDBBVoWWXkVls8Vwqn0EK6srPT9GqBQEQijKI0eLdXUtN9PNRCO1TJCclW9QXsIt7Z2PteM10NYcmF0rJYRsfoHS+7x6ur2QHjNGhc4jx4de30AAPKenwrh+nrp3/+WTjsttX0EDYQPPjj4PiLfh7XSI49IJ57I13gAAAWhurpaf/nLX2St1fjx4ztV9mbjonLeupmuEK6qqorbKxkoJgTCKEpeIGyt+0lnywjJhbhBK4Qj1/HEaxkhuTA6XsuIeAYObA+EvUCcQBgAULD8BMLPPecqb884I7V9+A2E33tP2ntv9xOU9z6am6VXX5U++ki64org2wEAIEdGjRql3/zmN3rllVf0u9/9rsNzzc3NqqyszGhYm62WEbSLQKkgEEZRGj3ahamffOKC3sbG9LWMkLITCMerEE7UanDQIAJhAEAR8dMyYvZsqX9/6fjjU9tHkArhVKqDpfb30dQkzZrlJv+vfjW1bQEAkCOXX365zjjjDE2aNEnLly/f/bhXWRtEWVlZxgPhoPtobm7mgnIoGQTCKEpeCFpT0x6QxguEvbYQkRXC1iZuGdGtW/xAONb8EXkeGCleD2EpdiCcrEI4OhDu2TO1dooAAOSFZBXCu3ZJzzwjfeUrrvl+KrwTy0QnjNa6QDiV/sFSxyvSPvWUa2/BxeQAAAXGGKMHH3xQ5eXl+sY3viFrraTUKmtTrRAu9y4Gm4F9UCGMUkIgjKIUGQh//LG7ve++sdctK3PBb2SF8K5d7twvUYVw0B7CketEviZIy4hkFcLDhkkrV0otLe69jx4t0f4IAFCwklUIP/ecm8DPPz/1fXgnlokqhFetcvvpaoXw/PnS6tXSueemth0AAHJs33331V133aV//vOfeuqppySlVllbXl7eqRdxIqm2jAiyDyqEUUoIhFGUhg51QWtNjWvVV1YmHXlk/PV79+4YCHuVuYkuKtfVlhGhkAtug7SMSFYhfNRRrj3GokXuvY8ZE39dAADyXrIK4alT3ddjTjop9X34aRnRlQvKSe3v44kn3P7OOiu17QAAkAeuuuoqHXTQQbrhhhvU0tKScoVwW1vb7irjZOghDKQXgTCKUlmZdMAB0pIlLhA+9NDEQWp1dceWEV5lbrwK4UQtI/wGwonW9/Ydq0I40fs49li3nDtXWr6c/sEAgAKXqEJ49Wrp2Welr39dqqxMfR/ZCIS997F2rZus+/dPbTsAAOSBiooK/eIXv9BHH32khx56KOUKYUm+K3izEQg3NzcTCKNkEAijaB19tPTPf7pA2AtK44lXIRzkonLeuWqsc9JYgbD3+kQtIyIrhEMhV/2bqGXEvvu6nsG/+IU7rz3qqPjrAgCQ9xJVCE+Z4vo7feMbXduHn0D4vfekffaR9twztX1Enlyedlpq2wAAII+ceeaZOvHEE/WTn/xEtbW1KVUIS/Id2HrrBQ2EJf+hcyoXxwMKFYEwitZ3v+sC1fp66ZhjEq8bXSHsp2VErH7A3brF7tmbqELYb8sIL7BOVCFsjAu/N2yQ9t9fOuOM+OsCAJD34gXCW7dKv/2t9OUvSyNHdm0ffiuEU60OljoGwqeemvp2AADIE8YY/fKXv1Rtba1efPHFjAfCqVQIe+v63QctI1BKCIRRtD71Ken0093tZIFwdIVwspYR8XoIx5s7UqkQjm4Z4Y0v2UXJvff6gx+0XycHAICCVFHhltEtI371K2n7dum227q+D+/EMt7JYlub9MEH0iGHpL6PyGqjRBc1AACggHzmM5/ReeedJ0lqbGwM9NqgYa0XCJcHOMkNGjpzUTmUEgJhFLXJk6U775RGjUq8XrwK4aA9hJMFwpHns97r470mumWEnwphSbr8cunmm6UJExKvBwBA3jPGTZSRn6guXizdfbd04YXSYYd1fR/eiWW8CuHly13Ppq5UCHvBdvRtAAAK3E033SRJ+uCDDwK9LhsVwqm0jKBCGKWCv0hR1MaMkW64Ifl61dWxewgnahmRjQrhpiZXsFRe3h5YJ6sQ3nvv9BRMAQCQF6qq2j9Rra+XLrnETdyTJ6dn+8laRnT1gnKSC7b/+7+ls85KfRsAAOShww8/XOPHj9d+++0X6HXZuqicRIUwEAuBMKDUWkbE6yEcS6o9hCUXTkeOL1mFMAAARcWrEF69WrrgAmnRIukvf5EGDUrP9pMFwl7F09ixXdvPXXd17fUAAOSpGTNmBH5NNiuE6SEMdEbLCEDtLSOsdfeTtYzIRoWwV53sjcWrECYQBgCUlKoq6Q9/kPbbT1q4UJo9O72VtskC4SVLXPjct2/69gkAQIkLGtZ662UyEK6vrycQRsmgQhiQO8dra3Pha2Tv3ngtI1LtIRwrEI73msgKYUnascMt+/SJ/z4AACg6l10mLVggHXqodN110ogR6d2+n0B4zJj07hMAgBKXbxXCr7/+ulavXq2jjjrK9/aBQkYgDEjaYw+33LrVhcB+WkY0N7uKYmPcY01NHS8iHsl7PJWWEd5Ytm51y379Er8XAACKSqZbLSQLhGtqpPAV1AEAQHp4wW4mA+Eg+7jvvvtUXV2tK664wvf2gUJGywhA7YHwtm1u2dDgLgIeL+D1QtzogDeTLSO8sXljBQAAaeCdWMY6Wdy82f1QIQwAQFqlWiHsvS6d+6itrdWsWbM0fvx49eEruSgRBMKA2qtuIwPheNXBUnvAG9k2ItMtI7Ztc49x0VMAANLIO7GMVSFcU+OWBMIAAKSVF9a2xfuGTpSutIxIto8pU6aoublZ1157re9tA4WOQBhQx5YRkmvTkCgQ9qp60xEI+20ZsW0b1cEAAKRdopYRS5a45ejR2RsPAAAlIF96CLe2tur3v/+9TjnlFB100EG+tw0UOgJhQLFbRsS7oJwUu2VEc3P8QLiszLWgCNJDOLplxNatBMIAAKRdokC4psZN4CNHZndMAAAUuXwJhOfMmaM1a9ZQHYySQyAMKHjLiKAVwpJ7LpUK4ciWEVxQDgCANEtWITxqlFRZmd0xAQBQ5IIGwt566Q6Ep06dqiFDhuiss87yvV2gGBAIA5L69nVLvy0jgvYQ9l7T3Nx+328PYVpGAACQQckCYfoHAwCQdtmoEPbWjbePtWvX6rnnntOECRMCXawOKAYEwoDct0Grq9srhLdvlxJdXDRdFcIVFe3Xsonm7X/7drfcupUKYQAA0i5eIBwKSUuXEggDAJABycLaaF4gHCS4TRY6T58+XW1tbZo4caLvbQLFgkAYCOvXr71CeNMmae+9468bq4dw0EC4qSl+uwjJVQj37OnGIlEhDABARngnltEni6tWucmaC8oBAJB22ewh3BbjW0BtbW16+OGHdeKJJ2q//fbzvU2gWBAIA2F77NFeIVxbK/XvH3/ddFUIJwqEJTeG2lpXtLR9O4EwAABpF69CuKbGLakQBgAg7RKFtbF46xljAu8jVug8b948LVu2TFdeeaXv7QHFhEAYCPMC4aYmaceOxBXC0T2ErU0tEE60vuTGUFsr1dW581RaRgAAkGbxAuElS9ySCmEAANIulQrhINXByfYxdepU9e3bV1/96lcDbRMoFgTCQNgee7iWEZs3u/t+WkZ4gXBLi1umu0J4771dywivcpkKYQAA0ixRhXCfPtLAgdkfEwAARS5oIBwKhdIWCG/fvl2zZ8/WxRdfrB49egTaJlAsCISBsH79XPBaW+vu+2kZ4QW83rKqKv5rqqqC9RCW2iuECYQBAMiQRBXCY8ZIAb6aCgAA/MllhfDMmTO1a9cu2kWgpBEIA2FeywgvEA7SMsILetPdMsLrIexd7I6WEQAApFmiQJh2EQAAZEQ2AmFv/eh9TJ06VYceeqiOPPLIQNsDigmBMBC2xx7uwm0bNrj7QVpGpBoI+6kQbmiQ1q1rHyMAAEgj7+Qy8mSxoUFavZoLygEAkCHxwtp42tradofIfsUKnd9991298cYb+vrXvx7oAnVAsSEQBsK86tuPP3ZLPy0jshEIS9LSpR3HCAAA0sQ7uYysEP7oI7ekQhgAgIzIZsuItog5fvr06aqsrNSll14aaFtAsSEQBsK86tuPPnLtAvfcM/660YFwc7NbZqJlhDemyDECAIA0idUyoqbGLakQBgAgI2KFtYmko4dwS0uLHnvsMZ111lnqn6gCDCgBBMJAWGQgvNde7QVDsVRUSJWVUn29u++3QtgLjiX32l69Eo/JqxD2Quo+fRKvDwAAAooVCC9Z4pYHHJD98QAAUAJycVG55557Ths3btTEiRMDbQcoRgTCQJgXvn7wQeJ2EZ5evYIHwpEVwkEC4Q8+cBXLAec/AACQTLwK4aFDk0/UAAAgJUED4VAo1OVAePr06RowYIBOP/30QNsBihHxEhB25JFSdbVUV5f4gnKebAbCdXXSSSclHxMAAAgoXoUw7SIAAMiYbFcIb9q0SU8//bQuu+wyVVZWBhssUIQIhIGw7t0l74NCP4FwdXXXA+Hq6sT76Nu3vXXFOeckHxMAAAgoOhC21gXCXFAOAICM8cLdIIFwQdmfxgAADShJREFUeaK+jkn2MXPmTLW0tOiKK64INlCgSBEIAxHOPdctg7aMaGhofyzR+g0N7jwzFHLhcLIK4bIyN5bycunMM5OPCQAABOSdXHonpLW10vbtVAgDAJBB2a4QnjZtmo444ggdeuihwQYKFKmKXA8AyCdf+pKrFB42LPm6kYHwzp1umajit7rahcGNjVJra/s2khk+XDrsMKlfv+TrAgCAgKIrhGtq3JIKYQAAMsYLa9siWzYl0JVAeOHChVq4cKF++9vfBhskUMQIhIEI/fpJ77wjDRmSfN1evaTNm91tv4Gwt673IaifQPiPf0zcigIAAHRBvECYCmEAADImmxXCM2bMUGVlpS6++OJggwSKGIEwEMVvQVCvXtKqVe52kEC4rs5VCnvbSGboUH/jAQAAKYgVCFdWuq/oAACAjAgaCIdCoZQD4XXr1um8885Tfz+9IYESQSAMpChWy4hEAW9khXDkNgAAQA5FB8JLlkj779/eWxgAAKRdNiuEJWnChAmBXgsUOwJhIEW9erWHuzt3Sj16JD53JBAGACAPxaoQpn8wAAAZlY1A2Ft/wIABOv3004MNEChywX6bAOwWXSGcqF2E1DEQ9l5HIAwAQI55n+a2tbkm/0uXEggDAJBhXlgbJBAuD/jtnW7hi/FcdtllqqysDDZAoMhRIQykqLpa2rXLnTsGDYS9DzaTvQYAAGSYNymHQu7iAM3NBMIAAGRYNiqEBw4cqBkzZujLX/5y4PEBxY5AGEiRV93b0BA8EPY+2KRCGACAHItsGVFT424TCAMAkFFeINzmtWxKIpVAWJLGjx8f+DVAKSAQBlLkhbn19cED4YqKjtsAAAA5QiAMAEDWeYFwU1OTr/VTDYQBxMZvE5CirgTC9BAGACBPRAfCvXtLAwfmdkwAABS53r17a/To0frTn/7ka/1QKEQgDKQRv01AioIGwt27u3NOAmEAAPJIdCA8erRkTG7HBABAkTPG6Dvf+Y7mz5+v+fPnJ12fCmEgvfhtAlIUNBA2xq3jBcJVVe2tIwAAQI5EB8JjxuR2PAAAlIgrrrhCffr00W9+85uk67a1te1uMwGg6wiEgRQFDYSljoEw1cEAAOQB7+Syvl5auZL+wQAAZEnv3r115ZVX6sknn9TatWsTrkuFMJBe/DYBKYrsCRw0EPa7PgAAyDDv5PKjjyRrCYQBAMiia6+9VqFQSPfff3/C9QiEgfTitwlIkVfhu2OH1NhIhTAAAAXJO7n88EO3JBAGACBrRo0apbPPPlsPPvigdu7cGXc9AmEgvfhtAlLkBbobN7olgTAAAAUoOhA+4IDcjQUAgBJ0ww03aPPmzXrggQfirhMKhQiEgTTitwlIkRfobtjglgTCAAAUIO/ksr5eGjRI6tMnt+MBAKDEHHPMMTr11FP1y1/+UvX19THXoUIYSC9+m4AURQfCfgLeXr0IhAEAyCvGtN+mXQQAADlx6623auPGjZo8eXLM5wmEgfTitwlIUWWl+6FCGACAAlde7pYEwgAA5MSxxx6rr3zlK7rzzju1fv36Ts+3tbWp3JuvAXQZgTDQBb16EQgDAFDwvIojAmEAAHLmzjvv1K5du3T99dd3eo4KYSC9+G0CuqC6OvVA2M/6AAAgCwiEAQDIudGjR+vGG2/UY489pr/97W8dniMQBtKL3yagC1KpEG5tlbZsoUIYAIC8QSAMAEBe+J//+R+NHTtWEydO1Lp163Y/TiAMpBe/TUAX9OolWetu+w2EJfcaAmEAAPJEWZn7GTUq1yMBAKCkdevWTU8++aR27typr3zlK6qrq5MkhUIhAmEgjfhtArpg2LD220ECYUnad9/0jwcAAKSgrEwaMULq1i3XIwEAoOSNHTtWjz32mN58802deuqpWrt2LRXCQJrx2wR0wa23tt/2EwhXVLTfvvTStA8HAACkoqxMGjMm16MAAABh5557rmbPnq1FixZp7NixWrZsGYEwkEb8NgFdcNhh0llnudtVVcnX9841n3zS3/oAACALhg6Vxo3L9SgAAECEc889V++++65OPvlkNTY2avDgwbkeElA0KpKvAiCROXOkhgbJmOTrHn20tHMn/YMBAMgrb70llZfnehQAACDKfvvtpz//+c9qaGhQN1o7AWlDIAx0kTHBAl7CYAAA8gxf2wEAIK/17Nkz10MAigotIwAAAAAAAACgRBAIAwAAAAAAAECJIBAGAAAAAAAAgBJBIAwAAAAAAAAAJYJAGAAAAAAAAABKhLHW+l/ZmFpJK9Ow3/6SNqVhO6WAY+Ufx8ofjpN/HCt/0nmchltr907TtkpWGudrid8DvzhO/nGs/OE4+cex8i9dx4r5Og2Yr3OGY+UPx8k/jpU/HCf/Mn6OHSgQThdjzBvW2s9kfccFiGPlH8fKH46TfxwrfzhOxY3/vv5wnPzjWPnDcfKPY+Ufx6p48d/WP46VPxwn/zhW/nCc/MvGsaJlBAAAAAAAAACUCAJhAAAAAAAAACgRuQqEH8rRfgsRx8o/jpU/HCf/OFb+cJyKG/99/eE4+cex8ofj5B/Hyj+OVfHiv61/HCt/OE7+caz84Tj5l/FjlZMewgAAAAAAAACA7KNlBAAAAAAAAACUiIwGwsaY040xS4wxS40xk2I8380Y80T4+fnGmBGZHE8+83Gs/ssY84ExZpEx5kVjzPBcjDPXkh2niPXON8ZYY0zJXsHSz7EyxlwQ/v/qfWPM49keYz7w8bs3zBjzkjFmYfj370u5GGc+MMY8bIzZaIx5L87zxhjz2/CxXGSMOSLbY0RqmK/9Y772jznbH+Zr/5iz/WG+Lm7M2f4wX/vHfO0fc7Y/zNf+5Hy+ttZm5EdSuaSPJY2SVCXpHUljo9b5lqQHwrcvkvREpsaTzz8+j9VJknqGb19TisfKz3EKr9db0jxJr0n6TK7Hna/HStIBkhZK6he+PyDX487T4/SQpGvCt8dKWpHrcefweB0v6QhJ78V5/kuS/ibJSBonaX6ux8yPr/+uzNfpPVYlP1/7PVbh9Up6zma+TvuxYs62zNfF/MOcndbjxHzt81iF1yvp+drvsWLOZr4OeKxyOl9nskL4KElLrbXLrLXNkmZJOidqnXMkzQjfni3pC8YYk8Ex5aukx8pa+5K1tiF89zVJQ7M8xnzg5/8pSfqppF9I2pXNweUZP8fqKkm/s9ZulSRr7cYsjzEf+DlOVlKf8O2+ktZlcXx5xVo7T9KWBKucI+kR67wmaQ9jzD7ZGR26gPnaP+Zr/5iz/WG+9o852yfm66LGnO0P87V/zNf+MWf7w3ztU67n60wGwkMkrY64vyb8WMx1rLWtkrZL2iuDY8pXfo5VpCvlPiUoNUmPkzHmcEn7WmufyebA8pCf/6dGSxptjHnFGPOaMeb0rI0uf/g5TrdKuswYs0bSs5K+k52hFaSg/5YhPzBf+8d87R9ztj/M1/4xZ6cP83XhYs72h/naP+Zr/5iz/WG+Tp+MztcV6dpQDLE+hbQprFMKfB8HY8xlkj4j6YSMjig/JTxOxpgySb+WNCFbA8pjfv6fqpD7SsuJcp+Iv2yMOcRauy3DY8snfo7TxZKmW2vvMcYcI+nR8HFqy/zwCg7/phcm5mv/mK/9Y872h/naP+bs9OHf9MLFnO0P87V/zNf+MWf7w3ydPhn99zyTFcJrJO0bcX+oOpeB717HGFMhVyqeqFy6WPk5VjLGnCLpfySdba1tytLY8kmy49Rb0iGS/mWMWSHXY2VOiTa99/v79xdrbYu1drmkJXKTVynxc5yulPQHSbLWviqpu6T+WRld4fH1bxnyDvO1f8zX/jFn+8N87R9zdvowXxcu5mx/mK/9Y772jznbH+br9MnofJ3JQPh1SQcYY0YaY6rkGtrPiVpnjqQrwrfPl/RPG+6cXGKSHqvw1zQelJusSrEPjZTkOFlrt1tr+1trR1hrR8j1gjrbWvtGboabU35+/56Su5iCjDH95b7esiyro8w9P8dplaQv6P+3d4c4VkRBGEb/Ejgkm2AHrGAEAo0i4DFkgibBsgESAhaBGzd7wIIahUEQBAYFKUQ/T5HMm57hnrOCyjXfS73u20mq6m62WH270ilvjrMkjw5fQ72X5Ed3f917KP5Kr+f0ek6zZ/R6TrMvj17fXJo9o9dzej2n2TN6fXmO2uujXRnR3b+q6mmS82xfGXzX3Z+q6mWSj919luRttkfDL7L9a/nwWPNcZ8OzepXkdpIPh28CfOnuB7sNvYPhOZHxWZ0nOamqz0l+J3ne3d/3m/rqDc/pNMmbqnqW7fWMxwv+qE6SVNX7bK8/3Tnc9/Qiya0k6e7X2e5/up/kIsnPJE/2mZR/oddzej2n2TN6PafZc3r9/9LsGb2e0+s5zZ7R67m9e10LnjkAAAAAwJKOeWUEAAAAAADXiIUwAAAAAMAiLIQBAAAAABZhIQwAAAAAsAgLYQAAAACARVgIAwAAAAAswkIYAAAAAGARFsIAAAAAAIv4A+lfedlUOdENAAAAAElFTkSuQmCC\n",
"text/plain": [
""
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
}
],
"source": [
"N = 300 # number of discrete samples per signal\n",
"\n",
"## Define a set of signals\n",
"t = np.linspace(0.,1.,N)\n",
"\n",
"s_template = []\n",
"s_template.append(signal_gabor(t)) \n",
"s_template.append(signal_sawtooth(t)) \n",
"s_template.append(signal_square(t))\n",
"\n",
"## generate confounds\n",
"num_classes = len(s_template) # number of classes\n",
"\n",
"## Plotting\n",
"fig, ax = plt.subplots(1, 3, sharex=False, sharey=False, figsize=(25,3))\n",
"\n",
"#c = ['b*', 'ro', 'kx'][i]\n",
"ax[0].plot(t,s_template[0],'b')\n",
"ax[0].set_title(\"Gabor wave\")\n",
"ax[0].set_yticks([])\n",
"\n",
"ax[1].plot(t,s_template[1],'r')\n",
"ax[1].set_title(\"Sawtooth wave\")\n",
"ax[1].set_yticks([])\n",
"\n",
"ax[2].plot(t,s_template[2],'k')\n",
"ax[2].set_title(\"Square wave\")\n",
"ax[2].set_yticks([])\n",
"\n",
"plt.show()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Generate confounds (apply 4-th degree polynomial)"
]
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {},
"outputs": [],
"source": [
"Lp = 4\n",
"p0 = np.linspace(-0.5,0.15,Lp) #[-0.5, 0.2]\n",
"p1 = np.linspace(0.85,1.25,Lp) #[0.75, 1.5]\n",
"p2 = np.linspace(0.5, 1.25,Lp) #[0.5, 1.5]\n",
"p3 = np.linspace(0.5, 1.25,Lp) #[0.5, 1.5]\n",
"p4 = np.linspace(0.5, 1.25,Lp) #[0.5, 1.5]\n",
"\n",
"s_conf = []\n",
"y_conf = []\n",
"coeff = np.zeros(5)\n",
"for k in range(num_classes):\n",
" Lc = 0\n",
" for a in range(Lp):\n",
" coeff[0]=p0[a]\n",
" for b in range(Lp):\n",
" coeff[1]=p1[b]\n",
" for c in range(Lp):\n",
" coeff[2]=p2[c]\n",
" for d in range(Lp):\n",
" coeff[3]=p3[d]\n",
" for e in range(Lp):\n",
" coeff[4]=p4[e]\n",
" g = coeff[4]*(t**4)+coeff[3]*(t**3)+coeff[2]*(t**2)+coeff[1]*(t**1)+coeff[0]*(t**0)\n",
" g_prime = 4*coeff[4]*(t**3)+3*coeff[3]*(t**2)+2*coeff[2]*(t**1)+coeff[1]\n",
" if k==0:\n",
" sc = g_prime*signal_gabor(g) \n",
" elif k==1:\n",
" sc = g_prime*signal_sawtooth(g) \n",
" elif k==2:\n",
" sc = g_prime*signal_square(g) \n",
" s_conf.append(sc)\n",
" y_conf.append(k)\n",
" Lc = Lc+1\n",
"s_conf = np.stack(s_conf, axis=0)\n",
"y_conf = np.stack(y_conf, axis=0)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Split data into train and test sets"
]
},
{
"cell_type": "code",
"execution_count": 6,
"metadata": {},
"outputs": [],
"source": [
"def random_index_matrix(max_index, n_samples_perclass, num_classes, repeat, y_train):\n",
" seed = int('{}{}{}'.format(n_samples_perclass, num_classes, repeat))\n",
" np.random.seed(seed)\n",
" tr_index = []\n",
" te_index = []\n",
" for classidx in range(num_classes):\n",
" max_samples = (y_train == classidx).sum()\n",
" tr_index.append(np.random.randint(0, max_samples, (n_samples_perclass)))\n",
" te_index.append(~tr_index[classidx])\n",
" return tr_index, te_index"
]
},
{
"cell_type": "code",
"execution_count": 7,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"(1536, 300)\n",
"(1536, 300)\n"
]
}
],
"source": [
"data_train = []\n",
"data_test = []\n",
"y_train = []\n",
"y_test = []\n",
"\n",
"max_index = len(y_conf) // num_classes\n",
"train_index, test_index = random_index_matrix(max_index, max_index//2, num_classes, num_classes, y_conf)\n",
"for i in range(num_classes):\n",
" data_class = s_conf[np.where(y_conf==i)]\n",
" y_class = y_conf[np.where(y_conf==i)]\n",
" data_train.append(data_class[train_index[i],:])\n",
" y_train.append(y_class[train_index[i]])\n",
" data_test.append(data_class[~train_index[i],:])\n",
" y_test.append(y_class[~train_index[i]])\n",
"data_train = np.concatenate(data_train,axis=0)\n",
"data_test = np.concatenate(data_test,axis=0)\n",
"y_train = np.concatenate(y_train)\n",
"y_test = np.concatenate(y_test)\n",
"print(data_train.shape)\n",
"print(data_test.shape)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"In this example we have used 256 randomly chosen samples per class to train the model. We have used another function *take_train_samples* function from *utils.py* script for this. User can use their own script."
]
},
{
"cell_type": "code",
"execution_count": 8,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"(768, 300)\n"
]
}
],
"source": [
"n_samples_perclass = 256\n",
"x_train_sub, y_train_sub = take_train_samples(data_train, y_train, n_samples_perclass, num_classes, repeat=0)\n",
"print(x_train_sub.shape)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Create an instance for SCDT_NS class"
]
},
{
"cell_type": "code",
"execution_count": 9,
"metadata": {},
"outputs": [],
"source": [
"scdt_ns_obj = SCDT_NS(num_classes, rm_edge=True)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Training Phase"
]
},
{
"cell_type": "code",
"execution_count": 10,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"\n",
"Calculating SCDTs for training data ...\n",
"Generating basis vectors for each class ...\n"
]
}
],
"source": [
"scdt_ns_obj.fit(x_train_sub, y_train_sub)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Testing Phase"
]
},
{
"cell_type": "code",
"execution_count": 11,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"\n",
"Calculating SCDTs for testing samples ...\n",
"Finding nearest subspace for each test sample ...\n"
]
}
],
"source": [
"preds = scdt_ns_obj.predict(data_test, use_gpu=False)"
]
},
{
"cell_type": "code",
"execution_count": 12,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"\n",
"Test accuracy: 100.0%\n"
]
}
],
"source": [
"from sklearn.metrics import accuracy_score\n",
"\n",
"print('\\nTest accuracy: {}%'.format(100*accuracy_score(y_test, preds)))"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": []
}
],
"metadata": {
"kernelspec": {
"display_name": "Python 3",
"language": "python",
"name": "python3"
},
"language_info": {
"codemirror_mode": {
"name": "ipython",
"version": 3
},
"file_extension": ".py",
"mimetype": "text/x-python",
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.7.7"
}
},
"nbformat": 4,
"nbformat_minor": 4
}