sec6_physical_interpretation/dzRH_contour_plot.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "a758a875-1c32-45f6-b133-702feb45b76f",
   "metadata": {},
   "source": [
    "**Modification to always satisfy the RH-constraint**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "56e4d405-d113-4c27-965b-bb01807ad169",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import sympy as sp\n",
    "import pandas as pd\n",
    "import matplotlib\n",
    "import matplotlib.pyplot as plt\n",
    "import json\n",
    "import sys\n",
    "import os\n",
    "import gc\n",
    "\n",
    "sys.path.insert(0, '../sec2_data')\n",
    "from functions import append_dict_to_json\n",
    "from sklearn import tree\n",
    "\n",
    "np.random.seed(10)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "67f91be4-799d-4915-8f85-4a3327ce92ee",
   "metadata": {},
   "outputs": [],
   "source": [
    "def round_expr(expr):\n",
    "    d = {}\n",
    "    for n in expr.atoms(sp.Number):\n",
    "        d[n] = sp.Number('%.6g'%n)\n",
    "    return expr.xreplace(d)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "dd38bf0e-4407-4709-a661-41216c023ceb",
   "metadata": {},
   "outputs": [],
   "source": [
    "no_of_regimes = 2\n",
    "regime = 1"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "af2f9173-f118-4d30-803c-db6b13cc05f5",
   "metadata": {},
   "source": [
    "**Read data**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "9cb7228f-5656-4737-b862-c0072551e8bb",
   "metadata": {},
   "outputs": [],
   "source": [
    "all_possible_features = ['hus', 'clw', 'cli', 'ta', 'pa', 'zg', 'fr_land', 'U', 'rh', 'ps', 'hus_z', 'hus_zz', 'clw_z', 'clw_zz', 'cli_z',\\\n",
    "            'cli_zz', 'ta_z', 'ta_zz', 'pa_z', 'pa_zz', 'U_z', 'U_zz', 'rh_z', 'rh_zz']\n",
    "\n",
    "loc_all = {}\n",
    "for i in range(len(all_possible_features)):\n",
    "    loc_all[all_possible_features[i]] = i\n",
    "    \n",
    "# Features\n",
    "features = ['rh', 'ta', 'clw', 'cli', 'rh_z']\n",
    "no_features = len(features)\n",
    "\n",
    "loc = {}\n",
    "for i in range(len(features)):\n",
    "    loc[features[i]] = i"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "b8462c5d-87c6-4d8d-983f-7daf8582855a",
   "metadata": {},
   "outputs": [],
   "source": [
    "path_data = os.path.join(os.environ['HOME'] + '/my_work/icon-ml_data/cloud_cover_parameterization/neighborhood_based_SR_DYAMOND')\n",
    "\n",
    "# Load the input data and pick the five best features (rh, ta, clw, cli, rh_z)\n",
    "input_data = np.load(path_data + '/cloud_cover_input_dyamond.npy')\n",
    "input_data = np.concatenate([np.expand_dims(input_data[:, loc_all[sel_var]], axis=1) for sel_var in features], axis = 1)\n",
    "\n",
    "output_data = np.load(path_data + '/cloud_area_output_dyamond.npy')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "6831a5fb-6fcb-4785-86a9-700a0fe82730",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(285179494, 5)"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "(samples_total, no_of_features) = input_data.shape\n",
    "(samples_total, no_of_features)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "9b0eaa09-11a3-485e-b9ff-6befe674d5b0",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "# Construct training and validation data\n",
    "training_folds = []\n",
    "validation_folds = []\n",
    "two_week_incr = samples_total//6\n",
    "\n",
    "for i in range(3):\n",
    "    # Note that this is a temporal split since time was the first dimension in the original tensor\n",
    "    first_incr = np.arange(samples_total//6*i, samples_total//6*(i+1))\n",
    "    second_incr = np.arange(samples_total//6*(i+3), samples_total//6*(i+4))\n",
    "\n",
    "    validation_folds.append(np.append(first_incr, second_incr))\n",
    "    training_folds.append(np.arange(samples_total))\n",
    "    training_folds[i] = np.delete(training_folds[i], validation_folds[i])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "ff49c846-d3de-4b64-b8d5-42f1db80e37d",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "844"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# The second fold yields the best model\n",
    "flattened_input_train = input_data[training_folds[1]]\n",
    "flattened_input_valid = input_data[validation_folds[1]]\n",
    "flattened_output_train = output_data[training_folds[1]]\n",
    "flattened_output_valid = output_data[validation_folds[1]]\n",
    "    \n",
    "# Remove input_data, output_data\n",
    "del input_data, output_data, training_folds, validation_folds\n",
    "gc.collect()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "18a94f27-da39-4411-9bf0-d139e8fc4c4c",
   "metadata": {},
   "outputs": [],
   "source": [
    "if no_of_regimes > 1:\n",
    "    # Already remove the regime with clw + cli = 0\n",
    "    reg_not_0_train = np.where(flattened_input_train[:, loc['clw']] + flattened_input_train[:, loc['cli']] > 1e-20)[0]\n",
    "    flattened_input_train = flattened_input_train[reg_not_0_train]\n",
    "    flattened_output_train = flattened_output_train[reg_not_0_train]\n",
    "\n",
    "    reg_not_0_valid = np.where(flattened_input_valid[:, loc['clw']] + flattened_input_valid[:, loc['cli']] > 1e-20)[0]\n",
    "    flattened_input_valid = flattened_input_valid[reg_not_0_valid]\n",
    "    flattened_output_valid = flattened_output_valid[reg_not_0_valid]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "dd5711b1-f92c-47db-bef5-f162a9559355",
   "metadata": {},
   "source": [
    "**Normalize the features**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "f7d119dc-7b3a-4fa0-974c-a1b755b7b40b",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "# Scale the data\n",
    "mean_all = [4.12205844e-03,2.25493498e-05,3.38180032e-06,2.57065512e+02,6.00030443e+04,5.64080139e+03,2.35046400e-01,1.32776682e+01,6.02512234e-01,9.86270417e+04,-1.27545273e-06,-4.02484958e-10,1.65204582e-08,-4.34660202e-11,4.29441131e-10,-1.82817316e-12,-4.68742483e-03,-7.54899040e-07,-7.51544542e+00,-1.06989723e-04,1.65615172e-03,-9.27604679e-06,-4.76200071e-05,-1.32246548e-07]\n",
    "std_all = [5.07648249e-03,5.69702638e-05,1.01308124e-05,3.00533874e+01,3.12514292e+04,5.66963918e+03,4.11184302e-01,1.11389888e+01,3.32494615e-01,6.24039256e+03,2.03179260e-06,1.17041141e-08,1.33311867e-07,1.42840744e-09,6.73384546e-09,5.07424672e-11,5.82875686e-03,6.34826092e-05,3.53136052e+00,1.13215264e-02,6.62892130e-03,6.08144307e-05,2.58065098e-04,2.49552692e-06]\n",
    "\n",
    "mean = np.concatenate([np.expand_dims(mean_all[loc_all[sel_var]], axis=0) for sel_var in features], axis = 0)\n",
    "std = np.concatenate([np.expand_dims(std_all[loc_all[sel_var]], axis=0) for sel_var in features], axis = 0)\n",
    "\n",
    "# Work with scaled training folds\n",
    "train_data_scaled = (flattened_input_train - mean)/std\n",
    "valid_data_scaled = (flattened_input_valid - mean)/std"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7886e0e1-49bd-4c87-bdb2-09d260e12710",
   "metadata": {},
   "source": [
    "**The condition for the artificial increase of RH**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "3fa28bf5-52b9-4b97-afad-8e94d26782f4",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "If the following term is greater than 0, then RH is increased:\n"
     ]
    },
    {
     "data": {
      "text/latex": [
       "$\\displaystyle - 3.00757 rh - 32.192 \\left(0.00389006 ta - 1\\right)^{2} + 0.952388$"
      ],
      "text/plain": [
       "-3.00757*rh - 32.192*(0.00389006*ta - 1)**2 + 0.952388"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "x0, x1 = sp.symbols('x0 x1')\n",
    "rh, ta = sp.symbols('rh ta')\n",
    "\n",
    "(a,b,c,d) = (38.6562122, 43.53500518, 19.78403208, 1.13637902)\n",
    "\n",
    "cond = 1/(2*c*d)*(-c*x1**2-a) - x0\n",
    "\n",
    "X0 = (rh - mean_all[loc_all['rh']])/std_all[loc_all['rh']]\n",
    "X1 = (ta - mean_all[loc_all['ta']])/std_all[loc_all['ta']]\n",
    "\n",
    "subs_set = []\n",
    "for k in range(2):\n",
    "    subs_set.append((locals()['x%d'%k], globals()['X%d'%k]))\n",
    "subs_set = set(subs_set)\n",
    "cond_orig_inputs = cond.subs(subs_set)\n",
    "\n",
    "print('If the following term is greater than 0, then RH is increased:')\n",
    "round_expr(cond_orig_inputs)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "id": "3a9cde1c-09ee-4204-ab7b-8af76b6a2144",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "That means whenever RH < 0.3269 - 10.5359*(0.00389006*T-1)^2\n"
     ]
    }
   ],
   "source": [
    "print('That means whenever RH < 0.3269 - 10.5359*(0.00389006*T-1)^2')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "id": "c17d0c8d-8b35-438f-8896-80355bf2d8c1",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle - 0.000159435210586049 T^{2} + 0.081970566308 T - 10.209$"
      ],
      "text/plain": [
       "-0.000159435210586049*T**2 + 0.081970566308*T - 10.209"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "sp.expand('0.3269 - 10.5359*(0.00389006*T-1)**2')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "id": "8dd2d168-726a-45b7-831e-340d47e982e9",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "3.2459051125995657"
      ]
     },
     "execution_count": 15,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.sqrt(10.5359)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "id": "4f7ff279-8f32-4ed8-b26f-9de93018ddd5",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.012626765642319066"
      ]
     },
     "execution_count": 16,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.sqrt(10.5359)*0.00389006"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "id": "f5341966-6833-4bff-a763-1496180ca7f8",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Increase the general font size\n",
    "matplotlib.rcParams['legend.fontsize'] = 16\n",
    "matplotlib.rcParams['axes.labelsize'] = 16 # For an axes xlabel and ylabel\n",
    "matplotlib.rcParams['axes.titlesize'] = 17\n",
    "matplotlib.rcParams['xtick.labelsize'] = 16\n",
    "matplotlib.rcParams['ytick.labelsize'] = 16"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "id": "c0d3d1e4-dd98-4ab1-b977-4e54c1e488e7",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZgAAAFdCAYAAAA+HntSAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8qNh9FAAAACXBIWXMAAAsTAAALEwEAmpwYAABVC0lEQVR4nO2deZycVZW/n2/1kj0kISsJECDsooCRRVSCoKCg4s64gRvIuM24DoP7xqjjz2VckFEGVNRxxAV3dhTZFRWQxYQkEJaQjexLd9f5/XHf7lSq79v9vl31VlV3nyef91NV973LudWVOnXvOfccmRmO4ziOU29KzRbAcRzHGZm4gnEcx3EKwRWM4ziOUwiuYBzHcZxCcAXjOI7jFIIrGMdxHKcQXME4juM4heAKxnEcxykEVzDOiELS+yQtk9TZbFmyIOntku6VtE1SWdKEZsvkOPXCFYwz0rgQmASc02xBBkPSPwFfBW4iyPsqM9vcXKkcp37IQ8U4Iw1JXwJeCuxjZuWMbSYA04HHzGxHgeJVjnk1MMvMDmvEeI7TaHwF44xEfgjsBTw3S2VJ7wE2AcuAVZLeV5xofWMKOBq4vuixHKdZuIJxhgWS9pb01cResVnSBklXS3pmpPqtwHrgZRm7/xkwC5gKfAn4vKQz6iF3DElfA8rAROAdkkzSE0WN5zjNwrfInGGBpFcAHycog+WE7aw3A3OBhWZ2d1X9IW8/SboF6DazZ9Uqd0r/xxOU37uAjwBLgbVm9usixnOcZuEKxhkWSBpvZluqyqYB9wE/N7O3Vt27CHgT0GE5P+SSvgi8ycx2q1HsgcY4D/gMMMXM1hc1juM0E98ic4YFlcpF0jhJuxM+v7cBT480WQu0AVMG6jfp638kbZL0kKSzgW6go27Cx3kasMyVizOScQXjDAskdUr6tKSHgC3AamAVcCpxJaKMXX+asF31auAM4G3Ay4FHK8aenthJNknaImmppNdE7s+pkvn7ki5IGfdpwF8zyug4wxJXMM5w4cvAvxFsMGcAJwPPA64l/jmeSliJPJnWoaTxwNnAV8zsV2Z2E/AKYE/g3oqqRwCPmNlEMxsPXABcLGlsxf11ZvZY1RCHAn+JjDsW2B/4W/p0HWf4095sARwnI2cA3zGzd1UWSvpESv39gPsGsb8sACYQlFQvbYT/F7+rKDscuLPi9Q3AGMKBzm3J/Xuq5GoDDiSiYICnJOP4CsYZ0fgKxhkulKn6vEp6NnBMdcXkjMnTgRsz9AnhDEwv7wU2AN+vKDsC+HPS9+7Ap4AbzGxVxf1dFAxBefUA/4iM+7Tk0RWMM6LxFYwzXPg5cKakTYRVwcHAWwhf7JOq6h4N7Ab8dJA+7wceAt4q6V7gjYQts7PNbG1FvcOBFycHMicStulOr7o/P3Gl7qUd+FtKJIGnEZTakkHkc5xhjSsYZ7jwbsJ21MsIiuAugr3kNcCiqrqvBh4Grh6oQzPrkvRy4GLCqmUZ8Foz+0FvHUnjgAMIW24PAa8F/jORpfL+S9l1G+080ncIngbcndd92nGGG34OxhlRJDHFHgI+bmZfydGuw8y6IuVHA1dWnomRtAR4v5n9JLl/IzDZzLZW1LkK+LGZfbOG6ThOXZF0HPBRwqp7LLAY+KqZXVzEeG6DcUYabwM2EqIqZyamXBKOIKyWKvkN8OKK+3dXKRcBC4kb+B2nKUh6KmFV3wG8leCOfzvwbUnnFjKmr2AcJx1JFwJlM/vnirLTCNtqs4GvAyUzO7vi/v6ECAOTqqMPOE6zkPQZ4H3ANDPbVFF+C2Bmdmy9x2z4CkbSPEn/Jenm5NCaSZpfVeeSpDx23TdA3+cldQbzHnKcTJjZ2yqVS1L2SzObaWbl5P7ZVff/YWZtrlycFqMT6AK2VpU/SUG6oBlbZAuAVwHrgD+k1PkkcGzV9U/JvStiDSTtC5wPeFRax3Gc/lySPH5F0h6Spkh6K3Ai8MUiBmz4FpmkUq/rpqS3AP9NSAy1bJB2HwY+ATzFzKrPHCDpdwQvoAOB9qIi4TqO4wxXJD2D4L4/NynqAs41s28XMV7D3ZSzZhiM8AbgTynK5TXAkYRVzk/ydDpl2m42e89Zu8qYUjeujOO1ozVTlLlFasfK0vqItk8bK9o+TqxuOc8cMsqaWjfXHPL0G62a+b3JJVeesVL7jZZmrptrrOxDRQdL/72asW6esVKHyiFXnh/Y1VW3dWM7yllj3uVC08caO3J+VW7suofEfT7hIjO7qK/PYBu8nHB27G2ErbKXABdK2mZml9UqdzXD4hxM4lq3gJA/o/reVMLy7gNmtjY48GRn9p6zuOi3X9+lrJyiA7vL3f3LrH9Zat1IWVofaXV39PR3duop92SqB7Cju3/5jp6UsWJ1I2Vp48Xq9lh/WdPlyj6HrpQ5dHX1L+/qjssQK+8p9/8sdHWlta9trO60upE59JTjX4xZZehOmUO5u/98Y2WhvH8fqXUj8hY2Vq5++5dbynvbr/zWAnfju8ro2FmD16vArlyxzcwWDlDlM4QVy2kVXpPXJNEpvizpBzUsAKIMFzflNxDemB9E7n0eeICd+4uO4zjDHpWU68rAYcBfIy75twG7AzPrPYeWX8FIGkNwCvilma2uuvdsgvI5Ms+p6CTnx9kAs+bW/T11HMepkcxKIw+PA4dL6jSzHRXlRxO21tbGmw2d4bCCeQkh38elkXvfBL4NrEg8IqYQlGZb8npMrEMzu8jMFprZwim7F5a00HEcZ2iokBXMV4F9gF9Ieomk50v6KsF2/Y0qpVMXWn4FA5xJSC4Vy1d+cHK9LXJvHfCvwJcG6lwqMbZt7C5lafaPkvrr41I5rqNLEd3dXoq/3VF7jeIyxPqItW8rtUXbt6l/eVspbuvobOs/1o72eKLHqA2lLbtdJS5X9jl0tqfYdiIy1GqvaSvF/+Ztbf3/k6fV7WjvP4c0e017pG6avaYt8kUTk6ErRa6utv7vQbk9Xrfc3b88lw0lImtPyhdl7AtUpfhYpfbabDBZ6xZj3g+I+Jxrwcx+LOmFwAeBbxFCxSwB3k74sV53WlrBSJoFPB/4ekoojxMiZV8i5Np4JyHOjuM4zvBCoJQfAbVgZr8hhDpqCE1RMBVhzXtzqb9A0ipglZndUFH1tQQZY9tjmNn1kb6fJJyD6XfPcRxneFCIDabhNGsF839Vr3v9hG9g19DrZxICCf65EUI5juO0BKr/FlkzaIqCMbNM75yZPW3wWv3aLMotkOM4jlN3WtoG4ziOMxoRUPIVzPBHiM5S5y5lMW8xiHuMtSvFM6yU/XR+zOMsjwx52ne29fcCy3NiPs3jrC0yXsx7qS3ieZRWt7Oc4rFWs3davG7USzDyd4x5gAF0dUc8w1K822In7lO906JeYPEvn9iXUkdH/89ozGMuyJXduy3mcZbPC6x/v6VUj7X+nl09OTzD0jzOYjKkbU1Vl1vOqCF58S0yx3Ecp/64DcZxHMcpBvcicxzHcYpAoMjB3eGGKxjHcZwWo4iT/M1g1CuYkkRn265G/nZLMdxHwrekhevPaowPMuQIQROpmzV8TFp5mkNArUb6WN2YkwFkD/dflFxpdTsjoXHS5Io5BJRSDPcxY3x63ZjzQfzLp6M9YtCPOBSktY+PlX0OXRGjOUC5I2Kk78oeaiaPMT7maJDVcA9Qao/Hza0Od9NdpAJwG4zjOI5TDCokVEyjcQXjOI7TgvgKxnEcx6k/vkXmOI7jFIEb+R3HcZxi8BXMyECU+oWKKZOSyKhGb69c4V9S6sY83PJ4htU6h1joFYgnB8sTlibm2ZWacCwSziSPXGl1s3qM5Qnp0hYJpxLqZk9OFivv6EgJV9MVCb9So8daW1uaXNnn0FPu/38q5nHWneKFlifhWMzjLJbcDOIeZ1aOe5FVy6BCQ8X4QUvHcRynIFzBOI7jOPXHt8gcx3GcIhgpRv7hf5LHcRzHaUl8BQOUqoy+srjejYaVqDFvS1ofaYb3qEE/YoRNNdwXkHsmrY+y9TfCphruc+RtyZp7BtKcB7IbouNy1T6HrLlnQnnkvY3IGsaLGelrcyjIk3smzXkgJm98rHj7InLPpNVNM/JXj1ewjd9XMI7jOE4RiFKplOvK3LP0Qkm/l7RJ0gZJd0h6bhGz8BWM4zhOiyHibuA19yudA3w1uT5JWGQcDoyv+2C4gnEcx2k9lL7dOOQupfnAl4D3m9mXKm79rq4DVeAKxnEcp8UQSrWT1cCbgDJwYb07TsNtMI7jOK2GoNSmXFcGngXcB5whaYmkbkmLJb29qGmM+hWMEO1Vb0OPUkJQEPFOiZRBimdXHULFxMqjnmVpoUAaGO4mloytkcnNQt2IZ1ck1AxkDwtTVHKzVM+uiAdVV1d2j7OikpvlSU7W3R0LYVNMcrNYuJn0hGP9/69XJxZL7aNAL69gg6n77/89kuvzwL8DS4BXAl+V1G5mX673gKNewTiO47Qc0lBsMNMl3VHx+iIzu6jidQmYBJxlZj9Jyq5NbDPnSfqKmcV9tIeIKxjHcZwWY4heZKvNbOEA99cA+wNXVZVfCZwCzAEezTvoQLiCcRzHaUHq7UUG3AMcEynv1WRx20ANuJHfcRynxZAKOWj50+Tx5Kryk4EVZvZ4XSeBr2CQRJt2zQMi4ka+skXKU1axsXAzaYbGXIb3jKFeYoZ/IPobJeaQkCZDreFuUuWKdNuWkvMkT46XWKiWWKiZULd/eczwn0Y+uWqbQ5qzRFd3yvvbr888OW3in9s8jgqxcDNdKcb0+FjZw93Ews3EQs1APNxM2v/TauN/0YFcCjho+WvgOuCbkqYDDwKvAJ4PvLHeg4ErGMdxnJZDBRy0NDOTdDpwAfBxYCrBbfm1Zvb9ug6W4ArGcRyn5SjkoCVmtgF4e3IVjisYx3GcFkPJQcvhjisYx3GcFqSIFUyjGf4zcBzHcVoSX8Eg2lUVKibFPyQaKsbidaPhZtLOyEa6SPXsyuitVWuomVA5XhytGvGaKyLUTF5yhWSJJierLdRMHuqRNC3uBZY91Eweag31UmuomTykhbuJhZuJhZqBRoeKGdJJ/pbDFYzjOE6roWLywTQaVzCO4zgthq9gHMdxnGLwFYzjOI5TBKKQWGQNZ9QrGAGl6hAdOYzxSrGEx8LNlFOcB2JOBdGwNEB7KeOfLCVsXdSwmidPTcp8owb5mAwp/2eKCJcD0FWuzfAeI0+omTz5YGp1Ekgj5qiQ+jcvZc8tlDUsTR7yhJrJk6cmT1iaNKpnW+j6QsUctGw0o17BOI7jtBpiZBy0bKiKlPQKSZdLWi5pq6T7JV0gaVJVvamSviVptaTNkq6WdFikv89IulLSGkkm6ayGTcZxHKcoFFZzea5WpNFSvQ/oIaTrPAX4BnAucJUU1uySBFyR3H8n8HKgA7hO0ryq/t4JjAN+2RDpHcdxGkCvF1mdw/U3nEZvkb3IzFZVvL5B0lrgUmARcC3wYuBZwHPN7DoASTcDS4EPAO+qaL+bmZUlLQDe0AD5HcdxGoJ7keWkSrn0cnvyODd5fDHwaK9ySdqtl/QL4CVUKBgzqzkDmxDtVW9Dj1JO8kdO7ZeVYjyMOQqkfV4idZVihFXE+B/LX5FmxI2e2k97F4s4yZ9iuI/KledHWY451JpnZkdPMcb4NGrNM1OU80B8/NryzHR11W6Mz0MteWbS8sbUgyLC9TeDVpjB8cnjvcnjocDdkXr3AHtJmtgQqRzHcZqGRoQNpqleZJLmAp8ArjazO5LiacCySPW1yeNUYFPx0jmO4zQHKb4aHW40TcEkK5GfE9zLK9N1inwbTEMZ+2zgbIB5e1b7DTiO4zQXodTzVsOJpigYSWMJnmL7Aseb2YqK22sJq5hqpiaP62od38wuAi4COOLIw9OOVTqO4zQNX8EMAUkdwOXAUcBJZnZXVZV7gOdHmh4CPGRmvj3mOM6IRn6SPz/JWZfLgBOBU83slki1K4A3SjrezG5I2k0GXgR8vwCpagoVky+sTJ48MyneNFllSNtQLMAzDHJ4pxXkGdboPDO10sg8M7WS6hlWUJ6ZWmlUnhmleJs6O2n0CuZrwCuBTwObJR1TcW9FslV2BXAz8D1J7ydsiZ1H+Mr8XGVnko4HZgCzk6KFkjYBmNmPi5yI4zhOkYyELbJGr8FekDyeT1AilddboO9sy2nAVcDXgZ8STv+fYGYPV/X3ceD/gP9KXr89ef1/xU3BcRynWHqN/HmuIY0j/TYJs/WpOk8BaPxBy/kZ660F3pRcA9VbVLtUjuM4rUUj3JQl/RPwtCLHGP5WJMdxnBFHsQctJU0Bvgi8pwjpexn14foVM/KnkcOYHssTk8twX6sMOZwPSm0NDCtTUFiaXIb7OsiQlbT/+C1ruI/kr2l0aJwYeQz3ReSpiZHmtFOfvgtfwXwOuMfMfiCpAOepwKhXMI7jOC2HijtoKelZhODAhW6PgSsYx3GclqOoFUxyDvGbwH+a2f11H6AKVzCO4zgthhjSQcvpku6oeH1RErWkkg8Scmh9uhb5suIKxnEcp9UYmhfZajNbmNqltBfhiMhbgDGSxlTcHpMY/jeapRmL8+NeZI7jOC1GQedg9gXGAt8jHGDvvSBkG14H9EtNXwu+ggHaqvVsrV5dgEWziOXoowjPspS6luJW1V7K8fEoICxManiPmMdYPTzDapxDbEuj2d5iafSkJH5rVUrl1gvLUnSkmAJsMH8BToiUX0dQOt8GFtdzQFcwjuM4LUYRwS7N7Eng+thYwHIz63evVlIVjKSHhtCfEYJYxjJSOo7jOBkIXmTDa5UZY6AVzDzg18CqjH2VgNcBnbUK5TiOM6qRGhbs0swK2+wbbIvsE2Z2W5aOJLUDr69dJMdxnNGNYMRntPwwUB29OBUz65b0YeCRmqVqMKLE5u4d3LZ6OYtmLUj/w+bQ81HDeUF5ZlIEyNze0jqN9NFIw380VA3QHvvYpvVbRE6aHGN1tqiVc0dKeatuP8Ty+jSbYvPBNG4FUySpH38zy3wQR9JEM9uUp00rsbW7iwsf+CMPb3mSVds38/K9DqPkyYQcx3FqIvVngaQ3Z+kgOZxzVb0EagbXr1zMw1ueBOCmVUv56cN3YZa2BHAcxymWEK6/uGjKjWKgBfxFkraZ2WVpFSRNIyiXg+ouWQN5/h4Hsnr7Jv60dgUANz6xlA618aJ5h3haVMdxGo5G+hYZ4eDN/yRK5vLqm5JmANcQToe+pCD5GkKbSrx2n4WUzbhzXTAhXbdyMR2lEi+Ye3CTpXMcZzQy0o38ZxFsft+X9Aoz+0XvDUlzCMplHvBCM/t9oVI2gJLE6/ZZSLeVuevJxwC48rEHaC+18bw5BzRZOsdxRhNqoJtykQxk5DdJrwN+BPxI0kvM7MokYNo1wAzgZDO7uUGyFob1BHtLCfGGvRdycfk27t2wEoBfP3IvHWrjhNn7E/1zp+yglaORYuJ2nbZYzw0MK5MW2S4mb9pvqszeZXUI6RLzLot6lqX128ikaSl1h5N3WUt4lnX1/5uX0+wODXpvi94+L42AFcyAM0iiar4auBr4qaSzgN8D04ATR4Jyqaa91MYbFxzNAZNm9JX9fMXd3LxqaROlchxnNBEy7ZZyXa3IoFKZWTfwcuAPhGBo44BFZvangmVrGp2lNt684Bj2nbh7X9mPlv+Fu9Y91kSpHMcZTYwEBTNQLLLvVBVtJGzGPAi8v2p5aGZ2Zv3Fax5j2to5e/9j+dr9N/Lwlicx4LsP3sE5BzyT/SbtPmh7x3GcoSJGxhbZQLuVz6H/jv1yYHZyVTIiD42Mbevg7P2fyZfvu4HV2zfTZWW+vfhW3nnQs5gzbnKzxXMcZ6QijWwFY2bzGyhH87CdRv4YE0udnLPfM/nKA79nY/d2tvZ08c0HbubdBx3PtDHj44Z/yGV4L0cSyCmlg5IiIxaRTyalj5jzQqja/0YrhJUpRTpJ+4+bKwRNjBrDyrSV4h4BbeX+5Y10EkgLKxOTK20ODXUUiKULKqV8cDuGPkzRR+RGgoIZ/jNoANPHTODs/Y5lTPKFub5rGxf+449s7t7eZMkcxxmJCFHK+a8VaU2pWpB546fw5gXH9B1+emLbJi5efBvdkV9xjuM4tTISjPwDxSLrkXRU1o4ktSVtjqyPaK3HAZNn8Lp9nt73esmm1fz4ob963DLHcepKr5F/uCuYgXZyBcxODlZm7WvEB+46Yto8Vm/fzK8e+TsAt65+iFljJ3HC7AVNlsxxnJHDCDfyJ/w0Z3/D8qd8vxVImiG7HPTnSTMP4ImtG7l9bUiX84sV9zBz7CSeMmVOX908p/7jQsWLeyLn7mMOAaW0wWJOAnlkyBFhoByRNe0/Taw81+n8NHIY3mOOAiXL4SQQc2qoQ+SCGDviPg1RI3vcSaAG6/YQiDkK5HESKIyYQ0BbXIZyVU6aNEecuqCRYeQfSMG8cYh9jvgj75J49V5HsGb7Fh7cvAYDvvPg7bz7oOOZO363ZovnOM4wp9fIP9wZyE350kYKMtxoL7Xxpn2P5ov3X8+aHVvYUe7hW4tv5r0Hn8DEjjHNFs9xnGHOSFjBDP8ZNJGJHWN46/7HMrYt6Ol1O7by3aV3UHajv+M4jiuYWpk9bjKv32dh3+v7NzzBbx69t4kSOY4z3CnCi0zSKyRdLmm5pK2S7pd0gaRJRc3DFUwdOHTKHJ4/58C+11c9dj93P+mBMR3HGSqFRFN+HyE7x78DpwDfAM4FrpKK2Y9r0awUDWSQUDGVlCMuQb2GuFNmH8zyTeu4f+MTAFy29E7ee/BuzBg7saJufByLeKOk5Zpos4y5Y1Kn1N+zy1J+Z0iRfDAp238WkbeUVVaI/tRJDf+Sw7MrVh5rD7XntOmOuCSlypVjDp1t9ff42tHdFS2P5XZPy/feFplDUd5pMS+0npRDzqWIvKWUUDEdNXwDDsNQMS8ys1UVr2+QtBa4FFgEXFvvAX0FUydKEm/YZyFTO8cBsK2ni4uX3MqOnhR/UsdxnBSk+q9gqpRLL7cnj3PrOoGEXApGUknSUyQdL2lCEQINZya0j+GN+xzdF07msa0buGLF3U2WynGc4UiDYpEdnzwWYjjOLJWktwOPA38lLKUOTMp/JuldRQg3HNlrwlReNu+pfa9vXLXU7TGO4+SiEaFiJM0FPgFcbWZ31HsOkFHBSHor8GXgZ4QUypW7j38gZLx0Ep45fT5PrTjV/4Nlf2b9jq1NlMhxnOHFkLbIpku6o+I6O7V3aSLwc0Isg6Eeqh+UrCau9wBfMLMPSv1ijdwHvL++YjUOA/rZrVNDxUTaRwzhAK+efwTL71nH+q5tbO7ewWVL/8TbDjiOUswymMNIX1Yk/ErE2JomlyJG+pgxP/Sbw0gf6aInIqssLc9NbYb7eoSViYaKySNXntwzOcLKxJwHUg3vMSN95PORbrjv/zfPk3tmBzmcB3LJVYxDQVfERFpKCVfTz3mgYCv/EFYlq81s4WCVJI0FrgD2BY43sxVDEC8TWWewD/C7lHubgSl1kWYEMaF9DK/bZ2Hf9/EDG1dx/cp/NFUmx3GGByrGTRlJHcDlwFHAC83sriLnkVXBrAbmp9w7EHikLtKMMPafPIPnzj6g7/WvHvk7j25Z30SJHMcZFqiQg5Yl4DLgROAlZnZL0dPIqmB+AXxE0r4VZSZpOvCvBNvMoGQ5SSrpREnfk7QkqbNE0jckzYz0t4+kH0t6UtJmSddJGnSJ2EheuMfB7DV+KgA9Zvxg2Z/pMU9S5jhOOqIQL7KvAa8E/hPYLOmYimteEfPIqmA+BGwH7gauJuy4f4Xg2tZD8ETIQpaTpG8Ddgc+ldS5AHgxcEtimAJA0u7AjcBTgHOAM5Jb10k6OKM8hdNWKvGafY7sc11+eMuT3LBycZOlchyntSlki+wFyeP5wM1V11uKmEUm852ZrUlWBv8CnAwsSdp+FfiimW3IOF6Wk6T/HKnzAHAD8Crg4qT8XGAWwUi1GEDStcCDwMeTui3B7HGTOWWPg/qSlP36kXs5dMocZo0tLASQ4zjDnHqf5Dez+XXtMAODKhhJncBnge+b2SeBTw51sCwnSXOcNj0G+Eevcknabpb0B+A0Se1mluEYvUG5ygUqxTkk6gDV31EqygkzFvDXdY+wYst6uq3M/y67k3cc+GxKUjRxUVpooKgXWMSFKy0ZUimWcCwtuVkOL7DoHCIL5JisaeWZQ7cwQFiZHJ5dUY+xPJ5hsf9OeZKb5fBOa0tJHJc94VhcrhixMC15xspTNy3hWGdUruwea+W2+HtbjshQTvFkK1eFmynSh0wUnNCsQQyqIs1sB2ELalxBMmQ5SRqr00P8s7+dIOt+tYtWP9pU4oz5R/a5KT+4aQ1/XDXic7M5jjMkhFTKdbUiWaW6Ezis3oNnOUmaOAB8iaBcflZx635g/8QW01u3RHC/A5g2wLhn9x5GWrNmTU1zyMO88VM4qcqrbGPX9oaN7zjO8KGEcl2tSFYF817gfZJOU1qY35xkOUkqqR34AWFr7IyqLa8LCfJ/R9J+kuYQHA/2Se6numqZ2UVmttDMFu6+++5p1Qrh+XMOZMaY4KuwraeLXz5yT0PHdxzHaRRZFcz/ETy7fg5sk/SwpIcqruV5Bq06SXpy7CRpshq5FDgJON3M/lZ538weBF4LPB1YDDwKHAt8ManSkgHA2kttvGyvnbHKbl29nOWb1jZRIsdxWpF865fW3CLLauq7hgEyjOSh6iTpSQOcJL2QEPfsFWZ2TayCmV0u6WfAAcAOM1si6RvAw2b2UCaBDKzacJ5iyI6/AylvS6RY5dDvwZNm8ZTd5nD3+qADL3/oLv714EW7hJGpNXdMNG8M0BPxSkh1CIiU15o7JpY3JtTtL285xYMiV/iWAnLHpDkfRA33KblnippDzMAdC6myI8X9JZ/hvrZQL1FZiYd/iRlb8zgU1CN3TL8QMgXuSom4889wI6ub8ln1GKzqJOmpaSdJJX2B4Jd9ppn9bBDZekiM/5L2ICilz9dD3iJ56bzDuG/DSrqtzMNbnuTW1cs5dsb8ZovlOE4roPSkg8OJRme07D1J+mmSk6QV91aY2QpJHyQE17wY+EdVnVVmtgT6VkKfI5yP2QAcCpwH3AN8ofCZ1MjuYyZw4qwD+N3j9wHwq0fu4YhpcxlbQPZCx3GGH6267ZWHTApG0hsGq2Nm38nQVeVJ0vOr7n0c+FhFnTclVyWXAmf1DgnsD7yGEGxzBUEpfSZxrW55Tpy9P7etXc66HVvZ1L2D61cu4ZQ9Dmq2WI7jtACjZosMuCSlvHKzclAFk+UkqZktyiJQ4lF2Wpa6rUpnqZ1T5hzMD5b/GYDrHv8Hz5qxDxM7xjRZMsdxmolQLVkqW4Y84fqrr4WEVcc/gKMLkW4U8IxpezE7CRmzvdzNVY/f32SJHMdpBSTlulqRrEb+mBvycuDPybmY9xC2qoYl1c5Oqg4d03cj0raUEjolY3IyAS+cewgXL7kVgBufWMqimQuYMmZsvN+I91AsfExUAIpLThYLIZM1fEy40V/eopKTpf0yzNxvysejqORkeeYQDSETqVqrt1eo23+stFAvecLC1C5X9rqxEDKx8DHQP4RM0VtYI8EGU48Z/AE4tQ79jFoOmzKHvSf0hvQv89vH7muyRI7jNBfl/teK1EPBHANsqkM/oxZJnDb30L7Xt695iHXbtzRRIsdxms2o2SKT9JFIcSchF8uphLD9Tg3sP3kG+03cnSWb1lA24/qVS3jpXnUP/+Y4zjAgRFMe/ltkWb3IPhYp206ww3yakBTMqZETZx/AksU3A3DL6uU8f48DmdAeC1TuOM5Ip1UDWOYhq5F/+KvSgagy6qdFiondSK0aM7yn+Q5UhJCZM3Yyj23bwI5yDzc+sZRT9tg1OWcshEwsfExaXJ+44T17jpe0D30shEzW8DEQT6uT9gsuljum1vAvqXUjZeWUlNfNDmGTVrfH+r+7baWUfDI5jPHxutlD2MSM8T3leHigmPNCPRwKYiFkYuFjQnldomVlRCNiBZNpBpKeU5muuOreREnPqa9YoxNJu4Tz/8MTD7K9J0PONMdxRhwjwQaTVUVeBxyScu/A5L5TBw6fOpdpneMB2Ny9g9vXZIvZ6TjOyKE3o+Vo8SIbSPoxZE4c7AxGm0ocP3NB3+ubVi2Lb7c5jjOC0cgO1y9pPiFfSy8LI9tk4wjxwvxndh15xrQ9+eWj99BV7uHRretZvnkd8yemJuh0HMdpSQYy8p8JfJRgLzbgv9h1JWPJ627g7UUJOBoZ397JkdPmcevqEEDhj6sedAXjOKOMIrzIJO1JSMr4PML399XAv2TOn5WTgRTMJcD1iRDXEpTI36vqbAceMLNhnZKx3w5UqrdXLCRLSqdtkRspkTFioVqeOX1+n4L5y9pHOH3PpzKhvTPuBZY1fAxgOXYz8/Qb8xjLGj4mlEdC2OSQNZcHVo3hW9Lap3mXZZar1hA2KX3EPLB60jOK92+fyzMsxTstYyKzPN5eaWT1WEurGwsfA/1DyBRpVxfx/3819SmNJ3yXbycsIAz4FHCdpKea2ea6DsgACiaJP7Y8EewE4E9m5if2G8ReE6Yyb/xurNiyni4rc/uah1g0a8HgDR3HGQEUYrh/K8HscaCZLQaQ9DdCwOJzgP9X7wEzqUgzu8GVS2ORxDNn7NP3+k9rHm6iNI7jNJoCvMheDNzSq1wAzGwp8EfgJUXMIfMaTNLJkn4q6e+SHqy6lhQh3Gjn8KlzaUvW4Q9veZLV21zHO86oQGGLLM+VgUOBuyPl95B+DKUmsh60fCHwa2A8cBBwH8FzbE+CZeH3RQg32hnf3slBk2f1vb5z3SNNlMZxnEZSwApmGrAuUr4WmFpX4ROyxiL7MPA14F+BLuBDZvZnSQcAvwN+U4RwjcCsf1iXlFQq8RAyaXFlYv4AQ9hSPXLaPO5Z/zgAf167gufNObBfnazhY4IMkfAvKQbQnhqPN8UM92n5ZGIhZCzlDYs5D1hqcJzs1GrkzxP+JY9DQHSsGmXIY7jv6an9mFvU8F6jMT6tblq4mXi/EQeIFIeC/iFkirPyy5SaD2kApku6o+L1RWZ2UVWd2H+UwiaSVcEcBHyEsFqx3nZm9oCkjxEU0I+KEHC085Qpc+gotdFV7uGxrRt4fOsGZo+b3GyxHMcpEMOGcsB6tZktHOD+OsIqppqpxFc2NZPVBlMGui3MeBWwV8W9R4H96i2YExjT1s6hu83ue/3XdY82URrHcRpF2F3JfmXgHoIdpppD6H8EpS5kVTD3A/OT53cA/yJpjqQZwHuBZfUXzenlsClz+p7ft2FlEyVxHKdhlC3fNThXAMdI6ovQkkRsOS65V3eybpFdBvTGjf8o4fTniuR1D/CaOsvlVHDg5JmIsDe5bNNatnTvYLzniXGckUvENlwH/ht4B/BzSR8Ko/BJ4GHgm/UeDLLng/laxfM/SToMOIXgVXa1mRWyvBrWRH5RpOeZ6V/UmyMGYGLbGPYcP5WHtqzDgAc2rOLwaXPrIyf58sGk1Y2Fteipg+E9RtR5IGWsok7iZ6WoHC9FyApxo3fsxH2om+10PtTuMFKcXMW8j3WhzqKZ2WZJzyWEivkuwbh/DSFUTCFnIAZVMJI6gXOBa8zs7kTQFcC3ihDIiXPQ5Jk8tCXY4e7bsLKuCsZxnNajiCjqScyxl9e94xQGtcGY2Q7gP4h7HzgN4uCK8zD3rl/pIfwdZyST08Dfql8HWY3897Jr6H6nwew1YSpjS2HBub5rG+t2bG2yRI7jFEr9jfwNJ6uC+Qjw4cT24jSBNpXYa8LOw7bLNg/rANaO4wxA7zmYPFcrktWL7IPAROBOScuAx9jVNG1mdnydZXOqmD9hGg9sXAXA8k1rOXLavCZL5DiOk05WBdNDQQdxWpHUHwNRb6+UyrF8MDWy94SdZrDKFUzWHDEQz+eSJ+9KreTzWIvPoSh5s+ZoyeMZVpS3Vy7vtIw5YqA4r6qoF1iLenulhdGpzhNTZD4YaF27Sh6yuikvKlgOJwN7V2yRrdiynu5yD+2R/7iO4wxzjJa1q+ShvinTnEKZ2D6G6WMmANBjZR7burHJEjmOUxSjyYsMSXMl/T9Jd0haKukpSfm/SDq6OBGdSuZUBLp8YpsrGMcZsYwWLzJJhwJ3Aa8nBLfcC+iNVbI38O5CpHP6MXPspL7nK13BOM7IJAkVM1q8yL5AOAtzMrAN2FFx7ybgs3WWq4FEtH8BBvoBidgvLSVvyqyxE/ueP9GEDJd5nAdiuV8a/f+g2aFi8pBquI/8Dmy2rJA9xwtAT09z5U3NfxMJz5Mnn0xRGK277ZWHrArmWcA/mdkmqZ/7yUpgdqSNUwC+gnGc0UDrbnvlIauCGejnx3TAj5U3iJkVK5hV2zZRNivcXdJxnMYzElYwWY38twFvTLn3KuCP9RHHGYzx7Z2Mb+sAoNvKbO7eMUgLx3GGHaPMBvNJ4GpJVwLfJ2wRniTp3cBLgecUJJ8TYXLHWLb0dAGwsWsbEzs6miyR4zh1p/lmtprJtIIxsxuA04F9gIsJeQT+A3g2cLqZ3VqUgE5/JnWM7Xu+oWtbEyVxHKcoRtMKBjP7FfArSQuAmcAaM7u/MMlGIGmfgbw2lMkdY/qeb+jaXoNEAxPzXqo1cVRRpIWgsYKSntVKI8PK5CGfZ1hzPwvDzTMsFyPkJH9mBdOLmS0GFhcgi5ORyhXMRl/BOM6IpEUXJbnIc5J/f0mXSnpA0ubk8ZJkRZO1j1dIulzScklbJd0v6QJJkyrqzJdkKdeUinofG6DeiP7Wndje2ffcjfyOM/II52BGyRaZpEXArwnuyL8inH2ZBbwIeLWkUxI7zWC8D3gI+HdgBXAE8DHgBEnPNNtlj+AC4Iqq9pUHP74F/Lbq/oSkrLrdiKKjtPPP1mXDbOnvOE42mr9jWjN5TvLfCZxsZn3Hx5OVx5XJ/YUZ+nmRma2qeH2DpLXApcAi4NqKew+a2S1pHZnZCoKS6kPS6wlzujSDLMOWjoq9567htrfsOM7gtPCqJA9ZFcwhwKsrlQuAmW2U9FngB1k6qVIuvdyePM7NKMtAnElYXf2uDn21LO0VwRS6U/Jn5EEpXgat+gGPyduiokYdJQDKI+HnqeMMQlYbzAp2BresphN4pAYZejNh3ltVfoGkbknrJV0xWLpmSfOAE4DLzKy7Bnlans6KHDC+ReY4IxMrW66rCCQdIOnLkv4maZOkx5Lv46dlaZ9VwXwW+LikXVYZyeuPAp/JJ/Yu7T8BXG1mdyTF24FvAucQFMb7gMOAmyQdPEB3ryfMZ0RvjwG0+xaZ44x4ymXLdRXE8wnfw5cSbO7/DMwAbpX09MEaZ90iOx6YBCyRdAs7jfzHJM8XJY4AAGZmZw7WoaSJwM+BbirC0JjZY8DbKqr+QdJvgXuA84HXpXT5BuBOM/tbhrHPBs4GmDd3+OW1L1fsB6VF4HUcZ/hiRmGrkpz8EPiaVeyXS7oWWEZI0/KGgRrniabcAzxGyP+yd1L+WPL47Iq6g74rksYSPL32BY5PDPapmNnDkm4EnpHS31HAQcC/DDZ20t9FwEUAhz/18Jb4K+ahu8LZrsMVjOOMSMotYFg0s9WRsvWSHiCD3TyTgjGzfYYgWxRJHcDlwFHASWZ2V9ampCuvMwkroe/XLmHrU7kt1l6qzp7gOM7wpzi7Sq1ImgY8BfifwermPslfCwrZqi4DTgROHcgNuardXsBxwE8j9zqBM4Bfp3ipjTgqPcc6UsJl5KFVvcXSGE7yureYMxTMKNKuUiv/RfjB/6XBKuZSMJL2BPYExlbfM7Nr+7fox9eAVwKfBjZLOqbi3gozWyHpCwRj/c3AKuBA4DzCsaOYM8FpwDRGgXG/l0rPsfZ++d8cxxkJDGEFM13SHRWvL0rMAX1IOgm4KkNfN5jZoupCSecBrwHenIQNG5CsJ/n3Jaw8juotSh6NnVtXWb7pXpA8np9clXyccKr/HuBc4CyCY8FqwgHMj6cE1zwTWAv8MsP4I4LtPTu9sMf4FpnjjDyGtoJZbWaDHXi/CRjIG7eXLdUFkt5G+JH/ITO7OItAWVcw3wL2IhjR7wOGFADLzOZnqHMxISVA1j5fMhRZhjOVAS4ndvRbTDqOMwIowgZjZlsI3+G5SKKkfB34gpl9Omu7rArmGcBZZnZ5XsGc+lMZor8ydL/jOCMDw1rCiwxA0ksJBv1vmdn78rTNqmBWMMRVi7OTvHlf0qhcwUwucAUznAzUrZr3JY1WyP0SoycSeihW1gqkydWq8uaiRc7BSHoOIRTY34BLquzm283szoHaZ1UwnwE+KOlaM9s8NFGderGhe+cKZqKvYBxnRNIKCgZ4LjCGEPn+j1X3lgPzB2qc9RzMdyUdBCxLTvKv619l8NP7Tu2UzXZJkzyp3RWM44w0jNZwUzazjxGcr4ZEVi+yswiuwj3AkfTfLmv+OzFKWN+1te+g5fi2Dsa1deBvv+OMMFpki6xWsm6RfZxwyPHNZvZkceI4g/HE1p0ZE2aNnYSkYWd/cBxnMAoNYNkwsiqY3YGvu3JpPiu37UzqOXPcpAFqOo4znBlOESvSyKpgbiQczrmmQFmahKBUJ/euoRKJ9pKWBOyJCgUza+zEoiRKxVK8nyzicdYK/0Fi3lqt6sGVJleryuseZ8XR4qFiMpNVwbwb+JGkdYSc99VGfiztm8epK09s27lFNnOsr2Acx2ldsiqY3myT30m5bzn6coaImbFiy5N9r2e7gnGcEctoMvJ/AndVajqrtm9iS08XABPaO9l9zIQmS+Q4TiGMpi2yxBfaaTLLNu/cmdx7wtRUO43jOMMba+F8MHnIva2VpDreHXjUzLrqL1LzSfvejpY30EFg+ea1fc/3njCt73nMTTmPMb6RpLlUR+fQYFmzOgS0gjE+jwyt4CTQU5Ekb8Ay61/WaFIdAnp2LS/ah2UkKJjM2aoknSbpz8B64EHgsKT8W5JeU5B8TgWVCmZ+hYJxHGeEkWyR5blakUwKRtLpwM8JuVk+yM58MABLCTlZnALZ2tPFo1s3AOHN33PClKbK4zhOsZhZrqsVybqC+SjwP2b2fPqnybybkJ/ZKZAHNjxBOdlGmjd+CuPbO5sskeM4RdEbi2y4r2Cy2mAOBj6QPK+eyTqCTcYpkPs2PNH3/KDdZjVREsdxCsdGl5F/AzA95d58YFVdpBlJRIz/aV5f0fKKtaWZcd/GlX2vD548s2bxKslneI/XLTfQiz0aNSBNriYbvXMZ41OcGoqQN82YnsfwHq9bzHtbu1zNdx7IS6uuSvKQdYvsKuA8SVMqykzSGOAdwG/qLZizk5XbNrJux1YAxra17+JB5jjOyMOSaMp5rlYkdQUj6UHgpWb2V+B84DbgfuDXhG2yfwOeCuwGnF64pKOYe9fvXL0cMGkmbaXMzn+O4wxTyi0aJy0PA31TzSdkMsPMlhHywPwSeB4hL8xzgFuAo83s0UKlHOX8Zd0jfc8PcfuL44x8LJ+Bv1W30zIftDSzFcCbC5TFibBm+2aWJyf4SxKHTd2jyRI5jlM0xsg4aDmYghn+Mxzm3Ll25+rloMkzmeDuyY4zKhgJW2SDKZiPS1qdoR8zs2F52FKKeHGlhH+JeXvlCiszBO5cu6Lv+RHT5kU9qGIeXKmeYZEDWWneS7USzxGT4ikVDXeT3butHsTeh1yhYjK2rwe1ytDo/ChZc8fkkasessY8zlLf237jFfj728B6hv/v+8EUzOHA9gz9DP93ogVZuW0jj2xdD0C7Shw2ZU6TJXIcpxEYNipWMKeb2W0NkcTpx22rl/c9P3i3WYxt62h6sErHcUYvkv4J+D7wiJnNG6y+JwlrUbrLPdxaoWCOnr53E6VxHKehtGA+mOQc5BeBx7O2cQXTovztycfY1L0DgCkd4zjY3ZMdZ1TRgl5knwP+CjwGnJSlgSsY+hvkUw30teZ+STl1FHMeuGnV0r7nx8yYT5tC46y5X+qxlZan35hBPk+omVrlrUd+lFrzweShCLnSymsNk5LP8F5b+Ja09nmo2XmgJ63urp/dIgMYm7WWDUbSccDrCIfrP5S1XaqCMTM/Lt4kVm7byOKNwXmvhDjGt8ccZ9TRKisYSR3ARcDnzWxxnky6voJpQW5c9WDf80OnzGZK57gmSuM4TjNooRXMBwlRXS7I29BXKS3Gpu7t3FJh3D9uxj5NlMZxnGZgQ8toOV3SHRXX2dX9SjpJkmW4rk/qLyDEonyHmW3LOw9fwbQYf3jiQbqSvem543bjwDqH5nccZzgwpAjJq81s4SB1biLk9xqMLcnjV4BrgVsqoul3AkpebzezrWmduIJpMdpLJcaU2tle7ubEOQek5pBxHGcEY8VskZnZFuC+HE0OAfYmJJasZh3wZeBf0hq7goF+3mHpicEGbztQH6nKomKj8nl7HMhxM/fh9jUP8dSps/t5V+UJCxMjn2dXLKxMdi+wmBdansRkqR5rMbkKSuyVx2OsqLAyRSQnS/X2yuHZFUsuVlSol3xyZa8bfW9T5GqkTaQ3ZXILcAYwtqrs34CnA68EVvRrUYErmCZiFr4q2wiKp2xGSWJ8eyfPmbVvc4VzHKeJGNYCRn4zu6W6TNJZhK2x6wdr7wqmSfz60b+zavtm1ndtZcHkGRw5bR5zxk1utliO47QCLXiSfyi4gmkC337wFh7Zsp4DJs1gbFsHd6x5iBtWLubVex/B4dPm9h2qdBxn9GIpBz6bjZmdlbWuK5gG85d1j/Do1g28cd+j2XP8FACWblnL759YwveX/YlN3Ts4bsZ82mqNGuA4zrDFfAUzcuhne8/x5Z4eViZWV2zu2UGPlZncMbbP6L/vpN3ZfcwExrd38vMVdzG2rZ1nTJ/XZ5OpJGv4lrRcKjHjcFqOljz9ZnUeqDXUDDQ2zEpRTgJF1Y0buIsxxkfDv9RoeK+PXDnqRlYJ1SFh0vooKi9Rb++uYJzcdFuZ7nKZtkRx9FiZNkrs1jmWF+xxMNt7uvnx8r8wd/wk9hi/G2bmrsqOM8owoyWM/LXS0M1+Sa+QdLmk5ZK2Srpf0gWSJlXVO1TSTyQ9KmmzpHskvVdSqkKUdF5yAvXG4mcydI6YMhfDuOKRewBoU6nP3XNSxxhOnL0/M8ZO4pcr7mVHuceVi+OMUoZwkr/laLQ1+X1AD/DvwCnAN4BzgaukYNmWtAdwPbAv4QDPi4CfAZ8HPh3rVNK+hHAGTxQpfD2Y0D6G42bsyz3rH+Oqx+8HdlUye4zfjadMmc3j2zbSXYfIso7jOM2i0VtkLzKzVRWvb5C0FrgUWEQISXAaMB04zsweSOpdK2k/4A2EwGvVfAO4DDiQFt/2ay+VOHr3vXl0y3puWr2UdpU4Yc7+tKnUZ3OZN34Kf3jiQbb2dDG+vbPZIjuO03Ba4xxMrTR0BVOlXHq5PXmcmzz2fqNuqKr3JBF5Jb0GOBI4rw4iFo6ZMX3MBE6dewgzxkzkupX/4Bcr7gagJFE246HN65jSOY7OUkvrSsdxCmKIwS5bjlb4Bjs+ebw3efw/4KPAVyW9H1gDnAi8Hvh4ZUNJUwkpPD9gZmuHZK9Q/xAuad2oLVv4l17DfD8DfQlEKN9jwm68bK+nct3KxVy/cjEPbFjFzLETaZP409pHeNX8pzGho7NfaJV4wrH+ZT3Et9fyhIqJh6VJ8wLLFhYm3butv7x18SIrICRLd7m7pvapdeswh5i31I7urv71UpKQxRODNc6za0dPf1nT5co+h9T3NiJDeqiYqs9jwd/prXoOJg9NVTCS5gKfAK42szsAzGylpGOBnwO9iVEM+JiZfa6qi88DDwCXNEbigTEzNnZvZ0vPDua074YBgn7uxr1KZ/a4ybx43qE8fdpcfv/Eg6zdsYXJHWM5c9+FHDLVUyQ7zqjFz8HUhqSJBCXSDbyxonwG8BNgM/AKwgrmucCHJG03s88m9Z5NsMkcaWk/i9PHPhs4G2DPeXvWPhmgu1zm4uW3snr7Ztbs2Mwe43bjqVPmcPyMBYxpa6fHyrTT1le/d3UzoX0M+0+ewf6TZ/T1014q0VOHlMeO4wxXzFcwQ0XSWOAKgqfY8WZWGZHzA8B8YG8z6w0Rfb2kNuCTkr5tZquBbwLfBlZU5CloB9qS11vNbHtsfDO7iJAClCMOP6Lmnwk9VuZrD95Iu0o8b+YBTO4Yyw2rl3DT6mXct2EVZ+93DGPbOsKZlwpjfiILaOcqp83dkh1n1GMG1jP8VzAND3qV5He+HDgKeKGZ3VVV5TBgcYVy6eU2oANYkLw+GHgbISdB73UccEzy/NxCJhBh2ea1rO/axgtnH8Izpu3FgZNm8sZ9j2bRzAWs2b6ZLz3we7b3dPdTLr999D5uWb1slzI/9+I4DgblnnKuqxVp6AomOetyGcFof2osFDTwOPBMSVOrlMzRyeMjyeMJkbZfAtqAdwKLM0rVPzRMnhwvbWJDOdhdZo6biNpEd7mHMR3tPHvWfkzoGMMvVtzNN5fcxNsPfDbtpaBkNnRt43eP3svUzvEcMmUWkzrG7NJvrvAteQz3sZAdaQ4BOcK3ZA0LkyfHS5oxPVbe0FAxKe27rb9cuRwCcsiVz0jfvyxm+A/tY8b42sbKUzdtrJjxP5dDQcoXcCwsTNZ+i11fDCmjZcvR6C2yrxGS1Hwa2CzpmIp7K5KtsguB1wJXSvo8wQaziHBI86dm9jBALBeBpCeB9ix5CurJnLGT6S6XuWv9Yxw7fT7tpTbKZrSXShw5bR5be3ZwzeMP8JtH7+W0uYdQkpjSOY73HHwC3dbTT7k4jjO6sWQFM9xp9BbZC5LH84Gbq663QF+Cm2cDqwjpOK8AXkbwNnttg+XNxO5jJnDApBncsmY5yzeHRVfvmZb2Uoljps9n7wnT+Pv6x9lR8QttzwlT2Gfi7s0S23GcliUY+fNcrUijD1rONzOlXB+rqHeLmb3QzOaY2QQzO9TMPmVmWwfpf5GZPavwiVTRWWrj5Xs+jZXbNvL7VUtYuW0jEJRMd7nMmLZ2Tt7jIB7fuoEVW55stHiO4ww3EiN/nqsV8cxWdWL6mAm8c/9nc+e6FVy98gEe2bIeCKFhADZ372BSx1gmeOgXx3EyUC6Xc12tSCuc5G8qov/J/bw5XnqZN2EK7zloEV+5//es79rKs2fux1On7sGTO7by9ycfZ0ypjTFtbf1OrceMnbGT7RA3vBeV4yUmQ8xJAOKOAo3M8RIzsENOh4CMJ+br4VAQkyttDlEjf40n5tPax4zpaQ4BNRvea4wEkCZXV3f/97GrO/65jZ/kT/mMVtfNd/wuFyPFTXnUK5h6s+f4qbznoBP44UN/5n+W3MrUzvF0lNrY1L2dNy04it06xzZbRMdxWh3zg5ZOBDNj1thJnLP/M1m6aS1LN61h2pjxLJg0neljxzdbPMdxhgllX8E41fx+1YPcv2Elr913IYdOmc2hU2b33Uvb9nIcx9kFz2jpVPPgpjX8fMVd/H3DSr7w9+tYvX1zs0VyHMepGUlzJV0s6XFJ2yUtlXTBYO18BVMnNnRt45Klt/WFqJ/Q3sluHW5vcRwnP4a1zEFLSfOBPwJLgXcBKwnxIhektwq4glH/PC+xvC8ApVJkwVcKwS4vXXo7G7q2AUG5nLXfM2gvaRevrVh+FMie4wXi3loxz6607bgicryEujHvtuyhRIrw9korz+NxVqu3Vz3mUKtnVzwfTDE5XmqVK08+mHrkeIl5l8W80AC6unatW6ATWd85mBbhQkKIrhPMrPcPdEOWhq5g6sAvVtzDkk2rgeD2/Pp9FjJtjBv0HccZOq3gRZakqj8ZeEOFcsmM22Bq5JZVy7h+5c64mi/Y4xAO2s2ThTmOM3TM8qVLLjA52XHJ41ZJVyX2l3WSviNp0DhXrmBqYPHG1fzoob/0vX7KlDmcNOeA5gnkOM6IoUVike2RPF5MyB78AuCDwKnA75II+an4FtkQWbVtExc/eAvlZCN2j3G78fp9Fu6SGtlxHGdIDM0GM13SHRWvL0qSK/Yh6STgqgx93WBmi9i5CLnezN6ePL9W0nrgh4Tts9+kdeIKhv55XlKTfiVv9ZbuHfz3gzezJTFITuoYw1v2P5rOtlKfwTyW8jgtzEquHC0Zc7+kGeOjDgGp4V9yzCEWKiYmVwON8Wl95MrRksMYn8chIFY3zcAdN9LnybuSwxifw/Beu1zZ88F09UTCv3SlGOOjhvs0uWK5Y1IcbPrlgynQCG9D8iJbbWYLB6lzEyFZ42BsSR7XJI/VSunK5PEIXMHUj55ymUuW3MYT2zYB0KESb97vaKZ2jmuyZI7jjBSMYg5amtkW4L4cTe6pECnGgEK6DSYHZsYPl9/JAxtX9ZX90z5HsvfEqU2UynGcEYeFUDF5roK4hZBl+JSq8t7Xtw/U2FcwOfjlI3/n9jUP9b0+ZY+DOGLa3CZK5DjOyKQ1gl2aWbekfwMukXQh8BPCActPA9cD1w7U3hVMRm54YgnXPP5A3+tjpu/NyXMOoujM3I7jjEJa6KClmV0qqUzwHnsjsBb4HnCepZ0IT3AFk4E/r1vBzx65q+/1obvN5pV7H46kYg19juOMSoz06APNwMy+C3w3bztXMIOEirl/wxN8f/mf+l7PnzCVN+z7dEoKIVZinlaQz1srq2dYWh9Zw8ekt8+RnCy1bn95i/L2yhVWJk8SsIwy1GMOeTyoYnXTQ7JEvNNinmG5kojVQ65sMsS8xSDuMZbLMywtVExXjlAx3Q0MFYOlfk6HE65gBmD5prVcvPhWepJP0qyxE3nL/sfQ2eZvm+M4xWHWWiuYoeLflCms2PIkF/7jJrYnv0CndIzjnAOeyYT2ziZL5jjOyMcop6wchxOuYCI8tnUD33jgj2xNlvUT2js554Bn+lkXx3EagpEe8Xo44edgqnhi2ya+fv+NbO7eAcC4tg7OPeA45oyb3GTJHMdxhhe+gmGnUX/19s18ffGNbOzeDsCYUjvn7H8sc8dPTjXopxvuYwby2nK85KmbJ/xLmuE+nmcmR5iUGkOn5AnpkqvfgvLBxOqmhU6pNSRLzJifWrfGkC755MpeN0/4l55ItOBUw30sVEzEmB/qZnce6G/kLzJUjNtgRhTrtm/h6/ffyPokaVhnqY2z9z/GT+k7jtMEzBXMSGHN9s187f4bWbsjxHfrUIm3LDiafScNmu7AcRyn7hjp7vTDCVcwwN/WPdqnXNok3rjgKPafPKPJUjmOM2ox9yIbMSyatYBN3du5YeUS3rTgaA7ebWazRXIcZxTTaif5h4orGEL+l9PmHspR0/dm1thJqUZvx3GcxuAn+UcMRhkEM8dOyB3+Jd0DKxKuIiWJWJ6EY1k9xvKEf0lPIhbx1soRvqWo8C+5+q1RhlrDv9QaOiXUzRb+Ja08T/iXfHJlr5vVYyxP+JdUb68awr8M1G+3e5HlxhWM4zhOi+FbZI7jOE5BWOpZpeGEKxjHcZwWw4NdOo7jOAXhBy1HDNUG8Xy5VLLneKk1/Eu6DHkM9zFHhbixsohcKHnapxr5CwpBk3UOaeFQag+zUlv4l7Tx8oSKKSL8C2Q36Ocy8hcQ/gX6G/N3jrdrvwXb+N2LzHEcxykAGxkrmIZGU5Z0sqRrJT0uabukFZJ+JOmQqnpTJX1L0mpJmyVdLemwSH/7SPqxpCeTetdJWti4GTmO49SfsILpyXW1Io0O1z8N+BPwDuD5wHnAocAtkvYGkCTgCuAU4J3Ay4EO4DpJ83o7krQ7cCPwFOAc4Izk1nWSDm7IbBzHcQohrGDyXK1IQ7fIzOwHwA8qyyTdBtwHvAL4AvBi4FnAc83suqTOzcBS4APAu5Km5wKzgOPNbHFS71rgQeDjwKuKno/jOE4hjBAvslZIOLYmeey1GL4YeLRXuQCY2XrgF8BLKtodA/yjV7kk9TYDfwBOk+T2JcdxhiW9GS3zXEUhaXdJX5b0oKStkpZK+qqkQSMCN+VLWFIb0AbsDfwH8Djww+T2ocDdkWb3AG+QNNHMNgE9wI5Ive3AOGA/4P7BZDGs3/5lXby9YnVTvLVqTRgWTwyW0j5St7DEXg309irK46yo0CmxL4Raw7/kqZsuV/ZwN3kShmX1GEsN6RLxGMuVcCxj+JcwVrY5pHlfjiQqTBYHAB8B7gUOAT4JPF3SM22AmDnN+pV/K/D05PliwnbYE8nracCySJu1yeNUYBNBeTxP0u5mtgZAUgk4qqIfx3GcYUjLeJHtDzwTOMfMLkrKrpdUBr5BUDypP+SbtUX2esIW12uADcBVkuYn9wTRnwaqen0hQf7vSNpP0hzgK8A+yf3Uv46ksyXdIemONatXD30WjuM4BWBJPpg8V0F0Jo8bqsqfTB4H1CFNUTBmdq+Z3ZoY/U8EJgL/ltxeS3z10Zu7eF3Sx4PAawkrocXAo8CxwBeTeo8NMP5FZrbQzBbuPn16rdNxHMepO2Ur57oK4h7g98CHJS2UNFHSUYTtst+Y2b0DNW66kd/MniQoiAVJ0T0EO0w1hwAPJfaX3raXA3OTewvM7OkEZfWwmT1UpNyO4zhF0RtNudluyol95YWEbbDbgY0EE8eDhCMkA6JCcxpkQNIsYAlwmZmdI+l04KfAIjO7IakzmeCm/H0ze+cAfe1BcBD4vJldkHH8VcByYDowUvfLRurcfF7Dj5Ewt9457G1mheRWl/TbZJw8jAW2Vby+qMJu0tvvScBVGfq6wcwWJW2+DywiHP+4Fzg4ef4n4EVm6cunhioYST8F/gz8jbCndwDwr8Bs4CgzeyAx1N8I7Am8n7Aldh7wVOBpZvZw0lcH8DnghqSvQ5N6S4ATzSzmYTaQbHeY2YiMAjBS5+bzGn6MhLkN5zlIGg/slaHqFjN7SNKpwC+Bk8zsmop+ngdcCZxuZj9P66TRXmS3EA5AvpdgPHoYuB64wMyWAZhZWdJpwH8CXydo5ZuBE3qVS4IRPBxeA0wBVgAXA5/Jq1wcx3FGA2a2hXCwPSu9Ibpuryq/LXk8GGgNBWNmnwU+m6HeWuBNyZVWpxs4rX7SOY7jOFU8njweBVxdUX508vjIQI2bbuRvIS4avMqwZaTOzec1/BgJcxsJc8jKTwgeut+RdK6kEySdC3yHsAP104EaN93I7ziO47QukvYEPgY8F5hDOAJyNfAxMxtwBeMKxnEcxymEEb9FJmnPJGfMekkbJP1EUhYvCiSNlfR5SY8lQd5ulvScomXOQi3zqurnPEkm6cYi5MxLjX+vvSRdKukhSVskPSDpU5ImFC13BtnmSfqv5DO0JXnP52dot1DSRZLuS9o9JOkySfsM1rZRDHVuFe0PlvR/Cvmftkq6X9K7CxQ5JsOQ5yDPS5XKiFYwiUvetcBBwJmEEDX7E3LGZPnS+TbwVsKp1dMIS8PfSTq8EIEzUod59fazL3A+8MRgdRtBLfNK7l8NPAf4MHAq8C2Cx+LFBYqdlQUED8p1hIjfWTmD4IL/FeAFhIgXRwJ3JFsXrcBQ50byRXwrMAZ4C+FQ3xcIwXAbyZDmIM9LNTBmNmIv4N2EqMsLKsr2AbqB9wzS9mkEV+g3VpS1E060XjFc51XVz++AbxJcxW8c5n+v5yd/r+dXlf9H0n58k+dWqnj+lkTW+RnazYiU7U2ItfeJZv/NapxbiRC546fDeA4fSj5flZ/ZCcBK4EfNnlezrxG9giHklrnFds0ZsxT4I7vmlklr2wX8b0XbbkJagZMljam/uJmpZV4ASHoN4ZfweYVIODRqmddAQflK9A+W2lBsgNPOg7RbFSlbDqwihElqOkOdG+F0+CHA/6ufNEOjhjl4XqoBGOkKZqDcModkaLvUwsGk6rad7Iyd1gxqmReSphKCgn7AwpmjVqGWeV0N/AP4rKRDkqB8zyWsii5M/tOPCJKtl5mEsB3DmWclj2Ml3SKpS9ITkr4iaVxTJctOlrxUo5aRrmCmkURfrmItO6MzD6Vt7/1mUcu8AD4PPABcUkeZ6sGQ52Vm2whfWL3bLhuBawhhLt5RXzGbR/KL+ELCCubbTRanVvZIHv+XEHbkeYTwT28Bvt8soXJyP7B/YosBPC9VJaNh+ZYlt0yMrHlpmsWQZJP0bOANwJGWbBi3GEOd11jCF9VMgnPAQ4T/5B8h7JGfW0cZm8lXCQmgTjWzmDIeTvT+wP2emX0keX69Qsbb/5B0iJn9vUmyZeVC4F2Eg4jvArYQHGcGzUs1GhjpK5h1pOeWGew/52B5aZq5tVTLvL5J+OW7QtIUSVMIPzTaktfNtC3VMq83E/b0X2hm3zOz35vZfxK8yN4m6Wl1lbQJSLoAOBt4k5ld2Wx56sCa5LE6um/v3A5vnChDw2rISzUaGOkKZqDcMoP9MroH2Cdxna1uu4PwYWoWtczrYOBthC/s3us4grFyHc39pV/LvA4D1pnZkqryyqB8wxZJ5xNclN9tZt9ttjx14p7ksXrV2rtiHRa//s3zUqUy0hXMFcAxyXkPAJLDU8cl9wZr2wG8sqJtO/Bq4Eoz2153abNTy7xOiFx/JRjXTwB+XIC8WallXo8DUyVVO19kCsrXyiRbL58Czjez/2q2PHXkNwRj+ClV5Scnj3c0VpyhY2Y9FjL1LlHIS/VqQs760U2z/aSLvAj+6IuBuwhuri8mfJk+CEysqLc3YZ/+I1Xtf0j4Vf8WQmrnHxMS+hw5nOcV6e96WuMczJDnBcwnuCg/QDikeQIhn9AGwhdVqdHziczvFcn1DcKv9nOT18cPMK8zCL/kf0NYZVZehzR7TrXMLSn/aFL+GeAkwiptK3DJcJgD4UfoF4HTCbG63knYJvsD0Nnsv0uzr6YL0IAPzV7A5ckXzUbgZ1QdoEq+nIwQvK2yfBzBR//xRLHcSsi0OaznFemrJRRMHf5ehwA/IkR53Zoom/8EpjZ7Xol8lnJdnzYvgqffgO1a4RrK3JJyAe8h/LDYQcgu+wmgYzjMgWC//CXhYOV2QsLDT9Hkg72tcnmwS8dxHKcQRroNxnEcx2kSrmAcx3GcQnAF4ziO4xSCKxjHcRynEFzBOI7jOIXgCsZxHMcpBFcwDpLOSlLE9l47JC2R9JkkiGTe/hYl/SwaQtuPJWH2q8svkbQsb3+1kMhisZwekhYk985qpEzJ2Jnei4q/6/yKsmWSLql4vSiZp38XOHXHP1ROJa8kBOo7lZDt8jxCaP9G8lHCiehqPgm8tMGytCq1vBcvTdr3sojwnvt3gVN3RkO4fic7f7GdmfmukrQ/8GZJ77ahZ/yrC9Y/iOWopZb3wszurKcsjjMQ/qvFGYg/E8LlTO8tkDRe0mclLU220pZKOn+wLRZJz5f0a0mPSdoi6W5J701yf/TW6Q0rcX7Fdt3Hknt920KSxkhaK+kLkXFenbQ7vKLseEnXSNooabOk30l6ytDfltQ5RreuJF0v6fqK171biKdL+mYyl3WSviipTdIzJN2YyHqPpJMHG0fSvpJ+lby3qyR9GeiXeqFyiyx5bz+a3OqqeM/HJH18MdK+d9vtoLzvjzP6cAXjDMR8YD1J3o7EFvE7QvDPLwMvAL4FfJjBt9L2JWSYfBNhC+5S4GPApyvqHJs8XpI8PzbpfxcsRLL+EfCaSgWV8DrgbjP7SyLzqcm4m5J7rwEmAX+QtOcgMvfSJqm98gKqxx0KXwI2EyLvfhX4l6TsO8DFwMsIeYd+Iml6tAdAUichp8oRwNuBswgJrz40yPjfYmdWzGeRvOfJ+/s/wJkRG9w5wA1mdl+G+TmjnWYHQ/Or+RfhC8mAAwnbplMJiqAbeEdFvdcn9Z5T1f58QqDCmcnrRUm9RSnjKRnnfEK06lLFPQM+FWlzCbCs4vVxSd2TK8pmAF3AByrKFgPXVPU1GVgNfGmQ9+VjpAdA7L3OSpOxovx6KgJTVrw/F1fV+3NS/qyKsqcmZWcO8F68NalzTEVZb+pooyJYKLCMikjFFXNsr5JlH0K++ddHZDmj2Z9Zv4bH5SsYp5L7CF/Qawm/bL9pZl+tuH8KIdrtTVW/5q8khC0/Jq1jSXOS7aDlBGXURYg6O4WQ5jgXZvZHQuTa11cUn0H4Yr0sGXN/YD/gsip5twA3A8/JONwxwDOqrno4HPym6vV9wGYzu7GqDGCg1daxhORWt/QWWLCZ/WiogpnZUsJq9ZyK4nOAVcBPhtqvM7pwI79TyUuBFYSVwHuAf5Z0q5l9J7k/k5AToyul/e6xwsQ+cwWwB+EX832EcPqnE1YxuV2hE74HvF/SRDPbRFA215pZb3KxXsX1bXZuBVWSNdvgn8ysu7JA0pNDkLeadVWvdwC79GtmOyTBwO/RHEK4+GpiZXn4OvCLxF61lLDFeKGZ7aixX2eU4ArGqeRuS7zIJF0L/A34vKTLzWwzwRazFHhVSvtlKeX7AQsJ2y3f6y2U9KIa5f0uwUj9Ukm3ElYWZ1bc7835fh5wdaR9vb8otwGdkfLdK2QpgseIp5qeVWO/vyb8Tc8hJH6bBFxUY5/OKMIVjBPFzLZLej/wc+CfCUb83wIvBzZZPiPv+OSxb+UjqQN4baTuDoLnWhYZl0i6mbByOYBgMK/cvrmf8AV5qJn9Rw55h8pyYJak6Wa2GkDSfgTb1k0Fjnsz8EZJx/RukyWrxrQfApX0pv4eR0jw1oeZlSV9k5Bl8tnA1ebu4k4O3AbjpGJmVwC3A++TNI5g27gJuEbSeySdKOkFkt4h6UpJ41O6upfw5ftpSa+Q9BKC11OMvwOnSnqepIUK+c0H4juEdNZvAX6abJX1ym8Er6ozJP2vpJcnLsuvkvQlSe/J+FZk5f8IRvDLJJ0s6bUEBb26zuNUcykhrfRPEjfiFxIygU7O0PbvyeN7JR0taWHV/W8TtueeBlxYJ3mdUYIrGGcwPkSwZbzNzLqAk4H/Bs4mbKFcRtiWuomULadkz/50Qurp7wBfA34PxFYV7yCsRH5BUG5nDyLf/xK83WYTtsyqx/41wZg/geCW+zvgc0n9mwfpOxfJ9uIrgLmEL/gPEGxZD9RznMi4O4DnAX8h2E0uJWxlfipD818mbf6Z8H7cXtX3KuAGwjbcFXUT2hkVeMpkx3FSkTSV4AzxJTP7cLPlcYYXboNxHKcfkmYQbEfvJux0fL25EjnDEd8icxwnxqnAH4CjCIc8H2uyPM4wxLfIHMdxnELwFYzjOI5TCK5gHMdxnEJwBeM4juMUgisYx3EcpxBcwTiO4ziF4ArGcRzHKYT/D/CuUjR+PTdvAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 576x360 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "from pylab import meshgrid,cm,imshow,contour,clabel,colorbar,axis,title,show\n",
    "from matplotlib import colors\n",
    "\n",
    "fig = plt.figure(figsize=(8, 5))\n",
    "\n",
    "# Define the bounds\n",
    "rh_min = np.min(flattened_input_train[:, loc['rh']])\n",
    "rh_max = np.max(flattened_input_train[:, loc['rh']])\n",
    "ta_min = np.min(flattened_input_train[:, loc['ta']])\n",
    "ta_max = np.max(flattened_input_train[:, loc['ta']])\n",
    "\n",
    "# Number of pixels in x- and y-direction\n",
    "N_x = 50\n",
    "N_y = 50\n",
    " \n",
    "x = np.linspace(rh_min, rh_max, N_x)\n",
    "y = np.linspace(ta_min, ta_max, N_y)\n",
    "X,Y = meshgrid(x, y)\n",
    "\n",
    "Z = -(-3.00757*X-31.6875*(0.00389006*Y-1)**2+0.983277)\n",
    "\n",
    "ax = plt.subplot(111, xticks=[0,10,20,30,40,49], yticks=[0,10,20,30,40,49])\n",
    "ax.set_xticklabels(['%.1f'%val for val in [x[0],x[10],x[20],x[30],x[40],x[49]]])\n",
    "ax.set_yticklabels(['%d'%val for val in [y[0],y[10],y[20],y[30],y[40],y[49]]])\n",
    "\n",
    "# plt.xticks(np.arange(0, N_x, N_x/5), ['%.1f'%val for val in np.linspace(min(x), max(x), 5)])\n",
    "# plt.yticks(np.arange(0, N_y, N_y/5), ['%d'%val for val in np.linspace(min(y), max(y), 5)])\n",
    "\n",
    "plt.title('a) $\\partial_{RH}\\,f$')\n",
    "\n",
    "im = ax.imshow(Z, cmap='PRGn', vmin=-8, vmax=8) \n",
    "cset = ax.contour(Z, [0], linewidths=3,cmap=cm.Set2)\n",
    "plt.clabel(cset,inline=True,fmt='%1.1f',fontsize=16)\n",
    "colorbar(im)\n",
    "\n",
    "plt.xlabel('Relative Humidity')\n",
    "plt.ylabel('Temperature [K]')\n",
    "plt.savefig('derivative_of_f_wrt_rh.pdf')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 121,
   "id": "5279512a-b7c2-4892-a543-792fb6bbbd20",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.983277\n",
      "0.3269340364480295\n"
     ]
    }
   ],
   "source": [
    "T_opt = 1/0.00389006\n",
    "cond_max = -31.6875*(0.00389006*T_opt-1)**2+0.983277\n",
    "rh_max_increase = cond_max/3.00757\n",
    "\n",
    "print(cond_max)\n",
    "print(rh_max_increase)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 130,
   "id": "157113d7-c53d-4152-9c37-e4e40629eb72",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The condition can only trigger for T between 211.861 and 302.278.\n",
      "The condition sets RH to be at least 0.3269 - 10.5359(0.00389T - 1)^2 (which decreases cloud cover)\n",
      "The condition is most active for T = 257.0654 and RH = 0. It then sets RH to 0.3269.\n"
     ]
    }
   ],
   "source": [
    "print('The condition can only trigger for T between %.3f and %.3f.'%((-np.sqrt(0.98/31.6875) + 1)/0.00389, (np.sqrt(0.98/31.6875) + 1)/0.00389))\n",
    "print('The condition sets RH to be at least 0.3269 - 10.5359(0.00389T - 1)^2 (which decreases cloud cover)')\n",
    "print('The condition is most active for T = 257.0654 and RH = 0. It then sets RH to 0.3269.')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "51702441-2dbb-4e54-ab31-796a587b2fb7",
   "metadata": {},
   "source": [
    "Optimize coefficients"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 38,
   "id": "26e734a3-580a-48e4-b2c4-a891a42b5685",
   "metadata": {},
   "outputs": [],
   "source": [
    "# See ~/workspace_icon-ml/symbolic_regression/finding_symmetries/pysr_results_dyamond_on_regimes/no_of_regimes_2/notes.txt\n",
    "def func(X, a,b,c,d,e,f,g,h,i,j):\n",
    "    x0 = X[:, 0] \n",
    "    x1 = X[:, 1] \n",
    "    x2 = X[:, 2] \n",
    "    x3 = X[:, 3]\n",
    "    x4 = X[:, 4]\n",
    "    \n",
    "    # Modified to always satisfy RH-constraint\n",
    "    x0 = np.maximum(x0, 1/(2*c*d)*(-c*x1**2-a))\n",
    "    \n",
    "    return a*x0 - b*x1 + c*x0*(d*x0 + x1**2) + x4**2*(e*x4 + f) + g - h/(x2 + i*x3 + j)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 39,
   "id": "5e1685e5-d5ee-4ea7-b4e5-af0566792377",
   "metadata": {},
   "outputs": [],
   "source": [
    "import scipy as sci\n",
    "from scipy.optimize import minimize"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 40,
   "id": "6f1fdcef-6f16-45f7-a22c-17e90674665e",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Warning: Desired error not necessarily achieved due to precision loss.\n",
      "         Current function value: 158.700476\n",
      "         Iterations: 30\n",
      "         Function evaluations: 484\n",
      "         Gradient evaluations: 44\n"
     ]
    }
   ],
   "source": [
    "def objective(P, X,Y):\n",
    "    '''\n",
    "        The objective function.\n",
    "    '''\n",
    "    a,b,c,d,e,f,g,h,i,j = P\n",
    "    train_preds = np.minimum(np.maximum(func(X, a,b,c,d,e,f,g,h,i,j), 0), 100) \n",
    "    train_mse = np.mean((train_preds - Y)**2, dtype=np.float64)\n",
    "\n",
    "    return train_mse\n",
    "\n",
    "T_subset = 10**6\n",
    "inds = np.random.randint(0, train_data_scaled.shape[0], T_subset)\n",
    "\n",
    "(a,b,c,d,e,f,g,h,i,j) = (36.69, 36.69, 16.48, 1, 0.21, 2.76, 50.29, 2.86, 1, 0.76)\n",
    "\n",
    "# Nelder-Mead gives me the same result\n",
    "res = minimize(objective, (a,b,c,d,e,f,g,h,i,j), args=(train_data_scaled[inds], flattened_output_train[inds]), \\\n",
    "               method='BFGS', options={'disp': True})"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 41,
   "id": "373a1caf-7c49-4937-8e50-11f244f3a26c",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[38.15388776 43.63539083 19.93732831  1.15446879  0.35551696  4.06188388\n",
      " 44.23799326  2.03600412  0.66480254  0.63943517]\n"
     ]
    }
   ],
   "source": [
    "print(res.x)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7878720f-5901-44d6-b696-efa29667c02b",
   "metadata": {},
   "source": [
    "New values"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 42,
   "id": "65f8d6d2-042e-4e35-938d-9c791591f698",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "159.23523622241254"
      ]
     },
     "execution_count": 42,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "objective(res.x, train_data_scaled, flattened_output_train)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 43,
   "id": "5fa3fb5b-c1bf-4d1b-bc40-1b72ceab13af",
   "metadata": {
    "tags": []
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "157.60467889204403"
      ]
     },
     "execution_count": 43,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "objective(res.x, valid_data_scaled, flattened_output_valid)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 48,
   "id": "0bb7874c-7ed5-4cb6-ae61-1a660393c97b",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "On the entire dataset\n",
      "Valid MSE: 103.86753\n"
     ]
    }
   ],
   "source": [
    "# From finding_symmetries/pysr_results_dyamond_on_regimes/save_optimized_eqns.ipynb\n",
    "mse_reg_0 = 0.0353\n",
    "n_0 = 32419018\n",
    "n_21 = 62640812\n",
    "N = n_0 + n_21\n",
    "\n",
    "valid_reg_mse = 157.60467889204403\n",
    "\n",
    "print('On the entire dataset')\n",
    "print('Valid MSE: %.5f'%((n_0*mse_reg_0 + n_21*valid_reg_mse)/N))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 49,
   "id": "f26f558e-c50a-41c4-82bf-9d8a4b0a161e",
   "metadata": {
    "tags": []
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "157.63718734237756"
      ]
     },
     "execution_count": 49,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "objective([38.6562122, 43.53500518, 19.78403208, 1.13637902, 0.35299939, 4.04888686, 44.21730274, 2.03128527, 0.66971589, 0.6409019 ], \\\n",
    "          valid_data_scaled, flattened_output_valid)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 50,
   "id": "101e7cad-761a-414a-b5a1-88daee663d9b",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "On the entire dataset\n",
      "Valid MSE: 103.88895\n"
     ]
    }
   ],
   "source": [
    "# From finding_symmetries/pysr_results_dyamond_on_regimes/save_optimized_eqns.ipynb\n",
    "mse_reg_0 = 0.0353\n",
    "n_0 = 32419018\n",
    "n_21 = 62640812\n",
    "N = n_0 + n_21\n",
    "\n",
    "valid_reg_mse = 157.63718734237756\n",
    "\n",
    "print('On the entire dataset')\n",
    "print('Valid MSE: %.5f'%((n_0*mse_reg_0 + n_21*valid_reg_mse)/N))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1993d671-c8c0-49e4-8064-7e275269ad7c",
   "metadata": {},
   "source": [
    "How often do we need to artificially increase the RH?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 45,
   "id": "b555d65d-8fac-4c62-8be3-191176a599ab",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.009560125771782007"
      ]
     },
     "execution_count": 45,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Optimal params without the RH-modification\n",
    "(a,b,c,d) = (38.6562122, 43.53500518, 19.78403208, 1.13637902)\n",
    "\n",
    "# In roughly 1% of all cases\n",
    "len(np.where(train_data_scaled[:, 0] < 1/(2*c*d)*(-c*train_data_scaled[:, 1]**2-a))[0])/len(train_data_scaled[:, 0])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 47,
   "id": "2b65f161-7fd5-4e71-a874-fc3eb19d4b20",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.010437487368024214"
      ]
     },
     "execution_count": 47,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Optimal params with the RH-modification\n",
    "(a,b,c,d) = (38.15388776, 43.63539083, 19.93732831, 1.15446879)\n",
    "\n",
    "# In roughly 1% of all cases\n",
    "len(np.where(train_data_scaled[:, 0] < 1/(2*c*d)*(-c*train_data_scaled[:, 1]**2-a))[0])/len(train_data_scaled[:, 0])"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (based on the module python3/2022.01)",
   "language": "python",
   "name": "python3_2022_01"
  },
  "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.9.9"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}