{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"You can run this notebook in a [live session](https://binder.pangeo.io/v2/gh/pangeo-data/climpred/main?urlpath=lab/tree/docs/source/examples/decadal/verify_dim_implications.ipynb) [](https://binder.pangeo.io/v2/gh/pangeo-data/climpred/main?urlpath=lab/tree/docs/source/examples/decadal/verify_dim_implications.ipynb) or view it [on Github](https://github.com/pangeo-data/climpred/blob/main/docs/source/examples/decadal/verify_dim_implications.ipynb)."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Implications of `verify(dim)`\n",
"\n",
"This demo demonstrates how setting `dim` in {py:meth}`.PerfectModelEnsemble.verify` and {py:meth}`.HindcastEnsemble.verify` alters the research question and prediction skill obtained, especially for spatial fields.\n",
"\n",
"**What's the purpose of `dim`?**\n",
"\n",
"`dim` is designed the same way as in {py:meth}`xarray.Dataset.mean`. A given metric is calculated over dimension `dim`, i.e. the result does not contain `dim` anymore. `dim` can be a string or list of strings.\n",
"\n",
"**3 ways to calculate aggregated prediction skill**\n",
"- Spatially aggregate first, then `verify`\n",
"- `verify` on each grid point, then aggregate spatially\n",
"- `verify` over spatial, `member` and `init`ialization directly"
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {},
"outputs": [
{
"data": {
"application/javascript": [
"\n",
" setTimeout(function() {\n",
" var nbb_cell_id = 1;\n",
" var nbb_unformatted_code = \"# linting\\n%load_ext nb_black\\n%load_ext lab_black\";\n",
" var nbb_formatted_code = \"# linting\\n%load_ext nb_black\\n%load_ext lab_black\";\n",
" var nbb_cells = Jupyter.notebook.get_cells();\n",
" for (var i = 0; i < nbb_cells.length; ++i) {\n",
" if (nbb_cells[i].input_prompt_number == nbb_cell_id) {\n",
" if (nbb_cells[i].get_text() == nbb_unformatted_code) {\n",
" nbb_cells[i].set_text(nbb_formatted_code);\n",
" }\n",
" break;\n",
" }\n",
" }\n",
" }, 500);\n",
" "
],
"text/plain": [
""
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"# linting\n",
"%load_ext nb_black\n",
"%load_ext lab_black"
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {
"ExecuteTime": {
"end_time": "2020-01-06T18:19:48.753233Z",
"start_time": "2020-01-06T18:19:46.236835Z"
}
},
"outputs": [
{
"data": {
"text/plain": [
""
]
},
"execution_count": 2,
"metadata": {},
"output_type": "execute_result"
},
{
"data": {
"application/javascript": [
"\n",
" setTimeout(function() {\n",
" var nbb_cell_id = 2;\n",
" var nbb_unformatted_code = \"%matplotlib inline\\nimport matplotlib.pyplot as plt\\nimport numpy as np\\nimport xarray as xr\\n\\nimport climpred\\n\\nclimpred.set_options(PerfectModel_persistence_from_initialized_lead_0=True)\";\n",
" var nbb_formatted_code = \"%matplotlib inline\\nimport matplotlib.pyplot as plt\\nimport numpy as np\\nimport xarray as xr\\n\\nimport climpred\\n\\nclimpred.set_options(PerfectModel_persistence_from_initialized_lead_0=True)\";\n",
" var nbb_cells = Jupyter.notebook.get_cells();\n",
" for (var i = 0; i < nbb_cells.length; ++i) {\n",
" if (nbb_cells[i].input_prompt_number == nbb_cell_id) {\n",
" if (nbb_cells[i].get_text() == nbb_unformatted_code) {\n",
" nbb_cells[i].set_text(nbb_formatted_code);\n",
" }\n",
" break;\n",
" }\n",
" }\n",
" }, 500);\n",
" "
],
"text/plain": [
""
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"%matplotlib inline\n",
"import matplotlib.pyplot as plt\n",
"import numpy as np\n",
"import xarray as xr\n",
"\n",
"import climpred\n",
"\n",
"climpred.set_options(PerfectModel_persistence_from_initialized_lead_0=True)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### sample spatial initialized ensemble data and nomenclature"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We also have some sample output that contains gridded data on MPIOM grid:\n",
"- `initialized3d`: The initialized ensemble dataset with\n",
" - `member`s (`1`, ... , `4`) denoted by $m$\n",
" - `init`s (initialization years: `3014`, `3061`, `3175`, `3237`) denoted by $i$\n",
" - `lead` years (`1`, ..., `5`) denoted by $l$.\n",
" - `x` and `j` (think as longitude and latitude): spatial dimensions denotes by $s$.\n",
"\n",
"- `control3d`: The control dataset spanning over `time` (`3000`, ..., `3049`) [not used here].\n",
"\n",
"In {py:meth}`.PerfectModelEnsemble.verify`, we calculate a $metric$ on $data(s,i,m,l)$ containing spatial dimensions $s$ and ensemble dimension $i,m,l$ over dimensions $dim$ denoted by $metric(data(s,i,m,l), dim)$."
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {
"ExecuteTime": {
"end_time": "2020-01-06T18:20:43.341389Z",
"start_time": "2020-01-06T18:20:43.339156Z"
}
},
"outputs": [
{
"data": {
"application/javascript": [
"\n",
" setTimeout(function() {\n",
" var nbb_cell_id = 3;\n",
" var nbb_unformatted_code = \"# Sea surface temperature\\ninitialized3d = climpred.tutorial.load_dataset(\\\"MPI-PM-DP-3D\\\")\\ncontrol3d = climpred.tutorial.load_dataset(\\\"MPI-control-3D\\\")\\nv = \\\"tos\\\"\\ninitialized3d[\\\"lead\\\"].attrs = {\\\"units\\\": \\\"years\\\"}\";\n",
" var nbb_formatted_code = \"# Sea surface temperature\\ninitialized3d = climpred.tutorial.load_dataset(\\\"MPI-PM-DP-3D\\\")\\ncontrol3d = climpred.tutorial.load_dataset(\\\"MPI-control-3D\\\")\\nv = \\\"tos\\\"\\ninitialized3d[\\\"lead\\\"].attrs = {\\\"units\\\": \\\"years\\\"}\";\n",
" var nbb_cells = Jupyter.notebook.get_cells();\n",
" for (var i = 0; i < nbb_cells.length; ++i) {\n",
" if (nbb_cells[i].input_prompt_number == nbb_cell_id) {\n",
" if (nbb_cells[i].get_text() == nbb_unformatted_code) {\n",
" nbb_cells[i].set_text(nbb_formatted_code);\n",
" }\n",
" break;\n",
" }\n",
" }\n",
" }, 500);\n",
" "
],
"text/plain": [
""
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"# Sea surface temperature\n",
"initialized3d = climpred.tutorial.load_dataset(\"MPI-PM-DP-3D\")\n",
"control3d = climpred.tutorial.load_dataset(\"MPI-control-3D\")\n",
"v = \"tos\"\n",
"initialized3d[\"lead\"].attrs = {\"units\": \"years\"}"
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {
"ExecuteTime": {
"end_time": "2020-01-06T18:20:43.701617Z",
"start_time": "2020-01-06T18:20:43.683697Z"
}
},
"outputs": [
{
"name": "stderr",
"output_type": "stream",
"text": [
"/Users/aaron.spring/Coding/climpred/climpred/utils.py:195: UserWarning: Assuming annual resolution starting Jan 1st due to numeric inits. Please change ``init`` to a datetime if it is another resolution. We recommend using xr.CFTimeIndex as ``init``, see https://climpred.readthedocs.io/en/stable/setting-up-data.html.\n",
" warnings.warn(\n",
"/Users/aaron.spring/mambaforge/envs/climpred-dev/lib/python3.10/site-packages/xarray/coding/times.py:351: FutureWarning: Index.ravel returning ndarray is deprecated; in a future version this will return a view on self.\n",
" sample = dates.ravel()[0]\n",
"/Users/aaron.spring/Coding/climpred/climpred/utils.py:195: UserWarning: Assuming annual resolution starting Jan 1st due to numeric inits. Please change ``init`` to a datetime if it is another resolution. We recommend using xr.CFTimeIndex as ``init``, see https://climpred.readthedocs.io/en/stable/setting-up-data.html.\n",
" warnings.warn(\n",
"/Users/aaron.spring/mambaforge/envs/climpred-dev/lib/python3.10/site-packages/xarray/coding/times.py:351: FutureWarning: Index.ravel returning ndarray is deprecated; in a future version this will return a view on self.\n",
" sample = dates.ravel()[0]\n"
]
},
{
"data": {
"application/javascript": [
"\n",
" setTimeout(function() {\n",
" var nbb_cell_id = 4;\n",
" var nbb_unformatted_code = \"# Create climpred PerfectModelEnsemble object.\\npm = (\\n climpred.PerfectModelEnsemble(initialized3d).add_control(control3d)\\n # .generate_uninitialized() # needed for uninitialized in reference\\n)\";\n",
" var nbb_formatted_code = \"# Create climpred PerfectModelEnsemble object.\\npm = (\\n climpred.PerfectModelEnsemble(initialized3d).add_control(control3d)\\n # .generate_uninitialized() # needed for uninitialized in reference\\n)\";\n",
" var nbb_cells = Jupyter.notebook.get_cells();\n",
" for (var i = 0; i < nbb_cells.length; ++i) {\n",
" if (nbb_cells[i].input_prompt_number == nbb_cell_id) {\n",
" if (nbb_cells[i].get_text() == nbb_unformatted_code) {\n",
" nbb_cells[i].set_text(nbb_formatted_code);\n",
" }\n",
" break;\n",
" }\n",
" }\n",
" }, 500);\n",
" "
],
"text/plain": [
""
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"# Create climpred PerfectModelEnsemble object.\n",
"pm = (\n",
" climpred.PerfectModelEnsemble(initialized3d).add_control(control3d)\n",
" # .generate_uninitialized() # needed for uninitialized in reference\n",
")"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"#### subselecting North Atlantic"
]
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
""
]
},
"execution_count": 5,
"metadata": {},
"output_type": "execute_result"
},
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAXkAAAEHCAYAAABLKzaMAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAAqKElEQVR4nO3de7BkZXku8Oddq7t37z0XZoYZYAQiGNFIWYJmRJGUEoiIHk7UEzWQEtFYISfxXlRU1MQkdRJJvCScc1JWRiCSIxIRMVJGEQ7qMReDDhcRBEUFcWBkZmCEYWb27N693vNH90yv71m71+qe3dfl86O66K+/dfn6Mt9ee+2n32XuDhERKado3AMQEZHh0SQvIlJimuRFREpMk7yISIlpkhcRKTFN8iIiJVYZx07N7GwAlwKIAVzm7pfkLV+rzPls9bCDbY/pZ1NknT7qcjNalvuLBlvQH2wsv205/cZRVl43oQcS3naqn7eV0MI8Duo/4eSnopv7brs/d1MWhS9wc3U9aDdW0ApR+kXgjXUdRpf+7nFg/hgUbcuWsa3Muplts87y/W47og8VH7VZqj+7bELL5u9rObyvf0jL88j3du109w2D3OYH3nW4r14V4d1/vmN0T2RAbNQ5eTOLAfwAwEsBbAXwbQDnufv3uq1z2OxGf9HxbzrYbh42G/Qvrqx27s/FQV+zGr4nSc3y++nHnqd+gGQ+p/TSRU36R7RI/Q3v2o4Xwr54f/gPMN4Xbiyab4bt/Y3OMBthH/bNh+1GuC2fD/tv2LkZ3bxs9vxw3SaNY8Vc0H7irGcF7UdOodd/VWd9r4SvgUX8Auf3Z9vp++HrGRWsyxNierKsxPnbimlfsRX0p9rZiThsV+Pw9a7HiwXtzudiNnW/1V4I2jP0geWxLEdSeDSVs26fJxw+evI1t7r7pkPeITGzE07dVP/B/Lzj9rv2P8Xdtw1q26MwjiP5UwD80N1/DABm9k8AXgmg6yQvIjIuv3XOyh+85U2H4dFdCW78+p6H0d/v92M3jnPyRwP4aaq9tf2YiMhE+fcvHusLDcdLXjSHV79iBb7/wwbM7MRxj6sf4ziSX+qnYOb3QjO7EMCFAFCvrB72mEREAmZmLzl1Fn93yREH2vjL9x+OS/7XrrsxRUfz4ziS3wrg2FT7GAAP80LuvtndN7n7plpljrtFRIbqs5dtTJ759Cqe9YzawcdO3TSLmZrhq9cdMzVFv8ZxJP9tACeY2fEAHgJwLoDfGcM4RESWZGbV5z57Bl+86imZvr983+F449sfwTd/yyJ3T5ZYfaKMfJJ390UzeyuAr6AVobzC3e8e9ThERLr5n3+xYeHRXU0cdUR2inz68TVsOqmOP3zTYU1MwWmbseTk3f1LAL7U+woAmp0fmLZIPzwTWjZvU5FRm/pj6g8TmbmacX4c02rUTix1n/qa4Y7jhXBjHLGs7J3p3H8yjMpl3mTfF+6rHmbZX37028LlFzvROpsNl41mwiflR6wL2ns2hi+wV/s48CnIzWcikzHFD1P9mWU5qsjtnFgkRyhjevMy/RSZrBjn03Oy7NTmiORcpUH9YXtFZf/B+xyR5HY9Ctddjib/w6L3Li9S2cx8mYUiqwOMdi7FzFY/58Qa/vX6Y7su84F3rcNZv70V57/FZtx9f9cFJ8BYJnkRkUl18dvXPv7UY6tYuaL7nyzXHx7j3FetwutfY/OY8KN5lTUQEUn50s178KZzixN9b3vzGnzmn3ePYETLo0leRCTlO3cv3BZVHEnBfzOzwJbv7L913OMtokleRIQ0Penp1gszu8LMtpvZXanHPmxm95rZnWb2eTNbM6znokleRIQsotnTrUefBHA2PXYTgGe7+3PQquV18eBGH9IkLyJCmu493Xrh7t8A8Bg9dqO7H4g3/SdaXwodiilJ1ziQrnZIEUpLxys580j1Uzm5xYtnq1Aufb+17W7j7aaw5mx3NPBoMRxMOv0W768GfbUnw9hj9UmqHLnApYi5WmaqSuL+MHaXUGnh+aPCCqEL9PervooRFkUm6f2IKbqYXj6vqmSrTdvi2GNq2xXq48qQ1ShsZ7bFFTHTEUrKAPOyc5WwciRHJrnS5FyqPWNhH0coq1w2tQ8J/eOoWvgaFEYqUw9E9PlLxhBeSQZYZrkHvwvgM8Pa+JRM8iIio9PgCzZ0t97MtqTam929e61uYmbvB7AI4Ko+htcXTfIiIqTXUzEAdh5q7XozuwDAOQDO9CFe2EOTvIgIGXZBmvbV8d4D4CXuvneY+9IkLyJCmgM8J29mVwM4Ha1TO1sBfBCtNM0MgJus9Qeh/3T3/z6wnaZokhcRIY0Bnjxx9/OWePjywe0hnyZ5ERHSnOxyNH2ZjknePayEyBePXkz92KWLaXMckLNySUHVyXQ7E//jOCa3CyKXQX/BtjJxQj5pmKT7woXnDw+flC2G7ewFx8N2nDqsqeyl6GEStufX0YXUqfJmHo4xcknRTOwx6h6ZBMLoIm87W3WSK0d2jzIWRSa5n/fFVSjT2+bIZKaqJF18u+hi3Ok2xxo5Msn9faF1OVIZZT6wofTyHCNN6M3jOOZyLhLefTwD3+TYTMckLyIyQgsl+p6oJnkRETKM3w7GRZO8iAjROXkRkRJr9HNJuAmnSV5EhOhIXkSkxDIF1abYdEzyiQP75g82bWYm6LaFTnwrWgyfUkFyqzgGad37MtcbzlwUPL8/ve/MugURyn6+kFd43WNOmVKSLlrs7DxesNxlG3NhOwkLYhZenLsfRbHIdOXIfqtOxrR8OhZZo4tp1/qMTPLy6f5sVUmOTOZXkuSLcadjkZkIZUFksij2GKdyu5mLbxdsmyOWzdQbkumjDwn3z2c+ZMuXKF0jIlJeCzonLyJSXuOoYT8smuRFREjm9NMU0yQvIkIaXp6psTzPRERkQJr6xquISHnpdM0ymdkDAHYDaAJYPNTLZ4mIDINO1wzGr7v7zl4W9MVFNB/bdbAdR1TGtNKJO8Wz4VNqzobLZqK/RRny9G9tRbn4PrPu6eULM/h95Mv7ztgXlDFOUpHnZvgVhQzu92q4M4+WUcOVx5kpF9y9XZSDz/RTO11OuF4Js+l5pYOX6ucs/Ewqdz8bde9r9Ye5+X7KB8f0xveTg19K+vsAVTRylsxqWvdywdm+7pn61jgGf7E+na4RESkxfRlq+RzAjdY6FPt7d988pnGIiGSorMHynebuD5vZEWhdyPZed/9GegEzuxDAhQBQx9xS2xARGYoyVaEcy48rd3+4/f/tAD4P4JQlltns7pvcfVMVBSeCRUQGqImop9s0GPkozWyFma06cB/AWQDuGvU4RES6Sdx6uk2DcZyuORLA5631F/IKgE+7+w1jGIeIyJIUoVwGd/8xgJP6Xq/ZiYI1H90V9MXVTqnReDYsOxrNhufWojpHKvN/GqdDZ3lliFsbo3ZBDDJoFx0UFO07Hccs2G9h6eG8eCaXT6Zlkxr1V7iOcff1+z4uKig1nI5BVgsiknmRSQCop6KMmYhkv5FJijmmywlz31xcEJk0LnNMzyv1CR5G1LCbmPYVUXyzinDcSSoWyX/0zItbAsUlkw+FLhoiIlJinM2fZprkRUSIjuRFREqskZRnaizP7yQiIgOSwHq69cLMrjCz7WZ2V+qxdWZ2k5nd1/7/2mE9F03yIiKk6VFPtx59EsDZ9Nh7Adzs7icAuLndHgpN8iIipOFxT7detL/N/xg9/EoAV7bvXwngVQMbPJnKE0++GEbS/PEnDt6PZsIMX7UeRioX5yhSuRhGuyzvCw79VoosWj+nb6Dfs1jmtvMOWDwOX7+kSv28bqbtS98HYAVtjkxGmRhkp12JuVpj2K5Rf4360zFJjkjW4vyqlLMUg+SYZLrNfUURySLBHxDpzVhupDJdpZIjk5lxF1S0TEss/9iTj6CHcT3WEXzR6Uh33wYA7r6tXeJlKKZykhcRGaY+qlCuN7MtqfbmSSu4qEleRIQ0kp4n+Z2HeNGjR8xsY/sofiOA7YewjZ7onLyICEk86um2DNcDuKB9/wIAX1j2oLvQkbyICBnkl6HM7GoAp6N1amcrgA8CuATANWb2ZgAPAnjtwHZINMmLiJBB/uHV3c/r0nXmwHaSQ5O8iAhZLNFFQ0oxyTd37z54v1IPLzASU4SyuoIu9F233HY6Elh4CephxiAL9pW7X4oa9hv1DJbn2ChXpaRIJVeh9EoYpbNUv1EEMqJ1Y1o3jjnKSLHHVJsviN1PZBIIY5EVjl/mRCJ7aQcX215ORBJAk96Q4OLdvO0+I5UcWQ36uMokRSj5ouB5+yo6181VKYdBF/IWESkxVaEUESmxRU3yIiLlNS2X9uuFJnkREaLTNSIiJaYjeRGREtM5eRGREtOR/ATz+f1B257cF7SrT9SD9uIsXRm+at3b9Gpl8+i9j7PQMrbFcWbngfICnIvPK6GcWdZz25nywTG3O3npqCAHH3N2nbLuMe0rXU6431w8Z9/TGfGiXPxsHJbC5sx4pp0aC+fNuaxuUZldXj9J5dGpCjTiwiy75/ann0c9Cp9z3cLyynFOxh4Ic+lF3xWIRlByS5O8iEiJ6XSNiEiJ6UheRKTENMmLiJTYYu8XDZl4muRFREgmrDDFNMmLiJBhXBx8XIY2yZvZFQDOAbDd3Z/dfmwdgM8AOA7AAwBe5+67BrnfdNlhAIir4VOMqRRxbZZKD9fCX9MWZztvdlIL95WJKvJglvM54Y3xtvL6uY8Tk/2OOyg1nB+/zJYipsgkrR+l+rmUMEckK32UFgbCWGSFYnl5pYSXXj5VtngZpYSBbESQY49pRZMNRyq5nHC68DCX9822w3HUjGKkNM50bJIjk/ycC6XKBxeV+eUoZ1E881CU6Zz8ME88fRLA2fTYewHc7O4nALi53RYRmSjNJOrpNg2GNkp3/waAx+jhVwK4sn3/SgCvGtb+RUQOlbv1dJsGoz4nf6S7bwMAd99mZkd0W9DMLgRwIQDUMTei4YmIlOt0zcT+4dXdNwPYDACrbd3gT7qJiHRRpsv/jfqk0iNmthEA2v/fPuL9i4gUKtPpmlFP8tcDuKB9/wIAXxjx/kVECiVuPd2mwTAjlFcDOB3AejPbCuCDAC4BcI2ZvRnAgwBeO6z9H9B8LExoxtRfrYaPJFWq8lerpvrCNzVTrZF/ZBbFIA91WSxVaTJn1X6rUuacHONVM1UnM1UmiyKUnThctspkf+1q1L3SZGFkMlOVMmynY5H9Vp3sJzLJCpct2HZQKdKoUmTUX7VMji7mRShZE/nVNPvBkckI+VUrD4WX6ATx0CZ5dz+vS9eZw9qniMggJFMSj+zFxP7hVURkXKblVEwvNMmLiJAyna4pz+8kIiIDkiRRT7demNm7zOxuM7vLzK42s3rxWoOjSV5EhHiPtyJmdjSAtwPY1K7hFQM4dxhj7kana0REyIAz8BUAs2bWADAH4OFBbryXnU+8Z/zq03DTls8ebL80OvTkJUcqK9Xw8sa1mEOWnZIKbnQpZPogNOmXsIQX502nfo/qO9WYV5WyqGJlUVXKvCqU3EfPqajqZF6kMuJ4JbU5ipiJXOa0OTLJFRi5P6+yJPdxtHA5MnHLgotax/Rm8vLp2ORh8d6wL3Px7bCdeb1zLuRdiD5jDZp60t8w5bgli5EfUR0ETwYzybv7Q2b2EbQi4/sA3OjuNw5k4z3S6RoREeLe2w2t7wFtSd0uTG/HzNaiVZjxeABPAbDCzF4/yucyFUfyIiKj1Mfpmp3uvimn/zcA3O/uOwDAzK4D8CIAn1reCHunSV5EhA3unPyDAF5oZnNona45E8CWQW28F5rkRUSID+g0v7vfYmbXArgNwCKA29GurjsqmuRFRMgg0zXu/kG0aneNhSZ5ERFWom+8apIXESGDilBOgqmc5G9KPhu0l5ObT3b9PGhHs7NBO12KuB6Hb7x5+PI15sL+xXBTaNZo36kcvVPfcuTm3HuRU06Yc/BeCU9eFpUWzrQ5pJ83LBpXROvm5uS5lLDllymuZvo7bd5vkUxZXd52Km9eVO6Xc/DcX7Mwwz8XdUoAr4n3BH0rqDxwldbtBz/HhNLZYQJ/qX11/i01C15evmpTdSjzsSZ5EZHy0ukaEZES0+kaEZHyKlOpYU3yIiJMk7yISInpylAiIuU1hMKWY1OKST4dqew3Tpks0FXmtz0SNOPFTtSr7muDPkvC2sKVFWHd3cW65bdTkcvGXNCFhCOVXNKXU3npPlqVl3WKMXItUi6JHKyfiUjSziw/MsnS3yzk86D9nheNuOxuqp2NW3KkspnbH+WUE+bSuBwB5BK9cSYK2unnSORMQTngbLng8PO8Kpo/eH9NvC/oq4HjmfkveIM+SPPeffpo9hlBTL++Ma3b4P3QpgdZ6vkgHcmLiJRYic7Jq568iAhLeryNgJmdZmYr2vdfb2YfM7On9rq+JnkREebW2200Pg5gr5mdBODdAH4C4B97XVmTvIgIMe/tNiKL7u5oXWHqUne/FMCqXlcunOTN7K3tS1iJiPxi8B5vo7HbzC4GcD6AfzGzGEC1YJ2DejmSPwrAt83sGjM724zLRImIlMuEHcn/NoD9AH7X3X8G4GgAH+515cJ0jbt/wMz+GMBZAN4E4H+b2TUALnf3H3Vbz8yuAHAOgO3u/uz2Y38K4PcA7Ggv9j53/1Kvgx2FZP982N760MH70WO7gr7ZtWvCleszQdNr4Q/b5mFhWcqFwzvL7zkyzC3uXxP+LG2GaU0k9M6lY4/84XN+oDAyycunqlBydJO3vYxDAL5QQ7ZNw+LIpHVv87IxLctVJ/Mk/Z6LpZxplJMN5SqTHJlcRTHIFdH+/OVT/Suor0qvQVxwaDpP/enoKFehzFTeXAZ+TdgcxUYHYoIilO7+MzO7CsDzzewcAN9y98Gek2+fD/pZ+7YIYC2Aa83sr3NW+ySAs5d4/G/c/eT2baImeBERABN1usbMXgfgWwBeC+B1AG4xs9f0un7hkbyZvR3ABQB2ArgMwB+5e8PMIgD3ofXX3gx3/4aZHdfrQEREJsWEfeP1/QCe7+7bAcDMNgD4vwCu7WXlXr4MtR7Af3P3n6QfdPek/atDv95qZm9A64rlF7n7rqIVRERGarK+DBUdmODbHkUfycjCBd39T3iCT/Xd0+uO2j4O4JcBnAxgG4CPdlvQzC40sy1mtmXHjh3dFhMRGbwJOl0D4Mtm9hUze6OZvRHAvwDo+VT3SHPy7v6IuzfdPQHwCQCn5Cy72d03ufumDRs2jG6QIvILzxLr6TYiDuDvATwHwEkANvez8khr15jZRnff1m6+GsBdo9y/iEhPJut0zUvd/T0ArjvwgJn9GYD39LLy0CZ5M7sawOkA1pvZVgAfBHC6mZ2M1kv4AIDfH/R++SLfbFkX/d67N2j7QhhJs2r4clolbFceXxm0471rUq3wC2xJXPDWhGnN4GplHIlknA7jyGQm/ZZu9xuZpP7+LtxNET8aZzWmi1jH3S96nYlX9vmvmC9MncYXlmYcTeQIZpxThTKhF5ArLlapkiS308tzZcgqjYsviM3VIOscO/VOdHE3jStJ8itW5kUs+TXgt4r7V0SDj1COMAPffQxmfwDgDwE8zczuTHWtAvDvvW5naJO8u5+3xMOXD2t/IiIDMwGTPIBPA/gygA8BeG/q8d3u/livG1GpYRERMgkRSnd/HMDjAJY6YO6ZJnkRETYZR/IDoSqUIiJkkLVrzGyNmV1rZvea2T1mdupwRx/SkbyICBvskfylAG5w99eYWQ3AXNEKg6RJXkSEDCpdY2arAbwYwBsBwN0XAAyholp3Ol0jIsIG943Xp6FVdfcfzOx2M7vswKX8RuUX7ki+KEf/stnzwwc8Sd0N31WLKFhM/b64GPbPh2WMbU+nXX0irCU8syr8+esVaueVC+acPH0zz2ic3M+HMUFl3EzIvuCTzmWPc9bP9BHOzRdl3ytBTj6MSxSVr+27nHDOutyO6JIMe5u1rtvijP0MlQvmksmZsaTy6lyGOKZ8eSYHb5yrj6i/8xrGCLfNmgmV2abcfPp7CBFl7rmc8pooLLdc51z9IPR+JL/ezLak2pvdPf2N1AqA5wF4m7vfYmaXohWH/OOBjLMHv3CTvIhIkT4ilDvdfVNO/1YAW939lnb7WoSZ96HT6RoRETKodE37Sk4/NbNnth86E8D3hjj0DB3Ji4iwwaZr3gbgqnay5sdoXWFvZDTJi4iQQX7j1d3vAJB3SmeoNMmLiLASfeNVk7yICNMkX15f2fd/uvadve73wgc4QtkII5PepGgXRxdTEctoX7hudU81XHSGYnlUFzYdoSxK/3FUMeEYJEcqvct9AE5/u7co/D2X98V/rAr78/9lZcdNbeS30xapHnNCmdRs5LIzNo41Mi7py2J6vRPrvjyX1X28mf9lyYUofF611PprorBUdt0p4pspY8yRSXqN0q8/RTiaCLfd8DCu2aRtNVK5X456rrBw3VUU/ZwbQjWxSShQNiia5EVEyCTUkx8UTfIiIkyTvIhIiWmSFxEpL52uEREpMU3yIiJlpkn+F9MNj30it//lR/5B+ABFKjORy3TXQrhsvC/McNWeDCNnmQhlXhWigsKRRpHJhKpYpgpxZtblVKhn4oA0Tq5wmRo3/7tKV1AEgAa19za6V28EgIW480RqcRhF5IqVmYqWlKFLV7SsROG2iipBsuy+u+f1djfD6qQcqeQ2i1ORVo528nOuUj9HJitU3rSSWjzmD1lEn2eEFVjrFvYvpLZdQ/ic5uj1XsG7KoisHgpFKEVEykxH8iIi5aVz8iIiZaZJXkSkvHROXkSkzHQkLyJSXjon3wMzOxbAPwI4Cq3ydpvd/VIzWwfgMwCOA/AAgNe5+65hjWOUvvzIx4P22esvDNpmOTlHyiZGjfD3xXg/VQjcy1UoO23+fHJVyqSW358ZZjpC2aSqkpxey1SwpHEWtINNUd9iMxzYPg8rdTYpYrkvFZusUgyP44OVqHtkEggjmNzHkcrMtgsuIp5+FlwNk1+DBlXPnKfXoOocqexcBLtK0cQaV53kKpMFVweNU8vHyI9brozD12QdRSwbORcC53FwdcxkCIfdmYvdT7FhXuN1EcBF7v4sAC8E8BYzOxGti9je7O4nALgZI76orYhIIe/xNgWGNsm7+zZ3v619fzeAewAcDeCVAK5sL3YlgFcNawwiIodiUBfyngQjOSdvZscBeC6AWwAc6e7bgNYPAjM7YhRjEBHpVZnSNcM8XQMAMLOVAD4H4J3u/kQf611oZlvMbMuOHTuGN0AREabTNb0xsypaE/xV7n5d++FHzGxju38jgO1Lrevum919k7tv2rBhwzCHKSISKNPpmqFN8mZmAC4HcI+7fyzVdT2AC9r3LwDwhWGNQUTkkJToSH6Y5+RPA3A+gO+a2R3tx94H4BIA15jZmwE8COC1QxyDiEjfyhShHNok7+7/hkyR24POHNZ+J4lV+nh5uYYvta1JOfoFztV37scVyqLzMOj3t0yZ4kz12s72uExuwm/x/nBjmcx+TqniRc7QV8KBNCknb/T78v4ofKLpMruVmLPtlBGnUsS1mMpEp0QVeg3oBYyMc/PhvmYoI15N9c/QfnlZ3lZMWfeYM/ipsdQy4wqambx5I5O5D6XLC3Muvmr5n/0ZC/P9De/+evO4EnrOCX+oBmBaTsX0Qt94FRFhmuRFRMqr4FosU2XoEUoRkWkz6HSNmcVmdruZfXF4o16ajuRFRNjgz/O/A61v/a8e9IaL6EheRIRY0tutp22ZHQPgvwC4bJhj7kZH8iIiZMDpmr8F8G4Aqwa61R5pkh8mjlA26a856V8JudTwYp+fstTifISR+SNSUX/36r+Z0sGURIRn2lQWtkExyWqnv0nRxGaFookUg+SSyEaDiSJf8j4AVCmeuVjJ/6U2XT64RnHLiPZbi/JjkLNxo2t/PeK+sF03jl82c9u11PJcApk16DPYRLgvPnCtpj4oc/SZ4dLDRWWL8yKXHK8cSVmZ3k/XrDezLan2ZnfffKBhZucA2O7ut5rZ6YMbYO80yYuIkD6O5He6+6ac/tMA/KaZvQJAHcBqM/uUu79+mUPsmc7Ji4iQQZ2Td/eL3f0Ydz8OwLkAvjrKCR7QkbyISJbKGoiIlNgQ5nh3/zqArw9+y/k0yYuIEBUoExEpMRUok97M1sP2/oWwHaUq9/EfcQqOJDymdqqkYKaqJOEPcOYPSHl/UCpY1yliyU8jalKEMlVd02NauBou24ypCiUvH3WPXFZqVNEyCbcdZ8bdPUfKUcQ6RSJXxuH7PEttjlSmY5PZyGTYzuyblq9SxJKrhqY16DnOF4QTm5StraeWn8l8MHI31ReOVy5SdczmUM6tDH6T46JJXkSE2BDKF4+LJnkREcLXb5hmmuRFRFh55nhN8iIiGTpdIyJSXopQioiUWK9lhKeBJvllOPV3Phq0v/npi4L2l+/766D98qf/UddtWZLkthnHJNMX687GK6nN8bblxN0K1s3kjbmdjjLStjLXkY7yj674wt7piGVMFSy5CmW9QtUdqdJkPXWBbY5MclXJ1ZV9QXsuE6Hk2GNnX3xhbq4qmXeh7qXW5+XzcKSySGyd5asWfuiKqk6y/d55TfhC3byt3Un4+u0dxqkVna4RESkvna4RESkzHcmLiJSYzsmLiJSXvvEqIlJmBcGHaaJJXkSElWeOH97l/8zsWDP7mpndY2Z3m9k72o//qZk9ZGZ3tG+vGNYYREQOhbn3dJsGwzySXwRwkbvfZmarANxqZje1+/7G3T8yxH0PxUv+64eDdo36f/2svwra8f4ww1xZORuu0Oz0c0GkaD7MbcfzlB1uUBg+JeFcPC9a9KM9r5wwb7soB8/Z60zt4Zx16WjKEg7SUy6ecvRxpbOBWjV8PWcyOXnKvlepfHC1k3XnXPxcph3m4lfG80E7L9vOuXYuFcz92Vx8uHwz9SWJBn1hokFvZpPWjWnfddr3XKoE8IxV0Y90Lh4AnvTOa9ak95Vz8I8mM0H7+b/0AG19AHWOdbqmmLtvA7CtfX+3md0D4Ohh7U9EZGCm5Ci9F0M7XZNmZscBeC6AW9oPvdXM7jSzK8xs7SjGICLSs6TH2xQY+iRvZisBfA7AO939CQAfB/DLAE5G60j/o13Wu9DMtpjZlh07dgx7mCIiB5XpnPxQJ3kzq6I1wV/l7tcBgLs/4u5Nd08AfALAKUut6+6b3X2Tu2/asGHDMIcpIhJqJr3dpsAw0zUG4HIA97j7x1KPb0wt9moAdw1rDCIih8S9t9sUGGa65jQA5wP4rpnd0X7sfQDOM7OT0cpSPADg94c4BhGR/k3JBN6LYaZr/g1LZ5m+NKx9DltCUUQuo/u1G98TtM8840NBu7E2jFBa6te9aCGM1RnHL/eGkbOZOHxpm/VOhK1ZC/ucf6vkEr4cVcwpRcyvAbjN/zY49bic3x0Ltm207Sj1BlWo1HCdIpIrqmHscY4ilauqnRjkqkoYiVxb2RO011d2h9uK9gdtjjnmyZYeppLIFMesWzjudOSySS/+ApUHrlFt53pErxl9MGb7jE2mRfTmzaQ+SAm9Po8n4bjmkxF8h3NKTsX0Qt94FRFhmSOj6aVJXkSEleh0zUhy8iIiU2VA6Zpu5V1GSUfyIiJscEfyS5Z3cffvDWoHRTTJi4iwAU3yOeVdNMmLiIzNEAqULVHeZSQ0yZMXvzKsNJmkoorN2fBPGN+8+qLcbd381YsPeRxnnfLnuf2VPWFUrr4zHUkL39bFehhXSyj5lmlXuvdn4pVclZK6OTLpFW77kvcBADG1qd8oFhlRu5qqNFmvhNHDWYpIcmRyrhJGKtOxybWVvUEfRyY3VJ4I143CyGVC8cFm6kXkmGONIpN1as9FC7n9eYqrTIbjXBWFH5SqHfr0EVPedc46NV2blGxZAFel5PqvQ9D7JL/ezLak2pvdfTMvtER5l5HRJC8iwrgkdnc73X1T3gJLlXcZJU3yIiLEB5ST71beZZQUoRQRYYMrUHagvMsZ47oano7kRUTYgP7wmlPeZWQ0yYuIsBJ941WTvIgIcV3jdXq89Nf+ImgnFa7eGFbi89mwnV4+qY7ut65kNv+t4Qt/V/Z2onP1XfQcqSplcyZsN2apvap75JIvCt5PRHKpdlLv/GOyCv3Dom0bVUWMa1Spky/Onao0yRfqXklVJw+r7QvaXGlyffXJzrJxGKHkyOSGOIxUrqCYI5tPvWjZCGX4nNZQRctV9BrNUWXJZip+2Cj4Y2KdIpGro9kuSy5flPPnwCZfnJwilC87fgTfI1IVShGRElMVShGR8vLec/ITT5O8iAjxJl9ZZ3ppkhcRSVmJw57nSROWqeERcneswtpfHdGwDpm+DCUikrIa67CjVTgy1zb8BOuwYQQjWh4dyYuIpDyM+4/ejV0PrfejEPGFg9uavogHcR+exONrRzy8vulIXkQkxd0fXo+NeBj3d13mQdyHp+A4uPvPRzeyQzOVR/IveENY5ydqdP4SXnsi/IOJrQyf4tdveE/QPu01HwnazXr4cy9JxY45Ez5MjdXV3P54Lz3PVBogWqD4VxKeW4wW8+sFL1JuPv28uSwx5+YTzsFXKSdfC8dmqXZEOXc+iMrk5CvdSwsDwEy1890BzsWvroa5+DXUXh2H7bWVPZ1lKSefbYcZ+xWUdeePUQOdcTY8fO25/O+6aHBZ9oTy6HnZ9VHiEsarovC9+H8PPCNo/zyZG/gY7sc9q1fisCeO8l9CxcIP/YLP4xE8hCfx85mB73gIJuNdFRGZIO6++xg8DT/BDzJ9P8Y9OA7PhLvnf8ttQmiSFxFZwr24vbYT27DfO79J7PHd2I2f4y7cMjVz59QMVERklNy9cTxOxI9Sl2P9Ib6Lx/Hoi92np4LZVJ6TFxEZhTvxH9FabEie9MexiFb9I3f/1zEPqy+a5EVEunB3N7MX3Ifv3tJEAz/Ho78y7jH1S6drRERyuPu3KqhiJQ6Du39/3OPpl03DqSUz2wHgJwDWA9g55uEsRePq36SOTePqzySO66nuPvlfRR2RqZjkDzCzLUVXRh8Hjat/kzo2jas/kzou6dDpGhGREtMkLyJSYtM2yW8e9wC60Lj6N6lj07j6M6njkrapOicvIiL9mbYjeRER6cNUTPJmdraZfd/Mfmhm7x3zWK4ws+1mdlfqsXVmdpOZ3df+/8hrTJvZsWb2NTO7x8zuNrN3TMLYzKxuZt8ys++0x/VnkzCu1PhiM7vdzL44YeN6wMy+a2Z3mNmWSRmbma0xs2vN7N72Z+3USRiXdDfxk7yZxQD+DsDLAZwI4DwzO3GMQ/okgLPpsfcCuNndTwBwc7s9aosALnL3ZwF4IYC3tF+ncY9tP4Az3P0kACcDONvMXjgB4zrgHQDuSbUnZVwA8OvufnIqojgJY7sUwA3u/isATkLrtZuEcUk37j7RNwCnAvhKqn0xgIvHPKbjANyVan8fwMb2/Y0Avj8Br9sXALx0ksYGYA7AbQBeMAnjAnAMWpPSGQC+OEnvJYAHAKynx8Y6NgCrAdyP9t/yJmVcuuXfJv5IHsDRAH6aam9tPzZJjnT3bQDQ/v8R4xyMmR0H4LkAbsEEjK19SuQOANsB3OTuEzEuAH8L4N1AcPWMSRgXADiAG83sVjO7cELG9jQAOwD8Q/sU12VmtmICxiU5pmGSX+qS6YoEdWFmKwF8DsA73f2JcY8HANy96e4no3XkfIqZPXvMQ4KZnQNgu7vfOu6xdHGauz8PrdOUbzGzF497QGgVNHwegI+7+3MB7IFOzUy8aZjktwI4NtU+BsDDYxpLN4+Y2UYAaP9/+zgGYWZVtCb4q9z9ukkaGwB463qYX0frbxrjHtdpAH7TzB4A8E8AzjCzT03AuAC0rjPa/v92AJ8HcMoEjG0rgK3t38QA4Fq0Jv1xj0tyTMMk/20AJ5jZ8WZWA3AugOvHPCZ2PYAL2vcvQOt8+EiZmQG4HMA97p6+CO5Yx2ZmG8xsTfv+LIDfAHDvuMfl7he7+zHufhxan6mvuvvrxz0uADCzFWa26sB9AGcBuGvcY3P3nwH4qZk9s/3QmQC+N+5xSYFx/1GglxuAVwD4AYAfAXj/mMdyNYBtABpoHdm8GcDhaP0B7772/9eNYVy/htZprDsB3NG+vWLcYwPwHAC3t8d1F4A/aT8+9tcsNcbT0fnD69jHhda57++0b3cf+MxPyNhOBrCl/X7+M4C1kzAu3brf9I1XEZESm4bTNSIicog0yYuIlJgmeRGREtMkLyJSYprkRURKTJO8iEiJaZIXESkxTfIyVczs+WZ2Z7tO/Yp2jfqx18IRmVT6MpRMHTP7HwDqAGbRqqXyoTEPSWRiaZKXqdOuYfRtAPMAXuTuzTEPSWRi6XSNTKN1AFYCWIXWEb2IdKEjeZk6ZnY9WuWBj0frikRvHfOQRCZWZdwDEOmHmb0BwKK7f7p9/d//MLMz3P2r4x6byCTSkbyISInpnLyISIlpkhcRKTFN8iIiJaZJXkSkxDTJi4iUmCZ5EZES0yQvIlJimuRFRErs/wM0MwcGcenRgQAAAABJRU5ErkJggg==\n",
"text/plain": [
"