diff --git a/examples/ShapeAnalysis.ipynb b/examples/ShapeAnalysis.ipynb new file mode 100644 index 0000000..0371d17 --- /dev/null +++ b/examples/ShapeAnalysis.ipynb @@ -0,0 +1,406 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Shape Analysis\n", + "\n", + "After finding atlas correspondence points across a shape population we can find insights into how shape informs other metrics and classes. In this example we train a Distance-Weighted Discrimination (DWD) classifier to predict whether a given mouse femur is healthy/unhealthy based on its shape features.\n", + "\n", + "DWD classification is designed to work with High-Dimensional Low Sample Size (HDLSS) data where the number of features for a given sample significantly outnumbers the total number of samples in the data set. In this case we have 28 mouse femur samples, each of which is represented by approximately 4000 points in three-dimensional space. We use DWD to get a hyperplane separating the feature space for healthy and unhealthy femur classes and analyze prediction accuracy and distance to the hyperplane in order to assess performance.\n", + "\n", + "This notebook assumes that a population of shapes in correspondence is available for training the DWD classifier. See the `TemplateGenerationIterative` notebook for a procedure to create a representative atlas from a shape population and the `MeshToMeshRegistration` notebook for examples of getting correspondence points on individual samples with the generated atlas." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "scrolled": true + }, + "outputs": [], + "source": [ + "import sys\n", + "!{sys.executable} -m pip install itk dwd sklearn seaborn matplotlib pandas" + ] + }, + { + "cell_type": "code", + "execution_count": 2, + "metadata": {}, + "outputs": [], + "source": [ + "import os\n", + "import glob\n", + "\n", + "from dwd.dwd import DWD\n", + "import itk\n", + "import numpy as np\n", + "import sklearn.model_selection\n", + "import pandas as pd\n", + "\n", + "module_path = os.path.abspath(os.path.join('..'))\n", + "\n", + "if module_path not in sys.path:\n", + " sys.path.append(module_path)\n", + "\n", + "from src.hasi.hasi import classify" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Load Correspondence Meshes" + ] + }, + { + "cell_type": "code", + "execution_count": 3, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "28\n" + ] + } + ], + "source": [ + "CORRESPONDENCE_INPUT = 'Output/correspondence/'\n", + "\n", + "paths = glob.glob(CORRESPONDENCE_INPUT+'*')\n", + "print(len(paths))" + ] + }, + { + "cell_type": "code", + "execution_count": 4, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "28 shapes found\n", + "14 healthy samples, 14 unhealthy samples\n", + "Meshes each have 3817 points\n" + ] + } + ], + "source": [ + "meshes = [itk.meshread(path, itk.F) for path in paths]\n", + "\n", + "# Get femur class from filename: right femurs are healthy, left feurs are unhealthy\n", + "labels = np.array(['Healthy' if '-R' in path else 'Unhealthy' for path in paths])\n", + "\n", + "\n", + "print(f'{len(meshes)} shapes found')\n", + "print(f'{len(labels[labels == \"Healthy\"])} healthy samples, '\n", + " f'{len(labels[labels == \"Unhealthy\"])} unhealthy samples')\n", + "\n", + "assert(not any(mesh.GetNumberOfPoints() != meshes[0].GetNumberOfPoints()\n", + " for mesh in meshes[1:]))\n", + "print(f'Meshes each have {meshes[0].GetNumberOfPoints()} points')" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Prepare Data\n", + "\n", + "DWD expects input of size `n x d` where `n` = number of samples and `d` = number of features. We can generate the shape features for a given mesh by flattening the list of 3D points to one dimension, however for our given population this exceeds DWD's memory availability. For this example we use a step size of 10 to select every tenth point for inclusion in the feature array, leading to a feature length of (int(3817 / 10) * 3) ~= 1146 for each sample." + ] + }, + { + "cell_type": "code", + "execution_count": 5, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "(28, 1146)\n" + ] + } + ], + "source": [ + "features = classify.make_point_features(meshes, step=10)\n", + "print(features.shape)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "We use `sklearn` to split the data set into samples for training and testing." + ] + }, + { + "cell_type": "code", + "execution_count": 6, + "metadata": {}, + "outputs": [], + "source": [ + "X_train, X_test, y_train, y_test = sklearn.model_selection.train_test_split(features,labels, train_size=0.6, random_state=73)" + ] + }, + { + "cell_type": "code", + "execution_count": 7, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "(16, 1146)\n", + "(16,)\n", + "(12, 1146)\n", + "(12,)\n" + ] + } + ], + "source": [ + "print(X_train.shape)\n", + "print(y_train.shape)\n", + "print(X_test.shape)\n", + "print(y_test.shape)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Train Classifier\n", + "\n", + "Here we fit the DWD classifier to the training set and then verify with the testing set." + ] + }, + { + "cell_type": "code", + "execution_count": 8, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "DWD(C=9.989272474981702)" + ] + }, + "execution_count": 8, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "classifier = DWD(C='auto')\n", + "classifier.fit(X_train,np.squeeze(y_train))" + ] + }, + { + "cell_type": "code", + "execution_count": 9, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "83.33% correct\n", + "[ True True True False False True True True True True True True]\n" + ] + } + ], + "source": [ + "predict = classifier.predict(X_test)\n", + "correct = np.squeeze(predict) == y_test\n", + "\n", + "print(f'{float(sum(correct) / len(correct)):0.2%} correct')\n", + "print(correct)" + ] + }, + { + "cell_type": "code", + "execution_count": 10, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "[ 0.51246972 1.37529028 -0.88142607 1.36234602 -0.78152306 -1.09616038\n", + " -0.1393085 0.32118888 -0.90293256 -0.69205907 -1.07170148 -0.70416412]\n" + ] + } + ], + "source": [ + "# Print distances to separating hyperplane\n", + "print(classify.get_distances(classifier, X_test))" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Visualize Results\n", + "\n", + "We can examine the fitness of the hyperplane by comparing the distance to the hyperplane for each test sample." + ] + }, + { + "cell_type": "code", + "execution_count": 11, + "metadata": {}, + "outputs": [ + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYIAAAD4CAYAAADhNOGaAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/Z1A+gAAAACXBIWXMAAAsTAAALEwEAmpwYAABFiklEQVR4nO3dd3gU5drA4d+bnpCQQAKE3qS3QAKEAKEqRVAQpFpAEBSxocdPOKjYDlawHFGKCKI0BQ6CgHSlBRLpvQdCJ4GE9PZ+f8wmBAwQYDeTZJ/7uvaaLbMzz26y88y8VWmtEUIIYb8czA5ACCGEuSQRCCGEnZNEIIQQdk4SgRBC2DlJBEIIYeeczA7gbnXu3FmvWLHC7DBEEdF2RlsA1g9ab2ocQuQDdasXCt0VweXLl80OQQghipRClwiEEEJYV6ErGhLCmsaGjjU7BCFMJ4lA2LWO1TqaHYIQppNEIOzazvM7AQjwDzA1DnuUlpZGVFQUycnJZodSpLi5uVGhQgWcnZ3z/B5JBMKuvbLiFUBaDZkhKioKLy8vqlSpglK3bNAi7oLWmujoaKKioqhatWqe3yeVxUIIUyQnJ+Pr6ytJwIqUUvj6+t71VZYkAiGEaSQJWN+9fKeSCIQQws5JIhBC2C1PT88bHs+YMYORI0fe07bWr19Pt27dsu9v3rw5+7VBgwbx66+/3nugNiaVxcKu/afDf8wOQRRB69evx9PTk5CQELNDyRO5IhB2LaRiCCEVC8ePVeSvS5cu0atXL5o2bUrTpk3ZtGkTANu2baNFixY0btyYkJAQDh06dMP7Tp48yXfffcfEiRMJCAhgw4YNAPz111+EhIRQrVq17KuDp556iv/973/Z7x04cCCLFy/Onw+Yg1wRCLu2+bRx+S7JwFzvLtnH/rNxVt1m3XLFead7vduuk5SUREBAQPbjmJgYHnnkEQBefvllXn31VVq1asWpU6fo1KkTBw4coHbt2mzYsAEnJydWr17NmDFjWLBgQfY2qlSpwnPPPYenpyevv/46AN9//z3nzp1j48aNHDx4kEceeYTevXszZMgQJk6cSI8ePYiNjWXz5s3MnDnTqt9DXkgiEHZtzJoxgPQjsFfu7u7s3Lkz+/GMGTOIiIgAYPXq1ezfvz/7tbi4OOLj44mNjeXpp5/myJEjKKVIS0vL07569OiBg4MDdevW5cKFCwC0adOGESNGcOnSJRYsWECvXr1wcsr/w7IkAiGE6e505m6GzMxMwsLCcHNzu+H5kSNH0q5dOxYtWsTJkydp27Ztnrbn6uqafV9rnX3/qaee4qeffmLu3Ln88MMPVon9bkkdgRBC5OKhhx7i66+/zn6cdeUQGxtL+fLlAeMKIjdeXl5cu3YtT/sZNGgQX3zxBQB169a953jvhyQCIYTIxVdffUVERAQNGzakbt26fPfddwC88cYbjB49msaNG5Oenp7re7t3786iRYtuqCy+lTJlylCnTh0GDx5s9c+QVyrnJUphEBQUpLPK8IS4XzJDmXkOHDhAnTp1zA7DdImJiTRo0IDt27fj7e1tlW3e4ru9ZZdjqSOwstTUVE6cOMGJEyc4efIkkZGRxMfHk5KSQmpqKpmZmfj4+ODn54efnx9ly5alfv36VKtWDUdHR7PDtztfdP7C7BCEHVu9ejVDhgzh1VdftVoSuBeSCO5TQkICf/75Jxs3bmTjxo1s27aNlJSU7NednZ3x8vLCxcUFFxcXHBwcuHr1KlevXr1hO+7u7tSrV49mzZrRsWNH2rVrh4+PT/5+GDskw08LM3Xs2JHIyEizw5BEcC8SExNZtmwZ8+fPZ+nSpSQlJeHk5ESTJk0YMWIEjRs3pmrVqlSpUoVy5crh4PDPqpi0tDRiYmI4deoUe/fuZc+ePezevZuZM2cyadIkHBwcaNasGb169aJ///7ZlVPCulYfXw3IBDXCvkkdwV04evQoX375JTNmzCA+Pp4yZcrQu3dvevbsSYsWLfDw8LjvfaSmphIWFsaqVatYsWIFERERKKVo3749TzzxBH379sXd3d0Kn0aA1BGYSeoIbOdu6wik1VAebNy4kZ49e1KzZk2mTJnCY489xpo1azhz5gz//e9/6dChg1WSAICLiwuhoaG8//77hIeHc/jwYd5++20iIyMZPHgwlSpV4q233uL8+fNW2Z8QQkgiuI19+/bRrVs3WrduzYYNGxg7diyRkZHMnDmT9u3b50vlbo0aNRg3bhyHDx9m/fr1tGrVig8//JBKlSrxzDPPcOLECZvHIIQo2iQR5OLs2bMMHTqUhg0bsnHjRj755BNOnTrFe++9h7+/vykxKaVo06YNixYt4vDhwwwfPpw5c+ZQq1YtRowYwdmzZ02JS4jC6uTJk9SvX/+G58aNG8dnn312y/fczzDVN8s5NPUXX3xBYmJi9ms3D49ta5IIctBaM2XKFOrUqcOPP/7ISy+9xLFjx/jXv/5ltaIfa3jggQf4+uuvOXbsGEOHDmXq1KlUr16dMWPGEB8fb3Z4Qoi7dHMiyG+SCCyOHTtGhw4dGD58OIGBgezfv5+JEyfi6+trdmi3VK5cOSZNmsShQ4fo3bs348ePp06dOsyfP5/C1gjALJO7TWZyt8lmhyEKmLZt2/J///d/NGvWjJo1a97QO/js2bN07tyZGjVq8MYbb2Q/v3LlSlq0aEGTJk14/PHHs0/K3nvvPZo2bUr9+vUZNmzYP36bX331FWfPnqVdu3a0a9cu+/l///vfNGrUiODgYC5cuMC1a9eoWrVq9iB3cXFxNzy+H3bffFRrzeTJkxk1ahROTk5MnjyZZ599tlDNpVqtWjVmzZrF888/z8iRI+nbty/fffcdkyZNonbt2maHV6DV8qtldggCYPmbcH6Pdbfp3wC6fHTPb09PT2fbtm0sW7aMd999l9WrjabGO3fuZMeOHbi6ulKrVi1efPFF3N3d+eCDD1i9ejXFihXj448/ZsKECbz99tuMHDmSt99+G4Ann3ySpUuX0r179+z9vPTSS0yYMIF169bh5+cHGP2TgoOD+fDDD3njjTeYOnUqY8eOpW3btvz+++/06NGDuXPn8thjj+Hs7HwfX5LBrq8I4uLi6N+/P88//zytW7dm//79DBs2rFAlgZxCQkIIDw/n22+/ZefOnQQEBPDxxx/fcjwUAUsOLWHJoSVmhyFMcKvfedbzjz32GACBgYGcPHky+/UOHTrg7e2Nm5sbdevWJTIykrCwMPbv30/Lli0JCAhg5syZ2R3F1q1bR/PmzWnQoAFr165l3759d4zNxcUle9rLnPsfOnRo9gilP/zwg9XGJ7LbK4KdO3fSp08fjh8/zvjx43njjTdy7fhV2Dg6OvLcc8/Ro0cPXnjhBd58800WLFjA9OnT/1ExJuDzLZ8D0L1W9zusKWzqPs7c75Wvry9Xrly54bmYmBiqVq0KXB822tHR8YaTqZzDSWe9prXmwQcfZM6cOTdsLzk5mREjRhAREUHFihUZN24cycnJd4zN2dk5OyHl3H/Lli05efIk69evJyMjw2q/6cJ/5LsHs2fPJjg4mISEBNatW8ebb75ZJJJATv7+/vz666/MmzePEydOEBgYyIQJE8jMzDQ7NCEKBE9PT8qWLcvatWsBIwmsWLGCVq1a3fW2goOD2bRpE0ePHgWMop3Dhw9nH/T9/PyIj4+/5QT2dzNs9VNPPcWAAQOsOlpp0Tr63YHWmnfffZeBAwcSHBzMjh07aN26tdlh2YxSij59+rB//366dOnCa6+9RteuXaUzmhAWP/74I++//z4BAQG0b9+ed955h+rVq9/1dkqVKsWMGTPo378/DRs2pEWLFhw8eBAfHx+effZZ6tevT6dOnWjatGmu7x82bBidO3e+obL4VgYOHMiVK1fo37//Xcd5S1rrQnULDAzU9yIpKUkPGDBAA3rQoEE6JSXlnrZTWGVmZupvv/1Wu7m56VKlSunff//d7JAKhDY/tNFtfmhjdhh2af/+/WaHUCj98ssv+oknnrjtOrf4bm95XLWLK4IrV67QsWNHZs+ezfjx45k+fTouLi5mh5WvlFI899xzRERE4O/vz8MPP8y///1vMjIyzA5NCJFHL774Im+++SZvvfWWVbdb5CuLL1y4QKdOnThw4ADz5s2jT58+Zodkqnr16rF161Zeeukl/vOf/xAWFsacOXMoXbq02aGZYlbPWWaHIESe5Zw605psekWglOqslDqklDqqlHrzNuv1UkpppVSQNfd/6tQpQkNDOXLkCEuXLrX7JJDF3d2dqVOnMn36dDZv3kyTJk3YsmWL2WGZoqJ3RSp6VzQ7DLulpeOj1d3Ld2qzRKCUcgS+AboAdYH+Sql/zMyslPICXga2WnP/R44coVWrVpw/f56VK1fy4IMPWnPzRcLgwYPZsmULrq6utGnThunTp5sdUr6bt3ce8/bOMzsMu+Tm5kZ0dLQkAyvSWhMdHY2bm9tdvc+WRUPNgKNa6+MASqm5wKPA/pvWex/4GPiXtXZ85MgR2rRpQ1paGuvXr6dx48bW2nSRExAQQHh4OH379mXIkCHs3r2bzz77DCenIl9qCMC3Ed8C0Ld+X5MjsT8VKlQgKiqKS5cumR1KkeLm5kaFChXu6j22/LWXB07neBwFNM+5glKqCVBRa/27UsoqieDEiRO0b98+OwnUq1fPGpst0kqWLMny5cv517/+xRdffMG+ffuYN28eJUuWNDs0UYQ5Oztnd94S5jKt1ZBSygGYALyWh3WHKaUilFIRtzt7OH36NB06dCAhIYFVq1ZJErgLTk5OTJw4kenTp/PXX38REhIicx0IYSdsmQjOADlr4SpYnsviBdQH1iulTgLBwG+5VRhrradorYO01kGlSpXKdWfnz5+nQ4cOREdH88cffxAQEGClj2FfBg8ezKpVq7h48SLBwcGEh4ebHZIQwsZsmQjCgRpKqapKKRegH/Bb1ota61ittZ/WuorWugoQBjyitb7rCYljY2Pp1KkTZ8+eZfny5bfsvSfyJjQ0lM2bN+Ph4UGbNm1YvHix2SEJIWzIZolAa50OjAT+AA4A87XW+5RS7ymlHrHWflJSUujRowf79+9n0aJFhISEWGvTdq127dqEhYVRr149evbsyX//+1+zQ7KJX/v8yq99ch//RQh7oQpb062goCAdEWFcNGRmZtKvXz9++eUXfvrpJwYOHGhydEVPQkICAwYM4LfffmPUqFF8+umnRW6APiHsxC3H1y+0v2itNa+88gq//PILn376qSQBGylWrBgLFy5k5MiRTJgwgT59+pCUlGR2WFYzY+cMZuycYXYYQpiq0DYW/+KLL/j666959dVXee21OzY8EvfB0dGRr776iqpVq/L6669z7tw5lixZUiSal2YlgUEBg0yNQwgzFcpE8Pvvv/Paa6/x2GOP8dlnn5k6o5jWmrikdK6lpBGfkk58sjGBRDFXJzxdnfByc8Lb3bnQznqWRSnFqFGjqFSpEgMHDiQ0NJQ//viD8uXLmx2aEOI+FbpEkJSURL9+/WjcuDE//vhjvpZXp2VksjvqKttOXOHIxWscuxjPsUsJxKfcfirIYi6OVCvlSbVSxXiglCdBVUrSuJIPbs6O+RS59fTu3ZuSJUvy6KOP0qpVK1auXEmNGjXMDksIcR8KXWWxq6ur9vX1JTw8PF/ORi9dS2HZnnP8dfgSW0/EZB/0/Yu7Ub20cWCvWNKD4m7OeLo5UczVyK0JKenEp6RzLTmd0zGJHLsUz/FLCZy5apSvuzo5EFi5BK1q+NG9YTkqlvSw+WexpoiICLp06YKDgwMrVqwotMN4tJ3RFoD1g9abGocQ+eCWxRKFLhE4ODjobdu2ERRk1YFKb5CclsHqAxdYuP0Mfx6+REampoqvBy0f8KPlA34EV/OlZLF7m88gNimNbSdi2HIsmi3HozlwLg6AwMoleDSgHN0blqPEPW47vx08eJCHHnqI2NhYlixZQmhoqNkh3TVJBMKOFJ1EUL16dX3s2DGbbDs2MY1ZYSf5YdNJohNS8S/uRo/G5XmsSXlqlvGyyT5PxyTy266zLN55hsMX4nFzdqBXkwo806oq1Ut52mSf1nT69GkeeughTp48ybx583jkEat1EckXiWmJAHg4F64rMiHuQdFJBDn7EVjLpWspTN1wnJ/DIklIzaBdrVI806oqIdX9cHTIv0re/WfjmLn5JIt2nCE1I5OOdUrzYvsaNKrok28x3IvLly/TtWtXtm/fzvfff8/TTz9tdkhCiH+SRJCb5LQMpm86waR1x0hMTad7o3IMD61O3XLFrbL9e3XpWgqzwiKZteUkVxLTeLhBWV7vVIuqfsVMjet2rl27Rs+ePVmzZg0TJkzg1VdfNTukPJkUPgmAEU1HmByJEDYniSAnrTXL9pxn/PIDRF1JomOdMozpWptqBawoJj4lnSl/HWfahuOkpGfSv1lFXnuwVoGtQ0hJSWHAgAEsXLiQd955h3feeafAN5uVOgJhR275Yyx0zUfv1/nYZEYv3M26Q5eo7e/Fz0Ob0/IBP7PDypWnqxOjHqzJE8GV+HrNUWZvO8WyPef5d9c6PNakfIE7yLq6ujJv3jyeffZZ3n33Xa5evcqECRNkSAohCji7SQRaa375O4r3l+4nLSOTt7vV5emQKvlaB3CvSnu58X6P+gxoXokxi/bw2i+7WLA9ig961C9wVzFOTk58//33eHt78+WXXxIbG8vUqVPtZsYzIQoju/h1Rsen8Povu1h36BLNqpbkk14NqVKAy9tvpU7Z4ix4LoTZ207x8YqDdPlyA/9+uA5PBlcuUFcHDg4OTJw4kRIlSjBu3Dji4uKYPXs2rq6uZocmhMhFkb9m33Yihq5fbWDTsWjGda/L3GeDC2USyOLgoHgiuDJrRrUhuJovby/exzMzwrl0LcXs0G6glOKdd95h4sSJLFy4kO7du5OQkGB2WEKIXBTZyuLMTM3kv47z2cpDVCrpwTcDmpjeGsjatNb8uCWSD5cdwMvVic8eb0S72qXNDusffvjhB4YOHUpwcDC///47Pj4+ZockhD0qesNQ3058SjrDZv3NxysO0rm+P7+NbFnkkgAYZ91Ph1Rh6YutKOXlyuAZ4UxYdZjMzIKV3AcPHsz8+fMJDw+nbdu2XLhwweyQhBA5FLkrgjNXkxgyI5wjF+N56+E6PB1SpUCVn9tKcloGY/+3l1//jqJtrVJ80TcAH4+C1cz0jz/+oGfPnlSoUIHVq1dTqVIls0Pis82fAfB6yOsmRyKEzdnHFcGOU1d49L+bOHMliR8GNWVQy6p2kQQA3Jwd+bR3Qz7oUZ9NRy/T/b8bs8cxKig6derEqlWruHjxIq1ateLQoUNmh8TSw0tZenip2WEIYaoikwiW7zlHvylheLg4snBECKE1S5kdUr5TyqhInj+8BWnpmt7fbmbdwYtmh3WDli1bsn79epKTk2ndujU7d+40OyQh7F6RSARztp3ihdnbqV/em/+90JIaNhogrrBoXKkEi0e2pGqpYgyZGc7MzSfNDukGAQEBbNiwATc3N9q2bcumTZvMDkkIu1boE8G3648xeuEeQmuW4qchze95eOiipkxxN+YPb0H72mV457d9jPttHxkFqBK5Vq1abNy4kdKlS/PQQw+xcuVKs0MSwm4V2kSgtWb8sgN8vOIgjwaUY+pTQbi7FL4Zv2zJw8WJyU8GMrRVVWZsPsmLc7aTkp5hdljZKlWqxIYNG6hRowbdunVjwYIF+R6Du7M77s7u+b5fIQqSQtlqKDw8nHd+28ePWyJ5qkVlxnWvh0MhGCrCTNM2HOeD3w8QUt2XKU8F4elacDqVX716lYcffpiwsDCmTZvG4MGDzQ5JiKKoaLUaykoCw0Kr8e4jkgTyYmjranz+eCO2noih/5QwouMLTk9kHx8fVq5cSceOHXnmmWf48ssvzQ5JCLtS6BLB2atJ/LglkmdbV2V0l9p20zzUGnoFVmDKk4EcvnCNxydv4XxsstkhZStWrBi//fYbvXr14pVXXuHdd98lP65W3//zfd7/832b70eIgqzQJYLohFSebV2VMV3rSBK4Bx3qlOGnoc25GJdC3ylbOHM1yeyQsrm6ujJ37lwGDRrEuHHjGDVqFJmZmTbd55oTa1hzYo1N9yFEQVfoEoGfp6skgfvUtEpJZg1pRkxCKn2+28LpmESzQ8qWNYz1yy+/zBdffMHQoUNJT083OywhirRClwjKertJErCCxpVKMHtoMAmp6fSZvIUTlwvOyKBZw1iPGzeOH374gX79+pGSUnDqNIQoagpdIhDW06CCN7OHBhvTYE4JIzK64CSDnMNYL1iwgEceeUSGsRbCRiQR2Lm65Yoz+9nmpKRnMGDqVqKuFJxiIoBXXnmF6dOns3r1ah566CFiYmKsun1fD198PXytuk0hCptC2Y/gfievF/+090wsA6aG4e3hzPzhLSjrXbA6WS1YsICBAwdStWpVVqxYQeXKlc0OSYjCpmj1IxDWV7+8Nz8Oac6VhDQGTN3KxbiC07QUoFevXqxcuZLz58/TokULGaxOCCuSRCCyBVT0YeYzTbkQl8xT07cRm5hmdkg3CA0NZePGjTg5OREaGsqqVavue5ujV49m9OrRVohOiMJLEoG4QWDlkkx5Mohjl+IZMjOcpNSCMzYRQL169diyZQtVqlSha9euzJo16762tyVqC1uitlgpOiEKJ0kE4h9a1fDji76N+fvUFUbO3k5ahm07dd2t8uXLs2HDBkJDQ3nqqaf46KOP8qUXshBFlSQCkauHG5bl/Ufrs+bgRf5vwe4CNw+yt7c3y5cvZ8CAAYwePZoXX3yRjIyCdfUiRGFh0yEolVKdgS8BR2Ca1vqjm15/DngByADigWFa6/22jEnk3RPBlYlJSGXCqsOU9HDh3w8XrB7dLi4uzJo1i/Lly/Ppp58SGRnJ7Nmz8fKy74mJhLhbNrsiUEo5At8AXYC6QH+lVN2bVputtW6gtQ4APgEm2CoecW9ebP8Ag0KqMG3jCb7985jZ4fyDg4MDn3zyCd988w3Lly+nVatWREZG5vn9FYpXoELxCjaMUIiCz5ZXBM2Ao1rr4wBKqbnAo0D2Gb/WOufs6sWAglX+IFBK8Xa3usQkpPLJikOU9HChX7NKZof1DyNGjKBGjRo8/vjjNGvWjMWLFxMcHHzH9/302E/5EJ0QBZst6wjKA6dzPI6yPHcDpdQLSqljGFcEL+W2IaXUMKVUhFIq4tKlSzYJVtyag4Pis8cbEVqzFGMW7WHF3nNmh5SrBx98kLCwMLy8vGjbti2zZ882OyQhCgXTK4u11t9orasD/weMvcU6U7TWQVrroFKlSuVvgAIAFycHvnuiCY0q+vDS3J1EnLTuUA/WUrt2bbZu3Urz5s0ZOHAgb7/99m2Hsn5lxSu8suKV/AtQiALIlongDFAxx+MKluduZS7Qw4bxiPvk4eLE9KebUt7HnWd/jChQI5bm5Ovry6pVq3jmmWd4//336devH4mJuY+htPP8Tnae35m/AQpRwNgyEYQDNZRSVZVSLkA/4LecKyilauR4+DBwxIbxCCsoUcyFHwY1BWDwD9uISUg1OaLcubi4MG3aND799FN+/fVXQkNDOXXqlNlhCVEg2SwRaK3TgZHAH8ABYL7Wep9S6j2l1COW1UYqpfYppXYCo4CnbRWPsJ4qfsWY9nQQZ2OTGfZjBMlpBbP9vlKK119/ncWLF3PkyBGaNGnC6tWrzQ5LiALHpnUEWutlWuuaWuvqWusPLc+9rbX+zXL/Za11Pa11gNa6ndZ6ny3jEdYTWLkkE/sEEBF5hdd/2VXgOpzl1L17d8LDw/H396dTp07SE1mIm5heWSwKr4cblmV0l9os3X2OT1ceMjuc26pZsyZhYWE8/vjjjB49ml69ehEXF0dN35rU9K1pdnhCmMqmPYtF0TcstBqnYhL5dv0xKpbwYEDzgtfHIIunpydz5swhODiY119/naZNm7Jo0SLq1r25n6MQ9kWuCMR9UUrx7iP1aFurFG8t3sv6QxfNDum2lFK88sorrF27ltjYWJo1a8b8+fPNDksIU+UpESilFiqlHlZKSeIQ/+Dk6MB/BzShVhkvXvh5O/vPxt35TSYLDQ1l+/bteFb0pG/fvrz00kskJxesyXiEyC95PbBPAgYAR5RSHymlatkwJlEIebo6MX1QU7zcnBk6M5yL1wr+QbVcuXLUfK0m5R8qz9dff02zZs3Yv1/GPBT2J0+JQGu9Wms9EGgCnARWK6U2K6UGK6WcbRmgKDz8vd2Y9nQQVxLTGD7r7wLbrDQnBycHHuj/AMuWLeP8+fMEBgby3XffSasiYVfyXNSjlPIFBgFDgR0Yw0s3Ae5/vkBRZNQv783Evo3Yceoqby7YXWgOqF26dGH37t2Ehoby/PPP06tXL6Kjo80OS4h8kdc6gkXABsAD6K61fkRrPU9r/SLgacsAReHTuX5Z/tWpFv/beZZv1h01O5w88/f3Z/ny5Xz++ecsXbqURo0asX79erPDEsLm8npFMFVrXVdrPV5rfQ5AKeUKoLUOsll0otAa0bY6PQLK8dnKwyzfUzBHKwUI8A8gwD8g+7GDgwOjRo0iLCyMYsWK0b59e8aMGUNKSop5QQphYyovl+5Kqe1a6yZ3ei4/BAUF6YiIiPzerbgHyWkZ9J8axoFzcfz6XAj1y3ubHdJdiY+P56WXXuKHH36gfv36/PDDDwQFyXmPKLRuOb3gba8IlFL+SqlAwF0p1Vgp1cRya4tRTCTELbk5OzLlySBKergwdGYEF+IKfkuinDw9PZk+fTpLliwhJiaG5s2b8+abb0ozU1Hk3KloqBPwGcYQ0hOAzy23UcAY24YmioJSXq58P6gpcclpPPtjBEmpBasl0RMLn+CJhU/cdp1u3bqxb98+Bg8ezMcff0zjxo3ZsmVLPkUohO3dNhForWdqrdsBgyyDwmXdHtFaL8ynGEUhV6dscb7s15g9Z2J5/dddBaolUVRcFFFxUXdcz8fHh2nTpvHHH3+QmJhIy5Ytee211245z4EQhcmdioayTpWqKKVG3XzLh/hEEfFg3TK82bk2v+8+x5drCu+0Ew899BB79uxh+PDhTJgwgUaNGrFu3TqzwxLivtypaKiYZekJeOVyEyLPhoVWo3dgBb5YfYQlu86aHc49K168ON9++y1r1qwhMzOT9u3b069fP86cud0EfEIUXHlqNVSQSKuhwi0lPYMnpm1ld1Qs84a3IKCij6nxtJ3RFoD1g9bf0/uTkpL45JNP+Oijj3B0dOTtt9/mlVdewcXFxXpBCmEd99ZqKPvdSn2ilCqulHJWSq1RSl3KUWwkRJ65Ojny3ROBlPJy5dkfIzgXm2RqPC0qtKBFhRb3/H53d3feeecd9u/fT4cOHfi///s/GjVqxJo1a6wYpRC2ldcOZQ9preOAbhhjDT0A/MtWQYmizdfTle+fbkpSagZDZ0aQmJpuWizjO45nfMfx972dqlWrsnjxYpYuXUpqaiodO3akT58+nD592gpRCmFbeU0EWRPYPAz8orWOtVE8wk7U8vfi6/6NOXAujlHzCvZUl3fj4YcfZt++fbz33nssWbKEWrVqMWbMGGJj5ScjCq68JoKlSqmDQCCwRilVCpBeNeK+tKtdmn8/XJcV+87z+SpzprrsNb8Xveb3suo23dzceOuttzhw4AA9e/Zk/PjxVKtWjYkTJ8pQFaJAyusw1G8CIUCQ1joNSAAetWVgwj4807IK/ZtV5Jt1x1i0487t+a0tOjGa6ETbjDJapUoVfv75Z7Zv305QUBCjRo2iVq1a/PTTT2RmZtpkn0Lci7uZcaw20Fcp9RTQG3jINiEJe2JMdVmf4Gol+b9f9/B3ZIzZIVld48aN+eOPP1i1ahW+vr48+eSTNGnShBUrVhSoznXCfuW11dAsjKEmWgFNLTcZfUtYhYuTA98ODKScjxvDfvybqCtFs7dux44dCQ8PZ86cOVy7do0uXbrQunVrSQjCdHm9IggCWmqtR2itX7TcXrJlYMK+lCjmwrSnm5KakcmQGRHEp5jXksiWHBwc6NevHwcOHOCbb77h1KlTdOnShWbNmrF48WIpMhKmyGsi2Av42zIQIR4o7cmkgU04eimel+fsICMfWhJ1qNqBDlU72Hw/N3NxcWHEiBEcPXqUadOmERMTQ48ePQgICGDevHlkZBSswflE0ZbX+QjWAQHANiC72YPW+hGbRXYL0rO46Ju15SRvLd7HsNBqjOlax+xw8kV6ejpz587lww8/5ODBg9SsWZPRo0fTv39/XF1dzQ5PFA237Fmc10TQJrfntdZ/3kdQ90QSgX14Z/FeZm6J5ONeDejbtJLZ4eSbzMxMFi5cyAcffMCuXbsoU6YMzz//PM899xxlypQxOzxRuN3fEBOWA/5JwNlyPxzYbpXQhMjFW93q0rqGH2P/t5ew47abRL7Lz13o8nMXm23/bjk4ONC7d2927NjBypUrCQoKYty4cVSqVImnn36a7dvlZyesL6+thp4FfgUmW54qD/zPRjEJgZOjA/8d0IRKJT14/qe/iYxOsMl+ktKSSEozd7yj3CilePDBB1m6dCmHDh1i2LBhLFiwgMDAQEJDQ1m4cCHp6UWzQl3kv7xWFr8AtATiALTWR4DStgpKCABvd2e+f7opGhgyM4K45DSzQzJFzZo1+frrr4mKiuLzzz/n9OnT9OrVi8qVKzN27FiOHz9udoiikMtrIkjRWqdmPVBKOQHS8FnYXBW/Ynw7MJCTlxMYOXsH6Rn227zSx8eHUaNGcfToURYtWkTjxo0ZP3481atXp0OHDsydO1fmUxb3JK+J4E+l1BiMSewfBH4BltguLCGua1Hdlw971uevw5f44PcDZodjOkdHR3r06MHSpUuJjIzk/fff5/jx4/Tv359y5crx8ssvs2vXLrPDFIVIXhPBm8AlYA8wHFgGjLVVUELcrG/TSgxtVZUZm08yKyzSatvtVrMb3Wp2s9r28luFChUYO3Ysx44dY/Xq1XTq1InvvvuOgIAA6tevz3/+8x8pOhJ3lOcZyiwjjqK1vmTTiO5Amo/ar4xMzbM/RvDn4UvMHNyMVjX8zA6pQIqOjmb+/PnMnj2bjRs3AhAcHEz//v3p06cP/v7SN9RO3Vs/AqWUAt4BRnL96iED+Fpr/Z41I8wrSQT27VpyGr2/3cK52CQWvdCS6qU8zQ6pQIuMjGTevHnMnj2bXbt24eDgQPv27enduzePPvqoJAX7cs+JYBTQBRimtT5hea4a8C2wQms90cqB3pEkAnE6JpEe32zC082Jhc+H4Ot57z1v73fO4sJk//79zJkzhzlz5nDs2DGUUgQHB9OzZ0969OhBjRo1zA5R2NY9dyh7EuiflQQAtNbHgSeAp+64V6U6K6UOKaWOKqXezOX1UUqp/Uqp3Za5kCvfaZtCVCzpwZSngjgfm8yzP0aQnCbj8uRF3bp1ef/99zly5Ai7d+/m3XffJSUlhTfeeIOaNWtSv359xo4dS3h4uAx+Z2fulAictdaXb37SUk/gfLs3KqUcgW8wrijqAv2VUnVvWm0HxmQ3DTE6rH2S18CFfQusXIIv+gaw4/RVXp23s8hMdZkflFI0aNCAt956i7///puTJ0/y5ZdfUrp0aT766COaNWtGmTJleOKJJ/jpp5+4ePGi2SELG7tTIki9x9cAmgFHtdbHLX0Q5nLTrGZa63Va66zB58OACnfYphDZujQoy7+71mH53vOMXy7NSu9V5cqVeemll1i7di0XLlzg559/pnPnzqxcuZInn3ySMmXKEBQUxNixY9m4cSNpafbZsa8oc7rD642UUnG5PK8Atzu8tzxwOsfjKKD5bdYfAizP7QWl1DBgGEClSvYzAJm4syGtqnI6JpGpG05QsaQHT7WoYnZIhZqvry8DBgxgwIABZGZmsmPHDlasWMGKFSv46KOP+PDDDylWrBitWrWiXbt2tG3blsDAQJyc7nQoEQVZnpuP3vWGleoNdNZaD7U8fhJorrUemcu6T2C0TGqjtb7t7N5SWSxulpGpGT7rb9YevMCUJ4PoWDfvo3ROCp8EwIimI2wVXpFx9epV1qxZw7p161i3bh379+8HwMvLi9atW9O2bVvatm1LQEAAzs63LTkW5ri/YajvaY9KtQDGaa07WR6PBtBaj79pvY7A1xhJ4I6FkZIIRG4SU9PpNyWMIxfimTc8mIYVfMwOqci7cOEC69evZ/369axbt45Dhw4B4O7uTrNmzWjZsiUhISG0aNGCkiVLmhytwKRE4AQcBjoAZzCGrh6gtd6XY53GGJXEnS0D2d2RJAJxK5eupdBz0iaS0zJZNCKEiiU97viexDSjisrD+c7rits7e/YsGzZsYPPmzWzevJkdO3Zkz7RWp04dQkJCaN68OYGBgdSvXx8XFxeTI7Y7+Z8IAJRSXYEvAEdgutb6Q6XUe0CE1vo3pdRqoAFwzvKWU3ea9UwSgbidoxev8dikzZQu7savz7XAx+P2Bxt76keQ3xISEggPD2fz5s1s2rSJLVu2cOXKFcCYqrNRo0YEBQURFBREYGAgdevWlSIl2zInEdiCJAJxJ2HHo3nq+200qODNT0Oa4+7ieMt1JRHkH601x48fJyIigoiICP7++2/+/vtv4uKM9ihubm4EBARkJ4fGjRtTu3ZtuXKwHkkEwr4s23OOF2Zvp12t0kx+MhBnx9xbSksiMFdmZiZHjx7NTgwRERFs376d+Ph4AJycnKhTpw4NGzakUaNGNGzYkIYNG+Lv748xAo64C5IIhP35KSySsf/bS68mFfjs8Ya5HjgkERQ8GRkZHD58mF27drF79252797Nrl27iIqKyl7Hz88vOyk0bNiQBg0aULduXTw8pK7nNm6ZCKTxryiyngiuTHR8KhNXH8bP04XRXeuYHZLIA0dHR+rUqUOdOnXo169f9vMxMTHs2bMnOzns3r2byZMnk5RkTDWqlKJ69eo0aNCA+vXrU79+fRo0aECNGjWkn8MdyLcjirSXOjxAdEIKk/86jq+nC8NCq9/w+qCAQeYEJu5ayZIladOmDW3atMl+LiMjg+PHj7N792727t3L3r172bNnD4sXL84eL8nFxYXatWv/I0FUqlRJipcspGhIFHkZmZqX5u7g993n+PzxRvQKlJFMirrk5GQOHjzInj17bkgQp09fH+zAy8uLevXq/SNBlCpVysTIbUrqCIR9S0nP4JkZ4YQdj2HqU4G0r230Pr6caIyp6Ochk9zYg9jYWPbt25edGLKW0dHR2euULl06OylkJYh69erh5eVlYuRWIYlAiPiUdPpPCePIxWv8NKQ5QVVKSmWxQGvNhQsXbrhy2Lt3L/v27SMhISF7vSpVqvwjQRSy5q2SCIQAuByfwuPfbeFyfAqzhwbz4mpjQFxJBOJmmZmZnDx58h8J4uDBg6SnpwPg7OxMgwYNsjvFBQUFFeRe05IIhMhy5moSfb7bQkJqOh6lp+LidkkSgciz1NRUjhw5wp49e9i5c2d2/4erV68CRuV0w4YNCQwMJDg4mJCQEGrUqFEQKqYlEQiR06noRPpM3sLFhGj8K81l8/CFZockCjGtNSdOnLih13RERER2r2lfX19atGhBSEgIISEhNG3a1Iw+D5IIhLjZ8UvxdPxiGUpp1rzSjSp+xcwOSRQhmZmZHDx4MHsQvs2bN2eP0Ork5ETjxo1p164d7dq1o1WrVnh6eto6JEkEQuTmy42/MGWVO95u7swb3iJPI5YKca+io6MJCwtj06ZNbNiwga1bt5KWloaTkxPNmjWjXbt2tG/fnhYtWuDu7m7t3UsiEOJW9p2Npf+UMLw9nJk/vAVlva3+AxQiVwkJCWzevJm1a9eybt06IiIiyMjIwNXVldDQULp27UrXrl2pWbOmNXYniUCI3JyONToYxcR5MXDaVkp7uTJ3eDClve40E6sQ1hcXF8eGDRtYu3Yty5cv58ABYy7uBx54IDsptGnTBje3e/r/lEQgRG5y9iOIOBnDU9O3Ud7HndnPBlPKy9Xc4ITdO3HiBMuXL2fZsmWsXbuWpKQk3N3d6dixI7169eKRRx6hRIkSed3cLRNB7mPzCmGHgqqU5PunmxJ1JYm+U7ZwPjbZ7JCEnatatSojRoxg6dKlREdHs3z5coYMGcLOnTsZNGgQpUuXpnPnzkydOpVLly7d834kEQiRQ4vqvvw4pBkX41LoM3kLp2MSzQ5JCMCYC7pz5858/fXXREZGsnXrVkaNGsWRI0cYNmwY/v7+tG/fnkmTJnHx4h2nf7+BJAIhbtK0Skl+Htqc2KQ0+kzewvFL8WaHJMQNlFI0a9aMjz/+mKNHj7Jjxw5Gjx7NuXPneOGFFyhXrhzdunVj3rx52cN0344kAiFy0aiiD3OHBZOankmfyWEcOn/N7JCEyJVSioCAAD744AMOHDjAnj17eP3119m5cyf9+vXD39+foUOH3n4bUlks7NmSQ0sA6F6re66vH70Yz8BpYaSmZzJrSHPql/fOz/CEuGcZGRmsX7+eWbNmsWDBAq5duyathoS4V5HRCQyYupW45DRmDG5GYOU8t9IQokBISEigWLFi0mpIiNwcunyIQ5cP3Xadyr7F+OW5Fvh5uvLk91vZfPRyPkUnhHUUK3b74VMkEQi7NnzpcIYvHX7H9cr5uDNveDAVSrgz6Idwft99Lh+iEyJ/SCIQIo9Ke7kxf3gLGlX0ZuSc7czYdMLskISwCkkEQtwFHw8XZg1pTsc6ZRi3ZD+frDhIYatnE+JmTmYHIERh4+bsyLcDm/DW4n1MWn+Mc7HJfNSrAa5OjmaHZshIh7QESE+FjBRIT4GM1OvLjFTQmdfXvzmRKQdwdAFHZ8vSBZxcrt939gBndzB/ohVhJZIIhLgHTo4O/Kdnfcp5u/H5qsOcuZLE5CcDKVHMSlMUpiZA/EVIioHEK5ZlNCTGGPeTrkByHKTGQ0o8pF6zLOMhPR+GxlCO4OIJrp43Lb3Awxc8ShpL95I3PvYsbawjChRpPirs2urjqwHoWK3jPW9j8c4z/OvX3ZT3cWf6oKZUvd0ENxlpEHcGrp6G2Ci4ds444MefN5bXzkP8BeOAnisFbt7GgdW1uHFQdfW66WBcHFw8LGfyrjnO6l3B0dU403e4+eolx9m9zjDizEgzrigy0q5fSaSnGEkqNd5Y3pyEkuOuJ63M9Nw/gnMx8PIHr7LgVcaytDz2qQw+FaFYaXCQkmsrk34EQthSxMkYhs36m8zMTCY9XIoQ72i4egpiT18/6MeeNg78OYtlAFy8jDNlL3/wLGNZljYOhlln0+4ljaWbdy4H8QJIa0i5ZiSEpBjjSiYx2khy184b38O1C5bleUi/aRgER1cjIfhUun4rURV8q0PJ6kbSE3dLEoEQudl5ficAAf4BeX9TagLEHIfoYxBzDKKPQ8wxIi9eZUjsYE7osrzlNIunHVeiHJ2geHnjQOZdEbwrGAc4b8uteFlwsfMpMrWGlDiIPWNJnKduvMWehoSbRtb09LckhWrG0rcGlK4DJaoUjkRpDkkEQuQm53wEN8jMgCsn4dIhuHz4hgM+127qQ+BZxjhL9a3GNa8HeHVfdVafcaRPo5K83ysIVxfn/PgoRVtqAsScsPwdjl7/W0QfvTFJOLlDqVpGUihdB0rXhVK1jQQsldu3/AKksljYNSedSfm0FNj3P8tB/5BlecQoH8/i4WeceVZrB77VLAd+yxlpjspPL2BKW83E1Yf5eu1RjlzZxqSBTWT6y/vlUgz86xu3myXHGn+viwfg0kG4uB+Or4ddc66v41rcSAhZyaFsQ/BvKEVMFnJFIOxDZiZcPQnn98KFvcbB4tIh0i8fvvFsyKeycUZZqpZx4PCrBX41wN3nrne5bM85/vXLLlydHfmqX2Na1fCz0ocReZIYY0kMB3Lc9hktrgBQxt+2bCMoGwDlAsC/gVEPUzRJ0ZCwI2lJxoH+/B7jwH9+D1zYZ7RuAaOdfMlqUKo2s86FE+nsxtjHfjQOClYurz96MZ7nf/qbo5fiebVjTUa2ewAHB7svojCP1kbl9LldcG6nsTy7E66dvb5OyepGcigXYFk2AbfiJgVsVZIIRBF17YLlQL/n+oE/+sj1ljkuXkZxQpn6xtmef32jaMDZKKq5ZR2BFSWkpDNm0R4W7zxLaM1SfP54I5kPuaCJv3g9OZzdadyPPW15URlXhxWbQoWmUKEZ+NUsjM1bzUkESqnOwJeAIzBNa/3RTa+HAl8ADYF+Wutf77RNSQR2KiPdOMCf33P9dmHvjRWF3pUs5cgNrh/4fSrf9ge7+fRmAEIqhtg0fK01P289xftL9+Pl5sTnfQJoU7OUTfcp7lNCNJzbAVF/Q1S4cUu+arzmWhzKB0LFZkZyKB9oNO8t2PI/ESilHIHDwINAFBAO9Nda78+xThWgOPA68JskAgEYlX9ZZfnndxv3Lx64Xnnr6GJU+pVpcP0sv0w9cC/48wQcOn+NF+ds5/CFeJ5tXZV/daqNi1OhO7O0T5mZRkulqHA4vQ2iIow6h6yrT98aULE5VG4BlVoYxY8Fq6WSKYmgBTBOa93J8ng0gNZ6fC7rzgCWSiKwM1rD1cibyvL3GG3Hs3j45TjLtxz4/WoYvWOtIL+uCHJKTsvgw98PMCssktr+XkzoE0DdckWiDNr+pMTD2e2W5BAOp7caHejAaFZcqQVUDjGWZeqZ3cfBlETQG+istR5qefwk0FxrPTKXdWdwm0SglBoGDAOoVKlSYGRkpE1iFjaUlmSc1V/Ye/3Af2Gv0ZEIAAW+D1w/w/dvaBTvePnb9KwqP+oIbmXNgQu8uXAPVxNTeblDDZ5rUx0nR7k6KNQyM41+J5Gb4NQWiNwCcVHGa67eRlFS5RZQKQTKNzGG/cg/hbsfgdZ6CjAFjCsCk8MRdxJ/8XqRTlZZ/uUjxhg2YIyHU6YeNOxjKctvaBT1uHiYG3c+61CnDCtfKcHbv+3js5WHWbX/Ah/3bkhtf7k6KLQcHKB0bePWdIjx3NVTRkI4tdlYHl1lPO/kZiSGKqFQtbXROsnJSoMW3iVbJoIzQMUcjytYnhNFRUa60bMzZ1n++T2QcPH6OsUrGGf5dbpfr8QtUbUwtriwiRLFXPi6f2M61/PnrcV76fbVRoa3qcaL7Wvg5ixDJRQJWWMlNeprPE6ItlwtbIITG2DdB7AOY3jvSsFQpTVUbWM0XXXMn3N1W+4lHKihlKqKkQD6AQNsuD9hS8mxRlv885aD/gVLBW7WkMcOzsZZUI0Hr7fYKVOvMLSkKBAebliWkOq+fLjsAN+sO8bvu8/xn54NCHlAOqEVOcV8oU434wZGYojcaCSFkxtgzbvG8y5eRv1C1dZGcvBvYLM6Bls3H+2K0TzUEZiutf5QKfUeEKG1/k0p1RRYBJQAkoHzWut6t9umVBbbmNbGpWx2Wb7ldjVHvYx7SUtZfoPrZ/l+NU27rL0fZtYR3Mqmo5cZs2gPkdGJPNywLGO61qG8jwxRYTfiLxoJISsxRB81nnfzhsqtjMRQNdToD3N39WfSoUzkIi3Z6IKfVY6fVYmbEmtZQRnj6WR3xrLcvMoWtGZx9+yeRh/NB8lpGUz+8ziT1h9FKRjR9gGGhVaT4iJ7FHcWTm6EE38aySHrpKxYaajWFqq3M5bFy91pS5II7F78pRt7357fY7RuyKrAdS5mFOVk98JtCGXqyhDJJou6ksh/lh1g2Z7zlPdxZ9SDNenRuDyOMkyF/bp6Ck78BcfWGYPrJV42nverdT0xVGmV20xwkgjsRmaGMU5+Vjl+1oE//vz1dYqXv3HIBf+GdluBa40ZyvLD5qOX+c/yA+w9E0dtfy/e7FKbNjVLoYrIlZm4R5mZRqe2rKQQudmY5MfBCcoHXb9aKB8Ijs6SCIqk1ARLBe5uOGc58F/Yf322JwcnY4yUnEMu+DeQCtwcCmIdwa1kZmqW7D7LZysPcTomiWZVSvJyxxqEVPeVhCAMackQtc1ICsfWwdkdgDYqnsdESSIo9LLb5u8xDvrn91gqkSx/Pzeff1bglqqV3x1WCp3ClAiypKZnMmfbKSatP8qFuBSCKpfg5Y41aPWAnyQEcaPEGKPC+dg66P6FJIJCIzMTrpy4fpaf1WonZ9GOdyXLxBpZB/6GMgPTPSqMiSBLcloG8yNOM2ndMc7HJVO/fHGebV2Nrg3K4iw9lMU/Fe6exUVWesr1cfOzDvoX9kJqvPG6cjSKdqq3u37A969fKAZXE7bn5uzIUy2q0LdpRRZuP8PUDcd5ee5OPllxiEEhVXg8qAI+HoWvSa/If3JFkF+Srlw/u8866F8+BJnpxusunjcW7fg3NJKAs5u5cRdxhfmK4GaZmZq1By8y5a/jbDsZg6uTA90aluOJ4EoEVPSRYiMhRUP5RmtjQoubD/qxOUbU9PT/Z9GOnbbaMduhy4cAqOVXy+RIrGv/2Th+3hrJ/3acISE1g9r+XjzWpDyPBpSnTHE5ubBTkghsImuylHO7LAd9yzJr8oqsETWzD/qWpWdpM6MWdiQ+JZ3/7TjDgu1R7Dh1FQcFLR/wo3ujcjxUt4wUHdkXSQT3LSPd6IWbczq783uuN9V0crN0yMo6y28kHbIKgSWHlgDQvVZ3kyOxveOX4lm04wyLdpwh6koSjg6KkOq+dKrnT/vapSknw1iYSmtNWoYmI1OTlplJeoYmPSOT9EzjOUcHlX1zyl464OigcHZUeSn6k0RwVzLSjAHVck5ufWHv9QHWXDyNs/ucE1z71si3kQKF9RSlOoK80lqz50wsy/eeZ8Xe85y4nABAzTKetKlZitCapQisXAIPF/l/ziutNQmpGcTEp3I1KZUriWlcTUzlamKacUtKJTYxjfiUdBJTM0hITSchJZ2ElAwSU41lakbmPe/fQUExFyeKuTpRzNURT9es+054uTnh5+nKmK51JBHcUnqq0XIn56TVF/ZdnxbRxevGA37ZAGP8HXNnGhJWYo+JICetNUcvxrP+0CXWH77IthMxpGVonBwUDSp407yqL0GVS9CwojelveyvbiE9I5PohFQuXUvh0rUULl5LznHfWF6KT+FiXApJaRm33I6XqxPF3Z3xcjMOzh4ujjccuD1cnHBzdsDJQeHkaFla7js7KhyUIlPr7KuDrFt6pnHVkJyWSXyKJbmkphOfkmFJNOlcS07ncnwKhz7oIs1HAaO55oV9OQ76O42euJlpxuuu3kZ5fvNhxgG/bIAx76hU4ooiSilFjTJe1CjjxbOh1UhISSf8ZAxbT8Sw7UQM3288znd/GieLZb3daFjBmzpli1PL8p4qvh6FblY1rTWxSWlcjk/h0rVUy/LmA3syl+NTiE5IJbdz5eJuTpTycqW0lxuNKvhQysuVUl6u+BZzwcfDhRIezvh4OOPj4YK3u7Pp/TrudMJfdBNBWtI/D/oXD1xvrunmY5zhtxhhOeg3kpY7wu4Vc3Wiba3StK1lNGhISs1g79lYdp2+yu6oWHZHXWXl/gvZB0cXRwcqlHCnsq8HlX2LUaGEO/7ebpQp7oZ/cTdKFnPBw8XRpk1X0zIyiUtKIy45ndikNGKT0riSYDnAx6dw2XKwz7pFx6eSnvnPA6Ozo6KUpyulirtRoYQHjSuVoLTlAJ91K+3lip+na6EbBfZO33/RSASpiUYZflbRTtZBP2tkTfcSxsE+5EVjWS4AfCpLT1wh7sDdxZGmVUrStMr18amSUjM4ejGeQxeuceTiNU5FJxIZnUj4ySvEp6T/YxuODgovNyeKuznfsPR0c8LJwSj2UErhoMAha+mgSM/QpKRnkJyWmb1MTssgJd1YxlkO+gmpty6ScXF0wM/TBV9P4yBet2xx/CwHcz9PF0p5uuJnOcB7uzvbbV+LwldHENhERyyadP0s/9wuozWPtlS0ePheP9hnLb0rykFf5Op07GkAKnpXvMOa4k6yilwuxKVwPi6ZC7HJXElMJS45jWvJ6cQlGctryenEJRsVp5mZmkwNGVqjtXE/Uxvl304OCjdnR1ydHLKXrjkeF3dzxts962aUwWc9LlHMBT9PV4q7OdntwT0XRaiyuJyTjhhmaZJZrPSNlbjlAowhluUPL4QQNytClcVeZaDfVOOgX4RmyhLmmLd3HgB96/c1ORIhzFMIE0FZqN3V7ChEEfFtxLeAJAJh36SJjBBC2DlJBEIIYeckEQghhJ2TRCCEEHau8FUWC2FFv/b51ewQhDCdJAJh1/w8/MwOQQjTSdGQsGszds5gxs4ZZochhKkkEQi7JolACEkEQghh9yQRCCGEnZNEIIQQdk4SgRBC2DlpPirs2rKBy8wOQQjTSSIQds3D2cPsEIQwXZEqGpq46nCu9wvUPtaNz/N+bPUZigprfD+TwicxKXySFaKxsTz+3xT4fQiryv4N3Offrkglgi/XHMn1foHax58f5Xk/tvoMRYU1vp/5++Yzf998K0RjY3n8vynw+xBWlf0buM+/XZFKBEIIIe6eTROBUqqzUuqQUuqoUurNXF53VUrNs7y+VSlVxZbxCCGE+CebJQKllCPwDdAFqAv0V0rVvWm1IcAVrfUDwETgY1vFI4QQIne2vCJoBhzVWh/XWqcCc4FHb1rnUWCm5f6vQAelZDZ6IYTIT0prbZsNK9Ub6Ky1Hmp5/CTQXGs9Msc6ey3rRFkeH7Osc/mmbQ0Dhlke1gIOWTFUP+DyHdcquuz984N8B/b++cE+voPLWuvOub1QKPoRaK2nAFNssW2lVITWOsgW2y4M7P3zg3wH9v75Qb4DWxYNnQEq5nhcwfJcrusopZwAbyDahjEJIYS4iS0TQThQQylVVSnlAvQDfrtpnd+Apy33ewNrta3KqoQQQuTKZkVDWut0pdRI4A/AEZiutd6nlHoPiNBa/wZ8D8xSSh0FYjCSRX6zSZFTIWLvnx/kO7D3zw92/h3YrLJYCCFE4SA9i4UQws5JIhBCCDtnd4lAKfW4UmqfUipTKXXL5mJ3Gh6jsFJKlVRKrVJKHbEsS9xivQyl1E7L7eZK/kJHhjvJ03cwSCl1KcfffagZcdqKUmq6Uuqipf9Sbq8rpdRXlu9nt1KqSX7HaBa7SwTAXuAx4K9brZDH4TEKqzeBNVrrGsAay+PcJGmtAyy3R/IvPOuT4U7u6n96Xo6/+7R8DdL2ZgC5dqiy6ALUsNyGAd/mQ0wFgt0lAq31Aa31nXom52V4jMIq57AeM4Ee5oWSb2S4k6L9P50nWuu/MFon3sqjwI/aEAb4KKXK5k905rK7RJBH5YHTOR5HWZ4rCsporc9Z7p8HytxiPTelVIRSKkwp1SN/QrOZvPw9s9fRWqcDsYBvvkSXP/L6P93LUizyq1KqYi6vF2VF+Xd/W4ViiIm7pZRaDfjn8tK/tdaL8zue/Ha7z5/zgdZaK6Vu1X64stb6jFKqGrBWKbVHa33M2rGKAmUJMEdrnaKUGo5xhdTe5JhEPiiSiUBr3fE+N5GX4TEKrNt9fqXUBaVUWa31Octl78VbbOOMZXlcKbUeaAwU1kRwN8OdRBXR4U7u+B1orXN+3mnAJ/kQV0FSqH/390OKhnKXl+ExCqucw3o8DfzjCkkpVUIp5Wq57we0BPbnW4TWJ8Od5OE7uKk8/BHgQD7GVxD8BjxlaT0UDMTmKEYt2rTWdnUDemKU/aUAF4A/LM+XA5blWK8rcBjjLPjfZsdtxc/vi9Fa6AiwGihpeT4ImGa5HwLsAXZZlkPMjtsKn/sff0/gPeARy3034BfgKLANqGZ2zCZ8B+OBfZa/+zqgttkxW/nzzwHOAWmWY8AQ4DngOcvrCqNl1THL/32Q2THn102GmBBCCDsnRUNCCGHnJBEIIYSdk0QghBB2ThKBEELYOUkEQghh5yQRCCGEnZNEIIQQdu7/AapDXHLWrI6PAAAAAElFTkSuQmCC\n", + "text/plain": [ + "
" + ] + }, + "metadata": { + "needs_background": "light" + }, + "output_type": "display_data" + } + ], + "source": [ + "# Density plot\n", + "test_distances = classify.get_distances(classifier, X_test)\n", + "classify.densityplot(test_distances, np.squeeze(y_test))" + ] + }, + { + "cell_type": "code", + "execution_count": 12, + "metadata": {}, + "outputs": [ + { + "data": { + "text/html": [ + "
\n", + "\n", + "\n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + "
DistanceLabel
00.512470Unhealthy
11.375290Unhealthy
2-0.881426Healthy
31.362346Healthy
4-0.781523Unhealthy
\n", + "
" + ], + "text/plain": [ + " Distance Label\n", + "0 0.512470 Unhealthy\n", + "1 1.375290 Unhealthy\n", + "2 -0.881426 Healthy\n", + "3 1.362346 Healthy\n", + "4 -0.781523 Unhealthy" + ] + }, + "execution_count": 12, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "X = pd.DataFrame(test_distances, columns=['Distance'])\n", + "y = pd.DataFrame(y_test, columns=['Label'])\n", + "X['Label'] = y\n", + "\n", + "X.head()" + ] + }, + { + "cell_type": "code", + "execution_count": 13, + "metadata": {}, + "outputs": [ + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAiUAAAB6CAYAAABgKiPgAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/Z1A+gAAAACXBIWXMAAAsTAAALEwEAmpwYAAAdA0lEQVR4nO3deXBc1aHn8e/pvaXWYlte5AWvbMY4JtgEm80kL8T1YogDnpohJEA9EirFEoYM4U0lhCUP3gxMBijgDxKGmGwVZgYm8YsDfsTBDpCEgBMcAgaDARmMjW3tarV6P/PHud2SwZsayX0l/T5Vt2737Vb30bF1+9fnnHuOsdYiIiIiUm2BahdAREREBBRKRERExCcUSkRERMQXFEpERETEFxRKRERExBdCg3nyihUr7Pr164erLCIyzIwxAOiqOxGpInOwBwbVUtLa2vrxiyIiIiJyAOq+EREREV9QKBERERFfUCgRERERX1AoEREREV9QKBERERFfUCgRERERXxjUPCUiIh+XtZaihXyxSKFoKRTdfT48dYqBgIGAMQQDhoAxhAKGQOCgUxyIyAinUCIiH5HNF+lJ5+hO50mm8/Skc/Rk3O1kxm293pbM5EllC6SyBfpyBdLZApl8kXS+QDZfJJsvkisUyRUsuUKRfNF+JGwEDGDADJhTyVqLBYrWUixCwboAEzAQCgYIBw2RYIBwMEAsHCQSChALBYhFgsTDQWojIWqiQRLREIloiLpYmEQsRF00RF3M3a+Ph6iPhWmsCRMPB8uTy4lIdSiUiIxi6VyBjlSW9t4snalc+fgDT79JRyrnHc/S2Zejuy9Hd1+enkyOXN5SGw1SGw1RGw0RDwepiQSJhYPEQgGi4QDRkLs/rjZCc0OwfCwackEh4gWHUoAIBQKEgv1BpBIDW1nyA0JOLl8kW3ABKFPeCqRzbr+vJ8POjhTpXJG+XIG+bIHebIFU1gtZ6TyFoqU+HqYx7kLK+NoI42oiNCUiTEhEGV/r9hNqIzQlokxIRAgH1QMuMpQUSkRGiHyhSGefCxJtySwdqSxtvVk6erO0JTO0Jl346Ei5rTOVw1qo81oD6mL9f+7b9yapjYaYWBdldlOtCx8R16pQEwkRCwd82WpgjCFoIBgIEh3is1c2Xyy3AJVahXrSOdp6c+xoT3n383T25ejytkQ05EJKXZRJdVEm18eYXB9lUl2MSfWl+zESQ11YkVFKfykiVZLJF+jozdHWm6G9N1sOG229GVp7srSnMrSVg0aOZDpPbSxIQzxcDhl1sRC1EdcVMXNCDQumNVAf6++eiIb2Dxe/9vb/cckx1fmlfSwSCjA+FGF8beSInl+0lp50nq6+HJ2prLfP8bf3uuhOt5YDYntvloAxTPSCy9SGOM2NMaY2xplSH2NKQ4zmhjgTaiMaLyNj3qBDSUtLCytXruSVV14pH7v11ltJJBLccMMNB/yZRx55hM2bN/PAAw9UXlLP5ZdfzsqVK1m9ejX33nsvV155JTU1NQAkEgmSyeTHfg+RwbLW0pstuFaLUutFb5b2XteC0Zr0gkcyS7vXkpHJFV3AiIe9IOHGPCSiIepjIZobGqjzHquPhUlEQ/rQ8pGAMTTEwzTEwxwzvuagz7PWksr2d6O5fY539vWWW7vakll6s3kmJqJMaXCBZVpjnGnj4jQ3xJnaGGNaY5yGeNiXLVgiQ2VEt5Tce++9fPnLXy6HEpGhUizacldJuxcuSmGjdUBXSWnrTOUIBCiHjLpof0tGIhqiqTbK7Am17jEvZNRENLByLDDGlMfmTB938HNVNl90ISWZKQeVd1p7XQtab5a9PWkKRcuU+v7QMmN8DVMb+0PLlIYY0VDwKP52IkNrSEPJ8uXL+dSnPsXGjRvp7Ozk4Ycf5qyzzgJg165drFixgrfeeosvfvGL3HXXXQA89dRT3HLLLWQyGebOncuaNWtIJBJ873vf49e//jV9fX0sW7aMH/zgB/udwO+77z527drFueeeS1NTExs3bgTgO9/5DuvWrSMej7N27VpqampYuHAhb7zxBuFwmO7ubj7xiU+U78vYkC8UywM7Syf99gFjMVp7M64Vw/sm253OUxPp7yqpj7twURsNURd134xPmuq6SkpBQx8G8nFEQoHyGJSDSWXz5S6+fT1Z3tmX5C87OmhNZsr/r+tj4f1aWqaPc7dLxxpr1Noi/jXkLSX5fJ4XXniBJ554gttuu40NGzYAsGXLFl566SWi0SjHH3881157LfF4nNtvv50NGzZQW1vLnXfeyd13383NN9/MNddcw8033wzAV77yFdatW8f5559ffp9vfOMb3H333WzcuJGmpiYAent7Of3007njjju48cYbeeihh7jppptYvnw5v/nNb1i1ahWPPvooF154oQLJCFcsWrr63HiM1uTAsRj93SVtXujoSOVIZvIkoiGvud11ldSVu0rCnDC5nvpZoQFjNcIE1VUiPlMTCVEzPsSMg3QXlVr4XGue+1t4/u02r4Uvy96eDPli8cCtLV630ZSGGLGwArZUx6BDycESdun4hRdeCMCpp55KS0tL+fHPfOYzNDQ0ADB//nx27NhBZ2cnW7du5YwzzgAgm82ydOlSADZu3Mhdd91FKpWivb2dk046ab9QciCRSISVK1eW3/+3v/0tAF/96le56667WLVqFWvWrOGhhx4a7K8tR0E65y7dLIWJNq8FY1+P21qT/QM/O/ty5ZaMcmtGLOTmoYiFOWFKfXkOCtedovEYMvoFAobxtW6w7nGT6w74nIO1trQlM7R6f3uJaGj/sS2NcZob4zQ3xGhuiDGpLkYkpMuhZegNOpRMmDCBjo6O/Y61t7cze/ZsAKLRKADBYJB8Pl9+Tun4wMestXz2s5/lF7/4xX6vl06nueqqq9i8eTMzZszg1ltvJZ1OH7Zs4XB/s+TA9z/jjDNoaWlh06ZNFAoFFixYMNhfWypQKFo6vYF8pUBR+ga3r8e77YWQjlSWfMHSUOPmiCh1mdRFQ9THIhwzvoYFUxtoqOkPICHNESEyaIdtbbFeK6QXXNqSWf7+fhe/f2Mf7d7fcmcqR308zOT6KM0NLqxMbYwzuT7mXVEULV8Kra4iGYxBh5JEIkFzczNPP/00n/70p2lvb2f9+vVcd911rFmzZlCvdfrpp3P11Vezfft25s2bR29vL++//z6TJk0CoKmpiWQyyWOPPcbq1as/8vN1dXX09PSUu28O5dJLL+VLX/oS3/3udwdVRulXKFo3B0b5yhJvn8yyrydd/pZVGvzZ3ecuYW2MR8pjL+rLM2mGy1cT1Mddt4pm1BSpvoAxjKtxE8fNI3HA5xS87lN3DsjQ0Ztl2wfd/Pnttv3m0rFYJiai7nLo+hhT6qNMro8zqc4dm1gXpcmbmE7dpQIVjin5yU9+wtVXX803v/lNAG655Rbmzp076NeZOHEijzzyCBdffDGZTAaA22+/neOOO46vfe1rLFiwgClTprBkyZID/vyVV17JihUrmDp1anmg68Fccskl3HTTTVx88cWDLudolMkX6Erl6Ei5ORYG7jt6XbdJR2+Wjt4c7Sk362dvpnDAeTISURc25k1McOox46iLuZChcRkio1NwQDfRwYILQJ93KXT/OSbH1l1d/MmbfK7Tm9slmclTHwsxIRGlKeFmzC3tx9eWZtONlGfZbYjr3DJaGWs/vArWwS1evNhu3rx5GIszfB577DHWrl3LT3/602oX5WOz1tKXK5Snx+4pbzl60nm6vX1XnzsRdA+YgbI02VPR2vIgz0TMXVWS8Gb1rI26sRn10TB1cXe1SV08RCKicRkj3bJ5rlXxj9tbq1wSkX75YrF8bupK5ehO9y950JPOkyyd07x9KlOgNhqk0QsojTWl5QEi7nZNZL8r41yXrze4PRbS8gDVd9APkhE9T8mRuvbaa3nyySd54oknhv29rLUfXYMj5y1QlnNrcaTzBTI5b/GyXJFUtkA6V6A3kyeVzdObcYub9Wbz3jHvvne7L1sgFDTURIKufzgSJB4JlvfxcNBbqyTElPoYcycm9ptGvDYa+shMnyIi1RIKBMpdRkw4/POLRUsym6c37S0J4H1B683m2dXZx/a9SfpyBVIDz59Z73am//xZ/jIWdefR0u3SlzN3TvX2YbfWUzzi1n+Khd16T7FQ/z4SChAJBQgFjM6vFRpUS0m0+VjbfNm9w1caERlWO+50V6fN/Od1VS6JiIxVLf/98wdNbIMKJcaYfcCOQzylCaisXdgEAuGJM082gaC/Wm9cBVmLtVisd9MC3r58vGhL9+HIKxUopnsjgVhtdljKP4qp3iqjequM6q0yY7PejMGYgAGvycQEMBjvuAF32wy4faAOjUKqi2BNw9EuPDafS+dad7w6jG/Raq1dcaAHBhVKDscYs9lau3jIXnCMUL1VRvVWGdVbZVRvlVG9VW4s1p1G+4iIiIgvKJSIiIiILwx1KPnhEL/eWKF6q4zqrTKqt8qo3iqjeqvcmKu7IR1TIiIiIlIpdd+IiIiILyiUiIiIiC8olIiIiIgvKJSIiIiILyiUiIiIiC8olIiIiIgvKJSIiIiILyiUiIiIiC8olIiIiIgvKJSIiIiILyiUiIiIiC8olIiIiIgvKJSIiIiILyiUiIiIiC8olIiIiIgvKJSIiIiILyiUiIiIiC8olIiIiIgvKJSIiIiILyiUiIiIiC8olIiIiIgvKJSIiIiILyiUiIiIiC8olIiIiIgvKJSIiIiIL4QG8+QVK1bY9evXD1dZRD4eayHTDX0dbsv0QCYJ+T7IZ6GYc88BCIQgGIZQDCI1EKmDWAPEG6FmgntsFDLGAGBL9SAicvSZgz0wqFDS2tr68YsiUilroecDaHsT2t+GtrehowW6d7rjvftc2Ig1QjQBkVoIRSEYdSHDBMH7UKZYAFtwYSWfdsEl0wuZLhdmwrVQNwnqmqHhGBg3y23jZ8OEuRAfV716EBEZpQYVSkSOmmwK9rwKu7fA7pdhzyvQ+gYEI9AwA+qmQGISNM2DmUuhpsm1coRiH/+9bdEFk1Q7pNpc2GndBi3PQc8u6NoJ4ThMOBYmn+S2KQth8nwXhEREpCJmMM24ixcvtps3bx7G4siYZK1r+Xjvz/Dun2Dni9D+DjTOhPFzYNxMaJzl9tG6apfWlTfVBl3vQee7butogc4d0DAdpp4C05fAtFNhysmutcYn1H0jIj5w0O4bhRI5+qyFva+5lod3fg/vPg+BIEyaD03HQtPxrpskGKl2SQenmHcBpfVNaH/L7bt2wqQTYMbpMHMZHHO6a+GpEoUSEfEBhRKpss734O2NsP130PKs62aZcrILIpMXVPWDeljl+lw42fea637au9V1Nc08A2af5faNM45acRRKRMQHFErkKMv1wY4/wJu/dVuqDaZ90gWQ5kWjN4QcTrHgunn2vOqCyu6XXZfUrLNgznIXVOqnDtvbK5SIiA8olMhR0LED3nwKtj3hxoeMnwvNn4Cpn3RXrBhNi/MR1rounw9edq0ou1+G2gkw+xwXUmad5e4PEYUSEfEBhRIZBoW8Cx/bnoQ31rvWkOmnuhDSfIq7LFcGxxbdIN8P/uZaUz54xQ3wnXMuzD0Xjln6sepVoUREfEChRIZIqh3eehpe/w289TuonQTTF8O0xW6QqlpDhlYx78aifPCyCyit21wX2JxzYc457iqfQVzdo1AiIj6gUCIVshb2bYM3/h22/cZ9e29e6FpDpi2G2qZql3Bsyadhz1YXUva86sanTDvVhZTZZ8PURYecjVahRER8QKFEBiHX5y7X3fakCyPFXH9ryJSFvpp3Y8zLJuGDV2HP311I6dnl/p1mnwOzznRzpoT6L61WKBERHxi6UPL666+TTCbLxx555BE2b97MAw88MOhSbdq0ie9///usW7eOTZs2EYlEWLZsGQCXX345K1euZPXq1YN+XalA21uwfQNsWw87vUGqpUnAGmf2T88u/pbudrPf7t0Ke1+Frvfd1U6zzoSZSzHzPgMolIhIVQ3N2jfDadOmTSQSiXIokWGW7oJ3nnVB5K3fudaRqZ+EGafBp66EiAapjkixejdJ20zv7yibdBPV7X3NDUYuWXe9m9Bt+mI3a65Cp4j4wJCOSty3bx8XXXQRS5YsYcmSJfzhD38A4IUXXmDp0qWccsopLFu2jG3btu33cy0tLTz44IPcc889LFq0iGeffRaAZ555hmXLljFnzhwee+wxAC699FJ+9atflX/2kksuYe3atUP5a4xOuT54+/ew4Xvww+XwP0+A5+5xA1PPvhEu+hEsu9Z9o1YgGT0iCdfa9clL4XP/OuABA1t+DmtWwJ0z4SdfgKfvcF123burVlwRGdsG3VLS19fHokWLyvfb29u54IILALjuuuu4/vrrOfPMM3n33Xf53Oc+x2uvvcYJJ5zAs88+SygUYsOGDXz729/m8ccfL7/GrFmz+PrXv04ikeCGG24A4OGHH2b37t0899xzvP7661xwwQWsXr2aK664gnvuuYdVq1bR1dXFH//4R3784x9/zGoYhTJJt4ZMyx+g5Rk3MHL8XLd43ILVMOnEkTeNuwydk77YfzvV7q7qadvuZt0tLXxYmuiueaG7PX4OBH3TuCoio9CgzzDxeJwtW7aU75fGlABs2LCBrVu3lh/r7u4mmUzS1dXFZZddxptvvokxhlwud0TvtWrVKgKBAPPnz2fPnj0AnHPOOVx11VXs27ePxx9/nIsuuohQaIyfKK11V2Hs3OzWkXn3eWjf7laxnXg8HP+PcNZ/gXBNtUsqflQz3s1/csxSd99atzJy+9vQ8Q68+L/cxHipVrcm0aT5MPFE93+r6Th3TIOfRWQIDOmnebFY5PnnnycW23/5+GuuuYZzzz2XX/7yl7S0tLB8+fIjer1otP9EN3Bg3qWXXsrPfvYzHn30UdasWTMkZR8xrIXuXa7lY9dLrjVk1xbXDVP6kPjkV2DCPLWESGWMccsAJCa5BQRLcn39KyPve93NV9O9E3o+gMRkGDfLzdw7fo673TADGo+BmgkasyIiR2RIQ8l5553H/fffz7e+9S0AtmzZwqJFi+jq6mLatGmAa1k5kLq6Orq7u4/ofS6//HJOO+00pkyZwvz584ek7L6U6YF9b7g1Uva86qYg3/uqCyZNx8K42W6w4ilfcYu86cQvwykcd6G36bj9jxfz0LMHena7bddL7lLy3r3ueD4DdVOgvhnqprq1feqmuIn3apvcVjMB4uNca57+H4uMWUMaSu677z6uvvpqFi5cSD6f5+yzz+bBBx/kxhtv5LLLLuP222/n85///AF/9vzzz2f16tWsXbuW+++//5DvM3nyZE488URWrVo1lMWvjlyf++bZ0eKay1vfcKvKtr4J6U63gmzjTKifDvP+AU77GsTH68Qt/hEIQcM0tx1Irs91/fS2QV+7G8PS/Rd3BVh563b/3wFiDRCtdwsVlrZIwk2vH67xtphbaToY8bYwBMIQCLrymIC3Gfa/+tC6UG+LYAtuXyx4W37/rZBzc/QU8m6fz7hjheyAx7zbpePFnPez+f7bpde2hf73ssUBmwWKB6g4r+yl3yUQABPs/x1LWzDsbREIRt2+VD+lLVLjQmW4xrtf238/Uuvdr3H1XHpM5xipghE5eVoqleLkk0/mr3/9Kw0NDdUuzsHl+iC5F5J7XJdLz27o2gmd70HXu+52uss1k9c1uy0xGRqmQ/00d1zTtssQMnPOBsC+/UyVS3IA1kIhA9ne/i2XcrPY5vrcPp+G/IBQsF+QKHgf/EWgCBbKH/aWAdkk4H3eex/2BNwHvTH9gSYQHBAASrfDHw0EgdCHjpWeG9o/IAWC3vsF+8PSQYNTuUK8wFIKUgNCVDlMDQg9xdwhglSmP0Dlvdv5zIA6TUMu7eo71+d+Nlzjwkwk0b9FE15Q9EJjrB4ipfBYek5d/z6acI9rgLTsz//zlBypDRs2cMUVV3D99dcfnUBSLLi5HjJJt093Q2bAN7y+TvfNL9Xmvg2m2tzW2+r+6GvGu6bpmvGueTo+wXW9zDoDaie6Vo9AcPh/DxG/M6b/m33N0K2MLBUoFiDf5wLKfluq/3i6040nKoea1P7Py/V5wbLPhbNI7YCQU7t/0Bl4P1LT34ITjkO4tB+wheJea1DctRKpVWfUGFxLydSg3Xyl5rAQGanMbW7clr2lvsolEZEx69auoZlm3hizD9hxiKc0Aa1HXrIjMy5Gw5xxgXlD/boHYsEWLUVrKRa9bbjfs73PhsbHTX6432e0Ub1VRvVWGdVbZYa73oxxnWHGuI4yYzDGEDCH6CLwu460bXu7w7YwTJ+pPtBqrV1xoAcGFUoOxxiz2Vq7eMhecIxQvVVG9VYZ1VtlVG+VUb1VbizWnUZRioiIiC8olIiIiIgvDHUo+eEQv95YoXqrjOqtMqq3yqjeKqN6q9yYq7shHVMiIiIiUil134iIiIgvKJSIiIiILwx5KDHG/Isx5mVjzBZjzFPGmKlD/R6jkTHmfxhjXvfq7pfGmMZql2kkMMb8B2PMq8aYojFmTF06VwljzApjzDZjzHZjzH+tdnlGAmPMj4wxe40xr1S7LCOJMWaGMWajMWar9zd6XbXLNBIYY2LGmBeMMX/z6u22apfpaBryMSXGmHprbbd3+xvAfGvt14f0TUYhY8x5wNPW2rwx5k4Aa+0/V7lYvmeMORG3wMkPgBustdVfnMmnjDFB4A3gs8BO4EXgYmvt1qoWzOeMMWcDSeAn1toF1S7PSGGMaQaarbV/NcbUAX8BVun/26EZYwxQa61NGmPCwHPAddba56tctKNiyFtKSoHEU4u3LJYcmrX2KWttadbD54Hp1SzPSGGtfc1au63a5RghTgO2W2vfttZmgUeBL1S5TL5nrX0GaK92OUYaa+1ua+1fvds9wGvAQZaSlhLrJL27YW8bM5+jwzKmxBhzhzHmPeAS4ObheI9R7p+AJ6tdCBl1pgHvDbi/E31IyFFgjJkFnAL8ucpFGRGMMUFjzBZgL/Bba+2YqbeKQokxZoMx5pUDbF8AsNZ+x1o7A/g5cM1QFngkO1y9ec/5DpDH1Z1wZPUmIv5kjEkAjwP/+UMt6XIQ1tqCtXYRrsX8NGPMmOk2DFXyQ9bafzjCp/4ceAK4pZL3GW0OV2/GmMuBlcBnrCaQKRvE/zc5tPeBGQPuT/eOiQwLb0zE48DPrbX/r9rlGWmstZ3GmI3ACmBMDLQejqtvjh1w9wvA60P9HqORMWYFcCNwgbU2Ve3yyKj0InCsMWa2MSYC/Cfg36pcJhmlvAGbDwOvWWvvrnZ5RgpjzMTS1ZfGmDhuYPqY+RwdjqtvHgeOx10RsQP4urVW38YOwxizHYgCbd6h53XV0uEZY74I3A9MBDqBLdbaz1W1UD5mjPlH4F4gCPzIWntHdUvkf8aYXwDLccvI7wFusdY+XNVCjQDGmDOBZ4G/4z4PAL5trX2ieqXyP2PMQuDHuL/RAPB/rLXfq26pjh5NMy8iIiK+oBldRURExBcUSkRERMQXFEpERETEFxRKRERExBcUSkRERMQXFEpE5KCMMcnDP6v83FuNMTcM1+uLyOinUCIiIiK+oFAiIoNijDnfGPNnY8xL3rpEkwc8/AljzJ+MMW8aY7424Ge+ZYx50RjzsjHmtioUW0RGAIUSERms54DTrbWnAI/ilkcoWQh8GlgK3GyMmWqMOQ84FjgNWAScaow5++gWWURGgooW5BORMW068L+NMc1ABHhnwGNrrbV9QJ+3kNhpwJnAecBL3nMSuJDyzNErsoiMBAolIjJY9wN3W2v/zRizHLh1wGMfXrfCAgb4b9baHxyV0onIiKXuGxEZrAagtMjmZR967AvGmJgxZgJuEbsXgX8H/skYkwAwxkwzxkw6WoUVkZFDLSUicig1xpidA+7fjWsZ+b/GmA7gaWD2gMdfBjbiVtT9F2vtLmCXMeZE4E9uNXuSwJeBvcNffBEZSbRKsIiIiPiCum9ERETEFxRKRERExBcUSkRERMQXFEpERETEFxRKRERExBcUSkRERMQXFEpERETEF/4/wPslX8783sEAAAAASUVORK5CYII=\n", + "text/plain": [ + "
" + ] + }, + "metadata": { + "needs_background": "light" + }, + "output_type": "display_data" + } + ], + "source": [ + "# Ridge plots\n", + "classify.density_by(X,'Label')" + ] + } + ], + "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.8.5" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +} diff --git a/setup.py b/setup.py index e181044..1c4dc49 100644 --- a/setup.py +++ b/setup.py @@ -13,7 +13,7 @@ setup( name='itk-hasi', - version='0.2.0', + version='0.2.1', author='Kitware Medical', author_email='itk+community@discourse.itk.org', packages=['itk'], @@ -47,5 +47,6 @@ r'itk>=5.2.0', r'itk-boneenhancement', r'itk-ioscanco', + r'dwd' ] ) diff --git a/src/hasi/hasi/classify.py b/src/hasi/hasi/classify.py new file mode 100644 index 0000000..8e8cba4 --- /dev/null +++ b/src/hasi/hasi/classify.py @@ -0,0 +1,114 @@ +#!/usr/bin/env python3 + +# Copyright NumFOCUS +# +# Licensed under the Apache License, Version 2.0 (the "License"); +# you may not use this file except in compliance with the License. +# You may obtain a copy of the License at +# +# http://www.apache.org/licenses/LICENSE-2.0.txt +# +# Unless required by applicable law or agreed to in writing, software +# distributed under the License is distributed on an "AS IS" BASIS, +# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +# See the License for the specific language governing permissions and +# limitations under the License. + +# Purpose: Python functions for shape classification and distance +# analysis with DWD classifier +import itk +import numpy as np +from dwd.dwd import DWD + +import seaborn as sns +import matplotlib.pyplot as plt + +# Generates an nxd feature matrix with +# - n = number of meshes (samples) +# - d = number of features, equal to +# number of mesh points * (1 / step) * mesh dimension +def make_point_features(meshes:list, + step:int=1) -> np.ndarray: + assert(step >= 1) + + features = None + for mesh in meshes: + points = [mesh.GetPoint(idx) + for idx in range(mesh.GetNumberOfPoints()) + if idx % step == 0] + points_array = np.expand_dims(np.asarray(points).flatten(),0) + + if features is None: + features = points_array + else: + features = np.append(features, points_array, axis=0) + return features + +def get_distances(classifier:DWD, features:np.ndarray) -> np.ndarray: + direction = classifier.coef_.reshape(-1) + intercept = float(classifier.intercept_) + distance = features.dot(direction) + intercept + + return distance + +def densityplot(X, y, ax=None): + ax = ax or plt.gca() + + xlim = X.min() * 1.05, X.max() * 1.05 + ax.set_xlim(*xlim) + + ax.axvline(x=0, color='green', linestyle='dashed') + + order = sorted(np.unique(y)) + sns.kdeplot(ax=ax, x=X, + color='black', + common_norm=True, + common_grid=True, + bw_adjust=1.0,) + + sns.kdeplot(ax=ax, x=X, + hue=y, hue_order=order, + common_norm=True, + common_grid=True, + bw_adjust=1.0) + + sns.rugplot(ax=ax, x=X, + hue=y, hue_order=order,) + + sns.despine() + +# Ridge Plots +def density_by(X, category, categories=None): + # adapted from https://seaborn.pydata.org/examples/kde_ridgeplot.html + + plt.close() + + variable = 'Distance' + if not categories: + categories = X[category].unique().tolist() + + title = '{} density by {}'.format(variable, category) + + g = sns.FacetGrid(X, + row=category, + hue=category, + row_order=categories, + hue_order=categories, + aspect=8, + height=1) + + g.map(sns.kdeplot, variable, fill=True) + g.map(plt.axhline, y=0, lw=2) + g.map(plt.axvline, x=0, lw=2, color='k') + + def label(x, color, label): + ax = plt.gca() + ax.text(0, .3, label, + ha='left', va='center', + transform=ax.transAxes) + + g.map(label, category) + + g.set_titles("") + g.set(yticks=[]) + g.despine()