{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "cbde3d9d-abc3-4e48-83c7-ff04da4cb99f",
   "metadata": {},
   "source": [
    "# PLSC (behavior/seed PLSC)\n",
    "\n",
    "Partial least squares correlation (PLSC; [Abdi & Williams, 2012](https://doi.org/10.1007/978-1-62703-059-5_23)), also known as behaviour PLSC or seed PLSC ([Krishnan et al., 2011](https://doi.org/10.1016/j.neuroimage.2010.07.034)), is a method for identifying latent patterns of shared variance between two sets of variables by decomposing correlation matrices that are stratified by experimental condition."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f54f0f3d-477b-4447-b797-fb00ab00cc44",
   "metadata": {},
   "source": [
    "## Setting up simulated data\n",
    "\n",
    "We'll simulate data in which a single latent variable explains the relatinoship between a set of 50 observed variables in the data array and a set of 4 covariates. The latent variable will have a sinusoidal pattern of loadings onto the observed variables, and positive loadings onto all of the covariates. First, we will import the necessary libraries and seed the random number generator for reproducibility:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "c70a33cd-981f-4459-aa41-4731cbbab0d3",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import pandas as pd\n",
    "from pyplsc import PLSC\n",
    "from matplotlib import pyplot as plt\n",
    "\n",
    "np.random.seed(123)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cdf5abba-f8cf-40c0-b78a-c53396a29e0f",
   "metadata": {},
   "source": [
    "Next, we will simulate the underlying latent variable, which will be normally distributed."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "d3a043c0-819f-40d9-bd28-a2e4599bc9b4",
   "metadata": {},
   "outputs": [],
   "source": [
    "n_subj = 40\n",
    "latent_var = np.random.normal(size=(n_subj, 1))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ce8ec4a8-d875-4294-abf5-19a78cac9bf4",
   "metadata": {},
   "source": [
    "Next, we will create a set of sinusoidal loadings:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "9a90f0fe-88c3-44b0-82ea-a7c97de257aa",
   "metadata": {},
   "outputs": [],
   "source": [
    "n_var = 50\n",
    "data_loadings = np.sin(0.2*np.arange(n_var)).reshape(1, n_var) # Sinusoidal loading pattern"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1457b8bc-4f23-4639-9d16-26a894fdf8dd",
   "metadata": {},
   "source": [
    "The following displays the pattern of loadings onto the observed variables:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "3291e3e6-260b-4b44-9946-a300c1d87548",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Text(0, 0.5, 'Correlation weight')"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAksAAAGwCAYAAAC5ACFFAAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjgsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvwVt1zgAAAAlwSFlzAAAPYQAAD2EBqD+naQAAcBdJREFUeJzt3Qd4VNXWBuAvvZHeIaH33iO9SlVBsV68gCheC/bKf0XsXBtW7CJyRVG8oNiQ3nuXTkILkB7SSZ//WXvmjAkkIWUm0773eQ6ZlpkzJ8PMmr3XXstJp9PpQEREREQVcq74YiIiIiJisERERER0FRxZIiIiIqoCgyUiIiKiKjBYIiIiIqoCgyUiIiKiKjBYIiIiIqqCa1VXUvWUlpbiwoUL8PX1hZOTEw8bERGRDZBSk9nZ2WjYsCGcnSsfP2KwZAISKEVHR5viroiIiKiexcfHIyoqqtLrGSyZgIwoaQfbz8/PFHdJREREZpaVlaUGO7TP8cowWDIBbepNAiUGS0RERLblaik0TPAmIiIiqgKDJSIiIqIqMFgiIiIiqgKDJSIiIqIqMFgiIiIiqgKDJSIiIqIqMFgiIiIiqgKDJSIiIqIqMFgiIiIiqgKDJSIiIiJ7CZY2bNiA66+/XnUHltLkP/3001V/Z926dejevTs8PDzQsmVLzJ8//4rbzJ07F02bNoWnpydiYmKwY8cOMz0DIiIisjU2FSzl5uaiS5cuKripjlOnTmHs2LEYMmQI9u3bh0cffRT33HMP/vzzT+Ntvv/+ezz++OOYNWsW9uzZo+5/5MiRSE5ONuMzISIiIlvhpNPpdLBBMrK0dOlSjB8/vtLbPPPMM/jtt99w8OBB42W33347MjIysHz5cnVeRpJ69eqFDz/8UJ0vLS1VHYgfeughPPvss9XuWuzv74/MzEw20rWA3IJieLq5wMW56kaIREREtfn8doUd27p1K4YPH17uMhk1khEmUVhYiN27d2PGjBnG652dndXvyO9WpqCgQG1lDzbVn+z8Imw/mY4tcWnYEpeKo4nZiPT3xD0DmuOO3tHwdrfrlzURWdCesxfx1p/HcDo1Fy3CGqBVmC9ahzdAq3BftApvAD9PN/597JBdf6okJiYiPDy83GVyXoKbS5cu4eLFiygpKanwNkePHq30fmfPno0XX3zRbPtN5V0qLMHuMxdVYCQB0l/nM1FSWn5ANCEzHy//ehgfrDmBKX2bYnKfpgj0ceehJCKTiE/Pwxt/HsMv+y8YL7uQmY+NJ1LL3U6+uEng1DqsAbo2DsDojpEc9bYDdh0smYuMREmek0aCL5m6I9OKTc7BC8sOYcepdBSWlJa7rlmID/q0CEbfFsHo0SQQ646l4NP1cTidlod3V53AZxtO4o7ejXHPgGaI9Pfin4aIaiUrvwgfrY3DvM2nUFhcCicn4ObuUbipexTOpOXieFIOTiRn43hSNpKyCtQXN9k2HE9Rvz+mUwLm3NpVpQqQ7bLrYCkiIgJJSUnlLpPzMi/p5eUFFxcXtVV0G/ndysjKOtnIfI4mZmHi59uRllto/LYmwVG/FiHqZ8OA8gGQBEa39ozGHwcT8PG6OBy6kIUvN53Cgq2ncWO3RvjXoBZoEdqAfzIiqpbiklJ8t+Ms3ll1AumG96F+LYPxf2PaoUNDf3Ve3ovKyswrMgROOeo9bNGOePz+VyJSs3fg80k94e/NKTpbZdfBUp8+ffD777+Xu2zlypXqcuHu7o4ePXpg9erVxkRxSfCW89OnT7fIPhNw6EIm7vxiOy7mFaFDQz+8d3s3tAj1UUn9VZEE7+s6N8TYTpHYcCIVH6+LxbaT6fhh1zks3n0OE7pHYfZNneDmYlOLQImoHsmaJxmpfvX3I2p0W8j7z7/HtsOQNmFVvg9JMNSzaZDaxKiOEfjXgt3YcTodN3+yBV9P7X3FFz2yDTYVLOXk5CA2NrZcaQApCRAUFITGjRur6bHz589jwYIF6vr77rtPrXJ7+umnMXXqVKxZswY//PCDWiGnkem0yZMno2fPnujduzfeffddVaLgrrvusshzdHQHzmXgn1/uQOalInSJ8seCqTE1/jYmb2aDWoeqTZIxZaRp5eEk/Lj7HDxcnfHK+I5XDbyIyDHzIx9YuBtrj+mn0IJ83PHY8Fa4vXfjWn3J6tsiBD/c1wdTvtqBE8k5uPGjzZh/V2+0i6x81RVZJ5sqHSAFJqVm0uUk2JFik1OmTMHp06fV7cr+zmOPPYbDhw8jKioKM2fOVLcrSwKqN998UyWEd+3aFe+//74qKVBdLB1gGhLYTP5yB7ILitG9cQDmT+1tspUlKw4l4l/f7Ia82mde1x53929mkvslIvsgH4WPfr8PP++7AHcXZ9zVvykeHNLSJO9B5zMuYco8fcDk6+GKTyf1UIEUWV51P79tKliyVgyW6m7n6XTc9dVO5BQUo3fTIMy7qxcaeJh24POLjSfxym9HVILmF5N6Yli78qsgichxzdt0Ci/9ehiuzk5YeE8MYpqXz0eqK8lnmrZgl5qSc3Nxwtu3dsUNXRqa9DHIfJ/fTN4gi9sal4bJ83aoQKlP82DMn2r6QEnIaJIkgsvXg4e+24vDF1gfi4iAbSfTVI6SkNwkUwdKQtIJFtzdG2M6RaCoRIeHv9urvsCRbWCwRBa16UQq7pq/A3mFJRjQKgTzpvQyW1FJyVN6aVwHtaJFHu+er3ciOSvfLI9FRLYhMTMf07/do2q3je/aUNVpMxcpH/DBHd2NjyEj3S/9chill9WNI+vDYIksZt2xZEz9eifyi0oxpE2oWlrr5W7eWiSSpPnRP3qgeaiPKignw+KS1ElEjqeguAT3L9yN1JxClXQ9+6bOZl/8Iat2Z13fHjNGt1XnpX7Tu6uOm/Uxqe4YLJFFbDyRgnsX7FZF3q5tH45P/tmj3oq2yXD4V1N6IdDbDfvPZeKJxfv4zY7IAcmozt6zGfDzdMWnd/Yw+5c1jQRkUvvt9Qmd1PmP1sXhWGJ2vTw21Q6DJbJI49unFh9QVblHd4zARxO7w8O1fqvbNgn2waf/7KkSLaVo3Nsrj9Xr4xORZf2wMx4Lt59VCz7eu6MbGgd71/s+SCFd+bJYXKrDs0sO8EubFWOwRPXu/dUnkJiVj+ggL7xzW1eLFYns3SwI/7mpszo9d20c/rf7nEX2g4jq1/74DDz380F1+vHhrVWxSUvQ8ihlQYuMcC3cfsYi+0FXx2CJ6tWJpGzVhkS8cH0Hi/dLmtAjCg8OaaFOyzc76UNHRPYrLacA93+jTwEY3i5c1VKyJOld+fSoNur068uPISHzkkX3hyrGYInqjZT0mvnzQTXkLG9S1lLn6Ilr2xiX8/7rv7vU6hgiss9+b1I2RBZ3NA/xwZzbusDZ2fLV/CfGNEG3xgGqfMqsnw9ZeneoAgyWqN4s239B9WqTliOyGsRayJvl27d0RcdGfqof3TsruTKFyB69+ecxbIlLg7e7Cz79Zw+TdQgwxQo5SQmQgpgrDidh+cEES+8SXYbBEtWL7PwivPqbvujb9CEtER1U/8mUVZFVMC/e0FGdXrw7Xk0XEpH92HA8BZ9u0BeBfOuWLmgV7gtr0ibCF/cN0qcEPP/zIWTlF1l6l6gMBktUL95ddQLJ2QVoGuyNaQObW+VR79EkECM7hEPqw72+/Kild4eITESKPmr/pyf3aYIxnSKt8thOH9oSzUJ81Hvl63/wPciaMFgiszuamIX5W06r0y/cYPmk7qo8PaqtGhJfdSSZyd5EdmL5oUQcupClVp09Mrw1rJW8N752o772kpQ1kJ6ZZB0YLJH5k7p/OqhaCYzqEIHBFlqiW10tQhvgtl7R6vTsP46o/Sci2yXvPXMMeYhT+zdDkI87rFmfFsG4raf+PWjGkr9UlXGyPAZLZFZL9pzHztMX4eXmgplWlNRdlUeHtVL7K3VPlh9MtPTuEFEd/LT3PGKTc+Dv5YZ7BjSziWM5Y0xbhDRwV/v9yTo227UGDJbIbDIvFanRGfHQsJZoFOBlE0c7zM8T0wxvqrJ6pqik1NK7RES1ILWU3l2tH1W6f3ALq1n9djUB3u6YdX0HdXru2lgVNJFlMVgis5mz4phqUClNa+/pb51J3ZWRJPRgH3ecTM3F9zvjLb07RFQL3++KR3z6JYT6emByn6Y2dQyv6xypGoxLW6j/W/IXW6FYGIMlMouD5zPx32360v0vj+sId1fbeqn5errhoaEtjSv5pJ8dEdmO/KISfLjmhLFcSX01yTVlK5SXx3dUNaF2nE5XgR9Zjm19gpHNLNOVSt2yBF++HfVrGQJb9I+YJmgS7I3UnAJ8sVHfooWIbMN/t55BUlaBmv6/vbc+YdrWRAV64/Fr9av33lt1Qk0rkmUwWCKT+3H3OZUc7ePugufG2kZSd0VkNOzJEfqeTZ9tiFNBExHZRhHcj9bFqtOPDGsFD1fbGlUq6599miDM10M1H/9l/wVL747DYrBEJiXffN5eeUydfnR4a0T4e9r0ER7bKRKdo/yRW1iC91frh/SJyLrN23RatS6S/m83dW8EWyaB3pR++nyrzzeeZDkTC2GwRCb1x8EENfQtCZWT+jax+aMrfeOeHd1Wnf52+1mcSs219C4RURUy8grxxUb9cvvHrm0NVxfb/5ib2LuJyl06mpiNDSdSLb07Dsn2X0VkNaSA45eb9Lk9k65pYtND32X1bRGCQa1DUVyqw1sr9KNmRGSdPll/EtkFxWgX6adGhu2Bv7cbbu/VWJ3+3NDfjuoXgyUymV1nLuLAuUx4uDpj4jW2P6pUlowuOTkBvx1IwP74DEvvDhFVIDk7H/O36L+wPXFtazUybC+m9m+qWjFtik1Vq42pfjFYIpP50rBiTHIErL2lQE3Jt9Qbu+lzH9gGhcg6zV0Ti/yiUnRrHIBh7ay7tVJtVsZpI2XaNCPVHwZLZBLx6XlYcVjfGuSufrbRUqCmnhjRRq2Q23YyHeuPp1h6d4iojHMX8/DtjrPq9FMj2qg6RfZm2gB9cd9fDiTgfMYlS++OQ2GwRCYxf8tpVVdpQKsQtA73tcujKvVa7ozRTy/O23za0rtDRGXIatWiEh36tghGXxut7XY1naL80ad5sGoO/JUhP5TqB4MlMklNE60lyN397XNUSTPZsMJvw/EUnEnjyjgia3AyJQf/23NenX5ypL42mr26d5B+dOm7HWdV/02qHwyWqM5+2HUOOQXFaBHqg4GtQu36iDYJ9sHA1qHGUgJEZHmfrj+pRluGtQ1D98aBsGeDW4eidXgDVftNAiaqHwyWqE7kDUpbfTK1fzO7Wn1SmTtj9Et4f9gVr/pPEZFlR7aXGSpb/2tQC7v/U0gulpa79NXmU2yBUk8YLFGdrDycpLp6B3i74aZuUQ5xNIe2DUNDf09VIViKcBKR5fy07wIuFZWgZVgD9Gpq36NKmhu6NlQtUKQAsBYoknkxWKI6mWdIMpwY09jmunrXllQEvqN3Y2OzTiKyXCFcbTr8H70b2+UKuIpIwV9t1bEUqZTjQOZlc8HS3Llz0bRpU3h6eiImJgY7duyo9LaDBw9W/3ku38aOHWu8zZQpU664ftSoUfX0bGzbX+cyseN0OlydnfDPa/S9ixzFbb2j1fPeczYDhy9kWXp3iBzS/nOZOJKQpUp62HoPuJr6R0xj1az8WBJboNQHmwqWvv/+ezz++OOYNWsW9uzZgy5dumDkyJFITk6u8PZLlixBQkKCcTt48CBcXFxwyy23lLudBEdlb/fdd9/V0zOybfM260eVruscafMNc2sqzNcTIztEqNPfbOfoEpElfGcYVbquUyQCvO2rEO7V+Hu54TZDC5TPNsRZenfsnk0FS3PmzMG0adNw1113oX379vjkk0/g7e2NefPmVXj7oKAgREREGLeVK1eq218eLHl4eJS7XWCgY8x710VSVj5+McyV391fn2zoaCZeo3+j+mnveZVkSkSWSey+w7DowtFoLVA2x6axBYqZ2UywVFhYiN27d2P48OHGy5ydndX5rVu3Vus+vvzyS9x+++3w8fEpd/m6desQFhaGNm3a4P7770daWlqV91NQUICsrKxym6NZsPW0aizbu2mQKpTmiKQ4nJRLyCssUQETEVkmsbtnE8f8glu2BcrnbIFiVjYTLKWmpqKkpATh4eHlLpfziYn6NhtVkdwmmYa75557rpiCW7BgAVavXo3XX38d69evx+jRo9VjVWb27Nnw9/c3btHR0XAklwpLsNAw/C3fbByV5LfdaWgY/M22s0yyJKonjprYXZF7B+pH9n9lCxSzsplgqa5kVKlTp07o3bt3uctlpOmGG25Q140fPx6//vordu7cqUabKjNjxgxkZmYat/h4ffVqR7F073lk5BUhOsgL17bX5+04qpu6R8HLTZ9kufP0RUvvDpFDcOTE7st1bOSvWrywBYp52UywFBISopKzk5KSyl0u5yXPqCq5ublYtGgR7r777qs+TvPmzdVjxcbGVnobyXHy8/MrtznSNzotsXtK32ZqvtyRSZLlDV0aqtPfbGOiN1F9JnaPdcDE7qpGlxbtjFcj/+TAwZK7uzt69Oihpss0paWl6nyfPn2q/N3FixerPKM777zzqo9z7tw5lbMUGamfB6by1h9PQWxyDhp4uOLWno5RhPJqtKk4KVCZmlNg6d0hcpjEblk+T1BtpqICvVTbqZVHyg8okIMFS0LKBnz++ef4+uuvceTIEZWMLaNGsjpOTJo0SU2RVTQFJ1NswcHB5S7PycnBU089hW3btuH06dMq8Bo3bhxatmypShLQlb40FKG8tWc0fD3deIgMncC7RPmrjufSAoWIzIeJ3VeSNlM3dtNPRy7Zc44vP0cPlm677Ta89dZbeP7559G1a1fs27cPy5cvNyZ9nz17VtVJKuvYsWPYtGlThVNwMq134MABlbPUunVrdRsZvdq4caOaaqPy4tPzsPFEKiSX8q5+jpvYXdXokiSdSu4AEZk3sfsOB0/svpwWLMl7dHJ2vqV3x+64wsZMnz5dbRWpKClbygFUVgrey8sLf/75p8n30V79vE+/PL5fixBEB3lbenesyvVdGuKV347g3MVLWH88GUPbll+1SUSmTeye4OCJ3ZdrHtoA3RoHYO/ZDCzbdwH3GJrtkgOOLJHlSMApq+DEuK76hGb6m6ebC27uEWUsI0BEpsfE7qrdZBhd0t6ryXQYLFG1HDyfhbiUXHi4OmNUR8cuF1AZaSYs1h5LVlOWRGQ6TOy+uus6N4SbixMOXcjC0UTHK5ZsTgyWqFp+MkzBXds+nIndVQyD92sZDJn1/W4HR5eITImJ3VcX6OOOIW3C1Omlezi6ZEoMluiqiktKjUt1x3dlnkBV7ozRJ3rLqriCYtY7ITIFJnbXrFCu9gWXi01Mh8ESXdWWuDSkZBcg0NsNA1uH8ohVYXj7cIT7eSA1pxDLD169DQ8RXR0Tu6tvSNtQVSw3KasAW+JS+fIyEQZLVO0puLGdI9UqFKqcm4szbu+lz11avIv1TohMgYnd1efh6oLru+iLKnMqznT4yUdVyissxp+GERKtjgdVTetVJd/qWNGbyHSJ3VJbiao/FffHwUTkFhTzkJkAgyWq0srDScgtLFFNc7s3DuTRqoYmwT7oHOUPqU35x1/li6QSUe0Tu3s15XtQdXSLDkCzEB913JgOYBoMlqhKP+/7O7Gb1XKr7/rO+lpUvxxgsERUF4sNLYRu7xXN96BqkvdqbSaANZdMg8ESVSotp0A1zhXjuAquRiS/S+w8nY7ETLYeIKqNs2l5OHAuE85OwHimAdSIFixtjktFQuYlvgDriMESVeq3vxLU0lOZUpIhcKq+hgFe6NkkUNVckuNIRDX3+0H9/50+LYIR0oD9OmtCWlL1bhqk3oO0GQKqPQZLVKm/25swsbs2rjOMLv1iSE4lopr5zTCNPbYTWyzVZbHJkj3nKu2RStXDYIkqdDo1VzVklOFvbRkq1cyYzpHq+O2Lz2D7E6IaOpOWi7/OZ8LF2QkjO7AxdW2M7qQv93I8KUe1QKHaY7BEFdKGbfu1DEGYryePUi3IcYtpFqxO/8pEb6Ia0aav+zQPRjCn4GpFilNKiyqxhO1P6oTBEl1Bhmu1QpSsrVQ313fRTx/8eoBTcUS1moIzTGdT7dxkSPRetv+8al1FtcNgiSpsLXAqNRdebi4Y2SGCR6gORnWMUNMIMgR+MiWHx5KoGuT9R/7P6Kfg+B5UF9KiKtjHXbVg2niC7U9qi8ESXeEnQ2K3DN/6eLjyCNVBkI87+rcMUac5FUdUPb8bpuD6tghW/4eobi2YbuiqH+H+3x62YKotBktUTlFJqXH1FqfgTIOr4ohqRvtiof3fobq5qVuUsSNDVn4RD2ctMFiicjbFpiItt1AN2/ZvpR8RoboZ0SEC7i7OOJGcg2OJ2TycRFWIS8nBkYQsuDo7YUR7TsGZQsdGfmgV1gAFxaVswVRLDJaonJ8NU3DyjU6Gb8k0K1Ikb0Aw0Zuoar8bRpVkJW4gp+BM1/7EUHPpf1wVVyv8NCQj6U7956EkdZqtBUxLq1UlU5wsDkd09ZIBXAVnWvr+nsCOU+k4n8H2JzXFYImMZD5bulQ3DfZG1+gAHhkTGt4uHJ5uzjidlsficESViE3OwdHEbLi5OGEkp+BM3oKpV5MgdXrFoUS+BmuIwRJV2N5Ehm3JdGRV4dC2Yeo0258QVb0KTlaQ+nu78TCZ2AhDJfQ/GSzVGIMlUlKyC1Ryt+AUnHlc31krUJnAqTiiKgtRshecOWg1q2QqLj23kK/BGmCwRMYpuJJSHTpH+aNZiA+PihkMaRsGH3cXlS+w52wGjzFRGSeSsnEsST8Fp7XoINOKDvJG+0g/lOqAVUf0+alUPQyWSFlxWD+HzWq55uPp5mL8EOCqOKKKE7sHtApVK0jJPLT3eOYt1QyDJUJ2fhG2xKapI8Hu3uZ1nWF6QaYbZCSPiC6bguvEQpTmNLKj/gvbhhOpagU0VQ+DJcL64ykoLClF8xAftAhtwCNiRgNah8DP0xXJ2QXYeTqdx5oIwPGkbFW0VYq3DucUnFm1CfdFk2BvFBaXqvd+qh4GS4QVhtpK13YI5yo4M/Nw/bs5MVfFEZVvbzKwdQin4MxMVjpr70FcFVd9DJYcnHy7WHssWZ1ma4H6cX0X/VTcHwcTUVxSWk+PSmSdpEirVjKAhSjrh5ZuseZosvoMoKtjsOTgtp9KQ3Z+MUIaeKAbC1HWC62Tuizd3RKnzxUjclTHk3JUMUp3V2dVvJXMr1t0IEJ9PdR7/9aTfA+yy2Bp7ty5aNq0KTw9PRETE4MdO3ZUetv58+erIceym/ze5d9qnn/+eURGRsLLywvDhw/HiRMn4HBTcO3D4OzMQpT1wdXFGaM76ofBuSqOHN1vBy6on4Nah8LXk6vg6oO812srczkVZ4fB0vfff4/HH38cs2bNwp49e9ClSxeMHDkSycn6aaSK+Pn5ISEhwbidOXOm3PVvvPEG3n//fXzyySfYvn07fHx81H3m5+fD3pWW6lR9JcEpOMusipNefJyKI0clX1Z/NUzBSfNuqj9a3pJ8BshnAdlRsDRnzhxMmzYNd911F9q3b68CHG9vb8ybN6/S35HRpIiICOMWHh5e7j/qu+++i+eeew7jxo1D586dsWDBAly4cAE//fQT7N1f5zORmJWvCiX2aRFs6d1xKL2aBiLA2w2Zl4qw+8xFS+8OkUVIH7iTKblqCm4Yp+DqVZ/mwfD1cFXdG/bG8z3IboKlwsJC7N69W02TaZydndX5rVu3Vvp7OTk5aNKkCaKjo1VAdOjQIeN1p06dQmJiYrn79Pf3V9N7Vd1nQUEBsrKyym22XIhycJswVTCR6ncqbmgbfa84VtIlR6+tNKRNKBp4uFp6dxyKBKjSVUAb4SY7CZZSU1NRUlJSbmRIyHkJeCrSpk0bNer0888/45tvvkFpaSn69u2Lc+fOqeu136vJfYrZs2eroErbJBCzRcYpOMPKCKpf2jfp1Ucqn0Ymslcysq9V7WYvOMsoW0JA/h5kB8FSbfTp0weTJk1C165dMWjQICxZsgShoaH49NNP63S/M2bMQGZmpnGLj4+HrTmVmqtWobg6O6mRJap/UlNG+mCdTM1FXEoO/wTkUOQ1L+9DUohyqGGEg+rX4DahaoTpTFqe6stHdhAshYSEwMXFBUlJ5YcL5bzkIlWHm5sbunXrhtjYWHVe+72a3qeHh4dKHC+72ZqVhim4a5oHswichcjKHzn+YjWbWpKD0UZUr2kRzCk4C/HxcMWAliHq9J8HORVnF8GSu7s7evTogdWrVxsvk2k1OS8jSNUh03h//fWXKhMgmjVrpoKisvcp+UeyKq6692nrJQM4BWdZWl2ZVYc5FUeOGSwNb8dRJUtiNW87C5aElA34/PPP8fXXX+PIkSO4//77kZubq1bHCZlykykyzUsvvYQVK1bg5MmTqtTAnXfeqUoH3HPPPcaVco8++iheeeUVLFu2TAVSch8NGzbE+PHjYa9k9cPus/rVDywCZ1nDDB8Uu86k42JuoYX3hqh+yGtdXvOCU3CWfw+SEnuHE7IQn55n4b2xXja1/OC2225DSkqKKiIpCdiSi7R8+XJjgvbZs2fVCjnNxYsXVakBuW1gYKAamdqyZYsqO6B5+umnVcB17733IiMjA/3791f3eXnxSnsiUz6Sy9c5yh8NA7wsvTsOLSrQG20jfNUSamk7c1P3KEvvEpHZSQNXKe0jr335P0CWE9zAA72aBmH7qXSV6H3PgOb8c1TASccU+DqTqTtZFSfJ3raQvzR1/k7VE+jJEa0xfWgrS++Ow3t7xTF8sCYWYztFYu7E7g5/PMj+Tf92j2qe++CQFnhqZFtL747Dm7fpFF769TB6Nw3CD/fZdwpKbT+/bWoajuout6AYm2JT1ekRhmWjZB0lBOTbNptakr0rKilVr3XBQpTWQctd3XkmHak5BZbeHavEYMnBbDB8IDcJ9karsAaW3h0C0LmRv2pqmVNQrBobE9mznafTVQPXYB93dIkKsPTukCEdoGMjP5WescpQf4/KY7DkYFYYe8GFqwR3so6mlsMMdWZYoJLsnfYal+rRLmzebTVGtv+7QCVdicGSgw1/a/V8OAVnXbRViVJVnWmEZM8kX1KwZIB1GdlRHyxtjk1Ddn6RpXfH6jBYciA7TqUjyzD83b1xoKV3h8ro1zIEHq7OOJ9xiZV0ySGqdvdvFWrp3aEyJC2jWYgPCktKse6YPqeM/sZgyYGsMAyvyigGh7+ti5e7Cwa00lfSZc4A2SttZDumeRCrdlsZScvQEr05FXclBksOQqZ22DjXumkrg1axsS7ZfdVuNu+25mreMrJUUFxi6d2xKgyWHMShC1m4kJkPLzcXNeVD1kdL8t4Xn4Hk7HxL7w6RSWXmFWHXGX3nAFbttk5dowIQZliZuyWWK3PLYrDkYFNwg1qHwtPNxdK7QxUI8/NElyh/dXqtIQmWyF6sO56MklId2oT7IjqIVbutdWXu8Pb6UT/pKEB/Y7DkaCUDDHPSZO2r4vhGRfY5Baf1QyTrNKRNmHHVIlfm/o3BkgM4k5areo9JUjeHv20jb2lTbAryi5gzQPZTtmSdYaSCwZJ169siWK1WPHfxEuJSci29O1aDwZID0BKGpe9PgLe7pXeHqtAu0heNAryQX1SKzYa2NES2btfpi6psSZCPO7pGs2yJNfPxcFWrFYUW4BKDJYfAb3S2tXxX++bNVXFkL9Yc1acBDG4TyrIlNmCwYSqO9Zb+xpElB2icu/1kern/AGQbeUtSk6a0VGfp3SGqM5YMsC1D2ugLhkqvSvkMIQZLdm9rXJqqyBod5IUWoT6W3h2qBhkC93F3QXJ2AQ5eyOQxI5t2MiUHJ1Nz4ebiZCy8StZNKnlLs/WiEh3TAWo7sjR06FBkZGRccXlWVpa6jqxvua4Y3DqMjXNthIerCwYZvtmxmjfZSy+4mGbB8PV0s/TuUDXTAbRVcWvZ+qR2wdK6detQWFh4xeX5+fnYuHFjTe+OzEiWfa49mmLMFSDbm4pj3hLZulWGFidcBWdbtM8MyXnV6ZgO4FrdA3fgwAHj6cOHDyMxUV/kUJSUlGD58uVo1KiRqf9eVMemldKY1d3VGX1aBPNY2hD5VufsBBxOyFJ/Q1khR2SLVbt3ntZX7R7WljXebMk1zYPh6eaMhMx81dy7bYQfHFm1g6WuXbuqoTnZKppu8/LywgcffGDq/aM60FYyxDQLgrd7tf/UZAUCfdzRs0kQdpxOx5ojSfhnn6aW3iWiGlt/IkVV7ZaO9o2DWbXblkinh74tQtQ06tqjKQ4fLFV7Gu7UqVOIi4tTw3E7duxQ57Xt/PnzKmdp6tSp5v3rUY1o5eq5Cs42adMWK9lYl2yUrOgsW2yVbHNV3FrWW6r+yFKTJk3Uz9LSUvP9ZchkZLnnjlPp5V7wZFukR9PsP45iW5x++a4UiyOyFcWqard+dHs4W5zYJP0X7UPYfeYiMi8Vwd/LcRP0a/Xue+LECaxduxbJyclXBE/PP/+8qfaN6kCqP8uyz8ZB3moZKNmeFqEN1PLdM2l5qgSE1uCSyBZoH7CB3m7o1phVu22RNDxuGdYAsck52HQiFWM7R8JR1ThY+vzzz3H//fcjJCQEERER5Zajy2kGS9Zh3fEU46hS2b8R2ZZBrUOxYOsZVQKCwRLZktWGkgGyWEH6UpJtks+Q2OQcNRXHYKkGXnnlFbz66qt45plnzPfXoTqRvLL1huFv5ivZSbB0LEX9XRn4ku2VDOCIqC2TYPfzjafUe5B0FHB20MC3xnWWLl68iFtuucU8e0MmcSL575IBsvyTbFefMh3AT6WyAzjZhtOpuTiZkgtXZycMaM2q3basZ1N9R4HUnAIcupAFR1XjYEkCpRUrVphnb8ikjXP7NA+Gl7sLj6oNk5IPvZvpO4CvN0ytElk77bXas2kg/Fi126bJl+5+LfUBryOviqtWztL7779vPN2yZUvMnDkT27ZtQ6dOneDmVj47/uGHHzb9XlKNaCtQWLXbfqbiNsWmqg+gu/o1s/TuEF3VBkOwxDQA+zCkbRhWHE5SwdLDw1rBEVUrWHrnnXfKnW/QoAHWr1+vtrIkn4LBkmXlFBRj52l9yQC+UdkH6RP36u9HsO1kGvKLSlSxOCJrVVBcgi1xaer0wFYsW2IPBhvKz+yLz0B6biGCfNzhaKoVLEnhSbKtkgFNg1kywF5I9eNIf0/VdmD7qXQ10kRkrXadvohLRSUI9fVAu0hfS+8OmUCkvxfaRvjiaGI2Np5IwbiujtfarMY5S2QrU3D66s9k+2TEVguQtFWORNaerySvWa7etK+pOLHWUBLC0dS4ztLjjz9e4eXyn8LT01PlNI0bNw5BQfqkVKo/srRcS+5mvpJ9kQ+eRTvjsf64/H3bW3p3iK6ar8QRUPsrIfDxujgVDEu/P0ernVXjkaW9e/fiyy+/xGeffWbMW5JClXLZ6tWrVTAlAdPhw4fNssNz585F06ZNVWAWExOj+tRVRvZrwIABCAwMVNvw4cOvuP2UKVOMDYK1bdSoUbBFx5Ny1FSNB0sG2J1+rULUm1NcSi7i0/MsvTtEFUrMzFdTNVIHt79hBRXZh+6NA+Dr6YqLeUXYfy4DjqbGwZKMGknQceHCBezevVtt586dw7XXXos77rhDNdUdOHAgHnvsMZPv7Pfff6+CsVmzZmHPnj3o0qULRo4cqdquVGTdunVqn6Q1y9atWxEdHY0RI0aofSxLgqOEhATj9t1338GmSwa0CGYSsJ2R5dc9DC0jNpzgVBxZ96hSl6gABDpgErA9c3VxxkBDOsA6B5yKq3Gw9Oabb+Lll1+Gn5+f8TJ/f3+88MILeOONN+Dt7a1ankgQZWpz5szBtGnTcNddd6F9+/b45JNP1OPNmzevwtsvXLgQDzzwALp27Yq2bdviiy++UL3sZASsLA8PD9W6RdtkFMoWaTUwZLiU7HNVXNm8NCJrzlci+zPE8Nmy1gHfg2ocLGVmZlY4kpOSkoKsLH11z4CAABQWFsKU5P4kAJNRLY2zs7M6L6NG1ZGXl4eioqIr8qlkBCosLAxt2rRRfe/S0vTLXitTUFCgnmvZzdKy84vUKhTBfCX7pH0AbYlNRWFx+QbWRJZWXFKqVkqVDezJPt+D/jqfieTsfDiSWk3DTZ06FUuXLlXTb7LJ6bvvvhvjx49Xt5G8oNatW5t0R1NTU1FSUoLw8PJ9huR8YmJite5D+tk1bNiwXMAlU3ALFixQo02vv/66ysEaPXq0eqzKzJ49W42maZtM71na5tg0FJfq0CzEB02CfSy9O2QG7SP9ENLAHbmFJaqjO5E12X8uE1n5xfD3clPTcGR/Qn090DnK3yFX5tY4WPr0008xbNgw3H777WjSpIna5LRcJtNiQpvysib/+c9/sGjRIhXYSXK4Rvb9hhtuUNXIJdj79ddfsXPnTjXaVJkZM2aoETZti4+Ph6VxFZz9kwaWWpE/tj4ha6O9JvsbFiOQfRpsmIpztHSAGgdLUr1bVpnJVJWsjJNNTsvqOB8f/YiG5AjJZkohISFwcXFBUpK+k7VGzkueUVXeeustFSxJT7vOnTtXedvmzZurx4qNja30NpLjJDlbZTfLlwxgfSVHoE1vMFgia8N8JccwxPAeJAtNikocJx2g1kUpJWiSwEM2OW1u7u7u6NGjR7nkbC1Zu0+fPpX+niSdS0L68uXL0bNnz6s+jkwrSvAXGRkJW3EsKRuJWfnwdHNGjKHpKtmnAa2k0B9wJCELSVmOlTNA1ktaYBwwLCdncrd96xwVoNqdZOcXY48DpQNUqyjlTTfdhPnz56sRFDldlSVLlsBcpGzA5MmTVdDTu3dvvPvuu8jNzVWr48SkSZPQqFEjlVMkJAdJVuZ9++23qjaTltskwZ1sOTk5ePHFFzFhwgQ1OhUXF4enn35a1YmSkgS2Yu1R/ahS3xYhLBlg5+RNSt6s9sdnqG/yt/a0fL4ckTR61umgWmKE+/2d5kD2x8VZ31Fg6d7zalVcTPNgOIJqjSxJErNWtr5sYnNFmznddtttakpNAiCZ5tu3b58aMdKSvs+ePavqJGk+/vhjtYru5ptvViNF2ib3IWRa78CBAypnSRLSJUldRq82btyoptpsBfOVHIux9YkhR4TI0rRkX44qOYZBhvcgbfWjI3DSScIL1YmUDpBAUZK96zt/KSu/CN1eWqnKz294aggaB3vX6+NT/ZOVcBM+3qJWHe1+brgqFkdkKaWlOvR+bTVScwrw7T0x6MvK3XYvNacAPV9ZpU7v/PdwtUrO3j+/a/UuW1xcjFWrVqmVcdnZ2eoyqegt01pUvzafSFWBUvNQHwZKDqJLlL8KlDIvSduBTEvvDjm4I4lZ6sPT290FPZraZkFfqpmQBh7o0FAfWGyKdYzRpRoHS2fOnFHL7KXe0oMPPqiKUWr5QU8++aQ59pGqYFwF15pVux2FjCTJ8mzBqTiyNO012LdFMDxcXSy9O1RPBmpTccdTHeKY1zhYeuSRR1SC9cWLF+Hl5WW8/MYbb7yijQiZl8ygan3CWDHXQfOWDC1uiCzdD0778CTHMMDwhW3DiVQ1FWvvqrUarixJft6yZYtayl+WrDa7vEEtmVdcSg4SMvPh7sqSAY5msOGD6cD5TKTlFCC4ge3mDJDtyikoNrZZYnK3Y+nZJEhNvcoU7NHEbLQ3TMvZqxqPLElto4pagUh9Il9fX1PtF1XDBsPwp9RW8nTj8LcjCfPzRLtIP7VcW5ZtE1mC9CmUNktNg73ZZsnBuLs6o4+hbIA2w2HPahwsjRgxQtU30khJAUnsnjVrFsaMGWPq/aMqaC9QrQUGOepUnP2/UZF1YtVuxzZAm4pzgDImNQ6W3n77bWzevBnt27dHfn4+/vGPfxin4CTJm+pHQXEJtp1MU6cHtNa/YMkxgyUJmh0hZ4CsL2dSC5aYr+SYBhreg2QqNq+wGPasxjlLUVFR2L9/v2pKKwUdZVRJijlOnDixXMI3mZe8OPOLShHm64E24Zz+dEQ9mgTCR+UMFOLQhSx0MnQDJ6oPp1Jzce7iJbi7OOMaB6niTOU1C/FBVKCXeh1sP5mOIW3td1W2a61+ydUVd955p+n3hmo8BafvFcYO346aMyAFAFceTsL648kMlqheaaNKvZoFwsejVh8lZOOcnJzUZ9B3O86q14M9B0s1noZr3Lix6sH25Zdf4uTJk+bZK6p2cvdATsE5tMGGDuCst0T1jflKJAYZPoPsvfVJjYOl1157DZ6enio/SRrORkdHq1Gmzz//HCdOnDDPXlI5ydn5quu86M/WAg5NS+7fczZDVfQmqg/5RX/nTDJfybH1aRGimuvGpeTifMYl2KsaB0sSGH322Wc4fvy4Sup+88031eUPPPAA2rZta459pMtsNiwV79jIj/V1HFx0kDdahPqoljeyjJuoPuw8na5yJsP9mDPp6Py93NA1OkCd3mjHq+JqNdGcl5eHTZs2Yd26dVi7di327t2Ljh07YvDgwabfQ6p8Co4lA8iQtybf6qSS7uhOkTwmZHZauQpZkcmcSRrYKlQ1+JZc2tt7N7bLA1LjkaW+ffsiODgYzz77rCodID8TEhJUwPTOO++YZy/JSJaIbzyRavyQJDKWEDieopZzE9VfvpL9JvRS9Q0w5C1tOpGK4pJSuzx0NQ6Wjh49Ch8fHzXlJlu7du0QGMhO0xbp8N2Ex52AmOZBcHNxUvkCp9PyeEjIrC5kXMKJ5Bw4OzFnkvS6RAWo6bis/GLVgske1ThYSktLw5o1a3DNNdfgzz//RL9+/dCoUSNVnFKSvMm8tFElKTMvS8eJvN1dVZ8m/evDfnMGyDpo1ZolT8Xf283Su0NWwMXZybjYyF6redf401bmpzt37oyHH34YP/74I/744w9ce+21WLx4Me677z7z7CUZscM3VTUMruWzEZn7CxtXwZEjtT6pcbC0Z88ezJkzBzfccIPKXerTp4+q5P3QQw9hyZIl5tlLUqScvNbhW3thEpVN9t8al4oiO80ZIMuTVZda42bmTFJZAwy5k/vi7bOMSY1Xw/Xu3RvdunXDoEGDMG3aNAwcOBD+/myzUB+2n0pHYUkpGgV4qTLzRJr2kX4I9nFHWm4h9py5iBi2nyAz+Ot8pvog9PV0RRe216Ey5HNJypjIylwpY2JvK3NrHCylp6fDz8/PPHtD1Z6C43JdKstZcgZaheDnfRfUNAmDJTIHrY5OvxYhcHVhziSVJ59N9lrGpMavdgZKVpArwCk4qoA2LcIkbzL3e5CWI0dU1kA7LmPCrwY2tFw31rBcV5qnEl1Oy2OTpbsXcwt5gMiksvOLsOesPmeSBXGpIjHNguDu4qzKmJxMzYU9YbBkI7TRArVc14vLdelK4X6eaBPuC/lCtzmOq+LItLadTEdxqQ5Ng71Vmx2iisqY9GoWaJetTxgs2VqLE8MwJ1FVo0sbWUKAzPSFje9BVJ10AMlbsicMlmwAl+tSTZfvygebveUMkGWxzRLVrIxJGgqKS+Cwq+FKSkowf/58rF69GsnJySgtLV/TRap7k2kdOKevW8HlunQ1vZsGqcruFzLzEZeSg5ZhvjxoVGfx6Xk4lZoLV2cnXNNcXy2eqCLtIn0R0sBDteWS5rp9W4Q45sjSI488ojYJmjp27IguXbqU28h83+iknDyX61JVvNxdVJKlYDVvMhXpJi+6Nw6ErydzJqlyUtZGW7FtT+9BNR5ZWrRoEX744QeMGTPGPHtEleYKsGIuVTdvSQJsed1M7d+MB43qTMuBY+cAqg7Ja1uy97x6D3p2dFs45MiSu7s7WrZsaZ69oStkqeW6Geo036ioOrSgWlYv2VPOAFlGcUmpcXWllhNHVBUpkCsOXchCSnYBHDJYeuKJJ/Dee+8xebSeSJKcJHg3D/Hhcl2qlrYR+pyBS0UlKmeAqC72n8tEdn6xKlnSqRFbW9HVyftPh4b6Th+bYlMccxpu06ZNWLt2Lf744w906NABbm7l56/ZTNd8LU6IapIzoB8GT7WbBEuybBqA5Ey6SFVcomqQzywZWZK8pRu7RcHhRpYCAgJw4403qka6ISEhqolu2Y3MtVyXH3hUfVo7CrY+Ib4HkUVrvp1IRWmpzvFGlr766itY0ty5c/Hmm28iMTFRrb774IMP0Lt370pvv3jxYsycOROnT59Gq1at8Prrr5dLTpdaNLNmzcLnn3+OjIwM9OvXDx9//LG6raWdTs3F2fQ8uLnIct1gS+8O2ZD+LfUjkQfPZ6klvDIsTlRTUrJkX3xGuTwUouro0SQQ3u4u6v3naGI22hum5RyuKGVKSoqakpNNTteH77//Ho8//rgKbvbs2aOCpZEjR6p6TxXZsmUL7rjjDtx9993Yu3cvxo8fr7aDBw8ab/PGG2/g/fffxyeffILt27fDx8dH3Wd+fj4sTRsVkBedj0eN41pyYKG+HmgfqX9z2hxrP8t3yUI5k6E+iApkixOqPg9XF+OXfHsY4a5xsJSbm4upU6ciMjISAwcOVFvDhg1VQJKXlwdzmjNnDqZNm4a77roL7du3VwGOt7c35s2bV+HtJRF91KhReOqpp9CuXTu8/PLL6N69Oz788EPjqNK7776L5557DuPGjUPnzp2xYMECXLhwAT/99FOl+1FQUICsrKxymzlo5eJZMoDqMhVnT7VOyEItTgwrLIlqOxXncMGSjOysX78ev/zyi5q2ku3nn39Wl8lKOXMpLCzE7t27MXz4cONlzs7O6vzWrVsr/B25vOzthYwaabc/deqUms4rexvJu4qJian0PsXs2bPL5WlFR0fDHMt1t8WlqdODmNxNtaB9wLH1CdUWcyapLgYaPrt2nE7HpcISxwqW/ve//+HLL7/E6NGj4efnpzbJAZKcnx9//NE8ewkgNTVVVQ0PDw8vd7mcl4CnInJ5VbfXftbkPsWMGTOQmZlp3OLj42FqUql71RODMOfWLsbpFKKakOlbTzdnJGcX4HhSDg8e1ciZNOZMUt00D/FBowAvFBaXYvsp/Zd/hwmWZKrt8uBChIWFmX0azlp4eHgYA0VtM4dwP0/c1D0KzlyuS7Xg6SatT+wnZ4Dql5YGIC1OmDNJtS1jYi9TcTUOlvr06aMSrMsmQF+6dAkvvviius5cpEyBi4sLkpKSyl0u5yMiIir8Hbm8qttrP2tyn0S2RHuj0j74iKprI2u8kUneg0L170GG15PDBEuSNL1582ZERUVh2LBhapOcHVl5JteZi7RZ6dGjB1avXm28rLS0VJ2vLEiTy8veXqxcudJ4+2bNmqmgqOxtJFlbVsWZM/Ajqi9avtv2k2nIL7LtnAGqP0UlpWolnGCNN6qLfi2DIZMjJ5JzkJB5CbaqxuvRO3bsiBMnTmDhwoU4evSoukyW50+cOBFeXl4wJ0kunzx5Mnr27KlqK8lKNlmdJ6vjxKRJk9CoUSOVgC0eeeQRVTzz7bffxtixY1UT4F27duGzzz4zDhE++uijeOWVV1RdJQmepCaTrO6TEgNEtq5lWANE+HkiMSsfO0+nc2UlVcv++AxkFxQj0NsNHRuy2DDVXoC3OzpHBah6XTIVd2tP0y+Iqg+1Kt4jy/VlCX99u+2221RNp+eff14lYHft2hXLly835lCdPXtWrZDT9O3bF99++60qDfB///d/KiCSkgAS8GmefvppFXDde++9amVf//791X16enrW+/MjMlfOwOLd59QbFctQUHVoUyb9W4UyZ5LqTNovSbAkrytbDZacdFJs6CqWLVumVr9JHzg5XZUbbrgBjkam7qSEgKyMM1eyN1FtLdt/AQ9/t1c12F3+6EAeSLqq8XM3qw+3N27ubLMfbmQ9dp1Ox82fbEWAtxt2P3etVfUYrO7nd7VGlmRKSkZyZMVbVdNT8i1WlvcTkfWQBqhOTlAtB5Kz8hHmx1FTqlxGXiEOnNO3OGG+EplCl+gA+Hq4IiOvCIcuZKppObtM8JZEagmUtNOVbQyUiKxPkI+7Me/E1pfvkvltiUuD9D1tFdYAkf7mzUMlx+Dm4ow+LYJt+j2oxqvhpB2ItPuoqMK2XEdE1megofUJ6y3R1WivEea3kSkNMKzMXW+jJQRqHCzJyjOZ27tcdna2cVUaEVkX7YNPvtWVyrABUQUkhVXrJaj1FiQyhUGG96A9Zy4ip6AYdh8syX8myU263Llz51SSFBFZH1WF2d0FabmFOJxgnsbPZPtOpebifMYluLs4I6ZZkKV3h+xI42BvNAn2RnGpztj31C5LB3Tr1k0FSbJJIUpX179/VXKVpCntqFGjzLWfRFQH7q6SMxCCVUeS1DB4x0b8YkNX0vJJejYNhLd7rSrLEFVKFgycSTuLDSdSMLz9lW3TrFm1/zdoq+D27duHkSNHokGDBuWqazdt2hQTJkwwz14SUZ0Naq0PlqTWyYNDWvKIUqX1lZivROYwoFUovtl21iaTvKsdLEk/OCFBkRSHZNFGItsy0JBguduQM9DAgyMH9LeC4hJsPZlWbkEAkSn1bRGsaizJdG98eh6ig7ztN2dJ2o0wUCKyPU2CfYw5A1rfLyLN7tMXkVdYgpAGHmgXweK6ZHq+nm7o3lhfY8nWRpdqHCxJftJbb72lerNJE9qgoKByGxFZr4F20gGcTG+9oWSAjCo5W1GFZbIvA2z0PajGwdKLL76IOXPmqKk4KSEgzW1vuukm1ZPthRdeMM9eEpFJp+IkwZKoLK1kwCDDa4TInO9Bm+NSUVxSar/B0sKFC/H555/jiSeeUCvi7rjjDnzxxReque22bdvMs5dEZBJSRdfV2Qln0vJwJi2XR5UUaYNzJCFLtcWR9jhE5tKpkT/8vdyQnV+M/eeurNloN8GS9Ijr1KmTOi0r4rQClddddx1+++030+8hEZmMJHX3aBJok8PgZD4bDPkj8kEW3MCDh5rMxsXZyRiQ29J7UI2DpaioKCQkJKjTLVq0wIoVK9TpnTt3wsOD/8mIbGUYfL1h2oVI+9DSctqIzGlAK9trv1TjYOnGG2/E6tWr1emHHnoIM2fORKtWrTBp0iRMnTrVHPtIRCak5aRsjUtFYbHt5AyQeZSU6owfWlogTVQffeL2xWcg81KRTRzsGhda+c9//mM8LUnejRs3xtatW1XAdP3115t6/4jIxNpH+iHYx121Ptlz9iKuaa7vBk6O6eD5TFzMK1JTtN0My7qJzKlRgBdahPogLiUXW2JTMbpTpP2NLF2uT58+akUcAyUi2yDLwrVhcFvKGSDz0F4D/VoGw82lzh8JRDUrIWAj9ZaqNbK0bNmyat/hDTfcUJf9IaJ6INMtP+27oEoIPD2qLY+5A9PKSHAKjurTwNYhmL/ltArWdTqd6jtr88GS1hfuauTJStFKIrKNb3UHz2chNadAVW0mx5OVX4Q9ZzPUaSZ3U326prmMZDrhfMYlnE7LQ7MQH6v+A1RrzLW0tLRaGwMlItsQ6uuhcpfEJhsZBifTk3wRSfBuHupjU326yPZ5u7uiZ5Mgm0kHqNMEdX5+vun2hIgsU83bBt6oyDy08hEcVSJLGNDadkoI1Ko33Msvv4xGjRqpopQnT55Ul0sJgS+//NIc+0hEZqB1lpcEy9JSHY+xg5E8ES1QZosTsoSBhnSALXFpVl/GpMbB0quvvor58+fjjTfegLu7u/Hyjh07qrYnRGQbpJK3l5uLylk6kphl6d2henYyNVfli7i7OCOmOZugU/1rH+mHkAbuyCsswa7T6fYVLC1YsACfffYZJk6cCBcXF+PlXbp0wdGjR029f0RkJh6uLqpXXNkmquQ41h/Tjyr1bhak8keILFHGZKAhHWCdlacD1DhYOn/+PFq2bHnF5ZLgXVRkG5U4iUhvIOstOay/SwawcS5ZzuA2YernumPJ9hUstW/fHhs3brzi8h9//BHdunUz1X4RUT3QvtXtOpOO3IJiHnMHkV9Ugm0n09Rp1lciSxrQMgTOTsDxpBxcyLhktX+MGo+9Pv/885g8ebIaYZLRpCVLluDYsWNqeu7XX381z14SkVlIbZOoQC+cu3hJfXgOaxfOI+0Adp2+iPyiUoT7eaBNuK+ld4ccWKCPO7pEB2Dv2QysP56CO3o3hl2MLI0bNw6//PILVq1aBR8fHxU8HTlyRF127bXXmmcvicgspJAsSwg4nvXHk42rkay9cjLZv8GtrX8qrkbBUnFxMV566SU0a9YMK1euRHJyMvLy8rBp0yaMGDHCfHtJRGZfvmsrPZqo7rSEfk7BkTUY3Eb/HrQ51npLCNQoWHJ1dVUlAyRoIiL70LdlMFycnXAqNRfx6XmW3h0ys8TMfBxLyoYMKPVvyeRusrxOjfwR7OOOnIJi7D5zEXYxDTds2DCsX7/ePHtDRPXOz9MN3RsHqNOSM0D2TStE2SUqQOWLEFlTCYH1VvoeVONgafTo0Xj22Wfx5JNP4rvvvsOyZcvKbeaSnp6uajv5+fkhICAAd999N3Jycqq8/UMPPYQ2bdrAy8sLjRs3xsMPP4zMzMxyt5P5+su3RYsWme15EFnzVJwttB2gullvLBmg/5sTWdNU3DorzVuq8Wq4Bx54QP2cM2fOFddJoGGuZroSKCUkJKhcKanndNddd+Hee+/Ft99+W+HtL1y4oLa33npLlTs4c+YM7rvvPnWZlDko66uvvsKoUaOM5yUYI3Ik8sH59srj2BKbhqKSUri51KltJFkpaZqrNU4exPpKZEUGqMUGwNHEbDVVHOHvCZsOlqRcQH2T1XbLly/Hzp070bNnT3XZBx98gDFjxqhgqGHDhlf8jrRf+d///mc836JFC9Wq5c4771Q5V5J/VTY4ioiIqPb+FBQUqE2TlcVWEWTbOjbyR6C3Gy7mFWFffAZ6NWX7C3t04FwGMi8Vwc/TVU3DEVmLIB93dI4KwP54KSGQjNt6WVcJgRp9fZQRHQkyDh48iPq0detWFdBogZIYPnw4nJ2dsX379mrfj0zByTRe2UBJPPjggwgJCUHv3r0xb9481WCyKrNnz4a/v79xi46OrsWzIrIekuDdX1sVZ6U5A1R3Wj5I/1YhcOXoIVmZwVrrE0MrHpsNltzc3FTuj7mm2iqTmJiIsDB9HQaNBDxBQUHquupITU3Fyy+/rKbuypJSCD/88IOa3pswYYKaZpRRq6rMmDFDBV7aFh8fX4tnRWRdBlnxGxWZhhYIazlqRNaYt7TpRKpKB7AmNU5M+Pe//43/+7//UwnUdSWJ4hUlWJfdTNGcV6bJxo4dq3KXXnjhhXLXzZw5E/369VOtWp555hk8/fTTePPNN6u8Pw8PDzVCVXYjspdg6a/zmUjOyrf07pCJZRqmWAWTu8kadZYVmt5uyC4oVhW9bTpn6cMPP0RsbKzKE2rSpImq4l3Wnj17qn1fTzzxBKZMmVLlbZo3b67yiaQAZlmSdyQB29VyjbKzs1Xytq+vL5YuXapGx6oSExOjRqAkJ0mCIiJHEerrgS5R/th/LlONLt3ai9PL9mRzXCpKdUCrsAZoGOBl6d0hqjAdQBK9l+2/oFbF9W4WZLvB0vjx40324KGhoWq7mj59+iAjIwO7d+9Gjx491GVr1qxRyeYS3FQ1ojRy5EgV9EhZA0/Pq2fX79u3D4GBgQyUyGE7gEuwtOZoMoMlO7PeML3KUSWy9qm4ZSpYSsHTo9rCZoOlWbNmob61a9dOjQ5NmzYNn3zyiUo0nz59Om6//XbjSjhp7CsFM6WhryRqS6AkLVikHcs333yjzmur1iRAc3FxUf3skpKScM0116hASvKWXnvtNVVDisgRDW0bhvdWn8Cm2FTVdsDdlSUE7IEsWlmn9YNjfSWyYgMNr8/DCVkqHSDMz9M2gyWNjPLIkn7RoUMHlfNjTgsXLlQBkgREsgpOkrHff/994/USQB07dkwFR9p0oLZSrmXLluXu69SpU2jatKmakps7dy4ee+wx9WYit5P6URKUETlq24GQBh5IzSnAztPp6Md2GHbh0IUsJGUVwNvdBTFWNLVBdDl5/+kc5Y8D5zLV6s1bekbbZrAkuUMyorNu3Tpj8UaZIhsyZIiqfF2dabXakJVvlRWgFBL8lF3yP3jw4KuWAJDRqrLFKIkcnbQdkGHwH3efU1NxDJbsw9qj+lEl+Xt6urlYeneIrlpCQIKldVYULNV4jF1aiEjS9KFDh1SCtWxSd0mmuKSdCBHZ/lScWGulbQeo5lYbgqVhhr8tkTUb1Eb/Ot14PAXFVlJCoMbBklTS/uijj1QekUaW5Mt01h9//GHq/SOieqYKFjo74WRKLs6k5fL42ziZUt1/Tr8MewiDJbIBXaMD4O/lhqz8YmO5C5sLlmQFWkXL7+UyS7RCISLT8vN0M7Y7kak4sm2yqkgyEjo28kO4lSTLEl29hECIVRXJrXGwNHToUDzyyCOqIa1GVqJJkrQkXxOR7RvSVp97yGDJ9q05mqR+Dm0bbuldIapRGROhreK0uWBJilJKfpIkVEtzWtmaNWumLrtamxAisq28pe0n05FbUGzp3aFakvIPG4+nlvubEtlSR4GD57OQkv1343qbWQ0nTWNlWf6qVauMrUgkf0ka2xKRfWgR2gDRQV6IT7+EzbGpGNGh6kr5ZJ12nU5XrSNCGrijcyN/S+8OUY06CsjUsQRL0tNwQo8o2FydJenZdu2116qNiOyP/B8f2iYMX289o1bFMViyTdo0qkxpSFkIIlsyuHWYCpbWWUGwVO1pOGkvIqvetCrYZWVmZqrClBs3bjT1/hGRhWgrp9YeTblqzTKy7mCJJQPIFg1qo5+K23giBSXS2NAWgqV3331XVbb28/O74jp/f3/861//UtWvicg+XNM8GF5uLkjMyseRhGxL7w7V0KnUXJxMzYWbi5MqB0Fka7pFB8DP0xUZeUUWLyFQ7WBp//79VVa7lj5s0gKFiOyDVHru1zJYnWaBStsdVZLO7b6eV5Z7IbJ2ri7OGNBKP7okrU9sIliShrMV1VfSuLq6IiXFOuohEJFpl++yhIDtlgwYYvgbEtnyVNx6C3cUqHaw1KhRI9XWpDIHDhxAZGSkqfaLiKwob2nv2Yu4mFto6d2hasrOL8KOU+nq9LB2rK9Ett0nThw4n4m0nALrD5bGjBmDmTNnIj8//4rrLl26hFmzZuG6664z9f4RkQU1CvBC2whfSG6lpYfBqfo2nUhFUYkOzUN80CzEh4eObFaYnyfaRfqhob8XzqbnWX/pgOeeew5LlixB69atMX36dLRp00ZdLrWWpC9cSUkJ/v3vf5tzX4nIQqNLRxOz1VTc+G6N+Dewoca57AVH9uDbe2IQ4O2mSppYfbAUHh6OLVu24P7778eMGTOMS4ll50eOHKkCJrkNEdkXqfz88bo4NbIky3elbxNZr9JSHdYZ8jtYMoDsQaCPu6V3oWZFKZs0aYLff/8dFy9eRGxsrAqYWrVqhcDAQPPtIRFZfPmudADPvFSkcpd6GprsknWS3I7UnEI08HDl34rIUr3hhARHvXr1Qu/evRkoETnA8l2tTxNXxVk/7W80sHUI3F1r9RZPRJfh/yQiuqohbRks2QqWDCAyPQZLRHRVg1qHQXIrJdH7QsYlHjErlZSVr3ppyd9Kq5FFRHXHYImIrirIx13lLglW87Zeaw1TcJ2jAlTXdiIyDQZLRFTtVXFlP5DJeksGcBUckWkxWCKiatFq9myOTUN+UQmPmpWRv8nm2NRygS0RmQaDJSKqlvaRfojw88SlohJsN7TSIOshf5O8whKE+3mgQ0M/S+8OkV1hsERE1SIFaLVVcZyKsz7a30RGlSxZ6ZjIHjFYIqJq01ZYrT6aZKziT5Ynfwv5m4ghXAVHZHIMloio2ga0CoGnmzPi0y/hcEIWj5yViE3OUX8TKULZr2WIpXeHyO4wWCKiavN2dzVW815+MJFHzsqqdl/TPBg+HjXqYkVE1cBgiYhqZFTHCPWTwZL1YMkAIvNisERENTK0bTjcXJxwIjlHTf+QZaXnFmL3mYuGvw1LBhCZA4MlIqoRfy83Y17M8oMJPHoWtvJwIkpKdaq0Q3SQt6V3h8guMVgiohob1cEwFXeIeUuW9vtf+r/BmE76vwkROXCwlJ6ejokTJ8LPzw8BAQG4++67kZNT9RTA4MGDVb2Rstt9991X7jZnz57F2LFj4e3tjbCwMDz11FMoLi4287Mhsm3Xtg+HsxNU09b49DxL747DyswrMlbtHt0p0tK7Q2S3bCZYkkDp0KFDWLlyJX799Vds2LAB995771V/b9q0aUhISDBub7zxhvG6kpISFSgVFhZiy5Yt+PrrrzF//nw8//zzZn42RLYtuIEHYpoFq9N/cnTJYlYeSUJxqQ5twn3RIrSB5XaEyM7ZRLB05MgRLF++HF988QViYmLQv39/fPDBB1i0aBEuXLhQ5e/KiFFERIRxk5EpzYoVK3D48GF888036Nq1K0aPHo2XX34Zc+fOVQFUZQoKCpCVlVVuI3LUVXF/sISAxfzxlz5nbDSn4IjMyiaCpa1bt6qpt549exovGz58OJydnbF9+/Yqf3fhwoUICQlBx44dMWPGDOTl5ZW7306dOiE8PNx42ciRI1XwI6NYlZk9ezb8/f2NW3R0dJ2fI5GtGWnIW5KVWElZ+ZbeHYeTlV+EjSf0U3BjOAVHZFY2ESwlJiaqfKKyXF1dERQUpK6rzD/+8Q81arR27VoVKP33v//FnXfeWe5+ywZKQjtf1f3KfWVmZhq3+Pj4Ojw7ItsU4e+Jbo0D1OkVnIqrd2uOJKOwpBQtQn3QKoxTcETmZNFSr88++yxef/31q07B1VbZnCYZQYqMjMSwYcMQFxeHFi1a1Pp+PTw81Ebk6EZ3jMDesxlqKu6ffZpaenccyu+GKTgZVWLjXCI7DpaeeOIJTJkypcrbNG/eXOUaJSfry/lrZMWarJCT66pL8p1EbGysCpbkd3fs2FHuNklJ+maUNblfIkc1qkMkXvv9KLafSlfFEYN83C29Sw4hp6AY646nqNOjO3IVHJFdB0uhoaFqu5o+ffogIyMDu3fvRo8ePdRla9asQWlpqTEAqo59+/apnzLCpN3vq6++qgIxbZpPVttJEnj79u1r+ayIHEfjYG9VDFGa6q46nIRbezF/rz6sPZqMwuJSNA32RrtI33p5TCJHZhM5S+3atcOoUaNUGQAZCdq8eTOmT5+O22+/HQ0bNlS3OX/+PNq2bWscKZKpNlnZJgHW6dOnsWzZMkyaNAkDBw5E586d1W1GjBihgqJ//vOf2L9/P/78808899xzePDBBznNRlSDqTjxB6t51xvtWEttJU7BEZmfTQRL2qo2CYYk52jMmDGqfMBnn31mvL6oqAjHjh0zrnZzd3fHqlWrVEAkvydTfhMmTMAvv/xi/B0XFxdVs0l+yiiTJH9LQPXSSy9Z5DkS2SJt2frm2DS1QovMK6+wGGuP6qfgxnAKjqheOOl0Ol39PJT9klIDUkJAVsaVreNE5CiGvb0OcSm5eO/2rhjXtZGld8fuayvdv3APogK9sPHpIRxZIqqHz2+bGVkiIuulJRn/YehTRubzu6EIKFfBEdUfBktEZLJq3uuOJ+NSYQmPqJnkF5VgzZGkcrliRGR+DJaIqM46NPRT00L5RaVYf7x8mQ8ynQ3HU5BbWIKG/p7oGq0vCEpE5sdgiYjqTFZk/b0qjlNx5qId21EduQqOqD4xWCIik07FSRuOgmJOxZmaHNNV2hQcG+cS1SsGS0RkEt2iAxHu54HsgmJsiU3jUTUxOabZ+cUI8/VAj8aBPL5E9YjBEhGZ5s3E2QkjO+hHl5ZzKs5sveBkBE+ONRHVHwZLRGQyowzB0orDiSguKeWRNZGiklKsOKytgmMvOKL6xmCJiEymd7MgBHq74WJeEXacSueRNZGtcWnIvFSEYB93dYyJqH4xWCIik3F1cca17cPV6eWHuCrO1L3gRnaMgAun4IjqHYMlIjIpbZpI8pZKS9lNqa5kOvPPQ/opOPaCI7IMBktEZFJ9WwbD19MVydkF2HaSq+LqSqYz03ML1fRmTHNOwRFZAoMlIjIpD1cXXN+loTq9ePc5Ht06+t0wBTeifQTcXPiWTWQJ/J9HRCZ3S48oY65Ndn4Rj3AtlZTqsPwgC1ESWRqDJSIyOelb1iLUR/WK++2AfmSEarcKLjWnAH6erujbIoSHkMhCGCwRkVl6xd3SM1qd/pFTcbX2/a549VOmNd1d+XZNZCn830dEZnFTt0aQVe67zlzEyZQcHuUaysgrxJ+G8gu39dIHnkRkGQyWiMgswvw8Mah1qDr9vz1M9K6pn/aeR2FxKdpG+KJTI38z/IWIqLoYLBGR2WhTcf/bfV4lK1P1/bDrnHFUSaY1ichyGCwRkdkMaxeGAG83JGblY1NsKo90NR08n4nDCVlwd3HG+K6NeNyILIzBEhGZtebSOEPNJSZ6V9/3O/WJ3SM6hCPQx91Mfx0iqi4GS0RkVjf30E/FSbJyZh5rLl1NflEJftp3Xp1mYjeRdWCwRERm1bGRn0pSlmTlXw5c4NG+Cumpl51fjEYBXujH2kpEVoHBEhGZlSQn32yo6M32J9WfgrulZxScpfYCEVkcgyUiMrvx3RrB1dkJ++MzcCIpm0e8EmfScrH1ZBpk8Zu2kpCILI/BEhGZXUgDDwxuE6ZOM9G7cosN5QL6twxR03BEZB0YLBFRvZBpJbFk73kUl5TyqF9G6lBpgeStHFUisioMloioXgxtG4ZgH3ekZBdgw4kUHvXLyDGRelRSl0pKBhCR9WCwRET1ws3FGeMMBRY5FXelHwyJ3VKEUupTEZH1YLBERPU+FbfqcDIu5hbyyBuk5RRg1ZEkdZq1lYisD4MlIqo37SL90KGhHwpLSvGzofAiAUv3nkdRiQ6do/zVMSIi62IzwVJ6ejomTpwIPz8/BAQE4O6770ZOTk6ltz99+rSq71LRtnjxYuPtKrp+0aJF9fSsiBzPLYaaSz/u0SczOzqdTmesrcTEbiLrZDPBkgRKhw4dwsqVK/Hrr79iw4YNuPfeeyu9fXR0NBISEsptL774Iho0aIDRo0eXu+1XX31V7nbjx4+vh2dE5Jhu6NoIbi5OOHg+C0cSsuDo9krtqeQceLo544au+j56RGRdXGEDjhw5guXLl2Pnzp3o2bOnuuyDDz7AmDFj8NZbb6FhwyvfYFxcXBAREVHusqVLl+LWW29VAVNZMlJ1+W2rUlBQoDZNVhbf8ImqK8jHHcPbheOPg4kq0Xvmde0d+uBpid1jOkbCz9PN0rtDRLY6srR161YV0GiBkhg+fDicnZ2xffv2at3H7t27sW/fPjV9d7kHH3wQISEh6N27N+bNm6eGxasye/Zs+Pv7GzcZxSKi6tPan/ykcnUct+ZSbkExftmv75d3ay++jxBZK5sIlhITExEWpq/+q3F1dUVQUJC6rjq+/PJLtGvXDn379i13+UsvvYQffvhBTe9NmDABDzzwgBq1qsqMGTOQmZlp3OLj9d8Miah6BrUOVVW903IL1QiTo/rtrwTkFpagabA3YpoFWXp3iMgag6Vnn3220iRsbTt69GidH+fSpUv49ttvKxxVmjlzJvr164du3brhmWeewdNPP40333yzyvvz8PBQieZlNyKqPlcXZ/zzmibq9Mfr4q46mmuvFu/SmuZGq/c7IrJOFs1ZeuKJJzBlypQqb9O8eXOVT5ScnFzu8uLiYrVCrjq5Rj/++CPy8vIwadKkq942JiYGL7/8sspJkqCIiMxjct8m+GxDnEryXncsBUPalh89tndxKTnYefoinJ3+npYkIutk0WApNDRUbVfTp08fZGRkqLyjHj16qMvWrFmD0tJSFdxUZwruhhtuqNZjSV5TYGAgAyUiMwvwdsfEayRgOokP18ZicJtQhxpd+WLjSfVzSJswhPt5Wnp3iMjWc5Yk12jUqFGYNm0aduzYgc2bN2P69Om4/fbbjSvhzp8/j7Zt26rry4qNjVVlBu65554r7veXX37BF198gYMHD6rbffzxx3jttdfw0EMP1dtzI3Jk9/RvBncXZ+w+cxE7TqXDUZy7mGds+XLf4BaW3h0isodgSSxcuFAFQ8OGDVMlA/r374/PPvvMeH1RURGOHTumptvKktVtUVFRGDFixBX36ebmhrlz56qRq65du+LTTz/FnDlzMGvWrHp5TkSOLszP09gCZe66ODiKT9bHqYrdfVsEo1dTJnYTWTsnnaNmVpqQ1FmSEgKyMo7J3kQ1czYtD0PeXoeSUh1+md4fnaL87foQJmRewqA31qmWL4vuvQbXNA+29C4ROaysan5+28zIEhHZp8bB3rihi346/aN1sbB3n6yLU4FS72ZBDJSIbASDJSKyuPsNeTvLDyUiNrnyno+2LikrH98ZKnY/OqyVpXeHiKqJwRIRWVzrcF+MaB8OSQqQfB57Jc+tsLgUPZsEok8LTr8R2QoGS0RkFR4Y0tLYAkVWi9mb5Ox8fLv9rDr9yPBWDlUmgcjWMVgiIqvQNToA/VoGo7hUh8836GsQ2ZPP1p9EQXEpujUOQP+WIZbeHSKqAQZLRGQ1HhysH11atDMeKdkFsBepOQX4ZvsZdfqRYRxVIrI1DJaIyGpIHo+MMMkIzLzNp2AvPt94EvlFpegS5a+aCBORbWGwRERWQ/J4HjTkLv136xlkXiqCrUvPLVTPRTBXicg2MVgiIqsyrG0Y2oT7IqegGP/dehr20AMur7AEHRv5qT5wRGR7GCwRkVVxdnbCA0P0dZfmbT6NvMJi2KqLuYX4eos+4Ht4KHOViGwVgyUisjpjO0WicZC3msJatENfxNEWSd5VbmEJ2kX64dr24ZbeHSKqJQZLRGR1XF2ccd8g/ejSZxtOqkKOtiYzrwjzN+tHlR4Z1pJ1lYhsGIMlIrJKE3o0QpivBxKlRcgOfTFHWxtVyi4oRtsIqU4eYendIaI6YLBERFbJw9UF04fqV8a9vvwozqTlwlZk5RcZSx88NLSVysMiItvFYImIrNadMU0Q0yxIrSZ7cvF+lJTqYAs+XheH7PxitAprgNEdOapEZOsYLBGR1ZIRmbdu6QIfdxfsPH0R8zZZf6HKbSfTjM2AnxjRhqNKRHaAwRIRWbXoIG88d117dfrNFcdwPCkb1lwq4LHv90GnA27pEYVRHFUisgsMlojI6t3eKxqD24SqVXGP/7APRSXWtzpOp9Phmf8dQEJmPpqH+OCFGzpYepeIyEQYLBGRTbRBeX1CZ/h7ueHg+Sx8uCYW1mbh9rNYcTgJbi5OeP+ObvDxcLX0LhGRiTBYIiKbEO7niZfHd1SnP1wbiwPnMmAtZGrw5V8Pq9PPjGqLjo38Lb1LRGRCDJaIyGbc0KUhxnaOVKviHv9hP/KLSiy9S2ofHvp2LwqKSzGodSim9mtm6V0iIhNjsERENuXlcR0R0sADsck5eHvFMUvvDl77/QiOJWWrfZKVe6ypRGR/GCwRkU0J8nHH6xM6qdNfbDqF7SfTLLYvKw4lYsHWM+r027d2Qaivh8X2hYjMh8ESEdmcYe3CcWvPKLVE/8kf9yOnoLje9yExMx9P/++AOj1tQDM1BUdE9onBEhHZpJnXtUejAC/Ep1/Cq78dqdfHlpypR7/fi4y8InRq5I+nRrat18cnovrFYImIbJKvpxvevKWzOi2Ndv88lFhvjy0VuredTIe3u4sqE+DuyrdSInvG/+FEZLP6tgjBlL5N1en7v9mNj9bFotTM/eN2n0nHnJXH1ekXb+iAZiE+Zn08IrI8BktEZNNmjGmLCd2jIDHSG8uP4d7/7kJmXpHJH0eCsAVbT+POL3aoaTgpY3BzjyiTPw4RWR8GS0Rk0zxcXfDWLZ0x+6ZOajps1ZFkXPfhRhw8n2myx7iQcQmT5u3A8z8fwqWiEvRrGYxXb+yoKosTkf1jsERENk+Cljt6N8aS+/siKlCf9H3Tx1uwaMdZ1bOttuR3f9x9DiPf2YBNsanwdHNWU2//nRqjcqaIyDE46eryTkJKVlYW/P39kZmZCT8/Px4VIguSKThptrv6aLI6L1NlUsjSy92lRveTkl2AGUv+wqojSep898YBePvWrsxRInLAz2+bGVl69dVX0bdvX3h7eyMgIKBavyNx4PPPP4/IyEh4eXlh+PDhOHHiRLnbpKenY+LEieogyf3efffdyMnJMdOzICJz8/d2w+eTeuKpkW3g7AQ1MnTjR5txOjW32vfx+18JGPHOehUoubs4q35vi+/ry0CJyEHZzMjSrFmzVDBz7tw5fPnll8jIuHoTzddffx2zZ8/G119/jWbNmmHmzJn466+/cPjwYXh6eqrbjB49GgkJCfj0009RVFSEu+66C7169cK3335b7X3jyBKRddoSm4qHF+1Fak4hfD1cMfGaJvD3coOPhwu83Fzg4+Gqlv97u+t/Ss7T3LWx+HnfBfX77SP9MOe2LmgbwRFjIntU3c9vmwmWNPPnz8ejjz561WBJnlbDhg3xxBNP4Mknn1SXycEIDw9X93H77bfjyJEjaN++PXbu3ImePXuq2yxfvhxjxoxRQZn8fnUwWCKyXlJpe/q3e7DrzMVq/46LsxMeHNwC04e2Yg0lIjtW3c9vV9ipU6dOITExUU29aeSAxMTEYOvWrSpYkp8yWqUFSkJu7+zsjO3bt+PGG2+s8L4LCgrUVvZgE5F1ivD3xHf3XoNvt59FXEoO8gpLkFdYrP9ZUIK8omL9z8IS5BYWo2mwD14e3xFdo6s33U9E9s9ugyUJlISMJJUl57Xr5GdYWFi5611dXREUFGS8TUVkau/FF180y34Tkem5uThjsqF4JRFRTVk0wfvZZ59VS36r2o4ePQprM2PGDDVkp23x8fGW3iUiIiKyx5ElySeaMmVKlbdp3rx5re47IiJC/UxKSlKr4TRyvmvXrsbbJCfrlxdriouL1Qo57fcr4uHhoTYiIiKyfxYNlkJDQ9VmDrL6TQKe1atXG4MjyS2SXKT7779fne/Tp49KFN+9ezd69OihLluzZg1KS0tVbhMRERGRzdRZOnv2LPbt26d+lpSUqNOyla2J1LZtWyxdulSdlik8WTX3yiuvYNmyZapkwKRJk9QKt/Hjx6vbtGvXDqNGjcK0adOwY8cObN68GdOnT1fJ39VdCUdERET2zWYSvKW4pNRL0nTr1k39XLt2LQYPHqxOHzt2TOUQaZ5++mnk5ubi3nvvVSNI/fv3V6UBtBpLYuHChSpAGjZsmFoFN2HCBLz//vv1+tyIiIjIetlcnSVrxDpLREREtsfu2p0QERERWQKDJSIiIqIqMFgiIiIiqgKDJSIiIqIqMFgiIiIiqgKDJSIiIqIqMFgiIiIiqgKDJSIiIiJ7qOBtzbS6nlLcioiIiGyD9rl9tfrcDJZMIDs7W/2Mjo42xd0RERFRPX+OSyXvyrDdiQmUlpbiwoUL8PX1VQ18TRnxSgAWHx9fZRl24vG2RXx983jbM76+beN4y4iSBEoNGzZU/WErw5ElE5ADHBUVBXORPzyDpfrD412/eLx5vO0ZX9/Wf7yrGlHSMMGbiIiIqAoMloiIiIiqwGDJinl4eGDWrFnqJ/F42xu+vnm87Rlf3/Z1vJngTURERFQFjiwRERERVYHBEhEREVEVGCwRERERVYHBEhEREVEVGCxZsblz56Jp06bw9PRETEwMduzYYeldsgsbNmzA9ddfryq2SsX1n3766YqKrs8//zwiIyPh5eWF4cOH48SJExbbX1s2e/Zs9OrVS1W3DwsLw/jx43Hs2LFyt8nPz8eDDz6I4OBgNGjQABMmTEBSUpLF9tnWffzxx+jcubOxOF+fPn3wxx9/GK/n8Taf//znP+o95dFHH+XxNpMXXnhBHeOyW9u2bc3++mawZKW+//57PP7442op5J49e9ClSxeMHDkSycnJlt41m5ebm6uOpwSjFXnjjTfw/vvv45NPPsH27dvh4+Ojjr38J6SaWb9+vXrj2rZtG1auXImioiKMGDFC/Q00jz32GH755RcsXrxY3V5aB91000081LUk3QTkQ3v37t3YtWsXhg4dinHjxuHQoUM83ma0c+dOfPrppypQLYuvb9Pr0KEDEhISjNumTZvMf7x1ZJV69+6te/DBB43nS0pKdA0bNtTNnj3bovtlb+S/wNKlS43nS0tLdREREbo333zTeFlGRobOw8ND991331loL+1HcnKyOubr1683Hls3Nzfd4sWLjbc5cuSIus3WrVstuKf2JTAwUPfFF1/weJtJdna2rlWrVrqVK1fqBg0apHvkkUfU5Xx9m96sWbN0Xbp0qfA6cx5vjixZocLCQvWtUKZ/yvafk/Nbt2616L7Zu1OnTiExMbHcsZe+QTINymNfd5mZmepnUFCQ+imvcxltKnu8ZUi9cePGPN4mUFJSgkWLFqmRPJmO4/E2Dxk9HTt2bLnXseDxNg9Ji5A0iubNm2PixIk4e/as2Y83G+laodTUVPUmFx4eXu5yOX/06FGL7ZcjkEBJVHTsteuodkpLS1UuR79+/dCxY0fj8XZ3d0dAQACPtwn99ddfKjiSqWPJ21i6dCnat2+Pffv28XibmASjkioh03CX4+vb9OSL6/z589GmTRs1Bffiiy9iwIABOHjwoFmPN4MlIqq3b9/yhlY2v4DMQz5IJDCSkbwff/wRkydPVvkbZFrx8fF45JFHVD6eLMQh8xs9erTxtOSHSfDUpEkT/PDDD2pBjrlwGs4KhYSEwMXF5YoMfjkfERFhsf1yBNrx5bE3renTp+PXX3/F2rVrVQJy2eMt084ZGRnlbs/Xet3It+uWLVuiR48eakWiLGh47733eLxNTKZ9ZNFN9+7d4erqqjYJSmWBiJyWEQ2+vs1LRpFat26N2NhYs76+GSxZ6RudvMmtXr263BSGnJehdTKfZs2aqf9UZY99VlaWWhXHY19zkkMvgZJMA61Zs0Yd37Lkde7m5lbueEtpAclB4PE2HXn/KCgo4PE2sWHDhqkpTxnF07aePXuqPBrtNF/f5pWTk4O4uDhV6sWs7yd1Sg8ns1m0aJFagTV//nzd4cOHdffee68uICBAl5iYyKNugpUre/fuVZv8F5gzZ446febMGXX9f/7zH3Wsf/75Z92BAwd048aN0zVr1kx36dIlHvsauv/++3X+/v66devW6RISEoxbXl6e8Tb33XefrnHjxro1a9bodu3apevTp4/aqHaeffZZtdrw1KlT6vUr552cnHQrVqzg8a4HZVfDCb6+TeuJJ55Q7yfy+t68ebNu+PDhupCQELXS1pzHm8GSFfvggw/UH93d3V2VEti2bZuld8kurF27VgVJl2+TJ082lg+YOXOmLjw8XAWsw4YN0x07dszSu22TKjrOsn311VfG20gQ+sADD6jl7d7e3robb7xRBVRUO1OnTtU1adJEvW+Ehoaq168WKPF413+wxNe3ad122226yMhI9fpu1KiROh8bG2v24+0k/9R9IIyIiIjIPjFniYiIiKgKDJaIiIiIqsBgiYiIiKgKDJaIiIiIqsBgiYiIiKgKDJaIiIiIqsBgiYiIiKgKDJaIiIiIqsBgiYhMrmnTpnj33Xft5sha8vk4OTnhp59+qvbtX3jhBXTt2rXK20yZMgXjx483wd4ROQYGS0RUbfHx8Zg6dSoaNmyoGj43adIEjzzyCNLS0ngUzSQhIQGjR4/m8SWyIAZLRFQtJ0+eVF3UT5w4ge+++w6xsbH45JNPVIdv6eidnp5usSNZUlKC0tJS2JPCwkL1MyIiAh4eHpbeHSKHxmCJiKrlwQcfVKNJK1aswKBBg9C4cWM14rFq1SqcP38e//73v8vdPjs7G3fccQd8fHzQqFEjzJ0713idtKSU6SK5DwkEZKTq4YcfNl5fUFCAJ598Uv2e/H5MTAzWrVtnvH7+/PkICAjAsmXL0L59e3UfX3zxBTw9PZGRkVFuP2Tka+jQocbzmzZtwoABA+Dl5YXo6Gj1uLm5ucbrk5OTcf3116vrmzVrhoULF1Z5XOR4XO1xZeRNjoU8H29vb3Tq1EkFnGUNHjwY06dPx6OPPoqQkBCMHDmywmm4Z555Bq1bt1b307x5c8ycORNFRUVX7Nenn36qnp/c7tZbb0VmZmalz0ECzdmzZ6vnK8+7S5cu+PHHH6t83kQOpc6teInI7qWlpemcnJx0r732WoXXT5s2TXX5Li0tVeel672vr69u9uzZumPHjunef/99nYuLi27FihXq+sWLF+v8/Px0v//+u+7MmTO67du36z777DPj/d1zzz26vn376jZs2KA6ir/55ps6Dw8P3fHjx9X1X331lc7NzU3dZvPmzbqjR4/qcnJydOHh4bovvvjCeD/FxcXlLpP78vHx0b3zzjvqvuR3u3XrppsyZYrxd0aPHq3r0qWLbuvWrbpdu3apx/Dy8lK/U5HLH6Oiy86dO6eew969e3VxcXHG4yHPu2y3+gYNGuieeuop9XxkE/I2vXTpUuPtXn75ZbXfp06d0i1btkw9zuuvv268ftasWeo5Dh06VD3e+vXrdS1bttT94x//MN5m8uTJunHjxhnPv/LKK7q2bdvqli9frvZPjq8c73Xr1lX5uiByFAyWiOiqtm3bdsWHdllz5sxR1yclJRmDpVGjRpW7zW233aYCEfH222/rWrdurSssLLziviR4kkDi/Pnz5S4fNmyYbsaMGeq0fJjL4+3bt6/cbR555BEVJGj+/PNP9aF/8eJFdf7uu+/W3XvvveV+Z+PGjTpnZ2fdpUuXVGAn97tjxw7j9UeOHFGXVRYsVedxKzJ27FjdE088US5YksDtclUddyFBWI8ePcoFS3L8JEDT/PHHH+o5JiQkXBEs5efn67y9vXVbtmwpd79yrO64445KH5fIkbhaemSLiGyH/rO7eiSP6fLz2oqyW265RZ2WaaRRo0ZhzJgxaurL1dUVf/31l8pBkqmmsmRqLjg42HhepgQ7d+5c7jYTJ07ENddcgwsXLqipPZlCGzt2rJqyE/v378eBAwfKTa3Jc5JpqFOnTuH48eNqH3r06GG8vm3btsbfr8zVHleez2uvvYYffvhBTVlKPpI8H5kiK6vs41bm+++/x/vvv4+4uDjk5OSguLgYfn5+5W4j05sy5Vf22MtzPHbsmMqBKktyz/Ly8nDttdeWu1z2sVu3blfdHyJHwGCJiK6qZcuWKnfmyJEjuPHGG6+4Xi4PDAxEaGhotY6m5NLIB7fkO61cuRIPPPAA3nzzTaxfv14FAC4uLti9e7f6WVaDBg2MpyW3RvaprF69eqFFixZYtGgR7r//fixdulTlN2nkvv/1r3+Vy48qG2BIsFQbV3tceW7vvfeeChAlX0nysCQ3SUvi1sjlVdm6dasKzF588UWV0+Tv768e8+2330ZtyTERv/32W7kASzCxnEiPwRIRXZWM6MjIw0cffYTHHntMBSqaxMRENZIyadKkcsHLtm3byt2HnG/Xrp3xvNyHjCbJJsnjMoIjo0oymiEjMZJoLYnYNSXBhOxPVFQUnJ2d1QiPpnv37jh8+LAK/ioi+yAjNRKoSQAkJKi7PHm7po+7efNmjBs3Dnfeeac6L6M8EphJcnpNbNmyRZVrKJtMf+bMmStud/bsWeMol3bsZZ/atGlzxW21BHn5HUncJ6IrcTUcEVXLhx9+qKaOZERjw4YNqubS8uXLVRAlIxKvvvpqudtLgPDGG2+ooEBWwi1evFitEBMy6vLll1/i4MGDqiTBN998o4InCQRk+k0CDwm+lixZoqbHduzYoVZryejH1cjv7tmzR+3PzTffXG50RFaSScAhq8727dunyiD8/PPP6ryQYEKmBWX0afv27Spouueee8oFh7V53FatWqkRNHlsGYWT+09KSqrxK0/uR4IaGU2SaTiZjpNRrMvJ6rzJkyeraceNGzeqkTRZEXf5FJzw9fVVKw8lCP7666/V/crz+OCDD9R5IuJqOCKqgdOnT6vkYFmBJavRoqOjdQ899JAuNTW13O0kwfvFF1/U3XLLLSp5OCIiQvfee+8Zr5eE5ZiYGLUiTlZuXXPNNbpVq1YZr5fE7+eff17XtGlT9TiRkZG6G2+8UXfgwAFjgre/v3+l+9m7d2+VGL1mzZorrpPk7WuvvVatPJPH7ty5s+7VV181Xi9J0JJ8LQnajRs31i1YsEA9n6oSvK/2uLKaUBKq5THDwsJ0zz33nG7SpEnlVqRJgrckil8twVtWywUHB6v7kqR52a+yx0ISvGU130cffaRr2LChztPTU3fzzTfr0tPTK10NJ6sY3333XV2bNm3U8Q4NDdWNHDlSraQjIp3OSQ4Co0YiIiKiinEajoiIiKgKDJaIiIiIqsBgiYiIiKgKDJaIiIiIqsBgiYiIiKgKDJaIiIiIqsBgiYiIiKgKDJaIiIiIqsBgiYiIiKgKDJaIiIiIqsBgiYiIiAiV+39GuD8Eou1DrwAAAABJRU5ErkJggg==",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.plot(data_loadings.T)\n",
    "plt.xlabel('Observed variable')\n",
    "plt.ylabel('Correlation weight')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "700c5222-3c6d-47d6-84d1-d37eb04c818b",
   "metadata": {},
   "source": [
    "Next, we can compute the data array as the latent variable multiplied by the loadings, plus some noise:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "361dadf9-ec0a-4dee-bbfe-a9d65a18116e",
   "metadata": {},
   "outputs": [],
   "source": [
    "data_noise = 0.5*np.random.normal(size=(n_subj, n_var))\n",
    "data = latent_var @ data_loadings + data_noise"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "502182f0-6bd6-478b-bcb0-2acd5d1c3f05",
   "metadata": {},
   "source": [
    "Next, we will simulate the covariates:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "0efe5ea7-56f6-4fcf-a177-bbe9d88440f7",
   "metadata": {},
   "outputs": [
    {
     "ename": "NameError",
     "evalue": "name 'np' is not defined",
     "output_type": "error",
     "traceback": [
      "\u001b[31m---------------------------------------------------------------------------\u001b[39m",
      "\u001b[31mNameError\u001b[39m                                 Traceback (most recent call last)",
      "\u001b[36mCell\u001b[39m\u001b[36m \u001b[39m\u001b[32mIn[1]\u001b[39m\u001b[32m, line 2\u001b[39m\n\u001b[32m      1\u001b[39m n_cov = \u001b[32m4\u001b[39m\n\u001b[32m----> \u001b[39m\u001b[32m2\u001b[39m cov_noise = \u001b[43mnp\u001b[49m.random.normal(size=(n_subj, n_cov))\n\u001b[32m      3\u001b[39m covariates = latent_var + cov_noise \u001b[38;5;66;03m# Effectively all loadings of 1\u001b[39;00m\n\u001b[32m      4\u001b[39m covariates = pd.DataFrame(covariates, columns=[\u001b[33m'\u001b[39m\u001b[33mCov. \u001b[39m\u001b[38;5;132;01m%s\u001b[39;00m\u001b[33m'\u001b[39m % cov \u001b[38;5;28;01mfor\u001b[39;00m cov \u001b[38;5;129;01min\u001b[39;00m \u001b[38;5;28mrange\u001b[39m(\u001b[32m1\u001b[39m, n_cov + \u001b[32m1\u001b[39m)])\n",
      "\u001b[31mNameError\u001b[39m: name 'np' is not defined"
     ]
    }
   ],
   "source": [
    "n_cov = 4\n",
    "cov_noise = np.random.normal(size=(n_subj, n_cov))\n",
    "covariates = latent_var + cov_noise # Effectively all loadings of 1\n",
    "covariates = pd.DataFrame(covariates, columns=['Cov. %s' % cov for cov in range(1, n_cov + 1)])"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "31fa5bd9-4aec-4be6-9000-bf6aaea11727",
   "metadata": {},
   "source": [
    "## Fitting and evaluating the model\n",
    "\n",
    "Since the data are not stratified by any experimental condition, we can fit the model by simply providing the data array and the covariates withut a design matrix specifying conditions:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "id": "4cd6e7d7-bd3e-4e05-9fbf-b5aad9f80699",
   "metadata": {},
   "outputs": [],
   "source": [
    "mod = PLSC(random_state=123) # Seed for reproducibility\n",
    "mod.fit(data=data,\n",
    "        covariates=covariates)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "71aa6949-f5bf-41b1-af27-862abc9d37fb",
   "metadata": {},
   "source": [
    "We can use permutation testing to evaluate the significance of the latent variables identified by the model. This yields one $p$ value per singular value:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "fbf57e51-c870-453d-a3dd-ac6cf8a3d64f",
   "metadata": {},
   "outputs": [
    {
     "ename": "NameError",
     "evalue": "name 'mod' is not defined",
     "output_type": "error",
     "traceback": [
      "\u001b[31m---------------------------------------------------------------------------\u001b[39m",
      "\u001b[31mNameError\u001b[39m                                 Traceback (most recent call last)",
      "\u001b[36mCell\u001b[39m\u001b[36m \u001b[39m\u001b[32mIn[2]\u001b[39m\u001b[32m, line 1\u001b[39m\n\u001b[32m----> \u001b[39m\u001b[32m1\u001b[39m \u001b[43mmod\u001b[49m.permute(\u001b[32m1000\u001b[39m, return_null_dist=\u001b[38;5;28;01mFalse\u001b[39;00m)\n\u001b[32m      2\u001b[39m \u001b[38;5;28mprint\u001b[39m(mod.pvals_)\n\u001b[32m      3\u001b[39m is_sig = mod.pvals_ < \u001b[32m0.05\u001b[39m\n",
      "\u001b[31mNameError\u001b[39m: name 'mod' is not defined"
     ]
    }
   ],
   "source": [
    "mod.permute(1000, return_null_dist=False)\n",
    "print(mod.pvals_)\n",
    "is_sig = mod.pvals_ < 0.05\n",
    "sig_lvs = np.where(is_sig)[0]\n",
    "print(sig_lvs)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8b632b84-70eb-49c5-b0e2-b897964b7e54",
   "metadata": {},
   "source": [
    "Thus there is one significant pair of latent variables, in line with how we simulated the data. Next, we can perform bootstrap resampling to assess the reliability of the data saliences and to evaluate how reliably the design scores correlate with each covariate."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "id": "304d945f-dd17-4c35-9999-396c0ecc4127",
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Getting resamples: 100%|█████████████████████████████████████████████████████████| 1000/1000 [00:00<00:00, 3664.26it/s]\n",
      "Resampling: 100%|█████████████████████████████████████████████████████████████████| 1000/1000 [00:01<00:00, 961.61it/s]\n"
     ]
    }
   ],
   "source": [
    "mod.bootstrap(1000, return_boot_stat_dist=False)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "896ff78b-3185-4c0b-87e2-5829300d9b7c",
   "metadata": {},
   "source": [
    "## Visualizing the model\n",
    "\n",
    "Next we will visualize the bootstrap ratios ($z$ scores) of the data saliences and the correlation of the data scores with the covariates:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "1f6b7fd7-ae86-4ef3-8fa1-acb12da80224",
   "metadata": {},
   "outputs": [
    {
     "ename": "AttributeError",
     "evalue": "'PLSC' object has no attribute 'design_sals_z_'",
     "output_type": "error",
     "traceback": [
      "\u001b[31m---------------------------------------------------------------------------\u001b[39m",
      "\u001b[31mAttributeError\u001b[39m                            Traceback (most recent call last)",
      "\u001b[36mCell\u001b[39m\u001b[36m \u001b[39m\u001b[32mIn[10]\u001b[39m\u001b[32m, line 20\u001b[39m\n\u001b[32m     16\u001b[39m ax[\u001b[32m0\u001b[39m].bar(x=labels,\n\u001b[32m     17\u001b[39m           height=mod.boot_stat_val_[:, lv_idx],\n\u001b[32m     18\u001b[39m           yerr=mod.get_boot_stat_yerr(lv_idx))\n\u001b[32m     19\u001b[39m ax[\u001b[32m0\u001b[39m].set_ylabel(\u001b[33m'\u001b[39m\u001b[33mCorrelation with data score\u001b[39m\u001b[33m'\u001b[39m)\n\u001b[32m---> \u001b[39m\u001b[32m20\u001b[39m ax[\u001b[32m1\u001b[39m].plot(\u001b[43mmod\u001b[49m\u001b[43m.\u001b[49m\u001b[43mdesign_sals_z_\u001b[49m[:, lv_idx])\n\u001b[32m     21\u001b[39m ax[\u001b[32m1\u001b[39m].set_ylabel(\u001b[33m'\u001b[39m\u001b[33mBootstrap ratio (z score)\u001b[39m\u001b[33m'\u001b[39m)\n\u001b[32m     22\u001b[39m ax[\u001b[32m1\u001b[39m].set_xlabel(\u001b[33m'\u001b[39m\u001b[33mObserved variable\u001b[39m\u001b[33m'\u001b[39m)\n",
      "\u001b[31mAttributeError\u001b[39m: 'PLSC' object has no attribute 'design_sals_z_'"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAosAAAHrCAYAAACn9tfQAAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjgsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvwVt1zgAAAAlwSFlzAAAPYQAAD2EBqD+naQAATodJREFUeJzt3Qm8TfX+//GPecqYWYYMmUPEpYFQRGi4hVxTohSKJjIVGZIkUqIk90+UJBcpEQ1EoUIoQ3HNkjnz/j/e3/tb57HPOXudQWcfZ5/zej4ei7PXXmvv71pr77U/6/v5fr8rXSAQCBgAAAAQQvpQMwEAAACCRQAAAMSJmkUAAAD4IlgEAACAL4JFAAAA+CJYBAAAgC+CRQAAAPgiWAQAAIAvgkUAAAD4IlhEqtWsWTPr2rXr5S4G/k+6dOnsueeeY39cogYNGrgp0vb93yk3ksaiRYvsiiuusIMHD7JLcUkIFpFoU6dOdT8+33//fbT5c+bMcfPfeust33UXL17slhk3blyc73HmzBl75plnrGjRopYtWzarU6eOWzehvvnmG/vss8/cawTbunWr/fOf/7S8efNa9uzZ7cYbb7Qvvvgi5Gu89tprVrFiRcuSJYsVK1bM+vTpYydPnkzQ+/fu3duuu+46y5cvn3sfvY5+rE+cOBFtuY0bN9q9995rpUuXdsvlz5/fbr75ZvvPf/4T6zXnzp1rFSpUsNy5c1uLFi1sz549sZZp2bKldevWLUFlBJA4u3fvtvvuu8/y5MljuXLlslatWtn27dsTvP6KFSvcOUff9cKFC1uvXr1inRMSc/7TOa5Lly5WpUoVy5Ahg5UqVSrk+zZt2tTKli1rI0aMSOQWA/9H94YGEuOdd97R/cQD3333XbT5p0+fDuTOnTtwyy23+K7bqVOnQIYMGQL79++P8z3atGkTyJgxY+DJJ58MvPnmm4G6deu6x1999VWCytiqVavAbbfdFm3ezp07A/nz5w8UKlQoMGzYsMDYsWMD1apVc6+7fPnyaMs+/fTTbhv/+c9/Bt54441Az5493XIxX9PPDTfcEOjVq1dg3LhxgUmTJgW6d+8eyJIli5t/4cKFqOUWLFgQaNKkSeC5555zy6lMN910k3tvbbdn27ZtgcyZMwfat28feP311wPXXHNNrLIsWrTI7f8DBw4EUqK//vorcO7cuctdjIhVv359N0Xavv875U5Jjh8/HihXrlygYMGCgRdffDEwZsyYQPHixQNXXXVV4NChQ/Guv27dukDWrFkDNWrUcOeU/v37u3NC06ZNL/n817FjR/ea9erVc+UoWbKk7/vrvJE9e/bAsWPHLnEPIC0jWESSBYvSpUuXQPr06QO7d+8O+YOlYCbUyTHYqlWr3Ou/9NJL0dYtU6aMO2nGR4GoTqxvvfVWtPmPPPKIm7958+aoeSdPnnQn/Ouuuy5q3p49e9xyCsyCjR8/3pVr3rx5gUsxevRot/7KlSvjXO78+fMuiC1fvnzUPP24lC5dOnDx4kX3+IsvvgikS5fO7RdRIFCxYsXAyy+/HEhJFBh7ZUTaDLoitdwxKUDU93f16tVR8zZt2uQufvv16xfv+rfffnugSJEigaNHj0bNmzx5snvNTz/99JLOfzrPnj171v3dvHnzOINFnRdV1rfffjsRWw38D2loJKl//etfdvHiRZs5c2as5xYsWGBHjx61du3axfkas2fPdimV4HRq1qxZXbpl5cqVtmvXrjjX1/ucP3/eGjduHG3+V199ZTVq1LDy5ctHzVM6SKnbtWvX2q+//urm6T20fps2baKt7z0OtW0J4aWIjhw5Eudy2vbixYtHW+6vv/5yqS+l8EXpbV3sab6XMr9w4YL17NkzweVR6uqWW26JNV/HT2l3pes9o0ePtnr16tmVV17p0mI1a9Z0xykmla9Hjx42ffp0q1y5skvhq71UqHZzv//+uz3yyCPueOg19dpKyf/2228hmz2oaYGaAhQoUMBy5Mhhd911V8g2WJ988onVr1/fcubM6VKF119/vc2YMSPaMqtWrXKpOaX09RnQ8nr9hFCKcPDgwS6tp+3TsXr66afdfE/Hjh3dZ3bTpk3R1m3SpIlrAuE1IfC27csvv7SHHnrI7QOVuUOHDvbnn3/GWY6zZ8/aoEGD3LHQdmif3HTTTSGbVcTc9/pb89Qso1OnTu6zpdfo3LmznTp1Ktb6/+///T/3PjpO+uzpuxDqezhp0iQrU6aMW6527druO/d3qGxqa6dUr/adtlGp2SFDhrjPf3LS512fJU0eNQtp1KiRvf/++3Gue+zYMZdG1vlRx9ej46ztC14/Mec/7YtMmTIlqPwFCxa0a6+91j7++OMEbzPgIVhEklJ7u6uuuirWj7Nonn6Y77zzzjhfY926dXbNNddEO6mKfnzkhx9+iLddkH50S5YsGW2+fsz1IxaTyiRr1qyJWk5iLhtzufgo4Dx06JALDNS2aMCAAS6A8bYjmNpCatlt27bZK6+84gIe/Qh59AOl/fLee+/Zjh07bNiwYS5YUeChgOn555+3MWPGJPiHQ1q3bu2ClH379kWb//XXX7syBwfLr776qgu09SM9fPhwy5gxowvsFJjHtHTpUtdmU6+v9fzaUX333XfuWOl91Ib14YcftiVLlrjOEKECFgXCP/74owvUunfv7tp1KjANpuCrefPmdvjwYevXr5+NHDnSqlevHhWweuXT51Q/4HotbY8C84YNG9rq1avj3GcKpHVxoeBZ7UbHjx/vPs86Ztre4P2loFZBo4J4efPNN93nQOvoRz6YtkOBpYI4BRAKtvW6cQVEKr/aB2t/vfjii25dfRYUVMX3HfGo/d3x48ddWzb9rf2nz1IwfdZUpnLlyrnP2OOPP+6Ok/Zh8AXN22+/7QJetcUbNWqU3XDDDW5fxXdxFx/tPwX2hQoVcq+roFXHTVN81B5Q36v4Jl3Exnfcf/rpJ6tVq1as5/R91vdW+9HP+vXr3fkg5vqZM2d2n099t5Pq/BcX7Tt954BE+78aRiBJ0tDy1FNPuee3bNkSNU+pF7Wtadu2bbyvX7ly5UDDhg1jzd+4caN73YkTJ8a5/o033hioWbNmrPktWrQI5MmTJ1abHaV29LpKE8uaNWvc46FDh8ZqE6j5V1xxRSAhlG7W8t6ktLLSx6E89NBDUcspja+2kocPH462jNpAesvky5cvsHTpUje/a9eu8ab2Q9Hx0WspvR4zXa9tPHXqVNS84L9Fqa8qVarEOk5e+XWsYtJzgwcP9n3N4H02bdq0WJ+3xo0bR6XhpXfv3i6tduTIEfdY/+fMmTNQp06dWKlvbz39r3Znaica/Foqy9VXXx249dZb49xn//73v932xWw7ps+kyvjNN99EzVNqUfNeeOGFwPbt290+vfPOO6Ot522bPq9eOlFGjRrl5n/88ce+6Vw1Vzhz5ky01/vzzz9dm9wHHnggzn2vvzUv5nJ33XVX4Morr4x6/Ntvv7l9rDa+wdavX++aanjzVXa15atevXq0Mqkdrt7nUtPQapOn9dVm2KPjppSr2vAePHgwQevHN8VXPr2PlhsyZEis5yZMmOCeC27eEtMHH3zglvnyyy9jPXfvvfcGChcu/LfPf/GloWX48OHuNeJrMw7ERM0ikpxSLRJcu/jhhx/a6dOn401Bi1KrSu/FpFSM93xc/vjjD1fjFpNqo1QTohogXb3/8ssvrpbE69Xtva56Mav3oWpr3nnnHZcWVU2fak1Ucxff+3sqVarkUk/qxaw0pVJooXo+isqhZd999127/fbbXW2K0ozBVFul1K1SqPpfKWTVMkybNs3VbKl2RPteKWTVNsVMgcak2gvVasyaNStqnt5XaTDVmgXXrAb/rfSo3kspT6XvY1JKV9sen+DXPHfunDtuqi1VSjTU6yot56XhRe+v8mpfiPafanf69u0b9VnxeOtpf6m5wf333+/ez6tZUs2uanJV06paJD8ffPCB69mu9GNwzZRqJSU4BXzbbbe5z4xqY++++25XJtUuhqJtC64V1mdVtbcLFy70LYtSlaqZEpVZtale7VWo/ReKanODaZ9qv6jW0hvhQK+tWsfg7VXtoWoave3Vd+jAgQPu9bwyeWlkpbf/ruAaZK+pg74fn3/+eZzr6Xunz0V808svvxzn63jf+Us9L8W3fvC6f/f8FxfvvKhjCCRGxkQtDSSA2sWoPZxSpl47KQWOGhZGKbKEBBHB7b88Cja95+MTKn2nIEwpQAUTCghFwYnSbPpRUduh4OBWQeUDDzwQ9cOs9nLLly+3LVu2WEIojeS1m9QQG9oH+l8/5NWqVYu2rIIPTaKUnwINBWwKDIMDpBIlSrjJo6E39AOtdRUoKuWnNkkKOrX+5s2bXdDhR9v47LPPuiFBFGQuW7bM/egHp1Rl/vz59sILL7hgK/jYBJfNc/XVVydo/+hHT+lPBeR6/+BjFiotGLzdwT98Xts+pQJFnz0/XrtUpYf96L1DXWx46ysIV4o5FO27YEpX63hov+n4q91YKAq8gumzWKRIkVjtN2PScVago+OsgDuxxyCufarPr7ZXxyVm+TxegOsF7DGX0/MaFurvSJ8+fazX0IWOxLd/dNGSkAuX+HjnnEs9L8W3fswLs797/vPjfcdCfW+BuBAsIiwUuCgoU42D2jCqBkK1LHEFLh79SCp4iGnv3r3u/5jtvWJSe0W/zgGqkVAjfrU/8toLqa1V8A+QKHBS2z39WKpNn34EVZui9w5eLjFUu9S+fXvXQSZmsBiTOpdof6n2M7hDTjDVCCpwmTdvnqthUyN5tYlTzZI6l0yePNm+/fZbN66bHwWFatunGjPVbuo1VBOkNmIedVJQ2zO1UXv99dfd8VEQoCAvVNvUhP6YqQ2iXkPvW7duXfe++hFTG8ZQtXsK2ENJTEcH73Vfeukld+xDCb5oCLV+1apVXdu9UNTZJZhqsL0AUu3W2rZta0lFnU5Uc6e2jU899ZQLRLWPFIB7gXN84tun2l4dE9Wsh1o2rn2VEijwT0hNnM4F6rjjR8+pts87ByX2vKTvTPCyMdcPXvfvnv/i4p0XdeEOJAbBIsJCP4oKQhRMqKOJgpmEpKBFP+IKLpUKC27krVo27/m4qJZNNYN+lA5WcOJRKksBjhrkx6Qg0ast+fnnn90JWz/Ql0K1Bfrxja8xvXg/cH7LqgOIAoShQ4e6tO3+/ftdzZL3Q6LtUS1RqB+dYKqBUsN5BZ4KpJV2VPARnAbTvlQK7NNPP402X4He36F0t2r4glOAqj2Jr7e4H/XClQ0bNrga47iWCa71Tex7qJONUtbx1c4ota0LE9VsqSe5OmeoB3dwb1qPLkqCe6aruYI+a7oLUVz7TzVu3mD4noR0/EjM9ipw1OckroskrzOZtsNLyYs+k+qQFd/FUVz0nVFv6OD310WU+HWe8jz22GOu9jU+ajqhWvW4ajd1kRDzRgTeeUnHQZ3X/Ki2WxfKWl8pfY9S6ap1Dp73d89/cdGxUKDoVzMO+KHNIsJC6S21f1IQohoQ/djoBzMhVKum4FLDcAQHWgpO1JYwZu1NTAoEdQWdkDsrqGegfmw1LEVcbav0g6VUtXpEB7fz0o+hUoDBNQYKdoJTgh7vzjbBPSJjpi2911Q7RAV8fik0tadUMOjdzlC1qfoxUlm8NknqGava0PiodlE1kFOmTHHrxUxBq0ZJwYjXq9dL/6kt5t+h141ZK6hmAsHvkxhK3esHWzVrXsrO472PeoMqAFJ6OFT70fhuh6YfdQXgqrUNFeAH3+FHd+DYuXOnC1ZUE6nARsFxqBSjPuvBn5k33njDtT9U0wk/Xk1f8D5UQKHhVZKKasP1PuohHfNY6bHaN3qfaQUgEydOjNbWVr2rLzX4D6ahoYLfV49Vux08YkA42yx65yX14A8OGNUkRb3rNTJAMH0Pdew9Orfo4kTnwuBe0//+97/d5zB4/b97/ouLRnIIvlAGEoqaRVwyBRfBQ5IEX83rR1upaDXc1zAs/fv3T/Dr6oSok6dqJhVMqZZIP7gKULyUcVw0dIoCJ9UYBo9VpnZV+rFXSlVBlG61px83tbHU8Ckxt0EBh67i9SOuGlINq6JyBLfzUuCgDg8KAvTDKKqhUFtCnfRVK6kfT6VyFZTqR9XrACRKNasGQSlepb6V8tawKfqx0Q9YqDSffoSURtWwNV7AoO1Ve0ildPX8Rx995GoZE/LDoH3y5JNPuknptpg1btqfCnaUmlbHEB2TCRMmuOOidP6luuOOO9yPpX5IFRQryNExU+B7KVQLo44+Dz74oKu9U1kVUKsmUDWxOnaqIVLQriBMqXrV/Gm/6ziqNkevEepWix41I1CqXhcMWl610fph1/HSfNW+6hgrgFDKXrV8XvtY/dir49HAgQNdLWMwfUYU+OhYKADRumo+oM9qXPtPnynVVuoYqdZIn2ftS7+OVImlwFptVfVd1PdPtc76buu99BnT90ufGwVuWk6fZ9Us6oJDy2ibQ7VZ1H5Q+9+ENCFQrbbOM/qO6dyglLg++2prG18NWVK1WRSNCaqLBO1rb5v1vdCQPk888US0ZXVOiFlbqbbRumDWfO23//73v+47rouc4GYfiTn/6funZiiiMTOVidBxENXmqt2yR6+l5R999NEk2R9IY2L1jwbi4Q334Tft2rXLLaehX3Q7K837+eefE7VfNfSJbnWlISX0Gtdff70buiahWrZsGWjUqFG0eSqPbgOo19SwGxoq5Zlnngl5+ytto+6ikiNHDjcci17LG6om2I4dO9z2aYgOz9atWwMdOnRwd1zJli2bGzJIw2FouJITJ05EW/+9995zQ8JouBMNRZI3b173OHjIlFBDbdx9992x5ms4DA0PpPLqjjTff/99gveXbkOo7XjwwQdDPq+7PmjIGR2LChUquP3jDb8STI8fffTRkK8Rc/gWDfPSuXNndwtGDSuj4Ww0/IiG/wjen35DNWkYIs2PORyR7rCj259p3+fKlStQu3Ztt59j3npN+1DDxGib9J733XdfYMmSJfHuKw0To7t56JhqXR0zDX3z/PPPuyGi9HnS6+kYxLzFnob70dA73l18vG3T7Sa7devmXkv7ol27doE//vgj2roxh87REDIaCkXvpXLoNnLz5893+y7mECp+Q+fEHHrGK48+18E+/PBDNySVvg+a9BnQcQ4eHsu7pZy+VypPrVq13FAxoe7gov0VPFyMH22L3k+3u9TtLXW7On1XVP7g22YmF53bNKyVPlc6TnfccUfg119/jbWc33A8GnJJn02dEwoUKOD2YajzT0LPf3Gdi4O/Q95doLjdHy5VOv1zuQNWIKmpJk+1F6rx8evJCVxuqo1W7abSm6EGfE6NlIZVDfbYsWPjreVS+2C1zUyqmtK0TIPq65yo2ncgsWiziFRJ7SWV3omZ7gNweWksS6X+vfa2CD+l8dX5SKlt4FLQZhGplto2AUhZ1OZPE5KP2kRSO4u/g5pFAAAA+KLNIgAAAHxRswgAAABfBIsAAADwRbAIAAAAXwSLAAAA8EWwCAAAAF8EiwAAAPBFsAgAAABfBIsAAABIPcHihAkTrFSpUpY1a1arU6eOrV69Os7lP/jgA6tQoYJbvmrVqrZw4cJkKyuAtHv/4xYtWljRokUtXbp0Nnfu3HjXWbZsmV133XWWJUsWK1u2rE2dOjVZygoAqSpYnDVrlvXp08cGDx5sa9eutWrVqlmTJk3swIEDIZdfsWKFtW3b1rp06WLr1q2zO++8000bNmxI9rIDSDtOnjzpzk+6uE2IHTt2uPsl33LLLfbDDz/Y448/bg8++KB9+umnYS8rAKSq2/2pJvH666+31157zT2+ePGiFS9e3Hr27Gl9+/aNtXzr1q3dSXv+/PlR8/7xj39Y9erVbeLEicladgBpk2oWP/roI3eh6ueZZ56xBQsWRLuQbdOmjR05csQWLVqUTCUFgNAyWoQ4e/asrVmzxvr16xc1L3369Na4cWNbuXJlyHU0XzWRwVQTGVdK6MyZM27yKCA9fPiwXXnlle6kDyBt0nX18ePHXWpZ556kpHOVzmUxz1WqYfTDuQpAcp2rIiZYPHTokF24cMEKFSoUbb4eb968OeQ6+/btC7m85vsZMWKEPf/880lUagCpza5du+yqq65K0tf0O1cdO3bM/vrrL8uWLVusdThXAUiuc1XEBIvJRTWXwbWRR48etRIlSridnitXrgS9RpXBkd/OaMPzTRK9Dtuddo53ajjWid1uBW5q9pIzZ05LLecqAKnPsTCcqyImWMyfP79lyJDB9u/fH22+HhcuXDjkOpqfmOVFPRE1xaSTb0JPwOmzZLdIdyk/Nmx32jneqeFYX+rnPBzNUfzOVSpfqFrFpDpXAUi90iXhuSpiekNnzpzZatasaUuWLInWnlCP69atG3IdzQ9eXhYvXuy7PABcDpyrAKRkERMsilIukydPtnfffdc2bdpk3bt3d72dO3fu7J7v0KFDtA4wjz32mOtJ+PLLL7t2jc8995x9//331qNHj8u4FQBSuxMnTrghcDR5Q+Po7507d7rHOk/pfOV5+OGHbfv27fb000+7c9Xrr79u77//vvXu3fuybQMARFwa2hsK5+DBgzZo0CDXIFxD4CgY9BqG60Qc3POnXr16NmPGDBswYIA9++yzVq5cOdcTukqVKpdxKwCkdroo1ZiJHq9tYceOHd1g23v37o0KHOXqq692Q+coOHz11Vddo/S33nrL9YgGgMstooJFUa2gX82g7oAQ07333usmAEguDRo0cMNX+Al1dxato5sHAEBKE1FpaAAAACQvgkUAAAD4IlgEAACAL4JFAAAA+CJYBAAAgC+CRQAAAPgiWAQAAIAvgkUAAAD4IlgEAACAL4JFAAAA+CJYBAAAgC+CRQAAAPgiWAQAAIAvgkUAAAD4IlgEAACAL4JFAAAA+CJYBAAAgC+CRQAAAPgiWAQAAIAvgkUAAAD4IlgEAACAL4JFAAAA+CJYBAAAgC+CRQAAAPgiWAQAAIAvgkUAAAD4IlgEAACAL4JFAAAA+CJYBAAAgC+CRQAAAPgiWAQAAIAvgkUAAAD4IlgEAACAL4JFAAAA+CJYBAAAgC+CRQAAAPgiWAQAAIAvgkUAAAD4IlgEAACAL4JFAAAAECwCAAAg8ahZBAAAgC+CRQAAAPgiWAQAAIAvgkUAAAD4IlgEAACAL4JFAAAA+CJYBAAAgC+CRQAAAPgiWAQAAIAvgkUAAAD4IlgEAACAL4JFAAAA+CJYBAAAgC+CRQAAAPgiWAQAAIAvgkUAAAD4IlgEAACAL4JFAAAA+CJYBAAAgC+CRQAAAPgiWAQAAIAvgkUAAAD4IlgEAACAL4JFAAAA+CJYBAAAgC+CRQAAAPgiWAQAAIAvgkUAAAD4IlgEAACAL4JFAAAA+CJYBAAAQOQHi4cPH7Z27dpZrly5LE+ePNalSxc7ceJEnMv37NnTypcvb9myZbMSJUpYr1697OjRo8labgAAgEgWMcGiAsWNGzfa4sWLbf78+fbll19at27dfJffs2ePm0aPHm0bNmywqVOn2qJFi1yQCQDhNmHCBCtVqpRlzZrV6tSpY6tXr45z+bFjx0Zd3BYvXtx69+5tp0+f5kABuOwyWgTYtGmTC/S+++47q1Wrlps3fvx4a9asmQsGixYtGmudKlWq2Icffhj1uEyZMjZs2DD717/+ZefPn7eMGSNi0wFEoFmzZlmfPn1s4sSJLlBUINikSRPbsmWLFSxYMNbyM2bMsL59+9qUKVOsXr169ssvv1inTp0sXbp0NmbMmMuyDQAQUTWLK1eudKlnL1CUxo0bW/r06W3VqlUJfh2loJXGjitQPHPmjB07dizaBACJoQCva9eu1rlzZ6tUqZILGrNnz+6CwVBWrFhhN9xwg91///2uNvK2226ztm3bxlkbybkKQHKJiGBx3759sa7GFfDly5fPPZcQhw4dsqFDh8aZupYRI0ZY7ty5oyalgwAgoc6ePWtr1qxxF7QeXdjqsS58Q1FtotbxgsPt27fbwoULXfaEcxWANB0sKu2iNEtc0+bNm//2+6h2sHnz5u4K/7nnnotz2X79+rkaSG/atWvX335/AGmHLkwvXLhghQoVijZfj/0ublWjOGTIELvxxhstU6ZMrtlMgwYN7Nlnn/V9H85VAJLLZW2498QTT7h2OXEpXbq0FS5c2A4cOBBtvtodqseznovL8ePHrWnTppYzZ0776KOP3Ik4LlmyZHETACSXZcuW2fDhw+311193bRy3bt1qjz32mMuGDBw4kHMVgLQbLBYoUMBN8albt64dOXLEpWlq1qzp5i1dutQuXrzoTqxx1SiqUbmCv3nz5rleiQAQTvnz57cMGTLY/v37o83XY7+LWwWE7du3twcffNA9rlq1qp08edI1m+nfv79LYwPA5RIRZ6CKFSu62kE1GFebnm+++cZ69Ohhbdq0ieoJvXv3bqtQoUJUmx8FimokrhPu22+/7R4rBaRJKSIACIfMmTO7i9olS5ZEzdOFrR7rwjeUU6dOxQoIFXBKIBDgQAG4rCJm/Jjp06e7ALFRo0bupHrPPffYuHHjop4/d+6cG5ZCJ11Zu3ZtVE/psmXLRnutHTt2uB6HABAOGjanY8eObgSH2rVru6FzdOGq3tHSoUMHK1asmOtQJy1atHA9qGvUqBGVhlZto+Z7QSMAXC4REyyq57PGIvOj4C/4ClyNw7kiB3A5tG7d2g4ePGiDBg1y2Yzq1au7sWK9Ti87d+6MVpM4YMAA16FP/ytLouY5ChQ1NiwAXG4REywCQCRRJkSTX4eWmEOBDR482E0AkNJERJtFAAAAXB4EiwAAAPBFsAgAAABfBIsAAADwRbAIAAAAXwSLAAAA8EWwCAAAAF8EiwAAAPBFsAgAAABfBIsAAADwRbAIAAAAXwSLAAAA8EWwCAAAAF8EiwAAAPBFsAgAAABfBIsAAADwRbAIAAAAXwSLAAAA8EWwCAAAAF8EiwAAAPBFsAgAAABfBIsAAADwRbAIAAAAXwSLAAAA8EWwCAAAgKQNFrdt22YDBgywtm3b2oEDB9y8Tz75xDZu3HgpLwcAAIDUEiwuX77cqlataqtWrbI5c+bYiRMn3Pwff/zRBg8eHI4yAgAAIFKCxb59+9oLL7xgixcvtsyZM0fNb9iwoX377bdJXT4AAABEUrC4fv16u+uuu2LNL1iwoB06dCipygUAAIBIDBbz5Mlje/fujTV/3bp1VqxYsaQqFwAAACIxWGzTpo0988wztm/fPkuXLp1dvHjRvvnmG3vyySetQ4cO4SklAAAAIiNYHD58uFWoUMGKFy/uOrdUqlTJbr75ZqtXr57rIQ0AAIDUI2NiFg4EAq5Gcdy4cTZo0CDXflEBY40aNaxcuXLhKyUAAAAiI1gsW7asG09RwaFqFwEAAJB6JSoNnT59ehck/vHHH+ErEQAAACK3zeLIkSPtqaeesg0bNoSnRAAAAIjMNLSox/OpU6esWrVqblDubNmyRXv+8OHDSVk+AAAARFKwOHbs2PCUBAAAAJEfLHbs2DE8JQEAAEDkB4ty4cIFmzt3rm3atMk9rly5srVs2dIyZMiQ1OUDAABAJAWLW7dutWbNmtnu3butfPnybt6IESPcMDoLFiywMmXKhKOcAAAAiITe0L169XIB4a5du2zt2rVu2rlzp1199dXuOQAAAKThmsXly5fbt99+a/ny5Yuad+WVV7ohdW644YakLh8AAAAiqWYxS5Ysdvz48Vjzdds/DaUDAACANBws3nHHHdatWzdbtWqVu/2fJtU0Pvzww66TCwAAANJwsDhu3DjXZrFu3bqWNWtWNyn9rHtGv/rqq+EpJQAAACKjzWKePHns448/dr2ivaFzKlas6IJFAAAApC6XNM6iKDgkQAQAAEjdEp2Gvueee+zFF1+MNX/UqFF27733JlW5AAAAEInB4pdffukG5Y7p9ttvd88BAAAgDQeLfkPkZMqUyY4dO5ZU5QIAAEAkBotVq1a1WbNmxZo/c+ZMq1SpUlKVCwAAAJHYwWXgwIF2991327Zt26xhw4Zu3pIlS+y9996zDz74IBxlBAAAQKQEiy1atLC5c+fa8OHDbfbs2ZYtWza79tpr7fPPP7f69euHp5QAAACInKFzmjdv7iYAAACkbolus7hr1y7773//G/V49erV9vjjj9ukSZOSumwAAACItGDx/vvvty+++ML9vW/fPmvcuLELGPv3729DhgwJRxkBAAAQKcHihg0brHbt2u7v999/3/WOXrFihU2fPt2mTp0ajjICAAAgUoLFc+fOWZYsWdzf6tTSsmVL93eFChVs7969SV9CAAAARE6wWLlyZZs4caJ99dVXtnjxYmvatKmbv2fPHrvyyivDUUYAAABESrCo+0K/+eab1qBBA2vbtq1Vq1bNzZ83b15UehoAAABpdOgcBYmHDh1yt/bLmzdv1Pxu3bpZ9uzZk7p8AAAAiLRxFjNkyBAtUJRSpUolVZkAAAAQqWloAAAApB0EiwAAAPBFsAgAAABfBIsAAABI2g4uJ0+etOXLl9vOnTvt7Nmz0Z7r1avXpbwkAAAAUkOwuG7dOmvWrJmdOnXKBY358uVzQ+lo2JyCBQsSLAIAAKTlNHTv3r2tRYsW9ueff1q2bNns22+/td9//91q1qxpo0ePDk8pAQAAEBnB4g8//GBPPPGEpU+f3o23eObMGStevLiNGjXKnn322fCUEgAAAJERLGbKlMkFiqK0s9otSu7cuW3Xrl0WLocPH7Z27dpZrly5LE+ePNalSxc7ceJEgtYNBAJ2++23W7p06Wzu3LlhKyMAAICl9WCxRo0a9t1337m/69evb4MGDbLp06fb448/blWqVLFwUaC4ceNGW7x4sc2fP9++/PJLd4vBhBg7dqwLFAEguUyYMMHd2Spr1qxWp04dW716dZzLHzlyxB599FErUqSIZcmSxa655hpbuHBhspUXAJIsWBw+fLg7mcmwYcPcbf+6d+9uBw8etDfffNPCYdOmTbZo0SJ766233En3xhtvtPHjx9vMmTNtz5498abNX375ZZsyZUpYygYAMc2aNcv69OljgwcPtrVr11q1atWsSZMmduDAgZA7S6NK3Hrrrfbbb7/Z7NmzbcuWLTZ58mQrVqwYOxdA5PWGrlWrVtTfSkMriAu3lStXutRz8Hs3btzYpcNXrVpld911V8j11GP7/vvvd1f4hQsXTtB7qQ2mJs+xY8eSYAsApCVjxoyxrl27WufOnd3jiRMn2oIFC9xFa9++fWMtr/lqarNixQrX1EdUKwkAEVmz2LBhQ5cuiUlBlZ4Lh3379rnANFjGjBndsD16Lq6e2/Xq1bNWrVol+L1GjBjh2l96kzrvAEBCqZZwzZo17oLWowtbPdaFbyjz5s2zunXrujR0oUKFXJMeZXEuXLjg+z66qNV5N3gCgBQRLC5btizWQNxy+vRp++qrrxL1WrrCVlvCuKbNmzfbpdDJd+nSpa69YmL069fPjh49GjWFs9MOgNRH484qyFPQF0yP/S5ut2/f7tLPWk/tFAcOHOiaz7zwwgu+78OFLYAUl4b+6aefov7++eefo530dIJTOjqx7Ws0BE+nTp3iXKZ06dIuhRyzrc/58+dd2sYvvaxAcdu2bS59Heyee+6xm266yQW9oahhuSYASC4XL1502ZNJkya5Ick0bu3u3bvtpZdecu0e/S5s1S7So5pFMiEALmuwWL169ajavlDpZg3QrU4niVGgQAE3xUfpGaW+ldrRSdQLBnWCVYcXv1rLBx98MNq8qlWr2iuvvOIGFQeAcMifP78L+Pbv3x9tvh77Xdyq06DaKmo9T8WKFd1FuTI5mTNnjrUOF7YAUlywuGPHDjdeoWr6NAREcJCnE5muioNPdElJJ82mTZu6BuNqKH7u3Dnr0aOHtWnTxooWLeqW0VV4o0aNbNq0aVa7dm13Ug51Yi5RooRdffXVYSknAOh8qIvaJUuW2J133ul2iC5s9VjnrVBuuOEGmzFjhlvOG8f2l19+cUFkqEARAFJksFiyZEn3v05ml4PGctSJVgGhTqZKJ48bNy7qeQWQGm5CPaAB4HJSerhjx45uBAddvKrt9MmTJ6N6R3fo0ME121G7Q9HwY6+99po99thj1rNnT/v1119dB5devXpxIAFE3tA5we0WdfeWmJ1dWrZsaeGgns+68vajYSZU8xmX+J4HgKTQunVrN/asblqgVLKa8ahdt9fpRedOrwZR1Nbw008/dSM4XHvttS6QVOD4zDPPcEAARF6wqF57Gtdw/fr1rv2iF4B5d0iJa6gHAEgrlAnxSzuH6mCnttnffvttMpQMAMI8dI6udtXmT72Ts2fP7m7Bp1vvKd3i18MYAAAAaaRmUYPKqieyevwpjaJJt99T2xu1r1m3bl14SgoAAICUX7OoNHPOnDnd3woYvXszqwOMOpgAAAAgDdcs6jZUP/74o0tFa4zDUaNGuaEdNJishtUBAABAGg4WBwwY4IaAkCFDhtgdd9zh7ohy5ZVX2qxZs8JRRgAAAERKsNikSZOov8uWLevu3azb7uXNmzeqRzQAAADS+DiLMcdABAAAQBoNFu++++4Ev+CcOXP+TnkAAAAQab2hc+fOHTXlypXL3eP0+++/j3p+zZo1bp6eBwAAQBqrWXznnXei/tbtp+677z6bOHGiZciQIWo4nUceecQFkgAAAEjD4yxOmTLFnnzyyahAUfR3nz593HMAAABIw8Hi+fPnXQ/omDTv4sWLSVUuAAAARGJv6M6dO1uXLl1s27ZtVrt2bTdv1apVNnLkSPccAAAA0nCwOHr0aCtcuLC9/PLLtnfvXjevSJEi9tRTT9kTTzwRjjICAAAgUoLF9OnT29NPP+2mY8eOuXl0bAEAAEid/tag3ASJAAAAqVuiO7gAAAAg7SBYBAAAgC+CRQAAAPgiWAQAAEDSdnDRfaA1HThwINZA3NzFBQAAIA0Hi88//7wNGTLEatWq5cZXTJcuXXhKBgAAgMgLFidOnGhTp0619u3bh6dEAAAAiNw2i2fPnrV69eqFpzQAAACI7GDxwQcftBkzZoSnNAAAAIi8NHSfPn2i/laHlkmTJtnnn39u1157rWXKlCnasmPGjEn6UgIAACDlBovr1q2L9rh69eru/w0bNoSnVAAAAIicYPGLL74If0kAAAAQ+W0WH3jgATt+/His+SdPnnTPAQAAIA0Hi++++6799ddfseZr3rRp05KqXAAAAIikcRaPHTtmgUDATapZzJo1a9RzFy5csIULF1rBggXDVU4AAACk5GAxT5487m4tmq655ppYz2u+7u4CAACANBgsqpOLahUbNmxoH374oeXLly/qucyZM1vJkiWtaNGi4SonAAAAUnKwWL9+fff/jh07rESJEtwTGgAAIA1IULD4008/WZUqVSx9+vR29OhRW79+ve+yGqgbAAAAaShY1CDc+/btcx1Y9LfaJyolHZPmq7MLAAAA0lCwqNRzgQIFov4GAABA2pCgYFGdV0L9DQAAgNQtwR1cPOrc0qBBA9fhRf+XKVMmPCUDAABA5N3BZfjw4W5A7hdffNHKlStnxYsXt3/96182efJk+/XXX8NTSgAAAERGzaICQ02yd+9eW758uc2fP98eeeQRu3jxIh1cAAAA0nKwKKdOnbKvv/7ali1b5gbrXrdunRtaR2lpAAAApOFgsV69ei44rFixogsO+/btazfffLPlzZs3PCUEAABA5LRZ3Lx5s+XIkcMqVKjgJgWNBIoAAACpU6KDxT/++MOWLl1q//jHP+zTTz+1G264wYoVK2b333+/6+QCAACANBws6i4tuqVfr169bPbs2fbJJ5/Yrbfeah988IE9/PDD4SklAAAAIqPN4tq1a13HFk3q5HL8+HGrWrWq9ezZ0429CAAAgDQcLNauXdtq1KjhAsOuXbu6zi25c+cOT+kAAAAQWcHi4cOHLVeuXOEpDQAAACK7zSKBIgAAQNqR6GARAAAAaQfBIgAAAHwRLAIAAMAXwSIAAACSrjf0hQsXbOrUqbZkyRI7cOCAXbx4MdrzursLAAAA0miw+Nhjj7lgsXnz5lalShV3RxcAAACkTokOFmfOnGnvv/++NWvWLDwlAgAAQOS2WcycObOVLVs2PKUBAABAZAeLTzzxhL366qsWCATCUyIAAABEbhr666+/ti+++MI++eQTq1y5smXKlCna83PmzEnK8gEAACCSgsU8efLYXXfdFZ7SAAAAILKDxXfeeSc8JQEAAEDkB4uegwcP2pYtW9zf5cuXtwIFCiRluQAAABCJHVxOnjxpDzzwgBUpUsRuvvlmNxUtWtS6dOlip06dCk8pAQAAEBnBYp8+fWz58uX2n//8x44cOeKmjz/+2M1TT2kAAACk4TT0hx9+aLNnz7YGDRpEzdMA3dmyZbP77rvP3njjjaQuIwAAACKlZlGp5kKFCsWaX7BgQdLQAAAAaT1YrFu3rg0ePNhOnz4dNe+vv/6y559/3j0HAACANJyG1t1bmjRpYldddZVVq1bNzfvxxx8ta9as9umnn4ajjAAAAIiUYLFKlSr266+/2vTp023z5s1uXtu2ba1du3au3SIAAADS+DiL2bNnt65duyZ9aQAAABB5weK8efPs9ttvd/eB1t9xadmypYXD4cOHrWfPnm7InvTp09s999zjUuJXXHFFnOutXLnS+vfvb6tWrbIMGTJY9erVXbqcWlAAAIAkChbvvPNO27dvn+vxrL/9pEuXzi5cuGDhoDT33r17bfHixXbu3Dnr3LmzdevWzWbMmBFnoNi0aVPr16+fjR8/3jJmzOjaVyrYBIBwmjBhgr300kvu3Kn23ToH1a5dO971Zs6c6Zr2tGrVyubOnctBAhAZweLFixdD/p1cNm3aZIsWLbLvvvvOatWq5ebpxKvxHUePHu3uIBNK7969rVevXta3b9+oebo1YVzOnDnjJs+xY8eSbDsApA2zZs1yNzCYOHGi1alTx8aOHes6BuoWqbro9vPbb7/Zk08+aTfddFOylhcA4pLoKrZp06ZFC6Y8Z8+edc+Fg2oI8+TJExUoSuPGjV0NodLLoRw4cMA9pxNzvXr13NiQ9evXt6+//jrO9xoxYoTlzp07aipevHiSbw+A1G3MmDGuXbcyIJUqVXJBo9p6T5kyxXcdZWWUQdEwZKVLl07W8gJAkgaLOvkdPXo01vzjx4+758LBS4EHU0o5X7587rlQtm/f7v5/7rnn3ElbNZPXXXedNWrUyPXm9qOUtbbPm3bt2pXEWwMgNdOF85o1a9wFrUcXtnqsC18/Q4YMcee5Ll26JOh9dNGuzEfwBAApIlgMBAKubWJM//3vf11NXGIoPazXimvyhudJLC9d/tBDD7kgtkaNGvbKK6+4NHRcV/dZsmSxXLlyRZsAIKEOHTrkaglj3ulKj/0ubpXxePvtt23y5MkJfh+yIABS3NA5Cra8AE61c6rZ8+jEuGPHDteZJDGeeOIJ69SpU5zLKB1TuHBhl1YOdv78eddDWs+FUqRIEfe/UkDBKlasaDt37kxUOQEgXJSVad++vQsU8+fPn+D1lAVRu0iPahZpNgPgsgaLXi/oH374wTXUDh6yJnPmzFaqVCk3nE1iFChQwE3x0W0Ejxw54lI7NWvWdPOWLl3qag/VeDwUlUcdX9SgPNgvv/zihgECgHBQwKdhuvbv3x9tvh6Hurjdtm2b69jSokWLWJkRXZTrHFamTJmQWRBNAJBigkXdD9oLwlq3bu1u75dcVBuoWku1PVRDcQ2d06NHD2vTpk1UT+jdu3e7Gk91stHwFKoBfeqpp1y5NWyFxld89913XVp79uzZyVZ2AGmLLp51UbtkyZKoi2wFf3qs81ZMFSpUsPXr10ebN2DAAFfjqLFkqS0EEHF3cOnYsaNdDrq9oE60Cgi9QbnHjRsX9bwCSF2Bnzp1Kmre448/bqdPn3ZD6ChlraBR4zSGukoHgKSi9LDOlRrBQRevGjrn5MmTUZ0AO3ToYMWKFXPtDnXhrduoBtPoDxJzPgBERLCo9onqKPL++++7tn/q+RdMQVk4qOdzXANwq8ZTnW9CdaIJHmcRAMJN2ZeDBw/aoEGDXKcWZTY0IoPX6UXnTm4OACDVBosaA+ytt95ynVOUKtGt9NTeRnca0IkRAGAuExIq7SzLli2LcxdNnTqVXQggcofOUTpYvfYULKrxtW5LpeBRgeK3334bnlICAAAgMoJFpVSqVq3q/laPaG+A7jvuuMMWLFiQ9CUEAABA5ASLV111le3du9f9rY4in332mftb921mGAeEw8Wzp+33F+9wk/4GAAApOFi866673BAQ0rNnTxs4cKCVK1fO9e574IEHwlFGAAAAREoHl5EjR0br8VeiRAl3v1MFjMGDygIAACANBouh7q6iCQAAAGk0WJw3b16CX7Bly5Z/pzwAAACItGDRu2VVfHSLPQ3aDQAAgDQULHo3tQcAAEDakuje0MF032UAAACkXokOFpVmHjp0qBUrVswNyr19+3Y3X0PovP322+EoIwAAACIlWBw2bJi7b+moUaMsc+bMUfOrVKnibvsHAACANBwsTps2zSZNmmTt2rWzDBkyRM2vVq2abd68OanLBwAAgEgKFnfv3m1ly5YN2Qnm3LlzSVUuAAAARGKwWKlSJfvqq69izZ89e7bVqFEjqcoFAACASLyDy6BBg6xjx46uhlG1iXPmzLEtW7a49PT8+fPDU0oAAABERs1iq1at7D//+Y99/vnnliNHDhc8btq0yc279dZbw1NKAAAApPyaxfPnz9vw4cPtgQcesMWLF4evVAAAAIi8msWMGTO6IXMUNAIAACD1S3QaulGjRrZ8+fLwlAYAAACR3cHl9ttvt759+9r69eutZs2art1isJYtWyZl+QAAABBJweIjjzzi/h8zZkys59KlS+duBwgAAIA0GixquBwAAACkDYlqs6g7tKiTy4YNG8JXIgAAAERmsJgpUyYrUaIEqWYAAIA0ItG9ofv372/PPvusHT58ODwlAgAAQOS2WXzttdds69atVrRoUStZsmSs3tBr165NyvIBAAAgkoLFO++8MzwlAQAAQOQHi4MHDw5PSQAAAJDiJDpY9KxZs8Y2bdrk/q5cubLVqFEjKcsFAACASAwWDxw4YG3atLFly5ZZnjx53LwjR47YLbfcYjNnzrQCBQqEo5wAAACIhN7QPXv2tOPHj9vGjRtdj2hNGnfx2LFj1qtXr/CUEgAAAJFRs7ho0SL7/PPPrWLFilHzKlWqZBMmTLDbbrstqcsHAACASKpZ1O3+NDh3TJrHrQABAADSeLDYsGFDe+yxx2zPnj1R83bv3m29e/e2Ro0aJXX5AAAAEEnBogblVvvEUqVKWZkyZdx09dVXu3njx48PTykBAAAQGW0Wixcv7u7SonaLmzdvdvPUfrFx48bhKB8AAAAibZzFdOnS2a233uomAAAApF4JTkMvXbrU9XpWujmmo0ePuoG5v/rqq6QuHwAAACIhWBw7dqx17drVcuXKFeu53Llz20MPPWRjxoxJ6vIhyMWzp+33F+9wk/4GAABIMcHijz/+aE2bNvV9XmMs6haAAAAASIPB4v79+0OOr+jJmDGjHTx4MKnKBQAAgEgKFosVK+Zu6+fnp59+siJFiiRVuQAAABBJwWKzZs1s4MCBdvp07LZyf/31lw0ePNjuuOOOpC4fAAAAImHonAEDBticOXPsmmuusR49elj58uXdfI21qPtCX7hwwfr37x/OsgIAACClBouFChWyFStWWPfu3a1fv34WCASixlxs0qSJCxi1DICkoR7vu175p/u7eO/Zlj5zVnYtACBlD8pdsmRJW7hwof3555+2detWFzCWK1fO8ubNG74SAgAAILLu4KLg8Prrr0/60gAAACAyO7gAAAAg7SFYBAAAgC+CRQAAAPgiWAQAAIAvgkUAAAD4IlgEAABA0g6dA4Ty28jmYdkxJ0+etCte+d/fm4Y2tRw5cnAAAABIJtQsAgAAwBfBIgAAAHwRLAIAAMAXwSIAAAB8ESwCAADAF8EiAAAAfBEsAgAAwBfBIgAAAHwRLAIAAMAXwSIAAAB8ESwCSFEunj1tv794h5v0NwDg8iJYBAAAgC+CRQAAAPgiWAQAAIAvgkUAAABEfrB4+PBha9euneXKlcvy5MljXbp0sRMnTsS5zr59+6x9+/ZWuHBhy5Ejh1133XX24YcfJluZAQAAIl3EBIsKFDdu3GiLFy+2+fPn25dffmndunWLc50OHTrYli1bbN68ebZ+/Xq7++677b777rN169YlW7kBAAAiWUQEi5s2bbJFixbZW2+9ZXXq1LEbb7zRxo8fbzNnzrQ9e/b4rrdixQrr2bOn1a5d20qXLm0DBgxwtZJr1qxJ1vIDSHsmTJhgpUqVsqxZs7rz1urVq32XnTx5st10002WN29eNzVu3DjO5QEgOUVEsLhy5UoX5NWqVStqnk6m6dOnt1WrVvmuV69ePZs1a5ZLYV+8eNEFl6dPn7YGDRr4rnPmzBk7duxYtAkAEkPnnT59+tjgwYNt7dq1Vq1aNWvSpIkdOHAg5PLLli2ztm3b2hdffOHOd8WLF7fbbrvNdu/ezY4HcNlFRLCotocFCxaMNi9jxoyWL18+95yf999/386dO2dXXnmlZcmSxR566CH76KOPrGzZsr7rjBgxwnLnzh016aQNAIkxZswY69q1q3Xu3NkqVapkEydOtOzZs9uUKVNCLj99+nR75JFHrHr16lahQgWXRdEF7pIlS9jxANJ2sNi3b19Lly5dnNPmzZsv+fUHDhxoR44csc8//9y+//57d6WvNotqv+inX79+dvTo0ahp165dl/z+ANKes2fPuqYuyn54lAXRY9UaJsSpU6fcha4uiP2QBQGQXDLaZfTEE09Yp06d4lxGbQ3Vmzlm+ub8+fMuvaznQtm2bZu99tprtmHDBqtcubKbp1TQV1995doS6Uo/FNVAagKAS3Ho0CG7cOGCFSpUKNp8PU7oxe8zzzxjRYsWjRZwhsqCPP/88xwkAKk7WCxQoICb4lO3bl1XQ6ir9Zo1a7p5S5cudWkaNRz3uzL3ruiDZciQwa0XTr+NbB6W1z158qRd8cr//t40tKkbDghA6jJy5EjXvlrtGNU5Jq4siLIlHrWvptkMgFQXLCZUxYoVrWnTpq4NkGoElZ7p0aOHtWnTxl19ixqCN2rUyKZNm+Z6P6vdj9omqp3i6NGjXbvFuXPnRg29AyQVLg4QLH/+/O6idP/+/dHm67FfJsSjc5WCRTWdufbaa+NcliwIgOQSER1cvAbgCgAVEDZr1swNnzNp0qSo5xVAakxFr0YxU6ZMtnDhQldz2aJFC3fiVSD57rvvuvUBIBwyZ87sMiDBnVO8zirKkvgZNWqUDR061A0TFjzyAwBcbhFRsyhq6D1jxgzf5zWeWSAQiDavXLly3LEFQLJTerhjx44u6FOmY+zYsa4ZiXpHezcMKFasmGt3KC+++KINGjTIneN0LvNGebjiiivcBACXU8QEiwAQKVq3bm0HDx50AaACPw2JoxpDr9PLzp07o7WnfuONN1wv6n/+85/RXkfjND733HPJXn4ACEawCABhoHbVmkJR55Vgv/32G8cAQIoVMW0WAQAAkPwIFgEAAOCLYBEAAAC+CBYBAADgi2ARAAAAvggWAQAA4ItgEQAAAL4IFgEAAOCLYBEAAAC+CBYBAADgi9v9AUi030Y2D9teO3nypF3xyv/+3jS0qeXIkSNs7wUAiB81iwAAAPBFsAgAAABfBIsAAADwRbAIAAAAXwSLAAAA8EWwCAAAAF8EiwAAAPBFsAgAAABfBIsAAADwRbAIAAAAXwSLAAAA8EWwCAAAAF8EiwAAAPBFsAgAAABfBIsAAADwRbAIAAAAXwSLAAAA8EWwCAAAAF8Z/Z8CcDnlyJHDAoEABwEAcFlRswgAAABfBIsAAADwRbAIAAAAXwSLAAAA8EWwCAAAAF8EiwAAAPBFsAgAAABfBIsAAADwRbAIAAAAXwSLAAAA8EWwCAAAAF8EiwAAAPBFsAgAAABfGf2fAlKGHDlyWCAQuNzFAAAgTaJmEQAAAL4IFgEAAOCLYBEAAAC+CBYBAADgi2ARAAAAvggWAQAA4ItgEQAAAL4IFgEAAOCLYBEAAAC+CBYBAADgi9v9RRBuewcAAJIbNYsAAADwRc0igBSFGnQASFmoWQQAAIAvgkUAAAD4IlgEAACAL4JFAAAA+CJYBAAAgC+CRQAAAPgiWAQAAIAvgkUAAAD4IlgEAACAL4JFAAAA+CJYBAAAgC+CRQAAAPgiWAQAAEDkB4vDhg2zevXqWfbs2S1PnjwJWicQCNigQYOsSJEili1bNmvcuLH9+uuvYS8rAABAahExweLZs2ft3nvvte7duyd4nVGjRtm4ceNs4sSJtmrVKsuRI4c1adLETp8+HdayAgAApBYREyw+//zz1rt3b6tatWqCaxXHjh1rAwYMsFatWtm1115r06ZNsz179tjcuXPDXl4AaduECROsVKlSljVrVqtTp46tXr06zuU/+OADq1Chglte57mFCxcmW1kBIFUEi4m1Y8cO27dvn0s9e3Lnzu1O2itXrvRd78yZM3bs2LFoEwAkxqxZs6xPnz42ePBgW7t2rVWrVs1lNQ4cOBBy+RUrVljbtm2tS5cutm7dOrvzzjvdtGHDBnY8gMsu1QaLChSlUKFC0ebrsfdcKCNGjHBBpTcVL1487GUFkLqMGTPGunbtap07d7ZKlSq5pjBqbz1lypSQy7/66qvWtGlTe+qpp6xixYo2dOhQu+666+y1115L9rIDQEwZ7TLq27evvfjii3Eus2nTJpeaSS79+vVzNQKeo0ePWokSJahhBNI4L8ugJi7xta9es2aNO5d40qdP77IcflkNzQ8+74hqIuNqMqMsiKbgc1VwOQGkTccSeK6KmGDxiSeesE6dOsW5TOnSpS/ptQsXLuz+379/v+sN7dHj6tWr+66XJUsWN8Xc6dQwApDjx4+7rIOfQ4cO2YULF0JmNTZv3hxyHWU7LiULorbcMXGuAiB//PFHnOeqiAkWCxQo4KZwuPrqq13AuGTJkqjgUIGfekUnpkd10aJFbdeuXZYzZ05Lly6dXW7aBv0YqEy5cuWytILt5nhfbrpKV6Coc0JKEDMLcuTIEStZsqTt3LkzyX4gLqfU9J1PTdsibE/K5mVE8+XLl2SveVmDxcTQCfDw4cPuf121//DDD25+2bJl7YorrnB/K12tq+277rrLBXaPP/64vfDCC1auXDkXPA4cONCd6NVwPKGUPrrqqqsspdEJJzWcdBKL7U5bUtrxTkgQlj9/fsuQIYPLYgTTYy/jEZPmJ2b5UFmQ4DKmpH2W2j4Df0dq2hZhe1I2xS9J9loWITS4do0aNVzvwhMnTri/NX3//fdRy2zZsiWq3Y48/fTT1rNnT+vWrZtdf/31br1Fixa5oSkAIBwyZ85sNWvWdFkNz8WLF93junXrhlxH84OXl8WLF/suDwDJKWJqFqdOneqmuMRszKnaxSFDhrgJAJKL0sMdO3a0WrVqWe3atd2YrydPnnS9o6VDhw5WrFgxlwmRxx57zOrXr28vv/yyNW/e3GbOnOkuhCdNmsRBA3DZRUywiP9R2km1q6HST6kZ283xjiStW7e2gwcPuoyIOqmo3bSyGl4nFjWnCU4R6VamM2bMcDcRePbZZ13TGfWErlKlSpr9jqSm7UlN2yJsT9o7PukCSdm3GgAAAKlKxLRZBAAAQPIjWAQAAIAvgkUAAAD4IlgEAACAL4JFAAAA+CJYTCYaPkMDhOte1+rOrls/tWjRItZAvOG0ceNGu+eee6xUqVJuDEqN/ZYWtnvy5Ml20003Wd68ed3UuHFjW716darf7jlz5rhx/vLkyWM5cuRww7f8+9//TvXbHUzjFeqznpi7NqVkEyZMcN9f3VigTp068X6OP/jgA3dnKy1ftWpVW7hwoUXq9lyO73E4j01K/4wmdnt0u8lHH33UihQp4r7711xzTYr6vCV2e/T7WL58ecuWLZs7j/Xu3dtOnz5tl9uXX37pzqm6G50+NxpiKz7Lli2z6667zh0X3fUuvjGrQ9LQOQivHTt2BIoWLRqoVKlSYPbs2YEtW7YENmzYEHj55ZcD5cuXT7bdv3r16sCTTz4ZeO+99wKFCxcOvPLKK2liu++///7AhAkTAuvWrQts2rQp0KlTp0Du3LkD//3vf1P1dn/xxReBOXPmBH7++efA1q1bA2PHjg1kyJAhsGjRolS93cHlKVasWOCmm24KtGrVKhDpZs6cGcicOXNgypQpgY0bNwa6du0ayJMnT2D//v0hl//mm2/c8R41apT7DAwYMCCQKVOmwPr16wORuD3J/T0O57ak9M9oYrfnzJkzgVq1agWaNWsW+Prrr912LVu2LPDDDz8EInF7pk+fHsiSJYv7X9vy6aefBooUKRLo3bt34HJbuHBhoH///u7crhDuo48+inP57du3B7Jnzx7o06ePOw+MHz/+kn4HCBaTwe233+5OCCdOnIj13J9//hn19++//x5o2bJlIEeOHIGcOXMG7r333sC+ffvcc/rh1QdDJ8lgY8aMCZQuXTrRZSpZsmTYg8WUuN1y/vx59z7vvvtuIC1tt9SoUcMFDal9u3WM69WrF3jrrbcCHTt2TFE/xJeqdu3agUcffTTq8YULF1xwPmLEiJDL33fffYHmzZtHm1enTp3AQw89FIjE7Unu73G4tyUlf0YTuz1vvPGG+36ePXs2kBIldnu0bMOGDaPN69OnT+CGG24IpCQJCRaffvrpQOXKlaPNa926daBJkyaJei/S0GF2+PBhd+cGVc8rFRiTUoTevWNbtWrlll++fLm7L+z27dvdnSBEVfpKKU6fPj3a+np8//33W0qTkrf71KlTdu7cOcuXL5+lle3WeUWpYN0//eabb7bUvt26xWfBggWtS5culhqcPXvW1qxZ41KvHt0BRo9XrlwZch3ND15emjRp4rt8St+e5PweJ8e2pNTP6KVsz7x589x9zPX9112KdOeh4cOH24ULFywSt0d3VNI6Xqp6+/btLqXerFkzizRJdh64pHAWCbZq1SoX/avKOC6fffaZqxreuXNn1DxVl2tdpY9FNYFlypSJet6vFiYl1Cym1O2W7t27u6vgv/76K5Dat/vIkSOuBi9jxowurfL2228HwiElbfdXX33lajgPHjzoHqe0WptLsXv3brcPVqxYEW3+U0895WpNQlHKecaMGdHmKY1bsGDBQCRuT3J+j8O9LSn5M3op26NmJjq/PPDAA4Hvv//epX3z5csXeO655wKR+ll79dVX3XdI504zCzz88MOBlCYhNYvlypULDB8+PNq8BQsWuHVPnTqV4PeiZjHMEno3xU2bNrlGtJo8lSpVcjUyek7atGljv/32m3377bdRtS1qtKoG7ClNSt3ukSNHugblH330kWvonNq3O2fOnPbDDz/Yd999Z8OGDbM+ffq4xs6pdbuPHz9u7du3d50h8ufPf8nbg5Qt3N/jcEqNn1FlDFRLOmnSJKtZs6bLFPTv398mTpxokUjnSNWMvv7667Z27VrXWXDBggU2dOhQS6syXu4CpHblypVzPZY2b978t1+rcOHC1rBhQ5sxY4b94x//cP93797dUqKUuN2jR492PzKff/65XXvttZYWtlvpFvV+E/WGVkA2YsQIa9CggaXG7d62bZsLNNVbMPiHTDJmzOjS8GXKlLFIo6AiQ4YMtn///mjz9Vj7KxTNT8zyKX17kvN7HM5tSemf0Us5NuoBnSlTJreep2LFim50BKWBM2fObJG0PQMHDnQB/YMPPugeV61a1U6ePGndunVzQbDOq5HC7zyQK1cu19M7oSJniyOU2tOofYC67evDFmq4Ae+LtWvXLjd5fv75Z/e8al487dq1s1mzZrn2BmpHoVqYlCilbfeoUaPcVaHa1alNXFrZ7pj0o3TmzBlLrdut2sf169e72lRvatmypd1yyy3u7+AazUiiH1vV2AQPQaRjqcdqKxaK5sccskhtRP2WT+nbk5zf43BuS0r/jF7Ksbnhhhts69atUUGv/PLLLy6IvJyB4qVuj9rDxgwIM/xfIJzQLEpKkWTngUtOliPBtm3b5oaq8YYU+eWXX1wXdrWJqFChglvm4sWLgerVq7shFNasWePagNWsWTNQv379aK917NixQLZs2QLVqlULNGrUKGq+ho9QuxGt50fDG2jYCU0aBkDD6OjvX3/9NVVv98iRI92wCSrD3r17o6bjx4+n6u1WOxW1EVR59P6jR4927W8mT56cqrc7ppTUHuzvUDswtQubOnWq26/dunVzw394Pcnbt28f6Nu3b7Shc3S8ddzV3nPw4MEpbuicxGxPcn+Pw7ktKf0zmtjtURtk9Uzv0aOHa2M8f/581zb2hRdeCETi9ui7ou3RMHMaeuazzz5z7ag1wsDlps+79zuuEE4jRehvjTIh2g5tT8yhc9RGU+cBtVtm6JwUbM+ePa47vjqW6ISnxs0aPkRj4SVkSJFg+sDqQ6IxozwaC0rzgl8vJm+ZmFPMH+rUtt1671DbrRNCat5ujcVVtmzZQNasWQN58+YN1K1b1500wyklbHdK/yH+OzRGWokSJdy+VeP8b7/9Nuo5fY+1rcHef//9wDXXXOOW1/AZatgeqdtzOb7H4Tw2Kf0zmtjtUQcSDc2koEwdj4YNG+aGB4rE7Tl37pzrnKMAUefP4sWLBx555JFoQ4BdLjr3hfoeeOXX/zF/07WOLtK17To277zzTqLfN53+SdpKTwAAAKQWtFkEAACAL4JFAAAA+CJYBAAAgC+CRQAAAPgiWAQAAIAvgkUAAAD4IlgEAACAL4JFAAAA+CJYBAAAgC+CRQAAAPgiWAQAAID5+f92hRmkgN8IegAAAABJRU5ErkJggg==",
      "text/plain": [
       "<Figure size 640x480 with 2 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "labels = mod.design_sal_labels_['covariate']\n",
    "n_sig = len(sig_lvs)\n",
    "\n",
    "fig = plt.figure(constrained_layout=True)\n",
    "subfigs = fig.subfigures(nrows=n_sig)\n",
    "if n_sig == 1:\n",
    "    subfigs = [subfigs]\n",
    "\n",
    "for plot_idx, lv_idx in enumerate(sig_lvs):\n",
    "    subfig = subfigs[plot_idx]\n",
    "    subfig.suptitle('LV %s (%.2f%% variance explained, p = %.3f)' % (\n",
    "        lv_idx,\n",
    "        100*mod.variance_explained_[lv_idx],\n",
    "        mod.pvals_[lv_idx]))\n",
    "    ax = subfig.subplots(ncols=2)\n",
    "    ax[0].bar(x=labels,\n",
    "              height=mod.boot_stat_val_[:, lv_idx],\n",
    "              yerr=mod.get_boot_stat_yerr(lv_idx))\n",
    "    ax[0].set_ylabel('Correlation with data score')\n",
    "    ax[1].plot(mod.data_sals_z_[:, lv_idx])\n",
    "    ax[1].set_ylabel('Bootstrap ratio (z score)')\n",
    "    ax[1].set_xlabel('Observed variable')\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f439ef29-2c57-4099-9ba6-a303bef3f561",
   "metadata": {},
   "source": [
    "As we can see, the boostrap ratios are roughly sinusoidal. `mod.flip_signs(lv_idx=0)` can be used to obtain an equivalent solution in which the covariates have positive correlations with the latent variable."
   ]
  }
 ],
 "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.13.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}