sampler_example.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "02fb2313",
   "metadata": {},
   "outputs": [],
   "source": [
    "from mdsampler.sampler import Sampler\n",
    "from pymatgen.core import Lattice, Structure"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7caca452",
   "metadata": {},
   "source": [
    "## Initialize the MD sampler"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "2865a270",
   "metadata": {},
   "outputs": [],
   "source": [
    "sampler = Sampler()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "296717d3",
   "metadata": {},
   "source": [
    "## Create dummy structures NaCl and KCl (lattice params are incorrect)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "8c0e1fde",
   "metadata": {},
   "outputs": [],
   "source": [
    "lattice = Lattice.cubic(4)\n",
    "species = [\"Na\", \"Cl\"]\n",
    "coords = [[0, 0, 0], [0.5, 0.5, 0.5]]\n",
    "NaCl = Structure(lattice, species, coords)\n",
    "\n",
    "lattice = Lattice.cubic(5)\n",
    "species = [\"K\", \"Cl\"]\n",
    "coords = [[0, 0, 0], [0.5, 0.5, 0.5]]\n",
    "KCl = Structure(lattice, species, coords)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2a7b1030",
   "metadata": {},
   "source": [
    "### create 20 slightly-perturbed structures for NaCl and KCl, one sevrely-strained KCL to structure list"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "d7d6f26c",
   "metadata": {},
   "outputs": [],
   "source": [
    "structures = []\n",
    "for i in range(20):\n",
    "    struc = NaCl.copy()\n",
    "    struc.perturb(distance=0.1)\n",
    "    structures.append(struc)\n",
    "for i in range(20):\n",
    "    struc = KCl.copy()\n",
    "    struc.perturb(distance=0.1)\n",
    "    structures.append(struc)\n",
    "\n",
    "struc = KCl.copy()\n",
    "struc.apply_strain([0.5, 0.5, 0.5])\n",
    "structures.append(struc)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b965cef6",
   "metadata": {},
   "source": [
    "## Sample 3 structures using the MD sampler with only one line"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "bc1c61f2",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "wcss = 0.07800547033548355\n",
      "returned the [10, 40, 25] structures\n",
      "We have sampled 3 structures\n"
     ]
    }
   ],
   "source": [
    "sampled_structures, wcss = sampler.sample_structures(structures=structures, num_out=3)\n",
    "print(f\"We have sampled {len(sampled_structures)} structures\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "24d4f62e",
   "metadata": {},
   "source": [
    "### WCSS stands for \"Within Cluster Sum of Squares\". \n",
    "###  It's a measure used in k-means clustering and other similar algorithms to evaluate the quality of the clustering. Here, the smaller means the better coverage of selected structures as a statistical representative for the original input distribution.\n",
    "### In practice, you should use WCSS to check whether the number of selected structures can present a good coverage of the input structure list"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "id": "56a3ff58",
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[5.530797481536865, 2.4744789600372314, 0.09148959815502167, 0.04608744755387306, 0.04174330085515976, 0.005896228365600109, 0.0025841835886240005]\n"
     ]
    }
   ],
   "source": [
    "test_num_outs, wcss_list = sampler.get_wcss_list(\n",
    "    structures=structures, test_num_outs=[1, 2, 3, 4, 5, 10, 15]\n",
    ")\n",
    "print(wcss_list)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "id": "73fbfe26",
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "import matplotlib\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "font = {\"family\": \"arial\", \"weight\": \"normal\", \"size\": 18}\n",
    "\n",
    "matplotlib.rc(\"font\", **font)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 35,
   "id": "df75bbf1",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x16028f850>]"
      ]
     },
     "execution_count": 35,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAeUAAAEsCAYAAAD916seAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAA4vklEQVR4nO3deZwcVbn/8c8zkz0zyUDCDEsSAjJBRUhAxF25iqKC4gqoKIiK+wIqv6teMdflulxXBFRwQURFLpu7sososgqyCQkSEgIJSci+Z+b5/XFOZSqV6pmeme6unu7v+/XqV8+cOn3q1Knqemo5ddrcHRERESleS9EVEBERkUBBWUREpE4oKIuIiNQJBWUREZE6oaAsIiJSJxSUK8zMrOg6iEj9G6n7inqsdz3WaajKCspm5mY27GenzGyOmf3ZzNbG1/nDLbNemNlTzewaYEoq7fDYdr8tsGpJXebGuny86LrI0NXTNjVcZvZhM5tvZpvN7AkzO6qgeriZravh/HbaV1S4/KlxmRZUoexnA3+udLlDVe22LMKoWs0oHsn8CpgB3A38C7ilVvOvgd8D+xRdCZGRwMxeCnwb2AZcBWwC5hVaqdoZyfuKm4ANRVciZSS3Za6aBWVgd0JAXgE809231nDetZB31eEW4GlAzY7CRUaIw+L7Oe7+kUJrUnsj+bZhvV0mHsltmauWQXlsfF/SgAE5l7tvIFwREJEdJfuDRYXWQqTeuPuAL8BD1h3Szo/p+wDvB+4hXIJ6DPgusFsq7/VJGelXprzXA9cBa4CNwJ3Ax4AxmXxz4+ePAi6KeZ8A3gUcHqfNBZ4LXAOsJZyd/wLoJBzpfYgQLDcA9wLvzllmA04gXFpbDmyN5VwFHJXKl8xzp2VLTfttTvmDXd7DgeOAW2O9lwM/A/Ypcx0m5Xw8ZznfBdwMrCec1f8deCdgOeUcE9t1Saz3g8A30us7lfdk4K+x3dYTblvMBdrKqO+tsb6HlJh+D9AL7B3/nwKcCdwX65Wsq2PKbJ/0tnMI4bLYqtge1wIvLpE/b91OjdMWlCh/yNtmer7A84C/xLxLgB8C00os32zg/wjflc1xvX0OmNhPO5wBrCRso+eW0YZdcR08HOexFLgYmJ3KM5Oc7wtw/gBltwIfj9vF6th2twGnAqNz8ncSLo8ndXkM+FGyveTs39blpO8L/BhYHMtYEMucWqKOewNnx3luAB4gfDemDLSvGOx6SuU/GrghrqMn4vyfkt3+BmjbveNyziPsw5cSbjW+MJXnpJy6L8hMeydwFuE78yTw2dT6vqfEvNdl22C4bcnwvpu52/wgt6cJwBeAu+LyrSJ8T08mZ5+60+fLXGn9BeUr4vvNwGVxZThwWyrvp4HLY/oq4ELgwtT0b8Zpm4ArY96knBuA8am8c2P6g4Qv5+XAfODQVONeA2yJK/Jywg7LCQHirDjtujivnjjtXZnl+1FMXwn8LjWfZAN4fcz3tLg862L6/yXLVmrjGOLyJu18F3Bp3CgceBQYV8Y6TMr5eCqtNZblhJ3cb+MrWZZfAi2p/K8nBMKNsd6XAo+k1seEVN5TY/rq2H5XEL7sTthA+904gY/GvF/Omfb0OO3G+P84woGEEwLapYRAmqzbnQ66cspM1tUf43pZGNfL3TF9C6kDhFLrtowv/rC2zVQ5C2I9HyFsc/+K6Y8DMzP1OYawI3FCILs0ft6BO4DJOeU/GOf/h9i2Hxug/fZPLcu8WKdb4v+b6fu+TCV8X/4Zp90e/z9lgPK/HfMvBX4N/CZuWw78LJN3P8L3wgnf2cvifJxwEDQnZ/+2LpP2HMK+ygkHgJcC98f/HwZmZPIfSjhQdsLB1CWxDZNtsoN+9hWDXU8x/4fo2zb/SPiebSQcmO6w/fXTrlOBh2L+f8Q6/S3+vw14Rcz3wlh3J5ygXAh8M047KbXNbCIE9H8Cb2QIQXm4bcnQv5u52zyD2J4IB9ZXxGmPEL7jf4jt4sAXB1wnA2VIbbTZhjufvsByZCq9K7UA6SOt3JVD2NEnCzArlT6ZsHPyZOXH9Lmp+e4f01oyjeuEL7Gl6rQ29blnp8p7d0z/eyrteTHtn8CkVHoL8NU47arMciyI6VNTaTttHMNY3l7gHan0NsKXyIG3lbEOk3LSQfk0+r6Mu6fS9yLsiBz4SCr9IcIXcv9U2mjgTzHvyTFtLOHMeBnQmcrbTjiocOAlA9R3d8JO4aGcaf8dy3hf/P9t8f+fZPIdHtttYRntk952vk/qigVwXkz/YX/rtswv/nC3zXQ5lwFjU9vmWTH9V5l2XE3YWae/p6OA78X855Uo/8T0tt9P2xnhSo8DX2THA7k3xG1mHbBXf9tjP+XPiHnvJ3XGCOxBOIt14Cmp9Jtj2mcydXk7fTv21lS6kwrKhO13Ydx2Tsos56fJfP8JB7fJwdvHMuk/jOlfG2BfMdj1NDNuL6uBg1Pp+xMOzHbY/vpp28/EvJ/LpCdtdUMmfYe2imknpbaZw9PbDIMMyhVqy8MZ+ndzp21+MNsT4eDFCScF6W3s6XFZN5E6ecltk4FWWmpFeCbt/Jj+g5z838pp1NyVQ+he78Crc8rZPS7EBuKXkb4v829y8ieNu57U2Wac9rs47buZ9F1j+vJU2pGEHd7ROfNILg09kEkva+MYxvJenZP/o3Had8pYh0k56aCcnOUelJP/4DhtYSptQ3ztmsn7DEIAeWr8vzP5LDFopPI+H3gHmTONEnW+KpbzzEz6/YQd/dT4/+kx31dyyngX4Yi9dYB5Jetqbc62MydOu7u/dZua1t8Xf7jbZlLORmCPTP4xhB1yTzKNviDyuZx6jiOc3W4GdsmUv5p+AnGmnP+g7+wu75bHt7N1yNse+yn/sJj3pmydgFcSDsqSbeH55ASTVP4r4vT0LagdAg3w1ph2QYky/hGnHxD/f1GpeRJuq8wHfpRKW8DO+4rBrqfPxvz/lZP/pOz210/bnkPqADeVbsCHgVdl0ndoq8z87s4pfyaDC8qVaMtkGx7sd3OnbX6w2xNwbPz/lzl5jwOOp8StiORViZ5rN+ekLYnvE/v7oJmNIlwm2kg4xd+Buy8hXs4FnpWZfE8/Rd/j7hszacvj+52Z9JXxfVxqvn9y99e7+/ZnQc1snJkdQtjBQ9gBDsowl3fI7VyiLjMIZyDz3P2fOXX5B+Ey5HQzmxmTk7rdbmafNrM5Me897n6eu/8r/v8E4ehxOnCbmZ1qZrPitL+6+4/dfWEZ1fxZfD82Ve/ZwFOBK909Wac3xPfTzewiM3uTmU2O8/uBu1/i7j3ltAvwz5xtZ8jtnGNY22bKn9398XSCu28hHMi0EK72ALw4vl+XLcDdNxEOEscQ7nOn3efuvSWWIeuF8f0yj3ufjF/G9xeVWV7WPYTbO88B/mJm743bL+7+B3f/aWpbKLm80Z8y+fIMVMaVmXzJcv0um9HdV7j7fu5+cj/z63eeJdZT0uZ/yuYnXD4uV/LdOdPMzjWzo8xsggdnuvvvB1FWf/vkclWiLYcqb5sf7PZ0E+GE4Vgz+4OZvd3MOgHc/ZfufpG7r++vEpUIyqty0raVWf4Uwoa22N23lcizIL7vnklfSWl505KdxYodEvN3IpjZBDP7kJldaWYLCWc4twNfSrL0M/9ShrO8q3LyltvOefaI74/0kydbl1MIZwkzCR0Z/mFmi83sbDN7euazxwP/JpxFfwN4wMweMrOvJjvUMlxGuHLwplRaEqCTgI27/53QCWgL4Wj0YmC5mV1nZu82s9Flzg8q385Zw942owUl0hfH9z3j+7T4fm0yCFD6RV977lVGPUsZaFtK6prdpsvi4SmGYwkdmZ5H6Ej6iJndbWZnmNnUVPZkec8osbznxOnZ5U1LyvhRiTJOz5SRLP9wepIPdj0l81xMhruvJHRUGpC7X0To49JCuNr1W+BJM/udmR07yJGyBrPNlFKJthyqvPoPanty90WEKwfrgFcAPwGWmNnNZnaamQ14YF+JR6L623EMJFnh/ZWR7Ag3Z9L7O4ovFfDKYmZ7EToj7UPYUd5C6Ol9B+Fo8O6hFh3fh7K8w2nnitTF3Rea2TMJl3teR7jMP4vQ+/4UMzve3S+Nee8ys6fGPMcALyf0Zv0E8AEzO8Ldb+qvgu6+Jo5c9UYzO9TdbyPsnDaQORtw96+b2YWEKxlHEY64D4+vk83scHfPtmnubMvIM5DWfqYNa9tM2VQiPVmvyWOHSV0uI1yhKSV75aLcs+T0PEu1Xaltumzufk28YvMa4NXAEYQDvmcAHzSz57n7fPqW96+UPnCBcKm9lKSMKwn9Ikq5P5N/ONvOUNdTKWVvZ+5+mpmdSTj4fRXhwOdV8XU8oR9MOQazzSSy35VKtOVg5peWV/9Bb0/u/nMz+z1hH3k08BLCLZjDgPeZ2XPcfUV+MbV9TjnPCsLZzTQzG1Xi7DEZreWJ2lWLL8T5nkPo6LS9XmbWPYxy62l5H8vML89OdYlnb9fFF3FH+Z/Ae4D/IfQYTfJupa9HN2Z2APB5wsY6lxCwB/IzQqB9g5n1EnpC/iLvEpC7LyU8RnG2mY2J5X+PcNnz9YRHjyol+QLnfck7KjifUvYokb53fH80vj9OOHD6UjyoqYaBtqWKbNPxsv8vgV/GM7hDCVdhXkA42HsPYXkBrnD3rw1xVkkZ57n7JWXkT25vTMubaGZvJhxE/dZLj9Ew2PW0mNB5aG/62j+Z33jCNri2jHIAcPcFwP8C/xvP5l5H2P+9zsyeO9ABdD9Kfk/MbCw735qpRFtW8rs5pO3J3VcRHjP7sZm1Ek4SziKss1Pou+K6k0JHQ4mNejPhPuUrs9PNrIvwhVtPOEutlWfH96/kBM6Xxfds2w14ZFdPyxvv6S4E9jOzg3LqcjDhzPYhd19sZnub2Z2WGXM5fpk/RPgiTI+ffbaZ3Wdm383kvZew8yTJW4bkeeGj6Tti/1k6g5l9zsweM7MXpOa1xd1/Q3i0bTDzK1cySltnzrTnVHheeV6YvSwfd6YvJ5wlJzvRG+P7Tttb/MyVZnajmT1jGHVJ5vG6Epc7k34YN+RMG5CZHWNm88zsk0lavOd5K+HgDvrWb1KXV5Qo64tmdpuZvaWfWQ7UZj8xs7+b2Uti0t/i+04HmXGd/IiwQ072JXn7isGup2vj+2tzsh9J+b9r8EMzW2Zm278f7r7e3S+k72rUcL47yfdkt5xtI+97Uom2rOR3c1Dbk4X+DgvM7K1JHnfvcffrCAeQMEB71sMQZWfG97PSZ6FmNomw8x0D/Dinc0w1JWcZr04nmtkRhLNB2PkIL7mcOHmAsutpeb8T338SDwiSuuxB6F0P4f4dhAA+BXilme3QLoT7uC2EgR0gPCe5D/B2M8t2WHtzfL+VMsTOS5cQLlOeRLjacGUm20LCmePn0/dszKydEMzLnt8gPEC46nFw+mDAzPYh9Iyttj2Ab5tZS5zvaMKjXLsSntdMLo+dR9g2P2lm23csFpxBOMicRt+l2KG4jnBL55mEdbB9v2JmxxBub6wHLhhi+fcRrpB8xMz2TZXdQtj2oG/9XkN4tvWlZvZfmbocSRig52D63x4uIpzVv8PMTkxPMLN3EHp7H0DfgfPVhGdcX2pmp6TythLu144DLkr1EcjbVwx2Pf2EcN/4I6mDg6QD5zco3+OEHslfTh/kmdnu9D1SeHsq/yZgXLwSNSB3f5JwVj8F2H4gZKEfwNdzPlKJtqzkd3Ow29N8wtWLT1uqr0NsrzfEf/vfF/XXNdv7unI7pR+JemNO/o/HaXNTaTMp0TWeEByc0MB/ItxXWRHTbmDHQSnmUuJRCvrvCt9ffXfo5k+4J+nx9TfCg+nJ87WPEToEbCX1uA9hQAMnfGkuIfTUza1PBZf3jXHa+WWsw53KIVzeuYK+R4GSQRmS52YvZsfn8l5F34AWN8d2SZ7hWw8cmsr73pjeE5fpEvoGjFhKZoCLAeqetKOTeWwoTh8d5+GEe4C/ia9kAIKLBzGPsh6jiOlnxvQthN70ySAB1xEO7BaUWf5gts2knNsIZwsPxPXw75j+T6AjU8ZbY14nBJJL4+eS9X5YOfUcoP0OoG9wmHlx20kGdNkIvGmg7XGA8r9M33fm6rgM82LaA8RHhWLeAwlB1QkDfVxOuHKQbEPv76+NY9pL4jbthIOCy+h7Fnsr8JpM/kMJ+4WkjS+hb1COu9hxQKCd9hWDXU8x/7GxLj2EM+dfEc4S747vC8po18mpeSwi7A/+QN8+4KuZ/Ml3+DbgpzHtpJh2Vol5JOMh9MZ6/ppwQHFvXE6vQltW5Ls5xO3pFzF9TZzv5fSN3XEjmVEbd5pfmV8Iz2m486lQUI7Tj6NvuLh1caV/CBhV7pe5v8YdoL55X8qjCTf3VxIun95FHE6ScF/LSQ3hCHTHBk+GeDxwgPpUYnmHFZRjeishgN5K2AmtInRyyx2QhNDB5o9xGbcSjoJ/SmpAkcxO48+E5/+2EDpKfBfYs5ztLlWOEXYYTmpAmkyeSYQd978IHYrWxWX6AAM8o1zGtlMqKLcSLscn83yEcCVlXFzWsr74g9k2U+V8i3Cp+jbCzubRuH1OLrF8hxF28snwjfMJlwL3K7cdymjDPQn39B+J81hMODs+sNztsZ+yW+J2egt9gzA8CHyFzEFIzL8X4TLnw6m6/I7U4Bb9ff9j+v6EM9JkmM1HCAcbB5eo475xXT5O+G4sJJzdTc7k22lfMdj1lMr/AsL3cRVhX/VzwoA0q7Lbaz9tuzvh/vHDhO/pasK+6S05eZ9DODjZHOu4CwME5fi5d8bPbYrtc0787PVkYksl2pIKfTeHuD2NJfSzuSvWawPhQOnTlDH6YjKqkIiIiBSsHu4pi4iICArKIiIidUNBWUREpE4oKIuIiNSJokf0kn5MnTrVZ86cWXQ1RKTB3H777cvdfbei6yE7U1CuYzNnzuS226o1MqKINCsz6++HaKRAunwtIiJSJxSURURE6oSCsoiISJ1QUBYREakTCsoiIiJ1QkFZRESkTigoi4iI1AkF5QZy5b1LOO+GfxddDRERGSIF5Qby5weX8Z1r56Gf4xQRGZkUlBvIrK521mzaxhNrNxddFRERGQIF5QbS3dkGwLyl6wquiYiIDIWCcgPp7moHYN4TawuuiYiIDIWCcgOZ2jaGjgmjeVBnyiIiI5KCcgMxM2Z1tjNfZ8oiIiOSgnKD2a+rjQeXrlMPbBGREUhBucF0d7axeuNWlq1TD2wRkZFGQbnBzEo6e+m+sojIiKOg3GD6HovSfWURkZFGQbnB7NY+lsnjRzPvCZ0pi4iMNArKDcbM6O5s0+VrEZERSEG5AXV3tfPgE2vVA1tEZIRRUG5A3Z1trNqwleXrthRdFRERGQQF5QY0S8NtioiMSArKDai7Sz9MISIyEikoN6DO9rG0jxulM2URkRFGQbkBmRmzutp1piwiMsIoKDeo7s42PassIjLCKCg3qO6udp5cv4UVGgNbRGTEUFBuUMlwm/ptZRGRkUNBuUElj0Xpt5VFREYOBeUG1TVpLO1jR+lMWURkBFFQblBmxn5dbXosSkRkBFFQbmCzOvVYlIjISKKg3MC6u9pYoR7YIiIjhoJyA+ve3tlLZ8siIiOBgnID2/5YlIKyiMiIoKDcwPaYPI62saOYv1SdvURERgIF5QZmZuzX2abHokRERggF5QY3q0tjYIuIjBQKyg2uu7Od5es2s3L9lqKrIiIiA1BQbnD7dYXOXjpbFhGpfwrKDS4ZA1sje4mI1D8F5Qa35+RxTBzTqpG9RERGAAXlBhfGwG7XmbKIyAigoNwEuvVYlIjIiKCg3ARmdbWxbO1mVm1QD2wRkXqmoNwEujuTzl46WxYRqWcKyk1gvzgGtjp7iYjUNwXlJrBXx3gmjGlVZy8RkTqnoNwEWlrCGNg6UxYRqW8Kyk2iu1OPRYmI1DsF5SbR3dXG0jWbWb1xa9FVERGREhSUm8SsOAb2fJ0ti4jULQXlJpE8FqVBRERE6peCcpPYq2M840drDGwRkXqmoNwktvfA1uVrEZG6paDcRLr1WJSISF1TUK4xMzvIzLaY2dxaz7u7q50lazaxZpN6YIuI1CMF5Roys1HA+cDoIubfreE2RUTqmoJybX0SOKComc/qCj2w9ViUiEh9UlCuETM7EPgv4PNF1WHaLuMZN7pFj0WJiNQpBeUaiJetfwxcDVxYVD36emArKIuI1CMF5dr4f0A38J6iK9Ld2c68pbp8LSJSjxSUq8zMDgDOAD7u7o+Wkf8UM7vNzG5btmxZxevT3dXG46s3sVY9sEVE6o6CchWZWSvhsvWN7n5eOZ9x93Pd/VB3P3S33XareJ2S4Tbn6xK2iEjdUVCurk8As4H/NLOpZjYV2CVOmxDTaroO9FiUiEj9UlCurlcAY4BbgGXxdUec9on4/4xaVmj6rhMYO6pFw22KiNShUUVXoMF9jL4z40QXoQf2T4ELgCW1rFBri/GU3dr0WJSISB1SUK4id789m2ZmM+Of/3b3q2tbo2BWVxu3LlhZxKxFRKQfunzdhLq72lm8aiPrNm8ruioiIpKioNyEks5e6oEtIlJfdPm6xtx9AWBF1qE7joE9b+la5kzvKLIqIiKSojPlJjRj1wmMGdWi4TZFROqMgnITSnpga7hNEZH6oqDcpLo79ViUiEi9UVBuUrO62li8aiPr1QNbRKRuKCg3qf00BraISN1RUG5Ss7riGNgKyiIidUNBuUnN2HUCY1o1BraISD1RUG5So1pb2He3ifq1KBGROqKg3MS6u9p1piwiUkcUlJtYd2cbi57cyIYt6oEtIlIPFJSbWNLZ66En1hdcExERAQXlAZnZWDM72Mz2L7oulZY8FvWgRvYSEakLCsqAme1iZl8ys+sy6ccDjwK3AfeZ2X1m9pxCKlkFM6dMYHSr6bEoEZE60fS/EmVmk4CbgFnAo2Y2yt23mdmhwIWEA5fvAfcAJwJXm9kcd59fWKUrZFRrC/tObWO+OnuJiNQFnSnD6cC+wLHuPsPdk15PnyH8xOKZ7v5+dz8HeDGwGPhkMVWtvO4ujYEtIlIvFJThtcAF7n5JkmBmE4Aj479nJ+nuvolw9vyyWlawmro721m0cgMbt/QUXRURkaanoAwzCfeM014IjAEedvd5mWmLgM4a1KsmurvacIeHlulsWUSkaArK0AO0ZtKOiO9X5eTfDVhT1RrVUN8Y2LqvLCJSNAVluA94bvKPmRnwBsCBX+fkPwa4vzZVq769p0xkdKvpvrKISB1QUIafAMea2UfM7OnANwiXtP8N/Cmd0cw+RQjgF9e6ktUyurWFfaZqDGwRkXrQ9I9EAd8HXgB8k3B2bMBK4C3u3gtgZicTelzvC9wIfLeYqlZHd2c79zy2uuhqiIg0vaY/U/bgBOD5hMej3gvs7+63prJNIwTrLwBHJsG6UXR3tbHwyQ1s2qoe2CIiRdKZcuTuNxEGEcmb9jngc7WtUe10d7Zv74F9wJ6Ti66OiEjTavoz5XLEYTjHFV2PatneA1v3lUVECqWgDJjZaDN7j5n9OJP+QjO7F1gOrDOzP5jZU4qpZfXsPWUio1pMj0WJiBSs6S9fm9lY4GrCPeUtZvbuOPZ1N6H39Vjgj8C9wBuBv5nZQe6+tLBKV9iYUS3MnDpRj0WJiBRMZ8rwEeB5hE5eu6TGvp4LjAN+4e5HufvpwKHANuBTRVS0mmZ1tTFfvxYlIlIoBWU4FrjE3b/m7hsBzGwM8BrCI1JfSzK6+5PA+cDRBdSzqvbrbOeRFevVA1tEpEAKytAN/CWT9lxgIvC4u9+ZmTYf2LMG9aqpWV1t9Dr8e9n6oqsiItK0FJTDuNfbMmkvje9X5+SfDGyoao0K0N3ZDmgMbBGRIikohzPfOZm01xEuXf82J/+R8TMNZZ+pE2ltMT0WJSJSIAVluAh4u5kdY2YTzOxU4ADgCTI/SGFmbwVeDvyq9tWsrjGjWpg5ZYLOlEVECtT0j0QRfoDiSOBy+sa+3gKc7O5bAMzsdcAHgcOBBwjjZDec7s52HlyqoCwiUpSmD8ruvsXMjgCOI3TwWgtc6O73pbI9k/Ac84XAaUkv7UYzq6uNK+9bwuZtPYwdlf2JaRERqbamD8oA7t4D/Dy+8vwPcEaj/RBF1n5d7dt7YD9tj0lFV0dEpOnonnJ52giXtRva9jGwNYiIiEghFJQjM/ugmd1tZnlXD74FPBY7gTWspAf2fN1XFhEpRNMHZQsuAM4kDAqyd062fwO9wNfM7Be1rF8tjR3Vyt5TJmgMbBGRgjR9UAbeBZwAnAPs5e4PZTO4+38B+wA/BY41s7fXtoq1093ZpseiREQKoqAcgvIN7v5Bd99UKlOcdjJwF/CeWlWu1mZ1tbNgxQY2b9MY2CIitaagHAYKKWswkNj7+hLgoKrWqED7dbbR0+ssWN5wI4mKiNQ9BeUw7nXJM+Qcywn3lxtSMga2BhEREak9BWWYR/id5HI9C1hYpboUbt/dJtJieixKRKQICsph7Ou3mtkBA2WMed4K/L7qtSrIuNGt7D1lIvPV2UtEpOYUlOH7wCPA9Wb2VjPbaXxJM2sxszcDVxGG4fxWbatYW92dbXosSkSkAE0flN19HfAaYAVwAfCkmV1rZj8zs4vM7FpgJWHc683Ake7+eHE1rr5ZXe0sWL6eTVvVA1tEpJY09jXg7g+Y2WzgA8DxwAvoa5stwE3AZcC57r65mFrWzoHTJrOt17nv8TUcMmOXoqsjItI0mj4om9nPgGuB6939G4SfcsTMpgI97r6yyPoVYc70DgDuXLhKQVlEpIaaPigTfrLxeAAzewy4Lnm5+4IC61WYrknj2H3SOO56dFXRVRERaSoKyjAZeDbwvPg6ijDsppvZo4QAfT0hSD9SVCVrbc70Du5atKroaoiINJWmD8ruvp5w+fraJM3MngY8H3hufL0tpi+k7yz6p7Wvbe3Mnt7BH+9dwsr1W9hl4piiqyMi0hSavvd1Hne/391/4O7vdPenA53ARwg9tE8Czi+wejWR3FfWJWwRkdpp+jPlPGY2mnAp+3BCT+xnAxOBHkJP7OsKq1yNHDhtMmZw56JVHL5/Z9HVERFpCgrKhN9UBp4JvBR4CeHS9XjCGNf/AL5LuLx9Y7zc3fDaxo5iVme77iuLiNRQ0wdlM7uMcEY8OSb9kzDK13WEn3RcU1DVCjd7+mSuum8p7k44bhERkWrSPWV4LTAB+Akwy90PdvePuftvmzkgQ+jstXLDVhY9ubHoqoiINAUFZfgpsJTQgetfZnabmX3VzF5hZm2VmIGZHWlmfzGzDWa2zsyuNrPnVKLsato+iIg6e4mI1ETTB2V3P9Hd9wZmAR8E5gNvJ/wS1JNm9jcz+4KZvdTMxg22fDN7MfAHoAP4NPDfwFOAP5vZYRVajKrYv6udcaNbuHPhqqKrIiLSFJr+nnLC3ecTAvL3AczsQEKnr/8A3g98CthsZjcD17j758ss+lvAIuDZ7r4hln0BcD/wReBlFVyMihrV2sKBe03WY1EiIjXS9GfKpbj73e7+bXd/LbAH4ez5DuBFwNxyyjCzXYDZwMVJQI5lLwX+THjsqq7NntbBPYtXs7Wnt+iqiIg0PJ0p5zCzpxCeTU5es4ExwAbCbyr/ucyi1gD7A3mPUU0Ftg27slU2Z0YHP7jxYR5YspZn7DV54A+IiMiQNX1Qjmezh9EXgA8DdgUMWA3cCJxBCMS3u3vZPzIc887LmedBhGeh/zTc+lfb7GkdAPxj0SoFZRGRKmv6oEwYOtMJQXg58BdCAL4BuNPdvZIziz26L4j/fjln+inAKQAzZsyo5KyHZNou45kycQx3LVrF256zd9HVERFpaArK8H+EX4G6wd3vreaMzGwC8GvC5fAvuftOl8Hd/VzgXIBDDz20ogcEQ2FmzJnewZ0a2UtEpOqaPii7+3G1mI+ZdQC/JVy2/hHh8agRYfb0Dq594AnWbNrKpHGji66OiEjDUu/rGjCzTsKwnc8nnAW/q9KXxatpzvQO3OGeR1cXXRURkYamoFxlZtZO6NA1B/imu79nJAVk2LGzl4iIVI+CcvWdTQjI33b30wquy5BMnjCafadO1C9GiYhUWdPfU64mM3sa8DbCo1V3mtkJ2TzufmHNKzYEs6d38Nf5y/WLUSIiVaSgXF0vju+TgR+XyDMygvK0yVz+j8UsWbOJPSaPL7o6IiINSZevq8jdv+fu1t+r6DqWa86MXQD04xQiIlWkoCxledoe7YxpbdHPOIqIVJGCspRl7KhWnrbnJHX2EhGpIgVlKducaZO5+9HV9PSOqCe6RERGDAVlKducGR2s39LD/CfWFV0VEZGGpKAsZUsGEblz0cpiKyIi0qAUlKVs+0ydyKRxo7hzkYbbFBGpBgVlKZuZMVu/GCUiUjUKyjIoc6Z38ODStWzYsq3oqoiINBwFZRmUOdM76Ol17lm8puiqiIg0HAVlGZTZ0zsA9LyyiEgVKCjLoExtG8u0XcZrZC8RkSpQUJZBmz29Q2Ngi4hUgYKyDNrB0ztYvGojy9ZuLroqIiINRUFZBk33lUVEqkNBWQbtGXtOprXFuEv3lUVEKkpBWQZt/JhW9u9q1yAiIiIVpqAsQzJ7egd3LVpFr34xSkSkYhSUZUgOnt7Bmk3bWLBifdFVERFpGArKMiRJZy9dwhYRqRwFZRmS/TrbmDimVT2wRUQqSEFZhqS1xThw2mSdKYuIVJCCsgzZ7Okd3Pf4GjZv6ym6KiIiDUFBWYbs4OkdbO1x7ntMvxglIlIJCsoyZBrZS0SkshSUZcj2mDyerkljdV9ZRKRCFJRlWGZP6+CuR1cXXQ0RkYagoCzDMmdGBw8vX8+qDVuKroqIyIinoCzDMmdaB4DOlkVEKkBBWYblwGmTMVNnLxGRSlBQlmFpHzea/XZrU2cvEZEKUFCWYUt+McpdvxglIjIcCsoybHOmd7Bi/RYeXbmx6KqIiIxoCsoybHP0i1EiIhWhoCzDtv/u7Ywd1aLOXiIiw6SgLMM2urWFZ+ylX4wSERkuBWWpiNnTOrjnsdVs7ektuioiIiOWgrJUxJwZHWza2ssDS9YWXRURkRFLQVkqom9kr1WF1kNEZCRTUJaKmL7reHadOIY7F64quioiIiOWgrJUhJkxe9pknSmLiAyDgrJUzJzpuzDviXWs3bS16KqIiIxICspSMbOnT8Yd7l6sX4wSERkKBWWpmNlJZ69FCsoiIkOhoCwVs8vEMcycMoE7F60suioiIiOSgrJUVPjFKJ0pi4gMhYKyVNSc6R0sWbOJJas3FV0VEZERR0FZKmq2fjFKRGTIFJSlop6+xyRGt5qCsojIECgoS0WNG93K0/aYpJ9xFBEZAgVlqbjZ0zq4e/Fqenq96KqIiIwoCspScXOmd7Bu8zYeWrau6KqIiIwoCspScersJSIyNArKUnH7Tp1I+7hRCsoiIoOkoCwV19JizJ7Woc5eIiKDpKAsVTFnegf/WrKWjVt6iq6KiMiIoaAsVTF7egc9vc69j2nITRGRcikoS1XMnj4ZUGcvEZHBUFCuATPbx8wuM7Mn4+sCM9ut6HpVU2f7OPbqGK+gLCIyCKOKrkCjM7MpwHXAGOArhDb/BHCQmR3m7luKrF81zZ4+mbseXVV0NURERgwF5eo7DZgGHOju9wOY2c3AVcCJwHkF1q2q5kzv4Pd3L+GQz1/F1LYxTG0bu/01pW0Mu7WNZWr7mB3Sxo5qLbraIiKFUVCuvuOB65OADODuV5vZA3FawwblYw+dzqatvSxZs4kV6zazfN0W7np0FcvXbmZ9iV7Zk8aNYmp7CNK7tY3dHsx3bRvD6JYWzKDFjJaW+L79BRbfW1tCWpK3taXv7+T/1hZjVMw3qjWmxWmjWvv+7svbQksL4T3OS0Sk0hSUq8jMdgH2BS7JmXwHcFRta1RbHRPG8OGXdudO27ilh+XrNrNs3WaWrw0Be/m6zX2vtVu4f8kalq/dzJpN22pc84FtD9gWgntrbiA3WuJ7a0sLrS2EdwvBvVTevs/E8lvjwUOmnJ3y9nOgkT64yPv8UA5OWi0s9/bPx/JEZOgUlKtrr/i+OGfa48AkM5vs7k333ND4Ma1M33UC03edMGDezdt6WLVhK9t6nd5exx163Ol1x93pdeh1pydO602l9fam/nantzd8tqe3l55etr9v6+2lp9e3v7b1hvzbeuJ7Zlr28z29vTt8pqfX6Ymf6+3d8fPp1+ZtPfR4/HxqXnmf2SHd+9LriRnbA3lytWL7tB3y2c7p5eQdTHmZ/OkppctJpw9uvjvMyXb+e7Dllah6bv63HDaDk1+wT25dZGRRUK6u9vi+IWfaxvg+EdgelM3sFOAUgBkzZlS1ciPF2FGtdE3SveY8yUFJ9qBixwOIvoOEnt7sgUbfQUnvAAcn2z+/w0FD/sFJTzwg6qtnqs471D9Jy8+bXdb+yuivnFL5KZV/GGVuz18ybznLMXD+9D9T28cijUFBubqSR876O53pTf/j7ucC5wIceuih9XUaJHXHzGg1aG3RQYtII9BzytW1Nr6Pz5k2PpNHRESanIJydS2M73vkTNsTWOXu62tYHxERqWMKylXk7quAh4FDciYfDNxW0wqJiEhdU1CuvkuBI8zsqUmCmR0B7A9cVFitRESk7qijV/V9FXg7cI2ZfR0YB5wO3A5cWGTFRESkvuhMucrcfRnwIuAu4HPAR4ErgFe6++biaiYiIvVGZ8o14O4PAK8quh4iIlLfdKYsIiJSJ8xLDZ8jhTOzZcAjRdejAU0FlhddiSagdq6dwbb13u7e0L/pPlIpKEvTMbPb3P3QouvR6NTOtaO2bhy6fC0iIlInFJRFRETqhIKyNKNzi65Ak1A7147aukHonrKIiEid0JmyiIhInVBQFhERqRMKytIUzOwWM/Oc1yVF161RmNl5ZnZ9Tvo+ZnaZmT0ZXxeYmZ6RHYZ+2lrb+QinYTal4ZmZAU8jjDl+aWayBmepADN7J/Au4M+Z9CnAdcAY4CuEfc4ngIPM7DB331Lruo50/bS1tvMGoKAszWAm0Ab8yt31y1wVZGatwKeBuSWynAZMAw509/vjZ24GrgJOBM6rQTUbQhltPRNt5yOeLl9LMzggvt9faC0ajJmNA+4A/hv4KbA4J9vxwPVJQAZw96uBB+I0KUOZba3tvAEoKEsz2GFnZWYTC6xLIxkHTAKOc/cTgW3piWa2C7Av4bfDs+4ANCxk+fpt60jbeQNQUJZm8AxgLfANM1sLrDOzh8xMZ2rDswbodveLS0zfK77nndU9Dkwys8lVqVnjGaitQdt5Q1BQlmZwANAOdABvB04m7Lx+YWZvK7BeI5q797p73hlboj2+b8iZtjG+62yuDGW0NWg7bwjq6CXN4Fyg1d3PThLM7CLgHuB/zezn7t5TWO0aV3LQ39+wgb21qEiT0HbeAHSmLA3P3b+X3lHFtI2EDjNdwNMLqVjjWxvfx+dMG5/JI8Ok7bwxKChLM3sivrcVWovGtTC+75EzbU9glbuvr2F9mpW28xFEQVkampntZWb3mtkZOZOfGt8frmWdmoW7ryK07SE5kw8GbqtphRqYtvPGoaAsDc3dFwOTgXeb2aQk3cymAycB17n7koKq1wwuBY4wsyQwYGZHAPsDFxVWqwaj7bxx6KcbpeGZ2WuBy4F7CSNItQMfJAz9+Pz0wBYydGa2AFjg7oen0nYjdDTaBnyd8Lzt6cB8Qttvrn1NR74Sbf1atJ2PeDpTlobn7lcArwXWE8Zf/hhwE/A87aiqy92XAS8C7gI+B3yUMDbzKxWQK0vbeWPQmbKIiEid0JmyiIhInVBQFhERqRMKyiIiInVCQVlERKROKCiLiIjUCQVlERGROqGgLCIiUicUlKXumJmb2flF16M/ZjbWzH5kZmvi69U1mOcCM7u+iuV3mlnFft/YzK6PI09VoqwxZrZXJcoa5HxbzGxmrecrzUtBWWRo3g28gzA61WmM8B9XMLNXAg8AuxVdlywz2xu4G3hZjec7Cfg7YexokZoYVXQFREaog+L7B9y9EX4T+NlAR9GVKGEfYFYB890VeBbw+wLmLU1KZ8oiQzMGoEECsojUCQVlAbbfr/yemZ0Qf5d1k5nNM7MP5OS7vsTnr8/8f5aZvcvMHjSzjWZ2q5kdZma7m9nFZrbWzBab2RfNbKdt0cw+FadvMLNrzexZOXmONrO/xTwrzexSM5uVyeNm9nkz+7WZbTaz+8ys5FUiMzsmlrnRzFbFzx2ULg84MVX2Tu2Ryjsj1unx2Kb3mdnp2eU1s+ea2VWxTdaa2ZVmdlipcgf7OTN7tpn9PrbRCjP7nZkdGKedD3w2Zn04sx6fbmaXx3bYYGZ/NbMjc8o/IrbZejN7yMyOG6ju8XNjzexbZvbvuG4WmdnZZrZLnH4ScF3M/uPY9pjZ3NierzOzJXHZ3xnTPXsfOC/dzCaZ2TfNbGFctrvN7F1x2uH0/f7wZ5PPmtlJ8e/DM+XvkJ76/w1m9nAsf26cNs7MvhDTt8Rl/5yZjSmnzaTBubteegEsAB4BngT+G3gfcCfgwKsy+a4v8fnrM/8vAh4j/FTffwJrY9o9hN/SPQX4U5zHianPOrAOWAF8Jn72ccKv3xyQyncS0AtcRfiJus/E+a0EZuWUdzXwHuDUftrhAzH/rcCpwBnA0vj5Z8U8JwA3xHwnAC8rUdZo4L7YFp8g3Ie+PH7uU6l8LwO2ALcQfkXp/wH/AjYBL+ynjcv93AuBzfHzn4rLtSC270zgucBlsV4fTZYHOBBYTbjXfHr83C1AD3BcqvwjgK2Enww8Ffif2F5rCD8v2N92d15cr/8DvBP4Zizryjh9X+CLsW7fB06I6XNjPVbE5f4s4Tea58a8MzPz2SGdcKXjTsJPSp4Tt4tfxTwfBrpiW3hsmxOAiYRtzoHDM+XvkJ76fw3whbj+XwC0ErbDjcC3gfcCP4jL8hvijwTp1byvwiugV3284k66FzgolbZ7TPtZJt/1JT5/fU55B6bSvhp3VBel0iYSAkZ6Hh53lunP7kcIQJfG/ycRAsYvMvXYnXBgcXmmvFVAxwBtMCUGiJuBMan0mUl6Ku18wAco77A47zem0gz4I/CT+H8L8BBwI9CaaZd5wD/y2niQn7uZcLAyJZU2KwaCr8b/55IJZsD1hN89nphKG0U4IFmStBGhk9tCYFIq33/E8hYM0EYbgLMyaZ8nHBS1xf8Pj2WdlMqT1Hdu5rM7LUdeOuGg04G3ZNbNDYQDwNa43neYB4MPyueXyHdkJv2UmH5MLb7vetXvS5evJe0Bd/9n8o+7LyGcJe4+xPIecve7U/8/GN8vT81jPfAEsEfms39Mf9bd5wN/AI40s1bCWeIk4Aozm5q8CMH82pgvfYn6ZndfNUB9XwpMAL7u7ltS814A/BQ4zMyy9ezPY8SzYjM70szGePAKdz8x5jmYcDZ4BbBLajnGE86c5pjZtJyyy/qcmXUSOiv93N1XpJbpQeBQwu/u7sTMpgAvJnRyGp8qv4Ow/rqAZ8Xyn0k4OFqTKv864J87FbyzR4Hj4uXejvjZz7j7s9x9XRmf/1MZefIcDSwDfpEkuLsDbyNcWegdYrlZ2fq9Ic739sx2+3vCQdLRFZqvjFDqfS1py3LSNhPOGoZiaeb/bfH9iUx6Dzv3b/hXTnkPAa8hPLbzlJh2UT/z341w1pM3zzz7xPcHcqYlPxK/d6rMfrn7o2Z2OvAlwtnxOjO7BvglcLG799C3HP8bX3mmE4JXWrmf20Y4A5yXU79/9FP9pPwPxVeeGYSrFxDWTda/CL26+/M+4GLgx8B5ZnYTIej/yN1XD/BZKG+95plJOGjc4Qfl3f2R5G8zG2LRO8jW7ymE7TLvuwahTaWJKShL2nDODvIC97acNAhnjwPJy5ME7p7U/E6hr0NO1srU3z1lzLO/vXAy7y395NmJu3/NzH4OvB54FfBy4Bjg7cAr6VuOzxCeic2Td4BS7uf2j38Pdt0m5Z9NOBvPcy+QDOgxLmf6gFfi3P0aM5sBvJpwlvhy4BvAqWb2THcvFbwS5axX2Hn7bKVyZ8N55Sey9WslHCC9v0T+lSXSpUkoKMtg9QBj0wnxMvFU8s+WhmpmTlo34T7ycsL9VYBl7n51pj6HE3Z+mwc5z6TMpwJ3ZaYlwS17xlqSme0KzAb+5u5nAWdZGDHrfOCNsfdzMs91OcvxLMKzshv7qetAn1sYk/fLqd9XgJXu/uV+yt+WU/7TCVcVNsR8Tv5zxPvmpKXLGQvMAR5194uAi2Kv9NMIZ//HA9/pr4wcSRAcm0nP3oJZSN+z5uk6vTLO9/Rhll/KAsJtg2vdfftBgZmNJhy4LSqzHGlQuqcsg7UE2N/MxqfSXkP+mdJwvNJSwyqa2TOAI4Ffx0uOVxF6GX8i7tCSfHsRetF+OXtpsgxJmaelH0+J93RPAG5x98FcLn054f729iE44z30e+K/PYROUo8DHzazttQ8J9F3WTfvikNZn3P3xwgHGG+O05J8+wAfIdwbTuoCcZ/g7o/HeZxkZnumPjca+BFwCTDK3ZcTOkedYGZdqXzPBQ4ZoH12BW4CPplqn15CJ690nXao2wCWxPc5qbq0E65SpP0e6DKz12XSTwWOIhz45c03r/xRhHvF5fg1Ybnfl0l/L+FWzBFlliMNSmfKMli/IJy9/NHMLiScgZ1CeJyqkjYBfzGzMwk9ik8lXNr7LwB3X25mnyJc6rwp1mU04ZGmccDHBztDd1+RKvOvZvYzoJ1wqbGF8KjMYPyGcH/6h2b2TEJP5qcSHt+61t3vAzCzDxEC6R1m9oO47O8m3L9+q7vvFJTdfesgPncqocPRrTFfL+E+8Sr6Onoll4k/YWZ/cPdfx+W9ltAp6RzC40dvJtwn/mSq49jHgL8Afzezs+lbX8v7axx3fzy28fvjFYS/EXrAf5DQH+HiTN1OsHCj9yf9FHsFcCbhqsTehKslpxAe0Ur7PnAy4ez8bMJ6OorQgfBkd+8xsxWxrV5jZo8QHo26nhCYz4gHpksJtyLKHTP8B4Rn3L9jZocQHjE7kPBI1h2EgylpZkV3/9arPl6U/6hTC+HZ3YWEIHAz4ZGV37LzI1HXZ8o6ifzHSbLzcMJzq18k7JDXEc4wunPq9ybCPdUNhKBxJfD8TJ6dHk0ZoC2OJewsN8YyLyX1eFbMcz4DPBIV882IeRcSAsQjwNdJPT4U872UEADXEi7R/xU4eqB1VM7nYr7nA9fEtlwWl+kpqekd9F0puD+Vfgjh4GIV4bGwO0g9U57K96xU+YsIz/heyMCPRI0nPBf/QGzv5YSOcPtl8p1JeOZ3LaGz1FxyHn2KeQ+P28Qmwu2GM4DjsvkJZ6znEK44bABuB96UKes/4zawgb7HnWYTnjXeQOjI9S3gReQ/EnV4Tv3aCZfnF9D3/PiZpB5Z06t5X+Y+2Ct8IiIiUg26pywiIlInFJRFRETqhIKyiIhInVBQFhERqRMKyiIiInVCQVlERKROKCiLiIjUCQVlERGROqGgLCIiUif+PxAMvkVwRqFlAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 360x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig = plt.figure(figsize=((5, 4)))\n",
    "ax = fig.add_subplot()\n",
    "ax.set_title(r\"Information loss vs number of selected structures\")\n",
    "ax.set_xlabel(\"number of selected structure\")\n",
    "ax.set_ylabel(\"WCSS\")\n",
    "\n",
    "ax.plot(test_num_outs, wcss_list)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6dd13042",
   "metadata": {},
   "source": [
    "### Note that the information loss is almost zero starting at x=3. This is because we have selected the 3 most representative structures: NaCl, KCl and strained KCl"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python (tf)",
   "language": "python",
   "name": "myenv"
  },
  "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.9.10"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}