{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "header",
   "metadata": {},
   "source": [
    "# The QIF Mean-Field Model\n",
    "\n",
    "An exact mean-field reduction of spiking neural populations using the Montbrio-Pazo-Roxin framework."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "learning-objectives",
   "metadata": {},
   "source": [
    "## Learning Objectives\n",
    "\n",
    "By the end of this tutorial, you will be able to:\n",
    "\n",
    "- Understand the QIF mean-field reduction and why it's exact\n",
    "- Explain the role of the Lorentzian distribution of excitabilities\n",
    "- Simulate population firing rate and mean membrane potential dynamics\n",
    "- Analyze bifurcations in the QIF model as a function of external input"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "background",
   "metadata": {},
   "source": [
    "## Background / Theory\n",
    "\n",
    "### From Spiking Neurons to Mean Fields\n",
    "\n",
    "Most neural mass models (Wilson-Cowan, Jansen-Rit) are phenomenological - they capture qualitative features but aren't derived from spiking neuron dynamics.\n",
    "\n",
    "The **Quadratic Integrate-and-Fire (QIF)** mean-field model is different: it's an **exact** reduction of a population of QIF neurons with heterogeneous excitabilities.\n",
    "\n",
    "### The Montbrio-Pazo-Roxin Reduction\n",
    "\n",
    "Consider $N$ all-to-all coupled QIF neurons:\n",
    "\n",
    "$$\n",
    "\\tau \\dot{V}_j = V_j^2 + \\eta_j + I(t) + J\\tau r(t)\n",
    "$$\n",
    "\n",
    "where $\\eta_j$ is drawn from a Lorentzian distribution with mean $\\bar{\\eta}$ and width $\\Delta$.\n",
    "\n",
    "Remarkably, in the limit $N \\to \\infty$, this reduces exactly to:\n",
    "\n",
    "$$\n",
    "\\begin{aligned}\n",
    "\\tau\\dot{r} &= \\frac{\\Delta}{\\pi} + 2\\tau r v \\\\\n",
    "\\tau\\dot{v} &= v^2 + \\bar{\\eta} + I(t) + J\\tau r - (\\pi\\tau r)^2\n",
    "\\end{aligned}\n",
    "$$\n",
    "\n",
    "where:\n",
    "- $r$: Population firing rate\n",
    "- $v$: Mean membrane potential\n",
    "\n",
    "### Key Parameters\n",
    "\n",
    "| Parameter | Symbol | Description |\n",
    "|-----------|--------|-------------|\n",
    "| Time constant | $\\tau$ | Membrane time constant |\n",
    "| Mean excitability | $\\bar{\\eta}$ | Center of Lorentzian distribution |\n",
    "| Heterogeneity | $\\Delta$ | Width of Lorentzian (HWHM) |\n",
    "| Coupling | $J$ | Recurrent synaptic strength |\n",
    "\n",
    "### Why Is This Important?\n",
    "\n",
    "1. **Exact reduction**: No approximations needed (unlike firing rate models)\n",
    "2. **Captures heterogeneity**: The $\\Delta$ parameter models neuron-to-neuron variability\n",
    "3. **Biophysically grounded**: Directly relates to spiking neuron properties\n",
    "4. **Low-dimensional**: Only 2 variables describe the entire population"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "implementation",
   "metadata": {},
   "source": [
    "## Implementation\n",
    "\n",
    "### Step 1: Setup and Imports"
   ]
  },
  {
   "cell_type": "code",
   "id": "imports",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2026-03-21T13:55:36.019457100Z",
     "start_time": "2026-03-21T13:55:30.297721700Z"
    }
   },
   "source": [
    "import brainstate\n",
    "import brainunit as u\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "\n",
    "import brainmass\n",
    "\n",
    "# Set simulation time step\n",
    "brainstate.environ.set(dt=0.01 * u.ms)  # Fine timestep for accuracy"
   ],
   "outputs": [],
   "execution_count": 1
  },
  {
   "cell_type": "markdown",
   "id": "single-pop",
   "metadata": {},
   "source": [
    "### Step 2: Single Population Simulation\n",
    "\n",
    "Create a QIF population and observe its dynamics:"
   ]
  },
  {
   "cell_type": "code",
   "id": "create-model",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2026-03-21T13:55:36.720975Z",
     "start_time": "2026-03-21T13:55:36.028457500Z"
    }
   },
   "source": [
    "# Create QIF model\n",
    "model = brainmass.MontbrioPazoRoxinStep(\n",
    "    in_size=1,  # single population\n",
    "    tau=1.0 * u.ms,  # time constant\n",
    "    eta=-5.0,  # mean excitability (negative = subthreshold)\n",
    "    delta=1.0 * u.Hz,  # heterogeneity (Lorentzian width)\n",
    "    J=15.0,  # recurrent coupling strength\n",
    ")\n",
    "\n",
    "# Initialize states\n",
    "brainstate.nn.init_all_states(model)\n",
    "\n",
    "# Set non-zero initial firing rate\n",
    "model.r.value = np.array([0.1]) * u.Hz\n",
    "model.v.value = np.array([-2.0])\n",
    "\n",
    "print(f\"tau = {model.tau.value()}\")\n",
    "print(f\"eta = {model.eta.value()}\")\n",
    "print(f\"delta = {model.delta.value()}\")\n",
    "print(f\"J = {model.J.value()}\")"
   ],
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "tau = 1. ms\n",
      "eta = -5.0\n",
      "delta = 1. Hz\n",
      "J = 15.0\n"
     ]
    }
   ],
   "execution_count": 2
  },
  {
   "cell_type": "code",
   "id": "simulate",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2026-03-21T13:55:37.037155Z",
     "start_time": "2026-03-21T13:55:36.824897300Z"
    }
   },
   "source": [
    "# Define simulation step\n",
    "def step_run(i):\n",
    "    with brainstate.environ.context(i=i, t=i * brainstate.environ.get_dt()):\n",
    "        r = model.update()\n",
    "        return model.r.value, model.v.value\n",
    "\n",
    "\n",
    "# Run simulation (100 ms)\n",
    "n_steps = int(100 * u.ms / brainstate.environ.get_dt())\n",
    "indices = np.arange(n_steps)\n",
    "r_trace, v_trace = brainstate.transform.for_loop(step_run, indices)"
   ],
   "outputs": [],
   "execution_count": 3
  },
  {
   "cell_type": "code",
   "id": "plot-dynamics",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2026-03-21T13:55:37.590828600Z",
     "start_time": "2026-03-21T13:55:37.055955200Z"
    }
   },
   "source": [
    "t_ms = indices * brainstate.environ.get_dt()\n",
    "\n",
    "fig, axes = plt.subplots(1, 2, figsize=(12, 4))\n",
    "\n",
    "# Firing rate\n",
    "axes[0].plot(t_ms, r_trace[:, 0], 'b-', linewidth=1.5)\n",
    "axes[0].set_xlabel('Time (ms)')\n",
    "axes[0].set_ylabel('Firing rate r (Hz)')\n",
    "axes[0].set_title('Population Firing Rate')\n",
    "\n",
    "# Mean membrane potential\n",
    "axes[1].plot(t_ms, v_trace[:, 0], 'r-', linewidth=1.5)\n",
    "axes[1].set_xlabel('Time (ms)')\n",
    "axes[1].set_ylabel('Mean potential v')\n",
    "axes[1].set_title('Mean Membrane Potential')\n",
    "\n",
    "plt.tight_layout()\n",
    "plt.show()"
   ],
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<Figure size 1200x400 with 2 Axes>"
      ],
      "image/png": "iVBORw0KGgoAAAANSUhEUgAABKUAAAGGCAYAAACqvTJ0AAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjgsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvwVt1zgAAAAlwSFlzAAAPYQAAD2EBqD+naQAAdytJREFUeJzt3XlYFeX///HXAQVUBDcEN8Ct1DQxTEL75EZC2UKamVkuEWa5Y5aYaWZFi2tqbqVlZZplZmaWqVkmam6VmaTmFgpqJigqKMzvD76cnydAOUfgLD4f1zUX58zcM/Oe+1hzX++573tMhmEYAgAAAAAAAEqRm70DAAAAAAAAwPWHpBQAAAAAAABKHUkpAAAAAAAAlDqSUgAAAAAAACh1JKUAAAAAAABQ6khKAQAAAAAAoNSRlAIAAAAAAECpIykFAAAAAACAUkdSCgAAAAAAAKWOpBSAUtOuXTu1a9euWI/54osvymQyFesxi4s11/v999/LZDLp+++/L9GYAAAAbHHw4EGZTCZNmDDB3qHgKvJ+q/fee8/qfd977z2ZTCYdPHiw2OMCCkJSCnAReTeQvMXLy0s33HCDBg4cqNTUVHuHd03OnTunF1980eESNpfX9+VLQECAvUOzyn//7ZQpU0a1atVSnz59lJycbNMxHfU3AwC4lsvvYRs2bMi33TAM1alTRyaTSffcc48dIiy64OBgmUwmRUREFLh97ty55mvdunVrKUfnuvIecOYt5cuXV5MmTTR69Gilp6dbfbxXX31Vy5YtK/5AC7Bw4UJNmTKlVM4FlJQy9g4AQPF66aWXVLduXV24cEEbNmzQzJkztXLlSu3atUvly5e3d3g2OXfunMaNGydJ+XoejR49WiNHjrRDVLnuvPNO9erVy2JduXLlJEnffvttkY9zxx136Pz58/Lw8CjW+Kxx+b+dTZs26b333tOGDRu0a9cueXl5WXWsK/1mAAAUNy8vLy1cuFC33367xfr169fr77//lqenp50is46Xl5fWrVunlJSUfA+5PvroI3l5eenChQt2is61zZw5U97e3jp79qy+/fZbvfLKK1q7dq1++uknq3rlv/rqq3rwwQcVHR1dcsH+n4ULF2rXrl0aOnSoxfqgoCCdP39eZcuWLfEYgGtFUgpwMXfddZdatmwpSXriiSdUtWpVTZo0SV988YV69Ohh5+iKX5kyZVSmjP3+V3bDDTfo0UcfLXBbURJMFy5ckIeHh9zc3KxO/BS3//7bqVatml5//XUtX75cDz30kF1jAwDgSu6++24tWbJEb731lkW7YOHChQoNDdXJkyftGF3RtWnTRj///LMWL16sIUOGmNf//fff+vHHH/XAAw/os88+s2OEtsvIyFCFChXsHUahHnzwQVWrVk2S1L9/f3Xt2lVLly7Vpk2bFB4ebuforJM3agJwBgzfA1xchw4dJEkHDhyQJF26dEnjx49X/fr15enpqeDgYI0aNUqZmZkW+wUHB+uee+7Rt99+q5CQEHl5ealJkyZaunSpRbnC5nQqynj0rKwsjRkzRqGhofL19VWFChX0v//9T+vWrTOXOXjwoPz8/CRJ48aNM3etfvHFFws9v7XXuGHDBrVq1UpeXl6qV6+eFixYcIUaLbr/zimVN2/UokWLNHr0aNWqVUvly5dXenp6gXNKtWvXTk2bNtXu3bvVvn17lS9fXrVq1dIbb7yR71yHDh3SfffdpwoVKqh69eoaNmyYvvnmm2uap+p///ufJGn//v3mdcXxm0nSnj179OCDD6pKlSry8vJSy5YttXz5cpviBACgR48e+ueff7R69WrzuqysLH366ad65JFHCtwnJydHU6ZM0U033SQvLy/5+/vrySef1L///mtR7osvvlDnzp1Vs2ZNeXp6qn79+ho/fryys7Mtyllz3y6Ml5eXunTpooULF1qs//jjj1W5cmVFRkYWuF9R7qt5bbMNGzZo8ODB8vPzU6VKlfTkk08qKytLp0+fVq9evVS5cmVVrlxZzz77rAzDKPB8kydPVlBQkMqVK6e2bdtq165dFtv79Okjb29v7d+/X3fffbcqVqyonj17SpJ+/PFHdevWTYGBgfL09FSdOnU0bNgwnT9/vsBjJCcnKzo6Wt7e3vLz89MzzzyTr+6L+lta479t6IyMDA0fPlx16tSRp6enbrzxRk2YMMGijkwmkzIyMvT++++b2z99+vQxb09OTtbjjz8uf39/eXp66qabbtK8efMszpvXJvzkk0/0yiuvqHbt2vLy8lLHjh21b98+c7l27drpq6++0qFDh8znCg4OllTwnFK//vqr+vTpo3r16snLy0sBAQF6/PHH9c8//9hcR0BxoKcU4OLyEgpVq1aVlNsD5v3339eDDz6o4cOHa/PmzUpISNAff/yhzz//3GLfvXv3qnv37urfv7969+6t+fPnq1u3blq1apXuvPPOa44tPT1d77zzjnr06KHY2FidOXNG7777riIjI7VlyxaFhITIz89PM2fO1FNPPaUHHnhAXbp0kSTdfPPNhR7Xmmvct2+fHnzwQcXExKh3796aN2+e+vTpo9DQUN10001XvYYLFy7ke/pasWLFKw4TGD9+vDw8PPTMM88oMzPzij2q/v33X0VFRalLly566KGH9Omnn+q5555Ts2bNdNddd0nKbSR16NBBx44d05AhQxQQEKCFCxdaJIpskZdQrFy5snldcfxmv//+u9q0aaNatWpp5MiRqlChgj755BNFR0frs88+0wMPPHBNcQMArj/BwcEKDw/Xxx9/bL4/fv3110pLS9PDDz+st956K98+Tz75pN577z317dtXgwcP1oEDBzR9+nTt2LFDP/30k3no03vvvSdvb2/FxcXJ29tba9eu1ZgxY5Senq4333zT4phFuW9fzSOPPKJOnTpp//79ql+/vqTcHl8PPvhggcOxrL2vDho0SAEBARo3bpw2bdqkOXPmqFKlStq4caMCAwP16quvauXKlXrzzTfVtGnTfNMULFiwQGfOnNGAAQN04cIFTZ06VR06dNBvv/0mf39/c7lLly4pMjJSt99+uyZMmGCeRmLJkiU6d+6cnnrqKVWtWlVbtmzRtGnT9Pfff2vJkiUW58rOzlZkZKTCwsI0YcIEfffdd5o4caLq16+vp556yurf0hqXt6ENw9B9992ndevWKSYmRiEhIfrmm280YsQIJScna/LkyZKkDz74QE888YRatWqlfv36SZL5N0xNTdVtt90mk8mkgQMHys/PT19//bViYmKUnp6ebwjea6+9Jjc3Nz3zzDNKS0vTG2+8oZ49e2rz5s2SpOeff15paWn6+++/zef39vYu9HpWr16tv/76S3379lVAQIB+//13zZkzR7///rs2bdrksC8OwnXAAOAS5s+fb0gyvvvuO+PEiRPGkSNHjEWLFhlVq1Y1ypUrZ/z999/Gzp07DUnGE088YbHvM888Y0gy1q5da14XFBRkSDI+++wz87q0tDSjRo0aRosWLczrxo4daxT0v5K8eA4cOGBe17ZtW6Nt27bm75cuXTIyMzMt9vv3338Nf39/4/HHHzevO3HihCHJGDt2bL7z/Pf8tlzjDz/8YF53/Phxw9PT0xg+fHi+c/2XpAKX+fPnF3i969atMyQZ9erVM86dO2dxrLxt69atM69r27atIclYsGCBeV1mZqYREBBgdO3a1bxu4sSJhiRj2bJl5nXnz583GjVqlO+YBSno386nn35q+Pn5GZ6ensaRI0fMZYvjN+vYsaPRrFkz48KFC+Z1OTk5RuvWrY2GDRteMVYAAC6Xdw/7+eefjenTpxsVK1Y032O7detmtG/f3jCM3Ht+586dzfv9+OOPhiTjo48+sjjeqlWr8q3/7z3bMAzjySefNMqXL29xLyvqfbsweTFeunTJCAgIMMaPH28YhmHs3r3bkGSsX7/e4nrzFPW+mrdvZGSkkZOTY14fHh5umEwmo3///uZ1ly5dMmrXrm3Rjjlw4IAhydyuzLN582ZDkjFs2DDzut69exuSjJEjR+a7zoLqMyEhwTCZTMahQ4fyHeOll16yKNuiRQsjNDTU/N2a37IgeW3JpKQk48SJE8aBAweM2bNnG56enoa/v7+RkZFhLFu2zJBkvPzyyxb7Pvjgg4bJZDL27dtnXlehQgWjd+/e+c4TExNj1KhRwzh58qTF+ocfftjw9fU110tem7Bx48YWba6pU6cakozffvvNvK5z585GUFBQvnPl/VZ5bVLDKLjeP/7443xt4YLa8EBJYvge4GIiIiLk5+enOnXq6OGHH5a3t7c+//xz1apVSytXrpQkxcXFWewzfPhwSdJXX31lsb5mzZoWT9d8fHzUq1cv7dixQykpKdccq7u7u7mXUE5Ojk6dOqVLly6pZcuW2r59u03HtPYamzRpYh6mJkl+fn668cYb9ddffxXpfPfff79Wr15tsRTWtT5P7969zZOhX423t7fFnFUeHh5q1aqVRXyrVq1SrVq1dN9995nXeXl5KTY2tkjnyHP5v50HH3xQFSpU0PLly1W7dm1zmWv9zU6dOqW1a9fqoYce0pkzZ3Ty5EmdPHlS//zzjyIjI7V3716b3/gHALi+PfTQQzp//rxWrFihM2fOaMWKFYUO3VuyZIl8fX115513mu9FJ0+eVGhoqLy9vS16G19+z867d/3vf//TuXPntGfPHovjFuW+fTXu7u566KGH9PHHH0vKneC8Tp06Fu2VPLbcV2NiYix6xYSFhckwDMXExFjE0LJlywLjjo6OVq1atczfW7VqpbCwMHMb7HKX92bKc3l9ZmRk6OTJk2rdurUMw9COHTvyle/fv7/F9//9738WcVnzW17JjTfeKD8/P9WtW1dPPvmkGjRooK+++krly5fXypUr5e7ursGDB1vsM3z4cBmGoa+//vqKxzYMQ5999pnuvfdeGYZhEWdkZKTS0tLytaP69u1r0Zs+7/e35t/S5S6v97ye/rfddpsk2dzuBooDw/cAFzNjxgzdcMMNKlOmjPz9/XXjjTfKzS03/3zo0CG5ubmpQYMGFvsEBASoUqVKOnTokMX6Bg0a5OvKe8MNN0jKHdr137fC2OL999/XxIkTtWfPHl28eNG8vm7dujYdz9prDAwMzHeMypUrF3kOgtq1axf66ubCWHNttWvXzvcbVK5cWb/++qv5+6FDh1S/fv185f5bB1eT928nLS1N8+bN0w8//FDgMMRr+c327dsnwzD0wgsv6IUXXiiwzPHjxy0auwAAFIWfn58iIiK0cOFCnTt3TtnZ2XrwwQcLLLt3716lpaWpevXqBW4/fvy4+fPvv/+u0aNHa+3atUpPT7col5aWZvG9KPftonjkkUf01ltv6ZdfftHChQv18MMPFzi8ypb76n/bPr6+vpKkOnXq5FtfUHuoYcOG+dbdcMMN+uSTTyzWlSlTxuLBVp7Dhw9rzJgxWr58eb7j/7c+vby8zPNU5vlvO82a3/JKPvvsM/n4+Khs2bKqXbu2edidlNvWqlmzpipWrGixT+PGjc3br+TEiRM6ffq05syZozlz5hQpzv/+TnnTKdg6T9apU6c0btw4LVq0KN+5/lvvQGkiKQW4mFatWpnfoFaY4hwzXtix/jsBZUE+/PBD9enTR9HR0RoxYoSqV68ud3d3JSQkWEyuXZxx/Ze7u3uB641CJvYsDkXtJSWVbnyX/9uJjo7W7bffrkceeURJSUnmOQqu9TfLycmRJD3zzDOF9iizNpkGAECeRx55RLGxsUpJSdFdd92lSpUqFVguJydH1atX10cffVTg9rxEyOnTp9W2bVv5+PjopZdeUv369eXl5aXt27frueeeM9/X8hTXfTssLEz169fX0KFDdeDAgStO1i5Zd18tLMaC1l9Le8PT09P8YDRPdna27rzzTp06dUrPPfecGjVqpAoVKig5OVl9+vQpcn1erqi/5dXccccd5rfvFbe863r00UfVu3fvAsv8d77U4m4DPvTQQ9q4caNGjBihkJAQeXt7KycnR1FRUfnqHShNJKWA60hQUJBycnK0d+9e85MdKXfixdOnTysoKMiifN7Tt8sTPH/++ackmd/ukffU5vTp0xYNv6s9MZKkTz/9VPXq1dPSpUstzjF27FiLctYk0ay9RlcQFBSk3bt35/utLn9Di7XyEk3t27fX9OnTNXLkSEnX/pvVq1dPklS2bFmre5gBAHA1DzzwgJ588klt2rRJixcvLrRc/fr19d1336lNmzZXfFj0/fff659//tHSpUt1xx13mNfnvZGtJPXo0UMvv/yyGjdurJCQkALL2OO+unfv3nzr/vzzT3Pb8Ep+++03/fnnn3r//fctJlC//K2J1irqb3ktgoKC9N133+nMmTMWvaXyhm9e3r4sqA3k5+enihUrKjs7u1h/p6K2kf/991+tWbNG48aN05gxY8zrC/otgdLGnFLAdeTuu++WJE2ZMsVi/aRJkyRJnTt3tlh/9OhRi7fVpaena8GCBQoJCTEP3cvr2vzDDz+Yy+W9Cvdq8p4AXf7EZ/PmzUpMTLQol/e2ltOnT1/1mNZeoyuIjIxUcnKyxaufL1y4oLlz517Tcdu1a6dWrVppypQpunDhgqRr/82qV6+udu3aafbs2Tp27Fi+c544ceKaYgYAXN+8vb01c+ZMvfjii7r33nsLLffQQw8pOztb48ePz7ft0qVL5vtXQfe9rKwsvf3228UbeAGeeOIJjR07VhMnTiy0jD3uq8uWLbOYp2rLli3avHlzkd4uWFB9GoahqVOn2hxPUX/La3H33XcrOztb06dPt1g/efJkmUwmi2uvUKFCvnO6u7ura9eu+uyzz7Rr1658x7f1d6pQoUKRht4VVO9S/vYyYA/0lAKuI82bN1fv3r01Z84cc3f0LVu26P3331d0dLTat29vUf6GG25QTEyMfv75Z/n7+2vevHlKTU3V/PnzzWU6deqkwMBAxcTEaMSIEXJ3d9e8efPk5+enw4cPXzGee+65R0uXLtUDDzygzp0768CBA5o1a5aaNGmis2fPmsuVK1dOTZo00eLFi3XDDTeoSpUqatq0qZo2bXrN1+gKnnzySU2fPl09evTQkCFDVKNGDX300Ufy8vKSdG3DNUeMGKFu3brpvffeU//+/YvlN5sxY4Zuv/12NWvWTLGxsapXr55SU1OVmJiov//+W7/88ss11wkA4PpV2PCoy7Vt21ZPPvmkEhIStHPnTnXq1Elly5bV3r17tWTJEk2dOlUPPvigWrdurcqVK6t3794aPHiwTCaTPvjggxId5p8nKChIL7744lXLlfZ9tUGDBrr99tv11FNPKTMzU1OmTFHVqlX17LPPXnXfRo0aqX79+nrmmWeUnJwsHx8fffbZZzbPkyQV/be8Fvfee6/at2+v559/XgcPHlTz5s317bff6osvvtDQoUMt5p8KDQ3Vd999p0mTJqlmzZqqW7euwsLC9Nprr2ndunUKCwtTbGysmjRpolOnTmn79u367rvvdOrUKavjCg0N1eLFixUXF6dbb71V3t7eBSZjfXx8dMcdd+iNN97QxYsXVatWLX377bel0uMPuBqSUsB15p133lG9evX03nvv6fPPP1dAQIDi4+PzDb+ScieynDZtmkaMGKGkpCTVrVtXixcvtpizoGzZsvr888/19NNP64UXXlBAQICGDh2qypUrq2/fvleMpU+fPkpJSdHs2bP1zTffqEmTJvrwww+1ZMkSff/99/niHjRokIYNG6asrCyNHTu2wKSUtdfoCry9vbV27VoNGjRIU6dOlbe3t3r16qXWrVura9eu5uSULbp06aL69etrwoQJio2NLZbfrEmTJtq6davGjRun9957T//884+qV6+uFi1aWHQpBwCgJM2aNUuhoaGaPXu2Ro0apTJlyig4OFiPPvqo2rRpI0mqWrWqVqxYoeHDh2v06NGqXLmyHn30UXXs2PGqb9stLaV9X+3Vq5fc3Nw0ZcoUHT9+XK1atdL06dNVo0aNq+5btmxZffnllxo8eLASEhLk5eWlBx54QAMHDlTz5s1tjqkov+W1cHNz0/LlyzVmzBgtXrxY8+fPV3BwsN58803zG57zTJo0Sf369dPo0aN1/vx59e7dW2FhYfL399eWLVv00ksvaenSpXr77bdVtWpV3XTTTXr99ddtiuvpp5/Wzp07NX/+fE2ePFlBQUGF9hBcuHChBg0apBkzZsgwDHXq1Elff/21atasadO5geJiMkojzQ/A6QQHB6tp06ZasWKFvUOBjaZMmaJhw4bp77//5m12AAAAABwOc0oBgAs4f/68xfcLFy5o9uzZatiwIQkpAAAAAA6J4XsA4AK6dOmiwMBAhYSEKC0tTR9++KH27NlT6OuRAQAAAMDeSEoBgAuIjIzUO++8o48++kjZ2dlq0qSJFi1apO7du9s7NAAAAAAoEHNKAQAAAAAAoNQxpxQAAAAAAABKHUkpAAAAF3Pw4EHFxMSobt26KleunOrXr6+xY8cqKyvrivtduHBBAwYMUNWqVeXt7a2uXbsqNTXVoszhw4fVuXNnlS9fXtWrV9eIESN06dKlkrwcAADgophTykY5OTk6evSoKlasKJPJZO9wAABAMTMMQ2fOnFHNmjXl5uZcz/H27NmjnJwczZ49Ww0aNNCuXbsUGxurjIwMTZgwodD9hg0bpq+++kpLliyRr6+vBg4cqC5duuinn36SJGVnZ6tz584KCAjQxo0bdezYMfXq1Utly5bVq6++WuT4aEcBAODaityOMmCTI0eOGJJYWFhYWFhYXHw5cuSIvZsdxeKNN94w6tatW+j206dPG2XLljWWLFliXvfHH38YkozExETDMAxj5cqVhpubm5GSkmIuM3PmTMPHx8fIzMwsciy0o1hYWFhYWK6P5WrtKHpK2ahixYqSpCNHjsjHx8fO0QAAgOKWnp6uOnXqmO/5zi4tLU1VqlQpdPu2bdt08eJFRUREmNc1atRIgYGBSkxM1G233abExEQ1a9ZM/v7+5jKRkZF66qmn9Pvvv6tFixYFHjszM1OZmZnm78b/vWeHdhQAAK6pqO0ouyelZsyYoTfffFMpKSlq3ry5pk2bplatWhVY9vfff9eYMWO0bds2HTp0SJMnT9bQoUOtPuaFCxc0fPhwLVq0SJmZmYqMjNTbb79t0cC6mryu5j4+PjSmAABwYa4wvGzfvn2aNm3aFYfupaSkyMPDQ5UqVbJY7+/vr5SUFHOZ/7aX8r7nlSlIQkKCxo0bl2897SgAAFzb1dpRdp0gYfHixYqLi9PYsWO1fft2NW/eXJGRkTp+/HiB5c+dO6d69erptddeU0BAgM3HHDZsmL788kstWbJE69ev19GjR9WlS5cSuUYAAIDiMnLkSJlMpisue/bssdgnOTlZUVFR6tatm2JjY+0Sd3x8vNLS0szLkSNH7BIHAABwLHbtKTVp0iTFxsaqb9++kqRZs2bpq6++0rx58zRy5Mh85W+99VbdeuutklTg9qIcMy0tTe+++64WLlyoDh06SJLmz5+vxo0ba9OmTbrttttK4lIBAACu2fDhw9WnT58rlqlXr57589GjR9W+fXu1bt1ac+bMueJ+AQEBysrK0unTpy16S6WmppofBgYEBGjLli0W++W9na+wB4aS5OnpKU9PzyueHwAAXH/s1lMqKytL27Zts5i3wM3NTREREUpMTCyxY15tvoTCZGZmKj093WIBAAAoTX5+fmrUqNEVFw8PD0m5PaTatWun0NBQzZ8//6pvEAwNDVXZsmW1Zs0a87qkpCQdPnxY4eHhkqTw8HD99ttvFj3QV69eLR8fHzVp0qQErhgAALgyuyWlTp48qezs7ALnJbjSnATXesyizJdQkISEBPn6+pqXOnXq2BQjAABASctLSAUGBmrChAk6ceKEUlJSLNo6ycnJatSokbnnk6+vr2JiYhQXF6d169Zp27Zt6tu3r8LDw809yTt16qQmTZroscce0y+//KJvvvlGo0eP1oABA+gJBQAArGb3ic6dRXx8vOLi4szf82aSBwAAcDSrV6/Wvn37tG/fPtWuXdtiW96b7y5evKikpCSdO3fOvG3y5Mlyc3NT165dLV4Gk8fd3V0rVqzQU089pfDwcFWoUEG9e/fWSy+9VDoXBgAAXIrdklLVqlWTu7u7eR6CPJfPW1ASxyzKfAkFYS4EAADgLPr06XPVuaeCg4PNCao8Xl5emjFjhmbMmFHofkFBQVq5cmVxhAkAAK5zdhu+5+HhodDQUIt5C3JycrRmzRrzvAUlccyizJcAAAAAAACAkmXX4XtxcXHq3bu3WrZsqVatWmnKlCnKyMgwvzmvV69eqlWrlhISEiTlTmS+e/du8+fk5GTt3LlT3t7eatCgQZGOefl8CVWqVJGPj48GDRpkMV8CAAAAAAAASpZdk1Ldu3fXiRMnNGbMGKWkpCgkJESrVq0yT1R++PBhizfFHD16VC1atDB/nzBhgiZMmKC2bdvq+++/L9IxpavPlwAAAAAAAICSZTL+O5kAiiQ9PV2+vr5KS0uTj4+PvcMBAADFjHt9yaFuAQBwbUW919ttTikU7JdfpOBg6fbb7R0JAACAk/nxx9yG1N132zsSAABQBHYdvof8Ll6UDh2S6L8GAABgpQsXchtSvr72jgQAABQBPaUcTN4UWjk59o0DAADA6ZhMuX95ugcAgFMgKeVg8tpSJKUAAACsRFIKAACnQlLKweT1lKItBQAAYCUaUgAAOBWSUg6G4XsAAAA2oss5AABOhaSUgyEpBQAAYCOG7wEA4FRISjkYHvABAADYiOF7AAA4FZJSDoa2FAAAgI14ugcAgFMhKeVgGL4HAABgI4bvAQDgVEhKORiSUgAAADYiKQUAgFMhKeVg6HUOAABgI+ZBAADAqZCUcjC0pQAAAGzE0z0AAJwKSSkHw/A9AAAAGzF8DwAAp0JSysGQlAIAALARSSkAAJwKSSkHQ69zAAAAGzEPAgAAToWklIOhLQUAAGAjnu4BAOBUSEo5GIbvAQAA2IjhewAAOBWSUg6GpBQAAICNSEoBAOBUSEo5mLy2lER7CgAAwCo83QMAwKmQlHIwbpf9IiSlAAAArEBPKQAAnApJKQdzeVKKh3wAAABWICkFAIBTISnlYEhKAQAA2IjXGAMA4FRISjmYy+eUIikFAABghbyGFI0oAACcAkkpB8OcUgAAADZi+B4AAE6FpJSDYfgeAACAjUhKAQDgVEhKORiG7wEAANiIOaUAAHAqJKUcDD2lAAAAbMScUgAAOBWSUg6GpBQAAICNGL4HAIBTISnlYJjoHAAAwEYkpQAAcCokpRwMc0oBAADYiDmlAABwKiSlHAxJKQAAABsxpxQAAE6FpJSDMZloTwEAANiE4XsAADgVklIOiJ7nAAAANqARBQCAUyEp5YDoKQUAAGADGlEAADgVklIOKO8hH+0pAAAAKzB8DwAAp0JSygGRlAIAALABSSkAAJwKSSkHxHQIAAAANqARBQCAU7F7UmrGjBkKDg6Wl5eXwsLCtGXLliuWX7JkiRo1aiQvLy81a9ZMK1eutNiempqqPn36qGbNmipfvryioqK0d+9eizLt2rWTyWSyWPr371/s12YrpkMAAACwAT2lAABwKnZNSi1evFhxcXEaO3astm/frubNmysyMlLHjx8vsPzGjRvVo0cPxcTEaMeOHYqOjlZ0dLR27dolSTIMQ9HR0frrr7/0xRdfaMeOHQoKClJERIQyMjIsjhUbG6tjx46ZlzfeeKPEr7eoGL4HAABgg7yklERiCgAAJ2DXpNSkSZMUGxurvn37qkmTJpo1a5bKly+vefPmFVh+6tSpioqK0ogRI9S4cWONHz9et9xyi6ZPny5J2rt3rzZt2qSZM2fq1ltv1Y033qiZM2fq/Pnz+vjjjy2OVb58eQUEBJgXHx+fEr/eoiIpBQAAYAOSUgAAOBW7JaWysrK0bds2RURE/P9g3NwUERGhxMTEAvdJTEy0KC9JkZGR5vKZmZmSJC8vL4tjenp6asOGDRb7ffTRR6pWrZqaNm2q+Ph4nTt3rliuqzgwHQIAAIAN3C5r2tKQAgDA4ZWx14lPnjyp7Oxs+fv7W6z39/fXnj17CtwnJSWlwPIpKSmSpEaNGikwMFDx8fGaPXu2KlSooMmTJ+vvv//WsWPHzPs88sgjCgoKUs2aNfXrr7/queeeU1JSkpYuXVpovJmZmeaklySlp6dbfc1FxZxSAAAANri8p1ROjuTubr9YAADAVdktKVUSypYtq6VLlyomJkZVqlSRu7u7IiIidNddd8m47GlZv379zJ+bNWumGjVqqGPHjtq/f7/q169f4LETEhI0bty4Er8GieF7AAAANmH4HgAATsVuw/eqVasmd3d3paamWqxPTU1VQEBAgfsEBARctXxoaKh27typ06dP69ixY1q1apX++ecf1atXr9BYwsLCJEn79u0rtEx8fLzS0tLMy5EjR656jbYiKQUAAGADhu8BAOBU7JaU8vDwUGhoqNasWWNel5OTozVr1ig8PLzAfcLDwy3KS9Lq1asLLO/r6ys/Pz/t3btXW7du1f33319oLDt37pQk1ahRo9Aynp6e8vHxsVhKCnNKAQAA2OC/w/cAAIBDs+vb9+Li4jR37ly9//77+uOPP/TUU08pIyNDffv2lST16tVL8fHx5vJDhgzRqlWrNHHiRO3Zs0cvvviitm7dqoEDB5rLLFmyRN9//73++usvffHFF7rzzjsVHR2tTp06SZL279+v8ePHa9u2bTp48KCWL1+uXr166Y477tDNN99cuhVQCOaUAgAA1+LgwYOKiYlR3bp1Va5cOdWvX19jx45VVlbWFfe7cOGCBgwYoKpVq8rb21tdu3bN10vdZDLlWxYtWlSSl1N0DN8DAMCp2HVOqe7du+vEiRMaM2aMUlJSFBISolWrVpknMz98+LDcLuuG3bp1ay1cuFCjR4/WqFGj1LBhQy1btkxNmzY1lzl27Jji4uKUmpqqGjVqqFevXnrhhRfM2z08PPTdd99pypQpysjIUJ06ddS1a1eNHj269C78Khi+BwAArsWePXuUk5Oj2bNnq0GDBtq1a5diY2OVkZGhCRMmFLrfsGHD9NVXX2nJkiXy9fXVwIED1aVLF/30008W5ebPn6+oqCjz90qVKpXUpViHpBQAAE7FZBjcsW2Rnp4uX19fpaWlFftQvsBA6cgR6eefpZYti/XQAACgiEryXm8Pb775pmbOnKm//vqrwO1paWny8/PTwoUL9eCDD0rKTW41btxYiYmJuu222yTl9pT6/PPPFR0dbXMsJVa3589L5cvnfj5zRvL2Lr5jAwCAIivqvd6uw/dQMOaUAgAAxS0tLU1VqlQpdPu2bdt08eJFRUREmNc1atRIgYGBSkxMtCg7YMAAVatWTa1atdK8efN0tWecmZmZSk9Pt1hKBHNKAQDgVOw6fA8FY04pAABQnPbt26dp06ZdceheSkqKPDw88g3F8/f3V0pKivn7Sy+9pA4dOqh8+fL69ttv9fTTT+vs2bMaPHhwocdOSEjQuHHjrvk6rorhewAAOBV6Sjkg5pQCAAAFGTlyZIETjV++7Nmzx2Kf5ORkRUVFqVu3boqNjb3mGF544QW1adNGLVq00HPPPadnn31Wb7755hX3iY+PV1pamnk5cuTINcdRIJJSAAA4FXpKOSCSUgAAoCDDhw9Xnz59rlimXr165s9Hjx5V+/bt1bp1a82ZM+eK+wUEBCgrK0unT5+26C2VmpqqgICAQvcLCwvT+PHjlZmZKU9PzwLLeHp6FrqtWF32ghySUgAAOD6SUg6IOaUAAEBB/Pz85OfnV6SyycnJat++vUJDQzV//nyLNxoXJDQ0VGXLltWaNWvUtWtXSVJSUpIOHz6s8PDwQvfbuXOnKleuXDpJp6thTikAAJwKSSkHxJxSAADgWiQnJ6tdu3YKCgrShAkTdOLECfO2vF5PycnJ6tixoxYsWKBWrVrJ19dXMTExiouLU5UqVeTj46NBgwYpPDzc/Oa9L7/8Uqmpqbrtttvk5eWl1atX69VXX9Uzzzxjl+vMh+F7AAA4FZJSDojhewAA4FqsXr1a+/bt0759+1S7dm2LbXlvyrt48aKSkpJ07tw587bJkyfLzc1NXbt2VWZmpiIjI/X222+bt5ctW1YzZszQsGHDZBiGGjRooEmTJhXLXFXFgqQUAABOxWRc7R2+KFB6erp8fX2VlpYmHx+fYj1206bS779La9ZIHToU66EBAEARleS9/npXonWbl5hKSZH8/Yv32AAAoEiKeq/n7XsOiJ5SAAAANspLSvHcFQAAh0dSygHRlgIAALARDSkAAJwGSSkHRE8pAAAAG/EaYwAAnAZJKQdEUgoAAMBGvMYYAACnQVLKAZGUAgAAsBHD9wAAcBokpRwQbSkAAAAb0ZACAMBpkJRyQPSUAgAAsBFzSgEA4DRISjkgklIAAAA2Yk4pAACcBkkpB0RSCgAAwEYM3wMAwGmQlHJAtKUAAABsxPA9AACcBkkpB0RPKQAAABsxfA8AAKdBUsoBkZQCAACwEV3OAQBwGiSlHBBJKQAAABuRlAIAwGmQlHJAtKUAAABsxJxSAAA4DZJSDoieUgAAADZiTikAAJwGSSkHRFIKAADARnQ5BwDAaZCUckAkpQAAAGxEUgoAAKdBUsoB0ZYCAACwEXNKAQDgNEhKOSB6SgEAANiIOaUAAHAaJKUcEEkpAAAAG9HlHAAAp0FSygGRlAIAALARw/cAAHAaJKUcEA/4AAAAbMTwPQAAnAZJKQdETykAAAAb8XQPAACnQVLKAZGUAgAAsBFJKQAAnAZJKQdEUgoAAMBGzCkFAIDTICnlgHjABwCA69iwYYO9Q7i+MKcUAABOg6SUA6KnFAAArqNDhw6qW7euRo0apd27d9s7HNfH0z0AAJwGSSkHRFIKAADXcfToUQ0fPlzr169X06ZNFRISojfffFN///23vUNzTSSlAABwGlYlpU6fPq358+fr8ccfV8eOHRUeHq777rtPY8eO1caNG0sqxusOSSkAAFxHtWrVNHDgQP3000/av3+/unXrpvfff1/BwcHq0KGDvcNzPcwpBQCA0yhSUuro0aN64oknVKNGDb388ss6f/68QkJC1LFjR9WuXVvr1q3TnXfeqSZNmmjx4sUlHbPLYyoEAABcU926dTVy5Ei99tpratasmdavX2/vkFwPDSkAAJxGmaIUatGihXr37q1t27apSZMmBZY5f/68li1bpilTpujIkSN65plnijXQ6wkP+AAAcD0//fSTPvroI3366ae6cOGC7r//fiUkJNg7LNfD8D0AAJxGkXpK7d69W2+88UahCSlJKleunHr06KHExET17du3yAHMmDFDwcHB8vLyUlhYmLZs2XLF8kuWLFGjRo3k5eWlZs2aaeXKlRbbU1NT1adPH9WsWVPly5dXVFSU9u7da1HmwoULGjBggKpWrSpvb2917dpVqampRY65pDF8DwAA1xEfH6+6deuqQ4cOOnz4sKZOnaqUlBR98MEHioqKsnd4roenewAAOI0iJaWqVq1q1UGLWn7x4sWKi4vT2LFjtX37djVv3lyRkZE6fvx4geU3btyoHj16KCYmRjt27FB0dLSio6O1a9cuSZJhGIqOjtZff/2lL774Qjt27FBQUJAiIiKUkZFhPs6wYcP05ZdfasmSJVq/fr2OHj2qLl26WHWNJYmkFAAAruOHH37QiBEjlJycrBUrVqhHjx4qX768vcNyXQzfAwDAaVj99j13d3e1b99ep06dslifmpoqd3d3q441adIkxcbGqm/fvmrSpIlmzZql8uXLa968eQWWnzp1qqKiojRixAg1btxY48eP1y233KLp06dLkvbu3atNmzZp5syZuvXWW3XjjTdq5syZOn/+vD7++GNJUlpamt59911NmjRJHTp0UGhoqObPn6+NGzdq06ZN1lZHiaAtBQCA6/jpp5/09NNPq1q1avYO5frA8D0AAJyG1UkpwzCUmZmpli1b6vfff8+3raiysrK0bds2RURE/P9g3NwUERGhxMTEAvdJTEy0KC9JkZGR5vKZmZmSJC8vL4tjenp6asOGDZKkbdu26eLFixbHadSokQIDAws9b96x09PTLZaSQq9zAAAAG5GUAgDAaVidlDKZTPrss8907733Kjw8XF988YXFtqI6efKksrOz5e/vb7He399fKSkpBe6TkpJyxfJ5yaX4+Hj9+++/ysrK0uuvv66///5bx44dMx/Dw8NDlSpVKvJ5JSkhIUG+vr7mpU6dOkW+VmsxfA8AAMBGPN0DAMBp2NRTyt3dXVOnTtWECRPUvXt3vfzyy1b1kiopZcuW1dKlS/Xnn3+qSpUqKl++vNatW6e77rpLbm5WX6qF+Ph4paWlmZcjR44UU9T5kZQCAACwEfMgAADgNMpcy879+vVTw4YN1a1bN/3www9W7VutWjW5u7vne+tdamqqAgICCtwnICDgquVDQ0O1c+dOpaWlKSsrS35+fgoLC1PLli3Nx8jKytLp06ctektd6byS5OnpKU9PT6uu0Va0pQAAAGzE8D0AAJyG1d2HgoKCLCY0b9++vTZt2mR1zyEPDw+FhoZqzZo15nU5OTlas2aNwsPDC9wnPDzcorwkrV69usDyvr6+8vPz0969e7V161bdf//9knKTVmXLlrU4TlJSkg4fPlzoeUsbvc4BAABsRFIKAACnYXVPqQMHDuRb16BBA+3YsSNfL6ariYuLU+/evdWyZUu1atVKU6ZMUUZGhvr27StJ6tWrl2rVqqWEhARJ0pAhQ9S2bVtNnDhRnTt31qJFi7R161bNmTPHfMwlS5bIz89PgYGB+u233zRkyBBFR0erU6dOknKTVTExMYqLi1OVKlXk4+OjQYMGKTw8XLfddpu11VEiGL4HAIBzq1y5cpHn2vzvG41xjXi6BwCA07im4XuX8/LyUlBQkFX7dO/eXSdOnNCYMWOUkpKikJAQrVq1yjyZ+eHDhy3mgmrdurUWLlyo0aNHa9SoUWrYsKGWLVumpk2bmsscO3ZMcXFxSk1NVY0aNdSrVy+98MILFuedPHmy3Nzc1LVrV2VmZioyMlJvv/32NVx98SIpBQCAc5syZYq9Q7h+MQ8CAABOw2QUcYbyoj7xu16e9qWnp8vX11dpaWny8fEp1mOPGCFNmCANH577FwAAlL6SvNdf70q0blu1kn7+WfryS+mee4r32AAAoEiKeq8vck+py5/4GYahp556Si+99JKqV69+TYEiP3qdAwDgmi5cuKCsrCyLdSS8ihlzSgEA4DSKnJTq3bu3xfdBgwapa9euqlevXrEHdb1j+B4AAK4jIyNDzz33nD755BP9888/+bZnZ2fbISoXRkMKAACnYfXb91DyaEsBAOA6nn32Wa1du1YzZ86Up6en3nnnHY0bN041a9bUggUL7B2e66GnFAAATqPYJjpH8WF+TgAAXMeXX36pBQsWqF27durbt6/+97//qUGDBgoKCtJHH32knj172jtE10JSCgAAp0FPKQfEnFIAALiOU6dOmac78PHxMb8U5vbbb9cPP/xgz9BcEw0pAACcRpF7SsXFxVl8z8rK0iuvvCJfX1+L9ZMmTSqeyK5jDN8DAMB11KtXTwcOHFBgYKAaNWqkTz75RK1atdKXX36pSpUq2Ts810OXcwAAnEaRk1I7duyw+N66dWv99ddfFutMeY0AXBOSUgAAuI6+ffvql19+Udu2bTVy5Ejde++9mj59ui5evMjDvJLA8D0AAJxGkZNS69atK8k4cBke8AEA4DqGDRtm/hwREaE9e/Zo27ZtatCggW6++WY7RuaiSEoBAOA0mFPKATEVAgAArisoKEhdunQp0YTUwYMHFRMTo7p166pcuXKqX7++xo4dq6ysrCvuN2fOHLVr104+Pj4ymUw6ffp0vjKnTp1Sz5495ePjo0qVKikmJkZnz54toSuxAQ0pAACcRpF6Sr322msaPHiwypcvf9Wymzdv1smTJ9W5c+drDu56ldeWys62bxwAAMA2b731lvr16ycvLy+99dZbVyw7ePDgYj//nj17lJOTo9mzZ6tBgwbatWuXYmNjlZGRoQkTJhS637lz5xQVFaWoqCjFx8cXWKZnz546duyYVq9erYsXL6pv377q16+fFi5cWOzXYRO6nAMA4DSKlJTavXu3goKC1K1bN917771q2bKl/Pz8JEmXLl3S7t27tWHDBn344Yc6evSoFixYUKJBuzrmlAIAwLlNnjxZPXv2lJeXlyZPnlxoOZPJVCJJqbzEUp569eopKSlJM2fOvGJSaujQoZKk77//vsDtf/zxh1atWqWff/5ZLVu2lCRNmzZNd999tyZMmKCaNWsW2zXYjOF7AAA4jSIlpRYsWKBffvlF06dP1yOPPKL09HS5u7vL09NT586dkyS1aNFCTzzxhPr06SMvL68SDdrVubvn/iUpBQCAczpw4ECBn+0pLS1NVapUuaZjJCYmqlKlSuaElJQ7T5abm5s2b96sBx544FrDvHYkpQAAcBpFnui8efPmmjt3rmbPnq1ff/1Vhw4d0vnz51WtWjWFhISoWrVqJRnndYXhewAAuI6XXnpJzzzzTL5pEM6fP68333xTY8aMKfEY9u3bp2nTpl2xl1RRpKSkqHr16hbrypQpoypVqiglJaXQ/TIzM5WZmWn+np6efk1xXBFdzgEAcBpWT3Tu5uamkJAQ3X///Xr44YcVERFBQqqY0VMKAADXMW7cuAInAj937pzGjRtn1bFGjhwpk8l0xWXPnj0W+yQnJysqKkrdunVTbGzsNV2LrRISEuTr62te6tSpU3Ino6cUAABOo8g9pVB68pJS9JQCAMD5GYYhU16i5DK//PKL1cPphg8frj59+lyxTL169cyfjx49qvbt26t169aaM2eOVecqSEBAgI4fP26x7tKlSzp16pQCAgIK3S8+Pl5xcXHm7+np6SWXmKKnFAAAToOklANi+B4AAM6vcuXK5t5LN9xwg0ViKjs7W2fPnlX//v2tOqafn5/5ZTNXk5ycrPbt2ys0NFTz58+Xm5vVHeTzCQ8P1+nTp7Vt2zaFhoZKktauXaucnByFhYUVup+np6c8PT2v+fxFQlIKAACnQVLKATF8DwAA5zdlyhQZhqHHH39c48aNk6+vr3mbh4eHgoODFR4eXiLnTk5OVrt27RQUFKQJEyboxIkT5m15PZqSk5PVsWNHLViwQK1atZKUO2dUSkqK9u3bJ0n67bffVLFiRQUGBqpKlSpq3LixoqKiFBsbq1mzZunixYsaOHCgHn74Ycd4855EQwoAACdCUsoB0VMKAADn17t3b0lS3bp11bp1a5UtW7bUzr169Wrt27dP+/btU+3atS22Gf8319LFixeVlJRkfpOyJM2aNctinqs77rhDkjR//nzzsMGPPvpIAwcOVMeOHeXm5qauXbvqrbfeKuErsgI9pQAAcBpWJaUuXryocuXKaefOnWratGlJxXTd4wEfAACuo23btsrJydGff/6p48ePK+c/N/i8xE9x6tOnz1XnngoODjYnqPK8+OKLevHFF6+4X5UqVbRw4cJrjLAEkZQCAMBpWJWUKlu2rAIDA5VNF54SRU8pAABcx6ZNm/TII4/o0KFD+ZJAJpOJdlVxIykFAIDTsHrGy+eff16jRo3SqVOnSiIeiLfvAQDgSvr376+WLVtq165dOnXqlP7991/zQnuqBJCUAgDAaVg9p9T06dO1b98+1axZU0FBQapQoYLF9u3btxdbcNcrhu8BAOA69u7dq08//VQNGjSwdyjXB5JSAAA4DauTUtHR0SUQBi7H8D0AAFxHWFiY9u3bR1KqtJCUAgDAaVidlBo7dmxJxIHL0FMKAADXMWjQIA0fPlwpKSlq1qxZvrfw3XzzzXaKzEWRlAIAwGlYnZRCyaOnFAAArqNr166SpMcff9y8zmQyyTAMJjovCTSkAABwGiSlHBA9pQAAcB0HDhywdwjXFxpSAAA4DZJSDoi37wEA4DqCgoLsHcL1heF7AAA4DTd7B4D86HUOAIBr+eCDD9SmTRvVrFlThw4dkiRNmTJFX3zxhZ0jc0EkpQAAcBo2J6WysrKUlJSkS5cuFWc8EL3OAQBwJTNnzlRcXJzuvvtunT592jyHVKVKlTRlyhT7BueKSEoBAOA0rE5KnTt3TjExMSpfvrxuuukmHT58WFLum2Vee+21Yg/wekRPKQAAXMe0adM0d+5cPf/883LPe/IkqWXLlvrtt9/sGJmLIikFAIDTsDopFR8fr19++UXff/+9vLy8zOsjIiK0ePHiYg3uekVPKQAAXMeBAwfUokWLfOs9PT2VkZFhh4hcHEkpAACchtVJqWXLlmn69Om6/fbbZTKZzOtvuukm7d+/v1iDu14x0TkAAK6jbt262rlzZ771q1atUuPGjUs/IFdHUgoAAKdh9dv3Tpw4oerVq+dbn5GRYZGkgu1oSwEA4Dri4uI0YMAAXbhwQYZhaMuWLfr444+VkJCgd955x97huR4aUgAAOA2rk1ItW7bUV199pUGDBkmSORH1zjvvKDw8vHiju07RUwoAANfxxBNPqFy5cho9erTOnTunRx55RDVr1tTUqVP18MMP2zs818PknAAAOA2rk1Kvvvqq7rrrLu3evVuXLl3S1KlTtXv3bm3cuFHr168viRivO7SlAABwLT179lTPnj117tw5nT17tsBe5ygmTM4JAIDTsHpOqdtvv107d+7UpUuX1KxZM3377beqXr26EhMTFRoaWhIxXndoSwEA4Do6dOig06dPS5LKly9vTkilp6erQ4cOdozMRTF8DwAAp2F1TylJql+/vubOnVvcseD/MHwPAADX8f333ysrKyvf+gsXLujHH3+0Q0QujqQUAABOw+qklLu7u44dO5av2/k///yj6tWrK5tMyjWjLQUAgPP79ddfzZ93796tlJQU8/fs7GytWrVKtWrVskdoro2GFAAATsPq4XuGYRS4PjMzUx4eHlYHMGPGDAUHB8vLy0thYWHasmXLFcsvWbJEjRo1kpeXl5o1a6aVK1dabD979qwGDhyo2rVrq1y5cmrSpIlmzZplUaZdu3YymUwWS//+/a2OvaTQUwoAAOcXEhKiFi1ayGQyqUOHDgoJCTEvoaGhevnllzVmzBh7h+l6SEoBAOA0itxT6q233pKU+7a9d955R97e3uZt2dnZ+uGHH9SoUSOrTr548WLFxcVp1qxZCgsL05QpUxQZGamkpKQCJwDduHGjevTooYSEBN1zzz1auHChoqOjtX37djVt2lRS7muX165dqw8//FDBwcH69ttv9fTTT6tmzZq67777zMeKjY3VSy+9ZP5evnx5q2IvSbSlAABwfgcOHJBhGKpXr562bNkiPz8/8zYPDw9Vr15d7nlPolB8aEgBAOA0ipyUmjx5sqTcnlKzZs2yaER5eHgoODg4X4+kq5k0aZJiY2PVt29fSdKsWbP01Vdfad68eRo5cmS+8lOnTlVUVJRGjBghSRo/frxWr16t6dOnm8+9ceNG9e7dW+3atZMk9evXT7Nnz9aWLVssklLly5dXQECAVfGWFnpKAQDg/IKCgiRJOSRHShdJKQAAnEaRh+8dOHBABw4cUNu2bfXLL7+Yvx84cEBJSUn65ptvFBYWVuQTZ2Vladu2bYqIiPj/wbi5KSIiQomJiQXuk5iYaFFekiIjIy3Kt27dWsuXL1dycrIMw9C6dev0559/qlOnThb7ffTRR6pWrZqaNm2q+Ph4nTt37orxZmZmKj093WIpKXltKZJSAAC4hv3792vQoEGKiIhQRESEBg8erP3799s7LNdEUgoAAKdh9UTn69atK5YTnzx5UtnZ2fL397dY7+/vrz179hS4T0pKSoHlL584dNq0aerXr59q166tMmXKyM3NTXPnztUdd9xhLvPII48oKChINWvW1K+//qrnnntOSUlJWrp0aaHxJiQkaNy4cbZcqtXyekrRlgIAwPl98803uu+++xQSEqI2bdpIkn766SfddNNN+vLLL3XnnXfaOUIXw9M9AACchtVJKUn6+++/tXz5ch0+fDjfK44nTZpULIHZatq0adq0aZOWL1+uoKAg/fDDDxowYIBq1qxp7mXVr18/c/lmzZqpRo0a6tixo/bv36/69esXeNz4+HjFxcWZv6enp6tOnTolcg0M3wMAwHWMHDlSw4YN02uvvZZv/XPPPUdSqrjRUwoAAKdhdVJqzZo1uu+++1SvXj3t2bNHTZs21cGDB2UYhm655ZYiH6datWpyd3dXamqqxfrU1NRC53oKCAi4Yvnz589r1KhR+vzzz9W5c2dJ0s0336ydO3dqwoQJ+Yb+5ckbdrhv375Ck1Kenp7y9PQs8vVdC9pSAAC4jj/++EOffPJJvvWPP/64pkyZUvoBuTq6nAMA4DSKPKdUnvj4eD3zzDP67bff5OXlpc8++0xHjhxR27Zt1a1btyIfx8PDQ6GhoVqzZo15XU5OjtasWaPw8PAC9wkPD7coL0mrV682l7948aIuXrwoNzfLy3J3d7/iJKM7d+6UJNWoUaPI8ZckekoBAOA6/Pz8zG2Ny+3cubPAtw3jGvF0DwAAp2F1T6k//vhDH3/8ce7OZcro/Pnz8vb21ksvvaT7779fTz31VJGPFRcXp969e6tly5Zq1aqVpkyZooyMDPPb+Hr16qVatWopISFBkjRkyBC1bdtWEydOVOfOnbVo0SJt3bpVc+bMkST5+Piobdu2GjFihMqVK6egoCCtX79eCxYsMA8r3L9/vxYuXKi7775bVatW1a+//qphw4bpjjvu0M0332xtdZQIpkIAAMB1xMbGql+/fvrrr7/UunVrSblzSr3++usWUwOgmJCUAgDAaVidlKpQoYJ5HqkaNWpo//79uummmyTlTl5uje7du+vEiRMaM2aMUlJSFBISolWrVpknMz98+LBFr6fWrVtr4cKFGj16tEaNGqWGDRtq2bJlatq0qbnMokWLFB8fr549e+rUqVMKCgrSK6+8ov79+0vK7aH13XffmRNgderUUdeuXTV69Ghrq6LE0OscAADX8cILL6hixYqaOHGi4uPjJUk1a9bUiy++qMGDB9s5OhdEUgoAAKdhMgzDsGaH6Ohode7cWbGxsXrmmWf0xRdfqE+fPlq6dKkqV66s7777rqRidSjp6eny9fVVWlqafHx8ivXYJ09Kfn65n7Oz/3/bCgAAlJ6SuNefOXNGklSxYsViOZ6zKsl2lN56SxoyRHr4Yen/evcDAIDSVdR7vdU9pSZNmqSzZ89KksaNG6ezZ89q8eLFatiwod3fvOcqLk9C5eSQlAIAwBUcP35cSUlJkqRGjRrJL+8JFIoXPaUAAHAaViWlsrOz9ffff5vnXqpQoYJmzZpVIoFdz/KG70m5PaXKWJ06BAAAjuLMmTN6+umn9fHHH5tfvOLu7q7u3btrxowZ8vX1tXOELoakFAAATsOqPjju7u7q1KmT/v3335KKB8rfUwoAADivJ554Qps3b9ZXX32l06dP6/Tp01qxYoW2bt2qJ5980t7huR6SUgAAOA2r++A0bdpUf/31l+rWrVsS8UD5e0oBAADntWLFCn3zzTe6/fbbzesiIyM1d+5cRUVF2TEyF8VrjAEAcBpWz1b08ssv65lnntGKFSt07NgxpaenWyy4diSlAABwHVWrVi1wiJ6vr68qV65sh4hcHK8xBgDAaVjdU+ruu++WJN13330ymUzm9YZhyGQyKZssyjVj+B4AAK5j9OjRiouL0wcffKCAgABJUkpKikaMGKEXXnjBztG5IIbvAQDgNKxOSq1bt64k4sBl6CkFAIDrmDlzpvbt26fAwEAFBgZKkg4fPixPT0+dOHFCs2fPNpfdvn27vcJ0HSSlAABwGlYnpdq2bVsSceAyl3VAoz0FAICTi46OtncI1xeSUgAAOA2rk1IoeSZTbnsqJ4eeUgAAOLuxY8faO4TrC0kpAACchtUTnaN0MEcnAACADUhKAQDgNEhKOSjeZgwAAGADklIAADgNklIOKq+nFEkpAAAAK5CUAgDAaZCUclC0pwAAAGxAd3MAAJyG1ROdt2jRQqbLXw/3f0wmk7y8vNSgQQP16dNH7du3L5YAr1f0lAIAALABE3MCAOA0rE5KRUVFaebMmWrWrJlatWolSfr555/166+/qk+fPtq9e7ciIiK0dOlS3X///cUe8PWCnlIAALiG7Oxsvffee1qzZo2OHz+unP/c3NeuXWunyFwUjSgAAJyG1UmpkydPavjw4XrhhRcs1r/88ss6dOiQvv32W40dO1bjx48nKXUN6CkFAIBrGDJkiN577z117txZTZs2LbDHOYoRSSkAAJyG1UmpTz75RNu2bcu3/uGHH1ZoaKjmzp2rHj16aNKkScUS4PWKnucAALiGRYsW6ZNPPtHdd99t71CuDySlAABwGlZPdO7l5aWNGzfmW79x40Z5eXlJknJycsyfYRvm6AQAwDV4eHioQYMG9g7j+kFSCgAAp2F1T6lBgwapf//+2rZtm2699VZJuXNKvfPOOxo1apQk6ZtvvlFISEixBnq9YfgeAACuYfjw4Zo6daqmT5/O0L3SQFIKAACnYXVSavTo0apbt66mT5+uDz74QJJ04403au7cuXrkkUckSf3799dTTz1VvJFeZ2hPAQDgGjZs2KB169bp66+/1k033aSyZctabF+6dKmdInNRNKIAAHAaVg/fk6SePXsqMTFRp06d0qlTp5SYmGhOSElSuXLlGL53jegpBQCAa6hUqZIeeOABtW3bVtWqVZOvr6/FUhIOHjyomJgY1a1bV+XKlVP9+vU1duxYZWVlXXG/OXPmqF27dvLx8ZHJZNLp06fzlQkODpbJZLJYXnvttRK5DpuQlAIAwGlY3VMqT1ZWVoGvNQ4MDLzmoMBE5wAAuIr58+eX+jn37NmjnJwczZ49Ww0aNNCuXbsUGxurjIwMTZgwodD9zp07p6ioKEVFRSk+Pr7Qci+99JJiY2PN3ytWrFis8V8TJuYEAMBpWJ2U2rt3rx5//PF8k50bhiGTyaRsGgDFgvYUAACwVV5iKU+9evWUlJSkmTNnXjEpNXToUEnS999/f8XjV6xYUQEBAcURavGjpxQAAE7D6qRUnz59VKZMGa1YsUI1atRgws4SQk8pAABcx6effqpPPvlEhw8fzjeEbvv27aUSQ1pamqpUqVIsx3rttdc0fvx4BQYG6pFHHtGwYcNUpozNHfCLF40oAACchtWth507d2rbtm1q1KhRScSD/0NPKQAAXMNbb72l559/Xn369NEXX3yhvn37av/+/fr55581YMCAUolh3759mjZt2hV7SRXV4MGDdcstt6hKlSrauHGj4uPjdezYMU2aNKnQfTIzM5WZmWn+np6efs1xFIqeUgAAOA2rJzpv0qSJTp48WRKx4DJ5DxtJSgEA4NzefvttzZkzR9OmTZOHh4eeffZZrV69WoMHD1ZaWppVxxo5cmS+Scb/u+zZs8din+TkZEVFRalbt24W80DZKi4uTu3atdPNN9+s/v37a+LEiZo2bZpF0um/EhISLCZ3r1OnzjXHUSiSUgAAOA2rk1Kvv/66nn32WX3//ff6559/lJ6ebrGgeOQlpS5dsm8cAADg2hw+fFitW7eWlPuG4jNnzkiSHnvsMX388cdWHWv48OH6448/rrjUq1fPXP7o0aNq3769WrdurTlz5hTfRV0mLCxMly5d0sGDBwstEx8fr7S0NPNy5MiREolFEkkpAACciNXD9yIiIiRJHTt2tFjPROfFi6QUAACuISAgQKdOnVJQUJACAwO1adMmNW/eXAcOHJBhGFYdy8/PT35+fkUqm5ycrPbt2ys0NFTz58+Xm5vVzyKLZOfOnXJzc1P16tULLePp6SlPT88SOX8+JKUAAHAaViel1q1bVxJx4D9ISgEA4Bo6dOig5cuXq0WLFurbt6+GDRumTz/9VFu3blWXLl1K5JzJyclq166dgoKCNGHCBJ04ccK8Le+tecnJyerYsaMWLFigVq1aSZJSUlKUkpKiffv2SZJ+++03VaxYUYGBgapSpYoSExO1efNmtW/fXhUrVlRiYqKGDRumRx99VJUrVy6Ra7EaSSkAAJyG1Umptm3blkQc+A+SUgAAuIY5c+Yo5/8SJAMGDFDVqlW1ceNG3XfffXryySdL5JyrV6/Wvn37tG/fPtWuXdtiW17vrIsXLyopKUnnzp0zb5s1a5bGjRtn/n7HHXdIkubPn68+ffrI09NTixYt0osvvqjMzEzVrVtXw4YNU1xcXIlch01ISgEA4DRMRhH6jf/6669q2rSp3Nzc9Ouvv16x7M0331xswTmy9PR0+fr6Ki0tTT4+PsV+/IgIac0aaeFCqUePYj88AAC4ipK+11/PSrRuf/tNuvlmyd9fSkkp3mMDAIAiKeq9vkg9pUJCQpSSkqLq1asrJCREJpOpwDkQmFOq+NBTCgAA1/Hjjz9q9uzZ2r9/vz799FPVqlVLH3zwgerWravbb7/d3uG5lryeUrRJAQBweEVKSh04cMA8qeaBAwdKNCDkIikFAIBr+Oyzz/TYY4+pZ8+e2rFjhzIzMyVJaWlpevXVV7Vy5Uo7R+hi3N1z/5KUAgDA4RXpNSxBQUEymUy6ePGixo0bp5ycHAUFBRW4oHjktadISgEA4NxefvllzZo1S3PnzlXZsmXN69u0aaPt27fbMTIXRVIKAACnYdW7gcuWLavPPvuspGLBZegpBQCAa0hKSjJPGH45X19fnT59uvQDcnU0ogAAcBpWJaUkKTo6WsuWLSuBUHC5vPYUD/kAAHBuAQEB2rdvX771GzZsUL169ewQkYujpxQAAE6jSHNKXa5hw4Z66aWX9NNPPyk0NFQVKlSw2D548OBiC+56xkM+AABcQ2xsrIYMGaJ58+bJZDLp6NGjSkxM1DPPPKMXXnjB3uG5Hp7sAQDgNKxOSr377ruqVKmStm3bpm3btllsM5lMJKWKCUkpAABcw8iRI5WTk6OOHTvq3LlzuuOOO+Tp6alnnnlGgwYNsnd4roeJOQEAcBpWD987cOBAoctff/1ldQAzZsxQcHCwvLy8FBYWpi1btlyx/JIlS9SoUSN5eXmpWbNm+d5Yc/bsWQ0cOFC1a9dWuXLl1KRJE82aNcuizIULFzRgwABVrVpV3t7e6tq1q1JTU62OvSSRlAIAwDWYTCY9//zzOnXqlHbt2qVNmzbpxIkTGj9+vL1Dc01lLnvmmpNjvzgAAMBVWZ2UKk6LFy9WXFycxo4dq+3bt6t58+aKjIzU8ePHCyy/ceNG9ejRQzExMdqxY4eio6MVHR2tXbt2mcvExcVp1apV+vDDD/XHH39o6NChGjhwoJYvX24uM2zYMH355ZdasmSJ1q9fr6NHj6pLly4lfr3WICkFAIBr8fDwUJMmTdSqVSt5e3vbOxzXlddTSqIhBQCAgzMZhmFcrVBcXJzGjx+vChUqKC4u7oplJ02aVOSTh4WF6dZbb9X06dMlSTk5OapTp44GDRqkkSNH5ivfvXt3ZWRkaMWKFeZ1t912m0JCQsy9oZo2baru3btbzNEQGhqqu+66Sy+//LLS0tLk5+enhQsX6sEHH5Qk7dmzR40bN1ZiYqJuu+22IsWenp4uX19fpaWlycfHp8jXXFQDBkhvvy2NHSu9+GKxHx4AAFzFtd7rH3/88SKVmzdvntXHdnYl2o46e1aqWDH387lzUrlyxXt8AABwVUW91xdpTqkdO3bo4sWL5s+FMZlMRQ4wKytL27ZtU3x8vHmdm5ubIiIilJiYWOA+iYmJ+ZJikZGRFm8DbN26tZYvX67HH39cNWvW1Pfff68///xTkydPliRt27ZNFy9eVEREhHmfRo0aKTAw0KqkVEmjpxQAAM7tvffeU1BQkFq0aKEiPANEcaGnFAAATqNISal169bpr7/+kq+vr9atW1csJz558qSys7Pl7+9vsd7f31979uwpcJ+UlJQCy6ekpJi/T5s2Tf369VPt2rVVpkwZubm5ae7cubrjjjvMx/Dw8FClSpWueJz/yszMVGZmpvl7enp6ka7TViSlAABwbk899ZQ+/vhjHThwQH379tWjjz6qKlWq2Dss13f5nFK8gQ8AAIdW5DmlGjZsqBMnTpi/d+/e3eEmB5dyk1KbNm3S8uXLtW3bNk2cOFEDBgzQd999d03HTUhIkK+vr3mpU6dOMUVcMJJSAAA4txkzZujYsWN69tln9eWXX6pOnTp66KGH9M0339BzqiTRUwoAAKdR5KTUfxtPK1euVEZGhs0nrlatmtzd3fMltlJTUxUQEFDgPgEBAVcsf/78eY0aNUqTJk3Svffeq5tvvlkDBw5U9+7dNWHCBPMxsrKydPr06SKfV5Li4+OVlpZmXo4cOWLtJVuFtxkDAOD8PD091aNHD61evVq7d+/WTTfdpKefflrBwcE6e/asvcNzTW6XNW/pKQUAgEOz29v3PDw8FBoaqjVr1pjX5eTkaM2aNQoPDy9wn/DwcIvykrR69Wpz+YsXL+rixYtyc7O8LHd3d+X83yuBQ0NDVbZsWYvjJCUl6fDhw4WeV8ptVPr4+FgsJYmeUgAAuBY3NzeZTCYZhqFskiUli4YUAABOoUhzSkm5k5j/dyJzayY2L0hcXJx69+6tli1bqlWrVpoyZYoyMjLUt29fSVKvXr1Uq1YtJSQkSJKGDBmitm3bauLEiercubMWLVqkrVu3as6cOZIkHx8ftW3bViNGjFC5cuUUFBSk9evXa8GCBea3Avr6+iomJkZxcXGqUqWKfHx8NGjQIIWHhzvMJOcSbSkAAFxBZmamli5dqnnz5mnDhg265557NH36dEVFReV7iIZi5O6e24gi+QcAgEMrclLKMAz16dNHnp6ekqQLFy6of//+qlChgkW5pUuXFvnk3bt314kTJzRmzBilpKQoJCREq1atMk9mfvjwYYsGW+vWrbVw4UKNHj1ao0aNUsOGDbVs2TI1bdrUXGbRokWKj49Xz549derUKQUFBemVV15R//79zWUmT54sNzc3de3aVZmZmYqMjNTbb79d5LhLQ15SirYUAADO6emnn9aiRYtUp04dPf744/r4449VrVo1e4d1fShTRsrM5OkeAAAOzmQUcabNvN5LVzN//vxrCshZpKeny9fXV2lpaSUylO+NN6TnnpP69JGukyoFAMChXOu93s3NTYGBgWrRosUVe5db80DPVZR0O0q+vlJ6uvTnn1LDhsV/fAAAcEVFvdcXuafU9ZJschQM3wMAwLn16tXrmqc6gI3ocg4AgFMoclIKpYukFAAAzu29996zdwjXL15jDACAU2CGTQdFUgoAAMBG9JQCAMApkJRyUCSlAAAAbERPKQAAnAJJKQdFUgoAAMBG9JQCAMApkJRyUCSlAAAAbERPKQAAnAJJKQdFUgoAAMBG9JQCAMApkJRyUDzgAwAAsBENKQAAnAJJKQdFTykAAAAb5SWl6CkFAIBDIynloOh1DgAAYCOe7gEA4BRISjko2lIAAAA2oqcUAABOgaSUgyIpBQAAYCO6nAMA4BRISjkoklIAAAA2YqJzAACcAkkpB0VSCgAAwEb0lAIAwCmQlHJQeW2pixftGwcAAIDToacUAABOgaSUgypbNvcvbSkAAAAr0VMKAACnQFLKQXl45P7NyrJvHAAAAE6HnlIAADgFklIOKi8plZlp3zgAAACcDj2lAABwCiSlHBQ9pQAAAGxETykAAJwCSSkHRVIKAADARnlJKXpKAQDg0EhKOShPz9y/WVmSYdg3FgAAAKeSN3yPnlIAADg0klIOKq+nlGHwkA8AAMAqea8xvnjRvnEAAIArIinloPKSUhJD+AAAAKzCPAgAADgFklIOiqQUAACAjfIaUvSUAgDAoZGUclB5UyFIJKUAAACskjd8j0YUAAAOjaSUgzKZ6HkOAABgE3pKAQDgFEhKObC89lRmpn3jAAAAcCr0lAIAwCmQlHJg9JQCAACwAY0oAACcAkkpB0Z7CgAAwAZ5PaUYvgcAgEMjKeXASEoBAADYgEYUAABOgaSUA/P0zP1LewoAAMAKTHQOAIBTICnlwHjIBwAAYAMmOgcAwCmQlHJgJKUAAABsQE8pAACcAkkpB0ZSCgAA2OLgwYOKiYlR3bp1Va5cOdWvX19jx45V1hUaFadOndKgQYN04403qly5cgoMDNTgwYOVlpZmUe7w4cPq3Lmzypcvr+rVq2vEiBG6dOlSSV+SdegpBQCAUyhj7wBQOJJSAADAFnv27FFOTo5mz56tBg0aaNeuXYqNjVVGRoYmTJhQ4D5Hjx7V0aNHNWHCBDVp0kSHDh1S//79dfToUX366aeSpOzsbHXu3FkBAQHauHGjjh07pl69eqls2bJ69dVXS/MSr4xGFAAAToGklAPLa09lZto3DgAA4FyioqIUFRVl/l6vXj0lJSVp5syZhSalmjZtqs8++8z8vX79+nrllVf06KOP6tKlSypTpoy+/fZb7d69W9999538/f0VEhKi8ePH67nnntOLL74oj7zGi70xfA8AAKfA8D0HxkM+AABQXNLS0lSlShWr9/Hx8VGZMrnPMRMTE9WsWTP5+/uby0RGRio9PV2///57ocfJzMxUenq6xVKiGL4HAIBTICnlwEhKAQCA4rBv3z5NmzZNTz75ZJH3OXnypMaPH69+/fqZ16WkpFgkpCSZv6ekpBR6rISEBPn6+pqXOnXqWHkFVqKnFAAAToGklAMjKQUAAC43cuRImUymKy579uyx2Cc5OVlRUVHq1q2bYmNji3Se9PR0de7cWU2aNNGLL754zXHHx8crLS3NvBw5cuSaj3lF9JQCAMApOERSasaMGQoODpaXl5fCwsK0ZcuWK5ZfsmSJGjVqJC8vLzVr1kwrV6602F5YI+3NN980lwkODs63/bXXXiuR67OVp2fuX9pTAABAkoYPH64//vjjiku9evXM5Y8ePar27durdevWmjNnTpHOcebMGUVFRalixYr6/PPPVTYvwSMpICBAqampFuXzvgcEBBR6TE9PT/n4+FgsJYonewAAOAW7T3S+ePFixcXFadasWQoLC9OUKVMUGRmppKQkVa9ePV/5jRs3qkePHkpISNA999yjhQsXKjo6Wtu3b1fTpk0lSceOHbPY5+uvv1ZMTIy6du1qsf6ll16yeGJYsWLFErhC2+UlpS5csG8cAADAMfj5+cnPz69IZZOTk9W+fXuFhoZq/vz5cnO7+rPI9PR0RUZGytPTU8uXL5eXl5fF9vDwcL3yyis6fvy4uZ22evVq+fj4qEmTJtZfUEnJS6QxfA8AAIdm955SkyZNUmxsrPr27asmTZpo1qxZKl++vObNm1dg+alTpyoqKkojRoxQ48aNNX78eN1yyy2aPn26uUxAQIDF8sUXX6h9+/YWTw6l3CTU5eUqVKhQotdqrXLlcv+eP2/fOAAAgHNJTk5Wu3btFBgYqAkTJujEiRNKSUmxmPcpOTlZjRo1MvdQT09PV6dOnZSRkaF3331X6enp5n2ys7MlSZ06dVKTJk302GOP6ZdfftE333yj0aNHa8CAAfLMe5rmCOgpBQCAU7BrUiorK0vbtm1TRESEeZ2bm5siIiKUmJhY4D6JiYkW5aXct74UVj41NVVfffWVYmJi8m177bXXVLVqVbVo0UJvvvmmLl26VGispf7WGEnly+f+PXeuxE8FAABcyOrVq7Vv3z6tWbNGtWvXVo0aNcxLnosXLyopKUnn/q+hsX37dm3evFm//fabGjRoYLFP3hxQ7u7uWrFihdzd3RUeHq5HH31UvXr10ksvvWSX6ywUE50DAOAU7Dp87+TJk8rOzi7wLS7/naQzT2FvfSnsjS/vv/++KlasqC5dulisHzx4sG655RZVqVJFGzduVHx8vI4dO6ZJkyYVeJyEhASNGzeuqJdWLOgpBQAAbNGnTx/16dPnimWCg4NlGIb5e7t27Sy+FyYoKCjffJ4OJ2/4XmamfeMAAABXZPc5pUravHnz1LNnz3xzIsTFxZk/33zzzfLw8NCTTz6phISEArufx8fHW+yTnp5e4q8zpqcUAACADfLafSSlAABwaHZNSlWrVk3u7u4FvsWlsDe4FPbWl4LK//jjj0pKStLixYuvGktYWJguXbqkgwcP6sYbb8y33dPTs9TnSqCnFAAAgA0ub0QZhmQy2TceAABQILvOKeXh4aHQ0FCtWbPGvC4nJ0dr1qxReHh4gfuEh4dblJdy500oqPy7776r0NBQNW/e/Kqx7Ny5U25ubgW+8c9e6CkFAABgg7yklERvKQAAHJjdh+/FxcWpd+/eatmypVq1aqUpU6YoIyNDffv2lST16tVLtWrVUkJCgiRpyJAhatu2rSZOnKjOnTtr0aJF2rp1q+bMmWNx3PT0dC1ZskQTJ07Md87ExERt3rxZ7du3V8WKFZWYmKhhw4bp0UcfVeXKlUv+oouInlIAAAA2uDwpdf78/x/OBwAAHIrdk1Ldu3fXiRMnNGbMGKWkpCgkJESrVq0yT2Z++PBhubn9/w5drVu31sKFCzV69GiNGjVKDRs21LJly9S0aVOL4y5atEiGYahHjx75zunp6alFixbpxRdfVGZmpurWrathw4ZZzBnlCOgpBQAAYIOyZSV3dyk7Ozcp5UAPHQEAwP9nMorymhXkk56eLl9fX6WlpcnHx6dEzvHtt1JkpBQSIu3YUSKnAAAAhSiNe/31qlTqtmJF6exZad8+qX79kjkHAAAoUFHv9XadUwpXltfznJ5SAAAAVmIeBAAAHB5JKQeWN3yPthQAAICVSEoBAODwSEo5MHpKAQAA2IikFAAADo+klAOjpxQAAICNSEoBAODwSEo5sMt7SjEdPQAAgBVISgEA4PBISjkwb+///5khfAAAAFYgKQUAgMMjKeXAypeX3N1zP6el2TcWAAAAp0JSCgAAh0dSyoGZTJKvb+5nklIAAABWYHJOAAAcHkkpB+fjk/uXpBQAAIAVKlTI/Xv2rH3jAAAAhSIp5eDoKQUAAGCDvCd76en2jQMAABSKpJSDIykFAABgAxpRAAA4PJJSDo72FAAAgA2YAwEAAIdHUsrB5SWl6HkOAABgBZ7sAQDg8EhKOTge8gEAANiApBQAAA6PpJSDq1Qp9+/p0/aMAgAAwMnQ3RwAAIdHUsrB+fnl/j1+3L5xAAAAOBW6mwMA4PBISjk4f//cvykp9o0DAADAqTB8DwAAh0dSysHlJaVSU+0bBwAAgFOpWjX377//Spcu2TcWAABQIJJSDo6kFAAAgA2qVZPc3CTDkE6csHc0AACgACSlHFxAQO7ff/+VsrLsGwsAAIDTcHeXqlfP/XzsmH1jAQAABSIp5eAqV5bKlMn9zLxSAAAAVqhRI/cvjSgAABwSSSkH5+YmBQfnfv7rL7uGAgAA4FzyupyTlAIAwCGRlHICDRvm/t27175xAAAAOJW8nlJHjtg3DgAAUCCSUk6ApBQAAIANbrgh9++ff9o3DgAAUCCSUk4grz21e7d94wAAAHAqN96Y+zcpyb5xAACAApGUcgK33pr7NzFRysmxbywAAABOo1Gj3L979kiXLtk3FgAAkA9JKSfQooVUrpx06pT0xx/2jgYAAMBJNGwoVaokZWRIO3faOxoAAPAfJKWcQNmyUrt2uZ8XL7ZrKAAAAM7D3V26447czytX2jcWAACQD0kpJ/HYY7l/335bOn7cvrEAAAA4jQcfzP07Z4509qx9YwEAABZMhmEY9g7CGaWnp8vX11dpaWny8fEp8fNdupQ7jG/Xrtye6IMH506AXqWKVKaM5OaW+zDQzU0ymUo8HAAAHMoNN+TeA4tTad/rryelWrcXLuT+AzlyRLrtNunpp6W6daUKFSQPj/wNp4IaUjSuAACuLG8OxmJU1Ht9mWI/M0pEmTLSkiVSx47S3r3SoEH2jggAAMdx5ozk7W3vKOCQvLykhQulu++WNm3KXQAAQC53d7u+DISklBNp1Ej6/ffc3ufr1kl//y39+6+UnZ37Vr68vwAAALjM7bdLv/6aOw/Cpk1SSkru5OdZWf+/zOWDBwr7DACAq3F3t+vpGb5nI7r0AwDg2rjXlxzqFgAA11bUez0TnQMAAAAAAKDUkZQCAAAAAABAqSMpBQAAAAAAgFJHUgoAAAAAAACljqQUAAAAAAAASp1DJKVmzJih4OBgeXl5KSwsTFu2bLli+SVLlqhRo0by8vJSs2bNtHLlSovtJpOpwOXNN980lzl16pR69uwpHx8fVapUSTExMTp79myJXB8AAAAAAAAs2T0ptXjxYsXFxWns2LHavn27mjdvrsjISB0/frzA8hs3blSPHj0UExOjHTt2KDo6WtHR0dq1a5e5zLFjxyyWefPmyWQyqWvXruYyPXv21O+//67Vq1drxYoV+uGHH9SvX78Sv14AAAAAAABIJsMwDHsGEBYWpltvvVXTp0+XJOXk5KhOnToaNGiQRo4cma989+7dlZGRoRUrVpjX3XbbbQoJCdGsWbMKPEd0dLTOnDmjNWvWSJL++OMPNWnSRD///LNatmwpSVq1apXuvvtu/f3336pZs+ZV405PT5evr6/S0tLk4+Nj9XUDAADHxr2+5FC3AAC4tqLe6+3aUyorK0vbtm1TRESEeZ2bm5siIiKUmJhY4D6JiYkW5SUpMjKy0PKpqan66quvFBMTY3GMSpUqmRNSkhQRESE3Nzdt3rz5Wi4JAAAAAAAARVDGnic/efKksrOz5e/vb7He399fe/bsKXCflJSUAsunpKQUWP79999XxYoV1aVLF4tjVK9e3aJcmTJlVKVKlUKPk5mZqczMTPP39PT0wi8MAAAAAAAAV2T3OaVK2rx589SzZ095eXld03ESEhLk6+trXurUqVNMEQIAAAAAAFx/7NpTqlq1anJ3d1dqaqrF+tTUVAUEBBS4T0BAQJHL//jjj0pKStLixYvzHeO/E6lfunRJp06dKvS88fHxiouLM39PS0tTYGAgPaYAAHBRefd4O0+/6ZLy6pR2FAAArqmo7Si7JqU8PDwUGhqqNWvWKDo6WlLuROdr1qzRwIEDC9wnPDxca9as0dChQ83rVq9erfDw8Hxl3333XYWGhqp58+b5jnH69Glt27ZNoaGhkqS1a9cqJydHYWFhBZ7X09NTnp6e5u95FUyPKQAAXNuZM2fk6+tr7zBcypkzZyTRjgIAwNVdrR1l97fvLV68WL1799bs2bPVqlUrTZkyRZ988on27Nkjf39/9erVS7Vq1VJCQoIkaePGjWrbtq1ee+01de7cWYsWLdKrr76q7du3q2nTpubjpqenq0aNGpo4caL69++f77x33XWXUlNTNWvWLF28eFF9+/ZVy5YttXDhwiLFnZOTo6NHj6pixYoymUzFUxmXxV6nTh0dOXKEN9KUAuq7dFHfpY86L13Ud+kqyfo2DENnzpxRzZo15ebm8jMelCraUa6D+i5d1Hfpo85LF/VduhyhHWXXnlKS1L17d504cUJjxoxRSkqKQkJCtGrVKvNk5ocPH7a4gNatW2vhwoUaPXq0Ro0apYYNG2rZsmUWCSlJWrRokQzDUI8ePQo870cffaSBAweqY8eOcnNzU9euXfXWW28VOW43NzfVrl3bhisuOh8fH/5DLEXUd+mivksfdV66qO/SVVL1TQ+pkkE7yvVQ36WL+i591Hnpor5Llz3bUXbvKYX80tPT5evrq7S0NP5DLAXUd+mivksfdV66qO/SRX3jv/g3Ubqo79JFfZc+6rx0Ud+lyxHqm77oAAAAAAAAKHUkpRyQp6enxo4dazGxOkoO9V26qO/SR52XLuq7dFHf+C/+TZQu6rt0Ud+ljzovXdR36XKE+mb4HgAAAAAAAEodPaUAAAAAAABQ6khKAQAAAAAAoNSRlAIAAAAAAECpIynlYGbMmKHg4GB5eXkpLCxMW7ZssXdILiEhIUG33nqrKlasqOrVqys6OlpJSUkWZS5cuKABAwaoatWq8vb2VteuXZWammqniF3La6+9JpPJpKFDh5rXUd/FLzk5WY8++qiqVq2qcuXKqVmzZtq6dat5u2EYGjNmjGrUqKFy5copIiJCe/futWPEzis7O1svvPCC6tatq3Llyql+/foaP368Lp+mkfq23Q8//KB7771XNWvWlMlk0rJlyyy2F6VuT506pZ49e8rHx0eVKlVSTEyMzp49W4pXAXugHVUyaEfZF+2o0kE7qvTQjipZztaOIinlQBYvXqy4uDiNHTtW27dvV/PmzRUZGanjx4/bOzSnt379eg0YMECbNm3S6tWrdfHiRXXq1EkZGRnmMsOGDdOXX36pJUuWaP369Tp69Ki6dOlix6hdw88//6zZs2fr5ptvtlhPfRevf//9V23atFHZsmX19ddfa/fu3Zo4caIqV65sLvPGG2/orbfe0qxZs7R582ZVqFBBkZGRunDhgh0jd06vv/66Zs6cqenTp+uPP/7Q66+/rjfeeEPTpk0zl6G+bZeRkaHmzZtrxowZBW4vSt327NlTv//+u1avXq0VK1bohx9+UL9+/UrrEmAHtKNKDu0o+6EdVTpoR5Uu2lEly+naUQYcRqtWrYwBAwaYv2dnZxs1a9Y0EhIS7BiVazp+/LghyVi/fr1hGIZx+vRpo2zZssaSJUvMZf744w9DkpGYmGivMJ3emTNnjIYNGxqrV6822rZtawwZMsQwDOq7JDz33HPG7bffXuj2nJwcIyAgwHjzzTfN606fPm14enoaH3/8cWmE6FI6d+5sPP744xbrunTpYvTs2dMwDOq7OEkyPv/8c/P3otTt7t27DUnGzz//bC7z9ddfGyaTyUhOTi612FG6aEeVHtpRpYN2VOmhHVW6aEeVHmdoR9FTykFkZWVp27ZtioiIMK9zc3NTRESEEhMT7RiZa0pLS5MkValSRZK0bds2Xbx40aL+GzVqpMDAQOr/GgwYMECdO3e2qFeJ+i4Jy5cvV8uWLdWtWzdVr15dLVq00Ny5c83bDxw4oJSUFIs69/X1VVhYGHVug9atW2vNmjX6888/JUm//PKLNmzYoLvuuksS9V2SilK3iYmJqlSpklq2bGkuExERITc3N23evLnUY0bJox1VumhHlQ7aUaWHdlTpoh1lP47YjipT7EeETU6ePKns7Gz5+/tbrPf399eePXvsFJVrysnJ0dChQ9WmTRs1bdpUkpSSkiIPDw9VqlTJoqy/v79SUlLsEKXzW7RokbZv366ff/453zbqu/j99ddfmjlzpuLi4jRq1Cj9/PPPGjx4sDw8PNS7d29zvRb0/xjq3HojR45Uenq6GjVqJHd3d2VnZ+uVV15Rz549JYn6LkFFqduUlBRVr17dYnuZMmVUpUoV6t9F0Y4qPbSjSgftqNJFO6p00Y6yH0dsR5GUwnVnwIAB2rVrlzZs2GDvUFzWkSNHNGTIEK1evVpeXl72Due6kJOTo5YtW+rVV1+VJLVo0UK7du3SrFmz1Lt3bztH53o++eQTffTRR1q4cKFuuukm7dy5U0OHDlXNmjWpbwAujXZUyaMdVfpoR5Uu2lG4HMP3HES1atXk7u6e760ZqampCggIsFNUrmfgwIFasWKF1q1bp9q1a5vXBwQEKCsrS6dPn7YoT/3bZtu2bTp+/LhuueUWlSlTRmXKlNH69ev11ltvqUyZMvL396e+i1mNGjXUpEkTi3WNGzfW4cOHJclcr/w/pniMGDFCI0eO1MMPP6xmzZrpscce07Bhw5SQkCCJ+i5JRanbgICAfJNbX7p0SadOnaL+XRTtqNJBO6p00I4qfbSjShftKPtxxHYUSSkH4eHhodDQUK1Zs8a8LicnR2vWrFF4eLgdI3MNhmFo4MCB+vzzz7V27VrVrVvXYntoaKjKli1rUf9JSUk6fPgw9W+Djh076rffftPOnTvNS8uWLdWzZ0/zZ+q7eLVp0ybf67n//PNPBQUFSZLq1q2rgIAAizpPT0/X5s2bqXMbnDt3Tm5ulrdQd3d35eTkSKK+S1JR6jY8PFynT5/Wtm3bzGXWrl2rnJwchYWFlXrMKHm0o0oW7ajSRTuq9NGOKl20o+zHIdtRxT51Omy2aNEiw9PT03jvvfeM3bt3G/369TMqVapkpKSk2Ds0p/fUU08Zvr6+xvfff28cO3bMvJw7d85cpn///kZgYKCxdu1aY+vWrUZ4eLgRHh5ux6hdy+VvjTEM6ru4bdmyxShTpozxyiuvGHv37jU++ugjo3z58saHH35oLvPaa68ZlSpVMr744gvj119/Ne6//36jbt26xvnz5+0YuXPq3bu3UatWLWPFihXGgQMHjKVLlxrVqlUznn32WXMZ6tt2Z86cMXbs2GHs2LHDkGRMmjTJ2LFjh3Ho0CHDMIpWt1FRUUaLFi2MzZs3Gxs2bDAaNmxo9OjRw16XhFJAO6rk0I6yP9pRJYt2VOmiHVWynK0dRVLKwUybNs0IDAw0PDw8jFatWhmbNm2yd0guQVKBy/z5881lzp8/bzz99NNG5cqVjfLlyxsPPPCAcezYMfsF7WL+25iivovfl19+aTRt2tTw9PQ0GjVqZMyZM8die05OjvHCCy8Y/v7+hqenp9GxY0cjKSnJTtE6t/T0dGPIkCFGYGCg4eXlZdSrV894/vnnjczMTHMZ6tt269atK/D/2b179zYMo2h1+88//xg9evQwvL29DR8fH6Nv377GmTNn7HA1KE20o0oG7Sj7ox1V8mhHlR7aUSXL2dpRJsMwjOLvfwUAAAAAAAAUjjmlAAAAAAAAUOpISgEAAAAAAKDUkZQCAAAAAABAqSMpBQAAAAAAgFJHUgoAAAAAAACljqQUAAAAAAAASh1JKQAAAAAAAJQ6klIAAAAAAAAodSSlADi9Pn36KDo62m7nf+yxx/Tqq6+W2PF3796t2rVrKyMjo8TOAQAArk+0owDYk8kwDMPeQQBAYUwm0xW3jx07VsOGDZNhGKpUqVLpBHWZX375RR06dNChQ4fk7e1dYud58MEH1bx5c73wwgsldg4AAOBaaEfloh0FOC6SUgAcWkpKivnz4sWLNWbMGCUlJZnXeXt7l2gj5mqeeOIJlSlTRrNmzSrR83z11VeKjY3V4cOHVaZMmRI9FwAAcA20o3LRjgIcF8P3ADi0gIAA8+Lr6yuTyWSxztvbO1+383bt2mnQoEEaOnSoKleuLH9/f82dO1cZGRnq27evKlasqAYNGujrr7+2ONeuXbt01113ydvbW/7+/nrsscd08uTJQmPLzs7Wp59+qnvvvddifXBwsF5++WX16tVL3t7eCgoK0vLly3XixAndf//98vb21s0336ytW7ea9zl06JDuvfdeVa5cWRUqVNBNN92klStXmrffeeedOnXqlNavX3+NNQoAAK4XtKNy0Y4CHBdJKQAu6f3331e1atW0ZcsWDRo0SE899ZS6deum1q1ba/v27erUqZMee+wxnTt3TpJ0+vRpdejQQS1atNDWrVu1atUqpaam6qGHHir0HL/++qvS0tLUsmXLfNsmT56sNm3aaMeOHercubMee+wx9erVS48++qi2b9+u+vXrq1evXsrrrDpgwABlZmbqhx9+0G+//abXX3/d4smlh4eHQkJC9OOPPxZzTQEAAFiiHQWgtJCUAuCSmjdvrtGjR6thw4aKj4+Xl5eXqlWrptjYWDVs2FBjxozRP//8o19//VWSNH36dLVo0UKvvvqqGjVqpBYtWmjevHlat26d/vzzzwLPcejQIbm7u6t69er5tt1999168sknzedKT0/Xrbfeqm7duumGG27Qc889pz/++EOpqamSpMOHD6tNmzZq1qyZ6tWrp3vuuUd33HGHxTFr1qypQ4cOFXNNAQAAWKIdBaC0kJQC4JJuvvlm82d3d3dVrVpVzZo1M6/z9/eXJB0/flxS7kSb69atM8+t4O3trUaNGkmS9u/fX+A5zp8/L09PzwInEb38/HnnutL5Bw8erJdffllt2rTR2LFjzY28y5UrV878RBIAAKCk0I4CUFpISgFwSWXLlrX4bjKZLNblNYBycnIkSWfPntW9996rnTt3Wix79+7N96QtT7Vq1XTu3DllZWVd8fx557rS+Z944gn99ddfeuyxx/Tbb7+pZcuWmjZtmsUxT506JT8/v6JVAAAAgI1oRwEoLSSlAEDSLbfcot9//13BwcFq0KCBxVKhQoUC9wkJCZEk7d69u1hiqFOnjvr376+lS5dq+PDhmjt3rsX2Xbt2qUWLFsVyLgAAgOJCOwqArUhKAYByJ8g8deqUevTooZ9//ln79+/XN998o759+yo7O7vAffz8/HTLLbdow4YN13z+oUOH6ptvvtGBAwe0fft2rVu3To0bNzZvP3jwoJKTkxUREXHN5wIAAChOtKMA2IqkFAAod/LLn376SdnZ2erUqZOaNWumoUOHqlKlSnJzK/x/lU888YQ++uijaz5/dna2BgwYoMaNGysqKko33HCD3n77bfP2jz/+WJ06dVJQUNA1nwsAAKA40Y4CYCuTkfceTQCA1c6fP68bb7xRixcvVnh4eImcIysrSw0bNtTChQvVpk2bEjkHAABAaaMdBYCeUgBwDcqVK6cFCxbo5MmTJXaOw4cPa9SoUTSkAACAS6EdBYCeUgAAAAAAACh19JQCAAAAAABAqSMpBQAAAAAAgFJHUgoAAAAAAACljqQUAAAAAAAASh1JKQAAAAAAAJQ6klIAAAAAAAAodSSlAAAAAAAAUOpISgEAAAAAAKDUkZQCAAAAAABAqSMpBQAAAAAAgFL3/wDWpdbC+2pW9QAAAABJRU5ErkJggg=="
     },
     "metadata": {},
     "output_type": "display_data",
     "jetTransient": {
      "display_id": null
     }
    }
   ],
   "execution_count": 4
  },
  {
   "cell_type": "markdown",
   "id": "phase-portrait",
   "metadata": {},
   "source": [
    "### Step 3: Phase Portrait\n",
    "\n",
    "The (r, v) phase space reveals the attractor structure:"
   ]
  },
  {
   "cell_type": "code",
   "id": "plot-phase",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2026-03-21T13:55:38.905722700Z",
     "start_time": "2026-03-21T13:55:37.673455200Z"
    }
   },
   "source": [
    "# Run multiple trajectories with different initial conditions\n",
    "fig, ax = plt.subplots(figsize=(8, 6))\n",
    "\n",
    "init_conditions = [\n",
    "    (0.01, -3.0),\n",
    "    (0.5, -1.0),\n",
    "    (1.0, 0.0),\n",
    "    (0.1, -5.0),\n",
    "    (2.0, -2.0),\n",
    "]\n",
    "\n",
    "for r0, v0 in init_conditions:\n",
    "    model = brainmass.MontbrioPazoRoxinStep(in_size=1, tau=1.0 * u.ms, eta=-5.0, delta=1.0 * u.Hz, J=15.0)\n",
    "    brainstate.nn.init_all_states(model)\n",
    "    model.r.value = np.array([r0]) * u.Hz\n",
    "    model.v.value = np.array([v0])\n",
    "\n",
    "\n",
    "    def step(i):\n",
    "        with brainstate.environ.context(i=i, t=i * brainstate.environ.get_dt()):\n",
    "            model.update()\n",
    "            return model.r.value, model.v.value\n",
    "\n",
    "\n",
    "    r, v = brainstate.transform.for_loop(step, np.arange(int(50 * u.ms / brainstate.environ.get_dt())))\n",
    "\n",
    "    ax.plot(v[:, 0], r[:, 0], 'b-', linewidth=0.8, alpha=0.7)\n",
    "    ax.plot(v[0, 0], float(r[0, 0] / u.Hz), 'go', markersize=6)\n",
    "    ax.plot(v[-1, 0], float(r[-1, 0] / u.Hz), 'r*', markersize=10)\n",
    "\n",
    "ax.set_xlabel('Mean potential v', fontsize=12)\n",
    "ax.set_ylabel('Firing rate r (Hz)', fontsize=12)\n",
    "ax.set_title('Phase Portrait: Trajectories Converge to Fixed Point', fontsize=14)\n",
    "ax.set_xlim([-6, 2])\n",
    "ax.set_ylim([0, 3])\n",
    "plt.tight_layout()\n",
    "plt.show()"
   ],
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<Figure size 800x600 with 1 Axes>"
      ],
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAxYAAAJOCAYAAAAqFJGJAAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjgsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvwVt1zgAAAAlwSFlzAAAPYQAAD2EBqD+naQAAcXdJREFUeJzt3XlcVGX///H3gAKigPuOiEummXuamoJm4pJpppndpVmZlWuWlWYulVlppbeZWv3SMs0tl8rSzNwqtbKs1PTOBbfEXTYVFK7fH/NlcmRQmIWB4fV8POYBc86ZOZ9zgJl5c13XuSzGGCMAAAAAcIGftwsAAAAAkP8RLAAAAAC4jGABAAAAwGUECwAAAAAuI1gAAAAAcBnBAgAAAIDLCBYAAAAAXEawAAAAAOAyggUAAAAAlxEs4NMeeughWSwWxcbGersU5NCcOXNksVg0Z84cb5fidr58bJLvHx98g8ViUXR0tLfLkJR33qtiY2NlsVj00EMPebUO5F8EC+QrGS96V94CAgIUHh6u+++/X3/88Ye3S/S4q4+/UKFCqlChgrp166aNGzd6bL955Q3HXR8GMj78Zvfm7eN2xfr162WxWDRu3Dhvl+IVly9f1uzZs9WpUyeVL19eAQEBCgsL0y233KLRo0fr4MGD3i6xQKhataqqVq2aK/u51t/y8uXLPV6Dp3nrfSC7xo0bJ4vFovXr13u7FOSyQt4uAHBG9erV9cADD0iSkpKStGXLFn366adaunSp1q5dq5YtW3q5Qs8qVaqUBg0aJEm6ePGitm/frhUrVujzzz/XwoUL1bNnTy9X6Lq7775bt956qypUqOCR52/QoIHGjh1rtyw2NlYfffSR6tevr27dumXa3l08fWzelpeO7+DBg+ratat+//13lStXTnfccYfCw8OVnJysX3/9Va+99pomT56sHTt2qEaNGt4uF27i7++v0aNHO1x34403SpL++usvBQcH52ZZbuWJ94FKlSrpr7/+UlhYmLvLRQFBsEC+VKNGjUz/fR09erQmTJigF154wef/S1K6dOlMx//BBx+of//+evbZZ30iWISFhXn0za1BgwaZwsL69ev10UcfqUGDBh79776nj83b8srxJSYmKiYmRnv27NGIESP08ssvKzAw0G6bvXv3avjw4UpKSvJSlfCEQoUKXfdvOCNg5FeeeB8oXLhwvj8v8DID5CMHDhwwkkxMTEymdXFxcUaSCQ4Oti3r27evkWT2799vpk6damrVqmUCAgJMlSpVzLhx40xaWprdc5w7d8689tprpnXr1qZChQqmcOHCpkKFCubBBx80e/fuzbTPCxcumMmTJ5t69eqZ0NBQExwcbCIiIkzPnj3N9u3bM22/fPly07ZtW1O8eHETGBhobrrpJjNp0iRz+fLlbJ8DSaZWrVqZlqelpZmiRYsaSebEiRO25R9++KFp2rSpKVq0qClatKhp2rSpmT17dqbHr1u3zkgyY8eONT/88IO54447TFhYmJFkZs+ebSQ5vK1bt84YY8zYsWNt92fPnm0aNmxoihQpYqKiopw6txn7zKg1oz5HN0fH44yMffTt29duuaeP7Ur79+83jzzyiAkPDzcBAQGmfPnypm/fviY2NtZhzfv27TP9+/c3VatWNQEBAaZMmTImKirK9twZtTu6HThwwPY8J0+eNEOHDrV7np49e5o///wz0z4z/q727dtnJk+ebGrXrm0CAgJs581dx7dt2zZzzz332LYtXbq0adKkiXnllVccnourjRkzxkgyDzzwwHW3TUlJsbv//fffm06dOpkSJUqYwMBAU6tWLTNmzBiTnJyc6bGSTFRUlImLizN9+vQxpUqVMkFBQaZZs2a2v48Mbdu2NRaLJcuf5+DBg40k880339gt37Bhg7nzzjtNqVKlTEBAgKlRo4Z54YUXMtVzrb/jDCdPnjT9+/c3ZcqUMUWKFDFNmjQxS5cuvebP7ffffze9evUy5cuXN4ULFzZVqlQxgwYNMqdOnbrWaTXG/Pu67eg2duxYu22z+3p1LRERESYwMPC622X83DL8/fffplixYqZSpUqZjiurdQkJCWbMmDGmTp06JigoyISFhZn27dubTZs2Odznjh07TOfOnU2xYsVMaGio6dixo/nzzz9tf1NX/k1er3ZPvA9k/Kyufg2MiooykkxqaqoZO3asiYiIMAEBAaZmzZpm+vTpDre9+hYREZGtY0P+RosFfI7FYsm0bMSIEdqwYYPuvPNOxcTEaPny5Ro3bpxSU1M1YcIE23Z//fWXxowZozZt2ujuu+9W0aJFtXv3bs2fP18rV67Ur7/+qoiICNv2ffv21aJFi1SvXj3169dPgYGBOnz4sNatW6eff/5Z9evXt207cuRIvfbaa6pUqZK6d++usLAwbdq0SSNGjNDWrVu1ePFit5+DIUOGaNq0aapUqZIeeeQRSdJnn32mfv366bffftPUqVMzPfbHH3/Uq6++qjZt2uixxx7ToUOH1KBBAw0dOlRTp07N1E3o6j7TkyZN0rp169S1a1e1b99e/v7+knJ+bq9WtWpVjR07VuPHj1dERITdmIcrWx6qVq2qgwcP6sCBA27vz+2pY8uwdetWxcTEKDk5WXfeeadq1qyp2NhYzZs3T19//bU2b96satWq2bb//vvv1blzZ9t/5u+77z6dPXvW9rN96KGHFB0dbeviFRUVZTc+pXjx4pKkkydPqnnz5tq3b5+io6N133336cCBA1qyZIlWrlyp1atX67bbbstU7+DBg7VlyxZ17txZXbp0UdmyZd12fNu3b1eLFi3k7++vrl27KiIiQufOndOuXbv03nvv6YUXXrju+fzwww8lSWPGjLnutgEBAbbvFy9erN69eyswMFC9evVS2bJl9c033+ill17S6tWrtX79egUFBdk9/ty5c7rtttsUFhamBx98UCdOnNDChQsVExOjbdu2qW7dupKkBx98UN99953mzZunUaNG2T3H5cuXtWDBAlWsWFG33367bfmMGTM0cOBAFS9e3Haef/nlF02YMEHr1q3TunXr7OqXHP8dS9auo1FRUdq1a5datGih1q1b68iRI7rvvvsUExPj8Nx8/vnnuvfee+Xn56euXbsqPDxcu3bt0jvvvKPVq1dr69atKlGiRJbntnjx4ho7dqymTJkiSRo2bJht3ZW/j868XrlTjRo1NG3aNPXr10+PPvqoli1bJkm6dOmSevfurfPnz+vzzz9XqVKlJElnzpxR69attXPnTrVs2VKPP/64EhIStGLFCrVp00aLFy+2e63csWOHWrZsqaSkJHXv3l01a9bUTz/9pJYtW9q9V7iLK+8DWendu7d++ukndezYUf7+/lq0aJEGDhyowoULq3///pJke23esGGD+vbta3sdzni9gY/zdrIBcuJaLRYZ/51s06aNbVnGf4EiIyPNP//8Y1t+8uRJU7x4cRMSEmL3n8pz586Z06dPZ3ru7777zvj5+ZlHH33UbluLxWIaN26cqcXh8uXL5uzZs7b733zzja3upKQk2/L09HTz+OOPG0lmyZIl2ToHyuI/VR9++KHtWI2x/odTkqldu7Y5d+6cbbszZ86YG264wUgyGzdutC2/skXgww8/zPT8Wf0nK0PGf8aLFi1q/vjjj0zrc3Jujcn6v9666r+MV4uIiMjRf/6udL0WC08eW2pqqqlataoJCQkxv/76q932mzZtMv7+/ubOO++0Lbt48aKpVKmS8fPzM19//XWmfR8+fDjTcV393+EM/fr1M5LMyJEj7ZavXLnSSDI1atSwa93L+LuqXLmyOXjwYKbnc8fxDR8+3Egyy5cvz/T82fkveWxsrK3GnIiPjzdhYWEmMDDQ/P7777blaWlpplevXkaSeemll+wek/F38+STT9qdpw8++MBIMgMGDLAtS0hIMEWKFDF16tTJtO8vvvjCSDLPPPOMbdnOnTtNoUKFTP369TMd98SJE40kM3nyZNuy6/0djx492kgyjz32mN3yb7/91mEL4KlTp0xoaKipVKlSplaWTz/91EgygwYNyrQfRyIiIrL8r3VOX6+utx9/f38zduzYTLdPP/3Utl1WryX33XefkWTeffddY4wxI0aMcPj3cf/99xtJ5v3337dbfvz4cRMeHm7KlCljLly4YFue8Z/8Tz75xG77kSNHOmxFvBZPvQ9cr8WiWbNmJj4+3rZ89+7dplChQplqubKVFwULwQL5SsaLXvXq1W1vFM8884xp1aqVkWSCgoLMjz/+aNs+4wOQozfYjHWOPig6cvPNN5uqVava7sfHxxtJpmXLliY9Pf2aj73rrruMJIcfwjICyj333JOtOiSZUqVK2Y7/ueeeMx06dDCSjJ+fny2gPPzww0aSWbhwYabnmDdvnpFkHn74YduyjA8kjRo1crjf7AaLp556KlvHcaWrz60xzgeLvXv3mr/++sukpqbmuI7rBQtPHtvSpUsdfmjN0L17d+Pn52d7U1+4cKGRZPr06XPdGq4VLFJSUkxQUJApVaqUw24+d9xxR6YPHxl/O1OnTnW4P3ccX0awWL169XWPz5EtW7YYSebWW2/N0eM+/vhjI8k88cQTmdYdPHjQFCpUyFSrVs1ueUboTExMtFt+6dIlU6hQoUx/U7179zaSzLZt2+yW33vvvUaSXTfKIUOGZPmhOi0tzZQpU8Y0btzYtux6f8cZXd3i4uIyrWvfvn2mn9tbb71lJJmPP/7Y4fM1atTIlC5d2uG6q10rWOT09ep6+8n4oH71rWvXrrbtsnotOXfunKlataopUqSI+e9//2ssFotp2rSpuXTpkm2bkydPGn9/f9O2bVuHNfz3v/81kswXX3xhjLH+7kgy9erVy7RtYmKiKV68eI6DhSfeB64XLL777rtMz5OxLiEhwbaMYFFw0RUK+dK+ffs0fvx4SdbBZuXKldP999+v559/XjfffHOm7Rs3bpxpWeXKlSVZuzBcaf369ZoyZYq2bt2qU6dO6fLly7Z1V3Y3CA0NVadOnfTVV1+pUaNG6tmzp6Kjo3XLLbeocOHCds+5ZcsWFS1a1NY142pFihTR7t27s3fwkk6fPm07fn9/f5UuXVpdu3bV008/rVatWkmSfvvtN0n2XQ0ytGnTRpK1u8nVbrnllmzX4UjTpk2zXJfdc+uK6tWru+V5HPHksW3ZskWStGfPHoeDTuPi4pSenq7//e9/atKkiX766SdJUvv27XN4FPZ2796tixcvqk2bNg6vkNOmTRutWbNG27dvt/1uZbjW+bhaTo/v3nvv1ZQpU3T33XerV69euuOOO9S6dWtVqlQpZweYQ9f6u6lSpYqqVaum//3vf0pMTFRISIht3Q033KBixYrZbV+oUCGVK1cu02vMgw8+qE8//VRz585Vo0aNJEkJCQn64osvdPPNN9t1i8k4b6tXr9batWsz1VS4cGGHrx2O/o4TEhIUGxurOnXqqFy5cpnWt2zZUt98843dsoz9b926Vfv27cv0mIsXL+rUqVM6deqUSpcunWl9djn7epWVwMBAXbx40alawsLCNG/ePLVu3VpDhgxRSEiI5s+fr0KF/v3I9PPPPystLU0pKSkOf5///vtvSda/rzvvvFO///67JDnsUlisWDE1aNAgxxcd8eT7QFau91565d8ECiaCBfKlmJgYrVq1Ktvbh4aGZlqW8SaRlpZmW7Z48WL16tVLxYoVU0xMjKpWrarg4GDbZF9XX+9+8eLFevXVVzV//nxbn+/Q0FD169dPr776qu2D2pkzZ3T58mXbm4AjycnJ2T6eWrVqXTeIJCQkyM/PT2XKlMm0rly5crJYLEpISHC4zhVZPT6n5zYv8uSxnTlzRpI0b968a26X8XsSHx8vSS5/0M74Hcjq2DIuGevq70pOj69Zs2Zav3697e9r9uzZkqwfmF9//XXbh6KslC9fXpJ09OjRbNcoZe98/O9//1NCQoLdhyhHrzGS9XXmytcYyRoGy5UrpwULFmjy5Mny9/fXkiVLdOHCBT344IN222actyvHgmWHo/ozji2rsTCOHpOx/+nTp19zf8nJyS4FC2dfrzylUaNGioiI0P79+9WxY8dM/7DIOC8//PCDfvjhhyyf5+q/15yc++vx5PtAVrL7XoqCi2ABXGHcuHEKCgrStm3bVLNmTbt1CxYsyLR9cHCwXnnlFb3yyis6cOCA1q1bp5kzZ2rq1Km6cOGCZs2aJcn6YmyxWHTq1KlcOY6Mfaanp+vkyZOZ3sxOnDghY4zDNwlHg99zIqvH5/Tc5kWePLaMn8UXX3yhO++887rbZwyEzOkH56z2e/z4cYfr4+Li7La7Uk5+V3J6fJLUqlUrff3117pw4YK2bt2qL774Qu+++646d+6sHTt22A1kv1pERIQqVaqkw4cP6++//870c7lenc6cj+zy9/dX7969NWXKFH377beKiYnR3Llz5efnp/vvv99hPVcHmetx9LPJeK4TJ044fIyjY854zJ9//mkbgO4Jzr5eecqIESO0f/9+lSpVSosWLVLfvn3VqVMnu3ol6emnn9bkyZOv+3wZl1/Oybl3h7x2XuH7mHkbuMK+fftUu3btTB9Cjh07pv3791/zsZGRkXr44Ye1YcMGFStWTJ9//rltXbNmzXT69Glb83huaNiwoSQ5bF7PWJaTSd8yroDk7H+lXDm3V/Lz88tz/xlzx7E1a9ZMkrR58+ZsbZ/RDenqriuOXOtnd+ONNyooKEg///yzzp8/n2m9M78rjuT0+K5UpEgRRUdH680339SoUaN04cIFrVmz5rqPy7gCziuvvHLdbVNTUyVd++/m8OHD2rdvn6pVq+Zyl4+MlolPPvlEhw8f1oYNG9SmTZtMLVAZ5y2jS5IrQkNDVbVqVe3du9fhB9wff/wx0zJXfm5X8/f3z/Jv192vV65YuXKl3nnnHUVFRemXX35RiRIl1K9fP7sP/7fccossFku2z0tG97bvv/8+07qkpKQcdUfKCW+dV1ffL5B/ESyAK0RERGjv3r12byAXL17UE088oUuXLtlte/LkSe3YsSPTc5w9e1YpKSl2l6McMmSIJOnhhx/W6dOnMz0mLi5Of/31l7sOQ5L1UriSNH78eLum7vj4eFuXrIxtsqNEiRKyWCw6fPiwU/Xk5NxeS8mSJXXkyJEs1+/bt0+7d+/O0XO6yh3H1rVrV1WpUkVvvfWWNm7cmGn9pUuX7D6U3HXXXapcubI++eQTrV69OtP2V7ZklCxZUpIc/uwCAgLUu3dvnTp1ShMnTrRbt2rVKq1evVo1atRweTb7nB7f5s2bHfaRzzjHV1/u1ZFnnnlGtWrV0scff6xRo0YpJSUl0zYHDhxQt27dtGvXLludYWFhmj17tnbu3Gnbzhij5557TpcvX7a71LGzGjVqpDp16mjZsmWaNWuWjDGZukFJ0pNPPqlChQpp8ODBtkvGXuncuXO2fvTZ8Z///EepqamZZp1fv369w9+jfv36KSQkRC+88ILd+chw/vz5bIeekiVL6tSpUw5/ru5+vXJWXFyc+vXrpxIlSuiTTz5R1apV9d577+nEiRPq27evjDGSrF3t7r33Xv3444+aNGmSbfmVtm7dagvrVapUUevWrfXHH39k6g746quvZhqH4y7eOq/Xes2Bb6MrFHCFwYMHa/DgwWrYsKF69Oihy5cva82aNTLGqH79+rYBeJL1g1vDhg1Vv3591atXT5UqVdLp06e1YsUKXbp0Sc8884xt2w4dOujFF1/Uyy+/rBo1aqhDhw6KiIjQ6dOntXfvXm3atEmvvPKKateu7bZjad26tQYPHqxp06apbt26uueee2SM0WeffaYjR45oyJAhat26dbafr1ixYrrlllu0ceNGPfjgg6pZs6b8/Pz04IMPZmuOhpyc22tp27atFi1apG7duqlhw4by9/fXXXfdpXr16kmSbr/9do/NY5EVdxxbYGCglixZoo4dOyoqKkpt27bVzTffLIvFooMHD2rTpk0qVaqUrU91YGCgFi1apA4dOqhjx47q0KGD6tevr4SEBG3fvl3nz5+3feC88cYbVbFiRS1YsECBgYGqXLmyLBaLBg8erLCwML3++uvasGGDXnnlFf34449q1qyZYmNjtXjxYgUHB2v27Nny83Pt/1A5Pb7XX39d69atU+vWrRUZGamgoCD9+uuvWrt2rapVq6a77777uvsMCQnR6tWr1bVrV02cOFGzZ89W+/btVblyZdv5+eGHH1SoUCFbd5bQ0FC9//776t27t5o1a6ZevXqpTJky+vbbb7Vt2zY1bdpUI0aMcOlcZHjwwQc1cuRIvfHGGwoODtY999yTaZu6devq3Xff1RNPPKFatWqpU6dOql69uhITE7V//35t2LBBDz30kGbOnJmtfT733HP67LPPNHPmTO3YsUOtWrXSkSNHtGjRInXp0kVffPGF3c+6TJky+vTTT9WzZ0/Vr19fHTp00I033qiUlBTFxsZqw4YNatGiRbbGvLVt21a//PKLOnbsqFatWikgIECtW7e23dz5euUMY4z69OmjkydPasmSJbZByT169NAjjzyi//f//p/eeustPf3005Kkd999V3v27NGzzz6ruXPnqnnz5ipevLgOHz6sX375RX///beOHTtmG2s3ffp0tWzZUn369NHy5ctt81j8/PPPatWqlTZt2uT2Y/LWeW3Tpo0sFotGjRqlnTt3KiwsTMWLF9egQYPcvi/kMd64FBXgrGvNY+HItWYzdXQ5vPT0dDNz5kxz0003maCgIFO+fHnzyCOPmBMnTtguqZfh7NmzZty4cbbZlgMCAkzFihVNhw4dHM4rYIwxa9asMV26dDFlypQxhQsXNuXLlzfNmzc3L7/8sjl06FC2jklZXL88Kx9++KG55ZZbTHBwsAkODja33HKLw8vvXm+uA2OM2bNnj+nUqZMpXry4sVgsdufvepcXzMm5NSbry80eO3bM3HvvvaZ06dLGz88v0zaenMfC08dmjDFHjhwxQ4cONTVr1jSBgYEmNDTU1K5d2zz66KNm7dq1mbbfu3eveeSRR0zlypVN4cKFTdmyZU10dHSmy4Nu2bLFREVFmZCQEIfXzD958qQZMmSIiYiIMIULFzalS5c2PXr0uObM21mdY3cc36pVq0yfPn1MrVq1TEhIiClWrJipU6eOGTVqlDl58qTD/WYlNTXVfPjhh6ZDhw6mXLlypnDhwiYkJMQ0atTIjBo1yuHf3saNG03Hjh1N8eLFTUBAgLnhhhvMiy++aDcPTQZd4xLI17rE6qFDh2y/w717977mMfz000/mvvvuMxUrVrT9fBo1amSef/5589dff9m2y87f8YkTJ8wjjzxiSpcubYKCgkzjxo3N0qVLzeTJk40ks2zZskyP2b17t3nkkUdsMy6XKFHC3HzzzWbIkCHmp59+umbtGRITE03//v1NhQoVjL+/v8M6s/t6dS3Ozrw9adIkIynTvDPGGJOUlGRuuOEGExAQYDcPy/nz580bb7xhGjdubIoWLWqKFCliIiMjTbdu3czHH39sd4laY4z5888/TadOnUyxYsVMSEiI22fezkp2z+v1LjfrSFa1z5kzx9x8880mMDDQSMy8XVBYjHHQfgcA8JiZM2fqiSee0Pz589W7d29vlwNIkh544AHNmzdPu3btcmvrKYCCgzEWAJDL9u7dK+nf678DuenYsWOZlm3YsEELFixQrVq1CBUAnEaLBQDkkjVr1mj58uWaPXu2SpQoodjY2EyTKQKe1rBhQxUpUkQNGjRQ0aJFtWvXLq1atUr+/v5auXKl7rjjDm+XCCCfosUCAHLJypUrNW/ePDVv3lzffPMNoQJe0bdvX126dEkLFizQlClT9NNPP6lLly7atGkToQKAS/JUi8WMGTM0Y8YMxcbGSpJuuukmjRkzRh07dszyMYsXL9aLL76o2NhY1axZU6+//rrdJDYAAAAAPC9PtVhUrlxZr732mrZt26ZffvlFbdu2VdeuXR1eO1uyTubTu3dvPfLII/rtt9/UrVs3devWzeHcAgAAAAA8J0+1WDhSsmRJTZo0yTaD6pV69eql5ORkffnll7Zlt956qxo0aJDta3oDAAAAcF2enSAvLS1NixcvVnJyspo3b+5wm82bN2v48OF2y2JiYrR8+fIsnzclJcVu9tX09HSdOXNGpUqVksVicUvtAAAAQH5kjFFiYqIqVqyY48lR81yw+PPPP9W8eXNdvHhRxYoV07Jly1SnTh2H28bFxalcuXJ2y8qVK6e4uLgsn3/ixIm2aewBAAAAZHb48OEcXxY9zwWLWrVqafv27YqPj9eSJUvUt29fbdiwIctwkVMjR460a+WIj49XlSpVdPjwYYWGhrplHwAAAEB+lJCQoPDwcIWEhOT4sXkuWAQEBKhGjRqSpMaNG+vnn3/W1KlTNWvWrEzbli9fXsePH7dbdvz4cZUvXz7L5w8MDFRgYGCm5aGhoQQLAAAAQHJqiECeuiqUI+np6XZjIq7UvHlzrV271m7ZmjVrshyTAQAAAMAz8lSLxciRI9WxY0dVqVJFiYmJmj9/vtavX6/Vq1dLkvr06aNKlSpp4sSJkqShQ4cqKipKb775pjp37qwFCxbol19+0XvvvefNwwAAAAAKnDwVLE6cOKE+ffro2LFjCgsLU7169bR69WrbTKCHDh2yG53eokULzZ8/X6NHj9aoUaNUs2ZNLV++XHXr1vXWIQAAAAAFUp6fx8LTEhISFBYWpvj4eMZYAAAAoEBz5bNxnh9jAQAAACDvI1gAAAAAcBnBAgAAAIDLCBYAAAAAXEawAAAAAOAyggUAAAAAlxEsAAAAALiMYAEAAADAZQQLAAAAAC4jWAAAAABwGcECAAAAgMsIFgAAAABcRrAAAAAA4DKCBQAAAACXESwAAAAAuIxgAQAAAMBlBAsAAAAALiNYAAAAAHAZwQIAAACAywgWAAAAAFxGsAAAAADgMoIFAAAAAJcRLAAAAAC4jGABAAAAwGUECwAAAAAuI1gAAAAAcBnBAgAAAIDLCBYAAAAAXEawAAAAAOAyggUAAAAAlxEsAAAAALiMYAEAAADAZQQLAAAAAC4jWAAAAABwGcECAAAAgMsIFgAAAABcRrAAAAAA4DKCBQAAAACXESwAAAAAuIxgAQAAAMBlBAsAAAAALiNYAAAAAHAZwQIAAACAywgWAAAAAFxGsAAAAADgMoIFAAAAAJcRLAAAAAC4jGABAAAAwGUECwAAAAAuI1gAAAAAcBnBAgAAAIDLCBYAAAAAXEawAAAAAOAyggUAAAAAlxEsAAAAALiMYAEAAADAZQQLAAAAAC4jWAAAAABwGcECAAAAgMsIFgAAAABcRrAAAAAA4DKCBQAAAACXESwAAAAAuIxgAQAAAMBlBAsAAAAALiNYAAAAAHAZwQIAAACAywgWAAAAAFxGsAAAAADgMoIFAAAAAJcRLAAAAAC4jGABAAAAwGV5KlhMnDhRt9xyi0JCQlS2bFl169ZNe/bsueZj5syZI4vFYncLCgrKpYoBAAAASHksWGzYsEEDBw7Uli1btGbNGl26dEnt27dXcnLyNR8XGhqqY8eO2W4HDx7MpYoBAAAASFIhbxdwpVWrVtndnzNnjsqWLatt27apdevWWT7OYrGofPnyni4PAAAAQBbyVIvF1eLj4yVJJUuWvOZ2SUlJioiIUHh4uLp27aqdO3fmRnkAAAAA/k+eDRbp6ekaNmyYWrZsqbp162a5Xa1atfThhx9qxYoV+uSTT5Senq4WLVroyJEjDrdPSUlRQkKC3Q0AAACAayzGGOPtIhx54okn9PXXX+v7779X5cqVs/24S5cuqXbt2urdu7defvnlTOvHjRun8ePHZ1oeHx+v0NBQl2oGAAAA8rOEhASFhYU59dk4T7ZYDBo0SF9++aXWrVuXo1AhSYULF1bDhg21d+9eh+tHjhyp+Ph42+3w4cPuKBkAAAAo0PLU4G1jjAYPHqxly5Zp/fr1ioyMzPFzpKWl6c8//1SnTp0crg8MDFRgYKCrpQIAAAC4Qp4KFgMHDtT8+fO1YsUKhYSEKC4uTpIUFhamIkWKSJL69OmjSpUqaeLEiZKkl156Sbfeeqtq1Kihc+fOadKkSTp48KAeffRRrx0HAAAAUNDkqWAxY8YMSVJ0dLTd8tmzZ+uhhx6SJB06dEh+fv/24Dp79qz69++vuLg4lShRQo0bN9aPP/6oOnXq5FbZAAAAQIGXZwdv5xZXBqgAAAAAvsTnBm8DAAAAyF8IFgAAAABcRrAAAAAA4DKCBQAAAACXESwAAAAAuIxgAQAAAMBlBAsAAAAALiNYAAAAAHAZwQIAAACAywgWAAAAAFxGsAAAAADgMoIFAAAAAJcRLAAAAAC4jGABAAAAwGUECwAAAAAuI1gAAAAAcBnBAgAAAIDLCBYAAAAAXEawAAAAAOAyggUAAAAAlxEsAAAAALiMYAEAAADAZQQLAAAAAC4jWAAAAABwGcECAAAAgMsIFgAAAABcRrAAAAAA4DKCBQAAAACXESwAAAAAuIxgAQAAAMBlBAsAAAAALiNYAAAAAHAZwQIAAACAywgWAAAAAFxGsAAAAADgMoIFAAAAAJcRLAAAAAC4jGABAAAAwGUECwAAAAAuI1gAAAAAcBnBAgAAAIDLCBYAAAAAXEawAAAAAOAyggUAAAAAlxEsAAAAALiMYAEAAADAZQQLAAAAAC4jWAAAAABwGcECAAAAgMsIFgAAAABcRrAAAAAA4DKCBQAAAACXESwAAAAAuIxgAQAAAMBlBAsAAAAALiNYAAAAAHAZwQIAAACAywgWAAAAAFxGsAAAAADgMoIFAAAAAJcRLAAAAAC4jGABAAAAwGUECwAAAAAuI1gAAAAAcBnBAgAAAIDLCBYAAAAAXEawAAAAAOAyggUAAAAAlxEsAAAAALiMYAEAAADAZQQLAAAAAC4jWAAAAABwWZ4KFhMnTtQtt9yikJAQlS1bVt26ddOePXuu+7jFixfrxhtvVFBQkG6++WZ99dVXuVAtAAAAgAx5Klhs2LBBAwcO1JYtW7RmzRpdunRJ7du3V3JycpaP+fHHH9W7d2898sgj+u2339StWzd169ZNO3bsyMXKAQAAgILNYowx3i4iKydPnlTZsmW1YcMGtW7d2uE2vXr1UnJysr788kvbsltvvVUNGjTQzJkzr7uPhIQEhYWFKT4+XqGhoW6rHQAAAMhvXPlsnKdaLK4WHx8vSSpZsmSW22zevFnt2rWzWxYTE6PNmzd7tDYAAAAA/yrk7QKykp6ermHDhqlly5aqW7dultvFxcWpXLlydsvKlSunuLg4h9unpKQoJSXFdj8hIcE9BQMAAAAFWJ5tsRg4cKB27NihBQsWuPV5J06cqLCwMNstPDzcrc8PAAAAFER5MlgMGjRIX375pdatW6fKlStfc9vy5cvr+PHjdsuOHz+u8uXLO9x+5MiRio+Pt90OHz7stroBAACAgipPBQtjjAYNGqRly5bpu+++U2Rk5HUf07x5c61du9Zu2Zo1a9S8eXOH2wcGBio0NNTuBgAAAMA1eWqMxcCBAzV//nytWLFCISEhtnESYWFhKlKkiCSpT58+qlSpkiZOnChJGjp0qKKiovTmm2+qc+fOWrBggX755Re99957XjsOAAAAoKDJUy0WM2bMUHx8vKKjo1WhQgXbbeHChbZtDh06pGPHjtnut2jRQvPnz9d7772n+vXra8mSJVq+fPk1B3wDAAAAcK88PY9FbmAeCwAAAMDKZ+exAAAAAJA/5KkxFgCAgictPU2bDm3SscRjqhBSQa2qtJK/n7+3ywIA5BDBAgDgNUv/Wqqhq4bqSMIR27LKoZU1tcNUda/d3YuVAQByiq5QAACvWPrXUvVY1MMuVEjS0YSj6rGoh5b+tdRLlQEAnEGwAADkurT0NA1dNVRGma8fkrFs2KphSktPy+3SAABOIlgAAHLdpkObMrVUXMnI6HDCYW06tCkXqwIAuIJgAQDIdccSj11/oxxsBwDwPoIFACDXVQip4NbtAADeR7AAAOS6VlVaqXJoZVlkcbjeIovCQ8PVqkqrXK4MAOAsggUAINf5+/lraoepkpQpXGTcn9JhCvNZAEA+QrAAAHhF99rdteTeJSpbpJLd8sqhlbXk3iXMYwEA+YzTE+Tt2rVLu3bt0qlTp2SxWFS6dGnVrl1bderUcWd9AAAf1r12d1VM7KrXP92kex9m5m0AyM9yFCzWr1+vOXPm6IsvvtC5c+dkjP31xy0Wi8LCwtSlSxf169dP0dHR7qwVAOCDLiT7q1ZgtHrf7O1KAACuyFawWLVqlV588UVt27ZNdevW1UMPPaTGjRurWrVqKlGihIwxOnv2rA4cOKBt27ZpzZo1mjt3rho1aqQJEyYoJibG08cBAMinkpKkYsW8XQUAwFXZChY9evTQo48+qrlz5+rGG2/McrvmzZvr/vvvlyTt3r1bM2fOVM+ePZWQkOCeagEAPic5mWABAL7AYq7uz+TAmTNnVLJkSad24Mpjc0NCQoLCwsIUHx+v0NBQb5cDAAXOBx9Yvz76qHfrAAC49tk4W1eFciUY5OVQAQDwPrpCAYBvcOpysw8//LBGjhypCxcuOFy/ZcsWPfzwwy4VBgAoGJKSpKJFvV0FAMBVTgWLOXPm6I033lCLFi0UGxubaf2+ffv00UcfuVobAKAAoMUCAHyD0xPkDR06VGfPnlWTJk20Zs0ad9YEAChACBYA4BucDhZNmjTRtm3b1LBhQ3Xq1EmvvfaaO+sCABQQBAsA8A1OBwtJKlWqlFavXq1nnnlGL7zwgnr27Knk5GR31QYAKAAIFgDgG1wKFpLk5+eniRMn6rPPPtOaNWt066236n//+587agMA+LjLl6WUFIIFAPgCl4NFhm7duumnn35Senq6Xn75ZXc9LQDAhyUlWb8SLAAg/8vWzNtXGzt2rOrVq5dp+Q033KCffvpJ48aN06lTp1wuDgDg25KSpIAAqXBhb1cCAHCV08EiK0WLFtWkSZOcLggAUHAwvgIAfIfbukIBAJBTTI4HAL4j2y0WISEhslgs2X5ii8Wi+Ph4p4oCABQMtFgAgO/IdrC455577IJFSkqKFixYoPbt26tChQoeKQ4A4NuSkwkWAOArsh0s5syZY3f/1KlTWrBggZ599lm1bdvW3XUBAAoAWiwAwHc4PcYiJ92iAABwhGABAL6DwdsAAK8hWACA7yBYAAC8hmABAL6DYAEA8BqCBQD4jmwP3l66dKnd/cTERFksFn3//fc6d+6cw8d0797dpeIAAL6NYAEAviPbwaJHjx6yWCwyxtgtHzdunMPtLRaL0tLSXCoOAODbCBYA4DuyHSzWrVvnyToAAAUQwQIAfEe2g0VUVJQn6wAAFDCXL0sXLxIsAMBXMHgbAOAVycnWrwQLAPAN2QoWAwYM0IEDB3L85Pv27dOAAQNy/DgAgO9LSpIKF5YCArxdCQDAHbIVLA4fPqxatWqpY8eOmjNnjg4fPpzltrGxsfrggw/Uvn173XjjjTpy5IjbigUA+I6kJKloUW9XAQBwl2yNsfjqq6/0ww8/aPLkyXrssceUlpamUqVKqWrVqipRooSMMTp79qwOHDigs2fPyt/fX506ddK6det02223efoYAAD5EAO3AcC3ZHvwdsuWLdWyZUudPHlSX375pTZv3qzdu3fbWiRKlSql7t27q3nz5urcubPKli3rsaIBAPlfcjLBAgB8SbaDRYYyZcqoX79+6tevnyfqAQAUELRYAIBv4apQAACvIFgAgG8hWAAAvIJgAQC+hWABAPAKggUA+BaCBQDAKxi8DQC+hWABAPAKWiwAwLfkOFhcvHhR//3vf7Vx40ZP1AMAKCAIFgDgW3IcLIKCgvTcc89pz549nqgHAFBAMPM2APgWp7pC1a1bV7GxsW4uBQBQkNBiAQC+xalgMWHCBM2aNUvffvutu+sBABQA6enS+fMECwDwJTmeeVuS3nnnHZUsWVIxMTGKjIxUZGSkihQpYreNxWLRihUr3FIkAMC3JCdbvxIsAMB3OBUs/vjjD1ksFlWpUkVpaWnau3dvpm0sFovLxQEAfFNSklSokBQY6O1KAADu4lSwYHwFAMAVGQO3+R8UAPgO5rEAAOQ6Bm4DgO8hWAAAch3BAgB8D8ECAJDrCBYA4HsIFgCAXJeczOR4AOBrCBYAgFxHiwUA+B6CBQAg1xEsAMD3OHW52QxHjx7Vxo0bdeLECd1zzz2qXLmy0tLSFB8fr7CwMPn7+7urTgCAD0lKkipV8nYVAAB3cqrFwhij4cOHKzIyUv/5z380fPhw/e9//5MkJSUlqWrVqpo2bZpbCwUA+A5aLADA9zgVLCZNmqSpU6fqmWee0Zo1a2SMsa0LCwtT9+7d9dlnn7mtSACAbyFYAIDvcSpYvP/+++rTp49effVVNWjQINP6evXq2VowAAC4GsECAHyPU8Hi8OHDatGiRZbrixYtqoSEBKeLAgD4NoIFAPgep4JF2bJldfjw4SzXb9u2TVWqVHG6KACA7zJGOn+eYAEAvsapYNG9e3fNnDlT+/fvty2zWCySpG+++UZz5sxRz5493VMhAMCnJCdbwwXBAgB8i1PBYvz48apQoYIaNGigPn36yGKx6PXXX9dtt92mjh07ql69eho1apS7awUA+ICkJMnPTwoK8nYlAAB3cipYhIWFacuWLXr22Wd19OhRBQUFacOGDTp37pzGjh2rTZs2KTg42N21AgB8QMb4iv9r6AYA+AinJ8grUqSIRo8erdGjR7uzHgCAj0tKkooW9XYVAAB3c6rFom3btlq7dm2W69etW6e2bds6XRQAwHdxRSgA8E1OBYv169fr+PHjWa4/ceKENmzYkOPn3bhxo7p06aKKFSvKYrFo+fLl163DYrFkusXFxeV43wCA3EGwAADf5FSwkP69CpQje/fuVUhISI6fMzk5WfXr19f06dNz9Lg9e/bo2LFjtlvZsmVzvG8AQO5ITiZYAIAvyvYYi48++kgfffSR7f4rr7yi999/P9N2586d0x9//KFOnTrluJiOHTuqY8eOOX5c2bJlVbx48Rw/DgCQ+2ixAADflO1gcf78eZ08edJ2PzExUX5+9g0eFotFRYsW1eOPP64xY8a4r8rraNCggVJSUlS3bl2NGzdOLVu2zHLblJQUpaSk2O4zQzgA5K6kJMmJRm0AQB6X7WDxxBNP6IknnpAkRUZGaurUqbrrrrs8Vlh2VKhQQTNnzlSTJk2UkpKiDz74QNHR0dq6dasaNWrk8DETJ07U+PHjc7lSAECGpCSpQgVvVwEAcDeLMcZ4uwhHLBaLli1bpm7duuXocVFRUapSpYrmzp3rcL2jFovw8HDFx8crNDTUlZIBANnw4otSq1ZS+/bergQAcLWEhASFhYU59dnY6XksMiQmJio+Pl7p6emZ1lWpUsXVp8+xpk2b6vvvv89yfWBgoAIDA3OxIgDAlRi8DQC+yelgMWPGDL311lvav39/ltukpaU5+/RO2759uyrQxg4AeRaDtwHANzkVLGbOnKmBAwcqJiZGDz/8sF544QU99dRTCgoK0pw5c1SuXDkNGTIkx8+blJSkvXv32u4fOHBA27dvV8mSJVWlShWNHDlSR48e1ccffyxJmjJliiIjI3XTTTfp4sWL+uCDD/Tdd9/pm2++ceawAAC5gGABAL7JqXkspk2bppiYGH399dd67LHHJEmdO3fWhAkTtGvXLiUmJur06dM5ft5ffvlFDRs2VMOGDSVJw4cPV8OGDW1XmDp27JgOHTpk2z41NVVPP/20br75ZkVFRen333/Xt99+q9tvv92ZwwIAeJgxBAsA8FVOtVjs27dPAwcOlCQVLlxYkvVDviSFhYXp0Ucf1bvvvqunn346R88bHR2ta40lnzNnjt39Z599Vs8++2yO9gEA8J7z563homhRb1cCAHA3p1oswsLCdPnyZUlSaGiogoODdfjwYdv6kJAQxcXFuadCAIDPSE6WLBYpONjblQAA3M2pYFG3bl39/vvvtvu33nqrZsyYoaNHj+rw4cOaNWuWbrjhBrcVCQDwDRndoCwWb1cCAHA3p7pCPfDAA5o5c6ZSUlIUGBio8ePHq127drbLyxYuXFifffaZWwsFAOR/jK8AAN/lVLDo16+f+vXrZ7vfsmVL7dy5U1988YX8/f3Vvn17WiwAAJkkJTG+AgB8VY6DxcWLF/Xee++pQYMGat26tW15tWrVNHToULcWBwDwLbRYAIDvyvEYi6CgID333HPas2ePJ+oBAPgwggUA+C6nB2/Hxsa6uRQAgK9LTiZYAICvcipYTJgwQbNmzdK3337r7noAAD6MFgsA8F1ODd5+5513VLJkScXExCgyMlKRkZEqUqSI3TYWi0UrVqxwS5EAAN+QlCSVKePtKgAAnuBUsPjjjz9ksVhUpUoVpaWlae/evZm2sXCRcgDAVWixAADf5VSwYHwFAMAZBAsA8F1OjbEAAMAZBAsA8F0ECwBAriFYAIDvIlgAAHKFMQQLAPBlBAsAQK64eFFKTydYAICvIlgAAHJFUpJksUhFi3q7EgCAJxAsAAC5IilJCg62hgsAgO8hWAAAcgXjKwDAtzk1j0VkZOQ1J8CzWCwKCgpS5cqV1aZNGw0YMEAlSpRwukgAQP6XlEQ3KADwZU61WERFRalYsWKKjY1VSEiIGjZsqIYNGyokJESxsbEqVqyY6tSpoxMnTmjUqFG6+eabdeDAAXfXDgDIR2ixAADf5lSw6Natm44ePaoNGzbo999/12effabPPvtMv//+u9atW6ejR4/qoYce0m+//abvvvtOZ8+e1ciRI91dOwAgHyFYAIBvcypYjBkzRoMHD1arVq0yrYuKitLAgQM1atQoSVJ0dLQGDBigb7/91rVKAQD5WnIywQIAfJlTweLvv/++5piJkiVL6u+//7bdr127tpKTk53ZFQDAR9BiAQC+zalgUa1aNX300Ue6cOFCpnXnz5/X7NmzFRkZaVv2zz//qEyZMs5XCQDI9wgWAODbnLoq1Lhx43TffffpxhtvVN++fVW9enVJ0t69e/Xxxx/r6NGj+vTTTyVJaWlp+uSTT9SyZUv3VQ0AyHcIFgDg25wKFj179lRwcLBGjhypV155xW5d3bp1NX36dN15552SJGOMvv32Wy43CwAFHMECAHybU8FCkjp37qzOnTvr2LFjOnjwoCQpIiJCFSpUsN9BoUKKiIhwrUoAQL5HsAAA3+Z0sMhQoUKFTGECAICrESwAwLc5HSzS0tK0evVq7d+/X2fPnpUxxm69xWLRiy++6HKBAID8zxguNwsAvs6pYPHLL7/onnvu0ZEjRzIFigwECwBAhpQU6fJlggUA+DKnLjf75JNP6sKFC1q+fLnOnDmj9PT0TLe0tDR31woAyKeSkqxfixb1bh0AAM9xqsXijz/+0IQJE9SlSxd31wMA8EFJSVJwsOTn1L+zAAD5gVMv8ZUrV86yCxQAAFdj4DYA+D6ngsVzzz2n999/XwkJCe6uBwDgg5KT6QYFAL7Oqa5QiYmJKlasmGrUqKH77rtP4eHh8vf3t9vGYrHoqaeeckuRAID8jRYLAPB9FuNEnya/bHSStVgs+WIAd0JCgsLCwhQfH6/Q0FBvlwMAPmnFCmnnTmnUKG9XAgC4Flc+GzvVYnHgwAFnHgYAKKDoCgUAvs+pYBEREeHuOgAAPoyuUADg+7jwHwDA42ixAADfl60Wi8jISPn5+Wn37t0qXLiwIiMjZbFYrvkYi8Wiffv2uaVIAED+RosFAPi+bAWLqKgoWSwW26DtjPsAAGQHLRYA4PuyFSzmzJlzzfsAAFxLcjItFgDg63I8xuL8+fPq3r275s2b54l6AAA+iK5QAOD7chwsgoOD9e233+r8+fOeqAcA4IMIFgDg+5y6KtRtt92mzZs3u7sWAIAPunxZuniRMRYA4OucChbvvPOONm3apNGjR+vIkSPurgkA4EMyGrhpsQAA32YxxpicPigkJESXL19WamqqJKlQoUIKDAy0f2KLRfHx8e6p0oNcmbYcAHB9//wjDRwoLV0qcUFBAMjbXPls7NTM2/fccw+XmwUAZEvGFaF42wAA3+ZUsOByswCA7EpKYnwFABQETo2xAAAgu5jDAgAKhmy1WHz88ceSpAcffFAWi8V2/3r69OnjfGUAAJ/ApWYBoGDI1uBtPz8/WSwWXbhwQQEBAfLzu35Dh8ViUVpamluK9CQGbwOAZy1ZIh04II0Y4e1KAADX4/HB2wcOHJAkBQQE2N0HAOB66AoFAAVDtoJFRESERo0apfvuu0/16tVTRESEp+sCAPiIpCSJBmEA8H3ZHrz92muvaceOHbb7p0+flr+/v7777juPFAYA8A1cFQoACgaXrgrlxNx6AIAChq5QAFAwcLlZAIBH0WIBAAUDwQIA4FHJyQQLACgIcjTzdmxsrH799VdJUnx8vCTp77//VvHixR1u36hRI9eqAwDke8xjAQAFQ7bmsZD+ncviSsaYTMuuXM48FgBQsBkj3X23NGOGVKGCt6sBAFyPx+exkKTZs2fnuDAAQMGWkiKlpdFiAQAFQbaDRd++fT1ZBwDAByUlWb8yxgIAfB+DtwEAHpOUJAUHS3682wCAz+OlHgDgMVwRCgAKDoIFAMBjuCIUABQcBAsAgMfQYgEABQfBAgDgMcnJtFgAQEFBsAAAeExSEi0WAFBQECwAAB5DiwUAFBwECwCAxzB4GwAKDoIFAMBj6AoFAAVHtmfeBgAgp+gKBV+Vlp6mTYc26VjiMVUIqaBWVVrJ38/f22UBXpWnWiw2btyoLl26qGLFirJYLFq+fPl1H7N+/Xo1atRIgYGBqlGjhubMmePxOgEA2UNXKPiipX8tVdWpVdXmoza6f+n9avNRG1WdWlVL/1rq7dIAr8pTwSI5OVn169fX9OnTs7X9gQMH1LlzZ7Vp00bbt2/XsGHD9Oijj2r16tUerhQAkB10hYKvWfrXUvVY1ENHEo7YLT+acFQ9FvUgXKBAy1NdoTp27KiOHTtme/uZM2cqMjJSb775piSpdu3a+v777/X2228rJibGU2UCALKJrlDIC4yRUlKsQTcpSUpMlC5elFJTrbdLl6TLl63b+vn9+7VwYalIESkoyHoLDErT4K+Gyshk3oeMLLJo2Kph6lqrK92iUCDlqWCRU5s3b1a7du3slsXExGjYsGFZPiYlJUUpKSm2+wkJCZ4qDwAKtLQ06cIFWizgecZIJ09K//xjvcXFSadP/3s7c8YaICTr72OxYtbAEBBgvRUuLBUq9O9zGWP9/b10yfo7nJJiDSIHtEn/RB7Jug4ZHU44rE2HNim6arTnDxzIY/J1sIiLi1O5cuXslpUrV04JCQm6cOGCihQpkukxEydO1Pjx43OrRAAosM6ft36lxQLuFB8vHThgfzt61NriUK6cVLGiVL68VLWq1LixVLq0VLKkFBZmDRV+LnQC//TPY9qYjZ5OxxKPOb8TIB/L18HCGSNHjtTw4cNt9xMSEhQeHu7FigDANyUlWf8LHBDg7UqQX6WlSfv3S7t3S3v2WG9xcdbgEBlpvTVvLlWpYg0VhTz8qaZCSAW3bgf4mnwdLMqXL6/jx4/bLTt+/LhCQ0MdtlZIUmBgoAIDA3OjPAAo0DIGblss3q4E+YUxUmys9Mcf1tuOHdbfnxtvlGrVkm6/XbrhBu91r2tVpZUqh1bW0YSjDsdZWGRR5dDKalWllReqA7wvXweL5s2b66uvvrJbtmbNGjVv3txLFQEAMjBwG9lx8aK0fbv000/Stm3WLnQ33STVry/df7+1VcKV7kvu5O/nr6kdpqrHoh6yyGIXLiyyJugpHaYwcBsFVp4KFklJSdq7d6/t/oEDB7R9+3aVLFlSVapU0ciRI3X06FF9/PHHkqTHH39c77zzjp599lk9/PDD+u6777Ro0SKtXLnSW4cAAPg/zGGBrCQlSVu3St9/bw0V5cpJTZpIw4dbQ4WnuzS5onvt7lpy7xINXTXU7pKzlUMra0qHKepeu7sXqwO8K0/96f7yyy9q06aN7X7GWIi+fftqzpw5OnbsmA4dOmRbHxkZqZUrV+qpp57S1KlTVblyZX3wwQdcahYA8gDmsMCVLl6UfvxR2rhR+v136+Dqli2l/v2tA67zk+61u6trra7MvA1cxWKMydxJsABJSEhQWFiY4uPjFRoa6u1yAMBnfPaZtG+f9Oyz3q4E3pKebh0r8d131lBRoYIUHS3ddpu1lQJA3uPKZ+M81WIBAPAddIUquM6ckb75xnq7dMkaJiZNso6XAOC7CBYAAI9ITqYrVEFijLV1YuVK6eefrYOvBwywjp3wp4cQUCAQLAAAHpGURHeXgiA1VdqwQVqxQjp7VoqJkR59VCpb1tuVAchtBAsAgEfQYuHbEhOtrRNffmmd1bprV2uXJyZEBAouggUAwCMYY+GbzpyRli+Xvv7aOlnd009LDRowESIAggUAwEOYIM+3nD4tLV4srVkjNWwovfKKdTZsAMhAsAAAeATzWPiGc+ekJUusLRRNmkhvvmmdgwIArkawAAC4nTG0WOR3Fy9a5yJZvtx6hadJk6Rq1bxdFYC8jGABAHC71FTp8mWCRX6Unm6d0G7uXOtVvejyBCC7CBYAALdLSrJ+DQ72bh3ImT/+kP7f/7O2NvXvL7VsyaBsANlHsAAAuF1yslSkCBOj5RdHj0qzZ0t//in16iXdeSeXjQWQcwQLAIDbMXA7f7h4UfrkE+vA7DvukN57zzonBQA4g2ABAHA7Bm7nfT/9JM2YIZUvL02ZIoWHe7siAPkdwQIA4HZMjpd3nTkjzZplHU/x8MNSu3aMowDgHgQLAIDb0RUq7zHG2uXpo4+kpk2lmTPp9gTAvQgWAAC3oytU3nLwoDRtmnWyu+eft86cDQDuRrAAALgdXaHyhvR0adkyaf58qUsXqXdvKTDQ21UB8FUECwCA2yUlWSdXg/ecOCG99ZZ09qz06qtMcgfA8/y8XQAAwPfQFcp7jJHWrpUGD5aqVJGmTiVUAMgdtFgAANyOwdveER8vTZ8u/fWXNGKE1KSJtysCUJDQYgEAcDtaLHLftm3SoEHW76dPJ1QAyH20WAAA3I4Wi9yTlibNmyd98YX0+ONS27bMSwHAOwgWAAC3o8Uid5w5I02aZO0C9dZbzJ4NwLvoCgUAcKv0dOn8eVosPO2PP6QhQ6QyZQgVAPIGWiwAAG6VnGz9SouFZxgjLVwoLVkiPfaYdMcddH0CkDcQLAAAbpWUJPn7MxGbJyQmWrs+xcVJb7whVavm7YoA4F90hQIAuFVysrUbFP9Fd6+DB6WnnpKCgqQpUwgVAPIeWiwAAG6VlEQ3KHfbulV6803p7rul++4jtAHImwgWAAC3ymixgOuMkRYvtt6eekpq0cLbFQFA1ggWAAC3osXCPVJSpKlTpd27pddfp+sTgLyPYAEAcCtaLFx36pT0yivWAfBvvy2FhXm7IgC4PgZvAwDcisnxXLN/v/T009YWigkTCBUA8g9aLAAAbkVXKOdt3y5NnCjdc4/UsyeDtAHkLwQLAIBbJSVZZ4NGzqxdK82YIQ0cKLVp4+1qACDnCBYAALeiK1TOZMykvWyZ9OKLUv363q4IAJxDsAAAuBVdobLv8mXp3XelX3+1XvmpalVvVwQAziNYAADcKimJq0JlR2qq9Npr0smT0uTJUunS3q4IAFxDsAAAuBVdoa7v/Hnp5ZeltDRruCCIAfAFXG4WAOA2xtBicT0JCdLo0VJAgPTSS5wrAL6DFgsAgNtcumQdN0CLhWOnT1sHaEdEWOeqKMS7MAAfQosFAMBtkpKsX4ODvVtHXnTsmPTss9KNN0ojRhAqAPgeXtYAAG6TnCwFBfGh+WqHD0svvCBFR0v9+jHxHQDfxEs/AMBtuNRsZgcPWkNFp05S796ECgC+i2ABAHAbgoW92FhrqOjSRbrvPm9XAwCexRgLAIDbcEWofx04II0aJXXtSqgAUDAQLAAAbsMcFlb791tbKu6+W7r3Xm9XAwC5g65QAAC3ocXCGipGj5a6d5d69PB2NQCQewgWAAC3SU4u2MHi4EFrqLjnHusNAAoSukIBANymIA/ePnrUGiruuotQAaBgIlgAANymoLZYHD9uDRXt2km9enm7GgDwDoIFAMBtCmKLxalT1oHazZtLffowTwWAgotgAQBwm4J2Vahz56wtFfXrS/37EyoAFGwECwCA2xSkq0IlJ0svvijVqCENHEioAACCBQDAbQpKi0VqqvTyy1KZMtKwYZIf76YAQLAAALhHenrBCBZpadIbb1iP97nnpEJcuB0AJDGPBQDATc6ft3715a5QxkjTplmvAjVxohQY6O2KACDvIFgAANwiOdnaJSgoyNuVeM5HH0l//ilNmuT7LTMAkFMECwCAW2QM3PbVQcwrVkjffmvtBlWypLerAYC8h2ABAHALXx5f8eOP0iefSK++KlWs6O1qACBvYvA2AMAtfHVyvN27pbfflkaMkGrW9HY1AJB3ESwAAG6RnOx7A7ePHbNeVvahh6SmTb1dDQDkbQQLAIBb+FqLRWKiNG6c1Lat1Lmzt6sBgLyPYAEAcAtfChYZE+BVrSo9/LC3qwGA/IFgAQBwi4yrQuV3xljHVKSnS08/7btXuQIAdyNYAADcwleuCvXxx9LevdKLL0oBAd6uBgDyDy43CwBwC19osVizRlq1Spo8WQoL83Y1AJC/0GIBAHCL/N5i8ddf0qxZ0siRUqVK3q4GAPIfggUAwC3yc4vFyZPShAlSv35SvXrergYA8ieCBQDALfJri8XFi9Irr0jNm0udOnm7GgDIvwgWAAC3yI8tFsZIU6ZIwcHSgAFcAQoAXMHgbQCAy1JTpUuX8l+LxcKF0t9/Wy8vW4h3RABwCS+jAACXJSVZv+anFosff5SWLpXeeEMKDfV2NQCQ/+XJrlDTp09X1apVFRQUpGbNmumnn37Kcts5c+bIYrHY3YKCgnKxWgBAcrIUGJh//usfG2ttpRg+3Dq7NgDAdXnuLWDhwoUaPny4Zs6cqWbNmmnKlCmKiYnRnj17VLZsWYePCQ0N1Z49e2z3LXSSBYBclZSUf7pBJSdLr74q3X23dOut3q7Gc9LS07Tp0CYdSzymCiEV1KpKK/n7+Xu7LAA+LM8Fi7feekv9+/dXv379JEkzZ87UypUr9eGHH+r55593+BiLxaLy5cvnZpkAgCskJ+ePblDGSG+9ZZ2nondvb1fjOUv/Wqqhq4bqSMIR27LKoZU1tcNUda/d3YuVAfBleaorVGpqqrZt26Z27drZlvn5+aldu3bavHlzlo9LSkpSRESEwsPD1bVrV+3cuTPLbVNSUpSQkGB3AwC4Jr+0WCxeLB08aO0C5auN20v/Wqoei3rYhQpJOppwVD0W9dDSv5Z6qTIAvi5PBYtTp04pLS1N5cqVs1terlw5xcXFOXxMrVq19OGHH2rFihX65JNPlJ6erhYtWujIkSMOt584caLCwsJst/DwcLcfBwAUNPmhxeLXX63BYtQoKSTE29V4Rlp6moauGiojk2ldxrJhq4YpLT0tt0sDUADkqWDhjObNm6tPnz5q0KCBoqKitHTpUpUpU0azZs1yuP3IkSMVHx9vux0+fDiXKwYA35PXWyxOnJAmT5aeeEKqVs3b1XjOpkObMrVUXMnI6HDCYW06tCkXqwJQUOSpMRalS5eWv7+/jh8/brf8+PHj2R5DUbhwYTVs2FB79+51uD4wMFCBgYEu1woA+FdennU7NdU6WLtVK6ltW29X4zlHjkiLVx3L1rbHErO3HQDkRJ5qsQgICFDjxo21du1a27L09HStXbtWzZs3z9ZzpKWl6c8//1SFChU8VSYA4Cp5edbtWbOkwoWl/v29XYn7xcdLy5ZJQ4dKQ4ZIqaez995XIYT3SADul6daLCRp+PDh6tu3r5o0aaKmTZtqypQpSk5Otl0lqk+fPqpUqZImTpwoSXrppZd06623qkaNGjp37pwmTZqkgwcP6tFHH/XmYQBAgZKcLEVEeLuKzNavlzZvlv773/wzx8b1pKVZx4usWSP9/LNUp450113WS+cGFWmlVVMr62jCUYfjLCyyqHJoZbWq0soLlQPwdXnuZbZXr146efKkxowZo7i4ODVo0ECrVq2yDeg+dOiQ/Pz+bWg5e/as+vfvr7i4OJUoUUKNGzfWjz/+qDp16njrEACgwMmLYyyOHpXefVd69lmpdGlvV+O6f/6xhonvvpP8/KTbb5cefliy7ynsr6kdpqrHoh6yyGIXLiyyXgZrSocpzGcBwCMsxpjM/9IoQBISEhQWFqb4+HiFhoZ6uxwAyJeGDrXOC5FXJpxLTZWeflpq1Ej6vwbvfCktzdoq8eWX0s6d1vPbrp3UsKE1XGTF0TwW4aHhmtJhCvNYALgmVz4b57kWCwBA/pPXBm+//74UFCQ9+KC3K3FOfLz0zTfSV19Z73fsKI0YIYWFZe/x3Wt3V9daXZl5G0CuIlgAAFyWlwZvb9woff99/hxXsW+f9Pnn0qZN1rETjz0mNW0q+TuRB/z9/BVdNdrtNQJAVvLZSy4AIK8xRjp/Pm+0WPzzj/TOO9ZuUGXKeLua7DHGOhh76VJpzx7r2ImpUyXmbwWQ3xAsAAAuOX/e+uHY2y0WqanS669L7dtLzZp5t5bsuHTJ2rqydKm161OXLtLzz/vurOAAfB/BAgDgkuRkyWKRihTxbh3/7/9Zuz499JB367ieixelr7+Wli+XgoOlu++WoqOlgABvVwYAriFYAABcknGpWYvFezVs2SJt2GDtQpRXx1WcPy+tXGkNFOXKSU8+aR0/4c3zBgDulEdffgEA+YW357A4c8Y6UPvJJ60f2POapCTrgOzPP5eqVLGO/2jYkEABwPcQLAAALvHmFaGMkd5+W2rSRGrd2js1ZOXiRWnFCusYiho1pFGjpJtvJlAA8F0ECwCAS7w5h8WKFdKxY9LIkd7ZvyOXLkmrVkkLF1pnxR492hooAMDXESwAAC7xVovF/v3SJ59IL79sHQTtbWlp0rp10vz51vMxZIh0yy20UAAoOAgWAACXeKPFIiVFmjxZuuceqXbt3N331YyRNm+W5s6VLl+W+vSxdsvy8/NuXQCQ2wgWAACXeKPFYvZs6z7vvTd393u13bul99+XTp6U7rvPOodGXr0qFQB4Gi9/AACXJCdLJUvm3v5+/tna5WjqVMnfP/f2e6VTp6zhZutWqUcPqVs3KSjIO7UAQF5BsAAAuCQ3LzebmGi9tOxjj1kHRue2ixelzz6Tli2TWrSQZs6USpfO/ToAIC8iWAAAXJKbwWLGDKlWLalt29zZXwZjpPXrpY8+sgaJCROsdQAA/kWwAAC4JDk5d8ZYbNokbd8uvftu7l5paf9+6z5Pn5b69bMOzOZKTwCQGcECAOCS3GixOHvW2loxcKBUvLhn95XhwgVp3jzp66+tYyjuvVcKDMydfQNAfkSwAAC4xNOXmzVGmjZNathQatnSc/u5cn8//ii9955UqZI0ZYoUHu75/QJAfkewAAA4LTXVevNksPjuO2nfPumddzy3jwzHjkmzZkl790qPPCJFR9PtCQCyi2ABAHBaUpL1q6fGWJw8aW05GDFCCgnxzD4k66zZy5ZJCxZYB4Y/80zuT/oHAPkdwQIA4LTkZOv8DZ6YFM4Y66Vlb7tNatLE/c+f4cAB65wYqalc7QkAXEGwAAA4zZOzbn/9tfTPP9LIkZ55/kuXpEWLpKVLpbvvlnr1kgoX9sy+AKAgIFgAAJzmqStCnTwpzZkjjRolBQe7//n37LG2UhQuLE2aJFWr5v59AEBBQ7AAADjNE1eEMkaaPt3aBapBA/c+9+XL1kvIfv65dN991pYKT3TjAoCCiJdTAIDTPNEVav1667iHESPc+7yHDklvvmn9/u23pSpV3Pv8AFDQESwAAE5zd1eoc+ek99+XhgxxX2AxxtpCMXeudNdd0v3300oBAJ7ASysAwGnJye5tsZg1S6pfX7r1Vvc836lT1taJuDjppZekOnXc87wAgMz8vF0AACD/cmeLxebN0u+/SwMGuOf5NmyQBg2SypWzztxNqAAAz6LFAgDgtKQkKTLSPc8zY4bUv79UvLhrz3XxovTuu9K2bdKwYe5r/QAAXBvBAgDgNHd1hfrwQ+slX6OjXXue/fulN96QSpWytlKULOl6bQCA7CFYAACc5o6uUH/8IX3/vfUSsxaLc89hjHVCvQ8/lHr0kO69V/Kjsy8A5CqCBQDAaa4Gi0uXrN2WHnhAKlPGuedITpb++19p925p3Dipbl3n6wEAOI//5wAAnOZqV6glS6SgIKlzZ+ce//ff1kvTpqZawwWhAgC8hxYLAIBT0tOl8+edb7E4etQaLF57TfL3z/njv/lGeu89qXdvqXt357tRAQDcg2ABAHBKcrL1qzPBwhhrF6j27aWaNXP22EuXrPNdbNkijRkj1auX8/0DANyPYAEAcEpSkrWlITAw549dv146ckQaNSpnjzt1Snr1VWvrxJQpUunSOd83AMAzGGMBAHBKxsDtnHZBSkyU/t//s85ZkZPxGX/8IQ0dKlWvbu0+RagAgLyFFgsAgFOSk53rBvXRR1KNGlLLltnb3hhp+XJp3jzrrNx33JHzfQIAPI9gAQBwSlJSzq8I9ddf1m5Q2Z2z4tIl6Z13pN9/lyZOzPl4DABA7iFYAACcktM5LC5ftgaKXr2kcuWuv318vDRhgvVxb73FLNoAkNcxxgIA4JScdoX66itrSLj77utvGxsrDR9unTTvtdcIFQCQH9BiAQBwSk66Qp07Zx0j8fzzUqHrvPP89JM0ebJ1bopevZifAgDyC4IFAMApSUlSSEj2tp0zR2rQQGrYMOttjJGWLZM+/VQaNiz7g7sBAHkDwQIA4JTkZKl8+etvt2eP9P331gnxspKWZp30butW6yDtGjXcVycAIHcQLAAATsnO4O30dGnmTKlHD6lsWcfbXLwoTZokHTtm7QJVpoz7awUAeB7BAgDglOwEizVrpIQE63gJR+LjpZdekgICpDfecG5eDABA3sBVoQAATklOvvbg7cRE62R4jz1mDQ5XO3ZMGjHCeunZl14iVABAfkeLBQDAKddrsfjkE+mGG6SmTTOv27PHGiZuv13q148rPwGALyBYAAByzJhrB4sDB6Rvv5X++9/MoWH7duvEd336SF26eLxUAEAuIVgAAHIsJcU6MNtRVyhjpA8+kDp3lipVsl+3ebP05pvSwIFSmza5UysAIHcwxgIAkGNJSdavjoLFzz9bZ87u1ct++dq10ltvWcdVECoAwPfQYgEAyLGMWbf9rvr31OXL0ocfSv/5j33o+OILae5c6cUXpXr1crdWAEDuIFgAAHIsI1hc7euvrWEjJsZ63xhp0SJpxQrp5ZelWrVyt04AQO6hKxQAIMeSkzMP3E5MlObPlx5+WPL3t4aKjz+WVq60zqZNqAAA30aLBQAgxxxdEWrhQqlGDalx439DxXffWUPF1YO4AQC+hxYLAECOXd0V6p9/rN2gHnnEej8jVLz6KqECAAoKWiyQ56Wlp2nToU06lnhMFUIqqFWVVvL38/d2WUCBdnWLxezZUtu2UkQEoQIACiqCBfK0pX8t1dBVQ3Uk4YhtWeXQypraYaq61+7uxcqAgu3KMRY7dki//y7NmkWoAICCjK5QyLOW/rVUPRb1sAsVknQ04ah6LOqhpX8t9VJlAEImj9Fdj5SSuaO9Pv5Y6t5d+vJLQgUAFGQEC+RJaelpGrpqqIxMpnUZy4atGqa09LTcLg2ApJv/XqrSOqP0b9fqn3+sy1atkiZMIFQAQEFFsECetOnQpkwtFVcyMjqccFhf79qky5dzsTAA0sKFqq2/JEl+StftKV9q2TJp/HipcmUv1wYA8BrGWCBPOpZ4LFvbjZl0TLPOSAEBUnCw9So1jr5ea13G14AADx8U4CsefVSWK1oTuy7rpya/nVSNGl6sCQDgdQQL5EkVQipka7vXR1dQk9LS+fPWwaRXf834/sSJrLfJaPEoVMhxEClSxHoLDpaCgv69f+XyIkXs1xUuLFksHjxBQG6ZMkX64Qf7ZUlJdndLpJ1SyZd62m/TsqU0bJhHSwMA5C0WY0zmTuwFSEJCgsLCwhQfH6/Q0FBvl4P/k5aepqpTq+powlGH4ywssqhyaGUdGHrA5UvPpqZaQ0ZWwePCBeniRevXCxesyzO+v3Ld+fPWScEk66zDVweQnN6uDCyFC7t0iIDzXEnIBfvtBQDyJVc+G9NigTzJ389fUztMVY9FPWSRxS5cWGT9oDOlwxS3zGcREGC9FS/u2vMYYw0pV4aOa92Sk6VTp+yXnT9vH2IyPpcVKmQNGYGB1qCR8TXj+6uXX7ksq+VXfvVjtBWyUr68FBcnIylHEaN8eQ8VBADIq2ixoMUiT3M0j0V4aLimdJji8/NYGCOlpNgHj5QUa/C4eNH6fcbtymVXrstq2cWLUnr6v/sqXNhx8MgIXQEB1m0CA61fXVmecSPMeJ8x1t+F8+etvZvi46Vz5/69nT0rHTkilV/yXw2JHa5CSsteuOjfX3rvPY/WDgDwDFc+GxMsCBZ5HjNve8bly9cOKikp1haYjNulS/b3r7f86nWXLtnv398/c/goVMh6K1z43+8d3XKy3t/fevPzs94yvne0LDvr/fyy7h105XJH21y5zBgpLc0a8K7+6mjZ1V+vPrdXf3/pkvVnmPE1o7vflbcrW8UyAl9GoE1OtgYNPz8pMlKKLv6rHp3RRH4Ouiba+e03qUGD6/7+AQDyJoKFCwgWQO4wJusPwxm3y5ett0uXrB+gM77PWH7l+quXXWt9Vh/Ws/t9xs2TchpwHLUKOWodytgm46IE/v7/tkycOiWdPCnFxUnHj1u3j4yUatSQqleXqlWTwsOtAU2SVLKktRkjKyVKSGfOePZEAQA8ijEWAPI8i+XfD7v5kTHW25UB48p/yzj6F8211l/dAuLOq4hdvmwNDP/8Ix07Zr0dOSIdOmQNEyEh1sAQHi41bvzv96VLX6eOc+euvePrrQcA+DSCBQBkQ8aH/7wyNiQ11drSkBEcrrydOGGts3x5qUIF6+3WW6WePa0BIizMyZ1enY78/a1NO1mtBwAUKAQLAMiB3BrzY4y1V9GJE9ZuSidO2AeJ06etXZwqVJAqVrR+bdHi3+9LlXJzCNqzx/5+xgDtxx6T3n/ffrtatdy4YwBAfkGwAIBscnSVssqhlTW1w9QcX6XMGOtwhYzQcPXXEyesjQElS0ply0rlylm/3nTTv60QxYvn4kSMe/dav1os0q+//jtA+733rCGjWTPrQe3fT7AAgAIqTw7enj59uiZNmqS4uDjVr19f06ZNU9OmTbPcfvHixXrxxRcVGxurmjVr6vXXX1enTp2ytS8GbwPIjqV/LVWPRT0yTdiYMa/KknuX2IWL7ASHy5etwSEjNJQrZ/996dJ5bHLE8+etI8CdXQ8AyPN8avD2woULNXz4cM2cOVPNmjXTlClTFBMToz179qhs2bKZtv/xxx/Vu3dvTZw4UXfeeafmz5+vbt266ddff1XdunW9cAQAfE1aepqGrhrqcBZ4839Txz22dJiOhnXVieP+1wwONWv+e7906Xw2mP16oYFQAQAFWp5rsWjWrJluueUWvfPOO5Kk9PR0hYeHa/DgwXr++eczbd+rVy8lJyfryy+/tC279dZb1aBBA82cOfO6+6PFAsD1rI9drzYftbnuds+XX6eoiGhbcChTJp8FBwBAgefKZ+M8cn0Tq9TUVG3btk3t2rWzLfPz81O7du20efNmh4/ZvHmz3faSFBMTk+X2AJBTxxKPZWu7ei2OqUMHqWFDqVIlQgUAoGDJU12hTp06pbS0NJUrV85uebly5bR7926Hj4mLi3O4fVxcnMPtU1JSlJKSYrsfHx8vyZrOAMCRUIVKF7O3Ha8lAID8LON9zJlOTXkqWOSGiRMnavz48ZmWh4eHe6EaAL7kztfu9HYJAAC4xenTpxWWw4mP8lSwKF26tPz9/XX8+HG75cePH1f58uUdPqZ8+fI52n7kyJEaPny47f65c+cUERGhQ4cO5fjkIXsSEhIUHh6uw4cPM47FQzjHnsc59izOr+dxjj2Pc+x5nGPPi4+PV5UqVVSyZMkcPzZPBYuAgAA1btxYa9euVbdu3SRZB2+vXbtWgwYNcviY5s2ba+3atRo2bJht2Zo1a9S8eXOH2wcGBiowMDDT8rCwMH5BPSw0NJRz7GGcY8/jHHsW59fzOMeexzn2PM6x5/k5MctqngoWkjR8+HD17dtXTZo0UdOmTTVlyhQlJyerX79+kqQ+ffqoUqVKmjhxoiRp6NChioqK0ptvvqnOnTtrwYIF+uWXX/Tee+958zAAAACAAiXPBYtevXrp5MmTGjNmjOLi4tSgQQOtWrXKNkD70KFDdgmqRYsWmj9/vkaPHq1Ro0apZs2aWr58OXNYAAAAALkozwULSRo0aFCWXZ/Wr1+faVnPnj3Vs2dPp/YVGBiosWPHOuweBffgHHse59jzOMeexfn1PM6x53GOPY9z7HmunOM8N0EeAAAAgPwnT02QBwAAACB/IlgAAAAAcBnBAgAAAIDLCBZXWblypZo1a6YiRYqoRIkStvk04B5Vq1aVxWKxu7322mveLssnpaSkqEGDBrJYLNq+fbu3y/Epd911l6pUqaKgoCBVqFBBDz74oP755x9vl+UzYmNj9cgjjygyMlJFihRR9erVNXbsWKWmpnq7NJ8xYcIEtWjRQsHBwSpevLi3y/EZ06dPV9WqVRUUFKRmzZrpp59+8nZJPmPjxo3q0qWLKlasKIvFouXLl3u7JJ8yceJE3XLLLQoJCVHZsmXVrVs37dmzJ8fPQ7C4wmeffaYHH3xQ/fr10++//64ffvhB999/v7fL8jkvvfSSjh07ZrsNHjzY2yX5pGeffVYVK1b0dhk+qU2bNlq0aJH27Nmjzz77TPv27VOPHj28XZbP2L17t9LT0zVr1izt3LlTb7/9tmbOnKlRo0Z5uzSfkZqaqp49e+qJJ57wdik+Y+HChRo+fLjGjh2rX3/9VfXr11dMTIxOnDjh7dJ8QnJysurXr6/p06d7uxSftGHDBg0cOFBbtmzRmjVrdOnSJbVv317Jyck5eyIDY4wxly5dMpUqVTIffPCBt0vxaREREebtt9/2dhk+76uvvjI33nij2blzp5FkfvvtN2+X5NNWrFhhLBaLSU1N9XYpPuuNN94wkZGR3i7D58yePduEhYV5uwyf0LRpUzNw4EDb/bS0NFOxYkUzceJEL1blmySZZcuWebsMn3bixAkjyWzYsCFHj6PF4v/8+uuvOnr0qPz8/NSwYUNVqFBBHTt21I4dO7xdms957bXXVKpUKTVs2FCTJk3S5cuXvV2STzl+/Lj69++vuXPnKjg42Nvl+LwzZ85o3rx5atGihQoXLuztcnxWfHy8SpYs6e0yAIdSU1O1bds2tWvXzrbMz89P7dq10+bNm71YGeCc+Ph4Scrx6y7B4v/s379fkjRu3DiNHj1aX375pUqUKKHo6GidOXPGy9X5jiFDhmjBggVat26dBgwYoFdffVXPPvust8vyGcYYPfTQQ3r88cfVpEkTb5fj05577jkVLVpUpUqV0qFDh7RixQpvl+Sz9u7dq2nTpmnAgAHeLgVw6NSpU0pLS1O5cuXslpcrV05xcXFeqgpwTnp6uoYNG6aWLVuqbt26OXqszweL559/PtNg4atvGf15JemFF17QPffco8aNG2v27NmyWCxavHixl48ib8vuOZak4cOHKzo6WvXq1dPjjz+uN998U9OmTVNKSoqXjyJvy+45njZtmhITEzVy5Ehvl5zv5OT3WJJGjBih3377Td988438/f3Vp08fGeYbvaacnmNJOnr0qDp06KCePXuqf//+Xqo8f3Dm/ALA1QYOHKgdO3ZowYIFOX6sz8+8ffLkSZ0+ffqa21SrVk0//PCD2rZtq02bNum2226zrWvWrJnatWunCRMmeLrUfCu75zggICDT8p07d6pu3bravXu3atWq5akS873snuN7771XX3zxhSwWi215Wlqa/P399Z///EcfffSRp0vNt1z5PT5y5IjCw8P1448/qnnz5p4qMd/L6Tn+559/FB0drVtvvVVz5syRn5/P/y/MJc78Ds+ZM0fDhg3TuXPnPFydb0tNTVVwcLCWLFlidzXJvn376ty5c7RoupnFYtGyZcu4cqcHDBo0SCtWrNDGjRsVGRmZ48cX8kBNeUqZMmVUpkyZ627XuHFjBQYGas+ePbZgcenSJcXGxioiIsLTZeZr2T3Hjmzfvl1+fn4qW7asm6vyLdk9x//973/1yiuv2O7/888/iomJ0cKFC9WsWTNPlpjvufJ7nNHiScvbteXkHB89elRt2rSxtR4TKq7Pld9huCYgIECNGzfW2rVrbR9209PTtXbtWg0aNMi7xQHZYIzR4MGDtWzZMq1fv96pUCEVgGCRXaGhoXr88cc1duxYhYeHKyIiQpMmTZIk9ezZ08vV+YbNmzdr69atatOmjUJCQrR582Y99dRTeuCBB1SiRAlvl+cTqlSpYne/WLFikqTq1aurcuXK3ijJ52zdulU///yzbrvtNpUoUUL79u3Tiy++qOrVq9Na4SZHjx5VdHS0IiIiNHnyZJ08edK2rnz58l6szHccOnRIZ86c0aFDh5SWlmab66ZGjRq21w3kzPDhw9W3b181adJETZs21ZQpU5ScnKx+/fp5uzSfkJSUpL1799ruHzhwQNu3b1fJkiUzvfch5wYOHKj58+drxYoVCgkJsY0NCgsLU5EiRbL/RO6/QFX+lZqaap5++mlTtmxZExISYtq1a2d27Njh7bJ8xrZt20yzZs1MWFiYCQoKMrVr1zavvvqquXjxordL81kHDhzgcrNu9scff5g2bdqYkiVLmsDAQFO1alXz+OOPmyNHjni7NJ8xe/ZsI8nhDe7Rt29fh+d33bp13i4tX5s2bZqpUqWKCQgIME2bNjVbtmzxdkk+Y926dQ5/Z/v27evt0nxCVq+5s2fPztHz+PwYCwAAAACeR6dVAAAAAC4jWAAAAABwGcECAAAAgMsIFgAAAABcRrAAAAAA4DKCBQAAAACXESwAAAAAuIxgAQAAAMBlBAsAALIwbtw4WSwWpx4bHR2t6Oho9xYEAHkYwQIAvGTOnDmyWCyyWCz6/vvvM603xig8PFwWi0V33nmnFyr0nl27dmncuHGKjY31+L7Onz+vcePGaf369R7fFwD4MoIFAHhZUFCQ5s+fn2n5hg0bdOTIEQUGBnqhKu/atWuXxo8fn2vBYvz48Q6DxejRo3XhwgWP1wAAvoBgAQBe1qlTJy1evFiXL1+2Wz5//nw1btxY5cuX91JlKFSokIKCgrxdBgDkCwQLAPCy3r176/Tp01qzZo1tWWpqqpYsWaL777/f4WPS09M1ZcoU3XTTTQoKClK5cuU0YMAAnT171m67FStWqHPnzqpYsaICAwNVvXp1vfzyy0pLS7PbLjo6WnXr1tWuXbvUpk0bBQcHq1KlSnrjjTeydQwWi0WDBg3SvHnzVKtWLQUFBalx48bauHFjpm1/++03dezYUaGhoSpWrJhuv/12bdmyxbZ+zpw56tmzpySpTZs2tu5iV7YofP3112rVqpWKFi2qkJAQde7cWTt37rTbz0MPPaRixYrp6NGj6tatm4oVK6YyZcromWeesR1/bGysypQpI0kaP368bV/jxo2T5HiMxezZs9W2bVuVLVtWgYGBqlOnjmbMmJGt83S1unXrqk2bNpmWp6enq1KlSurRo4dTzwsA3kCwAAAvq1q1qpo3b65PP/3Utuzrr79WfHy87rvvPoePGTBggEaMGKGWLVtq6tSp6tevn+bNm6eYmBhdunTJtt2cOXNUrFgxDR8+XFOnTlXjxo01ZswYPf/885me8+zZs+rQoYPq16+vN998UzfeeKOee+45ff3119k6jg0bNmjYsGF64IEH9NJLL+n06dPq0KGDduzYYdtm586datWqlX7//Xc9++yzevHFF3XgwAFFR0dr69atkqTWrVtryJAhkqRRo0Zp7ty5mjt3rmrXri1Jmjt3rjp37qxixYrp9ddf14svvqhdu3bptttuy9R1Ki0tTTExMSpVqpQmT56sqKgovfnmm3rvvfckSWXKlLGFgrvvvtu2r+7du2d5nDNmzFBERIRGjRqlN998U+Hh4XryySc1ffr0bJ2nK/Xq1UsbN25UXFyc3fLvv/9e//zzT5Y/fwDIkwwAwCtmz55tJJmff/7ZvPPOOyYkJMScP3/eGGNMz549TZs2bYwxxkRERJjOnTvbHrdp0yYjycybN8/u+VatWpVpecbzXWnAgAEmODjYXLx40bYsKirKSDIff/yxbVlKSoopX768ueeee657LJKMJPPLL7/Ylh08eNAEBQWZu+++27asW7duJiAgwOzbt8+27J9//jEhISGmdevWtmWLFy82ksy6devs9pOYmGiKFy9u+vfvb7c8Li7OhIWF2S3v27evkWReeuklu20bNmxoGjdubLt/8uRJI8mMHTs203GNHTvWXP1W6eicxsTEmGrVqtkti4qKMlFRUZm2vdKePXuMJDNt2jS75U8++aQpVqyYw30BQF5FiwUA5AH33nuvLly4oC+//FKJiYn68ssvs+wGtXjxYoWFhemOO+7QqVOnbLfGjRurWLFiWrdunW3bIkWK2L5PTEzUqVOn1KpVK50/f167d++2e95ixYrpgQcesN0PCAhQ06ZNtX///mwdQ/PmzdW4cWPb/SpVqqhr165avXq10tLSlJaWpm+++UbdunVTtWrVbNtVqFBB999/v77//nslJCRccx9r1qzRuXPn1Lt3b7tj9/f3V7NmzeyOPcPjjz9ud79Vq1bZPiZHrjyn8fHxOnXqlKKiorR//37Fx8fn6LluuOEGNWjQQAsXLrQtS0tL05IlS9SlSxe7fQFAXlfI2wUAAKxdctq1a6f58+fr/PnzSktLy7J//d9//634+HiVLVvW4foTJ07Yvt+5c6dGjx6t7777LtOH9qs/BFeuXDnTeIISJUrojz/+yNYx1KxZM9OyG264QefPn9fJkyclWa/AVKtWrUzb1a5dW+np6Tp8+LBuuummLPfx999/S5Latm3rcH1oaKjd/aCgINsYigwlSpTINBYlJ3744QeNHTtWmzdv1vnz5+3WxcfHKywsLEfP16tXL40aNUpHjx5VpUqVtH79ep04cUK9evVyukYA8AaCBQDkEffff7/69++vuLg4dezYUcWLF3e4XXp6usqWLat58+Y5XJ/xQfrcuXOKiopSaGioXnrpJVWvXl1BQUH69ddf9dxzzyk9Pd3ucf7+/g6fzxjj/EG5WUbNc+fOdXi1rEKF7N/WsjomZ+3bt0+33367brzxRr311lsKDw9XQECAvvrqK7399tuZzml29OrVSyNHjtTixYs1bNgwLVq0SGFhYerQoYNbawcATyNYAEAecffdd2vAgAHasmWLXdeYq1WvXl3ffvutWrZsec2uMuvXr9fp06e1dOlStW7d2rb8wIEDbq07Q0ZrwpX+97//KTg42BZ2goODtWfPnkzb7d69W35+fgoPD5ekLGe7rl69uiSpbNmyateunVvqzsnM2l988YVSUlL0+eefq0qVKrbljrpgZVdkZKSaNm2qhQsXatCgQVq6dKm6detWIOcvAZC/McYCAPKIYsWKacaMGRo3bpy6dOmS5Xb33nuv0tLS9PLLL2dad/nyZZ07d07Sv/+tv7LFITU1Ve+++657C/8/mzdv1q+//mq7f/jwYa1YsULt27eXv7+//P391b59e61YscLu6k3Hjx/X/Pnzddttt9m6MhUtWlSSbMeSISYmRqGhoXr11Vftrn6VIaPLVU4EBwc73Jcjjs5pfHy8Zs+eneP9XqlXr17asmWLPvzwQ506dYpuUADyJVosACAP6du373W3iYqK0oABAzRx4kRt375d7du3V+HChfX3339r8eLFmjp1qnr06KEWLVqoRIkS6tu3r4YMGSKLxaK5c+d6rGtT3bp1FRMToyFDhigwMNAWYMaPH2/b5pVXXtGaNWt022236cknn1ShQoU0a9YspaSk2M2Z0aBBA/n7++v1119XfHy8AgMDbXNHzJgxQw8++KAaNWqk++67T2XKlNGhQ4e0cuVKtWzZUu+8806O6i5SpIjq1KmjhQsX6oYbblDJkiVVt25d1a1bN9O27du3V0BAgLp06aIBAwYoKSlJ77//vsqWLatjx445eeasYfGZZ57RM888o5IlS7qtNQYAchMtFgCQD82cOVPvvfeeTpw4oVGjRmnkyJH67rvv9MADD6hly5aSpFKlSunLL79UhQoVNHr0aE2ePFl33HFHtie9y6moqChNmTJFc+fO1ZgxY1SyZEl9/fXXqlevnm2bm266SZs2bVLdunU1ceJEjR8/XhEREVq3bp2aNWtm2658+fKaOXOmTpw4oUceeUS9e/fWrl27JFnHoqxdu1aVKlXSpEmTNHToUC1YsEANGjRQv379nKr9gw8+UKVKlfTUU0+pd+/eWrJkicPtatWqpSVLlshiseiZZ57RzJkz9dhjj2no0KFO7TdD5cqV1aJFCyUmJqp79+4qXLiwS88HAN5gMXlpVB4AIF+yWCwaOHBgjlsLAAC+gxYLAAAAAC4jWAAAAABwGcECAAAAgMu4KhQAwGUM1wMA0GIBAAAAwGUECwAAAAAuI1gAAAAAcBnBAgAAAIDLCBYAAAAAXEawAAAAAOAyggUAAAAAlxEsAAAAALiMYAEAAADAZf8fivJGXze2AWwAAAAASUVORK5CYII="
     },
     "metadata": {},
     "output_type": "display_data",
     "jetTransient": {
      "display_id": null
     }
    }
   ],
   "execution_count": 5
  },
  {
   "cell_type": "markdown",
   "id": "external-input",
   "metadata": {},
   "source": [
    "### Step 4: Response to External Input\n",
    "\n",
    "The QIF model responds to external input $I(t)$:"
   ]
  },
  {
   "cell_type": "code",
   "id": "input-response",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2026-03-21T13:55:39.499103400Z",
     "start_time": "2026-03-21T13:55:38.954454500Z"
    }
   },
   "source": [
    "# Create model with noise for realism\n",
    "model = brainmass.MontbrioPazoRoxinStep(\n",
    "    in_size=1,\n",
    "    tau=1.0 * u.ms,\n",
    "    eta=-5.0,\n",
    "    delta=1.0 * u.Hz,\n",
    "    J=15.0,\n",
    "    noise_r=brainmass.OUProcess(1, sigma=0.1 * u.Hz, tau=1.0 * u.ms),\n",
    ")\n",
    "model.init_all_states()\n",
    "\n",
    "\n",
    "# External input: step function\n",
    "def get_input(t):\n",
    "    \"\"\"Step input at t=50ms\"\"\"\n",
    "    return u.math.where(t > 50. * u.ms, 5.0, 0.0)\n",
    "\n",
    "\n",
    "def step_run(i):\n",
    "    t = i * brainstate.environ.get_dt()\n",
    "    v_inp = get_input(t)\n",
    "    with brainstate.environ.context(i=i, t=t):\n",
    "        model.update(v_inp=v_inp)\n",
    "        return model.r.value, model.v.value\n",
    "\n",
    "\n",
    "# Simulate\n",
    "n_steps = int(150 * u.ms / brainstate.environ.get_dt())\n",
    "r_trace, v_trace = brainstate.transform.for_loop(step_run, np.arange(n_steps))\n",
    "t_ms = np.arange(n_steps) * brainstate.environ.get_dt()"
   ],
   "outputs": [],
   "execution_count": 6
  },
  {
   "cell_type": "code",
   "id": "plot-response",
   "metadata": {
    "ExecuteTime": {
     "start_time": "2026-03-21T13:55:39.636082100Z"
    }
   },
   "source": [
    "fig, axes = plt.subplots(3, 1, figsize=(10, 8), sharex=True)\n",
    "\n",
    "# External input\n",
    "input_trace = [get_input(t) for t in t_ms]\n",
    "axes[0].plot(t_ms, input_trace, 'k-', linewidth=2)\n",
    "axes[0].set_ylabel('External Input')\n",
    "axes[0].set_title('Step Input Response')\n",
    "axes[0].axvline(50, color='gray', linestyle='--', alpha=0.5)\n",
    "\n",
    "# Firing rate\n",
    "axes[1].plot(t_ms, r_trace[:, 0], 'b-', linewidth=1)\n",
    "axes[1].set_ylabel('Firing rate (Hz)')\n",
    "axes[1].axvline(50, color='gray', linestyle='--', alpha=0.5)\n",
    "\n",
    "# Mean potential\n",
    "axes[2].plot(t_ms, v_trace[:, 0], 'r-', linewidth=1)\n",
    "axes[2].set_xlabel('Time (ms)')\n",
    "axes[2].set_ylabel('Mean potential v')\n",
    "axes[2].axvline(50, color='gray', linestyle='--', alpha=0.5)\n",
    "\n",
    "plt.tight_layout()\n",
    "plt.show()"
   ],
   "outputs": [],
   "execution_count": null
  },
  {
   "cell_type": "markdown",
   "id": "bifurcation",
   "metadata": {},
   "source": [
    "### Step 5: Bifurcation Analysis\n",
    "\n",
    "Let's explore how the steady-state firing rate depends on mean excitability $\\bar{\\eta}$:"
   ]
  },
  {
   "cell_type": "code",
   "id": "bifurcation-sim",
   "metadata": {},
   "source": [
    "# Sweep eta values\n",
    "eta_values = np.linspace(-10, 5, 50)\n",
    "steady_r = []\n",
    "steady_v = []\n",
    "\n",
    "for eta_val in eta_values:\n",
    "    model = brainmass.MontbrioPazoRoxinStep(\n",
    "        in_size=1,\n",
    "        tau=1.0 * u.ms,\n",
    "        eta=eta_val,\n",
    "        delta=1.0 * u.Hz,\n",
    "        J=15.0,\n",
    "    )\n",
    "    brainstate.nn.init_all_states(model)\n",
    "    model.r.value = np.array([0.1]) * u.Hz\n",
    "    model.v.value = np.array([-1.0])\n",
    "\n",
    "\n",
    "    # Run to steady state\n",
    "    def step(i):\n",
    "        with brainstate.environ.context(i=i, t=i * brainstate.environ.get_dt()):\n",
    "            model.update()\n",
    "            return model.r.value, model.v.value\n",
    "\n",
    "\n",
    "    r, v = brainstate.transform.for_loop(step, np.arange(int(100 * u.ms / brainstate.environ.get_dt())))\n",
    "\n",
    "    # Take final value (steady state)\n",
    "    steady_r.append(float(r[-1, 0] / u.Hz))\n",
    "    steady_v.append(v[-1, 0])\n",
    "\n",
    "steady_r = np.array(steady_r)\n",
    "steady_v = np.array(steady_v)"
   ],
   "outputs": [],
   "execution_count": null
  },
  {
   "cell_type": "code",
   "id": "plot-bifurcation",
   "metadata": {},
   "source": [
    "fig, axes = plt.subplots(1, 2, figsize=(12, 4))\n",
    "\n",
    "# r vs eta\n",
    "axes[0].plot(eta_values, steady_r, 'b-', linewidth=2)\n",
    "axes[0].axhline(0, color='gray', linestyle='--', alpha=0.5)\n",
    "axes[0].axvline(0, color='gray', linestyle='--', alpha=0.5)\n",
    "axes[0].set_xlabel(r'Mean excitability $\\bar{\\eta}$', fontsize=12)\n",
    "axes[0].set_ylabel('Steady-state firing rate (Hz)', fontsize=12)\n",
    "axes[0].set_title('Firing Rate vs Excitability')\n",
    "axes[0].fill_between(eta_values, 0, steady_r, alpha=0.3)\n",
    "\n",
    "# v vs eta\n",
    "axes[1].plot(eta_values, steady_v, 'r-', linewidth=2)\n",
    "axes[1].axhline(0, color='gray', linestyle='--', alpha=0.5)\n",
    "axes[1].axvline(0, color='gray', linestyle='--', alpha=0.5)\n",
    "axes[1].set_xlabel(r'Mean excitability $\\bar{\\eta}$', fontsize=12)\n",
    "axes[1].set_ylabel('Steady-state mean potential', fontsize=12)\n",
    "axes[1].set_title('Mean Potential vs Excitability')\n",
    "\n",
    "plt.suptitle('QIF Bifurcation Diagram', y=1.02, fontsize=14)\n",
    "plt.tight_layout()\n",
    "plt.show()"
   ],
   "outputs": [],
   "execution_count": null
  },
  {
   "cell_type": "markdown",
   "id": "heterogeneity",
   "metadata": {},
   "source": [
    "### Step 6: Effect of Heterogeneity ($\\Delta$)\n",
    "\n",
    "The parameter $\\Delta$ controls population heterogeneity - how diverse the neurons are:"
   ]
  },
  {
   "cell_type": "code",
   "id": "delta-sweep",
   "metadata": {},
   "source": [
    "delta_values = [0.1, 0.5, 1.0, 2.0, 5.0]\n",
    "\n",
    "fig, axes = plt.subplots(1, len(delta_values), figsize=(15, 3))\n",
    "\n",
    "for idx, delta_val in enumerate(delta_values):\n",
    "    model = brainmass.MontbrioPazoRoxinStep(\n",
    "        in_size=1,\n",
    "        tau=1.0 * u.ms,\n",
    "        eta=-5.0,\n",
    "        delta=delta_val * u.Hz,\n",
    "        J=15.0,\n",
    "        noise_r=brainmass.OUProcess(1, sigma=0.1 * u.Hz, tau=1.0 * u.ms),\n",
    "    )\n",
    "    brainstate.nn.init_all_states(model)\n",
    "    model.r.value = np.array([1.0]) * u.Hz\n",
    "    model.v.value = np.array([-1.0])\n",
    "\n",
    "\n",
    "    def step(i):\n",
    "        with brainstate.environ.context(i=i, t=i * brainstate.environ.get_dt()):\n",
    "            model.update()\n",
    "            return model.r.value\n",
    "\n",
    "\n",
    "    r = brainstate.transform.for_loop(step, np.arange(int(100 * u.ms / brainstate.environ.get_dt())))\n",
    "    t = np.arange(len(r)) * brainstate.environ.get_dt()\n",
    "\n",
    "    axes[idx].plot(t, r[:, 0], 'b-', linewidth=0.8)\n",
    "    axes[idx].set_title(f'$\\Delta$ = {delta_val} Hz')\n",
    "    axes[idx].set_xlabel('Time (ms)')\n",
    "    if idx == 0:\n",
    "        axes[idx].set_ylabel('Firing rate (Hz)')\n",
    "\n",
    "plt.suptitle('Effect of Population Heterogeneity', y=1.02, fontsize=14)\n",
    "plt.tight_layout()\n",
    "plt.show()"
   ],
   "outputs": [],
   "execution_count": null
  },
  {
   "cell_type": "markdown",
   "id": "exercises",
   "metadata": {},
   "source": [
    "## Exercises\n",
    "\n",
    "### Exercise 1: Oscillatory Regime\n",
    "\n",
    "The QIF model can exhibit oscillations under certain parameter regimes. Try:\n",
    "\n",
    "```python\n",
    "model = brainmass.MontbrioPazoRoxinStep(\n",
    "    in_size=1,\n",
    "    tau=1.0 * u.ms,\n",
    "    eta=5.0,     # Increase eta\n",
    "    delta=0.5 * u.Hz,  # Reduce heterogeneity\n",
    "    J=20.0,      # Strong coupling\n",
    ")\n",
    "```\n",
    "\n",
    "Questions:\n",
    "1. Can you find oscillatory dynamics?\n",
    "2. What parameter combinations lead to oscillations?\n",
    "\n",
    "### Exercise 2: Multiple Populations\n",
    "\n",
    "Create E and I populations:\n",
    "\n",
    "```python\n",
    "# Excitatory population\n",
    "E = brainmass.MontbrioPazoRoxinStep(1, eta=-3.0, J=15.0)\n",
    "# Inhibitory population  \n",
    "I = brainmass.MontbrioPazoRoxinStep(1, eta=-3.0, J=10.0)\n",
    "```\n",
    "\n",
    "Couple them and observe E-I dynamics.\n",
    "\n",
    "### Exercise 3: Lorentzian Distribution\n",
    "\n",
    "Visualize the Lorentzian distribution of excitabilities:\n",
    "\n",
    "```python\n",
    "def lorentzian(eta, eta_bar, delta):\n",
    "    return delta / (np.pi * ((eta - eta_bar)**2 + delta**2))\n",
    "```\n",
    "\n",
    "How does $\\Delta$ affect the distribution shape?"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "summary",
   "metadata": {},
   "source": [
    "## Summary\n",
    "\n",
    "In this tutorial, you learned:\n",
    "\n",
    "1. **QIF mean-field model**: An exact reduction of spiking neuron populations\n",
    "2. **Lorentzian heterogeneity**: The $\\Delta$ parameter captures neuron diversity\n",
    "3. **Two state variables**: Population rate $r$ and mean potential $v$\n",
    "4. **Bifurcation behavior**: Transition from low to high firing states\n",
    "\n",
    "### Key Insights\n",
    "\n",
    "- **Exact reduction**: Unlike phenomenological models, QIF is mathematically exact\n",
    "- **Heterogeneity matters**: $\\Delta$ captures realistic neuronal diversity\n",
    "- **Computational efficiency**: 2 ODEs vs thousands of spiking neurons\n",
    "\n",
    "### Comparison with Other Models\n",
    "\n",
    "| Model | Type | Basis | Heterogeneity |\n",
    "|-------|------|-------|---------------|\n",
    "| Wilson-Cowan | Phenomenological | Firing rate | None |\n",
    "| Jansen-Rit | Phenomenological | PSP | None |\n",
    "| QIF | Exact mean-field | Spiking neurons | Lorentzian |\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "references",
   "metadata": {},
   "source": [
    "## References\n",
    "\n",
    "1. Montbrio, E., Pazo, D., & Roxin, A. (2015). Macroscopic description for networks of spiking neurons. *Physical Review X*, 5(2), 021028.\n",
    "\n",
    "2. Gast, R., Schmidt, H., & Knosche, T. R. (2020). A mean-field description of bursting dynamics in spiking neural networks with short-term adaptation. *Neural Computation*, 32(9), 1615-1634.\n",
    "\n",
    "3. Coombes, S., & Byrne, A. (2019). Next generation neural mass models. *Lecture Notes in Nonlinear Dynamics in Computational Neuroscience*."
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "name": "python",
   "version": "3.10.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}