{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "d5cf2b3d-0f34-4d13-9ba3-dd2e3fb66337",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Cluster analysis for 2D diurnal histogram of rainfall\n",
    "\n",
    "from docx import Document\n",
    "import pandas as pd\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import matplotlib.ticker as mticker\n",
    "from itertools import chain\n",
    "from matplotlib.colors import LinearSegmentedColormap, Normalize\n",
    "import matplotlib.colors as mcolors\n",
    "from matplotlib.colors import BoundaryNorm\n",
    "\n",
    "from sklearn.utils import resample\n",
    "from collections import defaultdict\n",
    "from scipy.cluster.hierarchy import linkage, fcluster, dendrogram\n",
    "from sklearn.metrics import davies_bouldin_score, silhouette_score\n",
    "\n",
    "\n",
    "# Load the Word document\n",
    "file_path = \"/Users/cdz/Desktop/Khafid_etal/BMKG_SurfaceStationData/CombinedStations31.docx\"\n",
    "doc = Document(file_path)\n",
    "\n",
    "len(doc.tables)\n",
    "\n",
    "table = doc.tables[0]\n",
    "\n",
    "data = []\n",
    "for row in table.rows:\n",
    "    data.append([cell.text.strip() for cell in row.cells])\n",
    "\n",
    "df = pd.DataFrame(data)\n",
    "df.head()\n",
    "\n",
    "# Define column names\n",
    "df.columns = [\n",
    "    \"station_id\",\n",
    "    \"station_name\",\n",
    "    \"region\",\n",
    "    \"longitude\",\n",
    "    \"latitude\",\n",
    "    \"elevation_m\",\n",
    "    \"value\"\n",
    "]\n",
    "\n",
    "id = pd.to_numeric(df[\"station_id\"], errors=\"coerce\")\n",
    "nid = id.shape[0]\n",
    "\n",
    "y_labs = ['Time', 'T (˚C)', 'p (hPa)', 'V (m/s)', 'q (g/kg)', 'WDA (degree)', 'Rain (mm/3hrs)'] \n",
    "x_labs = \"Local Time\"\n",
    "labs = ['(a)', '(b)', '(c)', '(d)', '(e)', '(f)', '(g)', '(h)']\n",
    "\n",
    "cmap = LinearSegmentedColormap.from_list(\"white_to_dark\", [\"white\", \"black\"])\n",
    "norm = BoundaryNorm([0, 1, np.inf], ncolors=cmap.N, clip=True)\n",
    "\n",
    "count_std = np.zeros(6) \n",
    "\n",
    "y1 = [2017, 20, 1000,  0, 13,   0,  0]\n",
    "y2 = [2023, 25, 1021, 10, 25, 360, 30]\n",
    "\n",
    "#var = ['Time', 'rainfall3h', 'Temp', 'Specific humidity', 'Pres', 'Windspeed', 'Windir', 'RH'] \n",
    "vsh = ['T', 'p', 'V', 'q', 'WDA', 'R3h']\n",
    "v = vsh[5]\n",
    "#print(v)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "45b4897f-49d7-4b46-b423-6043b7eced58",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(31, 8, 24)\n"
     ]
    }
   ],
   "source": [
    "# Make a combined PDF including all stations\n",
    "\n",
    "PR_list = []\n",
    "\n",
    "# loop through stations\n",
    "for n in range(nid): \n",
    "#    for n in range(2): \n",
    "    sid = id[n]\n",
    "    \n",
    "    data_in_path = f\"/Users/cdz/Desktop/Khafid_etal/Manuscript/2Dhist_stations_Data/hist2d_data_{sid}_{v}.npz\"\n",
    "#        print(data_in_path)\n",
    "    \n",
    "    data = np.load(data_in_path)\n",
    "    HR = data['H']\n",
    "    PR = 8*100*HR/HR.sum()\n",
    "\n",
    "    PR_list.append(PR)\n",
    "    \n",
    "PT = np.stack(PR_list, axis=0)\n",
    "\n",
    "print(PT.shape)\n",
    "\n",
    "xedges = data[\"xedges\"]\n",
    "yedges = data[\"yedges\"]\n",
    "\n",
    "data_out_path = f\"/Users/cdz/Desktop/Khafid_etal/Manuscript/2Dpdf_data_AllStationsStacked_{v}.npz\"\n",
    "\n",
    "np.savez(data_out_path, H=PT, xedges=xedges, yedges=yedges)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "97aea6d0-29ec-4466-9c05-5c2f7dae6ae4",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Flatten data\n",
    "\n",
    "linkage_method = 'complete'  # \"average\" and \"ward\" produce the same result\n",
    "\n",
    "N = PT.shape[0]\n",
    "X = PT.reshape(N, -1)\n",
    "\n",
    "# Generate the linkage matrix for the entire dataset\n",
    "Z = linkage(X, method=linkage_method) \n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "id": "ffe99be5-daad-4bb0-a337-fc7e447f5f86",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "k     | Stability  | Silhouette   | Davies-Bouldin \n",
      "-------------------------------------------------------\n",
      "2     | 1.000      | 0.292        | 1.148          \n",
      "3     | 1.000      | 0.266        | 1.112          \n",
      "4     | 1.000      | 0.177        | 1.369          \n",
      "5     | 1.000      | 0.152        | 1.267          \n",
      "6     | 1.000      | 0.144        | 1.261          \n",
      "7     | 0.993      | 0.134        | 1.068          \n",
      "8     | 0.964      | 0.133        | 0.995          \n",
      "9     | 0.924      | 0.130        | 0.967          \n",
      "10    | 0.897      | 0.141        | 0.954          \n",
      "11    | 0.943      | 0.165        | 0.935          \n",
      "12    | 0.943      | 0.164        | 0.851          \n",
      "13    | 0.914      | 0.152        | 0.792          \n",
      "14    | 0.910      | 0.148        | 0.779          \n",
      "15    | 0.957      | 0.142        | 0.741          \n"
     ]
    }
   ],
   "source": [
    "# Cluster analysis\n",
    "\n",
    "n_bootstrap = 100 # not change for larger n\n",
    "\n",
    "\n",
    "k_list = [2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15]\n",
    "results = {}\n",
    "\n",
    "for k in k_list:\n",
    "    n_clusters = k\n",
    "\n",
    "    # --- reference clustering for this k ---\n",
    "    Z = linkage(X, method=linkage_method)\n",
    "    labels_ref = fcluster(Z, n_clusters, criterion='maxclust')\n",
    "\n",
    "    \n",
    "    sil_avg = silhouette_score(X, labels_ref) # Calculate silhouette score \n",
    "    db_index = davies_bouldin_score(X, labels_ref) # Calculate Davis_Bouldin \n",
    "    \n",
    "    # --- bootstrap co-occurrence for this k (same logic as before) ---\n",
    "    co_occurrence = np.zeros((N, N))\n",
    "    counts = np.zeros((N, N))\n",
    "\n",
    "    for b in range(n_bootstrap):\n",
    "        indices = np.array(resample(range(N), replace=True, n_samples=N))\n",
    "        X_boot = X[indices]\n",
    "        Z_boot = linkage(X_boot, method=linkage_method)\n",
    "        labels_boot = fcluster(Z_boot, n_clusters, criterion='maxclust')\n",
    "\n",
    "        for c in range(1, n_clusters + 1):\n",
    "            mask_boot = (labels_boot == c)\n",
    "            if mask_boot.sum() < 2:\n",
    "                continue\n",
    "\n",
    "            orig_ids = indices[mask_boot]\n",
    "            ii = np.repeat(orig_ids[:, None], len(orig_ids), axis=1)\n",
    "            jj = np.repeat(orig_ids[None, :], len(orig_ids), axis=0)\n",
    "            counts[ii, jj] += 1.0\n",
    "            co_occurrence[ii, jj] += 1.0\n",
    "\n",
    "    stability_matrix = np.divide(\n",
    "        co_occurrence, counts,\n",
    "        out=np.zeros_like(co_occurrence),\n",
    "        where=counts > 0\n",
    "    )\n",
    "\n",
    "    # cluster-level scores for this k\n",
    "    cluster_stability = defaultdict(list)\n",
    "    for i in range(N):\n",
    "        for j in range(i+1, N):\n",
    "            if labels_ref[i] == labels_ref[j]:\n",
    "                cluster_stability[labels_ref[i]].append(stability_matrix[i, j])\n",
    "\n",
    "    cluster_scores = {\n",
    "        c: np.mean(vals) if len(vals) > 0 else 0.0\n",
    "        for c, vals in cluster_stability.items()\n",
    "    }\n",
    "\n",
    "    # store for later comparison\n",
    "    results[k] = {\n",
    "        \"labels_ref\": labels_ref,\n",
    "        \"cluster_scores\": cluster_scores, \n",
    "        \"mean_stability\": np.mean(list(cluster_scores.values())), \n",
    "        \"silhouette_score\": sil_avg,\n",
    "        \"db_index\": db_index  # Added to results\n",
    "    }\n",
    "\n",
    "\n",
    "print(f\"{'k':<5} | {'Stability':<10} | {'Silhouette':<12} | {'Davies-Bouldin':<15}\")\n",
    "print(\"-\" * 55)\n",
    "for k in k_list:\n",
    "    print(f\"{k:<5} | {results[k]['mean_stability']:<10.3f} | \"\n",
    "          f\"{results[k]['silhouette_score']:<12.3f} | {results[k]['db_index']:<15.3f}\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "5dcc736a-45f2-4466-811f-4dfb5e974061",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA+4AAAJiCAYAAABZ4ZipAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjguMCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy81sbWrAAAACXBIWXMAAA9hAAAPYQGoP6dpAABWI0lEQVR4nO3deXQUVd7G8ac6ayckgQRICIRNo4CgbIosShhkU0QHEBBkUBlFQRwERDIuBEdhREUUFIdFYWQRN9wXQFlEEBFwARFcEBAImxBAmgTIff9g6JcmrEklVUl/P+f0Oemq23V/le6u7qduVbVljDECAAAAAACu5HG6AAAAAAAAcHoEdwAAAAAAXIzgDgAAAACAixHcAQAAAABwMYI7AAAAAAAuRnAHAAAAAMDFCO4AAAAAALhYqNMFuEVubq62bt2qmJgYWZbldDkAAAAAgBLOGKP9+/crOTlZHs/px9UJ7v+zdetWpaSkOF0GAAAAACDIbN68WZUqVTrtfIL7/8TExEg69g+LjY11uBoAAAAAQEm3b98+paSk+PPo6RDc/+f44fGxsbEEdwAAAABAkTnb6dpcnA4AAAAAABcjuAMAAAAA4GIEdwAAAAAAXIzgDgAAAACAixHcAQAAAABwMYI7AAAAAAAuRnAHAAAAAMDFCO4AAAAAALgYwR0AAAAAABcjuAMAAAAA4GIEdwAAAAAAXIzgDgAAAACAixHcAQAAAABwMYI7AAAAAAAuRnAHAAAAAMDFCO4AAAAAALgYwR0AAAAAABcjuAMAAAAA4GIEdwAAAAAAXIzgDgAAAACAixHcAQAAAABwsVCnC0DJYoyR7/BRp8sA4ALesBBZluV0GQAAAMUewR22Mcao84tLtWLjHqdLAeACDauU0et3NSa8AwAAFBCHysM2vsNHCe0A/L7euIcjcAAAAGzAiDsKxdcPXaOo8BCnywDggIM5R9XwsXlOlwEAAFBiOB7cFy1apCeffFIrVqzQtm3bNHv2bN14442nbNunTx9NmDBBzzzzjAYMGOCfnp2drcGDB2vmzJny+Xxq2bKlXnjhBVWqVKloVgJ5RIWHKCrc8ZcXAAAAABR7jh8q/+eff+qyyy7TuHHjztju7bff1rJly5ScnJxn3oABAzR79my9+uqrWrx4sQ4cOKD27dvr6FEO0QQAAAAAFG+OD4m2a9dO7dq1O2ObLVu26J577tEnn3yi6667LmBeVlaWJk+erFdeeUXXXHONJGnatGlKSUnRvHnz1KZNm0KrHQAAAACAwub4iPvZ5ObmqmfPnrr//vt1ySWX5Jm/YsUKHT58WK1bt/ZPS05OVu3atbVkyZLTLjc7O1v79u0LuAEAAAAA4DauD+5PPPGEQkNDde+9955yfmZmpsLDw1WmTJmA6YmJicrMzDztckeOHKm4uDj/LSUlxda6AQAAAACwg6uD+4oVK/Tss89qypQp5/07wMaYMz4mPT1dWVlZ/tvmzZsLWi4AAAAAALZzdXD//PPPtWPHDlWuXFmhoaEKDQ3Vxo0bNWjQIFWtWlWSlJSUpJycHO3ZE/j74Tt27FBiYuJplx0REaHY2NiAGwAAAAAAbuPq4N6zZ0999913+uabb/y35ORk3X///frkk08kSQ0aNFBYWJjmzp3rf9y2bdu0evVqNWnSxKnSAQAAAACwheNXlT9w4IB+/vln//0NGzbom2++UXx8vCpXrqyEhISA9mFhYUpKStLFF18sSYqLi1Pv3r01aNAgJSQkKD4+XoMHD1adOnX8V5kHAAAAAKC4cjy4f/3112rRooX//sCBAyVJvXr10pQpU85pGc8884xCQ0PVpUsX+Xw+tWzZUlOmTFFISEhhlAwAAAAAQJFxPLinpaXJGHPO7X/77bc80yIjIzV27FiNHTvWxsoAAAAAAHCeq89xBwAAAAAg2BHcAQAAAABwMYI7AAAAAAAuRnAHAAAAAMDFCO4AAAAAALgYwR0AAAAAABcjuAMAAAAA4GIEdwAAAAAAXIzgDgAAAACAixHcAQAAAABwMYI7AAAAAAAuRnAHAAAAAMDFCO4AAAAAALgYwR0AAAAAABcjuAMAAAAA4GIEdwAAAAAAXIzgDgAAAACAixHcAQAAAABwMYI7AAAAAAAuRnAHAAAAAMDFCO4AAAAAALgYwR0AAAAAABcjuAMAAAAA4GIEdwAAAAAAXIzgDgAAAACAixHcAQAAAABwMYI7AAAAAAAuRnAHAAAAAMDFCO4AAAAAALgYwR0AAAAAABcjuAMAAAAA4GIEdwAAAAAAXIzgDgAAAACAixHcAQAAAABwMYI7AAAAAAAuRnAHAAAAAMDFCO4AAAAAALgYwR0AAAAAABcjuAMAAAAA4GIEdwAAAAAAXIzgDgAAAACAixHcAQAAAABwMYI7AAAAAAAuRnAHAAAAAMDFCO4AAAAAALgYwR0AAAAAABcjuAMAAAAA4GIEdwAAAAAAXIzgDgAAAACAixHcAQAAAABwMYI7AAAAAAAuRnAHAAAAAMDFCO4AAAAAALgYwR0AAAAAABcjuAMAAAAA4GIEdwAAAAAAXIzgDgAAAACAixHcAQAAAABwMYI7AAAAAAAuRnAHAAAAAMDFCO4AAAAAALiY48F90aJFuv7665WcnCzLsvT222/75x0+fFgPPPCA6tSpo+joaCUnJ+tvf/ubtm7dGrCM7Oxs9e/fX2XLllV0dLQ6dOig33//vYjXBAAAAAAA+zke3P/8809ddtllGjduXJ55Bw8e1MqVK/Xwww9r5cqVeuutt7R+/Xp16NAhoN2AAQM0e/Zsvfrqq1q8eLEOHDig9u3b6+jRo0W1GgAAAAAAFIpQpwto166d2rVrd8p5cXFxmjt3bsC0sWPH6oorrtCmTZtUuXJlZWVlafLkyXrllVd0zTXXSJKmTZumlJQUzZs3T23atCn0dQAAAAAAoLA4PuJ+vrKysmRZlkqXLi1JWrFihQ4fPqzWrVv72yQnJ6t27dpasmTJaZeTnZ2tffv2BdwAAAAAAHCbYhXcDx06pKFDh6p79+6KjY2VJGVmZio8PFxlypQJaJuYmKjMzMzTLmvkyJGKi4vz31JSUgq1dgAAAAAA8qPYBPfDhw+rW7duys3N1QsvvHDW9sYYWZZ12vnp6enKysry3zZv3mxnuQAAAAAA2KJYBPfDhw+rS5cu2rBhg+bOnesfbZekpKQk5eTkaM+ePQGP2bFjhxITE0+7zIiICMXGxgbcAAAAAABwG9cH9+Oh/aefftK8efOUkJAQML9BgwYKCwsLuIjdtm3btHr1ajVp0qSoywUAAAAAwFaOX1X+wIED+vnnn/33N2zYoG+++Ubx8fFKTk5W586dtXLlSr3//vs6evSo/7z1+Ph4hYeHKy4uTr1799agQYOUkJCg+Ph4DR48WHXq1PFfZR4AAAAAgOLK8eD+9ddfq0WLFv77AwcOlCT16tVLGRkZevfddyVJdevWDXjc/PnzlZaWJkl65plnFBoaqi5dusjn86lly5aaMmWKQkJCimQdAAAAAAAoLI4H97S0NBljTjv/TPOOi4yM1NixYzV27Fg7SwMAAAAAwHGuP8cdAAAAAIBgRnAHAAAAAMDFCO4AAAAAALgYwR0AAAAAABcjuAMAAAAA4GIEdwAAAAAAXIzgDgAAAACAixHcAQAAAABwMYI7AAAAAAAuRnAHAAAAAMDFCO4AAAAAALgYwR0AAAAAABcjuAMAAAAA4GIEdwAAAAAAXIzgDgAAAACAixHcAQAAAABwMYI7AAAAAAAuRnAHAAAAAMDFCO4AAAAAALhYqNMFACg5jDHyHT7qdBlw2MGcI6f8G8HLGxYiy7KcLgMAgGKL4A7AFsYYdX5xqVZs3ON0KXCRho996nQJcIGGVcro9bsaE94BAMgnDpUHYAvf4aOEdgCn9PXGPRyNAwBAATDiDsB2Xz90jaLCQ5wuA4DDDuYcVcPH5jldBgAAxR7BHYDtosJDFBXO5gUAAACwA4fKAwAAAADgYgR3AAAAAABcjOAOAAAAAICLEdwBAAAAAHAxgjsAAAAAAC5GcAcAAAAAwMUI7gAAAAAAuBjBHQAAAAAAFyO4AwAAAADgYgR3AAAAAABcjOAOAAAAAICLEdwBAAAAAHAxgjsAAAAAAC5GcAcAAAAAwMUI7gAAAAAAuBjBHQAAAAAAFyO4AwAAAADgYgR3AAAAAABcjOAOAAAAAICLEdwBAAAAAHAxgjsAAAAAAC5GcAcAAAAAwMUI7gAAAAAAuBjBHQAAAAAAFyO4AwAAAADgYgR3AAAAAABcjOAOAAAAAICLEdwBAAAAAHAxgjsAAAAAAC5GcAcAAAAAwMUI7gAAAAAAuBjBHQAAAAAAFyO4AwAAAADgYgR3AAAAAABcjOAOAAAAAICLEdwBAAAAAHAxgjsAAAAAAC5GcAcAAAAAwMUcD+6LFi3S9ddfr+TkZFmWpbfffjtgvjFGGRkZSk5OltfrVVpamtasWRPQJjs7W/3791fZsmUVHR2tDh066Pfffy/CtQAAAAAAoHA4Htz//PNPXXbZZRo3btwp548aNUqjR4/WuHHjtHz5ciUlJalVq1bav3+/v82AAQM0e/Zsvfrqq1q8eLEOHDig9u3b6+jRo0W1GgAAAAAAFIpQpwto166d2rVrd8p5xhiNGTNGDz74oDp27ChJmjp1qhITEzVjxgz16dNHWVlZmjx5sl555RVdc801kqRp06YpJSVF8+bNU5s2bYpsXQAAAAAAsJvjI+5nsmHDBmVmZqp169b+aREREWrevLmWLFkiSVqxYoUOHz4c0CY5OVm1a9f2tzmV7Oxs7du3L+AGAAAAAIDbuDq4Z2ZmSpISExMDpicmJvrnZWZmKjw8XGXKlDltm1MZOXKk4uLi/LeUlBSbqwcAAAAAoOBcHdyPsywr4L4xJs+0k52tTXp6urKysvy3zZs321IrAAAAAAB2cvwc9zNJSkqSdGxUvUKFCv7pO3bs8I/CJyUlKScnR3v27AkYdd+xY4eaNGly2mVHREQoIiKikCoHgKJljJHx+ZwuAwiQm/P/F4nNPehT7pEQB6sBTs3yes86IAQATnN1cK9WrZqSkpI0d+5c1atXT5KUk5OjhQsX6oknnpAkNWjQQGFhYZo7d666dOkiSdq2bZtWr16tUaNGOVY7ABQVY4w2du8h36pVTpcCBDgUEi5dP0KS9FPTZoo8muNwRUBe3vr1VWX6NMI7AFdzPLgfOHBAP//8s//+hg0b9M033yg+Pl6VK1fWgAEDNGLECKWmpio1NVUjRoxQVFSUunfvLkmKi4tT7969NWjQICUkJCg+Pl6DBw9WnTp1/FeZB4CSzPh8hHa4UuTRHH309mCnywDOyLdypYzPJysqyulSAOC0HA/uX3/9tVq0aOG/P3DgQElSr169NGXKFA0ZMkQ+n099+/bVnj171KhRI82ZM0cxMTH+xzzzzDMKDQ1Vly5d5PP51LJlS02ZMkUhIRySByC4pH6xWB6v1+kyAMD1cn0+/dS0mdNlAMA5cTy4p6WlyRhz2vmWZSkjI0MZGRmnbRMZGamxY8dq7NixhVAhABQfHq9XHkaNAAAASpRicVV5AAAAAACCFcEdAAAAAAAXI7gDAAAAAOBiBHcAAAAAAFyM4A4AAAAAgIsR3AEAAAAAcDGCOwAAAAAALkZwBwAAAADAxQjuAAAAAAC4GMEdAAAAAAAXI7gDAAAAAOBiBHcAAAAAAFyM4A4AAAAAgIsR3AEAAAAAcDGCOwAAAAAALkZwBwAAAADAxQjuAAAAAAC4GMEdAAAAAAAXI7gDAAAAAOBiBHcAAAAAAFyM4A4AAAAAgIsR3AEAAAAAcDGCOwAAAAAALkZwBwAAAADAxQjuAAAAAAC4GMEdAAAAAAAXI7gDAAAAAOBiBHcAAAAAAFyM4A4AAAAAgIsR3AEAAAAAcDGCOwAAAAAALkZwBwAAAADAxQjuAAAAAAC4GMEdAAAAAAAXI7gDAAAAAOBiBHcAAAAAAFyM4A4AAAAAgIsR3AEAAAAAcDGCOwAAAAAALkZwBwAAAADAxQjuAAAAAAC4GMEdAAAAAAAXI7gDAAAAAOBiBHcAAAAAAFyM4A4AAAAAgIsR3AEAAAAAcDGCOwAAAAAALkZwBwAAAADAxQjuAAAAAAC4GMEdAAAAAAAXI7gDAAAAAOBiBHcAAAAAAFws1OkCcAbGSIcPOl3Fucs5esLfByWFOFbKeQuLkizL6SoAAAAAIA9bgvu6deu0a9cu1a1bV9HR0XYsEsZIL7WRNi9zupJzZyIkvXzs7ycvlKxsR8s5LylXSrd/THgHAAAA4DoFOlT+v//9rypVqqRatWrp6quv1rp16yRJXbp00cSJE20pMGgdPli8QrukKCtbv0V212+R3RVVnEK7JG3+sngd3QAAAAAgaOR7xP3111/Xrbfeqvbt26tdu3bq16+ff179+vX12muv6Y477rClyKA3+GcpPMrpKkqmnIPSUxc6XQUAAAAAnFa+g/vIkSN12223afLkyTp69GhAcK9Zs6bGjh1rS4HQsdAezikIAAAAABCM8n2o/Nq1a9WtW7dTzouPj9fu3bvzXRQAAAAAADgm38E9KipKWVlZp5y3ZcsWlSlTJt9FAQAAAACAY/Id3Js2bapx48bJGJNn3pQpU5SWllaQugAAAAAAgApwjvsjjzyiZs2a6YorrlD37t1lWZbeeustDRs2TIsWLdJXX31lZ50AAAAAAASlfI+4N2zYUB999JEOHDigQYMGyRijESNGaP369frwww9Vu3ZtWwo8cuSIHnroIVWrVk1er1fVq1fXo48+qtzcXH8bY4wyMjKUnJwsr9ertLQ0rVmzxpb+AQAAAABwUr5H3CWpRYsWWrt2rX755Rdt375dZcuW1UUXXWRXbZKkJ554Qi+++KKmTp2qSy65RF9//bVuu+02xcXF6R//+IckadSoURo9erSmTJmiiy66SI899phatWqldevWKSYmxtZ6AAAAAAAoSgUK7sddcMEFuuCCC+xYVB5Lly7VDTfcoOuuu06SVLVqVc2cOVNff/21pGOj7WPGjNGDDz6ojh07SpKmTp2qxMREzZgxQ3369CmUugAAAAAAKAr5PlT+iSeeUP/+/U85r3///nrqqafyXdSJmjVrpk8//VTr16+XJH377bdavHixrr32WknShg0blJmZqdatW/sfExERoebNm2vJkiWnXW52drb27dsXcAMAAAAAwG3yHdynTp162vPYL7vsMk2dOjXfRZ3ogQce0M0336waNWooLCxM9erV04ABA3TzzTdLkjIzMyVJiYmJAY9LTEz0zzuVkSNHKi4uzn9LSUmxpV4AAAAAAOyU7+C+cePG057PfuGFF+q3337L76IDzJo1S9OmTdOMGTO0cuVKTZ06VU899VSeHQOWZQXcN8bkmXai9PR0ZWVl+W+bN2+2pV4AAAAAAOyU73Pcw8LCtGPHjlPO2759+xlD8/m4//77NXToUHXr1k2SVKdOHW3cuFEjR45Ur169lJSUJOnYyHuFChX8j9uxY0eeUfgTRUREKCIiwpYaAQAAAAAoLAX6ObiJEyeect7EiRPVsGHDfBd1ooMHD8rjCSwzJCTE/3Nw1apVU1JSkubOneufn5OTo4ULF6pJkya21AAAAAAAgFPyPeI+ePBgXXfddUpLS1Pfvn1VsWJF/f7773rxxRe1aNEiffjhh7YUeP311+vxxx9X5cqVdckll2jVqlUaPXq0br/9dknHDpEfMGCARowYodTUVKWmpmrEiBGKiopS9+7dbakBAAAAAACn5Du4t23bVhMmTNCgQYPUrVs3WZYlY4zi4uI0ceJEtWnTxpYCx44dq4cfflh9+/bVjh07lJycrD59+uiRRx7xtxkyZIh8Pp/69u2rPXv2qFGjRpozZw6/4Q4AAAAAKPYK9DvuvXv3Vrdu3bRkyRLt3LlT5cqVU5MmTRQdHW1XfYqJidGYMWM0ZsyY07axLEsZGRnKyMiwrV8AAAAAANygQMFdkqKjo9WqVSs7agEAAAAAACcpUHA3xmj58uXauHGjfD5fnvl/+9vfCrJ4AAAAAACCXr6D+/r169WhQwf99NNPMsbkmW9ZFsEdAAAAAIACyndw79evnw4dOqRZs2bp0ksv5TfRAQAAAAAoBPkO7l999ZUmTpyozp0721kPAAAAAAA4gSe/DyxVqpRiY2PtrAUAAAAAAJwk38H9tttu04wZM+ysBQAAAAAAnCTfh8rXrl1bM2fOVIcOHXT99dcrISEhT5uOHTsWqDgAAAAAAIJdvoN79+7dJUkbNmzQ+++/n2e+ZVk6evRo/isDAAAAAAD5D+7z58+3sw4AAAAAAHAK+Q7uzZs3t7MOAAAAAABwCvm+OB0AAAAAACh8+R5xl6SffvpJ//nPf7R27Vr5fL6AeZZl6dNPPy1QcQAAAAAABLt8B/fVq1fryiuvVMWKFfXzzz/r0ksv1a5du7RlyxalpKToggsusLNOAAAAAACCUr4Plf/nP/+pNm3aaM2aNTLGaPLkydq8ebPee+89HTp0SI899piddQIAAAAAEJTyHdxXrlypXr16yeM5tojc3FxJ0nXXXafBgwcrPT3dngoBAAAAAAhi+Q7ue/bsUXx8vDwej8LCwrRnzx7/vIYNG2rlypW2FAgAAAAAQDDLd3CvWLGidu3aJUm68MILtWjRIv+87777TqVKlSp4dQAAAAAABLl8X5yuWbNmWrJkiW688Ub16NFDw4YN07Zt2xQeHq4pU6bolltusbNOAAAAAACCUr6D+4MPPqitW7dKkh544AFlZmZq+vTpsixLXbp00ZNPPmlbkQAAAAAABKt8B/cLLrjA/5NvISEheu655/Tcc8/ZVhgAAAAAACjAOe633367NmzYcMp5Gzdu1O23357vogAAAAAAwDH5Du5TpkzRzp07Tzlv165dmjp1ar6LAgAAAAAAx+Q7uJ/JH3/8oYiIiMJYNAAAAAAAQeW8znFftGiRFixY4L8/adIkffzxxwFtfD6f3nnnHdWqVcuWAgEAAAAACGbnFdznz5+v4cOHS5Isy9KkSZNO2a5KlSp6/vnnC14dAAAAAABB7ryC+5AhQ3TPPffIGKPy5cvrk08+Uf369QPaREREqFSpUrYWCQAAAABAsDqv4O71euX1eiVJGzZsUHJyssLCwgqlMAAAAAAAUICL0yUmJsrn8wVMe+211zR06FDNmzevwIUBAAAAAIACBPeePXvq3nvv9d9/7rnn1K1bN40aNUpt2rTRhx9+aEuBAAAAAAAEs3wH96+++kpt27b133/uued0yy23aO/everYsaOeeuopWwoEAAAAACCY5Tu479y5UxUrVpR07Hz3X3/9Vf3791dsbKx69+6t1atX21YkAAAAAADBKt/BPSoqSllZWZKkzz//XKVKlVLDhg0lSZGRkTpw4IA9FQIAAAAAEMTO66ryJ6pTp46ef/55ValSRS+88IJatGghy7IkSZs2bVJSUpJtRQIAAJQkxhiZky7yi6KVe8L/P5fnwnGW1+vPEgDyyndwf/jhh9W+fXvVrVtX4eHhAVeS/+CDD/L8vjsAAACOhfaN3XvIt2qV06Xgf35q2szpEoKet359VZk+jfAOnEa+g/tf/vIXrV27VitWrFDdunVVvXr1gHl169a1oz4AAIASxfh8hHbgJL6VK2V8PllRUU6XArhSvoO7JFWpUkVVqlTJM71Pnz4FWSwAAEBQSP1isTxer9NlAI7J9fk44gE4BwUK7gAAAMg/j9crDyOMAICzOK+ryoeEhOirr7469kCPRyEhIae9hYayTwAAAAAAgII6r3T9yCOPqFKlSv6/uXgEAAAAAACF67yC+7Bhw/x/Z2Rk2F0LAAAAAAA4Sb6OZ9+5c6f+85//aNGiRdq6daskKTk5WS1atNCdd96phIQEW4sEAAAAACBYnXdw//TTT9WpUyft27dPISEhKlu2rIwxWrdunebNm6ennnpKs2fP1tVXX10Y9QIoZMYY+Y74zvtxBw8fPeFvn2SFnPcyvKFeTsEBAAQlY4yM7/w/f4u73BPWOTcI11+SLC/ff3B25xXcd+7cqa5duyouLk6TJk3Stddeq6j/XQn14MGDev/99zV48GB17txZa9euZeQdKGaMMfrbR3/TNzu/Of/H5oZJ+pckKe215rI8h897GfXK19PUtlP58AIABBVjjDZ27yHfqlVOl+KoYP1ZOG/9+qoyfRrff3BG53VV+cmTJ+vo0aP64osv1LlzZ39ol6SoqCh16dJFixcv1uHDhzV58mTbiwVQuHxHfPkK7ZJkeQ4rpuZQxdQcmq/QLkmrdqzK12g/AADFmfH5gj60BzPfypVBebQFzs95jbjPmTNHt99+u//K8qdSuXJl3Xbbbfr44481ZMiQAhcIwBkLuiyQN9RbJH35jviU9lpakfQFAICbpX6xWB5v0Xz+wlm5Pl/QHmWA83dewX3t2rXq37//WdtdddVVmjlzZr6LAuA8b6hXUWFRZ28IAABs4/F65Yni8xdAoPM6VH7v3r0qX778WduVL19ee/fuzW9NAAAAAADgf84ruGdnZyssLOys7UJDQ5WTk5PvogAAAAAAwDHn/XNw69atU2jomR/2448/5rsgAAAAAADw/847uN96661nbWOM4ecMAAAAAACwwXkF95dffrmw6gAAAAAAAKdwXsG9V69ehVUHAAAAAAA4hfO6OB0AAAAAAChaBHcAAAAAAFyM4A4AAAAAgIud91XlAaAkMcbI+HxOl1EguSfUn1vM18XyevlVEgAAgJMQ3AEELWOMNnbvId+qVU6XYpufmjZzuoQC8davryrTpxHeAQAATsCh8gCClvH5SlRoLwl8K1cW+yMgAAAA7MaIOwBISv1isTxer9NlBK1cn6/YHy0AAABQWAjuACDJ4/XKExXldBkAAABAHhwqDwAAAACAixWL4L5lyxbdcsstSkhIUFRUlOrWrasVK1b45xtjlJGRoeTkZHm9XqWlpWnNmjUOVgwAAAAAgD1cH9z37Nmjpk2bKiwsTB999JF++OEHPf300ypdurS/zahRozR69GiNGzdOy5cvV1JSklq1aqX9+/c7VzgAAAAAADZw/TnuTzzxhFJSUvTyyy/7p1WtWtX/tzFGY8aM0YMPPqiOHTtKkqZOnarExETNmDFDffr0KeqSAQAAAACwjetH3N999101bNhQN910k8qXL6969epp4sSJ/vkbNmxQZmamWrdu7Z8WERGh5s2ba8mSJaddbnZ2tvbt2xdwAwAAAADAbVwf3H/99VeNHz9eqamp+uSTT3TXXXfp3nvv1X//+19JUmZmpiQpMTEx4HGJiYn+eacycuRIxcXF+W8pKSmFtxIAAAAAAOST64N7bm6u6tevrxEjRqhevXrq06eP7rjjDo0fPz6gnWVZAfeNMXmmnSg9PV1ZWVn+2+bNmwulfgAAAAAACsL1wb1ChQqqVatWwLSaNWtq06ZNkqSkpCRJyjO6vmPHjjyj8CeKiIhQbGxswA0AAAAAALdxfXBv2rSp1q1bFzBt/fr1qlKliiSpWrVqSkpK0ty5c/3zc3JytHDhQjVp0qRIawUAAAAAwG6uv6r8fffdpyZNmmjEiBHq0qWLvvrqK02YMEETJkyQdOwQ+QEDBmjEiBFKTU1VamqqRowYoaioKHXv3t3h6gEAAAAAKBjXB/fLL79cs2fPVnp6uh599FFVq1ZNY8aMUY8ePfxthgwZIp/Pp759+2rPnj1q1KiR5syZo5iYGAcrBwAAAACg4Fwf3CWpffv2at++/WnnW5aljIwMZWRkFF1RAAAAAAAUAdef4w4AAAAAQDAjuAMAAAAA4GLF4lB5AMC5McbI+HxOl3Heck+oObcY1i9Jltcry7KcLgMAAJRABHcAKCGMMdrYvYd8q1Y5XUqB/NS0mdMl5Iu3fn1VmT6N8A4Ahay47qQ+WUnYaX0idmAXLoI7AJQQxucr9qG9OPOtXCnj88mKinK6FAAosUrKTuqTFded1idiB3bhIrgDQAmU+sViebxep8sICrk+X4n4wgUAxQE7qd2LHdiFi+AOACWQx+uVhw9OAEAJxk5qd2AHdtEguAMAAAAodthJjWDCz8EBAAAAAOBiBHcAAAAAAFyM4A4AAAAAgIsR3AEAAAAAcDGCOwAAAAAALkZwBwAAAADAxQjuAAAAAAC4GMEdAAAAAAAXI7gDAAAAAOBiBHcAAAAAAFws1OkCAAAoLMYYGZ+vUPvIPWH5uYXclyRZXq8syyr0fgAAgHsQ3AEAJZIxRhu795Bv1aoi6/Onps0KvQ9v/fqqMn0a4R0AgCDCofIAgBLJ+HxFGtqLim/lykI/igAAALgLI+4AgBIv9YvF8ni9TpdRILk+X5GM6AMAAPchuAMASjyP1ytPVJTTZQAAAOQLh8oDAAAAAOBiBHcAAAAAAFyM4A4AAAAAgIsR3AEAAAAAcDGCOwAAAAAALkZwBwAAAADAxfg5OBQfxkiHD9q7zJyDp/7bTmFRkmUVzrIBAAAAlHgEdxQPxkgvtZE2Lyu8Pp66sHCWm3KldPvHhHcAAAAA+cKh8igeDh8s3NBemDZ/af+RAgAAAACCBiPuKH4G/yyFRzldxdnlHCy8UXwAAAAAQYPgjuInPEoKj3a6CgAAAAAoEhwqDwAAAACAixHcAQAAAABwMYI7AAAAAAAuRnAHAAAAAMDFCO4AAAAAALgYwR0AAAAAABcjuAMAAAAA4GIEdwAAAAAAXIzgDgAAAACAixHcAQAAAABwMYI7AAAAAAAuRnAHAAAAAMDFCO4AAAAAALhYqNMFAAAA5IcxRsbnc7qM85Z7Qs25xbB+SbK8XlmW5XQZABA0CO4AAKDYMcZoY/ce8q1a5XQpBfJT02ZOl5Av3vr1VWX6NMI7ABQRDpUHAADFjvH5in1oL858K1cWy6MdAKC4YsQdAAAUa6lfLJbH63W6jKCQ6/MV26MEAKA4I7gDAIBizeP1yhMV5XQZAAAUGoI7AACACxSHi+0VpwvrcQE9oGi2K0W9XQjW9zbBHQAAwGHF8WJ7bj9kngvoIdg5sV0piu1CsL63uTgdAACAw7jYnv24gB6CXUndrgTre5sRdwAAABfhYnsFwwX0gLxKwnYl2N/bBHcAAAAX4WJ7AOzGdqX441B5AAAAAABcjOAOAAAAAICLEdwBAAAAAHAxgjsAAAAAAC5W7IL7yJEjZVmWBgwY4J9mjFFGRoaSk5Pl9XqVlpamNWvWOFckAAAAAAA2KVbBffny5ZowYYIuvfTSgOmjRo3S6NGjNW7cOC1fvlxJSUlq1aqV9u/f71ClAAAAAADYo9gE9wMHDqhHjx6aOHGiypQp459ujNGYMWP04IMPqmPHjqpdu7amTp2qgwcPasaMGQ5WDAAAAABAwRWb4N6vXz9dd911uuaaawKmb9iwQZmZmWrdurV/WkREhJo3b64lS5acdnnZ2dnat29fwA0AAAAAALcJdbqAc/Hqq69q5cqVWr58eZ55mZmZkqTExMSA6YmJidq4ceNplzly5EgNHz7c3kIBAAAAALCZ60fcN2/erH/84x+aNm2aIiMjT9vOsqyA+8aYPNNOlJ6erqysLP9t8+bNttUMAAAAAIBdXD/ivmLFCu3YsUMNGjTwTzt69KgWLVqkcePGad26dZKOjbxXqFDB32bHjh15RuFPFBERoYiIiMIrHAAAAAAAG7g+uLds2VLff/99wLTbbrtNNWrU0AMPPKDq1asrKSlJc+fOVb169SRJOTk5WrhwoZ544gknSgaAszLGyPh8ti4z94Tl5dq87OMsr/eMRzMBAADAfq4P7jExMapdu3bAtOjoaCUkJPinDxgwQCNGjFBqaqpSU1M1YsQIRUVFqXv37k6UDABnZIzRxu495Fu1qtD6+Klps0JZrrd+fVWZPo3wDgAAUIRcH9zPxZAhQ+Tz+dS3b1/t2bNHjRo10pw5cxQTE+N0aQCQh/H5CjW0FybfypUyPp+sqCinSwEAAAgaxTK4L1iwIOC+ZVnKyMhQRkaGI/XIGOnwQXuXmXPw1H/bKSxKYtQMcFTqF4vl8XqdLuOscn2+QhvFBwAAwJkVy+DuKsZIL7WRNi8rvD6eurBwlptypXT7x4T3s8nvjhk7dr6wc6XE83i98jB6DQAAgDMguBfU4YOFG9oL0+Yvj9UfHu10Je5l146Z/O58YecKAAAAEPQI7nYa/LMUXgxGznIOFt4ofknj9I4Zdq4AAAAAQY/gbqfwKAJWSVaUO2bYuQIAAADgfwjuwLlixwwAAAAAB3icLgAAAAAAAJwewR0AAAAAABcjuAMAAAAA4GIEdwAAAAAAXIzgDgAAAACAixHcAQAAAABwMYI7AAAAAAAuRnAHAAAAAMDFCO4AAAAAALgYwR0AAAAAABcLdboAAAAAAMWPMUbG5yvSPnNP6C+3iPuWJMvrlWVZRd4vQHAHAAAAcF6MMdrYvYd8q1Y5VsNPTZsVeZ/e+vVVZfo0wjuKHMEdAACbFcYoVFGMMjGSBOBcGZ/P0dDuFN/KlTI+n6yoKKdLQZAhuANAEMlvoLQjNAZLKCyKUajCGmViJAlAfqR+sVger9fpMgpVrs/nyAg/cBzBHQCChF2BMr9fXIIlFBbnUShGkgDkh8frlYftBlCoCO4AECScDpTBGAqLyygUI0kAALgbwR0AglBRBspgDoWMQgEAADsQ3AEgCBEoAQAAig+P0wUAAAAAAIDTY8QdAAAAAGAbfhbVfgR3AAAAAIAt+FnUwsGh8gAAAAAAWzj9KzYFcfwXcNyIEXcAAAAAgO34WVT7ENwBAADgOvk9R9aO82DdfJ4rUJzwKzb2IbgDAADAVew6Rza/I2huPs8VQHDiHHcAAAC4itPnyLr5PFcAwYkR9+LMGOnwwfN/XM7BU/99PsKiJPZCAwCAQlaU58gWh/NcAQQngntxZYz0Uhtp87KCLeepC/P3uJQrpds/JrwDAIBCxTmyAEBwL74OHyx4aC+IzV8eqyE82rkaAAAAgPPARQ9RXBHcS4LBP0vhRbQnOudg/kfpAQAAAIdw0UMUZwT3kiA8ipFvAAAjSQBwBm656KHFqR/IB4I7AAAlACNJAHDuuOghihuCOwAAJQAjSQBw7rjoIYobgjsAACUMI0kAAJQsBHcAAEoYRpIAnCuujQEUDwR3AACA85DfoHMmdoSgsyEk4WRcGwMoPgjuAAAA58iuoHMmhXXqASEJJ+PaGEDxQXAHAAA4R04HnYIgJOFMuDYG4G4EdwAAgHwoyqBTEIQknAuujQG4G8EdAAAgHwg6AICi4nG6AAAAAAAAcHoEdwAAAAAAXIzgDgAAAACAixHcAQAAAABwMYI7AAAAAAAuxlXlAQAAigFjjIzPd96Pyz3hMbn5eLwkWV6vLMvK12MBAAVHcAcAAHA5Y4w2du8h36pVBVpOfn/P3Vu/vqpMn0Z4BwCHcKg8AACAyxmfr8ChvSB8K1fma7QfAGAPRtwBAACKkdQvFsvj9RZJX7k+X75H6QEA9iG4AwAAFCMer1eeqCinywAAFCEOlQcAAAAAwMUI7gAAAAAAuBjBHQAAAAAAFyO4AwAAAADgYgR3AAAAAABcjOAOAAAAAICLEdwBAAAAAHAx1wf3kSNH6vLLL1dMTIzKly+vG2+8UevWrQtoY4xRRkaGkpOT5fV6lZaWpjVr1jhUMQAAAAAA9nF9cF+4cKH69eunL7/8UnPnztWRI0fUunVr/fnnn/42o0aN0ujRozVu3DgtX75cSUlJatWqlfbv3+9g5QAAAAAAFFyo0wWczccffxxw/+WXX1b58uW1YsUKXX311TLGaMyYMXrwwQfVsWNHSdLUqVOVmJioGTNmqE+fPk6UDQAAAACALVw/4n6yrKwsSVJ8fLwkacOGDcrMzFTr1q39bSIiItS8eXMtWbLktMvJzs7Wvn37Am4AAAAAALhNsQruxhgNHDhQzZo1U+3atSVJmZmZkqTExMSAtomJif55pzJy5EjFxcX5bykpKYVXOAAAAAAA+VSsgvs999yj7777TjNnzswzz7KsgPvGmDzTTpSenq6srCz/bfPmzbbXCwAAAABAQbn+HPfj+vfvr3fffVeLFi1SpUqV/NOTkpIkHRt5r1Chgn/6jh078ozCnygiIkIRERGFVzAAAAAAADZw/Yi7MUb33HOP3nrrLX322WeqVq1awPxq1aopKSlJc+fO9U/LycnRwoUL1aRJk6IuFwAAAAAAW7l+xL1fv36aMWOG3nnnHcXExPjPW4+Li5PX65VlWRowYIBGjBih1NRUpaamasSIEYqKilL37t0drh4AAAAAgIJxfXAfP368JCktLS1g+ssvv6xbb71VkjRkyBD5fD717dtXe/bsUaNGjTRnzhzFxMQUcbUAAAAAANjL9cHdGHPWNpZlKSMjQxkZGYVfEAAAAADAdYwxMj7feT8u94TH5Obj8ZJk/e9o8MLi+uAOAAAAAMCZGGO0sXsP+VatKtByfmraLF+P89avryrTpxVaeHf9xekAAAAAADgT4/MVOLQXhG/lynyN9p8rRtwBAAAAACVG6heL5fF6i6SvXJ8v36P054PgDgAAAAAoMTxerzxRUU6XYSsOlQcAAAAAwMUI7gAAAAAAuBjBHQAAAAAAFyO4AwAAAADgYgR3AAAAAABcjOAOAAAAAICLEdwBAAAAAHAxgjsAAAAAAC5GcAcAAAAAwMUI7gAAAAAAuBjBHQAAAAAAFyO4AwAAAADgYgR3AAAAAABcjOAOAAAAAICLEdwBAAAAAHAxgjsAAAAAAC5GcAcAAAAAwMUI7gAAAAAAuBjBHQAAAAAAFyO4AwAAAADgYgR3AAAAAABcjOAOAAAAAICLEdwBAAAAAHAxgjsAAAAAAC5GcAcAAAAAwMUI7gAAAAAAuBjBHQAAAAAAFyO4AwAAAADgYgR3AAAAAABcjOAOAAAAAICLEdwBAAAAAHAxgjsAAAAAAC5GcAcAAAAAwMUI7gAAAAAAuBjBHQAAAAAAFyO4AwAAAADgYgR3AAAAAABcjOAOAAAAAICLEdwBAAAAAHAxgjsAAAAAAC5GcAcAAAAAwMUI7gAAAAAAuBjBHQAAAAAAFyO4AwAAAADgYgR3AAAAAABcjOAOAAAAAICLEdwBAAAAAHAxgjsAAAAAAC5GcAcAAAAAwMUI7gAAAAAAuBjBHQAAAAAAFyO4AwAAAADgYgR3AAAAAABcjOAOAAAAAICLEdwBAAAAAHAxgjsAAAAAAC5GcAcAAAAAwMVKVHB/4YUXVK1aNUVGRqpBgwb6/PPPnS4JAAAAAIACKTHBfdasWRowYIAefPBBrVq1SldddZXatWunTZs2OV0aAAAAAAD5VmKC++jRo9W7d2/9/e9/V82aNTVmzBilpKRo/PjxTpcGAAAAAEC+hTpdgB1ycnK0YsUKDR06NGB669attWTJklM+Jjs7W9nZ2f77WVlZkqR9+/adZ+d/StlG/3uwFH70/B6fX/RLv4Xg4OGDOuo7+r9u9+lI2JES3W/uwYM6cPT/+/UcoV/6pV/6pV/6pV/6pV/6Lbp+j+dPY8wZ21nmbC2Kga1bt6pixYr64osv1KRJE//0ESNGaOrUqVq3bl2ex2RkZGj48OFFWSYAAAAAAHls3rxZlSpVOu38EjHifpxlWQH3jTF5ph2Xnp6ugQMH+u/n5ubqjz/+UEJCwmkfAwAAAACAXYwx2r9/v5KTk8/YrkQE97JlyyokJESZmZkB03fs2KHExMRTPiYiIkIREREB00qXLl1YJQIAAAAAkEdcXNxZ25SIi9OFh4erQYMGmjt3bsD0uXPnBhw6DwAAAABAcVMiRtwlaeDAgerZs6caNmyoxo0ba8KECdq0aZPuuusup0sDAAAAACDfSkxw79q1q3bv3q1HH31U27ZtU+3atfXhhx+qSpUqTpcGAAAAAEC+lYirygMAAAAAUFKViHPcAQAAAAAoqQjuAAAAAAC4GMEdAAAAAAAXI7gDAAAAAOBiBHcAAAAAAFyM4A4AAAAAQD4dPHhQHTt21A8//FBoffBzcMB5OnTokHbs2KHKlSs7XQpwVtu3b9f333+v3bt3y+PxqGLFimrQoIEiIiKcLq1IfPbZZ1q5cqU8Ho+uvPJKNWnSxOmSUADGGC1btkyrV6/W7t27ZVmW4uPjVbt2bTVq1EiWZTldIvLhq6++0uWXX+7Y85eVlaXw8HB5vV7/tM8//1yrV69WpUqV1K5dO4WGhjpSW1HJzs6WZVkKDw93uhSUILt379avv/6q6tWrKyEhwelyClVWVpbKlCmjhQsX6qqrriqcTgxs8+effzpdQpEqqvXNzs42mZmZZvv27SY7O7tI+jyTN954w3g8HqfLsJ3b/s/B4PDhw+add94xu3fvtn3ZX375pWnWrJnxeDx5bjExMWbQoEFFvs3auXNnoayrMcY89thj5l//+pf//r59+8xVV11lPB6PsSzLWJZlPB6PueGGG0xOTk6h1OAWhfm6ctLMmTNNpUqVAp7TE5/bSpUqmZkzZzpdZolVmO9fy7JM5cqVzbBhw8xvv/1WKH2cysGDB82NN95oPB6PCQ0NNQMGDDDGGNO7d++A11adOnXMH3/8UWR1FZZly5bl2e7PmTPH1K1b1//50KBBA/PZZ585VCHsVJTf6x5//HFTpUoVc8EFF5iJEycaY4x5+umnTXh4uP/9NXToUNv7nT9/vmnVqpWpUaOG6dy5s1m1alWeNl9++aVt39tjYmLOeLMsy0RFRZmYmBgTGxtrS58nIrifp9mzZ5sbb7zRdOvWzSxdutQYc2yjV716dePxeExCQoIZMWKErX3WrFnTPPjgg+aXX36xdbnnwon1NcaYXbt2maFDh5qLL77YhISE+D9QQkJCzMUXX2z++c9/ml27dtne77kojOB+6NAhM336dPPvf//bvPfee6ds88svv5jbbrvN1n7d/H9euHChadGiha3LPHLkiJk0aZJp166dqVixoomMjDRRUVEmNTXV/P3vfzfffPONrf2dzd69e43H4zGLFi2ydblffPGFiYyMNOXKlTOdOnUy3bt3NxdccIEJDQ016enp5s477zSlS5c2jRo1MgcPHrS176VLl5p58+YFTJs0aZKpXLmy//WVmppq3nrrLVv7TU1NNZMnT/bfv+uuu0xUVJR59tlnzfr16826devM6NGjjdfrNQ899JCtfTv1/j2dwnpdnc6ff/5pxowZY7p162a6d+9uxo4da/vr6tVXXzWWZZnWrVubmTNnmjVr1pht27aZbdu2mTVr1piZM2eatm3bGo/HY1577TVb+3Zqu7F+/XozbNgw89hjj5nff//dGGPMDz/8YG688UZTq1Ytc8MNN5iVK1fa2qdT71/Lsky5cuWMZVkmJCTEtGnTxrz22muFvpPt0UcfNaGhoaZPnz5myJAhpnTp0mbAgAEmLi7OvPzyy2b16tVmwoQJJiYmxgwZMqRQazkVuz8HPR6PWbZsmf/+4sWLTWhoqElJSTF33323ueuuu0ylSpVMeHi4WbFihW39GuPM6/lUimJ75eT3die+102bNs1YlmWuvPJK0759exMeHm5efPFFExISYu6++27z8ssvm549exqPx2NmzJhhW78rVqwwYWFhJikpyfzlL38xpUuXNuHh4eaFF14IaGdncLcsy1SqVMnceuuteW7du3c3lmWZtm3b+qfZjeB+Hj744ANjWZYpXbq0SUpKMlFRUebTTz810dHRJi0tzQwePNg0bdrUeDweW/f6H9/rGxISYlq2bGlmzpxZJCOiTq3vr7/+aipWrGjCw8NNmzZtzP33329GjhxpRo4cae6//37Ttm1bExERYVJSUsyvv/5qW7/Dhw8/p1vXrl1tDe579+41tWvXDti736hRozyjDnZueIxx7v98ruzeQbJnzx5z+eWXB/yfLcsyoaGhpmHDhiY2NtaEhoaaUaNG2danMcb079//tLc+ffoYy7JMx44dTf/+/c29995rS58tWrQwTZo0Mfv37/dPO3LkiOndu7dp2rSpMcaYjRs3mnLlyplhw4bZ0udxTZs2DRj5njp1qrEsyzRo0MD885//NOnp6aZu3bomJCTEfPjhh7b16/V6zYIFC/z3y5cvb5566qk87UaOHGmqVq1qW79OvX+deF0Zc+wL6Xfffee//8cff5hatWoZy7JMqVKlTHR0tLEsy1x66aVm3759tvVbt25dc+edd5613R133GHq1q1rW79ObTfWrFljSpUq5e+3QoUK5scffzRJSUkmMTHRNG/e3MTHx5tSpUqZtWvX2tavU+9fy7LMsmXLzIIFC0zPnj1NVFSU8Xg8ply5cmbgwIFmzZo1tvV1oho1apjhw4f777/33nvG4/Hk2XY8/vjjpkaNGoVSw5nY/Tl4/P98XJs2bcwll1xi9u7d65/2xx9/+Ecv7eLU69mp7ZVT39ud+l7XuHFj0717d//9559/3kRGRpq+ffsGtOvatautO6JuuOEG07hxY3PgwAFjjDFZWVmmR48exuPxmH//+9/+dnZ+/k6cONGULl3a3HDDDWbLli0B8/bs2WMsyzILFy60pa9TIbifh+bNm5vWrVubQ4cOGWOMue+++0zZsmUDXqzGGNOpUyfTpEkT2/q1LMs8++yz5tZbbzWlSpUyHo/HxMfHm3vvvdd8++23tvVzMqfWt2PHjqZ27dpm48aNp22zceNGU6dOHdOpUyfb+j3xS9nZbnZ+kA4dOtTEx8ebN954w2zatMm8/PLLJikpyVSoUMGsXr3a387uL/5O/Z83btx4Trfx48fbur533323iYuLMzNnzjT79+83hw4dMp9++qm56KKLTHp6usnJyTHPPvusCQkJOe2oaX6c7XV18qHcdihVqpR5++2380zftGmT8Xg8ZvPmzcYYY5566ilz0UUX2dLncaVLlw74Qn/JJZeYLl265Gn317/+1dbtRpkyZcz777/vvx8aGnrKD8958+aZ8PBw2/p16v3rxOvqeL8nfvG/5557THR0tHnzzTdNbm6uyc3NNbNmzTJer9ekp6fb1m9kZGTAjpnTmT9/vomMjLStX6e2G926dTMXX3yxWbt2rdm5c6e54YYbTI0aNczll1/u3yG3c+dOU7NmTdOrVy/b+nXq/Xvy6yorK8uMHz/eNGzY0P8abty4sZk8ebL/S7odoqKizPz58/339+/fbyzLMp9//nlAu88++8xER0fb1q9Tn4Mn/59jY2PNlClT8rSbMGGCSUxMtK1fp17PTm2vnPre7tT3uvj4ePPBBx/47+/cudNYlmU++uijgHZvvvmmiY+Pt63fihUrmjfffDPP9Iceesh4PB7/Tki7P3+3bNli2rdvb+Li4syLL77on753716Cu5skJCQEfDHcunWrsSzLfPzxxwHt3nvvPVO2bFnb+j1xw7N//34zYcIE06hRI/+H2RVXXGH+85//2Lq30Bjn1rd06dJm9uzZZ2331ltvmdKlS9vWb7ly5cydd95pdu3adcbbSy+9ZOsGoGbNmuaZZ54JmPb777+bSy+91JQrV86/kbd7w+PU//n46/Zcb3apUKGCee655/JMP35Y+fEvD3//+99N8+bNbev3qquuMsnJyacM0oW1dzYmJsa8++67eaZv2bLFWJblH72aN2+erSHHmGMj3yeuT2hoqJk7d26edh9++KGJioqyrd9rr73W9OjRw3+/bt265oknnsjT7rHHHjNVqlSxrV+n3r9OvK6MyftFuFKlSubhhx/O0+74oZp2SU5O9p83eSYTJ040ycnJtvXr1HajcuXKAYFq7dq1xrKsPIeZjh8/3lSuXNm2fp16/578ujrRd999Z+69916TkJBgLMuy9bzR5ORk88Ybb/jv//rrr8ayrDyfja+99pqtryunPgdP/j+HhoaaxYsX52m3YMECExYWZlu/Tr2endpeOfW93anvdScf8XbkyBFjWVae0y0WLlxo647z6Ojo054SlpGRYSzLMsOGDbP98/e4V155xZQtW9Y0a9bM/Pjjj0US3Ev2JTJtdvjwYUVGRvrvH7/66MlX4IyIiNDhw4cLpYZSpUrpjjvu0B133KEffvhBkyZN0rRp03TXXXdp4MCB6tKli1566SVb+nJqfY8cORJwZdfT8Xq9OnLkiG391qtXT+vXrz/rVS9jY2Nt61OSNm7cqPr16wdMq1ixohYsWKBrrrlGf/nLXzR37lxb+5Sc+z+Hh4eradOmuvbaa8/Y7rvvvtO0adNs6/ePP/5QrVq18ky/5JJLlJ2drd9++021a9fWtddeq1tvvdW2fhctWqRx48apZ8+eatOmjcaNG6fExERJKrQrKDdu3FijR49W69atA64e//jjjysmJkapqamSjl1FuFSpUrb2feGFF2rp0qW6+uqrJUkJCQnat29fnnb79u1TWFiYbf2mp6crLS1NZcqU0ZAhQ/TEE0+oS5cusixLrVu3liR99NFHevzxxzVgwADb+nXq/evE6+pUtm3bpubNm+eZnpaWpjFjxtjWz0033aQHHnhAsbGx6ty5szyewF+zzc3N1ZtvvqmhQ4eqZ8+etvXr1HZj+/btql69uv/+8b8vuOCCgHYXX3yxMjMzbevXqffvmdSpU0fPPvusnnzySc2ePdu27ziSdPnll+vxxx9XgwYNFBsbq6FDh6pGjRoaPXq02rZtq8jISB08eFDPPPOMLr30Utv6depzUJIWLFig33//XZJUtmxZ7d69O0+bXbt22frZ4NTr+WRFtb06UVF+b3fqe125cuW0ZcsW/32Px6M+ffqofPnyAe22b9+u0qVL29ZvSkqKfvjhh1NewX3YsGGSpOHDh2v58uW29XmiW265Ra1bt1bfvn1Vt25d9e/fv/A/fwttl0AJ1KBBA9OvXz///cmTJ5vIyMg8FywZMGCArefYnWlPtDHG5OTkmNdff920bdvWhIaG2tavU+vbsmVLk5aWdsbD4Q4cOGCaN29uWrVqZVu/999/vylTpsxZ23300Ue2niNbqVKlUx7qY8yxkbP69eub+Ph4M3bsWFv3GDr1f27QoIHp2LHjWdvZfW5fjRo1Tnko3Ouvv248Ho/ZuXOnMebYIZGlSpWyrd/jNmzYYFq1amXKlCljJk2aZIwpvMOqli1bZiIjI02FChXMzTffbG699VZTo0YN4/F4As7FHT58uK2jhMYcu4psmTJlzFdffWWMMeaf//ynqV27tv/wfGOOHapXo0YNc/3119va9/Tp0/2HJZYtW9Z4vd6AUSvLskyHDh38p//Ywan373FF+boy5tjn0fHn1hhjkpKS8lzMzJhj20k7Dy0+cOCAadWqlX/EtVGjRqZdu3bm2muvNY0aNTKxsbHG4/GY1q1b2/prCU5tN8qWLWveeecd//2jR4+ayMjIPFdM/uCDD2wdgXbq/Xu27zmFZeXKlQHbifj4ePPjjz+aCy+80H/Bq8TERBMSEnJOp2qcK6c+B091Ws19992Xp93QoUNN/fr1bevXqdezU9srp763O/W97rrrrjN33333WdsNHjzY1u8ct99+u2nZsuUZ2zz66KO2nzJ2Km+88YZJSkpixN1N+vTpoz59+uj7779XXFycPv74Y40bN07/+Mc/dOjQITVs2FBffvmlJkyYoCeeeKLI6goLC1Pnzp3VuXNn/15UOzi1vqNGjVKLFi1UvXp1de7cWXXq1FF8fLwsy9Lu3bv1/fff66233tLBgwe1YMEC2/p95JFH1K9fv7O2a9u2rTZs2GBbv7Vr19a8efPUsWPHPPNKly6tTz/9VC1btrR1lFBy7v9cr149ffrpp+fU1hhjW7933323Bg4cKJ/Pp3bt2ikiIkJffPGFnnzySbVr105ly5aVJP3yyy+qUqWKbf0eV7VqVc2ZM0eTJk3S4MGDNX369ELbTlxxxRVatGiRHn74Yb3//vvKyclRrVq1NGnSJN12223+dp06ddLNN99sa9/33nuvPvvsMzVp0kTNmzdX7dq1tXnzZqWmpqpGjRqSpLVr1yo6OlpvvvmmrX13795daWlpmjhxoj7//HNt2bJFubm5SkhI0KWXXqpOnTqpVatWtvbp1Pv3uKJ8XR3XokUL/4i3z+fT6tWr1bJly4A2P//8s/8IADtER0drzpw5+vDDDzV79mytWbNGv/zyi6Rjo8LdunXTX//6V7Vt29a2PiXnthsXXnihvv32W3Xo0EHSsREsn8+Xp91PP/2klJQU2/p16v3bq1cvlStXzrblnat69eppxYoVmj59usLDw9WrVy9VqVJF8+bN0wMPPKBvv/1Wl1xyiQYMGHDKkdqC9OvE5+D8+fPzTIuLi8szbdOmTeratatt/Tr1epac2V6dTWF9b3fqe93jjz+uvXv3nrXdvn37Ar6DFNQtt9yi8ePHa9euXf5t8ckefvhhRUVF6f3337et31Pp1KmT2rZtq127dikpKanwOiq0XQIl1MMPP2wqVapkqlWrZsaMGWOMMeall14yUVFR/r2XnTp1svXqkU7tiTbGmfU1xpgff/zR3HTTTQH9HL9FRUWZrl27mh9//NHWPp3y9NNPm9jY2DP+RuzxKxvbvcfQif/z4sWLA67iezo7d+60dYTDGGMGDhxoQkNDA0Zgmzdv7h81M+bYObITJkywtd+T/f7776Z9+/YmLCzMeDyeQt0764SjR4+aZ555xlSrVi3P6yo6Otp0797dkZ/JKQxOvn9PVhSvq1P9BM7xz4YTNW3a1NaLHznJie3Gc889d06/+HDFFVec09X2z0cwvX+d4uTnoBOcej07tb1y8nt7MH1/DkaWMTbuygtie/fu1bp165ScnGz73kI3Kqr1PXr0qH755Rf/OVgJCQm64IILFBISUmh9FrXc3Fz5fD55vd48522e6NChQ9q+fXuhjAYHw//5uK1bt2rp0qU6fPiwatWqZet5i+frjTfe0Jo1a3TbbbepcuXKjtVRmH7//feAke+S9rpyw/v3ZG54XW3atEmxsbG2ns/opK1bt+rLL7/0H73i5HbjRPv27VNkZGSea8/YpaS/f+Euhf16Pp2Str2Sgut7XTAhuBeCQ4cOaceOHUX+hcmpfovS7t279euvv6p69epnvYhcSRMMz28wC7bnN9jWFyUbn/soDDy/RSPY/s/Btr4lyemHB5BvH3zwgapVq1Yi+l2wYIFat26tmjVr6qabbtI333yTp82yZcts34M3YsQIVa1aVRdeeKEmTZokSRo9erSSk5N15ZVXKikpSenp6bb26XYl6fkNtn7PhVPbDaeUpPUNttezm99HTuFzH4WB7WTRCLb3UUl6XQUbgjtOa+XKlWrdurW+//57JScna968eWrUqJHGjx9fqP1Onz5dDz30kCpUqKCaNWuqX79++s9//qMhQ4aod+/emjx5sm6++WaNGjVKM2fOLNRaSjKnnt9g6xclW7C9nnkflWw8vygMwfa6Crb1RdHhqvLn4dFHHz2ndj/88EOJ6bdhw4aaO3euoqOjtW/fPvXt21f33HOP9u3bpwceeMDW/o57/vnndfPNN2v69OmSpBdeeEEDBgxQnz599Pzzz0uSbr31VuXk5GjixIm2XxW7qAXb8xuM/Z4Lu59fpwTj+gbb69mJfp3C537Jfn6dwnayZH/+Btv6ogg5e2284uX47wCe6ncwT77Z/bubTvRbsWLFU/4+8UMPPWQ8Ho/517/+ZYwx5ssvv7S13/j4ePPBBx/47+/cudNYlmU++uijgHZvvvmmiY+Pt61fpwTb8xts/Tr1/Dol2NY32F7PTvXrFD73jympz69T2E4eU1I/f4NtfVF0GHE/D2XLltVf//pXjRgx4ozt3n33Xf39738v9v3u3bv3lL+r+q9//UuhoaF65JFHdOTIEbVr1862PqVjv7MZHR3tv1+mTBlJUvny5QPalS1bVgcOHLC1bycE2/MbbP069fw6JdjWN9hez0716xQ+948pqc+vU9hOHlNSP3+DbX1RdAju56FevXpav379Wa9mHhsbWyL6TUlJ0Q8//KCrrroqz7xhw4ZJkoYPH67ly5fb2m+5cuW0ZcsW/32Px6M+ffrkCe7bt28vET/dEWzPb7D169Tz65RgW99gez071a9T+Nz/fyXx+XUK28n/VxI/f4NtfVF0uDjdebjsssv07bffnrVddHS0rT+x4FS/TZo00euvv37a+cOGDdPw4cP10Ucf2danJNWpU0eLFy/237csS+PHj1elSpUC2n311VeqWbOmrX07Idie32Dr16nn1ynBtr7B9np2ql+n8LkfqKQ9v05hOxmopH3+Btv6ogg5fax+cbJ//37z22+/BU2/n332mbnpppvMzp07z9juqaeeMmlpabb1+80335gFCxactd2dd95ppkyZYlu/Tgm25zfY+nXq+XVKsK1vsL2enerXKXzun1pJeX6dwnby1ErK52+wrS+KjmWMMU7vPAAAAAAAAKfGofIAAAAAALgYwR0AAAAAABcjuAMAAAAA4GIEdwAAAAAAXIzgDgBAEVu2bJn++te/qnLlyoqIiFBiYqIaN26sQYMGBbR74YUXNGXKlAL1NWLECL399tt5pi9YsECWZWnBggUFWn5+3HrrrSpVqpSty5wyZYosy9Jvv/1m63IBAHADgjsAAEXogw8+UJMmTbRv3z6NGjVKc+bM0bPPPqumTZtq1qxZAW0LM7jXr19fS5cuVf369Qu0fAAAUPhCnS4AAIBgMmrUKFWrVk2ffPKJQkP//2O4W7duGjVqVJHVERsbqyuvvLLI+gMAAPnHiDsAAEVo9+7dKlu2bEBoP87j+f+P5apVq2rNmjVauHChLMuSZVmqWrWqJOnQoUMaNGiQ6tatq7i4OMXHx6tx48Z65513ApZnWZb+/PNPTZ061b+MtLQ0Sac/VP7dd99V48aNFRUVpZiYGLVq1UpLly4NaJORkSHLsrRmzRrdfPPNiouLU2Jiom6//XZlZWXl6/9StWpVtW/fXh9//LHq168vr9erGjVq6KWXXsrT9ssvv1TTpk0VGRmp5ORkpaen6/Dhw6dc7qxZs9S4cWNFR0erVKlSatOmjVatWuWfv3jxYoWFhWnw4MEBjzt+6P3kyZPztT4AANiJ4A4AQBFq3Lixli1bpnvvvVfLli07beCcPXu2qlevrnr16mnp0qVaunSpZs+eLUnKzs7WH3/8ocGDB+vtt9/WzJkz1axZM3Xs2FH//e9//ctYunSpvF6vrr32Wv8yXnjhhdPWNmPGDN1www2KjY3VzJkzNXnyZO3Zs0dpaWlavHhxnvadOnXSRRddpDfffFNDhw7VjBkzdN999+X7f/Ptt99q0KBBuu+++/TOO+/o0ksvVe/evbVo0SJ/mx9++EEtW7bU3r17NWXKFL344otatWqVHnvssTzLGzFihG6++WbVqlVLr732ml555RXt379fV111lX744QdJUrNmzfTYY4/p6aef1rvvvitJWrNmjfr166dbbrlFvXv3zvf6AABgGwMAAIrMrl27TLNmzYwkI8mEhYWZJk2amJEjR5r9+/cHtL3kkktM8+bNz7rMI0eOmMOHD5vevXubevXqBcyLjo42vXr1yvOY+fPnG0lm/vz5xhhjjh49apKTk02dOnXM0aNH/e32799vypcvb5o0aeKfNmzYMCPJjBo1KmCZffv2NZGRkSY3N/eM9fbq1ctER0cHTKtSpYqJjIw0Gzdu9E/z+XwmPj7e9OnTxz+ta9euxuv1mszMzID1r1GjhpFkNmzYYIwxZtOmTSY0NNT0798/oJ/9+/ebpKQk06VLF/+03Nxcc+2115rSpUub1atXm1q1apkaNWqYAwcOnHE9AAAoKoy4AwBQhBISEvT5559r+fLl+ve//60bbrhB69evV3p6uurUqaNdu3ad03Jef/11NW3aVKVKlVJoaKjCwsI0efJkrV27Nl91rVu3Tlu3blXPnj0DDtkvVaqUOnXqpC+//FIHDx4MeEyHDh0C7l966aU6dOiQduzYka8a6tatq8qVK/vvR0ZG6qKLLtLGjRv90+bPn6+WLVsqMTHRPy0kJERdu3YNWNYnn3yiI0eO6G9/+5uOHDniv0VGRqp58+YBpwhYlqX//ve/iomJUcOGDbVhwwa99tprio6Oztd6AABgN4I7AAAOaNiwoR544AG9/vrr2rp1q+677z799ttv53SBurfeektdunRRxYoVNW3aNC1dulTLly/X7bffrkOHDuWrnt27d0uSKlSokGdecnKycnNztWfPnoDpCQkJAfcjIiIkST6fL181nLy848s8cXm7d+9WUlJSnnYnT9u+fbsk6fLLL1dYWFjAbdasWXl2kCQkJKhDhw46dOiQ2rZtqzp16uRrHQAAKAxcVR4AAIeFhYVp2LBheuaZZ7R69eqztp82bZqqVaumWbNmybIs//Ts7Ox813A8NG/bti3PvK1bt8rj8ahMmTL5Xr5dEhISlJmZmWf6ydPKli0rSXrjjTdUpUqVsy537ty5Gj9+vK644grNnj1bb775pjp16mRP0QAAFBAj7gAAFKFTBWNJ/kPck5OT/dNOHm0+zrIshYeHB4T2zMzMPFeVP9MyTnbxxRerYsWKmjFjhowx/ul//vmn3nzzTf+V5p3WokULffrpp/4RdUk6evSoZs2aFdCuTZs2Cg0N1S+//KKGDRue8nbctm3bdMstt6h58+ZasmSJOnTooN69e2vDhg1Ftl4AAJwJI+4AABShNm3aqFKlSrr++utVo0YN5ebm6ptvvtHTTz+tUqVK6R//+Ie/bZ06dfTqq69q1qxZql69uiIjI1WnTh21b99eb731lvr27avOnTtr8+bN+te//qUKFSrop59+CuivTp06WrBggd577z1VqFBBMTExuvjii/PU5fF4NGrUKPXo0UPt27dXnz59lJ2drSeffFJ79+7Vv//970L/35yLhx56SO+++67+8pe/6JFHHlFUVJSef/55/fnnnwHtqlatqkcffVQPPvigfv31V7Vt21ZlypTR9u3b9dVXXyk6OlrDhw/X0aNHdfPNN8uyLM2YMUMhISGaMmWK6tatq65du2rx4sUKDw93aG0BADiG4A4AQBF66KGH9M477+iZZ57Rtm3blJ2drQoVKuiaa65Renq6atas6W87fPhwbdu2TXfccYf279+vKlWq6LffftNtt92mHTt26MUXX9RLL72k6tWra+jQofr99981fPjwgP6effZZ9evXT926ddPBgwfzXJjtRN27d1d0dLRGjhyprl27KiQkRFdeeaXmz5+vJk2aFOa/5ZzVrl1b8+bN06BBg9SrVy+VKVNGPXv2VKdOnXTnnXcGtE1PT1etWrX07LPPaubMmcrOzlZSUpIuv/xy3XXXXZKkYcOG6fPPP9fcuXP958mXKVNGr776qq6++moNGTJEY8aMKerVBAAggGVOPB4OAAAAAAC4Cue4AwAAAADgYgR3AAAAAABcjOAOAAAAAICLEdwBAAAAAHAxgjsAAAAAAC5GcAcAAAAAwMUI7gAAAAAAuBjBHQAAAAAAFyO4AwAAAADgYgR3AAAAAABcjOAOAAAAAICL/R/Aau0kXWzXQgAAAABJRU5ErkJggg==",
      "text/plain": [
       "<Figure size 1200x700 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# PLOT THE DENDROGRAM\n",
    "\n",
    "plt.figure(figsize=(12, 7))\n",
    "#plt.title(\"Hierarchical Clustering Dendrogram\")\n",
    "plt.xlabel(\"Station Index\", fontsize=12,)\n",
    "plt.ylabel(\"Distance\", fontsize=12)\n",
    "\n",
    "dendrogram(\n",
    "    Z,\n",
    "    leaf_rotation=90.,  # rotates the x axis labels\n",
    "    leaf_font_size=12., # font size for the x axis labels\n",
    "    color_threshold=0.7 * max(Z[:, 2]) # colorizes clusters\n",
    ")\n",
    "\n",
    "#plt.axhline(y=15, color='r', linestyle='--') # Example \"cut\" line\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "id": "db655b50-45f2-46f3-a3b9-64f08f7d9553",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "station ID, combined cluster, original 3\n",
      "96011 1 3\n",
      "96035 2 1\n",
      "96163 1 3\n",
      "96109 1 3\n",
      "96253 1 3\n",
      "96195 1 3\n",
      "96087 1 3\n",
      "96221 1 3\n",
      "96295 1 3\n",
      "96237 1 3\n",
      "96749 1 3\n",
      "96839 1 3\n",
      "97230 2 1\n",
      "97240 1 3\n",
      "97300 1 3\n",
      "96581 1 3\n",
      "96685 1 3\n",
      "96633 2 1\n",
      "96509 2 1\n",
      "97180 1 3\n",
      "97072 1 2\n",
      "97048 1 2\n",
      "97014 1 3\n",
      "97430 1 3\n",
      "97724 1 3\n",
      "97900 1 3\n",
      "97372 1 3\n",
      "97502 1 3\n",
      "97560 2 1\n",
      "97980 1 3\n",
      "97690 2 1\n"
     ]
    }
   ],
   "source": [
    "print(\"station ID, combined cluster, original 3\")\n",
    "for n in range(nid): \n",
    "    if labels_ref[n] == 1:\n",
    "        c = 2\n",
    "    if labels_ref[n] == 2:\n",
    "        c = 1\n",
    "    if labels_ref[n] == 3:\n",
    "        c = 1\n",
    "    \n",
    "    print(id[n], c, labels_ref[n])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4862fafa-d440-450b-b92f-1c217f9ba929",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "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.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}