From aeda321f97d392ecde05a5631215552823056b00 Mon Sep 17 00:00:00 2001 From: Lev Kozlov Date: Mon, 16 Sep 2024 23:22:07 +0900 Subject: [PATCH] feat: add mppi like sampling of parameters --- examples/nonlinear_residuals.ipynb | 439 ++++++++++++++++++++++++++--- 1 file changed, 396 insertions(+), 43 deletions(-) diff --git a/examples/nonlinear_residuals.ipynb b/examples/nonlinear_residuals.ipynb index da292db..3d64e93 100644 --- a/examples/nonlinear_residuals.ipynb +++ b/examples/nonlinear_residuals.ipynb @@ -14,7 +14,9 @@ "from robot_descriptions.skydio_x2_mj_description import MJCF_PATH\n", "\n", "from mujoco_sysid.mjx.convert import logchol2theta, theta2logchol\n", - "from mujoco_sysid.mjx.parameters import get_dynamic_parameters, set_dynamic_parameters\n" + "from mujoco_sysid.mjx.parameters import get_dynamic_parameters, set_dynamic_parameters\n", + "\n", + "import matplotlib.pyplot as plt" ] }, { @@ -65,22 +67,114 @@ "cell_type": "code", "execution_count": 3, "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "((1001, 7), (1001, 6), (1001, 4))" + ] + }, + "execution_count": 3, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "log_qpos = jnp.array(log.data(\"qpos\"))\n", + "log_qvel = jnp.array(log.data(\"qvel\"))\n", + "log_ctrl = jnp.array(log.data(\"ctrl\"))\n", + "\n", + "log_qpos.shape, log_qvel.shape, log_ctrl.shape" + ] + }, + { + "cell_type": "code", + "execution_count": 4, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "(Array([[ 0.0000000e+00, 0.0000000e+00, 0.0000000e+00, -3.1506685e-14,\n", + " -3.7471137e-12, -2.4804585e-15],\n", + " [-5.7232217e-04, 5.5040164e-07, -4.9597770e-04, 1.0275190e-05,\n", + " 1.0678568e-02, -1.0522764e-06],\n", + " [-1.4207972e-03, 1.3654253e-06, -1.2307838e-03, 2.5685131e-05,\n", + " 2.6694240e-02, -2.6480805e-06],\n", + " ...,\n", + " [ 3.9635710e-02, -2.3915130e-01, 3.0338511e-02, -1.9615057e-01,\n", + " -6.3093685e-02, 1.5923199e-01],\n", + " [ 4.6400912e-02, -2.3677342e-01, 3.1514063e-02, -1.9395384e-01,\n", + " -6.8520784e-02, 1.5893112e-01],\n", + " [ 5.3117447e-02, -2.3419258e-01, 3.2680660e-02, -1.9162498e-01,\n", + " -7.3916078e-02, 1.5861784e-01]], dtype=float32),\n", + " Array([[ 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,\n", + " 0.00000000e+00, 0.00000000e+00, 0.00000000e+00],\n", + " [-2.84185503e-02, 2.73143050e-05, -2.46452875e-02,\n", + " 5.13718231e-04, 5.33897042e-01, 5.28622550e-05],\n", + " [-5.60338274e-02, 5.37765154e-05, -4.85276319e-02,\n", + " 1.02718384e-03, 1.06760776e+00, 1.07214546e-04],\n", + " ...,\n", + " [ 7.09591210e-01, 2.17365265e-01, 1.43795148e-01,\n", + " 2.47991294e-01, -1.05928528e+00, -3.08974367e-02],\n", + " [ 7.05333948e-01, 2.37760365e-01, 1.42507792e-01,\n", + " 2.60988653e-01, -1.05597985e+00, -2.98040118e-02],\n", + " [ 7.00467050e-01, 2.58056104e-01, 1.41612679e-01,\n", + " 2.73828089e-01, -1.05192769e+00, -2.87376605e-02]], dtype=float32))" + ] + }, + "execution_count": 4, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "import jaxlie\n", + "\n", + "\n", + "@jax.jit\n", + "def diff_qpos(qpos1, qpos2):\n", + " # qpos = [x, y, z, qw, qx, qy, qz]\n", + " quat1 = qpos1[3:][jnp.array([3, 0, 1, 2])]\n", + " quat2 = qpos2[3:][jnp.array([3, 0, 1, 2])]\n", + " q1 = jaxlie.SO3.from_quaternion_xyzw(quat1)\n", + " q2 = jaxlie.SO3.from_quaternion_xyzw(quat2)\n", + "\n", + " return jnp.concatenate([qpos1[:3] - qpos2[:3], jaxlie.SO3.log(q1.inverse() @ q2)])\n", + "\n", + "\n", + "@jax.jit\n", + "def diff_qvel(qvel1, qvel2):\n", + " return qvel1 - qvel2\n", + "\n", + "\n", + "diff_qpos(log_qpos[0], log_qpos[0]), diff_qvel(log_qvel[0], log_qvel[0])\n", + "\n", + "vmap_diff_qpos = jax.vmap(diff_qpos, in_axes=(0, None))\n", + "vmap_diff_qvel = jax.vmap(diff_qvel, in_axes=(0, None))\n", + "\n", + "vmap_diff_qpos(log_qpos, log_qpos[0]), vmap_diff_qvel(log_qvel, log_qvel[0])" + ] + }, + { + "cell_type": "code", + "execution_count": 5, + "metadata": {}, "outputs": [], "source": [ "mjx_model = mjx.put_model(model)\n", "\n", "\n", - "def smart_step(acc, vel, pos, ctrl, parameters):\n", + "def smart_step(vel, pos, ctrl, parameters):\n", " # update the parameters\n", " new_model = set_dynamic_parameters(mjx_model, 1, parameters)\n", "\n", " mjx_data = mjx.make_data(new_model)\n", " # set initial data for the step\n", - " mjx_data = mjx_data.replace(qacc=acc, qvel=vel, qpos=pos, ctrl=ctrl)\n", + " mjx_data = mjx_data.replace(qvel=vel, qpos=pos, ctrl=ctrl)\n", " # step the simulation\n", " mjx_data = mjx.step(new_model, mjx_data)\n", " return mjx_data.qpos, mjx_data.qvel\n", - " # return mjx_data.qpos.at[0], mjx_data.qvel.at[0]\n", "\n", "\n", "smart_step = jax.jit(smart_step)" @@ -88,101 +182,360 @@ }, { "cell_type": "code", - "execution_count": 4, + "execution_count": 6, "metadata": {}, "outputs": [ { - "name": "stdout", - "output_type": "stream", - "text": [ - "[-2.8813668e-04 2.7725861e-07 9.9750474e-02 9.9999642e-01\n", - " 2.5690015e-06 2.6697947e-03 2.6182667e-07] [-2.8813669e-02 2.7725859e-05 -2.4952721e-02 5.1380089e-04\n", - " 5.3395963e-01 5.2365394e-05] [-2.8813670e+00 2.7725860e-03 -2.4952722e+00 5.1380090e-02\n", - " 5.3395962e+01 5.2365395e-03] [ 4.8408070e+00 4.8512077e+00 -5.7623768e-09 -5.7623768e-09]\n" - ] + "data": { + "text/plain": [ + "Array([ 1.325 , 0. , 0. , 0.0715 , 0.04051 , 0. ,\n", + " 0.02927 , -0.0021 , 0. , 0.060528], dtype=float32)" + ] + }, + "execution_count": 6, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "true_parameters" + ] + }, + { + "cell_type": "code", + "execution_count": 7, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "(1000, 10)" + ] + }, + "execution_count": 7, + "metadata": {}, + "output_type": "execute_result" } ], "source": [ - "pos = jnp.array(log.data(\"qpos\")[0])\n", - "vel = jnp.array(log.data(\"qvel\")[0])\n", - "acc = jnp.array(log.data(\"qacc\")[0])\n", - "ctrl = jnp.array(log.data(\"ctrl\")[0])\n", + "N_SAMPLES = 1_000\n", + "key = jax.random.PRNGKey(0)\n", + "\n", + "\n", + "vmapped_logchol2theta = jax.vmap(logchol2theta)\n", + "\n", + "\n", + "@jax.jit\n", + "def sample_parameters(theta):\n", + " \"\"\"\n", + " Get theta estimate (10,) and return (N_SAMPLES, 10)\n", + " \"\"\"\n", + " logchol = theta2logchol(theta)\n", + " distrib = jax.random.normal(key, shape=(N_SAMPLES, logchol.shape[0])) * 0.1\n", + " logchol = jnp.array([logchol] * N_SAMPLES)\n", + "\n", + " # cast to theta backward\n", + " new_logchol = logchol + distrib\n", + " return new_logchol\n", + " # return vmapped_logchol2theta(new_logchol)\n", + "\n", "\n", - "print(pos, vel, acc, ctrl)" + "sample_parameters(true_parameters).shape" ] }, { "cell_type": "code", - "execution_count": 5, + "execution_count": 8, + "metadata": {}, + "outputs": [], + "source": [ + "vmapped_step = jax.vmap(smart_step, in_axes=(0, 0, None, 0))\n", + "# vmapped_step(log_qvel[0], log_qpos[0], log_ctrl[0], sample_parameters(true_parameters))" + ] + }, + { + "cell_type": "code", + "execution_count": 9, + "metadata": {}, + "outputs": [], + "source": [ + "# def compute_error(logchols):\n", + "# thetas = vmapped_logchol2theta(logchols)\n", + "# qpos = log_qpos[0]\n", + "# qvel = log_qvel[0]\n", + "\n", + "# qpos = jnp.array([qpos] * N_SAMPLES)\n", + "# qvel = jnp.array([qvel] * N_SAMPLES)\n", + "\n", + "# errors = jnp.zeros(N_SAMPLES)\n", + "\n", + "# for i in range(len(log_qpos) - 1):\n", + "# ctrl = log_ctrl[i]\n", + "# qpos, qvel = vmapped_step(qvel, qpos, ctrl, thetas)\n", + "\n", + "# next_qpos = log_qpos[i + 1]\n", + "# next_qvel = log_qvel[i + 1]\n", + "\n", + "# errors += jnp.linalg.norm(vmap_diff_qpos(qpos, next_qpos), axis=1) ** 2\n", + "# errors += jnp.linalg.norm(vmap_diff_qvel(qvel, next_qvel), axis=1) ** 2\n", + "\n", + "# # replace nan with 1e10\n", + "# errors = jnp.where(jnp.isnan(errors), 1e10, errors)\n", + "\n", + "# return errors\n", + "\n", + "\n", + "# # compute_error = jax.jit(compute_error)\n" + ] + }, + { + "cell_type": "code", + "execution_count": 10, + "metadata": {}, + "outputs": [], + "source": [ + "import jax\n", + "import jax.numpy as jnp\n", + "\n", + "\n", + "def compute_error(logchols):\n", + " thetas = vmapped_logchol2theta(logchols)\n", + " qpos0 = log_qpos[0]\n", + " qvel0 = log_qvel[0]\n", + "\n", + " qpos = jnp.tile(qpos0[None, :], (N_SAMPLES, 1))\n", + " qvel = jnp.tile(qvel0[None, :], (N_SAMPLES, 1))\n", + "\n", + " errors = jnp.zeros(N_SAMPLES)\n", + "\n", + " def step_fn(carry, inputs):\n", + " qpos, qvel, errors = carry\n", + " ctrl, next_qpos, next_qvel = inputs\n", + "\n", + " # Update the positions and velocities\n", + " qpos, qvel = vmapped_step(qvel, qpos, ctrl, thetas)\n", + "\n", + " # # Compute the differences\n", + " # next_qpos = jnp.broadcast_to(next_qpos, qpos.shape)\n", + " # next_qvel = jnp.broadcast_to(next_qvel, qvel.shape)\n", + "\n", + " diff_qpos = vmap_diff_qpos(qpos, next_qpos)\n", + " diff_qvel = vmap_diff_qvel(qvel, next_qvel)\n", + "\n", + " # Accumulate the squared errors\n", + " errors += jnp.linalg.norm(diff_qpos, axis=1) ** 2\n", + " errors += jnp.linalg.norm(diff_qvel, axis=1) ** 2\n", + "\n", + " return (qpos, qvel, errors), None\n", + "\n", + " # Prepare inputs for scanning\n", + " inputs = (log_ctrl[:-1], log_qpos[1:], log_qvel[1:])\n", + "\n", + " # Run the scan\n", + " (qpos, qvel, errors), _ = jax.lax.scan(step_fn, (qpos, qvel, errors), inputs)\n", + "\n", + " # Replace NaNs with a large constant\n", + " errors = jnp.where(jnp.isnan(errors), 1e10, errors)\n", + "\n", + " return errors\n", + "\n", + "\n", + "# JIT-compile the compute_error function\n", + "compute_error = jax.jit(compute_error)\n" + ] + }, + { + "cell_type": "code", + "execution_count": 11, "metadata": {}, "outputs": [ { "data": { "text/plain": [ - "Array([ 1.325 , 0. , 0. , 0.0715 , 0.04051 , 0. ,\n", - " 0.02927 , -0.0021 , 0. , 0.060528], dtype=float32)" + "(0.0, 10000000.0)" ] }, - "execution_count": 5, + "execution_count": 11, "metadata": {}, "output_type": "execute_result" + }, + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAiMAAAGsCAYAAAAPJKchAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjkuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8hTgPZAAAACXBIWXMAAA9hAAAPYQGoP6dpAACooklEQVR4nO2dd5xVxfn/P+ferbAsnaUXKypVEETEEolo1Gg0/ogaQWOMMZiofJNYohKjBtL8mkI0atQ0S4oaEw1+FcWKohQVFTSKguhSRDpsu+f3x91795TpM6fc3XnnZdh7zply5syZec7zPPOM47quC4vFYrFYLJaEyCRdAYvFYrFYLB0bK4xYLBaLxWJJFCuMWCwWi8ViSRQrjFgsFovFYkkUK4xYLBaLxWJJFCuMWCwWi8ViSRQrjFgsFovFYkkUK4xYLBaLxWJJFCuMWCwWi8ViSRQrjFgsFovFYkmUkhJGnn32WZxyyino378/HMfBww8/LJX+hz/8IRzHCf3XuXPnaCpssVgsFouFS0kJI7t27cLo0aMxf/58pfTf/e538cknn/j+O/jgg3HmmWcarqnFYrFYLBZRSkoYOfHEE3HjjTfiS1/6EvF8Q0MDvvvd72LAgAHo3LkzJk6ciEWLFhXP19TUoG/fvsX/NmzYgLfeegsXXHBBTHdgsVgsFoslSEkJIzwuueQSLF68GPfffz9ef/11nHnmmTjhhBPw7rvvEq+/8847ccABB2DKlCkx19RisVgsFkuBdiOMrF27FnfffTf+9re/YcqUKdh3333x3e9+F0ceeSTuvvvu0PV79+7FX/7yF6sVsVgsFoslYcqSroAp3njjDbS0tOCAAw7wHW9oaEDPnj1D1z/00EPYsWMHZs6cGVcVLRaLxWKxEGg3wsjOnTuRzWaxdOlSZLNZ37mamprQ9XfeeSdOPvlk1NXVxVVFi8VisVgsBNqNMDJ27Fi0tLRg48aNXB+QNWvW4Omnn8YjjzwSU+0sFovFYrHQKClhZOfOnfjvf/9b/L1mzRqsWLECPXr0wAEHHIBzzjkHM2bMwC9+8QuMHTsWmzZtwsKFCzFq1CicdNJJxXR33XUX+vXrhxNPPDGJ27BYLBaLxeLBcV3XTboSoixatAjHHnts6PjMmTNxzz33oKmpCTfeeCP++Mc/Yv369ejVqxcOP/xwXH/99Rg5ciQAIJfLYciQIZgxYwZuuummuG/BYrFYLBZLAGlh5Nlnn8XPfvYzLF26FJ988gkeeughnHbaacw0ixYtwuzZs/Hmm29i0KBBuOaaa3DeeedpVNtisVgsFkt7QXppr2wU1DVr1uCkk07CscceixUrVuCyyy7D17/+dTz++OPSlbVYLBaLxdL+0DLTOI7D1YxcccUVePTRR7Fy5crisa985SvYunUrFixYoFq0xWKxWCyWdkLkDqyLFy/G1KlTfcemTZuGyy67jJqmoaEBDQ0Nxd+5XA5btmxBz5494ThOVFW1WCwWi8ViENd1sWPHDvTv3x+ZDN0YE7kwUl9fH4rlUVdXh+3bt2PPnj2orq4OpZk7dy6uv/76qKtmsVgsFoslBtatW4eBAwdSz6dyae9VV12F2bNnF39v27YNgwcPxrp161BbW5tgzcwyYk7eb+boA3ph/jnjir/7dKnEU989hptu3ODu+MMFE/CrJ9/F7c+9DwBY9N1j0KtLJf780geY95/VAICV108rpinwneP2wzeO2hdXP/gGHnntY9+5XjUVWPS98KolHrc/+x5+tbBt6fX8c8bit0+/hzc/3g4AOHfSEPxp8YehdKT6/e/0Mfj8weSAdH9c/AF+uiB/b89fcSy6daoonjty3kJs3dPsy9vLFX9/DY++UR8695MFq4p1e/Q7R2JIz878G6ZQv20Ppt78LADgn7OOwKnzX/SddxzgjR9OIyUtEmyPQn2n3fIM1n+2N1T/Z9/ZhG/9ZVm+zEsm483123D1Q3nT6OOXTcGA7p2YZQTbCQDOu2sJXv3wM+J5Utqbn3gHdz2/BgDw2pzjkc043DRezrr9Jbyxfhv1PIlgO/3u3HGYvF8vobSm8dZl7KBuWL5uK/G6H5w0HGdNGCKV5xmHDsA/lq0vHue1z86GZhz+44UAgB+ecjC+PH5QqK1I7x0p/4eWfYRr//kmAOC3Xz0U44d0x4Sb8nlfe8pBmD5+sNC9AMDpv30B72zYSbyHQl36d6vCph0NaGpxfdd569q7pgJPE8ao8Tc+gb1NueLv4PhAg9YOBf717ckY1iscPDPIyvVb8ZXbX/YdW3Hd51GWDWsEeO8Di7++ug4/+tdboeNPf/do9O5ShSferMflf32Nmf+rH2zBeXe/wqyjCbZv345BgwahS5cuzOsiF0YKu+N62bBhA2pra4laEQCorKxEZWVl6HhtbW27EkYylfkJoqK6BrW1tcXfmxuBbz7wFu467zBUlWep6cqrO6O2thZVnWuKx7rU1qK2SyU6de5SPObNu0B15y6ora1FZaea0LlsVaVSO1d7ygSATp27oLy6MzKVeeGgsjpcFq1+nWtqqHXwltOlSy1qO7cNNtmqzsjkmnx5e/Her/dcled4ly61qK1VF0Z2ueXFvGq6hO/NccL1CkJrp7KqzshUZkL171yz11P/LqjunPP3idpwft4ySPXJP7sG4nlSWn8bdgkNbmLlNVHPkwj3my6JjRHeuuTvpZF4XeHdk8kz+J7y0jt7m4rXd2ptk2BbkY6R8q+uaXvfOnf259VJ4l4AtPbfHPEeCnmWVVUj05hFplUYKVznrSttjMpWdkYm0+K7j1oBYYTWDgVqutSitpYvjNRszxHbmTTRyzzPIPW7HWKd82NXFaprdnHzr+nS7BsLoxJGCvBcLCLfKG/SpElYuHCh79gTTzyBSZMmRV10yUB6Ri++9yn+sewjofQu1EPFmHTB4flC69RTuA6K571V160lL68oIvsE29b7K65IQnE8X34d0o/KK5dkNKjge106kaniryupvFJorjTUUVoY2blzJ1asWIEVK1YAaIuCunbtWgB5E8uMGTOK13/zm9/E+++/j+9///tYtWoVfvvb3+Kvf/0rLr/8cjN30I7xqhtFKUwIpe7o632p71+yFi/8d7PnnPlXxzuRllAcQCKum9A9GBToLH4eeHWd1PUm298vXKf7yUY17JXKkKBazTTcn7Qw8uqrr2Ls2LEYO3YsAGD27NkYO3YsrrvuOgDAJ598UhRMAGDYsGF49NFH8cQTT2D06NH4xS9+gTvvvBPTpsnZyNo30QoOcU1MwWJ4v9uO0+u3cv02XPngGzjnzpeJ54MpeS0ZR1PIaiUam3P40b/ewrPvbIqsTnGQhDYmVIc0jKo8SuxDIZdAkzoRj4lREZcmlEcpvAZBpH1GjjnmGOYLf8899xDTLF++XLaokmZPYwvWfbYbB9SxnXaiwjvepaVjykwUhSs37tgrV4bU1Z50EX3V874kP966B39c/CHuemEN7nphDT6YdxLzemo5gfon8chLQhCwSOPXGpaW1stcXUvprsXwzREpuL9UrqZpD5w2/wWs3rADd503Hp8bTl4VUkD3Q4k3B8TVzaJQEVZkww68aUd0Ut6xtwlHzHsqokpEky2zyBSo85MfUvOw3uk4vvlNyoVJaEZUSVqfQnr343wXCsWnQbiQJXIH1o7K6g07AAAPLf+YcyUdrsmB0eG8aVmTo8mXN2SWCZ5XyLOirK2LNreEfWgi8R8xmCUrr3Vb9kjnJ6q+TmIwSoOZxhIBrkv6M1o43Zwm7EXlK5e2/mz6NtNwf1YYSQFG+5UBb26T9fELReRrWHZWrzDSSBBGZKGvpnEFrhIsQzC5yoBCEzJ8Wgm3A9up01AHDrG4jJgUqIO/Y3jQppoobtNhbKvnqGNpCbwAFKwwkmJEBy1i9/MkTsvSTvpESk9Xnm27jwbC6iJTtxZVE6Vhgk7L87eUJjmPnaajPuFSu+80jDuyWGEkBUT9pZSY/V4jPkGhzt4onkY0IwKViOtFVnnuImYaN6EnHtTOBIlDI6By5299vB2XP7AC67bsjqBGYUptpUgpzWtJt2zSY4dq+WkQXqwDaztFxDzyybY9+NtSscBqIvAC/lBNJIJ5FjQjcao9tfONaSgPq9JjKbZd8IVfPQcAWF2/A49dOiXh2pjBZL8LOrAaXW3mupHGRDKmOdXIqJTGqySxmpEUQPtSUorUWEgrkPjEXz6nUIJ8XYq/ZTQjrv9fAGhsaaFep1on1WuY6TkagqjJL7/0qtZjEohcdplJrSLZvrcJp81/AXe27t9E478bdxqrR6lpP1j4nmuggXXuc8uuRkye9xTm/WdV6FyJhWLxkA4pQbYWaTCxWmGkHcBd2ks5v3V3E/mEaj0IxxzuFez6e1+SQkTaKF6cYCwFY/ky8opywkr6yynp8r3c+dwarFi3FTc++nbSVYmNuPqwzrv4++ffx8fb9uK2Z95TziNECoWYKMYrqplGMb80vK9WGIkYkXcjimVq3okuNqmXtzeNlM9IOI0JnxFqeUZXH7Rllost+m0Cyy+DdUimWB+ke9/d0Bw+aGHyyGsfY8uu/EZ/fi2buf4VRT8NjpimVraJjqFpmNSB+FcRmcAKIxETZZcQ7W/JheYWvE6wlU7/7YvYutu/E2pogzgD92rW3k7PK0pVtM9XR9WpTfZ6r3lKrchISFNdgNIwQXznvuWY/rvFAJLT7rUX4hx/ZYQQn1+h+apIY4WRFEDXjMjndeui90Jp4+povHJMvJS3PfO+Xj6UtGad8tr+Tix6ZQISKG+zwTg2byTddQl+JGpj4pbfbfWh8fZhk23JjlLL7itxP1PxDyuLKlYYiZh44hu1vQL3vPhB+HxKRmORgF1tx1ziuYbmsBOrF9X5Liqn09jMNL6/zeh2ZJuSd6uk/PY0tmB3Y7RmlDQ453nhteuL723Gg8s+wj0vrImlPjxC7RdDc5qSW+f9ZxWWrf2MX56Z4ojE2fvUl/Ym/47Ypb0Ro/OIdV6QJFRw4f6sb0IhTSQu9YeAM69Aa+i+l97kOYZqJNIB0Ii5Sv160bRf/M3z+HDLbiy/9vOSpVHqQNobJIlxVuPhnn0HeYdqGUxOLi7rhUsZQe3bP5Z9hH8s+0h5A8oCaTeJ65KGalvNSAdA9gWJSpuehg5PJ5raJXHPvD2C4ihXJOiZ67p4d+NONDbn8MJ/N0dYr3T1vFLwGfESVfvJ+pvUVrV9O2cMN2KksU5S1v/SihVGIkZoNQ3tKs4L4ob+oCSNazENLxw8zV+DtXSQoPnQcZQUqYOuWt87+Jh2YBWNvGhiAJSvnlyZ3iquXL9NujTRGtipQI+gz0gcZi9e36O9V9F9SImupomnt/GEJ+n4Syl4Saww0k5JZGlvgPAXung9ikHPDNYnCdgvucJ0H6MjnZagRxKQGVFn6rfvlSxNoA6MY0niwMG2PU34wi+fw/yn/xtJGVE5ZKcFWpXSqHSKovloQk/OdbGnke1Xl1asMJIGaIoRQ3nGNZioBl9jCSlEH4Aogp4ZdGD1+YwkNJL77ycmJ1pJ7ZK3XpEGgEuhSHv3C2vw1ifb8bPHVyddFS7BOCOm0NFimO7TvKrofATEOQR89fcv46DrFsgL9yl4Raww0k5JxIGV9Nsz4shMzHEHGYpqwmIt7TVppgk5C8tnHS5L8npuGwYy9LaNOfU634E1DTb8huboAvgBpleERZc3CZH8U/AIE4Vmplm3ZQ8A4JEVH0vllwaB3QojKUB3HOZ1o11JRqAUiAzKDjctd71EdSLDH2ck/pfcnICmXq7QhMIo4U8vfShZOqucNppbcjj518/joj+9Gr7QoHKGmZUT72TqQlMAi8yBlXGOMNmarMXK9dvw9T+8inc37DCYK6WpEpjnS81JGrDCSIdgyk+fjqUc/rJa83maIqo4I6xJINqlvdGo1oXLJxyTCdV97cMr1crl+IwsX7cVb368HY+/uYFwoVKRSsStndEpzmTQM1OrVqg+I4L5n/Kb5/Hk2xsw464lremMVCs1iNyPk4Apn4UVRlKA6WVlaVBDA+EOTtMSsLUf7HuRVS/Go45sKyMXrTa+rURJfw0R5M00nPxCS3sF83VdZrwWgRyKfzUpmkcamlvw8vufosnQ/khR98Lgxo96rWfeBAiAOWPyfTj0alFI/sk2Md8KcZ+ReHzceJRimH4rjKQY7uZNxSil/s7ektPr/qodmb+0lyKMpEB48i8EMVefNDiwKuehUWbhmX68dU9oPyEgP7mv37qn+JvV1y/846s48ZfPCQkCPLOe6maL3//765h++0u4ydDOv3rClTw675hPM2Lw3aA98vVb92AHx7TsAli29jNMuOlJPPKanH+EJUzyI7AVRozR1JLDg8s+8g2wopiWYZtaNL8aTHZNzywjEggreC3J+VBnYBWKM6J5+36fEfp1UQVaiisWRKjcwKqLz3Y14oh5T2HMj54IXXvqb17A1Juf8Ryht8WTb2/E6g078PpHW9Xq5RVGFDUj/2x1CCRttyCLgxgG/+B7Q7pEoKP/7pn3UqHCB+C7CdcFLvzDq9i4owHfuW958TjvjVq0eiN+/FhYoDSmSeCYCeNCdmhJwwehDQdviN8/vwbz/rMKlWUZrL7xRCN5qr4gjc25VDgKBCdEkpaANzAHz4V+J/8OMUlib5rQOcLJqx96g5un7t40q+rpzoGsc3T4NSL6jHjNNCxBPUbNNtVkGUF/ceEqx1+Z+59VuOjofXxp0jBxASBqyniT8Hl3v6JUVhpWm8hQekYaqxkhsrepBf967WN8tiusXqbx3LubAKgt2WO9QI3NOVzx99fxLwlVpKoqWhfe60rWjJBvvnBpXEGsTMZS8KZnCSMqg7rIIOOC30b3vrxWKB9ViM86oSHSq50y5fPBg7kjreNQn09U1hvuHk+stEztnlp9gmkn3PSkVFpTAlHGcJckaqAk0j+8fD0Wrd7IvY5bbckH862/LMPjb9ZLpTGNFUYIzPvPKnz7vuU4586X8e6GHdgiIJREMdA6DvDAK2vxwKvr8G2PKjJI8L3c29SSSkmeViPm0l7OoBPWlCjed1STAOveVPIzfJ1J/Ct45PYOlvX+p9aBE2dE1UxjGlYETSP5+8qiCfViZbFWZpn6MNi4o4FfD8rfOhT3uOH65xkqkMG6Lbtx2QMrhLQ3pqvz8potuOhPSw3nKocVRgj8+/W8FuKtT7bj8//7LA69IWzzFkVokGWcE3lJg0z56dP412ufSKdrq4+iA2vgjRUJNkUrie7sqlIzT3qRawx56gMxmmlCbe37JZmXi7n/eRtL1myRS8c5n9TySa+AkpTW0AvLNKnTXV5btxUn3PIsnn1nk9D1okWFgp7JVSsaXJpWVa6Tmd5wT0fo27xTfqynIaLxSYm1rYgVRgwhuoGZVJ5Qn8yej3AXVFVIKmiH4zSiq/YUJYqvLoCn9ZHPT8hM48ppJYI89+5m/O6Z95XT5yuhl1y5WE6HictMw4NuplFvuBl3LcGq+h2YcdcSjnDKrgPrOpP+IjraZFO1EJVFtEyWGmlplKJPCA8rjETAl299Uep61qoKlg05ZYJteAAMnSdpRuj33tSS49oxZQdHkrbm369/jHcMRmL0CgKmNSNxPHPVLzR/rJMwzIibSiWK4a2LCTPNxwor5oJEYUbdtqeJUha5vB/+602hfJNanh4HBc2Iqf5n2jQZJWl7qlYYiYBXP/ys+Leug5dYWO3SgChYOXTnupl3LcHdL3wQPmfwhp97dzMuuXc53vx4eyT5sx0SzRUUdMDVuQdj9x/IhyV0mxqYiYoRzw2Z0Ixc8Y/XtdKz3uuoJn5StiJOzLy8ZIUq72PWWX5K+wiRzVPUgTWKcAJRyiKlqDmxwoghIosZUTKihsBqGsIxB/SX9cX3PiXkodse/vReIcQUPp8RhjQSlZlGtxzVCTHoNCmV1o1OqPY+gsYW/oTGo14gaifPBEHrFlpmXW9470BQDp22TaNmhFYj2fejqBkxJQxrPT9zc4hIXml7rFYYISLfKaKSRJPoMC2uK7SCKAjPYZVopqF6sDLK8WoBNNuH/GVkrtGT2iiPNdFyVygZqDJp+mO9I3SHzrYzov4yrLy9mhHV+9SdM/LpyYXr9BedVWmq+SaFqTrF4zNC0+JEp7+wmpEOjJY5hnrcEQobbXqw2LSjAYfe8ARWKwWmokOLPZHkWEfyptdeseNJ/9y7dEfiKB1lWffAi9CrrBnxqdHl7o9WpBHByJOJCZ8RE8v4aXsWGbOQcfy3pPJC8Lmqfwx4L9dtRRNzecZwoJG0yG3hfaAIQnrKpEwrjBgiMs1IRPmK8OCyj6Su5ymDZSY5Vl4m3yETg9Hephbsouyl8egb6kusTRFsrmbO7n2qgbe4yZhNTU7tm7wUZx+qZkQpNzOTIK1/uxqyEnVFH/TemagCsUn7jPj+NuUzUnBgZSfUGnMS8BkpRawwYggtlRtdNZJKe60qRM2IIyehm9RaAECW0PayRUy46UkcMudx7GlsaU0vlkN0j5ZdA55mxISfku4EWMzHjM2oiN9Mk9y7RfUZiSz6nkZSDYfVIMZWrRhqJuMRWCUqxpoyuPlw6h0UrohLu9lZxI4VRgioyBWRSLlusvZa2aJ5daXtTSObl8kmMaEZ2b43rxV5b9NOqXQmJ55gTqz2a+asKFH+Evb5TEr6jBhqCvKg6zXTJD8EO4g+HLxPi+Dq9TUTQtsbH23DnH+u9C0/1vmAI9Vo5fpt2LBdblm6eB002k/wOm9NuE3OOx800wjWIUnsRnkElFY5aClGaPuzsHeoLZyK6mvK9JcjKTfH0fMZ0bFZF8rXzdN0el14dv1m3oxn6AZk+iXVZ8REPShCUpKPiWqmiXFpr0pa14Vf6BTM45TfPK9eAUI9SAWf/Gv5MgrfIlEuLafhMGLx6vaC4O3k+xVfW5IkVjMSIw3NLfjdM+9hVb3YclLXjc5eGweuG+j+JDONbJ6FjA1BUoyoTggqy1llEf2SYzuwRuUzouPYSJuc9fPx5eFQjkuguwqCHWdEK+sioZVsGnm1J1NxENPh4EkNHUnztUOHEyuMEFDrn/xEdz63BnP/swon3PKcUHk5TfWqLqZfIuKgRhmYr3yQvMU9r05c7WYgg6zmYBTcHC5qgvVfv3UPzrztRSxY2eYoy6tFM89nxJRmJJANaxKna0bk6kLKR3QyFe0JItfxuhV176UI+pDrsjWs3PS+vPTrYwJT7SQqjIjcd0vOxVOrwjvuUuuq4zMiCSm3tMWwsmYaDVas24qde5tx5P69hASY19ZtlcrfhSv0pRTVACHvMxL8GvP/lvEZYZbDKEMW4tJemboYUn+rcvWDb+CVDz4L581IE5lmRKctBPKU7Su3LnoPm3Y00PNGsL+KYeJjOpKgZ2C0o3q27GXCCUknIa2rIoKb9gpx35K1+NNLHxrISd9MExzX0iJEsrDCiAanzX8BAPDy1cdp5cNy4kzS41/XH0P0fJISuu7EQvpqjPKRBbULmwi7Ouft+nSNTVRxRoJ1CPvnsK83we3Pvo/TDx0IAPjJglUAgCE9Oxkt08zSXspxQwKdfwWMrtAcjZ+N/NJe8/UQ1owIXEPbR4vW9sGSfRF04xgSUyagWDONATZs38tanauM7iCii2khgbxrryP9UhiNM6LpwKozUJto3xbKJzYrZ16cEVVMLgEl5Skyb6yq3xHaMI4WAyZJaAJfZHvTxLAapBQxubSXZ/4MwhTOdZdQh1bTpP8pWmGEgGzYaQcOM+AQtzxK2oeWfZSo85iuQ2YoPclMIz0YBFTrmtqbLGE0kloF4quLK5XexKOlCRZsB9ZoNCMhs0cgH3Y4eH6ZotVqaG5RSieKbgRWlgCuU1WFmHJCePtD3v9EPS9TuK4r7UhM0jK37U2jH/RM5cMgKsIRWMPXpOAx+rBmGkW8Hc9xxAYo2Yl32dqtbKc/uewSh6gZQbT3wRu0dOOMJD0wkwbA8J5A/vO8OCMm7iksmHIypU7O8pUJyme5wIdDW53818W5QIF2XyLbP0iXRTCZyaanntPIV0eoy2uN5Uon+s8bfOgtkvVh3b/uO2hiy4K4sZoRRbwdz3HonVrI855x1cYd/B1C0wJv4iBunqYwGujErwj+JsoiMmYagh3blH1eBFrMEFYbkdK8/UnbcnMTDqwkpZiKYC1rpgHCmh36F6vajYrUg3WNg2iEWJXVSiKkZWlvKN6JJCRhQXw1Db9A6ruYkvZLO1YYUcT79cUy05jk98+vwc3/tzr6glrRfYmCyUkWBd2BmZc0NCkGfuvGGdAdIINc/OdleP2jrcLXEzUjnLqQVtOcf/crxb9NTT7BXNgRWMNlPv5mvbQdHsjXP7hpX7EOKflgpMcZ0fDtoNwzYM5nxNXMy8sb67dqpZf9kCG1bdveNPrIarWC1Wdp7WQRMtOkTEayZhpFgpqRqPB2mBv+/VZ0BZHKlr2ek4C4tFfWo96Ve4k+3dmAqx96A9MPG4zPH1wXOk8MeiZXJeX0pPtY8GY9FlC88knQNSP0ckgT/OadbatylIO++cp0qed4aQtc9KeluOjofaTrkcv575k2wauvptF/4U3XiQVPOBXKIAIee0O8n5tAx0wj0gRUDRwlsa4PlQzWgbUd4+14GcfRstFpLy+NygNfMlvVWkS5b8aPH3sbT769ERf+8VXied1w8P5rJW3YEa2m4QlsvK9vIz4jhGMqS3sf8+x6/MXfvIDlaz/jlp1z/S1L35RODRPfHrSydTQj3r6sGkOFhK9OHLNnFPzs8VXYsbdJ+20hvStFzaiBh0oTRqKAN98ExzWyA2u6BBQrjCiSCziwRuUvxBRU0qZnC0DyGwij4jMiDm/jLG0zjddnxC38G99z4TmjkuDVLoqgZ6ZC7HuFSlqeLQEzTVp8HrzQI7BGU5ZWBNYEBBAv859+Dz9+7O3Qcdl7IpppBGdAkaJkl8zr+PikTZAwQYc209y3ZC0+3roHp44ZgP361BSPi8xPQWeoqCw1SY6jJkJxe6GZaaK8R95yU5E0zGvpH41SaVUh+4y4zHsgtYFfk6Cu42orA4QG0V8Z1tDMH/DdkGaEpj5Xu08jQc9o2pqUaDm9eNvv3Y07sOSDLQZqJMeKddtCx2RviSRkr1y/HdP+91ns2Ksfi0Ymqu6G7XuxkxH/RrcXhDbKI5WRMnmmQwsjf311HZav3YqRA7r6hBER/GvvxezIKqYckQ4TVZ/S7awifgOyS3uDJgjde9d20hU8ZopgD6LuwOtro6BAxi7DTP3lzARUTUFIvcY45znuPedbhi9RJx2YyzYRvc+IyYnGm9Udz60JlBPPjNZC0DrIFk1zMF29YQc3rcgHioxm5Ot/8JuN8yEI/L9ZyM4lpbCip0ObaQqPU+UxeftdznXFAqVRSkrCw/9nj/NX5expamFK72HYLWnEgZX3tCRVyrpfDKwVDNy0cpcT06gEWqKde+WDLdi2u0k51oXOCg7RK13K315aXNf/sSBcCzF0X9eg5sZLVG4HWivWUjCRNbe42g8yanMdTRYhPe031vs1PeF3J3qSf6p+OrYw0joTqvTRFsFJyITnfVI8uGw9Rsx5HHubWvgXU3CpP/JkHEc+eJHBSY7s2CVTF29eLjVPctn6wwFRM+Ky8yadasm5OPO2xTjxl89KB28SKYObpcE2C64WYmlQkoKu0IpCL6WXbwpkETQRNSNyFYtK0Ptg8y6s27Jba5sF01VjKBRTS4c203h9z2Xxfj3mAio2L8Gw8byayBBXB/t46x7s05tvxiJK9xxHQv2vTE6G3ElR00xj0GQkgkh7Pbh8Pf6+9KPi7+A9s74QP962V3lVQFBrESxGZS8OtlBFPpePM0Ivy1OoEiY+MKh119g2iL3hpnq+aZjIWgjL0WXrpSX8U5LuamjGMT9fBADo3qmcUq5A9kGTtmaji+SXBo2Xl44tjLS+vUqaEc+A7brxhpIOkrI+RYVUTcdx5JfSqtggWPmFjslUSL0SUT02ryASpKklh9X1bBt5FEsUXbBNmVTtBeM3rZbNObYDry5RvuuipgTSNgcGXwvhOsU19FB9oyQwofELUr+9LUI2b88nFv9ZWY9/v/5x24GoVSUppIMLI/mXWaWfezt2jjAwBMtg14N+Lm3SKwueycOEz4gs/BD15vIvhUf1rb8swxNvbWBeo6wZYZgu43TMzOXENCPBvtHQnMOjr3+Ck0b1Y6bTjwuk77+Qd5pnnQ9+GZe2mcaEgKkjz9CSeqMZUzetFMj/2/ctD6Qxp7GlXqNVgnk6ts9I67/BBy/iqew107hgqUhdNDbnsHztZ5FI5mki9AIJTEgqcT5cxq/QtZwml/UZacnlnycpfeH+xX1GxK4zCU8QAcjh4kXway3CAoFK0DOWUywtTUuOPpT7VyyEz8+6dxm9khLwujWt7lpBz2hlaU5taYjTQoynI1mtKDYhbGpuy1Nl6wIa2maa0CiZ/DPk0cE1I/l/dR1YeS/r//ztNfzrtY+p59lhgUsb32RNFEb08t++txl/ePEDnDiyL/p0qeLXR+QlZTT6WXe8hLc/2Y6XrjoOnSvLdL9ftFJHxdY9TbGXSfUZUWijlpxYkC/V1tfdEZUkqBXPiQqykuejkifiklNIZhrZovX2/SEfb2xpc+5P00Z5QprJlA0/HVwz0rqaRiGt32eE7jTiOA5TEDFBWqRevhYifEFG1mck8PvKf7yOOY+8ia/e+bLQ9bosWbMFO/Y246X3P83nT/Bgjepp3L9kLd7fvCui3Nv4bFejUrqg4BlsB5Ut04PjO0kTFaSFsXQ2yMrAEkshTAQ9o9QwCi1EPkv1fNOhGQkLcGlYTSMShE8FbfMxW0mdSjq2MFLUjMg/Kq95MC+LJLc3TVoJRgIltrIjL0x5n9c7G3b6/mVdm/8dPE9II1CfQsRGv2lCDtlud+WDb0iWoMZnu/U1I7I+I7TTrEGVlmdzzoUrMEe4rouv/p4sxEZJ3meEck44D9l3RupyY2lN0ZzLaU+oWpoRSukiTqsqpeoHY+Snb1Q0x0aFNdMo4jfTRCdQpGEgEB4gOefj3EhKFNUa7WgNBud9PoXBLg3PzItsfVT2uwFIpgH/ERWfkbDjEb8eOUFnRxfAVgXBy8irTjXTyHeezTsbcM4dL/vMBCa/jNPQn3Mu4Og6sEbiMyKyPYF8vsY1I4EDuZyLi/60VLMUs3RszUjBTKPiMxKMM2KgHiREpPk0DBYkgvENiGZLV77+OmYdEVWvSP479ja15u8XSmVI6WMzopaXd6Cj2dvZ+ZJoEVxNkxSs1TQqPiO/XvguM6S5K5EviTQs7QX0BSyt1TSUtFFpF/QdWNm/d1MCWSa5erNjCyMFM43CK+V9QVlBz2iIPnTWC5S2AVdp5UqM6maTZewsbKzlFbaKmhGxSqbt+RVQrVdwaW8wGxMaBZ5DNMBeTSOSnof20l7QJ1LRCdNbd7LfgqwwyChLOWW0yD6/KHxfxFaemS+X2wcVi0xSea0kjMyfPx9Dhw5FVVUVJk6ciCVLljCvv+WWW3DggQeiuroagwYNwuWXX469e/cy08SJzBLEAj6Tg6uycZFYeSUVZ4TwBriEydqfRsHXQuZaBQFJhMKePd7kaX1UshNRFEvQXVdty3SVmrS4oqtp1O5TdzUNQH+vo3jfVbSP/vTp7NjS/VpjpqWlbIzMgdVsmwtr3EpJM/LAAw9g9uzZmDNnDpYtW4bRo0dj2rRp2LhxI/H6e++9F1deeSXmzJmDt99+G7///e/xwAMP4Oqrr9auvC46e9P4w8GrbPgmWE5KBwIReKpCwPwELlKm/zxZQOJR0IwQfUYE0ufTpvPZmupzOua0tjyCX/jev8mpmJoR8SpR0Q965lK/QIU1I547EalPqTuwmiCK+xAx0yiVq22mUdOMlZRm5Oabb8aFF16I888/HwcffDBuu+02dOrUCXfddRfx+hdffBGTJ0/G2WefjaFDh+L444/HWWedxdWmxEHhHVYZfGXijJDw71nDKEegd0Tdf3Y1NIt9VZDMMJ6/Se2Uj7kgfgesGA2i6X2/FfPa1ZgXRny7w7aTQVt1fxSWCYQrFNI0BYzrmGYakWeR4POiC0vRCII6+TJ9RhLs9HGaaWj3KeLAqlSebnrZr7JWkvz4lRJGGhsbsXTpUkydOrUtg0wGU6dOxeLFi4lpjjjiCCxdurQofLz//vt47LHH8IUvfIFaTkNDA7Zv3+77LwoyRZ8ReUIOrBFpRtIwwX3xNy/g1PnP+441Nufw8dY93LTel5jkzS59f5zrzfgl8CtVuBWSsKXihJgmEgnSlC+YUBf5vFhBz3gO1SIY8RnRdWD1XceukL4An05k6xWF+VFIM2K8VD5pfWYspJb2bt68GS0tLairq/Mdr6urw6pVq4hpzj77bGzevBlHHnkkXNdFc3MzvvnNbzLNNHPnzsX1118vUzUlHIcsjYiMNb4vYkYq1k6aoXpwykmSlev9AuEZt76IN9Zvwz8uPgLjhnQHQPqK9f8m7nbvqviMML7USAUw6kR0qhWpQ0EY8WnIBBJyyo4C+S9IxXLgEv/O14GdqYoAR0vCMtOolBnEjM+I3HHp/M1kAyBZ1b1JohCyG5rSuZomOG+IZlcymhEVFi1ahB//+Mf47W9/i2XLluHBBx/Eo48+ihtuuIGa5qqrrsK2bduK/61bty6SutH2piE9jmBHbsnRz/HykkHIOpJA/3mjNXLlg8voO8QG0V3OyONPiz8oOpaKolp028qZ8LGFb/P3f0kzJr4gpc001HzkzWp5B1byuX9I9NfIcOnvQlQRWLWyTckHEQsZbaZS/pTjUfmMJOXAmqTgKaUZ6dWrF7LZLDZs8A+2GzZsQN++fYlprr32Wpx77rn4+te/DgAYOXIkdu3ahW984xv4wQ9+gEwmLA9VVlaisrJSpmpKyOxNk18N0PbbZ6bJqTiwij31tGhGaPg3HmPXleR5ruZvQz5+7T/fDF/L+y1hGvD5K5DOt/575/NryBmErk/ns1X9ggwn8x8wHRiQVs8W4aBn5GuaW3JYVb8DB/WrRZaweZKJ+9BdPST7iNrj0l7/snr+5VEEXYwqkKPusK/qM1QympGKigqMGzcOCxcuLB7L5XJYuHAhJk2aREyze/fukMCRzWYBJOv8lEd9b5pQnBFmCWHEJVWW1iXp9mPb4IO2auJmVy4pJaM8qdqZzc+3mpugGUmr4ChbKxPjK8lkp7KEnVUVpplGQ6M49z+rcPKvn8cvn3yHn4kCLlxqX3nirXozZQiYI0VJb79mfxwE0boPSlIxv36FDy7pFHLpdX2WokDaTDN79mzccccd+MMf/oC3334bF198MXbt2oXzzz8fADBjxgxcddVVxetPOeUU3Hrrrbj//vuxZs0aPPHEE7j22mtxyimnFIWSpKBpRkhjZvAZ+TbKA32gjcNMU8pELVBxXy7iedqLGh78dCKwpkCWJKI6aAdNVkYGNoU8RH1GaPy+VbP1q6f+SzzP8vGSuYbEn19aK3Sd9NJepdq0pk1pPyWZSEWvN1eHZMqVLVTUbJqkgkB6b5rp06dj06ZNuO6661BfX48xY8ZgwYIFRafWtWvX+jQh11xzDRzHwTXXXIP169ejd+/eOOWUU3DTTTeZuwtFaD4jJPIPqe2t9y/pdLWCnrEQmxjSMVqQvsZ4battz+bViVc+4byICr3NgdV7TO5G0vHUwqju4cH8UnXVdu1lakYoJ1mraUTzBoDBPTpx86DxCWOlWdR9Hgj3a51JJq1RoF3K3zSi2CgvKq0R73nJBmAlaSpJlIzPSIFLLrkEl1xyCfHcokWL/AWUlWHOnDmYM2eOSlGRIuUzEvjtX9pL/zqhmmkEv2xYdVu5fnvipi4Tyxyl0xhUt5Kyas65+MmCVThyv16YvF8v7GpoRlnW8U2mpABn6f2CdH3/8jA1IMlkQxvsmXWmCSOCWhleewzpSRZGRLr8uxvJu0gD+WrragRl+prr6pWW9BjjJZtxiH4aIkJBFP4dUS0w4KXhm2Hky8ynS+5Zd+y9aSg+IyIqVm/Hbs65+GSbXHh7U898wUozNmYTkG6J+1K5kpOWrPZBoZ3vfXktbl30Hs6582XsamjGIXMexxFzn6JoUbxCafR100Fku3MA0iuSCvi1RIFzYMfi0XXo9NLSYsZE1LW6nF9WzsUtT76DF/+7WSpvk8+eaFYWELpFSZEsAoI/MQCxOuoFPTOfZ5RwNWPtRTPSXigOjq6L5Ws/Q3k2gxEDuhKvDT5Lbye8/IEV0lK392qduAXvbaJ/gSWNWItE2/t5XuWk0r1ftavq8/FVPt3VSJxsvenT7t/TrBpalcOti97DscN7B47KCQT0FUyMNJS+0+IKrqbhXJKhSE/eww8vX49bnnyXW1aw3OAk1ticwzf+9Kp4HlIlqqVoS0lPG7cTff5DUU7TVyCKDXbFtOoKT4uThGum4QijSprIiOnYmpHWJ7ptTxO+9NsXcfKvn6fam4MPz9uxVdR/ph56JuMk+uXiFaRoG+GxyCnYz+XU/5zzhAsamtu21/avoPHm64aPSfuMxPvgRDUjsvxkwSqccMtzXP8Owyt7mWWZeCdoX+FePvx0l1LewSHjPys/waLVm4TT+7aT4NRTtz2k30/XxR9e/AAvvienLRIh45DvV6SOWhvllZpmRNVMY7YaUnRszUjr8PjprsbiMdGvR93AUD7NiKMunJSJjJgpxhX8ijVXXuA3oWzvluz+VVOev922o7S82fUwtNpEpKzWf5uj+DT0sGTNlrYyXTlhi/qlxvoqZ2pN9KFqRjTzdeGGJsa9TS2UqxXLENPKi+Ulef3i9z/FnEfyMX8+mHeSRslhMk5+1HYD9RLpa1EFlDNxTSiNrk8R73wKhSurGYH/wThwiD4jITONpk7e1DPPEoLGJUXollz+qgbZZnBh1oZNyssbnI22ER5pHxqZFzmqYEksSHFeoiJYkgu2L5aamYZ+3MiSS5pTuoGoZ8HnL70aTyKt7jsj+6G0bstu9cI4eAVE/7vHT6ulGaEcj241jW56N/A7cJ6SzvqMJERhUPHHj2A/jfpte/Hcu5t8qnwlfAKQOmUJm2l8EOohIqFHWX9tYcgl/qntM9KiucJBhSQEoAKqz0GpxoJty5tIaJO896hK/Ug+IzryDdGx2uDSXp75LU4cp3XcDk22/IoktZpGBdPZivaHkooz0p4ovP80vwAvheMn/vJZfLa7CV0q9ZrO1Fb2pHDVcaK9tFdhBYq5bdb5khAp6irQNliQ/EjEyha+VBuS4BRHmVL3qCCNsKK2CvkQ8IQRB9i2uwk3PvpW6HhbHfjlEMsOzGI0kxCNKP2sdBNH2bcLZppQmQJp9XxG4p28dfOlObBu2L4Xyz78DOOGdhdKFycdWxhp7dUyqrbPdjcBAHYoLn8s4DMNaUzoafIZIa5U4TStabNLKH+OepJXtD/svzffVjONJ4dczhVeFisamIvFytbNCkXRNS3KEIxt4YKtAVTyGWHXgHkWAL7/99eZ5zMOcP2/38SDy9Zz85LBRXjMkba2+jSrHDONpvZRJunfl36Ee178QL0wDloOrBEMNCJzRxITPO29+dzPF2FXYwsuPW5/4vkkfUY6tjDS+i9J/R7E+C6KhvJJWjOijcIXnqn3RSQf38tJMkN5jl37zzeJm/XR8tW9jZN//bzU9XGPMzLClmk/IBG5a+mHnzHPZxwHbxIFPo/fguJT1PUZ4SHiqB0F3/3ba5Hmn9eMtC3vLSLyLicV9Exlaa9CXXzpKc2zqzHvXvD06o3EdEn6jKTH+zEBCj4j73tiddA6znfuW0HcdVYV/9I89YGoLOvEvkTUi89+TlANimgeoqw/z1bKmwRpZpii6UOx6hGF/KCQr2ScXz0hDRRPNaJSBuPDwcStOo56ADgmbng1jc6u39ylvRB4EVnpGY3pCl5nCtpYKTKG6GhGaClT68Aays9/hFZvG2ckIQrd+qX3/UsSSTz59gb89dV1xsqWNRfQSNNqGhVU7lvKFYHlcwD+y8ddTaP45HIikpphYhVGFDReJvIppDEj4DrYsbcpfNTAFgjBL1BZnxGVMqNOG0f3yjggCrUiX/RRaEZE7lmtXfTqyh3XKB9DVjOSFCTbI+PyzzzxSHTKAAJf3BpvcTbiQUwG0m3w7m13Y4t0BEsdZIVAmpXGJZyXIQr7NY+4B5qgVonpM0JpD5YApRqDRAaSZoSlDRQldgdWLc2I4HXqRQhDdWCNeDUNKfv1W/fg5fc/Vc5TtjwvfG0Y+3caNSMd3Gck/ERZg59SX2aokk1AWOWWGKrLCZ94a4N4Ga65YGEieXmFBv8+NHplR22eIhHnQONC7v6omhGFRHmNl3DRTEjP2beaRiFP1w0LozquXyKhwXX6mqjQF0f/ojqwCqQ1HYB48rynhK6LolW4TR04H7yeKoyoV0mbDq0ZkfXKjuprNo6vlqQwXT3Z/CTf2fB5TwPf9Ojb3hPef6TJ5eJ7dkks7YUbuD/XfNAzVl5mBD1yHiacTYMmA3mfEc55zmQkVVaKNCOO4xDbX8gZXesLIuYPB+Pp/UdoWiK7miYhiO+/Sx8YlCR/2iATGizUOkGa9kYgVcV4/dzi/4ldznG+4/uMtP394PK2JZ4FwVR5NUUiZpq4NTH+3+ylvQr5S5StQlTN5bqu/nYSEg7wLjR9RkSFkRi6F81FLmoHVh1UxnZtM00gg407Gnw7vNPyj9ex3k+HFkZIdlpWpzY5mLuUv3XySRsuzHdu2ck/ZDuVbDDaMy/cl7pmJL69adrKNJPP259s514TnADjVFWbC4pHPl4YNh5/sx5PvU1eIimbt85Hu8jGnnGYUOIwO2YcJ1Ub5cmyq6EZawXC5Zvem+bCP7xaXNZLOm+qXB06tDAi26lNOgAa9Xswk5US3q+yOOrhBtX/EedF/YIoakbUSGKZrakyz/39y/wy3fAEyPyaM9kerpl7Zb1Z2/Y04aI/LVXMN4xsfU2bK1XxthEvbosJ6A6s/LRxBv3z4i11+u2LsXK9gDAfqOruxmbc/H/v4IQRfTF+aA/p9F5BBKALZkkq2q3PSADWszDZmYMxK0phUlPBdP2k8/NcHhTclqzZwt1jiKoZKfqMqJrX4tdqmXoUm3eKrSrzaUZ4Qp9CPehRW6MdVKOIP6Lz1c6NU6QzwEC8j599B19I1aW4N02ApOKMiPD759YU/xYRREjMf/q/uPP5NfjybYuFrue1B29cS4IOrRkhWbFZL57JjZZMft2nBZLTnOnOrSGLhDRbX+V84bM2ISzkpXp3SWxaF682xj8B5pf2sh1YpZ04Gbdj4k6Z+etMbISk0v1aUqOnoz8NflUnCW0JdJo1I4++8QnmS6YJ3s/7m3YZqw8pf97xOLCakQBMzYhZTbLnb/WMkwxSIwKtfoN7dFLKT+erV2RjPC9dqsoYPiP54/9+7RPlusS9pj/uoGdyS3vN1e32Z9/HpzsbtPOh1ciB+Z2yoxZOdeq7bU848FtSZByyI7TI7WlpRmIeZ4Pvg0lBHUinZqRjCyOEY2yfEZOaETN55dz4HSHpuKFftPtUjasg/wXZlkB2vK8qz1LTvL95F1zXxV0vrCFfwKElRl+fRJb2wv+sVH1zmGkY5659eKV8hsH86dKI3lJZwjF5nxH29W7g7+hWBkWTL42MQ5ZGRLQeLQmuFJHFdfPOrptbhWrZ5eSqwkiSdGgzDVkzQn9IZoUR/9/KWaeoT8monzOK0kje70Pmi9v7t+SAzzEz3fbM+1L5eUliCV28Qc/kuqbpqn28ba92HrT+4jDOqaLjCyXy1RzVk//Z46vxxdH9MUhR0ymLQ3FgFUFn/E5imcDo6/8PzTkXYwd3w94m/4DB36mZ5zMSPnbR0ftg7ODu0vU0RQfXjJDtNLSXO43SZBrr5IVWvzJVYUSjLrJ+fPkN1+gpbnnyHeW6JKHR0rUEfPTZbuwSdNwMRrd1ob+nC6mMUoRU7yitNKSVTSbh+V6ZJENzYBW4vWbTIVgjprm1UyxfuzW0nF5GM0Y8T2gw0ztHy9KxhRFJn5E/v7TWWNmmxoacpG3eNKzQ2K7rUgdZ1Y3BZMPBe6+VVoW77ElCVaACoLwsVIXCwKPjwLdm8y4c+ZOnMfknYiGwAVmfEXkil0Uo+TuOns9IM+E5yPqMyE42UTbVh5/y42aYgurAKnCHaf9w8xJ1VXnbHCRBxxZGCMfi6q/+pb3q+5QkG2WED1Uzko1HMxJaQi0jyIA9gJVl1V+f9Vv3IG4bm87X9zOr88G9tu4Wc2YM+inwvszTqOVguIxoTWxBwSMKZ+bwh4HR7MPlxfT8aHvTiAhzaQh6JlweT/Oh6YNF6r9Jb7nasYURxfXqJjCpGUkLxAGJUr8sLa4zB9n4HEHNiNTXeoSakUL+cRLrhE94Ttx4GCmDt5WAKiRzgbRmhHP5pfcvb7tW2oNHnriWqtN8RkSEw6TCwaugW1VeeyQRWoBHBxdGwsdMPyOaHc6l/C3LrU//VyN19NBeCtVpSWdC3bq7SfIlZ3+xZjWFkbgo3IFO315Vv0OhYJf0J+9SlexjRXen7JaA9zJP6OXWh/A2rduyJ1RGlMQ10dN8RkTaT8dMmbapm2um45wntUXS3wodWxiRDHqmAjVKpKFyPt62FyvWbTWSFw+uqp1wLAoBXLXppvz0aeyRCODEM+voCiNxD3A6Wr/7X1kndf19S9aGNSOM69NobmSZaXSeHslnxPTSXt+1ktpEFeL60qb5jAhpRlKoDaDB9wnSO0/0GbEOrMmhuuGSCXyaEUlfhiBbdomF59ZF5QWgDRKqt6v7fD7auod/UaEssAewUtGMFIhzLP6/tzYoOxoLp4l4iqXVSdeBNeQzAjXn6iivl4UkYEVBJkOeMoV8RnSW9sashou6PKLPiNWMJAet7U1KiFQzjcG+Flcfkq2yrpBFQsXvw/9brkJMB9YS8xlJcjWB6yY/2MnCai0tnxHiahr1/MTijEQ8ucWmGVH/iNSp46X3r8BDyz9STi8L3wyjZ6ghtVfSr6cNehbA9HhN7zSuwDWlBekuqPem2NAvvvepUjqVYnmCi76ZJp7nXriNNH/dqbhYRn07G7eTA6c5mmWTVtNEuWuvi+Q0I9c+vBKr6sU3h+veqRw1VWUhn5cCeTNN+L0T04wIV4PI5Q+8hi+NHaiXiSC6z0tlNU3SXwsdXBhJx2qaP7+0FtPHD1bOSzVmhyz5ycVfFnPzM9DjjMSFTrwFnvq8XGNpbxLEHfXVZ4rk9AUVQSnqrvXymi2MsnV8RgIOrND7ahd5+6NuK1r9//TSh1L5ZDMOczzLOA55aa9A/0lqozwVPvqME7uFZzLn5J/GmCulNZoahhZnxOTcLmryOeU3z6uXEZNAy3eq0nfMM00wLLjU17rL9rPQ1oy0czONjIksfUMjndUbdmgJdqSveNl5UrYfR60VM+Uzwgv3ThvrRO6vlBxYL71/BfO8yljshezAmiwdWhihRIOPZZIwWURc8RtUnOZoaZIaFmTGI576XNdnJC4KX/Fxj8Xe4hqac+zlwSoOrAkJum9+vB1PtwaBU6EpYC9w3dKKgUHC1ESfdRzmeJahxhnh513qbSyDyp0m7dPVoYWROJb20ihFB9aHln8Uah9WBy71V58nmJba0t64S/T2lZv/j72Pj0rNkuxfj6z4WDktaeKWHXfStprGlDBCc1DlnRcpv5TMNDz4EY3l87RLexNEdm8ak5j0TYlLor3iH29g0TubqOeDLwBLq1ASHykuewArU4wiW8w+5kaIXTPiKW/hqg2ca0uhQ7ShY/IKpnXhRmpCcMFegWZCw2fKTENbuls8T1lW3d4isPLgh1mQv1erGUmQOHxGaJjVjMTXi1Z9Ih6Js9Tf/bwDK/28piwSH633ELfPSENzW4C5KIpOsn+Z3Jsmn596XYSW9jLyV90nqkDGMde3MhwzjeM4xLbi+fA0NufQpLN+usRQMtMYr4UcHXw1Delo6Zlp4nZd+OPiD4jHg19frK+9UljOzPMZiWsVkyni1ow0NLcN/lyHu2irYhyd9zfs2Btx0DOG7xZQ0PCpT9TZjEPcb0eFjMMezzKOmqP8+BufwPa9zbrVK5K0Jo9beqm9UOjgwghpMokvAqtJM018k+Jbn2zHv15rs5ezSk56JY0ueZ+R6ExNv1r4rl4GksQ9gPqEkQhs3EmOuMY1I5KS4me7G/HXV9fhy+MGCr3/rNz1NSOOOZ+RjMPU9GYch2hu4T0Pk4JIvjyj2RmnFB1YO7QwQjTTUI6bplTn6Q2UQFAAyWck4spETNRLe9/btEsrvSyxm2maPGYazrUqwnmS75BO0aTJVNaf4dT5LwAA7nnxA5w9kR+jiCUM6vo+ZTNkAUEFWhyR4vkMWXCLuy8krhmJxGfEOrAmB3H3x9KbQZOWaGlEqVWIC1Z/qK7IxlgTdQp3EHfQM79mhH1tqfSHAjqagOBkyloCz2NnA/+L3wVHM6IpVOc1I2Y6F29pL81nJG6zb9IfWrz7LbHXCUAHF0ZoDqxxPEijDqwJCiPesmVW05QCvKihpRRECYh/gPIKIzxUukqSra/TtUnvhU5fEnr9I9TwZRwY8xlxHPb9OCC3Vdy+qWn3eVNa2mtX0yRHknvTsDrzH742AWMGdRMvI6XvBWuASGudvfAiV8a1U6kpklxNEwUlu5om6MAK+b1pZHBdztJeTZ8Rk2aabMZhrlLLOA6xrZKOLhw3XE2jgrCUdJyRDu4zQt6bJmmfkQlDe0hJqYnazj1lB1+ApO2quvDU26ZU01FTeAxxP4+9TaXRPiroaDKMb5ugu7RXUzOSzRh0YOUMfI5Dvpf4N4GMtThprGakxIhDM0LdtJaTLKVuIMx6hZYsRlqTGHDZk44p1XTUfPX3L+N/n3gngaW94poRtY3yklxNo542vGuvpj+PQF1Yl+hu+Gh8NQ1zaS9tNY2R4oVJvZkm6Qoo0LGFkRjKoHUK1uDrOHKezXF+FQSrxaomO85I+uGpz0vJZ+SXC9+NXZUt0z4qNUv261TDTEMMehalmYbdVro+I2Y1I+yxj+YzEndfSPrV5475KpoRtaoYo2MLI3FoRiiYLCbOF4NlVwxWo+QdWDmDeFPSI5IkcVdXprzG5hw+2BzvUmcdtFbThMLBa4YqFzHTMH1G9DUjxsLBc3ftJZ+N32ck3e++ks+INdMkB6lju3CNrremdVpeX5apQVonfVa10v4yA2gNB8/SjJSWT0Tse+FITFCf7mqUDkyVrAOretq4v+xdTv5en5EB3aql86fF/lDBAXtSpJ275uGVRsoXRfVuTU0tUXQX68CaMsyvppE/Y66M6PEJboGKpFVIEqUl5+KBV9ZRz5eKz0iB+FXZpdU+MuhMvsGkrqu3UZ7IJMLK3UScEZOaEVZOtKpu29NkpHxRXMXvEAfpDR1hNSMJQgwHH1PZXM1IiaymYcGbjEZLLF9OClasjFLyGQHiFw6iltWSdWBNkZlGgLWf0k1gukt7c5rClA9unJGkPRvyqPY9U1p3XncpxQ+BDi2MkH1GzD5EWnb81TSl4cDqq0fgrnjj0x0zxhmoUXKUnjASc3kRF1iq4eCJQc8ilNw++mw37nhuDfV81hPYQ3WuNLUjbt6BlXE+JTOWatc2sakpL24MUBoLBIKk5NEmA21vGpPEopKLoYwCTGEkaKbhvLGdK0rbSlhqQc/i9hkpNWFNBr2N8vy/XVdPM8ITIF754DPmea9iREUYcV2gyVQEVrDDwSe/5iOP6rtkQrMjUrSamcbuTZMYia6mYZTjOJB65+JUyQWLUt2113WTt1HqUmqTbVz9pPD1V4qqYlF0Hj1JSI9ai8TCu7RXZbJ0XaDZkDN3JsMx06RkzEjSgTXnupHMU0k3bccWRmLQjdBX07DLkVpNk5JJMby0l319Wuy/qphSTcdFXN2k4IsVtbCWbORhDc0IITpglD4jvLrqxhkB8kuzTZDXjNDPmzBzmEBV0DYjjOjnQSJpQa9jCyOExs+5ZiVEWr8x2Z9kl0TqINNhmZoRuIl3fl1KTTMS1+RdFEbasWbEZJwR3fx00RVGXNc1Z6Zx2CHhRT9gulaXG6kPFcXbNfEBlnOjcd1Oejju0MIIiTSYaYDkpVQaoZfJt2tvYJVA+52LAJSez0hcZpNC343cgbVEw8EH28WFwdUoBHg5eyd/JZ8RAM2GtISOo6cZ6d+1Cs9+71j0qqkwUp8ghTFO9XGZ6LO8YIyqWJ+RBCEGPYsp0AgzHHzr/0odZjh46zMSO3E5sBa+tKNunlITBguQNEZ6cUb08PuMyJN3YDVlpuFoRjiDRs+aSgzu2clIXUgUHpOqUGHiFWyvvlgdWxghHDO/mobiM2K4nDQQ9hlpj3fZhimnvQL79O5sNL8g7c1npNSEwQLBbuO60e9NwyLr04woOLDCRaMhMw1P88GrXuF8VF/5heek+rhMtFK+Dub7S9Ifhx1bGElwNQ1v8Em6Y5iAH/Le/E1+b9qBxvOkYXoy/O05hxrNL0h7W01TahFwC4R27SUck4E3VvA0YhltnxHTZhp1n5GoTQ1rt+wGoNG3jWhG9PMgkfSU07GFEcIx03ZoWp9ldag0CyJScUaYDqzR3OfA7vJ7a6hi2kwQtWkuNs1IxmpGWBAdWBO8FZ9mRCG9C7NBz1iyEVdzYqQWdI77xTOo37ZXQzNiwmckmqW9SU88HVsYIapG4imbu7Q3xQIJDVkzTQneog/TA0LUzzw2nxGn4DMSsWakxDYqLBAKB++6Wpsu8oRY3lPI+qKeyZfvujBmpuHdjbCZxlBtSLz+0VbltCbk51xUDqzms5SiQwsjJCnbBYw+FWo4+HYegwNgv3iuK7c78rEH9ha6LmmPcB2irnlsS3utZoQJqdqRhqyR8BlRLcCUmSYfDl7dgTWOtz+n4eNj4oMgKiE/6aGzQwsjpNY3v5iGnGG7dO4kfPGxkOn7LA971TxVMREkikTUg0HcPiMiwsJ5RwxVLqdkV9MQfEbSE4FVDWOraThmGr5mJI4RQN1MYsqBNcll7VGhJIzMnz8fQ4cORVVVFSZOnIglS5Ywr9+6dStmzZqFfv36obKyEgcccAAee+wxpQqbJA6fERpMnxEkL6XSeO7dzb7fLA0ONwKrxD2KDjJxtJv+lySNduIzUjTTsK/756zJOH/yUOVySlYzQqh3pBFYOedFBX1q/q7ZvWlY7wFPY1wQZKIcB/Q0IwYqEJkDa7KTjvROZQ888ABmz56N2267DRMnTsQtt9yCadOmYfXq1ejTp0/o+sbGRnz+859Hnz598Pe//x0DBgzAhx9+iG7dupmovxZxrKahm2lKcyBlIeMzkndgFe/8osoI3YFVhGzGAVrM59tefEZEl/bqPquS1YwQggMmqRkp8/iMqGgWcq5rzoE1AziMpuA7sEb//ruuvDwga5b20rNzBT7d1Vj8HZnPSMIfwNLCyM0334wLL7wQ559/PgDgtttuw6OPPoq77roLV155Zej6u+66C1u2bMGLL76I8vJ8iN6hQ4fq1doQpI5rPs4ImRIdR6Uw+eWaJjNNWVRmmkhybSM2M02rvpX3te84epOHKT+FuCGvpklwaa/mahrAnJmmoSnnE46CcIeBGAYAlY3q8tfLP+MenStQWeY3YNhw8MhrOZYuXYqpU6e2ZZDJYOrUqVi8eDExzSOPPIJJkyZh1qxZqKurw4gRI/DjH/8YLS30T8uGhgZs377d918UkDUj8QzY/NU0SXcNMbzVDO3/xbpFN5yeRUawp5LyK2cMbiroxmWgEfUzj91MI6AZ0bnlUtWMhIKewS3GTLnoqH0waZ+eRsvjrqbxvFuq4eBNmWl27G1mCqi8jxKn+G9071JeGJG7Xxdq71/GCY831oEVwObNm9HS0oK6ujrf8bq6OtTX1xPTvP/++/j73/+OlpYWPPbYY7j22mvxi1/8AjfeeCO1nLlz56Jr167F/wYNGiRTTWHIPiNmUYszkrT1To2gvC7y0ojep/hEHX3LlapmJC4fC9GlvaICZpDCfZSqzwhp195CW508qj8G9YgvVg6g7wNlMhz89r1NzEmRqxiJaeBUMdOofeg6IQEssr1pEp51Il9Nk8vl0KdPH9x+++0YN24cpk+fjh/84Ae47bbbqGmuuuoqbNu2rfjfunXrIqkbTTNi9pGQe0179IYOwvMZkUHYTBPD+xSdZiSSbIs0xmTWyAjuTaOqGSnEFylZzQhjb5pMRn5S4F3Nm7j0I7Ca8xnZvqdJa2+aOCbUnOtKaydcqH/oBlfvyZQt9WhLyWekV69eyGaz2LBhg+/4hg0b0LdvX2Kafv36oby8HNlstnjsoIMOQn19PRobG1FREd5dsbKyEpWVlTJVU4LoMxLT+GZypUlaCEdg5adxHEeo0UVfqjiarTwyzUi0tW9sjkkYEXY2VjNNNbe4qCwDWkrUZ4SgGClOMFEsG+d9+JT5lvbKl2/STLN9bzPzJeZ1l4K2LdLVNDn5eUJ1XmnO5UL3knPFP2ary7PY1RiBt30ESGlGKioqMG7cOCxcuLB4LJfLYeHChZg0aRIxzeTJk/Hf//4XOY+h9J133kG/fv2IgkisxLCaZsfeZuJx9q698XLWhMHKaVl1FVkhIHqv4pqR6FuvoiwahWLUVTf19cpD5lmp3HJBI1KqmhESBc1IVkVbpNlxvJoRb1Z/vWgSxg7uxs/ANfcsdjY0syOw8vamiUkzIjtPqMYGaWjKhd6nnMRynuqKLP+iVpL+/pUeVWfPno077rgDf/jDH/D222/j4osvxq5du4qra2bMmIGrrrqqeP3FF1+MLVu24NJLL8U777yDRx99FD/+8Y8xa9Ysc3ehSBw+Iy+v2aKULs6OoRPrwYvU3jSt50THUeHriMfMtmZUwkjUxKcZEWtv1aW9pe4zEsT1TObZjJ5TLy1/FjSfkQnDegi9OS6gFc7ey1cOG8TsF7q7+prAq8mSSkdJMrxvF2qaxpZc6PnI+J5UlkkIIwmr46WX9k6fPh2bNm3Cddddh/r6eowZMwYLFiwoOrWuXbsWGY9n2qBBg/D444/j8ssvx6hRozBgwABceumluOKKK8zdhSKkTq+7muaMQwfiH8s+4l7H37WX3zGm7N8rFIRMhYisDmJmGjgQEQG1fEYM358VRtiIOqbmzTTy+ReW9P5x8YfyiVNKziOMxI1umfm9dfQFw5+fORonjeyHKx98nXpNGszXKnMEKwlrrG/JuRQzjRgyTuJJN620MAIAl1xyCS655BLiuUWLFoWOTZo0CS+99JJKUZESxT55+/WpQZeqMqp5pkCa9vjSkYh9S3sD51gvbfGMhH+B2HUxmGmyyZhpqsoz2Nuk3nEaUqYZcVr/J0tzzsW6Lbuxfuse6bRpxIVbXGGTb7s4XOjb0HEYLeRvQhj5/MF1qK7I6jmwxvD+q0Rg1VmwEBQWX/lgi7AzusxKqaQFvdL8xDNEHBFYabA6cz4YFB9TE69OPt7bCAofr320jZvetM9IHOJ9dD4j9Mq/cOXnUKa6FraVhuZ4HNlktFhqmhEXu0vEKU+UwsdJFGYa3jyoG4HVdc2EsxfZcZfrwBrD+68a9IyWhlfl4Pv0g4dW4tHXPxEqV2allBVGEoT8VaYetjefWmxpMK8vi1TBVOfRyUZ1CCq8mOJBz5JWIrZRHpVmhHK8b20VBnSr1paz4tKMiKr9M4oTb3MuF4pKWcp4J3MlkwlnZuR9lbMisJKqc9jQ7qH8WwyspnFCf4QRDnoW4cyqFA4e9DT8FULq9xKHptgU7eeNViAqzYiQatPEl4R2DnmS7LCianrRQGMxuIwgm3EiCXxGewzFojSLbNAw8cggs7S3tqpcOv/mnFtSgyyPllybz0U2E3/oqSxlNQ0g9n56HXB1KDxTVpm82sRhpnEV4oyoRG0toDPUSJlp2nvQs1IjLv98XgRWkZnH1Iunk42qTFVUyQqWLfy1HcNg5CAa7QhtMDD1nNNmpsk4DqrKs+hSKee61tziRrrLbdycOv+F4t8qS3t1W0JWG0PqpyZClBfum6XJEdWMRIkLtTgjtDGf97x1IuRaM02JQBrkdd8p0fS864TMNGJF8fPREUY8A4dMhE/Zj31hzUgML5TjOML73Ywa2FUiX/LxYiAn4ZzIpM6BtfWynjVy8Yaac7nYNv2Lgy2eHVlVVPK8plBd2gtASNXowoxmpCjkMD/UOHkI+J3oksspaDk0QrDqfGBFZFGOhBKqqnnIHiP64eBF+o7JLwldTGgTdjU044PNuwAA804fKZxO9Ks/K+i8GYeq0YG4E+sjlxyJGZOGCOdLoqi+1nxOqfMZab2fXjVy0Zabcy4xoN7XJg+TyieN5M00Zvswb6Rhm2kE8je0tNfhyyIiuWjXg4fM0toCLuhBz3jPW8dvXW41jTXTJAbNZySOby7urr0CeaTBTFNorLc+2Y6cC/TrWoXjDyFvDeAvU65QUQk/Hs0IpFa2CFeJ6jNSEEaEiyTS0BSPmUbYKVlRGGnJuUSVd4p8nJVRUcnr7nPli8CqOJmbDEDHGhvT4CtkejUND517lnksSbdsxxZGSHvTGMmXj4l3Nw0OrIXb2NWQj6vSs6ZCalIQvVRcM0I4Zvwtc1BeppbpT788ipErxWck8K8qaTPTFPpJJ4mQ1UA+rD1p8kvTiitVoljay9PCegWgYNmkZxk84rpmhJFCWaychM00EXcFkxvlcX1GNPp1KUUpVgp61l4gPWPzu/aS4WpGRHxGDFVUz4HVbf03/zu/E6uA823oDzbCq1di0oyoOrCO6E/2IfnN2WPpz6E4wOrdXFx7uQg/qtb7kRUi8pqR8L2k4KNZG5XVNNo+I4z2F2lTU0HPimYaRlalqhnJMeYV3h3phZoQJ+mm7diakQgcWGn5BuHu2iuymsbQzGviBS9MDo7jiE1Gko5mol8HcfiMZB1HKgqr91GTmvr0Qwfg5FH9uYNV8sOwGOI+I63XS/a/5hayMJKGiUoXUtN1rS7HQf1qqWliX00TdGB1zaxuKmSroxmJg5wrH5qBtZqGh6CvPKVc8UKbDe28rEqHFkaImhEDhhoRFR7tGpmXLRVBz1pvI1fUjIQnBVI9ixOs4E2orqbJOOYFlPKyjJRmxPuoVSbMQhulYSAWQXajPFkHveZ27DPiEDSLvbtU4pD+dGFEF29XDjYh+d31H8xp+EP4y2o10zCjU6uNAybJ78ArB8uBlVdZPZ8R8Zru2NukXI4JOrQwQvoiMPFSiQSXMrG0Nw1fgoXbKHT6jOOE6lVOmG1kJ1hxzYifjELcBh7lGUfqa9I7CJGSFQZ32kDrEP5KM7LCiKwaurmFvLQ3De9DVLDuTHfM8rV/8EMixj4npBmJoyIclJxRdZb2akjZMvXc0cDeTy1qOrQwQt61V0/t6bou9goEl+LJ1kLjqjGfEf2MCl8zDsJ1N7ETqahmJPjiRjFBlWUd5T5CamtebIS4nPJMISw4qpppKEt7k16aGBVR3xUrHHycOALSSBoETpUIrEwHVk5aneFTppa8zV2jpkMLI8SlvdAP1S6SPO7VNPv27qxfIIGwmSasGSEJEoUjW3eLqQazgmYRETWzLmXZjFQf8Ztp6NfR6lrUnAiXmCyyS3tlB9t80DNSfnL5pJVQrA+H3abNEsEGSXjbjVR2kKjkgaKZhjGFipYdpUZHJc6IitNrAb2NTMUL3b7HmmkSg6wZMeE1woe3g6OQA6tEJ/3+CcPp+QjnEmbJB59i4dsbPA6s4UmhTMcDqxXRr2eRpYm6lGccqYHFeylrqSQ9HLx4WWlA1oFV3kxDDrKVhq/m08cO0M6D1A9YzoV3Pr/GeHnM6yNuZta7lfwTLphpVBxY1fwEtcw0EteeNKqfcjkmsMJIAF1BRDgcPK8kgf4noyqM6iVeuX47LvjDq/h46x4AZM0IKUaI7ICmujGd45i/97JsxidgnTqmv1R96Cc10qYIWZ8RWTMeTXvZXjQjQRw42BbhVyvrcSVh+mqPS3v1zDT0K04ayREgBOvZu0sljti3l9jFEdHBhRHCQQ1HIxlMOMqa2PkXMDPJfbx1L4D8yohgfmQzjVyh4n4I8fiMdK1u2232xBHsAYG3mobnE1J09k3FdyEf2b1ppH2KXBCXkrYXnxHSbXy2uzF80FR5lL9Jv/PHom1n5odaoGjatgxRdoW89lxWM6K+ay9LsXz4vj3Z5QqW8fCsyeIViogOLYwQ44wYWdzLh+SA50XkXcpJmIqjHqgLO8IWgp55iyNNNtKaEUFTT/Cqvl2r5AoSoDyTQW1VmzCiuxyb5xNSNOOUyFwrG/RM9r6+/4/Xsash7CSehq9mnXn6V2eNJWbhOMDWKO35Pp8REfNwdFUB+JqRsYO7AcgHChzYvTrayhDIuXJjL8BegcNrc1a/5j0KEe15RVkGA7rF345BOrQwQo7AqpenaHL6dtLib7opscnEl05hOXMxqqbnPkR3uGWhohlxHOB3547TLjtIeTaD2uq24MX8mrU9JxWhkLfaJm3ICgUqQsRdBD+JUjfTfHE03dy3TdDRWwVv+3eu9AflTkK+4y3tve/Cw7Hgsik4aWQ/IU2OaeL6YC3A8hnhPZ9S2ty6QwsjxDgjMGf+YMFf2st/rTSd6D2F6WdR2Pek0KTepjWxtFfYgdXz918umIh9e9dolx2kLOugb22bxkVGwCDGGRH0FSkVM4Ssw53K5nCbdzaEyy2R9iHhFdijNNOcMrp/yLThLe6qE4djxIBa/OSMkaFzccF0YHWAqvIshvetJQaIK14XUd0AtQiseqtp1NIB5j5Y46BD701DCwcfi5lG339V6oWIelDxmmmAQtvm60fa4VZ+1145p0idsniUZx2cPXEoFq7aiKkH1XHbVjsCa8noRPLIDp4qgy1pn50SlkUCGr3wjXxp7ED8Y9lHRsqqyGbQ6Nk00Vte39oq/PvbU4yUo056xjUSKqHd8/MKZTUNJy1z7yBOaiFzUkrklQ6tGaGFg9dRjLguMKwXP6YHTZAQfbmOG95HOvAODRODeFAz4s1SVJDoVVOBkynLy2jh10XiIpimLJNBTWUZHvrWZMw6dj8pVSnRIZCTvm0JrFQ1E0N+rxP5G0vr0l5VwZEdf8bBj049BL/8yhi1SgXzY5Qdx9J4HjLDWiJmJFd+lQMrxZUn0sMuAKWjEdWlgwsj4YdsIhjZAxcdjqkH1TGv0Q0HP/+cQ6XqGnV/ZvmMkJxPSYPvIf27YvTAbsT8vap/nwkogRc16APDFUa4PiNO63UUFB09k0J2AlMx45GFEelsUgMrCqqDvC/HqWP0Y5i4rhsqgL20l3Qs2oZmDWvCgkqEdXxv067ix5corNU0IwZ0xTs3nkhNyxrj+B9CKVF7CNChzTTEQVNh2VaQPl2q8OVxA/Dk2xuo1/C0GiKqO2OaEQN5tJlp4PsXoCztpfhO0F4uWpyRTMbhSpCmh6UyiU3ygIBmRKMypWKukXdglS+DtLRXJzhU0jjUHxGXFTgitlt4tLAm0OCZYH3jENiffHsD3t+0UyoNb6SmLVEGOFozzXLTRIfWjBBlEeh5ILcJMhxbnmYvyTpyUUB56I7jbWYakmZEv5t5v569tx38aohDrRwUjGSEBFackZoK8reBE7gu7cj2JVNmmlJWZzOXbxq8LRdhoY2dfwJmGolrk3rk72/eJXU9e2kvO62OkC3ywZoWJ9cOLYzQNsozgar6THTlREZSM8KLsqg7kAeFEW92okHPSBMMKQ/WXi8MK4gxQv4rXDNNG6xxJZNxcO+FE0PHS22OlTW7qJhpSPuxJGGyM4W36lFrwEhmIOq1gZPnHTHUcG3CmBiD09YTdIKe6QiqIkWmxZLTsYURwt1/vG0PPvpsj3KeogMJ7/mL5MKavOOmoSlvpim8HF5pXmclDC+PYBrfz4hGpKAPjEwxxF17PX9XlmWp5+MYYE3M57KCrZKZpp2tpmF9/c6cNNRcQYQhI+OXhKg8f8WxmHPKwdEHPZO41rQ2LLINRaFuMtEJjZCeGYJPhxZGvF9ShT9/98z7WnkWVF687sOLwCpUlowDK6NGDvQnOqaZRtBnhLX/DO2FzIYEg+hnpKBmhDcgyviMsBwGTQy8vIHNhJmLN3YO7F6Na046SKtM0tLeNLiMqFokfTK058fEYT1w5viBWnUKwhTgGfUa2L1TLKYwps9IxJ/xpzACz+nANNNwxiz28xEfe9JOhxZGvC+WaV8DXn7UCKyFziVQHVMOrCYoxC0oaka8ZhrCCE26PdZESWtPkYHV9PAZWk3Dud63moZwnvdhalIzwttwkDehf2nsAPz0jFHMa3gCz8/PHI2vT9nHU2bp+4xMHNYDw/t2wRd4G5dRoK2mGd63i9H7cuGGl8N7HVgFzJ4pkPmK0OpCa7Krv8BeRhvVkOpqxG1lmh8V3QHSSIcWRrxjpqmvqsKz5w3INEGiUSKsqpQwwvQZ0Vdx7yUGPctDnAAJBdJiiXjz5R2PY6AMNrtM25HNNA7zfLEvGbg5njBCm/hOP3QAvjftQNxw2gjUcfb7icNMQ9KMJOkzcv2ph2DBZUehqtxvZvvKYYOE0lMjiUZyTzKaEZJWM9p2Zg1rodU0ElUZ2L0aZxzK1jJ1qgibSU2g5cCq0d4iM0RaxJUOLoywJwGtvDWlGxFzQ4pcRtDU0mqeImhGiBvlEfJgCXBUM03QlzSG+eizwD4h3GelKTMW7t3ErfH6Oe3ssJ6dMevY/VBTWcZfdi5ZhqkluUn6jBSF8MBx0Tp5m8Dbv0zv1uu6ssKz0eKFMOWYH6RPl0ruxH7upCE4fJ8euIoTiEyWqMLBc7WyVjNSGng7pml7M29A5sYZEaiPqY7mtP7PBKI+IyRYO/PSsgi3s1fdHc1IWh34+pUa3EnHOOmLwkgMMwOtCJmiRXdYbsvbbN9LgqIpjdEfg9Ain3pDtUeh7QkJgxQTUVKY8IXbtCO8d1FwN3ESnSrKcP83JuGCI4eJV0KAfAh5xdU0Wkt7lZPGTocOeubtmKYGssKzJ63U8WLC34MU+IkG6+5MjndtQc/aMs0S96YJp2UJLbQXUi5mgjykmGrTDvFH1+V+nXiv5VxMOk/76laBlwdtcA/uhsyitoo9rAQnbFMTbpJBz9qcjMXTlGczIcdvwG+qlXnHRQhqRk4Z3V8+AqvRGoVh+lYETjVRzNqkFZEO4nFwJ/HgsvXUMZ9XI/bSXvGP3j5dKrGRIKSlhQ4tjPhDjJvtpLwBlhsOXqAMoU2QRDF0+6Q4I0GHTxokoaV4TtBnxPRGedmMg1xL28Pq0blCIQKreIRY0mDZphmRKlYJkTJ4A3qXqnKpMof27CR1PY0kv+oL7SZjpqmgCCOFZfJANMv3vWX99IxR+HBLWwCv4LuS9qi/W3ZJmrEEb8e0FvKuF9Yop2VvlMfGu2qTdktpMeV0cDNN29+mB3qeAys3zohAfeRsq2q+ArKQ9qYR9RlhakYo9Q/mLXsfl089ABOG9RAul1gPSqEFDYFPM8JxCCRl3+YzEoOZhnZcouguXM2I//ekfXvi25/bT7wACqRnc9HR+2DFdZ/Xzlu07PBKFTo0c9YejzASxYo5n/DrcPoVqbsb6IaVjPDnbAdW/8mgrMYMYeBEv/x7YPdq6TQm5gIahwzoWvw7DRtJsujQwoj3a7vwoES/4qm0vkk8lfG9L6/VK8cwpvqp+N404WNMnxFKTxWKwMpgzOBuzC/zoLBDUop4B8Ap+/fCL84cjVnH7otHLjkSgP7eNFnKRBcF1FUdjOWfQWQ1I47j4OtH7sO/kANx0RYcdOtUoZ23eNniZkNv3/L2770eYaS5xbCZBm5Y+yEni1CPyvATxvJwOflLxlQtHmla9Q47U7Z0YMG7X95GebSPqYqyDH591ljp+iRFhxZG/HFGwsd0ULWDF7y4Rb6CTalwzfqMkDQjYnFGVDQjNZX+l5+seaBmyyVYLum5eg8N7N4JZ4wbiO9NG46hveSjObI0I17OO2Iozpk4OHS8e6dy7pbkzPIF6sVrztpqjmaEcCwYvE4F0gdAXB+DhX4X1oyI9Wnv33ub2uyvccQSyjCerUlzhTcrU/uhSDWPIy5kqN62bLrBPToxN8kD+BqN288dRzz+oy8egrpa9jL8NNGhhRFSYK5Gya2haaiG8C3Uw7iZRqk28hRv2yEc80CedBk+I4FMfnLGSBzSvxZXf+Eg33HZwcABe0AL1p044TGuByTDWzN8RrycPKofvj9tOMYN6Y4bThtRPF6ezeCbR+/LKoBXAaVkXmQ1I4D4iisWpBx0c3141mScJBDIjOYzwrqe5kC/N0KfkaADa/7veEYHbykyphgWkrJIosu/SdwxYzz3GpaW3XFA1fyJmgzT4THS4YWRtsdTYWBnWcCzmkax18t8IBrb1A/mdgAW9Rkh0btLJSNf/+/phw3Go9+Zgr6cAFwAe6h1HPbLGAzERvR/kVhpQnQ5CU0OfkhLex0H6NqpHP+4+Aice/iQ4vFYBhbOPVZkM9yvvSA6+28UYO2IrEqvmgoh022G8IxY5Wcc/3J672VenxHDVppi2W3lspe7EgU8xTb1ylWs8YYlfwXTyTpfio7LcUXz5RXzvWkHsjfX1FhpkzasMNKK7ODJzVsxu7ZBjX9tVMGBdGgz03iOCe7aSzI7BPPlwZvYSfVgNWOwX/DMNKT74q+mYVeU1KY0dIVKWhHBCYyZh8OOZEl6LiY0IzICzf98/gCh60R3tCYoBIm/C7TkXGpf9ZlpTGtGEG5/f/BH/7moxg3WXems7mDVNw5BXVYAYF195H69MOvY/Zj9uowx0URpcouCDi2MOJ67Dw6GfTVtbao+IzIaFVPjFE87IAMxzojACpT9+9SEQmn78yW3S3BilB4MeF/5AWFE10yjsilWQTljYsLmQQ9LLpEH5Pu/iYGS3M3I+Yr6qIiq9qmraQTvi2amaTa6fj9MfjWNZBoD5bquiyP360U+x0oncS0pccrnYx8FjRyrD7EDRZbQzaKDCyPehxVUx8tGkSxQEOpVAzBlSU4XFNK0UV4BkrpapC147w0tC1G7KLVcsG3UITMNTzOiMADwNCuFPuGPb0IrR69PCGmThK5hXRTNIEnc94dhJhHLU+zaQv8k7ZlDzZtSn1PHtO0ca1oWcd1wxFVW/zNppvHVA8Dd5x+GYQQnb5lhTepauJEsj79s6v7Fv6XHH4EELMGeZUIUfU5pmUY6uDDS9nfwC1j3K1RdM5L/VyS5ibDJ+XNhenSuwMxJQwhnOOUUNSNtx4gTeOA3b8BX+cIscEj/roQrKRUJEPQlUhIyJZb2ks4X7knEDCHTJ86fPDRcvkA63jWO4xCXQEcNeWkvGdGn6MARi03Res22PU3s62jJPWVce/LBxb9NR2AFSKYk1oQWnZ2mPJvBEMKyepk7ljHpyO7L89qc4/G1ycOo58syDlbdcAIum9pm8pNvLnFBl1wH+c1F00oHF0baHlZw0pGNslmgkKWqQ55MB4oiOqO3HocxgoGx0nn/BSixOSTV2aLNQrru5umj6ddzBoPK8qADKzmXAlyhinOMdL4gGMsIyK9ddzy3/K9PCcf2MLV7LC82QhQQnyWlMBnNiMjXdCG/hia/KoPpHOrz1Wj722uuNP+Ou75OwNuvJarpjB3ynX4u5MAqWa7M+Nq1uhw9a+gxaqrLsyHTsmzfFjIBsnxGBDQjE4bmx/HjA9tYpA0rjLRiSjNiykwjklrGTMN6+WnqbRXJupAiExjwwtfx1cFeqD4jAgNpny50/x8nb6ehEhRS+WaacB68pYq89BmCmYZ234WSunbiL6/lCUa04yLRfJNwlpPya5EQbkWc0QvZHTu8d+C4qEaPfDzqCKyAvy1CfZWobdJ/toXbIt1eVJ9YLuSFK9YYuKOhWas+PEgrE4MUNCO/Pmss9u9Tg/83fmDomvu/cThW3XACenaOPvifDh1cGGn7O2h7C/oKyKJqppHZh8SYAyvluIo81ebI15ZYRLDj3W9wp9y2dHJCDQlWMxpxYPWZabhGjtCRQl8SaUfdfSZo1eMJTEF09tNQhXTr9L4tqBmB2GqaQn6VZVmcNaFtVZjoMECrj/EIrG5AIwP/77j8B1hj1/VfPIQhbLNVI6z2dl1XSmDl5Ue8XrJ3s64mfdgFKWhGThndH0/MPhr79akpnitqqTMOqsqzdjVNmonEgRWFcPBqdZLpMLL2UvYFgXpI1qUAKRw8OTZHMB25rG8dsy/e+OHxwmYvlTqz2lFMM0JWt7fl77mWUAZ3b5rWvmgiFoe/XOJRI/kmuIGuD7oDq3h6kWv9Y4l/shcth0QUmpHgnlzeooOlyUQ07ltbJdyuLG3h2MHd8c6NJ6KbgHYvtLqGtwGp5PjAupq0Gkhe2OEnYL335YGJhrVtQ9p9SDq0MOJ9NmlxYC3uQ2J4NQ3rSlpVVe6BpFoUWeVAK6pbp3JmRE+WytkEspoRUvG6oa8Lz0Ek+JakzBmCqhmh/E1D1UwpwwmH9MVXD2drIWjvkWj9HAgO4p5LWE6FviQ+oYBcho7PCG2ZeXg5fNvfwTGF2KaU5uhTW4l5jD1nfPXg3FZ5NkO8ZmB3v8Or1AeZ8JVtsB79SaPCkXlle73I9TJLe72Xhjf5lKlZ/Mjv6tOO8D7kkAOrqmqjFdWv2EI1TJtpeC9tcMIUtZUHITqwyo3lgePsxCGhhl9UCBkzje7SXqJPSEBtHqTQl7zh8mmlcL8M2afpz4GjvfHn4XAcWM2MivPOGIn7lqxj14UmXAk7sIop3v1bS4i1lX/iIF+js5qmsizri+ZKKtdxHNR1qcLwvl1Qns2gS2ivJzlEv74Ld8W63CsY/eXrE/HWx9tx1P5+bYRM66jIdbKrEE1SaBvlpb2B32nXjHRoYcRLWuKMyAzU5jQjpK998R0uvZCWJhPNNMFXRcFBlXy9fJ1ZzSgUDp6hGg3mT6qftwjS+cIgYmT/Fp/gI66Cl2pWJ56Br6o86xOiZeZtYTMNxPqU9369z0lnSTqgF4G1sjwTEkZIHyWZjIPHvjMFAMEHS0LblD8niMjD8lwyeb9emEwJkhYlosKkUALFy1l9lbWfl8zuzGmgQ5tpvJSX+Z9UYg6sRTONnzevnxYKsR21s5nKhFIQwnw7kgr5jOT/Daph+Ss3Ai+caEU9aDuwer80FcoX1YzE4cDauZL8fSLbyszNvSjH7z7/sOLf804fyS2jsoysyhcpy3TQM+8lomEBaHvTeJEJohakqkzM6RvIPy/RbRtoOJDXprI3zJNPz3xUCu8Ga/wxEjVY811ijQnWZ6REEfkClkHdTENOV57NhDQh3TtL7I4aeA8PG9qdebmqE2Khv3sHJRHBrHDFtEPq8N3j24IIHX8we2181O9X2IGVfT3phecNgay9QQBvBFb+FzevrNrqNmGDlEV3yg6g4NQxeKmKLH/sgX2Kf4t8DAjtGUO5RFwzIhb0jK4ZESsnKs1IEJXlrSFoGTgOc3Ityzg49/AhGNSjGl86NLwENYhp510lnxHJc7JtK6YZYZlpAg6s3o/AEtOMWDNNK8FJR8RZkEShw6tKoazdP4Pj0h0zxuOEW54TrFdb4pXXT8Pi9z7FKx+8yq6LkplGTDNCS+c4Di753P646Oh9sXNvM7pz1sYHc5bx3wD4H0tBj36ekKny2Hmb0LX5jOiPJieN7I/6bXswfig5oB1tTybZktm7iQqkFxRmuJsQUs1/4nYaWv8d3rcLVtXvaM2v7bioidffV8nX6PiM0DQjMtoLeTMp/dxPvzwKpx86sHWJLT9jEWFExjlcpSnZZhq+xlePwntPv4IZ9Czw22pGSoRgx0rKgZWVLPiVNLxvLf50wQShfL0voqimQm1iDR/jLYcFoazybIYriJCQXefvwiVOaNMOqcPEYT1CUUp529STSucNgn6fEcL5os+IQJ/krlJwcMtXxuKrhw8h1rWOJoxImKIc6A98os9Rda721m/qQX2o1zkOSxFAFiLLfY7GYvdBa68Wjb1piJoRV9L0QriUqS1gPPcs5UOLhsizlXLiV9CNsPqxzPYDOjBX0wQq4f0VHC7sapoSgTQZqlCIcqf64GmCggPyl4LKoE0aYEkvvjnNSPi6sEZDscFCGallE2TmEUNxxL5hhzm+AytRHGGWRQsFXiBD8Bmh3SZ36S7nfL+uNM0I7x79sIOeqZtXgqjqDUw4sPoGfs+PrKCZxnuKdp2OqYLuMyKTC/9Dgn21WDoSpn3iVDYdlHZgNZh/AdbHY9A/yf/REPzgS7c0YjUjrQQHJxUzzZfGDsBXD89vLqf64NvMNP7jjuMQvwJUihFJ4ziOWjj4ojDSdozl8V1MJ11SIR1bw0Ii6C9DGvNo9851YFW4Ee/kRQr0NKV1OaPotvdMOBXsSxFG5IoQ3FyOk4cI/CBX5OPiDqzsvVtI+XnHjtEDuwqVQ4/ASp5BD+5Xy80z6HwNACeN7Cf1rsn0Z55GjDSpsh6fkCZDQmBREeyYwpWAoMZ7D2RXagWxDqztkOCDUjGz/O/0McQvW5V6iGo8RGvpfQ3zy3ZF6iKYOSGNP84IewIn/VZFJJs/fm2i/4BEGPEzx4Ud71hfI4CAmcY3kflfySdnH40D6rrkzwkIdVwfCuqPPKNaJ89g//dF7eTUgTcpiWBKpUx7jxxeQ3jO0M5mKOY179dq58oyXHPSQZQ68DUoNDMEqSv0qqnEyZ5AXN48zzh0IO678HB8edxAqZdN9jGwspZVNouYYELxkRg1VtK0SPo+BQ+RVlYJuysVx1L6NSEHVl96+Q+1JLHCSCvBB6e7tFcVWrG0fiT+Bdn2JpKSkN5TtTgjBc2O4zlGyDv0W+1NIWmQpNK45C8wmtPiMQeG/Qu8dadFvRSuTwDvXhMi6n+dsgCgR+cKvHbd8Vhx3eep6YRUyxJfbMRrRM2P3jgjhLvX14zQr6WZroL3Xl1BMZf4yiGXQYvASqrTRUftg5v/3xji9X27VmLSvj2Ryehvc8fWFjDOSY4nQg6sEWtGWIKAiAMrSXNxQJ8ubdcL1GH73ibqOdZ7FjyVdp8RK4y0EnxOqqtpdGmbzP3Hae+xim3dEUiXH4TF8vZC2puGbNoIfHmr7uXDOic4EQ3u0ZlwTTjdkJ6dwgfBn6h52orggNKliuzKZSTomS+2RTi/6vIsunJC8HPLYEzgMnmwmDAsvxrIRJwRtl8HXYtIMyP496ZxhCZM2QispMsdhz5u1VSW+64TRfYxMh1YI/AZkREv1Jb2ik/2JEjCwoF924SRxma+I8ukfXpRx54grKCG1kxTIphyYNWlzUzjR8SJjoX3xRbtlCqmKtLeNEQzTei3mRdFJJdgdb79uf1wzsTB+NrkYZI5ha8kbpTHSR98HrRYHyLBtHgDOOuRVmQz1DJ4QkyoHE3BiZf8D+fnV5EFzY+m8i/m6dB7Aq2tg6ueqI/EK7BLakZklyx3rmzTzmhrRpjaAvo50xs9yqLkM8IRVHnHSB8QPTpX4JgDe+PAui4Y2L2aW4euncrxzPeOFbrW8fUp/7mjD+zNTZ8kVhhpJTgYiEZRNI3sC0sbfL559L6BI2wzTThfNUm6aOf0Bj0TuCdVoV3FLhq8pHNlGW760khM3q9n8ZjMY9B1YA1O3NUUvyP/ahpyQTynP1ZdqwJLQX94ysGei+npQmXAEdqPiJkHp5CiBsBVM9PQ8r/25IOxv8c0xsqELoyIaV280PrbPr3CWjvW9TS8vmwy7zV5wqVcyxkzRHwsZDG6czkBFc3rGZ6AbiTn/bKMg3vOn4AFl01hzjPB7Em7BDMJZNCnSxVe/+HxwlqWuLHCSCvBl7siZWYaWb4wsi+OOqBNEubtjxJE1Mk1iGjQs6gcVmmT9IGtTqCjBnb13b93fOLtNCxSCxUH1mDz0IoWEep0tlDvVOE3D53n0RTJPi7doGe8Swr585cyk3OiNeUFRw7DrV89tC09QzNCIzTBCJla/L8e/c6ROHlUP9wxYzwxLTneDb2mXgFXJpCatJmGcY7Uf3VX74bSMyqgslWCbJwROMCPTx9R/EnSjBRWxcn60Fxz8sE4hqPd8GlpCY1RW1VuxNwbBUrCyPz58zF06FBUVVVh4sSJWLJkiVC6+++/H47j4LTTTlMpNlJqq/028uQcWOU6CusFO33sAAAIf+kJ5qNiOilU33sfIrZiVXumqG/NPV87DN/53H64c8Z4310F46+01UepOkrpgu1DG6RE/Jh0BneaoyUgN3A6jgkzDTt94bTq0l6/bd1P0KRJaxda0ZWeJbWirRBsrkP6d8Vvzj4UQymaEZHH4b0PrzDS0CQecEP2KTLNNBH4LMg5sOb//ePXJuC0Mf3FEjGrTNYaVXriu5DGc1VhoKayDF8/ch/+ha3QiklrvBHpGfeBBx7A7NmzMWfOHCxbtgyjR4/GtGnTsHHjRma6Dz74AN/97ncxZcoU5cpGwQ2njcBJo/rh1NEDfMeDwshZEwbHUp+2XW/1O8ypY/rj4VmT8fCsyaGBkydo6JpNvAMP2TnVX0DUwnq/rtWYffyB6FNbJThByamnSX8X4PqMMJbRehGJ18IrzFfXwLmu1XSnVd8Xl8CzYq5CEJji+A7WusKO2HUO8rE5SNA+BHoEIgfTHglrHxFuvQT8sLx4Baq9zS2MK/nlnEhpDwc8jVgMkyCj/xd8Ro46oDfOaY0HxUPWTDMsIDySXlmh95iSv/cYMfq2gHY32Pe9GvQkkRZGbr75Zlx44YU4//zzcfDBB+O2225Dp06dcNddd1HTtLS04JxzzsH111+PffYRl+zi4NzDh2D+2YeGAgR5v0K/c9z+mBvYRfTu8w7D9PGDjNenGPRMMx/XzXfGMYO6oXNlmZK9VCd8Mm+gDWs0FDUjQWW3iAlA4CWVWnHg+VtlG4Fg+9AGdBVfBF5ZXoKTqBe/EMMvnKXhE+lXopOzb2mvRHdl5e9zinUcDOrRCQ/PmoyvHCb2vveskd/GQLbbiAhTXofNKkXNyCH9w8HVjj+4jii4OhybVtIOrGp707DMNG3n/nrRJHzlsEG46kR/TBnSuyKqGeF9EJF8y7x5UzUjgXx/OX2MUH2iRuoVaGxsxNKlSzF16tS2DDIZTJ06FYsXL6am+9GPfoQ+ffrgggsuECqnoaEB27dv9/0XNcHBwKsZqSREMjx2eB/MOnY/4/VgqTILA0AXyjbvLEITgBM8Hz6t8vKSlvaSQ6jnKTiNnjtJ7EuFmpHcKcr1al+q3gGLFPWSv7Q3mB/5OpEN2LgOrL5y/PlRd+xl1Il2rarZ7YxDB2J43y6YvF8voTJ5gbHoq07If4eua/13zKBuGDeku+8c7bF6hbrte5uo1/mehWRvFRHwm1tcHLFvT+zTuzNGDmiLBNsgoRk549CB+NGph+DR7xzpKSf/kSNarwIJWb6LeIUzVfOZF++pCcN6YN4Zo9CVEEE5iI5Q5g+UF87HO2+JvoIqe4BFgdSstnnzZrS0tKCuzr+te11dHVatWkVM8/zzz+P3v/89VqxYIVzO3Llzcf3118tUTZugECDiM6KiOeBRfJkJHen+bxyOX/zfavzP8Qd66kCrW+C3ZFV3NogPWF4K8SmyGfbEXjj0h/Mn4JNtezGohxkPb50lwiY0I8Hdn8XKDWp3yIWTBJ0gqj4UANCjM8tMI6eVYWoeGHX8xf8bXdzVNeM42s6WtFp46zesVw2ADdw8Q+8U5e3z+gx8uquRKoyKCkS8tMVjgd9NuRweuOgI5Fz/+yijGclkHMyYNDR0nGp6YuRFXPZuegMaBuZ37RUok9BSOg6kvD7jHSNEBPE0EamsumPHDpx77rm444470KuX+LKkq666Ctu2bSv+t27dughrmSc4eJYLDPxRUKwG4cU5qF8t7px5GA4S2JeCWw7n/OadDUovb0Gty4pI6aUsm9ESREQdWGl4BwtVG7730vIy+Tc9bKYhX0fb+EwGlrAm+oXE9eeAY+TrTyQHE3vTDOvVCX/82gQ89p28P5tX/c1S04u8HzWCWkxZTdKQnmTHVqBtL6NzDx8Cxwk/iwaBQFuqxO3AGq4A/ZRPMyJYFXbQM7X70Xo3OOVX+ALuUfJIqTQipRnp1asXstksNmzwf0Vs2LABffv2DV3/3nvv4YMPPsApp5xSPJZr3TqxrKwMq1evxr77BuNhAJWVlaisrJSpmjZBJ0LfQ6U8O95gtOQHx+H/3tyAax5eKV0f0QA9tG4V/OJQ+f5QCRJU2PXVZ6YhOtuZeSGCuejkqvql6r2XimxYYOAv7Q0KI+TCSVvCh8rinGfdV09BnxERTER7lP3yJLUzXTPivcbxOfEN6tEJFx+zL2o5UWhZ78fvZ47Hi+99ipNG9sOfX/qQeA1vGwESD3zjcDy0fD2uOGE47n15rT+/1ga7c+Z4rNm8q7icPUgjZfM9E7DNNMlOgiJ73QRh9kHF29HTjLCFDa9GX+Qd7C5gVooLKWGkoqIC48aNw8KFC4vLc3O5HBYuXIhLLrkkdP3w4cPxxhtv+I5dc8012LFjB375y19i0CDzDqCmEHFE5PXtPl2qMHZwN24+vzhzNP7nb6/58zasvQwKJ6al41vPORQDu3dqW03jNdNEqGQKBz3TMdPITw75Mtv+Ji2/5ZnzgpYdqmZEYBNG7kZ5jPbp25Uf4VGEvM+IiXwccN8yxffEH2smnMkVJwznlkWLjgoAxx1Uh+MOqhOuomi/nbhPT0zcpyfzmsqyLIb31deesiCHApBffRIvKhFY6ZVWvZ2soCmXt5qGJGyI+Ix4382HZ00WqkscSHtCzp49GzNnzsT48eMxYcIE3HLLLdi1axfOP/98AMCMGTMwYMAAzJ07F1VVVRgxYoQvfbdu3QAgdDxtlPtiBZCfqsiAy9MAHNSvFieM6BsSRkSleFMyC2kOk8k7uNzPFw4+wqBnoXwlr/fHlNDJKQ/JvCcbiIw2AHrNNKrtx0rWv1WrRUzn+yLjF25EMyJwDV8TpP7u8mhW+dT2wJtYdPJLCubeNAlrRvzCo1hdWFepPrNOlI+Kh751BH6yYBVeen+LUH14PiPUDR49h1nmvriRFkamT5+OTZs24brrrkN9fT3GjBmDBQsWFJ1a165di0yUn8IxIRJgarCArwOvKfLOeuHjKiYSX76hcvy/RV6jvrX0yYmHX8sQnTASMtMomlqC16tqRioVHFhDQc8o1wXDtZPQMdP0ZQkj3JL915p4viKDPV8TRDuuX8HmFrF3VORVFnFO5pECWSR5nxEGmrJjCNXb6dWF7IIwdnB33P+NSRh65aNCZZLKTypYpwnk14gCuOSSS4hmGQBYtGgRM+0999yjUmTsiKi7HMfBTV8agR88RPcJ4X1F0oSOOL3MafTtWoU/XzARs+5dhm176NtYk/Db5MMY8xkxOL7pBKEqoKIZCX4xUn1GBBxYZXaxDfYx1k69Uj40DnsrAdGuLeQzInG/XrxNrvqqiWpGRK4ihQ+IimMP7I2nV2/CSaPIwct0kA16Zvq+WV3GVXBgZaE6RrB8s/iwNZQVAj4jad29V0kY6QiIPi7RsNU0cm7ee3/isB7Y29SCAd3ydnvVpbU0VJchH7l/L4wZ1A3PvLNJKp13giW2QVRmGskXjb43jVqZpC8Tmc3rALo2TUQzwoqimi+LfGPEaI7edJJLe01oHqIcMk3UT1d76cWEZkS00/7yrLF4etVGTD2ojn8xBaKzsMOuAknbOOeUQ/DfjTvx9SnRB8M0/X2n2oN61YgtziC1pbcNSeOEdzUf9SNaqPT4scIIBZ86jHUdLx/O+UJMhfu/cXhrufkUm3c2cOuYTy90Weg6EWfIAiqm3riWj5nSsOTzakNuZ9M2VOKMBDUjtHuq9DyzoIngzxdMxE2PvY2fnDEymMyfN+W2eMvF5VYXAS2CJgwWQmYabmX4PiOqNW0WXJVCjzPiFWLjmyJqq8px6pgB/AsZ0ARs5lhJeBaDenTCou8dq1UXUeQ9Rjh9Q/GRsSId8+D5boloRtKKFUaoqH0lh3LhpC2MU8EXVVQY4eVL+z1+SHecPKof9unN30RPpVPznEFNvSa671sfj/3WxGZ9FYQ4I7JLe2nV8GpGgrEijty/F/5zKX/fJ9qOxbw79/pHibTS7ia6Zk9YSydQEE87QTfTyD/rYL1Zq2lE8AuxJmLIJE8Uq09MYVKTBaiPFzpaMG+JPJ+RKP2losAKIxREnxf/OvYFtNdj0w49YYRXTibj4DdnH0q8NohK543Lc151zf7d5x2G9Vv3+DQCPgdWiXy9XygqDmSicUa8Xz17GZM9C9nW+sfFR2D91j0Y4Qknzi3DAfY0NkuWFEZ2sCe9S/TljR6hTHGOajK4mkYlWF7acMD2FUr6S90f9KytLgf3q8Vbn5C3HGH57qnczYtXfk4hladM7xhFqIB/NQ0/jzRhhREKPgmU0e10d7+lSes1VWXYaFggUUXXTENeL2/mjShT9B4/dnif0DFR0xwrHdlnhE2wfUW2/laNoinb7OOGdA/tySLC7kaGZiRWB1ZyJia6n65mpHNF2/CrYt4LEuckQ2t3tgNrRJURhPa4fnTqIfj365/gy+MGSuWnMob17yYey4e8UZ5/ZgriHX9E9kRKE6W7DihihDsa5zK+zwj5+K/PGitWPj3nQDnqA6fKF413GR8pdVweJTLwliOLlKKighVdTeNFXTOi3/K86jmOwxRGhMuJ8HqfZkTQbBR8hUSFEdqrV1PVJoyYWFVi0n9KsQKJBz1jh+8nP4hunSrwwy8eIqX9A8zEqpFFRjPSRPHbSquZxgojFHi2OdJ1xPOcB0/TjBzSvytuPUfMjCKCzjecStgY3otq8n0wNSioRmD1LvEkmml4PiNBB1ZG40wY2gOVZRlf+HIZ4hqH9hgQRqKMMxJnKKSeNWSHxS4eYaSU40N4SbeZpu1v0fGdhQnNnQ7kXXvbjtGE5XSKIlYYoSLuM8Ix03DSszqrzrsbylfjpVCRpDO+pb0kdaM5OlWQrY2y+y74qilRwSbPqgqSup335S3qwArkd25+bc7x3CW89LLIx01/Le1u0vcZERrsKX8X86CkM+EzIsoXR/fHORMH41cBbadXGDES9CwhM83nD84vEf76kcMSN9OwhFPqqibl0qK9oU4VYadmnmak3CNlN+XIptyJ+/TQrlsUWJ8RCqIqT75mhH2e/WWn3tn36+NfJaMaZwQgf9EcWNcFqzfskEoTFdUVWexsaJv85p99KDbvbBBaKeTFt4Zfov5edShpUuEGPQs5sNKvzWQcVGXUV16YUeWT89i3d2dccGQ+XgRLMyLeE0U0I5wcqA6swpXQpiybwU1fyi+5/s59y4vHayrbBEojcUZixDue3PbVcajfvhcDulVj5fpt1DRJm5F8mhGvfxjTtEPPT2g7EIVbnnv6SNy3ZC2+d8KB4fw4Qc+8H4G0CMHf/tz+6F1TiWMODPvNJUlpvQExItpZufZzT4c5dUz/0HmW2VlkpUhQyDjqgN548crPoVsnv2pY9uvv7vMPK/5NqsbhHOmaG4HVoLBSHYiZctKofph5xFCFnNh+LjSG9uyEAd2qMbxvF+IqItkQ7ZEKcp6sTWsELjp6X5w9cTCAvJmRhqj/kshgzw0oR3VgTV5Z3V7MNNmMUwzW6O27v/zKGHzrmLZd2TtX6i9fjgLVniDSh1TesbMmDMYjlxyJPl3C2zN4zYu84vtRtneoKs/ivMnDMLRXevalAaxmRBuZMe2HpxyCV9Zswcfb9haPsQbTow7ojZEDuuKQ/uI7cHarLid6bMu+E97tx0mT456mFvSqqcDmnY3E9KrRTFUgqTN1kREIyrIZPPO9Y5SFiFDQs3hkkUj52ZdH4dZn3sPdL3ygnIdsO0jtJ+T5O6mNF2oNm2lMx9FQwfvMjti3F04dMwCHDesB13WZ2w2YQjSgmGxEYRKJOLAK1PuRSyZjy65GDBLYOy1NWGGEgOOIL/OUXdqbDURaHNKDLp1WlGXwr28fycxfFNYmaCTKMuxOn18tQb93XpwOk++xTDRZUWQFC9YSY54mIOwzEt0oZyRMu0AWfWqrMOeUQ4rCSEU2g0bBiKUFxBxY2/4+fJ+eGDWwK17/yGMqSF4BQqXW4/djYmlvzvROcAyoy0YJPg3HxmQOOObA3pj9+QOk06maj5IwO4lo7EcN7BZPZQxTurrBCCnPZoQlZ5mxPXjtKaP743+/MkaucpJlFDjmgPyLeueM8UL5eL/WiZoRzmqJIb6InYRKGXyPo9CMEKusOJHzpojgPBSlmcabs/drXCYcuUrtKgX21VEpxzv/lmczeOQSv/BOy8P3LBPSKHi/XE1oRnTjnsggUlLcq2fuOX9CyDxNwz+p069jaa6TsPR5i0xCMxMlVjNCoCKb4QayKsCyjQcJTsr6sUQgvDbecRx857j9hfMt8xgnSZ1+T1ML82U87qA+uOKE4RgxgGxiMvlVMbC7eCAhUeJ80YNCTpRle23OXavLcfUXhsOBE7kKvbIsix3IOxmLTpkiwt/h+/TAfUvWauUhirfe//72kTj5189r5Te0Z5tW1IRWI0ZZhIpXrkt6KW/UJCKM+DT27at9rTBCoKLML4ywVKj79anBPy4+wrfHSZzUVpXjtTnHY/T1/2c0X685acyg7vjrqx/5zvOCWjmOg4tbndfWbdlNOG+gkq1c/vkD8PibG7Bvb3MOWUYnMenVNFFqRvx5f+OofSlXytONsdxYZMfhICLN8MXR/VGezWAkRSinakY8f6vM4bIBskh4l56XGdgoLw0+Iz7a11wZIhlhwKuxTqD4CLHCCIGKgJmGt2+EcLjsiDqPaswJFl6fkemHDUKL62LC0B6YdsuzAPJmmrS8C/26VuP5K44NrarRIc4XPVYHVgN50wQ11pb0Q3p2wkef7QFgNhy84zj4wsh+3OuenH0Upt78rFjBETJ9/CA88Oo6HH9wHRzHwW/PORSr6nfg0MGCYwiDWIURgaJKZbJkmmkiigOlimospFLACiMEysscYTMNj2CHTdvHCw3vBJnNODj38CG+87pBrUy/R6bNDCa/emQfeSmqt887YijRaflPF0zAXc+vwU1fGokj5j0llaeJdihksV+fLsTjgPg7KeK0yeP6Uw/BtBF1OHyfngCAL4zsJyRMiSDpHxwJJs00hQ3sThzRV7NWbFS1oEm8p36fkdIbJ1hYYYRAeTbjk+r1hBH9JWRS5RnKJ2g6CLKnsUX4ZSiVd0bUqc00QcfDKIuOO77GlP17Y8r+iqHrDZQfx+2WS8SWryrP4nPD6RokHdJmptFt+z9eMAH/eeMTnDp2gJkKeRBdLSmaR1z45pP4i48Uu5qGQN5HxAn81sdB6UzMvKW5uxvZDqw80tgOkTnfcSaJ5oAwEuUXjwnVuX4WYpOmquA0c1KbFo8a9MzAXQxuXQ1zzIFqwpZpYl3aK/AMdftxr5pKnDtpKGojdq5WNtMw8rzo6Hwk4mtPPlitUgK0N82IFUYIBB1Yg5oRqeW8huqUFr44Oh9F9sIp+winIe9Nk+6WMSqLcM4HnTuj3MQtqnbX2RWaBu0ZdOYs5f7yuEHcPLyo1v2Biw7HVScOx8++PFopvWlaYtSMiBSV5Fx5xQnDmculeWHVRWAJy1edeBBe/+HxmHZIdCamdiaLWDMNiXycEe/vgIOhRt4p06RK8/MzR+O8yUMxemA3PPDKOuV80vgimVDdkqA988qyDC7//AEY2D0YKTHC1TRGHFj10gs7sBKO/eXrE7kRicUcXz31EatOiH5dq3HR0eZWI+kSp2aEhldjkuSX+8XH7IsLpwzDfj/4D/dadpwRtXQAItHoeAXntH/QyWKFEQL796nxSb3lGgGJ4n4fo/YJqCjLSHv+k2qURmHEuyrJ5EBKU2lPP2wQvkmYzKJchRBV1rKbEopAegYjBnTlrh6TfXal/oFQIAWyiI+kX3FWVGQjPiOK6UzhtDO7hhVGPPz1okl4aPlHuPLEg/DZrrY9V0z5jADpnISjplTuuVdNJX591lhUlWe5PjNREmmcEcNZ33P+YVixbmtxgzyTqNaVFjJ75qQh+MPiD5Xro7PzdRzEaqYRuKZkfBoUq5nE/XnbvURaVxgrjHiYMKwHJgzL70a7dXebMFLwGXGc/FeU6legi/bzFQboTmzpfJVOGR3eWVkX6pJQyvWRakZMLJf11Hz80B6RbUVOUkOLVN+3SaPn+KVTDyAKI+3llYx3bxpyWd7DaZNFjjqgzdFYNJIpy58o6fsrGWFPkHam6DGHL+hZq8/Ivy45EieP6offzxTb3yWYT6kIIi9ffZzQdaKvguqk0lGJdKO8FOQnHg5eIXP4hbmO1M/i3JtGhLiXkbM4dHA36vYb6v0s2ftLUfMawWpGKJBW04wY0BW/OftQ5Xzi+AQz0T/rauV2+FWhnb1HTGSF0GjNNKXT8uRVWGrpWGlFV9Ok/WMiDbJI985iG9XFzZfHDaL6GpXOG9G+9/6xwogAOjtqBrtLO+s/QpDuuSO1g+wcEWk4eBN5xLSaRtVc1Z43E2MRZ9AzWkkDulXjp2eMQm0EW1TowIqErSqgd6QxLA6sMCKATgRWryCTzTrRf13F+IKIvsT2nZUjWp+R6PI2DVmI5d+Az2ekA5ls0hKB9f8dNoh/UYpgdQtWiyatmUi6fNNYYYRCk2ejBx3NSLdOFbjyxOHIOEBNZQzNnY7xiEtH+mKVDapVSqtpZPI7cURffPDpbhw6uJtY3op9xOczYiC/tHL6oQPw4LL1xd9x+oykRO5RRtTRdkR/+u7MScsCSZdvGiuMUPCG6A4GPZOFFEuiQ2HopSmsZmr3pHw1jSq3fnUcXNcVrgNJQyTkMwKyZqS9Mff0kfh/4wfhK7e/BCA9mpE0wuoGLCH14P61uPfCibj4z8uwbU+TcLroSEdQuSiwq2koeDUjMhthJU6J9E+V96hEbk2bSDUjkeUsWL7UXgqGbfmU46U6h1eWZYu7/wJALsZde0u0yYr4YsZwutkR+/bC0J7BKMnRmlNFSPpdNk0JzbLx0qXSE40z6V6XUsoENUbeL4jCviJRxPOIGtVeUFVO3kuFNjGnPc5IXBA1IyJxRnwJzd1v2ifgOIOelQr79u4MADh2OD0WjvIrkcCr1LumbaVjCb3KQlgzDYXBPTvhR6cegm6d0rlULQ38+qyxOP/uV3DFCcOZ13lfmgWXHYXmnIthvTpLl+ckbKdRLfnnZ47GN/+8FN/+3H7+/Cj3knbNiG9flwgfR1ZVM+L9W8CBNe2RVUWJdW+alAg+Rx3QG8++swknjexHPL/gsqOwu7EltKzXFVeMUEnCTNO1Uzn+OWsyKsszJfVhIYIVRhjMmDQ06SqkmlEDu+HVa6ZKvRRV5Vn07lKpVF6pvnr79anBk7OPFr4+0qBnBrL21i/KKSlLUI2ITAC0CKxeUjKXGqUj+ozMP3ssnlq1EVMPqiOeL89m0LWabQBQfd9ENZgZx2wMmNGDupnLLEVYYaSdEbe0LvIi075USw3j0UsTMNOUktObqBkwiD8Cq0HBJOWTfUucipH4imLSpaocp44ZoJWHsmZENLRBh/G818P6jFhiRWcqTHoeNT2c0Mw0UQqUpSOKAFmC47hQH6BcQ5s8hhCcE0uJutq8pnHaIWTtgCWM9B46hItEPxpOH5sXlg4b2l0sQQfFakbaGUlP2BZ9IvWXNmGm8fwtG0NFhjLFhhAx0wDAfRcejtc+2opph/RVKictPH7ZUXhnw85YJ7v29KGvKvyLpvvRqSMw5YDeOHr/3vyLOzBWGGlnHFjXJekqhDA1buVf/vhHwa7V5di2pwlH7tfLaL5UM02E0kgpBf4i+YyIQDNFBY9O2rcnJu3bk3htKdGtU0Vxt3GLGF6n5ahX01RXZPHFElw9GDdWGGknPHLJZDz37mbMPGJo0lVhouWcmdA8+uh3jsQTb23AdMNhrqPUKtAoJc0ZKdigSP19mhtz1bFYfNiID2axwkg7YdTAbhg1sFvS1SBiYhldkgzs3gnnTx4WW3lROpmayNu3tFc7NzoknxERvPfoFfhKSRBLO6W+HFrWZ4R0SXtbWps01oHVEjlGVKIoTUGGRTJBz6LL2zQknxERM5PjGdVKe8pML9ZnpP2NR0ljhRFLyVBKE6kO+/WpiSxvE03o3fCxQmNHax6qPiN+B1vvcW98lHY0m1qk8T591XGllJbJlwLWTGOJHjvuS/G54X3wo1MPwSH9a43nbWL87FJVjrvOG4+M41BD3ZuAqBkRCQdPc2DVvHfbjdsn6nFGjFajw2OFEUuslNJqjqRwHCfC6L9m2v9zw6OPaWFmNY0VIaKg1M00fl8iOyalAWumiZljDsyvNe+jGBK9FDG7tNeiQymNuyTNiMgkGNfeOR2ZUm9Wn5lGMQ9rpjGL1YzEzNVfOAgH1nXB1IM7aLREHQdW++5rU0oDKGk1jYiyJKrVPlawaZ+ovhIl9CqVBFYYiZlOFWU4t4NtwGdqELfvvj5JxDZRxbs3zbmHD8EJI/qiTMBh1r+0N5KqWUoc/9Jeu5omDVgzjQUA8JuzxwIAfvblUcbz7t2lEv26VmFAt2p0qbTyb5K0xLnNvCZen5Fjh/fGZMEIuH6PkdK531KilITaqCglLWMpYGcGCwDg5FH9cfzBfVFRZl4+zWYcPPf9YwHohTq3jmb6NMW5tasm5Z6+IuMvFNVqGkt7Qu49IPUd25/MYoURS5EoBJECIup1Hvbd16e0NCOePiPx8K0DqyUO7MeRWayZxmJp53hXpTTncgnWRA6vz4jMsO+dJHJWGrFYSgIrjFgs7RzvB1wpTc5en5E0fIVaP4k2bFNYTGOFEUvpkPx8VPKUks8IKc6IXn5tw13X6nKjeVtKC1lhyg490WN9Riwlgx0Q1Mg7f+ZH39LyGVEz03jxTjrZjIMHvnE4Glty6NapQq9ylpJG9i0onbemdLHCiMXS3vHM5MP7dkmuHpKU+cw0ankEl/ZO3Kencn3shNSGXTJtMY0VRiwlQxr8BkoRB8DzVxyLzTsbsU/v6HYENo13NY3qVgDWtyEaSr1dS73+7RHrM2IpGdqLLHLu4UNQWZbB16cMi63Mgd07YcygbrGVZwIjmhE76VgI7N9HTihvJ0NPqrGaEUvJ0F4GhBtOG4E5pxxsJPaKCKUqxJnwGbFEw42njcD021/C/3z+gKSrokT3zhV4/opj0anCToFpwT4JiyUB4hJEgNLd7dgbZ0T1FqxiJBom7tMT79x4YqSBEqNmYPdOSukev+wodKmyU6dpbItaLO2cUtWMlBF27ZXFxgaJjlIWRHQ4sIScwEuJjtmbLCWJdWBVo1RbzTvZKTuwmqoMrP+JxRIlVhixlAylOqkmTakKcVXlHmFEOdCImbpYOjal+g6VElYYsVgsqaS6PFv8W10WsdKIxVIKKAkj8+fPx9ChQ1FVVYWJEydiyZIl1GvvuOMOTJkyBd27d0f37t0xdepU5vUWCw37caJGqTZblVcYScHDT0EVLAlRax1WI0daGHnggQcwe/ZszJkzB8uWLcPo0aMxbdo0bNy4kXj9okWLcNZZZ+Hpp5/G4sWLMWjQIBx//PFYv369duUtFosAJTqJ+oURtTxM+nmcMW4gBnSrxlcPH2wuU0tJcMNpIzBmUDf85uyxSVel3SIt7t1888248MILcf755wMAbrvtNjz66KO46667cOWVV4au/8tf/uL7feedd+If//gHFi5ciBkzZihW29IxKdFZNWFKtdV8PiOKeZjciqe2qhzPX3FsKrQ0lngZ2L0THp41OelqtGukNCONjY1YunQppk6d2pZBJoOpU6di8eLFQnns3r0bTU1N6NGjB/WahoYGbN++3fefxXLUAb0AAL1q7CZnHQGvz0ha9vezgojFEg1SmpHNmzejpaUFdXV1vuN1dXVYtWqVUB5XXHEF+vfv7xNogsydOxfXX3+9TNUsHYDrv3gIDu5XixNH9ku6KiVFqU6g1RVtwkhjc04pj9pqa+u3WEqBWFfTzJs3D/fffz8eeughVFVVUa+76qqrsG3btuJ/69ati7GWlrTSpaocX5+yDwZ0q066KiVBz855DdIkjZ1qk6SqrE0Y2dvUIpV23ukj8eVxA3HiCCu4WiylgNRnQ69evZDNZrFhwwbf8Q0bNqBv377MtD//+c8xb948PPnkkxg1ahTz2srKSlRWVspUzWKxBHh41mQ8vHw9ZkwamnRVlMh49qZpkNSMfGXCYHxlgnU0tVhKBSnNSEVFBcaNG4eFCxcWj+VyOSxcuBCTJk2ipvvpT3+KG264AQsWLMD48ePVa2uxWIQZ1KMTvn3c/ujaqTzpqmgjqxmxWCylhbRBdfbs2Zg5cybGjx+PCRMm4JZbbsGuXbuKq2tmzJiBAQMGYO7cuQCAn/zkJ7juuutw7733YujQoaivrwcA1NTUoKZGbhtni8XSMdnbbIURi6U9Iy2MTJ8+HZs2bcJ1112H+vp6jBkzBgsWLCg6ta5duxYZzwZXt956KxobG/HlL3/Zl8+cOXPwwx/+UK/2FoulXXPyqH5Y/N6nOHlU/6SrYrFYIsRxS2Bby+3bt6Nr167Ytm0bamtrk66OxWKJCdd10ZJzUZa1O1dYLKWI6Pxt171ZLJbU4jgOyrKluTTZYrGIYz83LBaLxWKxJIoVRiwWi8VisSSKFUYsFovFYrEkihVGLBaLxWKxJIoVRiwWi8VisSSKFUYsFovFYrEkihVGLBaLxWKxJIoVRiwWi8VisSSKFUYsFovFYrEkihVGLBaLxWKxJIoVRiwWi8VisSSKFUYsFovFYrEkihVGLBaLxWKxJIoVRiwWi8VisSSKFUYsFovFYrEkihVGLBaLxWKxJIoVRiwWi8VisSSKFUYsFovFYrEkihVGLBaLxWKxJIoVRiwWi8VisSSKFUYsFovFYrEkihVGLBaLxWKxJIoVRiwWi8VisSSKFUYsFovFYrEkihVGLBaLxWKxJIoVRiwWi8VisSSKFUYsFovFYrEkihVGLBaLxWKxJIoVRiwWi8VisSSKFUYsFovFYrEkihVGLBaLxWKxJIoVRiwWi8VisSSKFUYsFovFYrEkihVGLBaLxWKxJIoVRiwWi8VisSSKFUYsFovFYrEkihVGLBaLxWKxJIoVRiwWi8VisSSKFUYsFovFYrEkihVGLBaLxWKxJIoVRiwWi8VisSSKFUYsFovFYrEkihVGLBaLxWKxJIoVRiwWi8VisSSKFUYsFovFYrEkihVGLBaLxWKxJIoVRiwWi8VisSSKFUYsFovFYrEkihVGLBaLxWKxJIoVRiwWi8VisSSKFUYsFovFYrEkihVGLBaLxWKxJIoVRiwWi8VisSSKFUYsFovFYrEkihVGLBaLxWKxJIoVRiwWi8VisSSKFUYsFovFYrEkihVGLBaLxWKxJIoVRiwWi8VisSSKFUYsFovFYrEkihVGLBaLxWKxJIoVRiwWi8VisSSKkjAyf/58DB06FFVVVZg4cSKWLFnCvP5vf/sbhg8fjqqqKowcORKPPfaYUmUtFovFYrG0P6SFkQceeACzZ8/GnDlzsGzZMowePRrTpk3Dxo0bide/+OKLOOuss3DBBRdg+fLlOO2003Daaadh5cqV2pW3WCwWi8VS+jiu67oyCSZOnIjDDjsMv/nNbwAAuVwOgwYNwre//W1ceeWVoeunT5+OXbt24d///nfx2OGHH44xY8bgtttuEypz+/bt6Nq1K7Zt24ba2lqZ6losFovFYkkI0fm7TCbTxsZGLF26FFdddVXxWCaTwdSpU7F48WJimsWLF2P27Nm+Y9OmTcPDDz9MLaehoQENDQ3F39u2bQOQvymLxWKxWCylQWHe5uk9pISRzZs3o6WlBXV1db7jdXV1WLVqFTFNfX098fr6+npqOXPnzsX1118fOj5o0CCZ6losFovFYkkBO3bsQNeuXannpYSRuLjqqqt82pRcLoctW7agZ8+ecBzHWDnbt2/HoEGDsG7dOmv+iRjb1vFg2zkebDvHg23n+IiqrV3XxY4dO9C/f3/mdVLCSK9evZDNZrFhwwbf8Q0bNqBv377ENH379pW6HgAqKytRWVnpO9atWzeZqkpRW1trO3pM2LaOB9vO8WDbOR5sO8dHFG3N0ogUkFpNU1FRgXHjxmHhwoXFY7lcDgsXLsSkSZOIaSZNmuS7HgCeeOIJ6vUWi8VisVg6FtJmmtmzZ2PmzJkYP348JkyYgFtuuQW7du3C+eefDwCYMWMGBgwYgLlz5wIALr30Uhx99NH4xS9+gZNOOgn3338/Xn31Vdx+++1m78RisVgsFktJIi2MTJ8+HZs2bcJ1112H+vp6jBkzBgsWLCg6qa5duxaZTJvC5YgjjsC9996La665BldffTX2339/PPzwwxgxYoS5u1CksrISc+bMCZmELOaxbR0Ptp3jwbZzPNh2jo+k21o6zojFYrFYLBaLSezeNBaLxWKxWBLFCiMWi8VisVgSxQojFovFYrFYEsUKIxaLxWKxWBKlQwsj8+fPx9ChQ1FVVYWJEydiyZIlSVepZJg7dy4OO+wwdOnSBX369MFpp52G1atX+67Zu3cvZs2ahZ49e6KmpgZnnHFGKADe2rVrcdJJJ6FTp07o06cPvve976G5uTnOWykp5s2bB8dxcNlllxWP2XY2x/r16/HVr34VPXv2RHV1NUaOHIlXX321eN51XVx33XXo168fqqurMXXqVLz77ru+PLZs2YJzzjkHtbW16NatGy644ALs3Lkz7ltJLS0tLbj22msxbNgwVFdXY99998UNN9zg27vEtrMazz77LE455RT0798fjuOE9oAz1a6vv/46pkyZgqqqKgwaNAg//elP9SvvdlDuv/9+t6Kiwr3rrrvcN998073wwgvdbt26uRs2bEi6aiXBtGnT3LvvvttduXKlu2LFCvcLX/iCO3jwYHfnzp3Fa775zW+6gwYNchcuXOi++uqr7uGHH+4eccQRxfPNzc3uiBEj3KlTp7rLly93H3vsMbdXr17uVVddlcQtpZ4lS5a4Q4cOdUeNGuVeeumlxeO2nc2wZcsWd8iQIe55553nvvzyy+7777/vPv744+5///vf4jXz5s1zu3bt6j788MPua6+95n7xi190hw0b5u7Zs6d4zQknnOCOHj3afemll9znnnvO3W+//dyzzjoriVtKJTfddJPbs2dP99///re7Zs0a929/+5tbU1Pj/vKXvyxeY9tZjccee8z9wQ9+4D744IMuAPehhx7ynTfRrtu2bXPr6urcc845x125cqV73333udXV1e7vfvc7rbp3WGFkwoQJ7qxZs4q/W1pa3P79+7tz585NsFaly8aNG10A7jPPPOO6rutu3brVLS8vd//2t78Vr3n77bddAO7ixYtd182/OJlMxq2vry9ec+utt7q1tbVuQ0NDvDeQcnbs2OHuv//+7hNPPOEeffTRRWHEtrM5rrjiCvfII4+kns/lcm7fvn3dn/3sZ8VjW7dudSsrK9377rvPdV3Xfeutt1wA7iuvvFK85j//+Y/rOI67fv366CpfQpx00knu1772Nd+x008/3T3nnHNc17XtbIqgMGKqXX/729+63bt3940dV1xxhXvggQdq1bdDmmkaGxuxdOlSTJ06tXgsk8lg6tSpWLx4cYI1K122bdsGAOjRowcAYOnSpWhqavK18fDhwzF48OBiGy9evBgjR4707eo8bdo0bN++HW+++WaMtU8/s2bNwkknneRrT8C2s0keeeQRjB8/HmeeeSb69OmDsWPH4o477iieX7NmDerr631t3bVrV0ycONHX1t26dcP48eOL10ydOhWZTAYvv/xyfDeTYo444ggsXLgQ77zzDgDgtddew/PPP48TTzwRgG3nqDDVrosXL8ZRRx2FioqK4jXTpk3D6tWr8dlnnynXL5W79kbN5s2b0dLS4hucAaCurg6rVq1KqFalSy6Xw2WXXYbJkycXI+vW19ejoqIitMFhXV0d6uvri9eQnkHhnCXP/fffj2XLluGVV14JnbPtbI73338ft956K2bPno2rr74ar7zyCr7zne+goqICM2fOLLYVqS29bd2nTx/f+bKyMvTo0cO2dStXXnkltm/fjuHDhyObzaKlpQU33XQTzjnnHACw7RwRptq1vr4ew4YNC+VRONe9e3el+nVIYcRillmzZmHlypV4/vnnk65Ku2PdunW49NJL8cQTT6Cqqirp6rRrcrkcxo8fjx//+McAgLFjx2LlypW47bbbMHPmzIRr137461//ir/85S+49957ccghh2DFihW47LLL0L9/f9vOHZgOaabp1asXstlsaMXBhg0b0Ldv34RqVZpccskl+Pe//42nn34aAwcOLB7v27cvGhsbsXXrVt/13jbu27cv8RkUzlnyZpiNGzfi0EMPRVlZGcrKyvDMM8/gV7/6FcrKylBXV2fb2RD9+vXDwQcf7Dt20EEHYe3atQDa2oo1bvTt2xcbN270nW9ubsaWLVtsW7fyve99D1deeSW+8pWvYOTIkTj33HNx+eWXFzdXte0cDabaNarxpEMKIxUVFRg3bhwWLlxYPJbL5bBw4UJMmjQpwZqVDq7r4pJLLsFDDz2Ep556KqS2GzduHMrLy31tvHr1aqxdu7bYxpMmTcIbb7zh6/xPPPEEamtrQ5NCR+W4447DG2+8gRUrVhT/Gz9+PM4555zi37adzTB58uTQ8vR33nkHQ4YMAQAMGzYMffv29bX19u3b8fLLL/vaeuvWrVi6dGnxmqeeegq5XA4TJ06M4S7Sz+7du32bqQJANptFLpcDYNs5Kky166RJk/Dss8+iqampeM0TTzyBAw88UNlEA6BjL+2trKx077nnHvett95yv/GNb7jdunXzrTiw0Ln44ovdrl27uosWLXI/+eST4n+7d+8uXvPNb37THTx4sPvUU0+5r776qjtp0iR30qRJxfOFJafHH3+8u2LFCnfBggVu79697ZJTDt7VNK5r29kUS5YsccvKytybbrrJfffdd92//OUvbqdOndw///nPxWvmzZvnduvWzf3nP//pvv766+6pp55KXBo5duxY9+WXX3aff/55d//99+/wS069zJw50x0wYEBxae+DDz7o9urVy/3+979fvMa2sxo7duxwly9f7i5fvtwF4N58883u8uXL3Q8//NB1XTPtunXrVreurs4999xz3ZUrV7r333+/26lTJ7u0V4df//rX7uDBg92Kigp3woQJ7ksvvZR0lUoGAMT/7r777uI1e/bscb/1rW+53bt3dzt16uR+6Utfcj/55BNfPh988IF74oknutXV1W6vXr3c//mf/3GbmppivpvSIiiM2HY2x7/+9S93xIgRbmVlpTt8+HD39ttv953P5XLutdde69bV1bmVlZXucccd565evdp3zaeffuqeddZZbk1NjVtbW+uef/757o4dO+K8jVSzfft299JLL3UHDx7sVlVVufvss4/7gx/8wLdU1LazGk8//TRxXJ45c6bruuba9bXXXnOPPPJIt7Ky0h0wYIA7b9487bo7rusJe2exWCwWi8USMx3SZ8RisVgsFkt6sMKIxWKxWCyWRLHCiMVisVgslkSxwojFYrFYLJZEscKIxWKxWCyWRLHCiMVisVgslkSxwojFYrFYLJZEscKIxWKxWCyWRLHCiMVisVgslkSxwojFYrFYLJZEscKIxWKxWCyWRLHCiMVisVgslkT5/2cPJIT3rJ2MAAAAAElFTkSuQmCC", + "text/plain": [ + "
" + ] + }, + "metadata": {}, + "output_type": "display_data" } ], "source": [ - "true_parameters" + "params = sample_parameters(true_parameters)\n", + "\n", + "errors = compute_error(params)\n", + "plt.plot(errors)\n", + "plt.ylim(0, 1e7)\n" ] }, { "cell_type": "code", - "execution_count": 6, + "execution_count": 12, "metadata": {}, "outputs": [ { "data": { "text/plain": [ - "Array([ 0.47510388, -0.28072184, 0.3476282 , 0.26952767, 0.85689604,\n", - " 0.33116183, 1.5397867 , 0.31192046, -0.16261268, 1.8731921 ], dtype=float32)" + "Array([ 9.20532346e-01, 1.13579944e-01, -5.36158355e-03, 5.26336469e-02,\n", + " 3.09643317e-02, 7.71308027e-04, 4.48253453e-02, -9.11160186e-03,\n", + " 4.05511179e-04, 6.89987987e-02], dtype=float32)" ] }, - "execution_count": 6, + "execution_count": 12, "metadata": {}, "output_type": "execute_result" } ], "source": [ - "# sample random parameters in log-cholesky space\n", - "key = jax.random.PRNGKey(0)\n", + "# find the best parameters\n", + "idx = jnp.argmin(errors)\n", + "best_parameters = params[idx]\n", "\n", - "# sample 10 random parameters\n", - "random_parameters = jax.random.normal(key, (10,))\n", - "theta_random = logchol2theta(random_parameters)\n", + "logchol2theta(best_parameters)" + ] + }, + { + "cell_type": "code", + "execution_count": 13, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "Array([ 1.3327169e+00, 6.5090909e-04, -3.8447324e-03, 7.5462282e-02,\n", + " 4.1322153e-02, -6.5983908e-04, 3.0134678e-02, -2.2140213e-03,\n", + " 1.5658194e-04, 6.1365250e-02], dtype=float32)" + ] + }, + "execution_count": 13, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "lambda_ = 1e-10\n", + "weights = jnp.zeros(N_SAMPLES)\n", + "for i in range(N_SAMPLES):\n", + " error = errors[i]\n", + " weights = weights.at[i].set(jnp.exp(-error * lambda_))\n", + "\n", + "weights /= jnp.sum(weights)\n", + "# use the weights to compute the new parameters\n", + "new_parameters = jnp.sum(weights[:, None] * params, axis=0)\n", "\n", - "theta_random" + "logchol2theta(new_parameters)" ] }, { "cell_type": "code", - "execution_count": 7, + "execution_count": 14, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ - "pos [-8.60458879e-04 8.27660234e-07 9.92544927e-02 9.99967927e-01\n", - " 7.70652416e-06 8.00899567e-03 7.87957120e-07]\n", - "next_pos [-5.1760540e-04 1.3108459e-05 1.0056238e-01 9.9998802e-01\n", - " 1.2252582e-04 4.8988331e-03 -7.4357864e-05]\n", - "diff l2 0.0033942016\n" + "Step: 0, Error: 58176900.0 Theta: [ 1.3327945e+00 1.4981302e-03 -4.2728218e-03 7.5407363e-02\n", + " 4.1317686e-02 -7.3245692e-04 3.0044770e-02 -2.2448506e-03\n", + " 2.1679589e-04 6.1285086e-02]\n", + "Step: 1, Error: 26420402.0 Theta: [ 1.3334370e+00 6.4618990e-04 -3.7520514e-03 7.5623877e-02\n", + " 4.1360307e-02 -7.0125330e-04 3.0167455e-02 -2.2152588e-03\n", + " 1.4892367e-04 6.1402701e-02]\n", + "Step: 2, Error: 26410714.0 Theta: [ 1.3328096e+00 8.2522689e-04 -3.8017740e-03 7.5434141e-02\n", + " 4.1328318e-02 -6.8677368e-04 3.0142022e-02 -2.2241846e-03\n", + " 1.5410635e-04 6.1385911e-02]\n", + "Step: 3, Error: 26408272.0 Theta: [ 1.3328104e+00 8.2509348e-04 -3.8020690e-03 7.5434171e-02\n", + " 4.1328348e-02 -6.8674842e-04 3.0142015e-02 -2.2241736e-03\n", + " 1.5412638e-04 6.1385933e-02]\n", + "Step: 4, Error: 16411237.0 Theta: [ 1.3322150e+00 7.1831804e-04 -3.8392865e-03 7.5388931e-02\n", + " 4.1304972e-02 -7.3728012e-04 3.0114530e-02 -2.2147393e-03\n", + " 1.5573596e-04 6.1342273e-02]\n" ] } ], "source": [ - "next_pos, next_vel = smart_step(acc, vel, pos, ctrl, theta_random)\n", + "def optimize(theta0):\n", + " ITERATIONS = 5\n", + "\n", + " theta = theta0\n", + " for idx in range(ITERATIONS):\n", + " logchols = sample_parameters(theta)\n", + " errors = compute_error(logchols)\n", + "\n", + " weights = jnp.zeros(N_SAMPLES)\n", + " lambda_ = 1e-9\n", + " for i in range(N_SAMPLES):\n", + " error = errors[i]\n", + " if jnp.isinf(error):\n", + " weights = weights.at[i].set(0)\n", + " else:\n", + " weights = weights.at[i].set(jnp.exp(-error * lambda_))\n", + "\n", + " weights /= jnp.sum(weights)\n", + " new_parameters = jnp.sum(weights[:, None] * params, axis=0)\n", + " theta = logchol2theta(new_parameters)\n", "\n", - "print(\"pos\", log.data(\"qpos\")[1])\n", - "print(\"next_pos\", next_pos)\n", - "print(\"diff l2\", jnp.linalg.norm(log.data(\"qpos\")[1] - next_pos))" + " print(f\"Step: {idx}, Error: {errors.mean()} Theta: {theta}\")\n", + "\n", + " return theta\n", + "\n", + "\n", + "true_logchol = theta2logchol(true_parameters)\n", + "true_logchol += jax.random.normal(key, shape=true_logchol.shape) * 1.0\n", + "test_theta = logchol2theta(true_logchol)\n", + "\n", + "last_theta = optimize(test_theta)" + ] + }, + { + "cell_type": "code", + "execution_count": 15, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "Array([ 1.325 , 0. , 0. , 0.0715 , 0.04051 , 0. ,\n", + " 0.02927 , -0.0021 , 0. , 0.060528], dtype=float32)" + ] + }, + "execution_count": 15, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "true_parameters\n" + ] + }, + { + "cell_type": "code", + "execution_count": 17, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "Array(0.006688, dtype=float32)" + ] + }, + "execution_count": 17, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "jnp.linalg.norm(theta2logchol(last_theta) - theta2logchol(true_parameters))" ] } ],