diff --git a/docs/_config.yml b/docs/_config.yml index 1cb5dc24..c0711748 100644 --- a/docs/_config.yml +++ b/docs/_config.yml @@ -12,6 +12,7 @@ execute: execute_notebooks: cache exclude_patterns: - clay-v0-*.ipynb + - partial-inputs.ipynb # Define the name of the latex output file for PDF builds latex: diff --git a/docs/_toc.yml b/docs/_toc.yml index 299721a6..e05334c7 100644 --- a/docs/_toc.yml +++ b/docs/_toc.yml @@ -40,6 +40,8 @@ parts: file: clay-v0-location-embeddings - title: Interpolating images in embedding space file: clay-v0-interpolation + - title: Generate embeddings from partial inputs + file: partial-inputs - caption: About Clay chapters: - title: GitHub diff --git a/docs/partial-inputs.ipynb b/docs/partial-inputs.ipynb new file mode 100644 index 00000000..54cf3bc3 --- /dev/null +++ b/docs/partial-inputs.ipynb @@ -0,0 +1,451 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "id": "8dd35554-a9dd-49cf-b9fa-24fa8ae6cecf", + "metadata": {}, + "source": [ + "# Burn scar analysis using embeddings from partial inputs\n", + "This notebook contains a complete example for how to run Clay. It\n", + "combines the following three different aspects\n", + "\n", + "1. Create single-chip datacubes with time series data for a location and a date range\n", + "2. Run the model with partial inputs, in this case RGB + NIR\n", + "3. Study burn scares through the embeddings generated for that datacube\n", + "\n", + "## Let's start with importing and creating constants" + ] + }, + { + "cell_type": "code", + "execution_count": 1, + "id": "b7bcff1e-bdb5-47f8-aa0e-d68d6fdd3476", + "metadata": {}, + "outputs": [], + "source": [ + "# Ensure working directory is the repo home\n", + "import os\n", + "\n", + "os.chdir(\"..\")" + ] + }, + { + "cell_type": "code", + "execution_count": 2, + "id": "15d65ec9-86aa-4275-89ba-ec79fdbad361", + "metadata": {}, + "outputs": [], + "source": [ + "import warnings\n", + "from pathlib import Path\n", + "\n", + "import geopandas as gpd\n", + "import matplotlib.pyplot as plt\n", + "import numpy\n", + "import pandas as pd\n", + "import pystac_client\n", + "import rasterio\n", + "import rioxarray # noqa: F401\n", + "import stackstac\n", + "import torch\n", + "from rasterio.enums import Resampling\n", + "from shapely import Point\n", + "from sklearn import decomposition\n", + "\n", + "from src.datamodule import ClayDataModule\n", + "from src.model_clay import CLAYModule\n", + "\n", + "warnings.filterwarnings(\"ignore\")\n", + "\n", + "BAND_GROUPS = {\n", + " \"rgb\": [\"red\", \"green\", \"blue\"],\n", + " \"rededge\": [\"rededge1\", \"rededge2\", \"rededge3\", \"nir08\"],\n", + " \"nir\": [\n", + " \"nir\",\n", + " ],\n", + " \"swir\": [\"swir16\", \"swir22\"],\n", + " \"sar\": [\"vv\", \"vh\"],\n", + "}\n", + "\n", + "STAC_API = \"https://earth-search.aws.element84.com/v1\"\n", + "COLLECTION = \"sentinel-2-l2a\"" + ] + }, + { + "cell_type": "markdown", + "id": "a6341305-9c44-4a1e-847c-80d77b01c0bf", + "metadata": {}, + "source": [ + "## Search for imagery over an area of interest\n", + "In this example we use a location and date range to visualize a forest fire that happened in [Monchique in 2018](https://pt.wikipedia.org/wiki/Inc%C3%AAndio_de_Monchique_de_2018)" + ] + }, + { + "cell_type": "code", + "execution_count": 3, + "id": "a1886f5a-8669-40e7-8fae-e45619570e3c", + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Found 12 items\n" + ] + } + ], + "source": [ + "# Point over Monchique Portugal\n", + "poi = 37.30939, -8.57207\n", + "\n", + "# Dates of a large forest fire\n", + "start = \"2018-07-01\"\n", + "end = \"2018-09-01\"\n", + "\n", + "catalog = pystac_client.Client.open(STAC_API)\n", + "\n", + "search = catalog.search(\n", + " collections=[COLLECTION],\n", + " datetime=f\"{start}/{end}\",\n", + " bbox=(poi[1] - 1e-5, poi[0] - 1e-5, poi[1] + 1e-5, poi[0] + 1e-5),\n", + " max_items=100,\n", + " query={\"eo:cloud_cover\": {\"lt\": 80}},\n", + ")\n", + "\n", + "items = search.get_all_items()\n", + "\n", + "print(f\"Found {len(items)} items\")" + ] + }, + { + "cell_type": "markdown", + "id": "c4ba5c36-90a6-427c-80c5-2a83ad11a1b0", + "metadata": {}, + "source": [ + "## Download the data\n", + "Get the data into a numpy array and visualize the imagery. The burn scar is visible in the last five images." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "c371501c-3ef0-4507-9073-0521a1c733be", + "metadata": {}, + "outputs": [], + "source": [ + "# Extract coordinate system from first item\n", + "epsg = items[0].properties[\"proj:epsg\"]\n", + "\n", + "# Convert point into the image projection\n", + "poidf = gpd.GeoDataFrame(\n", + " pd.DataFrame(),\n", + " crs=\"EPSG:4326\",\n", + " geometry=[Point(poi[1], poi[0])],\n", + ").to_crs(epsg)\n", + "\n", + "coords = poidf.iloc[0].geometry.coords[0]\n", + "\n", + "# Create bounds of the correct size, the model\n", + "# requires 512x512 pixels at 10m resolution.\n", + "bounds = (\n", + " coords[0] - 2560,\n", + " coords[1] - 2560,\n", + " coords[0] + 2560,\n", + " coords[1] + 2560,\n", + ")\n", + "\n", + "# Retrieve the pixel values, for the bounding box in\n", + "# the target projection. In this example we use only\n", + "# the RGB and NIR band groups.\n", + "stack = stackstac.stack(\n", + " items,\n", + " bounds=bounds,\n", + " snap_bounds=False,\n", + " epsg=epsg,\n", + " resolution=10,\n", + " dtype=\"float32\",\n", + " rescale=False,\n", + " fill_value=0,\n", + " assets=BAND_GROUPS[\"rgb\"] + BAND_GROUPS[\"nir\"],\n", + " resampling=Resampling.nearest,\n", + ")\n", + "\n", + "stack = stack.compute()\n", + "\n", + "stack.sel(band=[\"red\", \"green\", \"blue\"]).plot.imshow(\n", + " row=\"time\", rgb=\"band\", vmin=0, vmax=2000, col_wrap=6\n", + ")" + ] + }, + { + "cell_type": "markdown", + "id": "ce633fb1-fc82-4c88-8204-cda47aa9c874", + "metadata": {}, + "source": [ + "![Minicube visualization](https://github.com/Clay-foundation/model/assets/901647/c6e924e5-6ba1-4924-b99a-df8b90731a5f)" + ] + }, + { + "cell_type": "markdown", + "id": "77e7c22c-1bfd-4281-bb12-8330c3eedc25", + "metadata": {}, + "source": [ + "## Write data to tif files\n", + "To use the mini datacube in the Clay dataloader, we need to write the\n", + "images to tif files on disk. These tif files are then used by the Clay\n", + "data loader for creating embeddings below." + ] + }, + { + "cell_type": "code", + "execution_count": 5, + "id": "6509c3b2-a67c-447d-a7a1-e5fbcc1e35b5", + "metadata": {}, + "outputs": [], + "source": [ + "outdir = Path(\"data/minicubes\")\n", + "assert outdir.exists()\n", + "\n", + "# Write tile to output dir\n", + "for tile in stack:\n", + " # Grid code like MGRS-29SNB\n", + " mgrs = str(tile.coords[\"grid:code\"].values).split(\"-\")[1]\n", + " date = str(tile.time.values)[:10]\n", + "\n", + " name = \"{dir}/claytile_{mgrs}_{date}.tif\".format(\n", + " dir=outdir,\n", + " mgrs=mgrs,\n", + " date=date.replace(\"-\", \"\"),\n", + " )\n", + " tile.rio.to_raster(name, compress=\"deflate\")\n", + "\n", + " with rasterio.open(name, \"r+\") as rst:\n", + " rst.update_tags(date=date)" + ] + }, + { + "cell_type": "markdown", + "id": "ebc4b6ee-db58-4005-9689-a7d0acdc6a79", + "metadata": { + "scrolled": true + }, + "source": [ + "## Create embeddings\n", + "Now switch gears and load the tiles to create embeddings and analyze them. \n", + "\n", + "The model checkpoint can be loaded directly from huggingface, and the data\n", + "directory points to the directory we created in the steps above.\n", + "\n", + "Note that the normalization parameters for the data module need to be \n", + "adapted based on the band groups that were selected as partial input. The\n", + "full set of normalization parameters can be found [here](https://github.com/Clay-foundation/model/blob/main/src/datamodule.py#L108)." + ] + }, + { + "cell_type": "markdown", + "id": "d89e0135-9473-4f76-9f09-e4e295dd51c9", + "metadata": {}, + "source": [ + "### Load the model and set up the data module" + ] + }, + { + "cell_type": "code", + "execution_count": 6, + "id": "301ee2db-c5fc-4628-b837-12e6ea477415", + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Total number of chips: 12\n" + ] + } + ], + "source": [ + "DATA_DIR = \"data/minicubes\"\n", + "CKPT_PATH = \"https://huggingface.co/made-with-clay/Clay1/resolve/main/Clay_v0.1_epoch-24_val-loss-0.46.ckpt\"\n", + "\n", + "# Load model\n", + "rgb_model = CLAYModule.load_from_checkpoint(\n", + " CKPT_PATH,\n", + " mask_ratio=0.0,\n", + " band_groups={\"rgb\": (2, 1, 0), \"nir\": (3,)},\n", + " bands=4,\n", + " strict=False, # ignore the extra parameters in the checkpoint\n", + ")\n", + "# Set the model to evaluation mode\n", + "rgb_model.eval()\n", + "\n", + "\n", + "# Load the datamodule, with the reduced set of\n", + "class ClayDataModuleRGB(ClayDataModule):\n", + " MEAN = [\n", + " 1369.03, # red\n", + " 1597.68, # green\n", + " 1741.10, # blue\n", + " 2858.43, # nir\n", + " ]\n", + " STD = [\n", + " 2026.96, # red\n", + " 2011.88, # green\n", + " 2146.35, # blue\n", + " 2016.38, # nir\n", + " ]\n", + "\n", + "\n", + "data_dir = Path(DATA_DIR)\n", + "\n", + "dm = ClayDataModuleRGB(data_dir=str(data_dir.absolute()), batch_size=20)\n", + "dm.setup(stage=\"predict\")\n", + "trn_dl = iter(dm.predict_dataloader())" + ] + }, + { + "cell_type": "markdown", + "id": "db3f3e5e-8668-4830-9c77-cc1d8cb35234", + "metadata": {}, + "source": [ + "### Create the embeddings for the images over the forest fire\n", + "This will loop through the images returned by the data loader\n", + "and evaluate the model for each one of the images. The raw\n", + "embeddings are reduced to mean values to simplify the data." + ] + }, + { + "cell_type": "code", + "execution_count": 7, + "id": "c5762240-9d22-4ebd-8e39-83fc6594a459", + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Average embeddings have shape (12, 768)\n" + ] + } + ], + "source": [ + "embeddings = []\n", + "\n", + "for batch in trn_dl:\n", + " with torch.inference_mode():\n", + " # Move data from to the device of model\n", + " batch[\"pixels\"] = batch[\"pixels\"].to(rgb_model.device)\n", + " # Pass just the specific band through the model\n", + " batch[\"timestep\"] = batch[\"timestep\"].to(rgb_model.device)\n", + " batch[\"latlon\"] = batch[\"latlon\"].to(rgb_model.device)\n", + "\n", + " # Pass pixels, latlon, timestep through the encoder to create encoded patches\n", + " (\n", + " unmasked_patches,\n", + " unmasked_indices,\n", + " masked_indices,\n", + " masked_matrix,\n", + " ) = rgb_model.model.encoder(batch)\n", + "\n", + " embeddings.append(unmasked_patches.detach().cpu().numpy())\n", + "\n", + "embeddings = numpy.vstack(embeddings)\n", + "\n", + "embeddings_mean = embeddings[:, :-2, :].mean(axis=1)\n", + "\n", + "print(f\"Average embeddings have shape {embeddings_mean.shape}\")" + ] + }, + { + "cell_type": "markdown", + "id": "72db5745-21c6-4b8e-b8f7-cb48c0f9c9ef", + "metadata": {}, + "source": [ + "## Analyze embeddings\n", + "Now we can make a simple analysis of the embeddings. We reduce all the\n", + "embeddings to a single number using Principle Component Analysis. Then\n", + "we can plot the principal components. The effect of the fire on the\n", + "embeddings is clearly visible. We use the following color code in the graph:\n", + "\n", + "| Color | Interpretation |\n", + "|---|---|\n", + "| Green | Cloudy Images |\n", + "| Blue | Before the fire |\n", + "| Red | After the fire |" + ] + }, + { + "cell_type": "code", + "execution_count": 8, + "id": "88f3b2dc-8f2a-447b-a6af-b04e0d1ff61c", + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "" + ] + }, + "execution_count": 8, + "metadata": {}, + "output_type": "execute_result" + }, + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAjgAAAHDCAYAAAApyGCxAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjguMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8g+/7EAAAACXBIWXMAAA9hAAAPYQGoP6dpAAA7oUlEQVR4nO3deXgUZbr+8bsJJAE0jYqEAJHFg+wjEIQEjYJKZFVEhuCSGRVxOI6jyM8NmXFcJy5HB1xwGR05OgyikqgzLIoKgkMQgQR3REDDkoAgdANKCJ3n90dOWpokkEA63V39/VxXXdpvv115n4Ki71S9VeUyMxMAAICDNAj1AAAAAOoaAQcAADgOAQcAADgOAQcAADgOAQcAADgOAQcAADgOAQcAADhOw1APIBTKysq0detWnXjiiXK5XKEeDgAAqAEz0549e9SqVSs1aHDkYzRRGXC2bt2q5OTkUA8DAAAcg02bNqlNmzZH7BOVAefEE0+UVL6BEhISQjwaAABQE16vV8nJyf7v8SOJyoBTcVoqISGBgAMAQISpyfQSJhkDAADHIeAAAADHIeAAAADHIeAAAADHIeAAAADHIeAAAADHIeAAAADHIeAAAADHIeAAAADHico7GSO8HSj1afrcpVq/rUinJybphmHpim0UE+phAQAiCAEHYeX2l3L0+Jc3y3fC5vKGYunW/7TRpK7T9Mg1o0I7OABAxOAUFcLG7S/l6NHvR8vXdHNAu6/pFj36/Wjd/lJOiEYGAIg0QQ04S5Ys0YgRI9SqVSu5XC69+eabR+yfk5OjQYMG6dRTT1VCQoLS0tL0zjvvBPSZMWOGXC5XpWX//v1BrATBdqDUp8e/vFmSSYc/Q81lkqTHv5yoA6W+eh8bACDyBDXg7Nu3T2eeeaaeeuqpGvVfsmSJBg0apHnz5mnVqlUaOHCgRowYofz8/IB+CQkJKioqClji4+ODUQLqyfS5S8tPS1X3gFiXyXfCJk2fu7RexwUAiExBnYMzZMgQDRkypMb9p06dGvD6L3/5i9566y3961//Uq9evfztLpdLLVu2rKthIgys31ZUp/0AANEtrCcZl5WVac+ePTr55JMD2vfu3au2bdvK5/OpZ8+euv/++wMC0OFKSkpUUlLif+31eoMyXq7+OXanJyZJxTXsBwDAUYT1JOPHHntM+/bt05gxY/xtnTt31owZM/T2229r1qxZio+P19lnn61169ZVu57s7Gy53W7/kpycXOdjvf2lHDW5q51uWTNQTxVfoVvWDFSTu9oxMbaGbhiWrpi9bSSr5hyVuRSzN1k3DEuv34EBACKSy8ysXn6Qy6Xc3FyNHDmyRv1nzZql6667Tm+99ZYuvPDCavuVlZWpd+/eOvfcc/XEE09U2aeqIzjJycnyeDxKSEioVR1Vqbj6p9IE2f/7sr6t7Rtc4lwDv2xH+ScWS2I7AgAklX9/u93uGn1/h+URnNmzZ2vcuHF67bXXjhhuJKlBgwY666yzjngEJy4uTgkJCQFLXeHqn7rzyDWjdFvbNxSzr3VAe8y+NoQbAECthN0cnFmzZunaa6/VrFmzNGzYsKP2NzMVFBSoR48e9TC6yvxX/1TnkKt/Jo4cUG/jilSPXDNKD5RewlwmAMBxCWrA2bt3r7799lv/640bN6qgoEAnn3yyTjvtNE2ePFlbtmzRyy+/LKk83PzmN7/RtGnTlJqaquLi8lmnjRs3ltvtliTde++9Sk1NVceOHeX1evXEE0+ooKBATz/9dDBLqRZX/9S92EYxhEEAwHEJ6imqlStXqlevXv4rnCZNmqRevXrp7rvvliQVFRWpsLDQ3/+5557TwYMH9fvf/15JSUn+5eabb/b32b17t66//np16dJFGRkZ2rJli5YsWaK+ffsGs5Rq1fSqHq7+AQCg/tTbJONwUptJSkdzoNSnJne1k6/plsCJsRXMpZh9bfTTXzZymgUAgOMQ8ZOMI0lsoxhN6jqt/MXhlzj/3+tJXacSbgAAqEcEnDrA1T8AAIQXTlHV8SXjXP0DAEBw1Ob7O+wuE49kXP0DAEB44BQVAABwHAIOAABwHAIOAABwHAIOAABwHAIOAABwHAIOAABwHAIOAABwHAIOAABwHAIOAABwHAIOAABwHAIOAABwHAIOAABwHAIOAABwHAIOAABwHAIOAABwHAIOAABwHAIOAABwHAIOAABwHAIOAABwHAIOAABwHAIOAABwHAIOAABwHAIOAABwHAIOAABwHAIOAABwHAIOAABwHAIOAABwHAIOAABwHAIOAABwHAIOAABwHAIOAABwHAIOAABwHAIOAABwnKAGnCVLlmjEiBFq1aqVXC6X3nzzzaN+5sMPP1RKSori4+PVoUMHPfvss5X6zJkzR127dlVcXJy6du2q3NzcIIweAABEqqAGnH379unMM8/UU089VaP+Gzdu1NChQ5Wenq78/HzddddduummmzRnzhx/n7y8PGVmZiorK0tr1qxRVlaWxowZo48//jhYZQAAgAjjMjOrlx/kcik3N1cjR46sts8dd9yht99+W1999ZW/bcKECVqzZo3y8vIkSZmZmfJ6vZo/f76/z+DBg3XSSSdp1qxZNRqL1+uV2+2Wx+NRQkLCsRUEAADqVW2+v8NqDk5eXp4yMjIC2i666CKtXLlSpaWlR+yzbNmyatdbUlIir9cbsAAAAOcKq4BTXFysxMTEgLbExEQdPHhQO3bsOGKf4uLiatebnZ0tt9vtX5KTk+t+8AAAIGyEVcCRyk9lHariDNqh7VX1ObztUJMnT5bH4/EvmzZtqsMRAwCAcNMw1AM4VMuWLSsdidm+fbsaNmyoU0455Yh9Dj+qc6i4uDjFxcXV/YABAEBYCqsjOGlpaVq4cGFA27vvvqs+ffqoUaNGR+zTv3//ehsnAAAIb0E9grN37159++23/tcbN25UQUGBTj75ZJ122mmaPHmytmzZopdffllS+RVTTz31lCZNmqTx48crLy9PL774YsDVUTfffLPOPfdcPfzww7rkkkv01ltv6b333tNHH30UzFIAAEAECeoRnJUrV6pXr17q1auXJGnSpEnq1auX7r77bklSUVGRCgsL/f3bt2+vefPmafHixerZs6fuv/9+PfHEE7rsssv8ffr3769XX31VL730kn71q19pxowZmj17tvr16xfMUgAAQASpt/vghBPugwMAQOSJ2PvgAAAA1AUCDgAAcBwCDgAAcBwCDgAAcBwCDgAAcBwCDgAAcBwCDgAAcBwCDgAAcBwCDgAAcBwCDgAAcBwCDgAAcBwCDgAAcBwCDgAAcBwCDgAAcBwCDgAAcBwCDgAAcBwCDgAAcBwCDgAAcBwCDgAAcBwCDgAAcBwCDgAAcBwCDgAAcBwCDgAAcBwCDgAAcBwCDgAAcBwCDgAAcBwCDgAAcBwCDgAAcBwCDgAAcBwCDgAAcBwCDgAAcBwCDgAAcBwCDgAAcBwCDgAAcBwCDgAAcBwCDgAAcBwCDgAAcJx6CTjTp09X+/btFR8fr5SUFC1durTavldffbVcLlelpVu3bv4+M2bMqLLP/v3766McAAAQ5oIecGbPnq2JEydqypQpys/PV3p6uoYMGaLCwsIq+0+bNk1FRUX+ZdOmTTr55JP161//OqBfQkJCQL+ioiLFx8cHuxwAABABGgb7Bzz++OMaN26crrvuOknS1KlT9c477+iZZ55RdnZ2pf5ut1tut9v/+s0339SuXbt0zTXXBPRzuVxq2bJljcZQUlKikpIS/2uv13sspQAAgAgR1CM4Bw4c0KpVq5SRkRHQnpGRoWXLltVoHS+++KIuvPBCtW3bNqB97969atu2rdq0aaPhw4crPz+/2nVkZ2f7g5Pb7VZycnLtiwEAABEjqAFnx44d8vl8SkxMDGhPTExUcXHxUT9fVFSk+fPn+4/+VOjcubNmzJiht99+W7NmzVJ8fLzOPvtsrVu3rsr1TJ48WR6Px79s2rTp2IsCAABhL+inqKTy00mHMrNKbVWZMWOGmjVrppEjRwa0p6amKjU11f/67LPPVu/evfXkk0/qiSeeqLSeuLg4xcXFHdvgAQBAxAnqEZzmzZsrJiam0tGa7du3Vzqqczgz09///ndlZWUpNjb2iH0bNGigs846q9ojOAAAILoENeDExsYqJSVFCxcuDGhfuHCh+vfvf8TPfvjhh/r22281bty4o/4cM1NBQYGSkpKOa7wAAMAZgn6KatKkScrKylKfPn2Ulpam559/XoWFhZowYYKk8vkxW7Zs0csvvxzwuRdffFH9+vVT9+7dK63z3nvvVWpqqjp27Civ16snnnhCBQUFevrpp4NdDgAAiABBDziZmZnauXOn7rvvPhUVFal79+6aN2+e/6qooqKiSvfE8Xg8mjNnjqZNm1blOnfv3q3rr79excXFcrvd6tWrl5YsWaK+ffsGuxwAABABXGZmoR5EffN6vXK73fJ4PEpISAj1cAAAQA3U5vubZ1EBAADHIeAAAADHIeAAAADHIeAAAADHIeAAAADHIeAAAADHIeAAAADHIeAAAADHIeAAAADHIeAAAADHIeAAAADHIeAAAADHIeAAAADHIeAAAADHIeAAAADHIeAAAADHIeAAAADHIeAAAADHIeAAAADHIeAAAADHIeAAAADHIeAAAADHIeAAAADHIeAAAADHIeAAAADHIeAAAADHIeAAAADHIeAAAADHIeAAAADHIeAAAADHIeAAAADHIeAAAADHIeAAAADHIeAAAADHIeAAAADHIeAAAADHqZeAM336dLVv317x8fFKSUnR0qVLq+27ePFiuVyuSsvXX38d0G/OnDnq2rWr4uLi1LVrV+Xm5ga7DAAAECGCHnBmz56tiRMnasqUKcrPz1d6erqGDBmiwsLCI35u7dq1Kioq8i8dO3b0v5eXl6fMzExlZWVpzZo1ysrK0pgxY/Txxx8HuxwAABABXGZmwfwB/fr1U+/evfXMM8/427p06aKRI0cqOzu7Uv/Fixdr4MCB2rVrl5o1a1blOjMzM+X1ejV//nx/2+DBg3XSSSdp1qxZRx2T1+uV2+2Wx+NRQkJC7YsCAAD1rjbf30E9gnPgwAGtWrVKGRkZAe0ZGRlatmzZET/bq1cvJSUl6YILLtCiRYsC3svLy6u0zosuuqjadZaUlMjr9QYsAADAuYIacHbs2CGfz6fExMSA9sTERBUXF1f5maSkJD3//POaM2eOcnJy1KlTJ11wwQVasmSJv09xcXGt1pmdnS232+1fkpOTj7MyAAAQzhrWxw9xuVwBr82sUluFTp06qVOnTv7XaWlp2rRpk/7nf/5H55577jGtc/LkyZo0aZL/tdfrJeQAAOBgQT2C07x5c8XExFQ6srJ9+/ZKR2COJDU1VevWrfO/btmyZa3WGRcXp4SEhIAFAAA4V1ADTmxsrFJSUrRw4cKA9oULF6p///41Xk9+fr6SkpL8r9PS0iqt8913363VOgEAgHMF/RTVpEmTlJWVpT59+igtLU3PP/+8CgsLNWHCBEnlp4+2bNmil19+WZI0depUtWvXTt26ddOBAwf0j3/8Q3PmzNGcOXP867z55pt17rnn6uGHH9Yll1yit956S++9954++uijYJcDAAAiQNADTmZmpnbu3Kn77rtPRUVF6t69u+bNm6e2bdtKkoqKigLuiXPgwAHdeuut2rJlixo3bqxu3bpp7ty5Gjp0qL9P//799eqrr+qPf/yj/vSnP+n000/X7Nmz1a9fv2CXAwAAIkDQ74MTjrgPDgAAkSds7oMDAAAQCgQcAADgOAQcAADgOAQcAADgOAQcAADgOAQcAADgOAQcAADgOAQcAADgOAQcAADgOAQcAADgOAQcAADgOAQcAADgOAQcAADgOAQcAADgOAQcAADgOAQcAADgOAQcAADgOAQcAADgOAQcAADgOAQcAADgOAQcAADgOAQcAADgOAQcAADgOAQcAADgOAQcAADgOAQcAADgOAQcAADgOAQcAADgOAQcAADgOAQcAADgOAQcAADgOAQcAADgOAQcAADgOAQcAADgOAQcAADgOAQcAADgOAQcAADgOPUScKZPn6727dsrPj5eKSkpWrp0abV9c3JyNGjQIJ166qlKSEhQWlqa3nnnnYA+M2bMkMvlqrTs378/2KUAAIAIEPSAM3v2bE2cOFFTpkxRfn6+0tPTNWTIEBUWFlbZf8mSJRo0aJDmzZunVatWaeDAgRoxYoTy8/MD+iUkJKioqChgiY+PD3Y5AAAgArjMzIL5A/r166fevXvrmWee8bd16dJFI0eOVHZ2do3W0a1bN2VmZuruu++WVH4EZ+LEidq9e/cxjcnr9crtdsvj8SghIeGY1gEAAOpXbb6/g3oE58CBA1q1apUyMjIC2jMyMrRs2bIaraOsrEx79uzRySefHNC+d+9etW3bVm3atNHw4cMrHeE5VElJibxeb8ACAACcK6gBZ8eOHfL5fEpMTAxoT0xMVHFxcY3W8dhjj2nfvn0aM2aMv61z586aMWOG3n77bc2aNUvx8fE6++yztW7duirXkZ2dLbfb7V+Sk5OPvSgAABD26mWSscvlCnhtZpXaqjJr1izdc889mj17tlq0aOFvT01N1VVXXaUzzzxT6enpeu2113TGGWfoySefrHI9kydPlsfj8S+bNm06voIAAEBYaxjMlTdv3lwxMTGVjtZs37690lGdw82ePVvjxo3T66+/rgsvvPCIfRs0aKCzzjqr2iM4cXFxiouLq93gAQBAxArqEZzY2FilpKRo4cKFAe0LFy5U//79q/3crFmzdPXVV+uf//ynhg0bdtSfY2YqKChQUlLScY8ZAABEvqAewZGkSZMmKSsrS3369FFaWpqef/55FRYWasKECZLKTx9t2bJFL7/8sqTycPOb3/xG06ZNU2pqqv/oT+PGjeV2uyVJ9957r1JTU9WxY0d5vV498cQTKigo0NNPPx3scgAAQAQIesDJzMzUzp07dd9996moqEjdu3fXvHnz1LZtW0lSUVFRwD1xnnvuOR08eFC///3v9fvf/97f/tvf/lYzZsyQJO3evVvXX3+9iouL5Xa71atXLy1ZskR9+/YNdjkAACACBP0+OOGI++AAABB5wuY+OAAAAKFAwAEAAI5DwAEAAI5DwAEAAI5DwAEAAI5DwAEAAI5DwAEAAI5DwAEAAI5DwAEAAI5DwAEAAI5DwAEAAI5DwAEAAI5DwAEAAI5DwAEAAI5DwAEAAI5DwAEAAI5DwAEAAI5DwAEAAI5DwAEAAI5DwAEAAI5DwAEAAI5DwAEAAI5DwAEAAI5DwAEAAI5DwAEAAI5DwAEAAI5DwAEAAI5DwAEAAI5DwAEAAI5DwAEAAI5DwAEAAI5DwAEAAI5DwAEAAI5DwAEAAI5DwAEAAI5DwAEAAI5DwAEAAI5TLwFn+vTpat++veLj45WSkqKlS5cesf+HH36olJQUxcfHq0OHDnr22Wcr9ZkzZ466du2quLg4de3aVbm5ucEavqP4fNLixdKsWeX/9flCPSIAAOpe0APO7NmzNXHiRE2ZMkX5+flKT0/XkCFDVFhYWGX/jRs3aujQoUpPT1d+fr7uuusu3XTTTZozZ46/T15enjIzM5WVlaU1a9YoKytLY8aM0ccffxzsciJaTo7Urp00cKB0xRXl/23XrrwdAAAncZmZBfMH9OvXT71799Yzzzzjb+vSpYtGjhyp7OzsSv3vuOMOvf322/rqq6/8bRMmTNCaNWuUl5cnScrMzJTX69X8+fP9fQYPHqyTTjpJs2bNqrTOkpISlZSU+F97vV4lJyfL4/EoISGhTuoMdzk50ujR0uF/2i5X+X/feEMaNar+xwUAQE15vV653e4afX8H9QjOgQMHtGrVKmVkZAS0Z2RkaNmyZVV+Ji8vr1L/iy66SCtXrlRpaekR+1S3zuzsbLndbv+SnJx8rCVFJJ9PuvnmyuFG+qVt4kROVwEAnCOoAWfHjh3y+XxKTEwMaE9MTFRxcXGVnykuLq6y/8GDB7Vjx44j9qlunZMnT5bH4/EvmzZtOtaSItLSpdLmzdW/byZt2lTeDwAAJ2hYHz/EVXEe5P+YWaW2o/U/vL0264yLi1NcXFytxuwkRUV12w8AgHAX1IDTvHlzxcTEVDqysn379kpHYCq0bNmyyv4NGzbUKaeccsQ+1a0z2iUl1W0/lJ/OW7q0PBQmJUnp6VJMTKhHFXnYjgCCJainqGJjY5WSkqKFCxcGtC9cuFD9+/ev8jNpaWmV+r/77rvq06ePGjVqdMQ+1a0z2qWnS23a/DKh+HAul5ScXN4PR8fVaHWD7QggqCzIXn31VWvUqJG9+OKL9uWXX9rEiROtadOm9t1335mZ2Z133mlZWVn+/hs2bLAmTZrYLbfcYl9++aW9+OKL1qhRI3vjjTf8ff7zn/9YTEyMPfTQQ/bVV1/ZQw89ZA0bNrTly5fXaEwej8ckmcfjqdtiw9icOWYuV/lSPuumfKlomzMn1COMDBXb8dBtGK7b8eBBs0WLzP75z/L/HjwY6hH9IpK2I4DwUZvv76AHHDOzp59+2tq2bWuxsbHWu3dv+/DDD/3v/fa3v7XzzjsvoP/ixYutV69eFhsba+3atbNnnnmm0jpff/1169SpkzVq1Mg6d+5sc2rxL2I0Bhyz8i+NNm0Cv1CSk/kyqamDBytvv8O/nJOTwyNIVPVn3aZNePxZR9J2BBBeavP9HfT74ISj2lxH7zTMeTh2ixeXn0Y5mkWLpAEDgj2a6oX7PY8iZTsCCD+1+f6ul6uoED5iYvjSOFaRcDXa0e555HKV3/PokktCF2wjYTsCiHw8bBOooUi4Gi0S7nkUCdsRQOQj4AA1FAlXo0XC0ZFI2I4AIh8BB6ihmBhp2rTy/z/8y7ni9dSpoZ3TFAlHRyJhOwKIfAQcoBZGjSqfpNu6dWB7mzahn7wrRc7RkXDfjgAiH1dRRdlVVKgb4Xw1WsVVVFLgZONwuYrqUOG8HQGEn9p8fxNwCDhwoJyc8qupDp1wnJxcfuonXMINANQWl4kDUW7UqPJLwTk6AiBaEXAAh+KeRwCiGZOMAQCA4xBwAACA4xBwAACA4xBwAACA4xBwAACA4xBwAACA4xBwAACA4xBwAACA4xBwAACA4xBwAACA4xBwAACA4xBwAACA4xBwAACA4/A0cQA4Ap9PWrpUKiqSkpKk9PTyJ7UDCG8EHACoRk6OdPPN0ubNv7S1aSNNmyaNGhW6cQE4Ok5RAUAVcnKk0aMDw40kbdlS3p6TE5pxAagZAg4AHMbnKz9yY1b5vYq2iRPL+wEITwQcADjM0qWVj9wcykzatKm8H4DwRMABgMMUFdVtPwD1j4ADAIdJSqrbfgDqHwEHAA6Tnl5+tZTLVfX7LpeUnFzeD0B4IuAAwGFiYsovBZcqh5yK11Oncj8cIJwRcACgCqNGSW+8IbVuHdjepk15O/fBAcIbN/oDgGqMGiVdcgl3MgYiEQEHAI4gJkYaMCDUowBQW5yiAgAAjkPAAQAAjhPUgLNr1y5lZWXJ7XbL7XYrKytLu3fvrrZ/aWmp7rjjDvXo0UNNmzZVq1at9Jvf/EZbt24N6DdgwAC5XK6AZezYscEsBQAARJCgBpwrrrhCBQUFWrBggRYsWKCCggJlZWVV2/+nn37S6tWr9ac//UmrV69WTk6OvvnmG1188cWV+o4fP15FRUX+5bnnngtmKQAAIIIEbZLxV199pQULFmj58uXq16+fJOlvf/ub0tLStHbtWnXq1KnSZ9xutxYuXBjQ9uSTT6pv374qLCzUaaed5m9v0qSJWrZsGazhAwCACBa0Izh5eXlyu93+cCNJqampcrvdWrZsWY3X4/F45HK51KxZs4D2mTNnqnnz5urWrZtuvfVW7dmzp9p1lJSUyOv1BiwAAMC5gnYEp7i4WC1atKjU3qJFCxUXF9doHfv379edd96pK664QgkJCf72K6+8Uu3bt1fLli31+eefa/LkyVqzZk2loz8VsrOzde+99x5bIQAAIOLU+gjOPffcU2mC7+HLypUrJUmuKh7kYmZVth+utLRUY8eOVVlZmaZPnx7w3vjx43XhhReqe/fuGjt2rN544w299957Wr16dZXrmjx5sjwej3/ZtGlTbcsGAAARpNZHcG688cajXrHUrl07ffrpp9q2bVul93744QclJiYe8fOlpaUaM2aMNm7cqA8++CDg6E1VevfurUaNGmndunXq3bt3pffj4uIUFxd3xHUAAADnqHXAad68uZo3b37UfmlpafJ4PFqxYoX69u0rSfr444/l8XjUv3//aj9XEW7WrVunRYsW6ZRTTjnqz/riiy9UWlqqpKSkmhcCAAAcK2iTjLt06aLBgwdr/PjxWr58uZYvX67x48dr+PDhAVdQde7cWbm5uZKkgwcPavTo0Vq5cqVmzpwpn8+n4uJiFRcX68CBA5Kk9evX67777tPKlSv13Xffad68efr1r3+tXr166eyzzw5WOQAAIIIE9T44M2fOVI8ePZSRkaGMjAz96le/0iuvvBLQZ+3atfJ4PJKkzZs36+2339bmzZvVs2dPJSUl+ZeKK69iY2P1/vvv66KLLlKnTp100003KSMjQ++9955ieAIeAACQ5DIzC/Ug6pvX65Xb7ZbH4znq/B4AABAeavP9zbOoAACA4xBwAACA4wTtRn8AgPrhO+DTZ9OX6qf1RWpyepJ63JCumFjmJCK6EXAAIIItvz1Hpz1+s3r6Nvvbtt7aRoWTpin1kVEhHBkQWpyiAoAItfz2HPV9dLRaHhJuJKmlb4v6Pjpay2/PCdHIgNAj4ABABPId8Om0x2+WZJX+IW+g8otjkx+fKN8BX72PDQgHBBwAiECfTV+qVr7N1f4j3kCm1r5N+mz60nodFxAuCDgAEIF+Wl9Up/0ApyHgAEAEanJ6zZ69V9N+gNMQcAAgAvW4IV1bY9qoTK4q3y+TS1tiktXjhvR6HhkQHgg4ABCBYmJjVDhpmiRVCjkVrzdNmsr9cBC1CDgAEKFSHxmlFbe9oeKY1gHtRTFttOK2N7gPDqIaD9vkYZsAIhx3Mka0qM33N3cyBoAIFxMbo54TB4R6GEBY4RQVAABwHAIOAABwHAIOAABwHAIOAABwHAIOAABwHAIOAABwHAIOAABwHO6DAwAA6ky43HiSgAMAAOrE8ttzdNrjN6unb7O/beutbVQ4aVq9PzqEU1QAAOC4Lb89R30fHa2Wh4QbSWrp26K+j47W8ttz6nU8BBwAAHBcfAd8Ou3xmyVZpWDRQOWPvEx+fKJ8B3z1NiYCDgAg6HwHfCqYuljL/jBLBVMX1+sXHYLvs+lL1cq3udpQ0UCm1r5N+mz60nobE3NwAABBFU7zMhAcP60vqtN+dYEjOACAoAm3eRkIjianJ9Vpv7rgMjOrt58WJrxer9xutzwejxISEkI9HABwJN8Bn7Y1aaeW1Zy6KJNLRTFt1PKnjSG5jBh155c/6y3+OTeHqqs/69p8f3MEBwAQFOE4LwPBERMbo8JJ0ySVh5lDVbzeNGlqvQZZAg4AICjCcV4Ggif1kVFacdsbKo5pHdBeFNNGK257o97nWzHJGAAQFOE4LwPBlfrIKPkeuEQFh93JuHUITkEyB4c5OAAQFPU1L6OuhMsjBlA95uAAAEIuHOdlVGf57Tna1qSdet4yUP2fukI9bxmobU3acZVXBCPgAACCJtzmZVSFS9mdiVNUnKICgKAL19M/XMoeWWrz/c0kYwBA0MXExqjnxAGhHkYln01fGnCH5cNVXMpeMH1pWIw/XINiOArqKapdu3YpKytLbrdbbrdbWVlZ2r179xE/c/XVV8vlcgUsqampAX1KSkr0hz/8Qc2bN1fTpk118cUXa/Pm6v+CAgBQlUi6lJ15QrUT1IBzxRVXqKCgQAsWLNCCBQtUUFCgrKyso35u8ODBKioq8i/z5s0LeH/ixInKzc3Vq6++qo8++kh79+7V8OHD5fPx8DYAQM1FyqXszBOqvaDNwfnqq6/UtWtXLV++XP369ZMkLV++XGlpafr666/VqVOnKj939dVXa/fu3XrzzTerfN/j8ejUU0/VK6+8oszMTEnS1q1blZycrHnz5umiiy466tiYgwMAkCLjUnbmCf0iLC4Tz8vLk9vt9ocbSUpNTZXb7dayZcuO+NnFixerRYsWOuOMMzR+/Hht377d/96qVatUWlqqjIwMf1urVq3UvXv3atdbUlIir9cbsAAAEAmXsvPIi2MTtIBTXFysFi1aVGpv0aKFiouLq/3ckCFDNHPmTH3wwQd67LHH9Mknn+j8889XSUmJf72xsbE66aSTAj6XmJhY7Xqzs7P984DcbreSk5OPozIAgJOE+6XskTRPKJzUOuDcc889lSYBH76sXLlSkuRyuSp93syqbK+QmZmpYcOGqXv37hoxYoTmz5+vb775RnPnzj3iuI603smTJ8vj8fiXTZs21aJiAIDTpT4ySok/faeCvy7Sshv/qYK/LlLLnzaGPNxIkTNPKNzU+jLxG2+8UWPHjj1in3bt2unTTz/Vtm3bKr33ww8/KDExscY/LykpSW3bttW6deskSS1bttSBAwe0a9eugKM427dvV//+/atcR1xcnOLi4mr8MwEA0SdcL2XvcUO6tt7a5qjzhHrckB6C0YWvWgec5s2bq3nz5kftl5aWJo/HoxUrVqhv376SpI8//lgej6faIFKVnTt3atOmTUpKKk+mKSkpatSokRYuXKgxY8ZIkoqKivT555/rkUceqW05AACEtYp5Qi0fHa0yuQJCzqHzhELxQMtwFrQ5OF26dNHgwYM1fvx4LV++XMuXL9f48eM1fPjwgCuoOnfurNzcXEnS3r17deuttyovL0/fffedFi9erBEjRqh58+a69NJLJUlut1vjxo3T//t//0/vv/++8vPzddVVV6lHjx668MILg1UOAAAhE+7zhMJRUO9kPHPmTN10003+K54uvvhiPfXUUwF91q5dK4/HI0mKiYnRZ599ppdfflm7d+9WUlKSBg4cqNmzZ+vEE0/0f+avf/2rGjZsqDFjxujnn3/WBRdcoBkzZigmhvQKAHCm1EdGyffAJSo47E7GHLmpGs+i4j44AABEhLC4Dw4AAECoEHAAAIDjEHAAAIDjEHAAAIDjEHAAAIDjEHAAAIDjEHAAAIDjEHAAAIDjEHAAAIDjBPVRDeGq4ubNXq83xCMBAAA1VfG9XZOHMERlwNmzZ48kKTk5OcQjAQAAtbVnzx653e4j9onKZ1GVlZVp69atOvHEE+VyuUI9nOPm9XqVnJysTZs2OfbZWtFQoxQ9dQZbNGxHanSOaKmzLpiZ9uzZo1atWqlBgyPPsonKIzgNGjRQmzZtQj2MOpeQkOD4nSMaapSip85gi4btSI3OES11Hq+jHbmpwCRjAADgOAQcAADgOAQcB4iLi9Of//xnxcXFhXooQRMNNUrRU2ewRcN2pEbniJY661tUTjIGAADOxhEcAADgOAQcAADgOAQcAADgOAQcAADgOAScOrRv374aPR8jkkVDjVL01Bls0bAdqdE5oqXOaEHAqSM+n0/jxo3T0KFD9fzzz4d6OEERDTVK0VNnsEXDdqRG54iWOqMJAaeO7N69Wz179tSgQYP00EMP6ZxzztG6detCPaw6FQ01StFTZ7BFw3akRueIljqjCffBCYK9e/dq9OjRSk1N1T333BPq4QRFNNQoRU+dwRYN25EanSNa6nQ6juDUkYqcePDgQZ1wwglKT0/XSy+9pL179wa8H8mioUYpeuoMtmjYjtTojBql6KkzmhBw6oCZyeVySZIaNix/QPvy5cvVqVMnnXDCCZLkf3/Hjh2hGeRxioYapeipM9iiYTtSozNqlKKnzmjTMNQDcIKKv/hbt27V+++/r1dffVUfffSR5s+fr7KyMjVo0EBz587Va6+9pvXr10uSHn/8cfXt2zeUw66VaKhRip46gy0atiM1OqNGKXrqjDbMwTkOXq9Xn376qRYsWKBPPvlEX3zxhWJjYzV69GiNGDFC6enpkqSFCxdqzJgxSklJUWZmpr7++mu98sorevXVV3X++eeHuIoji4YapeipM9iiYTtSozNqlKKnzqhlOCYlJSXWrFkza9asmY0dO9Yee+wx+/DDD83n8wX027BhgyUnJ9u1115ru3fvNjOzgwcP2vDhw+3WW281M7OysrJ6H39NREONZtFTZ7BFw3akxl9Eco1m0VNnNCPgHKMffvjB+vfvby6Xy1588cWA9w7dQW6++WZr3bq1eTweM/tlRxg7dqwNHjw44HPbtm0L8qhrJxpqNIueOoMtGrYjNTqjRrPoqTOaEXCO0+uvv26nnnqqde7c2ebOnRvw3rZt26xJkyb+nefgwYNmZlZUVGSpqak2efJk/460dOlSGz9+vI0fP96+/PLL+i3iKKKhRrPoqTPYomE7UqMzajSLnjqjEVdRHafRo0dr+/btuuyyy3TppZcqJydHZWVlkqR//OMfatGihUaOHCkzU0xMjCRp7ty5MjN17dpVDRqU/xGceOKJ6tSpkwoLC9WtWze9+eaboSqpkmioUYqeOoMtGrYjNTqjRil66oxKoUhVTrVnzx4rLCw0s/LDmPfee6+lpqaa1+v19/niiy/skksusWHDhtnevXv9fc3Kzwnfc8891rhx47A91BkNNZpFT53BFg3bkRrLRXqNZtFTZ7TgCE4dOuGEE5ScnCyp/LLD7t2768cff1RJSYm/z0MPPaSdO3dq/Pjxatq0qXw+n1wul8xMGzZs0JNPPqm7775bLVq08P8WEU6ioUYpeuoMtmjYjtRYLtJrlKKnzqgRwnDleIWFhZaammo9e/a0yZMn21lnnWVJSUn2wgsv+PtUJP/9+/fbuHHjrH379gHrqDjn6/V6bePGjWH3W0E01GgWPXUGWzRsR2osF+k1mkVPnU5FwKkHU6ZMseHDh9stt9xia9as8bcfOlN/0aJF1rBhQ8vNzfW/V7HjbN261Xr06GG9e/e20047zW688cawuywxGmo0i546gy0atiM1lov0Gs2ip06nIeDUk9LS0mrf2717tw0cONDOP//8gPaKHWDUqFGWkpJib7/9tuXl5dnZZ59t559/vv3www9BHXNtRUONZtFTZ7BFw3akRmfUaBY9dToJASdEVqxYYbm5uebxeOyVV16xBg0a2BdffGFm5cm/4jeDsrIyGzNmjN15553+z3755Zd2xhln2Lx580Iy9pqKhhrNoqfOYIuG7UiNzqjRLHrqjGQ8iypEDh48qHHjxqlZs2bau3evJkyYoK5du/qfe1Jh+fLl6ty5sxYsWCCPxyO3262WLVtq3bp12r9/fwgrOLpoqFGKnjqDLRq2IzU6o0YpeuqMZFxFFSJpaWnavn27Lr30Uv3www/67rvv9PXXX/t3jN27d+umm27S6NGjtWbNGjVu3FgtWrTQJZdcogEDBqhnz55h/6C3aKhRip46gy0atiM1OqNGKXrqjGihPoQEs40bN1pKSor17dvXNmzYYGblhz8bNmxoixYt8vebPn26nXnmmbZq1aqA+zJEgmio0Sx66gy2aNiO1LjI3y+SazSLnjojDQEnjKxfv95/46i//e1v1rhx44AbSf3888+WmJhor7/+esDnDh48WGlGvtfrtQ0bNti3335bP4OvoWio0Sx66gy2aNiO1OiMGs2ip85IwSmqMNKhQwc1bdpUknTeeefpjDPO8N/u2+VyafPmzYqPj1eTJk0kSX/961/l8XgUExMjl8slSfL5fFq+fLn69euniy++WGeddZbGjh2rffv2haSmw0VDjVL01Bls0bAdqdEZNUrRU2fECHXCQtXKysrsvvvus2bNmtn48ePtoYcesjZt2tjw4cNt5cqVVlpaasOHD7dmzZrZAw884P/czJkzLSUlxQYNGmT/+c9/bNWqVXbeeefZxIkTw+6+C9FQo1n01Bls0bAdqdEZNZpFT53hjIAT5j777DMbOXKkjRo1yn77299WOlz55ptvWp8+fWzWrFn2/fff23nnnWdjx44NuL/Cgw8+aF26dLE9e/bU9/BrJBpqNIueOoOtLrdjxemDcEONzqjRLHrqDEcEnAixb9++at/bv3+/mZndfffdlpaWZv/+97/N7Je7bL722mvWsWNH/+S3/Px8u+666+zzzz8P8qhrJxpqNIueOoOtLrbj+vXrAz73888/B2m0x4YanVGjWfTUGU6YgxMhKs7ZViUuLk7r1q1TXl6ezjrrLA0ZMkSSZGaSpHfffVdNmzZV+/btJUmdO3dWkyZN1KNHD1199dVhcy+GaKhRip46g+14t2OTJk3UoUMHSdJnn32ma6+9VldccYVuuummgIcrhhI1OqNGKXrqDCshDFeoQ1988YW1aNHC5s+fb2a/PODt22+/tdjYWHv55ZcD2s3KD50OGDDA4uLibPr06fU/6FqKhhrNoqfOYDvadvzf//1fMzNbu3atXXTRRXb66afbk08+aeeff761bt3a/7lwRo3OqNEseuqsTwQch1i8eLHFxMRYUVGRv83n89nw4cOtb9++tnPnzoD+h3453nbbbeZyuezKK6+st/Eei2io0Sx66gy2o23HHTt2mFn5F4bL5bJPPvnE3+93v/udZWZm1vuYa4sanVGjWfTUWZ8IOA7h8XgsLS3Npk6damblqf+BBx6wRo0a2X/+85+AvhXndfft22evvPKKnXLKKXbllVfahx9+GPB+uImGGs2ip85gq247NmzY0D766CN/v++++84uv/xye+edd/xtU6dOtQ4dOoT9PUio0Rk1mkVPnfWJgOMgU6dOtbi4OMvIyLDWrVtbt27d7Jlnnqm2/+9+9zs744wz7KqrrqrHUR6faKjRLHrqDLaqtuOTTz5pZmZfffWV3XXXXZaUlGQpKSnWoEEDu/HGG23atGl29tln23nnnRfawdcQNTqjRrPoqbO+EHAcZvPmzfbAAw/Y3//+d/vmm2/87WVlZf7f5levXm1XX321NWrUyF544QX/pYeR8tt+NNRoFj11Btuh23Ht2rX+9htvvNHOO+88mzlzpu3cudPef/99O/300+3iiy+2F154wd/30FOA4YoanVGjWfTUWR8IOFGg4i98WVmZffrpp9ahQwc766yzLDc3t8q+FTeTOvymUuF875VjrbHCvn37bNmyZbZu3br6GO4xO9Y6KwLP+vXr7aWXXrJPP/203sYcbsrKymzPnj02cuRIGzNmTMB7f/7zn23YsGEBbZG4T1DjkWusEAn7/fHUGe37PZeJR4GYmBhJ0mOPPaarrrpKycnJeuWVVzRy5Eh/nx9//DHgluGlpaX+W4d/+OGHysrK0qhRo5SZmakffvghFGUc0bHU6PP5/O998803+vvf/66RI0fq8ssv1+7du+u5gpo51jornnA8depUzZw5U0OGDNFll12mnTt3hqKMkHK5XDrhhBOUmJion376yX8priT17dtXq1at0q5du+T1eiN2n6DG6muMtP3+eOqM+v0+tPkK9SU/P99iY2NtwoQJtmvXrkrvL1iwwGJjY23KlCkB7c8++6x17tzZUlJS7Pnnn7eLL77Y2rVrZ19//XU9jbzmjrVGs/LfelavXm2vvvqquVwue/DBB62kpKQeRl17x1Onz+ezoqIiW7Nmjf3617+2c845x7Zu3VoPow4/H3zwgbVr186uueYa++KLL2zWrFnWvXt3GzlypJmZzZ07N+L3CWqsvkazyNrvj6fOaN3vCThR5IcffjCv12tmVc/R+Ne//mVdunSxAQMG2I8//mgbN260uLg4u/POO/2XKHq9XktNTfVPfAs3Na1x4MCBtm3btkrv33jjjdapUyf77LPPgj7W41GbOouLi6vst3btWjvjjDNszpw5wR9wmFq3bp1dcMEF9l//9V/WqVMnS01NDbgM3wn7BDU6Z7+vTZ3s9wQcHGb//v2Wl5dnZWVlNnr0aOvbt69t3749oE+vXr3stttuC9EIj9/+/fttxYoV5vF4zOyXeS3r1683l8tl06dPt9LS0oD3du3aZZs3bw7NgI9RRZ0VQcis8gTE5s2b27Rp0/yvPR5PwDNwosXnn39uhYWF/m116HY6ln3iwIEDYTd3pa5rDEdHq9Ep+31N6mS/J+DgMBVpf82aNeZyuWzu3LkB73/++efWv39/e/DBB/1tK1eutBkzZtjjjz8esFNFiooJeRdeeKGdc845tmXLlkp9srOz7ZRTTrHLL7/c/w9kJDl8cuWGDRvs1ltvtYSEhIDnWP3ud7+z9PR0mzBhQkT+WQbDsewT7777rv33f/+3/fGPf4yI7ch+z37vxP2eScYIUDEp7bnnnlN6err69OkT8P6yZcv0448/qlevXpKk+fPn68orr1R2drb+/e9/q3Xr1po6dWp9D/uY+Xw+uVwuLViwQB988IHuuOMOtWrVyv+eJH366afKy8tTamqqioqKdOqpp+qee+4J4ahrr7S0VGvXrtULL7ygrKwspaen69///rdeeeUVdevWTZJ08OBBnXnmmbrhhhu0YcMGnXHGGZozZ06IRx56Nd0nKtp9Pp8aN26s5s2ba8GCBXK73VqzZk29j7s22O/Z7x2534c6YSH8lJWV2U033WQXXXRRQPsnn3xiGRkZNmLECH/bhAkTbPDgwbZmzRozM5s9e7b96le/soULF9brmGurrKzM/wRfM7MOHTrYVVdd5Z+0e+hvPrfffrt16dLFf4llbm6unX766XbppZeG7T0nKsa1YsUKu+WWW+yCCy6wZs2a2Zlnnml33HGH5eTk+O+Zc/hveRWmTJliQ4cODdsa61Nt9olDXXPNNdahQ4eImNDJfs9+b+as/Z6AgyplZ2dbamqq/y/5/v377brrrrPu3bsHPAMlNzfXBg8eHPDZjh072l133VWv462tb7/91kaMGGHz58+3Rx991E466SRbsWJFpX4bNmyw22+/3Vq3bm0jRoywDRs2+N+ruHdGON9U79xzzzWXy2VXXXWVfffdd0fsW1ZWFvCP3nvvvWcnnniiffzxx8EeZkQ40j6xevVqMzMrLS31b8NPP/3UGjZsaDNnzgzZmGuL/b4c+70z9nsCDqq0Zs0aa9eunQ0bNsweffRR69Kli/Xu3dv/pOpt27bZww8/bEOHDjW3223jx4+3Xbt2WW5urnXp0sUef/zxan9DCBe33HKLNWrUyFwul1177bUBv9kdavv27fbFF1/YuHHjrE+fPlZQUFDPIz0+jz76qMXExNill15q33//fbX9Dv1y/vnnn+2xxx6zZs2a2b59++prqGGtun3iueeeM7PKvxEPGjTIBgwYUO3fq3DEfv8L9vvI3+8JOKjW9u3b7fLLL7fhw4fbNddcY99++63/KoO0tDRLTU21Bx980N566y3r06ePNW7c2AYNGmSTJk2KmIe+7dixw6655hpzuVz2hz/8odKTug+1ceNG69Gjhz388MP1OMK6UVxcbCNGjDCXy2X33ntvpTudVtixY4fl5ubakCFD7JRTTrHnn38+FMMNW1XtExXb8NBtmZOTYzExMQEPSYwU7PeB2O8jFwEHR/Xzzz8HvM7Pz7dGjRpZXl6ev2379u2Wmppq8+fPr+/h1YmlS5fapZdeah988IGZlR+iPvSqiop/EDIzMy0jIyMkY6wLixYtssmTJ1thYaGZld9PZ8WKFfbHP/7RzjvvPGvdurV1797drr32WluwYEGIRxu+Dt8nDrV//37r2LGjjRs3rh5HVPfY79nvIx0BB7WWn59vrVu3tn/+858B7eecc4795S9/CdGo6tbVV19t9913X0Dbvn37rGPHjvanP/0pRKOqW16v17p06WIul8suv/xyu//++23ZsmWhHlZEWr16tS1evNjMzB5++GE76aSTAuZtOAH7Pft9pOEycdRa586dlZaWptdff11btmyRz+dTbm6uXC6Xfvrpp1AP77jY/z3npVevXnriiSd0xRVXqKCgQB988IGuueYaNWrUSOecc06IR1k3TjzxRA0dOlQxMTGKj4/Xf//3fystLU1S+eWlqLnCwkINHDhQ559/vu6//35NmTJF7du3D/Ww6hT7Pft9xAl1wkJkKi4utqFDh9pJJ51kQ4YMsSZNmtjo0aP9h0Cd4Pvvv7dhw4ZZ586drWPHjpaYmGhz586NqEmjNVFYWGjp6enWpEkTu+eee8J+kmi42rhxo1122WXmcrns+uuvt6KiolAPqc6x3ztHNOz3LrNDHk0K1NJHH32k1atX6/TTT1dKSopatmwZ6iHVuW+++Ubx8fFq2rSpTjnllFAPJ2jeeecdZWZm6g9/+IP+/Oc/q2HDhqEeUkT64IMPdOWVV6pfv3569tlnHblPsN87h5P3ewIOgADbtm1TYmJiqIcR8VavXq2ePXv67xIMhDMn7vcEHAAA4Dj8agEAAByHgAMAAByHgAMAAByHgAMAAByHgAMAAByHgAMAAByHgAMAAByHgAMAAByHgAMAAByHgAMAAByHgAMAABzn/wOZYdIZzKl5qgAAAABJRU5ErkJggg==", + "text/plain": [ + "
" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "pca = decomposition.PCA(n_components=1)\n", + "pca_result = pca.fit_transform(embeddings_mean)\n", + "\n", + "plt.xticks(rotation=-30)\n", + "# All points\n", + "plt.scatter(stack.time, pca_result, color=\"blue\")\n", + "\n", + "# Cloudy images\n", + "plt.scatter(stack.time[0], pca_result[0], color=\"green\")\n", + "plt.scatter(stack.time[2], pca_result[2], color=\"green\")\n", + "\n", + "# After fire\n", + "plt.scatter(stack.time[-5:], pca_result[-5:], color=\"red\")" + ] + }, + { + "cell_type": "markdown", + "id": "a16fbdb8-1c2d-4c84-8526-283fa14faa53", + "metadata": {}, + "source": [ + "In the plot above, each image embedding is one point. One can clearly \n", + "distinguish the two cloudy images and the values after the fire are\n", + "consistently low." + ] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python 3 (ipykernel)", + "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.11.6" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +} diff --git a/src/datamodule.py b/src/datamodule.py index 8a7d6d5c..e653bb49 100644 --- a/src/datamodule.py +++ b/src/datamodule.py @@ -157,7 +157,7 @@ def setup(self, stage: Literal["fit", "predict"] | None = None) -> None: dp = torchdata.datapipes.iter.IterableWrapper(iterable=[self.data_dir]) chips_path = list(dp.list_files_by_s3(masks="*.tif")) else: # if self.data_dir is a local data path - chips_path = list(Path(self.data_dir).glob("**/*.tif")) + chips_path = sorted(list(Path(self.data_dir).glob("**/*.tif"))) print(f"Total number of chips: {len(chips_path)}") if stage == "fit": diff --git a/src/model_clay.py b/src/model_clay.py index cffea57e..7b9dc134 100644 --- a/src/model_clay.py +++ b/src/model_clay.py @@ -791,6 +791,7 @@ def __init__( # noqa: PLR0913 b1=0.9, b2=0.95, embeddings_level: Literal["mean", "patch", "group"] = "mean", + band_groups=None, ): super().__init__() self.save_hyperparameters(logger=True) @@ -801,11 +802,14 @@ def __init__( # noqa: PLR0913 "large": clay_large, } if model_size in model_map: - self.model = model_map[model_size]( - mask_ratio=mask_ratio, - image_size=image_size, - patch_size=patch_size, - ) + model_args = { + "mask_ratio": mask_ratio, + "image_size": image_size, + "patch_size": patch_size, + } + if band_groups: + model_args["band_groups"] = band_groups + self.model = model_map[model_size](**model_args) else: raise ValueError( f"Invalid model size {model_size}. Expected one of {model_map.keys()}"