Last active
January 13, 2022 12:25
-
-
Save jamm1985/79b43a6dc1655a32e44ae64ee437d5db to your computer and use it in GitHub Desktop.
magnitude_det_rate_estimation_mle.ipynb
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
{ | |
"nbformat": 4, | |
"nbformat_minor": 0, | |
"metadata": { | |
"colab": { | |
"name": "magnitude_det_rate_estimation_mle.ipynb", | |
"provenance": [], | |
"collapsed_sections": [], | |
"mount_file_id": "1bbYs7f1hAQ9UMgj9Ps_6HmLdc_WLn6b_", | |
"authorship_tag": "ABX9TyMMVraKNqNI/3i6HW2axdHa", | |
"include_colab_link": true | |
}, | |
"kernelspec": { | |
"name": "python3", | |
"display_name": "Python 3" | |
}, | |
"language_info": { | |
"name": "python" | |
} | |
}, | |
"cells": [ | |
{ | |
"cell_type": "markdown", | |
"metadata": { | |
"id": "view-in-github", | |
"colab_type": "text" | |
}, | |
"source": [ | |
"<a href=\"https://colab.research.google.com/gist/jamm1985/79b43a6dc1655a32e44ae64ee437d5db/magnitude_det_rate_estimation_mle.ipynb\" target=\"_parent\"><img src=\"https://colab.research.google.com/assets/colab-badge.svg\" alt=\"Open In Colab\"/></a>" | |
] | |
}, | |
{ | |
"cell_type": "code", | |
"metadata": { | |
"colab": { | |
"base_uri": "https://localhost:8080/" | |
}, | |
"id": "HWLL0dmj3hs9", | |
"outputId": "7d4250ea-0ff7-4a63-ae52-ebc996b5ee82" | |
}, | |
"source": [ | |
"%pip install scipy==1.7.3 " | |
], | |
"execution_count": null, | |
"outputs": [ | |
{ | |
"output_type": "stream", | |
"name": "stdout", | |
"text": [ | |
"Collecting scipy==1.7.3\n", | |
" Downloading scipy-1.7.3-cp37-cp37m-manylinux_2_12_x86_64.manylinux2010_x86_64.whl (38.1 MB)\n", | |
"\u001b[K |████████████████████████████████| 38.1 MB 1.3 MB/s \n", | |
"\u001b[?25hRequirement already satisfied: numpy<1.23.0,>=1.16.5 in /usr/local/lib/python3.7/dist-packages (from scipy==1.7.3) (1.19.5)\n", | |
"Installing collected packages: scipy\n", | |
" Attempting uninstall: scipy\n", | |
" Found existing installation: scipy 1.4.1\n", | |
" Uninstalling scipy-1.4.1:\n", | |
" Successfully uninstalled scipy-1.4.1\n", | |
"\u001b[31mERROR: pip's dependency resolver does not currently take into account all the packages that are installed. This behaviour is the source of the following dependency conflicts.\n", | |
"albumentations 0.1.12 requires imgaug<0.2.7,>=0.2.5, but you have imgaug 0.2.9 which is incompatible.\u001b[0m\n", | |
"Successfully installed scipy-1.7.3\n" | |
] | |
} | |
] | |
}, | |
{ | |
"cell_type": "code", | |
"metadata": { | |
"id": "6XqC8G4J3u8R" | |
}, | |
"source": [ | |
"import numpy as np\n", | |
"import matplotlib.pyplot as plt\n", | |
"import pandas as pd\n", | |
"\n", | |
"from scipy import stats\n", | |
"from scipy import special\n", | |
"from scipy.stats import norm\n", | |
"from scipy.optimize import minimize" | |
], | |
"execution_count": null, | |
"outputs": [] | |
}, | |
{ | |
"cell_type": "markdown", | |
"metadata": { | |
"id": "HNypu_ESwwBV" | |
}, | |
"source": [ | |
"### get $M_1,..,M_n$" | |
] | |
}, | |
{ | |
"cell_type": "code", | |
"metadata": { | |
"id": "jN_KpGnJ39CH" | |
}, | |
"source": [ | |
"catalog = pd.read_csv(\"/content/drive/MyDrive/Colab Notebooks/Magnitude detection rate estimation (MLE)/catalog.txt\",\n", | |
" sep = \"\\t\",\n", | |
" names=[\"date\",\"lat\",\"lon\",\"depth\",\"ML\"])\n", | |
"catalog['date'] = pd.to_datetime(catalog['date'], format='%Y %m %d %H %M %S')\n", | |
"M_observation = catalog['ML'].to_numpy()" | |
], | |
"execution_count": null, | |
"outputs": [] | |
}, | |
{ | |
"cell_type": "code", | |
"metadata": { | |
"colab": { | |
"base_uri": "https://localhost:8080/", | |
"height": 300 | |
}, | |
"id": "95qjbM01HdpA", | |
"outputId": "a2d0eb41-a565-4b18-8431-1fab58c38786" | |
}, | |
"source": [ | |
"catalog.describe()" | |
], | |
"execution_count": null, | |
"outputs": [ | |
{ | |
"output_type": "execute_result", | |
"data": { | |
"text/html": [ | |
"\n", | |
" <div id=\"df-e9703547-02d9-4362-aea2-d41aca48fc33\">\n", | |
" <div class=\"colab-df-container\">\n", | |
" <div>\n", | |
"<style scoped>\n", | |
" .dataframe tbody tr th:only-of-type {\n", | |
" vertical-align: middle;\n", | |
" }\n", | |
"\n", | |
" .dataframe tbody tr th {\n", | |
" vertical-align: top;\n", | |
" }\n", | |
"\n", | |
" .dataframe thead th {\n", | |
" text-align: right;\n", | |
" }\n", | |
"</style>\n", | |
"<table border=\"1\" class=\"dataframe\">\n", | |
" <thead>\n", | |
" <tr style=\"text-align: right;\">\n", | |
" <th></th>\n", | |
" <th>lat</th>\n", | |
" <th>lon</th>\n", | |
" <th>depth</th>\n", | |
" <th>ML</th>\n", | |
" </tr>\n", | |
" </thead>\n", | |
" <tbody>\n", | |
" <tr>\n", | |
" <th>count</th>\n", | |
" <td>437.000000</td>\n", | |
" <td>437.000000</td>\n", | |
" <td>437.000000</td>\n", | |
" <td>437.000000</td>\n", | |
" </tr>\n", | |
" <tr>\n", | |
" <th>mean</th>\n", | |
" <td>51.472320</td>\n", | |
" <td>143.194384</td>\n", | |
" <td>11.853547</td>\n", | |
" <td>1.870023</td>\n", | |
" </tr>\n", | |
" <tr>\n", | |
" <th>std</th>\n", | |
" <td>0.124103</td>\n", | |
" <td>0.185499</td>\n", | |
" <td>5.655157</td>\n", | |
" <td>0.739913</td>\n", | |
" </tr>\n", | |
" <tr>\n", | |
" <th>min</th>\n", | |
" <td>51.100000</td>\n", | |
" <td>143.000000</td>\n", | |
" <td>0.000000</td>\n", | |
" <td>0.100000</td>\n", | |
" </tr>\n", | |
" <tr>\n", | |
" <th>25%</th>\n", | |
" <td>51.393000</td>\n", | |
" <td>143.082000</td>\n", | |
" <td>8.000000</td>\n", | |
" <td>1.300000</td>\n", | |
" </tr>\n", | |
" <tr>\n", | |
" <th>50%</th>\n", | |
" <td>51.537000</td>\n", | |
" <td>143.104000</td>\n", | |
" <td>12.300000</td>\n", | |
" <td>1.800000</td>\n", | |
" </tr>\n", | |
" <tr>\n", | |
" <th>75%</th>\n", | |
" <td>51.570000</td>\n", | |
" <td>143.269000</td>\n", | |
" <td>15.100000</td>\n", | |
" <td>2.200000</td>\n", | |
" </tr>\n", | |
" <tr>\n", | |
" <th>max</th>\n", | |
" <td>51.600000</td>\n", | |
" <td>144.000000</td>\n", | |
" <td>31.700000</td>\n", | |
" <td>4.500000</td>\n", | |
" </tr>\n", | |
" </tbody>\n", | |
"</table>\n", | |
"</div>\n", | |
" <button class=\"colab-df-convert\" onclick=\"convertToInteractive('df-e9703547-02d9-4362-aea2-d41aca48fc33')\"\n", | |
" title=\"Convert this dataframe to an interactive table.\"\n", | |
" style=\"display:none;\">\n", | |
" \n", | |
" <svg xmlns=\"http://www.w3.org/2000/svg\" height=\"24px\"viewBox=\"0 0 24 24\"\n", | |
" width=\"24px\">\n", | |
" <path d=\"M0 0h24v24H0V0z\" fill=\"none\"/>\n", | |
" <path d=\"M18.56 5.44l.94 2.06.94-2.06 2.06-.94-2.06-.94-.94-2.06-.94 2.06-2.06.94zm-11 1L8.5 8.5l.94-2.06 2.06-.94-2.06-.94L8.5 2.5l-.94 2.06-2.06.94zm10 10l.94 2.06.94-2.06 2.06-.94-2.06-.94-.94-2.06-.94 2.06-2.06.94z\"/><path d=\"M17.41 7.96l-1.37-1.37c-.4-.4-.92-.59-1.43-.59-.52 0-1.04.2-1.43.59L10.3 9.45l-7.72 7.72c-.78.78-.78 2.05 0 2.83L4 21.41c.39.39.9.59 1.41.59.51 0 1.02-.2 1.41-.59l7.78-7.78 2.81-2.81c.8-.78.8-2.07 0-2.86zM5.41 20L4 18.59l7.72-7.72 1.47 1.35L5.41 20z\"/>\n", | |
" </svg>\n", | |
" </button>\n", | |
" \n", | |
" <style>\n", | |
" .colab-df-container {\n", | |
" display:flex;\n", | |
" flex-wrap:wrap;\n", | |
" gap: 12px;\n", | |
" }\n", | |
"\n", | |
" .colab-df-convert {\n", | |
" background-color: #E8F0FE;\n", | |
" border: none;\n", | |
" border-radius: 50%;\n", | |
" cursor: pointer;\n", | |
" display: none;\n", | |
" fill: #1967D2;\n", | |
" height: 32px;\n", | |
" padding: 0 0 0 0;\n", | |
" width: 32px;\n", | |
" }\n", | |
"\n", | |
" .colab-df-convert:hover {\n", | |
" background-color: #E2EBFA;\n", | |
" box-shadow: 0px 1px 2px rgba(60, 64, 67, 0.3), 0px 1px 3px 1px rgba(60, 64, 67, 0.15);\n", | |
" fill: #174EA6;\n", | |
" }\n", | |
"\n", | |
" [theme=dark] .colab-df-convert {\n", | |
" background-color: #3B4455;\n", | |
" fill: #D2E3FC;\n", | |
" }\n", | |
"\n", | |
" [theme=dark] .colab-df-convert:hover {\n", | |
" background-color: #434B5C;\n", | |
" box-shadow: 0px 1px 3px 1px rgba(0, 0, 0, 0.15);\n", | |
" filter: drop-shadow(0px 1px 2px rgba(0, 0, 0, 0.3));\n", | |
" fill: #FFFFFF;\n", | |
" }\n", | |
" </style>\n", | |
"\n", | |
" <script>\n", | |
" const buttonEl =\n", | |
" document.querySelector('#df-e9703547-02d9-4362-aea2-d41aca48fc33 button.colab-df-convert');\n", | |
" buttonEl.style.display =\n", | |
" google.colab.kernel.accessAllowed ? 'block' : 'none';\n", | |
"\n", | |
" async function convertToInteractive(key) {\n", | |
" const element = document.querySelector('#df-e9703547-02d9-4362-aea2-d41aca48fc33');\n", | |
" const dataTable =\n", | |
" await google.colab.kernel.invokeFunction('convertToInteractive',\n", | |
" [key], {});\n", | |
" if (!dataTable) return;\n", | |
"\n", | |
" const docLinkHtml = 'Like what you see? Visit the ' +\n", | |
" '<a target=\"_blank\" href=https://colab.research.google.com/notebooks/data_table.ipynb>data table notebook</a>'\n", | |
" + ' to learn more about interactive tables.';\n", | |
" element.innerHTML = '';\n", | |
" dataTable['output_type'] = 'display_data';\n", | |
" await google.colab.output.renderOutput(dataTable, element);\n", | |
" const docLink = document.createElement('div');\n", | |
" docLink.innerHTML = docLinkHtml;\n", | |
" element.appendChild(docLink);\n", | |
" }\n", | |
" </script>\n", | |
" </div>\n", | |
" </div>\n", | |
" " | |
], | |
"text/plain": [ | |
" lat lon depth ML\n", | |
"count 437.000000 437.000000 437.000000 437.000000\n", | |
"mean 51.472320 143.194384 11.853547 1.870023\n", | |
"std 0.124103 0.185499 5.655157 0.739913\n", | |
"min 51.100000 143.000000 0.000000 0.100000\n", | |
"25% 51.393000 143.082000 8.000000 1.300000\n", | |
"50% 51.537000 143.104000 12.300000 1.800000\n", | |
"75% 51.570000 143.269000 15.100000 2.200000\n", | |
"max 51.600000 144.000000 31.700000 4.500000" | |
] | |
}, | |
"metadata": {}, | |
"execution_count": 5 | |
} | |
] | |
}, | |
{ | |
"cell_type": "code", | |
"metadata": { | |
"colab": { | |
"base_uri": "https://localhost:8080/", | |
"height": 424 | |
}, | |
"id": "bFoNzGIattM1", | |
"outputId": "f792de53-4097-448c-de99-6a93f1342f05" | |
}, | |
"source": [ | |
"catalog" | |
], | |
"execution_count": null, | |
"outputs": [ | |
{ | |
"output_type": "execute_result", | |
"data": { | |
"text/html": [ | |
"\n", | |
" <div id=\"df-820923d3-2fc7-49b3-b3e3-2efaf9947906\">\n", | |
" <div class=\"colab-df-container\">\n", | |
" <div>\n", | |
"<style scoped>\n", | |
" .dataframe tbody tr th:only-of-type {\n", | |
" vertical-align: middle;\n", | |
" }\n", | |
"\n", | |
" .dataframe tbody tr th {\n", | |
" vertical-align: top;\n", | |
" }\n", | |
"\n", | |
" .dataframe thead th {\n", | |
" text-align: right;\n", | |
" }\n", | |
"</style>\n", | |
"<table border=\"1\" class=\"dataframe\">\n", | |
" <thead>\n", | |
" <tr style=\"text-align: right;\">\n", | |
" <th></th>\n", | |
" <th>date</th>\n", | |
" <th>lat</th>\n", | |
" <th>lon</th>\n", | |
" <th>depth</th>\n", | |
" <th>ML</th>\n", | |
" </tr>\n", | |
" </thead>\n", | |
" <tbody>\n", | |
" <tr>\n", | |
" <th>0</th>\n", | |
" <td>1966-09-14 05:42:29</td>\n", | |
" <td>51.600</td>\n", | |
" <td>143.900</td>\n", | |
" <td>15.0</td>\n", | |
" <td>3.6</td>\n", | |
" </tr>\n", | |
" <tr>\n", | |
" <th>1</th>\n", | |
" <td>1967-01-09 20:26:09</td>\n", | |
" <td>51.500</td>\n", | |
" <td>143.500</td>\n", | |
" <td>10.0</td>\n", | |
" <td>3.6</td>\n", | |
" </tr>\n", | |
" <tr>\n", | |
" <th>2</th>\n", | |
" <td>1967-05-21 03:52:25</td>\n", | |
" <td>51.400</td>\n", | |
" <td>144.000</td>\n", | |
" <td>15.0</td>\n", | |
" <td>3.4</td>\n", | |
" </tr>\n", | |
" <tr>\n", | |
" <th>3</th>\n", | |
" <td>1968-01-25 23:53:50</td>\n", | |
" <td>51.300</td>\n", | |
" <td>143.400</td>\n", | |
" <td>15.0</td>\n", | |
" <td>3.6</td>\n", | |
" </tr>\n", | |
" <tr>\n", | |
" <th>4</th>\n", | |
" <td>1970-11-04 11:34:53</td>\n", | |
" <td>51.200</td>\n", | |
" <td>143.400</td>\n", | |
" <td>15.0</td>\n", | |
" <td>3.3</td>\n", | |
" </tr>\n", | |
" <tr>\n", | |
" <th>...</th>\n", | |
" <td>...</td>\n", | |
" <td>...</td>\n", | |
" <td>...</td>\n", | |
" <td>...</td>\n", | |
" <td>...</td>\n", | |
" </tr>\n", | |
" <tr>\n", | |
" <th>432</th>\n", | |
" <td>2020-03-30 07:30:20</td>\n", | |
" <td>51.398</td>\n", | |
" <td>143.183</td>\n", | |
" <td>7.5</td>\n", | |
" <td>2.6</td>\n", | |
" </tr>\n", | |
" <tr>\n", | |
" <th>433</th>\n", | |
" <td>2020-04-27 20:51:13</td>\n", | |
" <td>51.581</td>\n", | |
" <td>143.081</td>\n", | |
" <td>11.4</td>\n", | |
" <td>2.1</td>\n", | |
" </tr>\n", | |
" <tr>\n", | |
" <th>434</th>\n", | |
" <td>2020-05-13 17:58:42</td>\n", | |
" <td>51.560</td>\n", | |
" <td>143.103</td>\n", | |
" <td>17.1</td>\n", | |
" <td>2.7</td>\n", | |
" </tr>\n", | |
" <tr>\n", | |
" <th>435</th>\n", | |
" <td>2020-05-14 19:02:06</td>\n", | |
" <td>51.561</td>\n", | |
" <td>143.096</td>\n", | |
" <td>15.7</td>\n", | |
" <td>1.8</td>\n", | |
" </tr>\n", | |
" <tr>\n", | |
" <th>436</th>\n", | |
" <td>2020-05-28 12:50:16</td>\n", | |
" <td>51.561</td>\n", | |
" <td>143.093</td>\n", | |
" <td>14.2</td>\n", | |
" <td>1.2</td>\n", | |
" </tr>\n", | |
" </tbody>\n", | |
"</table>\n", | |
"<p>437 rows × 5 columns</p>\n", | |
"</div>\n", | |
" <button class=\"colab-df-convert\" onclick=\"convertToInteractive('df-820923d3-2fc7-49b3-b3e3-2efaf9947906')\"\n", | |
" title=\"Convert this dataframe to an interactive table.\"\n", | |
" style=\"display:none;\">\n", | |
" \n", | |
" <svg xmlns=\"http://www.w3.org/2000/svg\" height=\"24px\"viewBox=\"0 0 24 24\"\n", | |
" width=\"24px\">\n", | |
" <path d=\"M0 0h24v24H0V0z\" fill=\"none\"/>\n", | |
" <path d=\"M18.56 5.44l.94 2.06.94-2.06 2.06-.94-2.06-.94-.94-2.06-.94 2.06-2.06.94zm-11 1L8.5 8.5l.94-2.06 2.06-.94-2.06-.94L8.5 2.5l-.94 2.06-2.06.94zm10 10l.94 2.06.94-2.06 2.06-.94-2.06-.94-.94-2.06-.94 2.06-2.06.94z\"/><path d=\"M17.41 7.96l-1.37-1.37c-.4-.4-.92-.59-1.43-.59-.52 0-1.04.2-1.43.59L10.3 9.45l-7.72 7.72c-.78.78-.78 2.05 0 2.83L4 21.41c.39.39.9.59 1.41.59.51 0 1.02-.2 1.41-.59l7.78-7.78 2.81-2.81c.8-.78.8-2.07 0-2.86zM5.41 20L4 18.59l7.72-7.72 1.47 1.35L5.41 20z\"/>\n", | |
" </svg>\n", | |
" </button>\n", | |
" \n", | |
" <style>\n", | |
" .colab-df-container {\n", | |
" display:flex;\n", | |
" flex-wrap:wrap;\n", | |
" gap: 12px;\n", | |
" }\n", | |
"\n", | |
" .colab-df-convert {\n", | |
" background-color: #E8F0FE;\n", | |
" border: none;\n", | |
" border-radius: 50%;\n", | |
" cursor: pointer;\n", | |
" display: none;\n", | |
" fill: #1967D2;\n", | |
" height: 32px;\n", | |
" padding: 0 0 0 0;\n", | |
" width: 32px;\n", | |
" }\n", | |
"\n", | |
" .colab-df-convert:hover {\n", | |
" background-color: #E2EBFA;\n", | |
" box-shadow: 0px 1px 2px rgba(60, 64, 67, 0.3), 0px 1px 3px 1px rgba(60, 64, 67, 0.15);\n", | |
" fill: #174EA6;\n", | |
" }\n", | |
"\n", | |
" [theme=dark] .colab-df-convert {\n", | |
" background-color: #3B4455;\n", | |
" fill: #D2E3FC;\n", | |
" }\n", | |
"\n", | |
" [theme=dark] .colab-df-convert:hover {\n", | |
" background-color: #434B5C;\n", | |
" box-shadow: 0px 1px 3px 1px rgba(0, 0, 0, 0.15);\n", | |
" filter: drop-shadow(0px 1px 2px rgba(0, 0, 0, 0.3));\n", | |
" fill: #FFFFFF;\n", | |
" }\n", | |
" </style>\n", | |
"\n", | |
" <script>\n", | |
" const buttonEl =\n", | |
" document.querySelector('#df-820923d3-2fc7-49b3-b3e3-2efaf9947906 button.colab-df-convert');\n", | |
" buttonEl.style.display =\n", | |
" google.colab.kernel.accessAllowed ? 'block' : 'none';\n", | |
"\n", | |
" async function convertToInteractive(key) {\n", | |
" const element = document.querySelector('#df-820923d3-2fc7-49b3-b3e3-2efaf9947906');\n", | |
" const dataTable =\n", | |
" await google.colab.kernel.invokeFunction('convertToInteractive',\n", | |
" [key], {});\n", | |
" if (!dataTable) return;\n", | |
"\n", | |
" const docLinkHtml = 'Like what you see? Visit the ' +\n", | |
" '<a target=\"_blank\" href=https://colab.research.google.com/notebooks/data_table.ipynb>data table notebook</a>'\n", | |
" + ' to learn more about interactive tables.';\n", | |
" element.innerHTML = '';\n", | |
" dataTable['output_type'] = 'display_data';\n", | |
" await google.colab.output.renderOutput(dataTable, element);\n", | |
" const docLink = document.createElement('div');\n", | |
" docLink.innerHTML = docLinkHtml;\n", | |
" element.appendChild(docLink);\n", | |
" }\n", | |
" </script>\n", | |
" </div>\n", | |
" </div>\n", | |
" " | |
], | |
"text/plain": [ | |
" date lat lon depth ML\n", | |
"0 1966-09-14 05:42:29 51.600 143.900 15.0 3.6\n", | |
"1 1967-01-09 20:26:09 51.500 143.500 10.0 3.6\n", | |
"2 1967-05-21 03:52:25 51.400 144.000 15.0 3.4\n", | |
"3 1968-01-25 23:53:50 51.300 143.400 15.0 3.6\n", | |
"4 1970-11-04 11:34:53 51.200 143.400 15.0 3.3\n", | |
".. ... ... ... ... ...\n", | |
"432 2020-03-30 07:30:20 51.398 143.183 7.5 2.6\n", | |
"433 2020-04-27 20:51:13 51.581 143.081 11.4 2.1\n", | |
"434 2020-05-13 17:58:42 51.560 143.103 17.1 2.7\n", | |
"435 2020-05-14 19:02:06 51.561 143.096 15.7 1.8\n", | |
"436 2020-05-28 12:50:16 51.561 143.093 14.2 1.2\n", | |
"\n", | |
"[437 rows x 5 columns]" | |
] | |
}, | |
"metadata": {}, | |
"execution_count": 6 | |
} | |
] | |
}, | |
{ | |
"cell_type": "markdown", | |
"metadata": { | |
"id": "xXGnsppa8tbw" | |
}, | |
"source": [ | |
"## log-likelihood for constant rate\n", | |
"\n", | |
"Based of GR law, the distribution of $M$ is \n", | |
"\n", | |
"$$w(M)=A \\exp(-\\beta M) \\approx \\exp(-\\beta M)$$\n", | |
"\n", | |
"where $\\beta=b \\log 10$\n", | |
"\n", | |
"The detection rate function (Ogata, Kastura, 1992):\n", | |
"\n", | |
"$$q(M)=\\mathrm{erf}(M|\\mu,\\sigma)=\\int_{-\\infty}^M\\frac{1}{\\sqrt{2\\pi}\\sigma}e^{\\frac{-(x-\\mu)^2}{2\\sigma}}dx$$\n", | |
"\n", | |
"Given the observer magnitudes $M_1,...,M_n$ the likelihood of the point process is (Ogata, Kastura, 1992):\n", | |
"\n", | |
"$$\\log L(\\theta)=\\Sigma_{i=1}^n \\log f(M_i|\\theta)=$$\n", | |
"\n", | |
"$$N\\log(\\beta)-\\beta\\Sigma_{i=1}^N[M_i + \\log \\mathrm{erf}\\{M_i|\\mu,\\sigma\\}] + \\exp \\left[ N\\beta\\mu - \\frac{N}{2}\\beta^2\\mu^2\\right]$$\n", | |
"\n", | |
"where\n", | |
"\n", | |
"$$f(M_i|\\theta)=\\frac{w(M)q(M)}{\\int_{-\\infty}^{+\\infty}w(M)q(M)dm}$$\n", | |
"\n", | |
"$$\\theta=[\\beta,\\mu,\\sigma]^T$$\n", | |
"\n", | |
"and erf function\n", | |
"\n", | |
"$$F_X(x)=\\Phi(\\frac{x-\\mu}{\\sigma})=\\frac{1}{2} \\left[1+ \\mathrm{erf}(\\frac{x-\\mu}{2\\sqrt{\\sigma}}) \\right] $$" | |
] | |
}, | |
{ | |
"cell_type": "markdown", | |
"metadata": { | |
"id": "Ry5fe9EFCult" | |
}, | |
"source": [ | |
"## evaluation function from observation and given parameters" | |
] | |
}, | |
{ | |
"cell_type": "code", | |
"metadata": { | |
"id": "cRW0jBC82spP" | |
}, | |
"source": [ | |
"def density(M,theta):\n", | |
" beta = theta[0]\n", | |
" mu = theta[1]\n", | |
" sigma = theta[2]\n", | |
" constant = (np.exp(-beta*M)*norm.cdf(M,mu,sigma)).sum()\n", | |
" return np.exp(-beta*M)*norm.cdf(M,mu,sigma)/constant" | |
], | |
"execution_count": null, | |
"outputs": [] | |
}, | |
{ | |
"cell_type": "code", | |
"metadata": { | |
"id": "6X7c8A0f3Mrz" | |
}, | |
"source": [ | |
"def log_l(theta):\n", | |
" beta = theta[0]\n", | |
" mu = theta[1]\n", | |
" sigma = theta[2]\n", | |
" N = M_observation.size\n", | |
" return (N * np.log(beta) - beta * (M_observation + np.log(norm.cdf(M_observation,mu,sigma))).sum() + np.exp(N*beta*mu - N/2 * beta**2 * mu**2))" | |
], | |
"execution_count": null, | |
"outputs": [] | |
}, | |
{ | |
"cell_type": "markdown", | |
"metadata": { | |
"id": "bp9n0HcmQhZm" | |
}, | |
"source": [ | |
"## first derivatives of $\\log L$\n", | |
"\n", | |
"$\\log L (\\theta) = N\\log(\\beta)-\\beta\\Sigma_{i=1}^N[M_i + \\log \\mathrm{erf}\\{M_i|\\mu,\\sigma\\}] + \\exp \\left[ N\\beta\\mu - \\frac{N}{2}\\beta^2\\mu^2\\right]$\n", | |
"\n", | |
"$\\frac{\\partial \\log L}{\\partial \\beta} = \\frac{N}{\\beta} - \\Sigma_{i=1}^N \\left[M_i + \\log\\Phi \\left( \\frac{M_i - \\mu}{\\sigma} \\right) \\right] + \\exp \\left[ N\\beta\\mu - \\frac{N}{2}\\beta^2\\mu^2\\right](N\\mu - N\\beta\\mu^2)$\n", | |
"\n", | |
"$\\frac{\\partial \\log L}{\\partial \\mu} = \\beta N \\frac{1}{\\sigma} + \\exp \\left[ N\\beta\\mu - \\frac{N}{2}\\beta^2\\mu^2\\right] (N \\beta - N\\beta^2\\mu)$ \n", | |
"\n", | |
"$\\frac{\\partial \\log L}{\\partial \\sigma} = \\beta \\Sigma_{i=1}^N \\left[ \\frac{M_i - \\mu}{\\sigma} \\right]$" | |
] | |
}, | |
{ | |
"cell_type": "code", | |
"metadata": { | |
"id": "96b6iajScm-G" | |
}, | |
"source": [ | |
"def der_log_1(theta):\n", | |
" beta = theta[0]\n", | |
" mu = theta[1]\n", | |
" sigma = theta[2]\n", | |
" N = M_observation.size\n", | |
" der = np.zeros_like(theta)\n", | |
" exp_term = np.exp(N*beta*mu - N/2 * beta**2 * mu**2)\n", | |
" exp_term_beta_prime = N*mu - N*beta*mu**2\n", | |
" exp_term_mu_prime = N*beta - N*beta**2*mu\n", | |
" der[0] = N/beta - (M_observation + np.log(norm.cdf(M_observation,mu,sigma))).sum() + exp_term*exp_term_beta_prime\n", | |
" der[1] = N*beta*(1/sigma) + exp_term * exp_term_mu_prime\n", | |
" der[2] = beta* (((M_observation - mu)/sigma**2).sum())\n", | |
" return der" | |
], | |
"execution_count": null, | |
"outputs": [] | |
}, | |
{ | |
"cell_type": "code", | |
"metadata": { | |
"id": "J25gMxwLYq8O" | |
}, | |
"source": [ | |
"def hess_log_1(theta):\n", | |
" beta = theta[0]\n", | |
" mu = theta[1]\n", | |
" sigma = theta[2]\n", | |
" N = M_observation.size\n", | |
" hess = np.zeros((3,3))\n", | |
" exp_term = np.exp(N*beta*mu - N/2 * beta**2 * mu**2)\n", | |
" exp_term_beta_prime = N*mu - N*beta*mu**2\n", | |
" exp_term_mu_prime = N*beta - N*beta**2*mu\n", | |
" \n", | |
" hess[0,0] = -N/beta**2 + exp_term_beta_prime**2*exp_term - N*mu**2*exp_term\n", | |
" hess[0,1] = N*(1/sigma) + exp_term_beta_prime*exp_term*exp_term_mu_prime + exp_term*(N-2*N*beta*mu)\n", | |
" hess[0,2] = ((M_observation - mu)/sigma**2).sum()\n", | |
"\n", | |
" hess[1,0] = N/sigma**2 + exp_term_mu_prime*exp_term*exp_term_beta_prime + exp_term*(N-2*N*beta*mu)\n", | |
" hess[1,1] = exp_term_mu_prime**2*exp_term - N*beta**2*exp_term\n", | |
" hess[1,2] = -beta*N*(1/sigma**2)\n", | |
"\n", | |
" hess[2,0] = ((M_observation - mu)/sigma).sum()\n", | |
" hess[2,1] = -N*beta*(1/sigma)\n", | |
" hess[2,2] = -beta*((M_observation - mu)/sigma**2).sum()\n", | |
"\n", | |
" return hess" | |
], | |
"execution_count": null, | |
"outputs": [] | |
}, | |
{ | |
"cell_type": "code", | |
"metadata": { | |
"id": "BCb_vT_8Vucb" | |
}, | |
"source": [ | |
"# estimate b-value \n", | |
"#this is modified version of \n", | |
"#source: https://github.com/qingkaikong/Learning_Obspy/blob/master/04_catalogAnalysis_20140421/calc_rec.py\n", | |
"\n", | |
"def calc_recurrence(magnitudes, min_mag = None, max_mag = None, interval = 0.05):\n", | |
"\n", | |
" \"\"\"\n", | |
" This function reads an earthquake catalogue file and calculates the\n", | |
" Gutenberg-Richter recurrence parameters using maximum likelihood (Aki 1965) \n", | |
" approach.\n", | |
"\n", | |
" Funtion arguments:\n", | |
"\n", | |
" magnitudes: list of magnitudes values\n", | |
" min_mag: minimum magnitude for which data will be used - i.e. catalogue\n", | |
" completeness\n", | |
" max_mag: maximum magnitude used in bounded G-R curve. If not specified,\n", | |
" defined as the maximum magnitude in the catlogue + 0.1 magnitude\n", | |
" units.\n", | |
" interval: Width of magnitude bins for generating cumulative histogram\n", | |
" of earthquake recurrence. Default avlue 0.1 magnitude units. \n", | |
" \n", | |
" \"\"\"\n", | |
" \n", | |
" # If minimum magnitude is not specified, read all magnitudes\n", | |
" if min_mag is not None:\n", | |
" pass\n", | |
" else:\n", | |
" min_mag = -1.0\n", | |
" \n", | |
" # If minimum magnitude is not specified default value to minimum in catalogue\n", | |
" if min_mag == -1.0:\n", | |
" min_mag = min(magnitudes)\n", | |
" # If maximum magnitude is not specified default value to maximum in catalogue\n", | |
" if max_mag is not None:\n", | |
" pass\n", | |
" else:\n", | |
" max_mag = max(magnitudes) + 0.1\n", | |
"\n", | |
" num_eq = len(magnitudes)\n", | |
" print('Minimum magnitude: {}'.format(min_mag))\n", | |
" print('Total number of earthquakes: {}'.format(num_eq))\n", | |
" #num_years = max(years)-min(years)\n", | |
" #annual_num_eq = num_eq/num_years\n", | |
" #print 'Annual number of earthquakes greater than Mw', min_mag,':', \\\n", | |
" #annual_num_eq\n", | |
" print('Maximum catalog magnitude: {}'.format(max(magnitudes)))\n", | |
" print('Mmax = {}'.format(max_mag))\n", | |
" max_mag_bin = max(magnitudes) + 0.15\n", | |
" \n", | |
" # Magnitude bins\n", | |
" bins = np.arange(min_mag, max_mag_bin, interval)\n", | |
" # Magnitude bins for plotting - we will re-arrange bins later\n", | |
" plot_bins = np.arange(min_mag, max_mag, interval)\n", | |
"\n", | |
" ###########################################################################\n", | |
" # Generate distribution\n", | |
" ###########################################################################\n", | |
" # Generate histogram\n", | |
" hist = np.histogram(magnitudes,bins=bins)\n", | |
"\n", | |
" # Reverse array order\n", | |
" hist = hist[0][::-1]\n", | |
" bins = bins[::-1]\n", | |
"\n", | |
" # Calculate cumulative sum\n", | |
" cum_hist = hist.cumsum()\n", | |
" # Ensure bins have the same length has the cumulative histogram.\n", | |
" # Remove the upper bound for the highest interval.\n", | |
" bins = bins[1:]\n", | |
"\n", | |
" # Maximum Likelihood Estimator fitting\n", | |
" # b value\n", | |
" b_mle = np.log10(np.exp(1)) / (np.mean(magnitudes) - min_mag)\n", | |
" beta_mle = np.log(10) * b_mle\n", | |
" print('Maximum Likelihood: b value {}'.format( b_mle))\n", | |
"\n", | |
" return cum_hist[::-1], bins[::-1]" | |
], | |
"execution_count": null, | |
"outputs": [] | |
}, | |
{ | |
"cell_type": "markdown", | |
"metadata": { | |
"id": "1e0HSbKVVfH6" | |
}, | |
"source": [ | |
"## before 2006" | |
] | |
}, | |
{ | |
"cell_type": "code", | |
"metadata": { | |
"colab": { | |
"base_uri": "https://localhost:8080/" | |
}, | |
"id": "I87zkEpPt92J", | |
"outputId": "d2d376af-bec6-42e5-a438-2a07df696d38" | |
}, | |
"source": [ | |
"M_observation_before_2006 = catalog[catalog['date'] < '20060917']['ML']\n", | |
"print(M_observation_before_2006.describe())\n", | |
"M_observation_before_2006 = M_observation_before_2006.to_numpy()" | |
], | |
"execution_count": null, | |
"outputs": [ | |
{ | |
"output_type": "stream", | |
"name": "stdout", | |
"text": [ | |
"count 22.000000\n", | |
"mean 3.477273\n", | |
"std 0.387829\n", | |
"min 3.000000\n", | |
"25% 3.225000\n", | |
"50% 3.400000\n", | |
"75% 3.600000\n", | |
"max 4.500000\n", | |
"Name: ML, dtype: float64\n" | |
] | |
} | |
] | |
}, | |
{ | |
"cell_type": "code", | |
"metadata": { | |
"colab": { | |
"base_uri": "https://localhost:8080/" | |
}, | |
"id": "tJZabdRGccx4", | |
"outputId": "d09c2d70-9256-4320-9c12-291532ae0092" | |
}, | |
"source": [ | |
"calc_recurrence(M_observation_before_2006.tolist(), min_mag = 2.0, max_mag = None, interval = 0.05)" | |
], | |
"execution_count": null, | |
"outputs": [ | |
{ | |
"output_type": "stream", | |
"name": "stdout", | |
"text": [ | |
"Minimum magnitude: 2.0\n", | |
"Total number of earthquakes: 22\n", | |
"Maximum catalog magnitude: 4.5\n", | |
"Mmax = 4.6\n", | |
"Maximum Likelihood: b value 0.2939839569806628\n" | |
] | |
}, | |
{ | |
"output_type": "execute_result", | |
"data": { | |
"text/plain": [ | |
"(array([22, 22, 22, 22, 22, 22, 22, 22, 22, 22, 22, 22, 22, 22, 22, 22, 22,\n", | |
" 22, 22, 22, 22, 19, 19, 17, 17, 16, 16, 13, 13, 10, 10, 9, 9, 4,\n", | |
" 4, 4, 4, 4, 4, 2, 2, 2, 2, 2, 2, 1, 1, 1, 1, 1, 1,\n", | |
" 0, 0]),\n", | |
" array([2. , 2.05, 2.1 , 2.15, 2.2 , 2.25, 2.3 , 2.35, 2.4 , 2.45, 2.5 ,\n", | |
" 2.55, 2.6 , 2.65, 2.7 , 2.75, 2.8 , 2.85, 2.9 , 2.95, 3. , 3.05,\n", | |
" 3.1 , 3.15, 3.2 , 3.25, 3.3 , 3.35, 3.4 , 3.45, 3.5 , 3.55, 3.6 ,\n", | |
" 3.65, 3.7 , 3.75, 3.8 , 3.85, 3.9 , 3.95, 4. , 4.05, 4.1 , 4.15,\n", | |
" 4.2 , 4.25, 4.3 , 4.35, 4.4 , 4.45, 4.5 , 4.55, 4.6 ]))" | |
] | |
}, | |
"metadata": {}, | |
"execution_count": 13 | |
} | |
] | |
}, | |
{ | |
"cell_type": "code", | |
"metadata": { | |
"id": "LazULtNsqXkM" | |
}, | |
"source": [ | |
"M_observation=M_observation_before_2006\n", | |
"beta=np.log(10)*0.2939839569806628\n", | |
"mu_0 = 3.5\n", | |
"std_0=0.38\n", | |
"theta0=[beta,mu_0,std_0]\n", | |
"bounds=((0.3*np.log(10),0.8*np.log(10)),(3.0,5.0),(0.01,0.8))" | |
], | |
"execution_count": null, | |
"outputs": [] | |
}, | |
{ | |
"cell_type": "markdown", | |
"metadata": { | |
"id": "NdHFqO92kByp" | |
}, | |
"source": [ | |
"### trust regeon method\n", | |
"\n", | |
"Conn, A. R., Gould, N. I., and Toint, P. L. Trust region methods. 2000. Siam. pp. 169-200." | |
] | |
}, | |
{ | |
"cell_type": "code", | |
"metadata": { | |
"colab": { | |
"base_uri": "https://localhost:8080/" | |
}, | |
"id": "YdG6Pd_ihcOK", | |
"outputId": "1614f102-5760-489b-850e-c5450d2c8722" | |
}, | |
"source": [ | |
"res3 = minimize(log_l, theta0, method='trust-exact', jac=der_log_1, hess=hess_log_1, options={'gtol': 1e-8, 'maxiter': 10000, 'disp': True}, bounds=bounds)\n", | |
"res3.x" | |
], | |
"execution_count": null, | |
"outputs": [ | |
{ | |
"output_type": "stream", | |
"name": "stdout", | |
"text": [ | |
"Warning: Maximum number of iterations has been exceeded.\n", | |
" Current function value: -45.180518\n", | |
" Iterations: 10000\n", | |
" Function evaluations: 2\n", | |
" Gradient evaluations: 1\n", | |
" Hessian evaluations: 2\n" | |
] | |
}, | |
{ | |
"output_type": "execute_result", | |
"data": { | |
"text/plain": [ | |
"array([0.67692308, 3.5 , 0.38 ])" | |
] | |
}, | |
"metadata": {}, | |
"execution_count": 20 | |
} | |
] | |
}, | |
{ | |
"cell_type": "code", | |
"metadata": { | |
"colab": { | |
"base_uri": "https://localhost:8080/" | |
}, | |
"id": "sSlZCXbHlydb", | |
"outputId": "805b3447-e1ef-46ec-8366-02d7db81ea76" | |
}, | |
"source": [ | |
"res3 = minimize(log_l, theta0, method='trust-constr', jac=der_log_1, hess=hess_log_1, options={'xtol': 1e-8, 'verbose': 1, 'maxiter': 10000}, bounds=bounds)\n", | |
"res3.x" | |
], | |
"execution_count": null, | |
"outputs": [ | |
{ | |
"output_type": "stream", | |
"name": "stdout", | |
"text": [ | |
"The maximum number of function evaluations is exceeded.\n", | |
"Number of iterations: 10000, function evaluations: 4, CG iterations: 19996, optimality: 1.14e+01, constraint violation: 0.00e+00, execution time: 1.5e+01 s.\n" | |
] | |
}, | |
{ | |
"output_type": "execute_result", | |
"data": { | |
"text/plain": [ | |
"array([1.52698598, 3.00211599, 0.26584321])" | |
] | |
}, | |
"metadata": {}, | |
"execution_count": 21 | |
} | |
] | |
}, | |
{ | |
"cell_type": "code", | |
"metadata": { | |
"colab": { | |
"base_uri": "https://localhost:8080/" | |
}, | |
"id": "MbCP4HJRXcF7", | |
"outputId": "ca6bc327-0bfd-446a-bd50-843e265ee069" | |
}, | |
"source": [ | |
"log_l([1.52698598, 3.00211599, 0.26584321])" | |
], | |
"execution_count": null, | |
"outputs": [ | |
{ | |
"output_type": "execute_result", | |
"data": { | |
"text/plain": [ | |
"-101.45002121112489" | |
] | |
}, | |
"metadata": {}, | |
"execution_count": 22 | |
} | |
] | |
}, | |
{ | |
"cell_type": "markdown", | |
"metadata": { | |
"id": "AKyazUxskLta" | |
}, | |
"source": [ | |
"### other methods" | |
] | |
}, | |
{ | |
"cell_type": "code", | |
"metadata": { | |
"colab": { | |
"base_uri": "https://localhost:8080/" | |
}, | |
"id": "32CyshzYKrRg", | |
"outputId": "95399347-abee-41c4-c204-b6287017fa7e" | |
}, | |
"source": [ | |
"res = minimize(log_l, theta0, method='BFGS', jac=der_log_1, options={'gtol': 1e-8,'disp': True})\n", | |
"res.x" | |
], | |
"execution_count": null, | |
"outputs": [ | |
{ | |
"output_type": "stream", | |
"name": "stdout", | |
"text": [ | |
"Warning: Desired error not necessarily achieved due to precision loss.\n", | |
" Current function value: nan\n", | |
" Iterations: 9\n", | |
" Function evaluations: 685\n", | |
" Gradient evaluations: 685\n" | |
] | |
}, | |
{ | |
"output_type": "stream", | |
"name": "stderr", | |
"text": [ | |
"/usr/local/lib/python3.7/dist-packages/ipykernel_launcher.py:6: RuntimeWarning: overflow encountered in double_scalars\n", | |
" \n", | |
"/usr/local/lib/python3.7/dist-packages/ipykernel_launcher.py:7: RuntimeWarning: overflow encountered in double_scalars\n", | |
" import sys\n", | |
"/usr/local/lib/python3.7/dist-packages/ipykernel_launcher.py:6: RuntimeWarning: invalid value encountered in log\n", | |
" \n" | |
] | |
}, | |
{ | |
"output_type": "execute_result", | |
"data": { | |
"text/plain": [ | |
"array([-5.64089017e+40, 1.02607809e+41, -5.86270257e+39])" | |
] | |
}, | |
"metadata": {}, | |
"execution_count": 491 | |
} | |
] | |
}, | |
{ | |
"cell_type": "code", | |
"metadata": { | |
"colab": { | |
"base_uri": "https://localhost:8080/" | |
}, | |
"id": "IjnRHYl_dDxn", | |
"outputId": "d3661e76-cb40-47d5-c96c-38f2df71959a" | |
}, | |
"source": [ | |
"res3 = minimize(log_l, theta0, method='Newton-CG', jac=der_log_1, hess=hess_log_1, options={'xtol': 1e-8, 'disp': True})\n", | |
"res3.x" | |
], | |
"execution_count": null, | |
"outputs": [ | |
{ | |
"output_type": "stream", | |
"name": "stdout", | |
"text": [ | |
"Warning: CG iterations didn't converge. The Hessian is not positive definite.\n", | |
" Current function value: nan\n", | |
" Iterations: 2\n", | |
" Function evaluations: 31\n", | |
" Gradient evaluations: 27\n", | |
" Hessian evaluations: 3\n" | |
] | |
}, | |
{ | |
"output_type": "execute_result", | |
"data": { | |
"text/plain": [ | |
"array([ 115.76497681, -9.09424982, -458.53593901])" | |
] | |
}, | |
"metadata": {}, | |
"execution_count": 492 | |
} | |
] | |
}, | |
{ | |
"cell_type": "code", | |
"metadata": { | |
"colab": { | |
"base_uri": "https://localhost:8080/" | |
}, | |
"id": "6hNju2ejgEkt", | |
"outputId": "a3494ab7-a501-4e9f-f80a-cc8f90bfb26f" | |
}, | |
"source": [ | |
"res3 = minimize(log_l, theta0, method='trust-ncg', jac=der_log_1, hess=hess_log_1, options={'gtol': 1e-8, 'maxiter': 10000, 'disp': True}, bounds=bounds)\n", | |
"res3.x" | |
], | |
"execution_count": null, | |
"outputs": [ | |
{ | |
"output_type": "stream", | |
"name": "stderr", | |
"text": [ | |
"/usr/local/lib/python3.7/dist-packages/scipy/optimize/_minimize.py:539: RuntimeWarning: Method trust-ncg cannot handle constraints nor bounds.\n", | |
" RuntimeWarning)\n" | |
] | |
}, | |
{ | |
"output_type": "stream", | |
"name": "stdout", | |
"text": [ | |
"Warning: Maximum number of iterations has been exceeded.\n", | |
" Current function value: -83.749301\n", | |
" Iterations: 10000\n", | |
" Function evaluations: 3\n", | |
" Gradient evaluations: 2\n", | |
" Hessian evaluations: 2\n" | |
] | |
}, | |
{ | |
"output_type": "execute_result", | |
"data": { | |
"text/plain": [ | |
"array([1.1584505 , 2.62513242, 0.43232642])" | |
] | |
}, | |
"metadata": {}, | |
"execution_count": 493 | |
} | |
] | |
}, | |
{ | |
"cell_type": "code", | |
"metadata": { | |
"colab": { | |
"base_uri": "https://localhost:8080/" | |
}, | |
"id": "AIC2JKPOMP0r", | |
"outputId": "0b58c14c-41b5-4354-ee23-b69f645cc7c3" | |
}, | |
"source": [ | |
"res2 = minimize(log_l, theta0, method='nelder-mead', bounds=bounds, options={ 'xatol': 1e-8,'disp': True, 'maxiter': 10000})\n", | |
"res2.x" | |
], | |
"execution_count": null, | |
"outputs": [ | |
{ | |
"output_type": "stream", | |
"name": "stdout", | |
"text": [ | |
"Optimization terminated successfully.\n", | |
" Current function value: -123.648179\n", | |
" Iterations: 20\n", | |
" Function evaluations: 35\n" | |
] | |
}, | |
{ | |
"output_type": "stream", | |
"name": "stderr", | |
"text": [ | |
"/usr/local/lib/python3.7/dist-packages/ipykernel_launcher.py:1: OptimizeWarning: Initial guess is not within the specified bounds\n", | |
" \"\"\"Entry point for launching an IPython kernel.\n" | |
] | |
}, | |
{ | |
"output_type": "execute_result", | |
"data": { | |
"text/plain": [ | |
"array([1.84206807, 3. , 0.01 ])" | |
] | |
}, | |
"metadata": {}, | |
"execution_count": 497 | |
} | |
] | |
}, | |
{ | |
"cell_type": "markdown", | |
"metadata": { | |
"id": "ZFlSol8DtDDQ" | |
}, | |
"source": [ | |
"## after 2006" | |
] | |
}, | |
{ | |
"cell_type": "code", | |
"metadata": { | |
"colab": { | |
"base_uri": "https://localhost:8080/" | |
}, | |
"id": "aAInXcUJVHXl", | |
"outputId": "9c41154f-7577-4785-e861-d2a3fb1ea1fe" | |
}, | |
"source": [ | |
"M_observation_after_2006 = catalog[catalog['date'] > '20060917']['ML']\n", | |
"print(M_observation_after_2006.describe())\n", | |
"M_observation_after_2006 = M_observation_after_2006.to_numpy()" | |
], | |
"execution_count": null, | |
"outputs": [ | |
{ | |
"output_type": "stream", | |
"name": "stdout", | |
"text": [ | |
"count 415.000000\n", | |
"mean 1.784819\n", | |
"std 0.651446\n", | |
"min 0.100000\n", | |
"25% 1.300000\n", | |
"50% 1.700000\n", | |
"75% 2.100000\n", | |
"max 4.500000\n", | |
"Name: ML, dtype: float64\n" | |
] | |
} | |
] | |
}, | |
{ | |
"cell_type": "code", | |
"metadata": { | |
"colab": { | |
"base_uri": "https://localhost:8080/" | |
}, | |
"id": "T0JJ_EjAtABI", | |
"outputId": "1f1a44ac-cd9a-4edb-9fa8-bab375353096" | |
}, | |
"source": [ | |
"calc_recurrence(M_observation_after_2006.tolist(), min_mag = 1.0, max_mag = None, interval = 0.05)" | |
], | |
"execution_count": null, | |
"outputs": [ | |
{ | |
"output_type": "stream", | |
"name": "stdout", | |
"text": [ | |
"Minimum magnitude: 1.0\n", | |
"Total number of earthquakes: 415\n", | |
"Maximum catalog magnitude: 4.5\n", | |
"Mmax = 4.6\n", | |
"Maximum Likelihood: b value 0.5533687749151043\n" | |
] | |
}, | |
{ | |
"output_type": "execute_result", | |
"data": { | |
"text/plain": [ | |
"(array([392, 382, 382, 353, 334, 334, 300, 300, 270, 270, 248, 248, 225,\n", | |
" 225, 199, 199, 165, 165, 144, 144, 124, 124, 100, 100, 85, 85,\n", | |
" 71, 71, 61, 61, 56, 56, 45, 45, 32, 32, 26, 26, 20,\n", | |
" 20, 16, 16, 13, 13, 9, 9, 7, 7, 7, 7, 7, 7,\n", | |
" 5, 5, 5, 5, 5, 5, 4, 4, 2, 2, 2, 2, 2,\n", | |
" 2, 1, 1, 1, 1, 0, 0]),\n", | |
" array([1. , 1.05, 1.1 , 1.15, 1.2 , 1.25, 1.3 , 1.35, 1.4 , 1.45, 1.5 ,\n", | |
" 1.55, 1.6 , 1.65, 1.7 , 1.75, 1.8 , 1.85, 1.9 , 1.95, 2. , 2.05,\n", | |
" 2.1 , 2.15, 2.2 , 2.25, 2.3 , 2.35, 2.4 , 2.45, 2.5 , 2.55, 2.6 ,\n", | |
" 2.65, 2.7 , 2.75, 2.8 , 2.85, 2.9 , 2.95, 3. , 3.05, 3.1 , 3.15,\n", | |
" 3.2 , 3.25, 3.3 , 3.35, 3.4 , 3.45, 3.5 , 3.55, 3.6 , 3.65, 3.7 ,\n", | |
" 3.75, 3.8 , 3.85, 3.9 , 3.95, 4. , 4.05, 4.1 , 4.15, 4.2 , 4.25,\n", | |
" 4.3 , 4.35, 4.4 , 4.45, 4.5 , 4.55]))" | |
] | |
}, | |
"metadata": {}, | |
"execution_count": 24 | |
} | |
] | |
}, | |
{ | |
"cell_type": "code", | |
"metadata": { | |
"id": "u-HKsU56tSNi" | |
}, | |
"source": [ | |
"M_observation=M_observation_after_2006\n", | |
"beta=np.log(10)*0.5533687749151043\n", | |
"mu_0 = 1.78\n", | |
"std_0=0.65\n", | |
"theta0=[beta,mu_0,std_0]\n", | |
"bounds=((0.3*np.log(10),0.8*np.log(10)),(1.0,5.0),(0.01,0.9))" | |
], | |
"execution_count": null, | |
"outputs": [] | |
}, | |
{ | |
"cell_type": "markdown", | |
"metadata": { | |
"id": "C1jgu6b2kvyR" | |
}, | |
"source": [ | |
"### trust regeon method\n", | |
"\n", | |
"Conn, A. R., Gould, N. I., and Toint, P. L. Trust region methods. 2000. Siam. pp. 169-200." | |
] | |
}, | |
{ | |
"cell_type": "code", | |
"metadata": { | |
"colab": { | |
"base_uri": "https://localhost:8080/" | |
}, | |
"id": "r9KD0WDhuFO9", | |
"outputId": "df3ee79c-dbfe-4568-8828-b12d92049682" | |
}, | |
"source": [ | |
"res3 = minimize(log_l, theta0, method='trust-exact', jac=der_log_1, hess=hess_log_1, options={'gtol': 1e-8, 'disp': True}, bounds=bounds)\n", | |
"res3.x" | |
], | |
"execution_count": null, | |
"outputs": [ | |
{ | |
"output_type": "stream", | |
"name": "stdout", | |
"text": [ | |
"Warning: A bad approximation caused failure to predict improvement.\n", | |
" Current function value: -555.888176\n", | |
" Iterations: 41\n", | |
" Function evaluations: 43\n", | |
" Gradient evaluations: 36\n", | |
" Hessian evaluations: 43\n" | |
] | |
}, | |
{ | |
"output_type": "execute_result", | |
"data": { | |
"text/plain": [ | |
"array([1.41156477, 1.41704815, 0.46302158])" | |
] | |
}, | |
"metadata": {}, | |
"execution_count": 26 | |
} | |
] | |
}, | |
{ | |
"cell_type": "code", | |
"metadata": { | |
"colab": { | |
"base_uri": "https://localhost:8080/" | |
}, | |
"id": "HKd8q0jzuJbi", | |
"outputId": "a8331ad6-4c17-4f1d-9b1b-614107678b79" | |
}, | |
"source": [ | |
"res3 = minimize(log_l, theta0, method='trust-constr', jac=der_log_1, hess=hess_log_1, options={'xtol': 1e-8, 'verbose': 1}, bounds=bounds)\n", | |
"res3.x" | |
], | |
"execution_count": null, | |
"outputs": [ | |
{ | |
"output_type": "stream", | |
"name": "stdout", | |
"text": [ | |
"`xtol` termination condition is satisfied.\n", | |
"Number of iterations: 324, function evaluations: 324, CG iterations: 366, optimality: 1.72e+02, constraint violation: 0.00e+00, execution time: 0.61 s.\n" | |
] | |
}, | |
{ | |
"output_type": "execute_result", | |
"data": { | |
"text/plain": [ | |
"array([1.41298784, 1.41457276, 0.38406817])" | |
] | |
}, | |
"metadata": {}, | |
"execution_count": 27 | |
} | |
] | |
}, | |
{ | |
"cell_type": "code", | |
"metadata": { | |
"colab": { | |
"base_uri": "https://localhost:8080/" | |
}, | |
"id": "e8fCxNZVYj0Z", | |
"outputId": "6012d8fe-6b66-4c10-d89a-be3689e2adf6" | |
}, | |
"source": [ | |
"log_l([1.41298784, 1.41457276, 0.38406817])" | |
], | |
"execution_count": null, | |
"outputs": [ | |
{ | |
"output_type": "execute_result", | |
"data": { | |
"text/plain": [ | |
"-529.1446723322614" | |
] | |
}, | |
"metadata": {}, | |
"execution_count": 28 | |
} | |
] | |
}, | |
{ | |
"cell_type": "markdown", | |
"metadata": { | |
"id": "4ljvHeeJlNip" | |
}, | |
"source": [ | |
"### other methods" | |
] | |
}, | |
{ | |
"cell_type": "code", | |
"metadata": { | |
"colab": { | |
"base_uri": "https://localhost:8080/" | |
}, | |
"id": "aPZ8-lrQt3FS", | |
"outputId": "d61ab2e9-d15d-4a6d-fede-2610e195e4c3" | |
}, | |
"source": [ | |
"res = minimize(log_l, theta0, method='BFGS', jac=der_log_1, options={'gtol': 1e-8,'disp': True})\n", | |
"res.x" | |
], | |
"execution_count": null, | |
"outputs": [ | |
{ | |
"output_type": "stream", | |
"name": "stdout", | |
"text": [ | |
"Warning: Desired error not necessarily achieved due to precision loss.\n", | |
" Current function value: -333.101804\n", | |
" Iterations: 0\n", | |
" Function evaluations: 7\n", | |
" Gradient evaluations: 4\n" | |
] | |
}, | |
{ | |
"output_type": "execute_result", | |
"data": { | |
"text/plain": [ | |
"array([1.27417869, 1.78 , 0.65 ])" | |
] | |
}, | |
"metadata": {}, | |
"execution_count": 483 | |
} | |
] | |
}, | |
{ | |
"cell_type": "code", | |
"metadata": { | |
"colab": { | |
"base_uri": "https://localhost:8080/" | |
}, | |
"id": "nsfPTxPGt_kE", | |
"outputId": "e3fbcf18-e3d4-4245-b42d-509c627cfcea" | |
}, | |
"source": [ | |
"res3 = minimize(log_l, theta0, method='Newton-CG', jac=der_log_1, hess=hess_log_1, options={'xtol': 1e-8, 'disp': True})\n", | |
"res3.x" | |
], | |
"execution_count": null, | |
"outputs": [ | |
{ | |
"output_type": "stream", | |
"name": "stdout", | |
"text": [ | |
"Warning: CG iterations didn't converge. The Hessian is not positive definite.\n", | |
" Current function value: nan\n", | |
" Iterations: 1\n", | |
" Function evaluations: 16\n", | |
" Gradient evaluations: 16\n", | |
" Hessian evaluations: 2\n" | |
] | |
}, | |
{ | |
"output_type": "execute_result", | |
"data": { | |
"text/plain": [ | |
"array([ 337.84035575, -18697.41039866, -137.99089267])" | |
] | |
}, | |
"metadata": {}, | |
"execution_count": 484 | |
} | |
] | |
}, | |
{ | |
"cell_type": "code", | |
"metadata": { | |
"colab": { | |
"base_uri": "https://localhost:8080/" | |
}, | |
"id": "DK5R-U1huBIL", | |
"outputId": "a3db9fbc-5e36-4466-c418-ace74d084dae" | |
}, | |
"source": [ | |
"res3 = minimize(log_l, theta0, method='trust-ncg', jac=der_log_1, hess=hess_log_1, options={'gtol': 1e-8, 'disp': True}, bounds=bounds)\n", | |
"res3.x" | |
], | |
"execution_count": null, | |
"outputs": [ | |
{ | |
"output_type": "stream", | |
"name": "stdout", | |
"text": [ | |
"Warning: A bad approximation caused failure to predict improvement.\n", | |
" Current function value: -578.053951\n", | |
" Iterations: 94\n", | |
" Function evaluations: 96\n", | |
" Gradient evaluations: 88\n", | |
" Hessian evaluations: 88\n" | |
] | |
}, | |
{ | |
"output_type": "stream", | |
"name": "stderr", | |
"text": [ | |
"/usr/local/lib/python3.7/dist-packages/scipy/optimize/_minimize.py:539: RuntimeWarning: Method trust-ncg cannot handle constraints nor bounds.\n", | |
" RuntimeWarning)\n" | |
] | |
}, | |
{ | |
"output_type": "execute_result", | |
"data": { | |
"text/plain": [ | |
"array([1.42948191, 1.39936953, 0.48862331])" | |
] | |
}, | |
"metadata": {}, | |
"execution_count": 485 | |
} | |
] | |
}, | |
{ | |
"cell_type": "code", | |
"metadata": { | |
"colab": { | |
"base_uri": "https://localhost:8080/" | |
}, | |
"id": "pvlYvHUluUPO", | |
"outputId": "3f0dbbaa-c9c1-40f7-aa20-829a1fceab09" | |
}, | |
"source": [ | |
"res2 = minimize(log_l, theta0, method='nelder-mead', bounds=bounds, options={ 'xatol': 1e-8,'disp': True, 'maxiter': 10000})\n", | |
"res2.x" | |
], | |
"execution_count": null, | |
"outputs": [ | |
{ | |
"output_type": "stream", | |
"name": "stdout", | |
"text": [ | |
"Optimization terminated successfully.\n", | |
" Current function value: -908.755727\n", | |
" Iterations: 220\n", | |
" Function evaluations: 398\n" | |
] | |
}, | |
{ | |
"output_type": "execute_result", | |
"data": { | |
"text/plain": [ | |
"array([1.84206807, 1.08603177, 0.36621488])" | |
] | |
}, | |
"metadata": {}, | |
"execution_count": 443 | |
} | |
] | |
}, | |
{ | |
"cell_type": "code", | |
"metadata": { | |
"colab": { | |
"base_uri": "https://localhost:8080/" | |
}, | |
"id": "e3x3voNxk8SE", | |
"outputId": "750380fe-91a4-4402-b8c9-1765444274da" | |
}, | |
"source": [ | |
"1/hess_log_1([1.37903484, 1.45075722, 0.49417915])" | |
], | |
"execution_count": null, | |
"outputs": [ | |
{ | |
"output_type": "execute_result", | |
"data": { | |
"text/plain": [ | |
"array([[ 3.30381416e-06, 3.46633159e-06, 1.33290164e-03],\n", | |
" [ 3.45549010e-06, 3.65363392e-06, -4.05239864e-04],\n", | |
" [ 2.69720331e-03, -8.20026227e-04, -9.66546749e-04]])" | |
] | |
}, | |
"metadata": {}, | |
"execution_count": 118 | |
} | |
] | |
}, | |
{ | |
"cell_type": "code", | |
"metadata": { | |
"colab": { | |
"base_uri": "https://localhost:8080/" | |
}, | |
"id": "DNCJdir9x-6y", | |
"outputId": "b1299595-6ed2-4626-f6c9-8a6459e85c68" | |
}, | |
"source": [ | |
"1.37643205/np.log(10)" | |
], | |
"execution_count": null, | |
"outputs": [ | |
{ | |
"output_type": "execute_result", | |
"data": { | |
"text/plain": [ | |
"0.5977768440297808" | |
] | |
}, | |
"metadata": {}, | |
"execution_count": 109 | |
} | |
] | |
}, | |
{ | |
"cell_type": "markdown", | |
"metadata": { | |
"id": "Ael13lUsiFc1" | |
}, | |
"source": [ | |
"## ≥ M2.8" | |
] | |
}, | |
{ | |
"cell_type": "code", | |
"metadata": { | |
"colab": { | |
"base_uri": "https://localhost:8080/" | |
}, | |
"id": "hgceJvTphh8p", | |
"outputId": "fbd5ddc6-877a-41cc-f076-c64d26a7ee9a" | |
}, | |
"source": [ | |
"M_observation_M3 = catalog[catalog['ML'] >= 3.0]['ML']\n", | |
"print(M_observation_M3.describe())\n", | |
"M_observation_M3 = M_observation_M3.to_numpy()" | |
], | |
"execution_count": null, | |
"outputs": [ | |
{ | |
"output_type": "stream", | |
"name": "stdout", | |
"text": [ | |
"count 42.000000\n", | |
"mean 3.454762\n", | |
"std 0.422653\n", | |
"min 3.000000\n", | |
"25% 3.100000\n", | |
"50% 3.300000\n", | |
"75% 3.600000\n", | |
"max 4.500000\n", | |
"Name: ML, dtype: float64\n" | |
] | |
} | |
] | |
}, | |
{ | |
"cell_type": "code", | |
"metadata": { | |
"colab": { | |
"base_uri": "https://localhost:8080/" | |
}, | |
"id": "Wa1uZRyriQzE", | |
"outputId": "63ddfa28-d332-4109-e9bb-d13c541dc0f3" | |
}, | |
"source": [ | |
"calc_recurrence(M_observation_M3.tolist(), min_mag = 2.0, max_mag = None, interval = 0.05)" | |
], | |
"execution_count": null, | |
"outputs": [ | |
{ | |
"output_type": "stream", | |
"name": "stdout", | |
"text": [ | |
"Minimum magnitude: 2.0\n", | |
"Total number of earthquakes: 42\n", | |
"Maximum catalog magnitude: 4.5\n", | |
"Mmax = 4.6\n", | |
"Maximum Likelihood: b value 0.2985330317501895\n" | |
] | |
}, | |
{ | |
"output_type": "execute_result", | |
"data": { | |
"text/plain": [ | |
"(array([42, 42, 42, 42, 42, 42, 42, 42, 42, 42, 42, 42, 42, 42, 42, 42, 42,\n", | |
" 42, 42, 42, 42, 35, 35, 30, 30, 25, 25, 20, 20, 17, 17, 16, 16, 9,\n", | |
" 9, 9, 9, 9, 9, 6, 6, 4, 4, 4, 4, 3, 3, 2, 2, 2, 2,\n", | |
" 0, 0]),\n", | |
" array([2. , 2.05, 2.1 , 2.15, 2.2 , 2.25, 2.3 , 2.35, 2.4 , 2.45, 2.5 ,\n", | |
" 2.55, 2.6 , 2.65, 2.7 , 2.75, 2.8 , 2.85, 2.9 , 2.95, 3. , 3.05,\n", | |
" 3.1 , 3.15, 3.2 , 3.25, 3.3 , 3.35, 3.4 , 3.45, 3.5 , 3.55, 3.6 ,\n", | |
" 3.65, 3.7 , 3.75, 3.8 , 3.85, 3.9 , 3.95, 4. , 4.05, 4.1 , 4.15,\n", | |
" 4.2 , 4.25, 4.3 , 4.35, 4.4 , 4.45, 4.5 , 4.55, 4.6 ]))" | |
] | |
}, | |
"metadata": {}, | |
"execution_count": 499 | |
} | |
] | |
}, | |
{ | |
"cell_type": "code", | |
"metadata": { | |
"id": "yWvIbfbOiwtA" | |
}, | |
"source": [ | |
"M_observation=M_observation_M3\n", | |
"beta=np.log(10)*0.2985330317501895\n", | |
"mu_0 = 3.4\n", | |
"std_0=0.42\n", | |
"theta0=[beta,mu_0,std_0]\n", | |
"bounds=((0.1*np.log(10),0.8*np.log(10)),(3.0,5.0),(0.01,0.8))" | |
], | |
"execution_count": null, | |
"outputs": [] | |
}, | |
{ | |
"cell_type": "code", | |
"metadata": { | |
"colab": { | |
"base_uri": "https://localhost:8080/" | |
}, | |
"id": "7aHPmUink8Ff", | |
"outputId": "83efc994-719a-4930-a628-2b104ca16514" | |
}, | |
"source": [ | |
"res3 = minimize(log_l, theta0, method='trust-exact', jac=der_log_1, hess=hess_log_1, options={'gtol': 1e-8, 'disp': True}, bounds=bounds)\n", | |
"res3.x" | |
], | |
"execution_count": null, | |
"outputs": [ | |
{ | |
"output_type": "stream", | |
"name": "stdout", | |
"text": [ | |
"Warning: Maximum number of iterations has been exceeded.\n", | |
" Current function value: -90.958199\n", | |
" Iterations: 600\n", | |
" Function evaluations: 2\n", | |
" Gradient evaluations: 1\n", | |
" Hessian evaluations: 2\n" | |
] | |
}, | |
{ | |
"output_type": "execute_result", | |
"data": { | |
"text/plain": [ | |
"array([0.68739771, 3.4 , 0.42 ])" | |
] | |
}, | |
"metadata": {}, | |
"execution_count": 625 | |
} | |
] | |
}, | |
{ | |
"cell_type": "code", | |
"metadata": { | |
"colab": { | |
"base_uri": "https://localhost:8080/" | |
}, | |
"id": "kRSCG-yckRfo", | |
"outputId": "af421d9a-4af6-4377-99bf-aebb3a0d913d" | |
}, | |
"source": [ | |
"res3 = minimize(log_l, theta0, method='trust-constr', jac=der_log_1, hess=hess_log_1, options={'xtol': 1e-10, 'verbose': 1}, bounds=bounds)\n", | |
"res3.x" | |
], | |
"execution_count": null, | |
"outputs": [ | |
{ | |
"output_type": "stream", | |
"name": "stdout", | |
"text": [ | |
"The maximum number of function evaluations is exceeded.\n", | |
"Number of iterations: 1000, function evaluations: 4, CG iterations: 1996, optimality: 2.30e+01, constraint violation: 0.00e+00, execution time: 1.4 s.\n" | |
] | |
}, | |
{ | |
"output_type": "execute_result", | |
"data": { | |
"text/plain": [ | |
"array([1.56896081, 3.00187221, 0.22273732])" | |
] | |
}, | |
"metadata": {}, | |
"execution_count": 626 | |
} | |
] | |
}, | |
{ | |
"cell_type": "code", | |
"metadata": { | |
"colab": { | |
"base_uri": "https://localhost:8080/" | |
}, | |
"id": "UBhMiFYKly4F", | |
"outputId": "2b5453d3-d6d7-4964-e8d0-a80b834f0629" | |
}, | |
"source": [ | |
"log_l([1.56896081, 3.00187221, 0.22273732])" | |
], | |
"execution_count": null, | |
"outputs": [ | |
{ | |
"output_type": "execute_result", | |
"data": { | |
"text/plain": [ | |
"-195.31041653986026" | |
] | |
}, | |
"metadata": {}, | |
"execution_count": 507 | |
} | |
] | |
}, | |
{ | |
"cell_type": "markdown", | |
"metadata": { | |
"id": "PWT4SgJYrmD8" | |
}, | |
"source": [ | |
"# Plot results" | |
] | |
}, | |
{ | |
"cell_type": "code", | |
"metadata": { | |
"colab": { | |
"base_uri": "https://localhost:8080/", | |
"height": 578 | |
}, | |
"id": "8gbkC0hP3XlL", | |
"outputId": "ac7ea829-12de-44ae-a487-ca9151366039" | |
}, | |
"source": [ | |
"plt.rcParams.update(plt.rcParamsDefault)\n", | |
"plt.rcParams['figure.figsize'] = [10, 6]\n", | |
"plt.rcParams['ytick.major.size'] = 10\n", | |
"#plt.rcParams['ytick.major.width'] = 4\n", | |
"plt.rcParams['ytick.minor.size'] = 5\n", | |
"#plt.rcParams['xtick.minor.width'] = 2\n", | |
"\n", | |
"plt.subplot(121)\n", | |
"plt.title(\"Prior 2006-09\",fontsize = 15)\n", | |
"plt.ylabel('Probability', fontsize = 15)\n", | |
"plt.xlabel('Magnitude', fontsize = 15)\n", | |
"plt.yscale('log')\n", | |
"plt.yticks(fontsize=14)\n", | |
"plt.xticks(fontsize=15)\n", | |
"plt.ylim([0.001, 1])\n", | |
"plt.xlim([0.0, 5.0])\n", | |
"plt.text(0.4, 0.5,r'$b=0.7$'+'\\n'+r'$\\mu=3.0$'+'\\n'+r'$\\sigma=0.3$', ha='left', va='center', fontsize = 15)\n", | |
"# ECDF\n", | |
"x = np.sort(M_observation_before_2006)\n", | |
"y = np.arange(len(x))/float(len(x))\n", | |
"# EPMF\n", | |
"y[1:] -= y[:-1].copy()\n", | |
"u, inv = np.unique(x, return_inverse=True)\n", | |
"sums = np.zeros(len(u), dtype=y.dtype) \n", | |
"np.add.at(sums, inv, y)\n", | |
"plt.scatter(u, sums, marker='+')\n", | |
"# estimated density\n", | |
"x3=np.linspace(0,5,len(x))\n", | |
"plt.plot(x2,density(x2,[1.52698598, 3.00211599, 0.26584321]))\n", | |
"\n", | |
"plt.subplot(122)\n", | |
"plt.title(\"After 2006-09\", fontsize = 15)\n", | |
"#plt.ylabel('Probability', fontsize = 15)\n", | |
"#plt.xlabel('Magnitude', fontsize = 15)\n", | |
"plt.yscale('log')\n", | |
"plt.yticks([])\n", | |
"plt.xticks(fontsize=15)\n", | |
"plt.ylim([0.001, 1])\n", | |
"plt.xlim([0.0, 5.0])\n", | |
"plt.text(0.4, 0.5,r'$b=0.6$'+'\\n'+r'$\\mu=1.4$'+'\\n'+r'$\\sigma=0.4$', ha='left', va='center', fontsize = 15)\n", | |
"# ECDF\n", | |
"x = np.sort(M_observation_after_2006)\n", | |
"y = np.arange(len(x))/float(len(x))\n", | |
"# EPMF\n", | |
"y[1:] -= y[:-1].copy()\n", | |
"u, inv = np.unique(x, return_inverse=True)\n", | |
"sums = np.zeros(len(u), dtype=y.dtype) \n", | |
"np.add.at(sums, inv, y)\n", | |
"plt.scatter(u, sums, marker='+')\n", | |
"# estimated density\n", | |
"x2 = np.arange(0,5,0.1)\n", | |
"plt.plot(x2,density(x2,[1.41298784, 1.41457276, 0.38406817]))\n", | |
"#plt.savefig('fig_1.png',dpi=600)\n", | |
"plt.show()\n" | |
], | |
"execution_count": null, | |
"outputs": [ | |
{ | |
"output_type": "display_data", | |
"data": { | |
"image/png": "iVBORw0KGgoAAAANSUhEUgAAA3sAAAIxCAYAAAD0eaFdAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAAPYQAAD2EBqD+naQAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAgAElEQVR4nOzdeVhTV/oH8G9YEkA2BSwoArKpdUXUahWDrSO4FbUqThdFxRmn2lrt8rPWZbR2bEfR6VjajrVV61atWpdqtdWKjljcEHBBUSrgLigEWQPh/P5gkpomIIFAMHw/z8Mzcs5d3pvpc1/e3HPukQghBIiIiIiIiMisWJg6ACIiIiIiIjI+FntERERERERmiMUeERERERGRGWKxR0REREREZIZY7BEREREREZkhFntERERERERmiMUeERERERGRGWKxR0REREREZIZY7BEREREREZkhFntEtSSRSLR+LCws4OzsjJCQEKxevRpCCIOOFxcXB4lEgqioqPoJuIYyMzOxcuVKhIeHw93dHdbW1nB1dUV4eDh2795d7b65ubmYMWMGvL29IZPJ4O3tjTfffBN5eXlV7qNSqbBixQp07twZtra2cHNzw9ixY5GamlrtuYQQWLt2Lfr3748WLVrA1tYWvr6+eOmll3DhwoVaXfvatWvRq1cv2Nvbo0WLFhgyZAiOHz9ebewrV65EcHAwmjVrBicnJ/Tv3x87duyo1fmJiIzl5MmTmvy0aNGiardNT0/HyJEj4erqCgsLC0gkEsTFxTVMoAa4e/cuvvrqK4wcORKenp6QSqVwdnaGXC7HunXrqs27xcXFmD9/PgIDA2FjY4NWrVph0qRJuHnzZrXnNDQvqO3cuRPh4eFwc3ODjY0N2rRpg5EjR+LYsWMGXzcA7NmzB3K5HI6OjnB0dERoaCj27t1b7T4bN25E37594eDgAHt7e/Ts2RNffvmlwX+f0BNOEFGtABAAxIQJE8SECRPEK6+8Ivr06SMkEokAIMaNG2fQ8Q4fPqw5nin17dtXABAymUzI5XIRGRkpevbsqbnemTNn6t0vOztb+Pv7CwDC19dXjB07VnTs2FEAEIGBgeL+/fs6+6hUKjFy5EgBQDg7O4sXX3xRyOVyIZFIhJ2dnThx4oTecxUXF4vw8HABQLRo0UIMGzZMjBkzRgQHBwtLS0uxfv16g697xowZAoCwtbUVERERIiwsTFhZWQlLS0vx/fff62xfXl4uhg0bJgAIe3t7MWjQIDFw4EBhZ2cnAIgFCxYYHAMRkbFMnz5dc98ODAyscjuVSiW6desmAIjevXuLV199VUyYMEGkpqaKCRMmCADi8OHDDRd4NV5++WUBQFhZWYnevXuLyMhI0a9fP2FhYSEAiNGjR4vy8nKd/YqLi0Xv3r0FAOHh4SHGjh0revXqJQAINzc3kZ6ervd8huYFISo/z0mTJgkAolmzZiIsLExERkaKPn36CKlUKj744AODr3vFihWa6w4PDxcRERHC1tZWABArV67Uu8/UqVMFACGVSoVcLhdDhgwRzs7OjeLvDGpYLPaIakmdRP/op59+ElZWVgKA2LNnT42PV1hYKFJTU8WtW7eMGabBIiMjxcqVK0V+fr5W+w8//KC5rgMHDujsp07Co0aNEmVlZZr2119/vcrk8uWXXwoAIiAgQNy5c0fTvm3bNgFA+Pv7ax1LTf0HyJQpU0RRUZFW361bt0RmZqZB1/zzzz8LAMLFxUWkpaVp2o8fPy6kUqlwdnYWubm5WvssW7ZMABA+Pj7i6tWrmvbU1FTRqlUrAUAcP37coDiIiIxBqVQKV1dXAUC4u7sLACIhIUHvtunp6QKACAkJ0elrbMXeG2+8IT788ENx7949rfaTJ08KR0dHAUD85z//0dnv/fffFwBEnz59xMOHDzXtMTExAoCQy+U6+9QmLwghxIIFCwQAMXz4cJ0vOR88eKB1rJq4dOmSsLS0FDKZTCunXL58Wbi4uAgrKytx5coVrX3UObR58+bi9OnTmvZbt26JTp06CQBi06ZNBsVBTy4We0S1VFWxJ4QQEydOFADE5MmTGziq+vWXv/xFABBRUVFa7bdu3RIWFhZCKpVqFW1CCFFSUiLc3NyEpaWluHv3rlZfhw4dBAC935C+8MILAoDYtm2bVvuJEycEANGrVy9RUVFhlOsaPHiwACBWrFih0/fGG28IAGLZsmVa7X5+fgKA2Lhxo84+q1atEgDEiBEjjBIfEZEhdu/eLQCIvn37ikWLFgkAYtq0aXq3PXLkSJVfyDW2Yq86//jHPwQAERoaqtVeWloqnJycBACRmJios1+XLl0EAK2iSIja5YXr168LqVQqvLy8dL6IrK2//e1vAoCYMWOGTt/y5csFADF9+nSt9ueff14AEB9++KHOPj/99JMAILp162aU+Kjx45w9onoQFBQEALh+/bqmTSKRwMfHB0qlEosWLUL79u0hk8kwYsQIANXP2SsvL9fMDbO3t4e9vT169eqFzz//HCqVSmf70NBQSCQSZGRkYNOmTejduzccHBzg7Oxcp+vq2rUrAODWrVta7fv370dFRQVCQkLw1FNPafXJZDIMHz4cKpUK+/bt07Rfu3YNqampsLW1xdChQ3XONXr0aACV8xQe9eWXXwIApk+fDolEUqfrASrncfzyyy9a53xcHAqFAunp6QAqP+s/GjBgAADgwIEDKC0trXOMRESG2LBhAwDglVdewSuvvAIA2LJlC8rKyrS2k0gkkMvlAIB169Zp5vipc8i6desAVN7THp2jnpGRoXWc/fv3Y+jQoXBzc4NMJoOvry9mzZqF+/fv68QWFRWlmRN44MABDBgwAM7OzpBIJNXO736cqvJTfHw8FAoF/Pz8NLn5Ufru8bXJC0DlZ6hUKhEdHQ1bW9taX8uj1PPyDInjzJkzAPTnJ7lcDgsLCyQlJSErK8soMVLjZmXqAIjM0cOHDwFUFjqPqqiowIgRI3D06FHI5XJ06dIFLi4u1R5LpVIhIiIC+/btg6OjI/70pz9BCIFffvkFr732Gn7++Wds27YNFha6390sWbIEq1evRt++fTFs2DCt4rM2fvvtNwCAu7u7VntycjIAoHv37nr36969O77++mukpKTo7NOpUydYW1vr3QeA1j4ANAn42WefRXp6OjZv3ozr16/Dzc0N4eHh6Nevn0HXdPnyZZSWlsLNzQ2enp41iqOwsFDz7+bNm+vso/7/tLi4GGlpaejcubNBMRER1ZZCocDu3bshlUoxduxYtGjRAs8++yyOHz+O/fv3Y/jw4ZptJ0yYgDt37uDAgQPw8/PT3D/bt28PHx8fHDt2DOnp6QgLC9O679vb22v+PXv2bHz88ceQSqXo2bMnPDw8kJycjBUrVmD37t2Ij4/X+RIQADZt2oTVq1ejR48eGDx4MNLT0+v0BV5d8hOgfY+vTV4AtPPT7du3sXHjRly9ehVOTk4YMGAAwsLCDLrGvLw8TUGmr1Bt06YNXF1dkZmZifz8fDg6OgL4PUfpy09SqRT29vbIz89HcnIyvLy8ahwPPaFM/WiR6EmFKoZxVlRUiD59+ggA4v3339fZ3t/fX9y4cUNnv6pe0KKeG9axY0etIZK3bt0S7dq10ztBWy6XCwDCxsZGxMXF1fFKK+Xm5go3NzcBQGzfvl2rT/2SlU8++UTvvjt37tTM51P75JNPBAAxcuRIvfvk5eVpXsCiVlxcrPkcV61aJWQymeZ39U9kZKQoLS2t8XXt2rVLABBBQUFVbqOe1K6ex1hcXCwsLS0FAJGamqqz/ZkzZzTxGDJvk4iorlavXi0AiIiICE3bZ599JgCIMWPG6Gxf3cvBHjeMc+vWrQKA6NSpk9a8sYqKCjF//nzNPVnfMQGIb7/9tnYX+QdKpVIzLSAmJkarb+bMmdW+XCwpKUkAEN27d9e01SYvCCE08yP//e9/a4aOPvoTGhqqd55fVZKTkzVz76qifrlOSkqKpk09b/zHH3/U2f7+/fuaeKp6uQuZFw7jJDISlUqFK1euYNKkSfj1118hk8kwceJEne2WLFmC1q1b1/i4//73vwEAy5cv1/p21MPDA0uXLgUAfPLJJ3r3nTx5smaITl1NnToV2dnZ6N27N0aOHKnVV1BQAACws7PTu2+zZs0A/P7Es7b7PDrE57XXXsPQoUORmpqKvLw87NixA66urtiyZQvef//9Gl/X4+LQF4uNjQ169uwJoPK13H/09ddfa/79aPxERPVt/fr1AKAZvgkAY8eOhbW1Nfbs2QOFQmG0c3344YcAgM2bN8Pf31/TLpFI8Pe//x3dunXDtm3bkJOTo7Pv0KFDERkZaZQ45s2bh9TUVLRt2xZTp07V6quP/FTVfrm5uQCAWbNmoUuXLkhMTER+fj4OHjyItm3bIi4uDlOmTKnxddU2jv79+wNgfqJKLPaI6kg9h8HKygqBgYFYu3YtHBwcsHnzZvj5+els++gQmsfJyspCVlYW3NzcMGjQIJ3+YcOGwdnZGVevXsWdO3d0+l944QXDL0iPjz/+GFu2bEGLFi2wceNGo8yVq42KigrNv9u3b4/vvvsO7du3h5OTE0aOHKmZX/Lpp58iPz+/XmOZPXs2ACAmJgbLli3DnTt3cOvWLXz44Yf44osvYGVVOUpe3/BaIqL6kJWVhaNHj8LZ2Vkr17i4uGDIkCEoKSnBd999Z5Rz3bt3D8nJyQgICECnTp10+iUSCfr27QuVSqWZQ/YoY+Wnb7/9Fv/85z9hY2ODTZs2VVsY1Td1jmrevDl+/PFHBAUFwcHBAc8//zx2794NiUSCbdu2IS0trV7jePvtt2FlZYUtW7bg3XffRVZWFnJycvCf//wH8+fPZ35qYjhnj6iOJkyYAKDypuno6IjOnTtj1KhResfKt2zZUmceX3XUE829vb319kskEnh7eyMvLw83b97UmatgjLH4GzZswHvvvYdmzZph79698PX11dlGPX+jqKhI7zHU8wccHByMsg8AjB8/XidRDRkyBC1btsS9e/dw8uRJDBw4EJcuXcJHH32kc/zZs2ejffv2j42jqlgiIiLw8ccfY86cOXjnnXfwzjvvaPr+8pe/4OzZszh16pTe/w6IiOrDxo0bIYTA6NGjdXLNK6+8gl27dmHDhg2Ijo6u87nUL2m5cuXKY78A1Pdkzxj56ZdffkFUVBQsLCywefNm9O7dW2eb+shP1e2Xm5uLMWPGaJ64qXXq1Ak9e/bEyZMncfToUQQGBiInJwdvv/22zrGjo6PRr1+/WscRHByMNWvWYMqUKVi6dKlmFBBQ+UTV2toaO3fuZH5qIljsEdWRvmESVbGxsTH6+atLsnU93w8//ICJEyfC2toaO3bs0JtIgd+T9o0bN/T2q9sfLVprs4+joyOaN2+O3Nxc+Pj46N3Px8cH9+7dw7179wAAd+7c0Tzxe1RUVBTat2//2DgKCwuRl5eH5s2bayVTAHj33XcxcuRIbNu2DRkZGXBycsLQoUMhl8s1k/o7duyo97hERMamHsIZFxen87IqpVIJADh69CgyMzOr/BKxptRPsdzd3REWFlbttvrOVdf8dOrUKURERECpVOKrr77SvNn6j+ojP1WVF7y9vR+bn06ePKnJTwUFBXrzU2hoKPr166eJIzc3F4WFhToFZFXxA5XF/YABA7B161akpaXBxsYGzz//PIYOHYqQkBAAzE9NBYs9okasVatWAIDMzMwqt1H3GTIPsCaOHDmCMWPGQAiBTZs26R1GqqZ+5XViYqLefnV7ly5ddPY5f/48ysrKdN7IqW8fAOjWrRsOHz6smRvxRw8ePADw+zezoaGhEEJUGXu7du0gk8mQnZ2Nmzdv6nyOVcWhFhAQgPfee0+rLSsrCzdv3oS/v7/R/38hItLnzJkzSE1NBQBcvXoVV69e1budEAIbN27EnDlz6nQ+9Rdarq6uBn3paQwXL17E4MGDUVBQgBUrVuidH69Wm/xU27wQFBSEpKSkGucnHx+favOTs7MzvLy8kJWVhbNnz+oU8NevX0dOTg68vb01b+J8VOvWrTFz5kyttuLiYiQlJcHBwaHKN5SSeeFgXaJGzMvLC15eXsjOzsahQ4d0+vfu3Yvc3Fz4+/vrDOGsi8TERLzwwgsoLS3F6tWr8eKLL1a7fXh4OCwsLPDf//5X842lWmlpKfbs2QNLS0sMGTJE0962bVt06NABxcXFmnWEHrVt2zYA0JnjqJ7nERcXp7NPVlaWZmiRvtdU62Nra4vnnnsOAPTOZakqjuqsXLkSQOVwTiKihqBeW+/tt9+GEELvj/q+qd72caRSKYDKtV7/yNPTE+3bt8fFixfrfQ7aozIyMjBo0CDcv38ff//73/Hmm29Wu33fvn3h5OSE9PR0JCUl6fTru8fXNi+o89ORI0d09ikoKNAUiTXNTwA069Cqz1mTOKrz9ddfo7CwEK+++qrR1gKkRs5EbwEleuKhiqUXqtve29u7yv6qXn+9dOlSAUB06dJF3Lt3T9N++/Ztzaumq1p64dq1azWOT+3SpUuaJRb+/e9/13i/l19+WQAQL774oigrK9O0v/HGG1W+1vvLL78UAERAQIC4e/eupn379u2aZSoePZYQQigUCuHq6iosLCzErl27NO2FhYVi6NChAoAYMmSIAVcsxM8//ywACBcXF5GWlqZpP378uJDJZMLZ2VnnddkFBQXi4sWLOsf64osvhKWlpWjXrp0oKSkxKA4iotooLy8XTz31lAAgzpw5U+V2KpVKtG7dWgAQp0+fFkJUv/TCggULBADx6aef6j3exo0bNUsvnD17Vqc/JydHrFq1Sqvtccs5VOfu3bsiICBAABBvvfVWjfd7//33BQDx7LPPioKCAk17TEyMACDkcrnOPrXJC+Xl5Zq8HBsbq9U+ZcoUzWdVUVFR49gvXbokLC0thUwmE7/++qumPS0tTbi4uAgrKyutZS/UTp06pdO2c+dOYWdnJ1xdXUV2dnaNY6AnG4s9olpqqGKvvLxcDB48WAAQTk5OYuTIkWLEiBHCwcFBABAjRowQKpVKa5+6FHvqNXvc3NzEhAkT9P4sWbJEZ7/s7Gzh5+cnAAg/Pz8RGRkpOnXqpCnm7t+/r7OPSqXSrNHXvHlzMXr0aBEaGiokEomwtbUVCQkJemPcv3+/sLa2FhKJRPTu3VuMHDlSs66Qj4+P3nUMH2fGjBkCgLCzsxMRERFi8ODBwsrKSlhaWorvv/9eZ/tr165p1j8cMWKEGDNmjPD19dXE8NtvvxkcAxFRbezbt08AEIGBgY/ddtasWQKAmDFjhhCi+mLv9OnTQiKRCBsbGxERESEmT54sJk+eLHJycjTbzJkzRwAQFhYWonv37mLMmDFi9OjRIigoSFhaWgonJyetY9al2BsxYoTmPl1VftJXBBYXF4tnnnlGABAeHh5i7Nixmt/d3NxEenq63vMZmheEEOLs2bPC0dFRABBdu3YVL774oiY3uLi4aK2HV1PLly8XAISVlZUYPHiwiIiIELa2ttV+KavOxcOGDRORkZGiY8eOmhj0FYJkvljsGeBf//qX8PLyEjKZTPTt21ckJSWZOiQyoYYq9oQQoqysTHzyySciKChI2NnZCTs7O9GjRw8RGxsrysvLdbavS7Hn7e2tsxDsH3/0fQsqROVira+//rpo06aNkEqlok2bNuKNN96odhHZ8vJyERMTIzp27ChsbGyEi4uLGD16tLhw4UK1cSYlJYlRo0YJV1dXYW1tLXx9fcXMmTPr9G3lmjVrRHBwsLCzsxPOzs4iPDxcxMfH691WoVCIqVOniqefflo4ODgIOzs70bFjR7FgwQLx8OHDWsdARGSoP//5zwKAWLBgwWO3PXXqlAAgWrZsKcrKyqrNPUJUPr3r3r27prjQl1uOHDkixowZI1q1aiWsra2Fi4uL6NKli5g+fbo4cuSI1rZ1KfbUua26n6rybFFRkZg3b57w8/MTUqlUuLu7i6ioKHH9+vVqz2lIXlD77bffxPjx44W7u7uwtrYWnp6eIjo6WmRkZBh8zWq7d+8WISEhwt7eXtjb24uQkBCxZ8+eKrefOXOmCAoKEs7OzkImk4mAgAAxa9YsrVE01DRIhKhmZihpbNq0CZMmTcKqVasQHByMpUuXYv/+/UhLS9M7KZaIiIiIiMiUWOzVUM+ePdG3b1/861//AlA5Wdnd3R2LFy/G1KlTTRwdERERERGRNrN5G+eGDRvw17/+FT169IBMJoNEInnsq4BPnTqFIUOGwNnZGc2aNUPv3r2xdetWne2USiXOnj2LgQMHatqsrKwQGhqKX3/91diXQkREREREVGdms87e3LlzkZmZCVdXV3h4eFS7LhkAHD58GGFhYbCxscG4cePg4OCA7du3IzIyEtevX8dbb72l2TYnJwcqlQpPPfWU1jFatmyJ9PT0erkeIiIiIiKiujCbJ3urV69GRkYGsrOzHzussry8HFOmTIGFhQWOHj2KVatWISYmBsnJyQgMDMScOXMeWywSERERERE1ZmZT7A0cOBDe3t412vaXX35Beno6XnrpJXTr1k3T7uTkhDlz5kCpVGLdunWadldXV1haWuLu3btax7l3755RF7ImIiIiIiIyFrMp9gwRFxcHABg0aJBOX1hYGADgyJEjmjapVIqgoCAcOnRI01ZeXo64uDj06dOnfoMlIiIiIiKqBbOZs2eIK1euAAACAgJ0+tzd3WFvb6/ZRm3mzJmYPHkygoOD0b17dyxbtgxWVlZ46aWX9J6jtLQUpaWlBsVVUVGBBw8ewMXFBRKJxKB9iYioboQQePjwIVq1agULiyb5XaheFRUVuHXrFhwcHJibiIgakDHyUpMs9hQKBYDKYZv6ODo6arZRe+mll5CdnY05c+bg7t276NGjBw4cOFDlGntLlizBwoULjRs4ERHVu+vXr8PT09PUYTQat27dQps2bUwdBhFRk1WXvNQki73amjFjBmbMmFGjbd977z3MmjXLoOMrFAp4eXnh+vXrXKidiKiB5efno02bNnBwcDB1KI2K+vNgbiIialjGyEtNsthTP9H749M7tfz8fDRv3rxO55DJZJDJZLXa19HRkQmViMhEOFRRm/rzYG4iIjKNuuSlJjkpQT1X74/z8gDgzp07KCgo0Dufj4iIiIiI6EnRJIs9uVwOAPjpp590+g4cOKC1DRERERER0ZOoSRZ7zz//PHx9fbFp0yYkJSVp2hUKBf7xj39AKpVi/PjxJoyQiIiIiIiobsxmzt7q1atx7NgxAMC5c+c0beo19fr164fo6GgAgJWVFVavXo2wsDD0798f48aNg4ODA7Zv347MzEwsW7YMPj4+prgMIiIiIiIiozCbYu/YsWNYt26dVlt8fDzi4+M1v6uLPQAYMGAAjh07hgULFmDLli0oKytD586d8fHHHyMyMrJeY42NjUVsbKxOu0qlqtfzEhERVYW5iYjI/EiEEMLUQVCl/Px8ODk5QaFQ8I1nREQNjPdg/fi5EBGZhjHuv01yzh4REREREZG5Y7FHRERERERkhljskZbDhw9DIpFg0qRJpg5FS3FxMebPn4/AwEDY2NigVatWmDRpEm7evFnjY8TFxUEikTz2Z9GiRfV4JUREZChzzk1/lJ2djbfffhvt2rWDra0tWrRoge7du+Odd94xYuRE1FSYzQtayDjOnj0LAOjRo4eJI/ldSUkJnnvuOSQkJMDDwwMRERHIyMjAmjVr8MMPPyAhIQG+vr6PPY67uzsmTJigt0+lUmHDhg0AgJCQEKPGT0REdWPOuelRZ86cQVhYGO7fv4+OHTsiIiIC+fn5uHjxIlasWIGlS5fW09UQkblisUdaEhMTAQDBwcEmjuR3ixcvRkJCAvr06YOffvoJ9vb2AIDly5fjrbfewqRJkzRLbFSnffv2WLt2rd6+H3/8ERs2bECbNm0QGhpqvOCJiKjOzDk3qWVnZyM8PBzFxcXYtWsXXnjhBa3+kydPGjN8ImoiOIyTtJw9exZWVlbo2rWrqUMBACiVSnz66acAKl8Lrk6mADBr1ix06dIFR44cwZkzZ+p0HvVTvZdffhkSiaROxyIiIuNqCrlpwYIFyMnJwdKlS3UKPQDo1atX3QMnoiaHxZ4JxMbG4umnn9b56dmzp0njKioqwuXLl9GxY0colUrMnj0bvr6+kMlkaN++PVatWtXgMcXHx0OhUMDPzw9BQUE6/aNHjwYA7Nmzp9bnKCwsxK5duwAAr776aq2PQ0T0JGNuqjlj56bi4mJs2LABzZo1w8SJE40aKxE1bRzGaQLTpk3DtGnTdNrVa2mYSkpKClQqFTw8PBAUFITS0lL06dMHHh4eOH78OP76179CqVRi+vTpDRZTcnIyAKB79+56+9XtKSkptT7Hjh07UFhYiKCgIDz99NO1Pg4R0ZOMuanmjJ2bTp8+jYcPH6Jfv36wtbXFjz/+iJ9//hklJSUIDAzE2LFj0apVK+MET0RNCp/skYZ6TsT+/fsxbNgwZGRk4LvvvkN8fDz+85//AAAWLlwIlUpV5TFCQ0Nr9MbLR3+qmkcHAFlZWQAAT09Pvf3q9szMzNpcMoDfh3DyqR4RUePTFHLTxYsXAQAtW7bEiBEjMGTIEKxYsQKff/45Zs6cCX9/f2zevLlGxyIiehSf7JGG+m1nQ4YMwSeffKLV95e//AVLly7F1atXkZqaik6dOuk9Rnh4OHx8fAw6r7+/f5V9BQUFAAA7Ozu9/c2aNQMAPHz40KBzqt2+fRuHDh2CpaUl/vznP9fqGEREVH+aQm7Kzc0FAOzevRuWlpaIjY3FmDFjUFRUhE8//RTLli3DhAkT0KFDB3Tr1s2QyyCiJo7FHmmovz1dsmSJ3v4OHTrg6tWryM/Pr/IYs2fPrpfY6svmzZuhUqkQHh4Od3d3U4dDRER/0BRyU0VFBQCgvLwcH374IV577TVN39KlS5GZmYnvvvsOS5cuxcaNG00VJhE9gTiMkwAAZWVlOH/+PDp06IAuXbro3aakpAQA4Obm1mBxqd9wVlRUpLe/sLAQAODg4FCr43MIJxFR49VUctOjb/PU94IWdduRI0cMipOIiE/2CABw4cIFKJXKKhKMIVsAACAASURBVCebCyGQmJgIJycntG3btsrjfPTRR7h06ZJB546Ojka/fv309nl5eQEAbty4obdf3e7t7W3QOQEgNTUVZ8+ehb29PUaMGGHw/kREVL+aSm5Sb2dnZ6e3aFUPQb13716NjkdEpMZijwD8Pieiqm8hf/75Z9y/fx/jxo2DlVXV/9ns37/f4G8eQ0NDq0yo6jWV1MN4/kjdXtU3vtVZv349AGDUqFFVzrsgIiLTaSq5Sb18Q3FxMUpLSyGTybT6Hzx4AED7CSARUU2w2CMAvyem9PR0nT6VSoX58+cDAN55551qjxMXF2fUuPr27QsnJyekp6cjKSlJZ2L6tm3bAADDhw836LhCCGzatAkAh3ASETVWTSU3eXl5oWvXrkhOTsaRI0cwaNAgrX51oapvTT8ioupwzp4JNMaFa9Xfnh46dEjr28+SkhJMnjwZJ06cwLRp06ocSlNfpFKpZu2kadOmaeZBAMDy5cuRkpICuVyO4OBgnX3Hjx+P9u3b4/vvv9fp++9//4vMzEy0bt0azz33XP1dABHRE4K5qebqIze9++67AIC3334bt2/f1rQnJSUhJiYGADB16lSjXwsRmTlBjYZCoRAAhEKhaNDzqlQqYW9vL3x8fERkZKSwtrYW4eHhYuzYsaJly5YCgBgyZIgoLS1t0LjUiouLxTPPPCMACA8PDzF27FjN725ubiI9PV3vfnK5XAAQa9as0embMmWKACDeeeedeo6eiJ4UproHN3bMTfrVR26aMGGCACCcnZ3FkCFDxIABA4RMJhMAxJQpU+r5ioiosTHG/ZdP9ghXrlxBQUEBnn76aaxbtw7Tp0/H2bNnsXv3brRu3RqfffYZdu/eDalUapL4bGxscPjwYcybNw92dnbYuXMnMjMzERUVhcTERPj6+hp0vNLSUs0Qm1deeaU+QiYiojpqarkJANasWYNVq1bBz88PcXFxOHnyJLp37461a9di1apV9XAVRGTuJEIIYeogqFJ+fj6cnJygUCjg6Oho6nCIiJoU3oP14+dCRGQaxrj/8skeERERERGRGWKxR0REREREZIZY7BEREREREZkhFntERERERERmiMUeERERERGRGWKxR0REREREZIasTB1AUxQbG4vY2FiddpVKZYJoiIiImJuIiMwR19lrRLiWERGR6fAerB8/FyIi0+A6e0RERERERKQXiz16IixfvhyjRo1CQEAAnJycIJPJ4O3tjfHjx+PcuXO1OmZxcTHmz5+PwMBA2NjYoFWrVpg0aRJu3rxp5OiJiMgcnTlzBh999BFGjRoFT09PSCQSSCQSox3//v37aNmyJSQSCfz9/Y12XCJqOjiMsxHhUJmqubq6orCwEF26dEHr1q0BABcuXEBaWhqsra2xY8cODBs2rMbHKykpwYABA5CQkAAPDw+EhIQgIyMDJ0+ehJubGxISEuDr61tfl0NEjRDvwfrxc6naiBEjsGvXLp12Y/1pFRUVhW+++QZCCPj5+eHq1atGOS4RPRk4jJOajF27diE3NxcnTpzAjh07sGPHDly+fBmxsbEoKytDdHQ0ysvLa3y8xYsXIyEhAX369EFaWhq2bNmCEydOICYmBtnZ2Zg0aVI9Xg0REZmDPn36YN68edi9ezdu374NmUxmtGMfOnQI69atw5QpU4x2TCJqevhkrxHht6e14+/vj/T0dCQnJ6NLly6P3V6pVKJly5ZQKBRITExEUFCQVn/Xrl2RkpKC06dPIzg4uL7CJqJGhvdg/fi51JyNjQ1KS0vr/GSvuLgYnTt3hkwmw86dOxEYGMgne0RNEJ/skVHt2bMHEokEr7/+ut7+6OhoSCQSHDp0qIEjq561tTUAQCqV1mj7+Ph4KBQK+Pn56RR6ADB69GgAlZ8HERGZ1pOam+pi4cKF+O233/DFF19ochwRUW2w2CONs2fPAgC6detWq35TWL9+PS5fvoyAgAAEBATUaJ/k5GQAQPfu3fX2q9tTUlKMEyQREdXak5ib6iIlJQUxMTGYOHEiQkJCTB0OET3hWOyRhjph6nvaVVZWhvPnz8PT0xMuLi5VHiM0NFTzNrKa/qxdu7bGMS5duhRRUVEYM2YMOnXqhPHjx8PDwwObN2+GpaVljY6RlZUFAPD09NTbr27PzMyscVxERFQ/noTcZCwVFRWIjo6Gs7Mz/vnPfzb4+YnI/FiZOgBqPM6ePQtra2t07NhRp+/ixYtQKpWP/eY0PDwcPj4+Bp3XkNdJHzhwQGuojre3N7755huD5tYVFBQAAOzs7PT2N2vWDADw8OHDGh+TiIjqx5OQm4xl5cqVOHXqFNasWVNt8UpEVFMs9ggA8ODBA2RmZmomhP9RUlISgMcPk5k9e3a9xKd28OBBAEBeXh7OnTuHRYsWQS6XY/HixXj//ffr9dxERNSwnpTcZAxZWVmYO3cu5HI5oqKiTB0OEZkJFnsmEBsbi9jYWJ12lUplgmgqVTdM5tH+xjInwtnZGSEhIdi3b5/m1deDBg1Cz549H7uvvb09AKCoqEhvf2FhIQDAwcHBeAETETVyzE2mNW3aNCiVSnzxxRemDoWIzAiLPROYNm0apk2bptOufr2qKRhrAvxHH32ES5cuGXTu6Oho9OvXz6B91KytrREZGYkzZ85gz549NSr2vLy8AAA3btzQ269u9/b2rlVMRERPIuYmbXXJTbXxww8/wNnZGVOnTtVqLykpAQDcvHkToaGhAIBvv/0W7u7uDRYbET25WOwRgN8TZufOnXX6Hj58iISEBDg6OsLX17fa4+zfvx9Hjhwx6NyhoaF1Sqiurq4AgOzs7Bpt37VrVwBAYmKi3n51e03W7CMiovrzJOem2sjLy6syzpKSEk2fugAkInocFnsE4PcCR99LS9atWwelUolevXpBIpFUe5y4uLj6CK9a6uTn5+dXo+379u0LJycnpKenIykpSecb4W3btgEAhg8fbtxAiYjIIE9ybjJUVQuxZ2RkoG3btlxUnYhqhUsvEAoLC5GWlgYA2LBhg1bC2bdvH/7v//4PABAYGGiS+OLj47F//35UVFRotZeVlWHlypVYv349bG1tERkZqdU/fvx4tG/fHt9//71Wu1QqxfTp0wFUDltSz9EDgOXLlyMlJQVyudygN3wSEZFxNfbcVFtV5SYiovrAJ3uElJQUVFRUwNfXF59//jkOHToEf39/XLt2DampqejZsydOnTqFnTt3oqKiAmvWrGnQ+K5cuYKJEyfC1dUVwcHBcHFxQU5ODs6dO4fbt2/DxsYGa9euRZs2bbT2y8rKwuXLl6FQKHSOOXfuXBw8eBDHjx9HQEAAQkJCkJmZiRMnTsDNzQ1ff/11Q10eERHp0dhzEwDs3bsXH3zwgeZ3pVIJAOjdu7embd68eRg6dKjm9+pyExGRsbHYI82ciAkTJuCpp57CkiVLcOjQIbRr1w6bNm1CcHAwBg4ciJKSEnTv3r3B45PL5ZgzZw6OHDmClJQU5OTkQCqVwsfHB6NHj8Ybb7xh8HpINjY2OHz4MJYsWYJNmzZh586daNGiBaKiovDBBx9UueA6ERE1jMaem4DKueInTpzQaX+0rabzyYmI6oNEVDVInBqc+o1nCoUCjo6ODXbe6OhofPXVV9izZw+GDRvWYOclImpMTHUPbuyYm4iITMMY91/O2SPNt6eco0ZERI0FcxMRUd2x2GviysrKcP78ebi7u8PDw8PU4RARETE3EREZCYu9Ju7ixYtQKpUmm+9ARET0R8xNRETGwRe0NHFdu3atcm0fIiIiU2BuIiIyDj7ZIyIiIiIiMkMs9oiIiIiIiMwQiz0iIiIiIiIzxDl7JhAbG4vY2FiddpVKZYJoiIiImJuIiMwRF1VvRLigb9WKi4uxZMkSfPvtt8jKykKLFi0QHh6ODz74AK1btzboWMuXL8exY8dw7tw53Lt3DyUlJXB3d4dcLsc777yDzp0719NVEFFjxnuwfvxcqmbM3PRH9+/fR4cOHZCdnQ0/Pz9cvXrVSFET0ZPCGPdfFnuNCBOqfiUlJRgwYAASEhLg4eGBkJAQZGRk4OTJk3Bzc0NCQgJ8fX1rfDxXV1cUFhaiS5cummR84cIFpKWlwdraGjt27MCwYcPq63KIqJHiPVg/fi76GTs3/VFUVBS++eYbCCFY7BE1Uca4/3LOHjV6ixcvRkJCAvr06YO0tDRs2bIFJ06cQExMDLKzszFp0iSDjrdr1y7k5ubixIkT2LFjB3bs2IHLly8jNjYWZWVliI6ORnl5eT1dDRERmQNj56ZHHTp0COvWrcOUKVOMGDERNUV8steI8NtTXUqlEi1btoRCoUBiYiKCgoK0+rt27YqUlBScPn0awcHBdT6fv78/0tPTkZycjC5dutT5eET05OA9WD9+LrrqMzcVFxejc+fOkMlk2LlzJwIDA/lkj6iJ4pM9qhcbN27EwIED4erqCktLS0gkEp2f0aNHN0gs8fHxUCgU8PPz00mmADRx7Nmzxyjns7a2BgBIpVKjHI+IiIyjqeSmhQsX4rfffsMXX3yhyUlERLXFt3GSRkVFBf785z9j69atkMlkkMvlaN68OeLj43Hjxg2tbbt169YgMSUnJwMAunfvrrdf3Z6SklLnc61fvx6XL19GQEAAAgIC6nw8IiKqu6aUm1JSUhATE4OJEydq5gASEdUFiz3S+OCDD7B161Z06NABe/fuRdu2bQFUDikZMGAATpw4gYULF+Ldd9+t8tvG0NBQHDlyxKDzrlmzBlFRUXr7srKyAACenp56+9XtmZmZBp0TAJYuXYoLFy6gsLAQqampuHDhAlq1aoXNmzfD0tLS4OMREZHxNZXcVFFRgejoaDg7O+Of//ynQbESEVWFxR4BAO7evYslS5ZAKpVi69atmmQKALa2tpg0aRJOnDiB+Ph42NjYVHmc8PBw+Pj4GHRuf3//KvsKCgoAAHZ2dnr7mzVrBgB4+PChQecEgAMHDuDQoUOa3729vfHNN98YZe4fERHVXVPKTStXrsSpU6ewZs0auLi4GBApEVHVWOwRgMohjKWlpYiOjkanTp10+jt06AAAyM3NrfY4s2fPrpf46sPBgwcBAHl5eTh37hwWLVoEuVyOxYsX4/333zdxdERE1FRyU1ZWFubOnQu5XF7l00QiotrgC1oIwO+Fz6hRo/T25+XlAQDc3NwaLCYAsLe3BwAUFRXp7S8sLAQAODg41Poczs7OCAkJwb59+xAcHIx58+bh1KlTtT4eEREZR1PJTdOmTYNSqcQXX3xhnACJiP6HT/YIAHDu3DkAQK9evfT2q4ufHj16VHucjz76CJcuXTLo3NHR0ejXr5/ePi8vLwDQmYSvpm739vY26Jz6WFtbIzIyEmfOnMGePXvQs2fPOh+TiIhqr6nkph9++AHOzs6YOnWqVntJSQkA4ObNmwgNDQUAfPvtt3B3d6/RcYmIWOwRACAnJwcAqlzD4/vvvwcADBs2rNrj7N+/3+BJ8KGhoVUm1K5duwIAEhMT9far2421Jp6rqysAIDs72yjHIyKi2mtKuSkvL6/KGEtKSjR96gKQiKgmOIyTAPyeSPUt2vr111/j/PnzePbZZx/7tCsuLg5CCIN+qpuf0LdvXzg5OSE9PR1JSUk6/du2bQMADB8+3ICrrZo6mfr5+RnleEREVHtNJTdVFcO1a9cAVOYkdZuhL5ohoqaNxR4BAPr37w+g8hXXFRUVmvaffvoJb775Juzs7PDZZ581eFxSqRTTp08HUDmnQT0PAgCWL1+OlJQUyOVyvW/QHD9+PNq3b6/55heoXAh3//79WtcIAGVlZVi5ciXWr18PW1tbREZG1tMVERFRTTWV3EREVF84jJMAAH//+9+xf/9+bN68GWfPnkVQUBCuXbuGhIQEODk54bvvvtMMW2loc+fOxcGDB3H8+HEEBAQgJCQEmZmZOHHiBNzc3PD111/r3S8rKwuXL1+GQqHQtF25cgUTJ06Eq6srgoOD4eLigpycHJw7dw63b9+GjY0N1q5dizZt2jTU5RERURWaSm4iIqovfLJHAIDOnTvj2LFjGDZsGO7evYvt27cjOzsbr7/+Os6dO4ewsDCTxWZjY4PDhw9j3rx5sLOzw86dO5GZmYmoqCgkJibC19e3xseSy+WYM2cO2rVrh5SUFHz33XeIj49HixYtNNc6duzYerwaIiKqqaaSm4iI6otECCFMHQRVys/Ph5OTExQKRZWT0YmIqH7wHqwfPxciItMwxv2XT/aIiIiIiIjMEOfsmUBsbCxiY2N12lUqlQmiISIiYm4iIjJHHMbZiHCoDBGR6fAerB8/FyIi0+AwTiIiIiIiItKLxR4REREREZEZYrFHRERERERkhljsERERERERmSEWe0RERERERGaIxR4REREREZEZYrFHRERERERkhljsERERERERmSEWe0RERERERGaIxR4REREREZEZYrFHRERERERkhljsERERERERmSEWe0RERERERGaIxR4RkYkUKcvhM3svfGbvRZGy3NThEBERkZlhsUdERERERGSGrEwdABFRU6N+ilekVD3S9vu/7aS8NRMREVHd8S8KIqIG9vT8AzptPRYf1Pw746OhDRkOERERmSkO4yQiIiIiIjJDfLJHRNTALi4KA1A5dFP9RO/03IGwk1qaMiwiIiIyMyz2iIgamL45eXZSS87VIyIiIqPiME4iIiIiIiIzxK+RiYhMxE5qxZexEBERUb3hkz0iIiIiIiIzxGKPiIiIiIjIDLHYIyIiIiIiMkOcs2cCsbGxiI2N1WlXqVQmiIaIiIi5iYjIHEmEEMLUQVCl/Px8ODk5QaFQwNHR0dThEBE1KbwH68fPhYjINIxx/+UwTiIiIiIiIjPEYo+IiIiIiMgMsdgjIiIiIiIyQyz2iIiIiIiIzBCLPSIiIiIiIjPEYo+IiIiIiMgMsdgjIiIiIiIyQyz2iIiIiIiIzBCLPSIiIiIiIjPEYo+IiIiIiMgMsdgjIiIiIiIyQyz2iIiIiIiIzBCLPSIiIiIiIjPEYo+IiIiIiMgMsdgjIiIiIiIyQyz2iIiIiIiIzBCLPSIiIiIiIjPEYo+IiIiIiMgMsdgjIiIiIiIyQyz2iIiIiIiIzBCLPSKiGihSlsNn9l74zN6LImV5jfuIiIiITIXFHhERERERkRmyMnUARESNmfpJXZFS9UibStNnJ7XS2wcAdlLeYomIiMh0+JcIEVE1np5/QKetx+KDVW7/aF/GR0PrJSYiIiKimuAwTiIiIiIiIjPEJ3tERNW4uCgMQOXwTPVTu9NzB8JOaqk1jPOPfURERESmxmKPiKga+ubd2UktYSe1qraPiIiIyNQ4jJOIiIiIiMgM8etnIqIasJNaVfnCler6iIiIiEyFT/aIiIiIiIjMEIs9IiIiIiIiM8Rij4iIiIiIyAyx2CMiIiIiIjJDLPaIiIiIiIjMEIs9IiIiIiIiM8Rij4iIiIiIyAyx2KuhHTt24E9/+hNatGgBiUSCjIwMU4dERERERERUJRZ7NVRYWIj+/ftj0aJFpg6FiIiIiIjosaxMHcCT4tVXXwUAnD9/3sSREBERERERPV6jfrK3YcMG/PWvf0WPHj0gk8kgkUiwdu3aavc5deoUhgwZAmdnZzRr1gy9e/fG1q1bGyZgIiIiIiKiRqJRP9mbO3cuMjMz4erqCg8PD2RmZla7/eHDhxEWFgYbGxuMGzcODg4O2L59OyIjI3H9+nW89dZbDRQ5ERERERGRaTXqJ3urV69GRkYGsrOzMXXq1Gq3LS8vx5QpU2BhYYGjR49i1apViImJQXJyMgIDAzFnzhydYnH27NmQSCTV/hARERERET2JGvWTvYEDB9Z4219++QXp6emYOHEiunXrpml3cnLCnDlzEBUVhXXr1mH+/PmavrfeegtRUVHGDJmIiIiIiKhRaNTFniHi4uIAAIMGDdLpCwsLAwAcOXJEq93NzQ1ubm71HhsREREREVFDM5ti78qVKwCAgIAAnT53d3fY29trtqmNBw8eICsrC+np6QCAixcvIi8vD15eXmjRooXO9qWlpSgtLTXoHPn5+bWOj4iIiIiI6FFmU+wpFAoAlcM29XF0dNRsUxu7d+/GxIkTNb8PHToUALBmzRq9Q0GXLFmChQsX1vp8REREREREdWE2xV59i4qKMmh+33vvvYdZs2YZdI78/Hy0adPGwMiIiIiIiIh0mU2xp36iV9XTu/z8fDRv3rzB4pHJZJDJZA12PiL6XZGyHE/PPwAAuLgoDHZSqxr1mQNzvz4iIiKquUa99IIh1HP19M3Lu3PnDgoKCvTO5yMiIiIiIjJHZlPsyeVyAMBPP/2k03fgwAGtbYjIPBUpy//3o3qkTYUiZTlyCkqq7CtSlpsiXKOq7trN4fqIiIjIcGYzvuf555+Hr68vNm3ahDfeeEOz1p5CocA//vEPSKVSjB8/3sRRElF9Ug9ffFSPxQer3P7RvoyPhtZLTA3lcdf+pF8fERERGa5RF3urV6/GsWPHAADnzp3TtKnX1OvXrx+io6MBAFZWVli9ejXCwsLQv39/jBs3Dg4ODti+fTsyMzOxbNky+Pj4mOIyiIiIiIiIGlyjLvaOHTuGdevWabXFx8cjPj5e87u62AOAAQMG4NixY1iwYAG2bNmCsrIydO7cGR9//DEiIyMbLO7HiY2NRWxsrE67SqXSszUR1dTFRWEAKocvqp9qnZ47EHZSSxQpy2EntdLbZw6qu3aimmBuIiIyPxIhhDB1EFQpPz8fTk5OUCgUcHR0NHU4RE8svo3TfK+vPvEerB8/FyIi0zDG/ddsXtBCRERE9CQrUpbDZ/Ze+Mzea/IXKzWmWIio9viVLxGZHTupVZUvJNHXV66qQPINBf57JRtnMnPR2tkW4Z3c8ayfK6RWT9Z3YtVdOxERETUtLPaIqEm6/qAIR69k479pOYhPz8HDEu1vrr89dR2ONlYY2OEphHdyR/9AN9hYc/4bERmf+snZH5dOUWvI4diNKRYiqjvO2WtEOC+CqP6l3s7HjG/PIu1ugVa7o40V+gW44pm2Lki7+xAHLtxFTkGppt9OaokB7Vti4rM+6OHToqHDpgbAe7B+/Fzqn8/svdX2N+TT+sYUC1FTZ4z7L7+eIaImI/V2Pl76MgG5RWWwspCgu1dzhAS4IiTQDZ1bO8HSQqLZdlFEJ5zJzMWP52/jwPk7uKUowd6U2/jx3G28E9YeU+W+kEgk1ZyNiIiIyLT4ZK8R4benRPXn8p2H+POXCXhQqERXTyesndgLzZtJa7SvEALJNxRYE38Nu5JuAQD+9PRTiBnbFY421vUZNjUg3oP14+dS/x4dOqlv6RRTDeM0dSxETR2f7D2huJYRUcNKu/sQL/2v0Ovi6YRvJj8DJ9uaF2kSiQTd2jjjX5Hd8ExbF/x99wX8fPEuXlh5DJ+/EowOHvwDmJ58zE2mo6+AspNamqSwakyxEFHd8cleI8JvT4mM78rdyid6OQVKdGrtiI2Te8PJrm5P41Ju5OFvGxJxM68YNtYW+HBEZ7wY7GmkiMlUeA/Wj59Lw2lM62Q2pliImipj3H9Z7DUiTKhExnX13kOMW3UCOQWleNrDEZumPANnu5oN3Xyc3EIlZmxJwtG0bADAy894Yf7wpyGz4hs7n1S8B+vHz6VxYPFF1PRwUXUioipcvVegKfQ6eDhiY7TxCj0AaN5MijVRPfHmwABIJMDGE1l4dfVJFCs55I2IiIgaBxZ7RGR2buYV46UvE5BTUIr27g7YGP1MjV/GYghLCwneHBiINVE94WBjhZMZD/DaxjMoU1UY/VxE1DQVKcv/96O97p26nYioOhwDQERm5+MfL+Hew1K0e8oBm6b0Rot6KPQeFdquJdZE9cQrX53A4cvZePu7ZKwY2w0WFlyagagpMuaQS/VxHqV+SybAde+IqHp8skdEZuX8TQV2J1cuj7A8smu9F3pqPXxa4PNXgmFlIcGupFtYuOcCOCWaiIiITIlP9ojIrHy8/xIAIKJbK3Rs5dSg5x7QriVixnbFm1uSsO7XTDRvJsWbAwMbNAYiMp1H16j7ve33f9fmCd/FRWGa4+hb946IqDos9ojIbMRfzcF/r+TA2lKCt/7UziQxRHRrjbyiMizYfQH/OngFze2kmPCsj0liIaKGVR9DLrnuHRHVBe8UJsCFa4mMTwihear38jPe8HKxM1ksE571QW6REv86eAULdl+As501Irq1Nlk8RDXB3EREZH64zl4jwrWMiGpvb8ptTNuUiGZSSxx5dwBc7WUmjUcIgYV7LmLt8QxYWUjw5YQeGNCupUljourxHqwfP5eae3QYp74hl3waR0SGMPk6ewsXLsSNGzfqcggiojorU1Vg2U+XAQDRIb4mL/QAQCKRYP6wpxHRrRXKKwSmbUzElbsPTR0WEdUjO6nV/34sH2mz1LQTETW0Ohd7bdu2xfDhw7F7925UVHBtKSJqeFtOXce1nEK4NJNiSn9fU4ejYWEhwbIxXdHH1wVFShX+uv4MHpaUmTosImoihBBQFJUh7e5DKIp07z1FynL4zN4Ln9l7uWYfkZmq09dMixcvxldffYW9e/di3759cHd3x6RJkzB58mT4+PgYKUQioqoVKcvxyaErAIDXn/OHvaxxfXtubWmBlS8FYfjKY/gtpxBvf5eML14JhkTCNfiIzJWd1KpB178rKVPhwIU7+C27ELcVxbitKMGtvMr/ffRtoL5uzdDdqzmCvJzR3as5PJvbNliMRGQaRpmzd/DgQaxevRo7d+6EUqmEhYUFBg4ciClTpiAiIgJWVo3rj6/GivMiiAwXe/gqlh64jDYtbHFoViikVo1z+dCk63kY+8WvUKoq8G54O7wW6m/qkOgPeA/Wj59L41WsVGHjiUz85+hvyH5YWuV2DjZWeFii++TOTmqpKQaPz34OznbW/2vn321EFjWFzAAAIABJREFUjYEx7r9GfUHL/fv3sW7dOnz11VdITU2FRCKBm5sboqKiMHnyZAQEBBjrVGaJCZXIMLmFSvT/52E8LC3HJ+O6Nfo3Xm4+mYX3dpyDhQRYN6kXQgLcTB0SPYL3YP34uTQ+BaXlWP9rJlb/9zfcL1QCAFo52UDezg0eTrbwcLJBK+fK//VwsoWt1BIPCpVIup6LxMw8fHr4arXHb8inkkRUtUZX7D3q+PHj+Oyzz7Bp0ybNcKXQ0FBMnz4dI0eOrI9TPvGYUIkMs/iHi1h97Bqe9nDED6/3g4VF4x8a+X/bUrDl9HU0t7PGntf7wbO56ZaIIG28B+vHz6XxyC8pw7r4DHwVfw15/5uD59XCDq+F+mFUd88aj2zwmb232v6PRnXGi8GesLZsnCMliJoKk7+Nsyrp6enYs2cPDh06pGnz9PTE4cOHMXr0aPTq1QvXr1+vj1MTURNxM68Y3/yaCQB4N7zdE1HoAcDCiI7o4umE3KIy/G1DIkrKuIYZET3e4Uv30O+jXxDzcxryisrg69oMMWO64pe35BjXy8ugIewXF4Xh4qIwnJ47UNP2dlggXO2lAIDZO87huZg4bD19HeUq4718jy+EIWp4Riv2ysrK8O233+L5559HYGAgPv74Y5SXl2PWrFm4dOkSMjMzER8fj8GDB+P06dOYPn26sU5NRE3QhoRMKFUV6O3bAvLAJ2c4pI21JT5/JRgtmklx7qYC83f9P3t3HhZ1vf0B/D0zMMCwr7LLJiqigJqhpuaSpqaplaZ2b93qljcszeqmRXYzU9uvFnlb7OYvc8vUNNfrmnuioCguKJuyyr4Jwyy/P2BGDRSGGfjO8n49zzwPfGeYOcL4nTnzOZ9zzoHjTonoXjacuo7n/y8RFbUKdPFywNIno/G/OUPwWB9/WLVh9a25ERHPDgzG4TeH4Z1HIuDhIMW1kpv454azGPflEVy9UWXIfw4RdSC9d+BeuHAB3377LX788UeUlJRArVZjwIABmDFjBp544gnY2Nyad9W/f3/89ttviI2NxcGDB/V9aJOVkJCAhISEJseVSn7CT9QaKpUaW5JzAQB/iQ0yuc6Wfi52+GJqDP6y4gTWJ15HdIArpt0fKHRYZOH42mR81Go1vv49HUt2XAQATIrxw4eP92q38kpbawmeeyAYU/sFYNXxLCw/cBUX8iow7ovDWDSxJybEtG1f9O3D5m8du/U1G8IQtR+99uw98MADOHbsGNRqNZycnPDUU09hxowZiIyMvOfPPffcc/jhhx/4AvIn3BdB1Dp/ZJRg8tfH4GBjhcT4EbC1lrT8Q0Zo+YGr+HDnRUglYmyOG4gIX/6/FxLPwc3j70UYKpUaC7ddwPdHMgAALw4OwZsPd9O5ZL1GrkDE/F0AGso3dUmsCitqMWttMo6lFwMApvQNwL/G94CdVLdzbkt7BNkQhqh5gu/ZO3r0KGJiYvDNN98gNzcXX375ZYuJHgA8//zz+P777/V5aCKyYJuTcwAAD0d6m2yiBwAzhoRgRPdOkCtVmLnmNPewEBEAQK5QYfa6ZG2iFz+2O+aN6d7he5O9nGyx6vn7MXtEF4hEwLrEaxj/5WFcLqjs0DiIqO30Wjc/efIk+vTpo/PP9e/fH/3799fnoYnIQskVKmxPyQMATDDyUQstEYlE+PjxXhi99BDSb1TjvS2p+PDxXkKHRUQCqqpTYMaPp3D4ShGsxCJ88kRUm8onDVU6KRGLMHtEOPoFu2HW2mSkFVZh/JeHseDRSDzRx79VZfSpC0ZpH7/vwj0AgMT4EXfsGSSi9qHXyt62bduwZcuWFm+3detWLFiwQJ+HIiICABy8fANlNfXwcrRB/1B3ocPRm6u9FJ9PidZ+ar71TK7QIRGRQEqr5Zj6zXEcvlIEmVSCFc/c16pEr7kulxHzdyFi/i5tcgUAfRfu0R7X1YBQD+yYNQiDunigtl6Ff244i39uOIv6VnTrbK4hjEwq0R4novajV7L3r3/9C5s3b27xdlu2bMF7772nz0MREQG4VcI5LsoXEhMZt9CS/qHumDk0DADw1sYUXCupETgiIupoKpUar65PRkpOOdzspVjz91ij6zTs4WCDlX/rhzdGdYVELMLPp67jH6tOcYQMkRHrkI9TlEolxGIO5iQi/VTVKbAntQAA8Gi0r8DRGNas4V1w9GoxTmWV4pW1SVj/Yn8ONCayIN8dTseBSzdgYyXGqufub1XDpnuVaibGD4dMamXw0kmxWIS4oWHo7uOIf6w6jT0XCvH093/gu6f7wtHW+p4/K5NasRkLUQfrkGTv/PnzcHV17YiHIiIztutcPuoUKoR42KOnn7PQ4RiUlUSMpU9GY/TSQ0jKLsPn/7uMfz7cTeiwiKgDnM4uxUc7LwEA5o+LaHVn3ubKMW8v2/xzYqUpndRXjVyBZ39IBADY20hwIqME0787gR/+1g9u9lK975+IDEfn//HPPvvsHd8fPny4yTENhUKBS5cuITExERMmTGhbhEREjTQlnI9G+5ncbL3W8HeV4cPHeuGln05j+cGrGBjmgYFhHkKHRUTtqPxmPV5ZkwSFSo2xvXwwrZ9pzdz84W/34cUfT+Ps9XJM/voYVj13P7ydbYUOi4ga6Txn7/ZyTJFIhNb8eK9evbBx40aEhIToHqEF4SwjorsrrKxF7KK9UKmBA68/iCAPe6FDajfzNqZgzR/Z8HK0wY5Zg+DuYCN0SBaB5+Dm8ffSftRqNV766TR2nMtHoJsMv73yAJxaKIW83e1lnM2Vahq6+cndHi+//CaeX3kK+RW18He1w6rn7jfrczRRRzHE+Vfns8D+/fsBNJyghg0bhocffhhvvvlms7eVSqXw9fVF586d2xQcEZHGb2fyoFID0QEuZv8mYv4jETiZWYIrhVV4Y8NZrHi6r1muZBJZulXHs7DjXD6sJSJ8MTVGp0QPaD6ZM1SpZnNaKhsN9rBHRlE1Hv/PMfz4XD909+GHA0RC0/lsMGTIEO3XTz/9NAYNGnTHMSKi9vCrtoTTvBqzNMdOKsGX02Iw/ssj2HexEKuOZ+Ev/YOEDouIDOh8bjne33YBAPDmw90QFeAicET6W/9ifzy14gQu5Vdi9NJD+HXmQET5m/6/i8iU6VzGSfpLSEhAQkJCk+NKpRKXL19mqQzRn2QUVWPoJwcgEYtwfN5weDpaRlnjf49k4L2tqbCxEmPbKw8gzMtR6JDMmqWXK/K1qeNU1ykw7ovDSC+qxvBuXvjORFbvW1M2ml9+E7GL9wEAvBxt8Ms/BiDATSZMwEQmTpAyTtJfXFwc4uLimhzX/EGJ6E6aVb2BYR4Wk+gBwNP9g7DvYiEOpRVh9rpkbPzHQEitOI6B2gdfmzrOO5vPIb2oGj7OtvjkiSiTSPSAe5eN1sgVqJErYHXbyJjCyjpM/+4EfnyuHzq7m3f5PZGx0inZCwkJgUgkwp49exAcHKxTwxWRSISrV6/qHCARWTa1Wo1fk3MBABPauYSzRq7Q7klJXTDKYC3K23qfYrEInzwRhVH//h3ncirw7z3CjGMwtt8LkSn75dQ1bExq+ADrw8d6wdVMRhU0t58PALJLajDk4wM48+5IONvptieRiPSn06trZmYmAKC+vv6O74mI2svZ6+XIKKqGrbUYI3t4Cx1Oh+vkZIslk3pixqqGcQwPdvVCv2A3ocMiojaorVfik92Xtd/3DTLNGcRtGY7+7A8n8eNz/fjBDlEH06keSKVSQaVSITw8/I7vW3shItKVZrbeQxHecLBpnzcJmvKjGrnytmNK7XGh7/PhSB880ccfajXw6rpkVNTWtykmXRn774XIlNTIFfjuUDryymtvO2Y+z/3UBaOQumAUEuNHaI8lxo/A5rgBcLK1wqmsUvxj1WnIFXw/SNSR2KDFiFh6cwCiP1MoVYhdvA9FVXX47q99MSKiU7s8TtDcbfe8XtdPsNvjPqvqFBiz9BCyS2owMcYPn0+J1jkmXZnC78WQeA5uHn8vhmHMz31Daq5E+1RWKZ767gRu1isxtpcPlj0ZA4nYNPYpEgnJEOdf7vQnIqN19Goxiqrq4CKzxuBwT6HDEZSDjRU+nxIFsQjYlJSDLWdyhQ6JiKhV+nR2xdd/6QNriQjbzubhnV/PgWsNRB2DhdNEZLT2XigAAIyO9G7XLpSpC0YBuHs7cWO5zz6d3TBzWBcs25uG+E0p6NvZFb4udm2+v5aYyu+FyNhdK6mBtUSEeqUaS5+Mxqy1yQDM87l/t/18g8M9sfTJGMxcfRqrT2Qj0E2GGUNCBYiQyLLo9O5JIpG0+WJlxbySiHRz5GoxAGBwl/Zd1ZNJrRovktuOSbTHjeU+AeDlYWGICnBBRa0Cr60/A5Wq/T4dN6XfC5Ex++x/l1GvVGNAqDtGdPfSHre05/6Ynj54d1wPAMCHOy9i57l8gSMiMn86nV0CAgJMZhYMEZm2wopaXCmsgkgExIa4Cx2O0bCWiPHvKdEYs/QQjqUX4/sjGXh+UOvH4BBRxzqXU45NjaMW5o3ubvHvo54eEIT0G1VYeSwLs9cl4WeXAejpzzmORO2FDVqMCDfBE93ya3IOZq1NRg9fJ2x7ZZDQ4Ridn05k4e1N5yC1EmPrzAfQ1dtR6JBMHs/BzePvpe3UajWeWnECR64U49FoXyx9MkbokIyCQqnC8/+XiAOXbsDL0Qa/zhwIH+f2K0knMlVs0EJEZuvIlSIAwMAwD4EjMU7T+gViaFdPyBUqzF6XjDqFsuUfIqIO9XtaEY5cKYZUIsbrI7sKHY7RsJKI8cXUGHTt5IjCyjo890MiqutMf/wEkTFiskdERulo4369/qEs4WyOSCTCh4/3gqvMGhfyKvDvPWlCh0REt1Gq1Fi8/QIA4K/9OyPATSZwRMbF0dYaK57pCw8HG6TmVWDW2iQo23EPMpGl0mnPXnZ2NgDAz88PEolE+31rBQYG6nR7IrJM2cU1uF56E1ZiEfoFuQkdjtHycrTF4kk9MWPVafzn4FUM6+aF+/j7IjIKm5JycDG/Ek62Vpg5LEzocIySv6sM3/61D5785jj2XCjEou0X8M4jEUKHRWRWdEr2goKCIBaLkZqaivDwcAQFBbV6o7FIJIJCwSV6ImrZ0asNJZzRAS6wt7GMLnVt9XCkDx7r7Y9fTl/HnPXJ2DFrMBz4OyMSVG29Ep/uvgQAiBsaBheZVOCIjFdMoCs+nRyFmauTsOJwBoI97PFUbGehwyIyGzq9Ixg8eDBEIhFkMtkd35NuEhISkJCQ0OS4Usk9N0TArZELA7hfr1XeHR+B4+nFuFZyE+9vTcWHj/cSOiQyQXxtMpxVx7OQV14LPxc7PD0gSOhwjN4jvXyRWVSNT3ZfxrtbziPU04El/EQGwm6cRoQdz4gautfd98EeFFXJsfaFWI5daKXj6cWY+u1xqNXAN3/pg5E9vIUOyeTwHNw8/l50o1SpMfij/cgpu4klk3riyX7cwtIaarUar65LxubkXLjZS7H15Qfg58IOnWTZ2I2TiMxOWmEViqrksLUWIybQRehwTEZsiDteaJy3N29jCm5U1gkcEZFl2nOhADllN+Eqs8aEGD+hwzEZIpEIiyf1Qg9fJ5RUy/Hij4moreeqMpG+DJ7slZaWorS0FFwwJKK20IxcuC/IDTZWEoGjMS1zRoajm7cjiqvlmLfxLM/DRAL44UgmAGBqv0DYWvMc1lo1cgW6z9+J87kVcJVZ41xOBeZtTOF5jEhPBkn2tmzZgpEjR8LBwQEeHh7w8PCAo6MjRo4ciV9//dUQD0FEFoIjF9rOxkqCz6dEQyoRY8+FQqxPvCZ0SEQW5WJ+BY6lF0MiFrHJiB4+nxINiViETUk5+L4xeSaittEr2VOr1Xj22WcxceJE7NmzBzU1NXB2doazszNqamqwZ88eTJo0Cc888ww/mSGiFimUKhxPb0j2BoayOUtbdPdxwpyR4QCABVtTkV1cI3BERJZj5dFMAMCoHp3gy/1mrVIjVzRebpVsRvo545+jGobQL9p+AUcbKz6ISHd6JXtLly7FDz/8AB8fHyxfvhxlZWUoKSlBSUkJysvL8Z///Ac+Pj748ccfsXTpUkPFTERm6nxuBSprFXC0tUKkn7PQ4Zisvw8KQb8gN1TLlXjt52QOKibqAKXVcmxKygEAPDMgWOBoTEfE/F2ImL8LfRfu0R7ru3APFu+4CKCh4U3c6tO4VsIProjaQq9k75tvvoFMJsOhQ4fw4osv3tElxtHRES+88AIOHToEOzs7fPPNN3oHS0Tm7UjjfL3YEHdIxBzr0lYSsQifTo6CvVSCk5ml+PZQutAhEZm9dYnXUFuvQoSPE+4LchU6HLPR088ZpTX1GPTRfgTN3YYaOWc2E+lCr8m7GRkZGDlyJIKD7/4JVnBwMIYPH47du3fr81BEZAGOaebrcb+e3gLcZJg/LgJv/pKCz3ZfxpBwT3T3Ydt8ovagUKrw47EsAMAzA4M4g1gHqQtGAQBq5Ert6l5i/AjIpA3Nbcpq6vHIF4dRUi0HAG4LItKRXit7np6ekEqlLd7O2toaHh7cf0NEd1enUOJkZgkAYCCHqRvE5L4BGNHdC3KlCq+uS0adgm3MidqDZtyCm70U46N8hQ7HpMikVo0XyW3HJJBJG9YjXGTWWDKpp/a6749kavf5EVHL9Er2Jk6ciH379qG0tPSutykpKcG+ffswYcIEfR6KiMxcUnYZautV8HCwQRcvB6HDMQuauVXu9lJczK/E5/9LEzokIrP0X+24hQCOWzAgzX6+F348pT326e7L2uNE1DK9kr2FCxciJCQEw4YNw759+5pcv3//fjz00EMIDQ3FokWL9HkoIjJzR28r4WQJlOF4OtpgUeOn4l//flW7ekpEhpGaW4ETGSUct6AnmdQKmUvGInPJWO2qHhHpT6f/TcOGDWtyTCqV4tSpU3jooYfg5uaGzp0bTnTZ2dkoLm548xYbG4sJEyZg7969BgiZiMyRprU29+sZ3qge3ni8jz82nLqOOeuTsWPWYDjY8M0UkSFoxi08HOkNH2eOWzCk5vbzBXvYI6OoGv2CXKFQqmAlMcjI6DvUyBXalcPUBaOYfJJJ0+nZe+DAgbtep1arUVxcrE3wbnfs2DF+Uk9Ed1Vdp0DytTIA3K/XXt4dF4FjV4txreQmFv6WiiWP9RI6JCKTV1Itx+bkhnELfxsQJGwwZqi5JOuLqdGY8vVx/JFZio93X8K80d0FiIzIdOiU7GVkZLRXHERkwf7ILIFCpYa/qx0C3GRCh2OWHG2t8enkKEz99jjWnryGhyI6YXj3TkKHRWTS1p7MRp1ChUg/J/TpzHELHSHE0wEfPR6FuNWn8fXBdMQEuOLhSG+D3Lem6cvtA95v/5orfGSKdHrWako0iYgMSTNyYWAoV/XaU2yIO55/IBjfHsrAm7+kYNdsF7g72AgdFpFJun3cwtP9OW6hPWn282mM7eWD09nBWHE4A2/8fAZdvR0R7GGv9+M01/Tl9mHvt8dAZCoMX+hMRKSjI5r9emHcr9feXhvZFV07OaKoqg5vbzrHmVVEbbQ7tQB55bVws5diHMctdLi5o7vhviBXVNYp8I9VpziKgeguDLoeXVZWhsrKyru+eQgMDDTkwxGRGSitliM1rwIA0J/NWdqdrbUEn02JwoSEI9h5Ph8bT+fgsT7+QodFZHLWJ14DADx5H8ctCMFaIkbCtN4Ys+wwLuZX4u1N5/DZ5Ci9VlhbGvBOZIr0XtnLz8/H888/Dy8vL7i7uyMoKAjBwcFNLiEhIYaIl4jMzPH0YqjVQBcvB3g52godjkXo4euM2SPCAQD/2nIeOWU3BY6IyLQUV9XhUFpDRcKk3vywRCheTrZImBYDiViETUk5WP1Htl73d68B79yvR6ZKr2duXl4e7rvvPuTm5sLPzw+enp4oLCxE//79kZ6ejoKCAohEIvTv3x/W1taGitnkJSQkICEhoclxpVLZzK2JzNsfjXPfuKrXsV4cHIK9FwpwOrsMr68/g5+evx9iMfccWTK+NrXe9nP5UKrUiPRzQpiXg9DhWLT7Q9zxz1FdsXjHRby3JRVR/i6I9HMWOiwio6H3UPXc3FwsWLAA165dw+jRoyESiXDkyBHk5eXhwIED6NatG0QiEXbs2GGomE1eXFwcUlNTm1xOnjwpdGhEHS7lejkAIDrAReBILIuVRIzPJkfDzlqCY+nF+G/jrDCyXHxtar0tjeMWHo3yEzgSAoAXBodgRPdOkCtVeOmn0yi/Wa/X/XHAO5kTvZK9nTt3Ijg4GPHx8c1eP3jwYOzevRtJSUl4//339XkoIjJDSpUa53Mb9uv18ucnsR0tyMMeb49tmFH14c6LSCuoFDgiIuN3vbQGJzNLIRIBj0T5CB0OARCJRPj0iSj4u9ohu6QGb/x8hs2niBrplezl5OQgOjpa+71E0lDjXFdXpz3m5+eHoUOHYv369fo8FBGZoas3qnCzXgl7qQTBHiyFEsL0+wMxJNwTcoUKr65PRr1SJXRIREZt65k8AMD9wW7wcbYTOBrScJZZ46vpvSGViLE7tQArDnM2NBGgZ7Ln5OR0x/cuLg1lWDk5OXcct7W1bXKMiOhsYwlnDz9nSLhfTBAikQgfPd4LLjJrnMupwBd704QOiciobTmTCwB4NJolnMaml78L3nmkoVphyY6LOJVVInBERMLTK9kLDAxEdvatzkeRkZEAgO3bt2uP1dTU4MiRI/DxYakDEd0p5XoZAKAnN9MLqpOTLRZOaDh/Jxy4iqTsUoEjIjJOlwsqcSGvAtYSEUZHegsdDjXjqdjOGBflC4VKjZmrk1BSLW/2djVyBYLmbkPQ3G2c0UdmTa9kb9iwYTh79ixu3LgBABg/fjzs7e3xxhtvYO7cufjiiy8wdOhQFBQUYPTo0QYJmIjMR0pOw8oe9+sJ75Fevng02hdKlRpz1p/BTTk7MBL92ZbkhlW9IeGecJFJBY6GmiMSibB4Uk+EeNojr7wWs9clQ6Xi/j2yXHole9OnT8ekSZOQmpoKAHBzc8PXX38NtVqNjz76CLNnz8bJkycRERGBDz74wCABE5F5UChV2uYsXNkzDgvGR8LbyRYZRdVYvOOC0OEQGRW1Wq0t4RzPEk6j5mBjheXT+8DWWozfL99Awv4r2utq5IrGi/K2Y0rtcSJzI1K3Q7ui7OxsbN++HaWlpQgPD8f48eM5Z68VKioq4OzsjPLy8ib7IYnMzYW8CoxeegiONlY48+5IzngzEofSbuAvK/4AAPzfs/0wONxT4Ig6Ds/BzePvpUFSdikmfnUUMqkEifEj2JLfBGw4dR2v/3wGYhGw6rn7MSDMA0Fzt93zZzKXjO2g6IhaZojzr14re3cTGBiIGTNmYN68eXjssceY6BFREyna5ixOTPSMyKAunni6f2cAwBsbzqC8Rr95VUTm4tfGEs6REZ2Y6JmIx/v4Y0rfAKjUwCtrk1FYWSt0SEQdzuBnq9LSho39Li4uEIn4Bo6Imndrvx6HqRubuaO741BaEdKLqvHOr+ewbGqM0CERCUqhVOG3sw0jF9iF07T8a3wPJF8rw6WCSsxak4yUf42ERCxCjVyJvgv3AEDjSq1E4EiJ2odBVva2bNmCkSNHwsHBAR4eHvDw8ICjoyNGjhyJX3/91RAPQURm5mxjssf9esbHTirBZ1OiIRGLsOVMLrY27lMislTH0otRVFUHV5k1HujiwU6ORqSlv4WdVIKE6b0hk0pwLL0Y3x7KgExqdUdyJ5NKGo9xxZbMj17JnlqtxrPPPouJEydiz549qKmpgbOzM5ydnVFTU4M9e/Zg0qRJeOaZZ9AOWwOJyETJFSpcyGtozsJOnMYpOsAFcQ+GAgDiN59DQQXLn8hyaUo4x/T0gbWkXXbAUDsK83LAook9AQBf7EvD4bQigSMi6jh6nbGWLl2KH374AT4+Pli+fDnKyspQUlKCkpISlJeX4z//+Q98fHzw448/YunSpYaKmYhM3OWCSsgVKjjZWiHQTSZ0OHQXLw/vgp5+zii/WY83Npzlh3ZkkWrrldh1Lh8AMKqHNzs5Ggldu2pOiPHD1H4BUKuB2euSUFWrQOaSschcMpYremTW9OrGGRERgezsbKSkpCA4OLjZ22RkZKBnz54IDAzUjmig5rHjGVmKtX9kY+7GFAwMc8dPz8cKHQ7dw5XCSoxZdhhyhQoLJ0TiqdjOQofUbngObp6l/152nsvDjFWn4etsi9zye69ws5Njx2lLV83aeiUmfnUUF/IqcH+wG356/n5YcaWWjJjg3TgzMjIwfPjwuyZ6ABAcHIzhw4cjIyNDn4ciIjNya78em7MYuzAvR7z5cDcAwAfbLiCjqFrgiIg6lqaEc1y0r8CRkL5srSVImBYDe6kEJzJK8O89aUKHRNTu9Fq39vT0hFQqbfF21tbW8PDw0OehiMiMaMYucL+eafjbgCDsSS3AsfRivLY+Getf7M9Pw8kiVNTWY+/FQgDAo1F+mDW8CwCwk6MRSF0wCoDuf4sQTwcsfqwXXlmThIQDV3BfsBuGWNA8UbI8er1aT5w4Efv27dOOW2hOSUkJ9u3bhwkTJujzUERkJuoUSlzMb2jOwk6cpkEsFuGTyVFwtLHC6ewyfP17utAhEXWIXefyIVeo0MXLAd19HLUdG9nJUXj6/C3GR/li+v2BUKuBV9clI7+F8lwiU6ZXsrdw4UKEhIRg2LBh2LdvX5Pr9+/fj4ceegihoaFYtGiRPg9FRGbicn4V6pVquMis4e9qJ3Q41Ep+Lnb41/geAIDP/3cZ5xpLcYnM2c7Gxizjonw5O9jMvPNIBCJ8nFBSLccra5KgUKqEDomoXej0MdSwYcOaHJNKpTh16hQeeughuLm5oXPnhs372dnZKC4uBgDExsZiwoQJ2Lt3rwFCJiJTdjanDEDDqh7fPJmWSb39sDs1H7vLCDOLAAAgAElEQVTOF2DO+mRsmfkAbK1ZvkbmqbZeiSNXG1r0j+zR6Y7rZFIrNmMxEm39W9haN8zfG/fFYfyR2bB/7/VRXdshQiJh6ZTsHThw4K7XqdVqFBcXaxO82x07doxv6ogIAPfrmTKRSIRFE3viVFYZLhdU4dPdl/D22AihwyJqF8euFqO2XgVfZ1t07eQodDjUDoI97LF4Uk+8zP17ZMZ0SvbYUZOI9HX2OjtxmjJ3BxssmdQTz/9fIr47nIHh3TshNsRd6LCIDG7vxQIAwLDuXvzA2oyNi/LF8fRi/HQiG3PWJWP7rEHo5GQrdFhEBqNTsqcp0SQiaovaeiUuF1QCAHpyZc9kjYjohCl9A7Au8RpeW38GO2cPgqOttdBhERmMWq3G/os3AADDunkJHA21t3ceicDp7DJcyKvAK2uSOH+PzAqfyUTUYS7mV0KhUsPdXgpfZ35yasreGRcBf1c75JTdxPu/pQodDpFBXSqoRE7ZTdhaizEglKOjzJ2ttQRfTe+tnb+3dC/n75H5MEiyV1BQgMWLF2PMmDGIiopCVFQUxowZgyVLlqCgoMAQD0FEZiDlemNzFn82ZzF1DjZW+GxyNEQiYH3idew+ny90SEQGs69xtt6AUA82IbIQwR72WDSpJwDgy/1XcCjthsARERmG3sneL7/8gvDwcMTHx2Pnzp1ISUlBSkoKdu7cibfffhtdu3bFL7/8YohYicjEafbr9eJ8PbPQL9gNLwwKAQDM25iCoqo6gSMiMox9FxqSPZZwWpZHo/0wrXH+3uy1ySio4Pw9Mn16JXuJiYmYOnUqqqurMXHiRGzatAlJSUlITk7G5s2bMWnSJFRVVWHatGlITEw0VMxEZKJSGmezRTLZMxtzRoajm7cjiqvlmLcxBWq1WuiQiFqtRq5A0NxtCJq7DTVyBQCgtFqO09mlAIChTPYszvxHIrTntFfWJEGp4jmNTJtODVr+bPHixVAqldiwYQMmTpx4x3W9evXC+PHjsWnTJjz22GNYsmQJNmzYoFew5iIhIQEJCQlNjiuVSgGiIeoYN+VKpBVWAQB6+bMTp7mwsZLgs8nReDThMP6XWoANp67jib4BQodFbcDXpgYHL9+ASg1083aEn4ud0OFQB9PM3xv/xWHt/r05D4ULHRZRm4nUenwM26lTJ4SHh+PQoUP3vN2gQYNw+fJl7t9rQUVFBZydnVFeXg4nJyehwyEyqFNZpXhs+VF4Otrgj7eGc8+emfnqwBV8tPMSHGyssGPWIAS4yYQOSWc8BzfPHH8vmlW8GrkSfRfuAQAkxo+ATCrBGz+fwbaUfMQNDcUbo7oJGSYJ6NfkHMxamwyRCFj13P0YGMZGPdTxDHH+1Wtlr7y8HIGBgS3eLjAwECdPntTnoYjIxGmas/TyY3MWc/Ti4FDsvVCIU1mleP3nM1jz91iIxfw7k3GKmL+ryTFN0qfB/XqW7dFoPxxPL8aaP65h1tpkbJ/1ALwc2UWaTI9ee/a8vb2RlJTU4u2Sk5Ph7e2tz0MRkYlLyakAwPl65koiFuGzyVGQNbYu//5IhtAhEbWZq8wa0QGuQodBAnt3XA9083ZEUVUdZq9N5v49Mkl6JXujRo3CpUuX8NZbbzVb069WqxEfH4+LFy/i4Ycf1uehiMjEpeQ0jl1gcxaz1dndHvFjIwAAH+26hMsFlQJHRNS81AWjkLpgFBLjR2iPJcaPwHMPBAEAHuzqBQlXpi2erbUEX07rDZlUgqNXi/HlvitCh0SkM72SvXfeeQdubm748MMPERYWhjfffBPLly/H8uXLMXfuXISFhWHx4sVwd3dHfHy8oWImIhNTXafAlcbmLEz2zNvUfgEY2tUTcoUKr65LhlyhEjokoiZkUqvGi+S2YxL8frkIAEs46ZYwLwcsnBAJAFi69zKOXi0SOCIi3ei1Z8/f3x/79u3D9OnTce7cOXz88cfavTiavi89e/bETz/9BH9/f/2jJSKTlJpXAZUa8HayhZcT9zyYM5FIhA8f64VR//4d53MrsGxvGl4f1VXosIhadL20BmmFVZCIRRgc7il0OGREJvX2x/H0YqxPvN6wf++VQfB0tBE6LKJW0SvZAxqSubNnz+LAgQM4dOgQcnNzAQC+vr4YNGgQHnzwQX0fgohMXErjMHXu17MMXk62+GBiT7z002l8deAKhnbzQp/O3P9ExkcmtULmkrEAgJVHMwEAfTu7wtnOWsCoyBi9Nz4SydfKcLmgCnPWJ2Pl3/qxCRWZBL2SvUmTJsHHxwcJCQl48MEHmdgRUbNS8xqas/TwNY+27dSyMT19MCHaF5uTc/Ha+mRsnzUIMqneny8StZt9FwsBsISTmmcnlSBhWm+M//IIDqUV4asDVzBzWBehwyJqkV579rZv347i4mJDxUJEZkozTL1rJ0eBI6GO9N6jkfBxtkVmcQ0Wbb8gdDhEd1Vdp8Cxqw3vZ4Z3Z7JHzevSyRELHu0BAPjsf5dxIp3vgcn46ZXsBQcHo7q62lCxEJEZUqvVuNLYlbFLJweBo6GO5GxnjU+eiAIArDqejYOXbwgcEVHzjlwpglypQqCbDKGePE/R3T3RNwCTevtBpQZeWZuE4qo6oUMiuie9kr2pU6fi4MGDyM/PN1Q8RGRm8sprUS1XwkosQmd3e6HDoQ42MMwDzwwIAgC88fMZlNXIhQ2IqBn7L90q4dQ0miO6m/cfjUSopz0KKuowZ/0ZqDh/j4yYXsnevHnzMGjQIAwZMgSbNm1CfX29oeIiIjOhKeEM9rCHtUSvUw6ZqDcf7oYQT3sUVtYhfvM5ocMhuoNareZ+PdKJvY0VEqb3ho2VGAcv38DXv6cLHRLRXen1zqtr1644f/48rly5gscffxx2dnbw9fVFSEhIk0toaKihYiYiE5LGEk6LZyeV4PPJ0ZCIRfjtbB62nMkVOiQirfO5FSioqINMKsH9IW5Ch0Mmopu3E94b37B/75Pdl5CYWSJwRETN06s1WmZm5h3fq9VqlnQS0R00w9TDvNicxZJFBbhg5tAwLN2bhvhNKegX5AZvZ85cJOFpVvUeCPOAjZWkhVsT3TLlvgAcSy/Gr8m5eHlNEra/Mgiu9lKhwyK6g14reyqVSqcLEVkeTRlnFy+u7Fm6mcPC0MvfGRW1Cryx4QzUau5zIeEdSmtoHDSUJZykI5FIhA8m9kSIhz3yymvx+s88r5Hx4QYaImo3arWaZZykZS0R47PJ0bCxEuNQWhFWHc8SOiSycDflSiRfKwMADAh1FzgaMkUONlb4clpvSK3E2HuxEN8dyhA6JKI7tCnZ2759O1544QWMHj0aEyZMwPz585GRwSc3Ed3pRmUdKmoVEIsaGrQQhXk5YO7obgCAD7ZfQPqNKoEjIkuWmFWCeqUavs62CHSTCR0OmagIXyfMfyQCAPDhzos4nV0qcEREt+ic7E2fPh3jxo3DihUrsGvXLmzZsgUffPABevTogS1btrRHjERkojQlnEHu9twLQ1pP9w/CwDB31NarMGf9GSiULPMnYWgGqceGunPkArVKjVyBoLnbEDR3G2rkCu3x6fcH4uFIbyhUakz66ijyym8KGCXRLToleytWrMCaNWsgkUjwzDPPYNmyZfjggw8QGxuL2tpa/PWvf0V5eXl7xUpEJkZTwhnG/Xp0G7FYhI8fj4KjrRWSr5Vh+YGrQodEFupYekOy1z+EJZykH5FIhAWN3TkB4O1N57h/j4yCTsneypUrIRaLsWPHDqxYsQIzZ87EvHnzcOTIETz99NOorKzExo0b2ytWIjIx2uYs3K9Hf+LrYocFjza8MVq6Nw3ncvhBIXWsqjoFzl5veN715349akGNXNF4Ud52TIkauQJFVbWokSsgFt9aHd53sRDf/J5+x+ofkRB0Gr2QkpKC2NhYDB8+vMl1b731FlauXImUlBSDBUdEpu1WJ06OXaCmJkT7Yff5Auw4l49X1yVj68sPwNaa5b7UMU5mlkCpUiPAzQ7+rtyvR/cWMX9Xk2N9F+65588s3nERi3dcROaSse0VFlGLdFrZq6iouOtwdM3xiooK/aMiIrNwa8YeV/aoKU3bcg8HG6QVVuHjXZeEDoksyPGrLOGkjlF+s17oEMiC6ZTsqdVqSCTNf+oqFjfcFefpEREAFFfVoaRaDpEICPVkskfNc7OX4qPHewIAVhzOwNGrRQJHRJZCu1+PJZzUCqkLRiF1wSgkxo/QHkuMH9F4bHiT6/a/PgT+rnYAgLm/nOX+PRIM5+wRUbvQlHAGuMpgJ2VpHt3dsG6dMLVfAADgjZ/PoqKWn4JT+6qordfuE+0f4iFwNGQKZFKrxovktmMSyKRW8HCwbXJdJydbJEzrDWuJCDvO5eNHzhUlgeic7K1cuRISiaTZi0gkuuv1VlY6bQ8kIhOnSfbC2ZyFWiF+bAQC3WTIKbuJBVtThQ6HzNzJjBKo1A3zP72dbYUOh8xUVIAL5o7uDgBY+NsFNqIiQeic7KnV6jZdWN5JZFmuaMcusDkLtczexgqfTo6CSARsOHUdu87nCx0SmTHtfD3u1yMdyaRWyFwyFplLxkImtWrxumcHBmFE906QK1WIW30alaxcoA6mU7KnUqn0uhCR5bjViZMre9Q69wW54cXBDc2+5m1MwY3KOoEjInOl2a8XG+ImcCRk7kQiET55ohf8XOyQVVyDeRtTuH+POhT37BFRu7hcwBl7pLtXH+qCbt6OKKmWY95GNjUgwyurkSM1r6FzODtxUkdwkUmxbGoMrMQi/HY2D6v/yBY6JLIgTPaIyOBKq+UoqmpYlWEnTtKFjZUEn0+JhlQixp4Lhfg58brQIZGZOZFRArUaCPW0h5cT9+tRx+jT2RVvjOoKAHhvaypSczmqjDoGkz0iMrgrNxpW9fxc7GBvw+ZMpJvuPk6YMzIcAPDe1vO4VlIjcERkymrkCgTN3YagudtQI1do9+tx5IL5+/PfXmh/HxSCoV09IVeoMHP1aVTVCR8TmT8me0RkcGks4SQ9/X1QCO4LckW1XInXfj4DpYrlnGQYxzXz9ThygTqYWCzCp5Oj4e1ki/SiasRv4v49an9M9ojI4NIKGzpxsjkLtZVELMKnT0RDJpXgj4wSfH84Q+iQyMTUyBWNF6X2WE7pTVzMbzg/sTmL+Wrub18jV2qPC8nNXoovpsVAIhZhc3IuS9Wp3bG+iogM7oq2EyfHLlDbBbrL8M4jEZi3MQUf77qEweGe6OrN5xS1TsT8XU2OPfT579qv3R1sOjIc6kDN/e37Ltyj/TpzydiODKeJ+4LcMOehcHy86xLmbzmH6EAXhHfiuY3aB1f2iMjgNGWcYSzjJD09eV8Ahnfzglypwux1yZArOMaHiEzfP4aEYnC4J2rrVXjpp9OCrziS+WKyR0QGVVFbj/yKWgBAGMs4SU8ikQiLH+sJV5k1LuRVYOney0KHRCYidcEopC4YhcT4EdpjwR72AIBlT0YLFRZ1gOb+9onxI7THjYFYLMJnk6Pg5WiDK4VVmP/reaFDIjPFZI+IDEpTwuntZAsnW2uBoyFz4OVoi0UTewIAlh+4ilNZJQJHRKZAJrVqvEi0xzKKqiESAYPDPQWMjNpbc397mVSiPW4sPBxssGxqDMQiYMOp6/jlFPfvkeEx2SMig7rCTpzUDkb39MGkGD+o1MCc9WdQzZbl1EbdvZ3gIpMKHQYRACA2xB2zRzSMmonffA5XGhucERkKk71WWrx4Mfr27QtHR0d06tQJkydPRmZmptBhERkdTSdOlnCSob07vgd8nG2RVVyDRdsvCB0OmQiZ1AqZS8Ziar9AAJyvZ0k0f/vMJWONakXvz+KGhmFgmDtu1isxc3USauuVLf8QUSsx2WulgwcP4uWXX8aJEyewc+dOlJSUYPTo0VAo+Oky0e3S2ImT2omznTU+eSIKAPDTiWzsv1QocERkSm7N12OyR8ZFIhbh8ynR8HCwwcX8Sry3lfv3yHCY7LXSzp078fTTTyMiIgIxMTH49ttvcfHiRaSmpgodGpFR4UB1ak8Dwzzwt4FBAIA3N5xFabVc2IDIJOSX1yKjqBpiEdCP8/XICHk52mLpk9EQiYA1f1zDr8k5QodEZsKok71Vq1bhxRdfRN++fWFjYwORSIQffvjhnj9z8uRJjBkzBi4uLrC3t0dsbCzWr19v8NjKy8sBAG5ufNEg0qiuUyCn7CYAIMyTyR61jzcf7oZQT3sUVtYh/tdzUKvVQodERu5YehEAINLPmY2jyGgNDPPAy0PDAABvbUxB+o0qgSMic2DUyV58fDy++eYbZGVlwcfHp8Xb79+/HwMHDsThw4cxefJkzJgxA/n5+ZgyZQo+/fRTg8WlVCrx+uuvY8yYMfD39zfY/RKZuquNL0weDjZwtWcDBGofttYSfD4lGlZiEbadzcOWM7lCh0RG7vjVhg6uLOEkYzdrRDjuD3ZDtVyJOO7fIwMw6mTvu+++Q2ZmJm7cuIEZM2bc87YKhQJ///vfIRaL8fvvv+Obb77Bp59+ijNnziA8PBxvvfUWsrKy7viZuXPnQiQS3fPyZ2q1GjNmzEB2dnaLq4xElkZbwsnmLNTOevm74OVhXQAA72w+h7zymwJHRMbsZOO4jn7BrMYh4yYRi7Bsagzc7KW4kFeBD7axGRXpx6iTvREjRqBz586tuu2+fftw9epVTJs2DdHRt4alOjs746233oJcLsfKlSvv+JnXXnsNFy5cuOfldmq1Gi+99BL27NmDvXv3wtOTc3qIbqdtzsL9etQBXhoaiih/Z1TUKvDPDWehUrGck5oqq5Ej/UY1ACAm0FXgaIha1snJFp9NbmhG9ePxLGxPyRM4IjJlxtuHVkcHDhwAAIwcObLJdaNGjQLQ0FHzdp6enq1O2NRqNeLi4rBt2zYcPHgQAQEB+gVMZIY084G4skcdwVoixmdTojF22SEcSivCj8ez8PSAIKHDIiOTdK0MABDsYQ83HcrLa+QKRMzfBQBIXTDKqFv3k/HT9fn0YFcv/OPBUCw/cBUv/XS61T9H9Gdm84xJS0sDAHTp0qXJdd7e3nBwcNDepi3i4uKwZs0abN26FXZ2dsjPzwfQ0KBFKm364lFXV4e6ujqdHqOioqLN8REZA83KXhjHLlAHCfV0wLzR3fHulvNYvOMCHujigVA2B6LbJGWVAgBiAl0EjoRIN689FI4T6cU4nd3wgYVcoYKM2+FJR0ZdxqkLTXdMZ2fnZq93cnLS3qYtli9fjrKyMgwaNAg+Pj7ay9GjR5u9/eLFi+Hs7KzThauFZMpq65XILqkBwDJO6lh/ie2MB8I8UFuvwpz1Z6BQqoQOiYyI5o1y71aWcNbIFY0X5W3HlNrjRLpo6/OpRq6AXKnCgkcjtceW7LjI5yHpzGxW9tqbrq29582bhzlz5uj0MxUVFUz4yGRdvVEFtRpwlVnDnZ04qQOJxSJ8/EQvjPr8d5y5VoaE/Vcxa0TTKg+yPEqVGsnXdEv2NKV2t+u7cI/268wlYw0THFmEtj6fmvu5tSevYe3Ja/f8OaI/M5tkT7Oid7fVu4qKCri6dtzGbBsbG9jY2HTY4xEJ7YqmOYuXY7OdbInak4+zHd6fEIlZa5OxbF8ahnbzRC9/lu1ZurTCSlTVKSCTStDVm+XlRGR5zCbZ0+zVS0tLQ58+fe64Lj8/H1VVVejXr58QoRFZBM3YhTCWcJJAxkf5YndqAbadzcOr65Kx7ZVBsLWWCB0WCeh0VsOqXpS/CyTi1n0IlbqgoalbjVypXYFJjB8BmZTPJdJdW59Pzf1cpJ8TzuVUoKefM+QKFaRWZrMbi9qR2TxLhgwZAgDYvXt3k+t27dp1x22IyPDSixqSPTbHIKGIRCIsfDQSXo42uHqjGh/uvCh0SCSwpOyG5iy9O7d+lVcmtWq8SG47JtEeJ9JFW59Pzf3c55Oj4WRrhZSccnzE8xu1ktkke8OHD0dISAhWr16N5ORk7fHy8nIsWrQIUqkUf/3rXwWMkMi8ZRY1NGcJcpcJHAlZMld7KT58vBcA4L9HMnH0SpHAEZGQTmuSPc7XIzPg52qHj59omL/33eEM7EktEDgiMgVG/RHVd999h8OHDwMAUlJStMc0M/UeeOABPP/88wAAKysrfPfddxg1ahQGDx6MJ598Eo6Ojvjll1+QlZWFTz75BEFBQUL8M4jMnlqt1nbi7MxkjwQ2tKsXpt0fiNUnsvH6z2ewY/ZgONtZCx0WdbCyGjmu6jFMXSa1YhMMMpi2Pp/+/HOjenjjbwOD8N8jmXh9wxlsf2UQfF3sDBkqmRmjTvYOHz6MlStX3nHsyJEjOHLkiPZ7TbIHAEOHDsXhw4fx7rvvYt26daivr0fPnj3x4YcfYsqUKR0Wd0sSEhKQkJDQ5LhSqWzm1kTGr6Rajqo6BUQiwN+VyR4J7+0x3XHkShGyimvw3tbz+GxytNAhGT1ze21q6zB1ImM3b3R3nMoqxdnr5Xh5TRLWvhALa4nZFOuRgYnUus4UoHZTUVEBZ2dnlJeXw8nJSehwiFrtdHYpJn11FL7Otjg6b7jQ4RABAE5lleCJ/xyDSg3856neeDjS55635zm4eab6e/ls9yUs23cFk3r7Mdkns5NdXIOxyw6hsk6BGUNCMXd0N6FDonZgiPMvPwYgIr1lFzeUcAayhJOMSJ/ObpgxJBQAMG9jCgorawWOiDqSZmWvLSWcRMYu0F2m3Z/8n4NXsf9SocARkbFiskdEesssbtgX09nNXuBIiO40e0Q4uvs4obSmHvN+SQGLWSyDSqVGcrZmmDrnLZJ5GtPTB3+J7QwAeG39GeSX8wMtaorJHhHpjSt7ZKykVmL8e0o0pBIx9l4sxLqT14QOiTpAWmEVKjXD1DtxmDqZr7fHdkeEjxNKquV4ZU0SFEqV0CGRkWGyR0R6y2InTjJiXb0d8fqocADA+7+laj+cIPOlGbkQ5e8CKzauIDNmay1BwvTesJdK8EdmCZbuTRM6JDIyPAMSkd6yijUz9ljGScbpuQdC0C/IDdVyJd7enCJ0ONTOTmfpPkydyFQFe9hj8WMN+/e+3H8Fh9M4X5RuYbJHRHqpqlOgqKoOAMs4yXhJxCJ8OjkKA0Ld8e64HkKHQ+1Ms7IXE8DmLGQZxkf5Ymq/QKjVwOx1SSis4P49amDUc/bMlbnNMiLLpimJc5VZw8mWg6vJeAW4ybD677FCh2G0zOW1qbym/rZh6lzZI8vx7rgIJGWX4mJ+JWatTcaq5++HRCwSOiwSGJM9AcTFxSEuLq7Jcc0sDSJTkl3S8KYqkCWcRCbNXF6bkq41rOoFucvg7mAjcDREHUezf2/cF4dxLL0YX+xLw+wR4UKHRQJjGScR6eXWfj2WcBJR+6mRKxA0dxuC5m5DjVxx19ud1o5cYAknmYbWPrdbI9TTAR9MjAQALN2bhqNXuX/P0jHZIyK9ZDYme53dmOwRkfCSNPv1OjPZI8s0McYfk/v6Q60GZq1N1u6rJ8vEZI+I9MIyTiJqTzVyReNFedsxpfb47ThMnUyJLs9tXf1rfA908XLAjco6vLouGSqVWt9wyURxzx4R6UVTxskZe0TUHiLm72pyrO/CPdqvM5eM1X595QaHqZPp0OW5rSuZ1ApfTe+NcV8exqG0Iiw/eBVxQ8PafH9kuriyR0RtJleokFt2EwDLOIlIeJr5er38nTlMnSxel06OWPBow/69T3dfwh8ZJQJHRELgyh4RtVlO2U2o1ICdtQSejux6R0SGl7pgFICG8jbNqkdi/AjIpJImt9XM12NzFjIFujy32+qJPv44frUYG5Ny8MqaJGyfNQhu9lKD3T8ZP37sRURtllncsF+vs7sMIhFn+RCR4cmkVo0XyW3HJNrjt2MnTjIlujy320okEuH9CZEI9bRHfkUt5qzn/j1Lw5U9AZjL4FoizUD1QJZwEpk8U39tKq+px5XCKgAcpk50O3sbK3w5rTcmJBzBgUs38M2hdMwYEip0WNRBRGq1mum9kdAMri0vL4eTk5PQ4RC1aMHWVHx/JAN/HxSMt8dGCB0OkV54Dm6eqfxeDlwqxDP/PYkgdxkOvDFU6HCIjM7qE9l4a1MKJGIR1r8Yiz6d3YQOiVpgiPMvyziJqM00Yxc6c+wCEQksqbGEM4YlnETNmtovAOOifKFUqfHy6iSUVsuFDok6AJM9ImqzTI5dICIjkXSN8/WI7kUkEmHRxEgEucuQW16LNzacAQv8zB+TPSJqE5VKjeySxmTPjSt7RCQctVqNlOsNyV5UAJM9ortxtLXGl9N6QyoRY8+FQqw4nCF0SNTOmOwRUZsUVNZCrlDBSiyCr4ut0OEQkQXLLa9FaU09rMQihHOYOtE9Rfo5451HugMAPtx5EcmNq+JknpjsEVGbZDWWcPq52nF4MREJ6lxOOYCGIdK21oabUUZkrp6K7YwxPb1Rr1Rj5urTKL9ZL3RI1E74Do2I2iRbu1+PJZxEJCxNstfTz3i7hRIZE5FIhCWP9UKAmx2ul97EP7l/z2wx2SOiNtEOVOeMPSISmCbZi/RzFjgSItPhZGuNhGm9YS0RYdf5AvzfsSyhQ6J2wGSPiNokq4SdOIlIeGq1Gik5FQCY7BHpqpe/C+aNbti/98G2C9oPTsh8WAkdgCVKSEhAQkJCk+NKpVKAaIjaRlPGGciVPSKzYKqvTYWVdSiqqoNYBHT3Zhknka7+NjAIx9OLsTu1AHGrT+O3lx+Ao6210GGRgYjULNA1GhUVFXB2dkZ5eTmcnPiCRcZLrVaj13u7UVmrwK7Zg9HVm93vyPTxHNw8Y/+97L1QgOdWJqJrJ0fsenWw0OEQmaTymnqMWXYIOWU3MbaXD76cGgORSCR0WBbPEOdflnESkc7KaupRWasAwJU9IhJWSmPZWQ82ZyFqM2eZNb6cFgMrsQjbzuZh9R/ZQodEBsJkj4h0ptmv18nJBnZStjknIjBKX7MAACAASURBVOGc0+zX8+V+PSJ9xAS64s2HuwEA3tuaitTcCoEjIkNgskdEOsvSduLk2AUiEpZ27II/kz0ifT0/KBjDu3lBrlBh5urTqKpTCB0S6YnJHhHpTNuchZ04iUhANyrrkF9RC5EIiPBhGSeRvkQiET55Igo+zrZIL6pG/KYUzt8zcUz2iEhnmjLOICZ7RCSg87kNq3ohHvawt2GDcSJDcLWX4oupMZCIRdicnIv1ideEDon0wGSPiHSmKeMMdGcZJxEJh8PUidpH3yA3vDYyHADw7pbzuJRfKXBE1FZM9ohIZ1mNZZyd2YmTiATE5ixE7WfG4FAMCfdEbb0KcatPo0bO/XumiMkeEenkplyJwso6AEBnlnESkYBSuLJH1G7EYhE+mxyFTk42uFJYhfm/nhc6JGoDJntEpJPsxv16TrZWcJFJBY6GiCxVabUcOWU3AXDGHlF7cXewwdInYyAWARtOXccvp64LHRLpiMkeEekks3G/XpAH9+sRkXDONTZnCXKXwcnWWuBoiIRRI1cgaO42BM3dZrAyyz/fZ2yIO2aPaNi/F7/5HK4UGnb/Xnv8G+gWtq4SQEJCAhISEpocVyqVAkRDpBvt2AXu1yMyK6b22qTZr9eDJZxE7S5uaBhOZBTjyJVixP2UhM1xA2EnlQgdFrUCkz0BxMXFIS4ursnxiooKODvzRYuMW1ZJ40B17tcjMium9tqkWdnryWSPLJBmBaxGrrzt2K2vZVLd3+K3dJ+fT4nGmKWHcamgEu9tPY8lj/XS+TF0eby2/BuoKf4WiUgntzpxsoyTiISjHbvATpxkgSLm72pyrO/CPdqvM5eMbZf7XPpkNJ5acQJrT15D/1B3PBrtp/Pj6PJ4pD/u2SMinWgatHBlj4iEUn6zXvvBUw9fNmch6igDwzzw8tAwAMBbG1OQfqNK4IioJVzZI6JWq1eqcL20oftdZw5UJyKBnG8s4fR3tYOrPbsCk+VJXTAKQEPZo2Y1LDF+BGR67KNr7X3OGhGOExklOJFRgrjVSdj00gDYWuv+uO3xb6CmuLJHRK2WW3YTSpUaNlZieDnaCB0OEVmo8xymThZOJrVqvEhuOybRHm/P+5SIRVg2NQbu9lJcyKvA+7+lGs2/gZpiskdErZZ1WydOsVgkcDREZKm0zVn8mewRCaGTky0+mxINAPjpRDZ+O5srcER0N0ybiajVsrhfj4iMQEpjcxbu1yNLJ5NaGbyRSWvvc0i4J156MBRfHbiKub+kINLXuU0zeNvj30C3cGWPiFotq0gzdoH79YhIGFV1CmQ0nosiOXaBSFBzHgrHfUGuqKpTYOaa06hTGOdcTkvGZI+IWu1aKQeqE5GwUnMroFYDPs628HDg3mEiIVlJxFg2NQauMmucy6nA4u0XhQ6J/oTJHhG1Wm5ZLQDAz8VO4EiIyFKd05ZwclWPyBj4ONvh08lRAIAfjmZi57k8gSOi2zHZI6JWyylrGLvg58pkj4iEoW3OwhJOIqMxrFsnvDA4BADwxoazuNa4x5+Ex2SPiFrlplyJkmo5AMCXK3tEJBDNyl6kH5uzEBmTN0Z1RUygCyprFZi5+jTkCpXQIRGY7BFRK+WWN6zqOdhYwcmWjXyJqOPdlCtxpbAKAFf2iIyNtUSML6bGwNnOGmeul+PDndy/ZwyY7BFRq+Q2lnD6uthCJOKMPSLqeKl5FVCpAU9HG3g52QodDhH9ib+rDB8/3gsAsOJwBv6XWiBwRMSP5wWQkJCAhISEJseVSrarJeN1K9ljCSeROTKF16bzjfv1Ijlfj8hojezhjWcHBuP7Ixl4/ecz2D5rEBu7CYjJngDi4uIQFxfX5HhFRQWcnVmWQsYpp7SxOQtP2ERmyRRem1KuszkLkSmYO7obTmWV4Mz1cry8+jTWvdgf1hIWFAqBv3UiapWcxrELXNkjIqGcz60AAERw7AKRUZNaifHltN5wtLXC6ewyfLL7ktAhWSwme0TUKpoyTq7sEZEQFEqVtjlLhA/LOImMXYDbrf17Xx9Mx/6LhQJHZJmY7BFRq2i6cXJlj4iEkFFUDblSBXupBP6c9UlkEh6O9MHT/TsDAOasT0Ze43sJ6jhM9oioRSqVGnnaMk52wCOijncxvxIAEO7tCLGYHYGJTMVbY7sj0s8JpTX1mLUmGQol5+91JCZ7RNSioqo6yJUqiEVAJ7Y7JyIBXGpM9rp5OwocCRHpwsZKgi+n9oaDjRX+yCzBv/ekCR2SRWGyR0Qtymncr+ftZMtuWkQkCM3KXtdOTPaITE2Qhz0WT+oJAEg4cAWH0m4IHJHl4Ls2ImpRLjtxEpHALuY3dOLs6s3mLESmaFyUL6bdHwi1Gpi9NhmFFbVCh2QRmOwRUYs4UJ2IhFRVp8D1xlmfLOMkMl3zH4lAN29HFFfLMWttMpQqtdAhmT0me0TUohwme0QkIM1+vU5ONnC1lwocDRG1la21BAnTe0MmleBYejGW7eX+vfbGZI+IWpSjnbHH5ixE1PE0yR5LOIlMX6inAxZNbNi/t2xfGo5eKRI4IvPGZI+IWqQdqM7ZVkQkgEuN+/VYwklkHibE+GFK3wCo1cCsdcm4UVkndEhmi8keEbWIe/aISEgXOHaByOz8a3wPhHdywI3/b+/Ow6Kq+/6Bvw8DDDvIJqgIyeKK4r4kLt3eYWpqLqG0YWlZdNuT3pVZpmm/8i41e5Tql5pbapp7VtqGK+4LopkLmwsg+77Def7AGeVmBhXhfIeZ9+u6/KNz5sCbc8V8+Mz5LvmlmL75LKo4f69RsNkjojoVlVUgu6gcAJs9IlKeLMt3DeNks0dkLKwtVYgM6wZrCxUOXsnAl/uuio5klNjsEVGdNNsu2KvN4WBlITgNEZmaW3mlyC0uh8pMgp+7neg4RNSA/JvbY96ojgCAxb9dxrH4TMGJjA+bPSKqE4dwEpFImv31HnG1hdpcJTgNETW08T28MKZbS1TJwLTvzyCzgPP3GhKbPSKq000uzkJEAv3NIZxERm/+qE7wdbPFrbxSzPghhvP3GhCbPSKq050ne9x2gYiUp5mv157NHpHRslWbI/KZblCbm2HfpXR8czBedCSjYS46gCmKjIxEZGRkreOVlZUC0hDVjRuqE5kGQ61Nf3OPPSKT0M7DAXNHdsS722Lx2d5L6OnTDN29nUXHavLY7AkQERGBiIiIWsfz8vLg6OgoIBGRfto99tjsERk1Q6xN5ZVViEsrAMBtF4hMwYSeXjgSl4ldMcn414Yz+GlaMJrZWoqO1aRxGCcR1UmzGief7BGR0hIyClFWWQVbSxU/cCIyAZIk4eMxgXjE1RbJuSV4a0sMZJnz9x4Gmz0i0quqSkZKLp/sEZEYdy/OYmYmCU5DREqwU5tjWVhXWJqb4feLaVh5KEF0pCaNzR4R6ZVeUIryShkqMwnu9mrRcYjIxFy6ve0C5+sRmZaOLRwxe0QHAMB/9vyNs9dzBCdqutjsEZFemsVZPBysYK7i2wURKUuzEifn6xGZnmd7t8awQA+UV8p4fcNp5BaVi47UJPGvNyLSi9suEJFIF1O4xx6RqZIkCQvGdkZrZxvcyC7G21s5f68+2OwRkV7J3HaBiATJLynXji7gkz0i0+RgZYFlYV1hoZKw98ItrIlOFB2pyWGzR0R63czm4ixEJMblW9VP9TwcrOBkw6XXiUxV51ZOmDWsPQDg45//RuyNXMGJmhY2e0Sk101uu0BEgty9EicRmbbwfj4I6dgcZZVViNhwGnklnL93v9jsEZFe3FCdiETh4ixEpCFJEj4d2wWtmlnjWlYR3t0ay/l794nNHhHplZzLOXtEpJyisgr4zPwJPjN/woVkzbYL92727r6uqKyisWMS0QNqiN9RRxsLLJ3YFeZmEn6KTcH6Y9ca/PsZ43sJmz0i0qmwtAI5t5c55mqcRKS0K7c0T/a4xx4RVevauhlmPtEOADBv91+4kMz5e/fCZo+IdNIM4XSwMoe9lYXgNERkzIrKKm7/q9QeyyupgMpMqvPDJl3XFZVVao8TkViN8Tv6Uv9HMKS9O8oqqvD6hjMoKL3zder7/Yz5vUSSOeDVYOTl5cHR0RG5ublwcOAnmSTWvktpCF91Au087LHnfwaIjkPU6PgerJsS98Vn5k91nk9cMLxBryMiZTTW72hOURmGfXEQybklGBXUAktCgyBJktG9lzTE+y+f7BGRTsm3V+Lk4ixERERkSJxsLLE0rCtUZhJ2nk3GphPXRUcyWOaiAxCRYeKG6kSklL/mhQCoHjbV46Pftcff+Ic/XhnY5oGuO/n+ENhYqhoxLRHdr8b8He3u7Yy3QtpiwS9/Y86uCwhq7VTv72fM7yV8skdEOrHZIyKl2Fia3/5X8w+rwJaOsLHU/7m0rutsLFXa40QkVmP/jr4c3AaD2rqhtKIKEetPQ5ZRr+9nzO8lbPaISKcb2maPK3ESkRjcUJ2I6mJmJmHR+C5o7qBGXHohZu88LzqSwWnarSoRNRrNk71Wzfhkj4iUYWNpjr3/MwAhSw7ATm1+3+8/NpbmXIyFyIA15u+oi50a/zuhKyYuP4ptp2+ibxsXjO/hVa/vZ4zvJXyyR0S1VFbJSM2tXqCFwziJSEl/p1Zvph7Q3A6SJAlOQ0RNQe82Lpj+zwAAwAc7L2j36SQ2e0SkQ3p+KSqqZKjMJLjbcxgnESnnUurtzdQ9uf0FEd2/1wb5IdjfFcXllYjYcBrFd+2ZZ8rY7BFRLTdvD+H0cLCCyoyfrBORcrTNHufrEdEDMDOTsPjpILjZq3H5VgHm7rogOpJBYLNHRLVomj3usUdESrt0e/hVQHM2e0T0YNzs1fgiNAiSBGw6eR07ztwUHUk4NntEVItmcZaWXJyFiBRUXFap/bDJ391OcBoiaor6+bli2mP+AIBZ22MRl14gOJFYbPaIqJZkbrtARALEpRdAloFmNhZwsVOLjkNETdS0f/ijbxsXFJVVImL9aZSUm+78PTZ7RFQLN1QnIhE0n8D78akeET0ElZmELyYEwdXOEn+n5mPe7r9ERxKG++wJEBkZicjIyFrHKytN91MHMiw3c7jtApGpMYTadDWNzR4RNQx3Byt8HhqE5789jg3HrqFvGxc82aWF6FiKY7MnQEREBCIiImodz8vLg6Ojo4BERDUlc4EWIpNjCLVJ0+z5urHZI6KHF+zvhohBflgWdRXvbotFYEtH+Ljaio6lKA7jJKIaCkorkFtcDoBP9ohIWXyyR0QN7X+G+KOXjzMKSisQscH05u+x2SOiGjRP9RytLWCn5sN/IlJGRWUVEjMLAbDZI6KGY64ywxcTg9DMxgIXkvPw8c8XRUdSFJs9IqrhJhdnISIBkrKKUF4pw9pChRaOfP8hoobj6WiNxaFBAIC1R5Lwc2yK4ETKYbNHRDXcma/HbReISDna+XrutjAzkwSnISJjM7itO14Z2AYA8M6Wc7iWWSQ4kTLY7BFRDdx2gYhE0M7X4+IsRNRI/v14W3T3bob80gq8vvE0yiqqREdqdGz2iKiGm9lciZOIlBfHxVmIqJFZqMzwvxO7wtHaAudu5GLBL3+LjtTo2OwRUQ0pudV77Hk4chgnESnnKjdUJyIFtHSyxqLxXQAA3x5OwK8XUgUnalxs9oiohrT8UgCAhwObPSJShizLfLJHRIoZ0qE5Xur/CADg3z/E4Ea28c7fY7NHRFqyLCP19pO95mz2iEghKbklKCyrhLmZBG8X09rwmIjEeGdoO3TxckJeSQX+tfEMyiuNc/4emz0i0sovrUDx7c1G2ewRkVI0i7N4u9jAQsU/TYio8Vmam2HZxK6wtzLHmWs5WLj3kuhIjYLvqESklZZX/VTPwcoc1pYqwWmIyFTE3Z6v58uVOIlIQV7ONvhsXPX8vf9/IB5//n1LcKKGx2aPiLRSc6vn6/GpHhEp6Srn6xGRIEM7eSC8nw8AYMbmGKTkFosN1MDY7BGR1q08rsRJRMpjs0dEIr07rB0CWzoiu6gc0zaeQYURzd9js0dEWrfyq5s9d3s2e0SknDhuu0BEAqnNVVgW1hV2anOcSMzG579fFh2pwbDZIyKtW9qVONWCkxCRqcgpKkNGQRkAztkjInG8XWyxYGwgACAyKg77L6cLTtQw2OwRkdatvNt77HEYJxEpRDOEs4WjFWzV5oLTEJEpG9G5BZ7p3RoAMH3TWe30lqaMzR4RaaXmcRgnESlL0+z5cggnERmA2SM6oL2nAzILy/DG92dQWSWLjvRQ2OwRkVYaF2ghIoVxcRYiMiRWFipEhnWFraUKR+Oz8MUfV0RHeihs9ogIAFBVJSMtX7P1AufsEZEyrnJxFiIyMG3c7PDxmOr5e0v/vILDVzMEJ6o/NntEBADILCxDRZUMSQJc7djsEZEytE/2uDgLERmQUUEtMaGnF2QZeOP7s0jLb5rz99jsERGAO3vsudqpYaHiWwMRNb7iskrczKnewJhP9ojI0Mx5siPaNrdHRkEp3tx0tknO3+NfdEQEANpPrDiEk4iUEpdeAFkGmtlYwIUjCojIwFhbqhD5TFdYW6hw+GomIqOuio70wNjsEREAIDX39nw9rsRJRArhZupEZOj83O0xf3QnAMCS3y/jaHym4EQPhs0eEQG4M4yzOVfiJCKFcCVOImoKxnVvhbHdWqFKBqZtPIOMglLRke4bmz0iAnDXME4+2SMihWj32OPiLERk4OaP7gg/dzuk5Zdi+uYYVDWR+Xts9ogIAJCayzl7RKQsPtkjoqbCxtIcy8K6Qm1uhgOX0/H1gTjRke4Lmz0iAgDcyrs9Z4/DOIlIARWVVUjMLATAZo+ImoZ2Hg74cGRHAMCiXy/jRGKW4ET3xmaPiABwGCcRKSspqwjllTKsLVRo4WgtOg4R0X0J7emFUUEtUFklY9rGM8guLBMdqU5s9ogIZRVVyCiofrPiME4iUkKcZr6euy3MzCTBaYiI7o8kSfh/TwWijastUnJLMOOHGMiy4c7fY7NHREi/vaqUhUqCs62l4DREZAqupnNxFiJqmuzU5lgW1g2W5mb48+80rDiYIDqSXmz2iEi7OIu7vRUkiZ+wE1Hj0y7OwmaPiJqgDi0c8MGIDgCA/+z5G6evZQtOpBubPSJC2u099jy4OAsRKSSOK3ESURP3TO/WGN7ZExVVMv614Qxyi8pFR6qFzR4R3dlQnfP1iEgBsiwjLp0rcRJR0yZJEhaMCYS3iw1u5hTjrS2GN3+PzR4RIfX2tgvuXImTiBSQmleCgtIKqMwkeLvYio5DRFRv9lYWiAzrBkuVGX796xZWRyeKjlQDmz0i4jBOIlKUZr6et4sNLM35pwgRNW2dWjriveHtAQAf/3wR527kCE50B99hiQi38jmMk4iUw8VZiMjYPN/XG0M7eqC8UsbrG84gr8Qw5u+x2btPn3/+OTp27Ag7Ozs4OTnhsccew7Fjx0THImoQmtU4uaE6ESnhKhdnISIjI0kS/jOuM1o1s8a1rCLM3HrOIObvsdm7T97e3li8eDFiYmIQHR0NPz8/hISEIDMzU3Q0ooeWdnvOXnMO4yQiBbDZIyJj5GhtgWVh3WChkvBzbCq+O5okOhKbvfs1ZswYhISEwNfXFx06dMDChQuRm5uL8+fPi45G9FAKSyuQX1oBAGjuwGaPiBofV+IkImMV5OWEd4a2AwDM330R52/mCs1j0M3ed999h1deeQU9evSAWq2GJElYvXp1ndecOHECw4YNg5OTE2xtbdGnTx9s3ry5QXOVlZXhm2++QbNmzRAYGNigX5tIaZptF2wtVbBTmwtOQ0TGLr+kHBkF1aMJHnHlSpxEZHxe6v8IhrR3R1llFV7fcBr5AufvGXSz9/777+Obb75BUlISPD097/n6qKgoPProozh06BCefvppTJ06FampqQgNDcWiRYseOs/BgwdhZ2cHa2trfP755/jtt9/g7Oz80F+XSKRbHMJJRApKzCgCALjaWcLeykJwGiKihidJEhaO74IWjlZIzCzCrO3nhc3fM+hmb8WKFUhMTER6ejqmTp1a52srKiowZcoUmJmZ4cCBA/jmm2+waNEixMTEICAgALNmzUJSUs1xszNnzoQkSXX+u1uPHj1w9uxZREdH44knnsDTTz+NjIyMBv+5iZSUls/FWYhIOQmZ1UM4+VSPiIyZk40lloZ1hcpMwo8xyfj+xHUhOQy62RsyZAi8vb3v67V//vkn4uLiEBYWhqCgIO1xR0dHzJo1C2VlZVizZk2Na2bMmIGLFy/W+e9u1tbW8PPzQ+/evbFixQqYmZlh1apVD/+DEgmkXYmT2y4QkQISM6qbPR9upk5ERq67tzPeCmkLAJi76wIupuQpnsFoJujs27cPAPD444/XOhcSEgIA2L9/f43jbm5ucHNzq/f3lGUZpaWl9b6eyBBwGCcRKSlB0+zxyR4RmYCXg9vgaHwm9l1KR8SG0/jx9f6wVXCNBKNp9q5cuQIA8Pf3r3XOw8MDdnZ22tfUxzvvvIORI0eiVatWyMrKwpdffokbN25g7NixOl9fWlr6wI1gbm71aj15ecp3/WS6rt/KQFVpERzMyvn/Hpk0zf//hrAvkiHR3I+Gen+4fD0NVaVFaG5VxfccIjIJc0MewbiEVFy9kY63NxzFx2MCa00X06Uh6pLRNHuaRsnR0VHneQcHB+1r6iM5ORkTJkxAWloanJ2d0bNnTxw8eBDt27fX+fpPPvkEH374Yb2+l5eXV71zEtXX60uA10WHIDIAmZmZemuJKdLsJ9vQtWn8kgb9ckRETcKXAL588cGueZi6ZDTNXmNbt27dA73+3XffxfTp0x/ompycHHh7e+PatWv8Q+O/5OXlwcvLC9evX4eDg4PoOAaD90U/3hvdeF/0y83NRevWrbnK8n/R3A/Wppr4u6Qf741uvC/68d7o1hB1yWiaPU0B0vf0Li8vD82aNVMsj1qthlpdvwUvHB0d+T+6Hg4ODrw3OvC+6Md7oxvvi35mZga9dpniNPeDtUk3/i7px3ujG++Lfrw3uj1MXTKaiqaZq6drXl5qaioKCgp0zucjIiIiIiIyRkbT7A0cOBAA8Ouvv9Y6t3fv3hqvISIiIiIiMnaquXPnzhUd4n4cOnQIf/zxB0aPHl1jHz0NHx8frF+/HtHR0RgxYgQ8PDwAVA/rfP7555Gfn48VK1bAyclJ6egPRKVSYdCgQTA3N5oRtg2G90Y33hf9eG90433Rj/dGN94X3Xhf9OO90Y33RT/eG90e9r5IsgGvMb1ixQocOnQIABAbG4vTp0/j0UcfhZ+fHwCgf//+mDx5svb1UVFRCAkJgZWVFSZMmAB7e3ts3boVSUlJWLhwIWbMmCHk5yAiIiIiIlKaQTd74eHhWLNmjd7zL7zwAlavXl3j2PHjxzFnzhxER0ejvLwcgYGBmD59OkJDQxs5LRERERERkeEw6GaPiIiIiIiI6sdoFmghIiIiIiKiO9jsGYDi4mJ88MEHCAgIgJWVFVq0aIEXX3wRN2/eFB1NmFOnTmHBggUYM2YMWrVqBUmSIEmS6FjCFRUVYceOHXjppZfQtm1bWFlZwdbWFl26dMG8efNQUFAgOqJQixcvxpgxY+Dv7w9HR0eo1Wp4e3vj+eefR2xsrOh4BiEzMxPu7u6QJEk7/9mUDRo0SPv+ouvfnj17REcUhrWpNtYm3Vib9GNduj+sTXc0dF3iME7BSkpKMHjwYBw9ehSenp4IDg5GYmIijh8/Djc3Nxw9ehRt2rQRHVNxo0ePxs6dO2sdN/X/XVesWIEpU6YAANq3b49OnTohLy8P0dHRyM/PR7t27bB//364u7sLTiqGq6srCgsL0blzZ7Rs2RIAcOHCBVy+fBkWFhbYtm0bRowYITilWOHh4Vi7di1kWYavry+uXr0qOpJQgwYNwv79+zF27FjY2dnVOj9jxgwEBgYKSCYWa5NurE26sTbpx7p0f1ib7mjwuiSTUO+9954MQO7bt6+cn5+vPb5o0SIZgDxw4EBx4QRasGCBPHv2bHnXrl1ySkqKrFarZf7vKsurV6+WX375Zfmvv/6qcTw5OVnu2rWrDECeOHGioHTiHTp0SC4uLq51PDIyUgYgN2/eXC4vLxeQzDD8/vvvMgD55ZdflgHIvr6+oiMJN3DgQBmAnJCQIDqKQWFt0o21STfWJv1Yl+6Ntammhq5LfLInUFlZGdzd3ZGbm4vTp0+ja9euNc536dIF586dw8mTJ9G9e3dBKQ2DlZUVSktLTf7T07ocOXIE/fr1g1qtRl5eHiwtLUVHMih+fn6Ii4tDTEwMOnfuLDqO4oqLixEYGAi1Wo0dO3YgICDA5D89Be58gpqQkAAfHx/RcQwCa9P9Y226N9Ym/Uy9LgGsTbo0dF3inD2BDh8+jNzcXPj6+tYqpgAwbtw4AMCPP/6odDRqgrp06QIAKC0tRWZmpuA0hsfCwgIATPYPjQ8//BDx8fH4+uuvtfeCSBfWJmpIrE36mXpdAliblMAt6gWKiYkBAHTr1k3nec3xc+fOKZaJmq74+HgA1cXD2dlZcBrDsm7dOly6dAn+/v7w9/cXHUdx586dw6JFizBp0iTt3CuqaeXKlcjMzISZmRkCAgIwevRotG7dWnQsIVibqCGxNulm6nUJYG26l4aqS2z2BLp27RoAoFWrVjrPa44nJSUplomari+++AIAMHToUKjVasFpxPrss89w4cIFFBYW4uLFi7hw4QJatGiBjRs3QqVSiY6nqKqqKkyePBlOTk749NNPRccxWB999FGN//73v/+N2bNnY/bs2YISicPaRA2Jtaka61JNrE331lB1ic2eQJqliG1sbHSet7W1BQDk5+crlomapp9//hkrV66EhYUF5s+fLzqOcHv37sUff/yh/W9vb2+sXbvWJOcXLV26FCdOnMCqVavg4uIiOo7BGTBgACZPnox+/frB09MT169fsgAm9AAADlNJREFUx5YtW/DRRx/hgw8+gIODA9544w3RMRXF2kQNhbXpDtalmlib9GvwutQgy7xQvUyZMkUGIL/33ns6z1+5ckUGIPv7+yuczPBwxTP9Ll68KDdr1kwGIC9ZskR0HIOSnZ0tHzhwQB4yZIgMQP7oo49ER1JUUlKSbGdnV2vlxISEBK54dg979+6VAchOTk5yUVGR6DiKYm26f6xN+rE26WbqdUmWWZvqq751iQu0CKTZO6OoqEjn+cLCQgCAvb29Ypmoabl58yaGDh2K7OxsTJ8+3eSeQNyLk5MTgoOD8fPPP6N79+6YPXs2Tpw4ITqWYiIiIlBWVoavv/5adJQm5/HHH0ePHj2Qk5ODY8eOiY6jKNYmelisTfqZel0CWJvqq751icM4BdJMsrxx44bO85rj3t7eimWipiMrKwuPP/44kpKSMGnSJCxcuFB0JINlYWGB0NBQnDp1Cj/++CN69uwpOpIidu/eDScnJ0ydOrXG8ZKSEgDVf5ANGjQIAPD999/Dw8ND6YgGzd/fHydPnkRKSoroKIpibaKHwdp0f0y1LgGsTQ+jPnWJzZ5AmuWIT58+rfO85rip7r1C+hUUFOCJJ57AX3/9hTFjxmD58uWQJEl0LIPm6uoKAEhPTxecRFk5OTnYv3+/znMlJSXac5oiS3dkZ2cDuDNHzVSwNlF9sTY9GFOtSwBrU33Vpy5xGKdAjz76KBwdHREXF4ezZ8/WOr9lyxYAwJNPPql0NDJgpaWlGDVqFI4fP46QkBCTXcnrQWkKh6+vr+AkypFlWee/hIQEANX3QnOMG4rXlJ6ejoMHDwLQvwWBsWJtovpgbXpwpliXANam+qpvXWKzJ5ClpSVef/11ANXjlzXzIABg8eLFOHfuHAYOHGiyKzVRbZWVlZg4cSL+/PNPBAcHY9u2bSa9GevdDh8+jD179qCqqqrG8fLycixduhTr1q2DtbU1QkNDBSUkQxMdHY0dO3agsrKyxvHExEQ89dRTKCwsxMiRI/VuQWCsWJvoQbE26ca6RA+qMeoSh3EK9v777+P3339HdHQ0/P39ERwcjKSkJBw7dgxubm749ttvRUcU4qeffqqxTHNZWRkAoE+fPtpjs2fPxvDhwxXPJtKyZcuwfft2ANXDP1577TWdr1u4cKF2eIipuHLlCiZNmgRXV1d0794dLi4uyMjIQGxsLFJSUmBlZYXVq1fDy8tLdFQyEJcvX8akSZPg4eGBbt26wcnJCUlJSTh16hRKSkrQsWNHLF++XHRMIVibdGNt0o21STfWJXpQjVGX2OwJZmVlhaioKHzyySfYsGEDduzYAWdnZ4SHh2P+/Pkm94myRnp6us6Vhu4+Zopj3DVjtQFoC6suc+fONamCCgADBw7ErFmzsH//fpw7dw4ZGRmwtLSEj48Pxo0bh2nTpsHPz090TDIgvXv3xquvvopjx47hxIkTyM7Ohq2tLYKCgjB+/Hi8+uqrsLa2Fh1TCNYm3VibdGNt0o11iR5UY9QlSZZluZHyEhERERERkSCcs0dERERERGSE2OwREREREREZITZ7RERERERERojNHhERERERkRFis0dERERERGSE2OwREREREREZITZ7RERERERERojNHhERERERkRFis0dERERERGSE2OwRCSJJkvbfkSNH9L5u8+bN2tf5+PgoF7ABGErm1atXQ5IkzJ07V3QUIiIiIsWw2SMyAOvXr9d77rvvvlMwSePbt28fJElCeHi46ChERERERo3NHpFAKpUKgYGB2LRpEyoqKmqdz8zMxJ49e9CtWzcB6R7exYsX8ccff4iOQURERGSS2OwRCfbMM88gIyMDe/furXVu06ZNKC8vx7PPPisg2cNr164dfH19RccgIiIiMkls9ogECwsLgyRJOodrfvfdd7Czs8OoUaN0XivLMjZu3IgJEyYgICAAtra2sLe3R69evfDll1+iqqpK53WFhYWYOXMmfHx8YGVlBT8/P8yfPx/l5eXw8fGBJEk1Xn/30MusrCy8+uqr8PT0hFqtRqdOnfDtt9/q/D7/PWcvPDwcgwcPBgCsWbOmxrxFzXy6ew3zDA8PhyRJ2LdvX61zhw8fxpAhQ2Bvbw8nJyeEhITg2LFjOr+ORkVFBb766iv07dsXDg4OsLa2RlBQEJYsWaLzaSsRERFRU2EuOgCRqfPy8sKAAQOwa9cuFBQUwM7ODgAQHx+PI0eO4LnnnoONjY3Oa0tLSxEWFgYXFxd06NAB3bp1Q2ZmJqKjoxEREYHjx49j9erVta4ZMmQIjh49CmdnZ4wYMQKlpaVYsGABTp8+XWfWnJwc9O3bFwUFBQgODkZGRgYOHDiAl156CVVVVZg8eXKd1/fv3x+pqanYu3cvfH190b9/f+25oKCg+7hb+u3evRtPPfUUKioq0KtXL7Rp0wYxMTEYMGCA3saxuLgYw4cPR1RUFJydndGnTx9YWVnh2LFjePPNNxEVFYXt27fDzIyfixEREVETJBOREABklUoly7IsL1++XAYgr1mzRnt+3rx5MgB57969ckpKigxA9vb2rvE1ysvL5e3bt8tlZWU1jqelpck9evSQAcj79++vcW7+/PkyALlXr15ydna29nhCQoLs5eUlA5D/+60hKipKe3zChAlySUmJ9tz27dtlAHLr1q11/oz/nVnztV544QWd9+Ve51944QUZgBwVFaU9lpeXJ7u5uckA5G+//VZ7vKqqSn7nnXe02efMmVPja7322msyADk0NFTOycmp8fWGDRsmA5C/+uornTmIiIiIDB0/riYyAOPGjYNara6xKuf69evh6emJf/zjH3qvMzc3x+jRo2FhYVHjuJubGz755BMAwM6dO2uc+/rrrwEAixYtgpOTk/a4j48PPvjggzpzOjg4YNmyZVCr1dpjo0ePRqdOnXDt2jUkJibW/YM2ki1btiA9PR0DBgzApEmTtMclScL8+fPRqlWrWtekpaVh+fLl8PLywqpVq+Do6Kg9Z29vj5UrV8LS0hJfffWVIj8DERERUUPjME4iA+Dk5IThw4dj586dSE1NxfXr13Hp0iW8+eabUKlU97z+7Nmz+PXXX5GUlISioiLIsoz8/HwAwJUrV7SvS0pKws2bN+Hh4VFjCKVGaGgopkyZovf7dO/eHS4uLrWOBwQE4Pz580hJSRGyr97BgwcBABMmTKh1zsLCAuPGjcOSJUtqHN+3bx/Ky8sxdOhQWFtb17rOw8MD/v7+iI2NRXFxsc7XEBERERkyNntEBuLZZ5/Ftm3b8P333yMhIUF7rC5lZWUIDw/Hxo0b9b5G0/QBQEpKCoDqeYK6aBY2ycnJ0Xle1xMyzXVA9XxAEZKTkwEA3t7eOs/rakA1TyGXL1+O5cuX1/n1s7Ky0LJly4fKSERERKQ0NntEBmLYsGFwcnLC2rVrkZycjPbt299zf73Fixdj48aNCAwMxKeffopu3bqhWbNmsLCwwOXLl9G2bVvIstxgGQ1hoRJ9K4zW9+sEBQWhS5cudb727mGrRERERE0Fmz0iA6FWqzF+/HjtU6Zp06bd85rt27cDADZu3IiOHTvWOBcfH1/r9Z6engCA69ev6/x6+fn5ep/qKcXS0hIAUFBQoPO8ruyanyspKUnnNbqOa55S9u/fH0uXLq1XViIiIiJDJv5jeiLSeu655+Di4gJXV1c888wz93x9dnY2AN3DKzdv3lzrmLe3N1q2bInU1FRER0fXOv/DDz/UI/WD0TRz+vaw0zRuly9frnUuKytL5/YQwcHBAHT/zBUVFdi6dWut44MHD4ZKpcLu3btRXl5+/z8AERERURPBZo/IgGj2rktPT9c7/+xuAQEBAO6ssKmxZcsWrF27Vuc1U6dOBQDMmDEDubm52uNJSUmYN29efaPftxYtWgAALl26pPP8I488gtatWyM2NrbGSqKFhYV4+eWXkZeXV+ua8ePHw8XFBfv27cOaNWu0x2VZxpw5c3Dt2rVa17Rs2RIvvvgiEhMTMXHiRNy6davWa65evaqzUSQiIiJqCtjsETVhb7/9NlQqFWbOnIkePXogLCwMPXv2xPjx4/Hmm2/qvOatt95Cnz59cPToUfj6+mL8+PEYOXIkOnbsiC5duqB169a1tnJoSD4+PujcuTNOnjyJXr16YdKkSZg8eTJ27dqlfc2cOXMAAGPHjsVjjz2GkSNHwtfXF+fOncOoUaNqfU3NVgkqlQrh4eHo06cPwsLC0KlTJ3z22Wd6Vxj94osv8M9//hNbt27VbvIeFhaGUaNGwd/fH/7+/li3bl3j3AgiIiKiRsZmj6gJGzBgAA4dOoTHHnsM8fHx2L17NywtLbF161ZERETovEatVuO3337D22+/DVtbW+zatQsXLlzAjBkzsGnTJty6dUvn9goNaevWrRg9ejTi4+Oxdu1arFy5ssbwzBdffBGrVq1C+/btcfjwYRw/fhxPPvkkjhw5UmNvwLuNGjUKUVFRGDx4MM6fP4+ffvoJnp6e2L9/P/r166fzGmtra/zyyy9Ys2YNevfujYsXL2LLli04efIk3Nzc8OGHH+LTTz9tlHtARERE1NgkuSGX6iOiJu3o0aPo27cvhg4dil9++UV0HCIiIiJ6CHyyR2SCzpw5U2sLg/j4eLzyyisA7r2/HxEREREZPj7ZIzJB7dq1Q25uLgIDA+Hi4oJr167h1KlTKC0txciRI7Fjxw5IkiQ6JhERERE9BDZ7RCYoMjISmzdvxqVLl5CVlQUrKyt06NABzz77LKZOnQpzc27BSURERNTUsdkjIiIiIiIyQpyzR0REREREZITY7BERERERERkhNntERERERERGiM0eERERERGREWKzR0REREREZITY7BERERERERkhNntERERERERGiM0eERERERGREfo/n1gnLVNm6gEAAAAASUVORK5CYII=\n", | |
"text/plain": [ | |
"<Figure size 1000x600 with 2 Axes>" | |
] | |
}, | |
"metadata": {} | |
} | |
] | |
}, | |
{ | |
"cell_type": "code", | |
"source": [ | |
"plt.clf()" | |
], | |
"metadata": { | |
"id": "R4hJGr4SjuaD" | |
}, | |
"execution_count": null, | |
"outputs": [] | |
} | |
] | |
} |
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
1966 09 14 05 42 29 51.6 143.9 15 3.6 | |
1967 01 09 20 26 09 51.5 143.5 10 3.6 | |
1967 05 21 03 52 25 51.4 144 15 3.4 | |
1968 01 25 23 53 50 51.3 143.4 15 3.6 | |
1970 11 04 11 34 53 51.2 143.4 15 3.3 | |
1970 11 14 17 48 45 51.4 143.1 15 3.4 | |
1972 09 09 23 06 11 51.2 143.5 10 3.3 | |
1975 01 04 03 58 48 51.4 143.5 15 3.9 | |
1976 04 16 08 30 33 51.1 143.3 10 3.9 | |
1977 10 23 05 07 35 51.6 143 15 3.4 | |
1988 07 23 19 31 20 51.46 143.81 7 3.5 | |
1988 08 24 05 44 49 51.36 143.48 9 4.5 | |
1989 04 17 03 50 34 51.55 143.01 8 4.2 | |
1992 01 06 14 21 33 51.14 143.36 9 3.1 | |
1994 04 06 05 40 16 51.3 143.8 6 3.1 | |
1995 10 05 19 27 26 51.35 143.45 10 3 | |
1997 01 22 20 51 23 51.32 143.61 10 3.6 | |
1997 07 09 08 00 43 51.56 143.28 9 3 | |
1997 07 10 22 17 02 51.54 143.03 10 3.3 | |
1998 04 30 21 31 39 51.43 143 9 3 | |
2000 05 09 05 03 16 51.59 144 9 3.2 | |
2002 01 03 10 49 58 51.22 143.83 14 3.6 | |
2006 10 03 14 37 19 51.422 143.058 4.1 1.4 | |
2006 10 13 20 04 11 51.569 143.085 10 3.2 | |
2006 11 01 19 36 08 51.444 143.2 13.7 1.2 | |
2006 11 28 07 09 50 51.106 143.014 18.3 1.3 | |
2007 02 15 12 15 23 51.327 143.369 4.1 1.9 | |
2007 02 20 17 14 18 51.28 143.55 2 1.5 | |
2007 03 16 08 41 52 51.494 143.008 0 0.7 | |
2007 04 01 17 57 10 51.531 143.16 0 2.1 | |
2007 04 05 16 41 29 51.458 143.114 24.2 1.4 | |
2007 06 12 10 49 58 51.397 143.039 0 1.2 | |
2007 06 16 13 44 23 51.272 143.374 26.5 1.1 | |
2007 06 25 15 57 01 51.469 143.111 5.2 1.9 | |
2007 08 25 04 51 10 51.287 143.368 23.7 2.5 | |
2007 08 31 12 02 32 51.237 143.339 3 1.4 | |
2007 09 19 17 34 21 51.357 143.084 4.2 1.4 | |
2008 08 08 20 35 25 51.29 143.035 10 0.6 | |
2009 02 24 10 21 12 51.229 143.166 12.1 1.1 | |
2009 04 01 12 29 50 51.175 143.586 2 1.3 | |
2009 06 18 17 20 49 51.396 143.352 3.4 1.6 | |
2009 10 15 09 52 47 51.276 143.022 8.3 2.2 | |
2009 11 18 01 49 05 51.426 143.01 31.7 2.9 | |
2009 11 18 06 57 05 51.364 143.161 13.9 1.8 | |
2009 12 16 08 57 57 51.435 143.072 15.1 1.8 | |
2009 12 16 11 29 45 51.433 143.088 18 1.7 | |
2010 01 11 15 26 34 51.465 143.018 10 0.1 | |
2010 08 23 13 36 11 51.33 143.05 3.5 0.8 | |
2010 10 12 17 34 56 51.212 143.043 10 0.5 | |
2011 03 15 10 53 07 51.58 143.212 14.3 2 | |
2011 03 15 21 12 24 51.599 143.178 14.4 2.8 | |
2011 04 26 22 05 08 51.29 143.021 15 1.3 | |
2011 06 11 07 19 41 51.536 143.009 11.7 4 | |
2011 07 21 17 22 57 51.468 143.124 15 2.2 | |
2011 10 19 04 47 09 51.439 143.054 4 2.1 | |
2011 11 10 12 37 04 51.531 143.189 14 1.8 | |
2011 12 25 03 02 43 51.21 143.227 4.2 2.2 | |
2012 01 27 19 02 37 51.42 143.194 17 1.3 | |
2012 03 18 18 17 33 51.307 143.406 13 1.5 | |
2012 05 28 10 26 26 51.396 143.092 9.3 1.8 | |
2012 09 28 21 19 45 51.244 143.413 4.1 2.1 | |
2012 10 05 03 55 07 51.325 143.272 8.3 1.5 | |
2012 10 07 02 24 47 51.565 143.191 22.6 2.5 | |
2012 10 07 03 03 13 51.593 143.173 24 2.2 | |
2012 10 21 22 11 50 51.334 143.239 4.3 1.3 | |
2012 12 28 06 53 11 51.351 143.091 5 2.5 | |
2013 01 03 10 52 12 51.301 143.362 10 1.6 | |
2013 01 16 13 41 36 51.221 143.017 4.7 1.7 | |
2013 02 09 14 16 56 51.456 143.194 8 2.6 | |
2013 02 11 04 18 22 51.471 143.207 5 2.9 | |
2013 03 27 03 36 33 51.459 143.177 7.5 1.6 | |
2013 04 05 14 32 25 51.506 143.073 15.1 1.6 | |
2013 04 14 03 23 55 51.537 143.18 7.5 2.1 | |
2013 05 18 21 23 00 51.133 143.368 15 2.4 | |
2013 08 06 12 46 28 51.343 143.402 13.4 2.7 | |
2013 08 06 15 44 24 51.309 143.404 4.2 1.1 | |
2013 08 06 16 47 33 51.565 143.107 14.8 2.1 | |
2013 08 09 17 46 33 51.566 143.102 14.6 1.5 | |
2013 08 11 15 00 44 51.566 143.098 12.5 1.3 | |
2013 09 20 03 10 35 51.404 143.459 9.4 3 | |
2013 09 20 12 27 12 51.429 143.415 9.1 2.1 | |
2013 09 20 21 16 10 51.424 143.397 5 1 | |
2013 09 20 21 17 28 51.431 143.38 5 0.5 | |
2013 10 14 02 31 04 51.424 143.441 9.8 3.3 | |
2013 10 14 02 31 24 51.425 143.417 9.1 3.6 | |
2013 10 14 02 32 02 51.424 143.448 9.3 3.2 | |
2013 10 14 07 34 14 51.43 143.426 8.6 2.9 | |
2013 10 14 08 36 24 51.42 143.463 8.9 4.3 | |
2013 10 14 08 50 12 51.413 143.396 9.1 2 | |
2013 10 14 16 09 42 51.431 143.439 8.4 2 | |
2013 10 15 09 39 04 51.313 143.011 14.3 3.1 | |
2013 10 15 16 23 17 51.313 143.032 13.6 1.8 | |
2013 10 15 20 13 29 51.309 143.035 10.4 1.9 | |
2013 10 16 02 12 48 51.411 143.359 9.8 2.2 | |
2013 10 17 12 29 24 51.403 143.349 9.8 1.7 | |
2013 12 19 14 30 38 51.101 143.66 10 2 | |
2014 01 14 14 15 52 51.418 143.435 2 1.5 | |
2014 01 20 13 52 32 51.235 143.269 4.1 1.2 | |
2014 03 29 09 02 20 51.369 143.133 11.4 1.4 | |
2014 06 12 07 50 40 51.564 143.109 18.7 1.4 | |
2014 06 17 10 50 20 51.56 143.075 15.1 1.9 | |
2014 06 17 14 30 32 51.428 143.401 9.3 1.6 | |
2014 06 17 16 41 32 51.56 143.116 17.2 1.6 | |
2014 06 25 13 19 57 51.4 143.122 15.1 1.3 | |
2014 06 27 00 05 13 51.566 143.113 19.8 2.7 | |
2014 06 27 15 38 46 51.567 143.106 18.7 3 | |
2014 06 27 16 42 35 51.564 143.104 17.5 1.8 | |
2014 06 27 20 41 27 51.567 143.103 18.2 1.7 | |
2014 06 28 00 18 05 51.57 143.104 20.7 2 | |
2014 06 29 16 04 43 51.566 143.096 17.3 2.1 | |
2014 06 29 18 04 53 51.568 143.081 15 1.4 | |
2014 06 29 22 46 45 51.569 143.097 20.9 2.1 | |
2014 06 30 14 42 29 51.563 143.115 20.4 3.9 | |
2014 06 30 14 48 12 51.563 143.108 19.3 2.6 | |
2014 06 30 15 11 27 51.572 143.094 19.8 3.6 | |
2014 06 30 16 19 49 51.572 143.095 18.2 2.3 | |
2014 06 30 16 57 21 51.571 143.086 18.5 1.5 | |
2014 06 30 17 40 39 51.577 143.096 17.8 3 | |
2014 06 30 18 45 23 51.575 143.091 15 1.8 | |
2014 06 30 19 28 22 51.541 143.005 3.9 1.4 | |
2014 06 30 20 58 12 51.57 143.101 16.6 4.5 | |
2014 06 30 21 06 12 51.565 143.086 19.7 2.9 | |
2014 06 30 21 09 57 51.569 143.088 19.8 2.7 | |
2014 06 30 21 18 52 51.504 143.095 21.9 1.8 | |
2014 06 30 21 24 19 51.566 143.078 20 2.2 | |
2014 06 30 21 25 08 51.577 143.08 18.8 2.7 | |
2014 06 30 21 32 10 51.575 143.076 18.5 2.2 | |
2014 06 30 21 37 43 51.575 143.083 19.9 2.6 | |
2014 06 30 22 12 39 51.572 143.075 14.4 1.7 | |
2014 06 30 22 23 45 51.571 143.074 20 2.3 | |
2014 06 30 22 57 15 51.567 143.102 18.9 2.7 | |
2014 06 30 22 57 15 51.571 143.1 20.5 2.7 | |
2014 06 30 22 57 16 51.568 143.102 18.6 2.6 | |
2014 06 30 23 25 13 51.562 143.107 19.1 2.7 | |
2014 06 30 23 57 02 51.569 143.095 15.1 1.8 | |
2014 07 01 00 00 24 51.565 143.102 18.9 1.7 | |
2014 07 01 00 18 31 51.581 143.072 14.8 2 | |
2014 07 01 00 48 38 51.561 143.1 18.8 3 | |
2014 07 01 01 03 20 51.574 143.093 15.1 2.4 | |
2014 07 01 01 14 11 51.579 143.081 13.5 3.2 | |
2014 07 01 05 57 10 51.581 143.086 19.6 2.2 | |
2014 07 01 05 57 10 51.579 143.081 18.7 2.2 | |
2014 07 01 07 09 55 51.582 143.088 15 2.4 | |
2014 07 01 07 42 07 51.574 143.076 10.1 1.5 | |
2014 07 01 07 55 43 51.582 143.072 12.8 1.7 | |
2014 07 01 09 20 37 51.579 143.074 18.3 2.1 | |
2014 07 01 09 57 48 51.565 143.109 16.1 1.6 | |
2014 07 01 10 21 27 51.563 143.083 18.8 1.9 | |
2014 07 01 11 21 37 51.571 143.106 18.2 2.8 | |
2014 07 01 12 30 49 51.565 143.084 18.5 2 | |
2014 07 01 13 47 43 51.549 143.003 4.1 1.2 | |
2014 07 01 14 32 18 51.577 143.092 17.8 2 | |
2014 07 01 17 53 51 51.569 143.078 13.4 1.8 | |
2014 07 01 18 03 06 51.568 143.079 17.3 1.8 | |
2014 07 01 23 29 01 51.576 143.077 18 2.1 | |
2014 07 02 00 16 15 51.566 143.105 17.3 2.4 | |
2014 07 02 03 40 34 51.573 143.068 15 1.9 | |
2014 07 02 10 53 33 51.569 143.074 19.1 2.4 | |
2014 07 02 11 39 06 51.567 143.072 18.5 2.9 | |
2014 07 02 16 13 19 51.561 143.096 18.4 1.4 | |
2014 07 02 17 38 05 51.563 143.096 19.2 1.7 | |
2014 07 02 19 02 47 51.562 143.1 19.3 2 | |
2014 07 02 22 56 44 51.58 143.078 15 2.2 | |
2014 07 03 00 26 54 51.586 143.093 19.2 4 | |
2014 07 03 13 38 09 51.579 143.073 18.7 2.5 | |
2014 07 03 14 28 43 51.56 143.108 18.6 1.8 | |
2014 07 03 20 01 09 51.569 143.096 18.6 2.2 | |
2014 07 03 20 29 56 51.564 143.109 18.2 2.7 | |
2014 07 03 21 01 41 51.571 143.077 12.2 1.6 | |
2014 07 03 22 42 17 51.583 143.104 13.5 1.9 | |
2014 07 04 04 20 00 51.583 143.096 13.3 2.1 | |
2014 07 04 04 27 34 51.585 143.1 14.4 1.8 | |
2014 07 04 18 44 52 51.592 143.062 15 1.9 | |
2014 07 04 21 10 58 51.583 143.1 14.1 1.9 | |
2014 07 05 10 31 30 51.579 143.079 12.9 1.4 | |
2014 07 05 14 27 36 51.589 143.077 15 2.7 | |
2014 07 05 19 13 37 51.567 143.093 15.1 1.6 | |
2014 07 06 08 45 08 51.6 143.085 17.1 1.8 | |
2014 07 06 14 40 13 51.6 143.091 17.1 2.4 | |
2014 07 07 01 20 09 51.561 143.075 15.1 1.9 | |
2014 07 07 02 04 37 51.594 143.077 13.3 2.4 | |
2014 07 07 13 48 04 51.59 143.079 15 1.8 | |
2014 07 09 17 55 56 51.593 143.091 20.5 2.6 | |
2014 07 09 20 13 12 51.594 143.065 15 2.3 | |
2014 07 10 04 12 29 51.598 143.077 14.1 2.3 | |
2014 07 10 12 26 27 51.599 143.089 18.2 1.9 | |
2014 07 11 01 50 18 51.594 143.078 19.2 2.9 | |
2014 07 11 03 38 51 51.592 143.068 11.3 1.9 | |
2014 07 13 17 53 22 51.57 143.079 12.1 1.1 | |
2014 07 14 08 44 55 51.589 143.085 18.2 2.1 | |
2014 07 14 21 40 32 51.588 143.106 16.9 2.3 | |
2014 07 14 21 42 05 51.588 143.104 15.9 2.1 | |
2014 07 15 04 27 50 51.538 143.122 21.9 1.8 | |
2014 07 15 10 04 35 51.587 143.1 13 1.8 | |
2014 07 15 13 01 25 51.569 143.101 15 2.1 | |
2014 07 15 19 58 33 51.561 143.103 16.8 1.4 | |
2014 07 15 21 37 14 51.598 143.088 15 2.6 | |
2014 07 15 22 40 57 51.596 143.094 16.3 1.9 | |
2014 07 16 01 06 39 51.6 143.076 18.2 2 | |
2014 07 16 09 50 17 51.595 143.087 14.1 2.1 | |
2014 07 16 16 04 33 51.575 143.099 12.6 1.3 | |
2014 07 17 08 19 37 51.595 143.087 17.9 2.1 | |
2014 07 17 23 58 27 51.579 143.084 16.8 2.3 | |
2014 07 18 01 14 00 51.577 143.069 12.7 1.5 | |
2014 07 18 23 41 41 51.599 143.087 18.8 3.1 | |
2014 07 19 15 08 56 51.597 143.092 15.1 2.8 | |
2014 07 19 15 30 33 51.6 143.094 17.4 2.8 | |
2014 07 19 15 30 33 51.595 143.096 15.4 2.7 | |
2014 07 19 15 36 01 51.597 143.102 16.7 2.2 | |
2014 07 19 15 45 05 51.596 143.089 15 1.8 | |
2014 07 19 17 28 03 51.599 143.093 16.5 1.5 | |
2014 07 19 17 29 54 51.597 143.097 17.3 2.3 | |
2014 07 19 17 38 41 51.596 143.089 13.5 1.5 | |
2014 07 22 05 31 19 51.564 143.093 12.3 2.2 | |
2014 07 25 10 18 45 51.568 143.062 12.3 1.7 | |
2014 08 02 22 45 36 51.594 143.082 18 1.8 | |
2014 08 03 08 10 10 51.58 143.072 15.7 1.8 | |
2014 08 03 10 38 38 51.585 143.077 17.7 2 | |
2014 08 03 10 57 47 51.582 143.061 22.1 1.4 | |
2014 08 03 16 04 11 51.558 143.11 17.3 1.4 | |
2014 08 06 18 40 31 51.567 143.1 13.6 1.3 | |
2014 08 11 18 03 35 51.584 143.081 10.5 1.1 | |
2014 08 17 17 31 08 51.569 143.09 17.2 1.6 | |
2014 08 21 13 17 52 51.299 143.404 15 1.8 | |
2014 09 04 12 56 11 51.546 143.189 15 1.7 | |
2014 09 14 15 22 38 51.551 143.17 15 2 | |
2014 09 14 15 25 00 51.543 143.199 15 1.6 | |
2014 09 16 15 13 21 51.572 143.192 15.1 1.5 | |
2014 09 19 08 16 42 51.351 143.112 10.6 1.7 | |
2014 10 08 00 54 35 51.587 143.102 13.2 1.8 | |
2014 10 08 02 39 18 51.561 143.091 12.1 2 | |
2014 10 30 16 10 36 51.408 143.5 5.4 2 | |
2014 10 31 12 13 28 51.582 143.11 12.8 1.6 | |
2014 10 31 23 26 03 51.581 143.113 12.9 2 | |
2014 11 01 16 59 11 51.573 143.107 13.3 2.4 | |
2014 11 02 12 50 41 51.551 143.04 4.1 1.3 | |
2014 11 06 20 05 57 51.57 143.124 13.1 1.7 | |
2014 11 12 18 21 16 51.257 143.49 4.1 2 | |
2014 12 10 20 54 22 51.542 143.089 2.4 1.3 | |
2014 12 21 10 59 41 51.481 143.601 12.3 2.3 | |
2014 12 22 17 13 40 51.494 143.622 7.6 2.3 | |
2014 12 29 19 06 07 51.488 143.607 12.2 1.7 | |
2015 01 01 14 23 57 51.372 143.137 4 1.1 | |
2015 01 09 09 43 24 51.567 143.091 12.6 1.6 | |
2015 01 14 08 17 40 51.519 143.599 7.5 1.8 | |
2015 01 18 05 55 30 51.578 143.095 11 1.3 | |
2015 01 18 06 05 47 51.578 143.094 11 1.7 | |
2015 01 18 08 54 22 51.58 143.097 11 1.4 | |
2015 01 18 08 54 22 51.577 143.094 11.1 1.4 | |
2015 01 18 15 52 37 51.578 143.096 14.4 2.2 | |
2015 01 19 12 43 24 51.417 143.118 12.3 1.2 | |
2015 01 21 13 59 41 51.181 143.103 4.2 1.2 | |
2015 01 29 13 57 52 51.575 143.086 12.8 1.3 | |
2015 01 31 07 04 21 51.6 143.095 15.8 1.2 | |
2015 02 21 04 44 01 51.58 143.083 12.4 1.4 | |
2015 03 07 03 07 24 51.427 143.117 13.3 1.3 | |
2015 03 07 20 07 57 51.41 143.334 14.2 0.8 | |
2015 03 14 23 49 11 51.564 143.131 4.4 1.8 | |
2015 04 04 10 44 30 51.56 143.112 13.2 2.1 | |
2015 04 04 20 00 48 51.562 143.113 12 1.1 | |
2015 04 04 20 14 16 51.561 143.111 11.1 1.1 | |
2015 04 05 18 46 20 51.558 143.11 11.9 1.8 | |
2015 04 08 15 11 38 51.563 143.085 11.4 1 | |
2015 05 14 16 44 23 51.445 143.066 11.3 1.2 | |
2015 06 29 17 34 29 51.343 143.487 10.4 2 | |
2015 06 29 17 55 32 51.356 143.506 9.7 1.7 | |
2015 07 15 13 57 07 51.166 143.004 5.9 1.5 | |
2015 08 07 13 57 02 51.591 143.169 4.1 1.2 | |
2015 08 08 18 45 18 51.552 143.181 15.1 1.8 | |
2015 08 08 18 50 02 51.545 143.195 15.1 2.6 | |
2015 08 08 21 53 21 51.551 143.176 15 1.9 | |
2015 08 14 12 02 30 51.556 143.204 15.1 1.8 | |
2015 08 14 19 39 41 51.597 143.151 7.3 1.3 | |
2015 08 21 02 34 35 51.578 143.204 14.4 2.4 | |
2015 08 30 04 03 07 51.54 143.216 15 2.1 | |
2015 08 30 04 38 21 51.558 143.175 4.2 1.6 | |
2015 08 30 16 54 51 51.212 143.091 10 1.6 | |
2015 09 09 11 56 53 51.539 143.269 12 1.3 | |
2015 09 23 09 06 58 51.429 143.108 15 1.9 | |
2015 09 26 00 54 31 51.395 143.13 4.3 1.3 | |
2015 10 08 22 46 56 51.588 143.137 12.2 1.9 | |
2015 10 09 01 32 04 51.583 143.148 10.6 2.1 | |
2015 11 24 18 22 25 51.506 143.473 4.1 1.5 | |
2015 11 27 16 30 30 51.462 143.081 16.2 1.7 | |
2015 12 18 07 49 06 51.392 143.345 27.6 1.6 | |
2015 12 31 04 17 08 51.524 143.171 10.5 2.8 | |
2015 12 31 19 42 27 51.265 143.022 15.1 1 | |
2016 01 02 21 43 14 51.53 143.183 4.3 0.7 | |
2016 01 13 09 25 05 51.417 143.155 15.9 1.5 | |
2016 01 13 10 05 29 51.299 143.078 4.1 1.1 | |
2016 01 13 16 22 07 51.171 143.002 5 1.2 | |
2016 01 15 20 03 22 51.337 143.498 12 1.3 | |
2016 02 07 23 59 52 51.353 143.42 11.3 2 | |
2016 02 29 21 31 25 51.488 143.081 4.2 1.1 | |
2016 03 01 09 55 41 51.544 143.076 12.3 1.1 | |
2016 03 01 12 53 18 51.595 143.111 17 1.8 | |
2016 03 02 08 43 51 51.316 143.234 5 1 | |
2016 03 05 09 32 45 51.489 143.017 3.3 2.7 | |
2016 03 05 11 30 51 51.498 143.006 4 0.8 | |
2016 03 05 11 46 22 51.5 143.008 8 1.2 | |
2016 03 18 09 59 06 51.301 143.272 5 0.7 | |
2016 03 19 00 51 05 51.472 143.088 12.5 2.1 | |
2016 03 24 19 06 48 51.55 143.091 9 2.3 | |
2016 03 24 23 58 25 51.537 143.074 0.6 1.3 | |
2016 04 11 19 04 41 51.369 143.042 10 1.6 | |
2016 04 18 18 54 45 51.565 143.086 15.2 1.1 | |
2016 04 20 21 50 09 51.245 143.338 5 1.5 | |
2016 04 22 05 49 26 51.224 143.181 2 1.9 | |
2016 04 25 18 45 17 51.241 143.081 10 1.7 | |
2016 04 29 15 26 03 51.492 143.608 10 1.1 | |
2016 05 23 12 57 16 51.57 143.136 2.8 1.7 | |
2016 07 20 14 57 15 51.314 143.187 10 1.8 | |
2016 07 29 13 36 29 51.309 143.355 5 1.3 | |
2016 07 30 17 35 14 51.4 143.16 15.8 1.3 | |
2016 08 01 17 28 39 51.446 143.437 5 0.4 | |
2016 08 01 19 25 56 51.441 143.415 6.5 3.1 | |
2016 08 01 19 33 40 51.439 143.424 6.8 1.7 | |
2016 08 03 09 44 56 51.459 143.428 5 0.8 | |
2016 08 04 02 06 13 51.292 143.272 5 1.1 | |
2016 08 05 08 08 25 51.411 143.353 4.1 1.4 | |
2016 08 08 12 49 39 51.426 143.349 10 1.3 | |
2016 08 11 17 34 20 51.527 143.053 6.8 1 | |
2016 08 31 17 45 17 51.459 143.032 10 1.1 | |
2016 10 07 20 25 53 51.344 143.055 11.3 1.7 | |
2016 10 13 22 20 14 51.566 143.134 10.2 1.6 | |
2016 10 26 13 00 49 51.377 143.071 13.5 1.3 | |
2016 11 26 10 31 43 51.362 143.011 16 1 | |
2016 11 27 08 59 17 51.523 143.201 4.3 1.3 | |
2016 11 29 00 08 04 51.6 143.097 13.7 1.8 | |
2016 12 11 07 38 02 51.123 143.201 10 1.5 | |
2016 12 11 17 50 20 51.416 143.596 9.6 1.6 | |
2016 12 20 13 27 53 51.294 143.54 13.5 2.4 | |
2016 12 24 00 30 03 51.557 143.121 11.9 1.7 | |
2016 12 26 22 40 38 51.369 143.043 11.6 1.3 | |
2017 01 09 12 09 36 51.372 143.279 14.5 1 | |
2017 02 21 09 24 04 51.573 143.123 2.4 1.1 | |
2017 02 25 16 27 24 51.51 143.457 21.3 1.2 | |
2017 02 27 14 27 39 51.354 143.031 15 1.1 | |
2017 04 30 08 56 14 51.467 143.016 19 2.8 | |
2017 04 30 08 56 30 51.519 143.039 19 2.6 | |
2017 05 09 23 15 32 51.592 143.315 4.2 2.5 | |
2017 06 10 14 38 37 51.401 143.521 5.5 1.4 | |
2017 07 08 15 39 27 51.245 143.526 12 2.7 | |
2017 07 12 10 57 52 51.523 143.055 1.2 1.1 | |
2017 07 29 22 53 47 51.409 143.371 3.9 1.1 | |
2017 08 28 13 57 55 51.311 143.534 10.5 1.7 | |
2017 08 31 18 42 21 51.409 143.695 10.4 2.3 | |
2017 09 30 00 49 59 51.35 143.408 8.4 1.4 | |
2018 01 21 23 48 06 51.524 143.086 11 1.2 | |
2018 01 22 01 17 18 51.513 143.088 1.5 1.3 | |
2018 02 09 19 37 01 51.338 143.046 4.3 1.1 | |
2018 04 06 20 42 10 51.515 143.535 0.5 2 | |
2018 04 09 16 06 16 51.376 143.329 5 0.9 | |
2018 04 19 11 36 27 51.573 143.133 2.7 0.7 | |
2018 04 24 01 16 45 51.557 143.049 4.3 1.3 | |
2018 04 26 02 29 55 51.559 143.114 2.6 1.4 | |
2018 04 29 04 27 18 51.447 143.548 13 1.9 | |
2018 05 24 17 41 44 51.581 143.079 10.6 2.1 | |
2018 05 28 11 50 26 51.57 143.055 4.2 1.3 | |
2018 06 21 15 33 07 51.468 143.475 4.1 1.4 | |
2018 07 05 22 49 44 51.524 143.187 5 1.8 | |
2018 07 11 12 07 16 51.368 143.489 2.1 1.2 | |
2018 07 15 19 40 42 51.475 143.044 4.2 0.9 | |
2018 07 17 19 40 53 51.36 143.111 8 1.4 | |
2018 07 23 12 54 00 51.254 143.649 10.9 1.8 | |
2018 07 31 10 05 59 51.345 143.267 20 1.6 | |
2018 08 19 11 48 06 51.379 143.537 11.3 1.7 | |
2018 08 28 14 53 17 51.168 143.012 5 1.3 | |
2018 09 11 08 45 33 51.347 143.431 0 1.5 | |
2018 09 18 08 20 40 51.336 143.046 8.5 1.9 | |
2018 09 19 07 55 26 51.336 143.053 4.2 1.4 | |
2018 09 20 07 36 24 51.199 143.201 6.6 1.5 | |
2018 09 25 12 31 14 51.33 143.045 10.4 1.1 | |
2018 10 05 22 01 08 51.375 143.531 11.5 1.7 | |
2018 10 20 14 02 30 51.463 143.055 4.2 1.2 | |
2018 10 28 22 48 13 51.513 143.011 4.2 1.6 | |
2018 11 04 14 47 22 51.381 143.543 9.9 1.8 | |
2018 11 07 19 47 08 51.544 143.178 4.3 1.3 | |
2018 12 03 12 44 10 51.565 143.108 16.2 1.2 | |
2018 12 12 04 59 39 51.223 143.197 2 1.6 | |
2018 12 20 12 17 47 51.374 143.521 10.1 1.3 | |
2018 12 25 15 02 52 51.486 143.475 11.7 3.3 | |
2018 12 25 15 41 01 51.487 143.449 11 3.2 | |
2018 12 30 15 28 31 51.423 143.109 11.9 2 | |
2019 02 19 17 54 51 51.38 143.525 11 0.8 | |
2019 02 27 22 54 45 51.311 143.406 9.6 2.3 | |
2019 03 10 14 46 54 51.456 143.433 14.1 1.4 | |
2019 03 12 15 31 58 51.171 143.937 14.8 1.4 | |
2019 03 19 13 31 31 51.317 143.021 4.2 0.7 | |
2019 04 01 17 31 39 51.563 143.091 2.8 1.1 | |
2019 04 09 19 50 10 51.221 143.056 1.2 1.1 | |
2019 05 01 01 25 30 51.369 143.145 11.8 1.5 | |
2019 05 02 04 13 33 51.589 143.101 15 1.3 | |
2019 05 02 04 43 57 51.232 143.483 0.6 2.1 | |
2019 05 08 13 54 40 51.414 143.147 15.7 1.4 | |
2019 05 13 15 58 32 51.518 143.468 11.8 2.3 | |
2019 05 24 15 48 11 51.268 143.553 7.7 1.5 | |
2019 06 06 23 58 59 51.498 143.527 3.8 1.7 | |
2019 06 17 01 32 22 51.398 143.126 12.1 1 | |
2019 07 03 15 45 47 51.454 143.04 5 0.5 | |
2019 08 05 14 08 03 51.401 143.172 4.1 1.5 | |
2019 08 31 16 32 52 51.468 143.442 0 1 | |
2019 09 03 11 41 07 51.57 143.131 17.5 2.1 | |
2019 09 12 14 20 22 51.458 143.544 5 1.1 | |
2019 10 15 10 02 00 51.568 143.12 12.2 1.9 | |
2019 11 24 16 06 43 51.18 143.155 2.8 1.5 | |
2019 12 05 03 57 08 51.553 143.007 27.6 1.4 | |
2019 12 30 18 22 33 51.554 143.104 11.7 2.6 | |
2019 12 30 18 24 27 51.557 143.102 15.8 1.3 | |
2020 01 01 12 54 47 51.561 143.103 13.1 0.8 | |
2020 01 05 09 39 52 51.563 143.097 14.8 0.9 | |
2020 01 05 13 05 54 51.559 143.091 15 1.1 | |
2020 01 05 13 09 14 51.564 143.1 14.6 0.8 | |
2020 01 05 14 40 37 51.557 143.099 11.4 2.2 | |
2020 01 05 16 14 50 51.562 143.097 12.4 2.3 | |
2020 01 05 16 16 20 51.562 143.105 15 1.1 | |
2020 01 05 16 41 29 51.56 143.101 15 2.6 | |
2020 01 05 16 52 01 51.562 143.1 18.4 1.1 | |
2020 01 06 13 45 28 51.566 143.088 13.8 1.1 | |
2020 01 06 16 38 21 51.565 143.079 16 1 | |
2020 01 06 17 36 38 51.562 143.076 16 0.8 | |
2020 01 09 17 43 21 51.562 143.095 10 1.2 | |
2020 01 11 04 25 35 51.505 143.06 16 1.3 | |
2020 01 12 06 09 02 51.192 143.202 4.1 1.4 | |
2020 01 25 14 43 25 51.369 143.104 11.9 0.9 | |
2020 01 30 14 11 26 51.363 143.107 9.9 1.2 | |
2020 02 05 22 31 07 51.572 143.086 2.3 1.7 | |
2020 02 25 18 05 34 51.323 143.553 1.1 1.4 | |
2020 02 25 18 39 16 51.325 143.547 1.2 1.4 | |
2020 03 02 18 13 25 51.317 143.56 5 1.1 | |
2020 03 22 14 45 13 51.595 143.141 15 1.6 | |
2020 03 28 20 59 14 51.393 143.181 12.9 1.9 | |
2020 03 30 07 30 20 51.398 143.183 7.5 2.6 | |
2020 04 27 20 51 13 51.581 143.081 11.4 2.1 | |
2020 05 13 17 58 42 51.56 143.103 17.1 2.7 | |
2020 05 14 19 02 06 51.561 143.096 15.7 1.8 | |
2020 05 28 12 50 16 51.561 143.093 14.2 1.2 |
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
""" | |
File: magnitude_det_rate_estimation_mle.py | |
Author: Andrey Stepnov | |
Email: [email protected] | |
Github: https://github.com/jamm1985 | |
Description: Estimation the detection bound of local magnitude using a widely used statistical model developed by Ogata and Katsura. **Constant** rate is assumed. | |
Ogata, Y., Katsura, K. (1993). https://doi.org/10.1111/j.1365-246x.1993.tb04663.x | |
Input: Earthquake catalog produced by REPORT command from SEISAN distribution | |
https://seis.geus.net/software/seisan/node124.html | |
Output: Estimated mu (50% detection rate), b-value and sigma. Fig_1.png with the magnitude-frequency relations: ‘+’ are observed magnitudes, line is estimated frequency-densities curves for catalog. X-axes is for local magnitude and Y-axes given in logarithmic scale. | |
Dependencies: scipy >= 1.7.3 | |
""" | |
import numpy as np | |
import matplotlib.pyplot as plt | |
import pandas as pd | |
from scipy import stats | |
from scipy import special | |
from scipy.stats import norm | |
from scipy.optimize import minimize | |
import argparse | |
parser = argparse.ArgumentParser(description='Process input') | |
parser.add_argument('report', type=str, help='Input report.out') | |
parser.add_argument( | |
'--min_mag', | |
type=float, | |
default=1.0, | |
help='min_mag for b-value estimation (default: 1.0)' | |
) | |
args = parser.parse_args() | |
input_report = args.report | |
min_mag = args.min_mag | |
def density(M,theta): | |
"""evaluation function from observation and given parameters | |
:M: M_1,...,M_n | |
:theta: parameter vector (mu, betta, and sigma) | |
:returns: probability for every observation | |
""" | |
beta = theta[0] | |
mu = theta[1] | |
sigma = theta[2] | |
constant = (np.exp(-beta*M)*norm.cdf(M,mu,sigma)).sum() | |
return np.exp(-beta*M)*norm.cdf(M,mu,sigma)/constant | |
def log_l(theta): | |
"""log-likelihood function from observation and given parameters | |
:M_observation: M_1,...,M_n (global context!) | |
:theta: parameter vector (mu, beta, and sigma) | |
:returns: log-likelihood value for all observations | |
""" | |
beta = theta[0] | |
mu = theta[1] | |
sigma = theta[2] | |
N = M_observation.size | |
return (N * np.log(beta) - beta * (M_observation + np.log(norm.cdf(M_observation,mu,sigma))).sum() + np.exp(N*beta*mu - N/2 * beta**2 * mu**2)) | |
def der_log_1(theta): | |
"""first derivatives of log_l | |
:M_observation: M_1,...,M_n (global context!) | |
:theta: parameter vector (mu, beta, and sigma) | |
:returns: gradient vector with values of first derivatives | |
""" | |
beta = theta[0] | |
mu = theta[1] | |
sigma = theta[2] | |
N = M_observation.size | |
der = np.zeros_like(theta) | |
exp_term = np.exp(N*beta*mu - N/2 * beta**2 * mu**2) | |
exp_term_beta_prime = N*mu - N*beta*mu**2 | |
exp_term_mu_prime = N*beta - N*beta**2*mu | |
der[0] = N/beta - (M_observation + np.log(norm.cdf(M_observation,mu,sigma))).sum() + exp_term*exp_term_beta_prime | |
der[1] = N*beta*(1/sigma) + exp_term * exp_term_mu_prime | |
der[2] = beta* (((M_observation - mu)/sigma**2).sum()) | |
return der | |
def hess_log_1(theta): | |
"""second derivatives of log_l | |
:M_observation: M_1,...,M_n (global context!) | |
:theta: parameter vector (mu, beta, and sigma) | |
:returns: hessian matrix with values of second derivatives | |
""" | |
beta = theta[0] | |
mu = theta[1] | |
sigma = theta[2] | |
N = M_observation.size | |
hess = np.zeros((3,3)) | |
exp_term = np.exp(N*beta*mu - N/2 * beta**2 * mu**2) | |
exp_term_beta_prime = N*mu - N*beta*mu**2 | |
exp_term_mu_prime = N*beta - N*beta**2*mu | |
hess[0,0] = -N/beta**2 + exp_term_beta_prime**2*exp_term - N*mu**2*exp_term | |
hess[0,1] = N*(1/sigma) + exp_term_beta_prime*exp_term*exp_term_mu_prime + exp_term*(N-2*N*beta*mu) | |
hess[0,2] = ((M_observation - mu)/sigma**2).sum() | |
hess[1,0] = N/sigma**2 + exp_term_mu_prime*exp_term*exp_term_beta_prime + exp_term*(N-2*N*beta*mu) | |
hess[1,1] = exp_term_mu_prime**2*exp_term - N*beta**2*exp_term | |
hess[1,2] = -beta*N*(1/sigma**2) | |
hess[2,0] = ((M_observation - mu)/sigma).sum() | |
hess[2,1] = -N*beta*(1/sigma) | |
hess[2,2] = -beta*((M_observation - mu)/sigma**2).sum() | |
return hess | |
# estimate b-value | |
# this is modified version of | |
# source: https://github.com/qingkaikong/Learning_Obspy/blob/master/04_catalogAnalysis_20140421/calc_rec.py | |
def calc_recurrence(magnitudes, min_mag = None, max_mag = None, interval = 0.05): | |
""" | |
This function reads an earthquake catalogue file and calculates the | |
Gutenberg-Richter recurrence parameters using maximum likelihood (Aki 1965) | |
approach. | |
Funtion arguments: | |
magnitudes: list of magnitudes values | |
min_mag: minimum magnitude for which data will be used - i.e. catalogue | |
completeness | |
max_mag: maximum magnitude used in bounded G-R curve. If not specified, | |
defined as the maximum magnitude in the catlogue + 0.1 magnitude | |
units. | |
interval: Width of magnitude bins for generating cumulative histogram | |
of earthquake recurrence. Default avlue 0.1 magnitude units. | |
""" | |
# If minimum magnitude is not specified, read all magnitudes | |
if min_mag is not None: | |
pass | |
else: | |
min_mag = -1.0 | |
# If minimum magnitude is not specified default value to minimum in catalogue | |
if min_mag == -1.0: | |
min_mag = min(magnitudes) | |
# If maximum magnitude is not specified default value to maximum in catalogue | |
if max_mag is not None: | |
pass | |
else: | |
max_mag = max(magnitudes) + 0.1 | |
num_eq = len(magnitudes) | |
print('Minimum magnitude: {}'.format(min_mag)) | |
print('Total number of earthquakes: {}'.format(num_eq)) | |
#num_years = max(years)-min(years) | |
#annual_num_eq = num_eq/num_years | |
#print 'Annual number of earthquakes greater than Mw', min_mag,':', \ | |
#annual_num_eq | |
print('Maximum catalog magnitude: {}'.format(max(magnitudes))) | |
print('Mmax = {}'.format(max_mag)) | |
max_mag_bin = max(magnitudes) + 0.15 | |
# Magnitude bins | |
bins = np.arange(min_mag, max_mag_bin, interval) | |
# Magnitude bins for plotting - we will re-arrange bins later | |
plot_bins = np.arange(min_mag, max_mag, interval) | |
########################################################################### | |
# Generate distribution | |
########################################################################### | |
# Generate histogram | |
hist = np.histogram(magnitudes,bins=bins) | |
# Reverse array order | |
hist = hist[0][::-1] | |
bins = bins[::-1] | |
# Calculate cumulative sum | |
cum_hist = hist.cumsum() | |
# Ensure bins have the same length has the cumulative histogram. | |
# Remove the upper bound for the highest interval. | |
bins = bins[1:] | |
# Maximum Likelihood Estimator fitting | |
# b value | |
b_mle = np.log10(np.exp(1)) / (np.mean(magnitudes) - min_mag) | |
beta_mle = np.log(10) * b_mle | |
print('Maximum Likelihood: b value {}'.format( b_mle)) | |
return b_mle, cum_hist[::-1], bins[::-1] | |
# get M_1,...,M_n from catalog | |
catalog = pd.read_fwf(input_report) | |
M_observation = catalog['Ml'].dropna().to_numpy() | |
# obtain initial theta0 | |
mu_0 = M_observation.mean() | |
b_value_0, _, _ = calc_recurrence(M_observation.tolist(), min_mag = min_mag, max_mag = None, interval = 0.05) | |
beta_0 = np.log(10)*b_value_0 | |
std_0 = M_observation.std() | |
theta0=[beta_0,mu_0,std_0] | |
bounds=((0.3*np.log(10),0.8*np.log(10)),(3.0,5.0),(0.01,0.8)) | |
print('theta0 beta, mu, sigma is {}'.format(theta0)) | |
print('bounds is {}'.format(bounds)) | |
# estimate numerically using trust region method | |
res_exact = minimize(log_l, theta0, method='trust-exact', jac=der_log_1, hess=hess_log_1, options={'gtol': 1e-8, 'maxiter': 10000, 'disp': True}, bounds=bounds) | |
print('estimated theta is {}'.format(res_exact.x)) | |
# estimate using constrains | |
res_constrains = minimize(log_l, theta0, method='trust-constr', jac=der_log_1, hess=hess_log_1, options={'xtol': 1e-8, 'verbose': 1, 'maxiter': 10000}, bounds=bounds) | |
print('Current function value:{}'.format(log_l(res_constrains.x))) | |
print('estimated theta is {}'.format(res_constrains.x)) | |
# choose final theta | |
if log_l(res_exact.x) < log_l(res_constrains.x): | |
final_theta = res_exact.x | |
print('exact method is better, final theta is {}'.format(final_theta)) | |
else: | |
final_theta = res_constrains.x | |
print('bounded method is better, final theta is {}'.format(final_theta)) | |
# plot results | |
plt.cla() | |
plt.rcParams.update(plt.rcParamsDefault) | |
plt.rcParams['figure.figsize'] = [10, 6] | |
plt.rcParams['ytick.major.size'] = 10 | |
#plt.rcParams['ytick.major.width'] = 4 | |
plt.rcParams['ytick.minor.size'] = 5 | |
#plt.rcParams['xtick.minor.width'] = 2 | |
#plt.title("Prior 2006-09",fontsize = 15) | |
plt.ylabel('Probability', fontsize = 15) | |
plt.xlabel('Magnitude', fontsize = 15) | |
plt.yscale('log') | |
plt.yticks(fontsize=14) | |
plt.xticks(fontsize=15) | |
plt.ylim([0.001, 1]) | |
plt.xlim([M_observation.min()-0.5, M_observation.max()+0.5]) | |
plt.text(0.5, 0.5,r'$b$={}'.format(np.round(final_theta[0]/np.log(10),decimals=1))+'\n'+r'$\mu$={}'.format(np.round(final_theta[1],decimals=1))+'\n'+r'$\sigma$={}'.format(np.round(final_theta[2],decimals=1)), ha='left', va='center', fontsize = 15) | |
# ECDF | |
x = np.sort(M_observation) | |
y = np.arange(len(x))/float(len(x)) | |
# EPMF | |
y[1:] -= y[:-1].copy() | |
u, inv = np.unique(x, return_inverse=True) | |
sums = np.zeros(len(u), dtype=y.dtype) | |
np.add.at(sums, inv, y) | |
plt.scatter(u, sums, marker='+') | |
# estimated density | |
#x3=np.linspace(0,5,len(x)) | |
plt.plot(u,density(u,final_theta)) | |
plt.savefig('fig_1.png',dpi=600) |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment