{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "1e249f07",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2026-06-19T08:10:41.690003Z",
     "iopub.status.busy": "2026-06-19T08:10:41.689701Z",
     "iopub.status.idle": "2026-06-19T08:10:47.160937Z",
     "shell.execute_reply": "2026-06-19T08:10:47.159958Z"
    },
    "tags": [
     "remove-cell"
    ]
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "An NVIDIA GPU may be present on this machine, but a CUDA-enabled jaxlib is not installed. Falling back to cpu.\n"
     ]
    }
   ],
   "source": [
    "%matplotlib inline\n",
    "import brainmass\n",
    "import brainstate\n",
    "import brainunit as u\n",
    "import jax.numpy as jnp\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "brainstate.random.seed(0)\n",
    "brainstate.environ.set(dt=0.1 * u.ms)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "136266f7",
   "metadata": {},
   "source": [
    "# Hopf Oscillator\n",
    "\n",
    "The **Hopf (Stuart-Landau normal-form) oscillator** is the canonical model of rhythm onset. Each node obeys the supercritical Hopf normal form\n",
    "\n",
    "$$\\dot x = (a - x^2 - y^2)\\,x - \\omega\\,y,\\qquad \\dot y = (a - x^2 - y^2)\\,y + \\omega\\,x,$$\n",
    "\n",
    "where the bifurcation parameter $a$ controls the regime: for $a<0$ the origin is a stable fixed point (the node is silent), while for $a>0$ a stable limit cycle of radius $\\sqrt{a}$ appears and the node oscillates at angular frequency $\\omega$. It is the workhorse phenomenological node for whole-brain resting-state models.\n",
    "\n",
    "**Reference:** Deco, Kringelbach, Jirsa & Ritter (2017), *The dynamics of resting fluctuations in the brain*, Scientific Reports 7:3095."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "947d858f",
   "metadata": {},
   "source": [
    "## Build the model\n",
    "\n",
    "We pick `a = 0.25 > 0` so the node sits on a limit cycle."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "281bb69d",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2026-06-19T08:10:47.164101Z",
     "iopub.status.busy": "2026-06-19T08:10:47.163512Z",
     "iopub.status.idle": "2026-06-19T08:10:47.183306Z",
     "shell.execute_reply": "2026-06-19T08:10:47.182517Z"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "HopfStep(\n",
       "  in_size=(1,),\n",
       "  out_size=(1,),\n",
       "  init_x=Constant(value=0.0),\n",
       "  init_y=Constant(value=0.0),\n",
       "  method=exp_euler,\n",
       "  a=Const(\n",
       "    fit=False,\n",
       "    t=IdentityT(),\n",
       "    reg=None,\n",
       "    val=Array(0.25, dtype=float32)\n",
       "  ),\n",
       "  w=Const(\n",
       "    fit=False,\n",
       "    t=IdentityT(),\n",
       "    reg=None,\n",
       "    val=Array(0.3, dtype=float32)\n",
       "  ),\n",
       "  beta=Const(\n",
       "    fit=False,\n",
       "    t=IdentityT(),\n",
       "    reg=None,\n",
       "    val=Array(1., dtype=float32)\n",
       "  )\n",
       ")"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "node = brainmass.HopfStep(in_size=1, a=0.25, w=0.3)\n",
    "node"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "59c305b9",
   "metadata": {},
   "source": [
    "## Run a simulation\n",
    "\n",
    "The high-level `Simulator` compiles the whole rollout into one XLA program — no hand-rolled Python loop. We discard a short transient."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "8588a579",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2026-06-19T08:10:47.185474Z",
     "iopub.status.busy": "2026-06-19T08:10:47.185276Z",
     "iopub.status.idle": "2026-06-19T08:10:47.377733Z",
     "shell.execute_reply": "2026-06-19T08:10:47.376981Z"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "((1800, 1), Quantity(200., \"ms\"))"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "sim = brainmass.Simulator(node, dt=0.1 * u.ms)\n",
    "res = sim.run(200. * u.ms, monitors=['x', 'y'], transient=20. * u.ms)\n",
    "res['x'].shape, res['ts'][-1]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "06055409",
   "metadata": {},
   "source": [
    "## Visualize\n",
    "\n",
    "The timeseries shows the sustained oscillation; the phase portrait shows the trajectory settling onto the circular limit cycle of radius $\\sqrt{a} = 0.5$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "eaa3b8c9",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2026-06-19T08:10:47.380670Z",
     "iopub.status.busy": "2026-06-19T08:10:47.380435Z",
     "iopub.status.idle": "2026-06-19T08:10:47.580969Z",
     "shell.execute_reply": "2026-06-19T08:10:47.580094Z"
    }
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA90AAAGGCAYAAABmGOKbAAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjgsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvwVt1zgAAAAlwSFlzAAAPYQAAD2EBqD+naQAARqpJREFUeJzt3XlclWX+//H3AWRRBFIRREHclyRtdERssUa+rk2uaeSCadqi1uSSWSbp1NDkVLa41Mw3NbPRLLW+5mjm0jhJaqjlPma4C7gEuCSiXL8//HHyyOIBzw0cej0fj/NQrnPd5/5c94Fz3e9zn3PfNmOMEQAAAAAAcDmPsi4AAAAAAICKitANAAAAAIBFCN0AAAAAAFiE0A0AAAAAgEUI3QAAAAAAWITQDQAAAACARQjdAAAAAABYhNANAAAAAIBFCN0AAAAAAFiE0A24uf3796tTp04KDAyUzWbTsmXLCux38OBB2Ww2zZ0719724osvymazubSeIUOGKDIy0qWPaaWCtosrHDlyRL6+vvrmm2+KvWy7du30zDPPuLQeAEDpmjt3rmw2m7777ruyLuU3KW/7Hzx4sKxLAQjdQEndaDK955571KJFC8vriI+P144dO/Tyyy9r/vz5atOmjeXrLI4LFy7oxRdf1Pr168u6lFI1depURUdH64477ij2shMmTNCMGTOUmppqQWUAgJuRN//n3Xx9fdW4cWONGjVKaWlpZV2eW/nLX/5S6MECK8ycOdPlb7IDziB0A27sl19+UVJSkoYNG6ZRo0Zp4MCBqlOnjtPLT5o0Sb/88otLa/r73/+uffv22X++cOGCpkyZ8psK3SdPntS8efP02GOPlWj5Hj16KCAgQDNnznRxZQAAV5k6darmz5+vd955R+3bt9esWbMUExOjCxculHVpbsPK0D1o0CD98ssvqlu3rr2N0I2yQugG3NjJkyclSUFBQSVa3svLS76+vi6sSKpUqZJ8fHxc+pju5sMPP5SXl5f++Mc/lmh5Dw8P9e3bVx988IGMMS6uDgDgCl27dtXAgQP1yCOPaO7cufrTn/6klJQUffbZZ2VdWrlmjCnRG/7nz58vVn9PT0/5+vq6/Gt0QEkQuoFSdPnyZf35z39WgwYN5OPjo8jISD333HPKzs526BcZGan77rtPX375pVq1aiVfX181b95cS5Yssfd58cUX7e/ejh8/XjabrdjfpS7oO902m02jRo3S4sWL1bx5c/n5+SkmJkY7duyQJL377rtq2LChfH19dc899+T7rtS13+k+ePCggoODJUlTpkyxfxTvxRdfLLKujIwMPf3004qMjJSPj4/q1KmjwYMH69SpUzp37pyqVKmip556Kt9yR48elaenpxITE516rKLs3btXffv2VbVq1eTr66s2bdro888/L3KZPMuWLVN0dLT8/f3tbXv27JGfn58GDx7s0Pc///mPPD09NWHCBIf2//mf/9GhQ4e0fft2p9YJAChbf/jDHyRJKSkpDu3Z2dkaM2aMgoODVaVKFfXq1cv+pnmezz77TN27d1dYWJh8fHzUoEED/fnPf9aVK1cc+u3fv199+vRRaGiofH19VadOHT344IPKzMx06Pfhhx+qdevW8vPzU7Vq1fTggw/qyJEjNxxD3n7B3r171a9fPwUEBKh69ep66qmndPHiRYe+xd2nWbVqldq0aSM/Pz+9++67stlsOn/+vObNm2ffPxgyZIhDHbt379ZDDz2kW265RXfeeack6YcfftCQIUNUv359+fr6KjQ0VEOHDtXp06cd1nv9d7ojIyO1a9cuff311/b13XPPPTfcJoAreJV1AYC7y8zMLDDA5eTk5Gt75JFHNG/ePPXt21djx47Vpk2blJiYqD179mjp0qUOfffv36/+/fvrscceU3x8vObMmaMHHnhAK1eu1P/8z/+od+/eCgoK0tNPP624uDh169bNIeTdjA0bNujzzz/XyJEjJUmJiYm677779Mwzz2jmzJl64okn9PPPP+vVV1/V0KFDtXbt2gIfJzg4WLNmzdLjjz+uXr16qXfv3pKk2267rdB1nzt3TnfddZf27NmjoUOH6ne/+51OnTqlzz//XEePHlWrVq3Uq1cvLVq0SK+//ro8PT3ty/7zn/+UMUYDBgxw6rFq1KhRYA27du3SHXfcodq1a+vZZ59VlSpV9PHHH6tnz5769NNP1atXr0Lrz8nJ0ZYtW/T44487tDdr1kx//vOfNX78ePXt21f333+/zp8/ryFDhqhp06aaOnWqQ//WrVtLkr755hvdfvvtha4PAFA+HDhwQJJUvXp1h/bRo0frlltuUUJCgg4ePKjp06dr1KhRWrRokb3P3Llz5e/vrzFjxsjf319r167V5MmTlZWVpWnTpkmSLl26pM6dOys7O1ujR49WaGiojh07puXLlysjI0OBgYGSpJdfflkvvPCC+vXrp0ceeUQnT57U22+/rbvvvlvbtm1z6tNx/fr1U2RkpBITE/Xtt9/qrbfe0s8//6wPPvjA3qc4+zT79u1TXFycHn30UQ0fPlxNmjTR/Pnz9cgjj6ht27YaMWKEJKlBgwYOyz3wwANq1KiR/vKXv9g/+bV69Wr99NNPevjhhxUaGqpdu3bpvffe065du/Ttt98WemR7+vTpGj16tPz9/fX8889LkkJCQm64LQCXMABKZM6cOUZSkbdbb73V3n/79u1GknnkkUccHmfcuHFGklm7dq29rW7dukaS+fTTT+1tmZmZplatWub222+3t6WkpBhJZtq0aTesN6/vnDlz7G0JCQnm+pcBScbHx8ekpKTY2959910jyYSGhpqsrCx7+8SJE40kh77x8fGmbt269p9PnjxpJJmEhIQb1miMMZMnTzaSzJIlS/Ldl5uba4wxZtWqVUaS+de//uVw/2233WY6dOhQrMcqaLt07NjRREVFmYsXLzr0b9++vWnUqFGR9f/4449Gknn77bfz3XflyhVz5513mpCQEHPq1CkzcuRI4+XlZbZs2VLgY3l7e5vHH3+8yPUBAEpX3vz/1VdfmZMnT5ojR46YhQsXmurVqxs/Pz9z9OhRh36xsbH2OccYY55++mnj6elpMjIy7G0XLlzIt55HH33UVK5c2T4Xbdu2zUgyixcvLrS2gwcPGk9PT/Pyyy87tO/YscN4eXnla79e3n7B/fff79D+xBNPGEnm+++/N8aUbJ9m5cqV+dZXpUoVEx8fX2gdcXFx+e4raFv985//NJLMv//9b3tb3va/dh/l1ltvddhPAEoLHy8HbtKMGTO0evXqfLfrj+auWLFCkjRmzBiH9rFjx0qSvvjiC4f2sLAwhyOqAQEBGjx4sLZt22b5Wa07duzo8FH16OhoSVKfPn1UtWrVfO0//fSTy9b96aefqmXLlgUeTc579zo2NlZhYWFasGCB/b6dO3fqhx9+0MCBA4v1WNc7c+aM1q5dq379+uns2bM6deqUTp06pdOnT6tz587av3+/jh07Vmj9eR9vu+WWW/Ld5+Hhoblz5+rcuXPq2rWrZs6cqYkTJxZ6xvlbbrnlhh+DBwCUjdjYWAUHBys8PFwPPvig/P39tXTpUtWuXduh34gRIxzmnLvuuktXrlzRoUOH7G1+fn72/+fNPXfddZcuXLigvXv3SpL9SPaqVasKPVnbkiVLlJubq379+tnnr1OnTik0NFSNGjXSunXrnBpb3ifd8owePVrSr/syxd2nqVevnjp37uzUuq9V0AlJr91WFy9e1KlTp9SuXTtJ0tatW4u9DqA08PFy4Ca1bdu2wNB0fWA6dOiQPDw81LBhQ4d+oaGhCgoKcph8Jalhw4b5gmHjxo0lXf2udGhoqKuGkE9ERITDz3kTfXh4eIHtP//8s8vWfeDAAfXp06fIPh4eHhowYIBmzZqlCxcuqHLlylqwYIF8fX31wAMPFOuxrvfjjz/KGKMXXnhBL7zwQoF90tPT8+1UXc8UcgK0Bg0a6MUXX9T48ePVokWLQteR9xicAAYAyqcZM2aocePG8vLyUkhIiJo0aSIPj/zHs66fU/PelL127ty1a5cmTZqktWvXKisry6F/3ve169WrpzFjxuj111/XggULdNddd+n+++/XwIED7fPx/v37ZYxRo0aNCqy5UqVKTo3t+uUbNGggDw8P+/eji7tPU69ePafWe72Cljtz5oymTJmihQsXKj093eG+67/bDpQXhG6glLlDiLr2e9LOtBcWMK00ePBgTZs2TcuWLVNcXJw++ugj3XffffYdj5LKzc2VJI0bN67Qd+Wv38m4Vt53+Yp6I+LLL7+UJB0/flynT58u9A2UjIyMQr93DgAoW4W96X69G82dGRkZ6tChgwICAjR16lQ1aNBAvr6+2rp1qyZMmGCflyTptdde05AhQ/TZZ5/pyy+/1JNPPmn/3nWdOnWUm5srm82mf/3rXwWut6Tnfils38XZfZprj04XR0HL9evXTxs3btT48ePVqlUr+fv7Kzc3V126dHHYVkB5QugGSkndunWVm5ur/fv3q1mzZvb2tLQ0ZWRkOFxHUvr1iOu1E9p///tfSSr2WcrLUnHfZGjQoIF27tx5w34tWrTQ7bffrgULFqhOnTo6fPiw3n777RI91rXq168v6erRgNjY2GItK109ouHn55fv7LV5Zs+erdWrV+vll19WYmKiHn300QIvL3Ps2DFdunTJ4XcFAFDxrF+/XqdPn9aSJUt0991329sLm0eioqIUFRWlSZMmaePGjbrjjjs0e/ZsvfTSS2rQoIGMMapXr57903ElsX//foejzD/++KNyc3Pt+x/F3acpTHH3EX7++WetWbNGU6ZM0eTJkx3qtWJ9gKvwnW6glHTr1k3S1bNnXuv111+XJHXv3t2h/fjx4w5n/8zKytIHH3ygVq1aWfrRclerXLmypKvv5DujT58++v777/Od+VTKf0R90KBB+vLLLzV9+nRVr15dXbt2LfFj5alZs6buuecevfvuuzpx4kS++6+/zMv1KlWqpDZt2ui7777Ld19KSorGjx+vPn366LnnntPf/vY3ff755w5ng82TnJwsSWrfvn2R6wMAuLe8I9LXzkuXLl3SzJkzHfplZWXp8uXLDm1RUVHy8PCwX6ard+/e8vT01JQpU/LNc8aYfJfVKsyMGTMcfs57Uztvni3uPk1hqlSp4vT+gVTwtiqoDletD3AVjnQDpaRly5aKj4/Xe++9Z/8o2ebNmzVv3jz17NlT9957r0P/xo0ba9iwYdqyZYtCQkL0/vvvKy0tTXPmzCmjEZSMn5+fmjdvrkWLFqlx48aqVq2aWrRooRYtWhTYf/z48frkk0/0wAMPaOjQoWrdurXOnDmjzz//XLNnz1bLli3tfR966CE988wzWrp0qR5//PF831UrzmNda8aMGbrzzjsVFRWl4cOHq379+kpLS1NSUpKOHj2q77//vsgx9+jRQ88//7yysrIUEBAg6eoOwtChQ+Xn56dZs2ZJkh599FF9+umneuqpp+wnh8uzevVqRUREcLkwAKjg2rdvr1tuuUXx8fF68sknZbPZNH/+/HzBcu3atRo1apQeeOABNW7cWJcvX9b8+fPl6elpP39JgwYN9NJLL2nixIk6ePCgevbsqapVqyolJUVLly7ViBEjNG7cuBvWlJKSovvvv19dunRRUlKSPvzwQz300EP2ebO4+zSFad26tb766iu9/vrrCgsLU7169ewnaS1IQECA7r77br366qvKyclR7dq19eWXXxb6qYCC1jdr1iy99NJLatiwoWrWrGm/vjpgqdI/YTpQMeRdiqKwyz116NDB4ZJhxhiTk5NjpkyZYurVq2cqVapkwsPDzcSJEx0uTWXM1ctrdO/e3axatcrcdtttxsfHxzRt2jTfZUKsumTYyJEjnVrPunXr8l2+5PpLhhljzMaNG03r1q2Nt7e3U5cPO336tBk1apSpXbu28fb2NnXq1DHx8fHm1KlT+fp269bNSDIbN24s0WMVtF2MMebAgQNm8ODBJjQ01FSqVMnUrl3b3HfffeaTTz4psnZjjElLSzNeXl5m/vz59rY333wz32XgjDHm8OHDJiAgwHTr1s3eduXKFVOrVi0zadKkG64LAFC6bjT/36hf3ty5bt06e9s333xj2rVrZ/z8/ExYWJh55pln7JfHzOv3008/maFDh5oGDRoYX19fU61aNXPvvfear776Kt+6P/30U3PnnXeaKlWqmCpVqpimTZuakSNHmn379hVZc95+we7du03fvn1N1apVzS233GJGjRplfvnlF4e+xd2nKcjevXvN3Xffbfz8/Iwk++XD8uo4efJkvmWOHj1qevXqZYKCgkxgYKB54IEHzPHjx/PtXxR0ybDU1FTTvXt3U7VqVSOJy4eh1NiMKYMzIAEoUmRkpFq0aKHly5eXdSnlXq9evbRjxw79+OOPZV2Kg2HDhum///2vNmzYUOxlly1bpoceekgHDhxQrVq1LKgOAID8XnzxRU2ZMkUnT57kRJ6AC/GdbgBu68SJE/riiy80aNCgsi4ln4SEBG3ZskXffPNNsZf961//qlGjRhG4AQAAKgC+0w3A7aSkpOibb77RP/7xD1WqVEmPPvpoWZeUT0REhC5evFiiZZOSklxcDQAAAMoKR7oBuJ2vv/5agwYNUkpKiubNm+dWZ3MHAADAbwvf6QYAAAAAwCIc6QYAAAAAwCKEbgAAAAAALMKJ1FwgNzdXx48fV9WqVWWz2cq6HACAmzPG6OzZswoLC5OHB++PuwrzNQDAlZydrwndLnD8+HGFh4eXdRkAgArmyJEjqlOnTlmXUWEwXwMArHCj+ZrQ7QJVq1aVdHVjBwQElHE1AAB3l5WVpfDwcPv8AtdgvgYAuJKz8zWh2wXyPqIWEBDAJA4AcBk+Au1azNcAACvcaL7mi2IAAAAAAFiE0A0AAAAAgEUI3QAAAAAAWITQDQAAAACARQjdAAAAAABYhNANAAAAAIBFCN0AAAAAAFiE0A0AAAAAgEUI3QAAAAAAWITQDQAAAACARQjdAAAAAABYhNANAAAAAIBFCN0AAAAAAFiE0A0AAAAAgEUI3QAAAAAAWITQDQAAAACARQjdAAAAAABYhNANAAAAAIBFCN0AAAAAAFiE0A0AAAAAgEUI3QAAAAAAWITQDQAAAACARQjdAAAAAABYhNANAAAAAIBFCN0AAAAAAFiE0A0AAAAAgEUI3QAAAAAAWITQDQAAAACARQjdAAAAAABYhNANAAAAAIBFCN0AAAAAAFiE0A0AAAAAgEUI3QAAAAAAWITQDQAAAACARdwudM+YMUORkZHy9fVVdHS0Nm/eXGT/xYsXq2nTpvL19VVUVJRWrFhRaN/HHntMNptN06dPd3HVAAD89jBnAwDgZqF70aJFGjNmjBISErR161a1bNlSnTt3Vnp6eoH9N27cqLi4OA0bNkzbtm1Tz5491bNnT+3cuTNf36VLl+rbb79VWFiY1cMAAKDCY84GAOAqtwrdr7/+uoYPH66HH35YzZs31+zZs1W5cmW9//77BfZ/88031aVLF40fP17NmjXTn//8Z/3ud7/TO++849Dv2LFjGj16tBYsWKBKlSqVxlAAAKjQmLMBALjKbUL3pUuXlJycrNjYWHubh4eHYmNjlZSUVOAySUlJDv0lqXPnzg79c3NzNWjQII0fP1633nqrNcUDAPAbwpwNAMCvvMq6AGedOnVKV65cUUhIiEN7SEiI9u7dW+AyqampBfZPTU21//zXv/5VXl5eevLJJ52uJTs7W9nZ2fafs7KynF4WAICKrrzM2czXAIDywG2OdFshOTlZb775pubOnSubzeb0comJiQoMDLTfwsPDLawSAACUZM5mvgYAlAduE7pr1KghT09PpaWlObSnpaUpNDS0wGVCQ0OL7L9hwwalp6crIiJCXl5e8vLy0qFDhzR27FhFRkYWWsvEiROVmZlpvx05cuTmBgcAQAVSXuZs5msAQHngNqHb29tbrVu31po1a+xtubm5WrNmjWJiYgpcJiYmxqG/JK1evdref9CgQfrhhx+0fft2+y0sLEzjx4/XqlWrCq3Fx8dHAQEBDjcAAHBVeZmzma8BAOWB23ynW5LGjBmj+Ph4tWnTRm3bttX06dN1/vx5Pfzww5KkwYMHq3bt2kpMTJQkPfXUU+rQoYNee+01de/eXQsXLtR3332n9957T5JUvXp1Va9e3WEdlSpVUmhoqJo0aVK6gwMAoAJhzgYA4Cq3Ct39+/fXyZMnNXnyZKWmpqpVq1ZauXKl/cQrhw8flofHrwfv27dvr48++kiTJk3Sc889p0aNGmnZsmVq0aJFWQ0BAIDfBOZsAACushljTFkX4e6ysrIUGBiozMxMProGALhpzCvWYLsCAFzJ2XnFbb7TDQAAAACAuyF0AwAAAABgEUI3AAAAAAAWIXQDAAAAAGARQjcAAAAAABYhdAMAAAAAYBFCNwAAAAAAFiF0AwAAAABgEUI3AAAAAAAWIXQDAAAAAGARQjcAAAAAABYhdAMAAAAAYBFCNwAAAAAAFiF0AwAAAABgEUI3AAAAAAAWIXQDAAAAAGARQjcAAAAAABYhdAMAAAAAYBFCNwAAAAAAFiF0AwAAAABgEUI3AAAAAAAWIXQDAAAAAGARQjcAAAAAABYhdAMAAAAAYBFCNwAAAAAAFiF0AwAAAABgEUI3AAAAAAAWIXQDAAAAAGARQjcAAAAAABYhdAMAAAAAYBFCNwAAAAAAFiF0AwAAAABgEUI3AAAAAAAWIXQDAAAAAGARQjcAAAAAABYhdAMAAAAAYBFCNwAAAAAAFiF0AwAAAABgEUI3AAAAAAAWIXQDAAAAAGARQjcAAAAAABYhdAMAAAAAYBFCNwAAAAAAFiF0AwAAAABgEUI3AAAAAAAWcbvQPWPGDEVGRsrX11fR0dHavHlzkf0XL16spk2bytfXV1FRUVqxYoX9vpycHE2YMEFRUVGqUqWKwsLCNHjwYB0/ftzqYQAAUOExZwMA4Gahe9GiRRozZowSEhK0detWtWzZUp07d1Z6enqB/Tdu3Ki4uDgNGzZM27ZtU8+ePdWzZ0/t3LlTknThwgVt3bpVL7zwgrZu3aolS5Zo3759uv/++0tzWAAAVDjM2QAAXGUzxpiyLsJZ0dHR+v3vf6933nlHkpSbm6vw8HCNHj1azz77bL7+/fv31/nz57V8+XJ7W7t27dSqVSvNnj27wHVs2bJFbdu21aFDhxQREeFUXVlZWQoMDFRmZqYCAgJKMDIAAH5VEeaV8jhnV4TtCgAoP5ydV9zmSPelS5eUnJys2NhYe5uHh4diY2OVlJRU4DJJSUkO/SWpc+fOhfaXpMzMTNlsNgUFBRXaJzs7W1lZWQ43AABwVXmZs5mvAQDlgduE7lOnTunKlSsKCQlxaA8JCVFqamqBy6Smphar/8WLFzVhwgTFxcUV+U5FYmKiAgMD7bfw8PBijgYAgIqrvMzZzNcAgPLAbUK31XJyctSvXz8ZYzRr1qwi+06cOFGZmZn225EjR0qpSgAA4OyczXwNACgPvMq6AGfVqFFDnp6eSktLc2hPS0tTaGhogcuEhoY61T9v8j506JDWrl17w+95+fj4yMfHpwSjAACg4isvczbzNQCgPHCbI93e3t5q3bq11qxZY2/Lzc3VmjVrFBMTU+AyMTExDv0lafXq1Q798ybv/fv366uvvlL16tWtGQAAAL8RzNkAAPzKbY50S9KYMWMUHx+vNm3aqG3btpo+fbrOnz+vhx9+WJI0ePBg1a5dW4mJiZKkp556Sh06dNBrr72m7t27a+HChfruu+/03nvvSbo6efft21dbt27V8uXLdeXKFft3x6pVqyZvb++yGSgAAG6OORsAgKvcKnT3799fJ0+e1OTJk5WamqpWrVpp5cqV9hOvHD58WB4evx68b9++vT766CNNmjRJzz33nBo1aqRly5apRYsWkqRjx47p888/lyS1atXKYV3r1q3TPffcUyrjAgCgomHOBgDgKre6Tnd5xXU/AQCuxLxiDbYrAMCVKtx1ugEAAAAAcDeEbgAAAAAALELoBgAAAADAIoRuAAAAAAAsQugGAAAAAMAihG4AAAAAACxC6AYAAAAAwCKEbgAAAAAALELoBgAAAADAIoRuAAAAAAAsQugGAAAAAMAihG4AAAAAACxC6AYAAAAAwCKEbgAAAAAALELoBgAAAADAIoRuAAAAAAAsQugGAAAAAMAihG4AAAAAACxC6AYAAAAAwCKEbgAAAAAALELoBgAAAADAIoRuAAAAAAAsQugGAAAAAMAihG4AAAAAACxC6AYAAAAAwCKEbgAAAAAALELoBgAAAADAIoRuAAAAAAAsQugGAAAAAMAihG4AAAAAACxC6AYAAAAAwCJezna8/fbbZbPZnOq7devWEhcEAAAAAEBF4XTo7tmzp4VlAAAAAABQ8TgduhMSEqysAwAAAACACofvdAMAAAAAYBGnj3Rf68qVK3rjjTf08ccf6/Dhw7p06ZLD/WfOnHFJcQAAAAAAuLMSHemeMmWKXn/9dfXv31+ZmZkaM2aMevfuLQ8PD7344osuLhEAAAAAAPdUotC9YMEC/f3vf9fYsWPl5eWluLg4/eMf/9DkyZP17bffurpGAAAAAADcUolCd2pqqqKioiRJ/v7+yszMlCTdd999+uKLL1xXHQAAAAAAbqxEobtOnTo6ceKEJKlBgwb68ssvJUlbtmyRj4+P66oDAAAAAMCNlSh09+rVS2vWrJEkjR49Wi+88IIaNWqkwYMHa+jQoS4tEAAAAAAAd1Wis5e/8sor9v/3799fERERSkpKUqNGjfTHP/7RZcUBAAAAAODOShS6rxcTE6OYmBhXPBQAAAAAABVGiUP3/v37tW7dOqWnpys3N9fhvsmTJ990YQAAAAAAuLsSfaf773//u5o1a6bJkyfrk08+0dKlS+23ZcuWubhERzNmzFBkZKR8fX0VHR2tzZs3F9l/8eLFatq0qXx9fRUVFaUVK1Y43G+M0eTJk1WrVi35+fkpNjZW+/fvt3IIAAD8JjBnAwBQwtD90ksv6eWXX1Zqaqq2b9+ubdu22W9bt251dY12ixYt0pgxY5SQkKCtW7eqZcuW6ty5s9LT0wvsv3HjRsXFxWnYsGHatm2bevbsqZ49e2rnzp32Pq+++qreeustzZ49W5s2bVKVKlXUuXNnXbx40bJxAABQ0TFnAwBwlc0YY4q7UEBAgLZv36769etbUVOhoqOj9fvf/17vvPOOJCk3N1fh4eEaPXq0nn322Xz9+/fvr/Pnz2v58uX2tnbt2qlVq1aaPXu2jDEKCwvT2LFjNW7cOElSZmamQkJCNHfuXD344INO1ZWVlaXAwEBlZmYqICDABSMFAPyWuWpeiY+P17Bhw3T33Xe7sDrnlMc5m/kaAOBKzs4rJfpO9wMPPKAvv/xSjz32WIkLLK5Lly4pOTlZEydOtLd5eHgoNjZWSUlJBS6TlJSkMWPGOLR17tzZ/hH4lJQUpaamKjY21n5/YGCgoqOjlZSU5HTodoWLOVf0ztofS219AABrPH5PA1Xxccl5Sm9aZmamYmNjVbduXT388MOKj49X7dq1LV9vRZ+zAQAojhLtFTRs2FAvvPCCvv32W0VFRalSpUoO9z/55JMuKe5ap06d0pUrVxQSEuLQHhISor179xa4TGpqaoH9U1NT7ffntRXWpyDZ2dnKzs62/5yVleX8QAp7zMu5emcdoRsA3F18+8hyE7qXLVumkydPav78+Zo3b54SEhIUGxurYcOGqUePHvnmb1cpL3O2FfM1AADFVaK9gvfee0/+/v76+uuv9fXXXzvcZ7PZLAnd5UliYqKmTJni0sf09vTQkPaRLn1MAEDp861UotOlWCY4OFhjxozRmDFjtHXrVs2ZM0eDBg2Sv7+/Bg4cqCeeeEKNGjUq6zItYcV8DQBAcZUodKekpLi6jhuqUaOGPD09lZaW5tCelpam0NDQApcJDQ0tsn/ev2lpaapVq5ZDn1atWhVay8SJEx0+ApeVlaXw8PBijed6ft6eevH+W2/qMQAAKMyJEye0evVqrV69Wp6enurWrZt27Nih5s2b69VXX9XTTz/tsnWVlznbivkaAIDiKl9vxxfB29tbrVu31po1a+xtubm5WrNmjWJiYgpcJiYmxqG/JK1evdrev169egoNDXXok5WVpU2bNhX6mJLk4+OjgIAAhxsAAOVNTk6OPv30U913332qW7euFi9erD/96U86fvy45s2bp6+++koff/yxpk6d6tL1lpc5m/kaAFAelOhI9/UnOsljs9nk6+urhg0bqkePHqpWrdpNFVfQeuPj49WmTRu1bdtW06dP1/nz5/Xwww9LkgYPHqzatWsrMTFRkvTUU0+pQ4cOeu2119S9e3ctXLhQ3333nd577z17vX/605/00ksvqVGjRqpXr55eeOEFhYWFqWfPni6tHQCA0larVi3l5uYqLi5OmzdvLvCI8L333qugoCCXr5s5GwCAq0oUuvOux33lyhU1adJEkvTf//5Xnp6eatq0qWbOnKmxY8fqP//5j5o3b+6yYvv376+TJ09q8uTJSk1NVatWrbRy5Ur7SVUOHz4sD49fD963b99eH330kSZNmqTnnntOjRo10rJly9SiRQt7n2eeeUbnz5/XiBEjlJGRoTvvvFMrV66Ur6+vy+oGAKAsvPHGG3rggQeKnNOCgoIs+doYczYAAFeV6Drd06dP14YNGzRnzhz7R7UyMzP1yCOP6M4779Tw4cP10EMP6ZdfftGqVatcXnR5w3U/AQCuxLxiDbYrAMCVnJ1XShS6a9eurdWrV+c7ir1r1y516tRJx44d09atW9WpUyedOnWq+NW7GSZxAIArMa9Yg+0KAHAlZ+eVEp1ILTMzU+np6fnaT548ab8GZlBQkC5dulSShwcAAAAAoEIoUeju0aOHhg4dqqVLl+ro0aM6evSoli5dqmHDhtlPZrJ582Y1btzYlbUCAAAAAOBWSnQitXfffVdPP/20HnzwQV2+fPnqA3l5KT4+Xm+88YYkqWnTpvrHP/7hukoBAAAAAHAzJfpOd55z587pp59+kiTVr19f/v7+LivMnfAdMQCAKzGvWIPtCgBwJWfnlRId6c7j7++v22677WYeAgAAAACACsvp0N27d2/NnTtXAQEB6t27d5F9lyxZctOFAQAAAADg7pwO3YGBgbLZbPb/AwAAAACAojkduufMmWP//8yZM5Wbm6sqVapIkg4ePKhly5apWbNm6ty5s+urBAAAAADADZX4kmHz58+XJGVkZKhdu3Z67bXX1LNnT82aNculBQIAAAAA4K5KFLq3bt2qu+66S5L0ySefKCQkRIcOHdIHH3ygt956y6UFAgAAAADgrkoUui9cuKCqVatKkr788kv17t1bHh4eateunQ4dOuTSAgEAAAAAcFclCt0NGzbUsmXLdOTIEa1atUqdOnWSJKWnp3PdSwAAAAAA/r8She7Jkydr3LhxioyMVHR0tGJiYiRdPep9++23u7RAAAAAAADcldNnL79W3759deedd+rEiRNq2bKlvb1jx47q1auXy4oDAAAAAMCdlSh0S1JoaKhCQ0Md2tq2bXvTBQEAAAAAUFGU6OPlAAAAAADgxgjdAAAAAABYhNANAAAAAIBFCN0AAAAAAFiE0A0AAAAAgEUI3QAAAAAAWITQDQAAAACARQjdAAAAAABYhNANAAAAAIBFCN0AAAAAAFiE0A0AAAAAgEUI3QAAAAAAWITQDQAAAACARQjdAAAAAABYhNANAAAAAIBFCN0AAAAAAFiE0A0AAAAAgEUI3QAAAAAAWITQDQAAAACARQjdAAAAAABYhNANAAAAAIBFCN0AAAAAAFiE0A0AAAAAgEUI3QAAAAAAWITQDQAAAACARQjdAAAAAABYhNANAAAAAIBFCN0AAAAAAFjEbUL3mTNnNGDAAAUEBCgoKEjDhg3TuXPnilzm4sWLGjlypKpXry5/f3/16dNHaWlp9vu///57xcXFKTw8XH5+fmrWrJnefPNNq4cCAECFxpwNAMCv3CZ0DxgwQLt27dLq1au1fPly/fvf/9aIESOKXObpp5/W//3f/2nx4sX6+uuvdfz4cfXu3dt+f3JysmrWrKkPP/xQu3bt0vPPP6+JEyfqnXfesXo4AABUWMzZAAD8ymaMMWVdxI3s2bNHzZs315YtW9SmTRtJ0sqVK9WtWzcdPXpUYWFh+ZbJzMxUcHCwPvroI/Xt21eStHfvXjVr1kxJSUlq165dgesaOXKk9uzZo7Vr1zpdX1ZWlgIDA5WZmamAgIASjBAAgF+587xSnudsd96uAIDyx9l5xS2OdCclJSkoKMg+eUtSbGysPDw8tGnTpgKXSU5OVk5OjmJjY+1tTZs2VUREhJKSkgpdV2ZmpqpVq+a64gEA+A1hzgYAwJFXWRfgjNTUVNWsWdOhzcvLS9WqVVNqamqhy3h7eysoKMihPSQkpNBlNm7cqEWLFumLL74osp7s7GxlZ2fbf87KynJiFAAAVHzlac5mvgYAlAdleqT72Weflc1mK/K2d+/eUqll586d6tGjhxISEtSpU6ci+yYmJiowMNB+Cw8PL5UaAQAoK+44ZzNfAwDKgzI90j127FgNGTKkyD7169dXaGio0tPTHdovX76sM2fOKDQ0tMDlQkNDdenSJWVkZDi8c56WlpZvmd27d6tjx44aMWKEJk2adMO6J06cqDFjxth/zsrKYiIHAFRo7jhnM18DAMqDMg3dwcHBCg4OvmG/mJgYZWRkKDk5Wa1bt5YkrV27Vrm5uYqOji5wmdatW6tSpUpas2aN+vTpI0nat2+fDh8+rJiYGHu/Xbt26Q9/+IPi4+P18ssvO1W3j4+PfHx8nOoLAEBF4I5zNvM1AKA8cIuzl0tS165dlZaWptmzZysnJ0cPP/yw2rRpo48++kiSdOzYMXXs2FEffPCB2rZtK0l6/PHHtWLFCs2dO1cBAQEaPXq0pKvfA5OufjztD3/4gzp37qxp06bZ1+Xp6enUjkUezoYKAHAld59Xyuuc7e7bFQBQvjg7r7jFidQkacGCBRo1apQ6duwoDw8P9enTR2+99Zb9/pycHO3bt08XLlywt73xxhv2vtnZ2ercubNmzpxpv/+TTz7RyZMn9eGHH+rDDz+0t9etW1cHDx4slXEBAFDRMGcDAPArtznSXZ7xzjkAwJWYV6zBdgUAuFKFuk43AAAAAADuiNANAAAAAIBFCN0AAAAAAFiE0A0AAAAAgEUI3QAAAAAAWITQDQAAAACARQjdAAAAAABYhNANAAAAAIBFCN0AAAAAAFiE0A0AAAAAgEUI3QAAAAAAWITQDQAAAACARQjdAAAAAABYhNANAAAAAIBFCN0AAAAAAFiE0A0AAAAAgEUI3QAAAAAAWITQDQAAAACARQjdAAAAAABYhNANAAAAAIBFCN0AAAAAAFiE0A0AAAAAgEUI3QAAAAAAWITQDQAAAACARQjdAAAAAABYhNANAAAAAIBFCN0AAAAAAFiE0A0AAAAAgEUI3QAAAAAAWITQDQAAAACARQjdAAAAAABYhNANAAAAAIBFCN0AAAAAAFiE0A0AAAAAgEUI3QAAAAAAWITQDQAAAACARQjdAAAAAABYhNANAAAAAIBFCN0AAAAAAFiE0A0AAAAAgEUI3QAAAAAAWITQDQAAAACARQjdAAAAAABYhNANAAAAAIBFCN0AAAAAAFjEbUL3mTNnNGDAAAUEBCgoKEjDhg3TuXPnilzm4sWLGjlypKpXry5/f3/16dNHaWlpBfY9ffq06tSpI5vNpoyMDAtGAADAbwNzNgAAv3Kb0D1gwADt2rVLq1ev1vLly/Xvf/9bI0aMKHKZp59+Wv/3f/+nxYsX6+uvv9bx48fVu3fvAvsOGzZMt912mxWlAwDwm8KcDQDAr2zGGFPWRdzInj171Lx5c23ZskVt2rSRJK1cuVLdunXT0aNHFRYWlm+ZzMxMBQcH66OPPlLfvn0lSXv37lWzZs2UlJSkdu3a2fvOmjVLixYt0uTJk9WxY0f9/PPPCgoKcrq+rKwsBQYGKjMzUwEBATc3WADAb547zyvlec525+0KACh/nJ1X3OJId1JSkoKCguyTtyTFxsbKw8NDmzZtKnCZ5ORk5eTkKDY21t7WtGlTRUREKCkpyd62e/duTZ06VR988IE8PJzbHNnZ2crKynK4AQCA8jVnM18DAMoDtwjdqampqlmzpkObl5eXqlWrptTU1EKX8fb2zvfud0hIiH2Z7OxsxcXFadq0aYqIiHC6nsTERAUGBtpv4eHhxRsQAAAVVHmas5mvAQDlQZmG7meffVY2m63I2969ey1b/8SJE9WsWTMNHDiw2MtlZmbab0eOHLGoQgAAygd3nLOZrwEA5YFXWa587NixGjJkSJF96tevr9DQUKWnpzu0X758WWfOnFFoaGiBy4WGhurSpUvKyMhweOc8LS3NvszatWu1Y8cOffLJJ5KkvK+316hRQ88//7ymTJlS4GP7+PjIx8fHmSECAFAhuOOczXwNACgPyjR0BwcHKzg4+Ib9YmJilJGRoeTkZLVu3VrS1ck3NzdX0dHRBS7TunVrVapUSWvWrFGfPn0kSfv27dPhw4cVExMjSfr000/1yy+/2JfZsmWLhg4dqg0bNqhBgwY3OzwAACoM5mwAAEqmTEO3s5o1a6YuXbpo+PDhmj17tnJycjRq1Cg9+OCD9rOgHjt2TB07dtQHH3ygtm3bKjAwUMOGDdOYMWNUrVo1BQQEaPTo0YqJibGfBfX6SfrUqVP29RXn7OUAAOAq5mwAABy5ReiWpAULFmjUqFHq2LGjPDw81KdPH7311lv2+3NycrRv3z5duHDB3vbGG2/Y+2ZnZ6tz586aOXNmWZQPAMBvBnM2AAC/covrdJd3XPcTAOBKzCvWYLsCAFypQl2nGwAAAAAAd0ToBgAAAADAIoRuAAAAAAAsQugGAAAAAMAihG4AAAAAACxC6AYAAAAAwCKEbgAAAAAALELoBgAAAADAIoRuAAAAAAAsQugGAAAAAMAihG4AAAAAACxC6AYAAAAAwCKEbgAAAAAALELoBgAAAADAIoRuAAAAAAAsQugGAAAAAMAihG4AAAAAACxC6AYAAAAAwCKEbgAAAAAALELoBgAAAADAIoRuAAAAAAAsQugGAAAAAMAihG4AAAAAACxC6AYAAAAAwCKEbgAAAAAALELoBgAAAADAIoRuAAAAAAAsQugGAAAAAMAihG4AAAAAACxC6AYAAAAAwCKEbgAAAAAALELoBgAAAADAIoRuAAAAAAAsQugGAAAAAMAiXmVdQEVgjJEkZWVllXElAICKIG8+yZtf4BrM1wAAV3J2viZ0u8DZs2clSeHh4WVcCQCgIjl79qwCAwPLuowKg/kaAGCFG83XNsPb6DctNzdXx48fV9WqVWWz2SxfX1ZWlsLDw3XkyBEFBARYvj4rVaSxSBVrPIylfKpIY5Eq1nhcORZjjM6ePauwsDB5ePBNMFcp7fm6NFWkvyUrsZ2cw3ZyDtvJORV5Ozk7X3Ok2wU8PDxUp06dUl9vQEBAhfnFrUhjkSrWeBhL+VSRxiJVrPG4aiwc4Xa9spqvS1NF+luyEtvJOWwn57CdnFNRt5Mz8zVvnwMAAAAAYBFCNwAAAAAAFiF0uyEfHx8lJCTIx8enrEu5aRVpLFLFGg9jKZ8q0likijWeijQWuB9+/5zDdnIO28k5bCfnsJ04kRoAAAAAAJbhSDcAAAAAABYhdAMAAAAAYBFCNwAAAAAAFiF0l1OJiYn6/e9/r6pVq6pmzZrq2bOn9u3b59Dnnnvukc1mc7g99thjZVRx0V588cV8tTZt2tR+/8WLFzVy5EhVr15d/v7+6tOnj9LS0sqw4sJFRkbmG4vNZtPIkSMlle/n5d///rf++Mc/KiwsTDabTcuWLXO43xijyZMnq1atWvLz81NsbKz279/v0OfMmTMaMGCAAgICFBQUpGHDhuncuXOlOIpfFTWenJwcTZgwQVFRUapSpYrCwsI0ePBgHT9+3OExCno+X3nllVIeyY2fmyFDhuSrs0uXLg59ystzc6OxFPT3Y7PZNG3aNHuf8vK8OPNa7Mzr1+HDh9W9e3dVrlxZNWvW1Pjx43X58uXSHAoqgJL8jRdnfj19+rTq1Kkjm82mjIwMC0ZQOqzYTt9//73i4uIUHh4uPz8/NWvWTG+++abVQ3GpGTNmKDIyUr6+voqOjtbmzZuL7L948WI1bdpUvr6+ioqK0ooVKxzud2afwd24chs5ux/ijlz9u3Stxx57TDabTdOnT3dx1WXMoFzq3LmzmTNnjtm5c6fZvn276datm4mIiDDnzp2z9+nQoYMZPny4OXHihP2WmZlZhlUXLiEhwdx6660OtZ48edJ+/2OPPWbCw8PNmjVrzHfffWfatWtn2rdvX4YVFy49Pd1hHKtXrzaSzLp164wx5ft5WbFihXn++efNkiVLjCSzdOlSh/tfeeUVExgYaJYtW2a+//57c//995t69eqZX375xd6nS5cupmXLlubbb781GzZsMA0bNjRxcXGlPJKrihpPRkaGiY2NNYsWLTJ79+41SUlJpm3btqZ169YOj1G3bl0zdepUh+fr2r+z0nKj5yY+Pt506dLFoc4zZ8449Ckvz82NxnLtGE6cOGHef/99Y7PZzIEDB+x9ysvz4sxr8Y1evy5fvmxatGhhYmNjzbZt28yKFStMjRo1zMSJE0t9PHBvJfkbL8782qNHD9O1a1cjyfz8888WjKB0WLGd/vd//9c8+eSTZv369ebAgQNm/vz5xs/Pz7z99ttWD8clFi5caLy9vc37779vdu3aZYYPH26CgoJMWlpagf2/+eYb4+npaV599VWze/duM2nSJFOpUiWzY8cOex9n9hnciau3kbP7Ie7Git+lPEuWLDEtW7Y0YWFh5o033rB4JKWL0O0m0tPTjSTz9ddf29s6dOhgnnrqqbIrqhgSEhJMy5YtC7wvIyPDVKpUySxevNjetmfPHiPJJCUllVKFJffUU0+ZBg0amNzcXGOM+zwv14eh3NxcExoaaqZNm2Zvy8jIMD4+Puaf//ynMcaY3bt3G0lmy5Yt9j7/+te/jM1mM8eOHSu12gtSULi73ubNm40kc+jQIXtb3bp1y90Le2Ghu0ePHoUuU16fG2eelx49epg//OEPDm3l8XkxJv9rsTOvXytWrDAeHh4mNTXV3mfWrFkmICDAZGdnl+4A4LZK8jdenPl15syZpkOHDmbNmjVuHbqt3k7XeuKJJ8y9997ruuIt1LZtWzNy5Ej7z1euXDFhYWEmMTGxwP79+vUz3bt3d2iLjo42jz76qDHGuX0Gd+PqbVSQgvZD3I1V2+no0aOmdu3aZufOneV2H+Bm8PFyN5GZmSlJqlatmkP7ggULVKNGDbVo0UITJ07UhQsXyqI8p+zfv19hYWGqX7++BgwYoMOHD0uSkpOTlZOTo9jYWHvfpk2bKiIiQklJSWVVrlMuXbqkDz/8UEOHDpXNZrO3u9PzkiclJUWpqakOz0NgYKCio6Ptz0NSUpKCgoLUpk0be5/Y2Fh5eHho06ZNpV5zcWVmZspmsykoKMih/ZVXXlH16tV1++23a9q0aeX2Y7/r169XzZo11aRJEz3++OM6ffq0/T53fW7S0tL0xRdfaNiwYfnuK4/Py/Wvxc68fiUlJSkqKkohISH2Pp07d1ZWVpZ27dpVitXDnZXkb9zZ+XX37t2aOnWqPvjgA3l4uPeuoZXb6XqZmZn59svKo0uXLik5OdlhfB4eHoqNjS10fElJSQ79pauvW3n9ndlncCdWbKOCFLYf4i6s2k65ubkaNGiQxo8fr1tvvdWa4suYV1kXgBvLzc3Vn/70J91xxx1q0aKFvf2hhx5S3bp1FRYWph9++EETJkzQvn37tGTJkjKstmDR0dGaO3eumjRpohMnTmjKlCm66667tHPnTqWmpsrb2zvfC1BISIhSU1PLpmAnLVu2TBkZGRoyZIi9zZ2el2vlbetrg0Hez3n3paamqmbNmg73e3l5qVq1auX+ubp48aImTJiguLg4BQQE2NuffPJJ/e53v1O1atW0ceNGTZw4USdOnNDrr79ehtXm16VLF/Xu3Vv16tXTgQMH9Nxzz6lr165KSkqSp6en2z438+bNU9WqVdW7d2+H9vL4vBT0WuzM61dqamqBf1d59wHOKMnfuDO/n9nZ2YqLi9O0adMUERGhn376yZL6S4tV2+l6Gzdu1KJFi/TFF1+4pG4rnTp1SleuXCnwdWjv3r0FLlPY69a1r2t5bYX1cSdWbKPrFbYf4k6s2k5//etf5eXlpSeffNL1RZcThG43MHLkSO3cuVP/+c9/HNpHjBhh/39UVJRq1aqljh076sCBA2rQoEFpl1mkrl272v9/2223KTo6WnXr1tXHH38sPz+/Mqzs5vzv//6vunbtqrCwMHubOz0vvxU5OTnq16+fjDGaNWuWw31jxoyx//+2226Tt7e3Hn30USUmJsrHx6e0Sy3Ugw8+aP9/VFSUbrvtNjVo0EDr169Xx44dy7Cym/P+++9rwIAB8vX1dWgvj89LYa/FwM149tln9de//rXIPnv27LFs/RMnTlSzZs00cOBAy9bhCmW9na61c+dO9ejRQwkJCerUqVOprBPuraj9kN+65ORkvfnmm9q6davDp0YrGvf+DNFvwKhRo7R8+XKtW7dOderUKbJvdHS0JOnHH38sjdJuSlBQkBo3bqwff/xRoaGhunTpUr4zpaalpSk0NLRsCnTCoUOH9NVXX+mRRx4psp+7PC952/r6s9pe+zyEhoYqPT3d4f7Lly/rzJkz5fa5ypvoDh06pNWrV9/w3eXo6GhdvnxZBw8eLJ0CS6h+/fqqUaOG/ffKHZ+bDRs2aN++fTf8G5LK/nkp7LXYmdev0NDQAv+u8u7Db9vYsWO1Z8+eIm/169cv0d+4M7+fa9eu1eLFi+Xl5SUvLy/7m3g1atRQQkKC6wdcQmW9nfLs3r1bHTt21IgRIzRp0iSXjtEqNWrUkKenZ5Hz+/UKe9269nUtr83ZxyzPrNhGeYq7H1KeWbGdNmzYoPT0dEVERNhfhw4dOqSxY8cqMjLSknGUBUJ3OWWM0ahRo7R06VKtXbtW9erVu+Ey27dvlyTVqlXL4upu3rlz53TgwAHVqlVLrVu3VqVKlbRmzRr7/fv27dPhw4cVExNThlUWbc6cOapZs6a6d+9eZD93eV7q1aun0NBQh+chKytLmzZtsj8PMTExysjIUHJysr3P2rVrlZuba39zoTzJm+j279+vr776StWrV7/hMtu3b5eHh0e+jyeWN0ePHtXp06ftv1fu9txIVz8p0rp1a7Vs2fKGfcvqebnRa7Ezr18xMTHasWOHQxDI2/Fq3rx56QwE5VZwcLCaNm1a5M3b27tEf+PO/H5++umn+v7777V9+3Zt375d//jHPyRd3RHOuxRmeVDW20mSdu3apXvvvVfx8fF6+eWXrRusi3l7e6t169YO48vNzdWaNWsK3c+KiYlx6C9dfd3K6+/MPoM7sWIbSSXbDynPrNhOgwYN0g8//GB/Ddq+fbvCwsI0fvx4rVq1yrrBlLayPY8bCvP444+bwMBAs379eodL5ly4cMEYY8yPP/5opk6dar777juTkpJiPvvsM1O/fn1z9913l3HlBRs7dqxZv369SUlJMd98842JjY01NWrUMOnp6caYq5fqiIiIMGvXrjXfffediYmJMTExMWVcdeGuXLliIiIizIQJExzay/vzcvbsWbNt2zazbds2I8m8/vrrZtu2bfazaL7yyismKCjIfPbZZ+aHH34wPXr0KPCSYbfffrvZtGmT+c9//mMaNWpUZpcMK2o8ly5dMvfff7+pU6eO2b59u8PfUd4Zozdu3GjeeOMNs337dnPgwAHz4YcfmuDgYDN48OByNZazZ8+acePGmaSkJJOSkmK++uor87vf/c40atTIXLx40f4Y5eW5udHvmTHGZGZmmsqVK5tZs2blW748PS83ei025savX3mXDOvUqZPZvn27WblypQkODuaSYSi2G/2NHz161DRp0sRs2rTJ3lbc+XXdunVuffZyY6zZTjt27DDBwcFm4MCBDq8Fefsx5d3ChQuNj4+PmTt3rtm9e7cZMWKECQoKsl9VYdCgQebZZ5+19//mm2+Ml5eX+dvf/mb27NljEhISCrxk2I32GdyJq7eRM/sh7siK36XrVcSzlxO6yylJBd7mzJljjDHm8OHD5u677zbVqlUzPj4+pmHDhmb8+PHl5nrQ1+vfv7+pVauW8fb2NrVr1zb9+/c3P/74o/3+X375xTzxxBPmlltuMZUrVza9evUyJ06cKMOKi7Zq1Sojyezbt8+hvbw/L3k7U9ff4uPjjTFXLwHywgsvmJCQEOPj42M6duyYb4ynT582cXFxxt/f3wQEBJiHH37YnD17tgxGU/R4UlJSCv07yrumenJysomOjjaBgYHG19fXNGvWzPzlL39xCLLlYSwXLlwwnTp1MsHBwaZSpUqmbt26Zvjw4Q6XoDKm/Dw3N/o9M8aYd9991/j5+ZmMjIx8y5en5+VGr8XGOPf6dfDgQdO1a1fj5+dnatSoYcaOHWtycnJKeTRwdzf6G8973ct7jTOm+PNrRQjdVmynhISEAl8L6tatW4ojuzlvv/22iYiIMN7e3qZt27bm22+/td/XoUMHh9doY4z5+OOPTePGjY23t7e59dZbzRdffOFwvzP7DO7GldvImf0Qd+Xq36XrVcTQbTPGGFcdNQcAAAAAAL/iO90AAAAAAFiE0A0AAAAAgEUI3QAAAAAAWITQDQAAAACARQjdAAAAAABYhNANAAAAAIBFCN0AAAAAAFiE0A0AAAAAgEUI3QBu2vr162Wz2ZSRkVHWpQAAAADlis0YY8q6CADu5Z577lGrVq00ffp0SdKlS5d05swZhYSEyGazlW1xAAAAQDniVdYFAHB/3t7eCg0NLesyAAAAgHKHj5cDKJYhQ4bo66+/1ptvvimbzSabzaa5c+c6fLx87ty5CgoK0vLly9WkSRNVrlxZffv21YULFzRv3jxFRkbqlltu0ZNPPqkrV67YHzs7O1vjxo1T7dq1VaVKFUVHR2v9+vVlM1AAAH5DTp48qdDQUP3lL3+xt23cuFHe3t5as2ZNGVYGuD+OdAMoljfffFP//e9/1aJFC02dOlWStGvXrnz9Lly4oLfeeksLFy7U2bNn1bt3b/Xq1UtBQUFasWKFfvrpJ/Xp00d33HGH+vfvL0kaNWqUdu/erYULFyosLExLly5Vly5dtGPHDjVq1KhUxwkAwG9JcHCw3n//ffXs2VOdOnVSkyZNNGjQII0aNUodO3Ys6/IAt0boBlAsgYGB8vb2VuXKle0fKd+7d2++fjk5OZo1a5YaNGggSerbt6/mz5+vtLQ0+fv7q3nz5rr33nu1bt069e/fX4cPH9acOXN0+PBhhYWFSZLGjRunlStXas6cOQ7vvAMAANfr1q2bhg8frgEDBqhNmzaqUqWKEhMTy7oswO0RugFYonLlyvbALUkhISGKjIyUv7+/Q1t6erokaceOHbpy5YoaN27s8DjZ2dmqXr166RQNAMBv3N/+9je1aNFCixcvVnJysnx8fMq6JMDtEboBWKJSpUoOP9tstgLbcnNzJUnnzp2Tp6enkpOT5enp6dDv2qAOAACsc+DAAR0/fly5ubk6ePCgoqKiyrokwO0RugEUm7e3t8MJ0Fzh9ttv15UrV5Senq677rrLpY8NAABu7NKlSxo4cKD69++vJk2a6JFHHtGOHTtUs2bNsi4NcGucvRxAsUVGRmrTpk06ePCgTp06ZT9afTMaN26sAQMGaPDgwVqyZIlSUlK0efNmJSYm6osvvnBB1QAAoCjPP/+8MjMz9dZbb2nChAlq3Lixhg4dWtZlAW6P0A2g2MaNGydPT081b95cwcHBOnz4sEsed86cORo8eLDGjh2rJk2aqGfPntqyZYsiIiJc8vgAAKBg69ev1/Tp0zV//nwFBATIw8ND8+fP14YNGzRr1qyyLg9wazZjjCnrIgAAAAAAqIg40g0AAAAAgEUI3QAAAAAAWITQDQAAAACARQjdAAAAAABYhNANAAAAAIBFCN0AAAAAAFiE0A0AAAAAgEUI3QAAAAAAWITQDQAAAACARQjdAAAAAABYhNANAAAAAIBFCN0AAAAAAFjk/wEm9NASio16+AAAAABJRU5ErkJggg==",
      "text/plain": [
       "<Figure size 1000x400 with 2 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, axes = plt.subplots(1, 2, figsize=(10, 4))\n",
    "brainmass.viz.plot_timeseries(res['x'], ts=res['ts'], ax=axes[0])\n",
    "axes[0].set_title('Hopf limit cycle (x)')\n",
    "brainmass.viz.plot_phase_portrait(res['x'], res['y'], ax=axes[1])\n",
    "axes[1].set_title('Phase portrait')\n",
    "plt.tight_layout()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3e4de3f8",
   "metadata": {},
   "source": [
    "## Try it: cross the bifurcation\n",
    "\n",
    "Sweep the bifurcation parameter `a` across zero, giving each node a small kick off the fixed point. For `a <= 0` the node decays back to silence; for `a > 0` the settled amplitude grows as $\\sqrt{a}$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "fad1c5fd",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2026-06-19T08:10:47.583360Z",
     "iopub.status.busy": "2026-06-19T08:10:47.583170Z",
     "iopub.status.idle": "2026-06-19T08:10:48.198884Z",
     "shell.execute_reply": "2026-06-19T08:10:48.197982Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "a = -0.20  ->  settled amplitude = 0.000  (sqrt(a) = 0.000)\n",
      "a = +0.00  ->  settled amplitude = 0.069  (sqrt(a) = 0.000)\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "a = +0.25  ->  settled amplitude = 0.502  (sqrt(a) = 0.500)\n",
      "a = +1.00  ->  settled amplitude = 1.001  (sqrt(a) = 1.000)\n"
     ]
    }
   ],
   "source": [
    "for a in [-0.2, 0.0, 0.25, 1.0]:\n",
    "    m = brainmass.HopfStep(in_size=1, a=a, w=0.3)\n",
    "    brainstate.nn.init_all_states(m)\n",
    "    m.x.value = m.x.value + 0.1  # small kick off the fixed point\n",
    "    r = brainmass.Simulator(m, dt=0.1 * u.ms).run(\n",
    "        300. * u.ms, monitors=['x'], transient=150. * u.ms, init_states=False)\n",
    "    amp = float(np.sqrt(np.mean(u.get_magnitude(r['x']) ** 2)) * np.sqrt(2))\n",
    "    print(f'a = {a:+.2f}  ->  settled amplitude = {amp:.3f}  '\n",
    "          f'(sqrt(a) = {np.sqrt(max(a, 0)):.3f})')"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.13.11"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}