diff --git a/examples/example-04.ipynb b/examples/example-04.ipynb index fb04f86..800cff4 100644 --- a/examples/example-04.ipynb +++ b/examples/example-04.ipynb @@ -1,159 +1,335 @@ { - "nbformat": 4, - "nbformat_minor": 4, - "metadata": {}, - "cells": [ - { - "metadata": {}, - "source": [ - "# Example 4: Simulated site profile\n", - "\n", - "Generate simulated shear-wave velocity profiles." - ], - "cell_type": "markdown" - }, - { - "metadata": {}, - "source": [ - "import matplotlib.pyplot as plt\n", - "\n", - "import pystrata\n", - "\n", - "%matplotlib inline" - ], - "cell_type": "code", - "outputs": [], - "execution_count": null - }, - { - "metadata": {}, - "source": [ - "# Increased figure sizes\n", - "plt.rcParams[\"figure.dpi\"] = 120" - ], - "cell_type": "code", - "outputs": [], - "execution_count": null - }, - { - "metadata": {}, - "source": [ - "Create a simple site profile" - ], - "cell_type": "markdown" - }, - { - "metadata": {}, - "source": [ - "profile = pystrata.site.Profile(\n", - " [\n", - " pystrata.site.Layer(\n", - " pystrata.site.SoilType(\"Soil-1\", 18.0, None, 0.05), 30, 400\n", - " ),\n", - " pystrata.site.Layer(\n", - " pystrata.site.SoilType(\"Soil-2\", 19.0, None, 0.05), 20, 600\n", - " ),\n", - " pystrata.site.Layer(pystrata.site.SoilType(\"Rock\", 24.0, None, 0.01), 0, 1200),\n", - " ]\n", - ")" - ], - "cell_type": "code", - "outputs": [], - "execution_count": null - }, - { - "metadata": {}, - "source": [ - "Initialize the variations." - ], - "cell_type": "markdown" - }, - { - "metadata": {}, - "source": [ - "toro_thickness = pystrata.variation.ToroThicknessVariation()\n", - "toro_velocity = pystrata.variation.ToroVelocityVariation.generic_model(\"USGS B\")" - ], - "cell_type": "code", - "outputs": [], - "execution_count": null - }, - { - "metadata": {}, - "source": [ - "Create the varied thickness and velocity." - ], - "cell_type": "markdown" - }, - { - "metadata": {}, - "source": [ - "%pdb" - ], - "cell_type": "code", - "outputs": [ - { - "name": "stdout", - "output_type": "stream", - "text": [ - "Automatic pdb calling has been turned ON\n" - ] - } - ], - "execution_count": null - }, - { - "metadata": {}, - "source": [ - "count = 10\n", - "# Create realizations of the profile with varied thickness\n", - "varied_thick = [toro_thickness(profile) for _ in range(count)]\n", - "\n", - "# For eaach realization of varied thickness, vary the shear-wave velocity\n", - "varied_vel = [toro_velocity(rt) for rt in varied_thick]" - ], - "cell_type": "code", - "outputs": [], - "execution_count": null - }, - { - "metadata": {}, - "source": [ - "Create a plot of the varied velocity models." - ], - "cell_type": "markdown" - }, - { - "metadata": {}, - "source": [ - "fig, ax = plt.subplots()\n", - "\n", - "for profile in varied_vel:\n", - " ax.plot(\n", - " [layer.initial_shear_vel for layer in profile],\n", - " [layer.depth for layer in profile],\n", - " drawstyle=\"steps-pre\",\n", - " )\n", - "\n", - "ax.set(xlabel=\"$V_s$ (m/s)\", xscale=\"log\", ylabel=\"Depth (m)\", ylim=(55, 0))\n", - "\n", - "ax.grid()\n", - "fig.tight_layout();" - ], - "cell_type": "code", - "outputs": [ - { - "data": { - "image/png": "iVBORw0KGgoAAAANSUhEUgAAAvQAAAI0CAYAAAByNL2iAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjkuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/GU6VOAAAACXBIWXMAABJ0AAASdAHeZh94AABlzUlEQVR4nO3deXhU5d3/8c9kJslkmSSEEBKWBBQkoEQlgFhFwZXNXdRiRWOtFNSnVq1SaxFaa9UHtP60Kq021n0D2keKWJeodQMhWqIQBAwBDCFAMtnIJJnJ/P4IRGMSMhNm5sxJ3q/r4jow55z7/p4lwyf3nDnH4vV6vQIAAABgShFGFwAAAACg+wj0AAAAgIkR6AEAAAATI9ADAAAAJkagBwAAAEyMQA8AAACYGIEeAAAAMDECPQAAAGBiBHoAAADAxAj0AAAAgIkR6AEAAAATM22gr62t1c0336wBAwbIbrfrhBNO0EsvvWR0WQAAAEBI2YwuoLsuvvhiffbZZ7rvvvt0zDHH6IUXXtCPf/xjNTc3a9asWUaXBwAAAISExev1eo0uwl+rVq3S9OnTW0P8Ieecc46++uor7dixQ1ar1cAKAQAAgNAw5SU3K1asUHx8vGbOnNnm9dzcXJWWlmrNmjUGVQYAAACElikD/ZdffqmRI0fKZmt7xVB2dnbrfAAAAKA3MOU19Pv379dRRx3V7vXk5OTW+YdTXl6uvXv3tnmturpaX3/9tUaPHq3o6OjAFQsAAAD4qKGhQTt37tTpp5+upKQkn9YxZaCXJIvF0q15kvTYY49p0aJFgS4JAAAACIh//OMfuuCCC3xa1pSBvm/fvh2OwldUVEj6bqS+M/PmzWt3/f3GjRt12WWX6YUXXuhw9D9Y7KtvkW3Pf+Xuf7xcUx4MWb9mVf/6Dnl218uaHqOY8zKMLqdTq1ev1p49e9S/f39NmTKlzbzi4gWqO7BRcbGjNHTo74Jax8cvF2v/rgPqOyhWP7p8aFD7Qnu/ff1rfbm7RselO/T7847pcJn5xeUqPNCo0bFRum9oarv5f1j/B21ybtLIpJGa+y+PmjZuUuSokUpesCDY5aMXez/vCe3d/o36DTlKp+f+3OhyEGL19fUqLCzU6NGjFRMTY3Q5vc4333yjWbNmafDgwT6vY8pAP3r0aL344otyu91trqMvLCyUJB133HGHXT81NVWpqe3/45RarsM/9thjA1dsVzYmSV6rlJEknXRS6Po1qfIvYtTorlLU4ESlnpRtdDmd2rhxo7xerzIyMnTSD46rLbKvnM4oJSX1Vc6Y4B7z0g8jFXXAqQEDk3TSSWOC2le4Kaiu00Pb96jW4zGshpKTR8nrcqvEbtN9UQkdLvPtwBjZ3M1KSorTSScObze/b2Vf2ffY1bd/X41K9uhAdLRik/sqk/cLBNH21StkrdyrQen9272Hoeerrq5WVVWVxowZo4SEjt+7EDzx8fGS5Ncl4KYM9BdddJH++te/atmyZbr88stbX//73/+uAQMG8OYDQA9t36O39lcbW0SsVYq1qkrSJ866wy4az612AQDdZMpAP3XqVJ199tmaO3euqqurNWzYML344otavXq1nnvuOe5BD6B1ZD7BFqFj4435yHhjabWqXW4l2G0aNaDzUa54q1W3DEkLYWUAgJ7ElIFekpYvX67f/OY3WrBggSoqKpSVlaUXX3xRV1xxhdGlAQgjx8bHaEUHl7KEwuVrP9Ga4gqdMDRZL083pgYAQM9n2kAfHx+vhx9+WA8//LDRpQAAAACGMW2gR+/WtLtW5Us3hKSviGirEs7MUNRgR0j6AwAA8AeBHqYSEd3y/Qivy6PG4qqQ9p1yTQjvfgQAAOAjAj1MJeHMlnvPNzeE5laETbtr5XV5QtYfAACAvwj0MJWowY6QjpSXL90Q8k8CAAAA/BFhdAEAAAAAuo9ADwAAAJgYgR4AAAAwMQI9AAAAYGIEegAAAMDECPQAAACAiRHoAQAAABMj0AMAAAAmRqAHAAAATIxADwAAAJgYgR4AAAAwMQI9AAAAYGI2owsA4JuC6jo9tH2Paj0en9fZN9yixkEORcVY9MznW4JYXfj5qrbe6BIAAAgJAj3gg6bdtSpfusH35ffXHVyvrt16TQNrpVipqbRW5Z/53uZ9A5r1QbzPi7dwWCRHpCRpi7POz5V7hshddSpf2/l+joi2KuHMDEUNdoSwKgAAAodADxxGRLRVkuR1edRYXOXzes1RbilCana5263X3NctxXY873BqU2Ik2RTf5NUxNb6P0vdmcW7pp9sa1Fjd3OWyKdccG4KKAAAIPAI9cBgJZ2ZIkpob/AvQEfttUqMUYbcpamBi23l2W+s0amhiR6t33Ka9JZRmNVn0t31RPq2zb1eNGuvdioqxKWVQLx2B7hsl9e14VtPuWnldHr+PLwAA4YRADxxG1GBHt0ZuI/PWSyX7FJkep9Tc7DbzdhbES04pckC8Umdkd9xAR21+vkVy1ilyQJxSpw/3aZ2PlhSodK9TA9IdGjXH9756i/KlG/z6lAQAgHBEoAfQ6/n7HQlfRERbNbjBqzUBbRUAgPYI9EAPtae4WqVbnEaXEda6+x0JX50dY9FrAW8VAIC2CPRAD7VuVXHr36PsVgMrCV/d/Y5EVw5dmx/t9Qa0XQAAOkKgB3qoRtd3IXXs9KEGVhK+uvsdia5wbT4AIJR4UizQww0YnqT+QxKMLgMAAAQJI/QAJEmNO2tU/c4ObuEYAE27a40uAQDQixDoAUiSqt/ZIVdRhdFl9CgNFovRJQAAegECPQBJ330x1GK3KjI93uBqzC8i2qq3qqukeqMrAQD0dAR6wAeFewu1dMNS1TXV+bR8/4r+ssuuzRWblbs6t828s21F6h8hba4o0qM/mHc4myMvkSIGHWzz3i6XH1FxthxK0+aKIuWufqTL5a+tmKqhStc30bv0t8w3fK4LndvorVZMRpO2R0Uqd3X3vsewuWJzgKsCAPQ0BHrAB0s3LNX7u973efnTmk6TXXbVNNVo3Z51beZN6OdSf7tU01ijdXvXddJCezWp50iH1ivver0BTRPkUFqHNXRkZtOpktJ9Xh4+iJBscdIBSev2HFlTcZFxkqoDURUAoIch0AM+ODQy74h0aETyiC6Xd+xzSK6W5cf2H9t2nu2/kqrkiHJobP/jfa5hTaRDFdLB9cZ2ubzja0drzT4tX+rf8ujaxt3Vqq5vUkJMpEald/9OQ3GRcfr58T+X9EDgigMA9BgEesAPI5JHKG9KXpfL5e3OU0lNiUYkj9B9U+5rM299wSw5nWs0IjlLs8Z03dYhF32+RZ8461pqOLHr9VZ8VaDS/U6NSM7S/Cmzuly+vGSDGquqNCI5S3lTLvO5LnTu8qWfaM2OCmUPTVbelJOPuL2SANQEAOh5uA89AAAAYGIEegAAAMDECPQAAACAiRHoAQAAABMj0AMAAAAmxl1uALTRtLtW5Us3hKSviGirEs7MUNRgR0j6AwCgJyLQA5DUEq4lyevyqLG4KqR9p1xzbEj7AwCgJyHQA5AkJZyZIUlqbvCEpL+m3bXyujwh6w8AgJ6KQA9AkhQ12BHSkfLypRtC/kkAAAA9EYEeXdpTXK11q4rV6Oq9I6kjKs7WgKYJcnzt0IqvCrpcft+B2pbprlqtWNJ2+dhhNbI5pH07a7Qiv+u2WtscbpEcFu3bWasV7/pQw65an9sGAADmRaBHl9atKtb2wv1Gl2Eoh9LkUJokqXS/s8vlG5PdUpTUWO9W6bdtl88Y6JbNcXDelq7bam1zkENyRB5cr8bn9aLsVp+XBQAA5kOgR5cOjcxHxdiUMije4GqMsbmiSDVNNXJEOjQiOavL5Q8csKnJ07LPBgxPajMvKsbWOv3hvMOJirH4vV6U3aqx04f63AcAADAfAj18ljIoXhfdOsboMgyRu/oRrduzTmP7j9X8KbO6XD4v77+qK9nfss9y2+6z9QUOOZ1SymCHci7wfX8+8/kWbXHWKWVwvC46f7i/mwAAAHooHiwFAAAAmBgj9AarL3Vp3wfJaraWSu/ONrqcDrnsUyRrmlxFRSq56k9Gl2OIKyqKNL3JLUfkVyp5vuvj1DAgXYqJUUNRkUquaru8a8omKV1ybdqkkod8P+YNU2dK6YPVsKlIJQ/+3u9tCDcRKdNliR6ghqJNKrlqsdHlBMU1u6t1cX2TEtZFquTDhCNuz1VUFICqAAA9DYHeYPs+dqq21C7JJe38zOhyOtR8wo+kpDQ111TrwBfhWWOwZbT+rUYHtnW9DzxnTJZiYuSpqdGBz9ou33xKk5QuNXcw77BtnjpFSleHbZpRzKmnyRbdsj31PWB7OjLke38/sCtw7UbExQWuMQCA6RHoDdbc2CxJioiOkD07x+BqOhZhbxlZjHAkKHbcOIOrMUbR974Um+XDl2KtDkfr9If7LMJRKKlaEQ6HYseN9rmGw7VpRhE9bHs6snF3tarrm5QQE6lR6Uc+Qi+1hPmUG+YFpC0AQM9AoA8T9tQoZT77jNFldKhgSYG0xSl7VpYyb+36C6E90cLVuQe/FHus8qbkdbl8dF6eVFKi6KwsZd5/f5t5+wpmyeVcI/vIkcq80vdjHv35FslZp+iRWcqcFZ7nij8OPVgqOmukUudcbnQ5QXH70k+0prhCJw1N1stzTja6HABAD8WXYgEAAAATI9ADAAAAJkagBwAAAEyMQA8AAACYGIEeAAAAMDECPQAAAGBi3LYyXDTWSXnTja6iY2WXSRoslRVKeb81uhr/RMdLp90uDQrPe/wDAAAcKQK90SKsLdNmt1Ty4WEXLYyK0tI+iaqzWEJQ2HdGNB+QQ9Lm5gPKbfg6pH0fsQZJK6+Q7IlH1MxmNUoWHfylxodfvMpGSEroePn+JZJdUtkG/36JS/upZD/qYJs3+75euCqbJSnj4Pb8uvPl+KUMAIDDItAbLXGwpKKWwJk58rCLLrWU632LKzR1fc8Aa4QckmqsEVoXYw95/4HREJBW4uqd0p4tPiyZJilBaqhq/4taYqJkj5Rc1V3+EtdG0qWSXfrK1k8XJV3qR9XhqSlhkJoVowgNUh/vdP2y5O8aU1PU+QqzXg5dcQAAmAiB3mjRjpZp2mgp9/BP/6xbnSvtWSdHpEMjkkeEoLgWjq9banREOjS2/9iQ9XvEGmqkqp1SsycgzcXJop/b+0uZPuz7ssSW3yGiE6W0U9vOs5dIqpfsCVLmaJ/7j4+KkSRV2xz6JOlE3wsPe4mSfiTFJOvZ8mfbziorbPmlqKHWkMoAADADAr0JjUgeobwpeSHrb8VXBSrd79SI5CzNnzIrZP2aWl6eVFJy8Be1B9vOK5glOddIadnStBd8bvKX1XXS9j2q9QTmFxSjNZXWqdnl1pZEm2qsUu2AHGn6FW0Xypvu36cYAAD0QgR6wCTGJMTp2eyjjC4jYMrXblBjcb1+PtGhdbFGVwMAgHlx20oAAADAxAj0AAAAgIkR6AEAAAATI9ADAAAAJkagBwAAAEyMQA8AAACYGIEeAAAAMDECPQAAAGBiBHoAAADAxAj0AAAAgIkR6AEAAAATsxldwA/V1NTo97//vb744gt9/vnn2rdvn+6++24tXLiw3bIFBQW6/fbb9emnn8pms+mMM87Q4sWLddRRR4W+8BDaXLFZuatzQ9bfiIqz5VCaNlcUKXf1IyHr18z6V/SXXfYOj9XZtiL1j5A2VxTp0RAex3BzbcVUDVW6DrgPSIo9uK/ubbuQZY+UltoyNeG+2h5VrZiMJm2PilTu6gSjyzFEXGSc5mTP0eh+o40uBQB6rLAL9Pv379df/vIXHX/88brwwgv15JNPdrhcUVGRJk2apBNOOEGvvPKKXC6XFixYoIkTJ+qLL75Qv379Qlx58MVFxkmSappqtG7PupD1O6BpghxKC3m/ZnZa02myy97hPpvQz6X+dqmmsUbr9vbe/Tmz6VRJ6fI0eyQd3B/lP9gfFkkxdkkNkhnPvQjJFicdkLRuj9HFGOvRMx81ugQA6LHCLtBnZmaqsrJSFotF+/bt6zTQL1iwQNHR0Vq5cqUSElpGvnJycjR8+HAtXrxY999/fyjLDok52XMkSXVNdSHt1/G1o2Ua6dDY/mND2rdZOfY5JFfH+8xh+6+kKjmiHBrb/3hjCgwDjtKW88oaYW35d1QH51dZoeSqkuyJUpr5Rng37q5WdX2TEmIiNSq9943Qb67YrJqmmpC/ZwFAbxN2gd5isXS5jNvt1sqVKzV79uzWMC+1/DIwefJkrVixokcG+tH9RhsyyrXiqwKV7ndqRHKW5k+ZFfL+zShvd55Kako0InmE7ptyX5t56wtmyelcoxHJWZo1Js+gCo1XXrJBjVVVirXFSpJGJI9Q3ok/2B9506WyLVLmMdIU8+2ry5d+ojU7KpQ9NFl5U042upyQy12dy6d6ABACpvxS7LZt21RfX6/s7Ox287Kzs7V161a5XC4DKgMAAABCK+xG6H2xf/9+SVJycnK7ecnJyfJ6vaqsrFR6enqH65eXl2vv3r1tXtu6daskqba2VtXV1QGuuHNuj6d1Gsp+/eHxuFun4VpjuPEcPK6eDo6rx+1unfbm/XnovPJ6vZIs8rjb76tYj1s2SW6PWwdMuK8Odx70Br19+83M4/a0Tjl2vU9dXV2bKUKrtrbW73VMGegPOdzlOYeb99hjj2nRokUdzlu7dq3KysqOuDZfDXJWKlaS01mpL/PzQ9avPyqdMZJsqnQ6lR+mNYYbp9PZOv3hPouJccpq63heb3KM0yGHIuV2uyVFHtwfm9ssc4rTqRS17KuPTLivnE6rJEuvPdaVNZUt08rKXrn9ZlZ58D2M9/3ebe3atUaX0Cvt2LHD73VMGej79u0r6buR+u+rqKiQxWJRUlJSp+vPmzdPM2fObPPa1q1bdeGFF2r8+PEaOXJkQOs9nL0vv6JGFSspqY+GT54csn798e/NRdpTUas+SUmaPHmC0eWYwt69e1VbW6ukpCRN/sFxLSp6XjW1UlJSkiZMCM9jHgp1326Vp7pONlvL21BSUpImnzSqzTKxex+TDu6rH+5HM3j22w1SddXB+ttfItjTLfvPMm3ft119+vTR5InmO3692cp1H2p3+e6D7/scu96mrq5Oa9eu1fjx4xUXF2d0Ob3Opk2b/F7HlIH+6KOPVkxMjAoLC9vNKyws1LBhw2S32ztdPzU1VampqR3Oi4+Pb/NF22CrtFrVKMlmtYa0X39YrbbWabjWGG6sVmvr9If7zHowwFptvXt/uqw2efTdp2lWWwc/AwfPPZtJz73DnQe9QW/ffjOz2qytU45d7xUXF8fxN0B8fLzf65jyS7E2m03nnXeeli9frpqamtbXd+zYofz8fF188cUGVgcAAACETliO0L/xxhuqq6trDesbN27Ua6+9JkmaNm2aYmNjtWjRIo0bN04zZszQ/PnzWx8slZKSoltvvdXI8gEAAICQCctAP3fuXJWUlLT++9VXX9Wrr74qSSouLtaQIUOUlZWl9957T3fccYcuvfRS2Ww2nXHGGVq8eHGPfEosDLRrvfTBA1KDH986LxshKaHlwUh509vO618i2SWVbWg/rzcpmyUpQ15XoxQbpabt5Sr/vx98AavxfKl5urQjUVq6wZAyj8Tc3W79RLFy7Har3IT1H6lrK6ZqZtOpcpQ6VF4S2O2PiLYq4cwMRQ12BLRdADCjsAz027dv92m5nJwcvf3228EtBvjgAenr1X6ulCYpQWqokko+bDsrMVGyR0qu6vbzepGIxjMkZch78Mq/5uZoNdZndLxwvaTiqpDVFijDJEk2ySU1mrD+IzVU6ZJabh/cWBWc7U+55tigtAsAZhKWgR4IK4dG5qMTpbTRvq1Tlig1HFrn1Lbz7CWS6iV7gpTpY3s9UEL9Jml/siIijpFkU0REg6JiOrhVV4RVShosRZlvJHbj7irVuNxy2G0alZ5odDkht7miSDVNNXJEOjQiOStg7TbtrpXX5VFzgydgbQKAmRHoAV+ljZZy/+Xbsnl5UknJwXUebDuvYJbkXCOlZUvTXgh8nSYRJSlFUuTnWyRnnSKHpCr1olOMLiugblr6idYUV+uk9GS9PKf33bbyjtUPad2edRrbf6zyplwWsHbLl27olZ94AEBnTHmXGwAAAAAtCPQAAACAiRHoAQAAABMj0AMAAAAmRqAHAAAATIxADwAAAJgYgR4AAAAwMe5DDwCACezeulmfLntJja76oPe1d3tx0PsAEDgEegAATODTZS/pm4LPQtpnlD0mpP0B6B4CPQAAJnBoZD46Nk79hgwNen9R9hidfMmPg94PgCNHoAcAwET6DRmqy+++z+gyAIQRvhQLAAAAmBiBHgAAADAxAj0AAABgYgR6AAAAwMQI9AAAAICJEegBAAAAE+O2lfDZvl21WrGkIChtR9mtGjttqPoPTQhK+wAAAD0VgR5dirJbJUmN9W6VbnEGta/pNxwf1PYBAAB6GgI9ujR2WssTCRtdnqC0v29XrRrr3UFrHwAAoCcj0KNL/YcmBHXkfMWSgqCP/AMAAPRUfCkWAAAAMDECPQAAAGBiBHoAAADAxAj0AAAAgIkR6AEAAAATI9ADAAAAJkagBwAAAEyM+9AbzJVSo/1zm1SRUqh9BbOMLscQscNqlDHQragYm9YXOIwup73+JVJiomQvkXw8Rqn9dyshsUF2e7TWF7zVZl5t7cZgVAkAAHopAr3BKo7fqYYMr6RquZxrjC7HEDZHyx9JcjoNLaVjdkn2SEn1ko/HyG5v+SNJTueODpexWeMCUx8AAOjVCPQGa470SJIiGqxK6D/W4GqMsW9njRrrW0boUwaH4Qh92QbJVS3ZE6S0bN9WKdstl6tlhD4tLb3dfJs1TkOH3hToSgEAQC9EoA8TURVxypn6gtFlGGJFfoFKtzg1YHiSci4YY3Q57eVNl0o+lDJHS9N8O0Z5eXkqKSlRZmampk/LDXKBAACgN+NLsQAAAICJEegBAAAAEyPQAwAAACZGoAcAAABMjEAPAAAAmBh3uQGAINu4u1qXL/3E6DJCbntUtRQR2O2Pj7bprgarIgPSGgD0DAR6AAiS+OiWt9gal1triisMrib0YjKaZIuTquubtGZH4Lb/QsXqRNm0cXeVbupFvyhllVYrQdLG0sD/ghgfbdNNZw7XCYOTAtougNAg0IeJxuQ6rS+YZXQZhogdVqOMgS0PllpfEIYPlupfIiUmSvYSycdjlNp/txISWx4stb7grSAXaG61tVdIylBtzSatL7i7W23YrHEaMvRGJSYcH9jijtBNZw6XJNU2uA2uxBjboyJ1QFJCTKSyhyYfcXs//KWo5Rel6iNu1ywGuNxKkFQdxF8Qn7pmXFDaBRBcBHqDRTRZJUnN0R45nWsMrsYYNkfLH0lyOg0tpWN2SfZISfWSj8fIbm/5I0lO546gldYTuDVNskhud/UR/wwcf/xfA1RVYJwwOKlXB6Tc1Qlat0calZ6gvCknH3F7X+x06pF3tsix3SW5JIfdppPSj/wXBbNIqLRJLinBbtNJAfgF6ZCNu6tV43L32l88gZ6AQG+w5P8OlqfKKUuKQ/aRI40uxxD7dtaosb5lhD5lcBiO0JdtkFzVkj1BSsv2bZWy3XK5Wkbo09LSg1ygudlqEySPZLMlKCn+JL/Xr63dKLe7Rm5PXRCqQzg59AtS+dINaiyu0qj0RL08x7efyZ7g5UX/1C6nNGpAghbNOfJfkA65fOknvfKSMKAnIdAbzL7Pob5vRCp23GhlXvmM0eUYYkV+gUq3ODVgeJJyLhhjdDnt5U2XSj6UMkdL017wbZW8PJWUlCgzM1PTp+UGuUBzi/98i+SsU7xjpHJO9G3/ft/6glm99tMtAAAkAn3YcBUVqeSq2UaXYQiXfYpkTTu4D/5kdDntlZVKrr6SvVR617dj1DAgXYqJUUMvPq6+apg6U0ofrIZNRSp58Pd+r++asklKl1ybNqnkod6xryPi4pQyb65isnvP6DQAoHMEeoNFxMVJkppranTgs88MrsYYzSf8SEpKU3NNtQ58Ea77IFqSS9rpW32eMyZLMTHy9OLj6ivPqVOkdHV7XzWf0iSl986focFPPG50CQCAMECgN1jKvLmSpOa63nv9b4Q9oWXqSFDsuDD8AmFZoeSqkuyJUtpon1axOhyt07DcpjBypPsqwlEoqVoRDodix/l2fMzMVVSk5pqaXv2eAQBoi0BvsJjs7F4/ylawpEDa4pQ9K0uZt4bhrTvzpksl30iZI6Vc377nEJ2XJ5WUKDorS5n33x/kAs0t+uA19NEjs5Q5y//vkewrmCWXc43sI0f2iu+hlFw1u9d9EgEAOLwIowsAAAAA0H0EegAAAMDECPQAAACAiRHoAQAAABMj0AMAAAAmRqAHAAAATIxADwAAAJgYgR4AAAAwMQI9AAAAYGIEegAAAMDECPQAAACAiRHoAQAAABOzGV0AAARCbe1GrS+YZXQZQeeasknNpzQpwlGofSHc3rqmWpXW7laz19Pm9SZZVejJ0H5vQrt1NldsDlV5ANCrEegBmJrNGidJcrtr5HSuMbiaEEg/+EfVcoV4e/tZJFnav+5scOrNffZO14uLjAteUQAAAj0Acxsy9EZJkttTZ3AloeHatEnNNTWKcDhkHzkyZP1urihSTWONrBFWxdpiJUl9LLWKsnjUL9qhsf2P73C9uMg4/fz4n4esTgDojQj0CBv7dtVqxZKCbq8fZbdq7LSh6j+0/Uf/MFbh3kIt3bBUdU3tQ/fmyEukiEHaXLFZuavvPYJeoo9gXfO44oVIZWyL1I6jI/XkTaHb5s0VdtU0NWls/7HKOytPkrS+YJaczjUakZylWWPyQlYLAKAtAj0MF2W3SpIa690q3eI84vam39DxSCGMs3TDUr2/6/0O59WkniPZpZrGGq0rXxfiysxnepNbklTTVKN1e0K/v7h8BgDCD4Eehhs7bagkqdHl6WLJzu3bVavGevcRtYHgOTQy74h0aETyiDbz1kQ6VCHJEeXQ2P5jDajOXByRX0mqkSPSobH9jw1p31w+AwDhiUAPw/UfmnDEo+orlhQEZHQfwTUieYTyprS9NOOiz7foE2ddy7wTuWyjKyXPz9aBbZ8pKzmr3b4EAPRO3IceAAAAMDECPQAAAGBiBHoAAADAxAj0AAAAgIkR6AEAAAATI9ADAAAAJhZWgf7dd9/Vtddeq6ysLMXFxWngwIG64IILtH79+nbLFhQU6KyzzlJ8fLySkpJ08cUX65tvvjGgagAAAMA4YRXoH3/8cW3fvl2/+MUvtGrVKj388MMqLy/XhAkT9O6777YuV1RUpEmTJqmxsVGvvPKK/va3v+nrr7/WxIkTtXfvXgO3AAAAAAitsHqw1J///Gelpqa2eW3KlCkaNmyY7r33Xp1xxhmSpAULFig6OlorV65UQkKCJCknJ0fDhw/X4sWLdf/994e8dgAAAMAIYTVC/8MwL0nx8fEaNWqUdu7cKUlyu91auXKlLrnkktYwL0mZmZmaPHmyVqxYEbJ6AQAAAKOF1Qh9R6qqqlRQUNA6Or9t2zbV19crOzu73bLZ2dl666235HK5ZLfbQ11q2NpTXK11q4rV6PIYXUrQ7NtVa3QJAAAAhgj7QH/DDTeorq5Ov/nNbyRJ+/fvlyQlJye3WzY5OVler1eVlZVKT0/vtM3y8vJ219pv3bpVklRbW6vq6upAlR8WPvm/Lfp2U5XRZYSExeYN+PGL9bhlk+T2uHXAx7Y9Hk/rtKedT91xuP3hcXtap+yrrrkP7kt3GJxbHre7dRrKWjwed+vU6H0QSsH6WeH9Cj9UV1fXZorQqq31f5AyrAP9b3/7Wz3//PN65JFHlJOT02aexWLpdL3DzZOkxx57TIsWLepw3tq1a1VWVuZ/sWFs754YSTZZbF5FJvTcUfoIq9SUVKr8/F0BbfcUp1MpkpxOpz7Kz/dpHafT2TrN93GdnqyyprJlWlnZbn84Y/tLNvvBfbXZiPJMZZCzUrGSnM5KfWnwuRUT45TVFvrz/BinQw5Fyul0am0v+vmqPPi+Uhng/e10WiVZeL9CO2vXrjW6hF5px44dfq8TtoF+0aJFuueee/SHP/xBN954Y+vrffv2lfTdSP33VVRUyGKxKCkp6bBtz5s3TzNnzmzz2tatW3XhhRdq/PjxGjly5JFvQBj59+Yi7amoVWqGQ+f8PMvockwndu9jUq2UlJSkyZMn+7TO3r17VVtb69c6Pdmy/yzT9n3b1adPH02e2HZ/PLW5VKp1teyrk0YZVKF57H35FTWqWElJfTTc4HOrqOh51Rz82ZgwIXS11H27VZ7quoM/X2ND1q/RVq77ULvLd6tPgN9Xnv12g1RddXB/tr+cFb1PXV2d1q5dq/HjxysuLs7ocnqdTZs2+b1OWAb6RYsWaeHChVq4cKHuvPPONvOOPvpoxcTEqLCwsN16hYWFGjZsWJfXz6empnb4BVyp5Uu43/+ybU9gtdpapz1t20Li4P6z+bH/rFZr65R9fvj9YbXtOThlX/mi0mpVoyRbGJxbVputdRrKWlxWmzzqfe9pVpu1dRrI7eb9Cp2Ji4vjnDBAfHy83+uE1V1uJOn3v/+9Fi5cqLvuukt33313u/k2m03nnXeeli9frpqamtbXd+zYofz8fF188cWhLBcAAAAwVFiN0C9ZskQLFizQlClTNH36dH366adt5k+YMEFSywj+uHHjNGPGDM2fP18ul0sLFixQSkqKbr31ViNKBwAAAAwRVoH+9ddflyStXr1aq1evbjff6/VKkrKysvTee+/pjjvu0KWXXiqbzaYzzjhDixcvVr9+/UJaMwAAAGCksAr07733ns/L5uTk6O233w5eMQACoqC6Tuts58mZeo7WRDp00edb2sz/qrbeoMoAAOgZwirQA+h5Htq+R3utR0lWqULSJ86O72scf/CLeQAAwD8EegBBVXvwoTWW5jr1UaVGJI9ot0y81apbhqSFujQAAHoEAj2AkLA17tBJln8r78Q8o0sBAKBHIdAD6DWqqv+r7cWPyu0Jr8eZexvcclc2SM3eLpdtPu+AvNPSVGkt1e6V54egus7VRxdLVqmptFbln20IWb9Nu/1/LDoA9GQEegC9xvbiR7Vv/7tGl9Gxwz8P7zuxQa2ie+psaiyuCnm3EdF87wIAJAI9gF7k0Mi8zeZQfPwog6v5TlNprZpdbinCooiow4fU5gMH5PV4ZLFaFRFrfLqPaI5Rav1lihqaGNp+o61KODMjpH0CQLgi0APodeLjRylnzAtGl9Gq/LMNaiyuUtTQRKXOyT7ssiVXzdaBzz5T7Lhxynz2mRBVCAAIZxFGFwAAAACg+wj0AAAAgIkR6AEAAAAT4xp69BiNO2tU/c4ONTd4Attw2SypYbq0I1Fa6tut+Zr2t3z5sml3ncp9XKenahrULMVKR7kG6dqiqSovMW5/NA2slWJDf5vFrnAbRgDAkSDQo8eofmeHXEUVQWj54J006iX5eGu+5ii3FCE1u9yG3M4vnDSnxEixNsU1x2hoVboaq4zbH8193VJs+B4XbsMIAOgOAj16jEMj8xa7VZHp8YFruKxQclVJ9kQpbbRPq0Tst0mNUoTdpqiBob2dX7iJsDdLkuoi6lWcuFsjkrMMrMXWOg31bRa7wm0YAQDdRaBHjxOZHt/lrf/8kvdrqeRDKeNUKfdfvtWQt14q2afI9Dil5gawFhOK/HyL5KzTN/Zd+tuJ/1belMsMq2VnQbzklCIHxCt1Ru8+LgCAnoMvxQIAAAAmRqAHAAAATIxADwAAAJhYt6+hLy4u1qpVq/TRRx/p22+/VX19vVJSUjRq1CidccYZOvvssxUZGRnIWgEAAAD8gN8j9O+9956mTJmi4cOH66abbtJ//vMf1dbWKjIyUsXFxXriiSc0Y8YMDRo0SAsWLFB1dXUw6gYAAAAgPwP9RRddpHPOOUdRUVF68cUXtWfPHu3cuVPr16/XRx99pE2bNqmqqkrr16/XnDlz9Nxzz2n48OF6++23g1U/AAAA0Kv5dcmNw+FQUVGRjjrqqE6XsVqtOvHEE3XiiSdq4cKFevbZZ/Xtt98ecaE91Z7iaq1bVaxGV4Cfbvo9+3bxFEoAAICeyq9A/8wzz/jVeEREhK6++mq/1ult1q0q1vbC/SHpK8rOUygB9Exf7HTqkXe2qLbBbXQpQZNVWq0ESRtLq3X50k8C1u7G3VwaC5gdD5Yy2KGR+agYm1IGBfDppj8QZbdq7PShQWsfAIz0yDtb9E5RudFlBNUAl1sJkqpdbq0prgh4+/HRRALArI74p/err75SSUmJXC5Xu3kXX3zxkTbfa6QMitdFt44xugwAMKVDI/MOu02j0hMMriY4EiptkktKsNt00tDkgLYdH23T/5w5PKBtAgidbgf6bdu26dJLL9WGDRskSV6vt818i8Uijyd414UDAPBDo9IT9PKck40uIyheXvRP7XJKowYkaFEP3UYA3dPtQH/99derrKxMDz30kEaOHKmoqKhA1gUAAADAB90O9GvXrtVf//pXXXHFFYGsBwAAAIAf/H6w1CH9+vVTYmJiIGsBAAAA4KduB/q5c+fqr3/9ayBrAQAAAOCnbl9y86tf/Uq33nqrcnJyNHXqVCUnt/3GvcVi0S9/+csjLhAAAABA57od6NesWaO///3vqqio0Oeff95uPoEeAAAACL5uB/obb7xRKSkp+tvf/sZdbgAAAACDdDvQf/XVV3rppZd0/vnnB7IeAAAAAH7odqDPyMho9zApAG2VlZUpLy/viNqIjo7WaaedpkGDBgWoKgAA0JN0O9DPnz9fixcv1rnnniu73R7ImgDTi46OliQ1NDSopKQkIG3OmjUrIO0AAICepduBvqCgQN9++62OPvpoTZ48ucO73Dz88MNHXCBgRqeddpqklkB/JMrKytTQ0HDE7QAAgJ6r24H+0Ucfbf37Cy+80G4+gR692aBBgwIyop6XlxewEX4AANAzdTvQNzc3B7IOAAAAAN3Q7SfFAgAAADCeX4G+rq6uW510dz0AAAAAh+dXoB86dKgeeughVVdX+7T8Z599pvPPP18PPvhgt4oDAAAAcHh+XUO/ePFi/eY3v9Fdd92l8847T5MnT9aYMWOUmpoqu92uiooKbdu2TZ9++qn++c9/auPGjbrssst07bXXBqt+AAAAoFfzK9DPnj1bM2fO1NNPP60nnnhCr7zyiiwWS5tlvF6vYmJidOmll+rpp59WTk5OQAsGAAAA8B2/73ITExOjuXPnau7cufr222/18ccfq7S0VPX19UpJSVFWVpZOOukkRUZGBqNeAAAAAN/T7dtWStLAgQM1c+bMQNUCBETT7lqVL90QuAbLZkkN06UdiVIg2/VB0/6WL5Q37a4L7DaFUNOgZilWGtIwQAMbU4wuBwCAHueIAj0QTiKirZIkr8ujxuKqALac0TKplxTQdrvWHOWWIqRmlzvA2xQ6zSkxUqxNDk+cJpWcYHQ5AAD0OAR69BgJZ7YE7+YGT2AbLiuUXFWSPVFKGx3YtrsQsd8mNUoRdpuiBiaGtO9AibB/9xC6aDeX4gEAEGgEevQYUYMdSrnm2MA3nPdrqeRDKeNUKfdfgW//MCLz1ksl+xSZHqfU3OyQ9h0okZ9vkZw8iwIAgGDhSbEAAACAiRHoAQAAABPjkhvABMrKypSXl2d0Gd1SljZMsju031Kr9Y3btTXI2xEdHa3TTjtNgwYNCmo/AACEiyMK9DU1NXrjjTdUUlKi+vr6NvMsFot++9vfHlFxQG8XHR0tSWpoaFBJSYnB1XSPK2mwZHeo0eKW03tAzhBtx6xZs0LSDwAARut2oF+zZo2mT5+uioqKDucT6IEjd9ppp0lqCfRmZbe3/FIS5bUpyRKrxIx+QeurrKxMDQ0Npt5fAAD4q9uB/pe//KUGDhyo1atXKzs7W1FRUYGsC4CkQYMGmX6k+f8+36JiZ536euOVEzVEE3MvC1pfeXl5pv0kAwCA7up2oC8sLNQLL7ygsWPHBrIewDQad9ao+p0dgb/vfQ9z6EmxAAAgOLod6Pv1C97H5oAZVL+zQ66iji85w3cOPSlWkhpsTQZXAwBAz9PtQH/TTTfpiSee0IwZM2SxWAJZE2AKh0bmLXarItPjDa4mfB16UmyNtU7vZX6hs3SlwRUBANCz+BXoH3zwwTb/3rRpk0488URNnz5dffv2bTPPYrHol7/85ZFXCIS5yPR4pc4x51NcQ+HQk2K3R5dqmH2f0eUAANDj+BXob7vttg5f37BhQ7vXCPQAAABA8PkV6IuLi4NVBwAAAIBu8CvQZ2ZmBqsOAAAAAN0Q0d0VrVar1q5d2+G89evXy2q1drsoAAAAAL7pdqD3er2dzmtububONwAAAEAIdDvQS+o0tK9fv16JiYlH0jQAAAAAH/h1Df3DDz+shx9+WFJLmL/wwgsVHR3dZpn6+nqVl5fr0ksvDVyVAAAAADrkV6BPTU3VscceK0navn27jjrqKCUlJbVZJjo6WqNHj9YvfvGLgBUJAAAAoGN+Bfof//jH+vGPfyxJmjx5sh5//HFlZWUFpTAAPc/mis3KXZ0btPb7V/SXXfZO+znbVqT+EdLmiiI9GsQ6gumKiiJlSCqqKNJCk25DMGyPqlZMRpO2R0Uqd3WCX+vGRcZpTvYcje43OkjVAUBw+RXovy8/Pz+QdQDoBWqaarRuz7qgtX9a02myy95pPxP6udTfLtU01mjd3uDVEUzTm9ySgr8vTSdCssVJBySt29O9Jh4989GAlgQAodLtQC9J1dXV+vOf/6z8/Hzt379fffv21eTJkzV37tx2l+IA6N0So5M0tv/YoPbh2OeQXJIj0tFhXw7bfyVVyRHl0Nj+xwe1lmBxRH4lqebgNh5rdDlhY+PualXXNykhJlKj0n0fod9csVk1TTWqa6oLYnUAEFzdDvTFxcWaPHmyduzYoczMTKWlpWnLli16++239cQTTyg/P19HHXVUIGsFYGLD+wxT3ol5Qe0jb3eeSmpKNCJ5hO6bcl+7+esLZsnpXKMRyVmaNSa4tQRLyfOzdWDbZ8pKzlLeFHNuQzBcvvQTrdlRoeyhycqbcrLP6+WuzuWTDgCm1+3bVv7iF7+Qy+XSRx99pOLiYn3yyScqLi7Whx9+qIaGBt18880BLBMAAABAR7o9Qv/uu+/q4Ycf1skntx0J+dGPfqR77rmHQA8ACLmNu6t1+dJPfF5+e1S1FOH/ekbIKq1WgqSNpeFfK0IrPtqmm84crhMGJxldCgzS7UAfHR2twYMHdzgvIyOj3f3pAQAIlvjolv/OalxurSmu8Hm9mIwm2eKk6vomrdnh+3pGGOByK0FStZ/biN7jqWvGGV0CDNLtQH/BBRfo1Vdf1TnnnNNu3quvvqoZM2b43eYXX3yh3/zmNyosLNTevXsVExOjESNG6IYbbtBPfvKTNssWFBTo9ttv16effiqbzaYzzjhDixcv5rp9AOiFbjpzuCSptsHt13rboyJ1QFJCTKSyhyYHobLASai0SS4pwW7TSWFeK0Jn4+5q1bjcfp/76Fm6HehnzZqln/70p5o5c6ZmzZqltLQ0lZWV6fnnn9e6dev01FNPqaCgoHX5MWPGdNmm0+nU4MGD9eMf/1gDBw5UXV2dnn/+eV111VXavn277rrrLklSUVGRJk2apBNOOEGvvPKKXC6XFixYoIkTJ+qLL75Qv379urtZAAATOmFwUrdGJ3NXJ2jdHmlUeoJfX6Y1wsuL/qldTmnUgAQtmhPetSJ0Ll/6CZ/YoPuB/tDI/M6dO7V8+fLW171eb5v5Xq9XFotFHo+nyzYnTZqkSZMmtXltxowZKi4u1l/+8pfWQL9gwQJFR0dr5cqVSkhouT1ZTk6Ohg8frsWLF+v+++/v7mYBAAAAptLtQJ+XF7rbpaWkpKi8vFyS5Ha7tXLlSs2ePbs1zEtSZmamJk+erBUrVhDoAQAA0Gt0O9BfffXVgayjjebmZjU3N6uyslKvvvqq3nzzTT36aMsT/LZt26b6+nplZ2e3Wy87O1tvvfWWXC6X7HZ70OoDAAAAwsURPSn2kM2bN2vfvn064YQTFBcXd8TtzZs3T0uXLpUkRUVF6f/9v/+nOXPmSJL2798vSUpObv+FoOTkZHm9XlVWVio9Pb3T9svLy7V37942r23dulWSVFtbq+rq6iPeBl95PO7WaSj7he9iPW7ZJLk9bh343jHi2PnG4/a0ToO9nw5d2ldWVqYnn3yy3fy09N2KiZF2797d4fxwExUVpZNOOqnN+5n74Da6PcHfn73BoXPGY4L9GcqfJZhHMM7hurq6NlOEVm1trd/rHFGgf+aZZ3TnnXdq9+7dkqTPPvtMY8aM0WWXXaazzz5bP/vZz7rV7p133qnrrrtO5eXlev3113XjjTeqrq5Ot912W+syFoul0/UPN0+SHnvsMS1atKjDeWvXrlVZWVm36u6OSmeMJJvKd9To5fs+DVm/8F3igZmK9ExXkzNWVd98d4zGN0epr6wq31GttRy7TpXnJEnJUdpQeUDnvFUU1L6akoepue+QTufHR86VTY3y2iPkyQjIeEbQ/bVkgywlX7X+2zL7Imn2RfJ6I+TO/0LN7sN/GmmxSNaYZkWYY3MV7W3W9IYqDW1uDEl/lTWVLdPKSuXn54ekz+6qdDpbp+FeK0LH6bRKssgZhPNi7dq1AW0PvtmxY4ff63T7Lf7VV1/VNddcoxkzZmjq1Km64YYbWueNGTNGr7zySrcDfUZGhjIyMiRJ06ZNkyT9+te/1tVXX62+fftK+m6k/vsqKipksViUlJR02PbnzZunmTNntnlt69atuvDCCzV+/HiNHDmyW3V3x7vbt+jbiip53RY1Vpjkf9xeZq8O3gq1SZLru9e98RbJppZjV8ux60xkfcvUFRmh7clRQe6tq/ZTWiYWHcFzssNID72yMCUlRdcOSwtJX8v+s0zb921Xnz59NHni5JD02V0r132o3eW71ScpSZMnh3etCJ1nv90gVVcpKSlJkye3vxy5O+rq6rR27VqNHz8+IFdewD+bNm3ye51up5A//vGPys3N1VNPPSWPx9Mm0I8cOVKPPPJId5tuZ/z48XriiSf0zTffKCcnRzExMSosLGy3XGFhoYYNG9bl9fOpqalKTU3tcF58fHybL9sG28nnD9c6W7EaXV3fBQgGKSuUXFWSPVFKG936cpSzXmpqVlSMTQPSHQYWGN4urpFWVXnlCoMAbbHWyxK1VxY1G11KlzxyS16vZLHI+r23am+zRxZboyyWZnk9cfLUD+m0jaYGj7zNXkXF2JQyOD4EVR+Zr2rrVe1ulstiCdn7sNVqbZ2G8r2/O6w2a+s03GtF6ATzHI6Li+NcM0B8vP/v190O9Js2ber0bjLJyckdjqB3V35+viIiInTUUUfJZrPpvPPO0/Lly/XAAw/I4WgJUjt27FB+fr5++ctfBqzfUOg/NEHTbzje6DJwOHm/lUo+lDJPlXK/+zJ4+dINaiyuUsogh0bNCcyoSE9lrp/K8JCXl6eSkhJlZmYqNze39fWSq2Zrxykfq/EYr5L6nqScMS902saKJQUq3eLUgOFJuuj84aEo+4hc9PkWfeLkml0A8Fe3A31sbKyqqqo6nPftt9+qT58+frd5/fXXKyEhQePHj1f//v21b98+vfrqq3r55Zf1q1/9qvWBUYsWLdK4ceM0Y8YMzZ8/v/XBUikpKbr11lu7u0kAAACA6XT7Q/BTTjlFjz76aOuDpL7v6aefbveAKF+cfPLJWrt2rW644QadddZZuu6661RWVqZnn31WDzzwQOtyWVlZeu+99xQZGalLL71U11xzjYYNG6YPPviAp8QCAACgV+n2CP2CBQt06qmnavz48Zo1a5YsFouWL1+uu+++Wx988EG3vhmdm5vb5qPlw8nJydHbb7/tdx8AAABAT9LtEfqxY8fqjTfeUG1trW699VZ5vV7de++9+vrrr7Vq1Sodd9xxgawTAAAAQAeO6F57kydP1qZNm7Rt2zbt2bNHKSkpOuaYYwJVGwAAAIAuBOTm2UcffbSOPvroQDQFmE7T7lqVL91gdBnoYZr2t9ztpWl3XZvzKyJluiISd0raoabSWpV/1vm5N8JZr6HxVkU5601xjjYNapZipabSOpWvDU2911ZM1cymU+Uodai8JHh9RkRblXBmhqIGc4tbAIHXrUC/d+9eLV26VB988IFKS0slSQMGDNDkyZN1/fXXtz78CejJIqJb7v3rdXnUWNzxHZ+A7mqOcksRUrPL3eb8skQPkCUypmWZH8z7oQRJskVITc2mOEebU2KkWNvB7aoPSZ9DlS4pXZLU2Mmd2wIp5Zpjg94HgN7H70D/zjvv6JJLLlF1dbWsVqtSUlLk9Xq1efNmvf3221q8eLFWrFih0047LRj1AmEj4cyWpxk3N/BQMARexH6b1ChF2G2KGpjY+npD0SZ5m1rCboTdpqihiZ01oX27atRY7255sNSg8B8ZjrA3H5wefrsCaXNFkWqaauSIdGhEclZQ+mjaXSuvy8N7BYCg8SvQ7927V5dffrkSExP15JNPatq0aYqNjZUkHThwQCtXrtRtt92mSy+9VJs2bWKkHj1a1GAHo20Imsi89VLJPkWmxyk197sHl5VctVjNVTul/lLkgHilzuj8oWYfLSlQ6V6nBqSb4+FnkZ9vkZx1ihwQp9TpoXkQ1h2rH9K6Pes0tv9Y5U25LCh9HHoIHQAEi193uXnqqafk8Xj00Ucf6dJLL20N81LLg6Yuu+wyffjhh2pqatJTTz0V8GIBAAAAtOVXoP/3v/+ta6+9VoMGDep0mYyMDOXm5mr16tVHXBwAAACAw/Mr0G/atEmnnnpql8tNnDhRmzZt6nZRAAAAAHzjV6B3Op1KTU3tcrnU1FQ5nc7u1gQAAADAR34F+oaGBkVGRna5nM1mU2NjY7eLAgAAAOAbv29buXnzZtlsh1+tqKio2wUBAAAA8J3fgf6aa67pchmv1yuLxdKdegAAAAD4wa9An5eXF6w6AAAAAHSDX4H+6quvDlYdAAAAALrB70tuAKC327Vrlz744AM1NDQErY+ysrKgtQ0A6FkI9ADgpw8++EBff/11SPqKjo4OST8AAPMi0AOAnw6NzEdHRystLS1o/URHR+v000/vdL5r0yaVPDS78/n2KZI1Ta6iIpVc9acgVBhYDVNnSumD1bCpSCUP/j4kfV5RUaTpTW45Ir9SyfOd78sjEZEyXZboAWoo2qSSqxZ3ux1XU3XLdFORSq4KTq04vIi4OKXMm6uY7GyjSwHaINADQDelpaUpNzc3pH1GxMW1/r25pkYHPvus02WbT/iRlJSm5ppqHfii8+XChefUKVK65OliuwIpo/VvNTqwLTh9xpx6mmzRLdtVfwTb1Xz0ACk+puW4b9scwArhr8FPPG50CUAbBHoAMJGUeXO1c+NnkioV4XAodtzoTpeNsCe0TB0Jih03LkQVdp/V4Widhqreoooi1TTVyBHpUFZyVlD6iAjQdkU0VUte98HjHv7Hs6dxFRWpuaZGzXV1RpcCtEOgBwATicnOVrT7GNU718g+cqQyr3ym02ULlhRIW5yyZ2Up89ZZIayye6I/3yI56xQ9MkuZszrfrkBauDpX6/as09j+xypvSnBuzVy+dIMai6sUnTVSqXMu73Y79kXzpY1fyj4yS5l33xfACuGLkqtmh+yTI8BfEUYXAAAAAKD7CPQAAACAiRHoAQAAABMj0AMAAAAmRqAHAAAATIy73AAAYCJ7txfr5UXzjS6j13E1Vav56AGKaKpuueNQGIiyxygu4lhJMUaXAoMR6AEAMIEoe0toazhQp10bvzS4ml4qPkbyuqUw2v8D+lZKCWcbXQYMRqAHAMAEJlxyhSSp0VVvcCW9k2tTy4OlIhwO2UcG5yFk/ti7vVgNB+pk9TQaXQrCAIEeANDrba7YrNzVuUFp+9qKqRqqdG2uKNIdqx86ssaOD0xNRomLjNOc7Dka3a/zJxyHq5KrZuvAts2KHTcuLB7s9fKi+XxSg1YEegBArxUXGSdJqmmq0bo964LSx8ymUyWlB7UPs3n0zEeNLgHoUQj0ANDD7dtVqxVLCowuo0v7hlskh0X7dtZqxbvBrTfKbtXYaUM1J3uOJKmuqS5ofTlKHS3TSIfG9h8btH7C3eaKzappqgnqvgZ6KwI9APRQUXarJKmx3q3SLU5ji/FB4yCH5Ig8WG+Nz+s1u3fL7Vojr9e/a4lL/hupvgPjNFHxkuL9rNZ3yQ1RLdOaKE1Z0z9o/YS7oRVO1TTGyhHl1MtrwuMuMf4It7vc7N1ebHQJCCMEesBXZYVS3nSjq0A4KBshKcG4c6J/iWSXVLbhsP2PbU6Tkiao0RMVutq6I8IqJQ5WVEzLf0lRMTYNGJ7k8+plX69UY9M3fndbXyXtqvJ7Nb8NTxstxfRVQ12ddm3rvdc8x0mKk11Sk3aVmXQ/hOFdbjzWMP/5RkgQ6IGuRB8cuWuokko+NLYWhIk0SQnGnROJiZI9UnJVH7b//pKm218LXV1HYsgUPTP4Hm1x1illcLwuOn+4z6u+vOgVHXBK0bFx6jdkaJfL79tZq4Z6t6JjbEoZHLyR+UOiG+KkZik6Lk6DRh0X9P7C1eaKzapprJEjyqERySOMLsdv4XaXG6nlVqavRxwn7Te6EhiNQA905bTbW6YNtcbWgfBRlig1SIpOlNJODX3/9hJJ9ZI9Qco0391C2igrbPnFKAA/X/2GDNXlPtx9ZMWSApVucWrA8CRddOuYI+63K+VLN6ixuEr9Mo/S5XMuDHp/4Sp3da7W7dmqsf0Ha8EU4+8S469wu8vNIS8t/UTaX2F0GTAYgR7oyqAcadbLRleBcJKXJ5WUSGmjpdwHQ99/wSzJuUZKy5amvRD6/gMpbzqffAHAEYowugAAAAAA3UegBwAAAEyMQA8AAACYGIEeAAAAMDECPQAAAGBiBHoAAADAxAj0AAAAgIkR6AEAAAATI9ADAAAAJkagBwAAAEyMQA8AAACYGIEeAAAAMDECPQAAAGBiBHoAAADAxAj0AAAAgIkR6AEAAAATI9ADAAAAJkagBwAAAEyMQA8AAACYGIEeAAAAMDECPQAAAGBiBHoAAADAxAj0AAAAgIkR6AEAAAATI9ADAAAAJkagBwAAAEyMQA8AAACYmM3oAgAAQO+xuWKzclfnGl2G366oKFKGpKKKIi0Mo/q3R1UrJqNJ26Milbs6ISBtejweVdZUatl/lslqtQakTV/FRcZpTvYcje43OqT9mh2BHgAABF1cZJwkqaapRuv2rDO4Gv9Nb3JLCsP6IyRbnHRA0ro9gW16+77tgW3QD4+e+ahhfZsRgR4AAATdnOw5kqS6pjqDK+keR+RXkmrkiHRobP9jjS6n1cbd1aqub1JCTKRGpQdwhL6yUn369AnpCP3mis2qaaox7TliJAI9ABVU1+mh7XtU6/EYXYoplKUNkytpsOz2aP3f51tC3n9t7RVya5pstQmKN6D/gEr7qZR0qWRP1Fe19UZXgyAa3W+0qUddS56frQPbPlNWcpbypuQZXU6ry5d+ojU7KpQ9NFl5U04OSJvV1dXKz8/X5ImTlZAQmF8SfJG7Oje8Pv0wEQI9AD20fY/e2l9tdBnmYXe0/JFU7DRiJClDskjySDKk/wCyHyXZD/7d3SxJig/xNbsAYHYEegCtI/MJtggdGx9jcDXhr6ysTC5Xg+z2aKWlpYW8/9qaTXK7q2WzJSjeMTLk/QdUWaHkqpLsiVLaaMVbrbplSOj3KQCYGYEeQKtj42O04sThRpcR9vLyPlRJSYkyMzOVO3ViyPtfX3C3nM41Soo/STknvhDy/gMq72ap5EMp81Rp6r+MrgYATIn70AMAAAAmRqAHAAAATIxADwAAAJhY2Af6J598UhaLRfHx8e3mFRQU6KyzzlJ8fLySkpJ08cUX65tvvjGgSgAAAMAYYR3ov/32W912220aMGBAu3lFRUWaNGmSGhsb9corr+hvf/ubvv76a02cOFF79+41oFoAAAAg9ML6Ljc///nPddpppyk5OVmvvfZam3kLFixQdHS0Vq5c2frQg5ycHA0fPlyLFy/W/fffb0TJAIAQ2VNcrXWrirVvZ60kad/OWq1YUtDlevt21Qa7NAAIqbAN9M8995zef/99bdy4UXfddVebeW63WytXrtTs2bPbPMEsMzNTkydP1ooVKwj0ANDDrVtVrO2F+9VQ75YkNdS7VbrF6fP6UXYeYAWgZwjLQF9eXq6bb75Z9913nwYNGtRu/rZt21RfX6/s7Ox287Kzs/XWW2/J5XLJbre3mw8A6BkaXS0PRIuIsMgjKTrGpgHDk3xaN8pu1djpQ4NXHACEUFgG+nnz5mnEiBGaO3duh/P3798vSUpOTm43Lzk5WV6vV5WVlUpPT+9w/fLy8nbX2W/dulWSVFtbq+rq6iMpHzAdj9vTOuX875rn4JN1PR5j9pfH7W6dmv14xXrcsklye9w64Oe2eDwt+8EWFSFPo9QnPUZn/myYX22EYv8dqtPjMf/x6s3cB3/u3Qb93HcmGO9HdXV1baahYvR7a7iorfX/ssCwC/TLli3T66+/rs8//1wWi+Wwyx5u/uHmPfbYY1q0aFGH89auXauysjLfigV6CGdsf8lml9PpVH7+ZqPLCXtOp7N1mp+fH/L+Y2KcstqM6z+QTnE6lSLJu3uDGv9ytl/rNpf+TNJR8jQ1SpIqw3R/HON0yKFIOZ1OrQ3D+uCbQc5KxUpyOiv1ZRgdR6fTKskSlPeDtWvXBrS9rlTWVLZMKyvD8mc5VHbs2OH3OmEV6Gtra3XDDTfopptu0oABA1r/02xsbHmzdjqdioyMVN++fSV9N1L/fRUVFbJYLEpKSuq0n3nz5mnmzJltXtu6dasuvPBCjR8/XiNHjgzMBgEm8dTmUqnWpaSkJE0+aZTR5YS9vXv3qra2tmV/TZ4c8v6Lip5XTa2UlJSkCRNC338gxVQ9I9UWKdJzQCm1RX6tG+k5IEmyyS23pD4GHY+u1H27VZ7quoPny1ijy0E37X35FTWqWElJfTQ8jM6zZ7/dIFVXHTy/2l+K3B11dXVau3atxo8fr7i4uIC06Ytl/1mm7fu2q0+fPpo8MXz2caht2rTJ73XCKtDv27dPe/bs0ZIlS7RkyZJ28/v06aMLLrhAr732mmJiYlRYWNhumcLCQg0bNuyw18+npqYqNTW1w3nx8fFtvmgL9AZW256DUyvnvw+sVmvr1Ij9ZbXZWqemP15n3CnZbFJDN+48U9n2v7BwPX9dVps8kqzWHnC8erFKq1WNkmwG/dx3JpjvR3FxcSHdVqPfW8NFR89e6kpYBfq0tLQOP2K577779P777+uNN95QSkqKbDabzjvvPC1fvlwPPPCAHA6HpJaPKPLz8/XLX/4y1KUDQPDtWi998ICU9JVkl1S2QcqbbnRVAACDhVWgt9vtmjRpUrvXn376aVmt1jbzFi1apHHjxmnGjBmaP3++XC6XFixYoJSUFN16662hKxoAQuWDB6SvV0vZiZI9UnJVSyUfGl2VcZrPNboCAAgLYRXo/ZGVlaX33ntPd9xxhy699FLZbDadccYZWrx4sfr162d0eQAQeIcuS4k4eP90e4KUOdq4eoxWlyg1SbLZpQajiwEA45gi0D/99NN6+umn272ek5Ojt99+O/QFAYCRouIk1Utp2dK0F4yuxjhLCqQa53e/4ABAL2WKQA8A+J7GOske4ds19NHx0mm3S4NyQlMbACDkCPQAYBbRB+980OyRFOHfNfSzXg5aWQAAYxHoAcAsTru9ZWr/SlJ919fQlxVKDVXduyUkAMA0CPQAYBaDclpG2gtmSc41XV9Dnze9d98FBwB6iQijCwAAAADQfQR6AAAAwMQI9AAAAICJEegBAAAAEyPQAwAAACZGoAcAAABMjEAPAAAAmBiBHgAAADAxAj0AAABgYgR6AAAAwMQI9AAAAICJEegBAAAAEyPQAwAAACZGoAcAAABMjEAPAAAAmBiBHgAAADAxAj0AAABgYgR6AAAAwMQI9AAAAICJEegBAAAAEyPQAwAAACZGoAcAAABMjEAPAAAAmBiBHgAAADAxAj0AAABgYgR6AAAAwMQI9AAAAICJ2YwuAAAQOF/sdOqRd7aotsGtBfurdKykr3ZX6XdLPzG6tIDL3t2gJEkHGt2KkLSxtFqXh+F2zt3t1jBJG3dX6aYwrA++uWZ3tYZI2ri7WreH0XHcuLva6BIQBgj0ANCDPPLOFr1TVC5JqolySxFSjcutNcUVBlcWeEPqo5Qkq9wer6IkVYfpdv5EsZJsB48D4cusLq5vkiRV1zeF5XkWH02k6804+gDQg9Q2uCVJDrtNjiib1Njy95MGJhtcWeAlbGmQaptls1okt5Rgt+mkoeG3nY7dbsl18Dikh1998E3CusiWaUxk2J1n8dE2/c+Zw40uAwYi0ANADzQqPUHHRiVKJdKx6Yl6Ofdko0sKuBVLClS6xanYKJtcDdKoAQlaNCf8trN86QY1FldpVHqiXp6TbXQ56KaSDxN0YFfLz9bLYXieoXfjS7EAAACAiRHoAQAAABMj0AMAAAAmRqAHAAAATIxADwAAAJgYgR4AAAAwMQI9AAAAYGIEegAAAMDECPQAAACAiRHoAQAAABMj0AMAAAAmRqAHAAAATIxADwAAAJgYgR4AAAAwMQI9AAAAYGIEegAAAMDECPQAAACAiRHoAQAAABMj0AMAAAAmRqAHAAAATIxADwAAAJgYgR4AAAAwMQI9AAAAYGIEegAAAMDECPQAAACAiRHoAQAAABMj0AMAAAAmRqAHAAAATIxADwAAAJgYgR4AAAAwMQI9AAAAYGIEegAAAMDECPQAAACAiRHoAQAAABMj0AMAAAAmRqAHAAAATIxADwAAAJgYgR4AAAAwMQI9AAAAYGJhFejfe+89WSyWDv98+umnbZYtKCjQWWedpfj4eCUlJeniiy/WN998Y1DlAAAAgDFsRhfQkXvvvVeTJ09u89pxxx3X+veioiJNmjRJJ5xwgl555RW5XC4tWLBAEydO1BdffKF+/fqFumQAAADAEGEZ6IcPH64JEyZ0On/BggWKjo7WypUrlZCQIEnKycnR8OHDtXjxYt1///2hKhUAAAAwVFhdcuMLt9utlStX6pJLLmkN85KUmZmpyZMna8WKFQZWBwAAAIRWWAb6G264QTabTQkJCTr33HP14Ycfts7btm2b6uvrlZ2d3W697Oxsbd26VS6XK5TlAgAAAIaxeL1er9FFHPL555/r73//uyZNmqS+fftq69at+t///V99/fXX+te//qVzzz1XH3/8sU455RS9+OKLuuKKK9qs/8c//lF33nmnSktLlZ6e3mk/5eXl2rt3b5vXtm7dqgsvvFCffvqpRo4cGZTtA8LVTzaXam2tSw5rhEbGRBldTtgr37tXDQ0NioiwKDIy9PsrKmqPrBEueRWh5ua2/Td7vTr0rm5Vsyxq+YdXllCXGXQt22SR5D34Rwf/HV4sYVgTAN8MPLBPS8+4MqR9btq0SRMmTNCXX36pY4891qd1wuoa+hNPPFEnnnhi678nTpyoiy66SKNHj9btt9+uc889t3WexdL5G+Th5knSY489pkWLFnU4b+3atSorK/OzcsDc6mP6SZGxqvE0a20tn3B1KcbR8scw3+vbalwVANDjxUr5+fkh7XLHjh1+rxNWgb4jSUlJmjFjhp544gnV19erb9++kqT9+/e3W7aiokIWi0VJSUmHbXPevHmaOXNmm9cOjdCPHz+eEXr0Osl1Lj2226k6T7PRpZhCY2Ojqqur1WzQB5wREQ2y2ar03ah0W4fqssgrS/h8CBsUZvnkgVF6wJwGHtinyQaM0Psr7AO9JB26Kshisejoo49WTEyMCgsL2y1XWFioYcOGyW63H7a91NRUpaamdjgvPj6+zZdtgd5gYkKCJqZ3/DMBAOhdqqurlZ+fr8mTJ5OJDBAfH+/3OmH5pdjvq6ys1MqVK3XCCSfIbrfLZrPpvPPO0/Lly1VTU9O63I4dO5Sfn6+LL77YwGoBAACA0AqrEfpZs2YpIyNDY8eOVUpKirZs2aIlS5Zoz549evrpp1uXW7RokcaNG6cZM2Zo/vz5rQ+WSklJ0a233mrcBgAAAAAhFlYj9NnZ2XrzzTd13XXX6ayzztJvfvMbjRo1Sh9//LHOOuus1uWysrL03nvvKTIyUpdeeqmuueYaDRs2TB988AFPiQUAAECvElYj9PPnz9f8+fN9WjYnJ0dvv/12kCsCAAAAwltYjdADAAAA8A+BHgAAADAxAj0AAABgYgR6AAAAwMQI9AAAAICJEegBAAAAEyPQAwAAACZGoAcAAABMjEAPAAAAmBiBHgAAADAxAj0AAABgYgR6AAAAwMQI9AAAAICJEegBAAAAEyPQAwAAACZGoAcAAABMjEAPAAAAmBiBHgAAADAxAj0AAABgYgR6AAAAwMQI9AAAAICJEegBAAAAEyPQAwAAACZGoAcAAABMjEAPAAAAmBiBHgAAADAxAj0AAABgYgR6AAAAwMQI9AAAAICJEegBAAAAEyPQAwAAACZGoAcAAABMjEAPAAAAmBiBHgAAADAxAj0AAABgYjajCwgXDQ0NkqStW7caXAkAAIBxamtrtWPHDm3atEnx8fFGl9PrHMqih7KpLwj0B+3cuVOSdOGFFxpbCAAAAHq9nTt3asyYMT4ta/F6vd4g12MKTqdT77//vgYPHqycnBx9+eWXRpcEkzvuuOM4j8IUx+bIsP9a9OT9YPZtM1P94Vjr1q1bdeGFF+of//iHhg0bZnQ5vU5DQ4N27typ008/XUlJST6tQ6DvgMViEbsFR4rzKHxxbI4M+69FT94PZt82M9UfjrV+9dVXrb9oHHvssUaXAx/wpdgO3H333UaXgB6A8yh8cWyODPuvRU/eD2bfNjPVb6ZaEb4YoQcAAEArRujNhxF6AAAAwMQI9AAAAGjVr18/3X333erXr5/RpcBHXHIDAAAAmBgj9AAAAICJEegBAAAAEyPQAwAAACZGoA+xhoYGXXvttcrIyFBCQoImTJigjz/+2OiyYCKcQ+GN44NA4DwC4A8CfYi53W4NGTJEH374oZxOp+bOnavzzz9fBw4cMLo0mATnUHjj+CAQOI8A+IO73ISB5ORk5efn6/jjjze6FJgU51B44/ggEDiPAHTGFCP07777rq699lplZWUpLi5OAwcO1AUXXKD169cHve+amhrdfvvtOuecc9SvXz9ZLBYtXLiww2Vra2t18803a8CAAbLb7TrhhBP00ksvHbb9zZs3q76+XkcffXQQqschX3zxhaZPn66MjAzFxMQoOTlZJ598sp577rmg98055L8nn3xSFotF8fHxQe+L49PzfPjhh5o2bZr69OmjmJgYDR8+XL///e+D2ifnEXqbmTNnqn///kpISFB2drZWrlxpdEm9mikC/eOPP67t27frF7/4hVatWqWHH35Y5eXlmjBhgt59992g9r1//3795S9/UUNDgy688MLDLnvxxRfr73//u+6++2698cYbGjdunH784x/rhRde6HD5+vp6zZ49W3fddVdIgktv5nQ6NXjwYN17771atWqVnnnmGQ0ZMkRXXXWV7rnnnqD2zTnkn2+//Va33XabBgwYEJL+OD49ywsvvKDTTz9diYmJeuaZZ7Rq1SrdcccdCvaH0ZxH6G0WLlyonTt3qrq6Wk8++aSuvPJK7d+/3+iyei+vCezZs6fdazU1Nd7+/ft7zzzzzE7Xq6qq8q5atarT+a+//rq3trb2sH03Nzd7m5ubvV6v17t3716vJO/dd9/dbrl//etfXkneF154oc3rZ599tnfAgAFet9vd5vXGxkbv9OnTvbNnz25tH6F30kkneQcPHtzpfM6h0JsxY4b3vPPO81599dXeuLi4wy7L8cH37dq1yxsXF+edO3euX+txHgFH5rPPPvNGR0d7//vf/xpdSq9lihH61NTUdq/Fx8dr1KhR2rlzZ6frPfnkk5oxY4b+/ve/t5u3dOlSnX/++XrmmWcO27fFYpHFYumyxhUrVig+Pl4zZ85s83pubq5KS0u1Zs2a1team5s1e/ZsWa1WPfXUUz61j+BISUmRzWbrdD7nUGg999xzev/99/XYY4/5tDzHB9/35JNPqq6uTnfccYff63EeoScL1iVhV155pex2u8aNG6ezzjpLo0ePDuJW4LCM/o2iu5xOpzcxMdF70UUXHXa5W265xWuxWLxPPPFE62t/+tOfvJK8d955p199Hm7UZcKECd5x48a1e/3LL7/0SvIuXbq09bXrrrvOe/rpp3vr6+v96h9HzuPxeJuamrzl5eXeP//5z16bzdbm3OgI51Bo7Nmzx9u3b1/vn//8Z6/X6/VphN7r5fjgO2eccYY3OTnZu3r1au/xxx/vtVqt3n79+nnnzJnjraqqOuy6nEfoyYqLi72JiYne0047zXvdddd1en56vS2fFiUlJXmfeOIJ77vvvtu6/PPPP9/h8k1NTd4333zT++CDDwZxC9AV0wb6K6+80muz2bzr1q3rctk777zTK8n7pz/9yXvfffd5JXl/97vf+d3n4d6khw8f7j333HPbvV5aWuqV5L333nu9Xq/Xu337dq8kr91u98bFxbX++eCDD/yuB/6bM2eOV5JXkjcqKsr72GOP+bQe51DwXXLJJd4f/ehHrZcV+BrovV6OD1qMGDHCa7fbvQ6Hw3vvvfd68/PzvQ888IA3JibGe8opp3R5yQrnEXqqYF0S9n3Tp0/3/utf/wpo3fBd59cahLHf/va3ev755/XII48oJyeny+X/8Ic/yG636+abb5Yk/e///q9uu+22gNd1uI9DD83LzMwM+pez0Lk777xT1113ncrLy/X666/rxhtvVF1dXZfnA+dQcC1btkyvv/66Pv/8825dVsDxgdRyiYrL5dLdd9+t+fPnS5ImTZqkqKgo3XzzzXrnnXd01llndbo+5xF6Kl/fVw93SdisWbO0Zs0a/ehHP+pwXY/Ho61btx5xregeU1xD/32LFi3SPffcoz/84Q+68cYbfV6vqqqq9YR2Op0Br6tv374dfru7oqJCUsv9g2G8jIwMjR07VtOmTdPjjz+u66+/Xr/+9a+1d+/eLtflHAqO2tpa3XDDDbrppps0YMAAOZ1OOZ1ONTY2SmrZ13V1dV22w/FB3759JUnnnntum9enTp0qSSooKOiyDc4j9GZffvmlRo4c2e67ZdnZ2a3zJamsrEzLli1TXV2d3G63XnnlFeXn5+u0004Lec1oYaoR+kWLFmnhwoVauHCh7rzzTp/W8Xq9uummm/T444/rqaeeUkNDg+bNmyeXy6XFixcHrLbRo0frxRdflNvtbvODUFhYKEk67rjjAtYXAmf8+PF64okn9M0336hfv34dLsM5FFz79u3Tnj17tGTJEi1ZsqTd/D59+uiCCy7QP/7xjw7X5/jgkOzsbH366aftXj800h0R0fkYFucR0HL71aOOOqrd64d+0fz+L6N/+tOfdO2110qShg8frhdffFEnnHBCSOpEBwy72MdPv/vd77ySvHfddZfP63g8Hu9Pf/pTr81m87744outrz/99NNeq9XqnTdvnl+3ATvcdWerVq3ySvK+9NJLbV6fMmVKl9edwThXXXWVNyIiwlteXt7hfM6h4Kuvr/fm5+e3+3Puued67Xa7Nz8/31tYWNjhuhwffN+bb77pleT9wx/+0Ob1Bx980CvJ+5///KfD9TiP0Jt09R2PKVOmtHv90Hc8/vjHP4agQnSHKUbolyxZogULFmjKlCmaPn16uxGYCRMmdLjegw8+qGeffVavvvpqmwd9XH311bLb7frJT36ikSNHdnnpzhtvvKG6ujrV1NRIkjZu3KjXXntNkjRt2jTFxsZq6tSpOvvsszV37lxVV1dr2LBhevHFF7V69Wo999xzslqtR7AHcKSuv/56JSQkaPz48erfv7/27dunV199VS+//LJ+9atfdTo6zzkUfHa7XZMmTWr3+tNPPy2r1drhvEM4Pvi+c845R+edd55+97vfqbm5WRMmTNC6deu0aNEizZgxQ6eeemqH63EeAS24JMzEjP6Nwhenn356651JOvrTmQMHDng//vjjTud/+OGH3oaGhi77z8zM7LTv4uLi1uVqamq8//M//+NNS0vzRkVFebOzs9uM9sA4f/vb37wTJ070pqSkeG02mzcpKcl7+umne5999tnDrsc5ZBxf7nLD8cEPHThwwHvHHXd4Bw8e7LXZbN6MjAzvr3/9a6/L5TrsOpxH6C0ON0L/s5/9zBsfH+9tampq8/qLL77oleT96KOPQlQl/GXxevkaPQAAQG+wb98+9evXT3fffXe7h0u98cYbmjZtml566SVdfvnlra9PnTpVGzZs0I4dO/gUKUyZ4pIbAAAAdB+XhPVsjNADAAD0cEOGDFFJSUmH84qLizVkyBBJLbcS/s1vfqNXXnlFFRUVysrK0q9//WtdccUVIawW/iLQAwAAACZmugdLAQAAAPgOgR4AAAAwMQI9AAAAYGIEegAAAMDECPQAAACAiRHoAQAAABMj0AMAAAAmRqAHAAAATIxADwAAAJgYgR4AAAAwMQI9APQw1157rSIjI9XY2NjpMtOmTVNsbKx27NjRrT5+97vfadSoUWpubu5umR3yeDxKTU3VQw895PM6Tz31lAYOHKi6urqA1gIAZkGgB4AeJjs7W263W5s3b+5w/ptvvqk33nhD8+fPV0ZGht/tl5aW6oEHHtDvfvc7RUQE9r+RDz74QHv37tXFF1/s8zpXX3214uLi9MADDwS0FgAwCwI9APQw2dnZkqQvv/yy3Ty3261bbrlFQ4YM0e23396t9h9++GElJSX5Fbp99dprr2ns2LHKzMz0eR2bzaY5c+bo4Ycf1oEDBwJeEwCEOwI9APQwhwv0TzzxhDZu3KglS5bIbrf73XZjY6OeeuopzZo1q83o/MKFC2WxWLRhwwbNnDlTiYmJSk5O1i233NL6acGUKVPkcDg0ZMiQDkfTvV6vVqxYoUsuuaT1tb179+r666/X4MGDFR0drX79+umUU07R22+/3WbdK6+8UtXV1XrppZf83iYAMDub0QUAAAIrJSVFaWlp7QJ9ZWWlFi5cqDPPPLPbo+tr1qzR/v37NXny5A7nX3bZZfrJT36iOXPm6K233tIDDzygpqYmvf3225o3b55uu+02vfDCC7rjjjs0bNiwNnV8/PHH2r17d5tAf9VVV6mgoEB/+MMfdMwxx8jpdKqgoED79+9v029aWpqysrL0r3/9S9dee223tg0AzIpADwA9UHZ2tr766qs2ry1cuFBVVVV6+OGHu93uJ598IkkaM2ZMh/Ovv/563XLLLZKks846S//+97/16KOPavny5broooskSZMmTdLKlSv1/PPPtwn0r732mkaPHq3hw4e3vvbRRx/puuuu089+9rPW1y644IIO+x4zZky7kXsA6A245AYAeqDRo0eruLi49ZryzZs36/HHH9cNN9ygY489ttvtlpaWymKxKCUlpcP5M2bMaPPvkSNHymKxaOrUqa2v2Ww2DRs2TCUlJW2WXb58eZvReUkaP368nn76ad1zzz369NNP1dTU1GltqampKi8vl9vt9nezAMDUCPQA0ANlZ2erublZGzdulCTdcsstSkpK0sKFC9ssV1RUpLPPPlvJycnq06eP5s2bd9h26+vrFRkZKavV2uH85OTkNv+OiopSbGxsu+v1o6Ki5HK5Wv+9du1a7dixo12gf/nll3X11VfrySef1Mknn6zk5GTNnj1bZWVl7fq22+3yer1t2gWA3oBADwA90Pe/GPvmm29q1apVuvfee5WUlNRmuSuvvFK5ubnav3+/SkpKlJube9h2U1JS1NjYGPB7vi9btkzHHHOMjjvuuHb9/elPf9L27dtVUlKiP/7xj1q+fLmuueaadm1UVFQoOjpa8fHxAa0NAMIdgR4AeqCRI0fKZrPpiy++0C233KKcnJwOvyz6zTffqKGhQR6PRwkJCRo3btxh283KypIkbdu2LaD1Llu2rN3o/A9lZGToxhtv1Nlnn62CgoJ287/55huNGjUqoHUBgBkQ6AGgB4qOjtYxxxyjv/zlL9q0aZMeeeSRDh8C9dJLL+npp5/W4MGDdcsttxz2GnWp5QutkvTpp58GrNYvvvhC27Ztaxfoq6qqNGbMGC1evFgrV67U+++/r8WLF2v16tU6++yz2yzb3NystWvXdnr3HQDoybjLDQD0UNnZ2dq4caOuuuoqnXzyyR0uc+655+rcc8/Vrl27dMopp+iss87StGnTOm1z8ODBmjhxov75z3/q+uuvD0idy5YtU2ZmpnJyctq8brfbddJJJ+nZZ5/V9u3b1dTUpIyMDN1xxx3tHor13nvvqaqqSldeeWVAagIAM7F4vV6v0UUAAEJv+fLlOuGEE3TUUUfpv//9r84991x98sknGjp06GHXW7ZsmS6//HKVlJRo4MCBR1zHqFGjNHXqVC1ZsqTbbVx11VX65ptv9NFHHx1xPQBgNgR6AOilfvGLX+jll19WbW2tjjrqKN1zzz06//zzu1zP6/XqRz/6kXJycvToo4+GoNLD27Ztm0aOHKl3331Xp556qtHlAEDIEegBAH778ssv9X//93+aP39+h9fmh1J+fr62bNkSsEuAAMBsCPQAAACAiXGXGwAAAMDECPQAAACAiRHoAQAAABMj0AMAAAAmRqAHAAAATIxADwAAAJgYgR4AAAAwMQI9AAAAYGIEegAAAMDECPQAAACAiRHoAQAAABMj0AMAAAAmRqAHAAAATOz/A8C6W8CJiJIJAAAAAElFTkSuQmCC", - "text/plain": [ - "
" - ] - }, - "metadata": {}, - "output_type": "display_data" - } - ], - "execution_count": null - } - ] + "cells": [ + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Example 4: Simulated site profile\n", + "\n", + "Generate simulated shear-wave velocity profiles." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "tags": [] + }, + "outputs": [], + "source": [ + "import matplotlib.pyplot as plt\n", + "\n", + "import pystrata\n", + "\n", + "%matplotlib inline" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# Increased figure sizes\n", + "plt.rcParams[\"figure.dpi\"] = 120" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Create a simple site profile" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "profile = pystrata.site.Profile(\n", + " [\n", + " pystrata.site.Layer(\n", + " pystrata.site.SoilType(\"Soil-1\", 18.0, None, 0.05), 30, 400\n", + " ),\n", + " pystrata.site.Layer(\n", + " pystrata.site.SoilType(\"Soil-2\", 19.0, None, 0.05), 20, 600\n", + " ),\n", + " pystrata.site.Layer(pystrata.site.SoilType(\"Rock\", 24.0, None, 0.01), 0, 1200),\n", + " ]\n", + ")" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Initialize the variations." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "toro_thickness = pystrata.variation.ToroThicknessVariation()\n", + "toro_velocity = pystrata.variation.ToroVelocityVariation.generic_model(\"USGS B\")" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Create the varied thickness and velocity." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "%pdb" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "count = 10\n", + "# Create realizations of the profile with varied thickness\n", + "varied_thick = [toro_thickness(profile) for _ in range(count)]\n", + "\n", + "# For eaach realization of varied thickness, vary the shear-wave velocity\n", + "varied_vel = [toro_velocity(rt) for rt in varied_thick]" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Create a plot of the varied velocity models." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "fig, ax = plt.subplots()\n", + "\n", + "for profile in varied_vel:\n", + " ax.plot(\n", + " [layer.initial_shear_vel for layer in profile],\n", + " [layer.depth for layer in profile],\n", + " drawstyle=\"steps-pre\",\n", + " )\n", + "\n", + "ax.set(xlabel=\"$V_s$ (m/s)\", xscale=\"log\", ylabel=\"Depth (m)\", ylim=(55, 0))\n", + "\n", + "ax.grid()\n", + "fig.tight_layout();" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Example of varying bedrock depth" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Setup" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# Import to specify distribution of bedrock depth\n", + "from scipy.stats import uniform" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "tags": [] + }, + "outputs": [], + "source": [ + "# Parameters for distribution of bedrock depth\n", + "depth_bedrock = 50\n", + "depth_bedrock_min = 40\n", + "depth_bedrock_max = 65\n", + "uniform_dist_loc = depth_bedrock_min\n", + "uniform_dist_scale = depth_bedrock_max - depth_bedrock_min\n", + "uniform_dist = uniform(loc=uniform_dist_loc, scale=uniform_dist_scale)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "tags": [] + }, + "outputs": [], + "source": [ + "# Parameters for distribution of velocities\n", + "PARAMS = {\n", + " \"USGS B\": { # To be consistent with example above\n", + " \"ln_std\": 0.27,\n", + " \"rho_0\": 0.97,\n", + " \"delta\": 3.8,\n", + " \"rho_200\": 1.00,\n", + " \"h_0\": 0.0,\n", + " \"b\": 0.293,\n", + " }\n", + "}" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Initialize instances of classes for randomizing profiles" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "tags": [] + }, + "outputs": [], + "source": [ + "var_bedrock = pystrata.variation.HalfSpaceDepthVariation(dist=uniform_dist)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "tags": [] + }, + "outputs": [], + "source": [ + "var_thickness = pystrata.variation.ToroThicknessVariation()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "tags": [] + }, + "outputs": [], + "source": [ + "var_velocity = pystrata.variation.ToroVelocityVariation(\n", + " ln_std=PARAMS[\"USGS B\"][\"ln_std\"],\n", + " rho_0=PARAMS[\"USGS B\"][\"rho_0\"],\n", + " delta=PARAMS[\"USGS B\"][\"delta\"],\n", + " rho_200=PARAMS[\"USGS B\"][\"rho_200\"],\n", + " h_0=PARAMS[\"USGS B\"][\"h_0\"],\n", + " b=PARAMS[\"USGS B\"][\"b\"],\n", + " vary_bedrock=True # Enable variation of bedrock depth\n", + ")" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Randomize soil profiles" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "tags": [] + }, + "outputs": [], + "source": [ + "num_sim = 3\n", + "\n", + "# Create realizations of the profile with varied bedrock depths\n", + "varied_depths = [var_bedrock(profile) for _ in range(num_sim)]\n", + "\n", + "# For each realization of varied depth, vary the thicknesses\n", + "varied_thicknesses = [var_thickness(rd) for rd in varied_depths]\n", + "\n", + "# For each realization of varied thicknesses, vary the velocities\n", + "varied_velocities = [var_velocity(rt) for rt in varied_thicknesses]" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Plot randomized profiles" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "tags": [] + }, + "outputs": [], + "source": [ + "fig, ax = plt.subplots()\n", + "\n", + "for profile in varied_velocities:\n", + " ax.plot(\n", + " [layer.initial_shear_vel for layer in profile],\n", + " [layer.depth for layer in profile],\n", + " drawstyle=\"steps-pre\",\n", + " )\n", + "\n", + "ax.invert_yaxis()\n", + "ax.set(xlabel=\"$V_s$ (m/s)\", xscale=\"log\", ylabel=\"Depth (m)\")\n", + "\n", + "ax.grid()\n", + "fig.tight_layout();" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python 3 (ipykernel)", + "language": "python", + "name": "python3" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.11.5" + } + }, + "nbformat": 4, + "nbformat_minor": 4 }