{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Table of Contents \n",
"\n",
" **1 Imports and predefined routines** \n",
"\n",
"\n",
" **2 Loading the data** \n",
"\n",
" **3 Dimension reduction**\n",
"\n",
" **4 Clustering**\n",
"\n",
" **5 Clustering outputs**"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## 1 Imports and predefined routines"
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {},
"outputs": [
{
"name": "stderr",
"output_type": "stream",
"text": [
"WARNING:hyperspy.api:The ipywidgets GUI elements are not available, probably because the hyperspy_gui_ipywidgets package is not installed.\n",
"WARNING:hyperspy.api:The traitsui GUI elements are not available, probably because the hyperspy_gui_traitsui package is not installed.\n"
]
}
],
"source": [
"%matplotlib inline\n",
"\n",
"import hyperspy.api as hs\n",
"import numpy as np\n",
"import matplotlib.pyplot as plt\n",
"from mpl_toolkits.mplot3d import Axes3D\n",
"from skcmeans.algorithms import Probabilistic\n",
"from sklearn import preprocessing\n",
"import h5py\n",
"from sklearn.decomposition import PCA, NMF, TruncatedSVD\n",
"from sklearn.cluster import DBSCAN\n",
"from sklearn.mixture import GaussianMixture\n",
"from scipy import sparse\n",
"\n",
"from matplotlib.colors import to_rgb, LinearSegmentedColormap\n",
"from skcmeans.algorithms import Probabilistic, Possibilistic, GustafsonKesselMixin\n",
"from skimage.transform import hough_line, hough_line_peaks, rescale\n",
"from skimage.color import label2rgb\n",
"\n",
"LINEWIDTH = 20\n",
"def full_width_figure(aspect_ratio):\n",
" return plt.figure(figsize=(LINEWIDTH, aspect_ratio * LINEWIDTH))\n",
"def half_width_figure(aspect_ratio):\n",
" return plt.figure(figsize=(0.5 * LINEWIDTH, 0.5 * aspect_ratio * LINEWIDTH))\n",
"def quarter_width_figure(aspect_ratio):\n",
" return plt.figure(figsize=(0.25 * LINEWIDTH, 0.25 * aspect_ratio * LINEWIDTH))\n",
"\n",
"MPL_COLORS_RGB = [to_rgb('C{}'.format(i)) for i in range(10)]"
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {},
"outputs": [],
"source": [
"# for the scatter plots\n",
"def scatter_loadings(loadings, c='k', aspect_ratio=1):\n",
" fig = quarter_width_figure(aspect_ratio)\n",
" gridspec = plt.GridSpec(2, 2)\n",
" ax01 = fig.add_subplot(gridspec[0, 0], aspect='equal')\n",
" ax01.scatter(loadings[:, 0], loadings[:, 1], s=0.25, c=c, edgecolor='none')\n",
" ax01.set_xlabel('Loading 0')\n",
" \n",
" ax21 = fig.add_subplot(gridspec[0, 1], aspect='equal')\n",
" ax21.scatter(loadings[:, 2], loadings[:, 1], s=0.25, c=c, edgecolor='none')\n",
" ax21.set_xlabel('Loading 2')\n",
" ax21.set_ylabel('Loading 1')\n",
" \n",
" ax23 = fig.add_subplot(gridspec[1, 1], aspect='equal')\n",
" ax23.scatter(loadings[:, 2], loadings[:, 3], s=0.25, c=c, edgecolor='none')\n",
" ax23.set_ylabel('Loading 3')\n",
" \n",
" ax03 = fig.add_subplot(gridspec[1, 0], aspect='equal')\n",
" ax03.scatter(loadings[:, 0], loadings[:, 3], s=0.25, c=c, edgecolor='none')\n",
" \n",
" width_ratio = np.diff(ax21.get_xlim()) / np.diff(ax01.get_xlim())\n",
" height_ratio = np.diff(ax03.get_ylim()) / np.diff(ax01.get_ylim())\n",
" gridspec.set_width_ratios([1, width_ratio])\n",
" gridspec.set_height_ratios([1, height_ratio])\n",
" \n",
" plt.tight_layout()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
" **Return to Table of Contents**\n",
"## 2 Loading the data "
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {},
"outputs": [],
"source": [
"im = hs.load('dataset1_CBED.hspy')"
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {
"scrolled": false
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"\n"
]
}
],
"source": [
"# Check the data dimensions of the data\n",
"print(im)"
]
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {},
"outputs": [],
"source": [
"# Change the data over from hs-list to np-array\n",
"data1=[]\n",
"for i in range (0,im.data.shape[0]):\n",
" for j in range (0,im.data.shape[1]):\n",
" a=im.inav[j,i].data.ravel()\n",
" data1.append(a)\n",
"data=np.array(data1)\n",
"data = data.astype('float64') "
]
},
{
"cell_type": "code",
"execution_count": 6,
"metadata": {},
"outputs": [],
"source": [
"# standard scaling\n",
"data -= data.mean()\n",
"data /= data.std()\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
" **Return to Table of Contents**\n",
"\n",
"## 3 Perform PCA dimension reduction "
]
},
{
"cell_type": "code",
"execution_count": 7,
"metadata": {},
"outputs": [],
"source": [
"decomp = PCA(n_components=10).fit(data)\n",
"loadings=decomp.transform(data)"
]
},
{
"cell_type": "code",
"execution_count": 8,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"[]"
]
},
"execution_count": 8,
"metadata": {},
"output_type": "execute_result"
},
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAXQAAAD4CAYAAAD8Zh1EAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAAAag0lEQVR4nO3de5Sc9X3f8fd3Ljuz98vsLrqsbrOAhQAJ8AK7dRLqkBNDfZFxCYE2zolbl9Ia241JMLk4TUI5NanJpSkxh0NTn7ZpOBRjDsEYnKTk1CkIa+VIQkIgdINddFtpV9r7ZWZ+/WNmV7PLSju7mtUzzzOf1zl75rn8NPvVHOnzPPN9fvOMOecQERH/C3ldgIiIFIcCXUQkIBToIiIBoUAXEQkIBbqISEBEvPrFzc3Nbv369V79ehERX9qxY8cp51zLfPs8C/T169fT3d3t1a8XEfElM3vvfPvUchERCQgFuohIQCjQRUQCQoEuIhIQCnQRkYBQoIuIBIQCXUQkIHwX6PtPDPEfXnyL8am016WIiJQU3wV678AoT/39YX7y3oDXpYiIlBTfBfqN65sIh4xth057XYqISEnxXaDXxqNcs7qe1xXoIiKz+C7QATqTTezsOcPYpProIiLTfBnoXckEU2lH93v9XpciIlIyfBnoN65vIhIyXj+otouIyDRfBnp1LMLmNvXRRUTy+TLQAbraE+zuPcvIRMrrUkRESoJ/Az3ZTDrj2H5EfXQREfBxoH90XSPRsKntIiKS49tAr6wIc92aBrbpwqiICODjQIfs9MU3PzjL4PiU16WIiHjO14He2Z4g42D7YfXRRUR8Heg3rG2kIhLSfV1ERPB5oMejYa5f06ALoyIi+DzQITsffe/RQc6Oqo8uIuXN/4GeTOAcvHFYZ+kiUt58H+jXrW0gFgmp7SIiZc/3gR6LhOlY38i2Q5rpIiLlzfeBDtm2y75jgwyMTHpdioiIZwIR6J3JBKA+uoiUt0AE+ua2BiqjYd0fXUTKWiACvSISomN9oy6MikhZC0SgQ3Y++v4Tw5wanvC6FBERTxQU6GZ2m5m9Y2YHzOyhC4y70czSZnZn8UosTNd0H12zXUSkTC0Y6GYWBh4Hbgc2AfeY2abzjHsUeKXYRRbi2tX1VFeEef3QKS9+vYiI5wo5Q78JOOCcO+ScmwSeBrbOM+7LwHeBk0Wsr2CRcIgbNzTpwqiIlK1CAn010JO33pvbNsPMVgN3AE8Ur7TF60omONg3wsnBcS/LEBHxRCGBbvNsc3PW/xj4unMufcEnMrvXzLrNrLuvr6/AEgvX1Z7to2u2i4iUo0ICvRdYk7feBhydM6YDeNrMjgB3An9mZp+d+0TOuSedcx3OuY6WlpalVXwBV6+qpzYe0f3RRaQsRQoYsx24wsw2AB8AdwP/LH+Ac27D9LKZfQd40Tn3fPHKLEw4ZNy8oUn3dRGRsrTgGbpzLgXcT3b2yj7gGefcXjO7z8zuW+4CF6szmeDwqRGOn1UfXUTKSyFn6DjnXgJemrNt3gugzrlfufiylm76vi6vHzrFHde3eVmKiMglFZhPik7btLKO+sqopi+KSNkJXKCHcn10zXQRkXITuECH7PTFnv4xegdGvS5FROSSCWygA5rtIiJlJZCBfmVrLY1V6qOLSHkJZKCHQkZnMsG2Q6dxbu6HWkVEgimQgQ7ZtssHZ8bo6R/zuhQRkUsiuIGeNx9dRKQcBDbQL2+tobmmQhdGRaRsBDbQzYybkwleP6g+uoiUh8AGOmTbLscHxzlyWvPRRST4gh3o0/dH1/RFESkDgQ70ZHM1rbUx3QZARMpCoAPdzOhqVx9dRMpDoAMdsn30U8MTHOwb8boUEZFlFfhAP3d/dLVdRCTYAh/o6xJVrKyPs00XRkUk4AIf6GZGl+7rIiJlIPCBDtDZnuD0yCT7Twx7XYqIyLIpi0Cfvq/LNvXRRSTAyiLQ1zRV0dZYqQ8YiUiglUWgQ3a2y7bDp8lk1EcXkWAqm0DvSiY4MzrF28eHvC5FRGRZlE+gt2s+uogEW9kE+qqGStYlqtRHF5HAKptAh2zb5ceHT5NWH11EAqisAr0zmWBwPMW+Y4NelyIiUnRlFei6P7qIBFlZBfpldXGSzdW6MCoigVRWgQ7Z2wD8+HA/qXTG61JERIqq7AK9K5lgeCLFnqPqo4tIsJRdoHfqvi4iElBlF+gttTEub63RhVERCZyyC3TItl22H+lnSn10EQmQ8gz09gSjk2l29571uhQRkaIpKNDN7DYze8fMDpjZQ/Ps32pmu81sp5l1m9lPFb/U4lEfXUSCaMFAN7Mw8DhwO7AJuMfMNs0Z9rfAFufcdcC/AJ4qcp1F1VRdwcYVtQp0EQmUQs7QbwIOOOcOOecmgaeBrfkDnHPD7twXdlYDJX+zlM5kgu4jA0ym1EcXkWAoJNBXAz156725bbOY2R1m9jbwfbJn6R9iZvfmWjLdfX19S6m3aDqTCcam0uzqPeNpHSIixVJIoNs82z50Bu6c+55zbiPwWeDh+Z7IOfekc67DOdfR0tKyqEKLrTPZhJnu6yIiwVFIoPcCa/LW24Cj5xvsnPu/QLuZNV9kbcuqoaqCq1bUKdBFJDAKCfTtwBVmtsHMKoC7gRfyB5jZ5WZmueUbgAqg5JOyqz3BjvcHGJ9Ke12KiMhFWzDQnXMp4H7gFWAf8Ixzbq+Z3Wdm9+WG/VNgj5ntJDsj5hfzLpKWrK5kgslUhp09Z7wuRUTkokUKGeScewl4ac62J/KWHwUeLW5py+/GDU2Ecn306bnpIiJ+VZafFJ1WXxnl6lX1uj+6iARCWQc6ZPvoO98/oz66iPieAj2ZYDKdYcd7A16XIiJyUco+0G/c0EQ4ZJq+KCK+V/aBXhOLcO3qet3XRUR8r+wDHbK3AdjVe4bRyZTXpYiILJkCneyF0am0o/uI+ugi4l8KdKBjXSORkGn6ooj4mgIdqI5F2LKmQRdGRcTXFOg5XckEb35wluEJ9dFFxJ8U6DmdyQTpjGP7kX6vSxERWRIFes5H1zUSDRvb1HYREZ9SoOdUVoS5fk2jLoyKiG8p0PN0tifY88FZBsenvC5FRGTRFOh5upIJMg5+fEh9dBHxHwV6nuvXNlARCek2ACLiSwr0PPFomBvWNqiPLiK+pECfoyvZzFvHBjkzOul1KSIii6JAn6OrPYFz8MZh9dFFxF8U6HNsWVNPPBrSbQBExHcU6HPEImE61jXpwqiI+I4CfR6dySbePj5E/4j66CLiHwr0eXS1JwB4Q2fpIuIjCvR5bG5roKoirOmLIuIrCvR5RMMhOtY36cKoiPiKAv08upIJ3j05TN/QhNeliIgURIF+HjN99MM6SxcRf1Cgn8c1q+qoiUXUdhER31Cgn0ckHOLG9bo/uoj4hwL9ArraExzqG+HE4LjXpYiILEiBfgFdyWYAfWpURHxBgX4Bm1bVURtXH11E/EGBfgHhkHHzBt3XRUT8QYG+gM5kgiOnRzl2dszrUkRELkiBvoDp+ehqu4hIqSso0M3sNjN7x8wOmNlD8+z/52a2O/fzmpltKX6p3rhqRR0NVVEFuoiUvAUD3czCwOPA7cAm4B4z2zRn2GHgFufcZuBh4MliF+qVUK6PrvnoIlLqCjlDvwk44Jw75JybBJ4GtuYPcM695pwbyK1uA9qKW6a3upIJegfG6Okf9boUEZHzKiTQVwM9eeu9uW3n8y+BH8y3w8zuNbNuM+vu6+srvEqPdeb66JrtIiKlrJBAt3m2uXkHmn2cbKB/fb79zrknnXMdzrmOlpaWwqv02JWttTRVV6jtIiIlLVLAmF5gTd56G3B07iAz2ww8BdzunAtU8oVCRmeyiW0HT+Ocw2y+Y5yIiLcKOUPfDlxhZhvMrAK4G3ghf4CZrQWeAz7vnNtf/DK915VMcPTsOO+rjy4iJWrBM3TnXMrM7gdeAcLAnzvn9prZfbn9TwC/AySAP8udvaaccx3LV/allz8ffV2i2uNqREQ+rJCWC865l4CX5mx7Im/5i8AXi1taaWlvqaG5Jsa2Q6e5+6a1XpcjIvIh+qRogcyyffTXD2X76CIipUaBvghd7QlODE5w+NSI16WIiHyIAn0RupK5PrqmL4pICVKgL8KG5mouq4vpvi4iUpIU6ItgZnQlE2w71K8+uoiUHAX6InUmE5wanuBg37DXpYiIzKJAXyTdH11ESpUCfZHWNlWxqj6uC6MiUnIU6ItkZnS2Z/vomYz66CJSOhToS9CVTNA/Msn+k0NelyIiMkOBvgSdSfXRRaT0KNCXYE1TFW2NlfrCCxEpKQr0JepKJnj94GnOjk15XYqICKBAX7LPd61jZDLNwy++5XUpIiKAAn3JNrc18G9uaefZHb38zVsnvC5HRESBfjG+cusVbFxRy0PPvcnAyKTX5YhImVOgX4SKSIg/vOs6zo5N8jsv7PW6HBEpcwr0i7RpVR1fvfUK/mrXUV7c/aHvzhYRuWQU6EVw3y3tbGmr5xvP76FvaMLrckSkTCnQiyASDvHYXVsYmUzzG8+9qVvriognFOhFcnlrLQ9+4iP8zb4TPPeTD7wuR0TKkAK9iL7wsQ3cuL6R3/2rvRw7O+Z1OSJSZhToRRQOGd/6hS2k0o4Hn92t1ouIXFIK9CJbl6jmNz95FT969xT/68fve12OiJQRBfoy+KWb1/JTlzfzyPf38f7pUa/LEZEyoUBfBmbGo3duJmzGrz+7S1+EISKXhAJ9maxuqOQbn97EG4f7+c5rR7wuR0TKgAJ9Gf3CR9u4dWMrj778Ngf7hr0uR0QCToG+jMyM//i5a6msCPPAM7tIpTNelyQiAaZAX2atdXF+f+s17Ow5w5M/OuR1OSISYAr0S+DTm1fyyWtX8kd/vZ+3jw96XY6IBJQC/RIwMx7+7DXUV0Z54JldTKbUehGR4lOgXyJN1RU8cse17D06yH959YDX5YhIACnQL6FPXL2Cz92wmsdfPcCbvWe9LkdEAkaBfon9+09fTUtNjK89s5PxqbTX5YhIgBQU6GZ2m5m9Y2YHzOyhefZvNLPXzWzCzH6t+GUGR31llEfv3My7J4f5o7/e73U5IhIgCwa6mYWBx4HbgU3APWa2ac6wfuArwLeKXmEA3XJlC/fctJYnf3SIHe/1e12OiAREIWfoNwEHnHOHnHOTwNPA1vwBzrmTzrntwNQy1BhIv/XJq1jdUMkDz+xidDLldTkiEgCFBPpqoCdvvTe3bdHM7F4z6zaz7r6+vqU8RWDUxCL8pzu3cOT0KH/w8jtelyMiAVBIoNs825Z0+0Dn3JPOuQ7nXEdLS8tSniJQutoTfOFj6/nOa0d47cApr8sREZ8rJNB7gTV5623A0eUpp/w8+ImNJJur+fVndzM0ro6ViCxdIYG+HbjCzDaYWQVwN/DC8pZVPiorwnzrri0cOzvGI9/f53U5IuJjCwa6cy4F3A+8AuwDnnHO7TWz+8zsPgAzW2FmvcDXgN82s14zq1vOwoPkhrWN/Otb2nl6ew+vvn3S63JExKfMqy8y7ujocN3d3Z787lI0kUrzmT/9fwyMTvLDX/0ZGqoqvC5JREqQme1wznXMt0+fFC0RsUiYx+7aQv/IJL/7wl6vyxERH1Kgl5BrVtfz5Z+9gud3HuXlPce8LkdEfEaBXmL+7cfbuXZ1Pb/1vT2cGp7wuhwR8REFeomJhkM8dtcWhsZT/Pb39uDVNQ4R8R8Fegm68rJavvbzV/Ly3uO8sEtT/kWkMAr0EvWvfjrJDWsb+MbzezgxOO51OSLiAwr0EhUOGY/ddR2T6Qxf/+5utV5EZEEK9BK2obmah27byN+908cz3T0L/wERKWsK9BL3y13r6UomePjFffQOjHpdjoiUMAV6iQuFjD+4czPOOR58djeZjFovIjI/BboPrGmq4huf2sRrB0/zP7a953U5IlKiFOg+8Ys3ruEff6SFb/7gbQ6fGvG6HBEpQQp0nzAzvvm5zUTDxq/9712k1XoRkTkU6D6yoj7O72+9hh3vDfDUjw55XY6IlBgFus9svW4Vn7j6Mh774X72nxjyuhwRKSEKdJ8xMx6541pq4hEeeGYXU+mM1yWJSIlQoPtQc02MRz57DW9+cJZv/91Br8sRkRKhQPep269dydbrVvGf//Zd9nxw1utyRKQERLwuQJbu9z5zNa8fPM2n/vTvaa6poK2xijVNVbQ1VrKmsYo1TZW0NVaxqiFOLBL2ulwRWWYKdB9rqKrgL+/t5OU9x+npH6VnYJTdvWf4wZvHSOVNazSDy2rjrGnKBn1bYyVtecG/sj5OJKw3ayJ+p0D3ufaWGr708ctnbUtnHMcHx+ntH6VnYIye/lF6B8boGRjljcP9PL9zjPxp7OGQsbI+PuusfuaxsYrW2hihkF3iv5mILJYCPYDCIWN1QyWrGyq5eZ79k6kMx86OZUM+d2Y/vfzqO330Dc3+6ruKcIjVjZXZM/q8M/vp9UR1BWYKfBGvKdDLUEUkxLpENesS1fPuH59K0zswRu9A9gy/N+8Mf8+bxxgYnZo1vjIapq2xkpUNlbTWxrisLkZrbZzL6mK0zDzG1McXWWYKdPmQeDTM5a01XN5aM+/+4YlUNuz7x2YeewZGOTE4zv7jQ/QNT8x7a4KGqiiX1cZpzQV+a12My2pjtNbFcweCOC21MeJRBb/IUijQZdFqYhE2rqhj44q6efenM47+kUlODI7TNzTByaFxTgzmP05w8OQpTg5NzLp4O62+MkprbSwX+HFaco+tddnQb63NHhAqKxT8IvkU6FJ04ZDRUptts1xIJuMYGJ2cCfuTQxOcHMw+nsg9vnG4n5ND40ylPxz8tfHIzJl9a96Zfn1llNp4lNp4hJpYJPsYj1AXjxKLhNTvl8BSoItnQiEjURMjURNjE/Of7QM45zgzOsWJoXFODp4L+7684O9+b4CTQxNMpi58K4RIyGYCvjYWzT2eC/3aeJSaWIS63HpN7NyBoS6eGx+PENU0TylBCnQpeWZGY3UFjdUVbFxx/nHOOQbHUgyOTzE4PsXweIqh8RTDEymGJlIMzd02nt12fHCcd09Ob5ua993AXLFIiNq8A8C5dwN5B4DK7EGgNh6lrjK7ry4eoa4yO0YXiaXYFOgSGGZGfVWU+qrokp/DOcdEKjMT+MO50B+aWZ86dzDI2zY0nuL9kdGZg8TwRIqFblmfPShE84I/G/Z18bnr55Zrc/vqKqNUV4TVPpJZFOgiecyMeDRMPBqmuebC1wAuxDnHyGSaofGpmXcN08tD41MMjufeSczsSzE4NsXRM2PZfWNTTCzQPgoZM+8IznsAiEWIRUPEIiFikTAVkXPLsWiIinAotz9MLBKatT8aNh0wfEaBLrIMzIyaWLb1srJ+ac8xkUrnzvizAT80cxA4tzy9b/oA0dM/OrNveCKFu4gvtjLLfqgsFgkRi+YHfjgX+nnrMweN+ffHo9nnqMz9xKNhKitCMwfP/O2xSEifTF4iBbpIiYpFwsRqlv5OIZNxjEymmExlmEhlZh4nUum89TQTU5nZ6zPjLrx/MpVmeCLF6eFJJlJpJtOZmbHT45Z6QIlHQ7OCPnswCFFZESYeCROvmD4AhOaMCWfH5LbPOmBUhImGQ4QMQmaEQnZu2bLL4VD2Xcn0csgMMwhPjynxA40CXSSgQiGjNr706wkXyzlHKpO9JjE+lc77yTA2lWZsMs1Y3vaxyTTjqUz2cXrbVJqxqcysMWdGp2Y9z/S4i3k3shjzHRDCueAPhaaXjXAo72ARmn3guOemtXzxp5NFr02BLiLLwsyIho1oOERNbHmjZvpi9sT0wWLmAJFmfDI9sy2VdmScI51xOEd22TkyLvuOJjPfsnNkMufGudyfn7ucHT/9/Nl9+cvpvOdb6DMaS1XQq2xmtwF/AoSBp5xz35yz33L7/wkwCvyKc+4nRa5VRGRe+Rez6/HuXYnXFvx0hJmFgceB24FNwD1mtmnOsNuBK3I/9wLfLnKdIiKygEI+7nYTcMA5d8g5Nwk8DWydM2Yr8N9d1jagwcxWFrlWERG5gEICfTXQk7fem9u22DEiIrKMCgn0+ebpzL2eXMgYzOxeM+s2s+6+vr5C6hMRkQIVEui9wJq89Tbg6BLG4Jx70jnX4ZzraGlpWWytIiJyAYUE+nbgCjPbYGYVwN3AC3PGvAD8smV1Amedc8eKXKuIiFzAgtMWnXMpM7sfeIXstMU/d87tNbP7cvufAF4iO2XxANlpi19YvpJFRGQ+Bc1Dd869RDa087c9kbfsgC8VtzQREVkMc5fq87Jzf7FZH/DeEv94M3CqiOX4nV6P2fR6nKPXYrYgvB7rnHPzXoT0LNAvhpl1O+c6vK6jVOj1mE2vxzl6LWYL+uuh79ESEQkIBbqISED4NdCf9LqAEqPXYza9HufotZgt0K+HL3voIiLyYX49QxcRkTkU6CIiAeG7QDez28zsHTM7YGYPeV2Pl8xsjZm9amb7zGyvmX3V65q8ZmZhM/sHM3vR61q8ZmYNZvasmb2d+zfS5XVNXjGzX839H9ljZn9pZnGva1oOvgr0Ar9so5ykgAecc1cBncCXyvz1APgqsM/rIkrEnwAvO+c2Also09fFzFYDXwE6nHPXkL2Fyd3eVrU8fBXoFPZlG2XDOXds+qv+nHNDZP/Dlu196M2sDfgk8JTXtXjNzOqAnwH+K4BzbtI5d8bTorwVASrNLAJUMc/dYIPAb4GuL9I4DzNbD1wPvOFxKV76Y+BBIONxHaUgCfQB/y3XgnrKzKq9LsoLzrkPgG8B7wPHyN4N9ofeVrU8/BboBX2RRrkxsxrgu8C/c84Nel2PF8zsU8BJ59wOr2spERHgBuDbzrnrgRGgLK85mVkj2XfyG4BVQLWZ/ZK3VS0PvwV6QV+kUU7MLEo2zP/COfec1/V46GPAZ8zsCNlW3M+a2f/0tiRP9QK9zrnpd2zPkg34cvRzwGHnXJ9zbgp4DvhHHte0LPwW6IV82UbZMDMj2yPd55z7Q6/r8ZJz7jecc23OufVk/138H+dcIM/CCuGcOw70mNlHcptuBd7ysCQvvQ90mllV7v/MrQT0AnFB90MvFef7sg2Py/LSx4DPA2+a2c7ctt/M3b9e5MvAX+ROfg5Rpl8845x7w8yeBX5CdmbYPxDQWwDoo/8iIgHht5aLiIichwJdRCQgFOgiIgGhQBcRCQgFuohIQCjQRUQCQoEuIhIQ/x93BD9IwvYfpwAAAABJRU5ErkJggg==\n",
"text/plain": [
"