{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "0cc5e729-9116-4ec9-bf1e-8346cbccdf7b",
   "metadata": {},
   "source": [
    "## Run Clay v1\n",
    "\n",
    "This notebook shows how to run Clay v1 wall-to-wall, from downloading imagery\n",
    "to training a tiny fine tuning head. This will include the following steps:\n",
    "\n",
    "1. Set a location and date range of interest\n",
    "2. Download Sentinel-2 imagery for this specification\n",
    "3. Load the model checkpoint\n",
    "4. Prepare data into a format for the model\n",
    "5. Run the model on the imagery\n",
    "6. Analyise the model embeddings output using PCA\n",
    "7. Train a Support Vector Machines fine tuning head"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "add63cd9",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Add the repo root to the sys path for the model import below\n",
    "import sys\n",
    "\n",
    "sys.path.append(\"..\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "6a17b8a8-a9c6-4053-833e-de97287fae49",
   "metadata": {},
   "outputs": [],
   "source": [
    "import math\n",
    "\n",
    "import geopandas as gpd\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "import pystac_client\n",
    "import stackstac\n",
    "import torch\n",
    "import yaml\n",
    "from box import Box\n",
    "from matplotlib import pyplot as plt\n",
    "from rasterio.enums import Resampling\n",
    "from shapely import Point\n",
    "from sklearn import decomposition, svm\n",
    "from torchvision.transforms import v2\n",
    "\n",
    "from src.model import ClayMAEModule"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "beac6394-9762-422b-9f5d-82d226018c0c",
   "metadata": {},
   "source": [
    "### Specify location and date of interest\n",
    "In this example we will use a location in Portugal where a forest fire happened. We will run the model over the time period of the fire and analyse the model embeddings."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "08d7787d-1506-4de7-89dc-c1054910acf7",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Point over Monchique Portugal\n",
    "lat, lon = 37.30939, -8.57207\n",
    "\n",
    "# Dates of a large forest fire\n",
    "start = \"2018-07-01\"\n",
    "end = \"2018-09-01\""
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2bd226c9-003b-4867-a64a-8ae887e7e20a",
   "metadata": {},
   "source": [
    "### Get data from STAC catalog\n",
    "\n",
    "Based on the location and date we can obtain a stack of imagery using stackstac. Let's start with finding the STAC items we want to analyse."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "2e80743c-7c77-459b-9984-f6c26cdff549",
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/home/tam/apps/miniforge3/envs/claymodel/lib/python3.11/site-packages/pystac_client/item_search.py:850: FutureWarning: get_all_items() is deprecated, use item_collection() instead.\n",
      "  warnings.warn(\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Found 12 items\n"
     ]
    }
   ],
   "source": [
    "STAC_API = \"https://earth-search.aws.element84.com/v1\"\n",
    "COLLECTION = \"sentinel-2-l2a\"\n",
    "\n",
    "# Search the catalogue\n",
    "catalog = pystac_client.Client.open(STAC_API)\n",
    "search = catalog.search(\n",
    "    collections=[COLLECTION],\n",
    "    datetime=f\"{start}/{end}\",\n",
    "    bbox=(lon - 1e-5, lat - 1e-5, lon + 1e-5, lat + 1e-5),\n",
    "    max_items=100,\n",
    "    query={\"eo:cloud_cover\": {\"lt\": 80}},\n",
    ")\n",
    "\n",
    "all_items = search.get_all_items()\n",
    "\n",
    "# Reduce to one per date (there might be some duplicates\n",
    "# based on the location)\n",
    "items = []\n",
    "dates = []\n",
    "for item in all_items:\n",
    "    if item.datetime.date() not in dates:\n",
    "        items.append(item)\n",
    "        dates.append(item.datetime.date())\n",
    "\n",
    "print(f\"Found {len(items)} items\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5b7c68ae-7c8a-446a-8bc7-5afba70183c2",
   "metadata": {},
   "source": [
    "### Create a bounding box around the point of interest\n",
    "\n",
    "This is needed in the projection of the data so that we can generate image chips of the right size."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "0f3573b5-5a00-47d9-a648-5c4d7cd2c996",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Extract coordinate system from first item\n",
    "epsg = items[0].properties[\"proj:epsg\"]\n",
    "\n",
    "# Convert point of interest into the image projection\n",
    "# (assumes all images are in the same projection)\n",
    "poidf = gpd.GeoDataFrame(\n",
    "    pd.DataFrame(),\n",
    "    crs=\"EPSG:4326\",\n",
    "    geometry=[Point(lon, lat)],\n",
    ").to_crs(epsg)\n",
    "\n",
    "coords = poidf.iloc[0].geometry.coords[0]\n",
    "\n",
    "# Create bounds in projection\n",
    "size = 256\n",
    "gsd = 10\n",
    "bounds = (\n",
    "    coords[0] - (size * gsd) // 2,\n",
    "    coords[1] - (size * gsd) // 2,\n",
    "    coords[0] + (size * gsd) // 2,\n",
    "    coords[1] + (size * gsd) // 2,\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "bbbd3f67-5f2c-46dc-9ee1-2ef1f50fa032",
   "metadata": {},
   "source": [
    "### Retrieve the imagery data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "8b8d3824-e48c-4f9d-9c7b-181c0800f96f",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "<xarray.DataArray 'stackstac-91f577dfbc973b2eff4d609fa1d57243' (time: 12,\n",
      "                                                                band: 4,\n",
      "                                                                y: 256, x: 256)> Size: 13MB\n",
      "dask.array<fetch_raster_window, shape=(12, 4, 256, 256), dtype=float32, chunksize=(1, 1, 256, 256), chunktype=numpy.ndarray>\n",
      "Coordinates: (12/53)\n",
      "  * time                                     (time) datetime64[ns] 96B 2018-0...\n",
      "    id                                       (time) <U24 1kB 'S2B_29SNB_20180...\n",
      "  * band                                     (band) <U5 80B 'blue' ... 'nir'\n",
      "  * x                                        (x) float64 2kB 5.366e+05 ... 5....\n",
      "  * y                                        (y) float64 2kB 4.131e+06 ... 4....\n",
      "    platform                                 (time) <U11 528B 'sentinel-2b' ....\n",
      "    ...                                       ...\n",
      "    gsd                                      int64 8B 10\n",
      "    proj:transform                           object 8B {0, 4200000, 10, 49998...\n",
      "    common_name                              (band) <U5 80B 'blue' ... 'nir'\n",
      "    center_wavelength                        (band) float64 32B 0.49 ... 0.842\n",
      "    full_width_half_max                      (band) float64 32B 0.098 ... 0.145\n",
      "    epsg                                     int64 8B 32629\n",
      "Attributes:\n",
      "    spec:        RasterSpec(epsg=32629, bounds=(536640.79691545, 4128000.7407...\n",
      "    crs:         epsg:32629\n",
      "    transform:   | 10.00, 0.00, 536640.80|\\n| 0.00,-10.00, 4130560.74|\\n| 0.0...\n",
      "    resolution:  10\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/home/tam/apps/miniforge3/envs/claymodel/lib/python3.11/site-packages/stackstac/prepare.py:408: UserWarning: The argument 'infer_datetime_format' is deprecated and will be removed in a future version. A strict version of it is now the default, see https://pandas.pydata.org/pdeps/0004-consistent-to-datetime-parsing.html. You can safely remove this argument.\n",
      "  times = pd.to_datetime(\n"
     ]
    }
   ],
   "source": [
    "# Retrieve the pixel values, for the bounding box in\n",
    "# the target projection. In this example we use only\n",
    "# the RGB and NIR bands.\n",
    "stack = stackstac.stack(\n",
    "    items,\n",
    "    bounds=bounds,\n",
    "    snap_bounds=False,\n",
    "    epsg=epsg,\n",
    "    resolution=gsd,\n",
    "    dtype=\"float32\",\n",
    "    rescale=False,\n",
    "    fill_value=0,\n",
    "    assets=[\"blue\", \"green\", \"red\", \"nir\"],\n",
    "    resampling=Resampling.nearest,\n",
    ")\n",
    "\n",
    "print(stack)\n",
    "\n",
    "stack = stack.compute()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "77354bee-bea3-43e0-8936-5808b352e25f",
   "metadata": {},
   "source": [
    "### Let's have a look at the imagery we just downloaded\n",
    "\n",
    "The imagery will contain 7 dates before the fire, of which two are pretty cloudy images. There are also 5 images after the forest fire."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "008b9c7d",
   "metadata": {},
   "outputs": [],
   "source": [
    "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": "e3c5f64a",
   "metadata": {},
   "source": [
    "![Minicube visualization](https://github.com/Clay-foundation/model/assets/901647/c6e924e5-6ba1-4924-b99a-df8b90731a5f)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a7aa15e9-0285-4cac-816d-a58d2ceda389",
   "metadata": {},
   "source": [
    "### Load the model\n",
    "\n",
    "We now have the data to analyse, let's load the model."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "f3fb9b6f-00e4-45a1-b575-2484b5afd511",
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/home/tam/apps/miniforge3/envs/claymodel/lib/python3.11/site-packages/lightning/pytorch/utilities/migration/utils.py:55: The loaded checkpoint was produced with Lightning v2.2.4, which is newer than your current Lightning version: v2.1.4\n",
      "/home/tam/apps/miniforge3/envs/claymodel/lib/python3.11/site-packages/torch/nn/modules/transformer.py:282: UserWarning: enable_nested_tensor is True, but self.use_nested_tensor is False because encoder_layer.self_attn.batch_first was not True(use batch_first for better inference performance)\n",
      "  warnings.warn(f\"enable_nested_tensor is True, but self.use_nested_tensor is False because {why_not_sparsity_fast_path}\")\n"
     ]
    }
   ],
   "source": [
    "device = torch.device(\"cuda\") if torch.cuda.is_available() else torch.device(\"cpu\")\n",
    "ckpt = \"https://clay-model-ckpt.s3.amazonaws.com/v0.5.7/mae_v0.5.7_epoch-13_val-loss-0.3098.ckpt\"\n",
    "torch.set_default_device(device)\n",
    "\n",
    "model = ClayMAEModule.load_from_checkpoint(\n",
    "    ckpt, metadata_path=\"../configs/metadata.yaml\", shuffle=False, mask_ratio=0\n",
    ")\n",
    "model.eval()\n",
    "\n",
    "model = model.to(device)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9422d855-c73b-43c3-9a06-fd89d7f4eb08",
   "metadata": {},
   "source": [
    "### Prepare band metadata for passing it to the model\n",
    "\n",
    "This is the most technical part so far. We will take the information in the stack of imagery and convert it into the formate that the model requires. This includes converting the lat/lon and the date of the imagery into normalized values.\n",
    "\n",
    "The Clay model will accept any band combination in any order, from different platforms. But for this the model needs to know the wavelength of each band that is passed to it, and normalization parameters for each band as well. It will use that to normalize the data and to interpret each band based on its central wavelength.\n",
    "\n",
    "For Sentinel-2 we can use medata file of the model to extract those values. But this cloud also be something custom for a different platform."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "1e9cb385-1aa6-453f-9a5b-1b388a95c407",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Extract mean, std, and wavelengths from metadata\n",
    "platform = \"sentinel-2-l2a\"\n",
    "metadata = Box(yaml.safe_load(open(\"../configs/metadata.yaml\")))\n",
    "mean = []\n",
    "std = []\n",
    "waves = []\n",
    "# Use the band names to get the correct values in the correct order.\n",
    "for band in stack.band:\n",
    "    mean.append(metadata[platform].bands.mean[str(band.values)])\n",
    "    std.append(metadata[platform].bands.std[str(band.values)])\n",
    "    waves.append(metadata[platform].bands.wavelength[str(band.values)])\n",
    "\n",
    "# Prepare the normalization transform function using the mean and std values.\n",
    "transform = v2.Compose(\n",
    "    [\n",
    "        v2.Normalize(mean=mean, std=std),\n",
    "    ]\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c2d8e1f3-011d-4be5-8071-547f0ad91ad6",
   "metadata": {},
   "source": [
    "### Convert the band pixel data in to the format for the model\n",
    "\n",
    "We will take the information in the stack of imagery and convert it into the formate that the model requires. This includes converting the lat/lon and the date of the imagery into normalized values."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "a5bbe0c2-2cc3-428c-8d38-7e9516b6134a",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Prep datetimes embedding using a normalization function from the model code.\n",
    "def normalize_timestamp(date):\n",
    "    week = date.isocalendar().week * 2 * np.pi / 52\n",
    "    hour = date.hour * 2 * np.pi / 24\n",
    "\n",
    "    return (math.sin(week), math.cos(week)), (math.sin(hour), math.cos(hour))\n",
    "\n",
    "\n",
    "datetimes = stack.time.values.astype(\"datetime64[s]\").tolist()\n",
    "times = [normalize_timestamp(dat) for dat in datetimes]\n",
    "week_norm = [dat[0] for dat in times]\n",
    "hour_norm = [dat[1] for dat in times]\n",
    "\n",
    "\n",
    "# Prep lat/lon embedding using the\n",
    "def normalize_latlon(lat, lon):\n",
    "    lat = lat * np.pi / 180\n",
    "    lon = lon * np.pi / 180\n",
    "\n",
    "    return (math.sin(lat), math.cos(lat)), (math.sin(lon), math.cos(lon))\n",
    "\n",
    "\n",
    "latlons = [normalize_latlon(lat, lon)] * len(times)\n",
    "lat_norm = [dat[0] for dat in latlons]\n",
    "lon_norm = [dat[1] for dat in latlons]\n",
    "\n",
    "# Normalize pixels\n",
    "pixels = torch.from_numpy(stack.data.astype(np.float32))\n",
    "pixels = transform(pixels)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "fe34640d-48df-41b0-a749-4d6ada1a42d7",
   "metadata": {},
   "source": [
    "### Combine the metadata and the transformed pixels\n",
    "\n",
    "Now we can combine all of these inputs into a dictionary that combines everything."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "70fea0d4-f3ac-4430-a7d0-7b3c302a0754",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Prepare additional information\n",
    "datacube = {\n",
    "    \"platform\": platform,\n",
    "    \"time\": torch.tensor(\n",
    "        np.hstack((week_norm, hour_norm)),\n",
    "        dtype=torch.float32,\n",
    "        device=device,\n",
    "    ),\n",
    "    \"latlon\": torch.tensor(\n",
    "        np.hstack((lat_norm, lon_norm)), dtype=torch.float32, device=device\n",
    "    ),\n",
    "    \"pixels\": pixels.to(device),\n",
    "    \"gsd\": torch.tensor(stack.gsd.values, device=device),\n",
    "    \"waves\": torch.tensor(waves, device=device),\n",
    "}"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4a18e52b-b21e-4c58-a1f2-26b66d73ecbe",
   "metadata": {},
   "source": [
    "### Run the model\n",
    "\n",
    "Pass the datacube we prepared to the model to create embeddings. This will create one embedding vector for each of the images we downloaded."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "69335005-ebd1-4edf-b493-5e314637ef68",
   "metadata": {},
   "outputs": [],
   "source": [
    "with torch.no_grad():\n",
    "    unmsk_patch, unmsk_idx, msk_idx, msk_matrix = model.model.encoder(datacube)\n",
    "\n",
    "# The first embedding is the class token, which is the\n",
    "# overall single embedding. We extract that for PCA below.\n",
    "embeddings = unmsk_patch[:, 0, :].cpu().numpy()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "62b3ebb8-18a2-4918-b863-01ea36095d9f",
   "metadata": {},
   "source": [
    "### Analyse the embeddings\n",
    "\n",
    "A simple analysis of the embeddings is to reduce each one of them into a single number using Principal Component Analysis. For this we will fit a PCA on the 12 embeddings we have, and do the dimensionality reduction for them. We will se a separation into three groups, the previous images, the cloudy images, and the images after the fire, they all fall into a different range of the PCA space."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "id": "e8debf91-e38a-46d2-81c9-24b71a3adfc4",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.collections.PathCollection at 0x7f847cda1890>"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAi8AAAHSCAYAAAAkMCseAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjguMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/H5lhTAAAACXBIWXMAAA9hAAAPYQGoP6dpAAA4IUlEQVR4nO3de3hU5bn+8XtIYAiBhGMOmHAQ5SyigBAqiLKJglLZgIVNq9AqlYqiZAdtsFalP5tW1OIRNy0gFi1eJVBQsAUrAVqgEiSonItIQkg4iGYANYfJ+/sjzciQAAlmDu/M93Ndc8mseVfyPK9ZWXfWrLXGYYwxAgAAsESDQBcAAABQF4QXAABgFcILAACwCuEFAABYhfACAACsQngBAABWIbwAAACrEF4AAIBVIgNdQH2rqKjQkSNH1KxZMzkcjkCXAwAAasEYo1OnTqlt27Zq0ODCx1ZCLrwcOXJEycnJgS4DAABcgvz8fCUlJV1wTMiFl2bNmkmqbD4mJibA1QAAgNpwuVxKTk727McvJOTCS9VbRTExMYQXAAAsU5tTPjhhFwAAWIXwAgAArEJ4AQAAViG8AAAAqxBeAACAVQgvAADAKoQXAABgFcILAACwSsjdpM5XSsvcemXVRh04WqhO8Ym679ZBatQwItBlAQAQdggvtfDwwmV6bteDcjc9XLmgSEr/Z5LSuj+vp388OrDFAQAQZnjb6CIeXrhMsw+NlTv6sNdyd3SBZh8aq4cXLgtQZQAAhCfCywWUlrn13K4HJRnp3I9acBhJ0nO7HlJpmdvvtQEAEK58Gl7mzp2rXr16eT4kMSUlRe++++4F11m/fr369Omjxo0b6/LLL9err77qyxIv6JVVGyvfKjrfZ0Q5jNxN8/XKqo1+rQsAgHDm0/CSlJSk3/zmN8rJyVFOTo5uuukm3X777dq5c2eN4w8ePKgRI0Zo0KBB2r59u2bOnKlp06YpKyvLl2We14GjhfU6DgAAfHc+PWF35MiRXs+feuopzZ07V1u2bFGPHj2qjX/11VfVrl07zZkzR5LUrVs35eTk6JlnntGYMWN8WWqNOsUnSkW1HAcAAPzCb+e8uN1uLVmyRGfOnFFKSkqNYzZv3qzU1FSvZTfffLNycnJUVlZW4zolJSVyuVxej/py362DFHE6STLned/IOBRxOln33Tqo3r4nAAC4MJ+Hl48//lhNmzaV0+nUlClTtHz5cnXv3r3GsUVFRYqPj/daFh8fr/Lycp04caLGdTIzMxUbG+t5JCcn11vtjRpGKK3785VPzg0w/3me1n0O93sBAMCPfB5eunTpotzcXG3ZskU/+9nPNHHiRO3ateu84x0O75BgjKlxeZWMjAwVFxd7Hvn5+fVXvKSnfzxaM9ovVcSZy7yWR5xJ0oz2S7nPCwAAfubzm9Q1atRIV1xxhSSpb9++2rp1q55//nn93//9X7WxCQkJKiryPsnk2LFjioyMVKtWrWr8+k6nU06ns/4LP8vTPx6t/1d2O3fYBQAgCPj9DrvGGJWUlNT4WkpKit5++22vZWvWrFHfvn3VsGFDf5R3Xo0aRuihUUMCWgMAAPDx20YzZ87Uxo0b9dlnn+njjz/Wo48+quzsbP3whz+UVPmWz1133eUZP2XKFB06dEhpaWnavXu3FixYoPnz5ys9Pd2XZQIAAIv49MjL0aNHdeedd6qwsFCxsbHq1auX/vrXv2rYsGGSpMLCQuXl5XnGd+zYUatXr9b06dP18ssvq23btnrhhRcCcpk0AAAITg5TdUZsiHC5XIqNjVVxcbFiYmICXQ4AAKiFuuy/+WwjAABgFcILAACwit+vNkJ4Ky1zc8k5AOA7IbzAbx5euEzP7Xqw8pO6JalISv9nktK6P8/N/gAAtcbbRvCLhxcu0+xDY+WOPuy13B1doNmHxurhhcsCVBkAwDaEF/hcaZlbz+16UJKRzv2UB0flxW7P7XpIpWVuv9cGALAP4QU+98qqjZVvFZ3nw7nlMHI3zdcrqzb6tS4AgJ0IL/C5A0cL63UcACC8EV7gc53iE+t1HAAgvBFe4HP33TpIEaeTJHOe942MQxGnk3XfrYP8WxgAwEqEF/hco4YRSuv+fOWTcwPMf56ndZ/D/V4AALVCeIFfPP3j0ZrRfqkizlzmtTziTJJmtF/KfV4AALXGBzPCr7jDLgCgJnXZf3OHXfhVo4YRemjUkECXAQCwGG8bAQAAqxBeAACAVQgvAADAKoQXAABgFcILAACwCuEFAABYhfACAACsQngBAABWIbwAAACrEF4AAIBVCC8AAMAqhBcAAGAVwgsAALAK4QUAAFiF8AIAAKxCeAEAAFYhvAAAAKsQXgAAgFUILwAAwCo+DS+ZmZnq16+fmjVrpri4OI0aNUp79+694DrZ2dlyOBzVHnv27PFlqQAAwBI+DS/r16/X1KlTtWXLFq1du1bl5eVKTU3VmTNnLrru3r17VVhY6HlceeWVviwVAABYItKXX/yvf/2r1/OFCxcqLi5O27Zt0+DBgy+4blxcnJo3b+7D6gAAgI38es5LcXGxJKlly5YXHXvNNdcoMTFRQ4cO1bp16847rqSkRC6Xy+sBAABCl9/CizFGaWlpuv7669WzZ8/zjktMTNS8efOUlZWlZcuWqUuXLho6dKg2bNhQ4/jMzEzFxsZ6HsnJyb5qAQAABAGHMcb44xtNnTpVq1at0j/+8Q8lJSXVad2RI0fK4XBo5cqV1V4rKSlRSUmJ57nL5VJycrKKi4sVExPznesGAAC+53K5FBsbW6v9t1+OvDzwwANauXKl1q1bV+fgIkkDBgzQ/v37a3zN6XQqJibG6wEAAEKXT0/YNcbogQce0PLly5Wdna2OHTte0tfZvn27EhMT67k6AABgI5+Gl6lTp+rNN9/UihUr1KxZMxUVFUmSYmNjFRUVJUnKyMhQQUGBXn/9dUnSnDlz1KFDB/Xo0UOlpaVavHixsrKylJWV5ctSAQCAJXwaXubOnStJGjJkiNfyhQsXatKkSZKkwsJC5eXleV4rLS1Venq6CgoKFBUVpR49emjVqlUaMWKEL0sFAACW8NsJu/5SlxN+AABAcAi6E3YBAADqC+EFAABYhfACAACsQngBAABWIbwAAACrEF4AAIBVCC8AAMAqhBcAAGAVwgsAALAK4QUAAFiF8AIAAKxCeAEAAFYhvAAAAKsQXgAAgFUILwAAwCqEFwAAYBXCCwAAsArhBQAAWIXwAgAArEJ4AQAAViG8AAAAqxBeAACAVQgvAADAKoQXAABgFcILAACwCuEFAABYhfACAACsQngBAABWIbwAAACrEF4AAIBVCC8AAMAqhBcAAGAVn4aXzMxM9evXT82aNVNcXJxGjRqlvXv3XnS99evXq0+fPmrcuLEuv/xyvfrqq74sEwAAWMSn4WX9+vWaOnWqtmzZorVr16q8vFypqak6c+bMedc5ePCgRowYoUGDBmn79u2aOXOmpk2bpqysLF+WCgAALOEwxhh/fbPjx48rLi5O69ev1+DBg2sc88gjj2jlypXavXu3Z9mUKVO0Y8cObd68+aLfw+VyKTY2VsXFxYqJiam32gEAgO/UZf/t13NeiouLJUktW7Y875jNmzcrNTXVa9nNN9+snJwclZWVVRtfUlIil8vl9QAAAKHLb+HFGKO0tDRdf/316tmz53nHFRUVKT4+3mtZfHy8ysvLdeLEiWrjMzMzFRsb63kkJyfXe+0AACB4+C283H///froo4/0pz/96aJjHQ6H1/Oqd7bOXS5JGRkZKi4u9jzy8/Prp2AAABCUIv3xTR544AGtXLlSGzZsUFJS0gXHJiQkqKioyGvZsWPHFBkZqVatWlUb73Q65XQ667VeAAAQvHx65MUYo/vvv1/Lli3T+++/r44dO150nZSUFK1du9Zr2Zo1a9S3b181bNjQV6UCAABL+DS8TJ06VYsXL9abb76pZs2aqaioSEVFRfr66689YzIyMnTXXXd5nk+ZMkWHDh1SWlqadu/erQULFmj+/PlKT0/3ZakAAMASPg0vc+fOVXFxsYYMGaLExETP46233vKMKSwsVF5enud5x44dtXr1amVnZ6t379761a9+pRdeeEFjxozxZakAAMASfr3Piz9wnxcAAOwTtPd5AQAA+K4ILwAAwCqEFwAAYBXCCwAAsArhBQAAWIXwAgAArEJ4AQAAViG8AAAAqxBeAACAVQgvAADAKoQXAABgFcILAACwCuEFAABYhfACAACsQngBAABWIbwAAACrEF4AAIBVCC8AAMAqhBcAAGAVwgsAALAK4QUAAFiF8AIAAKxCeAEAAFYhvAAAAKsQXgAAgFUILwAAwCqEFwAAYBXCCwAAsArhBQAAWIXwAgAArEJ4AQAAViG8AAAAqxBeAACAVXwaXjZs2KCRI0eqbdu2cjgc+stf/nLB8dnZ2XI4HNUee/bs8WWZAADAIpG+/OJnzpzR1VdfrR//+McaM2ZMrdfbu3evYmJiPM/btGnji/IAAICFfBpehg8fruHDh9d5vbi4ODVv3rz+CwIAANYLynNerrnmGiUmJmro0KFat27dBceWlJTI5XJ5PQAAQOgKqvCSmJioefPmKSsrS8uWLVOXLl00dOhQbdiw4bzrZGZmKjY21vNITk72Y8UAAMDfHMYY45dv5HBo+fLlGjVqVJ3WGzlypBwOh1auXFnj6yUlJSopKfE8d7lcSk5OVnFxsdd5MwAAIHi5XC7FxsbWav8dVEdeajJgwADt37//vK87nU7FxMR4PQAAQOgK+vCyfft2JSYmBroMAAAQJHx6tdHp06f173//2/P84MGDys3NVcuWLdWuXTtlZGSooKBAr7/+uiRpzpw56tChg3r06KHS0lItXrxYWVlZysrK8mWZAADAIj4NLzk5Obrxxhs9z9PS0iRJEydO1GuvvabCwkLl5eV5Xi8tLVV6eroKCgoUFRWlHj16aNWqVRoxYoQvywQAABbx2wm7/lKXE34AAEBwCKkTdgEAAM5GeAEAAFYhvAAAAKsQXgAAgFUILwAAwCqEFwAAYBXCCwAAsArhBQAAWIXwAgAArEJ4AQAAViG8AAAAqxBeAACAVQgvAADAKoQXAABgFcILAACwCuEFAABYhfACAACsQngBAABWIbwAAACrEF4AAIBVCC8AAMAqhBcAAGAVwgsAALAK4QUAAFiF8AIAAKxCeAEAAFYhvAAAAKsQXgAAgFUILwAAwCqEFwAAYBXCCwAAsArhBQAAWMWn4WXDhg0aOXKk2rZtK4fDob/85S8XXWf9+vXq06ePGjdurMsvv1yvvvqqL0sEAACW8Wl4OXPmjK6++mq99NJLtRp/8OBBjRgxQoMGDdL27ds1c+ZMTZs2TVlZWb4sEwAAWCTSl198+PDhGj58eK3Hv/rqq2rXrp3mzJkjSerWrZtycnL0zDPPaMyYMT6qEgAA2CSoznnZvHmzUlNTvZbdfPPNysnJUVlZWYCqAgAAwcSnR17qqqioSPHx8V7L4uPjVV5erhMnTigxMbHaOiUlJSopKfE8d7lcPq8TAAAETlAdeZEkh8Ph9dwYU+PyKpmZmYqNjfU8kpOTfV4jAAAInKAKLwkJCSoqKvJaduzYMUVGRqpVq1Y1rpORkaHi4mLPIz8/3x+lAgCAAAmqt41SUlL09ttvey1bs2aN+vbtq4YNG9a4jtPplNPp9Ed5AAAgCPj0yMvp06eVm5ur3NxcSZWXQufm5iovL09S5VGTu+66yzN+ypQpOnTokNLS0rR7924tWLBA8+fPV3p6ui/LBAAAFvHpkZecnBzdeOONnudpaWmSpIkTJ+q1115TYWGhJ8hIUseOHbV69WpNnz5dL7/8stq2basXXniBy6QBAICHw1SdERsiXC6XYmNjVVxcrJiYmECXAwAAaqEu+++gOmEXAADgYggvAADAKoQXAABgFcILAACwCuEFAABYhfACAACsQngBAABWIbwAAACrEF4AAIBVCC8AAMAqhBcAAGAVwgsAALAK4QUAAFiF8AIAAKxCeAEAAFYhvAAAAKsQXgAAgFUILwAAwCqEFwAAYBXCCwAAsArhBQAAWIXwAgAArEJ4AQAAViG8AAAAqxBeAACAVQgvAADAKoQXAABgFcILAACwCuEFAABYhfACAACsQngBAABWIbwAAACrEF4AAIBV/BJeXnnlFXXs2FGNGzdWnz59tHHjxvOOzc7OlsPhqPbYs2ePP0oFAABBzufh5a233tJDDz2kRx99VNu3b9egQYM0fPhw5eXlXXC9vXv3qrCw0PO48sorfV0qAACwgM/Dy3PPPae7775b99xzj7p166Y5c+YoOTlZc+fOveB6cXFxSkhI8DwiIiJ8XSoAALCAT8NLaWmptm3bptTUVK/lqamp2rRp0wXXveaaa5SYmKihQ4dq3bp15x1XUlIil8vl9QAAAKHLp+HlxIkTcrvdio+P91oeHx+voqKiGtdJTEzUvHnzlJWVpWXLlqlLly4aOnSoNmzYUOP4zMxMxcbGeh7Jycn13gcAAAgekf74Jg6Hw+u5MabasipdunRRly5dPM9TUlKUn5+vZ555RoMHD642PiMjQ2lpaZ7nLpeLAAMAQAjz6ZGX1q1bKyIiotpRlmPHjlU7GnMhAwYM0P79+2t8zel0KiYmxusBAABCl0/DS6NGjdSnTx+tXbvWa/natWs1cODAWn+d7du3KzExsb7LAwAAFvL520ZpaWm688471bdvX6WkpGjevHnKy8vTlClTJFW+7VNQUKDXX39dkjRnzhx16NBBPXr0UGlpqRYvXqysrCxlZWX5ulQAAGABn4eXcePG6fPPP9esWbNUWFionj17avXq1Wrfvr0kqbCw0OueL6WlpUpPT1dBQYGioqLUo0cPrVq1SiNGjPB1qQAAwAIOY4wJdBH1yeVyKTY2VsXFxZz/AgCAJeqy/+azjQAAgFUILwAAwCqEFwAAYBXCCwAAsArhBQAAWIXwAgAArEJ4AQAAViG8AAAAqxBeAACAVQgvAADAKoQXAABgFcILAACwCuEFAABYhfACAACsQngBAABWIbwAAACrEF4AAIBVCC8AAMAqhBcAAGAVwgsAALAK4QUAAFiF8AIAAKxCeAEAAFYhvAAAAKsQXgAAgFUILwAAwCqEFwAAYBXCCwAAsArhBQAAWIXwAgAArEJ4AQAAViG8AAAAqxBeAACAVfwSXl555RV17NhRjRs3Vp8+fbRx48YLjl+/fr369Omjxo0b6/LLL9err77qjzIBAIAFfB5e3nrrLT300EN69NFHtX37dg0aNEjDhw9XXl5ejeMPHjyoESNGaNCgQdq+fbtmzpypadOmKSsry9elAgAACziMMcaX36B///669tprNXfuXM+ybt26adSoUcrMzKw2/pFHHtHKlSu1e/duz7IpU6Zox44d2rx580W/n8vlUmxsrIqLixUTE1M/TQAAAJ+qy/7bp0deSktLtW3bNqWmpnotT01N1aZNm2pcZ/PmzdXG33zzzcrJyVFZWVm18SUlJXK5XF4PAAAQunwaXk6cOCG32634+Hiv5fHx8SoqKqpxnaKiohrHl5eX68SJE9XGZ2ZmKjY21vNITk6uvwYAAEDQ8csJuw6Hw+u5MabasouNr2m5JGVkZKi4uNjzyM/Pr4eKAQBAsIr05Rdv3bq1IiIiqh1lOXbsWLWjK1USEhJqHB8ZGalWrVpVG+90OuV0OuuvaAAAENR8euSlUaNG6tOnj9auXeu1fO3atRo4cGCN66SkpFQbv2bNGvXt21cNGzb0Wa0AAMAOPn/bKC0tTX/4wx+0YMEC7d69W9OnT1deXp6mTJkiqfJtn7vuusszfsqUKTp06JDS0tK0e/duLViwQPPnz1d6erqvSwUAABbw6dtGkjRu3Dh9/vnnmjVrlgoLC9WzZ0+tXr1a7du3lyQVFhZ63fOlY8eOWr16taZPn66XX35Zbdu21QsvvKAxY8b4ulQAAGABn9/nxd+4zwsAAPYJmvu8AAAA1DfCCwAAsArhBQAAWIXwAgAArEJ4AQAAViG8AAAAqxBeAACAVQgvAADAKoQXAABgFcILAACwCuEFAABYhfACAACsQngBAABWIbwAAACrEF4AAIBVCC8AAMAqhBcAAGAVwgsAALAK4QUAAFiF8AIAAKxCeAEAAFaJDHQBAOrO7ZY2bpQKC6XERGnQICkiItBVAYB/EF5CCDu0+hHs87hsmfTgg9Lhw98uS0qSnn9eGj06cHWdK9jnUbKjRgDV8bZRiFi2TOrQQbrxRmnChMr/duhQuRy1F+zzuGyZNHasd3CRpIKCyuXBVGcwz6NkR40AauYwxphAF1GfXC6XYmNjVVxcrJiYmECX4xdVO7Rz/086HJX/Xbo0uP4iD1bBPo9ud+XO9dzgUsXhqDwCc/BgYI8eBPs8SnbUCISbuuy/CS+Ws2WHFuxsmMfs7MqjAxezbp00ZIivq6mZDfNoQ41AOKrL/pu3jSy3ceP5fwlLlX9Z5udXjsP52TCPhYX1O84XbJhHG2oEcGGEF8vZsEOzgQ3zmJhYv+N8wYZ5tKFGABdGeLGcDTs0G9gwj4MGVb6dUXVexrkcDik5uXJcoNgwjzbUCODCCC+Ws2GHZgMb5jEiovJy6Kp6zlb1fM6cwJ6nYcM82lAjgAsjvFjOhh2aDWyZx9GjK6+Euewy7+VJScFxhYwN82hDjQAujPASAoJ9h2YLW+Zx9Gjps88qryp6883K/x48GFz1Bfs82lAjgPPz6aXSX3zxhaZNm6aVK1dKkr7//e/rxRdfVPPmzc+7zqRJk7Ro0SKvZf3799eWLVtq9T3D7VLps3G30PrBPNYPG+bRhhqBcBE093kZPny4Dh8+rHnz5kmSfvrTn6pDhw56++23z7vOpEmTdPToUS1cuNCzrFGjRmrZsmWtvmc4hxcAAGxVl/23zz7baPfu3frrX/+qLVu2qH///pKk3//+90pJSdHevXvVpUuX867rdDqVkJDgq9IAAIDFfHbOy+bNmxUbG+sJLpI0YMAAxcbGatOmTRdcNzs7W3FxcercubMmT56sY8eOnXdsSUmJXC6X1wMAAIQun4WXoqIixcXFVVseFxenoqKi8643fPhwvfHGG3r//ff17LPPauvWrbrppptUUlJS4/jMzEzFxsZ6HsnJyfXWAwAACD51Di9PPPGEHA7HBR85OTmSJEcNN1IwxtS4vMq4ceN06623qmfPnho5cqTeffdd7du3T6tWrapxfEZGhoqLiz2P/Pz8urYEAAAsUudzXu6//36NHz/+gmM6dOigjz76SEePHq322vHjxxUfH1/r75eYmKj27dtr//79Nb7udDrldDpr/fUAAIDd6hxeWrdurdatW190XEpKioqLi/XBBx/ouuuukyT961//UnFxsQYOHFjr7/f5558rPz9fidyrGwAAyIfnvHTr1k233HKLJk+erC1btmjLli2aPHmybrvtNq8rjbp27arly5dLkk6fPq309HRt3rxZn332mbKzszVy5Ei1bt1a//3f/+2rUgEAgEV8eofdN954Q1dddZVSU1OVmpqqXr166Y9//KPXmL1796q4uFiSFBERoY8//li33367OnfurIkTJ6pz587avHmzmjVr5stSAQCAJXx6k7pA4CZ1AADYpy77bz7bCAAAWIXwAgAArEJ4AQAAViG8AAAAqxBeAACAVQgvAADAKoQXAABgFcILAACwSp0/2wgA4D/uUrc+fmWjvjpQqCadEnXVfYMU0Sgi0GUBAUV4AYAgteXhZWr33IPq7T7sWXYkPUl5ac9rwNOjA1gZEFi8bQQAQWjLw8t03eyxSjgruEhSgrtA180eqy0PLwtQZUDgEV4AIMi4S91q99yDkky1X9INVPlxdMnPPSR3qdvvtQHBgPACAEHm41c2qq378Hl/QTeQ0WXufH38yka/1gXfcpe6lTsnW5se+JNy52QTTi+Ac14AIMh8daCwXsch+HF+U90QXgAgyDTplFiv43yNK6K+m6rzm/SftwSrJLgLlDB7rLZoKQHmHA5jjLn4MHu4XC7FxsaquLhYMTExgS4HAOrMXerW0SYdlOAu8JzjcrYKOVQYkaSErw4GPCRUHTFoe/YRgwiOGNTWt/+va36bMJj+X/taXfbfnPMCAEEmolGE8tKel1S58zpb1fP8tDkB35nZdEVUsJ5PYtv5TcEyj4QXAAhCA54erQ9mLFVRxGVeywsjkvTBjMC/jWDTFVFbHl6mo006qPf0GzXwpQnqPf1GHW3SISjClU3nNwXTPBJeACBIDXh6tOK/+ky5v1unTfe/qdzfrVPCVwcDHlwke44YBPvRIVvObwq2eeScFwBAnW164E8a+NKEi4+7/00NfPF//FBRdTacT2LD+U3+mkfOeQEA+JQNRwxsODpkw/lNwTiPhBcAQJ1ddd8gHYlIqrbDrVIhhwoiknXVfYP8XNm3bDmfJNjPbwrGeeQ+LwCAOqs6YpAwe6wq5PB6y+PsIwaXBfCIgQ1Hh6oMeHq03P/vduWec7+cQM5flWCcR855AQBcspru81IQkaz8tDkBP2Jgw/kkNvDXPHLOCwDAL4L5iigbziexQTDOI0deAAAhLZiPDtnE1/NYl/034QUAEPL4/KX64ct5JLwQXgAAsArnvAAAgJBFeAEAAFYhvAAAAKsQXgAAgFUILwAAwCo+DS9PPfWUBg4cqCZNmqh58+a1WscYoyeeeEJt27ZVVFSUhgwZop07d/qyTAAAYBGfhpfS0lLdcccd+tnPflbrdZ5++mk999xzeumll7R161YlJCRo2LBhOnXqlA8rBQAAtvBpeHnyySc1ffp0XXXVVbUab4zRnDlz9Oijj2r06NHq2bOnFi1apK+++kpvvvmmL0sFAACWCKpzXg4ePKiioiKlpqZ6ljmdTt1www3atGlTjeuUlJTI5XJ5PQAAQOiKDHQBZysqKpIkxcfHey2Pj4/XoUOHalwnMzNTTz75ZLXlhBgAAOxRtd+uzY3/6xxennjiiRrDwtm2bt2qvn371vVLezgc3p9aaYyptqxKRkaG0tLSPM8LCgrUvXt3JScnX/L3BwAAgXHq1CnFxsZecEydw8v999+v8ePHX3BMhw4d6vplJUkJCQmSKo/AJCYmepYfO3as2tGYKk6nU06n0/O8adOmys/PV7Nmzc4beGzicrmUnJys/Pz8kP2spnDoUQqfPn0tHOaRHkNHuPRZH4wxOnXqlNq2bXvRsXUOL61bt1br1q0vqbCL6dixoxISErR27Vpdc801kiqvWFq/fr1++9vf1uprNGjQQElJST6pL5BiYmJC/gc/HHqUwqdPXwuHeaTH0BEufX5XFzviUsWnJ+zm5eUpNzdXeXl5crvdys3NVW5urk6fPu0Z07VrVy1fvlxS5dtFDz30kH79619r+fLl+uSTTzRp0iQ1adJEEyZM8GWpAADAEj49YfeXv/ylFi1a5HledTRl3bp1GjJkiCRp7969Ki4u9ox5+OGH9fXXX+u+++7TF198of79+2vNmjVq1qyZL0sFAACW8Gl4ee211/Taa69dcMy5ZxU7HA498cQTeuKJJ3xXmEWcTqcef/xxr/N6Qk049CiFT5++Fg7zSI+hI1z69DeHqc01SQAAAEEiqG5SBwAAcDGEFwAAYBXCCwAAsArhBQAAWIXwUk+OHj2qkydPBroMn6JH1Fa4zGM49EmPCEaEl3qwf/9+ff/739f//M//aNeuXYEuxyfoEbUVLvMYDn3SI4IV4aUeXHnllXrmmWfUtWtXDR48OCQ3AHpEbYXLPIZDn/SIYMV9Xr6jiooKNWjwbQZMS0vTzp07tWzZMkVHRwewsvpDj6HRoz+EyzyGQ5/0GBo9hiqOvFyiiooKGWPUoEEDlZWVqaKiQpLUr18/HT9+3PPcZvQYGj36Q7jMYzj0SY+h0WOo8+nHA4Sys9N6w4YNPf/OycnR0aNHq33sgY3oMTR69Idwmcdw6JMeQ6PHUEd4qaOcnBz97W9/U1lZmdxut4wxioiIUEFBgT799FP94x//0MKFC63+6HN6DI0e/SFc5jEc+qTH0OgxXHDOSx2cOXPG8+nWkydPVn5+vtxut+Lj41VRUaFWrVpp0qRJuvrqqz3JvqCgQG63WzExMWrevHkAq68degyNHv0hXOYxHPqkx9DoMawY1MnOnTtNbGysmT179gXHHT9+3IwcOdJ07drVtG/f3vTp08fs2LHDT1V+N/T4LZt79Idwmcdw6JMev2Vzj+GC8HIJtm7dapxOp3nssce8lldUVBhjjNm8ebPp2LGj6d27t1m4cKH529/+ZmbMmGHi4+PNoUOHAlFyndFjaPToD+Eyj+HQJz2GRo/hgPByibZt22YcDof59a9/bYwxpry83BhjzMcff2ySk5PNsGHDzBdffOEZf/LkSXPNNdeYd955JxDlXhJ6DI0e/SFc5jEc+qTH0Ogx1HGp9CW69tprtX37djmdTp06dUoRERGSpAkTJujyyy/XO++8o+bNm6u8vFyS1KJFCx08eFCFhYWBLLtO6DE0evSHcJnHcOiTHkOjx5AX6PRku7KyMs+/X3rpJdOpUyezZ8+eaq89+eSTJjk52Rw5csTvNX5X9BgaPfpDuMxjOPRJj6HRY6jiyMt3FBn57dXmhw8fVosWLdS+fXuv115//XW99957GjdunGJjYyVJX375pVasWKGioiL/F11H9BgaPfpDuMxjOPRJj6HRY6givNQD85+rzaOiopSYmKjGjRt7Xps3b57mz5+vFi1a6L777lOTJk0kSVu2bNGSJUs0duxY7d69OyB11wU9hkaP/hAu8xgOfdJjaPQYkgJ41Cfk7NmzxzRr1sxMnjzZ/P73vzdjx441vXr1MhMmTDD79+83xhjjdrs94zds2GCGDBliHA6HKSoqClTZdUKPodGjP4TLPIZDn/QYGj2GEsJLPcvNzTXDhg0zAwYMMCkpKWbRokWe90nP/sE3xph169aZhIQEc9tttwWi1EtGj6HRoz+EyzyGQ5/0GBo9hgrCiw988803pqSkxGtZ1T0Eqv79/vvvm/j4eHP77bfXOCbY0WNo9OgP4TKP4dAnPYZGj6GA8OJDNf0wu91us27dOhMXF+f1g+92u6384afH0OjRH8JlHsOhT3oMjR5tRnjxsVOnTpns7GzPIcfz/eDbjB4r2d6jP4TLPIZDn/RYyfYebcUHM/rY8ePH9V//9V/q37+/7rjjDv3oRz/SwIEDtXz5ckmVZ7o7HA7P+AULFmjjxo267LLLNGnSJF1xxRWBKr3W6DE0evSHcJnHcOiTHkOjR2sFLDaFka1bt5ro6GjjcDjM6NGjPcvLy8u9DjXOnz/fNGjQwHz/+983t912m4mPjzcfffRRIEquM3oMjR79IVzmMRz6pMfQ6NFG3OfFD/r27asPPvhALVq0UGJioiTJ7XYrIiLCK7Xn5+drwoQJWrp0qd5++23dc889Gj9+vI4fPx6o0muNHkOjR38Il3kMhz7pMTR6tFHkxYegPnTv3l1r1qxRv379dNVVV+mnP/2ppMof+DfeeEP79u1TcXGxoqKi1LBhQ0nSj370I61YsUKFhYVq06ZNIMuvFXoMjR79IVzmMRz6pMfQ6NE2HHnxoz59+ujf//63evXq5UnsjzzyiGbNmqV27dopISFBf/7zn/XrX/9ap0+f1pYtW7Rz584AV1039BgaPfpDuMxjOPRJj6HRo1UC/b5VONuwYYNxOBxm06ZNnmV/+ctfTEREhBkxYoRxOBzm8ccfD1yB9YAeQ6NHfwiXeQyHPukxNHoMZoSXAMrNzTXdunUz27ZtM2VlZcbtdpvi4mLTr18/s23bNpOfnx/oEr8zegyNHv0hXOYxHPqkx9DoMZjxtlEAJSYmKjIyUitWrFBkZKQaNGigTz/9VHv27NHp06eVlJQkqfLkMFvRY2j06A/hMo/h0Cc9hkaPwYz7vARYbm6uhg0bpmHDhsnpdGrTpk3q3r27Fi9erNOnT8sYo4SEBM/Z7VL1ewsEO3oMjR79IVzmMRz6pMfQ6DFYceQlwHr37q0NGzaoVatWOnTokMaPH6+MjAxFR0frjTfeUO/evXXo0CFFRESotLRUkjw/+BUVFVakenoMjR79IVzmMRz6pMfQ6DFo+f+dKtTE7XabsrIyr2VlZWVmypQp5rLLLjOfffaZMcaYr7/+2qxatcpMmDDBDB061IwfP97s3LkzECXXGT2GRo/+EC7zGA590mNo9BhsuM9LkGjQoIEaNPA+EBYZGalXXnlFTZs21datW9W6dWtlZGToH//4hxo0aKDbbrtN+/bt08CBA/Wvf/1LXbp0CVD1tUOPodGjP4TLPIZDn/QYGj0GnUCnJ9ROSUmJmTJliunXr5+ZNWuW122px40bZ9LT063/gDB6DI0e/SFc5jEc+qTH0OjR3zjyYgG3263nnntOH374ocaNG6dp06bJ4XDI/Odc6xMnTqhp06bVkr9N6DE0evSHcJnHcOiTHkOjx0Bgtizw9ddf66233tLVV1+tqVOnqmHDhnK73XI4HDp27JhOnTqldu3aqaKiQpJUVlYW4Irrjh5Do0d/qM08Jicnyxhj9QmR4dAnPYZGjwERoCM+qIPly5ebiIgI8+WXXxpjjNeJYT/5yU9Mq1atTFFRkTHGmBMnTpjJkydb92mm9BgaPfpDbeaxsLDQa52zD9PbIhz6pMfQ6DEQOPJigdatW6t79+46c+aMpMoTwSoqKjR16lQtWbJEixcvVnx8vNxut1q1aqXk5GT17t1bO3bsCHDltUePodGjP9RmHhMSElRUVKSlS5dq0qRJmjVrlvbv3x/gyusmHPqkx9DoMSACnZ5wcXl5eSYxMdHMnDnTfPLJJ2bVqlVm1KhRpmnTpmbp0qU1pvSZM2eaJk2amO3bt/u/4EtAj6HRoz+cbx6jo6PN0qVLTXl5uTGm8i/ejh07mpSUFPODH/zAtG7d2qojWeHQJz2GRo+BQHixxLZt20znzp1Nhw4dTIcOHcy1115rPvzww2rjzj5jffbs2SY6Otrs3r3bn6VeMnr8ls09+sO589i7d2+Tk5PjNeYXv/iFGTp0qOf5Y489Znr16mWOHj3q73IvWTj0SY+VbO/R3/h4AIscP35chYWFcjqdateunaKiojy3mjbn3HL6+PHj2rVrl2688UbFx8drx44diouLC2D1tUOPodGjP5w9j0lJSYqOjtaOHTuUnZ2tf/7zn4qOjlaDBg00f/58SdLWrVv1ox/9SH/+85/Vq1evAFdfe+HQJz2GRo/+xKXSFmnTpo3atGnjeV5eXu55/7RBgwYqLy9XeXm5nnnmGa1bt07Z2dm68847de2111qzw6PH0OjRH86dxyNHjmjWrFn68ssvdfnllys+Pl7PPPOMrrzySk2cOFEfffSRXC6XTp06FcCq6y4c+qTH0OjRrwJ30Af1pby83Hz66afmvvvuMz169DCxsbFmypQp5m9/+1u1cbaiR+9xqC4nJ8e0adPGvP76655l77zzjmnSpIm54YYbTLdu3cyPf/zjAFZYP8KhT3oMjR59iSMvlnv99de1bds2/eEPf1DPnj1100036bHHHlPz5s3VsGFDr7FVn2pqG3oMjR59rbS0VC1atFDfvn0lVX5671VXXaWOHTvqqaeeUlJSktq3by9JXp/ya5tw6JMeQ6NHX+JSaYsdP35cM2bM0I4dO/SrX/1Kf//73zVnzhy1adPGa4dXddMzG9FjaPToD506dZLT6dSiRYskVX56b8OGDVVcXKyvvvrKsyMoKyuzekcQDn3SY2j06EucsGu5kydPqqysTPHx8Z5l5j8nfebl5Sk6OlqtWrWqltzNOSeGBrNL7fHsccHuu/TIX2Xfys3N1bBhw3TzzTcrOjpa2dnZatKkiXJycnT8+HE1bNjQ+u1BCo8+L7VHyZ4+v0uP4b7dc+TFci1btvTs8KpyaNVG+/7776tTp046dOiQIiIiVFpa6rVuRUWFFbegr0uPZ99eu6ioSO+9956Kior8X3Qd1aXHs4/ALF++XLNnz9ahQ4f8X3QQ6t27tzZu3KioqCgdOXJEKSkp2rBhgyIiIrRq1aqQ2B6k8OizLj3aut3XpUe2+3P4/Swb+E1FRYW59957TfPmzc3hw4c9y5csWWJ+9rOfme9973vmzjvvNDt37gxgld9NVY8tWrQwBQUFXq9t2LDBjB071gwYMCBke3zkkUfMmDFjTKdOnbih1VlqOqnZ7XaH3PYQDn1eqMdQ2e7r2iPbPTepC3lut9s89thjZvny5cYYYzIyMkxUVJS54447zNSpU83dd99tmjdvbvbs2RPYQr8Dt9ttnnjiCbNs2TKv5RUVFebAgQPme9/7nmnZsqU5fvx4gCr87s7t8ey78Z48edL8/Oc/N+3atTMHDx4MUIXB6+wb/oXy9hAOfZ7bYyhu9xfqke3+W4SXEHf2D/vDDz9snE6n+cMf/mBOnjzpWT5y5Ejz85//PBDl1auSkhLPv0tLS40xxvz97383jRs3Nj/5yU9Mfn5+oEqrN2f3eO7yIUOGmJdfftnPFdklXLaHcOnTGLb7cN3uuVQ6xFWdNzF79mw9++yz+uMf/6jx48dXO5nts88+C0B19atRo0YyxqiiokINGzbUe++9p+HDh3s+6CwxMTHQJV6SqpvXGWPUqFEjSd/e2E6SZ7nL5dKnn34ayFKDXrhsD+HSp8R2H7bbfSCTE/xjz549plu3buapp54y33zzjddrubm5plOnTmbevHkBqq5+VR1yfe+990xkZKS55557zJEjRzyvn/1x9LYpKioyZ86cMadOnTIul8ucOXPGfPXVV2b37t1m9uzZpkGDBmbNmjWBLjPohcv2EC59GsN2H47bPUdewkBBQYEKCwuVmpoqp9PpWf7ZZ5/ptddeU+vWrdWjR48AVvjdfPPNN/rkk0/Us2dPNW7cWO+//75uueWWan95ud1uz18tf/7zn3Xdddd57qUQ7ObNm6dp06apc+fOcrlcioqK0jfffKPS0lI1btxYp06d0ty5c3XDDTdIqux17dq16t27txISEgJcfXCpzfZw9dVXe5bv27dPcXFxat68eQCqvXRs92z3oYzwEgbKysrUvn17XXnllZ5le/fu1YIFC7R8+XLNmDFDAwcOlCQ9/fTTKigoUFRUlO655x5dccUVgSq71r7++muNHj1aY8eO1S233KJbb71VkyZN0pNPPun1C6zqngiLFi3Sb37zGxUWFmr9+vVeO6pg1b9/f0VGRqpFixZatGiRIiMjderUKVVUVOiyyy5TdHS01+ceLViwQHPnzlVkZKQWLFignj17BrD64HKx7eGRRx7RddddJ0navXu3Fi5cqP379+u3v/2tOnfuHKiy64ztnu0+pLf7QB/6ge8dPnzYJCQkmPT0dLNr1y6zYsUKM3bsWNO1a1fz7LPPesb99re/NV27djWpqalmzJgxpk2bNtZchvevf/3LxMbGGofDUe2Q8dkWLFhgBg8ebCZPnmwmTpxooqOjTW5urp+rvTQ7duwwjRo1Mr/85S89y869MuFse/bsMf/7v/9rWrdubd1VJb50oe3hd7/7XbXxK1euNBMmTDCNGjUyBw4c8H/Bl4jt/lts96GH8BImPvzwQ9OpUyfTo0cPExkZaX7wgx+YxYsXe43JzMw0o0aNMp9//rkxxpgZM2aYHj16WHOp4a5du0zz5s3ND3/4wxrvm5Cbm2smTJhgWrRo4dmo09PTTePGja35ZZ2Tk2Oio6PNI488ct4x5/4yu/vuu83EiROtft+/vtW0PSxZssTzenl5udc83nXXXcbhcJj3338/EOVeMrZ7tvtQ3e4JL2Hk2LFjZt++fWbfvn3mzJkznuXp6enmwQcfNGPGjPH6FNMDBw6Yrl27mg8++CAQ5V6Sbdu2maZNm5oXX3yx2i+yb775xmzdutUMHz7cDBs2zLNRP/7442bGjBmBKPeSfPDBB6Z58+bm73//u9fy0tJSzy+wqktGjam8VHbYsGF+rdEGZ28PZ1+Keu7Pzb333muio6PN6tWr/V1ivWC7Z7sPRYSXMHfDDTeYpKQkM2PGDDN9+nQTGRlpnnzySeN2u82sWbNMmzZtrLlLZZV9+/aZtWvXGmOq/zVijDFHjx41gwcPNk899ZS/S6s3R44cMYcOHbrgmIqKCnP8+HGTkpJihg8f7qfK7Hb2/VGM+Ta4vPPOOwGqyDfY7u3Edv8tTtgNY2vXrtWHH36o1atX6/rrr5ckfe9739Mdd9yhTZs2KTc3V/fee6/atm0b4Err5sorr/ScpPj444+rf//+uu222yRV3juhTZs2iouL04EDB7zWM5Z8mJskzwmJJ0+e1PTp0+V2u3X69GkZY9SgQQNVVFSosLBQJ0+elDFGS5YskWRXj/525MgRxcXFea5Muffee/XGG2/orbfe0q233hrg6uoP2z3bfSggvISxVq1aKT4+3nMDJEkaOnSounfvrl/84hfq3LmzmjVrpqioqABWeek+//xzvfvuuzpz5oxuu+02zwYuSSdOnKj24XQ2btynTp3SH//4RzkcDj377LM6fvy4SktL1bJlSzmdTrVp00bjxo3z/D+2sUd/+Oabb/TUU09Jkn73u99p2rRpeuONN7RkyZKQCi4S2z3bfWggvISxyy67TE2bNtWqVas8l4bu2LFDpaWlatKkidcleDZq1aqVFi5cqOuvv17R0dGaOHGirrjiCmVkZGj//v1asGCBJLv/Kmnfvr0++eQTXXfddTLGeHbA56q6Wydq1rhxY919993q16+f3n33XX3xxRchd8SlCts9230oILyEsfj4eC1cuFA33XSTDh8+rNOnT+vLL79UmzZt1KZNm0CXVy+uuuoqrV+/XpMmTdKKFSt08OBBOZ1Opaenq2/fvpLs/6uke/fuys7O1uDBg/XFF19o1qxZ1caE6i+w+nTttddqx44dGjJkiK677jrdcsstgS7JJ9ju2e5DgcMYYwJdBAJr165devnll3XgwAH16tVL48eP17XXXhvosurVkSNHtHfvXu3atUs333yzEhIS1LRp00CXVa8+/PBD9e3bVy+++KKmTp0a6HKstWPHDl1zzTXKyMjQr371q5DdAbDdh4Zw3e4JL5D07Qd+nf3BX7DPrl27dOLECQ0aNMj6vywD6eOPP9aKFSs0bdo0xcTEBLocn2G7Dw3huN0TXuDF5veBgfoUTjt0tnvYhvACAACsEppv5gIAgJBFeAEAAFYhvAAAAKsQXgAAgFUILwAAwCqEFwAAYBXCCwAAsArhBQAAWIXwAgAArEJ4AQAAViG8AAAAq/x/0DbJkIHviEUAAAAASUVORK5CYII=",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Run PCA\n",
    "pca = decomposition.PCA(n_components=1)\n",
    "pca_result = pca.fit_transform(embeddings)\n",
    "\n",
    "plt.xticks(rotation=-45)\n",
    "\n",
    "# Plot all points in blue first\n",
    "plt.scatter(stack.time, pca_result, color=\"blue\")\n",
    "\n",
    "# Re-plot cloudy images in green\n",
    "plt.scatter(stack.time[0], pca_result[0], color=\"green\")\n",
    "plt.scatter(stack.time[2], pca_result[2], color=\"green\")\n",
    "\n",
    "# Color all images after fire in red\n",
    "plt.scatter(stack.time[-5:], pca_result[-5:], color=\"red\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b38b70a6-2156-41f8-967e-a490cc8e2778",
   "metadata": {},
   "source": [
    "### And finally, some finetuning\n",
    "\n",
    "We are going to train a classifier head on the embeddings and use it to detect fires."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "id": "1da07de0-b8f2-46c9-bd2a-58b15ca2224f",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Matched 5 out of 5 correctly\n"
     ]
    }
   ],
   "source": [
    "# Label the images we downloaded\n",
    "# 0 = Cloud\n",
    "# 1 = Forest\n",
    "# 2 = Fire\n",
    "labels = np.array([0, 1, 0, 1, 1, 1, 1, 2, 2, 2, 2, 2])\n",
    "\n",
    "# Split into fit and test manually, ensuring we have all 3 classes in both sets\n",
    "fit = [0, 1, 3, 4, 7, 8, 9]\n",
    "test = [2, 5, 6, 10, 11]\n",
    "\n",
    "# Train a support vector machine model\n",
    "clf = svm.SVC()\n",
    "clf.fit(embeddings[fit] + 100, labels[fit])\n",
    "\n",
    "# Predict classes on test set\n",
    "prediction = clf.predict(embeddings[test] + 100)\n",
    "\n",
    "# Perfect match for SVM\n",
    "match = np.sum(labels[test] == prediction)\n",
    "print(f\"Matched {match} out of {len(test)} correctly\")"
   ]
  }
 ],
 "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.8"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}