{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Chron.jl coupled eruption/deposition age and age-depth model\n",
"\n",
"This Jupyter notebook demonstrates the `Chron.jl` Julia package for integrated Bayesian eruption age and stratigraphic age-depth modelling, based in part on the work of [Keller, Schoene, and Samperton (2018)](https://doi.org/10.7185/geochemlet.1826), with age-depth modelling capabilities extended for use in [Schoene et al. (2019)](https://doi.org/10.1126/science.aau2422), [Deino et al. (2019a)](https://doi.org/10.1016/j.quascirev.2019.05.009) and [Deino et al. (2019b)](https://doi.org/10.1016/j.palaeo.2019.109258). For more information, see [github.com/brenhinkeller/Chron.jl](https://github.com/brenhinkeller/Chron.jl) and and [doi.org/10.17605/osf.io/TQX3F](https://doi.org/10.17605/osf.io/TQX3F)\n",
"\n",
" \n",
"
If running this notebook as an online Binder notebook and the webpage times out, click the badge at left to relaunch (refreshing will not work). Note that any changes will be lost!
\n", "\n", "Hint: `shift`-`enter` to run a single cell, or from the `Cell` menu select `Run All` to run the whole file. Any code from this notebook can be copied and pasted into the Julia REPL or a `.jl` script.\n", "***\n", "\n", "## Load required Julia packages" ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [], "source": [ "# Load (and install if necessary) the Chron.jl package\n", "try\n", " using Chron\n", "catch\n", " using Pkg\n", " Pkg.add(\"Chron\")\n", " using Chron\n", "end\n", "\n", "using Statistics, StatsBase, DelimitedFiles, SpecialFunctions\n", "using Plots; gr(); default(fmt = :png)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "***\n", "## Enter sample information\n", "First, let's enter some basic information about your samples. How many are there, what are the sample names (needs to be a valid filename, BTW), and what are the stratigraphic heights and uncertainties? Then, we'll enter the ages as .csv files." ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [], "source": [ "nSamples = 5 # The number of samples you have data for\n", "smpl = NewChronAgeData(nSamples)\n", "smpl.Name = (\"KJ08-157\", \"KJ04-75\", \"KJ09-66\", \"KJ04-72\", \"KJ04-70\",) # These have to match the names used above\n", "smpl.Height .= [ -52.0, 44.0, 54.0, 82.0, 93.0,]\n", "smpl.Height_sigma .= [ 3.0, 1.0, 3.0, 3.0, 3.0,]\n", "smpl.Age_Sidedness .= zeros(nSamples) # Sidedness (zeros by default: geochron constraints are two-sided). Use -1 for a maximum age and +1 for a minimum age, 0 for two-sided\n", "smpl.Path = \"MyData/\" # Where are the data files? Must match where you put the csv files below.\n", "smpl.inputSigmaLevel = 2 # i.e., are the data files 1-sigma or 2-sigma. Integer.\n", "\n", "smpl.Age_Unit = \"Ma\" # Unit of measurement for ages and errors in the data files\n", "smpl.Height_Unit = \"cm\"; # Unit of measurement for Height and Height_sigma" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Note that smpl.Height *must* increase with increasing stratigraphic height -- i.e., stratigraphically younger samples must be more positive. For this reason, it is convenient to represent depths below surface as negative numbers.\n", "***\n", "Now let's see what's in our current directory (we'll use `;` to activate Julia's command-line shell mode, followed by a unix command, in this case `ls`" ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Chron1.0Coupled.ipynb\n", "Chron1.0Coupled.jl\n", "Chron1.0DistOnly.ipynb\n", "Chron1.0DistOnly.jl\n", "Chron1.0Radiocarbon.ipynb\n", "Chron1.0Radiocarbon.jl\n", "Chron1.0StratOnly.ipynb\n", "Chron1.0StratOnly.jl\n", "DenverUPbExampleData\n", "EruptionDepositionAgeDemonstration.ipynb\n", "EruptionDepositionAgeDemonstration.jl\n", "MyData\n", "ffsend\n" ] } ], "source": [ ";ls" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Equivalently, we can also run unix commands using the `run()` function" ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Chron1.0Coupled.ipynb\n", "Chron1.0Coupled.jl\n", "Chron1.0DistOnly.ipynb\n", "Chron1.0DistOnly.jl\n", "Chron1.0Radiocarbon.ipynb\n", "Chron1.0Radiocarbon.jl\n", "Chron1.0StratOnly.ipynb\n", "Chron1.0StratOnly.jl\n", "DenverUPbExampleData\n", "EruptionDepositionAgeDemonstration.ipynb\n", "EruptionDepositionAgeDemonstration.jl\n", "MyData\n", "ffsend\n" ] } ], "source": [ "run(`ls`);" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now that we know how to access the command line, let's make a new folder for our example data. This can be called whatever you want, just make sure it matches `smpl.Path` above" ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ "mkdir: MyData/: File exists\n" ] } ], "source": [ ";mkdir MyData/" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now, let's make some files and paste in csv data for each one. For now, I'm pasting in example data from Clyde et al. (2016), [10.1016/j.epsl.2016.07.041](https://doi.org/10.1016/j.epsl.2016.07.041)" ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [], "source": [ "# You can just paste your data in here, in the following two-column format.\n", "# The first column is age and the second column is two-sigma analytical uncertainty.\n", "# You should generally exclude all systematic uncertainties here.\n", "data = [\n", "66.12 0.14\n", "66.115 0.048\n", "66.11 0.1\n", "66.11 0.17\n", "66.096 0.056\n", "66.088 0.081\n", "66.085 0.076\n", "66.073 0.084\n", "66.07 0.11\n", "66.055 0.043\n", "66.05 0.16\n", "65.97 0.12\n", "]\n", "\n", "# Now, let's write this data to a file, delimited by commas (',')\n", "# In this example the filename is KJ08-157.csv, in the folder called MyData\n", "writedlm(\"MyData/KJ08-157.csv\", data, ',')" ] }, { "cell_type": "code", "execution_count": 7, "metadata": {}, "outputs": [], "source": [ "data = [\n", "66.24 0.25\n", "66.232 0.046\n", "66.112 0.085\n", "66.09 0.1\n", "66.04 0.18\n", "66.03 0.12\n", "66.016 0.08\n", "66.003 0.038\n", "65.982 0.071\n", "65.98 0.19\n", "65.977 0.042\n", "65.975 0.066\n", "65.971 0.082\n", "65.963 0.074\n", "65.92 0.12\n", "65.916 0.088\n", "]\n", "writedlm(\"MyData/KJ04-75.csv\",data,',')" ] }, { "cell_type": "code", "execution_count": 8, "metadata": {}, "outputs": [], "source": [ "data = [\n", "66.13 0.15\n", "66.066 0.052\n", "65.999 0.045\n", "65.989 0.057\n", "65.98 0.11\n", "65.961 0.057\n", "65.957 0.091\n", "65.951 0.066\n", "65.95 0.11\n", "65.929 0.059\n", "]\n", "writedlm(\"MyData/KJ09-66.csv\",data,',')" ] }, { "cell_type": "code", "execution_count": 9, "metadata": {}, "outputs": [], "source": [ "data = [\n", "66.11 0.2\n", "66.003 0.063\n", "66.003 0.058\n", "65.98 0.06\n", "65.976 0.089\n", "65.973 0.084\n", "65.97 0.15\n", "65.963 0.055\n", "65.959 0.049\n", "65.94 0.18\n", "65.928 0.066\n", "65.92 0.057\n", "65.91 0.14\n", "]\n", "writedlm(\"MyData/KJ04-72.csv\",data,',')" ] }, { "cell_type": "code", "execution_count": 10, "metadata": {}, "outputs": [], "source": [ "data = [\n", "66.22 0.27\n", "66.06 0.11\n", "65.933 0.066\n", "65.918 0.087\n", "65.92 0.34\n", "65.916 0.067\n", "65.91 0.18\n", "65.892 0.09\n", "65.89 0.063\n", "65.89 0.15\n", "65.88 0.2\n", "65.812 0.069\n", "65.76 0.15\n", "]\n", "writedlm(\"MyData/KJ04-70.csv\",data,',')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Alternatively, you could download .csv files that you have posted somewhere online using the Julia `download()` function, or using the unix command `curl` throught the command-line interface" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "For each sample in `smpl.Name`, we must have a `.csv` file in `smpl.Path` which contains each individual mineral age and uncertainty.\n", "\n", "***\n", "## Configure and run eruption/deposition age model\n", "To learn more about the eruption/deposition age estimation model, see also [Keller, Schoene, and Sameperton (2018)](https://doi.org/10.7185/geochemlet.1826) and the [BayeZirChron demo notebook](http://brenh.in/BayeZirChron). It is important to note that this model (like most if not all others) has no knowledge of open-system behaviour, so *e.g.*, Pb-loss will lead to erroneous results.\n", "\n", "#### Boostrap relative pre-eruptive (or pre-depositional) mineral age distribution" ] }, { "cell_type": "code", "execution_count": 11, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAlgAAAGQCAIAAAD9V4nPAAAABmJLR0QA/wD/AP+gvaeTAAAgAElEQVR4nOzdd1gUV9sG8DMzu/Qq0rvSRBRFAUWl2IK9t9g1SqLp0USNMXmTLzFGY15jYqJRY++CNfbeKwJWBAQEpIj03WXZmfn+2LyEoMKi2/f+XblyLbOzs88OMveeM2fOUDzPEwAAAENFa7oAAAAATUIQAgCAQUMQAgCAQUMQAgCAQUMQAgCAQUMQAgCAQUMQAgCAQUMQAgCAQUMQAgCAQUMQAgCAQdO9IExOTlZ8WjiWZVVaDLwQx3GaLsHgcByH6RLVD0cYjVD6EUb3gjAsLKy6ulrBlUUikUqLgedxHCcWizVdhcGRSCT4/qF+OMKoH8/zSt/tuheEAAAASoQgBAAAg4YgBAAAg4YgBAAAg4YgBAAAg4YgBAAAg4YgBAAAg4YgBAAAgybQdAF6Qiwjv9/nlt3mxCxvKaRsjIilkFgKqV6u1DutaAG+bwAAaCsE4euSsGTlPe6HZK6TA5XQi3Exoypq+DIpKa8hZVJ+xV1u5X1uWWemhwul6UoBAOAFEISvjuXJirvcomQutDn1VywT3OzvqHM0rc08arAnvSeLm36ObWdHLQmnvS0RhwAA2gV9dq+IJ2TaOXbHI25/byah1z8p+LzBnvSd4YIOzamwvbJFSZgNEgBAu6BF+Io+vsymlvFH+gjMFdiFJgyZ146e4Ev1OczyhMwJxvcPAABtgSPyq/jqJnsyj9/fW6EUrOVmTp3oK1iXyq28j3YhAIC2QIuwyZbf4bak82f7C2yNm/xaB1PyVywTdYC1FpLRLfEtBABA8xCETbPhIbckhTvbn3EyfcUttLCkDsUyvf6SWRlRfd0xdgYAQMPQKGmC/dnc3Gvc8T6Mp8VrBViQLZXQSzD5rOxCAW4pDgCgYQhCRZVJyTsXuB09GF9rJTTjOjlQm6IFw47LsiqRhQAAmoQgVNTn19n+7lQXR6V1ZvZypT4MYqacZZGEAAAahCBUyPWnfHwmtzCUUe5mZ7elK2rIagwiBdBN5ubmmi4BlABB2DiWJ3Hn2cVhzCsME20YQ5H1Ucz8G2w2OkgBdBBN4xCqD/BbbNzyO5ylkLzpo5J91cqGer81OkgBADQGQdiIJyLy3S329y6M6i50+KwtXSYlf6aigxQAQAMQhI2YeZGdGcgE2Kjwgj8BTdZEMp9dZR9XoVkIAKBuCMKGHHrM3y7hP1P91KBtm1HvBjJvn2dV/UYAAFAPgvClWJ58cJldEcGYKHms6IvNa0fnVpFt6eggBQBQKwThS+3M4JxMSU9XNc2CJqTJ8ghm3nVOiigEAFAjBOFLLU7hPgtWS2Pwf7o5Uf7W5A9cVggAoEYIwhc7nMPLOKL+SbEXhTH/l8hW1Kj5bQEADBeC8MW+T2I/C6bVf2+Its2oGBf65ztoFAIAqAmC8AWuFvHZlWSkt2Z2zrcd6WW32eJqjbw5gEFYuXLlnDlz5syZ89133+3fv5/jmvDVs6SkJDU19WXPFhcXp6WlKaPGprly5Ur//v1f9qxEImnZsmV1dXX37t2TkpLqPVtTU3Pz5s2XvbZnz543b96cNm3a7t27Gy2j3sf38fF59uyZAuVrGILwBb69xc1uSws0tG+8Lanh3vT3t3ApBYCqbNmy5cGDB7a2tjU1NbNnz548ebLirz158uTbb7/9smf/+uuvDz/8UBk1No1YLM7JyXnZszzPZ2Rk8Dw/atQoBweHes8WFxd36tTpZa8dOXKko6NjQUFBRUVFo2XU+/hvvfWWiYmJAuVrmDpuzFteXn7z5s2HDx/26dPHzc2tdjnLsufOnbt69SrHcVFRUZ07dyaEHDlyJCsrS76Cubn52LFj1VBhXfdK+auF3LYYoZrft64v2jNtdte8H0S7m+POvQAq0bt373feeYcQ0rNnz969e69fv54QkpGRcfjwYUJInz59vL295WseOXIkOTnZ09NzyJAhMpns0qVLhYWFO3fuFAqFgwcPzsjI2L9/f1VVlY+PT0xMzNWrV588ebJz505TU9O+ffvu3r27Z8+eO3bsMDU1HT9+/LVr1y5dulRTU9OlSxf5ES8+Pj4yMjI+Pl4ikQwePNjDw0MikRw8eDAqKmrr1q0CgWDkyJF2dnaEEI7j9u/f/+DBAz8/v0GDBlEURQi5c+fO4cOHHR0dnZycnv+MT548iY+PpyhqyJAh8iWurq7GxsYVFRU7duwoKCiwt7fv16/f6dOnOY7buXOnfLckJib6+PhcvXr1/v37s2fPlr9E/vKEhIS0tLSoqKiwsDBCyIEDB7p162ZtbU0ISUlJEYlEdT9+//79vby8GIYhhJw/f/7KlSsODg5Dhw41Nzd/+vTpjRs3goKCtm/fLj/OW1hYqPK33Qh1tHo8PDw+/PDDDz744M6dO3WXr1279v333y8qKiorK+vfv//SpUsJIb/99tuWLVtu3Lhx48aN5ORkNZRXz/dJ3PtBjKk6viG8lLMZiWtFf30TZwoBVIvjuDNnzvj7+xNCLl68GB4enpeXl5ubGxYWdunSJULIu+++O2fOHJqm161b17179+rq6oKCAolEkpGRkZmZmZGRERERUVFRYW1tffLkyUePHhUVFYlEooyMjKysLJlMNnLkyIEDB2ZlZUkkEolE8tVXX1VXVxNCJk6cuGrVKkLIuHHj+vbtW1hYmJ2dHRYWlpaWVlZWNmrUqKFDh1ZXV9+4caNjx44lJSUsy77xxhvr1q0zNjZetWqVvIVw5syZ6OhokUiUnJz80Ucf1ftoT548CQ0NvX//fnl5+YgRI+QLP/7447S0tMjIyCtXrtjZ2T169Oj06dO5ubmEkIyMjIyMjJqamv/85z99+/Y9dOhQTU0Nx3GzZs2SdwUvWrToyJEjPM8PHjx4+/bthJB33nnn8ePH8i1v3bp169atdT++/NOVl5cvXLhw8uTJHMcdPHgwNDS0srLy3r17kyZNmjRpklQqPXDgQL9+/dTxy24Ar3oSiYTneQ8Pj8OHD9ddXlZWVvt4x44drq6uPM8PGjRo7dq1DWzN2NhYLBYr+Nbl5eVNKjW7kmu+UVpS3aQXqURpNe+wSXq3hNN0IU3GsmxlZaWmqzA4VVVVMplM01U0wbxrskFH1fTfnef+jiIjI83MzGxtbY2MjJo1a3bx4kWe57t377506VL5Cj/88EPPnj1TU1PNzc0LCgp4nmdZNiAgYPv27bt27YqJiZGvtmvXrujo6Lpb3rBhQ79+/eSP5Zm3b9++5z/+lStXWrduzfO8qanpmjVr5As/+uijuLi4/Px8Qsjx48flC4cMGbJw4cKNGzdGRkbKl9TU1Hh7eyclJcXExPz+++/yhZ9++mlwcHDdt5gzZ86UKVPkj+Wn98Risa+v78mTJ42MjCoqKmrXfPLkiVAorP0xOjr6vffeq/0xICDgwoULAwYMGD16tHzJgQMHWrZsyfO8m5tbSkqKfOHcuXM/++yzuh+f53mGYR4+fGhqanr37t3ajS9duvTs2bMmJiaFhYU8z4vFYlNT07y8vOf30gtxHFe3eKVQR8Ontlldj5WVVe1jmqZrV7ty5Up1dbW/v390dLS8+a82S5K5qf60jZE63/PFrI3IrDbMf25y27qr9VpGAPXo40aH2qtpcl1H0xccRpYsWfLOO+9IJJKEhISBAwfevXs3OTl58eLF8mdjYmIWL158+/btVq1ayU+q0TQdGRmZlJQUEhJSu5Ho6OgFCxb4+PgMHDhw9OjR8g7DeuRdoIQQjuMWLFiwe/dugUBA03R2drZ8edeuXeUPunXr9sMPPxBCKIqqu1De31hQUDBy5Ej5QrFYfO/eveTk5OXLl9euduTIkbrvm5ycPHz48Npna5ebm5uPHDnSw8Ojf//+AwcOHDp06PM1v/CUYW1JkZGR6enplZWVz6/zvPT0dGtr61atWsl/jImJSUpK6tixo6enp729PSHExMTEwcGhsLDQ2dlZkQ2qgkZ7AP9HJBItWLBA3rR3dXUViUQpKSmLFi0KCAg4cOCAvIu5Fsdxn3zySd2FzZo1mzNnzgu3LJFIhEJFz/aVSKmND6nEgbxEohXX8U1pQX5Iou8W1bSw1HQpTcFxnEQiqfdbA1WT97to526naVooFNb7UtvViSJE86fATUxMxowZ88knn5w9e9bY2FgqlcqXS6VSExOTuksIIdXV1fW+1tvZ2aWkpFy7di0hISE2Nnbz5s0vfAv5g23btp04ceLq1auWlpZ3794NDQ2tfa/aB/Lt8zxf+1j+QCgURkREzJ8/v3az9vb29Qqu9751n5W3TWtt3LgxLS3twIEDX3zxxaVLl2bPnv2ymuuquzX575RhGJb9e1ifWCx+YZun3j6s/VwCwT/pQ1EUzzf+rYjjOKlUyvO8RCKp+/KGGRkZNXrbSM0HYXV19bBhw9q3bz9jxgxCyK+//ipfXlFRERgYuGvXrlGjRtVdn6IoNze3unvB2tr6ZX//DMMofmjYnkn6uBEXrRmfYsWQ6QHklwfMsnBNl9IUFEU1abeDUjD/o+lCXoCiKDV37TRJYmJiUVGRl5dX165dt2/fLm8Mbdu2rWvXrh06dMjIyLh9+3ZQUFBZWdnhw4e3bt3Ksmxpaan8tRUVFZaWluHh4eHh4fn5+bdu3WrdunXts/Xk5OT4+PhYWlrKt1+7fPfu3UFBQYSQXbt21ba64uPjJ06cKJPJ9uzZM3XqVHd397i4uKVLl9rY2BBCJBIJTdNdunTZvXt3+/bt5a+t93ZdunSJj4+fPn06RVF1r3yQSqXV1dU+Pj4ffvihg4PDH3/8YWVlJZPJqqqqzM3NG9hRe/bsef/99xmG2b17d2hoqLGxsbxrNDg4uLq6+ujRowMGDLC2tq738Vu2bCkQCI4dO9arVy+JRLJnz56XtVsaJT+2yL/wKf5PXZF/exoOwpqamlGjRpmbm69bt65eaFtaWoaGhj5/vQ5FUR999JGCQ3KFQqHiLcL16bKlnRihUIv+Yj8IIoG7av7TUWivAyOQ/8ZxnEwmU3y3g1LI/6lrZxBqp++//37NmjXl5eWFhYVff/11hw4dFi5c2KdPnzt37nAcl5ube+jQIUdHx6VLl/bu3VvePzly5MiYmJiqqiqpVNq+fXtPT89+/fotWbIkNDRUJBLdunXryy+/tLW1nT17docOHXx9fTds2FD3HYcPH/7jjz8OHTq0qqpKHodyN27cGDRoUGlpaXl5+cqVK2tqagQCQXx8/MGDB9PT05s1azZhwgQjI6MJEya0adOmW7dulZWViYmJV69e/fbbb3v27JmcnFxRUfH8qMu33347Pj4+IiLCxcWl7oWSRUVFHh4eXbt2NTc3P3369H//+18zM7MxY8a0adOmRYsWq1evftkeMzMz69q1q5ub29mzZ/fu3UsI+eSTT6ZOnZqQkJCdnS2/IiAyMrL248vD3sTEZOXKlePHj4+MjJR/pRgzZszFixdf4VdGUZT8XGaTDuwKbVmRBqlSeHp6rlq16o033iCE5OTkSKVST0/PcePGlZWVJSQkyBvL8t+WPBHLyspat279008/1Y53kjMxMSktLVUwCOXf1xRZ88ZTfsQJNm2kQAPTyTTonQusoyn5KkRnDnAcx4nF4oa/WoLSiUQiY2NjBKGC8vLyJBIJIcTY2NjZ2bn2W3h1dXVycjJFUW3btjUy+nuwQGFh4b179zw9Pb28vORLeJ7Pz8+vrq728vLKyspKT083NTXt0KGD/CXyZ6VSqYeHx6NHj7y8vGq3X1JSkpSU1Lx584CAgMePH3t7e5uZmT18+LCgoEAqlYaEhBgZGRUUFLi7u0skkitXrhgZGbVr16721yqvRL7Q1NSUEFJVVXXz5k17e3tPT8/i4uK616cRQmQy2c2bNymKCgkJycrK8vb2zs3Ntbe3F4lEt2/f5jguKChIfm2GvLaSkhI3N7enT59aW1vX/gnn5OTY29uXlJSYm5sXFxdnZma2a9dO3jAlhDx+/Dg1NVXeKiWENGvWrPbje3p6Pnr0yNPTk6bpkpKSlJQUR0dH+QBdiUTy9OnT2mqzs7OdnJxqd3jDeJ6vqqpS7uUW6gjCd99998GDB+fPn2/VqpWdnd3KlStXrFiRmZkZExPz7rvvRkZG1n7+rVu3dujQITo62sTE5NChQ8HBwXv27Kn3t62iIJx5kXU2pea317oZBjIq+E57ZY9GC801342tEAShRiAIdZQ8CF1dXWuXyIPw+XN+IKeKIFTHwXXChAnl5eWfffaZ/EcHB4e4uDiJRGJlZSX/dlDL1tY2ISEhJSWF47iJEydGRESooTxCiIQl29O5m0O0MWpaWFKRzvSfqdy7gVoX0gDwmqZMmVLvW6Opqelbb72lqXoMk/q6RpVFFS3CLencpjTurze0MQgJIdeK+FEn2dQRAk3N+tYkaBFqBFqEYCBU0SLUhSOr6q15wE3x095dEWpPeViQ3ZmYaAYAQPm09+ivNpkVfMozfoCHVu+K2W2Zhbc4HWu8AwDoAq0++qvH2lRunA9trN1dSn3dKZYnJ/MQhQAASmboQcjxZP1DfrIW94vKUYR80oZenIx7MwEAKJm2B4CqHcvlnUxJm2ZadvHgi7zpQyc/4++VolEIAKBMhh6Ea1O5Kf66sROMaDLNn15xF0NmAACUSTcyQEWKq8nRHG5UC53ZCW+3Yrakc+VaMSU4AICe0JkMUIVt6Vx/D6246ZKCnM1Idxd6cxoahQBaod5dHUBHGXQQbs/gxrTUsT0wM5D+5Q6uowDQCpgITT/oWAwoUb6Y3C3he7rqwDCZuqKdKYYmZ58gCgEAlMNwg3BnBjfAkzbSwR3wdgD9K4bMAAAoiQ7mgJLsfMSN8NbJjz/Blz6Zx+VWoVEIAKAEOpkEr09H+0XlLIRkdEv6jwdoFAIAKIGBBqHu9ovKzQykV93nahCFAACvTWej4PXobr+oXCsbKsCaSsD9KAAAXpsOh8Er0+l+0VozAzFkBgBACQwxCHc94vp76HC/qNwgTzqjgtwuwZAZAIDXouNp8Ep2ZHAjdWdatZcR0GSKH/XHfTQKAQBei87nQVPpR7+o3FsB9OY0TiTTdB0AALrM4IJQ18eL1uVuToXaY8gMAMBr0YtAaApdHy9az7QAXFAIAPBa9CcSFKFP/aJyAz3oh2UEd+sFAHhlhhWE+tQvKiegyXhfai0ahQAAr0qPMkEB8ZncMC99+8hv+dMb0zgpohAA4JXoWyo04Fk1ufmU7+6iP/2icj5WVGtbam8WkhAA4FUYUBAeeszFuNBmAk3XoQLT/GlcUAgA8GoMKAj3Z/MDPPStOSg31JtOesanl2PIDABAkxlKENZw5Ggu18ddP4PQiCZjW9J/pqJRCADQZIYShOfyeR8rysVMP4OQEDI9gF6bihszAQA0maEE4YFsboCHPn/YABuqhSV16DGSEACgafQ5G+o68Jjvr6cnCGtNC6BXP8BpQgCApjGIILxfyotkpJ2dngfhSG/6UiGXJ0IWAgA0gUEE4f5sfqAHpecxSIipgAz1ojc+RBACADSBgQShnp8grDXFj16byiEJAQAUp//x8KyaJD/jY/RuQpkXCnegTBhyIR9RCACgKP0Pwr8eczHOtAmj6TrUZaIvvQYXFAIAKEz/g3B/tv6PF61roh+9L4srr9F0HQAAOkLPg7CGI8dzuX6GcYJQzs6YxLjQOzLQKAQAUIieJ8SFItrHinIy1XQd6jXFj8YdCgEAFKTnQXg4jzaQ8aJ1xbpRuSJyF7etBwBQgJ6HxLEnjEGdIJSjKTLeB7etBwBQiL4HYU9psL5PKPNCk/3oTWmYgxsAoHF6HoTNjHhDjEFCWlpRrWyoA9lIQgCARuh5EBqyKf70GvSOAgA0BkGot4Z70ZcL+dwqDJkBAGgIglBvmQrIiBb0hjQEIQBAQxCE+myKH736PubgBgBoCIJQn4XaU5ZCcg5zcAMAvJw6gpDn+ezs7KKiouefqqioOHXqVGJiYt2FeXl5x48fz8zMVENtem+SH4bMAAA0ROVBuHDhQltb25YtW86ePbveUykpKX5+ft9///2oUaNGjRrF8zwhZNOmTW3btl22bFlYWNjPP/+s6vL03nhfel8WVyrVdB0AANpK5UE4cODA5OTkWbNmPf/U/Pnz33rrrSNHjty4cePKlSvHjh2rrq7+5JNPtm3btn///hMnTsybN6+kpETVFeo3O2PSyxVzcAMAvJTKg7B169YeHh7PL5dIJAcPHpwwYQIhxNLScvDgwfHx8efPnxcIBD169CCEtGnTxs/P7/Dhw6quUO/hgkIAgAYINPXG+fn5LMvWZqSnp+eJEyceP37s6elJUX/PBuPh4fH48eN6L+R5PiEhwcjIqHaJubl57969X/guHMdxnKFnQE9nEicmt56ybZupY5od7n/U8F5QS77Pa/92QD3wT139eJ5v0m6n6cbbexoLQolEQlGUQPB3AcbGxmKxWCKR1C6RL5RIJPVeyPP85s2bGeafW87b29tHRka+8F2qq6uFQqGya9c9b3rRa+9T34ewangvjuMkEkndXxCogUQi4Xkeu13NcIRRP57n6yVFw0xMTBrNQo0FoZOTE8/zz549s7e3J4QUFRU5Ozs7OzsXFxfXriNfWO+FNE3v2rXLxMREkXdhWdbMzEyJZeuouNZ82F7Z4s7Gxqo/TsrbJdjt6mdsbIwgVDMcYdSP53me55W72zV2HaGNjU1AQMCZM2fkP545c6Zz584dO3ZMS0vLz88nhIjF4mvXrnXq1ElTFeoTL0uqTTNqP+bgBgB4jspbhNeuXTt58uTly5fLy8sXLVoUHh5uY2MTHh5eUVExa9asjz/+WCqVJicn37t3b9y4cdbW1qNHjx4zZszMmTM3btzYtWvXNm3aqLpCAzHFj16byg33xhQKAAD/ovIgrK6uLikpCQ8PJ4SUlJSIxWJ/f/+5c+cyDDN16lRra+v9+/c3b9784sWL1tbWhJBVq1b98ssve/fuDQsL++CDD1RdnuEY6kV/eJl9XMW7m2M8BQDAPyj5Zew6xMTEpLS0VMFzhBUVFZaWlqouSVe8e5F1MKUWtFdto5DjOLFYbG5urtJ3gXpEIhHOEaofjjDqx/N8VVWVhYWFEreJjjIDMj2AXn2fY3Xsmw8AgGopGoRffPFFcnKySksBVWvbjHIyI8dykYQAAP9QNAi3bt0aHBwcERGxbt06sVis0ppAdab503/cx9hRAIB/KBqEd+7c2bFjh7m5+ZQpU5ycnOLi4pKSklRaGajCmJb06SdcngiNQgCAvykahMbGxiNGjDh27Ni9e/feeeed+Pj4du3adezYcdWqVSKRSKUlghJZCMkwb3rDQwQhAMDfmjxYxt/f//vvv8/Ozv74449v3LgRFxfn7u4+f/583CZCV8h7R3HfegAAuSZfR1hVVbVt27ZVq1ZdvXrVy8tr2rRpeXl5y5Yt27ZtW2JiIkYSa79Qe8raiJx+wnd3wQWFAABNaRHeunVrxowZLi4u06dPt7e3P3DgQHp6+rx583755ZeUlJS8vLyjR4+qrlBQoqn+9B+4MRMAACFE8RZhr169jh8/7uDgMGPGjLi4OC8vr7rPenl5eXl5oXdUV4zzob+4UVMkYewVmpYAAECfKRqEVlZWW7ZsGTZsWN0bAdZ19OhR+RxpoP2sjcgAD3pTGvdREGZUAABDp+hxsEWLFmFhYfVSMCMjIy4uTv7Yzc0NJwh1yDR/etV9jJgBAFA4CDdu3FhQUFBvYUFBwapVq5RdEqhDVyeKociFfEQhABi61+oZKyoqsrW1VVYpoGZv+dMrMcsMABi8Rs4RHj9+fOfOnYSQ8vLyH374wdHRsfYpmUx27NixkJAQ1RYIKjPJj/46saZQzDiYaroUAADNaSQIHz58KA/C6urqkydPCgT/rG9vbx8aGrpw4ULVFggqY2NEhnrRf6ZynwVjyAwAGC5F70fo5OQUHx8fERGh6oIahfsRKtGtYn7wMTZ9lIBR3rX1uB+hRuB+hBqBI4z6qeJ+hA21CC9fvrxhw4bBgwf37t170qRJmzZt2rRp0/OrrVixQokFgTq1s6McTcnhHL6fO2aZAQAD1VAQPn78+ODBg0FBQb179z5+/HhRUZHaygK1eacV/dtdtp97kyfbAwDQD4p2jWoPdI0qVzVLPLfVXBoo8LZUTqMQXaMaga5RjcARRv1U0TWq6CiJwsLC/Pz82h9379796aefysfRgE4zZsh4H1xHAQCGS9Eg7N27d+25wBUrVgwfPnz58uUjR4786quvVFUaqMvbreg/UzkJq+k6AAA0QaEglEgkycnJffr0kf+4dOnSIUOGVFVV/fbbbz/++KNYLFZlhaByLa2oEDtq1yM0CgHAECkUhCUlJTzPOzs7E0IePnyYnp4eFxdH0/TYsWMrKyszMzNVWyOo3jut6N/uIQgBwBApFITy20oUFhYSQhISEoyMjLp27Vr7rEQiUVFxoDb9POg8Ebn5VMdGTgEAvD6FgtDMzCwkJOTzzz/fs2fP77//3qNHD/mYwNTUVEKIq6uramsE1WMoMg1TjwKAQVJ0sMyKFStu3749ZMiQ6urqRYsWyRdu3brVx8fHwcFBZeWB+kz1p3c+4p5Va7oOAAD1asJ1hBzH5eXlOTo6CoVC+ZL79+8LBAIfHx+VlfcCuI5QdaacZX2tqbmvN/UoriPUCFxHqBE4wqifuqdYq0sikfz2228JCQlPnjzhuH91oKWnpyuxINCgj9vQsYfZT9rQRpiFGwAMhqJBOHHixB07dkRGRkZHR9M0DpP6KciWCrAmOzO4sT74FQOAoVAoCKVS6Z49e7777ru5c+equiDQrI/aMAtusAhCADAcCh3vSktLpVJpbGysqqsBjevrToll5Fw+rqMAAEOhUBDa29u3bNny7t27qq4GNI4i5L3W9E+3cR0FABgKhYKQoqg1a9Z88803p06d0rm7VUBTTe/ctdUAACAASURBVPSlLxRw6eX4RQOAQVD08gl/f/+srKzq6moTExNTU9O6Tz179kw1tb0YLp9Qg7nXWDFL/tvpVcbi4/IJjcDlExqBI4z6afLyiYEDB5aXlyvxjUGbvdeabrNb9p8QxtpI06UAAKiYokG4ePFildYBWsXFjIp1o9c84D5ug+GjAKDnmnCYO3fu3LBhw/z9/Vu3bi1fsnz58rVr16qmMNCwD4Po5Xc5FicKAUDfKRqEe/bsiYmJefToUVBQUFlZmXyhsbHxggULMHxGL4XaU54WZEcGho8CgJ5TNAg//vjj8ePHX79+/b333qtdGB0dnZubm5ubq5raQMPmBjP/l8hx+J4DAHpNoSAsLCx89OjRe++9R9M0RVG1y+U3YMrPz1dVdaBRb7hRZgJy8DEahQCgzxQKQvnkovXm2iaEyNuCyh3GClplbjv6/xIRhACgzxQKwubNm7ds2XL9+vWEkLotwp9//tnBwcHX11dV1YGmDfGiRTJyMg/dowCgtxS9fOLrr78eO3ZsYWGhn5+fVCrdtm3bjh07EhISfv31V1zDq8coQj4Npr+7xXZ3UfSfCgCAbmnCjXnXr18/Z86c2jOCVlZWX3311UcffaSy2l4MM8uoGcuTgJ2y9VFMhCPV+NqYWUZDMLOMRuAIo36qmFmmCUFICJHJZLdu3SoqKrKysgoJCak315p6IAjV77d73OEcfm8vhQ6yCEKNQBBqBI4w6qexKdZEItHmzZtPnTqVl5fH87yDg0NkZGRQUJBGghDUb4of/e0tWfIzvm0zhRqFAAA6pPEW4e3bt/v165ednU0IsbGxoSiqpKSEEGJnZxcfHx8ZGamOMutAi1AjlqRwN5/yW2Iab3CgRagRaBFqBI4w6qeKFmEjo0bFYvGQIUPKyspWrFjx9OnTkpKSZ8+elZWVbdy40cjIaNiwYYWFhUqsBrTW2wH0iTzuYRmGjwKAvmkkCHfu3Jmenn7w4MF33nnHzs5OvtDKymrcuHGnT5+uqqpavXq1Im8jk8lqampet1jQHAsheS+Q+QbXFAKA3mkkCI8dOxYdHd2lS5fnn/Lz8xs+fPjRo0cb3gLHcTNnzrSzs2vevPnUqVNlMlntU19//XWzOpycnGQy2dixY2uX1M7uDdrgwyD6WC6X9AyNQgDQK40E4b179zp37vyyZzt37nzv3r2Gt7B58+aTJ09mZWXl5uYmJibWbUHOnj07/X/GjBnTpUsXgUBQVVX1zTffyBdevHixSR8GVMpCSD4LZr68gUYhAOiVRoKwrKzM1tb2Zc82b968tLS04S2sX78+Li7OxsbGwsJi5syZ8ulp5ExNTW1tbW1tbS0tLRMSEqZOnSpfbmZmJl9ubW2t8AcBdZgRSCc/4y8VolEIAPqjkSCsqalpYBwawzBSqbThLaSlpbVq1Ur+ODAwMD09/fl1Dhw4QNP0G2+8If9x9uzZ1tbWYWFhhw8ffuE2i4qKCusoLi5uuAZQFiOafN6OnnOV1XQhAABK0/h1hJcuXXrZUPjExMRGX15WVlY7ztXCwkJ+6UU9a9eunTRpkjxx586d6+XlZWZmtm3btiFDhty4cSMwMLDuyizLtm/fvu6Up25ubmfPnn3hu1dWVjZaITTJMGeyOMl4X5ooxvHFcchxnEQieX6KdlApsVgslUpx+YSa4QijfjzPi0QixaeCMTMza/TvovEg3Llz586dOxV8y+c1b9689ka+paWlDg4O9VbIz88/cuTI0qVL5T+Gh4fLH0ybNm3Xrl2HDh2qF4QMw+Tk5Ch4HSEhBFf5KN03odw3SfQAH8ELr67nOE4gEOA6QjVjGAbXEWoEjjBqxvM8TdNqnVlm//79jXZ+Nqx169aJiYl9+/YlhCQmJtZLNULIunXrIiIifHx8nn+tVCoVCDDXs9YZ7k0vSeYSMrmhXore2BkAQGs1EjNt2rR5zTeIi4ubOnVqbGyskZHR0qVLFy9eTAiZM2dOYGDghAkTCCHr16+fN2+efGWxWLx69eoePXqYmJhs3779+vXra9asec0CQOkoQr7uwHxwiR3oQQsQhQCg41Te3urTp8/8+fMnT57Mcdz7778/YsQIUuemhunp6V5eXsOGDZP/SFHU+fPnf/vtN47jAgICTpw40aJFC1VXCK/gDTfKxYxsSuMm+SEJAUC3Ne3uE9oAc41qicuF/MgT7L0RAvN/f5vCXKMagblGNQJHGPXTwFyjAC/TyYGKcqYW3sKlFACg2xCE8OqWhDOr7nOpmIkbAHQZghBenaMpmdWWmX0VlwwCgA5TNAi7du26aNEi3HQJ6vkwiH5Qyv/1GI1CANBVigahl5fX/PnzPTw83nzzzdOnT+vcEBtQESOaLI9gPrjEVuNcIQDoJkWDcNOmTdnZ2f/5z38uX74cExPj7++/aNGioqIilRYHOqGXK9XKhvr5DjpIAUAnNeEcobOz82effZaWlnbw4MHAwEB5A3HcuHHnzp1TXX2gE37qRP+QzOaJ0E8AALqnyYNlaJru27fvqlWrZsyYIZFINm/eHBkZGRoaev78eVXUBzqhpRU1PYCeg1EzAKCDmhaEPM8fP3585MiR7u7ua9asmT59+vXr1/ft2ycQCHr06PHo0SMVVQnab1475mw+fzwXjUIA0DGKTrFWUFCwbt261atXp6WlBQYG/vjjj+PHj5ffOLdDhw6xsbEuLi4XL1709vZWZbWgvcwFZE0k89Y59tZgGhOlA4AOUfSQFRIS8vTp0yFDhqxevToqKqres0KhMCQkRCgUKrs80CU9XKjuztSca9yS9pouBQBAYYoG4bx584YPH+7o6PiyFY4cOaKkkkCHLe3EtI2X9XWmB7TUdCkAAIpR9Bzhhg0bcnNz6y1MTExs2RIHPPiHtRFZEUG/d01QUaPpUgAAFKNoEGZlZUkkknoLRSJRZmamkisCHdfHjXS15+ZdwwX2AKAbXmuu0ZSUlAY6S8FgLe7A7svmzzzBCFIA0AGNnCNcu3btt99+Swh5+vTpiBEj6t4FsLq6Oi8vb/LkyaotEHSQlZBf0YWZfJZNHiqwwAgqANBujQShh4dHz549CSEbN24MDQ2t2/4zMTEJCgqaMGGCagsE3dTPndrpRH10mf2jG24VCwBarZEg7NmzpzwIOY6bNWuWv7+/WqoCffBrFyZ0j2xjGjfeB3f7AgDt1dARateuXS1atPj9998JIUlJSX369GnxIuoqFXSMuYDs6MF8cpm9V4qThQCgvRpqEbq6uvbq1cvLy4sQEhUVVV5erqaiQF8E2VILQ5mRJ9grgwRmmG8GALQSpXN3FjQxMSktLa07bKcBFRUVlpaWqi4J6uI4TiwWm5ub1y6ZeIYV0mQ1ThaqkkgkMjY2ZhjsZLXCEUb9eJ6vqqqysLBQ4jYb+pa+a9euTz/9tNFNZGRkKK8e0EMrujBhe2QbHnITfHGyEAC0TuNdo2orBfSV/GRhzEFZR3sq0IbSdDkAAP/SUBB27ty5c+fOaisF9FhrW2pxODP0GHtxoKCZsaarAQCoQ9FRo2FhYS8cMopRo6Cgib70QE9qyDFZNSZfAwBtglGjoD6Lwpg3T7FvX2D/jMSYDgDQFhg1Ckr2/KjRusQy0v0v2QAPel47DJxRJowa1QgcYdRP3aNG63n8+PHy5ctv3bqVm5vr5OQUFBQ0Y8YMzDUDTWIqIAm9BJ32yjwtyFjMOAMAWkDRI9GFCxeCgoL++9//VlZWBgQESKXSP/74Izg4eO/evSqtD/SPkyn5K5aZdYW9VKhjvREAoJcU6hrled7f39/c3Hzv3r0eHh7yhYWFhWPHjk1MTMzJyVGwo1Ip0DWq5RruGq11OIefclZ2vK8AF1QoBbpGNQJHGPVTRdeoQi3CgoKChw8fLlu2rDYFCSEODg6///57cXHx3bt3lVgQGIhYN2pZZ6bXX5iJFAA0TKFzhBYWFjRN29jY1FsuX2JlZaX8usAAjPCmRTISe5g91Y9pYYl2IQBohkItQgsLi6FDhy5btqze8mXLloWHh/v4+KigMDAIE33pL0Po7gfZzAq0CwFAMxpqET548ODMmTPyx+Hh4f/3f/938+bNwYMHOzs7P3369NChQ9evX583b55a6gS9NcWPrqwhvQ6xZ/ozLmZoFwKAujU0WObPP/+cMmVKo5tQ85WIGCyj5RQcLFPPD8ncn6ncqX4CJ1MV1aXnMFhGI3CEUT91X0c4fPjwqKgoJb4ZwMt82pbmeNJ1v+xILNPSCu1CAFCfhoLQ0tKy7pcdmUx2/fr19PT0qqqquqtNnz5dVdWBIZkTTDubkaiD7P7eTHs7ZCEAqImiM8ukpqb269cvLS3t+acQhKAsE31pGyPyxiHZlhhBT1dkIQCog6Izy0yfPl0gEFy6dGnMmDEff/yx/LJCZ2dnzCwDyjXIk97ZQzD2tGzXI07TtQCAQVAoCFmWvXTp0qJFizp16mRsbGxsbOzj4/P+++8vWbJk5syZLIvb6oAyRTlTh2MFH1ziVt1HFgKAyikUhEVFRVKpNCAggBBibm5eVlYmXz5gwICcnJz79++rsEAwSO3tqDP9mZ9ucx9cYmVIQwBQJYWC0M7OjmGYp0+fEkLc3NwSExPly588eaLC0sCw+VhRVwcJHlWQvkdkJdWargYA9JdCQSgUCjt06HDu3DlCyNChQ69cuRIXF7dy5crRo0e7uLj4+vqquEgwUJZCktCLCWlOhe2V3cWUpACgGoqOGv32228LCwsJIX5+fj/++OPXX3+9atUqT0/PLVu2GBkZqbJCMGgMRb4PZQKsuZiDsvVRglg3DCXVRoVi8rjqX99UzATE14oS4I6ToAte/Q71JSUltra2yq1GEZhZRsu92swyjbpYwI88yb7lT33RnmGQhs9R88wyeSL+7BP+1jM+qZhPesZXs8T739OmV9SQJyI+pDkV4UB1dqQ6O9DN1XevNvXBEUb9VDGzTBOCMC0tbdmyZYmJibV3qH/vvffatm2rxGoUgSDUcioKQkJIgZiMPSUjhGyOEThiJrZ/U08QVtSQ+ExuUxqX+JSPdqHbNaOC7UhwM8rD4gXfTcpryOVC/lIBf6mQu1zIt7OjPgqiB3jQtB59j8ERRv00GYQnT54cMGAAx3Fdu3Z1dHQsLi6+cOGCRCLZtGnTyJEjlVhQoxCEWk51QUgIYXny9U12bSq/KZqJctajA+prU3UQHs3l16Vyhx5zUc70eB+qnwdt0pS3knEkPpNbepsrlpAPg+hJfrS5oqdltBqOMOqnsSDkeb5ly5YODg579+51dHSULywtLR03btyFCxfy8vJMTdX3/RxBqOVUGoRyx3L5iWdkMwOZucF61bx4HaoLwkOP+a9ushKWvNOKHtmCbmb8Wlu7WMAvvc2decLNDKTnBjPGOj5JOI4w6qexIMzPz3d2dj5//nyXLl3qLs/MzPT29r5x40ZISIgSa2oYglDLqSEICSF5Iv7NUyxFyIZoxt0cYaiSIDySw391k62sIV+G0EO9lPmdI6OCn3WFSy3j10YyYfY6/OvDEUb9VBGECg3qsrS0ZBjm+UObfMnzd66vSyaTff311+Hh4b179z579mzdpw4cODCyjvz8fEJIUVHRlClTOnbsOGbMmOzs7KZ9GjAYLmbUib6CN9zojntk29Jxyb2SXSviu+yXfXyZ/TCIThoqGO6t5JZ3C0sqviczvx096KhszjVWgsmpQKMUCkJzc/NRo0YtWbKE4/51xPnhhx+6devWokWLBl67ePHiPXv2/P777+PHjx8wYEBubm7tU6mpqUVFRSP+R57wY8eO5Thuw4YNzs7OgwYNeqUPBQaBocicYPrQG4L/JHLjTrNlUk0XpBdKpWTmRXbQMdn0ADplmGBUCxV2Po9uSScNFaaXk5AE2eVCXCcKGtNQ12jdO9RXVVV99913dnZ2Q4YMcXJykt+h/sGDB/PmzZs7d+7LtsDzvIeHxx9//BEbG0sIGTp0aEhIyPz58+XPLl26NCUl5c8//6xdPzU1tW3btoWFhVZWVizLOjk57dmzp15/LLpGtZx6ukbrEsnI7KvsX4/5PyOZaEMdQaOUrtHNadynV7lBntS3HRnb1zsX2CQ7H3HvX2S/6sDEBejYhYc4wqifum/Me/Hixbi4uLpLnj59+v3339dd0nAQPnv2LCcnJywsTP5jWFjYzZs3665w/vz52NhYV1fXadOmderUKSUlJSAgwMrKihDCMExISEhSUlK9IASox0xAfo1gDj3mx59mh3hRC0MZ/RiRqE4Py/i3L7ClUpLQSwMn7UZ40yF2VN8j7KMKfmEoLhMFdWvogDF27NjBgwe/ztblk9HUnkS0tbUtKCiofTY4OPjzzz93cXG5du1a9+7dDx06VFBQUPeMo62trXwLdbEsGxgYSFH//LG4urr+9ddfLyygqqqq7pqgBhzHSSSSV56o4ZV1syWXY6nPbjJtdrG/hcu62BvWiUORSFRTU/MKLUKekD8eMt/dZj4NZOP8WIYilZWqKLARjjQ51oMafU4w4qjs9/CaJl2boUE4wqgfz/MikUjx9c3MzGi6kZ6GhoLQyMio3vRpWVlZt27dysnJcXZ2DgoK8vPza3jr1tbWhBCRSCRvxlZWVtbNuR49esgf9O7du7i4+Pfff+/fv3/dT1hZWSnfQl0Mw+zbt8/Y+J+OG6FQ+LJmMs/zym1BQ6M4jnvh0Co1sCBkUw9yOIeffp6OdaN+DGcshOqvQjNomn6FrtEnIjLtnOxpNbkwkPG31vBciRaEnOpPppxlB5xh9vUW2OvCTDQ4wqgfz/MURWlg1CghRCqVTpo0ydvbe/Dgwe++++6wYcP8/f0HDRpUXl7ewKscHBzMzMwePnwo/zEtLc3Ly+uFa7q4uJSWlnp5eWVkZNTe4PBl6/v4+PjW8bJtgmGKdaOShgqqWRIcLzueiyEYL7XhIdcuoaaLE31hgMDfWiuaNcYM2RTDRDpRXffLMivwuwM1UTQIZ8+evWnTpvfee+/ChQsZGRlXrlyZP3/+sWPH3nrrrQZeJRAIRo0atXz5ckJIXl7e7t2733zzTULIwoUL09PTr1+/Ls+8x48fr169OiYmpnPnztbW1ps3byaEHD169OnTp/JRNgBNYmNE1kUxv0Qw086zk86wxbiL078VV5Nhx9kfU7ijfQRzg2mtOilHEbIojJkZSPc8xOaJkIWgFrwCpFKpubn5d999V2/5+vXraZouLCxs4LV5eXmhoaHu7u42NjYLFiyQL3RwcDh+/PiAAQPMzMzc3d3NzMxmzJghlUp5nj9z5oyLi0vLli3t7e337dv3/AaNjY3FYrEiZfM8X15eruCaoCwsy1ZWVmq6ir9VSPkPL8mcN0u3pLGarkW1qqqqZDKZImueyOXct9TMuiyTKLS6xnx/i229q6ZI0b91zcARRv04jquoqFDuNhWaWSYvL8/V1TU5OblNmzZ1l5eUlDRr1uzq1auhoaENbyE/P9/S0vL580aVlZWlpaXOzs51z22wLJufn+/g4CAUvuAMDy6f0HLqv3yiUVeL+GnnWDdz8msE42WpTc0f5VHk8okajnxxg92cxv8ZyfR01YH98Pl19kgOf7KfwEpbz/XiCKN+vKZmlrG2tmYY5v79+/WW37t3jxDSvHnzRrfg5OT0wiOjhYWFm5tbvb9ehmFcXV1fmIIAryDMnro+WNDViQ7dK1uYxNUY1njSv6WW8RH7ZPdLSeIQgU6kICHk245MZweq/xGZSKbpUkCvKTqzTO/evd99991Dhw7VtiAvXLgwadKk4OBgb29vVVYIoARCmswNpq8NElws4NrFy848MayTT6sfcF33y6b603t6Mbp1X8CfI5gWltTQ4zKpQX59AfVQ9DZMWVlZPXv2TEtLs7KycnZ2LioqevbsmZOT05EjR9R8S0J0jWo5LewarWdPFvfhJS7KmVocxjjoy30NX9Y1WlxNpp9jH1Xwm2OYVja60RCsh+XJqJOskCZbYrRqWA8hOMJogsa6Rgkhnp6et27dWrVq1eDBg729vWNjY3/++ed79+6p/8a8AK9psCd9Z7jAyZQE7a756bY+95Qey+XbxctaWpHLgwQ6moKEEIYim6OZzAr+20T9/VWBRinUIiwpKRk5cuSCBQu6deumhpoahhahltP+FmGtB2X8h5fY7Ery385MLx05bfYy9VqEEpZ8fp3dmcGvi2K6u+j2R5N7IiLhe2XLI+hBnlo0HymOMOqnsRYhy7LHjx/H6BXQM/7W1KFYwfdh9IwL7NDj7CN9uYL7ahHfIUH2uJLcGirQjxQkhDibkfhezPTz7O0SPfk1gfZQKAibN28eGBh49epVVVcDoH4DPOjbwwSh9lToHtmsK+wzXb76XsKSOdfYQUdlC0LoHT2Y17ybvLbp2Jz6qRMz+BhmSAAlU7STYc2aNT/99NOaNWuenwUbQNcZM2RuMH1/hFBAk8BdNYuSOF28VezVIhKSILtTwt8YIhjVQov6D5XozZb0CG9qyDEMIgVlUnTUqJOTU90bR9Sl4BaUBecItZwOnSN8oful/JxrXNIz/tuO9GhV3pZWiSpqyPwrkh1ZzLLOzEg9jcBaHE8GHZN5WVDLIzR/iwocYdRPFecIFb1v2/z586uqqpT4xgDaKcCG2tOLOZfPf3aVXXiLWxBCD/PS3jjkCdmSxn12jevuSN0cRDtb6HkKEkJoimyOEYTukW1L50a31P/PC2qgUBBmZmb6+fk5ODi0bt0aQ2bAEHRzoi4OFJzP57+4wf7nJvdpW3qsj3ZNTk0ISXrGv3eRrZKRHd2ZdpZSY2NDuR+xlZDs6sH0+EvWvjmlJffNAJ3WyPcpqVQ6bNgwb2/vN954o3379r6+vrdv31ZPZQAa19WJOtVPsCScWXGPa58g25ymLRcd5lbxMy6wsYdk433oa4MEEY4GFwZtmlFfhjBjT7HVOng2F7RNI0H4888/x8fHDxw48Jdffpk1a1ZRUdHkyZPVUxmAloh1oy4PFPwQxqx7yHlvl313i3sq0Vgxj6v4dy+ybeNlZgJyd7hwWoD2dtuq2sxA2tOC+uwakhBeVyN9KceOHevRo8fevXvlP/r7+0+bNq2kpMTW1lb1tQFokVg3KtZNcLuEX3ab89tZM8KbnhlIt22mvhTKrOC/T+J2PuLe8qfvDRfqzeRwr2NNJBOSIOvuzA3UpqvsQec08q8nKyure/futT/KH2dmZqq0JgCtFWRL/dGNeTBC6G5ODTzKBsfLFidzOVUqHDhdw5H92dyw42zHPTI7E/JghHCRHk2R+ppsjMiWGGb6eTa7ElfZw6trpEUokUhMTf/5mzMzMyOEiMVi1RYFoN3sTcj89vTn7enz+fzmNK59AtvGlnrTh+7tSnlYKK2NmPyMX5fKbUnnfK2pib70n1FCrb0tnwZ1cqA+acO8eYo93U8gQLMQXknjw8xSU1OPHz8uf/zs2TNCyPXr10UiUe0KPXv2VFFxANqMIqSbE9XNiVnWmTmUw+3I4OdfZy2FVA8XqocrFeNMv8INj7Ir+bP5/Nl8/uwTXsKSCb7U+QECHytDPQ2omFlt6ZN53P/dYr8K0fyVhaCLGrmg3svLKysrq+FN4IJ6qEvXL6h/HTwht5/xJ/L4k3n8mSectRHlZ038rCk/a8rfmnIwJSYMMRUQQoiFgGJ5kivi86r4rEqSJ+KzKsmlQl4i4yOd6W5OVJQT1dqWUnwgjCJ3qNdjT0QkJKFmX29BqL1avzTgCKN+Grig/tdff63b+AOABlCEtGlGtWlGfRhEOJ7JruRTy0lqGX+/lD/0mCuuJmIZkU/eVinjaUJczSlXc8rDnLiYU/3cyZchNK6KezXOZmR5BDP+NHtziMDMUC6nBKVRdIo17YEWoZYz5BahBhl4i1Bu/GnWxoioc+o1HGHUT5M35gUA0HK/RDD7s/lDj3Xsyz1oHIIQAPSEtRFZG8lMP6/b99IC9UMQAoD+6O5CDfOm3r2I6WagCRCEAKBXFnZkbhXz2zO0Y1pY0AUIQgDQK6YCsj6K+eASW4CZP0AxCEIA0Deh9tQUP3r6eXSQgkIQhACgh74MYTLK0UEKCkEQAoAeMmbI+ijmQ3SQggIQhACgn0KaUxN86Th0kEJjEIQAoLe+7sCklfO7HqGDFBqCIAQAvWXMkDXdmHcvsoXoIIWXQxACgD4Ld6DG+9AfXEYHKbwUghAA9NzXHZibT/m9WegghRdDEAKAnjMVkHVRzMyLXAnmIIUXQRACgP7r7EAN8aRmXUEHKbwAghAADML3YcyZfP5IDm7SBPUhCAHAIJgLyKquzNsX2IoaTZcCWgZBCACGorsLFeNMzbuGDlL4FwQhABiQnzoxe7P4c/noIIV/IAgBwIBYG5HlEfS0c6xYpulSQGsgCAHAsAzypNvZUV/dRAcp/A1BCAAG55cIZmMad7EAHaRACIIQAAxQcxPyUyfmrXOsBM1CQBACgGEa1YJuZUN9dwtJCAhCADBUv3ZhVt7nbjxFB6mhQxACgIFyMiWLw5ipZ9kaTMdt2BCEAGC4JvjSXpbUoiQkoUFDEAKAQfslgv75DnunBB2khgtBCAAGzc2c+rYjM+kMK0Oz0FAJ1PAeT548Wb16dXFx8cCBA7t37167XCwWHz58+OrVqzzPR0dHx8bGEkLi4+MfPnwoX8HS0nLGjBlqqBAADNm0ADohi1uYxH3RHm0DQ6Ty33pFRUWnTp1ycnJ8fHxGjx69e/fu2qc2btz4008/WVpaNmvWbMqUKV999RUhZMOGDWfPni0pKSkpKSktLVV1eQAAhJDV3Zhf7rIYQWqYVN4i3LRpk4eHx8qVKwkhVlZWCxcuHDZsmPypiRMnTp8+Xf44ICAgLi5OnoXDhw+fPHmyqgsDAKjlYkYtDmMmnGZvDBGYMJquBtRL5S3Cc+fO9ezZU/64Z8+eN27cqKqqkv9obGxc2CYq8AAAIABJREFUu5pIJLKyspI/PnXq1MKFC+Pj42UyTIsLAGoywZduZUv9XyIusTc4Km8R5ufnd+3aVf7YwcGBoqi8vDxfX9+665SVlX3++efz588nhPj7+9fU1FRVVX355ZdLliw5ffq0kZFR3ZU5jpsyZQrD/POdzc7O7rvvvnvhu4vF4rprghpwHCcWiymK0nQhhkUkErEsi3/tr2lpCNXpL0FvR2lHO4X6SHGEUT+e50UiEU0r2oozMTFpdGWVB6FQKKxt2MlkMp7n6zYECSEikWjQoEGxsbHy7tBFixbJl3/++eetW7fesWPHuHHj6q5PUVRkZKRQKKxdYmlpWW+btaRS6cueAhXhOI7jOOx2NWNZ1tjYGAfl1+RqTH6J4OOuCK8NoEwVODriCKN+PM/LZDLFd7siX8pVHoQuLi65ubnyxzk5OQzDODk51T4rFosHDBjQokWL5cuX13uhqalpu3btMjMz6y2nKGrSpEkmJiaKvDvDMDg0qBlFUdjt6sf8j6YL0XlDvcn2R+xXt8iS8MZ3Jva5+vE8r/TdrvJzhAMHDty7d291dTUhZMeOHbGxsUZGRg8fPkxOTpZKpSNGjGjevPkff/whb7qyLCuRSOQvLCwsvHDhQps2bVRdIQBAXSu6MNsz+FNPMILUUKi8RThgwIBff/21S5cuvr6+x48fP3z4MCFk5cqVmZmZ3bp1O3jwYLt27cLDw+UrHzp0KCAgoEuXLiYmJqdOnerTp8/AgQNVXSEAQF12xmRtJDPxNHtrqKAZOj4NAMXzKv/Ww7LsmTNnysrKunbtam9vTwjJzc2trq42NzfPycmpu2ZISEhmZmZKSgrHcf7+/q1atXp+ayYmJqWlpQp2jVZUVFhaWirlU4CC5INlzM3NNV2IYRGJRDhHqFwfXWazK8nung3tUhxh1I/n+aqqKgsLCyVuUx1BqFwIQi2HINQIBKHSVbMkfK/sozb0RN+XnkLCEUb9VBGEmE8IAOAFjBmypTsz+wr7sEzHWgvQVAhCAIAXC7ShPm/HjD2NGxbqOQQhAMBLvR9EO5iQ724hCfUZghAA4KUoQlZHClbeZ8/lo4NUbyEIAQAa4mRKNkYLRp9k88WaLgVUA0EIANCIHi7UZD9q7CkZi2ahPkIQAgA07usODE3hZKF+QhACADSOpsimaMHK+9yxXLQK9Q2CEABAIY6mZFM0M+kMmydCFuoVBCEAgKKinam4VvSbp1icLNQnCEIAgCaY3442psnca7iRvf5AEAIANAFNka3dBQmZ/NpUDJzREwhCAICmaWZM9vdm5l1jrxfjEKoP8FsEAGiyABtqXZRg3EWjnCqcLdR5CEIAgFcR60ZN95ENOsaKZJouBV4PghAA4BV93ErWyoaadAZjSHUbghAA4NX90ZXJrOS/T8LAGR2GIAQAeHWmArKnF7PqPrf+IbJQVwk0XQAAgG5zMaOO9WGiDrBWQjLEC60L3YPfGQDA6/Kxovb0YuLOs+dx20IdhCAEAFCCUHtqa3fBsBOy5GfIQh2DIAQAUI4eLtR/OzH9j7BZlchCXYIgBABQmjEt6Y/a0H0OswW4nb3uQBACACjTR0H0mJZ0zEHZE5GmSwHFIAgBAJTsi/b0ZD86+qAME7DpBAQhAIDyzW5LTw+gux1gH1UgC7UdriMEAFCJT9rQ5gLS/S/2RF+mhSWl6XLgpRCEAACq8nYrmuVJj7/YY30YHytkoZZCEAIAqNDMQNpUQLrtl+3qKejiiCzURjhHCACgWlP86E0xgqHHZZvSMB+pNkIQAgCoXA8X6mRfwYIb3Fc3cc8mrYMgBABQh9a21NVBguO5/OQzrBQtQ22CIAQAUJPmJuRoH0F5DenxFy4x1CIIQgAA9TETkN09meFedNhe2YFsZKFWQBACAKgVRcgHQfTeXoIPLrEfXEI3qeYhCAEANCDUnro5RPBERLrsk6WXo2moSQhCAADNsDYi23sw43zoiP2y3+5xHNJQQxCEAAAaI+8mPd9fsPsRF7pXdvMpwlADEIQAABrma00d6yv4oDXd74jsg0tslUzTBRkYBCEAgOZRhEzwpROHCgslpF28LCETHaXqgyAEANAWTqZkawzzaxfmm0Su017ZiTykoTogCAEAtEtvV+rGEMGstvSMC2zX/bJz+YhD1UIQAgBoHYqQEd70nWGCSX702FNs/yOy47k88lBFEIQAAFpKQJO3/OmHIwVDvelPrrBtd8v+uM+JMZRG2RCEAABazZghU/zopKGCnyOYg495r+01c66xt0vQPlQa3JgXAEA3xDhTMc5Mejm96j7X9zBra0zG+tCjW1AeFrjf72tBixAAQJe0tKIWhTGZowU/d2bSy/mQBFnUAdmPKdwdtBFfFVqEAAC6h6ZIlDMV5cwsj2CO5fIHs7kBRzmWJ7FuVB93KtqZtjHSdIm6Qx1BeP/+/ZSUFF9f33bt2tV7qrCw8Pz58zY2NtHR0TRN113fz88vODhYDeUBAOguI5r0c6f6uTOEkAdl/KHH/O/3uAmnWXdzqrMjFeFIdXagAmwodJ42QOVdoytXroyKitq/f/+AAQO++eabuk9dvXo1MDBwx44ds2bN6tu3L8uyhJDffvtNvn7//v3rrQ8AAA3wt6Y+DKIPxwqejRduimHa21Gn8vgBR1nbDTVd98vePs+uuMudzeefVWu6UC1D8aq8NEUsFru5ue3bt69Lly5paWlt27bNzMx0cHCQPxsbGxsVFTV37lyJRNKmTZuffvqpR48erq6uBw4ciIiIePjwYXBwcFZWlr29fd1tmpiYlJaWmpiYKFJARUWFpaWl8j8YvBzHcWKx2NzcXNOFGBaRSGRsbMwwjKYLMSy6coR5Vk1SnvG3S/iUZ3xKCX+nhGco0sKSamFFtbAkLa0od3PK1Zw4m1F2xpqutTE8z1dVVVlYWChxm6rtGj1//ryZmVmXLl0IIT4+PkFBQYcPH54wYQIhRCwWHzt2bMWKFYQQExOTQYMG7d2718jIyMLCIiIighDi6+vbunXrI0eOjBs3TqVFAgDot2bG8hOK//SPFklIRjmfUcFnVJBLBfzOKi5XRJ6IeJGMOJtRTqbEzoTYGVO1/7c1IlZGlKWQWAmJlRGxElKmAmKiL9+7VBuEubm57u7utT+6u7vn5OTIHz958oTjODc3N/mPbm5ud+/ezcnJedn6tXieX7t2rVAorF1iaWk5YsSIFxbAsqy8xxXUhuM47Hb1wz7XCN3d7c2EpJkd6WhXu0CekZSEJXkivkBMFVfzxRK+uJoUV/OZFaSkmlTU8OU1pFxKKmVUeQ0vlhEJS6yNiBFNLIUURRFrIU8IsRIShiY0IdZGf+eutRGh/xfBxgwx+3fs0IRYCUmTDHHmfE0V3e00TVONnSFVbRDW1NTU7asRCAQ1NTW1T1EUVfusUCiUSqUNrF+L5/lLly7VXc3Ozm7w4MEvK+D5LYBKcRyH3a5+NTU1NE1zHKfpQgyL/v1TZwhxNyHuCp16IoSQMimp4amKGp7jSXkNRQgpryEsRzhCymv+Waf2VhrVLBGz/4oljpBiSdPO0ImlTdjtRv/f3r2HNXGljwM/gBBAbgkQLolSqdxUShBEWAhdQK2sFlDsRSuKgKWgS+0F12dtfaxuu7WrbgGVFalbBDFF5FFu1gX0UQki5WIRFSRIuUMgBJBbEsj8/ji/nW+eIIgsl4a8n79mzpzMvDkMeTOZM+doaMxxIjQzM+vs7CRX+Xy+j48PXjY1NSUIoqury8TEBG8yNzcfW9/MzExun6qqqufOnZvkPUKJRDLJmmC6SKVSgiCg2WeZVCqFe4SzDz5hZv/dEwQxMECZ3maf2V6jLi4uv/32W2NjI0Lo+fPnJSUlbDabIIihoSF9fX0HB4f8/HxcMz8/n81mu7q61tfXNzU1IYT6+vp++eUXNps9oxECAABQcjN7RUin00NDQwMDAz/88EMOh7N+/Xo7O7vy8nInJ6fh4eEDBw7s27evp6ensrKytbV169atOjo6ISEhgYGBu3fv5nA4vr6+NjY2MxohAAAAJTezj08ghEZHRy9cuFBRUWFraxsaGkqhUDo7O1NTU/fu3aumppafn5+dnW1kZBQWFmZqaorrJyUlPXjwwNbWNiwsTENDfnQEeHzidw4en5gT8PjEnIBPmNk3E49PzHginHavlAiTk5ODgoJmOiQgq6mpqba21tvbe64DUS537txZtGjRkiVL5joQ5ZKamvree+/B94/Z1N7eXllZuW7dumnc53wedJsgCPzMIphN9+/fj4mJmesolM6ZM2fu3r0711EonbCwsOHh4bmOQrn8+uuvx44dm959zudECAAAALwUJEIAAABKTfHuEdra2lpaWpJTVUyssLDQw8NjpkMCsgQCQUdHx7Jly+Y6EOVSU1NDpVLJgXzB7CgqKnJ1dZ3kxxGYFj09PU1NTfb29pOsHxcX99J754qXCOvq6p48eTLXUQAAAFAAnp6eenp6E9dRvEQIAAAATCO4ogcAAKDUIBECAABQajM7xNosEIlE8fHxT548WbZsWUREhNxINOnp6QUFBUZGRnv37sWjew8PD8fHx1dXVy9fvvyjjz4aO3INmAyy2e3s7CIiIiiU/z+bp1QqvXPnTl5enlAotLGxCQ0N1dHRqaysvH79OvnanTt34lGEwKsqLS1NTk5WUVHZuXOno6MjWd7Q0MDhcMjVgIAAPDbh5cuXCwoKjI2NyfMfTEFmZub169epVOqePXsYDAZZnpKS0tLSQq6am5sHBQXFxMSQTxba2dn5+fnNdriKb3h4uLy8vKKiQkVFJTIyUm5rVVXVDz/8MDo6unXrVjc3N1yYm5ubnZ2tr68fERGxePHiVz2iwl8R7tixIzMzk81mX716ddeuXbKbTp8+HR0dvXr16o6ODnd3d5FIhOtnZWWx2eyMjIyQkJA5ilrhBQcHX7t2jc1mZ2ZmBgcHk+V8Pj88PFxdXX3lypVZWVmenp5isfiXX345f/688L8UdP62OVdRUeHt7W1hYcFgMP74xz8+evSI3FRXV3fy5EmyhcViMULo1KlT+/fvd3V1bW9v9/DwwOc/eFVJSUmRkZHOzs79/f1ubm79/f3kpufPn5Ntfvbs2cLCQoTQV199VV9fjwsHBgbmLnAFlpWVtXv37vT09L///e9ym3g8nru7u5GRkZWV1fr16+/du4cQ4nA4ISEhK1euFIlErq6uPT09r3xIQpHV1dVRKBSBQEAQRFdXF4VCqa+vx5tGR0ctLCxycnIIgpBKpSwWKyUlhcfjaWpqdnd3EwTR2dlJoVAaGhrmLnxF9ezZMwqF0tXVRRCEQCCgUCjPnj3Dm/A8pXgZjzhaXFycmJi4adOmOQt3vggKCoqOjsbLH3/88e7du8lNBQUFLBZLtvLIyIiFhUVubi5BEFKp1MHBITU1dTajnTfs7OzS0tLwsru7+9mzZ8fWwSOO3r9/nyAIKpXK4/FmNcR5Ki8vj8lkyhVGRUWFhobi5a+++iowMJAgCBaLlZycjAu9vb1jY2Nf9ViKfUV47949BwcHGo2GEDI0NFyxYgX+goAQamlpaWhowCNeqqioeHl5cblcXJ9KpSKEjIyMli9fTtYHk1dcXGxvb29oaIgQotFob7zxBtmMqqqq5DNVg4ODIpEI/3WePXt28ODBuLi41tbWuQpb0XG5XHIEV29vby6XK7tVIBAcOnToxIkTT58+RQg1Nzc3NjZ6eXkhmfN/9mNWdN3d3U+ePMHNiF7U7BiHw1m8eLGLiwtejY+PP3LkyI0bN2YvUKUx9r9gYGDgwYMHL/0bTUyxE2F7e7uRkRG5SqfT29ra8HJbW5uOjg45NrexsXFra+sE9cHkTbIZ9+7dGxgYaGVlZWRk5O7urqury+Vy7ezsHjx4MIvBzh9tbW1ksxsbG8u2uY6Ozpo1a7S1taurq1euXJmbm9ve3j72/J+DoBVcW1ubqqoq/jKHxm/G8+fPh4aG4uW33npLX19fLBaHhoZGRETMXqzKQfbDx9jYmM/nNzc3I4RkC6dwqit2ZxkKhSKRSMhVkUhE/udramrKbhKLxVpaWhPUB5M3mWb8y1/+8vTp04KCAoSQv7+/v78/Lt+zZ8/f/va39PT0WYt23pBtdrFYLNvmLi4u5OWIjY3Nl19++cMPP4w9/2cz2vlBU1NTKpWOjIzgXnUvbMaampqysrJr167h1UuXLuGFXbt2WVtb79+/H6YEmUZy/wXq6ura2toIIYlEgrvsTe1UV+wrQgaDgb8OYM3NzWSfLnNzc5FI1NnZiVebmprMzc0nqA8m76XNeOjQoZ9//vnGjRv6+vpyr3V2dm5qapqNKOcdJpNJNntTU9N4py5uYQaDMTw83NXVRdY3NzefpUDnETMzM1VVVbJr6AubMTEx0c/Pz9jYWK789ddfp1KpcLZPL9kPH3yem5qaqqmpyRZO4VRX7ETo4+PT3NxcXl6OECorK2tvb/fx8REIBHiyXzabnZqaihDq6+vLzc3dtGnTmjVrGhsbKyoqEEKlpaV8Ph+mzZsCb2/v1tbWsrIyhFB5eXlraytu9qysLITQ8ePH09LSbty4gW8iIoR6e3vxwsjISEZGhoODw1xFrtACAgLw+UwQBIfDCQgIQAhxuVwej0e2MEEQly9fdnBwMDY29vDwwPV7e3vx+T+HwSsobW3tdevW4Yu8wcHBzMxM3IzZ2dn4S8bIyEhKSgrZ/7y/v5/sFJ2fn9/f329nZzdHsc83169f5/P5AQEBHA5HKpUihC5duhQQEKCurr5x40b8NxoaGrp69epUTvX/pVfP70FsbKyJicm7775rYmJy+vRpgiBu3bqloaFBEERhYaGRkdHmzZttbGy2bduG68fExJD1z5w5M5ehK7JTp07R6XTcjKdOnSII4vbt2wsWLODxeAghS0tLp//Ky8vz8vJycXEJCAiwtLRksVhtbW1zHb5C4vP5NjY2Xl5enp6eK1aswJ2f2Wz2sWPHgoODHRwcAgICli9fbmlpWVVVRRDE3bt38flvbW39wQcfzHX4iqq0tNTY2Bi3bUBAAO4UrampWVBQQBBERkYGg8EYGRnBlbOysszMzDZu3Ojt7a2rq5uQkDCXoSus2tpaJycnKysrDQ0NJyen7du3EwRBpVJzcnL6+vocHR3d3NzWrl1raWnZ0tJCEMSvv/5Kp9P9/f3t7e3/9Kc/kX+OyZsPY43yeDz8QP3rr7+OEBocHGxsbLS1tUUIdXV13bt3z9zc3MnJabz6YGrq6uoeP34s2+wNDQ1Lly6V+y0IP8ddUVHR3d3NZDJZLBYM1T9lw8PDhYWFqqqq7u7u+I5IQ0ODjo6Orq5uRUVFR0eHqakpi8Uih4l44fkPXpVQKORyuXQ6fdWqVSoqKgih6urqxYsXa2trd3d3SyQS2cEKampqeDyelpbWG2+8IdunDEyeWCyWvflCoVAYDEZNTQ2DwdDR0ZFIJFwuVywWs9ls8nZgT08Pl8s1NDR0cXGZwifMfEiEAAAAwJTBd3MAAABKDRIhAAAApQaJEAAAgFKDRAgAAECpQSIEAACg1CARAmV36dIlPMbCzCEIIiMjYwpDICYnJ483LvzTp0+TkpJkx1FTIM+ePUtKShoaGpr8S/Ly8p48eTJzIQFlBokQzGdisdhyQhKJJDw8fKbHPr148WJ4eLiuru6rvjA6Ovry5csv3HTr1q3g4OCBgYFHjx4lJCQoVkYsLCwMDg7u6emprq5OSEggZ7KdQGlp6ebNm0dGRmYhPKBsFHvQbQAmpqamJjv98pkzZxBCsnNeq6mpffnll6tWrZq5GCQSycGDBz/77LMpJMIJuLi4HD16VEtL6+bNm1FRUdu2bVNXV5/G/c8oR0fHo0eP6urq5uTkhIeHb9q06aXD3//5z3/+7rvvfvzxx7CwsNkJEigPSIRgPlNTU/viiy/I1StXrqioqMiWIISio6PlXiUSifr6+mSHUR4YGCAIQkdHR65mf39/f3+/sbGxmpraeDFcu3atubl5+/btsoXPnz8nCEJPT++FL+Hz+dra2rKHE4vFQqHQyMiIPJCjo6Ojo+N4ByUJhUKEEJ6DE++nt7d37AjR2MDAQH9/v+w4KdOlp6eHQqGQ44DY29vb29tPUH9kZEQgEOjp6ZEv0dHR2bJlS1xcHCRCMO3gp1Gg7Kytrb/99luEEI1GO378eFRUlL6+Pp1Ot7Gxqays7Ojo2Lhxo66urr6+/oYNG3BeQQgVFxf/4Q9/0NPTMzMzo9Ppx48fH2//Fy5ccHV1ZTKZeDU0NJROp+vp6enr6zOZzBMnTuDyqKgoJyennJycJUuWmJiYfPrpp7hcKpVGR0cbGBiYmpoyGAw8jjZCKCUlxczM7MiRIwcOHEAILVq0iEaj0Wi03t7e2NhYGo1WUlLCYrFoNBoenvvixYsODg6ampp0Ol1HR+edd97Bw0YXFhbSaLTMzExfX189PT1LS0sOh0Oj0R4/fiz7Lj755JOlS5eKxWK5d3f9+nUajYbHmJUr+emnn2g0GpfL9fT0pFKpurq6np6e+EZpenq6mZnZ119//fHHH+M/AQ6+o6OjpaXl7bff1tTUNDU11dbWtrKyKikpwXsODAysrKyE+SzBtINECJRde3v78+fPEUJCofD48eNdXV35+fm5ubkikWjbtm3vvPOOs7NzcXFxYmLizZs3jx49ihB68OABHlW5oKDg4cOHn3766YEDB2JjY8fufHR09O7du25ubmSJSCSKj49/+PBhaWnpu++++/nnn1+8eBEhNDg4yOPxIiMjv/jii6KiInKi15SUlJKSkry8vPv3769evTooKOjWrVsIoYGBgfb2dj8/v127diGELly4kJaWlpaWtnDhwuHhYaFQ+N57723fvr2oqOjw4cP4bYaFhXG53MePH58+ffrOnTtBQUEIIYlEIhQKw8PDbW1tb968+dNPP/n7+6uoqJw7d46MeXBw8Mcff3z77bfJUUxJ+FKVnHJBtkQkEgmFwqCgID8/v5KSksTExIqKiv379+Mdtre3r1279qOPPkIIJSYm4uANDAzCw8Orq6tzcnLq6uru378fGRlJDgPp6uqqqqp68+bN/+0PDsAY0zlmOAC/bywWy9HRUa5QV1f3r3/9K0EQCCEWi4WnFyAI4l//+hdCCH8QY0FBQa+99hpBEBs2bLC2th4eHiY3ffjhh0wmc+wR8aXS2bNnxwtp7dq1vr6+BEHgzJeXlye71cTERE9PTygU4lWxWMxkMtetW0eGJxQKcQLGv7Vix44dQwjhWUHGk5SUhBASCAQ4rwQHB8tu3bdvn4GBAf5BmCCIxMREhBCe1ELO1atXEULV1dVjS/AhvvvuO3LTJ598snDhQvLora2tON3y+XzZt3zw4MHxwjY1Nd2xY8cE7wuAKYB7hAD8n3Xr1pFD11tbW+MScqu1tTWHw5FIJAUFBW+99dbdu3fJTTQarbm5WSgUknfjMPzzI41GI0skEklGRkZVVRWfz0cItbW1kb83Lly4cM2aNXIh+fj4GBgY4GV1dXU/P7/x+pHK8ff3lyspKyv7+eef29raJBIJnrO6rq4Ob/Lz85OtuWfPnpiYmIyMDHxrMyEhwdPTc/ny5ZM5rhxfX19yedmyZQMDAwKBYIL6jo6OeH60LVu2ODg4yM0kYGhoSM62DcB0gUQIwP8hUw5CCP8MKJvYNDQ0JBJJT0/P8PDwf/7znzt37si+lkql8vl8uUS4YMEChBDZ6b+zs9PDw6Ozs3P9+vV0Ol1TU1NLS4v8ZH9hLxVTU1PZVTMzs87OTpFI9NL3Ire3qKio06dPe3l5WVtbU6lUvIfe3l7c+0buKEuXLvX29k5ISNi+fXtlZWVJSQn+/XYK5BoQITT2RqOsf//739HR0bGxsd98842JiUlISMihQ4fIPqUSiUSBOscCRQGJEIBXs3DhQjU1tfDw8H/+858vrYyzEXkNlJKSUl9fz+PxFi9ejEs++OCDxsZGvPzCedTwNSWps7PTwMAAT0Y4Mdm9dXR0xMXF/eMf//j8889xSW5uLv59EsPT7MmKiIjYsmXLo0ePzp49a2houHnz5hceRS7TI4T6+vpeGtsETE1Nk5OTxWJxSUlJWloa7sf0zTff4K0CgWAmOrUCJQedZQB4NRQKxdXVNSsrazKPgTOZTHNz84cPH+LV+vp6Op1OZsHBwUHc82UCt2/fJkdgkUqlN27ckHvwAD+eOPEoLb/99htCSHZ63tzc3ImP6+/vz2Aw4uLiUlNTg4ODx3vOj8FgIJmfWHHAE+9Z1njBa2hoeHh4xMbGurm5FRcX48LW1laBQODi4jL5/QMwGZAIAXhlX3/9dUNDw6ZNm8rLy4eGhpqbm69evfrZZ5+9sLKPjw+Xy8XLLBarpaUlPj5eJBLxeLz3338fd1idQH9//65du1paWvh8flRUVE1Nzb59+2Qr4Ft3J0+evHfvXllZmWwHTpK1tbWWltbJkyfb29t7enpiYmLIxzDGs2DBgrCwsISEhN7e3t27d49Xzc7ObtGiRYcPH3706FFLS8uJEyeuXLky8Z5lLVu2DCH0/fffFxUVlZWVDQ0NhYSE3L59u7u7WywW5+XlVVVVkfm7sLBQRUXFx8dn8vsHYDIgEQLwyt58882srKy6ujonJydtbe1Fixbt2LFjvNG/QkNDq6qqqqqqEEI7d+7cunVrZGSkpqamjY0NjUbDDz9MIDg4eHR0lMlkmpiYnDt37tixY3K/Uq5aterw4cMpKSlsNtvZ2fmFmZVKpZ47d+727dtmZmZUKjUhIeHUqVMvfZthYWEqKipvvvmmjY3NeHUoFMr58+ebmppWrFjBZDKvXLmC+6xOkr29/bfffpuenu7p6ens7CwQCMrLy729vQ0NDSkUiq+v78aNG48cOYIrczgcb2/vJUuWTH7/AEyGCvHfZ3QAmPfRcRpAAAABRUlEQVSkUikacytOKpWqqKiMvUk2GTwer7u7m0ajWVhYTNCJw9nZ2cPD4/vvv8erra2tLS0tFhYWdDp9kgdqampqa2uzsrIi+54QBCGVSicY0Was/v7+2tpabW1ta2vrybzfsrIyZ2fnS5cuvf/++xPXFIvF1dXVmpqauKvtS00cvFAobG5uHhkZee2118j329bWtmTJkitXrmzYsGEyhwBg8iARAjDjioqK1qxZU1tbi++oKQSpVBoYGFhRUVFbW/t76KgZFRVVVVUFT9ODmQCJEIDZ0NDQYGhoOHa00t+niIiI7Ozs1tbWy5cvj9dfdJY1NTXp6+uPNzorAP8LSIQAAHnZ2dmdnZ2rV6/GnVkAmN8gEQIAAFBq0GsUAACAUoNECAAAQKlBIgQAAKDUIBECAABQapAIAQAAKDVIhAAAAJTa/wPmShtIHaEYtAAAAABJRU5ErkJggg==" }, "metadata": {}, "output_type": "display_data" } ], "source": [ "# Bootstrap a KDE of the pre-eruptive (or pre-deposition) zircon distribution\n", "# shape from individual sample datafiles using a KDE of stacked sample data\n", "BootstrappedDistribution = BootstrapCrystDistributionKDE(smpl)\n", "h = plot(range(0,1,length=length(BootstrappedDistribution)), BootstrappedDistribution,\n", " label=\"Bootstrapped distribution\", xlabel=\"Time (arbitrary units)\", ylabel=\"Probability Density\", \n", " fg_color_legend=:white)\n", "savefig(h, joinpath(smpl.Path,\"BootstrappedDistribution.pdf\"))\n", "display(h)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Run eruption/deposition age model" ] }, { "cell_type": "code", "execution_count": 12, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Estimating eruption/deposition age distributions...\n", "1: KJ08-157\n", "2: KJ04-75\n", "3: KJ09-66\n", "4: KJ04-72\n", "5: KJ04-70\n" ] } ], "source": [ "# Number of steps to run in distribution MCMC\n", "distSteps = 10^6\n", "distBurnin = distSteps÷100\n", "\n", "# Choose the form of the prior distribution to use\n", "# Some pre-defined possiblilities include UniformDisribution,\n", "# TriangularDistribution, and MeltsVolcanicZirconDistribution\n", "# or you can define your own as with BootstrappedDistribution above\n", "dist = BootstrappedDistribution\n", "\n", "# Run MCMC to estimate saturation and eruption/deposition age distributions\n", "smpl = tMinDistMetropolis(smpl,distSteps,distBurnin,dist)\n", "results = vcat([\"Sample\" \"Age\" \"2.5% CI\" \"97.5% CI\" \"sigma\"], hcat(collect(smpl.Name),smpl.Age,smpl.Age_025CI,smpl.Age_975CI,smpl.Age_sigma))\n", "writedlm(smpl.Path*\"results.csv\", results, ',');" ] }, { "cell_type": "code", "execution_count": 13, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "BootstrappedDistribution.pdf\n", "KJ04-70.csv\n", "KJ04-70_distribution.pdf\n", "KJ04-70_distribution.svg\n", "KJ04-70_rankorder.pdf\n", "KJ04-70_rankorder.svg\n", "KJ04-72.csv\n", "KJ04-72_distribution.pdf\n", "KJ04-72_distribution.svg\n", "KJ04-72_rankorder.pdf\n", "KJ04-72_rankorder.svg\n", "KJ04-75.csv\n", "KJ04-75_distribution.pdf\n", "KJ04-75_distribution.svg\n", "KJ04-75_rankorder.pdf\n", "KJ04-75_rankorder.svg\n", "KJ08-157.csv\n", "KJ08-157_distribution.pdf\n", "KJ08-157_distribution.svg\n", "KJ08-157_rankorder.pdf\n", "KJ08-157_rankorder.svg\n", "KJ09-66.csv\n", "KJ09-66_distribution.pdf\n", "KJ09-66_distribution.svg\n", "KJ09-66_rankorder.pdf\n", "KJ09-66_rankorder.svg\n", "results.csv\n" ] }, { "data": { "text/plain": [ "6×5 Array{Any,2}:\n", " \"Sample\" \"Age\" \"2.5% CI\" \"97.5% CI\" \"sigma\"\n", " \"KJ08-157\" 66.0703 66.0343 66.0935 0.0149756\n", " \"KJ04-75\" 65.9318 65.886 65.9707 0.0218798\n", " \"KJ09-66\" 65.9368 65.8915 65.9779 0.0219181\n", " \"KJ04-72\" 65.9551 65.9219 65.9766 0.0138594\n", " \"KJ04-70\" 65.8307 65.7535 65.8949 0.0365122" ] }, "execution_count": 13, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# Let's see what that did\n", "run(`ls $(smpl.Path)`); sleep(0.5)\n", "results = readdlm(smpl.Path*\"results.csv\",',')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "To see what one of these plots looks like, try pasting this into a markdown cell like the one below\n", "```\n", "