notebook.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import pandas as pd\n",
    "import matplotlib.pyplot as plt\n",
    "import scanpy as sc\n",
    "\n",
    "# from chiral import CHIRAL\n",
    "# from helper_fn import optimal_shift\n",
    "\n",
    "from pyCHIRAL import CHIRAL, optimal_shift, ccg"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/Users/salati/miniforge3/envs/numpyro/lib/python3.12/site-packages/anndata/_core/anndata.py:1908: UserWarning: Variable names are not unique. To make them unique, call `.var_names_make_unique`.\n",
      "  utils.warn_names_duplicates(\"var\")\n"
     ]
    }
   ],
   "source": [
    "path = \"/Users/salati/Documents/CODE/github/scCircadianMeta/data/BULK/zhang.h5ad\"\n",
    "adata = sc.read_h5ad(path)\n",
    "\n",
    "organ = \"liver\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "# log - transform the data\n",
    "cc = np.median(adata.X.sum(axis=1))\n",
    "eps = 1/cc\n",
    "adata.layers[\"s_log\"] = np.log(adata.layers[\"s_norm\"] + eps)\n",
    "\n",
    "adata2 = adata[adata.obs.organ == organ]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array(['Arntl', 'Npas2', 'Cry1', 'Cry2', 'Per1', 'Per2', 'Nr1d1', 'Nr1d2',\n",
       "       'Tef', 'Dbp', 'Ciart', 'Per3', 'Bmal1'], dtype='<U5')"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "ccg"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/Users/salati/miniforge3/envs/numpyro/lib/python3.12/site-packages/anndata/_core/anndata.py:1908: UserWarning: Variable names are not unique. To make them unique, call `.var_names_make_unique`.\n",
      "  utils.warn_names_duplicates(\"var\")\n",
      "  0%0.00s - Debugger warning: It seems that frozen modules are being used, which may\n",
      "0.00s - make the debugger miss breakpoints. Please pass -Xfrozen_modules=off\n",
      "0.00s - to python to disable frozen modules.\n",
      "0.00s - Note: Debugging will proceed. Set PYDEVD_DISABLE_FILE_VALIDATION=1 to disable this validation.\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1[{\"variableName\": \"ID_TO_MEANING\", \"type\": \"dictionary\", \"supportedEngines\": [\"pandas\"], \"isLocalVariable\": true}]\n",
      "\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "ERROR:tornado.general:SEND Error: Host unreachable\n",
      "/Users/salati/Documents/CODE/github/scCircadianMeta/rhythm_module/pyCHIRAL/pyCHIRAL/em.py:304: FutureWarning: The behavior of DataFrame concatenation with empty or all-NA entries is deprecated. In a future version, this will no longer exclude empty or all-NA columns when determining the result dtypes. To retain the old behavior, exclude the relevant entries before the concat operation.\n",
      "  Q_hist = pd.concat([Q_hist, Q_temp])  # Append the new iteration results\n",
      "  2%"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Algorithm has converged\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "\n"
     ]
    }
   ],
   "source": [
    "res = CHIRAL(E=adata2, clockgenes=ccg ,layer=\"s_log\", iterations=500)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.10404895698537764"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "w =  2 * np.pi / 24\n",
    "true_phase = (adata2.obs.ZT.values % 24) * w\n",
    "phi = res[\"phi\"]\n",
    "\n",
    "phase_aligned, mad = optimal_shift(phi, true_phase)\n",
    "mad"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.lines.AxLine at 0x16dbb7ad0>"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAioAAAGwCAYAAACHJU4LAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjguMCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy81sbWrAAAACXBIWXMAAA9hAAAPYQGoP6dpAABB9ElEQVR4nO3deXhU1f3H8fcMISAwAWUH2bdiQdYqIIKIFFwqoJWtIoigUhXEpRA3rApIVVaxWEVUQBRtVfrDBUVAClKBAKKgVYMIEcKaBQgkIef3xxmjkcUMTHLuzHxez/N9ntybm5lvYu18POfec3yAQURERMSD/K4bEBERETkZBRURERHxLAUVERER8SwFFREREfEsBRURERHxLAUVERER8SwFFREREfGsONcNnKkaNWqQmZnpug0REREJQSAQ4IcffvjV6yI6qNSoUYOUlBTXbYiIiMhpqFmz5q+GlYgOKj+OpNSsWVOjKiIiIh7lA0YAY7HBIzkQoH5KSqE+uyM6qPwoMzNTQUVERMSDKgIvA1cEj+cD9wCFnQ/RzbQiIiJSJDoCG7AhJQsYBgwAQhlaUFARERGRsPIB9wHLgHOBL4ELgedP47WiYupHREREvKEKMAf4ffD4ZeDPwKHTfD0FFREREQmLS4BXgOrAYeA24MUzfE1N/YiIiMgZ8QMPAR9iQ8oXQFvOPKSARlRERETkDFQD5gGXBo9nAXdgb54NBwUVEREROS2XAXOBqsBB4FZsaAknTf2IiIhISEoAjwLvY0PKRqAN4Q8poBEVERERCUFN7A2znYLHM4FRwJEiej8FFRERESmUy7GPG1cCMrALuC0o4vfU1I+IiIicUhwwEXgHG1LWAa0p+pDy43uLiIiInFAt4FWgQ/B4Onavnuxien8FFRERETmhP2DXQjkHSANuAv5VzD1o6kdEREQKKAk8BSzEhpRPgVYUf0gBBRURERH5mbrAf4C7gseTsLsgf+eoH039iIiICAC9gReACsB+YDDwb4f9gEZUREREYl4pYBp2aqcCsApoifuQAgoqIiIiMa0BNpjcETyeCHQGtjvrqCBN/YiIiMSoPsBzQAKwF7gBeNdpR8fTiIqIiEiMKQ38HXgNG1I+xk71eC2kgIKKiIhITGkMrMbudJwHPAZcCqS4bOoUNPUjIiISI/6E3USwHJAKXA986LSjX6cRFRERkSh3FvA8MBcbUj7CTvV4PaSAgoqIiEhUawqswS5/nweMBboBu1w2FQJN/YiIiESpwcAMoAywExgALHPYz+lwPqJSo0YN5syZw969ezl06BDr16+ndevWrtsSERGJWGWBl4DZ2JCyGDvVs8xdS6fN6YhKhQoVWLlyJUuXLuXyyy9n9+7dNGjQgLS0NJdtiYiIRKzmwALgN8Ax4EHgccC4bOoMOA0qo0ePZvv27QwZMiT/3LZt2xx2JCIiErmGAVOxN8/uAPpjNxiMZE6nfq6++mrWrl3LggULSE1NJSkpiaFDh570+vj4eAKBQIESERGJdQHgFeAf2JDyDnaqJ9JDyo+Mq8rKyjJZWVlm3LhxpmXLlubmm282hw8fNgMHDjzh9WPHjjUnEggEnP0OKpVKpVK5rJZg/gfGgMkBcw8Ynwf6OlUFAoFQPr/dNXr06FGzcuXKAuemTp1qVq1adcLr4+PjTSAQyK8aNWooqKhUKpUqZuvPYI5gQ8p3YNp5oKfCVChBxenUz86dO9m8eXOBc1u2bKF27donvD47O5vMzMwCJSIiEmvKY2+YnQGUAt4GWmGXxo82ToPKypUradKkSYFzjRs31g21IiIiJ9EWSAKuA7KBO4FewAF3LRU5Z0M/bdu2NdnZ2SYxMdE0aNDA9O/f3xw8eNAMGDAg7ENHKpVKpVJFeo0EcxQ71ZMMpq0Hejqdiph7VABz5ZVXms8++8xkZWWZzZs3m6FDhxbVL6pSqVQqVUTW2WDexAYUA+YNMOU90NfpViif377gFxEpEAiQkZFBQkKC7lcREZGI5PP7qd+6BQmVK5GxZy/JSRsxeXn5378QeA2oAxwF7gKecdNq2ITy+a29fkRERBxp3rUzvcaMokK1qvnn0nal8tbjk/l8yXLuBsYDJYFvgD7AejetOqOgIiIi4kDzrp0ZNGkCv5zYKF+lMrc+OJpmGzbRad9+AF4FbgZice7A+aaEIiIiscbn99NrzCjA4PMX/Cg+d/1nDOwziE779nMEG1D6E5shBTSiIiIiUuzqt25RYLoHgLw8Lpj1Mh1mPI//2DH2163NwDJleGfzl26a9AiNqIiIiBSzhMqVChyftW8/1wy/i47TnsV/7Bibr+rBvNdms7NOLUcdeodGVERERIpZxp69+V+fuyaJK0aPpdyeveSULsVH993DF72uBJ+vwHWxSkFFRESkmCUnbSQ95Qe6v/0O7Z6djT8vj70N6rHoycfY17A+Ji+PtF2pJCdtdN2qcwoqIiIixaxKXh5X9L+JCw+kAfB5ryv5KPFucsucFVxDxcfbE6cUWE8lVimoiIiIFKOuwDyg6oE0skrGsXTsGL7seWX+940xLHtpHpuWLHfWo5foZloREZFiUAJ4FFgMVAX+V64sr77+Mlv+cHmB63w+H10GX0/zrp0ddOk9CioiIiJFrAbwEfAA9oP3WeDVV19gf706x62jYo8NPUffedz3YpH+AiIiIkWoB7AB6ARkAP2AJ9q2olyd2icNIj6/n7OrV6N+6xbF1qdXKaiIiIgUgTjgceBdoDKQBLTBbjD4y3VUTqaw10Uz3UwrIiISZrWA+cBFweOngXuwux8DhV4fReuoaERFREQkrK7CTvVcBKQB1wJ38FNIAbuOStqu1JM+fmzy8jiwc5fWUUFBRUREJCxKAk8C/wbOAT4FWgP/OsG1Ji+Ptx6fDPiOCytaR6UgBRUREZEzVBdYAdwdPJ4MdAS2nuJnNi1Zzkt3JZK+e0+B82mpu3nprkStoxLkA4zrJk5XIBAgIyODhIQEMjNjdQNsERFxqTfwAlABOAAMBhaG8PM+v5/6rVuQULkSGXv2kpy0MepHUkL5/NbNtCIiIqchHjvVc0fw+BPso8ffh/g6Ji+Pb9euD2drUUVTPyIiIiFqAKzip5AyEbtOSqghRX6dRlRERERCcB3wPJAA7AVuwK6VIkVDIyoiIiKFUBr4O7AAG1JWAC1RSClqCioiIiK/ojGwGrgVyAPGAV2AFJdNxQhN/YiIiJzCAOwmguWA3cD1wAdOO4otGlERERE5gbOA54B52JCyFGiBQkpxU1ARERH5habYlWWHYqd6HgYuA3Y57ClWaepHRETkZwYBM4CywE7gT9jRFHFDIyoiIiLYYPJisMpip3haopDimoKKiIjEvGbAGuxoyjHgfqA79uZZcUtTPyIiEtOGAtOwN8+mAP2xa6SINyioiIhITApgHzvuHzx+BzuistdZR3IimvoREZGY0xJYhw0pucBfgKtQSPEijaiIiEhMGQ5MBkphNxHsh935WLxJQUVERGJCAnYzweuCxwuBwcABVw1JoWjqR0REol4bYD02pGQDo4CeKKREAo2oiIhIVBsBPAHEA1uBvthHkSUyKKiIiEhUqgC8APQOHv8TuAlId9WQnBZN/YiISNS5EDvV0xs4CtwO/BGFlEikoCIiIlHDB9yNXbCtLvAN0AG7d49EJk39iIhIVKiI3afnquDxq8DNQKarhiQsNKIiIiIR7yLsVM9VwBHgFuxibgopkU9BRUREIpYPGAMsA2oBX2HvT/mHw54kvDT1IyIiEakyMAe7yzHBr4cDh5x1JEXB6YjK2LFjMcYUqJ07d7psSUREIkBnYAM2pBwGhgA3oJASjZyPqHz++edcdtll+cfHjh1z2I2IiHiZH7gfGAuUADZjV5vd7LIpKVLOg0pubi6pqamu2xAREY+rCswDugaPXwDuwI6oSPRyfjNto0aNSElJITk5mfnz51OvXr2TXhsfH08gEChQIiIS/bpip3q6Yqd3BmJXmVVIiX5Og8p///tfbrjhBrp3786wYcOoVq0aq1at4pxzzjnh9YmJiWRkZORXSkpKMXcsIiLFqQTwCLAYqAZ8ht1gcK7LpqTYGa9UmTJlzM6dO82oUaNO+P34+HgTCATyq0aNGsYYYwKBgPPeVSqVShXeqgFmGRgTrGfBlPZAX6ozr0AgUOjPb+f3qPzc4cOH2bRpE40aNTrh97Ozs8nOzi7mrkREpLh1xz5uXBm7aNvN2JVmJfY4v0fl5+Lj42natKkeURYRiVFxwATgPWxIWQ+0RiElljkNKk888QSdOnWibt26XHDBBbzxxhskJCTw0ksvuWxLREQcOBe7wuyY4PHTQHvsxoISu5xO/Zx77rnMnz+fSpUqsWfPHlavXk27du34/vvvXbYlIiLF7CrshoIVgXTsEz3/dNmQeIbToNK/f3+Xby8iIo6VxE713B08XgP0BbY660i8xlM304qISOyoA7yG3UQQYDJ22kePTMjPKaiIiEix64VdWfZs4AAwGFjosB/xLk899SMiItEtHpgKvIkNKZ8ALVFIkZNTUBERkWJRH1gJjAgePwF0AvT4hJyKpn5ERKTIXQc8DyQAe4FBwDtOO5JIoREVEREpMqWAZ4AF2JCyAjvVo5AihaWgIiIiRaIRsBoYHjweD3QBtJ2shEJTPyIiEnYDgGeBcsBuYCB2B2SRUGlERUREwuYs4DlgHjakLMVO9SikyOlSUBERkbD4DfBfYCiQB/wVuAzQNrNyJjT1IyIiZ+wG7E2zZYFd2KmfpU47kmihERURETltZbCbCb6EDSkfAC1QSJHwUVAREZHT0gxYi10T5RjwANADe/OsSLho6kdEREJ2EzAde/NsCnaq52OnHUm0UlAREZFCK4d97HhA8Phd7P0pe511JNFOUz8iIlIoLYB12JCSC4wGrkQhRYqWRlRERORX3QpMBkpjNxHsh935WKSoKaiIiMhJJWAXcOsTPF4I3Ajsd9aRxBpN/YiIyAm1AZKwISUHuAvoiUKKFC+NqIiIyHHuAJ4E4oHvgL7Apy4bkpiloCIiIvkqAC8AvYPH/8I+ipzmqB8RTf2IiAgAFwDrsSHlKHZU5VoUUsQtBRURkRjnw95/8h+gLvAt0AF42mFPIj/S1I+ISAw7B7tPz1XB49eAm4EMZx2JFKQRFRGRGNUB2IANKUewa6X0QyFFvEVBRUQkxviwq8ouB2oBXwEXYpfGF/EaTf2IiMSQysDL2F2OAeYCw4GDzjoSOTUFFRGRGNEJmA/UAA5jn+p5wWlHIr9OUz8iIlHODzwAfIQNKZuxjyIrpEgk0IiKiEgUq4qd3rkseDwbuB07oiISCRRURESi1KXAPKAacAh7L8ocpx2JhE5TPyIiUaYE8FfgA2xI2QS0RSFFIpNGVEREokh17A2znYPH/wBGYtdJEYlECioiIlGiO3bUpDKQCdyCDS0ikUxTPyIiEa4EMB54DxtSNgBtUEiR6KARFRGRCHYuNpB0DB7PAO7G7n4sEg0UVEREItSV2A0FKwLpwFDgDacdiYSfpn5ERCJMSeAJ4P+wIWUt0BqFFIlOGlEREYkgdYBXgXbB4ynYDQazXTUkUsQUVEREIkRP7MqyZwMHgBuBt512JFL0NPUjIuJx8diRk7ewIWU10AqFFIkNCioiIh5WH1iJXbQN4EnsLsjbnHUkUrw8E1TGjBmDMYbJkye7bkVExBP+CCRhl7/fB1wF3AvkuGxKpJh5Iqi0bduWm2++mY0bN7puRUTEuVLY9VBeB8oD/wFaAosc9iTiivOgUrZsWebNm8ewYcM4cOCA63ZERJxqhL0H5c/B4/FAF2CHs45E3DqtoNKxY0fmzJnDqlWrqFGjBgDXX389F110UcivNWPGDBYtWsSSJUt+9dr4+HgCgUCBEhGJFv2BddjRk93YvXvuB3Id9iTiWshB5ZprruH9998nKyuLVq1aUapUKQACgQD33XdfSK/Vt29fWrduTWJiYqGuT0xMJCMjI79SUlJCbV9ExHPOwu5y/AoQAJZhw8pidy2JeIoJpZKSkszAgQMNYDIyMky9evUMYFq0aGF27txZ6Nc599xzza5du8z555+ff27p0qVm8uTJJ/2Z+Ph4EwgE8qtGjRrGGGMCgUBIv4NKpVJ5pX4D5jMwBswxMA+DKeGBvlSqoqxAIBDK53doL37o0CFTp04dAwWDSr169UxWVlahX6dnz57GGGNycnLyyxhjjh07ZnJycozf7w/3L6pSqVSeqhvAHMSGlJ1gLvVATypVcVQon98hr0y7c+dOGjZsyLZtBZ/i79ixI8nJyYV+nSVLltCsWbMC52bPns2XX37JxIkTycvLC7U1EZGIUAb7VM/g4PGHwPVAqquGRDws5KDy7LPPMnXqVIYMGYIxhho1atC+fXuefPJJHnnkkUK/zsGDB/niiy8KnDt06BD79u077ryISLT4LbAAOA84BjyMfbJH/2kmcmIhB5UnnniC8uXLs3TpUkqXLs3HH3/M0aNHefLJJ5kxY0ZR9CgiEhVuAqZjb55NAQYAHzvtSMT7fNg5oJCdddZZnHfeefj9fjZv3syhQ4fC3NqvCwQCZGRkkJCQQGZmZrG/v4hIYZQDZgJ/Ch6/BwwE9jrrSMStUD6/T3vBt6ysLNatW8eXX37JZZddxm9+85vTfSkRkajVAliLDSm5wBjgChRSRAor5KDy2muvcdtttwFQunRp1qxZw4IFC/jss8+45pprwt6giEikuhW7ymwTYDvQGZjIaQ5ji8SokINKp06dWLFiBQC9e/fG7/dToUIFRowYwQMPPBD2BkVEIk0C8Crwd6A08G/sAm6rHPYkEqlCDirly5dn//79APTo0YN//vOfZGVlsWjRIho1ahT2BkVEIklr7I7HfbG7HN8FXA3sd9mUSAQLOahs376d9u3bU6ZMGXr06MHixXaR57PPPpsjR46EvUERkUhxO3bUpAHwHdARmFzE7+nz+2nQthWtLu9Gg7at8Pmd7zUrElYhP548ZcoU5s2bx8GDB9m2bRvLli0D7JTQpk2bwt2fiIjnVQBmAT/epfcmMARIK+L3bd61M73GjKJCtar559J2pfLW45PZtGR5Eb+7SPE4rceTW7duTe3atfnggw/yH0u+4oorSEtLY9Wq4puF1ePJIuLaBdj7UeoB2cA92LVSilrzrp0ZNGkCYAqMopi8PMDHS3clKqyIZ4Xy+X3a66h4gYKKiLh0F/A4UBL4FntfyrpieF+f388D7/+L8lUqn3Cqx+TlkZa6m3E9rg0GFxFvCeXzO+SpH4CaNWty9dVXU7t2beLj4wt87+677z6dlxQRiRjnAC8CfwgeLwCGARnF9P71W7coMN3zSz6/n7OrV6N+6xZ8u3Z9MXUlUjRCDiqXXnopCxcuZOvWrTRp0oTPP/+cunXr4vP5SEpKKooeRUQ8owMwH6gNHAFGYVedLU4JlSuF9ToRLwv59vAJEybw1FNP0bx5c44cOcK1115LrVq1WL58Oa+//npR9Cgi4pwPGA0sx4aU/wHtKP6QApCxp3Dr2hb2OhEvCzmoNG3alJdeegmA3NxczjrrLA4dOsRDDz3E6NGjw96giIhrlYBF2PtR4oB5QBtgo6N+kpM2krYr9aT3n5i8PA7s3EVykqsORcIn5KBy6NAhSpUqBcAPP/xAgwYN8r9XqZKGGUUkulwMbAAuB7KwOyBfDxx02JPJy+OtxycDvuPCyo9P/bw9cYpupJWoEHJQWb16NRdddBEAixYt4qmnnuK+++7jhRdeYPXq1WFvUETEBT9wP7AUqAlsAX4HvOCyqZ/ZtGQ5L92VSPruPQXOp6Xu1qPJElVCfjy5Xr16lCtXjk2bNnHWWWfx5JNP0rFjR7755htGjRrF999/X0StHk+PJ4tIUaiCnd65LHj8InAbcNhVQ6fgj4vjor7XUKn2uez9fgcrX/sXebm5rtsSOSWtoyIicpouxYaUasAh4M/Ay047OjmtTCuRqliCSsmSJalSpQr+Xyw2tH379tN5udOioCIi4eIHHgIeDH69CbuA2xaXTZ2CVqaVSFakC741atSIWbNm0aFDhwLnfT4fxhji4k5rDTkREWeqY0dRugSPnwNGYm+e9SKf30+vMaP4ZUj58XsmL4+eo+/k86UrdEOtRLyQU8Xs2bPJzc3lqquuYufOnRgTsTNHIiL8HpiDvS8lE7gFu6Cbl2llWoklIQeVli1b0qZNG7766qui6EdEpFiUAB4FEoPHG4A+wNeuGgqBVqaVWBLy48mbN2/WeikiEtHOBZbxU0h5BrvKbCSEFNDKtBJbChVUAoFAfo0ePZq//e1vdO7cmXPOOafA9wKBQFH3KyJyRq7Ajp50xG4i2Af76PFRhz2FSivTSiwp1NRPWlpagXtRfD4fS5YsKXCNbqYVES+LA8YD9waP12Kf6kl21tHp+3Fl2kGTJmDy8k741I9WppVoUahU0aVLl1+/SETEo2oDr2GndwCmAn8Bsp11dOZ+XJn2uHVUUnfz9sQpejRZooYWfBORqNYTmA2cDRwAhgBvuWwozHx+P/VbtyChciUy9uwlOWmjRlLE84p0HRWAChUqcNNNN9G0aVOMMWzZsoXZs2dz4MCB02pYRCTcSgJ/A+4MHv8XO9WzzVVDRcTk5ekRZIlqIT/106lTJ7777jtGjBjB2WefzTnnnMOIESPYunUrnTp1KooeRURCUg9YyU8h5UnsLsjRFlJEYoUJpTZt2mSeffZZ4/f788/5/X4zc+ZMs2nTppBe60wrEAgYY4wJBALF+r4qlcq7dS2YNDAGzF4wV3mgJ5VKVbBC/PwO7cUPHz5sGjdufNz5xo0bm8OHD3v5F1WpVFFcpcA8jQ0oBsx/wJzrgb5UKtXxFcrnd8hTP0lJSTRt2vS4802bNmXDhg2hvpyIyBlrCHyCXQ8FYAJwCbDDVUMiEjYh30w7bdo0pk6dSsOGDVm9ejUA7dq147bbbmPMmDE0b948/9pNmzaFr1MRkRPoB/wDCAB7gIHA+047EpFwCvnx5GPHjp3y+8aYYlv8TY8ni8Su0tj1UG4OHi8HBgA/OOtIRAqrSB9Prlev3mk3JiISDk2A14HmQB7wGPAIcOr/jBKRSBRyUPn++++Log8RkUIZCPwdKAvsAq4HlpzyJ0QkkoV8M62IiAtlgBeAl7EhZQnQEoUUkWinoCIinnce8ClwI3Z650Hg90Cqy6ZEpFhoq2MR8bQhwHTsiMoP2Btmtd2eSOxQUBERTyqHvRfl+uDx+9j7U/Y460hEXNDUj4h4zvnAWmxIyQXGAJejkCISiwo1orJ//36MKdxyKxUrVjyjhkQktt0CTMGuk7Id6I/dYFBEYlOhgsqdd96Z/3XFihV54IEHeP/99/nkk08AaN++Pd27d+fRRx8tkiZFJPoFgOeAvsHj/wMGA/tcNSQinhHSRkJvvPGGue222447f9ttt5k333zTs5saqVQq71ZrMF9jNxPMBnM3GJ8H+lKpVEVTRbp7cmZmpmnQoMFx5xs2bGgyMzO9/IuqVCoP1u1gjmBDylYwF3qgJ5VKVbRVpLsn79u3j969ex93vlevXuzbF9og7a233srGjRtJT08nPT2dVatW0aNHj1BbEpEIVB54A/vocSngTaAV8F+XTYmI54T8ePLYsWOZNWsWl1xySf49Ku3ataNHjx4MHTo0pNfasWMHY8aM4ZtvvgFg0KBBvP3227Rq1YrNmzeH2pqIRIjfAa8B9YBs4B5sYBEROZGQh2wuuOACM3fuXLNu3TqTlJRk5s6day644IKwDAft27fPDBky5ITfi4+PN4FAIL9q1KihqR+VKsJqFPY+FAPmGzBtPNCTSqUq3irSe1SKqvx+v+nbt685cuSIadq06QmvGTt2rDkRBRWVyvt1Npi3sQHFgFkAJsEDfalUquKvIg8q9evXN48++qiZN2+eqVy5sgFM9+7dzXnnnRfyazVr1sxkZmaanJwcc+DAAXP55Zef9FqNqKhUkVntwWzDBpQsMLd6oCeVSuWuijSodOrUyRw6dMgsXrzYHDlyxNSrV88A5t577zWvv/56yM2WLFnSNGjQwLRp08aMHz/e7N69+6QjKmf4i6pUqmIuH5i/gMnBhpSvwLTwQF8qlcptFWlQWbVqlRk1apQBTEZGRn5Qadu2rdmxY8cZN//BBx+YmTNnFsUvqlKpirEqgVnET1M988CU80BfKpXKfRXp48nNmzfnzTffPO78nj17wrJ8vs/no1SpUmf8OiLizsXABuAKIAsYBvwJOOiwJxGJTCE/npyWlkb16tX57rvvCpxv1aoVKSkpIb3WuHHjePfdd9m+fTuBQIB+/fpxySWXaC0VkQjlBxKBvwIlgC1AH+Bzl02JSEQLOai88sorTJw4keuuuw5jDH6/nw4dOvDkk0/y8ssvh/RaVatWZc6cOVSvXp309HQ+++wzevTowYcffhhqWyLiWBVgLtAtePwScBtwyFlHIhItQppXiouLM3PnzjW5ubnm2LFj5ujRoyY3N9e8/PLLxu/3e3aOS6VSFV11AfMD9l6Ug2AGeaAnlUrl3Qrl89sX/CJk9evXp1WrVvj9ftavX5+/umxxCgQCZGRkkJCQQGZmZrG/v0is8wMPAQ8Gv/4cO9WzxWVTIuJ5oXx+hzz186Pk5GSSk5NP98dFJMJVB+YBXYLHzwMjsDfPioiES8hBxe/3M3jwYLp27UqVKlXw+ws+ONS1a9ewNSci3tQNez9KFeyTPLcArzjtSESiVchBZerUqQwePJhFixbx+eefY8xpzRyJSAQqATwC3Bc83oCd6vnaVUMiEvVCDir9+vWjT58+vPvuu0XRj4gUM5/fT/3WLUioXImMPXtJTtqIycs77rqawHzsGikAfwfuAo4UX6siEoNCDirZ2dlObpwVkfBr3rUzvcaMokK1qvnn0nal8tbjk9m0ZHn+ucuBl4FKQAYwFHi9mHsVkdgU8sq0Tz31FCNHjiyKXkSkGDXv2plBkyZQvkrlAufLV6nMoEkTaN61M3HAROAdbEhZB7RGIUVEik/IIyodO3akS5cuXH755XzxxRfk5OQU+P61114btuZEpGj4/H56jRkFGHy/uCHe5/dj8vK4YfAALlqynPbB89OAe4HsYu5VRGLbaS2hf6K9fkQkctRv3aLAdM8vNVz2H7o/+BilgTRgCKB/60XEhZCDypAhQ4qiDxEpRgmVK53wvD8nh06Tnqb13AUAbCqfwNXpGXxXjL2JiPxcyPeoiEjky9iz97hz5ben0G/gLfkhZe0N/enToJ5Ciog4FfKISlJS0gnXTjHGcOTIEb755htefPFFli1bFo7+RKQIbN2wibxjx/D5/fh8Phot/ojfjx1PqYOHOJIQ4N3HHuTbTh34X9tLXLcqIjEu5BGV9957j/r163Po0CGWLl3KsmXLOHjwIA0aNGDNmjVUr16dDz/8kKuvvroo+hWRMKjXsjn+EiWIy87m0see4A9330+pg4dIadmcOW+8zNYuF+MvUYJ6LZu7blVEYlzIIyqVKlXiqaee4rHHHitw/v7776dOnTp0796dhx9+mAcffJCFCxeGrVERCZ+EypWosG07V977AFW3/A+AT4cMZNXtN5NXMq7AdSIiLoW8e3JaWhpt2rTh22+/LXC+QYMGrFu3jgoVKtCkSRPWrFlDQkJCOHs9jnZPFjk9I+rX5Yldu4k/fJjDZ1fgvXEP8d3F7Y+77pkb/8y3a9c76FBEolkon98hT/0cOXKEDh06HHe+Q4cOHDliF9P2+/0cPXo01JcWkSJWGpgJTE3+jvjDh9nRugVzX3/puJBi8vI4sHMXyUkbnfQpIvKjkKd+pk+fzsyZM2nTpg1r1qzBGMMFF1zA0KFDGT9+PADdu3dn/Xr9V5iIlzQBFgDnA3nAc3Vrc+i56Zi4Evh+dp3d58fH2xOnnHDPHxGR4hTy1A/AgAEDuP3222nSpAkAX331FdOnT2f+/PkAlC5dGmNMkY+qaOpHpHCux24iWA5IDR5/yIn3+jmwcxdvT5xSYK8fEZFwCuXz+7SCilcoqIicWhlgOnZlWYCPgD8Bu352TWF3TxYRCZdQPr9DnvoRkchwHnaq57fAMeCvwDjstM/Pmbw83TArIp5VqKCyb98+GjduzL59+9i/f/8JF3z7UcWKFcPWnIicnhuBp7EjKj8AAwBN5IhIJCpUUBk1alT+0Mydd95ZlP2IyBkoi70XZWDw+P3g13ucdSQicmZ0j4pIlDgfO9XTBMgFHgQmEsH/gotI1Ar7PSqBQKDQb67AIFL8bgamYtdJ2QH0A1Y67UhEJDwKFVTS0tJOeV8KgM/nwxhDXJzuzxUpLgHgH9hgArAIGATsc9aRiEh4FSpVdOnSpaj7EJEQtcJO9TQEcoBEYBKa6hGR6FKooPLxxx8X6sVatGhxRs2ISOHcBjwFlAK2YUdUVjvtSESkaIS8188vJSQkMHz4cNatW8e6devC0ZOInER54A3so8elgLexIysKKSISrU47qHTp0oU5c+awc+dO7rjjDt555x3atm0bzt5E5Gd+B6wHrgWygZFAL+CAw55ERIpaSHe+1qxZk8GDBzNkyBDKli3LggULKFmyJNdeey1btmwpqh5FYt6d2EeN44FkoC+w1mVDIiLFpNAjKosWLWLz5s2cd9553HHHHdSoUYMRI0YUZW8iMe9s4C1gMjakvAG0RiFFRGJHoUdUfv/73zNt2jT+/ve/88033xRlTyICtAdeBWoDR4FR2FVnRURiSaFHVC6++GICgQBr165l9erV3HbbbVSqVKkoexOJST7gXuBjbEj5GmiHQoqIxC4TSp111lnmxhtvNCtWrDBHjx41ubm5ZsSIEaZcuXIhvU44KhAIGGOMCQQCxf7eKlVRVCUwi8CYYL0CJuCBvlQqlSqcFeLn9+m/UePGjc3EiRPNDz/8YA4fPmzefvttL/+iKpWnqyOYHdiAchjMUA/0pFKpVEVRoXx+n9E6Kv/73/8YPXo05557Lv379z+TlxKJWT7gPmAZUBP4ErgQeN5hTyIiXqHdk0UcqgLMAX4fPH4Z+DNwyFlHIiJFL+y7J4tI+HUB5gHVgcPYZfFfdNmQiIgHnfES+iISGj8wFvgQG1K+ANqikCIiciIaUREpRtWwoyiXBo9nAXcAWc46EhHxNgUVkWLSDZiLvS/lIHArNrSIiMjJaepHpIiVAB4D3sOGlI1AGxRSREQKw2lQGTNmDJ9++ikZGRmkpqby5ptv0rhxY5ctiYRVTeAj4H7sv2wzsavM/s9lUyIiEcRpUOncuTMzZsygXbt2dOvWjbi4OBYvXkyZMmVctiUSFpcDG4BOQAZ2x+PhwBGHPYmIRCLnK9T9WJUqVTLGGHPxxReHfWU7laq4Kg7MRH5aBn8tmAYe6EulUqm8UqF8fnvqZtry5csDsH///hN+Pz4+nlKlSuUfBwKBYulLpLBqYXc87hA8ng7cA2Q760hEJLJ56mbaSZMmsWLFCr744osTfj8xMZGMjIz8SklJKeYORU7uD9ipng5AGnAtMAKFFBGRM+V8CAgwTz/9tNm6daupWbPmSa+Jj483gUAgv2rUqKGpH5XzKglmEj9N9fwXTF0P9KVSqVRerYib+pk2bRpXX301nTp1OuUoSXZ2NtnZ+u9T8Y66wGvABcHjScAYIMdVQyIiUcZ5UJk+fTq9e/fmkksu4bvvvnPdjkihXYNdWbYCsB8YDPzbYT8iItHIaVCZMWMGAwYMoGfPnmRmZlK1alUA0tPTOXJED3GKN5UCngRuDx6vAvoB2511JCIS3ZzNUZ3MoEGDwj7HpVKFoxqAWcdP96M8jn0c2XVfKpVKFUkVMfeo+Hw+l28vEpK+wD+ABGAvcAPwrtOORESin6ceTxbxotLYpe9fxYaUj4GWKKSIiBQHBRWRU2gMrAZuAfKwmwteCmgFHxGR4uH8qR8Rr/oTdiSlHJAKXA986LQjEZHYoxEVkV84C/vY8VxsSPkIO9WjkCIiUvwUVER+5jxgDTAEO9UzFugG7HLZlIhIDNPUj0jQYGAGUAbYCQwAljnsR0RENKIiQlngJWA2NqQsxk71LHPXkoiIBCmoSExrDqzFrolyDLgP6AHsdtmUiIjk09SPxKxhwDTsOik7gP7Af5x2JCIiv6SgIjEnADyLDSYA72BHVPY560hERE5GUz8SU1oC67AhJRe4F7gKhRQREa/SiIrEjD8Dk7C7H2/D7ni82mlHIiLyaxRUJOqVB54H/hg8fhu4ETjgrCMRESksTf1IVGsLJGFDSjZwJ9ALhRQRkUihERWJWiOBvwHxwFagD/ZRZBERiRwKKhJ1zsYu3tYzePxP4CYg3VlHIiJyujT1I1GlHbAeG1KOArdhp30UUkREIpOCikQFH3AP8DFQB/gGaA8847IpERE5Y5r6kYhXEbtXz5XB41eBm4FMZx2JiEi4aERFIlpHYAM2pBzBBpT+KKSIiEQLBRWJSD4gEVgKnAt8CVwAPOeyKRERCTtN/UjEqQzMBX4fPJ4DDAcOOetIRESKioKKRJRLgFeA6sBh7FM9LzrsR0REipamfiQi+IGHgA+xIeUL4HcopIiIRDuNqIjnVcNO9XQNHr8A3IEdURERkeimoCKedhk2pFQFDmLvRZnrtCMRESlOmvoRTyoBPAq8jw0pn2E3GFRIERGJLRpREc+pib1htlPweCYwCrtOioiIxBYFFfGUHtjHjSsBGdgF3F5z2pGIiLikqR/xhDjgceBdbEhJAtqgkCIiEus0oiLO1cLuz9MhePw0doPBo846EhERr1BQEaf+gF0L5RwgDbgJ+JfDfkRExFs09SNOlASeAhZiQ8qnQGsUUkREpCCNqEixq4ud6rkweDwZGA3kuGpIREQ8S0FFzojP76d+6xYkVK5Exp69JCdtxOTlnfT63tiVZSsAB4DB2FEVERGRE1FQkdPWvGtneo0ZRYVqVfPPpe1K5a3HJ7NpyfIC18YDT2KXvgf4BOgHfF9MvYqISGTSPSpyWpp37cygSRMoX6VygfPlq1Rm0KQJNO/aOf9cA2AVP4WUidjF3BRSRETk1yioSMh8fj+9xowCDD6//7jvgaHn6Dvx+f304ac1UfYCVwBjgNxi7llERCKTgoqErH7rFlSoVvW4kPIjn99P5bMrMLfi2bwGJAArgJbYBd1EREQKS0FFQpZQudIpv3/21m30/9MwBuzZRx4wDugCpBRHcyIiElV0M62ELGPP3pN+r+m/36Pro38jPiuLfXFx9M/N5YNi7E1ERKKLgoqELDlpI2m7UilfpXL+9E9c1hEunfAUzd78P3tNi2Z0/OwLdrpsVEREIp6mfiRkJi+Ptx6fDPgweXmc8+1WBgy4iWZv/h/G52PVrUO4puI57DTGdasiIhIFjKu6+OKLzcKFC01KSooxxpiePXuG9POBQMAYY0wgEHD2O8RyNe/a2fzznjtMdulSxoDJrFTRzJr4V9O8a2fnvalUKpXKuxXK57fTqZ+yZcuyceNGZs+ezb/+pV1eIklZ4O4ly7kmuLDb6ornMLJmddYk/vWUK9OKiIiEwmlQee+993jvvfcKfX18fDylSpXKPw4EAkXRlvyK5sBrQFPgGPAQMGHffsy+/U77EhGR6BNR96gkJiaSkZGRXykpeuC1uA0D/osNKSnYx47HY8fnREREwi2igsqECRNISEjIr5o1a7puKWYEgFeAfwBnAe9gF3Bb4bAnERGJfhH1eHJ2djbZ2dmu24g5LYEFQCPs0vf3YTcY1CiKiIgUtYgKKlL8hgOTgVLYTQT7YXc+FhERKQ4KKnJCCcDzwHXB44XAYOCAq4ZERCQmOX88uWHDhvnH9erVo0WLFuzfv5/t27c77Cy2tcU+1VMfyAZGA1NcNiQiIjHN2YIvnTt3Nicye/bssC8YoypcjQBzFIwBkwzmdx7oSaVSqVTRVRGz4Nvy5cvx+XwuW5Cgs4EXgF7B438CNwHprhoSEREhwh5PlqJxIbAeG1KOArcDf0QhRURE3FNQiWE+4B7sWih1gG+ADsAMl02JiIj8jJ76iVEVgReBq4LHrwI3A5muGhIRETkBjajEoIuADdiQcgS4BeiPQoqIiHiPgkoM8QFjgGXAucBX2PtT/uGwJxERkVPR1E+MqAzMAboHj+dgV5095KwjERGRX6egEgM6YzcUrAEcxj7VM9tpRyIiIoWjqZ8o5gceBJZgQ8pm4HcopIiISOTQiEqUqgrMA7oGj18A7sCOqIiIiEQKBZUo1BUbUqpi70G5FZjrtCMREZHTo6mfKFICeARYjA0pnwFtUEgREZHIpRGVKFEDe8Ns5+DxP4CR2HVSREREIpWCShTojn3cuDJ20babsSvNioiIRDpN/USwOGAC8B42pKwHWqOQIiIi0UMjKhGqFjAfuxw+wNPYDQaPOutIREQk/BRUItBV2A0FKwLpwE3AP102JCIiUkQ09RNBSgJPAv/GhpQ1QCsUUkREJHopqESIusAK4O7g8WSgI7DVVUMiIiLFQFM/EaAXdtn7CsABYDCw0F07IiIixUYjKh4WD0wF3sSGlE+AliikiIhI7FBQ8agGwCpgRPD4CaAT8L2zjkRERIqfpn486DrgeSAB2AsMAt5x2pGIiIgbGlHxkNLAM8ACbEhZgZ3qUUgREZFYpaDiEY2B1cDw4PF4oAuQ4qwjERER9zT14wEDgGeBcsBuYCB2B2QREZFYpxEVh84CngPmYUPKUuxUj0KKiIiIpaDiSFPgU2AokAf8FbgM2OmyKREREY/R1I8Dg4AZQFlgF3bqZ6nTjkRERLxJIyrFqCx2M8EXg19/ALRAIUVERORkFFSKSTPsJoKDgGPAA0AP7M2zIiIicmKa+ikGQ4Fp2JtnU7BTPR877UhERCQyKKgUoXLYx44HBI/fBW7ArjYrIiIiv05TP0WkJZCEDSm5wGjgShRSREREQqERlSIwHJiEXRL/e6AfdudjERERCY2CShglYDcTvC54vBC4EdjvrCMREZHIpqmfMGmDneq5DsgB7gJ6opAiIiJyJjSiEgYjgCeAeOA7oC921VkRERE5MwoqZ6AC8ALQO3j8L+AmIM1RPy74/H7qt25BQuVKZOzZS3LSRkxenuu2REQkSiionKYLgVeBusBR4B7gaZcNOdC8a2d6jRlFhWpV88+l7Urlrccns2nJcoediYhItNA9KiHyAXcDK7Ah5VugA7EZUgZNmkD5KpULnC9fpTKDJk2gedfOjjoTEZFooqASgnOwT/I8CZQEXgNaY2+ijSU+v59eY0YBBp/ff9z3wNBz9J3HfU9ERCRU+iQppIuADcBVwBHgVuz6KBkOe3KlfusWVKhW9aRBxOf3c3b1atRv3aKYOxMRkWjjPKgMHz6c5ORksrKyWLt2LR07dnTdUgE+YAywDKgFfIW9P+VZhz25llC5UlivExERORmnQaVPnz5MmTKFcePG0apVK1asWMG7775LrVq1XLaVrzLwDjABe9fxXKAt8JnLpjwgY0/hNgIo7HUiIiIn4zSo3HXXXcyaNYtZs2bx5ZdfMmrUKLZv387w4cNdtgVAJ+xUTw/gMPax44HAQYc9eUVy0kbSdqWe9DFkk5fHgZ27SE7aWMydiYhItHEWVEqWLEmbNm1YvHhxgfOLFy+mQ4cOJ/yZ+Ph4AoFAgQo3P/Ag8BFQA9gMXIBdL0Usk5fHW49PBnzHhRV77OPtiVO0noqIiJwxZ0GlUqVKxMXFkZqaWuB8amoq1apVO+HPJCYmkpGRkV8pKSlh7akq8D7wCFACmA38DvgirO8SHTYtWc5LdyWSvntPgfNpqbt56a5EraMiIiJh4XzBN2NMgWOfz3fcuR9NmDCBSZMm5R8HAoGwhZVLgXlANeAQdgfkOWF55ei1aclyPl+6QivTiohIkXEWVPbu3Utubu5xoydVqlQ5bpTlR9nZ2WRnZ4e1jxLAQ8AD2OGlTUAf4Muwvkv0Mnl5fLt2ves2REQkSjmb+snJyWHdunV069atwPlu3bqxatWqYumhBrAEG1T8wD+w96MopIiIiHiD06mfSZMmMWfOHNauXcsnn3zCzTffTO3atZk5c2aRv3d37NROZSATuAWYX+TvKiIiIqFwGlQWLFhAxYoVeeihh6hevTqff/45V1xxBd9//32RvWcc9mbZxODxBuxUz9dF9o4iIiJyunzAie9cjQCBQICMjAwSEhLIzMz81evPxe54fFHweAZ2g8GjRdeiiIiI/EIon9/On/opLlcCLwEVgXRgKPCG045ERETk1zjf66eolcTudvx/2JCyFrvjsUKKiIiI90X1iEod7FRPu+DxFGA0EN4HnEVERKSoRG1Q6YldWfZs4ABwI/C2045EREQkVFE39ROPHTl5CxtSVgOtUEgRERGJRFEVVOoDK4GRweMnsbsgb3PWkYiIiJyJqJn6+SPwPFAe2AcMAhY57UhERETOVFQElaeAYcGv/wP0B3a4a0dERETCJCqmfn4MKeOBLiikiIiIRIuoGFHZFwgwFLvB4FmumxEREZFTCgQChb42ooPKj79oxZQU3nTci4iIiIQmEAj86hL6Eb3XD0CNGjUKtc9PqAKBACkpKdSsWbNIXj+S6W9zavr7nJr+Pqemv8+p6e9zcpH2twkEAvzwww+/el1Ej6gAhfolz0RmZmZE/AN3QX+bU9Pf59T09zk1/X1OTX+fk4uUv01he4yKm2lFREQkOimoiIiIiGcpqJzE0aNHefjhhzl69KjrVjxHf5tT09/n1PT3OTX9fU5Nf5+Ti9a/TcTfTCsiIiLRSyMqIiIi4lkKKiIiIuJZCioiIiLiWQoqIiIi4lkKKicwfPhwkpOTycrKYu3atXTs2NF1S55x8cUXs3DhQlJSUjDG0LNnT9ctecaYMWP49NNPycjIIDU1lTfffJPGjRu7bsszbr31VjZu3Eh6ejrp6emsWrWKHj16uG7Lk8aMGYMxhsmTJ7tuxRPGjh2LMaZA7dy503VbnlKjRg3mzJnD3r17OXToEOvXr6d169au2woLBZVf6NOnD1OmTGHcuHG0atWKFStW8O6771KrVi3XrXlC2bJl2bhxI7fffrvrVjync+fOzJgxg3bt2tGtWzfi4uJYvHgxZcqUcd2aJ+zYsYMxY8bQtm1b2rZty0cffcTbb7/Neeed57o1T2nbti0333wzGzdudN2Kp3z++edUq1Ytv5o3b+66Jc+oUKECK1euJCcnh8svv5zzzjuPu+++m7S0NNethY1R/VSrV682zzzzTIFzmzdvNuPHj3fem9fKGGN69uzpvA+vVqVKlYwxxlx88cXOe/Fq7du3zwwZMsR5H16psmXLmq+++sp07drVLF261EyePNl5T16osWPHmvXr1zvvw6s1YcIE8/HHHzvvo6hKIyo/U7JkSdq0acPixYsLnF+8eDEdOnRw1JVEqvLlywOwf/9+x514j9/vp2/fvpQtW5ZPPvnEdTueMWPGDBYtWsSSJUtct+I5jRo1IiUlheTkZObPn0+9evVct+QZV199NWvXrmXBggWkpqaSlJTE0KFDXbcVNgoqP1OpUiXi4uJITU0tcD41NZVq1ao56koi1aRJk1ixYgVffPGF61Y8o1mzZmRmZnL06FFmzpxJ79692bJli+u2PKFv3760bt2axMRE1614zn//+19uuOEGunfvzrBhw6hWrRqrVq3inHPOcd2aJ9SvX5/hw4fz9ddf0717d2bOnMm0adMYOHCg69bCxvmwjleqevXqxhhj2rVrV+D8fffdZ7Zs2eK8P6+Vpn5OXk8//bTZunWrqVmzpvNevFQlS5Y0DRo0MG3atDHjx483u3fvNk2bNnXel+s699xzza5du8z555+ff05TPyevMmXKmJ07d5pRo0Y578ULdfToUbNy5coC56ZOnWpWrVrlvLdwlEZUfmbv3r3k5uYeN3pSpUqV40ZZRE5m2rRpXH311XTp0oWUlBTX7XhKTk4O3377LevWreO+++5j48aNjBw50nVbzrVp04aqVauybt06cnJyyMnJ4ZJLLmHEiBHk5OTg9+v/qn/u8OHDbNq0iUaNGrluxRN27tzJ5s2bC5zbsmULtWvXdtRReOl//T+Tk5PDunXr6NatW4Hz3bp1Y9WqVY66kkgyffp0rrnmGi699FK+++471+14ns/no1SpUq7bcG7JkiU0a9aMli1b5teaNWuYN28eLVu2JC8vz3WLnhIfH0/Tpk31iHLQypUradKkSYFzjRs3Ztu2bY46Cj/nwzpeqj59+pijR4+aG2+80fzmN78xkyZNMpmZmaZ27drOe/NClS1b1rRo0cK0aNHCGGPMnXfeaVq0aGFq1arlvDfXNWPGDHPgwAHTqVMnU7Vq1fwqXbq08968UOPGjTMdO3Y0derUMc2aNTOPPfaYyc3NNZdddpnz3rxYmvr5qZ544gnTqVMnU7duXXPBBReYhQsXmvT0dP3/crDatm1rsrOzTWJiomnQoIHp37+/OXjwoBkwYIDz3sJUzhvwXA0fPtxs3brVHDlyxKxdu1aPl/6sOnfubE5k9uzZzntzXSczaNAg5715oZ5//vn8f69SU1PNBx98oJByilJQ+anmz59vUlJSzNGjR82OHTvMG2+8oXubflFXXnml+eyzz0xWVpbZvHmzGTp0qPOewlW+4BciIiIinqN7VERERMSzFFRERETEsxRURERExLMUVERERMSzFFRERETEsxRURERExLMUVERERMSzFFRERETEsxRURCTibN26VZsZisQIBRURwRhzypo9e7brFkUkRsW5bkBE3KtWrVr+13379uWRRx4psBtrVlZWgevj4uLIzc0ttv5EJHZpREVESE1Nza/09HSMMfnHpUuXJj09neuuu46lS5eSlZXF9ddfz9ixY1m/fn2B1xk5ciRbt24tcG7w4MFs3ryZrKwstmzZwvDhw0/Zy9KlS5k+fTrTp0/nwIED7N27l0cfffS468qUKcOsWbPIyMhg27ZtDBs2rMD3H3/8cb766isOHTrEt99+yyOPPEJc3E//bXb++efz0UcfkZGRQXp6OmvXrqVNmzb532/fvj3Lly/n8OHDfP/990ydOpUyZcoU+m8qIuHjfGdElUrlnRo0aJA5cOBA/nGdOnWMMcYkJyeb3r17m7p165rq1aubsWPHmvXr1xf42ZEjR5qtW7fmHw8dOtSkpKTk/1zv3r3N3r17zQ033HDS91+6dKnJyMgwkydPNo0bNzYDBgwwBw8eLLAb7NatW83evXvN8OHDTYMGDczo0aNNbm6uadKkSf41999/v2nfvr2pU6eOueqqq8zOnTvNvffem//9TZs2mZdfftk0adLENGzY0Pzxj380559/vgFMs2bNTEZGhhk5cqRp2LChad++vVm3bp154YUXnP/zUalisJw3oFKpPFQnCyojRowocF1hgsq2bdtMv379Clxz//33m5UrV570/ZcuXWq++OKLAucmTJhQ4NzWrVvNyy+/XOCaXbt2mVtuueWkr3vPPfeYNWvW5B+np6efNDC99NJLZubMmQXOXXTRRSY3N9eUKlXK+T8jlSqWSveoiEihrF27NqTrK1WqRO3atZk1axbPPfdc/vm4uDjS09NP+bOrV68ucPzJJ59w99134/f7ycvLA+Czzz4rcM2uXbuoUqVK/vG1117LnXfeScOGDSlXrhxxcXFkZGTkf3/SpEk8//zzDBw4kA8//JDXX3+d5ORkANq0aUPDhg3505/+lH+9z+ejRIkS1KtXjy+//DKkv4WInD7doyIihXLo0KECx3l5efh8vgLnSpYsmf+132//72XYsGG0bNkyv5o1a0a7du3OuJ+cnJwCx8aY/Pe88MILefXVV3n33Xe56qqraNWqFePGjSM+Pj7/+r/+9a/89re/ZdGiRVx66aVs3ryZXr165ff+7LPPFui7RYsWNGzYkG+//faMexeRwtOIioiclj179hR4WgigZcuW+V/v3r2bHTt2UL9+fV555ZWQXvuXQaZdu3Z8/fXX+aMpv+aiiy5i27ZtjB8/Pv9cnTp1jrvu66+/ZsqUKUyZMoVXXnmFG2+8kbfeeoukpCR++9vfKpSIeIBGVETktCxbtozKlSvzl7/8hfr16/PnP/+Zyy+/vMA1Dz/8MImJiYwYMYJGjRrRrFkzBg8ezKhRo0752rVq1eKpp56icePG9OvXjzvuuIOpU6cWurdvvvmG2rVr07dvX+rXr88dd9xB7969879funRppk+fTufOnalduzYdOnTgd7/7HVu2bAFg4sSJtG/fnqeffjp/JOUPf/gD06ZNC+EvJCLh4vxGGZVK5Z062c20LVq0OO7aW265xWzbts1kZmaaF1980SQmJha4mRYw/fv3N0lJSebIkSNm3759ZtmyZaZXr14nff+lS5eap59+2jzzzDMmLS3N7Nu3z4wfP77ANVu3bjUjR44scG79+vVm7Nix+ccTJ040e/bsMRkZGWb+/Plm5MiR+b9XyZIlzSuvvGK2bdtmjhw5Ynbs2GGmTZtW4EbZtm3bmvfff99kZGSYzMxMs2HDBpOYmOj8n49KFWvlC34hIuIJS5cuZcOGDb866iIisUFTPyIiIuJZCioiIiLiWZr6EREREc/SiIqIiIh4loKKiIiIeJaCioiIiHiWgoqIiIh4loKKiIiIeJaCioiIiHiWgoqIiIh4loKKiIiIeNb/A1RT+NeJErJwAAAAAElFTkSuQmCC",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.scatter(true_phase, phase_aligned)\n",
    "plt.xlabel(\"True phase\")\n",
    "plt.ylabel(\"Aligned phase\")\n",
    "plt.axline([0, 0], [2*np.pi, 2*np.pi], color=\"red\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "numpyro",
   "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.12.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}