{ "cells": [ { "cell_type": "markdown", "id": "1ec21eec", "metadata": {}, "source": [ "# LQR: bringing a double integrator to rest\n", "\n", "This notebook solves a finite-horizon **linear quadratic regulator** (LQR) problem with\n", "`diffilqrax.lqr.solve_lqr`. It then checks the solution against the KKT optimality conditions\n", "and against the parallel solver `diffilqrax.plqr.solve_plqr`.\n", "\n", "**The problem.** The state is a position and velocity, $x = (p, v)$, and the input $u$ is an\n", "acceleration. The discretised dynamics with time step $\\delta t$ are\n", "\n", "$$\n", "x_{t+1} = A x_t + B u_t, \\qquad\n", "A = \\begin{pmatrix} 1 & \\delta t \\\\ 0 & 1 \\end{pmatrix}, \\quad\n", "B = \\begin{pmatrix} \\tfrac12 \\delta t^2 \\\\ \\delta t \\end{pmatrix}.\n", "$$\n", "\n", "We want to bring the state to rest at the origin using little control effort, so we minimise\n", "\n", "$$\n", "J = \\tfrac12 x_T^\\top Q_f x_T + \\sum_{t=0}^{T-1} \\tfrac12 \\left( x_t^\\top Q x_t + u_t^\\top R u_t \\right).\n", "$$" ] }, { "cell_type": "code", "execution_count": 1, "id": "161c6442", "metadata": { "execution": { "iopub.execute_input": "2026-10-14T22:41:11.394161Z", "iopub.status.busy": "2026-10-14T22:41:11.393860Z", "iopub.status.idle": "2026-10-14T22:41:12.039298Z", "shell.execute_reply": "2026-10-14T22:41:12.038432Z" } }, "outputs": [], "source": [ "from functools import partial\n", "\n", "import jax\n", "import jax.numpy as jnp\n", "import jax.random as jr\n", "import matplotlib.pyplot as plt\n", "\n", "from diffilqrax.lqr import solve_lqr, kkt\n", "from diffilqrax.plqr import solve_plqr\n", "from diffilqrax.typs import LQR, LQRParams, ModelDims\n", "\n", "jax.config.update(\"jax_enable_x64\", True)" ] }, { "cell_type": "markdown", "id": "d6f7392e", "metadata": {}, "source": [ "## Setting up the problem\n", "\n", "diffilqrax uses the general time-varying LQR cost\n", "\n", "$$\n", "\\ell_t(x, u) = \\tfrac12 x^\\top Q_t x + q_t^\\top x + \\tfrac12 u^\\top R_t u + r_t^\\top u + x^\\top S_t u,\n", "\\qquad x_{t+1} = A_t x_t + B_t u_t + a_t ,\n", "$$\n", "\n", "with terminal cost $\\tfrac12 x^\\top Q_f x + q_f^\\top x$. Every time-varying term is stacked along\n", "a leading time axis of length $T$. This problem is time-invariant, so we tile each matrix\n", "$T$ times." ] }, { "cell_type": "code", "execution_count": 2, "id": "5165ce90", "metadata": { "execution": { "iopub.execute_input": "2026-10-14T22:41:12.043265Z", "iopub.status.busy": "2026-10-14T22:41:12.042922Z", "iopub.status.idle": "2026-10-14T22:41:13.084892Z", "shell.execute_reply": "2026-10-14T22:41:13.084147Z" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "initial state: [ 1.81608667 -0.75488484]\n" ] } ], "source": [ "dims = ModelDims(n=2, m=1, horizon=100, dt=0.1)\n", "T, n, m = dims.horizon, dims.n, dims.m\n", "\n", "A = jnp.array([[1.0, dims.dt], [0.0, 1.0]])\n", "B = jnp.array([[0.5 * dims.dt**2], [dims.dt]])\n", "\n", "span_mat = partial(jnp.tile, reps=(T, 1, 1)) # (n, k) -> (T, n, k)\n", "span_vec = partial(jnp.tile, reps=(T, 1)) # (n,) -> (T, n)\n", "\n", "lqr = LQR(\n", " A=span_mat(A),\n", " B=span_mat(B),\n", " a=span_vec(jnp.zeros(n)),\n", " Q=span_mat(jnp.eye(n)),\n", " q=span_vec(jnp.zeros(n)),\n", " R=span_mat(jnp.eye(m)),\n", " r=span_vec(jnp.zeros(m)),\n", " S=span_mat(jnp.zeros((n, m))),\n", " Qf=jnp.eye(n),\n", " qf=jnp.zeros(n),\n", ")\n", "\n", "x0 = jr.normal(jr.PRNGKey(0), (n,))\n", "params = LQRParams(x0=x0, lqr=lqr)\n", "print(\"initial state:\", x0)" ] }, { "cell_type": "markdown", "id": "7153fc3c", "metadata": {}, "source": [ "## Solving\n", "\n", "`solve_lqr` runs a Riccati backward pass and then a forward rollout. It returns the optimal states\n", "$x_{0:T}$, the inputs $u_{0:T-1}$ and the costates (Lagrange multipliers) $\\lambda_{0:T}$." ] }, { "cell_type": "code", "execution_count": 3, "id": "1a31b041", "metadata": { "execution": { "iopub.execute_input": "2026-10-14T22:41:13.087831Z", "iopub.status.busy": "2026-10-14T22:41:13.087636Z", "iopub.status.idle": "2026-10-14T22:41:13.935907Z", "shell.execute_reply": "2026-10-14T22:41:13.935091Z" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "(101, 2) (100, 1) (101, 2)\n" ] } ], "source": [ "Xs, Us, lambdas = solve_lqr(params)\n", "print(Xs.shape, Us.shape, lambdas.shape)" ] }, { "cell_type": "code", "execution_count": 4, "id": "e7696b83", "metadata": { "execution": { "iopub.execute_input": "2026-10-14T22:41:13.939158Z", "iopub.status.busy": "2026-10-14T22:41:13.938909Z", "iopub.status.idle": "2026-10-14T22:41:15.479939Z", "shell.execute_reply": "2026-10-14T22:41:15.479349Z" }, "tags": [ "thumbnail" ] }, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAscAAAH/CAYAAAC7N3BzAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjguMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8g+/7EAAAACXBIWXMAAA9hAAAPYQGoP6dpAACCiElEQVR4nOzdeXxU1f3/8dedmcxk38geAgHZEdlBwAUFBRcUQcVKBZfa9ldwQ1ux39atWrStLUWtS2tdqtQdF1wpIggqqyyy7wlLEpKQTNZJMnN/f0wyEEggQMhNwvv5eNzeO+cu85mM1ndOzj3XME3TREREREREsFldgIiIiIhIc6FwLCIiIiJSTeFYRERERKSawrGIiIiISDWFYxERERGRagrHIiIiIiLVFI5FRERERKo5rC6gufP5fOzbt4+IiAgMw7C6HBERERE5CaZpUlRUREpKCjZb/f3DCsfHsW/fPtLS0qwuQ0REREQaQWZmJm3btq13v8LxcURERAD+H2RkZKTF1YiIiIjIyXC73aSlpQWyXX0Ujo+jZihFZGSkwrGIiIhIC3e8YbK6IU9EREREpJrCsYiIiIhINYXjZiYjr5RXluzE5zOtLkVERETkjKMxx82IaZo8MGctS7bl8fHa/Tw5/hw6JYRbXZaIiIgcxufzUVFRYXUZcoSgoCDsdvspX0fhuBkxTRh9djKrMwpYufsgl8/6hrtGdObnF3QkyK5OfhEREatVVFSwc+dOfD6f1aVIHaKjo0lKSjqlZ1MYpmnq7/fH4Ha7iYqKorCwsMlmq9hbUMZv31/Hwi0HAOiZEsmT48/h7NSoJnl/EREROZppmmRkZFBZWXncB0lI0zJNk9LSUnJycoiOjiY5OfmoYxqa6RSOj8OKcAz+L3nOD3t5dO4GCkorsdsMfnlhR+4c0RmX49T/ZCAiIiInprKykm3btpGSkkJUlDqsmqO8vDxycnLo0qXLUUMsGprp9CtPM2UYBuP6tWXePRdyea8kvD6TZxds5+pnlrBhn9vq8kRERM44Xq8XAKfTaXElUp/Q0FDA/4vMyVI4bubiI1z8Y2J/npvYjzZhTjZlFXH1s4t55qutVHk13klERKSpncp4Vjm9GuO7UThuIS7rlcwX91zApT0SqfSa/OXLLYx//ju25RRbXZqIiIhIq6Fw3ILEhbt44ab+/PX63kQEO1iTWcAVs77h34s1L7KIiIg0va+//hrDMCgoKDjmcenp6cycObNJajpVCsctTM1Y5C/vuYDzO8fhqfLx6NwNTH55GVmF5VaXJyIiImeQoUOHsn///sANiq+88grR0dFHHbd8+XJ+/vOfN3F1J0fhuIVKjgrhtVsH8YerexIcZOObrbmMmrmIT9but7o0EREROUM4nc4GzSscHx8fuFmuuVM4bsEMw+CmIenMveN8eqVGUVhWyZTZq5j29mrc5Sd/l6aIiIi0HsOHD2fq1KlMnTqVqKgo4uLi+P3vf0/NbL4HDx5k0qRJxMTEEBoaymWXXcbWrVsD5+/evZsxY8YQExNDWFgYPXv25NNPPwVqD6v4+uuvueWWWygsLMQwDAzD4OGHHwaOHlaRkZHB1VdfTXh4OJGRkVx//fVkZ2cH9j/88MP06dOH//znP6SnpxMVFcUNN9xAUVHRaf95KRy3Ap0Swnnv/w1l6kWdsBnw/qq9XDbzG5btzLe6NBERkVbLNE1KK6osWU70MRWvvvoqDoeDZcuW8fe//52//vWv/Otf/wLg5ptvZsWKFXz00Ud89913mKbJ5ZdfHpgObcqUKXg8HhYtWsS6det48sknCQ8PP+o9hg4dysyZM4mMjGT//v3s37+f++6776jjfD4fV199Nfn5+SxcuJB58+axY8cOJkyYUOu47du388EHHzB37lzmzp3LwoULeeKJJ07oc58MPT66lXA6bNw3qivDu8Zzz9urycwv44YXv+NXwztx18jOevy0iIhIIyur9NLjwS8see8Nj44i1NnwGJeWlsbf/vY3DMOga9eurFu3jr/97W8MHz6cjz76iCVLljB06FAA3njjDdLS0vjggw+47rrryMjIYPz48fTq1QuAjh071vkeTqeTqKgoDMMgKSmp3lrmz5/PunXr2LlzJ2lpaQC89tpr9OzZk+XLlzNw4EDAH6JfeeUVIiIiALjpppuYP38+jz/+eIM/98lQYmplBqTH8umd5zO+X1t8JjyzYBvXPvctO3NLrC5NRERELHLuuefWGhc8ZMgQtm7dyoYNG3A4HAwePDiwr02bNnTt2pWNGzcCcOedd/LYY48xbNgwHnroIdauXXtKtWzcuJG0tLRAMAbo0aMH0dHRgfcE/1CMmmAMkJycTE5Ozim9d0Oo57gViggO4qnre3NRt3h++/461uwp5IpZ3/DwmJ5cN6CtJi8XERFpBCFBdjY8Osqy924qP/vZzxg1ahSffPIJX375JTNmzOCpp57ijjvuOK3vGxQUVOu1YRj4fKf/AWjqOW7Frjwnhc/vvoDBHWIprfDym/fW8qs3VlFQWmF1aSIiIi2eYRiEOh2WLCfa0bV06dJar7///ns6d+5Mjx49qKqqqrU/Ly+PzZs306NHj0BbWloav/zlL3n//fe59957+ec//1nn+zidzsBjtuvTvXt3MjMzyczMDLRt2LCBgoKCWu9pFYXjVi4lOoTZt5/L/aO74bAZfPZjFqNnfsO323KtLk1ERESaSEZGBtOmTWPz5s3897//5emnn+auu+6ic+fOXH311dx+++0sXryYNWvW8NOf/pTU1FSuvvpqAO6++26++OILdu7cyapVq1iwYAHdu3ev833S09MpLi5m/vz55ObmUlpaetQxI0eOpFevXkycOJFVq1axbNkyJk2axIUXXsiAAQNO68+hIRSOzwB2m8H/G34Wc341jI5xYWS5y5n40lJmfLqRiqrT/+cJERERsdakSZMoKytj0KBBTJkyhbvuuivwUI6XX36Z/v37c+WVVzJkyBBM0+TTTz8NDGvwer1MmTKF7t27M3r0aLp06cI//vGPOt9n6NCh/PKXv2TChAnEx8fzpz/96ahjDMPgww8/JCYmhgsuuICRI0fSsWNH3nrrrdP3AzgBhnmic4GcYdxuN1FRURQWFhIZGWl1OaestKKKP8zdyH+XZQBwdmokMyf0pVPC0VOyiIiIyCHl5eXs3LmTDh06EBwcbHU5DTZ8+HD69OnTYh7ffCqO9R01NNOp5/gME+p0MGNcL164qT8xoUH8uNfNlU9/wxtLd5/wnIkiIiIirY3C8RlqVM8kPr/7As7vHEd5pY//m/Mjt7+2krxij9WliYiIiFhG4fgMlhgZzKu3DOJ3V3THabfxv43ZjJr5DQs2n/45BEVERKRpfP3112fEkIrGonB8hrPZDH52fkc+mDKMLonh5BZ7uOXl5Tz04Y+UVx57KhYRERGR1kbhWADokRLJR1PP4+ah6QC8+t1urnx6Mev3FVpbmIiIiEgTUjiWgOAgOw9f1ZNXbx1EfISLbTnFjH12Cc99vR2vTzfriYiISOuncCxHubBLPJ/fdT6X9kik0mvy5Oeb+MmL35OZf/RE3iIiIiKticKx1KlNuIsXburPn649hzCnnWW78rns79/wzopMTfkmIiIirZbCsdTLMAyuH5DG53dfwMD0GIo9Vfz63bX88nVN+SYiIiKtk8KxHFdabChv/nwIvxndlSC7wRfrsxk1cxFfrM+yujQRERFpAoZh8MEHHzTb6zUmhWNpELvN4FfDOzHnV8PomhhBbnEFv/jPSqa9vZrCskqryxMREZEWZP/+/Vx22WUA7Nq1C8MwWL16tbVFVVM4lhNydmoUH90xjP83/CxsBry/ai+j/raIRVsOWF2aiIiItBBJSUm4XC6ry6iTwrGcMJfDzv2ju/HOL4fSIS6MLHc5k/69jP+bs45iT5XV5YmIiMhhXnzxRVJSUvD5fLXar776am699VYAPvzwQ/r160dwcDAdO3bkkUceoaqq/v+mr1u3josvvpiQkBDatGnDz3/+c4qLi2sd8+9//5uePXvicrlITk5m6tSpgX2HD6vo0KEDAH379sUwDIYPH86iRYsICgoiK6v2EM67776b888//6R/Fg2hcCwnrX/7GD6589CDQ95YmsGovy1i8dZcawsTERFpCqYJFSXWLCcwc9R1111HXl4eCxYsCLTl5+fz+eefM3HiRL755hsmTZrEXXfdxYYNG3jhhRd45ZVXePzxx+u8XklJCaNGjSImJobly5fzzjvv8L///a9W+H3uueeYMmUKP//5z1m3bh0fffQRnTp1qvN6y5YtA+B///sf+/fv5/333+eCCy6gY8eO/Oc//wkcV1lZyRtvvBEI9KeL47ReXVq9UKeDh6/qyaU9EvnNe2vZc7CMn760lJ8MSuOBy7sTGRxkdYkiIiKnR2Up/DHFmvf+7T5whjXo0JiYGC677DJmz57NiBEjAHj33XeJi4vjoosu4tJLL2X69OlMnjwZgI4dO/KHP/yB3/zmNzz00ENHXW/27NmUl5fz2muvERbmr+GZZ55hzJgxPPnkkyQmJvLYY49x7733ctdddwXOGzhwYJ31xcfHA9CmTRuSkpIC7bfddhsvv/wyv/71rwH4+OOPKS8v5/rrr2/Q5z5Z6jmWRjG0Uxxf3H0Bk4e0B+C/yzIZ9bdFfL05x+LKREREZOLEibz33nt4PP6pWN944w1uuOEGbDYba9as4dFHHyU8PDyw3H777ezfv5/S0qMfALZx40Z69+4dCMYAw4YNw+fzsXnzZnJycti3b18giJ+sm2++mW3btvH9998D8Morr3D99dfXet/TQT3H0mjCXA4eufpsLuuVzP3vrWV3Xik3v7yc8f3a8vsruxMd6rS6RBERkcYTFOrvwbXqvU/AmDFjME2TTz75hIEDB/LNN9/wt7/9DYDi4mIeeeQRxo0bd9R5wcHBJ1xaSEjICZ9Tl4SEBMaMGcPLL79Mhw4d+Oyzz/j6668b5drH0urD8aJFi/jzn//MypUr2b9/P3PmzGHs2LFWl9WqnduxDZ/ddT5/+WILL3+7k/dW7WHhlhweHNOTMeckYxiG1SWKiIicOsNo8NAGqwUHBzNu3DjeeOMNtm3bRteuXenXrx8A/fr1Y/PmzfWOCT5S9+7deeWVVygpKQn04i5ZsgSbzUbXrl2JiIggPT2d+fPnc9FFFx33ek6nv/PM6/Uete9nP/sZP/nJT2jbti1nnXUWw4YNa+hHPmmtflhFSUkJvXv35tlnn7W6lDNKqNPBg2N68O4vh9I5IZzc4gru/O8P3PbqCvYWlFldnoiIyBln4sSJfPLJJ/z73/9m4sSJgfYHH3yQ1157jUceeYT169ezceNG3nzzTX73u9/Ve53g4GAmT57Mjz/+yIIFC7jjjju46aabSExMBODhhx/mqaeeYtasWWzdupVVq1bx9NNP13m9hIQEQkJC+Pzzz8nOzqawsDCwb9SoUURGRvLYY49xyy23NOJPo36tPhxfdtllPPbYY1xzzTVWl3JG6t8+hrl3nsc9I7vgtNv4alMOl/51Ia8s2YnX1/A7bUVEROTUXHzxxcTGxrJ582ZuvPHGQPuoUaOYO3cuX375JQMHDuTcc8/lb3/7G+3bt6/zOqGhoXzxxRfk5+czcOBArr32WkaMGMEzzzwTOGby5MnMnDmTf/zjH/Ts2ZMrr7ySrVu31nk9h8PBrFmzeOGFF0hJSeHqq68O7LPZbNx88814vV4mTZrUSD+JYzNM8wTmAmnhDMM47rAKj8cTGKwO4Ha7SUtLo7CwkMjIyCaosvXallPE9PfWsWL3QQB6t43i8Wt6cXZqlMWViYiIHF95eTk7d+6kQ4cOJzUWV07ObbfdxoEDB/joo4+Oe+yxviO3201UVNRxM12r7zk+UTNmzCAqKiqwpKWlWV1Sq9EpIYK3fzGEP4w9mwiXgzV7CrnqmcU8/NF63OV6BLWIiIgcUlhYyOLFi5k9ezZ33HFHk72vwvERHnjgAQoLCwNLZmam1SW1KjabwU3ntmf+vRdyVe8UfCa88u0uRj61kI/X7OMM+kOGiIiIHMPVV1/NpZdeyi9/+UsuueSSJnvfVj9bxYlyuVzN9lnfrUlCZDCzftKX6wek8fsPf2Rnbgl3/PcH3l6RyUNjetIpIdzqEkVERMRCTTFtW13UcyyWOq9zHJ/ddb7/hj2HjW+25jJ65iL++OlGijTUQkRERJpYqw/HxcXFrF69mtWrVwOwc+dOVq9eTUZGhrWFSUBwkJ27Rnbmy7svYES3BKp8Ji8u2sHFTy3kvZV78GlWCxEREWkirT4cr1ixgr59+9K3b18Apk2bRt++fXnwwQctrkyOlB4Xxks3D+TlmweS3iaUA0Ue7n1nDdc+/y1r9xRYXZ6IiAiA7o9pxnw+3ylf44yayu1kNHTaD2lcniov/168i6e/2kpphf+JOeP6pfLrUV1Jjmqcx1KKiIicCK/Xy9atWwkNDSU+Pl5PfG1GTNOkoqKCAwcO4PV66dy5MzZb7T7ghmY6hePjUDi2VlZhOU98tpEPVvufXR8cZOPn53fkFxeeRZhL95OKiEjTKi4uZs+ePeo9bqZCQ0NJTk4OPJL6cArHjUThuHlYnVnA459sYPku/wNE4iNc3HdpF67tn4bdpt/cRUSk6Xi9XiorddN4c2O323E4HPX26CscNxKF4+bDNE0+/zGLJz7fxO68UgC6JIbz61HdGNk9QX/eEhERkXq1+CfkFRQU8K9//YsHHniA/Px8AFatWsXevXstrkysYhgGl/VK5st7LuB3V3QnMtjBluxibn9tBdc+/x3LduZbXaKIiIi0cM2y53jt2rWMHDmSqKgodu3axebNm+nYsSO/+93vyMjI4LXXXmuyWtRz3HwVllby/KLtvLxkJ+WV/rtTL+oaz29Gd6N7sr4rEREROaRF9xxPmzaNm2++ma1btxIcHBxov/zyy1m0aJGFlUlzEhUaxP2ju7Hw1xcxcXA77DaDBZsPcPmsb5gyexVbs4usLlFERERamGYZjpcvX84vfvGLo9pTU1PJysqyoCJpzhIjg3n8ml78b9qFXHlOMqYJn6zdz6UzF3HHf39gW45CsoiIiDRMswzHLpcLt9t9VPuWLVuIj4+3oCJpCTrEhfHMjf347K7zGd0zCdOEj9fs45K/LeKuN39g+4Fiq0sUERGRZq5ZhuOrrrqKRx99NDBNimEYZGRkcP/99zN+/HiLq5PmrntyJM/f1J9P7jyPS3skYprw4ep9jPzrQqa8sYof9xZaXaKIiIg0U83yhrzCwkKuvfZaVqxYQVFRESkpKWRlZTFkyBA+/fRTwsLCmqwW3ZDX8v24t5CZ/9vK/zZmB9ou7BLPlIs6MahDrIWViYiISFNpFfMcL1myhDVr1lBcXEy/fv0YOXJkk9egcNx6bMpy89zX2/l4zT581f/UD2gfw/8bfhYXdU3ApoeJiIiItFotOhy/9tprTJgwAZfLVau9oqKCN998k0mTJjVZLQrHrc/uvBKeX7iD91buocLrnwKuY3wYt53XgXF92xLitFtcoYiIiDS2Fh2O7XY7+/fvJyEhoVZ7Xl4eCQkJeL3eJqtF4bj1ynaX89Linfx3aQZFnioAYkKDuOnc9vx0SHsSIoKPcwURERFpKVp0OLbZbGRnZx81M8WaNWu46KKLAk/MawoKx61fsaeKt5Zn8u/FO9lbUAaA027jynOSuWlIe/qkRevR1CIiIi1ciwzHffv2xTAM1qxZQ8+ePXE4HIF9Xq+XnTt3Mnr0aN5+++0mq0nh+MxR5fXx5YZs/vnNDn7IKAi090qNYtKQ9ozpnUJwkIZciIiItEQNzXSOevdYYOzYsQCsXr2aUaNGER4eHtjndDpJT0/XVG5y2jjsNi7vlczlvZJZnVnAa9/tYu6a/azbW8iv313L459u5PoBadwwMI2O8eHHv6CIiIi0OM2q57jGq6++yoQJE2o9Otoq6jk+s+UVe3h7xR5e/353YMgFwKAOsdwwMI3Lzk7WDXwiIiItQIscVtEcKRwLgNdnsmBTDrOXZfD15pzAVHARwQ7G9knl+gFpnJ0aqbHJIiIizVSLDsder5e//e1vvP3222RkZFBRUVFrv27IEyvtLyzj3RV7eGtFJnsOHupN7pwQzjX9UhnbJ5WU6BALKxQREZEjNTTTNcvHRz/yyCP89a9/ZcKECRQWFjJt2jTGjRuHzWbj4Ycftro8OcMlR4Vwx4jOLPr1Rbx+22CuOCcZp8PG1pxi/vT5ZoY9+RU3/vN73lmRSVF5pdXlioiIyAlolj3HZ511FrNmzeKKK64gIiKC1atXB9q+//57Zs+e3WS1qOdYGqKwrJLP1u3n/R/2smznob9sOB02hneJ58reKYzolkCYq1ndAysiInLGaNHDKsLCwti4cSPt2rUjOTmZTz75hH79+rFjxw769u1LYWFhk9WicCwnKjO/lA9+2Muc1XvZcaAk0B4cZOPibglceU4Kw7vGE+pUUBYREWkqLXIqtxpt27Zl//79tGvXjrPOOosvv/ySfv36sXz58qMeKS3S3KTFhnLHiM5MvbgTm7KKmLt2H3PX7md3Ximfrsvi03VZuBw2zu8cx6U9kxjRLYE24frnWkREpDlolj3H06dPJzIykt/+9re89dZb/PSnPyU9PZ2MjAzuuecennjiiSarRT3H0hhM02T9Pjcfr93HZ+uyyMgvDeyzGTAgPZZLeyRyUbcEOsaFadYLERGRRtaih1Uc6fvvv+fbb7+lc+fOjBkzpknfW+FYGptpmmzOLuLL9dl8sT6L9fvctfa3iw3l4m4JDO8az7kd2+ipfCIiIo2gRYfjRYsWMXTo0FqPjwaoqqri22+/5YILLmiyWhSO5XTbc7CUL9dn89WmHJbuzKPSe+hfyeAgG+d2bMN5neI4r3McXRMj1KssIiJyElp0OLbb7ezfv5+EhIRa7Xl5eSQkJOD1epusFoVjaUrFniq+3ZbLgs05LNh0gCx3ea39ceEuzuvUhmGd4hhyVhvaxoRaVKmIiEjL0qJvyDNNs87esby8PMLCwiyoSKRphLscXNoziUt7JmGaJpuyili8NZfF23JZujOP3GIPH6zexwer9wGQGh3C4I6xnNuhDYM6xNK+Tah6lkVERE5Bs+o5HjduHAAffvgho0ePrjUzhdfrZe3atXTt2pXPP/+8yWpSz7E0F54qL6t2F7Bkmz8sr9tbiNdX+1/fxEgX/dvH0K9dDP3bx9AzJQqno1k+60dERKRJtcie46ioKMDfcxwREUFIyKFH8DqdTs4991xuv/12q8oTsZTLYWfIWW0YclYb7hvVlRJPFasyDrJ0Rz5Ld+axOrOAbLcnMF0c+B9Cck5qFH3bRXNO22h6t40mLTZEvcsiIiL1aFbh+OWXXwYgPj6ehx9+mNBQ/3jKXbt28cEHH9C9e3fi4uKsLFGk2QhzOTi/czznd44HoLzSy+rMAlZlHGTV7oOs3H2Qg6WVrNh9kBW7DwbOiwkNolfbaM5JjeLs1Ch6pkTSNkaBWUREBJrZsIoal1xyCePHj+eXv/wlBQUFdOvWjaCgIHJzc/nrX//K//t//6/JatGwCmmpTNNkV14pK3bls2ZPAev2FLJxfxEVXt9Rx0YEO+iRHEmPlEh6JEfSLSmSTgnhhDg1jZyIiLQOLXq2iri4OBYuXEjPnj3517/+xdNPP80PP/zAe++9x4MPPsjGjRubrBaFY2lNPFVeNmcVsWZPIWszC9iw382W7KJa08fVMAz/nMtdEiPomhhB58RwzooPp2N8mB59LSIiLU6LHHNco7S0lIiICAC+/PJLxo0bh81m49xzz2X37t0WVyfScrkcds5p6x9/zLntAaio8rH9QDHr97nZsM/Nhv2FbMkuJr+kgt15pezOK2Xehuxa10mNDqFjfFggLKe38S8p0cE47LoBUEREWq5mGY47derEBx98wDXXXMMXX3zBPffcA0BOTo56b0UamdNho3tyJN2TI6H/ofbcYg9bsorYnF3E5qwituUUs/1AMQdLK9lbUMbegjK+2Zpb61oOm0FabCjt24SSFhNKWmwIbWMObUeFBGlss4iINGvNMhw/+OCD3Hjjjdxzzz2MGDGCIUOGAP5e5L59+1pcnciZIS7cRVwnF0M71b4JNr+kgh0H/EF5+4ESduaWsDuvhN15pXiqfOzM9bfVJdzlIDkqmJToEFKig0mJCiE5OoSUqGASIoNJjHQRERzUFB9PRESkTs1yzDFAVlYW+/fvp3fv3ths/j/TLlu2jMjISLp169ZkdWjMsUjD+Hwm2UXl7MotZXdeCZkHS8nML2PPwVIyD5ZxoMjToOuEOe0kRgaTEOkiISKYuHAX8RH+JS7cSXyEizZhLmLDnJrDWUREGqxF35DX2J599ln+/Oc/k5WVRe/evXn66acZNGhQg85VOBZpHOWVXvYWlLG/oJx9BWXsKyzzrwvKyXKXk+0up6i86oSuGRHsoE2Yk9gwJ7FhLqJDg4gJDSI61ElMqJPo0CCiQ4KIDAkiqnod4XJgs2loh4jImaZF35DXmN566y2mTZvG888/z+DBg5k5cyajRo1i8+bNJCQkWF2eyBkjOMjOWfH+GS/qU+KpIqfIQ1ZhOTlF5Rwo8viXYk9gO7fYQ35JBT4TisqrKCqvYldeaYPrMAyIcDmICA4iIthBZHAQ4cEOIoIdhLv8S1j1Eu6yB7ZDg/zbIU47oU47oU4HIUF2guyGxlGLiLQirb7nePDgwQwcOJBnnnkGAJ/PR1paGnfccQfTp08/7vnqORZpfnw+k8KySvJLK8gvqSCv2L8uKKugoLSSgyUVHCytpKC0gsKySgrLKnGXV1JeefQcz6fKbjMICbITHGQnOMgW2HY5bEetnQ7bYWs7zurtILsNp90IbB9aDILsNhzV6yC7DYfNwGE3/GubrXrbht1mBBZH9dpm1KxpmgBvmuDzguk9Yu079PqobZ9/Xe9iHt3GkW1mdZt56Pia18ddc8Q2h/YHtqnj9WFtR7XXsf/In9OJOqHvr55ja13DOM51DUxMf6nVx5imiXnYtWu2a30a87CVYRz6qIZRff4R7YHrmf73rDm/pi7TDGybh1/7iFoP/4yHjqujvpp3q/W5qPPnUOsfCQ47xqTWzwHjGMfWc71a7Ud+hnra/O1m3e0N+Ay1r92wOupj1l1Gvdeufwek9RqGKzisQe97qtRzDFRUVLBy5UoeeOCBQJvNZmPkyJF89913dZ7j8XjweA6NjXS73ae9ThE5MTabQUyYk5gwJ2fFN/w8T5UXd1kVhWWVFJVXBnqei8orKfZU4S6vosTjX4o9NdteSiqqKKvwUlpxaLvK5/9/e6/PpLj6+NPBwEcwFYRQgYtKXEb1umYxKnFRgZMqnFTiNKpwUYmTKoKq24KMKlyGFydVuIwqnIbXv696HYSXIMOLo3rbUbMYVThMH3a82Kvb7IHFh930b9vwBdqk9TBoaFQSOXn7E5aR3L6r1WXU0qrDcW5uLl6vl8TExFrtiYmJbNq0qc5zZsyYwSOPPNIU5YlIE3M57MRH2ImPcJ3ytTxVXsorfZRXeinzVFFeVkRlaSFVpYV4y4vxlbvxeUowPUXgKYaKEozKUmxVpdiryrBXlWL3lmL3lhPkLcfhKyfIV06Qz4PTLMdpVuA0PQRxekL3UU6g4/Nk+UyjOkrb8FGzbdR6zRHtpmngw8Dk0PrwbV91hDv6GGq117TVbFOzbR7qXTx0Xu3jOKzt8PWRP6Ijj6+vrTEY9Xw59bcf+9j6zmvI8UYdx9TbH20c/z3rep+63qP2tlHHN1bH8cbRbbWPPdb7HH11A/Oo5uP9bOquz6ijrb5j69bQf9JO5hqn8w9PwY7mN0NRqw7HJ+OBBx5g2rRpgddut5u0tDQLKxKR087ng/ICKDtYvS6ovS4vhHI3eNz+dXkhLo8bl6eIKE8xVBRV/yn/NLMFQVAIOFzgqFkHg8MJdhem3RlYfHYnps2FzxaEaQvCZ3fiswXhMxyYNic+mwOv4d/ntTnwGnZMw7/fZziq91e32ez4sOMzqhccgW0z0GYLvMaw4zP8/cnmYeuaP63DEaGyrr/L16mOcGLUvbdmGIlxxHGHwpFR6yTjsPOOOgfjsG3//9S01TrPOPL9ap93eE2GUfu61NFW8/rQdu3r1lfzkdc4/JjjfabD36u+a9f6uR3n8x1+zOGfs6aGI69R8561Xx9qE2kKrTocx8XFYbfbyc6u/XSv7OxskpKS6jzH5XLhcp16r5KIWMQ0wVMEJQeql1wozatjyfeH4bJ8fwhujG5Sww6ucHBGVK/DwBlevYSBMxSCwo7YDvUH3qAj1o4QCAo+bB0MNvux355DgeLYR4qISH1adTh2Op3079+f+fPnM3bsWMB/Q978+fOZOnWqtcWJyImp8kDRfijKguJsKMr2r4uzoTjHvy7J9Qdib8PmVD5KUBiExEBINARHH1oHR0FwJLgij9iuXjvDwRXhD7Xq4RIRadFadTgGmDZtGpMnT2bAgAEMGjSImTNnUlJSwi233GJ1aSJSw1ME7n1QuAfcew9t14Rh9z5/D++JCAqDsDgIi4fQNv7t0FgIjfO/Do2FkNjqMFy9OJyn5/OJiEiL0erD8YQJEzhw4AAPPvggWVlZ9OnTh88///yom/RE5DQqOwgHd8HB3VCYCQWZh60z/GN6G8Lugogk/xKeAOGJhy0JEJZQHYjj/EMXRERETlCrn+f4VGmeY5EG8Pn8Pb75Ow4tB3f5l4LdDQu/riiISoXIFIhMrV5SIDIZIqqXkBgNWxARkZOieY5FpHGZpv9GttytkLe1er0N8rb7Q/DxxvmGJUB0u+olDaLSILp99XZb/5hdERERiykci0htPh+498CBzdXLJv86d4t/WrP62BwQkw6xHf1LTAeIae9vi26nYQ4iItIiKByLnMmKcyBnA2RvgJz1kLMRcjZBZUk9Jxj+Ht+4TtCmM7TpBG3O8i+RbcGu/0sREZGWTf8lEzkTeCv9wyCyf4SstZD1I2Stg9Lcuo+3BfmDb3zXQ0tcV38IDgpp2tpFRESakMKxSGtTWe7vBd6/xr/sW+3vEa5zTLABsR0goQck9oSE7v7t2I5gb36P9BQRETndFI5FWjJvpT/47l0J+1bB3h/gwEbwVR19rDPCH4CTekHS2ZDYyx+GnaFNX7eIiEgzpXAs0lKYpn9u4D3LYc8KfyDevwaqyo8+NiQWUvpAcm9I7gPJ50B0OthsTVy0iIhIy6JwLNJcVZbBvh8gc1l1IF7uf0TykVxR/iCc2h9S+/nDcFRbzQcsIiJyEhSORZqLklzI+B4yv/ev960GX2XtY2wO/7CItgMhdYA/EMd2VI+wiIhII1E4FrFKQSbs/hZ2L/Gv87YefUx4IqQNgraD/IE4pY9mixARETmNFI5FmoJp+p8it2uxf9n9LRRmHH1cfHdoNxjSzoV25/ofoKHhESIiIk1G4VjkdCnIhF3fwM5v/OvCzNr7Dbu/J7j9UGg/DNIGQ2isJaWKiIiIn8KxSGMpyYNdi2DHQti5EPJ31N5vc/jHCKef51/aDgJXuDW1ioiISJ0UjkVOVmUZZHwH2xfAjq/9T5zDPLTfsENKX38Q7nC+f6iEwrCIiEizpnAs0lA+n//xyzsW+ANxxndHzzGc0AM6XAgdL/QPlwiOsqZWEREROSkKxyLHUnzAH4a3zYft86HkQO39ESlw1sXQcTh0uAAiEi0pU0RERBqHwrHI4bxVsGcZbPufPxDvX117f1CYf5jEWRfDWRdBXBfNJiEiItKKKByLFO719wpvnee/mc5TWHt/Ui84awR0GumfUcLhtKZOEREROe0UjuXM462EzKWw9UvY+j/IWV97f0isv2e48yXQ8SINlRARETmDKBzLmcG9H7bNq+4d/ho87sN2Gv4p1jpfAp0u8c89bLNbVKiIiIhYSeFYWidvFexZXt07PA+y19XeH9rGP0yi0yX+XuKwNtbUKSIiIs2KwrG0HkXZ1TfSzYPtX0H54WOHDUjtB50v9fcQJ/cFm82yUkVERKR5UjiWlstbBXtX+HuGt82D/Wtq7w+JOdQ73GkEhMVZU6eIiIi0GArH0rLUzCyxbb5//uHyI2aWSO7j7xnufKl/HLHGDouIiMgJUDiW5q3mEc3b5vuHSuRsqL0/JMY/Zrimdzg8wZo6RUREpFVQOJbm5biPaK6eWaLTSH8YVu+wiIiINCKFY7GWaUL+Dtj1Dexc5F+OfERzeFL1vMMj/fMOh8ZaU6uIiIi0egrHzU1JLtidEBxpdSWnh2lCwW7Y/S3srA7E7j21jwkKg/Rh/iB81kUQ302PaBYREZEmoXDc3Cz+Gyz/F3S7As6Z4O8xtQdZXdXJ83khez1kfA8Z3/rXRftrH2MLgrYDocMF0OF8aDtIj2gWERERSygcNzf71/jH2P74nn8JjYOzx0PvCZDSr3n3oJomuPfC3pXVyyrYtxoqimofZwvyP4Uu/Tx/IE4bDM4wKyoWERERqcUwTdO0uojmzO12ExUVRWFhIZGRTTDUwTRh3w+w9i1Y9y6U5h7aF90e0s+H9kP9S0y6dWG5ygO5WyB7A+Ss96/3r4GSnKOPdUZA2iBoPwTaDfGHfGdo09csIiIiZ6yGZjqF4+No8nB8OG+lf8aGtW/Bpk+gqqz2/shUf0hOOgdiO/qXmPTGC55VHijIhIM74eAu/5K/E/K3Q+5WML1Hn2PYIbGnfxaJ1H7+dXw3zSghIiIillI4biSWhuPDeYr843V3L/HfzLZ3Ffgq6z42ItkfkkNiwBUBznD/2hUOjhD/ed7qpWa7vBBK8/w3BJbmQkkeeArrvn4NVxQk9oCEHv51Yi9IPgeCQhr944uIiIicioZmOo05bilcEdVPfrvE/7qi1P/o5N3f+Yc35O/w9+iWF/pveDvypreTFRTmD9ox6RDb4dB2Qnd/z3VzHgMtIiIicoIUjlsqZ2j17A4X1G4vzfcPfSjYDR63v8fZUwSeYv/rKo9/9gubwz9lXM22KwLC4vw3AB6+DolRABYREZEzhsJxaxMa61/a9re6EhEREZEWx2Z1ASIiIiIizYXCsYiIiIhINYVjEREREZFqGnN8HDUz3bndbosrEREREZGTVZPljjeLscLxcRQV+R99nJaWZnElIiIiInKqioqKiIqKqne/HgJyHD6fj3379hEREYHRBFOaud1u0tLSyMzMtPahI3LK9F22DvoeWwd9j62HvsvWwYrv0TRNioqKSElJwWarf2Sxeo6Pw2az0bZt2yZ/38jISP1L30rou2wd9D22DvoeWw99l61DU3+Px+oxrqEb8kREREREqikci4iIiIhUUzhuZlwuFw899BAul8vqUuQU6btsHfQ9tg76HlsPfZetQ3P+HnVDnoiIiIhINfUci4iIiIhUUzgWEREREammcCwiIiIiUk3hWERERESkmsKxiIiIiEg1PSHvOJr68dEiIiIi0vj0+OhGsm/fPtLS0qwuQ0REREQaQWZmJm3btq13v8LxcURERAD+H6Se4S4iIiLSMrndbtLS0gLZrj4Kx8dRM5QiMjJS4VhERESkhTveMFndkCciIiIiUk3hWERERESkmoZVSItV4a3gYPlBCisKKfQcWooriwmyBeGyu3DanTjtTlx2F+FB4SSFJZEYmkiQPcjq8kVERKQZUjiWFqHQU8im/E1szNvIxvyNbMrfxC73Lnym74SvZWAQHxJPUngSyWHJpEem0y22G11ju5IanorN0B9UREREzlQKx9IsmabJ1oKtfLnrS+btnseOwh11Hmc37ES5ovyL078OCwqjyldFhbcCj9dDha+CCm8FBZ4CskuyqfBVkFOWQ05ZDmsPrK11vbCgMLrGdKVrbFf6xPdhQNIAEkITmuIji4iISDNgmKZpWl1Ec+Z2u4mKiqKwsFCzVTSBbQe38cXuL/hi1xfsLNxZa19qeCrdY7vTLbYb3dv41/Eh8Sf0cBbTNMkvzyerJIv9JfvZV7yP7YXb2ZS/iW0Ht1HhqzjqnHYR7RiQNIABif4lOTz5lD+niIiINK2GZjqF4+NQOD79TNPk233f8o81/6jVkxtkC2JY6jAubX8pF7S9gChX1Gmto9JXya7CXWw+uJn1uetZmb2STfmbMKn9r0in6E6c3/Z8Lki9gD4JfXDY9AcYERGR5k7huJEoHJ8+dYXiIFsQw1KGcWn6pQxPG06E89gTdZ9u7go3q3NWsyJrBSuyV7A+b32tcc4RQREMTR3K8LThXNj2QsvrFRERkbopHDcShePGV1codtldXN/1em49+1biQuIsrrB+hZ5Cluxdwjd7v2Hx3sUUeAoC+4JsQQxNGcol7S9heNrw097TLSIiIg2ncNxIFI4bV15ZHo9+9yhfZX4FQLA9mOu7Xs8tZ9/SrENxXbw+L+ty17FozyL+l/G/WmOkHYaDwSmDGZ0+mpHtRhLuDLewUhEREVE4biQKx41nQcYCHv7uYfLL83HYHPyk20+afU/xidhesJ0vd/tn19h6cGug3WV3cWHbC7mi4xWcn3q+5lgWERGxgMJxI1E4PnXFFcX8afmfmLNtDuC/oW3G+TPoFtvN4spOn52FO/ly15d8svOTWj3Kkc5IRqWP4qqzrqJ3fO8TmmlDRERETp7CcSNROD41K7JW8Lslv2Nv8V4MDG7ueTNT+07FaXdaXVqTME2Tjfkb+WTHJ3y28zMOlB0I7EuPTOfqTlczpuMYEsMSLaxSRESk9VM4biQKxyfv7c1v8/jSx/GZPlLDU3n8vMfpn9jf6rIs4/V5WZ69nI+3f8y83fMoqyoDwGbYGJI8hKs7Xc3F7S7GZXdZXKmIiEjro3DcSBSOT5xpmjy/9nn+sfofAIzpOIb/O/f/CAsKs7iy5qOksoQvd33Jh9s/ZGX2ykB7lCuKMR3HcE3na+gS08XCCkVERFoXheNGonB8Yrw+LzOWzeCtzW8B8ItzfsGUPlM0tvYYMt2ZfLj9Qz7Y9gHZpdmB9l5xvRjXeRyXdbhMv1iIiIicIoXjRqJw3HAV3gqmfzOdebvnYWAwfdB0bux+o9VltRhen5dv933LnG1zWJCxgCqzCoAQRwiXd7ica7tcS882PfWLhoiIyElQOG4kCscNU1xRzF0L7mJZ1jIcNgczzp/B6PTRVpfVYuWW5TJ3+1ze2/oeu9y7Au1dY7pybZdruaLjFXoan4iIyAlQOG4kCsfHV1pZys++/BnrctcR6gjl7xf/nXOTz7W6rFbBNE1WZq/k3a3vMm/XPCp8FYD/4SmjO4zm2i7Xck7cOepNFhEROQ6F40aicHxsXp+Xu7++m68zvybaFc3zlzxPzzY9rS6rVSr0FPLx9o95d8u7bC/cHmjvEtOFa7tcy5Udr1RvsoiISD0UjhuJwvGxPbnsSV7f+DpOm5OXRr1En4Q+VpfU6pmmyeoDq3l3y7t8sesLPF4PoN5kERGRY1E4biQKx/V7Y+MbPLHsCQD+cuFfGJU+yuKKzjyFnkLm7pjLO5vfqbM3+YqOVxDp1D+3IiIiCseNROG4bl9nfs1dC+7CZ/q4u9/d3NbrNqtLOqMdqzd5VPooru1yrR5XLSIiZ7SGZjpbE9bUKJ599lnS09MJDg5m8ODBLFu2rN5j//nPf3L++ecTExNDTEwMI0eOPObx0jDr89bzm0W/wWf6GN95PLeefavVJZ3xDMOgb0JfHj/vceZfN5/pg6bTKboT5d5yPtz+ITd9dhPjPhrH6xtep9BTaHW5IiIizVaL6jl+6623mDRpEs8//zyDBw9m5syZvPPOO2zevJmEhISjjp84cSLDhg1j6NChBAcH8+STTzJnzhzWr19Pampqg95TPce17S/ez42f3khuWS5DU4byzIhnCLIFWV2W1ME0TdYcWMM7W97hy11fUu4tB8BpczKi/Qiu7XwtA5MGqjdZRETOCK1yWMXgwYMZOHAgzzzzDAA+n4+0tDTuuOMOpk+fftzzvV4vMTExPPPMM0yaNKlB76lwfEiVr4rJn09m7YG1dIruxGuXvabZEVoId4WbT3d8yntb32NT/qZAe1pEGtd0uoarzrqKxLBECysUERE5vVrdsIqKigpWrlzJyJEjA202m42RI0fy3XffNegapaWlVFZWEhsbW+8xHo8Ht9tdaxG/f//4b9YeWEt4UDjPjHhGwbgFiXRGckO3G3j7yrd584o3ua7LdYQFhZFZlMmsH2Zx6XuXMmX+FObvnk+lr9LqckVERCzjsLqAhsrNzcXr9ZKYWLt3KzExkU2bNtVzVm33338/KSkptQL2kWbMmMEjjzxySrW2Ruvz1vPc6ucA+O3g35Ia3rBhKdK8GIZBz7ie9IzryX0D7uPL3V8yZ+scVuWsYtGeRSzas4jY4Fiu7HglYzuNpXNMZ6tLFhERaVItpuf4VD3xxBO8+eabzJkzh+Dg4HqPe+CBBygsLAwsmZmZTVhl81ReVc5vv/ktVWYVl7S/hCs7Xml1SdIIQoNCGdtpLK9e9iofjf2IW86+hTbBbcgvz+e1Da8x7qNxTJg7gdkbZ+smPhEROWO0mJ7juLg47HY72dnZtdqzs7NJSko65rl/+ctfeOKJJ/jf//7HOeecc8xjXS4XLpfrlOttTf6+6u/sKNxBXEgcD577oG7gaoU6RHVgWv9p3NH3DpbsXcKH2z7k6z1fsyFvAxvyNvCXFX9heNpwrjrrKoalDtNNmCIi0mq1mHDsdDrp378/8+fPZ+zYsYD/hrz58+czderUes/705/+xOOPP84XX3zBgAEDmqja1uP7/d/z+sbXAXh06KNEB0dbW5CcVkG2IIanDWd42nAOlh/k052f8uG2D9mYv5F5u+cxb/c8YlwxjO4wmqvOuoqebXrqlyUREWlVWtRsFW+99RaTJ0/mhRdeYNCgQcycOZO3336bTZs2kZiYyKRJk0hNTWXGjBkAPPnkkzz44IPMnj2bYcOGBa4THh5OeHh4g97zTJ6twl3hZtyH48guzea6Ltfx4JAHrS5JLLI5fzMfbf+IT3Z8Ql55XqA9PTKdMWeN4bIOl5EWkWZhhSIiIsfWKqdyA3jmmWf485//TFZWFn369GHWrFkMHjwYgOHDh5Oens4rr7wCQHp6Ort37z7qGg899BAPP/xwg97vTA7HD3zzAHN3zCUtIo13x7xLaFCo1SWJxap8VXy//3s+2v4RCzIWBOZOBjgn/hwu73A5o9JHERcSZ2GVIiIiR2u14bipnanh+OvMr7njqzuwGTZeHf0qfRL6WF2SNDPFFcX8L+N/fLrjU5ZmLcVn+gCwGTYGJw3msg6XcXG7i4lyRVlcqYiIiMJxozkTw3Glt5JrPrqG3e7d3NLzFqYNmGZ1SdLM5Zbl8sWuL/h0x6eszV0baHfYHAxJHsLoDqO5KO0izY0tIiKWUThuJGdiOP7Phv/wp+V/ok1wGz4Z9wlhQWFWlyQtSIY7g892fsYXu79g68GtgfYgWxDDUoYxsv1IhqcNV4+yiIg0KYXjRnKmheOC8gIun3M5RRVFPDzkYcZ3GW91SdKC7SjYwRe7v+CLnV+wvXB7oN1hOBiYNJCR7UdycbuLNUZZREROO4XjRnKmheMZS2cwe9NsusR04e0r38Zus1tdkrQS2w5uY97uefwv439sObgl0G5g0CehDxelXcRFaReRHpVuXZEiItJqKRw3kjMpHO8o3MG4D8fhNb3889J/cm7yuVaXJK3Ubvdu5mfMZ/7u+bXGKIN/eriL0i7ionYXcU7cOfoFTUREGoXCcSM5k8LxlPlTWLRnEcPbDufpEU9bXY6cIbJKsvg682sWZC5gWdYyqnxVgX3RrmiGpQ7j/NTzGZYyTA+hERGRk6Zw3EjOlHD87b5v+cW8X+AwHMy5eo7+tC2WKK4oZvG+xXyd+TWL9iyiqKIosM9m2OgV18sflFOH0T22u3qVRUSkwRSOG8mZEI69Pi/Xzb2OrQe38tPuP+X+QfdbXZIIVb4q1hxYwzd7vuGbvd/UGqcM/l7lc5PPZWjKUIakDCEpLMmiSkVEpCVQOG4kZ0I4fmfLOzz63aNEOiP5dNynmmJLmqWskiwW713M4r2LWbp/KcWVxbX2d4jqwOCkwQxOHsyAxAEagiEiIrUoHDeS1h6OSytLuez9y8gvz+f+gffz0x4/tbokkeOq8lXxY+6PfLvvW5bsW8KPuT8GntAH/hkwusV2Y2DSQAYlDaJvYl8ina3v318REWk4heNG0trD8esbXufJ5U+SFpHGh2M/JMgWZHVJIies0FPIiuwVLN2/lGX7l9WaUxn8YblLTBcGJA2gf2J/+iX0o01IG4uqFRERKygcN5LWHI4rvZVcPudyskqyeHDIg1zX5TqrSxJpFLlluSzbv4ylWUtZmb2S3e7dRx3TPrI9veN70yehD33i+3BW9FnYDJsF1YqISFNQOG4krTkcf7jtQ3635HfEhcTx+fjPcdldVpckclrkluWyMntlYDny5j6A8KBwzok/h15xvegV14uecT315D4RkVakoZnO0YQ1STPiM338+8d/A3BTj5sUjKVViwuJY1T6KEaljwL8wzDWHljLmgNrWH1gNWsPrKW4sphv933Lt/u+DZyXEpbC2XFn0zOuJz3a9KB7bHfdsCoi0sqp5/g4WmvP8VcZX3HXgruICIrgy2u/JNwZbnVJIpap8lWxrWAba3LWsC53HT/m/siOwh2YHP1/j6nhqYGg3DW2K11jupIQmoBhGBZULiIiDaWeY6mXaZq8tO4lACZ0m6BgLGc8h81Bt9hudIvtxgQmAP4HkmzI28Da3LVsyNvAxryN7Cnew97ivewt3su83fMC50e7oukS04UuMV3oGtuVTtGd6BjVkdCgUKs+koiInCT1HB9Ha+w5Xp61nFu/uBWnzckX136hcZUiDVToKWRz/mY25G1gQ/4GtuRvYZd7F17TW+fxqeGpdIruRKfoTpwVfRYdojqQHpmuX0hFRCygnmOpV02v8TWdr1EwFjkBUa4oBiUPYlDyoECbx+thW8E2tuRvYcvBLWw+uJntBdvJL88P9DIv3LOw1nUSQhL8QTkqnfTIdNpFtqNdRDtSI1I1naKIiMUUjs8wG/M2smTfEmyGjck9J1tdjkiL57K76NmmJz3b9KzVnl+ez/aC7Wwr2Mb2gu1sL9jOLvcucstyySnLIacsh6VZS2udYzfspISn0C6yHWnhabSNaEvb8La0jWhLaniqepxFRJrASYXj1157jQkTJuBy1Z7hoKKigjfffJNJkyY1SnHS+GpmqBiVPoq0iDSLqxFpvWKDY4lNimVg0sBa7e4KN7sKd7GzcCc7CneQ4c5gd9FuMt2ZlHvLySzKJLMos85rRrmiSAlLITksmZTw2uvEsERig2M1V7OIyCk6qTHHdrud/fv3k5CQUKs9Ly+PhIQEvN66x9+1RK1pzHGGO4MxH4zBZ/p4d8y7dI3tanVJIlLNNE1ySnPIKMogw53hv/mvaC97ivewp2gPBz0Hj3sNh81BYmhiYEkITSA+NJ74kHjiQ+P9r0PidaOgiJyRTuuYY9M065y2aM+ePURFaQ7Q5uqV9a/gM32cl3qegrFIM2MYBolhiSSGJR7V2wxQUlnCnqI9ZJVksa9kH/uL9wfWWSVZHCg7QJWvKjDO+VhCHCHEhcTRJriNfx3ShjbBbYgJjiEmOMbf6x0cS0xwDFHOKOw2++n62CIizc4JheO+fftiGAaGYTBixAgcjkOne71edu7cyejRoxu9yMM9++yz/PnPfyYrK4vevXvz9NNPM2jQoHqPf+edd/j973/Prl276Ny5M08++SSXX375aa2xOSquKGbujrkA3Hr2rRZXIyInKiwozD+vcj2/2Fb6KsktzSW7NJus0iyyS7LJKc3hQNkBDpQe4EDZAXJKcyirKqOsquyYwzcOZ2AQ4Ywg2hXtX4L960hnJJGuSCKdkUS5ovyvq5dwZzgRzgiC7cGa/1lEWpwTCsdjx44FYPXq1YwaNYrw8EM3hzidTtLT0xk/fnyjFni4t956i2nTpvH8888zePBgZs6cyahRo9i8efNRQzwAvv32W37yk58wY8YMrrzySmbPns3YsWNZtWoVZ5999mmrszn6dOenlFWV0SGqAwMSB1hdjog0siBbEMnhySSHJx/zuJLKEvLK8sgtyyWvvHpdlkdeeR4Hyw9ysPwg+eX5HPQcpNBTiImJu8KNu8JNRlHGCdXksDmICIog3BlOeFA4YUFh/rXTvw4NCiXMEUZoUCihjlDCgvzbIY6QOheX3aWwLSKn3UmNOX711VeZMGECwcHBp6Omeg0ePJiBAwfyzDPPAODz+UhLS+OOO+5g+vTpRx0/YcIESkpKmDt3bqDt3HPPpU+fPjz//PMNes/WMub4hrk3sD5vPfcNuE+zVIhIg1T5qij0FFLoKeSg5yAFngIKPYWBtbvCjdvjprCiELfHH6CLK4spqijCZ/oavR4Dg2BHMCGOEILtwQQ7gnHZXf7F4SLYfui10+6stV3z2mFz+F/bnATZggL7HDYHQbagwFLz2mFzHLVtN+wE2YKw2+y6AVKkBTmtY44nT276cFVRUcHKlSt54IEHAm02m42RI0fy3Xff1XnOd999x7Rp02q1jRo1ig8++KDe9/F4PHg8nsBrt9t9aoU3A5vyN7E+bz0Om4MxZ42xuhwRaSEcNod/PHJImxM6zzRNyqrK/GG5opiiyiJKKksoriympMK/Lq4sprSylJLKEkqrSmtt1wz9KKsqo6yyjApfhf+6mIH25sJm2LAb9kBottvs/teGIxCeD2+veX34uq7FMAz/fvzbgXaMwP7ANkbgmJrtWuvqbeCo/TVtQO1j6mir3ggwDntx+HagrY5e/vr64+p6VHu97Wbd+w+/dq326u269te0HeuYuuoPHI951Gc6/Fp1bR/rWkcdG1jVXWddn//wmkxO/Fr1Hn+S16zZntxzMhe3u5iW4KTCsc1mO+aftk7HbBW5ubl4vV4SExNrtScmJrJp06Y6z8nKyqrz+KysrHrfZ8aMGTzyyCOnXnAz8u6WdwG4OO1iYoNjLa5GRFo7wzD8QyWCQiHs1K9X5auivKqccm+5f129XVZVRnlVOR6vJ7DUvK7wVgTWFb6KwOtKX2WgrcpbFdhX5aui0ld5aPFWUmVW+de+KqrMqjp7w32mD5/po9JXeeofVKQVu6LjFVaX0GAnFY7ff//9WuG4srKSH374gVdffbXFB8sHHnigVm+z2+0mLa3lzgdcVlXGpzs+BWB8l9M3HlxE5HRx2Bz+cctY+xAUn+nzB+XqsOz1eanyVeE1vVT6KvGZPn9bdZCu2VfT7jW9gdc1bT7T52/Dh8/nw4cP0zTxml5M0wzsr3l/n+nDxN9e08NYs67Zd3i7Dx+Yh3r7jjwHjujh9B8ccCI9s4c7iRGb9Xa6BXqxD++pNurptT5Oz3ZNr/mx3ru+XvMje88Pb6+vF76h1zzy/Q+v81SudeRnrvcvAnUdW8fP6kT+unDkz797m+60FCcVjmtuzDvctddeS8+ePXnrrbe47bbbTrWuo8TFxWG328nOzq7Vnp2dTVJSUp3nJCUlndDxAC6X66iHm7Rk83bPo6iyiNTwVM5NPtfqckREWiybYQuMURaR1qtR7yQ499xzmT9/fmNeMsDpdNK/f/9a1/f5fMyfP58hQ4bUec6QIUOOqmfevHn1Ht8avbflPQCu6XSNbhwREREROY6T6jmuS1lZGbNmzSI1NbWxLnmUadOmMXnyZAYMGMCgQYOYOXMmJSUl3HLLLQBMmjSJ1NRUZsyYAcBdd93FhRdeyFNPPcUVV1zBm2++yYoVK3jxxRdPW43NyY7CHazKWYXNsDG201iryxERERFp9k4qHMfExNQa02KaJkVFRYSGhvL66683WnFHmjBhAgcOHODBBx8kKyuLPn368PnnnwduusvIyMBmO9Q7OnToUGbPns3vfvc7fvvb39K5c2c++OCDM2aO4/e3vA/A+annkxiWeJyjRUREROSk5zk+nM1mIz4+nsGDBxMTE9NoxTUHLXWe40pvJSPeGcFBz0FmXTSLi9pdZHVJIiIiIpZpdfMcy4n5KvMrDnoOEh8Sz/ltz7e6HBEREZEW4aTHHB88eJCXXnqJjRs3AtCjRw9uueUWYmM1j25z8P5W/5CKsZ3G4rA12tByERERkVbtpKYvWLRoEenp6cyaNYuDBw9y8OBBZs2aRYcOHVi0aFFj1ygnaG/xXr7b539q4DWdr7G4GhEREZGW46S6FKdMmcKECRN47rnnsNvtgP+peL/61a+YMmUK69ata9Qi5cR8uO1DTEzOTT6XtIiW+wATERERkaZ2Uj3H27Zt49577w0EYwC73c60adPYtm1boxUnJ+eLXV8AcNVZV1lciYiIiEjLclLhuF+/foGxxofbuHEjvXv3PuWi5ORtO7iNHYU7CLIFMTxtuNXliIiIiLQoJzWs4s477+Suu+5i27ZtnHuu/5HE33//Pc8++yxPPPEEa9euDRx7zjnnNE6l0iBf7Pb3Gg9NGUqEM8LiakRERERalpOa5/jwB23UeVHDwDRNDMPA6/WedHHNQUub5/jqD65mR+EO/njeHxlz1hiryxERERFpFk7rPMc7d+486cLk9NGQChEREZFTc1LhuH379o1dhzSCL3d/CWhIhYiIiMjJOumnQ2zdupUFCxaQk5ODz+erte/BBx885cLkxH25yx+OL02/1OJKRERERFqmkwrH//znP/l//+//ERcXR1JSEoZhBPYZhqFwbIFtB7exvXC7hlSIiIiInIKTCsePPfYYjz/+OPfff39j1yMn6fAhFZHO5n/joIiIiEhzdFLzHB88eJDrrruusWuRU6AhFSIiIiKn7qTC8XXXXceXX37Z2LXISdpesJ3thdtx2BwaUiEiIiJyCk5qWEWnTp34/e9/z/fff0+vXr0ICgqqtf/OO+9slOKkYWp6jTWkQkREROTUnNRDQDp06FD/BQ2DHTt2nFJRzUlLeAjI2A/Gsr1wO4+f9zhXnXWV1eWIiIiINDt6CMgZQkMqRERERBpPg8PxtGnT+MMf/kBYWBjTpk2r9zjDMHjqqacapTg5Pg2pEBEREWk8DQ7HP/zwA5WVlYHt+hw+57GcfjVTuF3aXrNUiIiIiJyqBofjBQsW1Lkt1tlRsINtBdtw2Bxc1O4iq8sRERERafFOaio3aR4W7VkEwOCkwRpSISIiItIIWkw4zs/PZ+LEiURGRhIdHc1tt91GcXHxMY+/44476Nq1KyEhIbRr144777yTwsLCJqz69Fq8bzEA56WeZ3ElIiIiIq1DiwnHEydOZP369cybN4+5c+eyaNEifv7zn9d7/L59+9i3bx9/+ctf+PHHH3nllVf4/PPPue2225qw6tOntLKUVdmrAIVjERERkcZyUvMcN7WNGzfSo0cPli9fzoABAwD4/PPPufzyy9mzZw8pKSkNus4777zDT3/6U0pKSnA4GjbcurnOc7wgYwF3LriT1PBUPhv3mW6EFBERETmGhma6FtFz/N133xEdHR0IxgAjR47EZrOxdOnSBl+n5odxrGDs8Xhwu921luZoyb4lgL/XWMFYREREpHG0iHCclZVFQkJCrTaHw0FsbCxZWVkNukZubi5/+MMfjjkUA2DGjBlERUUFlrS0tJOu+3QxTZPFe/3jjc9PPd/iakRERERaD0vD8fTp0zEM45jLpk2bTvl93G43V1xxBT169ODhhx8+5rEPPPAAhYWFgSUzM/OU37+x7XTvZG/xXoJsQQxMGmh1OSIiIiKtxkk9Prqx3Hvvvdx8883HPKZjx44kJSWRk5NTq72qqor8/HySkpKOeX5RURGjR48mIiKCOXPmEBQUdMzjXS4XLperQfVbZfEef69x/8T+hAaFWlyNiIiISOthaTiOj48nPj7+uMcNGTKEgoICVq5cSf/+/QH46quv8Pl8DB48uN7z3G43o0aNwuVy8dFHHxEcHNxotVvp8PHGIiIiItJ4WsSY4+7duzN69Ghuv/12li1bxpIlS5g6dSo33HBDYKaKvXv30q1bN5YtWwb4g/Gll15KSUkJL730Em63m6ysLLKysvB6vVZ+nFNSVlXGiqwVgMYbi4iIiDQ2S3uOT8Qbb7zB1KlTGTFiBDabjfHjxzNr1qzA/srKSjZv3kxpaSkAq1atCsxk0alTp1rX2rlzJ+np6U1We2NanrWcCl8FyWHJdIjqYHU5IiIiIq1KiwnHsbGxzJ49u9796enpHD5l8/Dhw2kBUzifsJpZKjSFm4iIiEjjaxHDKuSQw8OxiIiIiDQuheMWZLd7N5lFmThsDgYn138jooiIiIicHIXjFqSm17hfQj/CgsIsrkZERESk9VE4bkE0pEJERETk9FI4biHKq8pZnrUcgGGpwyyuRkRERKR1UjhuIVZmr8Tj9ZAQmkDn6M5WlyMiIiLSKikctxA1QyrOTz1fU7iJiIiInCYKxy1ETTjWkAoRERGR00fhuAXILslml3sXNsOmKdxERERETiOF4xZgZfZKALrGdCXSGWlxNSIiIiKtl8JxC1ATjvsn9re4EhEREZHWTeG4BViRvQKAAUkDLK5EREREpHVTOG7m8sry2FG4A/A/GU9ERERETh+F42ZuVc4qADpFdyImOMbiakRERERaN4XjZk7jjUVERESajsJxM7ciS+ONRURERJqKwnEzVugpZMvBLQAMSFQ4FhERETndFI6bsR9yfsDEJD0ynbiQOKvLEREREWn1FI6bMY03FhEREWlaCsfNWM14Y4VjERERkaahcNxMlVSWsDF/IwADkwZaXI2IiIjImUHhuJlanbMar+klNTyVpLAkq8sREREROSMoHDdTGm8sIiIi0vRaTDjOz89n4sSJREZGEh0dzW233UZxcXGDzjVNk8suuwzDMPjggw9Ob6GNZEV29fzGmsJNREREpMm0mHA8ceJE1q9fz7x585g7dy6LFi3i5z//eYPOnTlzJoZhnOYKG095VTnrctcBCsciIiIiTclhdQENsXHjRj7//HOWL1/OgAH+sPj0009z+eWX85e//IWUlJR6z129ejVPPfUUK1asIDk5+bjv5fF48Hg8gddut/vUP8AJWntgLVW+KhJCEmgb0bbJ319ERETkTNUieo6/++47oqOjA8EYYOTIkdhsNpYuXVrveaWlpdx44408++yzJCU17Ka2GTNmEBUVFVjS0tJOuf4TVTOkon9S/xbV4y0iIiLS0rWIcJyVlUVCQkKtNofDQWxsLFlZWfWed8899zB06FCuvvrqBr/XAw88QGFhYWDJzMw86bpPVs3NeBpSISIiItK0LA3H06dPxzCMYy6bNm06qWt/9NFHfPXVV8ycOfOEznO5XERGRtZamlKFt4I1B9YACsciIiIiTc3SMcf33nsvN9988zGP6dixI0lJSeTk5NRqr6qqIj8/v97hEl999RXbt28nOjq6Vvv48eM5//zz+frrr0+h8tNnfd56PF4PscGxdIjqYHU5IiIiImcUS8NxfHw88fHxxz1uyJAhFBQUsHLlSvr398/7+9VXX+Hz+Rg8eHCd50yfPp2f/exntdp69erF3/72N8aMGXPqxZ8mhz8yWuONRURERJpWi5itonv37owePZrbb7+d559/nsrKSqZOncoNN9wQmKli7969jBgxgtdee41BgwaRlJRUZ69yu3bt6NCh+fbI6uEfIiIiItZpEeEY4I033mDq1KmMGDECm83G+PHjmTVrVmB/ZWUlmzdvprS01MIqT90N3W6gfWR7hqQMsboUERERkTOOYZqmaXURzZnb7SYqKorCwsImvzlPRERERBpHQzNdi5jKTURERESkKSgci4iIiIhUUzgWEREREanWYm7Is0rNkGy3221xJSIiIiJysmqy3PFut1M4Po6ioiIA0tLSLK5ERERERE5VUVERUVFR9e7XbBXH4fP52LdvHxEREU3yUA63201aWhqZmZmaHaOF03fZOuh7bB30PbYe+i5bByu+R9M0KSoqIiUlBZut/pHF6jk+DpvNRtu2bZv8fSMjI/UvfSuh77J10PfYOuh7bD30XbYOTf09HqvHuIZuyBMRERERqaZwLCIiIiJSTeG4mXG5XDz00EO4XC6rS5FTpO+yddD32Droe2w99F22Ds35e9QNeSIiIiIi1dRzLCIiIiJSTeFYRERERKSawrGIiIiISDWFYxERERGRagrHIiIiIiLV9IS842jqx0eLiIiISOPT46Mbyb59+0hLS7O6DBERERFpBJmZmbRt27be/QrHxxEREQH4f5B6hruIiIhIy+R2u0lLSwtku/ooHB9HzVCKyMhIhWMRERGRFu54w2R1Q56IiIiISDWFYxERERGRagrHzcyKXflc+9y3HCypsLoUERERkTOOwnEz4vWZ/HbOOlbsPsgD76/DNE2rSxIRERE5oygcNyN2m8Ffr++Dw2bw+fos3l25x+qSRERERM4oCsfNzNmpUUy7tAsAD3+0noy8UosrEhERETlzKBw3Q7+44CwGpcdSUuHlnrdXU+X1WV2SiIiIyBlB4bgZstsMnrq+NxEuByt3H+S5r7dbXZKIiIjIGUHhuJlKiw3l0bE9AZg5fyurMwusLUhERETkDKBw3IyN7ZPKleck4/WZ3PPWakorqqwuSURERKRVUzhuxgzD4PGxvUiOCmZnbgl/mLvR6pJEREREWjWF42YuKjSIp67rDcB/l2Xw+Y9ZFlckIiIi0nopHLcAQzvF8YsLOwLwm3fXsOegpncTEREROR0UjluI+y7tSp+0aNzlVdz53x+o1PRuIiIiIo1O4biFCLLbePonfYlwOViVUcDM/22xuiQRERGRVkfhuAVJiw3lifHnAPCPr7ezeGuuxRWJiIiItC4Kxy3MFeck85NB7TBNuOft1Rwo8lhdkoiIiEiroXDcAj14ZQ+6JIZzoMjDve+sweczrS5JREREpFVQOG6BQpx2nrmxH8FBNhZtOcCL3+ywuiQRERGRVkHhuIXqkhjBQ2P8j5f+8xebWbojz+KKRERERFo+heMW7IaBaYztk4LXZzL1vz+QU1RudUkiIiIiLVqLDsczZsxg4MCBREREkJCQwNixY9m8eXOtY8rLy5kyZQpt2rQhPDyc8ePHk52dbVHFjcswDP44rldg/PHU2T9QpfmPRURERE5aiw7HCxcuZMqUKXz//ffMmzePyspKLr30UkpKSgLH3HPPPXz88ce88847LFy4kH379jFu3DgLq25coU4Hz/20P+EuB8t25vPnLzcf/yQRERERqZNhmmarmergwIEDJCQksHDhQi644AIKCwuJj49n9uzZXHvttQBs2rSJ7t27891333Huuece95put5uoqCgKCwuJjIw83R/hpH26bj+/emMVAC/c1J9RPZMsrkhERESk+WhopmvRPcdHKiwsBCA2NhaAlStXUllZyciRIwPHdOvWjXbt2vHdd9/VeQ2Px4Pb7a61tASX90rmtvM6AHDf22vYlVtynDNERERE5EitJhz7fD7uvvtuhg0bxtlnnw1AVlYWTqeT6OjoWscmJiaSlZVV53VmzJhBVFRUYElLSzvdpTea6Zd1Y0D7GIo8Vfzy9ZWUVXitLklERESkRWk14XjKlCn8+OOPvPnmm6d0nQceeIDCwsLAkpmZ2UgVnn5BdhvP3NiPuHAnm7KK+O2cdbSiUTMiIiIip12rCMdTp05l7ty5LFiwgLZt2wbak5KSqKiooKCgoNbx2dnZJCXVPSbX5XIRGRlZa2lJkqKCefon/bDbDOb8sJeXFu+0uiQRERGRFqNFh2PTNJk6dSpz5szhq6++okOHDrX29+/fn6CgIObPnx9o27x5MxkZGQwZMqSpy20yQ85qw++v6A7AHz/dyKItByyuSERERKRlaNHheMqUKbz++uvMnj2biIgIsrKyyMrKoqysDICoqChuu+02pk2bxoIFC1i5ciW33HILQ4YMadBMFS3Z5KHpXD+gLT4Tps5epRv0RERERBqgRU/lZhhGne0vv/wyN998M+B/CMi9997Lf//7XzweD6NGjeIf//hHvcMqjtRSpnKri6fKyw0vfs8PGQV0TghnzpRhhLscVpclIiIi0uQamuladDhuCi05HAPkuMsZ88xist0eLumRyAs/7Y/NVvcvFSIiIiKt1Rk5z7EcLSEymBduGoDTYWPehmxmzt9qdUkiIiIizZbC8RmgT1o0M67pBcCs+Vv5eM0+iysSERERaZ4Ujs8Q4/u35fbz/bN53PvOGlbuzre4IhEREZHmR+H4DDL9su5c0iORiioft7+2koy8UqtLEhEREWlWFI7PIHabwd9v6EOv1CjySyq4+ZVlFJZWWl2WiIiISLOhcHyGCXU6eGnyAFKigtlxoIRfvL6Ciiqf1WWJiIiINAsKx2eghMhgXrp5IOEuB9/vyOe3c9ahGf1EREREFI7PWN2TI3nmxr7YbQbvrtzDM19ts7okEREREcspHJ/BhndN4OGregLw1LwtvL080+KKRERERKylcHyGu+nc9vzywrMAmP7+WuZtyLa4IhERERHrKBwL94/uyrX92+IzYersVSzbqTmQRURE5MykcCwYhsET43oxolsCnioft726nE1ZbqvLEhEREWlyCscCgMNu45kb+zGgfQxF5VVMemkZmfl6SIiIiIicWSwLx9u3b+d3v/sdP/nJT8jJyQHgs88+Y/369VaVdMYLcdp5afJAuiSGk1PkYfK/l5FX7LG6LBEREZEmY0k4XrhwIb169WLp0qW8//77FBcXA7BmzRoeeughK0qSalGhQbx262BSo0PYkVvCpH/rKXoiIiJy5rAkHE+fPp3HHnuMefPm4XQ6A+0XX3wx33//vRUlyWGSooJ57bZBtAlzsn6fm8kvL6OoXAFZREREWj9LwvG6deu45pprjmpPSEggNzfXgorkSGfFh/P6zwYTHRrE6swCbn1lOaUVVVaXJSIiInJaWRKOo6Oj2b9//1HtP/zwA6mpqRZUJHXpnhzJf24dTESwg+W7DvKzV1dQXum1uiwRERGR08aScHzDDTdw//33k5WVhWEY+Hw+lixZwn333cekSZOsKEnq0attFK/eOogwp51vt+fxy9dX4qlSQBYREZHWyZJw/Mc//pFu3bqRlpZGcXExPXr04IILLmDo0KH87ne/s6IkOYZ+7WL4980DCQ6y8fXmA0yd/QOVXp/VZYmIiIg0OsM0TdOqN8/MzGTdunUUFxfTt29fOnfubFUp9XK73URFRVFYWEhkZKTV5Vhq8dZcbn11ORVVPi7pkcgzN/bF5bBbXZaIiIjIcTU001nSc/zoo49SWlpKWloal19+Oddffz2dO3emrKyMRx991IqSpAHO6xzHCzf1x+mwMW9DNr/4z0qNQRYREZFWxZKeY7vdzv79+0lISKjVnpeXR0JCAl5v8wlc6jk+2uKtufzsteWUV/oY1qkN/5w0gFCnw+qyREREROrVrHuOTdPEMIyj2tesWUNsbKwFFcmJOK9zHK/e4r9Jb8m2PG7+93KKPZrmTURERFq+Jg3HMTExxMbGYhgGXbp0ITY2NrBERUVxySWXcP311zdlSXKSBndsw2u3DSbC5WDZrnxuemkphWV6UIiIiIi0bE06rOLVV1/FNE1uvfVWZs6cSVRUVGCf0+kkPT2dIUOGNFU5DaJhFce2bk8hP60OxmenRvLKLYOIC3dZXZaIiIhILQ3NdJaMOV64cCFDhw4lKCioqd/6hCkcH9/G/W5++q+l5JVUkN4mlNduHUy7NqFWlyUiIiIS0KzHHF944YWBYFxeXo7b7a61SMvSPTmSd345hNToEHbllTLuuW9Zv6/Q6rJERERETpgl4bi0tJSpU6eSkJBAWFgYMTExtRZpeTrGh/P+r4bSLSmC3GIPE174nm+351pdloiIiMgJsSQc//rXv+arr77iueeew+Vy8a9//YtHHnmElJQUXnvtNStKkkaQGBnMW78YwqAOsRR7qrj538v5ZO1+q8sSERERaTBLxhy3a9eO1157jeHDhxMZGcmqVavo1KkT//nPf/jvf//Lp59+2tQl1Utjjk9ceaWXu99czefrszAMeHhMTyYPTbe6LBERETmDNesxx/n5+XTs2BGAyMhI8vPzATjvvPNYtGiRFSVJIwoOsvPsxH5MHNwO04SHPlrPwx+tp8rrs7o0ERERkWOyJBx37NiRnTt3AtCtWzfefvttAD7++GOio6OtKEkamd1m8NjYs/n1qK4AvPLtLm57dQXucs2FLCIiIs2XJeH4lltuYc2aNQBMnz6dZ599luDgYO655x5+/etfW1GSnAaGYTDlok48/9N+BAfZWLjlANc+9y2Z+aVWlyYiIiJSJ0vC8T333MOdd94JwMiRI9m0aROzZ8/mhx9+4K677mrwdRYtWsSYMWNISUnBMAw++OCDWvtN0+TBBx8kOTmZkJAQRo4cydatWxvzo0gDjD47mXd+MZTESBdbsosZ++wSVu7Ot7osERERkaNYEo5fe+01PB5P4HX79u0ZN24c3bp1O6HZKkpKSujduzfPPvtsnfv/9Kc/MWvWLJ5//nmWLl1KWFgYo0aNory8/JQ/g5yYXm2j+HDKeZydGkleSQU/eXEp76zItLosERERkVosma3Cbrezf/9+EhISarXn5eWRkJCA1+s94WsahsGcOXMYO3Ys4O81TklJ4d577+W+++4DoLCwkMTERF555RVuuOGGBl1Xs1U0rtKKKu55azVfrM8GYOLgdjw4pgcuh93iykRERKQ1a9azVZimiWEYR7Xv2bOHqKioRnmPnTt3kpWVxciRIwNtUVFRDB48mO+++67e8zwej57YdxqFOh08N7E/94zsgmHAG0szmPDC9+wvLLO6NBEREREcTflmffv2xTAMDMNgxIgROByH3t7r9bJz505Gjx7dKO+VlZUFQGJiYq32xMTEwL66zJgxg0ceeaRRapC62WwGd43szDlpUdz95mpWZxZw5azFPP2TvgztFGd1eSIiInIGa9JwXDPkYfXq1YwaNYrw8PDAPqfTSXp6OuPHj2/Kko7ywAMPMG3atMBrt9tNWlqahRW1Xhd1TeDjqefxy9dXsmG/m5++tJTfjO7GLy7oWOdfFkREREROtyYNxw899BAA6enp3HDDDbhcrtP2XklJSQBkZ2eTnJwcaM/OzqZPnz71nudyuU5rXVJbuzahvP+rofzfnB95b9UenvhsE0t35PHn63oTF67vQURERJqWJWOOL774Yg4cOBB4vWzZMu6++25efPHFRnuPDh06kJSUxPz58wNtbrebpUuXMmTIkEZ7Hzl1wUF2/nLdOTx+zdk4HTYWbD7AZX//hsVbc60uTURERM4wloTjG2+8kQULFgAEbppbtmwZ//d//8ejjz7a4OsUFxezevVqVq9eDfhvwlu9ejUZGRkYhsHdd9/NY489xkcffcS6deuYNGkSKSkpgeEd0nwYhsHEwe35aOowOieEc6DIw03/XsoTn22iUo+dFhERkSZiSTj+8ccfGTRoEABvv/02vXr14ttvv+WNN97glVdeafB1VqxYQd++fenbty8A06ZNo2/fvjz44IMA/OY3v+GOO+7g5z//OQMHDqS4uJjPP/+c4ODgRv9M0ji6JUXy0dTzuHFwO0wTnl+4nWuf+5bdeSVWlyYiIiJnAEvmOQ4PD+fHH38kPT2dq666imHDhnH//feTkZFB165dKStrPtN6aZ5j63z+435+8+5a3OVVhDntPHB5dyYObqeb9UREROSENet5jnv27Mnzzz/PN998w7x58wLTt+3bt482bdpYUZI0Q6PPTuazuy9gUIdYSiq8/O6DH7nppWXsOVhqdWkiIiLSSlkSjp988kleeOEFhg8fzk9+8hN69+4NwEcffRQYbiECkBodwpu3n8uDV/YgOMjG4m25jJ75DW8uy8CCP3qIiIhIK2fJsArwP/TD7XYTExMTaNu1axehoaFHPVbaShpW0XzsOFDMr99dy8rdBwG4oEs8T4zrRUp0iMWViYiISHPX0ExnWTgGOHDgAJs3bwaga9euxMfHW1VKvRSOmxevz+TlJTv50xebqajyEea0c88lXbh5aDoOuyV/CBEREZEWoFmPOS4pKeHWW28lOTmZCy64gAsuuICUlBRuu+02Sks1nlTqZ7cZ/Oz8jnx65/kMaB9DSYWXxz7ZyJhnlvBDxkGryxMREZEWzpJwPG3aNBYuXMjHH39MQUEBBQUFfPjhhyxcuJB7773XipKkhemUEM7bvxjCE+N6ERUSxMb9bsY99y2/+2AdhWWVVpcnIiIiLZQlwyri4uJ49913GT58eK32BQsWcP3119d6ep7VNKyi+csr9vDHTzfx3qo9AMSFu7h/dFfG92uLzaZp30RERKSZD6soLS0lMTHxqPaEhAQNq5AT1ibcxVPX9+a/t59Lx/gwcos9/PrdtVz17GKW7sizujwRERFpQSzpOR4xYgRt2rThtddeCzytrqysjMmTJ5Ofn8///ve/pi6pXuo5blk8VV5e+3Y3s+ZvpchTBcDlvZJ44LLupMWGWlydiIiIWKVZz1axbt06Ro8ejcfjCcxxvGbNGlwuF19++SU9e/Zs6pLqpXDcMuUWe/jbvC38d1kGPhOcdhu3nJfOry7sRFRokNXliYiISBNr1uEY/EMr3njjDTZt2gRA9+7dmThxIiEhzWvOWoXjlm1TlpvH5m5k8bZcACKCHfzigo7cMqwDYS6HxdWJiIhIU2nW4XjGjBkkJiZy66231mr/97//zYEDB7j//vubuqR6KRy3fKZp8tWmHP78xWY2ZRUB0CbMya8u6sTEwe0IDrJbXKGIiIicbs36hrwXXniBbt26HdXes2dPnn/+eQsqktbMMAxGdE/k0zvP5+839CG9TSh5JRX8Ye4GLvrL17yxdDeeKq/VZYqIiEgzYEnPcXBwMBs3bqRDhw612nfs2EGPHj0oLy9v6pLqpZ7j1qfS6+O9lXv4+/yt7C/0/7OWGOni9vM78pNB7TTcQkREpBVq1j3HaWlpLFmy5Kj2JUuWkJKSYkFFciYJstu4YVA7Ftw3nIfG9CA5Kphst4fHPtnIeU9+xaz5Wyks1YNEREREzkSWdJHdfvvt3H333VRWVnLxxRcDMH/+fH7zm9/oCXnZG2D+ozD+n+CKsLqaVi04yM4twzowcXB75vywh+e+3s6uvFL+Om8LLy7awYSBadw8NF1TwImIiJxBLBlWYZom06dPZ9asWVRUVAD+oRb3338/Dz74YFOXc0xNOqzC54VnB0HeNuh6OUx4A2yWdO6fkbw+k0/W7ecfC7YFbtyzGTCyeyK3nteBwR1iMQw9cU9ERKQlatazVdQoLi5m48aNhISE0LlzZ1wul1Wl1KvJxxzvWQEvXw5eD5w3DUY+dPrfU2oxTZOvtxzg5SW7WLTl0KPMeyRHcsuwdMb0TtEMFyIiIi1MiwjHLYElN+SteQvm/Ny/Pe5fcM51TfO+cpRtOUW8vGQX763aQ3mlD4DIYAfX9E3lJ4Pb0S1JN2mKiIi0BArHjcSy2SrmPQRLZoLdBbd8Bm37N917y1EKSiv477JM3li6mz0HywLtfdKiuXFQO67snUyoU7NciIiINFcKx43EsnDs88KbN8KWzyE8CX6+ACI1k4fVfD6TxdtyeXN5Bl+uz6bK5//XJ8xpZ9TZSVzTN5WhZ8Vht2lssoiISHOicNxILJ3nuNwNL10KBzZCSl9/D3JQ83q89pnsQJGHd1fu4c3lGezOKw20J0S4uLpPCmP7ptIjOVI38YmIiDQDCseNxPKHgOTvhH9eDGX5cPZ4/xhkzWDRrJimyaqMg7y/ai9z1+6nsOzQHMkd48O4/OxkLuuVpKAsIiJiIYXjRmJ5OAbY+Q38Zyz4qmDg7XD5n0Ehq1mqqPLx9eYc5vywl/kbc6jw+gL72rcJZfTZSVx2djLnpEZh09ALERGRJqNw3EiaRTiG6hksfgGYMPQOuOQPCsjNnLu8kq825vDZj/v5evMBPFWHgnJ8hIuLuyZwcfcEzusUp0dWi4iInGYKx42k2YRjgBUvw9y7/dsXToeLHrC0HGm4Ek8VCzbn8NmPWSzYlENphTewz2m3MbhjLBd1TeD8znF0SgjX8AsREZFGpnDcSJpVOAb4/jn4fLp/e+TDcN49lpYjJ85T5WXpjny+2pTDV5tyyMgvrbU/IcLFeZ3iOK9zHMM6xZEYGWxRpSIiIq2HwnEjaXbhGOCbv8L8R/zbl/0JBv/C2nrkpJmmyfYDJXy1KZtvtuaybGd+reEX4L+pb1B6LAPTYxnUIZa2MSHqWRYRETlBCseNpFmGY4CvHodFf/JvX/k3GHCrtfVIoyiv9LJy90EWb8tlybZc1u0t5Mh/Q5MigxnYIZZ+7aLpkxZNj5RIXA49zlpERORYFI4bSbMNx6YJX/4OvnvG//qi/4MLfq2b9FqZgtIKVuw6yPJd+Szblc+6PYWBB4/UcNptdE+JpE/bKPq0i6ZnShQd48Jw2DXln4iISA2F40bSbMMx+APy/Edh8V/9r/tNhiv+CnbNfNBalVV4+SHzICt2HWR1ZgGrMwvIL6k46jiXw0a35Eh6JEfSMyWSHimRdE4IJyI4yIKqRURErKdw3EiadTiuseyf8NlvwPRBl9Fw7b/BGWZ1VdIETNMkM7+M1XsKWJ1RwLq9BWzY56bksNkwDpcaHULnxHC6JkbQOTGCTgnhdIgLIypEoVlERFo3heNG0iLCMcDGufDebVBVDqn94ca3ISzO6qrEAj6fye78UtbvK2T9Pjfr97nZnOUm2+2p95y4cCcd4/xBuUN8GO1jQ2nXJpS02FAi1dssIiKtgMJxI2kx4RggYyn8dwKUHYTYjnDjOxDXyeqqpJkoLK1kS04Rm7OK2JpdxObsInYcKCGnqP7QDBATGkS72FDaxoaSGh1CSlQwKdEhpESHkBodQnRokGbPEBGRZk/h+DDPPvssf/7zn8nKyqJ37948/fTTDBo0qEHntqhwDJC7FV4fBwUZ4AyHq2bB2eOtrkqasaLySnbllrIjt5gdB0rYmVtC5sFSMvNLyS0+ejzzkVwOG4mRwSRGuqrX/u34CBdtwlzEhbuIi3DSJsyFXY/MFhERiygcV3vrrbeYNGkSzz//PIMHD2bmzJm88847bN68mYSEhOOe3+LCMUBRNrx7K+xe7H894FYYNQOC9DAJOTElnioyD5aSkVdKRn4p+wvL2VdQxr6CMvYWlJNbfOxe58MZBsSGOokJc1avg4gNcxId6iQmNIioEP8SGRJEdIiTqNAgIoIdhDsd2BSqRUTkFCkcVxs8eDADBw7kmWf8U575fD7S0tK44447mD59+nHPb5HhGMBbBV/PgG+eAkxI6gXXvQptzrK6MmlFyiu9HCjykO0uJ8tdTrbbQ071dm6xh9yiCnKLPeSXVhw1X/OJCHc5/EG5eh3mchDmdBDqshPuchDqdBDmtBNSvYQ67YQE2QlxOgh22AgOshMcZMcV2LbhdNhw2m2a8q458fnAVwWmF3ze6m3fYdvV7abXf6zpPbQ/sO3zr2te1yy+mtfmoTbM2scEFvPQcbWO8W+bplm9eP1rnwmYmKYP0+d/iI9ZfRyBY32B7UAbpv/ypg8TEyPQzmHHgGH6MPHfgOt/n+r9NT+3w69VvYbqbf/GofOorqG62TjsOmatY2v/Cxuo6bD3DLQfftzhb3lEW62a6jmv1vsd+Z619h/9C3Pta5hHbBlHXffwqU9r13jsX8aP+sx1NdYx1Mz/HnW0Yx51fP3H1tlc/d3XU/eR1677qGOq+9r11Ffvz6/u9h5j7yUyus1JVHXiGprpWvWcXxUVFaxcuZIHHngg0Gaz2Rg5ciTfffddned4PB48nkO9YW63+7TXeVrYHTDi99B+CLz/c8haBy9cWD3MYpzV1UkrERxkJy3Wf+PesXh9Jvkl/qB8sLSCgyWV5JdWcLCkgvySCgrLKiksq6SgtGa7CndZJRVef9Ao9lRR7Kk6LZ/BbjNw2qvDcnVgDrIbBFUHZ6fdwGG34bAZOOwGDlvtbZvNwGEzsBnVa5uB3QY2w99mtxnYDLBVH2Mz/PuM6m0DA8Pw/2fDMAiM367575nBka+PFgg3h23YfB7svgrsPg92XyV2nweb1//a5qvE5qvA5qvA7qvA5qsMrG2+Smxmpf9Y87DXvipsZlV1m3+7VptZhd2swjC92M3q/YFtLza8/nXNghf74e34sJ3Uf7abnkH9//kXkROzv3Byk4XjhmrV4Tg3Nxev10tiYmKt9sTERDZt2lTnOTNmzOCRRx5pivKaRqeR8MvF8O5tkPEtvHsLbPwIRj8JEYnHP1+kEdhtBvER/nHIJ6K80ktRuT8YF5VXUlxehbu8itKKKko8VZRUeCn1VFHs8VJaUUVZpZfSCi/llV7KKqq3q7x4Kn2UV/rbPVW+Wg9S8fpMynxeyirrnv6ucZk4qSKUcsIoJ9TwEEY5IYaHEDyE4iHEqF7jIcSowEWFf5sKQgwPwVT4F6N6TSXBVOAyKnFVv3YZlU3wWZpWlWmjJkIfvvZiw8TAhxFoN82abQPfYWuTw9sNzOrtmnYTA9Os2UfgHLPWazADayNwfs1rDts+cj/V21TvN4zDz/H/GmQa/tf+bf85VJ97qAew5jelQ68PXavmWP9RNdc97H8C73Mo4BuHLnnUtQ+deujs2p2Rgfc54jcGo9b2od/ujvzFovZ5Rq1244hfmIyjDztUk2ke9jnqr6XOGgJb9bzfUXXW7D/6F7pjvncdDMyjDzJNjDrPNOt+1lc9f5qru7+3vl9C6+t7bvgvrfV/1vqv0SUsqsHXbyqtOhyfjAceeIBp06YFXrvdbtLS0iysqBFEpsDkj+HrP8LimbB+Dmz7Ci55xP/gEJv+rCzNU81wiBMN1cdT5fVR4fXhqfSvK6p8eKr86yqfj0qvj4oqk0qvf7uqwoOtogjDUwieIuwVbmwVRdgrS7BXFuOoLMZR5V+CKksI8pbg8JbirColyFuK01uC01eG01uGjaYI4Yf4MPAaQVQZTqpsTqoO2/YaQXht/te+w7a9Nkfgtc8IwmcLwlfd5rMF4TPsmNVtpuHAtAVh2hz4Dts2bQ4I7LeDLQhs/vOw2TENO9iDMGx2DNthxxg2bDYHpt2OYbODzY7NsPt72o/ofcfgiN74Qz3xNT32UHu/zTCwGxBU3bNf02tfc22Dw651xHk117dV9/DXrAPXOPya9RwrIs1fqw7HcXFx2O12srOza7VnZ2eTlJRU5zkulwuXq3H/Q9ws2B0w4kHoMRY+vhP2/QBz74Y1b8KYv0NCN6srFDm9Kkr90xyW5eMoO4ijrIDQ8oLqtgIoL4DywrqXqvLGr8cRDEGh/lllnKH+7aDQQ9vOMAgKqV5C/WtHiP/G2lrrkOprBfvXgcUFjmBs9iBshoFmqxYRaZhWHY6dTif9+/dn/vz5jB07FvDfkDd//nymTp1qbXFWST4HfjYflr0I8/8Amd/D8+fBsDth2F0Q3Pz+vCFylIpSKM2F0jwoyfOvD1/K8qG0einL9wfgxgi4zghwRUBwJLgiq9cR/oDrqt52hVe/rmmvfu0M9wfemsVmP/V6RESk0bXqcAwwbdo0Jk+ezIABAxg0aBAzZ86kpKSEW265xerSrGOzw7n/D7pdCZ/eB1s+989qsfwlOO9uGPQLf++VSFPx+fwBtjgbSnKg+IB/XXKgesmtvV1ZenLvY3NASMyhJTgaQqIPbQdH+V8HR9VeaoKvAq2ISKvX6qdyA3jmmWcCDwHp06cPs2bNYvDgwQ06t8VO5dZQpgmbP4X5j8KB6psUwxPhgl/7xyM7nNbWJy1blQeKsvyht2i/fw7u4iOXHP9inuBYXLsTQuMgtA2EtfGva5aQWAitXkJi/eE3NNbfe6txnyIiZyTNc9xIWn04ruHzwrp3YMEfoWC3vy2qHZx3F5xzg/9PwyI1fD5/L27RPnDv96+Lsqq3a5Ys/5CGExES4//lLCwewhMgLAHC4vyvw+L926Ft/NuuCAVdERFpMIXjRnLGhOMaVRXww2uw8M9QnOVvc0VB34kw8Gd6iMiZwOf19+gW7gV3zbLvsHV1+PU1cLowuxMikiA8yT99YHiSPwCHJ1S3J/hfh8bpLxUiInLaKBw3kjMuHNeoKIUf/gNLX4D87dWNBnS+FAb/HDpepPGXLZFp+m9Sc++BwsMW997q7b3+4NugIQ6GP9RGpkBEMkQm+9cRyf7QW7MOiVEPr4iIWE7huJGcseG4hs8H27+Cpc/DtnmH2sMToec10HMctB2ouZKbiyrPYUF3DxRkQmFm7SBcVXb86xh2f+iNTIHI1LrX4Yn+KQJFRERaAIXjRnLGh+PD5W2HZf+ENf/1zwlbIyrNH5R7jIWUPupRPl1M0z+jw+FhtyCjdvAtzqZBTzMKS4CoVIhqC5Ft/euo1OrtVH/w1fcoIiKtiMJxI1E4rkNVBez4Gn58DzZ9AhVFh/aFxECHC/zDLjoOh9gOVlXZ8lSW+cf0FmYeGu9bs10TiBsyhZkjpDrs1ixptV9HpvofGCEiInIGUThuJArHx1FZBlvn+YPy9q/A4669PyYd2g2F1H6Q0hcSzz4zg1lFyaFZHdz7/WN+3ftq3/RWmtewa4Ul+ENudNoRwbd6O7SNxviKiIgcQeG4kSgcnwBvFexbBdsXwI4FsGc5+KpqH2NzQEIP//CL+G7QppN/iW4H9hb2gFvT9A8vKcmtnsc3q3rO3qxD8/kW7feHYU9hw64ZFHZoiEOtIQ/q9RURETkVCseNROH4FHiKYNcS2LsS9v3gD8719Y7aHBDd3j8Mo2aWg/DEQ9uhbQ49rrexQ7TPCxXF/no9RVDu9veAlxUcevRwac06r/YT2xo6nRn4H0BRM6tDZFv/TW1RqdU3uFXf5KaZHURERE6LhmY63Woup48rArqO9i/g72kt3OMPyvvXQN42/01++dv9Y2nztx82bdwxOIIPPc7X4fKHZVuQfz5du8MftE0TTF/txVsJVeX+oSCBtQcqS07xc0b6H05x+Dy+NevDpzcL1i9XIiIizZ3CsTQdw/CPk41Ogx5XHWo3Tf/wg7xtcHBX9ZCELP8whZpHD5fmHwqxVeX+pSSnceuzBfkDrCvCv4TEVC81jyGu3g4/7KltoXEa5iAiItKKKByL9Qzj0Jy6HS6o/zhvlX9mjHL3oSEQXo+/R9hbUb2u9I9zNmzVi+Gfksyw+XuUHcEQFFJ77Qw/1AutIQ0iIiJnNIVjaTnsjkO9uSIiIiKngR5rJiIiIiJSTeFYRERERKSawrGIiIiISDWNOT6Ommmg3W73cY4UERERkeaqJssd7xEfCsfHUVRUBEBaWprFlYiIiIjIqSoqKiIqKqre/XpC3nH4fD727dtHREQERhNM8+V2u0lLSyMzM1NP5Gvh9F22DvoeWwd9j62HvsvWwYrv0TRNioqKSElJwWarf2Sxeo6Pw2az0bZt2yZ/38jISP1L30rou2wd9D22DvoeWw99l61DU3+Px+oxrqEb8kREREREqikci4iIiIhUUzhuZlwuFw899BAul8vqUuQU6btsHfQ9tg76HlsPfZetQ3P+HnVDnoiIiIhINfUci4iIiIhUUzgWEREREammcCwiIiIiUk3hWERERESkmsJxM/Pss8+Snp5OcHAwgwcPZtmyZVaXJCdgxowZDBw4kIiICBISEhg7diybN2+2uiw5RU888QSGYXD33XdbXYqchL179/LTn/6UNm3aEBISQq9evVixYoXVZckJ8Hq9/P73v6dDhw6EhIRw1lln8Yc//AHNKdD8LVq0iDFjxpCSkoJhGHzwwQe19pumyYMPPkhycjIhISGMHDmSrVu3WlNsNYXjZuStt95i2rRpPPTQQ6xatYrevXszatQocnJyrC5NGmjhwoVMmTKF77//nnnz5lFZWcmll15KSUmJ1aXJSVq+fDkvvPAC55xzjtWlyEk4ePAgw4YNIygoiM8++4wNGzbw1FNPERMTY3VpcgKefPJJnnvuOZ555hk2btzIk08+yZ/+9Ceefvppq0uT4ygpKaF37948++yzde7/05/+xKxZs3j++edZunQpYWFhjBo1ivLy8iau9BBN5daMDB48mIEDB/LMM88A4PP5SEtL44477mD69OkWVycn48CBAyQkJLBw4UIuuOACq8uRE1RcXEy/fv34xz/+wWOPPUafPn2YOXOm1WXJCZg+fTpLlizhm2++sboUOQVXXnkliYmJvPTSS4G28ePHExISwuuvv25hZXIiDMNgzpw5jB07FvD3GqekpHDvvfdy3333AVBYWEhiYiKvvPIKN9xwgyV1que4maioqGDlypWMHDky0Gaz2Rg5ciTfffedhZXJqSgsLAQgNjbW4krkZEyZMoUrrrii1r+X0rJ89NFHDBgwgOuuu46EhAT69u3LP//5T6vLkhM0dOhQ5s+fz5YtWwBYs2YNixcv5rLLLrO4MjkVO3fuJCsrq9b/x0ZFRTF48GBLs4/DsneWWnJzc/F6vSQmJtZqT0xMZNOmTRZVJafC5/Nx9913M2zYMM4++2yry5ET9Oabb7Jq1SqWL19udSlyCnbs2MFzzz3HtGnT+O1vf8vy5cu58847cTqdTJ482erypIGmT5+O2+2mW7du2O12vF4vjz/+OBMnTrS6NDkFWVlZAHVmn5p9VlA4FjlNpkyZwo8//sjixYutLkVOUGZmJnfddRfz5s0jODjY6nLkFPh8PgYMGMAf//hHAPr27cuPP/7I888/r3Dcgrz99tu88cYbzJ49m549e7J69WruvvtuUlJS9D1Ko9OwimYiLi4Ou91OdnZ2rfbs7GySkpIsqkpO1tSpU5k7dy4LFiygbdu2VpcjJ2jlypXk5OTQr18/HA4HDoeDhQsXMmvWLBwOB16v1+oSpYGSk5Pp0aNHrbbu3buTkZFhUUVyMn79618zffp0brjhBnr16sVNN93EPffcw4wZM6wuTU5BTb5pbtlH4biZcDqd9O/fn/nz5wfafD4f8+fPZ8iQIRZWJifCNE2mTp3KnDlz+Oqrr+jQoYPVJclJGDFiBOvWrWP16tWBZcCAAUycOJHVq1djt9utLlEaaNiwYUdNp7hlyxbat29vUUVyMkpLS7HZakcWu92Oz+ezqCJpDB06dCApKalW9nG73SxdutTS7KNhFc3ItGnTmDx5MgMGDGDQoEHMnDmTkpISbrnlFqtLkwaaMmUKs2fP5sMPPyQiIiIwZioqKoqQkBCLq5OGioiIOGqceFhYGG3atNH48RbmnnvuYejQofzxj3/k+uuvZ9myZbz44ou8+OKLVpcmJ2DMmDE8/vjjtGvXjp49e/LDDz/w17/+lVtvvdXq0uQ4iouL2bZtW+D1zp07Wb16NbGxsbRr1467776bxx57jM6dO9OhQwd+//vfk5KSEpjRwhKmNCtPP/202a5dO9PpdJqDBg0yv//+e6tLkhMA1Lm8/PLLVpcmp+jCCy8077rrLqvLkJPw8ccfm2effbbpcrnMbt26mS+++KLVJckJcrvd5l133WW2a9fODA4ONjt27Gj+3//9n+nxeKwuTY5jwYIFdf53cfLkyaZpmqbP5zN///vfm4mJiabL5TJHjBhhbt682dKaNc+xiIiIiEg1jTkWEREREammcCwiIiIiUk3hWERERESkmsKxiIiIiEg1hWMRERERkWoKxyIiIiIi1RSORURERESqKRyLiIiIiFRTOBYRacW+/vprDMOgoKDA6lJERFoEPSFPRKQVGT58OH369GHmzJkAVFRUkJ+fT2JiIoZhWFuciEgL4LC6ABEROX2cTidJSUlWlyEi0mJoWIWISCtx8803s3DhQv7+979jGAaGYfDKK6/UGlbxyiuvEB0dzdy5c+natSuhoaFce+21lJaW8uqrr5Kenk5MTAx33nknXq83cG2Px8N9991HamoqYWFhDB48mK+//tqaDyoichqp51hEpJX4+9//zpYtWzj77LN59NFHAVi/fv1Rx5WWljJr1izefPNNioqKGDduHNdccw3R0dF8+umn7Nixg/HjxzNs2DAmTJgAwNSpU9mwYQNvvvkmKSkpzJkzh9GjR7Nu3To6d+7cpJ9TROR0UjgWEWkloqKicDqdhIaGBoZSbNq06ajjKisree655zjrrLMAuPbaa/nPf/5DdnY24eHh9OjRg4suuogFCxYwYcIEMjIyePnll8nIyCAlJQWA++67j88//5yXX36ZP/7xj033IUVETjOFYxGRM0xoaGggGAMkJiaSnp5OeHh4rbacnBwA1q1bh9frpUuXLrWu4/F4aNOmTdMULSLSRBSORUTOMEFBQbVeG4ZRZ5vP5wOguLgYu93OypUrsdvttY47PFCLiLQGCsciIq2I0+msdSNdY+jbty9er5ecnBzOP//8Rr22iEhzo9kqRERakfT0dJYuXcquXbvIzc0N9P6eii5dujBx4kQmTZrE+++/z86dO1m2bBkzZszgk08+aYSqRUSaD4VjEZFW5L777sNut9OjRw/i4+PJyMholOu+/PLLTJo0iXvvvZeuXbsyduxYli9fTrt27f5/e3doBAAAwkBs/60xPwKiIhkA/VfDy32AFT7kAQBALMcAABBxDAAAEccAABBxDAAAEccAABBxDAAAEccAABBxDAAAEccAABBxDAAAEccAAJAD2yXCrBq+docAAAAASUVORK5CYII=", "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "time = jnp.arange(T + 1) * dims.dt\n", "fig, axes = plt.subplots(3, 1, figsize=(7, 5), sharex=True, layout=\"constrained\")\n", "axes[0].plot(time, Xs, label=[\"position\", \"velocity\"])\n", "axes[0].legend()\n", "axes[1].plot(time[:-1], Us, c=\"C2\")\n", "axes[2].plot(time, lambdas)\n", "for ax, label in zip(axes, [\"state\", \"input\", \"costate\"]):\n", " ax.set_ylabel(label)\n", "axes[-1].set_xlabel(\"time\");" ] }, { "cell_type": "markdown", "id": "1d7be5c7", "metadata": {}, "source": [ "## Checking optimality\n", "\n", "At the optimum, the gradients of the Lagrangian with respect to $x$, $u$ and $\\lambda$ are zero.\n", "`diffilqrax.lqr.kkt` returns these residuals. They should be at the level of float64\n", "round-off." ] }, { "cell_type": "code", "execution_count": 5, "id": "175eea0b", "metadata": { "execution": { "iopub.execute_input": "2026-10-14T22:41:15.483235Z", "iopub.status.busy": "2026-10-14T22:41:15.483016Z", "iopub.status.idle": "2026-10-14T22:41:15.888947Z", "shell.execute_reply": "2026-10-14T22:41:15.888130Z" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "max|dL/dx| = 0.00e+00\n", "max|dL/du| = 1.89e-15\n", "max|dL/dλ| = 2.78e-17\n" ] } ], "source": [ "dLdX, dLdU, dLdLambda = kkt(params, Xs, Us, lambdas)\n", "for name, r in [(\"dL/dx\", dLdX), (\"dL/du\", dLdU), (\"dL/dλ\", dLdLambda)]:\n", " print(f\"max|{name}| = {jnp.max(jnp.abs(r)):.2e}\")" ] }, { "cell_type": "markdown", "id": "a31549c1", "metadata": {}, "source": [ "The parallel solver uses an associative scan instead of the sequential recursion. It returns the\n", "same solution:" ] }, { "cell_type": "code", "execution_count": 6, "id": "6c26c325", "metadata": { "execution": { "iopub.execute_input": "2026-10-14T22:41:15.892530Z", "iopub.status.busy": "2026-10-14T22:41:15.892249Z", "iopub.status.idle": "2026-10-14T22:41:17.542967Z", "shell.execute_reply": "2026-10-14T22:41:17.542249Z" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "max|Xs - Xp| = 4.44e-16, max|Us - Up| = 1.28e-15\n" ] } ], "source": [ "Xp, Up, _ = solve_plqr(params)\n", "print(f\"max|Xs - Xp| = {jnp.max(jnp.abs(Xs - Xp)):.2e}, max|Us - Up| = {jnp.max(jnp.abs(Us - Up)):.2e}\")" ] }, { "cell_type": "markdown", "id": "31fd73cc", "metadata": {}, "source": [ "## How to adapt\n", "\n", "- **Your own linear system:** replace `A`, `B` and `dims`. If the system is time-varying, build\n", " arrays of shape `(T, n, n)` and `(T, n, m)` directly instead of tiling.\n", "- **Costs:** a larger `R` relative to `Q` gives gentler inputs and slower convergence.\n", " `Qf` sets how strongly the final state is pinned. `Q`, `Qf` must be positive semi-definite and\n", " `R` positive definite.\n", "- **Offsets and cross terms:** use `a` for affine dynamics (for example gravity) and `S` for\n", " state–input coupling in the cost.\n", "- **Tracking a target:** use the linear terms `q` and `qf`, as shown in the\n", " [tracking example](lqr_tracking_integrator.ipynb).\n", "- **Long horizons on a GPU:** switch to `solve_plqr`. It has the same inputs and outputs." ] } ], "metadata": { "kernelspec": { "display_name": "Python 3", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.11.11" } }, "nbformat": 4, "nbformat_minor": 5 }