{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "c382ea91",
   "metadata": {},
   "source": [
    "# Monte Carlo View: Nominal vs Robust Grating\n",
    "\n",
    "> The robust adjoint design trades a sliver of peak efficiency for tighter fabrication yield. Building on the fabrication-aware optimizer from the previous notebook, we now quantify how much that robustness actually helps under process variation.\n",
    "\n",
    "> This notebook compares the nominal adjoint design against the robustness-optimized variant using a matched Monte Carlo experiment, highlighting the yield benefits of carrying fabrication awareness into the optimization loop."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "62860f03",
   "metadata": {
    "lines_to_next_cell": 2
   },
   "outputs": [],
   "source": [
    "import json\n",
    "from pathlib import Path\n",
    "\n",
    "import autograd.numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import pandas as pd\n",
    "import tidy3d as td\n",
    "from setup import (\n",
    "    center_wavelength,\n",
    "    default_spacer_thickness,\n",
    "    get_mode_monitor_power,\n",
    "    make_simulation,\n",
    ")\n",
    "from tidy3d import web"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "c5669f17",
   "metadata": {},
   "outputs": [],
   "source": [
    "design_paths = {\n",
    "    \"nominal\": Path(\"./results\") / \"gc_adjoint_best.json\",\n",
    "    \"robust\": Path(\"./results\") / \"gc_adjoint_robust_best.json\",\n",
    "}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "3965c290",
   "metadata": {},
   "outputs": [],
   "source": [
    "def load_nominal_parameters(path):\n",
    "    \"\"\"Load a design JSON (Bayes or adjoint) into numpy-friendly fields.\"\"\"\n",
    "    data = json.loads(Path(path).read_text(encoding=\"utf-8\"))\n",
    "    return {\n",
    "        \"widths_si\": np.array(data[\"widths_si\"], dtype=float),\n",
    "        \"gaps_si\": np.array(data[\"gaps_si\"], dtype=float),\n",
    "        \"widths_sin\": np.array(data[\"widths_sin\"], dtype=float),\n",
    "        \"gaps_sin\": np.array(data[\"gaps_sin\"], dtype=float),\n",
    "        \"first_gap_si\": float(data[\"first_gap_si\"]),\n",
    "        \"first_gap_sin\": float(data[\"first_gap_sin\"]),\n",
    "        \"spacer_thickness\": default_spacer_thickness,\n",
    "    }\n",
    "\n",
    "\n",
    "def make_variation_builder(nominal):\n",
    "    \"\"\"Return a closure that maps process deltas to a tidy3d Simulation.\"\"\"\n",
    "    base_widths_si = np.array(nominal[\"widths_si\"])\n",
    "    base_gaps_si = np.array(nominal[\"gaps_si\"])\n",
    "\n",
    "    def builder(*, overlay_delta=0.0, spacer_delta=0.0, etch_bias=0.0):\n",
    "        # Etch bias widens features when positive and narrows them when\n",
    "        # negative, so widths grow with the bias while gaps shrink, mirroring\n",
    "        # the fabrication effect of over/under etching.\n",
    "        pert_widths_si = base_widths_si + etch_bias\n",
    "        pert_gaps_si = base_gaps_si - etch_bias\n",
    "\n",
    "        return make_simulation(\n",
    "            pert_widths_si,\n",
    "            pert_gaps_si,\n",
    "            nominal[\"widths_sin\"],\n",
    "            nominal[\"gaps_sin\"],\n",
    "            first_gap_si=nominal[\"first_gap_si\"] + overlay_delta,\n",
    "            first_gap_sin=nominal[\"first_gap_sin\"],\n",
    "            spacer_thickness=nominal[\"spacer_thickness\"] + spacer_delta,\n",
    "        )\n",
    "\n",
    "    return builder"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "113c27f0",
   "metadata": {},
   "source": [
    "## Shared Monte Carlo Draws\n",
    "We reuse the exact same random perturbations for both designs so any differences in the statistics stem from the geometry rather than sampling noise."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "d773fdf8",
   "metadata": {},
   "outputs": [],
   "source": [
    "sigma_overlay = 0.025  # microns\n",
    "sigma_spacer = 0.02\n",
    "sigma_widths_si = 0.01\n",
    "\n",
    "seed = 42\n",
    "num_mc_samples = 100\n",
    "\n",
    "sigma_vector = np.array([sigma_overlay, sigma_spacer, sigma_widths_si], dtype=float)\n",
    "rng = np.random.default_rng(seed)\n",
    "perturbations = rng.standard_normal(size=(num_mc_samples, sigma_vector.size)) * sigma_vector"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "a20de9b6",
   "metadata": {},
   "outputs": [],
   "source": [
    "def run_monte_carlo(nominal_params, *, label=\"design\"):\n",
    "    \"\"\"Evaluate a design under shared Monte Carlo perturbations.\"\"\"\n",
    "    builder = make_variation_builder(nominal_params)\n",
    "\n",
    "    sims = {\n",
    "        f\"{label}_nominal\": builder(),\n",
    "    }\n",
    "    for idx, (overlay_delta, spacer_delta, etch_bias) in enumerate(perturbations):\n",
    "        sims[f\"{label}_mc_{idx:03d}\"] = builder(\n",
    "            overlay_delta=overlay_delta,\n",
    "            spacer_delta=spacer_delta,\n",
    "            etch_bias=etch_bias,\n",
    "        )\n",
    "\n",
    "    batch = web.run_async(sims, verbose=False)\n",
    "    freq0 = td.C_0 / center_wavelength\n",
    "\n",
    "    nominal_value = None\n",
    "    sample_values = []\n",
    "\n",
    "    for key in sims:\n",
    "        sim_data = batch[key]\n",
    "        power_da = get_mode_monitor_power(sim_data)\n",
    "        center_power = power_da.sel(f=freq0, method=\"nearest\").item()\n",
    "        if key.endswith(\"nominal\"):\n",
    "            nominal_value = float(center_power)\n",
    "        else:\n",
    "            sample_values.append(float(center_power))\n",
    "\n",
    "    return {\n",
    "        \"nominal\": nominal_value,\n",
    "        \"samples\": np.array(sample_values, dtype=float),\n",
    "    }"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "a816af68",
   "metadata": {},
   "outputs": [],
   "source": [
    "design_results = {}\n",
    "for label, path in design_paths.items():\n",
    "    params = load_nominal_parameters(path)\n",
    "    design_results[label] = run_monte_carlo(params, label=label)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "7d575288",
   "metadata": {},
   "outputs": [],
   "source": [
    "def linear_to_loss_db(values):\n",
    "    \"\"\"Convert linear transmission to loss in dB (positive = loss).\"\"\"\n",
    "    linear = np.clip(np.array(values, dtype=float), 1e-12, None)\n",
    "    return -10.0 * np.log10(linear)\n",
    "\n",
    "\n",
    "def summarize(center_samples):\n",
    "    \"\"\"Compute summary statistics in linear and dB scales.\"\"\"\n",
    "    linear = np.array(center_samples, dtype=float)\n",
    "    stats = {\n",
    "        \"mean_linear\": np.mean(linear),\n",
    "        \"std_linear\": np.std(linear, ddof=0),\n",
    "        \"p10_linear\": np.percentile(linear, 10),\n",
    "        \"p90_linear\": np.percentile(linear, 90),\n",
    "    }\n",
    "\n",
    "    loss_db = linear_to_loss_db(linear)\n",
    "    stats.update(\n",
    "        {\n",
    "            \"mean_db\": np.mean(loss_db),\n",
    "            \"p10_db\": np.percentile(loss_db, 10),\n",
    "            \"p90_db\": np.percentile(loss_db, 90),\n",
    "        }\n",
    "    )\n",
    "    return stats"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "254e0dfb",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>mean_linear</th>\n",
       "      <th>std_linear</th>\n",
       "      <th>p10_linear</th>\n",
       "      <th>p90_linear</th>\n",
       "      <th>mean_db</th>\n",
       "      <th>p10_db</th>\n",
       "      <th>p90_db</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>nominal</th>\n",
       "      <td>0.554774</td>\n",
       "      <td>0.027361</td>\n",
       "      <td>0.517067</td>\n",
       "      <td>0.587782</td>\n",
       "      <td>2.564235</td>\n",
       "      <td>2.307842</td>\n",
       "      <td>2.864534</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>robust</th>\n",
       "      <td>0.561130</td>\n",
       "      <td>0.028289</td>\n",
       "      <td>0.522749</td>\n",
       "      <td>0.598293</td>\n",
       "      <td>2.514948</td>\n",
       "      <td>2.230860</td>\n",
       "      <td>2.817066</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "         mean_linear  std_linear  p10_linear  p90_linear   mean_db    p10_db  \\\n",
       "nominal     0.554774    0.027361    0.517067    0.587782  2.564235  2.307842   \n",
       "robust      0.561130    0.028289    0.522749    0.598293  2.514948  2.230860   \n",
       "\n",
       "           p90_db  \n",
       "nominal  2.864534  \n",
       "robust   2.817066  "
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "summary = {label: summarize(result[\"samples\"]) for label, result in design_results.items()}\n",
    "summary_df = pd.DataFrame(summary).T\n",
    "summary_df"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2d4819ee",
   "metadata": {},
   "source": [
    "## Distribution of Center-Wavelength Loss\n",
    "Both designs now face identical process draws. The plot below overlays the center wavelength loss distributions in dB. Dashed vertical lines mark the nominal (unperturbed) efficiency for each design."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "828f2601",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAhMAAAGJCAYAAAAwtrGcAAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjcsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvTLEjVAAAAAlwSFlzAAAPYQAAD2EBqD+naQAAbzdJREFUeJzt3Xd4FFXbBvB7tqY3SJVAQuhNlBp6CVVBiiCIQhAQNUgTFPSTJhoUEUQpwqvAqyJKERVpghSp0hFBagJIh/Se3TnfH5B5s2QTdrOb7AbuH9deZM+cmXnm7MzsszNnZiQhhAARERFRMakcHQARERGVbUwmiIiIyCZMJoiIiMgmTCaIiIjIJkwmiIiIyCZMJoiIiMgmTCaIiIjIJkwmiIiIyCZMJoiIiMgmTCaoTAoLC0N0dLSjw3gosC0Lio+PhyRJ+Pjjjx0dikXatGmDNm3aODqMUjVlyhRIkoTbt2+X+Ly2b98OSZKwffv2Ep9XWcVkooQtXboUkiRBkiTs2rWrwHAhBEJDQyFJEp5++ukSjWXPnj2YMmUKkpKSSmT658+fx/Dhw1G5cmW4uLjAy8sLzZs3x6efforMzMwSmScRAKxfvx5TpkxxdBiUz8mTJzFlyhTEx8c7OhSLzZ8/H0uXLnV0GGUSk4lS4uLiguXLlxco37FjB/7991/o9foSj2HPnj2YOnVqiSQTv/76K+rWrYsffvgB3bp1w2effYbY2FhUrFgR48ePx6hRo+w+T7KP06dPY/HixY4Owybr16/H1KlTHR0G5XPy5ElMnTr1oUgmWrVqhczMTLRq1cohcZUFGkcH8Kjo2rUrVq5ciblz50Kj+V+zL1++HA0aNCiVQ3UlJS4uDv369UOlSpXw+++/Izg4WBkWExODc+fO4ddff7V5PkIIZGVlwdXV1eZpPeryt2VpJLJkXnp6Otzd3R0dhl1lZWVBp9OV2PQd0WYqlQouLi6lOs+yhkcmSkn//v1x584d/Pbbb0pZTk4OVq1aheeff97sOOnp6XjjjTcQGhoKvV6P6tWr4+OPP8b9D3qVJAkjRozA2rVrUadOHej1etSuXRsbN25U6kyZMgXjx48HAISHhyunXvL/avjmm2/QoEEDuLq6ws/PD/369cPly5cfuGwfffQR0tLS8OWXX5okEnmqVKlicmRiyZIlaNeuHQICAqDX61GrVi0sWLCgwHhhYWF4+umnsWnTJjRs2BCurq744osvCo3jwoUL6NOnD/z8/ODm5oamTZtalcR88803aNy4Mdzc3ODr64tWrVph8+bNJnXmz5+P2rVrQ6/XIyQkBDExMQWO9LRp0wZ16tTB8ePH0bp1a7i5uaFKlSpYtWoVcO9oVJMmTeDq6orq1atjy5YtJuPnnQv+559/0LdvX3h5eaFcuXIYNWoUsrKyTOraoy3v7zORm5uLqVOnomrVqnBxcUG5cuXQokULk3UXAH7//Xe0bNkS7u7u8PHxwTPPPINTp06ZXZZz584hOjoaPj4+8Pb2xuDBg5GRkfHAz+SPP/5Anz59ULFiRej1eoSGhmLMmDEmp82io6Mxb9484N62kPcqysGDB9GpUyeUL18erq6uCA8Px0svvWS27qJFixAREQG9Xo9GjRrhwIEDJsOPHz+O6Oho5fReUFAQXnrpJdy5c8dsW5w8eRLPP/88fH190aJFC2W4pdtfXjyurq5o3Lgx/vjjjwe2Y568fcW3336L6tWrw8XFBQ0aNMDOnTsL1L1y5QpeeuklBAYGKvuUr776yqROXl+CFStW4P/+7//w2GOPwc3NDXPnzkWfPn0AAG3btlU+k7w+B5IkmT0tdf+6mHeaeMeOHXjttdcQEBCAChUqmIxz+/Ztu2wnYWFh+Pvvv7Fjxw4l3rx+KIX1mVi5cqXymZUvXx4vvPACrly5YlInOjoaHh4euHLlCnr06AEPDw/4+/tj3LhxMBqNJnVXrFiBBg0awNPTE15eXqhbty4+/fRTs5+ls+GRiVISFhaGyMhIfPfdd+jSpQsAYMOGDUhOTka/fv0wd+5ck/pCCHTv3h3btm3DkCFDUL9+fWzatAnjx4/HlStXMHv2bJP6u3btwpo1a/Daa6/B09MTc+fORe/evXHp0iWUK1cOvXr1wpkzZ/Ddd99h9uzZKF++PADA398fAPD+++/j3XffRd++fTF06FDcunULn332GVq1aoUjR47Ax8en0GX75ZdfULlyZTRr1syitliwYAFq166N7t27Q6PR4JdffsFrr70GWZYRExNjUvf06dPo378/hg8fjmHDhqF69epmp3njxg00a9YMGRkZGDlyJMqVK4dly5ahe/fuWLVqFXr27FlkTFOnTsWUKVPQrFkzTJs2DTqdDvv378fvv/+Ojh07Ave+DKZOnYqoqCi8+uqrOH36NBYsWIADBw5g9+7d0Gq1yvQSExPx9NNPo1+/fujTpw8WLFiAfv364dtvv8Xo0aPxyiuv4Pnnn8fMmTPx7LPP4vLly/D09DSJqW/fvggLC0NsbCz27duHuXPnIjExEf/9739LtC2nTJmC2NhYDB06FI0bN0ZKSgoOHjyIw4cPo0OHDgCALVu2oEuXLqhcuTKmTJmCzMxMfPbZZ2jevDkOHz6MsLCwAssSHh6O2NhYHD58GP/5z38QEBCADz/8sMjPZeXKlcjIyMCrr76KcuXK4c8//8Rnn32Gf//9FytXrgQADB8+HFevXsVvv/2Gr7/+usjpAcDNmzfRsWNH+Pv7Y8KECfDx8UF8fDzWrFlToO7y5cuRmpqK4cOHQ5IkfPTRR+jVqxcuXLigfN6//fYbLly4gMGDByMoKAh///03Fi1ahL///hv79u0rkNj06dMHVatWxQcffKD8MLB0+/vyyy8xfPhwNGvWDKNHj8aFCxfQvXt3+Pn5ITQ09IHLjnvJ7Pfff4+RI0dCr9dj/vz56Ny5M/7880/UqVMHuLc9NW3aVEk+/P39sWHDBgwZMgQpKSkYPXq0yTTfe+896HQ6jBs3DtnZ2ejYsSNGjhyJuXPn4u2330bNmjUBQPnfWq+99hr8/f0xadIkpKenmwyz13YyZ84cvP766/Dw8MA777wDAAgMDCw0pqVLl2Lw4MFo1KgRYmNjcePGDXz66afYvXt3gX2m0WhEp06d0KRJE3z88cfYsmULZs2ahYiICLz66qvAvfWof//+aN++vbJdnDp1Crt37y4bp4kFlaglS5YIAOLAgQPi888/F56eniIjI0MIIUSfPn1E27ZthRBCVKpUSTz11FPKeGvXrhUAxPTp002m9+yzzwpJksS5c+eUMgBCp9OZlB07dkwAEJ999plSNnPmTAFAxMXFmUwzPj5eqNVq8f7775uU//XXX0Kj0RQozy85OVkAEM8884zFbZK3/Pl16tRJVK5c2aSsUqVKAoDYuHFjgfqVKlUSgwYNUt6PHj1aABB//PGHUpaamirCw8NFWFiYMBqNhcZz9uxZoVKpRM+ePQvUk2VZCCHEzZs3hU6nEx07djSp8/nnnwsA4quvvlLKWrduLQCI5cuXK2X//POPACBUKpXYt2+fUr5p0yYBQCxZskQpmzx5sgAgunfvbhLLa6+9JgCIY8eOKWUl0ZaPP/64ybpoTv369UVAQIC4c+eOUnbs2DGhUqnEwIEDCyzLSy+9ZDJ+z549Rbly5YqcR2HLFxsbKyRJEhcvXlTKYmJihKW7sx9//FHZJgsTFxcnAIhy5cqJhIQEpfynn34SAMQvv/xSZIzfffedACB27typlOW1Rf/+/U3qWrr95eTkiICAAFG/fn2RnZ2t1Fu0aJEAIFq3bv3AZQcgAIiDBw8qZRcvXhQuLi6iZ8+eStmQIUNEcHCwuH37tsn4/fr1E97e3soyb9u2TQAQlStXLtAOK1euFADEtm3bzMYxefLkAuX3r4t5+88WLVoIg8FgUrcktpPatWubbce85cxblrzPok6dOiIzM1Opt27dOgFATJo0SSkbNGiQACCmTZtmMs0nnnhCNGjQQHk/atQo4eXlVWA5ywqe5ihFffv2RWZmJtatW4fU1FSsW7eu0FMc69evh1qtxsiRI03K33jjDQghsGHDBpPyqKgoREREKO/r1asHLy8vXLhw4YFxrVmzBrIso2/fvrh9+7byCgoKQtWqVbFt27ZCx01JSQGAAr+qi5K/z0NycjJu376N1q1b48KFC0hOTjapGx4ejk6dOj1wmuvXr0fjxo1NDht7eHjg5ZdfRnx8PE6ePFnouGvXroUsy5g0aRJUKtNNIu9X5ZYtW5CTk4PRo0eb1Bk2bBi8vLwKnE7x8PBAv379lPfVq1eHj48PatasiSZNmijleX+b+5zuP7Lw+uuvK8uapyTa0sfHB3///TfOnj1rdvi1a9dw9OhRREdHw8/PTymvV68eOnToYBJfnldeecXkfcuWLXHnzh1l/SlM/uVLT0/H7du30axZMwghcOTIkQcuizl5vxjXrVuH3NzcIus+99xz8PX1NYkb931e+WPMysrC7du30bRpUwDA4cOHC0zz/rawdPs7ePAgbt68iVdeecWkT0J0dDS8vb0tXv7IyEg0aNBAeV+xYkU888wz2LRpE4xGI4QQWL16Nbp16wYhhElMnTp1QnJycoHlGjRoUIn2ZRo2bBjUarXZYfbeTiyR91m89tprJn0pnnrqKdSoUcPs6VVz20D+9cjHxwfp6ekFTieWFUwmSpG/vz+ioqKwfPlyrFmzBkajEc8++6zZuhcvXkRISEiBL+m8w4QXL140Ka9YsWKBafj6+iIxMfGBcZ09exZCCFStWhX+/v4mr1OnTuHmzZuFjuvl5QUASE1NfeB88uzevRtRUVHKuXZ/f3+8/fbbwL0NPb/w8HCLpnnx4kWzh+0La6/8zp8/D5VKhVq1ahU5fdxLCvLT6XSoXLlygelXqFChwOFtb2/vAoei874EzH1OVatWNXkfEREBlUpl0s+lJNpy2rRpSEpKQrVq1VC3bl2MHz8ex48ff2Bb4F573759u8Ch6PvXz7wv6Aetn5cuXVKSlrxzza1btza7fJZq3bo1evfujalTp6J8+fJ45plnsGTJEmRnZxeoa0ncCQkJGDVqFAIDA+Hq6gp/f3+lrc3FeP/nYOn2l9fu968XWq0WlStXtnj57x8fAKpVq4aMjAzcunULt27dQlJSEhYtWlQgnsGDBwP3ThUVtUz2VtT07b2dWKKobaBGjRoF9gcuLi7KKeU89++fX3vtNVSrVg1dunRBhQoV8NJLL5n0e3N27DNRyp5//nkMGzYM169fR5cuXYrsi2CNwrL2+ztrmiPLMiRJwoYNG8xOx8PDo9Bxvby8EBISghMnTlgU5/nz59G+fXvUqFEDn3zyCUJDQ6HT6bB+/XrMnj0bsiyb1C+rV24U9nnY8jndn5yUVFu2atUK58+fx08//YTNmzfjP//5D2bPno2FCxdi6NChFk3jfsVZbqPRiA4dOiAhIQFvvfUWatSoAXd3d1y5cgXR0dEFls9SkiRh1apV2LdvH3755Rds2rQJL730EmbNmoV9+/aZrO+WxN23b1/s2bMH48ePR/369eHh4QFZltG5c2ezMd7/Odiy/ZWEvJhfeOEFDBo0yGydevXqmby313Z6f4fE4kzf1u2kJBS2HuUXEBCAo0ePYtOmTdiwYQM2bNiAJUuWYODAgVi2bFmJx2grJhOlrGfPnhg+fDj27duH77//vtB6lSpVwpYtW5CammpydOKff/5RhlursB7uEREREEIgPDwc1apVs3q6Tz/9NBYtWoS9e/ciMjKyyLq//PILsrOz8fPPP5v86ivqVIolKlWqhNOnTxcot6S9IiIiIMsyTp48ifr16xc6fdzrxJj/V2BOTg7i4uIQFRVlU/zmnD171uQX2blz5yDLstK5saTaEgD8/PwwePBgDB48GGlpaWjVqhWmTJmCoUOHmrTF/f755x+UL1/eLpfu/fXXXzhz5gyWLVuGgQMHKuXmDgM/6OoNc5o2bYqmTZvi/fffx/LlyzFgwACsWLHCqoQpMTERW7duxdSpUzFp0iSlvLBTROZYuv3ltfvZs2fRrl07pTw3NxdxcXF4/PHHLZqfudjOnDkDNzc35dezp6cnjEajTet1UZ+Jr69vgaugcnJycO3aNavnY8/txNL1KP82kP+zyCsrzv4Z9450duvWDd26dYMsy3jttdfwxRdf4N1330WVKlWKNc3SwtMcpczDwwMLFizAlClT0K1bt0Lrde3aFUajEZ9//rlJ+ezZsyFJknJFiDXydvD3b8S9evWCWq3G1KlTC/xSFEIUuMTtfm+++Sbc3d0xdOhQ3Lhxo8Dw8+fPK5c35WXo+eeTnJyMJUuWWL08+XXt2hV//vkn9u7dq5Slp6dj0aJFCAsLK/IURo8ePaBSqTBt2rQCv1Ly4oyKioJOp8PcuXNNYv/yyy+RnJyMp556yqb4zcm73DHPZ599BgDKZ19SbXn/5+3h4YEqVaoopwGCg4NRv359LFu2zGRdOnHiBDZv3oyuXbvaNP885pZPCGH2UrnC1m1zEhMTC6zneUmkuVMd1saIe1cGWMrS7a9hw4bw9/fHwoULkZOTo9RZunSpVTei27t3r0mfh8uXL+Onn35Cx44doVaroVar0bt3b6xevdrsEcdbt25ZNJ+iPpOIiIgCl6MuWrSo0CMTRbHnduLu7m5RWzZs2BABAQFYuHChyTqzYcMGnDp1qlj7g/u3O5VKpRwBsna9dAQemXCAwg4d5tetWze0bdsW77zzDuLj4/H4449j8+bN+OmnnzB69GiTzpaWyut09c4776Bfv37QarXo1q0bIiIiMH36dEycOBHx8fHo0aMHPD09ERcXhx9//BEvv/wyxo0bV+h0IyIisHz5cjz33HOoWbMmBg4ciDp16iAnJwd79uzBypUrlWvHO3bsqGTfw4cPR1paGhYvXoyAgIBi/SrJM2HCBOWy25EjR8LPzw/Lli1DXFwcVq9eXaBjZX5VqlTBO++8g/feew8tW7ZEr169oNfrceDAAYSEhCA2Nhb+/v6YOHEipk6dis6dO6N79+44ffo05s+fj0aNGuGFF14oduyFiYuLQ/fu3dG5c2fs3bsX33zzDZ5//nnlF2hJtWWtWrXQpk0bNGjQAH5+fjh48CBWrVqFESNGKHVmzpyJLl26IDIyEkOGDFEuDfX29rbbba1r1KiBiIgIjBs3DleuXIGXlxdWr15ttp9F3ro9cuRIdOrUCWq12qQDbH7Lli3D/Pnz0bNnT0RERCA1NRWLFy+Gl5eX1YmQl5cXWrVqhY8++gi5ubl47LHHsHnzZsTFxVk8DUu3P61Wi+nTp2P48OFo164dnnvuOcTFxWHJkiVW9ZmoU6cOOnXqZHJpKO5dHp1nxowZ2LZtG5o0aYJhw4ahVq1aSEhIwOHDh7FlyxYkJCQ8cD7169eHWq3Ghx9+iOTkZOj1euVeD0OHDsUrr7yC3r17o0OHDjh27Bg2bdqkXLJuDXtuJw0aNMCCBQswffp0VKlSBQEBAQWOPOBeP5UPP/wQgwcPRuvWrdG/f3/l0tCwsDCMGTPG6uUYOnQoEhIS0K5dO1SoUAEXL17EZ599hvr16xf7ktpS5ejLSR52+S8NLcr9l4aKe5c2jhkzRoSEhAitViuqVq0qZs6cqVyumAeAiImJMTvN/JdZCSHEe++9Jx577DGhUqkKXCa6evVq0aJFC+Hu7i7c3d1FjRo1RExMjDh9+rRFy3rmzBkxbNgwERYWJnQ6nfD09BTNmzcXn332mcjKylLq/fzzz6JevXrCxcVFhIWFiQ8//FB89dVXBeIx1yZFLdv58+fFs88+K3x8fISLi4to3LixWLdunUWxCyHEV199JZ544gmh1+uFr6+vaN26tfjtt99M6nz++eeiRo0aQqvVisDAQPHqq6+KxMREkzqtW7cWtWvXNhuzueW5//PLu+Tt5MmT4tlnnxWenp7C19dXjBgxwuQyNFFCbTl9+nTRuHFj4ePjI1xdXUWNGjXE+++/L3JyckzG27Jli2jevLlwdXUVXl5eolu3buLkyZMmdfKW5datWybledvF/Zcp3+/kyZMiKipKeHh4iPLly4thw4Yplz3nv5zWYDCI119/Xfj7+wtJkoq8TPTw4cOif//+omLFikKv14uAgADx9NNPm1wumXdp6MyZMwuMf/9ljf/++6/o2bOn8PHxEd7e3qJPnz7i6tWrBeoV1hZ5LN3+5s+fL8LDw4VerxcNGzYUO3fuFK1bt7b40tCYmBjxzTffiKpVqwq9Xi+eeOIJs5dv3rhxQ8TExIjQ0FCh1WpFUFCQaN++vVi0aJFSJ++SyZUrV5qd3+LFi0XlypWFWq02ubTSaDSKt956S5QvX164ubmJTp06iXPnzhV6aai5/WdJbCfXr18XTz31lPD09DS53Pb+S0PzfP/998o+w8/PTwwYMED8+++/JnUGDRok3N3dC40/z6pVq0THjh1FQECA0Ol0omLFimL48OHi2rVrZtvW2UjCkp5fRFRq8m6OdevWrWL9UiMqjCRJiImJKXD6lMhW7DNBRERENmEyQURERDZhMkFEREQ2YZ8JIiIisgmPTBAREZFNmEwQERGRTR76m1bJsoyrV6/C09OzWLfcJSIielQJIZCamoqQkJAib/7n0JtWzZ8/X9StW1d4enoKT09P0bRpU7F+/XpleGZmpnjttdeEn5+fcHd3F7169RLXr1+3ah6XL18WAPjiiy+++OKLr2K+Ll++XOR3rUM7YP7yyy9Qq9WoWrUqhBBYtmwZZs6ciSNHjqB27dp49dVX8euvv2Lp0qXw9vbGiBEjoFKpsHv3bovnkZycDB8fH1y+fFl5XPbDQJZlJCYmwtfXt+hskdhWVmBbWY5tZTm2leWcra1SUlIQGhqKpKQkeHt7F1rP6a7m8PPzw8yZM/Hss8/C398fy5cvx7PPPgvceyJhzZo1sXfvXjRt2tSi6aWkpMDb2xvJyckPXTKRkJAAPz8/p1jhnBnbynJsK8uxrSzHtrKcs7WVpd+hTtNnwmg0YuXKlUhPT0dkZCQOHTqE3Nxck0fg1qhRAxUrViwymcjOzjZ5wlpKSgpw7wMqjefWlxZZliGEeKiWqaQUt61O3jmJXDkXWpUWtcr976mjf11JRq5RhlatQt3HCs/UyyKuV5ZjW1mObWU5Z2srS+NweDLx119/ITIyEllZWfDw8MCPP/6IWrVq4ejRo9DpdPDx8TGpHxgYiOvXrxc6vdjYWJOn3+VJTEyEwWAokWVwBFmWkZqaCiGEU2Svzqy4bfX676/jdtZtlHcpj5XtVyrlw5Ydxs20XAR4aLH+lSdLKGrH4HplObaV5dhWlnO2tkpNTbWonsOTierVq+Po0aNITk7GqlWrMGjQIOzYsaPY05s4cSLGjh2rvM873+Pr6/vQneaQJMlpzqs5s+K2VV5dlUoFPz+/B5Y/DLheWY5tZTm2leWcra00GsvSBIcnEzqdDlWqVAHuPUv+wIED+PTTT/Hcc88hJycHSUlJJkcnbty4gaCgoEKnp9frodfrC5SrVCqn+GDsSZKkh3K5SoKtbWU6nqT8/zC2Pdcry5VUWwkhYDAYYDQa7TpdR5FlGQaDATk5OVyvHqC020qtVkOj0RR66wRLY3B4MnE/WZaRnZ2NBg0aQKvVYuvWrejduzcA4PTp07h06RIiIyMdHSYRUYnIycnBtWvXkJGR4ehQ7CavD0BiYiLv9/MAjmgrNzc3BAcHQ6fTFXsaDk0mJk6ciC5duqBixYpITU3F8uXLsX37dmzatAne3t4YMmQIxo4dCz8/P3h5eeH1119HZGSkxVdyEBGVJbIsIy4uDmq1GiEhIdDpdA/Fl2/ekZaifgHTXaXZVkII5OTk4NatW4iLi0PVqlWLfTTEocnEzZs3MXDgQFy7dg3e3t6oV68eNm3ahA4dOgAAZs+eDZVKhd69eyM7OxudOnXC/PnzHRkyEVGJycnJgSzLCA0NhZubm6PDsRsmE5Yr7bZydXWFVqvFxYsXkZOTAxcXl2JNx6HJxJdfflnkcBcXF8ybNw/z5s0rtZiIiByN/QqoNNljfeMaS0RERDZhMkFEREQ2YTJBRERl3tKlSwvc5NBW8fHxkCQJR48etWq8sLAwzJkzxy4xlMRylQSnuzSUCACSM3OQmmWfO5YKISMtNQuZqgxIkuX584I2KyCEgCRJ+Dfxf5fp/felRhD37jaRvzyPp4sG3q7Fv8SKyBx7bhMPUpx1ODo6GsuWLUNsbCwmTJiglK9duxY9e/ZETk5OCUT6P8899xy6du1aovNwhLKyXEwmyCmlZhnw89GrSLHHzlMIuIlMZEipQLF7R9+xqJaXiwbd64cwmSC7s+s2UQRb1mEXFxd8+OGHGD58OHx9fUskvsK4urrC1dW1VOdZGsrKcvE0BzmtlCwDkjNzbX9l5SI924DkLDtM6wGvkt7R06PNbttECa3DUVFRCAoKQmxsbJH1Vq9ejdq1a0Ov1yMsLAyzZs0yGR4WFobp06dj4MCB8PDwQKVKlfDzzz/j1q1beOaZZ+Dh4YF69erh4MGDyjj3nw6YMmUK6tevj6+//hphYWHw9vZGv379TJ41sXHjRrRo0QI+Pj4oV64cnn76aZw/f96qZb558ya6desGV1dXhIeH49tvvy1QJykpCUOHDoW/vz+8vLzQrl07HDt2TBl+7NgxtG3bFp6envD29kaTJk2UZTN3mmP69OkICAiAp6cnhg4digkTJqB+/frK8OjoaPTo0QMff/wxgoODUa5cOcTExCA3N9eqZbMGkwkiIrILtVqNDz74AJ999hn+/fdfs3UOHTqEvn37ol+/fvjrr78wZcoUvPvuu1i6dKlJvdmzZ6N58+Y4cuQInnrqKbz44osYOHAgXnjhBRw+fBgREREYOHAghBCFxnP+/HmsXbsW69atw7p167Bjxw7MmDFDGZ6eno6xY8fi4MGD2Lp1K1QqFXr27GnVEzujo6Nx+fJlbNu2DatWrcL8+fNx8+ZNkzp9+vTBzZs3sWHDBhw6dAhPPvkk2rdvj4SEBADAgAEDUKFCBRw4cAAHDx7E+PHjodVqzc7v22+/xfvvv48PP/wQhw4dQsWKFbFgwYIC9bZt24bz589j27ZtWLZsGZYuXVqgje2JpzmICnEh51fkikxoJVdU1j2llP9zLUV5BHmN4Ifn4XFE9tCzZ0/Ur18fkydPNnsvoU8++QTt27fHu+++CwCoVq0aTp48iZkzZyI6Olqp17VrVwwfPhwAMGnSJCxYsACNGjVCnz59AABvvfUWIiMji3xekyzLWLp0KTw9PQEAL774IrZu3Yr3338fAJRHNeT56quv4O/vj5MnT6JOnToPXNYzZ85gw4YN+PPPP9GoUSPg3v2TatasqdTZtWsX/vzzT9y8eVN5btTHH3+MtWvXYtWqVXj55Zdx6dIljB8/HjVq1IAQAuHh4YU+YOuzzz7DkCFDMHjwYKVtNm/ejLS0NJN6vr6++Pzzz6FWq1GjRg089dRT2Lp1K4YNG/bA5SoOHpkgKsSFnA04m7MGF3I2mJT/cz0VJ66m4J/rlj2al+hR8+GHH2LZsmU4depUgWGnTp1C8+bNTcqaN2+Os2fPmjzYrF69esrfgYGBAIC6desWKLv/KEB+YWFhSiIBAMHBwSb1z549i/79+6Ny5crw8vJCWFgYAODSpUsWLeepU6eg0WjQoEEDpaxGjRompyWOHTuGtLQ0lCtXDh4eHsorLi5OOaUyduxYDB06FFFRUZgxY0aRp1pOnz6Nxo0bm5Td/x4AateuDbVaXeiy2xuTCSIisqtWrVqhU6dOmDhxYrGnkf8wf95tpc2VFXVK4v5TBZIkmdTv1q0bEhISsHjxYuzfvx/79+8H7t3W3F7S0tIQHByMo0ePmrxOnz6N8ePHA/f6d/z999946qmnsG3bNjz++OP48ccfbZrvg5bd3niag4iI7G7GjBmoX78+qlevblJes2ZN7N6926Rs9+7dqFatmskv6ZJ2584dnD59GosXL0bLli2Be6ckrFGjRg0YDAYcOnRIOc1x+vRpJCUlKXWefPJJXL9+HRqNRjnyYU61atVQrVo1jB49Gv369cPSpUvRq1evAvWqV6+OAwcOYODAgUrZgQMHrIq7JDCZICIqI7xcSn6Xba951K1bFwMGDMDcuXNNyt944w00atQI7733Hp577jns3bsXn3/+eak/xNHX1xflypXDokWLEBwcjEuXLpncH8MS1atXR+fOnTF8+HAsWLAAGo0Go0ePNrmUMyoqCpGRkejRowc++ugjVKtWDVevXsWvv/6Knj17onbt2hg/fjyeffZZhIeH4/Llyzh06JDZRAIAXn/9dQwbNgwNGzZEs2bN8P333+P48eOoXLmyzW1iCyYTRERlgOe9+z+U1rzsYdq0afj+++9Nyp588kn88MMPmDRpEt577z0EBwdj2rRpJp0vS4NKpcKKFSswcuRI1KlTB9WrV8fcuXPRpk0bq6azZMkSDB06FK1bt0ZgYCCmT5+udC7FvdML69evxzvvvIPBgwfj1q1bCAoKQqtWrRAYGAi1Wo07d+5g4MCBuHHjBsqXL48ePXpg6tSpZuc3YMAAXLhwAePGjUNWVhb69u2L6Oho/Pnnnza3iS0kUdR1NQ+BlJQUeHt7Izk5GV5eD0/Pe1mWkZCQAD8/v4fyCYP/Jmbgm32XkJxpj+uiBbyQhRS43LtvpWW2pI1AlkiAi+SHKI/PlfK1R64gM9cIV60aPZ54zGQcb1ctXmhaERV8y+bjox/29cqeSqKtsrKyEBcXh/Dw8GI/CtoZ8RHklitOW3Xo0AFBQUH4+uuvizXPotY7S79DeWSCiIiojMjIyMDChQvRqVMnqNVqfPfdd9iyZQt+++03h8bFZIKIiKiMyDtt8v777yMrKwvVq1fH6tWrERUV5dC4mEwQERGVEa6urtiyZYujwyiAyQRRIbxVYXARftBLpucJ/dx1yMo1wkVbepexERE5MyYTRIVo5DbObHmrav6lHgsRkTNjd20iIiKyCZMJIiIisgmTCSIiIrIJ+0wQFeJAxsfIFinQS14m/Sd2nrmldMBk/wkiIiYTRIVKluOVO2Dml5Ceo9wBk6i0pGSnIC03rVTm5aH1gJe+9O8YvH37drRt2xaJiYkmj/F2FtHR0UhKSsLatWttnlZ8fDzCw8Nx5MgR1K9f3y7xORKTCSKiMiAtNw3r49YjNSe1ROfjqfNE1/CuVicT0dHRWLZsGQBAo9GgQoUK6NOnD6ZNmwa9Xl9C0VpvypQpWLt2LY4ePerQOEJDQ3Ht2jWUL1/eoXHYC5MJIqIyIjUnFSk5KY4Oo1CdO3fGkiVLkJubi0OHDmHQoEGQJAkzZsxwdGhOR61WIygoyNFh2A07YBIRkV3o9XoEBQUhNDQUPXr0QFRUlMkzI7KzszFy5EgEBATAxcUFLVq0wIEDBwpMZ/fu3ahXrx5cXFzQtGlTnDhxQhk2ZcqUAqcF5syZg7CwMOX99u3b0bhxY7i7u8PHxwfNmzfHxYsXsXTpUkydOhXHjh2DJEmQJAlLly41uyxGoxFjx46Fj48PypUrhzfffBP3PxdTlmXExsYiPDwcrq6uePzxx7Fq1SpleGJiIgYMGAB/f3+4urqiatWqWLJkCXDvNIckSSZHSH7++WdUq1YNnp6eaNeuHZYtWwZJkpCUlAQAWLp0KXx8fLBp0ybUrFkTHh4e6Ny5M65du2bV51QSmEwQEZHdnThxAnv27IFOp1PK3nzzTaxevRrLli3D4cOHUaVKFXTq1AkJCQkm444fPx6zZs3CgQMH4O/vj27duiE317InCBsMBvTo0QOtW7fG8ePHsXfvXrz88suQJAnPPfcc3njjDdSuXRvXrl3DtWvX8Nxzz5mdzqxZs7B06VJ89dVX2LVrFxISEvDjjz+a1ImNjcV///tfLFy4EH///TfGjBmDF154ATt27AAAvPvuuzh58iQ2bNiAU6dOYcGCBYWe1oiLi8Ozzz6LZ555BgcPHsTLL7+Md955p0C9jIwMfPzxx/j666+xc+dOXLp0CePGmb/BXmniaQ4iIrKLdevWwcPDAwaDAdnZ2VCpVPj8888BAOnp6Vi4cCGWLl2KLl26AAAWL16M3377DV9++SXGjx+vTGfy5Mno0KEDAGDZsmWoUKECfvzxR/Tt2/eBMaSkpCA5ORlPP/00IiIiAAA1a9ZUhnt4eECj0TzwFMOcOXMwceJE9OrVCwCwcOFCbNq0SRmenZ2NDz74AFu2bEFkZCQAoHLlyti1axe++OILtG7dGpcuXcITTzyBhg0bAoDJ0ZP7ffHFF6hevTpmzpwJg8GA2rVr4++//8b7779vUi83NxcLFy5Ulm3EiBGYNm3aA9ulpDGZICIiu2jbti0WLFiA9PR0zJ49GxqNBr1794YQAufPn0dubi6aN2+u1NdqtWjcuDFOnTplMp28L2cA8PPzQ/Xq1QvUKYyfnx+io6PRqVMndOjQAVFRUejbty+Cg4MtXo7k5GRcu3YNTZo0Uco0Gg0aNmyonOo4d+4cMjIylKQnT05ODp544gkAwKuvvorevXvj8OHD6NixI3r06IFmzZqZnefp06fRqFEjk7LGjRsXqOfm5qYkEgAQHByMmzdvWrxsJYWnOYiIyC7c3d1RpUoVPP744/jqq6+wf/9+fPnll3adh0qlKtB34f5TIEuWLMHevXvRrFkzfP/996hWrRr27dtn1zjS0u5epvvrr7/i6NGjyuvkyZNKv4kuXbrg4sWLGDNmDK5evYr27dvbfEpCq9WavJckqUB7OAKTCSIisjuVSoW3334b//d//4fMzExERERAp9Nh9+7dSp3c3FwcOHAAtWrVMhk3/xd/YmIizpw5o5yq8Pf3x/Xr102+QM1d5vnEE09g4sSJ2LNnD+rUqYPly5cDAHQ6HYxGY5Gxe3t7Izg4GPv371fKDAYDDh06pLyvVasW9Ho9Ll26hCpVqpi8QkNDlXr+/v4YNGgQvvnmG8yZMweLFi0yO8/q1avj4MGDJmXmOqc6K57mICpEZV0X5IpMaCVXk/IaQZ7INcrQqpmLU+ny1HmWqXn06dMH48ePx7x58zB69Gi88sorGD9+PPz8/FCxYkV89NFHyMjIwJAhQ0zGmzZtGsqVK4fAwEC88847KF++PHr06AEAaNOmDW7duoWPPvoIzz77LDZu3IgNGzbAy+vufTHi4uKwaNEidO/eHSEhITh9+jTOnj2LgQMHAvf6LcTFxeHo0aOoUKECPD09zd4HY9SoUZgxYwaqVq2KGjVq4JNPPlGuqgAAT09PjBs3DmPGjIEsy2jRogWSk5Oxe/dueHl5YdCgQZg0aRIaNGiA2rVrIzs7G+vWrTPpv5Hf8OHD8cknn+Ctt97CoEGDcOLECeVKE0mS7PaZlBQmE0SFqKx7ymx5jeDSvzMgkYfWA13Du5bavOxBo9FgxIgRmDlzJoYNG4YZM2ZACIEXX3wRqampaNiwITZt2gRfX1+T8WbMmIFRo0bh7NmzqF+/Pn755RflqpCaNWti/vz5+OCDD/Dee++hd+/eGDdunPKL383NDf/88w+WLVuGO3fuIDg4GDExMRg+fDgAoHfv3lizZg3atm2LpKQkLFmyBNHR0QVif+ONN3Dt2jUMGjQIKpUKL730Enr27Ink5GSlznvvvQd/f3/ExsbiwoUL8PHxwZNPPom3334buHcUZOLEiYiPj4erqytatmyJFStWmG2r8PBwrFq1Cm+88Qbmzp2LyMhIvPPOO3j11Ved6qZfhZGEM5xsKUEpKSnw9vZGcnKykrk+DGRZRkJCAvz8/KBSPXy/kP9NzMA3+y4hOdOyy8GKJuCFLKTABUDJZvjerlq80LQiKvi6leh8SsrDvl7ZU0m0VVZWFuLi4hAeHg4XFxe7TNMZCCFgMBig0WjKxK9sR8rfVh988AEWLlyIy5cvl+g8i1rvLP0O5ZEJIiIiJzB//nw0bNgQ3t7e2L9/P2bOnIkRI0Y4OiyLMJkgKoRBZEJAQIIETb5+E7lGWfmb/SaIyF7Onj2L6dOnIyEhARUrVsQbb7yBiRMnOjosizCZICrE9vTxylNDozw+V8p/PX5NeWpojycec2iMRPTwmD17Nj755JMyeUqIP6uIiIjIJkwmiIiczEPeL56cjD3WN4cmE7GxsWjUqBE8PT0REBCAHj164PTp0yZ12rRpozzdLe/1yiuvOCxmIqKSknd3w4yMDEeHQo+QvPXt/rtrWsOhfSZ27NiBmJgYNGrUCAaDAW+//TY6duyIkydPwt3dXak3bNgwkweZuLmVzcvuiIiKolar4ePjozxrwc3NrUydNy8MLw21XGm2lRACGRkZuHnzJnx8fKBWq4s9LYcmExs3bjR5v3TpUgQEBODQoUNo1aqVUu7m5vbAJ7wRET0M8vZ1zvDwJnsRQkCWZahUKiYTD+CItvLx8bH5O9aprubIu7OYn5+fSfm3336Lb775BkFBQejWrRvefffdQo9OZGdnIzs7W3mfkpIC3LvBjCzLZscpi2RZVla6h5EQMiAEAHucOxb5XtaOhweMe1+5EBCi7K5rD/t6ZU8l2VaBgYEoX758gQdYlVWyLCMlJQVeXl68GdoDlHZbabVaqNVqCCHM9p2wdP12mmRClmWMHj0azZs3R506dZTy559/HpUqVUJISAiOHz+Ot956C6dPn8aaNWvMTic2NhZTp04tUJ6YmAiDwVCiy1CaZFlGamoqhBAP5caZlpoFN5EJAXt8ZgJuyLn3t+WZvpTvfy9kKeWqewmE6t6dNfNzEwakJSchQTYtLyse9vXKnthWlpNlGZmZmdBoNGyrB3C2tkpNTbWontMkEzExMThx4gR27dplUv7yyy8rf9etWxfBwcFo3749zp8/b/JM9zwTJ07E2LFjlfcpKSkIDQ2Fr6/vQ3c7bUmS4Ovr6xQrnL1lqjKQIaUiBfa5nTYAq2+nnf+4xN1x75LvTUOGZFIOAJKkhYe3D/x8yma/nod9vbIntpXl2FaWc7a20mgsSxOcIpkYMWIE1q1bh507d6JChQpF1m3SpAkA4Ny5c2aTCb1eb/ahKCqVyik+GHuSJOmhXC4AkCQVIEl2fJaGlO9lzTh4wHj3lUsSJKlsfyYP83plb2wry7GtLOdMbWVpDA5NJoQQeP311/Hjjz9i+/btCA8Pf+A4ec+tDw4OLoUI6VHW0HUsBAyQ7ttMWlXzhywLqFTsSEZEBEcnEzExMVi+fDl++ukneHp64vr16wAAb29vuLq64vz581i+fDm6du2KcuXK4fjx4xgzZgxatWqFevXqOTJ0egT4qCubLfdz15V6LEREzsyhycSCBQuAezemyi/v+fI6nQ5btmzBnDlzkJ6ejtDQUPTu3Rv/93//56CIiYiI6H4OP81RlNDQUOzYsaPU4iEiIiLrOUUHTCJndMNwGEaRA7WkQ6DmSaX8SmImjEJALUl4zNe1yGkQET0KmEwQFeKvrK+UR5AHevwvmTgQn6A8gvwxXz6CnIjI8dedEBERUZnGZIKIiIhswmSCiIiIbMJkgoiIiGzCZIKIiIhswmSCiIiIbMJkgoiIiGzCZIKIiIhswmSCqBBqSQ8NXKGWTB9pr1FL0KgkaNR8aigREXgHTKLCtXWfZbb86XohpR4LEZEzYzJBZZKrVg2dxsIDa0LATRggSVpAKtmjCV4uGqhLeB5ERM6GyQSVOa5aNVrX9ICQMi2qLwRgTAfU7hKkEv6i16gEMsVtXE1TF1nPQ+sBL71XicZCRFRamExQmaPTqCCkTPx87lckZKY8sL4EoBw8cAdpKPqh97bTqiWEl3eHq67wZMJT54mu4V2ZTBDRQ4PJBJVZCZkpuJWR+MB6EiSoVQK35GQIK9KJROM5yCIXKkkLX3UVpfz61WowGjVQqw0ICjljMo5OrYJvtgG53LSI6BHCqzmICpEh30C6uIYM+YZJeXJyEJISKyA5OchhsRERORMmE0RERGQTJhNERERkEyYTREREZBMmE0RERGQTJhNERERkEyYTREREZBMmE0RERGQTJhNERERkE96mj6gQrqpyMIpcqCWtSbmn5y0YjVqo1bkOi42IyJkwmSAqhJ+6htnykAqnSj0WIiJnxtMcREREZBMmE0RERGQTJhNERERkE/aZICrEdcMBGEUO1JIOQZpGSvn5s01gMOih0WQjoup+h8ZIROQMmEwQFcIocmBENiBMyw0GPQy5Lo4Ki4jI6fA0B5EDSJAcHQIRkd3wyARRKdOr9VBJKlxNu+roUAoQskBqZiqy07IhqQpPeDy0HvDSe5VqbETkvJhMEJUynUqHjNwMbL28Fak5qY4Ox5QAXLJdkKXPQmEHTzx1nuga3pXJBBEpmEwQOUhqTipSclIcHYYpAcgGGWlSWqHJBBHR/dhngoiIiGzCZIKIiIhswmSCiIiIbMJkgoiIiGzi0GQiNjYWjRo1gqenJwICAtCjRw+cPn3apE5WVhZiYmJQrlw5eHh4oHfv3rhx44bDYqZHh486An7qGvBRR5iUBwadQchjfyMw6IzDYiMiciYOTSZ27NiBmJgY7Nu3D7/99htyc3PRsWNHpKenK3XGjBmDX375BStXrsSOHTtw9epV9OrVy5Fh0yPCXRUED1UI3FVBJuU+vtfhW+4KfHyvOyw2IiJn4tBLQzdu3GjyfunSpQgICMChQ4fQqlUrJCcn48svv8Ty5cvRrl07AMCSJUtQs2ZN7Nu3D02bNi0wzezsbGRnZyvvU1LuXnonyzJkWS7xZSotsixDCPFQLVN+QsiAEChwL+u7AyHE3SsXLbmTpHSvVmncdVKCdDdkM2ErBCCEeHA9RxCwLH754V33LPWwb4P2xLaynLO1laVxONV9JpKTkwEAfn5+AIBDhw4hNzcXUVFRSp0aNWqgYsWK2Lt3r9lkIjY2FlOnTi1QnpiYCIPBUKLxlyZZlpGamgohBFSqh6/rS1pqFtxEJgQKfmZuwgBjOlAOHlCrHvxtLEGCn8odgARRwt/eerUK3kZXiIzC6wihRmpSGlTpeqhzXUs0HqsJQCv08AAKvc+Ei3BBanIq9Dn6Ug7OuTzs26A9sa0s52xtlZpq2Y31nCaZkGUZo0ePRvPmzVGnTh0AwPXr16HT6eDj42NSNzAwENevmz/EPHHiRIwdO1Z5n5KSgtDQUPj6+sLL6+G5Y58sy5AkCb6+vk6xwtlbpioDGVIqUpBbYJgkaaF2l3AHabglJz9wWnePSAhcl1OsSiZyRToEBCRI0EruSnl2ltu9UgG9i2nW4CGp4SLl4kpSJnIM5udV1U+LSHcXnEm7g2tpiRbHUxr0Ggm1vQJh1GYWmkyodCp4envCz8OvtMNzKg/7NmhPbCvLOVtbaTSWpQlWJxM7d+5Es2bNCszAYDBgz549aNWqlbWTBADExMTgxIkT2LVrV7HGz6PX66HXF/zFpFKpnOKDsSdJkh7K5QIASVIBkmT+57EkQZKke0fjLUsO8upak0zcMByBEdlQQ4/HtM2V8ri4BjDkukCjzUL1mjsLzAcSkG2QkWM0f3gwV75bKUeWkW00WhxP6VDnnT8q/A6YEiCppIdyvbPWw7wN2hvbynLO1FaWxmB1pG3btkVCQkKB8uTkZLRt29bayQEARowYgXXr1mHbtm2oUKGCUh4UFIScnBwkJSWZ1L9x4waCgoLMTImIiIhKm9XJhBACklTwJ8udO3fg7u5udpyipjVixAj8+OOP+P333xEeHm4yvEGDBtBqtdi6datSdvr0aVy6dAmRkZHWhk5EREQlwOLTHHmXY0qShOjoaJNTCUajEcePH0ezZs2smnlMTAyWL1+On376CZ6enko/CG9vb7i6usLb2xtDhgzB2LFj4efnBy8vL7z++uuIjIw02/mSiIiISp/FyYS3tzdw72iCp6cnXF3/1wtdp9OhadOmGDZsmFUzX7BgAQCgTZs2JuVLlixBdHQ0AGD27NlQqVTo3bs3srOz0alTJ8yfP9+q+RAREVHJsTiZWLJkCQAgLCwM48aNs/qUhjlCPLgznIuLC+bNm4d58+bZPD8iIiKyP6uv5pg8eXLJREJERERlktUdMG/cuIEXX3wRISEh0Gg0UKvVJi8iIiJ6tFh9ZCI6OhqXLl3Cu+++i+DgYLNXdhAREdGjw+pkYteuXfjjjz9Qv379komIiIiIyhSrk4nQ0FCLOk4SlXVBmobK7bTzq1xlHyAkQOJ2QESE4vSZmDNnDiZMmID4+PiSiYjISaglPTSSC9SS6e3ZtdocaHXZ0GpzHBYbEZEzsfrIxHPPPYeMjAxERETAzc0NWq3WZLi5W20TERHRw8vqZGLOnDklEwkRERGVSVYnE4MGDSqZSIicTJp8BbIwQiWp4aF6TClPuPMYZFkDlcoAv3JXHBojEZEzsDqZuHTpUpHDK1asaEs8RE4j2RivPII8fzJx62aE8ghyJhNERMVIJsLCwoq8t4TRaLQ1JiIiIipDrE4mjhw5YvI+NzcXR44cwSeffIL333/fnrERERFRGWB1MvH4448XKGvYsCFCQkIwc+ZM5VHlRERE9Giw+j4ThalevToOHDhgr8kRERFRGWH1kYmUlBST90IIXLt2DVOmTEHVqlXtGRsRERGVAVYnEz4+PgU6YAohEBoaihUrVtgzNiIiIioDrE4mtm3bZvJepVLB398fVapUgUZj9eSIiIiojLP6279169YlEwkRERGVScU6lHD+/HnMmTMHp06dAgDUqlULo0aNQkREhL3jI3IYjeQKldBAJZk+f0avy4BaZYBGwwd9ERGhOMnEpk2b0L17d9SvXx/NmzcHAOzevRu1a9fGL7/8gg4dOpREnESlLlDzpNnysIiDpR4LEZEzszqZmDBhAsaMGYMZM2YUKH/rrbeYTBARET1irL7PxKlTpzBkyJAC5S+99BJOnjxpr7iIiIiojLA6mfD398fRo0cLlB89ehQBAQH2iouIiIjKCKtPcwwbNgwvv/wyLly4gGbNmgH3+kx8+OGHGDt2bEnESOQQtw1/Q0YuVNCivKa2Uv7vpbowGLTQaHJRoeJfDo2RiMgZWJ1MvPvuu/D09MSsWbMwceJEAEBISAimTJmCkSNHlkSMRA6RLZKUR5Dnl57uqzyCnIiIipFMSJKEMWPGYMyYMUhNTQUAeHp6lkRsREREVAZYnUzExcXBYDCgatWqJknE2bNnodVqERYWZu8YiYiIyIlZ3QEzOjoae/bsKVC+f/9+REdH2ysuIiIiKiOsTiaOHDmi3Kwqv6ZNm5q9yoOIiIgeblYnE5IkKX0l8ktOTobRaLRXXERERFRGWJ1MtGrVCrGxsSaJg9FoRGxsLFq0aGHv+IiIiMjJWd0B88MPP0SrVq1QvXp1tGzZEgDwxx9/ICUlBb///ntJxEhEREROzOojE7Vq1cLx48fRt29f3Lx5E6mpqRg4cCD++ecf1KlTp2SiJCIiIqdVrEeQh4SE4IMPPrB/NEROxEMVAhkGqO7bTHz9/oVs1EClNjgsNiIiZ1KsZILoUeCtDjdbHhB4odRjISJyZlaf5iAiIiLKj8kEERER2YTJBBEREdmkWH0mDAYDtm/fjvPnz+P555+Hp6cnrl69Ci8vL3h4eNg/SiIHuJK7W3lq6GPa/9319fSpVspTQ6vX3OnQGImInIHVRyYuXryIunXr4plnnkFMTAxu3boF3Lv/xLhx46ya1s6dO9GtWzeEhIRAkiSsXbvWZHh0dDQkSTJ5de7c2dqQiYiIqARZnUyMGjUKDRs2RGJiIlxdXZXynj17YuvWrVZNKz09HY8//jjmzZtXaJ3OnTvj2rVryuu7776zNmQiIiIqQVaf5vjjjz+wZ88e6HQ6k/KwsDBcuXLFqml16dIFXbp0KbKOXq9HUFCQtWESERFRKbE6mZBl2ewDvf799194enraKy7F9u3bERAQAF9fX7Rr1w7Tp09HuXLlCq2fnZ2N7Oxs5X1KSooStyzLdo/PUWRZhhDioVqm/FQQ8NKrASEKDPPSqyFBhgRAgvTAaUn3allSt6hpWFIuAYDIm2ch8xOAEAKSsC2mkpAXPwo2+/8IQMgP77pnqYd9G7QntpXlnK2tLI3D6mSiY8eOmDNnDhYtWgTce4poWloaJk+ejK5du1ofaRE6d+6MXr16ITw8HOfPn8fbb7+NLl26YO/evVCr1WbHiY2NxdSpUwuUJyYmwmB4eO5YKMsyUlNTIYSASmX9RTnpOenINGaWSGz2kGOQUdM/B0YzK7JWrYabDPhLHlCrivrWu0uCBD+VOwAJoshvSVNXIcEIQAUJQSpvpfwsVDAAUENlUg4ArpIa3kY9AlVa5ArzG6G37IKslAz4wQOyyjl2GHlc1Cp4yK5Iy8jLLAqScnW4cycZaQmAbCbZczRXnRoeem2Jz8fWbfBRwraynLO1lbmnhJtjdTIxa9YsdOrUCbVq1UJWVhaef/55nD17FuXLl7d7f4Z+/fopf9etWxf16tVDREQEtm/fjvbt25sdZ+LEiRg7dqzyPiUlBaGhofD19YWXl5dd43MkWZYhSRJ8fX2LtcJlp2VjV/wupOZYtqKUtswcI+JvpyNXLvhlFe7zGHrWbI0kKR3X5eQHTuvur3+B63KKVcmEfK+uDGEyHyNk5f/75+8h1NCpXXFDzkBOIRm9l8oHLl5uSECaRfGXJk9JjQC44FTKTWQbzLdVeVdvNBAa/BWXgZQs50rQvVw06FY/GH4+biU+L1u3wUcJ28pyztZWGo1laYLVyUSFChVw7NgxrFixAsePH0daWhqGDBmCAQMGmHTILAmVK1dG+fLlce7cuUKTCb1eD71eX6BcpVI5xQdjT5IkFXu5JJWE1NxUpOSmlEhstkrPMeBWZhpyjAW/kH1dvSFJEoQEi5ODu0fuhVXJhOn45se7v1zg7i/6Iuclwer4S0te/NkGgWwzpzMBIEcWMAqBlGwZyU6WTECSIEmlt63bsg0+athWlnOmtrI0hmLdZ0Kj0eCFF14ozqg2+ffff3Hnzh0EBweX+ryJiIjIPIuSiZ9//tniCXbv3t3iumlpaTh37pzyPi4uDkePHoWfnx/8/PwwdepU9O7dG0FBQTh//jzefPNNVKlSBZ06dbJ4HkRERFSyLEomevToYdHEJEkye6VHYQ4ePIi2bdsq7/P6OgwaNAgLFizA8ePHsWzZMiQlJSEkJAQdO3bEe++9Z/Y0BpG9lVPXgoAM6b7bsVQI/QtCqCBJztV5kojIUSxKJkrqEpU2bdpAFNEbfNOmTSUyXyJLuKh8zZa7eySWeixERM7M8b07iIiIqEwrVjKxdetWPP3004iIiEBERASefvppbNmyxf7RERERkdOzOpmYP38+OnfuDE9PT4waNQqjRo2Cl5cXunbtWuQzNojKmiw5EZnyHWTJpqc10tN8kZZaDulp5k+DEBE9aqy+NPSDDz7A7NmzMWLECKVs5MiRaN68OT744APExMTYO0Yih7hjPPm/R5Cr/vcI8n8v1+UjyImI8rH6yERSUpLZx4B37NgRycnOdTc/IiIiKnlWJxPdu3fHjz/+WKD8p59+wtNPP22vuIiIiKiMsPo0R61atfD+++9j+/btiIyMBADs27cPu3fvxhtvvIG5c+cqdUeOHGnfaImIiMjpWJ1MfPnll/D19cXJkydx8uRJpdzHxwdffvml8l6SJCYTREREjwCrk4m4uLiSiYSIiIjKJN60ioiIiGxi9ZEJIQRWrVqFbdu24ebNmwVutb1mzRp7xkdEREROzupkYvTo0fjiiy/Qtm1bBAYGQpKkkomMiIiIygSrk4mvv/4aa9asQdeuXUsmIiIiG/D3DVHpszqZ8Pb2RuXKlUsmmjIkJTsFablpDpu/kAVSM1ORnZYNSWXd3lMtqWGQDSUW28PiMW1zs+W866XzctGqoJIk/JuYUWgdSQLU6mxkGGzbfoUQSE9NQ7Iqw65HaNUqCVp10d3ZPLQe8NJ72W2eRLayOpmYMmUKpk6diq+++gqurq4lE1UZkJabhvVx65Gak+qYAATgku2CLH0WYOV+LMg9CE2DmpZUZEQOo1OrkJ5twOa/byAly3zC7OWiQb0wgZ/PrUNiVkqx5yUJwA8eSEAahJ1yCa1ahQq+rnDVqQut46nzRNfwrkwmyKlYnUz07dsX3333HQICAhAWFgatVmsy/PDhw/aMz6ml5qQiJaf4OyObCEA2yEiT0qxOJjy1niUVFZFTSMkyIDkzt9DhBlnCzfRk3MpILLTOg0iQIKtkXJeTISCKPZ38dGoVPNw8kGv9rpnIoaxeYwcNGoRDhw7hhRdeYAdMIiIisj6Z+PXXX7Fp0ya0aNGiZCIichLJxjjIMEAFDbzV4Ur5zRuVIRs1UKkNCAi84NAYiYicgdXJRGhoKLy8eK6OHn5p8lXlEeT5k4nEhArKI8iZTBARFeMOmLNmzcKbb76J+Pj4komIiIiIyhSrj0y88MILyMjIQEREBNzc3Ap0wExISLBnfEREROTkrE4m5syZUzKREBERUZlUrKs5iIiIiPLYdDFzVlYWcnJyTMrYOdM55Bhk5Bpls8OydEYYhUBWjoz0bOe7E6YECbJ9LtsnIqJSYHUykZ6ejrfeegs//PAD7ty5U2C40Wi0V2xkg1yjjIt3MpBjKJhQSIYMZBuMuJyYgWupjrsleGHcdGoEebs4OgwiIrKQ1cnEm2++iW3btmHBggV48cUXMW/ePFy5cgVffPEFZsyYUTJRUrHkGGTkmDk6kWsEhLibcJgb7mg62eqLjIiIyIGsTiZ++eUX/Pe//0WbNm0wePBgtGzZElWqVEGlSpXw7bffYsCAASUTKRERETklq38CJiQkKE8N9fLyUi4FbdGiBXbu5NMU6eGhl3zgIvlBL/mYlLu7J8Ld4zbc3Yv/XAciooeJ1UcmKleujLi4OFSsWBE1atTADz/8gMaNG+OXX36Bj4+PBVMgKhvKa2qbLa9Q8a9Sj4WIyJlZfWRi8ODBOHbsGABgwoQJmDdvHlxcXDBmzBiMHz++JGIkIiIiJ2b1kYkxY8Yof0dFReHUqVM4fPgwqlSpgnr16tk7PiIiInJyNt1nAgDCwsIQFhZmn2iIiIiozLH4NMfevXuxbt06k7L//ve/CA8PR0BAAF5++WVkZ2eXRIxEDnHDcBjXcvfjhuGwSXn8+YY4d7oZ4s83dFhsRETOxOJkYtq0afj777+V93/99ReGDBmCqKgoTJgwAb/88gtiY2NLKk6iUmcQmchFOgwi06Q8O8cN2dkeyM5xc1hsRETOxOJk4ujRo2jfvr3yfsWKFWjSpAkWL16MsWPHYu7cufjhhx9KKk4iIiJyUhYnE4mJiQgMDFTe79ixA126dFHeN2rUCJcvX7Z/hEREROTULE4mAgMDERcXBwDIycnB4cOH0bRpU2V4amoqtFptyURJRERETsviZKJr166YMGEC/vjjD0ycOBFubm5o2bKlMvz48eOIiIgoqTiJiIjISVl8aeh7772HXr16oXXr1vDw8MCyZcug0+mU4V999RU6duxYUnESERGRk7I4mShfvjx27tyJ5ORkeHh4QK1WmwxfuXIlPDw8SiJGIiIicmJW307b29u7QCIBAH5+fiZHKiyxc+dOdOvWDSEhIZAkCWvXrjUZLoTApEmTEBwcDFdXV0RFReHs2bPWhkxEREQlyOpkwp7S09Px+OOPY968eWaHf/TRR5g7dy4WLlyI/fv3w93dHZ06dUJWVlapx0pERETm2Xw7bVt06dLF5PLS/IQQmDNnDv7v//4PzzzzDHDvjpuBgYFYu3Yt+vXrV8rR0qPGWx0GWRihkkyPxPkHnIcsa6BSGRwWGxGRM3FoMlGUuLg4XL9+HVFRUUqZt7c3mjRpgr179xaaTGRnZ5vc1jslJQUAIMsyZFm2W3xCFoDA3ZcjCBQ9fwFI9/6ZGyaEgHSvjrORYN/486Zi7bJ6qiqYLS9X7ur90VoeezHiL03/i7/w2CQAakmCl14FCEdtAOZ56NRQSwJeenWhsXnp1ZAgF2udyK+469WDpvnA/Yq4u/+x5/6spMmyDCHKVsyO4mxtZWkcTptMXL9+Hbh3f4v8AgMDlWHmxMbGYurUqQXKExMTYTDY75dkamYqXLJdIBsc9IELwNXgevdvM/sytdGIQJUGuaJgfN6yC7JSMuAHD8gq51hh83OV1PA26hGo0tolfgkS/FTuACSIEs7+HhQ7nLz93SQ13GU3BKi8C40/SF0OnsKImv7ZMDrJDi+PTp2DlJQM1PTPLTQ2rVoNNxnwlzygVhV/fSiJ9UqrUsHL6ArX3IL90vK4CBekJqdCn6O3yzxLgyzLSE1NhRACKpVDz647PWdrq9TUVIvqOW0yUVwTJ07E2LFjlfcpKSkIDQ2Fr68vvLy87Daf7LRsZOmzkCal2W2aVrm370rTpplNJtJlA27I6cgxs0P1UvnAxcsNCUjDdTm5FIK1jodQQ6d2xQ05wy7x3/3lKHBdTinxZOJBscPJ299TqFFepcdNORnZhcWv9YHaQ40N57fhVoZzxe+mVSPI2wWXEzKQK5v/rMN9HkPPmq2RJKXb1P4lsV7pJBU81QYYtYXvmlU6FTy9PeHn4WeXeZYGWZYhSRJ8fX2d4gvSmTlbW2k0lqUJTptMBAUFAQBu3LiB4OBgpfzGjRuoX79+oePp9Xro9QUzdpVKZdcPRlJJeceCHUdC4TFIgLj3z9wwSZIg7tVxNgL2j//ukeNCplcIo8iGgIAECWrpf+tUbq4OEBIgCWi1OdbFXsz4S8v/4i8itnvxJ2an4GZGYilHWDQPvQYurq64lZmOHKP5ZMjX1dtu7V+c9aro6YkH71eku/sfZ/iisYYkSXbfDz+snKmtLI3B8ZEWIjw8HEFBQdi6datSlpKSgv379yMyMtKhsdGj4brhIK4a9uC64aBJ+YVzTXHmn9a4cK5poeMSET1KHHpkIi0tDefOnVPex8XF4ejRo/Dz80PFihUxevRoTJ8+HVWrVkV4eDjeffddhISEoEePHo4Mm4iIiPJxaDJx8OBBtG3bVnmf19dh0KBBWLp0Kd58802kp6fj5ZdfRlJSElq0aIGNGzfCxcXFgVETERFRfg5NJtq0aQNRxKVlkiRh2rRpmDZtWqnGRURERJZz2j4TREREVDYwmSAiIiKbMJkgIiIimzCZICIiIpswmSAiIiKbMJkgIiIimzjt7bSJHC1AU1+5nXZ+YeEH75U6162wiYgchckEUSG0krvZcr1LRqnHQkTkzHiag4iIiGzCZIKIiIhswtMcRIVIl69DQIYEFdxVQUp5UmIQhKyGpDLCx/e6Q2MkInIGTCaICpFkPA8jsqGG3iSZuHG9Ggy5LtBos5hMEBHxNAcRERHZiskEERER2YTJBBEREdmEyQQRERHZhMkEERER2YTJBBEREdmEyQQRERHZhMkEERER2YQ3rSIqhFrSAeLe//loNNkm/xMRPeqYTBAVIkjTyGx5RNX9pR4LEZEz42kOIiIisgmTCSIiIrIJkwkiIiKyCftMEBUiwfgPjCIXakkLP3UNpfzqvzVhNGqhVucipMIph8ZIROQMmEwQFSJTvnP3EeRCD6j/V56a6q88ghxgMkFExNMcREREZBMmE0RERGQTJhNERERkEyYTREREZBMmE0RERGQTJhNERERkE14aSkTkJNQqCZIkIT3bUGgdLYzIyjXiWlImjEKUanyW8HTRwNtVZ0FNepgwmSAichJqlQSDUca/iZnIMchm65R30+BqUhaOx19BSlbhSYcjeLlo0L1+CJOJRxCTCaJCuKkCIYtcqCStSbm393UYjRqo1c61I6eHR45BRo7RfDKRaxQwyDJSsgSSM3NLPTYic5hMEBXCV13FbHlQyJlSj4WIyJmxAyYRERHZhMkEERER2cSpk4kpU6ZAkiSTV40aNSwYk4iIiEqL0/eZqF27NrZs2aK812icPmR6SFzN3Xf3qaHQI0TbVCk/e7o5DLl6aLTZqFp9t0NjJCJyBk7/zazRaBAUFOToMOgRJGBUXvnJshqyrIEs82oOIiKUhWTi7NmzCAkJgYuLCyIjIxEbG4uKFSsWWj87OxvZ2dnK+5SUFACALMuQZfOXWhWHkAUgcPflCAJFz18A0r1/5oYJISDdq+NsJNg3/ryp2LKshY17f/kDY4dzt///4i8itjIRf8m3vz3Wq4LTfHD8Eu7GD+HIHVAhhIAQBfe1sixDCGHXffDDytnaytI4nDqZaNKkCZYuXYrq1avj2rVrmDp1Klq2bIkTJ07A09PT7DixsbGYOnVqgfLExEQYDPb7JZmamQqXbBfIhdxYpsQJwNXgevdvM/sctdGIQJUGuaJgfN6yC7JSMuAHD8gq51hh83OV1PA26hGo0tolfgkS/FTuACQIK3a+VyHBCEAFCUEqb6X8LFQwAFBDZVJuSezFib80uUlquMtuCFB5l8n4S7P9i7teFcWS+H3gAWN6BtwEIOBcR8fchAFpyUlIkLNMymVZRmpqKoQQUKmcuquewzlbW6WmplpUz6mTiS5duih/16tXD02aNEGlSpXwww8/YMiQIWbHmThxIsaOHau8T0lJQWhoKHx9feHl5WW32LLTspGlz0KalGa3aVrl3r4rTZtmNplIlw24Iacjx0xW6aXygYuXGxKQhutycikEax0PoYZO7YobcoZd4r/7C0/gupxi1U5fvldXhjCZjxGy8v/9839Q7MWJvzR5CjXKq/S4KScjuwzGX5rtX9z1qiiWxG+EBLW7GzIkgRQ4102rJEkLD28f+Pm4mZTLsgxJkuDr6+sUX5DOzNnaytJ+ik6dTNzPx8cH1apVw7lz5wqto9frodfrC5SrVCq7fjCSSso7Fuw4EgqPQQLEvX/mhkmSBHGvjrMRsH/8dw8IFzI9i8Y3P9795Q+MHc7d/v+Lv4jYykT8pdP+tq5X5qb3oPgF7sZ/d7t3rtNMkCRIkvl9rSRJdt8PP6ycqa0sjcHxkVohLS0N58+fR3BwsKNDISIionucOpkYN24cduzYgfj4eOzZswc9e/aEWq1G//79HR0aERER3ePUpzn+/fdf9O/fH3fu3IG/vz9atGiBffv2wd/f39GhERER0T1OnUysWLHC0SEQERHRAzh1MkHkSH7q6hCQId13NjDksZOQZTVUKmOh4xIRPUqYTBAVwlVV3my5p9ftUo+FiMiZMZkgIipjNCoVvFyc7LJQAF4uGmhUzhcXlTwmE0REZYi71hWBXq6oF5YJg+xcX9x6jUCmfBsXEk3LhRBIT01DsioDkiTBTeMBo1F/947gTsbTRQNvV52jwyhzmEwQFSJHpNx9hoMkQSf97+6pmRmeEEIFSZLh6mbZrWaJ7EWv0SHbmImfz63DzXTnugOpm06NYG8XXEvOQq7xf3fwlATgBw8kIA0+rl7oXuVpHI+XkJLlXLcD93LRoHv9ECYTxcBkgqgQtwx/KY8gf0zbXCm/dPEJGHJdoNFmoXrNnQ6NkR5diVkpuJWRaEHN0uNh1MDVxRXXUtORkz+ZgARZdff28zmygEGWkZIlkJzpXLcDp+Jz6ptWERERkfNjMkFEREQ2YTJBRERENmEyQURERDZhMkFEREQ2YTJBRERENmEyQURERDZhMkFEREQ2YTJBRERENuEdMIkKEaxpYra8SrXdpR4LEZEzYzJBVAiVZH7zUKuNpR4LEZEz42kOIiIisgmTCSIiIrIJT3MQFSLFeAkyDFBBAy91RaX89q1KkI0aqNQGlPe/6NAYiYicAZMJokKkypeVR5DnTybu3K6kPIKcyQQREU9zEBERkY2YTBAREZFNmEwQERGRTZhMEBERkU2YTBAREZFNmEwQERGRTZhMEBERkU2YTBAREZFNeNMqokLoJE8YoYcaOpNyV5cUGLRZ0KhzHBYbEZEzYTJBVAh/TT2z5RXDj5Z6LEREzoynOYiIiMgmTCaIiIjIJjzNQUREpUqjUsHLRXJ0GAV4uWigUTlfXGUBkwmiQtwyHIcROVBDZ9J/4lJcfRiMOmjUOew/QWQld60rAr1cUS8sEwbZub649RqBTPk2LiQWXc9N4wGjUQ8h7B+DEDLSUrOQqcqAJFl/8sDTRQNvV50FNe2LyQRRIXJEqvII8vwys7yUR5ATkXX0Gh2yjZn4+dw63ExPdnQ4Jtx0agR7u+BachZyjbLZOr4uXuhe5Wkcj5eQkmWwfxBCwE1kIkNKBSTrki0vFw261w9hMkFERI+GxKwU3Mp4wCGAUuZh1MDVxRXXUtORU0gykWsUMMgyUrIEkjNzSyAKAQEDUpALwLmO3BSFHTCJiIjIJkwmiIiIyCZMJoiIiMgmZSKZmDdvHsLCwuDi4oImTZrgzz//dHRIREREdI/TJxPff/89xo4di8mTJ+Pw4cN4/PHH0alTJ9y8edPRoREREVFZSCY++eQTDBs2DIMHD0atWrWwcOFCuLm54auvvnJ0aEREROTsl4bm5OTg0KFDmDhxolKmUqkQFRWFvXv3mh0nOzsb2dnZyvvk5LvXMSclJUGWzV/qUxwpaSnQZGugyyn963kBAAKQciToZJ3Zq4eMuWqUV8nIFQWX2UPWIC01Fd7CBbkqt9KJ1wouUMMlV4/yKmGX+CUBeAodjJIrhBVXWl3JEjAKI9SSgL/+f/M5nZMFOUcGRA7875v/g2IvTvylyVWooclWobzkihxV2Yu/NNu/uOtVUcr6+lNY/PnbqizGn5+30CMrLR1agwS9XDL3mVCLTOglWH2fCa3BgJTkJCRJ9nuicUpKyr2wHnCHLuHErly5IgCIPXv2mJSPHz9eNG7c2Ow4kydPFne/avniiy+++OKLL3u8Ll++XOT3tVMfmSiOiRMnYuzYscp7WZaRkJCAcuXKQbIyy3NmKSkpCA0NxeXLl+Hl5eXocJwa28pybCvLsa0sx7aynLO1lRACqampCAkJKbKeUycT5cuXh1qtxo0bN0zKb9y4gaCgILPj6PV66PWmtz/28fEp0TgdycvLyylWuLKAbWU5tpXl2FaWY1tZzpnaytvb+4F1nLoDpk6nQ4MGDbB161alTJZlbN26FZGRkQ6NjYiIiO5y6iMTADB27FgMGjQIDRs2ROPGjTFnzhykp6dj8ODBjg6NiIiIykIy8dxzz+HWrVuYNGkSrl+/jvr162Pjxo0IDAx0dGgOpdfrMXny5AKndKggtpXl2FaWY1tZjm1lubLaVpJ44PUeRERERIVz6j4TRERE5PyYTBAREZFNmEwQERGRTZhMEBERkU2YTDih2NhYNGrUCJ6enggICECPHj1w+vTpIsdZvHgxWrZsCV9fX/j6+iIqKuqReFR7cdoqvxUrVkCSJPTo0aNE43QGxW2rpKQkxMTEIDg4GHq9HtWqVcP69etLJWZHKW5bzZkzB9WrV4erqytCQ0MxZswYZGVllUrMjrJgwQLUq1dPuclSZGQkNmzYUOQ4K1euRI0aNeDi4oK6des+9OtTHmvbqizt15lMOKEdO3YgJiYG+/btw2+//Ybc3Fx07NgR6enphY6zfft29O/fH9u2bcPevXsRGhqKjh074sqVK6Uae2krTlvliY+Px7hx49CyZctSidXRitNWOTk56NChA+Lj47Fq1SqcPn0aixcvxmOPPVaqsZe24rTV8uXLMWHCBEyePBmnTp3Cl19+ie+//x5vv/12qcZe2ipUqIAZM2bg0KFDOHjwINq1a4dnnnkGf//9t9n6e/bsQf/+/TFkyBAcOXIEPXr0QI8ePXDixIlSj720WdtWZWq/bs8Hc1HJuHnzpgAgduzYYfE4BoNBeHp6imXLlpVobM7G0rYyGAyiWbNm4j//+Y8YNGiQeOaZZ0otRmdhSVstWLBAVK5cWeTk5JRqbM7GkraKiYkR7dq1MykbO3asaN68eSlE6Fx8fX3Ff/7zH7PD+vbtK5566imTsiZNmojhw4eXUnTOpai2up8z79d5ZKIMyHuMup+fn8XjZGRkIDc316pxHgaWttW0adMQEBCAIUOGlFJkzseStvr5558RGRmJmJgYBAYGok6dOvjggw9gNBpLMVLHs6StmjVrhkOHDimHoS9cuID169eja9eupRanoxmNRqxYsQLp6emFPvJg7969iIqKMinr1KkT9u7dW0pROgdL2up+Tr1fd3Q2Q0UzGo3iqaeesvrXzauvvioqV64sMjMzSyw2Z2NpW/3xxx/iscceE7du3RJCiEfyyISlbVW9enWh1+vFSy+9JA4ePChWrFgh/Pz8xJQpU0otVkezZhv89NNPhVarFRqNRgAQr7zySqnE6GjHjx8X7u7uQq1WC29vb/Hrr78WWler1Yrly5eblM2bN08EBASUQqSOZ01b3c+Z9+tMJpzcK6+8IipVqvTAZ8nnFxsbK3x9fcWxY8dKNDZnY0lbpaSkiLCwMLF+/Xql7FFMJixdr6pWrSpCQ0OFwWBQymbNmiWCgoJKIUrnYGlbbdu2TQQGBorFixeL48ePizVr1ojQ0FAxbdq0UovVUbKzs8XZs2fFwYMHxYQJE0T58uXF33//bbbuo55MWNNW+Tn7fp3JhBOLiYkRFSpUEBcuXLB4nJkzZwpvb29x4MCBEo3N2VjaVkeOHBEAhFqtVl6SJAlJkoRarRbnzp0rtZgdxZr1qlWrVqJ9+/YmZevXrxcARHZ2dglG6RysaasWLVqIcePGmZR9/fXXwtXVVRiNxhKM0vm0b99evPzyy2aHhYaGitmzZ5uUTZo0SdSrV6+UonMuRbVVnrKwX2efCSckhMCIESPw448/4vfff0d4eLhF43300Ud47733sHHjRjRs2LDE43QG1rZVjRo18Ndff+Ho0aPKq3v37mjbti2OHj2K0NDQUou9tBVnvWrevDnOnTsHWZaVsjNnziA4OBg6na6EI3ac4rRVRkYGVCrTXaparVam9yiRZRnZ2dlmh0VGRmLr1q0mZb/99pvF/QYeNkW1FcrSft3R2QwV9Oqrrwpvb2+xfft2ce3aNeWVkZGh1HnxxRfFhAkTlPczZswQOp1OrFq1ymSc1NRUBy1F6ShOW93vUTnNUZy2unTpkvD09BQjRowQp0+fFuvWrRMBAQFi+vTpDlqK0lGctpo8ebLw9PQU3333nbhw4YLYvHmziIiIEH379nXQUpSOCRMmiB07doi4uDhx/PhxMWHCBCFJkti8ebMQZtpp9+7dQqPRiI8//licOnVKTJ48WWi1WvHXX385cClKh7VtVZb260wmnBAAs68lS5YodVq3bi0GDRqkvK9UqZLZcSZPnuygpSgdxWmr+z0qyURx22rPnj2iSZMmQq/Xi8qVK4v333/fpA/Fw6g4bZWbmyumTJkiIiIihIuLiwgNDRWvvfaaSExMdNBSlI6XXnpJVKpUSeh0OuHv7y/at2+vfDmKQtapH374QVSrVk3odDpRu3ZtqzohlmXWtlVZ2q/zEeRERERkE/aZICIiIpswmSAiIiKbMJkgIiIimzCZICIiIpswmSAiIiKbMJkgIiIimzCZICIiIpswmSAiIiKbMJkgohITFhaGOXPmODoMAMDSpUvh4+Nj1Tjbt2+HJElISkoqsbge5N1338XLL79cZJ02bdpg9OjRVk335MmTqFChAtLT022MkIjJBBGuX7+O119/HZUrV4Zer0doaCi6detW4GFEtirODp+Kx5mSGFtcv34dn376Kd555x2rxmvTpg0kSVJegYGB6NOnDy5evKjUqVWrFpo2bYpPPvmkBCKnRw2TCXqkxcfHo0GDBvj9998xc+ZM/PXXX9i4cSPatm2LmJgYR4dnVk5OjqNDoFLyn//8B82aNUOlSpWsHnfYsGG4du0arl69ip9++gmXL1/GCy+8YFJn8ODBWLBgAQwGgx2jpkcRkwl6pL322muQJAl//vknevfujWrVqqF27doYO3Ys9u3bp9RLSkrC0KFD4e/vDy8vL7Rr1w7Hjh1Thk+ZMgX169fH119/jbCwMHh7e6Nfv35ITU0FAERHR2PHjh349NNPlV+L8fHxAIATJ06gS5cu8PDwQGBgIF588UXcvn1bmXabNm0wYsQIjB49GuXLl0enTp0KLMeJEyegUqlw69YtAEBCQgJUKhX69eun1Jk+fTpatGgBADAajRgyZAjCw8Ph6uqK6tWr49NPP1Xqbt68GS4uLgUO748aNQrt2rVT3u/atQstW7aEq6srQkNDMXLkyCIPm9vajgCQmpqKAQMGwN3dHcHBwZg9e7bJUZ82bdrg4sWLGDNmjNLW+W3atAk1a9aEh4cHOnfujGvXrhUarzmrV69G7dq1odfrERYWhlmzZpkMnz9/PqpWrQoXFxcEBgbi2WefVYatWrUKdevWhaurK8qVK4eoqKgi22vFihXo1q2bSVl6ejoGDhwIDw8PBAcHF5h/Hjc3NwQFBSE4OBhNmzbFiBEjcPjwYZM6HTp0QEJCAnbs2GFVGxDdj8kEPbISEhKwceNGxMTEwN3dvcDw/OfX+/Tpg5s3b2LDhg04dOgQnnzySbRv3x4JCQlKnfPnz2Pt2rVYt24d1q1bhx07dmDGjBkAgE8//RSRkZHKr8Vr164hNDQUSUlJaNeuHZ544gkcPHgQGzduxI0bN9C3b1+TWJYtWwadTofdu3dj4cKFBWKtXbs2ypUrp3wp/PHHHybvAWDHjh1o06YNAECWZVSoUAErV67EyZMnMWnSJLz99tv44YcfAADt27eHj48PVq9erYxvNBrx/fffY8CAAcrydu7cGb1798bx48fx/fffY9euXRgxYkShbW5rOwLA2LFjsXv3bvz888/47bff8Mcff5h8Sa5ZswYVKlTAtGnTlLbOk5GRgY8//hhff/01du7ciUuXLmHcuHGFxnu/Q4cOoW/fvujXrx/++usvTJkyBe+++y6WLl0KADh48CBGjhyJadOm4fTp09i4cSNatWoFALh27Rr69++Pl156CadOncL27dvRq1cvFPasxYSEBJw8eRINGzY0KR8/fjx27NiBn376CZs3b8b27dsLJAnmpvXDDz+gSZMmJuU6nQ7169fHH3/8YXEbEJnl6MeWEjnK/v37BQCxZs2aIuv98ccfwsvLS2RlZZmUR0REiC+++EIIIcTkyZOFm5ubSElJUYaPHz9eNGnSRHnfunVrMWrUKJNpvPfee6Jjx44mZZcvXxYAxOnTp5XxnnjiiQcuT69evURMTIwQQojRo0eL8ePHC19fX3Hq1CmRk5Mj3NzcTB53fL+YmBjRu3dv5f2oUaNEu3btlPebNm0Ser1eeaT2kCFDxMsvv1ygrVQqlcjMzBTi3iOUZ8+ebbd2TElJEVqtVqxcuVIZnpSUJNzc3EzaNv988yxZskQAEOfOnVPK5s2bJwIDAwttk23btgkAyjI///zzokOHDiZ1xo8fL2rVqiWEEGL16tXCy8vLJP48hw4dEgBEfHx8ofPL78iRIwKAuHTpklKWmpoqdDqd+OGHH5SyO3fuCFdXV5Plb926tdBqtcLd3V24ubkJAKJatWoiLi6uwHx69uwpoqOjLYqJqDA8MkGPrMJ+Ed7v2LFjSEtLQ7ly5eDh4aG84uLicP78eaVeWFgYPD09lffBwcG4efPmA6e9bds2k+nWqFEDuPcLPU+DBg0eGGfr1q2xfft24N5RiHbt2qFVq1bYvn07Dhw4gNzcXDRv3lypP2/ePDRo0AD+/v7w8PDAokWLcOnSJWX4gAEDsH37dly9ehUA8O233+Kpp55SjtgcO3YMS5cuNYm9U6dOkGUZcXFxJdKOFy5cQG5uLho3bqwM9/b2RvXq1R/YPrh36D8iIsLstC1x6tQpkzYEgObNm+Ps2bMwGo3o0KEDKlWqhMqVK+PFF1/Et99+i4yMDADA448/jvbt26Nu3bro06cPFi9ejMTExELnlZmZCQBwcXFRys6fP4+cnByTIwx+fn5ml3/AgAE4evQojh07hl27dqFKlSro2LGjySkjAHB1dVViJCoujaMDIHKUqlWrQpIk/PPPP0XWS0tLQ3BwsPJFnV/+UyFardZkmCRJkGX5gdPu1q0bPvzwwwLDgoODlb/NnYa5X16/gbNnz+LkyZNo0aIF/vnnH2zfvh2JiYlo2LAh3NzcgHvn4seNG4dZs2YhMjISnp6emDlzJvbv369Mr1GjRoiIiMCKFSvw6quv4scff1QO5+fFPnz4cIwcObJALBUrVjS7rCXVjpYyN21Lk0pLeHp64vDhw9i+fTs2b96MSZMmYcqUKThw4AB8fHzw22+/Yc+ePdi8eTM+++wzvPPOO9i/fz/Cw8MLTKt8+fIAgMTERPj7+1sdi7e3N6pUqQIAqFKlCr788ksEBwfj+++/x9ChQ5V6CQkJJgkWUXEwmaBHlp+fHzp16oR58+Zh5MiRBb6wk5KS4OPjgyeffBLXr1+HRqNBWFhYseen0+lgNBpNyp588kmsXr0aYWFh0Ghs2xzr1q0LX19fTJ8+HfXr14eHhwfatGmDDz/8EImJiUp/CQDYvXs3mjVrhtdee00py390IM+AAQPw7bffokKFClCpVHjqqadMYj958qTyhfUg9mjHypUrQ6vV4sCBA0rCkpycjDNnzih9E1BIW9tDzZo1sXv3bpOy3bt3o1q1alCr1QAAjUaDqKgoREVFYfLkyfDx8cHvv/+OXr16QZIkNG/eHM2bN8ekSZNQqVIl/Pjjjxg7dmyBeUVERMDLywsnT55EtWrVlDKtVov9+/cry5+YmIgzZ86gdevWRcaeF1/eEY88J06cMOkkSlQcPM1Bj7R58+bBaDSicePGWL16Nc6ePYtTp05h7ty5iIyMBABERUUhMjISPXr0wObNmxEfH489e/bgnXfewcGDBy2eV1hYGPbv34/4+Hjcvn0bsiwjJiYGCQkJ6N+/Pw4cOIDz589j06ZNGDx4sNVfhpIkoVWrVvj222+VxKFevXrIzs7G1q1bTb5sqlatioMHD2LTpk04c+YM3n33XRw4cKDANAcMGIDDhw/j/fffx7PPPgu9Xq8Me+utt7Bnzx6MGDECR48exdmzZ/HTTz8V2gHTHu3o6emJQYMGYfz48di2bRv+/vtvDBkyBCqVyuSqjbCwMOzcuRNXrlwxuTLGVm+88Qa2bt2K9957D2fOnMGyZcvw+eefK504161bh7lz5+Lo0aO4ePEi/vvf/0KWZVSvXh379+/HBx98gIMHD+LSpUtYs2YNbt26hZo1a5qdl0qlQlRUFHbt2qWUeXh4YMiQIRg/fjx+//13nDhxAtHR0VCpCu7KMzIycP36dVy/fh3Hjh3Dq6++ChcXF3Ts2FGpEx8fjytXriAqKspubUSPJiYT9EirXLkyDh8+jLZt2+KNN95AnTp10KFDB2zduhULFiwA7n1Jr1+/Hq1atcLgwYNRrVo19OvXDxcvXkRgYKDF8xo3bhzUajVq1aoFf39/XLp0CSEhIdi9ezeMRiM6duyIunXrYvTo0fDx8TH7BfEgrVu3htFoVJIJlUqFVq1aKb+I8wwfPhy9evXCc889hyZNmuDOnTsmRynyVKlSBY0bN8bx48eVqzjy1KtXDzt27MCZM2fQsmVLPPHEE5g0aRJCQkLMxmavdvzkk08QGRmJp59+GlFRUWjevDlq1qxp0rdg2rRpiI+PR0RERLFOERTmySefxA8//IAVK1agTp06mDRpEqZNm4bo6Gjg3umaNWvWoF27dqhZsyYWLlyI7777DrVr14aXlxd27tyJrl27olq1avi///s/zJo1C126dCl0fkOHDsWKFStMTvPMnDkTLVu2RLdu3RAVFYUWLVqY7VOzePFiBAcHIzg4GG3btsXt27exfv16k/4V3333HTp27Fis+1gQ5ScJe54wJCIqZenp6Xjssccwa9YsDBkyxNHh2JUQAk2aNMGYMWPQv39/u047JycHVatWxfLlywt0KiWyFo9MEFGZcuTIEXz33Xc4f/48Dh8+rBwxeeaZZxwdmt1JkoRFixaVyB0qL126hLfffpuJBNkFj0wQUZly5MgRDB06FKdPn4ZOp0ODBg3wySefoG7duo4OjeiRxWSCiIiIbMLTHERERGQTJhNERERkEyYTREREZBMmE0RERGQTJhNERERkEyYTREREZBMmE0RERGQTJhNERERkk/8HLdWyuvlWeaUAAAAASUVORK5CYII=",
      "text/plain": [
       "<Figure size 600x400 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, ax = plt.subplots(figsize=(6, 4))\n",
    "bins = \"auto\"\n",
    "colors = {\n",
    "    \"nominal\": \"tab:blue\",\n",
    "    \"robust\": \"tab:green\",\n",
    "}\n",
    "\n",
    "for label, result in design_results.items():\n",
    "    losses_db = linear_to_loss_db(result[\"samples\"])\n",
    "    ax.hist(\n",
    "        losses_db,\n",
    "        bins=bins,\n",
    "        alpha=0.6,\n",
    "        label=f\"{label.capitalize()} design\",\n",
    "        color=colors.get(label, None),\n",
    "        edgecolor=\"white\",\n",
    "    )\n",
    "    nominal_loss = linear_to_loss_db([result[\"nominal\"]])[0]\n",
    "    ax.axvline(\n",
    "        nominal_loss,\n",
    "        color=colors.get(label, None),\n",
    "        linestyle=\"--\",\n",
    "        linewidth=2,\n",
    "    )\n",
    "\n",
    "ax.set_xlabel(\"Center wavelength loss (dB)\")\n",
    "ax.set_ylabel(\"Sample count\")\n",
    "ax.set_title(\"Monte Carlo comparison at shared perturbations\")\n",
    "ax.legend()\n",
    "ax.grid(alpha=0.25)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ff18ff1a",
   "metadata": {},
   "source": [
    "## What the numbers say\n",
    "\n",
    "Both designs were tested under identical Monte Carlo perturbations (N = 100, σₒᵥₑᵣₗₐᵧ = 25 nm, σₛₚₐcₑᵣ = 20 nm, σ_wᵢdₜₕ = 10 nm) and evaluated at the center wavelength.\n",
    "\n",
    "**Results:**\n",
    "\n",
    "* **Average loss:** Robust 2.51 dB vs nominal 2.56 dB (Δ = −0.05 dB). In linear scale, that’s 0.561 vs 0.555, or about **+1.1 % higher mean transmission.**\n",
    "* **Variability:** Standard deviation (linear) increases slightly (0.027 -> 0.028, **+3 %**), suggesting a comparable level of fluctuation between samples.\n",
    "* **Spread (10th–90th percentile):** 0.0707 -> 0.0755 (**+7 %**) - a slightly broader distribution.\n",
    "* **Tails:**\n",
    "  90th-percentile loss improves (2.86 -> 2.82 dB, **better worst-case**).\n",
    "  10th-percentile loss worsens (2.31 -> 2.23 dB, **slightly lower best-case**).\n",
    "\n",
    "**In short:**\n",
    "The robust design maintains essentially the same overall spread but shifts the entire distribution slightly toward lower loss. While variability remains comparable, the robust version delivers **a modest boost in average transmission and improved worst-case performance**, at the cost of a marginally weaker best-case - a balanced, realistic outcome consistent with fabrication-aware optimization."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "45517bad",
   "metadata": {},
   "source": [
    "At first glance, the numbers may not seem dramatic, but the difference is real and meaningful.  \n",
    "Even a few hundredths of a decibel can translate to higher wafer-level yield when scaled to thousands of devices.  \n",
    "It’s also worth remembering that the specific magnitude depends on many details of the experiment:\n",
    "\n",
    "- How and when robustness was introduced into the optimization (for example, from the start or as a final fine-tuning).  \n",
    "- The starting point, optimizer settings, and number of iterations used.  \n",
    "- The perturbation model and its assumed standard deviations or correlations.  \n",
    "- The type of device. Grating couplers are quite resonant and inherently sensitive to fabrication noise, so they tend to show smaller relative gains.\n",
    "\n",
    "This notebook is meant as a conceptual demonstration rather than an exhaustive benchmark.  \n",
    "There are many other ways to define and train for robustness, and exploring them is part of what makes photonic inverse design both challenging and exciting."
   ]
  }
 ],
 "metadata": {
  "jupytext": {
   "cell_metadata_filter": "-all",
   "encoding": "# coding: utf-8",
   "executable": "/usr/bin/env python",
   "main_language": "python",
   "notebook_metadata_filter": "-all"
  },
  "kernelspec": {
   "display_name": ".venv",
   "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.9"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}