diff --git a/examples/inverse_solver_2d.ipynb b/examples/inverse_solver_2d.ipynb new file mode 100644 index 0000000..589c078 --- /dev/null +++ b/examples/inverse_solver_2d.ipynb @@ -0,0 +1,274 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "#### Inverse solver example for body under uniform gravity loading\n", + "This notebook demonstrates an example for inverse solver for a body under uniform gravity loading and tries to estimate the youngs modulus of the body needed to get the observed stress. \n", + "\n", + "First we solved the problem of calculating the stress on the body under the given conditions and given youngs modulus of 1000. \n", + "\n", + "The following code runs the MPM solver for a body represented by 4 material points and a unit cell mesh. The body is modelled as linear elastic material with density=1, youngs modulus=1000. The particles are located at [(0,0),(0.5,0),(0,0.5),(0.5,0.5)].\n" + ] + }, + { + "cell_type": "code", + "execution_count": 1, + "metadata": {}, + "outputs": [ + { + "name": "stderr", + "output_type": "stream", + "text": [ + "No GPU/TPU found, falling back to CPU. (Set TF_CPP_MIN_LOG_LEVEL=0 and rerun for more info.)\n" + ] + } + ], + "source": [ + "import sys,os\n", + "sys.path.append(os.path.abspath('../'))\n", + "from diffmpm.material import SimpleMaterial, LinearElastic\n", + "from diffmpm.particle import Particles\n", + "from diffmpm.element import Quadrilateral4Node\n", + "from diffmpm.constraint import Constraint\n", + "from diffmpm.mesh import Mesh2D\n", + "from diffmpm.solver import MPMExplicit\n", + "import jax.numpy as jnp\n", + "import numpy as np\n", + "from jax import jit, value_and_grad\n", + "import optax\n", + "from tqdm import tqdm\n", + "import time\n", + "from bayes_opt import BayesianOptimization\n", + "import matplotlib.pyplot as plt\n", + "mesh_config = {}\n", + "density = 1\n", + "poisson_ratio = 0.0\n", + "youngs_modulus = 1000\n", + "# Linear Elastic Material with given parameters\n", + "material = LinearElastic(\n", + " {\n", + " \"id\": 0,\n", + " \"youngs_modulus\": youngs_modulus,\n", + " \"density\": density,\n", + " \"poisson_ratio\": poisson_ratio,\n", + " }\n", + ")\n", + "# location of the material points\n", + "particle_loc = jnp.array([[0.0, 0.0], [0.5, 0.0], [0.0, 0.5], [0.5, 0.5]]).reshape(\n", + " 4, 1, 2\n", + ")\n", + "# particle object for material points\n", + "particles = Particles(\n", + " particle_loc, material, jnp.zeros(particle_loc.shape[0], dtype=jnp.int32)\n", + ")\n", + "particles.velocity = particles.velocity.at[:].set(0.0)\n", + "# velocity constraint for material point at (0,0)\n", + "constraints = [(0, Constraint(1, 0.0))]\n", + "# gravity loading\n", + "gravity_loading = jnp.array([0.0, -9.8]).reshape(1, 2)\n", + "element = Quadrilateral4Node([1, 1], 1, [1, 1], constraints)\n", + "mesh_config[\"particles\"] = [particles]\n", + "mesh_config[\"elements\"] = element\n", + "mesh_config[\"particle_surface_traction\"] = []\n", + "mesh = Mesh2D(mesh_config)\n", + "solver = MPMExplicit(mesh, 0.01, sim_steps=10)\n", + "# getting the target value for youngs modulus 1000\n", + "target_ans = solver.solve_jit(gravity_loading)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "After getting the target stress. The following code tries to estimate the youngs_modulus for the body. The following code uses Adam optimizer for estimating the youngs modulus of the body by calculating error between target stress and estimated stress. The estimated stress is calculated by running the MPM solver for the estimated youngs modulus. \n", + "\n", + "The compute loss function generates a material object at each iteration for the youngs modulus at the given iteration and solves the generated configuration. \n", + "\n", + "The optax_adam function calculates gradient of the stress with respect to the youngs_modulus and updates the parameter. The inital guess for the youngs modulus is 1050" + ] + }, + { + "cell_type": "code", + "execution_count": 2, + "metadata": {}, + "outputs": [ + { + "name": "stderr", + "output_type": "stream", + "text": [ + "Youngs_modulus: 1000.3665161132812: 100%|██████████| 500/500 [00:51<00:00, 9.65it/s]\n" + ] + }, + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Time taken by adam: 51.79610347747803\n", + "parameter = 1000.3665161132812 and loss = 0.004986115265637636\n" + ] + }, + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAjcAAAHFCAYAAAAOmtghAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/bCgiHAAAACXBIWXMAAA9hAAAPYQGoP6dpAABeiklEQVR4nO3dd1RU19oG8GfoHaVbEMEGilRjjcYSe+/tqslNkasGa2KLscRri7GLJbElsSsaYw2JGjUSjVJEwY6CCCIovc/s7w+u82XEwujAgZnnt9as5ezZs8979uXCk7NPkQkhBIiIiIi0hJ7UBRARERFpEsMNERERaRWGGyIiItIqDDdERESkVRhuiIiISKsw3BAREZFWYbghIiIircJwQ0RERFqF4YaIiIi0CsMNkZZatWoVZDIZPD091fre1q1bIZPJcO/evbIprAKIjo7GnDlzXriPH3zwAWrXrq3R7ZXFmET0cgw3RFpq8+bNAIBr167hwoULEldTsURHR2Pu3LkvDDezZs3CgQMHNLq9shiTiF6O4YZIC126dAmRkZHo3r07AGDTpk0SV1R51KlTB76+vhV+zNfJzc0t1+0RVSQMN0Ra6FmYWbRoEVq2bIldu3YhJyenRL+//voLrVq1gomJCapXr47p06ejsLCwRL/du3ejU6dOqFatGkxNTeHh4YFp06YhOztbpd8HH3wACwsLXL9+HZ07d4a5uTmqVauGRYsWKbf37rvvwtzcHPXr18e2bdtKtT9PnjzBmDFjUKNGDRgZGcHNzQ0zZ85Efn6+Sj+ZTIZx48Zhw4YNqF+/PoyNjdGwYUPs2rVL2Wfr1q0YOHAgAKBdu3aQyWSQyWTYunWrch+eX0J6Nu6WLVvQoEEDmJqaokmTJvjrr78ghMA333wDV1dXWFhYoH379rh9+3aJefnnmHPmzFFu9/nXBx98oOxXUFCA+fPnw93dHcbGxrC3t8eHH36Ix48fq4xfu3Zt9OjRA8HBwfD19YWJiQnmzp1bqrkl0kqCiLRKTk6OsLa2Fu+8844QQojvv/9eABBbt25V6Xft2jVhZmYmGjZsKHbu3Cl+/vln0blzZ1GrVi0BQMTGxir7fv3112L58uXiyJEj4vTp02L9+vXC1dVVtGvXTmXMUaNGCSMjI+Hh4SFWrlwpQkJCxIcffigAiOnTp4v69euLTZs2iRMnTogePXoIAOLSpUuv3J/c3Fzh5eUlzM3NxdKlS8Wvv/4qZs2aJQwMDES3bt1U+gIQzs7Oyn06dOiQ6NKliwAg9u7dK4QQIjk5WSxYsEAAEGvXrhWhoaEiNDRUJCcnK/fBxcWlxLguLi6iZcuWIjg4WBw4cEDUr19f2NjYiIkTJ4revXuLw4cPi+3btwtHR0fh5eUlFAqFyrz8c8z4+Hjldp+9Pv/8cwFALFmyRAghhFwuF126dBHm5uZi7ty5IiQkRHz//feiRo0aomHDhiInJ0c5nouLi6hWrZpwc3MTmzdvFqdOnRIXL1585bwSaTOGGyIt88MPPwgAYv369UIIITIzM4WFhYVo3bq1Sr/BgwcLU1NTkZSUpGwrKioS7u7uJcLNPykUClFYWCj++OMPAUBERkYqPxs1apQAIPbv369sKywsFPb29gKACAsLU7anpqYKfX19MWnSpFfuz/r16wUAsWfPHpX2xYsXCwDi119/VbYBeOk+1a1bV9m2d+9eAUCcOnWqxPZeFm6cnJxEVlaWsu3gwYMCgPDx8VEJMitWrBAAxJUrV1455j+dPXtWmJiYiOHDhyvH2rlzZ4m5FEKIv//+WwAQQUFByjYXFxehr68vbty48dJtEOkSLksRaZlNmzbB1NQUQ4YMAQBYWFhg4MCBOHv2LG7duqXsd+rUKXTo0AGOjo7KNn19fQwePLjEmHfv3sWwYcPg5OQEfX19GBoa4r333gMAxMTEqPSVyWTo1q2b8r2BgQHq1q2LatWqqZx3YmNjAwcHB9y/f/+V+3Py5EmYm5tjwIABKu3Plm9+//13lfaX7dPt27fx4MGDV27rVdq1awdzc3Plew8PDwBA165dIZPJSrS/br+eiYmJQa9evdCyZUts3rxZOdbhw4dRpUoV9OzZE0VFRcqXj48PnJyccPr0aZVxvLy8UL9+/TfePyJtwnBDpEVu376NM2fOoHv37hBCIC0tDWlpacpg8OwKKgBITU2Fk5NTiTGeb8vKykLr1q1x4cIFzJ8/H6dPn8bff/+N4OBgACVPXDUzM4OJiYlKm5GREWxsbEpsy8jICHl5ea/cp2d1/jNAAICDgwMMDAyQmpr6yvr/2fZ8X3U8X7+RkdEr21+3XwDw8OFDdOnSBTVr1kRwcLDyuwDw6NEjpKWlwcjICIaGhiqvpKQkpKSkqIxVrVq1N9ovIm1kIHUBRKQ5mzdvhhAC+/btw759+0p8vm3bNsyfPx/6+vqwtbVFUlJSiT7Pt508eRIPHz7E6dOnlUdrACAtLU3j9b+Ira0tLly4ACGESsBJTk5GUVER7OzsVPq/ap9sbW3Ltlg1ZGRkoFu3blAoFDh69Cisra1VPrezs4OtrS2OHz/+wu9bWlqqvH8+/BHpMoYbIi0hl8uxbds21KlTB99//32Jzw8fPoxvv/0Wx44dQ48ePdCuXTscOnQIjx49Ui7jyOVy7N69W+V7z/5oGhsbq7Rv2LChjPZEVYcOHbBnzx4cPHgQffv2Vbb/8MMPys//6ffff3/hPtWpUwc1a9YE8P/7ItXl0gUFBejbty/u3buHc+fOKev6px49emDXrl2Qy+Vo1qyZBFUSVV4MN0Ra4tixY3j48CEWL16Mtm3blvjc09MTa9aswaZNm9CjRw98+eWXOHToENq3b4+vvvoKZmZmWLt2bYnLu1u2bImqVasiICAAs2fPhqGhIbZv347IyMhy2a+RI0di7dq1GDVqFO7du4fGjRvj3LlzWLBgAbp164b3339fpb+dnR3at2+PWbNmwdzcHEFBQbh+/brK5eDP7tq8ceNGWFpawsTEBK6uruV2ZGfixIk4efIkFixYgKysLPz111/Kz+zt7VGnTh0MGTIE27dvR7du3TB+/Hg0bdoUhoaGePDgAU6dOoXevXurhD0i+n8854ZIS2zatAlGRkb48MMPX/i5nZ0d+vbti8OHD+PRo0fw9PTEb7/9BisrK4waNQqffvopvLy8MGvWLJXv2dra4siRIzAzM8O//vUv/Pvf/4aFhUWJIzxlxcTEBKdOncLw4cPxzTffoGvXrti6dSumTJmiPO/nn3r16oVx48bhyy+/RP/+/XHv3j1s375d5URpV1dXrFixApGRkWjbti3eeecd/PLLL+WyP0DxXaMBYMaMGWjRooXK6+uvvwZQfCL0oUOHMGPGDAQHB6Nv377o06cPFi1aBBMTEzRu3Ljc6iWqbGRCCCF1EUREmiCTyTB27FisWbNG6lKISEI8ckNERERaheGGiIiItApPKCYircFVdiICeOSGiIiItAzDDREREWkVhhsiIiLSKjp3zo1CocDDhw9haWnJ25UTERFVEkIIZGZmonr16tDTe/WxGZ0LNw8fPoSzs7PUZRAREdEbiI+Pf+EjS/5J58LNs4fNxcfHw8rKSuJqiIiIqDQyMjLg7Oxc4qGxL6Jz4ebZUpSVlRXDDRERUSVTmlNKeEIxERERaRWGGyIiItIqDDdERESkVRhuiIiISKsw3BAREZFWYbghIiIircJwQ0RERFqF4YaIiIi0CsMNERERaRWGGyIiItIqDDdERESkVRhuiIiISKsw3BAREZHGnLuVgvC4p5LWwHBDREREb02uEFgechMjNl/AuB3heJpdIFktBpJtmYiIiLTC48x8TNgdjj9vpwIA2tS3h6mRvmT1MNwQERHRGwu9k4rAXeF4nJkPMyN9LOjbGH18a0haE8MNERERqU2hEFj3xx18++sNKARQ39ECQcP9UNfBUurSGG6IiIhIPU+yCzBxdwT+uPkYADDAvybm9W4EM6OKESsqRhVERERUKVy69wTjdoQjKSMPJoZ6mNfbE4OaOEtdlgqGGyIiInothULgu7N3seTEDcgVAm725gga7gd3JyupSyuB4YaIiIheKS2nAFP2RuK3mGQAQG+f6vhv38awMK6YMaJiVkVEREQVQnjcU4zbEY6EtFwYGehhds+GGNa0FmQymdSlvRTDDREREZUghMCWP+9h4bEYFMoFXGzNsHaYHzxrWEtd2msx3BAREZGKjLxCfLH3Co5fSwIAdGvshEX9vWBlYihxZaXDcENERERKVxPSMWZ7GOKe5MBQX4aZ3TwwqmXtCr0M9TyGGyIiIoIQAj9diMPXv0SjQK5AjSqmCBruB2/nKlKXpjaGGyIiIh2XlV+Eafuv4PCVRADA+x6O+HagN6zNKscy1PMYboiIiHRYTGIGxm4Pw92UbBjoyTC1izs+bu1aqZahnsdwQ0REpIOEENhzKR5f/XwN+UUKVLM2wZphfvB3qSp1aW+N4YaIiEjH5BQU4csDVxEcngAAaNvAHssG+cDG3EjiyjSD4YaIiEiH3HqUiTHbw3ArOQt6MmBypwb4z3t1oKdXeZehnsdwQ0REpCOCwx5g5oGryC2Uw8HSGKuH+qKZm63UZWkcww0REZGWyyuUY/bP17D7UjwA4N26dlgxxAd2FsYSV1Y2GG6IiIi02N3HWRizPQzXkzIhkwETOtTHuPZ1oa9Fy1DPY7ghIiLSUociH2L6/ivILpDDzsIIK4f4olVdO6nLKnMMN0RERFomr1CO+Uei8dNfcQCAZq42WD3UFw5WJhJXVj4YboiIiLTI/dRsjN0RhqsJGQCAce3qYsL79WCgrydxZeWH4YaIiEhLHL+aiM/3XkFmfhGqmhli+WAftG3gIHVZ5U7yGBcUFARXV1eYmJjA398fZ8+efWnf06dPQyaTlXhdv369HCsmIiKqWAqKFJj7yzUE/BSGzPwi+LtUxZHA1joZbACJj9zs3r0bEyZMQFBQEFq1aoUNGzaga9euiI6ORq1atV76vRs3bsDKykr53t7evjzKJSIiqnAePM3B2B3hiIxPAwCMbuOGKZ0bwFCHlqGeJxNCCKk23qxZM/j5+WHdunXKNg8PD/Tp0wcLFy4s0f/06dNo164dnj59iipVqrzRNjMyMmBtbY309HSVgERERFTZ/Bb9CJP3RiI9txDWpob4dqA33m/oKHVZZUKdv9+SxbqCggJcvnwZnTp1Umnv1KkTzp8//8rv+vr6olq1aujQoQNOnTr1yr75+fnIyMhQeREREVVmhXIFFh6Nwcc/XEJ6biG8navg8Gfvam2wUZdk4SYlJQVyuRyOjqr/Qzg6OiIpKemF36lWrRo2btyI/fv3Izg4GA0aNECHDh1w5syZl25n4cKFsLa2Vr6cnZ01uh9ERETlKTE9F0M2/oUNZ+4CAD5sVRt7R7eAs42ZxJVVHJJfLSWTqd4hUQhRou2ZBg0aoEGDBsr3LVq0QHx8PJYuXYo2bdq88DvTp0/HpEmTlO8zMjIYcIiIqFI6fSMZk/ZE4kl2ASyNDbBkgBe6Nq4mdVkVjmThxs7ODvr6+iWO0iQnJ5c4mvMqzZs3x08//fTSz42NjWFsrJ3PziAiIt1QJFdgxW+3sObUbQBAo+pWCBruBxdbc4krq5gkW5YyMjKCv78/QkJCVNpDQkLQsmXLUo8THh6OatWYWomISDslZ+Rh+PcXlMHmX81rYf9/WjLYvIKky1KTJk3CiBEj0KRJE7Ro0QIbN25EXFwcAgICABQvKSUkJOCHH34AAKxYsQK1a9dGo0aNUFBQgJ9++gn79+/H/v37pdwNIiKiMnH+dgoCd4UjJasA5kb6WNjfC728q0tdVoUnabgZPHgwUlNTMW/ePCQmJsLT0xNHjx6Fi4sLACAxMRFxcXHK/gUFBZgyZQoSEhJgamqKRo0a4ciRI+jWrZtUu0BERKRxcoXAmpO3seL3mxACcHeyxNrhfqhjbyF1aZWCpPe5kQLvc0NERBVZSlY+JuyKwLnbKQCAwU2cMadXI5ga6UtcmbTU+fst+dVSREREVOzC3VR8tjMcyZn5MDXUx/w+nujvX1PqsiodhhsiIiKJKRQC68/cwbe/3oRcIVDXwQLrhvuhnqOl1KVVSgw3REREEnqaXYBJeyJw6sZjAEA/3xqY39cTZkb8E/2mOHNEREQSuXz/CcbtCEdieh6MDfQwt1cjDH7H+aU3s6XSYbghIiIqZ0IIbDoXi0XHrqNIIeBqZ46g4X7wqMYLXTSB4YaIiKgcpecUYsq+SIREPwIA9PSujoX9GsPCmH+SNYUzSUREVE4i49MwdkcYHjzNhZG+Hmb1bIh/NavFZSgNY7ghIiIqY0IIbDt/D/89GoNCuUAtGzOsHeaHxjWtpS5NKzHcEBERlaGMvEJM238FR6OKHxTduZEjlgzwhrWpocSVaS+GGyIiojJyNSEd43aE4V5qDgz0ZJjRzQMftqrNZagyxnBDRESkYUIIbL8Qh3mHo1FQpECNKqZYM8wXvrWqSl2aTmC4ISIi0qDMvEJMD47C4SuJAID3PRywdKA3qpgZSVyZ7mC4ISIi0pBrD9Mxbkc4YlOyYaAnw9Qu7vi4tSuXocoZww0REdFbEkJgx8U4zP2leBmqurUJVg/zg78Ll6GkwHBDRET0FrLyizAjOAqHIh8CANq7O+Dbgd6oas5lKKkw3BAREb2hmMQMjN0ehrsp2dDXk+GLzg3wSWs36OlxGUpKDDdERERqEkJg19/xmHPoGvKLFKhmbYI1w3zh72IjdWkEhhsiIiK1ZOcXYcaBKPwcUbwM1a6BPZYN8uEyVAXCcENERFRKMYkZGLsjDHcfFy9Dfd65AT7lMlSFw3BDRET0GkII7P47HrP/twzlZFW8DNWkNpehKiKGGyIiolfIzi/Clwev4kB4AgDgvfr2WD7YBzZchqqwGG6IiIhe4kZSJsZsv4w7/1uGmtypPgLa1OEyVAXHcENERPQcIQT2XnqArw5dRV5h8TLUqqG+aOrKZajKgOGGiIjoH3IKipehgsOKl6Ha1LfH8kHesLUwlrgyKi2GGyIiov+5+SgTY7aH4XZyFvRkwORODfCf97gMVdkw3BAREQHYeykes34uXoZysDTG6qG+aOZmK3VZ9AYYboiISKflFBRh1sFr2B/2AADQup4dlg/2gR2XoSothhsiItJZt/63DHXrf8tQkzrWx5i2dbkMVckx3BARkU7af/kBvjx4FbmFcjhYGmPVUF805zKUVmC4ISIinZJbIMdXP1/F3svFy1Dv1i1ehrK35DKUtmC4ISIinXE7uXgZ6uaj4mWoCe/Xx9h2daHPZSitwnBDREQ6ITjsAWYeKF6Gsrc0xsohPmhZx07qsqgMMNwQEZFWyy2QY/ahq9hzqXgZqlVdW6wY7MtlKC3GcENERFrrdnIWxm4Pw41HmZDJgAkd6mNcey5DaTuGGyIi0koHwxMw40AUcgrksLMwxqohPmhZl8tQuoDhhoiItEpeoRxzDl3Drr/jAQAt69hixRAfOFiaSFwZlReGGyIi0hp3HhcvQ11PKl6GCmxfD4Ed6nEZSscw3BARkVb4OSIB04OfLUMZYeUQX7TiMpROYrghIqJKLa9Qjrm/RGPnxTgAQHM3G6wa4gsHKy5D6SqGGyIiqrTuPs7CmH8sQ33Wvh7GcxlK5zHcEBFRpXQo8iGm77+C7AI5bM2NsGKID1rXs5e6LKoAGG6IiKhSySuU4+vD0dh+oXgZqpmrDVYN9YUjl6HofxhuiIio0ohNycaY7WGIScyATAaMa1cX4zvUg4G+ntSlUQXCcENERJXCL5EPMT04Cln5RbAxN8KKwT5oU5/LUFQSww0REVVoeYVyzD8SjZ/+Kl6GaupafDWUkzWXoejFGG6IiKjCupeSjbE7wnDtYQaA4mWoCe9zGYpejeGGiIgqpCNXEjF1/xXlMtTywT54j8tQVAoMN0REVKHkF8nx3yMx+CH0PgDgndpVsWqoL6pZm0pcGVUWDDdERFRh3E8tXoa6mlC8DDWmbR1M6lify1CkFoYbIiKqEI5GJWLqvivIzC9CVTNDLBvsg3YNHKQuiyohjYSbtLQ0VKlSRRNDERGRjskvkmPBkRhs+98yVBOXqlg9jMtQ9ObUPs63ePFi7N69W/l+0KBBsLW1RY0aNRAZGal2AUFBQXB1dYWJiQn8/f1x9uzZUn3vzz//hIGBAXx8fNTeJhERVQxxqTkYsC5UGWwC3quDnZ82Z7Cht6J2uNmwYQOcnZ0BACEhIQgJCcGxY8fQtWtXfP7552qNtXv3bkyYMAEzZ85EeHg4Wrduja5duyIuLu6V30tPT8fIkSPRoUMHdcsnIqIK4vjVRHRffRZRCemoYmaILR+8g2ld3WHI82voLcmEEEKdL5iamuLmzZtwdnbG+PHjkZeXhw0bNuDmzZto1qwZnj59WuqxmjVrBj8/P6xbt07Z5uHhgT59+mDhwoUv/d6QIUNQr1496Ovr4+DBg4iIiCj1NjMyMmBtbY309HRYWVmV+ntERKQZ+UVyLDx6HVvP3wMA+LtUxeqhvqhehUdr6OXU+futdjyuWrUq4uPjAQDHjx/H+++/DwAQQkAul5d6nIKCAly+fBmdOnVSae/UqRPOnz//0u9t2bIFd+7cwezZs9UtnYiIJBb/JAeD1ocqg83oNm7Y9WlzBhvSKLVPKO7Xrx+GDRuGevXqITU1FV27dgUAREREoG7duqUeJyUlBXK5HI6Ojirtjo6OSEpKeuF3bt26hWnTpuHs2bMwMChd6fn5+cjPz1e+z8jIKHWNRESkOcevJuHzfZHIzCtCFTNDfDvQGx08HF//RSI1qR1uli9fjtq1ayM+Ph5LliyBhYUFACAxMRFjxoxRuwCZTKbyXghRog0A5HI5hg0bhrlz56J+/fqlHn/hwoWYO3eu2nUREZFmFBQpsOjYdWz+MxYA4FerClYP80MNHq2hMqL2OTeaUlBQADMzM+zduxd9+/ZVto8fPx4RERH4448/VPqnpaWhatWq0NfXV7YpFAoIIaCvr49ff/0V7du3L7GdFx25cXZ25jk3RETlIP5JDsbtDEdkfBoA4NM2bvi8cwOeNExqK9NzbrZt24YjR44o33/xxReoUqUKWrZsifv375d6HCMjI/j7+yMkJESlPSQkBC1btizR38rKClFRUYiIiFC+AgIC0KBBA0RERKBZs2Yv3I6xsTGsrKxUXkREVPZ+vZaE7qvOIjI+Ddamhvh+ZBPM6ObBYENlTu2fsAULFsDUtPhQYmhoKNasWYMlS5bAzs4OEydOVGusSZMm4fvvv8fmzZsRExODiRMnIi4uDgEBAQCA6dOnY+TIkcWF6unB09NT5eXg4AATExN4enrC3Nxc3V0hIqIyUFCkwNeHo/Hpj5eRkVcEH+cqOBL4Lt5vyPNrqHyofc5NfHy88sThgwcPYsCAAfj000/RqlUrtG3bVq2xBg8ejNTUVMybNw+JiYnw9PTE0aNH4eLiAqD4PJ7X3fOGiIgqjgdPczBuRzgi/rcM9fG7rviiizuMDHi0hsqP2ufcODg44MSJE/D19YWvry8mTpyIkSNH4s6dO/D29kZWVlZZ1aoRvM8NEVHZ+C36ESbvjUR6biGsTAywdKA3OjVykros0hLq/P1W+8hNx44d8fHHH8PX1xc3b95E9+7dAQDXrl1D7dq136hgIiKqvArlCiw5fh3fnS2+GsrbuQrWDPWFs42ZxJWRrlL7OOHatWvRokULPH78GPv374etrS0A4PLlyxg6dKjGCyQiooorIS0XgzaEKoPNR++6Yu/oFgw2JCnJLgWXCpeliIg04/eYR5i0p3gZyvJ/y1CduQxFZaRMLwUHgLNnz+Jf//oXWrZsiYSEBADAjz/+iHPnzr3JcEREVIkUyhVYcDQGH227hPTcQnjXtMbRwNYMNlRhqB1u9u/fj86dO8PU1BRhYWHKG+RlZmZiwYIFGi+QiIgqjodpuRi8IRQbz9wFAHzYqjb2BrTkMhRVKGqHm/nz52P9+vX47rvvYGhoqGxv2bIlwsLCNFocERFVHCevP0K3VWcRFpcGSxMDrP+XP2b3bMTLvKnCUftqqRs3bqBNmzYl2q2srJCWlqaJmoiIqAIplCuw9Ncb2PBH8dEar5rWWDPUD7VsebSGKia1w021atVw+/btEpd9nzt3Dm5ubpqqi4iIKoCHabkI3BmOS/efAgA+aFkb07u5w9hA/zXfJJKO2uFm9OjRGD9+PDZv3gyZTIaHDx8iNDQUU6ZMwVdffVUWNRIRkQRO3UjGpN0ReJpTCEtjAywZ4IWujatJXRbRa6kdbr744gukp6ejXbt2yMvLQ5s2bWBsbIwpU6Zg3LhxZVEjERGVoyK5At+G3MS603cAAI1rWGPNMF+42PIZflQ5qHWfG7lcjnPnzqFx48YwMTFBdHQ0FAoFGjZsCAsLi7KsU2N4nxsiopdLTC9ehvr7XvEy1KgWLpjR3YPLUCS5Mnv8gr6+Pjp37oyYmBjY2NigSZMmb1UoERFVHKdvJGPSnkg8yS6AhbEBFvf3QncvLkNR5aP2slTjxo1x9+5duLq6lkU9RERUzorkCiz/7SbWnipehmpU3Qprh/mhth2XoahyUjvc/Pe//8WUKVPw9ddfw9/fH+bmqj/8XOohIqo8ktLzELgzHBfvPQEAjGjugpndPWBiyGUoqrzUfraUnt7/36xJJpMp/y2EgEwmg1wu11x1ZYDn3BARFTtz8zEm7o5A6v+WoRb1b4weXtWlLovohcrsnBsAOHXq1BsXRkRE0iuSK7Dit1tYe/o2hAAaVrPC2uF+cOUyFGkJtcONq6srnJ2dVY7aAMVHbuLj4zVWGBERad6jjOJlqAuxxctQ/2peC192b8hlKNIqbxRuEhMT4eDgoNL+5MkTuLq6VvhlKSIiXXX21mNM2FW8DGVupI+F/b3Qy5vLUKR91A43z86teV5WVhZMTEw0UhQREWmOXCGw8rebWH2qeBnKo5oVgrgMRVqs1OFm0qRJAIpPIp41axbMzP7/gWlyuRwXLlyAj4+PxgskIqI3l5yRh8Bd4fjrbvEy1LBmtfBVDy5DkXYrdbgJDw8HUHzkJioqCkZGRsrPjIyM4O3tjSlTpmi+QiIieiPnbqVgwu5wpGQVL0Mt6NcYvX1qSF0WUZkrdbh5dpXUhx9+iJUrV/IyaiKiCkquEFj5+y2sPnkLQgDuTpZYO9wPdewrx2NyiN6W2ufcbNmyReV9RkYGTp48CXd3d7i7u2usMCIiUl9yZh7G74xA6N1UAMDQprUwuyeXoUi3qB1uBg0ahDZt2mDcuHHIzc1FkyZNcO/ePQghsGvXLvTv378s6iQiotc4fzsFgbsikJKVDzMjfSzo2xh9fLkMRbpH7/VdVJ05cwatW7cGABw4cABCCKSlpWHVqlWYP3++xgskIqJXkysEVvx2E8M3XUBKVj4aOFri0Lh3GWxIZ6kdbtLT02FjYwMAOH78OPr37w8zMzN0794dt27d0niBRET0co8z8zFy8wWs+K34/Joh7zjj4NhWqOvA82tId6m9LOXs7IzQ0FDY2Njg+PHj2LVrFwDg6dOnvM8NEVE5+vN2Csb/bxnK1FAfC/p5oq9vTanLIpKc2uFmwoQJGD58OCwsLODi4oK2bdsCKF6uaty4sabrIyKi5zx/NVQDR0usHe6Lug6WUpdGVCGoHW7GjBmDZs2aIS4uDh07dlQ+JdzNzY3n3BARlbFHGXkY/4+b8g15xxmzezaCqRGvhiJ6RiaEEFIXUZ7UeWQ6EVFFcubmY0zc/f/PhuJN+UiXqPP3W+0jN0REVL6K5AosC7mJoNN3ABQ/G2rtMF+48aZ8RC/EcENEVIElpucicGc4/r73FADwr+a18GV33pSP6FUYboiIKqhT15MxaU8EnuYUwsLYAIv6N0YPr+pSl0VU4THcEBFVMIVyBZaeuIENZ+4CADxrWGHtMD+42JpLXBlR5fBG4SYtLQ0XL15EcnIyFAqFymcjR47USGFERLooIS0Xn+0IQ1hcGgDgg5a1Mb2bO4wNuAxFVFpqh5tffvkFw4cPR3Z2NiwtLSGTyZSfyWQyhhsiojcUEv0IU/ZGIj23EJYmBvhmgBe6eFaTuiyiSkftcDN58mT8+9//xoIFC2BmZlYWNRER6ZSCIgUWH7+OTediAQDeNa2xZpgfnG34O5boTagdbhISEhAYGMhgQ0SkAfFPcjBuRxgiH6QDAD561xVTu7jDyEDtR/8R0f+oHW46d+6MS5cuwc3NrSzqISLSGcevJuLzfVeQmVcEa1NDLB3ojY4NHaUui6jSUzvcdO/eHZ9//jmio6PRuHFjGBoaqnzeq1cvjRVHRKSN8ovkWHAkBttC7wMA/GpVwephfqhRxVTiyoi0g9qPX3j2LKkXDiaTQS6Xv3VRZYmPXyAiKd1Lyca4nWG4mpABABj9nhumdGoAQ30uQxG9Spk+fuH5S7+JiKh0Dl95iGn7o5CVX4SqZoZYNsgH7dwdpC6LSOvwJn5ERGUsr1COrw9HY/uFOADAO7WrYtVQX1Sz5jIUUVkoVbhZtWoVPv30U5iYmGDVqlWv7BsYGKiRwoiItMGdx1kYuz0M15MyIZMBY9rWwcT368OAy1BEZaZU59y4urri0qVLsLW1haur68sHk8lw9+5djRaoaTznhojKy8HwBMw4EIWcAjlszY2wfLAP2tS3l7osokpJ4+fcxMbGvvDfRERUUm6BHHN/uYZdf8cDAJq72WDlEF84WplIXBmRbuA5N0REGnQ7ORNjt4fjxqPiZajA9vUQ2KEe9PVkr/8yEWkEww0RkYbsu/wAsw5eRW6hHHYWxlg1xAct69pJXRaRzmG4ISJ6SzkFRZh18Br2hz0AALxb1w7LB/vA3tJY4sqIdBPDDRHRW7iRlImxO8JwOzkLejJg4vv1MaZdXS5DEUmI4YaI6A0IIbDnUjy++vka8osUcLQyxsohvmjuZit1aUQ6r1Th5sqVK6Ue0MvLS60CgoKC8M033yAxMRGNGjXCihUr0Lp16xf2PXfuHKZOnYrr168jJycHLi4uGD16NCZOnKjWNomI3kZWfhG+PBCFgxEPAQDv1bfHskHesLXgMhRRRVCqcOPj4wOZTAYhBGSyVx9qVefZUrt378aECRMQFBSEVq1aYcOGDejatSuio6NRq1atEv3Nzc0xbtw4eHl5wdzcHOfOncPo0aNhbm6OTz/9tNTbJSJ6U9EPMzBuRxjupmRDX0+GKZ0aYHQbN+hxGYqowijVTfzu37+v/Hd4eDimTJmCzz//HC1atAAAhIaG4ttvv8WSJUvQp0+fUm+8WbNm8PPzw7p165RtHh4e6NOnDxYuXFiqMfr16wdzc3P8+OOPperPm/gR0ZsQQmD7hTjMOxyNgiIFqlmbYPVQXzSpbSN1aUQ6QeM38XNxcVH+e+DAgVi1ahW6deumbPPy8oKzszNmzZpV6nBTUFCAy5cvY9q0aSrtnTp1wvnz50s1Rnh4OM6fP4/58+eXqj8R0ZvIzCvEtOAoHLmSCABo7+6Abwd6o6q5kcSVEdGLqH1CcVRU1AsfweDq6oro6OhSj5OSkgK5XA5HR0eVdkdHRyQlJb3yuzVr1sTjx49RVFSEOXPm4OOPP35p3/z8fOTn5yvfZ2RklLpGIqKrCekYuyMM91NzYKAnw9Qu7vjoXVcuQxFVYGo/uc3DwwPz589HXl6esi0/Px/z58+Hh4eH2gU8fw5Pac7rOXv2LC5duoT169djxYoV2Llz50v7Lly4ENbW1sqXs7Oz2jUSke4RQmDb+XvoF3Qe91NzUKOKKfYEtMAnPL+GqMJT+8jN+vXr0bNnTzg7O8Pb2xsAEBkZCZlMhsOHD5d6HDs7O+jr65c4SpOcnFziaM7znh05aty4MR49eoQ5c+Zg6NChL+w7ffp0TJo0Sfk+IyODAYeIXik9txBT913B8WvFv586NnTE0gHesDYzlLgyIioNtcNN06ZNERsbi59++gnXr1+HEAKDBw/GsGHDYG5uXupxjIyM4O/vj5CQEPTt21fZHhISgt69e5d6HCGEyrLT84yNjWFszMsziah0IuPTMG5nGOKf5MJQX4bpXT3wYavarz2iTEQVxxvdxM/MzEwjl15PmjQJI0aMQJMmTdCiRQts3LgRcXFxCAgIAFB81CUhIQE//PADAGDt2rWoVasW3N3dARTf92bp0qX47LPP3roWItJtQghs/vMeFh2LQaFcwNnGFGuG+sHbuYrUpRGRmt4o3Pz444/YsGED7t69i9DQULi4uGD58uVwc3NT66jL4MGDkZqainnz5iExMRGenp44evSo8uqsxMRExMXFKfsrFApMnz4dsbGxMDAwQJ06dbBo0SKMHj36TXaDiAgAkJZTgCl7r+C3mEcAgK6eTljU3wvWplyGIqqMSnWfm39at24dvvrqK0yYMAHz58/HtWvX4Obmhq1bt2Lbtm04depUWdWqEbzPDRH90+X7TxG4MxwJabkw0tfDrB4e+FdzFy5DEVUw6vz9VvtqqdWrV+O7777DzJkzYWDw/wd+mjRpgqioKPWrJSKSgEIhsOGPOxi8IRQJabmobWuG4DEtMaIFz68hquzUXpaKjY2Fr69viXZjY2NkZ2drpCgiorL0JLsAk/dE4NSNxwCAnt7VsaCvJyxNuAxFpA3UDjeurq6IiIhQuWsxABw7dgwNGzbUWGFERGXh73tP8NmOcCRl5MHYQA+zezbC0KbOPFpDpEXUDjeff/45xo4di7y8PAghcPHiRezcuRMLFy7E999/XxY1EhG9NYVCYN0fd7As5CbkCgE3e3OsHeYHj2o8945I26gdbj788EMUFRXhiy++QE5ODoYNG4YaNWpg5cqVGDJkSFnUSET0VlKy8jFxdwTO3koBAPT1rYH5fTxhbvxGF4wSUQWn9tVS/5SSkgKFQgEHBwdN1lSmeLUUkW4JvZOK8bvCkZyZDxNDPczr5YmBTWpyGYqokinTq6Xat2+PtLQ0AMWPUHgWbDIyMtC+fXv1qyUiKgNyhcDK325h+Pd/ITkzH/UcLHBo3LsY9A7PryHSdmofkz19+jQKCgpKtOfl5eHs2bMaKYqI6G0kZ+Zhwq4InL+TCgAY6F8Tc3s3gpkRl6GIdEGp/59+5coV5b+jo6NVHngpl8tx/Phx1KhRQ7PVERGp6c/bKRi/KwIpWfkwNdTHf/t6op9fTanLIqJyVOpw4+PjA5lMBplM9sLlJ1NTU6xevVqjxRERlVbxMtRNrD51G0IA7k6WWDPMD3UdLKQujYjKWanDTWxsLIQQcHNzw8WLF2Fvb6/8zMjICA4ODtDX1y+TIomIXuVRRh4Cd4bjQuwTAMDQps6Y3bMRTAz5O4lIF5U63Dy7aZ9CoSizYoiI1PXHzceYuDsCT7ILYG6kjwX9GqO3D5fIiXTZG59dFx0djbi4uBInF/fq1eutiyIiep0iuQLfhtzEutN3AAANq1lhzTBfuNlzGYpI16kdbu7evYu+ffsiKioKMpkMz26T8+zSSrlcrtkKiYie8zAtF4E7w3Hp/lMAwIjmLpjZ3YPLUEQE4A3uczN+/Hi4urri0aNHMDMzw7Vr13DmzBk0adIEp0+fLoMSiYj+38nrj9Bt1Vlcuv8UlsYGWDvMD1/38WSwISIltY/chIaG4uTJk7C3t4eenh709PTw7rvvYuHChQgMDER4eHhZ1ElEOq5QrsA3J25g45m7AIDGNayxZpgvXGzNJa6MiCoatcONXC6HhUXxmradnR0ePnyIBg0awMXFBTdu3NB4gURED57m4LOd4QiPSwMAfNCyNqZ3c4exAY/WEFFJaocbT09PXLlyBW5ubmjWrBmWLFkCIyMjbNy4EW5ubmVRIxHpsF+vJWHK3khk5BXBysQASwZ4o4unk9RlEVEFpna4+fLLL5GdnQ0AmD9/Pnr06IHWrVvD1tYWu3fv1niBRKSbCooUWHgsBlv+vAcA8HaugjVDfeFsYyZtYURU4b3VU8GfefLkCapWrVopHkbHp4ITVXxxqTkYtzMMVx6kAwA+ae2Kzzu7w8hA7WsgiEhLqPP3WyNPkbOxsdHEMEREOBaViC/2XUFmfhGqmBli6QBvvN/QUeqyiKgSKVW46devX6kHDA4OfuNiiEh35RXKseBoDH4IvQ8A8HepilVDfVGjiqnElRFRZVOqcGNtbV3WdRCRDruXko2xO8Jw7WEGACDgvTqY3Kk+DPW5DEVE6itVuNmyZUtZ10FEOuqXyIeYHhyFrPwi2Jgb4dtB3mjXwEHqsoioEtPIOTdEROrKK5Rj3uFo7LgQBwBoWtsGq4b6wsnaROLKiKiyUzvcuLq6vvKqqLt3775VQUSk/e48zsLY7WG4npQJmQwY164uxneoBwMuQxGRBqgdbiZMmKDyvrCwEOHh4Th+/Dg+//xzTdVFRFrqQPgDzDxwFTkFcthZGGH5YB+0rmcvdVlEpEXUDjfjx49/YfvatWtx6dKlty6IiLRTboEcsw9dxZ5LDwAALdxssXKIDxysuAxFRJqlsWPAXbt2xf79+zU1HBFpkVuPMtF77TnsufQAMhkwvkM9/PRxMwYbIioTGjuheN++fbyZHxGVsPdSPL76+RpyC+WwtzTGysE+aFnXTuqyiEiLqR1ufH19VU4oFkIgKSkJjx8/RlBQkEaLI6LKKzu/CLN+vorgsAQAQOt6dlg2yAf2lsYSV0ZE2k7tcNOnTx+V93p6erC3t0fbtm3h7u6uqbqIqBK7npSBsdvDcOdxNvRkwKSO9TGmbV3o6VX8588RUeWnkQdnViZ8cCZR2RFCYPff8Zh96BryixRwtDLGqiG+aOZmK3VpRFTJlemDM48ePQp9fX107txZpf3EiRNQKBTo2rWrukMSkRbIyi/CjOAoHIp8CAB4r749lg3yhq0Fl6GIqHypfbXUtGnTIJfLS7QLITBt2jSNFEVElcu1h+noufocDkU+hL6eDFO7uGPLB+8w2BCRJNQ+cnPr1i00bNiwRLu7uztu376tkaKIqHIQQuCnC3H4+nA0CooUqG5tgtXDfOHvwisniUg6aocba2tr3L17F7Vr11Zpv337NszNzTVVFxFVcBl5hZgeHIUjVxIBAO97OOCbAd6oam4kcWVEpOvUXpbq1asXJkyYgDt37ijbbt++jcmTJ6NXr14aLY6IKqaoB+noseocjlxJhIGeDF9298B3I5sw2BBRhaD2kZtvvvkGXbp0gbu7O2rWrAkAePDgAVq3bo2lS5dqvEAiqjiEENh2/h4WHL2OArkCNaqYYs0wX/jWqip1aURESm+0LHX+/HmEhIQgMjISpqam8PLyQps2bcqiPiKqINJzCvHF/kicuPYIANCpoSO+GeANazNDiSsjIlLF+9wQ0WtFxKdh3I4wPHiaC0N9GWZ088AHLWur3K2ciKgsafw+N6tWrSr1xgMDA0vdl4gqNiEENp2LxaJj11GkEKhlY4Y1w3zhVbOK1KUREb1UqY7cuLq6lm4wmQx3795966LKEo/cEJVOWk4BpuyNxG8xyQCAbo2dsKi/F6xMuAxFROVP40duYmNjNVIYEVUOl+8/wWc7wvEwPQ9GBnqY1aMh/tWsFpehiKhSUPuEYiLSXgqFwMazd/HNiRuQKwRc7cyxZpgvGlW3lro0IqJSe6Nw8+DBAxw6dAhxcXEoKChQ+WzZsmUaKYyIyldqVj4m743E6RuPAQC9vKtjQb/GsDDmfwMRUeWi9m+t33//Hb169YKrqytu3LgBT09P3Lt3D0II+Pn5lUWNRFTGLsY+wWc7w/AoIx/GBnqY06sRhrzjzGUoIqqU1L5D8fTp0zF58mRcvXoVJiYm2L9/P+Lj4/Hee+9h4MCBZVEjEZURhUJgzclbGLIxFI8y8uFmb46DY1thaFOeX0NElZfa4SYmJgajRo0CABgYGCA3NxcWFhaYN28eFi9erPECiahspGTlY9SWi1j6600oBNDPtwZ+GfcuPKrxKkIiqtzUXpYyNzdHfn4+AKB69eq4c+cOGjVqBABISUnRbHVEVCbO30nBhF0RSM7Mh4mhHub19sRA/5o8WkNEWkHtcNO8eXP8+eefaNiwIbp3747JkycjKioKwcHBaN68eVnUSEQaIlcIrPr9FladvAUhgHoOFgga7od6jpZSl0ZEpDFqh5tly5YhKysLADBnzhxkZWVh9+7dqFu3LpYvX67xAolIM5LS8zB+VzguxD4BAAxu4ow5vRrB1Ehf4sqIiDRL7XNu3Nzc4OXlBQAwMzNDUFAQrly5guDgYLi4uKhdQFBQEFxdXWFiYgJ/f3+cPXv2pX2Dg4PRsWNH2Nvbw8rKCi1atMCJEyfU3iaRrjl1IxndVp3FhdgnMDfSx8ohPlg8wIvBhoi0ktrh5u+//8aFCxdKtF+4cAGXLl1Sa6zdu3djwoQJmDlzJsLDw9G6dWt07doVcXFxL+x/5swZdOzYEUePHsXly5fRrl079OzZE+Hh4eruBpFOKJQrsPBoDD7c8jeeZBegYTUrHA5sjd4+NaQujYiozKj9VPCmTZviiy++wIABA1Tag4ODsXjx4hcGn5dp1qwZ/Pz8sG7dOmWbh4cH+vTpg4ULF5ZqjEaNGmHw4MH46quvStWfz5YiXRH/JAeBu8IRHpcGABjVwgXTu3nAxJBHa4io8tH4s6X+KTo6+oU36/P19UV0dHSpxykoKMDly5cxbdo0lfZOnTrh/PnzpRpDoVAgMzMTNjY2L+2Tn5+vvLoLKJ4cIm13/GoSvtgXiYy8IliZGGDJAC908awmdVlEROVC7WUpY2NjPHr0qER7YmIiDAxKn5VSUlIgl8vh6Oio0u7o6IikpKRSjfHtt98iOzsbgwYNemmfhQsXwtraWvlydnYudY1ElU1eoRyzf76KgJ8uIyOvCD7OVXAksDWDDRHpFLXDTceOHTF9+nSkp6cr29LS0jBjxgx07NhR7QKev6+GEKJU99rYuXMn5syZg927d8PBweGl/Z7V+uwVHx+vdo1ElUFsSjb6rzuPbaH3AQCj27hhb0ALONuYSVwZEVH5UntZ6ttvv0WbNm3g4uICX19fAEBERAQcHR3x448/lnocOzs76OvrlzhKk5ycXOJozvN2796Njz76CHv37sX777//yr7GxsYwNjYudV1EldHPEQmYERyF7AI5bMyN8O0gb7Rr8PLQT0SkzdQONzVq1MCVK1ewfft2REZGwtTUFB9++CGGDh0KQ0PDUo9jZGQEf39/hISEoG/fvsr2kJAQ9O7d+6Xf27lzJ/79739j586d6N69u7rlE2mV3AI55hy6ht2Xio9INnW1waohvnCyNpG4MiIi6agdboDiRzB8+umnb73xSZMmYcSIEWjSpAlatGiBjRs3Ii4uDgEBAQCKl5QSEhLwww8/ACgONiNHjsTKlSvRvHlz5VEfU1NTWFtbv3U9RJXJzUeZGLcjDDcfZUEmAz5rXw+B7evCQF/t1WYiIq3yRuFGUwYPHozU1FTMmzcPiYmJ8PT0xNGjR5U3A0xMTFS5582GDRtQVFSEsWPHYuzYscr2UaNGYevWreVdPpEkhBDYe+kBvjp0FXmFCthbGmPlYB+0rGsndWlERBWC2ve5qex4nxuqzLLyizDzQBR+jngIAGhdzw7LBvnA3pLnlRGRdivT+9wQkTSuJqTjs53hiE3Jhr6eDJM71UdAmzrQ0+OTvImI/onhhqiCE0Lgx7/uY/7hGBTIFahubYJVQ33RpPbLb15JRKTL1A438fHxkMlkqFmzJgDg4sWL2LFjBxo2bKiRk4yJ6P+l5xZi6r4rOH6t+OT59z0csXSgF6qYGUlcGRFRxaX2ZRXDhg3DqVOnAABJSUno2LEjLl68iBkzZmDevHkaL5BIV4XFPUW3lWdx/FoSDPVl+KpHQ3w30p/BhojoNdQON1evXkXTpk0BAHv27IGnpyfOnz+PHTt28IolIg1QKAQ2nrmDQetDkZCWi1o2Ztj/n5b497uupbp7NxGRrlN7WaqwsFB5x9/ffvsNvXr1AgC4u7sjMTFRs9UR6ZjUrHxM3huJ0zceAwC6e1XDwn6NYWVS+htkEhHpOrWP3DRq1Ajr16/H2bNnERISgi5dugAAHj58CFtbW40XSKQr/rqbim6rzuL0jccwNtDDgr6NsWaoL4MNEZGa1D5ys3jxYvTt2xfffPMNRo0aBW9vbwDAoUOHlMtVRFR6coXAmpO3sfL3m1AIoI69OdYO94O7E+/DRET0Jt7oJn5yuRwZGRmoWrWqsu3evXswMzN75RO6KwLexI8qkuSMPIzfFYHQu6kAgAH+NTGvdyOYGfEuDURE/1SmN/HLzc2FEEIZbO7fv48DBw7Aw8MDnTt3frOKiXTQmZuPMXF3BFKzC2BmpI/5fTzRz6+m1GUREVV6aoeb3r17o1+/fggICEBaWhqaNWsGQ0NDpKSkYNmyZfjPf/5TFnUSaY1CuQLLQm5i3ek7AAB3J0usHe6HOvYWEldGRKQd1D6hOCwsDK1btwYA7Nu3D46Ojrh//z5++OEHrFq1SuMFEmmThLRcDNn4lzLY/Kt5LRwc24rBhohIg9Q+cpOTkwNLS0sAwK+//op+/fpBT08PzZs3x/379zVeIJG2CIl+hCl7I5GeWwhLYwMs6u+F7l7VpC6LiEjrqH3kpm7dujh48CDi4+Nx4sQJdOrUCQCQnJzME3SJXiC/SI65v1zDJz9cQnpuIbxrWuNIYGsGGyKiMqJ2uPnqq68wZcoU1K5dG02bNkWLFi0AFB/F8fX11XiBRJXZ/dRsDFgXii1/3gMAfPyuK/YGtEQtWzNpCyMi0mJvdCl4UlISEhMT4e3tDT294nx08eJFWFlZwd3dXeNFahIvBafy8kvkQ0wPjkJWfhGqmBni24He6ODhKHVZRESVUpleCg4ATk5OcHJywoMHDyCTyVCjRg3ewI/of/IK5Zj7SzR2XowDALxTuypWDvFF9SqmEldGRKQb1F6WUigUmDdvHqytreHi4oJatWqhSpUq+Prrr6FQKMqiRqJK43ZyJvqs/RM7L8ZBJgPGtauLnZ80Z7AhIipHah+5mTlzJjZt2oRFixahVatWEELgzz//xJw5c5CXl4f//ve/ZVEnUYUmhMDeyw8w++dryC2Uw87CCMsH+6B1PXupSyMi0jlqn3NTvXp1rF+/Xvk08Gd+/vlnjBkzBgkJCRotUNN4zg1pWlZ+Eb48EIWDEQ8BAO/WtcOywd5wsDSRuDIiIu1RpufcPHny5IUnDbu7u+PJkyfqDkdUqV1NSMdnO8MRm5INfT0ZJnWsj/+8Vwd6ejKpSyMi0llqn3Pj7e2NNWvWlGhfs2aN8gnhRNpOCIGtf8aiX9B5xKZko7q1CXZ/2hxj29VlsCEikpjaR26WLFmC7t2747fffkOLFi0gk8lw/vx5xMfH4+jRo2VRI1GFkpZTgC/2XcGv0Y8AAO97OGLpQC9UMTOSuDIiIgLe4MjNe++9h5s3b6Jv375IS0vDkydP0K9fP9y4cUP5zCkibXX5/hN0X3UOv0Y/gpG+Hmb3bIjvRvoz2BARVSBvdBO/F4mPj8fs2bOxefNmTQxXZnhCMb0JhUJg3R93sCzkJuQKgdq2ZlgzzA+eNaylLo2ISCeo8/db7SM3L/PkyRNs27ZNU8MRVRiPM/MxastFfHPiBuQKgd4+1XE4sDWDDRFRBfVGdygm0hXnbqVgwu4IpGTlw8RQD/N6eWJgk5qQyXjSMBFRRcVwQ/QCRXIFVvx2C2tP34YQQANHS6wZ5ot6jpZSl0ZERK/BcEP0nIdpuQjcGY5L958CAIY2rYXZPRvCxFBf4sqIiKg0Sh1u+vXr98rP09LS3rYWIsmFRD/ClL2RSM8thKWxARb0a4ye3tWlLouIiNRQ6nBjbf3qkyetra0xcuTIty6ISAr5RXIsOnYdW/68BwDwqmmNNUP9UMvWTNrCiIhIbaUON1u2bCnLOogkcy8lG+N2huFqQgYA4ON3XfFFF3cYGWjsYkIiIipHPOeGdNrPEQmYeeAqsvKLUNXMEEsHeqODh6PUZRER0VtguCGdlFsgx5xD17D7UjwAoKmrDVYO8UE1a1OJKyMiorfFcEM650ZSJsbtCMOt5CzIZMBn7eshsH1dGOhzGYqISBsw3JDOEEJg19/xmHPoGvKLFLC3NMbKwT5oWddO6tKIiEiDGG5IJ2TmFWJ6cBQOX0kEALSpb49lg7xhZ2EscWVERKRpDDek9a48SMO4HeGIe5IDAz0ZPu/cAJ+0doOeHh+hQESkjRhuSGsJIbDpXCwWH7+OQrlAjSqmWD3MF361qkpdGhERlSGGG9JKT7MLMGVvJH6/ngwA6NLICYv7e8HazFDiyoiIqKwx3JDWuRj7BIE7w5GUkQcjAz3M6tEQ/2pWi0/yJiLSEQw3pDXkCoG1p25jxW83oRCAm5051gzzQ8PqVlKXRkRE5YjhhrRCckYeJuyOwPk7qQCAfn418HVvT5gb80eciEjX8Dc/VXp/3HyMSbsjkJpdADMjfXzd2xP9/WtKXRYREUmE4YYqrUK5At/+ehPr/7gDAHB3ssTa4X6oY28hcWVERCQlhhuqlB48zUHgznCExaUBAEY0d8HM7h4wMdSXtjAiIpIcww1VOsevJuKLfVeQkVcESxMDLOnvha6Nq0ldFhERVRAMN1Rp5BXKseBoDH4IvQ8A8HGugtVDfeFsYyZxZUREVJEw3FClcPdxFsbtCEd0YgYAYPR7bpjSqQEM+SRvIiJ6DsMNVXjBYQ/w5cGryCmQw9bcCN8O8kbbBg5Sl0VERBUUww1VWNn5Rfjq52vYH/YAANDCzRYrhvjA0cpE4sqIiKgik/yYflBQEFxdXWFiYgJ/f3+cPXv2pX0TExMxbNgwNGjQAHp6epgwYUL5FUrlKiYxAz3XnMP+sAfQkwGTOtbHTx83Y7AhIqLXkjTc7N69GxMmTMDMmTMRHh6O1q1bo2vXroiLi3th//z8fNjb22PmzJnw9vYu52qpPAgh8ONf99F77Z+4+zgbTlYm2PlJcwR2qAd9PT4bioiIXk8mhBBSbbxZs2bw8/PDunXrlG0eHh7o06cPFi5c+Mrvtm3bFj4+PlixYoVa28zIyIC1tTXS09NhZcVnDlUk6bmFmB58BUejkgAA7d0dsHSgN2zMjSSujIiIpKbO32/JzrkpKCjA5cuXMW3aNJX2Tp064fz58xrbTn5+PvLz85XvMzIyNDY2aU5Y3FN8tiMcCWm5MNSXYWoXd3z0riuf5E1ERGqTLNykpKRALpfD0dFRpd3R0RFJSUka287ChQsxd+5cjY1HmqVQCGw4cxdLf70BuUKglo0ZVg/1hbdzFalLIyKiSkryE4qf/y9zIYRG/2t9+vTpSE9PV77i4+M1Nja9nceZ+Ri15SIWH78OuUKgp3d1HAl8l8GGiIjeimRHbuzs7KCvr1/iKE1ycnKJozlvw9jYGMbGxhobjzTj7K3HmLg7EilZ+TAx1MO8Xp4Y2KQml6GIiOitSXbkxsjICP7+/ggJCVFpDwkJQcuWLSWqispaoVyBxcevY+Tmi0jJyoe7kyV+GfcuBr3jzGBDREQaIelN/CZNmoQRI0agSZMmaNGiBTZu3Ii4uDgEBAQAKF5SSkhIwA8//KD8TkREBAAgKysLjx8/RkREBIyMjNCwYUMpdoHUEP8kB4G7whH+vyd5/6t5LXzZvSGf5E1ERBolabgZPHgwUlNTMW/ePCQmJsLT0xNHjx6Fi4sLgOKb9j1/zxtfX1/lvy9fvowdO3bAxcUF9+7dK8/SSU3HohLxxf4ryOSTvImIqIxJep8bKfA+N+Urr1COrw9HY/uF4pDqW6sKVg3hk7yJiEg9leI+N6T9bj3KxLgd4bjxKBMyGRDwXh1M6lifT/ImIqIyxXBDGieEwJ5L8Zh96BryChWwszDG8sHeaF3PXurSiIhIBzDckEZl5hVixoGr+CXyIQCgdT07LBvkA3tLXo5PRETlg+GGNCYyPg2f7QxH3JMcGOjJMKVzA3za2g16fOAlERGVI4YbemsKhcCmc7FYfPw6ihQCNauaYtVQX/jVqip1aUREpIMYbuitpGblY/LeSJy+8RgA0K2xExb284K1qaHElRERka5iuKE3dv52CibsjkByZj6MDfTwVc+GGNa0Fu80TEREkmK4IbUVyRVY+fstrDl1G0IAdR0ssGaYL9ydeN8gIiKSHsMNqSUhLRcTdoXj73tPAQBD3nHG7J6NYGrERygQEVHFwHBDpXbiWhK+2HcF6bmFsDA2wIJ+jdHLu7rUZREREalguKHXyiuUY+HRGGwLvQ8A8K5pjdVD/VDLlo9QICKiiofhhl7pzuMsjNsRjpjEDADAp23cMKVTAxgZ8BEKRERUMTHc0Evtu/wAX/18FTkFctiaG2HpIG+0a+AgdVlERESvxHBDJWTlF2HWwas4EJ4AAGhZxxbLB/vA0cpE4sqIiIhej+GGVFxNSMe4HWG4l5oDfT0ZJr5fD/9pWxf6fIQCERFVEgw3BKD4Sd5b/ryHRceuo0CuQHVrE6wc6ot3attIXRoREZFaGG4IT7ML8Pm+SPwWkwwA6NTQEUsGeKGKmZHElREREamP4UbHXbibivG7IpCUkQcjfT182cMDI5q78BEKRERUaTHc6Ci5QmD1yVtY9fstKATgZm+O1UN90ai6tdSlERERvRWGGx2UlJ6H8bvCcSH2CQBggH9NzO3VCObG/HEgIqLKj3/NdMzvMY8wZW8knuYUwtxIH/P7eqKvb02pyyIiItIYhhsdkV8kx+JjN7D5z1gAgGcNK6we6gdXO3OJKyMiItIshhsdcC8lG5/tDEdUQjoA4N+tXDG1awMYG/BJ3kREpH0YbrTczxEJmBEchewCOaqaGeKbAd54v6Gj1GURERGVGYYbLZVTUITZP1/D3ssPAABNXW2wcogPqlmbSlwZERFR2WK40ULRDzMwbmcY7j7Ohp4MCOxQD5+1r8dHKBARkU5guNEiQgj8+Nd9zD8Sg4IiBRytjLFyiC+au9lKXRoREVG5YbjREmk5BZi6/wpOXHsEAOjg7oBvBnrDxpyPUCAiIt3CcKMFLt17gvG7IpCQlgtDfRmmd/XAh61q8xEKRESkkxhuKjG5QmD9H3ewLOQm5AqB2rZmWD3UD41r8hEKRESkuxhuKqnkjDxM3BOBP2+nAgD6+FTH/L6NYcFHKBARkY7jX8JK6PSNZEzeE4nU7AKYGupjXu9GGOBfk8tQREREYLipVAqKFFj66w1sPHMXAOBRzQqrh/qiroOFxJURERFVHAw3lURcag4+2xWOyPg0AMDIFi6Y0c0DJoZ8hAIREdE/MdxUAr9EPsSM4Chk5hfB2tQQi/t7oYunk9RlERERVUgMNxVYboEc8w5fw86L8QCAJi5VsXKoL2pU4SMUiIiIXobhpoK6kZSJcTvCcCs5CzIZMLZtXUx4vx4M9PWkLo2IiKhCY7ipYIQQ2HkxHnN/uYb8IgXsLY2xYrAPWtW1k7o0IiKiSoHhpgJJzy3EjOAoHIlKBAC8V98e3w7yhp2FscSVERERVR4MNxVEWNxTBO4Mx4OnuTDQk+GLLg3w8btu0OOTvImIiNTCcCMxhUJg49m7WHriBooUAs42plg91A8+zlWkLo2IiKhSYriR0OPMfEzaE4Gzt1IAAN29qmFhv8awMjGUuDIiIqLKi+FGIudupWDC7gikZOXDxFAPc3o2wuB3nPkIBSIiorfEcFPOCuUKLA+5iXV/3IEQQH1HC6wZ5of6jpZSl0ZERKQVGG7K0YOnOQjcGY6wuDQAwLBmtfBVj4Z8hAIREZEGMdyUk2NRiZi6/woy8opgaWKARf280N2rmtRlERERaR2GmzKWVyjH/CPR+OmvOACAj3MVrB7qC2cbM4krIyIi0k4MN2XodnImxu0Ix/WkTABAwHt1MLlTfRjyEQpERERlhuGmDAghsPfSA8w+dA25hXLYWRhh2SAftKlvL3VpREREWo/hRsMy8gox88BV/BL5EADwbl07LBvsDQdLE4krIyIi0g0MNxoUEZ+Gz3aGIf5JLvT1ZJjcqT4C2tThIxSIiIjKkeQnfwQFBcHV1RUmJibw9/fH2bNnX9n/jz/+gL+/P0xMTODm5ob169eXU6WvFnonFQPWnUf8k1zUqGKKPaNbYEzbugw2RERE5UzScLN7925MmDABM2fORHh4OFq3bo2uXbsiLi7uhf1jY2PRrVs3tG7dGuHh4ZgxYwYCAwOxf//+cq68JH+XqvCoZoVujZ1wdHxr+LtUlbokIiIinSQTQgipNt6sWTP4+flh3bp1yjYPDw/06dMHCxcuLNF/6tSpOHToEGJiYpRtAQEBiIyMRGhoaKm2mZGRAWtra6Snp8PKyurtd+KfY+cVwtLYgI9QICIi0jB1/n5LduSmoKAAly9fRqdOnVTaO3XqhPPnz7/wO6GhoSX6d+7cGZcuXUJhYeELv5Ofn4+MjAyVV1mxMjFksCEiIpKYZOEmJSUFcrkcjo6OKu2Ojo5ISkp64XeSkpJe2L+oqAgpKSkv/M7ChQthbW2tfDk7O2tmB4iIiKhCkvyE4uePdAghXnn040X9X9T+zPTp05Genq58xcfHv2XFREREVJFJdim4nZ0d9PX1SxylSU5OLnF05hknJ6cX9jcwMICtre0Lv2NsbAxjY2PNFE1EREQVnmRHboyMjODv74+QkBCV9pCQELRs2fKF32nRokWJ/r/++iuaNGkCQ0PDMquViIiIKg9Jl6UmTZqE77//Hps3b0ZMTAwmTpyIuLg4BAQEACheUho5cqSyf0BAAO7fv49JkyYhJiYGmzdvxqZNmzBlyhSpdoGIiIgqGEnvUDx48GCkpqZi3rx5SExMhKenJ44ePQoXFxcAQGJioso9b1xdXXH06FFMnDgRa9euRfXq1bFq1Sr0799fql0gIiKiCkbS+9xIoSzvc0NERERlo1Lc54aIiIioLDDcEBERkVZhuCEiIiKtwnBDREREWoXhhoiIiLQKww0RERFpFUnvcyOFZ1e+l+XTwYmIiEiznv3dLs0dbHQu3GRmZgIAnw5ORERUCWVmZsLa2vqVfXTuJn4KhQIPHz6EpaXlK58+XloZGRlwdnZGfHw8bwpYxjjX5YdzXX441+WHc11+ymKuhRDIzMxE9erVoaf36rNqdO7IjZ6eHmrWrKnxca2srPh/lnLCuS4/nOvyw7kuP5zr8qPpuX7dEZtneEIxERERaRWGGyIiItIqDDdvydjYGLNnz4axsbHUpWg9znX54VyXH851+eFclx+p51rnTigmIiIi7cYjN0RERKRVGG6IiIhIqzDcEBERkVZhuCEiIiKtovPh5syZM+jZsyeqV68OmUyGgwcPqnwuhMCcOXNQvXp1mJqaom3btrh27ZpKn/z8fHz22Wews7ODubk5evXqhQcPHqj0efr0KUaMGAFra2tYW1tjxIgRSEtLK+O9q1g0MdcbN25E27ZtYWVlBZlM9sI55Fy//Vw/efIEn332GRo0aAAzMzPUqlULgYGBSE9PVxmHc11MEz/bo0ePRp06dWBqagp7e3v07t0b169fV+nD+dbMXP+zb9euXV84DudaM3Pdtm1byGQyldeQIUNU+pTFXOt8uMnOzoa3tzfWrFnzws+XLFmCZcuWYc2aNfj777/h5OSEjh07Kp9RBQATJkzAgQMHsGvXLpw7dw5ZWVno0aMH5HK5ss+wYcMQERGB48eP4/jx44iIiMCIESPKfP8qEk3MdU5ODrp06YIZM2a8dDuc67ef64cPH+Lhw4dYunQpoqKisHXrVhw/fhwfffSRyjic62Ka+Nn29/fHli1bEBMTgxMnTkAIgU6dOvH3yHM0MdfPrFix4qWP4eFca26uP/nkEyQmJipfGzZsUPm8TOZakBIAceDAAeV7hUIhnJycxKJFi5RteXl5wtraWqxfv14IIURaWpowNDQUu3btUvZJSEgQenp64vjx40IIIaKjowUA8ddffyn7hIaGCgDi+vXrZbxXFdObzPU/nTp1SgAQT58+VWnnXJf0tnP9zJ49e4SRkZEoLCwUQnCuX0ZT8x0ZGSkAiNu3bwshON8v8jZzHRERIWrWrCkSExNLjMO5LulN5/q9994T48ePf+m4ZTXXOn/k5lViY2ORlJSETp06KduMjY3x3nvv4fz58wCAy5cvo7CwUKVP9erV4enpqewTGhoKa2trNGvWTNmnefPmsLa2VvbRdaWZ69LgXL/em851eno6rKysYGBQ/Eg6znXpvMl8Z2dnY8uWLXB1dYWzszMAzndplHauc3JyMHToUKxZswZOTk4lxuFcv546P9fbt2+HnZ0dGjVqhClTpqgc2SmruWa4eYWkpCQAgKOjo0q7o6Oj8rOkpCQYGRmhatWqr+zj4OBQYnwHBwdlH11Xmrku7Tic61d7k7lOTU3F119/jdGjR6uMw7l+PXXmOygoCBYWFrCwsMDx48cREhICIyMj5Tic71cr7VxPnDgRLVu2RO/evV86Duf61Uo718OHD8fOnTtx+vRpzJo1C/v370e/fv1UximLuda5p4K/iefXZIUQL12nfVmfF/UvzTi65k3m+nVjvOk42q60c52RkYHu3bujYcOGmD179ivHeNU4uq408z18+HB07NgRiYmJWLp0KQYNGoQ///wTJiYmLxzjZePoulfN9aFDh3Dy5EmEh4erNcbz41Cx1/1cf/LJJ8p/e3p6ol69emjSpAnCwsLg5+f3wjFeNI66eOTmFZ4drnw+PSYnJyvTqpOTEwoKCvD06dNX9nn06FGJ8R8/flwi9eqq0sx1acfhXL+aOnOdmZmJLl26wMLCAgcOHIChoaHKOJzr11Nnvq2trVGvXj20adMG+/btw/Xr13HgwAHlOJzvVyvNXJ88eRJ37txBlSpVYGBgoFxm7d+/P9q2basch3P9am/6O9vPzw+Ghoa4deuWcpyymGuGm1dwdXWFk5MTQkJClG0FBQX4448/0LJlSwDFVzgYGhqq9ElMTMTVq1eVfVq0aIH09HRcvHhR2efChQtIT09X9tF1pZnr0uBcv15p5zojIwOdOnWCkZERDh06pDx68AznunTe5mdbCIH8/HwAnO/SKM1cT5s2DVeuXEFERITyBQDLly/Hli1bAHCuS+NNf66vXbuGwsJCVKtWDUAZzvUbn4qsJTIzM0V4eLgIDw8XAMSyZctEeHi4uH//vhBCiEWLFglra2sRHBwsoqKixNChQ0W1atVERkaGcoyAgABRs2ZN8dtvv4mwsDDRvn174e3tLYqKipR9unTpIry8vERoaKgIDQ0VjRs3Fj169Cj3/ZWSJuY6MTFRhIeHi++++04AEGfOnBHh4eEiNTVV2Ydz/fZznZGRIZo1ayYaN24sbt++LRITE5Uv/lyX9LbzfefOHbFgwQJx6dIlcf/+fXH+/HnRu3dvYWNjIx49eqTcDudbM79HnofnrgQSgnMtxNvP9e3bt8XcuXPF33//LWJjY8WRI0eEu7u78PX1LfPfIzofbp5dUvz8a9SoUUKI4svdZs+eLZycnISxsbFo06aNiIqKUhkjNzdXjBs3TtjY2AhTU1PRo0cPERcXp9InNTVVDB8+XFhaWgpLS0sxfPjwEpcxaztNzPXs2bNfOMaWLVuUfTjXbz/XL/s+ABEbG6vsx7ku9rbznZCQILp27SocHByEoaGhqFmzphg2bFiJS2E535r5PfK8F4UbzvXbz3VcXJxo06aNsLGxEUZGRqJOnToiMDBQ5T9GhSibuZYJIcSbH/chIiIiqlh4zg0RERFpFYYbIiIi0ioMN0RERKRVGG6IiIhIqzDcEBERkVZhuCEiIiKtwnBDREREWoXhhojoDZw+fRoymQxpaWml/s4HH3yAPn36lFlNRFSM4YZIhwgh8P7776Nz584lPgsKCoK1tTXi4uIkqIyISHMYboh0iEwmw5YtW3DhwgVs2LBB2R4bG4upU6di5cqVqFWrloQVEhG9PYYbIh3j7OyMlStXYsqUKYiNjYUQAh999BE6dOgAV1dXNG3aFMbGxqhWrRqmTZuGoqIi5Xdr166NFStWqIzn4+ODOXPmKN/LZDJ8//336Nu3L8zMzFCvXj0cOnRI5TuHDh1CvXr1YGpqinbt2mHbtm0qSzz3799Hz549UbVqVZibm6NRo0Y4evToa/ft2VLRiRMn4OvrC1NTU7Rv3x7Jyck4duwYPDw8YGVlhaFDhyInJ0f5vfz8fAQGBsLBwQEmJiZ499138ffff6uMffToUdSvX19Z871791Q+nzNnDnx8fFTaVqxYgdq1a7+03tLM55w5c1CrVi0YGxujevXqCAwMfO08EOk6hhsiHTRq1Ch06NABH374IdasWYOrV69i5cqV6NatG9555x1ERkZi3bp12LRpE+bPn6/2+HPnzsWgQYNw5coVdOvWDcOHD8eTJ08AAPfu3cOAAQPQp08fREREYPTo0Zg5c6bK98eOHYv8/HycOXMGUVFRWLx4MSwsLEq9/Tlz5mDNmjU4f/484uPjMWjQIKxYsQI7duzAkSNHEBISgtWrVyv7f/HFF9i/fz+2bduGsLAw1K1bF507d1bWHB8fj379+qFbt26IiIjAxx9/jGnTpqk9L+rat28fli9fjg0bNuDWrVs4ePAgGjduXObbJar03uqxm0RUaT169EjY29sLPT09ERwcLGbMmCEaNGggFAqFss/atWuFhYWFkMvlQgghXFxcxPLly1XG8fb2FrNnz1a+ByC+/PJL5fusrCwhk8nEsWPHhBBCTJ06VXh6eqqMMXPmTAFA+STgxo0bizlz5qi9T8+eYvzbb78p2xYuXCgAiDt37ijbRo8eLTp37qysz9DQUGzfvl35eUFBgahevbpYsmSJEEKI6dOnCw8PD5W5mTp1qkrNs2fPFt7e3ir1LF++XLi4uCjfjxo1SvTu3Vv5/nXz+e2334r69euLgoICdaeCSKfxyA2RjnJwcMCnn34KDw8P9O3bFzExMWjRogVkMpmyT6tWrZCVlYUHDx6oNbaXl5fy3+bm5rC0tERycjIA4MaNG3jnnXdU+jdt2lTlfWBgIObPn49WrVph9uzZuHLlyhtv39HREWZmZnBzc1Npe1bPnTt3UFhYiFatWik/NzQ0RNOmTRETEwMAiImJQfPmzVXmpkWLFmrV9CYGDhyI3NxcuLm54ZNPPsGBAwdUlgmJ6MUYboh0mIGBAQwMDAAUX0n1zz/ez9oAKNv19PSUbc8UFhaWGNfQ0FDlvUwmg0KheO12nvn4449x9+5djBgxAlFRUWjSpInKMtLr/HP7MpnstfX8cx//WdOztufre5HSzo0633F2dsaNGzewdu1amJqaYsyYMWjTps1rxyXSdQw3RAQAaNiwIc6fP6/yx/b8+fOwtLREjRo1AAD29vZITExUfp6RkYHY2Fi1tuPu7l7iZN1Lly6V6Ofs7IyAgAAEBwdj8uTJ+O6779TaTmnVrVsXRkZGOHfunLKtsLAQly5dgoeHB4Diufnrr79Uvvf8e3t7eyQlJanMX0RExCu3XZr5NDU1Ra9evbBq1SqcPn0aoaGhiIqKUmsfiXQNww0RAQDGjBmD+Ph4fPbZZ7h+/Tp+/vlnzJ49G5MmTYKeXvGvivbt2+PHH3/E2bNncfXqVYwaNQr6+vpqbWf06NG4fv06pk6dips3b2LPnj3YunUrgP8/ejJhwgScOHECsbGxCAsLw8mTJ5VBQ9PMzc3xn//8B59//jmOHz+O6OhofPLJJ8jJycFHH30EAAgICMCdO3cwadIk3LhxAzt27FDW/Ezbtm3x+PFjLFmyBHfu3MHatWtx7NixV277dfO5detWbNq0CVevXsXdu3fx448/wtTUFC4uLhqfByJtwnBDRACAGjVq4OjRo7h48SK8vb0REBCAjz76CF9++aWyz/Tp09GmTRv06NED3bp1Q58+fVCnTh21tuPq6op9+/YhODgYXl5eWLdunfJqKWNjYwCAXC7H2LFj4eHhgS5duqBBgwYICgrS3M4+Z9GiRejfvz9GjBgBPz8/3L59GydOnEDVqlUBALVq1cL+/fvxyy+/wNvbG+vXr8eCBQtUxvDw8EBQUBDWrl0Lb29vXLx4EVOmTHnldl83n1WqVMF3332HVq1awcvLC7///jt++eUX2Nraan4SiLSITJRmMZmIqAz997//xfr16xEfHy91KUSkBQykLoCIdE9QUBDeeecd2Nra4s8//8Q333yDcePGSV0WEWkJLksRUbm7desWevfujYYNG+Lrr7/G5MmTVe7K+zIBAQGwsLB44SsgIKDsCyeiSoHLUkRUaSQnJyMjI+OFn1lZWcHBwaGcKyKiiojhhoiIiLQKl6WIiIhIqzDcEBERkVZhuCEiIiKtwnBDREREWoXhhoiIiLQKww0RERFpFYYbIiIi0ioMN0RERKRV/g9D/cilf0xMkQAAAABJRU5ErkJggg==", + "text/plain": [ + "
" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "# using adam optimizer\n", + "@jit\n", + "def compute_loss(youngs_modulus, solver, target_stress):\n", + " # creating a particle object with the given youngs modulus in the function\n", + " material_props = solver.mesh.particles[0].material.properties\n", + " material_props[\"youngs_modulus\"] = youngs_modulus\n", + " #updating the mesh object\n", + " solver.mesh.particles[0].material = LinearElastic(material_props)\n", + " gravity_loading = jnp.array([0.0, -9.8]).reshape(1, 2)\n", + " solver.mesh.particles[0].velocity = solver.mesh.particles[0].velocity.at[:].set(0.0)\n", + " #solving for the given youngs_modulus and computing the loss with the expected stress\n", + " result = solver.solve_jit(gravity_loading)\n", + " #using stress to calculate loss\n", + " loss = jnp.linalg.norm(result[\"stress\"] - target_stress)\n", + " return loss\n", + "\n", + "#adam optimizer\n", + "def optax_adam(params, niter, mpm, target_stress):\n", + " start_alpha = 0.1\n", + " optimizer = optax.adam(start_alpha)\n", + " opt_state = optimizer.init(params)\n", + " param_list = []\n", + " loss_list = []\n", + " t = tqdm(range(niter), desc=f\"Youngs_modulus: {params}\")\n", + " for _ in t:\n", + " lo, grads = value_and_grad(compute_loss)(params, mpm, target_stress)\n", + " updates, opt_state = optimizer.update(grads, opt_state)\n", + " params = optax.apply_updates(params, updates)\n", + " t.set_description(f\"Youngs_modulus: {params}\")\n", + " param_list.append(params)\n", + " loss_list.append(lo)\n", + " return param_list, loss_list\n", + "\n", + "#initial guess for the youngs modulus\n", + "params = 1050.0\n", + "#measuring the time taken by the optimizer\n", + "start_t = time.time()\n", + "#calling the optimizer\n", + "parameter_list, loss_list = optax_adam(params, 500, solver, target_ans[\"stress\"])\n", + "end_time = time.time()\n", + "print(f\"Time taken by adam: {end_time-start_t}\")\n", + "print(f\"parameter = {parameter_list[-1]} and loss = {loss_list[-1]}\")\n", + "plt.plot(parameter_list,loss_list)\n", + "plt.xlabel(\"Youngs_modulus\")\n", + "plt.ylabel(\"Loss calculated in stress\")\n", + "plt.title(\"Adam optimizer\")\n", + "plt.show()\n" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "The ADAM optimizer based method is takes a lot of time in estimating the youngs modulus which is not ideal for large simulations. Hence a much better method of Bayesian Optimization can be used which is very fast compared to ADAM optimizer, and generates a much better estimate of youngs modulus. In Bayesian Optimization method we need to provide a prior, a range where our youngs modulus would be estimated to lie, in most practical applications this information would be known. We also need to provide the number if random points Bayesian Optimization can explore and the number of iterations for which our it should run. \n", + "\n", + "In the following code block we will be using Bayesian Optimization to estimate the youngs modulus of the material. The number if iterations for which our solver is run is 8 while in ADAM optimizer it was 500, as the solver function is very costly hence minimizing the number if iterations would speed up the simulation." + ] + }, + { + "cell_type": "code", + "execution_count": 3, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "| iter | target | youngs... |\n", + "-------------------------------------\n", + "| \u001b[0m1 \u001b[0m | \u001b[0m-0.1767 \u001b[0m | \u001b[0m983.4 \u001b[0m |\n", + "| \u001b[0m2 \u001b[0m | \u001b[0m-0.4804 \u001b[0m | \u001b[0m1.044e+03\u001b[0m |\n", + "| \u001b[0m3 \u001b[0m | \u001b[0m-1.019 \u001b[0m | \u001b[0m900.0 \u001b[0m |\n", + "| \u001b[95m4 \u001b[0m | \u001b[95m-0.1635 \u001b[0m | \u001b[95m984.7 \u001b[0m |\n", + "| \u001b[0m5 \u001b[0m | \u001b[0m-1.108 \u001b[0m | \u001b[0m1.1e+03 \u001b[0m |\n", + "| \u001b[95m6 \u001b[0m | \u001b[95m-0.1306 \u001b[0m | \u001b[95m1.012e+03\u001b[0m |\n", + "| \u001b[95m7 \u001b[0m | \u001b[95m-0.002043\u001b[0m | \u001b[95m1e+03 \u001b[0m |\n", + "| \u001b[0m8 \u001b[0m | \u001b[0m-0.5549 \u001b[0m | \u001b[0m947.0 \u001b[0m |\n", + "=====================================\n", + "Time taken by Bayesian Optimization: 3.57847261428833\n", + "{'target': -0.0020428309217095375, 'params': {'youngs_modulus': 1000.1903192810634}}\n" + ] + } + ], + "source": [ + "# Using Bayesian Optimization model based on gaussian processes maximizes the\n", + "#given function using probability distributions\n", + "\n", + "#target stress for the given youngs modulus\n", + "target_stress = target_ans[\"stress\"]\n", + "\n", + "#loss function for bayesian optimization\n", + "@jit\n", + "def loss_func(youngs_modulus, solver=solver, target_stress=target_stress):\n", + " material_props = solver.mesh.particles[0].material.properties\n", + " material_props[\"youngs_modulus\"] = youngs_modulus\n", + " solver.mesh.particles[0].material = LinearElastic(material_props)\n", + " external_loading_local = jnp.array([0.0, -9.8]).reshape(1, 2)\n", + " result = solver.solve_jit(external_loading_local)\n", + " stress = result[\"stress\"]\n", + " loss = jnp.linalg.norm(stress - target_stress)\n", + " #returning negative of the loss as bayesian optimizer maximizes the function\n", + " return -loss\n", + "\n", + "#giving bound to the value of the youngs_modulus for bayesian optimizer\n", + "pbounds = {\"youngs_modulus\": (900, 1100)}\n", + "optimizer = BayesianOptimization(\n", + " f=loss_func, pbounds=pbounds, random_state=1, verbose=2\n", + ")\n", + "#measuring the time taken by the optimizer\n", + "start_t = time.time()\n", + "#calling the optimizer, init_points is the number of random points to be sampled\n", + "#n_iter is the number of iterations to be performed\n", + "optimizer.maximize(init_points=3, n_iter=5)\n", + "end_time = time.time()\n", + "print(f\"Time taken by Bayesian Optimization: {end_time-start_t}\")\n", + "print(optimizer.max)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "After analysing the result of both the estimations we observe that both methods correctly estimates the youngs modulus to be 1000 as stated in the original MPM solver. We can also notice that the ADAM solver took almost 50 seconds for the simulation while the bayesian optimization took only 3 seconds, this shows the significant benefit of using probability based methods over Linear Regression based methods. \n", + "\n", + "The Inverse Solver has many practical applications, it can be used to estimate the necessary properties that a material should possess to have a certain behaviour under different loading conditions." + ] + } + ], + "metadata": { + "kernelspec": { + "display_name": "optaximpo", + "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.4" + }, + "orig_nbformat": 4 + }, + "nbformat": 4, + "nbformat_minor": 2 +} diff --git a/examples/inverse_solver_2d.py b/examples/inverse_solver_2d.py new file mode 100644 index 0000000..ab40d2b --- /dev/null +++ b/examples/inverse_solver_2d.py @@ -0,0 +1,129 @@ +import sys +import os + +sys.path.append(os.getcwd()) +from diffmpm.material import SimpleMaterial, LinearElastic +from diffmpm.particle import Particles +from diffmpm.element import Quadrilateral4Node +from diffmpm.constraint import Constraint +from diffmpm.mesh import Mesh2D +from diffmpm.solver import MPMExplicit +import jax.numpy as jnp +import numpy as np +from jax import jit, value_and_grad +import optax +from tqdm import tqdm +import time +from bayes_opt import BayesianOptimization + +mesh_config = {} +density = 1 +poisson_ratio = 0.0 +youngs_modulus = 1000 +# Linear Elastic Material with given parameters +material = LinearElastic( + { + "id": 0, + "youngs_modulus": youngs_modulus, + "density": density, + "poisson_ratio": poisson_ratio, + } +) +# location of the material points +particle_loc = jnp.array([[0.0, 0.0], [0.5, 0.0], [0.0, 0.5], [0.5, 0.5]]).reshape( + 4, 1, 2 +) +# particle object for material points +particles = Particles( + particle_loc, material, jnp.zeros(particle_loc.shape[0], dtype=jnp.int32) +) +particles.velocity = particles.velocity.at[:].set(0.0) +# velocity constraint for material point at (0,0) +constraints = [(0, Constraint(1, 0.0))] +# gravity loading +gravity_loading = jnp.array([0.0, -9.8]).reshape(1, 2) +element = Quadrilateral4Node([1, 1], 1, [1, 1], constraints) +mesh_config["particles"] = [particles] +mesh_config["elements"] = element +mesh_config["particle_surface_traction"] = [] +mesh = Mesh2D(mesh_config) +solver = MPMExplicit(mesh, 0.01, sim_steps=10) +# getting the target value for youngs modulus 1000 +target_ans = solver.solve_jit(gravity_loading) + + +# using adam optimizer +@jit +def compute_loss(youngs_modulus, solver, target_stress): + # creating a particle object with the given youngs modulus in the function + material_props = solver.mesh.particles[0].material.properties + material_props["youngs_modulus"] = youngs_modulus + #updating the mesh object + solver.mesh.particles[0].material = LinearElastic(material_props) + gravity_loading = jnp.array([0.0, -9.8]).reshape(1, 2) + solver.mesh.particles[0].velocity = solver.mesh.particles[0].velocity.at[:].set(0.0) + #solving for the given youngs_modulus and computing the loss with the expected stress + result = solver.solve_jit(gravity_loading) + #using stress to calculate loss + loss = jnp.linalg.norm(result["stress"] - target_stress) + return loss + +#adam optimizer +def optax_adam(params, niter, mpm, target_stress): + start_alpha = 0.1 + optimizer = optax.adam(start_alpha) + opt_state = optimizer.init(params) + param_list = [] + loss_list = [] + t = tqdm(range(niter), desc=f"E: {params}") + for _ in t: + lo, grads = value_and_grad(compute_loss)(params, mpm, target_stress) + updates, opt_state = optimizer.update(grads, opt_state) + params = optax.apply_updates(params, updates) + t.set_description(f"E: {params}") + param_list.append(params) + loss_list.append(lo) + return param_list, loss_list + +#initial guess for the youngs modulus +params = 1050.0 +#measuring the time taken by the optimizer +start_t = time.time() +#calling the optimizer +parameter_list, loss_list = optax_adam(params, 500, solver, target_ans["stress"]) +end_time = time.time() +print(f"Time taken by adam: {end_time-start_t}") +print(f"parameter = {parameter_list[-1]} and loss = {loss_list[-1]}") + +# Using Bayesian Optimization model based on gaussian processes maximizes the +#given function using probability distributions + +#target stress for the given youngs modulus +target_stress = target_ans["stress"] + +#loss function for bayesian optimization +@jit +def loss_func(youngs_modulus, solver=solver, target_stress=target_stress): + material_props = solver.mesh.particles[0].material.properties + material_props["youngs_modulus"] = youngs_modulus + solver.mesh.particles[0].material = LinearElastic(material_props) + external_loading_local = jnp.array([0.0, -9.8]).reshape(1, 2) + result = solver.solve_jit(external_loading_local) + stress = result["stress"] + loss = jnp.linalg.norm(stress - target_stress) + #returning negative of the loss as bayesian optimizer maximizes the function + return -loss + +#giving bound to the value of the youngs_modulus for bayesian optimizer +pbounds = {"youngs_modulus": (900, 1100)} +optimizer = BayesianOptimization( + f=loss_func, pbounds=pbounds, random_state=1, verbose=2 +) +#measuring the time taken by the optimizer +start_t = time.time() +#calling the optimizer, init_points is the number of random points to be sampled +#n_iter is the number of iterations to be performed +optimizer.maximize(init_points=3, n_iter=5) +end_time = time.time() +print(f"Time taken by Bayesian Optimization: {end_time-start_t}") +print(optimizer.max)