https://github.com/blab/mumps-wa-phylodynamics
Tip revision: 03472023d0a2c4cf89c490b889e31507054a5317 authored by Louise H Moncla on 08 October 2020, 00:05:57 UTC
adding in mtt output files
adding in mtt output files
Tip revision: 0347202
Enumerate-DTA-geographic-transitions.ipynb
{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Parse transitions along DTA trees\n",
"\n",
"August 15, 2019 \n",
"\n",
"To properly deal with the uncertainty of the number of transitions into Washington, I am going to go through the posterior distribution of trees in the .trees output file from running DTA in beast. I will use the counts of the number of transitions to generate a distribution over the number of transitions. "
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {},
"outputs": [
{
"name": "stderr",
"output_type": "stream",
"text": [
"WARNING (theano.configdefaults): install mkl with `conda install mkl-service`: No module named 'mkl'\n"
]
}
],
"source": [
"import sys, subprocess, glob, os, shutil, re, importlib\n",
"from subprocess import call\n",
"import imp\n",
"bt = imp.load_source('baltic', '../../baltic/baltic/baltic.py')\n",
"\n",
"%matplotlib inline\n",
"import matplotlib as mpl\n",
"from matplotlib import pyplot as plt\n",
"import matplotlib.patheffects as path_effects\n",
"import matplotlib.lines as mlines\n",
"from matplotlib.font_manager import FontProperties\n",
"import matplotlib.colors as clr\n",
"import pandas as pd\n",
"import textwrap as textwrap\n",
"from textwrap import wrap\n",
"import pymc3\n",
"import rpy2\n",
"import numpy as np\n",
"from scipy.special import binom\n",
"\n",
"%load_ext rpy2.ipython"
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {},
"outputs": [],
"source": [
"mpl.rcParams['font.serif'] = \"Helvetica\"\n",
"typeface='Helvetica'\n",
"mpl.rcParams['font.weight']=300\n",
"mpl.rcParams['axes.labelweight']=300\n",
"mpl.rcParams['font.family']=typeface\n",
"mpl.rcParams['font.size']=16"
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {},
"outputs": [],
"source": [
"def get_taxa_lines(tree_path): \n",
" \n",
" # write out a temp tree file\n",
" temp_tree = tree_path.replace(\".trees\",\".temp.tree\")\n",
" with open(temp_tree, \"w\") as outfile: \n",
" outfile.write(\"\")\n",
"\n",
" lines_to_write = \"\"\n",
" with open(tree_path, 'rU') as infile:\n",
" for line in infile: ## iterate through each line\n",
" if 'state' not in line.lower(): #going to grab all the interesting stuff in the .trees file prior to the newick tree strings\n",
" lines_to_write = lines_to_write + line\n",
"\n",
" return(lines_to_write)"
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {},
"outputs": [],
"source": [
"def get_burnin_value(tree_path, burnin_percent):\n",
" with open(tree_path, 'rU') as infile:\n",
" numtrees = 0\n",
" for line in infile: ## iterate through each line\n",
" if 'state' in line.lower(): #going to grab all the interesting stuff in the .trees file prior to the newick tree strings\n",
" numtrees += 1\n",
" \n",
" burnin = numtrees * burnin_percent\n",
" return(burnin)"
]
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {},
"outputs": [],
"source": [
"def retrieve_subtrees(tree):\n",
" \n",
" traitName='geo'\n",
"\n",
" tree.root.traits[traitName]='ancestor' ## give root node some trait value that's different from what the actual tree root has, so it registers as a switch\n",
"\n",
" tree_strings={division:[] for division in division_order}\n",
" subtype_trees={division:[] for division in division_order}\n",
"\n",
" for k in sorted(tree.Objects,key=lambda x:x.height):\n",
" kp=k.parent # kp is the parent node of k\n",
"\n",
" ## get current node's (k) and its parent's (kp) trait states\n",
" kloc=k.traits[traitName] # kloc = trait of k; kc = trait of k; they are the same thing\n",
" if traitName in k.parent.traits: # if parent has a trait block, use that trait, else assign to ancestor\n",
" kploc=kp.traits[traitName] # kploc = trait of parental node\n",
" else:\n",
" kploc='ancestor'\n",
"\n",
" ## if states do not match\n",
" if kloc!=kploc: # if node and parental node do not have the same trait\n",
" #N_children=len(k.leaves)\n",
" traverse_condition=lambda w:w.traits[traitName]==kloc # traverse tree for all nodes whose traitname = kc\n",
" #print('subtree resulting from %s>%s switch, traversing within %s')%(kploc,kloc,kloc)\n",
"\n",
" subtree=tree.subtree(k,traverse_condition=traverse_condition) #traitName = traitName ## this function returns a new baltic object that contains a trait-traversed subtree, starting from node k, for as long as the traversal stays within the starting trait value state\n",
"\n",
" if subtree != None:\n",
" subtree.traverse_tree()\n",
" subtree.sortBranches()\n",
" subtype_trees[kloc].append((kploc,subtree))\n",
" \n",
" return(subtype_trees)"
]
},
{
"cell_type": "code",
"execution_count": 6,
"metadata": {},
"outputs": [],
"source": [
"def count_transitions(tree_path,temp_tree):\n",
" \n",
" ## count transitions in the posterior set of trees \n",
" trees_per_file={}\n",
" transitions_dict = {}\n",
" treeCounter = 0\n",
"\n",
" # open .trees file and read in line by line\n",
" with open(tree_path, 'rU') as infile:\n",
" for line in infile: ## iterate through each line\n",
"\n",
" # write an empty tree to read file; this will get filled in as a new nexus file for each tree in the posterior\n",
" with open(temp_tree, \"w\") as outfile:\n",
" outfile.write(\"\")\n",
"\n",
" # if tree line...\n",
" if 'tree STATE_' in line:\n",
" treeCounter += 1\n",
" identifier = \" \".join(line.split(\" \")[0:2])\n",
"\n",
" if treeCounter > burnin:\n",
"\n",
" with open(temp_tree, \"a\") as outfile:\n",
" outfile.write(lines_to_write)\n",
" outfile.write(line)\n",
" tree = bt.loadNexus(temp_tree)\n",
"\n",
" # run this function to retrieve all of the transitions into Washington (i.e., all the Washington subtrees)\n",
" subtype_trees = retrieve_subtrees(tree)\n",
"\n",
" # output the number of washington subtrees to dictionary\n",
" transitions_dict[identifier] = len(subtype_trees['washington'])\n",
"\n",
" with open(dict_outfile, \"a\") as outfile: \n",
" outfile.write(identifier + \"\\t\" + str(len(subtype_trees['washington'])) + \"\\n\")"
]
},
{
"cell_type": "code",
"execution_count": 7,
"metadata": {},
"outputs": [],
"source": [
"division_order = [\"california\",\"new_york\",\"massachusetts\",\"virginia\",\"georgia\",\"manitoba\",\"arkansas\",\"north_dakota\",\n",
" \"louisiana\",\"ontario\",\"texas\",\"indiana\",\"michigan\",\"montana\",\"new_hampshire\",\"ohio\",\n",
" \"washington\",\"illinois\",\"kansas\",\"pennsylvania\",\"new_jersey\",\"british_columbia\",\"iowa\",\"missouri\",\n",
" \"north_carolina\",\"alabama\",\"wisconsin\"]"
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {},
"outputs": [],
"source": [
"# specify tree path\n",
"tree_path = \"/Users/lmoncla/Documents/Mumps/beast/DTA/2020-10-23-skygrid-25-exp-bssvs-prior/mumps-north-american-seqs-DTA-2020-01-16.aligned.no-metadata.trees\"\n",
"#tree_path = \"/Users/lmoncla/Documents/Mumps/beast/DTA/2020-10-23-skygrid-25-exp-bssvs-prior/offset-26/mumps-north-american-seqs-DTA-2020-01-16.aligned.no-metadata.trees\"\n",
"\n",
"# specify burnin percentage\n",
"dict_outfile = tree_path.replace(\".trees\",\".transitions.txt\")\n",
"temp_tree = tree_path.replace(\".trees\",\".temp.tree\")\n",
"burnin_percent = 0.10"
]
},
{
"cell_type": "code",
"execution_count": 9,
"metadata": {},
"outputs": [],
"source": [
"with open(dict_outfile, \"w\") as outfile:\n",
" outfile.write(\"state\\tnumber_transitions\\n\")"
]
},
{
"cell_type": "code",
"execution_count": 10,
"metadata": {
"scrolled": true
},
"outputs": [
{
"name": "stderr",
"output_type": "stream",
"text": [
"/Users/lmoncla/anaconda/envs/LHM-basics/lib/python3.7/site-packages/ipykernel_launcher.py:9: DeprecationWarning: 'U' mode is deprecated\n",
" if __name__ == '__main__':\n",
"/Users/lmoncla/anaconda/envs/LHM-basics/lib/python3.7/site-packages/ipykernel_launcher.py:2: DeprecationWarning: 'U' mode is deprecated\n",
" \n"
]
},
{
"name": "stdout",
"output_type": "stream",
"text": [
"1000.1\n"
]
}
],
"source": [
"# get taxa line and burnin number\n",
"lines_to_write = get_taxa_lines(tree_path)\n",
"burnin = get_burnin_value(tree_path, burnin_percent)\n",
"print(burnin)"
]
},
{
"cell_type": "code",
"execution_count": 11,
"metadata": {},
"outputs": [
{
"name": "stderr",
"output_type": "stream",
"text": [
"/Users/lmoncla/anaconda/envs/LHM-basics/lib/python3.7/site-packages/ipykernel_launcher.py:9: DeprecationWarning: 'U' mode is deprecated\n",
" if __name__ == '__main__':\n"
]
}
],
"source": [
"count_transitions(tree_path,temp_tree)"
]
},
{
"cell_type": "code",
"execution_count": 12,
"metadata": {},
"outputs": [
{
"name": "stderr",
"output_type": "stream",
"text": [
"/Users/lmoncla/anaconda/envs/LHM-basics/lib/python3.7/site-packages/ipykernel_launcher.py:1: FutureWarning: from_csv is deprecated. Please use read_csv(...) instead. Note that some of the default arguments are different, so please refer to the documentation for from_csv when changing your function calls\n",
" \"\"\"Entry point for launching an IPython kernel.\n"
]
},
{
"data": {
"text/html": [
"<div>\n",
"<style scoped>\n",
" .dataframe tbody tr th:only-of-type {\n",
" vertical-align: middle;\n",
" }\n",
"\n",
" .dataframe tbody tr th {\n",
" vertical-align: top;\n",
" }\n",
"\n",
" .dataframe thead th {\n",
" text-align: right;\n",
" }\n",
"</style>\n",
"<table border=\"1\" class=\"dataframe\">\n",
" <thead>\n",
" <tr style=\"text-align: right;\">\n",
" <th></th>\n",
" <th>state</th>\n",
" <th>number_transitions</th>\n",
" </tr>\n",
" </thead>\n",
" <tbody>\n",
" <tr>\n",
" <th>0</th>\n",
" <td>tree STATE_10000000</td>\n",
" <td>14</td>\n",
" </tr>\n",
" <tr>\n",
" <th>1</th>\n",
" <td>tree STATE_10010000</td>\n",
" <td>12</td>\n",
" </tr>\n",
" <tr>\n",
" <th>2</th>\n",
" <td>tree STATE_10020000</td>\n",
" <td>14</td>\n",
" </tr>\n",
" <tr>\n",
" <th>3</th>\n",
" <td>tree STATE_10030000</td>\n",
" <td>14</td>\n",
" </tr>\n",
" <tr>\n",
" <th>4</th>\n",
" <td>tree STATE_10040000</td>\n",
" <td>13</td>\n",
" </tr>\n",
" </tbody>\n",
"</table>\n",
"</div>"
],
"text/plain": [
" state number_transitions\n",
"0 tree STATE_10000000 14\n",
"1 tree STATE_10010000 12\n",
"2 tree STATE_10020000 14\n",
"3 tree STATE_10030000 14\n",
"4 tree STATE_10040000 13"
]
},
"execution_count": 12,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"df = pd.DataFrame.from_csv(tree_path.replace(\".trees\",\".transitions.txt\"), sep=\"\\t\")\n",
"df.reset_index(inplace=True)\n",
"df.head()"
]
},
{
"cell_type": "code",
"execution_count": 15,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"13.507054771692035"
]
},
"execution_count": 15,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"df['number_transitions'].mean()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Read back in the output file to dataframe and plot"
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {},
"outputs": [
{
"name": "stderr",
"output_type": "stream",
"text": [
"/Users/lmoncla/anaconda/envs/LHM-basics/lib/python3.7/site-packages/ipykernel_launcher.py:1: FutureWarning: from_csv is deprecated. Please use read_csv(...) instead. Note that some of the default arguments are different, so please refer to the documentation for from_csv when changing your function calls\n",
" \"\"\"Entry point for launching an IPython kernel.\n"
]
},
{
"data": {
"text/html": [
"<div>\n",
"<style scoped>\n",
" .dataframe tbody tr th:only-of-type {\n",
" vertical-align: middle;\n",
" }\n",
"\n",
" .dataframe tbody tr th {\n",
" vertical-align: top;\n",
" }\n",
"\n",
" .dataframe thead th {\n",
" text-align: right;\n",
" }\n",
"</style>\n",
"<table border=\"1\" class=\"dataframe\">\n",
" <thead>\n",
" <tr style=\"text-align: right;\">\n",
" <th></th>\n",
" <th>state</th>\n",
" <th>number_transitions</th>\n",
" </tr>\n",
" </thead>\n",
" <tbody>\n",
" <tr>\n",
" <th>0</th>\n",
" <td>tree STATE_10000000</td>\n",
" <td>14</td>\n",
" </tr>\n",
" <tr>\n",
" <th>1</th>\n",
" <td>tree STATE_10010000</td>\n",
" <td>12</td>\n",
" </tr>\n",
" <tr>\n",
" <th>2</th>\n",
" <td>tree STATE_10020000</td>\n",
" <td>14</td>\n",
" </tr>\n",
" <tr>\n",
" <th>3</th>\n",
" <td>tree STATE_10030000</td>\n",
" <td>14</td>\n",
" </tr>\n",
" <tr>\n",
" <th>4</th>\n",
" <td>tree STATE_10040000</td>\n",
" <td>13</td>\n",
" </tr>\n",
" </tbody>\n",
"</table>\n",
"</div>"
],
"text/plain": [
" state number_transitions\n",
"0 tree STATE_10000000 14\n",
"1 tree STATE_10010000 12\n",
"2 tree STATE_10020000 14\n",
"3 tree STATE_10030000 14\n",
"4 tree STATE_10040000 13"
]
},
"execution_count": 4,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"skygrid = pd.DataFrame.from_csv(\"../beast/output/no-metadata/mumps-north-american-seqs-DTA-2020-12-10-exp-mean-26.transitions.txt\", sep=\"\\t\")\n",
"skygrid.reset_index(inplace=True)\n",
"skygrid.head()"
]
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"13.327074769470059"
]
},
"execution_count": 5,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"skygrid['number_transitions'].mean()"
]
},
{
"cell_type": "code",
"execution_count": 6,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"0 13\n",
"dtype: int64"
]
},
"execution_count": 6,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"skygrid['number_transitions'].mode()"
]
},
{
"cell_type": "code",
"execution_count": 7,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"array([12, 15])"
]
},
"execution_count": 7,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"x = skygrid['number_transitions']\n",
"hpd = pymc3.stats.hpd(x)\n",
"hpd"
]
},
{
"cell_type": "code",
"execution_count": 8,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"<div>\n",
"<style scoped>\n",
" .dataframe tbody tr th:only-of-type {\n",
" vertical-align: middle;\n",
" }\n",
"\n",
" .dataframe tbody tr th {\n",
" vertical-align: top;\n",
" }\n",
"\n",
" .dataframe thead th {\n",
" text-align: right;\n",
" }\n",
"</style>\n",
"<table border=\"1\" class=\"dataframe\">\n",
" <thead>\n",
" <tr style=\"text-align: right;\">\n",
" <th></th>\n",
" <th>state</th>\n",
" <th>number_transitions</th>\n",
" </tr>\n",
" </thead>\n",
" <tbody>\n",
" <tr>\n",
" <th>0</th>\n",
" <td>tree STATE_10000000</td>\n",
" <td>14</td>\n",
" </tr>\n",
" <tr>\n",
" <th>1</th>\n",
" <td>tree STATE_10010000</td>\n",
" <td>12</td>\n",
" </tr>\n",
" <tr>\n",
" <th>2</th>\n",
" <td>tree STATE_10020000</td>\n",
" <td>14</td>\n",
" </tr>\n",
" <tr>\n",
" <th>3</th>\n",
" <td>tree STATE_10030000</td>\n",
" <td>14</td>\n",
" </tr>\n",
" <tr>\n",
" <th>4</th>\n",
" <td>tree STATE_10040000</td>\n",
" <td>13</td>\n",
" </tr>\n",
" </tbody>\n",
"</table>\n",
"</div>"
],
"text/plain": [
" state number_transitions\n",
"0 tree STATE_10000000 14\n",
"1 tree STATE_10010000 12\n",
"2 tree STATE_10020000 14\n",
"3 tree STATE_10030000 14\n",
"4 tree STATE_10040000 13"
]
},
"execution_count": 8,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"skygrid.head()"
]
},
{
"cell_type": "code",
"execution_count": 18,
"metadata": {},
"outputs": [],
"source": [
"min_value = min(skygrid['number_transitions'])\n",
"max_value = max(skygrid['number_transitions'])\n",
"print(min_value, max_value)"
]
},
{
"cell_type": "code",
"execution_count": 44,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"<div>\n",
"<style scoped>\n",
" .dataframe tbody tr th:only-of-type {\n",
" vertical-align: middle;\n",
" }\n",
"\n",
" .dataframe tbody tr th {\n",
" vertical-align: top;\n",
" }\n",
"\n",
" .dataframe thead th {\n",
" text-align: right;\n",
" }\n",
"</style>\n",
"<table border=\"1\" class=\"dataframe\">\n",
" <thead>\n",
" <tr style=\"text-align: right;\">\n",
" <th></th>\n",
" <th>number_transitions</th>\n",
" <th>number_trees</th>\n",
" <th>proportion_trees</th>\n",
" </tr>\n",
" </thead>\n",
" <tbody>\n",
" <tr>\n",
" <th>0</th>\n",
" <td>13</td>\n",
" <td>4193</td>\n",
" <td>0.465837</td>\n",
" </tr>\n",
" <tr>\n",
" <th>1</th>\n",
" <td>14</td>\n",
" <td>2238</td>\n",
" <td>0.248639</td>\n",
" </tr>\n",
" <tr>\n",
" <th>2</th>\n",
" <td>12</td>\n",
" <td>1247</td>\n",
" <td>0.138540</td>\n",
" </tr>\n",
" <tr>\n",
" <th>3</th>\n",
" <td>15</td>\n",
" <td>822</td>\n",
" <td>0.091323</td>\n",
" </tr>\n",
" <tr>\n",
" <th>4</th>\n",
" <td>11</td>\n",
" <td>225</td>\n",
" <td>0.024997</td>\n",
" </tr>\n",
" <tr>\n",
" <th>5</th>\n",
" <td>16</td>\n",
" <td>199</td>\n",
" <td>0.022109</td>\n",
" </tr>\n",
" <tr>\n",
" <th>6</th>\n",
" <td>17</td>\n",
" <td>48</td>\n",
" <td>0.005333</td>\n",
" </tr>\n",
" <tr>\n",
" <th>7</th>\n",
" <td>10</td>\n",
" <td>21</td>\n",
" <td>0.002333</td>\n",
" </tr>\n",
" <tr>\n",
" <th>8</th>\n",
" <td>18</td>\n",
" <td>6</td>\n",
" <td>0.000667</td>\n",
" </tr>\n",
" <tr>\n",
" <th>9</th>\n",
" <td>20</td>\n",
" <td>1</td>\n",
" <td>0.000111</td>\n",
" </tr>\n",
" <tr>\n",
" <th>10</th>\n",
" <td>9</td>\n",
" <td>1</td>\n",
" <td>0.000111</td>\n",
" </tr>\n",
" </tbody>\n",
"</table>\n",
"</div>"
],
"text/plain": [
" number_transitions number_trees proportion_trees\n",
"0 13 4193 0.465837\n",
"1 14 2238 0.248639\n",
"2 12 1247 0.138540\n",
"3 15 822 0.091323\n",
"4 11 225 0.024997\n",
"5 16 199 0.022109\n",
"6 17 48 0.005333\n",
"7 10 21 0.002333\n",
"8 18 6 0.000667\n",
"9 20 1 0.000111\n",
"10 9 1 0.000111"
]
},
"execution_count": 44,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"df = pd.DataFrame(skygrid[\"number_transitions\"].value_counts())\n",
"df.reset_index(inplace=True)\n",
"df.columns = ['number_transitions','number_trees']\n",
"df['proportion_trees'] = (df['number_trees'])/(sum(df['number_trees']))\n",
"df"
]
},
{
"cell_type": "code",
"execution_count": 46,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"<div>\n",
"<style scoped>\n",
" .dataframe tbody tr th:only-of-type {\n",
" vertical-align: middle;\n",
" }\n",
"\n",
" .dataframe tbody tr th {\n",
" vertical-align: top;\n",
" }\n",
"\n",
" .dataframe thead th {\n",
" text-align: right;\n",
" }\n",
"</style>\n",
"<table border=\"1\" class=\"dataframe\">\n",
" <thead>\n",
" <tr style=\"text-align: right;\">\n",
" <th></th>\n",
" <th>number_transitions</th>\n",
" <th>number_trees</th>\n",
" <th>proportion_trees</th>\n",
" </tr>\n",
" </thead>\n",
" <tbody>\n",
" <tr>\n",
" <th>0</th>\n",
" <td>13</td>\n",
" <td>4193</td>\n",
" <td>0.465837</td>\n",
" </tr>\n",
" <tr>\n",
" <th>1</th>\n",
" <td>14</td>\n",
" <td>2238</td>\n",
" <td>0.248639</td>\n",
" </tr>\n",
" <tr>\n",
" <th>2</th>\n",
" <td>12</td>\n",
" <td>1247</td>\n",
" <td>0.138540</td>\n",
" </tr>\n",
" <tr>\n",
" <th>3</th>\n",
" <td>15</td>\n",
" <td>822</td>\n",
" <td>0.091323</td>\n",
" </tr>\n",
" <tr>\n",
" <th>4</th>\n",
" <td>11</td>\n",
" <td>225</td>\n",
" <td>0.024997</td>\n",
" </tr>\n",
" <tr>\n",
" <th>5</th>\n",
" <td>16</td>\n",
" <td>199</td>\n",
" <td>0.022109</td>\n",
" </tr>\n",
" <tr>\n",
" <th>6</th>\n",
" <td>17</td>\n",
" <td>48</td>\n",
" <td>0.005333</td>\n",
" </tr>\n",
" <tr>\n",
" <th>7</th>\n",
" <td>10</td>\n",
" <td>21</td>\n",
" <td>0.002333</td>\n",
" </tr>\n",
" <tr>\n",
" <th>8</th>\n",
" <td>18</td>\n",
" <td>6</td>\n",
" <td>0.000667</td>\n",
" </tr>\n",
" <tr>\n",
" <th>9</th>\n",
" <td>20</td>\n",
" <td>1</td>\n",
" <td>0.000111</td>\n",
" </tr>\n",
" <tr>\n",
" <th>10</th>\n",
" <td>9</td>\n",
" <td>1</td>\n",
" <td>0.000111</td>\n",
" </tr>\n",
" <tr>\n",
" <th>0</th>\n",
" <td>19</td>\n",
" <td>0</td>\n",
" <td>0.000000</td>\n",
" </tr>\n",
" </tbody>\n",
"</table>\n",
"</div>"
],
"text/plain": [
" number_transitions number_trees proportion_trees\n",
"0 13 4193 0.465837\n",
"1 14 2238 0.248639\n",
"2 12 1247 0.138540\n",
"3 15 822 0.091323\n",
"4 11 225 0.024997\n",
"5 16 199 0.022109\n",
"6 17 48 0.005333\n",
"7 10 21 0.002333\n",
"8 18 6 0.000667\n",
"9 20 1 0.000111\n",
"10 9 1 0.000111\n",
"0 19 0 0.000000"
]
},
"execution_count": 46,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"# finally, add in a count of 0 for any transition value that is not observed \n",
"for i in range(min_value, max_value+1):\n",
" x = df[df['number_transitions'] == i]\n",
" if len(x) == 0:\n",
" new_row = pd.DataFrame.from_dict({'number_transitions':[i],\"number_trees\":0,\"proportion_trees\":0})\n",
" df = df.append(new_row)\n",
" \n",
"df"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": []
},
{
"cell_type": "code",
"execution_count": 50,
"metadata": {},
"outputs": [
{
"name": "stderr",
"output_type": "stream",
"text": [
"/Users/lmoncla/anaconda/envs/LHM-basics/lib/python3.7/site-packages/rpy2/rinterface/__init__.py:146: RRuntimeWarning: Scale for 'x' is already present. Adding another scale for 'x', which will\n",
"replace the existing scale.\n",
"\n",
" warnings.warn(x, RRuntimeWarning)\n"
]
},
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAA+gAAAH0CAYAAACuKActAAAEGWlDQ1BrQ0dDb2xvclNwYWNlR2VuZXJpY1JHQgAAOI2NVV1oHFUUPrtzZyMkzlNsNIV0qD8NJQ2TVjShtLp/3d02bpZJNtoi6GT27s6Yyc44M7v9oU9FUHwx6psUxL+3gCAo9Q/bPrQvlQol2tQgKD60+INQ6Ium65k7M5lpurHeZe58853vnnvuuWfvBei5qliWkRQBFpquLRcy4nOHj4g9K5CEh6AXBqFXUR0rXalMAjZPC3e1W99Dwntf2dXd/p+tt0YdFSBxH2Kz5qgLiI8B8KdVy3YBevqRHz/qWh72Yui3MUDEL3q44WPXw3M+fo1pZuQs4tOIBVVTaoiXEI/MxfhGDPsxsNZfoE1q66ro5aJim3XdoLFw72H+n23BaIXzbcOnz5mfPoTvYVz7KzUl5+FRxEuqkp9G/Ajia219thzg25abkRE/BpDc3pqvphHvRFys2weqvp+krbWKIX7nhDbzLOItiM8358pTwdirqpPFnMF2xLc1WvLyOwTAibpbmvHHcvttU57y5+XqNZrLe3lE/Pq8eUj2fXKfOe3pfOjzhJYtB/yll5SDFcSDiH+hRkH25+L+sdxKEAMZahrlSX8ukqMOWy/jXW2m6M9LDBc31B9LFuv6gVKg/0Szi3KAr1kGq1GMjU/aLbnq6/lRxc4XfJ98hTargX++DbMJBSiYMIe9Ck1YAxFkKEAG3xbYaKmDDgYyFK0UGYpfoWYXG+fAPPI6tJnNwb7ClP7IyF+D+bjOtCpkhz6CFrIa/I6sFtNl8auFXGMTP34sNwI/JhkgEtmDz14ySfaRcTIBInmKPE32kxyyE2Tv+thKbEVePDfW/byMM1Kmm0XdObS7oGD/MypMXFPXrCwOtoYjyyn7BV29/MZfsVzpLDdRtuIZnbpXzvlf+ev8MvYr/Gqk4H/kV/G3csdazLuyTMPsbFhzd1UabQbjFvDRmcWJxR3zcfHkVw9GfpbJmeev9F08WW8uDkaslwX6avlWGU6NRKz0g/SHtCy9J30o/ca9zX3Kfc19zn3BXQKRO8ud477hLnAfc1/G9mrzGlrfexZ5GLdn6ZZrrEohI2wVHhZywjbhUWEy8icMCGNCUdiBlq3r+xafL549HQ5jH+an+1y+LlYBifuxAvRN/lVVVOlwlCkdVm9NOL5BE4wkQ2SMlDZU97hX86EilU/lUmkQUztTE6mx1EEPh7OmdqBtAvv8HdWpbrJS6tJj3n0CWdM6busNzRV3S9KTYhqvNiqWmuroiKgYhshMjmhTh9ptWhsF7970j/SbMrsPE1suR5z7DMC+P/Hs+y7ijrQAlhyAgccjbhjPygfeBTjzhNqy28EdkUh8C+DU9+z2v/oyeH791OncxHOs5y2AtTc7nb/f73TWPkD/qwBnjX8BoJ98VQNcC+8AAAA4ZVhJZk1NACoAAAAIAAGHaQAEAAAAAQAAABoAAAAAAAKgAgAEAAAAAQAAA+igAwAEAAAAAQAAAfQAAAAA0zt2FQAAQABJREFUeAHs3QuUV1W9OPDvAPIQxRASBVNQfIWKYImmgSAkgpCaqKBe03x0uwtv5VLTNJS8lo/rI1zmI71/o64pKoqPUFylXGOlkqApgU8sIFBDAZGH4Pxnn3tnmgcwvxnmwMz8Pnut8XfOPvvsc/bnHIf5/s4+e5eUlqWQCBAgQIAAAQIECBAgQIAAga0q0GKrHt3BCRAgQIAAAQIECBAgQIAAgUxAgO5GIECAAAECBAgQIECAAAECjUBAgN4ILoJTIECAAAECBAgQIECAAAECAnT3AAECBAgQIECAAAECBAgQaAQCAvRGcBGcAgECBAgQIECAAAECBAgQEKC7BwgQIECAAAECBAgQIECAQCMQEKA3govgFAgQIECAAAECBAgQIECAgADdPUCAAAECBAgQIECAAAECBBqBgAC9EVwEp0CAAAECBAgQIECAAAECBATo7gECBAgQIECAAAECBAgQINAIBATojeAiOAUCBAgQIECAAAECBAgQICBAdw8QIECAAAECBAgQIECAAIFGICBAbwQXwSkQIECAAAECBAgQIECAAAEBunuAAAECBAgQIECAAAECBAg0AgEBeiO4CE6BAAECBAgQIECAAAECBAgI0N0DBAgQIECAAAECBAgQIECgEQgI0BvBRXAKBAgQIECAAAECBAgQIEBAgO4eIECAAAECBAgQIECAAAECjUBAgN4ILoJTIECAAAECBAgQIECAAAECAnT3AAECBAgQIECAAAECBAgQaAQCAvRGcBGcAgECBAgQIECAAAECBAgQEKC7BwgQIECAAAECBAgQIECAQCMQEKA3govgFAgQIECAAAECBAgQIECAgADdPUCAAAECBAgQIECAAAECBBqBgAC9EVwEp0CAAAECBAgQIECAAAECBATo7gECBAgQIECAAAECBAgQINAIBATom3kR1q1bt5k12J0AAQIECBAgQIAAAQIECEQI0DfzLhg5cmSsWrVqM2uxOwECBAgQIECAAAECBAgUu4AAvdjvAO0nQIAAAQIECBAgQIAAgUYhIEBvFJfBSRAgQIAAAQIECBAgQIBAsQsI0Iv9DtB+AgQIECBAgAABAgQIEGgUAgL0RnEZnAQBAgQIECBAgAABAgQIFLuAAL3Y7wDtJ0CAAAECBAgQIECAAIFGISBAbxSXwUkQIECAAAECBAgQIECAQLELCNCL/Q7QfgIECBAgQIAAAQIECBBoFAIC9EZxGZwEAQIECBAgQIAAAQIECBS7gAC92O8A7SdAgAABAgQIECBAgACBRiEgQG8Ul8FJECBAgAABAgQIECBAgECxC7RqTgALFiyIWbNmRfv27aNfv37ZZ13a99e//jVKS0s3uEu3bt2iVatmxbXBdsokQIAAAQIECBAgQIAAga0j0GwiznHjxsXVV18d69atyyRbtmyZrV900UUFyb733nux++67b7TsvHnzYu+9997odhsIECBAgAABAgQIECBAgMDmCDSLAH3atGkxfvz4OP744+Pyyy+PTz/9NH70ox/FxRdfHO3atYuxY8fWajR79uyszODBg2P//fevUb5jx4418mQQIECAAAECBAgQIECAAIGGEigp69K94T7dDXWEnOv55JNPolevXllQ/u6770Z6cp7S2rVrY5999smeqM+fP78if2Onc80118QPfvCDeOaZZ2LAgAEbK1Yjf9iwYfHggw9mXwTU2CiDAAECBAgQIECAAAECBAgUKNDkB4l79tlnIwXgp512WpUgvHXr1jFmzJhI76VPnTq1Vo70BL2kpCT69u1ba1kFCBAgQIAAAQIECBAgQIBAQws0+QD9hRdeyEwOOeSQGjbleTNnzqyxrXpGCtDTO+bpyfu9994bN954Yzz55JOxatWq6kWtEyBAgAABAgQIECBAgACBBhdo8u+gL1myJEPp1KlTDZwdd9wxy1u4cGGNbZUzUjf5119/PT7/+c9Hjx49YsWKFRWb99prr/jVr34V5cF+xQYLBAgQIECAAAECBAgQIECgAQWa/BP05cuXZxydO3euwVIeoK9cubLGtsoZr7zySnz22Wfx4YcfZoPMzZkzJ1577bW45JJL4u23344RI0bE0qVLK3ZZtmxZHHroodlPGv1dIkCAAAECBAgQIECAAAECmyvQ5J+gt23bNjNIAXb1tH79+iyrfOC46tvL1/fYY4+sW/sXvvCFOPzww8uzs2naUh3XXntt3HDDDXHVVVdl27bffvusfFo577zzKspbIECAAAECBAgQIECAAAEC9RVo8k/Qu3btmrW98hPucozyvB122KE8a4OfO+20U5xyyilVgvPygv/yL/+SLc6aNas8K1q0aJF1hU/d4Vu1avLfcVS0ywIBAgQIECBAgAABAgQIbD2BogjQu3XrVm/h9F56SuVd6etdkR0JECBAgAABAgQIECBAgMAmBJp8gL7ffvtlzUvTrVVP5Xm1DfCWRmxPc6an0durp7lz52ZZabtEgAABAgQIECBAgAABAgTyEmjyAfqAAQPigAMOiPvuu6/KU+40kFvKO+igg6J///6b9Nttt92yUdz/4z/+I0pLSyvKpuWrr746Wy/v6l6x0QIBAgQIECBAgAABAgQIEGhAgSYfoCeLNNr64sWLY+DAgfHAAw/EpEmTsuUPPvgg7rrrrirviacR20tKSqJ3794VjMcdd1xWPo3cPmjQoJg4cWJMnjw5jj766Gwu9LPPPrvWIL+iMgsECBAgQIAAAQIECBAgQKAeAs1ihLPRo0dn06SNHTs2Ro0alTF07Ngxbr/99ujbt2+tLGmU9wcffDAuvfTSuPPOO+OZZ57J9klzq6cR3C+88MJa61CAAAECBAgQIECAAAECBAhsjkBJWTfuf/bp3pyaGsG+qSlvvfVWrFmzJnr27Blt2rSp81mtXr063njjjUhTqXXv3r3W/YcNG5YF9+3atau1rAIECBAgQIAAAQIECBAgQGBjAs3iCXp541LX9RSYb05K86qnd9olAgQIECBAgAABAgQIECCwJQWaxTvoWxLMsQgQIECAAAECBAgQIECAQB4CAvQ8VNVJgAABAgQIECBAgAABAgTqKCBAryOY4gQIECBAgAABAgQIECBAIA8BAXoequokQIAAAQIECBAgQIAAAQJ1FBCg1xFMcQIECBAgQIAAAQIECBAgkIeAAD0PVXUSIECAAAECBAgQIECAAIE6CgjQ6wimOAECBAgQIECAAAECBAgQyENAgJ6HqjoJECBAgAABAgQIECBAgEAdBQTodQRTnAABAgQIECBAgAABAgQI5CEgQM9DVZ0ECBAgQIAAAQIECBAgQKCOAgL0OoIpToAAAQIECBAgQIAAAQIE8hAQoOehqk4CBAgQIECAAAECBAgQIFBHAQF6HcEUJ0CAAAECBAgQIECAAAECeQgI0PNQVScBAgQIECBAgAABAgQIEKijgAC9jmCKEyBAgAABAgQIECBAgACBPAQE6HmoqpMAAQIECBAgQIAAAQIECNRRoFUdyytOgAABAs1EYOjQoc2kJfVvxtSpU+u/sz0JECBAgAABAg0s4Al6A4OqjgABAgQIECBAgAABAgQI1EdAgF4fNfsQIECAAAECBAgQIECAAIEGFhCgNzCo6ggQIECAAAECBAgQIECAQH0EBOj1UbMPAQIECBAgQIAAAQIECBBoYAEBegODqo4AAQIECBAgQIAAAQIECNRHQIBeHzX7ECBAgAABAgQIECBAgACBBhYQoDcwqOoIECBAgAABAgQIECBAgEB9BATo9VGzDwECBAgQIECAAAECBAgQaGABAXoDg6qOAAECBAgQIECAAAECBAjUR0CAXh81+xAgQIAAAQIECBAgQIAAgQYWEKA3MKjqCBAgQIAAAQIECBAgQIBAfQQE6PVRsw8BAgQIECBAgAABAgQIEGhgAQF6A4OqjgABAgQIECBAgAABAgQI1EdAgF4fNfsQIECAAAECBAgQIECAAIEGFhCgNzCo6ggQIECAAAECBAgQIECAQH0EBOj1UbMPAQIECBAgQIAAAQIECBBoYAEBegODqo4AAQIECBAgQIAAAQIECNRHQIBeHzX7ECBAgAABAgQIECBAgACBBhYQoDcwqOoIECBAgAABAgQIECBAgEB9BATo9VGzDwECBAgQIECAAAECBAgQaGABAXoDg6qOAAECBAgQIECAAAECBAjUR0CAXh81+xAgQIAAAQIECBAgQIAAgQYWEKA3MKjqCBAgQIAAAQIECBAgQIBAfQQE6PVRsw8BAgQIECBAgAABAgQIEGhgAQF6A4OqjgABAgQIECBAgAABAgQI1EdAgF4fNfsQIECAAAECBAgQIECAAIEGFhCgNzCo6ggQIECAAAECBAgQIECAQH0EBOj1UbMPAQIECBAgQIAAAQIECBBoYAEBegODqo4AAQIECBAgQIAAAQIECNRHQIBeHzX7ECBAgAABAgQIECBAgACBBhYQoDcwqOoIECBAgAABAgQIECBAgEB9BATo9VGzDwECBAgQIECAAAECBAgQaGABAXoDg6qOAAECBAgQIECAAAECBAjUR0CAXh81+xAgQIAAAQIECBAgQIAAgQYWEKA3MKjqCBAgQIAAAQIECBAgQIBAfQQE6PVRsw8BAgQIECBAgAABAgQIEGhgAQF6A4OqjgABAgQIECBAgAABAgQI1EdAgF4fNfsQIECAAAECBAgQIECAAIEGFhCgNzCo6ggQIECAAAECBAgQIECAQH0EWtVnp8a6z4IFC2LWrFnRvn376NevX/a5Oef6u9/9Ltq2bRtf+cpXNqca+xIgQIAAAQIECBAgQIAAgVoFms0T9HHjxkWPHj1i5MiRcdRRR8UOO+wQ1157ba0AGyvwxBNPZPVceeWVGysinwABAgQIECBAgAABAgQINJhAswjQp02bFuPHj48RI0bESy+9FM8//3wMHjw4Lr744pgwYUKdsd5///0466yz6ryfHQgQIECAAAECBAgQIECAQH0FmnyA/sknn8S5554b3bp1i0mTJkWfPn3ikEMOiSlTpkT37t2zp+jr16+vk8/ZZ58dn332WZ32UZgAAQIECBAgQIAAAQIECGyOQJMP0J999tmYP39+nHbaadGyZcsKi9atW8eYMWMivZc+derUivzaFu64444suE+fKZWUlNS2i+0ECBAgQIAAAQIECBAgQGCzBZp8gP7CCy9kCOmpefVUnjdz5szqmza4/sYbb8T3v//9+Ld/+7cYOnToBsvIJECAAAECBAgQIECAAAECeQg0+QB9yZIlmUunTp1q+Oy4445Z3sKFC2tsq56xbt26OPXUU2PXXXfdrMHlqtdrnQABAgQIECBAgAABAgQIFCLQ5KdZW758edbOzp0712hveYC+cuXKGtuqZ6TR2tMUbTNmzIhtt902Vq9eXb1IxfqKFSvi5JNPztY/+OCDinwLBAgQIECAAAECBAgQIECgvgJNPkBP85SntKFB3coHh6v8bvqGoFJQ/pOf/CQuv/zy+PKXv7yhIlXy2rVrFz/+8Y+zvDRSvESAAAECBAgQIECAAAECBDZXoMkH6F27ds0Mli5dWsOiPC/Nib6xlJ6GpwHmDjzwwPje974XaVT4lMqfoKcgP+W1atUq0sBzKaXlgw8+OFsuz8tW/IcAAQIECBAgQIAAAQIECNRToCgC9DQF28ZS6tb+zjvvZJs3FMg//fTT0b59+zjllFPi3nvv3Vg18gkQIECAAAECBAgQIECAwGYJNPkAfb/99ssA0nRrxx9/fBWMlJdS+WjuVTb+30p6Aj927Ngam9KgcT//+c9jt912i69//evRt2/fGmVkECBAgAABAgQIECBAgACBhhIoKS1LDVXZ1qondU9///33Y968edGhQ4fsNJYtWxb77LNP7LLLLvHiiy9m3dLrcn6pi3t61/zoo4/e5Dzqw4YNiwcffDArW5f6lSVAgMDWFjCdZGzy9/vWvj6OT4AAAQIECBSfQJOfZi1dsksuuSQWL14cAwcOjAceeCAmTZqULacR1u+6664qwfkrr7wSJSUl0bt37+K72lpMgAABAgQIECBAgAABAo1WoMl3cU+yo0ePzkZxT13VR40alWF37Ngxbr/9dl3TG+2t58QIECBAgAABAgQIECBAoLJAswjQU4NOPfXUGDNmTLz11luxZs2a6NmzZ7Rp06ZyW7Pl1B2+kF79afq2QsrVOIAMAgQIECBAgAABAgQIECBQD4FmE6Cntqeu6ykwlwgQIECAAAECBAgQIECAQFMTaBbvoDc1dOdLgAABAgQIECBAgAABAgSqCwjQq4tYJ0CAAAECBAgQIECAAAECW0FAgL4V0B2SAAECBAgQIECAAAECBAhUFxCgVxexToAAAQIECBAgQIAAAQIEtoKAAH0roDskAQIECBAgQIAAAQIECBCoLiBAry5inQABAgQIECBAgAABAgQIbAUBAfpWQHdIAgQIECBAgAABAgQIECBQXUCAXl3EOgECBAgQIECAAAECBAgQ2AoCAvStgO6QBAgQIECAAAECBAgQIECguoAAvbqIdQIECBAgQIAAAQIECBAgsBUEBOhbAd0hCRAgQIAAAQIECBAgQIBAdQEBenUR6wQIECBAgAABAgQIECBAYCsICNC3ArpDEiBAgAABAgQIECBAgACB6gIC9Ooi1gkQIECAAAECBAgQIECAwFYQEKBvBXSHJECAAAECBAgQIECAAAEC1QUE6NVFrBMgQIAAAQIECBAgQIAAga0gIEDfCugOSYAAAQIECBAgQIAAAQIEqgsI0KuLWCdAgAABAgQIECBAgAABAltBQIC+FdAdkgABAgQIECBAgAABAgQIVBcQoFcXsU6AAAECBAgQIECAAAECBLaCwFYN0BctWhRLly7dCs12SAIECBAgQIAAAQIECBAg0LgEtkiAvnDhwvjBD34Qv/3tb7PWr127No4++ujo1q1b7LTTTvGNb3wjPv3008Yl42wIECBAgAABAgQIECBAgMAWFMg9QF+3bl0cd9xxcc0118Ts2bOzpl199dXx1FNPZQF6v3794qGHHorvfe97W7DZDkWAAAECBAgQIECAAAECBBqXQO4B+uOPPx4zZ86MSy+9NL7//e9nrf/lL38Zbdu2jZdffjn+8Ic/xOmnnx4pr7S0tHHpOBsCBAgQIECAAAECBAgQILCFBHIP0F977bVo3bp1XHLJJdGmTZuYO3duvPPOO3HkkUdGp06dsmYee+yxsWLFinjzzTe3ULMdhgABAgQIECBAgAABAgQINC6B3AP0999/P7bffvvYbrvtspaXv4c+dOjQColVq1Zly+nddIkAAQIECBAgQIAAAQIECBSjQO4B+h577BH/+Mc/Ys6cOVkX9vvvvz9zPuaYYyq8J0+eHC1atIjdd9+9Is8CAQIECBAgQIAAAQIECBAoJoHcA/QTTzwx69o+YMCAGDhwYPzxj3/MPvfee++YN29efPnLX45HHnkkTjrppIqn7MV0AbSVAAECBAgQIECAAAECBAgkgdwD9F122SUee+yxrJv7//zP/8RXv/rV+M1vfpPpp+7vaQC5FLjfcsstrggBAgQIECBAgAABAgQIEChagVZbouWDBw+Ot99+O1avXp2N3l5+zD59+mRTr/Xu3bs8yycBAgQIECBAgAABAgQIEChKgdyfoFdWTVOrLVy4MJ599tnsvfRtttkmvvjFL1YuYpkAAQIECBAgQIAAAQIECBSlwBYJ0NP85vfcc0907do1dt1112yKtdmzZ8crr7wSX/rSl+LFF18sSnyNJkCAAAECBAgQIECAAAEC5QJbJEC/7LLL4pvf/GZ89NFH0b9///Jjx/r167N50Q877LAsgK/YYIEAAQIECBAgQIAAAQIECBSZQO4BenpK/tOf/jTGjBkT7733XkyYMKGCuF+/fvHSSy9F9+7d44c//GEWsFdstECAAAECBAgQIECAAAECBIpIIPcA/ZlnnomSkpK47bbbNjiNWq9eveLcc8/N3k1fsGBBEdFrKgECBAgQIECAAAECBAgQ+KdA7gH6/PnzI021tv322//zqNWW0mjuKX344YfVtlglQIAAAQIECBAgQIAAAQLFIZB7gL7PPvtkT8cXLVq0UdEZM2ZEixYtYu+9995oGRsIECBAgAABAgQIECBAgEBzFsg9QB80aFCk6dRGjx4df/rTn2pYPvHEE9k76ocffnhsu+22NbbLIECAAAECBAgQIECAAAECxSCQe4C+1157xdVXXx3Tp0/PplQ78cQTM9frrrsuWx8+fHj2jvodd9xRDN7aSIAAAQIECBAgQIAAAQIENiiQe4CejnrBBRfElClTInV3f+ONN7ITefLJJ7MR3I8++uiYNWtW7Lvvvhs8QZkECBAgQIAAAQIECBAgQKAYBFptqUaOGDEi0k+aCz0F6W3bto2ePXtGu3btttQpOA4BAgQIECBAgAABAgQIEGi0AlvkCXrl1q9cuTI++eST6Nq1a7Rs2TI+/fTTypstEyBAgAABAgQIECBAgACBohTYIgF6aWlp3HPPPVlQvuuuu8aRRx4Zs2fPjldeeSV7D/3FF18sSnyNJkCAAAECBAgQIECAAAEC5QJbJEC/7LLL4pvf/GbWvb1///7lx47169fH3Llz47DDDssC+IoNFggQIECAAAECBAgQIECAQJEJ5B6gp6fkP/3pT2PMmDHx3nvvxYQJEyqI+/Xrlw0U17179/jhD3+YBewVGy0QIECAAAECBAgQIECAAIEiEsg9QH/mmWeyadRuu+222G677WrQ9urVK84999xYuHBhLFiwoMZ2GQQIECBAgAABAgQIECBAoBgEcg/Q58+fH7vssktsv/32G/Xs06dPtu3DDz/caBkbCBAgQIAAAQIECBAgQIBAcxbIPUBPc5+np+OLFi3aqOOMGTOiRYsWsffee2+0jA0ECBAgQIAAAQIECBAgQKA5C+QeoA8aNCi22WabGD16dPzpT3+qYfnEE09k76gffvjhse2229bYLoMAAQIECBAgQIAAAQIECBSDQO4B+l577RVXX311TJ8+PZtS7cQTT8xcr7vuumx9+PDh2Tvqd9xxRzF4ayMBAgQIECBAgAABAgQIENigQO4BejrqBRdcEFOmTInU3f2NN97ITuTJJ5/MRnA/+uijY9asWbHvvvtu8ARlEiBAgAABAgQIECBAgACBYhBotaUaOWLEiEg/H330URakt23bNnr27Bnt2rXbUqfgOAQIECBAgAABAgQIECBAoNEK5P4E/dVXX40hQ4ZkXdyTwuc+97n48pe/HAcccIDgvNHeFk6MAAECBAgQIECAAAECBLa0QO4B+iOPPBJPP/30Jkdx39KNdjwCBAgQIECAAAECBAgQINDYBHIP0Dt16pS1+eOPP25sbXc+BAgQIECAAAECBAgQIECg0Qjk/g76WWedFamb+4UXXhirV6+Ofv36RY8ePaJ169Y1EDp06FAjTwYBAgQIECBAgAABAgQIECgGgdwD9EcffTTSz4oVK2Ls2LGbNC0tLd3k9to2LliwIBsRvn379tkXAemzrmnevHkxZ86c6Nq1a/Tt2zebw72udShPgAABAgQIECBAgAABAgTqKpB7gN6xY8c48MADs5+6nlxdyo8bNy6bb33dunXZbi1btszWL7roooKqWbp0aXzzm9/Mvkwo3yGNMH/TTTfFueeeW57lkwABAgQIECBAgAABAgQI5CKQe4D+pS99Ka6//vpsDvSNteBvf/tbPP/88xvbXGv+tGnTYvz48XH88cfH5ZdfHp9++mn86Ec/iosvvjgbKb62J/fpAKecckqkes4555w4++yz4+9//3t23uedd16kYP9b3/pWreehAAECBAgQIECAAAECBAgQqK9ASVm38s3rV17LkX/5y1/GGWecEZs6TAqkr7322liyZEnstNNOtdRYdfMnn3wSvXr1yoLyd999NwumU4m1a9dmXwqkJ+rz58+vyK+69/+uzZw5M5v6LX2Z8OKLL1YUeeedd2LPPfeMww47LP7whz9U5FdeGDZsWDz44IOmjKuMYpkAgSYhMHTo0CZxnnme5NSpU/OsXt0ECBAgQIAAgToJ5PIE/b777ovyUdtnzJiRndBdd921wRNbuXJlTJ48OVq1ahXbb7/9BstsKvPZZ5/NAvAU5Kcn3eUpDUI3ZsyYrJt7+gNs+PDh5ZtqfG677bbZk/f+/ftX2ZYGs0s/f/nLX6rkWyFAgAABAgQIECBAgAABAg0tkEuA/te//jWqv/uduo1vKp188sn1egr9wgsvZNUecsghNaovz0tPyDcVoH/xi1/MushXr2DWrFlZ8H/CCSdU32SdAAECBAgQIECAAAECBAg0qEAuAfp3v/vdbBT01K39ueeei9tvvz0mTpxY48RLSkqyUdJTt/YjjjiixvZCMlK3+JTK51uvvM+OO+6YrS5cuLBy9iaX0znfc8898eSTT8bjjz+edZ+/7rrrNrmPjQQIECBAgAABAgQIECBAYHMFcgnQt9lmmzj11FOzc9t9992zKdZOO+20zT3XDe6/fPnyLL9z5841tpcH6KkbfaEpDQ535plnVhQfOXJkdOvWrWI9LaTu+9/+9rezvDT6u0SAAAECBAgQIECAAAECBDZXoMXmVlDb/l/96lfj17/+dW3F6r29bdu22b6fffZZjTrWr1+f5VV+N71GoWoZaVq41EU/DRaXRnC/5ppr4qCDDqp4pz4Vb9OmTTbwXRr8Lr2/LhEgQIAAAQIECBAgQIAAgc0VyD1A39wTrG3/rl27ZkU29CS7PG+HHXaorZqK7Wnu8y984QuRRnS/7bbb4rjjjos5c+ZkXd7LC6UeAkOGDMl+yr8gKN/mkwABAgQIECBAgAABAgQI1EegKAL06l3U6wJVPv95eh9dIkCAAAECBAgQIECAAAECeQk0+QB9v/32y2zSdGvVU3le+Wju1beXr6dB4FLX9t/97nflWRWfLVr8L9F2221XkWeBAAECBAgQIECAAAECBAg0tECDB+jvvfdepIHWtlQaMGBAHHDAAZHmXi8fMC4de9myZVleen+8+vzm1c9t3333jY8++igmTJhQfVPcfPPNWd5RRx1VY5sMAgQIECBAgAABAgQIECDQUAINHqD/13/9VzbF2qJFi7JznD17dowbN66hzneD9VxyySWxePHiGDhwYDzwwAMxadKkbPmDDz6Iu+66K1q1+udg9a+88kqk6d169+5dUdexxx4bxxxzTDz88MPxta99Lf77v/87Wx46dGg88cQTMWrUqPj6179eUd4CAQIECBAgQIAAAQIECBBoaIF/Rq4NVHP5lGbvv/9+FqingHj8+PFx5ZVXNtARalYzevToSKO4jx07NgumU4nUZT3Nv963b9+aO1TLSQH7b37zm7jsssvi1ltvjWnTpmUl0gjtP/7xj+Piiy+utodVAgQIECBAgAABAgQIECDQsAIlpWWpIat87LHHYsSIEZHeDR8+fHi8/fbb8dBDD2WBbm3HSQHy5qTUlLfeeivWrFkTPXv2zKZDq2t9q1atinnz5mXTp+25555R2xRtw4YNiwcffDDS6O8SAQIEmpJA6iVU7Gnq1KnFTqD9BAgQIECAQCMSaPAAPQW46Yn2I488UudmNvB3BXU+fn12EKDXR80+BAg0BgEBeoQAvTHcic6BAAECBAgQKBdo8C7u6Uny5MmTsyfnCxcujDQ92bXXXhtPP/10+TF9EiBAgAABAgQIECBAgAABAtUEGjxAT/Wnd7pT9/D0k0ZTf+qpp8Io6NXkrRIgQIAAAQIECBAgQIAAgUoCuQTolerP3kdP76SnlAaQmzt3brz55pvRtm3b2H333bPR1FNALxEgQIAAAQIECBAgQIAAgWIWyD1AT7jr1q3L5hO/4oor4uOPP67i3aNHj6xLfOVpz6oUsEKAAAECBAgQIECAAAECBIpAYIsE6BdddFHceOONsfPOO8fpp5+ePTlfsWJFzJkzJ6ZMmRKDBg3KpjYrZEq0IrgmmkiAAAECBAgQIECAAAECRSiQe4A+e/bsuOmmm7Ku7vfee2+0b9++CvOrr74aAwYMiBTEG0iuCo0VAgQIECBAgAABAgQIECgigRZ5t3X69OnZIe6+++4awXnasP/++8e4ceNixowZsXbt2rxPR/0ECBAgQIAAAQIECBAgQKBRCuQeoL/zzjtZ1/bOnTtvFKBXr16R5k9PA8hJBAgQIECAAAECBAgQIECgGAVyD9DTSO2LFy+OJUuWbNQ3dYNPqXv37tmn/xAgQIAAAQIECBAgQIAAgWITyD1AHzJkSLRo0SLOOOOMbE706sAzZ86M8ePHx8EHHxwdOnSovtk6AQIECBAgQIAAAQIECBAoCoHcB4lL3dfPP//8bBT3PfbYI4YPH56N4p7mRE+juD/11FPRqlWruPPOO4sCXCMJECBAgAABAgQIECBAgMCGBHIP0NNBb7jhhmwwuAsuuCAmTpxY5TwOPfTQuOWWW6JPnz5V8q0QIECAAAECBAgQIECAAIFiEtgiAXoCPeuss+LMM8+Md999N+bNm5eN6L7XXntFly5dislbWwkQIECAAAECBAgQIECAwAYFtliAno5eUlKSDQRnMLgNXguZBAgQIECAAAECBAgQIFDEArkPElfEtppOgAABAgQIECBAgAABAgQKFhCgF0ylIAECBAgQIECAAAECBAgQyE9AgJ6frZoJECBAgAABAgQIECBAgEDBAgL0gqkUJECAAAECBAgQIECAAAEC+QkI0POzVTMBAgQIECBAgAABAgQIEChYQIBeMJWCBAgQIECAAAECBAgQIEAgP4EtMs3a0qVL44477ogXXnghli9fHuvXr99gi37/+99vMF8mAQIECBAgQIAAAQIECBBo7gK5B+hr1qyJwYMHx6xZs5q7pfYRIECAAAECBAgQIECAAIF6C+Texf2hhx7KgvPhw4fHa6+9FitXrox169Zt8KferbAjAQIECBAgQIAAAQIECBBo4gK5P0GfPXt2RvSLX/widt555ybO5fQJECBAgAABAgQIECBAgEA+Ark/Qe/SpUu0atUq2rRpk08L1EqAAAECBAgQIECAAAECBJqBQO4B+pFHHpkNCnffffc1Ay5NIECAAAECBAgQIECAAAEC+Qjk3sW9b9++cf3118ePfvSjKCkpia985StZV/cWLWp+N9CpU6d8WqlWAgQIECBAgAABAgQIECDQyAVyD9AnTZqUBegffPBBfPvb394kR2lp6Sa320iAAAECBAgQIECAAAECBJqrQO4B+uc///k49NBDm6ufdhEgQIBAMxYYOnRoM25d4U2bOnVq4YWVJECAAAECBOotkHuAnt5BTz8SAQIECBAgQIAAAQIECBAgsHGB3AP0yodOc6DPnTs33nzzzWjbtm3svvvu0bt37+zd9MrlLBMgQIAAAQIECBAgQIAAgWIT2CIB+rp16+Lmm2+OK664Ij7++OMqxj169IjJkydngXqVDVYIECBAgAABAgQIECBAgEARCWyRAP2iiy6KG2+8MRu9/fTTT8+enK9YsSLmzJkTU6ZMiUGDBsW0adMijfguESBAgAABAgQIECBAgACBYhTIPUCfPXt23HTTTTFixIi49957o3379lWcX3311RgwYECkIP7pp5+uss0KAQIECBAgQIAAAQIECBAoFoGak5E3cMunT5+e1Xj33XfXCM7Thv333z/GjRsXM2bMiLVr1zbw0VVHgAABAgQIECBAgAABAgSahkDuAfo777yTdW3v3LnzRkV69eoVq1atygaQ22ghGwgQIECAAAECBAgQIECAQDMWyD1ATyO1L168OJYsWbJRxtQNPqXu3btnn/5DgAABAgQIECBAgAABAgSKTSD3AH3IkCHRokWLOOOMM2LZsmU1fGfOnBnjx4+Pgw8+ODp06FBjuwwCBAgQIECAAAECBAgQIFAMArkPEpe6r59//vnZKO577LFHDB8+PBvFPc2JnkZxf+qpp6JVq1Zx5513FoO3NhIgQIAAAQIECBAgQIAAgQ0K5B6gp6PecMMN2WBwF1xwQUycOLHKiRx66KFxyy23RJ8+farkWyFAgAABAgQIECBAgAABAsUksEUC9AR61llnxZlnnhnvvvtuzJs3LxvRfa+99oouXboUk7e2EiBAgAABAgQIECBAgACBDQpssQA9Hb2kpCQbCM5gcBu8FjIJECBAgAABAgQIECBAoIgFGjxAT3OZf/TRR9GmTZvYYYcdYvXq1bF8+fKCiHfaaaeCyilEgAABAgQIECBAgAABAgSam0CDB+hPPPFEHH/88dlgcI899ljcf//92QjuhcCVlpYWUkwZAgQIECBAgAABAgQIECDQ7AQaPEDv2rVrnHjiidG3b98MK82DntYlAgQIECBAgAABAgQIECBAYOMCDR6gH3LIITFp0qSKI6bR2a+66qrYZ599KvKqL/ztb3+L559/vnq2dQIECBAgQIAAAQIECBAgUDQCLfJu6cMPPxz77rvvJg+TplkbNWpUvPfee5ssZyMBAgQIECBAgAABAgQIEGiuAg3+BD1B3XffffHxxx9nZjNmzMg+77rrruyz+n9WrlwZkydPjlatWsX2229ffbN1AgQIECBAgAABAgQIECBQFAK5BOh//etf46KLLqoCePbZZ1dZr75y8sknR7t27apnWydAgAABAgQIECBAgAABAkUhkEuA/t3vfjfSYHFpVPbnnnsubr/99pg4cWIN0DQv+jbbbBNperUjjjiixnYZBAgQIECAAAECBAgQIECgWARyCdBT0H3qqadmhmkU9xUrVsRpp51WLKbaSYAAAQIECBAgQIAAAQIE6iyQ+yBxHTt2zAZ/mz59ep1Pzg4ECBAgQIAAAQIECBAgQKBYBHIP0B955JF4+umnY9GiRcViqp0ECBAgQIAAAQIECBAgQKDOArkH6J06dcpOqnxU9zqfoR0IECBAgAABAgQIECBAgEARCOTyDnplt7POOiteffXVuPDCC2P16tXRr1+/6NGjR7Ru3bpysWy5Q4cONfJkECBAgAABAgQIECBAgACBYhDIPUB/9NFHI/2kgeLGjh27SdM06rtEgAABAgQIECBAgAABAgSKUSD3AD0NEnfggQdmP8UIrM0ECBAgQIAAAQIECBAgQKAQgdwD9EGDBkX62RJpwYIFMWvWrGjfvn3WlT591jW9/fbbMXfu3Pj0009j3333jX322aeuVShPgAABAgQIECBAgAABAgTqLJB7gF75jFauXJkFv2+++Wa0bds20hzpvXv3jpKSksrF6rU8bty4uPrqq2PdunXZ/i1btszWL7roooLqW7x4cXz729+ONOp85TRw4MD4xS9+EXvssUflbMsECBAgQIAAAQIECBAgQKBBBXIfxT2dbQqa//M//zN23nnn+NKXvhSnnHJKHHfccdGnT5/Yc8894+WXX96sRk2bNi3Gjx8fI0aMiJdeeimef/75GDx4cFx88cUxYcKEWuv+7LPPsnNKwflJJ50UTzzxRDzzzDORBrhLnyNHjswGuKu1IgUIECBAgAABAgQIECBAgEA9BUrKBmbLfWS273//+3HjjTdmAfrxxx+fPTlPg8bNmTMnpkyZEjvssEOkILtv3751bsYnn3wSvXr1yrqkv/vuu5GenKe0du3arHt6+nJg/vz5FfkbOsCzzz4bRx55ZBx22GExY8aMKkWGDx+eBez3339/jBo1qsq2tDJs2LB48MEHo127djW2ySBAgEBjFhg6dGhjPr0tcm5Tp07d5HEY/S9PbU6bRLSRAAECBAgQKFgg9yfos2fPjptuuil7up26tt96663Zk+2rrroqHnrooUjbUyq0K3r1lqXgOgXgp512WpUgPE3jNmbMmEjvpdf2h0Xav3v37tkT8+r1n3766VlW+jJBIkCAAAECBAgQIECAAAECeQnkHqBPnz49O/e77747G7ytekP233//SO+PpyfX6al3XdMLL7yQ7XLIIYfU2LU8b+bMmTW2Vc4444wz4p133omzzz67cna2nAaNSyl1xZcIECBAgAABAgQIECBAgEBeArkPEpcC3/TueefOnTfahtRFfdWqVdkAcmlKtrqkJUuWZMU7depUY7cdd9wxy1u4cGGNbYVkfPDBB1nX/A4dOmTvtJfvkwa7u+SSS7LVjz76qDzbJwECBAgQIECAAAECBAgQqLdA7gF6Gqk9jZCeAukuXbps8ETLu7mnbuZ1TcuXL8922dAXAOUBegqo65rSPscee2ykID2N4p6+ZChP22yzTRxxxBHZ6iuvvFKe7ZMAAQIECBAgQIAAAQIECNRbIPcu7kOGDIkWLVpE6ka+bNmyGieaup+nEdgPPvjgSE+q65rSdG0ppZHYq6f169dnWeUDx1XfvrH1FJSn806jwZ9//vnxrW99q0rR9H57Gu09/Wy77bZVtlkhQIAAAQIECBAgQIAAAQL1Ecj9CXrqvp6C3DSKe5pLPI2Knp6qpyfUaeC1p556Klq1ahV33nlnfc4/unbtmu23dOnSGvuX56VR4gtNb731VqRRe9OAdj/84Q8jDWYnESBAgAABAgQIECBAgACBvAVyD9BTA2644YZIg8FdcMEFMXHixCptOvTQQ+OWW27J5kSvsqHAlUIC9G7duhVU26uvvhpf+9rX4v3334877rgjzjnnnIL2U4gAAQIECBAgQIAAAQIECGyuwBYJ0NNJnnXWWXHmmWdGmqt83rx52Yjue+2110bfSy+0Yfvtt19WNE23luZYr5xSXkrlo7lX3lZ9OXW1P/roo7P51B9//PEsUK9exjoBAgQIECBAgAABAgQIEMhLIPd30CufeElJSfbOdurS3rFjx02O7F55v00tDxgwIA444IC47777onzAuFQ+ve+e8g466KDo37//pqrIRpAfNWpU1u3+ySefFJxvUstGAgQIECBAgAABAgQIEMhDYIs9QU9dxq+88spYtGhRRTvSAGtf/epX47/+679il112qciv60Ka8mzMmDExcODAbPqz0tLS+MlPfpKNwP7EE09k77iX15lGXe/du3ek6dxefvnlLDuVnT9/fvY++zXXXFNetMpnGtF9Q/OkVylkhQABAgQIECBAgAABAgQI1FNgiwTo5513XvZOd3oXPHVz33XXXbOn3XPnzs0GiUtd0FO38rrOgV7e5tGjR2ejuI8dOzbSk/CU0hP622+/Pfr27VtebKOfjz76aLYtfXnwyCOPbLBcGthOIkCAAAECBAgQIECAAAECeQnkHqC/8cYbFQOuTZgwIdq0aVOlLekp9uDBg7MB5KZNm1ZlW11WTj311OwpehqFfc2aNdGzZ88ax0r1pS8B0hP2ymnWrFmVVy0TIECAAAECBAgQIECAAIEtLpD7O+jpyfjnPve5uPXWWzcYMKfu5mke9OnTp8eKFSs2CyC9454C8zS1W/UvAjarYjsTIECAAAECBAgQIECAAIGcBXIP0D/88MPYaaedqrwHXr1Ne+65Z6xdu7bK++nVy1gnQIAAAQIECBAgQIAAAQLNWSD3AD2NoP7666/Hc889t1HH9N53ly5dYu+9995oGRsIECBAgAABAgQIECBAgEBzFsg9QD/qqKPiu9/9bowcOTJ+9rOfZSOrJ9D0HngaOT1t+/nPfx7XX399LF26NP7xj39U/DRneG0jQIAAAQIECBAgQIAAAQKVBXIfJC7NRT5x4sRIXd3//d//Pftp3759fPrpp1m39vKTOf3008sXs8/WrVtng71VybRCgAABAgQIECBAgAABAgSaqUDuAXrqup66udc1tWqV+6nV9ZSUJ0CAAAECBAgQIECAAAECuQnkHgUfeeSRkX4kAgQIECBAgAABAgQIECBAYOMCuQfolQ+9cuXKmDt3brz55pvRtm3b2H333SNNs5amR5MIECBAgAABAgQIECBAgEAxC2yRAH3dunVx8803xxVXXBEff/xxFe8ePXrE5MmTs0C9ygYrBAgQIECAAAECBAgQIECgiAS2SIB+0UUXxY033hg777xzpMHg0pPzFStWxJw5c2LKlCkxaNCgmDZtWvTt27eI6DWVAAECBAgQIECAAAECBAj8UyD3AH327Nlx0003xYgRI+Lee++NNIJ75fTqq6/GgAEDIgXxTz/9dOVNlgkQIECAAAECBAgQIECAQNEI5D4P+vTp0zPMu+++u0Zwnjbsv//+MW7cuJgxY0aVadeK5gpoKAECBAgQIECAAAECBAgQKBPIPUB/5513sq7tnTt33ih4r169YtWqVdkAchstZAMBAgQIECBAgAABAgQIEGjGArkH6Ol988WLF8eSJUs2ypi6wafUvXv37NN/CBAgQIAAAQIECBAgQIBAsQnkHqAPGTIkWrRoEWeccUYsW7ashu/MmTNj/PjxcfDBB0eHDh1qbJdBgAABAgQIECBAgAABAgSKQSD3QeJS9/Xzzz8/G8V9jz32iOHDh2ejuKc50dMo7k899VS0atUq7rzzzmLw1kYCBAgQIECAAAECBAgQILBBgdwD9HTUG264IRsM7oILLoiJEydWOZFDDz00brnllujTp0+VfCsECBAgQIAAAQIECBAgQKCYBLZIgJ5AzzrrrDjzzDPj3XffjXnz5mUjuu+1117RpUuXYvLWVgIECBAgQIAAAQIECBAgsEGB3N9BT/Ocp/fQ03RrJSUl2UBwRx99dBxxxBGC8w1eEpkECBAgQIAAAQIECBAgUIwCuQfojzzySDz99NOxaNGiYvTVZgIECBAgQIAAAQIECBAgUJBA7gF6p06dshP5+OOPCzohhQgQIECAAAECBAgQIECAQDEK5P4Oenr3PHVzv/DCC2P16tXRr1+/6NGjR7Ru3bqGt2nWapDIIECAAAECBAgQIECAAIEiEcg9QH/00Ucj/axYsSLGjh27SdbS0tJNbreRAAECBAgQIECAAAECBAg0V4HcA/SOHTvGgQcemP00V0TtIkCAAAECBAgQIECAAAECmyuQe4A+aNCgSD8SAQIECBAgQIAAAQIECBAgsHGB3AP06od+77334s9//nPsvPPOse+++0bLli2rF7FOgAABAgQIECBAgAABAgSKTiD3UdzLRe+4447o1q1bNvf54MGDY//99480KNzQoUPj73//e3kxnwQIECBAgAABAgQIECBAoCgFtsgT9PPOOy/KA/Qzzzwzdt1111i+fHnMnTs3nnrqqTjkkEPi8ccf9556Ud6CGk2AAAECBAgQIECAAAECSSD3AP2NN97IgvNzzjknJkyYEG3atKki//LLL0d6on7BBRfEtGnTqmyzQoAAAQIECBAgQIAAAQIEikUg9y7u6cn45z73ubj11ltrBOcJuXfv3jF+/PiYPn16NhVbscBrJwECBAgQIECAAAECBAgQqCyQe4D+4Ycfxk477RStWm38Yf2ee+4Za9eujUWLFlU+N8sECBAgQIAAAQIECBAgQKBoBHIP0Pv37x+vv/56PPfccxtFfeSRR7LB4/bee++NlrGBAAECBAgQIECAAAECBAg0Z4HcA/Sjjjoqvvvd78bIkSPjZz/7WXzwwQeZZ2lpacyfPz/b9vOf/zyuv/76WLp0afzjH/+o+GnO8NpGgAABAgQIECBAgAABAgQqC2y833nlUpuxfN9998XEiRMjdXX/93//9+ynffv28emnn2bd2surPv3008sXs8/WrVvHmjVrquRZIUCAAAECBAgQIECAAAECzVUg9wC9S5cukbq51zVt6p31utalPAECBAgQIECAAAECBAgQaOwCuQfoRx55ZKQfiQABAgQIECBAgAABAgQIENi4QO4B+sYPbQsBAgTyExg6dGh+lTeRmqdOndpEztRpEiBAgAABAgQIJIHcB4nDTIAAAQIECBAgQIAAAQIECNQuIECv3UgJAgQIECBAgAABAgQIECCQu4AAPXdiByBAgAABAgQIECBAgAABArULCNBrN1KCAAECBAgQIECAAAECBAjkLiBAz53YAQgQIECAAAECBAgQIECAQO0CAvTajZQgQIAAAQIECBAgQIAAAQK5CwjQcyd2AAIECBAgQIAAAQIECBAgULuAAL12IyUIECBAgAABAgQIECBAgEDuAgL03IkdgAABAgQIECBAgAABAgQI1C4gQK/dSAkCBAgQIECAAAECBAgQIJC7gAA9d2IHIECAAAECBAgQIECAAAECtQsI0Gs3UoIAAQIECBAgQIAAAQIECOQuIEDPndgBCBAgQIAAAQIECBAgQIBA7QIC9NqNlCBAgAABAgQIECBAgAABArkLCNBzJ3YAAgQIECBAgAABAgQIECBQu4AAvXYjJQgQIECAAAECBAgQIECAQO4CAvTciR2AAAECBAgQIECAAAECBAjULiBAr91ICQIECBAgQIAAAQIECBAgkLuAAD13YgcgQIAAAQIECBAgQIAAAQK1CwjQazdSggABAgQIECBAgAABAgQI5C7QrAL0BQsWxKOPPhq/+93vYuXKlZuFN3/+/Pj1r3+9WXXYmQABAgQIECBAgAABAgQIFCrQbAL0cePGRY8ePWLkyJFx1FFHxQ477BDXXnttoQ5Vyi1fvjyGDRsW3/72t6vkWyFAgAABAgQIECBAgAABAnkJNIsAfdq0aTF+/PgYMWJEvPTSS/H888/H4MGD4+KLL44JEybUye7DDz+Mk046Kf7yl7/UaT+FCRAgQIAAAQIECBAgQIDA5gg0+QD9k08+iXPPPTe6desWkyZNij59+sQhhxwSU6ZMie7du2dP0devX1+Q0eTJk6NXr17x5JNPRuvWrQvaRyECBAgQIECAAAECBAgQINAQAk0+QH/22WcjvS9+2mmnRcuWLStMUoA9ZsyYSO+lT506tSJ/Ywu//e1v44QTToi1a9fGI488kgXqGysrnwABAgQIECBAgAABAgQINLRAkw/QX3jhhcwkPTWvnsrzZs6cWX1TjfVWrVrFZZddFq+//nr2HnuNAjIIECBAgAABAgQIECBAgECOAq1yrHuLVL1kyZLsOJ06dapxvB133DHLW7hwYY1t1TOGDBkS6aeQlLrVlw9AlwaUkwgQIECAAAECBAgQIECAwOYKNPkAvTxA7ty5cw2L8gB9c6dcq15xixYtYuedd86yK3err17OOgECBAgQIECAAAECBAgQKFSgyQfobdu2zdr62Wef1Whz+eBwDR1Ep2OWT8GWBqOTCBAgQIAAAQIECBAgQIDA5go0+XfQu3btmhksXbq0hkV5XpoTXSJAgAABAgQIECBAgAABAo1ZoCgC9DQFm0SAAAECBAgQIECAAAECBBqzQJMP0Pfbb7/MN023Vj2V55WP5l59u3UCBAgQIECAAAECBAgQINBYBJp8gD5gwIA44IAD4r777ovyAeMS7rJly7K8gw46KPr3799YvJ0HAQIECBAgQIAAAQIECBDYoECTHyQuteqSSy6JMWPGxMCBA7Pl0tLS+MlPfhIffPBBPPHEE5HmOC9Pr7zySvTu3TsOPPDAePnll8uzfRIgQIAAAQL1EBg6dGg99mp+u0ydOrX5NUqLCBAgQGCLC/wzct3ih264A44ePTrSKO5jx46NUaNGZRV37Ngxbr/99ujbt2/DHUhNBAgQIECAAAECBAgQIEAgJ4FmEaAnm1NPPTV7iv7WW2/FmjVromfPntGmTZsabOnJeXrCXlt66aWXaitiOwECBAgQIECAAAECBAgQaDCBZhOgJ5GSkpIsMG8wHRURIECAAAECBAgQIECAAIEtJNDkB4nbQk4OQ4AAAQIECBAgQIAAAQIEchUQoOfKq3ICBAgQIECAAAECBAgQIFCYgAC9MCelCBAgQIAAAQIECBAgQIBArgIC9Fx5VU6AAAECBAgQIECAAAECBAoTEKAX5qQUAQIECBAgQIAAAQwZe+UAACFJSURBVAIECBDIVUCAniuvygkQIECAAAECBAgQIECAQGECAvTCnJQiQIAAAQIECBAgQIAAAQK5CgjQc+VVOQECBAgQIECAAAECBAgQKExAgF6Yk1IECBAgQIAAAQIECBAgQCBXAQF6rrwqJ0CAAAECBAgQIECAAAEChQkI0AtzUooAAQIECBAgQIAAAQIECOQqIEDPlVflBAgQIECAAAECBAgQIECgMAEBemFOShEgQIAAAQIECBAgQIAAgVwFBOi58qqcAAECBAgQIECAAAECBAgUJiBAL8xJKQIECBAgQIAAAQIECBAgkKuAAD1XXpUTIECAAAECBAgQIECAAIHCBATohTkpRYAAAQIECBAgQIAAAQIEchUQoOfKq3ICBAgQIECAAAECBAgQIFCYgAC9MCelCBAgQIAAAQIECBAgQIBArgIC9Fx5VU6AAAECBAgQIECAAAECBAoTEKAX5qQUAQIECBAgQIAAAQIECBDIVUCAniuvygkQIECAAAECBAgQIECAQGECAvTCnJQiQIAAAQIECBAgQIAAAQK5CgjQc+VVOQECBAgQIECAAAECBAgQKExAgF6Yk1IECBAgQIAAAQIECBAgQCBXAQF6rrwqJ0CAAAECBAgQIECAAAEChQkI0AtzUooAAQIECBAgQIAAAQIECOQqIEDPlVflBAgQIECAAAECBAgQIECgMAEBemFOShEgQIAAAQIECBAgQIAAgVwFBOi58qqcAAECBAgQIECAAAECBAgUJiBAL8xJKQIECBAgQIAAAQIECBAgkKuAAD1XXpUTIECAAAECBAgQIECAAIHCBATohTkpRYAAAQIECBAgQIAAAQIEchUQoOfKq3ICBAgQIECAAAECBAgQIFCYgAC9MCelCBAgQIAAAQIECBAgQIBArgIC9Fx5VU6AAAECBAgQIECAAAECBAoTEKAX5qQUAQIECBAgQIAAAQIECBDIVUCAniuvygkQIECAAAECBAgQIECAQGECAvTCnJQiQIAAAQIECBAgQIAAAQK5CgjQc+VVOQECBAgQIECAAAECBAgQKExAgF6Yk1IECBAgQIAAAQIECBAgQCBXAQF6rrwqJ0CAAAECBAgQIECAAAEChQkI0AtzUooAAQIECBAgQIAAAQIECOQqIEDPlVflBAgQIECAAAECBAgQIECgMAEBemFOShEgQIAAAQIECBAgQIAAgVwFWuVau8oJEMhFYOjQobnU25QqnTp1alM6XedKgAABAgQIECBAoFYBT9BrJVKAAAECBAgQIECAAAECBAjkL+AJev7GjkCAAAECBAgUsYBeT/978fV8KuL/CTSdAIGCBTxBL5hKQQIECBAgQIAAAQIECBAgkJ+AAD0/WzUTIECAAAECBAgQIECAAIGCBQToBVMpSIAAAQIECBAgQIAAAQIE8hMQoOdnq2YCBAgQIECAAAECBAgQIFCwgAC9YCoFCRAgQIAAAQIECBAgQIBAfgIC9Pxs1UyAAAECBAgQIECAAAECBAoWEKAXTKUgAQIECBAgQIAAAQIECBDIT0CAnp+tmgkQIECAAAECBAgQIECAQMECrQou2QQKLliwIGbNmhXt27ePfv36ZZ91Pe2GqKOux1SeAAECBAgQIECAAAECBAg0myfo48aNix49esTIkSPjqKOOih122CGuvfbaOl3hhqijTgdUmAABAgQIECBAgAABAgQI/J9AswjQp02bFuPHj48RI0bESy+9FM8//3wMHjw4Lr744pgwYUJBF7sh6ijoQAoRIECAAAECBAgQIECAAIENCDT5Lu6ffPJJnHvuudGtW7eYNGlStGzZMmvmlClTYp999smeon/nO9+pyN+AQTREHRuqVx4BAgQIECBAgEDtAkOHDq29UBGUmDp1ahG0UhMJENiUQJN/gv7ss8/G/Pnz47TTTqsShLdu3TrGjBkT6Z3y2n7ZNUQdm0K2jQABAgQIECBAgAABAgQI1CbQ5J+gv/DCC1kbDznkkBptLc+bOXNmDB8+vMb28oyGqKO8Lp+bL+Bb9Kj1S6XNV1YDAQIECBAgQIAAAQKNTaDJB+hLlizJTDt16lTDdscdd8zyFi5cWGNb5Yy61rF69eq44447sio+/vjjylXVurxs2bJayzT3AmkAP4kAAQIECBAgQKBwgf/+7/8uvHAzLpl6yEr1FxCL/K9dY45HmnyAvnz58ky5c+fONe7U8gB95cqVNbZVzqhrHaWlpVEemKfluqTGfDPUpR15lq3tlYQ8j91U6mZU+5VixKh2gdpLuI8Y1S5Qewn3EaPaBWovITCt3UiJ2gXEIrUbbe0STT5Ab9u2bWb42Wef1bBcv359llc+cFyNAv+XUdc62rVrF5deemm293PPPbexauUTIECAAAECBAgQIECAAIGCBZr8IHFdu3bNGrt06dIajS7Pq+2booaoo8bBZRAgQIAAAQIECBAgQIAAgToIFEWAnqZg21QqJECvrY5N1W8bAQIECBAgQIAAAQIECBCoTaDJB+j77bdf1sY0VVr1VJ5XPpp79e3l6w1RR3ldPgkQIECAAAECBAgQIECAQH0ESsoGOavbKGf1OUrO+xx44IHx/vvvx7x586JDhw7Z0dIIhfvss0/ssssu8eKLL0arVpt+3b6+dQwbNiwefPDBSO+lSwQIECBAgAABAgQIECBAoL4CTf4Jemr4JZdcEosXL46BAwfGAw88EJMmTcqWP/jgg7jrrruqBOevvPJKlJSURO/evauY1aWOKjtaIUCAAAECBAgQIECAAAECDSCw6cfKDXCALVHF6NGjI43iPnbs2Bg1alR2yI4dO8btt98effv2LegUGqKOgg6kEAECBAgQIECAAAECBAgQ2IBAs+jiXt6u1Fv/rbfeijVr1kTPnj2jTZs25ZsK/qxrHbq4F0yrIAECBAgQIECAAAECBAhsQqBZPEEvb1/qup4C881JDVHH5hzfvgQIECBAgAABAgQIECBQnALN4h304rx0Wk2AAAECBAgQIECAAAECzUlAgN6crqa2ECBAgAABAgQIECBAgECTFRCgN9lL58QJECBAgAABAgQIECBAoDkJCNCb09XUFgIECBAgQIAAAQIECBBosgIC9CZ76Zw4AQIECBAgQIAAAQIECDQnAQF6c7qa2kKAAAECBAgQIECAAAECTVZAgN5kL50TJ0CAAAECBAgQIECAAIHmJNCs5kHfGhdm3bp18Y9//CPatWu3NQ7vmAQIECBAgAABAgQIECCQk0CnTp1yqnnD1ZaUlqUNb5JbiMDXvva12G677aJFi8bfGeGzzz6Lv/71r7H77rtHSUlJIc1ThsAGBebPnx9f+MIXomXLlhvcLpNAIQLp99Euu+wS22yzTSHFlSGwQYEFCxbE5z//+WjTps0Gt8skUIjAokWLomPHjh64FIKlzEYFFi9eHNtvv320b99+o2VsaHoCv/nNb6JVqy33XFuA3vTukXqf8ccffxwHH3xwzJkzR2BVb0U7JoH99tsv/vCHP8SOO+4IhEC9Bfr16xf3339/9qVhvSuxY9ELDB48OG644YY48MADi94CQP0FvvGNb8T3v//9OPzww+tfiT2LXuCb3/xmnHzyyXHMMccUvQWA+gs0/se+9W+bPQkQIECAAAECBAgQIECAQJMR8AS9yVyqhjnR1atXR9u2bRumMrUUrcCaNWt0Jy3aq99wDU/3UevWrb1y03CkRVnT2rVrs9ckvLpVlJe/wRqd7qPUhbUpvLLYYI1WUYMLpPsovf7nFcAGpy2qCgXoRXW5NZYAAQIECBAgQIAAAQIEGqvAlnvbvbEKFMl5vfvuu/Hyyy9Hjx494oADDiiSVmtmQwg8/PDDsddee0WvXr1qrS4NHpfeTT/11FNrLatAcQkUch/NmzcvGyOja9eu0bdvX4PHFdctUlBra/sdkwYe3NjYt926dduig/wU1CCFtorAxu6jpUuXxrJlyzZ5Tttuu2106dJlk2VsLA6Bjd1HlVvv7+/KGpYLFkijuEvNV6DsH5vSssEq0kj9FT+77bZb6S233NJ8G61lDSZwxx13ZPfN9ddfX2udZX/UlJYNHldaNqtBrWUVKC6B2u6jsqkqS0eMGFHxOyr9viqburL09ttvLy4ord2kQG2/Y5YsWVLlHqr8715aLvsCaJP121gcApu6j8oGidvkPZTuo7LBv4oDSis3KbCp+yjt6O/vTfLZWIuAJ+hlv22bc0ojSU6bNi1OOumkuOiii7L3PX/2s5/F2LFjIz1NOO6445pz87VtMwQeeeSR+Ld/+7eCavjwww9j9OjR8Ze//CWbdrCgnRQqCoFC7qNTTjkl+z11zjnnxNlnnx1///vfo+xLoTjvvPOy9/i+9a1vFYWVRm5coJDfMbNnz84qSKO677///jUqS1NoScUtUNt91L9//432wJgyZUq89dZbccQRRxQ3otZHbfdRIvL3txtlswRqCeBtbsICjz32WPZN8JAhQ2q0YtiwYdmTzr/97W81tskoboEPPvigtKyLenbvlM0rnH1u6gn6Qw89VFo2l3VWrmzAL0/Qi/v2qWh9offRiy++mN07X/rSlyr2TQtvv/12admgX6Vf+cpXquRbKT6BQn/H/PSnP83upWeeeab4kLS4VoFC76MNVfSnP/2ptGwAudKRI0eWfvbZZxsqIq9IBAq5j/z9XSQ3Q47NNM3aZn290bh3/uMf/5id4Pnnn1/jRM8444xI86I//vjjNbbJKG6Bsi9v4te//nWMGjUqyrombxLjt7/9bZxwwgmRRi1NT0oLeU99kxXa2GwECr2P0vucl19+efzkJz+p0vY0Xkb6Sb0ypOIVqMvvmPQEPY3knsYvkAhUFqjLfVR5v7Sc/n1Lc1vvsMMOcffdd5t1ojpQEa0Xeh/5+7uIboqcmipAzwm2MVRb9nQ8O43dd9+9xul06NAhy5s5c2aNbTKKW+Dggw/Ouhvff//98bnPfW6TGGlKmssuuyxef/31KHuysMmyNhaXQKH30Re/+MUYP358pG7JldOsWbMiDcBz1FFHVc62XGQCdfkdkwL0vffeOwuo7r333rjxxhvjySefjFWrVhWZmuZWF6jLfVR93yuuuCL+/Oc/x3XXXRedOnWqvtl6EQkUeh/5+7uIboqcmuod9JxgG0O13bt3z04jvTNVfeT2Rx99NNuWRiyVCFQWuPXWWyuvbnK57PWJSD8SgeoCdbmPyvct6y0W99xzTxZUpd49qUdG+qNYKl6BQn/HfPLJJ9kXhZ///OeznhcrVqyoQEuzUPzqV7+KQw45pCLPQnEJFHofVVd5880349prr430RWJ6ii4Vt0Ch95G/v4v7PmmI1nuC3hCKjbSOgQMHZmd21VVXRfrjpTytXLkyyt6hyVYr/xFTvt0nAQIEtoZAGhzuzDPPjN/85jeRfjelXhlpMEuJQG0Cr7zySpS9G5wN3pRemZgzZ0689tprcckll0TZeAZRNktA+EK6NkXbqwvceeedsX79+vjOd76ja3t1HOsbFfD390ZpbChQwBP0AqGaYrEBAwZkI2unrn59+vTJllM70vvFBx10UEydOjXat2/fFJvmnAkQaIYCaZTtNI912XRZ8Ytf/CKuueaamDx5cjz//PNmB2iG17shm7THHntE+rfuC1/4Qhx++OEVVV999dVZgJWegt5www2RvrCWCBQi8Omnn8b/+3//L7bffvv4l3/5l0J2UYZAJuDvbzfC5goI0DdXsJHvn7qLpoGWUnfTK6+8Mrp27ZoN/pUGjttzzz2zQU8aeROcHgECRSJQNvd5FmClIKtsVPcomx89HnjggazL+ze+8Y0iUdDM+gjstNNOkabr21BKwVUK0NO4BhKBQgXStGrvvfde9vQ8BekSgboI+Pu7LlrKVhfQxb26SDNb32abbeI//uM/sm5/ixcvjoULF8ZNN90UCxYsyFqa3s2TCBAg0BgFyuc/N9tEY7w6Teec0nvpKS1fvrzpnLQz3eoCqXt7Sql7u0SgrgL+/q6rmPKVBQTolTWa2XLqKpq6iaZ381Lq0qVLRQsffvjhbPmrX/1qRZ4FAgQIbGmBNAhc6tr+u9/9rsahW7T433+itttuuxrbZBCoLJBGbN9nn32ybu6V89Py3Llzs6y0XSJQiEB67/yZZ56Jnj17mj60EDBlqgj4+7sKh5V6CAjQ64HWVHZJTwvOOeec+N73vlfllNMvjl/+8pfxta99LQToVWisECCwhQX23Xff+Oijj2LChAk1jnzzzTdneaZaq0Ejo5rAbrvtlo3innqMpdkAylNaTu+hp+Q94nIVn7UJvPPOO7FmzZrYf//9aytqO4EaAv7+rkEio44C3kGvI1hTKp7+YUlzCz/99NNx9tlnxwknnBDpH53rr78+G400DX5SUlLSlJrkXAkQaGYCxx57bBxzzDGRevWkLw3TVEbbbrtt3Hbbbdm756NGjYqvf/3rzazVmtPQAscdd1ykkZN///vfx6BBg+Kss87KBhb8+c9/HtOmTcv+Dezfv39DH1Z9zVQgzQKQkgC9mV7gnJvl7++cgYugegF6M7/IaVTbf/3Xf4277ror+0nNTb840hP0XXbZpZm3XvMIEGjsAulLwjSt2mWXXZYNZpmCqZRSkP7jH/84Lr744sbeBOfXCARatmwZDz74YFx66aWR3h1O3ZNT6tSpUzZA3IUXXpit+w+BQgQE6IUoKbMpAX9/b0rHttoESsq6f/2zL1htpW1vsgJp/tc0F2wa6TaNkOzJeZO9lE6cQLMVWLVqVcybNy8LztMsEynokgjUVWD16tXxxhtvZNNjde/eva67K0+AAIEGE/D3d4NRFlVFAvSiutwaS4AAAQIECBAgQIAAAQKNVcAgcY31yjgvAgQIECBAgAABAgQIECgqAQF6UV1ujSVAgAABAgQIECBAgACBxiogQG+sV8Z5ESBAgAABAgQIECBAgEBRCQjQi+pyaywBAgQIECBAgAABAgQINFYBAXpjvTLOiwABAgQIECBAgAABAgSKSkCAXlSXW2MJECBAgAABAgQIECBAoLEKCNAb65VxXgQIECBAgAABAgQIECBQVAIC9KK63BpLgAABAgQIECBAgAABAo1VQIDeWK+M8yJAgAABAgQIECBAgACBohIQoBfV5dZYAgQIECBAgAABAgQIEGisAgL0xnplnBcBAgQIECBAgAABAgQIFJWAAL2oLrfGEiBAgAABAgQIECBAgEBjFRCgN9Yr47wIECBAgAABAgQIECBAoKgEBOhFdbk1lgABAgTyFnjhhRdi8uTJsXbt2rwP1WD1f/rpp9k5p3PfGqkpmm0NJ8ckQIAAgeYvIEBv/tdYCwkQIEBgCwpcd911ccIJJ8SyZcvqddQULKc67rvvvnrtX5+dPv744+ycr7nmmvrsXqd9Xn311Tj11FOr7LO5ZlUqs0KAAAECBJqwQKsmfO5OnQABAgQINDqBvn37xqpVq6J169b1Orf7778/LrroovjFL35Rr/0b+07HH398rF69usppbq5ZlcqsECBAgACBJiwgQG/CF8+pEyBAgEDjE7jkkksa30k18jNi1sgvkNMjQIAAgS0mIEDfYtQORIAAAQLFIPDss8/G/Pnz4+STT462bdvGhx9+GFOmTInDDjssdt5553jyySfjxRdfjC5dusSQIUPiwAMPrGD5n//5n0g/Kc2YMSNatWoV6Ylz6vb+2GOPxeGHHx6pO/qkSZOiT58+ccwxx0T79u2z8q+99lpMmzYt3n333ejRo0f0798/DjrooGxb9f+8/vrr8cQTT8RHH30UgwYNigMOOKB6kXj66adj4cKFcdppp0XLli0rtqd36++9997YbbfdYuDAgRX5aWHdunUxc+bM+P3vf5/Vndp20kknxTbbbBPLly/P3nNPn6k999xzT8V5VjcrrzQd66GHHoo///nPWVbv3r3j2GOPjW233ba8SJ18006lpaUxderU+OMf/5hZ7rvvvjFgwIDYe++9K+q0QIAAAQIEtppA2T9UEgECBAgQINBAAieeeGJp2T/qpe+9915W48svv5ytX3755aVf/OIXs+WygDX7LAvAS2+//faKI48ZMybLT/unn5KSktK5c+eWlgW92fqll15ausMOO1SUee6557J9v/e975WWBdFZ+V122SX7bNGiRWnKLwuGK+pPC9/5zney/dM5dO7cOVv+1re+lX2WvTtfUfboo4/O8j755JOKvLSwdOnSLP+4446rkl/2RUTp/vvvn23bfvvtS7fbbrtseb/99ist6/JfWvalQHZeldt2yimnZHVUN0uZL730UmnPnj0r6uvQoUO2nPKef/75imPXxbcs4C8dOnRoVk/y2mmnnbLlstcRSidMmFBRpwUCBAgQILC1BAwSV/aXgkSAAAECBPIWuOqqq6JXr17Zk/GyoDcefvjhKAui4wc/+EGk9ZR+/etfx69+9atsOb2D/tlnn8U+++yTraf/XHvttTF48OB4/PHH4/rrr8+eqN99991x4403Zvl///vfY9GiRVH25UAMHz48y7/55psr9k913nrrrdlT8fT0/P3338+elKdR5zcnrVmzJkaOHBlvvPFGTJw4MRsgL53DT3/60/jLX/4S//mf/xl77bVX1p6yADt23XXXbDk9id9QSu/wjx49OlJ7Upl0runnwQcfjCVLlmS9CtKT+MqpEN9UV3p6fuGFF2ZP3lNdqedBp06dsrx0DIkAAQIECGxNAQH61tR3bAIECBAoGoHUJTwF36mre+q6/vWvfz37SV3g33zzzYIcUrf4FMQPGzYsLrjggqy7dvrccccds1Hf0/aUyp6MZ4Ft6lJ/xRVXxMqVK7P8K6+8Mutaf8cdd1R0Ez/qqKNi3Lhx2fb6/mf69OlZ1/x//dd/zYL/sif/0a5du2ywu9SFPnVhT93aC0233XZbzJs3Ly6++OIoe8qefZGR6kyj46f2pC8h0pcSlVMhvv+/vfsHiWWH4jieBw9sRGysRBRsLP0DYmOpnZ22ioJYKIKlraL2IliJ2glWYicoloKIjdgpFgpaWIitknd+gQyz487evdd3s4XfFHfHSTbZ+VT3TJIT9amirQE2yx+ubVWD29vbCy889KKBggACCCCAQCMFCNAbqc/YCCCAAAI/RmBgYOBLZveurq7w/MXZ4DIU7cFuamrKqrXfXLO+mi23pe/ZfV1ob7r2r2vPui2Td6+vr+7x8dGNjo6G4DnfWPvlv1Our6/D1xVA54uC6tPTU3dychL2oefral3bsvVQXTyOTTe1J15Fe93zpR5fvYxQ0Wy/LfV3x8fH4eWFTObn58PLi3yfXCOAAAIIIJBagAA9tTjjIYAAAgj8SAHNZheLksipaCl7PUXJ3/Ilzgh3dnbmb2fX8b5m6GOitfb29qw+Xthe7IrAP96v9zMG1B0dHfV+pWY7PZeC+2r96bdqdr646qAeX83ma5m/kt5tb2+HQF3L25V4zva11/xNVCKAAAIIIJBCgAA9hTJjIIAAAgj8eAHtN/9uKZ6tHjO4xyXsxf7f39/DLb0IUCCqUq2tJcIpfUmgunyJfebvKUu7StxLn6/7k2s9l8bVXvRiUWZ3naMeX27E+np9LSGee35+DrPnCwsLYT+89vQr6/35+Xnsjk8EEEAAAQQaIvD9/y005GczKAIIIIAAAggo8ZrK7e1tVYx4XzPp2mutmWctdy8WLZUv7hGPAbCC4Xy5u7vL/xmulfhNpVqdEttp6fv9/X1oU88/tZ5Lv1/Be1wdUE9/sY1m3XW8nJ5Ns+aWuT3MxCuZXTw+LrblEwEEEEAAgUYIEKA3Qp0xEUAAAQQQKBGIs9HVZrqLX1FSuKGhobDHu7gn++bmJpy/rmXx2ruuZd06R13nm9sRZhVdbW5uVvytP5RoTuXs7Cx8xn+UUE0lP7OupHVakq6AN1+0t35tbS2MGZfW6/l+9Wx2hFvoZmNjo2Ic3VxfXw912l//u0XZ27VfXzPm+dLf3x/+zJ+vnq/nGgEEEEAAgVQC/6YaiHEQQAABBBBA4NcC2mOtsrW15Z6enpyWYdcqaqfM8Dp+zc5ad729vU7B+crKSkgcpyPUFDyr6Ei2wcFBNzIy4lZXV51mvo+Ojtzu7m4I4PPjTExMuJ2dnZA8TXvCFWDraDjt1Y4Z0GN7BbgzMzOhvbLTaxm5jkjTfm8F6TrqLSa30/Pp6LXp6Wk3PDwcvhf7iZ92BntIcKffPjY2FvrTEvb9/X2ne7Ozs25ycjI2r/tzcXExvLSYm5tzU1NTYXwluNNz6vfZOfR190VDBBBAAAEE/oqAvQGnIIAAAggggMD/JDA+Pq5N297OAQ89WgK18LcF2l9GsIA61NkxZFmdLTX3llXd21Fsoe7w8NDb7Hi4XlpaytrlL2xG3FsW89BGY9tSdm8Zy/3V1VW+Wbi+vLz0lizN20x2aG+z8N72Xvvm5mZvS9Er2tvMurdgPLSzIN/39fV5WybuLXmbt1nuirYfHx/egn5v+8ez39Ha2urVR75oLAvSQxs7Fz5UFc1005ace/nk++vu7vbLy8vekuplXf6u78HBgbfl8dlv1HPZ8n9/cXGR9ckFAggggAACjRL4RwP/lcifThFAAAEEEEDgjwWUIE0J2eKMej0dabb64eHB9fT0fDnSrfj9t7c3Zy8RXNzvXayPfyvDvGbQ29rasmXvsa7ap9prv/nn56ezgDqc+V6t3cvLS5jhj3vdq7XRPf03Rf0pQV61rO5l36t1X79RZ6nr+bWKoKWlpVZz6hBAAAEEEEgmQICejJqBEEAAAQQQQAABBBBAAAEEECgXIElcuQ01CCCAAAIIIIAAAggggAACCCQTIEBPRs1ACCCAAAIIIIAAAggggAACCJQLEKCX21CDAAIIIIAAAggggAACCCCAQDIBAvRk1AyEAAIIIIAAAggggAACCCCAQLkAAXq5DTUIIIAAAggggAACCCCAAAIIJBMgQE9GzUAIIIAAAggggAACCCCAAAIIlAsQoJfbUIMAAggggAACCCCAAAIIIIBAMgEC9GTUDIQAAggggAACCCCAAAIIIIBAuQABerkNNQgggAACCCCAAAIIIIAAAggkEyBAT0bNQAgggAACCCCAAAIIIIAAAgiUCxCgl9tQgwACCCCAAAIIIIAAAggggEAyAQL0ZNQMhAACCCCAAAIIIIAAAggggEC5AAF6uQ01CCCAAAIIIIAAAggggAACCCQTIEBPRs1ACCCAAAIIIIAAAggggAACCJQLEKCX21CDAAIIIIAAAggggAACCCCAQDIBAvRk1AyEAAIIIIAAAggggAACCCCAQLkAAXq5DTUIIIAAAggggAACCCCAAAIIJBMgQE9GzUAIIIAAAggggAACCCCAAAIIlAv8B8WEXYegXqZ3AAAAAElFTkSuQmCC\n"
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"%%R -w 1000 -h 500 -u px -i df,min_value,max_value # this sets the size of the plot...otherwise, it will go off the page\n",
"require(ggplot2)\n",
"library(ggplot2)\n",
"\n",
"array_of_transitions <- as.character(min_value:max_value)\n",
"\n",
"df$number_transitions = as.character(df$number_transitions)\n",
"df$number_transitions = factor(df$number_transitions, levels=array_of_transitions)\n",
"\n",
"p <- ggplot(data=df, aes(x=number_transitions, y=proportion_trees)) + \n",
" geom_col()+ #aes(y = (..count..)/sum(..count..))\n",
" labs(x=\"introductions\",y=\"proportion of trees\")+\n",
" #ggtitle(\"Re estimates\") + \n",
" theme(plot.title = element_text(size=20, hjust=0.5))+\n",
" scale_x_discrete(breaks=array_of_transitions)+\n",
" scale_x_discrete(breaks=c(\"9\",\"\",\"11\",\"\",\"13\",\"\",\"15\",\"\",\"17\",\"\",\"19\",\"\"))+\n",
" scale_y_continuous(limits=c(0,0.5), breaks=seq(0,0.5,0.1))+\n",
" #scale_color_manual(values=c(\"#4784C7\", \"#E39B39\"))+\n",
" #scale_fill_manual(values=c(\"#4784C7\", \"#E39B39\"))+\n",
" theme(panel.grid.major=element_line(colour=NA,size=NA))+\n",
" theme(panel.grid.minor=element_line(colour=NA,size=NA))+ \n",
" theme(strip.background = element_rect(colour=NA, fill=NA))+\n",
" theme(axis.line.x=element_line(colour=\"black\"))+\n",
" theme(axis.line.y=element_line(colour=\"black\"))+\n",
" theme(strip.text.x=element_text(size=16))+\n",
" theme(axis.title.y=element_text(size=20, vjust=4))+\n",
" theme(axis.title.x=element_text(size=20, vjust=-3))+\n",
" theme(axis.text=element_text(size=20, colour=\"black\"))+\n",
" theme(axis.text.x=element_text(size=20))+\n",
" theme(legend.title=element_blank())+\n",
" theme(panel.margin=unit(1, \"lines\"))+\n",
" theme(plot.margin=unit(c(1,1,1,1),\"cm\"))+\n",
" theme(legend.key.size=unit(0.7, \"cm\"))+\n",
" theme(panel.background=element_rect(fill=NA))+\n",
" theme(legend.key=element_rect(fill=NA))\n",
"\n",
"ggsave(\"2020-12-14-DTA-skygrid-25-WA-introductions.pdf\", p, width = 5, height = 3.5, path=\"/Users/lmoncla/Documents/Mumps/paper-and-figure-drafts/individual-PDFs/\")\n",
"p"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": []
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": []
}
],
"metadata": {
"kernelspec": {
"display_name": "LHM-basics (python3)",
"language": "python",
"name": "lhm-basics"
},
"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.7.3"
}
},
"nbformat": 4,
"nbformat_minor": 2
}
