{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Analyzing a lammps trajectory\n",
    "\n",
    "In this example, a lammps trajectory in dump-text format will be read in, and Steinhardt's parameters will be calculated."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "import pyscal.core as pc\n",
    "import os\n",
    "import pyscal.traj_process as ptp\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "First, we will use the `split_trajectory` method from `pyscal.traj_process` module to help split the trajectory into individual snapshots."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 41,
   "metadata": {},
   "outputs": [],
   "source": [
    "trajfile = \"traj.light\"\n",
    "files = ptp.split_trajectory(trajfile)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "`files` contain the individual time slices from the trajectory."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "10"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "len(files)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "'traj.light.snap.0.dat'"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "files[0]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we can make a small function which reads a single configuration and calculates $q_6$ values."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 39,
   "metadata": {},
   "outputs": [],
   "source": [
    "def calculate_q6(file, format=\"lammps-dump\"):\n",
    "    sys = pc.System()\n",
    "    sys.read_inputfile(file, format=format)\n",
    "    sys.find_neighbors(method=\"cutoff\", cutoff=0)\n",
    "    sys.calculate_q(6)\n",
    "    q6 = sys.get_qvals(6)\n",
    "    return q6"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "There are a couple of things of interest in the above function. The `find_neighbors` method finds the neighbors of the individual atoms. Here, an adaptive method is used, but, one can also use a fixed cutoff or Voronoi tessellation. Also only the unaveraged $q_6$ values are calculated above. The averaged ones can be calculate using the `averaged=True` keyword in both `calculate_q` and `get_qvals` method. Now we can simply call the function for each file.."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 42,
   "metadata": {},
   "outputs": [],
   "source": [
    "q6s = [calculate_q6(file) for file in files]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can now visualise the calculated values"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x7fbd87929b20>]"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXQAAAD4CAYAAAD8Zh1EAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAgAElEQVR4nO2dd5gUVdbG3zMzDDkIDAJDGKLkDJKDoICgGBd0zbKYV3f9XMGIGXVXXSNiWDOsAYUVEBERVOKARMkwwADCkDOT7vdHd8/UdFe4VV3VFeb8noeH6erqqnOrbr11wznnkhACDMMwjP9JctsAhmEYxh5Y0BmGYQICCzrDMExAYEFnGIYJCCzoDMMwASHFrRPXrFlTZGRkuHV6hmEYX7J8+fIDQog0te9cE/SMjAxkZma6dXqGYRhfQkQ7tL7jIReGYZiAwILOMAwTEFjQGYZhAgILOsMwTEBgQWcYhgkILOgMwzABgQWdYRgmILCgM0wpQAiB5TsOww/pshdsysHOg6fcNsOXsKAzjI3M27AfM9fsdduMGL5dvRdXvrUQX63Y7bYphtzw/lL0fXGe22b4EtciRRkmiNz8wTIAQNaEYS5bUpIdB08CALYfOOGyJYyTcAudYRgmILCgMwzDBAQW9ICxcMsBTF+1x20zGIZxAR5DDxjXvrsEAHBp+7ouW8IwTKLhFjrDlCJ84LXIxAELeing4ImzGDd1Dc7mF7htCsN4gjN5BXj3520oKAzWG44FvRTw3KwNmLx0J/63ynv+0UxiIXLbAm/w6tzNeHrGekxdke22KbYiJehENISINhLRFiIaq7FPfyJaSUTriGi+vWYy8VDI/WxfcfJsPk6czXfk2FwVQhw/E7q+p/OC1Ws1FHQiSgbwBoChAFoBuIaIWkXtUw3AmwAuFUK0BnC1A7YyjO8QQmDSgq3Yf+yM9G/ajJ+NNo/PdtAqJqjItNC7AdgihNgmhMgFMAXAiKh9rgUwVQixEwCEEPvtNZNh/Mnm/Sfw7MwNuOuzFdK/Ubaiz3tkFoa8ssA2e7SGXPILCnHsTJ5t52HcQUbQ0wHsUnzODm9T0hzAOUT0ExEtJ6Ib1A5ERGOIKJOIMnNycqxZzDA+Ir8gpM6RLr5ZzuYXYsMfx22zR2vI5R9frka78d8XfT6dW+CJnDR/HD3jaKKuoA1ByQi62js9+jKkAOgMYBiAwQAeJaLmMT8SYpIQoosQoktaWpppY5n48EOmPTWuf28JLvjXT26b4WvIYDZ06m8lk3aNn74Od366Ar/tPOykWYZ0f25u3Im6Dp3MxfIdJcth1+Rw12d+wINfrkZBocCyrEP2HDQOZAQ9G0B9xed6AKJDEbMBfCeEOCmEOABgAYD29pjIlHZ+3nwA23JOum2Gr5F9mR89nYfdR05j95HTAODY5GwiuWriQlz51sIS2+xq2+QcP4v/Zu7CWz9twdUTF2HR1oP2HNgiMoK+DEAzImpERKkARgGYHrXPNAB9iCiFiCoAOB/AentNDTZbc04gY+wMR9/yRq00J9h58BQOncxN+HkZawz813z0mvCj22bYil5jwK5HYvP+UBbLfSYmv53AUNCFEPkA7gYwGyGR/lwIsY6Ibiei28P7rAfwHYDVAJYCeFcIsdY5s4PHr1sOAACmr7Q/D8tUF3Ng931xHvq+UHpzW3vF71v2ZX7gxFmHLWGcRCqXixBiJoCZUdsmRn1+EcCL9pnGJJLHpq1FtfJl8PeLzrP92EHotlvFK9MWVudPvGK/LCfP5mPfsTNonFZJan+/lc8IjhSNg0VbD2JPeKzRD+g91B8t2oFXf9ySQGsYxn5u/s8yXPAv47hGu3tOXkkhwIIeB9e8sxiDXvJ/UGzO8eJudvdn5/rqJaXFbzsP4/Fpa33r2WM3VudPvDJkJMtSlzxNvl0dcvF0+3qxoFtg16FTuOk/SwEAp3L9Hzrc9Zkfiv7+49gZTItzHL/nc3Pxwncb4jUrLq6euAgfLtqBfJdaTvkFhQDcf8DjRfk+PJWbj5OlePhMBrfbDyzoFnhu1nr8tDHxgVGPfrMWvZ933gNBxIQZmGPP0TN486etNlnjT/72+Sq3TYgLtRdR68dno7VOSoLDJ3Nx35TfPDNnsvGP41i7+6jqd24Lr1OwoPuIjxfvQPbh01i564il37vhtug20Q9uXkEh/jjqvGvZ/xxYNSo3vxD3f77KtSExIxF8Y94WfLNyDyYv2Sl9zPd+2Y63HHr5D35lAYa/9osjx/YqLOgWINXg2cRx2Ru/WvrdhFkb8EXmLuMdA8y4qWvQ/bm5OJWb2FakHS3Cnzbux1crsvHYtHW6+y3ceiDh/tCdnpqDd3/Zbvp3T337O553YHjO6P4GtW3ja0Ffu/so1u89ZsuxhBD4ft0fKLRhzPVMXkHRGGqEwydzcSavAL9uOYDRH2aWmKwb+9Vqw4c0GiuhxgdOnMUDX6429Ru/E/3g/rB+HwDgbF6hyt7usXDrAfy4YZ8tx7r2nSW2tUxveH8pVkn0CJ0IHtsSDtaxwvnPzrXREnncflH4WtCHv/YLhv77Z0u/3bL/OPq+MK+oIn79226M+Xg5PlqUBSBUQTfts5YUqcWj3+HG8KRphI5PzcHVExdh9IeZ+GH9vhJ5mKcsM99qfnXuZlw9cZGlyNIdB43D6Af+6yc8N8tcsG/O8bP4arm3FgzQahknegjV6EG/9p0luOWDTNvOp/RcipfF20qGs8cjtGZQepC1fPQ7U66BVpOh+R1fC3o8TJy/DTsPnSpqse07FnoA9oa7qkNeWYCLXraetvTXLbE5HdbsPlok5NpCI3D/56uwZJt+TojN+0Mvm/3HzD+4/V78yXCfrTkn8fb8baaO+5ePMnH/F9qTgUIIvPfLdlWxeXv+VtMvED20BNTNBtRtH2fi29X2j60nGjfKcDqvwBdLKLo92RpYQd9z5LTUOOL8jTnIGDsDm8Ot8UMncnHkVC7267VwHFSF3PxCfLUiG9e/t1R3P7crjhp6izg88MUqNBo3E099+zvumRybG/y5WRtMv0DU2LL/BFbuOmL6+mw/cFLTI8IOhABmr9uHuz/7zbFzWCXSK/U6XqzzMgz45094/cfNCTmXVOi/H+kZTjCUNWGY6veRyvHduj8AACvCaUK/WJ6NL+IYNohO06kFUUh84l0ezokxOyceHOU1PXbaue5wpJuekmTuwgz4508AtOuLVZwZU7X3BhnN3/hUR6WYsWYvbuyZYdvx1O739gMn8c/vN+GKTvVQt1p5286lhi9b6F8uz0anp+ZI7Tt+un5ljVx/uyptdJpOPQa9ND9mWKdITA2EwA7RPXk2H+2f+N54R0lkTbJq+l8n/4bL35Tz8Ik8WHkF3pr89CNChOacIsh6edmd6MuJF8vS7YmLLO2ZgCyWvhT0h79eIz2r/sHCLNXt8QbPyJBXUIgzDi9Ce+enK3Q9c7IOaE+Abtl/AkdPxy47NjvcazFi7e6j2K44vtNd4umr9uC3neZ88Fs/Phtbc4on8Q6fCpV3u851cROlrVp1Z0d4BZ812Ucx4vVfcFojWnnG6r0xqQ9kb9HPmw8U/f3hwiwMeqm44SHz7Pxv1R50efoHZDoQii+EwOSlO9H28dmmvNL6vjDPsIHnd3wp6HYSack50TO+5LVf0OLR71S/s1P8zubHtkL3Hw+NZ/cPDyWYYXW2+lhyQaHALR8sK2rVDH/tl6KhCjNErvXPm0PzF7sOObfEGABs2BvrrTR3vT0ugnYzUJFYSivt8Ob9J3DsTB6e/HYdVmUfxRqNsf+7PluBeRvjX973D425kWkrdyNj7AxkH469f5E68ruKW3H24VNxTXDe+ekKjJu6BsfP5iOvUL4HtvPQKc0GXjSHT+bi9z0lbd979DSenbneFtdmpwicoE9euhMZY2fYdryfN1sP8bdzLUg91FpMP643fpDNVsuDJ87ixw37NRc8lu31RPb6MjyunrnD2W6vgMDnmbswbaVxXvh47rfd6E3Mtxv/PTZK1K8DJ3Kx69ApbA27GtrRcIkMuXwTXrbOjHtvbn4hej8/D3/770rT5430Nmat/UOxLeRWafeE9qVv/IKLXy3pEv23/67EpAXb8Nsud5fl08OXgq430TRu6hqLx1Q/6F2fxoqX1ullBCPCR4t2qG6XHEIvIZ5qrf1EtiEiD5rZXoedvaLN+47rpkT4x5erce+UYhHRMtXIu8hLHJPwtf5oURb6vDCvaM1QO+pFbkGB7rH0ns/8cIt63gbzL86b/7NMdfuoSYtNB1IZ1b1dh2LTK+QVWKvnicSXgm4LUTfFjjSr/zURIGRmbLGwUODhr629qOJlx8GT2H3kNK6ZtBjHzsSOtwPANyZeZErsfC4ufHmBpZQIB0+cdSTvildYu9ueSGolb8zbinkb9hcJW/QkqRDavbV4HrNMSQ8yGdTMmDBrg2r8x+SlO3XnotRYk30UGWNnJDzvTmDdFmUJVUZ7pMUOF7Wih0RxrC05J/BpVMIjowdD5sF5ac4mw32UQUgzVoe6ujnHz2LehuIhnQPHQxPUHm64xBC5PmM+Xi7tahrX+Wy4Olq9SDMRlHb1iuZvUrSwJQ+6//gZdHvGnZB8GSbO34qJ87di3RODi7YJITBu6hpUq1AGTSRXQQKAT5eEeuAlrlMCKLUt9KJHwMZ+v5mkXfE+3psV4ddqx5IRkAUmK5vymO9JJGJ6Z4F6oFCkN5SoqE21YJ5IWewKYz+TV4DvVbyD7EzkptWLnPDdBuQXFOLQyVys0ZjQLjqGbdYUc/RUbM9NWe4Fm3KQmXUIOw86O/ltF/dOia0vRxRllLmGbg3LlPoWuhWcTEMbeWjzC/RrREl3wdh9na5Qai+M6HM+M9O+UH4ly3ccQueG1Ys+yySP0kLNbVONDxdmYdbavZgypofq9099+3tML0qL/IJC7DlyBgICHyzMwsiu9aXtVWPVriNo+vAsNKheATsd9hiKoHwE7jOY4Lzh/dC8xBe3q1+7CEP//XPCPEj0nuB1Cu8WL4+Xq+FLQY+n1bM6+wja1K0ae8wEp0nTOltk8iqy0s78TTlYt8daq2uzjveBlYAbrZdM8aWTq/3R3j9aD83Pm3OQkhTbibzyrUUAgCs6peOlP3XACIvphLXIzS9EakrJ8z5u4L+sJaRqL74XZ2/E24rey39+zTJvpAkbnEJvUtTKEJNdmVNlzjN1RfG8z5m8ghKflaiVYt3uo2hfr1pMHfECvhR0sxQUCnyRuQuNalbEyEmL8cDg4pXtc8M+3AfjjGo7nVtgKsf23A1y/sE3vm/sdRGK5DuBDX8UPxBvL9iKR75Zq/mbV34wHj+P5vV59i0ivXb3UcPWsZHHydQVu/HSnzpYOn+uiu9+hH9+vxEPXdzS0nG12LSveGjnly0HdPb0D35drzU6Q+srP2zGxPn6i2wo23vj//c7tuacxFOXtTE8V6KTwZUKQf9syQ48Om0d2tcLtcxfnL0Rl3dML7GPlguYjGsYAHR75gcc11h667edzk+6RS9WreZ2pST7sP2z7wdOyOfEtitf98i3F6luzzMYsvrPr1l4/JLWqt+Z9WgAtHt4ar1Jq1GqXlpxavbaP9CklvokoRDq5VYLQLITqy8YZVoDrWMKUTJP06ps68N8TuK9PoMDRCY0DisnNkzc/IJCgQ8XZhW16tQeKy0xB4DL35TP76JEbbJJFX82lIpYuFU/VbAeSxKYi8MurCwsvveoO8vOabHHwjJ+f/uvN9dZPW0hPYeZOJFEUioEPSmcec9qZsPJS3fi8enr8M7P+l4bdnPNO4ul9rNSeaatjM/3Wi3fu1W+9NiiGF7k5Nl8Xw1xROqk3QtGR8cafLJYPUAvHpRirXnFPdRbUhIYQT+bX6DZdYpce+XzsMvEkEMkoEYrsKbRuJnSxzKDWh4MNXz0nBsyean8AsNOsmjrQbyr8QIHQpkEM8bOwJzfE5MT5raPlyfkPHag1LoXvtso9RvZl1V0NPDTM+L3pMrLN5/A7JTNLyq78KWgq70cH/tmXYmMcEqSKLaFbiaYJPKz2Wv/wJJtBz33cg6QnmPc1DXYvO943JPU8XLNO4t1xSLikfGhQbKnr1bY0/vYmuO97JBBaUgs1Yna1irjZsn4hURrRSAmRY+eysNcnQV2X50bWi0k3gqYdfAURk5ajFqVy8Z3oACi15o1y4UvL0C1CmVsO56TRIYW1J7bvUdPY5JGcFUQ0PPWMetarJXh0wx2vWCUljd/ZJbmfic92EoPhKC3f1J/kQYrk1B66C5P5wJuj61uP3BSOqhGliOyE8IOoDW0pkRGsIyCw8xit5fL0dN5OC5RVrPo5XLRwu5YgniQtbzL0z8Y7pNrcx0wIhCCLovVSVG3BdMIt9MzR68K7xe07uvibcVd8F2HTqF+9QqGx3JCGJ3mwpfme65x4idkvGMe1YkFcQJfjqFbxaruGfk0u41bLlIRvDi+K4PMhJqWF8XS7aGXWGSowG5vDqdQNmqcEnOvzTGZJR7zP1m8A0dP5WnGGrw8ZxN2HHTuefGloFu94FZ/l29iVRQ38Lh5nkUmwZjWq/LVH0NRs8fDgWeJ6MTZ0VMsKBDYsv84jpySDwKzgp1JyWTZpBMglCjW7TmG4a//jGVZ6k4X/567GTd/oJ7X3Q5K1ZBLUCnw+JCQn5ER0b1HT2NHAvKojLHBdXHN7qN496XtqF2lnA0WaWPUazST8leWIa/8bLyTwwhhHKWtl3YiXljQJfC6Xn7wq3FLk7FG9L0/fiavxMpHANDjOfXV3L049BCJrNVaJ9QujhukzMi1kBwuUViJgo3g9vCnL4dcEo3H9RzTA7zijtf4ZuUe/CiZWK20IkT8kchBxok8ShGkBJ2IhhDRRiLaQkRjVb7vT0RHiWhl+N9j9psaP1YngbzeQncbL7ZE7YJvPeMEaouh2IHhkAsRJQN4A8CFALIBLCOi6UKI36N2/VkIMdwBG0tQWChw0ma/ciPc7kYZ4eSYnAxBfuFFly3Rbmh+RDbZldaKVqWBHQ6t3iTTQu8GYIsQYpsQIhfAFAAjHLFGgk+W2J+MxxCPC9ZhF4Nwgo6AsDyB56V0t4lEZq1awLkVrdzE7caNjKCnA1AuZ58d3hZNDyJaRUSziEg10TQRjSGiTCLKzMmxtniqk+NPDBONEMDHi7LcNoNhpJARdLVmRvR7aAWAhkKI9gBeA/CN2oGEEJOEEF2EEF3S0tLMWapjjNNYjTBl/E9eQSEOcQ+IkURWKZzqvMkIejYA5Sq29QCUmMIWQhwTQpwI/z0TQBkiqmmblS7jdmg94x5W1l6NcP/n+osnM8FDdj7LqTaijKAvA9CMiBoRUSqAUQCmK3cgotoUHjAkom7h4/ozwYcKMhGFTDCJ52WuzAnDMInA0MtFCJFPRHcDmA0gGcD7Qoh1RHR7+PuJAK4CcAcR5QM4DWCU8HpGK4aRgGsx4yekIkXDwygzo7ZNVPz9OoDX7TVNg9LpOMC4hIDAl5m7jHdkGBM4NYbOof8Mo8PUFbvdNoFhpOHQf4ZhmIDAgs4wDBMQWNAZhmECgu8E3Y3E+QzDMH7Ad4LOMAzDqMOCzjAMExB8J+ilNIEdwzCMIb4TdIZhGEYd3wn6qbP6axUyDMOUVnwn6Hmc+pBhGJ+z+4gz6zr4TtA5WRLDMH7nk8XOrLzmO0H3/HpwDMMwBji1PKHvBL3Q3fWQGYZh4sYpZz3/CTqPuTAMw6jiO0EvYEFnGMbnuLmmqKdgPWcYxu84lZPKd4LOMAzjd7iFHoZD/xmGYdTxnaCnJLGiMwzjb9jLJUz3xjXcNoFhGCYu2A89TKWyvK41wzCMGr4TdIZhGL/DQy4MwzBBgb1cQrAbOsMwfodb6AzDMAGBJ0UZhmEYXXwn6Bz6zzCM3+FIUYZhGEYXFnSGYZgEw5OiDMMwjC4+FHQeRGcYxt+wlwvDMExA4CEXhmEYRhcpQSeiIUS0kYi2ENFYnf26ElEBEV1ln4klYbdFhmH8jlNLaRoKOhElA3gDwFAArQBcQ0StNPZ7HsBsu41kGIYJEoWFLgk6gG4AtgghtgkhcgFMATBCZb97AHwFYL+N9jEMwwQOp0YaZAQ9HcAuxefs8LYiiCgdwOUAJtpnmjo84sIwpZc26VXcNsEWjp/Nd+S4MoKuNiEbrauvAHhQCFGgeyCiMUSUSUSZOTk5sjYyDMMAAFKS2I9DD5mrkw2gvuJzPQB7ovbpAmAKEWUBuArAm0R0WfSBhBCThBBdhBBd0tLSLJrMMExppWcTXoJSDxlBXwagGRE1IqJUAKMATFfuIIRoJITIEEJkAPgSwJ1CiG9stxbs5cIwpZm7BjR12wRbqFkp1ZHjGi7QKYTIJ6K7EfJeSQbwvhBiHRHdHv7e8XFzhmEYAEhyKk1hgnHIycVY0AFACDETwMyobapCLoS4KX6zGIZhYgmInqPARbdFT1G9ojNdFYZhmETh1FyA7wS9B0+KMAzjc6pVKOPIcX0n6AzDlF6CMuTiFCzoDMMwCYfT5zKMJ2lfv5rbJjAMAJ8Kes1KZd02gWGK+O+Y7m6bUGoow5Giuvjy6giOLmI8RLkyyW6bwDAAfCroDOMVrj2/gdsmlCp4UlQfXwq6Xvt87RODDX/PlYKxg2oVyuDZy9u6bQbjQ5zSIH8Kus6QS6WyUsGvqky8rrPl3zrBjT0axn2MB4e0sMESxoj61cu7bUKpwKnFlYOCLwV9cOvacf2+blX1h++CFrXiOq7d2BEd7Eb9L42ThEPb1HHbBMZGhrV19n7yItEK+jWPL/VuDY1MZ8Jjy2fYYY8f2zN9mtXU/T69mvdaw9wTChZ+nRvxpaDH2+pMTvKHzNmVv6dyHMNQ5cskIyXB18sv3WqllX6pU0yw8amgx/fwJFqgrGKXd+aLV7e3/Nv3buxijxEmMHJLrVO1XIIs0cdb/TkGALpmnGP5t63rFi9vV9sjdcws/hR0C7/596gORX+r5VT2ikgoMUqC3y2jutRx4pko7tlUf/hDDTuFTu3dWz6V/b5LI50aGEfknlPBejbWKYq5nyZplTDv//rr7t+gegXL52IvFwVWktwrG31+EYRbezfS38EfHY24qF6Ro4L9xrlV3Ltn8Qhl5XIlMyA2qlkxTmsSjz8FPQ6rK6QmfkzYKnaNy8pU8vsGNYvr904x9Y6e7p08QUy4Ili+7E4FcmsNtSrniEp7ELkvBZ00mqYXtToXAPD4Ja1ivjPyGEmEZpntokUqcAWNHoWMzURylfy+Qc1NWCZPtFg9OrwVzju3suq+t/drErOtQY3Ya+aXSVNZmp1bybH82E4wXuX5UuKUpmrNrayRCCZ0ArVqWKWc9eFNO/CloGspWeQCp6bEFssLb+7mGkJmlZFd66tud7sH0jgt1FX919Xt0a1RyXH+IW1qo1wZ9WqntX39k0Nw36Bm6CCZ1fDJEa3xysgOxjvGiVaduqlnhuPndpNB4YaTGl50KU0Uc/7WF7+MvUBqX61Gabz4UtC1xtC9INp6dG8sN4kZjVppJ17XCVd0qqfaG4n5fYL1vVblcsiaMAxXdq5ny/HKpybjvkHNdYeFohnRoS6+vL2HLefXQuu6Vna5labGi1e1w7KHB+HlkSGPp3b1qjpyHiGEq89hIut69KmqViiDSqncQjdNvXNKtgKu6Jhu+Buvi71Vbu6lP3HqVEsgmrTKchNhQghc2kH/fmndKzO3kIjQRcUL6OrO9fDr2Avw2HDjF2EicfouCYTuUUaNUO8pnmGrii6LlhVqSdbPeElyuXfsS0GPdjG8SCIVQEQMtC63VgX/5cEBJizzBkrhK5+ajM4N1X1zM1TGpwFg8l/Mh+6beWHe0itDdbvsy8dor0va1dX87sWr2yO9Wnm0Uvgc6zH3/n6a39nXSCDfzAukpiThHMOF2p1pPaldo0Et5dJ1DGypPUykx4e3dEN7jd5MtD2Jajzp4UtBj71w8VWgfs3T8Mqo2DHXT249H/XOse5rGs1lEj0JJfoBNvKVRytft5aPuZ0LcauVQEu8DFMdaHwdPfGqJTjKBFrnN6oeM76vRpO0Sob7xI8/uo/PXN4GM//a23A/p3rDas9DRckYC6t5mvo1T8O0u43LDJgb7mE/dBP0aRqb66VCajIGtqiFt6+PjXz88JZu6KrSPbcrg97Gp4cga8IwyystqQlg3+bmA35ijiu9X+yeqx6/KO7zWyX6cjSpZV50iUhqqM4KKV5cVUdCZH9/Ut9b5M/nN0TTWvoT+2mVy2qeauxQa/lu9CZaZeuwHc+LEe63zwMj6CUvZYMaFZA1YViJbUkEvHdTV/Q2SPxkN03SKqJsir2BTE+OaI0KNoxjyrYSolvOE6/rjKrl7Xezi7w4tFrqXkuepsVf+jbCaKOgsBJ4QQpgS53SC9pTuqX2airfC4zUUyM3Xt1jOHCNx1/auuQ5wvbJxI9wtsU4GVRiDK34cvrNxWzWvX1wfff486QD1iv5kDbxpS/WwkiwI4JjdYIrurxOdXsrpKbgEZOTrs5Pirr7Mlz+yKASn6uVlw/Rj75PM/7aB701hgu16rQT97pf87QSL6nIKfS8nCLPbiWHPKECLeibnxla9HdKsrKoxZVb6ducmuyty6H2CLasU8XUBJodFdkLkz1AaNz7hava4fFLWhvvrILMtWiTHjtZqpUzR3Zt25Fd1OMFrHJ5x3TUqJhawm//jv5N8PGt3Uwd57v7+mjmK/l09PmWbOvdtKbqC7dG1HBjPDEZjWpWxJWd1YfL1F5cdrivzr2/H/5zc1fd86nVr8ZR6QMu75SOR4a1xD0XyLvgmsFbCmaRSCBNdEBRGQOBvq1fY4zoUOwRQUQlVglyYlhBFmVFMCunSqGJ1pxHhrWMw6r40dNAoxcHEeFPXepLT4TFHt/4fN/e0ydWkCRvwCANT4p/DDnP8Lct6oQE7rZ+jXU9awDg5ZEdsPzRC/HLg8VBLA8OaaGZe2RI69oYruL506J2Fc3f9LQ4MX5n/6b4v4tC5e3dtKbm8Mht/Rqrbo8WQCV69VqP82pXjrtJ0iStEgacpz+xGqlPynNFz80lE2F0n8aOLSzuP4dSxL4JB7SohTv6N8Ff+qhXEhNMUe8AABLPSURBVJUjAAA6Nzgn1vVI8blaHJnbZGhcsyK2HTip+t2EK9tZPu4FLWrhh/X7Vb8b3acxnp6xHoDd3VCNcW8LPX2zv7HsVWFQ/s9v02/ZRdedt6/vjPzCQvNmEFAuPM/StWF1ac8a2Un2idfLL6048bpOqKOxopcMSUlUomG1cOwFOJVbUPT517EX4NCJXM2698jwlrjlg0zp8+n1ViukJuNUbgGIEusWqlcdnR74CkQLPTmJ8OCQFqhu6B9rjVev6ejIcYHQeHCPxjWKIvjsoHvjGrimW2jFFTP1uKGGX3ppxeh6RA+5JCeR6gR4ZK/K5VLQNl3dpzkyWe/2PRjSpg7aR6VYaFSzoim7imI+KNQoqqvwUkmvVh5t61XVHuvWecuqiXL0PVD+/uWRHdCqThVUkGgNa/WuzEAqaho9BOT0ayUQgu4EtygiMCNuU23Tq1p2vdJi6cODMHlMd1zeMTZM3uno1ujK9cXtPUyPw2ofTY6FitwX0ZXfKB88AHz2l/OlX1qyrTS7L3vkPqYmJ+E1jcbBTT0zsPShgWimMbaslp8oGidaoa3rVsG8/+uP+Q8kLsBu2l29bDnO4Na1MfPePkhKcm4WSHlkivpfdX+HFb1UCnpa5ZBQqLlpnR8ONlHz5EhJJtWMgEpGaSTMiguTlaBO1fK4tXcGalUuiwt1EilFU6tyOfRpFt96rWapW618TCqHCDJ573s2qan74vvqjp74IjwpZvlZ0jh+Uwv+71oQEWpV0V5kxUweGz0iwzTdDYKqiAhZE4Zhxl/7qH4vs9iE/vG1v2tRpzLqVi2HS9uHxv3NtNpraeRiN6pKVoX2zgH6enCvQ1lMtSiVgv7o8FaYcEVbVV/YoW3rIPORQapRhDIt5odcnnS8a0ATXNy2NprWqoylDw9CrcraIqHX6ptwRVsMkUipYIx8e7fID93GJnLnhudormIj+wxrufzdM9C8yLod4V+/egXM+7/+eGCw8UStHs9dYX2OByh57SdeV3KMv2xKMhaOG4h/Xt0e13VvgIfDz5SMV9FDF1t7/qzWuSqKRTHUekgyvUw7KZWCXiE1BaO6NdDsokZPNsk+hD0am/MMuKO//tvdSt+/T7M0w673D3/vCwC4opN2NsRR3RoUT6ZJld+6CjcOTwKWjUqfa5/4hW2zWUxlV86KvBDcEPPoADsgNCaeEqeLrl7sjKw7Z4QSaZMVx01NScLTl7VFDZW5sd5NayKJgBuj4ki0vEcSMSkafYbujasn3BW6VAq6U0we0x2VUlPQt3maYVrbrAnDcLWkf7LdVbFprcrImjAMLevIJahyemr+9Ws74oObuxa5C0Z6Dg2ryy0B1qF+NaneRIzbouRDbuT+Kg+ptgTtur9eWYireFJUbqikX/PiYb7qOp5lyt/UqlIO254bJp0j3wg79D7aRfPNP3cGEeHpy9rEf3BJpGoqEQ0hoo1EtIWIxqp8P4KIVhPRSiLKJCK5bDY+Q0bXkpIIH93SDX00Ugx8Ovp8vPQnOY+WyKx5FRf94RNBlXJl0F/h41u9YireuaELJt0g526XmpKk65rXqGYlXN25Ht78c8l9ZJ/hug4v2mCXT3JtnTF4I8zMtdiNUqijPWyUmG35JxqtF9h13RuqBqw5gaEfOhElA3gDwIUAsgEsI6LpQojfFbvNBTBdCCGIqB2AzwHY6w6iQfNzE5ENzz56aYQsq1GlXBk8cWlrXNCiFvq8MM9BqwyIo/USLVayC3TbKTDJSYQXr7bPLdQ0OjqUnETyPSUVnr+yLT5YuAOA9WGF1eMvQnmTL5V4pdUjnQlHceMFJNNC7wZgixBimxAiF8AUACOUOwghTohi6ysigflAp92VuM6AmUpo5V6uHn8Rnggn/In4/d7YMwP1Ta5F6g7qV6feORXw9vWdsfShgZh+dy/LGSedwNDzQWN7ZE1U2QdW6ZcdzX0mJlbV1qQd2bUBZt1b0hPF7DJwVcqVsXFYCVJTFlbePU6Mg2u5kWox0GIa3kQhEymaDmCX4nM2gJhED0R0OYDnANQCEDsTE9pnDIAxANCgQQOztqoi2+KzA63Hd8B5aZi3MSfu41cpVwY39GiI67o3lMrY5i20xW1weHxbzy3PLlrF0drVIlq3z61aDhv3HTf83bS7emFZ1qGiz8o72rBGBTx7eVt0NzGRPvu+vsjN149EnX1fX9SqXBYdn5ojfVwryNROtzx6XhnZASt2HpbaV83Ty641ed1YtERG0NWsinl6hRBfA/iaiPoCeArAIJV9JgGYBABdunTx9oCYAqPbMumGLjEPmtV7SURI9puWe4TNzwyV9jxJBO3rV0P7+tWw79gZ1e/NDL8JEWq8GDVgzqsd6j08OaK1o7mI4h5ycfA+XdYx3fRiMhE+uqWbZn4bs7gx5CIj6NkAlO4Y9QDs0dpZCLGAiJoQUU0hxIF4DVTDQ88sgJAXhK1dVp8wqKV2zhg3sPseqD2O79/UpWjMWvZxdWMu74YeGQk93829MtCxgfpSh3bhhEBGH7Jv8/gD63S9exyePZB5ApYBaEZEjYgoFcAoANOVOxBRUwqXgog6AUgFcNBuY90icoN8NwriMJNUVn+Kh4bhBYydWpE+GvmUAcV/X9DiXMuPJJH1lq3XGjHRtE2vWhTZ6Ubu9ceGt8L/JJeKi4fbDCLFAXe9cQxb6EKIfCK6G8BsAMkA3hdCrCOi28PfTwRwJYAbiCgPwGkAI4XXfYxM0C69Km7p1Qg3ayxurEZwSq9NPCuct69XNSabZeeG5+CHv/dDkzR7urxGxN1aSuA99lp90pvjiWSO1IrQtYpey/cWUytEWadbo+pIq1zWcC7DLaTS5wohZgKYGbVtouLv5wE8b69p3iEpifCYQaBQkHGicai18K6d+VG8QolFEFy0w0708pb3aFID4y9phSs7a0ciew0z92XJuIH6x1LNCmnSIIuUvoHfBOFWF9ls+tXWdRMT8OBFrLotWj4fyCerohoTLVrKFjsR4aZejVC5nH8C4szcl6Qkiqt36iQs6Bb4VZHy1UtsfHqI6cUJvrqjJ9aMv0h3n09Gn48rOlnzGggC0a0rsy9rJ8P9vcCYvo1xcds6bpthivsvbK6Zm97P+FLQ3V7j0mzghtNM/kt3LHhggOriCkaUK5Ns2JLqmlEdL/2pg1XzfMeA80KeDkYugrKTf5HUDVZd6bzOQxe39J2X1z0Dm+F/9xQP+9mpKK9d0xE9GtdwfMUzNXy5BJ2fSERwQQ+J9R+n3dUL8zZ6x8XQCwzWSOj19GVt8deBzTRfdGYDTyqVTcG6JwajfJlk7Dp8CoBx/paKqck4qVi6LegseWgg8guDMSDVq2lNUzEGduKv16qPaFyzEm7p1QiTTKzn6CTt61fDfQ4m24/0WqwuLmyWd27ogovbxpevvVyZZCweNzAmjXFqShLqnVM8FxH9Tn728ra4sUdD9DWxGEjFsilISiI0qF4BDww+D+/eqO/y+dWdPaWPHQTOrVLOcz1fJ3C6fcctdIcobZ4xkRwtN/dKjPvYha3OtSWBV+2q5fDgkBZ466et0r+pVaUcnhhhLSUqEeGuAU0t/ZZhjGBBZxhJRifI15lJHANapOGKTulxr+BkRKLcFlnQfcLP/xiAvAJvBjMA3o9kNOKz0efDaAi3eW31BZxLM40TFATmFGVTkgM14e9LQfeSeDiyKLQKXk+h67VIRrP01JnEcqNskTH89GrlsfvI6cQbIMG39/TWXOCbcQeeFI2TCVfGt1guw6hRqWwKsiYMw6Ud6rptiiZt0mPTNzhBZP3TimV92f5MKHyFLPLs5W1xOq/0uJWVZtzsEXqoM+oa7etVxdihLXBVnKkE/n5hc7w0Z5NNVnkTFnSLXHu+PQt0MAyjDxHhdoksh0b8dWAzFnSGiTD1zp44V2PVIS/NazCM18ioWQG/7z2GCg6vsMaCzkjTSWcBA79PijKMk7xwVXtc0bEeGqc5m03Ul5Oi3BhkGMZPVCqbgkE2BMIZwS10Ji5SkkhzGIaxjwCtF8M4iC9b6Ix32PDUECz4x4Ciz0HsPd3auxGqlEtB//PiX2+SKYkd6RuYYriFzsRFis/SplqhZZ0qWD1+sKs2JCJrZ6JZ98RglE0Jfv1JJCzoDOMDgjjk4lagUEWHPU3chAWdYTxMABvmrvL6tR3RLr2a22Y4hi8FPYjdT7/z2CWtQAT0buZOYv+gEsCGuasMb+fdVAp24EtBZ7xHwxoV8e6NXd02g2FKNTwjwTAehjujjBlY0BnGw0QWXy4N3kRM/PCQC8N4mNv6NsHp3ALc1DPDbVMYH8CCzjAepnxqMsZd3NJtMxif4Mt+HA8rMgzDxOJLQWcYhmFiYUFnGIYJCCzoDMMwAYEFnWEYJiCwoDMMwwQEXwo6R88xDMPE4ktBZxiGYWLhwCKGYZg4aFarEtrWq+q2GQAkBZ2IhgD4N4BkAO8KISZEff9nAA+GP54AcIcQYpWdhjIMw3iROX/v57YJRRgOuRBRMoA3AAwF0ArANUTUKmq37QD6CSHaAXgKwCS7DWUYhmH0kRlD7wZgixBimxAiF8AUACOUOwghFgohDoc/LgZQz14zGYZhGCNkBD0dwC7F5+zwNi1uBTBL7QsiGkNEmUSUmZOTI29l7HEs/5ZhGCaoyAi6mnqqLoxFRAMQEvQH1b4XQkwSQnQRQnRJS0uTt5JhGIYxRGZSNBtAfcXnegD2RO9ERO0AvAtgqBDioD3mxc/c+/shiVv0DMOUAmQEfRmAZkTUCMBuAKMAXKvcgYgaAJgK4HohxCbbrYyDJmmV3DaBYRgmIRgKuhAin4juBjAbIbfF94UQ64jo9vD3EwE8BqAGgDfD49v5QoguzpnNMAzDRCPlhy6EmAlgZtS2iYq/RwMYba9pDMMwjBk49J9hGCYgsKAzDMMEBBZ0hmGYgMCCzjAMExB8LeipKb42n2EYxlZ8mz73kWEt0bc5R5syDMNE8K2gj+7T2G0TGIZhPAWPWTAMwwQEFnSGYZiAwILOMAwTEFjQGYZhAgILOsMwTEBgQWcYhgkILOgMwzABgQWdYRgmIJAQqsuDOn9iohwAOyz+vCaAAzaa4we4zKUDLnPpIJ4yNxRCqIbJuybo8UBEmaVtRSQuc+mAy1w6cKrMPOTCMAwTEFjQGYZhAoJfBX2S2wa4AJe5dMBlLh04UmZfjqEzDMMwsfi1hc4wDMNEwYLOMAwTEHwn6EQ0hIg2EtEWIhrrtj3xQETvE9F+Ilqr2FadiOYQ0ebw/+covhsXLvdGIhqs2N6ZiNaEv3uViCjRZZGBiOoT0TwiWk9E64jo3vD2IJe5HBEtJaJV4TI/Ed4e2DJHIKJkIvqNiL4Nfw50mYkoK2zrSiLKDG9LbJmFEL75ByAZwFYAjQGkAlgFoJXbdsVRnr4AOgFYq9j2AoCx4b/HAng+/HercHnLAmgUvg7J4e+WAugBgADMAjDU7bJplLcOgE7hvysD2BQuV5DLTAAqhf8uA2AJgO5BLrOi7H8H8BmAb4Net8O2ZgGoGbUtoWX2Wwu9G4AtQohtQohcAFMAjHDZJssIIRYAOBS1eQSAD8N/fwjgMsX2KUKIs0KI7QC2AOhGRHUAVBFCLBKh2vCR4jeeQgixVwixIvz3cQDrAaQj2GUWQogT4Y9lwv8EAlxmACCiegCGAXhXsTnQZdYgoWX2m6CnA9il+Jwd3hYkzhVC7AVCAgigVni7VtnTw39Hb/c0RJQBoCNCLdZAlzk89LASwH4Ac4QQgS8zgFcA/ANAoWJb0MssAHxPRMuJaEx4W0LL7LdFotXGkkqL36VW2X13TYioEoCvANwnhDimM0QYiDILIQoAdCCiagC+JqI2Orv7vsxENBzAfiHEciLqL/MTlW2+KnOYXkKIPURUC8AcItqgs68jZfZbCz0bQH3F53oA9rhki1PsC3e7EP5/f3i7Vtmzw39Hb/ckRFQGITH/VAgxNbw50GWOIIQ4AuAnAEMQ7DL3AnApEWUhNCx6ARF9gmCXGUKIPeH/9wP4GqEh4oSW2W+CvgxAMyJqRESpAEYBmO6yTXYzHcCN4b9vBDBNsX0UEZUlokYAmgFYGu7GHSei7uHZ8BsUv/EUYfveA7BeCPGS4qsglzkt3DIHEZUHMAjABgS4zEKIcUKIekKIDISe0R+FENchwGUmoopEVDnyN4CLAKxFosvs9sywhZnkixHyjtgK4GG37YmzLJMB7AWQh9Cb+VYANQDMBbA5/H91xf4Ph8u9EYqZbwBdwpVnK4DXEY4A9to/AL0R6j6uBrAy/O/igJe5HYDfwmVeC+Cx8PbAljmq/P1R7OUS2DIj5Hm3KvxvXUSbEl1mDv1nGIYJCH4bcmEYhmE0YEFnGIYJCCzoDMMwAYEFnWEYJiCwoDMMwwQEFnSGYZiAwILOMAwTEP4fIR7ERiy43KkAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.plot(np.hstack(q6s))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Adding a clustering condition"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We will now modify the above function to also find clusters which satisfy particular $q_6$ value. But first, for a single file."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "sys = pc.System()\n",
    "sys.read_inputfile(files[0])\n",
    "sys.find_neighbors(method=\"cutoff\", cutoff=0)\n",
    "sys.calculate_q(6)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now a clustering algorithm can be applied on top using the `cluster_atoms` method. `cluster_atoms` uses a `condition as argument` which should give a True/False value for each atom. Lets define a condition."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [],
   "source": [
    "def condition(atom):\n",
    "    return atom.get_q(6) > 0.5"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The above function returns `True` for any atom which has a $q_6$ value greater than 0.5 and `False` otherwise. Now we can call the `cluster_atoms` method."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "16"
      ]
     },
     "execution_count": 16,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "sys.cluster_atoms(condition)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The method returns 16, which here is the size of the largest cluster of atoms which have $q_6$ value of 0.5 or higher. If information about all clusters are required, that can also be accessed."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [],
   "source": [
    "atoms = sys.atoms"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "`atom.cluster` gives the number of the cluster that each atom belongs to. If the value is -1, the atom does not belong to any cluster, that is, the clustering condition was not met."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [],
   "source": [
    "clusters = [atom.cluster for atom in atoms if atom.cluster != -1]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we can see how many unique clusters are there, and what their sizes are."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [],
   "source": [
    "unique_clusters, counts = np.unique(clusters, return_counts=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "`counts` contain all the necessary information. `len(counts)` will give the number of unique clusters."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Text(0.5, 0, 'Cluster ID')"
      ]
     },
     "execution_count": 28,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAX4AAAEGCAYAAABiq/5QAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAWd0lEQVR4nO3de5RlZXnn8e+vAUUuCgwFQaRsYCkuRUUsMhozihgSEhDiDcOMLLwkzWSJAY0XiBrMGEc0akTN0vQsEbJC8IqIogISgUxWgjQdkLsodggXaRQNiMjFfuaPs3soiqrqXZd9DlX7+1mrVp29zz77fTa7++mX97z7eVNVSJL6Y8WoA5AkDZeJX5J6xsQvST1j4peknjHxS1LPbD7qANrYcccda+XKlaMOQ5KWlMsuu+zHVTU2df+SSPwrV65kzZo1ow5DkpaUJP8+3X6HeiSpZ0z8ktQzJn5J6hkTvyT1jIlfknrGxC9JPdNZ4k9ySpL1Sa6asv9NSa5PcnWSD3bVviRpel32+E8FDpq8I8mLgcOAZ1XVM4APddi+JGkanSX+qroYuHPK7j8GTqqq+5pj1nfVviRpesN+cvepwH9L8j7gl8Bbq+rS6Q5MsgpYBTA+Pj68CKWeWXn8OZs8Zt1JBw8hEg3LsL/c3RzYHnge8Dbg80ky3YFVtbqqJqpqYmzsEaUmJEnzNOzEfzNwZg18B9gA7DjkGCSp14ad+M8CDgBI8lTgMcCPhxyDJPVaZ2P8Sc4A9gd2THIzcCJwCnBKM8XzfuCocrV3SRqqzhJ/VR0xw1uv6apNSdKm+eSuJPWMiV+SesbEL0k9Y+KXpJ4x8UtSz5j4JalnTPyS1DMmfknqGRO/JPWMiV+SesbEL0k9Y+KXpJ4x8UtSz5j4JalnTPyS1DMmfknqmc4Sf5JTkqxvVtua+t5bk1QS19uVpCHrssd/KnDQ1J1JdgMOBG7qsG1J0gw6S/xVdTFw5zRv/TXwdsC1diVpBIY6xp/kUOCWqrpimO1Kkh7S2WLrUyXZCngn8Nstj18FrAIYHx/vMDJJ6pdh9vj3BHYHrkiyDngSsDbJr013cFWtrqqJqpoYGxsbYpiStLwNrcdfVVcCO23cbpL/RFX9eFgxSJK6nc55BvAvwF5Jbk7yhq7akiS111mPv6qO2MT7K7tqW5I0s1l7/ElWJDl8WMFIkro3a+Kvqg3AMUOKRZI0BG3G+M9vSizslmSHjT+dRyZJ6kSbMf7XN7/fOGlfAXssfjiSpK5tMvFX1e7DCESSNBybHOpJslWSdyVZ3Ww/Jckh3YcmSepCmzH+zwD3A7/RbN8M/GVnEUmSOtUm8e9ZVR8EHgCoqnuBdBqVJKkzbRL//UkeR1NGOcmewH2dRiVJ6kybWT3vAb4J7JbkdOAFwOu6DEqS1J02s3rOS3IZ8DwGQzzHWlhNkpauNrN6Lqiqn1TVOVX1tar6cZILhhGcJGnxzdjjT7IlsBWwY5LteegL3ccDTxxCbJKkDsw21HM0cByDJH8ZDyX+u4C/6TguSVJHZkz8VXUycHKSN1XVx4cYkySpQ22mc/4oybYAzRO8ZybZt+O4JEkdaZP4311Vdyf5TeB3gNOAT3YbliSpK20S/6+a3wcDn6yqrwCP2dSHkpySZH2Sqybt+6sk1yX5bpIvJ9lufmFLkuarTeK/JcnfAocDX0/y2JafOxU4aMq+84G9q+pZwPeAE+YQqyRpEbRJ4IcD5wIHVdXPgB2At23qQ1V1MXDnlH3nVdWDzea/Ak+aW7iSpIVqk/h3BNYA9yUZB7YArluEtl8PfGOmN5OsSrImyZo77rhjEZqTJEG7Wj3nMCjQFmBLYHfgeuAZ8200yTuBB4HTZzqmqlYDqwEmJiZqvm1Jkh6uTa2eZ07ebqZyHj3fBpMcBRwCvKSqTOiSNGRtevwPU1Vrk+w3n8aSHAS8A3hRVf1iPueQJC3MJhN/krdM2lwB7AtsctA9yRnA/gxq/dwMnMhgFs9jgfOTAPxrVf3PuYctSZqvNj3+bSe9fpDBmP+XNvWhqjpimt2fbhmXJKkjbcb4/2IYgUiShmO2ssxfpVlucTpVdWgnEUmSOjVbj/9DQ4tCkjQ0s5VlvgggydbAvVW1odnejMEXtJKkJajNk7sXMFiJa6PHAd/qJhxJUtfaJP4tq+rnGzea11vNcrwk6VGsTeK/Z/LCK0meC9zbXUiSpC61mcd/HPCFJLc227sAr+4uJElSl9rM4780ydOAvRgUaruuqh7oPDJJUida1eppEv1VmzxQkvSo12aMX5K0jJj4JalnWg31JNkVePLk45ulFSVJS0ybsswfYDCL5xrgV83uAkz8krQEtenx/z6wV1Xd13UwkqTutRnjv5HBAuuSpGWgTY//F8DlSS4A/n+vv6r+pLOoJEmdaZP4z25+5iTJKQwWVV9fVXs3+3YAPgesBNYBh1fVT+d6bknS/LV5cve0eZ77VOATwN9N2nc8cEFVnZTk+Gb7HfM8vyRpHmZbgevzVXV4kiuZZiWuqnrWbCeuqouTrJyy+zAGC7ADnAZciIlfkoZqth7/sc3vQxaxvZ2r6jaAqrotyU4zHZhkFbAKYHx8fBFDkKR+m20Fro0J+t+HF87D2l8NrAaYmJiYce1fSdLcDLtkw+1JdgFofq8fcvuS1HvDTvxnA0c1r48CvjLk9iWp9+aU+JNsn2TWL3UnHXsG8C/AXkluTvIG4CTgwCQ3AAc225KkIWpTq+dC4NDm2MuBO5JcVFVvme1zVXXEDG+9ZK5BSpIWT5se/xOq6i7g5cBnquq5wG91G5YkqSttEv/mzRexhwNf6zgeSVLH2iT+/wWcC3y/WX93D+CGbsOSJHWlTcmGLwBfmLR9I/CKLoOSJHWnzZe7uwNvYlBYbfIKXId2F5YkqSttqnOeBXwa+CqwodtwJElda5P4f1lVH+s8EknSULRJ/CcnORE4j4cvxLK2s6gkSZ1pk/ifCRwJHMBDQz3VbEuSlpg2if9lwB5VdX/XwUiSutdmHv8VwHZdByJJGo42Pf6dgeuSXMrDx/idzilJS1CbxH9i51FIkoamzZO7FyXZGdiv2fWdqnIBFUlaojY5xp/kcOA7wKsYFGq7JMkruw5MktSNNkM97wT229jLTzIGfAv4YpeBSZK60WZWz4opQzs/afk5SdKjUJse/zeTnAuc0Wy/GvjGQhpN8mbgDxk8CHYl8Lqq+uVCzilJameTPfeqehvwt8CzgGcDq6vq7fNtMMmuwJ8AE1W1N7AZ8AfzPZ8kaW7alGX+QFW9Azhzmn0LafdxSR4AtgJuXcC5JElz0Gao50BgapL/3Wn2tVJVtyT5EHATcC9wXlWdN/W4JKuAVQDj4+PzaUpSB1Yef84mj1l30sFDiETzNeNQT5I/TnIlsFeS7076+SHw3fk2mGR74DBgd+CJwNZJXjP1uKpaXVUTVTUxNjY23+YkSVPM1uP/BwZf4r4fOH7S/rur6s4FtPlbwA+r6g6AJGcCvwH8/QLOKUlqacbEX1X/CfwncARAkp2ALYFtkmxTVTfNs82bgOcl2YrBUM9LgDXzPJckaY7aPLn70iQ3AD8ELgLWsYDpnFV1CYOHv9YymMq5Alg93/NJkuamzYNYfwk8D/heVe3OoIf+zwtptKpOrKqnVdXeVXVkVd236U9JkhZDm8T/QFX9BFiRZEVVfRvYp+O4JEkdaTOd82dJtgEuBk5Psh54sNuwJEldadPjPwz4BfBm4JvAD4CXdhmUJKk7berx39O83ACc1m04kqSuWWVTknrGxC9JPTNbyYYLmt8fGF44kqSuzTbGv0uSFwGHJvkskMlvVtXaTiOTJHVitsT/5wxq9DwJ+MiU9wo4oKugJEndma1WzxeBLyZ5d1W9d4gxSZI61GY653uTHAq8sNl1YVV9rduwJEldaVOk7f3AscA1zc+xzT5J0hLUpmTDwcA+VbUBIMlpwL8BJ3QZmCSpG23n8W836fUTughEkjQcbXr87wf+Lcm3GUzpfCH29iVpyWrz5e4ZSS4E9mOQ+N9RVT/qOjBJUjfa9PipqtuAszuORZI0BCOp1ZNkuyRfTHJdkmuTPH8UcUhSH7Xq8XfgZOCbVfXKJI8BthpRHJLUO7P2+JOsSHLVYjaY5PEMviD+NEBV3V9VP1vMNiRJM5u1x19VG5JckWS8qm5apDb3AO4APpPk2cBlwLGTFnwBIMkqYBXA+Pj4vBtbefw5rY5bd9LB825DkpaSNmP8uwBXJ7kgydkbfxbQ5ubAvsAnq+o5wD0MisE9TFWtrqqJqpoYGxtbQHOSpMnajPH/xSK3eTNwc1Vd0mx/kWkSvySpG23m8V+U5MnAU6rqW0m2Ajabb4NV9aMk/5Fkr6q6HngJgxpAkqQh2GTiT/JHDMbadwD2BHYFPsUgYc/Xm4DTmxk9NwKvW8C5JElz0Gao543ArwOXAFTVDUl2WkijVXU5MLGQc0iS5qfNl7v3VdX9GzeSbM5gBS5J0hLUJvFflOTPgMclORD4AvDVbsOSJHWlTeI/nsG8+yuBo4GvA+/qMihJUnfazOrZ0Cy+cgmDIZ7rq8qhHklaotrM6jmYwSyeHzAoy7x7kqOr6htdBydJWnxtZvV8GHhxVX0fIMmewDmAiV+SlqA2Y/zrNyb9xo3A+o7ikSR1bMYef5KXNy+vTvJ14PMMxvhfBVw6hNgkSR2YbajnpZNe3w68qHl9B7B9ZxFJkjo1Y+KvKssoSNIy1GZWz+4MauusnHx8VR3aXViSpK60mdVzFoPVsr4KbOg2HElS19ok/l9W1cc6j0SSNBRtEv/JSU4EzgPu27izqtZ2FpUkqTNtEv8zgSOBA3hoqKeabUnSEtMm8b8M2GNyaWZJ0tLV5sndK4Dtug5EkjQcbXr8OwPXJbmUh4/xL2g6Z5LNgDXALVV1yELOJUlqr03iP7Gjto8FrgUe39H5JUnTaFOP/6LFbjTJk4CDgfcBb1ns80uSZtbmyd27eWiN3ccAWwD3VNVCeuofBd4ObDtLu6uAVQDj4+MLaErSUrHy+HNaHbfupIM7jmR52+SXu1W1bVU9vvnZEngF8In5NpjkEAalni/bRLurq2qiqibGxsbm25wkaYo2s3oepqrOYmFz+F8AHJpkHfBZ4IAkf7+A80mS5qDNUM/LJ22uACZ4aOhnzqrqBOCE5tz7A2+tqtfM93ySpLlpM6tncl3+B4F1wGGdRCNJ6lybWT2d1eWvqguBC7s6vyTpkWZbevHPZ/lcVdV7O4hHktSx2Xr890yzb2vgDcB/AUz8krQEzbb04oc3vk6yLYMnbV/HYCbOh2f6nCTp0W3WMf4kOzB4svZ/AKcB+1bVT4cRmCSpG7ON8f8V8HJgNfDMqvr50KKSJHVmtge4/hR4IvAu4NYkdzU/dye5azjhSZIW22xj/HN+qleS9OjX5gGuXmlTJMoCUZrKPzdaSuzVS1LPmPglqWdM/JLUMyZ+SeoZE78k9YyJX5J6xsQvST1j4peknjHxS1LPDD3xJ9ktybeTXJvk6iTHDjsGSeqzUZRseBD406pa29T5vyzJ+VV1zQhikaTeGXqPv6puq6q1zeu7gWuBXYcdhyT11UiLtCVZCTwHuGSa91YBqwDGx8eHGldfdF1YrMvztzn3Qs7fpaUcuxbPKP8cjOzL3STbAF8CjquqR9T3r6rVVTVRVRNjY2PDD1CSlqmRJP4kWzBI+qdX1ZmjiEGS+moUs3oCfBq4tqo+Muz2JanvRtHjfwFwJHBAksubn98bQRyS1EtD/3K3qv4vkGG3K0ka8MldSeoZE78k9YyJX5J6xsQvST1j4peknjHxS1LPmPglqWdGWqRtOZhrIbK5HD/XIk5dF31aykXXurxPXeviz8GjNfYuzr+Q+9rl39dRsscvST1j4peknjHxS1LPmPglqWdM/JLUMyZ+SeoZE78k9YyJX5J6xsQvST0zqsXWD0pyfZLvJzl+FDFIUl+NYrH1zYC/AX4XeDpwRJKnDzsOSeqrUfT4fx34flXdWFX3A58FDhtBHJLUS6mq4TaYvBI4qKr+sNk+EvivVXXMlONWAauazb2A6xcxjB2BHy/i+R7NvNblyWtdnhb7Wp9cVWNTd46iOmem2feIf32qajWwupMAkjVVNdHFuR9tvNblyWtdnoZ1raMY6rkZ2G3S9pOAW0cQhyT10igS/6XAU5LsnuQxwB8AZ48gDknqpaEP9VTVg0mOAc4FNgNOqaqrhxxGJ0NIj1Je6/LktS5PQ7nWoX+5K0kaLZ/claSeMfFLUs/0KvH3qVREknVJrkxyeZI1o45nsSU5Jcn6JFdN2rdDkvOT3ND83n6UMS6WGa71PUluae7v5Ul+b5QxLoYkuyX5dpJrk1yd5Nhm/7K7r7Nc61Dua2/G+JtSEd8DDmQwpfRS4IiqumakgXUkyTpgoqqW5YMvSV4I/Bz4u6rau9n3QeDOqjqp+Yd9+6p6xyjjXAwzXOt7gJ9X1YdGGdtiSrILsEtVrU2yLXAZ8PvAa1lm93WWaz2cIdzXPvX4LRWxjFTVxcCdU3YfBpzWvD6NwV+kJW+Ga112quq2qlrbvL4buBbYlWV4X2e51qHoU+LfFfiPSds3M8T/0CNQwHlJLmvKX/TBzlV1Gwz+YgE7jTierh2T5LvNUNCSH/6YLMlK4DnAJSzz+zrlWmEI97VPib9VqYhl5AVVtS+DKqhvbIYLtHx8EtgT2Ae4DfjwaMNZPEm2Ab4EHFdVd406ni5Nc61Dua99Svy9KhVRVbc2v9cDX2Yw1LXc3d6MnW4cQ10/4ng6U1W3V9WvqmoD8H9YJvc3yRYMEuHpVXVms3tZ3tfprnVY97VPib83pSKSbN18YUSSrYHfBq6a/VPLwtnAUc3ro4CvjDCWTm1MhI2XsQzub5IAnwauraqPTHpr2d3Xma51WPe1N7N6AJqpUR/loVIR7xtxSJ1IsgeDXj4MynL8w3K71iRnAPszKGN7O3AicBbweWAcuAl4VVUt+S9FZ7jW/RkMBxSwDjh64zj4UpXkN4F/Aq4ENjS7/4zB2Peyuq+zXOsRDOG+9irxS5L6NdQjScLEL0m9Y+KXpJ4x8UtSz5j4JalnTPxa9pL8WpLPJvlBkmuSfD3JU5OsnFzxco7nfG2SJy4wrtcm+UTzenJVxhuSnJnk6Qs5vzQTE7+WteZBmS8DF1bVnlX1dAbzpXde4KlfC8wp8SfZ1FKnf11V+1TVU4DPAf+YZGye8UkzMvFruXsx8EBVfWrjjqq6vKr+afJBk3vfzfbXkuyfZLMkpya5qlnf4M1JXglMAKc3PfTHJXlukouaonjnTioxcGGS/53kIuDYtkFX1eeA84D/vrDLlx5p6IutS0O2N4Na5/O1D7DrpDr421XVz5IcA7y1qtY0NVc+DhxWVXckeTXwPuD1zTm2q6oXzaPttcDTFhC7NC0TvzS7G4E9knwcOIdBL3yqvRj8A3P+YGSJzRhUVtzoc/Nse7qKstKCmfi13F0NvLLFcQ/y8KHPLQGq6qdJng38DvBGBiskvX7KZwNcXVXPn+Hc98wp4oc8B1h2y2Zq9Bzj13L3j8Bjk/zRxh1J9ksydehlHbBPkhVJdqMph5tkR2BFVX0JeDewb3P83cC2zevrgbEkz28+s0WSZywk6CSvYFBV9YyFnEeajj1+LWtVVUleBny0Wa/1lwyS/HFTDv1n4IcMqiVexWB8HQartH0mycZO0gnN71OBTyW5F3g+g/+r+FiSJzD4e/VRBv+3MRdvTvIaYOsmhgOq6o45nkPaJKtzSlLPONQjST1j4peknjHxS1LPmPglqWdM/JLUMyZ+SeoZE78k9cz/A8tonKZKR/+QAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.bar(range(len(counts)), counts)\n",
    "plt.ylabel(\"Number of atoms in cluster\")\n",
    "plt.xlabel(\"Cluster ID\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we can finally put all of these together into a single function and run it over our individual time slices."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 43,
   "metadata": {},
   "outputs": [],
   "source": [
    "def calculate_q6_cluster(file, cutoff_q6 = 0.5, format=\"lammps-dump\"):\n",
    "    sys = pc.System()\n",
    "    sys.read_inputfile(file, format=format)\n",
    "    sys.find_neighbors(method=\"cutoff\", cutoff=0)\n",
    "    sys.calculate_q(6)\n",
    "    def _condition(atom):\n",
    "        return atom.get_q(6) > cutoff_q6\n",
    "    sys.cluster_atoms(condition)\n",
    "    atoms = sys.atoms\n",
    "    clusters = [atom.cluster for atom in atoms if atom.cluster != -1]\n",
    "    unique_clusters, counts = np.unique(clusters, return_counts=True)\n",
    "    return counts"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {},
   "outputs": [],
   "source": [
    "q6clusters = [calculate_q6_cluster(file) for file in files]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can plot the number of clusters for each slice"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Text(0, 0.5, 'number of unique clusters')"
      ]
     },
     "execution_count": 33,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAX4AAAEGCAYAAABiq/5QAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAgAElEQVR4nO3deXyU9bX48c/JBkmAhISwhSUJIrJvkYRNEResdatLW6toBfSqvb3aWlptf7f7rba23t7aFUGl1bpU0bZWDYoba9i3sLgkbIEsEAIhZM/5/TETGpCQCcwzzyzn/Xo9L2aezMxzHJMzz3yf7zlfUVWMMcZEjii3AzDGGBNYlviNMSbCWOI3xpgIY4nfGGMijCV+Y4yJMDFuB+CLHj16aEZGhtthGGNMSFm3bt1BVU07dX9IJP6MjAzWrl3rdhjGGBNSRGT36fbbUI8xxkQYS/zGGBNhLPEbY0yEscRvjDERxhK/McZEGMdm9YhIZ+BDoJP3OC+r6g9EJAV4EcgAdgFfVNXDTsVhjDEd8dqGYh7L28n+yhr6Jsczd8YQrh+b7nZYfuXkGX8dMF1VRwNjgCtFJBd4CFiiqoOBJd77xhjjutc2FPPwoi0UV9agQHFlDQ8v2sJrG4rdDs2vHEv86nHMezfWuylwHbDQu38hcL1TMRhjTEc8lreTmoamk/bVNDTxWN5OlyJyhqNj/CISLSIbgTLgbVXNB3qp6gEA778923ju3SKyVkTWlpeXOxmmMcYAsL+ypkP7Q5WjiV9Vm1R1DNAPmCAiIzrw3Hmqmq2q2Wlpn6k4NsYYv+ubHH/a/VFRwtPLi6iuawxwRM4IyKweVa0E3geuBEpFpA+A99+yQMRgjDHtmTtjCCIn74uLFvp3j+dH/9zGxEeW8OibOyg5UutOgH7iWOIXkTQRSfbejgcuA3YA/wDu8D7sDuDvTsVgjDEdcWFmCqrQrXMMAqQnx/OLm0bz/txLWHTfJKYM7sG8Dz9l6i/e5ZsvbmTb/qNuh3xWnGzS1gdYKCLReD5gXlLV10VkJfCSiMwG9gA3OxiDMcb4bHFBCQCvfW0yWWldTvrZuAHd+f2t49lz6DhPLS/ipbV7WbShmCnn9WDO1EwuPj8NOfXrQpCSUFhsPTs7W607pzHGaV+et5KK6noWf+Pidh975HgDz63ezcIVuyg9Wsf5vbowZ0oW143tS6eY6ABE2z4RWaeq2afut8pdY4wBKqrrWV1UwYzhvX16fFJCLPdNO4+l357Or24eTZQI335lM5MffY8nlnzM4ep6hyM+eyHRj98YY5z2zrZSmhWfE3+LuJgobhzfjxvGpbP8k0M8ubSQX739Eb97/xNuHt+f2VMyyeiR6FDUZ8cSvzHGAHkFJaQnxzO8b7ezer6IMGVwD6YM7sFHpVXMX1rIi2v28mz+bi4f2ou7Lsoie2D3oLgOYInfGBPxjtU1svSTg9yWM9Avifn8Xl35xU2j+daMIfx5xW6ezd/N4m2ljO6fzF1TM7lyeG9iot0babcxfmNMxHt/Zxn1jc1cOaJjwzzt6dm1M9+aMYQVD03nJ9cN58jxev7zrxuY9sv3WbCsiGMuFYRZ4jfGRLy8glJSE+MYP7C7I6+fEBfDzIkZLHlwGn+aOZ4+SZ35yeuegrBH3tzOgSOBbQlhQz3GmIhW19jEezvKuHpUH6KjnB1/j44SZgzvzYzhvdm4t5Inlxby5IeFLFhaxDWj+zJnaibD+yYBzraHtsRvjIloKz45xLG6Rmb4eZinPWP6J/O7r4xjb8Vxnl6+ixfX7OHVDcVMGpTKiL7d+Muq3dQ0NAP/bg8N+CX521CPMSai5RWU0KVTDJMGpbpy/P4pCXz/mmGsePhSHvrcBRSWVzNvadGJpN/Cn+2hLfEbYyJWU7Py9rZSLrmgp+vVtknxsdxz8SA+/PYlbT7GX+2hLfEbYyLW2l0VHKqu58oOFm05KS4mivQ22kO31Ta6oyzxG2MiVl5BKXExUUwbElxrfsydMYT42JO/gcTHRjN3xhC/vL5d3DXGRCRVJa+ghKnn9SCxU3ClwpYLuDarxxhj/Khg/1GKK2u4/7LBbodyWtePTfdboj+VDfUYYyJSXkEJUQKXDe3ldigBZ4nfGBOR8gpKmJCZQkpinNuhBJwlfmNMxCksP8ZHpceCajZPIFniN8ZEnLyCUgCusMRvjDGRIa+ghFH9kvw2Lz7UWOI3xkSUkiO1bNxb2eGVtsKJJX5jTERZvK0E6PgSi+HEscQvIv1F5D0R2S4iBSJyv3f/GBFZJSIbRWStiExwKgZjjDlVXkEJg9ISOa9nF7dDcY2TBVyNwIOqul5EugLrRORt4BfAj1T1TRG5ynt/moNxGGPOwMm+78HmcHU9qworuOfiLLdDcZVjiV9VDwAHvLerRGQ7kA4o0LKacRKw36kYjDFn9tqGYh5etIWahibA/33fg82SHWU0NWtED/NAgMb4RSQDGAvkAw8Aj4nIXuCXwMOBiMEY81mP5e08kfRb+LPve7DJKyihb1JnRqYnuR2KqxxP/CLSBXgFeEBVjwL3At9Q1f7AN4AFbTzvbu81gLXl5eVOh2lMRGqrv7u/+r4Hk+P1jXz4UTlXDO+NiLNLLAY7RxO/iMTiSfrPqeoi7+47gJbbfwNOe3FXVeeparaqZqelBVfLVGPCRVvz2MNxfvsHO8upa2yO+GEecHZWj+A5m9+uqo+3+tF+4GLv7enAx07FYIw5s69NH/SZfZ1iovzW9z2YvFVQQveEWC7M6O52KK5zclbPZGAmsEVENnr3fRe4C/g/EYkBaoG7HYzBGHMGVTWNAKR17cTBqjoUGNs/Oewu7NY3NvPujjI+N6I3MdFWvtShxC8i3YH+qrq5vceq6jKgrYG08R05rjHG/xqamnlmxS4mZqXy/N25AHzv1S38bd0+Dh6ro0eXTi5H6D8rCw9RVdtowzxe7X70icj7ItJNRFKATcDTIvJ4e88zxgS3N7Yc4MCRWuZMzTyxb9aUTOobm3l21W4XI/O/t7aWkBgXzeTzergdSlDw5TtPknc2zg3A06o6HrjM2bCMMU5SVRYsKyKrRyKXDOl5Yv+gtC5cekFP/rJyN7WnTPMMVU3NytvbSpl2QU86n7KObaTyJfHHiEgf4IvA6w7HY4wJgLW7D7N53xFmTckkKurkEdnZUzM5VF3PPzaGR23lhj2HOXiszoZ5WvEl8f8IyAM+UdU1IpKFzcQxJqTNX1pIckIsN47r95mfTcxKZVifbsxfVoiquhCdf721tYS46CguGWLTwlucMfGLSDSei7mjVPU+AFUtVNUbAxKdMcbvdh+qZvG2Um7NGUB83GeHPkSE2VMy+aj0GEs/PuhChP6jquRtK2Hyeal07RzrdjhB44yJX1WbgGsDFIsxJgCeXr6LmCjh9okZbT7mmtF96dm1E/OXFQUuMAdsP1DF3ooaG+Y5hS9DPStE5LciMlVExrVsjkdmjPG7I8cbeGntXq4Z3Zde3Tq3+bi4mCjumJTBhx+Vs7OkKoAR+tdbBSVECVw2rJfboQQVXxL/JGA48GPgV97tl04GZYxxxvNr9nC8vonZUzLbfexXJgygc2wUT4XwWf/ighKyM1LCqibBH9pN/Kp6yWm26YEIzhjjPw1NzTyzfBeTBqUyvG/73Sm7J8Zx0/h+vLqxmPKqugBE6F+7Dlazo6TKhnlOw5cCrl4iskBE3vTeHyYis50PzRjjT29sOUDJ0ZMLttoza3LoFnTlFXiWWLzChnk+w5ehnmfwTOfs673/EZ6e+saYEHGiYCstkWnn92z/CV5ZaV24bGhPnl0VegVdeQUljEjvRv+UBLdDCTq+JP4eqvoS0Aygqo1AaP0GGBPh1uzyFmxN/mzBVntmTfEUdP19Y7FD0flf2dFa1u+pZMYwG+Y5HV8Sf7WIpOJZMhERyQWOOBqVMcavzlSw1Z4TBV1Li0KmoCtvWykAM0ZY4j8dXxL/N4F/AINEZDnwZ+C/HI3KGOM3uw5W8/b2Um7LGXjagq32iAhzpmbycdkxPgyRgq7FBSVk9UhkcM8ubocSlHxJ/AV4Fk6ZBPwHnqmdO5wMyhjjP08vL/IWbA0869e4epS3oGtpoR8jc8aR4w2s/PSQLbF4Br4k/pWq2qiqBaq6VVUbgJVOB2aMOXeegq19XDs6nZ5nKNhqT0tB19KPDwZ9QdeSHaU0Niszhttsnra0mfhFpLeIjAfiRWRsq6rdaYBdJjcmBPx19R5qGnwr2GrPrTkDiI+NZsGy4D7rzysooXe3zozul+x2KEHrTCtwzQC+CvTDU63b8p2pCs8SisaYINbQ1MzCFbuYfF4qw/p2O+fXS07wFHS9uGYvc2dcQFrX4KuGralv4oOPyvlidv8Oz16KJG2e8avqQlW9BPiqqk5vVbV7raouCmCMxpizcKJga0qW317zzskZNDQHb0HXBx+VU9vQbNW67fBljL+fd+lFEZH5IrJeRK5wPDJjzFlTVZ5cWkhWWiIXn++/PvRZ3hW6grWga3FBCckJsUzITHE7lKDmS+Kf5V168QqgJ3An8KijURljzsnqogq2Fh9l9mlW2DpXs6dkcai6ntc2BFdBV0NTM+9sL+XSC3oRG+1Laotcvrw7Lb81V+FZc3dTq33GmCA0f1kR3RNiuWFsxwu22pOblcLwvt2Yvyy4CrpWFR7iaG2jzebxgS+Jf52ILMaT+PNEpCve9g1nIiL9ReQ9EdkuIgUicn+rn31dRHZ69//i7MNv22sbipn86LtkPvQvJj/6btCdnRjjlKKD1byzvZTbcs+uYKs9LQVdn5Qd44OPyv3++mcrr6CE+NhoLvLj0Fa4OtOsnhazgTFAoaoe97ZvuNOH5zUCD6rqeu+HxToReRvoBVwHjFLVOhHxvWOUj17bUMzDi7ZQ4x2DLK6s4eFFWwC4fmy6vw9nTFB5enkRsVFRzDyHgq32fH5kXx59cwcLlhUxbYjf/4Q7rLlZWVxQyrQhaXSO9f+HXbjx5Yx/CtAFGCUiF+Gp3G13gqyqHlDV9d7bVcB2IB24F3hUVeu8Pys7y9jb9FjezhNJv0VNQxOP5e3096GMCSqVx+v529p9XDumLz27nn3BVntaF3TtKDnq2HF8tWFvJWVVdTabx0e+JP65rbb/Bv4J/LAjBxGRDGAskA+cD0wVkXwR+UBELmzjOXeLyFoRWVte3rGvk/srazq035hw4c+CrfZ8ZYK3oGup+yt0LS4oITZauOQC9799hAJfVuC6ptV2OTACKPX1ACLSBXgFeMA7OygG6A7k4vkweUlO01BDVeeparaqZqeldWzMrm9yfIf2GxMO6hs9BVtTzuvB0D7nXrDVnuSEOG7O7sffN+53dYUuVSWvoISJg3qQFB/rWhyh5GzmPO3Dk/zbJSKxeJL+c62KvvYBi9RjNZ4LxT3OIo42zZ0xhPhTxvniY6OZO2OIPw9jTFB5Y8sBSo/WBeRsv8WdkzNpaG7mLy4WdO0srWLXoeM2m6cD2r24KyJP4O3Fj+eDYgywyYfnCbAA2K6qj7f60WvAdOB9ETkfiAP82uu15QLuY3k7Ka6sITZaeOSGkXZh14QtVWX+skIG+blgqz2ZPRK59IJePLtqN/dNG+TKhdW8raWIwOW2xKLPfDnjXwus824rge+o6m0+PG8yMBOYLiIbvdtVwFNAlohsBV4A7lAHJgNfPzad5Q9N5+vTz6NZ4TL7pTBhLP9EwVZWwHvUzJmaSUV1Pa+6NGU6r6CE8QO6O3oxO9y0e8avqgvP5oVVdRltF3r58sHhFzmZqTzx7ies3VURFNPOjHHC/KXegq1xgf9Wm5OZwoj0bixYVsSXAtwcbW/FcbYdOMr3rhoasGOGgzO1Zd4iIptPs20Rkc2BDPJcjBuYTGy0kF9U4XYoxjiisPwYS3aUMjN3oCtDLSLCnClZnoKujwNb0JVXUAJg0zg76Exn/FcHLAoHJcTFMKpfMqsKD7kdijGOeHr5LmKjorjNwYKt9lw1so+noGtpEZcE8Jt1XkEJQ/t0Y0CqLRHSEWdqy7xbVXd7H1Pa6n4ZIdarJyczhS37jnC8vtHtUIzxq8rj9by8bh/XOVyw1Z6Wgq5lnxxk+4HAFHSVV9Wxdvdhm81zFny5uPs3Tu7N0+TdFzJyslJpbFbW7T7sdijG+NWJgq2pgZvC2ZaWgq6nlgWmoOvtbaWowpUjbJino3xJ/DGqWt9yx3s7zrmQ/C97YHeio8SGe0xYaV2wdUFv5wu22pOUEHuioKusqtbx471VUMLA1ASG9Orq+LHCjS+Jv1xErm25IyLX4ed5905L7BTDyPQk8gvtAq8JH//ast9TsBUEZ/stWgq6nl3pbEHX0doGVn56kBnDe3Oawn/TDl8S/z3Ad0Vkj4jsAb4D3O1sWP6Xk5XCpn2V1NQH36pBJrDCoWW3qjJ/aRHn9ezCxYODpw1xZo9ELhvai784vELXezvKaGhSm81zlnzp1fOpquYCw4DhqjpJVT91PjT/ys1KpaFJ2bDHxvkjWUvL7uLKGpR/t+wOteS/qrCCgv3OrLB1ruZMyeTw8QYWrXfuPX1rawk9u3ZibP92GwWb0/C5V4+qHvO2Vw5J2QO7EyXYOH+EC5eW3QuWFZKSGMcXgrANyYTMFEamJ7FgWSHNzf5foau2oYn3d5Zz+bBeQfehFyoiZmHKrp1jGZGexCor5IpYqhoWLbsLy4/xzvYybnOpYKs9LSt0fVpe7cgKXUs/PkhNQ5PN5jkHEZP4wTOff+PeSkfHHk3waWxq5vXN+7n+9yto6/yzV7fQ6fPy1PIi4qKjmJnrXsFWe64a2Yfe3Tozf1mh31/7ra0ldOscQ25Wqt9fO1K0m/hFJEFE/ltEnvTeHywiIVnVm5uVSn1jMxv3VrodigmAY3WNnqUBf/k+//nXDVQer+fGcel0jv3sr31tQyOflAX/SGbrgq20rp3cDqdNsdGegq7lnxzya0FXY1MzS3aUcunQXsRGR9R5q1/58s49DdQBE7339wE/dSwiB2VnpCA2zh/2Dhyp4ZE3tzPxkSX85PVt9O7WmT/eNp53H5zGr744hkdvGEV6cjwCpCfH8+AV5xMTHc2Nf1jJut3BPRT4XP4eahuag2oKZ1tOrNDlx4Ku1UUVVB5vsNk858iXxdYHqeqXROQWAFWtOd2KWaEgKT6WYX262Xz+MFWw/wjzlxbxz037aVblcyP6MGdqJmMHdD/pcdePTf/M2gzXj0nn9qdWc+v8fH57y7igbOPdUrA1dXBwFGy1Jykhli9m9+Ovq/fw7RlD6OmH4bS3CkroHBsV0DUHwpEvZ/z1IhKPdzEWERmE5xtASMrJTGX9nsPUNdo4fzhoblbe21HGV55cxed/s4y8ghJmThzIB3Mv4Xe3jvtM0m9L/5QEXr5nIkN6deXuv6zlhdV7HI68417fvJ+yqsCusHWu7pycSWOz+mWFruZmZXFBKRcNTiM+LvguaocSXxL/D4C3gP4i8hywBPi2o1E5KDcrhbrGZjbvO+J2KOYc1DY08eKaPVzx6w+585k1FJZX89DnLmDlw5fyg2uG0z+l490aU7t04q935XLR+Wk8tGgLv1nyMQ6sEXRWWgq2BvfsElJnuxk9Erl8qGeFrnMtntxcfISSo7U2m8cPfFmI5W0RWY9ncXQB7lfVkGrZ0NqETO84/6eHuDAjxe1wTAdVVNfz7Krd/HnlLg4eq2don248/sXRXD2qL3Ex536xL7FTDE/ens1Dr2zh8bc/ovRoLT++bgTRLs8XX1l4iG0HjvLoDSNDrkXBnKlZLN5WyqIN+7g15+xnIr21tYSYKOHSC4JvGC7U+LLm7kXemy1THoaJCKr6oXNhOSc5IY4hvbqSX1TB190OxvissPwYC5YV8cr6fdQ2NDNtSBp3T81i4qBUvyfC2OgofnnzKHp268Qf3v+Ug8fq+L8vj3V1zvyCpUWkJsaF5LrRF2Z0Z1S/JBYsK+KWCwecVdGVqrK4oITcrFSSEmIdiDKy+HJxd26r252BCXjW353uSEQBkJuVyotr9tLQ1GxTwoKYqrK6qIInlxaxZEcpsVFRfGFsOnOmZjLY4Y6MIsJ3rryAnl078ePXtzFzQT7zb7/QlaTzafkxluwo4/5LBwdlwVZ7RITZUzK5/4WNvP9RGdPP4oz9k7JjFB6s5s4Qur4RzHzp1XNNq+1yYARQ6nxozsnNSqGmocnG+YNUY1Mz/9y0n+t/t5wvzVvFut0VfP2S81j+0HR+ftMox5N+a3dOzuSJW8ayae8Rbv7TCg4cCXyF79Pegq3bgrhgqz1XjexDn6TOZz21862tniUWrwjC2VahyJcz/lPtw5P8Q9aETE/F36rCQ4wf6NusD+O8qtoGXlyzl6eX76K4sobMHon89PoR3Diun6uzOK4e1ZeUxDju/vM6bvj9Cv48a0LAPnwOV3sKtq4fG9wFW+1pKeh69M0dbNt/lGF9OzYdNW9bCWMHJIdUhXUw86Vy9wkR+Y13+y2wFNjkw/P6i8h7IrJdRApE5P5Tfv4tEVER6XH24Z+dlMQ4zu/VxRZgd8HpWiLvr6zhZ29sZ9Ij7/LTf20nPTmeeTPHs+SbF3Nb7sCgmLo3aVAPXvyPXBqblZv+uJK1uwLzu/PX1d6CrSlZATmek265cAAJcR0v6Np3+Dhbi49ypRVt+Y0vZ/xrW91uBJ5X1eU+PK8ReFBV14tIV2CdiLytqttEpD9wOeDaZOmczFQWrd9n4/wB1NISuaU7ZnFlDQ++tIlmVaKihM+N6M1dU7MYHaStdof3TWLRvZO4w1vo9ZtbxjpaQVrX2MQz3oKtIb1Df5UpT0FXf57L3813rvS9oCuvwDOybNW6/uPLGP/CVttzPiZ9VPWAqq733q4CtgMtUxL+F08tgGuTpHOzUqmub2JrsY3zB8rpWiI3qZLQKZoP5k7jt18ZF7RJv0X/lAT+ds9ELujTjXufXcdf8507d3l90wHKq+qYMzX0z/Zb3Dk5g8Zm5c8dWKErr6CEIb26ktEj0cHIIosvQz1bRGTzabYtIrLZl4OISAYwFsj3LuNYrKpnHC4SkbtFZK2IrC0v939r1wmZnjn8NtwTOG21Pj5e10S/7h0vuHJLapdOPH9XDhefn8Z3X93Cr9/5yO+FXqrK/GWegq2LBgd8NNQxA1MTuWJYL57N962g6+CxOtbuqmCGFW35lS9jHG/iqdy91bu9AbwMXA1c096TRaQL8ArwAJ7hn+8B32/veao6T1WzVTU7Lc3/lYppXTsxKC2RfGvYFjB9k+M7tD+YJcTFMO/2bG4a349fv/Mx3311K41NzX57/ZWferpazpmaGXIFW+2ZMzWLyuMNvLJ+X7uPfWdbKc0KM4bbbB5/8iXxT1bVb6vqFu/2EDBDVXer6hm/r4lILJ6k/5yqLgIGAZnAJhHZBfQD1ouIKx/nOVmprN112K9/sKZts6dkfGZffGw0c2cMCXwwfhAbHcVjN43ia5cM4vnVe7j3ufV+W+th/jJPwdZ1Y0KvYKs92QO7M7pfEk8tK2p3ha68ghL6dY9nWJ/gb0oXSnxJ/IkiMqXljohMAtodbPN28FwAbFfVxwG8Hxw9VTVDVTPwTA0dp6olZxX9OcrNSqWqrpFtfuwXbtp24EgtAvTq1ulES+RHbhgZktWoLUSEuTMu4EfXDued7aXcNj+fyuP15/San5Yf490dwbvC1rkSEWZNyaTwYDXvf1TW5uOqahtY/skhrhzeO+y+9bjNl1k9s4GnRCTJe78SmOXD8yYDM4EtIrLRu++7qvpGx8N0Rm7LOH9hBaP6BfdFxVB3rK6RF1bv5erRfXnilrFuh+N3d0zKoEeXTnzjxY3c/MeVLJw14ayHsJ5aVkRcTGgXbLXnqpF9ePTNHcxfWtRmJe97O8upb2q28X0H+DKrZ52qjgZGAaNVdUzLbJ12nrdMVUVVR3mfM+bUpO8983et4VvPbp3J7JFIfpGN8zvtpTV7qaprDKmWwh31+VF9WDhrAiVHarnh9yv4qLTjK3pVVNfzyvp9fGFMekgXbLUnNjqKr07KYMWnhyjYf/qZdXkFJfToEsc4H1trG9+1mfhF5Dbvv98UkW8Cc4DZre6HhZzMFFYXVdDUzlijOXtNzcpTy4u4MKM7Y4J8uua5mjgolZfumUizKjf9YQWrOzhr7K/5u0Nmha1z9eUJbRd01TY08f6OMi4f1tv1zqjh6Exn/C3j+F3b2MJCblYqR2sb/bouqDnZ4oIS9h2uCYvqU18M7dONRfdNokfXTty2IP9En5n21DU2sXDlbi46P43zA9iPyC1J8Z6Crn9u2k/p0dqTfrb8k4NU1zfZbB6HtJn4VfVP3n9/dLotcCE6KyfL5vM7bf6yIgakJHB5BDXY6tc9gZfvmcSwPt2477l1POvDClT/bCnYCuPhsFPN8q7Q9eeVu07an1dQQtdOMUwaFD41DMHElwKuNBH5rojME5GnWrZABBcIfZLiGZCSYPP5HbJ+z2HW7T7MrMkZEfeVPSUxjr/elcO0IT35f69t5fHFO9ss9PKssFXI+b26MDWMCrbaMyA1gRnDevNc/p4TBV2NTc28s72M6UN7+mVxHfNZvryrfweSgHeAf7XawkZOZgqrd1W0O6fYdNyCZUV07RzDzdn93Q7FFQlxMcybOZ4vZvfjN+9+wsOLtpy2bmTFp4fYUVLF7CnhV7DVntlTM08q6Fqz6zAV1fXWm8dBvkznTFDV7zgeiYtys1L527p97CytYqgVivjNvsPHeXPLAe6amkVip7PpAB4eYqKj+PmNo+jZtTO/fe8TDh6r44lbxp3UdXRBGBdstad1QddXJgwgr6CEuJiokFpbONT4csb/uohc5XgkLjoxzm/DPX61cMUuRIQ7JmW4HYrrRIRvzRjCT64bzpIdZdw6fxWHqz2FXp+UeQq2Zk4Mz4Kt9ogIs6dmUXiwmnd3lLG4oISLBqdF9MmC03x5Z+8HvisidUADngXXVVXD5tS4X/cE0pPjyS+q4KuTI+fCmpOqaht4YfVePj+yT0j24nHKzImeQq/7X9jIjF9/iAiUHq0DILVLnMvRuedzI3qTHB/Df/xlHU2qHG9o4pShJRUAABbbSURBVLUNxSFd1R3M2k38qhr+88rwDPe8t7MMVY24MVYnvLR2H1V1jcyJgPnoHfW5kX0o2H+E37736Un7f/avHXTtFBuRye5fmw9wrK6JJu/F78rjDTy8aAtARL4fTvNlVs9Fp9sCEVwg5WSlUFFdz8dlx9wOJeQ1NjXz9PIiJmSkWCuMNry6Yf9n9tU0NPFY3k4XonHfY3k7aTxlckUkvx9O82WoZ26r252BCcA6YLojEbkk17sOb37hoYgonnHS4m2l7Dtcw39fPcztUIJWW2sTtLU/3Nn7EVi+9Oq5ptV2OZ6F1kudDy2w+qfE0yepM6uskOuczV9ayMDUBC4bGjkFWx0VTmsT+IO9H4F1NtUR+/Ak/7AiIuRmpZJfeMjvqylFknW7D7N+TyV3Toq8gq2OmDtjCPGnzOAJ5bUJzpW9H4HV7lCPiDzBv9fGjQLGAGdcNjFU5WSm8OqGYj4tr+a8nl3cDickPRXhBVu+arlg+VjeTvZX1tA3OZ65M4ZE7IVMez8Cy5cx/rWtbjcCz/u64HqoycnyjvMXHbLEfxb2Vhznza0HuOuiyC7Y8tX1Y9MtsbVi70fg+DKdc2EgAgkGGakJ9OzaifzCCm7NCd9FMJyycMUuokT4qhVsGRPUrANSKy3j/KtsnL/DqmobeGHNXj4/qg99kuyCnDHBzBL/KXKyUiirqmPXoeNuhxJSXlyzl2NhvsKWMeHiTCtw/cX77/2BC8d9Oa3m8xvfeAq2djEh0wq2jAkFZzrjHy8iA4FZItJdRFJab4EKMNAGpSXSo0snW5ilA/IKSimurImoBUSMCWVnurj7R+AtIAtPpW7rSdnq3R92RIScrJQT4/zWt6d985d5CrYutYItY0LCmZZe/I2qDgWeUtUsVc1stYVl0m+Rm5nCgSO17K2wcvH2rNt9mA17Kpk1OdMKtowJEb60bLhXREaLyH96t1G+vLCI9BeR90Rku4gUtFwrEJHHRGSHiGwWkVdFJOgGhVvm868qsnH+9jy1rIhunWO4aXw/t0MxxvjIl+6c/wU8B/T0bs+JyNd9eO1G4EHvt4Zc4GsiMgx4GxihqqOAj4CHzzZ4pwzu2YWUxDjyC22c/0xaCra+kjPQCraMCSG+/LXOAXJUtRpARH4OrASeONOTVPUAcMB7u0pEtgPpqrq41cNWATedTeBOEhFyMj3j/KZtz3gLtu6YZMVuxoQSX+bxC9DU6n4TJ1/obf8FRDKAsUD+KT+aBbzZxnPuFpG1IrK2vLy8I4fzi5zMFIora9h32Obzn87R2gZeXLOXq61gy5iQ40vifxrIF5EfisgP8ZylL/D1ACLSBXgFeEBVj7ba/z08w0HPne55qjpPVbNVNTstLfCLLp/o22PDPaf10omCrbC+zm9MWPLl4u7jwJ1ABXAYuFNVf+3Li4tILJ6k/5yqLmq1/w7gauBWDdLeCEN6dSU5IZZ8u8D7GS0FWzmZKYzsl+R2OMaYDvLpipyqrgfWd+SFxTMBfgGw3fvh0bL/SuA7wMWqGrTjKFFRwoSMFFbZGf9nvFVQQnFlDT+4xlbYMiYUOdmrZzIwE5guIhu921XAb4GuwNvefX90MIZzkpOVyp6K4xw4YvP5W1uwrIgMK9gyJmQ5NgdPVZdx+ovAbzh1TH/LyfR0psgvrLA+4V4tBVs/vm64FWwZE6LOeMYvItEi8k6gggk2Q/t0o2vnGJvW2cqCZYUkxcdawZYxIeyMiV9Vm4DjIhKRV/Ciozzz+a1hm8feiuO8tbWEr+QMICHOCraMCVW+/PXWAltE5G2gumWnqv6XY1EFkZzMVN7ZXkbZ0Vp6duvsdjiuenq5t2BrYobboRhjzoEvif9f3i0i5WR5xvlXFVVw7ei+LkfjHk/B1h6uGd2X3kmR/QFoTKjzac1dEYkHBqjqzgDEFFSG9elG106ecf5ITvwvrt5LdX2TrbBlTBjwpUnbNcBGPL35EZExIvIPpwMLFjHRUWRndI/oFbkam5p5ZoWnYGtEekRe7jEmrPgyj/+HwASgEkBVNwIRddqXk5XKp+XVlFfVuR2KK1oKtuZMtfYMxoQDXxJ/o6oeOWVfULZZcErLfP7VETi7R1V5cqm3YOuCnm6HY4zxA18S/1YR+QoQLSKDReQJYIXDcQWVEelJJMZFR+R8/vV7DrNpbyWzp2QSZQVbxoQFXxL/14HhQB3wPHAUeMDJoIJNbHQU4zNSIrJh2/ylRSTFx3KjFWwZEzZ86c55XFW/B1wKXKKq31PVWudDCy45mSl8VHqMiup6t0MJmD2HjpNXUMKtVrBlTFjxZVbPhSKyBdiMp5Brk4iMdz604JKb1TLOHzln/U+vKCI6SrhjUobboRhj/MiXoZ4FwH2qmqGqGcDX8CzOElFGpicTHxsdMW2aj9Q08NKavVwzqi+9Irxi2Zhw40vir1LVpS13vF03q5wLKTjFxUQxfmD3iLnA++KaPVTXNzHLCraMCTttJn4RGSci44DVIvInEZkmIheLyO+B9wMWYRDJyUxhZ2kVlcfDe5y/samZZ5bvIjfLCraMCUdnumL3q1Pu/6DV7Yiax98iJysVVc98/iuG93Y7HMe8ubWE/Udq+fF1I9wOxRjjgDYTv6peEshAQsHo/kl0ioliVWH4Jn5VZf7SQjJ7JDLdCraMCUvtztETkWTgdiCj9eMjpS1za51iohk3oHtYz+dft/swm/Yd4SfXj7CCLWPClC8Xd9/Ak/S3AOtabREpJyuFbQeOcqSmwe1QHDF/aRHJCbHcOM6WmjQmXPlSldNZVb/peCQhIiczFdWPWburIuwWG999qJq8bSXcN22QFWwZE8Z8OeP/i4jcJSJ9RCSlZXM8siA1dkAycTFRYTmt8+nlu4iJEm63FbaMCWu+JP564DFgJf8e5lnb3pNEpL+IvCci20WkQETu9+5PEZG3ReRj77/dz+U/INA6x0Yzpn9y2K3De6SmgZfW7uWa0VawZUy48yXxfxM4z1u5m+ndfGnM3gg8qKpDgVzgayIyDHgIWKKqg4El3vshJTczha3FR6iqDZ9x/hfX7OG4rbBlTETwJfEXAMc7+sKqekBV13tvVwHbgXTgOmCh92ELges7+tpuy8lKpVlh7a7DbofiFw3egq2JWakM72sFW8aEO1+u4DUBG0XkPTytmYGOTecUkQxgLJAP9FLVA97XOCAip50sLiJ3A3cDDBgwwNdDBcS4Ad2JjRZWFR3ikjCY695SsPWT661gy5hI4Evif827nRUR6QK8AjygqkdFfJsbrqrzgHkA2dnZQVUpHB8Xzeh+yeSHQcO2loKtrB6JXDIk9D/EjDHtazfxq+rC9h7TFhGJxZP0n1PVRd7dpSLSx3u23wcoO9vXd1NOVgp//KCQ6rpGEjuF7tTHtbsPs3nfEX5qBVvGRAxf+vEXiUjhqZsPzxM8LZ23q+rjrX70D+AO7+07gL+fTeBuy8lMpalZWbs7tMf55y8t9BZs2QpbxkQKX05Vs1vd7gzcDPgyj38yMBPP4i0bvfu+CzwKvCQis4E93tcLOeMHdicmSsgvPMTF56e5Hc5Z2X2omsXbSvnatPOIj4t2OxxjTID4MtRzaqXSr0VkGfD9dp63DGhr7OBS38ILXomdYhjZLymk5/P/u2BroNuhGGMCyJcmbeNa3Y3C8w2gq2MRhZCczFQWLCvkeH1jyLU4aF2w1dMKtoyJKL5kq9Z9+RuBXcAXHYkmxHgu8H7K+t2VTBncw+1wOuSF1VawZUyk8mWox/rytyF7YHeio4T8okMhlfgbmpp5ZsUuJg2ygi1jIpEvQz2dgBv5bD/+HzsXVmjo2jmWEX27hdx8/je2HODAkVr+5wtWsGVMJPKlZcPf8bRZaASqW20GT/uGjXsrqW1ocjsUn6gqC5YVkZWWyLTzrWDLmEjkyxh/P1W90vFIQlRuVgrzPixk/Z7DTBoU/MM9a3Z5Crb+5wtWsGVMpPLljH+FiIx0PJIQlZ2RQpQQMsM985cW0j0hlhvGWsGWMZHKl8Q/BVgnIjtFZLOIbBGRzU4HFiq6dY5lWN9uQb8O72sbisn52Tss3lZKY5OSV1DidkjGGJf4MtTzOcejCHE5mak8u2o3tQ1NdI4NvgrY1zYU8/CiLdR4r0NU1TXy8KItAFw/1tbWNSbStHvGr6q7T7cFIrhQkZuVSl1jM5v2Vrodymk9lrfzRNJvUdPQxGN5O12KyBjjJl+Gekw7JmSkIELQtm8orqw57f79bew3xoQ3S/x+kJQQywW9g2+cv7lZeeTN7W3+vG9yfACjMcYEC0v8fpKTmcK63Yepb2x2OxTAU537rb9t4k8fFDJ5UCrxsSf/r46PjWbujCEuRWeMcZMlfj/JzUqltqGZzfvcH+evrmtk9sK1LNpQzLeuOJ9n5+TwyA2jSE+OR4D05HgeuWGkXdg1JkKFVkvJIDYh07NEQX5RBdkZvixX4IyDx+qY9cwaCvYf5ec3juRLF3rWK75+bLolemMMYGf8fpOSGMeQXl1ZVejeOP+eQ8e56Q8r+Ki0inkzx59I+sYY05olfj/KyfKM8zc0BX6cf2vxEW74w3Iqaxp4bk4ulw7tFfAYjDGhwRK/H+VmpXK8voktxUcCetxlHx/kS39aSaeYaF6+ZxLjB3YP6PGNMaHFEr8fnRjnD2Dfnr9vLObOZ1bTPyWBRfdN4ryeXQJ2bGNMaLLE70c9unTivJ5dAjaff/7SQu5/YSPjBnTnxf+YSC9bQtEY4wNL/H6Wk5nCmqIKGh0c529uVn72xnZ++q/tXDWyNwtnTSApPtax4xljwoslfj/LzUqlur6Jgv1HHXn9+sZmvvnSRuZ9WMjtEwfyxC3jgrIxnDEmeDmW+EXkKREpE5GtrfaNEZFVIrJRRNaKyASnju+WnKyW+fz+H+45VtfI7IVreG3jfubOGMKPrh1OtC2mYozpICfP+J8BTl256xfAj1R1DPB97/2w0rNrZ7J6JPr9Am95VR23zFvFik8P8YubRvG1S85DxJK+MabjHEv8qvohcGr2U6Cb93YSsN+p47spJyuF1UUVNDWrX15v96FqbvrjCj4uq+LJ28fzxez+fnldY0xkCvQY/wPAYyKyF/gl8HBbDxSRu73DQWvLy8sDFqA/5GalUlXXyPYD5z7Ov2XfEW78wwqO1jTw/F25TL/ACrOMMecm0In/XuAbqtof+AawoK0Hquo8Vc1W1ey0tLSABegPOZmpAOfcvmHpx+V8eZ63MOveSYwdYIVZxphzF+jEfwewyHv7b0DYXdwF6J3UmYGpCee0MMtrG4q58+k1DEhNZNF9kxiUZoVZxhj/CHTi3w9c7L09Hfg4wMcPmJxMzzh/81mM8z/5YSEPvLiRCzNSePE/cq0wyxjjV05O53weWAkMEZF9IjIbuAv4lYhsAn4G3O3U8d2Wm5XKkZoGdpRU+fyc5mblp69v43/e2M7nR/XhmVkX0q2zFWYZY/zLsX78qnpLGz8a79Qxg0lOlmecP7/oEMP6dmvn0Z7CrG/9bRP/2LSfr07K4PtXDyPK5ugbYxxglbsOSU+Op1/3eJ8u8B6ra2TWM2v4x6b9fOfKC/jBNZb0jTHOsRW4HJSblcqS7aU0N2ubiby8qo47n1nN9gNV/PLm0dw0vl+AozTGRBo743dQTmYKh4838HHZsdP+vOhgNTf+YQWfllUz/45sS/rGmICwxO+g3Fbj/KfavK+Sm/6wgmN1jTx/dy6XDOkZ6PCMMRHKEr+D+nWPp29S58+M83/wUTlfnreK+LhoXr5nImP6J7sUoTEmElnid5CIkJuVyuqiClQ98/kXrd/H7GfWkJGayKJ7J5FlhVnGmACzi7sOi40RDh6rJ+vhN+jaOYajtY1MGpTKn2aOp6vN0TfGuMDO+B302oZiXtvgaUCqwNHaRqIFbhibbknfGOMaS/wOeixvJ3WNJy/B2KTwv++EbacKY0wIsMTvoP2VNR3ab4wxgWCJ30F9k+M7tN8YYwLBEr+D5s4YQvwpC6HHx0Yzd8YQlyIyxhib1eOo68emA56x/v2VNfRNjmfujCEn9htjjBss8Tvs+rHpluiNMUHFhnqMMSbCWOI3xpgIY4nfGGMijCV+Y4yJMJb4jTEmwkhL18hgJiLlwO6zfHoP4KAfwwl19n78m70XJ7P342Th8H4MVNW0U3eGROI/FyKyVlWz3Y4jWNj78W/2XpzM3o+ThfP7YUM9xhgTYSzxG2NMhImExD/P7QCCjL0f/2bvxcns/ThZ2L4fYT/Gb4wx5mSRcMZvjDGmFUv8xhgTYcI68YvIlSKyU0Q+EZGH3I7HLSLSX0TeE5HtIlIgIve7HVMwEJFoEdkgIq+7HYvbRCRZRF4WkR3e35OJbsfkFhH5hvfvZKuIPC8ind2Oyd/CNvGLSDTwO+BzwDDgFhEZ5m5UrmkEHlTVoUAu8LUIfi9aux/Y7nYQQeL/gLdU9QJgNBH6vohIOvBfQLaqjgCigS+7G5X/hW3iByYAn6hqoarWAy8A17kckytU9YCqrvfersLzRx3RiwSISD/g88B8t2Nxm4h0Ay4CFgCoar2qVroblatigHgRiQESgP0ux+N34Zz404G9re7vI8KTHYCIZABjgXx3I3Hdr4FvA81uBxIEsoBy4Gnv0Nd8EUl0Oyg3qGox8EtgD3AAOKKqi92Nyv/COfHLafZF9NxVEekCvAI8oKpH3Y7HLSJyNVCmquvcjiVIxADjgD+o6ligGojIa2Ii0h3PyEAm0BdIFJHb3I3K/8I58e8D+re6348w/MrmKxGJxZP0n1PVRW7H47LJwLUisgvPEOB0EXnW3ZBctQ/Yp6ot3wJfxvNBEIkuA4pUtVxVG4BFwCSXY/K7cE78a4DBIpIpInF4LtD8w+WYXCEigmf8druqPu52PG5T1YdVtZ+qZuD5vXhXVcPurM5XqloC7BWRId5dlwLbXAzJTXuAXBFJ8P7dXEoYXugO28XWVbVRRP4TyMNzZf4pVS1wOSy3TAZmAltEZKN333dV9Q0XYzLB5evAc96TpELgTpfjcYWq5ovIy8B6PLPhNhCGrRusZYMxxkSYcB7qMcYYcxqW+I0xJsJY4jfGmAhjid8YYyKMJX5jjIkwlvhN2BKRVBHZ6N1KRKTYe/uYiPw+QDFMa+n+KSLXRnKXWBM8wnYevzGqeggYAyAiPwSOqeovXYznH0RoEaEJLnbGbyLOKWfhPxSRhSKyWER2icgNIvILEdkiIm95W10gIuNF5AMRWScieSLS5zSve7O3h/smEfnwND//qoj81nu7l4i86n3sJhGZ5N1/m4is9n4z+ZO3vbgxfmWJ3xgYhKdF83XAs8B7qjoSqAE+703+TwA3qep44Cngf07zOt8HZqjqaODado75G+AD72PHAQUiMhT4EjBZVccATcCt5/xfZ8wpbKjHGHhTVRtEZAue9h5vefdvATKAIcAI4G1P+xai8bTsPdVy4BkReQlPc68zmQ7cDqCqTcAREZkJjAfWeI8TD5Sd/X+WMadnid8YqANQ1WYRadB/9zFpxvM3IkCBqp5xOUJVvUdEcvB8e9goImM6GIcAC1X14Q4+z5gOsaEeY9q3E0hrWYdWRGJFZPipDxKRQaqar6rfBw5yclvwUy0B7vU+L9q7CtYS4CYR6endnyIiA/3832KMJX5j2uNduvMm4OcisgnYyOl7tD/mvSi8FfgQ2HSGl70fuMQ7vLQOGK6q24D/BywWkc3A28BnLiIbc66sO6cxxkQYO+M3xpgIY4nfGGMijCV+Y4yJMJb4jTEmwljiN8aYCGOJ3xhjIowlfmOMiTD/HwpAsxQAYDB6AAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.plot(range(len(q6clusters)), [len(x) for x in q6clusters], 'o-')\n",
    "plt.xlabel(\"Time slice\")\n",
    "plt.ylabel(\"number of unique clusters\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can also plot the biggest cluster size"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Text(0, 0.5, 'Largest cluster size')"
      ]
     },
     "execution_count": 34,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYgAAAEGCAYAAAB/+QKOAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAgAElEQVR4nO3dd3hU95Xw8e/RqEsDAo2QkOio2GBMMXFBxImdgo3rOs3e9ORN1m96WdYmySa2N8UJqU6ctk7i7Juss05is7FNjEuKY3AcgwSmCgRGAgmEBJJQb3PeP2YEQh5JI2lm7pTzeZ55NHPn3rlHA5oz91fOT1QVY4wxZrgkpwMwxhgTnSxBGGOMCcgShDHGmIAsQRhjjAnIEoQxxpiAkp0OIJQ8Ho/OmzfP6TCMMSZmbN++vUlV8wI9F1cJYt68eWzbts3pMIwxJmaISM1Iz1kTkzHGmIAsQRhjjAnIEoQxxpiALEEYY4wJyBKEMcaYgOJqFJMxobSxso4Nm6uob+miMCeDdWvKuHl5kdNhGRMxliCMCWBjZR3rH9lFV98AAHUtXax/ZBeAJQmTMKyJyZgANmyuOpscBnX1DbBhc5VDERkTeZYgjAmgvqVrXNuNiUeWIIwJoDAnY1zbjYlHliCMCWDdmjLSks//88hIcbFuTZlDERkTeWFLECIyW0T+LCL7RGSPiHzSv32DiOwXkZdF5FERyRnh+CMisktEdoiIFVgyEXXz8iLWXlRw9nFORgpfu2WJdVCbhBLOK4h+4LOqeiFwOfBREVkEPA1cpKoXAweA9aO8xlWqukxVV4YxTmMC6vMqBVPSSU9J4i2XzLLkYBJO2BKEqh5X1Qr//TZgH1Ckqk+par9/t78Ds8IVgzGTUVnbwiXzplEyw82BhjanwzEm4iLSByEi84DlwIvDnvoA8McRDlPgKRHZLiIfHuW1Pywi20RkW2NjYyjCNYYTrd3UtXRxyZxplOZbgjCJKewJQkSygd8Dn1LVM0O2fx5fM9SvRzi0XFVXANfia566MtBOqvpTVV2pqivz8gKueWHMuFXUNgOwYu40ygqyaTjTQ0tnr8NRGRNZYU0QIpKCLzn8WlUfGbL9vcD1wDtVVQMdq6r1/p8ngUeBS8MZqzFDVdQ0k5acxKKZUyjJdwNwoKHd4aiMiaxwjmIS4GfAPlX99pDt1wB3ADeqaucIx2aJiHvwPvBmYHe4YjVmuIraZpYUTSU1OYkyf4KosmYmk2DCeQVRDrwbuNo/VHWHiKwFfgC4gaf9234MICKFIrLJf2w+8LyI7AT+ATyhqk+GMVZjzurpH2B33RkumTsNgJlT03GnJXPQEoRJMGEr1qeqzwMS4KlNAbYNNimt9d8/DCwNV2zGjGZ33Rl6B7wsn+NLECJCaYGbqhOWIExisZnUxgxTebaD+twczsGRTCN0mRkTlyxBGDPM9ppmZk/PYIY7/ey20vxsmjv7aGzvcTAyYyLLEoQxQ6gqFbXNrPA3Lw0a7Kg+cMJGMpnEYQnCmCHqW7tpONPzqgRRWjA41NX6IUzisARhzBAVNf7+h2EJwpOdRm5WqiUIk1AsQRgzxPaaZjJSXFww0/2q50rys20uhEkoliCMGaKytpmLZ00lxfXqP42yfDcHG9ptJJNJGJYgjPHr7htgT/0ZVsydFvD50gI37T391Ld2RzgyY5wRtolyJvZsrKxjw+Yq6lu6KMzJYN2asoRaA+HlY630e/VV/Q+Dzo1kaqPIlh41CcCuIAzgSw7rH9lFXUsXCtS1dLH+kV1srKxzOrSIOVvBdU7ARQ7PFu2zfgiTKCxBGAA2bK6iq2/gvG1dfQNs2FzlUESRV1HTzLzcTHKz0wI+PzUjhYIp6RywkhsmQViCMADUt3SNa3u88U2QaxmxeWlQaYHbriBMwrAEYQAoHKFNfaTt8ebo6S6a2ntG7KAeVJafTfXJdga8NpLJxD9LEAaAdWvKcCWdX3w3I8XFujVlDkUUWef6H8a4gsh309PvpfZ0wKVMjIkrliAMAAvzshnwKtlprrPb7rlpccKMYqqobSYr1UVZwasnyA01+LyV/jaJwBKEQVW5+7E9eLJT2br+Dfzifa8BEqd5CXwJYunsnFddRQ1XPCMbsJpMJjFYgjD8YWc922qaWbemjCnpKVw6fzrJScLz1U1OhxYRnb397DvednYFudFkpiYzZ3qmdVSbhGAJIsF19vbztU37uahoCm+9ZDYAWWnJLJ+Tw5YESRA7j7YyMMoEueFK8922/KhJCJYgEtwP/3yIE2e6ueuGxec1r5QXe9hV10prZ5+D0UXGYAf18hEmyA1XVpDN4cYOevu94QzLGMdZgkhgtac6+enfDnPTskJWzpt+3nPlxR5U4YXD8X8VUVHTzIK8LHIyU4PavzTfTb9XeaWpI8yRGeOssCUIEZktIn8WkX0iskdEPunfPl1EnhaRg/6fAa/rReQaEakSkWoRuTNccSayr2zai0uEO6+94FXPLZudQ1aqK+77IVSVyqMtXBJk8xL4EgRYyQ0T/8J5BdEPfFZVLwQuBz4qIouAO4FnVbUEeNb/+Dwi4gLuB64FFgG3+Y81IbKluonNexr46FULmTn11aOVUlxJXLYgl63VpxyILnKOnOrkdEfvmBPkhlqQl4UrSazkhol7YUsQqnpcVSv899uAfUARcBPwS/9uvwRuDnD4pUC1qh5W1V7gN/7jTAj0D3i5+7E9zJ6ewf957YIR91u1MJfDTR3UxXG5jZFWkBtNWrKL+Z4sG+pq4l5E+iBEZB6wHHgRyFfV4+BLIsCMAIcUAUeHPD7m3xbotT8sIttEZFtjY2Mow45bv36xlgMN7Xx+7SLSU1wj7re6xAMQ16OZttc2405LpsQ/vyFYZfluSxAm7oU9QYhINvB74FOqeibYwwJsC1j8RlV/qqorVXVlXl7eRMNMGM0dvXz76QOUF+eyZnH+qPuW5bvxZKfGdYKoqGlm2ZwcksaYIDdcab6bmtOddPUOjL2zMTEqrAlCRFLwJYdfq+oj/s0NIjLT//xM4GSAQ48Bs4c8ngXUhzPWRPGtp6to7+nnSzcsRmT0D0URobzYw5bqU3G5zGZ7Tz8HGtrG1bw0qDQ/G1WoPtkehsiMiQ7hHMUkwM+Afar67SFP/QF4r//+e4H/DXD4S0CJiMwXkVTgVv9xZhL21p/hv1+s5d2Xzz07Emcs5Qs9NLX3cKAh/j4Idx5twauMq4N6UGmBjWQy8S+cVxDlwLuBq0Vkh/+2FrgXeJOIHATe5H+MiBSKyCYAVe0HPgZsxte5/bCq7gljrHFvsN7S1IwUPv3G0qCPK/f3Q8TjcNftNc2I+Ib0jtfc6ZmkJifZjGoT18K2JrWqPk/gvgSANwTYvx5YO+TxJmBTeKJLPJt2neDFV07z5ZsvYmpmStDHFeVkMN+TxdbqJj64en4YI4y8itpmSmZkMzUj+PdjULIrieK8bLuCMHHNZlIngK7eAb66aR8XFLi57dI54z5+1cJc/n74FH0D8VNawutVKoNYQW40ZQVumwth4poliATwk+cOUdfSxV03Lh6znHUgq4s9dPQOsPNoSxiic8bhpnZau/omlSBK8rOpb+3mTHf816syickSRJyra+nix389xHVLZnL5gtwJvcYVC3MRgS1xNKu6osaX7CbSQT2ozN/Rb/0QJl5ZgohzX920D1VYv/bV9ZaClZOZykWFU+NqPkRFbTNTM1JY4Mma8GsMjgSLxxFexoAliLj24uFTPPHycW5/3UJmTcuc1GuVF3uoqG2mo6c/RNE5q6K2meUTmCA3VFFOBlmpLlt+1MQtSxBxasCr3PXYXgqnpnP76xZO+vVWF3vo9yr/OHI6BNE5q7WrjwMN7eOq4BpIUpJQYiU3TByzBBGnfvNSLfuOn+Fz111IRurI9ZaCtXLeNFKTk9hyMPabmXYcnXz/w6DS/GxLECZuWYKIQ62dfXxzcxWXzZ/OdUtmhuQ101NcrJw7LS4mzFXUNJMksHQCE+SGK81309Tey6n2nhBEZkx0sQQRh77zzAFau/qCqrc0HuXFHvafaKMpxj8MK2qbKc13k502+XmiZQXWUW3ilyWIOHOgoY3/9/cabrt0DosKp4T0tcuLfWU3th6K3eGuXq+yo7aFS0LQvATnhrpaM5OJR5Yg4oiqcs9je8lKdfHZN5eF/PWXFE3FnZ4c0/0QB0+209bTP6kJckPludPIyUyxkhsmLgWVIERktYi8338/T0TiqyhPnHhqbwPPVzfxmTeVMj0rNeSv70oSVi3M5fnqppgt/11R619BLkRXECJC6QwruWHi05gJQkS+BNwBrPdvSgF+Fc6gzPh19w3w5Sf2Upqfzbsunxu286wu9lDX0kXt6c6wnSOcttc0Mz0rlXm5k5sXMlRpga9oX6wmTWNGEswVxD8BNwIdcLbqanCLCZiI+dnzr3D0dBdfvH4xya7wtRyuKo7t8t8Vtc2smJMT0s77snw3bd39NJyJ7c57Y4YL5pOkV31fjRRARCZem8CExYnWbu7/czVvXpR/dh3pcFngyWLm1HS2xmBdppbOXg43drA8RP0PgwZLblg/hIk3wSSIh0XkJ0COiHwIeAZ4ILxhmfH4+pP76fcqX7huUdjPdXYZ0kNNeL2x1aRSWeufIBemBGH9ECbejJkgVPWbwO/wrS1dBnxRVe8Ld2AmONtrmnm0so4PvXY+c0LYrj6a8uJcWjr72Hv8TETOFyrba5pxJQlLZ08N6etOy0olz51mVxAm7gTTSf3vwH5VXaeq/6qqT4vIhyMQmxmD1+tbRjR/ShofeX1xxM5bvjA2+yEqapu5cKabzNTQL6RYZjWZTBwKponp48BmEblqyLbbwxSPGYffbT/Gy8daWX/thWSFYFZwsGZMSac0Pzumyn8PeJWdRye3gtxoSvPdHGxoj7lmN2NGE0yCqAOuAe4VkXX+baEbAmIm5Ex3H9/YvJ9L5k7jpmWFET//qoUeXjpymu6+gYifeyKqTrTR0TsQtgRRVpBNV98Ax5q7wvL6xjghqPGQqloLvA5YJCK/BTLGOkZEfi4iJ0Vk95Bt/yMiO/y3IyKyY4Rjj4jILv9+24L8XRLK9589yKmOXu4Kcb2lYK0u9tDd5z078SzabffHGaoSG8OV2EgmE4eCSRDbAFS1W1XfD/wFCGaa7oP4rjzOUtV3qOoyVV2Gr9P7kVGOv8q/78ogzpVQDjW284stR3j7JbNZMiu0Ha7BumzBdFxJEjPDXStrmvFkpzFr2pjfbSakZEY2YDWZTHwJZhTTh4Y9vl9VFwRx3HNAwNVlxPeV9+3AQ0HGaYb4j8f3kpHiYt01oa+3FCx3egpLZ02NmY7qcEyQG8qdnkJRToYlCBNXRkwQIvKw/+cuEXl5+G2S530t0KCqB0d4XoGnRGT7WCOmROTDIrJNRLY1NjZOMqzo96f9DfylqpFPvrEET3aao7GsLvbw8rEWWrv6HI1jLE3tPRw51Rmy+ksjKStw2/KjJq6MdgXxSf/P64EbAtwm4zZGv3ooV9UVwLXAR0XkypF2VNWfqupKVV2Zl5c3ybCiW2+/l/94fB8L8rJ4zxXznA6H8mIPXvWtfR3NBifIhav/YVBpvpvDjR30DXjDeh5jImXEBKGqx/13m4CjqloDpAFLgfqJnlBEkoFbgP8Z5dz1/p8ngUeBSyd6vnjyiy2v8EpTB/9+/SJSk52v1L58zjQyUlxRP9y1oraZ5CRhSVF4+2tK87PpHfBSc6ojrOcxJlKC+ZR5DkgXkSLgWeD9+DqgJ+qN+CbeHQv0pIhkiYh78D7wZmB3oH0Tycm2br7/p2quvmAGV5XNcDocAFKTk7h0/vSo74eoqGlmceEU0lMmvzb3aM7WZDphq8uZ+BBMghBV7cT3rf/7qvpPwJhFf0TkIeAFoExEjonIB/1P3cqw5iURKRSRTf6H+cDzIrIT+AfwhKo+GdyvE782PFlFT/8A/359+OstjcfqYg+HGjs40drtdCgB9Q142XmsJeQF+gIpnpFNkthIJhM/gpl+KyJyBfBOYPBDfszjVPW2Eba/L8C2emCt//5hfM1Yxm/n0RZ+u/0Y/3LlAuZ7oquY7qriXAC2VDfxlktmORzNq+0/3kZ3nzfs/Q8A6Sku5uVmWYIwcSOYK4hP4lss6FFV3SMiC4A/hzcsM8jrVe56bA+e7DQ+dnXk6i0F68KCKUzPSo3afohQryA3ltJ8t02WM3EjmCuB5/D1Qww+Pgx8IpxBmXM27qijsraFDW+9GHd6itPhvEqSfxnSLYd8y5A6Mat7NNtrmsmfkkbh1PSInK80P5un9p6gu28g7H0exoSb80NhzIjae/q594/7WTo7h7esiL7mm0Griz00nOnhUGP0dc5W1DZzydxpEUtcpQVuvEpUvhfGjJcliCh2/5+rOdnWw103LCIpKbq+mQ9VPrgM6cHoamY62dbNseausBXoC6TMP5LpYIMlCBP7Rk0QIuISkU9HKhhzzpGmDn72t1e4ZUVRREbgTMbs6ZnMmZ7JlkPRNWGuosY3QS6S7988TxYpLrF+CBMXRk0QqjoA3BShWMwQX35iHyku4c5rLnA6lKCUF3v4+6FT9EfRLOKK2mZSXUlcVDQlYudMcSWxMC/blh81cSGYJqYtIvIDEXmtiKwYvIU9sgT23IFGntnXwMeuLmHGlMh0rk5WeXEubT39vFzX6nQoZ1XUNHNR0RTSkiPbWVxiI5lMnAhmHsQq/897hmxT4OrQh2P6Brzc8/he5uZm8oHV85wOJ2ir/MuQbjnYFNE2/5H09nt5ua6V91w+N+LnLsvP5rGd9XT09Ed0pT9jQi2YYa5XjbWPCZ3/eqGG6pPt/Od7Vkb8m+9kTM9KZXHhFLYcauLjbyhxOhz2Hj9Db783YvMfhhosuXHwZDvLZudE/PzGhMqYCUJE8oGvAoWqeq2ILAKuUNWfhT26CNhYWceGzVXUt3RRmJPBujVl3Ly8yJE47n1yPydau0lLTqK9O7pLaAdSXuzhwS1H6OztJzPV2W/O22v8E+QcuJopK/AliAMn2ixBmJgWzF/xg8AvgM/7Hx/AV4k15hPExso61j+yiy7/usp1LV3c8fuXae7sZc3igojFsXnPCe794356+n0dvD39Xj736G5ExJFkNVHlxR5++txhXjrSzOtKnS29XlHbTFFOBgURmiA31OxpmaSnJFk/hIl5wSQIj6o+LCLrAVS1X0RiY6X6MWzYXHU2OQzq6fdy92N7ufuxvQ5F5dPVN8CGzVUxlSBeM28aqa4ktlY3OZ4gKmuaHWleAt/s8pIZbqvJZGJeMAmiQ0Ry8XVMIyKXA9EzVGUS6lu6Rnzu629ZErE47vj9roDbR4svGmWmJrN8To7j5b+Pt3ZR39rN/3Gws7w0383fDsb/CocmvgWTID4D/AFYKCJbgDzgbWGNKkIKczKoC/AhXJSTwTteMydicdz3bHXAOApzMiIWQ6isLvbwracPcLqjl+lZqY7EMDhBLhIVXEdSVpDN7yuO0dLZS06mM++DMZMVzDyIPcDr8A13/RdgMbA/nEFFyro1ZWQMK6iWkeJi3ZqyhIwjFMpLfMNdX3BwVnVFbTNpyUlcODNyE+SGGxzJdMBKbpgYFkyCeEFV+1V1j6ruVtU+fAsBxbyblxfxtVuWUJSTgeC7cvjaLUsi3u4fLXGEwsVFU3GnJTvazFRR28zFs6Y6uizr4Egm66g2sWzEJiYRKQCKgAwRWQ4MVoubAmRGILaIuHl5UVR8EEdLHJOV7ErisgW5jq0P0d03wO66Vj5QPt+R8w8qmJKOOy3ZSm6YmDZaH8Qa4H3ALOBbnEsQbcDnwhuWiWWri3N5Zl8DR093Mnt6ZL9L7KlvpW9AHRvBNEhEKC2wkhsmto2YIFT1l8AvReQtqvr7CMZkYtxqfz/Eluombr00cp39cK6DOhrKfZTmu3ly9/GoXEjJmGAE00g7S0SmiM8DIlIhIm8Oe2QmZi3My2aGO82RfoiK2mZmT88gz50W8XMPV5afTXNnH43tPU6HYsyEBJMgPqCqZ4A3AzOA9wP3jnWQiPxcRE6KyO4h2+4SkToR2eG/rR3h2GtEpEpEqkXkziB/FxMlRITVxR5eOHQKr1cjdl5VZXtNc1RcPYBvdTmAAydsJJOJTcEkiMFr47XAL1R155Bto3kQuCbA9u+o6jL/bdOrTibiAu4HrgUWAbf56z+ZGFJe7OFURy/7I9hJW9fSxcm2HkfnPww1ONTV+iFMrAomQWwXkafwJYjNIuIGxlwVRlWfA05PIKZLgWpVPayqvcBvsEWLYs7gMqSRHM1UURs9/Q8Anuw0crNSbSSTiVnBJIgPAncCr1HVTiAVXzPTRH1MRF72N0EF+ksuAo4OeXzMvy0gEfmwiGwTkW2NjVbaIFoUTE1nYV4WWw5FMEHUNJOR4uICf9NONCjNd3PgpCUIE5uCSRCrgWzgYhG5Et9M6onWMP4RsBBYBhzHN3x2uEDNVyM2ZKvqT1V1paquzMtztkCcOd/qYg8vHj5Nb39kliGtqG1m6eypJLucmyA3XFmBmwMn2lCNXF+MMaESzF/SuiG3fwceA+6ayMlUtUFVB1TVC/wnvuak4Y4Bs4c8ngXUT+R8xlmrij109Q1QWdsc9nN19w2wt/5M1DQvDSrNd9PROxCw1pYx0W7MBKGqNwy5vQm4CGiYyMlEZOaQh/8E7A6w20tAiYjMF5FU4FZ8xQJNjLl8QS5JEpl+iJePtdLv1ShMENkAVvrbxKSJXIsfw5ckRiUiD+Gr2VQmIsdE5IPAN0Rkl4i8DFwFfNq/b6GIbALfehPAx4DNwD7gYVXdM4E4jcOmZqRw8awctkSgcN/gCnLL50TXCm4lgyOZbKiriUHBLDn6fc71ASTh6z/YOdZxqnpbgM0BV6FT1Xp8o6QGH28CXjUE1sSe8uJcfvzXw7R19+FOTwnbeSpqm5nvySI32/kJckNNzUhh5tR0DtoVhIlBwVxBbAO2+28vAHeo6rvCGpWJG+XFHga8youHJzLiOTiqSmVtc9RdPQwqzbeaTCY2jXkF4a/JZMyErJgzjfSUJLYcauKNi/LDco6jp7toau+Nuv6HQaX52bxw+BQDXsWVZDWZTOwYrdz3LgIPLxVAVfXisEVl4kZ6iovXzJse1o7q7bW+q5PoTRBuevu91JzqYEFettPhGBO00a4gro9YFCaulRd7uPeP+zl5ppsZU9JD/voVNS1kpbrOLtITbQbjOtDQZgnCxJQR+yBUtUZVa/z7NAx5fJLgajEZA/gmzAFsDdNoporaZpbNyYna5pviGdmI2PKjJvYE00n9W86vvTTg32ZMUBbNnEJOZkpYyn939PSz73j0TZAbKjM1mTnTM62j2sScYBJEsr9oHgD++6nhC8nEm6QkYdVC3zKkoS45sfNYC17F8RXkxlIyw21F+0zMCSZBNIrIjYMPROQmwLkV6U1MKi/2cLy1m1eaOkL6upWDFVxnR3eCKCvI5pWmjojVpTImFIJJELcDnxORWhGpBe4APhzesEy8WR2m8t8VNc0szMtiamb4JuGFQmm+m36vhjxBGhNOwdRiOqSql+NbvGexqq5S1UPhD83EkznTMynKyQhpP4SqUlEbPSvIjWZwJJP1Q5hYEnQtJlVtV1X7320mZOgypAMhWob0laYOmjv7omYFudEs8GSTnCTWD2FiSvQUzjdxr7zEw5nufnbXtYbk9c6uIBcDCSI1OYl5niy7gjAxZcwEISKvqn4WaJsxY1m1MBcgZM1MFbXNuNOTKY6RyWdl+W4r+21iSjBXEC8Euc2YUXmy07igwM3WEC1DWlHTzLLZOSRF6QS54Urz3dSe7qSrd8DpUIwJyogJQkQKROQSIENElovICv/t9UBmxCI0cWV1sYeXjjTT3Te5D8m27j6qGtpiov9hUFlBNqpQfdJmVJvYMFotpjXA+/At+fktzpXXaAM+F96wTLwqL/bwwPOvsO1IM6tLPBN+nZ1HW1GN3gJ9gZTmnxvJtGTWVIejiV0bK+vYsLmK+pYuCnMyWLemjJuXFzkdVlwaMUH4y3z/UkTeoqq/j2BMJo5dOn86yUnC89VNk0oQ22uaEYFlUboGRCBzc7NITU6yfohJ2FhZx/pHXqarzzfhsK6li/WP7AKwJBEGwfRBzBKRKeLzgIhUiMibwx6ZiUtZacmsmDNt0v0QFbXNlM5wMyWMq9SFmitJKM7LpsqGuk7YvX/cfzY5DOrqG2DD5iqHIopvwSSID6jqGeDNwAzg/cC9YY3KxLVVxbnsqmulpbN37J0D8Hp9K8itmBs7Vw+DygrctvzoBD25+wQnznQHfK6+pSvC0SSGYBLEYN/DWuAXqroTK/dtJmF1sQdVeGGC5b8PN7Vzpruf5THU/zCoNN9NfWs3Z7r7nA4lZrR29fGZh3dw+6+2k+IK/NFTmJMR4agSQzAJYruIPIUvQWwWETfnl/8OSER+LiInRWT3kG0bRGS/iLwsIo+KSMCvgCJyRER2icgOEdkW7C9jYsPS2TlkpbrYMsFmpu01zUBsdVAPKivwzdmwq4jgbKlu4trvPsf/7qjnE1cXc+8tS8hIcZ23jytJWLemzKEI49uYa1IDHwSWAYdVtVNEcvE1M43lQeAHwH8N2fY0sF5V+0Xk68B6fMX/ArlKVa1qbBxKcSVx2YJctlRP7AqioqaFnMwUFniyQhxZ+JXM8I9kOtHOJXOnOxxN9OrqHeDrT+7nwa1HWODJ4ne3X3H2itGVlHR2FFN2WjJtPf0xMxcm1gSTIBRfob7rgXuALGDMdSNV9TkRmTds21NDHv4deGuwgZr4Ul7s4U/7T1LX0kXROJsHKmqbWR5DE+SGKsrJICvVZSOZRrHzaAuffngHhxs7eN+qedxxzQVkpJ67arh5edHZEUt9A17e8ZMX+Pwju1g+O4fZ022KVigF08T0Q+AK4Db/4zbg/hCc+wPAH0d4ToGnRGS7iIxaWlxEPiwi20RkW2NjYwjCMpEw0fLfrV19HDzZHpPNS+BbPKnESm4E1Dfg5dtPH+CWH22lq3eAX33wMu66cfF5yWG4FFcS37t1OQCf+E0lfQO23kYoBZMgLpjM360AABhESURBVFPVjwLdAKrazCRXlBORzwP9wK9H2KVcVVcA1wIfFZErR3otVf2pqq5U1ZV5eXmTCctEUGl+Np7stHEniMpaf/9DDM2gHs5qMr3awYY2bvnhVu579iA3LS3kyU9dGfQ8mdnTM/nqLUuorG3hu88cCHOkiSWYBNEnIi583+oRkTyC6KQeiYi8F19z1Tt1hPUnVbXe//Mk8Chw6UTPZ6KTiFBe7OuHGM8ypBW1LSSJr6M7VpUWuGlq76WpvcfpUBzn9SoP/O0w133/eY41d/Kjd67g2+9YxtSM8c1vuWFpIe9YOZsf/uUQW8Ow9nmiCiZB3IfvQ3qGiHwFeB746kROJiLX4OuUvlFVO0fYJ8s/UgoRycI3/2J3oH1NbCsv9tDU3sOBhuBrE1XWNlNWMIXstGC6z6JTab5vJFOiX0Uca+7knx/4O19+Yh9XlnjY/OkruXbJzAm/3pduXMR8Txaf+p8dnLLkGxLBrCj3a+DfgK8Bx4GbVfW3Yx0nIg/hq/paJiLHROSD+EY1uYGn/UNYf+zft1BENvkPzQeeF5GdwD+AJ1T1yQn8bibKlfv7IYIt/z3gVSprW1gRQ+U1Ainz12RK1MWDVJXfbjvKNd/9G7uOtfKNt1zMf75nJTPcY459GVVmajLfv205LZ19rPvdy+O6MjWBjfk1TESmAyeBh4ZsS1HVUWf6qOptATb/bIR96/HNs0BVDwNLx4rLxL6inAzme7LYUt3EB1fPH3P/gyfbaO/pj6kKroHkudPIyUzhQAJWdW1q72H9I7t4em8Dl86fzrfetjSkI48WF07lc2sv4K7H9vLg1iO8v3zs/1dmZMFcp1cAs4FmfDOoc4DjInIS+JCqbg9jfCbOlRfn8mhFHX0DXlJco1/QVtT4V5CL0RFMg0SE0nx3wl1BPLn7BJ9/dBdt3f18fu2FfHD1/LAMVX7vqnn87WATX9u0n0vnT2dxoVXOnahg+iCeBNaqqkdVc/GNLHoY+Ai+IbDGTNjqYg8dvQPsPNoy5r4Vtc1Mz0plbm7sj3Uvy3dT1dCWEM0gZ7r7+OzDO7n9V9spmJrO459YzYeuXBC2eSwiwoa3LSUnM4WPP1RJZ29/WM6TCIJJECtVdfPgA/9ktytV9e+ALT1qJuXyBbmIBNcPUVHTzIo5OYjE3gS54Urzs2nr7h+x+Fy82FrdxDXfeY6NO+r4xNXFPPqR8rPrYoTT9KxUvvuOZbzS1MHdf9gb9vPFq2ASxGkRuUNE5vpv/wY0+4e+2qwUMyk5maksKZrK1jHKbjR39HK4qSOm5z8MdXbxoDhtZuruG+Dux/bwzw+8SFqKi9/dfgWfeXMZqcnBfOSExqpiDx95/UL+Z9tRHttZH7HzxpNg/rX+Gd+qchv9t9n+bS7g7eELzSSK8mIPFbXNdPSM3BRQeTR2C/QFMpggDo5jiG+s2Hm0hevu+xu/2HKE914xl02feK1jlXc/9cZSls/J4XOP7OLo6YAj680oRk0Q/quE76rqx1V1uf/2cVVtVNVeVa2OUJwmjpUv9NDvVf7xyukR99le04wrSbg4TpbqnJaVygx3GlVxNBeib8DLd/ylMjr9pTLuvumiUUtlhFuKK4n7rBTHhI2aIFR1AMgTkUmV1jBmNCvnTSM1OWnUfoiKmhYunOkmMzV2J8gNV1YQPyU3qk/6SmV879mD3DjOUhnhZqU4Ji6Yv7YjwBYR+QPQMbhRVb8drqBMYklPcfGaedNGrMvUP+Bl57EW3nbJrAhHFl4lM9z89z9q8Ho1JivTgq9UxoNbj/D1J/eTmerih+9cwdpJzIYOlxuWFvL8wSZ++JdDlC/0sKo4OpJXtAumD6IeeNy/r3vIzZiQWbXQw/4TbTS2vbpEQlVDG529A3HTQT2orCCb7j4vR5tjs228rqWLdz7wIvc8vpfVxb5SGdGYHAYNLcVxumNiy90mmjGvIFT17kgEYhLb6mIPGzZXsfVQEzctKzrvuYoYXkFuNIMd1Qca2pmbG/2LH22srDu7UM/UzBS6evpJdiXx9bcs4e0rZ0f98OPBUhz/dP9W1v12Jw+8d2XUx+y0Ma8gRCTPv1ToJhH50+AtEsGZxHFR0VSmpCcHHO5aUdtCnjuNWdPia93hkrMJIvr7ITZW1rH+kV3UtXShQEtnH31e5TNvLuUdr5kTMx+0iwunsn7tBTy7/yS/3HrE6XCiXjBNTL8G9gPzgbvx9Um8FMaYTAJyJQmrFnp4vrrpVbOLK2rjZ4LcUNlpycyalhETcyE2bK6iq2/gvG1ehZ8/f8SZgCbhfavmcfUFM/jqpv3srT/jdDhRLZgEkauqPwP6VPWvqvoB4PIwx2USUHlxLnUtXdScOtcm39TeQ82pzrhrXhpUGiOLB9W3dI1rezQTETa89WJ/KY4KK8UxiqAWDPL/PC4i14nIcnwT54wJqcHy31sOnRvNdLb/Ic46qAeV5rs51Nge9ePz89yBq+oUjnM98WiRm53Gd9+xjMNWimNUwSSIL4vIVOCzwL8CDwCfCmtUJiHN92RRODX9vOGuFbUtpLiEJUXxMUFuuLKCbPoGlJpTHWPv7BBVDbhAU0aKi3VryhyIKDSsFMfYglkw6HFVbVXV3ap6lapeAiyMQGwmwYgIq4o9bD10Cq/X1w9RUdvMosKppKc4Nxs3nM7VZIrekhtP7DrO4aYO3nrJLIpyMhB8a3l87ZYl3Ly8aMzjo5mV4hjdRCtnfSakURjjt7rYQ0tnH3uPn6FvwMvLx2J/BbnRLMzLJkmI2pIbXb0DfPWJfVw4cwpff8vFbLnzal659zq23Hl1zCcHsFIcY5logoiv4SQmaqwqzgV85b/3HT9Dd583bjuowTeLfF5uVtQuHvTjvx6ivrWbu25YhCtGZ3uPZWgpju89c9DpcKLKRAvbxP8qJ8YRM9zplOZns6W6iXR/aehYX2J0LNE6kqmupYsf//UQ1108k8sW5DodTljdsLSQvx1s5P6/VLOqOJdVC60UB4xyBSEibSJyJsCtDSiMYIwmwZQXe3jpyGleOHyKginpMTtSJlilBW6OnOqge9g8A6d9ddM+ROBzay90OpSIuOvGxcz3ZPFpK8Vx1ogJQlXdqjolwM2tqmNeeYjIz0XkpIjsHrJtuog8LSIH/T8DfjUUkWtEpEpEqkXkzon9aiZWuZKE7j4vm/c00NLVy8bKOqdDCquyfDdehUON0dNR/ffDp3ji5ePc/rqFFMV5gh40WIqjuaOPf/vdzoRYDnYs4Vze6UHgmmHb7gSeVdUS4Fn/4/P416C4H9/a14uA20RkURjjNFFkY2Udv/p7zdnH3X1e1j+yK66TRGl+NhA9JTcGvMrdj+2lKCeDf7kysQYsDpbieGafleKAMCYIVX0OGL4CzE3AL/33fwncHODQS4FqVT2sqr3Ab/zHmQSwYXMV3X3njyTp6htgw+YqhyIKv3meLFJcEjVDXR/6Ry37jp/hc2svdHSxH6dYKY5zIrdArE++qh4H8P+cEWCfIuDokMfH/NsCEpEPi8g2EdnW2NgY0mBN5MVTSYdgpbiSWJiXzcEouIJo7ezjW09Vcdn86axdUuB0OI6wUhznRDpBBCPQWLoRGwNV9aequlJVV+bl5YUxLBMJI3VIx31Hdb47KuZCfOeZA7R29XHXjYvjrjjieORmp/EdfymOex5L3FIckU4QDSIyE8D/82SAfY4Bs4c8noVv0SKTANatKSNj2KzpWC/pEIyyAjfHmrto73Hu22rViTb+399r+OfL5nDhzCmOxREtyos9/N/XLeQ3Lx3l8ZcT8yMo0gniD8B7/fffC/xvgH1eAkpEZL5/Lexb/ceZBHDz8iK+dsuSuCvpMJaSGb6OaqeamVSVex7fQ3ZaMp99U3wn4/H49JtKWTY7h/UJWoojbCvAi8hDwOsBj4gcA74E3As8LCIfBGqBt/n3LQQeUNW1qtovIh8DNgMu4OequidccZroc/PyorhPCMOVFZxbPGi5AzPHN+9pYEv1Ke66YRHTslIjfv5oleJK4vu3LWft9/7GJ39TycP/cgXJrmhsmQ+PsCUIVb1thKfeEGDfemDtkMebgE1hCs2YqDN7WibpKUkcaIj8SKbuvgG+smkvpfnZvOvyuRE/f7QbLMXx8Ycq+e4zB/nXOG/uHCpxUqExUSwpSRwrufGz51/h6OkuvnTD4oT6djweNywt5O0rZ3H/X6rZOmS9knhn/xuMiRKl+e6ILz96orWb+/9czZrF+WcXbDKBJWIpDksQxkSJ0vxsTrb10BzBD597/7iPfq/yheusWMFYMlOTue/WxCrFEbY+CGPM+AwuHnSgoS0i1VO315xm4456PnZVMbOnZ4b9fPHgoqKp3HntBdzz+F6W3fM0Z7r6KMzJYN2aMkcGVmysrGPD5irqW7rCEoddQRgTJc6OZDoZ/o5qr1e56w97KZiSzkeuSqx6S5M1LTOFJIHWrj4UX1l0J+qFbaysY/0ju6hr6QpbHHYFYUyUKJiSjjs9OSKLB/12+1F21bXyvVuXkZlqHwPj8c2nDuAd1rrU1TfAHb9/mUcimCRePHyKnv7AdctCdRVh/zOMiRIiQlkESm6c6e5jw+YqVs6dxo1LbWmX8RqpLlhPv5czXX0Ri2N4chgUyrplliCMiSIl+W7+uPs4qhq2Wkj3PXOQUx29PPj+SxO63tJEFeZkUBfgQ7goJ4ONHy2PWBzl9/4pYByhrFtmfRDGRJGy/GxaOvtobOsJy+tXn2znwa1HeMfK2VxUNDUs54h30VIvLBJx2BWEMVGk9GzJjXZmTEkP6WurKv/x+F4yUlwJNRs41Abb98M5eiha4rAEYUwUKfMPda1qaGN1SWgnrv1p/0n+eqCRL1x3IZ7stJC+dqKJlnph4Y7DmpiMiSK52Wl4slNDPpKpp3+A/3h8LwvysnjPFfNC+tomflmCMCbKlMwI/UimB7cc4cipTr54/SJSk+3P3gTH/qcYE2XKCtwcbGgLWSmHk23dfP9P1bzhghm8vizQKr/GBGYJwpgoU5rvpqN3IOAQxon4xpNV9PQP8IXrrd6SGR9LEMZEmbIC3+pyoSj9veNoC7/bfowPrJ7PfE/WpF/PJBZLEMZEmZLBkUwnJleTyVdvaQ957jQ+fnVJKEIzCcYShDFRZkp6CjOnpk/6CuLRyjp2HG3hjmsuIDvNRrSb8bMEYUwUmuziQe09/dz75H6Wzs7hligYr29ikyUIY6JQWYGb6sZ2BoaXDQ3SD/5UTWNbD3fdsIikJKu3ZCYm4glCRMpEZMeQ2xkR+dSwfV4vIq1D9vlipOM0xkml+W56+73UnOoY97FHmjr4+fOv8JYVs1g+Z1oYojOJIuINk6paBSwDEBEXUAc8GmDXv6nq9ZGMzZhoUTZkdbkFednjOvbLT+wlxSXccY3VWzKT43QT0xuAQ6pa43AcxkSV4hnZiIx/JNNfDzTyzL6TfOzqkpAX+zOJx+kEcSvw0AjPXSEiO0XkjyKyeKQXEJEPi8g2EdnW2NgYniiNibCMVBdzpmeOayRT34CXex7bw7zcTD6wel74gjMJw7EEISKpwI3AbwM8XQHMVdWlwPeBjSO9jqr+VFVXqurKvLy88ARrjANK893jShD/9UINhxo7+MJ1i0hLdo19gDFjcPIK4lqgQlUbhj+hqmdUtd1/fxOQIiKhrX1sTJQry3fzSlMHPf0DY+57qr2H7z5zgCtL83jDhVZvyYSGkwniNkZoXhKRAvGvhSgil+KL81QEYzPGcSX52fR7lVeaxh7J9M2nqujqHeCL1y+yZURNyDiSIEQkE3gT8MiQbbeLyO3+h28FdovITuA+4FYNVWlLY2JEWcFgyY3Rm5l217Xym5eO8t5V8yieMb4RT8aMxpH596raCeQO2/bjIfd/APwg0nEZE00WeLJJTpJR+yFUffWWpmem8ok3WL0lE1pOj2IyxowgNTmJ+Z4sDjSMPNT1Dzvr2VbTzLo1ZUzNSIlgdCYRWIIwJoqVFow8kqmzt5+vbdrPRUVTeNvK2RGOzCQCSxDGRLHSGW5qT3fS2dv/qud+9JdDnDjTzV03LMZl9ZZMGFiCMCaKlRVkowrVJ89vZjp6upOfPHeYm5YVsnLedIeiM/HOEoQxUaw0P/BIpq88sQ+XCHdee4ETYZkEYQnCmCg2NzeL1OQkDg65gtha3cSTe07wkdcvZObUDAejM/HOEoQxUcyVJJTMyD57BdE/4OXux/Yya1oGH7pygcPRmXhnCcKYKDe0JtN//6OWqoY2vnDdhaSnWL0lE16WIIyJcqX5bo63dlN7qpNvPXWAVQtzWbO4wOmwTAKwlcyNiXKnO3oAuHLDnwF4bYnH6i2ZiLArCGOi2MbKOv7rhfPX07rv2Wo2VtY5FJFJJJYgjIliGzZX0dPvPW9bV98AGzZXORSRSSSWIIyJYvUtXePabkwoWYIwJooV5gSe5zDSdmNCyRKEMVFs3ZoyMoYNZ81IcbFuTZlDEZlEYqOYjIliNy8vAnx9EfUtXRTmZLBuTdnZ7caEkyUIY6LczcuLLCEYR1gTkzHGmIAsQRhjjAnIEoQxxpiALEEYY4wJyBKEMcaYgERVnY4hZESkEagZc8fAPEBTCMOJZfZenM/ej/PZ+3FOPLwXc1U1L9ATcZUgJkNEtqnqSqfjiAb2XpzP3o/z2ftxTry/F9bEZIwxJiBLEMYYYwKyBHHOT50OIIrYe3E+ez/OZ+/HOXH9XlgfhDHGmIDsCsIYY0xAliCMMcYElPAJQkSuEZEqEakWkTudjsdJIjJbRP4sIvtEZI+IfNLpmJwmIi4RqRSRx52OxWkikiMivxOR/f7/I1c4HZOTROTT/r+T3SLykIikOx1TqCV0ghARF3A/cC2wCLhNRBY5G5Wj+oHPquqFwOXARxP8/QD4JLDP6SCixPeAJ1X1AmApCfy+iEgR8AlgpapeBLiAW52NKvQSOkEAlwLVqnpYVXuB3wA3ORyTY1T1uKpW+O+34fsASNiFCERkFnAd8IDTsThNRKYAVwI/A1DVXlVtcTYqxyUDGSKSDGQC9Q7HE3KJniCKgKNDHh8jgT8QhxKRecBy4EVnI3HUd4F/A7xOBxIFFgCNwC/8TW4PiEiW00E5RVXrgG8CtcBxoFVVn3I2qtBL9AQhAbYl/LhfEckGfg98SlXPOB2PE0TkeuCkqm53OpYokQysAH6kqsuBDiBh++xEZBq+1ob5QCGQJSLvcjaq0Ev0BHEMmD3k8Szi8DJxPEQkBV9y+LWqPuJ0PA4qB24UkSP4mh6vFpFfORuSo44Bx1R18Iryd/gSRqJ6I/CKqjaqah/wCLDK4ZhCLtETxEtAiYjMF5FUfJ1Mf3A4JseIiOBrY96nqt92Oh4nqep6VZ2lqvPw/b/4k6rG3TfEYKnqCeCoiJT5N70B2OtgSE6rBS4XkUz/380biMNO+2SnA3CSqvaLyMeAzfhGIfxcVfc4HJaTyoF3A7tEZId/2+dUdZODMZno8XHg1/4vU4eB9zscj2NU9UUR+R1QgW/0XyVxWHbDSm0YY4wJKNGbmIwxxozAEoQxxpiALEEYY4wJyBKEMcaYgCxBGGOMCcgShEl4IpIrIjv8txMiUue/3y4iP4xQDK8frBgrIjcmemVhEx0Seh6EMQCqegpYBiAidwHtqvpNB+P5Awk8YdNED7uCMGYEw77V3yUivxSRp0TkiIjcIiLfEJFdIvKkv0QJInKJiPxVRLaLyGYRmRngdd/mX0Ngp4g8F+D594nID/z380XkUf++O0VklX/7u0TkH/4rnZ/4S9cbE1KWIIwJ3kJ85b9vAn4F/FlVlwBdwHX+JPF94K2qegnwc+ArAV7ni8AaVV0K3DjGOe8D/urfdwWwR0QuBN4BlKvqMmAAeOekfztjhrEmJmOC90dV7RORXfhKszzp374LmAeUARcBT/vK8+DCVwp6uC3AgyLyML4ib6O5GngPgKoOAK0i8m7gEuAl/3kygJMT/7WMCcwShDHB6wFQVa+I9Om5OjVefH9LAuxR1VGX4lTV20XkMnxXIztEZNk44xDgl6q6fpzHGTMu1sRkTOhUAXmDazWLSIqILB6+k4gsVNUXVfWLQBPnl5wf7lng//qPc/lXdnsWeKuIzPBvny4ic0P8uxhjCcKYUPEvW/tW4OsishPYQeA1Ajb4O7d3A88BO0d52U8CV/mbtbYDi1V1L/AF4CkReRl4GnhVZ7gxk2XVXI0xxgRkVxDGGGMCsgRhjDEmIEsQxhhjArIEYYwxJiBLEMYYYwKyBGGMMSYgSxDGGGMC+v/9lbHHgF8dHwAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.plot(range(len(q6clusters)), [max(x) for x in q6clusters], 'o-')\n",
    "plt.xlabel(\"Time slice\")\n",
    "plt.ylabel(\"Largest cluster size\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The final thing to do is to remove the split files after use."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 44,
   "metadata": {},
   "outputs": [],
   "source": [
    "for file in files:\n",
    "    os.remove(file)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Using [ASE](https://wiki.fysik.dtu.dk/ase/)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The above example can also done using ASE. The ASE read method needs to be imported."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 36,
   "metadata": {},
   "outputs": [],
   "source": [
    "from ase.io import read"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 37,
   "metadata": {},
   "outputs": [],
   "source": [
    "traj = read(\"traj.light\", format=\"lammps-dump-text\", index=\":\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In the above function, `index=\":\"` tells ase to read the complete trajectory. The individual slices can now be accessed by indexing."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 38,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Atoms(symbols='H500', pbc=False, cell=[18.21922, 18.22509, 18.36899], momenta=...)"
      ]
     },
     "execution_count": 38,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "traj[0]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can use the same functions as above, but by specifying a different file format."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 46,
   "metadata": {},
   "outputs": [],
   "source": [
    "q6clusters_ase = [calculate_q6_cluster(x, format=\"ase\") for x in traj]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We will plot and compare with the results from before,"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 47,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Text(0, 0.5, 'Largest cluster size')"
      ]
     },
     "execution_count": 47,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYgAAAEGCAYAAAB/+QKOAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAgAElEQVR4nO3dd3hU95Xw8e/RqEsDAo2QkOio2GBMMXFBxImdgo3rOs3e9ORN1m96WdYmySa2N8UJqU6ctk7i7Juss05is7FNjEuKY3AcgwSmCgRGAgmEBJJQb3PeP2YEQh5JI2lm7pTzeZ55NHPn3rlHA5oz91fOT1QVY4wxZrgkpwMwxhgTnSxBGGOMCcgShDHGmIAsQRhjjAnIEoQxxpiAkp0OIJQ8Ho/OmzfP6TCMMSZmbN++vUlV8wI9F1cJYt68eWzbts3pMIwxJmaISM1Iz1kTkzHGmIAsQRhjjAnIEoQxxpiALEEYY4wJyBKEMcaYgOJqFJMxobSxso4Nm6uob+miMCeDdWvKuHl5kdNhGRMxliCMCWBjZR3rH9lFV98AAHUtXax/ZBeAJQmTMKyJyZgANmyuOpscBnX1DbBhc5VDERkTeZYgjAmgvqVrXNuNiUeWIIwJoDAnY1zbjYlHliCMCWDdmjLSks//88hIcbFuTZlDERkTeWFLECIyW0T+LCL7RGSPiHzSv32DiOwXkZdF5FERyRnh+CMisktEdoiIFVgyEXXz8iLWXlRw9nFORgpfu2WJdVCbhBLOK4h+4LOqeiFwOfBREVkEPA1cpKoXAweA9aO8xlWqukxVV4YxTmMC6vMqBVPSSU9J4i2XzLLkYBJO2BKEqh5X1Qr//TZgH1Ckqk+par9/t78Ds8IVgzGTUVnbwiXzplEyw82BhjanwzEm4iLSByEi84DlwIvDnvoA8McRDlPgKRHZLiIfHuW1Pywi20RkW2NjYyjCNYYTrd3UtXRxyZxplOZbgjCJKewJQkSygd8Dn1LVM0O2fx5fM9SvRzi0XFVXANfia566MtBOqvpTVV2pqivz8gKueWHMuFXUNgOwYu40ygqyaTjTQ0tnr8NRGRNZYU0QIpKCLzn8WlUfGbL9vcD1wDtVVQMdq6r1/p8ngUeBS8MZqzFDVdQ0k5acxKKZUyjJdwNwoKHd4aiMiaxwjmIS4GfAPlX99pDt1wB3ADeqaucIx2aJiHvwPvBmYHe4YjVmuIraZpYUTSU1OYkyf4KosmYmk2DCeQVRDrwbuNo/VHWHiKwFfgC4gaf9234MICKFIrLJf2w+8LyI7AT+ATyhqk+GMVZjzurpH2B33RkumTsNgJlT03GnJXPQEoRJMGEr1qeqzwMS4KlNAbYNNimt9d8/DCwNV2zGjGZ33Rl6B7wsn+NLECJCaYGbqhOWIExisZnUxgxTebaD+twczsGRTCN0mRkTlyxBGDPM9ppmZk/PYIY7/ey20vxsmjv7aGzvcTAyYyLLEoQxQ6gqFbXNrPA3Lw0a7Kg+cMJGMpnEYQnCmCHqW7tpONPzqgRRWjA41NX6IUzisARhzBAVNf7+h2EJwpOdRm5WqiUIk1AsQRgzxPaaZjJSXFww0/2q50rys20uhEkoliCMGaKytpmLZ00lxfXqP42yfDcHG9ptJJNJGJYgjPHr7htgT/0ZVsydFvD50gI37T391Ld2RzgyY5wRtolyJvZsrKxjw+Yq6lu6KMzJYN2asoRaA+HlY630e/VV/Q+Dzo1kaqPIlh41CcCuIAzgSw7rH9lFXUsXCtS1dLH+kV1srKxzOrSIOVvBdU7ARQ7PFu2zfgiTKCxBGAA2bK6iq2/gvG1dfQNs2FzlUESRV1HTzLzcTHKz0wI+PzUjhYIp6RywkhsmQViCMADUt3SNa3u88U2QaxmxeWlQaYHbriBMwrAEYQAoHKFNfaTt8ebo6S6a2ntG7KAeVJafTfXJdga8NpLJxD9LEAaAdWvKcCWdX3w3I8XFujVlDkUUWef6H8a4gsh309PvpfZ0wKVMjIkrliAMAAvzshnwKtlprrPb7rlpccKMYqqobSYr1UVZwasnyA01+LyV/jaJwBKEQVW5+7E9eLJT2br+Dfzifa8BEqd5CXwJYunsnFddRQ1XPCMbsJpMJjFYgjD8YWc922qaWbemjCnpKVw6fzrJScLz1U1OhxYRnb397DvednYFudFkpiYzZ3qmdVSbhGAJIsF19vbztU37uahoCm+9ZDYAWWnJLJ+Tw5YESRA7j7YyMMoEueFK8922/KhJCJYgEtwP/3yIE2e6ueuGxec1r5QXe9hV10prZ5+D0UXGYAf18hEmyA1XVpDN4cYOevu94QzLGMdZgkhgtac6+enfDnPTskJWzpt+3nPlxR5U4YXD8X8VUVHTzIK8LHIyU4PavzTfTb9XeaWpI8yRGeOssCUIEZktIn8WkX0iskdEPunfPl1EnhaRg/6fAa/rReQaEakSkWoRuTNccSayr2zai0uEO6+94FXPLZudQ1aqK+77IVSVyqMtXBJk8xL4EgRYyQ0T/8J5BdEPfFZVLwQuBz4qIouAO4FnVbUEeNb/+Dwi4gLuB64FFgG3+Y81IbKluonNexr46FULmTn11aOVUlxJXLYgl63VpxyILnKOnOrkdEfvmBPkhlqQl4UrSazkhol7YUsQqnpcVSv899uAfUARcBPwS/9uvwRuDnD4pUC1qh5W1V7gN/7jTAj0D3i5+7E9zJ6ewf957YIR91u1MJfDTR3UxXG5jZFWkBtNWrKL+Z4sG+pq4l5E+iBEZB6wHHgRyFfV4+BLIsCMAIcUAUeHPD7m3xbotT8sIttEZFtjY2Mow45bv36xlgMN7Xx+7SLSU1wj7re6xAMQ16OZttc2405LpsQ/vyFYZfluSxAm7oU9QYhINvB74FOqeibYwwJsC1j8RlV/qqorVXVlXl7eRMNMGM0dvXz76QOUF+eyZnH+qPuW5bvxZKfGdYKoqGlm2ZwcksaYIDdcab6bmtOddPUOjL2zMTEqrAlCRFLwJYdfq+oj/s0NIjLT//xM4GSAQ48Bs4c8ngXUhzPWRPGtp6to7+nnSzcsRmT0D0URobzYw5bqU3G5zGZ7Tz8HGtrG1bw0qDQ/G1WoPtkehsiMiQ7hHMUkwM+Afar67SFP/QF4r//+e4H/DXD4S0CJiMwXkVTgVv9xZhL21p/hv1+s5d2Xzz07Emcs5Qs9NLX3cKAh/j4Idx5twauMq4N6UGmBjWQy8S+cVxDlwLuBq0Vkh/+2FrgXeJOIHATe5H+MiBSKyCYAVe0HPgZsxte5/bCq7gljrHFvsN7S1IwUPv3G0qCPK/f3Q8TjcNftNc2I+Ib0jtfc6ZmkJifZjGoT18K2JrWqPk/gvgSANwTYvx5YO+TxJmBTeKJLPJt2neDFV07z5ZsvYmpmStDHFeVkMN+TxdbqJj64en4YI4y8itpmSmZkMzUj+PdjULIrieK8bLuCMHHNZlIngK7eAb66aR8XFLi57dI54z5+1cJc/n74FH0D8VNawutVKoNYQW40ZQVumwth4poliATwk+cOUdfSxV03Lh6znHUgq4s9dPQOsPNoSxiic8bhpnZau/omlSBK8rOpb+3mTHf816syickSRJyra+nix389xHVLZnL5gtwJvcYVC3MRgS1xNKu6osaX7CbSQT2ozN/Rb/0QJl5ZgohzX920D1VYv/bV9ZaClZOZykWFU+NqPkRFbTNTM1JY4Mma8GsMjgSLxxFexoAliLj24uFTPPHycW5/3UJmTcuc1GuVF3uoqG2mo6c/RNE5q6K2meUTmCA3VFFOBlmpLlt+1MQtSxBxasCr3PXYXgqnpnP76xZO+vVWF3vo9yr/OHI6BNE5q7WrjwMN7eOq4BpIUpJQYiU3TByzBBGnfvNSLfuOn+Fz111IRurI9ZaCtXLeNFKTk9hyMPabmXYcnXz/w6DS/GxLECZuWYKIQ62dfXxzcxWXzZ/OdUtmhuQ101NcrJw7LS4mzFXUNJMksHQCE+SGK81309Tey6n2nhBEZkx0sQQRh77zzAFau/qCqrc0HuXFHvafaKMpxj8MK2qbKc13k502+XmiZQXWUW3ilyWIOHOgoY3/9/cabrt0DosKp4T0tcuLfWU3th6K3eGuXq+yo7aFS0LQvATnhrpaM5OJR5Yg4oiqcs9je8lKdfHZN5eF/PWXFE3FnZ4c0/0QB0+209bTP6kJckPludPIyUyxkhsmLgWVIERktYi8338/T0TiqyhPnHhqbwPPVzfxmTeVMj0rNeSv70oSVi3M5fnqppgt/11R619BLkRXECJC6QwruWHi05gJQkS+BNwBrPdvSgF+Fc6gzPh19w3w5Sf2Upqfzbsunxu286wu9lDX0kXt6c6wnSOcttc0Mz0rlXm5k5sXMlRpga9oX6wmTWNGEswVxD8BNwIdcLbqanCLCZiI+dnzr3D0dBdfvH4xya7wtRyuKo7t8t8Vtc2smJMT0s77snw3bd39NJyJ7c57Y4YL5pOkV31fjRRARCZem8CExYnWbu7/czVvXpR/dh3pcFngyWLm1HS2xmBdppbOXg43drA8RP0PgwZLblg/hIk3wSSIh0XkJ0COiHwIeAZ4ILxhmfH4+pP76fcqX7huUdjPdXYZ0kNNeL2x1aRSWeufIBemBGH9ECbejJkgVPWbwO/wrS1dBnxRVe8Ld2AmONtrmnm0so4PvXY+c0LYrj6a8uJcWjr72Hv8TETOFyrba5pxJQlLZ08N6etOy0olz51mVxAm7gTTSf3vwH5VXaeq/6qqT4vIhyMQmxmD1+tbRjR/ShofeX1xxM5bvjA2+yEqapu5cKabzNTQL6RYZjWZTBwKponp48BmEblqyLbbwxSPGYffbT/Gy8daWX/thWSFYFZwsGZMSac0Pzumyn8PeJWdRye3gtxoSvPdHGxoj7lmN2NGE0yCqAOuAe4VkXX+baEbAmIm5Ex3H9/YvJ9L5k7jpmWFET//qoUeXjpymu6+gYifeyKqTrTR0TsQtgRRVpBNV98Ax5q7wvL6xjghqPGQqloLvA5YJCK/BTLGOkZEfi4iJ0Vk95Bt/yMiO/y3IyKyY4Rjj4jILv9+24L8XRLK9589yKmOXu4Kcb2lYK0u9tDd5z078SzabffHGaoSG8OV2EgmE4eCSRDbAFS1W1XfD/wFCGaa7oP4rjzOUtV3qOoyVV2Gr9P7kVGOv8q/78ogzpVQDjW284stR3j7JbNZMiu0Ha7BumzBdFxJEjPDXStrmvFkpzFr2pjfbSakZEY2YDWZTHwJZhTTh4Y9vl9VFwRx3HNAwNVlxPeV9+3AQ0HGaYb4j8f3kpHiYt01oa+3FCx3egpLZ02NmY7qcEyQG8qdnkJRToYlCBNXRkwQIvKw/+cuEXl5+G2S530t0KCqB0d4XoGnRGT7WCOmROTDIrJNRLY1NjZOMqzo96f9DfylqpFPvrEET3aao7GsLvbw8rEWWrv6HI1jLE3tPRw51Rmy+ksjKStw2/KjJq6MdgXxSf/P64EbAtwm4zZGv3ooV9UVwLXAR0XkypF2VNWfqupKVV2Zl5c3ybCiW2+/l/94fB8L8rJ4zxXznA6H8mIPXvWtfR3NBifIhav/YVBpvpvDjR30DXjDeh5jImXEBKGqx/13m4CjqloDpAFLgfqJnlBEkoFbgP8Z5dz1/p8ngUeBSyd6vnjyiy2v8EpTB/9+/SJSk52v1L58zjQyUlxRP9y1oraZ5CRhSVF4+2tK87PpHfBSc6ojrOcxJlKC+ZR5DkgXkSLgWeD9+DqgJ+qN+CbeHQv0pIhkiYh78D7wZmB3oH0Tycm2br7/p2quvmAGV5XNcDocAFKTk7h0/vSo74eoqGlmceEU0lMmvzb3aM7WZDphq8uZ+BBMghBV7cT3rf/7qvpPwJhFf0TkIeAFoExEjonIB/1P3cqw5iURKRSRTf6H+cDzIrIT+AfwhKo+GdyvE782PFlFT/8A/359+OstjcfqYg+HGjs40drtdCgB9Q142XmsJeQF+gIpnpFNkthIJhM/gpl+KyJyBfBOYPBDfszjVPW2Eba/L8C2emCt//5hfM1Yxm/n0RZ+u/0Y/3LlAuZ7oquY7qriXAC2VDfxlktmORzNq+0/3kZ3nzfs/Q8A6Sku5uVmWYIwcSOYK4hP4lss6FFV3SMiC4A/hzcsM8jrVe56bA+e7DQ+dnXk6i0F68KCKUzPSo3afohQryA3ltJ8t02WM3EjmCuB5/D1Qww+Pgx8IpxBmXM27qijsraFDW+9GHd6itPhvEqSfxnSLYd8y5A6Mat7NNtrmsmfkkbh1PSInK80P5un9p6gu28g7H0exoSb80NhzIjae/q594/7WTo7h7esiL7mm0Griz00nOnhUGP0dc5W1DZzydxpEUtcpQVuvEpUvhfGjJcliCh2/5+rOdnWw103LCIpKbq+mQ9VPrgM6cHoamY62dbNseausBXoC6TMP5LpYIMlCBP7Rk0QIuISkU9HKhhzzpGmDn72t1e4ZUVRREbgTMbs6ZnMmZ7JlkPRNWGuosY3QS6S7988TxYpLrF+CBMXRk0QqjoA3BShWMwQX35iHyku4c5rLnA6lKCUF3v4+6FT9EfRLOKK2mZSXUlcVDQlYudMcSWxMC/blh81cSGYJqYtIvIDEXmtiKwYvIU9sgT23IFGntnXwMeuLmHGlMh0rk5WeXEubT39vFzX6nQoZ1XUNHNR0RTSkiPbWVxiI5lMnAhmHsQq/897hmxT4OrQh2P6Brzc8/he5uZm8oHV85wOJ2ir/MuQbjnYFNE2/5H09nt5ua6V91w+N+LnLsvP5rGd9XT09Ed0pT9jQi2YYa5XjbWPCZ3/eqGG6pPt/Od7Vkb8m+9kTM9KZXHhFLYcauLjbyhxOhz2Hj9Db783YvMfhhosuXHwZDvLZudE/PzGhMqYCUJE8oGvAoWqeq2ILAKuUNWfhT26CNhYWceGzVXUt3RRmJPBujVl3Ly8yJE47n1yPydau0lLTqK9O7pLaAdSXuzhwS1H6OztJzPV2W/O22v8E+QcuJopK/AliAMn2ixBmJgWzF/xg8AvgM/7Hx/AV4k15hPExso61j+yiy7/usp1LV3c8fuXae7sZc3igojFsXnPCe794356+n0dvD39Xj736G5ExJFkNVHlxR5++txhXjrSzOtKnS29XlHbTFFOBgURmiA31OxpmaSnJFk/hIl5wSQIj6o+LCLrAVS1X0RiY6X6MWzYXHU2OQzq6fdy92N7ufuxvQ5F5dPVN8CGzVUxlSBeM28aqa4ktlY3OZ4gKmuaHWleAt/s8pIZbqvJZGJeMAmiQ0Ry8XVMIyKXA9EzVGUS6lu6Rnzu629ZErE47vj9roDbR4svGmWmJrN8To7j5b+Pt3ZR39rN/3Gws7w0383fDsb/CocmvgWTID4D/AFYKCJbgDzgbWGNKkIKczKoC/AhXJSTwTteMydicdz3bHXAOApzMiIWQ6isLvbwracPcLqjl+lZqY7EMDhBLhIVXEdSVpDN7yuO0dLZS06mM++DMZMVzDyIPcDr8A13/RdgMbA/nEFFyro1ZWQMK6iWkeJi3ZqyhIwjFMpLfMNdX3BwVnVFbTNpyUlcODNyE+SGGxzJdMBKbpgYFkyCeEFV+1V1j6ruVtU+fAsBxbyblxfxtVuWUJSTgeC7cvjaLUsi3u4fLXGEwsVFU3GnJTvazFRR28zFs6Y6uizr4Egm66g2sWzEJiYRKQCKgAwRWQ4MVoubAmRGILaIuHl5UVR8EEdLHJOV7ErisgW5jq0P0d03wO66Vj5QPt+R8w8qmJKOOy3ZSm6YmDZaH8Qa4H3ALOBbnEsQbcDnwhuWiWWri3N5Zl8DR093Mnt6ZL9L7KlvpW9AHRvBNEhEKC2wkhsmto2YIFT1l8AvReQtqvr7CMZkYtxqfz/Eluombr00cp39cK6DOhrKfZTmu3ly9/GoXEjJmGAE00g7S0SmiM8DIlIhIm8Oe2QmZi3My2aGO82RfoiK2mZmT88gz50W8XMPV5afTXNnH43tPU6HYsyEBJMgPqCqZ4A3AzOA9wP3jnWQiPxcRE6KyO4h2+4SkToR2eG/rR3h2GtEpEpEqkXkziB/FxMlRITVxR5eOHQKr1cjdl5VZXtNc1RcPYBvdTmAAydsJJOJTcEkiMFr47XAL1R155Bto3kQuCbA9u+o6jL/bdOrTibiAu4HrgUWAbf56z+ZGFJe7OFURy/7I9hJW9fSxcm2HkfnPww1ONTV+iFMrAomQWwXkafwJYjNIuIGxlwVRlWfA05PIKZLgWpVPayqvcBvsEWLYs7gMqSRHM1UURs9/Q8Anuw0crNSbSSTiVnBJIgPAncCr1HVTiAVXzPTRH1MRF72N0EF+ksuAo4OeXzMvy0gEfmwiGwTkW2NjVbaIFoUTE1nYV4WWw5FMEHUNJOR4uICf9NONCjNd3PgpCUIE5uCSRCrgWzgYhG5Et9M6onWMP4RsBBYBhzHN3x2uEDNVyM2ZKvqT1V1paquzMtztkCcOd/qYg8vHj5Nb39kliGtqG1m6eypJLucmyA3XFmBmwMn2lCNXF+MMaESzF/SuiG3fwceA+6ayMlUtUFVB1TVC/wnvuak4Y4Bs4c8ngXUT+R8xlmrij109Q1QWdsc9nN19w2wt/5M1DQvDSrNd9PROxCw1pYx0W7MBKGqNwy5vQm4CGiYyMlEZOaQh/8E7A6w20tAiYjMF5FU4FZ8xQJNjLl8QS5JEpl+iJePtdLv1ShMENkAVvrbxKSJXIsfw5ckRiUiD+Gr2VQmIsdE5IPAN0Rkl4i8DFwFfNq/b6GIbALfehPAx4DNwD7gYVXdM4E4jcOmZqRw8awctkSgcN/gCnLL50TXCm4lgyOZbKiriUHBLDn6fc71ASTh6z/YOdZxqnpbgM0BV6FT1Xp8o6QGH28CXjUE1sSe8uJcfvzXw7R19+FOTwnbeSpqm5nvySI32/kJckNNzUhh5tR0DtoVhIlBwVxBbAO2+28vAHeo6rvCGpWJG+XFHga8youHJzLiOTiqSmVtc9RdPQwqzbeaTCY2jXkF4a/JZMyErJgzjfSUJLYcauKNi/LDco6jp7toau+Nuv6HQaX52bxw+BQDXsWVZDWZTOwYrdz3LgIPLxVAVfXisEVl4kZ6iovXzJse1o7q7bW+q5PoTRBuevu91JzqYEFettPhGBO00a4gro9YFCaulRd7uPeP+zl5ppsZU9JD/voVNS1kpbrOLtITbQbjOtDQZgnCxJQR+yBUtUZVa/z7NAx5fJLgajEZA/gmzAFsDdNoporaZpbNyYna5pviGdmI2PKjJvYE00n9W86vvTTg32ZMUBbNnEJOZkpYyn939PSz73j0TZAbKjM1mTnTM62j2sScYBJEsr9oHgD++6nhC8nEm6QkYdVC3zKkoS45sfNYC17F8RXkxlIyw21F+0zMCSZBNIrIjYMPROQmwLkV6U1MKi/2cLy1m1eaOkL6upWDFVxnR3eCKCvI5pWmjojVpTImFIJJELcDnxORWhGpBe4APhzesEy8WR2m8t8VNc0szMtiamb4JuGFQmm+m36vhjxBGhNOwdRiOqSql+NbvGexqq5S1UPhD83EkznTMynKyQhpP4SqUlEbPSvIjWZwJJP1Q5hYEnQtJlVtV1X7320mZOgypAMhWob0laYOmjv7omYFudEs8GSTnCTWD2FiSvQUzjdxr7zEw5nufnbXtYbk9c6uIBcDCSI1OYl5niy7gjAxZcwEISKvqn4WaJsxY1m1MBcgZM1MFbXNuNOTKY6RyWdl+W4r+21iSjBXEC8Euc2YUXmy07igwM3WEC1DWlHTzLLZOSRF6QS54Urz3dSe7qSrd8DpUIwJyogJQkQKROQSIENElovICv/t9UBmxCI0cWV1sYeXjjTT3Te5D8m27j6qGtpiov9hUFlBNqpQfdJmVJvYMFotpjXA+/At+fktzpXXaAM+F96wTLwqL/bwwPOvsO1IM6tLPBN+nZ1HW1GN3gJ9gZTmnxvJtGTWVIejiV0bK+vYsLmK+pYuCnMyWLemjJuXFzkdVlwaMUH4y3z/UkTeoqq/j2BMJo5dOn86yUnC89VNk0oQ22uaEYFlUboGRCBzc7NITU6yfohJ2FhZx/pHXqarzzfhsK6li/WP7AKwJBEGwfRBzBKRKeLzgIhUiMibwx6ZiUtZacmsmDNt0v0QFbXNlM5wMyWMq9SFmitJKM7LpsqGuk7YvX/cfzY5DOrqG2DD5iqHIopvwSSID6jqGeDNwAzg/cC9YY3KxLVVxbnsqmulpbN37J0D8Hp9K8itmBs7Vw+DygrctvzoBD25+wQnznQHfK6+pSvC0SSGYBLEYN/DWuAXqroTK/dtJmF1sQdVeGGC5b8PN7Vzpruf5THU/zCoNN9NfWs3Z7r7nA4lZrR29fGZh3dw+6+2k+IK/NFTmJMR4agSQzAJYruIPIUvQWwWETfnl/8OSER+LiInRWT3kG0bRGS/iLwsIo+KSMCvgCJyRER2icgOEdkW7C9jYsPS2TlkpbrYMsFmpu01zUBsdVAPKivwzdmwq4jgbKlu4trvPsf/7qjnE1cXc+8tS8hIcZ23jytJWLemzKEI49uYa1IDHwSWAYdVtVNEcvE1M43lQeAHwH8N2fY0sF5V+0Xk68B6fMX/ArlKVa1qbBxKcSVx2YJctlRP7AqioqaFnMwUFniyQhxZ+JXM8I9kOtHOJXOnOxxN9OrqHeDrT+7nwa1HWODJ4ne3X3H2itGVlHR2FFN2WjJtPf0xMxcm1gSTIBRfob7rgXuALGDMdSNV9TkRmTds21NDHv4deGuwgZr4Ul7s4U/7T1LX0kXROJsHKmqbWR5DE+SGKsrJICvVZSOZRrHzaAuffngHhxs7eN+qedxxzQVkpJ67arh5edHZEUt9A17e8ZMX+Pwju1g+O4fZ022KVigF08T0Q+AK4Db/4zbg/hCc+wPAH0d4ToGnRGS7iIxaWlxEPiwi20RkW2NjYwjCMpEw0fLfrV19HDzZHpPNS+BbPKnESm4E1Dfg5dtPH+CWH22lq3eAX33wMu66cfF5yWG4FFcS37t1OQCf+E0lfQO23kYoBZMgLpjM360AABhESURBVFPVjwLdAKrazCRXlBORzwP9wK9H2KVcVVcA1wIfFZErR3otVf2pqq5U1ZV5eXmTCctEUGl+Np7stHEniMpaf/9DDM2gHs5qMr3awYY2bvnhVu579iA3LS3kyU9dGfQ8mdnTM/nqLUuorG3hu88cCHOkiSWYBNEnIi583+oRkTyC6KQeiYi8F19z1Tt1hPUnVbXe//Mk8Chw6UTPZ6KTiFBe7OuHGM8ypBW1LSSJr6M7VpUWuGlq76WpvcfpUBzn9SoP/O0w133/eY41d/Kjd67g2+9YxtSM8c1vuWFpIe9YOZsf/uUQW8Ow9nmiCiZB3IfvQ3qGiHwFeB746kROJiLX4OuUvlFVO0fYJ8s/UgoRycI3/2J3oH1NbCsv9tDU3sOBhuBrE1XWNlNWMIXstGC6z6JTab5vJFOiX0Uca+7knx/4O19+Yh9XlnjY/OkruXbJzAm/3pduXMR8Txaf+p8dnLLkGxLBrCj3a+DfgK8Bx4GbVfW3Yx0nIg/hq/paJiLHROSD+EY1uYGn/UNYf+zft1BENvkPzQeeF5GdwD+AJ1T1yQn8bibKlfv7IYIt/z3gVSprW1gRQ+U1Ainz12RK1MWDVJXfbjvKNd/9G7uOtfKNt1zMf75nJTPcY459GVVmajLfv205LZ19rPvdy+O6MjWBjfk1TESmAyeBh4ZsS1HVUWf6qOptATb/bIR96/HNs0BVDwNLx4rLxL6inAzme7LYUt3EB1fPH3P/gyfbaO/pj6kKroHkudPIyUzhQAJWdW1q72H9I7t4em8Dl86fzrfetjSkI48WF07lc2sv4K7H9vLg1iO8v3zs/1dmZMFcp1cAs4FmfDOoc4DjInIS+JCqbg9jfCbOlRfn8mhFHX0DXlJco1/QVtT4V5CL0RFMg0SE0nx3wl1BPLn7BJ9/dBdt3f18fu2FfHD1/LAMVX7vqnn87WATX9u0n0vnT2dxoVXOnahg+iCeBNaqqkdVc/GNLHoY+Ai+IbDGTNjqYg8dvQPsPNoy5r4Vtc1Mz0plbm7sj3Uvy3dT1dCWEM0gZ7r7+OzDO7n9V9spmJrO459YzYeuXBC2eSwiwoa3LSUnM4WPP1RJZ29/WM6TCIJJECtVdfPgA/9ktytV9e+ALT1qJuXyBbmIBNcPUVHTzIo5OYjE3gS54Urzs2nr7h+x+Fy82FrdxDXfeY6NO+r4xNXFPPqR8rPrYoTT9KxUvvuOZbzS1MHdf9gb9vPFq2ASxGkRuUNE5vpv/wY0+4e+2qwUMyk5maksKZrK1jHKbjR39HK4qSOm5z8MdXbxoDhtZuruG+Dux/bwzw+8SFqKi9/dfgWfeXMZqcnBfOSExqpiDx95/UL+Z9tRHttZH7HzxpNg/rX+Gd+qchv9t9n+bS7g7eELzSSK8mIPFbXNdPSM3BRQeTR2C/QFMpggDo5jiG+s2Hm0hevu+xu/2HKE914xl02feK1jlXc/9cZSls/J4XOP7OLo6YAj680oRk0Q/quE76rqx1V1uf/2cVVtVNVeVa2OUJwmjpUv9NDvVf7xyukR99le04wrSbg4TpbqnJaVygx3GlVxNBeib8DLd/ylMjr9pTLuvumiUUtlhFuKK4n7rBTHhI2aIFR1AMgTkUmV1jBmNCvnTSM1OWnUfoiKmhYunOkmMzV2J8gNV1YQPyU3qk/6SmV879mD3DjOUhnhZqU4Ji6Yv7YjwBYR+QPQMbhRVb8drqBMYklPcfGaedNGrMvUP+Bl57EW3nbJrAhHFl4lM9z89z9q8Ho1JivTgq9UxoNbj/D1J/eTmerih+9cwdpJzIYOlxuWFvL8wSZ++JdDlC/0sKo4OpJXtAumD6IeeNy/r3vIzZiQWbXQw/4TbTS2vbpEQlVDG529A3HTQT2orCCb7j4vR5tjs228rqWLdz7wIvc8vpfVxb5SGdGYHAYNLcVxumNiy90mmjGvIFT17kgEYhLb6mIPGzZXsfVQEzctKzrvuYoYXkFuNIMd1Qca2pmbG/2LH22srDu7UM/UzBS6evpJdiXx9bcs4e0rZ0f98OPBUhz/dP9W1v12Jw+8d2XUx+y0Ma8gRCTPv1ToJhH50+AtEsGZxHFR0VSmpCcHHO5aUdtCnjuNWdPia93hkrMJIvr7ITZW1rH+kV3UtXShQEtnH31e5TNvLuUdr5kTMx+0iwunsn7tBTy7/yS/3HrE6XCiXjBNTL8G9gPzgbvx9Um8FMaYTAJyJQmrFnp4vrrpVbOLK2rjZ4LcUNlpycyalhETcyE2bK6iq2/gvG1ehZ8/f8SZgCbhfavmcfUFM/jqpv3srT/jdDhRLZgEkauqPwP6VPWvqvoB4PIwx2USUHlxLnUtXdScOtcm39TeQ82pzrhrXhpUGiOLB9W3dI1rezQTETa89WJ/KY4KK8UxiqAWDPL/PC4i14nIcnwT54wJqcHy31sOnRvNdLb/Ic46qAeV5rs51Nge9ePz89yBq+oUjnM98WiRm53Gd9+xjMNWimNUwSSIL4vIVOCzwL8CDwCfCmtUJiHN92RRODX9vOGuFbUtpLiEJUXxMUFuuLKCbPoGlJpTHWPv7BBVDbhAU0aKi3VryhyIKDSsFMfYglkw6HFVbVXV3ap6lapeAiyMQGwmwYgIq4o9bD10Cq/X1w9RUdvMosKppKc4Nxs3nM7VZIrekhtP7DrO4aYO3nrJLIpyMhB8a3l87ZYl3Ly8aMzjo5mV4hjdRCtnfSakURjjt7rYQ0tnH3uPn6FvwMvLx2J/BbnRLMzLJkmI2pIbXb0DfPWJfVw4cwpff8vFbLnzal659zq23Hl1zCcHsFIcY5logoiv4SQmaqwqzgV85b/3HT9Dd583bjuowTeLfF5uVtQuHvTjvx6ivrWbu25YhCtGZ3uPZWgpju89c9DpcKLKRAvbxP8qJ8YRM9zplOZns6W6iXR/aehYX2J0LNE6kqmupYsf//UQ1108k8sW5DodTljdsLSQvx1s5P6/VLOqOJdVC60UB4xyBSEibSJyJsCtDSiMYIwmwZQXe3jpyGleOHyKginpMTtSJlilBW6OnOqge9g8A6d9ddM+ROBzay90OpSIuOvGxcz3ZPFpK8Vx1ogJQlXdqjolwM2tqmNeeYjIz0XkpIjsHrJtuog8LSIH/T8DfjUUkWtEpEpEqkXkzon9aiZWuZKE7j4vm/c00NLVy8bKOqdDCquyfDdehUON0dNR/ffDp3ji5ePc/rqFFMV5gh40WIqjuaOPf/vdzoRYDnYs4Vze6UHgmmHb7gSeVdUS4Fn/4/P416C4H9/a14uA20RkURjjNFFkY2Udv/p7zdnH3X1e1j+yK66TRGl+NhA9JTcGvMrdj+2lKCeDf7kysQYsDpbieGafleKAMCYIVX0OGL4CzE3AL/33fwncHODQS4FqVT2sqr3Ab/zHmQSwYXMV3X3njyTp6htgw+YqhyIKv3meLFJcEjVDXR/6Ry37jp/hc2svdHSxH6dYKY5zIrdArE++qh4H8P+cEWCfIuDokMfH/NsCEpEPi8g2EdnW2NgY0mBN5MVTSYdgpbiSWJiXzcEouIJo7ezjW09Vcdn86axdUuB0OI6wUhznRDpBBCPQWLoRGwNV9aequlJVV+bl5YUxLBMJI3VIx31Hdb47KuZCfOeZA7R29XHXjYvjrjjieORmp/EdfymOex5L3FIckU4QDSIyE8D/82SAfY4Bs4c8noVv0SKTANatKSNj2KzpWC/pEIyyAjfHmrto73Hu22rViTb+399r+OfL5nDhzCmOxREtyos9/N/XLeQ3Lx3l8ZcT8yMo0gniD8B7/fffC/xvgH1eAkpEZL5/Lexb/ceZBHDz8iK+dsuSuCvpMJaSGb6OaqeamVSVex7fQ3ZaMp99U3wn4/H49JtKWTY7h/UJWoojbCvAi8hDwOsBj4gcA74E3As8LCIfBGqBt/n3LQQeUNW1qtovIh8DNgMu4OequidccZroc/PyorhPCMOVFZxbPGi5AzPHN+9pYEv1Ke66YRHTslIjfv5oleJK4vu3LWft9/7GJ39TycP/cgXJrmhsmQ+PsCUIVb1thKfeEGDfemDtkMebgE1hCs2YqDN7WibpKUkcaIj8SKbuvgG+smkvpfnZvOvyuRE/f7QbLMXx8Ycq+e4zB/nXOG/uHCpxUqExUSwpSRwrufGz51/h6OkuvnTD4oT6djweNywt5O0rZ3H/X6rZOmS9knhn/xuMiRKl+e6ILz96orWb+/9czZrF+WcXbDKBJWIpDksQxkSJ0vxsTrb10BzBD597/7iPfq/yheusWMFYMlOTue/WxCrFEbY+CGPM+AwuHnSgoS0i1VO315xm4456PnZVMbOnZ4b9fPHgoqKp3HntBdzz+F6W3fM0Z7r6KMzJYN2aMkcGVmysrGPD5irqW7rCEoddQRgTJc6OZDoZ/o5qr1e56w97KZiSzkeuSqx6S5M1LTOFJIHWrj4UX1l0J+qFbaysY/0ju6hr6QpbHHYFYUyUKJiSjjs9OSKLB/12+1F21bXyvVuXkZlqHwPj8c2nDuAd1rrU1TfAHb9/mUcimCRePHyKnv7AdctCdRVh/zOMiRIiQlkESm6c6e5jw+YqVs6dxo1LbWmX8RqpLlhPv5czXX0Ri2N4chgUyrplliCMiSIl+W7+uPs4qhq2Wkj3PXOQUx29PPj+SxO63tJEFeZkUBfgQ7goJ4ONHy2PWBzl9/4pYByhrFtmfRDGRJGy/GxaOvtobOsJy+tXn2znwa1HeMfK2VxUNDUs54h30VIvLBJx2BWEMVGk9GzJjXZmTEkP6WurKv/x+F4yUlwJNRs41Abb98M5eiha4rAEYUwUKfMPda1qaGN1SWgnrv1p/0n+eqCRL1x3IZ7stJC+dqKJlnph4Y7DmpiMiSK52Wl4slNDPpKpp3+A/3h8LwvysnjPFfNC+tomflmCMCbKlMwI/UimB7cc4cipTr54/SJSk+3P3gTH/qcYE2XKCtwcbGgLWSmHk23dfP9P1bzhghm8vizQKr/GBGYJwpgoU5rvpqN3IOAQxon4xpNV9PQP8IXrrd6SGR9LEMZEmbIC3+pyoSj9veNoC7/bfowPrJ7PfE/WpF/PJBZLEMZEmZLBkUwnJleTyVdvaQ957jQ+fnVJKEIzCcYShDFRZkp6CjOnpk/6CuLRyjp2HG3hjmsuIDvNRrSb8bMEYUwUmuziQe09/dz75H6Wzs7hligYr29ikyUIY6JQWYGb6sZ2BoaXDQ3SD/5UTWNbD3fdsIikJKu3ZCYm4glCRMpEZMeQ2xkR+dSwfV4vIq1D9vlipOM0xkml+W56+73UnOoY97FHmjr4+fOv8JYVs1g+Z1oYojOJIuINk6paBSwDEBEXUAc8GmDXv6nq9ZGMzZhoUTZkdbkFednjOvbLT+wlxSXccY3VWzKT43QT0xuAQ6pa43AcxkSV4hnZiIx/JNNfDzTyzL6TfOzqkpAX+zOJx+kEcSvw0AjPXSEiO0XkjyKyeKQXEJEPi8g2EdnW2NgYniiNibCMVBdzpmeOayRT34CXex7bw7zcTD6wel74gjMJw7EEISKpwI3AbwM8XQHMVdWlwPeBjSO9jqr+VFVXqurKvLy88ARrjANK893jShD/9UINhxo7+MJ1i0hLdo19gDFjcPIK4lqgQlUbhj+hqmdUtd1/fxOQIiKhrX1sTJQry3fzSlMHPf0DY+57qr2H7z5zgCtL83jDhVZvyYSGkwniNkZoXhKRAvGvhSgil+KL81QEYzPGcSX52fR7lVeaxh7J9M2nqujqHeCL1y+yZURNyDiSIEQkE3gT8MiQbbeLyO3+h28FdovITuA+4FYNVWlLY2JEWcFgyY3Rm5l217Xym5eO8t5V8yieMb4RT8aMxpH596raCeQO2/bjIfd/APwg0nEZE00WeLJJTpJR+yFUffWWpmem8ok3WL0lE1pOj2IyxowgNTmJ+Z4sDjSMPNT1Dzvr2VbTzLo1ZUzNSIlgdCYRWIIwJoqVFow8kqmzt5+vbdrPRUVTeNvK2RGOzCQCSxDGRLHSGW5qT3fS2dv/qud+9JdDnDjTzV03LMZl9ZZMGFiCMCaKlRVkowrVJ89vZjp6upOfPHeYm5YVsnLedIeiM/HOEoQxUaw0P/BIpq88sQ+XCHdee4ETYZkEYQnCmCg2NzeL1OQkDg65gtha3cSTe07wkdcvZObUDAejM/HOEoQxUcyVJJTMyD57BdE/4OXux/Yya1oGH7pygcPRmXhnCcKYKDe0JtN//6OWqoY2vnDdhaSnWL0lE16WIIyJcqX5bo63dlN7qpNvPXWAVQtzWbO4wOmwTAKwlcyNiXKnO3oAuHLDnwF4bYnH6i2ZiLArCGOi2MbKOv7rhfPX07rv2Wo2VtY5FJFJJJYgjIliGzZX0dPvPW9bV98AGzZXORSRSSSWIIyJYvUtXePabkwoWYIwJooV5gSe5zDSdmNCyRKEMVFs3ZoyMoYNZ81IcbFuTZlDEZlEYqOYjIliNy8vAnx9EfUtXRTmZLBuTdnZ7caEkyUIY6LczcuLLCEYR1gTkzHGmIAsQRhjjAnIEoQxxpiALEEYY4wJyBKEMcaYgERVnY4hZESkEagZc8fAPEBTCMOJZfZenM/ej/PZ+3FOPLwXc1U1L9ATcZUgJkNEtqnqSqfjiAb2XpzP3o/z2ftxTry/F9bEZIwxJiBLEMYYYwKyBHHOT50OIIrYe3E+ez/OZ+/HOXH9XlgfhDHGmIDsCsIYY0xAliCMMcYElPAJQkSuEZEqEakWkTudjsdJIjJbRP4sIvtEZI+IfNLpmJwmIi4RqRSRx52OxWkikiMivxOR/f7/I1c4HZOTROTT/r+T3SLykIikOx1TqCV0ghARF3A/cC2wCLhNRBY5G5Wj+oHPquqFwOXARxP8/QD4JLDP6SCixPeAJ1X1AmApCfy+iEgR8AlgpapeBLiAW52NKvQSOkEAlwLVqnpYVXuB3wA3ORyTY1T1uKpW+O+34fsASNiFCERkFnAd8IDTsThNRKYAVwI/A1DVXlVtcTYqxyUDGSKSDGQC9Q7HE3KJniCKgKNDHh8jgT8QhxKRecBy4EVnI3HUd4F/A7xOBxIFFgCNwC/8TW4PiEiW00E5RVXrgG8CtcBxoFVVn3I2qtBL9AQhAbYl/LhfEckGfg98SlXPOB2PE0TkeuCkqm53OpYokQysAH6kqsuBDiBh++xEZBq+1ob5QCGQJSLvcjaq0Ev0BHEMmD3k8Szi8DJxPEQkBV9y+LWqPuJ0PA4qB24UkSP4mh6vFpFfORuSo44Bx1R18Iryd/gSRqJ6I/CKqjaqah/wCLDK4ZhCLtETxEtAiYjMF5FUfJ1Mf3A4JseIiOBrY96nqt92Oh4nqep6VZ2lqvPw/b/4k6rG3TfEYKnqCeCoiJT5N70B2OtgSE6rBS4XkUz/380biMNO+2SnA3CSqvaLyMeAzfhGIfxcVfc4HJaTyoF3A7tEZId/2+dUdZODMZno8XHg1/4vU4eB9zscj2NU9UUR+R1QgW/0XyVxWHbDSm0YY4wJKNGbmIwxxozAEoQxxpiALEEYY4wJyBKEMcaYgCxBGGOMCcgShEl4IpIrIjv8txMiUue/3y4iP4xQDK8frBgrIjcmemVhEx0Seh6EMQCqegpYBiAidwHtqvpNB+P5Awk8YdNED7uCMGYEw77V3yUivxSRp0TkiIjcIiLfEJFdIvKkv0QJInKJiPxVRLaLyGYRmRngdd/mX0Ngp4g8F+D594nID/z380XkUf++O0VklX/7u0TkH/4rnZ/4S9cbE1KWIIwJ3kJ85b9vAn4F/FlVlwBdwHX+JPF94K2qegnwc+ArAV7ni8AaVV0K3DjGOe8D/urfdwWwR0QuBN4BlKvqMmAAeOekfztjhrEmJmOC90dV7RORXfhKszzp374LmAeUARcBT/vK8+DCVwp6uC3AgyLyML4ib6O5GngPgKoOAK0i8m7gEuAl/3kygJMT/7WMCcwShDHB6wFQVa+I9Om5OjVefH9LAuxR1VGX4lTV20XkMnxXIztEZNk44xDgl6q6fpzHGTMu1sRkTOhUAXmDazWLSIqILB6+k4gsVNUXVfWLQBPnl5wf7lng//qPc/lXdnsWeKuIzPBvny4ic0P8uxhjCcKYUPEvW/tW4OsishPYQeA1Ajb4O7d3A88BO0d52U8CV/mbtbYDi1V1L/AF4CkReRl4GnhVZ7gxk2XVXI0xxgRkVxDGGGMCsgRhjDEmIEsQxhhjArIEYYwxJiBLEMYYYwKyBGGMMSYgSxDGGGMC+v/9lbHHgF8dHwAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.plot(range(len(q6clusters_ase)), [max(x) for x in q6clusters_ase], 'o-')\n",
    "plt.xlabel(\"Time slice\")\n",
    "plt.ylabel(\"Largest cluster size\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As expected, the results are identical for both calculations!"
   ]
  }
 ],
 "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.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}