{ "cells": [ { "cell_type": "markdown", "metadata": { "lines_to_next_cell": 0, "toc": true }, "source": [ "

Table of Contents

\n", "
" ] }, { "cell_type": "markdown", "metadata": { "lines_to_next_cell": 0 }, "source": [ "# Analyzing precip and temperature data\n", "\n", "## Learning objectives\n", "\n", "1. Document your project data by capturing important information\n", " as metadata in a json file\n", "\n", "1. Use the groupby method for a pandas dataframe to group observations\n", " by season\n", "\n", "1. Find the probability density function (histogram) of temperature and\n", " precipitation data and determine possible statistical models (normal distribution,\n", " exponential distribution) that best describe the variation in the data\n", "\n", "1. Find best fit parameters for statistical models that describe the precipitation\n", " and temperature data, assuming that temperature and precipitation are\n", " independent of each other.\n", "\n", "1. Use groupby to find the median temperature for each decade and each season\n", " and determine whether the seasonal temperatures have increased at YVR\n", " since the 1930's\n", "\n", "## Introduction\n", "\n", "We have about 30,000 days worth of YVR precipitation and temperature\n", "data generated by the 10-pandas1, 10-pandas2, and 10-pandas3 notebooks, these\n", "gave you a new csv file: data/processed/weather_YVR.csv. We will use\n", "that data set to estimate statistical distributions that could produce realistic\n", "simulated rainfall and temperature inputs to a model. It makes use of several\n", "new python modules/methods:\n", "\n", "1. The **context** module: We will write a new module called context.py that\n", " will be used to locate folders and libraries associated with a project, such as\n", " `weather_YVR.csv`.\n", "1. **The json module**. This reads and writes files written\n", " in \"javascript object notation\" i.e.\n", " [json](https://en.wikipedia.org/wiki/JSON). We will use it to\n", " save \"metadata\" or data about our data.\n", "\n", "1. **pandas.DataFrame.groupby**\n", " We will use groupby to split the dataframe into 4 seasons, each corresponding\n", " to 3 months of the year: winter (month initials djf), spring (mam),\n", " summer (jja), and fall (son). This is the tip of the iceberg of what\n", " you can do with pandas split/apply/combine functions. For more examples\n", " see [the Python Data Science Handbook](https://github.com/jakevdp/PythonDataScienceHandbook/blob/master/notebooks/03.08-Aggregation-and-Grouping.ipynb)\n", "\n", "1. The [scipy.stats](https://docs.scipy.org/doc/scipy/reference/stats.html) module. As you will see\n", " below, the temperature data is best fit with a a Gaussian (normal) distribution\n", " using the norm.fit:\n", " - [scipy.stats.norm.fit](https://docs.scipy.org/doc/scipy/reference/generated/scipy.stats.norm.html#scipy.stats.norm)\n", "\n", " Precipitation, on the other hand, follows an exponential distribution, which\n", " can be fit with:\n", " - [scipy.stats.expon.fit](https://docs.scipy.org/doc/scipy/reference/generated/scipy.stats.expon.html#scipy.stats.expon)" ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "in context.py, setting root_dir to /Users/phil/repos/eosc213_students/notebooks/pandas\n", "******************************\n", "context imported. Front of path:\n", "/Users/phil/repos/eosc213_students/notebooks/pandas\n", "/Users/phil/mini37/envs/e213/lib/python36.zip\n", "******************************\n", "\n" ] } ], "source": [ "import json\n", "\n", "import context\n", "import pandas as pd\n", "import scipy.stats as st\n", "from matplotlib import pyplot as plt\n", "import pprint\n", "import numpy as np" ] }, { "cell_type": "markdown", "metadata": { "lines_to_next_cell": 0 }, "source": [ "## The context module (repeated from 10-pandas1)\n", "\n", "As your project grows more complicated, it's good to have a central\n", "module that keeps track of important files and sets your scripts\n", "up so that they can import functions and classes from you modules.\n", "If you were planning to distribute your project using conda, then\n", "you would need to write an installation script, which is a fair\n", "amount of work. At this stage, it's easier and more flexible to\n", "store that information in a file that travels along with your notebook.\n", "By convention, this file is called:\n", "[context.py](https://github.com/phaustin/eosc213_students/blob/master/notebooks/pandas/context.py)\n", "\n", "Clicking on that link shows you the source -- the code is executed when\n", "you do `import context` and does the following:\n", "\n", "1. Defines Path objects that give the location of the raw and processed data folders\n", "\n", "1. Puts the notebooks/pandas folder on python's list of places to look for\n", " library modules (sys.path). We will use this when we start extracting\n", " code from notebooks into libraries" ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "/Users/phil/repos/eosc213_students/notebooks/pandas/data/processed\n" ] } ], "source": [ "#\n", "# Here is the path to the processed csv file\n", "#\n", "print(context.processed_dir)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Reading the processed data\n", "\n", "The next cell uses [\"globbing\"](https://en.wikipedia.org/wiki/Glob_(programming)) to\n", "find every csv file in all folders below `data/processed`. (See\n", "[this optional tutorial](https://realpython.com/python-pathlib/) for more information)." ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "yvr files: \n", "[PosixPath('/Users/phil/repos/eosc213_students/notebooks/pandas/data/processed/weather_daily_YVR_1938-2017_all.csv'),\n", " PosixPath('/Users/phil/repos/eosc213_students/notebooks/pandas/data/processed/weather_daily_YVR_1938-2019_all.csv'),\n", " PosixPath('/Users/phil/repos/eosc213_students/notebooks/pandas/data/processed/weather_YVR.csv')]\n" ] } ], "source": [ "yvr_files = list(context.processed_dir.glob(\"**/*YVR*csv\"))\n", "print(f\"yvr files: \\n{pprint.pformat(yvr_files)}\")" ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [], "source": [ "yvr_file = context.processed_dir / \"weather_YVR.csv\"" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Metadata is data about data\n", "\n", "Since we care about this dataset, it's a good idea to save\n", "details for future reference. We can't put this kind of information\n", "into a csv file, but there are many file formats that can store\n", "this type of descriptive metadata -- one of the most popular is\n", "[json](https://en.wikipedia.org/wiki/JSON). In the next cell we\n", "write the data into a nested dictionary called meta_dict, and\n", "dump that information into a new json file called metadata.json" ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [], "source": [ "meta_data = dict()\n", "\n", "file_info = dict()\n", "file_info[\"written_by\"] = \"10-pandas3-process-daily-data.ipynb\"\n", "file_info[\"contact\"] = \"paustin@eoas.ubc.ca\"\n", "file_info[\n", " \"description\"\n", "] = \"EC canada YVR data downloaded by 10-pandas2-download-daily-data.ipynb\"\n", "\n", "meta_data[\"weather_daily_YVR_1938-2017_all.csv\"] = file_info\n", "\n", "file_info = dict()\n", "file_info[\"written_by\"] = \"10-pandas3-process-daily-data.ipynb\"\n", "file_info[\"contact\"] = \"paustin@eoas.ubc.ca\"\n", "file_info[\n", " \"description\"\n", "] = \"Reduced YVR data subsetted by 10-pandas3-process-daily-data.ipynb\"\n", "\n", "meta_data[\"weather_YVR.csv\"] = file_info\n", "\n", "history_file = context.processed_dir / \"history_metadata.json\"\n", "with open(history_file, \"w\") as f:\n", " json.dump(meta_data, f, indent=4)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Take a minute and look at history_metadata.json with jupyter to see how the metadata is written out." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Read weather_YVR.csv into a dataframe\n", "\n", "Here's the dataframe produced by 10-pandas3\n", "\n" ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "there are 29632 days in the dataframe\n" ] }, { "data": { "text/html": [ "
\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
DateYearMonthDayT_mean (C)T_high (C)T_low (C)Rain (mm)Snow (cm)Total Precip (mm)
01938-01-011938114.49.4-0.6NaNNaN0.3
11938-01-021938124.57.21.7NaNNaN0.5
21938-01-031938131.77.2-3.90.00.00.0
31938-01-041938142.27.2-2.80.00.00.0
41938-01-051938152.27.2-2.80.00.00.0
\n", "
" ], "text/plain": [ " Date Year Month Day T_mean (C) T_high (C) T_low (C) Rain (mm) \\\n", "0 1938-01-01 1938 1 1 4.4 9.4 -0.6 NaN \n", "1 1938-01-02 1938 1 2 4.5 7.2 1.7 NaN \n", "2 1938-01-03 1938 1 3 1.7 7.2 -3.9 0.0 \n", "3 1938-01-04 1938 1 4 2.2 7.2 -2.8 0.0 \n", "4 1938-01-05 1938 1 5 2.2 7.2 -2.8 0.0 \n", "\n", " Snow (cm) Total Precip (mm) \n", "0 NaN 0.3 \n", "1 NaN 0.5 \n", "2 0.0 0.0 \n", "3 0.0 0.0 \n", "4 0.0 0.0 " ] }, "execution_count": 6, "metadata": {}, "output_type": "execute_result" } ], "source": [ "yvr_df = pd.read_csv(yvr_file)\n", "print(f\"there are {len(yvr_df)} days in the dataframe\")\n", "yvr_df.head()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Handling missing data\n", "\n", "What should we do about missing vlaues like the Rain and Snow totals for Jan./Feb. 1938?\n", "This is a judgement call -- we could through away the whole row, or do something\n", "more complicated like masking the missing values in calculations where they are referenced.\n", "For this notebook we'll fill the NaNs with zeros -- not a great idea since 0 is an\n", "actual temperature, but we can see if that makes a difference later" ] }, { "cell_type": "code", "execution_count": 7, "metadata": {}, "outputs": [], "source": [ "yvr_df.fillna(0.,inplace=True)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Use apply to tag the 29,190 days with their season\n", "\n", "We know that the seasons are quite different, and if we are interested in\n", "generating daily data we have to take that into account. In the next cell\n", "we set up a lookup table (dictionary) called season_table that maps the\n", "month number to one of four seasons. With this table we can write\n", "a function called find_season that takes a row of the dataframe and\n", "returns the season for that row" ] }, { "cell_type": "code", "execution_count": 8, "metadata": {}, "outputs": [], "source": [ "season_table = dict()\n", "for i in range(12):\n", " m = i + 1\n", " if m in [12, 1, 2]:\n", " # winter\n", " season_table[m] = \"djf\"\n", " elif m in [3, 4, 5]:\n", " # spring\n", " season_table[m] = \"mam\"\n", " elif m in [6, 7, 8]:\n", " # summer\n", " season_table[m] = \"jja\"\n", " else:\n", " # fall\n", " season_table[m] = \"son\"\n", "\n", "\n", "def find_season(row, season_table):\n", " month = row[\"Month\"]\n", " return season_table[month]" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### creating a new column with apply\n", "\n", "One of the most useful pandas dataframe methods is [apply](https://pandas.pydata.org/pandas-docs/stable/reference/api/pandas.DataFrame.apply.html). This method has lots of bells and whistles, all we need here is to apply the `find_season` function to every\n", "row of the dataframe, which means that we want to return a new column with the same number of rows\n", "as the other columns. We do that by tell apply that we want to return `axis=1` (rows are axis=0)" ] }, { "cell_type": "code", "execution_count": 9, "metadata": {}, "outputs": [], "source": [ "season = yvr_df.apply(find_season, args=(season_table,), axis=1)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now add that column to the dataframe" ] }, { "cell_type": "code", "execution_count": 10, "metadata": {}, "outputs": [ { "data": { "text/html": [ "
\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
DateYearMonthDayT_mean (C)T_high (C)T_low (C)Rain (mm)Snow (cm)Total Precip (mm)season
01938-01-011938114.49.4-0.60.00.00.3djf
11938-01-021938124.57.21.70.00.00.5djf
21938-01-031938131.77.2-3.90.00.00.0djf
31938-01-041938142.27.2-2.80.00.00.0djf
41938-01-051938152.27.2-2.80.00.00.0djf
\n", "
" ], "text/plain": [ " Date Year Month Day T_mean (C) T_high (C) T_low (C) Rain (mm) \\\n", "0 1938-01-01 1938 1 1 4.4 9.4 -0.6 0.0 \n", "1 1938-01-02 1938 1 2 4.5 7.2 1.7 0.0 \n", "2 1938-01-03 1938 1 3 1.7 7.2 -3.9 0.0 \n", "3 1938-01-04 1938 1 4 2.2 7.2 -2.8 0.0 \n", "4 1938-01-05 1938 1 5 2.2 7.2 -2.8 0.0 \n", "\n", " Snow (cm) Total Precip (mm) season \n", "0 0.0 0.3 djf \n", "1 0.0 0.5 djf \n", "2 0.0 0.0 djf \n", "3 0.0 0.0 djf \n", "4 0.0 0.0 djf " ] }, "execution_count": 10, "metadata": {}, "output_type": "execute_result" } ], "source": [ "yvr_df[\"season\"] = season\n", "yvr_df.head()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Grouping the data by season\n", "\n", "The next cell uses the groupby method to sort all of the\n", "rows into 4 seasons. The resulting group object (`seasons`) has\n", "4 dataframes inside it, keyed by the season marker djf, mam, etc." ] }, { "cell_type": "code", "execution_count": 11, "metadata": {}, "outputs": [], "source": [ "seasons = yvr_df.groupby(\"season\")" ] }, { "cell_type": "code", "execution_count": 12, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "\n" ] }, { "data": { "text/html": [ "
\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
DateYearMonthDayT_mean (C)T_high (C)T_low (C)Rain (mm)Snow (cm)Total Precip (mm)season
01938-01-011938114.49.4-0.60.00.00.3djf
11938-01-021938124.57.21.70.00.00.5djf
21938-01-031938131.77.2-3.90.00.00.0djf
31938-01-041938142.27.2-2.80.00.00.0djf
41938-01-051938152.27.2-2.80.00.00.0djf
591938-03-011938315.09.40.63.60.03.6mam
601938-03-021938324.78.31.10.50.00.5mam
611938-03-031938338.112.83.30.00.00.0mam
621938-03-041938348.613.33.90.00.00.0mam
631938-03-051938356.210.61.71.30.01.3mam
1511938-06-0119386111.718.94.40.00.00.0jja
1521938-06-0219386213.620.07.20.00.00.0jja
1531938-06-0319386314.721.18.30.00.00.0jja
1541938-06-0419386416.725.08.30.00.00.0jja
1551938-06-0519386516.422.210.60.00.00.0jja
2431938-09-0119389116.220.611.71.30.01.3son
2441938-09-0219389215.318.911.72.30.02.3son
2451938-09-0319389315.019.410.60.00.00.0son
2461938-09-0419389414.817.811.73.30.03.3son
2471938-09-0519389515.619.411.71.30.01.3son
\n", "
" ], "text/plain": [ " Date Year Month Day T_mean (C) T_high (C) T_low (C) \\\n", "0 1938-01-01 1938 1 1 4.4 9.4 -0.6 \n", "1 1938-01-02 1938 1 2 4.5 7.2 1.7 \n", "2 1938-01-03 1938 1 3 1.7 7.2 -3.9 \n", "3 1938-01-04 1938 1 4 2.2 7.2 -2.8 \n", "4 1938-01-05 1938 1 5 2.2 7.2 -2.8 \n", "59 1938-03-01 1938 3 1 5.0 9.4 0.6 \n", "60 1938-03-02 1938 3 2 4.7 8.3 1.1 \n", "61 1938-03-03 1938 3 3 8.1 12.8 3.3 \n", "62 1938-03-04 1938 3 4 8.6 13.3 3.9 \n", "63 1938-03-05 1938 3 5 6.2 10.6 1.7 \n", "151 1938-06-01 1938 6 1 11.7 18.9 4.4 \n", "152 1938-06-02 1938 6 2 13.6 20.0 7.2 \n", "153 1938-06-03 1938 6 3 14.7 21.1 8.3 \n", "154 1938-06-04 1938 6 4 16.7 25.0 8.3 \n", "155 1938-06-05 1938 6 5 16.4 22.2 10.6 \n", "243 1938-09-01 1938 9 1 16.2 20.6 11.7 \n", "244 1938-09-02 1938 9 2 15.3 18.9 11.7 \n", "245 1938-09-03 1938 9 3 15.0 19.4 10.6 \n", "246 1938-09-04 1938 9 4 14.8 17.8 11.7 \n", "247 1938-09-05 1938 9 5 15.6 19.4 11.7 \n", "\n", " Rain (mm) Snow (cm) Total Precip (mm) season \n", "0 0.0 0.0 0.3 djf \n", "1 0.0 0.0 0.5 djf \n", "2 0.0 0.0 0.0 djf \n", "3 0.0 0.0 0.0 djf \n", "4 0.0 0.0 0.0 djf \n", "59 3.6 0.0 3.6 mam \n", "60 0.5 0.0 0.5 mam \n", "61 0.0 0.0 0.0 mam \n", "62 0.0 0.0 0.0 mam \n", "63 1.3 0.0 1.3 mam \n", "151 0.0 0.0 0.0 jja \n", "152 0.0 0.0 0.0 jja \n", "153 0.0 0.0 0.0 jja \n", "154 0.0 0.0 0.0 jja \n", "155 0.0 0.0 0.0 jja \n", "243 1.3 0.0 1.3 son \n", "244 2.3 0.0 2.3 son \n", "245 0.0 0.0 0.0 son \n", "246 3.3 0.0 3.3 son \n", "247 1.3 0.0 1.3 son " ] }, "execution_count": 12, "metadata": {}, "output_type": "execute_result" } ], "source": [ "#\n", "# this creates a new object of type DataFrameGroupBy\n", "#\n", "print(type(seasons))\n", "seasons.head()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Use a dictionary comprehension put the dataframes into a dictionary\n", "\n", "In the next cell we create a new dictionary by iterating over the DataFrameGroupBy object. We\n", "will use a [dictionary comprehension](https://jakevdp.github.io/WhirlwindTourOfPython/11-list-comprehensions.html),\n", "to simplify the loop -- note that we wind up with 4 dataframes in the dictionary,\n", "one for each season. We do this because a dictionary of dataframes is simpler to work with than\n", "the DataFrameGroupBy object, which somewhat obscures the season keys." ] }, { "cell_type": "code", "execution_count": 13, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "dict_keys(['djf', 'jja', 'mam', 'son'])" ] }, "execution_count": 13, "metadata": {}, "output_type": "execute_result" } ], "source": [ "season_dict = {key: value for key, value in seasons}\n", "season_dict.keys()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Here is how you get the fall (son) dataframe from the dictionary" ] }, { "cell_type": "code", "execution_count": 14, "metadata": {}, "outputs": [ { "data": { "text/html": [ "
\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
DateYearMonthDayT_mean (C)T_high (C)T_low (C)Rain (mm)Snow (cm)Total Precip (mm)season
2431938-09-0119389116.220.611.71.30.01.3son
2441938-09-0219389215.318.911.72.30.02.3son
2451938-09-0319389315.019.410.60.00.00.0son
2461938-09-0419389414.817.811.73.30.03.3son
2471938-09-0519389515.619.411.71.30.01.3son
\n", "
" ], "text/plain": [ " Date Year Month Day T_mean (C) T_high (C) T_low (C) \\\n", "243 1938-09-01 1938 9 1 16.2 20.6 11.7 \n", "244 1938-09-02 1938 9 2 15.3 18.9 11.7 \n", "245 1938-09-03 1938 9 3 15.0 19.4 10.6 \n", "246 1938-09-04 1938 9 4 14.8 17.8 11.7 \n", "247 1938-09-05 1938 9 5 15.6 19.4 11.7 \n", "\n", " Rain (mm) Snow (cm) Total Precip (mm) season \n", "243 1.3 0.0 1.3 son \n", "244 2.3 0.0 2.3 son \n", "245 0.0 0.0 0.0 son \n", "246 3.3 0.0 3.3 son \n", "247 1.3 0.0 1.3 son " ] }, "execution_count": 14, "metadata": {}, "output_type": "execute_result" } ], "source": [ "season_dict[\"son\"].head()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Fitting the distributions\n", "\n", "What does it mean to \"fit\" a distribution? We want to be able to generate a series of random\n", "numbers that resemble temperatures and precip amounts that we might see at YVR during a \n", "particular season. We need to find two things:\n", "\n", "1. A probability density function (pdf) that has a similar shape to the histogram\n", " we get from the data.\n", " \n", "1. The best choice of parameters for that pdf that getting it \"closest\" in some statistical\n", " sense to the data.\n", " \n", "You'll see below that temperature is best fit with a Gaussin (normal) distribution, whilc\n", "precipitation more closely follows an exponential distribution. Both the normal and\n", "exponential distributions are functions of two parameters. For the normal distribution,\n", "the best estimate (maximum likelihood values) of those parameters is provided by the \n", "mean and the standard deviation of the temperature data. For the exponential distribution,\n", "they are given by the minimum value and the mean.\n", "\n", "### Daily average temperature\n", "\n", "The next cell shows histograms temperature for each of the seasons. We use\n", "[scipy.stats.norm.fit](https://docs.scipy.org/doc/scipy/reference/generated/scipy.stats.norm.html#scipy.stats.norm) to\n", "find the mean and standard deviation (called the `loc` and `scale` here) of the data\n", "and then shows that fitted distribution as a red line. Specifically, for temperatures x the red\n", "line is a plot of:\n", "\n", "\\begin{align*}\n", "f(y) &= \\frac{\\exp(-y^2/2)}{\\sqrt{2\\pi}} \\\\\n", "y &= (x - loc) / scale\n", "\\end{align*}\n", "\n", "The four plots show that a normal distribution give as pretty good representation of each of the seasons." ] }, { "cell_type": "code", "execution_count": 15, "metadata": { "lines_to_next_cell": 2 }, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAfsAAAHxCAYAAABqEBW0AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAIABJREFUeJzs3Xl8VNXZwPHfk4QEwr6GJUAi+07YVdC4gLhUFDesraK2qK++ttUuWq1ba6tvtdpa24p1KWoV0YooKKASCJRg2Pd9kRAUZA9LQpLn/eNOMIQsk2Tu3Fme7+czn5ncuXfmmZu558w959zniKpijDHGmMgV43UAxhhjjHGXVfbGGGNMhLPK3hhjjIlwVtkbY4wxEc4qe2OMMSbCWWVvjDHGRDir7I0rROR1EfmdiIwQkQ2llncTkWUickRE7vUyRmOMiRZxXgdgIpuqZgLdSi36JZChqmkehWSMMVHHzuxNsHUE1ngdhDHGRBOr7E1AiEiaiCz1Nc9PBur6lqeLSI7v8RfABcBfRSRPRLp6GLIxUUlEtovIL0RkpYgcFZFXRCRJRD7xHb+fiUhT37pTRORrETkkIvNEpFep13ldRP7m2y5PRBaISGsReV5EDojIehGxFrwQYZW9qTURiQemAm8AzYApwDVl11PVC4FM4B5VbaCqG4MaqDGmxDXASKAr8D3gE+DXQAuceqFkPM0nQBegFbAUeKvM61wPPOzbLh9Y6FuvBfAe8Cc3P4Txn1X2JhCGAXWA51X1pKq+B2R7HJMxpmIvqOo3qroL5wf4IlVdpqr5wAdAGoCqvqqqR3zLHwP6iUjjUq/zgaouUdUTvu1OqOokVS0CJpe8jvGeVfYmENoCu/T0WZV2eBWMMaZK35R6fLycvxuISKyIPCUiW0TkMLDd93yL6rxO4EI2tWGVvQmE3UA7EZFSyzp4FYwxJiC+D4wBLgYaAym+5VLRBiZ0WWVvAmEhUAjcKyJxIjIWGOJxTMaY2mmI0w+/D0gEfu9tOKY2rLI3taaqBcBYYDxwALgB+I+XMRljam0STnfcLmAtkOVtOKY25PRuVmMCS0QuBP6pqmd5HYsxxkQrO7M3busNbPM6CGOMiWaWLte4RkT+DFwJ3OJ1LMYYE82sGd8YY4yJcNaMb4wxxkQ4q+yNMcaYCBcxffYtWrTQlJQUr8Oo0NGjR6lfv77XYYQE2xff8WJfLFmy5FtVbRnUN62mUD+eK2Pfb4fth+DsA3+P54ip7FNSUli8eLHXYVQoIyOD9PR0r8MICbYvvuPFvhCRkE9lHOrHc2Xs++2w/RCcfeDv8WzN+MYYY0yEs8reGGOMiXBW2RtjjDERzip7Y4wxJsJZZW+MMcZEOKvsjTHGmAhnlb0xxhgT4SLmOntj/JXywPQq19n+1OVBiMSY8OHPcQN27IQqO7M3xlRJREaLyAYR2SwiD5Tz/HkislRECkXk2lLL+4vIQhFZIyIrReSG4EZujAGXK3srIIwJfyISC7wIXAr0BG4UkZ5lVvsKGA/8u8zyY8DNqtoLGA08LyJN3I3YGFOWa834pQqIkUAOkC0i01R1banVSgqIn5fZvKSA2CQibYElIjJTVQ+6Fa8xpkJDgM2quhVARN4BxgCnjmVV3e57rrj0hqq6sdTjXBHZA7QE7Fg2Jojc7LO3AsKYyNAO2Fnq7xxgaHVfRESGAPHAlgqenwBMAEhKSiIjI6PagYaCvLy8sI29Mvf3KfRrvZLPHqn7oTpCaR+4WdkHpYAwpjR/BxGZapFylmm1XkCkDfAGcIuqFpe3jqpOBCYCDBo0SMN1EpVInQBmvL8D9G5KByJ3P1RHKO0DNyt71wuIcDoTCKVfeF5zc1/4e/ZRlWD9r8Lke5EDtC/1dzKQ6+/GItIImA48rKpZAY7NGOMHNyt71wuIcDoTCKVfeF5zc1/4e/ZRlZKzE7eFyfciG+giIqnALmAc8H1/NhSReOADYJKqTnEvRGNMZdwcjX+qgPAd8OOAaf5saAWEMaFDVQuBe4CZwDrgXVVdIyJPiMiVACIyWERygOuAl0RkjW/z64HzgPEistx36+/BxzAmqrl2Zq+qhSJSUkDEAq+WFBDAYlWdJiKDcSr1psD3RORx3yU6JQVEcxEZ73vJ8aq63K14jTEVU9UZwIwyyx4p9Tgbp/Wu7HZvAm+6HqAxplKuZtCzAsIYY4zxnmXQM8YYYyKcVfbGGGOqrX7+MernH/M6DOMnq+yNMcZUjypvTv4N8176EcO3LfM6GuMHq+yNMcZUy9lfrSJt9wZitZhJ7z7CvQveRsrPlWRChFX2xhhjquXHX/6HvYlNSJ8wkam90rlv/lu8PuUxmh475HVopgJW2RtjjPFbl707uHDrYiYNuJyD9Rpx3+X38eAl9zDsq5VMf/0nkGVJEkORVfbGGGP89uPsDzgel8AbAy53Fojwdv/RXPODZyiKiYXzzoMXXgCtVnZ04zKr7I0xxvil1ZF9XLUmg3f7XszBeo1Oe251685cPv7PMHo03HsvPZ94Ao4c8ShSU5ZV9iYyqXJ79lT+OP15YoqLvI7GmIgwfulHxGoxrwy6qtznD9dtAFOnwtNP03LePBg8GFavDnKUpjyuZtAzxgsxxUU89tlEbl7mTIqzpF133uk/2uOojAlv9fOPcdOyT/i069l81bRNxSvGxMAvf8ny+HjSnnoK0tNh1y5ISAharOZMdmZvIkrCyXz+PvUP3LxsOv8YMpZFyb345bxJND5uzYnG1MYNK2fTOP8oLw8Z69f6h/r3h1degX37YO5cl6MzVbHK3kSOffv49zsPMXLTIh67aAJPXXAbj468k8Yn8vh55hteR2dM+Cos5LbFU1mU3Ivlbbv5v92FF0K9evDxx+7FZvxilb2JDNu2wbnn0vubLfzPVQ/w+qArAVjfKpVJAy7npmWf0OvrzR4HaUyYmjKF5MN7/T6rP6VePbjoIqeyt9H5nrLK3oS/pUvh7LNhzx5uGvc7Pu127mlPPzf8JvYnNuKJ2f+wLF/GVJcqPPMMW5ol83nnwdXf/oornB/j69cHPjbjN6vsTXibORPOP98Z/LNgAYuTe52xyuG6DXj6/PEMzF3P1WvmeBCkMWEsIwOWLuXlwVehUoMq47LLnHtryveUVfYmfL3+unPW0KkTLFwIPXpUuOp7fS5iWZtuPJjxGg3zjwYvRmPC3R//CK1a8UHvC2u2ffv20K8fTJ8e2LhMtVhlb8LTX/4Ct97qnNXPmwdt21a6ukoMj4y8k+ZHD/GT+f8OUpDGhLnVq+GTT+B//5f8uPiav87ll8P8+XDgQOBiM9Vilb0JPwUF8OijMGoUzJgBjRpVvQ2wqk0X3ul3CeOXfESXvTtcDtKYCPCnP0FiItx1V+1e54oroKgIZs0KTFym2qyyN+EnIwMOHoS774b46p1t/PG8H5KXkMjjn71ko4ONqUxuLrz5Jtx2GzRvXrvXGjIEWrSwfnsPWQY9ExZSHviuv+/3n/6VK+PrMXBeEfn/rV4/4IHExjxz3s38btbfuHz9fKb3GBHoUI2JDC+84JyN/+xntX+t2Fi49FKnJa6oyPnbBJWd2ZuwElNcxMhNWcw5a1CN+xD/3e8S1rQ6i4fmvEJiwfEARxiZRGS0iGwQkc0i8kA5z58nIktFpFBEri3z3C0issl3uyV4UZsaO3IE/v53uOYaOOuswLzmFVc42fQWLQrM65lqcbWytwLCBNqgXetoeewgn3Y9p8avURwTy29G3kXbI99y98J3AxhdZBKRWOBF4FKgJ3CjiPQss9pXwHjg32W2bQY8CgwFhgCPikhTt2M2tfTKK3DoENx/f+Bec9QoiIuzpnyPuFbZWwFh3HDphgXkx9ZhTqdBtXqdpck9eK/3Rfz4yw9I3b8rQNFFrCHAZlXdqqoFwDvAmNIrqOp2VV0JlM1adAkwW1X3q+oBYDZgsxKFMlX4859hxAgYOjRwr9ukCQwfbpW9R9w8s7cCwgSWKpdsXMi81AEci69X65d7+vzxnIiL59HPJtpgvcq1A3aW+jvHt8ztbY0XNm2C7dvhhz8M/GtfcQWsWgVffRX41zaVcnOAXnkHub8/E62AMGfot3sjbY98yzPnBaYQ2tugKX85dxwPz3mV3t9sYXXrzgF53Qgk5Szz99eR39uKyARgAkBSUhIZGRl+vkVoycvLC9vYAdpMn0434MuEBI6V+hz39yn0a/uSz17efqjXsiVDgY3PPUfumDFnbBtpQum74GZl73oBEU6FQyj9071Wk31xf59Czlkzn6LYWLpfOZD76/tX8FSlXofzYM6rPHQsm8V9Uk4tD9b/Kky+FzlA+1J/JwO51dg2vcy2GeWtqKoTgYkAgwYN0vT09PJWC3kZGRmEW+ylr3Z59tOlNE9szPVrmsHa0tkm/awuVjnb3N+niGfnl8lWqc3IaNKGrhs30jXM9lFNhNJ3wc3K3vUCIpwKh1D6p3utJvti/K8+5ooFWSxo35ffb20SwGhaMDSpEzELl/PsWTeeWrr9purFV1Nh8r3IBrqISCqwCxgHfN/PbWcCvy815mYU8GDgQzSBMnjnGhYn9wQp75yrlkT4otNgbvtiFhw75iTsMUHhZp/9qQJCROJxCohpfm47ExglIk19hcQo3zITpbrv3U7qgd18UmZGu0CYe9ZABu5aZznzK6CqhcA9OMfgOuBdVV0jIk+IyJUAIjJYRHKA64CXRGSNb9v9wG9xyoNs4AnfMhOCko58S4dD35BdzoRSgfJFp8Fw4gR88YVr72HO5FplbwWECaRLN/yXYoRZXYYF/LXnpQ4gTos5Z8eKgL92pFDVGaraVVU7qeqTvmWPqOo03+NsVU1W1fqq2lxVe5Xa9lVV7ey7vebVZzBVG7JzDQBfuljZf9m+NzRoYBPjBJmr19lbAWEC5ZKN/yW7fS/21Q9kE75jadvuHImvx/lblwb8tY0JJ4Nz1pIXX4+1SQFKpFOOgrg6zjX3H39sV8EEkWXQM6Fv40a6f7uDT2qRSKcyhbFx/LdjP87bttQKHxPVBuesYWnb7hTFuJzO9vLLIScHVq50933MKVbZm9D3n/8AMLPr2a69xbzUASQf3kOn/TmuvYcxoazRiTy67d1BdnLZ3GcuuOwy596a8oPGKnsT+t5/n+VturK7UUvX3mJu6gAA5+zemCg0KGctMSjZ7d3rrz+ldWsYPNiy6QWRVfYmtO3YAYsX80k3d5rwS+Q0ac2WZu2s395ErSE5ayiIiWNZm27BecPLL4esLNi7NzjvF+WqrOxFpHcwAjGmXB98AFCriW/8NS91AEN3ribhZL7r72VMqBm8cw2rWncmv05CcN7wiiucMTKffhqc94ty/pzZ/0NEvhSR/xGRwA+FNqYy778Pffuyo2lb199qbuoA6hXmMyRnjevvZUwoSTiZT5+vNwenCb9EWprTnG9N+UFRZWWvqsOBm3Cy4S0WkX+LyEjXIzPm669hwQJnTu0gWNS+D/mxdazf3kSdtN0biC8uZFH7IDbkxsQ4TfmffgonTwbvfaOUX332qroJeBj4FXA+8BcRWS8iY90MzkS5qVOdZr6xwfmaHY+vy5fJvazf3kSdwTvXUIywpF2P4L7xFVfA4cPOj3rjKn/67PuKyHM4WfAuBL6nqj18j59zOT4Tzd5/H7p2hV7Ba1qce9YAuu77CnburHplYyLE4Jy1bGjZkcN1GwT3jS++GOLjrSk/CPw5s/8rsBTop6p3q+pSAFXNxTnbNybw9u2DOXOcs3o3JuSowDzfJXjMtKkYTJQoLGTgrnV8Gcz++hINGkB6ul1vHwT+VPaXAf9W1eMAIhIjIokAqvqGm8GZKPbRR1BUFLT++hIbW3Rkd4PmVtmb6LFsGfVPnnB18ptKXXEFrF8PW7d68/5Rwp/K/jOgXqm/E33LjHHP++9Dhw4wcGBw31fEObufPRsKC4P73sZ4ITMTcHfym0pdeKFzP3++N+8fJfyp7Ouqal7JH77HNgmxcc/hwzBrVtCb8EvMPWsgHDoEX34Z9Pc2JugyM9nRpDV7Gjb35v27d4eGDWHRIm/eP0rE+bHOUREZUNJXLyIDgePuhmWi2owZUFAQ9Cb8EvNT+lMkMfz1Fy/w3IgDFa63/anLgxiVMS5QhfnzyU7u510MsbEwZIiTTc+4xp8z+58CU0QkU0Qygck489Qb447333eSbZzjfta88hyu24Dlbbpyvl1vbyLd+vXw7bfeNeGXGDYMVqyAY8e8jSOC+ZNUJxvoDtwF/A/QQ1WXuB2YiVL5+c6Z/VVXOUk3PDIvdQB9d2+iyfHDnsVgjOt8/fVBzZxXnmHDnAG5S6xqcYu/pelgoC+QBtwoIje7F5KJaosXO7/uL7nE0zDmnjWQGJQR25Z5GocxrsrMhFat2BaEdNSVGjrUubemfNdU2WcvIm8AnYDlQJFvsQKTXIzLRKt585z74cM9DWNl684cqNuQ87Yt46Oe53saizG1kfJAxdewz582i5WtO3syEPY0LVtCp042SM9F/gzQGwT0VFV1OxhjyMyEnj2hRQtPwyiOiWV+Sn/O277UGcTkdWFoTIC1ObyX5MN7eGXwGK9DcQwdCnPneh1FxPKnGX810NrtQIyhqMi51va887yOBHD67ZPy9tN973avQzEm4Ab7Znf0fHBeiWHDYNcuyMnxOpKI5E9l3wJYKyIzRWRayc3twEwUWrECjhyBESO8jgRwprwFOH+bDRoykWfIzjUcia/HulapXofiGDbMubd+e1f4U9k/BlwF/B54ttTNmMDyjQwOlcp+T8PmrGuZYlPe+ojIaBHZICKbReSBcp5PEJHJvucXiUiKb3kdEfmXiKwSkXUi8mCwYzdnGpyzhqXtelAcE+t1KI5+/SAhwSp7l1TZZ6+qc0WkI9BFVT/z5cX369shIqOBP/vW/6eqPlXm+QScgX4DgX3ADaq6XUTqAP8EBvhinKSqf6jG5zLhaN48SE2F9u29juSUuakDuG3xNBILjnMsvl7VG0QoEYkFXgRGAjlAtohMU9W1pVa7HTigqp1FZBzwNHADcB2QoKp9fOXHWhF5W1W3B/dTmBJNjh+m27dfMa2HN4NPKxo0+F6Lsyie/AnXx11gSasCzJ8pbn8MvAe85FvUDpjqx3YlhcOlQE+cS/Z6llntVOGAM13u077lpwoHnB8Cd5ScJZgIpeqc2YdIf32JeakDiC8uZNhXq7wOxWtDgM2qulVVC4B3gLIju8YA//I9fg+4SEQE5+qd+iIShzPPRgFgCQw8NChnHRAC19eXsbxNV/p+vZm4IpuXItD8aca/GzgX38GpqpuAVn5sZ4WD8d+GDbB3b8g04ZdYnNyLY3USrN/e+ZG/s9TfOb5l5a6jqoXAIaA5zrF9FNgNfAU8o6r73Q7YVGzIztXkx8axok1Xr0M5zbK23albWGCDYl3gz6V3+apaIL5Lj3wVsD+X4ZVXOAytaB1VLRSR0oXDGJzCIRH4mRUOEa7k+voQO7MviKvDwg59rd8eyrv2sGw5UNE6Q3BydLQFmgKZIvKZqp42p6mITAAmACQlJZGRkVHbmD2Rl5cXUrHf3+fMs+Qr31vNvk5duCctBnDnLDqpXvnvXZkGbTrDNPg5a0NqH9ZUKH0X/Kns54rIr4F6IjISJ2XuR35sZ4VDKaH0T/daefuix5QpNGnWjIU5Oc7lN2VUt9AIpLjd/UidlM3jrXZyOKnNqeWB+H+G0fciByg9mCIZyK1gnRzfSUFjYD/wfeBTVT0J7BGRBTj5O047nlV1IjARYNCgQZqenu7Cx3BfRkYGoRT7+DL94/UKTnDntq28PORqnl3lTxVQM/f3Kaz+62trLmvQjP2LNzH29XRX4gqmUPou+POfeACnb30VcAcwA2fwXFWscCgllP7pXit3X9x8M1x8MekXXFDuNmULrGBKTRjMHP7JxpmreCvtu6/09pvSa/3aYfS9yAa6iEgqsAsYh3OcljYNuAVYCFwLfKGqKiJfAReKyJs4LXXDgOeDFrk5TVrueuoUF4XO9fWlibCsbTfSctd7HUnE8WcinGJVfVlVr1PVa32P/WnGP1U4iEg8TuFQ9vr8ksIBShUOOP16F4qjPk7hYP/9SLVjB+zcGXL99SW2NW3LnvpNTyUhiUa+Pvh7gJnAOuBdVV0jIk+IyJW+1V4BmovIZuA+nBMFcAbqNsBJ0JUNvKaqK4P6AcwpQ3LWUIywtF0Pr0Mp17K23Ug9sBv27fM6lIjiT278bZTTR6+qZ1W2na8PvqRwiAVeLSkcgMWqOg2ncHjDVzjsx/lBAE7h8BpO4SBY4RDZQrS//hQRvkzuxeCda6teN4Kp6gyclr3Syx4p9fgEzpU0ZbfLK2+58cbgnDWsb5XC4boNvA6lXMvbdHMeLFoEl13mbTARxN/c+CXq4hy0zfx5cSscjF/mzYMmTaB3b68jqVB2+15csWE+7Q7tYVdjfy5GMSb0xBUVkpa7gXf7jPQ6lAqtbN2FIokhNivLKvsA8qcZf1+p2y5VfR64MAixmWiRmenMcufh/PVVyfb1b0ZzU74Jf72+2ULiyfxT3+dQdDy+LutbplgmvQDzJ6nOgFK3QSJyJ9AwCLGZaPDNN8419qHahO+zvmVHDscnMmSnVfYmfKXlbgBgSYj215dY1rab04xfXOx1KBHDn2b80nnwC4HtwPWuRGOiT0k+/BCv7ItjYlmS3MPO7E1YS8vdwO4Gzfm6kbdTSFdlWdvu/GD5J86JQI/Q/mESLvzJjV/+tVDGBMK8eZCYCAMGeB1JlbKTe3HB1iU0PXaIA4mNvQ7HmGrrv3uDc9Yc4pa39WX2y8qyyj5A/BmNf19lz6vqnwIXjok6mZlw9tlQp47XkVTpy/Yl/fZrmdX1bI+jMaZ6mh89SMeDX/Nm/9Af9La1WTtn0G5WFtx6q9fhRAR/RkQNAu7CSW3bDrgTZ2KbhljfvamNgwedOexDvAm/xMrWXcmPrWNN+SYs9d/t9Ncvaxf6Z/YqMTB0qA3SCyB/+uxbAANU9QiAiDwGTFHVH7kZmIkCCxY4s92FSWVfEFeH5W26WmVvwlJa7gYKJYbVSZ28DsU/w4bBb38LeXnQIDRzAoQTf87sO+DMOleiAEhxJRoTXTIzneb7oWXnRwpd2e170fvrLSQWHPc6FGOqJS13PetapXKiTl2vQ/HPsGHOaPzFi72OJCL4U9m/AXwpIo+JyKPAImCSu2GZqDBvHgweDPXqeR2J37KTexGnxacuYTImHMQUF9Fv9yaWte3udSj+GzLEubem/IDwJ6nOk8CtwAHgIHCrqv7e7cBMhDt2DLKzw6YJv8SSdj0okhi73t6ElS77dtKg4HhYjMQ/pVkz6NrVKvsA8TdlWSJwWFX/jDNDXaqLMZlokJUFhYVhV9nnJSSyrlWq9dubsJK2y5lHLKwqe3Ca8rOynLE9plb8yaD3KPAr4EHfojrAm24GZaJAZqaTHvecc7yOpNqyk3s6zfgnT3odijF+ScvdwIG6DdnetK3XoVTPsGFOls0dO7yOJOz5Mxr/aiANWAqgqrkiYpfcmdqZNw/69YPG4Zec5svkXty65COuuv0FlldxprT9qcuDFJUxFUvL9SXTEfE6lOoZNsy5z8qClBRPQwl3/jTjF/jmmFcA3/zyxtSYnDwJCxeGXRN+ieyS5DrWb2/CQMP8o3TetzP8mvAB+vRxBvAuWuR1JGHPnzP7d0XkJaCJiPwYuA142d2wTCRruGkTHD9+qrJPeWC6xxFVz7f1m7K1aVuG5Kzh5aFjvQ7HmEr13b2JGDS8RuKXiIuDQYNskF4A+DMa/xngPeB9oBvwiKq+4HZgJnI1XrHCeTB8uLeB1EJ2ci8G5axF1GblMqEtLXc9xQgrSvLNh5thw2DpUsjP9zqSsFZpZS8isSLymarOVtVfqOrPVXV2sIIzkanJypXQvTu0auV1KDWW3b4XTU8cofO3O70OxZhKpeVuYHPz9hxJCNMe2GHDoKAAli/3OpKwVmkzvqoWicgxEWmsqoeCFZSJYEVFNF61Cm66yetIaiU7uScAQ3LWsKllR4+jMaYCqqTlbmB25/DJUlmipHuv1ZE8vgQef+gVXhv07Wnr2ABY//nTZ38CWCUis4GjJQtV9V7XojKRa/Vq4o4eDdvBeSV2NGnDnvpNGZyzhrfSQn8WMROltmyh2fHD4Tk4z2dPw+bsatiS/rkbvQ4lrPlT2U/33YypvXnznPsRI7yNo7ZE+DK5F4N3rvU6EmMq5hvYtqxdGA7OK2VZ226k5a73OoywVmGfvYh87nvYU1X/VfYWpPhMpJk3jxNJSdChg9eR1Fp2+160O7KXdof2eB2KMeXLyiIvvh6bmrf3OpJaWd62Kx0OfUOLowe8DiVsVTZAr42InA9cKSJpIjKg9M2fFxeR0SKyQUQ2i8gD5TyfICKTfc8vEpGUUs/1FZGFIrJGRFaJSJhM1WQqpAqZmRzs29frSAIiO9l3vX0UpM61YzlMZWWxsnUXimNivY6kVkouG7QJqGqussr+EeABIBn4E/BsqdszVb2wiMQCLwKXAj2BG0WkZ5nVbgcOqGpn4Dngad+2cTgpee9U1V5AOmC5ScPdpk3wzTcc6tfP60gCYn3LjhyOT4z4SXHsWA5Tx47BihVh3V9fYnVSJwpi4hiwy5rya6rCyl5V31PVS4H/U9ULytwu9OO1hwCbVXWrqhYA7wBjyqwzBijpEngPuEhEBBgFrFTVFb5Y9qlqUTU/mwk1c+YARMyZfXFMLEuSe0TDmb0dy+Fo6VIoLAzPZDpl5NdJYFXrztFwrLnGn6Q6v63ha7cDSl+EnONbVu46qloIHAKaA10BFZGZIrJURH5ZwxhMKJkzB9q25XhysteRBEx2ci+67NtJ02MRfWWqHcvhyDc4b3m4JtMpI7t9L/ru3kTCSUuuUxP+jMavqfJmXCg7T2FF68QBw4HBwDHgcxFZoqqfn7axyARgAkBSUhIZGRm1jdk1eXl5IR2f61Q5Z9Ys9g8aRN7Ro6fti/v7FHoXVy21ie8O8+B3savY2mfYGc9X9T8Pk++F68cyhNfxXJlQ+Z/2/PhjGrZpwy3DGgLBP8aS6gX22G53shvxiwrhlTCVAAAgAElEQVT5feJ6crs742VCYT9XJlS+C+BuZZ8DlB4CmgzkVrBOjq9vrzGw37d8rqp+CyAiM4ABwGkFhKpOBCYCDBo0SNPT0wP/KQIkIyODUI7PdWvXwoEDtB43jgYNGpy2L8aHWW780uILe/C92DrkLljPswlnpv/dflN6pduHyffC9WMZwut4rkzI/E+3bIH0dJ5d5WYxX7H7+xQG9L2bFPbme8CWeRv420ln3E9Vx5fXQua7gH+z3iEiw0XkVt/jliKS6sdm2UAXEUkVkXhgHDCtzDrTgFt8j68FvvDNsDcT6Csiib6C43zALmgOZ77+ei64wNs4Aqwgrg7L23SN9L5EO5bDTU6Ocxt2ZmtTuDpYrxEbWnSI9GPNNVVW9iLyKPAr4EHfojo4o2sr5eu3uwfnYF8HvKuqa0TkCRG50rfaK0BzEdkM3Icz+h9VPYBzBUA2sBxYqqrhe/pnnMq+Y0dI9ed3YnjJbt+L3l9voV7BCa9DcYUdy2GoZErYCKrsARYn92TgrvXEFNsYz+ryp43laiANWAqgqrki0tCfF1fVGcCMMsseKfX4BHBdBdu+iR8/KkwYKC52Kvsrr6x63TCUndyLe/Rd0nLX89+U/l6H4wo7lsNMVhYkJED//vCfyJm7LDu5Fzct/5Ru3+5gXauzvA4nrPjTjF/ga45TABEJ06mTjGdWrYL9+yOuCb/EknY9KJIYhljzogkVWVkwYADEx3sdSUCVJLIalGM9QdXlT2X/roi8BDQRkR8DnwEvuxuWiSgR2l9fIi8hkXWtUq0v0YSGkydh8eKIa8IH2NWoJbkNWzDYKvtqq7IZX1WfEZGRwGGgG/CIzWlvqmXOHOjcGdqHd37uymQn9+SGlbOIKyqkMNab0c/GALByJZw4EZGVPSIsTu7J4J1rnPTbxm/+DND7GbBOVX+hqj+3it5US1ERzJ0bsWf1Jb5M7kXiyXx6f7PF61BMtPMl04nIyh7nWGuTt4/kwzYBVXX404zfCJgpIpkicreIJLkdlIkgy5bBoUMRX9lnt/dNihPhefJNGMjKgjZtIrYlbXGyMy2D9dtXjz/pch/3TWBxN9AWmCsin7kemYkMJf31IZJYwi3f1m/KlmbtOOerFV6HYqJdVhYMHQpSXlLD8LexRQcOJ9SP+AmoAs2vpDo+e4CvgX1AK3fCMRFnzhzo3t0504hwmSlpDN25mvhCm9TNeGTfPti8OWKb8MGZgGpxux52Zl9N/vTZ3yUiGTjpLVsAP1bVyJi2zLjr5EnIzIQL/ZkkMfxlpqaReDKfgbvWeR2KiVYRmkynrMXJPem67yvnx43xiz9n9h2Bn6pqL1V9VFXt55Txz5IlkJcX8f31JbLa9+FkTCwjti/1OhQTrbKyICYGBg3yOhJXZfv67VmwwNtAwkiFlb2INPI9/D/gKxFpVvoWnPBMWIuS/voSRxMSWdq2OyO2LfM6FBOtsrKgb1+oH9m5z1a26Up+bBzMn+91KGGjsjP7f/vulwCLffdLSv1tTOW++AL69IEWLbyOJGgyU9Po9c1WmkX2/PYmFBUXO834Ed6ED5AfF8/K1l2dbkLjlwore1W9wnefqqpn+e5LbpaU2FQuP99pYouSJvwSmSlpxKCcu32516GYaLN+PRw+HBWVPfguwVuyBI4d8zqUsODPAL0z5p0ub5kxp/nySzh+POoq+1WtO3OwbgNGbLemfBNk//2vcz90qLdxBMmX7Xs5g4Czs70OJSxUmNdTROoCiUALEWkKlFy02QjnentjKjZnjnOd7/nnex1JUBXHxDK/Y3+n397SeRqXpTzw3WzBf5n2BsPqN2XIa5tBIj+T45J2PZwHmZlRV87URGVn9nfg9M935/T++g+BF90PzYS1OXMgLQ2aNvU6kqDLTE2jTd4+Ou/b6XUoJkqIFnPu9uVkpvSP2GQ6ZR2u2wB697ZBen6qrM/+z6qaCvy8TJ99P1X9axBjNOHm+HGnSTHKmvBLzE9JA+A8G5VvgqTnN1tpfvwwmakDvA4luEaMcMqaoiKvIwl5/qTLfUFEeovI9SJyc8ktGMGZMLVwIRQURG1lv6txK7Y0S7br7U3QjPANCF3Qsb/HkQTZ8OFw5Igz05+pVJVzcYrIo0A60BOYAVwKzAcmuRqZCV9z5kBsLIwYcVqfYon7+xQyvpzlkWReahrjVsxyrkpISPA6HBPhRmxfyrqWKextEGXdZsOHO/fz5zvdhqZC/mTQuxa4CPhaVW8F+gFWepmKzZkDAwdCo0ZVrxuhMlPSqFeYbxm+jOvqnjzBoJy1ZKZEYWXXoYNzs377KvlT2R9X1WKg0JdVbw9g19mb8h096lx2FyX58CuS1aEPBTFxMGuW16GYCDd05xoSigqZnxJlTfglhg93RuTb1S+V8qeyXywiTYCXcUbjLwW+dDUqE74WLHCufY3S/voSx+LrsbRdd6vsjeuGb19Gfmwd57rzaDR8OOzeDdu2eR1JSPNngN7/qOpBVf0HMBK4xdecb8yZvvgC6tSBc8/1OhLPzUsdAMuWwZ49XocSECIyWkQ2iMhmEXmgnOcTRGSy7/lFIpJS5vkOIpInIj8PVszRYMS2ZWQn9+REnbpeh+KNESOce2vKr1RlE+EMKHsDmgFxvsdVssIhCs2ZA0OGRPxEHP441Yf62WfeBhIAIhKLk1/jUpzBujeKSM8yq90OHFDVzsBzwNNlnn8O+MTtWKNJy7z9dP92x6nLPaNSz57QpInlya9CZaPxn63kOQUq7ZQtVTiMBHKAbBGZVmaK3FOFg4iMwykcbij1vBUO4eTwYSdX9YMPeh1JSFiTdBY0a+Y05X//+16HU1tDgM2quhVARN4BxgClj+cxwGO+x+8BfxURUVUVkauArcDR4IUc+Yb7LrnLTI3iyj4mxmlJtDP7SlVY2atqbTtdrXCINpmZTnKLKO+vL1EcEwsXX+xU9qrhntmsHVA6JWAOUDYJ+6l1VLVQRA4BzUXkOPArnB/+1koXQCO2L+PbxMasbZXqdSjeGjECpk+HvXuhZUuvowlJ/lxnnwjcB3RQ1Qki0gXopqofV7GpFQ7RZs4c55rys8/2OpLQMWoUvPsurFnjpPYMX+X9Uik7/LmidR4HnlPVPKnkB4+ITAAmACQlJZGRkVGzSD2Wl5cXnNhVuSRnGTn9+nJf32Kg2P33rIakek5ODTeV7OdGiYkMAFa/9BLfllx7HwKC9l3wQ5WVPfAazij8c3x/5wBTgKoqeyscSgmlf7pbBk6bRmGPHqxYtOjUsvIO9mAUAqFiYVwDzgY2//3v5Fx33RnPh9H3IgdoX+rvZCC3gnVyRCQOaAzsx/mRf62I/B/QBCgWkRNl026r6kRgIsCgQYM0PT3djc/huoyMDIIS+8qVcOgg/2wygCmr/CnKg+v+PoU863Jc229Kdx6cfTb84hf0PngQQuh7E7Tvgh/8+U90UtUbRORGAFU9LpXVwN+xwqGUUPqnu2L/fti8GR577LTPWV6mvGAUAqFi+1M3wGOP0XnrVjqX8/8Po+9FNtBFRFKBXcA4oOxAhGnALcBCnGRcX6iqAiNKVhCRx4A8m18jAGbPBojOZDplJSQ4A4NtkF6F/LnOvkBE6uE7KxeRTkC+H9udKhxEJB6ncJhWZp2SwgFKFQ6qOkJVU1Q1BXge+L0VDiFu3jynX9r66880ahTMnQsnTngdSY2paiFwDzATWAe8q6prROQJEbnSt9orON1wm3G6/s64AscE0KxZbGrenq8btfA6ktAwfDgsXeok9jJn8KeyfxT4FGgvIm8BnwO/rGojKxyizEcfQcOGzq9rc7pRo5yZAMM8da6qzlDVrqraSVWf9C17RFWn+R6fUNXrVLWzqg4pGZxb5jUeU9Vngh17xDlxAubNi96seeUZPhwKC50MnuYMlbal+prr1wNjgWE4few/UdVv/XlxVZ2BM3lO6WWPlHp8AjizI/P09R/z572MhwoK4D//gauusklfynP++U6ioVmz4KKLvI7GRIL58+HECSdxk3Gcc45zxUtmprUwlqPSyt53CdxUVR0IRPY0ZabmZs2CgwfhhhuqXjcaNWjgXAc8axY8XTbPjDE1MHs21KnDovZhfYVHrZWdVfOjVmdR8NLbXHNs4GnLtz91eTDDCkn+NONnichg1yMx4WvyZGjaFEaO9DqS0DVqFCxfDt9843UkJhLMng3nnMOx+HpeRxJSPu16DgNz19P6sF+Nz1HFn8r+AmChiGwRkZUiskpEVrodmAkTx4/Dhx/C2LEQH+91NKFr1CjnPgJS5xqP7dnjzLlgP67PMKO7c439pRvDe3yMG/yp7C8FOuGkx/0ecIXv3hj45BM4csSa8KuSlgbNm9sseKb2Pv/cubfK/gzbmrVjXcsULltvlX1Z/sx6t6O8WzCCM2Fg8mQnPaUNiKlcTIxTOJekzjWmpmbNcrrNBg6set0oNKPbuQzetZakI9aUX5o/Z/bGlO/oUfj4Y7j2WoiLjiQ5tTJqFHz9Naxe7XUkJlypOv31F10EsbFeRxOSSpryR29c6HEkocUqe1NzH30Ex45ZE76/SppdrSnf1NS6dbBr13djQMwZtjRvz/oWHblsvc2CV5pV9qbmJk+Gtm2dZBamasnJ0KMHzJzpdSQmXPlS5Fp/feU+6XYug3PW0jJvv9ehhAyr7E3NHDoEM2bAdddZc2J1jBkDX3wBuWWniTDGD7NnQ+fOkJLidSQhbXr34cSgjN74X69DCRlW2Zua+fBDJ3OeNeFXz223QVER/OtfXkdiwk1BAWRkWBO+Hza36MDG5h243JryT7HK3tTM5MnQsSMMG+Z1JOGlSxcnfe4//wnFoTX/uAlxCxc6g2KtCd8vM7qfy5Cda2iZd8DrUEKCVfam+vbtcwaZXX+9k4vaVM+PfgRbtzpnacb4a9Ysp8vMLnP1y4xu5xKDcskmG5UPVtmbmvjgA2d2qXHjvI4kPF1zDTRp4pzdG+Ov2bNh6FBo3NjrSMLCxhYd2dws2Ubl+1hlb6pv8mRnkFBamteRhKd69eAHP4D33yfu0CGvozHhYN8+WLzY+uurQ4Tp3YczdOdqJ8VwlLNMKKZ6vvnGGU3+4IOkPDij6vWjXNlZuUr0ONqVTwoKKHr/Y1IWxtmsXKZyX3zhJNSx/vpqmdF9OD/57zvOFNx33ul1OJ6yM3tTPe+/7wwssyb8WlnX6ixWtO5CrzmWPtf4Yfp0aNQIhgzxOpKwsqFFR7Y0awfvved1KJ6zyt5Uz+TJ0LMn9I7uebQDYXK/UbTI+Yp+uzd6HYoJZfv3w7vvOpe5Wlrq6hFhRrfhMGcO7N3rdTSessre+G/XLsjMtGvrA2Raj/M5mZDAuBWWUc9U4pVXnKmk//d/vY4kLM3ofq7TGvnBB16H4in7mWj8N2WK0+RslX1A5CUksmnocK5cOM+ZJrhhQ69DMiEk5YHpxBQXMe+lZ8hp35txb30FfOV1WGFnXctUJ7/FlCkwYYLX4XjGzuyN/955B/r3h27dvI4kYqxJH0n9kyecZlpjyrh485ckH97D6wO/53Uo4UvESes9Zw58G73T3lplb/yzfTssWmRn9QH2dZdubGre3q65N+W6ZelH7GrYktldLFNlrVx7rZOmeupUryPxjKuVvYiMFpENIrJZRB4o5/kEEZnse36RiKT4lo8UkSUissp3f6GbcRo/lJx5WmUfWCK803cUZGWF9Dz3diwHX5e9Ozh3x0reSruUohibbKpW+veHTp2cpvwo5VplLyKxwIvApUBP4EYR6VlmtduBA6raGXgOeNq3/Fvge6raB7gFeMOtOI2f3nnHuewnNdXrSCLOB70vhDp1nIFYIciOZW/csvRj8mPr8Ha/S7wOJfyVNOV//rmToCgKuXlmPwTYrKpbVbUAeAcYU2adMUDJ9F/vAReJiKjqMlUtmQN0DVBXRBJcjNVUZtMmWLbMzupdsj+xMVx1FUyaBPn5XodTHjuWg+3gQcau+YJpPc7nQKKlxw2I665zmvI//NDrSDzhZmXfDthZ6u8c37Jy11HVQuAQ0LzMOtcAy1Q1JEvBqDB5snN//fXexhHJfvxj53rq0Lw8yI7lYHvtNRJP5vP6wCu8jiRypKU5LZNR2pTv5qV35U2HVjZVWKXriEgvnObAchNCi8gEYAJAUlISGSE8i1heXl5Ix1eR2GPHGPqnP3E0LY0VmzfD5s2nnru/T2GNXjOpXs23jTQl+yIjth7DkpI4/sc/sqJ1a6/DKsv1Y9m3Ttgcz5Wp9bFeXMzQZ57hQNcejB6ZwmjC81gJpeP8hbecs/lzevWn/yfTeOWlN8lvcPqlrn3aBb4FJZTKfTcr+xygfam/k4HcCtbJEZE4oDGwH0BEkoEPgJtVdUt5b6CqE4GJAIMGDdL09PRAxh9QGRkZhHJ8FXr4YThwgPhPPyW9TKrO8RXkfa/K/X0KeXaVpXiA7/bF9psuhLvvpu4jj5DeoQOcdZbXoZXm+rEM4XU8V6bWx/r06ZCbyy+u/AEfh/FxEorH+axm5/FR0Qds+nAxU/qePs/A9pvSA/5+oVTuu9mMnw10EZFUEYkHxgHTyqwzDWfQDsC1wBeqqiLSBJgOPKiqC1yM0VRm50549lm48UbLyR0M48dDTAy8+qrXkZRlx3Iw/eUv0LYtn3Y9x+tIIs6q1p3Z2TiJy6Nw2lvXKntfv909wExgHfCuqq4RkSdE5Erfaq8AzUVkM3AfUHJJzz1AZ+A3IrLcd2vlVqymAg895GTM+8MfvI4kOrRvD6NHw2uvQWFoNH+CHctBtWEDzJoFd95JYWxonRVHBBHe730h6duWkLZrvdfRBJWr19mr6gxV7aqqnVT1Sd+yR1R1mu/xCVW9TlU7q+oQVd3qW/47Va2vqv1L3WxC4mBasgTeeAN++lPo2NHraKLHj34Eubnw6adeR3IaO5aD5K9/hfj4qE7r6raJQ8byTYNmPPr5RESLvQ4naCyDnjmTKtx/P7RoAQ8+6HU00eWKKyApyTLqRaPDh+H1152rXpKSvI4mYh2Lr8dT54+n/+6NjF09x+twgsbaicwZfnztI7w8dy4Pj7yLN/8QfX1bnqpTx+m7f+YZJ7dBWprXEZlg+de/IC/PZrcLgqm90rl56XR+Nfd1Pu16NkcTEr0OyXV2Zm9Od/IkD2S8xuZmybzdf7TX0USnX/4SWraEW24J1SQ7JtCKi50m/CFDbDBsEKjE8PjFE2h19AB3Z0XHJFRW2ZvT/eMfdNq/i99fcJvl4/ZKs2ZOM/6qVfD4415HY4Jh9mzYuNHO6oNoedtuvN/7Qm7PnkqHA7u9Dsd1Vtmb7xw8CI8/zoKOffmi02Cvo4lul18Ot90GTz/tzDZoItsLL0CrVk5KVxM0T50/npOxdXhoTmjOSxFIVtmb7zz5JOzfz5MX/MiZOMJ467nnIDkZbr4Zjh3zOhrjls2bYcYMuPNOSLBpA4Jpb4NmvHj29VyyKQs++8zrcFxllb1xbNvmJPO45RbWJoVU9rbo1aiRk2Bn40Yn54GJTC++CLGxcMcdXkcSlV4dNIYdTVrDT34SUvktAs0qe+N48EGnwPnd77yOxJR20UVwzz3w/PMQIjm2TQAtXAh/+5szo2Tbtl5HE5Xy4+J58oLbYe1a+Mc/vA7HNVbZG8jKcma2+/nPoV3ZycyM5556Cjp3hltvhSNHvI7GBMrOnXD11U7mxD//2etootqsLsOcH9aPPBKx891bZR/tVOG++6B1a+eSLxN66td3rsH+6ivnB5kJf8eOwVVXOffTpkHzsrMBm6AScVrPDh2CRx/1OhpXWGUf7d57z2lK/O1voUEDr6MxFTnnHCer4cSJMHOm19GY2lB1rrRYtgzefht69vQ6IgPQuzfcdRf8/e/OZa8Rxir7aLZxI/zsZ9Cnj9NEbELbE084FcPtt8OBA15HY2rqySedbrOnnnIusTSh44knoEkTZ04QVa+jCShLlxutFi+GSy91mq8mTXIG55nQVreu878aOtQZOTxpktcRmer64AP4zW94v9cF3L+vJzww3euITGnNmjkV/j33wNSpzpiKCGGVfZRJeWA652xfzsQPnuRg3Yb88Ibfsu2dXfDOLq9Di2opfhT625+6HAYOhIcfdjLrXX11RBVGEW/lSvjhD1nWphu/Hv2/lssiVN1xh9OUf999cPbZznimCGCVfQTxp8K4bP18nvv4GbY1bcfN1z/BnoY2MCjsPPSQM6jrjjtg+HAnj74JbXv2wJVXQpMmTLjqIfLj4r2OyJRRuvwc1P9mJr37CIe69OaOsQ+xsk1XwPeDO0xZn30U+cGyGfz1w6dZ2bor19/0tFX04apOHacJXxVWrPA6GlOVggK49lr45huYOpW9DZp5HZGpwuLkXlzzgz9SGBvHlLd+xdjVn3sdUq1ZZR8NVLl3wdv8btbfmNNpED+84QkO17WR92Gtd2/YsQMuvtjrSExlVOHuuyEz05mrftAgryMyflrX6iyuvPlPLGnXgz9Nf47ffP5yWGfYs2b8MOFPE315RIt59LOJjF/6Me/3vpBfjb6Xwlj7t0eExMifgzvsvfCCM4PhQw85WfJMWDmQ2Jibr3+Ch+a8wu2LP4TRo50rKcIwL4Kd2UewOkUn+fNHzzB+6cdMHHw1P7/sp1bRGxMMq1fD2LHOVRNjxjgjvE1YKoyN4/GL7+AXl/7EaaEZPDgsr8O3kj8SqdLn6838au7rDN+xgj+kj+elodd6HZUxES99wkReWvMWxf/NJC++Hv8cfhMTO13NiV9/4nVoppam9B3JHx8Z5/yIO/tsZ9zM2LFeh+U3q+wjSKMTeYxZm8GNK2bSc882jscl8ItL72VK31Feh2ZMZNuxA377Wz579TU0vg7/GHYNE4eM5WC9Rl5HZgJp2DAnR8nYsXDNNfCb38Cvf+3kwAhxVtmHO1UG7VrLjStmctn6BdQrzGd1UiceHvU/fNjzfI4k1Pc6QmMi1+7dTka8iRNBhEkDrkBvvZr/22mXQ0astm2dGSjvustJM/7cc06CsquucjIiNm7sdYTlcrWyF5HRwJ+BWOCfqvpUmecTgEnAQGAfcIOqbvc99yBwO1AE3KuqEZsQvCaD75oeO8TY1V8wbuUsuuzbyZH4evyn9wW83W80q1t3diFKE+3sePbZt8+55HHGDGcu+sJCJ9f9ww/zxIsrub9JIez0OkjjhtPK6lbXMPz6Tly2YQEjZ8ym5ZQpFMTEEX/xhU6yqzFjvAu0HK5V9iISC7wIjARygGwRmaaqa0utdjtwQFU7i8g44GngBhHpCYwDegFtgc9EpKuqFrkVrxtKfzHu71PI+GpW6qLFJB3ZT+qBXFJ8t9QDuaTsdx7HFxeytG03fnHpvUzvPoJj8fUC/RFMCPE7y54LovJ4LiqCTZucir30bZcv26QI/OAHzixpnTr5NlrpWbgmyESYn5rG/NQ0Hi6+i/65G7lk00Lu2LLCOeu/6y7SevaEm292Ml+2bevcGjf2JHuim2f2Q4DNqroVQETeAcYApQuHMcBjvsfvAX8VEfEtf0dV84FtIrLZ93oLAxrhiRNw/Pjpy8qb/KBkmeqpxwOfmIXgPBYFUGK0mLjiYmK0mNjiIjppMTHFxcRqMS22n6TvbogvOkn9ghPUO3niu/uTx0ksyD/1uPmxQ3Q8sJuUA7upV5h/Koz82DrsaNKG7c3aMrvLUD7seT4bW6YEdJcYU4HQP56//to5yy4qOvO+9OPjx+Hw4dNuL364lAYFx2iQf4yGBcdplbePbnu/OnX8nYyJZUuzZNa16sK6LqNY1yqVNUmd2J/YGF5eD6wP6Ecx4aU4JpalyT1YmtyDO/5wGaxZA1OnEjNpEjzwwOkrJyZ+V/GXvvXoAZdd5lqMblb27Ti9MSsHGFrROqpaKCKHgOa+5Vlltm0XiKBKnx3dmfUeD8x9vUavs6QG29xYxfNH69TlWJ26HKrbgG3N2jI/pT/bm7Y9dctt1AIVu1rSeCLkjueyLR0b/3gV8cU1S3pyh8RwJKE+eQmJ5MXXY19iY97qP5p1rc5iXatUNjdvT0FcndqGbKJAyoMzfI/SuP/RPrz534N0PPg1SUf20SpvP63z9vHjTnUhN9cZ7Ldrl/MD9LzzwrayL6+douxpc0Xr+LMtIjIBmOD7M09ENlQnwAd9tyBpAXxb6RonTzi3Ywdhf05wovLAvf7siygR6H0hT/u1WseavHQ5y0LqeE6ozspnRFMMJ444txI7apaK2L7fDtsP3+2DL8ssn5Bdzsrz5tW0ed+v49nNyj4HaF/q72Qgt4J1ckQkDmgM7PdzW1R1IjAxgDG7RkQWq6rlysT2RWlhtC/sePZTGP1PXWX7IbT2gZttwtlAFxFJFZF4nAE608qsMw24xff4WuALVVXf8nEikiAiqUAXoOyPI2NM8NjxbEwYc+3M3tdndw8wE+dSnVdVdY2IPAEsVtVpwCvAG74BO/txChB8672LM/inELg75EfuGhPB7Hg2JryJljf63ASciEzwNVNGPdsX37F9EXnsf+qw/RBa+8Aqe2OMMSbC2XVcxhhjTISzyt5FInKdiKwRkWIRGVTmuQdFZLOIbBCRS7yKMZhEZLTv824WkQeq3iKyiMirIrJHRFaXWtZMRGaLyCbffVMvYzS1E43fcfteg4i0F5E5IrLOV+b/xLc8ZPaDVfbuWg2MBeaVXlgmfeho4G++dKQRq1S61UuBnsCNvv0QTV7H+X+X9gDwuap2AT73/W3CUBR/x1/HvteFwP2q2gMYBtzt+9+HzH6wyt5FqrpOVctLDHIqfaiqbgNK0odGslPpVlW1AChJtxo1VHUezij10sYA//I9/hdwVVCDMoEUld9x+16Dqu5W1aW+x0eAdThZIkNmP1hl743yUo8GJB1wCIvGz+yPJFXdDU6BAbTyOB5Tc/Yd/07Ufq9FJAVIAxYRQtMOcqQAACAASURBVPvB5rOvJRH5DGhdzlMPqeqHFW1WzrJIvywiGj+ziS72HY9yItIAeB/4qaoeFg9mt6uIVfa1pKoX12Azv9KHRpho/Mz++EZE2qjqbhFpA+zxOiBTY/Yd/07Ufa9FpA5ORf+Wqv7Htzhk9oM143sjGtOH+pNuNRqVTjF7C1BRa5AJffYd/05Ufa99Uzm/AqxT1T+Veipk9oMl1XGRiFwNvAC0BA4Cy1X1Et9zDwG34Yzi/KmqfuJZoEEiIpcBz/NdutUnPQ4pqETkbSAdZyasb4BHganAu0AH4CvgOlUtO9jJhIlo/I7b9xpEZDiQCawCin2Lf43Tbx8S+8Eqe2OMMSbCWTO+McYYE+GssjfGGGMinFX2xhhjTISzyt4YY4yJcFbZG2OMMRHOKntjjDEmwlllb4wxxkQ4q+xNwPjmcU4XkV+LyD+9jscYY4zDkuoYY4wxEc7O7I0xxpgIZ5W9CRgR2S4iF4vIYyLyZqnlU0TkaxE5JCLzRKSXl3EaY04nIr8SkV0ickRENojIRb6Jup4XkVzf7XkRSfCtny4iOSJyv4jsEZHdInKr15/DVMwqexMMn+DM7NcKWAq85W04xpgSItINuAcYrKoNgUuA7cBDwDCgP9APGAI8XGrT1kBjoB1wO/CiiDQNXuSmOqyyN65T1VdV9Yiq5gOPAf1EpLHHYRljHEVAAtBTROqo6nZV3QLcBDyhqntUdS/wOPDDUtud9D1/UlVnAHlAt2AHb/xjlb1xlYjEishTIrJFRA7jnDGAMx2mMcZjqroZ+CnOD/E9IvKOiLQF2gI7Sq26w7esxD5VLSz19zGggcvhmhqyyt647fvAGOBinCa/FN9y8SogY8zpVPXfqjoc6Ago8DSQ6/u7RAffMhOGrLI3bmsI5AP7gETg996GY4wpTUS6iciFvsF3J4DjOE37bwMPi0hLEWkBPAK8WclLmRBmlb1x2ySc5r9dwFogy9twjDFlJABPAd8CX+MMpP018DtgMbASWIUzuPZ3HsVoasmS6piAEZGvgB/gNNknq+ptHodkjDEGO7M3ASIiLYGWOGfxPYFt3kZkjDGmhFX2ptZEZDCwCXgBmAokAy97GpQxxphTrBnfGGOMiXB2Zm+MMcZEOKvsjTHGmAgX53UAgdKiRQtNSUmpcr2jR49Sv3599wMKYbYPonsfLFmy5FtVbel1HJXx93iuTCT/j+2zhSc3Ppu/x3PEVPYpKSksXry4yvUyMjJIT093P6AQZvsguveBiPx/e/cdHld5JX78e1Qs9yLZli03GWzLBTfc6DgUYwjghBLaBsgmMdkNJGyc/cUEFkjbkGQJIQnZQCD0UEICa3oxCIybXCRb7pZl2ZYl9yoX1fP7494x47HKVZm5U87nefRoNHPv3DOjmTlz33LeLU1v5S+v7+fGxPP/2B5bbArHY/P6frZmfGOMMSbOWbI3xhhj4pwle2OMMSbOWbI3xhhj4pwle2OMMSbOWbI3xhhj4pwle2OMMf46cACsdHtYxc08e2OMMTHoww9h2jTo0gVGjWJYRgasWAGjR8MZZ0Dv3n5HGBcs2RvfZc9++6S/M44c4GD7ztQkf/HyLHnoy5EOyxgTCX/+M/TsCTfdBIWF9Jo3D94O+kzo1ctJ+mPHwn33QUaGf7HGMEv2JqqkVVfyyV/uYHOPLL7z1R9T3jWqq7oaY1pjzx5480246y54+GEA5n/yCVNHjoRVq5yfwkLn9+9/7zT1/+53Pgcdm6zP3kSVIXu30bXyCGN3bOTNZ+9mytZCv0MyxoTLSy9BdTXcdtsX14lAZiZcfDF8//vw5JOwaBFcdRX8/e9QV+dfvDHMkr2JKsN3O2WeZ371Xg6278KLL9/L7Uvn2OAdY+LRM8/AmWfCmDFNb3vjjVBWBp9/Hvaw4pElexNVcnaXUJmcytwhk5lx62/55PRJPDj3Cbj9djh2zO/wjDFtZeVKWL7ceW97ceWV0KEDvPJKWMOKV2FN9iIyXUTWi0iRiMyu5/YLRGS5iNSIyHUhtw0UkQ9EZK2IrBGR7HDGaqLD8N0lFGUMoDYpmYq0jsy85l5+e94t8NxzcN55sHWr3yEaY9rCs89CaqozMM+Lzp2dhP/aa1BTE97Y4lDYkr2IJAOPAZcDI4GbRGRkyGZbgduBv9VzF88Bv1HVEcBkYFe4YjXRY9ieLazrNejE3ypJ/P7cm2DOHCgqggkTIDfXvwCNMa1XXQ0vvOD0w/fs6X2/G26AXbvsM6AFwnlmPxkoUtViVa0CXgZmBG+gqiWquhI4acSF+6UgRVU/dLerUNWjYYzVRIFuxw7Tp2If63tln3rjVVdBXp7zwXDJJfDoo9aPb0yseu89J2l7bcIPuOIK5wzfmvKbLZzJvh+wLejvUvc6L4YBB0TknyKSLyK/cVsKTBwbvrsEoP5kD5CTA4sXO4n/7rvhJz+JWGzGmDb0zDNOsZzp05u3X4cOcPXV8M9/Oq0DxrNwzrOXeq7zeiqWApwPjMdp6n8Fp7n/qZMOIDITmAmQmZlJroemnYqKCk/bxbNoew5mjXb638aUFwNw4Xn9mJhxcp/cSfHedRfjNm8m5YUXWDp1aouOGW3PgTEJI3hufWpq8/e/4Qb429/go4/g8svbPr44Fc5kXwoMCPq7P1DWjH3zVbUYQETeAM4iJNmr6hPAEwATJ07UqR4++HNzc/GyXTyLtufgdreC3i9WbuNgWid+uj0Tyk7+rlhyy9STd7roIvjzn5l6wQWQ1PwGqmh7DoxJGPXNrW+Oyy6Dbt2cpnxL9p6Fsxl/CTBURAaLSDvgRmBOM/btISKB8mkXAWvCEKOJIjm7t7CuV7ZTVKPJjXOcqXilpWGPyxjThpozt74+aWnwla/AG29AZWWbhhbPwpbsVbUGuBN4H1gLvKqqq0XkpyJyNYCITBKRUuB64HERWe3uWwv8EJgrIoU4XQJ/CVesJgqoMmz3lob760Pl5Di/168PW0jGmDbW3Ln1DbnhBjh4EN5/v03CSgRhrY2vqu8A74Rcd3/Q5SU4zfv17fsh0MKvfibW9Du0m65VR9kQNO0uWOhiOb0q9rEEuP/X/+C5uVUnrrcFc4yJYs2dW9+QSy6B9HSnKf/qq9smtjhnFfRMVBi2xymTu66BZB9qd6ceHG7XgdP2bQ9nWMaYttLSufX1SU2Fa65x6m9YZU1PLNmbqBCYdrehp7dkjwjF6f0t2RsTK1o6t74hN9wAFRXwzjtNb2ss2ZvokLO7hO1denGofWfP+xSn9+O0fTZAz5iY8Mwzztr0zZ1b35CpU537swI7nth69iYq5OzewnqPTfgBxen9+OqaXNpXH+d4avswRWYCRGQ68CiQDDypqg+F3J6GU+Z6ArAXuEFVS0QkFXgSOBPnM+c5Vf1lRIM3/grMrb/zTkhNPWUMTrBZo2u4ffbbTY+/SUmB66+Hp592zvA7ez9RSER2Zm98l1Jbw+l7SxscnNeQ4nRnbOfg/V7LN5iW8rjWxTeB/ao6BHgE+JV7/fVAmqqOxvkicIctbJVgAnPr26oJP+CGG5w++7featv7jUOW7I3vBu/bTru6GmeOfTMUZzjVl0/ba/32EdDkWhfu38+6l18DLhYRwamc2UlEUoAOQBVwKDJhm6jwzDMwfnzL59Y35LzzICvLmvI9sGZ847svauI378x+c48sAAbvt2QfAfWtdTGloW1UtUZEDgIZOIl/BlAOdAT+Q1X3hR6gJeWvGxPPJZFj6bF12rSJScuXs/HOO9nuxhwokV2fzA7O7V4f35CzzyZrzhzmv/UWtVHelO/n/82SvfFdzp4t1EgSm9IHNL1xkOOp7dnepZeNyI8ML2tdNLTNZKAWyAJ6APNE5KNAOewTG7ag/HVj4rkkckw9tjffhNRUhj7wAEPdKXe3N9Fn/3BhyqklsuuRPfttzuxwFv+s/gev//lj/nnGxfVuFy31N/z8v1kzvvFdzu4SNqf3oyql+Yti2Ij8iPGy1sWJbdwm+27APuBm4D1VrVbVXcB8YGLYIzb+C8ytv/LK1s+tb8DyrOGUdu3FlWvnheX+44Ule+O7nOaUyQ1RnNHP6bO3te3DzctaF3OAwOom1wEfq6rirFx5kTg64SxqtS5CcRs/tfXc+vqI8Pbw8zm/JJ9uxw6H7zgxzpK98dfhwww8uNNz5bxQxen96Vp1lF5HDrRxYCaYl7UucFalzBCRIuAHwGz3+seAzsAqnC8NT6vqyog+AOOPt95yVqgL8+p0bw0/n9S6Wi7bsDCsx4ll1mdv/LV6NUDLz+zT3RH5+0rZ3blHW0Vl6uFhrYvjONPsQverqO96kwAWLIBzzmnZuvXNUNhnCFu69+HKdfN4dey0sB4rVtmZvfFXYSFAs6fdBQTm2tsgPWOizMGDzpf5s88O/7FEeGv4+ZyzZQXpRw+G/3gxyM7sjb9WreJIantKu/Vu0e5lXXtyLCXNBukZE20WL3bG0pxzTot2b6zKXn3eGnE+3130dy5fP58Xx1/RomPGs7Ce2YvIdBFZLyJFIjK7ntsvEJHlIlIjItfVc3tXEdkuIn8MZ5zGR4WFbOw5EJWWvRRVkijp0dfO7I2JNgsWQFISTJ4ckcOt7TWYzT36ctGmJRE5XqwJ25l9UHnNS3Gm5CwRkTmquiZos63A7cAPG7ibnwGfhitG4zNVKCxkXdaZrbqbTen9GbVrUxsFZYxpjcAZ+XOvzKFnz0Fc8YvPInNgEZb0H8UlRXnOZ4vUV/YhcYXzzL7J8pqqWuKOyq0L3VlEJgCZwAdhjNH4aedO2LOn2ZXzQhWn92PAgZ2k1la3UWDGmNYQrWNc2XqW9RsR0ePmZw0n/dghBh0oj+hxY0E4++y9lNesl4gkAQ8DXwfqL4lEy8prxlKZyXCJluegx9KljAXGTOrPrDMaLp/ZlJyDfUlZWMcDvUo9P65oeQ6MiUdD92yla9VRlmcNj+hxl/dzjje+bD1b3HLaxhHOZO+lvGZD/h14R1W3SSNNMS0prxlTZSbDJGqeg/x8AH5ecTp7C1v+UhxzdCDTgM8W7eCJ/3enp32i5jkwJg5N2O7UTIr0mf3GjAFUtOvA+LJ1vDHqSxE9drQLZ7L3Ul6zIWcD54vIv+MU42gnIhWqesogPxPDCguhd2/2dureqrux6XfGRJcJ29eyp2M3tnbvE9Hj1iUls6LvUM7cbgUaQ4Wzz95Lec16qeotqjpQVbNxBu89Z4k+DhUWwujRrb6birSO7OrUw6bfGRMlxpetY3m/Eb4MksvPGs6IXZtpX3084seOZmFL9l7Ka4rIJBEpxamu9biIrA5XPCbK1NY6BTfOOKNN7s5ZEMfO7I3xW4+jBzl93/aI99cHLM8aTorWMXpHkS/Hj1ZhLarjobzmEpzm/cbu4xngmTCEZ/y0eTMcO+ac2W9s/d0Vp/dn+oYFrb8jY0yrjC9bD8Cyfv4k+4KsnBNxLBnQNicT8cDK5Rp/uGVy26IZH2BTej/Sjx2CvXvb5P6MMS0zYftaqpOSWdlnqC/H39exGyXd+3JmmfXbB7NyuSasGip5+b35/+BuhFEvb4N27Vt9nOIMt4Fo/foWl+c0xrTehO1rWZ15GpWpab7FkJ+Vw7lbVlhxnSCW7I0vcnaXsLV7H461QaKHL1a/s2RvjI+qqxmzYyMvj7nM1zDys3L46ppcsg7vpqxrb0919kse+nIEIvOPNeMbX+Ts3sKGVlbOC1baLZOqpBQn2Rtj/LFyJR2rK08Ut/HLcnd+//jt9nkQYMneRFxaTRXZ+8tY17Ptkn1tUjJbevS1ZG+MnxYuBCJfTCfUul7ZHE9px3jrtz/Bkr2JuCF7t5Gidaxv4Rr2DSlO7wcbNrTpfRpjmmHBAso7Z1DetZevYdQkp7CyzxAbpBfEkr2JuJzdJYDz7bstFaf3h6IiZw6/MSbyFi70/aw+ID9rOKN2bqJdjS2QBR6SvYjYREXTpnJ2b6EyOcVpdm9Dxen9oKoKSkra9H6NMR6UlUFJCfk+99cHLM8aTlptDSN3FfsdSlTwcmb/ZxHJE5F/F5HWFTE3BifZb8oYQE1y204G2ZQeNP3OGBNZUdJfH5B/oriONeWDh2SvqucBt+AsarNURP4mIpeGPTITt3J2l7R5Ez5AcUbQ9DtjTGQtXAhpaazOPM3vSADY1SWD7V16najol+g89dmr6kbgPuBHwIXA70VknYhcE87gTPzperyCvhV7Wd+G0+4CDnToChkZluyN8cOCBTBxItXJqX5HckJ+Vo6tgOfy0mc/RkQewVnM5iLgKlUd4V5+JMzxmTgz3B2ct75ndngOkJNjyd6YSKushGXLoq6gVX5WDv0P7aJXxT6/Q/GdlzP7PwLLgbGq+l1VXQ6gqmU4Z/vGeDZs9xaAsJzZOwcYZsnemEhbvtwZHHv22X5HcpLAYEFryveW7K8A/qaqxwBEJElEOgKo6vON7Sgi00VkvYgUicgp69GLyAUislxEakTkuqDrx4nIQhFZLSIrReSG5j0sE62G7y7hUFonyrv0DM8BcnKgvBwOHQrP/RtjTuUOzou2ZL8683SqklIs2eMt2X8EdAj6u6N7XaNEJBl4DLgcGAncJCIjQzbbCtwO/C3k+qPArao6CpgO/M5mAsSHnN1bWNdrUPgWp8hxRuBacR1jImjBAhg8GPr08TuSk1SmtGNN5mlWXAdvyb69qlYE/nAvd/Sw32SgSFWLVbUKeBmYEbyBqpao6kqgLuT6De6gwEB3wS7A35JMpvVUydmzpc0r550kkOytKd+YyFB1kn2U9dcH5GflMKZ8I8l1iV1sy0uyPyIiZwb+EJEJwDEP+/UDtgX9Xepe1ywiMhloB2xq7r4mumQd3k3XyiOsb8Oa+Kc4/XRISrJkb0ykbN3qdJ1FWRN+wPKs4XSoqTwxODhRealqcjfwdxEpc//uC3jpQ6+vnVa9BgYgIn2B54HbVLWunttnAjMBMjMzyc3NbfI+KyoqPG0XzyL5HMwaXXPi8sCVzuC88ZP7kTmipqFdWiV34UKm9OnD4XnzWNPIY7TXgTFtZMEC53e0ntkHDdJbnXm6z9H4p8lkr6pLRGQ4kIOTwNepqpdiw6U4hXgC+gNlDWx7ChHpCrwN3KeqixqI7QngCYCJEyfq1KlTm7zf3NxcvGwXzyL5HNwetI70v+TvYgbw6/392VnYttXzAkpumQrjxtFh+3Z6N/IY7XVgTBtZuBA6dYLRo/2OpF6lXXuzu1N3xpet44XxV/gdjm+8fuJOArLd7ceLCKr6XBP7LAGGishgYDtwI3Czl4OJSDvgdeA5Vf27xxhNlBt0oJzjKe3Y1Tk9vAfKyYFPPoG6OqdJ3xgTPgsWwOTJkBKeL/CtJkJ+1nDGJ3hxHS9FdZ4H/gc4DyfpTwImNrWfqtYAdwLv4xTkeVVVV4vIT0Xkave+J4lIKXA98LiIrHZ3/xpwAXC7iBS4P+Oa//BMNBl0YAdbu/VBJcwJOCcHjh2D0tLwHseYRHfkCBQURG0TfkB+Vg6n7S+j+7HEnZLr5avYRGCkqjarvx1AVd8B3gm57v6gy0twmvdD93sBeKG5xzPRbeD+crb0iMDUnOAR+QMHhv94CUBEpgOPAsnAk6r6UMjtacBzwARgL3CDqpa4t40BHge64sy8maSqxyMXvQmbpUudJaWjdHBewPIsp99+XNl6ck+f5HM0/vByirUKiK7Jkyb2qDLw4A62dm/bZW3rZdPv2pTHmhnfBPar6hCcMtq/cvdNwfni/h23bsZUwBYYjxeBwXlnneVvHE1Y2WcotZKU0MV1vJzZ9wTWiEgeUBm4UlWvDltUJu70OnKAjtWVbO0ege+NffpAly6W7NvOiZoZACISqJmxJmibGcCD7uXXgD+KiADTgJWqugJAVfdGKmgTAQsXOl+uMzL8jqRRx9q1Z12vbEv2TXgw3EGY+DfwQDkAWyKR7EVsQZy2VV/NjCkNbaOqNSJyEMgAhgEqIu/jFMZ6WVV/Hf6QTdipOsn+qqv8jsST/Kwcrl7zKaJ14R83FIW8TL37VEQGAUNV9SO3Ln5y+EMz8WSQm+wj0owPTrKfNy8yx4p/XmpmNLRNCl8M7j0KzBWRZao695SDtKBuRmPiuZZCNDy2DqWlTNmzh/UZGZQHxRJcW6MlMju0/j7qvd8Dw+ha8C4/Ty9hX/9Tx/JE4vn08//WZLIXkW/jvAHTgdNxvsH/Gbg4vKGZeDJo/w7qEEq7ZUbmgDk58OKLcPQodPRS3dk0wkvNjMA2pW4/fTdgn3v9p6q6B0BE3gHOBE5J9i2pm9GYeK6lEBWP7dlnAci5/XZyRo06cXVwbY2WmDW6hofDUIdjcN1ILgVW5Bbx6tjTTrm95JapbX7MUH7+37y0ZXwXOBc4BODWrO8dzqBM/Bl4oJyyrj2pSkmNzAEDg/SKiiJzvPh2omaGWwPjRmBOyDZzgNvcy9cBH7szeN4HxohIR/dLwIWc3NdvYtWCBdCtG4wY4XcknmzukcWB9p0Zn6CL4nhJ9pXuQjbAidG1zZ6GZxLboAPlkWvCBxuR34a81MwAngIyRKQI+AEw2913P/BbnC8MBcByVW3dqZ+JDosWwZQpsVO4SoT8rJyEXQHPS1vJpyLyY6CDiFwK/DvwZnjDMvFm4IEdfDgkdExXGA0d6vy2ZN8mPNTMOI5THKu+fa1uRrw5cgRWrYIZM5reNorkZw3nwuLldKk8wuG0Tn6HE1Fekv1snDm0hcAdOG/4J8MZlIkvnSqP0vPoQbb2iOCZfceOTkEdS/bGtInsoL74ydtW8WpdHd9Ym8Qnreyjj6TlWcNJQhlTvpH52YlVlNXLaPw64C/ujzHNNvDgDgC2RLIZH2z6nTFhMrZsA+AUq4klK7KGUYcwvmydJftQIrKZevroVfXU4YzG1GPg/kCyj3AhxpwceO45Zz6w1DczzBjTEmPLN7CtWyZ7O3X3O5RmOZzWiaKMAQlZXMdrbfyA9jj9cmFetszEkxNz7CPZjA9Ov/2hQ7BnD/TqFdljGxPHxpWvJ9+tNx9r8rNyuLRoccKdBDQ5jFJV9wb9bFfV3wEXRSA2EycGHShnf/sukR8QM2SI83vjxsge15g41qtiP/0P7aag7zC/Q2mRgqwc0o8dOnESkii8LHF7ZtDPRBH5DtDFy52LyHQRWS8iRSIyu57bLxCR5SJSIyLXhdx2m4hsdH9uC93XxI6B+3dEZrW7UIER+ZbsjWkzY3Y4/fUFWTk+R9IyBVnOl5RxCdaU76UZ/+GgyzVACc56840KWinrUpwqWktEZI6qBhfU2ArcDvwwZN904AGcLgQFlrn77vcQr4kygw6U+/PBkJ0NyclWWMeYNjSubAM1ksTqzNgctrWh5yCOpqYxrnwD/zfqS36HEzFeRuO39NlocqWsoPWu60L2vQz4UFX3ubd/CEwHXmphLMYnKbU1ZB3azf+NnBr5g6emOgnfzuyNaTNjyzewvlc2x1Pb+x1Ki9QmJVPYZ2jCDdLzMhr/B43drqq/beAmLytlNaS+fft53NdEkX6HdpGidZFZ2paT5wIDPFPXnYzcpVwVdH3JQ1+OSCzGxBvROsaWb+CtEef7HUqr5PcdxjeWzaFdTXXkSnj7zOto/El8UQv7KuAzTk7G9fGyUlar9m3JKlnRsFqU3yL5HPyg+3YAzprQm4HD234lq6b0yO/DsHlrmXVG9YmRt7m5ufY6MKYFBu8ro1vlkZgdnBdQkJVDWl4NI3cVx+zYg+bykux7Ameq6mEAEXkQ+LuqfquJ/byslNXYvlND9s0N3aglq2RFxWpRPovkc3Dff7/KDODX+/qxMwwrWTVlb20/Hjx2lGcWHzkxJ7jklqn2OjCmBcaWO4PzVsR6su/rJPhxZesTJtl7WcFgIFAV9HcVkO1hPy8rZTXkfWCaiPQQkR7ANPc6E2MGHtjB8ZR27OrsT2mGkh5ZAAzan1jTbIwJh3Hl66lo14GijAFNbxzFdnTtyY7O6YwrT5x+ey/J/nkgT0QeFJEHgMXAc03t5GWlLBGZJCKlOIV6HheR1e6++4Cf4XxhWAL8NDBYz8SWQQfK2dqtDyr+rIxV4hbyGbzfa6OSMaYhY8s3UNhnCHVJyX6H0moFWTmMc8v+JgIvo/F/ISLvAoERGd9Q1Xwvd+5hpawlOE309e37V+CvXo5jotfAAz7NsXeVdsukRpLItmRvTKu0q6lm5M7N/HVSbK1015CCvjlM37CQHkcPsr9jN7/DCTuvp1sdgUOq+ihQKiKDwxiTiReqDDywI7Lr2IeoSU5hW/dMS/bGtNLIXcW0q6shv2989HEHiusExiHEOy8V9B4AfgTc416Viq1NbbzYuZNO1ccjNu2uIVu6Z1myN6aV4mVwXsDKPkOplSTGJ0hTvpcz+68CVwNHAFS1DI/lck2C27QJ8GG1uxCb091kr15nfhpjQo0t38DOzuns6JLhdyht4mi7DmzoOTBhBul5SfZVqqq489xFJMKrmZiY5SZ7P5vxwRmR36XqGBlHD/oahzGxbFzZemd+fRytFJefleO0WCTAiYCXZP+qiDwOdBeRbwMfAX8Jb1gmLmzaRB1CabdMX8PY4n7ZsKZ8Y1po3z5O218WN034AQV9c+h+vCIhZut4GY3/PyJyKXAIyAHuV9UPwx6ZiX2bNlHWtafv5Sg3pztz7QfvL2NZ/5G+xmJMTFqyBHDOhONJ8Ap4oaW26xPLpbYbTfbuynXvq+olgCV40zzFxb434QNs79qbGkmywjrGtFReHnUIq/oM8TuSNlWUMYCKdh0YV76e18+4yO9wwqrRZnxVrQWOikj8T0I0bW/TJt8H58EX0+8SoanOmLDIy2NTRn8Op8XXkK26pGRW9hmaEMV1vBQrPw4UusvMHglcqarfC1tUJvYdPgy7drF1hP9n9uAM0htkyd6Y5lOFxYsp6DvW70jCoiBrGN/K66BioAAAIABJREFUe4O0mioqU9r5HU7YeEn2b7s/xnhXXAx8MTjObyU9sphYuiYhRt0a06a2bIHdu1kxPr4G5wUU9M2hXV0No3ZuYnm/EX6HEzYNJnsRmauqFwMjVfVHEYzJxIMomWMfEJh+1/PoAb9DMSa25OUBxPyytg0JDDocV7YhMZM90FdELgSuFpGXCVljXlWXhzUyE9sCc+x7RM+ZPdjqd8Y02+LFkJbGul7ZfkcSFrs7p7O9S6+4L67TWLK/H5iNs1DNb0NuUyC+hy6a1tm0CdLTo2ZAj61+Z0wL5eXBmWdSk+yl1zc2FWQNY3xZfCf7Bkfjq+prqno58GtV/VLIjyV607hNm+D00/2O4oTSbplUJyVbYR1jmqO6GpYtg8mT/Y4krAr65jDg4E4yjsRvN1+TFfRU9WctvXMRmS4i60WkSERm13N7moi84t6+WESy3etTReRZESkUkbUick/ovibKRVmyr01KprRbb7KtGd8Y71avhmPHYMoUvyMJqxPFdeK4Kd/rErfN5hbkeQy4HBgJ3CQioeXLvgnsV9UhwCPAr9zrrwfSVHU0MAG4I/BFwMSA6mrYujWqkj04/fZ2Zm9MM7iD8+L9zL6wzxBqJCmu59uHLdkDk4EiVS1W1SrgZWBGyDYzgGfdy68BF4uI4IwJ6CQiKUAHoAqnXK+JBVu3Qm1tVCb7Qbb6nTHeLV4MGRlw2ml+RxJWx1Pbs75XNuPiuN/e04gLETkPGKqqT4tIL6Czqm5uYrd+wLagv0uB0LagE9uoao2IHAQycBL/DKAc6Aj8h6ruqyeumcBMgMzMTHJzc5t8LBUVFZ62i2fhfg56LFnCWCD/0CFmjU4P23Ga67TyTLosO8b811+nol27hH8dGNOkvDznrD6OVrprSEHWMK5a8xmidaiE8zzYH00mexF5AJiIswjO00Aq8AJwblO71nNd6ClVQ9tMBmqBLKAHME9EPlLV4pM2VH0CeAJg4sSJOnXq1CZCgtzcXLxsF8/C/hysXQvA+Ouu46t/KAjfcZrpwmMDuBA4t3dvcmtqEv51YEyjDh92+uyvvdbvSCKioG8OtxS8x2l7t7Op5wC/w2lzXr6+fBW4GrdUrqqWAV087FcKBD9j/YHQDtMT27hN9t2AfcDNwHuqWq2qu4D5OF84TCzYtAnat4e+0THHPmCzO9eeoiJ/A4lRLR1wG3T7QBGpEJEfRipm0wrLljldXnE+OC8gUFxnfJwO0vOS7KtUVXHPykXE68TpJcBQERksIu2AG4E5IdvMAW5zL18HfOweaytwkTg6AWcB6zwe1/ht0yanjy8puprCtnfrTXVSMmzc6HcoMaeVA24DHgHeDXespo0EBudNmuRvHBGyKaM/h9p1jNt+ey+fxq+KyONAdxH5NvAR8JemdlLVGuBO4H1gLfCqqq4WkZ+KyNXuZk8BGSJSBPwAp4gPOB8qnYFVOF8anlbVlc14XMZPUTbtLqA2KZlt3TIt2bdMawbcIiJfAYqB1RGK17RWXp7zpb1nT78jiQiVJFb2Hcq48vgckd9kn72q/o+IXIozGj4HuF9VPa1tr6rvAO+EXHd/0OXjONPsQverqO96EwNUnUVwLr7Y70jqtaVHX06zZvyWaPGAWxE5BvwIuBSwJvxYsXgxnH++31FEVEFWDt9Z9Brtq49zPLW93+G0KS8D9P4D+LvXBG8S3M6dcORI1E7VKemRBes/sel3zdeaAbc/AR5R1QppZFR3S2bXNCaeZ96E+7G127OHc0pLKUpPp9Q9zqzRNWE7XrDMDpE7VqjsqiGkLKzj5502UJ4T2ktFTL8mvUy96wq8LyL7cJruXlPVneENy8QsdwGcaGzGB3eQXkUFqfv3+x1KrGnOgNvSkAG3U4DrROTXQHegTkSOq+ofg3duyeyaxsTzzJuwP7Y33gBgyM03M+SccwC4fXZkVjqfNbqGhwv9qcPfs3oEVwLrPi/iyaoxp9xecsvUVt2/n69JL834PwF+IiJjgBuAT0WkVFUvCXt0Jva469hHa7Lf4o7I77h9u8+RxJwTA26B7TgDbm8O2SYw4HYhJw+4PdEWLCIPAhWhid74Kzskkf/np39jZlIyZ/xjB5VzIpPko8GeTj0o7do7LivpNefr0y5gB7AX6B2ecEzM27TJKcCRne13JPUKTL/rUFrqcySxxe2DDwy4TQb+GhhwCyxV1Tk4A26fdwfc7sP5QmBi0NjyDazrlU1laprfoURcflYO48vib/KXlz77f8M5o++FM8L226q6JtyBmRi1aRMMGABp0fkhsb1bb0hJoYOd2TdbSwfchmz/YFiCM20mqa6WMeUbmTPyQr9D8UVBVg5XrZtHr4r97O7cw+9w2oyXM/tBwN2qGj2l0Ez0itJpdwG1SckweLCd2RvTgNE7iuhadZTFA87wOxRf5Pd1iuuMK1/Ph0PP8jmattPgPHsR6epe/DWwVUTSg38iE56JOVGe7AEYOtTO7I1pwLlbVgCwYNBYnyPxx+rM06hOSo674jqNndn/DbgSWIYzfSZ4zowC0Tm3yvjn8GHYtSs2kv0n7vS7BFjgw5jmOHdLAWt7ZbO3U3e/Q/FFZWoaa3sPjru17Rs8s1fVK93fg1X1NPd34McSvTlVlI/EP2HIEFKOHXNqAhhjTkirrmRi6VrmJ+hZfUBB3xzGlG8kqa7W71DaTJPlckVkrpfrjIn2OfYnDB3q/LayucacZML2taTVVjM/e5zfofhqeb/hdKk6xohdTa3kHjsa67Nv7/bN9xSRHkH99dk4S88ac7JYSfZDhji/rWyuMSc5d8sKqpOSyes/yu9QfBVo2Ti/JH7GpTd2Zn8HTn/9cPd34Of/cBaqMeZkxcWQng7duvkdSeMGDaIu2Va/MybUuVsKKOibw5G0jn6H4qvdndNZ2yub80uW+x1Km2msz/5RVR0M/DCkz36sVb8y9YqFkfgAKSkcz8qyM3tjgnQ9XsHoHZsSdhR+qM8Gn8nE0jV0qDrudyhtosk+e1X9g4icISJfE5FbAz9e7lxEpovIehEpEpHZ9dyeJiKvuLcvdrsIAreNEZGFIrJaRApFJL6WIIpHsZLsgWNZWXZmb0yQs7YWkqx1zM+2ZA8wL3s8abU1TNm2yu9Q2oSXAXoPAH9wf76EM+/+6kZ3cvZLxmnuvxwYCdwkIqHLCH0T2K+qQ4BHgF+5+6YALwDfUdVRwFSg2ttDMr6oroYtW2In2ffv7yR7W/3OGADO2bKCo6lp5Gfl+B1KVFjSfyTHU9pxfkm+36G0iSaTPc6CFhcDO1T1G8BYwEst1MlAkaoWq2oVzop5M0K2mQE8615+DbhYnDUwpwErVXUFgKruVdX4mQMRj7ZuhdramEn2R/v1c5bitel3xgBwXkkBef3PoDo51e9QokJlahp5/Udx/ubESfbHVLUOqHGr6u3CW0GdfsC2oL9L3evq3UZVa4CDQAYwDFAReV9ElovI//NwPOOnWBmJ7zrWz30pWlO+MWQe3sOQfaUJP78+1GeDxzNs71b6HNrjdyit5qU2/lIR6Q78BWc0fgWQ52G/+kqThbaZNrRNCnAeMAk4CswVkWWqetL8fhGZCcwEyMzMJDc3t8mgKioqPG0Xz8LxHGS9+y7DgAU7d1IVdN+zRte06XHaQm5uLnXpTsXndW++yY5aazQyie2cLSsBWGD99Sf5PHs8AOdtKeC10bG9qruX9ez/3b34ZxF5D+iqqis93HcpMCDo7/5AWQPblLr99N1wlsYsBT5V1T0AIvIOcCZwUrJX1SeAJwAmTpyoU6dObTKo3NxcvGwXz8LyHLz0EnTvzjnXXXdSCdrbZ0ffWtglt0zl09paSElheEoKwxP89WDMuVtWsK9DV9b0Hux3KFFlXa9sdnfqzvmb82M+2TdWVOfM0B8gHUhxLzdlCTBURAaLSDucta3nhGwzB7jNvXwd8LGqKs6a2WNEpKP7JeBCwJbVjWZ5eTB5cszUmtfkZDjtNGvGN0aVc0sKWDBwDCpeenYTiAifZY/nvJJ8ROv8jqZVGjuzf7iR2xS4qLE7VtUaEbkTJ3EnA39V1dUi8lNgqarOAZ4CnheRIpwz+hvdffeLyG9xvjAo8I6qRt8ponEcPQqFhXDPPX5H0jxDh9pce2M2bKBvxV5rwm/AvMFncu3qTxi5s9jvUFqlwWSvql9q7Z2r6jvAOyHX3R90+ThwfQP7voAz/c5Eu+XLnZH4U6b4HUnzDBkCubm2+p1JbHOd3lEbnFe/+YOcdQIuiPEpeE322YtIR+AHwEBVnSkiQ4EcVX0r7NGZqJbt9sd/M+91/guY+P4B9nweQw0wQ4c60+927IC+ff2Oxhh/zJ1LadfebOlu74H67O7cgzW9B8f8FDwvHTRPA1XAOe7fpcDPwxaRiTnjyjdQ2rU3ezr18DuU5gmsfmdN+SZR1dbCxx87Z/XWutWgz7LHM7F0jXNyEKO8JPvTVfXXuBXsVPUY9U+ZMwlqXPkGCvoO8zuM5gusfmeD9Eyiys+HAwesRG4T5g0+k3Z1NfDpp36H0mJekn2ViHTAnSMvIqcDlWGNysSMjCMHGHBwJytiMdkPHAipqZbsTeJy++tt8ZvGLXVL5/Lhh36H0mJekv0DwHvAABF5EWeuu1W0MwCM2eEkyoKsGEz2KSkweLA145vENXcunHFG7HXBRVhlSjvy+o+CDz7wO5QWazTZu3Xq1wHXALcDLwETVTU37JGZmDCubAO1ksSqzCF+h9IyQ4famb1JTMePw7x5cPHFfkcSEz4bPB7WrIHSUr9DaZFGk71b4OYNdyGat1X1rUBVO2MAxpZvYEPPgRxrF6MrEAfm2tvqdybRLFzoJHxL9p7MG+zWkovRpnwvzfiLRGRS2CMxsUeVsbE6OC9gyBBnhG15ud+RGBNZc+dCcjJceKHfkcSE9T0HQZ8+MduU7yXZfwlYKCKbRGSliBSKiJfa+CbODTpQTo/jh2NzcF7AWHdgUp6XtZ2MiSNz5zolrrt29TuS2CAC06Y5Z/Z1sVc610uyvxw4Hac87lXAle5vk+DGlm8AoCArx+dIWmHSJOjQwamkZ0yiOHjQ+YJrTfjNM20a7N3rTFmMMV5WvdsSiUBM7BlXtoGjqWls7DnQ71CaJXv228waXXNiRb4Xeg8j/aU5XNH+0hPblDz0Zb/CMyb8Pv3UOTu1ZN88l7gr333wAUyY4G8szWRLHJkWG1e+nsLMIdQmJfsdSqssHnAGw3eV0O3YYb9DMSYy5s51WrTOPtvvSGJLZqbT9ReD/faW7E3LVFUxamdxbPfXuxYNHE0SyuTS1X6HYkxkzJ0L558PaWl+RxJ7pk2D+fNjrnRuWJO9iEwXkfUiUiQis+u5PU1EXnFvXywi2SG3DxSRChH5YTjjNC1QWEhabXVs99e7VvTN4XhKO87aWuh3KMaEX3k5rF5tTfgtNW0aVFfHXOncJvvsW0pEkoHHgEtxFs9ZIiJzVHVN0GbfBPar6hARuRH4FXBD0O2PAO+GK0bTCu7o9Xg4s69KSWVZv+GW7E1cCqxOGTBj9Sc8Cnx5TRqrZ8fQKpXR4rzzoH17pyn/iiv8jsazcJ7ZTwaKVLVYVauAl4EZIdvMAJ51L78GXOxW7UNEvgIUA9a2Go3y8tjdsTvbu/byO5I2sWjAaEbs2kzX4xV+h2JMWJ27ZQX723dhTeZpfocSm9q3d2oTxFi/fTiTfT9gW9Dfpe519W6jqjXAQSBDRDoBPwJ+Esb4TGssXsyKvkPjZlnMxYF++2323bI+Le2SE5FLRWSZW59jmYhcFOnYTRBVzi1ZwcKBo1GxIVstNm0arF0L27Y1vW2UCFszPvUvgxtak7ShbX4CPKKqFdJIMhGRmcBMgMzMTHI9zJWuqKjwtF08a+1zkFxRwXnr1tH5mhuZNbqm7QKLoMwOnBR7cs5p1Py9Hd89uoKxoyck/GskWCu75PYAV6lqmYicAbzPqV/6TYScvq+Ufod3879nX+93KLFt2jTn94cfwr/+q7+xeBTOZF8KDAj6uz9Q1sA2pSKSAnQD9gFTgOtE5NdAd6BORI6r6h+Dd1bVJ4AnACZOnKhTp05tMqjc3Fy8bBfPWv0cfPwxqPJY8gjmFYbzJRQ+s0bX8PBJsacwru9wuuav5uGxKZTcMtWv0KLRiS45ABEJdMkFJ/sZwIPu5deAP4qIqGpw9ZHVQHsRSVNVWybbB9ev/JAaSeKDIVP8DiW2jRoFffs6TfkxkuzD2Y6zBBgqIoNFpB1wIzAnZJs5wG3u5euAj9Vxvqpmq2o28Dvgv0MTvfGROzhvZZ+hPgfSthYNHM3IncXWb3+qFnfJhWxzLZBvid4f7Wqqub7wIz4aOoVdXUL/NaZZRODSS+Gjj2KmdG7YTstUtUZE7sRptksG/qqqq0Xkp8BSVZ0DPAU8LyJFOGf0N4YrHtOGFi+GoUM52KGL35G0qcUDzrB++/q1pkvOuVFkFE7T/rQGD9KCbrnGxHOXXXMeW6C7aujCz8k4dojKr0yL6u630C62aBL8nPfu35+Re/ey7IknODx8uKf9/XxNhrUNVlXfAd4Jue7+oMvHgUY7j1T1wbAEZ1ouLw++9CW/o2hzBVk5VCanMmWbTcEL0ZouOUSkP/A6cKuqbmroIC3plmtMPHfZNeexBcpCv/TmB2ztlsndTEALo3dw3qldbNHjpO690aPh4YeZsHIlfOc7nvb38zUZvf9xE522b4eyMme1rDhTmdKO5Tbfvj4t7pITke7A28A9qjo/YhGbk5y2t5Sztxby8tjLbBR+W8nIgJtvhmefhQMH/I6mSfZfN80TWAp2SnwO8Fk0YDSjdhbHxJs3Utw++ECX3Frg1UCXnIhc7W72FM602SLgB0Bget6dwBDgv0SkwP3pHeGHkPBuWvEe1UnJ/H30pU1vbLy76y44ehSeftrvSJpkyd40z+LFkJr6xTrwcSZQJ5958/wOJaqo6juqOkxVT1fVX7jX3e+OvUFVj6vq9ao6RFUnB0buq+rPVbWTqo4L+tnl52NJNGk1VVxXOJcPhp7F7s49/A4nvowfD+eeC3/8I9TW+h1NoyzZm+bJy3MSffv2fkcSFoF+e1vf3sSLyzYsoMfxw7w0drrfocSn730Piovh3eiu7B6doyBMdKqthaVL4etf9zuSsKlMaUd+Vg5nWbI3ceKWgvfY0r0P87PjszUukkLXGQBIqU1jXucMNt51P7d+LpQ89GUfImuandkb79avh8OH47a/PmDRwNFQUGD99ib2rV3LlG2reGnsdBuYFyY1ySm8MP4KLijJ5/S90Vs+1/77xrvA4Lw4HIkfbNHA0U6hjM8/9zsUY1rniSeoSkrh76Mv8TuSuPby2MuoTE7h1uVv+R1KgyzZG+8WL4auXWFY7C9r25j8rOGQlmb99ia2HTsGzz7LB8POZm+n7n5HE9f2durOWyMu4LrCuXDwoN/h1MuSvfEuLw8mTYKk+H7ZVKa0c7oqLNmbWPbaa7B/Py+Os4F5kfD0hKvpVH0cnnnG71DqFd+f2qbtHDsGK1fGfRP+CVOnQn5+1H5LN6ZJjz8OQ4eycOAYvyNJCKv6DGFZ1nBnGl4U1su3ZG+8KSiAmpq4H5x3wtSp1m9vYtfq1TB/Psyc6SzaYiLi2QlXQVERvPee36GcwpK98WbxYud3opzZn3UWtGtnTfkmNj3+uPP6vf12vyNJKO/mnOMsffuHP/gdyiks2Rtv8vKgf3/nhZwIOnSwfnsTm44eheeeg2uvhZ49/Y4moVQnpzqL4rz3njNVOYpYsjfe5OUlzll9wNSpsHy59dub2PLqq85r9o47/I4kMd1xh1NS/LHH/I7kJGFN9iIyXUTWi0iRiMyu5/Y0EXnFvX2xiGS7118qIstEpND9fVE44zRN2LsXNm1KzGRv/fYm1jz+OOTkwAUX+B1JYsrMhBtucBbHOXTI72hOCFuyF5Fk4DHgcmAkcJOIjAzZ7JvAflUdAjwC/Mq9fg9wlaqOxlk28/lwxWk8WLLE+Z0og/MCAv32n37qdyTGeLNyJSxaZAPz/HbXXVBR4Sx/GyXCeWY/GShS1WJVrQJeBmaEbDMDCDwbrwEXi4ioar6qlrnXrwbai0haGGM1jVm82PngmDDB70giq2NH67c3seWJJ5yCULfd5nckiW3yZOezI4qm4YUz2fcDggsFl7rX1buNu2b2QSAjZJtrgXxVrQxTnKYxdXXwxhswejR06eJ3NJF34YWwbFlUNccZU68jR+D55+G66yAj9GPURNxdd8GGDfDBB35HAoR31bv62pC0OduIyCicpv1p9R5AZCYwEyAzM5NcD2dgFRUVnraLZ815Dnp//DEjCwpYe8897AzZZ9bomrYPLkIyOzQef+D56d69O+Pq6lj5pz+x76yzIhSdMS3wv//rfCm1gXnR4frrYdYsZxredP+rGIYz2ZcCA4L+7g+UNbBNqYikAN2AfQAi0h94HbhVVTfVdwBVfQJ4AmDixIk6derUJoPKzc3Fy3bxzPNzUFUF3/oWjBnDiJ/9jBHJySfdfHs9yz3Gilmja3i4sOGXf8ktU50LkyfDPfcwZt8+Z8CeMdFo6VL48Y/h6qvhvPP8jsaAM97nO9+Bn/wENm6EoUN9DSecyX4JMFREBgPbgRuBm0O2mYMzAG8hcB3wsaqqiHQH3gbuUdX5YYzRNObJJ51R+G+/DSGJPmFYv72JYtmz36ZL5RHeeub7pLTvzhXZN3Hwnnf8DssE3HEH/Pd/w49+BP/4h6+hhC3Zq2qNiNwJvA8kA39V1dUi8lNgqarOAZ4CnheRIpwz+hvd3e8EhgD/JSL/5V43TVV3hSte84Xs2W/TseoYnz5+L8UDzuCGT+vgs9g9i2+1qVPhl790mki7dvU7GmO+oMov3/0D/Q7u4ms3/4qDHRJwXE0069vXSfb/+Z/OZ8g55/gWSjjP7FHVd4B3Qq67P+jyceD6evb7OfDzcMZmGvfNJW/Q6+gBZl5zX0JO4ckO6qI4t6Q9L9bWMvP2X/PBsLNPXF/y0Jf9CM2YE24peJcr13/OQxfezvL+I/wOx9Rn1iynONd995H+i1/41h1oFfTMKdKPHuSOvH/y3rCzye833O9wfLe0/0i2dcvk3k+eomPVMb/DMcaxYgX3z/0LuYMn8PiUa/yOxjRExOkSHTeOkb/4hW9ldC3Zm1PcueAVOlRX8psLbvU7lKhQmdKOH3z5PxhwYCf3ffyk3+EYQ/LRo/C1r7G/Qxd+cOUPULGP8qjWsSO8/jp1KSkwY4YvJbjtFWJOtnkz/5L/Dq+OvoRNGQOa3j5BLBlwBo9PuZabV7zPxUWL/Q7HJDJVhj3yCBQV8f2r/pN9Hbv5HZHxYtAg1jz4oLME7te/HvFiO5bszcnuv5+6pCR+d17oxAnzyHm3sKb3YB569w9kHDngdzgmUT3zDJkffQQPPMDigaP9jsY0w4Fx4+B3v4M334QHH4zosS3Zmy+sWAEvvsjTE65mZxdbGjNUVUoqd185i66VFfzy/T+ChtaIMibM1qyB736X/ePHw733+h2NaYnvfhe+8Q342c/gn/+M2GEt2Zsv3HMPdOvG/551nd+RRK0NvbL59QW3MW3jImdVK2Mixe2np0sX1t57b+LWvoh1IvCnPzkFu269FVatishhLdkbR24uvPsu/PjHHGrf2e9ootpfJ81gwcAx8P3vQ3Gx3+GYRPG97zln9i+8QJXVvo9t7ds7Z/VdusBXvgL79oX9kJbsjdMc/aMfQf/+cOedfkcT9VSS+OGX74akJOebeW2t3yGZeKYKjz4KTz3ltL5deqnfEZm20K+fU1Vv61a46aawf46EtaiOiRGvvw55ec6HSYcOfkcTE8q69obHHnNG1f7mNzB7tt8hmXi0ahX827/B55/DZZc5ddZNVMtuZM2QWaNruH32218U5DrnHKdJ/9vfdtY2+NWvwhaXJftEV1PjvMhGjHDOUo13t9wCc+bA/fc7H8Tjx/sdkYkXFRU8fvFtfHPJGxxO68Qvp3+Pv4+5BL3vfeCLpGHiwLe+5VTY273bmY6XFJ4Gd0v2ie7pp52KTm+8ASn2cmgWEWdZ0c8/h3/5F2fd+/bt/Y7KxDJVp6Xt+9/njtJSXh4zjV9deBv7bS59fPvDH5wkH8bS5NZnn6jKypyazXffDWef7SyNaZovI8P5wrRmjdNCYkxLFRfDlVfCtddCejrX3PIbZl/+PUv0iSA5OexrkFiyTzTFxQx7+GEYPNgZ9HPNNfDKKwm52E2buewyZ2DjI4/AX/4CR474HZGJJZWVzpzrUaPgs8/gt7+FZctsYRvTpsLabisi04FHcZa4fVJVHwq5PQ14DpgA7AVuUNUS97Z7gG8CtcD3VPX9cMYa91atgoceovZvL5GZnMQLZ1zK41OuZVv3PvDYSmCl3xHGnOCBOO3bTeW1zHc5Y+ZMjn73Lj4+fTJvDT+PP794n1MXOw7Y+7mNVFdDQYHT/fP5506C37MHrr/e+cLYr5/fEZo4FLZkLyLJwGPApUApsERE5qjqmqDNvgnsV9UhInIj8CvgBhEZibO2/SggC/hIRIapqs1xaq4lS5z1lN94Azp14qmJM0i59Sp+Wtrb78jiyvHU9lx962+ZVLqGL6/7nMvXz+fKdfOg9++dLpKvfQ2mT4/ZPn17P7fC4cOwaNGJ5H503nw6VlcCsKV7H5b2H8PrF32JzwePhz8UAAX+xmvCrrER+8HachntcJ7ZTwaKVLUYQEReBmYAwR8OM4AH3cuvAX8UEXGvf1lVK4HNIlLk3t/CMMYbW1Sd5uLduxv+2bgR5s+HHj3ggQfgrrv4798sYlaPGufj2rSpuqRkFg8czeKBo/nJJTOZvG0VL3Xd4sylfeklp4DG1Vc761l37+783bWr8xO43KVLtFZGS+z3c10dVFU5Z+XV1c7l48edYih79sDevbBnD7+P28udAAAHlElEQVR7ZSE9jh2ix7HD9Dh6iF5H9jNk7zZStI5aSWJt78EsGTONJf1HsbTfCHZ1seI4JjLCmez7AduC/i4FpjS0jarWiMhBIMO9flHIvm3StlW4/eCJKSs3FbzHzLx/nLLN4IxObXGoU9VXS72h62prnQ+YwO/gy7W1zodNZWW9h6lMTmFfh27s6dSdN6d+gxfHXc6R4x3hN4vq3d60vdqkZBYOGks2Y0n5xhWctbWQK9fNY/o//o/uL77Y+M4dO0KnTl+Mzg38Dr4c+P3MM3D++ZF4SFH5fj7JyJHO+8IJgCnHj7PlaB0S9B4TcP9WREFQkrTuxPWi6l6npNZWk1pXS2ptDcnqbYWyu4H97buwr2NXDrTvwrbumbw/7GyW9h9JftZwKtLio0vHxJ5wJvv6RnyFZraGtvGyLyIyE5jp/lkhIus9xNUT2APwkPtziv0e7iWa1dZAxV7nZ+cmyD25hvv3gp6DRBXp52AT0ESK/8LRo86PFxdc0JJwBrVgn2h9PzfGn9f58cPOT0BRXpsfIp7fw/bYviDeaux4ej+HM9mXAsELovcHyhrYplREUoBuwD6P+6KqTwBPNCcoEVmqqhObs0+8sefAnoMWiMr3c2Pi+X9sjy02+fnYwjn1bgkwVEQGi0g7nAE6c0K2mQPc5l6+DvhYVdW9/kYRSRORwcBQoO2/HhtjvLL3szExLGxn9m6f3Z3A+zhTdf6qqqtF5KfAUlWdAzwFPO8O2NmH8wGCu92rOIN/aoDvJszIXWOikL2fjYltovUNEItjIjLTbS5MWPYc2HOQCOL5f2yPLTb5+dgSLtkbY4wxicbK5RpjjDFxLmGSvYhMF5H1IlIkIgmx+LiI/FVEdonIqqDr0kXkQxHZ6P7u4WeM4SYiA0TkExFZKyKrReT77vUJ9TwkEhG53v1f14nIxJDb7nE/A9aLyGV+xdga8fRZFq+fUdH4uZMQyT6o1OflwEjgJreEZ7x7Bpgect1sYK6qDgXmun/HsxpglqqOAM4Cvuv+7xPteUgkq4BrgM+Crwwp2zsd+JP72RAz4vCz7Bni8zMq6j53EiLZE1TqU1WrgECpz7imqp/hjIoONgN41r38LPCViAYVYaparqrL3cuHgbU41dsS6nlIJKq6VlXrK8hzomyvqm4GAmV7Y0lcfZbF62dUNH7uJEqyr6/UZ6IuLZWpquXgvCCBhFkRR0SygfHAYhL4eUhg8fA5EA+PoSlx9d6Mls+dsC5xG0U8les08UtEOgP/AO5W1UPO+iwmVonIR0Cfem66V1X/r6Hd6rku1j4H4uExJIxo+txJlGTvqVxngtgpIn1VtVxE+gK7/A4o3EQkFecN96Kq/tO9OuGeh3iiqpe0YLd4+ByIh8fQlLh4b0bb506iNON7KfWZKIJLmt4GNHQWFBfcJVafAtaq6m+Dbkqo58EA8VG2NxE+y2L+vRmNnzsJU1RHRK4AfscXpT5/4XNIYSciLwFTcVZa2gk8ALwBvAoMBLYC16tq6ACZuCEi5wHzgEIgsE7pj3H6zxLmeUgkIvJV4A9AL+AAUKCql7m33Qv8K85o6btV9V3fAm2hePosi9fPqGj83EmYZG+MMcYkqkRpxjfGGGMSliV7Y4wxJs5ZsjfGGGPinCV7Y4wxJs5ZsjfGGGPinCV7Y4wxJs5Zsk8gIpIhIgXuzw4R2R70dzu/4wslIr8TkQvcy6ki8pC7NOQqEckTkcvd2z6KxWUwjWkNez+b5rB59glKRB4EKlT1f/yOpT4ikg68o6pnuX8/BPQFZqpqpYhkAheq6qsichvQP5aLixjTGvZ+Nk2xM3tTLxHJFpF1IvKk+837RRG5RETmu9/GJ7vbdRKRv4rIEhHJF5EZQfvPE5Hl7s857vVTRSRXRF5z7/9FqX91iOuA99x9OgLfBu5S1UoAVd2pqq+6284BbgrvM2JM7LL3s7FkbxozBHgUGAMMB24GzgN+iFP6EeBe4GNVnQR8CfiNiHTCWeDhUlU9E7gB+H3Q/Y4H7gZGAqcB59Zz7HOBZUFxbFXVQ/UFqar7gTQRyWjh4zQmEdj7OYElyqp3pmU2q2ohgIisBuaqqopIIZDtbjMNuFpEfuj+3R6n7nMZ8EcRGQfUAsOC7jdPVUvd+y1w7+vzkGP3BXY3I9ZdQBawtxn7GJNI7P2cwCzZm8ZUBl2uC/q7ji9eOwJcq6rrg3d0+xB3AmNxWpCON3C/tdT/OjyG80EDUAQMFJEuqnq4gVjbu/sYY+pn7+cEZs34prXeB+4K9NOJyHj3+m5AuarWAV/HWaGrOdbiNPehqkdxlov8fWCUsYj0FZF/cS8L0Acoad1DMSbh2fs5TlmyN631MyAVWCkiq9y/Af4E3CYii3Ca/I40837fxln6MuA+nGbANe5x3uCLZsEJwCJVrWnRIzDGBNj7OU7Z1DsTtUTkc+BKVT3QxHaPAnNUdW5kIjPGNJe9n/1lZ/Ymms3CGRzUlFX2wWBM1LP3s4/szD7BudNb6ntjXayqNhLWmBhi72fTEEv2xhhjTJyzZnxjjDEmzlmyN8YYY+KcJXtjjDEmzlmyN8YYY+KcJXtjjDEmzv1/NB8LJOuam84AAAAASUVORK5CYII=\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "key_list = [\"djf\", \"mam\", \"jja\", \"son\"]\n", "df_list = [season_dict[key] for key in key_list]\n", "#\n", "# set up a 2 row, 2 column grid of figures, and\n", "# turn the ax_array 2-d array into a list\n", "#\n", "fig, ax_array = plt.subplots(2, 2, figsize=(8, 8))\n", "ax_list = ax_array.flatten()\n", "#\n", "# here is the variable we want to histogram\n", "#\n", "var = \"T_mean (C)\"\n", "#\n", "# store the fit parameters in a dictionary\n", "# called temp_params for future reference\n", "#\n", "temp_params = dict()\n", "#\n", "# loop over the four seasons in key_list\n", "# getting the figure axis and the dataframe\n", "# for that season. \n", "#\n", "for index, key in enumerate(key_list):\n", " the_ax = ax_list[index]\n", " the_df = df_list[index]\n", " the_temps = the_df[var]\n", " #\n", " # find the mean and standard deviation to\n", " # set the pdf mu, sigma parameters\n", " #\n", " mu, sigma = st.norm.fit(the_temps)\n", " #\n", " # save the parameter pair for that season\n", " #\n", " temp_params[key] = (mu, sigma)\n", " #\n", " # get the histogram values, the bin edges and the plotting\n", " # bars (patches) for the histogram, specifying 20 bins\n", " # note that we can specify the histogram variable by it's dataframe\n", " # column name when we pass the dataframe to matplotlib via\n", " # the data keyword\n", " #\n", " vals, bins, patches = the_ax.hist(var, bins=20, data=the_df, density=True)\n", " #\n", " # turn the bin edges into bin centers by averaging the right + left\n", " # edge of every bin\n", " #\n", " bin_centers = (bins[1:] + bins[0:-1]) / 2.0\n", " #\n", " # plot the pdf function on top of the histogram as a red line\n", " #\n", " the_ax.plot(bin_centers, st.norm.pdf(bin_centers, mu, sigma), \"r-\")\n", " the_ax.grid(True)\n", " the_ax.set(title=key)\n", " #\n", " # label the xaxis for the bottom two plots (2 and 3)\n", " # and the yaxis for the left plots (0 and 2)\n", " #\n", " if index > 1:\n", " the_ax.set(xlabel=var)\n", " if index in [0,2]:\n", " the_ax.set(ylabel='relative frequency')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Daily average total precipitation\n", "\n", "Precipitation data is a different story. Because negative precipitation\n", "is impossible, there is no way a normal distribution is going to work\n", "to represent the variablity. Instead we use\n", "[scipy.stats.expon.fit](https://docs.scipy.org/doc/scipy/reference/generated/scipy.stats.expon.html#scipy.stats.expon) to fit an exponential distribution.\n", "\n", "\\begin{align*}\n", "f(y) &= \\exp(-y)\\\\\n", "y &= (x - loc) / scale\n", "\\end{align*}\n", "\n", "\n", "where `loc` is the minimum of the data and `scale` is the distance\n", "between the minimum and the mean.\n", "\n", "The fits are not quite as good -- but do capture the one-sided nature\n", "of the variability. No comments on this one, but it's a copy and paste of\n", "the temperature grid plot above." ] }, { "cell_type": "code", "execution_count": 16, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAfsAAAHwCAYAAAChTMYRAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAIABJREFUeJzs3XuYHGWZ///3J5NkkszkSMIQCBCQiISzJOBKgAmCgrhEVxAUfwsuK+sK6wldcZdFxMMXDwgsooLCKggisApxQVlABkFFEg4SwmEJGCEgJCEhYXKe5P79UTWh05nM1Bxqqrvn87quvrq7uqr6np6pvqeeep77UURgZmZmtWtQ0QGYmZlZvpzszczMapyTvZmZWY1zsjczM6txTvZmZmY1zsnezMysxjnZWy4k/UjSVyQdJunpkuV7SnpE0uuSPlFkjGZmA8XgogOw2hYR9wF7liz6V6AlIg4sKCQzswHHZ/bW33YF5hcdhJnZQOJkb31C0oGSHk6b538GDEuXN0talD7+DTAT+I6kVklvLjBkswFJ0kJJn5P0mKRVkq6S1CTpV+nxe5eksem6N0l6WdIKSb+VtHfJfn4k6bvpdq2SfidpB0mXSFou6SlJbsGrEE721muShgK3ANcC44CbgPeXrxcRRwL3AWdFRGNE/F+/Bmpm7d4PHA28Gfhb4FfAvwHjSfJCe3+aXwFTgO2Bh4HryvbzAeDcdLt1wB/S9cYDNwPfzvOHsOyc7K0vvA0YAlwSERsi4mZgTsExmdm2XRYRr0TEiyT/gP8xIh6JiHXAL4ADASLi6oh4PV1+PrC/pNEl+/lFRDwUEWvT7dZGxDURsRH4Wft+rHhO9tYXdgRejC1nVfpLUcGYWZdeKXm8poPnjZLqJF0o6VlJK4GF6evju7OfvgvZesPJ3vrCX4GdJKlk2S5FBWNmfeJDwCzgKGA0MDldrm1tYJXLyd76wh+ANuATkgZL+jvg4IJjMrPeGUlyHf5VYATwtWLDsd5wsrdei4j1wN8BpwHLgZOAnxcZk5n12jUkl+NeBJ4AHig2HOsNbXmZ1axvSToS+GFE7F50LGZmA5XP7C1v+wB/LjoIM7OBzOVyLTeSLgWOB04tOhYzs4HMzfhmZmY1zs34ZmZmNc7J3szMrMbVzDX78ePHx+TJk7tcb9WqVTQ0NOQfUC84xr5RDTFC/8f50EMPLY2ICf32hj2Q5Xiu1N9vJcblmLKrxLg6iynz8RwRNXE76KCDIot77rkn03pFcox9oxpijOj/OIG50YNjDDgGeBpYAJzTwesfA+YBjwL3A1NLXvtCut3TwLu6eq8sx3Ol/n4rMS7HlF0lxtVZTFmPZzfjm1mXJNUBlwPHAlOBD0qaWrba9RGxb0QcAHyDdMazdL2Tgb1J/mH4bro/M+snTvZmlsXBwIKIeC6Siok3kNRN3ywiVpY8bQDah/rMAm6IiHUR8WeSM3yXUzbrRzVzzd7McrUT8ELJ80XAIeUrSToT+AwwFDiyZNvSUquL0mVm1k+c7M0si45mOtuqSEdEXA5cLulDwLkkBZUybSvpDOAMgKamJlpaWjoNqLW1tct1ilCJcTmm7Coxrr6IycnezLJYBOxc8nwS8FIn698AfK8720bElcCVANOmTYvm5uZOA2ppaaGrdYpQiXE5puwqMa6+iMnX7M0siznAFEm7SRpK0uFudukKkqaUPD0OeCZ9PBs4WVK9pN2AKcCD/RCzmaV8Zm9mXYqINklnAXcAdcDVETFf0gUkQ39mA2dJOgrYQDLV8anptvMl3UgyTWobcGZEbCzkBzEboAZcsp/34gpOO+e2TtdZeOFx/RSNWfWIiNuB28uWnVfy+JOdbPtV4Kt9GY+PZbPs3IxvZmZW45zszczMapyTvZmZWY1zsjczM6txTvZmZmY1zsnezMysxjnZm5mZ1TgnezMzsxqXa7KXdIykpyUtkHROB68fLulhSW2STujg9VGSXpT0nTzjNDMzq2W5JXtJdcDlwLHAVOCDkqaWrfY8cBpw/TZ282Xg3rxiNDMzGwjyPLM/GFgQEc9FxHqSWbBmla4QEQsj4jFgU/nGkg4CmoD/zTFGMzOzmpdnst8JeKHk+aJ0WZckDQIuAj6XQ1xmZmYDSp4T4aiDZZFx248Dt0fEC1JHu0nfQDoDOAOgqamJlpaWLnfcNBzO3ret03Wy7CdPra2thcfQFcfYd6olTjOrXnkm+0XAziXPJwEvZdz2b4DDJH0caASGSmqNiC06+UXElcCVANOmTYvm5uYud3zZdbdy0bzOf+yFp3S9nzy1tLSQ5WcpkmPsO9USp5lVrzyT/RxgiqTdgBeBk4EPZdkwIk5pfyzpNGBaeaI3MzOzbHK7Zh8RbcBZwB3Ak8CNETFf0gWSjgeQNF3SIuBE4ApJ8/OKx8zMbKDK88yeiLgduL1s2Xklj+eQNO93to8fAT/KITwzM7MBwRX0zMzMapyTvZmZWY1zsjczM6txTvZmZmY1zsnezMysxjnZm5mZ1bgBnewPeX4e59xzNUTWKr5mA1OG6ao/I+kJSY9JulvSriWvbZT0aHqb3b+RmxkM8GS/z8sL+NiDP2fUulVFh2JWsTJOV/0ISaXL/YCbgW+UvLYmIg5Ib8f3S9BmtoUBneyXNI4FYMKq5QVHYlbRskxXfU9ErE6fPkAXxbLMrH8N7GTf4GRvlkF3p6s+HfhVyfNhkuZKekDSe/MI0Mw6l2u53Eq3Odm3OtmbdSLzdNWSPgxMA44oWbxLRLwkaXfgN5LmRcSzHWzbrSmrK3W66kqcstgxZVeJcfVFTE72wIRVrxUciVlFyzRdtaSjgH8HjoiIde3LI+Kl9P45SS3AgcBWyb67U1ZX6nTVlThlsWPKrhLj6ouYBnQz/ophjawfNNjN+Gad2zxdtaShJNNVb9GrXtKBwBXA8RGxuGT5WEn16ePxwKHAE/0WuZkBA/zMHoklDWOd7M06ERFtktqnq64Drm6frhqYGxGzgW8CjcBNkgCeT3ve70UyffUmkpOLCyPCyd6snw3sZA8saRzjZG/WhQzTVR+1je1+D+ybb3Rm1pUB3YwP+MzezMxqXq7JPkPVrcMlPSypTdIJJcsPkPQHSfPTilwn5RWjk71Z9Ru0aSO3/vjTnDbXBfrMOpJbss9Ydet54DTg+rLlq4G/j4i9gWOASySNySPOJQ1jGbd6JYM2bcxj92bWDzYNqmOX117mTcsWFR2KWUXK88w+S9WthRHxGLCpbPn/RcQz6eOXgMXAhDyCXNIwlrrYxLg1K/PYvZn1k8WNY9m+dVnRYZhVpDyTfXerbnVI0sHAUDoYl9sXXEXPrDa80rgd27tAllmH8uyNn7nq1jZ3IE0ErgVOjYhNHbzerYpbsHXVrR3qR8Et8LEJr/L8vrsAxVTdKlWJFZzKOca+Uy1xVroljWPZ3c34Zh3KM9lnqrq1LZJGAbcB50bEAx2t092KW7B11a2dXxvPicBvH1/JfytZXkTVrVKVWMGpnGPsO9USZ6Vb3DAuaaGLAHV0rmE2cOXZjN9l1a1tSdf/BXBNRNyUY4wsHZH0+3Mzvll1W9w4lvqNbYxZ+3rRoZhVnNySfUS0Ae1Vt54EbmyvuiXpeABJ0yUtAk4kqbI1P938A8DhwGmSHk1vB+QR55qhw3h96HAne7Mqt7hhHABN7qRntpVcK+hlqLo1hw7mvY6InwA/yTO2Uh5rb1b9FjcmnW23b13G0xMmFxuMWYUZ8BX0wMnerBYsbkzO7N0j32xrTvakyd5fEGZVrb0Zf/tVbsY3K9dlspe0T38EUqQljT6zN6t2a4YOY+XQES6sY9aBLGf235f0oKSP51WytmhLGsYyet0q6tvWFx2KmfXCksZxbqUz60CXyT4iZgCnkIyZnyvpeklH5x5ZP1rSkPwPM37VawVHYma9sbhxrJvxzTqQ6Zp9Wqf+XODzwBHAf0p6StLf5Rlcf3HJXLPasLhhnIfemXUgyzX7/SRdTDJW/kjgbyNir/TxxTnH1y+c7M1qQzIZTlpFz8w2y3Jm/x3gYWD/iDgzIh6GzbPRnZtncP3Fyd6sNixuGMfwtnWMXL+66FDMKkqWojrvBtZExEYASYOAYRGxOiKuzTW6frJsxGg2IXfsMatyr4xsH2u/jNfrGwqOxqxyZDmzvwsYXvJ8RLqsZrTVDWbZiFE+szercksaXFjHrCNZkv2wiGhtf5I+HpFfSMVwFT2z6tdeMneCe+SbbSFLsl8l6a3tTyQdBKzJL6RiONmbVb83SuY62ZuVynLN/lPATZLa56KfCJyUX0jFWNIwht2Wv9T1imZWsV4fOoI1g+s9/M6sTJaiOnOAtwD/DHwc2CsiHso7sP62+czeQ3bMtiLpGElPS1og6ZwOXv+MpCckPSbpbkm7lrx2qqRn0tupOQf6xvA7M9ss6xS304HJ6foHSiIirsktqgIsaRjLsLb1HrJjVkZSHXA5cDSwCJgjaXZEPFGy2iPAtIhYLemfgW8AJ0kaB3wRmAYE8FC6bW7ZeHHDOFfRMyuTpajOtcC3gBkkSX86yYFbU5a0d+zxGYFZuYOBBRHxXESsB24AZpWuEBH3RET7f8oPAJPSx+8C7oyIZWmCvxM4Js9gX2kc5zN7szJZOuhNAw6NiI9HxL+kt09k2XmGpr/DJT0sqU3SCWWv9V/THy6sY9aJnYAXSp4vSpdty+nAr3q4ba8taRzLBF+zN9tClmb8x4EdgL92Z8cZm/6eB04DPlu2bb83/TnZm22TOljWYecWSR8mOW6P6MG2ZwBnADQ1NdHS0tJpUE3D4ex927ZavuefxzDqodV8/s2rutxHHlpbWwt53844puwqMa6+iClLsh8PPCHpQWBd+8KIOL6L7TY3/QFIam/625zsI2Jh+tqmsm03N/2lr7c3/f00Q7w94mRvtk2LSGa9bDcJ2GroiqSjgH8HjoiIdSXbNpdt29LRm0TElcCVANOmTYvm5uaOVtvssutu5aJ5W3+Fvb91PG8HfvrHlfz2Hz7Q6T7y0NLSQlex9zfHlF0lxtUXMWVJ9uf3cN8dNd8d0ottc236WzGskfWDBjvZm21tDjBF0m7Ai8DJwIdKV5B0IHAFcExELC556Q7ga5LGps/fCXwhz2DbC+u4k57ZG7pM9hFxbzqMZkpE3CVpBFCXYd+Zm+96um13m/1g201/AOvGjqF5yLLCm3AqsRmpnGPsO5UeZ0S0STqLJHHXAVdHxHxJFwBzI2I28E2gkaQmB8DzEXF8RCyT9GWSfxgALmhvsctLe2Gdpted7M3adZnsJX2UJKGOA95Ecob9feAdXWyaqemvk22by7ZtKV+pu81+sO2mP4DDho7ltUUrCm/CqcRmpHKOse9UQ5wRcTtwe9my80oeH9XJtlcDV+cX3ZYWN/jM3qxclt74ZwKHAisBIuIZYPsM221u+pM0lKTpb3bGuO4A3ilpbNr89850Wa5cMtes+i0fPor1gwZ7+J1ZiSzJfl06thYASYPJ0BwfEW1Ae9Pfk8CN7U1/ko5P9zVd0iLgROAKSfPTbZcB7U1/c+iHpj9wsjerCe1V9Hxmb7ZZlg5690r6N2C4pKNJSub+MsvOMzT9zeGN4hvl2/Zr0x8kyX671Stg40aoy9Itwcwq0ZKGcS6QZVYiy5n9OcASYB7wTyTJ+9w8gyrK0oYx1MUmWLq06FDMrBeS+vg+szdrl6U3/ibgB+mtprWPtefll6GpqdhgzKzHFjeOY/qiJ7pe0WyAyNIb/890cI0+InbPJaICbZHs99+/2GDMrMcWN4xl3JqVsG4d1NcXHY5Z4bJcsy+d9GYYSWe6cfmEU6wtkr2ZVa32sfa8/DLsumvnK5sNAFnms3+15PZiRFwCHNkPsfW7pQ1jkgdO9mZV7ZX2ZP/Xbk3pYVazsjTjv7Xk6SCSM/2RuUVUoNVDh9M6dDiNTvZmVW2Jk73ZFrI0419U8rgNWAj0/+wS/WRJwxgne7Mqt7jByd6sVJbe+DP7I5BKsaRhLLs52ZtVtVdHjGKjBlHnZG8GZGvG/0xnr0fEt/sunOItaRjra/ZmVW7ToDqWNoyhycneDMhWVGca8M8kE+DsBHwMmEpy3b7mrt072ZvVhsUNY92Mb5bKcs1+PPDWiHgdQNL5wE0R8Y95BlaUJQ1j4bXXYO1aGDas6HDMrIcWN46Dl7JOtGlW27Kc2e8CrC95vh6YnEs0FWDzWPtXXik2EDPrlVcax/nM3iyV5cz+WuBBSb8gqaT3PuCaXKMq0JLGksI6LsZhVrWWNIyDxYuhrQ0GZ/mqM6tdWXrjf1XSr4DD0kUfiYhH8g2rOK6iZ1YbFjeOhYgk4e+4Y9HhmBUqSzM+wAhgZURcCiyStFuOMRVqiavomdWExS6sY7ZZl8le0heBzwNfSBcNAX6SZ1BFenWEk71ZLVjc3krnZG+W6cz+fcDxwCqAiHiJjEPuJB0j6WlJCySd08Hr9ZJ+lr7+R0mT0+VDJP1Y0jxJT0r6Qvm2eWmrGwzjxzvZm1U5n9mbvSFLsl8fEUE6za2khiw7llQHXA4cSzIu/4OSppatdjqwPCL2AC4Gvp4uPxGoj4h9gYOAf2r/R6Bf7LCDk71Zlds8sZWH35llSvY3SroCGCPpo8BdwA8ybHcwsCAinouI9cANwKyydWYBP04f3wy8Q5JI/rFokDQYGE4y3G9lhvfsG072ZlVvQ90QmDDBZ/ZmZJvi9lskifi/gT2B8yLisgz73gl4oeT5onRZh+tERBuwAtgufb9VwF+B54FvRcSyDO/ZN5zszWrDxIlO9mZ0MfQubYq/IyKOAu7s5r7VwbLIuM7BwEZgR2AscJ+kuyLiubL4zgDOAGhqaqKlpaXLoJqGw9n7tnW6zvMLN7DTSy9x3z33gDoKMV+tra2ZfpYiOca+Uw1xSjoGuBSoA34YEReWvX44cAmwH3ByRNxc8tpGYF769PmIOL5/osbJ3izVabKPiI2SVksaHRErurnvRcDOJc8nAeUXz9rXWZQ22Y8GlgEfAn4dERuAxZJ+R1Kjf4tkHxFXAlcCTJs2LZqbm7sM6rLrbuWieZ2XF/iX6dPhZz+j+a1vhdGju9xnX2tpaSHLz1Ikx9h3Kj3Okv43R5Mcs3MkzY6IJ0pWex44DfhsB7tYExEH5B5oRyZOhPnzC3lrs0qS5Zr9WmCepKsk/Wf7LcN2c4ApknaTNBQ4GZhdts5s4NT08QnAb9LOgM8DRyrRALwNeCrLD9QndtghuXdTvhlk6H8TEQsj4jFgUxEBbtPEiclxvKmywjLrb1lqSN6W3rolItoknQXcQdL0d3VEzJd0ATA3ImYDVwHXSlpAckZ/crr55cB/AY+TNPX/V/pF0j9Kk/2ee/bb25pVqI763xzSje2HSZoLtAEXRsQtfRlcpyZOTMrlvvpq0lnPbIDaZrKXdHdEvAOYGhGf78nOI+J24PayZeeVPF5LMsyufLvWjpb3G5/Zm5XK0v+mM7tExEuSdgd+I2leRDy71Zt0sw9Olv438xcsY29gzq23smqPPboRcs9VYh8Mx5RdJcbVFzF1dmY/UdIRwPGSbqDsgI+Ih3v1zpXMyd6sVJb+N9uUFuIiIp6T1AIcCGyV7LvbBydT/5vjj4bzz2f6pEnQT/0iKrEPhmPKrhLj6ouYOjtSzgPOITmwv132WgBH9uqdK9nYsTBkiJO9WWJz/xvgRZLLbR/KsqGkscDqiFgnaTxwKPCN3CItN3Ficu8e+TbAbTPZp0Nnbpb0HxHx5X6MqXiDBkFTk5O9Gdn630iaDvyCZKjs30r6UkTsDewFXCFpE0mH4AvLevHny8neDMg2xe3ASvTtXFjHbLMM/W/mkLQClm/3e2Df3APclmHDYMwYJ3sb8LJOcTvwONmb1QYX1jFzst+mHXaAV14pOgoz6y0ne7NsyV7SDEkfSR9PSDvq1LYddoDFi2HjxqIjMbPecLI36zrZS/oi8HmgfU75IcBP8gyqIuywQ5LoX3216EjMrDd23DGZ5ja6UxrArLZkObN/H3A8ySx07WNmR+YZVEXwWHuz2jBxIqxbB6+9VnQkZoXJkuzXp/XqAyCtVV/7nOzNaoOH35llSvY3SroCGCPpo8BdwA/yDasCONmb1QYne7NM4+y/JeloYCWwJ3BeRHR3bvvq09SU3DvZm1U3J3uzrpO9pE8DNw2IBF+qsTG5OdmbVTcne7NMzfijgDsk3SfpTElNeQdVMVxYx6z6jRwJDQ1O9jagdZnsI6K9xvWZwI7AvZLuyj2ySuBkb1YbJk5Mht+ZDVDdqaC3GHgZeBXYPp9wKoyTvVlt2HFHn9nbgJalqM4/p3NQ3w2MBz4aEfvlHVhFcLI3qw2uomcDXJYz+12BT0XE3hHxxe5MTynpGElPS1og6ZwOXq+X9LP09T9Kmlzy2n6S/iBpvqR5koZlfd8+s8MOsHx5UpDDzKqXk70NcNtM9pJGpQ+/ATwvaVzprasdS6oDLgeOBaYCH5Q0tWy104HlEbEHcDHw9XTbwSQleT+W9hdoBjZ06yfrC+1j7T0hjll1mzgRWluTm9kA1NmZ/fXp/UPA3PT+oZLnXTkYWBARz0XEeuAGYFbZOrOAH6ePbwbeIUnAO4HHIuJPABHxakT0/4w0LqxjVhs8/M4GuG2Os4+I96T3PZ3hbifghZLni4BDtrVORLRJWgFsB7wZCEl3ABOAGyLiG+VvIOkM4AyApqYmWlpaugyqaTicvW9bp+u076dx0SKmAfPuvJNXV6/uct99pbW1NdPPUiTH2HeqJc6qVprsp0wpNhazAmQpqnN3RLyjq2UdbdrBsvJpp7a1zmBgBjAdWA3cLemhiLh7ixUjrgSuBJg2bVo0Nzd3ERJcdt2tXDSv8x974SnpftIvhX0nTIAM++4rLS0tZPlZiuQY+061xFnV2pO9h9/ZANXZNfth6bX58ZLGllyvn0wy3r4ri4CdS55PAsqPtM3rpNfpRwPL0uX3RsTSiFgN3A68NduP1Ie2T0cYuhnfrLrtmH5luRnfBqjOrtn/E8n1+bew5fX6W0k63nVlDjBF0m6ShgInA7PL1pkNnJo+PgH4TTrD3h3AfpJGpP8EHAFkHgXQZ4YMgfHjnezNqt2YMVBf72RvA1Zn1+wvBS6V9C8RcVl3d5xegz+LJHHXAVdHxHxJFwBzI2I2cBVwraQFJGf0J6fbLpf0bZJ/GAK4PSJu624MfcJj7c2qn+ThdzagZZn17jJJ+5AMnxtWsvyaDNveTtIEX7rsvJLHa4ETt7HtT0iG3xXLyd6sNjjZ2wCWpYLeF4HL0ttMknH3x+ccV+VwsjerDU72NoBlqaB3AvAO4OWI+AiwP1Cfa1SVpD3ZR/lAAjOrKk72NoBlSfZrImIT0JZW1VsM7J5vWBVkhx1gzRp4/fWiIzErTIbS14dLelhSm6QTyl47VdIz6e3U8m37zcSJSfnrtWsLC8GsKFmS/VxJY4AfkPTGfxh4MNeoKomr6NkAl7H09fPAabxRebN923HAF0kKah0MfFHS2Lxj7pCH39kAlmU++49HxGsR8X3gaODUtDl/YGhP9osWFRuHWXG6LH0dEQsj4jFgU9m27wLujIhlEbEcuBM4pj+C3opL5toAts3e+JK2WcRG0lsj4uF8Qqow+++f3D/4IBx5ZLGxmBUjS+nr7my7U0crdrf8dXdKXwM0LFrEdODxO+9k6fr1XUfeQ5VY/tgxZVeJcfVFTJ0Nvbuok9cCGBiZb/x42GsvuO8+OGerS5VmA0GW0te93ra75a+7VfoaYOpU+OhH2We77XItf12J5Y8dU3aVGFdfxNRZUZ2ZvdpzLZkxA268ETZtgkFZujmY1ZQspa8727a5bNuWPomqu8aPh8GD3YxvA1KWcfYjJJ0r6cr0+RRJ78k/tApy2GGwYgU8/njRkZgVIUvp6225A3hnOr/GWJLpq+/IKc7ODRoETU1O9jYgZTlN/S9gPfD29Pki4Cu5RVSJZsxI7u+/v9g4zAoQEW1Ae+nrJ4Eb20tfSzoeQNJ0SYtIKmJeIWl+uu0y4Msk/zDMAS5IlxVjxx3hxRcLe3uzonRZLhd4U0ScJOmDABGxRlJH1+Fq1+TJyZfE/ffDxz9edDRm/S5D6es5JE30HW17NXB1rgFmtd9+8POf+5KcDThZ/trXSxpO2qlG0puAdblGVWmkpCn/vvtcSc+smjU3J4V1Hnus6EjM+lWWZP9F4NfAzpKuA+4G/jXXqCrRjBnJWPvnny86EjPrqfYezffcU2gYZv2t02SfNtc/BfwdSXWsnwLTIqIl98gqja/bm1W/SZNgjz2gwsZRm+Wt02QfEQHcEhGvRsRtEfE/EbG0n2KrLPvuC6NGJU35Zla9Zs6Ee++FjRuLjsSs32Rpxn9A0vTcI6l0dXXw9rf7zN6s2s2cmQylffTRoiMx6zdZkv1M4A+SnpX0mKR5kjL1bskwU1a9pJ+lr/9R0uSy13eR1Crps1neL3czZsD8+bCsuJFDZtZL7dft3ZRvA0iWZH8s8CaS8rh/C7wnve9UxpmyTgeWR8QewMXA18tevxj4VYYY+8dhhyX3v/tdsXGYWc9NnAh77ulOejagZJn17i8d3TLsu8uZstLnP04f3wy8o30Mv6T3As8B87P+MLmbPh2GDHFTvlm1mzkTfvtbaOt8Ih2zWpGlqE5PZZkpa/M6EdEmaQWwnaQ1wOdJptTdZhN+d2fJgu7PlFXuwDe/GW67jUeOPbbL9+qpSpx1qZxj7DvVEmdNaW6G738fHnkk+SferMblmeyzzHa1rXW+BFwcEa2dFevr7ixZ0IOZssoddxxcfDHNhxwCw4d3+X49UYmzLpVzjH2nWuKsKaXj7Z3sbQDIs15klpmyNq8jaTAwGlhG0gLwDUkLgU8B/ybprBxjzW7GDNiwAebMKToSM+uppqZkyltft7cBIs9kn2WmrNnAqenjE4DfROKwiJgcEZOBS4CvRcR3cow1u7en8wH5ur1ZdWtuTo7jDRuKjsQsd7kl+ywzZQFXkVyjXwB8BthqeF7F2W472HtvF9cxq3YzZ0JrKzz0UNGRmOUuz2v2WWbKWksyJWZn+zg/l+B6Y8YM+OlPkwpcdXVFR2NmPXHegmM1AAAgAElEQVTEEcn9PffA295WbCxmOfMcjz0xYwasXAmPP150JGbWUxMmwD77uLiODQhO9j3RXlzHTflm1W3mzOS6/fr1RUdilisn+57YZZdk9ix30jOrbjNnwurVHl1jNc/Jviek5Oz+vvsgyksHmFnVOPzw5Hh2U77VOCf7npoxA156CRYuLDoSM+up7baD/fbzeHureU72PTVjRnLvpnyz6jZzZjK51bp1RUdilhsn+57aZx8YPdqd9MyqXXMzrF0LDz5YdCRmuXGy76lBg+DQQ31mb1bt2q/buynfaliuRXWq1eRzbutynYUXHpc05d9+OyxdCuPH90NkZtbnxo6FAw9Mkv1553W9vlkV8pl9b7SPt//d74qNwyxnko6R9LSkBZK2KmstqV7Sz9LX/yhpcrp8sqQ1kh5Nb9/v79gzaW6GP/whac43q0FO9r0xbRoMHeqmfKtpkuqAy4FjganAByVNLVvtdGB5ROwBXAx8veS1ZyPigPT2sX4Jurtmzkw66D3wQNGRmOXCyb43hg1L5sJ2srfadjCwICKei4j1wA3ArLJ1ZgE/Th/fDLxDkvoxxt457LCkH46v21uN8jX73jrsMPjWt5IqXCNGFB2NWR52Al4oeb4IOGRb60REm6QVwHbpa7tJegRYCZwbER0OYZF0BnAGQFNTEy1dFLppGg5n79vW6Tpd7aPUW6dMYdMtt/DozJmZt+lIa2trt963Pzim7Coxrr6Iycm+t2bMgAsvTIbtNDcXHY1ZHjo6Qy8vHbmtdf4K7BIRr0o6CLhF0t4RsXKrlSOuBK4EmDZtWjR3cTxddt2tXDSv86+whad0vo8tHH88XHopzYccAsOHZ9+uTEtLC13F3t8cU3aVGFdfxORm/N56+9uTYTtuyrfatQjYueT5JOClba0jaTAwGlgWEesi4lWAiHgIeBZ4c+4R98TMmcmEOL//fdGRmPU5J/veGjs2KbDj4jpWu+YAUyTtJmkocDIwu2yd2cCp6eMTgN9EREiakHbwQ9LuwBTguX6Ku3tmzIC6OtfJt5qUa7LvxXCdoyU9JGleen9knnH22owZydlAW+fXD82qUUS0AWcBdwBPAjdGxHxJF0g6Pl3tKmA7SQuAzwDtx/vhwGOS/kTSce9jEbGsf3+CjEaOTEbYuJOe1aDcrtmXDNc5mqSJb46k2RHxRMlqm4frSDqZZLjOScBS4G8j4iVJ+5B8yeyUV6y9NmMGfO97MG9eUpzDrMZExO3A7WXLzit5vBY4sYPt/hv479wD7CszZ8JFF8GqVdDQUHQ0Zn0mzzP7Hg/XiYhHIqL9muB8YJik+hxj7Z324jq//nWxcZhZ7zQ3w4YNvm5vNSfP3vi9Ga6ztGSd9wOPRMRWU1J1d6gOZBuuk0X5e+03bRojv/Y1HnzLW9gwdmyv9l2JQz/KOca+Uy1xDgiHHgqDBydN+UcfXXQ0Zn0mz2Tfm+E6yYvS3iRN++/s6A26O1QHsg3XyWKrIT3XXgv77suht90GP/xhr/ZdiUM/yjnGvlMtcVajLPNcQDrXBUBjY1Ioy9ftrcbkmey7M1xnUelwHQBJk4BfAH8fEc/mGGePdPQlcs5bj+djV13FrPVT+dOOe77xBWJm1eNd74IvfSkpnfu2txUdjVmfyPOafW+G64wBbgO+EBFVM8vMZW8/mVcax3HBnd9HsanocMysJz79aZg0CU4/PamXb1YDckv2vRyucxawB/AfJbNlbZ9XrH1lVf0IvjrzH9j/5Wf4wGN3Fh2OmfXEqFHw/e/DE0/A175WdDRmfSLXcrm9GK7zFeArecaWl9l7HcEpj/yKf733x7D8S0nRHTOrLu9+N3z4w0myf//7Yb/9io7IrFdcQa+vSZx/9D8xZm0rnHde1+ubWWW6+OLkn/XTT3fBLKt6TvY5eHL73bn2wHfDd78Lf/pT0eGYWU+MHw+XXQZz58IllxQdjVmvONnn5NuHfRjGjYOzzoIoH3FoZlXhAx+AWbPgP/4DFiwoOhqzHnOyz8nKYY3J1Lf33w/XX190OGbWE1LSQldfD//4j7DJo2ysOjnZ5+kjH0kKdHz2s7Byq+m7zawa7LgjfOtbcO+98IMfFB2NWY/k2ht/wBs0CC6/HA45BL78ZfjmN4uOyMw60GWlvdiBhUceCZ/7HBx3XDIO36yKONnnqP0L5MJ9j+b9376YY5a/iWfH77zVeq60Z1bhpOSsfp994GMfg1/+MllmViXcjN8PvnHEqaweMozz77rCnfXMqtXuu8NXvwq33QY//WnR0Zh1i5N9P1g2YjQXHfZhDvvLo7znqfuKDsfMeuoTn0guy33iE7BkSdHRmGXmZvx+ct2B7+bEeXdx6S+/xV6L/8ylh36I9YOHFB2WmWXUflluyt5/z21zPsmvjjiRTx7/uS3W8SU5q1Q+s+8nGwfV8aEPfo2b93kHZz5wE7/88SfZ96/PFB2WmXXTMxN25TtvP4lZT97LFT//CpNWvFJ0SGZdcrLvR6/XN/D5d3+S0044n9FrW/nFtWdz9m+v9cxaZlXmO3/zAb5+xKkctvAR7vrhP3PW72+gvm190WGZbZOb8QvQ8qZpvPP07/Ifd/+Qf/nDz3hq0gOcfdynmb/DHgCcvW8bp5UNBXLzoFnl2DSoju+97URumdrMv//mKj573084Yd7dnPaXP3HQew/c6vgt5WPZiuAz+4KsHNbI5477FB854YuMXfs6t17zGT59308YsnFD0aGZWUZ/HTWBs957Dqec9BU2DqrjRzd/ieO+/TUmvfZy0aGZbcHJvmD3vGk6R5/+XW7du5lP/v4GZv/40+zwzFMeomdWRX43+QCO+YfL+H/NpzFp/mPcddXH+cTvfuqmfasYbsavACuHNXL2cZ/h9j0P5f/9+juc+KVzOKJhLL/bdX9+N/kA7t/1gK4rfOHmQbMibagbwhWHnEDT+2aw/fd+zGfuv473P343P9/7SOZOmsojO+7J6qHDMx3L4OPZ+lauyV7SMcClQB3ww4i4sOz1euAa4CDgVeCkiFiYvvYF4HRgI/CJiLgjz1grwd17HMJR/7g33159P6vun8dhCx/hfU+0ALBg3CTun3wAv5t8AA/ssi+v1zcUG6wNOD6es1m13XjOmvV5rt//Xfzrb6/hk7/7KYMI2jSIJ7ffjbmTpjJ3p6nMmTSVxSO3KzpcGyByS/aS6oDLgaOBRcAcSbMj4omS1U4HlkfEHpJOBr4OnCRpKnAysDewI3CXpDdHxMa84q0UK4c18uT0o7ho3DEoNrHnkr9w6MJHmfGXR/nAvDs57eH/YaMGMb9pd14Y3cQrjdvx8sjteHnkeD5wyuO83Lgdr4zcjnWDh3a4f58tWE/4eO6+308+gPdOPoCR61Zx4ItPcdCLTzJ90ROc9Nj/8pGHfgnAC6ObmDNpKgu225klDWNY2jCWpSPGsLRhDG/+7C19UovDx7xBvmf2BwMLIuI5AEk3ALOA0i+HWcD56eObge9IUrr8hohYB/xZ0oJ0f3/IMd6KExrEU9vvxlPb78ZVB7+PIRs3cOBLT3Powkc56MUnecuSv3DEnx+mcf2arbZdPmwkL4/cjtahI1gzpJ41Q+pZPWQY1//6O5sfty9fXzeEtkF1tA0aTNugQVz64ekwZAgMHpzchgyBurpkYh8JBg1i1Lx5MHToFsuQtrzBtp+XK1+Wpe54F+uM+Mtf4Mknu95P1vfLyYjnn4ennurZxjvvDA390srj47mHXq9v4Le7H8Rvdz8IgMEb25i6+DmmLXqSaYvmM2Pho/zd/Hs63HZFfQNLGsaytGEMrw0fydrBQ1k7uD69H8q6wfWsHTJ08/J1g5NjeaPq2DhoEBsH1fEPJzzIrN3hw9+7n42DkuWbJIL0puS2eVn6PEiOiUgPjV996ojkQemxXHrfLsvxTTePz8723cd6dTz2xvDhsOuuue0+z2S/E/BCyfNFwCHbWici2iStALZLlz9Qtu1O+YVaHTbUDeHBnffhwZ332WJ547rVNL3+KhNfX8oOra/S9PqryX3rMhrWr2bkutVMWLWcERvWMmLDWoZtWMeIDeuoi23Mzf0/Xcfy1j74efJ2cNEBZNSrOO+4A975zr4KpTM+nvtIW91gHpv4Zh6b+Gaunj4LgPoN65iw+jXGr0puE1YtZ/yq5Yxf/cbz3Za9yLC29cltwzrqN25gWDc6AM7qbeD/1dsdbKlSj8+i4pqz01RO/PA3tljWl60yeSb7jv4FK+9ivq11smyLpDOAM9KnrZKezhDXeGBphvUK84kexDg/p1g6UfGfI9URI/Qmzne9qydb9eT0oRKP54r8/fbk+AX4vxxiKVGJn1UlxgRFxfXiE/D192yxSF/f/LCzmDIdz3km+0VA6Xyuk4CXtrHOIkmDgdHAsozbEhFXAld2JyhJcyNiWne26W+OsW9UQ4xQNXFW3PFcqZ9bJcblmLKrxLj6IqY8x9nPAaZI2k3SUJIOOrPL1pkNnJo+PgH4TUREuvxkSfWSdgOmAA/mGKuZdc7Hs1kVy+3MPr1mdxZwB8lQnasjYr6kC4C5ETEbuAq4Nu2ws4zkC4R0vRtJOv+0AWfWes9ds0rm49msuikGWKU2SWekzYUVyzH2jWqIEaonzkpTqZ9bJcblmLKrxLj6IqYBl+zNzMwGGtfGNzMzq3EDJtlLOkbS05IWSDqn6HgAJO0s6R5JT0qaL+mT6fJxku6U9Ex6P7YCYq2T9Iik/0mf7ybpj2mMP0s7bRUd4xhJN0t6Kv1M/6bSPktJn05/149L+qmkYZX4WVa6SjieJV0tabGkx0uWFfr3VqnfKenf+YOS/pTG9aV0eeF/+5X43SZpoaR5kh6VNDdd1qvf4YBI9nqj1OexwFTgg0pKeBatDTg7IvYC3gacmcZ1DnB3REwB7k6fF+2TQGm5q68DF6cxLicplVq0S4FfR8RbgP1J4q2Yz1LSTsAngGkRsQ9JR7f2srKV9llWrAo6nn8EHFO2rOi/t0r9TlkHHBkR+wMHAMdIehuV8bdfqd9tMyPigJIhd737HUZEzd+AvwHuKHn+BeALRcfVQZy3ktQefxqYmC6bCDxdcFyT0j+uI0nq64mkwMPgjj7fgmIcBfyZtB9KyfKK+Sx5o8LcOJKRMP8DvKvSPstKv1XS8QxMBh4veV4xf29pDBX3nQKMAB4mqcBY6N9+pX63AQuB8WXLevU7HBBn9nRc6rOiynVKmgwcCPwRaIqIvwKk99sXFxkAlwD/CrTX190OeC0i2tLnlfB57g4sAf4rbZL7oaQGKuizjIgXgW8BzwN/BVYAD1F5n2Wlq+TjuWL+3irtOyVtLn8UWAzcCTxL8X/7lfrdFsD/SnpISWVJ6OXvcKAk+0zlOosiqRH4b+BTEbGy6HhKSXoPsDgiHipd3MGqRX+eg0lK9n8vIg4EVlEZlz82S6+xzQJ2I5n9rYGkKbpc0Z9lpavEv7+KUonfKRGxMSIOIDmbPhjYq6PV+iueCv9uOzQi3kry/XCmpMN7u8OBkuwzlessgqQhJAfldRHx83TxK5Impq9PJPlPuCiHAsdLWgjcQNLcdQkwRklJVKiMz3MRsCgi/pg+v5kk+VfSZ3kU8OeIWBIRG4CfA2+n8j7LSlexxzMV8PdW6d8pEfEa0ELSp6DIv/2K/W6LiJfS+8XAL0j+OerV73CgJPsspT77nSSRVB17MiK+XfJSadnRU0muuxUiIr4QEZMiYjLJ5/abiDgFuIekJCoUHCNARLwMvCBpz3TRO0gqtlXMZ0nSfP82SSPS3317jBX1WVaBijyeU4X+vVXqd4qkCZLGpI+Hk/zj+yQF/u1X6nebpAZJI9sfA+8EHqe3v8P+7nhQ1A14N8nEUs8C/150PGlMM0iaiB4DHk1v7ya5bnQ38Ex6P67oWNN4m4H/SR/vTlLffAFwE1BfAfEdAMxNP89bgLGV9lkCXwKeSg/ea4H6SvwsK/1WCccz8FOSvhcbSFobTi/6761Sv1OA/YBH0rgeB85Ll1fE334lfbel7/+n9Da//e+7t79DV9AzMzOrcQOlGd/MzGzAcrI3MzOrcU72ZmZmNc7J3szMrMY52ZuZmdU4J3szM7Ma52RvZmZW45zsrc+k81Q3S/o3ST8sOh4zM0u4qI6ZmVmN85m9mZlZjXOytz4jaaGkoySdL+knJctvkvSypBWSfitp7yLjNLMtSfq8pBclvS7paUnvkFQv6RJJL6W3SyTVp+s3S1ok6WxJiyX9VdJHiv45bNuc7K0//AqYAmwPPAxcV2w4ZtYunSnyLGB6RIwE3gUsBP6dZBraA4D9SaZZPbdk0x2A0cBOJJMAXS5pbP9Fbt3hZG+5i4irI+L1iFgHnA/sL2l0wWGZWWIjyeyLUyUNiYiFEfEscApwQUQsjoglJDM2/n8l221IX98QEbcDrcCe5Tu3yuBkb7mSVCfpQknPSlpJcsYAML7AsMwsFRELgE+R/CO+WNINknYEdgT+UrLqX9Jl7V6NiLaS56uBxpzDtR5ysre8fQiYBRxF0uQ3OV2uogIysy1FxPURMQPYFQjg68BL6fN2u6TLrAo52VveRgLrgFeBEcDXig3HzEpJ2lPSkWnnu7XAGpKm/Z8C50qaIGk8cB7wk052ZRXMyd7ydg1J89+LwBPAA8WGY2Zl6oELgaXAyyQdaf8N+AowF3gMmEfSufYrBcVoveSiOtZnJD0PfJikyX5SRPxDwSGZmRk+s7c+ImkCMIHkLH4q8OdiIzIzs3ZO9tZrkqYDzwCXAbcAk4AfFBqUmZlt5mZ8MzOzGuczezMzsxrnZG9mZlbjBhcdQF8ZP358TJ48ucv1Vq1aRUNDQ/4BVTh/DomB+Dk89NBDSyNiQtFxdCbL8Vzpv7tKjw8qP0bH17Wsx3PNJPvJkyczd+7cLtdraWmhubk5/4AqnD+HxED8HCT9peu1ipXleK70312lxweVH6Pj61rW49nN+GZmZjXOyd7MzKzGOdmbmZnVOCd7MzOzGudkb2ZmVuOc7M3MzGqck72ZmVmNc7I3MzOrcTVTVCereS+u4LRzbut0nYUXHtdP0ZhVB0nHAJcCdcAPI+LCstcPBy4B9gNOjoibS17bBfghsDMQwLsjYmFvY/KxbJadz+zNrFOS6oDLgWOBqcAHJU0tW+154DTg+g52cQ3wzYjYCzgYWJxftGbWkQF3Zm9m3XYwsCAingOQdAMwC3iifYX2M3VJm0o3TP8pGBwRd6brtfZTzGZWwmf2ZtaVnYAXSp4vSpdl8WbgNUk/l/SIpG+mLQVm1o98Zm9mXVEHyyLjtoOBw4ADSZr6f0bS3H/VVm8inQGcAdDU1ERLS0unO24aDmfv29bpOl3tI0+tra2Fvn8WlR6j4+s7uSb7DJ16PgacCWwEWoEzIuKJ9LUvAKenr30iIu7IM1Yz26ZFJJ3r2k0CXurGto+UXAK4BXgbHST7iLgSuBJg2rRp0dVsYpdddysXzev8K2zhKZ3vI0+VMCNaVyo9RsfXd3Jrxs/Yqef6iNg3Ig4AvgF8O912KnAysDdwDPBdN/2ZFWYOMEXSbpKGkhybs7ux7VhJ7fNtH0nJtX4z6x95XrPf3KknItYD7Z16NouIlSVPG3ijaXAWcENErIuIPwML0v2ZWT+LiDbgLOAO4EngxoiYL+kCSccDSJouaRFwInCFpPnpthuBzwJ3S5pHckngB0X8HGYDWZ7N+B116jmkfCVJZwKfAYaS/Nffvu0DZdtm7RBkZn0sIm4Hbi9bdl7J4zkkzfsdbXsnyfh7MytInsk+U6eeiLgcuFzSh4BzgVOzbtvdDj1Q+Z16+ks1dSzJkz8HMxsI8kz23e3UcwPwve5s290OPVD5nXr6SzV1LMmTPwczGwjyvGbfZaceSVNKnh4HPJM+ng2cLKle0m7AFODBHGM1MzOrWbmd2UdEm6T2Tj11wNXtnXqAuRExGzhL0lHABmA5SRM+6Xo3kvTabQPOTDv6mJmZWTflOs4+Q6eeT3ay7VeBr+YXnZmZ2cDgcrlmZmY1zsnezMysxjnZm5mZ1TgnezMzsxrnZG9mZlbjnOzNzMxqnJO9mVU9xSYOXfgob3r1ha5XNhuAnOzNrOoF4r9uOp8T5t1ddChmFcnJ3syqn8TyEaMYu2Zl1+uaDUBO9mbWJUnHSHpa0gJJ53Tw+uGSHpbUJumEDl4fJelFSd/JK8Zlw0cxzsnerENO9mbWKUl1wOXAscBU4IOSppat9jxwGnD9NnbzZeDevGIEWDZiFGNXO9mbdcTJ3sy6cjCwICKei4j1JNNRzypdISIWRsRjwKbyjSUdBDQB/5tnkMuHj/aZvdk2ONmbWVd2Akq7uS9Kl3VJ0iDgIuBzOcS1hWXDfc3ebFtynfXOzGqCOlgWGbf9OHB7RLwgdbSbkjeRzgDOAGhqaqKlpaXT9ZuGw9n7tm1+/pb/a2Tso6/z2b3XEYPqALrcR55aW1sLff8sKj1Gx9d3nOzNrCuLgJ1Lnk8CXsq47d8Ah0n6ONAIDJXUGhFbdfKLiCuBKwGmTZsWzc3Nne74sutu5aJ5b3yFnfr6GA6O4Oo5a1k2YjQAC0/pfB95amlpoaufoWiVHqPj6ztO9mbWlTnAFEm7AS8CJwMfyrJhRJzS/ljSacC0jhJ9X1g+fBQAY1ev3JzszSwxoK/Z77H0eU6YdxdE1hZJs4EnItqAs4A7gCeBGyNivqQLJB0PIGm6pEXAicAVkub3d5yvpgl+3JoV/f3WZhVvQJ/Zv+PZB/lCy4+4bc8ZrBk6rOhwzCpWRNwO3F627LySx3NImvc728ePgB/lEB7wxpm9e+SbbW1An9mvqG8EYPTa1oIjMbPeWtae7D3W3mwrAzrZrxyWJvt1TvZm1W758JEAHn5n1oEBnexXDPOZvVmtWDeknlVDhjFuta/Zm5XLNdlnqKf9GUlPSHpM0t2Sdi15baOkR9Pb7Dzic7I3qy3LRoz2mb1ZB3LroFdST/toknG6cyTNjognSlZ7hGQozmpJ/wx8AzgpfW1NRByQV3zgZG9WazwZjlnH8jyzz1JP+56IWJ0+fYAuevP2NSd7s9qyfPgod9Az60Ceyb679bRPB35V8nyYpLmSHpD03jwCfL1+BJsQo5zszWrCshE+szfrSJ7j7DPX05b0YWAacETJ4l0i4iVJuwO/kTQvIp4t265btbRh63ra60eMoLlhJYNKllVLrePeqKaaznny51BblnsyHLMO5ZnsM9XTlnQU8O/AERGxrn15RLyU3j8nqQU4ENgi2Xe3ljZsXU971pBGnntx9RbLiqyn3V+qqaZznvw51JZXR4ymcf0a6tvWs27w0KLDMasYeTbjb66nLWkoST3tLXrVSzoQuAI4PiIWlywfK6k+fTweOBQo7djXZ1YMa2TUulV57NrM+ll7Fb0xPrs320JuyT5LPW3gmyQzYd1UNsRuL2CupD8B9wAXlvXi7zMrhjW6g55ZjVjmkrlmHcq1Nn6GetpHbWO73wP75hlbuxXDGmla+nx/vJWZ5Wz5CJfMNetIl2f2kvbpj0CK4jN7s9rxRn18V9EzK5WlGf/7kh6U9HFJY3KPqJ+trG9wsjerEZvntHczvtkWukz2ETEDOIWkZ/1cSddLOjr3yPrJymGN1G/cQP2GdV2vbDZAZSh9fbikhyW1STqhZPkBkv4gaX5aFvuk8m370mvDR7IJ+Zq9WZlMHfQi4hngXODzJGPh/1PSU5L+Ls/g+oOr6Jl1rqT09bHAVOCDkqaWrfY8cBpwfdny1cDfR8TewDHAJXm2EG4cVMeKYY0+szcrk+Wa/X6SLibpUX8k8LcRsVf6+OKc48udk71Zl7KUvl4YEY8Bm8qW/196stBeO2MxMCHPYJeNGMV27qBntoUsZ/bfAR4G9o+IMyPiYdh84J6bZ3D9YYXntDfrSndLX3dI0sHAUMqKY/W1pIqeO+iZlcoy9O7dJDPQbQSQNAgYFhGrI+LaXKPrBz6zN+tS5tLX29yBNBG4Fjg1IjZtY51ulb8uL33dblzTSEYtXczZ+7YVWgq5GkoxV3qMjq/vZEn2dwFHAe3ZcATwv8Db8wqqPznZm3UpU+nrbZE0CrgNODciHtjWet0tf11e+rrdhLYxNC9bwEXzBhda+roaSjFXeoyOr+9kacYfFhGbM2H6eER+IfWv9mQ/aq1L5pptQ5elr7clXf8XwDURcVOOMW62eTKc6Fbjg1lNy5LsV0l6a/sTSQcBa/ILqX+trG8AfGZvti1ZSl9Lmi5pEXAicIWk+enmHwAOB05LS2I/KumAPONdNnwU9RvbaFxfM19TZr2WpRn/UyS169ub7SYCuY6V7U+bBtWxcugIJ3uzTmQofT2HpHm/fLufAD/JPcAS7SVzPfzO7A1dJvuImCPpLcCeJB11noqIDblH1o9WDmt0b3yzGuGSuWZbyzoRznRgcrr+gZKIiGtyi6qfrRjWyCif2ZvVBJfMNdtal8le0rXAm4BHgY3p4gBqJtmvHOb6+Ga14tURowFPc2tWKsuZ/TRgakTtdm1dMayR3Za9WHQYZtYH3pjm1s34Zu2y9MZ/HNgh70CKtKLe09ya1YrXh45gw6A6n9mblchyZj8eeELSg8DmqeEi4vjcoupnyZz2HmdvVhOkZKy96+ObbZYl2Z+fdxBFWzGskeFt6xjatoH1g4cUHY6Z9dKy4aN8Zm9WIst89vcCC4Eh6eM5JBPj1AyXzDWrLctGjHZvfLMSWaa4/ShwM3BFumgn4JY8g+pvKzeXzHWyN6sFy4aPYpyb8c02y9JB70zgUGAlQDo39fZ5BtXfPM2tWW1ZPsLN+GalsiT7dRGxvv2JpMFknN5S0jGSnpa0QNI5Hbz+GUlPSHpM0t2Sdi157VRJz6S3U7O8X0+t8Jm9WU1ZNnwUY9a8Dhs3dr2y2QCQJdnfK+nfgOGSjgZuAn7Z1UaS6oDLgWOBqWrA460AACAASURBVMAHJU0tW+0RYFpE7EdyqeAb6bbjgC8ChwAHA1+UNDbbj9R9vmZvVluWDx/FIAKWLy86FLOKkCXZnwMsAeYB/0QyGca5GbY7GFgQEc+lLQM3ALNKV4iIeyJidfr0Ad6YSONdwJ0RsSwilgN3AsdkeM8ecbL//9u79zgpqjP/458vMzADjMDAIFFRQCEmXogiEldRUQdXjcFkowvRbDTramIubhKN0U3iLZpofhL1p27iNfH2ixqSVRKN/EQcNYoK4gVRVFTU8YpcBlGBDDz7R1Vj0/b01Mx0dXVXP+/Xq1/TXVWn+umarn6qTp06x7l0WRH2osd77yUbiHNlIspAOBuBq8NHV2wDvJ71upXgTL0jxwN/K1B2m9wCkk4ETgQYNmwYLS0tnQY1rC+csmv7ZtN6tdcBcNjA1QzdtT3SeirdmjVrquJzdsa3QzplBsPxZO9cIErf+K+Q5xq9mW3fWdE80/Je65f0NYJueffvSlkzuwq4CmD8+PE2adKkTkKCy26+g+kLcz92Lcf1rueZVz9k+sJalh7T+XoqXUtLC1G2V9r5dohG0iHApUANcI2ZXZAzfz/gEmAsMM3MZmTNO5aPawPPM7Pr444302Uuy5bF/VbOVYSofeNn1ANHAYMjlGsFts16PRx4M3chSc3AT4D9zWxdVtlJOWVbIrxnt3kves7ll9X+ZjLBvjlP0kwzezZrsdeA44BTc8pm2t+MJzhgfzwsG+vFdD+zd25zUTrVWZ71eMPMLgEOjLDuecAYSaMk9QGmATOzF5C0O8H9+1PM7N2sWbOAgyU1hg3zDg6nxabNx7R3riNR2t8sNbOngY05ZUva/iZjZd8tgiee7J0DolXjj8t62YvgCH2LzsqZWbuk7xIk6RrgOjNbJOlcYL6ZzQT+D9AA/FESwGtmNsXMVkj6OcEBA8C5ZraiKx+sq1bX+2A4znWgq+1vOiv7ifY3xba2dz0f9q6jnyd754Bo1fjTs563E3Sd+69RVm5mdxG03s+edmbW8+YCZa8DrovyPsXQVt/AtqveLtXbOVdJIre/6UnZrja4zdfYNlv7wAG8/cwzLE6oAWYlNP4s9xg9vuKJ0hr/gFIEkrS2ugZ28TN75/KJ1P6mQNlJOWVb8i3Y1Qa3+Rvbfmz/2oGMranhUwk1wKyExp/lHqPHVzxRqvF/WGi+mf26eOEkp62+v1fjO5ffpvY3wBsE7W+Ojlh2FvCLrE6xDgbOKH6In7Sy7wC/Zu9cKEqnOuOBkwius20DfIugR7wtiHDtvlK01TfQ/x9rqd3QcbWgc9XIzNqBTPub54DbMu1vJE0BkLSnpFaCu3WulLQoLLsCyLS/mUcJ2t9krPBk79wmUa7ZNwHjzOx9AElnA380s/+IM7BS29Q//jq//c65XBHa38zj4x4wc8uWtP1Nxsq+A+B1T/bOQbQz++2A9Vmv1wMjY4kmQd5lrnPpsrzfQHj/fVi3rvOFnUu5KGf2NwKPSfofgla0XwZuiDWqBHiydy5dNvWit3w5bL11ssE4l7AorfHPl/Q3YN9w0jfM7Il4wyq91XWe7J1Lk0296C1b5sneVb0o1fgA/YDVZnYp0Bq2yk0VP7N3Ll1Wepe5zm3SabKXdBbwYz6+XaY3cFOcQSVhdaaBnid751LB+8d37mNRzuy/DEwBPgAwszdJ0S13GX5m71y6bLpm78neuUjJfr2ZGWEXl5L6xxtSMtbX9uaj2jpP9s6lhFfjO/exKMn+NklXAoMknQDMBq6ON6xkeC96zqXHhl410NjoY9o7R7TW+BdJmgysBnYEzjSze2KPLAE+zK1zKdPU5Gf2ztFJspdUA8wKR6dLZYLP1lbfwIC13oOec6nhyd45oJNqfDPbAHwoaWCJ4klUm49p71y6eLJ3DojWg95aYKGkewhb5AOY2cmxRZWQ1fUNfPbdpUmH4ZwrlqFDYcGCpKNwLnFRkv2d4SP12uoa/D5759KkqSlooGcGUtLROJeYDpO9pHvN7CBgJzP7cQljSkxbfQMD1n8IGzZATU3S4TjneqqpCdavhzVrYIvUdQ/iXGSFrtlvJWl/YIqk3SWNy36UKsBSWl0fdiGwalWygThXhiQdIul5SUsknZ5nfp2kW8P5j0oaGU7vLel6SQslPSfpjNyysWlqCv76dXtX5QpV458JnE4wRvWvc+YZcGBcQSUl04seK1fCkCHJBuNcGQnvzLkCmAy0AvMkzTSzZ7MWOx5YaWajJU0DLgSmAkcBdWa2q6R+wLOS/mBmS2MPPDvZj0rdkB7ORdZhsjezGcAMST8zs5+XMKbEbJbsnXPZJgBLzOxlAEm3AEcA2cn+CODs8PkM4HJJIjg56C+pFugLrCfotyN+Q4cGf/3M3lW5TnvQ60mij1Dtt5+kBZLaJR2ZM2+DpCfDx8zuxtAVnuyd69A2wOtZr1vDaXmXMbN2oA0YQpD4PwDeAl4DLjKzFXEHDHg1vnOhKK3xuyVitd9rwHHAqXlW8ZGZ7RZXfPm01Xmyd64D+ZqyW8RlJgAbgK2BRuBBSbMztQSbCksnAicCDBs2jJaWloIBDesLp+zaXnCZvy9ezERgydy5tG67bcFli23NmjWdfoaklXuMHl/xxJbsiVDtl7lmJ2ljjHFE5mf2znWoFcjOlsOBNztYpjWssh8IrACOBu42s38A70p6CBgPbJbszewq4CqA8ePH26RJkwoGdNnNdzB9YeGfsKW/PAxqahg9aBCjO1lfsbW0tNDZZ0haucfo8RVPpGQvaSIwxsx+J2ko0GBmr3RSLF+13+e7EFu9pPlAO3CBmd2eJ64unQlA4bOBmvX1ALz8+OO8ViFHa91VSUekcfLtENk8YIykUcAbwDSCJJ5tJnAsMBc4EphjZibpNeBASTcB/YC9gEtKErXkveg5R4RkL+ksgqPwHYHfAb2Bm4B9OiuaZ1putV8h25nZm5K2B+ZIWmhmL222si6eCUBnZwO1HF/Tm+0HD2b7Cjla665KOiKNk2+HaMysXdJ3gVlADXCdmS2SdC4w38xmAtcCN0paQnBGPy0sfgXBb8czBL8LvzOzp0sWvCd75yKd2X8Z2B1YABAm4Ci9U0Sp9uuQmb0Z/n1ZUksYw0sFCxVBW30DW3o1vnOfYGZ3AXflTDsz6/lagtvscsutyTe9ZIYO9WTvql6U8ezXm5kRnpVL6h9x3Zuq/ST1ITjKj9SqXlKjpLrweRNBLcKzhUsVR1t9g1+zdy5NMl3mOlfFoiT72yRdCQySdAIwG7i6s0LhrTeZar/ngNsy1X6SpgBI2lNSK8FR/5WSFoXFPwvMl/QUcB/BNXtP9s65rvNqfOc6r8Y3s4skTSboBGNH4EwzizS2fYRqv3kE1fu55R4Gdo3yHsW2uq6/J3vn0qSpCVas8DEvXFWL0kDvB8Afoyb4Shec2b+adBjOuWJpaoKNG4MxL7wbbFelolTjDwBmSXpQ0nckDYs7qCR5Nb5zKeNd5joXqbvcc8xsZ+A7BD1g3S9pduyRJaStvgHa2oIzAedc5ct0meuN9FwVi3Jmn/Eu8DawHNgynnCSt7q+AcyChO+cq3zeP75znSd7SSeF97nfCzQBJ5jZ2LgDS4p3metcyniydy5SpzojgO+b2ZNxB1MOPNk7lzKe7J3rONlLGmBmq4Ffha8HZ88v2RCVJebJ3rmU6dcP+vb1ZO+qWqEz+/8HHA48TtB7XnZf9wZsH2NcidmU7FetSjYQ51zxeJe5rsp1mOzN7PDw76jShZM8H9PeuRTyLnNdlYvSQO/eKNPSwqvxnUsh7zLXVblC1+zrCcaebpLUyMfV+AMI7rdPpY9610Hv3p7snUuTpiZYsiTpKJxLTKEz+28SXK//TPg387iDYHzqdJJg0CBP9s5lkXSIpOclLZF0ep75dZJuDec/Kmlk1ryxkuZKWiRpYXgiUVp+Zu+qXKFr9pcCl0r6npldVsKYktfY6MneuZCkGoID/MlAKzBP0syckSiPB1aa2WhJ04ALgamSaoGbgH8zs6ckDQH+UeKPEDTQW70a1q+HPn1K/vbOJS3KqHeXSdoF2Amoz5p+Q5yBJcqTvXPZJgBLzOxlAEm3AEcA2cn+CODs8PkM4HJJAg4GnjazpwDMbHmpgt5M9r32W6f2KqRzHYrSQO8s4LLwcQDBffdTYo4rWZ7sncu2DfB61uvWcFreZcysHWgDhgCfBkzSLEkLJJ1Wgng/yTvWcVUuSg96RwKfA54ws2+Eo95dE29YCWtshBdfTDoK58qF8kyziMvUAhOBPYEPgXslPW5m+e7yORE4EWDYsGG0tLQUDGpYXzhl1/aCy2TWMai1ld2AJ2fPZtWK0vQHtmbNmk4/Q9LKPUaPr3iiJPuPzGyjpHZJAwgGxEllhzqb+Jm9c9lagW2zXg8H3uxgmdbwOv1AYEU4/X4zew9A0l3AOIKxNjZjZlcBVwGMHz/eJk2aVDCoy26+g+kLC/+ELT0mXEd4Zr/b8OHQyXqLpaWlhc4+Q9LKPUaPr3iijHo3X9Ig4GqC1vgLgMdijSppjY1BD3o+zK1zAPOAMZJGSeoDTANm5iwzEzg2fH4kMMfMDJgFjJXULzwI2J/Nr/WXhlfjuyoXpYHet8Onv5V0NzDAzJ6ON6yENTYGiX7NGhgwIOlonEuUmbVL+i5B4q4BrjOzRZLOBeab2UzgWuBGSUsIzuinhWVXSvo1wQGDAXeZ2Z0l/xBDhgR/vRc9V6UKdaozrtA8M1sQT0hloLEx+LtypSd75wAzuwu4K2famVnP1wJHdVD2JoLb75LTu3fQf4af2bsqVejMfnqBeQYc2NnKJR0CXEpwNnCNmV2QM38/4BJgLDDNzGZkzTsW+Gn48jwzu76z9yua7GQ/YkTJ3tY5FyPvWMdVsUKd6hzQkxVH7IjjNeA44NScsoOBs4DxBAcWj4dlS9NqLjvZO+fSwZO9q2JR7rPvJ+mnkq4KX4+RdHiEdW/qiMPM1gOZjjg2MbOl4fX/3JZw/wzcY2YrwgR/D3BIhPcsjkGDgr+e7J1LD0/2ropFaY3/O2A9sHf4uhU4L0K5KB1xxFG25/zM3rn08THtXRWLcp/9DmY2VdJXAczso7AbzM5E6YijR2W72gkHROuI48FnnmFf4KX583l9hx2ixFtxKqkziDj5dqgimTHtzYIBr5yrIlGS/XpJfQmTraQdgHURykXpiKNQ2Uk5ZVtyF+pqJxwQsSOOXx4GNTXsMHgwO1RIhwldVUmdQcTJt0MVaWqCdevggw+goSHpaJwrqSjV+GcBdwPbSrqZoOerKP1bR+mIoyOzgIMlNUpqJBhMY1bEsj3nw9w6lz7esY6rYgWTfVhdvxj4F4JW838AxptZS2crDgfDyHTE8RxwW6YjDklTwvXvKamV4P7cKyUtCsuuAH5OcMAwDzg3nFY63mWuc+niyd5VsYL12WZmkm43sz2ALvd6FaEjjnkEVfT5yl4HXNfV9ywaT/bOpcvQocFfT/auCkWpxn9E0p6xR1JuPNk7ly6ZM3vvMtdVoSgN9A4AvinpVeADgpbyZmZjY40saY2NsHRp0lE454rFq/FdFYuS7A+NPYpy5Gf2zqXLwIFQU+PJ3lWlKKPevVqKQMpOJtn7PbnOpUOvXsHod57sXRWKcs2+OjU2Qnt7cE+ucy4dvMtcV6U82XfE+8d3bjOSDpH0vKQlkk7PM79O0q3h/EcljcyZv52kNZJOzS1bMlttBa+8ktjbO5cUT/Yd8f7xndskaxTLQ4GdgK9K2ilnseOBlWY2GrgYuDBn/sXA3+KOtaB994Unn4QVpe22w7mkebLviCd757J1Oopl+Pr68PkM4KDMOBqSvgS8DCwqUbz5NTcH7XDuuy/RMJwrNU/2HfFk71y2KCNRblom7EGzDRgiqT/wY+CcEsRZ2IQJQb/4s2cnHYlzJRXl1rvq5MneuWxRRqLsaJlzgIvNbE2hATO7OopllBEs861jl113pd9f/sJjU6cWLNtTlTCiYrnH6PEVjyf7jniydy5blFEsM8u0SqoFBgIrgM8DR0r6FTAI2ChprZldnl24q6NYRhrB8pg865g6Fb7/fSaNHAkjRxYs3xOVMKJiucfo8RWPV+N3ZMCA4P56T/bOQbRRLGcCx4bPjwTmWGBfMxtpZiOBS4Bf5Cb6kmpuDv7ee29iIThXap7sO9Krlw9z61woyiiWwLUE1+iXAD8EPnF7XlnYaSf41Kf8ur2rKl6NX0hjI6xalXQUzpWFCKNYriUYrrrQOs6OJbiukIKz+1mzYOPG4MDeuZTzb3kh3j++c+nU3ByMfrdwYdKROFcSnuwL8WTvXDoddFDw16vyXZXwZF+IJ3vn0mn4cPjMZzzZu6rhyb4Qb6DnXHo1N8MDD8C6dUlH4lzsPNkXkj3MrXMuXZqb4cMP4ZFHko7Eudh5si+ksRHWr4ePPko6EudcsU2aFLTE96p8VwU82Rfiveg5l14DBwZ95Xuyd1Ug1mTf3fGvJY2U9JGkJ8PHb+OMs0Oe7J1Lt+ZmeOwxaGtLOhLnYhVbsi/C+Ncvmdlu4eNbccVZkCd759KtuTnoWKdCBjNxrrviPLPv0fjXZcGTvXPpttde0K+fV+W71Isz2Xd7/Otw3ihJT0i6X9K+McbZMU/2zqVbXR3st58ne5d6cfaN35Pxr98CtjOz5ZL2AG6XtLOZrd6scBfHv4aujYFdu3o1E4EX583jjREjOl13JamkcZjj5NvB0dwMp54Kra1BZzvOpVCcyb7b41+bmQHrAMzscUkvAZ8G5mcX7ur419DFMbA3bABgTFMTYypkzOKoKmkc5jj5dnCbDXl77LGFl3WuQsVZjd/t8a8lDQ0b+CFpe2AM8HKMseZXUxPcnuPV+M6l1667wtChXpXvUi22M3sza5eUGf+6BrguM/41MN/MZhKMf31jOP71CoIDAoD9gHMltQMbgG+Z2Yq4Yi3Iu8x1DkmHAJcS7MvXmNkFOfPrgBuAPYDlwFQzWyppMnAB0AdYD/zIzOaUKu6Rp98ZabmlBx0UJHuzYAhc51Im1vHsuzv+tZn9CfhTnLFF5oPhuCqXdRvtZIJLb/MkzTSzZ7MW23QbraRpBLfRTgXeA75oZm9K2oXg4D+3oW7ympvhllvg2Wdh552Tjsa5ovMe9Drjyd65bt9Ga2ZPmFmmrc4ioD6sBSgvmev2XpXvUsqTfWc82TvX09toM74CPGFm5TfM3IgRMHq0J3uXWrFW46eCJ3vnenIbbTBT2pmgav/gDt+ki7fSRrmNNqqWlhbG7LQTw+65h4dmz8Zqe/7TWAm3dZZ7jB5f8Xiy74wne+e6fRstgKThwP8AXzezlzp6k67eShvlNtqolh4zCZYvh5kz2b9vX9hnnx6vsxJu6yz3GD2+4vFq/M40NsLatcHDuerUk9toBwF3AmeY2UMli7g7DjggaInvVfkuhTzZd8a7zHVVLrwGn7mN9jngtsxttJKmhItdCwwJb6P9IZAZ5fK7wGjgZ1mjWG5Z4o8QzeDBsMcenuxdKnk1fmcyyX7VKthqq2RjcS4hPbiN9jzgvNgDLJbmZrjoInj/fdhii6Sjca5o/My+M5lk/847ycbhnItfczO0t8MDDyQdiXNF5cm+M3vsAQMGwC9+EfSu5ZxLr332gfp6r8p3qePJvjNNTXD++XDPPXDbbUlH45yLU309TJwId98NGzcmHY1zRePJPoqTTgrO8H/wA1i9uvPlnXOV6+ijYfFiOO20pCNxrmi8gV4UNTXw29/ChAnws5/BpZcmHZFzrog2GzDHhnL2uMM5bvp0fv5EG9fu+SUAll7whYSic67n/Mw+qvHjgzP8yy+HBQuSjsY5FxeJcw86gTt33IefzbmGLz57f9IROddjnuy74vzzg3GvTzoJNmxIOhrnXEw29qrhh4efwqPb7sL0Oy/mn159KumQnOsRT/ZdMWgQTJ8Ojz0GV1+ddDTOuRitq+3DCf/yU14ZvDVX/fk8eMoTvqtcnuy76uij4cAD4Ywz/N5751JudX0Dxx51Lu/X9YdDD4WlS5MOyblu8QZ6eWzWWCePHbY/insffBB+9CO44YYSReWcS8LbA5o49qhzmHHzaby3x7585Wu/YlXfAXmX9UZ8rlz5mX03vDRk2+C2nBtvhAoZ3tA5130vDh3B8UeeyfC2d7huxjnU/8MHxnKVxZN9d/3kJzBqFHz727B+fdLROOdiNn/4zpz8xR+x25svcNnMX1Gz0Rvpusrh1fjdNPKcOUza4+v8fsY5/GryCfz3P/3rJ5bxKj3n0mXWjntz5uRvcd49v+H8WVdwdvOJrO1dn3RYznXKk30PtOywJ3/79N6c/PAtzPzsfrQO+lTSITnnYnbTuC8wbM0Kvjf3VqY8dz9zdpjAXTvuw33bj9+svc8pu7ZzXJ72P34S4JIQa7KXdAhwKVADXGNmF+TMrwNuAPYAlgNTzWxpOO8M4HhgA3Cymc2KM9buOvegE9jvlQWcM/tKjv/KmSAlHZJzsaiG/Tmq6ft+jYdHfI7DFz/AP78wl8MXP8iHveu4b/vx3LXjRObssCcd/bx21gA4ww8KXDHFluwl1QBXAJOBVmCepJlm9mzWYscDK81stKRpwIXAVEk7AdOAnYGtgdmSPm1mZXeR7K0BQ7l44tH89L7r+O/bf8nCrcaweOhInh86gpE//mvRkr/v+C5J1bI/RyYxd8RY5o4Yy5mTT2LC689w2PMPccjzD/OF5x/io9o63th9HC98aiKLh47k3YbBtNU3+MmAS0ycZ/YTgCVm9jKApFuAI4DsH4cjgLPD5zOAyyUpnH6Lma0DXpG0JFzf3Bjj7bbf7zGF0ctbmbj0CQ574eFN01fX9ef5phE8PzTzGMk7DYNp71XLP3rV0F5TS3uvmvBRyz9qajDlbzMZ5WzADwhcjKpmf+6qDb1qmDvic8wd8TnOav4mE1oXcdjih/jKCw9x2byPP+K6mt682zCYd/s38m7DYN5pGMy7DYNZ1r+R1XX9WVfbm/U1vVlX24d1tX04+PgrWFfbZ9O09l41bFQvNqgXGyVMvXjhl1+AXr2Chx9IuALiTPbbAK9nvW4FPt/RMmbWLqkNGBJOfySn7Dbxhdoz7TW1nH7oyQBsse4DPr3sVT6zbCk7LnuVHZct5YvPPcDXnvwg0ro2hDuzbdpvhYU7sRH+DedlXmd88OtPrs86+AHYu5exZmP8Pw4NfWpif4+emLhhQzDQUaW6/XY46KBSvFPV7M89sbFXDY9sN5ZHthvL6p3/ndn3LGXbVW+z5QcrGbpmBVt+sJJha5Yzevnr7P3qUwxcF+13oUPTs947/K0w2Ow3I/u3hHBe5rdjn17GRxs3/13JyP196Ujub0wx9/ly3z+LGt9eewVDqcckzmSf75tiEZeJUhZJJwInhi/XSHo+QlxNwHsRluu2Z3pS2DYGj/jFvh0AKP+7EkuzHeLS3NydUiO6UaYc9+ey/t99r+TxGZhtehpR8WMs7j5f1v9jihnf7NndrZ2JtD/HmexbgW2zXg8H3uxgmVZJtcBAYEXEspjZVcBVXQlK0nwzG9+VMmnk2yHg2yGystufy/1/V+7xQfnH6PEVT5yd6swDxkgaJakPQQOdmTnLzASODZ8fCcwxMwunT5NUJ2kUMAZ4LMZYnXOF+f7sXAWL7cw+vGb3XWAWwa0615nZIknnAvPNbCZwLXBj2GBnBcEPCOFytxE0/mkHvlPRLXedq3C+PztX2WQW/eJOGkg6MawurGq+HQK+HSpXuf/vyj0+KP8YPb7iqbpk75xzzlUbHwjHOeecS7mqSfaSDpH0vKQlkk5POp5SkbStpPskPSdpkaT/DKcPlnSPpBfDv41Jx1oKkmokPSHpr+HrUZIeDbfDrWHjM1fmym1/rpT9rJy//5IGSZohaXG4Hf+pDLffD8L/7zOS/iCpvpy2YSFVkeyzuvo8FNgJ+GrYhWc1aAdOMbPPAnsB3wk/++nAvWY2Brg3fF0N/hN4Luv1hcDF4XZYSdDlqytjZbo/V8p+Vs7f/0uBu83sM8DnCOIsm+0naRvgZGC8me1C0FA10y10uWzDDlVFsierq08zWw9kuvpMPTN7y8wWhM/fJ9iBtiH4/NeHi10PfCmZCEtH0nDgC8A14WsBBxJ07QpVsh1SoOz250rYz8r5+y9pALAfwR0dmNl6M1tFGW2/UC3QN+xHoh/wFmWyDTtTLck+X1efqeyusxBJI4HdgUeBYWb2FgQ/VMCWyUVWMpcApwGZLgqHAKvMrD18XZXfiwpU1vtzGe9n5fz93x5YBvwuvMxwjaT+lNH2M7M3gIuA1wiSfBvwOOWzDQuqlmQfqbvONJPUAPwJ+L6ZrU46nlKTdDjwrpk9nj05z6JV9b2oUGX7fyvX/awCvv+1wDjgN2a2O/AByV/y2EzYXuAIYBTB6I39CS4l5SqL72KuWMezLyORuutMK0m9CX6AbjazP4eT35G0lZm9JWkr4N3kIiyJfYApkg4D6oEBBGc6gyTVhkfmVfW9qGBluT+X+X5W7t//VqDVzB4NX88gSPblsv0AmoFXzGwZgKQ/A3tTPtuwoGo5s4/S1WcqhdflrgWeM7PscfGyuzY9Frij1LGVkpmdYWbDzWwkwf9/jpkdA9xH0LUrVMF2SImy25/LfT8r9++/mb0NvC5px3DSQQQ9LpbF9gu9BuwlqV/4/87EWBbbsDNV06lOeER7CR939Xl+wiGVhKSJwIPAQj6+VvdfBNcTbwO2I/gSH2VmKxIJssQkTQJONbPDJW1P0MBrMPAE8LVw3HVXxsptf66k/axcv/+SdiNoPNgHeBn4BsEJadlsP0nnAFMJ7r54AvgPgmv0ZbENC6maZO+cc85Vq2qpxnfOOeeqlid755xzLuU82TvnnHMp58neOeecSzlP9s4551zKebKvIJKGSHoyfLwt6Y2s158YaSkc2aXUAQAAA9dJREFUMepbEdZbK2lVB9M3hOt/JhzRqW8RPsfnJV3cxTL9JbVIKvp3VtK9kgYWe73OFeL7s+/PpeS33lUoSWcDa8zsogLLjAZmmNlunayrFnjPzAZ1ND3sROIW4CEz+79Zy4jge7SRGCkYMrTdzK6IYd3HA01mdmGx1+1cFL4/F3Xdvj/n4Wf2KSHptPBo/RlJ3wsnXwDsGB7JXyBpgKQ5khZIejrsLzsSC44KHwRGSxodvs9vgQXAVpIOlTQ3XPetCgaxyBz1z5X0lIIxn/tJapZ0ezj/PEnXKxgL/EVJ/95BCMcQ9kwVlr9PwdjXL4br+LqkeeHnGhkud5OkK8JlX5K0X/heiyVdm7XuO4Cjo24L5+Lm+7Pvz0VnZv6owAdwNkEvWBAM+fkUwZCLWxAMrzkWGA08mVWmN7BF+HxL4MXweS3ByE2577Fpelj2r8AJ4Xo3Antmret+oF/4+icEvYfVA68A48LpAwl6PGsGbg+nnUfwA1MfrqeVYKSr7DjqgTezXjcDK4Bh4by3gTPDeacAF4XPbwJuCp9/hWCUqp0IDnKfBHbJWufLwKCk/6/+qM6H78++P8f98DP7dNgX+JOZfWjBWNq3AxPzLCfgQklPA/8f2FZSUyfr3kLSkwT9kb8E/D6c/pKZzQuf702w0z0cLnsMMBL4LPCafTzOd5uZbcjzHreb2Vozexd4ANgzZ/6WBD8G2R41s3fMbC3Bjj0rnL4wfO+Mv2RNf9PMnrWgivLZnOWWAVt1tBGcKyHfn31/LrpqGfUu7fINVZnP1wmOxseZWbukVoIj6ULet5xrhMFlPT7Ief+7zezfcpYbR7ThHnOXyX39UZ44s/ue3pj1eiObf6/X5Vkm33L14fs4lzTfn31/Ljo/s0+HB4AvS+qrYDztIwiux71PUA2YMZBgTOt2SZMJBnAohoeB/RUMqpFpaTsGWASMCH8kCK8x1uQp/yVJdeFZyb7A/OyZFgwpWa88LZSLQUGL4Cbg9TjW71wX+f7cA74/5+dn9ilgZo9J+gNB1RzAb8xsIYCk+ZIWAncCvwb+Imk+wXW1F4v0/u8oaAF7a9YO/F9m9qKkrwK/kZQ50j4wzyrmAX8jGKP8LDN7J88y9xJUL7YUI+YcE4C/d1Al6VxJ+f7cY74/5+G33rlESTqP4HagSzpZbk/g22b2jRhiuAK4zczuL/a6nasmvj+XL6/GdxUhbDz0d8XQCQfwhP8wOFc6vj+Xnp/ZO+eccynnZ/bOOedcynmyd84551LOk71zzjmXcp7snXPOuZTzZO+cc86lnCd755xzLuX+F8p0FIpNUXELAAAAAElFTkSuQmCC\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "fig, ax_array = plt.subplots(2, 2, figsize=(8, 8))\n", "ax_list = ax_array.flatten()\n", "var = \"Total Precip (mm)\"\n", "precip_params = dict()\n", "for index, key in enumerate(key_list):\n", " the_ax = ax_list[index]\n", " the_df = df_list[index]\n", " the_precip = the_df[var]\n", " loc, scale = st.expon.fit(the_precip)\n", " precip_params[key] = (loc, scale)\n", " vals, bins, patches = the_ax.hist(var, bins=20, data=the_df, density=True)\n", " bin_centers = (bins[1:] + bins[0:-1]) / 2.0\n", " the_ax.plot(bin_centers, st.expon.pdf(bin_centers, loc, scale), \"r-\")\n", " the_ax.grid(True)\n", " the_ax.set(title=key)\n", " #\n", " # label the xaxis for the bottom two plots (2 and 3)\n", " # and the yaxis for the left plots (0 and 2)\n", " #\n", " if index > 1:\n", " the_ax.set(xlabel=var)\n", " if index in [0,2]:\n", " the_ax.set(ylabel='relative frequency') " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Saving the fit parameters\n", "\n", "We have two new dictionaries: `temp_params` and `precip_params` that\n", "we should save for the future." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] }, { "cell_type": "code", "execution_count": 17, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "wrote the fit coefficients to /Users/phil/repos/eosc213_students/notebooks/pandas/data/processed/fit_metadata.json\n" ] } ], "source": [ "data_dict = dict()\n", "data_dict[\n", " \"metadata\"\n", "] = \"\"\"\n", " loc,scale tuples for daily average temperature (deg C)\n", " and precipitation (mm) produced by 11-pandas4 for YVR\n", " \"\"\"\n", "data_dict[\"temp\"] = temp_params\n", "data_dict[\"precip\"] = precip_params\n", "fit_file = context.processed_dir / \"fit_metadata.json\"\n", "with open(fit_file, \"w\") as f:\n", " json.dump(data_dict, f, indent=4)\n", " \n", "print(f\"wrote the fit coefficients to {fit_file}\")" ] }, { "cell_type": "code", "execution_count": 18, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "{\n", " \"metadata\": \"\\n loc,scale tuples for daily average temperature (deg C)\\n and precipitation (mm) produced by 11-pandas4 for YVR\\n \",\n", " \"temp\": {\n", " \"djf\": [\n", " 3.84256591465072,\n", " 3.595365405560527\n", " ],\n", " \"mam\": [\n", " 9.356482721671577,\n", " 3.4789191637037504\n", " ],\n", " \"jja\": [\n", " 16.88648212846009,\n", " 2.324995343255562\n", " ],\n", " \"son\": [\n", " 10.308878631550368,\n", " 4.372545963729551\n", " ]\n", " },\n", " \"precip\": {\n", " \"djf\": [\n", " 0.0,\n", " 4.8432998097309055\n", " ],\n", " \"mam\": [\n", " 0.0,\n", " 2.6093758371283147\n", " ],\n", " \"jja\": [\n", " 0.0,\n", " 1.2523918301531844\n", " ],\n", " \"son\": [\n", " 0.0,\n", " 3.770662503393972\n", " ]\n", " }\n", "}\n" ] } ], "source": [ "#\n", "# here's what the file looks like:\n", "#\n", "\n", "with open(fit_file,'r') as f:\n", " print(f.read())" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Grouping by decade\n", "\n", "The two functions below take a season dataframe and\n", "create a new column with the decade that each day belongs\n", "to. For example if a measurement in taken in 1941, then\n", "its decade is 194. We can use that see whether the temperature\n", "at the airport has been increasing since the 1930's and whether\n", "that increase depends on season" ] }, { "cell_type": "code", "execution_count": 19, "metadata": {}, "outputs": [], "source": [ "def find_decade(row):\n", " \"\"\"\n", " given a row from an Environment Canada dataframe\n", " return the first 3 digits of the year as an integer\n", " \n", " i.e. turn \"2010\" into the number 201\n", " \"\"\"\n", " year_string=f\"{row['Year']:4d}\"\n", " return int(year_string[:3])\n", "\n", "\n", "def decadal_groups(season_df):\n", " \"\"\"\n", " given a season dataframe produced by groupby('season'), add\n", " a column called 'decade' with the 3 digit decade\n", " and return a groupby dictionary of decade dataframes for that\n", " season with the decade as key\n", " \"\"\"\n", " #\n", " # add the decade column to the dataframe using apply\n", " #\n", " decade=season_df.apply(find_decade,axis=1)\n", " season_df['decade']=decade\n", " #\n", " # do the groupby and turn into a dictionary\n", " #\n", " decade_groups=season_df.groupby('decade')\n", " decade_dict={key:value for key,value in decade_groups}\n", " return decade_dict" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Finding the median temperature\n", "\n", "Once we have the dictionary that has the decade as key and\n", "the dataframe for that decade as value, we can find the median\n", "value of the daily average temperature for each decade using\n", "the `median_temps` function below." ] }, { "cell_type": "code", "execution_count": 20, "metadata": {}, "outputs": [], "source": [ "def median_temps(decade_dict):\n", " \"\"\"\n", " given a decade_dict produced by the decadal_temp function\n", " return a 2-column numpy array. The first column should be the\n", " integer decade (193,194,etc.) and the second column should be\n", " the median temperature for that decade\n", " \"\"\"\n", " values=list()\n", " for the_decade,the_df in decade_dict.items():\n", " T_median=the_df['T_mean (C)'].median()\n", " values.append((the_decade,T_median))\n", " result=np.array(values)\n", " return result\n", "\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Assignment\n", "\n", "Using the functions above, calculate the median temperature for\n", "each decade and save your (decade, temperature) lists into\n", "a dictionary with a key for each season. Use that dictionary\n", "to make a 2 x 2 plot of temperature vs. decade for each of the\n", "four seasons, using the 4 box plotting code above. (My solution\n", "is 14 lines long)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [] }, { "cell_type": "code", "execution_count": 21, "metadata": { "deletable": false, "nbgrader": { "checksum": "d937e270e59de55d9c0874ef073b9db7", "grade": true, "grade_id": "cell-55c843e9a4c0efe7", "locked": false, "points": 10, "schema_version": 1, "solution": true } }, "outputs": [ { "ename": "NotImplementedError", "evalue": "", "output_type": "error", "traceback": [ "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m", "\u001b[0;31mNotImplementedError\u001b[0m Traceback (most recent call last)", "\u001b[0;32m\u001b[0m in \u001b[0;36m\u001b[0;34m\u001b[0m\n\u001b[1;32m 1\u001b[0m \u001b[0;31m# YOUR CODE HERE\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m----> 2\u001b[0;31m \u001b[0;32mraise\u001b[0m \u001b[0mNotImplementedError\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m", "\u001b[0;31mNotImplementedError\u001b[0m: " ] } ], "source": [ "# YOUR CODE HERE\n", "raise NotImplementedError()" ] } ], "metadata": { "jupytext": { "cell_metadata_filter": "all", "notebook_metadata_filter": "all,-language_info", "text_representation": { "extension": ".py", "format_name": "percent", "format_version": "1.2", "jupytext_version": "1.0.3" } }, "kernelspec": { "display_name": "Python 3", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.6.7" }, "nbsphinx": { "execute": "never" }, "toc": { "base_numbering": 1, "nav_menu": {}, "number_sections": true, "sideBar": true, "skip_h1_title": true, "title_cell": "Table of Contents", "title_sidebar": "Contents", "toc_cell": true, "toc_position": {}, "toc_section_display": true, "toc_window_display": true } }, "nbformat": 4, "nbformat_minor": 2 }