{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "cell-title",
   "metadata": {},
   "source": [
    "# Flux Spectrum and E-Index Analysis\n",
    "\n",
    "This example demonstrates how to attach a neutron flux spectrum to a `Benchmark` using a Serpent detector file, visualize it with `FluxSpectrum.plot_spectrum`, and compute the E-index similarity matrix via `AssimilationSuite.e_index_matrix`. The E-index quantifies the cosine similarity between sensitivity profiles without the influence of a nuclear data covariance matrix."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cell-imports-md",
   "metadata": {},
   "source": [
    "### 1. Setup and Imports\n",
    "We import `andalus` for the core logic and `matplotlib` for visualization."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "cell-imports",
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "\n",
    "import andalus"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cell-reading-md",
   "metadata": {},
   "source": [
    "### 2. Reading Serpent Output\n",
    "We load the **HMF001** (Godiva) benchmark using `Benchmark.from_serpent`. The optional `flux_det` argument attaches a `FluxSpectrum` by reading a named detector from a Serpent `_det0.m` file. Here we pass the detector name `\"flux\"` and the path to the detector output file.\n",
    "\n",
    "The `pertlist` filters the sensitivity output to the most physically relevant reactions:\n",
    "* `mt 2 xs`: Elastic scattering\n",
    "* `mt 18 xs`: Fission\n",
    "* `mt 102 xs`: Radiative capture\n",
    "* `nubar prompt`: Average prompt neutron multiplicity\n",
    "* `chi prompt`: Prompt fission neutron spectrum"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "cell-read-benchmark",
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Reading ../data/hmf001.ser_res.m\n",
      "WARNING:  SERPENT Serpent 2.2.1 found in ../data/hmf001.ser_res.m, but version 2.1.31 is defined in settings\n",
      "WARNING:    Attempting to read anyway. Please report strange behaviors/failures to developers.\n",
      "  - done\n",
      "Reading ../data/hmf001.ser_sens0.m\n",
      "  - done\n",
      "Reading ../data/hmf001.ser_det0.m\n",
      "  - done\n"
     ]
    }
   ],
   "source": [
    "hmf001 = andalus.Benchmark.from_serpent(\n",
    "    title=\"HMF001\",\n",
    "    m=1.0,\n",
    "    dm=0.001,\n",
    "    sens0_path=\"../data/hmf001.ser_sens0.m\",\n",
    "    results_path=\"../data/hmf001.ser_res.m\",\n",
    "    pertlist=[\"mt 2 xs\", \"mt 18 xs\", \"mt 102 xs\", \"nubar prompt\", \"chi prompt\"],\n",
    "    flux_det=(\"flux\", \"../data/hmf001.ser_det0.m\"),\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cell-flux-inspect-md",
   "metadata": {},
   "source": [
    "The `flux` attribute is a `FluxSpectrum`, a pandas DataFrame subclass with a two-level `(E_min_eV, E_max_eV)` MultiIndex. We inspect the first few rows to verify the energy grid and flux values."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "cell-flux-head",
   "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></th>\n",
       "      <th>flux</th>\n",
       "      <th>flux_std</th>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>E_min_eV</th>\n",
       "      <th>E_max_eV</th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0.00001</th>\n",
       "      <th>0.10000</th>\n",
       "      <td>0.000000e+00</td>\n",
       "      <td>0.000000e+00</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>0.10000</th>\n",
       "      <th>0.54000</th>\n",
       "      <td>0.000000e+00</td>\n",
       "      <td>0.000000e+00</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>0.54000</th>\n",
       "      <th>4.00000</th>\n",
       "      <td>0.000000e+00</td>\n",
       "      <td>0.000000e+00</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4.00000</th>\n",
       "      <th>8.31529</th>\n",
       "      <td>5.409640e-09</td>\n",
       "      <td>5.409640e-09</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>8.31529</th>\n",
       "      <th>13.70960</th>\n",
       "      <td>0.000000e+00</td>\n",
       "      <td>0.000000e+00</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "                           flux      flux_std\n",
       "E_min_eV E_max_eV                            \n",
       "0.00001  0.10000   0.000000e+00  0.000000e+00\n",
       "0.10000  0.54000   0.000000e+00  0.000000e+00\n",
       "0.54000  4.00000   0.000000e+00  0.000000e+00\n",
       "4.00000  8.31529   5.409640e-09  5.409640e-09\n",
       "8.31529  13.70960  0.000000e+00  0.000000e+00"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "hmf001.flux.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cell-ealf-md",
   "metadata": {},
   "source": [
    "### 3. Flux Spectrum Characterization\n",
    "Two scalar metrics summarize where neutrons are most active in energy:\n",
    "\n",
    "$$\n",
    "\\text{EALF} = \\exp\\!\\left(\\frac{\\sum_i \\phi_i \\ln \\bar{E}_i}{\\sum_i \\phi_i}\\right)\n",
    "$$\n",
    "\n",
    "$$\n",
    "\\langle E \\rangle = \\frac{\\sum_i \\phi_i \\bar{E}_i}{\\sum_i \\phi_i}\n",
    "$$\n",
    "\n",
    "where $\\bar{E}_i = \\sqrt{E_{\\min,i} \\cdot E_{\\max,i}}$ is the geometric mean energy of bin $i$. EALF is the flux-weighted **geometric** mean (sensitive to the lethargy distribution), while $\\langle E \\rangle$ is the flux-weighted **arithmetic** mean. For a fast spectrum such as HMF001 we expect both to lie well above the thermal region."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "cell-ealf",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "EALF:        8.474e+05 eV\n",
      "Mean energy: 1.461e+06 eV\n"
     ]
    }
   ],
   "source": [
    "print(f\"EALF:        {hmf001.flux.ealf:.3e} eV\")\n",
    "print(f\"Mean energy: {hmf001.flux.mean_energy:.3e} eV\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cell-plot-md",
   "metadata": {},
   "source": [
    "### 4. Plotting the Flux Spectrum\n",
    "We plot the normalized flux per unit lethargy, $\\phi(E)/\\Delta u$, on a logarithmic energy axis. The shaded band represents the $1\\sigma$ statistical uncertainty from the Serpent tally. The sharp peak in the fast region is characteristic of the bare HEU Godiva configuration."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "cell-plot-flux",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAf8AAAGbCAYAAADKlJnyAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjkuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8hTgPZAAAACXBIWXMAAA9hAAAPYQGoP6dpAABLpUlEQVR4nO3de1xUdf4/8NcZdIa7CgwgioJCljcwUbLyViRqtVpp6lYqmpYblctqSZnkLbzFamWy2WLWZmlbuT2+7uIaSWrildTUNHVVvHERBBSUAc75/eGPyYFhmDlzBgbO6/l48Mj5fD7nfd7nfDq8mXPOnBEkSZJAREREqqFp6gSIiIiocbH4ExERqQyLPxERkcqw+BMREakMiz8REZHKsPgTERGpDIs/ERGRyrRq6gSaA1EUcfnyZXh5eUEQhKZOh4iIyCxJknD9+nUEBQVBo6n//T2LvxUuX76M4ODgpk6DiIjIKhcuXEDHjh3r7Wfxt4KXlxcA4Pz582jbtm2jrVcURRQUFECv11v8C07pGNYuY2mcnD5z7bXblNgncii1XlvjKDEXlvrtmQs526MUZz425M6FpT4eG/LHq+3YKC4uRufOnY11qz4s/laoOdXv7e0Nb2/vRluvKIq4desWvL297foFZ2sMa5exNE5On7n22m1K7BM5lFqvrXGUmAtL/fbMhZztUYozHxty58JSH48N+ePVeGwAaPASNW/4IyIiUhkWfyIiIpVh8SciIlIZXvNXUHV1NSorKxWLJ4oiKisrcevWLbuua9oaw9plLI2T02euvXab3H3SunVruLi4WD2eiKglY/FXgCRJyM3NRXFxseJxRVHE9evXZT9fQE4Ma5exNE5On7n22m327JO2bdsiMDCQz2ogItVj8VdATeH39/eHu7u7YsVFkiRUVVWhVatWdhV/W2NYu4ylcXL6zLXXbpO7PeXl5cjPzwcAtG/f3qrliIhaKhZ/O1VXVxsLv6+vr6KxWfyVKf4A4ObmBgDIz8+Hv78/LwEQkarxhj871Vzjd3d3b+JMqCE1c6TkfRlERM0Ri79CeB3Z+XGOiIhuY/EnIiJSGRZ/IiIileENf0REzYAkSSg3VJk8A+NmZbWxrfZra7m1duElMRVi8VexyZMno7i4GJs3bzZpz8zMxEMPPYSioiIcPnwYQ4cORdu2bXHlyhW4uroax+3fvx/R0dEAbv9iqll26NChddb15ptvYtGiRQCAI0eO4KWXXsKBAweg1+vx8ssv47XXXjMZ/9VXX+Gtt97CuXPnEBYWhqVLl+LRRx819n/zzTf429/+hoMHD6KoqAg///wzIiMjldgtRE5HkiRM33QSv1wpUzx2RMc22PzSA/wDQGV42p+s4uXlhW+//dak7e9//zs6depkdvyJEyeQk5ODy5cv48qVK5gzZw4AoLS0FLGxsejUqRMOHDiA5cuX4+2338ZHH31kXHb37t2YMGECpk6diuzsbPzhD3/AE088gaNHjxrHlJWV4cEHH8TSpUsdsLVEzuVmZbVDCj8AHL5YgovXbjokNjkvvvN3AEmScLOyWpE4VVVVaCVad6e6I0/fTZo0CWlpaZgwYQIA4ObNm9i4cSNeeuklvPPOO3XG+/v7w9PTs87n8T///HMYDAasXbsW7u7u6NmzJw4dOoSUlBRMmzYNAPDee+9h+PDhmD17NiRJwvz58/HDDz/ggw8+wJo1awAAzz33HARBwLlz5xyyvUTOat3kKAT7eECSRJQWFcLbxxeCoKnzuiG3Kqvx2Pu7AAAVVaKj0yYnw+LvADcrq9F93tZGX+/xBbFw1zpmSp977jksX74cOTk56NSpE77++muEhISgT58+NsXJysrCoEGDoNVqjW2xsbFYunQprl27Bi8vL2RlZSEhIcFkuWHDhuFf//qXIttC1JwF+3ggzN8ToigiXyqHv97TeM3/ztcNKTdUNUK25Kx42l/l/u///g+enp4mPyNHjqwzzt/fHyNGjMAnn3wCAEhLS0NcXFy9cYODg9GuXTt4eXnB09MThYWFAG4/Ctnf399kbEBAgLGv5r81bXeOqeknIiL78J2/A7i1dsHxBbF2x7H1UbZurW1/ZO3QoUONp9Jr7NmzB88991ydsVOmTMGrr76KZ599FllZWdi0aRMyMzPNxt2xYwfc3NyMubdr187m3IiIyDFY/B1AEARFTr9LkoQqDex6tn9DPDw8EBYWZtJ24cIFs2NHjBiB6dOnY+rUqXj88cctfpdBaGio2Wv+gYGBxi/YqZGXl2fsq/lvTdudY2r6iYjIPjztT1Zr1aoVJk6ciMzMTEyZMkVWjAEDBmDHjh0mz9fftm0bunXrZjw7MGDAAGRkZJgs9/3332PAgAHykyciIiMWf7LJwoULUVBQgNhYeZc1/vjHP0Kr1WL69Ok4duwYNm7ciFWrVpnc4PfKK68gPT0d7777Lk6cOIEFCxbgwIEDiI+PN44pKirCoUOHcPz4cQDAyZMncejQId4XQERkBRZ/solWq4Wfn5/syxBt2rTB1q1bce7cOURFReEvf/kL5s2bh+nTpxvH3H///diwYQM++ugjREZG4ptvvsG3336Lnj17Gsd899136NOnj/HBP+PHj0efPn2Qmppq3wYSEakAr/mrWM2d+7UNGTIEBoMBrVq1wpAhQ4xP7zNn1KhREEXRZFlJkow3K5rTu3dvbN++3eK9DGPHjsXYsWNNbnq80+TJky1+2oCIiOrHd/5EREQqw3f+REQqV1hWAc9S2z8qDACSKKGorAJS6S0Imt/P5HnpWsFD11qpFElhLP5ERCo37m97FI8Z6ueB/3v5Af4B4KR42p+ISIXcWrvg3k5tHRb/7NUyFFw3OCw+2Yfv/BVy501v5Jw4R0S/EwQBX8+4H4VlFagS67+ptyGSKKGosAA+vnoIGgE3DdUYuuJHALArLjkWi7+dtFotNBoNLl++DL1eD61Wq9jT+Gx9vK9SMaxdxtI4OX3m2mu3yd0eg8GAgoICaDQaky8VIlIzQRDg5+lqVwxRFCHc0sHf2xUajYZfGNRMsPjbSaPRIDQ0FFeuXMHly5cVjS1JEkRRhEajsav42xrD2mUsjZPTZ669dps9+8Td3R2dOnWy6hvPiIhaMhZ/BWi1WnTq1AlVVVWorq5WLK4oiigsLISvr6/sgiUnhrXLWBonp89ce+02ufvExcXFod+RQETUnLD4K0QQBLRu3RqtWyt3Z6soimjdujVcXV3tKv62xrB2GUvj5PSZa6/dpsQ+IXIkSZJws7Ia5YYqs/+PiqJYb399feUG5d5UEAEs/kREipEkCU//bQ8O5hQ3dSpEFvGtExGRQm5WVju08IfpPdDGje/ZyH78v4iIyAHSJkWhk69HnXZJElFaVAhvH18Igsaqvpr24KAA6L3suzufCGDxJyJyiGAfN4T5e9ZpF0UR+VI5/PWeZq/5m+urafdj4SeF8LQ/ERGRyrD4ExERqQyLPxERkco4ZfFfvXo1QkJC4OrqiujoaOzbt6/esd988w2ioqLQtm1beHh4IDIyEp999pnJGEmSMG/ePLRv3x5ubm6IiYnBqVOnHL0ZRERETsnpiv/GjRuRkJCApKQkZGdnIyIiArGxscjPzzc73sfHB2+++SaysrJw5MgRxMXFIS4uDlu3bjWOWbZsGd577z2kpqZi79698PDwQGxsLG7dutVYm0VEROQ0nK74p6SkYNq0aYiLi0P37t2RmpoKd3d3pKWlmR0/ZMgQPPHEE7jnnnvQtWtXvPrqq+jduzd27doF4Pa7/pUrV2Lu3LkYNWoUevfujU8//RSXL1/G5s2bG3HLiIiInINTfdTPYDDg4MGDSExMNLZpNBrExMQgKyurweUlScIPP/yAkydPYunSpQCAs2fPIjc3FzExMcZxbdq0QXR0NLKysjB+/Pg6cSoqKlBRUWF8XVpaCuD2x20a82thRVE0fpFNY8awdhlL4+T0mWuv3abEPpFDqfXaGkeJubDUb89cyNkepTjrsXFnuySZ/33R0o8Na/aBNXGUGq/GY8MaTlX8r169iurqagQEBJi0BwQE4MSJE/UuV1JSgg4dOqCiogIuLi748MMP8cgjjwAAcnNzjTFqx6zpqy05ORnz58+v015QUACDwWDTNtlDFEWUlJRAkiS7nu1vawxrl7E0Tk6fufbabUrsEzmUWq+tcZSYC0v99syFnO1RirMeGzcrf38Gf+m1IuSj7qXFln5smOyDokLkS+Wy4ig1Xm3HRklJiVXjnKr4y+Xl5YVDhw7hxo0byMjIQEJCArp06YIhQ4bIipeYmIiEhATj69LSUgQHB0Ov16Nt27bKJG0FURQhCAL0er1dv+BsjWHtMpbGyekz1167TYl9IodS67U1jhJzYanfnrmQsz1KcdZj487vsvdu5wN/f2+bYrSEY8NkH/j4wl9f90FH1sRRarzajg2tVmvVOKcq/n5+fnBxcUFeXp5Je15eHgIDA+tdTqPRICwsDAAQGRmJX3/9FcnJyRgyZIhxuby8PLRv394kZmRkpNl4Op0OOp3O7Hoa+5vkBEGwe71yYli7jKVxcvrMtdduU2KfyKHUem2No8RcWOq3Zy7kbI9SnPHYMN0vts+Fpb7mcmxYuw8aiqPkeDUdG9auy6lu+NNqtejbty8yMjKMbaIoIiMjAwMGDLA6jiiKxmv2oaGhCAwMNIlZWlqKvXv32hSTiIiopXCqd/4AkJCQgEmTJiEqKgr9+/fHypUrUVZWhri4OADAxIkT0aFDByQnJwO4fX0+KioKXbt2RUVFBf7973/js88+w5o1awDc/str5syZWLRoEcLDwxEaGoq33noLQUFBGD16dFNtJhERUZNxuuI/btw4FBQUYN68ecjNzUVkZCTS09ONN+zl5OSYnNYoKyvDn/70J1y8eBFubm64++678Y9//APjxo0zjnnttddQVlaG6dOno7i4GA8++CDS09Ph6sovySAiIvVxuuIPAPHx8YiPjzfbl5mZafJ60aJFWLRokcV4giBgwYIFWLBggVIpEhERNVtOdc2fiIiIHI/Fn4iISGVY/ImIiFSGxZ+IiEhlWPyJiIhUhsWfiIhIZVj8iYiIVIbFn4iISGVY/ImIiFSGxZ+IiEhlWPyJiIhUhsWfiIhIZVj8iYiIVIbFn4iISGVY/ImIiFSGxZ+IiEhlWPyJiIhUhsWfiIhIZVj8iYiIVIbFn4iISGVY/ImIiFSGxZ+IiEhlWPyJiIhUhsWfiIhIZVj8iYiIVIbFn4iISGVY/ImIiFSmVVMnQERELdOtymqUG6qsGiuKIm7+//EaTf3vS91au0AQBKVSVC0WfyIicojH3t+leMyIjm2w+aUHFI+rNjztT0REinFr7YKIjm0cFv/wxRJcvHbTYfHVgu/8iYhIMYIgYPNLD+DitZuoqBKtXk6SRJQWFcLbxxeCUPd96a3KauOZBFviknks/kREpChBEBDs427TMqIoIl8qh7/e0+w1f2vvHSDr8LQ/ERGRyrD4ExERqQyLPxERkcqw+BMREamMUxb/1atXIyQkBK6uroiOjsa+ffvqHbt27VoMHDgQ7dq1Q7t27RATE1Nn/OTJkyEIgsnP8OHDHb0ZRERETsnpiv/GjRuRkJCApKQkZGdnIyIiArGxscjPzzc7PjMzExMmTMD27duRlZWF4OBgDBs2DJcuXTIZN3z4cFy5csX488UXXzTG5hARETkdpyv+KSkpmDZtGuLi4tC9e3ekpqbC3d0daWlpZsd//vnn+NOf/oTIyEjcfffd+PjjjyGKIjIyMkzG6XQ6BAYGGn/atWvXGJtDRETkdJzqc/4GgwEHDx5EYmKisU2j0SAmJgZZWVlWxSgvL0dlZSV8fHxM2jMzM+Hv74927drhoYcewqJFi+Dr62s2RkVFBSoqKoyvS0tLAdz+HKooNt7DJURRhCRJdq1TTgxrl7E0Tk6fufbabUrsEzmUWq+tcZSYC0v99syFnO1RirMeG3e2S5L53xc8NuTFMbdveWzUZe36nKr4X716FdXV1QgICDBpDwgIwIkTJ6yK8frrryMoKAgxMTHGtuHDh+PJJ59EaGgozpw5gzfeeAMjRoxAVlYWXFxc6sRITk7G/Pnz67QXFBTAYDDYuFXyiaKIkpISSJJk8YsulI5h7TKWxsnpM9deu02JfSKHUuu1NY4Sc2Gp3565kLM9SnHWY+NmZbXx36XXipCPWzbF4LFh5b4tKkR+9Q0eG2aUlJRYNc6pir+9lixZgi+//BKZmZlwdXU1to8fP9747169eqF3797o2rUrMjMz8fDDD9eJk5iYiISEBOPr0tJSBAcHQ6/Xo23btg7dhjuJoghBEKDX6+36BWdrDGuXsTROTp+59tptSuwTOZRar61xlJgLS/32zIWc7VGKsx4bdz6FzrudD/z9vW2KwWPDyn3r4wt/X3ceG2ZotVqrxjlV8ffz84OLiwvy8vJM2vPy8hAYGGhx2RUrVmDJkiX4/vvv0bt3b4tju3TpAj8/P5w+fdps8dfpdNDpdHXaNRpNo04icPsxmfauV04Ma5exNE5On7n22m1K7BM5lFqvrXGUmAtL/fbMhZztUYpS6y0qr4S1J2YlUUJRuQG4YYCgqfu1sjcNv787FQTb58JSn9qPDdPtvj2Gx0Zd1q7LqYq/VqtF3759kZGRgdGjRwO4/ddTRkYG4uPj611u2bJlWLx4MbZu3YqoqKgG13Px4kUUFhaiffv2SqVORM2MJEmYtvEkjuaWOSR+K37nPDkxp7vbPyEhAWvXrsX69evx66+/YsaMGSgrK0NcXBwAYOLEiSY3BC5duhRvvfUW0tLSEBISgtzcXOTm5uLGjRsAgBs3bmD27NnYs2cPzp07h4yMDIwaNQphYWGIjY1tkm0koqZ3s7LaYYU/tJ0Ofl7WnX4lagpO9c4fAMaNG4eCggLMmzcPubm5iIyMRHp6uvEmwJycHJPTGmvWrIHBYMCYMWNM4iQlJeHtt9+Gi4sLjhw5gvXr16O4uBhBQUEYNmwYFi5caPbUPhGpz8YX7kNn34a/hU4SJRQVFsDHV2/2tH/NmIrSInjoWiudJpFinK74A0B8fHy9p/kzMzNNXp87d85iLDc3N2zdulWhzIioJfL10CHQ263BcaIoQrilg7+3a73XVkVRRP4tp/zVSmRk1/+hlZWVyM3NRXl5OfR6fZ3P1hMREZHzsfma//Xr17FmzRoMHjwY3t7eCAkJwT333AO9Xo/OnTtj2rRp2L9/vyNyJSIiIgXYVPxTUlIQEhKCdevWISYmBps3b8ahQ4fw22+/ISsrC0lJSaiqqsKwYcMwfPhwnDp1ylF5ExERkUw2nfbfv38/duzYgR49epjt79+/P6ZMmYLU1FSsW7cOO3fuRHh4uCKJEhERkTJsKv7WfhOeTqfDiy++KCshIiIiciyn+5w/EREROZZN7/zvfN59Q1JSUmxOhoiIiBzPpuL/888/m7zOzs5GVVUVunXrBgD47bff4OLigr59+yqXIRERESnKpuK/fft2479TUlLg5eWF9evXo127dgCAa9euIS4uDgMHDlQ2SyIiIlKM7Gv+7777LpKTk42FHwDatWuHRYsW4d1331UkOSIiIlKe7OJfWlqKgoKCOu0FBQW4fv26XUkRERGR48gu/k888QTi4uLwzTff4OLFi7h48SK+/vprTJ06FU8++aSSORIREZGCZD/bPzU1FbNmzcIf//hHVFZW3g7WqhWmTp2K5cuXK5YgERERKUt28Xd3d8eHH36I5cuX48yZMwCArl27wsPDQ7HkiIiISHl2PeRn586deOGFF/Diiy/C19cXHh4e+Oyzz7Br1y6l8iMiIiKFyS7+X3/9NWJjY+Hm5obs7GxUVFQAAEpKSvDOO+8oliAREREpS3bxX7RoEVJTU7F27Vq0bt3a2P7AAw8gOztbkeSIiIhIebKL/8mTJzFo0KA67W3atEFxcbE9OREREZEDyS7+gYGBOH36dJ32Xbt2oUuXLnYlRURERI4ju/hPmzYNr776Kvbu3QtBEHD58mV8/vnnmDVrFmbMmKFkjkRERKQg2R/1mzNnDkRRxMMPP4zy8nIMGjQIOp0Os2bNwssvv6xkjkRERKQg2cVfEAS8+eabmD17Nk6fPo0bN26ge/fu8PT0VDI/IiIiUpjs4l9Dq9Wie/fuSuRCREREjUD2Nf+HHnoI8+fPr9N+7do1PPTQQ3YlRURERI4j+51/ZmYmfvnlF/z888/4/PPPjY/1NRgM+PHHHxVLkIiIiJRl1+N9v//+e+Tm5uK+++7DuXPnFEqJiIiIHMmu4t++fXv8+OOP6NWrF/r164fMzEyF0iIiIiJHkV38BUEAAOh0OmzYsAGvvvoqhg8fjg8//FCx5IiIiEh5sq/5S5Jk8nru3Lm45557MGnSJLuTIiIiIseRXfzPnj0LPz8/k7annnoK3bp1w8GDB+1OjIiIiBxDdvHv3Lmz2faePXuiZ8+eshMiIiIix7Kp+CckJGDhwoXw8PBAQkKCxbEpKSl2JUZERESOYVPx//nnn1FZWWn8NxERETU/NhX/7du3m/03ERERNR82n/a3hiAIePfdd2UlRERERI5l82l/a9Q8A4CIiIicj+zT/o60evVqLF++HLm5uYiIiMD777+P/v37mx27du1afPrppzh69CgAoG/fvnjnnXdMxkuShKSkJKxduxbFxcV44IEHsGbNGoSHhzfK9hARETkTux7v6wgbN25EQkICkpKSkJ2djYiICMTGxiI/P9/s+MzMTEyYMAHbt29HVlYWgoODMWzYMFy6dMk4ZtmyZXjvvfeQmpqKvXv3wsPDA7Gxsbh161ZjbRYREZHTkP05/xrHjx9HTk4ODAaDSfsf/vAHWfFSUlIwbdo0xMXFAQBSU1OxZcsWpKWlYc6cOXXGf/755yavP/74Y3z99dfIyMjAxIkTIUkSVq5ciblz52LUqFEAgE8//RQBAQHYvHkzxo8fLytPIiKi5kp28f/f//6HJ554Ar/88gsEQTA+7rfmen91dbXNMQ0GAw4ePIjExERjm0ajQUxMDLKysqyKUV5ejsrKSvj4+AC4/STC3NxcxMTEGMe0adMG0dHRyMrKMlv8KyoqUFFRYXxdWloKABBFEaIo2rxdcomiCEmS7FqnnBjWLmNpnJw+c+2125TYJ3IotV5b4ygxF5b67ZkLOdujFKWOjRqSZN2xbc165c6FpT4eG+bni8dGXdauT3bxf/XVVxEaGoqMjAyEhoZi3759KCwsxF/+8hesWLFCVsyrV6+iuroaAQEBJu0BAQE4ceKEVTFef/11BAUFGYt9bm6uMUbtmDV9tSUnJ2P+/Pl12gsKCuqc4XAkURRRUlICSZKg0ci7QiMnhrXLWBonp89ce+02JfaJHEqt19Y4SsyFpX575kLO9ihFifWWVVQa/11aVIh8qVyR9cqdC0t9PDaAm5W/v6EsLSpEfvUNHhtmlJSUWDVOdvHPysrCDz/8AD8/P2g0Gmg0Gjz44INITk7GK6+80iQPAVqyZAm+/PJLZGZmwtXVVXacxMREk481lpaWIjg4GHq9Hm3btlUgU+uIoghBEKDX6+0q/rbGsHYZS+Pk9Jlrr92mxD6RQ6n12hpHibmw1G/PXMjZHqUosd4bt37/Q97bxxf+ek9F1it3Liz18dgAyg1Vxn97+/jC39edx4YZWq3WqnGyi391dTW8vLwAAH5+frh8+TK6deuGzp074+TJk7Ji+vn5wcXFBXl5eSbteXl5CAwMtLjsihUrsGTJEnz//ffo3bu3sb1muby8PLRv394kZmRkpNlYOp0OOp2uTnvNHzmNSRAEu9crJ4a1y1gaJ6fPXHvtNiX2iRxKrdfWOErMhaV+e+ZCzvYoxd71mm6DsvMhdy4s9an92DA3Xzw26rJ2XbIz6tmzJw4fPgwAiI6OxrJly/DTTz9hwYIF6NKli6yYWq0Wffv2RUZGhrFNFEVkZGRgwIAB9S63bNkyLFy4EOnp6YiKijLpCw0NRWBgoEnM0tJS7N2712JMIiKilkr2O/+5c+eirKwMALBgwQI89thjGDhwIHx9fbFx40bZCSUkJGDSpEmIiopC//79sXLlSpSVlRnv/p84cSI6dOiA5ORkAMDSpUsxb948bNiwASEhIcbr+J6envD09IQgCJg5cyYWLVqE8PBwhIaG4q233kJQUBBGjx4tO08iIqLmSnbxj42NNf47LCwMJ06cQFFREdq1a2fXE/7GjRuHgoICzJs3D7m5uYiMjER6errxhr2cnByT0xpr1qyBwWDAmDFjTOIkJSXh7bffBgC89tprKCsrw/Tp01FcXIwHH3wQ6enpdt0XQERE1FzZ9Tn/W7du4ciRI8jPz6/z8QK5n/MHgPj4eMTHx5vty8zMNHl97ty5BuMJgoAFCxZgwYIFsnMiIiJqKWQX//T0dDz33HMoLCys0ycIgqzP+RMREZHjyb7h7+WXX8bTTz+NK1euQBRFkx8WfiIiIuclu/jn5eUhISGhzsNziIiIyLnJLv5jxoypc/2diIiInJ/sa/4ffPABxo4di507d6JXr15o3bq1Sf8rr7xid3JERESkPNnF/4svvsB///tfuLq6IjMz0+TjfYIgsPgTERE5KdnF/80338T8+fMxZ86cRn+UJBEREcknu2obDAaMGzeOhZ+IiKiZkV25J02aZNdjfImIiKhp2PWtfsuWLcPWrVvRu3fvOjf8paSk2J0cERERKU928f/ll1/Qp08fAMDRo0dN+ux5tj8RERE5luziv337diXzICIiokYi65p/ZWUlHn74YZw6dUrpfIiIiMjBZBX/1q1b48iRI0rnQkRERI1A9t3+zz77LP7+978rmQsRERE1AtnX/KuqqpCWlobvv/8effv2hYeHh0k/7/YnIiJyTrKL/9GjR3HvvfcCAH777TeTPt7tT0RE5Lx4tz8REZHK8Nm8REREKiP7nT8AFBcX4+9//zt+/fVXAED37t0xdepUtGnTRpHkiIiISHmy3/kfOHAAXbt2xV//+lcUFRWhqKgIf/3rX9G1a1dkZ2crmSMREREpSPY7/z//+c/4wx/+gLVr16JVq9thqqqq8Pzzz2PmzJnYsWOHYkkSERGRcmQX/wMHDpgUfgBo1aoVXnvtNURFRSmSHBERESlPdvH39vZGTk4O7r77bpP2CxcuwMvLy+7EiIhqSJKEm5XVxteiKOJmZTXKDVXQaORdvSw3VDc8iKiFkl38x40bh6lTp2LFihW4//77AQA//fQTZs+ejQkTJiiWIBGpmyRJGJOahYPnrzV1KkQthuziv2LFCgiCgIkTJ6KqqgqSJEGr1WLGjBlYsmSJkjkSkYrdrKx2aOEP03ugjZtdH3wianZk/x+v1WqxatUqJCcn48yZMwCArl27wt3dXbHkiIjutG5yFIJ9PCBJIkqLCuHt4wtBkHfavyZGcFAA9F6uCmdK5Nzs+nM3IyMDGRkZyM/PhyiKJn1paWl2JUZEVFuwjwfC/D0hiiLypXL46z1lX/OvieHHwk8qJLv4z58/HwsWLEBUVBTat2/P5/kTERE1E7KLf2pqKj755BM899xzSuZDREREDib7CX8Gg8F4lz8RERE1H7KL//PPP48NGzYomQsRERE1Atmn/W/duoWPPvoI33//PXr37o3WrVub9KekpNidHBERESlPdvE/cuQIIiMjAQBHjx416ePNf0RERM5LdvHfvn27knkQERFRI5F9zZ+IiIiaJ5uKf05Ojk3BL126ZNN4AFi9ejVCQkLg6uqK6Oho7Nu3r96xx44dw1NPPYWQkBAIgoCVK1fWGfP2229DEASTn9pfRkRERKQmNhX/fv364YUXXsD+/fvrHVNSUoK1a9eiZ8+e+Prrr21KZuPGjUhISEBSUhKys7MRERGB2NhY5Ofnmx1fXl6OLl26YMmSJQgMDKw3bo8ePXDlyhXjz65du2zKi4iIqCWx6Zr/8ePHsXjxYjzyyCNwdXVF3759ERQUBFdXV1y7dg3Hjx/HsWPHcO+992LZsmUYOXKkTcmkpKRg2rRpiIuLA3D7QUJbtmxBWloa5syZU2d8v3790K9fPwAw21+jVatWFv84ICIiUhObir+vry9SUlKwePFibNmyBbt27cL58+dx8+ZN+Pn54ZlnnkFsbCx69uxpcyIGgwEHDx5EYmKisU2j0SAmJgZZWVk2x7vTqVOnjH+kDBgwAMnJyejUqVO94ysqKlBRUWF8XVpaCuD2s8Brf4eBI4miCEmS7FqnnBjWLmNpnJw+c+2125TYJ3IotV5b4ygxF5b67ZkLOdsjx52xJUk0HofOemzInQtLfTw2zP9/oPZjwxxr1yfrbn83NzeMGTMGY8aMkbO4WVevXkV1dTUCAgJM2gMCAnDixAnZcaOjo/HJJ5+gW7duuHLlCubPn4+BAwfi6NGj8PLyMrtMcnIy5s+fX6e9oKAABoNBdi62EkURJSUlkCTJri8vsTWGtctYGienz1x77TYl9okcSq3X1jhKzIWlfnvmQs72yHGzstr479KiQuRL5U59bMidC0t9PDbM/H9QfUP1x4Y5JSUlVo1r8V9iPWLECOO/e/fujejoaHTu3BmbNm3C1KlTzS6TmJiIhIQE4+vS0lIEBwdDr9ejbdu2jk7ZSBRFCIIAvV5v1y84W2NYu4ylcXL6zLXXblNin8ih1HptjaPEXFjqt2cu5GyPHOWGKuO/vX184a/3dOpjQ+5cWOrjsWHm/wNfd9UfG+ZotVqrxjlN8ffz84OLiwvy8vJM2vPy8hS9Xt+2bVvcddddOH36dL1jdDoddDpdnXaNRtOokwjcfmCSveuVE8PaZSyNk9Nnrr12mxL7RA6l1mtrHCXmwlK/PXMhZ3tsZbouZf8/cNSxIXcuLPWp/dgw9/+B2o8Nc6xdl9N8zl+r1aJv377IyMgwtomiiIyMDAwYMECx9dy4cQNnzpxB+/btFYtJRETUnDjNO38ASEhIwKRJkxAVFYX+/ftj5cqVKCsrM979P3HiRHTo0AHJyckAbt8kePz4ceO/L126hEOHDsHT0xNhYWEAgFmzZuHxxx9H586dcfnyZSQlJcHFxQUTJkxomo0kIiJqYrKLf1FREXx8fJTMBePGjUNBQQHmzZuH3NxcREZGIj093XgTYE5OjskpjcuXL6NPnz7G1ytWrMCKFSswePBgZGZmAgAuXryICRMmoLCwEHq9Hg8++CD27NkDvV6vaO5ERETNhezi7+fnhw4dOiAiIsLk56677rLri33i4+MRHx9vtq+moNcICQmBJEkW43355ZeycyEiImqJZBf/X375BYcOHcLhw4exf/9+fPTRRygqKoKrqyt69uyJvXv3KpknERERKUR28e/Rowd69OiBZ555BgAgSRLS09Px8ssv4+GHH1YsQSIiIlKWYnf7C4KAESNG4B//+Adyc3OVCktEREQKU/yjfvfddx+2b9+udFgiIiJSiOzT/p6enujVqxciIiLQu3dvRERE4O6778b+/ftx/fp1JXMkIiIyKiyrgIdWg6KyCkiltyBo6r/JXBIli+MkUULFHU8PVAvZxf+f//wnDh06hEOHDmHVqlU4c+YMJEmCIAhYuHChkjkSEREZjfvbHkXjhbbT4btX9PByq/tk15ZKdvEfPnw4hg8fbnxdXl6Os2fPwtfXl1+fS0REinJr7YJ7O7VFdk6x4rHPXqvA1esGFn853N3d0aNHDwDA0aNHZX2tLxERkTmCIODrGfejsKwCVaJ0+3R+YQF8fPUNn/avZ9xNQzWGrvgRAFDVwDNjWhrFiv/169fxxRdf4OOPP0Z2djaqqtR3DYWIiBxHEAT4eboC+P/fmndLB39v14a/1a+eceUqvNZfw+67/Xfs2IFJkyahffv2mDt3LoKDgxt86h4RERE1HVnFPzc3F0uWLEF4eDhGjhyJqqoqbNq0CZcvX8b8+fOVzpGIiIgUZPNp/8cffxwZGRkYOnQo3n77bYwePRoeHh7Gfnue609ERESOZ3Px37JlC/74xz9i5syZiIqKckRORERE5EA2n/bfvXs33Nzc8NBDD6Fbt25YsGABzpw544jciIiIyAFsLv733Xcf1q5diytXruD111/Hf//7X9x1112477778P777yMvL88ReRIREZFCZN/t7+HhgSlTpmDXrl04fvw4Bg0ahHfeeQcxMTFK5kdEREQKU+SLfbp164Zly5bh4sWL+Oabb/Doo48qEZaIiIgcQNFv9XNxccHo0aPx3XffKRmWiIiIFKT4V/oSERGRc7O5+M+bNw8HDx50RC5ERETUCGwu/hcvXsSIESPQsWNHzJgxA//5z39gMBgckRsRERE5gM3FPy0tDbm5ufjiiy/g5eWFmTNnws/PD0899RQ+/fRTFBUVOSJPIiIiUoisa/4ajQYDBw7EsmXLcPLkSezduxfR0dH429/+hqCgIAwaNAgrVqzApUuXlM6XiIiI7KTIDX/33HMPXnvtNfz000+4cOECJk2ahJ07d+KLL75QIjwREREpyOZn+zdEr9dj6tSpmDp1qtKhiYiISAH8qB8REZHK2PXOv7KyErm5uSgvL4der4ePj49SeREREZGD2PzO//r161izZg0GDx4Mb29vhISE4J577oFer0fnzp0xbdo07N+/3xG5EhERkQJsKv4pKSkICQnBunXrEBMTg82bN+PQoUP47bffkJWVhaSkJFRVVWHYsGEYPnw4Tp065ai8iYiISCabTvvv378fO3bsQI8ePcz29+/fH1OmTEFqairWrVuHnTt3Ijw8XJFEiYiISBk2FX9rP7qn0+nw4osvykqIiIiIHMvma/6TJ09GeXm5I3IhIiKiRmBz8f/ss89w48YN4+sZM2aguLjYZExVVZXdiREREZFj2Fz8JUkyef3555+bPM8/Ly8P3t7e9mdGREREDmH3Q35q/zEAALdu3bI3LBERETmIQ57wJwiCI8ISERGRAmQV/w0bNiA7OxuVlZVK54PVq1cjJCQErq6uiI6Oxr59++ode+zYMTz11FMICQmBIAhYuXKl3TGJiIhaOpuL/8CBA5GUlISoqCh4enqivLwcSUlJSE1NxZ49e0xuBrTVxo0bkZCQgKSkJGRnZyMiIgKxsbHIz883O768vBxdunTBkiVLEBgYqEhMIiKils7m4v/jjz+ipKQEJ0+exPr16/GXv/wFV65cwRtvvIH7778fd911l+xkUlJSMG3aNMTFxaF79+5ITU2Fu7s70tLSzI7v168fli9fjvHjx0On0ykSk4iIqKWT/cU+4eHhCA8Px/jx441tZ8+exYEDB/Dzzz/bHM9gMODgwYNITEw0tmk0GsTExCArK0tWjnJjVlRUoKKiwvi6tLQUACCKIkRRlJWLHKIoQpIku9YpJ4a1y1gaJ6fPXHvtNiX2iRxKrdfWOErMhaV+e+bCXJskSbhZWW3Vtlmr3PB7PEkSjcehsx4bcufCUh+PDfnjG9rfNWr+32poObnHRmOxdn02Ff+cnBx06tSp3v7Q0FCEhoZi7NixAIBLly6hQ4cOVsW+evUqqqurERAQYNIeEBCAEydO2JKm3TGTk5Mxf/78Ou0FBQUwGAyycpFDFEWUlJRAkiRoNPLuzZQTw9plLI2T02euvXabEvtEDqXWa2scJebCUr89c1G7TRAETN90Er9cKbN6f9iqtKgQ+VK5Ux8bcufCUh+PDfnjLY278w/V0mtFyMetBpeTc2w05nyUlJRYNc6m4t+vXz+MHj0azz//PPr161fvijdt2oRVq1Zh+vTpeOWVV2xZhVNITExEQkKC8XVpaSmCg4Oh1+vRtm3bRstDFEUIggC9Xm/XLzhbY1i7jKVxcvrMtdduU2KfyKHUem2No8RcWOq3Zy5qt92qEh1a+MP0HggOCoCfl6tTHxty58JSH48N+eMtjSs3/P5AOu92PvD3925wOTnHRmPOh1artWqcTcX/+PHjWLx4MR555BG4urqib9++CAoKgqurK65du4bjx4/j2LFjuPfee7Fs2TKMHDnS6th+fn5wcXFBXl6eSXteXl69N/M5KqZOpzN7D4FGo2nUSQRuf2zS3vXKiWHtMpbGyekz1167TYl9IodS67U1jhJzYanfnrkwbft92XWToxDs42HV9lmrjVsr6L1cG8zbFo46NuTOhaU+Hhvyx9c3znRfOvLYaLz5sHZdNhV/X19fpKSkYPHixdiyZQt27dqF8+fP4+bNm/Dz88MzzzyD2NhY9OzZ0+aEtVot+vbti4yMDIwePRrA7b+cMjIyEB8fb3M8R8UkooYF+3ggzN+zqdMgonrIuuHPzc0NY8aMwZgxYxRNJiEhAZMmTUJUVBT69++PlStXoqysDHFxcQCAiRMnokOHDkhOTgZw+4a+48ePG/996dIlHDp0CJ6enggLC7MqJhERkdrIvtsfAI4cOYKdO3dCq9XigQceQPfu3e1KZty4cSgoKMC8efOQm5uLyMhIpKenG2/Yy8nJMTmlcfnyZfTp08f4esWKFVixYgUGDx6MzMxMq2ISERGpjeziv2rVKvz5z3+Gt7c3XFxccO3aNfTq1Qvr169HZGSk7ITi4+PrPSVfU9BrhISEmP1uAVtiEhERqY1NdyGkpaUhOzsbFRUVWLx4MZYsWYJr166hsLAQ//vf/zBixAgMHDgQu3fvdlS+REREZCeb3vmvWLECp06dAnD7xrn9+/dj1apV6NOnDyIjI7FkyRIEBwdj1qxZ/AOAiIjISdn0zv/48eO4fv06du/ejdatW0Oj0eDLL7/EyJEj4ePjgy5duuDbb7/FwYMHsWXLFpw7d85BaRMREZFcNn/40NXVFf369cMDDzyAiIgI7NmzB9evX8cvv/yCRYsWISwsDJWVlZg4cSK6dOkCb2/vhoMSERFRo5F9w9+7776LIUOG4H//+x9efPFFREREIDg4GNnZ2QgKCsLFixdx8eJFHD16VMl8iYiIyE6yi39kZCQOHjyIF198Effdd5/xrvtWrVoZvzGvY8eO6NixozKZEhERkSLs+px/165dsW3bNuTl5WHPnj0wGAwYMGAACz4REZETs6v41wgICMCoUaOUCEVEREQO1rjf/kBERERNjsWfiIhIZVj8iYiIVIbFn4iISGVY/ImIiFSGxZ+IiEhlWPyJiIhUhsWfiIhIZVj8iYiIVIbFn4iISGVY/ImIiFSGxZ+IiEhlWPyJiIhUhsWfiIhIZVj8iYiIVIbFn4iISGVY/ImIiFSGxZ+IiEhlWPyJiIhUhsWfiIhIZVj8iYiIVIbFn4iISGVY/ImIiFSGxZ+IiEhlWPyJiIhUhsWfiIhIZVj8iYiIVKZVUydARE1DkiSUG6pws7Ia5YYqaDS/vxcQRbFOe0Ntt6qkJtkOIrIdiz+RCkmShDGpWTh4/lpTp0JETcApT/uvXr0aISEhcHV1RXR0NPbt22dx/FdffYW7774brq6u6NWrF/7973+b9E+ePBmCIJj8DB8+3JGbQOTUblZWO6zwh+k90caN7yuInJnTHaEbN25EQkICUlNTER0djZUrVyI2NhYnT56Ev79/nfG7d+/GhAkTkJycjMceewwbNmzA6NGjkZ2djZ49exrHDR8+HOvWrTO+1ul0jbI9RM7uw6e6IrxTEATh9/cCkiSitKgQ3j6+xnZr29q4tYLey7XxN4SIrOZ0xT8lJQXTpk1DXFwcACA1NRVbtmxBWloa5syZU2f8qlWrMHz4cMyePRsAsHDhQmzbtg0ffPABUlNTjeN0Oh0CAwOtyqGiogIVFRXG16WlpQBuX98URVH2ttlKFEVIkmTXOuXEsHYZS+Pk9Jlrr92mxD6RQ6n12hpHibkw13/nOH8PLUJ93etc8y8Qy6C/o93attrxHcGZjw1b58KaPh4b8sc3tL9rSJJocf9aare2rTFYuz6nKv4GgwEHDx5EYmKisU2j0SAmJgZZWVlml8nKykJCQoJJW2xsLDZv3mzSlpmZCX9/f7Rr1w4PPfQQFi1aBF9fX7Mxk5OTMX/+/DrtBQUFMBgMNm6VfKIooqSkBJIkmfxidXQMa5exNE5On7n22m1K7BM5lFqvrXGUmAtz/Tcrq4195WWlyM/Pt3ku5GyPUpz52LB1Lqzp47Ehf7ylcXceB6XXipCPWw0u5+zHRklJiVXjnKr4X716FdXV1QgICDBpDwgIwIkTJ8wuk5uba3Z8bm6u8fXw4cPx5JNPIjQ0FGfOnMEbb7yBESNGICsrCy4uLnViJiYmmvxBUVpaiuDgYOj1erRt29aOLbSNKIoQBAF6vd6uX3C2xrB2GUvj5PSZa6/dpsQ+kUOp9doaR4m5MNdfbqgy9rl7eMPf39/muZCzPUpx5mPD1rmwpo/HhvzxlsbdeRx4t/OBv793g8s5+7Gh1WqtGudUxd9Rxo8fb/x3r1690Lt3b3Tt2hWZmZl4+OGH64zX6XRm7wnQaDSNOokAIAiC3euVE8PaZSyNk9Nnrr12mxL7RA6l1mtrHCXmona/pf1rqd3atsbgzMeGLXNhbR+PDfnj6xtnui+t2+f1tTvLsWHtupzqbn8/Pz+4uLggLy/PpD0vL6/e6/WBgYE2jQeALl26wM/PD6dPn7Y/aSIiombGqYq/VqtF3759kZGRYWwTRREZGRkYMGCA2WUGDBhgMh4Atm3bVu94ALh48SIKCwvRvn17ZRInIiJqRpyq+ANAQkIC1q5di/Xr1+PXX3/FjBkzUFZWZrz7f+LEiSY3BL766qtIT0/Hu+++ixMnTuDtt9/GgQMHEB8fDwC4ceMGZs+ejT179uDcuXPIyMjAqFGjEBYWhtjY2CbZRiIioqbkdNf8x40bh4KCAsybNw+5ubmIjIxEenq68aa+nJwck2sa999/PzZs2IC5c+fijTfeQHh4ODZv3mz8jL+LiwuOHDmC9evXo7i4GEFBQRg2bBgWLlzIz/oTEZEqOV3xB4D4+HjjO/faMjMz67SNHTsWY8eONTvezc0NW7duVTI9IiKiZs3pTvsTERGRY7H4ExERqQyLPxERkcqw+BMREakMiz8REZHKsPgTERGpDIs/ERGRyrD4ExERqYxTPuSHiIioMd2qrDb5il9RFHHz/7fV/krf2u31tUmS1LgbYQMWfyIiUr3HP9iteMwege747mV/xeMqgaf9iYhIldxauyCiYxuHxT+WW45L1246LL49+M6fiIhUSRAEfDNjAI6euQj3Nu0gCL+/H5YkEaVFhfD28W2wvXbbrcpqPPb+LgBARbXYuBtlJRZ/IiJSLUEQEOCtg7/es861/Xyp3Kr22m133jvgrHjan4iISGX4zp/IyUmSVOeO49rquzO5vv5yQ7UjUyYiJ8fiT+TEJEnC9E0n8cuVsqZOhYhaEJ72J3JiNyurHVr4w/Qe8NS5OCw+ETknvvMnaibWTY5CsI+H2b767kxuqN9bp4F4s9RhORORc2LxJ2omgn08EObvabavvjuTG+oXRRH5LP5EqsPT/kRERCrD4k9ERKQyLP5EREQqw+JPRESkMiz+REREKsPiT0REpDIs/kRERCrD4k9ERKQyfMgPkYIkScLNyvq/NKehL+CpjV/AQ0SOwOJPpBBJkjAmNQsHz19r6lSIiCziaX8ihdysrHZY4Q/Te6CNG/9WJyJl8LcJkQPU9yU8DX0BT33jg4MCoPdydUSqRKRCLP5EDlDfl/A09AU89Y33Y+EnIgXxtD8REZHK8J0/qVJDd+XXZs1d+rwzn4iaCxZ/Uh3elU9EasfT/qQ6jrwrHwDC9J68M5+InJpT/oZavXo1li9fjtzcXEREROD9999H//796x3/1Vdf4a233sK5c+cQHh6OpUuXYuTIkcZ+SZKQlJSEtWvXori4GA888ADWrFmD8PDwxtgcskPN6fma0+22PiTHnDtPz9d3V37dPKy/S7+NWyvemU9ETs3piv/GjRuRkJCA1NRUREdHY+XKlYiNjcXJkyfh7+9fZ/zu3bsxYcIEJCcn47HHHsOGDRswevRoZGdno2fPngCAZcuW4b333sP69esRGhqKt956C7GxsTh+/DhcXa3/JV1uqILWUKXYtjZEiUInJ4a1y1gaJ6evdrskAWNTd+P4les2bLFt6rsrvzZb79InInJmTlf8U1JSMG3aNMTFxQEAUlNTsWXLFqSlpWHOnDl1xq9atQrDhw/H7NmzAQALFy7Etm3b8MEHHyA1NRWSJGHlypWYO3cuRo0aBQD49NNPERAQgM2bN2P8+PF1YlZUVKCiosL4urS0FABwX/J2aHTuim8zNY0wvQe8dbfPJjREFEVIkmTVWCXjWDu+oXH19dvSbm1bY1BivXJiWLOM3Lmw1GfNvm/OcyEnjjMfG3f2uWqERp0Ta9flVMXfYDDg4MGDSExMNLZpNBrExMQgKyvL7DJZWVlISEgwaYuNjcXmzZsBAGfPnkVubi5iYmKM/W3atEF0dDSysrLMFv/k5GTMnz9fgS0iJXT0boV3RobCQ9caEiSU3yiFu6c3BAh2xfXUuUC8WYr8m6UNjhVFESUlJZAkya53/rbGsXZ8Q+Pq67el3dq2xqDEeuXEsGYZuXNhqc+afd+c50JOHGc+NiRJQsaM3rhaVAyNoRT5+WW27Aq7lJSUWDXOqYr/1atXUV1djYCAAJP2gIAAnDhxwuwyubm5Zsfn5uYa+2va6htTW2JioskfFKWlpQgODsbmP0XDy7utTdtkj9vXmYvg7eNj1dPglIph7TKWxsnpM9cuSSLEshKEde5g/AVXUFAAvV7f6L/gBEGwe722xrF2fEPj6uu3pd3atsagxHrlxLBmGblzYanPmn3fnOdCTpzmcGy4a1s1+nxotVqrxjlV8XcWOp0OOp2uTnsXvTfatvVutDxEUUQ+bsHf39uuX3C2xrB2GUvj5PSZaxdFEfn5t6DRaIxtgiCYvG4sSq3X1jjWjm9oXH39trRb29YYlFivnBjWLCN3Liz1WbPvm/NcyInDY6Mua9flVHcu+fn5wcXFBXl5eSbteXl5CAwMNLtMYGCgxfE1/7UlJhERUUvmVMVfq9Wib9++yMjIMLaJooiMjAwMGDDA7DIDBgwwGQ8A27ZtM44PDQ1FYGCgyZjS0lLs3bu33phEREQtmdOd9k9ISMCkSZMQFRWF/v37Y+XKlSgrKzPe/T9x4kR06NABycnJAIBXX30VgwcPxrvvvotHH30UX375JQ4cOICPPvoIwO3TLjNnzsSiRYsQHh5u/KhfUFAQRo8e3VSbSURE1GScrviPGzcOBQUFmDdvHnJzcxEZGYn09HTjDXs5OTkm1zTuv/9+bNiwAXPnzsUbb7yB8PBwbN682fgZfwB47bXXUFZWhunTp6O4uBgPPvgg0tPTbfqMPxERUUvhdMUfAOLj4xEfH2+2LzMzs07b2LFjMXbs2HrjCYKABQsWYMGCBUqlSERE1Gw51TV/IiIicjwWfyIiIpVh8SciIlIZFn8iIiKVYfEnIiJSGRZ/IiIilWHxJyIiUhmn/Jy/s5EkCcDtxwI39rdlXb9+Ha6urnZ9sY+tMaxdxtI4OX3m2mu3KbFP5FBqvbbGUWIuLPXbMxdytkcpznxsyJ0LS308NuSPV9uxUVpaCuD3ulUfFn8rFBYWAgA6d+7cxJkQERE17Pr162jTpk29/Sz+VvDx8QFw+9HClnamI/Tr1w/79+9v9BjWLmNpnJw+c+13tpWWliI4OBgXLlyAt3fjfb1yfbk1Rhwl5sJSv9y5AJr/fDjq2JA7F5b6eGzIH6+mY0OSJPTt2xdBQUEWx7H4W6HmlE2bNm0a/aBycXGxe51yYli7jKVxcvrMtZtr8/b2bpZzISeOEnNhqd/euQCa73w46tiQOxeW+nhsyB+vtmNDq9U2eKmBN/w5uZdeeqlJYli7jKVxcvrMtSuxD5SgVB62xlFiLiz1N8e5AJz72JA7F5b6nHk+eGw4z1wA1uUiSA3dFUAoLS1FmzZtUFJS0uh/wZEpzoVz4Xw4D86Fc3H2+eA7fyvodDokJSVBp9M1dSqqx7lwLpwP58G5cC7OPh98509ERKQyfOdPRESkMiz+REREKsPiT0REpDIs/kRERCrD4k9ERKQyLP52KC4uRlRUFCIjI9GzZ0+sXbu2qVMiAOXl5ejcuTNmzZrV1KmoWkhICHr37o3IyEgMHTq0qdNRvbNnz2Lo0KHo3r07evXqhbKysqZOSZVOnjyJyMhI44+bmxs2b97c6Hnwo352qK6uRkVFBdzd3VFWVoaePXviwIED8PX1berUVO3NN9/E6dOnERwcjBUrVjR1OqoVEhKCo0ePwtPTs6lTIQCDBw/GokWLMHDgQBQVFcHb2xutWvEJ703pxo0bCAkJwfnz5+Hh4dGo6+Y7fzu4uLjA3d0dAFBRUQFJkhr8GkVyrFOnTuHEiRMYMWJEU6dC5DSOHTuG1q1bY+DAgQBuf1kZC3/T++677/Dwww83euEHVF78d+zYgccffxxBQUEQBMHsqZfVq1cjJCQErq6uiI6Oxr59+0z6i4uLERERgY4dO2L27Nnw8/NrpOxbHiXmY9asWUhOTm6kjFsuJeZCEAQMHjwY/fr1w+eff95ImbdM9s7HqVOn4Onpiccffxz33nsv3nnnnUbMvmVR4tiosWnTJowbN87BGZun6uJfVlaGiIgIrF692mz/xo0bkZCQgKSkJGRnZyMiIgKxsbHIz883jmnbti0OHz6Ms2fPYsOGDcjLy2us9Fsce+fjX//6F+666y7cddddjZl2i6TEsbFr1y4cPHgQ3333Hd555x0cOXKksdJvceydj6qqKuzcuRMffvghsrKysG3bNmzbtq0xN6HFUOLYAG4/+3/37t0YOXJkY6Rdl0SSJEkSAOnbb781aevfv7/00ksvGV9XV1dLQUFBUnJystkYM2bMkL766itHpqkacuZjzpw5UseOHaXOnTtLvr6+kre3tzR//vzGTLtFUuLYmDVrlrRu3ToHZqkecuZj9+7d0rBhw4z9y5Ytk5YtW9Yo+bZk9hwbn376qfTMM880RppmqfqdvyUGgwEHDx5ETEyMsU2j0SAmJgZZWVkAgLy8PFy/fh0AUFJSgh07dqBbt25Nkm9LZ818JCcn48KFCzh37hxWrFiBadOmYd68eU2VcotlzVyUlZUZj40bN27ghx9+QI8ePZok35bOmvno168f8vPzce3aNYiiiB07duCee+5pqpRbLGvmokZTnvIHAN7xUY+rV6+iuroaAQEBJu0BAQE4ceIEAOD8+fOYPn268Ua/l19+Gb169WqKdFs8a+aDGoc1c5GXl4cnnngCwO1PxUybNg39+vVr9FzVwJr5aNWqFd555x0MGjQIkiRh2LBheOyxx5oi3RbN2t9TJSUl2LdvH77++uvGTtGIxd8O/fv3x6FDh5o6DTJj8uTJTZ2CqnXp0gWHDx9u6jToDiNGjOCnYJxEmzZtmvz+MJ72r4efnx9cXFzqTFBeXh4CAwObKCv14nw4D86Fc+F8OI/mNBcs/vXQarXo27cvMjIyjG2iKCIjIwMDBgxowszUifPhPDgXzoXz4Tya01yo+rT/jRs3cPr0aePrs2fP4tChQ/Dx8UGnTp2QkJCASZMmISoqCv3798fKlStRVlaGuLi4Jsy65eJ8OA/OhXPhfDiPFjMXTfY5Ayewfft2CUCdn0mTJhnHvP/++1KnTp0krVYr9e/fX9qzZ0/TJdzCcT6cB+fCuXA+nEdLmQs+25+IiEhleM2fiIhIZVj8iYiIVIbFn4iISGVY/ImIiFSGxZ+IiEhlWPyJiIhUhsWfiIhIZVj8iYiIVIbFn4iISGVY/Imo2SksLIS/vz/OnTunaNzjx4+jY8eOKCsrUzQukbNh8SdqwSZPngxBEOr8DB8+vKlTs8vixYsxatQohISEWDX+8ccfr3ebd+7cCUEQcOTIEXTv3h333XcfUlJSFMyWyPnw2f5ELdjkyZORl5eHdevWmbTrdDq0a9fOYes1GAzQarUOiV1eXo727dtj69atuO+++6xaZvPmzXjqqadw/vx5dOzY0aRvypQp+OWXX7B//34AwJYtWzBt2jTk5OSgVStVf/EptWB850/Uwul0OgQGBpr83Fn4BUHAxx9/jCeeeALu7u4IDw/Hd999ZxLj6NGjGDFiBDw9PREQEIDnnnsOV69eNfYPGTIE8fHxmDlzJvz8/BAbGwsA+O677xAeHg5XV1cMHToU69evhyAIKC4uRllZGby9vfHPf/7TZF2bN2+Gh4cHrl+/bnZ7/v3vf0On09Up/JZyfOyxx6DX6/HJJ5+YLHPjxg189dVXmDp1qrHtkUceQVFREX788Ucr9zBR88PiT0SYP38+nn76aRw5cgQjR47EM888g6KiIgBAcXExHnroIfTp0wcHDhxAeno68vLy8PTTT5vEWL9+PbRaLX766Sekpqbi7NmzGDNmDEaPHo3Dhw/jhRdewJtvvmkc7+HhgfHjx9c5K7Fu3TqMGTMGXl5eZnPduXMn+vbta9LWUI6tWrXCxIkT8cknn+DOk51fffUVqqurMWHCBGObVqtFZGQkdu7cKWNPEjUTTfqFwkTkUJMmTZJcXFwkDw8Pk5/FixcbxwCQ5s6da3x948YNCYD0n//8R5IkSVq4cKE0bNgwk7gXLlyQAEgnT56UJEmSBg8eLPXp08dkzOuvvy717NnTpO3NN9+UAEjXrl2TJEmS9u7dK7m4uEiXL1+WJEmS8vLypFatWkmZmZn1btOoUaOkKVOmmLRZk+Ovv/4qAZC2b99uHDNw4EDp2WefrbOOJ554Qpo8eXK9ORA1d7ygRdTCDR06FGvWrDFp8/HxMXndu3dv4789PDzg7e2N/Px8AMDhw4exfft2eHp61ol95swZ3HXXXQBQ5934yZMn0a9fP5O2/v3713ndo0cPrF+/HnPmzME//vEPdO7cGYMGDap3e27evAlXV1eTNmtyvPvuu3H//fcjLS0NQ4YMwenTp7Fz504sWLCgzjJubm4oLy+vNwei5o7Fn6iF8/DwQFhYmMUxrVu3NnktCAJEUQRw+7r4448/jqVLl9ZZrn379ibrkeP555/H6tWrMWfOHKxbtw5xcXEQBKHe8X5+frh27ZpJm7U5Tp06FS+//DJWr16NdevWoWvXrhg8eHCdZYqKitC1a1dZ20PUHPCaPxFZdO+99+LYsWMICQlBWFiYyY+lgt+tWzccOHDApK3mjvo7Pfvsszh//jzee+89HD9+HJMmTbKYT58+fXD8+HFZOT799NPQaDTYsGEDPv30U0yZMsXsHxpHjx5Fnz59LOZB1Jyx+BO1cBUVFcjNzTX5ufNO/Ya89NJLKCoqwoQJE7B//36cOXMGW7duRVxcHKqrq+td7oUXXsCJEyfw+uuv47fffsOmTZuMd9vfWXDbtWuHJ598ErNnz8awYcPqfBSvttjYWBw7dszk3b+1OXp6emLcuHFITEzElStXMHny5Drxz507h0uXLiEmJsbKPUTU/LD4E7Vw6enpaN++vcnPgw8+aPXyQUFB+Omnn1BdXY1hw4ahV69emDlzJtq2bQuNpv5fIaGhofjnP/+Jb775Br1798aaNWuMd/vrdDqTsVOnToXBYMCUKVMazKdXr1649957sWnTJlk5Tp06FdeuXUNsbCyCgoLqxP/iiy8wbNgwdO7cucFciJorPuSHiBrN4sWLkZqaigsXLpi0f/bZZ/jzn/+My5cvW/VwoC1btmD27Nk4evSoxT9AbGUwGBAeHo4NGzbggQceUCwukbPhDX9E5DAffvgh+vXrB19fX/z0009Yvnw54uPjjf3l5eW4cuUKlixZghdeeMHqpwI++uijOHXqFC5duoTg4GDF8s3JycEbb7zBwk8tHt/5E5HD/PnPf8bGjRtRVFSETp064bnnnkNiYqLxsblvv/02Fi9ejEGDBuFf//qX2Y/qEZHyWPyJiIhUhjf8ERERqQyLPxERkcqw+BMREakMiz8REZHKsPgTERGpDIs/ERGRyrD4ExERqQyLPxERkcr8P0Ac+y4W0JLkAAAAAElFTkSuQmCC",
      "text/plain": [
       "<Figure size 500x400 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, ax = plt.subplots(figsize=(5, 4), layout=\"constrained\")\n",
    "\n",
    "hmf001.flux.plot_spectrum(ax=ax)\n",
    "\n",
    "ax.set(xlim=(1e3, 2e7))\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cell-eindex-intro-md",
   "metadata": {},
   "source": [
    "### 5. E-Index Similarity via AssimilationSuite\n",
    "The E-index measures the cosine similarity between two sensitivity vectors:\n",
    "\n",
    "$$\n",
    "E_{ij} = \\frac{S_i^T S_j}{\\sqrt{S_i^T S_i \\cdot S_j^T S_j}}\n",
    "$$\n",
    "\n",
    "Unlike $c_k$, the E-index does not involve a covariance matrix, so it captures purely geometric similarity of the sensitivity profiles. We construct a minimal `AssimilationSuite` with HMF001 as both the benchmark and the application to illustrate the API."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "cell-build-suite",
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Reading ../data/hmf001.ser_res.m\n",
      "WARNING:  SERPENT Serpent 2.2.1 found in ../data/hmf001.ser_res.m, but version 2.1.31 is defined in settings\n",
      "WARNING:    Attempting to read anyway. Please report strange behaviors/failures to developers.\n",
      "  - done\n",
      "Reading ../data/hmf001.ser_sens0.m\n",
      "  - done\n",
      "Reading ../data/hmf002-001.ser_res.m\n",
      "WARNING:  SERPENT Serpent 2.2.1 found in ../data/hmf002-001.ser_res.m, but version 2.1.31 is defined in settings\n",
      "WARNING:    Attempting to read anyway. Please report strange behaviors/failures to developers.\n",
      "  - done\n",
      "Reading ../data/hmf002-001.ser_sens0.m\n",
      "  - done\n",
      "Reading ../data/hmf002-002.ser_res.m\n",
      "WARNING:  SERPENT Serpent 2.2.1 found in ../data/hmf002-002.ser_res.m, but version 2.1.31 is defined in settings\n",
      "WARNING:    Attempting to read anyway. Please report strange behaviors/failures to developers.\n",
      "  - done\n",
      "Reading ../data/hmf002-002.ser_sens0.m\n",
      "  - done\n"
     ]
    }
   ],
   "source": [
    "suite = andalus.AssimilationSuite.from_yaml(\"config.yaml\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "cell-eindex",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "              HMF001  HMF002-001  HMF002-002\n",
      "HMF001      1.000000    0.960530    0.957067\n",
      "HMF002-001  0.960530    1.000000    0.999763\n",
      "HMF002-002  0.957067    0.999763    1.000000\n"
     ]
    }
   ],
   "source": [
    "e_index = suite.e_index_matrix()\n",
    "print(e_index)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cell-eindex-interp-md",
   "metadata": {},
   "source": [
    "The diagonal entries equal 1.0 by construction. The off-diagonal value between the benchmark and application reflects how closely aligned their sensitivity profiles are in the multivariate nuclear-data space. An E-index near 1.0 indicates that the two systems are driven by the same energy and reaction dependencies, making the benchmark a strong predictor of the application's behavior."
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "andalus",
   "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.12.9"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}