{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Pluton emplacement test\n", "\n", "This Jupyter notebook employs the stochastic sampling algorithm of [Keller, Schoene, and Samperton (2018)](https://doi.org/10.7185/geochemlet.1826) as implemented in [Chron.jl](https://github.com/brenhinkeller/Chron.jl).\n", "\n", "\"Launch \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 KernelDensity: kde\n", "using Plots; gr(); default(fmt = :png)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "***\n", "## Estimate the time of rheological lockup \n", "Given that the rheologically critical melt fraction ≈ 36% melt (see, e.g., [Caricchi 2007](https://doi.org/10.1016/j.epsl.2007.09.032) for experimental rheological data), which is also approximately the terminal porosity of a framework of close-packed spheres, a body of magma may be considered effectively \"emplaced\" by the time it cools to where F $\\leq$ 0.36. \n", "\n", "If we estimate the point within a MELTS zircon saturation distribution where melt fraction drops below 36%, we can use this to estimate the time by which a body of magma must have been emplaced.\n", "\n", "#### Input dataset (Try pasting in your own data here!)" ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [], "source": [ "# Age and one-sigma uncertainty.\n", "ages = [101.333, 102.175, 102.074, 102.814, 102.335, 101.212, 101.398, 101.772, 102.456, 102.302, 102.039, 101.579, 102.135, 101.978, 102.325, 102.36, 101.53, 102.409, 101.277, 102.766,]\n", "uncert = [0.137, 0.141, 0.134, 0.149, 0.132, 0.141, 0.129, 0.129, 0.135, 0.142, 0.144, 0.128, 0.15, 0.144, 0.132, 0.13, 0.133, 0.133, 0.148, 0.14,]\n", "age_unit=\"Ma\"\n", "\n", "# Sort by age (just to make rank-order plots prettier)\n", "t = sortperm(ages)\n", "ages = ages[t];\n", "uncert = uncert[t];" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Input MELTS-derived $\\ \\mathcal{\\vec{f}}_{xtal}(t_r)$" ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAlgAAAGQCAIAAAD9V4nPAAAABmJLR0QA/wD/AP+gvaeTAAAgAElEQVR4nOzdeXwM9/8H8M/MbO77TiQhkUMuonEECYJIlIbGrSjakNKiVVVUfX2rqlrU0Tb4Ube6676PuIVIELfEEUIuue/szOf3x7b5phHJYmc3u/t6PvyxmZndfe/andd+Zj7z+TCUUgIAAKCtWFUXAAAAoEoIQgAA0GoIQgAA0GoIQgAA0GoIQgAA0GoIQgAA0GoIQgAA0GoIQgAA0GoIQgAA0GoIQgAA0GqaEIQ5OTmpqamqrqJh4Xle1SVoOEqpIAiqrkLD4R1WAuwriGYE4Z49e2bOnKnqKhqWkpISVZeg4QRBKC8vV3UVGq6srAxZKDbsK4hmBCEAAMAbQxACAIBWQxACAIBWQxACAIBWQxACAIBWQxACAIBWQxACAIBWQxACAIBWk6i6AGi4eEoSX9DjafT4MyHxBeXp38ubGDO9GzMRTdhW1gyj0goBAN4eghBq97yE9DgklQok1JH5zIcNtGV1/zl8cDuP7kkVPozli6VkWRD3rjPSEEBdcRyn6hJUD0EItUgpoOGH+I+bsdP8azl43sGO6WDH/diGnE6nw2P5ga7MD204HRxlB1BDBgYGqi5B9bD3gpqu59DO+/kpLWpPweo62TMJkZJ7BSR4r/RhIa17YwBogBic30AQQg3XcmjYQemiduwYL7k+G1Z6ZFd3bogb23kfjywEAHWEQ6PwP7nlpN8xfnF7rr/ra/xCYgj53I/V40j3g/zp97hGhviBCQDqBC1C+Bsl5OMzfO/GzKCmb/KpGOvNfubDhh3kX2BuIgBQKwhC+NucRCGrjM5r++ZdyD73YyMaM+8ekhZVKrAuAABxIQiBEEKOptGY28KWrm/b+XNuG87HnBl7DnNeA4DaQBACKagkI05JN3VRzOm934O4xBf0j3uYWBwA1AOCEMh3CXxPZ7azg2I6uRhKyM5Qbuol/loOOpECgBpAEGq75AK69r4wu5UiR5fwNGPmteU+OMGXSBX4qAAAokAQarsJF/ip/pyDoYIfdpQn29qGmXABJwsB1FVhYSGltR/Xyc/PP3369OXLlxX1XCdPniwsLFTUo70uBKFW25sqPCwkE3xF+Rj81oE78YwefIIDpABiKSgosLS0bNOmTdWS7OxsW1vbwMBAQoitra11NYsXLz5z5oyvr69sy//85z/WL+nduzchZMGCBXZ2ds2aNTM2Nm7btm1ubm71J83NzW3RosXy5cuPHz8uZ53Hjh1bs2ZNHRuMHDnywYMHr/PSFQkX1GuvCoFMjhMWtxdrmFBjHbKqE/dhLJ/UT2KuK8pTAGg5Smlubm5eXl5cXJws/DZu3GhmZlZQUEAIyc7OPnv2rI+Pj2xjAwODs2fPvnjxQvbntGnTvvjiC0LI/PnzExMTN27cSAiRSCSxsbFz5849f/68p6enIAjHjh3T0dGp/qRxcXGOjo6y7eV0/fr1pKSkkSNHKuA1iwAtQu21/LbQzIzp4STiQDBdHJhezsyUSzhACiCiDz/8cO3atbLba9asqZ43JiYm5v/Q09Orfi99fX3Zcn19fYlEIrttbGx869YtT09PT09PQgjLsmFhYcbGxlX3io2N/frrr+/cudO9e/d169a9XMzZs2cHDx5cVlZGCLl9+3afPn3i4uKWL19++PDh7t27R0dHE0KioqJ8fHxcXFyGDBny/Plzxb8jrwlBqKV4ShbdEL5pKfoH4KdA7shTevgpDpCCJsstJ89KqBL+8bV9k4YOHfrXX3+VlpYmJSXxPN+yZcuqVbt27frjH0+fPpXntYSEhCQkJHz44Yd//vnns2fPaqzt0KHD5MmTfXx8tm7dOmDAgJfvHhwczLLs5MmTy8rKhgwZEhER0apVqw8//DAkJGTr1q0///wzIaRXr14XL168detWkyZNxo0bJ//7LBIcGtVS2x4IjY1JoK3o44Ka6pDlwdwn5/jrfSUmOvVvD6COpl3m96Qq49rZ8xESF5OaX1szM7POnTvv2bPnwoULH3/8cfVV169ff/Lkiex29VOJdfDx8bl06VJMTMx333139+7dvn37rl+/vmq2Jl1dXSMjIx0dHQsLi1c9QkxMTEBAwJUrV3x8fKKiogghBgYGenp6VXfp2bPn6dOn09LSLC0tT506Jd9LFxGCUEv9nKTgSybqEO7EdLZnfrjKz22DKUBBMy0L5pYRVX68R44cOX/+/Nu3b1+9evXSpUtVy2fOnNm8efOqP+Xs3tKiRYuYmBhCiOwQaExMzKRJk+QvxszMbPDgwT/88MOyZcteXpuenh4UFBQcHOzh4VFZWZmbm/uqvqlKg0Oj2uhYGq3giTJnlp/Xllt1V3hchAOkAKIICwu7e/duYGCgjY2NAh/Wy8urY8eOjx8/fq17Xb58eeXKlZMmTRozZkxFRQUhhOM4Qfi7xbxr167WrVuvXbt2xowZ7733ngKrfWNoEWqjn6/zk1uwypwtyc6AjPNhv4kXNoSgUQigeCzL3r59WyKpuUu/ePFiRkaG7LazszMhpKKi4tixY1UbdOzYscZdtm/fnpKSEh4ebm9vHxcXd+DAgfXr18tfSX5+/uDBgxctWjR48ODIyMjp06fPnz/f1dV15cqVcXFxJiYm1tbWSUlJT58+lUqlM2bMeJNXq2gIQq2TlENv5pEhbso+GPB1C85ru/RyFm1jgwkLARRDIpGEhobKLm+o6ttpbW3dvn17QkiPHj127NhRtXGfPn3at2/ftm3b+fPnVy1s2bKlm5sbz/+va7e3t/fRo0ejoqIKCwudnZ1jYmIiIiKqP6mtrW2rVq1eVdKGDRuGDx8+ZMgQQsiqVatGjRp1//79Xr16JSUlzZkzx8HBISYm5uLFi927d7eyspo6daqspRgSEmJqaqqQ9+QNMCo/OPv21qxZExsbW/fVmtqmsLDQxMSk1lXDY/kWlsxXLVRwVHzVXWHdfeHUe5rw84vn+YqKiqoeBCCGkpISfX19lsUZHBCXJuySQH5ZZWRfqvBrB9V03xzpyS65Kex6LLzfBLs2ALX37NmzuLi46ktMTU27deumqnreGIJQu2xKFiIas2YqGueFY8jPgdyE8/x7zqwEUQig5vLy8hISEqovsbOzQxBCQ7f2vjA/UJXdVcIcmUaGZF2y8JEnkhBAvfn4+MyePVvVVSiAknZGJSUlr1pVVlZW4zwlpbS0tLTGZjzPl5eXi1Kc1riRS3PKSYiC5h18Y9+35r5LEMox7BoANAyiB2FCQoKXl5eTk1OTJk1OnjxZfdWLFy/CwsLs7e2trKwWL15MCCkpKWnZsqWBgYG1tbWnp+e+fftkW3733XeWlpZ2dnaRkZFFRUVi16yp/rgrfOjBKPWyidp0sGN8LAimsAeABkLcIKSUDh8+fNy4cTk5OQsWLBgyZIjs4kqZb775xsLC4sWLFxcvXpw1a1ZSUpKOjs7KlStLSkqKi4u//fbbgQMHlpSUxMbGxsTE3Lx5Mysrq6SkZO7cuaLWrKmkAvkzRRjq3iAOSM5uxX2fKJRi2l4AaADE3S3Gx8enpaWNHTuWENK/f39jY+NDhw7JVkml0o0bN3711Vccx3l6ekZGRq5bt05HR6d169ay3tIRERGlpaW5ublr164dOnSok5OTjo7OpEmTqgZZh9dy4IngbsY0M1N1e5AQQkgra6atDbPsDhqFAKB64naWSU5Odnd3r5rLytvbOzk5WXY7IyOjqKjI29u7atWFCxdkt48cOfL48ePt27dPnDjR0dExJSUlKChItsrHxyctLa20tLTG9Vvl5eWZmZnVl5iamurr64v30tTO2vt0hEeDaA7KzG7Nhh6QRjVjMRI3AKiWuEGYn59vZGRU9aeJiUnVTMf5+fkMwxgaGsr+NDY2zsnJkd2+ePFicnLyo0ePQkJCCCF5eXlVIybIbuTl5dUIwiNHjlRNuyzz448/Dhw4UIwXpRZqnEnNqWCOp+kubVVRWNhQxk9oIiEdbXV/SSz9wlstj5DKLqiXStWyeHVRWlpaWVnZAC+of9VoFSAqSmlZWdnrfukMDQ05rp6u8uIGobW1tWyiZJnc3FxbW9uqVZTSwsJC2bA6eXl5dnZ2slUzZ84khGRlZbm7u3fr1s3a2jo/P1+2Ki8vj2VZa2vrGk8UERGBkWVqqP5dXXtLeK8JbWRhXMf2yvffNrTbfunkAANDNbyKByPLKAHHcRhZBqowDCPSN07cT5iPj8/du3eLi4sJIYIgXLt2zcfHR7bK2tra1ta26mLMxMTEqsOkMjY2Nra2tunp6T4+PomJibKFCQkJHh4eVcdaQU6bkoVhDaObTHU+5kyQPbsCZwoBVKewsFDVJaie6EHYqlWr6dOnp6enz50718TEpEuXLsnJyaGhoTzPR0VFzZw5MzU1dd++fYcOHRo1atS1a9e2bduWkpJy7969b7/9Ni8vr0OHDqNHj968efPJkycfPHjw/fffjxkzRtSaNc+TYppcQLs1ahDdZGqY0ZL96TpfhmsKAUB1RD8m9eeff06cOLFDhw7NmjXbs2dP9aMcM2fO/Prrr0NDQ83NzTdv3tykSZOSkpINGzbMmDGD4zh/f/+TJ09aW1tbW1svX7588uTJRUVF/fv3nzhxotg1a5g/U2hfF1anwTUICSGkpRXTyppZc0/4xLtB1gcAWgCzT2im6rNPtNolnR/IdVH1gDKvEpdJB57g7w+U6KpVFOIcoRJg9gklqGOmGu2BT5iGSymgz0toJ/sGmoKEkEBbxtOMbErGmUIAUA0EoYb7M4X2d2W5hpuDhBAyoyU395rAq/2xCQBQSwhCDbflgTC4aUP/X+7swFjqkb2paBQCgAo09F0kvI2kHJpfQdrbNez2ICGEkC/82EU3EIQAoAIIQk225YEwxI1RgxgkpJ8rm1pELmfh8CgAKBuCUJNtfUgHNfjjojIcQ8Z6s0tvolEIAMqmHntJeAPXcigvkABrtWgQEkLIaC92/xPh+SuncAYAEAWCUGPteUz7NFGbFCSEmOuSD9zY325hmBkAUCoEocbamypENFGz/9+JfuzyO0IJZnQAACVSsx0lyCmjjEkuoMHq0F+0OndTpr0tuxEX1wOAEiEINdPhZ2yYYwMdX7Run/uxS9BlBgCUSA33lCCHA8+43mp1grBK10YMw5DY57iOAgCUBEGogUql5EwGE+6krv+5Y73ZX2+hUQgASqKu+0qow4nn1N+SWumpuo43NdydPflMeFyERiEAKAOCUAPtfSy820iNW1TGOmS4B2auBwAlQRBqGkrI/if03UbqfTXeZz7syrsCZq4HACVAEGqahGyqzxF3E/U+ruhuygRYMVsfoFEIAKJDEGqavamCeg0o8yqf+XKLcR0FAIgPQahpTjyj6ttftLp3nZj8CnIxU72btgDQ8GnCHhOqlPMkIZu2s9WEFiHLkE99cHE9AIgOQahRrmRTL3PGREfVdShIVDP28FPhaTEahQAgIgShRjmbQYPUbXzROpjokCFuuI4CAMSFINQo5zQrCAkhE33ZFXdwHQUAiAhBqDkoIRcyhA6aFYQeZkxLK2YLrqMAANEgCDXHvXxqpMM4GWlUEBJCJvhyi28gCAFALAhCzXE2XdOOi8q868yUSMm5DHSZAQBRIAg1x/kM2kEjLpyogSFkHK6jAADRIAg1x9kMGmyvgUFICBnlyZ54JqRiPgoAEAGCUENkl5GMUuproZlBaKJDPvTAJIUAIAoEoYY4myG0t2U4zcxBQgiZ4MuuvicUVaq6DgDQOAhCDXEunQbZafL/ZhNjJsSBXXsfjUIAUDBN3nVqlXOae4Kwyud+7OKbgoAThQCgUAhCTVDOk6Rc2tZGw4MwyI6x1CP7n6BRCACKhCDUBJezqLc5YyhRdR3im+DLLsLF9QCgUAhCTRCXpSFTL9VrgCt7L59cz8HhUQBQGAShJriYqS1BqMOST31YjLgGAAqEINQE2hOEhJBoL3bXYyGjVNV1AICmQBCqvWcltEKgTU20JQgt9Eg/V/b/MEkhACgIglDtnc+g7Wy16/9xoi8bc1uoQBQCgCJo1w5UI8Vl0kBNv3CiBl8LxseCbH+IJAQABVBGj/uioqIjR45IpdLw8HAzM7PqqwRBOHHixLNnzzp16uTi4kII4Xn+0qVLKSkptra2Xbp00dHRKS8vP3PmTNVd3NzcXF1dlVC2uriYRf8boHU/aCb6crMS+A/ctO6FA4DCiR6E2dnZ7dq18/Ly0tfXnzRp0oULF5ydnavWRkZGPnny5J133pk0adLGjRvDw8O7dOlSUFDg7+9/48aN8vLy06dPl5WVhYeHd+3aVXaXkSNHIgirVArk6gvaWstahISQXo2ZyXHa1UsIAEQiehDGxMT4+Pjs2bOHEDJy5Mhffvll4cKFslWnTp2Kj4+/d++ekZFRUFDQN998Ex4evmrVKg8PD0KIIAjt2rVbtWrV0KFDdXV1jx49Knap6igph7oYM6Y6qq5D6WSTFC6+KbSz5VRdCwCoN9GPLO3du7dfv36y2/3799+7d2/1Ve+++66RkZFs1ZUrV549eyZLQUIIy7L29vbl5eWEEErp6dOnz507V1RUJHbB6kV7LqV/2ShP9uhT4TEmKQSAtyN6izAtLc3R0VF228nJKS0trfoqT09P2W1TU1MTE5O0tLRGjRrJlly6dOn06dO//PILIcTa2vqnn356/vz5kydPtmzZ0qVLlxrPkpKSsmTJkupLunTp4uXlJdKLajjOPSed7EllZc1uI5WVlZWVGj5lkT4hIz3IL9elP7dRwbPzPF9ZWSmRaMG4dqpTWVnJcRzL4kywiDR+XyGRSBimntaC6F9jnuc57u+DVxzH8TxPKZWVxfN89Y+4RCKRSqWy2w8fPuzfv/+SJUvc3NwEQUhNTZVtuXDhwpEjRz5+/LjGs+Tn59+6dav6khYtWvA8L97raiDisrlJPrW8UJ7nteHlj/MkbfaxU/0Ec11lPzX/D2U/sTaRvcOUotEvIo3/GHMcp/ogdHBwyMzMlN3OyMhwcHCoqqn6qvLy8ry8PFlz8MmTJ6GhoVOnTv3www8JIdXD8oMPPvjyyy+zs7Otra2rP0tAQMCyZcvEfi0NTU45ySytbGmnz770v1xZWamvr6+KopSqqT7p1Zhf+1Dytb+yGw2yn3Ha8CarkCAI+vr6aBGKSkv2FXUT/RMWEhJS1c/l6NGjnTt3JoQUFRWVlZV17tz5+PHjgiDIVrm6ujo7O2dkZISFhUVHR48bN+7lR0tISDAyMrKwsBC7bLVwMZO2tWVeTkGt8lULduktXFwPAG9O9Bbh+PHjW7dubW1tbWBgsHz58tOnTxNCBg4c2KpVq//85z///e9/Bw0aFBQUtGDBglmzZrEs+/777+fn56ekpERHRxNCOnToUFJScvnyZV9f3/T09FWrVn333XdVx1q1XFymoLU9Zao0t2R8zMmfKcIID7QbAOBNiB6ETZs2jY+P37BhA8/zFy5c8Pb2JoSMHz/exsZGIpGcPn36jz/+yMjIWL16dWhoKCFk4sSJBQUFVXd3cXFxdXWVnSa0tLQ8dOhQ27Ztxa5ZXcRl0U99sPcnXzbnvrrEf+ih5W1jAHhDjAaciF6zZk1sbOyaNWtUXYhSUUKs1lfe6a9ja1DL2sLCQhMTE6UXpRqUEP+d0p/bcuFOyotCnucrKioMDGp790FBSkpKcI5QbFq1r3gVfMLU1b18aqHL1JqC2oYh5KsW7IIkTe75BgDiQRCqq4uZNFDrTxBWGdyUvZdPErLV/vAGACgfglBdaeGkE3XQYckEX3Z+EjqPAsBrQxCqq7gstAj/JdqbPf5MSClAoxAAXg+CUC2VSsm9fPqOFYLwf4wk5GNPdslNNAoB4PUgCNVSwgvqY87o4XLKf5vox21MFrLLVF0HAKgVBKFaQk+ZWtkZkL6ubMxtNAoB4DXIG4Tffvvt9evXRS0F5Hcpi7ZFT5nafNmc/f0WXypVdR0AoD7kDcI///zT39+/Q4cOa9asKS0tFbUmqBdmZn+VZmZMGxt2fTIahQAgL3mD8ObNm1u3bjUyMvroo4/s7e2jo6OvXbsmamXwKumlpKiSupkiCGs3qTm7MEkQ0HsUAOQjbxDq6ekNGDDg6NGjt2/fHjt27M6dO1u2bNm6desVK1aUlJSIWiLUIBtrGzH4KiEOjKkuOfQUSQgAcnntzjLNmjX78ccfU1NTJ02adOXKlejoaGdn5xkzZuTm5opRH7wsLpMG2qKXU10m+LK/3MCIawAgl9fenxYXF69atSokJGThwoUuLi5z5swZMmTI4sWL27RpU1hYKEaJUENcFsaUqcegpuzdPHL1BRqFAFC/1wjCq1evjhs3rlGjRmPGjLGxsdm3b19KSsr06dN//fXXpKSkZ8+eHTlyRLxCQUag5Eo2bY0grJMOS8b5sItuoMsMANRP3vkIu3fvfuzYMVtb23HjxkVHR7u4uFRf6+Li4uLigqOjSnArj9oZMFZ6qq6jwYv2Yt22Vj4tZp2M8KMBAOoibxCamppu2rSpX79+urq6tW5w5MgRMzMzxRUGtbuMKwjlY6FHhrqxMbeFOa0xAA8A1EXeQ6NNmzZt27ZtjRR88OBBdHS07LaTkxNmd1SChGwaYI0glMuk5uz/3REKK1VdBwA0bPIG4fr16zMyMmoszMjIWLFihaJLgrokvKABGGtbPq4mTNdG7Mq7OFMIAHV5q174WVlZFhYWiioF6iVQciOH+iMI5TatJbswSahAFALAq9VzjvDYsWPbtm0jhBQUFPz00092dnZVq6RS6dGjRwMCAsQtEKq5l09tDBjz2s/SQi38LZlmZmTLA2G4O668BIDa1ROE9+/flwVheXn5iRMnJJL/bW9jY9OmTZu5c+eKWyBUk/ACcxC+tin+3OQ4fpg7izcOAGpVTxCOHTt27NixhBB7e/udO3d26NBBKVVB7RKzEYSvLcyR0WHJwSe0pzPeOgCoRV1BePHixXXr1r3//vthYWEjR47csGHDhg0bXt7s999/F608+JfEF3RyCxzie21fNmd/vs73dJb3YiEA0Cp17RqePHmyf/9+Pz+/sLCwY8eOZWVlKa0sqNXVF7QlWoSvb6Ar+028cDmLtsElmADwkrqCcMCAAQMGDJDdjo+PV0o98EqPi6gex9gbqLoONSRhySQ/9ufrwtZuuLgeAGqS9zhbZmZmenp61Z87duyYMmWKrB8NKEdCNg2wVnURauvjZuzpdOF+PobhBoCa5A3CsLCwqnOBv//+e//+/ZcuXTpw4MBZs2aJVRr8WyK6jL4FQwmJ9mJ/wTDcAPASuYKwrKzs+vXr7777ruzPhQsXRkZGFhcXx8TELFiwoLS0VMwK4W8Iwrf0mS+3+YGQjk8rAPybXEGYm5tLKXVwcCCE3L9/PyUlJTo6mmXZoUOHFhUVPXr0SNwagRBCSEI2wSijb8NGnwxuysbcwoS9APAvcgWhbFqJzMxMQshff/2lq6sbHBxctbasrEyk4qBKZikp5WljYwThW5ncgo25LRRhGG4AqEauIDQ0NAwICPjmm2927dq1bNmybt26GRkZEULu3btHCHF0dBS3RvhnrG3E4FtqasJ0dmBX38OZQgD4H3k7y/z+++83btyIjIwsLy+fN2+ebOGff/7p7u5ua2srWnnwN1xBqChTWrALbwiViEIA+Ie8QRgYGJiWlvbkyZNHjx41b95ctjAqKurgwYOi1Qb/cwXTECpIGxvGw5RsSkESAsDf5B10qqysLCYm5q+//nr+/Lkg/GsnkpKSIkJh8C9XX9DvWmFwNcWY3pL75Cw/3J3FONwAQOQPwhEjRmzdurVTp04hISEsiz2yUuVXkIxS6mmG3bZihDgwNgZkxyNhgCs+yQAgXxBWVFTs2rXrhx9+mDZtmtgFwcuuZFN/K4ZDDirONH9u+mW+vyvahAAg3znCvLy8ioqKHj16iF0N1CrhBW2FE4QK1dOZkbDk4BOMuAYA8gWhjY2Nm5vbrVu3xK4GanUlmwagy6iife3Pzk7ExfUAIF8QMgyzatWq2bNnnzx5klL8iFa2BHQZFUE/Fzavgpx8js8zgLaTt7PMmDFjHj9+3LVrV319fQODf00FlJOTI0Jh8LeCSpJWTL3QU0bRWIZM9WdnJ/BdemHCXgCtJu8uoHfv3gUFBW/wBJTSlStX7tu3z8bG5osvvvD19a2+9vjx4ytXruR5fuTIkT179hQEYceOHUeOHMnMzGzWrNnEiRNlw9bEx8cvWbKksLBw4MCBQ4YMeYMy1FdCNvW3YiTo3iiCYe7snKvC6XTayR6/MwC0l7xB+PPPP7/ZE/z6669Lly5dtGjR1atXQ0JC7t69a2lpKVt15cqVyMjIpUuX6unpDRs27K+//mrZsuXChQuHDx/euHHjbdu2BQUFJSUl5eTkdOvWbfbs2S4uLmPHjmUYZvDgwW9WjDpKwAlC0XAM+boFOzuRP/ouGoUA2us1vv9nzpxZtGjRjRs3JBLJzZs3CSFLly41MjL66KOPXnUXSumiRYt++eWXnj179uzZMzY2dt26dZ9//rls7dKlS6Ojo0eMGEEISU5OXrx48c6dOy9cuCBb27NnTzs7u0uXLp04caJnz54TJkwghLx48WLhwoVaFYRXsml3RwShWD70YOdcFc6m02A0CgG0lbxH3Hbt2tWlS5eHDx/6+fnl5+fLFurp6c2cObOO7jM5OTkPHjyomqoiODj48uXLVWsvX778qlWEkNzc3Pz8fEdHx/j4+OqbJSYm8rwWdfa7ko1rJ0Skw5Kp/uycq1r0iQKAGuRtEU6aNGn48OGrVq06ffp0XFycbGFISEh0dHRaWpqTk1Ot90pPT2dZ1tzcXPanlZVVenp61dqMjIyqw6RWVlYZGRmUUoZhCCGCIERFRQ0aNMjLy6vGZlKpNCsry97evvoTxcbGduzYsfqSzz//PDw8XM5X12AVS5knxbpOkmSXFLEAACAASURBVOKiote8Y3Exg8kq5DOgEZmTqHsqtaSV5WsMQMrzfEVFhVb9JlO+kpISqVSKoaxEpfH7CkNDw3o/QnIFYWZm5sOHD7dv386ybPW3TNaTJT09/VVBaGRkJAhCRUWFvr4+IaS0tNTY2Lh6fVVzGZaWlhoZGckenFL62WefZWdnb9y48eXNCCHVH0TG19f3q6++qr6kWbNmL2+mdq6m0+aWvIXpa78QSqkGvHylmdpS+Pk2uy/8Nc4UyIKwRg9qUCyWZfX19RGEosK+gsgZhLIPYo2xtgkhaWlppLZYquLg4KCjo/Po0SMvLy9CyKNHj5ydnavWOjs7V81uX33VpEmTEhISjhw5Ymho+PJmFhYWLz+jjY1NSEiIPK9FveBSeuX4yJP94aoQn01b4yg0gPaR66eWtbW1m5vb2rVrCSHVW4RLliyxtbX18PB41R319PR69+79xx9/EEJyc3P/+uuvgQMHEkL++OOPZ8+eDRw4cO3atVKplFK6evVq2arp06fHxsYeOHDA1NRU9iADBw7csmVLcXGx7I6yzbREfDZtbYNds+j0OPK1P/tdAuZmAtBG8h4L+u6774YOHZqZmenp6VlRUbF58+atW7f+9ddfv/32G8dxddxxzpw54eHhZ86cSU1Nfe+99zp16kQIGT9+/J49e0aPHr1r1y4/Pz8dHR0jI6Px48c/efJk7ty5RkZG7u7usrv//vvv/fv337x5s6+vr42NTVFR0bFjx97yNauRhGw6uQWOCynD6GbsvGvSy1m0DX55AGgZRv4h09auXTt16tSq3i6mpqazZs364osv6r1jRUXFjRs3LC0tXVxcZEtyc3NNTEwkEgkh5Pbt2zzP+/r6MgwjCEJVl1QZIyMjXV1dQkhycnJRUVHz5s1fzt01a9bExsauWbNGzheiLkqkxHZDZe6HOjqvH4WFhYUmJiYiFKXJlt4Ujj2ju7vX9cOuCs4RKkFJSQnOEYoN+wryWtcRjhgxYujQoVevXs3KyjI1NQ0ICJBzL6CrqxsQEFB9iYWFRdVtb2/vqtssy1ZfVV1VG1F7JL6gPhbMG6QgvJkxXuxP19EoBNA6cgVhSUnJxo0bT548+ezZM0qpra1tp06d/Pz88HNYVBhTRsn0ODKlBfv9VUHORiEAaIb6g/DGjRu9evVKTU0lhJibmzMMk5ubu3379v/+9787d+6UnfMDMSS8oO1tEYRKNdqLnXddiuk+ALRKPcfdSktLIyMj8/Pzf//99+zs7Nzc3JycnPz8/PXr1+vq6vbr1y8zM1M5hWqhBIwpo3T6HJnmz864gsvkAbRIPUG4bdu2lJSU/fv3jx071srKSrbQ1NR02LBhsbGxxcXFK1euFL9IbVQqJckF1M8CQahsY7zYe/mYpxBAi9QThEePHg0JCQkKCnp5laenZ//+/Y8cOSJOYdruWg71Mmf0cK5K6XRY8l0rduolHkkIoCXqCcLbt2+3b9/+VWvbt29/+/ZtRZcEhKCnjEoNcWOllOx+jOvrAbRCPUGYn5//qusZCCHW1tZ5eXmKLgkIISThBfprqAxDyJzW3NRLghRRCKAF6gnCysrKOgaO4TiuoqJC0SUBIegpo2o9nBhHI7IuGUkIoPnqv3ziwoULRkZGta5KTExUdD1ACCHlPLmbT5ujp4xKzWnNDTrBD2nKGmD6egCNVv9XfNu2bdu2bVNCKVAlKZe6mzLY/6pWO1umrQ2z6KYwzR+j+wBosnr2tXv37sXBT+XDBd0NxLy2bOBu6ShP1h5jKAFornqCsHnz5sqpA6pDl9EGoqkJM9ydnZ3I/9YBF7IAaCwc82mIEl6gp0xD8e073PaHws1cXFUIoLEQhA1OpUBu5dIWlgjCBsFCj3zdgpt2Gd1HATQWgrDBuZVHGxszxjqqrgP+8akPezOXYtA1AE2FIGxw0FOmodHjyNw27JQ4DLoGoJkQhA3OFVxK3/AMaMrqcmRLCg6QAmggeYMwODh43rx5mHRJCTCmTAPEEPJjG256vFCOCZoANI68Qeji4jJjxozGjRt/8MEHsbGxlOIokSh4SpJyqT96yjQ8He0ZPwtm2W00CgE0jbxBuGHDhtTU1P/+978XL17s0qVLs2bN5s2bl5WVJWpxWuhOHm1kyJjpqroOqM3cNuzca3w+RpgA0CyvcY7QwcHh66+/Tk5O3r9/v4+Pj6yBOGzYsDNnzohXn7aJz6ZtbNAcbKB8LZiIxuy8azg8CqBRXruzDMuyPXv2XLFixbhx48rKyjZu3NipU6c2bdqcPXtWjPq0DXrKNHCzAtj/uyukFuHUAIDmeL0gpJQeO3Zs4MCBzs7Oq1atGjNmTHx8/J49eyQSSbdu3R4+fChSldojPou2RhA2YI5GzERfbuIFnCkE0BzyBmFGRsa8efM8PT27d+9+8+bNBQsWpKWlLV++vFWrVhEREadPnzY1NT1//ryotWo8npIbubQlRhlt2Kb4s3fz6b4naBQCaAh5Z/oJCAjIzs6OjIxcuXJl586da6zV0dEJCAjQ0cFoKG/lRi51NmJM8C42bLosiQniRpzi43sRTEoBoAHkDcLp06f379/fzs7uVRscPnxYQSVpryvZtDV6yqiDzg5MB1vy0w32p/aqLgUA3pq8h0bXrVuXlpZWY2FiYqKbm5uiS9Je6CmjRua3ZdakMLfycIAUQO3JG4SPHz8uKyursbCkpOTRo0cKrkiLIQjViJ0B+aY5/eQsBiAFUHtvNdZoUlJSHQdL4bVUCuQmesqolSgPQaBkxR30IAVQb/WcI/zjjz/mzJlDCMnOzh4wYIC+vn7VqvLy8mfPno0aNUrcArXGjVzqasIYyXvSFlSPZcjKjlzn/dJezoyTEX7BAKireva7jRs3Dg0NJYSsX7++TZs21dt/+vr6fn5+H374obgFag1cQaiOvMyZsd7sJ2f5feH4CQOgrur59oaGhsqCUBCEyZMnN2vWTClVaSOcIFRT01tyAX9Jtz0UBrhiUjMAtVTXV3f79u1NmzZdtmwZIeTatWvvvvtu09ooq1QNhyBUU7osWdmRm3iBzylXdSkA8EbqahE6Ojp2797dxcWFENK5c+eCggIlFaV9KgRyOw+zL6mrdrbMwKbsl3H86k6cqmsBgNdWVxC2b9++ffu/Lxj++eeflVKPlrqeQz3MGAOcZlJbP7TmWuyUHn5Kw53wawZAzdS1692+ffuUKVPqfYgHDx4orh4thZ4y6s5QQlYEcx+f4ZP6SowxSB6AWqn/0KjSStFmV7JpAIJQzXVtxHS2Z769wv/SDgdIAdSJvIdGQVSXsugn3uhzqPYWteea75D2d6VBdvhZA6A25O012rZt21q7jKLX6NsrkZKUAtocPWXUn7kuWdiOjT7Ll2MSewD1oaReo5WVlSzLclwth4wopRUVFXp6em/84OouPpu2sGR00SDUCANc2c0pdM5V/rtWOEAKoB5E7zVaUVERFRW1a9cuhmGioqLmz5/PMP9r+syfP3/OnDmCIAQFBW3atMnc3PzQoUMxMTGJiYnt27ffsmULISQ9Pd3Hx6fqLtOnT588efKbFdMwXc6ibW3RHNQcvwVx/jsr32/C4rwvgFp4jWbIkydPpkyZEhYW5uvr261bt4kTJ969e7fee/322293795NT09PSUnZvXv3zp07q1bFx8fPnTv38uXLL1680NfX/89//kMIYVm2b9++kZGRRUVFss0EQSgpKUn5x2efffaar7Ghu5xF22CPqUHsDci8NtzHZ/hKDMcNoA7kDcJz5875+fktWrSoqKjIy8uroqLi//7v//z9/Xfv3l33HdetWzdhwgRDQ0Nra+uoqKi1a9dWX9W/f393d3eJRDJ58uT169dTSsPCwkaMGNG4cePqD8IwjMU/qg/8rRkuZdE2mI9Xs4z0ZO0NyIIkJCGAGpArCCmlo0aNatq0aXJy8vnz53fs2HHmzJlHjx517Njx448/fnmewuqSk5O9vb1lt318fFJSUl61Kjc398WLF7U+SEVFhb29vaOj44gRIzIzM1/eoLy8POvf6q6q4cgqI3kV1MMMQahplgVzC5P4O5i5F6DBk2ssk4yMjPv37586dap6Q83W1nbZsmXu7u63bt0KCAio9Y6VlZVFRUXGxsayP42NjXNzc6vW5uXlVV9FCMnNzbW2tq7xIGZmZqdOnQoICHj+/Pn48eOHDx9++PDhGtscOXKkKlNl5s2bN2jQIHlenWqdfsa+Y8EV/3McWFGKi4urn4sFheN5vqKigudf2T3UipAZzbkhJ+jx0Ar0hHozJSUlUqmUZfH2iUjj9xWGhob1foTkCkJjY2OWZc3NzWssly0xNTV91R11dHTMzc3z8/Nlf+bl5dna2lattba2rr6KEFJ9bRUjI6Pg4GBCiJub29KlSz08PKqHq0xERMSaNWvkeS0NTVIR38GBGBsruNMspbTGWwSKJQtCAwODOraZ4E/OveB/uG0wPxA9SN8Ey7L6+voIQlFhX0HkPDRqbGzct2/fxYsX11i+ePHiwMBAd3f3Ou7r6+ubkJAgu52YmFi93ebj45OYmFi1ytHR0czMrO5KKisrGYbRpC/GJfSU0WjLgrjtD+m+VBwgBWi46moR3r1799SpU7LbgYGB33//fUJCwvvvv+/g4JCdnX3w4MH4+Pjp06fX/QRjxoyZPXt2SEhIQUHBihUrNm/eTAgZN25cjx49Pv7441atWh06dMjHx2fWrFljxowhhKSnp9+4ceP+/fvZ2dnHjh1r1KhRRkZGWVmZn59fenr6l19+GRERYWhoqKCXr2KUkMtZdFVHzcl1qMFCj6wP4QadkCZE6tjX1XoEAJWpKwjPnz8fHR1dfcnVq1evXr1afcnMmTO//fbbOh5k+PDhT58+7dOnj66u7vfff9+lSxdCCKWUEOLh4bFu3boZM2bk5eX16dNn6tSphJCbN2/OmzePEGJqajpv3rzw8PAWLVr8+OOPqampFhYWoaGh33zzzZu+2AbnYSE14BgHDYl1qF1HeyaqGTvylPRgDwna/gANECPLpFoVFhZmZWXV+xAqH2VtzZo1sbGx6niOcMsDYesDuiNU8SeQCgsLTUxMFP6wUEWec4RVpALptE86qCk70Q+t/9dQUlKCc4Riw76C1N0iNDExqf4GSaXS+Pj4lJSU4uLi6pvJDmnCG7iURQMxpowWkLBkQxeu3W5pl0ZMCwwqC9DAyDsV7L1793r16pWcnPzyKgThG7uURb9rhV+7WqGpCTOvLffBST7+fYk++pACNCTy7oXHjBkjkUguXLgwZMiQSZMm3b9/f/HixQ4ODvWOLAOvwlNy7QVthS6jWmOUJ+trwXwTj5kpABoWuYKQ5/kLFy7MmzevXbt2enp6enp67u7uEyZMmD9//qefflrHNcVQhxu51NmIMcVs5trk9yBu6wN6LA1XUwA0IHIFYVZWVkVFhZeXFyHEyMio6ir4iIiIp0+f3rlzR8QCNdelTEw6oXWs9MiaztzI03xmqapLAYB/yBWEVlZWHMdlZ2cTQpycnKqugn/+/LmIpWm6i5k0EGNta59ujZhRnszgE1IezUKAhkGuINTR0WnVqtWZM2cIIX379o2Li4uOjl6+fPngwYMbNWrk4eEhcpGaCV1GtdZ/AziGIT9fx9wUAA2CvJ1l5syZ4+zsTAjx9PRcsGDBtm3bPvnkk5ycnE2bNunq6opZoWYqrCQPC6mfBYJQG7EMWdeZW3yDP5eBViGA6sl7+URoaGjV7c8///zzzz/Pzc21sLAQpyrNdyWb+lsxOrh0Qls5GjErO0mGnuTj35dYa9oMmwBq5jX2xMnJyePHjw8ODnZ1dW3fvv2UKVOuX78uXmWa7VIWbYsThNqtlzPzgRsz+IRUikOkAColbxCeOHHC399/5cqVBgYGQUFB5ubmW7Zsad269datW0WtT1PFZdJ2OEGo9b5vzelzZMolXIAEoEpyHRqllEZFRTVv3nz37t12dnayhXl5ecOGDYuOjo6IiJBzxEWoEpdF5wfiwKi2Yxmyqauk3W5p83vCKE98HgBUQ67vXkZGxsOHDxcsWFCVgoQQc3PzX3/9NS8v7/bt26KVp5mel5AKnrqaoEUIxFSH7Azlpl3mL2eh4wyAasgVhCYmJhzHGRkZ1VguW/LyzPVQt7gsoQ1OEMI/vMyZmCBuwHFcZQ+gGnIFoZGR0aBBg+bPny8I/zqt/9NPP3Xs2FHl0zCpnbhMGmiL42DwP5Eu7EhPJvKYtAIdZwCUTt4Z6lu3bv3DDz/4+PhERkba29vLZqi/e/duvTPUw8suZdEpLRCE8C//CeBu5fKfnOX/6ITJKQCUqq6JeVevXv3RRx/V+xB1PIJyqNfEvAIllusrUwbpWOmJ+CyYbFNsrzUxr5yKKkn7PdJPfdhPvPE7iRBMzKsU2FeQuluEQ4cOff/995VWipa4k09tDRhRUxDUlLEO2dWdC9or9TZnOjvgLDKAktQVhLq6ujWGT3v8+PHVq1efPn3q4ODg5+fn6ekpcnkaKA5jbcOruZkyW7tJBp+QxvaSeJrhcwKgDPIOsVZRUTFmzJh169ZVPxDau3fv9evXm5qailObZsKYMlC3TvbMnNZcz8P8hd4SG4y+BiA+eQ++f/XVVxs2bBg/fvy5c+cePHgQFxc3Y8aMo0ePRkVFiVqf5onLRBBCPUZ5sv1cmAHH0YkUQBnkahFWVlauWrVq9uzZ06ZNky1xdXVt27ath4fHqFGjsrKybGxsxCxScxRUkvsF9B1rBCHUY24bbsBxPvosvxqdSAFEJu8M9cXFxe+9916N5REREYIgPHr0SPF1aaiLmbSVNaOLTnBQH5Yh60O4pBz6w1W0CgHEJdcu2czMjOO4O3fu1FguG1zN2tpa8XVpqLPpQrAdmoMgF0MJ2RsmWX5H+DMFWQggInlHlgkLC/vss88OHjxY1Vnm3LlzI0eO9Pf3d3V1FbNCjXIugwbZoT0I8nIwJLu7cxMu8BczMRIpgFjk3SnHxMSYmpr27NnT3Nzcy8vLysoqODi4sLBw3bp1otanSaQCuZyF2Zfg9bS0YlZ15Pof51OLkIUAopD38okmTZpcvXp106ZNZ8+ezczMbNWqVbt27YYPH44Rt+V3NYc2MWYscCk9vKbeTdiHhaTHIf5MhARDMQAonFxBmJubO3DgwJkzZ44ePXr06NFi16SpzqbTjvZoDsKbmOjHZpTSdw9JT/SUGOuouhoAzSLXoVGe548dO6ajg+/fWzmXQYPQUwbe1Jw2XEsrps9RaTkmtAdQKLmC0Nra2sfH59KlS2JXo9nOZ9AOCEJ4UwwhMUGcuS4z4hQv4HQhgOLI21lm1apVv/zyy6pVqzIzM0UtSFM9KKSEEMxKD2+DY8jGLlx2Gf30PI8oBFAUeYPw/ffff/ToUVRUlJ2dHfNvotanMc6m02CcIIS3ps+RPWGSG7l0+mUcIQVQDHl7jc6YMaO4uFjUUjQbThCCohhKyL4wSch+qbmu8LU/LksFeFtyBeGjR488PT1tbW19fX3RZebNnM+gY7ywzwLFMNMlB8IlnfZJzfVIND5XAG+nnq9QRUVFv379XF1dw8PD33nnHQ8Pjxs3biinMk2SW05Si6i/JVqEoDAOhuTwu9z3icJmDMAG8HbqCcIlS5bs3Lmzd+/ev/766+TJk7OyskaNGqWcyjTJhUzaxoaR4Ic7KFRTE+ZgD+6Li/yBJ+g6A/Dm6jk0evTo0W7duu3evVv2Z7NmzUaPHp2bm2thYSF+bZrjTLrQyQExCIrnZ8HsDpNEHJFu6ybphN5YAG+knr3z48ePu3btWvWn7DbmXXpdZzCmDIimrQ2zqYtkwHFpHAbmBngj9QRhWVmZgYFB1Z+GhoaEkNLSUnGL0iylUnIthwZiVnoQTbdGzNrOkveOSGOfIwsBXlv9vUbv3bt37Ngx2e2cnBxCSHx8fElJSdUGoaGhIhWnGeKyaHMLxlDeC1UA3kQPJ2Z7N8mgE9I/u0q6OOBXF8BrqH/3vGzZsmXLllVfMnHixOp/Vs1Q+CoPHjw4d+6ck5NTSEhIjQvwc3Nzjxw5oqenFxYWJmtuEkJevHhx48YNFxeXJk2ayJaUl5cfPXq0sLAwNDTUxsam3pobFBwXBeXo7MBs6iIZfEK6paskBFkIILd6gvC3336r3vh7A3v37h01alSfPn0uX77s7e29ZcuWqlUPHz7s0KFDx44dCwoKpk+ffv78eXNz8z59+hw+fFgikcycOXPKlCmEkNLS0o4dOxoZGTk7O0+YMOHkyZN+fn5vU5KSnUkXJvhyqq4CtELXRszWrpKBJ6Sbuki6NUIWAsilniDs1avXWz7BtGnTFi1aNGzYsIKCAnd397i4uMDAQNmqn3/+OSIiYsWKFZTSHj16rFy5cvLkyYsWLXJycurfv3/VI2zatIll2RMnTnAcN23atO+//37z5s1vWZXSSAUSl0k7dMUuCZSkswOzM1QSeVS6qhMX0Rh9lQHqJ+735OHDh3fu3Onbty8hxNTUNDw8fO/evVVr9+3bN2DAAEIIwzD9+vWTrXJ1da0xeM3evXsjIyM5jiOE9O/ff9++faLWrFgJL6irCWOJyVRBiYLsmAM9JGPO8Hse41p7gPqJ24UjLS3NwsKi6uSfo6NjWlqa7LYgCM+fP3d0dJT96eTkVLXq5QepvllxcfHLFzKmpKT8+uuv1ZeEhIR4e3sr8LW8mdPPaLAd4Xllj4/M87zyn1Sr8P9QdSG1e8eC7OzGRh7nBUGIaKyuByRk73C9vRDgbTTkj7FCsCxb7+QQ4gYhz/Ms+79Gp0QikUqlstuCIAiCIGvnEUI4jqta9fKDVN+MEPLylvn5+UlJSdWX+Pn5VVZWKuJFvJXT6dxgV6GyUtnf5MrKyobw8jUYz/OVlZUSScPtDfyOOdndlelznHlRSoc2Vcssqays5Diu+j4EFE7j9xW6uroqDkJ7e/vc3FypVCrbX2RkZDg4OPz9xBKJjY1NZmZms2bNCCGZmZlVq2pwcHComgQxIyNDV1fX2tq6xjYBAQHLly8X62W8KUrIhazKZR119fWV/ZO8srJSX19fyU+qVWQ/8hr4mxzoQE6+R8MP8oRjP26mfnEiCIK+vj6CUFTYVxCxzxG6ubnZ2dmdPHmSEMLz/IkTJ0JCQgghubm5giB07tz5yJEjsi2PHDkiW/Wyzp07Hz16VHb76NGjnTp1UpdJEG/lUnNdppGhelQLGsnbnDnek5udKCy+gfOFALUTt0UokUimTp0aFRU1adKk2NhYS0vLHj16FBYWWlpaJiUlffXVV926ddPV1S0oKDh48GBiYiIhZP369WfPnr127VpGRkZKSsrIkSM/+uijRYsWffrpp40bN/7xxx+3b98uas0KdCadYvhHUDkPM+bUe1zPQ3xyAf2lHYfB3wFqEP0Mx6effurm5hYbGxsSEvLxxx/LjiYtW7bM0dHRwsLi7Nmz27ZtMzMzi4+Pd3Z2JoQ0adKktLS0VatWsrtbWVlZW1vHx8evXbu2oKDg2LFjVasavtPpNMwRQQiq18SYudhHMvQk3/OwdEtXiQW6MQNUw2hAj6w1a9bExsauWbNG1YXU5PynNLYX52aqgiwsLCw0MTFR/vNqD57nKyoqqo/E2/DxlHwVxx94QveGcR5mavATraSkBOcIxYZ9BRH7HKE2Sy6gDCEqSUGAWnEMWdiO+7I522mf9FyG2v8CBlAUBKFYTj6jXTDGFTQ8o73Y9SGSvsekmzC1PQAhBEEonpPPKQY+hoYp1JE5+q5k2mVh3jVkIQCCUDSxzwUEITRYLSyZ8xHclgfC6DN8JdIQtBuCUBR38qgOy7iaIAih4XI0Ys5ESLLKyLuHpHkVqq4GQHUQhKI4+ZxiEhxo+IwkZEco52/FdNgjTSlA9xnQUghCUcTiBCGoCY4hCwK5iX5s0F7p4afIQtBGCELFo4ScwglCUCvRXuyOUMlHp/mpl3mEIWgbBKHi3cqlJjpMY2MEIaiTIDvmUh8u9jkddJwv0uTZCABqQhAqHi6cADXlaMSc6iUx0yXt90jv56NlCNoCQah4J5/RrugpA+pJjyP/15H7qgUbvE/61yNcVwFaAUGoYAIlp9OFzmgRgjr70IPd013y+UVhVgIvoGUImg5BqGBJudRKD3MQgtoLtGUu9ZHEPqfvHZHmlKu6GgAxIQgVDEOMgsawMyDH3pV4mzNtdkkTX6BhCBoLQahgR9KE7piDEDSFhCULArkf27LhB6Uxt3GUFDQTglCRynhyLp12a4R3FTTKAFf2TITkj7tCxGFpRqmqqwFQNOyyFelMOm1hxZjpqroOAEVrZsac7y3xt2Le+atyXypahqBREISKdPipEO6ItxQ0kw5L5rTmtnaVTLjAjz3Hl0hVXRCAgmCvrUiHntIezjhBCJos2J652lciUNJql/RKNpqGoAkQhAqTVkwzS2mAFYIQNJypDlkezM0KYHsels65Kkhx2T2oOQShwhxOo6GOLIscBO0wqCl75X3JmXShw17prTw0DUGNIQgV5vBTGu6EGAQt4mTEHOwhiWrGhuyTzk8SMG8FqCkEoWLwlBxPE8LQUwa0DEPIGC/2Uh/JwSdC8F7pbTQNQQ1hx60Yl7OokxHjYKjqOgBUwcWEOdZTMtabDdkvnXqZr8BZQ1ArCELFwHFR0HIMIR96sJf7SK69oO12o0MpqBMEoWIceir0cMabCdqusTFzsIdkWku212Fp9Fm+EBP8gjrAvlsBXpST23k0yA4tQgBCCBngyt7op1PGkxY7pRiGBho+BKEC7E8VQh1ZXbyXAP+w1idrO3MrO3KT4/jIo/yTYsQhNFzYeSvArse0TxM0BwFq6taIudZXEmDNBPwl/fm6UIlONNAgIQjfVqmUnHgm9MIJQoDa6HHk23fYi70lp9OFFjulR9LQNIQGB7vvt3Xsw3XrDAAAH6tJREFUmdDKmrHUU3UdAA2YmymzN0wyP5D79Bzf9xj/sBBxCA0IgvBt7XpE+zTB2whQv17OzI1+ktbWTNvd0m+v8MWYvwIaBuzB3wpPyb4nAk4QAshJjyPTW7JXIyUPC4n3NunGZAHT3oPKIQjfyoUM6mjINDFGEAK8BkcjZkMIt7krt/SW0HoXThyCiiEI38quxwKOiwK8mQ52zIXekm9ashPO86EHMBgNqAx24m9l92P6vguagwBviCGknyt7o59kYFO2z1H+g5PoRwMqgCB8czdzqZQSf0sEIcBbkbBkjBd7d4DEx5xpu1v6+UX+eYmqawJtgiB8c7iOHkCBjCRkxjvsrf46LCF+OyonXOCfIQ5BKRCEb27HQyESJwgBFMpGnyxsx93qr6PLkrYHJGPPCffycbAUxIX9+BtKKaDPSmiwPVqEAIpnZ0DmB3JXI3g7A9JxnzTyKH82HXEIYpGI/QQVFRULFiw4f/68i4vLtGnTGjVqVH3t5s2bt27damxsPH78+DZt2lTf3tXVddq0aQ4ODnl5eV9//XXVXXr37t2rVy+xy67Xjkc00oXlkIMAorHWo7MC2Kkt2bX3hZGneWcjMjOA6+KAbx0omOgtwi+//PLAgQMTJkwQBKF79+48z1et2rp16xdffDFixIi2bduGhoY+fPiQEDJp0iTZ9lKptHv37oIglJSUrF69utU/HBwcxK5ZHtsfCv1c0J4GEJ2hhIz1Zu/0l4zyZD85y3faJz2K6w5BsaiY8vLyDA0Nb9++TSkVBKFJkyb79++vWtu+ffsVK1bIbg8bNmzq1Km5ubmGhoZ3796Vbd+4ceMDBw6kpaXp6+vX8SyrV68eMWKEiC/jJY8KBZv1FZW8Mp/z9RQUFKi6BA0nlUpLSkpUXYWGKy4u5vl/fc2kAl1/n/fZVtl2V+Vfj3hBVZVpEOwrKKXitmlu3bplaGjo5eVFCGEYJigo6PLly1UBfOXKlaCgINmfwcHBly5dunXrlrGxsaenZ43tpVLp+PHjv/jii3379olasJx2PKR9mrASNAgBlItjyDB3Nqmf5Gt/9vtEocUO6ep7Qhlf/x0B6iDuOcL09HRLS8uqP62trdPT02W3c3JyKioqqtZaWVmlp6fXur2ent7YsWMDAgLS09NHjx4dHR09a9asGk8UGxvbsWPH6ks+//zz8PBwMV4UIWRLss5UP76oqOHOrlZcXMwwOJUiIp7nKyoqqh/qB4UrKSmRSqUsW8tPzjBrEhZKTqSzv99lp15iR7nxUe6CvQEOmb42jd9XGBoa1voRqk7cIDQyMiovL6/6s6SkxNzcvGoVIaRqbWlpqbGx8cvbW1paWllZLVmyRLYkMDCwZ8+eM2bMkEj+Vbmvr+9XX31VfUmzZs2MjY1FeE3kaTFNKZK+56an04BbhJRSkV4+yMiC0MDAQNWFaDKWZfX19evYi/V2J73dyb18uvQmF3hIiGjMTmrOtsAYF68D+woidhA6OztnZGSUlpbK9hePHj2KiIiQrdLX17e2tn706FGTJk1kq5ydnZ2dnZ8/f15WVqavry9b6O/vX/0Bvby8ysrK8vPzraysqi+3sbEJCQkR9bVU2fmIRjRhG3IKAmgVTzNmaQfuu1bc8jvCu4d4PwsyqTkX5qTRzRxQKHF3597e3u7u7ps2bSKEJCcnnz9/vm/fvsXFxStWrCguLh44cOAff/xBCCktLd28efOAAQN8fHzc3Nxk29+/f//ChQuRkZHPnz+vrKwkhFBKly1b5u7uXiMFlQz9RQEaIAs9MtWffThY8oE7O/Uy77dd+n93hFJMeQjyELs3TmxsrK2tbXBwsLW19bx58yilT548IYSkpqampaV5e3sHBAS4urr27du3srKSUnry5Mmq7X/66SdK6ZIlSywtLdu1a+fi4uLu7n7x4sUaT6HMXqPPiqnFuooyqXKe7c2hJ5jY0GtUCV7uNSq/E8+EiMNSm/UVX1yQ3spF99JXwr6CUspQKvrp5eLi4tu3bzs5Odnb28uiNzc319zcnGVZnueTkpKMjIw8PDxetT0hJCMjIzU11cLCwsXFpcbZQULImjVrYmNj16xZI/YLIYQsuiFcy6GrO3FKeK63UVhYaGJiouoqNBnOESpBSUlJ3ecI6/WwkK66K6y+R5uaktHN2AGurIHog4ioGewrCCHKCEKxKTMIW++S/tSW69qooZ99wIdbbAhCJXj7IJSRCmT/E2HFHSEukw7zYKO9WG/zhv4VVhrsKwjGGn0tt/JoRikJwQhPAGpFwpI+Tdj94ZIrkRJjCQk9wAfvla65JxTjDCIQQhCEr2XDfeEDN4ZFDgKopybGzPetuceDJV/7s7se08Z/Vkaf5eOz1f6oGLwlHC+Xl0DJxhS6P7yhnx0EgLpJWBLRmI1oTNJLudX3hEHHeXM9MsaLHdyUNdNVdXGgCmgRyut0OrXUI34WaA8CaAh7AzLNn70/UPJjG+54GnXZXDnoBL8vlUob7phRIAq0COW1/r4w3B2/GwA0DcuQ7o5Md0cut5zb9lD48RofdYYOasoOc2fb2OCHr1bAnl0upVKy67EwxA1vF4DGstAjY7zYsxGS870l1vrMsFjea5v0+0ThYSFOImo47NnlsjdVaGPDOBiqug4AEF9TE+bbd9i7AyRrO3OZZbT9HmnwXunvt4T0UlVXBuLAoVG5rLkvDMNxUQAtE2jLBNpyCwO5I2n0zxRhxpVKPwumvyvb14VxMsJRU82BIKxfQja9nkN2hiIIAbSRhCU9nZmezlw5zx1NozseCd8l8B5mTD9Xtq8L09QEiaj2EIT1+y5RmNKC1cd1EwDaTY8j7zVm3mvMVQpc7HO646HQYQ/fyJDp68L2dWV8MFqN2kIQ1uPqC3o5i/7ZBTEIAH/TYf/uaPpbEHcug+58JPQ4KOhxpLsjE+rIdG3EmuN6RLWCIKzH7EThqxYYqBcAasExpJM908me+6UdScqhR9Po/90RRp3ifS2YcCc23IlpY8NwaCg2eNjB1+V6Dr2QKWwI0VF1IQDQoDGEtLBkWlgyXzZny3lyLoMefiqMPSc8KaJhTux7jZkeTqylnqqrhFdAENZldqIwuTmH5iAAyE+PI10bMV0bcfMIeVZC96fSrQ/o2LOVfpZMmCMb5sS0RTOxgcE+/pWScujZdGFtZzQHAeANNTJkRnsxo71IOc+dzaBH04RPzwkPC2lHe7ZrI6ZrI6a5JYNMVDkE4StNucRPb8kZ4h0CgLemx5FujZhujbgf25CsMnLymXDyOY25LeRX0FBHtrsjE+bIYsgOVcFuvnYHntDHReQTb1w7CAAKZqNPBjZlBzYlhJDUInokjR54Qr+8WGlrwHRpxHR1YDo5sDb6qq5SmyAIa1EpkC/j+IWBnA5yEADE1NiYiWrGRDUjAuWu5dATz+jqe8LHZ3g7A6a9LRNkxwTbM964QlFkCMJaxNwWXE3Iu8748AGAkrAMeceKeceK+bI5K1ByK4+ey6BnM+iP14QiKQ22Y4Ptmfa2TIA1o4sf6IqGIKzpRTmZc5WP7YV3BgBUg2WInwXjZ8FEexFCSFoxPZ1Oz2XQdfeF5ALa0ooJtmM62rNBdgxmElYI7O5rmnWFH9SUxbEIAGggHI2YIW7MEDdCCCmqJJey6Jl0ujCJH/L/7d17WBNX2gDwd0IAhSRAICFcRFQuKiAXaR9BwHJR1g+qKKvb27befdzWlaLddr2t27K6T+utCq6tN0StQNVvV11Q0CiUSxUQCygXtSLQcDUBAgQmyZzvj/HLsmpbtcUY8v4e/5g5Mzl5M85zXmbmzDlSMoZPBdtTU8RUkJhyt8JW6xlhIvwvRa3kZD1TFY+vTCCEXkQ8U/YlRQqAo2bgWgf5to1kNZINZUyvmrwspl4WUS+LOFPEFL6//+QwEf5HrwYW5GtTgk3wBEIIvfhMOexEUdQqAABoUcHVNuZqO9leqS1pJ46W1BQxNUVMvWRH+Qgp7Pr3EzAR/seHV7VBYmqOK54vCCHDIxkJs0ZzZo0GANASuKEgxW3kahtJucl83028hezFIvWyiHK3wrf4/wsmwgcu/EDONJDv5uIBQQgZPBPqwdinbHebXg2Ud5Cr7eRsA9lYxigGyEsiiv03zpzy5oOR50Vs9wEAOmlY/I12f6gJTp6CEBp+LLkQIqFCJA/yXXs/lLSTknZyoJb5rsNMoVZPtKZ8hJSfLeUnpHxtKb6RdZPARAgMgYV52lku1HQnI/+rCCFkFEQj4H9GUf8zigLgKJVKxpx/Q0EqFaS8gxy9zVTJiYMF5W9L+dtRfkLKWwijLId524iJENaVahU0yYjEQ4EQMkZWZhBsTwXbP8h2WgK1XeT6fXL9PtlexVTJSb8WvIUU+2rjRBvK24YaZiPAGXvrf/gWc+Iu+XY2FwdrQAghADChYKI1NdGaemPcg5L7A1ApJzcVpFJBMr9nqhSEywEfG8pHSE2wptytKDcBOFsacAcco06EBS3kw6vayzFcW3xfAiGEfoStObziQL3i8J9MJ+sjVQqokJOr7eTYbeZWN+mkwV1AeVhR7gLwsKLGCSg3AWUo82kYbyIsv0/mSzVHXuGOx0FkEELoaThaUI4WMGNQv4oeNdzqJnVdpK4LpDLyZQ1zR0l61TBOQI0TUOP4MFZAjeVTY/gwmkeZm+gx9scw0kR4UUbeuKTZO9UEO8gghNAvxzN9MGj44EKlGu50kzvd5I4Srt8n/1vP3FVCYy8RjaBc+TCGR43hw2g+5WJJufD0mSCNMRFmfM+sKtZ+HckNk2AWRAihocI3BT9byu+/s6OWwA+9pL4H6pXkrhIKWsi9HqahBxp7ia055coHVx7lwgNnS8qFB6MsKWdLym6I++YYXSIsv0/WXGFyZ3J9hJgFEULoeTOhwIVHufDgoUsRhkBzH7mrhPoe0tgLNzvJ+SZo6GEae4lKC648ytECHC0oR0uId+W8JPo1G3CjS4T+tlTFXK4N9o5BCKEXCYcCJ0vKyRJCHhnoplcD93qIrBdkfaS579f/aqNLhACAWRAhhAyIJZd9owNgaAaDw7fnEEIIGTVMhAghhIwaJkKEEEJGDRPhMNTf33/u3Dl9RzHM1dfXX7t2Td9RDHMlJSWNjY36jmKYy8rKomla31HoGSbCYej+/furVq3SdxTD3OXLl7/44gt9RzHMpaSkFBYW6juKYW7FihVKpVLfUegZJkKEEEJGDRMhQgghozYc3iO0s7OrqKiIjY3VdyAvCpqmNRoNHpAh1dra2tXVhQd5SNXU1Ny7d++rr77SdyDDGYfDefPNN7nc4ZALHmv37t1jxoz56X0oQsjziWZISaXSvr4hGG8AIYSQIQsLCxMIBD+9zzBJhAghhNCzwWeECCGEjBomQoQQQkZt2D4gNQY9PT379++vra21tbVdsGCBm5sbANTX1x84cKCjo8Pd3X3x4sVWVlYAUFtbu2/fvv7+/nnz5k2bNk3fgRuS27dvHzp0SC6Xe3p6LlmyhMfjAQBN06mpqWVlZQKB4PXXXw8ICAAAqVR68uRJCwuL5cuXs/8X6AmdO3fu3LlzDMOEh4fHxcVR1IOBlTUaza5duzw8PNhOSQzDpKWlFRcXjxo1auXKley5jZ5QVlbW+fPnASAiImLWrFmEkPz8/AsXLrDn9uLFi9lzW9dWzJ8/PywsTN9RPyd4RWjAYmJicnNzIyIiCCGBgYENDQ3Nzc2TJ0+maToiIiI/P3/GjBkA8MMPPwQFBVlYWHh5ec2ZMyc3N1ffgRuMe/fuBQYGUhQVERGRm5sbExMDAGq1evr06enp6f7+/mKxuKqqCgCysrLmzZvn4+NjZmY2ZcqUlpYWfcduMFJTUxcuXOjl5RUQEPD+++/v3LlTt+mzzz5LSkrKyMhgVzdu3Lh9+/bg4ODKysoZM2Zg/4Ynd+DAgSVLlvj4+Pj5+a1cuTI5Obm9vX358uVcLjcgIODs2bOhoaE0TTc1NQUFBVlaWnp5ec2ePfvChQv6Dvx5IcgwdXZ2AkBjYyO7Onny5LS0tMzMzAkTJrAl9+/fB4Dm5uYNGzbMmzePLdyxY0dUVJR+IjZAhw8fDggIYJfZsb4UCkVKSoqfn59Goxm857Rp03bv3s0uz5kz5+OPP37esRqsuXPn6g5XcnJyeHg4u1xTU+Pr65uYmPjWW28RQnp6eqysrK5du0YIUavVjo6OUqlUXzEbnNmzZ//tb39jl3fu3BkVFaXVarVaLVuiUqksLS2Li4vXrVv3u9/9ji3ctm0b+9eGMcArQkPF5/Pd3NwKCgoAoKmpqaGhwdfX18fHp7W1tb6+HgC++eYbFxcXW1vbwsLCyMhI9lMRERE4ZtWT8/X1bWxsbGpqAoCCggI3NzeBQJCTk/Paa69lZmb+9a9/ZS+vCSHFxcXh4eHsp/AgPxV/f/8rV66o1WqtVltUVMTeZ2YYZunSpcnJySNGjGB3q6qq4nA4/v7+AMDlcsPCwtiTHz0Jf3//q1evqtVqjUZTXFwcEBDA4XA4nAftv0qlGhgYsLGxMdq2Ap8RGioOh3Py5Mno6Og1a9bI5fIdO3ZMmjQJALZt2+bl5WVjY0PTdHZ2tqmpaUtLi52dHfspkUikUqk6Ozutra31Gr5h8PX1/eSTTzw8PIRCoVarzcnJ4XA49fX1169fj4mJcXd3X7Zs2dKlS5ctW0bT9OCDLJPJ9Bu5Afnzn/88d+5csVjM4XAmTZp08OBBANi5c+ekSZNCQkKys7PZ3QafxgAgEomam5v1E7EBWr9+/ezZs8ViMUVRAQEBaWlpg7e+9957c+bM8fT0fKit6O3t7e7u/tmX8IYBvCI0VH19fb/97W8TExNLS0vPnDmzadOmoqKisrKyDz/88NSpU6WlpevXr4+Pj1cqlebm5rrR5dkFc3NzvcZuMIqKijZt2nTmzJnS0tLExMT4+Pi+vj5TU1NfX9+UlJSEhISDBw/+/e9/Z4+nWq1mP0XT9MiRI/UauCHZtGlTV1fXt99+W1JSYmFhsXr16rt37+7duzcpKWnwbubm5rojDHiQn9L69ev7+vquXLlSUlJiamr6wQcf6DZ99NFHNTU1+/btA4CH2gqKoszMzPQT8fOFV4SGqri4uKuriz2hJRJJbGxsRkbGyJEjp0+fHh0dDQB//OMft2/fnpeX5+TkxN7cA4DGxkZbW1tsQZ5QRkZGbGwse7Pogw8+2Lp1a3FxsbOz87hx49gdPDw8lEqlVqsVCARNTU3Ozs4A0NjY6OjoqM+4DUpaWlpKSoqnpycArFu37je/+Y2np2dra+vkyZMBQC6XazQa9hFsW1sbTdNs09zU1KS7F41+1pEjR/bt2+fh4QEAa9eunTVr1ueffw4AmzZtysrKkkqlbBfch9oKOzs73a3p4Q2vCA2Vvb19Z2cnewuOEFJdXS2RSCQSSU1NDcMwANDe3t7e3i6RSOLi4jIzMzUaDQAcP348Li5Oz6EbDolEUl1dTQgBAJlM1tnZKZFI4uPjCwoK2IOcl5fn4uJibW0dFxd3/PhxAFCr1SdOnJgzZ46eQzcc9vb21dXV7PLNmzcdHBzeeeed8vLy3Nzc3Nzc119/PTIy8siRIz4+Po6Ojv/85z8BoLm5OS8vD8/kJ/foQQaArVu3pqen5+Tk6G6HGm1bgUOsGbAVK1acPn06Kirqxo0bNE3n5eWZmJiwne58fHykUmlUVNShQ4dUKtX06dMHBgbEYnFFRUV+fv7PDkGLWHK5fNq0aSNGjJg4ceKFCxdmz569Z88emqZjY2Plcrm7u7tUKk1NTZ05c+atW7deeeWVgIAAmUzG5/PPnz+P95+fEJvtQkJCuFyuVCo9ePDg4PZ33bp1DQ0NR44cAYDTp08vWrQoIiLiypUr8+bN27p1q/6iNjDnz59/4403wsLCOBzOpUuXUlNTvb29x40bN3bsWBsbG3afLVu2hISEREVFqdVqkUhUWVmZn5/v6uqq18CfE0yEhq2urq6urk4kEgUGBpqYmAAAwzClpaVtbW1ubm7jx49nd9NqtYWFhX19faGhoZaWlnoN2cBotdrS0tL29nZPT093d3e2kGGYq1evdnV1BQYG2trasoVKpbKgoIDP5wcFBbH/F+gJKRSK8vJyhmH8/PwG94iB/781KhaL2VWZTFZWVubq6urj46OPSA2YXC4vLy8nhLAHWa1Ws28E6YjFYh6Pp9FoioqKjK2twESIEELIqOEzQoQQQkYNEyFCCCGjhokQIYSQUcNEiBBCyKhhIkQIIWTUMBEihBAyajjEGkK/VGZm5kcfffRjW2NjY3ft2vU843k23d3d6enpkZGRugHkEDISmAgR+qU8PDwWLVrELre2tiYnJ8+cOTM4OJgt8fLy0l9oT6Gjo2P58uXp6emYCJGxwUSI0C/l5+fn5+fHLldWViYnJ0+fPv39998fvE9XV1d/fz87D86jNXR0dPB4PN0Ax4SQjo4OOzu7R3fu7+/v6uoSiUS6yeR0enp6enp6RCLRY8e16e7uHjxXlEajkcvlIpHosfE8lZ6enr6+Pt3gL1qttr29ncfj8Xi8X1gzQs8HPiNEaGhlZ2dPmjTJ2tpaIpE4OzvrpoLbuXOnUCgsKyvz9vYWiUR8Pj8hIYFhmKNHjzo6OorFYqFQeOzYMXZne3v7LVu2rFixgq1n9OjRZ86c0X1FUVHRlClTBAKBg4ODRCLZsWMHW56eni4UCvPy8oKDg62srEJDQwHg7NmzL730krm5ub29vYWFRUxMTGNjY3l5OTsj7uLFi4VCoVAo/Prrr/Pz84VC4fXr13VfpCv597//LRQKc3JywsPD+Xw+O1+uRqPZsGGDvb29g4ODQCAIDQ3VDfSM0IsMEyFCQyg7O/vVV1/19vYuLCwsLy+fP3/+O++8w06hoFKpFArF73//+8TExJKSktWrV3/++efLli379NNPv/zyy+Li4pCQkEWLFrHT4igUim3btjU0NOTl5RUVFXl4eMTHx1+7dg0ASktLIyMjhUKhVCqtqKhYuXLlmjVr/vGPfwDAwMCAQqF46623oqOjCwoK2Jl3Wlpa5s+fn5+fX11dffjw4Rs3bsTFxY0dO3b37t0AkJCQkJmZmZmZGRoaqlarFQoFOxcBS1dC07RCoViwYEFQUFB+fv7+/fsBYPny5Tt27Fi7dm1FRcXFixcHBgaioqIUCoU+DjxCT4MghH49FRUVALB9+3Z21dvbOywsjGEY3Q7R0dFTp04lhGzevBkAjh8/zpYzDOPq6mpqavr999+zJeyYyF988QUhxNTUVCKRqFQqdpNSqbSxsXnttdcIITNmzJgwYcLAwIDuKxYuXDhmzBhCSGpqKgBs3LjxJwJmp4Cvrq6+c+cOAKSnp+s2XbhwAQBKSkoeLTl16hQAJCQk6DZVVlYCQEpKiq5EJpONGDFi165dT3H4ENIHfEaI0FCRyWRVVVVLly69ePGirtDR0ZHNIix2FmUAoCjKw8PDwsJCN0mWs7OzpaWlboqAmJgY3UNEHo8XHR1dXl5O0/Tly5djY2Pz8/N1ddra2tbX1/f09LCrs2bNeiiwmzdvnj59WiaTDQwMKJVKALh9+/bEiROf9gcOrjknJwcAbGxs2GTJcnBwqKqqetpqEXrOMBEiNFRaW1sB4Kuvvjpx4sTgcg6H09fXxy7rZoMDADMzs8GrbAlN0+yyRCIZvMnBweHMmTMKhYKm6XPnzl26dGnwVmtra/bbAcDe3n7wpqSkpL/85S9BQUE+Pj42NjZcLhcAurq6nuEHDq6Z/bp33333oX26u7ufoWaEnidMhAgNFYFAAABJSUkJCQm/vLaOjo7Bq+3t7RKJhM/nczicd99999NPP330IwUFBQAwuH8pTdNJSUmJiYmfffYZW1JaWrpnz57HfiObIwc/I3woqw2u2crKCgBu3rz5UMJG6MWHnWUQGipjx44dNWrUiRMnyK8x62dubq5Wq2WXaZqWSqXe3t4WFhaBgYGnT58eGBh4kkrY26GTJ0/WlWRlZbEL7NsO/f39uk3Ozs4AwD47ZOXl5f1YzdOmTQOAzMzMJ/9FCL0gMBEiNFQoitqyZUthYeHbb79dXV2tUqnq6+uPHTuWlJT0DLW1tLT84Q9/aG1tlclkS5YsaW5uXrVqFQBs3rz59u3b8fHx5eXlKpWqsbHx1KlTf/rTnx5bibOzs62tbUpKSkNDg1KpPHToUEpKCrtJJBKJxeK0tLTLly+XlZXJ5XJXV1cPD4/Nmzd/9913MpksOTn56NGjPxbe1KlTX3311bVr16akpMhkst7e3oqKik8++YR9dojQiwwTIUJD6M033zxy5MilS5cmTpzIdoRZtWrVY194/1nvvfdeQ0ODRCJxcnI6efLknj17wsPDASAyMvJf//pXTU1NQECAhYWFi4vLwoULB9/PHIzL5aalpVVXV48ePVogEHz88cd79+5lN1EUlZqa2tbWNnPmzMDAwOzsbBMTk4MHD8rlcj8/Pycnp9TU1K1bt/5EhOnp6W+//fbq1audnJx4PJ6vr29GRgafz3+GH4vQ80T9KjdtEEI6Wq2Ww+EMHrGFYZi6urru7m6xWDxq1KhnSIRmZmYbNmzYsGFDfX19W1vbhAkTHk0wt27dUigUQqFw9OjRpqamP1GbSqWqra01MzMbP378oyPUPEStVtfW1nK53PHjxz9JqL29vXV1dYQQZ2dn3XAzCL3IMBEiZAB0iVDfgSA0DOGtUYQQQkYNX59AyADs2bOHHQsUIfSrw1ujCCGEjBreGkUIIWTUMBEihBAyapgIEUIIGTVMhAghhIwaJkKEEEJGDRMhQggho/Z/z1+DMNbAOqUAAAAASUVORK5CYII=" }, "execution_count": 3, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# Distribution, AKA MZrn from MELTS zircon saturation (meltstzirc) calculations\n", "dist = [0.00102226,0.00105785,0.00109886,0.00113827,0.00118283,0.00122527,0.00127,0.00131584,0.00136402,0.00141291,0.00146517,0.00151593,0.0015771,0.00163667,0.00169765,0.00176478,0.00183148,0.00189186,0.00196526,0.00203994,0.00211244,0.00219109,0.00226588,0.00234415,0.00242908,0.00250853,0.00259262,0.00267285,0.00276366,0.00285942,0.00295,0.00304668,0.0031507,0.00324811,0.00334256,0.00345641,0.00357056,0.00368579,0.00381205,0.00393925,0.0040682,0.00420839,0.00434517,0.00448705,0.00462617,0.00478174,0.00494046,0.00509347,0.005262,0.00543599,0.00559582,0.00576977,0.00596128,0.00615107,0.00635583,0.00656449,0.00678063,0.00700385,0.00723641,0.00747424,0.00772279,0.00797097,0.00824662,0.00852537,0.00880685,0.00910631,0.00941175,0.00970515,0.0100419,0.0103883,0.0107337,0.0111116,0.0114879,0.0118724,0.0122876,0.0127015,0.0131353,0.0135732,0.0140491,0.0145477,0.0150343,0.0155732,0.0161175,0.0166661,0.0172398,0.0178859,0.018519,0.0191786,0.0198861,0.0206008,0.0213348,0.0221012,0.0228635,0.0236548,0.0244204,0.025197,0.0260081,0.0268075,0.0276271,0.0284577,0.0292575,0.0299837,0.0306664,0.0313003,0.0318339,0.0323055,0.03268,0.0329388,0.0330462,0.0330452,0.0329532,0.0326285,0.032175,0.0315716,0.0308023,0.0298608,0.0288332,0.0276613,0.0261425,0.0245284,0.0228927,0.0210703,0.0192917,0.0175039,0.015764,0.01409,0.0124369,0.0109181,0.00945112,0.00806163,0.00679755,0.00566667,0.0045914,0.0036195,0.00290947,0.00229281,0.00177751,0.0013579,0.0010393,]\n", "# Corresponding MELTS temperatures. \n", "T = [809.92,810.521,811.122,811.723,812.325,812.926,813.527,814.128,814.729,815.331,815.932,816.533,817.134,817.735,818.337,818.938,819.539,820.14,820.741,821.343,821.944,822.545,823.146,823.747,824.349,824.95,825.551,826.152,826.754,827.355,827.956,828.557,829.158,829.76,830.361,830.962,831.563,832.164,832.766,833.367,833.968,834.569,835.17,835.772,836.373,836.974,837.575,838.176,838.778,839.379,839.98,840.581,841.182,841.784,842.385,842.986,843.587,844.188,844.79,845.391,845.992,846.593,847.194,847.796,848.397,848.998,849.599,850.2,850.802,851.403,852.004,852.605,853.206,853.808,854.409,855.01,855.611,856.212,856.814,857.415,858.016,858.617,859.218,859.82,860.421,861.022,861.623,862.224,862.826,863.427,864.028,864.629,865.23,865.832,866.433,867.034,867.635,868.236,868.838,869.439,870.04,870.641,871.242,871.844,872.445,873.046,873.647,874.248,874.85,875.451,876.052,876.653,877.255,877.856,878.457,879.058,879.659,880.261,880.862,881.463,882.064,882.665,883.267,883.868,884.469,885.07,885.671,886.273,886.874,887.475,888.076,888.677,889.279,889.88,890.481,891.082,891.683,892.285,892.886,]\n", "# Both of the above must be in order of increasing temperature (=increasing age)\n", "# We are implicitly assuming a constant cooling rate by treating a \n", "# scaled temperature distribution as equivalent to a scaled time distribution\n", "\n", "# Plot distribution\n", "plot(T, dist, label=\"MELTS f_xtal\", ylabel=\"Probability Density\", xlabel=\"Temperature\", xflip=true, fg_color_legend=:white)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Run MCMC to estimate eruption age" ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [], "source": [ "# Configure model\n", "nsteps = 400000; # Length of Markov chain\n", "burnin = 150000; # Number of steps to discard at beginning of Markov chain\n", "\n", "# Run MCMC\n", "(tminDist, tmaxDist, llDist, acceptanceDist) = metropolis_minmax(nsteps,dist,ages,uncert; burnin=burnin);" ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAlgAAAGQCAIAAAD9V4nPAAAABmJLR0QA/wD/AP+gvaeTAAAgAElEQVR4nO3dd1hTd/8//jdJ2HsjSxDEgQoKiDjqRsUtt7Rua7WO26rVqtXqR611tHVvrdaBExVHHShqq1UZ4mLIUATZewUIZJ3fH+d355sERFRISM7zcfXqRV45Oef15pg8OTlLg6IoAgAAwFQsZTcAAACgTAhCAABgNAQhAAAwGoIQAAAYDUEIAACMhiAEAABGU9UgFAgEPB5P2V0AAIDKU8kgpChq3rx5v/76q7IbAQAAlaeSQfjHH38MGzZM2V0AAIA6UL0gLCwsDA4OvnHjxp07d4qLi5XdDgAAqDaOshv4aJaWlv/++29KSsrNmzfNzc2lnyotLeVyuY6Ojp88c5FIxGazP7vHFkddxyUWi1ks1ftj7oPoCx9qaGgou5Gmp66rTF3fYuo6Ljmq+i/Szc1t4cKFcsXr16+vXLnyc2ZbXV39OS9vsdR1XDweTywWK7uLpicQCPh8vrK7aHoURanrMW7q+hZT13HJafYg5PF4I0aMcHR0tLa2HjBgQGRkZAMTi8XiDRs2jB071tvbOy4uji4KhcJFixbZ2tq6uLjs37+/uRsGAABGafYgZLFYs2bNioyMTEhIGDJkSEBAQG1tLSGksrJSerLa2lo+n09RVG5u7qhRo5KSkqqqquindu/e/eDBg2fPnoWGhv7f//1fREREc/cMAADMoaHI2zBVVlYaGhq+e/fO0dGxX79+X3311Zw5cwghtbW1gYGB/v7+CxYsoKc0MTEJCwvr0aMHIcTd3X3NmjVBQUGEkCVLlpSVlR05cqTe+Z88eXLz5s09e/aULk6ZMsXHx6eRHXK5XENDw08eYIulruOqqqrS1dVVv31O9B+F2traym6kiVEUVV1dra+vr+xGmp66vsXUYFxaWlof/IhQ0MEykZGRxcXFp0+fDgwMpA9mCQ4O7t+/v0gkmjVrVlBQkL6+/rx58+q+UCwWp6SkeHh40A89PT0PHDjQwIKMjY0lE9MsLCw0NTUb2aempmbjJ1Yh6j0u9QtCiqIoilK/VUYPSv3GRdT9LabsLj5LYw46U1AQ/vnnn8nJyampqRs2bKArDg4Od+/eHTBgwKFDh9zd3YODg+s9Nqm8vFwoFEr+JDEyMiosLGxgQc7Ozv/9738/uU82m62Wh0ip97jULwjZbDZFUeq3yuhBqd+4iLq/xZTdRbNTUBAeOnSIEJKamurh4dGtW7fOnTsTQmxsbNq2bXv//v1vv/32fb9rY2NjDocj2aFYUVFhYWGhmJ4BAIAJFPqntIuLi52d3Zs3bwghfD4/KCjI2Nj41atXO3bs2Lt3b/39sVht2rSJj4+nHyYkJLRt21ZxHQMAgLpr9i3C5OTk4uLibt26icXi4ODg3NxcX19fQsiECRMMDAxOnDjBZrPv3r3bv39/ExOTSZMmZWZmCgQCsVick5Pz9u1bR0fH6dOnb926dfDgwYWFhcePHz99+nRz9wwAAMzR7EFYWVn53XffJSYmcjgcDw+Pq1ev2traEkKWLVvm7e1NfyPq6Oj4999/03t6vv7667S0NEtLy6VLlxJC7t+/v3jx4jdv3jg4OGhray9ZsqR///7N3TMANN5Pv/y6ZeP6uvWN69ctWbJE8f0AfKxmD0IvL6+nT5/WrdPbhRKS66LduXOn7sRHjhx53ykTAKBcN26G8b89Q9r2li6y/1qvlhfHeR91PaJEXcclR/WuNQoALQ5Hh2gbyFTYmhUVFTk5OdI1+qxltTyPUFdXV9ktNAt1HZccBCEAND0qL+X3u7d3/XFcuiji1/Tq1fPuzWvK6qr5qOUV0on6jksOghAAmh5VkS8OWFEdsFym+uIvfsqfSuoI4L0QhACgOGXFhefOnZMr6urqjho1Sin9ABAEIQAoTlZcYlrWrF0XpWsUv4adHl1WiCAEpUEQAoCiCHjEtj3361MyxfI8o02NvSw+QHNQt4s0AgAAfBQEIQAAMBq+GgUAdfDu3buamhoFLEhHR6d169YKWBAoDIIQABqLoiixWFy3qoxe5A0dNS4zv5il2bx3MxYLau2tzJJinzXrUkDBEIQA0Fg9Bgx78iBcrkjpGJGeSmlHhkBMqmadI627Ne9iMp4LQuc07yLeg75i8+fcfufRo0ddunR5/fq1tbW1nZ1d41/47NkzKysreoPb1dX1oxb68OHDjh07mpmZfVyvioV9hADQWM8i7lN7SqmDNdL/sXXU8JJpH8XFxcXMzMzCwsLDw2PTpk1U4zaRw8LCTp069eHp/ufQoUMnT5781B4JIWTKlCmvX7/++eefw8Pl/5ohhCxevLigoKDeF27YsCEsLOzw4cOHDx/+4FIyMjJ++uknycN58+bFxcV9cs+KgSAEAPgsZWVlJ0+efP369YEDB7Zv3378+PEPv4aQ58+f//vvv41fyqZNm1asWPGpPf4/p06dmjhxYt368ePHKyoq6n1JcHDw1KlTGzn/wsLCM2fOSB4+evSod+/eDUzfEiAIAQA+l6GhoampqZ+fX0BAwOPHjysrK+fPn9+2bVsvL68DBw4QQgQCwYIFC1xdXR0dHf39/SMiIv78889r164NHjx4/vz5hJDExMRRo0Y5Ozv369cvMjKSEHL27NnNmzd/++239vb2t27d+uOPP06fPk1R1N69e7t169a2bduFCxdWVVURQmbPnn3y5Mm+ffs6OTnJNXb9+nVvb29XV9cdO3bQlQ0bNty9ezc6Orpv3752dnbu7u779+//6aefuFzutGnTBg8efP/+/V27du3evfvLL7+0s7OLi4vbuHHj7du3CSHl5eUTJkxwdHQcPXp0dnY2IeTHH3+8f/8+Pec7d+6sXr16wYIFeXl5gwcPHjx4MH0bvri4uKqqqoULF7Zt27Zbt2579+6lKKqiosLf3z80NLRz587t2rXbv3+/YtZUvbCPEACgaQgEgpcvX44dO3bBggUlJSXR0dHZ2dkBAQGWlpZcLjchIeH58+c6OjpPnz7t1KnTxIkTMzIytm3bxuFwSktLBw4cuGfPnhEjRjx8+HDs2LFxcXFpaWnr168/fvz4r7/+ymazw8LCjIyMQkJCtmzZcuPGDWtr6ylTpnz//feHDh2KioqKjo4+efKkubm5dD9v376dNGnS5cuXfX19f/rpp8zMTELIq1ev3Nzcli5dumbNmvHjxxcXF797927KlCn79+/fvXu3s7Ozvr7+uXPnQkJCzp49u3//fh0dncTERDpig4ODL126dPTo0bVr1wYFBT169Oj58+fduv3/+2Vzc3Nfvny5YcOG6dOnh4SEEEL09fWfPXtWXl6+ePHi7OzsqKio/Pz8gIAACwuLQYMGhYeHOzk5PXjwIDU1tX///kOHDnV2dlb0OiOEIAgBAD7f4sWLTUxMUlJSrKys5s6da29vHxsba2pqampqumDBghMnTowbN664uPjFixc9e/bs0aMHIURHR0dbW9vU1JQQcujQIU9Pz/79+3O5XCcnp06dOt29e5cQ0rdv3/Hjx0svKDg4+Pvvv+/QoQMhZPPmzb6+vvQW56JFi9zd3eW6CgkJGTlyZL9+/Qghv/zyy549eyRPsVis2NjY3r17t2rVio5PDQ0NIyMjuh9CyLhx4wYNGiQ3w6FDhw4ePJgQsm7dOjMzs7S0tLq/CkNDQxaLJZkPIYSiqODg4OjoaDMzMzMzs0WLFp04cYKe+YYNG0xNTb29vb28vGJjYxGEAACqauLEiZ6envb29q6urkVFRXw+X3KzcScnp3PnztHbf3PmzMnLy5syZcqWLVukX56RkREXFxcUFCSpcDgcQoiDg4PcgnJyciRnMTo5OfF4vNLS0nqnJITk5eXZ29vTP+vp6VlYWEieOnXq1Lp16zp06ODi4rJly5b+/fvLvbbeGUrmpq2tbW1tLXe/yfeprKzk8XiStlu3bk2/UENDw9LSki7q6+vTX/MqBYIQAOBzeXt79+nTh/7Z3NxcR0cnPT2d3m57+/atvb29pqbm6tWrV69enZSUNGbMmJ49e7LZbMlJmY6Oju3bt5c7mHPTpk0slvxhHHZ2dpLtsNTUVD09PfrMhHpvHGhra/vixQv658rKyqKiIslTnTt3vnDhQm1t7cGDB6dNm5aRkcFisaRPEq13hu/evaN/4PF4dMrq6elVV1fTRTrepMdFMzAw0NfXT0tL69KlCyEkLS1NEqgtBA6WAQBoSiwWa9q0aT/88EN2dnZ0dPTOnTtnzJgRHh7+6NGjqqoqLS0tiqJMTEycnZ1jYmKio6OTkpK+/PLLV69e7dixo7i4ODc39+zZsxkZGfXOfMaMGdu3b3/27FlWVtayZcu+/vrrBu6d+9VXX924cePmzZslJSXLly+XnnLv3r0ZGRlisVhTU5P+GtPZ2fnq1atPnz6lNzHrFR4efvXq1dLS0hUrVnTv3r1169Y+Pj7Hjh3Lzs6OiIigT66wt7cvKiq6c+fO06dPRSIRIURDQ+Prr79eunRpVlZWTEzM9u3bZ8yY8cm/3uaALUIAgM/St29fExMT6cr27dvXrl07fPhwAwODzZs3jxw58saNG2vXrn337p2Jicl3333n7+8vEAgSExN/+eUXR0fHPXv2PH78eM2aNUeOHGGz2T4+PgMGDHB2dtbV1ZXMs0OHDnp6emPHji0vL1+wYEFVVdXQoUPXrFlDCPH19a33jHVHR8cLFy6sW7eurKxs/vz5QUFBhoaG3bp1s7OzCw0NPXLkSHV1tbu7+9mzZwkhBw8e3Llz5+3bt3/66af27dtLH3fTtWtX+pvSX3755fLly0uXLu3atSt9X8mFCxempqYOGDCgS5cua9euTUpKsrCwOHz48KFDh0pLSy9dutS7d28zM7Pff/993bp1I0eO1NfX//nnn8eOHVtRUSG9D9LT09PGxqaJV0yjaTTy3E+VcPLkybCwsM8555TL5RoaGjZhSy2Euo6rqqpKV1e37tdHqo7P51MUpa3dvFcL+wSaOrrC7flEU1e6yF7uLPr6KGnfT7qo8WtfqnMAkbtD/aXV7HdPRIvCZIrleUabfMqL8j6zN9fOXqnDficOHp85nw/Iim1zfUlqPC6xplawRQgA6iBgxIhDu8coYEHmHp4KWAooEoIQANTBrk3rdm1ap+wuQCWp23dKAAAAHwVBCAAAjIavRgFAycQi4a1bt+rW+/btq6Ojo/h+gGkQhACgVNVlVbXCL1dslStXJfz7Ojmx7lWkAZocghAAlEokIIQqn3ddrmyw6tPvQAvwUbCPEAAAGA1BCAAAjIYgBAAARkMQAgB8Li6Xq+wWmoW6jksOghAAABgNQQgAAIyGIAQAAEZDEAIAAKPhhHoAkFdYWNhzUEDG60S5ulAgVEo/AM0KQQgA8nJycgoqePzfM+WfWGCljHYAmheCEADqw2IRbQNlNwGgCNhHCAAAjIYgBAAARkMQAgAAoyEIAQCA0RCEAADAaAhCAABgNAQhAAAwGoIQAAAYDUEIAACMhiAEAABGQxACAACjIQgBAIDREIQAAMBoCEIAAGA0BCEAADAa7kcIwHRPnz7l8/nSlTdv3ohFYmX1A6BgCEIARrt58+Z/ps/WNLOVLgq5JTy+SFktASgYghCA0TIzM4m7f/nE/TLVmAus0FVK6ghA0bCPEAAAGA1bhADQEomFguPHj5uZmcnVhw4d2rZtW6W0BOoKQQgALRGvhrfxzhuWrpFMNeG2paUlghCaFoIQAFoksYg/Yg0xby1dMzw6WVntgBrDPkIAAGA0BCEAADAaghAAABgNQQgAAIyGIAQAAEZDEAIAAKMhCAEAgNEQhAAAwGgIQgAAYDQEIQAAMBqCEAAAGA1BCAAAjIYgBAAARkMQAgAAoyEIAQCA0RCEAADAaAhCAABgNAQhAAAwGoIQAAAYDUEIAACMhiAEAABGQxACAACjIQgBAIDREIQAAMBoCEIAAGA0BCEAADAaghAAABgNQQgAAIyGIAQAAEZDEAIAAKMhCAEAgNEQhAAAwGgIQgAAYDQEIQAAMJqqBuHcuXPXrl2r7C4AAEDlcZTdwKc4fvz4gAEDEhISlN0IgCoRi8VFRUVyRS6Xq5RmAFoO1QvC4uLikydPdu7cOTY2trS01NTUVNkdAaiG7XsPrvhxOUdLR7rIr+Zq+IxXVksALYHqBaG5uXl4eHhKSsrNmzeRggCNFx39RDB+q6D3dJnqkemq9ykA0KRU9S3g5ubm5uam7C4AAEDlNfvBMgKBYPXq1T169HB1dQ0ICIiKimpgYrFYvHPnzqlTpw4ePDgxMVFSXL16tZubm6enZ3BwcHM3DAAtF0VxudyiOqqqqpTdGaiwZg/Cmpqa/Pz833///e7du3379h0yZEhBQQEhpLa2VnoyoVAoEonEYvGzZ898fHweP35cXl5OP3XgwIFLly7duHFjz549CxcujImJae6eAaBl4hdlzlvwvYNre+n/bB2clq/4SdmtgQpr9q9GDQ0NDx06RP+8fPny7du3v3jxwt/ff9iwYTNmzJg8eTIhRCgUTpo0qVevXgsWLDh+/DghZPXq1ZI5HDx4cOXKla6urq6urlOmTDl06JC3t/f7FpeQkLBkyRLpyvjx47t27drIbmtra7W0tD52jC2fGo+LxWKxWKp6FtD78Pl8iqKafLYisYiwm3yuTaDxQxVVlginHhT6BMlU7+ziC9Pl/rZWMDV+i6n6uLS0tDQ0NBqeRqH7CN+8eVNaWtqhQwdCyKFDhwYOHCgWiydNmjRt2rTKysrZs2fXfYlYLE5MTPTy8qIfenl5SWK1XlpaWhYWFtIVXV3dxn9QquWnKlH3canf0FgsFkVRTT6uD34cqC4NDQ3l/jNQy3+HRH3HJUdxQVhdXT1hwoRly5Y5ODgQQlxdXe/cuTNw4MADBw5YWlpeunSp3r87ysvLBQKBoaEh/dDY2Jj+ZvV92rZtu2LFik9uUlNTU1NT85Nf3mKp97jU741KURRFUU2+ylgaLfQXpfExG4X1YrFYyv0Xrt5vMWV30ewU9MaoqakZPXq0u7v7unXrJEVnZ2dPT88nT56MHj36fVvfRkZGbDa7srKSfsjlcs3NzRXRMQAAMIMitgj5fH5gYKCFhcWRI0ckf7yLRKJp06aJxeL4+PiAgAAOhzN16tS6r2Wz2a1bt05KSqJPlkhKSmrTpo0CegYAAIZo9iAUCATjx48XCoW7du2qqKgghOjr62tpaU2dOrW6uvrSpUuampq3bt0aOHCgrq7u+PHjKyoqRCIRRVFcLre0tNTY2HjatGnbt28fMmRISUnJ8ePHjxw50tw9AwAAczR7EBYWFsbHxxNCevToQVd+//33cePGzZw5s3fv3vS3z66urnfv3qU3FseNG5eWlmZhYTFnzhxCyP3793/44Yf4+Hhra2uKohYsWODv79/cPQMAAHM0exDa2tqmpqbWrffv31/6oaurK/3DnTt36k4cEhLC5/PZbDab3SKP/gYAAJWlMpdYU/VzWQAAoGVqoYdTAwAAKAaCEAAAGA1BCAAAjIYgBAAARkMQAgAAoyEIAQCA0RCEAADAaAhCAABgNAQhAAAwGoIQAAAYDUEIAACMhiAEAABGQxACAACjIQgBAIDREIQAAMBoKnM/QgBovIKCgri4OLlifkE+aaWUdgBaNAQhgBoa+dW0xJxStraedJGbFkvGjVFWSwAtFoIQQA1lZ77jTjlLbDtIFznruiqrH4CWDPsIAQCA0RCEAADAaAhCAABgNAQhAAAwGoIQAAAYDUEIAACMhiAEAABGQxACAACjIQgBAIDREIQAAMBoCEIAAGA0BCEAADAaghAAABgNQQgAAIyGIAQAAEZDEAIAAKPhxrwAoPJysrPv3LkjVzQ2Nvbx8VFKP6BaEIQAoOIyX95Iirj/tkS6JqquaG3Ijn8aqaymQIUgCAFAxfEqKM9R5RN2yhTfRgnDflBSQ6BisI8QAAAYDUEIAACMhiAEAABGQxACAACjIQgBAIDREIQAAMBoCEIAAGA0BCEAADDaB06o/+9//1tYWNjwNCEhIU3XDwAAgEJ9IAjj4uKys7PrfUosFqenpzd9RwDwMcLCwiIiIuSK3IoKpTQDoIo+EIQPHjyot3716tWffvqJEDJw4MCmbwoAGqeiomLEqDGiYcvk6hqVVUrpB0AVffS1Rh8/frxixYoHDx74+PhcvXp15MiRzdEWADQSW1NLNHK1fPXeXmX0AqCSPuJgmfj4+KCgoF69euXn54eEhERFRSEFAQBA1TUqCNPT02fPnu3h4REREXHw4MH4+Pjx48draGg0d3MAAADN7QNfjRYWFm7dunXHjh0GBgYbN25csGCBrq6uYjoDAABQgA8EoZ+fX2pqas+ePZcsWWJkZPTo0aO60wwaNKh5egMAAGh2HwhCoVBICHn8+PHjx4/fNw1FUU3cFAAAgKJ8IAhPnjxZU1OjmFYAAAAU7wNB2Lt3b8X0AQAAoBS41igAADAaghAAABgNQQgAAIyGIAQAAEZDEAIAAKMhCAEAgNEQhAAAwGgIQgAAYDQEIQAAMBqCEAAAGA1BCAAAjIYgBAAARkMQAgAAoyEIAQCA0RCEAADAaAhCAABgNAQhAAAwGoIQAAAYDUEIAACMhiAEAABGQxACAACjcZTdAABAM6DEpSXFCxf/IFdmaZDtW7copSNosRCEAKCOchILy7i7Mi3l6xd+RBCCHAQhgGrIyckZNGp8fk6mdJGixAKBQFkttXAaBubUkMXy1Qs/KqMXaNEQhACqIT4+PpvHqlj0j0y1PJ/83l85DQGoCwQhgMrQ0NIhZg6yJRzvBvC58C4CAABGQxACAACjIQgBAIDREIQAAMBoCEIAAGA0BCEAADAaghAAABgNQQgAAIyGIAQAAEZDEAIAAKMhCAEAgNEQhAAAwGiqGoQLFy5cu3atsrsAAACVp5JBeObMGW9vb2V3AQAA6kD1grCkpOTw4cPPnz9/+PBhaWmpstsBAADVpnr3IzQxMblw4UJKSoqRkZGxsbGy2wEAANWmekHIYrFMTU19fX19fX2V3QsAAKg8RQThrVu37t27l5aWNnny5FGjRjUwJUVRR48ejYmJKSoq2rBhQ9u2benib7/9dvbsWX19/SVLlowdO1YBPQMoV1JSkqampnQlIyODEouV1Q+AGlNEEIaHh2tqaqakpCQlJUmCUCgUcjj/b+lisZj+/7Vr17y9vY8dO7Z48WI6CI8cOXLkyJHz58/n5ORMmDChTZs2Hh4eCmgbQFkuXbo0bfZ32sYW0sXaimKRpYuyWgJQY4oIwi1bthBCXr16JV0cNmzY3Llzx40bRwgRi8XffPNNt27dvvvuu9DQUELIb7/9Jply3759q1at8vDw8PDwmDx58sGDB/ft2/e+ZSUnJ69Zs0a6MmrUqE6dOjWy1draWi0trUaPTGWo8bhYLBaLpXrHfDUsLi6O32dW7ei1MtUbv7LjbiinoWZGKXZxtbW1zTFPdX2Lqfq4tLS0NDQ0Gp5GaZ8gu3fvXrBgwZkzZ8Ri8YwZM7Kzs2fOnFl3MoqiEhISfHx86Ic+Pj5xcXENzFYsFgtkiUSiZhkAAACoBaUdLNO+ffubN28OHTr0yJEjLBbrypUrurq6dScrLy/n8/mSo0NNTEwKCgoamG2HDh02btz4yV3x+Xxtbe1PfnmLpa7jEgqF2tra6rdFyOFwCBEouwvF0VDsRmFzvBfU9S2mruOSo8xPEHd3dz8/v3v37k2fPr3eFCSEGBgYsNnsqqoq+mFlZaWpqakCewQAADWntC1CiqLmzZvH5XKfPXs2cuRIFov11Vdf1Z2Mw+HY29unpKTQB86kpKQ4OTkpulcAAFBfitgirK6uLi0tFQgEPB6vtLSUz+cTQmbNmpWWlnblyhVPT8+bN2/+8MMPV69eJYRUVFSUlpZSFMXlcktLS8Vi8eTJk3fv3i0SiUpKSk6cODF58mQF9AwAAAyhiCBcuXKli4tLZGTkzp07XVxc6ONCx4wZc+XKFR0dHUJIp06dwsLC6G2+gQMHuri4sNnsL7/80sXFJTMzc/ny5SwWy8bGxsXFZcyYMcOHD1dAzwAAwBCK+Gp0x44dO3bskCuOGDFC+qHkDIcnT57UncONGzfKysq0tbXftysRAADg06jMJdZMTEyU3QIAAKghdTvuHAAA4KMgCAEAgNEQhAAAwGgIQgAAYDQEIQAAMBqCEAAAGE1lTp8AAGgCOob9A8bI1Tgssm/77/Q1PYCBEIQAwCQ13H+cp8jV9K78WFJSopR2oCVAEAIAw3QdJVfQvPurUhqBFgL7CAEAgNEQhAAAwGgIQgAAYDQEIQAAMBoOlgFQpqSU17v27JMrxkRHUhb9ldIPAAMhCAGUaeZ/Fz2qbUVatZepphUTCyU1BMA8CEIAZfMcSboEyFRibyipFQAmwj5CAABgNAQhAAAwGoIQAAAYDUEIAACMhiAEAABGQxACAACjIQgBAIDREIQAAMBoCEIAAGA0BCEAADAaghAAABgNQQgAAIyGIAQAAEZDEAIAAKMhCAEAgNEQhAAKMjRwIkdLW+6/x48jlN0XANPhxrwACvLv3VuiTa+Jvpl0kf1/nUTKaggACCEIQgCFYmsRjrZsSUM5nQDA/+CrUQAAYDQEIQAAMBqCEAAAGA1BCAAAjIYgBAAARkMQAgAAo+H0CQBgOj6veuyEaVraMme2aGiQDauWT5w4UVldgcIgCAGA6WrLi3LH/0Zs3aWLWuFbCgsLldUSKBKCEACAECsX4tBFuqBhYKGsXkDBsI8QAAAYDUEIAACMhq9GAZpeQUEBl8uVK1KUWCnNAEDDEIQATSw/P79Ljy+qa4VydR6vRin9AEDDEIQATayiooInEFf+nCT/xDwjZbQDAB+AfYQAAMBoCEIAAGA0BCEAADAaghAAABgNQQgAAIyGIAQAAEZDEAIAAKMhCAEAgNEQhAAAwGgIQgAAYDRcYg0AoF9nrpgAACAASURBVB5UbdWfp0LuRz6Vqw/o4zd/3lyltATNBEEIAFAPUUZsrLVbrF5/merbKP69BwhCNYMgBPgsu3btKikpka6UlJTw+Xxl9QNNydmb+E2WqbA1Scl1JXUDzQVBCPDpgoODl249zPcYI1MtLNEQCJTUEQB8NAQhwGfRat2FP2q1TCnxLkm6p6R2AOCj4ahRAABgNAQhAAAwGoIQAAAYDUEIAACMhiAEAABGQxACAACj4fQJgEYpLy9fv3GTXPFVQoJQqK+UfgCgqSAIARrl++Urj0Zlkja+MtVX+RwLRyV1BABNA0EI0GgdB5H+sheZLHpHeCXvmRoAVAP2EQIAAKMhCAEAgNHw1SgAQKOV5YTfvq2tZyBXtndwTE1+pZSO4PMhCAEAGq34nbhND+Gsk7LFjPKDo5XUEDQBBCEAwMdgc4i27Bahlp6SWoGmgX2EAADAaAhCAABgNHw1CiDvz1Mhew8clCtmpKeTL9oppR8AaFYIQgB5m3797U2X6cS6rXSRlbpAWf0AQLNCEALUx8mbOHlJFzRwQASAmkIQAnNRFDVm0ozIhw/k6iXlFUrpBwCUAkEIzCUSif46F0z9kiBXZ23uo5R+QHWJxaKUlJS6dRcXFzabrfh+4KMgCIHxLNvUKWkooQ1QXVUl5ZXV3gNHyJdz3+ZkZ1tbWyulKWg8BCEAwOfhV1MsNneN/FcLusvsldIOfCycRwgAAIyGIAQAAEbDV6PAFNevX8/KypKuiMViQihl9QMALQSCEBghOjo6aMY80slfpioWUchBAMZDEAIj8Pl8TQv78gn7ZKrCGnL/uJI6AoCWAkEI6qa8vHz2giV8Pl+6WFxYUFtTo6yWAKAlQxCCuvnjjz8uvMgWeQfJVKv/ZgtzldQRMBVFFRYW1i2bmJhoa2srvh14HwQhqCGWbQeR32SZUm0lyX6ppHaAoWoplnfvfiyWzJVlBFXlF8+HjBo1SlldQV0IQlBhwecu/rRqtYbsdWB4lRVi36lK6gjg/xHzuLWbkomRzJVljA+NU1Y/8D4IQlBh50MvlQxcRjoNlake/xanx0KLRYlFycnJkZGRcnU7OzsHBweltAQIQlANRUVFaWlpcsWKigriYEHMHWWqmjqKawvgI9UUZq7edkDnyAXpYm1J7sJvJm/etFFZXTGcigVhTU3N9evXKysrBw8ebGtrq+x2oOkJBIJla355+vSZXP1FXDylbcDWkjnEgJuTRtwV2BzAZxPXVAmDttZ6yF6h+8ZmivCU1BGoVBDyeLxevXqZmpo6OTktXrz4zp07Xbt2VXZT8OnS0tL6B4wpLcyXLlJicSWPR80MlptY48EUat09YiZzFWPO6o7iZm8TQBGKi4sTExPlivr6+o6OjvVOD01IlYLwzJkzbDb79u3bbDbbwcFh/fr1oaGhym4K5JWXl8fFxckVxWLxu3fvzM3NpYuPHj0q1GlVveKazKTZ8RoHJxGP4XVmjFsjgfrKTTwWHhZyS+Ye0aKa6o6uTrf+uiQ3rYaGhqmpqQKbU3+qFITXr18fM2YMfZfLcePGbd68maIoDQ18PjYxsVhM1bnyWFhYWHp6ulzxzJWbuTnZchGVX1xKaelpGphIF6sykzX0TfVsXaWLNTnJIks3YmwjM9OynM9rH0AFVRSK+8zkjtsgU4w4+TzkB1snV7lpKT4v/PatuvPw8vLS19dvvh7VmCoFYU5OzsiRI+mf7ezsamtri4qKLC0tpadJTU3dsmWLdGXIkCHt27dv5CIEAoFAIJArnj59Oi8vr+6UpaWlrVq1ki7m5+dTRIPDkfmt1tTw0tPS/Pz85Obw6NEjT09PubtXJyQkmJlbmJjIpEh2dnYtn+/gIPMNSUV5+duMDCcnZ+miSCR68uy5lbVstBCSnv7OzMxMU0tLulhcUpKX+Y6jKdOtSMCv4nLFglq5ORA2h4iE8kUNFukeRDRkDtLUeHeL6h5ALGQa0yjcR1m3K3ftJ/PyKh6r6B0J3ylTLM2mhLXyRUIoiiIP/yS6xtJFcW01ib1BSjJlJi1Ko2q48nN4EyHmFskXxUJCqHqWJRaR6LPk9UOZYlUpSfqHCGUuWEOy48X8Kvk55L8htXWKtPuHiJaezGwFNeTFXyT/tcxkJVmEFS0/h9THVFmufJFXQUSiepYlFpOIk8TISqbGqyCv7hBeucyUuUliFkt+DlmxFK+8/iHc3UvYmjJDEArI01CSKXuaZnkeSY2Qn0PaE6o0U75YnksJ+e9Z48eInsx7QVxTReJukjLZayMUpIor66zct9FUZUn9Q6i7LJGIRIeQNxEyxcpiknyfiGW/fc98KeaVyc+hMI3i8+pbFkUeHCbaBjJD4NeQF3+RglSZCYszCCWWn0PyA0pbv3bgApmigEdu/ta3b1/5RekamxgZcjiy5yzW8ls5OFpaWkgXqyu5ViaGvt19pIs8Hi85OdnLy0turv/++6+3t7fcx9TzuHiegNLXl/mXXFJU5OLQqkOHDjLFkpK8vLyOHTtKF4VCYUxMTPfu3eWHQMjMmTPlPv0+H4fD+eD2kioFoVgsloyHXjFCofxHc3V19bt376QrXC5XJBI1chEikajuxLm5uVnZ2XJFAZ9fUlIidx2vgoICPSMTLdm84fF4+QWFmbL3PSCEZOfk2tjZy/0Lyy8sZmnqULJn4BaXlmmwWJXV1dLFispKAZ/PYcueKUCJtVnEzVl+pwK3uMC9nauOjsyRJqmvX3d2c5F7k+Tl5vJqapydZWKsqqoqNTW1S5cucrONjIjw7WYh928sQauTg22VkZHMWnjn01FTS8fWQqZY4mlbVKTtJlvkG/HjPD28ZIuEkGeendytCrS1y6SLiR7trKyF5uYyE+f4tBcIha1l58DVMsgwcHOXLVIUFeXZrUedZcV5dnRuVWFgIFNP83bX0+dYy05c1M2hrKzCVbZYYyBKqu3kWWe2MR5dPK3y5D6qkrp2tGlVa2IiM3Fm9/ZiisgNoayTWa65SwfZokgketq1a/c6y3rp6e7WqlRXV+b4i1SvjkbGGpayExd4t6msqm4jW6zW5bzR6NilzmyjPT28LbNZLJlV/qprR3v7ark1ntHdnc3RspNb4x7WeXmko2xRYMyP7VrPGn/u0amjdaG2tkxsv+7W0cxCLLfGc33c+AKB3K+rUlsvXbddpzqzjfDo6lfvGrflyq3xdJ+O2jpareTWeFf70lKDtrLFWkN+XKW7dz1rvLOHVb6mZrF0MalrR2sbvqmpzMTZvh1EYuIoO4cKTZMsU1e5X5dYLH7Ssb1vjx5yy3r54kVbt3Z6errSxbepbwwMja2sZLYWCvLzqyor8wqLpIvV1dV5BYU5eTI77Akh+YVFBSVlLJbM50x5WZmtnYOxsZF0UVDNLS0rk/ugo4PQ0Eh2SoEgKzvbNkf+ux8NQvh8fuM/rhuJzWarVRDa2NgUFBTQP+fl5XE4HCsrK7lpOnfuvHv37k9ehEAg0NGRP/j+xx9//OQZthBcLtfQ0FDZXTS9qqoqXV1duXepGuDz+RRFqd9VuCiKqq6uVsuv79T1Laau45KjSp8gAwYMCAsLo3++detW37595TanAAAAPpYqbRFOmzZt27ZtM2bMaNOmzdatWy9evKjsjgAAQOWp0hahiYlJTExMp06dRCLR33//PWDAAGV3BAAAKk+VgpAQYmlpuXjx4jVr1nh6ejb5zMVi8dWrV5t8ti3BnTt3Kisrld1F04uKisqucxyTGkhKSnr16pWyu2h6OTk5da+xqQaqq6tv376t7C6axbVr1+oek6h+VCwImxWfz//mm2+U3UWzWLp0aU6dY7TUwK5du9TygzU0NPT8+fPK7qLpPXnyZPv27cruounl5+cvXrxY2V00i1mzZvF46n/tNwQhAAAwGoIQAAAYTaPuxbRU182bN1etWiV3tZfGE4vFERERvXr1atquWoKYmBh3d3ddXd0PT6pSEhMTLS0tLSwsPjypSsnIyKAoqnXr1spupImVlJTk5ua6u6vbHUNqamri4uJ8fHw+PKmqiYiI6N69u0qfqLZ79265K4TUpVZBSAj5+++/q6qqlN0FAAC0CF988YWR7KVt6lK3IAQAAPgo2EcIAACMhiAEAABGU6VLrH2y6urqZ8+ePX/+XEdHZ9asWXSRz+cfPHgwPj6+Xbt28+bN09HREYlEL1++jImJKS0t/f777+VuIkH7559/oqKi5s+f3xIuHCwSiWJjY2NiYkpKShYuXCi5XPjFixfv3Lljbm4+f/58GxsbQkhKSsqTJ0+ysrK+/PJLJycnufk8ePAgNDRUIBB079592rRpCh5FvVJSUmJiYjIzM4OCgiQ7uiMjI0+fPs1ms2fMmNG5c2dCSE5OzpMnT1JSUnr16tWzZ0/pOVRWVh4+fDg5Odnc3Hz69OmurvI3dVOK3Nzc6OjolJSUnj17Sg7LSk9P379/P5fLHT169JAhQwghL168uH79ekZGhqOj44wZM+SO/7p8+fKtW7fMzMz++9//2traKmEYdfB4PPotxmaz586dSxcFAsEff/zx8uVLV1fX+fPn6+rqFhcXX7p0KTY2VlNTc8SIEf3795fMISMj49KlSykpKUZGRhMmTKh7qxOlqKqqiomJefHihYmJieStUVNTs3fv3pSUlM6dO8+ePVtTUzM3N/fy5csJCQl6enqBgYG+vr51Z3Xp0qV3794tWrRIsSOon1AofPHixdOnT8vKypYuXUpftp6iqLNnz96/f9/Gxmb+/PkWFhaVlZXXrl2LjIykKKpfv35jxoyRu43DlStX6DulDxkyZPjwujfTVh0UAxw/frxTp059+vRxc3OTFKdOndq7d+8TJ04MGjRo3LhxFEW9ePGiTZs2AQEBhJCKioq688nNzaU/T/Py8hTX/fslJCQ4OzuPGDGCEFJcXEwX9+/f37p16z///HP27Nlt2rTh8XgURTk5OQUEBBgYGNy+fVtuJnv27LG2tt64cePBgwe///57RY/hPegVYWhoePPmTboSGRlpaGi4Y8eOjRs3GhkZJScnUxQ1cuTIvn37Ojk5rV+/Xm4OvXv3Hj58eEhIyMqVK42NjenjMJVuzJgxX3zxhbOz89q1a+lKcXGxtbX1kiVLDh06ZGVldeHCBYqiOnfuvGTJkj/++GPq1KmWlpbZ2dmSORw+fNje3v7IkSPz5s1r3bp1ZWWlckYi68yZM+7u7l988UXr1q0lxVmzZvn5+Z04cWLYsGEjRoygKGrlypX/+c9/9uzZs2nTJlNT00OHDkkmnjBhwjfffHPgwIGVK1fq6ur+888/ih9FXXv27PHw8PDz8/P29pYUx44dO2jQoBMnTvTp02f69OkURc2dO3fSpEn79u1bt26dgYFBaGio3HxevHhhZWVlZGSk0O7fLzIy0tXVdejQoeR/tzqhKOrXX391c3M7evTo1KlT3d3d6b9j/P39t27dumvXLkdHx2XLlknPZOHChe3atdu2bduePXvqvgFVCyOCkHb58mVJEGZmZmpra9N5Vl5erqurS3+wUhSVlZX1viAMDAw8fPhwywlCWn5+viQIRSJRmzZt/vrrL/opLy+vEydOSKZ0cHCQC8L8/HxdXd3nz58rsuHGc3JykgRhUFDQqlWr6J/nzJkzf/58yWSBgYFy78OSkhJCSFZWFv2wa9euwcHBCmm5Ub788ktJEG7dutXf35/++ejRoz4+PhRFCQQCycSenp4HDx6kfxaLxW5ubhcvXqQf9ujR4/Dhw4rr+0Nu3rwpCcK8vDxtbe3MzEyKoiorK/X19ePi4qTHtXPnTj8/P8lD6afmzp379ddfK6jpRggODpYEYVJSkq6ubnl5OfW/MWZlZUk3v3r16pEjR0q/XCAQ+Pr67t27t+UEIe3169eSIOTz+TY2NvTfHyKRyM3N7fLly9Ljun79uqWlpeThkydPDA0NCwoKFN92c2DoPsKoqKh27dpZW1sTQoyMjLp16/b48eOGXxISEkIIGTt2rCL6+1S5ublv376VfOPUv3//R48eNTD9gwcP3NzchELhL7/8cujQoWrZe/+2KI8ePZJcZn3AgAENj8vY2NjFxeXhw4eEkMzMzMzMTA8PD0V0+fHkxhUTE1NbW8vh/P/7LMRicUVFhZmZGf2wqKgoJSVFsn4/+HtQoidPnjg5Odnb2xNC9PX1fX19Hz16JBkXIaSsrMzU1FTysIGnWpTHjx9369aNPhzf2tq6Xbt2kZGRDTf/66+/DhgwoFu3boru9WO8efOmpKSkd+/ehBAWi9WvX7+G11d4ePiQIUNiY2PXrVt36tSpJr+broIxNAjz8vKkz8K2srLKzc1tYPri4uJVq1bt3Lmz+Vv7LLm5ubq6upL9l5aWlg1fYjQ9Pb2goGDJkiXGxsbXr1/v2bMnn89XSKcfRywW5+fnS1aZpaVlw+uLxWKFhoYuWrTI3t6+Xbt2GzZsoPcptkC5ubnS46IoKi8vT/Ls77//rq+vP3r0aMnEHA7HxMREMn2LvYRsw2+xN2/ebN++feXKlXVfeO/evRs3bixcuFARXX68hsf17NmzY8eOLVu2TFJJSko6d+7c6tWrFdrlx8vLyzM1NZWcOC/3T6ukpGTlypWrVq2SVNLT06Oionbu3GlmZrZ///7AwEBFd9ykGHGwTF3a2toCgUDysLa2tu6N6aV99913y5Yts7Ozo79za7F0dHTobzPofdp8Pr/hq8lwOJyioqK4uDhzc/O5c+e6uLhcv369BW71slgsLS0tySrj8/kNr6/q6urAwMAlS5ZMnjw5ISFh0qRJnTt39vPzU0izH0dbW1vyxwf9g2RowcHBe/bs+eeffzQ1NekKfUiXSCSi/1T/4PpVogbeYllZWUOGDFm/fn3dqzg9e/ZswoQJZ86ccXR0VFyvH6OBcaWkpIwcOfLgwYOSS+eIxeKZM2du3769xa4mCblxSf/T4nK5w4cPHzt27JQpUyQTcDgciqIuXryoqak5efLkVq1aJSQkqO41gxi6RWhnZ0fvC6RlZWXZ2dk1MH1oaOjGjRtdXFy8vLwIIb6+vmFhYc3e5ceztbUVCoX0XkNCSGZmZsNHFdrb25uZmZmbmxNCOBxOmzZtWuwWhvQqy8zMbHh9PX78mMvl/vDDDzY2NgMHDhwxYsS5c+cU0uZHs7e3lx6XlpaWpaUlISQkJOTHH3+8ffu2i4uLZGL68FHJOvrg+lWi973FcnNzBw4cSO/llXtJbGzs8OHDDxw4MGzYMIX2+jHeN663b98OHjx4/fr1EyZMkH42Kirq22+/dXFx+c9//lNZWeni4pKUlKSEvj/Ezs6urKyMy+XSDyX/tKqqqoYPH+7h4bFt2zbp6e3t7Z2dnek/0UxNTS0sLFT6hmgMDcK+ffuWlZVFREQQQuLj41NTU/39/eud8t9//01NTU1MTLxz5054ePjFixcJISEhIX369FFox41jZmbWt2/fM2fOEEK4XG4Dm3cJCQkxMTH+/v41NTX0O7OsrCw2NraFHLZe15gxY06fPk3+d5D3mDFj6p0sJycnPDzc2tq6rKyM/s6KoqjExET6NJIWaMyYMaGhoTU1NYSQM2fOjBo1isViXbp0af78+X/99VeHDh3oyUpLS69cuWJoaDho0CB6/VZWVl69erUFbr7TevXqJRAI/vnnH0JIcnJyQkLCsGHDCgsLBw8ePGXKlKVLl0qmvHfvXkZGRnJyckBAwLZt21rsiGj+/v6pqanx8fGEkIiIiLKysn79+mVkZAwePHj58uUzZsyQTBkWFsZisZKTk8PDw8PDw3ft2qWnpxceHt6mTRvltf9erVu39vT0pP9eLC4uvn379tixY3k83qhRo1xcXPbt2yc5cSIqKioxMTEwMDAhIYH+hiwlJaWoqKhjx47KHMBnUuqhOgoSFxfn5eXl4uKio6Pj5eU1c+ZMiqIOHjxoZWUVFBRkY2Ozfft2iqKEQqGXlxedBJ6enr169aIoys/Pb8uWLZJZFRcXk5Z01Ki3tzd9GIiHh4evry9FUY8fP7a0tBw7dmz79u2/+uorerLJkyd7eXlpaWm1bdvWy8srJSVl8eLF//nPfyiK2r9/v42NzaRJk9q0aTNnzhxlDkbKtGnTpBtOSkrKzc11cXEZOHBgr169PD09y8rKKIravHmzl5eXqampra2tl5fXpUuXzp496+DgQFHU7NmzbW1tp06d6uXl1blz55KSEmWPiaIo6vfff5du+MKFCwKBwN/fv0uXLqNGjbKxsYmPj6coSl9fv1WrVl7/s2/fPvqPNrFYHB0dTa/fDh06BAYGisViZY+JoigqKSnJy8urbdu2WlpaXl5eU6dOpSjq6NGjlpaWQUFBrVq12rx5M0VR3333naampmRcw4YNoyiqQ4cOR44cGTRokJGRkeSpFvJPMTIy0svLy9nZWU9Pz8vLa+HChRRFbd261cbGJigoyMrKij4D5KuvvqInoNHvOwsLiytXrkhmFRER0XKOGq2qqvLy8urUqRMhpFu3bgMGDKAoKjw83MLCIjAw0MXFZdasWRRF7dixQ0NDo2vXrpKh1dTUjBgxYsWKFRRFLV261NnZmf5e9LffflPykD4PI641WltbK73ZrqOjQ2/1v337NiEhoX379m3btiWEUBSVlpYmmUxDQ8PZ2Tk9Pd3IyEhy2J5YLE5PT2/dunULuRz727dvJT/TDRNCiouLHz9+TH+S0n/HZWdn19bWSqa0t7cvLS0VCAT0QX3p6enx8fGurq7t27dX+AjqJ9ewnZ2dtrY2j8d7+PAhh8Pp1asXfbmDoqKiiooKyWT0wSb5+fn0Ck1JSUlJSbG0tPT29m4h66u4uLi8vFzy0MLCwsjISCwWR0ZGlpWV9erVy9jYmBCSlpYm/cY0MTHR1dVNS0uj/+guKSl5/PixtbW1t7e33AnOyiL3FtPW1qa/MExPT4+Li3Nzc2vXrh2ps744HI6jo2Nqaqq5uXl1dTW9WUzT1dX95NvINCEejyd9LIyenp7kChXJycmdOnWi33H5+fnS1/rX0tKyt7dPSUmxtbU1MDCgi7W1tbm5uXUvZ6EU9OeY5CGbzaZvdZKfnx8dHe3g4ODp6UkIKS8vp//0l3B2ds7KytLR0aG/wH/16tXbt2+7dOnSYvfpNhIjghAAAOB9GLqPEAAAgIYgBAAARkMQAgAAoyEIAQCA0RCEAADAaAhCAMV59OjRyZMnld1Fk7lz505iYuJnzkQoFJ4/f76wsLBJWgL4BAhCgIYkJia2adOmqdLr3Llzy5cvb5JZNezVq1cuLi6enp7NdxX1ly9fSu7FOmLECBcXl8mTJ8tNs2PHDhcXFxcXF+nT7ORwOJwTJ078+OOPzdQnwAchCAEaUltbm5aWJn0WvEo4evRoTk7Oy5cvr1271kyLWL58eWBgIH0RuKysrJycnLNnz0pfq5aiqH379uXk5Lx9+1YsFjcwq1WrVh07diwhIaGZWgVoGIIQ4LNQFFVQUFBWVlbvs1VVVXl5eQ3HQOMJBIL8/HzpuwTUSygUnjx5ctq0aV26dDl69Gi90xQWFlZWVtb7FI/Hy8vLa3hT8s2bN7dv354+fbqk4uvr26pVq+DgYEnl4cOHr1+/HjlypPQL6V+X3BB8fX07dOiwf//+hscF0EwQhACf7ujRo46OjtbW1qamph4eHvQ1pmmPHj3q2bOnkZFRq1atdHV1671h299//21vbz9v3jyRSOTv7y99mxtCiKTi5+c3d+7c9evXm5ub29jY2NjYHDx4sIGurl+/npeXN3Xq1ClTpoSFhcndUeTp06ddunSxsrIyMjIKCAg4d+6cmZlZbGwsISQ9PX3UqFF0z6ampvPnz5e+0J20Y8eOGRkZDRw4UFJhs9kTJ078888/pafp06cPfXE1Qkh8fPygQYN0dXWtra319fV9fHwiIyMlE48ZM+bkyZNCobCBcQE0EwQhwCc6derUjBkzevfuHRUVde/ePR0dnaFDh758+ZIQEh0dPXDgwJqamsuXLyckJFy5cqXu7S8uXLgQEBAQGBi4Z88eNpvN5XLldqRJKhUVFefPn798+fJff/319OnTYcOGzZkzJzQ09H2NHT161NXV1c/Pj85R6R2chYWFQ4cOFYlEt2/ffvnyZZcuXebPn19aWioSiYqKivr06ZOWlhYaGpqQkLBr166TJ0/Omzev3kXcu3fPx8dH7gqu33zzTUpKSlRUFCGkqqrq/Pnz0puMRUVFXl5eN2/eTExMvHbtGofDCQgIkFzK0s/Pr7y8/NmzZx/8tQM0PaVd7htAFTx//pwQsmfPnrpPubm5de7cWXL/h6KiIgMDgy+//JKiqAEDBlhZWZWXl8u95LvvvrO1taUoaseOHRwOZ/fu3ZKnevToMXbsWOmJJZWOHTtqaWllZ2fTdbFY3LFjRy8vr3obzsvL09TU/Pnnn+mHw4YNc3NzkzS5efNmFouVmJgomX706NGEkGfPnq1YsUJPTy8rK0vyFJ3Q+fn5dZeir68/b948yUMPDw/6Dga+vr5z586lKOr48eP6+voVFRX0bc0rKirk5pCfn89isYKDg+mHr1+/JoTQN3MAUDCG3qEe4DNVV1e/fv36559/ltz/wdzcfNCgQQ8ePKitrf33339nzZplZGRU94UikWjevHnHjh0LCQlp/L33evToIbkHr4aGRmBg4ObNm0UiUd27apw4cUIoFE6aNIl+OHXq1AkTJkRERPTs2ZMQEhsb27ZtW+nbjIwcOfLKlSuEkFu3brVp0yYxMVFyRoSWlpZIJEpMTLSyspIbe1VVleSWLNKmTZu2YsWKrVu3Hjt2LDAw0NDQUPrZvLy8kJCQtLS06upqev5v3ryR/PYIIQUFBY38hQA0IQQhwKfIyMigKEruVkG2trZFRUVlZWUCgcDBwaHeF5aXlx89erR///6jRo1q/OLkvllt1aoV7VIgqwAABT5JREFUfeBM3TvUHzt2rFWrVnfu3KEf1tTUaGhoHD16lA7CvLw8+gY6EpKHBQUFBQUFQUFB0s+ampoWFRXJLYLD4WhoaIhEorp9Tpw4ccmSJbt27bp//354eLj0U1euXPnqq6+cnZ179eplamrKYrHYbLbkrkz03kH6jucACoYgBPgU+vr6hBC508ALCwtNTEyMjIxYLFZeXl69LzQzMzt9+vSIESPGjRsXEhKira1N1zkcjtyxlBUVFZKglUujgoICNpttYWEhN/PIyMhXr14ZGhpKn5anq6sbEhKyY8cOfX19e3v7R48eSb9Ecr89IyMjNze3u3fvfnDsWlpaJiYmdQOSEGJsbDx69OhVq1Y5Ojr269dP+qkNGzb07NkzPDycxWIRQgQCwbZt2yTP0nOruycVQAFwsAzAp7C3t7ezs7t+/bqkUlVVde/evR49eujq6nbv3v3q1avvO+Syb9++N2/e/Oeff8aNGye5G62dnV1qaqpkmuzsbMnXhoSQ6Ojo0tJSycOwsLB27drRdyeWdvToUR0dnYyMjBIpt27dqqiouHjxIiHEy8srNTU1JiZG8pJz585JuoqMjMzIyGjM8H18fOLi4up9au7cuf369Vu+fDkdeBJpaWndunWTFMPDw6WDnz5stXv37o1ZOkDTQhACfNizZ8/OySorK1u6dOnDhw9XrFiRn5//9u3biRMnlpaW/vDDD4SQdevWvXv3LjAwMC4urrq6Ojk5+bfffpOeYe/eve/duxcVFTVs2DD6fL4hQ4YkJiZu27atsLDwyZMn48ePl845sVg8efLk9PT0kpKSn376KTIyctGiRXJN8ni8kJCQESNGmJiYSNd79erl7OxMn1A4ffp0JyenMWPGHDt2LCwsbPLkyZL0Xb58uba29siRI//+++/Kysr8/Px79+5Nnz693pMgBw8e/OLFC3pXn5wvvvgiPDx8zpw5cnUPD4+zZ88+f/68trb29u3b8+fPl2wNE0IeP35sb28vOdcCQKGUfbQOQItGHzVa15MnT8Ri8Zo1a3R0dOiKhYXFqVOnJC88f/689A68Hj16UFJHjUpmbmlp2adPn/LycoFAILlEmY6Ozq5du6SPGp00adLXX39Nb05xOJwVK1ZIDgSVOHHiBCEkNDS07ihWrFihoaHx5s0biqJSU1MDAgJ0dHRMTExmzpx5/Phx8r9dni9fvvTz85P0rKWlNXTo0LoLoiiqoKBAW1v79OnT9EPJUaN1SY4affXqlZubGz1nY2PjM2fO2NjYfP/99xRFiUQiOzu7NWvWNH69ADQhDYqimjxcAdRJvUeFSA7XrKqqevXqlZaWlru7O4cjs9NdLBYnJiZWVVXZ29vXPaqlXvn5+RkZGe3atZM+4tTd3d3b2/v48eO5ubmZmZnOzs5yB7x8jtWrV//222/V1dWSEeXk5GRlZRkaGrZu3VpPT+99L5wxY0Z6evq9e/cavyyhUPjmzRsej9ehQwfJHxCEkKtXr06YMOH169eN/C0BNC0EIUBLJwnCJpnbgwcPPD096aC9f//+6NGj/f39Q0JCPnY+OTk5bm5u165dkzso5hN07949ICBg7dq1nzkfgE+Do0YBmGXbtm03b95s3bp1RUVFfn6+j4/Prl27PmE+tra2r1+/lt7P92lEItHZs2ft7e0/cz4AnwxbhAAt3fnz5y0sLPr3798kc+Pz+dHR0e/evRMKhe3atfP19ZVcEwCAmRCEAADAaDh9AgAAGA1BCAAAjIYgBAAARkMQAgAAoyEIAQCA0RCEAADAaP8fUhEzQedz9bYAAAAASUVORK5CYII=" }, "metadata": {}, "output_type": "display_data" }, { "name": "stdout", "output_type": "stream", "text": [ "\n", "Estimated lockup age:\n", " 102.15909230599004 +0.2004912478051466/-0.2477371468776397 Ma (95% CI)\n", "\n", "Histogram N:\n", " [1, 0, 7, 0, 1, 0, 6, 4, 5, 4, 12, 5, 18, 8, 24, 29, 28, 34, 62, 88, 71, 105, 99, 112, 189, 203, 246, 260, 299, 359, 610, 596, 666, 1017, 1076, 1227, 1519, 1670, 2206, 2581, 3050, 3831, 3984, 4725, 5711, 6476, 7204, 8027, 9168, 10694, 11478, 12561, 14311, 15441, 16226, 16854, 17849, 18560, 18572, 18644, 18475, 18233, 17327, 16414, 14905, 13500, 12023, 10403, 8556, 7302, 5908, 4314, 3460, 2813, 2000, 1616, 1192, 913, 535, 389, 336, 176, 102, 109, 74, 19, 34, 13, 12, 4, 1, 2, 0, 22, 1, 1, 0, 1, 0, 0]\n", "\n", "Histogram bin centers:\n", " 101.44290659563755:0.012622372896352374:102.69252151237644" ] } ], "source": [ "# The fraction of the way through the MELTS dist that lockup occurs\n", "lockup_fraction_mu = 0.2814311218465038 \n", "lockup_fraction_sigma = 0; #0.07475159714446021\n", "\n", "lockup_fraction = lockup_fraction_mu .+ lockup_fraction_sigma .* randn(size(tminDist))\n", "lockup_dist = tmaxDist .* (1 .- lockup_fraction) .+ tminDist .* lockup_fraction\n", "\n", "# Print results\n", "AgeEst = mean(lockup_dist)\n", "AgeEst_sigma = std(lockup_dist)\n", "AgeEst_025CI = percentile(lockup_dist, 2.5)\n", "AgeEst_975CI = percentile(lockup_dist, 97.5)\n", "print(\"\\nEstimated lockup age:\\n $AgeEst +$(AgeEst_975CI-AgeEst)/-$(AgeEst-AgeEst_025CI) $age_unit (95% CI)\\n\\n\")\n", "\n", "# Plot results\n", "h = histogram(lockup_dist,nbins=100,label=\"Posterior distribution\",xlabel=\"Lockup Age ($age_unit)\",ylabel=\"N\",fg_color_legend=:white)\n", "# plot!(h,[wx,wx],collect(ylims()),line=(3),label=\"Weighted mean age\",legend=:topleft)\n", "savefig(h, \"histogram.pdf\")\n", "display(h)\n", "sleep(0.5)\n", "\n", "# Print out histogram\n", "edges = range(minimum(lockup_dist),maximum(lockup_dist),length=101) # Vector of bin edges\n", "hobj = fit(Histogram,lockup_dist,edges,closed=:left) # Fit histogram object\n", "print(\"Histogram N:\\n $(hobj.weights)\\n\\nHistogram bin centers:\\n $(cntr(edges))\")\n" ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAlgAAAGQCAIAAAD9V4nPAAAABmJLR0QA/wD/AP+gvaeTAAAgAElEQVR4nOzdd1wT9/8H8PddEvYWQUAQZbgnbhG1iqLgqLsOwF2r32rd1dqfbW2rtrXubRX3ah111VmtCwTFiRMFlS0ohJXk7n5/HI2IVhESAsnr+fCPu+NyeXMeeeVz4/NhBEEgAAAAQ8XqugAAAABdQhACAIBBQxACAIBBQxACAIBBQxACAIBBQxACAIBBQxACAIBBQxACAIBBQxACAIBBQxACAIBBq/BBmJube+/ePfUsx3E6LMbQCIKALvrKEg7vsoTDu4zp8PCu8EF469atQYMGqWdzcnJ0WIyhUalUCoVC11UYEBzeZUmhUKhUKl1XYUB0eHhX+CAEAAAoDQQhAAAYNAQhAAAYNAQhAAAYNAQhAAAYNAQhAAAYNAQhAAAYNAQhAAAYNKmuCwAAAAP18w1+x0OeiHiBUnONHc0KejCY2oAdUKPs2mkIQig5hmFYFicVAKCEptRnp9RniSg9nzx2Kp700k0kIQih5KRSHD8AUOHh6zwAABg0BCEAABg0BCEAABg0BCEAABg0BCEAABg0BCEAABg03P4OAAAfJlf1loUMQyaSMi9FExCEAADwAeRK6vpXQRKm5pGUIVtjIiJLGR3uUiEzpUIWDQAAumIho3+CCrJj4iXO3YKZWK9iX2Uri+pTU1OfPXtWeEliYuLff/+dnJxceMm5c+ciIyNzc3Pf3MLDhw/PnDmTmZmp9VqLp8sxYtermHVKZr2SXacccUozm1UqlXK5vPCSq1evJiUlaWbrH0ihUJw6paFfDACgHNNuEP7++++urq6Ojo49e/ZUL1yzZk39+vXnzZtXt27drVu3EtHWrVubNGkya9asMWPGuLu7nzt3rvBGpkyZ0rZt27lz53p5eV24cEGrBRcHu155LF4pCAIRkUAC0W+xSvY3vjTbfPjwYbdu3aysrFxdXatXr75p0yZx+cyZM3WVRhkZGV27dtXJWwMAlCXtBmGjRo0OHz68YcMG9ZKsrKwpU6YcPnz46NGju3btmjhxYn5+fu/evRMTE8+cORMVFTV69OivvvpKvf7t27fXrVsXFRV1/Pjx2bNnT548WasFv5dkPScmYBECz5mFve3ycTGkp6e3b9/ew8MjPT09IyNj7dq1EydO3Lx5c6kKBQCA4tFuEHp4eNSvX18ieXUj0V9//eXq6tq8eXMi+uijj8zMzP7++29TU1P1Cvb29oXX3717t7+/v5OTExENHTo0IiLi6dOnWq353Xh6WwwSEVGu6j9/9G5r1661srJavHixuB86der09ddff/3112Kj8+HDhx07dqxevfro0aOzs7OJaPXq1fXr13dxcWnRokVkZCQRXbp0yd/fv1q1at26dYuJiSGiJUuWLFu2bMCAAS4uLtu3bx85cqT67f7888+ZM2cS0ePHj/v161ejRo02bdqcOHFC/Okvv/xSq1atxo0bHzlypEideXl5Q4cOrVmzpoeHx7Bhw54/f05EUVFRfn5+1apVmzBhwldffXXs2DFBEFatWtWsWTNPT89x48aJJ3uXLFlSr149FxeXVq1a3bhxo2Q7CgBAG8r6ZpknT564u7urZ6tVqxYfH09EKSkps2fPfvHixePHj9evX//W9a2trW1tbePj46tWrVp4m9nZ2X/++ac4nZOTY2Zm5u7uXrdu3XeUUfLBg97aHiydK1eudO3atXBJ3bp1++KLL8Q9s379+kOHDrm6ugYHB0+bNm3y5Mlff/11RESEm5vbgwcPjIyMYmNje/bsuXPnztatW+/duzcoKOjWrVt37tzZtWvXjh07Vq5cSUTjxo2bPHly7dq1iWjJkiUDBgzIzc3t2LHjrFmztm7devXq1V69el24cCEqKmrVqlVHjx61sbEZPHjwG7+6MHDgwJUrVwqCMG3atOnTpy9evDgwMPCXX37p16/fjh07Ro4c6enpuXbt2o0bN+7evdvBwWHixImTJk2aMGHCggULIiIiqlSpcu/ePUtLS43vw4qL5z/gpDrP8x+0PpSGuKuxw99NEARB0MBeEjegjb1dnE/7sg7CvLw8IyMj9ayxsbF4d4yJiYmPj09aWlpkZOT58+fr1av37vULe/ny5fLly8VpjuMkEklgYGD16tXfUYaZmVkJsjA87Z0/LmlEJiUlNW7cuPASBwcHInrx4gURDR8+XAz17777ztfXd8qUKfn5+ZcvX3ZwcPDy8iKir776qmfPng0bNnz58mWzZs1MTEyuXLlCRH369OnUqZO4wYEDB27atOnHH3989uzZxYsXd+/efeTIEXt7+48//lgul1euXLlt27aHDh06efLkhAkTPDw8iGjWrFmnT58uXJWpqWmHDh3Onj2bmJhYuXLl7du3nz592tnZWYzM4ODgX375hYhWrlz51VdfWVtbv3jxIjQ0NCgoaOLEiXl5eREREQEBAbVq1SrhbtJHgiDk5OQUf/3c3NzC50tAq/Lz81mWlclkui6kXFOppAqFkJPDlXI7uQpGEKQf9OdQTMX5tC/rIKxSpYp4Sk2UlpYmnva0srIaPXo0Efn6+gYGBo4aNUosvfD6PM+np6eL6xfm7Ox89OhRcTorK0t7DY4W9lrZrLu7e5HzvU+ePCGiKlWqEJGrq6u40M3NLTMz08HBYdWqVYsXLw4NDfX391++fPmTJ08uXLjQv39/cTVnZ2eVSlX4hUQUEhLSp0+fuXPnhoWFffzxxzY2NvHx8XFxcepXEZGZmVliYmLhtytSZ2xsbPv27Tt37lyjRo3c3Nz09PSUlBSxSPVbE1F8fPwvv/yyatUqcaGPj4+np+fChQsXLFgwdOjQrl27Llu2TEx6YBjGwsKi+OsLgvBB60NpyGQyBOF7yWScsTFjYVHaq2yKfGIYha4O77IOwubNm48fPz47O9vc3Dw9PT0mJqZZs2aFV5BIJAzDFF7/m2++EacvX75sYWEhtlf0SZs2bebMmfPTTz+pr5WKt9E6OjoSUVxcnLjw8ePHNjY25ubmAwcOHDhwYEpKyrhx4+bOnevm5iaVSgufTyaiHTt2FN6NLVq0sLKyOnHixObNm5cuXUpEbm5uTk5Ox48fL/yqQ4cOPX78WP12RercuXNnYGCgeK716NGjv/32m7u7+7179wRBYBiG5/m7d++KW54yZUrv3r0LvzY4ODg4ODgpKWnUqFE//fTTTz/9VLp9BgCgMdq9Webx48fz588/ePBgUlLS/Pnzd+/eXbduXT8/v8GDBx84cGDw4ME9evRwd3efP3/+zz///Pvvvy9fvjwkJGT06NEsyw4ZMmT27Nm9evXKz88fP3783r17R48ePX78eGNjY63W/G5OZv/1E6a5Uwm/OYaGhlarVi0oKCg8PPzRo0cLFixYvnz5okWLxJ9u2LDh8uXLycnJX375ZUhISExMzP79+58/fy6VShmGsbGxGTly5N69e7du3ZqRkfH06dP169dnZWW99V0mT56cnZ390UcfEVG3bt1ycnL+7//+LyUlJTk5ed++fTExMcHBwYsWLbpx48azZ8/mzp1bZAt2dnZXr15NSUm5e/fud999R0Tt2rUzMzObNGnS+fPnv/jiixcvXjAMM3HixC+//PLixYuZmZm3b9/etGnT9evXDx06lJ6erq65ZDsKAEAbtBuEKpUqIyPD3d19yJAhGRkZ4g2Eu3btaty48ZYtW/z8/MLCwoioXbt2SUlJu3btunXr1oIFCxYsWEBEXbt2bdOmjZGR0dmzZ42NjXft2vXZZ5/Nnj1bqwW/V8IgmUmhlpaarVQIDyzhNo2MjI4dO9aqVatPP/20S5cu4eHhp0+fbtu2LRE1btz4m2+++fbbb9u3b+/l5TVv3jypVLpp0yY/P7927do5Ozt/9dVX1apVO3v27P79+1u3bt2rV69r165JpdJatWoVuUo6ZMgQJyen6dOni+ecTUxMzpw5k5SU5O/v36lTp71795qZmfXs2XPatGniedRPP/1UjEy1YcOGNW7c2M/P79NPP502bVq7du0kEsnJkycZhlm4cGHTpk0bNWpkZ2cXEhLy7bffzp49u2nTpmPHjs3NzZVKpWvXrvX19e3YsWPNmjWnTJlSwj0FAKAFjKCF2yDLUmRk5Keffio+RUBavkZYmPFGTqniiSETKZMTYrg91cnlcvG0/sOHD5s2bXrjxo0i9/SCBpXZ4Q2Em2WKR1NdrKXnk8dORUaw0ftX1QLD/QQvpfxQCRHu36Pp06dfunTJysrq1q1bP/74I1IQACocBCGUyvLlyxMTE7Oystzd3Qs/6AIAUFEgCKG0nJyc3nymBQCgoqjYY2cAAACUEoJQ/z1//jwhIaEEL4yNjRV7N32vzMzMO3fulOAtyolbt26hJy0Ag4Ug1L3g4OCmTZs2bdrU19d3/PjxJQutdwgLC5s1a1YJXti3b9/w8PDirHnhwoWQkJASvIWuXLt2bd++feK0SqWqV69eCfp2evnypfpxTwCouBCEunfnzp0BAwbs2rVr8eLFycnJgwYN0uz2g4ODxeffQS08PFwcC5OIpFLptWvXzMz+s6+E/5Kenj5nzhwNVwYAZQ5BWC5Urly5Ro0aPj4+48ePv3btGhFFRUWFhob6+vp+/PHHx44dUygUQ4YMSU1NFdfPzMz85JNPsrOzBUFYvXp1YGBgt27dtm3bRkSCICxcuLBjx45t27YdOXJkdnb25cuXDx8+/Pz588mTJ3fo0KFz586//vorx3FpaWnBwcF///13165du3Tpou6vtQie5xctWuTv79+rV68DBw6ICzMyMqZOndqhQ4du3br98ccfhdc/ffr0p59+mpmZOXbsWHX/cNu3b9+yZcvff/+9YMGChQsXtm3bNiQk5M0Rtf7666+BAwf6+voOHDjw6tWr6l+nffv2gwcPPnz48NSpU4koLS1twoQJ7du3Hzp06L17996sedOmTT169OjSpcu6devEJ2VXr17t7+/v6+sbEhISHR29Zs2aS5cu9e/f/3//+x/P8/Pnz1cqlePGjTtx4sTAgQPbtm27efPmFy9efPbZZ23btv3hhx8EQUhJSfniiy/at2/fpUuXZcuW8Tw/a9asnJyc/v379+/fPzk5mef5pUuXdu3aNSgoaM+ePSU8FACgzOGu0XLh9u3bJ06cyM/PX7p06cCBA4no5cuXgwYN8vDwuHv3bkhIyKlTpziOCwsLE7tl2bZt24sXL8zNzWfPnn3p0qUFCxaoVKoRI0ZIJBKpVLpt27YNGzaYm5tfvnxZEISYmJgbN2507dq1cePGo0aNys7OnjBhgkQi6dWr1/bt2zmOmzdv3u3bt/v37//o0aNKlSoVqe3rr78+fvz4ypUrxeCUSqX+/v5+fn4dOnRYtmxZbm5uYmKieuXdu3d//fXXe/futbKyOnTo0NixY8XlN2/eVCqVSqXy22+/nTx58m+//bZlyxZ/f/+bN28WHk4hJydn3LhxLi4u58+f79at2+3btzdu3Lhly5bffvtNLpePGTNGJpPl5+e3adNm+PDh48aNCw8P79Chw61btwp327Zw4cI9e/YsXrxYKpV+9tlnHMd5eXktWrRo27ZtNjY20dHRlSpVCggIuHTp0vTp042NjXme37Zt2+rVqw8fPnzlypVffvklKyurT58+27dvHzt27PDhwwcMGFC3bt369es3b958zJgxWVlZYld/Q4cO/fPPP6dPn05ENjY2kydPvn///rx583JycoYNG2ZiYhIUFKS1QwYANEeo4C5fvuzj46OezczMLOGGFiwQJJKy+LduXZF3btasWZ06dTp16tSuXbsaNWrs2rVLXP7ixYuIiIjjx48HBgYuXrz4zJkzHh4ePM8LgtCkSZO9e/cqlUpTU9P4+Hhx/T/++CMgIGDt2rVt2rR59OiRevu//PJLaGioIAgKheLGjRsnTpz45ptvAgICxOZaWlqauFqjRo1OnTpVuLDGjRufPHnS1tb24sWL/+6kBUFBQfv27atVqxbHceo1jxw50rx580WLFjVq1Ojp06fiQldX12vXronTM2fOnDp16m+//ab+FXie9/DwOHnyZJG9kZqaev78+ePHjzdq1Ojw4cN169Y9cuSI+KMVK1Y0bNhw+/bt7du3V6//8ccfb9mypfAWHBwcrl+/Lk6fPn26ZcuWe/bsadSokdg/uGj16tV9+/YVp5VKJRGJj0Lu27dPXOjv7y+OjSwIwvTp06dPny4IQn5+/vXr148fPz5z5sxevXrFxsZaW1uL68jlciMjI/XO3LRpk3r7GlTywxs+XF5enkKh0HUV5d2Ei6pfb3DvX+99nucJNmH5pd9OyaBF+K8pU+iLL8rijd42ntzUqVNDQ0OJKDY2tnbt2j4+PufPn589e3abNm0sLCyePHmSkpLi5+dnamr6999/W1hYJCYmBgYGPnv2LC8vr3379urteHl5DRo0KCIionHjxm5ubsHBwZMmTRJ/dPfu3cDAwNq1azs7Oz979iwlJYWIzM3N1U1Aa2vrzMzMIoXl5ORkZGTUrFlTnK1du/amTZtiY2Pr1q1bZIivu3fvTp06NTIy0sXF5R2/vaenpzgsBsMwXl5eRca4+O6778LCwnx9fY2NjTMyMlJTU58+faruNLVGjRpE9ODBgytXrhQehKRz587q6czMzJSUlJ49e6oH33B0dOzRo8eZM2datWrl4OAwePBgsQ33VupRqCwtLdXT4g6/ceNGr1696tev7+joGBcXV6Rn88ePH3Mc17x5c/WShg0bvmM/AED5gSD8F8OQVPd7o0aNGiYmJnfv3p0zZ87OnTtbtGhBRJ988on40zFjxqxZs8bS0nLkyJEymczBwYFhmIsXLxYZ3m/NmjUrV648efLkqFGjPD09xYVLliwZMmSIeHPHpk2bFi9eTETM2zoQL8zU1FQcp9DW1paIEhISHB0dnZyc3ry8V7NmzYkTJwYGBh49elQcSdjExCQvL0/86fPnz62srIgoOTlZ/ZLExERxqClRXl7e3LlzExISxGwWG/rOzs5xcXFiEoupWaVKlVatWv3XFU1LS0tzc/MjR46ow1v96y9atOjs2bNjx451cXFhWVb4wF52f/3111GjRs2YMYOIVq9evWnTJoZ51VWvOC7j1atXxV8TACoQ3CxTLqSmpsbGxsbExMyZM4dhmCZNmshkskePHhHRhQsX1LeoDB069NixYzt37hwxYgQRmZqaDhw4cMKECWLrJCkp6cKFC5GRkU+fPpVIJG3btrW3txfP+xGRTCaLi4vjeT45OfnXX38tZmEMw/Tt2/e7775TKBQpKSlLlizp379/QEDA48ePw8LCBEFQKpW3bt0SV/7kk09WrFjRtWvX69evE1H9+vXFe0auX7/++++/i+tcv359//79RLR///6EhIR27dqp34tlWZZlxRO2O3bsEG8aCgkJmT179p07d6KiopYtW0ZEPXr0uHLlyo4dO8QQCg8PL/zACcMwoaGhU6ZMycjIIKK0tLSzZ89eu3bt0aNHLMu2atXK0dFRpVI5Ozvfu3cvMTHx5cuXxdwV4g4UBCEhIUEc09HBwSEvL+/69esZGRm2traBgYETJ04UH8N49uxZREREMbcMALqFINQ9FxeXVatW+fv7Dxo06O7duydOnHB0dFyyZMlXX33l7Oz8448/jhs3TmyQWVtbBwUF+fr6VqtWTXztqlWrKleu3KBBA0dHx44dOz569OjBgwcfffSRk5NTnTp1OnTo0Lt370qVKrm4uEydOvXBgwdOTk4BAQGhoaEuLi5SqdTd3V1dhpOTk7m5eeHCatSoYW5uvmjRInHN5s2b9+nTZ8SIETY2NuLAvA4ODtWrVz969KiZmZl4RrR79+4rVqwYPXp0XFzcDz/8IF5inDFjxqhRo+zs7IjI399/x44djo6OM2fO/OOPPwoPSG1kZLRixYoePXpUrVr12LFjwcHBlpaWU6ZM6dmz5/Dhw7/99ttRo0ZZW1s7ODgcO3Zs7dq1Ytdus2bNUigUhcv+6aef6tSp07RpU0dHxzZt2sTExDx9+jQwMNDJycnb27t27drBwcGdO3f28/MLDAwUT6LWrl2bZVk3Nzf1aJeOjo7qtp2tra29vf2XX34ZHR3t5OTUvXv34cOHOzs7m5mZLV26dOzYsU2bNn369GlYWJhMJqtdu7ajo2PXrl2fPHmi6SMFALQCwzBVJCqVqkGDBj/99FNgYEkHP9SpDRs27N+/X/0ke3HwPK++GPn5558T0ZIlS7RSXEWg34d3eYNhmIoDwzBBmfr999+XL19uZ2fXrVs3XddSdq5fvz5y5MiGDRvGxsYmJycfO3ZM1xUBgL5BEFYYXl5eM2bM8PPze+8dLuVWhw4dvL29P+gljRo12rVrV2xsrJ2dXb169TDSEwBoHIKwwmjQoEGDBg10XUWpuLu7F74qWUw1atQQH5wAANAG3CwDAAA6xgnEk87OdSEIAQBAl1bG8J67lJkKct6m3BengwHREIQAAKAzF5KFcee5TAURUWIODTrNxcvL+lkGBCEAAOjMiYTXnuHLVdG5ZAQhAAAYDIs3btm0LPNHNxGEAACgM73dGatCyVfNgmnnVNbBhCAEAACdcbdkTgdKu7myEoaGerKnAiVWaBFCBaJSqYr08wkA8KGa2DOb20ssZcKm9pIaljp4iAIP1EPJiWNa6roKANCBlwqaGsFtvs8bSSgxR/iuqcSowjasEIQAAPDBxl/gtjzgiSiPowXXeSMJfefzllHHK4QKm+AAAKAjvEBFnnz//VEFPjmEIAQAgNKqqEMBEBGCEAAAPhTLUA+31+Kjl3sFjkJcIwQAgA+2oo3EVEpbH/AylsbUYv+vSUW9QEgIQgAAKAFrI1rXVmIhI42MUK9bFbt6AACAUkKLEADAIJxJFKLSCu7tVPCkfuyvvRPTxL4CX+ErPbQIAQAMgqWMqphRFTOyMqJvr3LidBUzMi/zLs3KG7QIAQAMQhP7gpbfSwVNvsQN8kBDqAB2BAAAGDQEIQAAGDQEIQAAGDQEIQAAGDQEIQAAGDTcNQoAUH7lqMjvoEqczlISy5C5lIjIXEpngvABrhnYjwAA5ZeZlCJ7FXxQTw7nnM2YyfVxJk/DsEMBAMCgIQgBAMCgIQgBAMCgIQgBAMCgIQgBAMCgIQgBAMCg4fEJAADN23Sf3xcnEJFAlJpLDqYFy0O8mJ7V0AIpXxCEAACa16Ma295JIKIcFTXZq7rTr+DD1tbYoIfALZ8QhAAAmmdjRDZGDBFlq4hhyM0C+Vd+IQgBAAooeZoWwYnT8nyGZclMxhGREUvzm0t0WhpoEYIQACq8XBXlcm9ZbiYlkw/JLwlDvd0LLuAtv8VZGxXMStCc02sIQgDQDU6gEWcL4itbRQyRmZSISMrSurYf1vzaHsvveMgTkSDQzQyhvl1BcI2oyQ6o8QF3prAMta1S8Np9j8je5NUs6DEEIQDoBsvQlAYFKbX4Jm8mpVG1WHH5hxruzQ73ZokojyObTcpjXfHJBh8AhwsA6AZDVM+2IPQqm5CFjFHPApQlPM4CAAAGDUEIAAAGDUEIAAAGDUEIAAAGDUEIAAAGDUEIAAAGTetBuHLlypCQEH9//8uXL4tLBEH48ccfvb2969Sps3z5ciLKzs6eNm2aj49PjRo1unfvHh0dXXgLvXv39v/X999/r+2CAQDAoGj9OcLo6OgmTZocPXo0LS1NXLJ9+/Z169YdOnQoNze3W7duXl5e9erVUygUK1ascHZ2XrNmTefOnWNjYy0sLMT1z5w5s379eldXVyKytbXVdsEAAFA24uRCUg4RUaaSOJ7CUwRxubsl42j6rhdqltaDcPXq1US0ePHiwksmTZpUq1YtIho7duyaNWv27NmzaNEi8afffvvtzz//HBMT06xZM/VL6tWr5+npqe1SAQCgLF1OFU4nFoza2MWZ3/SAF5f3q846mpZd7wo66Fnm1q1bTZo0EaebNGmyY8eOwj+NjIxkWbZI7A0ZMkQikbRo0WLmzJn29vZlVysAAGhN3+ps3+oF01lZOZaWJjopo6yDUBCEjIwMKysrcdba2jo1NVX904yMjODg4O+//77wKdB58+Y1atQoLy9v3rx5AQEBly5dkkpfKzs2NtbGxka9fYZhhg8f/u2332r/tzF0SqWS53mlUqnrQgxFdnY2w+hnJ2QKhUQhMHK5qpTbyeOIyFgul2uiJFKwjFyeX8rt5HAMkZFGSlIqpfn5glz+toE2PkS2khEEWbkqibR2eJuZmbHse+6GKesgZBjG2to6KytLnM3MzLSzsxOns7KyunbtGhAQMHHixMIvGTVqlDjh4+NTqVKlGzduNG7cuPAK7u7uJ0+eFKflcrmFhYWJiYmpaRmeYDZUYhAaGxvruhBDIQiC+tq5njEy4oxkjIVFaRsEUo6IlBrZS0ZGCiMjxsJCVsrtMCqNlSSTccbGjIVFaW9y5BTEMOWrJNLp4a2DU6MeHh537txp2bIlEd25c8fDw4OIsrOzu3Xr1qRJk4ULF/7XC01NTY2NjXNzc4ssZ1lW3YKUSqWWlpZaqx0AwNDJlRR0rKDtnpxLMlbYF8cTkYWUDnapkAM5aL3orKwslUrF87xcLhdPioaEhCxdurRv3775+flr1qz55ptvcnNzg4KCnJyc5s6d++LFCyKysLA4f/58RERE//79X7x40aBBA4VCMXfuXHNz84YNG2q7ZgAA+C8WMjrY+S3ZUXFP22s9CENDQ6OjoyUSyYwZM4jowIEDY8aMiY6OdnFxIaIRI0YMGDDgzp078fHx8fHx6jtF16xZEx8ff/bs2fbt2w8cOPDZs2dSqbRp06YHDhwwNzfXds0AAPrqcqrwVSQnV1G/k9zcpmxN65LEV6lPGJcvWg/C33///c2F69atW716NcMw4jXM2rVrP3z48M3Vhg0bRkSxsbEqlYpl2fde8AQAgHd4lCV0PKzKUhIR7XnEX0wRbjWmTXwAACAASURBVPaR2hjpuixd09n5XIlEUvyVi9wmCgA6dDFFUD/4rORJ9u8X1NaOTPPKFfbsmGHYFydkFbrL+1m2cCqB7+1u6G0MQ//9AeBDmUnJ3oTsTcjWmGZc5sRpexMyw/fVci/3jedTckr7xIo+wJELAB+moR3T0I4hIhVPI85yQzzxfbrC8HdhZkcRX9CeJzMp+VVBIx5BCGAwbmUIby5kGKpjg49CQ9GsMrO2rWRaBPc8j1zNmeVtWDcL/O8jCAEMg0C04HpBR473MwWWyMOKISKWoQ1+H3DBHiq64d5sH3fWdZsy/hN8/hfAjgAwCAxRWLuCwJsdxRmxzOzGOKVZYeRx9GM0t+k+b8wyRDSxHispXUOu4j7zpw0IQgCA8m5aBLf0ltigF6aEc0qeZjTE9xiNwa4EACjXBKKt/45PJNry+iyUEoIQAKBc4wVxVI1X3nwKAkoDQQgAUK5JGOpS9bXP6q6uuMSnSQhCAIDyblUbSZAbwxBJGBrkwc5rhht9NQk3ywAAlHcOpvRnZ+nES5yzGTOtgf40YLJO7cm9eoaISBBUL5/n2tiLyy0+6mvWuF2ZlYEgBADQlhwVrYrhFTytvcOHeLNGpYswCUOlfGqivLH8qK/lR32JiM/JSvwu1GHyUp2UoT/fLAAAypVcFbU8oJoSzql4Gn2O63RYpcLNnuUSghAAQCt2P+JvpL/q1u6fJOFEwlt6uQOdQxACAGhFnPzNJQjC8gjXCAEAXqPgad41fvNDRsaSkZQvcX9mjSq9scROvy7x6QsEIUC5ditDUI8akZvLmpoWXGWqZ8dg1Ig3RaQKnEBX0oQm9iXfOTMiuF9vFuznKeFcPkczG5Xk5FmQGzvUU9j8gCcihmhiPbaFA/7LyiOcGgUo114o6LG84N+nEUbq6Rf5uq6snBGIBp3m2h1UqXjy2acae557/2v+Q5EOzDaXtD8zhmhTe8m57lJjlm70kS5siYf/yim0CAG0IiFHEN52PcjF/MPaBG0cmTaODBEJRNMjOH16hkyzDsUL2x++SqxVMfxgD9b3w0ed5QXKfT1DSzmGe6NKjISlurZoC5ZfCEIArfjsfMEw4M+yBYGo6r/5t98fzQKtuJ5e9HvHtXShBEHIMtSlKvv7o1eZ2rUqMkzPIQgBtGLfv4H3zRWeF4RvfJB/b5etohkR3IZ7vIRhnucL3zeVmJRoV1W3LLqkhmUJA2xFa0k+Jxx6IjBEfauzC1rg/07PIQgBQJcmX+JW3ykYaW/hDYEh+rlEwdPbnfWx56PSCtqF7ZyYziVtyYn9mU04r7A3YWbjG4wBwPUGANCl3x/zr8+W8Ek7Ywmd6y5d3loiYWljO8nxrtJS9kYmZUmKD0jDgP9nAChH3nqHUTGZSGh4TVbKUIgXK8NnGxQbDhYA0KWe1V77FOrljjtToKzhGiEA6NKilhIJQ5vv8wzR8Jrsj01xTQ7KGoIQAHTJQkarfSWVjMlCxpSsAxeAUkIQArzCCaR8Wy8iUgb3TQDoLQQhwCub7/Pr7hYk4f2Xgpd1wfWqYd7siJpIQgD9hCAEeCXUmw31Lgg8dp0ycbCUxa0bAPoOX3IBAMCgIQgBAMCgFffU6I0bN06dOhUdHZ2WliaVSqtUqdK8efNOnTq5urpqtT4AAACtek8QchwXFha2bNmyq1evMgxTuXJlOzs7lUp1+vTpVatWsSzr7+8/adKkzp07l025AAAAmvWuU6O3b99u1KjRpEmTGjdufPDgwbS0tOTk5JiYmPv377948SIuLm7z5s3Gxsbdu3cPCAh4+fJlmRUNAACgKe8KwtjY2E8++eTp06fr168PDAy0s7Mr/FM3N7dBgwbt37//4cOHNWrUSE9P13KpAAAAmveuU6NBQUFBQUHv3UTVqlVXrFihuZIAAADKzgc8R3j27NmdO3c+fPgwLS2t8PLIyEhNVwUAAFBGihuEYWFhw4YN8/b2zsnJMTc3t7GxuXbtmkwm69Kli1brAwAA0KriPkc4Z86cIUOG3L59u2PHjn369Ll48eL9+/dr167t7e2t1foAoHza8ZB326HkBKq6XfX7o7f10ArlTHo+xWYJsVnCY7nAU8F0bJbwUqHrynStWC1CuVweFxe3c+dOlmUZhsnPzyciFxeXVatWtWjRYurUqdbW1lquEwDKkVsZQvAZTuyg/Fm2MPhvroEdo+6aFcqng/H8wXiBiASipvbsjIiCry8h3mygq0H/3xUrCPPy8gRBsLKyIqJKlSqprxHWrFlToVA8fPiwSZMmWqwRAMqZE8+EwsN05HN0KlFAEJZzwV5ssJeuiyiXinVqtFKlSpaWlo8fPyaimjVrHj9+XHxq8OjRo0RUuXJlbVYIAOWO6Rtfoc3QgT9UWMU6eBmG8ff337NnT0BAwCeffDJr1qxatWp5eHhERERUvF7WIiJIYfBnxDWEUalYQSCZTNeFaIXvAxWdk1KpGzlu93heIMrVQL++BSWVWrW7vJQlyi55SR/nCfviuCxlwayNEQXFSSmp5CW5x/CmUiJ5afcSy1Pr+5rZS9VvKq2NGcou7aYknOZKus3bmxC91Mc+ovNyjJ4+p3PnirWylxc5OmrwzRlBEIqzXm5urlKpFM+O3r17d9myZY8fP27UqNG0adMsLS01WNCHioyM/PTTT9WPcGRlZb2nnp496fnzsqjMAAiCIAgCy+rjnyXR+WShjaMGzvU9kZNAgpuFBjalqZIeywWGmGoWpdpItori5EJGPtkZUzULppQtwji5IGGYqual2ggR8QKFpwqtHDSwlx5lCTKWqpqXdlOcQBGaK8lIwriYlX5L5Q/PKZ7cN6pWq1grz5hBxXjGvfjeE4Tnzp3Ly8t79yY6deqkwYI+1AcHIWiOUqnked7Y2FjXhWgFu06pGiEr/XiE31zheUH4xkdSyu0IROw6pTCyVO3v5FwafFp1MkEgos4uzNYOUnuTkm9NxZPJBqVqhAZOCcy8zFnImJmNSvulKo8jm03KvGEaKGnyRYW9CfNl49JuKltFDluU2aGaKCmcczZjJtfXw6+efE5W4nehLj/+rpN3f8+3uCFDhsTFxb17nWK2KQFA58Zf4MQUJKJjz4QJF7mtHUqb0KBVuSoKOqYSp5NyBRnDHH7CE5GZlP7sjAuzmvH+/WhhYdG3b9/evXubmpqWQUEAoD0nnr32wN+JBJ4IQViumUrpj05v+aBmcIuu5rwnCDds2LB+/fqdO3fu3bu3f//+oaGhrVu3LpvKAEDjrIyYF4pXp3CsjfBpWgFYG+m6An33nnPNHTp02LJlS1JS0pIlSx49euTr61uzZs05c+aIj1IAQMUyzPu15Av10sOrTQAfqlh/BlZWVsHBwcePH4+Jiendu/e6des8PT0nTZqk7eIAQLNmN5YsbiVxtaBqFsyy1pIZDRGEAB8y+gQR1axZc8CAAdnZ2StXroyJidFSTQCgJRKGPq/LpuYJRiwzrg5SEICo+EGYkZGxe/fuNWvWREVF1alTZ+7cucOGDdNqZQAAFde1dCEmQyCifI44nnY8LLhNqVElppYNLs2WL+8JQo7jTp8+vWbNmv3795ubm/fr12/16tU+Pj5lUxwAQAX1UkFPcwqmZzeRqKc9rHRVEfyn9wRh48aN79y5061bt507dwYGBsr0tDMtgHJu72P+i0s8EVXfqfq1JdurGs5qlnd+VRi/Kmj5VQzvCcLMzEyWZc+ePXv27Nn/Wic9PV3TVQHAK3dfCp+c5vI5IqLHWcLAU9z13ow3hnoA0JD3BOGgQYOeo2dOAJ068UwQU1CUz9HJBAFBqCV5/+5qlUAqoWCWITJGxwP66z1B+MMPP5RNHQDwX978CDbBh7J25HPU6XBBf2ZZCmIZ+uuZiohMJHSiG/oz01vv+q89depUQkLCJ598IpG8689OLpcvXbp0yJAhFWw8JoAKomtVxtaYMvILZu2MqUtVNAdfc+SJcCqBJyJOIEspTQ0vaNZ1c2M7OH3AvjKW0LnuBZ+K+fn5LMvIZMg//feu/2NLS8upU6fOmTMnODi4d+/ederUKTzgjkKhuHz58s6dO7dt2+bs7Dx69GjtVwtQkXACLb7Jr73LC4JgacRMrMtKS3SPi4s5c7yr9KtI7uhToasrM9dH4myGIHxNTRsy+Xfndi90J5EHhqKBYnhXEDZr1uzevXu//PLLypUr/+///s/KysrLy8vW1pbjuOfPn9+5c0ehUNSpU2f+/PmhoaHvbjUCGKDvrnLfXCl4emxqOJepEL4t6WBMPvbM4QApu055uAsaKG9Rw5KpgcyDknrPF1RLS8s5c+bExcX9+eefI0aMsLGxSU5OzszMrFq16pdffnnu3Llbt26NGDECKQjwpo33XhuhLOw+BiwDKI+K9e3SyMgoKCgoSKMjAgPovRzVa8mXrUQQApRHWj/NEh8fHxkZmZaWNnz4cKm04O1OnTp16tSpKlWqhIaGWlhYcBz3999/h4eHcxzn5+fXrl27wlvIzs7euHFjQkJChw4dOnXqpO2CoSJqeUCl4omIspQkEFn92/FDeE+pRHdX0zpXZbc+4AvP6qwUAPhv2v3LjIqKatCgwfz588eMGaNQKMSF69evHzJkiLW19YkTJzp06MBx3ObNmydPnpyZmalQKPr37//dd9+pt8DzfMeOHY8ePWprazts2LDVq1drtWCooC71kEb2kkb2kg71ZPtXZ8TpyF66TEEiWtxK0t2t4E8syI1Z0gpXEIq691IQ/6XnC2l5BdP3X6LpDGVKuy3CRo0aZWRkxMfHu7u7i0t4nv/hhx9Wr17dvXv3SZMm1apV68iRI3379g0NDRVXaN68+bBhw2bPni3OHj9+PCkp6Z9//pHJZA0aNBg1atTIkSNxSRIqhErGdKCzZNZlhheEH5vjoC2KE2jOvzcTyZUMyxTMShna1B67C8qOdoPwzcSKi4t7/Phx586dxZ926tTp9OnTha8+CoJgamqqnj116lTHjh3FPk4/+uijpKSkBw8e1KxZU6tlA2iQkYR4tHDeRsLQtg4IPNC9sr4VOyEhwdra2tjYWJx1dHS8e/eu+qc5OTkzZ86cPHmyekliYqKLi4s4LZVKK1WqlJiYWCQIU1NTP/vsM3FaoVAYGRn5+fn16tVLu78JECmVSp7nBaEcfcyrVKyKE/LylJrYmCQvL48t9clVlYrlBQ2UJPxbUmkLIlKpGJZl8vL496/67u3wGiupHMrPz2dZluO4968KpSbk5xORNo4lIyOjwk/Av9UHBKFSqTx69OitW7cyMzPFrtfu3LljZmbm5uZW/I1IpVKef/Xnx3GcekQLhULRr1+/Bg0afP755/+1vkqlUt9xoyaTyby9vcXp/Px8Y2NjBwcHnD4tAzzPMwxTrnY1wxDLMqWs6Ek2zbtOAtHnlyUz6lM1i1KXxJS2JCoIwrecZSkBliWG0cCmBEZjJZVDEomEZVl9/e3KG55lSTvHEsO8/8tscYPw2bNnAQEBN2/eNDExqVSpkhiE69atu3DhwoULF4pfk7Ozc2Zmplwut7CwIKKEhAQnJyciUiqVAwcONDExCQsLK5zeTk5O8fHx4nRubm5GRoazs3ORbdrY2EycOFGczsrKsrTEg7Vlh+f5cjU4l0TCsyTIZCX/c8rIp/ZHVE+yBSJaf4+OPmOu9ZZWMi5VSbxQqpJEAhGRUiN7m2U5CcvIZKW9V47hNVZSOcTzPMuy+vrblTe8TEZEutrbxf1LGDVqVH5+flRU1JEjR9QL+/fvHx4e/vLly+K/n6ura/369f/44w8iysnJOXr0aGBgIMdxISEhubm527ZtU++Ie/fuxcfHBwUFHTt2TC6XE9GBAwe8vLyqV69e/LcD+FCHnvBiCoqeZQsH40t7ChEAyrNitQjlcvmxY8cOHDjQpEmTM2fOqJd7enryPP/kyRNra+u3vlChUAQGBoqnfYOCgiwtLffv3z937txhw4ZduHAhMjLSx8fHz89v7dq127dvb9u2rfqumf3798+YMcPDw+Onn35q06aNn59fixYt9uzZs3bt2uK0cwFK7IWi6BJ1b9cAoJeKFYSZmZkcx73ZFFOpVESkfkDwLVuXSqdPn154loi6d+9+6dKlf/755+OPPxYfkA8ICDh+/HjhFxoZGc2dO9fExISI9uzZc+LEiWfPnk2aNMnLy6u4vxlAibStwrDMq/s8WYYwzjiAfitWEFauXNnc3PzSpUu1a9cuvPzYsWNSqdTT0/O/Xsiy7Fv7gvH09Cz8KldX1zeHcKpTp456I+LjFgBloKEds6y1ZGo4l60iMynNby5pYo8gBNBnxQpCmUw2ZMiQGTNmODs7i08+cBx3+PDhyZMn9+vXz8rKSstFApSpsbXZUC/WfKMyebDMArdKAOi74t41+vPPP9+/fz8gIMDY2JjjOBsbG7lc3qRJk6VLl2q1PgCdMJUSEZlhyCMAA1DcP3QLC4vjx48fOnTo2LFjycnJ1tbW7dq1GzBgAO4tBgCACq24Qbh58+bc3Fwiql+/fv369YkoJydnw4YNtra27u7ujRs3fvM5dwAAgHdQJjxSpSUSkZCfSxyXe73gqXSZS3VpJacyK6O46TV16tTk5OT/+mn16tV37drVtGlTDVUFAAD6j8tMVybHEREJglGrbgXTRBJbeyqHQbh169aQkJD//e9/3bt3r1y5cmJi4vbt23fs2LF169aMjIzp06d//PHHDx48UHciCgAA8G4mtXxMavmI0zrsF6y4Qfj555/PmTNn5MiR4mzlypUbNGjAsuwPP/xw8OBBLy+vWrVqnT9//qOPPtJaqQAAAJpXrC7WUlNTb9++7efnV2R5u3btxI5mvL29q1at+vTpU80XCAAAoE3FCkKZTMYwzJUrV4osj4qKUp8L5Xne3Nxcw9UBAABoWbFOjdrY2HTu3Hn8+PFKpTIoKMjW1jYlJWXnzp3ffvutOLJ8fHx8YmKieiwkANBjR54Ix57xRCQQ2RjTF5cKRuzr5sr6u6AXHqh4inuNMCwsrFevXsHBwUQkkUjEwSqDgoJ+/vlnIkpPT58/f369evW0VygAlBN1bMni3yGc+ri/OqvkjgHQoGIqbhA6OjqeP3/+n3/+uX79elJSkqurq4+PT7NmzcSfNmrUqFGjRlorEgDKkWoWTCkHKwYoVz7gKXiWZdu1a9euXTv1ktzc3D/++GPw4MFaKAwAAKAslHCI6qioqAkTJri6ug4ZMkSzBQEAAJSlD+sXLTk5ecuWLRs2bLh165apqWnPnj0HDRqkpcoAAADKQLGCUKVSHT58+Lfffjt8+LBKpXJxcTEzM0tKStJVLwAAAACa8p5To3fu3JkxY4arq2vPnj2vXr06adKk+/fvf/PNNyzLIgUBAEAPvKdFGBAQkJKSMmTIkJCQkNatWzMMQ0RibzIAAAB64D1BaGlpGRcXFxkZWadOHS8vLwcHh7IpCwzBhnv88tu8OP0sR3A2Y8SHsYd7s5/VKeFtXPBfeIGa71eJ03IlEQn743gikjAU3hNjqIFBe88fQFRU1F9//bV58+apU6dOmTKlQ4cOQ4cOzcnJKZviQL8N82aHeRcEntFvyriBUiPEn9awDEX2QuABvMV7/jCMjIy6d+/evXv3Z8+ehYWFbdy4MSQkhGVZmUx29uxZX19flsVHFwAAVGDFjTEXF5eZM2fevXv37NmzwcHBMpmsXbt21atXnzFjhlbrAwAA0KoPa88xDNO2bdsNGzakpKTs2rXL29t7wYIFWqoMAIhoXxw/7Cw37Cw3/CxnbyKI08POcgfieF2XBqAnSnjNwNTUtF+/fv369bt7965mC4Ly73GWIH4Gq1TE82SkEMTl1S0ZDD2gca0dWA/Lgj08pnquerAzR1PsbADNKO3F85o1a2qkDqhAvoriFRwRUXwWqQS2hlVB02RrB4lM15eMczjKU+m4Bs1yMCWHfzMvSyZYWiL/ADQMd5HBB9vSXiJOfH+Ff6EQfmop0W09otQ8GnRadeKZQERXnwvbOkidzHRdEwBUBLr+Ag+gIRMucmIKEtHficK4C5xu6wGAigJBCHriZMJrN4+cSsC9JABQLAhC0BPWRsw7ZgEA/guCEPTECO/XDubh3ji2AaBYcLMM6ImpDVgbY5p/jeMEmtpAMra2joNw0U1evGJ5OU0QBPr1ZsGp2i/qIaEByhcEIegJlqExtdjUXMrjhPHloM9uB1MSBCKiT2rgJC1AuYYgBNCKQR66D2MAKA78rQIAgEFDEAIAgEFDEAIAgEFDEAIAgEFDEAIAgEFDEAIAgEHD4xMAr1xPF6KfC+rZzfd5cYjFhpWYhnZ4HBBAPyEIAV7J4+iFomD65xaSl8p/l+vXGIcAUBiCEOCV5pWZ5pXR8gMwLLhGCAAABg1BCAAABg2nRg3F/ZfCrRfCm8u9rZk6NjgZCACGC0FoKNLy6WZ6wfTPN/iJdVkpS0RkJaM6NjqsCwBAxxCEhqKVA9PKoaDl9300N70ha4r/fAAAXCMEAAADhyAEAACDhiAEAACDhiAEAACDhvslyjslT2956IHICN9hAAA0AUFY3gUcVeWqiIheKIjjqZJJwfJTgVITiQ7rAgDQEwjC8u5kt4L/o5+u8ym5wk8t9DD9ziQKPNGpBKFLVQbP9gNAGcP5NdAlgWjwaa79IRXHU9ejqo+Pc/xbTwQDAGgNghB06XSCsO0hr57dH8cfjOffsT4AgMYhCEGXbmUUbQDeeqGTQgDAcOEaoVZsfcD/erOgZfM0W3AxY8RrX0M82Yn18OXjFS/rotcEvax0UggAGC4EoVYM9mQHexYEnmWY8u8gqaVMtxWVU51dmICqzNGnBe1CvypMr2r4ogAAZQofOqBLLEOHukh3d5RIGNrcXnKym1SKQxIAyhY+daDkslWUrSrt8w4sQ32rsyxD/WuwSEEAKHs4NQolkZZHA0+pTiYQEV3PUO34SFLVHE8AAkCFhG/gUBJfXOJOJhRc2DufLIw9z+m2HgCAEtNNEAqCcP/+/cTERJ28O5TeqYTXHns4nSDgOXgAqKC0G4T5+fldunSxt7dnGObKlSviwpSUFB8fn+7duzdu3HjYsGE8zxPRZ5995ubmxjDM2rVri2ykUqVKNjY2dnZ2dnZ2I0aM0GrBUEw2Rq/PGqNrNACoqLQbhCzLDh8+PDw83NTUVL3w+++/r1mz5p07d+7du/fPP/8cPHiQiAICAg4ePNiyZcu3bicyMjI9PT09PX39+vVaLRiKaWSt146cEd46zsHo58LV58LV50JijpCUS1f/ndVtVQBQIWj3ZhmZTDZgwIAiC7dv375r1y4isrKyGjRo0Pbt23v06NGjRw8ikkrfXk9+fr5cLrewsNBqtVB8E+ux1kb0QzSn5IVpDaSf1tbxxeZVMTwnEBHlqEggWnG7oDeD1b4SFm1VAHinsr5rNCcnJzU11dPTU5z18PA4efLke1/Vvn373NzcatWqrV692tfXt8hP8/PzIyIi1Ns3MzOrUqWKm5ubZiuHwhii4d5sopx7oaBxdXR/y9UqXz0clAMAykZZB2FWVhYRmZgUjKpnZmaWmZn57pdERER4eHjwPL9gwYLevXs/fPjQ0tKy8AopKSmfffaZOM1xnEQi6d279+eff16C8q6kM5zwlhZEEzteUtKGhSAYZWdnM9LSnqZTKKRKJcnluaXcDhERGcvlcq7U//kKBaNQCHK5QlMlYbThd8vOzmZwNbas5Ofnsywrk6FTqDKipcPbzMyMZd/zyVLWQWhvby+VSjMyMuzt7YkoPT29SpUq736Jh4cHEbEsO2PGjPnz51+7dq1Io9DV1TUyMlKczsrKKhKTH+TPW1wuR0R096WQo6LGlQr+V1pVlZR4FFyGUZqbm1uU+q/JyIiXcYKFhcn7V30/pYWFhWmp//ONjJRGJFhYGL1/1eKVhCB8N0EQcIGgzMhkMgRhWdLh4V3WQSiRSBo0aHDx4kUvLy8iunjxYpMmTYr52qysrNzcXCsrLfbKrB72duEN/mm2sLAlTrgBAOg5rQfh9u3bs7KyVCrVH3/8ERkZOWTIkM8///zrr792cXF58uTJ3r17xccqjh8//ujRo6SkpH/++UcQhC5dupw5c2bv3r0zZ848ffq0j49PXl7ezz//3LRp07p162q7ZgAAMBxaD8IbN248f/582LBhqampqamp/fv3DwkJUSgU8+bNs7S0PHjwoHjmMzY29sqVKx999BERRUVFtWjRwsXFpWnTpvb29o8ePTp27JiJiUmnTp0+//xziQStNAAA0BitB+EPP/zw5sJRo0aNGjWq8JIxY8a8uVrHjh2JaOXKlVqqDQAAADcnAACAQUMQAgCAQUMQAgCAQUMQAgCAQUMQAgCAQUMQAgCAQUMQGhZeoD2PeJVAO2N5Fa/ragAAygEEoQHhBep+TNXvJKfiadhZrv0hlRJZCAAGD0FoQE4mCIefvBoE43yy8MdjJCEAGDoEoQG5/7LoUFD3X+qkEACAcgRBaEDq2hYd66uurU4KAQAoRxCEBqSdEzPU89X/eHc3tmc1HAAAYOjwOWhYNrWXnA2Sylg61EW6v7OExWjnAGDwEIQGp20VRsJQBycGIQgAQAhCAAAwcAhCAAAwaFofmBdKjxdo+W1+2W1ewQl2Jszk+qwRvsAAAGgIgrAC+PkGPz2CE6dnXuae5wk/t5DotiQAAL2BlkUFsPHea/2/hN1HdzAAABqDIKwAclSvzeaqiC/aRQwAAJQQTo1WAJ2rMmvvvIq+Ti6sfjz/F3afX3G7oHVra0S+f6rEXyvUmx1bG1/RAKCMIAgrgJ+aS57n0d44XhCoS1Vmla+eXCAM9mKHeL4l8PQj5gGgokAQVgDWRvR7J8kP0UxyrrC4lZ6kIBExRBJkHgDoGk5AVRgylvDUBACAxuGTVYvE4eAVPO2M5TEELgBA+YQg1BaBqOdxVb+T3qI0RgAADMBJREFUnIKjUf9wbQ+qFMhCAIDyB0GoLX8nCgfjX93qGZ4i7I5FEgIAlDsIQm15y3DwmXj6DwCg3MFdo9pS743h4N9cUkHxAomRzhMJRNy/+Y5bQAGgIkIQaktrR2a4N/vbv72jBboyH7vrSfu7/SFVPkdElK0kQWD+SS7o+eaf7lLc1woAFQ6CUIvW+0mG12Q7Hlbt7igNctOfgXDPBhUcNkqlkud5Y2Nj3dYDAFAa+AKvXW0cGRlL7TEcPABAeYUgBAAAg4YgBAAAg4ZrhEXxAq2I4Rff5PN4wc6YmdYQw8EDAOgzBGFRC2/yU8MLhoOfHcWl5QuLWupPP9cAAFAEGjtFhb0+HPym+xgEFwBAnyEIiyoyHHweh+HgAQD0GYKwqC5VX3vSoZMziw5TAAD0GIKwqPnNJf2qs+J+6ezCrNaX4eABAOCtcLNMUZYy2tVRMv8a8yRbWNYaKQgAoOfQInw7DAcPAGAg8GEPAAAGDUEIAAAGDUEIAAAGDUEIAAAGDUEIAAAGDY9PGIqTCcL+uILe46yMaEoEJ3YU0MmZ6VEN34cAwHAhCA1FbRsylRQE3sAar5LPxVxHBQEAlA8IQkPhbMY4m+m6CACA8gfnxAAAwKAhCAEAwKAhCAEAwKAhCAEAwKAhCAEAwKAhCAEAwKAhCAEAwKAhCAEAwKAhCAEAwKChZ5nybs4VTsEREUU/F3I5mnmZE5d/4yOR4WsMAECpIQjLO78qrEogImrv/NpyltFJOQAA+gZBWN595IzEAwDQIpxcAwAAg4YgBAAAg6bdIExNTV2xYsXIkSNDQkLUC1++fDl+/PimTZv27dv33r17RKRQKMLCwv73v//1798/MTGxyEaio6N79OjRvHnzqVOn5ubmarVgAAAwNNoNwvv37587d87Y2Hjfvn3qhWPHjk1ISNiwYUP9+vUDAgJUKpVcLt+7d6+9vf3u3bvlcnnhLWRnZ3fp0qVDhw7r1q2Ljo6eNm2aVgsGAACDI2hfeHi4lZWVOJ2UlGRkZPT06VNx1sPDY9++feK0SqUionv37hV+7fr165s1ayZO37hxw8LCQi6XF17h8uXLPj4+6tnMzEyN1PzLde6LiyqNbMpioyJToZEtlTsKhSIvL0/XVRgQTR3eUBx5eXkKhZ7+6ZZLOjy8y/oaYUxMjIODg4uLizjbvHnza9euvWP969evN2/eXJyuV6+eIAgPHz7UepUAAGAwyvrxiZSUFBsbG/Wsra1tSkrKO9ZPTU2tXr164fWTk5OLrBMXF+fp6SlOC4LAMMygQYNKeRI1P1+iVDJyuQYuSQqCUXZ2NiMVSr+p8kapVPI8r1QqdV2IocjOzmYYPE5TRvLz81mWlclkui7EUGjp8DYzM2PZ9zT5yjoIra2tc3Jy1LNyudzV1fUd61tZWRVePysrq3COipydnXfv3i1OZ2dnm5ub29raWlhYlKZOY2NephIsLExKsxERwyjNzc0t9PGvSQxCY2NjXRdiKARBKOWBDcUnk8kQhGVJh4d3WQdhtWrVEhISxLgiogcPHvj6+r5jfXd393PnzonTKSkpcrm8WrVqRdaRyWTe3t7idFZWlqWlpRYKBwAA/VTW1whr1apVt27dtWvXEtHly5evXbvWu3fvt665c+fO48ePDxw48PTp0zExMUS0bNmyjh07Ojg4lGnFJbL1Ad90n0r8Z2PEtD9UML3oJq/r0gAA4DXabRHGxcW5u7uL0wzDeHl53bt3b82aNX369Fm2bFlGRsaqVasqVapERJUrV05LSyMisW338uXLvXv31qhRw9/ff/78+a1bt7a1tS3yGEZ5NtiTHeyJzgoAACoA7QZhtWrVBKHoTSJNmjR5+PBhQkKCg4ODkZGRuDA1NbXIajt27BAnxo0bN2LEiLS0NBcXF9wpAAAAmqWbTrdZlq1atWrx1zcxMfmg9QEAAIoJp+8AAMCgIQgBAMCgIQgBAMCgYWDe1/zxmFfyRETRz4Xn+cLO2IKnHfq4s1J8ZwAA0EcIwtc8zaZ8joioji0jCEz8vyNh4Ok/AAB9hSB8zed10e4DADAs+NwHAACDhiAEAACDhiAEAACDpm9BeOjQIQyPV2bu3bt38+ZNXVdhQCpKX7v64ebNm/fu3dN1FYZCqVQeOnRIV++ub0E4ZsyY7OxsXVdhKPbt26fuEha0TRCE4OBgXVdhQLZt23bgwAFdV2EosrKyxo4dq6t317cgBAAA+CAIQgAAMGgV/jlCc3NzlUoVFBQkzkql0oEDB0qlFf73qhCePHmiUqnEYZOhDFhbW6sPddC2R48eyWSy8+fP67oQg6BSqSQSiTYO76VLl1avXv3d6zBvjhdY4Vy+fDk5OVnXVQAAQLnj5+dnZWX17nX0IQgBAABKDNcIAQDAoCEIAQDAoOnPTSUqlWrt2rXR0dEeHh7jx483MzPTdUX6LDo6+q+//lLPDhs2zMHBQYf16CW5XB4ZGXnt2jU7O7uhQ4eKC3Nzc5cvX37//v0GDRqMHj1aJpPptkh9EhcXFxERERsbGxgYWK9ePSKKjIw8efKkeoVRo0bZ2dnprkC9kpCQsHfv3piYGHNz8759+zZr1ox09zGuPy3C8ePHb9q0ydfX98yZM3379tV1OXouPDx848aNGf9SqVS6rkgPrV+//osvvtixY8eKFSvUCwcMGHDixAlfX9/t27fr8AFkvdSvX7+1a9cuXLgwKipKXHLu3LktW7aoj3OO43RboT6ZM2dOeHh4vXr1TE1NO3TocPDgQXr9Y7xfv35lV42gF5KTk42NjePi4gRByM7OtrCwiI6O1nVR+mzVqlX9+vXTdRUGYcOGDS1bthSnb968aWZmlpWVJQjCs2fPjI2NExISdFqdHmrRosXGjRvF6V9//TU4OFi39egrpVKpnp4xY0bv3r3f/Bi/du1a2RSjJy3CyMhINzc3Nzc3IjIzM2vZsiWe/tG2Bw8ezJo1a9myZYmJibquxVBcuHChWbNmFhYWROTs7Ozh4REREaHrovTcnTt3Zs2atXz58pSUFF3XolcKP+394sULW1tbHX6M60kQJiUl2dvbq2cd/r+9+wdJb43jOJ6/oCJriQpLwhNSGVFGJxCFHCualEAMotZoSxoq+jNEomu7RFBDECE01VBDNZRE4FImCJaQLUWD9oewOxxu/LjT7w63557nvF/TOU4f5Mvz8cjjY2Mjq/N/qrGx0ePx1NbWnpycdHZ2JpNJ0YkMgTn/YRaLxeVy1dTUHB0dORyOVColOpGELi4utre3Z2dnBY63JJtlKisrf//Tiff396qqKoF5pOf3+/1+v3Y9NTW1tra2s7MjNpIRMOc/LBgMBoNB7XpycjISiWxsbIiNJJnr62ufzxeLxRwOx+XlpajxluSJ0Gq15nK5r78PB8jlclarVWwk4+jv77+/vxedwhC0Of++Zc5/kqqqd3d3olNIJZ1ODw4ORqNRbV+MwGVckiL0eDylUun4+LisrCydTieTyZGREdGhZPby8qJdfH5+7u3tOZ1OsXkMYnh4+ObmRjvc9fT09PX1dWBgQHQomX3P+cfHRzwe7+3tFZtHJtlsdmhoaHFx8funQSKX8Z/Zk/MDNjc36+vrA4FAc3NzOBwWHUdyXq/X5XL5fL7W1ta+vr7Hx0fRiSR0dnamqqqiKGazWVXVUCj09fUVjUabmpoCgUBDQ0MsFhOdUSozMzOqqprNZkVRVFXVtia53W6fz2ez2Vwu19PTk+iM8hgdHdUGWzM+Pv4lbhmX6qzRbDabTCbb29s7OjpEZ5FcoVC4urp6enpqaWlxOp2/fkny1cL/SrFYzOfz37fV1dUWi6WsrCyVSt3e3nZ3dyuKIiycjPL5fLFY/L61WCylUunq6ur5+dlms/X09JhMJoHxJPOPd7uyslL7IlTIMi5VEQIA8G/xQR4AYGgUIQDA0ChCAIChUYQAAEOjCAEAhkYRAgAMjSIEdG9sbMxut4dCod9fPDg4sNvtiURCVCpALyhCQPceHh4ymcz6+rp2+pqmUChkMpm3tzeBwQBdoAgBGbS1tdXV1a2srIgOAugPRQjIoLa2dmFhYXd39/z8XHQWQGcoQkAS09PTiqIsLy+LDgLoDEUISKKiomJpaenw8PDo6Eh0FkBPKEJAHhMTE11dXXNzcxymD/w5ihCQR3l5+erqaiKRiMfjorMAukERAlLx+/1ut3t+fv7z81N0FkAfKEJANpFIJJVKbW1tiQ4C6ANFCMjG6/UODg7u7++LDgLoA0UISCgcDptMJtEpAH0wsbsMAGBkPBECAAyNIgQAGBpFCAAwNIoQAGBoFCEAwNAoQgCAoVGEAABDowgBAIZGEQIADI0iBAAY2l+VFbj4JDsYNQAAAABJRU5ErkJggg==" }, "execution_count": 6, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# Plot eruption age estimate relative to rank-order plot of raw zircon ages\n", "h = plot(ylabel=\"Age (Ma)\", xlabel=\"N\", legend=:topleft, fg_color_legend=:white)\n", "plot!(h,1:length(ages),ages,yerror=uncert*2,seriestype=:scatter, label=\"Observed ages\")\n", "plot!(h,[length(ages)],[AgeEst],yerror=2*AgeEst_sigma, label=\"Bayesian lockup age estimate\",color=:red)\n", "plot!(h,collect(xlims()),[AgeEst,AgeEst],color=:red, label=\"\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "***\n", "## Getting your data out\n", "Acccessing the command line with Julia's `;` command-line mode, we can use the unix command `ls` to see what's on the local filesystem and look for any files we have written" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "; ls" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "If you're running in a Binder notebook, you can use a utility called `ffsend` to upload and transfer any files that you want. Edit the cells below to specify " ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# Download prebuilt ffsend linux binary\n", "download(\"https://github.com/timvisee/ffsend/releases/download/v0.2.46/ffsend-v0.2.46-linux-x64-static\", \"ffsend\")\n", "\n", "# Make ffsend executable\n", "run(`chmod +x ffsend`);" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "; ./ffsend upload histogram.pdf" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "***" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] } ], "metadata": { "kernelspec": { "display_name": "Julia 1.5.0-rc1", "language": "julia", "name": "julia-1.5" }, "language_info": { "file_extension": ".jl", "mimetype": "application/julia", "name": "julia", "version": "1.5.0" } }, "nbformat": 4, "nbformat_minor": 2 }