Skip to content

Instantly share code, notes, and snippets.

@goerz
Last active August 5, 2019 07:42
Show Gist options
  • Save goerz/21e46ea7b45c9514e460007de14419bd to your computer and use it in GitHub Desktop.
Save goerz/21e46ea7b45c9514e460007de14419bd to your computer and use it in GitHub Desktop.
Time Discretization in Quantum Optimal Control
name: krotov
channels:
- defaults
- conda-forge
dependencies:
- python=3.6.7
- qutip=4.3.1
- pip:
- matplotlib
- watermark
- git+https://github.com/qucontrol/krotov.git@master#egg=krotov
Display the source blob
Display the rendered blob
Raw
{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Time Discretization in Quantum Optimal Control"
]
},
{
"cell_type": "markdown",
"metadata": {
"ExecuteTime": {
"end_time": "2019-02-14T00:48:17.761438Z",
"start_time": "2019-02-14T00:48:17.756117Z"
}
},
"source": [
"<a href=\"https://mybinder.org/v2/gist/goerz/21e46ea7b45c9514e460007de14419bd/master?filepath=Krotov_time_discretization.ipynb\"><img style=\"float: left;\" src=\"https://mybinder.org/badge_logo.svg\" alt=\"Binder\"/></a>"
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {
"attributes": {
"classes": [],
"id": "",
"n": "1"
}
},
"outputs": [],
"source": [
"import numpy as np\n",
"import matplotlib.pylab as plt"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"When implementing a gradient-based control method such as [Krotov's method](https://krotov.readthedocs.io), there are some subtleties in how controls are discretized onto the time grid. In particular, the requirements for gradient-based optimization differ from the choices that the [QuTiP package](http://qutip.org) makes when simulating quantum dynamics. This document explains these details in regards to the [`krotov` Python package](https://pypi.org/project/krotov/)."
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {
"ExecuteTime": {
"end_time": "2019-02-14T00:09:37.077916Z",
"start_time": "2019-02-14T00:09:36.921132Z"
}
},
"outputs": [
{
"data": {
"text/plain": [
"'0.1.0.post1+dev'"
]
},
"execution_count": 2,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"import krotov\n",
"krotov.__version__"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"QuTiP, in its propagation routine `mesolve`, takes the control values at the points of the time grid (`tlist`), and extends them in a piecewise-constant manner in the range ±dt/2 around each grid points. That is, controls switch at the midpoint between any two points in `tlist`.\n",
"\n",
"In contrast, Krotov's method (as well as any other gradient-based method, such as GRAPE) requires pulses to be constant on the intervals of the time grid, and switch *on* the time grid points. The implementation of Krotov's method handles this transparently: On input and output, the controls in the [`Objectives`](https://krotov.readthedocs.io/en/latest/API/krotov.objectives.html#krotov.objectives.Objective) follow the convention in QuTiP. Controls can be either functions that will be evaluated at the points in `tlist`, or an array of values of the same length as `tlist`. Internally, Krotov uses the functions [`control_onto_interval`](https://krotov.readthedocs.io/en/latest/API/krotov.structural_conversions.html#krotov.structural_conversions.control_onto_interval) and [`pulse_onto_tlist`](https://krotov.readthedocs.io/en/latest/API/krotov.structural_conversions.html#krotov.structural_conversions.pulse_onto_tlist) to convert the control values to be constant on the intervals of the time grid (resulting in an array of pulse values that has one value less than the points in `tlist`), and back:"
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {
"attributes": {
"classes": [],
"id": "",
"n": "2"
}
},
"outputs": [],
"source": [
"from krotov.structural_conversions import pulse_onto_tlist, control_onto_interval"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The names of these routines reflects a convention used in the `krotov` package, where \"controls\" are defined on `tlist`, and \"pulses\" are defined on the intervals of `tlist`.\n",
"\n",
"The transformation between the two is invertible. Since we generally care a lot about the boundaries of the controls (which often should be exactly zero), the initial and final value of the control is kept unchanged when converting to midpoints. For any other point, the \"control\" values at the time grid points will be the average of the \"pulse\" values of the intervals before and after that point on the time grid.\n",
"\n",
"We can illustrate both sampling conventions in a plot:"
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {
"attributes": {
"classes": [],
"id": "",
"n": "3"
}
},
"outputs": [],
"source": [
"def plot_sampling_points(func, N=30, T=10, peg='gridpoints'):\n",
" \"\"\"Show the numerical sampling for a control `func`\n",
"\n",
" Args:\n",
" func (callable): Callable control ``func(t, None)`` for every point on\n",
" the time grid\n",
" N (int): number of points on the time grid\n",
" T (float): final time\n",
" peg (str): If 'midpoints', make the sampling exact at the\n",
" midpoints of the time grid. If 'gridpoints', make the sampling\n",
" exact at the points of the time grid.\n",
" \"\"\"\n",
" tlist = np.linspace(0, T, N)\n",
" dt = tlist[1] - tlist[0]\n",
" tlist_midpoints = []\n",
" for i in range(len(tlist) - 1):\n",
" tlist_midpoints.append(0.5 * (tlist[i+1] + tlist[i]))\n",
" tlist_midpoints = np.array(tlist_midpoints)\n",
" tlist_exact = np.linspace(0, 10, 100)\n",
"\n",
" if peg == 'midpoints':\n",
" pulse = np.array([func(t, None) for t in tlist_midpoints])\n",
" control = pulse_onto_tlist(pulse)\n",
" elif peg == 'gridpoints':\n",
" control = np.array([func(t, None) for t in tlist])\n",
" pulse = control_onto_interval(control)\n",
" else:\n",
" raise ValueError(\"Invalid peg\")\n",
"\n",
" fig, ax = plt.subplots(figsize=(9, 7))\n",
" ax.plot(\n",
" tlist_exact, func(tlist_exact, None), '-', color='black',\n",
" label='exact')\n",
" ax.errorbar(\n",
" tlist_midpoints, pulse, fmt='o', xerr=dt/2, color='blue',\n",
" label=\"sampled on intervals\")\n",
" ax.errorbar(\n",
" tlist, control, fmt='o', xerr=dt/2, color='red',\n",
" label=\"sampled on time grid\")\n",
" ax.legend()\n",
" ax.set_xlabel('time')\n",
" ax.set_ylabel('pulse amplitude')\n",
" plt.show(fig)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"In the following we will use a Blackman shape between t=0 and T=10 as an example:"
]
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {
"attributes": {
"classes": [],
"id": "",
"n": "4"
}
},
"outputs": [],
"source": [
"func = krotov.shapes.qutip_callback(krotov.shapes.blackman, t_start=0, t_stop=10)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"If we performed an optimization using this Blackman shape as a guess control, with\n",
"N=10 time grid points, the sampling would look as follows:"
]
},
{
"cell_type": "code",
"execution_count": 6,
"metadata": {
"attributes": {
"classes": [],
"id": "",
"n": "5"
}
},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAiwAAAGtCAYAAAA1cy8JAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAIABJREFUeJzs3XlcVOX+B/DPwyKKoubSlUQBu6QopCK45paV5laJhCwqaGJqeW/en9pipV4tb4t1K8sLLrgAGm55kyx3c6HENksNFUHN3Mg1AZH5/v4AuSgoI8zwzPJ5v17zwjlzZuYzOMBnnnPOc5SIgIiIiMiSOegOQERERFQeFhYiIiKyeCwsREREZPFYWIiIiMjisbAQERGRxWNhISIiIovHwkJEREQWj4WFiIiILB4LCxEREVk8J90B7laDBg3Ey8tLdwwiIiIygb17954TkYblrWd1hcXLywtpaWm6YxAREZEJKKWyjFmPm4SIiIjI4rGwEBERkcVjYSEiIiKLZ3X7sBARkX75+fk4ceIEcnNzdUchK1G9enV4eHjA2dm5QvdnYSEiort24sQJuLm5wcvLC0op3XHIwokIsrOzceLECXh7e1foMbhJiIiI7lpubi7q16/PskJGUUqhfv36lRqRY2EhIqIKYVmhu1HZ9wsLCxEREVk8FhYiIqIStm7dil27dumOQbdgYSEiIiqBhcUysbAQEZFVWrp0Kdq3b482bdpg9OjRyMrKgo+PD86dOweDwYCuXbviq6++AgA8+eSTaNeuHVq1aoXY2Njix1i/fj0CAgLQunVr9OrVC5mZmZg7dy7ee+89tGnTBl9//bWul0e34GHNRERUKX//+9/xww8/mPQx27Rpg/fff/+2tx84cADLly/Hzp074ezsjLFjx2Lbtm2YPHkyxowZg/bt26Nly5Z47LHHAAALFixAvXr1kJOTg6CgIAQHB8NgMGDUqFHYvn07vL298ccff6BevXp49tlnUatWLfzf//2fSV8TVQ4LCxERWZ1NmzZh7969CAoKAgDk5OTg3nvvxdSpU5GcnIy5c+feVKI++OADrF69GgBw/PhxHDp0CGfPnkW3bt2K5wWpV69e1b8QMhoLCxFZlIQE4JVXgGPHgKZNgZkzgYgI3anoTu40EmIuIoLhw4fjzTffvGn51atXceLECQDAlStX4Obmhq1bt2Ljxo3YvXs3XF1d0aNHD87Qa4W4DwsRWYyEBCAmBsjKAkQKv8bEFC4nKqlXr15YsWIFzpw5AwD4448/kJWVhcmTJyMiIgLTp0/HqFGjAAAXL17EPffcA1dXVxw8eBCpqakAgI4dO2L79u04evRo8WMAgJubGy5fvqzhVdGdcISFiCqsRw/TPl5qKpCXd/Oyq1eBkSOBuDjTPc/WraZ7LNKjZcuWmDFjBh577DEYDAY4Oztj9uzZ2LNnD3bu3AlHR0esXLkSCxcuRHh4OObOnQtfX180b94cHTt2BAA0bNgQsbGxGDRoEAwGA+69915s2LABAwYMwODBg/HZZ5/hww8/RNeuXTW/WgIAJSK6M9yVwMBASUtL0x2DiGD6wrJt2+1v697ddM/DwlJ5Bw4cgK+vr+4YZGXKet8opfaKSGB59+UICxFVmKn/8Ht5FW4GupWnJ0sGkb3jPixEZDFmzgRcXW9e5upauJyI7BsLCxFZjIgIIDa2cERFqcKvsbE8SoiIuEmIiCxMRAQLChGVxhEWIiIisngsLERERGTxWFiIiIjK0aNHD9zNlBpbt25F//79TZ7jmWeewf79+++4zpo1a8pdxxTi4+Px3HPPmf15bmBhISIis0tIKDxs3cGh8CtnL66YefPmoWXLlndcpyKF5fr165WJVSVYWIiIyKzMccqFP//8E/369UPr1q3h5+eH5cuXAwCmT5+OoKAg+Pn5ISYmBjcmR+3RowdeeOEFBAYGwtfXF3v27MGgQYPg4+ODKVOmAAAyMzPRokULREREwNfXF4MHD8bVq1dLPfdXX32FTp06ISAgACEhIbhy5QoAYP369WjRogUCAgKwatWqMnPn5uYiOjoa/v7+aNu2LbZs2QKgcLRi0KBB6NOnD3x8fDBp0qQy719ypKdWrVp45ZVX0Lp1a3Ts2BGnT5/Grl27sHbtWkycOBFt2rTBkSNHcOTIEfTp0wft2rVD165dcfDgQQBAVFQUnn32WXTo0AGTJk2Cl5cXLly4UPxcPj4+OH36NP773/+iQ4cOaNu2LR555BGcPn26VK7k5GT4+fmhdevW6NatW/n/gRUhIlZ1adeunRARkV779+8vd53u3QsvLi4ihVXl5ouLS+HtFbFixQp55plniq9fuHBBRESys7OLl0VGRsratWuLsnSXSZMmiYjI+++/L+7u7nLy5EnJzc2Vxo0by7lz5+To0aMCQHbs2CEiItHR0fL2228X33/Pnj1y9uxZ6dq1q1y5ckVERGbNmiXTpk2TnJwc8fDwkPT0dDEYDBISEiL9+vUrlfudd96R6OhoERE5cOCANGnSRHJycmThwoXi7e0tFy5ckJycHGnatKkcO3asjO9pYQ4REQDFr2/ixInyz3/+U0REhg8fLsnJycX3efjhhyU9PV1ERFJTU6Vnz57F6/Xr10+uX78uIiLjx4+XBQsWFK/Xq1cvERH5448/xGAwiIhIXFycTJgwQUREFi5cKOPGjRMRET8/Pzlx4oSIiJw/f77s/zQp+30DIE2M+PvPERYiIjKrW88PVd5yY/j7+2PDhg2YPHkyvv76a9SpUwcAsGXLFnTo0AH+/v7YvHkzfvnll+L7DBw4sPi+rVq1gru7O1xcXNCsWTMcP34cANCkSRN06dIFABAZGYkdO3bc9LypqanYv38/unTpgjZt2mDRokXIysrCwYMH4e3tDR8fHyilEBkZWWbuHTt2FN/WokULeHp6Ij09HUDhCR3r1KmD6tWro2XLlsgqa9rnEqpVq1a8n0y7du2QmZlZap0rV65g165dCAkJQZs2bTB69Gj8/vvvxbeHhITA0dERABAaGlo8UrVs2TKEhoYCAE6cOIHevXvD398fb7/99k3f0xu6dOmCqKgoxMXFoaCg4I65K4rzsBARkVncOJ2COU658MADD+C7775DSkoKpkyZgl69emHSpEkYO3Ys0tLS0KRJE0ydOhW5ubnF93FxcQEAODg4FP/7xvUb+3AopW56nluviwgeffRRJCUl3bT8hx9+qNgLKaFkJkdHx3L3K3F2di7Od7v1DQYD6tate9t8NWvWLP53p06dcPjwYZw9exZr1qwp3lT2/PPPY8KECRg4cCC2bt2KqVOnlnqcuXPn4ptvvsG6devQrl077N27F/Xr1y/3Nd8NjrAQEZFZmeOUCydPnoSrqysiIyMxceJEfPfdd8XlpEGDBrhy5QpWrFhx14977Ngx7N69GwCQmJiIhx566KbbO3bsiJ07d+Lw4cMACvelSU9PR4sWLZCZmYkjR44AQKlCc0PXrl2RULTzTnp6Oo4dO4bmzZvfdc47cXNzw+XLlwEAtWvXhre3N5KTkwEUFq4ff/yxzPsppfDUU09hwoQJ8PX1LS4cFy9eROPGjQEAixYtKvO+R44cQYcOHTB9+nQ0bNiweMTKlFhYiIjIrMxxyoV9+/ahffv2aNOmDaZNm4YpU6agbt26GDVqFPz8/NC7d28EBQXd9eM2b94cc+bMga+vL86fP48xY8bcdHvDhg0RHx+PsLAwPPjgg+jUqRMOHjyI6tWrIzY2Fv369UNAQADuvffeMh9/7NixMBgM8Pf3R2hoKOLj428aWTGFIUOG4O2330bbtm1x5MgRJCQkYP78+WjdujVatWqFzz777Lb3DQ0NxdKlS4s3BwHA1KlTERISgnbt2qFBgwZl3m/ixInw9/eHn58fOnfujNatW5v0NQGAkqI9qE3+wEotANAfwBkR8SvjdgXg3wD6ArgKIEpEvivvcQMDA+VujoUnIiLTO3DgAHx9fXXHMKnMzEz0798fP//8s+4oNqus941Saq+IBJZ3X3OOsMQD6HOH2x8H4FN0iQHwiRmzEBERkRUz2063IrJdKeV1h1WeALC46JCmVKVUXaWUu4j8fof7EJENEhEcPnwYmzdvxtGjR2+6rXbt2ujWrRvat2+PatWqaUpI9sDLy4ujKxZM51FCjQGU3CvnRNGyUoVFKRWDwlEYNG3atErCEZF5iQg+//xzrFy5Eps3by7eSc/Z2RkODv8b/M0rOva1Zs2a6Nq1Kx577DGMHDkStWvX1pKbiPSwip1uRSRWRAJFJLBhw4a64xBRJYgI1q9fj6CgIAwcOBCff/45OnTogI8//hi//vor8vLykJubW3zJzs7GypUrERUVhczMTEyYMAHe3t546623ypyFlIhsk87C8huAJiWuexQtIyIbtWPHDnTr1g2PP/44srOzsXDhQpw6dQrJyckYM2YMHnjggVLzXtSrVw+DBg3CRx99hAMHDiAtLQ3t27fH5MmT0axZM3z44YdWcR4UIqocnYVlLYBhqlBHABe5/wqRbSooKMDrr7+Obt26ISMjo3g0JSoqCk5Od7dlul27dvjiiy/w9ddfo0WLFhg/fjx69ep10+ydRGR7zFZYlFJJAHYDaK6UOqGUGqmUelYp9WzRKikAMgAcBhAHYKy5shCRPufOnUPfvn0xffp0DBs2DIcOHcKYMWMqvQPtQw89hC1btmDJkiXYs2cP2rZti+3bt5soNdHNSp500Bhbt24tnjbflN5///2bNoX27dv3phMWVrXXXnsNGzduLLXcHK/fbIVFRMJExF1EnEXEQ0Tmi8hcEZlbdLuIyDgRuV9E/EWEk6sQ2ZhvvvkGAQEB2LZtG2JjY7Fw4UK43jrlaSXcOGfLt99+izp16uDhhx/GO++8A3PNL0WVkJBQOEe/g0Ph18qcqtmO3VpYUlJSULduXS1ZCgoKMH36dDzyyCNV8nxWsdMtEVmflJQUdOvWDY6Ojti5cydGjRpVav8UU/Hz88OePXvw5JNPYuLEiYiJiYHBYDDLc1EFJCQAMTGFJxQSKfwaE1Op0vLnn3+iX79+aN26Nfz8/IpP2jd9+nQEBQXBz88PMTExxeW1R48eeOGFFxAYGAhfX1/s2bMHgwYNgo+PT/E5czIzM9GiRQtERETA19cXgwcPLnPH7q+++gqdOnVCQEAAQkJCcOXKFQDA+vXr0aJFCwQEBGDVqlVl5s7NzUV0dDT8/f3Rtm1bbNmyBQAQHx+PQYMGoU+fPvDx8cGkSZNK3feDDz7AyZMn0bNnT/Ts2RNA4aHY586dK84eFRWFBx54ABEREdi4cSO6dOkCHx8ffPvtt8XftxEjRqB9+/Zo27ZtmbPeGgwGjB07Fi1atMCjjz6Kvn37Fp/mwMvLC5MnT0ZAQACSk5MRFRVVfJsxr79SjDmlsyVd2rVrd9vTVhORZVi/fr24uLhIQECAnDt3rsqe12AwyMsvvywA5NlnnxWDwVBlz21v9u/fX/5K3bsXXlxcRAqrys0XF5fC2ytgxYoV8swzzxRfv3DhgoiIZGdnFy+LjIyUtWvXFkXpLpMmTRIRkffff1/c3d3l5MmTkpubK40bN5Zz587J0aNHBYDs2LFDRESio6Pl7bffLr7/nj175OzZs9K1a1e5cuWKiIjMmjVLpk2bJjk5OeLh4SHp6eliMBgkJCRE+vXrVyr3O++8I9HR0SIicuDAAWnSpInk5OTIwoULxdvbWy5cuCA5OTnStGlTOXbsWKn7e3p6ytmzZ0tdP3r0qDg6OspPP/0kBQUFEhAQINHR0WIwGGTNmjXyxBNPiIjISy+9JEuWLBERkfPnz4uPj0/xa7khOTlZHn/8cSkoKJDff/9d6tatK8nJycXP969//at43eHDh0tycrLRr7+s9w2ANDHi7z9HWIjIpDZt2oQnn3wSvr6+2LBhg8nP2HonSinMmDEDkyZNwty5czF+/HhuHrIERXPpGL3cCP7+/tiwYQMmT56Mr7/+GnXq1AEAbNmyBR06dIC/vz82b96MX375pfg+AwcOLL5vq1at4O7uDhcXFzRr1qx4HqAmTZqgS5cuAIDIyEjs2LHjpudNTU3F/v370aVLF7Rp0waLFi1CVlYWDh48CG9vb/j4+BRvqizLjh07im9r0aIFPD09kZ6eDgDo1asX6tSpg+rVq6Nly5bIKusU13fg7e0Nf39/ODg4oFWrVujVqxeUUvD390dmZiaAwtGhWbNmoU2bNujRowdyc3Nx7NixUhlDQkLg4OCARo0aFY/m3FDyPEM3GPv6K0PnxHFEZGO2bt2KAQMGwMfHBxs2bEC9evWqPINSCrNmzcL169cxe/ZsODs749133zXb5ii6g61bC796eRVuBrqVp+f/1rlLDzzwAL777jukpKRgypQp6NWrFyZNmoSxY8ciLS0NTZo0wdSpU4vP4Ayg+CSDDg4ON51w0MHBofjQ+FvfJ7deFxE8+uijpc7G/MMPP1TodZRUMpOjo+NdH65/62sq+XpvPJaIYOXKlZU6Q3TNmjUrfN/K4AgLEZnE3r170b9/f3h7e2Pjxo23PatrVVBK4Z133sH48ePx3nvvYdq0adqyEICZM4Fbd7Z2dS1cXkEnT56Eq6srIiMjMXHiRHz33XfF5aRBgwa4cuVK8b4Vd+PYsWPYvXs3ACAxMREPPfTQTbd37NgRO3fuxOHDhwEU7hOSnp6OFi1aIDMzE0eOHAGAUoXmhq5duyKhaN+d9PR0HDt27K7Kg5ubGy5fvnzXr+uG3r1748MPPyweefz+++9LrdOlSxesXLkSBoMBp0+fxlYjSqWxr78yWFiIqNLOnDmDp556CvXr18emTZtw77336o4EpRTef/99REVFYdq0aVi9erXuSPYrIgKIjS0cUVGq8GtsbOHyCtq3bx/at2+PNm3aYNq0aZgyZQrq1q2LUaNGwc/PD71790ZQUNBdP27z5s0xZ84c+Pr64vz58xgzZsxNtzds2BDx8fEICwvDgw8+iE6dOuHgwYOoXr06YmNj0a9fPwQEBNz2Z2Ds2LEwGAzw9/dHaGgo4uPjbxoZKU9MTAz69OlTajONsV599VXk5+fjwQcfRKtWrfDqq6+WWic4OBgeHh5o2bIlIiMjERAQULzJ7XaMff2Voaxt+25gYKDczbHwRGRe169fx6OPPorU1FTs3LkTAQEBuiPdJDc3F927d8f+/fvx7bffljq1PVXMgQMHbO57mZmZif79+/MEiACuXLmCWrVqITs7G+3bt8fOnTvRqFGjSj9uWe8bpdReEQks774cYSGiu1diTo1L9erBfetWxMbGWlxZAQo/+a1cuRKurq6Y9/DDMDRtyrlAiMrRv39/tGnTBl27dsWrr75qkrJSWdzplojuzo05NYrmp6h3+TLinZxQzcFyP/94eHhga0wMms6Y8b9PaTfmAgEqtWmCbIeXlxdHV4oYs99KVWNhIbIHPXqY7rFSU0sdjlrt+nVg5EggLs50z2NivqmppRdevWr63Bb4i95cRIRHX5HRKrsLiuV+JCIiy2SGOTWqhLXmtlDVq1dHdnY257kho4gIsrOzUb169Qo/BkdYiOyBKT/1m2FOjSphrbktlIeHB06cOIGzZ8/qjkJWonr16vDw8Kjw/VlYiOiuHBk5Eo1eew03TR1VyTk1qsTMmTftewMA15ycUM3Sc1soZ2dneHt7645BdoSbhIjIaLm5uRiQlITJdeuiwMPDZHNqVIlb5gLJrlUL0devYxf/6BJZBRYWIjLaa6+9hgMHDmDAsmVwPH4cMBiAzEzLLys3REQU5jUYUO3kSez09ERUVFSZZ+QlIsvCwkJERtm9ezfeeecdxMTEoHfv3rrjVJqbmxsWLlyIQ4cO4eWXX9Ydh4jKwZluiahc169fR0BAAC5evIiff/4Zbm5uuiOZzLhx4/DJJ59g7969aNu2re44RHaHM90SkcnMnTsX+/btw+zZs22qrADAzJkz0aBBAzz//PM8RJfIgrGwENEdnT17Fq+++ip69eqFQYMG6Y5jcnXr1sWbb76JnTt3IjExUXccIroNFhYiuqNXXnkFV65cwQcffGCzs5pGR0cjKCgIEydOxOXLl3XHIaIysLAQ0W2lpaVh3rx5eP7559GyZUvdcczGwcEBH374IX7//XfMmDFDdxwiKgN3uiWiMhkMBnTp0gUZGRlIT09HnTp1dEcyu+joaCQkJGDfvn1o3ry57jhEdoE73RJRpSxduhSpqamYNWuWXZQVAJg1axZq1KiBv//977qjENEtWFiIqJScnBy8/PLLCAoKwvDhw3XHqTJ/+ctf8Prrr2P9+vXYuHGj7jhEVAILCxGV8vHHH+O3337Dv/71Lzg42NeviXHjxqFp06Z4+eWXeZgzkQWxr99ERFSuS5cu4c0338Sjjz6Knj176o5T5VxcXDB16lTs2bMHq1ev1h2HiIqwsBDRTd59911kZ2fjjTfe0B1Fm6FDh6JFixaYMmUKCgoKdMchIrCwEFEJZ8+exezZsxEcHIzAwHJ32rdZTk5OmDFjBg4cOIAlS5bojkNEYGEhohLefPNNXL16Ff/85z91R9Fu0KBBCAwMxOuvv468vDzdcYjsHgsLEQEAjh8/jo8//hjDhw+Hr6+v7jjaKaXwxhtv4NixY/jPf/6jOw6R3WNhISIAwLRp0yAimDp1qu4oFuORRx5Bz549MWPGDFy5ckV3HCK7xsJCRMjMzER8fDxGjx6Npk2b6o5jMW6Mspw9exZz587VHYfIrrGwEBHeeustODg4YNKkSbqjWJyOHTuiV69eePfdd5Gbm6s7DpHdYmEhsnO///47FixYgKioKHh4eOiOY5FeeeUVnDp1CgsWLNAdhchusbAQ2bl3330X+fn5mDx5su4oFqtHjx7o1KkT3nrrLeTn5+uOQ2SXWFiI7Fh2djbmzp2LsLAw3H///brjWCylFF555RVkZWUhISFBdxwiu8TCQmTH/v3vf+PPP//ESy+9pDuKxevbty9at26NN998k7PfEmnAwkJkpy5duoQPP/wQTz31FFq1aqU7jsVTSuHll19Geno6Vq5cqTsOkd1hYSGyU5988gkuXLiAV155RXcUqxEcHIzmzZvjjTfe4JmciaoYCwuRHcrJycHs2bPRu3dvtGvXTnccq+Ho6IgXX3wRP/74I7744gvdcYjsCgsLkR1KSEjAmTNneGRQBURERKBx48Z49913dUchsissLER2RkQwe/ZstG3bFj169NAdx+o4Oztj/Pjx2Lx5M3788UfdcYjsBgsLkQ1LSAC8vAAHh8KvCQnAl19+iQMHDmDChAlQSumOaJVGjRqFmjVr4r333iteVtb3mohMR1nbjmOBgYGSlpamOwaRxUtIAGJigKtX/7fM1RVo1mwWsrM/QGZmJqpVq6YvoJUbP3485s6di6ysLGze7F7m9zo2FoiI0JeRyBoopfaKSGB56zlVRRgiKp+pt86kpgJ5eTcvu3oV+PnnIfD2Fjz2WOXLytatlX4ILUzxvc7J+Rvy8z9Cx45zcPr0jDK/1yNHAnFxlX+uG6z1+01kCtwkRGSjbv0D+j9N4e4+uiqj2KQaNe5HgwZP4uTJT5CXV/ZI9e3/D4jobnGTEJGN8vICsrJKL69VKxuXL9ev8jy2aMeOHejatSvq1buEP/5wK3W7pyeQmVn1uYisibGbhDjCQmSjZs4s3I/iZn9i6tRrOuLYpC5duiAoKAjVqk2Fq+vNH/5cXQv/D4jINFhYiGxUREThTp+enoBSAgeH42jb9hP84x/uuqPZDKUUXnjhBZw6NRtjxnxf9L0u/J5zh1si02JhIbJhERGFmyTmz4+HwdAU77wToDuSzRk8eDAaN26MffteQmYmYDAUfs9ZVohMi4WFyA58/PHHaNmyJXr27Kk7is1xdnbG6NGj8dVXX+HQoUO64xDZLBYWIhv37bffIi0tDWPHjuVEcWYyatQoODs745NPPtEdhchmsbAQ2bg5c+agVq1aGDp0qO4oNqtRo0YIDg7GwoULcbXk7HFEZDIsLEQ27Ny5c1i+fDmGDRuG2rVr645j08aNG4cLFy4gMTFRdxQim8TCQmTD5s+fj7y8PIwdO1Z3FJvXpUsX+Pv7Y86cObC2+a2IrAELC5GNKigowNy5c9G9e3e0atVKdxybp5TCuHHj8MMPPyA1NVV3HCKbY9bCopTqo5T6VSl1WCn1Yhm3N1VKbVFKfa+U+kkp1deceYjsyRdffIHMzEyMGzdOdxS7ERERgdq1a2POnDm6oxDZHLMVFqWUI4A5AB4H0BJAmFKq5S2rTQHwqYi0BTAEwMfmykNkb+bMmQN3d3c8+eSTuqPYjVq1aiEqKgrJyck4c+aM7jhENsWcIyztARwWkQwRuQZgGYAnbllHANzYE7AOgJNmzENkN44cOYL169cjJiYGzs7OuuPYlTFjxuDatWuYN2+e7ihENsWchaUxgOMlrp8oWlbSVACRSqkTAFIAPG/GPER2Iy4uDo6Ojhg1apTuKHanRYsWePjhhxEXFweDwaA7DpHN0L3TbRiAeBHxANAXwBKlVKlMSqkYpVSaUirt7NmzVR6SyJpcu3YNCxcuRP/+/dG48a2fEagqxMTEIDMzExs3btQdhchmmLOw/AagSYnrHkXLShoJ4FMAEJHdAKoDaHDrA4lIrIgEikhgw4YNzRSXyDb897//xZkzZzi6otGTTz6J+vXrIy4uTncUIpthzsKyB4CPUspbKVUNhTvVrr1lnWMAegGAUsoXhYWFQyhElRAXFwcPDw/06dNHdxS75eLiguHDh2PNmjU4ffq07jhENsFshUVErgN4DsCXAA6g8GigX5RS05VSA4tW+weAUUqpHwEkAYgSzrhEVGGZmZn46quvMGLECDg6OuqOY9eeeeYZXL9+HYsWLdIdhcgmKGvrB4GBgZKWlqY7BpFFevXVVzFz5kxkZmaiadOmuuPYvW7duuHUqVP49ddfeeJJottQSu0VkcDy1tO90y0Rmcj169exYMEC9OnTh2XFQowaNQqHDh3Ctm3bdEchsnosLEQ2IiUlBSdPnuTOthZk8ODBqFu3LmJjY3VHIbJ6LCxENiIuLg6NGjVC//4uitL0AAAgAElEQVT9dUehIjVq1MDQoUOxcuVKZGdn645DZNVYWIhswIkTJ5CSkoLo6GjObGthRo0ahWvXrmHJkiW6oxBZNRYWIhuwcOFCGAwGjBw5UncUuoW/vz86dOiAuLg4WNtBDkSWhIWFyMoZDAbEx8ejZ8+euP/++3XHoTKMHDkS+/fvx549e3RHIbJaLCxEVu7rr79GRkYGoqOjdUeh2wgNDUWNGjWwcOFC3VGIrBYLC5GVW7hwIdzc3BAcHKw7Ct1G7dq1ERwcjKSkJOTk5OiOQ2SVWFiIrNjly5eRnJyM0NBQuLq66o5DdxAdHY2LFy9i9erVuqMQWSUWFiIrlpycjKtXr3JzkBXo0aMHvLy8uFmIqIJYWIis2MKFC9G8eXN06tRJdxQqh4ODA4YPH45Nmzbh2LFjuuMQWR0WFiIrdejQIezYsQNRUVE8T42VGD58OESEJ0QkqgAWFiIrFR8fDwcHBwwbNkx3FDKSt7c3evbsifj4eBgMBt1xiKwKCwuRFSooKMCiRYvQp08f3Hfffbrj0F0YMWIEMjIy8PXXX+uOQmRVWFiIrNCGDRvw22+/cWdbKzRo0CDUrl2bO98S3SUWFiIrtGjRItSrVw8DBgzQHYXukqurK0JDQ5GcnIwrV67ojkNkNVhYiKzMxYsXsWbNGoSFhcHFxUV3HKqA4cOH4+rVq1i1apXuKERWg4WFyMqsWLECubm53NnWinXu3BnNmjXD4sWLdUchshosLERWZvHixWjevDmCgoJ0R6EKUkph2LBh2Lx5M44fP647DpFVYGEhsiJHjx7F9u3bMWzYMM69YuWGDh0KEUFCQoLuKERWgYWFyIosXboUABAZGak5CVVWs2bN8NBDD2Hx4sUQEd1xiCweCwuRlRARLF68GD179kTTpk11xyETGDZsGA4cOIC9e/fqjkJk8VhYiKxEamoqDh8+zJ1tbUhISAhcXFy48y2REVhYiKzE4sWLUaNGDQQHB+uOQiZSt25dDBw4EElJSbh27ZruOEQWjYWFyArk5eVh+fLleOqpp+Dm5qY7DpnQsGHDcO7cOaxfv153FCKLxsJCZAXWrVuH8+fPc3OQDerduzcaNmyIJUuW6I5CZNFYWIiswOLFi+Hu7o5evXrpjkIm5uzsjPDwcKxduxbnz5/XHYfIYrGwEFm47OxspKSkICwsDE5OTrrjkBkMHToU165dw4oVK3RHIbJYLCxEFi45ORn5+fmce8WGBQQEoHnz5sXz7BBRaSwsRBYuISEBLVu2RJs2bXRHITNRSiEyMhLbt2/HsWPHdMchskgsLEQW7OjRo9ixYwciIyM5Fb+NCw8PBwAkJiZqTkJkmVhYiCzYjT9eN/6Yke1q1qwZOnfujKVLl3KqfqIysLAQWSgRwdKlS9G1a1d4enrqjkNVIDIyEr/88gt++ukn3VGILA4LC5GF+v7773Hw4EFERETojkJVJCQkBE5OTtz5lqgMLCxEFiohIQHOzs4ICQnRHYWqSIMGDfD4448jKSkJBQUFuuMQWRQWFiILVFBQgMTERPTt2xf16tXTHYeqUEREBH777Tds27ZNdxQii8LCQmSBNm/ejFOnTnHuFTs0YMAAuLm5cbMQ0S1YWIgsUEJCAmrXro3+/fvrjkJVzNXVFcHBwVi5ciVycnJ0xyGyGCwsRBYmJycHq1atwuDBg1G9enXdcUiDiIgIXLp0CZ9//rnuKEQWg4WFyMJ8/vnnuHz5Mo8OsmM9e/ZEo0aNkJSUpDsKkcVgYSGyMImJiXB3d0f37t11RyFNHB0dMWTIEKxbtw4XLlzQHYfIIrCwEFmQ8+fPIyUlBUOGDIGjo6PuOKRReHg4rl27hlWrVumOQmQRWFiILMiqVatw7do1TsVPCAwMxF//+leeW4ioCAsLkQVJTEyEj48P2rVrpzsKaaaUQnh4ODZv3ozff/9ddxwi7VhYiCzEyZMnsWXLFoSHh/PMzAQACAsLg4hg+fLluqMQaVduYVFKPaCU2qSU+rno+oNKqSnmj0ZkX5YvXw4RQVhYmO4oZCFatGiBgIAAbhYignEjLHEAXgKQDwAi8hOAIeYMRWSPEhMT0a5dOzRv3lx3FLIg4eHh2LNnDw4dOqQ7CpFWxhQWVxH59pZl180RhsguJSQgv3FjfJOWhs0ZGUBCgu5EZEFCQ0MRDqBhUBDg4AB4efE9QnbJyYh1ziml7gcgAKCUGgyAe4ARmUJCAhATA+erVwEAtc+fB2JiCm/jxHEEwGPbNsx3cED1ixcLF2Rl8T1CdkmJyJ1XUKoZgFgAnQGcB3AUQKSIZJo9XRkCAwMlLS1Nx1MTAT16mPbxUlOBvLzSy11cgI4dTfc8W7ea7rHozqzxPcL3B2mklNorIoHlrVfuCIuIZAB4RClVE4CDiFw2RUAiQtl/iO60nOwP3yNEAO5QWJRSE26zHAAgIrPNlInIcpn6k6iXV+EQ/608Pfmp11rxPUJkFnfa6dat6BIIYAyAxkWXZwEEmD8ake0zzJiBq7fOueLqCsycqScQWZ6ZMwvfEyXxPUJ26LaFRUSmicg0AB4AAkTkHyLyDwDtADStqoBEtmxH06Z4RgR/NmgAKFX4qTk2ljtT0v9ERACxsTA0aQIDgOxatfgeIbtkzGHNfwFwrcT1a0XLiKiSkpKS8JmrK5CZCRgMhV/5h4huFREBh2PHEB4aiuYuLsh/+mndiYiqnDGFZTGAb5VSU5VSUwF8A2CRWVMR2YH8/HwkJydj4MCBqFmzpu44ZAXCw8ORnZ2NjRs36o5CVOXKLSwiMhNANAoPaT4PIFpE3jDmwZVSfZRSvyqlDiulXrzNOk8rpfYrpX5RSnH+abIbGzZsQHZ2NqfiJ6P17t0bdevWRVJSku4oRFWu3MOalVJNAZwDsLrkMhE5Vs79HAHMAfAogBMA9iil1orI/hLr+KBw2v8uInJeKXVvxV4GkfVJSkpC3bp10bt3b91RyEq4uLggODgYy5cvR05ODmrUqKE7ElGVMWaT0DoAnxddNgHIAPCFEfdrD+CwiGSIyDUAywA8ccs6owDMEZHzACAiZ4wNTmTNrl69ijVr1iA4OBguLi6645AVCQsLw5UrV7Bu3TrdUYiqlDGbhPxF5MGiiw8Ki8huIx67MYDjJa6fKFpW0gMAHlBK7VRKpSql+hgbnMiarVu3DleuXEF4eLjuKGRlevTogUaNGnGzENkdY0ZYbiIi3wHoYKLndwLgA6AHgDAAcUqpureupJSKUUqlKaXSzp49a6KnJtInKSkJ7u7u6N69u+4oZGUcHR3x9NNPY926dbh44/xCRHag3MKilJpQ4vJ/RTvGnjTisX8D0KTEdY+iZSWdALBWRPJF5CiAdBQWmJuISKyIBIpIYMOGDY14aiLLdfHiRaSkpODpp5+Go6Oj7jhkhcLCwpCXl4fVq1eXvzKRjTBmhMWtxMUFhfu03LovSln2APBRSnkrpaoBGAJg7S3rrEHh6AqUUg1QuIkow6jkRFZq9erVyMvL49FBVGEdOnSAt7c3NwuRXSn3KCEA+0UkueQCpVQIgOTbrA8AEJHrSqnnAHwJwBHAAhH5RSk1HUCaiKwtuu0xpdR+AAUAJopIdkVeCJG1SEpKQrNmzdC+fXvdUchKKaUwZMgQvPXWWzhz5gzuvZcHWJLtM2aE5SUjl5UiIiki8oCI3F80nwtE5LWisgIpNEFEWhbt3LvM+OhE1ufMmTPYtGkThgwZUnwiUaKKCAsLQ0FBAVasWKE7ClGVuNPZmh8H0BdAY6XUByVuqg3gurmDEdmiFStWoKCgAEOGDNEdhaycv78/WrVqhaSkJIwdO1Z3HCKzu9MIy0kAaQByAewtcVkLgDNdEVVAUlISWrVqBX9/f91RyAaEhYVhx44dOHbsjvN4EtmEO52t+UcRWQTgfhFZVOKy6sZEb0RkvGPHjmHHjh3c2ZZMJjQ0FACwfPlyzUmIzO+2hUUp9WnRP79XSv1066WK8hHZjBt/VLg5iEzlr3/9K4KCgni0ENmFOx0l9Leir/2rIgiRrUtKSkL79u1x//33645CNiQsLAwTJkzAr7/+iubNm+uOQ2Q2d9ok9HvR16yyLlUXkcj6/frrr/j+++85ukImFxoaCqUUR1nI5t1pk9BlpdSlEpfLJb9WZUgia5eUlASlVPE+B0Smct9996F79+5ISkqCiOiOQ2Q2dxphcROR2iUubiW/VmVIImsmIli2bBm6d++O++67T3ccskFhYWFIT0/HDz/8oDsKkdkYdfJDpVSAUmq8Uup5pVRbc4cisiU//PADfv31Vx4dRGYTHBwMJycnbhYim2bMyQ9fA7AIQH0ADQDEK6WmmDsYka1ISkqCk5MTgoODdUchG1W/fn307t0by5Ytg8Fg0B2HyCyMGWGJABAkIq+LyOsAOgIYat5YRLbBYDBg2bJl6N27N+rXr687DtmwsLAwHD9+HLt27dIdhcgsjCksJwFUL3HdBcBv5olDZFt27tyJ48ePc3MQmd0TTzyBGjVqIDExUXcUIrMwprBcBPCLUipeKbUQwM8ALiilPrjlHENEdIukpCTUqFEDTzzxhO4oZONq1aqFAQMGIDk5Gfn5+brjEJmcMYVlNYCXAWwBsBXAKwA+w//OLUREZcjPz0dycjIGDBiAWrVq6Y5DdiAsLAznzp3Dpk2bdEchMrk7zXQLACg6nxAR3aVNmzbh3Llz3BxEVebxxx9HnTp1kJSUhD59+uiOQ2RSxhwl1F8p9b1S6g9OHEdkvKSkJNSpUwePP/647ihkJ1xcXDBo0CCsXr0aOTk5uuMQmZQxm4TeBzAcQH1OHEdknJycHKxevRrBwcFwcXHRHYfsSFhYGC5fvoyUlBTdUYhMypjCchzAz8I5n4mMlpKSgsuXL3NzEFW5nj174t577+UkcmRzyt2HBcAkAClKqW0A8m4sFJHZZktFZOWSkpLwl7/8BT179tQdheyMk5MTnn76acTFxeHSpUuoXZsD4mQbjBlhmQngKgrnYnErcSGiMly6dAmff/45nn76aTg6OuqOQ3YoLCwMeXl5WLNmje4oRCZjzAjLfSLiZ/YkRDZi9erVyMvL4+Yg0qZTp07w9PREYmIihg0bpjsOkUkYM8KSopR6zOxJiGxEUlISvLy80LFjR91RyE4ppRAWFoaNGzfizJkzuuMQmYQxhWUMgPVKqRwe1kx0Z6dPn8bGjRsRHh4OpZTuOGTHwsPDUVBQgOTkZN1RiEyi3MJSdBizg4jU4GHNRHeWnJyMgoIChIeH645Cds7f3x9+fn48txDZDGNGWKCUukcp1V4p1e3GxdzBiKxRYmIiHnzwQbRq1Up3FCKEh4dj165dOHr0qO4oRJVmzEy3zwDYDuBLANOKvk41bywi65ORkYHdu3dzdIUsxpAhQwAAy5Yt05yEqPKMGWH5G4AgAFki0hNAWwAXzJqKyArd+KNw448EkW7e3t7o3LkzNwuRTTCmsOSKSC4AKKVcROQggObmjUVkXUQECQkJeOihh+Dp6ak7DlGx8PBw/Pzzz9i3b5/uKESVYkxhOaGUqgtgDYANSqnPAGSZNxaRddm3bx/279+PiIgI3VGIbnJjAkOOspC1M+YooadE5IKITAXwKoD5AJ40dzAia5KQkAAnJycMHjxYdxSimzRs2BCPPfYYkpKSYDAYdMchqjCjjhK6QUS2ichaEblmrkBE1sZgMCApKQm9e/dGgwYNdMchKiU8PBxZWVnYvXu37ihEFXZXhYWIStu5cyeOHz/Oo4PIYj3xxBOoUaMGNwuRVWNhIaqkxMREuLq6YuDAgbqjEJXJzc0NAwcOxKeffor8/HzdcYgqxNiJ4zyVUo8U/buGUopnayYCcO3aNXz66ad44oknUKtWLd1xiG4rPDwc586dw4YNG3RHIaoQYyaOGwVgBYD/FC3yQOERQ0R2b/369fjjjz94dBBZvD59+qBevXpYunSp7ihEFWLMCMs4AF0AXAIAETkE4F5zhiKyFgkJCWjQoAEee4wnNCfLVq1aNTz99NNYs2YNLl++rDsO0V0zprDklTwqSCnlBEDMF4nIOly6dAlr165FaGgonJ2ddcchKldkZCRycnKwZg0Hycn6GFNYtimlXgZQQyn1KIBkAP81bywiy7dq1Srk5uYiMjJSdxQio3Tu3BleXl5ISEjQHYXorhlTWF4EcBbAPgCjAaQAmGLOUETWYOnSpbj//vvRoUMH3VGIjKKUQnh4ODZs2IBTp07pjkN0V4yZ6dYgInEiEgIgBsA3IsJNQmTXTp48ic2bNyM8PBxKKd1xiIwWEREBg8GA5cuX645CdFeMOUpoq1KqtlKqHoC9AOKUUu+ZPxqR5Vq2bBlEhEcHkdVp2bIl2rZty81CZHWM2SRUR0QuARgEYLGIdADQy7yxiCzb0qVLERgYiObNeeJysj6RkZHYs2cP0tPTdUchMpoxhcVJKeUO4GkAn5s5D5HF279/P77//nvubEtWa8iQIVBKcZSFrIoxhWU6gC8BHBaRPUqpZgAOmTcWkeVKSEiAg4MDQkNDdUchqpD77rsPvXr1wtKlS8FdEslaGLPTbbKIPCgiY4uuZ4hIsPmjEVkeg8GAhIQEPPLII2jUqJHuOEQVFhERgYyMDKSmpuqOQmQUp9vdoJT6EHeYIE5ExpslEZEF27FjB7KysjBjxgzdUYgqJTg4GGPHjsWSJUvQqVMn3XGIynXbwgIgrcpSEFmJxYsXo2bNmnjqqad0RyGqFDc3Nzz11FNYtmwZ3nvvPbi4uOiORHRHty0sIrKoKoMQWbqcnBwkJydj8ODBqFmzpu44RJU2dOhQJCYmIiUlhSWcLN6dRlgAAEqpLShj05CIPGyWREQWau3atbh06RKGDh2qOwqRSdzYF2vx4sUsLGTxyi0sAP6vxL+rAwgGcN08cYgs15IlS+Dh4YEePXrojkJkEk5OTggPD8eHH36I7Oxs1K9fX3ckotsy5iihvSUuO0VkAoAe5o9GZDlOnz6N9evXIyIiAo6OjrrjEJnMsGHDkJ+fz6n6yeIZMzV/vRKXBkqp3gDqVEE2IouxbNkyFBQUcHMQ2ZzWrVvD398fS5Ys0R2F6I6MmThuLwqPGNoLYDeAfwAYac5QRJZm8eLFCAgIQKtWrXRHITK5oUOHIjU1lVP1k0UzZpOQt4g0K/rqIyKPiciOqghHZAn279+P775rjqNHt8DBAfDyAjijOdmChITC9/Pkyf8HIBOTJ/+oOxLRbRmzSai6UmqCUmqVUmqlUurvSqnqVRGOyBK8+OJPAOJw/nxtiABZWUBMDEsLWbeEhML3cVYWIKIAeOKzz/ph6VKD7mhEZVLlnUdCKfUpgMsAlhYtCgdQV0RCzJytTIGBgZKWxjnt6PZMeRCPSAG2b/8dgEep21xcgI4dTfdcW7ea7rHINpnyvZ2aCuTllV7u7JyLzp1N95mU72sqj1Jqr4gElreeMYc1+4lIyxLXtyil9lc8GpH1OH9+M4BeZd5W1i97Imtxu/dvfj5nvCXLZExh+U4p1VFEUgFAKdUBRk7br5TqA+DfABwBzBORWbdZLxjACgBBIsLhE6oUU36ii4iIxy+/NIfB0LTUbZ6e/PRIVcuU7zcvr8LNQbdS6jjWravP2ZzJ4hhzlFA7ALuUUplKqUwUHikUpJTap5T66XZ3Uko5ApgD4HEALQGEKaValrGeG4C/AfimAvmJzObixYtYtWoVHn54E1xdb77N1RWYOVNPLiJTmDkTpd7XLi4FEHkRq1at0hOK6A6MKSx9AHgD6F508S5a1h/AgDvcrz2AwyKSISLXACwD8EQZ6/0TwL8A5N5FbiKz+/TTT5Gbm4s33vBDbGzhiIpShV9jY4GICN0JiSouIgKl3tfz5zugWbNvEB8frzseUSnlbhISkTIGDY3SGMDxEtdPAOhQcgWlVACAJiKyTik1sYLPQ2QW8fHxaNmyJQIDAxEUxIJCtici4tb3tUJGRhRee+01ZGVlwdPTU1c0olKMGWExC6WUA4DZKJyIrrx1Y5RSaUqptLNnz5o/HNm99PR07Nq1C1FRUVBK6Y5DVGWGDRsGoHCyRCJLYs7C8huAJiWuexQtu8ENgB+ArUX7xnQEsFYpVerQJhGJFZFAEQls2LChGSMTFVq0aBEcHBwQGRmpOwpRlfL09MTDDz+M+Ph4GAyck4UshzkLyx4APkopb6VUNQBDAKy9caOIXBSRBiLiJSJeAFIBDORRQqRbQUEBFi9ejD59+sDd3V13HKIqFxUVhYyMDOzYwUnNyXKYrbCIyHUAzwH4EsABAJ+KyC9KqelKqYHmel6iytq8eTNOnDiBqKgo3VGItBg0aBBq1arFnW/JopQ7062l4Uy3ZG4RERH44osvcPLkSVSvzrNQkH0aOXIkli9fjlOnTqFWrVq645ANM3amW2073RJZohtzr4SHh7OskF2Ljo7Gn3/+iZUrV+qOQgSAhYXoJklJScjNzeXmILJ7Xbp0gY+PDxYsWKA7ChEAFhaim8ybNw+tW7dGu3btdEch0kophREjRmD79u1IT0/XHYeIhYXohh9++AF79+7FyJEjOfcKEYDhw4fD0dGRoyxkEVhYiIrMnz8fLi4uiOCUtkQAAHd3d/Tr1w/x8fHIz8/XHYfsHAsLEYDc3FwkJCRg0KBBqFevnu44RBZj5MiROH36NFJSUnRHITvHwkIEYPXq1Th//jxGjhypOwqRRenbty/c3d0xf/583VHIzrGwEKFwZ1tvb2/07NlTdxQii+Lk5IThw4dj3bp1OHnypO44ZMdYWMjuZWRkYPPmzRgxYgQcHPgjQXSrESNGwGAwYNGiRbqjkB3jb2eyewsWLICDgwPnXiG6DR8fH3Tv3h3z58/nCRFJGxYWsmsFBQWIj49H79694eHhoTsOkcUaOXIkjhw5gu3bt+uOQnaKhYXs2vr16/Hbb79xZ1uicgQHB6NOnTqIi4vTHYXsFAsL2bW5c+eiUaNGGDiQJxAnuhNXV1cMHToUK1aswLlz53THITvEwkJ26/jx40hJScGIESPg7OysOw6RxRs9ejSuXbvGnW9JCxYWslvz5s2DiGDUqFG6oxBZBT8/P3Tu3BmxsbEQEd1xyM6wsJBdun79OubNm4c+ffrAy8tLdxwiqzF69Gikp6dj69atuqOQnWFhIbv0+eef4+TJkxg9erTuKERWJSQkBPfccw/mzp2rOwrZGRYWskv/+c9/0LhxY/Tr1093FCKrUqNGDQwfPhyrV6/GmTNndMchO8LCQnbn6NGj+PLLL/HMM8/AyclJdxwiqxMTE4P8/HwsXLhQdxSyIywsZHfi4uKglMIzzzyjOwqRVfL19UX37t0RGxvLmW+pyrCwkF25du0a5s+fj/79+3NmW6JKGD16NDIyMrBx40bdUchOsLCQXfnss89w5swZ7mxLVEmDBg1CgwYNuPMtVRkWFrIrc+bMgbe3N3r37q07CpFVc3FxwciRI7F27VocP35cdxyyAywsZDd++uknbNu2DWPHjoWjo6PuOERW79lnn4WIcJSFqgQLC9mNOXPmoHr16hgxYoTuKEQ2wcvLCwMGDEBsbCxyc3N1xyEbx8JCduH8+fNYunQpIiIiUK9ePd1xiGzGc889h3PnziE5OVl3FLJxLCxkF+Lj43H16lWMGzdOdxQim9KrVy80b94cH330ke4oZONYWMjmGQwGzJkzB126dEHbtm11xyGyKUopPPfcc/j222/x7bff6o5DNoyFhWzel19+iSNHjuC5557THYXIJg0bNgy1atXCnDlzdEchG8bCQjbvo48+QqNGjTBo0CDdUYhsUu3atTF8+HAsW7YMZ8+e1R2HbBQLC9m0w4cP44svvsDo0aNRrVo13XGIbNa4ceNw7do1zJs3T3cUslEsLGTT5syZA0dHR8TExOiOQmTTfH198cgjj+Djjz9Gfn6+7jhkg1hYyGZdunQJ8+fPR0hICO677z7dcYhs3vjx43HixAmsXLlSdxSyQSwsZLPmz5+Py5cvY8KECbqjENmFfv36wcfHB7Nnz4aI6I5DNoaFhWzS9evX8e9//xtdu3ZFYGCg7jhEdsHBwQEvvPAC9uzZg127dumOQzaGhYVs0urVq5GVlYUXXnhBdxQiuzJs2DDcc889mD17tu4oZGNYWMgmvffee2jWrBkGDhyoOwqRXalZsyaeffZZrFmzBhkZGbrjkA1hYSGbs3v3buzevRt///vfeVZmIg2ee+45ODo64oMPPtAdhWwICwvZjoQEwMsLHTp3RpZSGFWzpu5ERHbpvvvuw7/bt8eEDz6AODgAXl6FP59EleCkOwCRSSQkADExwNWrcADQVAR4/nnAxQWIiNCdjsi+JCQgZu9eON44Uigrq/DnE+DPI1WYsrZDzwIDAyUtLU13DKqsHj1M+3ipqUBeXunlLi5Ax46me56tW033WESWgj+PpJFSaq+IlHs4JzcJkW0o65fjnZYTkfnw55HMgJuESA9TfzLy8iocdr6Vpyc/hRGVhz+PZAU4wkI2Ie/113H11oWursDMmTriENm3mTMLf/5KuO7iwp9HqhQWFrIJcX/+iWcA5P7lL4BShZ/kYmO5gx+RDhERhT9/np4QpXDC0REzmjblzyNVCgsLWb38/Hy88847ONalC6qfOgUYDEBmJn85EukUEQFkZkIZDFjz/vuYdugQvv76a92pyIqxsJDVW758ObKysvDiiy/qjkJEZRgxYgQaNmyIWbNm6Y5CVoyFhayawWDArFmz4Ofnh379+umOQ0RlcHV1xd/+9jekpKTgp59+0h2HrBQLC1m1devW4ZdffsGLL74IpZTuOER0G2PHjkWtWrU4ykIVxsJCVktE8Oabb8LLywuhoaG64xDRHdxzzz149tlnsXz5chw5ckR3HLJCLCxktTZt2oTdu3dj0qRJcHLilEJElm7ChAmoVq0a3vyx+lsAABc9SURBVHjjDd1RyAqxsJBVEhFMnToVTZo0wYgRI3THISIjuLu7Y/To0Vi0aBEyMjJ0xyErw8JCVmnTpk3YuXMnXn75Zbi4uOiOQ0RGmjx5MpydnTGTk8jRXWJhIatTcnQlOjpadxwiugscZaGKYmEhq8PRFSLrxlEWqggWFrIqHF0hsn4cZaGKMGthUUr1UUr9qpQ6rJQqNQ2pUmqCUmq/UuonpdQmpZSnOfOQ9ePoCpFt4CgL3S2zFRallCOAOQAeB9ASQJhSquUtq30PIFBEHgSwAsBb5spD1k9E8Prrr3N0hcgGcJSF7pY5R1jaAzgsIhkicg3AMgBPlFxBRLaIyNWiq6kAPMyYh6xcSkoKdu3axdEVIhtxY5Tl9ddf1x2FrIA5C0tjAMdLXD9RtOx2RgL4wox5yIoVFBTgpZdewl//+leMHDlSdxwiMgF3d3f87W9/Q0JCAs8xROWyiJ1ulVKRAAIBvH2b22OUUmlKqbSzZ89WbTiyCImJidi3bx9mzJgBZ2dn3XGIyEQmT56MOnXq4KWXXtIdhSycOQvLbwCalLjuUbTsJkqpRwC8AmCgiOSV9UAiEisigSIS2LBhQ7OEJcuVl5eHV199FQEBAQgJCdEdh4hM6J577sFLL72ElJQUbN++XXccsmDmLCx7APgopbyVUtUADAGwtuQKSv1/e/ceHkV973H8/Q0BTAwVRWIEJBEwx/JECgWEltKHYoUWoXiAFhQDCB4sGCpesFy0B3ugsUXyeHnkEkBpbVoFFW+VClrFoASxgqQaNNISCggBeopBLmnY3/ljtx7kIgF29zebfF7/7O5kmP0wT/Kbz87MzlgnYD7hslIZwyySwObNm0dFRQX5+fkkJQVip6CIRFFeXh4tWrTgpz/9Kc4533EkoGI2+jvnaoA84GWgDFjinHvfzH5uZj+IzDYLSAOWmtkGM3v+JIuTeurTTz9lxowZ9O7dm6uvvtp3HBGJgdTUVKZPn05JSQnPPfec7zgSUJZobbZLly7unXfe8R1D4mT69Once++9vP3223Tt2tV3HBGJkZqaGnJyckhKSqK0tJQGDRr4jiRxYmZ/ds51OdV82r8ugbVr1y5mz57NkCFDVFZE6rjk5GRmzpxJWVkZixcv9h1HAkiFRQJr2rRpHDp0SFfCFKknBg0aRPfu3bn77rupqqryHUcCRoVFAmn9+vU8+uij/OQnPyE7O9t3HBGJAzPjwQcfZOfOnfziF7/wHUcCRoVFAsc5x6233kqzZs245557fMcRkTi68sorGTFiBAUFBbpkv3yBCosEztKlSykuLmbGjBk0bdrUdxwRibP8/HwaNmzInXfe6TuKBIgKiwTKwYMHmTRpEh06dOCmm27yHUdEPGjRogVTpkxh2bJlvPbaa77jSECosEigzJ49m61bt/LAAw/oa40i9djtt99OZmYmEydO5MiRI77jSACosEhgbN++nfz8fAYNGsR3vvMd33FExKOUlBTuv/9+Nm7cyIIFC3zHkQBQYZHAmDhxIqFQiFmzTngPTBGpZwYPHkyvXr2YOnUqlZW6e0t9p8IigfCHP/yBp556irvvvps2bdr4jiMiAWBmzJ07l/3793PHHXf4jiOeqbCId5999hm33HILX/3qV5k0aZLvOCISIJdffjmTJ0/mt7/9La+++qrvOOKRCot4d++991JRUcH8+fNp1KiR7zgiEjBTp06lXbt2jBs3jkOHDvmOI56osIhXGzdupKCggDFjxtCzZ0/fcUQkgM455xzmzp1LeXk5+fn5vuOIJ7pbs3gTCoXo0aMHmzdvZtOmTVxwwQW+I4lIgN1www0sWbKEjRs3cvnll/uOI1GiuzVL4M2dO5eSkhIKCgpUVkTklAoKCkhLS2Ps2LGEQiHfcSTOVFjEi82bN3PXXXfRp08fhg8f7juOiCSA9PR0Zs+eTXFxMQ8//LDvOBJnKiwSd0eOHGHUqFE0bNiQRYsWYWa+I4lIghg1ahT9+/dn8uTJfPjhh77jSBypsEjcPfDAA6xevZqHHnqIVq1a+Y4jIgnEzCgsLCQlJYWRI0dSU1PjO5LEiQqLxFVZWRnTpk1j4MCB5Obm+o4jIgno4osvZs6cOaxdu5b777/fdxyJExUWiZuamhpGjhxJWloa8+fP16EgETljQ4cOZciQIfzsZz+jtLTUdxyJAxUWiZv8/HzWrVvH3Llzueiii3zHEZEEZmbMmTOH888/nxEjRnD48GHfkSTGVFgkLoqLi5k+fTrXX389P/zhD33HEZE6oHnz5ixcuJANGzboth71gAqLxNyePXu47rrraNOmDfPmzfMdR0TqkAEDBjBx4kQefvhhli1b5juOxJAKi8RUKBRixIgR7N69myVLltCkSRPfkUSkjvnlL39Jly5dGD16NFu2bPEdR2JEhUViqqCggOXLl1NQUECnTp0oKoKsLEhKCj8WFflOKCKJ5thxZOnSRjz55JOEQiGGDRtGdXW174gSAyosEjMlJSVMmTKFwYMHM378eIqKYOxYqKgA58KPY8eqtIhI7Z1sHFmzpg2LFi1i7dq1TJs2zXdMiQHd/FBOqlevM/+31dW7ePfdLpgl07nzepKTm1JSAic6kb9xY+je/czf61ivvx69ZYnI2TmbceRETjWOlJffwo4dc2jffinNmw854/fROBI/uvmheBMKHeb99wfxr3/tpX37Z0hObgqceJD5sukiIsc61TjStm0BX/nKN9i0aST792+IXzCJOe1hkahyzjFmzBgee+wxlixZ8oWvMGdlhXffHiszE3SenIjURm3GkZ07d9K1a1eSkpJYt24d6enp8Ywop0l7WMSLhx56iMcee4x77rnnuOutzJwJqalfnD81NTxdRKQ2ajOOZGRk8Oyzz1JZWcmQIUN0Em4docIiUbNixQpuv/12rr32WqZPn37cz4cPh8LC8Cchs/BjYWF4uohIbdR2HOncuTOPPvooxcXF5OXlkWhHE+R4OiQkUVFaWkrPnj1p3bo1b731Fmlpab4jiYgwdepU8vPz+dWvfqWr4QZUbQ8JJccjjNRtW7ZsoW/fvqSlpfHCCy+orIhIYMyYMYPNmzdz1113kZ6ezsiRI31HkjOkwiJnZffu3fTt25eDBw9SXFxMZmam70giIp9LSkriN7/5DXv37mXMmDE0a9aM/v37+44lZ0DnsMgZq6qqol+/fmzdupUXX3yRnJwc35FERI7TuHFjli1bRseOHfnRj37EW2+95TuSnAEVFjkjhw4dYvDgwaxfv56lS5fSo0cP35FERE6qSZMmvPTSS7Rq1YprrrmG0tJS35HkNKmwyGk7ePAgAwcOZOXKlSxcuFC7V0UkIaSnp7NixQrOPfdcevfuzXvvvec7kpwGFRY5LQcOHGDAgAGsXLmSRYsWMWrUKN+RRERqLSsri9dff51zzjmH3r17s379et+RpJZUWKTW9u/fzzXXXMNrr73G4sWLGT16tO9IIiKnrV27dqxatYq0tDR69+6NLpWRGFRYpFb27dtHv379eOONN3j88ccZMWKE70giImesTZs2rFq1iqZNm/Ld736XNWvW+I4kp6DCIqe0bds2evbsyZo1a/jd737H9ddf7zuSiMhZy8rKYtWqVVx44YVcddVVPPvss74jyZdQYZHjFRWF7zCWlER1ixbcd8UVbNmyheXLlzN06FDf6UREoubfV+fu0KEDgwYN4o8jRnw+/pGVFR4PJRB0aX75oqIiGDsWDhz4fNIBM/bMnEnrKVM8BhMRiZ0DBw4w79vf5uY//5lzj/5BaqpuehZjtb00vwpLXdCrV/SWVVIChw8fP71xY+jePXrv8/rr0VuWiNRP0Rz7AFdSgmn8i7vaFhYdEpIvcCf6Y4UTlxgRkTrkhGUFNP4FhO4lVBdEqa2XlZXxlQ4daFlTc/wPMzP1qUBEgiXaY1JWFlRUHDd5f7NmpGn88057WATnHIsXL6Zr1678T0oKRxo3/uIMqakwc6afcCIi8TJzZni8O8qhpCTG7t3LqFGj2Ldvn6dgAios9d6OHTsYMGAAN954I507d+ZnmzbRYNGi8B4Vs/CjTjgTkfpg+PDweHfU+Ndo8WLa3XMPjz/+OFdccQUrV670nbLe0km39ZRzjqKiIiZMmMDhw4fJz89nwoQJJCWpw4qIHGvt2rWMHDmSDz/8kJtvvplZs2bRpEkT37HqBJ10Kyf10Ucf0b9/f3Jzc2nfvj0bNmzg1ltvVVkRETmJbt26sX79eu68804KCwvJycnh6aefJtE+9CcybaHqkX379jFp0iRycnIoLi6moKCAN954g+zsbN/RREQCLyUlhVmzZlFcXMx5553HkCFDuOqqq9i4caPvaPWCCks9UF1dTWFhIdnZ2cyePZvc3FzKy8u57bbbaNCgge94IiIJpUePHrz77rvMmTOH9957j06dOjFu3Dh27NjhO1qdpsJSh1VXVzN//nyys7O5+eabadeuHW+//TaLFi3ioosu8h1PRCRhJScnM27cOMrLy8nLy2PBggW0adOGCRMmsG3bNt/x6iQVljqoqqqKRx55hHbt2vHjH/+YjIwMXnrpJVavXk2XLqc8r0lERGrpggsu4MEHH+Sjjz4iNzeXefPm0bZtW8aPH095ebnveHWKCksd8sEHH5CXl0fLli3Jy8vjkksu4eWXX2bNmjV8//vfx8x8RxQRqZPatGnDggULKC8vZ9SoUSxcuJDs7Gz69u3L888/z5EjR3xHTHj6WnOC27t3L8888wxFRUWsWrWKRo0aMXToUMaPH0+3bt1UUkREPPjkk09YuHAh8+fPZ/v27bRu3Zrc3FyGDRtGTk6O73iBopsf1mGVlZUsX76cJ598kpUrV1JTU8Nll13GmDFjGD16NM2bN/cdUUREgJqaGl544QXmzZvHK6+8QigUon379gwbNoxrr72WnJycev/BMhCFxcy+BzwINAAWOufuO+bnjYHfAJ2BvcBQ59yWL1tmfSws+/fvp6SkhJUrV7JixQo2bNgAQGZmJsOGDWPo0KF07Nix3v/Si4gE2a5du3j66ad54oknWL16Nc45MjIyuPrqq+nTpw+9evWiZcuW9W4s915YzKwB8BFwNbANWAdc55z74Kh5xgMdnHM/NrNhwH8654Z+2XITrrAUFcG0abB1K7RuHb5XxUkuc++cY9euXZSXl1NaWsq6detYt24dZWVlhEIhGjZsSI8ePejTpw99+vRh06avM22a1WbRIiLi2dGbgxYtaujXbzVVVfN55ZVX2LNnDwAZGRl07dqVK6+8kk6dOpGdnU1WVhYNGzas/cITbIMQhMLyDWC6c65v5PUUAOdc/lHzvByZZ42ZJQM7gebuS0IlVGEpKsKNHYsdOPD5pJpGjSgZM4YN7dtTWVnJrl27qKyspKKigo8//piqqqrP523evPnnv7jdunXjW9/6Fmlpaf9eNGPHwlGLJjVVt/0REQmiLxuzr7suxIYNG3jzzTc//6C6adOmz+dLTk4mKyuLtm3bkpGRQXp6Ounp6Vx44YWkpKTQdu1aOs6ZQ/Lhw8cvPAE2CEEoLEOA7znnboq8zgW6OefyjprnL5F5tkVeb47Ms+dky41pYenVK7rLKymBo3+BIrYAlxL+itb51pALrBHNkxqRmZRK66QUWielcGmDVDKs8Ul3Dfb79Ak+cRnHTW/cGLp3j0583U1dROqrOG0OuNh28tJXhh03fb+r4eMjn7E1dJC/hw5SETrI9tAh/hGq5h+ummr+f9v9NyDrRG+aIBuE2haW5JgliCIzGwuMBWjdurXnNKfhRL+dQCbwapNvcJ41pMEZHqvc6dJP5y1FRMSjk43NJxvL0yyZjsnn0ZHzjvuZc44DHOF/3b847EJk7j/Jh/g6tkGIZWHZDlxy1OtWkWknmmdb5JDQeYRPvv0C51whUAjhPSwxSQvRb5BZWVBRcdxky8yk95a3zmrRrU+8aDIztWdERORsxWlzQOvMJDpuOcs3O9nC69gGIZYXjlsHXGZml5pZI2AY8Pwx8zwPjIw8HwL86cvOX0k4M2eGjyMeLTU1PD24ixYRkSiL6ZhdTzYIMSsszrkaIA94GSgDljjn3jezn5vZDyKzLQKamdnHwO3A5Fjl8WL48PBJT5mZYBZ+jNJJUDFctIiIRFlMx+x6skHQheNERETEm9qedKt7CYmIiEjgqbCIiIhI4KmwiIiISOCpsIiIiEjgqbCIiIhI4KmwiIiISOCpsIiIiEjgqbCIiIhI4KmwiIiISOCpsIiIiEjgqbCIiIhI4KmwiIiISOAl3M0PzWw3UOE7xxm4ENjjO0QdoPUYHVqP0aH1GB1aj9GTiOsy0znX/FQzJVxhSVRm9k5t7kYpX07rMTq0HqND6zE6tB6jpy6vSx0SEhERkcBTYREREZHAU2GJn0LfAeoIrcfo0HqMDq3H6NB6jJ46uy51DouIiIgEnvawiIiISOCpsMSYmX3PzD40s4/NbLLvPInIzC4xs9fM7AMze9/MbvWdKZGZWQMzW29mL/rOksjMrKmZPWVmm8yszMy+4TtTIjKz2yJ/138xs9+b2Tm+MyUCM3vUzCrN7C9HTbvAzFaaWXnk8XyfGaNNhSWGzKwB8AjwfaA9cJ2ZtfebKiHVAHc459oD3YFbtB7Pyq1Ame8QdcCDwB+dc5cDX0Pr9LSZWUvgJ0AX51wO0AAY5jdVwlgMfO+YaZOBV51zlwGvRl7XGSossXUl8LFz7q/OuWrgCWCg50wJxzn3iXPu3cjzKsIbhpZ+UyUmM2sFXAMs9J0lkZnZecC3gUUAzrlq59w//aZKWMlAipklA6nADs95EoJz7g3gH8dMHgj8OvL818C1cQ0VYyossdUS+PtRr7ehDe1ZMbMsoBOw1m+ShPUAcBcQ8h0kwV0K7AYeixxeW2hm5/oOlWicc9uB+4GtwCfAPufcCr+pEtpFzrlPIs93Ahf5DBNtKiySMMwsDXgamOic+9R3nkRjZv2BSufcn31nqQOSga8Dc51znYDPqGO73+Mhco7FQMIFsAVwrpnd4DdV3eDCXwGuU18DVmGJre3AJUe9bhWZJqfJzBoSLitFzrlnfOdJUD2AH5jZFsKHJ3ub2W/9RkpY24Btzrl/7+l7inCBkdPzXeBvzrndzrl/Ac8A3/ScKZHtMrOLASKPlZ7zRJUKS2ytAy4zs0vNrBHhk8me95wp4ZiZET5XoMw5V+A7T6Jyzk1xzrVyzmUR/l38k3NOn2bPgHNuJ/B3M/uPyKSrgA88RkpUW4HuZpYa+Tu/Cp28fDaeB0ZGno8EnvOYJeqSfQeoy5xzNWaWB7xM+Oz3R51z73uOlYh6ALlAqZltiEyb6px7yWMmkQlAUeTDyF+BGz3nSTjOubVm9hTwLuFvA66nDl+pNZrM7PdAL+BCM9sG/DdwH7DEzMYAFcCP/CWMPl3pVkRERAJPh4REREQk8FRYREREJPBUWERERCTwVFhEREQk8FRYREREJPBUWEQkbiJ3OB4fed4i8pVWEZFT0teaRSRuIveCejFyZ14RkVrTheNEJJ7uA9pGLgBYDnzVOZdjZqMI31n2XOAywjfEa0T4goGHgX7OuX+YWVvgEaA5cAD4L+fcpvj/N0Qk3nRISETiaTKw2TnXEZh0zM9ygEFAV2AmcCByY8E1wIjIPIXABOdcZ+BOYE5cUouId9rDIiJB8ZpzrgqoMrN9wAuR6aVAh8jdur8JLA3fdgaAxvGPKSI+qLCISFAcPup56KjXIcJjVRLwz8jeGRGpZ3RISETiqQpocib/0Dn3KfA3M/shhO/ibWZfi2Y4EQkuFRYRiRvn3F7gTTP7CzDrDBYxHBhjZu8B7wMDo5lPRIJLX2sWERGRwNMeFhEREQk8FRYREREJPBUWERERCTwVFhEREQk8FRYREREJPBUWERERCTwVFhEREQk8FRYREREJvP8DTijK3OSfdOQAAAAASUVORK5CYII=\n",
"text/plain": [
"<Figure size 648x504 with 1 Axes>"
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
}
],
"source": [
"plot_sampling_points(func, N=10, peg=\"gridpoints\")"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The red points indicate the sampling that QuTiP would use when propagating the guess control with `mesolve`. The blue points indicate the sampling that Krotov will use internally: Each iteration yields an update for the blue points, and once the optimization finishes, the blue points are transformed back to the sampling of the red points.\n",
"\n",
"An important detail is that if the guess controls are given as an (analytical) function, not an array, Krotov's [`optimize_pulses`](https://krotov.readthedocs.io/en/latest/API/krotov.optimize.html#krotov.optimize.optimize_pulses) will sample that function exactly onto `tgrid` before converting to the \"guess pulses\" defined on the midpoints. This is what was visualized above. If the controls were not discretized, but we were to choose the \"pulses\" based on the analytical formula, we would end up with the following sampling:"
]
},
{
"cell_type": "code",
"execution_count": 7,
"metadata": {
"attributes": {
"classes": [],
"id": "",
"n": "6"
}
},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAiwAAAGtCAYAAAA1cy8JAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAIABJREFUeJzs3XlclOXeBvDrBgVEUMSlcAPtoKCiCGiWqRiaG27IokLmipqe0krLN3Mr3+qYbcY5RamkkMqmkHlSyS1LU1xSUzM33BURTRRkmfv9Y4AXFHVYhnuW6/v5zGecZ5555ppxgN88v/u5HyGlBBEREZEhs1AdgIiIiOhxWLAQERGRwWPBQkRERAaPBQsREREZPBYsREREZPBYsBAREZHBY8FCREREBo8FCxERERk8FixERERk8GqoDlBeDRo0kC4uLqpjEBERURXYt2/fdSllw8etZ3QFi4uLC1JTU1XHICIioioghEjTZT22hIiIiMjgsWAhIiIig8eChYiIiAye0Y1hISIi9fLy8nDhwgXk5OSojkJGwsbGBk2bNkXNmjUr9HgWLEREVG4XLlyAvb09XFxcIIRQHYcMnJQSGRkZuHDhAlq0aFGhbbAlRERE5ZaTk4P69euzWCGdCCFQv379Su2RY8FCREQVwmKFyqOynxcWLERERGTwWLAQERGVsG3bNvz666+qY9B9WLAQERGVwILFMLFgISIioxQdHY3OnTvD09MTEydORFpaGlxdXXH9+nVoNBp069YNmzZtAgAMGTIE3t7eaNu2LSIjI4u38eOPP8LLywsdOnSAn58fzp49iy+//BKffPIJPD098fPPP6t6eXQfHtZMRESVMm3aNBw8eLBKt+np6YlPP/30ofcfO3YMa9aswS+//IKaNWvi5Zdfxvbt2/Hmm29i8uTJ6Ny5M9q0aYMXXngBALBs2TI4OjoiOzsbnTp1wrBhw6DRaDBhwgTs2LEDLVq0wI0bN+Do6IhJkybBzs4Ob7zxRpW+JqocvRUsQohlAPwBXJNStivjfgHgMwD9AdwFMFpKuV9feYiIyHT89NNP2LdvHzp16gQAyM7ORqNGjTBv3jzExcXhyy+/LFVEff7551i7di0A4Pz58/jrr7+Qnp6O7t27F88L4ujoWP0vhHSmzz0sUQC+ALDiIff3A+BaeHkawH8Kr4nIjMXEAG+/DZw7BzRvDixcCISGqk5Fj/KoPSH6IqXESy+9hPfff7/U8rt37+LChQsAgKysLNjb22Pbtm1ISUnBrl27YGtrC19fX87Qa4T0NoZFSrkDwI1HrDIYwAqptRuAgxDCSV95iMjwxcRIhIdLpKUBUgJpaUB4uERMjOpkZGj8/PwQHx+Pa9euAQBu3LiBtLQ0vPnmmwgNDcWCBQswYcIEAMCtW7dQr1492Nra4vjx49i9ezcAoEuXLtixYwfOnDlTvA0AsLe3x+3btxW8KnoUlWNYmgA4X+L2hcJll9XEIaLy8vWt/DZyc9Nx82YKbtzYhKtX3wXQtNT9d+8KjBp1FXPnfod69V6ArW2bSk9AtW1bpR5OBqBNmzZ477338MILL0Cj0aBmzZr4+OOPsXfvXvzyyy+wtLREQkICli9fjpEjR+LLL7+Eu7s7WrdujS5dugAAGjZsiMjISAQEBECj0aBRo0bYvHkzBg4ciMDAQCQlJWHJkiXo1q2b4ldLACCklPrbuBAuANY/ZAzLegAfSCl3Ft7+CcCbUsrUMtYNBxAOAM2bN/dOS0vTW2Yi0l1lCpa//96DM2f+BzdvbgEgUaOGI/Lz01H2jl8NAEsAgLW1M5yd38aTT46BEBX7zsWCpfKOHTsGd3d31THIyJT1uRFC7JNS+jzusSr3sFwE0KzE7aaFyx4gpYwEEAkAPj4++quwiKhcKvKH/9SpU3j77bexffuawkGSc9GvXz94e3vjqacsUNb3EWdnC+zYkYbNmzdj2bJl+PXXcFhafoIPP/wQ/v7+nCKeyAyonIclGcAoodUFwC0pJdtBRCYqNzcXM2bMgLu7O77//nvMmTMHJ0+exNy5c9G5c2dYWlpi4ULA1rb042xttQNvmzdvjnHjxmHnzp1Yu3YtCgoKMGjQIPTs2ROnTp1S86KIqNrorWARQqwCsAtAayHEBSHEOCHEJCHEpMJVNgA4DeAkgK8BvKyvLESk1rVr1+Dn54ePPvoIL730Ek6ePIn58+fD3t6+1HqhoUBkJODsDAihvY6MLH2UkBACQ4YMwZEjR/Cf//wHhw4dQqdOnZCSklLNr4qIqpNex7Dog4+Pj0xNfWCYCxEZqAMHDmDw4MG4fv06li1bhuHDh1fp9k+fPo1Bgwbh+PHjWLx4MV555RW2iKoBx7BQRVRmDAun5icivVmzZg26du0KKSV27txZ5cUKALRs2RK7du3CwIEDMW3aNIwdO5ZzbBCZIBYsRKQXUVFRGD58OLy8vJCamgovLy+9PZe9vT0SEhIwZ84cREVFISgoCHl5eXp7PiKqfixYiKjKrVu3DuPGjUOvXr3w008/4YknntD7c1pYWGD+/Pn497//jfXr12PMmDHQaDR6f14yD76+vijPcIRt27bB39+/ynOMHz8eR48efeQ669ate+w6VSEqKgpTp07V+/MUYcFCRFXqp59+QkhICDp16oS1a9fC2tq6Wp9/8uTJWLhwIWJiYvDqq6/C2MbpmaqYGMDFBbCw0F5z9uKK+eabb9CmTZtHrlORgiU/P78ysaoFCxYiqjJ79uzB4MGD0apVK2zYsAF2dnZKcsyaNQuvv/46vvjiC8ybN09JBvp/MTFAeDjuO+VC5YqWO3fuYMCAAejQoQPatWuHNWvWAAAWLFiATp06oV27dggPDy8uWH19fTF9+nT4+PjA3d0de/fuRUBAAFxdXTF79mwAwNmzZ+Hm5obQ0FC4u7sjMDAQd+/efeC5N23ahGeeeQZeXl4ICgpCVlYWAODHH3+Em5sbvLy8kJiYWGbunJwcjBkzBh4eHujYsSO2bt0KQLu3IiAgAH379oWrqytmzpxZ5uNL7umxs7PD22+/jQ4dOqBLly64evUqfv31VyQnJ2PGjBnw9PTEqVOncOrUKfTt2xfe3t7o1q0bjh8/DgAYPXo0Jk2ahKeffhozZ86Ei4sLbt68Wfxcrq6uuHr1Kr7//ns8/fTT6NixI3r16oWrV68+kCsuLg7t2rVDhw4d0L1798f/B1aElNKoLt7e3pKIDM9ff/0lHR0dZcuWLeXFixdVx5EajUaOHTtWApARERGq45ico0ePPnadHj20F2trKbWlSumLtbX2/oqIj4+X48ePL7598+ZNKaWUGRkZxcvCwsJkcnJyYZYecubMmVJKKT/99FPp5OQkL126JHNycmSTJk3k9evX5ZkzZyQAuXPnTimllGPGjJGLFi0qfvzevXtlenq67Natm8zKypJSSvnBBx/I+fPny+zsbNm0aVN54sQJqdFoZFBQkBwwYMADuT/66CM5ZswYKaWUx44dk82aNZPZ2dly+fLlskWLFvLmzZsyOztbNm/eXJ47d66M91SbQ0opARS/vhkzZsh3331XSinlSy+9JOPi4oof8/zzz8sTJ05IKaXcvXu37NmzZ/F6AwYMkPn5+VJKKV955RW5bNmy4vX8/PyklFLeuHFDajQaKaWUX3/9tXzttdeklFIuX75cTpkyRUopZbt27eSFCxeklFJmZmaW/Z8my/7cAEiVOvz95x4WIqq0nJwcBAUFAdB++2zcuLHiRNr5WiIjI+Hv749p06Zhz549qiOZrXv3yrdcFx4eHti8eTPefPNN/Pzzz6hbty4AYOvWrXj66afh4eGBLVu24I8//ih+zKBBg4of27ZtWzg5OcHa2hotW7bE+fPaU9s1a9YMXbt2BQCEhYVh586dpZ539+7dOHr0KLp27QpPT098++23SEtLw/Hjx9GiRQu4urpCCIGwsLAyc+/cubP4Pjc3Nzg7O+PEiRMAtCd0rFu3LmxsbNCmTRs87jQ0VlZWxeNkvL29cfbs2QfWycrKwq+//oqgoCB4enpi4sSJuHz5/+doDQoKgqWl9rQXISEhxXuqVq9ejZCQEADAhQsX0KdPH3h4eGDRokWl3tMiXbt2xejRo/H111+joKDgkbkrigULEVXatGnTcPDgQaxYsQJPPfWU6jjFLC0tsWLFCjRu3BjBwcHIzMxUHcmsbNumvTg7l32/s3PFz+vUqlUr7N+/Hx4eHpg9ezYWLFiAnJwcvPzyy4iPj8fhw4cxYcKEUoe4F42nsrCwKDW2ysLCongMx/1z+Nx/W0qJ3r174+DBgzh48CCOHj2KpUuXVuxF3KdkJktLy8eOK6lZs2Zxvoetr9Fo4ODgUJz34MGDOHbsWPH9tWvXLv73M888g5MnTyI9PR3r1q1DQEAAAOCf//wnpk6disOHD+Orr74qc9qAL7/8Eu+99x7Onz8Pb29vZGRklO/F64AFCxFVyqpVq/DVV19h5syZGDBggOo4D6hXrx5iY2Nx6dIljBkzhoNwFXjUKRcq6tKlS7C1tUVYWBhmzJiB/fv3F/8hbdCgAbKyshAfH1/u7Z47dw67du0CAHz33Xd47rnnSt3fpUsX/PLLLzh58iQA7ViaEydOwM3NDWfPni0+TcSqVavK3H63bt0QUzh458SJEzh37hxat25d7pyPYm9vj9u3bwMA6tSpgxYtWiAuLg6AtuD6/fffy3ycEAJDhw7Fa6+9Bnd3d9SvXx8AcOvWLTRp0gQA8O2335b52FOnTuHpp5/GggUL0LBhw+I9VlWJBQsRVdiff/6J8PBwdO3aFe+9957qOA/VuXNnLFq0CElJSfjkk09UxzE7upxyobwOHz6Mzp07w9PTE/Pnz8fs2bPh4OCACRMmoF27dujTpw86depU7u22bt0aERERcHd3R2ZmJiZPnlzq/oYNGyIqKgojRoxA+/bt8cwzz+D48eOwsbFBZGQkBgwYAC8vLzRq1KjM7b/88svQaDTw8PBASEgIoqKiqvxIuuHDh2PRokXo2LEjTp06hZiYGCxduhQdOnRA27ZtkZSU9NDHhoSEIDo6urgdBADz5s1DUFAQvL290aBBgzIfN2PGDHh4eKBdu3Z49tln0aFDhyp9TQCn5ieiCrp79y66dOmCy5cv48CBA2jatKnqSI8kpURgYCCSk5Px888/o0uXLqojGTVTnJr/7Nmz8Pf3x5EjR1RHMVmcmp+Iqt3s2bNx+PBhrFy50uCLFUC7u3vp0qVo3rw5QkNDyzxclYgMFwsWIiq3vXv34rPPPsPkyZPRt29f1XF05uDggGXLluH06dOcn4Ue4OLiwr0rBowFCxGVS15eHsaPH48nn3wS77//vuo45dajRw9MmDABixcvxv79+1XHISIdsWAhonL56KOPcOjQIfz73/8unvvC2PzrX/9Co0aNMH78eKOYkpyIWLAQUTn89ddfmD9/PoYNG4bBgwerjlNhDg4O+OKLL3DgwAEeNURkJFiwEJFOpJQIDw+HjY0NlixZojpOpQUEBGDw4MGYO3du8dwZRGS4WLAQkU6ioqKwbds2LFq0CE5OTqrjVJoQAhEREahZsyYmTZrECeXokUqedFAX27ZtK542vyp9+umnpY5w69+/f6kTFla3OXPmICUl5YHl+nj9LFiI6LFu376NWbNm4dlnn8W4ceNUx6kyTZo0wXvvvYeUlBSsX79edRzTFhMDuLgAFhba68qcqtmM3V+wbNiwAQ4ODkqyFBQUYMGCBejVq1e1PB8LFiJ6rH/961+4evUqPv74Y1hYmNavjUmTJqF169aYMWMG8vLyVMcxTTExQHg4kJamPVFzWpr2diWKljt37mDAgAHo0KED2rVrV3zSvgULFqBTp05o164dwsPDi/ec+fr6Yvr06fDx8YG7uzv27t2LgIAAuLq6Yvbs2QC0E8e5ubkhNDQU7u7uCAwMLHO+nk2bNuGZZ56Bl5cXgoKCkJWVBQD48ccf4ebmBi8vLyQmJpaZOycnB2PGjIGHhwc6duyIrVu3AtDuwQwICEDfvn3h6uqKmTNnPvDYzz//HJcuXULPnj3Rs2dPANpDsa9fv16cffTo0WjVqhVCQ0ORkpKCrl27wtXVtfjkn3fu3MHYsWPRuXNndOzYscxZbzUaDV5++WW4ubmhd+/e6N+/f/FpDlxcXPDmm2/Cy8sLcXFxGD16dPF9urz+StHllM6GdPH29n7oaauJqOqdP39e1qpVSw4fPlx1FL1JSkqSAGRERITqKEbj6NGjj1+pRw/txdpaSm2pUvpiba29vwLi4+Pl+PHji2/fvHlTSillRkZG8bKwsDCZnJxcGKWHnDlzppRSyk8//VQ6OTnJS5cuyZycHNmkSRN5/fp1eebMGQlA7ty5U0op5ZgxY+SiRYuKH793716Znp4uu3XrJrOysqSUUn7wwQdy/vz5Mjs7WzZt2lSeOHFCajQaGRQUJAcMGPBA7o8++kiOGTNGSinlsWPHZLNmzWR2drZcvny5bNGihbx586bMzs6WzZs3l+fOnXvg8c7OzjI9Pf2B22fOnJGWlpby0KFDsqCgQHp5eckxY8ZIjUYj161bJwcPHiyllHLWrFly5cqVUkopMzMzpaura/FrKRIXFyf79esnCwoK5OXLl6WDg4OMi4srfr4PP/yweN2XXnpJxsXF6fz6y/rcAEiVOvz9N62vSkRU5WbPng2NRmOUc67oauDAgejRowfmzp2LW7duqY5jeu7dK99yHXh4eGDz5s1488038fPPPxcfYr9161Y8/fTT8PDwwJYtW/DHH38UP2bQoEHFj23bti2cnJxgbW2Nli1bFp+sr1mzZujatSsAICwsDDt37iz1vLt378bRo0fRtWtXeHp64ttvv0VaWhqOHz+OFi1awNXVFUIIhIWFlZl7586dxfe5ubnB2dkZJ06cAAD4+fmhbt26sLGxQZs2bZCWllau96RFixbw8PCAhYUF2rZtCz8/Pwgh4OHhgbNnzwLQ7h364IMP4OnpCV9fX+Tk5ODcuXMPZAwKCoKFhQWefPLJ4r05RUqeZ6iIrq+/MmpU+RaJyGQcOHAAK1aswIwZM+Di4qI6jt4IIbB48WL4+Pjggw8+MOnirFpt26a9dnHRtoHu5+z8/+uUU6tWrbB//35s2LABs2fPhp+fH2bOnImXX34ZqampaNasGebNm1d8BmcAxScZtLCwKHXCQQsLi+L5eIQQpZ7n/ttSSvTu3fuBszEfPHiwQq+jpJKZLC0tyz1H0P2vqeTrLdqWlBIJCQmVOkN07dq1K/zYyuAeFiIqk5QSr7/+OhwdHTFr1qzSdxrrAMpH5Pb29saLL76ITz75pNzfbOkxFi4EbG1LL7O11S6voEuXLsHW1hZhYWGYMWMG9u/fX1ycNGjQAFlZWcVjK8rj3Llz2LVrFwDgu+++w3PPPVfq/i5duuCXX37ByZMnAWjHhJw4cQJubm44e/Zs8SHy9xc0Rbp164aYws/diRMncO7cuXIVD/b29rh9+3a5X1eRPn36YMmSJcVjew4cOPDAOl27dkVCQgI0Gg2uXr2KbToUlbq+/spgwUJEZVq/fj22bt2KefPmlT4KQQ8DKKuFDrkXLlwIIQT+53/+R2FQExQaCkRGaveoCKG9jozULq+gw4cPo3PnzvD09MT8+fMxe/ZsODg4YMKECWjXrh369OmDTp06lXu7rVu3RkREBNzd3ZGZmYnJkyeXur9hw4aIiorCiBEj0L59ezzzzDM4fvw4bGxsEBkZiQEDBsDLywuNGjUqc/svv/wyNBoNPDw8EBISgqioqFJ7Rh4nPDwcffv2faBNo6t33nkHeXl5aN++Pdq2bYt33nnngXWGDRuGpk2bok2bNggLC4OXl9djZ7XW9fVXhiiqsoyFj4+PLM+x8EQEwNe3XKtrpET71FTkSYkjPj6oWfLIoN27yx57YG0NdOlSuZz6pGPut8+cwf+eO4eD3t7oYGdX/uepYIvD2Bw7dgzu7u6qY1Sps2fPwt/fnydABJCVlQU7OztkZGSgc+fO+OWXX/Dkk09WertlfW6EEPuklD6Peyz3sBDRAxKuX8cfd+9inotL6WIF0MsAymqhY+4ZzZqhrqUlFrAtRGbM398fnp6e6NatG955550qKVYqi3tYiKgUjUaDDh06ID8/H0eOHIGlpWXpFR41gLLwSASDVI7cc+fOxYIFC/D777+jffv21RLP2JjiHhbSP+5hIaIqk5iYiCNHjuCdd955sFgB9DKAslqUI/e0adNQp04dLFiwoJrCGSdj+8JLalX288KChYiKaTQazJ8/H61bty5zrgUAehlAWS3KkbtevXqYNm0aEhIScOjQIQVhDZ+NjQ0yMjJYtJBOpJTIyMiAjY1NhbfBlhARFYuPj0dQUBBiYmIwcuRI1XGUyszMhIuLC3r37l2hw2NNXV5eHi5cuFBqnhOiR7GxsUHTpk1Rs2bNUst1bQmxYCEiANq9K56ensjNzcUff/xRdjvIzMyZMwfvvvsuDh06BA8PD9VxiEwSx7AQUbmsXbsWhw8ffvjYFTPEsSxEhoMFCxFBSol3330XrVq1wvDhw1XHMRiOjo549dVXER8fX+qcNERU/ViwEBE2b96M33//HW+99Rb3rtzn1Vdfha2tLRYvXqw6CpFZY8FCRFi0aBGcnJzMfqBtWerXr4+xY8ciOjoaly5dUh2HyGyxYCEycwcPHkRKSgpeffXVcp3TxJxMnz4dBQUFWLJkieooRGaLBQuRmVu8eDHs7OwwceJE1VEMVsuWLTFs2DD85z//qdSZcomo4liwEJmx8+fPY/Xq1Rg/fnzpMzLTA9544w3cunULS5cuVR2FyCyxYCEyY5999hmklJg2bZrqKAavc+fO6N69Oz755BPk5eWpjkNkdliwEJmpW7duITIyEsHBwXB2dlYdxyi88cYbOHfuHOLi4lRHITI7LFiIzFRkZCRu376NN954Q3UUozFgwAC4ubnho48+4jl0iKoZCxYiM5SXl4fPPvsMzz//PLy8vFTHMRoWFhZ4/fXXceDAAWzdulV1HCKzwoKFyAytXbsWFy9exPTp01VHMTphYWFo0KABD3EmqmYsWIjMUEREBFq2bIl+/fqpjmJ0bGxsMH78eCQnJ+PcuXOq4xCZDRYsRGbm0KFD2LFjByZPnsxp+Cto0qRJAIAvv/xScRIi88GChcjMREREwMbGBmPHjlUdxWg5Oztj0KBB+Prrr5GTk6M6DpFZYMFCZEZu3ryJ6OhohIaGwtHRUXUcozZlyhRcv34dsbGxqqMQmQUWLERmJCoqCnfv3sWUKVNURzF6fn5+aN26NSIiIlRHITILLFiIzIRGo0FERASeffZZdOzYUXUcoyeEwNSpU7Fnzx7s2bNHdRwik8eChchMbNq0CSdPnsTUqVNVRzEZo0aNgp2dHfeyEFUDFixEZiIiIgJPPPEEhg0bpjqKyahTpw5GjRqFNWvWID09XXUcIpPGgoXIDJw5cwY//PADJk6cCCsrK9VxTMqUKVNw7949nsWZSM9YsBCZgaVLl0IIgfHjx6uOYnLatGmD7t2745tvvoFGo1Edh8hksWAhMnH5+flYvnw5+vXrh2bNmqmOY5ImTJiAU6dOYdu2baqjEJksFixEJiwmBmjcOBeXLp3Hnj2xiIlRncg0DRs2DLVqjcPAgR6wsABcXMD3mqiK1VAdgIj0IyYGCA8H7t61BQCkp9siPFx7X2iowmAmKDGxFvLy/o38fO34oLQ08L0mqmJCSqk6Q7n4+PjI1NRU1TGIqpyvb9Vub/du4N69B5dbWwNdulTNcxhrB8QY32vAeN9vokcRQuyTUvo8bj22hIhMVFl/QB+1nCqO7zWR/rElRGQgqvrbs7OzxLlzoozl/KZe1a/fxUXbBrof32uiqqPXPSxCiL5CiD+FECeFEG+VcX9zIcRWIcQBIcQhIUR/feYhMichIYcA3Cm1zNYWWLhQTR5TtnCh9r0tie81UdXSW8EihLAEEAGgH4A2AEYIIdrct9psALFSyo4AhgP4t77yEJmbU6cWwN7+dTRvLiGE9tt+ZCQHgepDaKj2vXV2BgAJIc7h00/v8L0mqkL63MPSGcBJKeVpKWUugNUABt+3jgRQp/DfdQFc0mMeIrNx9epVJCcnY+JEe6SlCWg0wNmzLFb0KTRU+x6npu6HlM7IzY1SHYnIpOizYGkC4HyJ2xcKl5U0D0CYEOICgA0A/qnHPERmY+XKlcjPz8e4ceNURzE73t7e8PT0xLJly1RHITIpqo8SGgEgSkrZFEB/ACuFEA9kEkKECyFShRCpPMEY0aNJKREVFYVnnnkGbm5uquOYpTFjxmD//v04dOiQ6ihEJkOfBctFACXnAW9auKykcQBiAUBKuQuADYAG929IShkppfSRUvo0bNhQT3GJTMO+ffvwxx9/YPTo0aqjmK2RI0eiZs2a+Pbbb1VHITIZ+ixY9gJwFUK0EEJYQTuoNvm+dc4B8AMAIYQ7tAULd6EQVUJUVBRsbGwQHBysOorZatCgAfz9/REdHY28vDzVcYhMgt4KFillPoCpADYCOAbt0UB/CCEWCCEGFa72OoAJQojfAawCMFoa29S7RAbk3r17+O677zB06FA4ODiojmPWRo8ejWvXrmHjxo2qoxCZBL1OHCel3ADtYNqSy+aU+PdRAF31mYHInHz//ffIzMxkO8gA9OvXDw0bNkRUVBT8/f1VxyEyeqoH3RJRFYqKikKTJk3g5+enOorZq1mzJsLCwpCcnIyMjAzVcYiMHgsWIhNx5coV/Pjjj3jxxRdhaWmpOg5B2xbKy8vDqlWrVEchMnosWIhMRExMDAoKCvDSSy+pjkKF2rdvj44dOyIqKkp1FCKjx4KFyAQUzb3SpUsXzr1iYEaPHo19+/bhyJEjqqMQGTUWLEQmYP/+/Thy5AgH2xogzslCVDVYsBCZgBUrVsDa2ppzrxigBg0aYMCAAYiOjkZBQYHqOERGiwULkZErGtQ5cOBA1KtXT3UcKsOLL76IK1euYMuWLaqjEBktFixERi4lJQXp6ekICwtTHYUeon///nBwcEB0dLTqKERGiwULkZGLjo6Go6Mj+vWvCEarAAAgAElEQVTrpzoKPYSNjQ2CgoKQmJiIO3fuqI5DZJRYsBAZsdu3b2Pt2rUIDg6GlZWV6jj0CGFhYcjKykJy8v2nVCMiXbBgITJi69atQ3Z2NkJDQ1VHocd47rnn0KxZM7aFiCqIBQuREYuOjoaLiwueffZZ1VHoMSwsLBAaGoqNGzfi2rVrquMQGR0WLERG6sqVK0hJSUFoaCgsLPijbAzCwsJQUFCANWvWqI5CZHT4W47ISK1evRoajYbtICPStm1beHp6si1EVAEsWIiMVHR0NLy9veHu7q46CpVDWFgY9uzZgxMnTqiOQmRUWLAQGaFjx45h3759nHvFCI0YMQJCCMTExKiOQmRUWLAQGaGYmBhYWFhg+PDhqqNQOTVu3Bh+fn6Ijo6GlFJ1HCKjwYKFyMhIKbFq1Sr4+fnhySefVB2HKmDkyJE4ffo09uzZozoKkdFgwUJkZPbs2YPTp09jxIgRqqNQBQ0dOhRWVlZYtWqV6ihERoMFC5GRWb16NaysrDB06FDVUaiCHBwc0L9/f8TGxvIMzkQ6YsFCZESK5vAoOpkeGa8RI0bg8uXL2LFjh+ooREaBBQuREdmxYwcuX77MdpAJ8Pf3h52dHdtCRDpiwUJkRFatWgU7Ozv4+/urjkKVZGtri8GDByM+Ph65ubmq4xAZPBYsREYiNzcX8fHxGDx4MGxtbVXHoSowYsQIZGZmYtOmTaqjEBk8FixERmLTpk3IzMxkO8iE9O7dG46OjmwLEemABQuRkVi1ahUcHR3Ru3dv1VGoilhZWSEwMBBJSUm4e/eu6jhEBo0FC5ERuHv3LpKSkjBs2DBYWVmpjkNVaPjw4bhz5w6+//571VGIDBoLFiIjsH79ety5c4ftIBPUvXt3ODk5YfXq1aqjEBk0FixERmDVqlVwcnJC9+7dVUehKmZpaYmQkBBs2LABN2/eVB2HyGCxYCEycH///Tf++9//Ijg4GJaWlqrjkB4MHz4cubm5SEpKUh2FyGCxYCEycElJSbh37x5CQkJURyE96dy5M1xcXLBmzRrVUYgMFgsWIgMXGxuL5s2bo0uXLqqjkJ4IIRAcHIzNmzfjxo0bquMQGSQWLEQGLDMzExs3bkRQUBCEEKrjkB4FBwcjPz8fa9euVR2FyCCxYCEyYOvWrUNeXh7bQWbAy8sLTz31FNtCRA/BgoXIgK1ZswYtWrSAj4+P6iikZ0IIhISEYMuWLUhPT1cdh8jgsGAhMlAZGRlISUlBcHAw20FmIjg4GAUFBUhMTFQdhcjgsGAhMlCJiYkoKChgO8iMtG/fHq1bt2ZbiKgMLFiIDNSaNWvg6uoKT09P1VGomhS1hbZv344rV66ojkNkUFiwEBmga9euYevWrWwHmaHg4GBoNBokJCSojkJkUFiwEBmghIQEaDQatoPMUNu2bdG2bVu2hYjuw4KFyACtWbMG7u7uaNeuneoopEBISAh27tyJixcvqo5CZDBYsBAZmCtXrmDHjh1sB5mx4OBgSCnZFiIqgQULkYFJSEiAlBJBQUGqo5AirVu3hoeHB2JjY1VHITIYLFiIDExcXBzc3d3Rtm1b1VFIoaCgIPzyyy9sCxEVYsFCZEBKtoPIvBXtYWNbiEiLBQuRAUlMTGQ7iAAAbm5u8PDwQFxcnOooRAaBBQuRAYmNjWU7iIoVtYUuXbqkOgqRcixYiAxEUTuIe1eoSFBQEI8WIirEgoXIQLAdRPdzc3NDu3bt2BYiAgsWIoMRFxcHNzc3toOolKCgIOzcuZNtITJ7LFiIDMCVK1ewfft2ThZHD2BbiEiLBQuRAWA7iB6m6BQNbAuRuWPBQmQA2A6iRylqC12+fFl1FCJlWLAQKXbr3/9G1LZtOHr8OESLFkBMjOpIZGDG1aqF01LiycaNARcXfkbILNVQHYDIrMXEwHb6dNQtup2WBoSHa/8dGqoqFRmSmBg0mTfv/2/zM0JmSkgpVWcoFx8fH5mamqo6BpkrX9+q3d7u3cC9ew8ut7YGunSpuufZtq3qtkWPZoyfEX4+SCEhxD4ppc/j1mNLiEghWdYfIqDsP1BknvgZIQLAlhBR+VTxN9Gs+vVhf+PGg3c4O/Nbr7Gq6v83FxdtG+h+/IyQmeEeFiKFIho3xt37512xtQUWLlQTiAzPwoXaz0QJslYtfkbI7LBgIVIkIyMDs48dQ7K/v/bbshDa68hIDqak/xcaqv1MODtDCoGzALaMGMHPCJmdxxYsQohWQoifhBBHCm+3F0LM1mXjQoi+Qog/hRAnhRBvPWSdYCHEUSHEH0KI78oXn8h4JSUloaCgAK5z5wJnzwIajfaaf4jofqGh2s9GQQFecHXFh+fPq05EVO102cPyNYBZAPIAQEp5CMDwxz1ICGEJIAJAPwBtAIwQQrS5bx3Xwm13lVK2BTCtXOmJjFh8fDxatGgBLy8v1VHISAghEBgYiC1btuD69euq4xBVK10KFlsp5Z77luXr8LjOAE5KKU9LKXMBrAYw+L51JgCIkFJmAoCU8poO2yUyepmZmUhJSUFgYCDPHUTlEhgYiIKCAiQlJamOQlStdClYrgshngIgAUAIEQhAl/mhmwAoud/yQuGykloBaCWE+EUIsVsI0VeH7RIZveTkZOTl5SEwMFB1FDIyHTt2RMuWLREfH686ClG10qVgmQLgKwBuQoiL0LZtJlfR89cA4ArAF8AIAF8LIRzuX0kIES6ESBVCpKanp1fRUxOpEx8fj+bNm6NTp06qo5CRKWoLpaSk4EZZh8QTmajHFiyFLZ1eABoCcJNSPielPKvDti8CaFbidtPCZSVdAJAspcyTUp4BcALaAub+DJFSSh8ppU/Dhg11eGoiw3Xr1i1s2rSJ7SCqsMDAQOTn57MtRGbloRPHCSFee8hyAICU8uPHbHsvAFchRAtoC5XhAEbet846aPesLBdCNIC2RXRap+RERmr9+vXIzc1lO4gqzMfHB87OzkhISMCYMWNUxyGqFo/aw2JfePGBtgXUpPAyCcBjD2uQUuYDmApgI4BjAGKllH8IIRYIIQYVrrYRQIYQ4iiArQBmSCkzKvpiiIxBfHw8mjZtiqefflp1FDJSRW2hTZs24datW6rjEFWLx578UAixA8AAKeXtwtv2AH6QUnavhnwP4MkPyZjdvn0bDRs2xKRJk/Dpp5+qjkNGbPfu3XjmmWewcuVKhIWFqY5DVGFVefLDJwDklridW7iMiMrphx9+wL1799gOokrr3LkzmjZtiri4ONVRiKqFLic/XAFgjxBibeHtIQC+1V8kItMVHx8PJycnPPvss6qjkJGzsLDAsGHD8OWXX+Lvv/9GnTp1VEci0itdjhJaCGAMgMzCyxgp5f/qOxiRqblz5w42bNiAgIAAWFjwNF5UeYGBgbh37x5++OEH1VGI9E6Xcwk1B3AdwNrCS0bhMiIqh//+97/Izs5mO4iqzLPPPgsnJydOIkdmQZeW0A8onOUWQC0ALQD8CaCtvkIRmaL4+Hg0bNgQ3bp1Ux2FTISFhQUCAgKwdOlSZGVlwc7OTnUkIr3RpSXkIaVsX3hxhfYcQbv0H43IdGRnZ2P9+vUICAiApaWl6jhkQoKCgpCTk4P//ve/qqMQ6VW5G+lSyv0AOIEEUTls3LgRd+7cYTuIqtxzzz2HRo0asS1EJu+xLaH7Zry1gHbSuEt6S0RkguLj41G/fn306NFDdRQyMZaWlggICMDKlStx9+5d2Nraqo5EpBe67GGxL3GxhnZMy2B9hiIyJTk5OUhOTsaQIUNQs2ZN1XHIBAUGBuLOnTv48ccfVUch0htdBt0elVKWmplICBEEgLMVEelg8+bNuH37NttBpDc9evRA/fr1ERcXh4CAANVxiPRClz0ss3RcRkRliI+Ph4ODA55//nnVUchE1ahRA0OHDsX69euRk5OjOg6RXjy0YBFC9BNCLAHQRAjxeYlLFID8aktIZMRyc3ORlJSEIUOGwMrKSnUcMmFBQUHIysrCpk2bVEch0otH7WG5BCAVQA6AfSUuyQD66D8akfFLSUnBrVu32A4ivevZsyfq1avHcwuRyXroGBYp5e8AfhdCxEgpuUeFqALi4+NRp04d9OrVS3UUMnE1a9bEkCFDkJCQgHv37sHa2lp1JKIq9aiWUGzhPw8IIQ7df6mmfERGKy8vD+vWrcPgwYP5x4OqRVBQEP7++29s3rxZdRSiKveoo4ReLbz2r44gRKZmy5YtyMzMZDuIqo2fnx/q1q2L+Ph4+PvzVzeZlke1hC4XXqdVXxwi0xEfHw97e3u88MILqqOQmbCyssLgwYORlJSE3NxcDvQmk/KoltBtIcTfJS63S15XZ0giY5OXl4e1a9di4MCBsLGxUR2HzEhgYCBu3ryJn376SXUUoir10IJFSmkvpaxT4mJf8ro6QxIZm+3btyMjI4PtIKp2L7zwAuzt7XluITI5Op38UAjhJYR4RQjxTyFER32HIjJ28fHxqF27Nvr27as6CpkZa2trDBo0COvWrUNeXp7qOERV5rEFixBiDoBvAdQH0ABAlBBitr6DERmr/Px8JCYmwt/fH7Vq1VIdh8xQYGAgbty4ga1bt6qOQlRldNnDEgqgk5RyrpRyLoAuAF7Ubywi4/Xzzz8jPT0dQUFBqqOQmerTpw/s7OzYFiKTokvBcglAyVGD1gAu6icOkfGLi4uDra0t20GkTK1ateDv74/ExETk53PeTzINuhQstwD8IYSIEkIsB3AEwM2icwvpNx6RcSkoKEBiYiIGDBiA2rVrq45DZiwoKAgZGRnYtm2b6ihEVeJRE8cVWVt4KbJNP1GIjN/PP/+Mq1evsh1EyvXr1w+1a9dGXFwcTw1BJuGxBYuU8tvqCEJkCmJjY1GrVi30799fdRQycyXbQhEREahRQ5fvp0SGS5ejhPyFEAeEEDc4cRzRwxW1g/z9/dkOIoMQHByM69evsy1EJkGXMSyfAngJQH1OHEf0cGwHkaEp2RYiMna6FCznARyRUkp9hyEyZnFxcWwHkUHh0UJkSnQpWGYC2CCEmCWEeK3oou9gRMakoKAACQkJPDqIDE5QUBCuX7+O7du3q45CVCm6FCwLAdyFdi4W+xIXIipU1A4KDg5WHYWoFLaFyFToMmy8sZSynd6TEBkxtoPIUNna2ha3hb744gseLURGS5c9LBuEEC/oPQmRkWI7iAxdUFAQ0tPT2RYio6ZLwTIZwI9CiGwe1kz0IB4dRIauX79+sLW1ZVuIjNpjC5bCw5gtpJS1eFgz0YOK2kEDBgxQHYWoTLa2thg4cCCPFiKjpsseFggh6gkhOgshuhdd9B2MyBgUFBQgPj6e7SAyeGwLkbHTZabb8QB2ANgIYH7h9Tz9xiIyDtu3b8e1a9cQEhKiOgrRI/Xv3x+1a9fGmjVrVEchqhBd9rC8CqATgDQpZU8AHQHc1GsqIiMRGxuL2rVr8+ggMni1atXCoEGDkJiYiLy8PNVxiMpNl4IlR0qZAwBCCGsp5XEArfUbi8jw5efnIyEhAQMHDoStra3qOESPFRISgoyMDGzZskV1FKJy06VguSCEcACwDsBmIUQSgDT9xiIyfFu3bsX169c5WRwZjT59+qBOnTqIjY1VHYWo3HQ5SmiolPKmlHIegHcALAUwRN/BiAzdmjVrYGdnh379+qmOQqQTGxsbDB48GImJicjNzVUdh6hcdDpKqIiUcruUMllKyU86mbXc3FwkJiZi8ODBsLGxUR2HSGchISG4efMmNm/erDoKUbmUq2AhIq2ffvoJmZmZPDqIjE7v3r3h4ODAthAZHRYsRBWwZs0a1K1bFy+8wLNWkHGxsrLC0KFDsW7dOuTk5KiOQ6QzXSeOcxZC9Cr8dy0hBM/WTGbr3r17WLduHYYMGQJra2vVcYjKLTg4GH///Tc2bdqkOgqRznSZOG4CgHgAXxUuagrtEUNEZmnTpk24desWjw4io+Xn5wdHR0dOIkdGRZc9LFMAdAXwNwBIKf8C0EifoYgMWWxsLOrVq4devXqpjkJUITVr1kRAQACSk5ORnZ2tOg6RTnQpWO6VPCpICFEDgNRfJCLDlZ2djXXr1mHo0KGwsrJSHYeowkJCQpCVlYUNGzaojkKkE10Klu1CiP8BUEsI0RtAHIDv9RuLyDBt2LABWVlZGDFihOooRJXSs2dPPPHEE1i9erXqKEQ60aVgeQtAOoDDACYC2ABgtj5DERmqVatWoVGjRvD19VUdhahSLC0tERQUhPXr1+Pvv/9WHYfosXSZ6VYjpfxaShkEIBzAb1JKtoTIrMTEAM2ba5CQEIvs7GNYs6aG6khElVa//j+Rk3MMdevaw8VF+zknMlSP/a0rhNgGYFDhuvsAXBNC/CqlnK7nbEQGISYGCA8H7t7V1ve3bzsiPFx7X2iowmBElRATAyxa5ApAAADS0sDPNRk08bidJUKIA1LKjkKI8QCaSSnnCiEOSSnbV0/E0nx8fGRqaqqKpyYjUdXdmt27gXv3HlxubQ106VJ1z7NtW9Vti0xTVX62+bkmQyGE2Cel9HncerqMYakhhHACEAxgfaWTERmZsn6pP2o5kTHg55qMjS6N+AUANgLYKaXcK4RoCeAv/cYiqriq/kbn4qLdXX4/Z2d+e6TqVZWfN36uydjoMug2TkrZXkr5cuHt01LKYfqPRmQYFi4ELCxKn3PF1la7nMhYLVyo/RyXVKuWhp9rMlgP3cMihFiCR0wQJ6V8RS+JiAxMz56XoNHMQN26Efj7bwc0b679Zc+BiWTMij6/b78NnDsnIWUahg07htDQfmqDET3Eo1pCHNlKBO1U/MB32L37Hbi5OaiOQ1RlQkOLCheBjh2H4q+/rAGwYCHD9NCCRUr5bXUGITJUq1evhqenJ9zc3FRHIdKbESNG4M0338Tp06fRsmVL1XGIHqDL2Zq3CiG23H+pjnBEqp06dQq//fYbhg8frjoKkV6FhIQAAKfqJ4Oly1FCb5T4tw2AYQDy9ROHyLB89913AMBzB5HJc3Z2xnPPPYeYmBjMmjULQgjVkYhK0eUooX0lLr9IKV8D4KvLxoUQfYUQfwohTgoh3nrEesOEEFII8diJY4iqi5QSMTEx6N69O5o3b646DpHehYaG4ujRo/j9999VRyF6gC4tIccSlwZCiD4A6urwOEsAEdCO4GoDYIQQok0Z69kDeBXAb+VOT6RH+/fvx59//olQHg5EZiIoKAg1atRADE8qRAZIl5lu90F7xNA+ALsAvA5gnA6P6wzgZOG8LbkAVgMYXMZ67wL4EEBOGfcRKRMTE4OaNWsiMDBQdRSialG/fn3069cPq1atQkFBgeo4RKXo0hJqIaVsWXjtKqV8QUq5U4dtNwFwvsTtC4XLigkhvKA9P9EPj9qQECJcCJEqhEhNT0/X4amJKqegoACrV69G//794ejoqDoOUbUJDQ3FxYsXsWPHDtVRiErRpSVkI4R4TQiRKIRIEEJME0LYVPaJhRAWAD6Gdo/NI0kpI6WUPlJKn4YNG1b2qYkea+vWrbh8+TLbQWR2Bg4cCDs7O7aFyODo0hJaAaAtgCUAvij890odHncRQLMSt5sWLitiD6AdgG1CiLMAugBI5sBbMgQxMTGwt7eHv7+/6ihE1crW1hYBAQGIj49HTg479WQ4dClY2kkpx0kptxZeJkBbtDzOXgCuQogWQggrAMMBJBfdKaW8JaVsIKV0kVK6ANgNYJCUkjPsklLZ2dlISEjAsGHDUKtWLdVxiKpdaGgobt26hQ0bNqiOQlRMl4JlvxCiS9ENIcTT0GHafillPoCp0J7p+RiAWCnlH0KIBUKIQRUNTKRv69evx+3bt9kOIrP1/PPP44knnmBbiAyKLhPHeQP4VQhxrvB2cwB/CiEOA5BSyvYPe6CUcgOADfctm/OQdX11SkykZzExMXByckLPnj1VRyFSokaNGggJCcGXX36JmzdvwsGB59Ai9XTZw9IXQAsAPQovLQqX+QMYqL9oRNXvxo0b2LBhA4YPHw5LS0vVcYiUCQ0NRW5uLhISElRHIQKg22HNaY+6VEdIouoSGxuLvLw8hIWFqY5CpFSnTp3QqlUrrFypyzEWRPqnyx4WIrOxYsUKtGvXDh07dlQdhUgpIQRGjRqF7du34+zZs6rjELFgISry119/YdeuXRg1ahRP/EYEFO9pjI6OVpyEiAULUbGVK1fCwsKCRwcRFXJ2doavry9WrFgBKaXqOGTmWLAQAdBoNFi5ciV69eqFxo0bq45DZDBGjRqFv/76C7t371YdhcwcCxYiADt37sTZs2cxatQo1VGIDErRBIorVqxQHYXMHAsWImgH29rZ2WHIkCGqoxAZlDp16mDo0KFYvXo17t27pzoOmTEWLGT2srOzERsbi8DAQNSuXVt1HCKDM2rUKNy8eRPr169XHYXMGAsWMntJSUm4ffs220FED+Hn5wcnJye2hUgpFixk9lasWIFmzZqhR48eqqMQGaQaNWogNDQUGzZsQHp6uuo4ZKZYsJBZu3LlCjZu3IgXX3wRFhb8cSB6mFGjRiE/Px+rVq1SHYXMFH9Dk1mLjo6GRqNhO4joMTw8PNCxY0d8++23qqOQmWLBQmZLSolly5bh2WefRevWrVXHITJ4Y8eOxf79+3Hw4EHVUcgMsWAhs/Xbb7/h2LFjGDt2rOooREZh5MiRsLKywvLly1VHITPEgoXM1rJly2Bra4vg4GDVUYiMgqOjI4YMGYLo6GjOyULVjgULmaU7d+5g9erVCA4Ohr29veo4REZj7NixuHHjBpKTk1VHITPDgoXMUkJCAm7fvs12EFE59erVC02bNsWyZctURyEzw4KFzNKyZcvwj3/8A88995zqKERGxdLSEqNHj8bGjRtx/vx51XHIjLBgIbNz6tQpbN++HWPGjIEQQnUcIqMzevRoSCk58y1VKxYsZHaioqJgYWHBuVeIKuipp56Cr68vli9fDiml6jhkJliwkFkpKChAVFQU+vTpg6ZNm6qOQ2S0xo4di1OnTuHnn39WHYXMBAsWMispKSm4cOECB9sSVdKwYcNgb2+PpUuXqo5CZoIFC5mVyMhINGjQAAMHDlQdhcio2draYuTIkYiLi8PNmzdVxyEzwIKFzMbly5eRnJyM0aNHw9raWnUcIqMXHh6O7OxsREdHq45CZoAFC5mNqKgo5OfnY8KECaqjEJkELy8veHt7IzIykoNvSe9YsJBZ0Gg0+Prrr+Hr64tWrVqpjkNkMsLDw3H48GH89ttvqqOQiWPBQmYhJSUFZ86cwcSJE1VHITIpI0aMgJ2dHb766ivVUcjEsWAhsxAZGYn69etj6NChqqMQmRR7e3uMHDkSa9as4eBb0isWLGTyrly5gqSkJA62JdKTosG3MTExqqOQCWPBQiavaLBteHi46ihEJsnb2xve3t746quvOPiW9IYFC5k0DrYlqh4cfEv6xoKFTEdMDODiAlhYaK9jYvDTTz/h9OnT3LtCpGdFg28jIyO1C8r4eSSqjBqqAxBViZgYIDwcuHtXezstDQgPx5E2bVC/fn0EBASozUdk4ooG365cuRKfd+kCu+nTH/h5BACEhqoLSUZNGFu/0cfHR6ampqqOQZXl61u129u9G7h374HFaQD+06wZPmjZsmqeZ9u2qtkOkSGpop/H37Oy4LlvH27VqIE6+fkPrmBtDXTpUiXPBYA/jyZCCLFPSunzuPXYEiLTUEaxAgDNAExq3Lh6sxCZqQ52dniuTh3YlVWsAA/9OSXSBVtCpEZVfzNycdHudr7P9Vq14LJ7d9U+F5GpqcKfxymrV+PciBFwKetOZ2fuFaEK4x4WMg0LFwK2tqUW3QFw9ZVX1OQhMlMBAQH4V506yLG478+Lra3255SogliwkGkIDQUiI7Xf4ITAZSsrzHnySbT93/9VnYzIrFhZWaHhtGkYp9Egr3FjQAjtz2VkJAfcUqWwYCHTERoKnD2LfXv3onFuLpxnzYLF/d/yiEjvJk6ciNgaNTBrxAhAowHOnmWxQpXG3+ZkciIiIlC7dm289NJLqqMQmaXGjRtj6NChWLZsGe4WHdpMVEksWMikZGRk4LvvvsOoUaNQt25d1XGIzNbUqVORmZmJVatWqY5CJoIFC5mUpUuX4t69e5gyZYrqKERmrVu3bmjXrh2++OILnl+IqgQLFjIZ+fn5iIiIgK+vL9q2bas6DpFZE0Jg6tSpOHjwIHbu3Kk6DpkAFixkMhITE3Hu3DlMnz5ddRQiAvDiiy+ifv36+Pjjj1VHIRPAgoVMgpQSixcvxj/+8Q/4+/urjkNEAGxtbTFp0iQkJSXh5MmTquOQkWPBQiZh165d2LNnD6ZPn85DmYkMyJQpU1CjRg189tlnqqOQkeNvdjIJH3/8MerVq8dDmYkMjJOTE0aOHIlly5YhMzNTdRwyYixYyOidOXMGa9euxcSJE1G7dm3VcYjoPtOnT8fdu3cRGRmpOgoZMRYsZPQ+//xzWFhYYOrUqaqjEFEZOnToAD8/PyxZsgR5eXmq45CRYsFCRu3WrVv45ptvMHz4cDRp0kR1HCJ6iNdeew0XL15EXFyc6ihkpFiwkFH75ptvkJWVxUOZiQxc37594ebmhsWLF3MiOaoQFixktHJzc/HZZ5+hR48e8PLyUh2HiB7BwsIC06dPx/79+7F9+3bVccgIsWAho/Xdd9/h/PnzmDlzpuooRKSDF198EY0aNcL777+vOgoZIRYsZJQ0Gg0+/PBDdOjQAf369VMdh4h0UKtWLUyfPh2bNm3Cvn37VMchI8OChYzSunXrcPz4ccyaNQtCCNVxiEhHkydPRp06dfDBBx+ojkJGhgULGR0pJd5//3384x//QGBgoOo4RFQOdevWxZQpU5CQkIA///xTdRwyIixYyOikpKQgNTUVM2fOhKWlpeo4RFRO06ZNg7W1Nf71r3+pjkJGhAULGZ3330f1LZUAABe2SURBVH8fjRs3xqhRo1RHIaIKaNSoEcaNG4eVK1fiwoULquOQkWDBQkblt99+w9atW/Haa6/B2tpadRwiqqA33ngDGo0GixcvVh2FjIReCxYhRF8hxJ9CiJNCiLfKuP81IcRRIcQhIcRPQghnfeYh4/f++++jXr16CA8PVx2FiCrBxcUFI0eORGRkJK5fv646DhkBvRUsQghLABEA+gFoA2CEEKLNfasdAOAjpWwPIB4AG5r0UIcOHUJSUhL++c9/wt7eXnUcIqqkt956C9nZ2fjkk09URyEjoM89LJ0BnJRSnpZS5gJYDWBwyRWklFullHcLb+4G0FSPecjIzZs3D3Xr1sW0adNURyGiKtCmTRsEBwfj888/514Weix9FixNAJwvcftC4bKHGQfgv2XdIYQIF0KkCiFS09PTqzAiGYv9+/dj7dq1mD59OurVq6c6DhFVkTlz5uDOnTv46KOPVEchA2cQg26FEGEAfAAsKut+KWWklNJHSunTsGHD6g1HBmHevHlwcHDg3hUiE9OmTRuMGDECS5YswbVr11THIQOmz4LlIoBmJW43LVxWihCiF4C3AQySUt7TYx4yUqmpqfj+++/x+uuvo27duqrjEFEVmzNnDnJycrBoUZnfWYkA6Ldg2QvAVQjRQghhBWA4gOSSKwghOgL4CtpihaU1lWnu3LlwdHTEK6+8ojoKEelB69atERoaioiICFy5ckV1HDJQeitYpJT5AKYC2AjgGIBYKeUfQogFQohBhastAmAHIE4IcVAIkfyQzZGZ+u2337Bhwwa88cYbqFOnjuo4RKQn77zzDnJzczn7LT2UkFKqzlAuPj4+MjU1VXUMqiZ9+/bFvn37cPr0aR7KTGTixowZg9WrV+P06dNwcnJSHYeqiRBin5TS53HrGcSgW6Ky7Ny5Exs3bsSMGTNYrBCZgXfeeQf5+fl47733VEchA8SChQySlBIzZsxA48aNMXXqVNVxiKgatGzZEhMmTEBkZOT/tXfv0VFUCR7HvzdPDFGQYYw8ExSGx7ARmKxEYYLKco6CPARBSFTYBZEVRUVFA5wVZSKguCwwyhjFBxAHRowKLCAigoOHOAQRBpNVYCSIAiIKSAIhJHf/6JYhBJRHd2515/c5x5Pu6rLyo0761q+r68HWrVtdxxGPUWERT8rNzSUvL48nn3ySuLg413FEpJo8/vjjxMbGMnbsWNdRxGNUWMRzysrKyMzMpE2bNgwePNh1HBGpRgkJCTzyyCMsXLiQjz/+2HUc8RAVFvGcl156ia1btzJlyhSioqJcxxGRavbQQw+dKC6hdmKIBI8Ki3jKjz/+yIQJE0hLS6NHjx6u44iIA/Hx8UyYMIG//vWvLFmyxHUc8QgVFvGUZ599lm+//Zann34aY4zrOCLiyNChQ/nNb37Do48+yvHjx13HEQ9QYRHP2LNnD1OnTqV///507NjRdRwRcSg6OprJkydTWFjIK6+84jqOeIAKi3hGZmYmx44dIysry3UUEfGAPn360KlTJ8aPH8/BgwddxxHHVFjEE/Ly8nj11Vd58MEHadGihes4IuIBxhhmzJjBvn37mDBhgus44pgKizhXUVHBfffdR4MGDRg/frzrOCLiIR06dOCuu+5i5syZFBQUuI4jDqmwiHOvvPIK+fn5PPPMM7oEv4hUkZWVxSWXXMKoUaN0mnMNpsIiTh04cIDMzEw6d+5Menq66zgi4kH169dn4sSJvP/+++Tm5rqOI46osIhTjz/+OPv372fmzJk6jVlEzujuu+8mOTmZ0aNHU1JS4jqOOKDCIs5s2bKF5557jrvvvpt27dq5jiMiHhYVFcXMmTPZuXMnkydPdh1HHFBhESfKy8sZPnw4devWZeLEia7jiEgISEtLIz09nSlTplBYWOg6jlQzFRapVjk5kJQEUVERrFv3Ov37v82vfvUr17FEJERMmzaN6OjBtG9/KRERlqQk37gi4U+FRapNTg4MHw5FRQAGSGLOnE4abETkrL333mWUlT1PaenlWGsoKvKNKxpHwp8JtVPEUlJSbH5+vusYNcJ11wV2eXl5UFpadXpsLKSmBu73rF4duGWJyIXROCK/xBizwVqb8kvzaQ+LVJvTDTI/N11E5FQaR2quKNcBxLsC/QmjSZNydu2KrDI9MVGfZkTCVaDf20lJP32tXJnGkfCnPSxSbRo3fh4orjQtLg50r0MROVtZWb5xo7ISHn1UN0cMdyosUi1ef/118vJG0a/fuyQmgjG+T0TZ2ZCR4TqdiISKjAzfuPHTONKgwTEiI/+TZcvu0GX7w5wKiwRdUVER99xzD9deey3z5/dixw6oqIAdO1RWROTcZWRwYhz55psYpk5tz+LFi3nxxRddR5MgUmGRoCovL2fw4MGUl5czd+5coqJ02JSIBNaoUaPo1q0bDz74IF988YXrOBIkKiwSVM8++yxr1qxh5syZXHHFFa7jiEgYioiI4NVXX6VWrVrcfvvtlJWVuY4kQaDCIkGzceNGxo8fT79+/Rg8eLDrOCISxho2bMgLL7zA+vXrefLJJ13HkSBQYZGgOHToEIMGDaJ+/fq88MILuhOziATdrbfeypAhQ3jqqadYtWqV6zgSYCosEnDWWoYOHcq2bdt4/fXXda8gEak2M2fOpGXLlgwcOJCvv/7adRwJIBUWCbhp06axcOFCJk2axHWBvi63iMjPiI+P58033+TIkSP079+fY8eOuY4kAaLCIgH14YcfMmbMGPr27cvDDz/sOo6I1ECtW7dm9uzZrFu3jkceecR1HAkQFRYJmN27d3PbbbdxxRVX8PLLL+u4FRFxZsCAATzwwAPMmDGD+fPnu44jAaDCIgFRWlrKgAEDOHToELm5udSpU8d1JBGp4Z5++mk6d+7MsGHD2Lx5s+s4coFUWOSCWWsZNmwYa9euZfbs2bRt29Z1JBERoqOjWbBgAXXr1qVHjx588803riPJBVBhkQv2xBNPMG/ePCZOnMjAgQNdxxEROaFhw4YsWbKEH374gZ49e3L48GHXkeQ8qbDIBZkzZw5PPPEEQ4YMYdy4ca7jiIhU0a5dOxYsWMCnn35Keno65eXlriPJeVBhkapyciApCSIifD9zck472+rVqxk2bBg33HCDLg4nIp7Wo0cPZsyYweLFixk9evSZZzzL8U+qn+5EJ5Xl5MDw4VBS4nteVOR7DpVurbx582ZuueUWmjdvzptvvklMTIyDsCIiZ2/kyJFs376dadOm0bhx46qnPJ/l+CduGGut6wznJCUlxebn57uO4S2BvDhbXh6UlladHhsLqakAFBYX02XTJmKMYW379iTVqnXuv2f16gvLKSJyHmNfubWkFxbyl337+GPz5oxs1OifL57F+BcQGv8qMcZssNam/NJ82sMilZ3uzXrS9G1HjtB182YigFVXXXV+ZUVExJFIY5jXqhVHKyq4d9s2LoqI4D8aNPC9+Avjn7ilPSxSWVKSbzfoqRITKVqzhrS0NIqLi1m9erVOXxaRkFVaWkqvXr147733mDdvHunp6T87/rFjR3VHrDHOdg+LDrqVyrKyIC6u8rS4OPaPHk3Xrl05ePAgK1asUFkRkZAWGxvLW2+9RVpaGnfeeSdvvPHGGcc/srLchJRKVFiksowMyM72faIwBhIT+XrCBNpPncq+fftYvnw5HTp0cJ1SROSCxcXFsXjxYlJTUxk4cCAvlpRUGf/IztYBtx6hY1ikqoyME2/Q/Px8brrpJiIiIlizZg3t2rVzHE5EJHAuvvhi3n33Xfr378/w4cPZP2kSj375pS7T4EHawyJn9MEHH3D99dcTHx/P2rVrVVZEJCzVrl2bd955h/T0dDIzMxkzZgyhdnxnTaA9LHJa8+fPZ8iQITRv3pwVK1bQsGFD15FERIImOjqauXPnUq9ePaZOncrevXvJzs6mls6E9AztYZFKysvLyczMZNCgQVx99dV8+OGHKisiUiNEREQwY8YM/vCHPzB37ly6dOmiGyZ6iAqLnHDw4EF69erF5MmTGTFiBCtXrqRevXquY4mIVBtjDOPGjeOtt96ioKCAlJQU8vLyXMcSVFjEr7CwkI4dO7JixQpmzZrFrFmzdLl9Eamx+vTpw7p167jooovo0qULs2fP1nEtjqmw1HDWWp5//nk6dOjA999/z6pVqxgxYoTrWCIizrVt25b169eTlpbGsGHDGDhwID/88IPrWDWWCksNtnfvXnr27MnIkSPp0qULmzdv5ve//73rWCIinlGvXj2WL1/OU089RW5uLsnJyazWvYCcUGGpod555x2Sk5NZuXIl06dPZ+nSpVx++eWuY4mIeE5kZCSZmZknviK64YYbGDNmDEeOHHEdrUZRYalhioqK6N27N3369CEhIYH8/HxGjRpFRIT+FEREfk5KSgqffPIJw4YN45lnnuG3v/0tS5cudR2rxtBWqoY4duwYU6ZMoU2bNqxcuZIpU6awYcMG3RNIROQcxMfHk52dzapVq4iNjaVHjx707duXr776ynW0sKfCEuYqKip44403SE5O5rHHHqNbt24UFBQwZswYoqOjXccTEQlJ119/PZs2bSIrK4tly5bRqlUrxo8fz4EDB1xHC1sqLGHKWsuyZctISUlhwIABREZGsmjRIt5++20SExNdxxMRCXkxMTGMHTuWgoICbr75ZrKysmjWrBmTJ0+muLjYdbywo8ISZsrLy8nNzaVz5850796dgwcPMmfOHDZv3kzPnj1dxxMRCTvNmjVjwYIFbNy4kU6dOpGZmcmVV17JpEmT2L9/v+t4YUOFJdhyciApCSIifD9zcoKy6KZNK8jI+F+aN29Ov3792L17N7NmzaKwsJA77riDyMjIgP1eERGpql27dixZsoSPPvqI5ORkxo4dS5MmTRgxYgSFhYXB3BwEdVvjFSbUrtyXkpJi8/PzXcc4Ozk5MHw4lJT8c1pcHGRnQ0ZGABZtKSk5+RboxbRs+SyTJv0LvXr1UkkREXFoy5YtTJ8+nblz51Ja2peIiJepqPjnzRQDtDkI6ramOhhjNlhrU35xPhWWk1x3XWCXl5cHpaVVJh8zsRRcknpeiyyzFWwqP8QDxUsopnGV12NjIfX8Fl2Fro0kIjVVIDcHx47t429/i6G8vE6V1xLMHpZfMvCClt/mUB4xtuq2JlQ2CGdbWKKClsAX4kZgOhAJvGStnXzK67HAHOB3wH7gNmvtjmBmqlanKSsA0af7wzoDay27Ko6yofwAH5V9T97xHzhMOXD6Oyif4VeKiIgjMTG/prz89K/ttZcx5PBGOkfVo2PUpbSMjCfGnNvRGmfcpoTZBiFoe1iMMZHAF0A3YBewHhhkrS04aZ57gGRr7QhjzEDgFmvtbT+33JD6SigpCYqKqk5PTIQdO6pMPn78OF9++SWff/45mzZtIi8vj7y8PL777jsAGjZsSPfu3enevTv339+br76q+kd9hkWLiIhDZ9oc1KlzgBYtuvHTdi0mJoYOHTqQmppKSkoKrVq1omXLlsTHx5/7wkNkg+CFPSxXA9ustf/wB5oP9AYKTpqnNzDB/3gh8EdjjLGh9j3VmWRlYYcPx5z0veLxmBjW3XQTn0yfzp49e9i7dy979uxhx44dbNu2jbKyshPztm7dmp49e5Kamso111xD27ZtMcZ3zEpJyem/sszKqrZ/nYiInKWsrNOP2c89V5eMjPXs3buXtWvXnvig+qc//YmjR4+emLdRo0a0aNGCBg0akJCQQEJCApdddhlxcXFcecsttJ81i6iT96iE4QYhmHtYbgVutNYO8z+/A+horb33pHm2+OfZ5X++3T/Pd2dabkjtYQGOv/Yau4YMoSmwExgL/Nn/WlRUFAkJCVx++eU0adKEli1bnmjTrVu3pm7duj+77JwcGDcOdu6Epk19f5shcHyViEiNdC5jdllZGZ9//nml/7Zv337ig+7hw4crzT8IeApIBExiYkhtEJwfdBvIwmKMGQ4MB2jatOnvik6368ujrLWsWrWKWrVqUatWLWJjY4mNjaVevXpceumluoePiIics+LiYvbt28eRI0c4evQopaWlHD16lJYtW9KgQQPX8c6JF74S+hpoctLzxv5pp5tnlzEmCqiD7+DbSqy12UA2+PawBCVtkBhj6Nq1q+sYIiISRmrXrk3t2rVdx6hWwfx4vx5oYYxpZoyJAQYCi06ZZxEw2P/4VmBV2By/IiIiIgETtD0s1trjxph7gXfxndb8srX2M2PMk0C+tXYRMBuYa4zZBnyPr9SIiIiIVBLU67BYa5cCS0+Z9l8nPT4K9A9mBhEREQl9OuJTREREPE+FRURERDxPhUVEREQ8T4VFREREPE+FRURERDxPhUVEREQ8T4VFREREPE+FRURERDxPhUVEREQ8T4VFREREPE+FRURERDxPhUVEREQ8z1hrXWc4J8aYfUCR6xznoT7wnesQYUDrMTC0HgND6zEwtB4DJxTXZaK19te/NFPIFZZQZYzJt9amuM4R6rQeA0PrMTC0HgND6zFwwnld6ishERER8TwVFhEREfE8FZbqk+06QJjQegwMrcfA0HoMDK3HwAnbdaljWERERMTztIdFREREPE+FJciMMTcaYz43xmwzxjzmOk8oMsY0McZ8YIwpMMZ8Zoy533WmUGaMiTTGbDTGLHGdJZQZY+oaYxYaY/7PGFNojLnGdaZQZIx50P++3mKM+bMxppbrTKHAGPOyMeZbY8yWk6bVM8a8Z4zZ6v95qcuMgabCEkTGmEjgOeAmoA0wyBjTxm2qkHQceMha2wZIBUZqPV6Q+4FC1yHCwHRgubW2FXAVWqfnzBjTCBgFpFhr2wKRwEC3qULGq8CNp0x7DHjfWtsCeN//PGyosATX1cA2a+0/rLXHgPlAb8eZQo61dre19hP/4x/xbRgauU0VmowxjYEewEuus4QyY0wdIA2YDWCtPWatPeA2VciKAi4yxkQBccA3jvOEBGvth8D3p0zuDbzmf/wa0KdaQwWZCktwNQK+Oun5LrShvSDGmCSgPfCx2yQh63+AMUCF6yAhrhmwD3jF//XaS8aY2q5DhRpr7dfAVGAnsBs4aK1d4TZVSEuw1u72P94DJLgME2gqLBIyjDHxwJvAA9baQ67zhBpjzM3At9baDa6zhIEooAMwy1rbHigmzHa/Vwf/MRa98RXAhkBtY8ztblOFB+s7BTisTgNWYQmur4EmJz1v7J8m58gYE42vrORYa3Nd5wlRnYBexpgd+L6evMEYM89tpJC1C9hlrf1pT99CfAVGzs2/AV9aa/dZa8uAXOBax5lC2V5jTAMA/89vHecJKBWW4FoPtDDGNDPGxOA7mGyR40whxxhj8B0rUGit/W/XeUKVtTbTWtvYWpuE729xlbVWn2bPg7V2D/CVMaalf1JXoMBhpFC1E0g1xsT53+dd0cHLF2IRMNj/eDDwjsMsARflOkA4s9YeN8bcC7yL7+j3l621nzmOFYo6AXcAfzfGfOqfNtZau9RhJpH7gBz/h5F/AP/uOE/IsdZ+bIxZCHyC72zAjYTxlVoDyRjzZ+A6oL4xZhfwODAZ+IsxZihQBAxwlzDwdKVbERER8Tx9JSQiIiKep8IiIiIinqfCIiIiIp6nwiIiIiKep8IiIiIinqfCIiLVxn+H43v8jxv6T2kVEflFOq1ZRKqN/15QS/x35hUROWu6cJyIVKfJwJX+CwBuBVpba9saY4bgu7NsbaAFvhvixeC7YGAp0N1a+70x5krgOeDXQAlwl7X2/6r/nyEi1U1fCYlIdXoM2G6tbQc8csprbYG+wL8CWUCJ/8aC64A7/fNkA/dZa38HPAw8Xy2pRcQ57WEREa/4wFr7I/CjMeYgsNg//e9Asv9u3dcCb/huOwNAbPXHFBEXVFhExCtKT3pccdLzCnxjVQRwwL93RkRqGH0lJCLV6Ufg4vP5H621h4AvjTH9wXcXb2PMVYEMJyLepcIiItXGWrsf+MgYswV45jwWkQEMNcZsAj4Degcyn4h4l05rFhEREc/THhYRERHxPBUWERER8TwVFhEREfE8FRYRERHxPBUWERER8TwVFhEREfE8FRYRERHxPBUWERER8bz/Bx2V5L4X+UhbAAAAAElFTkSuQmCC\n",
"text/plain": [
"<Figure size 648x504 with 1 Axes>"
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
}
],
"source": [
"plot_sampling_points(func, N=10, peg=\"midpoints\")"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"This exhibits a serious problem: After the transformation back to the sampling points on the time grid at the end of the optimization, it does not preserve the boundary values of the control to be zero at t=0 and t=T. Hence the choice to sample the control fields on input. Note that the same also applies to the `shape` function in the pulse options (`pulse_options` argument in [`optimize_pulses`](https://krotov.readthedocs.io/en/latest/API/krotov.optimize.html#krotov.optimize.optimize_pulses)): these are also first sampled on the points of the time grid, and then converted to the interval representation, so that initial and final values of zero are preserved.\n",
"\n",
"Obviously, when using just a few grid points as in the above plots, the result of a propagation may depend significantly on the sampling (the blue/red points can deviate significantly from the \"exact\" curve). This is the \"discretization error\": In order for the optimization to be numerically useful, we must choose a time grid with a sufficiently small dt so that the difference between the two sampling choices becomes negligible:"
]
},
{
"cell_type": "code",
"execution_count": 8,
"metadata": {
"attributes": {
"classes": [],
"id": "",
"n": "7"
}
},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAiwAAAGtCAYAAAA1cy8JAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAIABJREFUeJzs3Xl8VNX5x/HPmUACyKYsokASVJRVtrAV0SAuVNT+xKLFBAFZRMS1BTe0wk/6s9KqtbVVEAFh3MCqaK27EQFREHdQBElCACFsSkgCSeb8/pgkJGGSTEhm7izf9+s1L5g7dybPzNyZPDn3nOcx1lpEREREQpnL6QBEREREqqOERUREREKeEhYREREJeUpYREREJOQpYREREZGQp4RFREREQp4SFhEREQl5SlhEREQk5ClhERERkZBXz+kAaqply5Y2MTHR6TBERESkDnz22Wd7rLWtqtsv7BKWxMRE1q1b53QYIiIiUgeMMRn+7KdTQiIiIhLylLCIiIhIyFPCIiIiIiEv7OawiIiI8woKCsjKyiI/P9/pUCRMNGjQgHbt2lG/fv3jur8SFhERqbGsrCyaNGlCYmIixhinw5EQZ61l7969ZGVl0aFDh+N6DJ0SEhGRGsvPz6dFixZKVsQvxhhatGhRqxE5JSwiInJclKxITdT2eFHCIiIiIiFPCYuIiEgZaWlprF692ukwpAIlLCIiImUoYQlNSlhERCQsLVmyhH79+tGzZ0+uv/56MjIy6NixI3v27MHj8TB48GDefvttAP7nf/6HPn360LVrV+bOnVv6GG+++Sa9e/emR48eDB06lPT0dJ544gkeeeQRevbsyUcffeTU05MKtKxZRERq5dZbb+WLL76o08fs2bMnjz76aKW3b9y4kRdeeIFVq1ZRv359pkyZwocffsgdd9zBDTfcQL9+/ejSpQsXXXQRAE8//TQnnXQSeXl59O3blyuvvBKPx8PEiRNZsWIFHTp0YN++fZx00klMnjyZxo0b84c//KFOn5PUjhIWEREJO++99x6fffYZffv2BSAvL4/WrVtz//33s3TpUp544olySdRjjz3Gyy+/DMC2bdv44YcfyM7O5txzzy2tC3LSSScF/4mI3wKWsBhjngYuBXZba7v5uN0AfwMuAXKBsdba9YGKR0TCg9sN99wDmZkQHw+zZ0NKitNRSVWqGgkJFGstY8aM4f/+7//Kbc/NzSUrKwuAnJwcmjRpQlpaGu+++y4ff/wxjRo1Ijk5WRV6w1Ag57AsBIZVcfuvgY7Fl0nAvwIYi4iEKrcbEhPB5SKnZSLvjnOTkQHWQkYGTJoEK6cc3YfERO99JKoNHTqUZcuWsXv3bgD27dtHRkYGd9xxBykpKcyaNYuJEycC8PPPP3PiiSfSqFEjvvvuO9asWQPAgAEDWLFiBVu3bi19DIAmTZpw8OBBB56VVCVgCYu1dgWwr4pdfgM8Y73WAM2NMacEKh4RcZ67Qt6xcorbm5EUZyiN92bwj4JJjOJoQvKbXDe9/nV0HzIyKLxuEje3dCt/iWJdunThgQce4KKLLuLss8/mwgsvJD09nbVr15YmLbGxsSxYsIBhw4ZRWFhI586dufPOOxkwYAAArVq1Yu7cuYwYMYIePXpw9dVXA3DZZZfx8ssva9JtiDHW2sA9uDGJwOuVnBJ6HXjQWruy+Pp7wB3W2nU+9p2EdxSG+Pj4PhkZGQGLWUQCw12cm+TmHt2WTjwJbDtm33QS6EA6AFtJJJFjP/Nl92nUCObO1amjYNq4cSOdO3d2OgwJM76OG2PMZ9bapOruGxbLmq21c621SdbapFatWjkdjojUQHKy9zJ+fPlkBaA9WT7vE0+mz/9Xtk9urvfxk5NrGayIhCwnE5btQPsy19sVbxORSFB8/uf9D108vyaREYePPW+TSbzPu5bdnmWq3wdgxGE3z69J1DwXkQjlZMKyHLjWeA0AfrbW7nQwHhGpK+6jc1NcWNoczmAeE8vNTQG4m9nkmkblthXGNuLhFrMxBhISIHPybO85nzIO0Yi7mV16fRRunmIibQ5XmK2rpEUkYgQsYTHGPAd8DJxljMkyxow3xkw2xkwu3uUN4EdgMzAPmBKoWEQkiCo5/3MCefwfd5Xb9mqjFNZPnuvNTIozlHpPz+WxPSl4PJCeDuf8M8U7QaV4n5wWCUytP5fnODph5f+4i0bklY+j5DyRiESEgNVhsdaOquZ2C9wYqJ8vIg46fNjn5ni2kRD3E5lH2pTWWDknJQWoZrZsSkrpjNrGwAVu+OAeyMzwEB+3m/jDx07crSoOEQk/YTHpVkRCX+mS5RVpbHP5nndiEhJIz29TOnpyvKt6UlK89/dYF+n5bTAJCT73y4pJ0JQWkQihhEVEas1dvpwKd3j+xCEalt+pUSPvkEogzPY1z6Uh04tma0qL1Ink5GTWrTum6kal0tLSuPTSS+s8jgkTJrBhw4Yq93nllVeq3acuLFy4kKlTpwb855RQwiIix62yJcvPkcJE5pHB0bkpAS2UknJ0nosHQwbxTGReuXkuWvrsrIpFA5U8Hp+nnnqKLl26VLnP8SQshYWFtQkrKJSwiEit+Zoq8hwpJJJOrc//+Kv4PNH553lIJKNcslJVnBJ4FUfg6mLE69ChQwwfPpwePXrQrVs3XnjhBQBmzZpF37596datG5MmTaKkOGpycjK33XYbSUlJdO7cmbVr1zJixAg6duzIjBkzAEhPT6dTp06kpKTQuXNnfvvb35JbsXgQ8PbbbzNw4EB69+7NyJEjycnJAeDNN9+kU6dO9O7dm3//+98+487Pz2fcuHF0796dXr168cEHHwDe0YoRI0YwbNgwOnbsyPTp033ev+xIT+PGjbnnnnvo0aMHAwYMYNeuXaxevZrly5czbdo0evbsyZYtW9iyZQvDhg2jT58+DB48mO+++w6AsWPHMnnyZPr378/06dNJTEzkwIEDpT+rY8eO7Nq1i9dee43+/fvTq1cvLrjgAnbt2nVMXEuXLqVbt2706NGDc889t/o38HhYa8Pq0qdPHysioeXkk/Os91dR+UtCgjPxJCQcGwtYGx/vcSagCLRhw4Zq9znvPO8lLs73+xEX5739eCxbtsxOmDCh9PqBAwestdbu3bu3dFtqaqpdvnx5cSzn2enTp1trrX300UftKaecYnfs2GHz8/Nt27Zt7Z49e+zWrVstYFeuXGmttXbcuHF2zpw5pfdfu3atzc7OtoMHD7Y5OTnWWmsffPBBO3PmTJuXl2fbtWtnN23aZD0ejx05cqQdPnz4MXH/5S9/sePGjbPWWrtx40bbvn17m5eXZxcsWGA7dOhgDxw4YPPy8mx8fLzNzMz08Zp647DWWqD0+U2bNs3+7//+r7XW2jFjxtilS5eW3uf888+3mzZtstZau2bNGjtkyJDS/YYPH24LCwuttdbefPPN9umnny7db+jQodZaa/ft22c9Hu9nZ968efb222+31lq7YMECe+ONN1prre3WrZvNysqy1lq7f/9+32+a9X3cAOusH7//NcIiIjVXZnz/yKmncvGeizGm/LLiQE5ZqY6PKS3AIaa1nabzEg6obGSrNiNe3bt355133uGOO+7go48+olmzZgB88MEH9O/fn+7du/P+++/z7bfflt7n8ssvL71v165dOeWUU4iLi+O0005j2zbvSrP27dszaNAgAFJTU1m5cmW5n7tmzRo2bNjAoEGD6NmzJ4sWLSIjI4PvvvuODh060LFjR4wxpKam+ox75cqVpbd16tSJhIQENm3aBHgbOjZr1owGDRrQpUsXqmtDExsbWzpPpk+fPqSnpx+zT05ODqtXr2bkyJH07NmT66+/np07j5Y8GzlyJDExMQBcffXVpSNVzz//fGlvpaysLC6++GK6d+/OnDlzyr2mJQYNGsTYsWOZN28eRUVFVcZ9vJSwiEjNVBjfj925k38VreCN1GfKllNxtLdPSvnSLcTHW2acNplxH/+1bs9LSJXS0ryXShZxkZDgvf14nHnmmaxfv57u3bszY8YMZs2aRX5+PlOmTGHZsmV8/fXXTJw4kfz8/NL7xMXFAeByuUr/X3K9ZA6HMabcz6l43VrLhRdeyBdffMEXX3zBhg0bmD9//vE9iQrKxhQTE1PtvJL69euXxlfZ/h6Ph+bNm5fG+8UXX7Bx48bS20844YTS/w8cOJDNmzeTnZ3NK6+8wogRIwC46aabmDp1Kl9//TVPPvlkude0xBNPPMEDDzzAtm3b6NOnD3v37q3Zk/eDEhYR8V8lReEaAcNevIX0xOSgTVmpTunS53OTyegwhJnbl3JCxZ00EzcofI141XYEbseOHTRq1IjU1FSmTZvG+vXrS3+RtmzZkpycHJYtW1bjx83MzOTjjz8G4Nlnn+Wcc84pd/uAAQNYtWoVmzdvBrxzaTZt2kSnTp1IT09ny5YtADz33HM+H3/w4MG4i5PkTZs2kZmZyVlnnVXjOKvSpEkTDh48CEDTpk3p0KEDS5cuBbwJ15dffunzfsYYrrjiCm6//XY6d+5MixYtAPj5559p27YtAIsWLfJ53y1bttC/f39mzZpFq1atSkes6pISFhGpmUCM7weBK0zjjgQVR7zqYgTu66+/pl+/fvTs2ZOZM2cyY8YMmjdvzsSJE+nWrRsXX3wxffv2rfHjnnXWWTz++ON07tyZ/fv3c8MNN5S7vVWrVixcuJBRo0Zx9tlnM3DgQL777jsaNGjA3LlzGT58OL1796Z169Y+H3/KlCl4PB66d+/O1VdfzcKFC8uNrNSF3/3ud8yZM4devXqxZcsW3G438+fPp0ePHnTt2pVXX3210vteffXVLFmypPR0EMD999/PyJEj6dOnDy1btvR5v2nTptG9e3e6devGr371K3r06FGnzwnA2OIZ1OEiKSnJ1mQtvIjUscRE7+mUihISvEMaoSpc4w5RGzdupHPnzk6HUafS09O59NJL+eabb5wOJWL5Om6MMZ9Za5Oqu69GWESkRjaNHcuhihudnGHrLx/nJY7Uqxf6cYsIoIRFRPxwdFGQpcsDE7i98cUUtWsXGjNs/VXhvMTexo0ZW1jInV+frYVDAkBiYqJGV0JYwJofikhkKFkU5J1naygqaseigv9w7oMxIZ+jHKNME8VGeXl8eMZ97HzodErOjJcsHCrZVURChxIWEfGpZOHMmjXHzks9fDiG8ePD95e697k1JDv7AawtP+GxZOFQuD43kUilU0IiUqVIXlxTUOB7dUYkPDeRSKMRFhHxqaSgV0KCJTPTHHN7ZcXAwkHJc6tq4ZCIhBaNsIhIlYYOfQ8qrAsKh0VB/vBd0MxGxHOTulW26aA/0tLSSsvm16VHH320XEPGSy65pFzDwmC77777ePfdd4/ZHojnr4RFRCq1bds2Xnzxf+ja9W/Ex9uwWhTkj7ILh8AC6Ywc+U5EPLeQU6b/lJZjHb+KCcsbb7xB8+bNHYmlqKiIWbNmccEFFwTl5ylhEZFK3XbbbXg8Hl57bRQZGSZkyu7XpdIS/h44//zxvPLKVezatcvpsCJLhf5TddHH6dChQwwfPpwePXrQrVu30qZ9s2bNom/fvnTr1o1JkyZRUhw1OTmZ2267jaSkJDp37szatWsZMWIEHTt2ZMaMGYC3cFynTp1ISUmhc+fO/Pa3vy2XHJR4++23GThwIL1792bkyJHk5OQA8Oabb9KpUyd69+7Nv//9b59x5+fnM27cOLp3706vXr344IMPAFi4cCEjRoxg2LBhdOzYkenTpx9z38cee4wdO3YwZMgQhgwZAniXYu/Zs6c09rFjx3LmmWeSkpLCu+++y6BBg+jYsSOffvpp6et23XXX0a9fP3r16uWz6q3H42HKlCl06tSJCy+8kEsuuaS0zUFiYiJ33HEHvXv3ZunSpYwdO7b0Nn+ef63409I5lC59+vSptG21iNSBJUusTUiwHmPsVrAvjxzpdERBs3HjRpvqctnsE06w1hhrExK8r4ccY8OGDdXvdN553ktcnLXeVKX8JS7Oe/txWLZsmZ0wYULp9QMHDlhrrd27d2/pttTUVLt8+fLiUM6z06dPt9Za++ijj9pTTjnF7tixw+bn59u2bdvaPXv22K1bt1rArly50lpr7bhx4+ycOXNK77927VqbnZ1tBw8ebHNycqy11j744IN25syZNi8vz7Zr185u2rTJejweO3LkSDt8+PBj4v7LX/5ix40bZ631Hm/t27e3eXl5dsGCBbZDhw72wIEDNi8vz8bHx9vMzMxj7p+QkGCzs7OPub5161YbExNjv/rqK1tUVGR79+5tx40bZz0ej33llVfsb37zG2uttXfddZddvHixtdba/fv3244dO5Y+lxJLly61v/71r21RUZHduXOnbd68uV26dGnpz/vzn/9cuu+YMWPs0qVL/X7+vo4bYJ314/e/RlhE5Kgyfwkba0kEfvOf/0TN8H2nzz5jvstFy0OH1NG5LgVgqVn37t155513uOOOO/joo49o1qwZAB988AH9+/ene/fuvP/++3z77bel97n88stL79u1a1dOOeUU4uLiOO2000qb9bVv355BgwYBkJqaysqVK8v93DVr1rBhwwYGDRpEz549WbRoERkZGXz33Xd06NCBjh07YowhNTXVZ9wrV64sva1Tp04kJCSwadMmAIYOHUqzZs1o0KABXbp0IcPXjPAqdOjQge7du+NyuejatStDhw7FGEP37t1JL24/8fbbb/Pggw/Ss2dPkpOTyc/PJzMz85gYR44cicvlok2bNqWjOSXK9hkq4e/zrw2tEhIRr+Rkn0VXTElhknnzji6viUTFzz+2sLD89mh5/oHgz3Ks43xNzzzzTNavX88bb7zBjBkzGDp0KNOnT2fKlCmsW7eO9u3bc//995d2cAZKmwy6XK5yDQddLheFxe+7MeVXxFW8bq3lwgsvPKYb8xdffHFcz6OssjHFxMSUxnQ89y/7HMs+P2stL730Uq06RJ9wwjF9z4NCIywiclQkF13xR7Q//0DxvRyrVkvNduzYQaNGjUhNTWXatGmsX7++NDlp2bIlOTk5pXMraiIzM5OPP/4YgGeffZZzzjmn3O0DBgxg1apVbN68GfDOCdm0aROdOnUiPT2dLVu2AByT0JQYPHgw7uIRu02bNpGZmVmj5KFJkyYcPHiwxs+rxMUXX8zf//730rk9n3/++TH7DBo0iJdeegmPx8OuXbtI8yOp9Pf514YSFhHxSkvDxsf7vq0WfwmHjbS0yguwRMPzD6QKfZzqYqnZ119/Tb9+/ejZsyczZ85kxowZNG/enIkTJ9KtWzcuvvhi+vbtW+PHPeuss3j88cfp3Lkz+/fv54Ybbih3e6tWrVi4cCGjRo3i7LPPZuDAgXz33Xc0aNCAuXPnMnz4cHr37k3r1q19Pv6UKVPweDx0796dq6++moULF5YbGanOpEmTGDZs2DGnafx17733UlBQwNlnn03Xrl259957j9nnyiuvpF27dnTp0oXU1FR69+5desqtMv4+/9owJVlWuEhKSrI1WQsvIv779NZb6fq3v1FuwLdRo8hZx1yd8o2TAMiPiaHBokXR8fxrYOPGjXTu3NnpMOpUeno6l156qRogAjk5OTRu3Ji9e/fSr18/Vq1aRZs2bWr9uL6OG2PMZ9bapOruqxEWEQGgoKCAlP/8h5lt23pHWiKt6Io/KowEHGjWjOuKivj4tNOcjkwkqC699FJ69uzJ4MGDuffee+skWaktTboViXJuN9xzD2Rk1APewfP7/Zi/9HI6LOeU6ehcLyeH9884g8/Gvk1+/gC2bTPEx3unXkRLDhdNEhMTNbpSzJ95K8GmERaRKFa2nhcYIJF//aunVvEWa9y4MZdcsoRNm/5AZqbRSucKwm1KgTirtseL5rCIRKHkZO+/PlYxAxAXB2VWg0alo6+R5fDhY5s/RvtrtHXrVpo0aUKLFi2OWforUpG1lr1793Lw4EE6dOhQ7jZ/57DolJBIFNMq3ur5Sla824McSIhp164dWVlZZGdnOx2KhIkGDRrQrl27476/EhaRKORPPa9op9eoavXr1z/mL2WRQNIcFpEoduON24FD5bbVsp5XxAlAzTMROQ5KWESi2OrVU4mLu4l27YqichWzP8qudAYPxmTy17/+otdIJMiUsIhEqXXr1vHKK69wzz0d2LYtBo8H0tOVrPiSkuJ9bTZu3IQxHdiy5X+dDkkk6ihhEYlSs2bN4sQTT+SWW25xOpSw0alTJ0aNGsU///lPdu/e7XQ4IlFFCYtIFFq/fj2vvfYat912G02bNnU6nLAyY8YM8vLy+Otf/+p0KCJRRQmLSDRxuyExkZ59+pBhDL8/5RSnIwo7nTp14rEBA7hxzhysy+VdRqQqciIBp2XNItGiTGM/FxBvLdxyCzRsqIkrNeF2M+Xzz3GVFN0sKX0Leh1FAkgjLCLRIDkZxo8v14UY8F4fP/5oWVepWvHr6KpY4lavo0jAKWERiRYqa1s39DqKOEIJi0g0SEurvDRrQsLRsq5SNb2OIo5RwiISJTInT65Q0xaVbD0ePkrfHqlfX6+jSIApYRGJEnd+9RU3xcVR1K4dKmtbC2VL3xpDdqNG3BATwy+XXeZ0ZCIRTQmLSBRIT0/nxRdf5MQbbyRm2zZU1raWSkrfejxkfPghT+fnM2/ePKejEoloSlhEosAjjzyCMYZbb73V6VAiTlJSEkOGDOGRRx7hyJEjTocjErGUsIhEMLcb2rcv4rHHHiEubicrVrR3OqSING3aNLZvP5dTTz2CasmJBIYKx4lEqKN14mIAOHSopeqbBci+fcMwJpm9exsCqiUnEgjGllRrDBNJSUl23bp1TochErJKapetWeO7NEhcHFSseybHR6+1SO0ZYz6z1iZVt59OCYlEKNU3Cx691iKBp1NCIhGmpHZZQoIlM9Mcc3tldc+k5kpe68RE72mgivRai9QdjbCIRKgrrlgLFUrFqU5cYPioJafXWqSOKWERiUDWWj755BZatbqH+HirOnEBVraWHFgggwce2KXXWqQOKWERiUCffPIJa9as4b77ziAjw6hOXBCU1JLbvn0n9et3JCPj/5wOSSSiKGERiUCPPfYYTZs2ZezYsU6HEnVOPfVUrrrqKp5++ml++eUXp8MRiRhKWEQizPbt21m6dCnjx4+ncePGTocTlW6++WYOHjzIokWLnA5FJGIoYRGJME888QRFRUVMnTrV6VCiVr9+/RgwYAB///vf8Xg8TocjEhGUsIhECrcbT3w8Mx94gJ8aNOC0jz92OqKo9tfevXn7hx8w9eqpVr9IHVAdFpFIUFyH35WbC0DrvDzVhneS283AhQsxANaqVr9IHVBpfpFw5k9t+AEDjlY4k8BLTtb7IVIDKs0vEk1UGz606P0QqXMBTViMMcOMMd8bYzYbY+70cXu8MeYDY8znxpivjDGXBDIekYiTlua9VFYDPiFBf80Hm94PkYAIWMJijIkBHgd+DXQBRhljulTYbQbworW2F/A74J+Bikckku29/fYKRfhRbXgn+ajVX1C/vt4PkVoI5AhLP2CztfZHa+0R4HngNxX2sUDT4v83A3YEMB6RiPXwrl1MAgrbtkV1+ENA2Vr9xvBTXBzTmjfHM2qU05GJhK1AJixtgW1lrmcVbyvrfiDVGJMFvAHc5OuBjDGTjDHrjDHrsrOzAxGrSNg6cuQITz31FDmXX069rCxUhz9ElNTq93j4cNEi/padzdtvv+10VCJhy+lJt6OAhdbadsAlwGJjzDExWWvnWmuTrLVJrVq1CnqQIqHs5ZdfZvfu3dxwww1OhyKVuOKKK2jdujVPPPGE06GIhK1AJizbgfZlrrcr3lbWeOBFAGvtx0ADoGUAYxKJOP/617/o0KEDF110kdOhSCViY2MZP348r732Gtu2bav+DiJyjEAmLGuBjsaYDsaYWLyTapdX2CcTGApgjOmMN2HROR+Rarjd3uKpLpflww8X0q/fo7hcTg+YSlUmTZqEx/M7undvgsul4rciNRWwbzhrbSEwFXgL2Ih3NdC3xphZxpjLi3f7PTDRGPMl8Bww1oZbJTuRICsuaktGBlhrgESWL79Mv/xC3KpVibhc8/n55+blit/qfRPxjyrdioQRFVENPypGLFI1VboViVAqohqe9L6J1I4SFpEwoiKq4UfFiEXqhhIWkTAzezY0aFBUbpuK2oY+H8Vv9b6J1IASFpEwk5IC/fs/jTGZGGNV1DZMlC1+Cx5crm384x9H9L6J+EkJi0iY+fnnn1m79lYmTHgAj8eoqG0YKSl++957aXg88dSr94LTIYmEDSUsImHmueeeIzc3l0mTJjkdihynIUOGcMYZZ/DUU085HYpI2FDCIhJm5s2bR8+ePenTp4/TochxMsYwYcIEVqxYwffff+90OCJhQQmLSBhZv34969evZ8KECRhjnA5HamHMmDHExMQwf/58p0MRCQtKWETCyPz582nQoAHXXHON06FILbVp04bLLruMRYsWceTIEafDEQl5SlhEwoHbjSc+nr//859si4nhxDfecDoiqQMTJkxg6O7dFLRtixoMiVStntMBiEg1ipsHuXJzAWh56JC3CQ1oeVCYG7ZvH0OModGePd4NJQ2GQO+tSAXqJSQSqvxpQpOfH9SQpI7ovRUppV5CIpFCTWgil95bEb/plJBIqCppMJOY6D1VUFFlzWkk9Om9FakxjbCIhLjCmTPJrbhRTWgigxoMifhNCYtIiHu1cWMmAHmtW4MxqHlQBCluMGTj4/EAuxs21HsrUgklLCIhbv78+XzUrh2xO3aAx4OaB0WYlBRMRgZ/uO022hYUkH3RRU5HJBKSlLCIhLDt27fz1ltvlVZFlch13XXXUVhYiFt1WER8UsIiEsIWL16Mx+Nh7NixTociAdatWzeSkpJYsGAB4VZuQiQYlLCIhChrLQsWLGDw4MGcccYZTocjQTBu3Di++uorPv/8c6dDEQk5SlhEQtTHH3/Mpk2bNLoSRUaNGkVcXBwLFy50OhSRkKOERSTEuN3e8hyDBg3EmAysHeV0SBIkJ554Ir16zeHxx6fhclm1FhIpQ4XjREJIcdsgvG1XZz23AAAgAElEQVSDDNbGc/PN0KCBFgZFA7cbPv98Ch6Pd4K1WguJHKURFpEQkZwM48eXJCtH5eZ6t5e0n5HIVPL+Hz5cfjWY3n8RLyUsIiFErWWim95/kcopYREJEWlplbeQSUg42n5GIpPef5GqKWERCSGzZ0P9+kfKbVNrmejhu7WQ1fsvghIWkZAyapSHZs2m0aDBLrUNikLFrYWKR1oskM60aT/o/RdBCYtISFmxYgV79jzGU0+9q7ZBUSolxfu+//JLDo0adeWnnx52OiSRkKCERSSELFq0iCZNmnDFFVc4HYo4rEmTJlx55ZW88MIL5OfnOx2OiOOUsIiEiEOHDrFs2TJGjhxJo4oTGSQqXXvttRw4cIDXXnvN6VBEHKeERSREvPzyy+Tk5DBmzBinQ5EQMWTIENq1a8eiRYucDkXEcUpYRELEM888Q2JiIuecc47ToUiIiImJITU1lTfffJNdu3Y5HY6Io5SwiDjN7aawXTvefOcdPt+/H9dzzzkdkYSQa6+9lquKimjYuTO4XKjBkEQr9RIScVJx86B6xfX4m//8s5rHSDmd169nvstFw/37vRvUYEiilLHWOh1DjSQlJdl169Y5HYZI7ZQ0hlmzxnfd9bg40MqQ6KZjRKKEMeYza21SdfvplJCIk9Q8RqqjY0QE0CkhEWeUNIZJTPQO8VdUWVMZiR46RkTK0QiLiIMKZ84kt+JGNQ+Ssnw3GNIxIlFHCYuIg15v1owJQG6rVqh5kPhU3GDIxsfjAbIbNdIxIlFJCYuIg5555hneP/lkYnfsQM2DpFIpKZiMDG6ZOpX2RUUcGD7c6YhEgk4Ji4hD9u3bx+uvv84111xDvXqaTibVGz16NIcPH2bZsmVOhyISdEpYRBzy4osvUlBQwOjRo50ORcJE3759OfPMM1m8eLHToYgEnRIWEYcsXryYrl270rNnT6dDkTBhjGH06NGsWLGC9PR0p8MRCSolLCIO2LJlC6tXr2b06NEYY5wOR8JIamoqAG6V55coo4RFxAFLlizBGEOKJthKDSUmJjJ48GAWL15MuFUqF6kNJSwiQeR2Q0KC5f777yU2dgcfftjO6ZAkDJ155v18//2bxMSoF6JEDy1NEAmS4j6H5OYawHD4cBv1sJMac7vh2WeHAAZr1QtRoodGWESCIDkZxo+H3AplbXNzvdtL+tyJVKXkOMrLKz/vSceRRAMlLCJBoh52Uhd0HEm0UsIiEgRpaZX3qktIONrnTqQqOo4kmilhEQmS2bMhJia/3Db1sJOa8tULsWFDq+NIIp4SFpEgueSS/cAkmjTZpz6HctyKeyGSkADGWCCdlJQPdRxJxFPCIhIky5Yto6hoMe+//6P6HEqtpKR4j5+iIujY8SK2bJnldEgiAaeERSRIlixZQqdOnejTp4/ToUiEMMaQmppKWloa27ZtczockYBSwiISBBkZGaxYsYLU1FSV4pc6lZKSgrWW5557zulQRAJKCYtIEDz77LMAXHPNNQ5HIpHm9NNPZ+DAgSxZssTpUEQCSgmLSIBZa1m8eDHnnHMOHTp0cDociUCpqal8/fXXfPXVV06HIhIwSlhEAsntpqBtW77ZuJE3NmxQ0xcJiKuuuopUY2g/eDC4XGowJBGp2oTFGHOmMeY9Y8w3xdfPNsbMCHxoImGuuHlQ7M6duIAm+/Z5m77oF4nUsZZvvcU8l4sTf/mFcg2GdKxJBDHVtSc3xnwITAOetNb2Kt72jbW2W7UPbsww4G9ADPCUtfZBH/tcBdwPWOBLa22VJ/mTkpLsunXrqvvRIs4paeiyZo3veulxcZCff+x2kZrSsSYRwBjzmbU2qbr9/OnW3Mha+2mFlQ2FfgQQAzwOXAhkAWuNMcuttRvK7NMRuAsYZK3db4xp7Uc8IuFBTV8kWHSsSRTwJ2HZY4w5He8ICMaY3wI7/bhfP2CztfbH4vs9D/wG2FBmn4nA49ba/QDW2t01iF0kNJU0dElM9A7NV1RZMxiRmtKxJlHEn0m3NwJPAp2MMduBW4Eb/LhfW6BsJaOs4m1lnQmcaYxZZYxZU3wK6RjGmEnGmHXGmHXZ2dl+/GgR5+Xfdx+5FTeqeZAEgq8GQzrWJMJUm7BYa3+01l4AtAI6WWvPsdam19HPrwd0BJKBUcA8Y0xzHzHMtdYmWWuTWrVqVUc/WiSwXm7YkAlA/skno+ZBElDFDYZsfDweYHfDhjrWJOJUekrIGHN7JdsBsNY+XM1jbwfal7nernhbWVnAJ9baAmCrMWYT3gRmbTWPLRLylixZwtft2xObnu5daioSSCkpmJQU7pg+nUceeYSdF19MS6djEqlDVX2LNim+JOE9BdS2+DIZ6O3HY68FOhpjOhhjYoHfAcsr7PMK3tEVjDEt8Z4i+rEG8YuEpN27d/PWW2+RkpKCS8mKBFFqaiqFhYUsXbrU6VBE6lSl36TW2pnW2pl4R0Z6W2t/b639PdAHiK/uga21hcBU4C1gI/CitfZbY8wsY8zlxbu9Bew1xmwAPgCmWWv31u4piTjvhRdeoKioiBQNyUuQnX322XTr1k2l+iXi+FOH5XvgbGvt4eLrccBX1tqzghDfMVSHRcJB//79OXz4MF988YXToUgU+vOf/8ydd97Jli1bOO2005wOR6RK/tZh8Wes+hngU2PM/caY+4FPgEW1jE8kYv3www98+umnpKamOh2KRKmSJpslTTdFIoE/q4RmA+OA/cWXcdbaPwU6MJFw5Xa7McYwatQop0ORKNW+fXvOO+88lixZQnWj6CLhwp9eQvHAHuDl4sve4m0iUobbDQkJlpkz7yM2dgdpaRXLDokEzxln3Mf3379JTIx6IUpk8KfS7X8ornILNAQ6AN8DXQMVlEi4Ke5zSG6uAQyHD7dh0iTvbZp3K8HmdsOzzw4BTLleiKDjUcJXtZNuj7mDMb2BKdbaCYEJqWqadCuhxJ/ecwMGHK2gLhJoyck6HiW81OWk23KsteuB/scVlUiEUu85CSU6HiUSVXtKqELFWxfeonE7AhaRSBjxp/ec/pqVYEpL0/EokcmfEZYmZS5xeOe0/CaQQYmEm9mzISam/J+v6j0nTvHVC7FhQ6vjUcKaP5NuN1hry9V4NsaMBFT3WaTY5ZcfxOW6iYYN/8qhQy2Ij/f+0tAER3FCyXF3zz2QmWmxNoOxY7NISTnH2cBEasGfEZa7/NwmErVefvllCgoW8d//bsTjgfR0JSvirJQU73FYWGiJjz+P9HSVz5LwVlW35l8DlwBtjTGPlbmpKVAY6MBEwsmSJUtITExk0KBBTociUo7L5SIlJYWHHnqIXbt2cfLJJzsdkshxqWqEZQewDsgHPitzWQ5cHPjQRMLDjh07eO+990hNTcUY43Q4IsdITU2lqKiIF154welQRI5bpSMs1tovgS+NMe7izssi4sPzzz+Px+NRZ2YJWV26dKFXr14sWbKEm2++2elwRI5LpSMsxpgXi//7uTHmq4qXIMUnEvKWLFlCUlISnTp1cjoUkUqlpqaydu1avv/+e6dDETkuVZ0SuqX430uBy3xcRKLet99+y+eff67OzBLyfve73+FyuXCrqZCEqUoTFmvtzuJ/M3xdgheiSAhyuyExkS7dupEOjK1f3+mIRKp06qmn8kCXLkz605+wLpc6IkrYqWqV0EGONj0EMMXXDWCttU0DHJtIaDra6RADJABMmwbNmmkts4Qut5tpmzZRr6jIe10dESXM1Lj5odPU/FAcpc5yEm786dCZnx/UkETK8rf5oT+Vbks6NJ+Dd4RlpbX281rGJxK+1FlOwpGOWwlz1Va6NcbcBywCWgAtgYXGmBmBDkwkJKWleTvI+aLOchKK0tKqP25FwoA/pflTgL7W2j9aa/8IDABGBzYskRA2ezaFcXHlt6nToYQ6Xx0RddxKGPEnYdkBNChzPQ7YHphwRMJASgp/796dbS4X1hjvX6hz52riooS2lBTvcZqQgAUygEOPPqrjVsKGPwnLz8C3xpiFxpgFwDfAAWPMYxV6DIlEhX379nHHl1/y8M03Y9TpUMJJcUfEtZ98QiLwnFpJSBjxZ9Lty8WXEmmBCUUkPLz44osUFBQwerTOjEp46tu3L2eeeSaLFy9mwoQJTocj4pdqExZr7aJgBCISLhYvXlzam0UkHBljGD16NPfeey/p6ekkJiY6HZJItfxZJXSpMeZzY8w+Y8wvxpiDxphfghGcSKjZsmULq1evZvTo0erMLGGtpJ2ESvVLuPBnDsujwBighbW2qbW2iarcSrRasmQJxhh1Zpawl5iYyODBg1m8eDHhVkBUopM/Ccs24BurI1qinLWWxYsXk5ycTPv27Z0OR6TWRo8ezffff4+qh0s48CdhmQ68YYy5yxhze8kl0IGJhBK3G0455TBbtmziq6+Wq2ecRIhrgHT69UtSL0QJef6sEpoN5OCtxRIb2HBEQs/RXofeckR79zZWzzgJe2433HrrCcAJgHohSuirtvmhMeYba223IMVTLTU/lGBRzziJVOrhKaHE3+aH/pwSesMYc1EdxCQSltQzTiKRjmsJN/6MsBzEO2Z4GCgADGCdWimkERYJtsRE73B5RQkJ3iK3IuFIx7WEijobYSlexuyy1jbUsmaJRnffnQMcKrdNPeMk3PnqhRgbW6DjWkKWP6eEMMacaIzpZ4w5t+QS6MBEQoW1bmAibdocRr0OJVKU6YWIMZZ69bZz2mkP6riWkOVPpdsJwArgLWBm8b/3BzYskdDxzDPP0KXLl+zYEYt6HUokKe6FiMdj+OMfF/Ddd/eRrvNBEqL8GWG5BegLZFhrhwC9gAMBjUokRPzwww+sXr2aMWPGqBS/RLSSZp6LFy92OBIR3/xJWPKttfkAxpg4a+13wFmBDUskNCxevBiXy6VS/BLxEhISSE5O5plnnlGpfglJ/iQsWcaY5sArwDvGmFcBH3PLRSKLx+PhmWee4YILLqBt27ZOhyMScNdeey2bN29mzZo1Tocicgx/VgldYa09YK29H7gXmA/8T6ADE3HaRx99REZGBtdee63ToYgExW9/+1saNmzIokWLnA5F5Bh+rRIqYa390Fq73Fp7JFABiTjO7YbERM5NTibDGH57RIe7RIcmTZowp1cv7p47F+tyoQZDEkr86SUkEj2ONg7CAPHWwtSpEBurpUES+dxuJn/2GTElc1jUYEhCSLWVbkONKt1KwKjBikQrfxpn6fiXAKnLXkIYYxKMMRcU/7+hMaZJbQMUCUlqsCLRTMe/hDB/CsdNBJYBTxZvaod3xZBIZElL85b99CUhQX9dSuRKS9PxLyHPnxGWG4FBwC8A1tofgNaBDErEMbNnU1C/fvltahwk0cJXgyEd/xIi/ElYDpddFWSMqQeE18QXET/Za67h7lat2BkbixoHSdQp02DIAulA5owZOv4lJPiTsHxojLkbaGiMuRBYCrwW2LBEnLF27Vr+smMHr//jH6hxkESl4gZDe7OzObN+fR7NznY6IhHAv4TlTiAb+Bq4HngDmBHIoEScsnDhQho0aMBVV13ldCgijmrZsiWXXXYZS5YsoaCgwOlwRPyqdOux1s6z1o4EJgGf2HBbCy3ih/z8fJ577jlGjBhBs2bNnA5HxHHjxo0jOzubN954w+lQRPxaJZRmjGlqjDkJ+AyYZ4x5JPChiQTXq6++yoEDBxg3bpzToYiEhGHDhnHyySezcOFCp0MR8euUUDNr7S/ACOAZa21/YGhgwxIJvoULF9K+fXuGDBnidCgiIaFevXqMHj2a119/nd27dzsdjkQ5fxKWesaYU4CrgNcDHI9I0Lnd0K5dIW+++R9++eVLnn8+xumQRELG2LFjKSwcSadODVB7IXGSP72EZgFvASuttWuNMacBPwQ2LJHgONo6yPtR+PnnE9U6RaSML77oiss1n/37GwJqLyTOUS8hiUr+tE7Jzw9qSCIhRZ8RCRZ/ewlVOsJijPk7VRSIs9befJyxiYQMtU4RqZo+IxIqqjolpGEMiVglbVESE71D3BVV1lJFJFroMyKhptKExVq7KJiBiDjhj388zHXXFQFH+6eodYrIUbNnl8zzOrpNnxFxQrWTbo0xH+Dj1JC19vyARCQSRPXrLwXeoHXr+WRnNyQ+3vtFrMmEIl4ln4W777ZkZloaNMhm7tyT9RmRoPNnldAfyvy/AXAlUBiYcESCa/78+Zx2WiY//BCHy59F/iJRKCUFUlIM99wzgwcffJDk5EygrdNhSZTxpzT/Z2Uuq6y1twPJ/jy4MWaYMeZ7Y8xmY8ydVex3pTHGGmOqnSUsUlc2b95MWloa48ePx6VsRaRa1113HR6PR5VvxRH+lOY/qcylpTHmYqDaRivGmBjgceDXQBdglDGmi4/9mgC3AJ/UOHqRWnj66adxuVyMGTPG6VBEwsLpp5/OkCFDmD9/Ph6Px+lwJMr482flZ3hXDH0GfAz8Hhjvx/36AZuttT9aa48AzwO/8bHf/wJ/BrSiX4KmsLCQhQsX8utf/5q2bTW0LeKv8ePHs3XrVtJKlhGJBIk/p4Q6WGtPK/63o7X2ImvtSj8euy2wrcz1LCqc9DTG9AbaW2v/U9UDGWMmGWPWGWPWZWdn+/GjRar25ptvsnPnTsaP9yf3FpESJd3M58+f73QoEmX8OSXUwBhzuzHm38aYl4wxtxpjGtT2BxtjXMDDeEdsqmStnWutTbLWJrVq1aq2P1qimdsNiYlcctllZLpcXHbwoNMRiYSVhg0b8kjfvvzp2Wexai4kQeTPKqFngIPA34uvXwMsBkZWc7/tQPsy19sVbyvRBOgGpBljANoAy40xl1trVbRO6t7RxkG4gPYeD9xwA8TEaB2ziL/cbq5duZIYAGvVXEiCptpeQsaYDdbaLtVt83G/esAmYCjeRGUtcI219ttK9k8D/lBdsqJeQnJckpOrbooyYMDR0p4i4ps+RxIA/vYS8mfS7XpjzIAyD9wfP8r2W2sLgal4Oz1vBF601n5rjJlljLncj58rUrfUFEWk9vQ5Eof4M8KyETgLyCzeFA98j7d4nLXWnh3QCCvQCIsct6qaoqSnBzsakfCkz5HUsVp3ay5jWB3EI+K82bM5PGYMcUVFR7epKYpIzfhoLmQbNsTocyQB5s+y5oyqLsEIUqQuHBg+nOtdLvY2bgzGeP8inDtXEwVFaiIlxfu5SUjAGkM68NG11+pzJAGneuQSNdxuN4sKCkhPSwOPxzt8rS9ZkZpLSfF+foqKuPjMM7n7m2+cjkiigBIWiQrWWp588kn69OlDnz59nA5HJCIYY5g0aRKrVq3i2299LgAVqTNKWCQqfPLJJ3z99ddMKqkXISJ1YsyYMcTGxjJ37lynQ5EIp4RFosKTTz5J48aNGTVqlNOhiESUli1bcuWVV/LMM8+Ql5fndDgSwZSwSMQ7cOAAL7zwAtdccw1NmjRxOhyRiDNp0iQOHDjA0qVLnQ5FIpgSFolobjecdpqLvLwcli9/TC1PRALgvPPOo02b25k06SLUXkgCxZ86LCJhyds6yJKb2xSAn36KU8sTkQB49lnD3r0PUlBQH1B7IQmMaivdhhpVupXqJCd7/62q5Ul+flBDEolI+qxJXajLXkIiYUktT0SCQ581CQadEpKIU9Istn37IrKyYo65PSEhuPGIRKqSz1pV7YVE6opGWCRiJSe/DRwqt02tg0Tq3uzZ3s9WWfqsSV1TwiIRyePx8Omnt3L66Q+RkKDWQSKBVKa9EGCBdO6660d91qROKWGRiPT++++zadMm7r+/I+npah0kEmgl7YX27/+ZRo26snWrhlekbilhkYj0+OOP07JlS0aOHOl0KCJRpXnz5qSmpvLss8+yb98+p8ORCKKERSJOZmYmy5cvZ8KECcTFxTkdjkjUmTJlCvn5+SxYsMDpUCSCKGGRiFPShG3y5MkORyISnXr06MGgQYP417/+hcfjcTociRBKWCRyuN3YhARmzZ7Nzrg4ElaudDoikaj1UI8evLtlC6ZePdXqlzqhOiwSGbx1+DG5uRigdV6eaoOLOMXtZuDChRgAa1WrX+qESvNLePOnNviAAUcrXIlIYCUn6/MoNaLS/BJdVBtcJHTo8ygBoIRFwltamvdSWQ3whAT9NScSTPo8SoAoYZGIkHP33RWK8KPa4CJO8VGrvzA2Vp9HqRUlLBIRHt+/n4nAkVNOUR1+EaeVqdVvjWFH/frce/LJ2GuucToyCWNKWCTsFRQU8I9//INd559P7I4dqsMvEgqKa/Ubj4c3/vlPHty2jQ8//NDpqCSMKWGRsPfyyy+TlZXFrbfe6nQoIuJDSkoKLVq04G9/+5vToUgYU8IiYe/RRx/l9NNPZ/jw4U6HIiI+NGzYkOuvv55XX32VH3/80elwJEwpYZGw5HZ7i2e6XJaPP36WX/3qH7hcOpxFQtWUKVOAFHr2bI7LpeK3UnP6hpewU1zUlowMsNYAiSxbdrG+/ERCWFpaW1yupzh48KRyxW/1uRV/qdKthBUV0RQJL/4Uo87PD2pIEmJU6VYilopoioQffW6ltpSwSFhREU2R8OJPMWoRfyhhkbAzezbExhaU26aitiKhzUfxW31upUaUsEjYGTXKw0kn3Un9+jswxqqorUgYKFP8FrBAOnfd9aM+t+I3JSwSdpYvX85PPz3M4sUf4fEYFbUVCRPFxW/5+eeDNG3ag2++udvpkCSMKGGRsDNnzhwSExO58sornQ5FRI5D06ZNmTx5MkuXLmXr1q1OhyNhQgmLhJVVq1axevVqbr/9durVq+d0OCJynG655RZiYmJ4+OGHnQ5FwoQSFgkPxaVtB55zDpkuFxNPOMHpiESkFk499VT+PmAA0x5/HKvSt+IH/Ykqoa+ktG1uLi6gvccDN93krTilySsi4cntZuLatbhKipeWlL4Ffa7FJ1W6ldClEpkikUklq6UMVbqVyKESmSKRR59rqSElLBK6VCJTJDKpZLUcByUsEvJyZ8wgt+JGlcgUCW8+St8Wxsbqcy2VUsIiIe9v2dlMAA63aQPGoNK2IhGgTOlbawzb69Xjntat8Ywa5XRkEqKUsEhIO3ToEA8//DA/X3IJcTt3gseDStuKRIji0rfG4+H9p5/moawsXn/9daejkhClhEVC2ty5c9mzZw8zZsxwOhQRCaBRo0bRoUMHHnjgAcJt9aoEhxIWCVn5+fnMmTOH888/n4EDBzodjogEUL169bjrrrtYu3Yt77zzjtPhSAhSwiIhp7ioLQ0bxrFz52oGDHjM6ZBEJAiuvfZaTjzxRi6//GxcLqvit1KOKt1KSClT1BYwQCKPPmrp0kXTVkQi3bJlceTkPEJBQX1AxW+lPFW6lZCh4pci0cmfotb6/EcuVbqVsKTilyLRS59/qYoSFgkZKn4pEp38KWqtz78oYZGQMmtWIcaUr2urorYi0cFH8Vvi4or0+RdACYuEmIKCRVg7gVatclXUViTKlCl+izEWl2sbp5/+Z33+BVDCIiHk8OHDzJo1i/79f2TXroYqaisShYqL3+LxGP7615fYsOEePvjgA6fDkhCghEVCxlNPPUVmZiYPPPAAxhinwxERh02ePJm2bdty7733qvqtKGGREOB244mP54apU9kZF8fQXbucjkhEQkCDBg2YMWMG8atWkd+mDbhcqJpc9FLhOHFWcaU4l7dSHG0OH1alKBEpNaFhQ641hoa7d3s3qJpc1FLhOHGGP5Wi8vODGpKIhBB9R0SNkCgcZ4wZZoz53hiz2Rhzp4/bbzfGbDDGfGWMec8YU8kqfIlYqhQlIlXRd4QUC9gpIWNMDPA4cCGQBaw1xiy31m4os9vnQJK1NtcYcwPwEHB1oGKSEFJcBaqofXtisrKOvb2yClIiEh1KKsUlJnpPA1Wk74ioE8gRln7AZmvtj9baI8DzwG/K7mCt/cBaW1IlbA3QLoDxSAhyd+3KoYobVSlOREr4qCZn9R0RlQKZsLQFtpW5nlW8rTLjgf/6usEYM8kYs84Ysy47O7sOQxQn/fDDD4x/7z2eGzKkpFKUKsWJSHllqslZY0gHXr/sMn1HRKGQWCVkjEkFkoDzfN1urZ0LzAXvpNsghiYBdNdddxEXF8elzz4Lbdo4HY6IhKqUFEhJwQBTL72Uj/77X7bs2UPLli2djkyCKJAjLNuB9mWutyveVo4x5gLgHuBya61mUUU4t9t7Strlsrz00l+4+OJFtFGyIiJ+euihh/jll0s57TSXyrJEmUAmLGuBjsaYDsaYWOB3wPKyOxhjegFP4k1WdgcwFgkBxSVXyMgAaw2QyH//O0JfNiLit88/70K9ek9z8OBJWHu0LIu+RyJfwBIWa20hMBV4C9gIvGit/dYYM8sYc3nxbnOAxsBSY8wXxpjllTychLnkZBg/HnLLN2ImL88wfvzRkgsiIpUp+R4pLIwrtz03F32PRIGAzmGx1r4BvFFh231l/n9BIH++hBaVUxCR2tL3SPRSLyEJirS0yssmJCQcLbkgIlIZfY9ENyUsEjTTpx8Ayp8TUjkFEakJH2VZiI0t1PdIFFDCIkHz6ae34nJdz6mnFqjkiogclzJlWTDGEhu7k7i4qQwffsDp0CTAlLBIUKxevZpFixYxfXo7tm+vj8cD6elKVkSk5lJSvN8fHo/h4493kpMzl/vvv9/psCTAlLBIYLnd2IQEBgwaxLaYGP7YsaPTEYlIBOnduzeTJ09mz2OPceTUU1FxlsgVEpVuJUIVF14xubkYoF1REdx0k7ctvIZWRKSOPNSzJy5rid2507uhpDgL6Lsmghhrw6vSfVJSkl23bp3TYUhVSoohrFnje61hXBzk5wc1JBGJQPquiQjGmM+stUnV7adTQhI4KpggIsGg75qooFNCUveKiyHkt2lDg127jr29skIKIiI1UVJ4JTHRexqoIn3XRBSNsEhAHD58mHtjYsg1pvwNKvUB99EAAA+nSURBVLwiInXNR3GWogYN9F0TYZSwSED86U9/4i87dvDDH/5QUjBBhVdEJDDKFGexxpAVE8NdLVpwZORIpyOTOqSEReqM2+0dmXW5LLNmXcegQY/T46GHSgomqPCKiAROcXEW4/Hw5auvMmf7ubRunatVzhFECYvUieIVzGRkgLUGSGD9+hv0JSEiQXfgwHBiYhbw88/NsfboKmd9H4U3JSxSayUt33PLtwkiL8+o5buIBFXJ91FRUVy57bm56PsozClhkTqhVYUiEir0fRSZlLBIraWlQXy87wKEavkuIsGUllb5amZ9H4U3JSxSJ371q9eBQ+W2aQWziDjBxypnXK48Zs0qdCYgqRNKWKTWvvzyS1566UqSkuYSH2+1gllEHFVmlTPGQMuWh/B4xvPjjw84HZrUghIWOX5uN574eLr37MmWoiLeG9+QjAyjFcwi4rjiVc54PJCdfQKjR9dj86xZHD7lFHV0DlNKWOT4FK9jdm3bhgto7/HQ9Pe/1xeAiISkJ849l7lA3E8/obXO4UndmqVm1B1VRMKJvrNCnro1S0BZrRsUkXCi76ywp27NUjNpaRQVFbGncWNO9vVXibqjikgoUUfniKERFqmxWbNmcVt+PgWxseVv0DpmEQlVPtY65xnD4T/+0aGApKaUsIhfyjc2HEfW4Ceo9/TT6sQsIuGhwlrnvNatGW8H0OKmy7RoKExo0q1Uq6SxYdleQQ0bWubNM8pPRCQsud0wbtwRCgqOjhQ3aqS/u5zg76RbJSxSKX8m1w8YoFLXIhJekpP1vRZKtEpI6szhw76TWk2uF5FwpUVD4UerhKRSJX9hnHTSQfbvb3rM7WokJiLhKC2t8kVD8fGWtDQT7JDEDxphkWMdnWFLbuvWDNs/jJiY8kuYtSBIRMKZrwaJcIhxceOwCQkq3x+ClLBIeSUzbDMywFoaZWcz3/UJ7173ghYEiUjEqNggMSEBHul7D9N+WITJzFT5/hCkSbdylGaiiUg0Kl5hYNeswej7L+g06VaOj2aiiUiU8pmsgL7/QoQm3cpRaWl44uNxbdt27G2aYSsikcqf8v36/nOcRlikXBXb6366i0PElN9BM2xFJBr4mImbC/znV9NK1iFoHq6DlLBEufJzbA2LCm5gsutpclpohq2IRJkKM3ELTj2VqXHnc9nz40rWIWgeroM06TZK+VPF1lczZhGRaFDyHbl6dfny/SX0HVl3NOlW/KIqtiIilfOVrIC+I52gSbdRKi0NPB4PzZrtJyenxTG3JyQEPyYRkVDhzzxcCS6NsESb4hm21uViX9OmXJZzGfXqHSm3i+bYioh4+aqIew0L+OqX9pqFG2RKWKJJmRm2xlpaHjrEwvrr+GDiUlWxFRHxoWJF3BuaLmIuE2m6P0uzcINMk26jharYiojUjr5HA0KTbuUYVlUcRURqR9+jjlHCEuFKi8Kt+IBM2vreSVUcRUSql5ZW6Wzb7EatSNj6gaa1BJASlghWsSjcXfyZQ1SYPaYZtiIi/vMxCzfPVZ9bch8kM9NoWksAKWGJQMnJ3sv48Zbc3KPbnyOFicwlA82wFRE5LmVm4Xow/BSXwGSzgOe4rtxuubn8f3t3H1tXXcdx/P1ZuzLuMKjbwsMeWiSLuhARN+VpUQP8gQ9hhqhhKY8uLDCeFB+yjhADyZBEQ8QIM5sMIW0GiCRMYkRlgGCQMEACYxAGa/cAk20E9lDp6Pr1j3u63Xt3u3aj957T088raXrO757efnd2es73/h6ZNy+lGHPKCUueJO0/q54cw4pnpnB+z4Hp/QpaaaET+vqgs9PJipnZoWpthc5OzvpaHxec1sm9ew+8j86lg9d6Wjz0eRh54ri86G//6e5mDHDcns0sYz4gVlD+x+QJj8zMPr6BJpebSwfLmM94uiHY30YE/pD4MbiGJQ+K7T+Utf8A4/kft9JWVuYuK2Zmw6uyW8st3FBMVkr1txH1L1Jkh8wJywi2bwTQk6vo69lT9ZipbKT5iC3usmJmViNlk8vRxzQ2VD0uenpo+fd9biU6TJ44boQqaQECYD0ttDDAghednXWNzcxsVBtgAaIuptJSkswUCv4QCZ44Lrf2jwAqbwFaxGIPWTYzy4IqQ593cyRt/KKszCOJDo0TlpGibATQZM7vaS97uX/IcqeHLJuZpavK0OfLWXbAAIjiSKJmjyQaIjcJjQSV7T8Us/VqfwBuATIzy56DjiTqN0rbiNwkNILt60yrPlrGbWHnpVd5BJCZ2Qg21JFEuy69ipZxW1zpUoUTlowpm06fMXT1HMv43h1Vj/UIIDOzkWGoI4kKvTvo6jnWU/xX4YQlC5IqlT6NYfaFzczpLr86NzCt6o+puZnOD4/1pLVmZiNAMkEufTGGMc3V7+uV9/s53R3MvrDF/VxwwlJ3+5p7kmvv6QX7q1TGEDSzgWXMZy77L0qPADIzy5mqI4kKLGL/fb2/n0szXZRWuTy9oKPsOTJachgnLLVWkqHsmtjCPy7rSFZPLl57U5e0Vemf0s0t3LBvfwWttE3or0t0+4+Z2YhX1kZUvK+3TVhaNpBioH4u05a0lT1H5s9PPvzmPItxwjKMDlZ7QgRHbe/itx+V155MZVPV95pWMbnQqbf31yW6/cfMLBday+/rp97eWlbpMlA/lykVz4053R2csmT/s4auLnp/MJ9rJ3bkKn+pacIi6VxJr0taJ2lhldePkHR/8vqzklpqGc+wG6T2ZNoQak8G6p/ydsM0V6aYmY0ilZUubzcMrZ9LtZqYxj3dXL/9hnzVwkRETb6ABuBN4DNAE/ASMKPimAXA75LtC4D7B3vfmTNnRs21t0c0N0dIxe/t7QcUPXVle0ShEFHMTyIgdlGIubTvK9qLyl7v/9qL9u3OpT12q/x9olAoxmBmZqNX+4HPmd0qf84cyrNmF+Xv9VFTIa6Z0F76qKv6/Ks1YHUMJa8YykGH8wWcDjxast8GtFUc8yhwerLdCGwjmcxuoK+aJyxVLpCexnFxUcO9ZddCJ1OrXiDrad63u57mQY8pFJLkp84XiJmZjQAVCcRTV7ZXPqKiS83D8jy6uOHe6GkcV/cP0ENNWGrZJDQZ2Fiyvykpq3pMRPQCHwATahjTwVVbpAdo6v2Qm/feWFY2lL4ni1hMt8p7gfc2FbhtwuKy5p7Zd7p/ipmZVVHRz2X2na2VfXXZcMXgI44G6g9TWn7T3htp6v2w/IAMLXg0IjrdSpovabWk1Vu3bq3tL+vpqVpc+Z89UN+T0vKHC628cEX5ldW4fCm/2dbq3MTMzA5LRQ5T/NBbksXsmtDM1WPLRxxt0uDPrIGSmoGei/VWy4RlMzC1ZH9KUlb1GEmNwNHA9so3ioilETErImZNmjSpRuECTzxR/A+vojJBce2JmZllRkkWc9S2Ts65u/WQa2EG+iA+0HOx3mqZsDwHTJd0gqQmip1qV1YcsxK4JNn+LrAqac9KT5XJfHqbCtw0tnySNteemJlZVh1OLcxNYxfT25TdSUprlrAkfVKuptixdi3wQESskXSzpPOSw+4CJkhaB1wPHDD0ue6qTObTuHzpAdmqa0/MzGxEGaQW5py7W2lcnt1JSpV2hcahmjVrVqxevTrtMMzMzGwYSHo+ImYNdtyI6HRrZmZmo5sTFjMzM8s8JyxmZmaWeU5YzMzMLPOcsJiZmVnmOWExMzOzzHPCYmZmZpnnhMXMzMwyzwmLmZmZZZ4TFjMzM8s8JyxmZmaWeU5YzMzMLPNG3OKHkrYCXXX6dROBbXX6XaOZz3Pt+RzXh89z7fkc10c9z3NzREwa7KARl7DUk6TVQ1lB0j4en+fa8zmuD5/n2vM5ro8snmc3CZmZmVnmOWExMzOzzHPCcnBL0w5glPB5rj2f4/rwea49n+P6yNx5dh8WMzMzyzzXsJiZmVnmOWEZgKRzJb0uaZ2khWnHk0eSpkp6XNKrktZIui7tmPJKUoOkFyU9knYseSTpk5IelPSapLWSTk87pjyS9KPkXvGKpBWSxqUdUx5IWi7pXUmvlJR9WtLfJb2RfP9UmjGCE5aqJDUAdwDfAGYAcyXNSDeqXOoFfhwRM4DTgKt8nmvmOmBt2kHk2O3AXyPic8DJ+FwPO0mTgWuBWRFxEtAAXJBuVLnxB+DcirKFwGMRMR14LNlPlROW6r4CrIuItyJiD3AfMCflmHInIt6JiBeS7Z0Ub/KT040qfyRNAb4F/D7tWPJI0tHAV4G7ACJiT0S8n25UudUIHCmpESgAb6ccTy5ExD+B9yqK5wD3JNv3AN+pa1BVOGGpbjKwsWR/E36Q1pSkFuAU4Nl0I8mlXwM/A/rSDiSnTgC2AncnzW6/lzQ+7aDyJiI2A78CNgDvAB9ExN/SjSrXjomId5LtLcAxaQYDTlgsAyQdBfwJ+GFE7Eg7njyR9G3g3Yh4Pu1YcqwR+BKwJCJOAXaTgerzvEn6UMyhmCAeD4yXdGG6UY0OURxOnPqQYics1W0GppbsT0nKbJhJGksxWemIiIfSjieHzgTOk9RJsWnzLEnt6YaUO5uATRHRXzv4IMUExobXOcD6iNgaER8BDwFnpBxTnv1X0nEAyfd3U47HCcsAngOmSzpBUhPFjl0rU44pdySJYrv/2oi4Le148igi2iJiSkS0ULyOV0WEP5UOo4jYAmyU9Nmk6Gzg1RRDyqsNwGmSCsm942zcubmWVgKXJNuXAA+nGAtQrMq0ChHRK+lq4FGKPdGXR8SalMPKozOBi4CXJf0nKVsUEX9JMSazw3EN0JF8wHkLuCzleHInIp6V9CDwAsURhi+SwdlYRyJJK4CvAxMlbQJ+DtwKPCBpHtAFfD+9CIs8062ZmZllnpuEzMzMLPOcsJiZmVnmOWExMzOzzHPCYmZmZpnnhMXMzMwyzwmLmdVNsqrxgmT7+GSYqpnZoDys2czqJlkz6pFktV0zsyHzxHFmVk+3AicmEwW+AXw+Ik6SdCnF1WDHA9MpLnLXRHFiwR7gmxHxnqQTgTuASUA3cHlEvFb/f4aZ1ZubhMysnhYCb0bEF4GfVrx2EnA+8GVgMdCdLCb4DHBxcsxS4JqImAn8BLizLlGbWepcw2JmWfF4ROwEdkr6APhzUv4y8IVkVe8zgD8Wl5IB4Ij6h2lmaXDCYmZZ0VOy3Vey30fxXjUGeD+pnTGzUcZNQmZWTzuBTxzOD0bEDmC9pO9BcbVvSScPZ3Bmll1OWMysbiJiO/AvSa8AvzyMt2gF5kl6CVgDzBnO+Mwsuzys2czMzDLPNSxmZmaWeU5YzMzMLPOcsJiZmVnmOWExMzOzzHPCYmZmZpnnhMXMzMwyzwmLmZmZZZ4TFjMzM8u8/wM0BSFVsjpLfQAAAABJRU5ErkJggg==\n",
"text/plain": [
"<Figure size 648x504 with 1 Axes>"
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
}
],
"source": [
"plot_sampling_points(func, N=50, peg=\"gridpoints\")"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The [`Objective`](https://krotov.readthedocs.io/en/latest/API/krotov.objectives.html#krotov.objectives.Objective) class has two methods `mesolve` and `propagate` that simulate the dynamics with the sampling on the grid points, and sampling on the intervals, respectively. You may compare the two to estimate the discretization error.\n",
"\n",
"You may play around with your own choice of `func` and sampling points:"
]
},
{
"cell_type": "code",
"execution_count": 9,
"metadata": {},
"outputs": [
{
"data": {
"application/vnd.jupyter.widget-view+json": {
"model_id": "19215f3b638a4a19957a97ce9914982b",
"version_major": 2,
"version_minor": 0
},
"text/plain": [
"interactive(children=(IntSlider(value=10, description='N', min=5, step=5), Dropdown(description='peg', options…"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"from ipywidgets import interact, interactive, fixed, interact_manual\n",
"import ipywidgets as widgets\n",
"\n",
"%matplotlib notebook\n",
"\n",
"interact(\n",
" plot_sampling_points,\n",
" N=widgets.IntSlider(min=5,max=100,step=5,value=10),\n",
" T=fixed(10),\n",
" peg=['gridpoints', 'midpoints'],\n",
" func=fixed(func));"
]
}
],
"metadata": {
"hide_input": false,
"kernelspec": {
"display_name": "Python 3",
"language": "python",
"name": "python3"
},
"language_info": {
"codemirror_mode": {
"name": "ipython",
"version": 3
},
"file_extension": ".py",
"mimetype": "text/x-python",
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.7.0"
},
"toc": {
"base_numbering": 1,
"nav_menu": {},
"number_sections": true,
"sideBar": true,
"skip_h1_title": false,
"title_cell": "Table of Contents",
"title_sidebar": "Contents",
"toc_cell": false,
"toc_position": {},
"toc_section_display": true,
"toc_window_display": false
},
"varInspector": {
"cols": {
"lenName": 16,
"lenType": 16,
"lenVar": 40
},
"kernels_config": {
"python": {
"delete_cmd_postfix": "",
"delete_cmd_prefix": "del ",
"library": "var_list.py",
"varRefreshCmd": "print(var_dic_list())"
},
"r": {
"delete_cmd_postfix": ") ",
"delete_cmd_prefix": "rm(",
"library": "var_list.r",
"varRefreshCmd": "cat(var_dic_list()) "
}
},
"types_to_exclude": [
"module",
"function",
"builtin_function_or_method",
"instance",
"_Feature"
],
"window_display": false
}
},
"nbformat": 4,
"nbformat_minor": 2
}
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment