{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "d2e9adfb-9962-40f7-8265-4861621e757e",
   "metadata": {},
   "source": [
    "# Adjusting correlation parameters for entropy scaling of viscosity"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4b54ad7f-d2ee-478c-8dd3-aebb2addec81",
   "metadata": {},
   "source": [
    "## Goal\n",
    "\n",
    "- Read in experimental data for viscosity of a pure substance.\n",
    "- Use the `DataSet`, `Loss` and `Estimator` objects to store and work with experimental data.\n",
    "- Define a `cost` function that will be used with `scipy.optimize.least_squares`.\n",
    "- Run the optimization."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "805bdfdb-d397-45a3-8354-082a8d13e4c1",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import pandas as pd\n",
    "from scipy.optimize import least_squares\n",
    "import matplotlib.pyplot as plt\n",
    "import seaborn as sns\n",
    "\n",
    "from feos.si import BAR, KELVIN, KILOGRAM, METER, MILLI, PASCAL, SECOND, RGAS\n",
    "from feos.eos import EquationOfState, State, Contributions\n",
    "from feos.pcsaft import Identifier, PcSaftParameters, PcSaftRecord, PureRecord\n",
    "from feos.eos.estimator import Estimator, Loss, DataSet, Phase\n",
    "\n",
    "sns.set_context('talk')\n",
    "sns.set_palette('Dark2')\n",
    "sns.set_style('ticks')\n",
    "colors = sns.palettes.color_palette('Dark2', 5)\n",
    "plt.rcParams['figure.figsize'] = [12,7]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "46c5a891-be8b-44d9-8f42-c84adf0cadb4",
   "metadata": {},
   "source": [
    "# `DataSet` objects\n",
    "\n",
    "FeO$_\\text{s}$ provides a range of data structures that facilitate the adjustment of parameters. One such structure is the `DataSet`, which serves as a storage unit for experimental data. When working with a specific model, a `DataSet` enables the evaluation of a `cost` function. How this cost function is defined depends on the property that you want to predict with your equation of state. In this notebook, we use the viscosity where the cost function is the relative difference between the experimental data and the model prediction.\n",
    "\n",
    "A DataSet encompasses the following components:\n",
    "\n",
    "- The `target`: This refers to the experimental data associated with the property for which we intend to adjust parameters. In our case, this is the viscosity.\n",
    "- Other properties: These are necessary for determining the thermodynamic conditions. For instance, temperature is required for vapor pressure calculations, while liquid density computations necessitate temperature and pressure. For the viscosity, we use pressure, temperature and the phase (liquid or vapor) as input.\n",
    "- Each property available in FeO$_\\text{s}$ is accompanied by a corresponding DataSet constructor."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "3c80cd76-4868-4b92-a654-287c7d65b875",
   "metadata": {},
   "outputs": [],
   "source": [
    "# read data from file into DataFrame\n",
    "data = pd.read_csv(\"data/hexane_viscosity.csv\")\n",
    "# map strings for phase to Phase objects\n",
    "phase = [Phase.Vapor if p == 'vapor' else Phase.Liquid for p in data.phase]\n",
    "\n",
    "# construct a `DataSet` for vapor pressure\n",
    "dataset = DataSet.viscosity(\n",
    "    target=data[\"viscosity / mPas\"].values * MILLI * PASCAL * SECOND,\n",
    "    temperature=data[\"temperature / K\"].values * KELVIN,\n",
    "    pressure=data[\"pressure / bar\"].values * BAR,\n",
    "    phase=phase\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "30f51ce6-a910-4e6f-93d2-7c02d070143e",
   "metadata": {},
   "source": [
    "For `DataSet.viscosity` we have to supply the experimental data for viscosity, temperature, and pressure.\n",
    "\n",
    "Optionally, a list of `Phase` objects can be supplied (via the `phase` argument) in which case not the stable phase is calculated for given temperature and pressure but the possibly instable phase according to the input phase."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "65c565d7-30a4-446a-8e22-0c225edeb88f",
   "metadata": {},
   "source": [
    "# `Estimator` and `Loss` objects"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "89312d88-4eba-4649-af06-a11d673a88b9",
   "metadata": {},
   "source": [
    "To collect the `DataSet` objects for the properties we wish to adjust parameters for, we can assemble them into an `Estimator` object. The `Estimator` requires the following inputs:\n",
    "\n",
    "- A list of `DataSet` objects.\n",
    "- A corresponding list of `weights`, with each weight assigned to a specific `DataSet`.\n",
    "- A list of `losses`, where each `DataSet` has an associated loss function.\n",
    "\n",
    "In our case, there is only a single `DataSet` so we can ignore the `weights`. Since our data comes from a correlation (from the NIST WebBook), we don't need a loss function."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "2e4c20fb-2415-43ac-87c4-68279dcb4e7d",
   "metadata": {},
   "outputs": [],
   "source": [
    "estimator = Estimator(data=[dataset],weights=[1],losses=[Loss.linear()])"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c49e1f8a-93f5-4981-b7ea-410b651058c8",
   "metadata": {},
   "source": [
    "# Defining a cost function\n",
    "\n",
    "When using `scipy`'s `least_squares` solver, it requires a function to calculate the cost or residuals. The first argument of this function must be the parameter vector that's being optimized. Therefore, it is necessary to define this function in Python before working with `least_squares`.\n",
    "\n",
    "To simplify the process, we create two functions:\n",
    "\n",
    "- `eos_from_parameters`: This function constructs the parameters and equation of state based on the current parameter vector. Since FeO$_\\text{s}$ does not allow mutation of an existing `EquationOfState` object, generating a new equation of state is necessary. We can use this function later to generate the equation of state for the final parameters.\n",
    "- `cost`: Taking the parameters and estimator as inputs, this function builds the equation of state, calculates the cost of the current model, and returns the result.\n",
    "\n",
    "These functions aid in the seamless integration of the solver, allowing for the generation of the equation of state and the calculation of the model's cost."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "8002fc42-49dd-4166-b85e-7c96f567fa8d",
   "metadata": {},
   "outputs": [],
   "source": [
    "# extract PC-SAFT parameters, identifier and molarweight, which are kept constant during optimization\n",
    "hexane = PcSaftParameters.from_json([\"hexane\"], \"../../../parameters/pcsaft/esper2023.json\").pure_records[0]\n",
    "identifier = hexane.identifier\n",
    "saft = hexane.model_record\n",
    "mw = hexane.molarweight\n",
    "m, sigma, epsilon_k = saft.m, saft.sigma, saft.epsilon_k\n",
    "\n",
    "def eos_from_parameters(c):\n",
    "    \"\"\"Returns equation of state (PC-SAFT) for current parameters.\"\"\"\n",
    "    global m, sigma, epsilon_k, identifier, mw\n",
    "    model_record = PcSaftRecord(m, sigma, epsilon_k, viscosity=c)\n",
    "    pure_record = PureRecord(identifier, mw, model_record)\n",
    "    parameters = PcSaftParameters.from_records([pure_record])\n",
    "    return EquationOfState.pcsaft(parameters)\n",
    "\n",
    "def cost(p, estimator):\n",
    "    \"\"\"Calculates cost function for current parameters.\"\"\"\n",
    "    return estimator.cost(eos_from_parameters(p))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b3d9f5f9-77ff-4523-a3ec-d7a3c05b829a",
   "metadata": {},
   "source": [
    "# Adjust parameters\n",
    "\n",
    "To begin parameter adjustment, there are two additional requirements:\n",
    "\n",
    "1. An initial guess for the parameters.\n",
    "2. Bounds for the parameters.\n",
    "\n",
    "It is crucial to note that trying multiple combinations of initial parameters is recommended to avoid getting stuck in a local minimum.\n",
    "\n",
    "Once these prerequisites are fulfilled, we can invoke `least_squares` with `args=(estimator, )` as an argument for the `cost` function. Please note that the tuple syntax is necessary in this case."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "1309b82d-dd20-4e9c-987c-03167cf4a7b8",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "   Iteration     Total nfev        Cost      Cost reduction    Step norm     Optimality   \n",
      "       0              1         4.7163e-03                                    6.45e-03    \n",
      "       1              2         4.5509e-04      4.26e-03       1.00e+00       2.23e-03    \n",
      "       2              3         1.6840e-05      4.38e-04       1.75e+00       2.57e-04    \n",
      "       3              4         1.6086e-07      1.67e-05       1.21e+00       1.39e-05    \n",
      "       4              5         8.4360e-08      7.65e-08       2.34e-01       7.53e-08    \n",
      "       5              6         8.4357e-08      3.66e-12       4.59e-03       2.38e-11    \n",
      "`gtol` termination condition is satisfied.\n",
      "Function evaluations 6, initial cost 4.7163e-03, final cost 8.4357e-08, first-order optimality 2.38e-11.\n",
      "CPU times: user 227 ms, sys: 12.2 ms, total: 239 ms\n",
      "Wall time: 75.9 ms\n"
     ]
    }
   ],
   "source": [
    "%%time\n",
    "initial_parameters = [0.0]*4 \n",
    "fitted_parameters = least_squares(cost, initial_parameters, args=(estimator,), verbose=2).x"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c7810262-d26d-4a5a-8d3b-5b1d217a9eb5",
   "metadata": {},
   "source": [
    "To provide statistics during the adjustment process, we can utilize the `verbose=2` option, which prints relevant information. However, our primary interest lies in obtaining the optimal parameters. These parameters can be accessed from the `x` field of the `OptimizeResult` object generated by the `least_squares` function.\n",
    "\n",
    "Once the adjustment is complete, we can investigate the results. A quick and straightforward method for analyzing the optimization results is to use the `Estimator.mean_absolute_relative_difference` method. This method returns the mean absolute relative difference (MARD) for each `DataSet`, without applying losses or weights. It is important to exercise caution when dealing with `NaN` values, as any predictions resulting in `NaN` are filtered out and thus not immediately visible in the reported values.\n",
    "\n",
    "It's essential to note that the resulting MARD **does not include** an evaluation of the loss functions or weights. Therefore, it may not equal the property that is minimized in the optimization process."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "b266f557-40f5-43d6-82eb-f5c393a3b464",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Adjusted correlation parameters\n",
      "c    = [-1.20010418 -2.15742129  0.24509735  0.21828781]\n",
      "MARD = 0.2038 %\n"
     ]
    }
   ],
   "source": [
    "print(\"Adjusted correlation parameters\")\n",
    "print(\"c    = {}\".format(fitted_parameters))\n",
    "\n",
    "mard = estimator.mean_absolute_relative_difference(eos_from_parameters(fitted_parameters)) * 100\n",
    "print(\"MARD = {:<5.4} %\".format(mard[0]))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "342179f3-ca31-4ecf-8fa6-b58680ffd697",
   "metadata": {},
   "source": [
    "# Plot resuts\n",
    "\n",
    "We can now use our results and calcuate the phase diagram and have a visual inspection of our parameters. Here, we calculate the reduced logarithmic viscosity to use the entropy scaled depiction."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "4569e7a0-54ec-4fb7-b735-1dbb4dca268a",
   "metadata": {},
   "outputs": [],
   "source": [
    "fitted_eos = eos_from_parameters(fitted_parameters)\n",
    "results = []\n",
    "for _, (t, p, viscosity, phase) in data.iterrows():\n",
    "    s = State(fitted_eos, temperature=t*KELVIN, pressure=p*BAR, density_initialization=phase)\n",
    "    results.append({\n",
    "        's_res': s.molar_entropy(Contributions.Residual) / RGAS / m,\n",
    "        'ln_eta_nist': np.log(viscosity * MILLI * PASCAL * SECOND / s.viscosity_reference()),\n",
    "        'ln_eta_saft': s.ln_viscosity_reduced(),\n",
    "        'temperature': t,\n",
    "        'pressure': p,\n",
    "    })\n",
    "    \n",
    "# collect into DataFrame\n",
    "results_df = pd.DataFrame(results)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "90dca269-a213-4b2b-aff2-f006b4f522e2",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAbYAAAEhCAYAAAATaoWLAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjguMCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy81sbWrAAAACXBIWXMAAAsTAAALEwEAmpwYAABPaUlEQVR4nO3deZxN9f/A8dedfd/MGDOYGaRP9m2EoqSoLEVKyRctSvFL+0YoUWmRb6lvCyKiJJUlJSRtSFlTHwkzxmzMvq/398e5M11jxuxzZ3k/Hw+PM875nHPe98yd+76f5XyOyWw2I4QQQjQWdrYOQAghhKhJktiEEEI0KpLYhBBCNCqS2IQQQjQqktiEEEI0Kg62DqCpUUrlY3yhSLV1LEII0YB4AYVa63LzliS2umcHmDw9Pb1tHYgQQjQUaWlpUMFWRklsdS/V09PTe+/evbaOQwghGozw8HDS0tIq1NIlfWxCCCEaFUlsQgghGhVJbEIIIRoVSWxCCCEaFRk80kCk52STkJPB15F/oJPiUL6BXBfSiWbO7ng4u9g6PCGEqDcksTUA6TnZ7I4/wT3bV5JbWFC8fv5v3/DeoP/Q1S8YNwdnSXBCCIE0RTYICTkZxUnNZPWYodzCAu7dvpKsgnzis1OJz0ghJz/PhpEKIYTtSWJrAL6O/IPcwgIuSY1m44//5Y19H9Ex5TRgJLfNkX/g7eRGVkE+cVlpJGVl2DhiIYSwHUlsDYBOigPAuTAf18I8OqdGs2j/Kp45soHA7BT+ToonOSeLlNws3B2cSMnLIj4jhZTsTBtHLoQQdU/62BqAi30CATjgE8JD3W5jyj/fcXF6HIPO/MXlCceIIYOEkEvw8mxGck4WdiYThUB0ZjI5BXk42zvi7eJm2xchhBB1RGpsDcB1oR1xsrMH4KBPa+7vOZ6X1PWcdfLAuTCfsN0f4fHK1RTsWk1CZipO9vYk52QR5OZDVkE+CTnp0v8mhGgyJLE1AC72jrx71X+Kk5vZZGJLi87c3fdeUgbeB44uOKSfxXHN4zi/eQNZ+gdyC/JJyc3mVHoSzVw8yCrI53RGCnGS4IQQjZw0RTYAvs5uhHj4snXkQ2yJ/JO/k+No7xPIkJAOpOfmkH7ZODy/WYD5t89xjfmLwjdHQefriLruMcy+rUjOySYlN4tQTz+ScjJJSk3A19kVH2c3nB0cbf3yhBCiRkmNrQFwdnAkzKsZjnb2dPNvxci23enkF0RuQQF+Lu4EtexI1m0LyZi6FlNoDwBcDn+N3+vX4/HtQiLPRuJkZ09yTjYZebl4ObmQVZDHidQEGWAihGh0pMbWQDg7ONLa049AVy/yzYUUFBaQmpvN8dSzBLh64ufijumSgcS26ob7oa8o2DAXx7QzuO94B4e9a8kfMYOIjtfS3N2L1NxsknIyjRpcbiY5BXl4S+1N2MC6det4+umncXFxYevWrQQEBJyz/cYbb8TLy4sVK1YAoJRiwoQJzJgxo7hMVFQUb731Fr/++itxcXF4eXkRFhZGnz59mDZtWvE5ynPppZcWn0c0bJLYGhgnBwecgMzcHJq5uuPl5Epidga/nYmkm38rgjz9+O2iAbg9uImWu1fCd+/inH4WVj+Ma1BH8m5+EXNId/ItiTEpJ5Mwr2aczkjB18kVX1d3W79E0QRlZ2ezZMkSnnrqqUrtFxERwc0334yzszOjR4+mVatWxMfHc+TIEd5//32mTZtG7969efnll88516xZsxg8eDCDBw8uXu/v719jr0fYliS2BsrNydn4wZxLC3cvBjhdxOHEaI4kxnBdaCec7Oz5rd8deHQdTqttb2LevwHvmCPw5o1k9RyF3eCHyXVyJb+wgJScLAJcPUjKzSSvMF+m5xJ1rkOHDnz88cfcc889NGvWrML7LVu2jKysLL744gtatmx5zraEhAQAWrduTevWrYvXp6amMmvWLJRS3HjjjTXzAkS9In1sDZyLoxPODo642DvQJzCMWy7qRWRaIutPHKSrfys6tA3nrxueI+7uZRDcEQDX3z/Hf+H15G97E/JzybWMmMzMz6UQOJWRKLcHiDo1efJkcnNzWbJkSaX2i4yMJDg4+LykBlQqQYrGRRJbI+Hm5IyzgyOuDg70a9GGsRf3JiItkS+OH6CrfytCuw4l7YEvKLz5RczufphyMvD4+lWaL7qR7ENf809KPA4mO6LSkwly8yFTbg8QdSgsLIxhw4axevVqEhMTK7xfy5YtOX36NHv27KnF6ERDI4mtkXFzdCYnPx8wo3wCuTK4PfvPRlFogiBPPxwvv4P0x7ZSOOAuzHb2mM6cwOfDyQSunEpe7DESstNJyc0mKj0JbycXMvJzZfRkI5STn8fp9CQ+OPIzj/64lg+O/Mzp9CSbfom5//77yc7O5oMPPqjwPuPHj8fBwYEJEyYwatQoXnjhBbZt20Z2dnYtRirqO0lsjZCHswvezm54O7kQ5t2My1q0JT4zlV/jI7G3tycoIATHm+aR+eBGzO37A+B/YjdOrw3G85sFRJ6NwGw2FzdPejm5kJCTQVxGCuk58oHR0OXk57Er7gQDPnuVmbvX88nfe5m5ez0DPnuVXXEnbJbc2rZty9ChQ1m5ciXJyckV2qd9+/Z88cUXjBgxgqioKJYvX86UKVO47LLLWLt2be0GLOotSWyNmIujE64OTng6OtPFv+V5Ca55u95sG/YsZ257HfxaYyrIw//nZfi8di0uBzbwT3IcDiY7ojNSiM1Mlf63RuJsdjp3bl1+zrP9wHhSxJ1bl3M223ZPh5gyZUqla21t2rThlVdeYdeuXaxfv55HHnkEe3t7ZsyYwS+//FKL0Yr6qkmMilRKBQEPAn2AcMADuEprvaOC+3cAXgf6A7nABuBRrfXZWgm4hrk4OgFgLjTTxb8lLbN82RN/grjMNIa17YqpTVdOdRiE149LcfzuXZwzEnBe8zjmll3IG/MyCT4tcXNw4nR6MmGefqTn55KYmkCQm5dMrtwAbYn887ykViS3sIAtkUe4s+NldRyVoV27dlx//fWsXLmSu+66q1L72tvbo5RCKUX37t2ZMGEC69evp1+/frUUraivmkqNTQFPAq2Ag5XaUalWwE6gHTAdeBUYAWxRSjWoO5rdnJyxw4SXkzOXB7XjyuD27DsTSWZBHqF+LTnZbyIx09Zj7joUAJ/Th3BaOAzvDc9jykwmITud6MxUotKT8HdxJyEnQ2pvDdDhxOgLbj9SzvbaNmXKFDIzM1m+fHmVj9G5c2cA4uPjayos0YA0lcT2G+CvtW4PvFLJfacDrsBArfUbWusXgDFAD2B8zYZZ+5wdHPF0dsXDwZlQLz/6tWhHfFYaG08eoot/Ky5u25NtVz9K3B3vQ2B7TOZC/PZ+SsDCobj8+in/JMdiNpuJzkwlNjOVApDBJQ1MZ7/gC27vWM722nbRRRdx7bXX8uGHH5KWlnbBsnv37iUv7/wvVt9//z1gNFOKpqdJNEVqrS/813Fho4H1WuvTVsfbqpQ6ipHgllY3Plv49wZv6B0YSjuvgH+bJ9t0wRTWhVPt+uG5axXOW9/AlJFI0Mbncd/7Kfmjnucf/zb4u3gUN08m5GRQaDbLzCUNwJCQDjz/66ZSmyOd7OwZEtLBBlGda8qUKXz99dekpaWVeo9akffff58//viDwYMHo5QC4MiRI3zxxRf4+PgwceLEugpZ1CNNpcZWJUqplkBzYG8pm/dg1NoatAs2T/oE4n7NAyQ+vJm8niMB8Ir9C793xuD7xWxMGUnFzZOxmalkF+bLqMkGwN/Fgw+umVj8GKQiTnb2fHDNRPxdPGwU2b8uvvhihgwZUm65yZMnM2TIEPbu3curr77K3Llz+eGHHxg2bBhr1649Z8YR0XSYzGazrWOoU0qpkcDnVGDwiFIqHPgVuF1rvbrEtpeBxwEHrXWB1frkckLw9vT0ZO/e0nKlbeXm55NbmE+a5Tlu0RnJXOIbhJujExuPH6RL6mmab34JU/QRAMyu3sReNZXEniPxd/XC1cGJtt7++LlIra2+y8nP42x2Blsij3AkMZqOfsEMCemIv4u7TIYt6qXw8HDS0tJStNY+5ZVtEk2R1eBqWeaUsi3bqkx63YRTu4wJlh2ww4SHnwv+Lh4cSIgiLjONoW264GjXjfUeLejx93f4bn8LU1YKQV+9gOv+9eSPeZkE31Y42ztIYmsAnB0caenhY7PRj0LUJklsF5ZlWTqXss2lRBkAyvs2YanReVc3sNpU1P9mj4kgj04kZmUW979d364bBWFdiOs2HO8tC3D87TN8og9jXjSKgoH3cazP7XTxL7tPRAghapv0sV1YjGUZVMq2ICDeuhmysSnqf3NzdKS7f2t6BrRm/5lT5JkLaRV8Md9fMYW4ie9T6NsKU34OgVv/S+slE8k5/YetQxdCNGGS2C7AMhLyDMZN3SVdCuyv04BswNnBEV8XdwJcPOjgF0SXZq34MymGr04eZmibzgR3H07CtPXk9jPufPCK/YvI2eEkrH8Bs9zfJoSwAUlsVpRS7ZRS7Uqs/gy4wTJCsqjc1cDFwKd1GZ8teTi74GRnj5+L23m1t9b+rdnZ9w7i716OfUBbzPm5JKybSeTz/ciO2G/r0IUQTUyT6WNTSj1j+bHoJp3xSqn+QLLWepFl3TbLMsxq1xeAW4DvlFJvYkzH9ThwAPiwVoOuZ5wdHI1/dg4EuHoUDy6Jz0pjeNuu+Dn1xe3Sm0hYN4ukLQvJidhH5Jw++A17imY3zMDk4GTrlyCEaAKaTGIDni/x/6KJ6CKARZRBa31KKXUlsAB4CWOuyI3AI1rr3NoItL4rGlzi4eRCW5+A87YHjH0Vj96jiVsyidyYv0hcP5f037+gxd1LcGlTWquuEELUnCZ3H5utKaWSPT09vevjfWw1rTA3m4Qv55D01StgLgSTHb7XP0azkbOxc3Ip/wBCCGFRmfvYpI9N1Bo7JxcCbnmBkFm/4NSqM5gLSfrqZSJn9yLrmDxORAhROySxiVrn0iac0Gd/xe/GmWDvQG7MX5yaN4D41Y9SmCOTJwshapYkNlEnTA5O+I96lpBZu3EO7QFmM8nfLCRiZncy//re1uEJIRoRSWyiTrmEdidk5i80u+l5TA5O5MX/Q9RLg4hb8QCF2Y1iZjIhGpSnnnqKQYMGFf8/KioKpRTr1q2zYVTVI4lN1DmTgyPNbphOyHN7cWl7KQAp297m5IyuZPyx1cbRiYZu9erVDfpDWVSfJDZhM84tO9F6xg/4j5mPycGZ/IQITr9yLXEfTKYgM8XW4YkG6uOPP+bzzz+3dRgNVsuWLTl48CA33nijrUOpMklswqZM9g74DX2M0Of34dL+cgBSvl9MxDNdyTi42cbRicYuM1MGL5VkMplwdnbG3t6+/ML1lCQ2US84BSlaP/0dAeMWYnJyIz8xitMLhhP7/h0UZCTZOjxRi2JiYnjiiSfo168fnTt3ZsSIEWzcuPGcMuvWrUMpxf79+5k3bx59+/ale/fuTJ06lcTExOJygwYN4q+//mLPnj0opVBKMX78+HOOsXfvXmbNmkWfPn0YPnx48b4rV67k+uuvp3PnzlxxxRW89NJLZGWd8/AOBg0axJQpU/j++++54YYb6NKlCyNGjOD77/8dALVr1y6UUnz77bfnvdY1a9aglOLYsWNlXo+8vDwWLVrEkCFD6NKlC3369GHs2LH89NNP55Tbt28fkyZNIjw8nB49ejBy5Eg+/fTfWf727t3LtGnTGDhwIJ07d+bKK6/khRdeIDv7wg8DLq2P7amnniI8PJyYmBjuu+8+evToQd++fZk/fz4FBefOA5+UlMTjjz9Oz549CQ8P58knn+Svv/6q0367pjTziKjnTHb2+A5+APduQ4lbei9Zf+0g9acVZBz+lsCJb+PRs+E2jYjSxcfHM2bMGBwdHZkwYQLe3t5s27aNRx99lNzcXG666aZzyj/33HP4+PjwwAMPEBUVxfLly5kzZw4LFy4EYPr06cybNw8XFxfuu+8+APz9/c85xuzZswkICGDatGnk5RkTdb/55pssWrSI/v37M27cOI4ePcqyZcs4evQoS5YswWQyFe9//PhxHn/8ccaOHcuoUaP49NNPmTJlCitWrKBnz5706dOHoKAgNmzYwODBg88594YNG+jYsSMXXXRRmddk0aJFLFmyhNtvv5327duTlpbGoUOH+OOPP7j8cqNVY+fOnUyZMoXAwEDuuOMOmjVrhtaaHTt2cMsttwDw9ddfk52dzdixY/Hx8eHgwYOsXLmS2NhY3njjjUr/rvLz87nrrrvo2bMnTzzxBD///DNLly6ldevW3H777QAUFhZy//33c/DgQW6//XbatGnDtm3bePLJJyt9vuqQxCbqHafm7Wj1xLek7HiPM2uepCAllug3bsKzz600/88b2Hv6l3+QJsCcn0d+0mlbh4GDb0tMVXzq9sKFC7Gzs+OLL77Ay8sLgNtvv51JkyaxYMECRo4ciZ3dvw1Lfn5+LF68uDjRFBYWsmLFCtLS0vD09OSaa67hzTffxMvLq8w+Ij8/P5YuXVp83MTERN59912uvPJK3n333eJjt2rVitdee43vvvvunFGDJ06c4H//+1/xuptuuokhQ4bw+uuvs2LFCkwmEyNGjGD58uWkp6fj4eEBGDXTX3/9tdwP+aLkNH369FK3FxQU8Oyzz9KiRQs+//xzPD09i7dZzyT12GOP4eLy7ww/t956K6GhoSxYsIDo6GiCg4MvGEdJWVlZjBw5ksmTJwMUJ/a1a9cWJ7atW7eyb98+Zs2axbhx44rL3XnnnZU6V3VJYhP1ksnODp9B9+He9Xrilk0m8/C3pO3+hMwj22k+/g08et9yzrfopsacn8fJ6Z3Ii//H1qHg2LwdYS/8UenkZjab+fbbbxk+fDj5+fnnNCkOGDCAH374gRMnTtCu3b8P3LjtttvO+b2Hh4ezbNkyTp8+zSWXXFKh844ZM+acZPnzzz+Tl5fHxIkTzzn2uHHjeOONN9ixY8c5iS04OPic/3t7ezN8+HBWrVpFZmYmbm5ujBw5kvfee48tW7YU1zo3btyIyWRi2LBhF4zPy8uLAwcOEBsbS4sWLc7b/scff3D69Glmzpx5TlIDzonfOqllZmaSnZ1Njx49MJvNHDlypNKJDYzkaK1Xr16sX7+++P8//PADTk5O3HzzzcXr7OzsGDduHLt27ar0+apKEpuo1xz9Q2n56GZSf/iAM6sfoyDtDDFvj8Wj1xqaT3gLB+9AW4coqigxMZHU1FRWrVrFqlWrSi2TlHRu/2pQ0LnP/C2q5aWmplb4vK1atTrn/9HR0QC0adPmnPXu7u40b968eHuRkJCQ844ZGhpKYWEhMTExtGvXjnbt2tGpUyc2btxYnNg2bNhA3759ad68+QXjmzZtGlOmTGHgwIF07NiRAQMGMGLEiOLmy6ioKADat29/weNER0fzxhtvsH37dlJSzh1lnJ5e+XtG3dzc8PHxOWedt7f3OceOjo4mMDAQZ2fnc8qVds1qkyQ2Ue+ZTCa8r7gLty7XEr/sPjIOfEX6b5+T+df3NB/3Op79xjW52pvJwZGwF/5o0E2RhYWFgNGUN2LEiFLLlPzwLmukXmUmc7euydSmkSNH8tJLL3H27FkSExPRWvPSSy+Vu1/v3r359ttv2b59Oz/99BOrV69m8eLFzJkzh9GjR1fo3AUFBdx5552kpKQwadIk2rZti5ubG3FxcTz11FPF174yGtIoSUlsosFw9G1J8EPrSfvlI+I/eojCjERi35tI2u41NL/jfzj6tiz/II2IycERx4AwW4dRZX5+fri7u2M2m7nssstq7LiV/ZJT1CR34sSJc5rnMjMziY+Pp3///ueUj4yMPO8YERER2NnZnVOjHD58OPPnz2fTpk2cPXsWFxeX8waTlMXHx4ebbrqJm266iczMTMaPH88bb7zB6NGjad26NQB///03ffr0KXX/o0ePcvLkSebPn8/IkSOL15ccWVnTgoOD2bNnDzk5OefU2kq7ZrVJhvuLBsVkMuF12X8Im3cYj16jAMg4sImI6V1I2bm0Ut/chW3Z29szePBgvvrqK44fP37edus+t8pwdXWtVNPkZZddhqOjIytWrDjn/bNq1Sry8vIYOHDgOeWjo6PZvn178f9TUlLYuHEj4eHhuLm5Fa/38/NjwIABbNiwgU2bNjFo0KDigSQXUrL51c3NjbCwMHJycgDo2LEjLVu2ZNmyZaSlpZ1Ttij+oj5E69djNpv58MPafTZy//79yc3NZe3atcXrCgsL+eijj2r1vCVJjU00SA4+LQj6v09J//VT4lc8QEHaWeKW3kPanjUE3vkejs3qtk1fVM2jjz7K7t27GT16NLfeeitt27YlKSmJQ4cOceTIkXMSSEV16tSJlStX8vbbbxMaGoqfnx/9+vUrs7yfnx+TJ09m0aJF3HvvvQwcOJCjR4+yZs0aLr/8cq666qpzyrdp04annnqKsWPH4uvry5o1a0hPT+fBBx8879g33HADDz/8MAAzZ86sUPzDhg2jd+/edO7cGR8fHw4fPsxXX31VPMrQ3t6e2bNnM2XKFEaOHMmoUaPw9/fn2LFjxMbGsmjRItq2bUtISAjz588nLi4ODw8Pvvnmm0ol/Kq45ppr6Nq1K/PmzePEiRO0adPmnD6+uuoykMQmGiyTyYTnpWNwveQqznz0oDFq8vC3nJzRhYAx8/EeeC8mO2mUqM+aN2/Op59+yqJFi9i8eTMJCQn4+PiglCo1UVTE/fffT1RUFIsXLyYjI4NLL730gokN4IEHHsDHx4ePPvqIF198EV9fXyZMmMCDDz543odx27Ztefrpp3n11Vc5efIkoaGhLFq0iPDw858Of/XVV+Pp6YmDgwMDBgyoUPzjx49n+/bt/Pzzz+Tm5hIcHMyDDz7I3XffXVzmyiuvZNmyZcX3vAGEhYUVJz9HR0feeecd5s6dy7vvvouzszODBw9m3LhxtTpVlr29Pe+++y7z5s1j3bp12NnZMXjwYKZOncrYsWPPG1RSW+QJ2nWsKT1Bu66l//YFcR9OpSAlFgDXDlcReOd7ODVva+PIRGMxaNAgLrnkEt5+++0Klc/Ly6N///4MGzaMWbNm1XJ09dfWrVuZOnUqq1atolevXlU6hjxBWzRJHr1GEjbvEF6XG1MoZf35HRHPdCPp2zcxV2EUmBDVtWXLFpKTkxv0hMKVVXLKroKCAlasWIGHhwedOnWqkxikKVI0KvYefrS4Zxkel44hftl95Ced5sxHD5G251Na3L0YpxYX2zpE0QQcOHAArTWLFi2ie/fudOvWzdYh1ZnnnnuOvLw8unfvTm5uLlu2bGHfvn088sgjdXarhdTYRKPk0W0oofMO4X3lJACy//6JiJk9SPzqVcyFBeXsLUT1rF69mmeffZaAgABeeOEFW4dTp/r27cs///zD66+/zoIFC0hLS2PWrFnFU3HVBeljq2PSx1b3Mv7YStzSe8lPiADApe2lBN69BOeWHW0cmRCioqSPTQgr7p2uIWzuAbyvngJA9vE9RM7uRcKGFzHn59k4OiFETZPEJpoEO1dPAse/SauntuPYvB3m/FwSPnuGyOf7kRN5wNbhCSFqkCQ20aS4XXIloc/vx+fah8BkIidiHxHPXcrZz5/FnJ9r6/CEEDVAEptocuyc3Wg+9jVaT9+JYwsFBfkkfvk8Ec/2JvuE9H0K0dBJYhNNlmv7ywid8zu+Q58Akx25UYeJfP4yznw6ncLc7PIPIISol5rEfWxKKWdgDjAe8AUOADO01tvK2e9ZYHYpm+K01uc/AVA0OHZOLgSMeRHP8JuIXTqJ3KjDJG2aT8bvXxJ492JcL7rwVExCiPqnqdTYlgEPAyuBB4FCYLNSqqKfWpMxkmLRv6m1EKOwIZe2vQmZvQe/G54BewdyY/7i1LwBnFn9GIU5mbYOTwhRCY2+xqaUuhS4DXhYa73Qsu5D4DAwH7iiAodZo7VOrq0YRf1g5+iM/03P4dFrFHFL7iYncj9J37xO+v4NBN71Pm6qIm8VIYStNYUa281AHrC4aIXWOhtYAvRXSgWVtaMVk1LKSynVtB7T3ES5hHYnZNYumt30PNg7khd3jKgXryJ+xTQKs9NtHZ4QohxNIbH1AP7SWpf8RNoDmIDuFThGJJACpCilliql/Go2RFHfmBwcaXbDdELn/IZzm94AJG97i5PPdCPjj602jk4IcSGNvikSCAJOl7I+xrIMLmVbkSTgTWAXkAsMwuhv66mU6qO1zim5g1IquZx4vMsLWNQfzi07EfLMjyR98zoJ62aTf/Ykp1+5Fu8rJ+F/68vYu8mvU4j6pikkNlfgvAQEZFttL5XW+r8lVq1VSh0G3gImAO/XSISiXjPZO+A39HE8etxA7JJJZB/7mZTvF5Nx6GsC73gH967X2zpEIYSVppDYsoDSHtvqYrW9Mt4BXgGuppTEVt4EnZYanXzNb4CcghStp+8geetbnF07nfzEKE4vGI7X5RMIuH0B9u6+tg5RCEHT6GOLwWiOLKloXXRlDqa1LsRo2pR+tibIZGeP75BphM49gKu6EoDUnz7k5PTOpP/+pY2jE0JA00hs+4FLlFIeJdb3sSwrNQOuUsoRaA2cqX5ooqFyat6OVk9upfmERZhcPChIiSX6jZuIeWccBWlnbR2eEE1aU0hsawFHYFLRCstMJHcCP2mtoy3rQpRSl1jvqJQKKOV4j2M0Y35TaxGLBsFkZ4fPoPsJm3sQt07XAJC262NOTu9M2p5PbRydEE1Xo+9j01rvVkp9CrxsuWftH2AiEArcYVX0Q+BKjFsAikQopT7GuJk7B7gKGA38CKyq/ehFQ+DoH0rLx74mdedSznz8GAVpZ4h5+zbS9qyh+fhFOHgH2jpEIZqUplBjA2ME438tyzcwanBDtdY/lbPfR0Bf4DlgAdAFeB4YorXOr71wRUNjMpnwvvJuQucdKh4lmb53HSendyb154+QJ9ULUXdM8gdXt5RSyZ6ent5798rjURors9lM2s8riV/1MIUZSQC4dx9O84lv4+jb0sbRCdEwhYeHk5aWllLeyHNoOjU2IeqMyWTC6/LxhM07jEevkQBk7N9IxPQupPzwgdTehKhlktiEqCUOPi0I+r+1BN2/CntPfwqzUohbMonTrw0lLyHS1uEJ0WhJYhOiFplMJjz73ErovMN4XjoGgMzDW4iY0ZXk797FXFho4wiFaHwksQlRBxy8AgiaspqgB9Zi7xVIYXYa8cunEPXKEHLjj9s6PCEaFUlsQtQhz16jCHvhMJ6X/QeArD+/I+KZbiR9+6bU3oSoIZLYhKhj9h5+BN27nOCH1uPg2xJzbiZnPnqIUy8OJDf2qK3DE6LBk8QmhI14dB9G6LxDeF1xNwDZf/9ExMweJG5+DXNhgY2jE6LhksQmhA3Zu3nT4q73aPnY1zg0C8Gcl83ZT57g1Nz+5Jw+YuvwhGiQJLEJUQ+4dx5M2NyDeA+6H4Ds43uInN2LhA0vYs7Ps3F0QjQsktiEqCfsXD0JnLCIVk9uwzGgLeb8XBI+e4bI5y8jJ7JSD6EQokmrlcSmlJIOAiGqyK3DQELn7sdnyINgMpET8TsRz13K2c+fw5yfa+vwhKj3aqvGZiq/iBCiLHbO7jS/fQGtp3+PY4uLoSCfxC/nEPHspWSf/M3W4QlRr9VWYqvQZHhKKUmAQlyAa/vLCZ3zO75DHweTHblRh4ic048zn06nMDfb1uEJUS/VSR+bUmq91c//s9q0rS7OL0RDZufkSsCYlwiZ+TNOLTtBYQFJm+YT+Ww4Wcd22To8Ieqduho84mn1s/VTqqXGJkQFubTtTcizv+I3YgbY2ZMb/Sen5vXnzOrHKMzJtHV4QtQbdZXYzBX4WQhRDjtHZ/xHzyFk9m6cQ7qD2UzSN68TMasHmXqnrcMTol6oq8TWTSn1lVJqc4mfu9bR+YVoVFxCexAyaxfNbpoD9o7kxR0j6sWriF8xjcLsdFuHJ4RNOdTReXrW0XmEaDJMDo40u2EGHj1vJHbJ3eSc2EvytrdIP7CJFne9h1vHq20dohA2UW5iU0q5AtcCTsAJ4KDWOqeS5xkCHNdaFw8WUUpdDbQBFlfyWEIIK86tOhPyzE8kffM6Cetmk3/2JFEvD8F74D34j5mPvZu3rUMUok5VpMb2KdAHI6m1BbyUUkeBfcD+oqXWOvECx5iote5vvUJrvU0p9SOS2ISoNpO9A35DH8ej+whil04i+9gvpOx4n4yDmwm8413cu15n6xCFqDMV6WO7Euistb5Ua+0PtAOeAjTQD3gPiC/nGGVNlyDTKAhRg5yCL6H19O8JGLsAk5Mr+YlRnF4wjNjFd1GQkWTr8ISoExWpsf0DFPdGa61PAaeAjUXrlFIe5RwjRykVorWOtNonFElsQtQ4k509vtc+iHv3YcQtvZcs/T2pPy4n4/AWAie+jUePG2wdohC1qiKJ7Q1gNvBEWQW01uUNw3oa2KCU+hIjKYYCI4AJFYxTCFFJToEX0erJraR89y5n1jxJQXIM0f8dhWff22g+7r/Ye/rbOkQhakVFmiLfAx5TSh1USr2olLpFKXVRZU6itd4PDAD+AnyBI8AArbVMWS5ELTLZ2eFz9f2EzT2IWydjlGTaro85OaMLab+utXF0QtSOitTYvIBuQA/LvyeBzkqpbK21z4V2LGVE5WdVGFEphKgmx4AwWj72Dak7l3Dm48cpSI0n5q1bSQu/iebjF+HgHWjrEIWoMeUmNq11JvCL5R8ASikHoGMFjl8TIyqFEDXAZDLhfeUk3LpcR/yy+8g4uJn0vevI/HMHzcctxLPf7ZhMMsudaPiqdIO21jofOFiBolcCF2mt4wCUUq0xan/dMUZU3g+EVTWOilJKOQNzgPEYTaEHgBnW99VdYN+WwOsY9+LZAduBh7XWJ2ovYiFqj6NfK4If3kDazyuJX/UwhRmJxL43gbQ9awic+D8cfINtHaIQ1VLbM4/UxIjKmrAMGA0sBI4BdwCblVJXaq1/KWsnS2zfYUziPA/IBx4GdiilumutZfy0aJBMJhNel4/HrdM1xH04lYzfvyRj/0ZO6h8IuP01vPrfIbU30WDV9lyRRSMqy1SBEZXVopS6FLgNeEJr/YTW+j1gEBAJzC9n9ynARcBQrfUrWuuimltLjAQnRIPm4BNE8AOfEXT/Kuw9/SnMSiFuySROvzaUvITI8g8gRD1U24mt2iMqa8DNQB5WM5xorbOBJUB/pVRQOfvu0lrvs9r3L4znyI2pnXCFqFsmkwnPPrcSOu8wnpcab+vMw1uImNGV5O/exWyWh3CIhqW2E5sXcDnwDuCPMaLysFIquZbPa60H8FcpNcM9GM+D617aTkopO4ynD+wtZfMe4GKllFsNximETTl4BRA0ZTVBD6zF3iuQwuw04pdPIerlIeTGH7d1eEJUWG31sZmg2iMqa0oQcLqU9TGWZVk95X6As1W5kvuaLMf+x3pDBZK2zEgr6jXPXqNwU1cQv+oR0n5eSdaf24l4phv+t7yIz9VTMNnV1dOuhKiaWnmHaq3LPK7WOl9rXZERlTXFFSjt3rlsq+1l7UcV9xWiQbP3aEbQvcsJfuhL7H2CMedmcuajB4l66SpyY4/aOjwhLqiunsdmS1kYNa+SXKy2l7Ufld23AjetJyO1NtFAeHQfjuu8/pz55HFSdy4l6+iPRMzsQbOb5uB77UOY7OxtHaIQ52kKbQoxGE2GJRWtiy5jv0SM2lpZ+5opvZlSiEbF3t2HFne9T8vHNuPQLARzXjZnP3mCU3MHkHP6iK3DE+I8TSGx7QcuKeV+uT6WZanzVWqtC4FDQHgpm/sAf1v6EIVoEtw7DyFs7kG8B90HQPbx3UTO7kXChhcx5+fZODoh/lXtpkilVD/g/4D2QDMsA0esmLXW7ap7nmpYCzwGTMK4QbtoJpI7gZ+01tGWdSGAm2U4v/W+LyqlehQN+VdKKYz74F6qs1cgRD1h5+pJ4IS38Ox9C3FL7yHvzHESPnuG9L3raDFpCc6tu9o6RCGqV2NTSk0AfsSY1cMF46bniBL/bHqXp9Z6N8aclS8rpeYrpe7FmBYrFOP2gyIfAn+W2P1t4DjwlVLqMaXUQ8C3GE2Qr9d27ELUV24dBhI6dz8+g6eByUROxO9EPNubs58/hzlfHrMobKu6NbYZGE/Svqao5lNPTQCetyx9Mea5HKq1/ulCO2mt05RSAzGS2EyMLwLfAQ9prRNqNWIh6jk7Z3eaj3sdz0tvJnbJJPJij5L45RzSf/ucFpOW4BLWy9YhiibKVJ1ZBZRS2cDjWus3ay6kxk0plezp6em9d29p930L0TAV5maR8MVzJG1+DcyFYGeP39DH8bthJnZOLuUfQIhyhIeHk5aWllLeyHOo/uCRKEofDi+EaELsnFwJGPMSrWf+hFNwRygsIHHjS0Q+G07WsV22Dk80MdVNbO8A45RScjOLEALXtpcS8txe/EbMADt7cqP/5NS8AZz5+HEKc2QQsagb1e1j+w1j4MgepdRbGA8ULShZSGu9s5rnEUI0EHaOzviPnoNH+CjiFt9NzqkDJH29gPR96wm8azFuaoCtQxSNXHUTm/WDOhdj3LRszWRZJzU6IZoYl9AehMzeTeKm+SSsn0te3DGiXhyIz9VT8b/lBexc6uJRjKIpqm5iu7NGohBCNEomB0ea3fgMHj1vJHbpJHJO7CV521tkHPyKwDvfw63jIFuHKBqhao2KFJUnoyJFU2UuyCfp6wUkfP4s5nxjbnHvgffgf+vL2Lt62Tg6Ud/V5ahIIYSoEJO9A37DniB0zu+4tOsLQMqO94mY0YWMg1/bODrRmFSqKdIy00ilaa0/rMp+QojGxyn4ElrP2Enyt29y9rNnyE+M4vSCYXj1n0jA2Newd/e1dYiigatsH9syjMEgJeeDvBAzxnRVQggBgMnOHt9rH8K9+3Dilt5Dlt5J6o/LyTi8hcCJb+PR4wZbhygasMomtqtqJQohRJPkFHgRrZ7cRsp373JmzZMUJMcQ/d9RePYdS/NxC7H39Ld1iKIBqlRi01p/X1uBCCGaJpOdHT5X34971+uJW3YvmX9sI23XajKPbKP5+Dfx7H2zrUMUDYwMHhFC1AuOAWG0fOwbAu98FztXLwpS44l561aiF40hPyXO1uGJBkQSmxCi3jCZTHhfOYnQuQdx63odAOl7P+PkjC6k/rIKuT1JVIQkNiFEvePYrDUtH95Ii3s+wM7Nh8L0BGLfHU/0f0eSn1Sfn5Al6gNJbEKIeslkMuF1+QTCXjiMu2WUZMb+jZyc3pmUHz6Q2psokyQ2IUS95uATRPC0dbS47yPsPJpRmJVC3JJJnF4wjLyESFuHJ+ohSWxCiHrPZDLh1fc2wl44jIdllGTmoW+ImNGV5O/eldqbOIckNiFEg+Hg1ZzgqZ8QNHUN9l7NKcxOI375FKJeHkLemRO2Dk/UE5LYhBANjmfv0YS9cBjPfuMAyPpzOydndCXp20WYCwttHJ2wNUlsQogGyd6jGUGTPyT4oS+x9wnGnJvJmY8eJOqlq8iN/dvW4QkbksQmhGjQPLoPJ2zeIbwGGI+HzDr6IxEzu5O4eQHmwgIbRydsQRKbEKLBs3f3ocXdi2n52GYcmoVgzsvm7CePc2ruAHJOH7F1eKKOSWITQjQa7p2HEDr3AN6D7gMg+/huImf3InHjS5gL8m0cnagrktiEEI2KvasXgRPeotWTW3EMaIM5P5eza2cQOacfOacO2jo8UQcksQkhGiW3DlcROvcAPoMfAJOJnIjfiXi2N2c/fw5zfq6twxO1SBKbEKLRsnN2p/m4hbR+egeOLS6GgnwSv5xDxHN9yD75m63DE7Wksg8abZCUUj7Ay8AowA3YDTyitd5fgX2XARNL2bRba9235qIUQtQW14v7EzrndxLWzSbpm9fJPXWQyDn98Bv6OH43zMTOycXWIYoa1OhrbEopO2ATcBvwJvAEEAjsUEq1q+BhMoHxJf7NrvlohRC1xc7JlYDbXqb1Mz/iFNwRCgtI3PgSkc+Gk3Vsl63DEzWoKdTYbgYuA0Zprb8AUEqtAY5iJKcJFThGntZ6Za1FKISoM67t+hDy3F4S188lcdN8cqP/5NS8Afhe+xDNbpqDnZOrrUMU1dToa2wYiS0a+LJohdb6DLAGGKmUcqzIQZRS9kopz9oJUQhRl+wcnfEf/Twhs3bh1LormAtJ+noBETN7kKl/sHV4opqaQmLrAfymtS45/fcewBO4qALH8ARSgVSl1Fml1AKllDTKC9HAuYT1JHT2bpqNehbsHcmL+5uol64ifuWDFGan2zo8UUVNoSkyCNheyvoYyzIY+PMC+8dgDDzZB9gDI4CHgQ7A9SULK6WSy4nHu5ztQog6ZHJwotmNM/HoOZLYJXeTc/I3krcuIuPAJgLvfA+3joNsHaKopAaV2CwDQZwqUlZrnW350RXIKaWI9fYLHefpEqtWK6WigMeVUoO11t9WJB4hRP3m3LoLITN/JmnzayR88Rx5Z04Q9fJgvAfei/+t87F39bJ1iKKCGlRiA64AvqtIQaVUgNb6LJAFOJdSpKgpMasKcbwGPA5cDZyT2LTWPuXElYzU2oSol0z2DvgNfxL3njcQt2QS2f/sImXHe2Qc3EzgHe/g3vU6W4coKqChJba/gDsrWDbNsozBaI4sqWhddGWD0FrHKaVyAb/K7iuEqP+cgzvQesZOkre8wdnPniE/8RSnFwzDq/9EAsa+hr27r61DFBfQoBKb1joWWFbJ3fYDlymlTCUGkPQB0oFjlY1DKdUKo0n0TGX3FUI0DCY7e3yvexj37sOJ++BesvROUn9cTsbhLQRO/B8ePUbYOkRRhqYwKnItxgCRG4tWKKX8gVuAL7XWeVbr21nftK2UciljiP9My/Kb2glZCFFfOLVoT6sntxHwnzcwObtTkBxD9H9HEvPOfyhIT7B1eKIUDarGVkVrgV3Ah0qpV4GzwBSMpP5sibLbLMswy7IFsE8ptQqjGdQOY1Tk1cAnWuudtRq5EKJeMNnZ4XvNVDy6DSXug8lkHtlG2q7VZB7ZRvPxb+LZ+2ZbhyisNPoam9a6ABiKcUP2NOAVjCbEq7TW5TVDJgMbgSHAS5Z/AcCjwLhaClkIUU85BrSh5ePf0PyOd7Bz8aQgNZ6Yt24letEY8lPjbR2esDCZzSXvWxa1SSmV7Onp6b13715bhyKEqIa8hFPELZtM5iGjR8LOoxnNxy3Es+9YTCaTjaNrfMLDw0lLS0spb+Q5NIEamxBC1AbHZq1p+cgmAictxc7Nh8L0BGLfHU/0G6PIT6r0YGtRgySxCSFEFZlMJrz7TyRs3iHcLaMkM/Zt4OSMLqT8sAxpEbMNSWxCCFFNDr7BBE/7nBb3rcTOoxmFmcnELbmb0wuGkZdwytbhNTmS2IQQogaYTCa8+o4l7IXDeFhGSWYe+oaIGV1I3vGe1N7qkCQ2IYSoQQ5ezQme+glBU9dg79Wcwuw04pfdT9TLQ8g7c8LW4TUJktiEEKIWePYeTdgLh/HsdzsAWX9u5+Qz3Uja+hbmwkIbR9e4SWITQohaYu/RjKDJKwh+8AvsfYIx52RwZuU0ol4aRG7s37YOr9GSxCaEELXMo8cIwuYdwmvAHQBkHf2BiJndSdy8AHNhgW2Da4QksQkhRB2wd/ehxd1LaPnoVzj4tcacl83ZTx7n1NwB5Jw+YuvwGhVJbEIIUYfcu1xL6LyDeA+8F4Ds47uJnN2LxI0vYS7It3F0jYMkNiGEqGP2rl4E3vE/Wj3xLY4BbTDn53J27Qwin7+MnFMHbR1egyeJTQghbMSt4yBCn9+Pz+AHwGQi5+RvRDx7KQlfzMGcn2vr8BosSWxCCGFDdi4eNB+3kNZP78AxsD0U5JHwxXNEPNeH7JO/2zq8BkkSmxBC1AOuF/cn9Pl9+F73KJjsyD11kMg5fTm79hkK83JsHV6DIolNCCHqCTsnVwJue5nWz/yIU3AHKCwgceOLRM7uRdY/u20dXoMhiU0IIeoZ13Z9CHnuN/yGPw129uRG/8mpuf058/HjFOZm2Tq8ek8SmxBC1EN2js743zyXkFm7cGrdFcyFJH29gIiZPcg6+qOtw6vXJLEJIUQ95hLWk9DZu2k2cjbYO5AX9zenXhxI/MoHKczJsHV49ZIkNiGEqOdMDk40GzmL0Gd/xTmsF5jNJG9dRMQz3cg8st3W4dU7ktiEEKKBcG7dlZCZP+N/8zxMDk7knTlB1MuDiVt2PwVZqbYOr96QxCaEEA2Iyd4Bv+FPETLnd1za9QUgZcd7RMzoSsahb2wcXf0giU0IIRog5+AOtJ6xk4DbXsXk6EJ+4ilOvzaU2CV3U5CRZOvwbEoSmxBCNFAmO3t8r3uY0Of343rxAABSf1jGyRldSN+3wcbR2Y4kNiGEaOCcWrSn1VPbCfjPG5ic3SlIjiH6vyOJeec/FKQn2Dq8OieJTQghGgGTnR2+10wlbO4BXDsMAiBt12pOTu9M2q+f2Ti6uiWJTQghGhHHgDa0emILze/4H3YunhSkxhPz1hiiF40hPzXe1uHVCUlsQgjRyJhMJnwG3kvovEO4dbkWgPS9n3FyemdSd63GbDbbOMLa5WDrAGqbUkoB9wF9gB6AC9BGa32yEse4DHgZ6AmkAp8AT2utM2s8YCGEqCGOzVrT8pFNpP64nDOrH6UwPYHYd/5D2u41BE58GwefIFuHWCuaQo2tHzAN8AL+rOzOSqnuwDaMhPgIsBiYjJHchBCiXjOZTHgPuIOweYdw7zECgIx96zk5vTMpPyxrlLW3ppDY1gM+WuvOwIdV2P8FIAEYqLV+R2v9DPB/wHCl1KAajFMIIWqNg28wwdM+p8XkFdi5+1GYmUzckrs5vWAYeQmnbB1ejWr0iU1rnai1TqvKvkopL2Aw8KHWOt1q04dAOjCmBkIUQog6YTKZ8Op3O2EvHMYjfDQAmYe+IWJGF5J3vEd2Xi5RaUksOfITj/zwKW8f3IFOiiUmI5mc/DwbR19xjb6PrZq6YFyjvdYrtda5Sqn9GH12QgjRoDh4BxL8f2tI+3Ut8SseoCA1nvhl95Ozczn3Nu/OKWfP4rKv7vuW1wfcQoiHHx39gnB2cLRh5BUjie3CinpWY0rZFoPRf3cOpVRyOcf0rmZMQghRIzx734zbJQOJ/+gh0natxvn4Lt45+Rvvt72CL4N7YDaZyC0s4OEfPmX1tZOIy0ol0NWr3ie3BpXYlFJ2gFNFymqts2vglK6WZU4p27KttgshRINk7+lP0H0r2RPUlYBNc/HPzWDasW0MPKN5WV1HtKsvuYUF/BJ7nJbuPrT3aU4H3xb1Ork1qMQGXAF8V5GCSqkArfXZap6v6BnszqVsc7HaXkxr7VNOXMlIrU0IUc/sbNaWTeF3cf/x77g+9jBdU6JYvHcZS9sMYF3LnkSmJeLn4s7SIz/xWM/B9brm1tAS21/AnRUsW6UBIyUUNUGWdrNHEBBdA+cQQgib6+wXzCeOLryirmdHwCU8evQbmuekMeWf77jijCYhrCOu9o5EZ6Sw7p/9XB7UjlYePvg5u9e7BNegEpvWOhZYVoenPAzkA+HAuqKVSiknoDuwqg5jEUKIWjMkpAPP/7qJ3MICfvVrw13hdzL5+A5GxBykc2o05lX3kzX4QX7NMBPi6cfxlDOMU304mhxP38A29Sq5Nfrh/pWhlLpEKRVS9H+tdQqwFRivlPKwKjoe8AA+reMQhRCiVvi7ePDB1RNxsrMHINPBmdcvvpanu91GrncQpvxc3Da/whv7PqJjbirRGSn8Enuc/x36nris1Hp1O0CDqrFVhVLKG3jA8t+iUYz/Z+nritBar7Aq/ifwPTDQat0M4Gdgh1JqMdAKeBTYrLXeWouhCyFEnXF2cKRvizbsHP0Y30Qe4UhCDO28/bmqlaIgK53fP5xKT72NDmkxXLJqKsda9yHqqvvwc3Gvd02TjT6xAb7A8yXWPWpZfg+s4AK01r8rpa4B5gOvY8wV+T7wdA3HKYQQNuXs4EgrD1/GX9yH6Mxk1h8/wKxd69kVdxJzi570cg1kXsROnJKiuOvkj6RsjOHD7mOItHesV02TjT6xWSY7NlWwbKnltNY/ApfXYFhCCFFvOTk4EOTmzch23SkEWnv6EeLpR98WbXnr4DZa7nyfm6P24n3mH+7fOp8TvW/jw9Z9+CX2OLtij9PGq5lNR01KH5sQQojzODs4EujqxcCWF5NfWMCP0ce4ZfN7bDh9jKXtryHizg9I8A7CwVxI+z2ruHvzHLL/2V3cNHkkMYb0nJq4nbjyTI1xZuf6TCmV7Onp6b13797yCwshhI3l5OdxJjudzRGH+TMxtrjmtuTIj+w8eZgPChLw+3k5duZCzCY7fmx/FbrfBLLtHJjWzZgn/njKGTo1C8bfxaPKtbjw8HDS0tJSyrtXGJpAU6QQQoiqK+p3m6j6kZqXxXdRR/n077108gvm7o79WXLkR45lFfCk3kzbjDMMOLqNnvF/sqr7GDZ6BzCkdUcyC/J4fs9XjGnfCxd7R5q5uuPl5FJrA02kKVIIIUS5nBwc8HR0ob1Pc/KsmiY3R/xBhHcw2Q9t4Mduo8gz2eGeHM09OxYS8s2rRJyN5OOjv3JHx8u4Z/tKAtw8WbhvG3viTnIkMaZWbhOQxCaEEKJCnB0c6eDbgid6Xcu1oR255aKePNZjMKuuncSmU38x108ROXk1Mc3aAND5zy0ELxpJj7P/sCv2OD0CWrPp5CH+r9sgHvlhLbmFBZzNzqjxOCWxCSGEqDDrpsnp4dfR0t2nuGly1bWTWJeZxfhON/F+u6swOzjjmhbPqK2v0G7TXFrbm4hMSyQiLYEeAa35JfY4B85G1XiM0scmhBCi0pwcHPDEhYu8A/gh+m9+jD7Ga/u2YsaMk70jl497lQ8ObmP8gU9xPPEr6u+dBJ/ajx7yKCe8/Wnu5klkWiL+Lu41HpskNiGEEFXi7OBIe+/mPNj9ajacOEiI1f1uy478zHVdr2ahpz8ujr7838mf8MxOpdemuXx77xriM9PoH3wRfi4e5Z+okqQpUgghRJV5OLvQyt2H2y++lGFhXXC0s2NX7HHu6HgZm04eYlvUUa6+7UUWD3+e5C5DORLWh/CQTuw7c4q+LdrQzb9ljcckNTYhhBDV4uzgSHMHRwY4XUSPgNbsij3OZ8d+o1dAKNO6DuJYyhlGX3oDSz38uS6sE8v+3MWCATfjbOeAfy3U2CSxCSGEqBFODg74OTgwNKwLg1t3ICMvl/1nTxGTkUx8VipTuw3kZMpZHu5xNV5Orvg5u9XKfWyS2IQQQtQ4R3sHfOwdGNhKMbCVKl7fpVnNNz2WJH1sQgghGhVJbEIIIRoVSWxCCCEaFUlsQgghGhVJbEIIIRoVGRVZ97zS0tIIDw+3dRxCCNFgpKWlAXhVpKwktrpXCNilpaWl2jqQGuJtWabYNIr6Qa7Fv+Ra/Euuxb+qcy28MD4/yyVP0BbVopRKBqjIU20bO7kW/5Jr8S+5Fv+qq2shfWxCCCEaFUlsQgghGhVJbEIIIRoVSWxCCCEaFUlsQgghGhVJbEIIIRoVSWxCCCEaFbmPTQghRKMiNTYhhBCNiiQ2IYQQjYokNiGEEI2KTIIsyqSUUsB9QB+gB+ACtNFan6zEMS4DXgZ6AqnAJ8DTWuvMEuWcgTnAeMAXOADM0Fpvq/4rqRlKKR+M1zIKcAN2A49orfdXYN8LdWZv1VoPtpQLA06UUe56rfXXlQi51lTzWiwDJpayabfWum+JsnbAY8D9QBBwFJintf6kGuHXqKpeC8trmwjcBHQH/DB+96uA17TWOVZlw6gn74vq/K0qpVoCrwNDMCpW24GHtdbnvTal1N0Yv/s2QCTwX631WxWJURKbuJB+wDTgT8u/HpXZWSnVHdgG/AE8ArTCeKO2BUaUKL4MGA0sBI4BdwCblVJXaq1/qWL8NcbyIbQJ6AK8CiQAU4AdSqleWut/yjnE+FLWhQMPAltK2bYS+KbEugOVCrqW1MC1AMgEJpdYd6aUcvOAp4D3gL3AjcDHSqkCrfXaKr6EGlPNa+EGLAV2Ae8A8Rh/c88Dg4BrStmnPrwvllGFv1WllAfwHeCJ8XvNBx7GuFbdtdZJVmUnY1yTT4EFwABgkVLKRWv9WnkBSmITF7Ie8NFapymlHqKSiQ14AeMPfaDWOh1AKXUSeF8pNUhrvd2y7lLgNoxvbgst6z4EDgPzgSuq/1Kq7WbgMmCU1voLAKXUGowaxGxgwoV21lqvLLlOKTUQMAOrS9nlt9L2qSeqdS0s8sp7fZZv949ifFN/yLJuMfA98KpSap3WukKPMalF1bkWucDlWuufrda9b/kbeU4pNVBrvaPEPjZ9X1Tzb3UKcBHQS2u9z7LvZsu+DwOzLOtcMRLfl1rrMZZ937d8iZitlFqstb7gY2+kj02USWudqLVOq8q+SikvYDDwYVFSs/gQSAfGWK27GcgDFludOxtYAvRXSgVVJYYadjMQDXxZtEJrfQZYA4xUSjlW5mCW5pzRwPda66gyyrgrpZyqHnKtqZFroZSyV0p5XqDIjYAj8LbVeczA/4BQ4NLKh17jqnwttNa5JZJakc8tyw6l7Wfj90V1/lZvBnYVJTXLvn9htOpYfx5cBTTD6vdu8RZGbe/68oKUxCZqSxeMFoG91iu11rnAfs6t/fUA/iqRAAH2ACaM/gdb64HxbblkX9kejD+2iyp5vKGAD/BRGdufx/gCkK2U+kUpVR9qrUVq4lp4YvS5piqlziqlFiilXEo5T6rW+mgp5ynabms1/b4AaGFZni1lm63fF1X6W7XUtrpS4vPAat+LlVJuVueglLK/YTxotNzfuyQ2UVuKvrnFlLItBgguUbascpQoays1HeM4IAco2U9UiNGH8hhwg2UZCmxVSg2o5DlqS3WvRQzGYIs7gdsx+hgf5t+aivV5YqtxnrpQG+/dJzCeMG3d91pf3hdVfb1+gPMF9jXx72dGEJCjtU60LmT5UpxwgXMUkz62JsLyjalCzReWpoXqcrUsc0rZlm21vahsWeUoUbbaqngtaixGSzPtMOArrXVyifNFAteVKP8xcAR4Cbi8ouepYCx1fi201k+XWLVaKRUFPK6UGqy1/rYmzlNZtn5fWGKYjjFoZLJ1P1Jdvy8uoKqvt7zPA+syrhj9j6Up+dlRKqmxNR1XAFkV+aeU8q+B82VZls6lbHOx2l5UtqxylChbE6pyLWoyxtGW/cpqhjyH1joaY4BJX6vmmppi62tRpGik29VW65rU+0IpdSswF3hXa/1eeeVr+X1Rlqq+3vI+D6zLlHWOorLlXlOpsTUdf2E0/VRElQaMlFDU5FBaZ3IQRoe7ddmyylGibE2oyrWoyRjHYTQ1bazEPqcwvoj6YAyVrym2vhYAaK3jlFK5GE1WRWIwhnnX2HnKYbNroZQajDGwagMwtYIxQO29L8pS1debiFFbK2tfM/9+ZsQATkopP+vmSMuAmWYXOEcxSWxNhNY6FuP+k7pyGOM+lXBgXdFKy5uzO8ZNqEX2Aw8qpTxKdEr3sSxr9D6dKl6L/cBlSilTiYECfTA6849V5CCWUWNXAcusb8CtgLZAAZBUXsHKsOW1sKaUaoXRDGh9L9t+YJJS6uISA0j6WG2vMba6FkqpPhj9i78Ct2mtCypx/lp5X1zAfqrwt6q1LlRKHcL4PCipD/C31aQN+y3LcM7tZwzHSOL7KYc0RYoaoZS6RCkVUvR/S//AVmC85cbMIuMBD4wbL4usxRjWPcnqeM4Y355/sjS52NpajE7rG4tWWJqjbsG43ybPan07pVS7Mo5zG8bfXanNkEqpgFLWXQSMBXZqrWu6+a0qqnwtlFIuZQzxn2lZWt98/CXG0PIpVvubMGbDicSY4cPWqvW+UEp1wLjB+yQwoqzfbz16X1Tob1UpFaKUuqSUffsqpXpY7aswbka3/jzYjlHDm3Lu7tyP8WVhc3lBSo1NlEkp5Q08YPlvP8vy/5RSyUCE1nqFVfE/MW6cHWi1bgbwM8bMAosxZh55FNistd5aVEhrvVsp9SnwsqVG8w/GVEOhGLMa1AdrMWaI+FAp9SrGUOwpGEnq2RJli6YWCivlOOMwmlJ2lHGel5VSbS3HiAHaYXyQgzESrj6ozrVoAexTSq3CaPqzw5iF5mrgE631zqIdtdZRSqmFwGOWWwH2AiMxmidvrQc3Z0M1roUlwX+DMS3VK8Aw43O+2EGt9UHLz/XifVGJv9UPgSsxRjsWeRu4B/hKKfUaRovOIxiv53Wrc2QppWYCb1ludt+C8Tv/D/BkyQFXpZHEJi7EF+O+GWuPWpbfAyu4AK3170qpazBmJHgd476l94GSo+LAmKHhecvSFzgIDNVa/1Tl6GuQ1rpAKTUU4wNoGsbIrD3ABK11RZshFdALWHCBD+UtGB9YD2D0myRZ1j2ntf6jWi+ihlTzWiRj9C0OwfggtMOYpeNR4L+llH8K4xpMxqgVHAVu11qvqfYLqQHVvBbNgNaWn18qZftzGH8HUL/eF1X6W7XMYDQQ47NgJsbv/jvgIa11Qomybyul8jDeFzdi9CU+qLV+oyIByoNGhRBCNCrSxyaEEKJRkcQmhBCiUZHEJoQQolGRxCaEEKJRkcQmhBCiUZHEJoQQolGRxCaEEKJRkcQmhBCiUZHEJkQtU0qZlFIyy48QdURmHhGiFiilTgKLgaEYTzMYBBwCXsSY79AVY7LXqVrrs5Z9HgYewnh0SxLG1FsLSxy3M8bs5tdbPZDTZpRSN2JMYNtJa/23reMRAmSuSCFq00TgBoz5DR2A5YA90APj2VlvYjwmZbhS6mJgHtBLa/2nZYb4kFKOuQBjFvVvASxPTkih4q0vzayfcVVdWusvLY8jmQ/cVFY5pdQiYBTQSmtttjxFPJlzJ8lNBTRGQv+4oseq/qsQjY0kNiFqz7ta6z8BlFJ+wBggsGjCV6XUdCBOKeWLMdO5CeiklIqw1OLOWh9MKdUPGIxR4yvigJFArd0PXIYx63uc1fqcmkxqVv4LLFdKdSptQl7Lo2ZGYjzGpSgR9cR4vZ9gTIpswpgQ+AFgtVIqV2u9roLHEuIcktiEqD2RVj+HYXx4/13i0SQ5QIjW+oBSajzGI0+WKqV+BZ7WWu+xKjsFI9l9VbTC8giPldYHtDRpZgP/1Vrn19irKds64H/8O/t8Sb2BlsAXVut6WpartNbri1YqpX7BaKIdh9UDass5lhDnkMQmRCVZmtEeAG7FeA6VCYgCvtNaT7Uqav1omkjAjJHEUks7rtZ6LbDW8uyxJzCe9RViOacDRk1lo/XDK0uJzRHoDOyviaRmecDlK8BwjEeUxAHvAnOKakxa63Sl1A/AzZSe2EZhNDt+Z7Wul2W5r0TZfyxL7zJCOu9YSqn5GNdLAf8HjLbE+hNwl9b6lOVLw0NAB4xHoDyptf6izBcuGjQZFSlEJVieFvwD8DjGQx8fw3i21M9A+7L201rHYSSqt5RSzS3Haq6Uutnys1JKDVFKuQK5GE8KLrA6RC+MJ4/v4cI6AU6cnzAqzTJQ5RDQF+MZWk9g9IM9C0wtUfwXoEUpT00GIxl9VSIh9wTOaK1PlSg73LL8uYywSjtWDyALo4bXDONZYR8B12Bc70UYv6ePgVmWMh+V9lRq0ThIjU2IyrkB6Apcq7XeUsl978JICrstg0PigQ0YCc8J48GSnTBqdn8AY6327WhZ/sOF9bAsf69kbOew1Bq/xBjMcb3WOtOyfhVG7fNuYJHVLkVxdcJ4MnbRcTpg1KSesVrnAVwM/GC5DgABGNd2tuWcr5US03nHsuiBMcr0Ba31Kqvy4RhP594EhBclQ8sDLBdi/B63IRodSWxCVI6vZXmpUmprWU/C1lqHlbIuHaPm8Fgp2w4B/S5w3qLaRXmDP4r6rqpbY3sAo19wcFFSA9Baxyil/gBalShf9ATk5iXWj8To7/vaal13jNaiK4EzVuvzMUZ9vqS1TiklpvOOpZRqBfgD662TmkUSkIfRHGldwytqCq6L/kdhA9IUKUTlrAUOYDR3RSul3ldKDVdK1fbfUtEIQNMFSxm1l3yMJsTquAX4RWt9vJRtbpyfYIviKjlScRSw1ZLUixT1rz2FMcpzOMbtAvZAP611UhkxlXasohrqJ6WU7wzs0FrHl1jfwbLUZZxHNHCS2ISoBMtw+V7AdRgfpoMxmhN/VEo51eKpi2o2fmUVsCTXbsCfWuvsqp7IcvtBOLC3lG1+QBuMfkZrRXGdsSrbynKcL0qULapVfqS13qq13qS1fgrjhvYBSqmrSzlvecfaVaJ8a4xa7i7O1xOI1lrHlrJNNAKS2ISoJK11gdb6G631g0A7jOH2/TCSSm05bFmWOUDFss2DavavYdSCTBgDWEq6HaMLo+QN1BdZloet1o3EqMGtL1G2F3BWax1VYv17luVYzlfWsXoAyaXULIsSXmnXokcZ60UjIYlNiApSSgVYbhAuprUuwBi9aAZO1+Lp92H0DfW9QJma6l8rat672vr1KuMGvLkYAy62l9inLxCntbZu3hsF/Ki1tq7FuQGXYDTnnkNrvReIBm4opWn3vGNZxVpakipq7vzNeqVSKgyjdimJrRGTwSNCVNyrQH+l1JfAMYwvhtdi9BG9orWOrq0Ta60LlFLrgJFKKWetdU4pxSo0ItIyj2Wo1rqs/rqeGCM2Q4H1SqnNQFtgMkbf2njrQTOWUY4DgKVW6/yAKzBuEbDWDaMvbX8Z594E3AP0B3Ze6FiW9SHAmjJeQ2m3E1yoJicaCamxCVFx2zCGso/BuK9rBsZAipFa65If4LXhf4AP/97rVVIPjJrjebWhEjwwakZl6YFR67sRo5/qNeA/GPeG9dZax5QoPxrjOrxrtW44xhfnz0uULUos+8s490bLclQFjnWhRN6L0muuNXI7hKjfZHZ/IRoQpdTXgLvWekAV9++Kkfju0lp/UMp2NyANeK2iyVop9TtwUmt9k9W6z4E2WuvuVYmzxPFr7FiiaZAamxANy6NAP6XUkCrufy1GYltexvZuGJ8LFbpdQCk1EmNY/ZMlNv0CTK9aiOepyWOJJkBqbEKIYkqpKcBbQHetdXlNmkLUS1JjE0JYK7rB+09bByJEVUmNTQghRKMiNTYhhBCNiiQ2IYQQjYokNiGEEI2KJDYhhBCNiiQ2IYQQjYokNiGEEI2KJDYhhBCNiiQ2IYQQjYokNiGEEI3K/wMWSemQTvBdkgAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "sns.scatterplot(\n",
    "    data=results_df, x='s_res', y='ln_eta_nist', \n",
    "    color=colors[0], label='NIST'\n",
    ")\n",
    "sns.lineplot(\n",
    "    data=results_df, x='s_res', y='ln_eta_saft', \n",
    "    color=colors[1], label='entropy scaling'\n",
    ")\n",
    "\n",
    "plt.xlabel(r\"$s^\\text{res}(T, \\rho) / R / m$\")\n",
    "plt.ylabel(r\"$\\ln\\frac{\\eta}{\\eta_\\text{CE}}$\")\n",
    "plt.legend(frameon=False);"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "96c715a3-dd27-4afd-a146-c5c516e43ce7",
   "metadata": {},
   "source": [
    "# Summary\n",
    "\n",
    "- The `Estimator` object in FeO$_\\text{s}$ allows the collection of `DataSet` objects for adjusting parameters.\n",
    "  - The `Estimator` takes a list of `DataSet` objects, weights, and `Loss` objects as inputs.\n",
    "  - For `DataSet.viscosity`, it calculates the cost of a model by evaluating the relative difference between the model's prediction and experimental data in each `DataSet`.\n",
    "- To work with `scipy`'s `least_squares` solver, a cost function is required.\n",
    "  - Two functions, `eos_from_parameters` and `cost`, are built for this purpose.\n",
    "  - `eos_from_parameters` constructs the parameters and equation of state for the current parameter vector.\n",
    "  - `cost` calculates and returns the cost of the current model based on the parameters and estimator.\n",
    "- Initial parameter guesses and parameter bounds are necessary for parameter adjustment.\n",
    "  - Checking multiple combinations of initial parameters is recommended to avoid local minima (not shown in this notebook).\n",
    "- For `DataSet.viscosity`, the `Estimator.mean_absolute_relative_difference` method provides the mean absolute relative difference (MARD) without applying weights or losses."
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.9.12"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}