diff --git a/data/processed/NPP_NorthAmerica_EqualArea_100km.tif b/data/processed/NPP_NorthAmerica_EqualArea_100km.tif new file mode 100644 index 0000000..b5c6f6d Binary files /dev/null and b/data/processed/NPP_NorthAmerica_EqualArea_100km.tif differ diff --git a/data/processed/Tavg_NorthAmerica_EqualArea_100km_month_01.tif b/data/processed/Tavg_NorthAmerica_EqualArea_100km_month_01.tif new file mode 100644 index 0000000..1538462 Binary files /dev/null and b/data/processed/Tavg_NorthAmerica_EqualArea_100km_month_01.tif differ diff --git a/data/processed/Tavg_NorthAmerica_EqualArea_100km_month_02.tif b/data/processed/Tavg_NorthAmerica_EqualArea_100km_month_02.tif new file mode 100644 index 0000000..5ca5289 Binary files /dev/null and b/data/processed/Tavg_NorthAmerica_EqualArea_100km_month_02.tif differ diff --git a/data/processed/Tavg_NorthAmerica_EqualArea_100km_month_03.tif b/data/processed/Tavg_NorthAmerica_EqualArea_100km_month_03.tif new file mode 100644 index 0000000..4df88c2 Binary files /dev/null and b/data/processed/Tavg_NorthAmerica_EqualArea_100km_month_03.tif differ diff --git a/data/processed/Tavg_NorthAmerica_EqualArea_100km_month_04.tif b/data/processed/Tavg_NorthAmerica_EqualArea_100km_month_04.tif new file mode 100644 index 0000000..bb86834 Binary files /dev/null and b/data/processed/Tavg_NorthAmerica_EqualArea_100km_month_04.tif differ diff --git a/data/processed/Tavg_NorthAmerica_EqualArea_100km_month_05.tif b/data/processed/Tavg_NorthAmerica_EqualArea_100km_month_05.tif new file mode 100644 index 0000000..854f88e Binary files /dev/null and b/data/processed/Tavg_NorthAmerica_EqualArea_100km_month_05.tif differ diff --git a/data/processed/Tavg_NorthAmerica_EqualArea_100km_month_06.tif b/data/processed/Tavg_NorthAmerica_EqualArea_100km_month_06.tif new file mode 100644 index 0000000..51e8093 Binary files /dev/null and b/data/processed/Tavg_NorthAmerica_EqualArea_100km_month_06.tif differ diff --git a/data/processed/Tavg_NorthAmerica_EqualArea_100km_month_07.tif b/data/processed/Tavg_NorthAmerica_EqualArea_100km_month_07.tif new file mode 100644 index 0000000..a4c2e29 Binary files /dev/null and b/data/processed/Tavg_NorthAmerica_EqualArea_100km_month_07.tif differ diff --git a/data/processed/Tavg_NorthAmerica_EqualArea_100km_month_08.tif b/data/processed/Tavg_NorthAmerica_EqualArea_100km_month_08.tif new file mode 100644 index 0000000..e17d230 Binary files /dev/null and b/data/processed/Tavg_NorthAmerica_EqualArea_100km_month_08.tif differ diff --git a/data/processed/Tavg_NorthAmerica_EqualArea_100km_month_09.tif b/data/processed/Tavg_NorthAmerica_EqualArea_100km_month_09.tif new file mode 100644 index 0000000..1da9611 Binary files /dev/null and b/data/processed/Tavg_NorthAmerica_EqualArea_100km_month_09.tif differ diff --git a/data/processed/Tavg_NorthAmerica_EqualArea_100km_month_10.tif b/data/processed/Tavg_NorthAmerica_EqualArea_100km_month_10.tif new file mode 100644 index 0000000..ad46a48 Binary files /dev/null and b/data/processed/Tavg_NorthAmerica_EqualArea_100km_month_10.tif differ diff --git a/data/processed/Tavg_NorthAmerica_EqualArea_100km_month_11.tif b/data/processed/Tavg_NorthAmerica_EqualArea_100km_month_11.tif new file mode 100644 index 0000000..61f7f10 Binary files /dev/null and b/data/processed/Tavg_NorthAmerica_EqualArea_100km_month_11.tif differ diff --git a/data/processed/Tavg_NorthAmerica_EqualArea_100km_month_12.tif b/data/processed/Tavg_NorthAmerica_EqualArea_100km_month_12.tif new file mode 100644 index 0000000..0cb03ab Binary files /dev/null and b/data/processed/Tavg_NorthAmerica_EqualArea_100km_month_12.tif differ diff --git a/data/raw/MOD17A3H_Y_NPP_2025-01-01_rgb_720x360.TIFF b/data/raw/MOD17A3H_Y_NPP_2025-01-01_rgb_720x360.TIFF new file mode 100644 index 0000000..6222b9c Binary files /dev/null and b/data/raw/MOD17A3H_Y_NPP_2025-01-01_rgb_720x360.TIFF differ diff --git a/data/raw/wc2.1_10m_tavg/wc2.1_10m_tavg_01.tif b/data/raw/wc2.1_10m_tavg/wc2.1_10m_tavg_01.tif new file mode 100644 index 0000000..3da58c1 Binary files /dev/null and b/data/raw/wc2.1_10m_tavg/wc2.1_10m_tavg_01.tif differ diff --git a/data/raw/wc2.1_10m_tavg/wc2.1_10m_tavg_02.tif b/data/raw/wc2.1_10m_tavg/wc2.1_10m_tavg_02.tif new file mode 100644 index 0000000..33ee99a Binary files /dev/null and b/data/raw/wc2.1_10m_tavg/wc2.1_10m_tavg_02.tif differ diff --git a/data/raw/wc2.1_10m_tavg/wc2.1_10m_tavg_03.tif b/data/raw/wc2.1_10m_tavg/wc2.1_10m_tavg_03.tif new file mode 100644 index 0000000..35cf8bc Binary files /dev/null and b/data/raw/wc2.1_10m_tavg/wc2.1_10m_tavg_03.tif differ diff --git a/data/raw/wc2.1_10m_tavg/wc2.1_10m_tavg_04.tif b/data/raw/wc2.1_10m_tavg/wc2.1_10m_tavg_04.tif new file mode 100644 index 0000000..f8e7d56 Binary files /dev/null and b/data/raw/wc2.1_10m_tavg/wc2.1_10m_tavg_04.tif differ diff --git a/data/raw/wc2.1_10m_tavg/wc2.1_10m_tavg_05.tif b/data/raw/wc2.1_10m_tavg/wc2.1_10m_tavg_05.tif new file mode 100644 index 0000000..2b552fd Binary files /dev/null and b/data/raw/wc2.1_10m_tavg/wc2.1_10m_tavg_05.tif differ diff --git a/data/raw/wc2.1_10m_tavg/wc2.1_10m_tavg_06.tif b/data/raw/wc2.1_10m_tavg/wc2.1_10m_tavg_06.tif new file mode 100644 index 0000000..0a2b3bc Binary files /dev/null and b/data/raw/wc2.1_10m_tavg/wc2.1_10m_tavg_06.tif differ diff --git a/data/raw/wc2.1_10m_tavg/wc2.1_10m_tavg_07.tif b/data/raw/wc2.1_10m_tavg/wc2.1_10m_tavg_07.tif new file mode 100644 index 0000000..a610d14 Binary files /dev/null and b/data/raw/wc2.1_10m_tavg/wc2.1_10m_tavg_07.tif differ diff --git a/data/raw/wc2.1_10m_tavg/wc2.1_10m_tavg_08.tif b/data/raw/wc2.1_10m_tavg/wc2.1_10m_tavg_08.tif new file mode 100644 index 0000000..5218640 Binary files /dev/null and b/data/raw/wc2.1_10m_tavg/wc2.1_10m_tavg_08.tif differ diff --git a/data/raw/wc2.1_10m_tavg/wc2.1_10m_tavg_09.tif b/data/raw/wc2.1_10m_tavg/wc2.1_10m_tavg_09.tif new file mode 100644 index 0000000..e4f1ce0 Binary files /dev/null and b/data/raw/wc2.1_10m_tavg/wc2.1_10m_tavg_09.tif differ diff --git a/data/raw/wc2.1_10m_tavg/wc2.1_10m_tavg_10.tif b/data/raw/wc2.1_10m_tavg/wc2.1_10m_tavg_10.tif new file mode 100644 index 0000000..ff20e2e Binary files /dev/null and b/data/raw/wc2.1_10m_tavg/wc2.1_10m_tavg_10.tif differ diff --git a/data/raw/wc2.1_10m_tavg/wc2.1_10m_tavg_11.tif b/data/raw/wc2.1_10m_tavg/wc2.1_10m_tavg_11.tif new file mode 100644 index 0000000..85ae418 Binary files /dev/null and b/data/raw/wc2.1_10m_tavg/wc2.1_10m_tavg_11.tif differ diff --git a/data/raw/wc2.1_10m_tavg/wc2.1_10m_tavg_12.tif b/data/raw/wc2.1_10m_tavg/wc2.1_10m_tavg_12.tif new file mode 100644 index 0000000..a602585 Binary files /dev/null and b/data/raw/wc2.1_10m_tavg/wc2.1_10m_tavg_12.tif differ diff --git a/experiments/notebooks/GEM.ipynb b/experiments/notebooks/GEM.ipynb new file mode 100644 index 0000000..1bdf1da --- /dev/null +++ b/experiments/notebooks/GEM.ipynb @@ -0,0 +1,254 @@ +{ + "cells": [ + { + "cell_type": "code", + "execution_count": 4, + "id": "7790bc03", + "metadata": {}, + "outputs": [], + "source": [ + "%load_ext autoreload\n", + "%autoreload 2\n", + "\n", + "#set the path to the root of the project dynamically\n", + "import sys\n", + "import os\n", + "sys.path.append(os.path.abspath(os.path.join(os.getcwd(), os.pardir, os.pardir)))\n", + "\n" + ] + }, + { + "cell_type": "markdown", + "id": "405058af", + "metadata": {}, + "source": [ + "# Data fetching" + ] + }, + { + "cell_type": "markdown", + "id": "fc462839", + "metadata": {}, + "source": [ + "You just run the script in `./src/utils/data_fetching.py`" + ] + }, + { + "cell_type": "markdown", + "id": "731ae5e6", + "metadata": {}, + "source": [ + "# Engine V2" + ] + }, + { + "cell_type": "code", + "execution_count": 10, + "id": "4ca3786d", + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Environment initialized with layers: ['carrying_capacity', 'npp', 'temperature']\n", + "Grid shape: (75, 120)\n" + ] + } + ], + "source": [ + "from src.gem.engine.ecosystem_engine import EcosystemEngine\n", + "from src.gem.engine.environment_state import EnvironmentState\n", + "from src.gem.engine.ecosystem_grid_state import EcosystemGridState\n", + "from src.gem.engine.species_registry import SpeciesRegistry\n", + "\n", + "from src.gem.engine.processes import apply_vegetation_growth, apply_atn_step, apply_dispersal\n", + "import numpy as np\n", + "import rasterio\n", + "\n", + "processed_dir = './data/processed'\n", + "npp_path = os.path.join(processed_dir, 'NPP_NorthAmerica_EqualArea_100km.tif')\n", + "temp_path = os.path.join(processed_dir, 'Tavg_NorthAmerica_EqualArea_100km_month_01.tif') # January example\n", + "\n", + "# 2. Load the actual raster data\n", + "with rasterio.open(npp_path) as src:\n", + " # src.read(1) returns a (75, 120) numpy array\n", + " npp_data = src.read(1)\n", + " \n", + "with rasterio.open(temp_path) as src:\n", + " temp_data = src.read(1)\n", + " \n", + " \n", + "\n", + "world_shape = (75, 120) # (12000 km x 7500 km with 100 km cell size)\n", + "env = EnvironmentState(world_shape)\n", + "env.add_layer(\"carrying_capacity\", np.full(world_shape, 100000.0))\n", + "env.add_layer(\"npp\", npp_data)\n", + "env.add_layer(\"temperature\", temp_data)\n", + "\n", + "print(f\"Environment initialized with layers: {list(env.layers.keys())}\")\n", + "print(f\"Grid shape: {env.shape}\")\n", + "\n", + "# Create 3 species: 1 plant, 2 animals\n", + "registry = SpeciesRegistry(num_species=20)\n", + "registry.add_species_to_group(\n", + " group_name=\"deciduous_trees\", \n", + " species_indices=list(range(0, 5)), \n", + " params={\n", + " \"base_growth_rate\": 0.15, \n", + " \"mortality_rate\": 0.05,\n", + " \"max_dispersal_rate\": 0.01\n", + " }\n", + ")\n", + "\n", + "registry.add_species_to_group(\n", + " group_name=\"coniferous_trees\", \n", + " species_indices=list(range(10, 15)), \n", + " params={\n", + " \"base_growth_rate\": 0.08, \n", + " \"mortality_rate\": 0.02,\n", + " \"max_dispersal_rate\": 0.02\n", + " }\n", + ")\n", + "\n", + "registry.add_species_to_group(\n", + " group_name=\"herbs\", \n", + " species_indices=list(range(5, 10)), \n", + " params={\n", + " \"base_growth_rate\": 0.2, \n", + " \"mortality_rate\": 0.1,\n", + " \"max_dispersal_rate\": 0.05\n", + " }\n", + ")\n", + "registry.add_species_to_group(\n", + " group_name=\"herbivores\", \n", + " species_indices=list(range(15, 18)), \n", + " params={\n", + " \"base_growth_rate\": 0.1, \n", + " \"mortality_rate\": 0.05,\n", + " \"max_dispersal_rate\": 0.1\n", + " }\n", + ")\n", + "registry.add_species_to_group(\n", + " group_name=\"carnivores\", \n", + " species_indices=list(range(18, 20)), \n", + " params={\n", + " \"base_growth_rate\": 0.05, \n", + " \"mortality_rate\": 0.02,\n", + " \"max_dispersal_rate\": 0.2\n", + " }\n", + ")\n", + "\n", + "#bigger groups for more interesting dynamics\n", + "registry.add_species_to_group(group_name=\"vertebrates\", species_indices=list(range(15, 20))) # All vertebrates for process targeting\n", + "registry.add_species_to_group(group_name=\"plants\", species_indices=list(range(0, 15))) # All plants for process targeting\n", + "\n", + "#add feeding links (plants -> herbivores -> carnivores)\n", + "registry.add_feeding_link(\"plants\", \"herbivores\")\n", + "registry.add_feeding_link(\"herbivores\", \"carnivores\")\n", + "\n", + "\n", + "#Grid state with 20 species (15 plants, 5 animals)\n", + "grid = EcosystemGridState(world_shape, registry)\n", + "\n", + "#ADDING LAYERS\n", + "grid.add_layer(\"metabolic_loss\")\n", + "grid.add_layer(\"net_growth_rate\")\n", + "\n", + "#ADDING DELTA LAYERS\n", + "grid.add_delta_layer(\"vegetation_delta\", source_layer=\"biomass\") # Register a delta for net growth rates\n", + "\n", + "# SEED INITIAL LAYERS\n", + "grid.layers[\"biomass\"][:, :, registry.get_group_indices(\"plants\")] = 50000.0 # Plants\n", + "grid.layers[\"biomass\"][:, :, registry.get_group_indices(\"animals\")] = 5000.0 # Animals\n", + "\n", + "# --- 2. Build the Engine ---\n", + "model = EcosystemEngine(grid, env)\n", + "\n", + "# The order here defines the pipeline sequence\n", + "model.add_process(apply_vegetation_growth)\n", + "model.add_process(apply_atn_step)\n", + "model.add_process(apply_dispersal)\n", + "\n", + "# --- 3. Run the Simulation ---\n", + "num_steps = 150\n", + "vertebrate_biomass_history = []\n", + "plant_biomass_history = []\n", + "\n", + "for tick in range(num_steps):\n", + " model.step()\n", + " \n", + " # Track Totals using the new matrix views\n", + " # np.sum() cleanly aggregates the whole grid\n", + " total_plant = np.sum(grid.get_layer_view(\"biomass\", \"plants\"))\n", + " total_vertebrates = np.sum(grid.get_layer_view(\"biomass\", \"vertebrates\"))\n", + " \n", + " plant_biomass_history.append(total_plant)\n", + " vertebrate_biomass_history.append(total_vertebrates)\n", + " \n" + ] + }, + { + "cell_type": "code", + "execution_count": 31, + "id": "57a478cf", + "metadata": {}, + "outputs": [ + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAABKUAAAHqCAYAAADVi/1VAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjYuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/P9b71AAAACXBIWXMAAA9hAAAPYQGoP6dpAACKY0lEQVR4nOzdeVyU5frH8e+wDYssIgqiqLjkvqXHNUvLVNzLyrQ0TU1PmRqVZptonswWs3JrcckW89dJzcpMyzRLTHGr1LQSxQXcWQRlGZ7fH8gcRxZBgQHm83695qXzzP08c90PKNdc3IvJMAxDAAAAAAAAQAlysncAAAAAAAAAcDwUpQAAAAAAAFDiKEoBAAAAAACgxFGUAgAAAAAAQImjKAUAAAAAAIASR1EKAAAAAAAAJY6iFAAAAAAAAEocRSkAAAAAAACUOIpSAAAAAAAAKHEUpYASsGTJEplMJuvDxcVF1atX1/Dhw3X8+HFru40bN8pkMmnjxo3FGs+8efO0ZMmSArevVauWTfzu7u6qW7euwsPDdebMGZu2ERERMplMRRxx2XH119rd3V1BQUHq0qWLZsyYoVOnTtk7xCI1bNgw1apVy95hAABKsbvuukseHh6Kj4/Ps80DDzwgV1dXnTx58obf78SJE4qIiNDu3buv+xrZP8+joqJuOJ6CWLNmjSIiIor8ulfnJSaTSZUrV1bnzp319ddf52hvMpmKJY6y4sqc18nJSb6+vmrYsKGGDh2qdevW2Tu8InX48GGZTKZCfSYAigNFKaAELV68WJGRkVq/fr1GjRqlZcuWqVOnTkpOTi7ROApblJKkjh07KjIyUpGRkfr22281evRovfvuu+rRo4dNu5EjRyoyMrIIoy2brvxaz507Vy1atNDMmTPVsGFDff/99/YOr8i88MILWrlypb3DAACUYiNGjNClS5f06aef5vp6QkKCVq5cqd69eyswMPCG3+/EiROaOnXqDRWlStqaNWs0derUYrt+dl6yZcsWvffee3J2dlafPn301Vdf2bSLjIzUyJEjiy2OsiA7592yZYu++OILjR07VtHR0erevbvuuecepaen2zvEIlG1alVFRkaqV69e9g4FDs7F3gEAjqRJkyZq3bq1JKlLly6yWCx66aWXtGrVKj3wwAN2ji5/fn5+ateunfV5ly5dlJSUpJdeekkHDx7UTTfdJEmqXr26qlevbq8wS40rv9aSNGDAAD3xxBO65ZZbdPfdd+uvv/4qksTb3urUqWPvEAAApVxYWJiCg4O1aNEiPfroozleX7ZsmS5evKgRI0bc0PtYLBZlZGTc0DWKSkpKijw9Pe0dhtXVeUmPHj1UsWJFLVu2TH369LEevzLXc1RX57xdu3bVY489poiICE2dOlXPP/+8Zs6caccIi4bZbObrjVKBkVKAHWX/IDhy5EiebaKionT//ferVq1a8vDwUK1atTRo0KAc52QPz/7xxx/173//WwEBAapUqZLuvvtunThxwtquVq1a2rt3rzZt2mQdnny90698fX0lSa6urtZjuU3fy8zM1KuvvqoGDRrIbDarSpUqGjp0qI4dO2bTrnPnzmrSpIkiIyPVoUMHa38XL14sSfrmm2908803y9PTU02bNtXatWttzv/77781fPhw1atXT56enqpWrZr69Omj33//PUc806dPV/369eXh4SE/Pz81a9ZMb731lrXN6dOn9cgjjygkJERms1mVK1dWx44db2iUU40aNfTGG28oKSlJ7777riTpo48+kslkynV02bRp0+Tq6mr9+mXfn+3bt6tTp07y9PRU7dq19corrygzM9N63qVLl/Tkk0+qRYsW8vX1lb+/v9q3b68vv/wyx3uYTCaNHTtWixcvtt6P1q1ba+vWrTIMQ6+99ppCQ0NVoUIF3X777fr7779tzs9t+l5mZqbeeecdtWjRwnp/27Vrp9WrV1vbbNiwQZ07d1alSpXk4eGhGjVqaMCAAUpJSbnu+wsAKJ2cnZ310EMPaceOHTl+JktZo3iqVq2qsLAwSVJcXJxGjx6t6tWry83NTaGhoZo6dapNwSl76tGrr76q6dOnKzQ0VGazWT/++KP+9a9/SZKGDx9uzXWunJIWFRWlvn37yt/fX+7u7mrZsqX+7//+L9fYz58/r+HDh8vf319eXl7q06ePDh06ZNMm++fzTz/9pA4dOsjT01MPP/ywJGn58uXq1q2bqlatKg8PDzVs2FDPPPOMzSj5YcOGae7cuZJkM83u8OHDkiTDMDRv3jzrz9WKFSvqnnvuyRFHYbi7u8vNzc0mh8t+/6un7/3xxx/q16+fKlasKHd3d7Vo0UIffvihTZvsJSg+/fRTTZo0SVWrVlWFChXUp08fnTx5UklJSXrkkUcUEBCggIAADR8+XBcuXLC5xty5c3XrrbeqSpUq8vLyUtOmTfXqq6/mGJm0a9cu9e7dW1WqVJHZbFZwcLB69eplk1d+/vnnatu2rXx9fa35UvbX5HpFRESocePGmjNnji5duiTDMFSvXj117949R9sLFy7I19dXjz32mM39WbZsmZ577jkFBwfLx8dHXbt21YEDB2zOXb9+vfr166fq1atbl8wYPXp0nktm/Pbbb7r33nutOV94eLgyMjJ04MAB9ejRQ97e3qpVq5ZeffVVm/Pzmr73559/atCgQQoMDJTZbFaNGjU0dOhQpaamSsoquD711FMKDQ2Vu7u7/P391bp1ay1btuyG7i8cFyOlADvK/oBfuXLlPNscPnxY9evX1/333y9/f3/FxsZq/vz5+te//qV9+/YpICDApv3IkSPVq1cvffrppzp69KiefvppPfjgg9qwYYMkaeXKlbrnnnvk6+urefPmScr6Tcm1GIZhTQYvXbqk7du3a/bs2erYsaNCQ0PzPfff//633nvvPY0dO1a9e/fW4cOH9cILL2jjxo3auXOnTR/i4uI0fPhwTZw4UdWrV9c777yjhx9+WEePHtV///tfPfvss/L19dW0adPUv39/HTp0SMHBwZKyhutXqlRJr7zyiipXrqxz587pww8/VNu2bbVr1y7Vr19fkvTqq68qIiJCzz//vG699Valp6frzz//tFnrYsiQIdq5c6f+85//6KabblJ8fLx27typs2fPXvNe5adnz55ydnbWTz/9JEkaOHCgJk6cqLlz56p9+/bWdhkZGXr33Xd11113WfuXfX8eeOABPfnkk5oyZYpWrlypyZMnKzg4WEOHDpUkpaam6ty5c3rqqadUrVo1paWl6fvvv9fdd9+txYsXW9tl+/rrr7Vr1y698sorMplMmjRpknr16qWHHnpIhw4d0pw5c5SQkKDw8HANGDBAu3fvznfdsGHDhunjjz/WiBEjNG3aNLm5uWnnzp3W5Prw4cPq1auXOnXqpEWLFsnPz0/Hjx/X2rVrlZaWVqp+swwAKBoPP/ywXnnlFS1atEhvvvmm9fi+ffu0bds2PfPMM3J2dlZcXJzatGkjJycnvfjii6pTp44iIyM1ffp0HT582PqLqmxvv/22brrpJr3++uvy8fFRYGCgFi9erOHDh+v555+3Tk3KHsX9448/qkePHmrbtq0WLFggX19fffbZZxo4cKBSUlI0bNgwm+uPGDFCd955pzWvev7559W5c2f99ttv8vPzs7aLjY3Vgw8+qIkTJ+rll1+Wk1PW7/7/+usv9ezZUxMmTJCXl5f+/PNPzZw5U9u2bbPmZi+88IKSk5P13//+1+aXVFWrVpUkjR49WkuWLNG4ceM0c+ZMnTt3TtOmTVOHDh20Z8+eAo28zh5FZhiGTp48qddee03JyckaPHhwvucdOHBAHTp0UJUqVfT222+rUqVK+vjjjzVs2DCdPHlSEydOtGn/7LPPqkuXLlqyZIkOHz6sp556SoMGDZKLi4uaN2+uZcuWadeuXXr22Wfl7e2tt99+23ruP//8o8GDBys0NFRubm7as2eP/vOf/+jPP//UokWLJEnJycm68847FRoaqrlz5yowMFBxcXH68ccflZSUJClrCuLAgQM1cOBARUREyN3dXUeOHLHe7xvRp08fvfLKK4qKitItt9yixx9/XBMmTNBff/2levXqWdstXbpUiYmJ1qLUlfenY8eO+uCDD5SYmKhJkyapT58+2r9/v5ydna33oX379ho5cqR8fX11+PBhzZo1S7fccot+//33HIXE++67Tw8++KBGjx6t9evXWwt533//vR599FE99dRT1mJh3bp1dffdd+fZvz179uiWW25RQECApk2bpnr16ik2NlarV69WWlqazGazwsPD9dFHH2n69Olq2bKlkpOT9ccff9xwjgwHZgAodosXLzYkGVu3bjXS09ONpKQk4+uvvzYqV65seHt7G3FxcYZhGMaPP/5oSDJ+/PHHPK+VkZFhXLhwwfDy8jLeeuutHO/x6KOP2rR/9dVXDUlGbGys9Vjjxo2N2267rcDx16xZ05CU49GmTRub6xqGYUyZMsW48r+W/fv35xrXr7/+akgynn32Weux2267zZBkREVFWY+dPXvWcHZ2Njw8PIzjx49bj+/evduQZLz99tt5xp2RkWGkpaUZ9erVM5544gnr8d69exstWrTIt88VKlQwJkyYkG+b3GR/HbZv355nm8DAQKNhw4bW51OmTDHc3NyMkydPWo8tX77ckGRs2rTJeiz7/vz6668212vUqJHRvXv3PN8vIyPDSE9PN0aMGGG0bNnS5jVJRlBQkHHhwgXrsVWrVhmSjBYtWhiZmZnW47NnzzYkGb/99pv12EMPPWTUrFnT+vynn34yJBnPPfdcnvH897//NSQZu3fvzrMNAKD8ue2224yAgAAjLS3NeuzJJ580JBkHDx40DMMwRo8ebVSoUME4cuSIzbmvv/66IcnYu3evYRiGER0dbUgy6tSpY3M9wzCM7du3G5KMxYsX54ihQYMGRsuWLY309HSb47179zaqVq1qWCwWwzD+9/P8rrvusmn3yy+/GJKM6dOn2/RLkvHDDz/k2//MzEwjPT3d2LRpkyHJ2LNnj/W1xx57zMjto1lkZKQhyXjjjTdsjh89etTw8PAwJk6cmO97Zvfj6ofZbDbmzZuXo70kY8qUKdbn999/v2E2m42YmBibdmFhYYanp6cRHx9vGMb/ctg+ffrYtJswYYIhyRg3bpzN8f79+xv+/v55xm2xWIz09HRj6dKlhrOzs3Hu3DnDMAwjKirKkGSsWrUqz3Ozv1eyYyuMmjVrGr169crz9fnz5xuSjOXLlxuGYRiJiYmGt7e3MX78eJt2jRo1Mrp06WJ9nn1/evbsadPu//7v/wxJRmRkZK7vl/09c+TIEUOS8eWXX1pfy865r/7eaNGihSHJWLFihfVYenq6UblyZePuu++2Hsv+N3Tlv5Pbb7/d8PPzM06dOpXnPWjSpInRv3//PF8HCsuhp+/99NNP6tOnj4KDg2UymbRq1apCnX/p0iUNGzZMTZs2lYuLi/r3759ru02bNqlVq1Zyd3dX7dq1tWDBghsPHmVSu3bt5OrqKm9vb/Xu3VtBQUH69ttv8/0N14ULF6y/2XBxcZGLi4sqVKig5ORk7d+/P0f7vn372jxv1qyZpPynCBbELbfcou3bt2v79u365ZdftHDhQp0+fVq33357juHEV/rxxx8lKcdvHtu0aaOGDRvqhx9+sDletWpVtWrVyvrc399fVapUUYsWLWxGDDVs2DBHvzIyMvTyyy+rUaNGcnNzk4uLi9zc3PTXX3/Z3Ks2bdpoz549evTRR/Xdd98pMTExR9xt2rTRkiVLNH36dG3durVIF7U0DMPm+b///W9J0vvvv289NmfOHDVt2lS33nqrTdugoCC1adPG5lizZs1yfH0///xzdezYURUqVJCLi4tcXV21cOHCXL9nunTpIi8vL+vz7HsbFhZmMyIqt3t+tW+//VaScvxm8EotWrSQm5ubHnnkEX344Yc3NP0AAFB2jBgxQmfOnLFO587IyNDHH3+sTp06WUeZfP311+rSpYuCg4OVkZFhfWRP7du0aZPNNfv27Ztj5Ehe/v77b/3555/WdTyvvH7Pnj0VGxubYyrV1Wt+dujQQTVr1rTmN9kqVqyo22+/Pcd7Hjp0SIMHD1ZQUJCcnZ3l6uqq2267TZJy/Zl8ta+//lomk0kPPvigTbxBQUFq3rx5gXdsXrp0qTWP+/bbb/XQQw/pscce05w5c/I9b8OGDbrjjjsUEhJic3zYsGFKSUnJsfxA7969bZ5n5w5XL6bdsGFDnTt3zmYK365du9S3b19VqlTJeq+GDh0qi8WigwcPSpLq1q2rihUratKkSVqwYIH27duXI+bs6Zv33Xef/u///s9mp+sbdXUO5+3treHDh2vJkiXWKZkbNmzQvn37NHbs2BznFyRPP3XqlMaMGaOQkBBrDlezZk1JuX/P5HbPTSaT9d+MJLm4uKhu3br55nApKSnatGmT7rvvvnxncbRp00bffvutnnnmGW3cuFEXL17Msy1QEA5dlEpOTlbz5s2v+Z9xXiwWizw8PDRu3Dh17do11zbR0dHq2bOnOnXqZB2qOm7cOH3xxRc3EjrKqOyEYNeuXTpx4oR+++03dezYMd9zBg8erDlz5mjkyJH67rvvtG3bNm3fvl2VK1fO9YdApUqVbJ5nT8270R8Yvr6+at26tVq3bq0OHTro4Ycf1qeffqr9+/frjTfeyPO87KG82UPQrxQcHJxjqK+/v3+Odm5ubjmOu7m5ScoqDmcLDw/XCy+8oP79++urr77Sr7/+qu3bt6t58+Y2/Z88ebJef/11bd26VWFhYapUqZLuuOMOm22fly9froceekgffPCB2rdvL39/fw0dOlRxcXH53aZrSk5O1tmzZ20KbIGBgRo4cKDeffddWSwW/fbbb9q8eXOuyczVX18p62t8Zf9WrFih++67T9WqVdPHH3+syMhIbd++XQ8//LDN/cqW170tyD2/2unTp+Xs7KygoKA829SpU0fff/+9qlSposcee0x16tRRnTp1bNb0AgCUP9nLB2RPwVuzZo1Onjxps8D5yZMn9dVXX8nV1dXm0bhxY0nK8Yuw3PKLvJw8eVKS9NRTT+W4fvYC7FdfP7efZ0FBQTnyl9ziuHDhgjp16qRff/1V06dP18aNG7V9+3atWLFCUsFys5MnT8owDAUGBuaIeevWrfn+YvBKDRs2tOZxPXr00Lvvvqtu3bpp4sSJNssXXO3s2bN55nDZr1/penOKmJgYderUScePH9dbb72lzZs3a/v27da1trLvla+vrzZt2qQWLVro2WefVePGjRUcHKwpU6ZYf4F46623atWqVcrIyNDQoUNVvXp1NWnSpEjWPMou6lyZxz3++ONKSkrSJ598IinrF4vVq1dXv379cpx/rTw9MzNT3bp104oVKzRx4kT98MMP2rZtm7Zu3WrT7kq53VtPT0+5u7vnOJ5fDnf+/HlZLJZrblj09ttva9KkSVq1apW6dOkif39/9e/fX3/99Ve+5wF5ceg1pcLCwmwqyFdLS0vT888/r08++UTx8fFq0qSJZs6cqc6dO0uSvLy8NH/+fEnSL7/8kut/6AsWLFCNGjU0e/ZsSVk/EKKiovT6669rwIABRd0llHLZCUFBJSQk6Ouvv9aUKVP0zDPPWI9nrxlkb9m/3dmzZ0+ebbJ/+MbGxub4IXfixIkca2LdiI8//lhDhw7Vyy+/bHP8zJkzNus+uLi4KDw8XOHh4YqPj9f333+vZ599Vt27d9fRo0fl6empgIAAzZ49W7Nnz1ZMTIxWr16tZ555RqdOncqxwHphfPPNN7JYLNb/R7KNHz9eH330kb788kutXbtWfn5+170j48cff6zQ0FAtX77cZqRT9gKVxaly5cqyWCyKi4vL94NCp06d1KlTJ1ksFkVFRemdd97RhAkTFBgYqPvvv7/Y4wQAlDwPDw8NGjRI77//vmJjY7Vo0SJ5e3vr3nvvtbYJCAhQs2bN9J///CfXa1xZDJCU7xqHV8vOOSZPnpznujrZ609my+2XUXFxcapbt+4149iwYYNOnDihjRs3WkdHScq3CJRbzCaTSZs3b851DdCCrAual2bNmum7777TwYMHc4zCzlapUiXFxsbmOJ69CUtR5XGrVq1ScnKyVqxYYR0VJEm7d+/O0bZp06b67LPPZBiGfvvtNy1ZskTTpk2Th4eHNV/u16+f+vXrp9TUVG3dulUzZszQ4MGDVatWLZs1PAvDMAx99dVX8vLyssnn69atq7CwMM2dO1dhYWFavXq1pk6dal0jqjD++OMP7dmzR0uWLNFDDz1kPX71RjPFwd/fX87Ozjk2Irqal5eXpk6dqqlTp+rkyZPWUVN9+vTRn3/+Wexxovxx6JFS1zJ8+HD98ssv+uyzz6y7GvTo0aNQVeDIyEh169bN5lj37t0VFRVVpNOBUD6ZTCYZhpEj4fjggw9ksViu+7pXj6y5XtmJQpUqVfJskz2U/eOPP7Y5vn37du3fv1933HHHDceRzWQy5bhX33zzTb7Dtv38/HTPPffoscce07lz56yLcV+pRo0aGjt2rO68807t3LnzuuOLiYnRU089JV9fX40ePdrmtVatWqlDhw6aOXOmPvnkEw0bNsxmSl1hmEwmubm52STIcXFxue6+V9SyC/3ZBftrcXZ2Vtu2ba2/Cb2R+wsAKP1GjBghi8Wi1157TWvWrNH9999vs8FF79699ccff6hOnTrWkT1XPq4uSuUmr1Hi9evXV7169bRnz55cr926dWt5e3vbnJM9+iXbli1bdOTIkRy/XMpN9s/hq3OT7B14CxJz7969ZRiGjh8/nmu8TZs2vWYcecnO4/KbqnXHHXdYi2tXWrp0qTw9Pa07Sd+o3O6VYRg2Sxvkdk7z5s315ptvys/PL9ccwmw267bbbtPMmTMlZU0RvF5Tp07Vvn37NH78+ByjkMaPH6/ffvtNDz30kJydnTVq1Kjreo/CfM8UNQ8PD9122236/PPPCzwCLzAwUMOGDdOgQYN04MABdlHGdXHokVL5+eeff7Rs2TIdO3bM+sPvqaee0tq1a7V48eIcIzHyEhcXl2O9oMDAQGVkZOjMmTOFGnIMx+Pj46Nbb71Vr732mgICAlSrVi1t2rRJCxcutBn5U1jZv2Favny5ateuLXd392smNfHx8dahw+np6dq/f79efvllmc3mfNcPql+/vh555BG98847cnJyUlhYmHX3vZCQED3xxBPX3Y+r9e7dW0uWLFGDBg3UrFkz7dixQ6+99lqOEVp9+vRRkyZN1Lp1a1WuXFlHjhzR7NmzVbNmTdWrV08JCQnq0qWLBg8erAYNGsjb21vbt2/X2rVr892x5Ep//PGHdd2HU6dOafPmzVq8eLGcnZ21cuXKXBPA8ePHa+DAgTKZTNZpBNd7H1asWKFHH31U99xzj44ePaqXXnpJVatWLfah1Z06ddKQIUM0ffp0nTx5Ur1795bZbNauXbvk6empxx9/XAsWLNCGDRvUq1cv1ahRQ5cuXbLuqpPXVGgAQPnQunVrNWvWTLNnz5ZhGDZT9yRp2rRpWr9+vTp06KBx48apfv36unTpkg4fPqw1a9ZowYIF15xeVKdOHXl4eOiTTz5Rw4YNVaFCBQUHBys4OFjvvvuuwsLC1L17dw0bNkzVqlXTuXPntH//fu3cuVOff/65zbWioqI0cuRI3XvvvTp69Kiee+45VatWrUA/pzt06KCKFStqzJgxmjJlilxdXfXJJ5/kOsI8Ow+bOXOmwsLC5OzsrGbNmqljx4565JFHNHz4cEVFRenWW2+Vl5eXYmNj9fPPP6tp06bWtSnzk52XSFlT7lasWKH169frrrvuyncX5SlTpljX+XrxxRfl7++vTz75RN98841effVV+fr6XvO9C+LOO++Um5ubBg0apIkTJ+rSpUuaP3++zp8/b9Pu66+/1rx589S/f3/Vrl1bhmFoxYoVio+P15133ilJevHFF3Xs2DHdcccdql69uuLj4/XWW2/ZrOeVnytz3uTkZB04cECfffaZNm/erPvuu09Tp07NNf5GjRrpxx9/1IMPPpjvL2zz06BBA9WpU0fPPPOMDMOQv7+/vvrqK61fv/66rldY2bv8tW3bVs8884zq1q2rkydPavXq1Xr33Xfl7e2ttm3bqnfv3mrWrJkqVqyo/fv366OPPlL79u3ZQRnXhaJUHnbu3CnDMHTTTTfZHE9NTc11TZf8XD2cN3uBvMIMN4bj+vTTTzV+/HhNnDhRGRkZ6tixo9avX59jwcjCmDp1qmJjYzVq1CglJSWpZs2auY4QutIvv/xiHe7s7OysatWqqU2bNnruuefUokWLfM+dP3++6tSpo4ULF2ru3Lny9fVVjx49NGPGjEL/e8pPdsIxY8YMXbhwQTfffLNWrFih559/3qZdly5d9MUXX1i34w0KCtKdd96pF154Qa6urnJ3d1fbtm310Ucf6fDhw0pPT1eNGjU0adKkHFsf52X48OGSsubv+/n5qWHDhpo0aZJGjhyZ528k+/fvL7PZrC5duthsK1xYw4cP16lTp7RgwQItWrRItWvX1jPPPKNjx47lmkgVtSVLlujmm2/WwoULtWTJEnl4eKhRo0Z69tlnJWUtdL5u3TpNmTJFcXFxqlChgpo0aaLVq1fnGFkKACh/RowYofHjx6tRo0Zq27atzWtVq1ZVVFSUXnrpJb322ms6duyYvL29FRoaqh49eqhixYrXvL6np6cWLVqkqVOnqlu3bkpPT9eUKVMUERGhLl26aNu2bfrPf/6jCRMm6Pz586pUqZIaNWqk++67L8e1Fi5cqI8++kj333+/UlNT1aVLF7311lu5roF5tUqVKumbb77Rk08+qQcffFBeXl7q16+fli9frptvvtmm7eDBg/XLL79o3rx5mjZtmgzDUHR0tGrVqqV3331X7dq107vvvqt58+YpMzNTwcHB6tixY57T7q6WnZdIWesyhYaGatasWdcsrtWvX19btmzRs88+q8cee0wXL15Uw4YNtXjx4hyb2NyIBg0a6IsvvtDzzz+vu+++W5UqVdLgwYMVHh5us9xKvXr15Ofnp1dffVUnTpyQm5ub6tevbzPdrW3btoqKitKkSZN0+vRp+fn5qXXr1tqwYYN1bbL8ZOe8JpNJXl5e1pz3+eefzzdPue+++xQREZHrmqAF5erqqq+++krjx4/X6NGj5eLioq5du+r7779XjRo1rvu6BdW8eXNt27ZNU6ZM0eTJk5WUlKSgoCDdfvvt1nXAbr/9dq1evVpvvvmmUlJSVK1aNQ0dOlTPPfdcsceH8slkXL2FgIMymUxauXKldQe95cuX64EHHtDevXtzzAeuUKFCjkUPhw0bpvj4+Bw7+N16661q2bKlzQK+K1eu1H333aeUlJQC7xYCoPz76quv1LdvX33zzTfq2bOnvcMBAABAAbVu3Vomk0nbt2+3dyhAmcJIqTy0bNlSFotFp06dUqdOna77Ou3bt9dXX31lc2zdunVq3bo1BSkAkqR9+/bpyJEjevLJJ9WiRYt8N2AAAABA6ZCYmKg//vhDX3/9tXbs2KGVK1faOySgzHHootSFCxdsdjKIjo7W7t275e/vr5tuukkPPPCAhg4dqjfeeEMtW7bUmTNntGHDBjVt2tQ6imHfvn1KS0vTuXPnlJSUZF0wMHs605gxYzRnzhyFh4dr1KhRioyM1MKFC4tkS1IA5cOjjz6qX375RTfffLM+/PBDpvYCAACUATt37lSXLl1UqVIlTZkyxTrrBkDBOfT0vY0bN6pLly45jj/00ENasmSJ0tPTNX36dC1dulTHjx9XpUqV1L59e02dOtW6GGGtWrV05MiRHNe48rZu2rRJTzzxhPbu3avg4GBNmjRJY8aMKb6OAQAAAAAAlHIOXZQCAAAAAACAfTjZOwAAAAAAAAA4HopSAAAAAAAAKHEOt9B5ZmamTpw4IW9vbxYTBgDAQRmGoaSkJAUHB8vJid/RFSVyLQAAUNBcy+GKUidOnFBISIi9wwAAAKXA0aNHVb16dXuHUa6QawEAgGzXyrUcrijl7e0tKevG+Pj42DkaAABgD4mJiQoJCbHmBSg65FoAAKCguZbDFaWyh5H7+PiQKAEA4OCYXlb0yLUAAEC2a+VaLKIAAAAAAACAEkdRCgAAAAAAACWOohQAAAAAAABKnMOtKVVQFotF6enp9g6jXHN1dZWzs7O9wwAAAAAAXAc+Nzuuovo8T1HqKoZhKC4uTvHx8fYOxSH4+fkpKCiIhWYBAAAAoIzgczOkovk8T1HqKtn/sKpUqSJPT0+KJcXEMAylpKTo1KlTkqSqVavaOSIAAAAAQEHwudmxFeXneYpSV7BYLNZ/WJUqVbJ3OOWeh4eHJOnUqVOqUqUKU/kAAAAAoJTjczOkovs8z0LnV8ieC+vp6WnnSBxH9r1mHjIAAAAAlH58bka2ovg8T1EqFww9LDncawAAAAAoe/gsh6L4HqAoBQAAAAAAgBJHUcoBbNy4USaTiZ0RAAAAAAAoBQ4fPiyTyaTdu3fbOxS7oihVTgwbNkwmk0kmk0murq6qXbu2nnrqKSUnJxfbe5pMJq1atarYrg8AAAAAQFHo06ePunbtmutrkZGRMplM2rlzZ6GvW1qLS7Vq1dLs2bPtHcY1UZQqR3r06KHY2FgdOnRI06dP17x58/TUU0/ZOywAAAAAAOxqxIgR2rBhg44cOZLjtUWLFqlFixa6+eabC3XNtLS0ogqvQAzDUEZGRom+Z3GjKFWOmM1mBQUFKSQkRIMHD9YDDzyQ60ims2fPatCgQapevbo8PT3VtGlTLVu2zKZN586dNW7cOE2cOFH+/v4KCgpSRESE9fVatWpJku666y6ZTCbr8z179qhLly7y9vaWj4+PWrVqpaioqGLqMQAAAAAA19a7d29VqVJFS5YssTmekpKi5cuXa8SIEdqyZYtuvfVWeXh4KCQkROPGjbOZfVSrVi1Nnz5dw4YNk6+vr0aNGqXQ0FBJUsuWLWUymdS5c2dr+8WLF6thw4Zyd3dXgwYNNG/evBxx/fnnn+rQoYPc3d3VuHFjbdy40fpa9lI83333nVq3bi2z2azNmzfrn3/+Ub9+/RQYGKgKFSroX//6l77//nvreZ07d9aRI0f0xBNPWGdUZbtWH+fNm6d69erJ3d1dgYGBuueee673lheIS7FevRwwDEMX0y12eW8PV+cbWs3ew8Mj160ZL126pFatWmnSpEny8fHRN998oyFDhqh27dpq27attd2HH36o8PBw/frrr4qMjNSwYcPUsWNH3Xnnndq+fbuqVKmixYsXq0ePHnJ2dpYkPfDAA2rZsqXmz58vZ2dn7d69W66urtfdB6A8MwxDhiEZkjIv/z3TMKTLzzONrDaZhiQj+5hhbS/j8nWs18t+blz1/H/vd+VxXXVegc7NcU5e7fO5Zh5xXh1PXrFcGc+NuPErXL5OUV2oCCIqqlhK071xd3VSs+p+N34hAABQbOz1ubkwn5ldXFw0dOhQLVmyRC+++KL1vM8//1xpaWlq3ry5unfvrpdeekkLFy7U6dOnNXbsWI0dO1aLFy+2Xue1117TCy+8oOeff16SNHbsWLVp00bff/+9GjduLDc3N0nS+++/rylTpmjOnDlq2bKldu3apVGjRsnLy0sPPfSQ9XpPP/20Zs+erUaNGmnWrFnq27evoqOjValSJWubiRMn6vXXX1ft2rXl5+enY8eOqWfPnpo+fbrc3d314Ycfqk+fPjpw4IBq1KihFStWqHnz5nrkkUc0atQo63V+//33fPsYFRWlcePG6aOPPlKHDh107tw5bd68+fq/QAVgMooisy9DEhMT5evrq4SEBPn4+Ni8dunSJUVHRys0NFTu7u6SpJS0DDV68Tt7hKp907rL061gdcNhw4YpPj7eOjJq27Zt6tmzp+644w79+9//VpcuXXT+/Hn5+fnlen6vXr3UsGFDvf7665KyKqsWi8XmG7BNmza6/fbb9corr0jKWlNq5cqV6t+/v7WNj4+P3nnnHZt/ZPnJ7Z7DMaRbMnUx3aJL6RalWwylZ2Qq3ZKp1Mt/plsMpVsylWbJvPzaFc+vOJZmyVTa5XMyMg1lWLIKNxmZmbJkSpYr/zSynxt5HLviYfzvWlces1j+VxTKKhpdLi7JtrCUfTzTyCquZBeYrK+rKIsZQPlVu7KXNjzZucivm18+gBvDvQWA8q00fW4uzGdmKWtUUsOGDbVhwwZ16dJFknTbbbepWrVqcnFxkYeHh959911r+59//lm33XabkpOT5e7urlq1aqlly5ZauXKltc3hw4cVGhqqXbt2qUWLFtbjNWrU0MyZMzVo0CDrsenTp2vNmjXasmWL9bxXXnlFkyZNkiRlZGQoNDRUjz/+uCZOnKiNGzeqS5cuWrVqlfr165dv3xo3bqx///vfGjt2rKSsUV0TJkzQhAkTrG2GDh2abx/XrFmj4cOH69ixY/L29r7m/czv83xB8wFGSpUjX3/9tSpUqKCMjAylp6erX79+euedd7Rv3z6bdhaLRa+88oqWL1+u48ePKzU1VampqfLy8rJp16xZM5vnVatW1alTp/KNITw8XCNHjtRHH32krl276t5771WdOnWKpoMocYZhKCXNogupGUq6lKELqRm6cClDF1LTlXTpf8cuplt0Mc2i1IysPy+mW3QxPVOXLhedLqZZdCnDootpWccupltkyaQicyNMJskkWX/DY7rieNZz619sXs+tjemqNldfM7frWNvkev3cY1Ie71eYmJTLOTeiKK4hXXG/b/Q6RdGnG79E1nWK6Obc6FWqV/QskjgAAAAaNGigDh06aNGiRerSpYv++ecfbd68WevWrdP48eP1999/65NPPrG2NwxDmZmZio6OVsOGDSVJrVu3vub7nD59WkePHtWIESNsRiplZGTI19fXpm379u2tf3dxcVHr1q21f/9+mzZXv2dycrKmTp2qr7/+WidOnFBGRoYuXryomJiYfOPasWNHvn288847VbNmTdWuXVs9evRQjx49dNddd8nTs/jyMYpS1+Dh6qx907rb7b0Lo0uXLpo/f75cXV0VHBxsnTZ3dVHqjTfe0JtvvqnZs2eradOm8vLy0oQJE3Is0nb1tDuTyaTMzMx8Y4iIiNDgwYP1zTff6Ntvv9WUKVP02Wef6a677ipUX1C0MiyZOpucpvMpaTqfnK74lDSdT0nX+ZQ0JVxM1/nkrOfxKWlKvJSuC5cylJSaVXAqidE8bi5OcnN2kquzSa7OTnJ1dpLZJetPV5f/Hctu4+Zy5XPbNs5OJrk4meRkuvzn5efO2cecc3/N2ckkZ9P//m597YpjzldcR5KcTCY5OWUVJJxMlwtFpqzXTZdfN10+nv337NdMVzx3MmVdw+Qk63lXnpt9/f8dK6qyAwAAAHDj7PW5ubCfmaWsBc/Hjh2ruXPnavHixapZs6buuOMOZWZmavTo0Ro3blyOc2rUqGH9+9WDOXKT/bn5/ffft1kiR5J16Zv8XJ3vX/2eTz/9tL777ju9/vrrqlu3rjw8PHTPPfdcc+H1a/XRzc1NO3fu1MaNG7Vu3Tq9+OKLioiI0Pbt2/OcdXWjKEpdg8lkKtRwQHvy8vJS3bp1r9lu8+bN6tevnx588EFJWd+Yf/31l7XyW1Curq6yWHLOG77pppt000036YknntCgQYO0ePFiilLFJDXDotj4S4pNuKTTF1J1KjHrz9NJto9zKWk3VFxydjKpgtlFFcwu8nbP+rNC9p9mF3m6ucjd1Ukers7ycHOW2dU56++uztbj7m7OcnfJej37uLurs8wuThRZAAAAgDKsLH1uvu+++zR+/Hh9+umn+vDDDzVq1CiZTCbdfPPN2rt3b4E+U18pew2pKz8bBwYGqlq1ajp06JAeeOCBfM/funWrbr31VklZI6l27NhhnYKXl82bN2vYsGHWz9kXLlzQ4cOHc8R19ef1gvTRxcVFXbt2VdeuXTVlyhT5+flpw4YNuvvuu/ON6XqVje8aFKm6devqiy++0JYtW1SxYkXNmjVLcXFxhS5K1apVSz/88IM6duwos9ksd3d3Pf3007rnnnsUGhqqY8eOafv27RowYEAx9aT8S0nL0LHzF3X8/EUdi8/683j8RR07n6Lj5y/q9IXUAhebnEySn6eb/DxdVdHTTRU9XeVn82fWa74erjaFJ2+zq9xdKRwBAAAAKPsqVKiggQMH6tlnn1VCQoKGDRsmSZo0aZLatWunxx57zLog+f79+7V+/Xq98847eV6vSpUq8vDw0Nq1a1W9enW5u7vL19dXERERGjdunHx8fBQWFqbU1FRFRUXp/PnzCg8Pt54/d+5c1atXTw0bNtSbb76p8+fP6+GHH863D3Xr1tWKFSvUp08fmUwmvfDCCzlmNdWqVUs//fST7r//fpnNZgUEBFyzj19//bUOHTqkW2+9VRUrVtSaNWuUmZmp+vXrX/8NvwaKUg7ohRdeUHR0tLp37y5PT0898sgj6t+/vxISEgp1nTfeeEPh4eF6//33Va1aNR08eFBnz57V0KFDdfLkSQUEBOjuu+/W1KlTi6kn5YMl01DMuRRFn7mgQ6eTdehMsqJPJyv6TLLiEi9d83x3VycF+3moirdZlb3dVbmCWZW9r3hcfu7v5SZnJwpLAAAAABzbiBEjtHDhQnXr1s06Na9Zs2batGmTnnvuOXXq1EmGYahOnToaOHBgvtdycXHR22+/rWnTpunFF19Up06dtHHjRo0cOVKenp567bXXNHHiRHl5ealp06Y2C49L0iuvvKKZM2dq165dqlOnjr788ksFBATk+55vvvmmHn74YXXo0MFabEpMTLRpM23aNI0ePVp16tRRamqqDMO4Zh/9/Py0YsUKRURE6NKlS6pXr56WLVumxo0bF/IOFxy7712BneBKnqPd87MXUvVnXFLWIzZRB04m6eDJJF1Kz3utLh93F1Wr6KnqFT1Uzc/D+me1y3/6e7kxigkACokd4ooP9xYAyjdH+wyHvLH7HlCKJV5K1x/HErT7WLz2HI3Xb8cSFJuQ+8gns4uTQgO8VLuyl2oHVLD5u6+na67nAAAAAABQllGUAorI0XMp2nrorH6NPqddMef1z+nkHG1MJqmmv6fqB3mrQZCPGgR5q36Qt2pW8mJqHQAAAADAoVCUAq5TdhFq66Fz2nrorI7HX8zRpnpFDzUP8VOL6n5qHuKnxsE+8jLzzw4AAAAAAD4dAwWUbsnUjiPnteHPU9rw5yn9feqCzesuTiY1q+6rdrUrqXWtimpe3U+VKpjtFC0AAAAAAKUbRSkgHwkp6dpw4KR+2H9KPx08rcRLGdbXXJxMah7ip7ah/mpXu5Ja1azIKCgAAAAAAAqIT9C5yMzMeyc0FK3SeK8vpln0w58ntXr3CW08cFpplv/F6O/lps71K+uOBoHqdFOAfNxZhBwAAACA4ymNn+VQsorie4Ci1BXc3Nzk5OSkEydOqHLlynJzc5PJxOLTxcEwDKWlpen06dNycnKSm5ubXePJsGTq57/PaPXuE/pub5yS0yzW1+oHeuvORoG6vWEVNa/ux4LkAAAAABwWn5tRlJ/nKUpdwcnJSaGhoYqNjdWJEyfsHY5D8PT0VI0aNeTk5GSX949LuKRPt8Vo2bYYnU5KtR6vXtFD/VoEq2/zaqof5G2X2AAAAACgtOFzM7IVxed5ilJXcXNzU40aNZSRkSGLxXLtE3DdnJ2d5eLiUuJVdcMw9Gv0OS2NPKzv9p6UJdOQlDU1r3ezqurXIlg316hItR8AAAAAcsHnZhTV53mKUrkwmUxydXWVqyvrBZUnl9It+u+OY1oaeVgHT/5v57w2tfw1pH1NdW8cJDcX+4zYAgAAAICyhM/NKAoUpVDuJV1K10dbj2jRz9E6cyFNkuTh6qz+LatpaPuaaljVx84RAgAAAADgeChKody6kJqhxT9H6/3Nh5R4KUOSVM3PQyNuCdWAVtXl60FFHwAAAAAAe6EohXLnUrpFn/wao3k//q2zyVkjo+pU9tKjneuqb4tguTozRQ8AAAAAAHujKIVywzAMfftHnF5es1/Hzl+UJNWq5Kkn7rxJfZoFy8mJhcsBAAAAACgtKEqhXNh7IkFTv9qnbdHnJEmBPmY90fUmDWhVnZFRAAAAAACUQhSlUKYlp2Zo9vcHteiXw7JkGnJ3ddLoW+to9G215enGtzcAAAAAAKUVn9pRZm3486ReWLVXx+Ozpur1alpVz/VqqGA/DztHBgAAAAAAroWiFMqcpEvpeunrffq/qGOSpOoVPfRS/ybqUr+KnSMDAAAAAAAFRVEKZUrkP2f11Od7dDz+okwmaUTHUIV3u4mpegAAAAAAlDF2XQH6p59+Up8+fRQcHCyTyaRVq1bl2/7nn39Wx44dValSJXl4eKhBgwZ68803SyZY2JUl09Cb6w9q8AdbdTz+oqpX9NBno9rp+d6NKEgBAAAAAFAG2fXTfHJyspo3b67hw4drwIAB12zv5eWlsWPHqlmzZvLy8tLPP/+s0aNHy8vLS4888kgJRAx7OJV0SeOX7VbkobOSpPtaV9eLfRqrgpliFAAAAAAAZZVdP9WHhYUpLCyswO1btmypli1bWp/XqlVLK1as0ObNmylKlVM7jpzT6I926syFVHm6Oes/dzXRXS2r2zssAAAAAABwg8r0UJNdu3Zpy5Ytmj59ep5tUlNTlZqaan2emJhYEqGhCPx3xzE9u+J3pVkydVNgBc17oJXqVqlg77AAAAAAAEARsOuaUterevXqMpvNat26tR577DGNHDkyz7YzZsyQr6+v9RESElKCkeJ6WDINvbxmv576fI/SLJnq3jhQKx/tSEEKAAAAAIBypEwWpTZv3qyoqCgtWLBAs2fP1rJly/JsO3nyZCUkJFgfR48eLcFIUViX0i0a++lOvffTIUnSuDvqaf4DreTF+lEAAAAAAJQrZfKTfmhoqCSpadOmOnnypCIiIjRo0KBc25rNZpnN5pIMD9cp8VK6Hlkapa2HzsnN2Umv39dcfZsH2zssAAAAAABQDMrkSKkrGYZhs2YUyqbTSaka+O5WbT10ThXMLloy/F8UpAAAuIZ58+YpNDRU7u7uatWqlTZv3pxv+02bNqlVq1Zyd3dX7dq1tWDBgjzbfvbZZzKZTOrfv38RRw0AAJDFriOlLly4oL///tv6PDo6Wrt375a/v79q1KihyZMn6/jx41q6dKkkae7cuapRo4YaNGggSfr555/1+uuv6/HHH7dL/Cgap5IuadB7W/XP6WQFVDBryfB/qUk1X3uHBQBAqbZ8+XJNmDBB8+bNU8eOHfXuu+8qLCxM+/btU40aNXK0j46OVs+ePTVq1Ch9/PHH+uWXX/Too4+qcuXKGjBggE3bI0eO6KmnnlKnTp1KqjsAAMAB2bUoFRUVpS5dulifh4eHS5IeeughLVmyRLGxsYqJibG+npmZqcmTJys6OlouLi6qU6eOXnnlFY0ePbrEY0fRuLIgFezrrk9HtVOtAC97hwUAQKk3a9YsjRgxwrrhy+zZs/Xdd99p/vz5mjFjRo72CxYsUI0aNTR79mxJUsOGDRUVFaXXX3/dpihlsVj0wAMPaOrUqdq8ebPi4+NLojsAAMAB2bUo1blzZxmGkefrS5YssXn++OOPMyqqHLm6ILXskXaqWYmCFAAA15KWlqYdO3bomWeesTnerVs3bdmyJddzIiMj1a1bN5tj3bt318KFC5Weni5XV1dJ0rRp01S5cmWNGDHimtMBAQAAbkSZXOgcZV/CxXQNXbiNghQAANfhzJkzslgsCgwMtDkeGBiouLi4XM+Ji4vLtX1GRobOnDmjqlWr6pdfftHChQu1e/fuAseSmppqs75nYmJiwTsCAAAcWplf6Bxlz6V0i0Z9GKU/45JU2dtMQQoAgOtkMplsnhuGkePYtdpnH09KStKDDz6o999/XwEBAQWOYcaMGfL19bU+QkJCCtEDAADgyBgphRKVYcnU48t2advhc/J2d9HSh9tQkAIAoJACAgLk7OycY1TUqVOncoyGyhYUFJRrexcXF1WqVEl79+7V4cOH1adPH+vrmZmZkiQXFxcdOHBAderUyXHdyZMnW9cFlbJGSlGYAgAABUFRCiXGMAy9uHqv1u87KTcXJ30wtLUaVvWxd1gAAJQ5bm5uatWqldavX6+77rrLenz9+vXq169frue0b99eX331lc2xdevWqXXr1nJ1dVWDBg30+++/27z+/PPPKykpSW+99VaehSaz2Syz2XyDPQIAAI6IohRKzIdbDuvTX2NkMklv399SbWtXsndIAACUWeHh4RoyZIhat26t9u3b67333lNMTIzGjBkjKWsE0/Hjx7V06VJJ0pgxYzRnzhyFh4dr1KhRioyM1MKFC7Vs2TJJkru7u5o0aWLzHn5+fpKU4zgAAEBRoCiFEvHzX2f00jf7JUmTwxqoR5MgO0cEAEDZNnDgQJ09e1bTpk1TbGysmjRpojVr1qhmzZqSpNjYWMXExFjbh4aGas2aNXriiSc0d+5cBQcH6+2339aAAQPs1QUAAODgTEb2CpcOIjExUb6+vkpISJCPD1PHSkL0mWT1n/uLEi6m6+6bq+mNe5vnuwgrAADFjXyg+HBvAQBAQfMBdt9DsUpJy9AjS6OUcDFdLUL89PJdTSlIAQAAAAAAilIoXi9+uVd/nbqgKt5mvTekldxdne0dEgAAAAAAKAUoSqHYfLHjmP6745icTNLbg1qqio+7vUMCAAAAAAClBEUpFIu/T13QC1/+IUkaf8dNasdOewAAAAAA4AoUpVDkUjMsGvvpTqWkWdShTiWNvb2uvUMCAAAAAAClDEUpFLm3vv9Lf8YlqZKXm2YPbCFnJxY2BwAAAAAAtihKoUjtPhqvBZv+kST9566mrCMFAAAAAAByRVEKReZSukVP/t9uZRpS/xbB6tEkyN4hAQAAAACAUoqiFIrMm+sP6p/TyarsbVZE38b2DgcAAAAAAJRiFKVQJPYcjdf7mw9Jkl6+q6n8PN3sHBEAAAAAACjNKErhhlkyDT2/6g9lGlK/FsG6s1GgvUMCAAAAAAClHEUp3LBl22L0+/EEeZtd9FyvhvYOBwAAAAAAlAEUpXBDzl5I1WvfHZAkPdntJlXxZrc9AAAAAABwbRSlcENe+fZPJVxMV6OqPnqwXU17hwMAAAAAAMoIilK4bjuOnNfnO45Jkl7q30Quznw7AQAAAACAgqGKgOtiGIZeXrNfknRPq+pqVbOinSMCAAAAAABlCUUpXJd1+05qx5Hzcnd10tPd69s7HAAAAAAAUMZQlEKhZVgy9eraPyVJI24JVaAPi5sDAAAAAIDCoSiFQvu/qGP653SyKnq6avRtdewdDgAAAAAAKIMoSqFQUtIy9Ob3ByVJj99eTz7urnaOCAAAAAAAlEUUpVAoi385rNNJqarh76kH29W0dzgAAAAAAKCMoiiFAktOzdAHmw9JksLvvEluLnz7AAAAAACA60NVAQX28dYjOp+SrtAAL/VpHmzvcAAAAAAAQBlGUQoFcjHNovcvj5J6tHMdOTuZ7BwRAAAAAAAoyyhKoUCWbYvRmQtpql7RQ/1bVrN3OAAAAAAAoIyjKIVrupRu0bs//SNJerRzXbk6820DAAAAAABuDNUFXNPnO47pZGKqqvq6a0ArRkkBAAAAAIAbR1EK+bJkGtYd90bfWltmF2c7RwQAAAAAAMoDilLI14Y/T+nI2RT5erjqvn+F2DscAAAAAABQTlCUQr4W/RwtSRrUpoY83VzsHA0AAAAAACgvKEohT/tOJCry0Fk5O5k0tH1Ne4cDAAAAAADKEYpSyNPiX7JGSYU1CVKwn4edowEAAAAAAOUJRSnk6syFVH25+4Qk6eFbQu0cDQAAAAAAKG8oSiFXn2yNUZolUy1C/HRzjYr2DgcAAAAAAJQzFKWQQ4YlU5/8ekSSNLxjLfsGAwAAAAAAyiWKUsjhxwOndSopVZW83BTWpKq9wwEAAAAAAOUQRSnk8Nm2GEnSgFbV5ebCtwgAAAAAACh6VBxgIy7hkn48cEqSNPBfIXaOBgAAAAAAlFcUpWDj86ijyjSkNrX8VadyBXuHAwAAAAAAyimKUrDKzDS0POqoJOn+NoySAgAAAAAAxYeiFKx++eeMjp2/KG93FxY4BwAAAAAAxcquRamffvpJffr0UXBwsEwmk1atWpVv+xUrVujOO+9U5cqV5ePjo/bt2+u7774rmWAdwGfbs0ZJ3dWymjzcnO0cDQAAAAAAKM/sWpRKTk5W8+bNNWfOnAK1/+mnn3TnnXdqzZo12rFjh7p06aI+ffpo165dxRxp+ZeQkq71e09KYoFzAAAAAABQ/Fzs+eZhYWEKCwsrcPvZs2fbPH/55Zf15Zdf6quvvlLLli2LODrH8u0fsUqzZKpBkLcaB/vaOxwAAAAAAFDO2bUodaMyMzOVlJQkf3//PNukpqYqNTXV+jwxMbEkQitzvtx9QpLUt0WwnSMBAAAAAACOoEwvdP7GG28oOTlZ9913X55tZsyYIV9fX+sjJISpaVc7mXhJW6PPSpL6NKMoBQAAAAAAil+ZLUotW7ZMERERWr58uapUqZJnu8mTJyshIcH6OHr0aAlGWTZ8teeEDENqXbOiQvw97R0OAAAAAABwAGVy+t7y5cs1YsQIff755+ratWu+bc1ms8xmcwlFVjat3pM1da8fU/cAAAAAAEAJKXMjpZYtW6Zhw4bp008/Va9evewdTpl36PQF/XYsQc5OJvVsWtXe4QAAAAAAAAdh15FSFy5c0N9//219Hh0drd27d8vf3181atTQ5MmTdfz4cS1dulRSVkFq6NCheuutt9SuXTvFxcVJkjw8POTry45x1yN7lNQtdQNUqQIjygAAAAAAQMmw60ipqKgotWzZUi1btpQkhYeHq2XLlnrxxRclSbGxsYqJibG2f/fdd5WRkaHHHntMVatWtT7Gjx9vl/jLOsMwmLoHAAAAAADswq4jpTp37izDMPJ8fcmSJTbPN27cWLwBOZj9sUk6dDpZZhcndWscZO9wAAAAAACAAylza0qh6KzblzX98dabKquCuUyueQ8AAAAAAMooilIO7Lu9JyVJ3RoF2jkSAAAAAADgaChKOaij51K0PzZRTiapa0OKUgAAAAAAoGRRlHJQ3+3NmrrXJtRfFb3c7BwNAAAAAABwNBSlHNS6fdlT91jgHAAAAAAAlDyKUg7o7IVURR0+J0nq1pipewAAAAAAoORRlHJAP+w/pUxDahzso+oVPe0dDgAAAAAAcEAUpRzQun1Z60kxdQ8AAAAAANgLRSkHk5yaoZ/+OiNJ6t6EqXsAAAAAAMA+KEo5mM1/nVFaRqZq+HuqfqC3vcMBAAAAAAAOiqKUg9l08LQk6fYGVWQymewcDQAAAAAAcFQUpRyIYRj66XJR6rabKts5GgAAcKPmzZun0NBQubu7q1WrVtq8eXO+7Tdt2qRWrVrJ3d1dtWvX1oIFC2xef//999WpUydVrFhRFStWVNeuXbVt27bi7AIAAHBgFKUcyD+nk3U8/qLcXJzUtra/vcMBAAA3YPny5ZowYYKee+457dq1S506dVJYWJhiYmJybR8dHa2ePXuqU6dO2rVrl5599lmNGzdOX3zxhbXNxo0bNWjQIP3444+KjIxUjRo11K1bNx0/frykugUAAByIyTAMw95BlKTExET5+voqISFBPj4+9g6nRC38OVovfb1PneoF6KMRbe0dDgAAdlMe8oG2bdvq5ptv1vz5863HGjZsqP79+2vGjBk52k+aNEmrV6/W/v37rcfGjBmjPXv2KDIyMtf3sFgsqlixoubMmaOhQ4cWKK7ycG8BAMCNKWg+wEgpB7KJqXsAAJQLaWlp2rFjh7p162ZzvFu3btqyZUuu50RGRuZo3717d0VFRSk9PT3Xc1JSUpSeni5/f0ZYAwCAoudi7wBQMi6lW/TrobOSKEoBAFDWnTlzRhaLRYGBgTbHAwMDFRcXl+s5cXFxubbPyMjQmTNnVLVq1RznPPPMM6pWrZq6du2aZyypqalKTU21Pk9MTCxMVwAAgANjpJSD2HrorFIzMhXs6666VSrYOxwAAFAErt5J1zCMfHfXza19bscl6dVXX9WyZcu0YsUKubu753nNGTNmyNfX1/oICQkpTBcAAIADoyjlIKxT9+pXzjdZBQAApV9AQICcnZ1zjIo6depUjtFQ2YKCgnJt7+LiokqVKtkcf/311/Xyyy9r3bp1atasWb6xTJ48WQkJCdbH0aNHr6NHAADAEVGUchDZRalb6zF1DwCAss7NzU2tWrXS+vXrbY6vX79eHTp0yPWc9u3b52i/bt06tW7dWq6urtZjr732ml566SWtXbtWrVu3vmYsZrNZPj4+Ng8AAICCoCjlAI6eS9Gh08lydjKpQ90Ae4cDAACKQHh4uD744AMtWrRI+/fv1xNPPKGYmBiNGTNGUtYIpit3zBszZoyOHDmi8PBw7d+/X4sWLdLChQv11FNPWdu8+uqrev7557Vo0SLVqlVLcXFxiouL04ULF0q8fwAAoPxjoXMHsPmvM5KkliF+8vVwvUZrAABQFgwcOFBnz57VtGnTFBsbqyZNmmjNmjWqWbOmJCk2NlYxMTHW9qGhoVqzZo2eeOIJzZ07V8HBwXr77bc1YMAAa5t58+YpLS1N99xzj817TZkyRRERESXSLwAA4DhMRvYKlw4iMTFRvr6+SkhIcJjh5Y8v26Wv9pzQhK71NKHrTfYOBwAAu3PEfKCkcG8BAEBB8wGm75VzhmEo8p+zkqT2tStdozUAAAAAAEDJoChVzv1zOllnLqTK7OKk5iF+9g4HAAAAAABAEkWpci/yUNYoqZtrVJS7q7OdowEAAAAAAMhCUaqc23q5KNW+DlP3AAAAAABA6UFRqhwzDEO/Xi5KtWM9KQAAAAAAUIpQlCrH/j51QWcupMnd1UnNQ3ztHQ4AAAAAAIAVRalyLHs9qVY1K8rswnpSAAAAAACg9KAoVY5lryfVLpSpewAAAAAAoHShKFVOGYahrYfOSWKRcwAAAAAAUPpQlCqnDp68oHPJafJwdVaz6n72DgcAAAAAAMAGRalyKnvqXutaFeXmwpcZAAAAAACULlQryqlth7Om7rUN9bdzJAAAAAAAADlRlCqndh45L0lqVZOiFAAAAAAAKH0oSpVDJ+IvKjbhkpydTGoe4mvvcAAAAAAAAHKgKFUO7YzJGiXVsKq3PN1c7BwNAAAAAABAThSlyqEd2VP3alS0cyQAAAAAAAC5oyhVDmWvJ3VzTYpSAAAAAACgdKIoVc5cTLNo74lESdLNjJQCAAAAAAClFEWpcua3Y/HKyDRUxdus6hU97B0OAAAAAABArihKlTM7Y+IlSa1qVpTJZLJvMAAAAAAAAHmgKFXOWBc5Zz0pAAAAAABQihW6KLV27Vr9/PPP1udz585VixYtNHjwYJ0/f75Ig0PhGIahnTFZX4OWrCcFAECxIy8CAAC4foUuSj399NNKTMxaSPv333/Xk08+qZ49e+rQoUMKDw8v8gBRcIfPpuhccprcnJ3UpJqPvcMBAKDcIy8CAAC4fi6FPSE6OlqNGjWSJH3xxRfq3bu3Xn75Ze3cuVM9e/Ys8gBRcDsvT91rWt1XZhdnO0cDAED5R14EAABw/Qo9UsrNzU0pKSmSpO+//17dunWTJPn7+1t/Uwj72BHDelIAAJQk8iIAAIDrV+iRUrfccovCw8PVsWNHbdu2TcuXL5ckHTx4UNWrVy/yAFFwuy/vvNcyxM+ucQAA4CjIiwAAAK5foUdKzZkzRy4uLvrvf/+r+fPnq1q1apKkb7/9Vj169CjUtX766Sf16dNHwcHBMplMWrVqVb7tY2NjNXjwYNWvX19OTk6aMGFCYcMvty6lW3TwZJIkqRlFKQAASkRR5kUAAACOptAjpWrUqKGvv/46x/E333yz0G+enJys5s2ba/jw4RowYMA126empqpy5cp67rnnruv9yrP9sYnKyDRUyctNwb7u9g4HAACHUJR5EQAAgKMpdFFq586dcnV1VdOmTSVJX375pRYvXqxGjRopIiJCbm5uBb5WWFiYwsLCCty+Vq1aeuuttyRJixYtKlzg5dzvxxMkZS1ybjKZ7BwNAACOoSjzIgAAAEdT6Ol7o0eP1sGDByVJhw4d0v333y9PT099/vnnmjhxYpEHeKNSU1OVmJho8yiPfjuWVZRqVs3XzpEAAOA4ylpeBAAAUJoUuih18OBBtWjRQpL0+eef69Zbb9Wnn36qJUuW6Isvvijq+G7YjBkz5Ovra32EhITYO6Ri8fux7JFSfvYNBAAAB1LW8iIAAIDSpNBFKcMwlJmZKSlr6+OePXtKkkJCQnTmzJmija4ITJ48WQkJCdbH0aNH7R1SkUtJy9Bfpy4vcl6dkVIAAJSUspYXAQAAlCaFXlOqdevWmj59urp27apNmzZp/vz5kqTo6GgFBgYWeYA3ymw2y2w22zuMYrXvRKIyDamKt1mBPixyDgBASSlreREAAEBpUuiRUrNnz9bOnTs1duxYPffcc6pbt64k6b///a86dOhQ5AHi2qzrSTFKCgCAEkVeBAAAcP0KPVKqWbNm+v3333Mcf+211+Ts7Fyoa124cEF///239Xl0dLR2794tf39/1ahRQ5MnT9bx48e1dOlSa5vdu3dbzz19+rR2794tNzc3NWrUqLBdKTesO+9V87NvIAAAOJiizIsAAAAcTaGLUnlxdy/8tLGoqCh16dLF+jw8PFyS9NBDD2nJkiWKjY1VTEyMzTktW7a0/n3Hjh369NNPVbNmTR0+fPj6Ai8HfjsWL4mRUgAAlBbXkxcBAAA4mkIXpSwWi95880393//9n2JiYpSWlmbz+rlz5wp8rc6dO8swjDxfX7JkSY5j+bV3REmX0nXoTLIkqUk1ilIAAJSkosyLAAAAHE2h15SaOnWqZs2apfvuu08JCQkKDw/X3XffLScnJ0VERBRDiMjP3hOJMgwp2Nddlb3L94LuAACUNuRFAAAA16/QRalPPvlE77//vp566im5uLho0KBB+uCDD/Tiiy9q69atxREj8vH75UXOmzJ1DwCAEkdeBAAAcP0KXZSKi4tT06ZNJUkVKlRQQkJWUaR379765ptvijY6XNNvx7N33vOzbyAAADgg8iIAAIDrV+iiVPXq1RUbGytJqlu3rtatWydJ2r59u8xmpo+VtD8uF6VYTwoAgJJHXgQAAHD9Cl2Uuuuuu/TDDz9IksaPH68XXnhB9erV09ChQ/Xwww8XeYDIW3Jqhg6fzVrkvHGwj52jAQDA8ZAXAQAAXL9C7773yiuvWP9+zz33qHr16tqyZYvq1q2rvn37FmlwyN+fcUkyDKmyt1kBFfhtLAAAJY28CAAA4PoVuih1tXbt2qldu3ZFEQsKaX9soiSpYVVGSQEAUBqQFwEAABTcdRWljh8/rl9++UWnTp1SZmamzWvjxo0rksBwbf8rSnnbORIAABwXeREAAMD1KXRRavHixRozZozc3NxUqVIlmUwm62smk4nkqwRlF6UaMVIKAAC7IC8CAAC4foUuSr344ot68cUXNXnyZDk5FXqddBSRzExDB+KSJDF9DwAAeyEvAgAAuH6Fzp5SUlJ0//33k3jZ2dHzKUpOs8jNxUm1A7zsHQ4AAA6JvAgAAOD6FTqDGjFihD7//PPiiAWFkD1176bACnJxJhEGAMAeyIsAAACuX6Gn782YMUO9e/fW2rVr1bRpU7m6utq8PmvWrCILDnnbF3t56l4QU/cAALAX8iIAAIDrV+ii1Msvv6zvvvtO9evXl6QcC3qiZPxv5z2KUgAA2At5EQAAwPUrdFFq1qxZWrRokYYNG1YM4aCgKEoBAGB/5EUAAADXr9CLEZnNZnXs2LE4YkEBJV5K17HzFyVJjShKAQBgN+RFAAAA16/QRanx48frnXfeKY5YUEB/Xl5PKtjXXb6ertdoDQAAigt5EQAAwPUr9PS9bdu2acOGDfr666/VuHHjHAt6rlixosiCQ+6YugcAQOlAXgQAAHD9Cl2U8vPz0913310csaCAKEoBAFA6kBcBAABcv0IXpRYvXlwccaAQKEoBAFA6kBcBAABcv0IXpbKdPn1aBw4ckMlk0k033aTKlSsXZVzIgyXT0IGTWWtKNazqbedoAACARF4EAABwPQq90HlycrIefvhhVa1aVbfeeqs6deqk4OBgjRgxQikpKcURI65w9FyKLqVnyuzipJqVvOwdDgAADo28CAAA4PoVuigVHh6uTZs26auvvlJ8fLzi4+P15ZdfatOmTXryySeLI0Zc4eDlUVJ1KleQs5PJztEAAODYyIsAAACuX6Gn733xxRf673//q86dO1uP9ezZUx4eHrrvvvs0f/78oowPV/nr1AVJ0k2BFewcCQAAIC8CAAC4foUeKZWSkqLAwMAcx6tUqcIw9RLw1+WRUvUCWU8KAAB7Iy8CAAC4foUuSrVv315TpkzRpUuXrMcuXryoqVOnqn379kUaHHI6eDJrpFS9KoyUAgDA3siLAAAArl+hi1JvvfWWtmzZourVq+uOO+5Q165dFRISoi1btuitt94qjhhxmSXT0D+ns6fvMVIKAAB7s3deNG/ePIWGhsrd3V2tWrXS5s2b822/adMmtWrVSu7u7qpdu7YWLFiQo80XX3yhRo0ayWw2q1GjRlq5cmVxhQ8AABxcoYtSTZo00V9//aUZM2aoRYsWatasmV555RX99ddfaty4cXHEiMuOnktRakbWznsh/p72DgcAAIdnz7xo+fLlmjBhgp577jnt2rVLnTp1UlhYmGJiYnJtHx0drZ49e6pTp07atWuXnn32WY0bN05ffPGFtU1kZKQGDhyoIUOGaM+ePRoyZIjuu+8+/frrr8XaFwAA4JhMhmEY9g6iJCUmJsrX11cJCQny8fGxdziFsm5vnB75aIcaVfXRmvGd7B0OAABlVlnOB7K1bdtWN998s81i6g0bNlT//v01Y8aMHO0nTZqk1atXa//+/dZjY8aM0Z49exQZGSlJGjhwoBITE/Xtt99a2/To0UMVK1bUsmXLChRXebi3AADgxhQ0HyjQ7nurV69WWFiYXF1dtXr16nzb9u3bt3CRosDYeQ8AAPsrDXlRWlqaduzYoWeeecbmeLdu3bRly5Zcz4mMjFS3bt1sjnXv3l0LFy5Uenq6XF1dFRkZqSeeeCJHm9mzZxdp/NfLMAxdTLfYOwwAAMoVD1dnmUwmu7x3gYpS/fv3V1xcnKpUqaL+/fvn2c5kMsliIVEoLuy8BwCA/ZWGvOjMmTOyWCw5dv4LDAxUXFxcrufExcXl2j4jI0NnzpxR1apV82yT1zUlKTU1VampqdbniYmJhe1OgV1Mt6jRi98V2/UBAHBE+6Z1l6dbgcpDRa5A75qZmZnr31Gy2HkPAAD7K0150dW/1TQMI9/fdObW/urjhb3mjBkzNHXq1ALHDAAAkM0+pTAUGjvvAQCAbAEBAXJ2ds4xgunUqVM5RjplCwoKyrW9i4uLKlWqlG+bvK4pSZMnT1Z4eLj1eWJiokJCQgrVn4LycHXWvmndi+XaAAA4Kg9XZ7u9d6GKUklJSTp48KDq16+vChUqaOfOnZo9e7YuXryo/v3764EHHiiuOB0eO+8BAFC62DMvcnNzU6tWrbR+/Xrddddd1uPr169Xv379cj2nffv2+uqrr2yOrVu3Tq1bt5arq6u1zfr1623WlVq3bp06dOiQZyxms1lms/lGulNgJpPJbtMLAABA0SvwT/WffvpJvXv31oULF6w7sNxzzz2qVq2anJ2dtWLFCqWkpGjUqFHFGa/DOnh5Pak6lSvI2ck+C5ABAIAspSEvCg8P15AhQ9S6dWu1b99e7733nmJiYjRmzBhJWSOYjh8/rqVLl0rK2mlvzpw5Cg8P16hRoxQZGamFCxfa7Ko3fvx43XrrrZo5c6b69eunL7/8Ut9//71+/vnnYusHAABwXE4Fbfj888/r3nvvVUxMjCZMmKCBAwdq7Nix2r9/v/744w9NnTpVc+fOLc5YHRo77wEAUHqUhrxo4MCBmj17tqZNm6YWLVrop59+0po1a1SzZk1JUmxsrGJiYqztQ0NDtWbNGm3cuFEtWrTQSy+9pLffflsDBgywtunQoYM+++wzLV68WM2aNdOSJUu0fPlytW3btlj7AgAAHJPJyF7h8hr8/Py0detWNWjQQGlpafLw8NDOnTvVvHlzSdLff/+tli1bKikpqVgDvlGJiYny9fVVQkKCfHx87B1OgU34bJdW7T6hp7vX12Nd6to7HAAAyrQbzQfKS15UHMpqrgUAAIpOQfOBAo+USkxMlL+/v6SsdQw8PT3l7f2/Bbe9vb2VkpJyAyEjP+y8BwBA6UFeBAAAcOMKXJQymUw5tgvOb3tgFB123gMAoHQhLwIAALhxBV7o3DAM3XHHHXJxyTolJSVFffr0kZubmyQpIyOjeCKEjp3P2nnPjZ33AAAoFciLAAAAblyBi1JTpkyxeZ7bdsNXLpSJopM9Sqp2gBc77wEAUAqQFwEAANy46y5KoeQcOp0sSapTmfWkAAAoDciLAAAAblyB15SC/fxzuShVu7KXnSMBAAAAAAAoGhSlyoBD2dP3KEoBAAAAAIBygqJUGXDozOWRUgFM3wMAAAAAAOUDRalSLulSuk4npUpipBQAAAAAACg/KEqVctmLnFf2Nsvb3dXO0QAAAAAAABSNAu2+9/bbbxf4guPGjbvuYJDToTOX15MKYJQUAAClAXkRAABA0ShQUerNN98s0MVMJlOhkq+ffvpJr732mnbs2KHY2FitXLlS/fv3z/ecTZs2KTw8XHv37lVwcLAmTpyoMWPGFPg9y5pD1p33WE8KAIDSoLjyIgAAAEdToKJUdHR0sbx5cnKymjdvruHDh2vAgAEFiqNnz54aNWqUPv74Y/3yyy969NFHVbly5QKdXxZlF6XqsJ4UAAClQnHlRQAAAI6mQEWp4hIWFqawsLACt1+wYIFq1Kih2bNnS5IaNmyoqKgovf766+W2KPXP6cvT9yhKAQAAAACAcuS6ilLHjh3T6tWrFRMTo7S0NJvXZs2aVSSB5SYyMlLdunWzOda9e3ctXLhQ6enpcnXNuRB4amqqUlNTrc8TExOLLb6ilplp6PDZy9P3Api+BwBAaWSvvAgAAKCsK3RR6ocfflDfvn0VGhqqAwcOqEmTJjp8+LAMw9DNN99cHDFaxcXFKTAw0OZYYGCgMjIydObMGVWtWjXHOTNmzNDUqVOLNa7iciLhoi6lZ8rV2aTqFT3sHQ4AALiKPfMiAACAss6psCdMnjxZTz75pP744w+5u7vriy++0NGjR3Xbbbfp3nvvLY4YbZhMJpvnhmHkevzKeBMSEqyPo0ePFnuMRSV7Pamalbzk4lzoLxUAAChm9s6LAAAAyrJCVzr279+vhx56SJLk4uKiixcvqkKFCpo2bZpmzpxZ5AFeKSgoSHFxcTbHTp06JRcXF1WqVCnXc8xms3x8fGweZcWh7PWkAlhPCgCA0sieeREAAEBZV+iilJeXl3WNpuDgYP3zzz/W186cOVN0keWiffv2Wr9+vc2xdevWqXXr1rmuJ1XWHTpzeT2pyqwnBQBAaWTPvAgAAKCsK/SaUu3atdMvv/yiRo0aqVevXnryySf1+++/a8WKFWrXrl2hrnXhwgX9/fff1ufR0dHavXu3/P39VaNGDU2ePFnHjx/X0qVLJUljxozRnDlzFB4erlGjRikyMlILFy7UsmXLCtuNMiF7+h477wEAUDoVZV4EAADgaApdlJo1a5YuXMiaVhYREaELFy5o+fLlqlu3rt58881CXSsqKkpdunSxPg8PD5ckPfTQQ1qyZIliY2MVExNjfT00NFRr1qzRE088oblz5yo4OFhvv/22BgwYUNhulAlM3wMAoHQryrwIAADA0ZiM7JXCHURiYqJ8fX2VkJBQqteXSknLUKMXv5Mk7XzhTvl7udk5IgAAyo+ykg+URdxbAABQ0Hyg0GtK1a5dW2fPns1xPD4+XrVr1y7s5ZCH6MvrSfl5ulKQAgCglCIvAgAAuH6FLkodPnxYFoslx/HU1FQdP368SILC/4pSTN0DAKD0Ii8CAAC4fgVeU2r16tXWv3/33Xfy9fW1PrdYLPrhhx9Uq1atIg3OkR05myJJqkVRCgCAUoe8CAAA4MYVuCjVv39/SZLJZNJDDz1k85qrq6tq1aqlN954o0iDc2SHL4+UqlWJohQAAKUNeREAAMCNK3BRKjMzU1LWDnjbt29XQEBAsQWF/42UqlnJ086RAACAq5EXAQAA3LgCF6WyRUdHF0ccuMrhs4yUAgCgtCMvAgAAuH6FXuhckjZt2qQ+ffqobt26qlevnvr27avNmzcXdWwOKyUtQ6eSUiVRlAIAoLQjLwIAALg+hS5Kffzxx+ratas8PT01btw4jR07Vh4eHrrjjjv06aefFkeMDid76p6fp6t8PV3tHA0AAMgLeREAAMD1MxmGYRTmhIYNG+qRRx7RE088YXN81qxZev/997V///4iDbCoJSYmytfXVwkJCfLx8bF3OLla+0esxny8U81D/PTlYx3tHQ4AAOVOUeUDZT0vKg5lIdcCAADFq6D5QKFHSh06dEh9+vTJcbxv376sq1BEDl8eKVWLRc4BACjVyIsAAACuX6GLUiEhIfrhhx9yHP/hhx8UEhJSJEE5uiOXFzmvyXpSAACUauRFAAAA16/Au+89/PDDeuutt/Tkk09q3Lhx2r17tzp06CCTyaSff/5ZS5Ys0VtvvVWcsTqMw2cYKQUAQGlGXgQAAHDjCrymlLOzs2JjY1WlShWtXLlSb7zxhnWdhIYNG+rpp59Wv379ijXYolAW1jnoMOMHnUi4pC/+3UGtala0dzgAAJQ7N5oPlJe8qDiUhVwLAAAUr4LmAwUeKXVl7equu+7SXXfddWMRIleX0i2KTbwkiZFSAACUVuRFAAAAN65Qa0qZTKbiigOXHTufIsOQKphd5O/lZu9wAABAHsiLAAAAbkyBR0pJ0k033XTNBOzcuXM3FJCjy15PqmYlT5JdAABKMfIiAACAG1OootTUqVPl6+tbXLFA0uHLO+/VYuc9AABKNfIiAACAG1OootT999+vKlWqFFcskHTk7P9GSgEAgNKLvAgAAODGFHhNKaaSlQxGSgEAUPqRFwEAANy4AhelrtxlBsWHkVIAAJR+5EUAAAA3rsDT9zIzM4szDkhKy8jUsfNZRalaAYyUAgCgtCIvAgAAuHEFHimF4nc8/qIyDcnd1UlVvM32DgcAAAAAAKDYUJQqRa5cT4q1KgAAAAAAQHlGUaoUOXImqyhVw5/1pAAAAAAAQPlGUaoUiTl3URKLnAMAAAAAgPKPolQpEnMua5FzRkoBAAAAAIDyjqJUKZK98151ilIAAAAAAKCcoyhVShiGwUgpAAAAAADgMChKlRLnktOUkmaRJFXz87BzNAAAAAAAAMWLolQpkT1KKsjHXe6uznaOBgAAAAAAoHhRlColjp7P2nkvxJ9RUgAAAAAAoPyjKFVKHL08UiqE9aQAAAAAAIADoChVShxlkXMAAAAAAOBAKEqVEtlrSoVUpCgFAAAAAADKP4pSpcTR85dHSlWiKAUAAAAAAMo/ilKlQIYlUyfiL0lipBQAAAAAAHAMFKVKgdiES7JkGnJzcVIVb7O9wwEAAAAAACh2FKVKgez1pKpX9JCTk8nO0QAAAAAAABQ/ilKlADvvAQAAAAAAR0NRqhRg5z0AAAAAAOBoKEqVAkfPX5TESCkAAAAAAOA4KEqVAtaRUv4edo4EAACUBefPn9eQIUPk6+srX19fDRkyRPHx8fmeYxiGIiIiFBwcLA8PD3Xu3Fl79+61vn7u3Dk9/vjjql+/vjw9PVWjRg2NGzdOCQkJxdwbAADgqChKlQLHrEUpRkoBAIBrGzx4sHbv3q21a9dq7dq12r17t4YMGZLvOa+++qpmzZqlOXPmaPv27QoKCtKdd96ppKQkSdKJEyd04sQJvf766/r999+1ZMkSrV27ViNGjCiJLgEAAAdkMgzDsHcQJSkxMVG+vr5KSEiQj4+PvcNRcmqGGk/5TpL0W0Q3+bi72jkiAADKv9KWDxTG/v371ahRI23dulVt27aVJG3dulXt27fXn3/+qfr16+c4xzAMBQcHa8KECZo0aZIkKTU1VYGBgZo5c6ZGjx6d63t9/vnnevDBB5WcnCwXF5cCxVeW7y0AACgaBc0HGCllZ0fPZ42S8vN0pSAFAACuKTIyUr6+vtaClCS1a9dOvr6+2rJlS67nREdHKy4uTt26dbMeM5vNuu222/I8R5I1kcyvIJWamqrExESbBwAAQEFQlLKzmLPsvAcAAAouLi5OVapUyXG8SpUqiouLy/McSQoMDLQ5HhgYmOc5Z8+e1UsvvZTnKKpsM2bMsK5t5evrq5CQkIJ0AwAAgKKUvbHzHgAAkKSIiAiZTKZ8H1FRUZIkk8mU43zDMHI9fqWrX8/rnMTERPXq1UuNGjXSlClT8r3m5MmTlZCQYH0cPXr0Wl0FAACQJBVscQAUm6OXFzmvzs57AAA4tLFjx+r+++/Pt02tWrX022+/6eTJkzleO336dI6RUNmCgoIkZY2Yqlq1qvX4qVOncpyTlJSkHj16qEKFClq5cqVcXfNfXsBsNstsNufbBgAAIDd2Hyk1b948hYaGyt3dXa1atdLmzZvzbT937lw1bNhQHh4eql+/vpYuXVpCkRaPY5dHSjF9DwAAxxYQEKAGDRrk+3B3d1f79u2VkJCgbdu2Wc/99ddflZCQoA4dOuR67dDQUAUFBWn9+vXWY2lpadq0aZPNOYmJierWrZvc3Ny0evVqubu7F1+HAQCAw7NrUWr58uWaMGGCnnvuOe3atUudOnVSWFiYYmJicm0/f/58TZ48WREREdq7d6+mTp2qxx57TF999VUJR150jl1e6LxaRUZKAQCAa2vYsKF69OihUaNGaevWrdq6datGjRql3r172+y816BBA61cuVJS1rS9CRMm6OWXX9bKlSv1xx9/aNiwYfL09NTgwYMlZY2Q6tatm5KTk7Vw4UIlJiYqLi5OcXFxslgsdukrAAAo3+w6fW/WrFkaMWKERo4cKUmaPXu2vvvuO82fP18zZszI0f6jjz7S6NGjNXDgQElS7dq1tXXrVs2cOVN9+vQp0diLynHrSCmKUgAAoGA++eQTjRs3zrqbXt++fTVnzhybNgcOHFBCQoL1+cSJE3Xx4kU9+uijOn/+vNq2bat169bJ29tbkrRjxw79+uuvkqS6devaXCs6Olq1atUqxh4BAABHZLeiVFpamnbs2KFnnnnG5ni3bt3y3Jo4NTU1xzByDw8Pbdu2Tenp6bmueZCamqrU1FTr89K0TXHCxXQlpWZIkoL9KEoBAICC8ff318cff5xvG8MwbJ6bTCZFREQoIiIi1/adO3fOcQ4AAEBxstv0vTNnzshisRRqa+Lu3bvrgw8+0I4dO2QYhqKiorRo0SKlp6frzJkzuZ5Tmrcpzp66V8nLTZ5urDkPAAAAAAAch90XOi/o1sSS9MILLygsLEzt2rWTq6ur+vXrp2HDhkmSnJ2dcz2nNG9TnD11j/WkAAAAAACAo7FbUSogIEDOzs45RkXltjVxNg8PDy1atEgpKSk6fPiwYmJiVKtWLXl7eysgICDXc8xms3x8fGwepUX2znvVKUoBAAAAAAAHY7eilJubm1q1amWzNbEkrV+/Ps/tjLO5urqqevXqcnZ21meffabevXvLycnug74K7Xj85ZFSrCcFAAAAAAAcjF0XMgoPD9eQIUPUunVrtW/fXu+9955iYmI0ZswYSVlT744fP66lS5dKkg4ePKht27apbdu2On/+vGbNmqU//vhDH374oT27cd2y15SqXtHTzpEAAAAAAACULLsWpQYOHKizZ89q2rRpio2NVZMmTbRmzRrVrFlTkhQbG6uYmBhre4vFojfeeEMHDhyQq6urunTpoi1btpTZLYoZKQUAAAAAAByVyXCwvX8TExPl6+urhIQEu68v1WLaOsWnpGvthE5qEFR61roCAKC8K035QHnDvQUAAAXNB8reQkzlxIXUDMWnpEtipBQAAAAAAHA8FKXs5Pjlnfd8PVzl7e5q52gAAAAAAABKFkUpO/nfIueMkgIAAAAAAI6HopSdsMg5AAAAAABwZBSl7OTY5el71St62jkSAAAAAACAkkdRyk6y15SqxvQ9AAAAAADggChK2QlrSgEAAAAAAEdGUcpOWFMKAAAAAAA4MopSdnAxzaIzF9IkSSGsKQUAAAAAABwQRSk7yB4lVcHsIh8PFztHAwAAAAAAUPIoStnBletJmUwmO0cDAAAAAABQ8ihK2QHrSQEAAAAAAEdHUcoOjp3PKkqx8x4AAAAAAHBUFKXsILsoVY2iFAAAAAAAcFAUpezg+OU1par5sfMeAAAAAABwTBSl7OBE/CVJjJQCAAAAAACOi6JUCUu3ZOpkUlZRKtjP3c7RAAAAAAAA2AdFqRIWl3BJhiG5OTspwMts73AAAAAAAADsgqJUCTsRn7XIeVU/dzk5mewcDQAAAAAAgH1QlCphJxKyilLBvqwnBQAAAAAAHBdFqRLGIucAAAAAAAAUpUrc8cvT94L9KEoBAAAAAADHRVGqhGWvKVWNnfcAAAAAAIADoyhVwk4wUgoAAAAAAICiVEkyDEPHz1OUAgAAAAAAoChVghIvZSg5zSKJ3fcAAAAAAIBjoyhVgrKn7vl7ucnDzdnO0QAAAAAAANgPRakS9L+peyxyDgAAAAAAHBtFqRJ0IuFyUYqpewAAAAAAwMFRlCpBx9l5DwAAAAAAQBJFqRJ1Iv6SJKkaRSkAAAAAAODgKEqVoBOMlAIAAAAAAJBEUapEZRelqlWkKAUAAAAAABwbRakSkm7J1MnErOl77L4HAAAAAAAcHUWpEnIy8ZIyDcnN2UkBXmZ7hwMAAAAAAGBXFKVKSPYi51X93OXkZLJzNAAAAAAAAPZFUaqEWBc592U9KQAAAAAAAIpSJeQ4O+8BAAAAAABYUZQqIdad91jkHAAAAAAAgKJUSTnBSCkAAAAAAAArilIlhOl7AAAAAAAA/0NRqoTEXt59L5jpewAAAAAAABSlSkLSpXQlpWZIkqqy+x4AAAAAAABFqZIQm5A1SsrH3UVeZhc7RwMAAAAAAGB/FKVKAIucAwAAAAAA2KIoVQLiLo+UqurLelIAAAAAAAASRakSceJyUSqI9aQAAAAAAAAkUZQqEbHZ0/cYKQUAAAAAACCpFBSl5s2bp9DQULm7u6tVq1bavHlzvu0/+eQTNW/eXJ6enqpataqGDx+us2fPllC01ycu8fL0PdaUAgAAAAAAkGTnotTy5cs1YcIEPffcc9q1a5c6deqksLAwxcTE5Nr+559/1tChQzVixAjt3btXn3/+ubZv366RI0eWcOSFc4KRUgAAAAAAADbsWpSaNWuWRowYoZEjR6phw4aaPXu2QkJCNH/+/Fzbb926VbVq1dK4ceMUGhqqW265RaNHj1ZUVFQJR15whmEo1rqmFEUpAAAAAAAAyY5FqbS0NO3YsUPdunWzOd6tWzdt2bIl13M6dOigY8eOac2aNTIMQydPntR///tf9erVK8/3SU1NVWJios2jJCVezFBKmkWSVJWFzgEAAAAAACTZsSh15swZWSwWBQYG2hwPDAxUXFxcrud06NBBn3zyiQYOHCg3NzcFBQXJz89P77zzTp7vM2PGDPn6+lofISEhRdqPa4lNzJq6V9HTVR5uziX63gAAAAAAAKWV3Rc6N5lMNs8Nw8hxLNu+ffs0btw4vfjii9qxY4fWrl2r6OhojRkzJs/rT548WQkJCdbH0aNHizT+a4mNv7zIOaOkAAAAAAAArFzs9cYBAQFydnbOMSrq1KlTOUZPZZsxY4Y6duyop59+WpLUrFkzeXl5qVOnTpo+fbqqVq2a4xyz2Syz2Vz0HSigEwlZI6Wqsp4UAAAAAACAld1GSrm5ualVq1Zav369zfH169erQ4cOuZ6TkpIiJyfbkJ2ds6bEGYZRPIHeIOtIKT+KUgAAAAAAANnsOn0vPDxcH3zwgRYtWqT9+/friSeeUExMjHU63uTJkzV06FBr+z59+mjFihWaP3++Dh06pF9++UXjxo1TmzZtFBwcbK9u5Ct75z2m7wEAgKJy/vx5DRkyxLpm5pAhQxQfH5/vOYZhKCIiQsHBwfLw8FDnzp21d+/ePNuGhYXJZDJp1apVRd8BAAAA2XH6niQNHDhQZ8+e1bRp0xQbG6smTZpozZo1qlmzpiQpNjZWMTEx1vbDhg1TUlKS5syZoyeffFJ+fn66/fbbNXPmTHt14ZpiL0/fC2akFAAAKCKDBw/WsWPHtHbtWknSI488oiFDhuirr77K85xXX31Vs2bN0pIlS3TTTTdp+vTpuvPOO3XgwAF5e3vbtJ09e3aea3wCAAAUFZNRWue9FZPExET5+voqISFBPj4+xf5+XV7fqOgzyVo2qp3a16lU7O8HAACuraTzgaK0f/9+NWrUSFu3blXbtm0lSVu3blX79u31559/qn79+jnOMQxDwcHBmjBhgiZNmiRJSk1NVWBgoGbOnKnRo0db2+7Zs0e9e/fW9u3bVbVqVa1cuVL9+/cvcHxl+d4CAICiUdB8wO6775VnhmHoRDwjpQAAQNGJjIyUr6+vtSAlSe3atZOvr6+2bNmS6znR0dGKi4tTt27drMfMZrNuu+02m3NSUlI0aNAgzZkzR0FBQQWKJzU1VYmJiTYPAACAgqAoVYziU9KVmpEpSQpi9z0AAFAE4uLiVKVKlRzHq1SpkmNX4yvPkZRjh+PAwECbc5544gl16NBB/fr1K3A8M2bMsK5t5evrq5CQkAKfCwAAHBtFqWJ04vJ6UgEV3GR2cbZzNAAAoDSLiIiQyWTK9xEVFSVJua73ZBjGNdeBuvr1K89ZvXq1NmzYoNmzZxcq7smTJyshIcH6OHr0aKHOBwAAjsuuC52Xd7HxWTvvMUoKAABcy9ixY3X//ffn26ZWrVr67bffdPLkyRyvnT59OsdIqGzZU/Hi4uJUtWpV6/FTp05Zz9mwYYP++ecf+fn52Zw7YMAAderUSRs3bsz12mazWWazOd+4AQAAckNRqhhl77xX1dfDzpEAAIDSLiAgQAEBAdds1759eyUkJGjbtm1q06aNJOnXX39VQkKCOnTokOs5oaGhCgoK0vr169WyZUtJUlpamjZt2mTdxfiZZ57RyJEjbc5r2rSp3nzzTfXp0+dGugYAAJArilLFKDYha6RUMCOlAABAEWnYsKF69OihUaNG6d1335UkPfLII+rdu7fNznsNGjTQjBkzdNddd8lkMmnChAl6+eWXVa9ePdWrV08vv/yyPD09NXjwYElZo6lyW9y8Ro0aCg0NLZnOAQAAh0JRqhhlF6Wq+jFSCgAAFJ1PPvlE48aNs+6m17dvX82ZM8emzYEDB5SQkGB9PnHiRF28eFGPPvqozp8/r7Zt22rdunXy9vYu0dgBAACyUZQqRifis6fvMVIKAAAUHX9/f3388cf5tjEMw+a5yWRSRESEIiIiCvw+V18DAACgKLH7XjGyjpRiTSkAAAAAAAAbFKWKiWEYirMWpRgpBQAAAAAAcCWKUsXkbHKa0iyZMpmkQB+KUgAAAAAAAFeiKFVMYuOzRkkFVDDLzYXbDAAAAAAAcCWqJcUkNiFrkfNgpu4BAAAAAADkQFGqmLDIOQAAAAAAQN4oShWTE5dHSgUxUgoAAAAAACAHilLFJHtNqWA/ilIAAAAAAABXoyhVTOKYvgcAAAAAAJAnilLFJHv6HiOlAAAAAAAAcqIoVQwyMw2dTMwaKRXESCkAAAAAAIAcKEoVgzMXUpVuMeRkkgK9zfYOBwAAAAAAoNShKFUMYi+vJ1XF210uztxiAAAAAACAq1ExKQaxl9eTqsp6UgAAAAAAALmiKFUMTsRn77xHUQoAAAAAACA3FKWKgXWkFIucAwAAAAAA5IqiVDHIXlOKkVIAAAAAAAC5oyhVDLKLUsF+jJQCAAAAAADIDUWpYhAbnzV9L4iRUgAAAAAAALmiKFXELJmGTialSpKCWVMKAAAAAAAgVxSlitjppFRZMg25OJlU2dts73AAAAAAAABKJYpSRezE5Z33An3c5exksnM0AAAAAAAApRNFqSIWG5+1yDnrSQEAAAAAAOSNolQRi708UqoqRSkAAAAAAIA8UZQqYrEJWSOlgv1Y5BwAAAAAACAvFKWKWPZIqSAfRkoBAAAAAADkhaJUETsRnz1SiqIUAAAAAABAXihKFbG4y9P3qvoyfQ8AAAAAACAvFKWKUIYlU6eSLhelGCkFAAAAAACQJ4pSRehkUqoyDcnV2aQAL7O9wwEAAAAAACi1KEoVodj4rEXOA33c5eRksnM0AAAAAAAApZeLvQMoT+oHeevTkW2Vasm0dygAAAAAAAClGkWpIuTt7qoOdQPsHQYAAAAAAECpx/Q9AAAAAAAAlDiKUgAAAAAAAChxFKUAAAAAAABQ4ihKAQAAAAAAoMRRlAIAAAAAAECJoygFAAAAAACAEkdRCgAAAAAAACXO7kWpefPmKTQ0VO7u7mrVqpU2b96cZ9thw4bJZDLleDRu3LgEIwYAAAAAAMCNsmtRavny5ZowYYKee+457dq1S506dVJYWJhiYmJybf/WW28pNjbW+jh69Kj8/f117733lnDkAAAAAAAAuBF2LUrNmjVLI0aM0MiRI9WwYUPNnj1bISEhmj9/fq7tfX19FRQUZH1ERUXp/PnzGj58eAlHDgAAAAAAgBtht6JUWlqaduzYoW7dutkc79atm7Zs2VKgayxcuFBdu3ZVzZo182yTmpqqxMREmwcAAAAAAADsy25FqTNnzshisSgwMNDmeGBgoOLi4q55fmxsrL799luNHDky33YzZsyQr6+v9RESEnJDcQMAAAAAAODG2X2hc5PJZPPcMIwcx3KzZMkS+fn5qX///vm2mzx5shISEqyPo0eP3ki4AAAAAAAAKAIu9nrjgIAAOTs75xgVderUqRyjp65mGIYWLVqkIUOGyM3NLd+2ZrNZZrP5huMFAAAAAABA0bHbSCk3Nze1atVK69evtzm+fv16dejQId9zN23apL///lsjRowozhABAAAAAABQTOw2UkqSwsPDNWTIELVu3Vrt27fXe++9p5iYGI0ZM0ZS1tS748ePa+nSpTbnLVy4UG3btlWTJk0K/Z6GYUgSC54DAODAsvOA7LwARYdcCwAAFDTXsmtRauDAgTp79qymTZum2NhYNWnSRGvWrLHuphcbG6uYmBibcxISEvTFF1/orbfeuq73TEpKkiQWPAcAAEpKSpKvr6+9wyhXyLUAAEC2a+VaJsPBfkWYmZmpEydOyNvbu0ALqhdWYmKiQkJCdPToUfn4+BT59Usz+k7f6bvjoO+O1/fy1m/DMJSUlKTg4GA5Odl935dyhVyr+Dhq3x213xJ9p+/03ZGUt74XNNey60gpe3ByclL16tWL/X18fHzKxTfS9aDv9N3R0Hf67kjKU78ZIVU8yLWKn6P23VH7LdF3+u546Hv56HtBci1+NQgAAAAAAIASR1EKAAAAAAAAJY6iVBEzm82aMmWKzGazvUMpcfSdvjsa+k7fHYmj9huljyN/Lzpq3x213xJ9p+/03ZE4at8dbqFzAAAAAAAA2B8jpQAAAAAAAFDiKEoBAAAAAACgxFGUAgAAAAAAQImjKFWE5s2bp9DQULm7u6tVq1bavHmzvUMqcjNmzNC//vUveXt7q0qVKurfv78OHDhg08YwDEVERCg4OFgeHh7q3Lmz9u7da6eIi8+MGTNkMpk0YcIE67Hy3Pfjx4/rwQcfVKVKleTp6akWLVpox44d1tfLa98zMjL0/PPPKzQ0VB4eHqpdu7amTZumzMxMa5vy0veffvpJffr0UXBwsEwmk1atWmXzekH6mZqaqscff1wBAQHy8vJS3759dezYsRLsxfXJr+/p6emaNGmSmjZtKi8vLwUHB2vo0KE6ceKEzTXKY9+vNnr0aJlMJs2ePdvmeFntO8oecq0s5eXnTn4cLc+SyLXItci1ymOuRZ51bRSlisjy5cs1YcIEPffcc9q1a5c6deqksLAwxcTE2Du0IrVp0yY99thj2rp1q9avX6+MjAx169ZNycnJ1javvvqqZs2apTlz5mj79u0KCgrSnXfeqaSkJDtGXrS2b9+u9957T82aNbM5Xl77fv78eXXs2FGurq769ttvtW/fPr3xxhvy8/OztimvfZ85c6YWLFigOXPmaP/+/Xr11Vf12muv6Z133rG2KS99T05OVvPmzTVnzpxcXy9IPydMmKCVK1fqs88+088//6wLFy6od+/eslgsJdWN65Jf31NSUrRz50698MIL2rlzp1asWKGDBw+qb9++Nu3KY9+vtGrVKv36668KDg7O8VpZ7TvKFnItx8m1HC3Pksi1yLWykGuVv1yLPKsADBSJNm3aGGPGjLE51qBBA+OZZ56xU0Ql49SpU4YkY9OmTYZhGEZmZqYRFBRkvPLKK9Y2ly5dMnx9fY0FCxbYK8wilZSUZNSrV89Yv369cdtttxnjx483DKN8933SpEnGLbfckufr5bnvvXr1Mh5++GGbY3fffbfx4IMPGoZRfvsuyVi5cqX1eUH6GR8fb7i6uhqfffaZtc3x48cNJycnY+3atSUW+426uu+52bZtmyHJOHLkiGEY5b/vx44dM6pVq2b88ccfRs2aNY0333zT+lp56TtKP3Itx8i1HDHPMgxyLXItcq2rlcdcizwrd4yUKgJpaWnasWOHunXrZnO8W7du2rJli52iKhkJCQmSJH9/f0lSdHS04uLibO6F2WzWbbfdVm7uxWOPPaZevXqpa9euNsfLc99Xr16t1q1b695771WVKlXUsmVLvf/++9bXy3Pfb7nlFv3www86ePCgJGnPnj36+eef1bNnT0nlu+9XKkg/d+zYofT0dJs2wcHBatKkSbm6F1LW/30mk8n6G+zy3PfMzEwNGTJETz/9tBo3bpzj9fLcd5Qe5FqOk2s5Yp4lkWuRa5FrXc1Rci3yLMnF3gGUB2fOnJHFYlFgYKDN8cDAQMXFxdkpquJnGIbCw8N1yy23qEmTJpJk7W9u9+LIkSMlHmNR++yzz7Rz505t3749x2vlue+HDh3S/PnzFR4ermeffVbbtm3TuHHjZDabNXTo0HLd90mTJikhIUENGjSQs7OzLBaL/vOf/2jQoEGSyvfX/UoF6WdcXJzc3NxUsWLFHG3K0/+Fly5d0jPPPKPBgwfLx8dHUvnu+8yZM+Xi4qJx48bl+np57jtKD3Itx8i1HDXPksi1yLXIta7kSLkWeRZFqSJlMplsnhuGkeNYeTJ27Fj99ttv+vnnn3O8Vh7vxdGjRzV+/HitW7dO7u7uebYrj33PzMxU69at9fLLL0uSWrZsqb1792r+/PkaOnSotd3/t3f/MVXVfxzHX5dfyg9FAeNCBCgNJHQmNIssndoPM1OmS2jOINaIFqmNcDlZ4h9Zf7lZy802wj/8EW25lmYqKr9arBZyEXOp1AUtYayfZig/vJ/vH309RYLZ9ysXODwf253ecz73nM/7MM597X0P59qx9vLycu3cuVO7d+9WSkqKXC6X1q5dq+joaGVnZ1vj7Fh7f/6XOu10LHp6epSVlSWPx6Nt27b94/iRXnt9fb22bt2q48eP/+s6RnrtGJ5Gy7n2mtGUtUZzzpLIWmStP5G1Rk/WImf9gT/fuwUiIiLk6+t7Xaeyo6Pjuk63Xbz44ov66KOPVFlZqZiYGGu50+mUJFsei/r6enV0dCgtLU1+fn7y8/NTdXW13nzzTfn5+Vn12bH2qKgo3XXXXX2WJScnWzeXtfPPvaioSK+88oqysrI0ffp0rVq1Si+99JJef/11Sfau/a9upk6n06nu7m79/PPPA44ZyXp6erRixQq53W5VVFRYn9xJ9q29trZWHR0dio2Ntc57ra2tKiwsVHx8vCT71o7hhaxl/6w1mnOWRNYia5G1pNGXtchZf6ApdQsEBAQoLS1NFRUVfZZXVFTo/vvvH6JZDQ5jjAoKCrR3714dO3ZMkydP7rN+8uTJcjqdfY5Fd3e3qqurR/yxWLBggZqamuRyuazHPffco5UrV8rlcmnKlCm2rX327NnXfR31mTNnFBcXJ8neP/fOzk75+PQ9Vfr6+lpfU2zn2v/qZupMS0uTv79/nzFtbW06efLkiD8W10LS2bNndeTIEYWHh/dZb9faV61apRMnTvQ570VHR6uoqEiHDh2SZN/aMbyQtf5k1/ed0ZyzJLIWWYusNRqzFjnrv7x5V3U7e++994y/v78pLS01p06dMmvXrjXBwcGmpaVlqKd2Sz3//PMmNDTUVFVVmba2NuvR2dlpjXnjjTdMaGio2bt3r2lqajJPPfWUiYqKMhcvXhzCmQ+Ov34rjDH2rf2LL74wfn5+5rXXXjNnz541u3btMkFBQWbnzp3WGLvWnp2dbW6//Xazf/9+43a7zd69e01ERIRZt26dNcYutf/222+moaHBNDQ0GElmy5YtpqGhwfrWk5upMz8/38TExJgjR46Y48ePm/nz55sZM2aY3t7eoSrrptyo9p6eHrNkyRITExNjXC5Xn3NfV1eXtQ071t6fv38rjDEjt3aMLGSt0Ze1RkvOMoasRdYia9k1a5Gz/hlNqVvo7bffNnFxcSYgIMCkpqZaX91rJ5L6fZSVlVljPB6P2bhxo3E6nWbMmDFmzpw5pqmpaegmPYj+HpbsXPu+ffvMtGnTzJgxY8zUqVPNO++802e9XWu/ePGiWbNmjYmNjTVjx441U6ZMMRs2bOjzBmmX2isrK/v9/c7OzjbG3Fydly9fNgUFBSYsLMwEBgaaxYsXm3Pnzg1BNf/OjWp3u90DnvsqKyutbdix9v70F5ZGau0Yechaf7DL+84/GU05yxiyFlmLrGXHrEXO+mcOY4y5NddcAQAAAAAAADeHe0oBAAAAAADA62hKAQAAAAAAwOtoSgEAAAAAAMDraEoBAAAAAADA62hKAQAAAAAAwOtoSgEAAAAAAMDraEoBAAAAAADA62hKAQAAAAAAwOtoSgEYciUlJbr77ruHehoAAAC2Q84CMJzRlAIwqBwOxw0fOTk5evnll3X06NEhmd8HH3yge++9V6GhoRo3bpxSUlJUWFhorSfIAQCA4YqcBWCk8xvqCQCwt7a2Nuv/5eXlevXVV3X69GlrWWBgoEJCQhQSEuL1uR05ckRZWVnavHmzlixZIofDoVOnTg1ZcAMAAPg3yFkARjqulAIwqJxOp/UIDQ2Vw+G4btnfPyXLyclRRkaGNm/erMjISE2YMEGbNm1Sb2+vioqKFBYWppiYGL377rt99vX9998rMzNTEydOVHh4uJYuXaqWlpYB57Z//3498MADKioqUlJSkhITE5WRkaG33npLkrRjxw5t2rRJjY2N1ieOO3bskCT9+uuvysvL02233abx48dr/vz5amxstLZ9rabt27frjjvuUFBQkJ588kn98ssv1piqqirNmjVLwcHBmjBhgmbPnq3W1tb/+5gDAIDRgZxFzgJGOppSAIalY8eO6cKFC6qpqdGWLVtUUlKixYsXa+LEifr888+Vn5+v/Px8nT9/XpLU2dmpefPmKSQkRDU1Nfr0008VEhKihQsXqru7u999OJ1OffXVVzp58mS/6zMzM1VYWKiUlBS1tbWpra1NmZmZMsbo8ccfV3t7uw4cOKD6+nqlpqZqwYIF+umnn6zXNzc36/3339e+fft08OBBuVwuvfDCC5Kk3t5eZWRkaO7cuTpx4oTq6uqUl5cnh8Nxi48kAABAX+QsAMOGAQAvKSsrM6Ghodct37hxo5kxY4b1PDs728TFxZmrV69ay5KSksyDDz5oPe/t7TXBwcFmz549xhhjSktLTVJSkvF4PNaYrq4uExgYaA4dOtTvfC5dumQWLVpkJJm4uDiTmZlpSktLzZUrVwacmzHGHD161IwfP77POGOMSUhIMNu3b7de5+vra86fP2+t/+STT4yPj49pa2szP/74o5FkqqqqBjhaAAAAN4+cRc4CRiKulAIwLKWkpMjH589TVGRkpKZPn2499/X1VXh4uDo6OiRJ9fX1am5u1rhx46x7J4SFhenKlSv65ptv+t1HcHCwPv74YzU3N6u4uFghISEqLCzUrFmz1NnZOeDc6uvrdenSJYWHh1v7CgkJkdvt7rOv2NhYxcTEWM/T09Pl8Xh0+vRphYWFKScnR48++qieeOIJbd26tc99IQAAAAYLOQvAcMGNzgEMS/7+/n2eOxyOfpd5PB5JksfjUVpamnbt2nXdtiZNmnTDfSUkJCghIUHPPvusNmzYoMTERJWXl+uZZ57pd7zH41FUVJSqqqquWzdhwoQB93PtkvFr/5aVlWn16tU6ePCgysvLVVxcrIqKCt133303nC8AAMD/g5wFYLigKQXAFlJTU1VeXm7dEPN/FR8fr6CgIP3++++SpICAAF29evW6fbW3t8vPz0/x8fEDbuvcuXO6cOGCoqOjJUl1dXXy8fFRYmKiNWbmzJmaOXOm1q9fr/T0dO3evZuwBAAAhhVyFoDBwp/vAbCFlStXKiIiQkuXLlVtba3cbreqq6u1Zs0afffdd/2+pqSkROvWrVNVVZXcbrcaGhqUm5urnp4ePfzww5L+CE9ut1sul0s//PCDurq69NBDDyk9PV0ZGRk6dOiQWlpa9Nlnn6m4uFhffvmltf2xY8cqOztbjY2Nqq2t1erVq7VixQo5nU653W6tX79edXV1am1t1eHDh3XmzBklJyd75XgBAADcLHIWgMFCUwqALQQFBammpkaxsbFatmyZkpOTlZubq8uXLw/4id7cuXP17bff6umnn9bUqVP12GOPqb29XYcPH1ZSUpIkafny5Vq4cKHmzZunSZMmac+ePXI4HDpw4IDmzJmj3NxcJSYmKisrSy0tLYqMjLS2f+edd2rZsmVatGiRHnnkEU2bNk3btm2z5vv1119r+fLlSkxMVF5engoKCvTcc88N/sECAAD4F8hZAAaLwxhjhnoSAGA3JSUl+vDDD+VyuYZ6KgAAALZCzgLsgyulAAAAAAAA4HU0pQAAAAAAAOB1/PkeAAAAAAAAvI4rpQAAAAAAAOB1NKUAAAAAAADgdTSlAAAAAAAA4HU0pQAAAAAAAOB1NKUAAAAAAADgdTSlAAAAAAAA4HU0pQAAAAAAAOB1NKUAAAAAAADgdTSlAAAAAAAA4HX/AUt3j2cEddqCAAAAAElFTkSuQmCC", + "text/plain": [ + "
" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "#plot the results\n", + "import matplotlib.pyplot as plt\n", + "\n", + "plt.figure(figsize=(12, 5))\n", + "\n", + "plt.subplot(1, 2, 1)\n", + "plt.plot(plant_biomass_history, label=\"Plants\")\n", + "plt.xlabel(\"Time Steps\")\n", + "plt.ylabel(\"Total Biomass\")\n", + "plt.title(\"Plant Biomass Dynamics\")\n", + "plt.legend()\n", + "\n", + "plt.subplot(1, 2, 2)\n", + "plt.plot(vertebrate_biomass_history, label=\"Vertebrates\")\n", + "plt.xlabel(\"Time Steps\")\n", + "plt.ylabel(\"Total Biomass\")\n", + "plt.title(\"Vertebrate Biomass Dynamics\")\n", + "plt.legend()\n", + "\n", + "plt.tight_layout()\n", + "plt.show()" + ] + } + ], + "metadata": { + "kernelspec": { + "display_name": ".venv", + "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.12.3" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +} diff --git a/experiments/tuesday_architecture_proposal/engine_v2/EcosystemEngine.py b/experiments/tuesday_architecture_proposal/engine_v2/EcosystemEngine.py deleted file mode 100644 index 5048527..0000000 --- a/experiments/tuesday_architecture_proposal/engine_v2/EcosystemEngine.py +++ /dev/null @@ -1,16 +0,0 @@ -from .EnvironmentState import EnvironmentState -from .EcosystemGridState import EcosystemGridState -class EcosystemEngine: - def __init__(self, grid_state: EcosystemGridState, env_state: EnvironmentState): - self.grid = grid_state - self.env = env_state - self.processes = [] - - def add_process(self, process_func): - """Registers a process to the simulation pipeline.""" - self.processes.append(process_func) - - def step(self): - """Executes one time step by running all processes in order.""" - for process in self.processes: - process(self.grid, self.env) \ No newline at end of file diff --git a/experiments/tuesday_architecture_proposal/engine_v2/EnvironmentState.py b/experiments/tuesday_architecture_proposal/engine_v2/EnvironmentState.py deleted file mode 100644 index 405f59b..0000000 --- a/experiments/tuesday_architecture_proposal/engine_v2/EnvironmentState.py +++ /dev/null @@ -1,22 +0,0 @@ -import numpy as np - -class EnvironmentState: - """ - Represents the physical world grid. - Passed through our modular functions to provide environmental context. - """ - def __init__(self, shape: tuple): - # shape could be (latitude_bins, longitude_bins) e.g., (180, 360) for 1-degree resolution - self.shape = shape - - # Dictionary to hold our vectorized spatial data (e.g., temperature, land/water mask) - self.layers = {} - - def add_layer(self, name: str, data: np.ndarray): - """Adds a new environmental variable layer.""" - if data.shape != self.shape: - raise ValueError(f"Layer shape {data.shape} does not match world shape {self.shape}") - self.layers[name] = data - - def get_layer(self, name: str) -> np.ndarray: - return self.layers[name] \ No newline at end of file diff --git a/experiments/tuesday_architecture_proposal/engine_v2/processes.py b/experiments/tuesday_architecture_proposal/engine_v2/processes.py deleted file mode 100644 index 09a97e2..0000000 --- a/experiments/tuesday_architecture_proposal/engine_v2/processes.py +++ /dev/null @@ -1,33 +0,0 @@ -from .EnvironmentState import EnvironmentState -import numpy as np - -from .EcosystemGridState import EcosystemGridState - -# --- Mock Processes (What your colleagues will write) --- - -def apply_vegetation_growth(grid: EcosystemGridState, env: EnvironmentState): - """Vegetation updates the basal biomass layer.""" - - r = grid.registry.get_group_parameter("plants", "base_growth_rate") - #broadcasting the carrying capacity across the species dimension for plants. - #inital shape of cc is (X, Y) and becomes (X, Y, 1) after adding the new axis, allowing it to broadcast correctly with the plants layer which has shape (X, Y, N_plants). - cc = env.get_layer("carrying_capacity")[..., np.newaxis] - - # The 'with' block handles the extraction and the write-back automatically! - with grid.edit_group_data("biomass", "plants") as plants: - # We do the math in-place on the yielded object - plants += (r * plants * (1 - (plants / cc))) - - - -def apply_atn_step(grid: EcosystemGridState, env: EnvironmentState): - """ATN handles feeding, metabolism, and sets up growth rates for dispersal.""" - # The ATN process modifies biomass based on trophic interactions - # and populates intermediate matrices like grid.net_growth_rate - pass - -def apply_dispersal(grid: EcosystemGridState, env: EnvironmentState): - """Dispersal moves biomass between adjacent grid cells.""" - # Dispersal reads grid.biomass and grid.net_growth_rate, - # calculates diffusion, and shifts biomass left/right/up/down. - pass \ No newline at end of file diff --git a/src/gem.egg-info/PKG-INFO b/src/gem.egg-info/PKG-INFO new file mode 100644 index 0000000..e8f6630 --- /dev/null +++ b/src/gem.egg-info/PKG-INFO @@ -0,0 +1,186 @@ +Metadata-Version: 2.4 +Name: gem +Version: 0.0.0 +Summary: General Ecosystem Model โ€” prototype by the GEM working group. +Author: GEM working group +License: MIT +Requires-Python: >=3.11 +Description-Content-Type: text/markdown +Requires-Dist: numpy +Requires-Dist: scipy +Requires-Dist: pandas +Requires-Dist: matplotlib +Requires-Dist: rasterio +Provides-Extra: dev +Requires-Dist: pytest; extra == "dev" + +# GEM Working Group + +This repository supports the GEM working group focused on designing and prototyping a new General Ecosystem Model. + +The project is inspired by Madingley, but the goal is not to reproduce it directly. The working group is exploring a new model structure that can better represent biodiversity change, species-level dynamics, trophic interactions, spatial structure, dispersal, vegetation, coexistence, and ecosystem processes. + +A second goal is to document how AI coding agents can support collaborative ecological modelling, while keeping scientific assumptions, model structure, tests, and validation explicit. + +## Objectives + +- Design a new process-based ecosystem model. +- Compare the effects of biodiversity change and climate change on ecosystem functioning. +- Explore species-level or hybrid species/guild representations. +- Develop modular model components for trophic dynamics, vegetation, dispersal, mortality, reproduction, and spatial processes. +- Define ecological tests, diagnostics, and validation targets. +- Document lessons from AI-assisted collaborative model development. + +## Resources and documentation + +[๐Ÿ“‚ SharePoint Documents](https://usherbrooke.sharepoint.com/sites/ielabworkinggroup) + +[๐Ÿ’ฌ Teams Chat](https://teams.microsoft.com/l/team/19%3A5UwBzgYI52ESK9znTZVksNqlEEprrzf8AeMWIcRYWpU1%40thread.tacv2/conversations?groupId=f8341ee7-d260-4f7a-9890-6f13bbc5ec80&tenantId=3a5a8744-5935-45f9-9423-b32c3a5de082) + +[๐Ÿ“… Logistics](https://usherbrooke.sharepoint.com/:f:/r/sites/ielabworkinggroup/Documents%20partages/logistics?csf=1&web=1&e=nsdYNf) + +## General Collaboration guidelines + +- This repo is where all notes, papers, and documents for the duration of the row group will live +- The `main` branch is protected, so please create a new branch for any changes and submit a pull request. + +## AI coding assistants + +[AGENTS.md](AGENTS.md) (mirrored as [CLAUDE.md](CLAUDE.md)) tells AI coding agents โ€” Claude Code, Codex, Cursor, โ€ฆ โ€” how to help on this project: how to communicate with you, what conventions to enforce, when to push back proactively, and what they must not do (no inventing process mechanics, no unfamiliar abstractions, no commits or pushes without asking). Read it before you let an agent touch the repo, so you know the behaviour to expect and what to correct. If you disagree with anything in there, open a PR against `AGENTS.md` โ€” it is the source of truth, and `CLAUDE.md` is kept in sync from it. + +## Project structure + +``` +gem-working-group/ +โ”œโ”€โ”€ data/ # Input data files (large rasters gitignored); see Input data files +โ”œโ”€โ”€ docs/ # Reference documents โ€” contracts, specifications, design notes +โ”œโ”€โ”€ experiments/ # Prototypes and experiments; one subfolder per experiment +โ”œโ”€โ”€ papers/ # Reference papers and bibliography +โ”œโ”€โ”€ src/ # The model code package (processes + engine), added when implementation starts +โ””โ”€โ”€ README.md +``` + +Everyday work happens in `experiments/` while the model is being prototyped. Once a process or engine component is stable, it migrates into `src/` so it can be imported by every experiment. + +## Python packaging and environment + +We use Python. The dependencies actually used across the prototype branches are: + +- `numpy` โ€” numerical arrays, the backbone of every process. +- `scipy` โ€” numerical integration for ODE-style processes (`scipy.integrate.solve_ivp`, `odeint`). +- `pandas` โ€” tabular data handling (species traits, parameter tables, diagnostics). +- `matplotlib` โ€” plotting for experiments and figure reproduction. +- `rasterio` โ€” reading and writing GeoTIFF spatial inputs. +- `pytest` โ€” the recommended test runner for unit tests (see *Model development*). + +Additional libraries (e.g. `xarray`, `geopandas`) are added only when a process actually needs them. + +Dependencies and packaging metadata live in a single `pyproject.toml` at the repo root. No `setup.py`, no `requirements.txt` โ€” `pyproject.toml` is the only source of truth, so there is nothing to keep in sync by hand. + +The package lives at [src/gem/](src/gem/) and imports as `import gem` (e.g. `from gem.vegetation import logistic_growth_delta`). The repo directory `gem-working-group/` is the *workspace*, not the package โ€” keep the distinction in mind when writing imports. + +**Getting started.** From the repo root: + +```bash +python -m venv .venv # create an isolated environment +source .venv/bin/activate # macOS / Linux +.venv\Scripts\activate # Windows +pip install -e . # install the project and its dependencies +pytest # confirm everything works +``` + +`pip install -e .` installs the project in **editable mode**: changes you make to the source code take effect immediately without reinstalling. The `.venv/` folder is local to your machine and is gitignored โ€” never commit it. Each team member recreates it from `pyproject.toml`. + +VS Code detects automatically a virtual environment in the project folder and recommends you to activate it when you start working with a virtual environment in your workspace. + +## Style guide and naming conventions + +We follow standard Python conventions (PEP 8) so the code is recognisable to anyone who has read a Python tutorial: + +- **Modules and packages** (folders and `.py` files): lowercase with underscores โ€” `vegetation.py`, `species_registry.py`. No CamelCase, no hyphens. +- **Functions and variables**: lowercase with underscores โ€” `logistic_growth_delta`, `body_mass`. +- **Constants**: uppercase with underscores โ€” `EXTINCTION_THRESHOLD`, `T0_K`. +- **Classes**: CamelCase โ€” `EcosystemGridState`, `EnvironmentState`. + +Names should describe **what** the thing is, not how it is implemented: prefer `metabolic_rate` over `m`, `carrying_capacity` over `K_arr`. Single-letter names are fine inside short math expressions where the meaning is local (`B`, `r`, `K` matching the equation you are encoding) but not in module-level APIs. + +Branch names: lowercase with hyphens โ€” `vegetation-logistic-growth`, `dispersal-density-dependent`. One feature per branch. + + +## Experiments and prototyping + +The `experiments` folder is where we will develop and share code for model prototyping, testing, and experiments. Each experiment should have its own subfolder with a README describing the purpose, methods, and results. + +Naming convention for experiment folders: `DAY_GROUPNAME_experimentNAME`. Example - `sunday_atn_bylot_experiment1`. + +We recommend using Jupyter notebooks for prototyping and documentation, but feel free to use other formats as needed. The key is to keep everything organized and well-documented for future reference. + +## Model development + +The model evolves the ecosystem by composing modular **processes** โ€” vegetation growth, ATN trophic dynamics, dispersal, metabolism, fire, and so on. Each process is implemented as a pure numpy function with a typed signature. + +The full contract is specified in [docs/processes_implementation_specification.md](docs/processes_implementation_specification.md). What every contributor needs to know: + +- **Two process categories.** Biomass-modifying processes (vegetation, ATN, dispersal) return a `biomass_delta` array โ€” the finite change in biomass over one time step `dt`. Dependency processes (metabolism, NPP, ...) return a shared intermediate quantity (rate, flux, factor) consumed by *multiple* biomass-modifying processes; they take no `dt`. +- **Numpy at the process's natural dimensionality.** The same science code runs on a single cell, a row, or the full `(X, Y, S)` grid via standard broadcasting. Per-cell python loops are the main performance trap and are not acceptable. +- **Typed signatures, runtime shape asserts.** All arrays are `NDArray[np.float64]`; scalars are `float`. A one-line shape `assert` at the top of each science function catches broadcast mismatches early. +- **ODE-style processes** (ATN and any other rate-based formulation) pick an integration strategy explicitly โ€” either integrate `dB/dt` to a delta (`scipy.solve_ivp` or a vectorised RK4 step) or rewrite the science to return `biomass_delta` directly. A continuous-time rate is never the public output of a process. + +Adding a new process is the same recipe every time: write a typed science function in its own module (`vegetation.py`, `atn.py`, `dispersal.py`, `metabolism.py`, ...) that imports nothing but `numpy`, add a runtime shape assert, and add a unit test against hand-built arrays. How the function is wired into a running simulation is the engine's concern, covered below. + +## Input data files + +Input data lives in `data/` with two subfolders: + +- `data/raw/`: untouched inputs as downloaded from their source (climate reanalyses, species traits, occurrences, ...). Never edit these by hand. +- `data/processed/`: inputs reprojected onto the engine grid (see *Geographic grid*), cleaned, or otherwise prepared for the engine to consume. + +**Formats.** GeoTIFF (`.tif`) for spatial raster data (temperature maps, carrying capacity, ...). CSV or JSON for tabular data (species traits, parameter sets, ...). Avoid project-specific binary formats โ€” they make data hard to inspect outside the engine. + +**File naming.** Include the date and a short descriptor: `data/raw/era5_temperature_2024-06-01.tif`, `data/processed/species_traits_bylot_v2.csv`. This keeps reruns reproducible because the filename itself records which version was used. + +**Large files.** Large rasters are gitignored. Each `data/raw/` subfolder should include a small script (`download.py` or `download.sh`) that reproduces the download, so the team can rebuild the dataset without committing gigabytes to git. + +**Reprojection** of raw inputs onto the engine's projected grid happens in `data/processed/`, not inside the engine or in processes. The engine consumes already-projected data. + +## Geographic grid + +For all experiments and simulations, we will use a common spatial grid to ensure comparability of results. The grid will cover North America with the following specifications: + +CRS: `ESRI:102008 North America Albers Equal Area Conic` +PROJ: `+proj=aea +lat_0=40 +lon_0=-96 +lat_1=20 +lat_2=60 +x_0=0 +y_0=0 +datum=NAD83 +units=m +no_defs` +Cell size: `100000 m ร— 100000 m` +Grid origin: `x = -7000000 m`, `y = -2000000 m` +Extent: `x = -7000000..5000000`, `y = -2000000..5500000` +Cell ID: `NA100_R{row}_C{col}` from upper-left or lower-left origin, documented explicitly + +## Simulation engine + +The engine is the runtime glue around the modular processes described in *Model development*. It is intentionally small โ€” most of the action is in the science modules โ€” and is built around three shared state objects ("the cart") that travel together through every process, plus a pipeline that runs the processes in order on each time step. + +- **State management.** Three objects hold all of the model's state, and they are the **single source of truth** โ€” processes read from them and write back to them rather than keeping their own copies: + - `EcosystemGridState`: the dynamic `(X, Y, Species)` biological state. Holds a `biomass` layer by default, plus any named layers processes register โ€” the shared `biomass_delta` layer that biomass-modifying processes accumulate into, dependency outputs like `metabolic_rate`, and so on. + - `EnvironmentState`: the `(X, Y)` environmental layers (temperature, carrying capacity, ...). Layers are added by name and shape-checked against the grid. + - `SpeciesRegistry`: the species list, the functional or trophic groups they belong to (`plants`, `herbivores`, ...), per-species traits stored as 1D arrays for fast vectorised math, and the feeding adjacency matrix. +- **Adapters (`processes.py`).** Science modules import nothing from the engine. All engine glue lives in a single file, `processes.py`, holding one `apply_(grid, env, dt)` per process. An adapter slices the right state arrays, fetches parameters, calls the science function, and writes the result back โ€” biomass-modifying adapters accumulate into the shared delta layer; dependency adapters write to a named shared layer. +- **Broadcasting.** Every process operates on numpy arrays at its natural shape โ€” `(S,)` for a single cell, `(Y, S)` for a row, `(X, Y, S)` for the full grid. Adapters reshape per-species parameters to `(1, 1, S)` and environmental layers to `(X, Y, 1)` so they broadcast cleanly against the `(X, Y, S)` biomass array. The same science code then runs unchanged from a unit test on a `(3,)` array to a global simulation. +- **Initialization.** Setting up a run means building the three state objects: define the grid extent and projection (see *Geographic grid*), load environmental layers from `data/processed/`, load the species list and traits, set initial biomass. Initialization is a regular Python function โ€” not a config file โ€” so it can be parametrised and reused across experiments. Each experiment owns its own initialization script. +- **Pipeline registration.** Each adapter is registered with the engine in pipeline order via `engine.add_process(apply_metabolism)`. Dependency-process adapters must be registered before any adapter that consumes their output. Removing a process is deleting one line. + +The end-of-step integration applies the accumulated `biomass_delta` to the biomass layer and zeroes the delta. This makes within-step computation order-independent (every process sees the same `state_t`) and matches how the underlying equations are usually written: a sum of contributions. + +## Running simulations + +Simulations live in `experiments/`, one subfolder per experiment (see *Experiments and prototyping*). Each experiment runs the engine for a defined set of conditions and stores its outputs locally. + +**Minimum reproducibility checklist.** Every run should record: + +- The **random seed** used (if any stochastic process is involved). +- The **configuration**: initial state, environmental data files, species list, parameter values, number of time steps, `dt`. A small JSON or YAML file alongside the outputs is sufficient. +- The **engine version**: the git commit hash the run was produced from. + +Without these three, a result cannot be reproduced. + +**Output storage.** Save outputs inside the experiment folder under a timestamped run name, e.g. `experiments/sunday_atn_bylot_experiment1/runs/2026-05-27_baseline/`. The folder should contain the biomass trajectory (NetCDF or `.npz` for arrays, CSV for tabular diagnostics), the configuration file, and any plots. Large outputs are gitignored; commit only what is small and informative (configuration, diagnostics, key plots). + +**Sharing runs.** Notebooks are useful for exploring outputs but should not be the canonical record. Keep the data files self-describing (column names, units in metadata) so a teammate can reopen them without needing your notebook. diff --git a/src/gem.egg-info/SOURCES.txt b/src/gem.egg-info/SOURCES.txt new file mode 100644 index 0000000..b9e5dc2 --- /dev/null +++ b/src/gem.egg-info/SOURCES.txt @@ -0,0 +1,15 @@ +README.md +pyproject.toml +src/gem/__init__.py +src/gem/vegetation.py +src/gem.egg-info/PKG-INFO +src/gem.egg-info/SOURCES.txt +src/gem.egg-info/dependency_links.txt +src/gem.egg-info/requires.txt +src/gem.egg-info/top_level.txt +src/gem/engine/ecosystem_engine.py +src/gem/engine/ecosystem_grid_state.py +src/gem/engine/environment_state.py +src/gem/engine/processes.py +src/gem/engine/species_registry.py +tests/test_vegetation.py \ No newline at end of file diff --git a/src/gem.egg-info/dependency_links.txt b/src/gem.egg-info/dependency_links.txt new file mode 100644 index 0000000..8b13789 --- /dev/null +++ b/src/gem.egg-info/dependency_links.txt @@ -0,0 +1 @@ + diff --git a/src/gem.egg-info/requires.txt b/src/gem.egg-info/requires.txt new file mode 100644 index 0000000..f7a8098 --- /dev/null +++ b/src/gem.egg-info/requires.txt @@ -0,0 +1,8 @@ +numpy +scipy +pandas +matplotlib +rasterio + +[dev] +pytest diff --git a/src/gem.egg-info/top_level.txt b/src/gem.egg-info/top_level.txt new file mode 100644 index 0000000..6e0ead2 --- /dev/null +++ b/src/gem.egg-info/top_level.txt @@ -0,0 +1 @@ +gem diff --git a/src/gem/engine/ecosystem_engine.py b/src/gem/engine/ecosystem_engine.py new file mode 100644 index 0000000..efa3fda --- /dev/null +++ b/src/gem/engine/ecosystem_engine.py @@ -0,0 +1,34 @@ +from .environment_state import EnvironmentState +from .ecosystem_grid_state import EcosystemGridState + +class EcosystemEngine: + def __init__(self, grid_state: EcosystemGridState, env_state: EnvironmentState): + self.grid = grid_state + self.env = env_state + self.processes = [] + + def add_process(self, process_func): + """Registers a process to the simulation pipeline.""" + self.processes.append(process_func) + + def step(self): + """ + Executes one time step: + 1. All processes accumulate their contributions into registered delta layers. + 2. Integrate all deltas: for each (source_layer, delta_layers) pair, + source_layer += ฮฃ(delta_layers). + 3. Zero all delta layers for the next step. + + This makes computation order-independent: all processes see the same state_t, + and multiple processes can contribute to the same source layer through separate deltas. + """ + # Run all processes; each accumulates into registered delta layers + for process in self.processes: + process(self.grid, self.env) + + # Integrate all registered deltas into their source layers + for source_layer, delta_layer_list in self.grid.delta_layers.items(): + for delta_layer in delta_layer_list: + self.grid.layers[source_layer] += self.grid.layers[delta_layer] + # Reset delta for next step + self.grid.layers[delta_layer][:] = 0.0 \ No newline at end of file diff --git a/experiments/tuesday_architecture_proposal/engine_v2/EcosystemGridState.py b/src/gem/engine/ecosystem_grid_state.py similarity index 62% rename from experiments/tuesday_architecture_proposal/engine_v2/EcosystemGridState.py rename to src/gem/engine/ecosystem_grid_state.py index efed161..67f9c49 100644 --- a/experiments/tuesday_architecture_proposal/engine_v2/EcosystemGridState.py +++ b/src/gem/engine/ecosystem_grid_state.py @@ -1,5 +1,5 @@ import numpy as np -from .SpeciesRegistry import SpeciesRegistry +from .species_registry import SpeciesRegistry from contextlib import contextmanager @@ -15,8 +15,14 @@ def __init__(self, shape: tuple, registry: SpeciesRegistry): # Dictionary to hold all 3D [X, Y, Species] matrices self.layers = {} + # Registry: maps source_layer_name -> list of delta_layer_names. + # Used by the engine to integrate all deltas into their source layers at step end. + self.delta_layers = {} + # The core currency requested by the ATN specs is always initialized self.add_layer("biomass") + # Register a default delta for biomass; processes can add more. + self.add_delta_layer("vegetation_delta", source_layer="biomass") def add_layer(self, layer_name: str, initial_value: float = 0.0): """Creates a new 3D tracking matrix for all cells and species.""" @@ -31,6 +37,28 @@ def add_layer(self, layer_name: str, initial_value: float = 0.0): initial_value, dtype=np.float32 ) + + def add_delta_layer(self, delta_name: str, source_layer: str): + """Register a delta layer that will update a source layer at step end. + + Multiple deltas can be registered for the same source layer. + The engine integrates all registered deltas into the source layer + and resets them to zero after each step. + + Args: + delta_name: Name of the delta layer (will be created if it doesn't exist). + source_layer: Name of the source layer this delta updates. + """ + # Create the delta layer if needed + if delta_name not in self.layers: + self.add_layer(delta_name) + + # Register in the delta registry + if source_layer not in self.delta_layers: + self.delta_layers[source_layer] = [] + + if delta_name not in self.delta_layers[source_layer]: + self.delta_layers[source_layer].append(delta_name) def get_layer_view(self, layer_name: str, group_name: str) -> np.ndarray: """ diff --git a/src/gem/engine/environment_state.py b/src/gem/engine/environment_state.py new file mode 100644 index 0000000..590910a --- /dev/null +++ b/src/gem/engine/environment_state.py @@ -0,0 +1,43 @@ +import numpy as np + +class EnvironmentState: + """ + Represents the physical world grid and its environmental layers (e.g., Temperature, NPP). + + Following the project's geographic guidelines, this class supports equal-area + projections (e.g., Albers North America) rather than assuming lat-lon. + This ensures that one cell equals one comparable unit of surface area. + """ + def __init__( + self, + shape: tuple, + crs: str = "ESRI:102008", + proj_string: str = "+proj=aea +lat_0=40 +lon_0=-96 +lat_1=20 +lat_2=60 +x_0=0 +y_0=0 +datum=NAD83 +units=m +no_defs", + cell_size: float = 100000.0, + origin: tuple = (-7000000.0, -2000000.0) + ): + """ + Initializes the environment state with metadata for an equal-area grid. + """ + self.shape = shape # (grid_x, grid_y) + self.crs = crs + self.proj_string = proj_string + self.cell_size = cell_size + self.origin = origin + + # Dictionary to hold our vectorized spatial data (e.g., temperature, land/water mask) + self.layers = {} + + def add_layer(self, name: str, data: np.ndarray): + """Adds a new environmental variable layer.""" + if data.shape != self.shape: + raise ValueError(f"Layer shape {data.shape} does not match world shape {self.shape}") + self.layers[name] = data + + def get_layer(self, name: str) -> np.ndarray: + return self.layers[name] + + @property + def cell_area(self) -> float: + """Returns the surface area of a single cell in square meters.""" + return float(self.cell_size ** 2) \ No newline at end of file diff --git a/src/gem/engine/processes.py b/src/gem/engine/processes.py new file mode 100644 index 0000000..6f912bf --- /dev/null +++ b/src/gem/engine/processes.py @@ -0,0 +1,152 @@ +from .environment_state import EnvironmentState +import numpy as np + +from .ecosystem_grid_state import EcosystemGridState +from ..vegetation import logistic_growth_delta + +# ============================================================================ +# PROCESS ARCHITECTURE: Pure Science + Thin Adapter +# ============================================================================ +# +# Every process has two parts: +# +# 1. PURE SCIENCE FUNCTION (e.g., src/gem/vegetation.py) +# - Input: explicit numpy arrays (biomass, parameters, environment data) +# - Output: delta array (change to be integrated) +# - No grid, no engine, no side effects +# - Works on any shape: (S,), (Y, S), or (X, Y, S) through broadcasting +# - Easy to test with hand-built data +# - Easy to reason about for ecology contributors +# - Example signature: +# +# def logistic_growth_delta(biomass, growth_rate, carrying_capacity, dt): +# """Logistic growth: dB/dt = r*B*(1 - B/K)""" +# assert biomass.shape == growth_rate.shape == carrying_capacity.shape +# return dt * growth_rate * biomass * (1.0 - biomass / carrying_capacity) +# +# 2. ADAPTER FUNCTION (e.g., apply_vegetation_growth below) +# - Input: grid and env (the engine's shared state) +# - Pulls the slices and parameters the science function needs +# - Calls the pure science function +# - Writes the delta back to the correct delta layer +# - Example pattern: +# +# def apply_vegetation_growth(grid, env): +# r = grid.registry.get_group_parameter("plants", "base_growth_rate") +# K = env.get_layer("carrying_capacity")[..., np.newaxis] +# B = grid.get_layer_view("biomass", "plants") +# +# delta = logistic_growth_delta(B, r, K, dt=1.0) +# +# with grid.edit_group_data("vegetation_delta", "plants") as d: +# d += delta +# +# WHY THIS PATTERN? +# - The science function signature documents what the process depends on +# (no hidden dependencies buried in grid/env reads) +# - Non-programmer ecologists can read, prototype, and test the science in a notebook +# without touching engine code +# - Adapters are boilerplate; the biology is in the science function +# - Delta layers make composition order-independent: every process sees frozen biomass_t +# +# ============================================================================ + + +# --- Mock Processes (What your colleagues will write) --- + +def apply_vegetation_growth(grid: EcosystemGridState, env: EnvironmentState): + """Adapter: vegetation growth process. + + Pulls data from grid/env, calls the pure science function, writes delta back. + The actual biology is in src.gem.vegetation.logistic_growth_delta(). + """ + # Extract parameters for the plants group + growth_rate = grid.registry.get_group_parameter("plants", "base_growth_rate") + carrying_capacity = env.get_layer("carrying_capacity")[..., np.newaxis] + biomass = grid.get_layer_view("biomass", "plants") + + # Broadcast growth_rate (scalar) and carrying_capacity to match biomass shape + growth_rate_broadcast = np.full_like(biomass, growth_rate) + + # Call pure science function + delta = logistic_growth_delta( + biomass=biomass, + growth_rate=growth_rate_broadcast, + carrying_capacity=carrying_capacity, + dt=1.0, # one time step + ) + + # Write delta back to the grid + with grid.edit_group_data("vegetation_delta", "plants") as plants_delta: + plants_delta += delta + + + +def apply_atn_step(grid: EcosystemGridState, env: EnvironmentState): + """ATN handles feeding, metabolism, and sets up growth rates for dispersal. + + Reads biomass and feeding links, writes trophic dynamics delta to biomass_delta. + """ + # The ATN process reads biomass to determine trophic interactions, + # computes feeding and metabolism fluxes, + # and accumulates the net trophic delta into biomass_delta. + # Also populates intermediate matrices like grid.net_growth_rate for other processes. + pass + +def apply_dispersal(grid: EcosystemGridState, env: EnvironmentState): + """Dispersal moves biomass between adjacent grid cells. + + Reads current biomass and net growth rates, writes dispersal delta to biomass_delta. + """ + # Dispersal reads grid.biomass (current state at start of step) + # and grid.net_growth_rate (computed by ATN), + # calculates diffusion and movement fluxes, + # and accumulates the dispersal delta into biomass_delta. + pass + + +# ============================================================================ +# CONTRIBUTING A NEW PROCESS +# ============================================================================ +# +# Step 1: Write the pure science function in a new module (e.g., src/gem/mortality.py): +# +# from numpy.typing import NDArray +# import numpy as np +# +# def mortality_delta( +# biomass: NDArray[np.float64], +# mortality_rate: NDArray[np.float64], +# dt: float, +# ) -> NDArray[np.float64]: +# """Loss due to background mortality. +# +# Args: +# biomass: Current biomass. Same shape for all inputs. +# mortality_rate: Per-capita mortality rate. Same shape as biomass. +# dt: Time step. +# +# Returns: +# Biomass delta (negative, since mortality removes biomass). +# """ +# assert biomass.shape == mortality_rate.shape +# return -dt * mortality_rate * biomass +# +# Step 2: Write tests in tests/test_mortality.py: +# - Use hand-built arrays, no engine +# - Test edge cases: zero mortality, all dead, etc. +# - Verify shape invariance: (S,), (Y, S), (X, Y, S) +# +# Step 3: Write the adapter here (e.g., apply_mortality below): +# - Register your delta layer (if new) +# - Extract parameters from grid/env +# - Call the science function +# - Write delta back +# +# Step 4: Register the adapter with the engine: +# model.add_process(apply_mortality) +# +# ============================================================================ + +# --- Placeholder adapters for ATN and dispersal --- +# (To be implemented by the respective teams; same pattern as apply_vegetation_growth) diff --git a/experiments/tuesday_architecture_proposal/engine_v2/SpeciesRegistry.py b/src/gem/engine/species_registry.py similarity index 100% rename from experiments/tuesday_architecture_proposal/engine_v2/SpeciesRegistry.py rename to src/gem/engine/species_registry.py diff --git a/src/gem/utils/data_fetching.py b/src/gem/utils/data_fetching.py new file mode 100644 index 0000000..47c7522 --- /dev/null +++ b/src/gem/utils/data_fetching.py @@ -0,0 +1,113 @@ +import os +import glob +import warnings +import numpy as np +import rasterio +from rasterio.transform import from_origin +from rasterio.warp import reproject, Resampling +from rasterio.crs import CRS + +warnings.filterwarnings('ignore') + +# โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€ +# MASTER GRID SPECIFICATIONS (North America Albers Equal Area Conic) +# โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€ +MASTER_PROJ = "+proj=aea +lat_0=40 +lon_0=-96 +lat_1=20 +lat_2=60 +x_0=0 +y_0=0 +datum=NAD83 +units=m +no_defs" +MASTER_CRS = CRS.from_proj4(MASTER_PROJ) +PIXEL_SIZE = 100_000 # 100km +X_MIN, X_MAX = -7_000_000.0, 5_000_000.0 +Y_MIN, Y_MAX = -2_000_000.0, 5_500_000.0 + +MASTER_HEIGHT = int((Y_MAX - Y_MIN) / PIXEL_SIZE) +MASTER_WIDTH = int((X_MAX - X_MIN) / PIXEL_SIZE) +MASTER_TRANSFORM = from_origin(X_MIN, Y_MAX, PIXEL_SIZE, PIXEL_SIZE) + +MASTER_PROFILE = { + 'driver': 'GTiff', + 'height': MASTER_HEIGHT, + 'width': MASTER_WIDTH, + 'count': 1, + 'dtype': 'float32', + 'crs': MASTER_CRS, + 'transform': MASTER_TRANSFORM, + 'nodata': np.nan, + 'compress': 'lzw' +} + +# โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€ +# DATA SOURCE DOCUMENTATION +# โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€โ”€ +# NASA NPP (MODIS): +# WorldClim Temperature: + +def check_required_files(file_patterns, source_name, download_url): + """Verifies existence of raw data files.""" + files = [] + for pattern in file_patterns: + found = glob.glob(pattern) + files.extend(found) + + if not files: + print(f"\n[MISSING DATA] {source_name} raw files not found.") + print(f"Please download from: {download_url}") + print(f"Search patterns attempted: {file_patterns}\n") + return None + return sorted(files) + +def reproject_to_master(src_path, dst_path, resampling=Resampling.bilinear): + """Reprojects an input raster to match the GEM Master Grid.""" + print(f" Processing {os.path.basename(src_path)} ...", end=' ', flush=True) + with rasterio.open(src_path) as src: + dst_array = np.full((1, MASTER_HEIGHT, MASTER_WIDTH), np.nan, dtype=np.float32) + + reproject( + source=rasterio.band(src, 1), + destination=dst_array, + src_transform=src.transform, + src_crs=src.crs, + src_nodata=src.nodata, + dst_transform=MASTER_TRANSFORM, + dst_crs=MASTER_CRS, + dst_nodata=np.nan, + resampling=resampling, + ) + + with rasterio.open(dst_path, 'w', **MASTER_PROFILE) as dst: + dst.write(dst_array[0], 1) + print("done.") + +def run_data_preparation(): + """Main entry point to prepare all environmental layers.""" + raw_dir = './data/raw' + proc_dir = './data/processed' + os.makedirs(proc_dir, exist_ok=True) + + print(f"--- Environmental Data Preparation ---") + print(f"Master Grid: {MASTER_HEIGHT} rows x {MASTER_WIDTH} cols ({PIXEL_SIZE/1000} km)\n") + + # 1. Process NASA NPP + npp_raw = check_required_files( + [os.path.join(raw_dir, 'MOD17A3H_Y_NPP_*.TIFF')], + "NASA MODIS NPP", + "https://neo.gsfc.nasa.gov/servlet/RenderData?si=2044462&cs=rgb&format=TIFF&width=720&height=360" + ) + if npp_raw: + dst = os.path.join(proc_dir, 'NPP_NorthAmerica_EqualArea_100km.tif') + reproject_to_master(npp_raw[0], dst) + + # 2. Process WorldClim Temperature + tavg_raw = check_required_files( + [os.path.join(raw_dir, 'wc2.1_10m_tavg', 'wc2.1_10m_tavg_*.tif')], + "WorldClim Tavg", + "https://geodata.ucdavis.edu/climate/worldclim/2_1/base/wc2.1_10m_tavg.zip" + ) + if tavg_raw: + for f in tavg_raw: + month = os.path.basename(f).split('_')[-1] + dst = os.path.join(proc_dir, f'Tavg_NorthAmerica_EqualArea_100km_month_{month}') + reproject_to_master(f, dst) + + print(f"\n--- Preparation Complete ---") + +if __name__ == "__main__": + run_data_preparation() \ No newline at end of file diff --git a/src/gem/vegetation.py b/src/gem/vegetation.py index 251cc43..ed9f42d 100644 --- a/src/gem/vegetation.py +++ b/src/gem/vegetation.py @@ -1,38 +1,80 @@ -"""Vegetation science functions. Pure numpy, no engine imports.""" +"""Vegetation science functions. Pure numpy, no engine imports. + +These functions implement the biological mechanisms of vegetation growth. +No grid plumbing, no side effects. Inputs and outputs are explicit arrays. + +Easy to test with hand-built data, easy to reason about for ecology contributors, +easy to prototype in a notebook with synthetic data. + +All functions work on any shape (S,), (Y, S), or (X, Y, S) through numpy broadcasting. +See docs/processes_implementation_specification.md for the shape contract. +""" import numpy as np from numpy.typing import NDArray + def logistic_growth_delta( biomass: NDArray[np.float64], growth_rate: NDArray[np.float64], carrying_capacity: NDArray[np.float64], dt: float, ) -> NDArray[np.float64]: - """Per-cell logistic growth applied to a biomass array. - - The function is shape-agnostic โ€” the same code runs on a 1-D species - vector or a full spatial grid. The team convention for which shapes - to actually use: - - - (S,) : one location, S species (smallest unit, handy for tests) - - (Y, S) : a row of Y locations, S species each - - (X, Y, S) : the full X-by-Y grid, S species each - - Convention on axes: the trailing axis is always species; leading axes - (if any) are spatial. - - All three array inputs must already have the **same shape**. The - adapter is responsible for reshaping per-species or per-cell - parameters (e.g. via ``np.broadcast_to``) before calling this - function. - - Returns the biomass delta over one time step ``dt``. + """Logistic (density-dependent) growth of vegetation biomass. + + Implements: dB/dt = rยทBยท(1 - B/K) + + where: + B = biomass per species (or per location and species) + r = intrinsic growth rate (species-dependent) + K = carrying capacity (environment and time-dependent) + + Typical ranges: + r โ‰ˆ 0.02โ€“0.3 per day (depends on growth strategy) + K โ‰ˆ 10^4โ€“10^5 kg/ha (depends on environment) + dt = 1 day (but parameterizable) + + Args: + biomass: Current vegetation biomass. Shape: (S,), (Y, S), or (X, Y, S). + growth_rate: Intrinsic growth rate r. Must broadcast with biomass. + carrying_capacity: Environmental capacity K. Must broadcast with biomass. + dt: Time step in days (default 1.0). + + Returns: + Biomass change dB over dt. Shape broadcasts to match all inputs. + + Shape contract (enforced by assertion): + All three array inputs must be broadcastable together. + Examples: + - (X, Y, 15) with (X, Y, 15) with (X, Y, 1) โœ“ broadcasts to (X, Y, 15) + - (X, Y, 15) with (15,) with (X, Y, 1) โœ“ broadcasts to (X, Y, 15) + - (X, Y, 15) with (X, Y, 10) with (X, Y, 1) โœ— 15 != 10, cannot broadcast + + Example: + >>> B = np.array([100.0, 200.0, 50.0]) # 3 species + >>> r = np.array([0.1, 0.15, 0.2]) # per-species rates + >>> K = np.array([500.0, 400.0, 300.0]) # per-species carrying capacity + >>> dt = 1.0 + >>> dB = logistic_growth_delta(B, r, K, dt) + >>> dB.shape + (3,) + >>> dB[0] # positive because B < K for species 0 + array(7.2) """ - # All inputs must already share the same shape โ€” catches broadcast - # mistakes in the adapter early, with a clear failure point. - assert biomass.shape == growth_rate.shape == carrying_capacity.shape + # Shape contract: all inputs must broadcast together. + # This allows (X, Y, 15) + (X, Y, 1) + (15,) etc., but catches mismatches. + try: + biomass_bc, rate_bc, capacity_bc = np.broadcast_arrays( + biomass, growth_rate, carrying_capacity + ) + except ValueError as e: + raise ValueError( + f"Inputs do not broadcast: biomass {biomass.shape}, " + f"growth_rate {growth_rate.shape}, carrying_capacity {carrying_capacity.shape}. " + f"Details: {e}" + ) + # Logistic equation: dB/dt = r * B * (1 - B/K) delta_biomass = dt * growth_rate * biomass * (1.0 - biomass / carrying_capacity) return delta_biomass diff --git a/tests/test_vegetation.py b/tests/test_vegetation.py index 97dd13e..c1bc1b3 100644 --- a/tests/test_vegetation.py +++ b/tests/test_vegetation.py @@ -1,21 +1,30 @@ -import numpy as np +"""Unit tests for vegetation science functions (no engine, no grid plumbing).""" +import numpy as np from gem.vegetation import logistic_growth_delta -def test_logistic_growth_zero_at_carrying_capacity(): - B = np.array([100.0, 100.0, 100.0]) - r = np.array([0.1, 0.2, 0.3]) - K = np.array([100.0, 100.0, 100.0]) +def test_logistic_growth_simple(): + """Biomass should grow when below carrying capacity, shrink when above.""" + B = np.array([50.0, 150.0]) # below and above K + r = np.array([0.1, 0.1]) + K = np.array([100.0, 100.0]) + delta = logistic_growth_delta(B, r, K, dt=1.0) - np.testing.assert_allclose(delta, 0.0) + + assert delta[0] > 0 # Growing (B < K) + assert delta[1] < 0 # Shrinking (B > K) -def test_logistic_growth_runs_on_grid(): - shape = (4, 5, 3) - B = np.full(shape, 50.0) - r = np.full(shape, 0.1) - K = np.full(shape, 100.0) +def test_logistic_growth_on_full_grid(): + """End-to-end: broadcast spatial grid with per-species parameters.""" + # Full spatial grid: (X, Y, S) + B = np.full((180, 360, 15), 50000.0) + r = np.full((180, 360, 15), 0.1) + K = np.full((180, 360, 1), 100000.0) # broadcast across species + delta = logistic_growth_delta(B, r, K, dt=1.0) - assert delta.shape == shape - assert np.all(delta > 0) + + assert delta.shape == (180, 360, 15) + assert np.all(delta > 0) # All species growing (B < K) +