optimize-jacobian-variance.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8903231c-3c64-41d7-8b84-c8ae23c24bb5",
   "metadata": {},
   "outputs": [],
   "source": [
    "%load_ext autoreload\n",
    "%autoreload 2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1dbab2e3-e238-414a-8383-c4e65d034052",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import torch\n",
    "import tqdm\n",
    "\n",
    "from utils import rotation_matrix\n",
    "#from methods import *"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b6d6a639-767f-42d4-8a7b-0eb74c72c6c4",
   "metadata": {},
   "outputs": [],
   "source": [
    "from models import Similitude, Similitude2, RobustSimilitude\n",
    "dtype=torch.float32\n",
    "model = Similitude2(scale=None,ncells=14,dtype=dtype,rectify=True,dropout=True) #146.2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5b5f4247-d3b2-48b5-add7-e6b18f00cca3",
   "metadata": {},
   "outputs": [],
   "source": [
    "if True:\n",
    "    # Overwrite Initial Phase Placements\n",
    "    start_phases = torch.normal(torch.zeros((model.ncells,2))+0.2,0.01)\n",
    "    #start_phases = torch.cat([torch.normal(torch.zeros((model.ncells//2,2)),0.01),torch.normal(torch.zeros((model.ncells//2,2))+0.2,0.01)])\n",
    "    \n",
    "    #model = torch.load('./saved-models/Similitude2-14-dirac')\n",
    "    #start_phases = model.phases.detach().numpy()\n",
    "    #start_phases = np.concatenate([start_phases,start_phases@rotation_matrix(10).T+0.2])\n",
    "    #start_phases = np.concatenate([start_phases,start_phases+0.2])\n",
    "    \n",
    "    model.set_phases(start_phases)\n",
    "start_phases = model.phases.detach().numpy()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "13128e90-d28f-4648-95c4-6d7dfb801851",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.collections.PathCollection at 0x7f8b603371c0>"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAi8AAAGdCAYAAADaPpOnAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy88F64QAAAACXBIWXMAAA9hAAAPYQGoP6dpAABVk0lEQVR4nO3dZ3xUZd7G8d9MOpAChDQIIB2kKUioNqJgxYaoWEDsYkNdCb3HgsqqKIr1WXVRsLGoQRZEBUIRpDfpLYUQSEJ6Zs7zIho3Ckggk3tmcn0/n/NiD+cwV25Y52L+55yxWZZlISIiIuIh7KYDiIiIiFSEyouIiIh4FJUXERER8SgqLyIiIuJRVF5ERETEo6i8iIiIiEdReRERERGPovIiIiIiHsXXdIDK5nQ6OXToEMHBwdhsNtNxRERE5DRYlkVOTg4xMTHY7af+bMXrysuhQ4eIjY01HUNERETOwP79+2nQoMEpj/G68hIcHAyU/vAhISGG04iIiMjpyM7OJjY2tux9/FS8rrz8PioKCQlReREREfEwp3PJhy7YFREREY+i8iIiIiIeReVFREREPEqVlJfp06fTuHFjAgMDiYuLY+XKlac8/tixYzz88MNER0cTEBBAixYt+Oabb6oiqoiIiLg5l1+w+8knnzBs2DBmzJhBXFwc06ZNo0+fPmzbto2IiIi/HF9UVMRll11GREQEc+bMoX79+uzdu5ewsDBXRxUREREPYLMsy3LlC8TFxXHBBRfw2muvAaUPkYuNjeWRRx5h+PDhfzl+xowZvPDCC2zduhU/P78Kv152djahoaFkZWXpbiMREREPUZH3b5eOjYqKili9ejXx8fF/vKDdTnx8PMnJySc8Z+7cuXTr1o2HH36YyMhI2rZty5QpU3A4HCc8vrCwkOzs7HKbiIiIeC+XlpeMjAwcDgeRkZHl9kdGRpKamnrCc3bt2sWcOXNwOBx88803jB49mhdffJFJkyad8PjExERCQ0PLNj1dV0RExLu53d1GTqeTiIgI3nrrLTp16sSAAQMYOXIkM2bMOOHxCQkJZGVllW379++v4sQiIiJSlVx6wW54eDg+Pj6kpaWV25+WlkZUVNQJz4mOjsbPzw8fH5+yfa1btyY1NZWioiL8/f3LHR8QEEBAQEDlhxcRERG35NJPXvz9/enUqRMLFy4s2+d0Olm4cCHdunU74Tk9evRgx44dOJ3Osn3bt28nOjr6L8VFREREqh+X3yo9bNgw7rrrLjp37kyXLl2YNm0aubm5DB48GIA777yT+vXrk5iYCMCDDz7Ia6+9xmOPPcYjjzzCr7/+ypQpU3j00UddHVW8nWVBcZ7pFCIC4FcDTuM7bEROxOXlZcCAARw+fJgxY8aQmppKx44dSUpKKruId9++fdjtf3wAFBsby/z583niiSdo37499evX57HHHuOZZ55xdVTxckePH+Lpjy7m8cxjtC0qMh1HpFqzEg5iC6hlOoZ4KJc/56Wq6TkvcjL/+H4Y3+5bQIvCImYfSnW/q9VFqpGvr1nNVZ2amY4hbqQi798u/+RFxB0s2LuAb/ctwMfmw/hrPsRet43pSCLV0mvf72D69zsImLeTC1o0ICI40HQk8UAqL+L1MgsymbS89DlBd7e9m7bRFxhOJFJ93R/fjm+2ZrM5JZuRX2zkrTs6YdO1L1JB+uRcvN7k5ZPJLMikee3mPNDhAdNxRKo1Px87L97cAT8fGws2p/HV2kOmI4kHUnkRr5a0J4nv9n6Hr82XST0m4e+j2+1FTGsdHcKjlzYHYOzcTaRlFxhOJJ5G5UW8VkZ+BpOXTwbgnvb30EbXuYi4jQcubkq7+qFk5Rcz4vMNeNm9I+JiKi/ilSzLYtLySRwrPEbL2i25r919piOJyP/w87EztX8H/H3sLNyazmdrDpqOJB5E5UW80re7v2XhvoX42nyZ3HMyfj5+piOJyJ+0jArm8ctKx0fj/7OJ1CyNj+T0qLyI1zmcd5jJK0rHRfd3uJ+WdVoaTiQiJ3NfryZ0iA0jp6CE4Z+v1/hITovKi3gVy7KYkDyB7KJsWtdpzZB2Q0xHEpFT8PWx82L/9vj72lm87TCzfz5gOpJ4AJUX8Srzds1j8YHF+Np9mdRzEn52jYtE3F2ziGCevKwFABPnbebQsXzDicTdqbyI10jPSydxZekXfD7U4SFa1G5hOJGInK57ejXhvIZh5BSW8MxnGh/Jqam8iFewLIvxyePJKcrh3LrnMrjtYNORRKQCfOw2pvbvQICvnZ9+zWDWqv2mI4kbU3kRr/DVzq/48cCP+Nn9mNxzMr52ffOFiKdpWq8WT/cpvcB+0rzNHDiaZziRuCuVF/F4qbmpPLfyOQAe7vgwTcOaGk4kImdqcI9z6NyoNrlFDo2P5KRUXsSjWZbFuORxHC8+Tvvw9tx17l2mI4nIWfCx23ihfwcC/ews3XGEj1bsMx1J3JDKi3i0L3Z8wdKDS/G3+zOx50SNi0S8wDnhNXmmbysApnyzhf2ZGh9JeSov4rFSjqfw/KrnAXj0/EdpEtrEcCIRqSx3dWtMl3PqkFfk4Ok563A6NT6SP6i8iEeyLIsxy8aQW5xLx3odub317aYjiUglstttTL2pA0F+Pizflcm/lu81HUnciMqLeKTZ22ezPGU5AT4BTOwxER+7j+lIIlLJGtatQcKVpeOjZ7/dyt4juYYTibtQeRGPc/D4QV78+UUAHjv/MRqHNjYbSERc5va4RnRrUpf8YgdPz16v8ZEAKi/iYZyWk7FLx5JXksf5EeczsPVA05FExIXsdhvP39Semv4+rNyTyfvL9piOJG5A5UU8yqfbPmVF6gqCfIOY2GMidpv+Cot4u9g6NUi4sjUAz8/fyq7Dxw0nEtP0X37xGPtz9vPS6peA0nFRw5CGhhOJSFUZGNeQns3CKSh28vSc9Tg0PqrWVF7EIzgtJ2OWjiG/JJ/OkZ25tdWtpiOJSBWy2Ww8e2M7agX4snrvUd5butt0JDFI5UU8wr+3/puf034myDeICT0maFwkUg01qF2DUVeVjo9emL+NHekaH1VXegcQt7cvex/TVk8D4MlOTxIbHGs2kIgYM+CCWC5sUY/CEidPzV6n8VE1pfIibs3hdDBq6SgKHAXERcXRv2V/05FExCCbzcazN7QjOMCXtfuPMfOnXaYjiQEqL+LWPtryEb+k/0IN3xqM7zFe4yIRISYsiNHXtAHgpe+282tajuFEUtX0TiBua0/WHl755RUAnrrgKerXqm84kYi4i/6dGnBJy3oUOUrHRyUOp+lIUoVUXsQt/T4uKnQU0i26Gzc1v8l0JBFxIzabjcQb2hMS6Mu6A1m8+aPGR9WJyou4pX9t/hfrDq+jll8tJvSYgM1mMx1JRNxMVGggY685F4Bp/93O1tRsw4mkqqi8iNvZdWwXr/7yKgBPX/A0UTWjDCcSEXd1w/n1iW8dQbHD4qnZ6yjW+KhaUHkRt1LiLGHU0lEUOYvoUb8H1ze73nQkEXFjNpuNKde3IzTIj40Hs5mxeKfpSFIFVF7ErXyw6QM2ZGwg2C+Ycd3GaVwkIn8rIiSQCf1Kx0evLPqVzYc0PvJ2Ki/iNnYc3cH0tdMBeKbLMxoXichpu7ZDDJe3iSwbHxWVaHzkzVRexC38Pi4qdhZzYYMLubbptaYjiYgHsdlsTL6+HbVr+LE5JZvp3+8wHUlcSOVF3MJ7G99j05FNBPsHM7bbWI2LRKTC6gUHMKFfWwCmf7+DjQezDCcSV1F5EeO2H93O6+teByChSwIRNSIMJxIRT3V1+2iubBdFiVPjI2+m8iJGFTuLGbVkFCXOEi6JvYSrm1xtOpKIeDCbzcbEfm2pW9Ofrak5vLroV9ORxAVUXsSotze8zZbMLYQGhDKm2xiNi0TkrNWtFcDE60rHR68v3sn6A8fMBpJKp/IixmzN3Mpb694CYESXEYQHhRtOJCLe4sp20VzdPhrHb+OjwhKH6UhSiVRexIhix2/jIquE+IbxXHHOFaYjiYiXmdCvLeG1/Nmedpx//lfjI2+i8iJGvLXhLbYd3UbtgNqM6jpK4yIRqXR1avoz6bp2AMz4YSdr9x8zG0gqjcqLVLnNRzYzc/1MAEZ2HUndoLqGE4mIt+rbNop+HWNwWvDkp2spKNb4yBuovEiVKnIUMXLJSByWg8sbXU6fxn1MRxIRLzfumnOpFxzAzsO5vLxgu+k4UglUXqRKzVg3gx3HdlAnsA4ju440HUdEqoHaNf2Zcn3p+GjmT7tYvfeo4URytlRepMpszNjIuxvfBWBU11HUCaxjOJGIVBeXtYnkhvPr47Tg6dnrND7ycFVSXqZPn07jxo0JDAwkLi6OlStXntZ5s2bNwmazcd1117k2oLhcoaOQUUtG4bAcXHHOFVzW6DLTkUSkmhl79blEhgSwKyOXqfO3mY4jZ8Hl5eWTTz5h2LBhjB07ljVr1tChQwf69OlDenr6Kc/bs2cPTz31FL169XJ1RKkCr699nZ1ZO6kbWJcRXUaYjiMi1VBoDT8SbygdH72zdDer9mQaTiRnyuXl5aWXXuLee+9l8ODBtGnThhkzZlCjRg3efffdk57jcDgYOHAg48ePp0mTJq6OKC627vA63t/0PgBjuo0hLDDMaB4Rqb4ubRVJ/04NsH4bH+UXaXzkiVxaXoqKili9ejXx8fF/vKDdTnx8PMnJySc9b8KECURERDBkyJC/fY3CwkKys7PLbeI+CkoKGL10NE7LydVNrubShpeajiQi1dyoq9sQHRrIniN5PD9/q+k4cgZcWl4yMjJwOBxERkaW2x8ZGUlqauoJz1myZAnvvPMOM2fOPK3XSExMJDQ0tGyLjY0969xSeaavnc7urN3UC6rH8C7DTccRESE0yI9nb2wPwHtL97B81xHDiaSi3Opuo5ycHO644w5mzpxJePjpfc9NQkICWVlZZdv+/ftdnFJO19r0tXyw6QMAxnYbS2hAqOFEIiKlLmpRj1suKP3H7j/mrCe3sMRwIqkIX1f+5uHh4fj4+JCWllZuf1paGlFRUX85fufOnezZs4drrrmmbJ/T6SwN6uvLtm3baNq0ablzAgICCAgIcEF6ORv5JfmMWjoKC4trm17LRbEXmY4kIlLOyKta8+P2w+zLzOO5pK1M6NfWdCQ5TS795MXf359OnTqxcOHCsn1Op5OFCxfSrVu3vxzfqlUrNmzYwNq1a8u2a6+9lksuuYS1a9dqJORBXv3lVfZm7yUiKIJnujxjOo6IyF8EB/rx/E0dAPi/5L0s25lhOJGcLpd+8gIwbNgw7rrrLjp37kyXLl2YNm0aubm5DB48GIA777yT+vXrk5iYSGBgIG3blm++YWFhAH/ZL+5rddpqPtz8IQDjuo8jxD/EcCIRkRPr2TycgXEN+WjFPv4xZz1Jj19IrQCXvzXKWXL5n9CAAQM4fPgwY8aMITU1lY4dO5KUlFR2Ee++ffuw293q0hs5C3nFeYxeOhoLixua30CvBnpOj4i4t4QrW7N422EOHM0n8ZstTP7tqwTEfdksy7JMh6hM2dnZhIaGkpWVRUiI/sVf1RJXJPLx1o+JqhnF59d+TrB/sOlIIiJ/a9mODG57ewUAHw6Jo2fz07tpRCpPRd6/9ZGHVJpVqav4eOvHAIzvNl7FRUQ8Rvdm4dzZrREAz3y2npyCYsOJ5FRUXqRS/D4uAripxU10r9/dcCIRkYp5pm8rGtapwcFj+Uz5ZovpOHIKKi9SKV5a/RIHjx8kpmYMT3V+ynQcEZEKqxngyws3lT687t8r9/PD9sOGE8nJqLzIWVuespxPtn0CwPge46npV9NwIhGRMxPXpC6DujcG4Jk568nK1/jIHam8yFnJLc5l7NKxAAxoOYCu0V0NJxIROTv/6NuSxnVrkJpdwKR5m03HkRNQeZGz8uLPL3Io9xD1a9VnWKdhpuOIiJy1Gv6+vNC/AzYbzF59gO+3ppuOJH+i8iJnbNnBZczePhuAiT0mUsOvhuFEIiKV44LGdRjS4xwAhn++nqw8jY/cicqLnJGcohzGJpeOi25rdRsXRF1gOJGISOV6qk9LmoTXJC27kPHzNpmOI/9D5UXOyNSfp5Kam0pscCyPnf+Y6TgiIpUu0M+HF/p3wG6Dz9cc5L+b0/7+JKkSKi9SYUsOLuHzXz/Hhk3jIhHxap0a1ebeXk0ASPhiA8fyigwnElB5kQrKLspm7LLScdHA1gPpFNnJcCIREdd64rIWNK1Xk8M5hYybq/GRO1B5kQp5fuXzpOel0yikEY+e/6jpOCIiLhfo58OLN3fEboMv1x4iaWOq6UjVnsqLnLYf9v/AVzu/woaNST0mEeQbZDqSiEiV6Bgbxv0XNQVg1JcbyMzV+MgklRc5LVmFWYxPHg/AnW3upGNER7OBRESq2OPxzWkRWYuM40WM+Wqj6TjVmsqLnJbnVj7H4fzDNA5pzNDzhpqOIyJS5QJ8fZjavwM+dhvz1qfwzYYU05GqLZUX+VuL9i3iP7v+g91mZ1LPSQT6BpqOJCJiRPsGYTx08e/jo41kHC80nKh6UnmRUzpWcIwJyRMAGHTuIDrU62A4kYiIWY9c2pxWUcFk5hYx+suNWJZlOlK1o/IipzRl5RSOFByhaWhTHur4kOk4IiLG+fvamdq/A752G99uTGXeeo2PqprKi5zUf/f+l293f4uPzYdJPScR4BNgOpKIiFtoWz+Uhy9pBsCYrzZyOEfjo6qk8iInlFmQycTlEwG4u+3dtA1vaziRiIh7efiSZrSJDuFoXjGjvtyg8VEVUnmRE5qyYgqZBZk0C2vGAx0eMB1HRMTt/O/4aP6mNOauO2Q6UrWh8iJ/MX/PfObvmY+PzYfJPSfj7+NvOpKIiFtqExPCo72bAzDmq02kZxcYTlQ9qLxIOUfyjzB5+WQA7ml3D23qtjGcSETEvT14cVPa1g8hK7+YEV9ofFQVVF6kjGVZTF4xmaOFR2lRuwX3t7/fdCQREbfn52Pnxf4d8fOx8d8t6Xzxy0HTkbyeyouUSdqTxIK9C/C1+TK552T8fPxMRxIR8Qgto4J5PL4FAOPmbiI1S+MjV1J5EQAy8jOYvKJ0XHRfh/toVaeV4UQiIp7l/gub0L5BKNkFJSR8vl7jIxdSeREsy2JC8gSyCrNoXac197S7x3QkERGP4+tj58X+HfD3sfP9tsPMXn3AdCSvpfIifL37a77f/z2+dl8m9ZyEn13jIhGRM9E8Mphhl5eOjyb+ZzMpWfmGE3knlZdqLj0vncQViQA82OFBWtRuYTiRiIhnu7dXE85rGEZOYQnPfKa7j1xB5aUa+31clF2UTZu6bbi77d2mI4mIeDwfu40XbuqAv6+dH7cf5pNV+01H8joqL9XY3J1z+eHAD/jZ/ZjcYzK+dl/TkUREvEKziFo8fXlLACZ9vYWDxzQ+qkwqL9VUWm4az618DoCHOj5Es9rNDCcSEfEud/c8h06NanO8sIRn5ujuo8qk8lINWZbFuORx5BTn0C68HYPOHWQ6koiI1ykdH7Un0M/Okh0ZfLxyn+lIXkPlpRr6cseXLDm4BH+7P5N6TNK4SETERZrUq8XTfUqfmzX56y3sz8wznMg7qLxUMynHU3h+1fMAPHLeIzQJa2I4kYiIdxvcvTFdGtchr8jBP+asx+nU+OhsqbxUI5ZlMXbZWI4XH6dDvQ7c0eYO05FERLye3W7j+ZvaE+TnQ/KuI3y4Yq/pSB5P5aUa+ezXz0hOSSbAJ4BJPSbhY/cxHUlEpFpoHF6T4VeUjo8Sv9nKviMaH50NlZdq4tDxQ7yw6gUAHj3vURqHNjYbSESkmrmjayO6NqlDfrGDp+as0/joLKi8VANOy8mYpWPIK8nj/IjzGdh6oOlIIiLVjv23h9fV8Pdh5e5MPkjeYzqSx1J5qQZmb5vNitQVBPoEMrHHRI2LREQMia1Tg4QrWwPwXNJWdmfkGk7kmVRevNyBnAO8uPpFAB7v9DgNQxoaTiQiUr0N7NKQHs3qUlDs5OnZ63BofFRhKi9ezGk5GbNsDPkl+XSK7MStrW41HUlEpNqz2208d2N7avr78PPeo7y3dLfpSB5H5cWLzdo6i1WpqwjyDWJij4nYbfrjFhFxBw1q12DkVW0AeGH+NnYePm44kWfRu5mX2p+9n2lrpgEwrNMwYoNjzQYSEZFybu0SS6/m4RSWaHxUUSovXshpORm1dBT5JfnERcVxc8ubTUcSEZE/sdlKx0fBAb6s2XeMd5bsMh3JY6i8eKGPt3zMmvQ11PCtwfge4zUuEhFxUzFhQYy+unR8NPW77exIzzGcyDPoXc3L7Mnawz/X/BOAJzs/Sf1a9Q0nEhGRU+nfuQEXt6xHUYmTJ2evp8ThNB3J7VVJeZk+fTqNGzcmMDCQuLg4Vq5cedJjZ86cSa9evahduza1a9cmPj7+lMfLHxxOB6OXjqbAUUDX6K70b9HfdCQREfkbNpuNxBvaERzoy7r9x3jrJ42P/o7Ly8snn3zCsGHDGDt2LGvWrKFDhw706dOH9PT0Ex6/ePFibr31Vr7//nuSk5OJjY3l8ssv5+DBg66O6vE+3PIhaw+vpaZfTSZ0n4DNZjMdSURETkN0aBBjrzkXgGkLfmVbqsZHp2KzLMullzfHxcVxwQUX8NprrwHgdDqJjY3lkUceYfjw4X97vsPhoHbt2rz22mvceeedf3t8dnY2oaGhZGVlERISctb5PcWurF3c/J+bKXQUMq7bOG5scaPpSCJShT5cvpep87cB8FSfltzetZHhRFJRlmVxzwc/s3BrOu3qh/L5Q93x86k+V3dU5P3bpatSVFTE6tWriY+P/+MF7Xbi4+NJTk4+rd8jLy+P4uJi6tSpc8JfLywsJDs7u9xW3TicDkYvGU2ho5AeMT24ofkNpiOJSBV7Y/FOjuUXcyy/mDcW7zQdR86AzWZjyg3tCA3yY8PBLN78QX+OJ+PS8pKRkYHD4SAyMrLc/sjISFJTU0/r93jmmWeIiYkpV4D+V2JiIqGhoWVbbGz1e57JB5s/YH3GeoL9ghnXfZzGRSJe6sPle+nx7CI+XL73L/+7U6Pa2IAgPzsPXtzUbFA5Y5EhgYy7tvTuo38u/JUtKdXvH+Snw60/j3r22WeZNWsWX3zxBYGBgSc8JiEhgaysrLJt//79VZzSrJ3HdvLaL6UjuX90+QdRNaMMJxIRV/hw+V7GfrWRg8fyGfvVxrIx0cFj+Uz+egtfrz+EBdSpGQBQruSIZ7muY30uaxNJscPiqdnrKNbdR3/h0vISHh6Oj48PaWlp5fanpaURFXXqN9mpU6fy7LPP8t1339G+ffuTHhcQEEBISEi5rboocZYwaskoip3FXNjgQvo17Wc6koi4wO/FxfHbFYoOi3KjoYJiBw4LfGzw4MVNeWPxTg4ey9f4yEPZbDYmX9+WsBp+bDqUzevf68/xz1xaXvz9/enUqRMLFy4s2+d0Olm4cCHdunU76XnPP/88EydOJCkpic6dO7syokd7f9P7bDyykWD/YMZ2G6txkYiXemPxzrJycm2HGOqHBfHgxU25sEU9fGzQrn4o9cOCGN+vLbd3bcSDFzctO0Y8U0RwIBP6tQXg1UW/sulQluFE7sXX1S8wbNgw7rrrLjp37kyXLl2YNm0aubm5DB48GIA777yT+vXrk5iYCMBzzz3HmDFj+Pjjj2ncuHHZtTG1atWiVq1aro7rMbYf3c70tdMBSOiSQESNCMOJRMRVfv805cGLm5a7i+j3UrPpUFZZcQG4vWsj3W3kBa5pH80361NI2pTKk5+uY+7Qnvj7uvXVHlXG5aswYMAApk6dypgxY+jYsSNr164lKSmp7CLeffv2kZKSUnb8G2+8QVFRETfddBPR0dFl29SpU10d1WMUO4sZtWQUJc4SLo69mKubXG06koi40O1dG7F0+KV/KSQPXtwUH9tfx0jiHWw2G5Oub0udmv5sTc3htUW/mo7kNlz+nJeqVh2e8/Lmujd5be1rhAaE8mW/LwkPCjcdSUQM+XD53hN+KnOy/eJ5vl6fwsMfr8HHbuOrh3vQtn6o6Ugu4TbPeZHKty1zGzPWzwBKx0UqLiLV28k+ldFFu97jqvbRXNU+GofT4slP11FY4jAdyTiVFw9S7Chm1NLScVHvhr258pwrTUcSETeli3a9y8R+bQmv5c+2tBxeWajxkcqLB5m5YSZbM7cSFhDGqK6jdHeRiJzUyT6REc9Up6Y/k64rvfvojcU7Wbf/mNlAhqm8eIgtR7Ywc/1MAEZ2HalxkYhINdO3bTTXdojBacGTs9dRUFx9x0cqLx6g2FHMyKUjKbFKuLzR5fRt3Nd0JBERMWD8tecSXiuAHenHmfbf6js+UnnxADPWz+DXo79SJ7AOI7uONB1HREQMqV3TnynXl46P3vpxJ2v2HTWcyAyVFze3KWMT72x4B4BRXUdRJ/DE364tIiLVw+XnRnH9efVxWvBUNR0fqby4sSJHEaOWjsJhObii8RVc1ugy05FERMQNjL2mDRHBAew6nMuL320zHafKqby4sTfWvcGOYzuoG1iXEXEjTMcRERE3EVbDn8Qb2gHw9pLdrN6baThR1VJ5cVMbDm/g3Y3vAjC622jCAsPMBhIREbfSu3UkN3VqgGXBU7PXk19UfcZHKi9uqNBRyMilI3FaTq5qchW9G/Y2HUlERNzQ6KvbEBUSyO6MXF6YX33GRyovbmj6L9PZnbWb8KBwErokmI4jIiJuKjTIj8QbS8dH7y3bzYpdRwwnqhoqL25mbfpa3t/0PgBju40lNMA7v4BLREQqxyUtIxjQORbLgqfnrCevqMR0JJdTeXEjBSUFjF46GguLa5tey8WxF5uOJCIiHmDk1a2JCQ1kX2Yezyd5//hI5cWNvPrLq+zJ3kNEUATPdHnGdBwREfEQIYF+PHdTewDeX7aH5J3ePT5SeXETa9LW8K/N/wJgbPexhPiHGE4kIiKepFfzetzapSEAT89ZR26h946PVF7cQH5Jftm46Ppm13NhgwtNRxIREQ808qrW1A8L4sDRfBK/3WI6jsuovLiBV9a8wr6cfUTWiOTpC542HUdERDxUrQBfnv9tfPTh8n0s3ZFhOJFrqLwYtip1FR9u+RCA8d3HE+wfbDiRiIh4sh7NwrmjayMA/jFnPTkFxYYTVT6VF4PyivMYs3QMADc2v5Ee9XsYTiQiIt5g+BWtiK0TxMFj+Uz5ZqvpOJVO5cWgl1e/zIHjB4iuGc1TnZ8yHUdERLxEzQBfnr+xAwD/XrmPH7cfNpyocqm8GLIyZSWzts0CSsdFtfxrGU4kIiLepFvTugzq3hiA4Z+tJ9uLxkcqLwbkFucyZlnpuGhAywF0i+lmOJGIiHijf/RtSaO6NTiUVcDked5z95HKiwEv/fwSB48fpH6t+gzrNMx0HBER8VI1/H154aYO2Gzwyc/7+X5buulIlULlpYotO7SMT7d/CsCE7hOo4VfDcCIREfFmXc6pw+Du5wCl46OsPM8fH6m8VKHjRccZu2wsALe2upUu0V0MJxIRkerg6T4tOSe8JmnZhUyYt9l0nLOm8lKFpv48ldTcVBrUasDj5z9uOo6IiFQTQf4+TO3fHpsNPltzgIVb0kxHOisqL1Vk6cGlfPbrZwBM7DFR4yIREalSnRrV4d5eTQBI+HwDx/KKDCc6cyovVSC7KLvs7qLbW99O56jOhhOJiEh1NOyyFjSpV5P0nELG/8dzx0cqL1XghVUvkJ6XTsPghjx6/qOm44iISDUV6OfD1P4dsNvgi18O8t2mVNORzojKi4v9eOBHvtzxJTZsTOo5iSDfINORRESkGju/YW3uu7ApACO+2MjRXM8bH6m8uFBWYRbjl40H4M42d3JexHmGE4mIiMDj8c1pHlGLjOOFjJ27yXScClN5caHnVz1Pen46jUMaM/S8oabjiIiIAH+Mj3zsNuauO8S3G1JMR6oQlRcX+X7f98zdORe7zc6knpMI9A00HUlERKRMh9gwHrio9O6jUV9u5MjxQsOJTp/KiwtkFWYxYfkEAO469y461OtgOJGIiMhfPdq7OS0jgzmSW8SYrzxnfKTy4gKJKxPJyM+gSWgTHu74sOk4IiIiJxTg68OLN5eOj77ekMK89YdMRzotKi+VbOHehXy96+vScVGPSQT4BJiOJCIiclJt64fy8CXNABj95UYO57j/+EjlpRIdLThaNi66u+3dtKvXznAiERGRvzf0kma0jg7haF4xo77cgGVZpiOdkspLJZqyYgqZBZk0C2vGgx0eNB1HRETktPj72pnavz2+dhvzN6Uxd517j49UXirJd3u+I2lPEj42Hyb1nIS/j7/pSCIiIqft3JhQHrm0OQBj524iPafAcKKTU3mpBEfyjzBp+SQA7ml3D+fWPddwIhERkYp76JKmnBsTwrG8YkZ+sdFtx0cqL2fJsiwmr5jM0cKjtKjdgvvb3286koiIyBnx87Hz4s0d8POxsWBzGl+uPWg60gmpvJyl+Xvms2DvAnxtvkzqMQk/Hz/TkURERM5Yq6gQHuv92/joq02kZbvf+Ejl5Sxk5GcwecVkAO5rfx+t67Y2nEhEROTsPXBRU9rVDyW7oIQRn7vf3UcqL2fIsiwmLZ/EscJjtKrTinva32M6koiISKXw/W185O9jZ+HWdD5b417jI5WXM/TN7m9YuG8hvvbfxkV2jYtERMR7tIgM5onLWgAw/j+bSMnKN5zoDyovZ+Bw3mGmrJgCwAPtH6BlnZaGE4mIiFS+e3udQ4fYMHIKShj+mfuMj6qkvEyfPp3GjRsTGBhIXFwcK1euPOXxs2fPplWrVgQGBtKuXTu++eabqoh5WizLYkLyBLKLsmlTtw13t7vbdCQRERGX8PWx82L/9vj72vlh+2E+/Xm/6UhAFZSXTz75hGHDhjF27FjWrFlDhw4d6NOnD+np6Sc8ftmyZdx6660MGTKEX375heuuu47rrruOjRs3ujrqaZm3ax6LDyzGz+6ncZGIiHi9ZhHBPHV56fho0rwtHDxmfnxks1z8GVBcXBwXXHABr732GgBOp5PY2FgeeeQRhg8f/pfjBwwYQG5uLvPmzSvb17VrVzp27MiMGTP+9vWys7MJDQ0lKyuLkJCQyvtBgLTcNK6fez05RTk8dv5j3NNOF+mKiIj3czgt+s9Yxpp9x+jVPJz/u7sLNputUl+jIu/fLv3kpaioiNWrVxMfH//HC9rtxMfHk5ycfMJzkpOTyx0P0KdPn5MeX1hYSHZ2drnNFRyFhcxPuIv4xVm0C2/HoHMHueR1RERE3I2P3cbU/h0I9IGQJf/Hex/P+/uTXMil5SUjIwOHw0FkZGS5/ZGRkaSmpp7wnNTU1Aodn5iYSGhoaNkWGxtbOeH/ZPnnr3PBd3vp/5OT8fUG42v3dcnriIiIuKMm9WoxNHw/TfL2kjrvfdIPHzWWxePvNkpISCArK6ts27/fNRcTdbl5KBlxzfF1gs+U17GKilzyOiIiIu4oY98eilclAdD6+juJqFfbWBaXlpfw8HB8fHxIS0srtz8tLY2oqKgTnhMVFVWh4wMCAggJCSm3uYKfjx/dXnoPn9q1Kdy6lYw333LJ64iIiLgbR0kJSW9Mw+kooUmnLlzX/1qjeVxaXvz9/enUqRMLFy4s2+d0Olm4cCHdunU74TndunUrdzzAggULTnp8VfKtW5eosWMAyHjzTfI3bTKcSERExPVWfTWHtF07CKxZi8vuHVrpF+tWlMvHRsOGDWPmzJl88MEHbNmyhQcffJDc3FwGDx4MwJ133klCQkLZ8Y899hhJSUm8+OKLbN26lXHjxvHzzz8zdOhQV0c9LSF9+xLcty+UlJCSMAKnxkciIuLF0vfsIvmzWQBcOvh+atWuYzhRFZSXAQMGMHXqVMaMGUPHjh1Zu3YtSUlJZRfl7tu3j5SUlLLju3fvzscff8xbb71Fhw4dmDNnDl9++SVt27Z1ddTTFjVmND516lC4fTsZr79uOo6IiIhLOEqKy8ZFzS7oSqueF5uOBFTBc16qmiuf81LudeZ/x8HHHgMfHxrPmkVQO/cpVyIiIpVh2eyPSJ7zbwKDQxg0dTo1w1x3ka7bPOfFm4X0uZyQq64Ch4NDCcM1PhIREa+StnsnK774FIDedz/g0uJSUSovZyFy1Eh8wsMp2rGTjFdfMx1HRESkUjhKikl6/WWcDgct4nrQslsv05HKUXk5C761axM9fhwAR955h/x168wGEhERqQTLP5tFxr49BAWH0HvIg8bvLvozlZezFNy7NyHXXgNOJ4cSRuAsLDQdSURE5Iyl7vyVFV/OBqD3kIeoERpmNtAJqLxUgqgRI/CtV4+iXbs4/MorpuOIiIickZLi0nGR5XTSslsvWnbraTrSCam8VAKfsDCiJowHIPPd98hb84vhRCIiIhWXPOdjjhzYR43QMC69+wHTcU5K5aWSBF9yCaHXXQeWRUpCAs78fNORRERETlvKjm2s+uozAOLveYgaIaGGE52cykslihyRgG9EBEV793J42j9NxxERETktJUVFJE1/Gcty0qrHRTTv0t10pFNSealEPiEhRE+aCEDm//0feT//bDiRiIjI31s2+yMyDx2gZlhtLh18v+k4f0vlpZLVuvBCQm+6ESyLQyNH4szLMx1JRETkpA5t38LP//kCgPh7hxIU7Lqn01cWlRcXiHzmGXyjoyneu4/0l6eZjiMiInJCxUWFJL0+Dcty0ubCS2nWOc50pNOi8uICPsHBRE8sHR8d/de/yF250nAiERGRv1o6618cTTlIrdp1uOSu+0zHOW0qLy5Sq2cPwm6+GYCUESNx5uYaTiQiIvKHA1s3sfqbrwC47L5HCKxVy3Ci06fy4kIR//gHfjExFB84QPqLL5mOIyIiAkBxYQHz35gGlsW5F8fT5PwLTEeqEJUXF/KpVZPoyZMAOPrxx+QuX244kYiICCz59/9xLDWFWnXDufjOe0zHqTCVFxer2a0bYbfeApSOjxzHNT4SERFz9m/ewJpv5wLQ575HCKzpOeOi36m8VIHIp57Cr359ig8dIv2FF0zHERGRaqq4oID5M0ofotru0stp3LGT4URnRuWlCthr1iR6yhQAjn3yCceXLjWcSEREqqMfP36frLRUguvW46I7PG9c9DuVlypSM64LtW+/HYCUUaNxHD9uOJGIiFQn+zauZ+38eQD0eeAxAmrUMJzozKm8VKGIYU/g17AhJSkppD/3nOk4IiJSTRTl55WNi9rH96VR+45mA50llZcqZK9Rg5gpk8Fm49jsORz/6SfTkUREpBr48aP3yD6cRki9CC66/W7Tcc6ayksVq9G5M3XuvAP4bXyUnW04kYiIeLO969eybsG3QOm4yD/Ic8dFv1N5MaDe44/j36gRJWlppD2r8ZGIiLhGYV4e898sHRd17HMVDdt2MJyocqi8GGAPCiI6MRFsNrI+/5ycxYtNRxIRES/0w4fvkJNxmNDIKHrdNsh0nEqj8mJIjfPPo86gQQCkjh6DIyvLbCAREfEqe9auZsPC+cBv46LAIMOJKo/Ki0H1HnsU/3POoeTwYdJ+ew6MiIjI2SrMy2X+W68CcN4V1xDbpp3hRJVL5cUge2AgMYlTwG4n66u55CxaZDqSiIh4gcX/9zbHj2QQFhVNr1vuMh2n0qm8GBbUsSN1h5TetpYydiwlR48aTiQiIp5s1y+r2Pj9ArDZ6PPg4/gFBpqOVOlUXtxA+NCh+DdriuNwBmmTNT4SEZEzU3D8OAveLB0XdbryWhq0OtdwItdQeXED9oAAYhITwceH7HnzyP7uO9ORRETEA33/wVscP5pJ7ej69Bhwh+k4LqPy4iaC2rWj7j2lX5KVOn6CxkciIlIhO1evYPOPi7DZ7PR96HH8ArxvXPQ7lRc3Ev7wQwQ0b47jyBHSJk40HUdERDxE/vEcFrz1GgCdrr6OmBatDSdyLZUXN2L39yf62d/GR998S3ZSkulIIiLiAb5/701yjx2lTkwDetx8u+k4Lqfy4maCzj2X8PvvA34bHx05YjiRiIi4s19XLmPLksWl46KHn8DX3990JJdTeXFD4Q88QEDLljiOHiV1/AQsyzIdSURE3FBedhb/fft1AC649gaim7U0nKhqqLy4IZu/PzHPJoKvLznffUf2N9+YjiQiIm5o0Xtvkpd1jLoNGtKt/0DTcaqMyoubCmzdmvAHHwAgbcJESg4fNpxIRETcyfblS9i27Edsdjt9H3oCXz8/05GqjMqLGwu/7z4C2rTGkZVFyrjxGh+JiAgAeVnHysZFXfr1J6ppc8OJqpbKixuz+fmVPrzOz4/jCxeSPW+e6UgiIuIGFr7zBvk52YQ3bEzXG28xHafKqby4ucCWLan38EMApE6aTHF6uuFEIiJi0rbkn9i+Yil2Hx/6Pvh4tRoX/U7lxQPUveceAs89F2dWFqljx2l8JCJSTeUeO8p/33kDgLjrbyaySTPDicxQefEANl9fYp5NxObnx/Hvvyfrq69MRxIRkSpmWRb/fft1CnKyqdfoHOKuv9l0JGNUXjxEQPPmhD/yCABpk6dQnJZmOJGIiFSlrUt/YMeq5NJx0UNP4ONb/cZFv1N58SB17x5MYPv2OHNySBk9WuMjEZFq4vjRTBa9OwOArjfeQkTjJoYTmaXy4kFsvr7EJE7B5u9P7o8/kfX556YjiYiIi1mWxYKZr1GQe5yIc5rSpV9/05GMU3nxMAFNm1LvsccASEt8luKUFMOJRETElbb89D27Vq/E7uPLFQ89gY+vr+lIxqm8eKA6g+4iqGNHnMePkzJK4yMREW+Vk5nBovffBKB7/9sIb9jYbCA34dLykpmZycCBAwkJCSEsLIwhQ4Zw/PjxUx7/yCOP0LJlS4KCgmjYsCGPPvooWVlZrozpcWw+PkRPmYItIIDcpUs5Nnu26UgiIlLJLMtiwVuvUZibS2ST5lxw7Y2mI7kNl5aXgQMHsmnTJhYsWMC8efP48ccfue+++056/KFDhzh06BBTp05l48aNvP/++yQlJTFkyBBXxvRIAU3Ood4TjwOQ/uxzFB88aDaQiIhUqk0/LGT3Lz/j4+vLFQ8/gd3Hx3Qkt2GzXDRz2LJlC23atGHVqlV07twZgKSkJK688koOHDhATEzMaf0+s2fP5vbbbyc3Nxff05jzZWdnExoaSlZWFiEhIWf1M7g7y+Fg7513kb96NTW6daXhu+9is9lMxxIRkbOUcySD9598iKL8PHrdNogu/W4yHcnlKvL+7bJPXpKTkwkLCysrLgDx8fHY7XZWrFhx2r/P7z/E6RSX6sbm40PMlMnYAgPJS17OsU8+MR1JRETOkmVZfPfmKxTl5xHdrCWdr7nedCS347LykpqaSkRERLl9vr6+1KlTh9TU1NP6PTIyMpg4ceIpR02FhYVkZ2eX26oT/0aNiBg2DIC051+g6MABw4lERORsbFj0HXvWrcHHz48+Dz2O3a5x0Z9VuLwMHz4cm812ym3r1q1nHSw7O5urrrqKNm3aMG7cuJMel5iYSGhoaNkWGxt71q/taWrfPpAanTtj5eWRMmIkltNpOpKIiJyB7Ix0fvjX2wD0HHAHdetXv/e001HhWcyTTz7JoEGDTnlMkyZNiIqKIv1P34BcUlJCZmYmUVFRpzw/JyeHvn37EhwczBdffIHfKb4xMyEhgWG/ffIApaWnuhUYm91OdOIUdvW7jryVKzn68b+pc/tA07FERKQCLMti/oxXKMrPJ6ZFa86/qp/pSG6rwuWlXr161KtX72+P69atG8eOHWP16tV06tQJgEWLFuF0OomLizvpednZ2fTp04eAgADmzp1LYGDgKV8nICCAgICAiv0QXsg/NpaIp54kbcJE0l98kVoX9sK/YUPTsURE5DSt/28S+zasxdc/gD4Palx0Ki675qV169b07duXe++9l5UrV7J06VKGDh3KLbfcUnan0cGDB2nVqhUrV64ESovL5ZdfTm5uLu+88w7Z2dmkpqaSmpqKw+FwVVSvUfuWW6gRF4eVn8+hESM0PhIR8RBZ6an88K93AOh5y53UialvOJF7c+lzXj766CNatWpF7969ufLKK+nZsydvvfVW2a8XFxezbds28vLyAFizZg0rVqxgw4YNNGvWjOjo6LJt//79rozqFWx2O9GTJ2OvUYP8n1dz9MMPTUcSEZG/YTmdzJ/xCsWFBdRvdS7nX3GN6Uhuz2XPeTGlOj3n5WSOzvqE1HHjsAUG0uTLL/Bv3Nh0JBEROYm1879m4btv4BsQwJ3Pv0rtqNN7Dpq3cYvnvIg5YQNupmb3blgFBRwaMRJLIzcREbd0LC2VHz56F4ALbxtUbYtLRam8eCGbzUb0pEnYa9Ykf80aMv/vX6YjiYjIn1hOJ/PfmEZJYSGxbdrR8fKrTEfyGCovXsovJoaI4c8AcHjaNAp37TKcSERE/tcv8+dxYMtG/AICufyBx7DZ9ZZ8urRSXizsppuo2bMnVmEhhxISND4SEXETR1MP8dPHHwBw4e13ExZ56uefSXkqL16sdHw0EXtwMAXr1pP53numI4mIVHtOp4Ok16dRUlRIw7Yd6BDf13Qkj6Py4uX8oqKIHD4cgMOvvErhjh2GE4mIVG+/fPsfDm3bjF9gEH00LjojWrFqIPSG66l50YVYRUUcShiBVVJiOpKISLWUeegAS/79fwBcfMcQQupF/M0ZciIqL9WAzWYjesIE7CEhFGzYwJF33jUdSUSk2nE6HSS9MY2S4iIatT+Pdr37mI7ksVReqgm/yEiiRo4A4PBrr1GwfbvhRCIi1cvqr78iZftW/INqcPn9j2Kz2UxH8lgqL9VIyLXXUuvSS6G4mJSEEVjFxaYjiYhUC0cO7GfpJ6XP3Lr4rnsICf/7LziWk1N5qUZsNhtR48ZiDw2lYNMmjrz9tulIIiJez+lwkPTGyziKizmnYyfaXnyZ6UgeT+WlmvGLiCBq1CgADr/+BgVbtxpOJCLi3Vb953NSd2wnoEZNLrvvEY2LKoHKSzUUcvVVBF8WD8XFpXcfFRWZjiQi4pUy9u8lefZHAFwy6D6C64YbTuQdVF6qIZvNRtTYsfiEhVG4ZQsZb75lOpKIiNdxlJSQ9PrLOEpKaHL+BbS58FLTkbyGyks15RseTtSY0QBkvPkmBZs3G04kIuJdVs39jLRdOwioWZPL7h2qcVElUnmpxoKvuILgPn2gpIRDwxM0PhIRqSSH9+4mec6/Abh08APUqlPXcCLvovJSjZWOj8bgU6cOhdu3c/iNN0xHEhHxeKXjomk4HSU07dyV1j0vNh3J66i8VHO+deoQNXYsAEfemkn+ho2GE4mIeLaVX84mfc9OAmsFc9m9D2tc5AIqL0JIn8sJufIKcDhIGZGAU+MjEZEzkr5nF8s/nwXApXc/QM2w2oYTeSeVFwEgcvRofOrWpfDXHWS8Nt10HBERj+MoKSZp+ks4HQ6ad+lOq+4Xmo7ktVReBADf2rWJGvfb+Ojtt8lfv95wIhERz7L88084vG8PQcEhxN/zkMZFLqTyImVCLruMkGuuAaeTQ8MTcBYWmo4kIuIR0nbtYMUXnwLQe8hD1AgNMxvIy6m8SDlRI0fgUy+col27OPzKK6bjiIi4vZLiYpJefxnL6aRFt1607NbTdCSvp/Ii5fiEhRE9fgIAme+9T94vvxhOJCLi3pZ/9m8y9u8lKCSU3nc/YDpOtaDyIn8RfOklhPbrB04nKQkjcBYUmI4kIuKWUndsZ+WXcwCIv+chaoSEGk5UPai8yAlFjkjANyKCoj17ODztn6bjiIi4nZKiIpLemIZlOWnV4yJaxPUwHanaUHmRE/IJDSV64m/jow8+IG/1asOJRETcy7I5H3PkwD5qhIZx6eD7TcepVlRe5KRqXXQRoTfcAJbFoREjcObnm44kIuIWDm3fys9zPwfgsnuHEhQcYjhR9aLyIqcUOfwZfKOiKN67j/SXXzYdR0TEuOKiwrJxUetel9Dsgq6mI1U7Ki9ySj4hIURPnAjA0f/7F7krVxpOJCJi1tJPPuTooQPUrF2HSwbdZzpOtaTyIn+rVq+ehPXvD0DKyFE4c3MNJxIRMePgti2s/vpL4LdxUa1gs4GqKZUXOS0Rz/wD35hoivfvJ/3Fl0zHERGpcsWFBcx/42WwLM69KJ6mnbqYjlRtqbzIafGpVYuYSZMAOPrxx+QuX2E4kYhI1Voy618cTTlErTp1ufiue0zHqdZUXuS01ezenbBbBgCQMnIkjuMaH4lI9XBgy0bWfDsXgMvve4TAmrUMJ6reVF6kQiKeehq/+vUpPniQ9KkvmI4jIuJyxQUFzH/jn2BZtL3kcs45r7PpSNWeyotUiE+tmkRPngzAsVmfcHzpUsOJRERc66d/f8CxtBSC69bj4juHmI4jqLzIGajZNY7at90GQMqo0TiOHzecSETENfZvWs8vSf8B4PL7HyGgRk3DiQRUXuQMRTw5DL/YWEpSUkh/7nnTcUREKl1RQT5Jb5R+t1v73n1p3OF8w4nkdyovckbsNWsSM+W38dHs2Rz/aYnhRCIilevHj94n+3AaIfUiuOiOu03Hkf+h8iJnrMYFF1D7zjsASBk1Ckd2tuFEIiKVY++Gtaz77msALr//UfyDahhOJP9L5UXOSsQTT+DXqCElaWmkPfuc6TgiImetMC+P+TNKx0UdLruSRu06mg0kf6HyImfFHhREzJQpYLOR9fnnHP/hB9ORRETOyo8fvktOxmFC6kVy4e2DTceRE1B5kbNWo1Mn6tx1FwApo8fgyMoynEhE5MzsWbeG9QuTAOj74GP4BwYZTiQnovIilaLe44/h37gxJenppE1JNB1HRKTCCvNy+e7NVwE4r+81xJ7b3nAiORmVF6kU9sBAYp5NBLudrK++ImfRItORREQqZPH/vUPOkcOERUbT69a7TMeRU1B5kUoT1LEjdQYPAiBl7Fgcx44ZzSMicrp2//IzG7//Dmw2+jz4GH6BgaYjySmovEilqvfoo/g3aYLjcAapk6eYjiMi8rcKco/z3ZuvAHD+FdfSoHVbw4nk76i8SKWyBwSUjY+y//MfshcsMB1JROSUFn/wNsePZlI7Ooaet9xhOo6cBpeWl8zMTAYOHEhISAhhYWEMGTKE46f5PTiWZXHFFVdgs9n48ssvXRlTKllQ+/bUveceAFLHjafk6FHDiURETmzn6pVs+uG/v42LnsAvQOMiT+DS8jJw4EA2bdrEggULmDdvHj/++CP33XffaZ07bdo0bDabK+OJC4UPfZiA5s1wHDlC2sSJpuOIiPxF/vEcFsx8DYBOV11H/ZatDSeS0+Wy8rJlyxaSkpJ4++23iYuLo2fPnrz66qvMmjWLQ4cOnfLctWvX8uKLL/Luu++6Kp64mN3fn+gpieDjQ/Y335KdNN90JBGRcr5//y1yj2ZSO6YBPQbcbjqOVIDLyktycjJhYWF07ty5bF98fDx2u50VK1ac9Ly8vDxuu+02pk+fTlRU1N++TmFhIdnZ2eU2cQ9B7dpS9757AUgdP56SI0cMJxIRKbVj1XK2/PQ9Npudvg8+jp9/gOlIUgEuKy+pqalERESU2+fr60udOnVITU096XlPPPEE3bt3p1+/fqf1OomJiYSGhpZtsbGxZ5VbKle9Bx8koGVLHEePkjp+ApZlmY4kItVcfk522bio87U3ENOileFEUlEVLi/Dhw/HZrOdctu6desZhZk7dy6LFi1i2rRpp31OQkICWVlZZdv+/fvP6LXFNWz+/sQkTgFfX3K++46cb781HUlEqrlF771JXtYx6jZoSPebbjMdR86Ab0VPePLJJxk0aNApj2nSpAlRUVGkp6eX219SUkJmZuZJx0GLFi1i586dhIWFldt/44030qtXLxYvXvyXcwICAggI0Md97iywTRvC77+fjOnTSZ0wkRpduuAbHm46lohUQ9tXLGXr0h+w2UvHRb7+/qYjyRmwWS76HH/Lli20adOGn3/+mU6dOgHw3Xff0bdvXw4cOEBMTMxfzklNTSUjI6Pcvnbt2vHPf/6Ta665hnPOOedvXzc7O5vQ0FCysrIICQmpnB9GzppVVMTuAbdQuGULwZfFU/+VV3Q3mYhUqbzsLN5/8iHys7OIu/5met5yp+lI8j8q8v7tsmteWrduTd++fbn33ntZuXIlS5cuZejQodxyyy1lxeXgwYO0atWKlStXAhAVFUXbtm3LbQANGzY8reIi7qtsfOTnR86C/5I972vTkUSkmln47gzys7MIj21E1xtvNR1HzoJLn/Py0Ucf0apVK3r37s2VV15Jz549eeutt8p+vbi4mG3btpGXl+fKGOImAlu1ot5DDwKQOmkSxX8aK4qIuMq25J/YnvxT6bjooSfw9fMzHUnOgsvGRqZobOTerOJi9gy4hYLNm6l1ySU0eH26xkci4lK5x47y/lMPU5CTTdcbb6HHzXqmiztyi7GRyInY/PyIfjYR/Pw4/v33ZM+dazqSiHgxy7JY+M4bFORkU69hY7reMMB0JKkEKi9S5QJbtKDe0KEApE6eQnFamuFEIuKtti37kV9XLsPu40Pfh4fh46txkTdQeREj6g65m8B27XBmZ5MyZoweXicilS732FEWvjsDgK433EJE4yaGE0llUXkRI2y+vsQkTsHm50fuDz+S9fkXpiOJiBexLIsFM6dTcDyHiMZN6XJdf9ORpBKpvIgxAc2aUe+xRwFIS0ykOCXFcCIR8RZblixm58/Lsfv40vehx/HxrfAzWcWNqbyIUXUGDyaoQwecx4+TMlrjIxE5e8czj7DovdJxUbebbqVeIz0nzNuovIhRNh8fohMTsQUEkLtkCcfmzDEdSUQ8WOm46DUKc3OJbNKMLv1uMh1JXEDlRYwLaHIO9R5/HID0Z5+j+OBBs4FExGNt/nERu9aswsfXl74PPYHdx8d0JHEBlRdxC3XuvIOg88/HmZtLyujRGh+JSIXlHMng+/dLn+Lerf9AwmMbGU4krqLyIm7B5uND9ORJ2AIDyV2WzLFPPjUdSUQ8iGVZfPfWqxTm5RLVrAUXXHOD6UjiQiov4jYCzjmHiGFPAJD2/PMUHThgOJGIeIqNixewZ+1qfPz86PugxkXeTuVF3Ert228nqHMnrLw8UkaOwnI6TUcSETeXnZHO4g/eBqDHgDuo2yDWcCJxNZUXcSs2u52YKVOwBQWRt2IFR//9b9ORRMSNWZbFd2++SlF+HtEtWtHpqn6mI0kVUHkRt+PfsCERTz4JQPrUFynat89wIhFxVxsWzmfv+l/w9fOn74OPY7drXFQdqLyIW6p9263U6NIFKz+flBEjNT4Skb/IPpzO4n+9A0DPW++kTkwDw4mkqqi8iFuy2e1ET5mMrUYN8n7+maMffmQ6koi4EcvpZP6MaRQX5FO/VRvOu+Ia05GkCqm8iNvyb9CAyH88DUD6Sy9RtGeP2UAi4jbW/TeJfRvX4+sfQB+Ni6odlRdxa2EDBlCjW1esggIOjRiJ5XCYjiQihh1LS+XHD98FoNdtd1E7KsZwIqlqKi/i1mw2GzGTJmGvUYP8NWvI/Ne/TEcSEYPKxkWFBTRo3Zbz+lxtOpIYoPIibs+vfn0ihj8DwOGXp1G4a7fhRCJiytrvvubA5o34BQTS58HHsdn1NlYd6U9dPEJY//7U7NEDq7CQlIQEjY9EqqGjqYf48eP3Abhw4GDCIqPMBhJjVF7EI9hsNqInTcReqxb569aR+f77piOJSBWynE7mvzGNksJCYs9tT4fLrjAdSQxSeRGP4RcdTWTCcAAO//MVCnfuNJxIRKrKmm//w8Gtm/ELDKLPA49pXFTN6U9fPEroDTdQ88JeWEVFHEoYgVVSYjqSiLhY5qGDLJn1fwBcdPvdhEZEGk4kpqm8iEex2WxET5yIPTiYgvXrOfLue6YjiYgLOZ2O0nFRUSEN23WkfXxf05HEDai8iMfxi4wkcsQIADJefZWC7dsNJxIRV1nz9Vcc2r4F/6Ag+jzwKDabzXQkcQMqL+KRQq/rR62LL8YqLiYlYQRWcbHpSCJSyY4c3M+ST0qf7XTRHfcQEh5hOJG4C5UX8Ug2m42o8eOxh4ZSsGkTR955x3QkEalETqeD+a9Pw1FcTOMO59Pu0stNRxI3ovIiHssvMoKoUSMBODz9dQq2bTOcSEQqy8//+YKUHdsIqFGTy+/XuEjKU3kRjxZy9dXUiu8NxcUcGp6g8ZGIFzhyYB/LPv0QgIvvupfguuGGE4m7UXkRj2az2YgeOxaf0FAKt2wh4823TEcSkbPgdDj4dvrLOEpKaHL+BZx7UW/TkcQNqbyIx/OtV4/IMaMByJgxg4ItWwwnEpEztWruZ6Tt+pWAmjWJv/dhjYvkhFRexCuEXHklwZdfDiUlpeOjoiLTkUSkgg7v28Oy2R8DcOmg+wmuo3GRnJjKi3gFm81G1Ngx+NSuTeG2bWTMmGE6kohUgKOkhKTXX8bpKKFp5zha97rEdCRxYyov4jV869YlauwYADLefIv8jZsMJxKR07Xyq9mk795JYM1axN+jcZGcmsqLeJWQvn0JvqIvOBykJAzHqfGRiNtL37OL5Z/NAuDSux+gVu06hhOJu1N5Ea8TNWYMPnXrUvjrDjKmv246joicgqOkmKQ3puF0OGh2QTda9bjIdCTxACov4nV8a9cmatxYAI7MnEn+hg2GE4nIyaz44lMO79lFYHAI8fc8pHGRnBaVF/FKIZddRsjVV4PTyaHhCTgLC01HEpE/Sdu1gxVffApA/JAHqRlW23Ai8RQqL+K1IkeOwCc8nKKdO8l49VXTcUTkf5QU/zEuahHXg5bdepmOJB5E5UW8lm/t2kSPHwfAkXffI3/tWqN5ROQPyz+bRca+PQSFhNL7nodMxxEPo/IiXi24d29C+11bOj5KGIGzoMB0JJFqL3Xnr6z8ajYA8fc8RI2QUMOJxNP4mg4g4mqRI0aQuyyZwt27SXn5NSKGPWE6kki1VVJczLfTX8JyOmnZ/UJaxPUwHUk8kMqLeD2f0FCiJk5g70OP8eXe8+GxH0xHEqm2ivN/wlGwn6CQUC4dfL/pOOKhNDaSaiH44osJufYa0zFEqjWnIxNHwc8AXDroQY2L5IzpkxepNmISnubK4aMJH/oQgc2am44jUu1YlsX25XVJ+XUrLbtrXCRnzmZZlmU6RGXKzs4mNDSUrKwsQkJCTMcRERGR01CR92+XjY0yMzMZOHAgISEhhIWFMWTIEI4fP/635yUnJ3PppZdSs2ZNQkJCuPDCC8nPz3dVTBEREfEwLisvAwcOZNOmTSxYsIB58+bx448/ct99953ynOTkZPr27cvll1/OypUrWbVqFUOHDsVu16U5IiIiUsolY6MtW7bQpk0bVq1aRefOnQFISkriyiuv5MCBA8TExJzwvK5du3LZZZcxceLEM35tjY1EREQ8j/GxUXJyMmFhYWXFBSA+Ph673c6KFStOeE56ejorVqwgIiKC7t27ExkZyUUXXcSSJUtO+VqFhYVkZ2eX20RERMR7uaS8pKamEhERUW6fr68vderUITU19YTn7Nq1C4Bx48Zx7733kpSUxPnnn0/v3r359ddfT/paiYmJhIaGlm2xsbGV94OIiIiI26lQeRk+fDg2m+2U29atW88oiNPpBOD+++9n8ODBnHfeebz88su0bNmSd99996TnJSQkkJWVVbbt37//jF5fREREPEOFnvPy5JNPMmjQoFMe06RJE6KiokhPTy+3v6SkhMzMTKKiok54XnR0NABt2rQpt79169bs27fvpK8XEBBAQEDAaaQXERERb1Ch8lKvXj3q1av3t8d169aNY8eOsXr1ajp16gTAokWLcDqdxMXFnfCcxo0bExMTw7Zt28rt3759O1dccUVFYoqIiIgXc8k1L61bt6Zv377ce++9rFy5kqVLlzJ06FBuueWWsjuNDh48SKtWrVi5ciUANpuNp59+mldeeYU5c+awY8cORo8ezdatWxkyZIgrYoqIiIgHctnXA3z00UcMHTqU3r17Y7fbufHGG3nllVfKfr24uJht27aRl5dXtu/xxx+noKCAJ554gszMTDp06MCCBQto2rSpq2KKiIiIh9HXA4iIiIhxxp/zIiIiIuIqKi8iIiLiUVx2zYspv0/B9KRdERERz/H7+/bpXM3ideUlJycHQE/aFRER8UA5OTmEhoae8hivu2DX6XRy6NAhgoODsdlspuOcluzsbGJjY9m/f78uMkbr8Wdajz9oLcrTepSn9fiDJ66FZVnk5OQQExOD3X7qq1q87pMXu91OgwYNTMc4IyEhIR7zl6wqaD3K03r8QWtRntajPK3HHzxtLf7uE5ff6YJdERER8SgqLyIiIuJRVF7cQEBAAGPHjtUXTP5G61Ge1uMPWovytB7laT3+4O1r4XUX7IqIiIh30ycvIiIi4lFUXkRERMSjqLyIiIiIR1F5EREREY+i8mJIZmYmAwcOJCQkhLCwMIYMGcLx48f/9rzk5GQuvfRSatasSUhICBdeeCH5+flVkNi1znQ9oPSpjFdccQU2m40vv/zStUGrQEXXIjMzk0ceeYSWLVsSFBREw4YNefTRR8nKyqrC1JVn+vTpNG7cmMDAQOLi4li5cuUpj589ezatWrUiMDCQdu3a8c0331RR0qpRkfWYOXMmvXr1onbt2tSuXZv4+Pi/XT9PU9G/H7+bNWsWNpuN6667zrUBq1BF1+LYsWM8/PDDREdHExAQQIsWLTz3/y+WGNG3b1+rQ4cO1vLly62ffvrJatasmXXrrbee8pxly5ZZISEhVmJiorVx40Zr69at1ieffGIVFBRUUWrXOZP1+N1LL71kXXHFFRZgffHFF64NWgUquhYbNmywbrjhBmvu3LnWjh07rIULF1rNmze3brzxxipMXTlmzZpl+fv7W++++661adMm695777XCwsKstLS0Ex6/dOlSy8fHx3r++eetzZs3W6NGjbL8/PysDRs2VHFy16joetx2223W9OnTrV9++cXasmWLNWjQICs0NNQ6cOBAFSd3jYqux+92795t1a9f3+rVq5fVr1+/qgnrYhVdi8LCQqtz587WlVdeaS1ZssTavXu3tXjxYmvt2rVVnLxyqLwYsHnzZguwVq1aVbbv22+/tWw2m3Xw4MGTnhcXF2eNGjWqKiJWqTNdD8uyrF9++cWqX7++lZKS4hXl5WzW4n99+umnlr+/v1VcXOyKmC7TpUsX6+GHHy773w6Hw4qJibESExNPePzNN99sXXXVVeX2xcXFWffff79Lc1aViq7Hn5WUlFjBwcHWBx984KqIVepM1qOkpMTq3r279fbbb1t33XWX15SXiq7FG2+8YTVp0sQqKiqqqogupbGRAcnJyYSFhdG5c+eyffHx8djtdlasWHHCc9LT01mxYgURERF0796dyMhILrroIpYsWVJVsV3mTNYDIC8vj9tuu43p06cTFRVVFVFd7kzX4s+ysrIICQnB19dzvr6sqKiI1atXEx8fX7bPbrcTHx9PcnLyCc9JTk4udzxAnz59Tnq8JzmT9fizvLw8iouLqVOnjqtiVpkzXY8JEyYQERHBkCFDqiJmlTiTtZg7dy7dunXj4YcfJjIykrZt2zJlyhQcDkdVxa5UKi8GpKamEhERUW6fr68vderUITU19YTn7Nq1C4Bx48Zx7733kpSUxPnnn0/v3r359ddfXZ7Zlc5kPQCeeOIJunfvTr9+/Vwdscqc6Vr8r4yMDCZOnMh9993nioguk5GRgcPhIDIystz+yMjIk/7sqampFTrek5zJevzZM888Q0xMzF8Knic6k/VYsmQJ77zzDjNnzqyKiFXmTNZi165dzJkzB4fDwTfffMPo0aN58cUXmTRpUlVErnQqL5Vo+PDh2Gy2U25bt249o9/b6XQCcP/99zN48GDOO+88Xn75ZVq2bMm7775bmT9GpXHlesydO5dFixYxbdq0yg3tIq5ci/+VnZ3NVVddRZs2bRg3btzZBxeP9eyzzzJr1iy++OILAgMDTcepcjk5Odxxxx3MnDmT8PBw03GMczqdRERE8NZbb9GpUycGDBjAyJEjmTFjhuloZ8RzPlP2AE8++SSDBg065TFNmjQhKiqK9PT0cvtLSkrIzMw86fgjOjoagDZt2pTb37p1a/bt23fmoV3IleuxaNEidu7cSVhYWLn9N954I7169WLx4sVnkbzyuXItfpeTk0Pfvn0JDg7miy++wM/P72xjV6nw8HB8fHxIS0srtz8tLe2kP3tUVFSFjvckZ7Iev5s6dSrPPvss//3vf2nfvr0rY1aZiq7Hzp072bNnD9dcc03Zvt//Eejr68u2bdto2rSpa0O7yJn83YiOjsbPzw8fH5+yfa1btyY1NZWioiL8/f1dmrnSmb7opjr6/aLMn3/+uWzf/PnzT3lRptPptGJiYv5ywW7Hjh2thIQEl+Z1tTNZj5SUFGvDhg3lNsD65z//ae3atauqole6M1kLy7KsrKwsq2vXrtZFF11k5ebmVkVUl+jSpYs1dOjQsv/tcDis+vXrn/KC3auvvrrcvm7dunnVBbsVWQ/LsqznnnvOCgkJsZKTk6siYpWqyHrk5+f/5b8R/fr1sy699FJrw4YNVmFhYVVGr3QV/buRkJBgNWrUyHI4HGX7pk2bZkVHR7s8qyuovBjSt29f67zzzrNWrFhhLVmyxGrevHm522EPHDhgtWzZ0lqxYkXZvpdfftkKCQmxZs+ebf3666/WqFGjrMDAQGvHjh0mfoRKdSbr8Wd4wd1GllXxtcjKyrLi4uKsdu3aWTt27LBSUlLKtpKSElM/xhmZNWuWFRAQYL3//vvW5s2brfvuu88KCwuzUlNTLcuyrDvuuMMaPnx42fFLly61fH19ralTp1pbtmyxxo4d63W3SldkPZ599lnL39/fmjNnTrm/Bzk5OaZ+hEpV0fX4M2+626iia7Fv3z4rODjYGjp0qLVt2zZr3rx5VkREhDVp0iRTP8JZUXkx5MiRI9att95q1apVywoJCbEGDx5c7j8wu3fvtgDr+++/L3deYmKi1aBBA6tGjRpWt27drJ9++qmKk7vGma7H//KW8lLRtfj+++8t4ITb7t27zfwQZ+HVV1+1GjZsaPn7+1tdunSxli9fXvZrF110kXXXXXeVO/7TTz+1WrRoYfn7+1vnnnuu9fXXX1dxYteqyHo0atTohH8Pxo4dW/XBXaSifz/+lzeVF8uq+FosW7bMiouLswICAqwmTZpYkydP9rh/4PzOZlmWVbWDKhEREZEzp7uNRERExKOovIiIiIhHUXkRERERj6LyIiIiIh5F5UVEREQ8isqLiIiIeBSVFxEREfEoKi8iIiLiUVReRERExKOovIiIiIhHUXkRERERj6LyIiIiIh7l/wGSxL9Ux81ksAAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig,ax=model.unit_cell.plot()\n",
    "ax.scatter(*model.phases.detach().T,s=1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "bce15e38-85d0-4e77-a7dd-cc47f5c572ac",
   "metadata": {
    "collapsed": true,
    "jupyter": {
     "outputs_hidden": true
    },
    "tags": []
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "100%|██████████████████████████████████████████████████████████| 2000/2000 [00:14<00:00, 142.73it/s]\n"
     ]
    }
   ],
   "source": [
    "# Train model\n",
    "nsamples_space = 2**10\n",
    "loss_history = []\n",
    "phases_history = [model.phases.clone().detach().numpy()]\n",
    "scale_history = []\n",
    "for epoch in tqdm.trange(2000):\n",
    "    rs = model.unit_cell.sample(nsamples_space)\n",
    "    rs = torch.tensor(rs,dtype=dtype)\n",
    "    loss = model.train_step(rs)\n",
    "    loss_history.append(loss)\n",
    "    phases_history.append(model.phases.clone().detach().numpy())\n",
    "    scale_history.append(model.scale.clone().detach().numpy())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "25de3ad7-99b5-472f-a403-2515c2b9a1f5",
   "metadata": {
    "collapsed": true,
    "jupyter": {
     "outputs_hidden": true
    },
    "tags": []
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "445 µs ± 105 µs per loop (mean ± std. dev. of 7 runs, 20 loops each)\n",
      "6.57 ms ± 76.6 µs per loop (mean ± std. dev. of 7 runs, 20 loops each)\n"
     ]
    }
   ],
   "source": [
    "# test time to sample and time to train to see where comp. bottleneck is.\n",
    "%timeit -n 20 model.unit_cell.sample(nsamples_space)\n",
    "%timeit -n 20 model.train_step(rs)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "526f3917-5385-4311-b326-ff227c5c862f",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "3.1031556253140648e+19"
      ]
     },
     "execution_count": 20,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAiMAAAGsCAYAAAAPJKchAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy88F64QAAAACXBIWXMAAA9hAAAPYQGoP6dpAAAoAUlEQVR4nO3dfVjVdZ7/8deB5CCrgEYcxFC0G81UVFwZauxmOyual+Xc7LLmhsOUbqXXOtGYQzeQNQndaO41Qzo5kbM7lczMVbbX6uIoyZrjmRxRKkspb3FMQHO4UQsQPr8/+nnqJChHD3048Hxc17ku+Hw/3+95f/rIOa++5/v5HocxxggAAMCSENsFAACAno0wAgAArCKMAAAAqwgjAADAKsIIAACwijACAACsIowAAACrCCMAAMAqwggAALCKMAIAAKwKqjCyefNmTZs2TfHx8XI4HFqzZo1f+1dUVOjWW2+Vy+VSeHi4hg4dqscee0zNzc3ePitXrtTEiRPVr18/9evXT263W9u2bQvwSAAAwFlBFUZOnTqlpKQkFRQUXNT+vXr1UkZGhv74xz+qoqJCy5Yt08qVK5Wbm+vtU1paqhkzZmjTpk3yeDxKSEjQpEmTdOTIkUANAwAAfI0jWL8oz+Fw6M0339T06dO9bY2NjXr00Uf1+uuvq7a2ViNHjtQzzzyjW265pd3jZGVl6S9/+YveeeedNre3tLSoX79++uUvf6mMjIwAjwIAAATVmZELmTdvnjwej1avXq33339f//RP/6TJkyfrk08+abP/3r17VVxcrJtvvrndY54+fVrNzc3q379/Z5UNAECP1m3OjFRWVmro0KGqrKxUfHy8t5/b7daECRO0ePFib9sNN9ygHTt2qLGxUXPmzNHy5csVEtJ2LnvggQe0fv16ffjhhwoPD+/UMQEA0BNdZruAQPnggw/U0tKia6+91qe9sbFRl19+uU9bUVGRGhoa9N5772nBggV6/vnn9fDDD59zzPz8fK1evVqlpaUEEQAAOkm3CSMnT55UaGioysrKFBoa6rOtT58+Pr8nJCRIkkaMGKGWlhbNmTNHDz30kM9+zz//vPLz87Vx40aNHj268wcAAEAP1W3CyNixY9XS0qKamhpNnDixw/u1traqublZra2t3jDy7LPP6umnn9b69es1fvz4zioZAAAoyMLIyZMntXfvXu/vBw4cUHl5ufr3769rr71WM2fOVEZGhpYsWaKxY8fq2LFjKikp0ejRozV16lS9+uqr6tWrl0aNGiWn06nt27crOztb6enp6tWrlyTpmWeeUU5Ojl577TUlJiaqqqpK0pdnV755hgUAAFy6oLqAtbS0VLfeeus57bNmzdKqVavU3Nysn//85/rP//xPHTlyRDExMfrOd76jRYsWadSoUSoqKtKzzz6rjz/+WMYYDR48WP/6r/+qBx980HtNSGJiog4dOnTOc+Tm5uqJJ57o7CECANDjBFUYAQAA3U+3us8IAAAIPoQRAABgVVBcwNra2qpPP/1Uffv2lcPhsF0OAADoAGOMGhoaFB8f3+7NRaUgCSOffvqp994gAAAguBw+fFhXXnllu9uDIoz07dtX0peDiYyMtFwNAADoiPr6eiUkJHjfx9sTFGHk7EczkZGRhBEAAILMhS6x4AJWAABgFWEEAABYRRgBAABWEUYAAIBVhBEAAGAVYQQAAFhFGAEAAFYRRgAAgFWEEQAAYBVhBAAAWEUYAQAAVhFGAACAVYQRIICMMfovz0GVHfqb7VIAIGgExbf2AsHi7T01evytDyVJB/OnWq4GAIIDZ0aAANp/7JTtEgAg6BBGAACAVYQRAABgFWEEAABYRRgBAABWEUYAAIBVhBEAAGAVYQQAAFhFGAEAAFYRRgAAgFWEEQAAYBVhBAAAWEUYAQAAVvkdRjZv3qxp06YpPj5eDodDa9asueA+paWlGjdunJxOp66++mqtWrXqIkoFAADdkd9h5NSpU0pKSlJBQUGH+h84cEBTp07VrbfeqvLycv3kJz/Rvffeq/Xr1/tdLAAA6H4u83eHKVOmaMqUKR3uv2LFCg0ZMkRLliyRJF133XXasmWLXnjhBaWlpfn79AAAoJvp9GtGPB6P3G63T1taWpo8Hk+7+zQ2Nqq+vt7nAQAAuqdODyNVVVVyuVw+bS6XS/X19fr888/b3CcvL09RUVHeR0JCQmeXCQAALOmSq2mys7NVV1fnfRw+fNh2SQAAoJP4fc2Iv+Li4lRdXe3TVl1drcjISPXu3bvNfZxOp5xOZ2eXBgAAuoBOPzOSmpqqkpISn7YNGzYoNTW1s58aAAAEAb/DyMmTJ1VeXq7y8nJJXy7dLS8vV2VlpaQvP2LJyMjw9r/vvvu0f/9+Pfzww9qzZ49efPFF/e53v9ODDz4YmBEAAICg5ncY2b59u8aOHauxY8dKkrKysjR27Fjl5ORIko4ePeoNJpI0ZMgQrV27Vhs2bFBSUpKWLFmiX//61yzrBQAAki7impFbbrlFxph2t7d1d9VbbrlFO3fu9PepAABAD9AlV9MAAICegzACAACsIowAAACrCCMAAMAqwggAALCKMAIAAKwijAAAAKsIIwAAwCrCCAAAsIowAgAArCKMAAAAqwgjAADAKsIIAACwijACAACsIowAAACrCCNAABkZ2yUAQNAhjAAAAKsII0AAOeSwXQIABB3CCAAAsIowAgAArCKMAAAAqwgjAADAKsIIAACwijACAACsIowAAACrCCMAAMAqwggAALCKMAIAAKwijAAAAKsIIwAAwCrCCAAAsIowAgAArCKMAAAAqwgjAADAKsIIAACwijACAACsIowAAACrCCMAAMAqwggAALCKMAIAAKwijAAAAKsIIwAAwCrCCAAAsIowAgAArCKMAAAAqwgjAADAKsIIAACwijACAACsIowAAACrCCMAAMAqwggAALCKMAIAAKy6qDBSUFCgxMREhYeHKyUlRdu2bTtv/2XLlmnYsGHq3bu3EhIS9OCDD+qLL764qIIBAED34ncYKSoqUlZWlnJzc7Vjxw4lJSUpLS1NNTU1bfZ/7bXX9LOf/Uy5ubnavXu3Xn75ZRUVFemRRx655OIBAEDw8zuMLF26VLNnz1ZmZqZGjBihFStWKCIiQoWFhW3237p1q2688UbdddddSkxM1KRJkzRjxowLnk0BAAA9g19hpKmpSWVlZXK73V8dICREbrdbHo+nzX1uuOEGlZWVecPH/v37tW7dOt1+++3tPk9jY6Pq6+t9HgAAoHu6zJ/Ox48fV0tLi1wul0+7y+XSnj172tznrrvu0vHjx/Xd735XxhidOXNG991333k/psnLy9OiRYv8KQ0AAASpTl9NU1paqsWLF+vFF1/Ujh079MYbb2jt2rV66qmn2t0nOztbdXV13sfhw4c7u0wAAGCJX2dGYmJiFBoaqurqap/26upqxcXFtbnP448/rrvvvlv33nuvJGnUqFE6deqU5syZo0cffVQhIefmIafTKafT6U9pAAAgSPl1ZiQsLEzJyckqKSnxtrW2tqqkpESpqalt7nP69OlzAkdoaKgkyRjjb70AAKCb8evMiCRlZWVp1qxZGj9+vCZMmKBly5bp1KlTyszMlCRlZGRo4MCBysvLkyRNmzZNS5cu1dixY5WSkqK9e/fq8ccf17Rp07yhBAAA9Fx+h5H09HQdO3ZMOTk5qqqq0pgxY1RcXOy9qLWystLnTMhjjz0mh8Ohxx57TEeOHNEVV1yhadOm6emnnw7cKAAAQNBymCD4rKS+vl5RUVGqq6tTZGSk7XKAdq3cvF9Pr9stSTqYP9VyNQBgV0ffv/luGgAAYBVhBAAAWEUYAQAAVhFGAACAVYQRIICMuvz14ADQ5RBGAACAVYQRIIAcctguAQCCDmEEAABYRRgBAABWEUYAAIBVhBEAAGAVYQQAAFhFGAEAAFYRRgAAgFWEEQAAYBVhBAAAWEUYAQAAVhFGAACAVYQRAABgFWEEAABYRRgBAABWEUYAAIBVhBEAAGAVYQQAAFhFGAEAAFYRRgAAgFWEEQAAYBVhBAggI2O7BAAIOoQRAABgFWEECCCHHLZLAICgQxgBAABWEUYAAIBVhBEAAGAVYQQIIFbTAID/CCMAAMAqwggQQKymAQD/EUYAAIBVhBEAAGAVYQQAAFhFGAEAAFYRRgAAgFWEEQAAYBVhBAAAWEUYAQAAVhFGAACAVYQRAABgFWEEAABYRRgBAABWEUYAAIBVhBEAAGAVYQQAAFhFGAEAAFYRRgAAgFUXFUYKCgqUmJio8PBwpaSkaNu2beftX1tbq7lz52rAgAFyOp269tprtW7duosqGAAAdC+X+btDUVGRsrKytGLFCqWkpGjZsmVKS0tTRUWFYmNjz+nf1NSkf/zHf1RsbKz+8Ic/aODAgTp06JCio6MDUT8AAAhyfoeRpUuXavbs2crMzJQkrVixQmvXrlVhYaF+9rOfndO/sLBQJ06c0NatW9WrVy9JUmJi4qVVDQAAug2/PqZpampSWVmZ3G73VwcICZHb7ZbH42lzn//+7/9Wamqq5s6dK5fLpZEjR2rx4sVqaWlp93kaGxtVX1/v8wCCgZGxXQIABB2/wsjx48fV0tIil8vl0+5yuVRVVdXmPvv379cf/vAHtbS0aN26dXr88ce1ZMkS/fznP2/3efLy8hQVFeV9JCQk+FMmAAAIIp2+mqa1tVWxsbF66aWXlJycrPT0dD366KNasWJFu/tkZ2errq7O+zh8+HBnlwkEhEMO2yUAQNDx65qRmJgYhYaGqrq62qe9urpacXFxbe4zYMAA9erVS6Ghod626667TlVVVWpqalJYWNg5+zidTjmdTn9KAwAAQcqvMyNhYWFKTk5WSUmJt621tVUlJSVKTU1tc58bb7xRe/fuVWtrq7ft448/1oABA9oMIgAAoGfx+2OarKwsrVy5Ur/5zW+0e/du3X///Tp16pR3dU1GRoays7O9/e+//36dOHFC8+fP18cff6y1a9dq8eLFmjt3buBGAQAAgpbfS3vT09N17Ngx5eTkqKqqSmPGjFFxcbH3otbKykqFhHyVcRISErR+/Xo9+OCDGj16tAYOHKj58+dr4cKFgRsF0EWwmgYA/OcwxnT5V8/6+npFRUWprq5OkZGRtssB2vXS5n1avG6PJOlg/lTL1QCAXR19/+a7aYAAYjUNAPiPMAIAAKwijAAAAKsIIwAAwCrCCAAAsIowAgAArCKMAAAAqwgjAADAKsIIAACwijACAACsIowAAACrCCMAAMAqwggAALCKMAIAAKwijAAAAKsII0AAGRnbJQBA0CGMAAAAqwgjQAA55LBdAgAEHcIIAACwijACAACsIowAAACrCCNAALGaBgD8RxgBAABWEUaAAGI1DQD4jzACAACsIowAAACrCCMAAMAqwggQQKymAQD/EUYAAIBVhBEggFhNAwD+I4wAAACrCCMAAMAqwggAALCKMAIAAKwijAAAAKsIIwAAwCrCCAAAsIowAgAArCKMAAAAqwgjAADAKsIIAACwijACAACsIowAAACrCCMAAMAqwggQQEbGdgkAEHQIIwAAwCrCCBBADjlslwAAQYcwAgAArCKMAAAAqwgjAADAKsIIEECspgEA/xFGAACAVYQRIIBYTQMA/ruoMFJQUKDExESFh4crJSVF27Zt69B+q1evlsPh0PTp0y/maQEAQDfkdxgpKipSVlaWcnNztWPHDiUlJSktLU01NTXn3e/gwYP66U9/qokTJ150sQAAoPvxO4wsXbpUs2fPVmZmpkaMGKEVK1YoIiJChYWF7e7T0tKimTNnatGiRRo6dOglFQwAALoXv8JIU1OTysrK5Ha7vzpASIjcbrc8Hk+7+z355JOKjY3VPffc06HnaWxsVH19vc8DCAaspgEA//kVRo4fP66Wlha5XC6fdpfLpaqqqjb32bJli15++WWtXLmyw8+Tl5enqKgo7yMhIcGfMgEAQBDp1NU0DQ0Nuvvuu7Vy5UrFxMR0eL/s7GzV1dV5H4cPH+7EKoHAYTUNAPjvMn86x8TEKDQ0VNXV1T7t1dXViouLO6f/vn37dPDgQU2bNs3b1tra+uUTX3aZKioqdNVVV52zn9PplNPp9Kc0AAAQpPw6MxIWFqbk5GSVlJR421pbW1VSUqLU1NRz+g8fPlwffPCBysvLvY877rhDt956q8rLy/n4BQAA+HdmRJKysrI0a9YsjR8/XhMmTNCyZct06tQpZWZmSpIyMjI0cOBA5eXlKTw8XCNHjvTZPzo6WpLOaQcAAD2T32EkPT1dx44dU05OjqqqqjRmzBgVFxd7L2qtrKxUSAg3dgWMMXI4uIYEAC7EYYzp8msR6+vrFRUVpbq6OkVGRtouB2jXys379fS63ZKkA3m3E0YA9Ggdff/mFAYAALCKMAIAAKwijAAAAKsII0AAff128F3/aiwA6BoIIwAAwCrCCBBA3A4eAPxHGAEAAFYRRgAAgFWEEQAAYBVhBAggn9U0FusAgGBCGAEAAFYRRoAAYjUNAPiPMAIAAKwijAAAAKsIIwAAwCrCCBBAvt9Nw3oaAOgIwggAALCKMAIEEKtpAMB/hBEAAGAVYQQAAFhFGAEAAFYRRoBOwloaAOgYwggAALCKMAIAAKwijAAAAKsIIwAAwCrCCAAAsIowAnQSvpoGADqGMAIAAKwijAAAAKsIIwAAwCrCCAAAsIowAgSQ+dpN4A03hAeADiGMAAAAqwgjQAA55LBdAgAEHcIIAACwijACAACsIowAAACrCCNAAPmspmExDQB0CGEEAABYRRgBAojVNADgP8IIAACwijACAACsIowAAACrCCNAAPF9NADgP8IIAACwijACBBCraQDAf4QRAABgFWEEAABYRRgBAABWEUaATsJ30wBAxxBGAACAVYQRAABg1UWFkYKCAiUmJio8PFwpKSnatm1bu31XrlypiRMnql+/furXr5/cbvd5+wMAgJ7F7zBSVFSkrKws5ebmaseOHUpKSlJaWppqamra7F9aWqoZM2Zo06ZN8ng8SkhI0KRJk3TkyJFLLh4AAAQ/v8PI0qVLNXv2bGVmZmrEiBFasWKFIiIiVFhY2Gb/V199VQ888IDGjBmj4cOH69e//rVaW1tVUlJyycUDXc3XbwfPreEBoGP8CiNNTU0qKyuT2+3+6gAhIXK73fJ4PB06xunTp9Xc3Kz+/fu326exsVH19fU+DwAA0D35FUaOHz+ulpYWuVwun3aXy6WqqqoOHWPhwoWKj4/3CTTflJeXp6ioKO8jISHBnzIBa7gdPAD471tdTZOfn6/Vq1frzTffVHh4eLv9srOzVVdX530cPnz4W6wSAAB8my7zp3NMTIxCQ0NVXV3t015dXa24uLjz7vv8888rPz9fGzdu1OjRo8/b1+l0yul0+lMaAAAIUn6dGQkLC1NycrLPxadnL0ZNTU1td79nn31WTz31lIqLizV+/PiLrxYAAHQ7fp0ZkaSsrCzNmjVL48eP14QJE7Rs2TKdOnVKmZmZkqSMjAwNHDhQeXl5kqRnnnlGOTk5eu2115SYmOi9tqRPnz7q06dPAIcC2OezmobFNADQIX6HkfT0dB07dkw5OTmqqqrSmDFjVFxc7L2otbKyUiEhX51wWb58uZqamvTDH/7Q5zi5ubl64oknLq16AAAQ9PwOI5I0b948zZs3r81tpaWlPr8fPHjwYp4CCEqspgEA//HdNAAAwCrCCAAAsIowAgAArCKMAAHk+900AICOIIwAAACrCCNAALGaBgD8RxgBAABWEUYAAIBVhBEAAGAVYQQIIN/vpmE9DQB0BGEEAABYRRgBAojVNADgP8IIAACwijACAACsIowAAACrCCNAJ2EtDQB0DGEEAABYRRgBAABWEUYAAIBVhBEAAGAVYQQIIN/bwVssBACCCGEEAABYRRgBAojbwQOA/wgjAADAKsIIAACwijACAACsIowAAfT11TTcDx4AOoYwAgAArCKMAAHEahoA8B9hBAAAWEUYAQAAVhFGAACAVYQRIIB8vpuG5TQA0CGEEQAAYBVhBAggVtMAgP8IIwAAwCrCCAAAsIowAgAArCKMAAHks5qGxTQA0CGEEQAAYBVhBAggVtMAgP8IIwAAwCrCCAAAsIowAgSQ7+3gAQAdQRgBAABWEUaAAOICVgDwH2EEAABYRRgBAABWEUYAAIBVhBEggHxvB896GgDoCMIIAACwijACBBCraQDAfxcVRgoKCpSYmKjw8HClpKRo27Zt5+3/+9//XsOHD1d4eLhGjRqldevWXVSxAACg+/E7jBQVFSkrK0u5ubnasWOHkpKSlJaWppqamjb7b926VTNmzNA999yjnTt3avr06Zo+fbp27dp1ycUDALqvhi+aufaqh/A7jCxdulSzZ89WZmamRowYoRUrVigiIkKFhYVt9v+P//gPTZ48WQsWLNB1112np556SuPGjdMvf/nLSy4eANA9fVzdoHFPbdBja/gf157gMn86NzU1qaysTNnZ2d62kJAQud1ueTyeNvfxeDzKysryaUtLS9OaNWvafZ7GxkY1NjZ6f6+vr/enzA57ecsB/fVvpzvl2OiZdh2p8/78/B8rFN4r1GI1QPBave2wmluMXn23UmGXcXnjt+HHNw5RQv8IK8/tVxg5fvy4Wlpa5HK5fNpdLpf27NnT5j5VVVVt9q+qqmr3efLy8rRo0SJ/Srsoa9//VDsqazv9edAzvb7tsO0SgG7hlT8dtF1CjzAtKT44wsi3JTs72+dsSn19vRISEgL+PD9IvlKpV10e8OOiZ9txqFYD+/WWK9JpuxQgqG355LiuHxilfhG9bJfSI7giw609t19hJCYmRqGhoaqurvZpr66uVlxcXJv7xMXF+dVfkpxOp5zOzn8hn5kyuNOfAwBwcRak2a4A3xa/PogLCwtTcnKySkpKvG2tra0qKSlRampqm/ukpqb69JekDRs2tNsfAAD0LH5/TJOVlaVZs2Zp/PjxmjBhgpYtW6ZTp04pMzNTkpSRkaGBAwcqLy9PkjR//nzdfPPNWrJkiaZOnarVq1dr+/bteumllwI7EgAAEJT8DiPp6ek6duyYcnJyVFVVpTFjxqi4uNh7kWplZaVCQr464XLDDTfotdde02OPPaZHHnlE11xzjdasWaORI0cGbhQAACBoOUwQ3FGmvr5eUVFRqqurU2RkpO1yAABAB3T0/ZvF2wAAwCrCCAAAsIowAgAArCKMAAAAqwgjAADAKsIIAACwijACAACsIowAAACrCCMAAMAqv28Hb8PZm8TW19dbrgQAAHTU2fftC93sPSjCSENDgyQpISHBciUAAMBfDQ0NioqKand7UHw3TWtrqz799FP17dtXDocjYMetr69XQkKCDh8+3G2/86a7j5HxBb/uPsbuPj6p+4+R8V08Y4waGhoUHx/v8yW63xQUZ0ZCQkJ05ZVXdtrxIyMju+U/sK/r7mNkfMGvu4+xu49P6v5jZHwX53xnRM7iAlYAAGAVYQQAAFjVo8OI0+lUbm6unE6n7VI6TXcfI+MLft19jN19fFL3HyPj63xBcQErAADovnr0mREAAGAfYQQAAFhFGAEAAFYRRgAAgFU9OowUFBQoMTFR4eHhSklJ0bZt22yX1CF5eXn6+7//e/Xt21exsbGaPn26KioqfPrccsstcjgcPo/77rvPp09lZaWmTp2qiIgIxcbGasGCBTpz5sy3OZQ2PfHEE+fUPnz4cO/2L774QnPnztXll1+uPn366Ac/+IGqq6t9jtFVxyZJiYmJ54zP4XBo7ty5koJz7jZv3qxp06YpPj5eDodDa9as8dlujFFOTo4GDBig3r17y+1265NPPvHpc+LECc2cOVORkZGKjo7WPffco5MnT/r0ef/99zVx4kSFh4crISFBzz77bGcPTdL5x9fc3KyFCxdq1KhR+ru/+zvFx8crIyNDn376qc8x2pr3/Px8nz62xiddeA5/9KMfnVP/5MmTffoE6xxKavNv0uFw6LnnnvP26cpz2JH3hUC9dpaWlmrcuHFyOp26+uqrtWrVqksfgOmhVq9ebcLCwkxhYaH58MMPzezZs010dLSprq62XdoFpaWlmVdeecXs2rXLlJeXm9tvv90MGjTInDx50tvn5ptvNrNnzzZHjx71Purq6rzbz5w5Y0aOHGncbrfZuXOnWbdunYmJiTHZ2dk2huQjNzfXXH/99T61Hzt2zLv9vvvuMwkJCaakpMRs377dfOc73zE33HCDd3tXHpsxxtTU1PiMbcOGDUaS2bRpkzEmOOdu3bp15tFHHzVvvPGGkWTefPNNn+35+fkmKirKrFmzxrz33nvmjjvuMEOGDDGff/65t8/kyZNNUlKS+fOf/2zeeecdc/XVV5sZM2Z4t9fV1RmXy2Vmzpxpdu3aZV5//XXTu3dv86tf/crq+Gpra43b7TZFRUVmz549xuPxmAkTJpjk5GSfYwwePNg8+eSTPvP69b9Zm+O70BiNMWbWrFlm8uTJPvWfOHHCp0+wzqExxmdcR48eNYWFhcbhcJh9+/Z5+3TlOezI+0IgXjv3799vIiIiTFZWlvnoo4/ML37xCxMaGmqKi4svqf4eG0YmTJhg5s6d6/29paXFxMfHm7y8PItVXZyamhojyfzf//2ft+3mm2828+fPb3efdevWmZCQEFNVVeVtW758uYmMjDSNjY2dWe4F5ebmmqSkpDa31dbWml69epnf//733rbdu3cbScbj8RhjuvbY2jJ//nxz1VVXmdbWVmNMcM+dMeacF/rW1lYTFxdnnnvuOW9bbW2tcTqd5vXXXzfGGPPRRx8ZSeYvf/mLt8///u//GofDYY4cOWKMMebFF180/fr18xnjwoULzbBhwzp5RL7aeiP7pm3bthlJ5tChQ962wYMHmxdeeKHdfbrK+Ixpe4yzZs0yd955Z7v7dLc5vPPOO80//MM/+LQF0xx+830hUK+dDz/8sLn++ut9nis9Pd2kpaVdUr098mOapqYmlZWVye12e9tCQkLkdrvl8XgsVnZx6urqJEn9+/f3aX/11VcVExOjkSNHKjs7W6dPn/Zu83g8GjVqlFwul7ctLS1N9fX1+vDDD7+dws/jk08+UXx8vIYOHaqZM2eqsrJSklRWVqbm5mafuRs+fLgGDRrknbuuPrava2pq0m9/+1v9+Mc/9vkSyGCeu286cOCAqqqqfOYsKipKKSkpPnMWHR2t8ePHe/u43W6FhITo3Xff9fa56aabFBYW5u2TlpamiooK/e1vf/uWRtMxdXV1cjgcio6O9mnPz8/X5ZdfrrFjx+q5557zOf0dDOMrLS1VbGyshg0bpvvvv1+fffaZd1t3msPq6mqtXbtW99xzzznbgmUOv/m+EKjXTo/H43OMs30u9b0zKL4oL9COHz+ulpYWn//gkuRyubRnzx5LVV2c1tZW/eQnP9GNN96okSNHetvvuusuDR48WPHx8Xr//fe1cOFCVVRU6I033pAkVVVVtTn+s9tsSklJ0apVqzRs2DAdPXpUixYt0sSJE7Vr1y5VVVUpLCzsnBd5l8vlrbsrj+2b1qxZo9raWv3oRz/ytgXz3LXlbE1t1fz1OYuNjfXZftlll6l///4+fYYMGXLOMc5u69evX6fU768vvvhCCxcu1IwZM3y+dOzf//3fNW7cOPXv319bt25Vdna2jh49qqVLl0rq+uObPHmyvv/972vIkCHat2+fHnnkEU2ZMkUej0ehoaHdag5/85vfqG/fvvr+97/v0x4sc9jW+0KgXjvb61NfX6/PP/9cvXv3vqiae2QY6U7mzp2rXbt2acuWLT7tc+bM8f48atQoDRgwQLfddpv27dunq6666tsu0y9Tpkzx/jx69GilpKRo8ODB+t3vfnfR/9C7qpdffllTpkxRfHy8ty2Y566na25u1j//8z/LGKPly5f7bMvKyvL+PHr0aIWFhenf/u3flJeXFxS3Gf+Xf/kX78+jRo3S6NGjddVVV6m0tFS33XabxcoCr7CwUDNnzlR4eLhPe7DMYXvvC11Zj/yYJiYmRqGhoedcRVxdXa24uDhLVflv3rx5+p//+R9t2rRJV1555Xn7pqSkSJL27t0rSYqLi2tz/Ge3dSXR0dG69tprtXfvXsXFxampqUm1tbU+fb4+d8EytkOHDmnjxo269957z9svmOdO+qqm8/29xcXFqaamxmf7mTNndOLEiaCZ17NB5NChQ9qwYcMFv4o9JSVFZ86c0cGDByV1/fF909ChQxUTE+Pz7zLY51CS3nnnHVVUVFzw71LqmnPY3vtCoF472+sTGRl5Sf+z2CPDSFhYmJKTk1VSUuJta21tVUlJiVJTUy1W1jHGGM2bN09vvvmm3n777XNOC7alvLxckjRgwABJUmpqqj744AOfF4+zL6AjRozolLov1smTJ7Vv3z4NGDBAycnJ6tWrl8/cVVRUqLKy0jt3wTK2V155RbGxsZo6dep5+wXz3EnSkCFDFBcX5zNn9fX1evfdd33mrLa2VmVlZd4+b7/9tlpbW71hLDU1VZs3b1Zzc7O3z4YNGzRs2DDrp/fPBpFPPvlEGzdu1OWXX37BfcrLyxUSEuL9aKMrj68tf/3rX/XZZ5/5/LsM5jk86+WXX1ZycrKSkpIu2LcrzeGF3hcC9dqZmprqc4yzfS75vfOSLn8NYqtXrzZOp9OsWrXKfPTRR2bOnDkmOjra5yrirur+++83UVFRprS01GeJ2enTp40xxuzdu9c8+eSTZvv27ebAgQPmrbfeMkOHDjU33XST9xhnl3BNmjTJlJeXm+LiYnPFFVd0ieWvDz30kCktLTUHDhwwf/rTn4zb7TYxMTGmpqbGGPPl8rRBgwaZt99+22zfvt2kpqaa1NRU7/5deWxntbS0mEGDBpmFCxf6tAfr3DU0NJidO3eanTt3Gklm6dKlZufOnd7VJPn5+SY6Otq89dZb5v333zd33nlnm0t7x44da959912zZcsWc8011/gsC62trTUul8vcfffdZteuXWb16tUmIiLiW1k2eb7xNTU1mTvuuMNceeWVpry83Odv8uwKhK1bt5oXXnjBlJeXm3379pnf/va35oorrjAZGRldYnwXGmNDQ4P56U9/ajwejzlw4IDZuHGjGTdunLnmmmvMF1984T1GsM7hWXV1dSYiIsIsX778nP27+hxe6H3BmMC8dp5d2rtgwQKze/duU1BQwNLeS/WLX/zCDBo0yISFhZkJEyaYP//5z7ZL6hBJbT5eeeUVY4wxlZWV5qabbjL9+/c3TqfTXH311WbBggU+96owxpiDBw+aKVOmmN69e5uYmBjz0EMPmebmZgsj8pWenm4GDBhgwsLCzMCBA016errZu3evd/vnn39uHnjgAdOvXz8TERFhvve975mjR4/6HKOrju2s9evXG0mmoqLCpz1Y527Tpk1t/pucNWuWMebL5b2PP/64cblcxul0mttuu+2csX/22WdmxowZpk+fPiYyMtJkZmaahoYGnz7vvfee+e53v2ucTqcZOHCgyc/Ptz6+AwcOtPs3efbeMWVlZSYlJcVERUWZ8PBwc91115nFixf7vJHbHN+Fxnj69GkzadIkc8UVV5hevXqZwYMHm9mzZ5/zP2/BOodn/epXvzK9e/c2tbW15+zf1efwQu8LxgTutXPTpk1mzJgxJiwszAwdOtTnOS6W4/8PAgAAwIoeec0IAADoOggjAADAKsIIAACwijACAACsIowAAACrCCMAAMAqwggAALCKMAIAAKwijAAAAKsIIwAAwCrCCAAAsIowAgAArPp/xJubLXY/5FQAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "#plt.semilogy(loss_history)\n",
    "plt.plot(loss_history)\n",
    "loss_history[-1] # -8.679537108440963"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "110744c6-c521-4c91-a657-9d5c9b14ef44",
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.semilogy(scale_history)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8ca853a8-29e3-487a-a882-dbe5ef6c8f2f",
   "metadata": {},
   "outputs": [],
   "source": [
    "J = model.jacobian(rs)\n",
    "metric_tensor = model.metric_tensor(J)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ff8e0ecf-e259-4d4a-8f79-90a349df257e",
   "metadata": {},
   "outputs": [],
   "source": [
    "fig,axs = plt.subplots(ncols=3,figsize=(9,3))\n",
    "axs[0].scatter(*rs.T,c=metric_tensor[...,0,0].detach())\n",
    "axs[1].scatter(*rs.T,c=metric_tensor[...,1,1].detach())\n",
    "axs[2].scatter(*rs.T,c=metric_tensor[...,1,0].detach())\n",
    "#plt.colorbar()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3c3f93fd-f3ba-49ce-b6b1-0e829904f9b2",
   "metadata": {},
   "outputs": [],
   "source": [
    "N = len(phases_history)\n",
    "fig, ax = plt.subplots()\n",
    "[ax.scatter(*phases_history[i].T, c='grey',alpha=(i+0.2)/(N+0.2)) for i in range(0,N,10)]\n",
    "ax.scatter(*phases_history[-1].T, c='red')\n",
    "ax.scatter(*start_phases.T,s=10)\n",
    "ax.axis('equal')\n",
    "model.unit_cell.plot(fig,ax)\n",
    "#model.unit_cell.plot(fig,ax,center=np.mean(phases_history[-1],axis=0))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "efe3ff21-5279-4eb0-b8ab-e3e031c02551",
   "metadata": {},
   "outputs": [],
   "source": [
    "torch.save(model,f'./saved-models/{type(model).__name__}-{model.ncells}-dirac')\n",
    "f'./saved-models/{type(model).__name__}-{model.ncells}-dirac'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "25370ce8-c4a7-4b92-9ce3-13db1c343b31",
   "metadata": {},
   "outputs": [],
   "source": [
    "# dropout?\n",
    "# figure out the geodesic between phases\n",
    "# - finds hexagons"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python [conda env:ml]",
   "language": "python",
   "name": "conda-env-ml-py"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.9.12"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}