{
"cells": [
{
"cell_type": "markdown",
"id": "d99e28b8-e335-4e56-952a-ddbef596b1d4",
"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 *Pb-loss-aware* 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, "id": "2cdabf6e-04fb-4082-a8e9-a24b6f60d2f2", "metadata": {}, "outputs": [], "source": [ "# Load (and install if necessary) the Chron.jl package\n", "using Chron\n", "using Plots, DelimitedFiles" ] }, { "cell_type": "markdown", "id": "cdee6490-0e43-4c38-b334-bb1e141f80dc", "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, "id": "6382e055-7c34-4878-8e90-71c5204abe6b", "metadata": {}, "outputs": [ { "data": { "text/plain": [ "\"m\"" ] }, "execution_count": 2, "metadata": {}, "output_type": "execute_result" } ], "source": [ "nsamples = 3 # The number of samples you have data for\n", "smpl = ChronAgeData(nsamples)\n", "smpl.Name = (\"KR18-04\", \"KR18-01\", \"KR18-05\")\n", "smpl.Height .= [ 0.0, 100.0, 200.0] # Arbitrary example heights\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?\n", "smpl.inputSigmaLevel = 1 # i.e., are the data files 1-sigma or 2-sigma. Integer.\n", "smpl.Age_Unit = \"Ma\" # Unit of measurement for ages and errors in the data files\n", "smpl.Height_Unit = \"m\" # Unit of measurement for Height and Height_sigma" ] }, { "cell_type": "markdown", "id": "89c5c1a2-dfdd-4745-8d87-686a538dd99f", "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, "id": "8d66999a-c19f-4bb4-8cbe-f483309d0338", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Chron1.0Coupled.ipynb\n", "Chron1.0Coupled.jl\n", "Chron1.0CoupledConcordia.ipynb\n", "Chron1.0CoupledConcordia.jl\n", "Chron1.0CoupledSystematic.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", "ConcordiaExampleData\n", "DenverUPbExampleData\n", "EruptionDepositionAgeDemonstration.ipynb\n", "EruptionDepositionAgeDemonstration.jl\n", "FCT\n", "Manifest.toml\n", "MyData\n", "PlutonEmplacement.ipynb\n", "Project.toml\n", "archive.tar.gz\n", "ffsend\n" ] } ], "source": [ ";ls" ] }, { "cell_type": "markdown", "id": "d34a15e4-4874-4761-84fa-0a9c1b554552", "metadata": {}, "source": [ "Equivalently, we can also run unix commands using the `run()` function, i.e.,\n", "```\n", "run(`ls`);\n", "```" ] }, { "cell_type": "markdown", "id": "a4cc88e6-4581-4949-a495-9326bdc966ca", "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": 4, "id": "914297ee-a0cf-49b8-851d-cc48046f62c4", "metadata": {}, "outputs": [], "source": [ ";mkdir -p MyData/" ] }, { "cell_type": "markdown", "id": "4be2117f-bf68-4bb0-b678-1d04397738e8", "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 MacLennan et al. (2020), [10.1126/sciadv.aay6647](https://doi.org/10.1126/sciadv.aay6647)" ] }, { "cell_type": "code", "execution_count": 5, "id": "b2d06f69-6298-440b-a747-f8e616695498", "metadata": {}, "outputs": [], "source": [ "# You can just paste your data in here, in the following five-column format:\n", "# ²⁰⁷Pb/²³⁵U, ²⁰⁷Pb/²³⁵U sigma, ²⁰⁶Pb/²³⁸U, ²⁰⁶Pb/²³⁸U sigma, correlation\n", "# You should generally exclude all systematic uncertainties here, and all uncertainties\n", "# (sigma) must be absolute uncertainties\n", "data = [\n", "1.1002 0.00060511 0.123908 0.00001982528 0.333\n", "1.1003 0.0005226425 0.123893 0.000020442345 0.421\n", "1.1 0.0011 0.123874 0.00002353606 0.281\n", "1.1002 0.00060511 0.123845 0.000025388225 0.418\n", "1.1007 0.0005338395 0.123833 0.000025385765 0.534\n", "1.0991 0.001154055 0.123797 0.000031568235 0.298\n", "1.09931 0.0004067447 0.123762 0.00003898503 0.709\n", "1.09947 0.0004617774 0.123752 0.00002598792 0.579\n", "1.0986 0.00093381 0.123738 0.00003650271 0.288\n", "1.09883 0.00047799105 0.123735 0.0000222723 0.506\n", "1.09904 0.000384664 0.123733 0.000021653275 0.404\n", "1.0758 0.0005379 0.121175 0.00002302325 0.427\n", "]\n", "\n", "# Now, let's write this data to a file, delimited by commas (',')\n", "# In this example the filename is KR18-01.csv, in the folder called MyData\n", "writedlm(\"MyData/KR18-01.csv\", data, ',')" ] }, { "cell_type": "code", "execution_count": 6, "id": "054f2fe6-9f67-44fd-85c4-058508428fe1", "metadata": {}, "outputs": [], "source": [ "data = [\n", "1.1009 0.000935765 0.123906 0.00002849838 0.319\n", "1.1003 0.00077021 0.123901 0.000035311785 0.415\n", "1.0995 0.000494775 0.123829 0.000025384945 0.434\n", "1.0992 0.00060456 0.123813 0.000036524835 0.616\n", "1.1006 0.00071539 0.123813 0.00002228634 0.321\n", "1.0998 0.00076986 0.123802 0.00002537941 0.418\n", "1.0992 0.00065952 0.123764 0.00003589156 0.509\n", "1.0981 0.0010981 0.123727 0.00003959264 0.232\n", "1.0973 0.000526704 0.123612 0.00002966688 0.47\n", "1.0985 0.0008788 0.123588 0.00002842524 0.341\n", "1.0936 0.0005468 0.123193 0.000032646145 0.575\n", "1.0814 0.000513665 0.121838 0.0000304595 0.587\n", "]\n", "writedlm(\"MyData/KR18-04.csv\",data,',')" ] }, { "cell_type": "code", "execution_count": 7, "id": "5181cfcd-d8cd-4ee4-843e-b8911f563ade", "metadata": {}, "outputs": [], "source": [ "data = [\n", "1.09741 0.00038958055 0.123676 0.00002226168 0.601\n", "1.097 0.00060335 0.123629 0.000024107655 0.575\n", "1.0967 0.00054835 0.123618 0.0000247236 0.41\n", "1.0952 0.00120472 0.123594 0.00003151647 0.366\n", "1.0969 0.000712985 0.123553 0.000035212605 0.581\n", "1.0956 0.0005478 0.123547 0.000032739955 0.562\n", "1.0958 0.00071227 0.123546 0.00002779785 0.489\n", "1.09621 0.00046588925 0.123539 0.000033973225 0.75\n", "1.09612 0.0004822928 0.123533 0.000032736245 0.747\n", "1.0969 0.00076783 0.123505 0.0000419917 0.552\n", "1.0958 0.00060269 0.123469 0.00002345911 0.342\n", "1.09517 0.00037783365 0.123447 0.000030244515 0.783\n", "1.0948 0.000525504 0.123352 0.00003885588 0.737\n", "1.09413 0.0003720042 0.123288 0.00002527404 0.508\n", "1.09321 0.00046461425 0.123179 0.00003202654 0.69\n", "1.08592 0.0004017904 0.122429 0.00002938296 0.634\n", "]\n", "writedlm(\"MyData/KR18-05.csv\",data,',')" ] }, { "cell_type": "markdown", "id": "36cc91c7-77e1-413c-8e87-13f03bd13dc4", "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": "markdown", "id": "78ade5fb-6487-4842-90b0-e77a3ff1300f", "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 Samperton (2018)](https://doi.org/10.7185/geochemlet.1826) and for the Pb-loss-aware version we are using in this notebook, [Isoplot.jl](https://github.com/JuliaGeochronology/Isoplot.jl)\n", "\n", "#### (optional) Boostrap relative pre-eruptive (or pre-depositional) mineral age distribution" ] }, { "cell_type": "code", "execution_count": 8, "id": "8918d67a-fb38-46a6-97a4-3fc584a4eab6", "metadata": {}, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ "\u001b[36m\u001b[1m┌ \u001b[22m\u001b[39m\u001b[36m\u001b[1mInfo: \u001b[22m\u001b[39mInterpreting the five columns of KR18-04.csv as:\n", "\u001b[36m\u001b[1m└ \u001b[22m\u001b[39m | ²⁰⁷Pb/²³⁵U | 1-sigma 𝑎𝑏𝑠𝑜𝑙𝑢𝑡𝑒 | ²⁰⁶Pb/²³⁸U | 1-sigma 𝑎𝑏𝑠𝑜𝑙𝑢𝑡𝑒 | correlation coefficient |\n", "\u001b[36m\u001b[1m┌ \u001b[22m\u001b[39m\u001b[36m\u001b[1mInfo: \u001b[22m\u001b[39mInterpreting the five columns of KR18-01.csv as:\n", "\u001b[36m\u001b[1m└ \u001b[22m\u001b[39m | ²⁰⁷Pb/²³⁵U | 1-sigma 𝑎𝑏𝑠𝑜𝑙𝑢𝑡𝑒 | ²⁰⁶Pb/²³⁸U | 1-sigma 𝑎𝑏𝑠𝑜𝑙𝑢𝑡𝑒 | correlation coefficient |\n", "\u001b[36m\u001b[1m┌ \u001b[22m\u001b[39m\u001b[36m\u001b[1mInfo: \u001b[22m\u001b[39mInterpreting the five columns of KR18-05.csv as:\n", "\u001b[36m\u001b[1m└ \u001b[22m\u001b[39m | ²⁰⁷Pb/²³⁵U | 1-sigma 𝑎𝑏𝑠𝑜𝑙𝑢𝑡𝑒 | ²⁰⁶Pb/²³⁸U | 1-sigma 𝑎𝑏𝑠𝑜𝑙𝑢𝑡𝑒 | correlation coefficient |\n" ] }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAlgAAAGQCAIAAAD9V4nPAAAABmJLR0QA/wD/AP+gvaeTAAAgAElEQVR4nOzdZ0AU194G8DMzS5PeeweRagF7xYhi7AXFXmMvCSa2a4yaRGNiiV2T6BujMfaOihA1GruIgI0mSpEiVTq7M/N+2HtRcZEiu8vuPr9Pu8Nx+LPuzDNnzpkZiud5AgAAoKpoeRcAAAAgTwhCAABQaQhCAABQaQhCAABQaQhCAABQaQhCAABQaQhCAABQaQhCAABQaQhCAABQaQhCAABQaQJ5F1C75cuXX7161czMrI7teZ6nKEqqJcH78LHLCz55ecEnLy91/+SLiopMTEz27dv34WYKEIS5ubk2NjaDBg2qS2OWZTmOU1NTk3ZVUE1FRYW6ujr2C7JXVlampaUl7ypUET55uWBZlmVZdXX1ujS+e/duVFRUrc0UIAh1dXXt7OyCgoLq0lgoFLIsq6mpKe2qoJqSkhItLS2axsl2WSsqKtLV1ZV3FaoIn7xciEQioVBYx0MQNTW1hISEWpthtwUAACoNQQgAACoNQQgAACoNQQgAACoNQQgAACoNQQgAACoNQQgAACpNAa4jBFAFJSISmcPHF/IJhXx8IYkv5LPKeDMtylyLWGtTZprEVof6xIryMcItCwAaGYIQQM6icvlfnnKHn3Fu+lQLA8pVnxrrQlz1aatm1KtyPquMpJfw2eXk2Wt+yyNOgyFBjlSQI+2NRJQ3hmHkXQI0DgQhgHwUC8lfSdyvcdyrcjKlOR07TGDVrHq2mWhS7gaEkP8u30LI3Vf8kWRuYDiryZA5HvT0FhjdkBvcX01pIAgB5OBgEhdym+1gRq/yZXpbU3Sde3dtTam2pszaduR2Nv9NJLvtMfetDzMM9/mSB9xZV2kgCAFkKrOMzL7OxhfyJ3oJ2ps1cE9KEdLBjArrK4hI5z+/KdiRKFrfgWmJk6UADYLzKgAywhPyy1PO55iwpTEVOaThKfi2XtbUtYCKoQ50n/OieTfZCvbjVwmgchCEALLwqpx8Eir6v3jucj/B8ta0euNteQKazPKgnwapZZSSbmdFqSV8o60aQDUgCAGk7lkR3/mMqIsF9e8AgaehVE5gGqiTw58wwx3p9qdElzOQhQD1gCAEkK7YPN4/lA3xolf5Mow0R/EoQr7yoQ/2FIy7wq6N5hCGAHWEIASQovB0vtd50ZaO9Ax3GW1r3SyomwOZY8+54EsYMgSoEzkHYVlZWXl5uXxrAJCS/Ync+CuiE70EA+1luqHZalPX+gsIIUMiRKqchUVFRfn5+fn5+SUlJfKuRYo2btxYVFT0gQanTp168OABISQmJiY0NLQBvyI6OvrcuXOEkJycnK1btzaszvv374eFhTXs30qbtLbPmzdv9u/f38bGpkOHDhIbhIWFeXl5mZiYmJqadu/e/fnz51KqBEAufnnKfR3JXeon6GQuh6saNBjyZw9GT41S5SwcNGhQ8+bN/fz8nJycHBwcGpABycnJo0aN+kCDS5cuLV269CNqbATLly8vKCj4QIO//vrr9u3bhJCYmJgzZ87U1GzgwIHZ2dkSf/TgwQPxp5ednf3999/XsbCcnJz+/ftXvY2MjDx//nwd/62MSSsINTQ0goODQ0JCXr9+LbGBurr6rl27ioqKcnNzbW1tZ82aJaVKAGQvNJX/JpK9GMi4G8jt2j4BTf70Z4w0qEHhonJVzcLly5cnJSVlZWUtXLhw6tSpb/+ooqKC47hq7VmWraysrHpbXFx87dq1am0KCwtZ9r8faEZGhrizVc3byVRTd624uPj9hRUVFUKhsNpCnuffP3MmEokqKiokrrmmXzF27NidO3dWvS0vLy8rK6t6e/ny5fd/S2lpaUVFxYQJE7Zt2/b28vLy8qoPoSYVFRWXLl2qevvZZ5/9/PPPbzd4+7dXEYlEEpdLlbSCsE2bNmPHjnV0dKypgb+/f+fOnWmaVldXHz58eFxcnJQqAZCxezn8pH9ExwMErvpyvsKdocje7oyJJjVYhbNQzMfHp6SkhOd5QkhaWtonn3zi4uJiZWX1xRdfiEQiQkhFRcWMGTOsra2dnJw+/fTTrKwsQsjw4cOzsrKcnZ2dnZ1zcnKOHj1qY2PTsWNHW1vbb775JiUlJSQk5J9//nF2dm7fvj0hZOjQoYsXL3Z1dW3evHlBQcGCBQtsbGx8fX1tbGyOHj0qrmTBggUzZ85s06aNt7e3t7f3kydPCCH5+flGRkarVq3y8vKytrZeuHChuNSKior58+fb29t7eHj07NkzPT1dvJLly5dbWVl5eXktXrxY4t8bFhZma2vbsmXL3r17FxYWihfu3btXfDTw6tWrrl27urm5+fj4uLi4iESi8ePHl5aWdu3a1dnZOTw8/PTp03379h02bJiLi8vWrVv37Nkzffp08Up4np83b56Xl5elpeXatWvFC0eOHHnkyBHx62PHjgUFBRFChg0bVl5eLv70kpKSdu7cOWfOHHGbXbt2WVlZubu7e3l53blzR7zQ0dHxhx9+cHNzs7GxGTZsWK1B24iaxJ1ljhw50rNnz5p+KhKJ8vLynj17Jn4rEAjs7OxkVRpA/SS95gdeFO3pJujYGNfLfzyGIr93Y8ZeYYdFiI73EmjI9jbRSa/55A+NXjUmy2bk/UtTTp8+nZ6eXl5efuXKlfXr14tvijZ9+nRPT8+IiIjCwsJu3brt3r17+vTpmzdvjo6OTkpK0tLSmjx58vz58w8ePHj06NG+ffsmJSWJ17Zw4cJjx46JMy8jI8PS0nLDhg1//vmnePyMEFJcXBwWFnbnzh0DAwNCyGeffbZu3TqKoh4+fOjv79+nTx9dXd3S0tLjx49HRkba2NisW7duwoQJd+7c4Xk+Pz//9evXCQkJubm5HTp06NSp0+DBg9euXZuSkhIfH6+pqblmzZq5c+ceP3784sWLe/fuffjwoZmZ2fLly9/vWb5+/XrMmDGHDh365JNPbty40a1bt8GDBxNCysvLxd3T3bt3N2/eXNzZzczMpGn6jz/+OHHixLVr18R716NHj4aFhYWGhh47dozjuJ07d1b1a7Oysuzs7BITE9PS0nx9fbt169axY8fi4uKqnnRlZaW4pGPHjrm6ulZ9eqdPnxYvv3fv3uLFi+/cuePq6rp3797hw4eL/8CCgoL4+Pj4+PiKigo/P7+zZ88OGjSoEb8hHyD/INy1a9f169fv3btXU4PHjx/funXr8OHD4rcMw1y8eNHMzExiY6FQyLKs+BAPZKm0tJRlWZpW6XnIuRVUnwi1xZ5sD6MKSee9pKKkpKTWm17u9CMTb6pNuMz+2kEovXzW1NQUCN7ZpRx/zl9Mr376UUpaG1M/tque8yYmJk5OTkVFRcbGxlevXp06dapQKAwLC3v+/DlFUQYGBtOnTz99+vT06dPPnDkze/ZsbW1tQshXX33Vrl07cZ/sbZaWlps3b546dWrnzp0tLS0lljFlyhRDQ0Pxaxsbm7179z579qyyspLn+bi4OD8/P0JIUFCQjY0NIWTu3LlLlizJyspSU1MjhHzxxReEEGNj4/Hjx585c2bw4MF//fXXpEmT/v33X0KIq6vrmjVreJ4/e/bsmDFjxPvAkJCQb7/9tloNN27csLGx+eSTTwghnTp1en+ihpWV1Z49e37//fe+fftaWFhI/EM8PT379u1LCKm2UQsEgnnz5on/uqCgoDNnznTs2FHiGmpy9uzZIUOGuLq6EkLGjx+/ZMmSBw8eiIucN28ewzDNmjXr2rXr+6cJeZ4vKSkRiUTi/Xxdfld5efn7/4/vk3MQ7t+//7vvvrt06ZKJiUlNbXx8fLp3775o0aK6rFD8AWlqajZejVAnFEVpaWmpchCWicjoy6JRLtS8luqy/L08z+vo6NTa7GAv0v2saGuS2pKWsvs/+sqH/spHnl+JTp06TZs2jRAyd+5cCwuLS5cueXt7syxbtcMxMTHJzc0lhOTl5VUtNDU1LS0tfX+k6vjx41u2bAkJCUlOTt60adOECRPe/42mpqbiF5WVlR07duzWrVu3bt1omv7zzz+rOlXGxsbiFxoaGjo6Orm5ueI0qlpuYmJy//59Qkh2dvadO3eqOlWjRo0SnyGrOiumr68vDtG35eXlVa3q7dVWGTduHMMwhw8fnjNnTq9evQ4dOqShoVHTH1KNjo6Ourp6VZ0ZGRnVGtQaPG+XR1GUiYlJTk6O+K2enp74hYaGxvsjoBRF6ejoiIOwjo/+0NTUrMu90eX5HT169OiiRYsuXLggPjQAUGgTr7LOutQq3yb6jDpNhpwIYLY/5k69kFEXrUkR7w2Li4tNTU319PRiY2PFy2NjY11cXAghTk5OVQtjYmKsrKyaNWumpqb2ds/D3Nz8u+++i4qK2r9/v/jQXE1NrabzT/Hx8fn5+du2bRs5cmRgYGB+fn7Vjx49eiR+kZaWVlJSYm9vX235w4cPnZ2dCSFubm79+vXb9RY1NTUnJ6fHjx+LWz59+vT9yTVOTk5xcXHiwnier1rt25/GmDFjTp06lZGR8fTp0wsXLoj/lrf/2Jryo6CgoGqo8uHDh+JPT19fv+oPTExMFL+otsIqzs7OVR91UVHRs2fPxCuRI2n1CAsLC+/evRsTE1NSUhIREWFoaOjr60sIGTJkyLx58/z9/c+fPz927NgffvghIyMjIyODpukPDBMCNHE7n3BPCvjbAwVNYmCwBlbNqFMBTN8wkaOuqjzp/saNGxoaGuXl5efOnTM2Nvb396coKiQkZMaMGT/++GNaWtqOHTvEc/o///zz4OBgW1tbQ0PDzz//fMGCBYQQOzu7srKyTZs2WVlZDRgw4Kuvvurfv7+xsfHly5c9PDwIIZ6enlFRUXv27DE1NR0wYMDbv9ra2rq4uPj//u//XFxcNm7c+PZJ43/++WfTpk2+vr6rVq2aPHmytra2uPezdOnSJUuWPHv2rOqChxUrVkyYMIHneU9Pz9TU1CdPnnz99ddTp05t3br1rl273N3d16xZ836PsH379ra2tnPnzh03btzhw4fz8vKqNdi9e7eampqbm1t2dnZ+fr64K+Lp6blhw4Zu3bp16tTpAx+ppqbmzJkzFy5cGBkZeeXKle3btxNC/P39N2zY0KJFi9TU1H379omDzcTExNDQcM2aNS4uLoGBgVVrGD9+/OrVq9esWdO5c+dNmzb16NGjRYsW9fyPbWTMihUrpLHeFy9erFu37tWrVzY2NpGRkfn5+T169CCEhIeH+/n52dnZRUZGchyXmpoaGRkZGRkZExMzdOhQiauKiIjQ0tLq0qVLXX4vx3E8z1cbqAAZEAqFampqqvmEtkf5/MR/2HN9BJbvPVlXBiorK98/r1UTy2aUoy41+Ro32pnWrr7/VDaVlZUcx+Xn55eWlnbu3Hnr1q26urqEkO7du6urqx86dCg7O3v9+vXiyS+Ojo7t2rU7ceLE3bt3P/vsM/EJVTU1tV69esXExKSmpnbp0uXFixfnz5+/ePGira3tunXrtLW1zczMvL29Y2Njc3JyunXrVllZ2bJlS3Nzc0KIlpZW586dDx48eOfOnc8++8zLy6tt27ZGRkahoaH+/v5CofDYsWNdu3b99ttvGYYpKyv78ccfDx06tH379pSUlE2bNvn4+BBCnJ2de/bsee7cudDQ0JycHPFkV319/YCAgOPHj9++fXvhwoUWFhb+/v5vjwdRFDVs2LBbt26dOXOma9eugwYN8vDwsLGxEQqFxsbGLVu2LCsru3DhwtmzZ1NSUlavXi0euQwICEhOTk5OTnZzczM2NjYwMBD3XsSfpKmpqY+Pj/g8/MiRI7dt25afn//LL7+ILw1o1aqVSCQ6cuSIQCBYtGiRmZlZ69atKYoKDAx8+PBhSkpKu3btBAKBubm5t7e3lpbWyJEjIyIiLl686OfnV3WUUF5e3r17d/E5T6FQ6OLiIvG6A47jOI57P/4levr0aXR0dHBw8IebUXUZSJSvxYsXGxoaYoywiSspKVHNMcJSEWl7UrS4FT3ORT5/e1FRkXj/XnfL7rFXMvi/P5X1JFIghMycOdPW1rbaZfjiYbOmvzduCuo1Rnjy5Mm9e/eeOHHiw81UbrcF0Ljm3mDbmlLySsGGWeXLmGlRX9xS7UsLAf4HpxABGu7wM+5aFh85WMG2I5oif3RnfE+KTr3gBsn2PqiwY8eO9xcaGRmhOyhH2AYAGijpNT/nBnuoJ6OrgINtOmpkfw9m5nU2o1TepQDIG4IQoCFEHAm+xH7ThmltrKjzg9qaUjPcmYn/iNATARWHIARoiPWxnLEmme2h2FvQ0pZ0iYhseqiKVxYCVFHszRhALhIK+Z9i2O2dFX7OpYAm+3owqx+w0XnoFoLqQhAC1A9PyKwb7LLWjJOuop4UfZujLrW+AzP6EluGG/SCqkIQAtTP7jjudSWZ66k82844F7qlMbXoLq6mqB9ZPicIpEp5NmYAGcgoJf+5x/7WlWGUoTf4xvbOzInn/L+ZOEFaD6WlmHGrJBCEAPUw9yY7rQXtrXQ36jRQJxs60DOus5WYNwOqB0EIUFdnU/jYPP4/rRR+joxEQY60ky61LgZJCCoHQQhQJ4WVZOZ19reujKZy5iAhhGzpRG98yCa9xglSUC0IQoA6WRXFBtpQXS2U7aTo2+x1qEUtmVnXMQcEVAuCEKB2ia/5PxK4b/2UtzP4P5970q/KyV9JOEEKKgRBCFC7z2+yS1oyFnV68ItiE9BkVxfmy9tcfoW8SwGQFQQhQC0i0vmnhQp/N7W6a2tKDXagluCyQlAZqrJtAzSMkCPzbrIbO9Aq9Qzb1X7M2VT+djZmzYBKQBACfMj2x5xVMzLATrW2FH118r0f/fktFkkIqkC1Nm+AesmrIGui2Z87qlJn8H/GudBCjhxLxqwZUH4IQoAafR3JjnSivQyV+ZKJmtAUWdeeWXiHq8BYISg7BCGAZI8L+KPJ3PI2qtgdFOthSXkZUlsfo1MISg5BCCDZl7fZZa0YYw151yFXP7aj10azubiUApQaghBAgmuZ/NMCMt1d1TeQFgbUSCd61X2cHgVlpurbOYBEi++y3/nR6tg+CFnpyxx6xsUXYgIpKC1s6ADVnXzBFQtJsBO2DkIIMdIgC7yZRXcwUghKC5s6wDtYniy7x61py9CqOFdUsnmedHQe/08GOoWgnBCEAO/Yl8AZa5BPbRGDb2gw5Ds/+j/3MFIIyglBCPBGJUe+jeLWtFXdSyZqEuxEvxaSsDR0CkEJIQgB3tj6iPMxojqZoztYHU2RZa3oZfdw0zVQQghCgP8qFpJ1sex3ftgoJAtyooUcOZeKKARlg20e4L9+jGH72NCeKnlDtbqgCFnehl4eiU4hKBsEIQAhhORWkO2PuRVtsEV8yBAHmufJ6Re4lAKUCjZ7AEIIWRfDjnCi7XXQHfwQipAVbeivIzkOvUJQIghCAJJbQX59yi1qic2hdgPtaS2GnHiOTiEoD2z5AOTHaDbYGd3BuvqmDfPNfXQKQXkgCEHV5ZSTPfHcYnQH6+xTW0pXjRzFM3tBWWDjB1W3Npod5UzbaKM7WA8rfZlVUegTgpJAEIJKe1VOfk/gFvlgQ6if3taUrho5g+mjoBSw/YNKW/OAHeNCW6M7WH8LvOnV0QhCUAYIQlBdmWVkbwL3lTe2goYY6kDnV5BrmTg/CgoPuwBQXT9EsxNc0R1sIJoiX3rTa6PxSApQeAIprZdl2adPn8bFxXl5eTVv3lxim5cvX544cYKm6aFDh5qbm0upEgCJMkrJ/gTu4XA1eReiwMa70ivuc7F5vLcRDiZAgUmrR9itW7eAgICpU6eGhoZKbJCUlOTj4xMTE3P37t2WLVumpaVJqRIAiX6KYce70hZa8q5DkWkwZL4X/WMMRgpBsUkrCENDQ1++fNmhQ4eaGmzcuHHo0KG7du3as2dPr169tm7dKqVKAN6XU07+SOAWYHTwo81wp8PSuOdFGCkEBSatHYGBgcGHG1y8eHHgwIHi1wMHDgwLC5NSJQDvWxfLBjtjdLAR6KmRKW70hofoFIICk9YYYa1evnxpaWkpfm1pafny5cuaWmZmZt68eTMnJ0f8VkND44svvtDR0ZHYWCgUsixLUdjByVpFRQVN0zStAH2svAry61P6Vj+uokIk71oaQUVFhbq6uhwLmNWcanWKLPJkTTRUq18o909eNYlEIqFQWMddjVAo5Pnav5Zy221RFFVVH8/zH44uTU1Ng//R1tZGzsHH2PqUGmzH22nLuw5lYabJD7UnO57Kuw6AhpJbj9DCwiIzM1P8OjMzs6p3KLGlu7v7okWL6rJamqZZltXQ0GicKqHORCKRhoZG0+8RFlaSX+KFtwcJNDSU5HCqsrJS7l/4Ra35zmdES9qoa8ttjyIHTeGTV0EMw9A0XcdPXk1NrS4dJ5nutkpKShITE8Wve/fuffbsWfHrs2fP9u7dW5aVgMra+pjrZ0s76SpJCjYRLnpUVwv6jwSMFIJCktbx2/bt269cuRIVFfXy5cubN2/OmzevS5cuV69eHT16dH5+PiHkiy++6NChA0VRQqHw4sWL9+7dk1IlAFVKRGTLI/af/qrUbZGV+Z701Gvs9BY0jWMMUDTS2iO0bdvW1NQ0KChI/NbOzo4Q0qZNm/3794uXuLi4xMTEiC+oX716tZmZmZQqAaiy7THX04p208euuvF1taD01Eh4Ot/HBh8vKBgpBmHbtm2rLTQ3N+/Xr1/VWysrq9mzZ0upAIBqylmy6SF3IZCRdyFKa44nvekR28cGHW5QME19agNAY9n1hOtoTuFmYNIT7ERH5fBxhap1EQUoAQQhqAQhRzY85JbgMfTSpMGQ6e705keYMgMKBvsFUAl74jlvQ+Jrgu6gdM1yZw4mcXkV8q4DoD4QhKD8WJ6si+GWtMLooNSZaZF+tvT/xaNTCIoEQQjK789Ezk6HdDZHd1AWPveitzziRIhCUBwIQlByHE9+jOH+g+6grLQxoWx1yOkUJCEoDAQhKLljzzkdNdLTCt1B2ZnvSW/C8yhAcSAIQZnxhKx+wH3dGt1BmRpsTz8vJlG5uI4CFAOCEJRZaApPEfKpLbqDMiWgyUx3ettjdApBMSAIQZmtiWaXtsLNL+Vgiht9/DmXj+soQBEgCEFpRaTz+RVkqAO+5HJgqkl6W9P7E9EpBAWAfQQore8fsEta4WEIcjPTnd75BEEICgBBCMrpZjafUkxGOeEbLjfdLSmKIlczMWUGmjrsJkA5rYhkl7aiBfiCy9VnbugUggKQvJ8ICQm5ceOGjEsBaCz3c/inhWScK2JQziY0py+kcdll8q4D4IMk7ykuXLjQuXNnHx+frVu3FhYWyrgmgI+04j63uCWtjhyUNwN1MsQetx6Fpk7yriI2NjY8PLxFixaff/65hYXFiBEjIiIiZFwZQMNE5/H3criJ6A42DTPc6V1POQ4DhdCESd5ZMAzTq1evw4cPv3jxYsWKFbdv3w4ICPDw8Fi7dm1eXp6MSwSol2+juIU+jBYek940tDWljDXIxXQkITRdtRw1W1tbL1q0KDExMSQk5MmTJ4sXL7a3t58/f35WVpZs6gOolycF/LVMbqobuoNNyHRcRwFNWy37i/z8/E2bNrVq1WrDhg2enp6bN2+ePXv2/v37/fz8CgoKZFMiQN19F8Ut8GZ01ORdB7xltDN9PYt7UYxOITRRNQbhv//+O378eCsrq4ULF/r4+Pzzzz8PHz6cO3fuDz/88PDhw/z8/PDwcFkWClCrhEI+4iU30x3dwaalmYCMcqZ3x6FTCE2U5IGUTp063bx5097e/uuvv54yZYq5ufnbP7W0tLSzs8NsUmhq1kRzczwYXXQHm57pLejAC+w3bQiDG/1A0yM5CF1dXZcuXdq3b1+Gkfz8mlu3bmloaEizMID6SS7iz6RwCSMQg02RpyFlo03C0ng8CQSaIMknkaZMmdK1a9dqKVhYWFh1EYWenh6CEJqU7x9wsz1oA3V51wE1mOJG78EFhdAkSQ7CESNGPHr0qNrCx48fBwQESL8kgHp7UcyffsHN98QDeJuuYCf60ksuE3eZgaanHtMKKioqNDU1pVcKQIOtus/N9mAMcZKiCdNRI0Ps8WAmaIreGSPMyspKS0sjhAiFwri4uLdPfubl5e3atcvR0VHWBQLUJuk1fzqFiw/C6GBTN8WNnnyVXeCNR2NB0/JOEB44cCAkJET8evLkydWaNmvW7LfffpNRXQB19t0Dbq4nuoMKoJM5xVDkZhbfyRxRCE3IO0E4bNgwb29vQsiIESNWrFjh4eFR9SNjY2NnZ2c9PT1ZFwjwQUmv+bOYLKo4JrvRu+O4TuYYzYUm5J0gtLOzs7OzI4Ts3bu3U6dOxsbGcqoKoK6+jeLmeTKYLKooJrjSbkeEG4WMHg5doMmQPFlmwIABSEFo+pJe86Gp3FxP3EpGYZhoEn8r+lASpsxAE/KmRxgeHr5+/frJkyePGDFi1KhR+fn5Ev/BhQsXZFUbQC1WRXHz0R1UNFOa0yuj2M9a4PAFmoo330WWZSsqKkQiESGksrKyogbyKxXgHQmF/IU0bp4X9qcKprcN9bKEPMzHPbihqXjTIwwMDAwMDBS/PnbsmJzqAair5fe5z70w1KR4GIpMbE7tjuM2dsCUGWgScDQNCik6j7+awc/D6KBimtyc/jORq8RAITQNkvcjR48ePX36tPh1Xl7eyJEjHRwcgoKCXr16JcPaAGq09C67tBWtjcfQKyYHXcrLkDqbgiSEJkFyEM6bN68q8xYvXnz69OkuXbrcunVr0qRJMqwNQLJ/M/nHBQSPoVdoE5rTexMwTAhNgoRdyevXrzMyMvz8/AghQqHw8OHDixcv3r9//9GjR8+dO5ebmyvzIgHe8XUku8qX1sAAkyIb7khfy8Q9uKFJkBCEJSUlhBB9fX1CyK1btwoLCwcOHEgIad26Nc/zKSkpMi4R4G3nUvlX5WS0M7qDik1bQNqwjpkAACAASURBVAba0X/hgkJoAiTsTczMzNTU1GJiYgghhw8fNjU1bdmyJSEkLy+PEIIHUIAc8YQsu8d+70fjQedKYEJzek8cghDkT8JkA4ZhgoKCpk6dGhAQcOzYsRkzZtA0TQiJjIxUU1MT34MNQC4OJXECmgy0R3dQGfSwpEpFJDqPb2mE4xqQJ8k7lJ07d44YMSI5OXnq1KnffvuteOHFixf79eunra1dl/WWlJTMnTvX09MzMDDw3r177zcoLy9funRp27Zt27Ztu3Tp0vLy8gb/DaAiWJ6sjOJ+aIveoJKgCBnjQu3FY+tB3iRPP9fV1d26dWu1hZs2bar7er/88svnz5+fPHkyIiIiMDAwOTlZV1f37QarVq36+++/9+3bRwiZOHEiz/Nr1qypZ/GgWnbHcbbapKcVclB5THClO58RrW3HqKGTD/IjlW9fcXHxvn371q1b5+rqOnPmTFdX10OHDlVr8+DBg5EjR7Zo0aJFixbBwcFRUVHSqASURrGQrLzPrWmLqaJKxVmPaq5PnU9FpxDkSXKPUCgU/v777ydPnnzx4kVZ2TsTnJOSkmpd6bNnz3ie9/T0FL9t27ZtbGxstTYjRozYtm1bQEAAIeTAgQOzZs1qSPmgMtZEs71tKF8TdAeVzXhXem8CP9Be3nWACpMchJ999tnevXtbtWrVpk0bLS2t+q701atX4qsvxAwMDBISEqq1GTp06MGDBz/55BNCSKtWrYYOHVrT2mJiYm7durVz507xW5qmw8PDzczMJDYWCoUsy4pvHQ6yVFpayrKseF5Vo0svJTsfq18PrCwulsbqFVtJSQlFKfDxQT9z6stbai9yy401FOz6ekX/5BWUSCQS7+fr0ri8vJzna/9eSQhCoVB48ODB5cuXr1y5st41EkII0dfXF1+MKFZcXGxoaFitzfjx452dnc+dO0cImT9//rhx406dOiVxbZ6ent7e3tOnTxe/1dTUtLKyqulXiz8gXOMhexRFaWlpSSkIV95l53tTLczwvCUJeJ7X0dGRdxUNp0NIPzv2TJZgjoeCjRMq+ievoMRBWMcemqamZl0OViR88/Ly8ioqKgYPHlzvAv/H3t6+pKQkIyND/DY+Pt7BwaFam+vXrw8fPpymaZqmhw8ffuPGjZrWxjCMkZGR0/98IAVBKd3O5q9l8gu8FWwvCXU3oTmNuaMgRxJ2Lqampvb29vHx8Q1eqampae/evX/++WdCSGxs7JUrV0aNGkUISUxMXL58ubiNl5fX8ePHeZ7nef748eNeXl4N/nWgxHhCQm6z3/vh/trKrJcVlVlGHuEJhSAnEoKQpuldu3atWLEiMjKywevdsmVLWFiYlZVV9+7df/75Z1tbW0JIenr67t27xQ22b99+584dKysrKyurW7dubd++vcG/C5TYoSSuTETGuqA7qMxoiox1ofYnolMI8iH5MHv58uXp6el+fn4mJiZ6enpv/6gus0YJIc7Ozg8ePMjLy9PR0VFX/+/QTvfu3dPT08Wv3d3db9++XVpaSghp1qxZw/8CUF7lLFlyj/ujO0NjRoKyG+dK9z7PfudHcLsEkD3JQditW7dWrVp9/NqNjIw+3AARCB+wIZbzM6G6WmDXqPw8DCgzTXI1k/e3xH83yJrkIPzpp59kXAdANWkl/MaH7O1BGBtUFeNc6X0JnL8l7pkAsoahF2ii5t/k5ngwTrroH6iKMc70qRdcKa4BBpmrMQhPnTrVpUsXIyMjGxsb8ZIff/xRPBEUQNoupPGx+fyiljhQUyFmWqSdKXUmBVNmQNYk72j27t07ePBgLS2tQYMGVS20sLBYs2ZNHa/nB2iwMhGZfZ3d3onRxEkyFTPWhd6XgCAEWZMQhDzPL126dP78+eHh4RMnTqxa3rlz5+zs7KppnwBSsiqK7WhO9bLGSVGVM8SBvpHNZ5XV3hKgEUkIwqysrJcvX06aNKnacgsLC0JIdna2LOoCVfUon/8tjvupHTqDqqiZgAywow89Q6cQZEpCEIov+6v20AlCyPPnzwkhb99NG6Bx8YTMucF+68tY4rIaVTXWhcaV9SBjEoLQyMjIw8Nj+/btPM9X3a6U5/m1a9fa2Ni4uLjItkJQIXviuBIRmdYCc2RU1ydWVEYpeVqA262B7Ei+SGvt2rWDBg1KT0/38PAoKyvbsmXLkSNHrl27tm/fPjx2BKTkVTn5zz32QqAA95FRZTRFRjlT+xO57/xwehxkRPKhd//+/U+dOvXy5cutW7fm5eXNmzcvOTl53759Y8eOlXF9oDpm/MtObE63MkYMqrqxLvSfSXV4iBxAI6nxth39+/fv379/WlpaVlaWnp6ei4sL+oIgPXsTuKcF/J/+uI8MEB8jSk+NXMvku+HueiATtex3bGxsqi6oB5CS1BJ+4R32Yl8BLhwEMfEFhd0s8IUAWZBwavTly5fLli3r2rWro6Ojg4NDhw4dFixYUMeHTgDUF0/IZ9fYEC+mpREO/+G/xrhQJ55zFbh7B8hE9SC8evWqh4fH999/HxMTo6enZ2homJCQsGHDBk9Pz6NHj8qlRFBumx9yRULypQ9misIbVs0obyPqXCquowBZeGfvU1hYOGLECD09vdDQ0Pz8/Ojo6KioqJycnKtXr7q4uEycOPHFixfyKhSU0tMC/rsH7O/dGDyFDqoZ40L/mYQZMyAL7wThgQMHcnNzz58//+mnn9L0f39EUVTXrl3Dw8MFAsEvv/wijyJBOYk4MuEf9js/xlUfMQjVBTnSf6dzBZXyrgNUwDtBeOXKld69e3t6er7fztLSMjg4+PLly7IqDJTfivussSYunwfJ9NVJTyv6WDLOjoLUvbMPSkhI8PX1rampr69vfHy89EsClXA2hd+XyP/eTYDOINRkjAt1IAlBCFL3ThC+fv3awMCgpqaGhoaFhYXSLwmUX+JrfvJV0QF/xkxL3qVAE9bPln6Qy6eVYKQQpOudIBSJRB+4ap6maZEIT4+Gj1UmIiP+Zr/1YzqbozcIH6LBkCEO9KFnCEKQruoX1J89ezYrK0tiU5wXhUYx8zrraUhNx9Ag1MEYF3rBLXaBN74tIEXVg/DSpUuXLl2SSymgCjY95B7k8jcG4lZqUCfdLahX5eRhPu9liPMHIC3v7I+io6M5DkPTIC03s/k10eyNgYJmyEGoG5oiwU7UX0nc93gYBUjNOzskPHQXpCfxNT88gv29u8BJF4f2UA9jXOghEex3fgTfG5ASnHkHWXhVTvqFsSt96UAb7M2gfloZUzoCciMLU2ZAWhCEIHWvhaTPedF4V3qqG75v0BCjnOk/EzFqA9KCHRNIVwVLBl8UdTan/tMKXzZooDEu1JFkTogoBOnAvgmkiOPJuCusvjr1c0fMdICGs9eh3PSpsDScHQWpQBCCtPCEzLnB5lXwB3vi4RLwsca40LjdGkiJ5CDctGlTenq6jEsBZcLxZNZ1NjKHPxEg0EBvED5akCN9PpUrFsq7DlBGkoPwp59+sre3Hzx48Llz53BlIdRXJUfGXGGfFvDhnwp01eRdDSgFE03S2Zw69QK7I2h8koPw4cOH27dvT0pK6tevn52d3eLFi1NSUmRcGSioUhEZHC4qEZJzfQR6SEFoPKNxdhSkQ3IQGhgYTJs2LTY29t69e/369du8ebOjo2NAQMCRI0dYlpVxiaBACipJnwsiYw3qeC9GC7ePgUY12J6+kcVnlcm7DlA6tUyW8fX13bVrV0pKyty5cyMiIkaMGOHq6rp58+ayMnwZobrMMuIfKvI1of7owQgwDQsaWzMB6WdHH8WjeqGx1bK7Yln2zJkzkyZN2rp1q6Gh4fz58/38/BYsWNCpU6fy8nLZlAgK4eYruv1pboQT/XMHTBEFaRntjLOj0PhqDMLU1NQVK1Y4ODgMHDgwKyvr119/TU9P//nnnw8fPhwZGfnkyZOwsDBZFgpN2S9PubH/Mr90oZe0RE8QpKi3NZX0mk8uwgWF0JgkD+OMGjXqyJEjGhoao0aNmjFjhp+f39s/9fHxcXR0fPXqlUwqhCatWEimXmMTXvN/BwjdTTEqCNIloMkwR/qvJH5pK5x3gEYjec+VkZGxfv36CRMmGBgYSGxw9OhRc3NzaRYGCuBxAT88gu1qQd0YIBCVV8i7HFAJo53p6f+yS3HHPmg8koNw79695ubmmpqaby8sLy9/+fKlk5MTIcTT01MW1UFTJeTI2mhu8yN2XXtmvCtNCBHJuyRQEZ3MqRIReZDLtzJGpxAah+Sjqvbt29+/f7/awqioKGdnZ+mXBE3d/Ry+w2nRzWzu7mCBOAUBZIYiZJQT9RemzEDjqcdeTCgUqqnV4wLp4uLiu3fvvnz58gNtcnJy7t69++LFi7qvFuSoVEQW32UHXGTne9KhfQT2OjgkBzkY60r/lcRzmDEDjeSdICwvL8/Pz8/Pz+c4rqioKP8tycnJhw8ftrGxqeN6r1y54uTkFBIS0rJly3Xr1kls880337i4uMyePbt79+579uz52D8FpEnIkV+fcm5HRDnl5OEwdARBnjwMKCMNci0TSQiN450xwh07doSEhIhfBwYGVmtKUVRNkfa+uXPnrlmzZsqUKfHx8a1btx49erSVldXbDf766699+/Y9fPhQHK4lJSUN/AtAyjieHEjiVt7nHHTJ0U+Y9mboBYL8jXah/0ziulvihu7QCN4Jwl69eu3atYsQ8tVXX82YMePtEUEjIyMvL68WLVrUZaWPHz9OTEwcM2YMIaR58+YdO3Y8fvz4nDlz3m6za9eukJAQU1PTrKwsc3NzbW3tRvhroFHxhJx8zn0dyemrk1+7Mj0sEYHQVIxyplofZ7d0ZPBsE/h47wSht7e3t7c3IaSysnLIkCHW1tYNW+mLFy+srKyqJp06Ojq+f8/u+Pj4u3fvbty4kaIoXV3d48ePOzo6SlxbaWlpbm5uRESE+C3DMN27d6dpnJqToiIh2ZvAbX3E6aiRH9sxn9oiAqFpsdWmvI2o82ncYHvsCuBjSb58olrvrb5KS0vV1dWr3mppab1/5jM/Pz8lJeXJkydqamrTp0//4osvTp48KXFtqampsbGxSUlJ4rcMwzg6OpqYmEhsLBQKWZYViTCZv4GeFVO/J9F7k5h2JtxPrTl/C44QUlxc+z8sLS1lWRYHKLJXUlJCUap4pDLMhvnjKd3LWG6PKFTZT16+RCKReD9fl8bl5eU8X/tY8psgvHXr1r59+wYNGtS7d+8lS5a8fv1a4j/Ytm1brSu1sLDIy8urepuTk+Pu7v5+m2HDhonzctSoUcHBwTWtzc3NrUOHDosWLar195L/BWG1KyChVgWV5Phzbn8C96SAn9qCjhlGW2vXbwunKEpLSwtBKHs8z+vo6Mi7CjkY04Isixay6jr66rU3lgaV/eTlSxyEWlpadWmsqalZl4OVN0GYmpp69uxZT0/P3r17h4eH13QHtboEoaenZ1FRUUJCgqurK8/zN27cGDduXLU27dq1y8nJEb9+9eqVoaFhrauFRlfBktBU7s9E/u+X3CdW9FxP+lNbGoMuoBAMNUhPK/r4c25Scxx+wUd5E4RBQUFBQUHi1/fu3fuYlRoYGEyaNGnKlClff/31yZMndXR0+vTpQwg5e/bssmXLHjx4QAgJCQkZOHBg8+bN9fX1ly1bNmvWrI/5jVAvBZXkfCp36gV/MZ1rZUSNdaV3d1MzkNNhNUCDjXamdj5BEMLHktZdkjdu3Lhhw4YtW7Y4OjqGh4eLz5jZ2dkNGTJE3KB9+/YHDx787bffWJZduXLlqFGjpFQJVEl8zZ9L5U+/4O6+4ntY0QPtqM0d1czqdIIBoCnqb0dP/5d9WcpbNcNYHTTcmyAUCoV1edyunp5eXdarrq6+ePHiagt9fHx8fHyq3vr7+/v7+9etTmigYiG5nMFdSOPD0vhylgTaUHM96QBruhkeFAGKT5Mhg+3pg0l8iDeCEBruze5w//79kydPrvUf1GUGDsgXy5N7r/iIl3x4OheVw7czo/rY0CcDaC9D7CxA2Yx2oRfeYUO8cXYUGu5NEHbq1Gn79u1yLAU+UuJrPjydj0jnL2dwNtpUgDW1yIfpbkmh8wdKzN+SyiojTwp4dwMc50EDvdlHurm5ubm5ybEUaIC8CnLpJReezoen8xUsCbCmhjpQ2zqrWWDkD1QDTZFgJ+pAEvetL6Y7QwOhs6B4OJ7cz+UvpPLnUrlH+XwXCyrAmp7nSXvizCeopNHOdNDf7Cpfgg0AGuZNEEZGRh49ejQwMLB79+6rV68uKiqS+A/WrFkjq9rgHa+F5EIqdzaFD0vnTDSovrbUt35MF3MKl/2BimtjQmkw5GYW38kcUQgN8SYI4+LiduzYYWZm1r17999//z07O1viP0AQylhqCX/mBX/qBXf7Fd/FnBpgT3/rhwcBArxjjAt9IInrZI6jQmiIN0E4evTo0aNHi1/Hx8fLqR74r+Qi/kgyf/gZl1LM97OjZ7jTxwNobZzJBpBkrAvV9iS7sQOjhtmjUH/YszYtL4rf5N9QB3pde6arBcWg+wfwQfY6lLsBdT6VG4iHUUD91RiEIpHo3LlzUVFR6enpFhYWXl5eAwcOxM2speS1kBxL5n6P554U8EMd6B/aMt0tkX8A9TDWhd6fyA+0l3cdoIAkB2FaWtqAAQMePHjAMIyxsXF+fr5QKHR2dj5z5sz7z5GABuN48vdLfm8CF5rC+VvRId70p7Y0zu0ANECQE73wjrCgksFdc6G+JO90x48fn5aWdvjw4fLy8qysrPLy8gsXLnAcN3ToUI7jZFyiUnpVTtZGcy6HRUvvsh3MqMSRasd7MYPskYIADWSgTnpZ00eTsYOCepOw383Ly7t8+fLOnTuDgoIEAgEhhKbpPn36HDhw4OnTp0+ePJF5kUolMoef/i/rdkQYmcPv68HcHSyY40Eba8i7LADFN9aF2p+IIIR6k3BqVHw3UU9Pz2rLxUvQI2wYEUcOJ3PrYrgylsxoQT8bicceATSyvrb0Z9fYF8U8ri+CepEQhMbGxn5+fmfOnGnRosXby0+fPm1jY4MxwvoqFpLdcdzGh5yNNlnpS/e3o7GNAkiDOk2GO9IHkvglLbGRQT28CcKioqKqp9KvWrVq0qRJL168GDp0qIWFRU5Ozrlz5/bs2bN582bxyVKoi9wKsiGW/eUp18uaPtaL8TXBxgkgXeNc6SlX2SUtMdgO9fAm1Y4ePVrtMUzbtm3btm3b20tGjRoVHBwso9IUWUEl2RDLbn/MBTnRdwcJHHQRgQCy0MGMEnIkMofHcSfU3Zsg9Pf3P3z4sBxLUQ7FQrLtMbc+lg20oW8PEjjrYWsEkB2KkDEu1P5EztcEt1uDunoThA4ODg4ODvKrROEJObLzCffdAzbQhr45EBEIIB9jXeiuZ0Q/tWMEOD8KdYMBv8ZxPpVfcJu10yGX+wk88IBQAPlx0aMcdanwdL6vLbZEqJMagzA6Ovro0aPPnj2r9hiK8PBw6VelSB4X8Atusc+Lyfr2zKfY8ACagLEu9P5Erq8tzo5CnUgOwmPHjgUHB5uYmHAcp6Wlpa6u/uzZMx0dnbZt28q4vqasRESWR7L7E7n/tGJmuuOmMABNxUhn+j/3hIWVjD6u1oU6kLzzXrJkSf/+/VNSUj799NOxY8fGx8c/fPjQzs6uT58+Mq6vyQpL472PiV6VkUfD1OZ5IgUBmhBjDfKJNX0Et1uDupGw/y4tLU1MTFy4cKGamhohpLKykhDSokWLX3/9dfny5cXFxbKusYnJryDT/2VnXmd3dGb+6MGY4IEcAE3PeBdqbwKCEOpEQhCWlJTwPG9kZEQIMTIyys3NFS/38fEpKytT8Wf2HkziPI4KddRI7DBBHxuMCAI0UX1t6YRCPvE1L+9CQAFICEITExNtbe2UlBRCiKur66VLl8rKyggh165dI4SIA1IFvRaScVfYVVHcmT6C9e0ZPCweoClTo8loZxqdQqgLCUFIUZS/v//JkycJIaNGjSooKGjZsuWgQYOGDBnSqVMne3tVfPDlnVe87wkRRcjdwQI/3LECQBFMdqP3JfAc+oRQG8lzPHbu3Dl79mxCiL6+fkREhJ+f36tXryZOnHjy5EmKUq0YYHmyNpobHC5a157+owc6ggAKw8uQMtQgVzKQhFALyft1a2tra2tr8WtfX98DBw7IsKQmJKOUjLwk0mDIvcECq2aqdQQAoAQmutJ7E7ieVrigED7kQ7P+y8vLHz16dPHixejo6JKSEpnV1ETczObbnRL1tqbDApGCAApptAt9+gVXJJR3HdC0SQ5CjuNWrFhhZmbm5eXVp0+fVq1amZiYfP755xUVFTKuT17+SOAGh4u2dKKXtabx/EAABWWqSbpb0sdwQSF8kORTo1999dXGjRuDgoKGDx9ubm4ufh7h9u3bs7Ozlf40aSVHvrrNhqXxV/oJ3HHXUAAFN8GV2vyIm9gc97yAGkkIwrKysu3bty9btmzVqlVVC4cOHdq2bduZM2f+9NNPVcOHyiezjAyLEFloUXcHC3TV5F0NAHy0fnb09H/Z5CLeEY8FhRpIOEoqKCgoLy8fNWpUteXBwcE8z2dmZsqkMDmIK+Q7nxb1tqaP9mKQggDKQZ0mwc70HwmYOwo1khCEZmZmZmZmz549q7b82bNnGhoazs7OMilM1m5m8z3Oipa1pr9pgzFBAKUy0ZX+IwHXE0KNJAQhwzDr16+fM2fOlStXqhZGRkZOmDDh22+/NTAwkF11snLqBTc4XPR/3QWTMJAAoHTamFA6auQfXFAINXgzRnjmzJmVK1dWvS0oKPD39zc2NhZPlsnOztbW1j506NBXX30ljzqlaHcct+I+d66PwBe3jAFQUlPc6N1xXA9LXFAIErwJQl1dXScnp6q3b79WYt9Esgef8f/0Z5wwkA6gvMa60CvuC/MrGEMNeZcCTc+bIOzRo0ePHj3kV4kcLL7LXkjl/x0gMMWjlACUmpEGCbShDyRxsz0w/AHVqeh3gick5BZ7MY3/ux9SEEAlTHGjf3mKK+tBghqD8NGjRxMmTBDfU8bb23vUqFG3b9+WZWXSwxPyxS32aiYf8anAGOdJAFRDTyuqRETu52DKDFQnOQivXbvWtm3bI0eOWFlZDRw40MHB4dy5c507dz5+/LiM62t0PCHzbrB3X/GX+gmMkIIAKoMiZKIrvTsenUKoTvIt1ubOnevu7h4aGmphYSFeUlBQMHLkyNmzZw8aNIhhap95JRKJ9uzZc//+fXd39+nTp2tq1nj+cceOHUZGRiNHjmzYH1AvPCHTrrFxhfyFQNw4BkDlTGpOtTzOrmvHaOF5avAWCT3CvLy86OjodevWVaUgIcTAwGDLli2ZmZlPnz6ty3rnzJmzZ8+e9u3bh4aGBgcH19Rs3759S5cu/e233xpQegOE3GKfFPDnkYIAKslam+pgRh19jk4hvEPCcVFlZSUhRE9Pr9py8ZK6PIAiKyvr999/T0hIsLW1DQoKsrCwePz4sYeHR7VmmZmZP/zww+eff/7vv/82sPz6WHGfvfSSv9JPgIfrAqisKW70pkfcOBcVnScIEkm+xZqFhcWWLVuqLd+6daumpqabm1utK71z546Dg4OtrS0hREdHp127dhKjbs6cOStXrjQ0NGxQ5fWz9TF3IIkP6yvAVUQAqmyAHR1fyD8twJQZeENC54im6eXLl8+aNevx48cjR460tLR89erVmTNn/v7776+//lpbW7vWlWZmZpqamla9NTU1zcjIqNbm8OHD5eXlw4cP37x584fXlpycHBoaeufOHfFbdXX1devW1RSfQqGQZVmOe+fUx4Fk+scY5mIvoR5fWVpaa/nQEKWlpTzP0zQOtGWtrKysLsP2UCXYnv7tMbuq1ceeIMUnLxcikUgoFPJ8nQ5lKioq6tJS8lnCmTNnqqmprVq16ssvvxQvMTMz+/HHHxcsWFCX362uri4SiareCoVCdXX1txvk5uYuW7bs0qVLdVmbkZGRt7f3oEGDxG9pmjYyMqq2wio0TbMsq6Hxpt93KoVfHs2H96FcDdAZlCKRSKShoYEglL3Kysq3v/BQq8/c+Z7n+e/aqql93LcVn7xcMAxD03QdP3k1NTWKqv2uYRKCUCQSXblyJTAwcOrUqcnJyXl5efr6+k5OTnXfx1lbW6elpVW9TU1NHTx48NsNLl++nJOTI16YnZ1dWFjYuXPn69evS1ybvr6+g4NDHaeVivuCVYdpl17ys26IwvoKPI1wBzXpYhhG/AWVdyEqR/zJy7sKReJuRFz1RRdeUoPtP+rrik9eLnie5ziujp98HfdIEoIwOzs7ICDg6tWrNjY2jo6Ojo6O9SuTkC5dupSVlV27dq1r167x8fGPHj3q27cvISQlJSUlJaVLly59+/a9d++euPHevXv//vvvffv21fe31OpJAT/qsuiAv6AlUhAA3jK9Bb3jMfeRQQhKQ8L3wNjYWEtLq6SkpMEr1dTU/OGHH4YPHz569Gh/f/9ly5aJhwxDQ0PnzJlDCNHW1nb6H/Gva0DcflhmGfk0jP2pHfOJFVIQAN4R5ERH5/HxhZgyA4RIDEINDY1p06b9/PPP4usoGmbq1Kk3btwYPnx4RETE0qVLxQtHjhx55MiRai2Dg4O3b9/e4F8kUZmIDA4XTWlOj3fFER8AVKdOk0nN6V9x61EghNQ0WcbBweHw4cNubm6BgYE2NjYCwZtmixYtquOqnZ2dqz3O3sjIyMjIqFozMzMzMzOz+tRcC5YnYy+zngbUstZIQQCQbKY77XtStNKXaYYLi1We5K/ADz/8kJWVRQjZuXNntR/VPQjlZWEkXVDJH/oE324AqJGdDtXOlDqSzE3AeSOVJ/kbkJmZyddAxvXV15lUcjWLOtVboI7vNgB80Ex3ZscTnB0FSUGYk5Nz48aNZ8+eVbssXSEEWJErfVg9T8qaLAAAIABJREFU3EoUAGrT15bKKsODmeDdIBSJRBMmTDAzM+vcubOzs7OHh0d8fLy8KmsYTYbgVqIAUBcMRaa1oNEphHeCcPv27X/88Ye/v//q1atnzJiRnJw8adIkeVUGACBtn7nRx55z+bU/SgCU2Tu9p/Pnz3fr1u3vv/8Wv/X29p49e3ZhYaG+vr48agMAkC4TTRJoQ+9L5OZ5YlqB6nrn//758+e9e/euetunTx9CSHJysqyLAgCQlZnu9I4nHMYJVdk7QVhWVtasWbOqt+IHTZTieQ0AoLy6WlBqNLmSgShUXdUnlqSlpUVGRopf5+bmEkLi4uLevs+3r6+vzIoDAJCBme70tsecvyXuoK2iqgfhhg0bNmzY8PaSyZMnv/226V9KCABQLxNc6W8ihc+KaCdd3JpYFb0ThKtXry4uLpZXKQAActFMQCY2p7c95ta3R6dQFb0ThKNHj5ZXHQAAcjTfk255XLS8NaMv+ZnfoMwwYxgAgFhrU72s6d/jcXG9KkIQAgAQQsiXPvTPjzgWsyBUD4IQAIAQQvxMKEstciYFnUKVgyAEAPivz73ojbEIQpWDIAQA+K9hjnRqCbn7CqdHVQuCEADgvxiKzPKgNz9Cp1C1IAgBAN6Y1oI+n8qllaBTqEIQhAAAb+ipkdEu9E48pFCVIAgBAN4xz5P+NY4rEcm7DpAVBCEAwDtc9KgelvQvT9EpVBUIQgCA6pa2on+KYctZedcBMoEgBACorqUR1caY+iMBnUKVgCAEAJDg69bM2mhOhChUAQhCAAAJ2ptRdjrk0DMkofJDEAIASPafVsz3DzgOlxQqOwQhAIBkvawpA3VyGrfhVnYIQgCAGi1qSa+6jz6hkkMQAgDUaKA9LeJJRDqiUJkhCAEAakQRsrglvfoBrihUZghCAIAPGelEp5WQq5noFCotBCEAwIcwFFnehl56F51CpYUgBACoxRhnulhIzqWiU6icEIQAALWgKbLSl158l8X8UaWEIAQAqN0ge1pHQA7jRjPKCEEIAFAnP7Rjlt/nhIhCpYMgBACok24WlL0O+T0eSahsEIQAAHW1pi2zMoorxcPrlQuCEACgrvxMqPam1I4n6BQqFQQhAEA9rG5L/xTDvhbKuw5oPFIMwu+//97c3NzIyGju3LkiUfVTCTdu3AgICDAyMjIxMRk9evSrV6+kVwkAQGNx06f62tI/xeD6euUhrSC8cOHCjh07bt26lZCQcP369Z07d1ZrkJubO23atMTExCdPnhQWFs6ePVtKlQAANK7vfOmdT7gXJZS8C4HGIa0g3L1799SpUx0dHY2NjUNCQvbs2VOtwYABA4KCgoyMjExNTadNm3bv3j0pVQIA0ListakvvJj/PFCTdyHQOKQVhHFxcV5eXuLX3t7ecXFxH2gcFhbWvn17KVUCANDovvShHxdSF9JwpxllIJDSevPz83V1dcWv9fT0SktLy8rKtLS03m954sSJI0eO3L9/v6ZVRUVFXbx4cfHixeK3WlpasbGxZmZmEhsLhUKWZYVCDGTLWmlpqUgkomlMv5K14uJieZegopa5VX5+Q/dGnwo1fOtlSCQSCYXC9+edSFRWVsZxtU/xlVYQGhsbv379Wvy6sLBQR0dHYgqGhYXNmDHj/Pnztra2Na2qdevWPXv2XLRoUV1+rzgINTU1G1Y2NBhN01paWghCuag66ARZGupcdDiL3puq/YUXvvayIw5CiYHyvjrulKT1/+fm5hYbGyt+HRsb6+bm9n6ba9eujR8//siRI35+flIqAwBAejZ1ZNY8YDNK5V0HfBxpBeGUKVN+++23uLi4zMzMdevWTZkyRbx85syZt27dIoRcv369b9++S5Ys0dbWjoyMjIqKklIlAABS4qJHTWxOfx2JSykUm7ROjfbu3XvBggUBAQEikWjMmDHTpk0TL09JSREPaURHR7do0WL//v379+8nhGhoaFy/fl1KxQAASMnXrRn3o6I7r/h2priaQlFJKwgJIQsWLFiwYEG1haGhoeIXs2bNmjVrlvR+OwCADOiqkdV+9Jwb7M2BAgZRqJgwxgsA8FHGudKG6mR9LG5AqqgQhAAAH4Ui5JeuzE8x7KN8XFaokBCEAAAfy16HWtmGmXKNZRGFCghBCADQCGZ60HpqZONDnCBVPAhCAIBGQBHySxdmbTT7uAC9QgWDIAQAaBwOutTy1syUqzhBqmAQhAAAjWa2B63BkM2PcIJUkSAIAQAaDU2R3V2Z1Q9wglSRIAgBABqTsx61qSMzNJwtwlNwFASCEACgkY12pjuZU9P/xT1IFQOCEACg8W3vzDwt4H+Lw2ChAkAQAgA0Pk2GHOrJ/Oceez8Hg4VNHYIQAEAqXPWpLR2ZkZfYwkp5lwIfhCAEAJCWEU50Hxtq/D+4sLBJQxACAEjR+vZMRin/YzQGC5suBCEAgBRpMORkALPjCbcvEVnYRCEIAQCky6oZdS6Q+fI2G5GOU6RNEYIQAEDqPAyowz0Foy+LYvKQhU0OghAAQBa6W1JbOjH9w9jUEmRh04IgBACQkZFO9GwP+tMLbAEuqGhKEIQAALKzqCXdzZIaEi4qEcm7FPgfBCEAgExt7cR4GlJ9zotwoX0TgSAEAJApipAtnZi2plTPc6KccnlXAwhCAADZowjZ2IHpZ0t1PyvKKJV3NSoPQQgAIB+rfJnxrrR/qCgd80jlCkEIACA3i1rSk5rTPULZJ3iivfwgCAEA5GlRS3pZa7pHqOj0C9yDTT4E8i4AAEDVTXClfYyooRHsjWx+tR9DU/IuSMWgRwgAIH+tjak7gwS3s/lB4bisQtYQhAAATYKpJgnrK7BqRnU+I3qMIUMZQhACADQV6jTZ1YUJ8aZ7nBWtjebwPF/ZQBACADQtk5vTdwcLIl5yvidED3IRhlKHIAQAaHLsdaiLfQWzPOhPzokW32WFmE8qTQhCAICmiCJkWgv6/hBBZA7f9qQoLA1dQ2lBEAIANF32OlR4X8G69sxXd9he50SROYjDxocgBABo6npZU/cHC4Kc6IEX2dGX2WdFiMPGhCAEAFAAAppMb0HHjxC4G1DtToqCL7G3shGHjQNBCACgMLQF5OvW9PNgtZ5W1OSrrN9J0R8JHKbSfCQEIQCAgtFRI9Na0A+HCb5pQ/+RwDkeEi26w+JCiwbDvUYBABQSTZEBdvQAO/pRPn8giRsWwaozZKQTFexEtzDA7UrrAUEIAKDYPA2p7/2Y7/3I7Wz+4DOu13lWT40EWFMB1nQPS0pHTd71NXlSDMLk5OTbt29bW1t37dpVYoOioqK///6bpulevXo1a9ZMepUAAKiC9mZUezNmfXsSlcuHp/M/P2RHXeZbG1OfWNHtTKm2ppSJprxLbJKkFYSnTp2aMmVKv3797ty506FDh//7v/+r1iAjI6Njx47e3t7/396dhzVx5g8Af3PInZADkIRL5CqKi4KKbo3AioDuakVRFLQKUlDbxwrVrdp6bR9P1KqtruL1iKuouKtVcOUsggp4IVQ8oARiOQLBQA6OkGTm98e7nV82HAYVWMj7+eudN9+ZeWdyfDMz78yrUqk2bNhw//59NpvdT41BEATRH2QS8LYgeVuQNnqS21QgX4jnCbHvn2EPRTjTkDTJkjTRguTBJH3EAA5mJDTkE+inRIjj+MaNG3/44YclS5aIxWJnZ+eSkhJPT0/NmEOHDk2dOjU5ORkA8Mknnxw7dmzLli390RgEQRC9ZUIFQbakIFsKAAAHoEKCPxThj5vw7DrslQQ0tONu5iQ3c5ITHdiZkuzNSA5mwN6MRNOzs6n9kggrKioqKytDQkIAACwWKyAg4MaNG1qJ8MaNG7t27YLlhQsXHj58GCVCBEGQ/kMCwNWc5GpOinD+T02bCrxswV9JcL4MFL/BfxJgr+VAIMfJJMAxIVkaAUsjkrUJsDICDEMSwwAwDP5TMDcAplSSCRUMj5TZL4mwtrbWwsLCyOg/Z6NtbW1ra2u7xtja2sKynZ1d1wCCWCx+9erV8ePH4SSFQomIiDA0NOw2WP27990GpI/gbsdx1IF7oKEP/GAZBnvekAQ8mcCTSVTA86QkSScQtuNNHaTGDry+DTR1AIEML+kELQrQ0om3dAKZkiRX4e0qIFMC2ghAJQOmAQkAYG6AkwCgjwAUMgAAMAxIcInmBoA4B2tCBYaU/2oGGQC6zgkVwzCeJfaxjU57HsN0usWyXxKhUqmkUP5/Q6lUamen9ojLKpWKiKFQKF0DCFKptLa29tGjR0TwnDlzGAxGT6tWq9Waa0cGhlKppFKpZDK6M3WgKZVKpVI52K3QR8N4z5uQwGgTMFq3LoxSJVDjpJZOHAAg6SThsAYDAICW33dPiwIQ/5Hb1STFf2cxDIA3Hbr+h8YwTNyO6bjndfx33i+JkMPhvHnzBsMw+LPY0NBgZ2fXNUYkEsFyY2Mjl8vtaWmjRo2aMGHC119/rcuqKRSKWq0mDkaRAQN3O0qEA0+pVKIP/KBAex6C+4AzUKtTqVRKJdBxz48YMYJEent3oH752XJzc6PT6ffu3QMAqFSq3Nzc6dOnAwDUanVrayuM4fF4mZmZsJyVlQUDEARBEGSA9csRoYGBwfr166OiouLj4zMzM62trQMCAgAAGRkZ4eHhzc3NAID4+Hgej8dgMFQqVXJyclFRUX+0BEEQBEF6118nstavX5+QkFBeXs7j8eBd8wCAsWPHJiQkwIA//OEP9+7da2trU6lUhYWFrq6uH2S9AoGgrKzsgywK6ZOioqI3b94Mdiv0UXp6OuqjNPAUCkVOTs5gt0If1dTUlJaWfthl9uOTZebNmzdv3jzNGnt7++joaGLSw8ODuIPiQ/npp58EAoG3t/eHXSzyVocPH16xYgW8ZwYZSKtWrfLz87OwsBjshuiXysrKjRs3fvLJJ4PdEL1z+/btR48e+fj4fMBlDreuDeiv8SBCOx9BkP7WH78zwy0RIgiCIEifoESIIAiC6DXS//7prNDQ0NLSUgcHB12CX79+3dHR8aG63iC6Kykp4XK5lpaWg90QvZOXlzd16tQRI4bFo66GjtbW1tLS0qlTpw52Q/ROTU2NTCZzd3fXJRjerf706dPew4ZAInz69GllZaW5ubkuwVKptLOzE3UcGHh1dXVsNrunR98h/aeqqsrR0XGwW6F3cBwXCASjRo0a7IboHblc3t7eruN/boVCYWJi4u/v33vYEEiECIIgCNJ/0DVCBEEQRK+hRIggCILoNZQIEQRBEL2GEiGCIAii1/rxEWv95MWLF4WFhaNHj/b19e02QCwWp6enGxgYBAcHm5qaEvUPHjwoKyvz9PT08vIaqMYOK48fPy4tLR07duzkyZO7vtrW1nb37l2hUOjs7PzHP/4RVnZ0dNy9e5eIcXZ2Rr3s+kqhUNy+fVsmkwUEBFhbW2u9imGY5hMv7e3tiXuHpFLp7du3AQDBwcF0On3AGjxsyGSy9PR0lUoVHBzcdQBUgUBQUVGhWcPj8QwNDV+8eEEMM04ikWbMmDFAzR0ucByvrKysrq6eMmWKmZlZtzF8Pj8vL4/L5QYEBBBDv6nV6qysLKFQ6Ovr2+ffGXxIOX/+vKWlZWxsrLu7e3R0dNcAPp8/cuTIsLCw2bNnu7m5icViWL9t2zYHB4fY2FhbW9t9+/YNbKuHg71799ra2q5atcrBwWHHjh1ar6pUKhqN5ufnt2LFitGjR8+dO1elUuE4Xl1dTaFQAn6XnJw8GG0fwtra2iZOnOjr67t8+XI2m11cXKwV0NHRAQDw9/eHe/jYsWOwvq6uzt7ePiQkJCQkxN7evr6+fsDbPrSJRCInJ6e//OUvCxcu5HK5r1+/1gpISUkhPtgeHh6mpqatra04jsfExDg7O8P6wMDAwWj7ECaXyxkMBpvNJpFIT58+7Tbm1q1bbDY7Ojray8tr7ty5sBLDsODgYG9v75UrV7LZ7MzMzD6tdyglQpVKZW9vn5qaiuN4U1MTjUZ7/vy5Vszq1atjYmJgOTg4eO/evTiOi0QiY2Pj8vJyHMdLS0vNzMwkEsnAtn1oa2lpMTMzKy0txXG8vLzcxMSkqalJMwDDsMrKSlhubm42NzfPycnBcby6uppGow18g4eNs2fPent7w38VW7dunT9/vlYATIRyuVyrftOmTWFhYbC8ePHiTZs2DUBrh5MdO3bMmTMHllesWLFu3bpegmNjY1esWAHLMTExu3fv7vf2DVMqlaqqqgrHcUNDw54S4YQJE86cOYPjuFwut7GxuXPnDo7jmZmZ9vb28L/IiRMnpkyZ0qf1DqVrhKWlpWKxOCgoCADAZrN9fX3T0tK0YlJTUxcsWADLCxYsSE1NBQBkZWW5urq6uLgAAMaNG8flcnNzcwe06UNcbm4ul8sdN24cAMDFxcXFxSUrK0szgEQijR49GpbNzc2NjY07OzvhJDx3d/fuXblcPsDNHgZSU1PnzZtHoVAAAKGhoWlpaXh3N/7eu3cvNze3paWFqLl582bXLwKiu9TU1NDQUFgODQ3tZQe2t7dfvnw5KiqKqKmpqfn3v/9dXl7e760cdigUSu9nNWtqaoqLi+Fn29TUNDg4GL41qamps2fPNjExAQCEhoYWFhbCZ8roaCglwtraWmtrayr1P9c1bWxsiHPxEIZhQqHQ1tZWK6C2tpao7HZGpHd92oEnT56k0WjTp0+HkxYWFocOHfryyy+dnJzu3LnT720dXmpra21sbGDZxsZGoVA0NTVpxVhbWx85cmTz5s2Ojo7Xrl3rdkb0ge8r3XdgSkqKlZXVtGnT4CSVSn369OmxY8emTp26aNEitVo9EM3VG3V1dXQ6nbjmrfkjT7xfLBbL2Ni4T5/5odRZRq1Wk0gkYpJCoahUKs0AHMcxDCNiiACtGalUqtaMSO9034FZWVlbtmxJS0szNjYGANjZ2VVVVcF59+3bFxUVVVlZOTBtHh7UajXRFwAeF2rteUNDw5qaGvhSUlJSZGTkn//8ZwMDA823rOs3BXkrrT2vVqtxHNf8FhDOnDkTFRVFvHTkyBH4dojFYm9v73PnzmkeLCLvqacsoPl+gb5/5ofSESGHwxGJRMSpoYaGBg6HoxlAoVAsLS2JI+KGhgYulwtnbGxsJMKIekRHOu7A/Pz88PDwq1evTpw4EdaQyWTiU7tkyRI+ny+RSAagwcOG5p5vaGigUqlWVlZaMfBnFwCwePFiqVRaVVUFfv+yEDOiD3xfae15DofTbRasqqq6f//+smXLiBri7WCxWEFBQcXFxQPQWv1hbW0tlUrhpXGgkQU036/W1la5XN6nz/xQSoTjxo0jk8kPHjwAACgUijt37sBHqXZ2dhJXR/z8/NLT02E5IyPDz88PAMDj8UpKSuDvQm1tbXl5+ccffzwomzBEffzxx+Xl5fBUg0gkKikpgWc+Ozo6pFIpjCkoKAgNDb106RKPx+t2IU+ePGEwGKgff5/4+fllZGTAcnp6Oo/Hg7+zEolEoVBoBRcXF5PJZHiCyN/fn5iR+CIguvP39+/6SwIAaG5uViqVRNipU6eCg4O7/c3FMOzp06d2dnb939jhr62tDXYycHBwcHR0zMzMBACo1ers7GyYBfz8/DIzMzEMAwBkZGS4ublpHSa9xXt08BkEO3fudHZ2PnjwYGBg4IwZM2Dl1atXuVwuLD969IhOp2/bti0+Pp7NZhOdnpcvX+7j43Po0CEvL681a9YMTuuHslWrVnl5eR06dMjHxycyMhJWHjx40MfHB8dxmUxGp9M9PT1jfpebm4vj+OHDh6Oiovbu3RsXF2dubv7jjz8O5jYMQWKx2NbW9rPPPtu9ezeDwcjIyID1EyZMOHr0KI7jFy5ciIiI2LNnz9dff21hYbFlyxYY8PLlS3Nz802bNm3evJnBYLx8+XLQtmFo4vP5TCZzw4YNW7duNTc3LykpgfVsNvvGjRuwrFKp7Ozsrl+/rjkjj8f79ttvd+/e7e/vP3r0aOIOLkRHGzdujImJoVAoCxYsiImJgT38v/rqK6LL9KlTp7hc7oEDB+bPnz9+/HilUonjeGdnp4eHR2ho6P79+zkcTlJSUp9WOvRGn7hx40ZBQcGoUaOWL19uZGQEAKiurr5//354eDgMKCsrS0lJMTAwiIiIIEYxVKvVFy5ceP78uaenZ1hYmObZZEQXGIZdunQJ3lAfHh4Oj0uePXvG5/Pnzp3b0dGRlJSkGT9t2rQxY8YIBIK0tLTffvuNyWQGBASgRxm8A6FQmJSUJJfL586dS5xzvnbt2kcffeTu7i4UCm/evFlVVUWj0Xg8HtFlAwBQUVGRnJxMIpGWLFni7Ow8SM0fwqqqqi5cuKBWq8PCwj766CNYef78eV9fX3t7ewCAWCy+evVqZGSk5mCQaWlpjx8/7uzsdHJyCgsLg/0YEd1dvHhRs4f50qVLTUxMioqK4DMlYGV2dnZOTs7IkSMjIyNpNBqslEgkZ8+eFYlEM2fO7OspkKGXCBEEQRDkA0IHRgiCIIheQ4kQQRAE0WsoESIIgiB6DSVCBEEQRK+hRIggCILoNZQIEQRBEL2GEiEyzGVnZxOPCBkAOI5fu3btnR9yLZfLExMTtUZ81TTAmzMALl269OTJk3eYsbCwED5nCkHeE0qEyNAWERExumfZ2dkHDx7ctWvXgLXn8uXLn332GXGTb181NzfHxsYWFRX1FHDgwIHdu3fDskgkSkxMHOojS6xZs+by5cuwXFZWlpiYqPkMs17U1dUFBwe/efOmP1uH6IWhNPoEgnQ1a9Ysd3d3WM7MzMzLy9u4caOpqSmssbe3//TTT4nBEfubSqXatGnTunXr+u+Rqp9++imRJ6qrq2NjY7OysogBaIaib775ZsKECbCck5Ozdu3a8PBwzWe19CQkJGTbtm179uxJSEjo5zYiwxxKhMjQtnTpUqLc2tqal5cXFxenOUQDHJC5K7lcrlKpGAwGUSMWi2k0Wrc/wW1tbVKplM1m9/4DffPmTYFAoNkkSKFQSKVSS0vLnmZUq9UikYjJZGqtVCaTjRw5UrNy8eLFvTRAk0qlampqYrFYBgYGsEYikZDJ5F6OVtvb2yUSyciRI7sdaeH9SaVSCoVC/E2BvvrqK13mFYlEI0aM0Hy/SCTS8uXLd+7cuWPHDvQkM+R9oFOjyDC3bNmy+fPnw3JBQQGLxbpx48bs2bPpdDqTyQwKCmpubn706NG4cePYbDaNRtuwYQN8hj3066+/zp4929zcnMPhMJnM+Pj4Xk7cJSUlTZ48WXOI7TNnznh4eBgZGVlZWdHp9PDw8ObmZuJVX1/flStXHjp0yMrKisPhnDt3DtZLpdJFixbR6XRra+sxY8YUFhYSsyxduhQOz52dnT1jxgwAQEhICIvFYrFYOTk5AIBJkyZ98cUXe/bssbS05HA4V65cAQBERERYWlrC0T8cHBx++OEHYoGvXr1isVgXL15ctGgRjUbjcDipqaksFis3N1dz03bv3s3hcIhhXjSFhITAJhHmzZtH1Ny6dQu2bebMmebm5nQ6fcqUKdXV1UTwmDFjdu7cCQDYu3fvxo0bAQB2dnZwi+ADl3fs2MFkMq2srJhMJoPB+O6774h5FyxY0NLScvPmzZ7eEQTRBUqEyDAnFouJy0gqlaq5uXnNmjWTJ08uLCw8ffp0Xl7eypUrIyIi1q1b9/Dhw7i4uP3791+/fh3G19fX83g8oVD4008/lZWV7d+//+TJk3Fxcd2uCMOw/Pz8KVOmaFY2NDR8/vnnhYWFZWVl33//fXp6enR0NPGqVCpNTU09ffr0iRMn7t27R4xgtXXrVjMzs4cPH+bk5FCp1ODgYOJCILE548ePhxcLt2zZcuXKlStXrnh6esJlpqSkpKSknDlz5u7duz4+PnDDExMTnz179vDhwzlz5qxdu/af//wnXKBarW5ubo6Li6PRaBkZGWlpaTNnzmQwGImJiZqbduLECX9/f80DMoJcLtd8SrJWTWdnZ3Nzc2RkpJ+fX1FR0T/+8Y/y8vK1a9cSwUKhEA7mtXDhwsjISABAUlIS3CJTU9Pk5OTvvvtu+/btL168ePbs2cmTJ1ksFjGvo6OjtbX1zz//3O07giC6+qADaCDIYILHEw0NDZqVs2fPnj59Oizn5eUBAD7//HPiVXgaMzk5GU5iGObo6Lhs2TI4CUeP0lxgQkLCiBEjWlpauq4dHuXA0ZF6cuzYMTKZ3NbWBifHjx9vZGRUU1NDBLx+/RoA4OXlRdTw+Xwqlbp582Y4OWvWLF9fX1iGfSazsrI0V+Hq6mpmZqa1E7RMmzaNGNSmrKwMADBr1izNgF27dhkYGDQ2NsJJeMgFh9bqKiAgIDAwULNmxowZRM21a9cAAN9++y3x6rZt28hkskKhgJNMJvOvf/0rLB85cgQAIJPJiOC4uLhRo0b1si2+vr5TpkzpJQBB3gpdI0T0TlBQEFGGVxADAwPhJIlEcnV1/e233+Bkenq6s7NzaWkpEW9kZKRUKl++fAmPtDQ1NTUBADSPV6CioqLMzEyhUKhUKuvq6jAM4/P5Y8eOha9Onjy5a1cX4lwuAMDR0dHLy6tPA51PmzZNayD7zs7Oq1evPn/+HA5P3dTUJJPJNAPmzp2rORkdHb1jx46kpCR4AS8xMdHNzQ2OxvxuZs2aRZTHjBmDYVhtba2jo+NbZ5wwYcL3338fFha2bNkyPz8/MzMzrQAWi1VSUvLODUMQgDrLIHpIs08K7EiiVUMkicbGRplMtmjRIq3ZYc7TQqVSAQAqlUqzMjo6+ty5c/7+/i4uLkwmE54wlEgkRIBWXxjI2tpac5LD4VRWVuq6eV2WWV9fP23aNIlEEhQUZGVlZWhoaGRkpHW1T2uNlpaW8+fPT0xMjI+Pr62tvXXrVkJCwvv0oOm6z3XsyhsREdHQ0HDs2LErV64YGhoGBgYmJCS4ubkRAUqlUpemdzggAAAESElEQVQupgjSC5QIEaRHdDrdx8cnNTVVl2CYfjRva6uurj59+vTRo0fXrFkDa1JSUi5evKg5V7djRGslWpFIxOFwdG+21jLPnj1bX1//66+/crlcWFNeXq51iNk1ya1evXr69Ol37ty5c+cOlUrt2hWWQKVSFQqFZo1UKtXqAfvOyGTy+vXr169fX1FRkZ6evmvXruDg4MrKSmIbxWJxt38mEER3qLMMgvRo+vTp+fn59fX1ugRbW1vb2dn98ssvRE1VVRUAwNvbm6i5deuWLou6ffs2UW5sbHzy5ImHh0fXMHiesL29vfelVVVVcblcIgtKpdL8/Py3toHH440bN+748eNnz55duHBhL/d+2NjY8Pl8/PchvkUi0fPnz9+6/G7BWzu63SIXF5cvvvjim2++qa6ubmhogJUYhj179mzy5MnvtjoEgVAiRJAebd68GQAwZ86cvLy81tZWoVCYmZkZFRXVU/yMGTPu379PTLq7uxsYGCQkJIhEIrFYvG/fvn/961+6rPfBgwd/+9vfmpub+Xx+REQEjuOrVq3qGubg4GBmZnbq1Kn8/PzHjx/DvpddjR8/ns/nnz17VqFQvHr1KiwsrKOjQ5dmxMbGXrlyRSAQxMTE9BIWGBgoEAh27tzZ2NhYXFy8cOHCbg9zdQEvnR44cKCgoODx48dqtXrfvn3nz59//fo1hmGVlZVXr161sbEhTuT+8ssvUqk0ICDg3VaHIBBKhAjSIxcXl59//plMJvv6+pqZmXE4nDlz5vTyTK+VK1e+ePGC6FxjbW3997///fbt21ZWVmw2++LFi4cPH9Zlvdu3b7958yaLxXJycnry5Mnly5ddXV27hpmYmJw+ffrly5cBAQETJ04sKCjodmkxMTHz58+PiooyMjIaO3asvb19eHi4Ls1YtmyZsbGxu7s7cV9Ht0JDQ6Ojo7ds2TJy5MipU6fOmjXrnQ/RJk2atH379gsXLvB4vIkTJ8pkMrlcvnr1agcHBwqF4uzs/ObNm+vXrxMnci9dumRvb48SIfKeSMQJDQQZ6nAcxzCMQqFoVsK749/5GAWqqampq6uj0WijRo0yNjbuJdLHx2fSpEk//vgjUSOTySoqKszMzLpNZr149eqVTCaD9+MTle+8OTU1NfX19Y6OjhYWFjrOIhAInJycDhw48OWXX741WCQSVVdXu7i4dHuvYS8wDCORSL30xFEqlXw+XyKRwJPPRKRCoXBycoqPj4+Pj+/TGhFEC0qECPIhFRQU/OlPfyovL7ezsxvstryvmJiYlJQUgUDQf49OfR9Hjx49cODAixcvDA0NB7styNCGEiGCfGACgYDFYr3zABT/CzZv3pycnFxdXX38+PHY2NjBbk73hEIhlUrV/QAXQXqCEiGCINqys7MFAsH48eO9vLwGuy0I0u9QIkQQBEH0Guo1iiAIgug1lAgRBEEQvYYSIYIgCKLX/g8dv+rQu4mAfAAAAABJRU5ErkJggg==", "image/svg+xml": [ "\n", "\n" ], "text/html": [ "\n", "\n" ] }, "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 upper intercepts \n", "# for each sample given an estimated time of Pb-loss\n", "BootstrappedDistribution = BootstrapCrystDistributionKDE(smpl, tpbloss=0)\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, framestyle=:box)\n", "savefig(h, joinpath(smpl.Path,\"BootstrappedDistribution.pdf\"))\n", "display(h)" ] }, { "cell_type": "markdown", "id": "d9d6edda-09be-4bbd-ab82-835dc5a9471e", "metadata": {}, "source": [ "#### Run eruption/deposition age model" ] }, { "cell_type": "code", "execution_count": 9, "id": "c5b1a1ef-6c98-46d2-9dc4-a1fa92cf20a0", "metadata": {}, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ "\u001b[36m\u001b[1m[ \u001b[22m\u001b[39m\u001b[36m\u001b[1mInfo: \u001b[22m\u001b[39mEstimating eruption/deposition age distributions...\n", "\u001b[36m\u001b[1m┌ \u001b[22m\u001b[39m\u001b[36m\u001b[1mInfo: \u001b[22m\u001b[39m1: KR18-04\n", "\u001b[36m\u001b[1m│ \u001b[22m\u001b[39mInterpreting the five columns of KR18-04.csv as:\n", "\u001b[36m\u001b[1m└ \u001b[22m\u001b[39m | ²⁰⁷Pb/²³⁵U | 1-sigma 𝑎𝑏𝑠𝑜𝑙𝑢𝑡𝑒 | ²⁰⁶Pb/²³⁸U | 1-sigma 𝑎𝑏𝑠𝑜𝑙𝑢𝑡𝑒 | correlation coefficient |\n", "\u001b[36m\u001b[1m┌ \u001b[22m\u001b[39m\u001b[36m\u001b[1mInfo: \u001b[22m\u001b[39m2: KR18-01\n", "\u001b[36m\u001b[1m│ \u001b[22m\u001b[39mInterpreting the five columns of KR18-01.csv as:\n", "\u001b[36m\u001b[1m└ \u001b[22m\u001b[39m | ²⁰⁷Pb/²³⁵U | 1-sigma 𝑎𝑏𝑠𝑜𝑙𝑢𝑡𝑒 | ²⁰⁶Pb/²³⁸U | 1-sigma 𝑎𝑏𝑠𝑜𝑙𝑢𝑡𝑒 | correlation coefficient |\n", "\u001b[36m\u001b[1m┌ \u001b[22m\u001b[39m\u001b[36m\u001b[1mInfo: \u001b[22m\u001b[39m3: KR18-05\n", "\u001b[36m\u001b[1m│ \u001b[22m\u001b[39mInterpreting the five columns of KR18-05.csv as:\n", "\u001b[36m\u001b[1m└ \u001b[22m\u001b[39m | ²⁰⁷Pb/²³⁵U | 1-sigma 𝑎𝑏𝑠𝑜𝑙𝑢𝑡𝑒 | ²⁰⁶Pb/²³⁸U | 1-sigma 𝑎𝑏𝑠𝑜𝑙𝑢𝑡𝑒 | correlation coefficient |\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 closure/crystallization distribution to use\n", "dist = HalfNormalDistribution\n", "## You might alternatively consider:\n", "# dist = BootstrappedDistribution \n", "# dist = UniformDistribution # A reasonable default\n", "# dist = MeltsVolcanicZirconDistribution # A single magmatic pulse, truncated by eruption\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(joinpath(smpl.Path, \"results.csv\"), results, ',');" ] }, { "cell_type": "code", "execution_count": 10, "id": "59fcf977-4a2a-4abb-8687-46098df7fc65", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "AgeDepthModel.pdf\n", "AgeDepthModel.svg\n", "BootstrappedDistribution.pdf\n", "DepositionRateModelCI.pdf\n", "DepositionRateModelCI.svg\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", "KR18-01.csv\n", "KR18-01_Concordia.pdf\n", "KR18-01_Concordia.svg\n", "KR18-01_Pbloss.pdf\n", "KR18-01_Pbloss.svg\n", "KR18-01_distribution.pdf\n", "KR18-01_distribution.svg\n", "KR18-04.csv\n", "KR18-04_Concordia.pdf\n", "KR18-04_Concordia.svg\n", "KR18-04_Pbloss.pdf\n", "KR18-04_Pbloss.svg\n", "KR18-04_distribution.pdf\n", "KR18-04_distribution.svg\n", "KR18-05.csv\n", "KR18-05_Concordia.pdf\n", "KR18-05_Concordia.svg\n", "KR18-05_Pbloss.pdf\n", "KR18-05_Pbloss.svg\n", "KR18-05_distribution.pdf\n", "KR18-05_distribution.svg\n", "agedist.csv\n", "distresults.csv\n", "height.csv\n", "lldist.csv\n", "results.csv\n" ] }, { "data": { "text/plain": [ "4×5 Matrix{Any}:\n", " \"Sample\" \"Age\" \"2.5% CI\" \"97.5% CI\" \"sigma\"\n", " \"KR18-04\" 751.96 751.274 752.449 0.304138\n", " \"KR18-01\" 752.126 751.698 752.494 0.200616\n", " \"KR18-05\" 750.723 750.381 750.964 0.148518" ] }, "execution_count": 10, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# Let's see what that did\n", "run(`ls $(smpl.Path)`);\n", "results = readdlm(joinpath(smpl.Path, \"results.csv\"),',')" ] }, { "cell_type": "markdown", "id": "1a1d4317-4424-42bd-b047-0fa6331c9598", "metadata": {}, "source": [ "To see what one of these plots looks like, try pasting this into a markdown cell like the one below\n", "```\n", "