{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "89c4679467d7",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2026-06-19T07:41:18.992251Z",
     "iopub.status.busy": "2026-06-19T07:41:18.991859Z",
     "iopub.status.idle": "2026-06-19T07:41:24.336226Z",
     "shell.execute_reply": "2026-06-19T07:41:24.335183Z"
    },
    "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 warnings\n",
    "import time\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import jax\n",
    "import jax.numpy as jnp\n",
    "import brainstate\n",
    "import braintools\n",
    "import brainunit as u\n",
    "import brainmass\n",
    "from brainmass import objectives\n",
    "from brainstate.nn import Param\n",
    "brainstate.environ.set(dt=0.1 * u.ms)\n",
    "brainstate.random.seed(0)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2612849f3317",
   "metadata": {},
   "source": [
    "# Compose a Custom Objective\n",
    "\n",
    "**Goal:** build the loss your fit actually needs — combine the\n",
    "`brainmass.objectives` builders, write your own, and plug the result into\n",
    "`brainmass.Fitter`.\n",
    "\n",
    "An *objective* in `brainmass` is a small `callable(prediction, target) -> scalar`\n",
    "that is `jit` / `grad` / `vmap`-safe and unit-aware. The builders in\n",
    "`brainmass.objectives` return such callables; you compose them or write your own\n",
    "with the same signature.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d833e36135e3",
   "metadata": {},
   "source": [
    "## The built-in toolkit\n",
    "\n",
    "Every entry in `brainmass.objectives` is a *builder* that returns a loss/score\n",
    "callable. They wrap `braintools.metric` — no metric maths is reimplemented.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "02d543e45f3c",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2026-06-19T07:41:24.339295Z",
     "iopub.status.busy": "2026-06-19T07:41:24.338837Z",
     "iopub.status.idle": "2026-06-19T07:41:27.988721Z",
     "shell.execute_reply": "2026-06-19T07:41:27.988084Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "timeseries_rmse : 1.433825135231018\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "fc_corr (score) : -0.19663743674755096\n",
      "fc_corr (loss)  : 1.1966373920440674\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "cosine_sim      : -0.025183098390698433\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "fcd (loss)      : 0.5260788798332214\n"
     ]
    }
   ],
   "source": [
    "import numpy as np\n",
    "rng = np.random.default_rng(0)\n",
    "pred = jnp.asarray(rng.standard_normal((200, 6)))    # (time, regions)\n",
    "target = jnp.asarray(rng.standard_normal((200, 6)))\n",
    "\n",
    "print(\"timeseries_rmse :\", float(objectives.timeseries_rmse()(pred, target)))\n",
    "print(\"fc_corr (score) :\", float(objectives.fc_corr()(pred, target)))\n",
    "print(\"fc_corr (loss)  :\", float(objectives.fc_corr(as_loss=True)(pred, target)))\n",
    "print(\"cosine_sim      :\", float(objectives.cosine_sim()(pred, target)))\n",
    "print(\"fcd (loss)      :\", float(objectives.fcd(as_loss=True)(pred, target)))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b7ada5ea8dfb",
   "metadata": {},
   "source": [
    "A few conventions worth knowing:\n",
    "\n",
    "- **Score vs loss.** `fc_corr` / `cosine_sim` / `fcd` return a *score* to maximise\n",
    "  by default; pass `as_loss=True` to get `1 - score`, a quantity to *minimise*.\n",
    "  `timeseries_rmse` / `fc_rmse` are already losses.\n",
    "- **Unit safety.** `timeseries_rmse` subtracts before stripping units, so mixing\n",
    "  mV and Hz raises. FC / cosine objectives are scale-invariant and operate on\n",
    "  magnitudes.\n",
    "- **Identity is exact.** Every loss returns `0` (and every score `1`) when\n",
    "  prediction equals target.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "06979b72e56d",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2026-06-19T07:41:27.990983Z",
     "iopub.status.busy": "2026-06-19T07:41:27.990703Z",
     "iopub.status.idle": "2026-06-19T07:41:29.098317Z",
     "shell.execute_reply": "2026-06-19T07:41:29.097375Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "rmse(x, x)    = 0.0\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "fc_corr(x, x) = 1.0\n"
     ]
    }
   ],
   "source": [
    "x = jnp.zeros((50, 4)) + jnp.asarray(rng.standard_normal((50, 4)))\n",
    "print(\"rmse(x, x)    =\", float(objectives.timeseries_rmse()(x, x)))\n",
    "print(\"fc_corr(x, x) =\", float(objectives.fc_corr()(x, x)))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "54ae9022176d",
   "metadata": {},
   "source": [
    "## Combine weighted objectives\n",
    "\n",
    "`objectives.combine` sums `(weight, objective)` pairs into one callable — the\n",
    "standard way to balance, say, a time-domain term against a connectivity term.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "52e6b0d6f4b5",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2026-06-19T07:41:29.101335Z",
     "iopub.status.busy": "2026-06-19T07:41:29.100932Z",
     "iopub.status.idle": "2026-06-19T07:41:29.362050Z",
     "shell.execute_reply": "2026-06-19T07:41:29.361262Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "combined loss: 2.0321438312530518\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(2.0 + 0.5) * 1.0 = 2.5\n"
     ]
    }
   ],
   "source": [
    "combined = objectives.combine(\n",
    "    (1.0, objectives.timeseries_rmse()),\n",
    "    (0.5, objectives.fc_corr(as_loss=True)),\n",
    ")\n",
    "print(\"combined loss:\", float(combined(pred, target)))\n",
    "\n",
    "# Sanity check the arithmetic on a controlled example:\n",
    "a = jnp.zeros((10, 3))\n",
    "check = objectives.combine(\n",
    "    (2.0, objectives.timeseries_rmse()),\n",
    "    (0.5, objectives.timeseries_rmse()),\n",
    ")\n",
    "print(\"(2.0 + 0.5) * 1.0 =\", float(check(a + 1.0, a)))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6c67323a8903",
   "metadata": {},
   "source": [
    "## Write your own objective\n",
    "\n",
    "A custom objective is any function with the `(prediction, target) -> scalar`\n",
    "signature, built from `jax.numpy` so it stays differentiable. Two rules:\n",
    "\n",
    "1. Reduce to a **scalar** (a gradient optimiser needs one number).\n",
    "2. Strip units with `u.get_magnitude` only where the quantity is scale-free; keep\n",
    "   them through a subtraction so unit mismatches are caught.\n",
    "\n",
    "Here is a **spectral-peak** objective: match the dominant oscillation frequency\n",
    "of each region. It reuses `braintools.metric.power_spectral_density` rather than\n",
    "hand-rolling an FFT. (`power_spectral_density(signal, dt_ms)` returns frequencies\n",
    "in cycles per `dt` unit; comparing two peaks in the *same* units is all we need.)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "593ad7b075af",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2026-06-19T07:41:29.364277Z",
     "iopub.status.busy": "2026-06-19T07:41:29.364056Z",
     "iopub.status.idle": "2026-06-19T07:41:30.640005Z",
     "shell.execute_reply": "2026-06-19T07:41:30.639296Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "peak-freq loss (pred vs fast): 6.103515625e-05\n",
      "peak-freq loss (pred vs pred): 0.0\n"
     ]
    }
   ],
   "source": [
    "def peak_frequency_mse(dt_ms=1.0):\n",
    "    \"\"\"Loss matching the dominant spectral peak of each region.\"\"\"\n",
    "    def loss(prediction, target):\n",
    "        def peak_freqs(sig):\n",
    "            sig = u.get_magnitude(sig)\n",
    "            def one(col):\n",
    "                freqs, psd = braintools.metric.power_spectral_density(col, dt_ms)\n",
    "                return freqs[jnp.argmax(psd)]\n",
    "            return jax.vmap(one, in_axes=1)(sig)  # (regions,)\n",
    "        diff = peak_freqs(prediction) - peak_freqs(target)\n",
    "        return jnp.mean(diff ** 2)\n",
    "    return loss\n",
    "\n",
    "# A 4-region bank of sinusoids at well-separated frequencies, sampled at 1 ms.\n",
    "samples = np.arange(1000)\n",
    "dt_s = 1e-3  # 1 ms sampling\n",
    "rng2 = np.random.default_rng(1)\n",
    "def sine_bank(freqs_hz):\n",
    "    cols = [np.sin(2 * np.pi * f * samples * dt_s) + 0.05 * rng2.standard_normal(1000)\n",
    "            for f in freqs_hz]\n",
    "    return jnp.asarray(np.stack(cols, axis=1))\n",
    "\n",
    "pred_sig = sine_bank([5, 10, 20, 30])\n",
    "fast_sig = sine_bank([8, 16, 28, 40])\n",
    "\n",
    "spec_loss = peak_frequency_mse(dt_ms=1.0)\n",
    "print(\"peak-freq loss (pred vs fast):\", float(spec_loss(pred_sig, fast_sig)))\n",
    "print(\"peak-freq loss (pred vs pred):\", float(spec_loss(pred_sig, pred_sig)))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0b6caa4c2e6d",
   "metadata": {},
   "source": [
    "Custom objectives compose with the built-ins through `combine`, exactly like\n",
    "the library ones — they share the same signature. Here a connectivity term and\n",
    "the spectral term are balanced into one loss (both inputs have ≥2 regions, so\n",
    "`fc_corr` is well-defined).\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "8bef94509cf7",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2026-06-19T07:41:30.642250Z",
     "iopub.status.busy": "2026-06-19T07:41:30.642023Z",
     "iopub.status.idle": "2026-06-19T07:41:30.748891Z",
     "shell.execute_reply": "2026-06-19T07:41:30.747611Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "mixed FC + spectral loss: 1.7818667888641357\n"
     ]
    }
   ],
   "source": [
    "mixed = objectives.combine(\n",
    "    (1.0, objectives.fc_corr(as_loss=True)),\n",
    "    (0.1, peak_frequency_mse(dt_ms=1.0)),\n",
    ")\n",
    "print(\"mixed FC + spectral loss:\", float(mixed(pred_sig, fast_sig)))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "823dd9bf8f46",
   "metadata": {},
   "source": [
    "## Plug an objective into `Fitter`\n",
    "\n",
    "`brainmass.Fitter` accepts an objective two ways:\n",
    "\n",
    "- **`objective=`** — a `(prediction, target)` callable; the Fitter builds the loss\n",
    "  as `objective(predict(model), target) + model.reg_loss()`, and you supply\n",
    "  `predict` (a `Simulator` closure) and the `target`.\n",
    "- **`loss_fn=`** — a full `loss_fn(model) -> (scalar, aux)` that you own entirely,\n",
    "  for losses against a *derived scalar summary* (the right choice for oscillators,\n",
    "  where a point-by-point `timeseries_rmse` is phase-degenerate).\n",
    "\n",
    "Here we fit the Hopf bifurcation parameter `a` so the settled amplitude matches a\n",
    "target, using a custom scalar objective inside `loss_fn`.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "77d1ffb40d2f",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2026-06-19T07:41:30.751126Z",
     "iopub.status.busy": "2026-06-19T07:41:30.750887Z",
     "iopub.status.idle": "2026-06-19T07:41:31.173385Z",
     "shell.execute_reply": "2026-06-19T07:41:31.172291Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "target amplitude : 1.2\n",
      "fitted a         : 1.232\n",
      "best loss        : 8.41e-03\n"
     ]
    }
   ],
   "source": [
    "def amplitude_loss(target_amp):\n",
    "    def loss(prediction):\n",
    "        x = u.get_magnitude(prediction)\n",
    "        amp = jnp.sqrt(jnp.mean(x ** 2)) * jnp.sqrt(2.0)\n",
    "        return (amp - target_amp) ** 2, amp\n",
    "    return loss\n",
    "\n",
    "target_amp = 1.2\n",
    "obj = amplitude_loss(target_amp)\n",
    "\n",
    "# Mark `a` trainable with a Param; a plain Param(fit=True) lives in physical space.\n",
    "node = brainmass.HopfStep(in_size=1, a=Param(0.1, fit=True), w=0.3,\n",
    "                          init_x=braintools.init.Constant(0.5))\n",
    "\n",
    "def loss_fn(m):\n",
    "    res = brainmass.Simulator(m, dt=0.1 * u.ms).run(\n",
    "        150 * u.ms, monitors=['x'], transient=50 * u.ms)\n",
    "    return obj(res['x'])   # (scalar loss, amplitude aux)\n",
    "\n",
    "fitter = brainmass.Fitter(node, braintools.optim.Adam(lr=0.05), loss_fn=loss_fn)\n",
    "result = fitter.fit(n_steps=40)\n",
    "print(f\"target amplitude : {target_amp}\")\n",
    "print(f\"fitted a         : {result.best_params['a']:.3f}\")\n",
    "print(f\"best loss        : {result.best_loss:.2e}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "913163788e73",
   "metadata": {},
   "source": [
    "The fit drives `a` from `0.1` toward the value whose limit cycle has amplitude\n",
    "`1.2`. Plot the loss curve from `result.history`:\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "fade74e41902",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2026-06-19T07:41:31.176242Z",
     "iopub.status.busy": "2026-06-19T07:41:31.176001Z",
     "iopub.status.idle": "2026-06-19T07:41:31.385853Z",
     "shell.execute_reply": "2026-06-19T07:41:31.384934Z"
    }
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAekAAAEiCAYAAADd4SrgAAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjgsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvwVt1zgAAAAlwSFlzAAAPYQAAD2EBqD+naQAATHpJREFUeJzt3XdUFNffBvBnl96LIE2QKooFIyIiVUUxltgiGhNFY4klMcZU09Qkb0w0xRixpai/JPaoKfaogCL2XpEqioCiIh3cve8fhI1LE3RxF3w+53AOOzvMfPfOss9OuXMlQggBIiIi0jhSdRdARERE1WNIExERaSiGNBERkYZiSBMREWkohjQREZGGYkgTERFpKIY0ERGRhmJIExERaSiGNBERkYZiSBM9ojFjxsDZ2VlpmkQiwezZs5/I+mfPng2JRPJE1vWglStXQiKR4NixYw+dNzQ0FKGhoQ1fVCXR0dGQSCSIjo5+4usmUiWGtAZLSkrCK6+8AldXV+jr68PU1BQBAQH47rvvUFRU1CDrXL16NRYsWNAgy34aHTx4ELNnz8bdu3fVXUqTtHjxYqxcuVLdZRA1GG11F0DV27p1K4YNGwY9PT2MHj0a7dq1Q2lpKQ4cOIC3334b58+fx/Lly1W+3tWrV+PcuXOYPn26ypf9NCgqKoK29n//VgcPHsScOXMwZswYmJubq68wNdm1a1eDLn/x4sWwsrLCmDFjlKYHBwejqKgIurq6Dbp+oobGkNZAKSkpGDFiBFq2bIm9e/fCzs5O8dzUqVORmJiIrVu3qrFCqom+vr66S9Ao6gpJqVTKbfGEFBYWwtDQUN1lNFk83K2B5s2bh/z8fPz0009KAV3B3d0dr7/+OgAgNTUVEomk2kN+lc+P5uXlYfr06XB2doaenh6aN2+OXr164cSJEwDKzx9u3boVaWlpkEgkkEgkSudcs7OzMW7cONjY2EBfXx/e3t5YtWqV0jor6vnqq68QFRUFV1dXGBoaonfv3khPT4cQAp9++ilatGgBAwMDDBw4ELdv335om5w5cwZjxoxRHPq3tbXFyy+/jJycHKX5Ks7TJiQk4KWXXoKZmRmsra3x0UcfQQiB9PR0DBw4EKamprC1tcXXX3+t9PcV5zLXrVuH999/H7a2tjAyMsJzzz2H9PT0h9b5YJvPnj0bb7/9NgDAxcVF0aapqan12m4AcODAAfj6+kJfXx9ubm5YtmxZjTX8+uuv8PHxgYGBASwtLTFixIg61Q4AJ0+exLPPPgtTU1MYGxujZ8+eOHToULXzFhYW4pVXXkGzZs1gamqK0aNH486dO0rzVHdOuqSkBLNmzYK7uzv09PTg6OiId955ByUlJdW+li5dusDQ0BAWFhYIDg5W7J07Ozvj/PnziImJUbRtxboqn5N+9dVXYWxsjMLCwirreOGFF2BrawuZTKaYtn37dgQFBcHIyAgmJibo168fzp8/r/R3ZWVluHTpEm7cuFFrm1a4dOkSIiIiYG1tDQMDA3h6euKDDz5Qmqcu7V9xTUBcXBxmzJgBa2trGBkZYfDgwbh586Zivv79+8PV1bXaWvz9/dG5c2elaXV534SGhqJdu3Y4fvw4goODYWhoiPfffx8AkJOTg1GjRsHU1BTm5uaIjIzE6dOnq32fX7p0Cc8//zwsLS2hr6+Pzp07488//3yk11lh+/btCAkJgYmJCUxNTeHr64vVq1crzXP48GH06dMHZmZmMDQ0REhICOLi4qptI40hSOM4ODgIV1fXOs2bkpIiAIgVK1ZUeQ6AmDVrluLxyJEjha6urpgxY4b48ccfxZdffikGDBggfv31VyGEELt27RIdO3YUVlZW4pdffhG//PKL2Lx5sxBCiMLCQtGmTRuho6Mj3njjDbFw4UIRFBQkAIgFCxZUqadjx47Cy8tLfPPNN+LDDz8Uurq6omvXruL9998X3bp1EwsXLhTTpk0TEolEjB079qGv86uvvhJBQUHik08+EcuXLxevv/66MDAwEF26dBFyuVwx36xZsxTrf+GFF8TixYtFv379BADxzTffCE9PTzF58mSxePFiERAQIACImJgYxd/v27dPABDt27cXHTp0EN9884147733hL6+vmjVqpUoLCxUzBsZGSlatmxZY5ufPn1avPDCCwKA+PbbbxVtmp+fX6/tdubMGWFgYCCcnJzE3LlzxaeffipsbGxEhw4dROV/4c8++0xIJBIxfPhwsXjxYjFnzhxhZWUlnJ2dxZ07d2pt43PnzgkjIyNhZ2cnPv30U/HFF18IFxcXoaenJw4dOqSYb8WKFYo2CgoKEgsXLhRTp04VUqlUBAcHK22PkJAQERISongsk8lE7969haGhoZg+fbpYtmyZePXVV4W2trYYOHCgUj2zZ88WAES3bt3E/PnzxXfffSdGjhwp3n33XSGEEJs3bxYtWrQQrVu3VrTtrl27lLbjvn37hBBCxMbGCgBi/fr1SusoKCgQRkZGYurUqYpp//vf/4REIhF9+vQR33//vfjyyy+Fs7OzMDc3FykpKYr5KrZhZGRkre0qRPl7wdTUVDRr1kzMnDlTLFu2TLzzzjuiffv2j9z+zzzzjOjRo4f4/vvvxZtvvim0tLRERESE0usAII4cOaJUS2pqqgAg5s+fr5hW1/dNSEiIsLW1FdbW1uK1114Ty5YtE1u2bBEymUz4+/sLLS0t8eqrr4pFixaJXr16CW9v7yrv83PnzgkzMzPh5eUlvvzyS7Fo0SIRHBwsJBKJ2LRpU71fZ8W8EolEtGvXTvzf//2fiIqKEuPHjxejRo1SzLNnzx6hq6sr/P39xddffy2+/fZb0aFDB6GrqysOHz780G2oLgxpDZObmysAVPnAqkl9PuzNzMyUPoyq069fvyrBI4QQCxYsEAAUgS6EEKWlpcLf318YGxuLe/fuKdVjbW0t7t69q5h35syZAoDw9vYWZWVliukvvPCC0NXVFcXFxbXW9WA4VlizZo0AIGJjYxXTKkJ64sSJimn3798XLVq0EBKJRHzxxReK6Xfu3BEGBgZKH7IVH+4ODg6K1ySEEOvXrxcAxHfffaeY9rCQFkKI+fPnCwBKH+5C1G+7DRo0SOjr64u0tDTFtAsXLggtLS2lkE5NTRVaWlri//7v/5SWd/bsWaGtrV1lemWDBg0Surq6IikpSTEtIyNDmJiYiODgYMW0ig9PHx8fUVpaqpg+b948AUD88ccfimmVQ/qXX34RUqlU7N+/X2ndS5cuFQBEXFycEEKIK1euCKlUKgYPHixkMpnSvA9+CWjbtq3S8itUDmm5XC4cHBzE0KFDlear2K4V76G8vDxhbm4uJkyYoDRfZmamMDMzU5pen5AODg4WJiYmStuw8mupb/uHhYUp/f0bb7whtLS0FP93ubm5Qk9PT7z55ptK65w3b56QSCSKWurzvgkJCREAxNKlS5Xm/f3336t8YZfJZKJHjx5V3uc9e/YU7du3V/qfl8vlolu3bsLDw6Per/Pu3bvCxMRE+Pn5iaKiomrbVy6XCw8PDxEeHq60rMLCQuHi4iJ69eolNBUPd2uYe/fuAQBMTExUvmxzc3McPnwYGRkZ9f7bbdu2wdbWFi+88IJimo6ODqZNm4b8/HzExMQozT9s2DCYmZkpHvv5+QEAXnrpJaULq/z8/FBaWorr16/Xun4DAwPF78XFxbh16xa6du0KAIrD9Q8aP3684nctLS107twZQgiMGzdOMd3c3Byenp5ITk6u8vejR49W2gbPP/887OzssG3btlrrVDWZTIadO3di0KBBcHJyUkxv06YNwsPDlebdtGkT5HI5IiIicOvWLcWPra0tPDw8sG/fvlrXs2vXLgwaNEjpEKmdnR1GjhyJAwcOKN6bFSZOnAgdHR3F48mTJ0NbW7vWNtqwYQPatGmD1q1bK9XYo0cPAFDUuGXLFsjlcnz88ceQSpU/ph6l25lEIsGwYcOwbds25OfnK6avW7cODg4OCAwMBADs3r0bd+/exQsvvKBUn5aWFvz8/JTa0NnZGUKIh15dfvPmTcTGxuLll19W2oYPvpZHbf8H2yIoKAgymQxpaWkAAFNTUzz77LNYv349hBBKr7lr166KWur7vtHT08PYsWOVpu3YsQM6OjqYMGGCYppUKsXUqVOV5rt9+zb27t2LiIgI5OXlKdaVk5OD8PBwXLlypcpnwcNe5+7du5GXl4f33nuvynUIFX936tQpXLlyBSNHjkROTo5ivQUFBejZsydiY2Mhl8uhiXjhmIYxNTUFUH7+WNXmzZuHyMhIODo6wsfHB3379sXo0aNrPG/1oLS0NHh4eFT5wGzTpo3i+QdV/jCqCGxHR8dqp1c+l1nZ7du3MWfOHKxduxbZ2dlKz+Xm5laZv7r16+vrw8rKqsr0yue1AcDDw0PpsUQigbu7O1JTU2utU9Vu3ryJoqKiKvUAgKenp1IgXrlyBUKIaucFoBSo1a2nsLAQnp6eVZ5r06YN5HI50tPT0bZtW8X0yusxNjaGnZ1drW105coVXLx4EdbW1tU+X7Ftk5KSIJVK4eXlVeOy6mv48OFYsGAB/vzzT4wcORL5+fnYtm0bXnnlFcWH+ZUrVwBA8aWhsor/z/qo+BLYrl27Gud5lPav/B63sLAAoPy/NHz4cGzZsgXx8fHo1q0bkpKScPz4caVulvV93zg4OFS5IDAtLQ12dnZVLiBzd3dXepyYmAghBD766CN89NFH1a4vOzsbDg4OdX6dSUlJAGpv34rtGhkZWeM8ubm5imVrEoa0hjE1NYW9vT3OnTtXp/lr2qt48CKYChEREQgKCsLmzZuxa9cuzJ8/H19++SU2bdqEZ5999rHqrkxLS6te0x/8pl+diIgIHDx4EG+//TY6duwIY2NjyOVy9OnTp9pvwNWt51HX3RDqs93qSi6XQyKRYPv27dW+VmNj40detqrI5XK0b98e33zzTbXPV/4Sp0pdu3aFs7Mz1q9fj5EjR+Kvv/5CUVERhg8frlQfAPzyyy+wtbWtsowHjwKpW13ezwMGDIChoSHWr1+Pbt26Yf369ZBKpRg2bJhinvq+bx48qlVfFe371ltvVTkSVKFysKvi/7ZivfPnz0fHjh2rnUcT/j+qoznvOFLo378/li9fjvj4ePj7+9c6b8U3v8o3y6i8Z1vBzs4OU6ZMwZQpU5CdnY1OnTrh//7v/xQhXVN4tGzZEmfOnIFcLlfam7506ZLi+YZy584d7NmzB3PmzMHHH3+smF7x7bghVF62EAKJiYno0KFDvZZTU3vWdbtVXAlc3Wu9fPmy0mM3NzcIIeDi4oJWrVrVq05ra2sYGhpWWSZQvo2lUmmVAL1y5Qq6d++ueJyfn48bN26gb9++Na7Hzc0Np0+fRs+ePWs9bO3m5ga5XI4LFy7U+KEK1P/Qd0REBL777jvcu3cP69atg7Ozs+K0ScV6AaB58+YICwur17JrUnGkqrYv3o/S/nVhZGSE/v37Y8OGDfjmm2+wbt06BAUFwd7eXjHP47xvKrRs2RL79u2r0h0rMTFRab6KttDR0VFZ+1Zss3PnzlUJ+MrzmJqaqmy9TwrPSWugd955B0ZGRhg/fjyysrKqPJ+UlITvvvsOQPmbzsrKCrGxsUrzLF68WOmxTCarcli4efPmsLe3V+r6YmRkVO3h4759+yIzMxPr1q1TTLt//z6+//57GBsbIyQkpP4vtI4qvklX/ubckHdG+9///qd0ymHjxo24ceNGvY84GBkZAagaxnXdblpaWggPD8eWLVtw9epVxfSLFy9i586dSvMOGTIEWlpamDNnTpW2EkJUe1j/wfX07t0bf/zxh9Lh6qysLKxevRqBgYFVDvUuX74cZWVlisdLlizB/fv3a22jiIgIXL9+HT/88EOV54qKilBQUAAAGDRoEKRSKT755JMqR0oefG1GRkb1upvb8OHDUVJSglWrVmHHjh2IiIhQej48PBympqb4/PPPlV5bhQe7/tS1C5a1tTWCg4Px888/K23DB1/Lo7R/XQ0fPhwZGRn48ccfcfr0aaUjB8DjvW8qhIeHo6ysTGm7yuVyREVFKc3XvHlzhIaGYtmyZdW2W3Vdqx6md+/eMDExwdy5c1FcXFylfgDw8fGBm5sbvvrqK6VrEh5nvU8K96Q1kJubG1avXo3hw4ejTZs2SnccO3jwIDZs2KB0h6Xx48fjiy++wPjx49G5c2fExsYiISFBaZl5eXlo0aIFnn/+eXh7e8PY2Bj//PMPjh49qtRX2MfHB+vWrcOMGTPg6+sLY2NjDBgwABMnTsSyZcswZswYHD9+HM7Ozti4cSPi4uKwYMGCBrnQrYKpqSmCg4Mxb948lJWVwcHBAbt27UJKSkqDrdPS0hKBgYEYO3YssrKysGDBAri7uytdGFMXPj4+AIAPPvgAI0aMgI6ODgYMGKD4Evaw7QYAc+bMwY4dOxAUFIQpU6Yovhy1bdsWZ86cUczn5uaGzz77DDNnzkRqaioGDRoEExMTpKSkYPPmzZg4cSLeeuutGmv97LPPsHv3bgQGBmLKlCnQ1tbGsmXLUFJSgnnz5lWZv7S0FD179kRERAQuX76MxYsXIzAwEM8991yN6xg1ahTWr1+PSZMmYd++fQgICIBMJsOlS5ewfv167Ny5E507d4a7uzs++OADfPrppwgKCsKQIUOgp6eHo0ePwt7eHnPnzlW075IlS/DZZ5/B3d0dzZs3r/F8MgB06tRJseySkpIqgWVqaoolS5Zg1KhR6NSpE0aMGAFra2tcvXoVW7duRUBAABYtWgQAuH79Otq0aYPIyMiHXjy2cOFCBAYGolOnTpg4cSJcXFyQmpqKrVu34tSpU4/U/nXVt29fmJiY4K233oKWlhaGDh2q9Pzjvm+A8i9VXbp0wZtvvonExES0bt0af/75p+IeCA8e8YiKikJgYCDat2+PCRMmwNXVFVlZWYiPj8e1a9dw+vTper0+U1NTfPvttxg/fjx8fX0xcuRIWFhY4PTp0ygsLMSqVasglUrx448/4tlnn0Xbtm0xduxYODg44Pr169i3bx9MTU3x119/1Wu9T8wTvZac6iUhIUFMmDBBODs7C11dXWFiYiICAgLE999/r9R9obCwUIwbN06YmZkJExMTERERIbKzs5W68pSUlIi3335beHt7CxMTE2FkZCS8vb3F4sWLldaZn58vRo4cKczNzQUApS5GWVlZYuzYscLKykro6uqK9u3bV+lCVNEt5cE+mEL81yVmw4YNStMrulkcPXq01ra4du2aGDx4sDA3NxdmZmZi2LBhIiMjo0p3pYouWDdv3lT6+8jISGFkZFRluSEhIaJt27ZV6lyzZo2YOXOmaN68uTAwMBD9+vWr0n2mLl2whBDi008/FQ4ODkIqlSp1x6rLdqsQExMjfHx8hK6urnB1dRVLly5VvNbKfv/9dxEYGCiMjIyEkZGRaN26tZg6daq4fPlyDa37nxMnTojw8HBhbGwsDA0NRffu3cXBgweV5qnYZjExMWLixInCwsJCGBsbixdffFHk5ORUad/KXaRKS0vFl19+Kdq2bSv09PSEhYWF8PHxEXPmzBG5ublK8/7888/imWeeUcwXEhIidu/erXg+MzNT9OvXT5iYmAgAinVV7oL1oA8++EAAEO7u7jW2w759+0R4eLgwMzMT+vr6ws3NTYwZM0YcO3ZMMU99umAJUd4/uOI9rK+vLzw9PcVHH32kNE992r/y/0xtr/nFF19UdGeqSV3eN5X/Xx508+ZNMXLkSGFiYiLMzMzEmDFjRFxcnAAg1q5dqzRvUlKSGD16tLC1tRU6OjrCwcFB9O/fX2zcuPGRX+eff/4punXrJgwMDISpqano0qWLWLNmjdI8J0+eFEOGDBHNmjUTenp6omXLliIiIkLs2bOnxnZRN4kQarhqhkhDRUdHo3v37tiwYQOef/55dZfT6AUFBUFPTw///POPukshNdiyZQsGDx6MAwcOICAgQN3lNEo8J01EDebGjRtVur1R01R5ZD6ZTIbvv/8epqam6NSpk5qqavx4TpqIVO7gwYPYtGkTkpKS8O6776q7HHoCXnvtNRQVFcHf3x8lJSXYtGkTDh48iM8///yxum097RjSRKRyP/zwA7Zv347p06dXuTsVNU09evTA119/jb///hvFxcVwd3fH999/j1dffVXdpTVqPCdNRESkoXhOmoiISEMxpImIiDQUz0nXQi6XIyMjAyYmJo808g4REVF1hBDIy8uDvb19lYGLHsSQrkVGRkaD3vCfiIiebunp6WjRokWNzzOka1Fxq8v09PRHvm8uERFRZffu3YOjo+NDb6nMkK5FxSFuU1NThjQREancw06l8sKxakRFRcHLywu+vr7qLoWIiJ5i7Cddi3v37sHMzAy5ubnckyYiIpWpa75wT5qIiEhDMaSJiIg0FEO6GjwnTUREmoDnpGuhynPSPx1IgYuVIXq0tlFRdURE1FjVNV/YBesJOHH1Dv5v6wXIBTAhyAVvh7eGrjYPYhARUe2YFE9AW3tTjPZ3BgD8sD8Fw5YexNWcQvUWRUREGo8h/QToaWth9nNtsWyUD8wMdHD6Wi76LdyPrWduqLs0IiLSYAzpJyi8rS22vR4En5YWyCu5j6mrT+D9zWdRXCZTd2lERKSBGNLVaMirux3MDbB2YldMCXWDRAKsPnwVg6LikJidp/J1ERFR48aru2vR0Hcc23/lJt5Ydwq38kthoKOFOQPbYphPCw6LSUTUxPGOY41AkIc1tr0ehAD3Zigqk+GdjWfwxrpTyCsuU3dpRESkARjSatbcRB//e9kPb4d7QksqwZZTGeizYD/2X7mp7tKIiEjNGNIaQEsqwdTu7lg7sSscLQ1w/W4RRv10BDM3nUV+yX11l0dERGrCkNYgvs6W2PF6MEb7twQArDlyFeHfxuLAlVtqroyIiNSBIV0Ndd6720hPG58MbIfVE/zQwqJ8r/qlnw7jg83cqyYietrw6u5aqHs86YKS+/hi+yX8cigNQHn3rfnPd0A3d6snXgsREakOr+5uAoz0tPHpoHZYPf6/veqRPx7Gh1vOooB71URETR5DuhHo5m6FHdOD8VJXJwDAr4euove3sdh3KVvNlRERUUNiSDcSxnra+GxQe/w23g8O5uV71WNXHsXU304g616xussjIqIGwJBuZALcrbB7RjAmBrtCSyrB1rM3EPZ1DP4XnwqZnJcXEBE1JQzpRshQVxvv922Dv14NREdHc+SV3MfHf5zHkMVxOJ+Rq+7yiIhIRRjSjZiXvSl+n9wNnw5sCxM9bZy+lovnFsXh/7Ze4IVlRERNAEO6GursJ11fWlIJRvk7Y8+bIejXwQ4yucAP+1PQ+9tY7LmYpe7yiIjoMbCfdC3U3U/6Uey7nI2PtpzDtTtFAICwNjb4uL8XnJoZqrkyIiKqwH7ST6nuns2x+40QTApxg7ZUgn8uZiHs2xh8szsBRaUydZdHRET1wD3pWjTGPekHJWbnYfafF3Agsfze3w7mBviovxfC29pwzGoiIjWqa74wpGvR2EMaAIQQ2HEuE5/+fQEZueX9qYM8rDD7ubZwszZWc3VERE8nhrQKNIWQrlBUKsPi6EQsi0lGqUwOHS0JXg5wwWs9PWCsp63u8oiInioMaRVoSiFdIfVWAT79+wL2/HtL0eYmeni/bxsM7GjPQ+BERE8IQ1oFmmJIV9h7KQtz/rqAtJxCAMAzTub4uL8XnnGyUHNlRERNH0NaBZpySANAcZkMPx1IQdS+RBT+e+X3oI72eKdPa9ibG6i5OiKipotdsB5DY7qZyePQ19HC1O7uiH4rFM/7tAAAbDmVgR5fR2PBP+yyRUSkbtyTrkVT35Ou7Oy1XHzy93kcTb0DALAz08e7fVrjOW97SKU8X01EpCo83K0CT1tIA+VdtradzcTn2y7i+t3yu5Z1dDTHxwO80Innq4mIVIIhrQJPY0hXqO589QBve7wT7glHS95ilIjocTCkVeBpDukK2feKMX/nZWw8cQ1CALpaUkR2a4lXu3vAzFBH3eURETVKDGkVYEj/53xGLj7fdhFxiTkAADMDHUzr6YFRXVtCV5vXHxIR1QdDWgUY0sqEEIhOuIm52y4iISsfAOBkaYh3+7RG3/a2vBkKEVEdMaRVgCFdvfsyOTYev4avdyfgZl4JAKCTkzk+6NcGPi0t1VwdEZHmY0irAEO6dgUl97E8NhnLY5NRVFZ+cVmftrZ4u48nB+8gIqoFQ1oFGNJ1k3WvGN/uTsD6Y+mQC0BLKsEIX0e83tMDzU311V0eEZHGYUirAEO6fi5n5mH+zkv452L54B0GOloYH+SCicGuMNHnleBERBUY0irAkH40h5Nz8MWOSzh59S4AwNJIF692d8eLXZ2gp62l3uKIiDQAQ/oxREVFISoqCjKZDAkJCQzpRyCEwM7zWZi38xKSbxYAABwtDfBmL0/eZpSInnoMaRXgnvTjuy+TY/2xa1jwTwKy/70S3MvOFO/08URIK2t22yKipxJDWgUY0qpTWHofK+JSsTQ6CXkl9wEAXVws8W4fT3bbIqKnDkNaBRjSqne7oBSL9yXif4fSUHpfDgAIa9Mcb4V7orUt25iIng4MaRVgSDecjLtF+O6fK9hwvLzblkQCDOrogDfCWsGpGQfwIKKmjSGtAgzphpeYnY9vdl/GtrOZAAAdLQlG+DrhtR7u7GNNRE0WQ1oFGNJPztlruZi38xL2X7kFANDXkWJsgAteCXaFuaGumqsjIlIthrQKMKSfvINJtzBvx2WcSr8LADDR18bEIFeMDXSBsZ62eosjIlIRhrQKMKTVQwiB3Rey8PWuBFzOygNQfkOUKaFueKlrS+jr8IYoRNS4MaRVgCGtXnK5wF9nMrDgnytIuVV+QxQbUz281sMDEZ0dOY41ETVaDGkVYEhrhvsyOX4/cQ3f/XMFGbnFAMrvXja9ZysMesYBWrx7GRE1MgxpFWBIa5aS+zKsOXwVi/Yl4VZ++d3L3Jsb442wVni2nS1vNUpEjQZDWgUY0pqpqFSGVfGpWBqThLuFZQCA1rYmmNGrFXp52fBWo0Sk8RjSKsCQ1mz3isvw84EU/LQ/RXGr0Q4tzPBGr1YI5X3BiUiDMaRVgCHdONwtLMXy2GSsPJiKwlIZAKCTkzlm9PJEgHszhjURaRyGtAowpBuXnPwSLI1Jwv/i01Dy733B/Vws8WZvT3Rx4SAeRKQ5GNIqwJBunLLvFWNxdBJWH76KUll5WAe6W+GNXh4ccYuINAJD+jFERUUhKioKMpkMCQkJDOlGKuNuERbtS8T6o+m4Ly9/mwe3ssYbYR54xslCzdUR0dOMIa0C3JNuGtJvF2LR3kRsPHENsn/DOtTTGm+EtYK3o7l6iyOipxJDWgUY0k3L1ZxCfL/3CjadvK4I656tm2N6WCu0b2Gm5uqI6GnCkFYBhnTTlHqrAAv3XsGWk9fxb1YjrI0Npod5oJ0Dw5qIGh5DWgUY0k1b8s18fL83EX+c+i+sw9vaYFpPD7S1Z1gTUcNhSKsAQ/rpkJidj4V7ruCvMxmo+G/o7WWD18MY1kTUMBjSKsCQfrokZudh4Z7EKmE9rScPgxORajGkVYAh/XRKzM7D93sT8efp/8K6l5cNXmdYE5GKMKRVgCH9dGNYE1FDYUirAEOagPJz1ov2XsGfpzMeuBq8Oab19ECHFuZqrY2IGieGtAowpOlBSTfzsajS1eDdPa0xrSfvYEZE9cOQVgGGNFUn+WY+ovYlYcup/26KEuRhhelhvDc4EdUNQ1oFGNJUm9RbBVgcnYjfT/wX1gHuzfB6z1YcdYuIasWQVgGGNNVF+u1CLI5OxIZj1xQDeXR1tcS0nh7wd+V41kRUFUNaBRjSVB/X7hRiSXQS1h9LR5ms/N+qc0sLTOvpgSAPK4Y1ESkwpFWAIU2PIuNuEZbGJGHt0XSU3i8fz9rb0RzTerijR+vmDGsiYkirAkOaHkfWvWIsi0nG6iNpKC4rD+u29qZ4rYcHenvZQCplWBM9rRjSKsCQJlW4mVeCH/cn45dDaSgslQEAPG1M8FpPdzzbzg5aDGuipw5DWgUY0qRKtwtK8dOBZKw6mIb8kvsAADdrI0zt7o7nvO2hrSVVc4VE9KQwpFWAIU0NIbewDD/HpWBFXAruFZeHtZOlIaaEumFIpxbQ1WZYEzV1DGkVYEhTQ7pXXIZf4tPw04EU3C4oBQDYm+ljUqgbIjo7Ql9HS80VElFDYUirAEOanoTC0vtYffgqlsUm42ZeCQDA2kQPrwS7YqSfEwx1tdVcIRGpWl3z5ak4rjZ48GBYWFjg+eefV3cpRFUY6mpjfJAr9r/THZ8MbAt7M33czCvBZ1svIvDLfYjal4i84jJ1l0lEavBU7ElHR0cjLy8Pq1atwsaNG+v8d9yTJnUovS/H5pPXsDg6CWk5hQAAU31tjOnmjLEBLrAw0lVzhUT0uLgn/YDQ0FCYmJiouwyiOtHVlmK4rxP2zAjBt8O94d7cGPeK72Ph3kQEfLkXn2+7iOy8YnWXSURPgNpDOjY2FgMGDIC9vT0kEgm2bNlSZZ6oqCg4OztDX18ffn5+OHLkyJMvlOgJ09aSYvAzLbBrejCWvNgJXnamKCyVYXlsMgK/3IeP/ziH63eL1F0mETUgtYd0QUEBvL29ERUVVe3z69atw4wZMzBr1iycOHEC3t7eCA8PR3Z2tmKejh07ol27dlV+MjIyntTLIGowUqkEz7a3w9ZpgVgxxhednMxRel+O/8WnIWTePryz8TRSbhWou0wiagAadU5aIpFg8+bNGDRokGKan58ffH19sWjRIgCAXC6Ho6MjXnvtNbz33nt1XnZ0dDQWLVpU6znpkpISlJSUKB7fu3cPjo6OPCdNGkUIgfjkHETtS0RcYg4AQCoB+nWwx9Tubmhty/cqkaZr0HPSq1atwtatWxWP33nnHZibm6Nbt25IS0t7lEVWq7S0FMePH0dYWJhimlQqRVhYGOLj41W2ngpz586FmZmZ4sfR0VHl6yB6XBKJBN3crPDb+K7YNKUberZuDrkA/jqdgT4L9mP8qmM4efWOusskIhV4pJD+/PPPYWBgAACIj49HVFQU5s2bBysrK7zxxhsqK+7WrVuQyWSwsbFRmm5jY4PMzMw6LycsLAzDhg3Dtm3b0KJFixoDfubMmcjNzVX8pKenP1b9RA2tk5MFfhrji63TAtGvvR0kEuCfi1kYvPggXvzxEA4m3YIGHSwjonp6pLskpKenw93dHQCwZcsWDB06FBMnTkRAQABCQ0NVWZ9K/PPPP3WaT09PD3p6eg1cDZHqtbU3Q9SLnZB0Mx9LopOw5eR1xCXmIC4xB52czDG1O4fJJGqMHmlP2tjYGDk55efCdu3ahV69egEA9PX1UVSkuqtNraysoKWlhaysLKXpWVlZsLW1Vdl6iJoKN2tjfDXMG/veCsWori2hqy3Fiat3MW7VMfRdeAB/n8mATM49a6LG4pFCulevXhg/fjzGjx+PhIQE9O3bFwBw/vx5ODs7q6w4XV1d+Pj4YM+ePYppcrkce/bsgb+/v8rWQ9TUOFoa4tNB7XDgne54JdgVRrpauHjjHl5dfRJh38Rg/dF0lN6Xq7tMInqIRwrpqKgo+Pv74+bNm/j999/RrFkzAMDx48fxwgsv1GtZ+fn5OHXqFE6dOgUASElJwalTp3D16lUAwIwZM/DDDz9g1apVuHjxIiZPnoyCggKMHTv2UUqvk6ioKHh5ecHX17fB1kH0JDQ31cfMvm0Q914PvN7TA2YGOki5VYB3fj+DkPn7sCIuBUX/jnFNRJpH7V2woqOj0b179yrTIyMjsXLlSgDAokWLMH/+fGRmZqJjx45YuHAh/Pz8Grw23haUmpr8kvtYfTgNP+xPUQzm0cxIFy8HumCUf0uY6uuouUKip0ODjoK1Y8cOGBsbIzAwEED5nucPP/wALy8vREVFwcLC4tEr1yAMaWqqistk2Hj8GpbGJOHanfLrSEz0tDG6W0uMDXCBlTEvoCRqSA3aT/rtt9/GvXv3AABnz57Fm2++ib59+yIlJQUzZsx4tIqJ6InR19HCS11bYt9bofh2uDc8mhsjr+Q+ovYlIfDLvZj953necpRIAzzSnrSxsTHOnTsHZ2dnzJ49G+fOncPGjRtx4sQJ9O3bt159mDVRVFQUoqKiIJPJkJCQwD1pavLkcoFdF7KwODoRZ67lAgC0pRIMesYBk0Lc4N7cWM0VEjUtDXq429LSEgcOHICXlxcCAwMxevRoTJw4EampqfDy8kJhYeFjFa8peLibnjZCCBxIvIXF+5IQn1zezVIiAfq0tcWUUHe0b2Gm5gqJmoa65ssj3cwkMDAQM2bMQEBAAI4cOYJ169YBABISEtCiRYtHq5iI1E4ikSDIwxpBHtY4cfUOFu9Lwj8Xs7D9XCa2n8tEkIcVpnZ3h5+LJW+MQvQEPNI56UWLFkFbWxsbN27EkiVL4ODgAADYvn07+vTpo9ICiUg9OjlZ4MfIztg5PRiDn3GAllSC/VduYcTyQxi65CD2XMziLUeJGpjau2BpMh7uJvpP+u1CLItNwvpj1xQ3Qmlta4LJoW7o194O2lpqH/mWqNFo0HPSACCTybBlyxZcvHgRANC2bVs899xz0NLSerSKNRBDmqiq7Lxi/HQgBb8duor8kvsAAEdLA0wMdsMwnxbQ12k6nwFEDaVBQzoxMRF9+/bF9evX4enpCQC4fPkyHB0dsXXrVri5uT165RqAV3cTPVxuYRl+OZSKn+NScbugFABgZayHcYEueKmrE0x4YxSiGjVoSPft2xdCCPz222+wtLQEAOTk5OCll16CVCpVGmu6MeOeNNHDFZXKsO7oVfywP0XRt9pEXxuj/XljFKKaNGhIGxkZ4dChQ2jfvr3S9NOnTyMgIAD5+fn1r1gDMaSJ6q5MJscfpzKwNCYJidnlnwF62lKM8HXE+CBXOFoaqrlCIs3RoHcc09PTQ15eXpXp+fn50NXVfZRFElEjp6MlxfM+LbBrejCWvuQD7xZmKLkvx6r4NIR+FY0Z604hIavq5wYR1eyRQrp///6YOHEiDh8+DCEEhBA4dOgQJk2ahOeee07VNRJRIyKVStCnnS22TA3Ab+P9EOhuBZlcYNPJ6+j9bSzGrzqG42l31F0mUaPwSIe77969i8jISPz111/Q0Sm/OKSsrAwDBw7EihUrYG5uruo61YKHu4lU48y1u1gSnYQd5zNR8Ynj52KJyaFuCGllzRuj0FOnwbtgAeVXeVd0wWrTpg3c3d0fdVEahVd3EzWMpJv5WB6TjE0nr6FMVv7R42Vnismhbujb3g5aUoY1PR1UHtL1Gd3qm2++qfO8mox70kQN40ZuEX7an4LVR66isFQGAGjZzBCvBLthSCcH9rWmJk/lId29e/c6rVgikWDv3r11q1LDMaSJGtadglKsik/FqoOpuFNYBgCwNinva/2iH/taU9P1RA53N3UMaaIno7D0PtYeSceP+5ORkVsMoLyv9aiu5X2trU3Y15qaFoa0CjCkiZ6s0vty/Hlaua+1rrYUEZ1bYGKQG5yasa81NQ0MaRVgSBOph1wu8M/FLCyOTsKp9LsAAKkE6N/BHpNC3OBlz/9HatwY0irAkCZSLyEEDqfcxpLoJMQk3FRMD/W0xqQQN45rTY0WQ1oFGNJEmuPc9Vwsi03G1jMZkP/7qfWMkzkmh7ghrI0NpOy+RY0IQ/oxsJ80keZKyynA8thkbDj+37jW7s2N8UqwKwZ2dICuNse1Js3HkFYB7kkTaa6beSVYEZeCX+LTkPfvuNZ2ZvoYH+SKEb6OMNLTVnOFRDVjSKsAQ5pI890rLsPqw1fx04EU3MwrAQCYGegg0r8lIrs5oxmHyiQNxJBWAYY0UeNRXCbDphPXsTw2Cak5hQAAfR0pRvg6YXyQC1pYsPsWaQ6GtAowpIkaH5lcYOf5TCyJTsLZ67kAAC2pBAM62GFSqBta2/J/mdSPIa0CDGmixksIgYNJOVgak4T9V24ppnf3tMbkUHf4Oluw+xapDUNaBRjSRE3D2Wu5WBqbhO1nbyi6b3VyMsckdt8iNWFIqwBDmqhpSb1VgOX7k7GR3bdIzRjSKsCQJmqasvOKsSIuFb9W6r41LtAFI7o4wZjdt6iBMaQfA29mQvR0yHug+1b2A923Rv/bfcuK3beogTCkVYB70kRPh5L7Mmw+cR3LY5ORfKsAAKCnLUVEZ0dMCHLl6FukcgxpFWBIEz1dZHKBXeczsTQmCaevlXffkkqAfh3sMSnEFW3tzdRcITUVDGkVYEgTPZ2EEIhPysHS2GTEPjD6VpCHFSaHuMHfrRm7b9FjYUirAEOaiKobfatDCzNMCnFDeFtbaLH7Fj0ChrQKMKSJqMLVnEL8sD8Z64+lo+Tf7lsuVkaYEOSKIZ0coK+jpeYKqTFhSKsAQ5qIKruVX4JVB1Pxv/g05BaVAQCsTfQwNsAZL3VtCVN9HTVXSI0BQ1oFGNJEVJOCkvtYezQdP+5Pxo3cYgCAsZ42XvRzwsuBLrAx1VdzhaTJGNIqwJAmoocpk8nx56kMLItNQkJWPgBAV0uKwc84YGKIK9ysjdVcIWkihrQKMKSJqK7kcoF9l7OxNCYJR1PvAAAkEqC3lw0mhbjhGScLNVdImoQhrQIMaSJ6FMfTbmNJdDL+uZilmObnYolJoW4IbWXN7lvEkFYFhjQRPY4rWXlYFpuMP05dR5ms/KO2ta0JXglxRf8O9tDR4oAeTyuG9GPgvbuJSJVu5Bbh5wMpWH34KgpKZQAAB3MDjA9ywXBfRxjqckCPpw1DWgW4J01EqpRbWIZfD6dhRVwKbuWXAgAsDHUw2t8Zkd2cYWmkq+YK6UlhSKsAQ5qIGkJxmQwbj1/DD/uTkZZTCADQ15FieGdHjA9yhaMlB/Ro6hjSKsCQJqKGJJMLbD93A0tjknDu+j0AgJZUgv4d7PBKsBu87Pm501QxpFWAIU1ET4IQAgeTcrA0Jgn7r9xSTA9uZY1JIa7wd+WAHk0NQ1oFGNJE9KSdu56LpTFJ2Hb2hmJAD+9/B/TozQE9mgyGtAowpIlIXdJyCvDD/mRsOHZNaUCPicGuGPwMB/Ro7BjSKsCQJiJ144AeTRNDWgUY0kSkKQpK7mPNkav46UAKB/RoAhjSKsCQJiJNU3pfjr9OKw/ooaMlKR/QI9gN7s05oEdjwJBWAYY0EWmqmgb06NXGBpNC3dCJA3poNIa0CjCkiagxOJ52G0tjkrH7wn8DenRxscSkEFd092zO7lsaiCGtAgxpImpMErPzsCwmGVseGNDD06Z8QI8B3hzQQ5MwpFWAIU1EjVFmbjF+jkvBb4fSFAN62JvpY1yQK0b4OsJIjwN6qBtDWgUY0kTUmOUWleHXQ2lYEZeKW/klAABzQx2M7toSkd2c0cxYT80VPr0Y0irAkCaipqC4TIZNJ65jeWwSUh8Y0COisyPGB7rCqRkH9HjSGNIqwJAmoqZEJhfYeT4TS2OScOZaLgBAKgH6dbDHpBBXtLU3U3OFTw+G9GOIiopCVFQUZDIZEhISGNJE1KQIIRCflIMllQb0CPKwwuQQN/i7cUCPhsaQVgHuSRNRU3c+IxfLYpLx95kMxYAeHf4d0COcA3o0GIa0CjCkiehpkX67ED/sT8a6o+lKA3pMCHLFkE4c0EPVGNIqwJAmoqdNzr8Deqx6YEAPK+P/BvQwM+CAHqrAkFYBhjQRPa0KSu5j3dF0/Lg/GRkPDOgx0s8J4zigx2NjSKsAQ5qInnZlsn8H9IhJxuWsPAAc0EMVGNIqwJAmIionxL8DekQn40jqbQAc0ONxMKRVgCFNRFTV8bTbWBKdjH8uKg/oMTnEDaGe1uy+VQcMaRVgSBMR1exKVh6WxSbjjwcG9Ghta4JJIW7o38EO2hzQo0YMaRVgSBMRPdyN3CL8tD8Fa45cVQzo4WBugAlBLhju6wQDXXbfqowhrQIMaSKiusstLMMvh1KxIi4VOQWlAAALQx1EdnNGpL8zLIx01Vyh5mBIqwBDmoio/orLZNhw/Bp+iE3G1dvlA3oY6GhhRBdHjA9yhYO5gZorVD+GtAowpImIHt19mRzbz2ViSXQSLty4BwDQlkrwnLc9Xglxg6etiZorVB+GtAowpImIHp8QAvuv3MLSmCQcTMpRTO/ZujkmhbrB19lSjdWpB0NaBRjSRESqdebaXSyNScL2c5moSJ/OLS0wKcQNPVo3h/QpGdCDIa0CDGkiooaRfDMfP+xPxu/Hr6NUVj6gRysbY7wS7IbnOtpDp4l332JIqwBDmoioYWXfK8bPcan47VAa8kruAwDszfQxLsgVI3wdYaSnreYKGwZDWgUY0kRET8a94jL8dugqfo5Lwc28EgCAuaEORvs7Y0w3Z1g2se5bDGkVYEgTET1ZxWUybDpxHctjk5CaU959S19HihG+Thgf5IIWFoZqrlA1GNIqwJAmIlIPmVxg5/ny7ltnr+cCALQU3bdc0dq2cX8mM6RVgCFNRKReQggcTMrBkugkHEi8pZjeo3VzTG7E3bcY0irAkCYi0hxnr+ViaUwStp27oei+5dPSApMbYfcthvS/0tPTMWrUKGRnZ0NbWxsfffQRhg0bVqe/ZUgTEWmelFsFWB6bjN+PX1PqvjUpxA0DvBtH9y2G9L9u3LiBrKwsdOzYEZmZmfDx8UFCQgKMjIwe+rcMaSIizVXRfevXQ2nI/7f7loO5AcYHuWC4ryMMdTW3+xZDugbe3t74+++/4ejo+NB5GdJERJovt6gMvx1Ow88HUnAr/7/Rt8Z0c8Fo/5YaOfpWXfNF7ccEYmNjMWDAANjb20MikWDLli1V5omKioKzszP09fXh5+eHI0eOPNK6jh8/DplMVqeAJiKixsHMQAdTQt1x4N0e+GxQOzhZGuJOYRm+/ScBAV/uxad/X8CN3CJ1l/lI1B7SBQUF8Pb2RlRUVLXPr1u3DjNmzMCsWbNw4sQJeHt7Izw8HNnZ2Yp5OnbsiHbt2lX5ycjIUMxz+/ZtjB49GsuXL2/w10RERE+evo4WXuraEnvfDMHCF56Bl50pCktl+OlACoLn7cPbG04jMTtf3WXWi0Yd7pZIJNi8eTMGDRqkmObn5wdfX18sWrQIACCXy+Ho6IjXXnsN7733Xp2WW1JSgl69emHChAkYNWpUrfOVlJQoHt+7dw+Ojo483E1E1AgJIRCTcBNLY5JwKPk2AEAiAXp72WByqDs6OpqrrbZGc7i7NqWlpTh+/DjCwsIU06RSKcLCwhAfH1+nZQghMGbMGPTo0aPWgAaAuXPnwszMTPHDw+JERI2XRCJBqGdzrJ3oj01TuqG3lw2EAHaez8KgqDi8sPwQYhNuQoP2VavQ6JC+desWZDIZbGxslKbb2NggMzOzTsuIi4vDunXrsGXLFnTs2BEdO3bE2bNnq5135syZyM3NVfykp6c/9msgIiL16+RkgeWjO2P3G8F43qcFtKUSxCfnYPTPR9D/+wP4+0wGZHLNC2vNvT5dRQIDAyGXy+s0r56eHvT09Bq4IiIiUhcPGxN8Ncwbb/RqhZ/2p2DNkas4n3EPr64+CedmlzEx2A1DfRygp62l7lIBaPietJWVFbS0tJCVlaU0PSsrC7a2tmqqioiIGjsHcwN8PMALB9/rgelhHjA31EFqTiHe33wWQV/uw/LYJEXfa3XS6JDW1dWFj48P9uzZo5gml8uxZ88e+Pv7q7EyIiJqCiyMdDE9rBXi3u2Bj/p7wc5MH9l5Jfh82yV0m7sHX+28jJz8kocvqIGo/XB3fn4+EhMTFY9TUlJw6tQpWFpawsnJCTNmzEBkZCQ6d+6MLl26YMGCBSgoKMDYsWMbrKaoqChERUVBJpM12DqIiEhzGOlpY1ygC0Z1bYktp65jaUwSkm8WYNG+RPx4IBnDOztiQrDrEx8qU+1dsKKjo9G9e/cq0yMjI7Fy5UoAwKJFizB//nxkZmaiY8eOWLhwIfz8/Bq8Nt5xjIjo6SSXC+y6kInF0Uk4c+2/oTIHettjcqgbPGxMHmv5vC2oCjCkiYiebkIIxCflYPEDQ2XOGuCFsQEuj7XcuuaL2g93ExERaSqJRIJu7lbo5m6FM9fuYtXBNAz3fXL30GBIV4PnpImIqLIOLczxdYT5E10nD3fXgoe7iYioITSJ24ISERE9zRjSREREGoohTUREpKEY0kRERBqKIV2NqKgoeHl5wdfXV92lEBHRU4xXd9ciNzcX5ubmSE9P59XdRESkMvfu3YOjoyPu3r0LMzOzGudjP+la5OXlAQAcHZ9cx3UiInp65OXl1RrS3JOuhVwuR0ZGBkxMTCCRSB5rWRXfmhrjXnljrh1o3PU35tqBxl1/Y64daNz1N+bagbrVL4RAXl4e7O3tIZXWfOaZe9K1kEqlaNGihUqXaWpq2ijfdEDjrh1o3PU35tqBxl1/Y64daNz1N+bagYfXX9sedAVeOEZERKShGNJEREQaiiH9hOjp6WHWrFnQ09NTdyn11phrBxp3/Y25dqBx19+Yawcad/2NuXZAtfXzwjEiIiINxT1pIiIiDcWQJiIi0lAMaSIiIg3FkH4CoqKi4OzsDH19ffj5+eHIkSPqLqlOZs+eDYlEovTTunVrdZdVo9jYWAwYMAD29vaQSCTYsmWL0vNCCHz88cews7ODgYEBwsLCcOXKFfUUW8nDah8zZkyVbdGnTx/1FFvJ3Llz4evrCxMTEzRv3hyDBg3C5cuXleYpLi7G1KlT0axZMxgbG2Po0KHIyspSU8XK6lJ/aGholfafNGmSmir+z5IlS9ChQwdFf1x/f39s375d8bwmtzvw8Po1td2r88UXX0AikWD69OmKaapof4Z0A1u3bh1mzJiBWbNm4cSJE/D29kZ4eDiys7PVXVqdtG3bFjdu3FD8HDhwQN0l1aigoADe3t6Iioqq9vl58+Zh4cKFWLp0KQ4fPgwjIyOEh4ejuLj4CVda1cNqB4A+ffoobYs1a9Y8wQprFhMTg6lTp+LQoUPYvXs3ysrK0Lt3bxQUFCjmeeONN/DXX39hw4YNiImJQUZGBoYMGaLGqv9Tl/oBYMKECUrtP2/ePDVV/J8WLVrgiy++wPHjx3Hs2DH06NEDAwcOxPnz5wFodrsDD68f0Mx2r+zo0aNYtmwZOnTooDRdJe0vqEF16dJFTJ06VfFYJpMJe3t7MXfuXDVWVTezZs0S3t7e6i7jkQAQmzdvVjyWy+XC1tZWzJ8/XzHt7t27Qk9PT6xZs0YNFdascu1CCBEZGSkGDhyolnrqKzs7WwAQMTExQojydtbR0REbNmxQzHPx4kUBQMTHx6urzBpVrl8IIUJCQsTrr7+uvqLqwcLCQvz444+Nrt0rVNQvRONo97y8POHh4SF2796tVK+q2p970g2otLQUx48fR1hYmGKaVCpFWFgY4uPj1VhZ3V25cgX29vZwdXXFiy++iKtXr6q7pEeSkpKCzMxMpW1hZmYGPz+/RrMtoqOj0bx5c3h6emLy5MnIyclRd0nVys3NBQBYWloCAI4fP46ysjKltm/dujWcnJw0su0r11/ht99+g5WVFdq1a4eZM2eisLBQHeXVSCaTYe3atSgoKIC/v3+ja/fK9VfQ9HafOnUq+vXrp9TOgOre97x3dwO6desWZDIZbGxslKbb2Njg0qVLaqqq7vz8/LBy5Up4enrixo0bmDNnDoKCgnDu3DmYmJiou7x6yczMBIBqt0XFc5qsT58+GDJkCFxcXJCUlIT3338fzz77LOLj46GlpaXu8hTkcjmmT5+OgIAAtGvXDkB52+vq6sLc3FxpXk1s++rqB4CRI0eiZcuWsLe3x5kzZ/Duu+/i8uXL2LRpkxqrLXf27Fn4+/ujuLgYxsbG2Lx5M7y8vHDq1KlG0e411Q9odrsDwNq1a3HixAkcPXq0ynOqet8zpKlGzz77rOL3Dh06wM/PDy1btsT69esxbtw4NVb29BkxYoTi9/bt26NDhw5wc3NDdHQ0evbsqcbKlE2dOhXnzp3T6GsXalNT/RMnTlT83r59e9jZ2aFnz55ISkqCm5vbky5TiaenJ06dOoXc3Fxs3LgRkZGRiImJUWtN9VFT/V5eXhrd7unp6Xj99dexe/du6OvrN9h6eLi7AVlZWUFLS6vK1XxZWVmwtbVVU1WPztzcHK1atUJiYqK6S6m3ivZuKtvC1dUVVlZWGrUtXn31Vfz999/Yt2+f0uhxtra2KC0txd27d5Xm17S2r6n+6vj5+QGARrS/rq4u3N3d4ePjg7lz58Lb2xvfffddo2n3muqvjia1+/Hjx5GdnY1OnTpBW1sb2traiImJwcKFC6GtrQ0bGxuVtD9DugHp6urCx8cHe/bsUUyTy+XYs2eP0jmXxiI/Px9JSUmws7NTdyn15uLiAltbW6Vtce/ePRw+fLhRbotr164hJydHI7aFEAKvvvoqNm/ejL1798LFxUXpeR8fH+jo6Ci1/eXLl3H16lWNaPuH1V+dU6dOAYBGtH9lcrkcJSUlGt/uNamovzqa1O49e/bE2bNncerUKcVP586d8eKLLyp+V0n7q/Y6N6ps7dq1Qk9PT6xcuVJcuHBBTJw4UZibm4vMzEx1l/ZQb775poiOjhYpKSkiLi5OhIWFCSsrK5Gdna3u0qqVl5cnTp48KU6ePCkAiG+++UacPHlSpKWlCSGE+OKLL4S5ubn4448/xJkzZ8TAgQOFi4uLKCoqUnPltdeel5cn3nrrLREfHy9SUlLEP//8Izp16iQ8PDxEcXGxuksXkydPFmZmZiI6OlrcuHFD8VNYWKiYZ9KkScLJyUns3btXHDt2TPj7+wt/f381Vv2fh9WfmJgoPvnkE3Hs2DGRkpIi/vjjD+Hq6iqCg4PVXLkQ7733noiJiREpKSnizJkz4r333hMSiUTs2rVLCKHZ7S5E7fVrcrvXpPLV6Kpof4b0E/D9998LJycnoaurK7p06SIOHTqk7pLqZPjw4cLOzk7o6uoKBwcHMXz4cJGYmKjusmq0b98+AaDKT2RkpBCivBvWRx99JGxsbISenp7o2bOnuHz5snqL/ldttRcWForevXsLa2troaOjI1q2bCkmTJigMV/0qqsbgFixYoVinqKiIjFlyhRhYWEhDA0NxeDBg8WNGzfUV/QDHlb/1atXRXBwsLC0tBR6enrC3d1dvP322yI3N1e9hQshXn75ZdGyZUuhq6srrK2tRc+ePRUBLYRmt7sQtdevye1ek8ohrYr25yhYREREGornpImIiDQUQ5qIiEhDMaSJiIg0FEOaiIhIQzGkiYiINBRDmoiISEMxpImIiDQUQ5qIiEhDMaSJGjFnZ2csWLDgsZYxe/ZsdOzYUSX11CQ1NRUSiURx72UiqhvecYyoEVi5ciWmT59eZUSdmzdvwsjICIaGho+87Pz8fJSUlKBZs2aPWWW5MWPG4O7du9iyZYtimkwmw82bN2FlZQVtbfWMkDt79mxs2bKFXxSoUeF40kSNmLW19WMvw9jYGMbGxiqopmZaWloaNTwiUWPBw91EDaykpATTpk1D8+bNoa+vj8DAQBw9elTxfHR0NCQSCbZu3YoOHTpAX18fXbt2xblz5xTPjx07Frm5uZBIJJBIJJg9ezaAqoe7JRIJli1bhv79+8PQ0BBt2rRBfHw8EhMTERoaCiMjI3Tr1g1JSUmKv6l8uLtiHQ/+ODs7AyjfIx43bhxcXFxgYGAAT09PpbF/Z8+ejVWrVuGPP/5Q/G10dHS1h7tjYmLQpUsX6Onpwc7ODu+99x7u37+veD40NBTTpk3DO++8A0tLS9ja2iped02io6PRpUsXGBkZwdzcHAEBAUhLS8PKlSsxZ84cnD59WlHXypUrAQB3797F+PHjYW1tDVNTU/To0QOnT5+u0j7Lli2Do6MjDA0NERERgdzc3FprIVIJFQ4AQkTVmDZtmrC3txfbtm0T58+fF5GRkcLCwkLk5OQIIf4bAatNmzZi165d4syZM6J///7C2dlZlJaWipKSErFgwQJhamqqGEYxLy9PCCFEy5YtxbfffqtYFwDh4OAg1q1bJy5fviwGDRoknJ2dRY8ePcSOHTvEhQsXRNeuXUWfPn0UfzNr1izh7e2tePzgcI2JiYnC3d1djBo1SgghRGlpqfj444/F0aNHRXJysvj111+FoaGhWLdunRCifMjNiIgI0adPH8UySkpKREpKigAgTp48KYQQ4tq1a8LQ0FBMmTJFXLx4UWzevFlYWVmJWbNmKeoICQkRpqamYvbs2SIhIUGsWrVKaRjGysrKyoSZmZl46623RGJiorhw4YJYuXKlSEtLE4WFheLNN98Ubdu2rTIUZVhYmBgwYIA4evSoSEhIEG+++aZo1qyZYvvMmjVLGBkZiR49eoiTJ0+KmJgY4e7uLkaOHPnobwqiOmJIEzWg/Px8oaOjI3777TfFtNLSUmFvby/mzZsnhPgvpNeuXauYJycnRxgYGCjCb8WKFcLMzKzK8qsL6Q8//FDxOD4+XgAQP/30k2LamjVrhL6+vuJx5ZCuIJfLxeDBg4WPj4/S2NCVTZ06VQwdOlTxODIyUgwcOFBpnsoh/f777wtPT08hl8sV80RFRQljY2Mhk8mEEOUhHRgYqLQcX19f8e6771ZbR05OjgAgoqOjq32+ute5f/9+YWpqWmVcbjc3N7Fs2TLF32lpaYlr164pnt++fbuQSqUaNewjNU08J03UgJKSklBWVoaAgADFNB0dHXTp0gUXL15Umtff31/xu6WlJTw9PavMUxcdOnRQ/G5jYwMAaN++vdK04uJi3Lt3D6ampjUu5/3330d8fDyOHTsGAwMDxfSoqCj8/PPPuHr1KoqKilBaWlrvq8MvXrwIf39/SCQSxbSAgADk5+fj2rVrcHJyqvJaAMDOzg7Z2dnVLtPS0hJjxoxBeHg4evXqhbCwMERERMDOzq7GOk6fPo38/PwqF80VFRUpnRJwcnKCg4OD4rG/vz/kcjkuX77Mc+3UoBjSRE2Mjo6O4veKEKxumlwur3EZv/76K7799ltER0crhdPatWvx1ltv4euvv4a/vz9MTEwwf/58HD58WNUvo0rdFbXXVveKFSswbdo07NixA+vWrcOHH36I3bt3o2vXrtXOn5+fDzs7O0RHR1d5ztzc/HFKJ1IJXjhG1IDc3Nygq6uLuLg4xbSysjIcPXoUXl5eSvMeOnRI8fudO3eQkJCANm3aAAB0dXUhk8meSM3x8fEYP348li1bViXc4uLi0K1bN0yZMgXPPPMM3N3dlfY461prxQVt4oEeoHFxcTAxMUGLFi0eq/5nnnkGM2fOxMGDB9GuXTusXr26xro6deqEzMxMaGtrw93dXenHyspKMd/Vq1eRkZGheHzo0CFIpVJ4eno+Vq1ED8OQJmpARkZGmDx5Mt5++23s2LEDFy5cwIQJE1BYWIhx48YpzfvJJ59gz549OHfuHMaMGQMrKysMGjQIQPlV3Pn5+dizZw9u3bqFwsLCBqk3MzMTgwcPxogRIxAeHo7MzExkZmbi5s2bAAAPDw8cO3YMO3fuREJCAj766COlK9Uraj1z5gwuX76MW7duoaysrMp6pkyZgvT0dLz22mu4dOkS/vjjD8yaNQszZsyAVPpoH0spKSmYOXMm4uPjkZaWhl27duHKlSuKLzrOzs5ISUnBqVOncOvWLZSUlCAsLAz+/v4YNGgQdu3ahdTUVBw8eBAffPABjh07pli2vr4+IiMjcfr0aezfvx/Tpk1DREQED3VTg2NIEzWwL774AkOHDsWoUaPQqVMnJCYmYufOnbCwsKgy3+uvvw4fHx9kZmbir7/+gq6uLgCgW7dumDRpEoYPHw5ra2vMmzevQWq9dOkSsrKysGrVKtjZ2Sl+fH19AQCvvPIKhgwZguHDh8PPzw85OTmYMmWK0jImTJgAT09PdO7cGdbW1kpHESo4ODhg27ZtOHLkCLy9vTFp0iSMGzcOH3744SPXbmhoiEuXLmHo0KFo1aoVJk6ciKlTp+KVV14BAAwdOhR9+vRB9+7dYW1tjTVr1kAikWDbtm0IDg7G2LFj0apVK4wYMQJpaWmK8/kA4O7ujiFDhqBv377o3bs3OnTogMWLFz9yrUR1xTuOEalZdHQ0unfvjjt37vA8qAbincpInbgnTUREpKEY0kRERBqKh7uJiIg0FPekiYiINBRDmoiISEMxpImIiDQUQ5qIiEhDMaSJiIg0FEOaiIhIQzGkiYiINBRDmoiISEMxpImIiDTU/wO5bhQjlku8fAAAAABJRU5ErkJggg==",
      "text/plain": [
       "<Figure size 500x300 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "hist = np.asarray(result.history)\n",
    "fig, ax = plt.subplots(figsize=(5, 3))\n",
    "ax.semilogy(hist)\n",
    "ax.set_xlabel('optimization step')\n",
    "ax.set_ylabel('loss')\n",
    "ax.set_title('Custom amplitude objective: convergence')\n",
    "fig.tight_layout()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8e126b33122a",
   "metadata": {},
   "source": [
    "## Using `objective=` with a `predict` closure\n",
    "\n",
    "When the loss *is* a `brainmass.objectives` builder over the trajectory (e.g. an\n",
    "FC fit), pass it via `objective=` and supply `predict` + `target`. The Fitter\n",
    "handles the `reg_loss` and the optimisation loop.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "890a224b56cb",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2026-06-19T07:41:31.387837Z",
     "iopub.status.busy": "2026-06-19T07:41:31.387593Z",
     "iopub.status.idle": "2026-06-19T07:41:35.086824Z",
     "shell.execute_reply": "2026-06-19T07:41:35.084810Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "fitted coupling k : 0.456  (true 1.0)\n",
      "best FC loss      : 0.0002\n"
     ]
    }
   ],
   "source": [
    "# 3-region drive so FC is non-degenerate.\n",
    "conn = brainmass.datasets.load_dataset('example_connectome')\n",
    "W3, D3 = conn.weights[:3, :3], conn.distances[:3, :3]\n",
    "brainstate.environ.set(dt=0.1 * u.ms)\n",
    "\n",
    "def make_net(k):\n",
    "    base = brainmass.HopfStep(in_size=3, a=0.3, w=0.3,\n",
    "                              init_x=braintools.init.Constant(0.3))\n",
    "    return brainmass.Network(base, conn=W3, distance=D3, speed=10 * u.mm / u.ms,\n",
    "                             coupling='diffusive', coupled_var='x', k=k)\n",
    "\n",
    "def simulate(net):\n",
    "    res = brainmass.Simulator(net, dt=0.1 * u.ms).run(\n",
    "        400 * u.ms, monitors=lambda m: m.node.x.value, transient=100 * u.ms)\n",
    "    return res['output']\n",
    "\n",
    "# Target FC from a \"ground-truth\" coupling of 1.0.\n",
    "target_signal = simulate(make_net(1.0))\n",
    "\n",
    "# Fit k with a trainable Param, FC-correlation objective.\n",
    "net = make_net(Param(0.3, fit=True))\n",
    "fitter = brainmass.Fitter(\n",
    "    net, braintools.optim.Adam(lr=0.05),\n",
    "    objective=objectives.fc_corr(as_loss=True),\n",
    "    predict=simulate,\n",
    ")\n",
    "res = fitter.fit(target=target_signal, n_steps=30)\n",
    "print(f\"fitted coupling k : {res.best_params['coupling.k']:.3f}  (true 1.0)\")\n",
    "print(f\"best FC loss      : {res.best_loss:.4f}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d8e3a5d9c8c9",
   "metadata": {},
   "source": [
    "## Recap\n",
    "\n",
    "- Objectives are `(prediction, target) -> scalar` callables; the builders in\n",
    "  `brainmass.objectives` produce them, and `combine` weights them together.\n",
    "- Write your own with `jax.numpy`, reduce to a scalar, and keep units through\n",
    "  subtractions.\n",
    "- Feed the Fitter a derived-scalar `loss_fn=` (oscillator amplitude / spectrum) or\n",
    "  an `objective=` + `predict=` (FC / FCD over the trajectory).\n",
    "\n",
    "## Next steps\n",
    "\n",
    "- {doc}`/tutorials/06_fitting_with_gradients` — the full gradient-fitting workflow.\n",
    "- {doc}`/tutorials/07_gradient_free_fitting` — the same objective, derivative-free.\n",
    "- {doc}`/howto/analyze_results` — the metrics these objectives are built on.\n",
    "- {doc}`/reference/orchestration` — the `Fitter` / `objectives` API.\n"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "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
}