{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "liquid-concert",
   "metadata": {},
   "outputs": [],
   "source": [
    "import pyscal.crystal_structures as pcs\n",
    "import pyscal.csystem as pcm\n",
    "import pyscal.core as pc\n",
    "from ase import atoms\n",
    "from ase.io import read, write\n",
    "from ase.build import bulk\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "from scipy.spatial import Voronoi, voronoi_plot_2d\n",
    "from pychromatic import Multiplot\n",
    "from pychromatic.colors import accent\n",
    "import copy"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "id": "9eb6edbb-6e56-42fd-837a-fd13d6a3347d",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1500\n"
     ]
    }
   ],
   "source": [
    "nx = 5\n",
    "print(nx**3*12)\n",
    "atoms, box = pcs.make_crystal(structure=\"bcc\", repetitions=(nx,nx,nx))\n",
    "sys = pc.System()\n",
    "sys.box = box\n",
    "sys.atoms = atoms"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "id": "40eb9a52-1426-43b3-ad20-96dca00446e3",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "135 µs ± 1.64 µs per loop (mean ± std. dev. of 7 runs, 10000 loops each)\n"
     ]
    }
   ],
   "source": [
    "%%timeit\n",
    "real_atomdict = {\"positions\":copy.copy(sys.positions), \n",
    " \"ghost\":copy.copy(sys.ghost)}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "id": "66dc942b-7723-4d35-84c2-58acfa7fc273",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1.47 s ± 20.6 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)\n"
     ]
    }
   ],
   "source": [
    "%%timeit\n",
    "sys.find_neighbors(method='voronoi', cutoff=0.1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "id": "1dfe6b57-1bec-42db-8f04-374680cd7866",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "96000"
      ]
     },
     "execution_count": 20,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "len(sys.atom.voronoi.vertex.unique_positions)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "1ce6f7cf-67a7-487d-abd8-65b30e7bcd61",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "2.57 s ± 60.1 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)\n"
     ]
    }
   ],
   "source": [
    "%%timeit\n",
    "sys.find_neighbors(method='voronoi', cutoff=0.1, clean_vertices=False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "id": "ae768c7d-5dc9-4cf4-bc09-3d49370c3420",
   "metadata": {},
   "outputs": [],
   "source": [
    "natoms = [(5**3)*2, (7**3)*2, (10**3)*2, (12**3)*2, (15**3)*2, (16**3)*2, (18**3)*2, (20**3)*2,]\n",
    "\n",
    "t_voroc = [22, 60, 171, 359, 672, 857, 1200, 1600]\n",
    "t_voro = [7, 21, 64, 130, 290, 351, 480, 637]\n",
    "t_scipy = [7, 21, 78, 151, 377, 465, 671, 1000]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "id": "373debf6-50ee-4732-b93c-12bdf470c1d9",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAcMAAAEWCAYAAAAadfxCAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAABTvElEQVR4nO3deXyU5bn4/8+dZSbrZLKRlSUBggIBAgERlC2IWlHBpaUim3DqaS31nNpWbX/nHM/3nH7bs9hW8bT11CUouFDr9m3dYNhcWAQEg6BhCQlZICEhmezJzNy/P57JmECAIJnMTHK9X6+8yNzPMzNXIOTK/Tz3fV1Ka40QQggxkAX5OgAhhBDC1yQZCiGEGPAkGQohhBjwJBkKIYQY8CQZCiGEGPBCfB1Ab0lISNDDhg3zdRhCCCH81N69e89orRO7O9ZvkuGwYcPYs2ePr8MQQgjhp5RSxRc6JpdJhRBCDHiSDIUQQgx4kgyFEEIMeP3mnmF32tvbKS0tpaWlxdehiD4WFhZGeno6oaGhvg5FCBEA+nUyLC0tJTo6mmHDhqGU8nU4oo9oramurqa0tJSMjAxfhyOECAD9+jJpS0sL8fHxkggHGKUU8fHxckVACNFjXkuGSqnnlFKVSqmD54yvVkp9pZT6Qin1n53GH1VKHXUfu7HT+CSlVIH72JPqMjObJMKBSf7dhQh8NpuN+5Ys4Ya5c7lvyRJsNpvX3subl0nzgaeAFzoGlFKzgduBcVrrVqXUIPf4aGARMAZIBTYppbK01k7gD8D3gJ3AO8BNwLtejFsIIYSP2Ww2nlmzhpUZGWRlZlJot/PMmjUA5OXl9fr7eW1mqLXeDtScM/x94Nda61b3OZXu8duBV7TWrVrrIuAoMEUplQJYtNY7tNF48QVggbdi7i+ioqJ8HcIl5efn88Mf/tDXYQgh/NT6/HxWZmQw2molJCiI0VYrKzMyWJ+f75X36+t7hlnA9UqpXUqpbUqpye7xNOBkp/NK3WNp7s/PHT9PdXU1eXl5no+1a9dednB9OSUPJE6ns0fnPfbYY+R76RtVCDFwaIeDkooKsiyWLuNZFgsnKyq88p59nQxDgFhgKvBTYIP7HmB3N3j0RcbPEx8fj81m83wsW7bssgLrmJLfk5DAs9Onc09CAs+sWXNFCfHEiRNcddVVLFu2jHHjxnHXXXfR1NSEzWZj4cKFnvM2btzIHXfcgdPpZPny5YwdO5bs7Gx++9vfAnD06FHmzp3L+PHjmThxIseOHaOhoYG8vDwmTpxIdnY2b7311kVjefjhh/n973/vefzYY4/x+OOPo7Xmpz/9qec9X331VQC2bt3K7Nmzueeee8jOzqalpYUVK1aQnZ1NTk4OW7Zs6fHfw3vvvcfEiRMZP358t5c3qqqquPPOO5k8eTKTJ0/m448/BmD37t1MmzaNnJwcpk2bxldffQUYs8o77riDm266iZEjR/Kzn/2sx7EIIfyX1pqmzz+n8qmnSA0Pp9Bu73K80G5ncEqKV967r7dWlAKvuy957lZKuYAE9/jgTuelA+Xu8fRuxntd5yk50GVKfiXXp7/66iueffZZpk+fzn333cfvf/97HnroIR544AGqqqpITEzk+eefZ8WKFezfv5+ysjIOHjTWHNXW1gKwePFiHnnkERYuXEhLSwsulwuTycQbb7yBxWLhzJkzTJ06ldtuu+2CC0cWLVrEP/zDP/CDH/wAgA0bNvDee+/x+uuvs3//fg4cOMCZM2eYPHkyM2bMAIxkdPDgQTIyMnj88ccBKCgo4Msvv2TevHkUFhYSFhZ20a+/qqqKv/u7v2P79u1kZGRQU3PulXN48MEH+cd//Eeuu+46SkpKuPHGGzl8+DBXXXUV27dvJyQkhE2bNvHzn/+cv/zlLwDs37+fzz77DLPZzKhRo1i9ejWDBw8+77WFEIGhtagI+8aNtFdUEJqczD333MOzr75q3DO0WCi023m2qIhVq1d75f37Ohm+CcwBtiqlsgATcAZ4G3hJKfUbjAU0I4HdWmunUqpeKTUV2AUsBdZ4I7CTFRVkZWZ2GcuyWDhZUHBFrzt48GCmT58OwL333suTTz7JT37yE5YsWcK6detYsWIFO3bs4IUXXqC+vp7jx4+zevVqbrnlFubNm0d9fT1lZWWemWRH8mlvb+fnP/8527dvJygoiLKyMk6fPk1ycnK3ceTk5FBZWUl5eTlVVVXExsYyZMgQfvvb3/Ld736X4OBgkpKSmDlzJp9++ikWi4UpU6Z49ul99NFHrHZ/E1511VUMHTqUwsJClFIsWbIEgFOnTmEymfjd734HGLPtnTt3MmPGDM/rxMXFnRfbpk2bOHTokOex3W6nvr6euro6li1bxpEjR1BK0d7e7jknLy+PmJgYAEaPHk1xcbEkQyECUHtlJfaNG2k9epTgmBisCxcSnp3NTUoROmgQ6/PzOVlQwOCUFFatXu2VxTPgxWSolHoZmAUkKKVKgX8BngOec2+3aAOWuWeJXyilNgCHAAfwgHslKRiLbvKBcIxVpF5ZSTo4JYVCu90zM4TemZKfO1PreLxixQpuvfVWwsLCuPvuuwkJCSE2NpYDBw7w/vvv8z//8z9s2LDBk1jOtX79eqqqqti7dy+hoaEMGzbskvvq7rrrLl577TVOnTrFokWLAOOyxIVERkZ6Pr/QednZ2ezfvx8wLr0OGzaM5cuXd3nepbY5uFwuduzYQXh4eJfx1atXM3v2bN544w1OnDjBrFmzPMfMZrPn8+DgYBwOx0XfQwjhX5x2O/Vbt9K0fz/KZMIydy6R11yDCvk6LXWsAekL3lxN+l2tdYrWOlRrna61flZr3aa1vldrPVZrPVFrvbnT+b/UWg/XWo/SWr/baXyP+/zhWusf6ov99L4Ci5cv59miIg7V1uJwuThUW8uzRUUs7vSD/ZsoKSlhx44dALz88stcd911AKSmppKamsq///u/e5LHmTNncLlc3Hnnnfzbv/0b+/btw2KxkJ6ezptvvglAa2srTU1N1NXVMWjQIEJDQ9myZQvFxRfsTOKxaNEiXnnlFV577TXuuusuAGbMmMGrr76K0+mkqqqK7du3M2XKlPOeO2PGDNavXw9AYWEhJSUljBo16pLvee2117Jt2zaKiooAur1MOm/ePJ566inP447kWldXR1qasV5KFuYI0T+4Wluxb95M5Zo1NB04QOQ115D0ox8RNX16l0TY1/p1ObbL0fHbR29Pya+++mrWrl3L/fffz8iRI/n+97/vObZ48WKqqqoYPXo0AGVlZaxYsQKXywXAr371KwBefPFF7r//fv75n/+Z0NBQ/vznP7N48WJuvfVWcnNzmTBhAlddddUlYxkzZgz19fWkpaWR4p7xLly4kB07djB+/HiUUvznf/4nycnJfPnll12e+4Mf/IC///u/Jzs7m5CQEPLz87vMzi4kMTGR//3f/+WOO+7A5XIxaNAgNm7c2OWcJ598kgceeIBx48bhcDiYMWMGf/zjH/nZz37GsmXL+M1vfsOcOXMu+V5CCP+lnU6a9u2jfutWXE1NhI8dS/ScOYTExvo6NACUlyZafS43N1ef29z38OHDXH311T6KyFhNOn/+fM+CmHP98Ic/JCcnh5UrV/ZxZAODr//9hRDGrZKWL7/EbrPhrK7GNHQolhtuwJTW7S45r1JK7dVa53Z3TGaGPjJp0iQiIyM9qzSFEKK/aTt5EvvGjbSdPElIQgJxixZhzsryy3KJkgy9aNiwYRecFe7du7ePoxFCiL7hqKnBbrPRcugQQZGRxMyfT0RODirIf3tDSDIUQgjRK5xNTTRs20bjnj2o4GCiZ84kcto0gkwmX4d2SZIMhRBCXBHd3k7Dzp00fPwxuq2NiJwcomfNIjg62teh9ZgkQyGEEN+Idrlo/vxz7Fu24LLbMWdlYZk7l9DERF+HdtkkGQohhLhsLceOYd+4Ecfp04SmphK7cCHmYcN8HdY3JsmwH4qKiqKhocHXYQgh+qH2U6eM8mnHjxNstRJ7552EjRnjlytEL4ckw05sNhv56/KpKK0gJT2F5fcu77NSQP7M6XQSHBzs6zCEED7krKvDvmULzQcOoMLDsdx4I5G5uT6tGtOb/Hedax+z2WyseWYNSd9JYtYfZ5H0nSTWPCMtnHqjhZMQInC5Wlqwb9rE6aeeovngQSKnTSNp9Wqipk7tN4kQZGbokb8unxHLRhB3tdFVIe7qOEYsG0H+OmnhdCUtnIQQgUk7nTTu2UPDtm24mpsJHzeO6NmzCenUzKA/kZmhW0VpBdaR1i5j1pFWKkqvrKvyuS2cPvroI0/bo3Xr1lFbW8uOHTu4+eabyczM9LRweu+997BYLN22cIqIiEBrzc9//nPGjRvH3LlzPS2cLqRzC6cDBw54Wjh99NFH3bZwAs5r4dTRqqlzCychRP+itab5iy+o/J//wf7ee4QkJ5Pwve8Ru3Bhv02EIDNDj5T0FGqP1HpmhgC1R2pJSZcWTpc6TwjRP7SWlGD/4APay8oIGTSIuHvuwTxiRMAvjukJmRm6Lb93OUfXHqXmcA0uh4uawzUcXXuU5fcuv6LXHegtnIQQ/s9x5gw1r7xC9fPP47Tbsd52G4n330/YyJEDIhGCzAw9Ou4L5q/LZ3/pflLSU1i9Slo4dfimLZyEEP7L2dBA/bZtNO3diwoNJXr2bCKvvZag0FBfh9bnpIWTF0kLJ9/y9b+/EP7K1dZG444dNHzyCdrhIGLSJKJnziS4062R/khaOPkhaeEkhOhr2uWiaf9+6rdswdXQQNjVV2PJyyMkPt7Xofmc15KhUuo5YD5QqbUee86xnwD/BSRqrc+4xx4FVgJO4Eda6/fd45OAfCAceAd4UAfIdFZaOAkh/IHWmtajR43yaVVVhKanE3v33ZiHDPF1aH7DmzPDfOAp4IXOg0qpwcANQEmnsdHAImAMkApsUkplaa2dwB+A7wE7MZLhTcC7XoxbCCH6jbbycqPB7okTBMfFEXv33YRdffWAWRjTU15Lhlrr7UqpYd0c+i3wM6BzyZTbgVe01q1AkVLqKDBFKXUCsGitdwAopV4AFiDJUAghLspRW0v95s00FxQQFBFBzM03EzFpEkpKK3arT+8ZKqVuA8q01gfO+a0kDWPm16HUPdbu/vzc8fNUV1d3Wfm5dOnSbrcICCFEf+Zqbqb+ww9p3L0blCLquuuImj6dIKkWdVF9lgyVUhHAL4B53R3uZkxfZPw88fHx59URPXz48GVGKYQQgUk7HDTu3k39hx+iW1oInzABy+zZBFssvg4tIPTlzHA4kAF0zArTgX1KqSkYM77Bnc5NB8rd4+ndjPc7e/bs4YUXXuDJJ5/0dShCiACitab54EHqN2/GWVuLecQIo8FuUpKvQwsofZYMtdYFwKCOx+77gbla6zNKqbeBl5RSv8FYQDMS2K21diql6pVSU4FdwFJgjbdi3GzbxPoXn6ek9DRD0pNYvGQFc/LmeuvtusjNzSU3t9vtL0II0a3WEyewb9xIe3k5IcnJxC9Zgjkz09dhBSSvlWNTSr0M7ABGKaVKlVIX3Fmutf4C2AAcAt4DHnCvJAX4PvAMcBQ4hpcWz2y2beK5p5/gh/OCePdfM/nhvCCee/oJNts2XdHrNjY2cssttzB+/HjGjh3Lq6++yqeffsq0adMYP348U6ZMob6+nq1btzJ//nzAaK+0ZMkS5syZw8iRI/nTn/4EwJIlS7q0alq8eDFvv/32FcUnhAg87VVVVL/8MtVr1+JqaMC6YAGJ3/ueJMIr4M3VpN+9xPFh5zz+JfDLbs7bA4w9d7y3rX/xeR5aYCUnMwKAnMwIHloAT734/BXNDt977z1SU1P529/+BkBdXR05OTm8+uqrTJ48GbvdTnh4+HnP+/zzz9m5cyeNjY3k5ORwyy23sGrVKn77299y++23U1dXxyeffMLatWu/cWxCiMDirK+nfutWmj77DGUyET13LlFTpqAGYPm03iYVaNxKSk+TPbTrb1XZQ8MpKT1+Ra+bnZ3NT37yEx5++GHmz5+P1WolJSWFyZMnA2C5wM3t22+/nfDwcMLDw5k9eza7d+9mwYIFPPDAA1RWVvL6669z5513EtKPmmsKIbrnam2l4ZNPaNyxA+10EjllClEzZhAcEeHr0PoN+UnqNiQ9iYLiZs/MEKCguJkh6Vd2EzorK4u9e/fyzjvv8OijjzJv3rwebXa9UOunJUuWsH79el555RWee+65K4pNCOHftMtF07591G/diquxkbAxY7DMmUNIXNylnywui7Rwclu8ZAWPv1nLZ8ebcDg1nx1v4vE3a1m8ZMUVvW55eTkRERHce++9/OQnP2Hnzp2Ul5d7GujW19fjcDjOe95bb71FS0sL1dXVbN261TOTXL58uafH4ZgxY64oNiGEf9Ja0/zll1T9/vfU/e1vhMTHk7ByJXF33SWJ0EtkZujWcV/wqRefp6T0OEPSk7jv/geveDVpQUEBP/3pTwkKCiI0NJQ//OEPaK1ZvXo1zc3NhIeHs2nT+Yt0pkyZwi233EJJSQn/9E//RGpqKgBJSUlcffXVLFiw4IriEkL4p7bSUqN8WkkJIQkJxC1ahDkrS8qneZm0cPJDjz32GFFRUfzkJz8571hTUxPZ2dns27ePmJgYH0QXOAL1318MTI6aGuybN9PyxRcERUYSPWsWERMnooLkAl5vkRZO/cSmTZu47777+PGPfyyJUIh+wtXURP327TR++ikqOJioGTOImjaNIGme3ackGfqhxx57rNvxuXPnUlJS0u0xIURg0e3tNOzaRcNHH6Hb2ojIySF61iyCo6N9HdqA1O+TodZarrUPQP3l8r/of7TWNH/+uVE+zW7HnJVllE9LTPR1aANav06GYWFhVFdXEx8fLwlxANFaU11dTZhU6Rd+pvX4ceo2bsRx6hShKSlYFy7EPGyYr8MS9PNkmJ6eTmlpKVVVVb4ORfSxsLAw0tPTL32iEH2g/fRp7Bs30nrsGMFWK9Y77iB87Fj5Jd2P9OtkGBoaSkZGhq/DEEIMUE67HfuWLTTv348KC8Mybx6RkyejpHKU35F/ESGE6GWu1lYaPvqIxp070VoTee21RF9/PUHd1CEW/kGSoRBC9BLtdNK0dy/127bhamoiPDub6DlzCLFafR2auARJhkIIcYW01rQcPozdZsNZU4Np2DAsN9yAyV05Svg/SYZCCHEF2k6epO6DD2gvLSUkMZG4e+7BPGKELI4JMJIMhRDiG3BUV2O32Wg5fJigqChibr2ViAkTpHxagJJkKIQQl8HZ2Ej9tm007d2LCgkhevZsIqdOJchk8nVo4gpIMhRCiB5wtbfTuGMHDR9/jG5vJ2LSJKJnziQ4KsrXoYle4LVkqJR6DpgPVGqtx7rH/gu4FWgDjgErtNa17mOPAisBJ/AjrfX77vFJQD4QDrwDPKil1pYQwotsNhvr8/M5WVFBenIyd8+ezcS6Olz19YRddRWWvDxCEhJ8HaboRd68uJ0P3HTO2EZgrNZ6HFAIPAqglBoNLALGuJ/ze6VUsPs5fwC+B4x0f5z7mkII0WtsNhvPrFnDPQkJPDt9OosTE3nh9dfZUVND/IoVxH3nO5II+yGvJUOt9Xag5pyxD7TWHW3ddwId9bJuB17RWrdqrYuAo8AUpVQKYNFa73DPBl8AFngrZiGEWJ+fz8qMDEZbrYQEBTHaamVVVhZvnziBecgQX4cnvMSX9wzvA151f56GkRw7lLrH2t2fnzt+nurqavLy8jyPly5dyrJly3ozXiFEP+dqbqakooKszMwu41kWC6UFBT6KSvQFnyRDpdQvAAewvmOom9P0RcbPEx8fj81m650AhRADina5aPrsM+o3byY1PJxCu53RnarGFNrtDE5J8V2Awuv6fEOMUmoZxsKaxZ0WwpQCgzudlg6Uu8fTuxkXQohe0XbyJGeeeYa6v/6VkMRE7l2+nGeLijhUW4vD5eJQbS3PFhWxePlyX4cqvKhPZ4ZKqZuAh4GZWuumTofeBl5SSv0GSMVYKLNba+1UStUrpaYCu4ClwJq+jFkI0T856+uxb9pE8+efExQdTeyddxI2ZgzzlCLYajVWkxYUMDglhVWrV3e5DSP6H29urXgZmAUkKKVKgX/BWD1qBja6SxXt1Fr/vdb6C6XUBuAQxuXTB7TWTvdLfZ+vt1a86/4QQohvRDudNO7cSf327Wink6jrriPq+uu7bJrPy8uT5DfAqP6yZS83N1fv2bPH12EIIfxYy9Gj2N97D0d1NeasLGJuvJGQuDhfhyX6iFJqr9Y6t7tjUoFGCNHvOc6exf7++7R89RXBcXHE3XMPYSNH+jos4UckGQoh+i1XWxsNH31EwyefoIKDiZ47l6ipU1HBwZd+shhQJBkKIfodrTUtX3xB3caNuOx2wseNwzJ3LsHR0b4OTfgpSYZCiH6l/fRp6t59l7biYkKTk4m9806pHCMuSZKhEKJfcDU3U79lC4179hAUFkbM/PlE5ORIf0HRI5IMhRABzVM9xmbD1dJCRG4ultmzCQoP93VoIoBIMhRCBKy2kyepe/dd2isqMA0dSszNNxOalOTrsEQAkmQohAg4F6oe4y7mIcRlk2QohAgY51WPuf56oq67rkv1GCG+CUmGQoiA0HL0KHXvvYdTqscIL5BkKITwa46aGuref5/WwkKC4+OleozwCkmGQgi/5Gpro+HDD2nYsQMVHIxl7lwipXqM8BJJhkIIvyLVY4QvSDIUQvgNqR4jfEWSoRDC51zNzdi3bKFJqscIH5FkKITwGe1y0bRvH/WbN+NqaSEyN5doqR4jfECSoRDCJ1pLSrC/+y7tp05J9Rjhc5IMhRB9qkv1GItFqscIv+C1ZKiUeg6YD1Rqrce6x+KAV4FhwAng21rrs+5jjwIrASfwI631++7xSUA+EA68AzyotdbeilsI4R3a4aBh1y4apHqM6KHNtk2sf/F5SkpPMyQ9icVLVjAnb65X3subM8N84CnghU5jjwA2rfWvlVKPuB8/rJQaDSwCxgCpwCalVJbW2gn8AfgesBMjGd4EvOvFuIUQvazlyBGjekxNDWGjRmGZN0+qx4iL2mzbxHNPP8FDC6xkD82koLiZx59+AsArCdFrS7W01tuBmnOGbwfWuj9fCyzoNP6K1rpVa10EHAWmKKVSAIvWeod7NvhCp+cIIfyco6aG6pdfpuall0Ap4hYvJm7RIkmE4pLWv/g8Dy2wkpMZQUiwIiczgocWWFn/4vNeeb++vmeYpLWuANBaVyilBrnH0zBmfh1K3WPt7s/PHT9PdXU1eXl5nsdLly5l2bJlvRi6EKKnzqsec8MNRF5zjVSPET3T3kBJ6Smyhw7vMpw9NJyS0uNeeUt/WUDT3Z1zfZHx88THx2Oz2Xo1KCHE5dFa03zwIPaNG3HV10v1GHF5mk9DxSY4s4MhiSYKipvJyYzwHC4obmZIundWHPd1MjytlEpxzwpTgEr3eCkwuNN56UC5ezy9m3EhhJ9pP3WKuvfeM6rHpKQQd/fdmAYPvvQTxcCmNdQfg4qNcPYAqGBInMriJbN4PH8tDy0wZoQFxc08/mYt993/oFfC6Otk+DawDPi1+8+3Oo2/pJT6DcYCmpHAbq21UylVr5SaCuwClgJr+jhmIcRFuJqbsW/eTNPevQSFh0v1GNEz2gk1+40k2FAEIZGQdjMkzQJTDHMyAbOVp158npLS4wxJT+K++x/07WpSpVQQMB4jUTUDX2itT1/iOS8Ds4AEpVQp8C8YSXCDUmolUALcDaC1/kIptQE4BDiAB9wrSQG+z9dbK95FVpIK4Rc6qsfYN29Gt7QQOXky0bNmSfUYcXHOFqj6BCps0HoGzIkw7LuQeC0Em7ucOidvrteS37nUxbbsKaWGAw8Dc4EjQBUQBmQBTcDTwFqttcv7oV5cbm6u3rNnj6/DEGJAaC0poe7dd3GcOoVp2DBibrpJqseIi2urg1Ob4fR2cDZB1HBIvQFix4Pqm6sISqm9Wuvc7o5damb47xj7/O4/d6O7eyXoPcASvt4uIYTox5x2u1E9pqDAqB5z112EjR4t1WPEhTWVGZdCz+wG7YK4HEiZC9HDL/3cPnTRZKi1/u5FjlUCv+vtgIQQ/kc7HDTs3GlUj3G5pHqMuDitwf4llH8AdYcgyASDrjeSYFiir6PrVk/vGd4NvKe1rldK/X/ARODftdb7vBqdEMLnzqsec+ONhMTG+jos4Y9cDqjeY8wEm0oh1AKDb4dBMyA0ytfRXVRPV5P+k9b6z0qp64Abgf/GuHx6jdciE0L4lKOmhrr336e1sJCQ+HjiFi8mbMQIX4cl/JGjGSq3G/cE22ohPAUyl0LCFAgK9XV0PdLTZNixsvMW4A9a67eUUo95JyQhhC9J9RjRY63VxqrQyo/A1QqWUZCxBKyj+2xRTG/paTIsU0o9jbGq9D+UUma8WNdUCNH3zqseM348lrw8qR4jztdwwqgUU73XeJyQCyk3QOQQn4Z1JXqaDL+N0S3iv7XWte7qMT/1XlhCiL7UfuoUde++S1tJiVSPEd3TLqg9COUbob4QgsMgJQ+S54A58Auv9ygZaq2blFJbgMFKqYnu4TPeC0sI0RdcTU3Yt2z5unrMrbcSMWGCVI8RX3O1Q9VOYybYcgpMsTDkLhh0HYT0nwILPV1N+m/AcuAYXxfK1sAc74QlhPAmqR4jLqm9AU5vg9NboL0eIgbDiJUQNwmC+t/948u5TDpca93mzWCEEN4n1WPERTWfhlM2o2Saqx2sYyFlHliyoB8XV+hpMjwIWPm6y4QQIsB0rh4TLNVjRGdaQ8Mx435gR+eIhGuMTfIRqb6Ork/0NBn+CvhMKXUQaO0Y1Frf5pWohBBXzGazsT4/n5MVFaTFxnJrcjLXJiURNWMGUdOnS/UYYSyKqfnsnM4RN0HSbDDF+Dq6PtXTZLgW+A+gAPB5UW4hxMXZbDaeWbOGlRkZZGVmUmi388yRI0TfcAPzZs/2dXjC17rtHLEIEqed1zlioOhpMjyjtX7Sq5EIIXrNuuefZ2VGBqOtVgBGW62sGjmSl/78Z+bdJhd0Bqy2Oji1xVgY09E5YuidEDsh4DbJ97aeJsO9SqlfYTTh7XyZVGqTCuFnHLW1nKyoIGt4164AWRYLJwsKfBSV8KmmMmNrxJndRlPduAnGJnk/6xzhSz1NhjnuP6d2GpOtFUL4mZbCQs6+8QapkZEU2u2emSFAod3O4JQU3wUn+panc8RGqPvCqBE66Dpjo3zYIF9H53d6uulebjII4ce0y0X95s00fPwxocnJ3LtqFc8+95xxz9BiodBu59miIlatXu3rUIW3uZxQ/WlAdo7wpYsmQ6XUvcBLF+pkr5QaDqRorT/yRnBCiEtz1tdz9i9/oa24mIiJE4m5+WbmhYQQHBlprCYtKGBwSgqrVq8mLy/P1+EKb+m2c8QSY4tEgHSO8KVLzQzjMbZU7AX2AlVAGDACmIlRku2Ry31TpdQ/AqswLrUWACuACOBVYBhwAvi21vqs+/xHgZUY3TN+pLV+/3LfU4j+qPXECc6+9hq6rQ3rggVEjB/vOZaXlyfJbyDo6BxR9bGxStQyCjLuBeuYAb8o5nIorfXFT1AqGOPe4HQgBWgGDgPvaq1LLvsNlUoDPgJGa62blVIbgHeA0UCN1vrXSqlHgFit9cNKqdHAy8AUIBXYBGRprZ2dXzc3N1fv2bPncsMRIiBprWn46CPqt2whJC6O2G9/m9BBch9oQGkoNi6FdnSOiM+F1MDuHOFtSqm9Wuvc7o5d8p6hO+lsdH/0lhAgXCnVjjEjLAceBWa5j68FtgIPA7cDr2itW4EipdRRjMS4oxfjESJguJqbOfvGG7QeOUL42LHEzJ9PkHlg7g0bcDo6R1RsBHv/6xzhSz1dTdprtNZlSqn/BkowZpkfaK0/UEolaa0r3OdUKKU6fs1NA3Z2eolS91gX1dXVXS4JLV26lGXLlnnryxDCJ9rKyjj75z/jrK8n5lvfIiI3V8qpDQSudjizy1gZ2o87R/hSnydDpVQsxmwvA6gF/uxeqHPBp3Qzdt613fj4eGw2W6/EKIS/0VrT9Omn1L3/PsHR0STcdx+mtPN+JxT9TbedI+6DuNx+2TnCl/o8GQJzgSKtdRWAUup1YBpwWimV4p4VpvB1UfBSoHOX0XSMy6pCDAiu1lbq/vpXmg8exDxyJLELF0qrpf6u284RNxiLY+RKgFf0tJ9hEvB/gVSt9c3uRS3Xaq2f/QbvWQJMVUpFYFwmzQP2AI3AMuDX7j/fcp//NvCSUuo3GAtoRgK7v8H7ChFw2isrObthA46aGqLz8oiaPl0ui/ZX0jnCp3o6M8wHngd+4X5ciLEN4rKTodZ6l1LqNWAf4AA+A/4XiAI2KKVWYiTMu93nf+FecXrIff4D564kFaI/ajpwgLq//hVlNhO/dCnmYcN8HZLwhnM7RwRHDNjOEb50ya0VAEqpT7XWk5VSn2mtc9xj+7XWE7wdYE/J1grRX2iHg7p336Vp3z5MQ4cSe+edBEdH+zos0ducre7OEZvcnSMSjFngAO4c4W1XtLXCrVEpFY974YpSaipQ10vxCSHcHDU1nP3zn2k/dYqo664jevZsVJBsnO5XzusckQFD7jSKZ8smeZ/paTL8Mca9u+FKqY+BROAur0UlxADU/OWX1L75JihF3He/S1hWlq9DEr2pqdy4FNrROSJ2PKTOk84RfqKnhbr3KaVmAqMwtjp8pbVu92pkQgwQ2unEbrPRuGMHoampxN59NyGduk2IANZt54jpxuVQ6RzhV3q6mjQY+BZG3dAQYJ5SCq31b7wYmxD9ntNuN4psl5QQkZtLzI03okJ8seNJ9CpP54hN0HQSQqMh/TZImimdI/xUT//X/T+gBaOodrcdLIQQl6f1+HHO/uUv6PZ2Yu+8k/CxY30dkrhSjmao/NDdOeKsdI4IID1Nhula63FejUSIAUJrTcP27dRv3UpIYiKxd99NaGKir8MSV6K1xtgkX/lRp84Ri6VzRADpaTJ8Vyk1T2v9gVejEaKfczY1Ufv667QeO0b4uHHE3HILQSaTr8MS35R0jug3epoMdwJvKKWCgHaMRTRaa23xWmRC9DNtpaXU/PnPuBobiZk/n4iJE6WaTCDSLqj9Aio+6NQ5Yg4k50nniADW02T4OHAtUKB7sktfCOGhtaZx1y7sGzcSbLGQsHIlppQUX4clLtd5nSOsxv7AQddL54h+oKfJ8AhwUBKhED1js9lYn5/PyYoKUi0WbktJYfb112O9/XYpsh0ANts2sf7F5ykpPc2QtEEsnj+ROZmnpHNEP9bTZFgBbFVKvQu0dgzK1gohzmez2XhmzRpWZmSQlZlJod3Os0ePEpOQwFxJhH5vs20Tzz39BA8tsJI9NJOC4mYef3Uj3D6eObf9A1iuks4R/VBPlzkVATbABER3+hBCnGPdc8+xMiOD0VYrIUFBjLZaWTliBC+tXevr0MSlOBpZn/9HHlpgJSczgpBgRU5mBA/dMYj1tjKIuVoSYT/V0wo0/+rtQIQIdNrppHHPHk6eOkXWiBFdjmVZLJwsKPBRZOKinG1Gy6Tq3VD7BSUVZ8ke2nUhTPbQcEpKj/soQNEXLpoMlVJPaa1/qJT6f3TTXV5rfZvXIhMiQGitaS0spG7jRpzV1aRZLBTa7YzuVFKt0G5nsCya8R/aCXWH4cynRvskVyuEWiF5DkPSaigobiYnM8JzekFxM0PSk3wX7wBls9nIX5dPRWkFKekpLL93OXl5eV55r0vNDJcCPwT+2yvvLkSAaz99mrr336etqIiQ+HjivvtdlkyfzrNPPWXcM3QnxmeLili1erWvwx3YtDb6BZ7ZZewLdNQbvQMTJkP8FLCMBBXE4qVWHn/6CR5aYMwIC4qbefzNWu67/0FffwUDis1mY80zaxixbASjRo6i9kgta55ZA+CVhHipZHgMQGu9rdffWYgA5mxooH7LFpo++wwVFoblppuIzM1FBQczNysLpZSxmrSggMEpKaxavdprv9GKS2iuMBLgmU+NvoEqBGLHGSXSrGPOK5M2J28uAE+9+DwlpccZkp7Effc/6BkXfSN/XT4jlo0g7mrjknXc1XGMWDaC/HX5PkmGiUqpH1/ooKwmFQONdjho2LmThg8/RDscRE6ZQvTMmedtl8jLy5Pk50utZ417gGc+NQployDmKki7BeJyLrkvcE7eXEl+fazJ2URFWwVlrWWUt5ZTfrKcUSNHdTnHOtLK/tL9Xnn/SyXDYCAKo+KMEAOW1pqWL77AvmkTzro6wkaNwnLDDYTEx/s6NNHB0QjV+4x+gfVHAA2Rw2Dot40yaaYYX0co3Jzayem200bSayunvLWcGkcNAApFQmgCsWmx1B6p9cwMAWqP1JKS7p1775dKhhVa6//T22+qlLICzwBjMRbm3Ad8BbyK0SbqBPBtrfVZ9/mPAisBJ/AjrfX7vR2TEBfSVlZG3fvv037yJCFJScTffjvmjAxfhyUAXG1w9nMjAdYeNBbGhCVB+nyInwzhsujF17TW1DnqKGtzz/hayznVdgonTgCigqNINaUyLmocqeZUUkwpmIPMZC7N9NwztI60UnuklqNrj7J6lXfuvV8qGXprRvgE8J7W+i6llAmIAH4O2LTWv1ZKPQI8AjyslBoNLALGAKnAJqVUltba6aXYhADAWVeH3WajuaCAoMhIYm69lYgJE1BB0oXAp7QT6r40EqBnJWgMJM82FsJEDpG9gD7U4mqhotV9udM962tyNQEQokJINiWTa8kl1ZRKmjmN6ODobmv0dtxmyF+Xz/7S/aSkp7B6lffuvauLVVhTSsVprWt69Q2VsgAHgMzO5d2UUl8Bs7TWFUqpFGCr1nqUe1aI1vpX7vPeBx7TWu/o/Lq5ubl6z549vRmqGKBcbW00fPQRDTt2gNZEXXstUdddR5DZ7OvQBi6toeEEVO+C6j1GWbTgcIibCAlTwJIlrZJ8wKVdVLZXemZ8ZW1lVLdXe47Hh8STZk4j1ZxKqjmVxNBEgpXvStgppfZqrXO7O3bRmWFvJ0K3TKAKeF4pNR7YCzwIJGmtK9zvW6GUGuQ+Pw2ja0aHUvdYF9XV1V1+Y1i6dCnLli3zQviiv9Ja07x/P/bNm3E1NBA+dizReXmEdNovKPpY8yn3StDd56wEnQLWsdIwt4/ZHfYuM75Tbado1+0AhAeFk2ZOY0zEGNLMaaSYUwgLCvNxxD3X09qkvf2eE4HVWutdSqknMC6JXkh31zvOm87Gx8djs9l6KUQx0LSeOIH9/fdpP3WK0PR04r7zHUzp6b4Oa2BqOwtn9hhJ8BusBBW9o83VRkVbhTHjcyfABmcDAMEEk2xKZkLUBGPWZ0rFGmIN6JZkvkiGpUCp1nqX+/FrGMnwtFIqpdNl0spO5w/u9Px0oLzPohX9mqOmBvvGjbR8+SXBMTFY77iD8LFjA/o/dUByNBr3/87sNnoEelaC3m0shJGVoF7l0i7OtJ/xrO4say3jTPsZtHveERsSy9CwoZ77fINMgwhRvkgf3tPnX43W+pRS6qRSapTW+isgDzjk/lgG/Nr951vup7wNvKSU+g3GApqRwO6+jlv0L66WFuq3baNx925USAjRc+YQNXUqKlQuu/UZVxucLei0EtQBYYMg/RZjIYysBPWaBmfD1zO+1nIq2ipo020AhAWFkWpKZVTEKONypymFiOCIS7xi4PNVal8NrHevJD0OrMDooLFBKbUSKAHuBtBaf6GU2oCRLB3AA7KSVHxT2umkae9e6rduxdXcTERODtFz5hAcFeXr0AYG7YS6r4wEePYzcLYYK0GTZhr3ASOHykrQXtbuaudU2ynPjK+8tRy70w5AEEEMMg0iOyqbVJOxyCUuJG5AXhm56GrSQCKrScXFaK1pPXoU+wcf4DhzBlNGBjHz5hGanOzr0Po/raHxhJEAq/dAux2CwzqtBB0lK0F7idaaGkfN11Vc2sqpbKvEhQuAmOAYUs3Gpc5UUypJpiRCB9AipG+8mlSI/qC9shL7Bx/QeuwYwXFxxC1ahNldP1R4UfMpdwL8FFoq3StBs90rQbNlJWgvaHI2ebY0dFzubHG1AGBSJlLNqUy1TPVsbYgKlisgFyLJUPRbzsZGo5j2vn0osxnLjTcSOXkyKth3+5z6vbazxuzvzG5oLAGUMfNLvcm9ErT/33vyFod2fF3CzD3rO+s4CxglzBJDE7kq4irPrC8+NJ4gmXH3mCRD0e9oh4OGXbuMYtptbUROnmwU046QH8Re4WiCmn3nrAQd6l4Jmgsmq68jDDhaa2odtV3u851uO+0pYRYdHE2qOZUJURNIM6eRbErGFGTycdSBTZKh6De01rQcPox940actbWYs7Kw3HADoQkJvg6t/3G1n1MT1L0SNO0W4zKorAS9LM3O5q87Nrg3tDe7mgEIVaEkm5KZbJls3O8zpREdEu3jiPsfSYaiX2grL8f+/vu0lZQQMmgQ8UuWYM7M9HVY/Yt2gf0rYzN8TcdKUIusBL1MTu2kqq3Kc5+vrLXM07EBICE0gazwrC4lzORyp/dJMhQBzWm3G8W0P//cKKY9fz4ROTlSTLu3aA2Nxe7u8LIS9HJprbE77V1Wd55qO4VDOwCIDIok1ZzaZWuDOUhq4PqCJEMRkFxtbTR88gmNH3+M1pqo6dOJuv56KabdW5pPu1eC7v56Jah1rNEdPnYsyP2pbrW6WqlorfBc6ixrK6PR2Qh83bFhYtREz6wvJjhGVjX7CUmGIqBorWk+cMAopl1fT9iYMVjmzpVi2r2hrbbTStBijJWgWZB6ozETlJWgXXSUMOu4z9dRwqxDXEgcGWEZnvt8g0yDfNqxQVycJEMRMFqLi41i2hUVhKalEXf33ZgGD770E8WFOZo7rQT9CmMl6BAYepe7JqjV1xH6hM1mI39dPhWlFaSkp7D83uVMmTmlS4PairaKLh0bUk2pXB1xtaeEWXiwFBQPJJIMhd9znD1rFNM+fJggiwXrwoWEZ2fL5aVvytVu1ASt3m38qR1gToS0b7lXgg7sqjw2m401f1rDiOUjGDVyFLVHavnt//6WzDOZJE9NJoggkkxJjI8a77nPFxsSK9+PAU6SofBbrpYW6rdvN4ppBwURPWsWkdOmESTFtHtks20T6198npLS0wxJT2LxnTcyZ7Q2ZoKelaAz3CtBhw3YlaAde/rKWssoayvjj2v/yIjlI4i7Og6AuKvjuGrFVZS8UsLDCx4myZTU7zo2CEmGwg9pl+vrYtpNTYRPmIBlzhyCo2VvVU9ttm3iuaef4KEFVrKHZlJQ3Mzja9fDvGTm5M02ukLEjIIBeA+ro09fWWuZ56PJ1QQYJczs5XasI61dnmMdaWV/2X7SzOf1FRf9hCRD4VdaOoppV1VhGjoUy403YkpJ8XVYgaW1mvX5v+ehBVZyMo1FLzmZETx0xyCe+qCdOd9b7tv4+pDWmrOOs10SX2V7padPX1xIHMPDh5NmTiPNnEZiaCIHBh+g9kitZ2YIUHuklpR0+T7szyQZCr/QXlVlFNM+epTg2Fhiv/1twq66Su7D9JSjGWr2QtVOqD9CSUUd2UO7Vt7JHhpOSelxHwXYNzq2NpS1fZ38Oiq5dBSunhYxzVO/s7s+fcvvXc6aZ9YwYtkIrCOt1B6p5ejao6xetbqvvxzRhyQZCp9yNjUZxbT37kWZTFjmzSNyyhQppt0TLifUHYSqXXD2gLskWhKk38aQtHUUFDd7ZoYABcXNDEnvP2XSOrcrKm0tpby1nKr2Ks+sLz40npHhIz2zvoTQhB5VcsnLywMgf10++0v3k5KewupVqz3jon+SZCh8QjscNO7eTf327ei2NiJyc4meNYtgKaZ9cR29Aat2GnsCHQ0QEgWDroPEqZ6FMIuXmnn86Sd4aIExIywobubxN2u57/4HffwFfHMtrhZj1tdxybOtzNOuyKzMpJpTyYrIIt2cfsVbG/Ly8iT5DTCSDEWf0lrT8uWXRjHts2cxjxxpFNNOTPR1aP6t5YxREu3MLmg57e4NOB4Sr4GYsRDUdSY9J28uAE+9+DwlpccZkp7Effc/6Bn3d1przrSf8WxmL20t7bKhPSE0gasiriLVlOqZ9ckldXElpNO96DNtFRVGMe3iYkISE7HceCNhw4f7Oiz/5WiC6r1wZifUHzXGokcaM8B+VhGmxdVCeWu553JneVu5Z9YXFhRGmimNVHOqMeszpxAWFObjiEUg8stO90qpYGAPUKa1nq+UigNeBYYBJ4Bva63Pus99FFgJOIEfaa3f90nQ4htx1tcbxbQPHCAoIoKYW24hYuJEKabdHZcDar8wEuDZz7++Dzj4dmM7RFjgt6PqUsastZzStlKq26uBrk1q083ppJpTiQ+Jl1mf8DpfXiZ9EDgMWNyPHwFsWutfK6UecT9+WCk1GlgEjAFSgU1KqSyttdMXQYuec7W30/jJJzR8/DHa5SJy2jSir7+eoDD5rb4LraHhhJEAqz8FR6NxHzDpekiYGvCtkZqdzV1Wd1a0VtCqWwGjjFmaOY2xkWOla4PwKZ8kQ6VUOnAL8Evgx+7h24FZ7s/XAluBh93jr2itW4EipdRRYAqwow9DFpdBa01zQQF2mw2X3U7Y6NFGMe3YWF+H5l8ueB9wKsSMOe8+YCBwaRdV7VVd9vV19OpTKAaFDmJ05GjSzemkmdOkjJnwG76aGf4O+BnQuaRIkta6AkBrXaGUGuQeTwN2djqv1D3WRXV1dZfVX0uXLmXZsmW9HLa4lNaSEqOYdnk5oampxN5xB+ahQ30dVp87rxTakhXG4hVHo/s+4K5O9wGzIHUexE2CkMAq7tzkbOqyurOitYI23QZARFAEaeY0xkWN8xSvNknrJ+Gn+jwZKqXmA5Va671KqVk9eUo3Y+et+omPj8dms11hdOKbcpw9i33TJloOHSIoOhrrggWEjxs3IH/r77YU2tO/g4rNzBlZ574PmGzcB0y4Bszxvg65R1zaRWV7ZZdZ31nHWcCY9SWZksiOyibNZOzrs4ZYB+S/vwhMvpgZTgduU0p9CwgDLEqpdcBppVSKe1aYAlS6zy8FOvfpSQfK+zRicUGu1lYaPvyQhp07jWLaM2caxbRNA3cGsP7F588vhbYAnvrrfub832VGAgyA+4CNzsYuia9zy6LIoEjSzGlMiJrgmfWFBkkBdRG4+jwZaq0fBR4FcM8Mf6K1vlcp9V/AMuDX7j/fcj/lbeAlpdRvMBbQjAR293HYAqO1zfr8fE5WVJCenMzdM2cysbbWKKY9frxRTNtiufQL9WfOFkpKT5E9tOuWkeyh4ZRUtsGw7/gosItzaieVbZVdLnnWOmoBurQsSjOnkWZKIyZEOrSL/sWfNt3/GtiglFoJlAB3A2itv1BKbQAOAQ7gAVlJ2vdsNhvPrFnDyowMsjIzKbTbeebNN2mbNIlv/f3fY0pN9XWIvtVUBqe3wZldDEk0+UUptO4a1HbcV29wNpw363NoBwBRwVGkmdOYGDWRNHMayaZkmfWJfk823YseuW/JEu5JSGC01eoZO1Rby/qqKp5ft853gfmSq93oDXh6G9QfM1aDxk9i81dhPLf2Zfc9w66l0PqqAozNZjuv2HRhfiHTvzOdqMlR1DnrAGPWl2xK9tTvTDOnYQm2yKxP9Et+ueleBJaTFRVkZWZ2GcuyWCgtKPBRRD7UUgWnt0PVJ0ZtUHMiDLkTEqdBaBRzRgBhg3xWCq3d1c4zLz7DiGVdG9RmLc9ix0s7WHrdUnLNuZ5ZnzSqFUKSoeiBxn37SAkPp9Bu7zIzLLTbGTxQeg1qJ5wtMGaBdYeAIIgbD4NmQMxVcE43hDl5c/s0+ZW1llHcWkxxSzHlreVUllYyZuSYLudZR1ppqGjgjsQ7+iQuIQKJJENxQa72dureeYfm/fu5Y9Iknv3sM1ZmZpJlsVBot/NsURGrVvfzHm9tZ6HyI+OjrRZMVki/FQZNB5Nvigg4tMNIfi1fJz8nThSKZFMyUyxT2J++XxrUCnEZJBmKbjlqaqjZsAHH6dNEzZjB7TNnEr1li7GatKCAwSkprFrdT3u8aRfUfWnMAs9+DrggZjQM+y7EZoPq28owDu2gvLXck/zKWsu6JL9cSy5Dw4Yy2DzYU8pML9HSoFaIyyALaMR5mg8fpvattyAoiNiFCwkbOdLXIfWN9gbjPmDlh9BSadQHTZwGSTMgrO9aTHUkv5KWEiP5tZV5Vnomm5IZYh5iJL+wwRft3nCx1aRCDEQXW0AjyVB4aKcTu81G444dRim1u+8mpNM9wn5Ja2g4ZiyIqd5rVIeJHgFJMyEuB/pgS4FTO43k12okv9LWUk/ySwpNYkjYEM/M70oa1gox0MlqUnFJzvp6zr72Gm0lJUTk5hJz442okH787eFoNuqDVm439ggGhxnd4pNmQMR5pW8v28VmZU7tpKKtwjPzK20t9VR2GRQ6iJyoHIaEDWGIeYgkPyH6SD/+aSd6qrWoiLN/+Qu6rQ3rHXcQkZ3t65C8p/Gke3P8bnC1QsRgyLgXEiYbCbEXdN7jN2rkKGqP1PLkn56ksKmQ2CmxnGw96Ul+iaGJjI8a70l+EcH9p2GvEIFEkuEAprWm4aOPqN+yhZC4OGKXLiV00KBLPzHQuNqMS6Cnt0FDEahQSMg1LoVGDuv1GqH56/LP2+M3cvlI3n3pXeZPms+4qHEMMQ9hSNgQIoMje/W9hRDfjCTDAcrV3MzZN9+ktbCQsDFjsN56K0HmftZUtfm0kQCrdoCzyegYP/RuSLwWQno3CbW52ihuKeZY8zHKT5YzauSoLsetI600VTTxvdTv9er7CiF6hyTDAaitooKzGzbgtNux3HQTkVOm9J/yWy4nnN1vJEH7V8Zm+NgcYxZoyerVWWBNew3Hmo9xrPkYxS3FOHESqkKJSYuRPX5CBBhJhgOI1pqmffuoe/ddgiMjSVixAlN6uq/D6h2tNcaWiMqPoN0OpjijX2DidDDF9MpbOLSDkpYSTwLs6OAeFxLHpOhJDA8fzuCwwUxcOlH2+AkRYCQZDhCu9nbq/vY3mg8cwDx8ONY77iA4IsAXa2gX1B6Cym1GqTQA61hjRah17Hkl0r6JOkedJ/mdaDlBu24nRIUw1DyU3OhchocPJza0ayWajlWj+evy2V+6n5T0FFav6qcFCoToJ2Sf4QDgqK42qslUVhI1cybRM2aggq48UfhMux0qPzG2RbRWQ2g0JF5nbI0IS7iil3ZqJ6WtpRxrPsbR5qOcaT8DQExwDCMiRjA8fDhDzUOlpZEQAUj2GQ5gzYcOUfvWW6jgYOIWLyZsxAhfh/TNaA31R4zN8TX7jMLZliwYcgfEToCgb/6t3OBo4FiLkfxONJ+gVbcSRBBDwoYwPmo8w8OHEx8S33/uqwohziPJsJ/STif2TZto3LmT0LQ0o5pMTO/cO+tTjiY4s9NIgs0VEBxhLIZJmgnhyd/oJV3aRXlrOUebj3Ks5Rin204DEB0czdWRVzM8fDjDwoZ56nwKIfo/SYb9kNNuN6rJnDxJ5OTJWG68ERXct8Wlr1jDCXeJtE+NfYKRwyBzKcRPhmDTZb9ck7OJ483HOdZ8jOMtx2l2NaNQpJvTmWWdxYjwESSGJsrsT4gBSpJhP9O5mkzsnXcSPnasr0PqOWcbVO82kmBjMQSZIGGK0TMwauhlvZTWmlNtpzyLX8raygCICIpgRLhx7y8jLEPKnQkhAB8kQ6XUYOAFIBlwAf+rtX5CKRUHvAoMA04A39Zan3U/51FgJeAEfqS1fr+v4/Z3XarJxMcTu2wZoYl912nhijSVG9siqnaAsxnCU2HYIkiYCiE9T1YtrhaKmos8CbDR1QhAqimV62OuZ0T4CJJNyTL7E0KcxxczQwfwkNZ6n1IqGtirlNoILAdsWutfK6UeAR4BHlZKjQYWAWOAVGCTUipLa+30Qex+ydXczNk33qD1yBHCx44l5tZbCTJd/qXEPuVyQM1nxiywvhBUCMRNNLZFRI/o0eZ4rTVV7VWelZ+lraVoNGFBYWSGZTI8fDiZ4ZlS8kwIcUl9ngy11hVAhfvzeqXUYSANuB2Y5T5tLbAVeNg9/orWuhUoUkodBaYAO/o2cv/UVl5uVJOpryfm5puJmDzZv2c+LWfcs8CPob0ezAkweCEMmgahlks+vc3VxomWE8bil+Zj1DvrAaPV0bWWaxkePpw0cxpBvbDHUAgxcPj0nqFSahiQA+wCktyJEq11hVKqo2J0GrCz09NK3WNdVFdXd9nUvHTpUpYtW+alyH1Pa03T3r3Uvfee/1eT0S6oLTBmgbVfGGOx44xZYMzoi26O11pT46jxzP5OtpzEiROTMpERnsHwsOEMDx9OdEh0H30xQoj+yGfJUCkVBfwF+Aettf0is5nuDpxXKSA+Ph6bzdaLEfoXm83G+vx8TlZUkJ6czMKxY8ltbcU8YgSxCxcS5I/VZNrqjPJolR9C21kIjYG0bxmb481xF3xau6udktYSz+yv1lELQEJoArkWo+rLYPNgglWArZAVQvgtnyRDpVQoRiJcr7V+3T18WimV4p4VpgCV7vFSYHCnp6cD5X0Xre/ZbDaeWbOGlRkZZGVmUmi388yHHxJ888186557/OuyqNZGgezT24yC2doFMVfDsG+DdTwEdZ/AattrOdpy1FP02qEdhKgQhoUN4xrLNQwPH441xNqnX4oQYuDwxWpSBTwLHNZa/6bTobeBZcCv3X++1Wn8JaXUbzAW0IwEdvddxL63Pj+flRkZjLZaARhttbIqK4uXdu3ilh/+0LfBdXA0GqtBT2+HltNGi6TkPBh0PYQnnXe6Uzs52XLSM/urdlQDEBsSy4SoCYwIH8GQsCGEKNn9I4TwPl/8pJkOLAEKlFL73WM/x0iCG5RSK4ES4G4ArfUXSqkNwCGMlagPDKSVpFprSioqyMrM7DKeZbFwsqCgT2PZbNvE+hefp6T0NEPSk1h87wrmXDPcmAVW7wHdDlGZMHwFxE+Cc+p31jvqPff+TrScoE23EUwwQ8KGkBOdw4jwEcSFXvjyqRBCeIsvVpN+RPf3AQG6Leuvtf4l8EuvBeWHtNa0fPUV9du2kRoeTqHd7pkZAhTa7QxO6bv+eJttm3ju6Sd4aIGV7KGZFBQ38/jT/w3H4piTk2A0zE2aAZFfX9F2aRdlrWWeBFjZblz5tgRbGBM5hhHhIxgaNhRTkJ9vAxFC9HtyDcrPaK1pLSykfutW2k+dIjgujkULF/Ls228b9wwtFgrtdp4tKmLV6r7rj7f+xed4aIGVnExjoU5OZgQPLUzkqXebmLPqPyE4DIBGZyPHm49ztPkoRS1FtLhaCCKIdHM6c6xzGB4+nITQBP+6zymEGPAkGfoJTxLcto32igqCY2Ox3n474ePG8a2gIMxDhxqrSQsKGJySwqrVfdAfr7XG6BNY+zklpafJHjq8y+HsoeGUVFRQ7qjhaL1x76+irQKAyOBIRoWPMopehw8jLCjMu7EKIcQVkGToY1prWo8cMZJgeXmXJNi552BeXp73k592GQWyaz83kmBTqTEeNohBseEUFDd7ZoYABcXNRMeZyD+VD0CaOY0ZMTMYET6CJFOSzP6EEAFDkqGPaK1pPXrUuBxaXk6w1Yr1ttuMJNiXHSacLUa3+NrP4exBcNQDQUZJtCF3Quw4nOYEmribX/2likfvTCR7qJEYf/WXKhxBodyecDsZYRlEBPvhXkchhOgBSYZ9TGtN67FjRhIsKyPYaiXm1luJGD++75Jgyxl38vsc7IVGo9zgCLCOhdhsdMxoamilqLmIovoPKa4qxn7Wzuj7ruaX75yg9lQ51uRwUu4cweHnv2RM5Ji+iVsIIbxEkmEf8STBbdtoLy0lOCaGmPnziZgwwftJULug/vjXCbDZuK9HWDIkz4HYcTRGJHOi9SQnmk9QdDofu9MOgDXEypjIMexJ34M5LoyJv5zmedmawzWkpPfdilYhhPAWSYZeprWm9fhxYybYl0nQ0Qy1B42aoLUHjU3xKgiis2DQdTisozlJK0UtRRQ1fsjps0a397CgMIaGDWVa2DQywjKIDY0FwLTExJpn1jBi2QisI63UHqnl6NqjrF7VdytahRDCWyQZeonWmraiIuq3bqXt5EmCLBZibrnFSIIhXvprbz5tJL+zn0P9EWNGGBIJ1my0NZuqyESOt1VQ1FLEyap1OLSDIII8C18ywjNIMaV02/GhY/FO/rp89pfuJyU9hdWr+mBFqxBC9AGl9Xk1rwNSbm6u3rNnj6/DMJLgiRNGEiwpIchiIfq664jIybmiJHhe9ZclK5gzezbUH/06AbYYszvCU41Ln5YRHAt2UtRSTFFLEU2uJsAoeJ0RlkFGWAZDwobIpnchxICglNqrtc7t7pjMDHtRa1ER9du20VZcTFB0tNFfcOLEK54Jdlv95Y+Pw7F1zBlnNhrjWrJwJF3PyXArR1xnKWoporruMACRQZFkhGd4EqC0OxJCiK4kGV6mzq2UBqeksHj5cq4bPtyYCRYXExQVheWmm4icNKl3Loe62lm/9k/dVH9J4Kl3ahl/+2oKQ4I42l5KWeteXG0uQlQIQ8xDmBA1gYywDBJDE2XPnxBCXIQkw8vQXSulP/32t9QOHcr0zMzeSYLaCQ0lYP/SaIVkP0px+Zluq78Unyrn6ZZdACSbkrnGcg3DwoYxOGywdHsQQojLID8xL0O3rZRGjmRdWRl3/O53qNDQi79Ad7QLmso7Jb9CYyM84AxPwR6fTURYYbfVXyLCglmQsIBhYcNkw7sQQlwBSYaX4eQFWimVFRT0PBFqDS2VRvKr+8pIgI4GANpMVs5EpXPCZOZgcDtnVCtQTWOLk1+9UcWjCztVf3mjiqYWJ6MjR/fyVymEEAOPJMMeaK+spH7rVlK+aSul1hoj6dV9ia77EtVeC0BLSDjlYdEUhkRxxBRMfXAICgeJobGkmpLJMSWRbErm0OBDmMeZ+OU7lZ7qLzHTUjF/3u69L1oIIQYQSYZu3S2MmZmTQ/3WrTQXFKBMJr4zbx7PbtrEyszMi7dSardD3Ve47F/iqjtESGsNAC1BIZwwhVEUHscJcxj2YDODzEbCu86UTJIpiUGmQefd71tx7wr3hverZMO7EEJ4gSRDLrAw5ne/o3bIEKalphI1bRpR06dzS0QEW0+W8MTB/TS3ugg3B3HV2AnkzbwWZ/Vemmv3E2w/Srg7+bUrRbEpjGJLLKXmKIIjBpNsTiHNlMwkUzLxofEEq0tXoZEN70II4V2y6R64b8kS7klI6HL581BtLetOnuTZ/HyCo6IAePLJJ9hu+xu/+E6S597dL189zfVjI3nw1kG0ozhpMlMWFkVjVDqmqBEkmVNJMiURFxLXbWUXIYQQfaNfbLpXSt0EPAEEA89orX/dW699sqKCXdrFmsKDNLa4iAwLYkrsIMpqa2kJacB+eidtDcfY9MFf+dd7krvs9/vFd5L4l5dOM+MHeURYRpMUlk5GSIzs6xNCiAASEMlQKRUM/A9wA1AKfKqUeltrfag3Xj8yIpzP6qt4bHFylxlfhFkRuf9fiHSf19TsIntoeJfnZg8Np6nZSU7at3sjFCGEED4QKNftpgBHtdbHtdZtwCvA7b314i5XK7/4ThI5mRGEBCvPjE8DxUlTODX8HppzfklkRAgFxc1dnltQ3ExkRED8TiGEEOICAiUZpgEnOz0udY95VFdXk5eX5/lYu3Ztj1+8qdnZ/YyvxcXQjJUkJ84k3JxA3g3f4pevnuaz4004nJrPjjfxy1dPk3fDt67gSxNCCOFrgTKl6e4GXJeVP/Hx8dhstm/04h0zvnMrvJw74/vRjx4E4LGX36GxyUFkRAh5N9ziGRdCCBGYAiUZlgKDOz1OB8p768WNGd/5q0TzbrjlvHN/9KMHJfkJIUQ/EyjJ8FNgpFIqAygDFgH39NaLy4xPCCEGtoDZZ6iU+hbwO4ytFc9prX/Z+bi/NPcVQgjhn/rFPkOt9TvAO76OQwghRP8TKKtJe93lrDb1J4EYdyDGDIEZdyDGDIEZdyDGDIEZd1/EPGCT4QsvvODrEL6RQIw7EGOGwIw7EGOGwIw7EGOGwIy7L2IesMlQCCGE6CDJUAghxIAXMAtoLmXv3r1nlFLFPT3fbDanKqV6ba9iXwnEuAMxZgjMuAMxZgjMuAMxZgjMuHsx5qEXOhAwWyuEEEIIb5HLpEIIIQY8SYZCCCEGPEmGQgghBrwBmQyVUjcppb5SSh1VSj3i41gGK6W2KKUOK6W+UEo96B6PU0ptVEodcf8Z2+k5j7pj/0opdWOn8UlKqQL3sSeVUt11++jN2IOVUp8ppf4aQDFblVKvKaW+dP+dX+vvcSul/tH9vXFQKfWyUirMH2NWSj2nlKpUSh3sNNZrcSqlzEqpV93ju5RSw7wU83+5vz8+V0q9oZSy+lPMF4q707GfKKW0UirBn+K+UMxKqdXuuL5QSv2nz2LWWg+oD4zapseATMAEHABG+zCeFGCi+/NooBAYDfwn8Ih7/BHgP9yfj3bHbAYy3F9LsPvYbuBajJZX7wI3ezn2HwMvAX91Pw6EmNcCq9yfmwCrP8eN0bezCAh3P94ALPfHmIEZwETgYKexXosT+AHwR/fni4BXvRTzPCDE/fl/+FvMF4rbPT4YeB8oBhL8Ke4L/F3PBjYBZvfjQb6K2Ws/dPz1w/2X+H6nx48Cj/o6rk7xvAXcAHwFpLjHUoCvuovX/Y1/rfucLzuNfxd42otxpgM2YA5fJ0N/j9mCkVjUOeN+GzdfN7aOw9gK9VeMH9Z+GTMw7Jwfdr0WZ8c57s9DgDPn/lv2RsznHFsIrPe3mC8UN/AaMB44wdfJ0G/i7ub7YwMwt5vz+jzmgXiZtOOHS4dS95jPuaf1OcAuIElrXQHg/nOQ+7QLxZ/m/vzccW/5HfAzwNVpzN9jzgSqgOeVcXn3GaVUpD/HrbUuA/4bKAEqgDqt9Qf+HPM5ejNOz3O01g6gDoj3WuSG+zBmH13e/5zY/CJmpdRtQJnW+sA5h/w57izgevdlzW1Kqcm+inkgJsPu7pP4fLOlUioK+AvwD1pr+8VO7WZMX2S81yml5gOVWuu9PX1KN2N9GrNbCMZlmj9orXOARoxLdxfi87jd99hux7hUlApEKqXuvdhTLhCbv33ff5M4+/RrUEr9AnAA6y/x/j6PWSkVAfwC+OfuDl8gBp/HjfF/MhaYCvwU2OC+B9jnMQ/EZFiKcV29Qzrg02oMSqlQjES4Xmv9unv4tFIqxX08Bah0j18o/lL35+eOe8N04Dal1AngFWCOUmqdn8fcEUep1nqX+/FrGMnRn+OeCxRprau01u3A68A0P4+5s96M0/McpVQIEAPUeCNopdQyYD6wWLuvu/l5zMMxfmE64P5/mQ7sU0ol+3ncpcDr2rAb40pTgi9iHojJ8FNgpFIqQyllwrjR+ravgnH/FvQscFhr/ZtOh94Glrk/X4ZxL7FjfJF75VQGMBLY7b4EVa+Umup+zaWdntOrtNaPaq3TtdbDMP7+Nmut7/XnmN1xnwJOKqVGuYfygEN+HncJMFUpFeF+rzzgsJ/H3Flvxtn5te7C+L7zxizrJuBh4DatddM5X4tfxqy1LtBaD9JaD3P/vyzFWJh3yp/jBt7EWHeAUioLY1HbGZ/E3Bs3cgPtA/gWxqrNY8AvfBzLdRhT+c+B/e6Pb2Fc67YBR9x/xnV6zi/csX9FpxWBQC5w0H3sKXrpRv0l4p/F1wto/D5mYAKwx/33/SbGJRq/jhv4V+BL9/u9iLHCzu9iBl7GuK/ZjvHDeGVvxgmEAX8GjmKsKMz0UsxHMe49dfx//KM/xXyhuM85fgL3Ahp/ifsCf9cmYJ07hn3AHF/FLLVJhRBCDHgD8TKpEEII0YUkQyGEEAOeJEMhhBADniRDIYQQA54kQyGEEAOeJEMh/JRS6ldKqVlKqQXqEt1V3OeM7qvYhOhvJBkK4b+uwahTOxP48BLnLsCo9C+E+AZkn6EQfkYp9V/AjXzdumY4RreN1zA2LX8PY7PyUWAJRiGBv2IUJq4D7sRoB/ZHIML9Gvdprc8qpbYCnwGTgESMCh6PAtkYLW/+P3fx8g0Ypa6CgX/TWr/q7a9bCF+SZCiEH1JKTcFIdD8Gtmqtp7vH47XW1e7P/x04rbVeo5TKx6gE9Jr72OfAaq31NqXU/wEsWut/cCfDXVrrh5XRSPphjMRYg5E0x2NUFbpJa/137teK0VrX9dXXLoQvyGVSIfxTDkYpsKsw6qd2GKuU+lApVQAsBsac+0SlVAxg1Vpvcw+txWis2qGjFm8B8IXWukJr3Qocxyh0XADMVUr9h1LqekmEYiAI8XUAQoivKaUmAPkYlyjPYFzmVEqp/RjNTfOBBVrrA0qp5RizuMvV6v7T1enzjschWutCpdQkjBq5v1JKfaC1/j/f4H2ECBgyMxTCj2it92utJ2AUkh8NbAZu1FpP0Fo3Y9wLrHC3/Vrc6an17mO4Z3JnlVLXu48tAbbRQ0qpVKBJa70Oo7HwxCv7qoTwfzIzFMLPKKUSgbNaa5dS6iqtdefLpP+EscK0GONyZrR7/BXgT0qpH2G0r1kG/NHd9PU4sOIyQsgG/ksp5cLoMPD9K/qChAgAsoBGCCHEgCeXSYUQQgx4kgyFEEIMeJIMhRBCDHiSDIUQQgx4kgyFEEIMeJIMhRBCDHiSDIUQQgx4/z9B+OOm3LJgnAAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 508.095x314.02 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "mlt = Multiplot()\n",
    "mlt[0,0].plot(natoms, t_voroc, '-', color=accent['pred'])\n",
    "mlt[0,0].plot(natoms, t_voro, '-', color=accent['pgreen'])\n",
    "mlt[0,0].plot(natoms, t_scipy, '-', color=accent['pyellow'])\n",
    "mlt[0,0].plot(natoms, t_voroc, 'o', label=\"pyscal voro+clean\", color=accent['pred'], markeredgecolor=accent[\"dgrey\"])\n",
    "mlt[0,0].plot(natoms, t_voro, 'o', label=\"pyscal voro\", color=accent['pgreen'], markeredgecolor=accent[\"dgrey\"])\n",
    "mlt[0,0].plot(natoms, t_scipy, 'o', label=\"scipy\", color=accent['pyellow'], markeredgecolor=accent[\"dgrey\"])\n",
    "mlt[0,0].set_xlabel(\"#atoms\")\n",
    "mlt[0,0].set_ylabel(\"Time (ms)\")\n",
    "mlt[0,0].legend()\n",
    "plt.savefig(\"voro_timing.png\", dpi=300, bbox_inches=\"tight\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "15379cfe-daf2-4092-b8ed-c643e8698d14",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "56b23cd3-ce49-4b9c-a66f-5fb928f5e008",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "752e1f62-e691-4bde-9dac-71e2b49efec0",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "8c20c942-8023-4425-b507-6e6a65e8e425",
   "metadata": {},
   "outputs": [],
   "source": [
    "def get_distance(x, box, dim=3):\n",
    "    box_2 = box/2\n",
    "    xdum = []\n",
    "    for i in range(dim):\n",
    "        xd = np.meshgrid(x[i], x[i])[1] - np.meshgrid(x[i], x[i])[0]\n",
    "        xdum.append(xd)\n",
    "    xd = np.array(xdum)\n",
    "    for i in range(dim):\n",
    "        xd[i] = np.where(xd[i] > box_2[i], xd[i]-box[i], xd[i])\n",
    "        xd[i] = np.where(xd[i] < -box_2[i], xd[i]+box[i], xd[i])\n",
    "    return xd"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "73b78891-ef91-45f0-aba1-1e8be0fc187f",
   "metadata": {},
   "outputs": [],
   "source": [
    "x = sys.atoms['positions']\n",
    "box = sys.boxdims"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "id": "6337c594-f8ca-447d-8277-f86987ce3e21",
   "metadata": {},
   "outputs": [],
   "source": [
    "x = np.array(x)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "id": "ed395d26-9344-4ca1-acd4-5ae7a8949748",
   "metadata": {},
   "outputs": [],
   "source": [
    "y = x.T"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "id": "9d26f239-c648-4454-b45d-4fbede636145",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "224 ms ± 857 µs per loop (mean ± std. dev. of 7 runs, 1 loop each)\n"
     ]
    }
   ],
   "source": [
    "%%timeit\n",
    "xd = get_distance(y, np.array(box))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "id": "0eb9ab1c-921e-48a2-a753-8eb9d04371c2",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(3, 2000, 2000)"
      ]
     },
     "execution_count": 27,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "xd.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "id": "5c35e8a0-4e7c-48fc-a7d2-e7a348d2896c",
   "metadata": {
    "collapsed": true,
    "jupyter": {
     "outputs_hidden": true
    },
    "tags": []
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[[ 0.    , -1.5635,  0.    , ...,  1.5635,  3.127 ,  1.5635],\n",
       "        [ 1.5635,  0.    ,  1.5635, ...,  3.127 ,  4.6905,  3.127 ],\n",
       "        [ 0.    , -1.5635,  0.    , ...,  1.5635,  3.127 ,  1.5635],\n",
       "        ...,\n",
       "        [-1.5635, -3.127 , -1.5635, ...,  0.    ,  1.5635,  0.    ],\n",
       "        [-3.127 , -4.6905, -3.127 , ..., -1.5635,  0.    , -1.5635],\n",
       "        [-1.5635, -3.127 , -1.5635, ...,  0.    ,  1.5635,  0.    ]],\n",
       "\n",
       "       [[ 0.    , -1.5635,  0.    , ...,  1.5635,  3.127 ,  1.5635],\n",
       "        [ 1.5635,  0.    ,  1.5635, ...,  3.127 ,  4.6905,  3.127 ],\n",
       "        [ 0.    , -1.5635,  0.    , ...,  1.5635,  3.127 ,  1.5635],\n",
       "        ...,\n",
       "        [-1.5635, -3.127 , -1.5635, ...,  0.    ,  1.5635,  0.    ],\n",
       "        [-3.127 , -4.6905, -3.127 , ..., -1.5635,  0.    , -1.5635],\n",
       "        [-1.5635, -3.127 , -1.5635, ...,  0.    ,  1.5635,  0.    ]],\n",
       "\n",
       "       [[ 0.    , -1.5635, -3.127 , ...,  4.6905,  3.127 ,  1.5635],\n",
       "        [ 1.5635,  0.    , -1.5635, ...,  6.254 ,  4.6905,  3.127 ],\n",
       "        [ 3.127 ,  1.5635,  0.    , ...,  7.8175,  6.254 ,  4.6905],\n",
       "        ...,\n",
       "        [-4.6905, -6.254 , -7.8175, ...,  0.    , -1.5635, -3.127 ],\n",
       "        [-3.127 , -4.6905, -6.254 , ...,  1.5635,  0.    , -1.5635],\n",
       "        [-1.5635, -3.127 , -4.6905, ...,  3.127 ,  1.5635,  0.    ]]])"
      ]
     },
     "execution_count": 23,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.where(xd[]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7669c585-0a71-47dd-9f02-a2d8c50f526f",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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.10.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}