Workbook.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Tutorial on OQMD \n",
    "### (Nove 23, 2021)\n",
    "\n",
    "\n",
    "## Part 1: Querying OQMD and Retrieving data using OPTIMADE"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Import libraries\n",
    "\n",
    "import requests\n",
    "from copy import deepcopy\n",
    "import os\n",
    "import shutil\n",
    "import json"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Set the base URL for OPTIMADE REST queries\n",
    "\n",
    "rest_base = \"http://oqmd.org/optimade/structures?\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Created Query: \n",
      "\n",
      "http://oqmd.org/optimade/structures?filter=stability<=0 AND elements HAS \"O\" AND nelements=3 AND _oqmd_band_gap>0&response_fields=id,_oqmd_entry_id,lattice_vectors,cartesian_site_positions,species_at_sites,_oqmd_band_gap&page_offset=0&page_limit=200\n"
     ]
    }
   ],
   "source": [
    "# Create the query URL with filter, response_fields and paging requirements\n",
    "\n",
    "# The following query filters for data of ternary non-metallic oxides\n",
    "# Crystal structure parameters and Band gap values are returned in response_fields\n",
    "# Crystal structures will be used to generate representational vectors (input features) for ML\n",
    "# Bandgap values will be used as targets for ML\n",
    "\n",
    "\n",
    "filter_   = 'stability<=0 AND elements HAS \"O\" AND nelements=3 AND _oqmd_band_gap>0'\n",
    "response_ = 'id,_oqmd_entry_id,lattice_vectors,cartesian_site_positions,species_at_sites,_oqmd_band_gap'\n",
    "page_     = [\"page_offset=0\", \n",
    "             \"page_limit=200\"]\n",
    "\n",
    "filter_   = 'filter=' + filter_\n",
    "response_ = 'response_fields=' + response_\n",
    "oqmd_optimade_query = rest_base + \"&\".join([filter_, response_]+page_)\n",
    "print(\"Created Query: \\n\\n{}\".format(oqmd_optimade_query))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Success!\n"
     ]
    }
   ],
   "source": [
    "# Now do a test query on URL using requests.get()\n",
    "\n",
    "response = requests.get(oqmd_optimade_query)\n",
    "if response.status_code == 200:\n",
    "    print(\"Success!\")\n",
    "    #print(response.json())\n",
    "else:\n",
    "    print(\"Query failed. Status: {}\".format(response.status_code))\n",
    "    print(\"Error Message: {}\".format(response.text))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "Querying: http://oqmd.org/optimade/structures?filter=stability<=0 AND elements HAS \"O\" AND nelements=3 AND _oqmd_band_gap>0&response_fields=id,_oqmd_entry_id,lattice_vectors,cartesian_site_positions,species_at_sites,_oqmd_band_gap&page_offset=0&page_limit=200\n",
      "Success!\n",
      "\n",
      "Querying: http://oqmd.org/optimade/structures?filter=stability%3C%3D0+AND+elements+HAS+%22O%22+AND+nelements%3D3+AND+_oqmd_band_gap%3E0&page_limit=200&page_offset=200&response_fields=id%2C_oqmd_entry_id%2Clattice_vectors%2Ccartesian_site_positions%2Cspecies_at_sites%2C_oqmd_band_gap\n",
      "Success!\n",
      "\n",
      "Querying: http://oqmd.org/optimade/structures?filter=stability%3C%3D0+AND+elements+HAS+%22O%22+AND+nelements%3D3+AND+_oqmd_band_gap%3E0&page_limit=200&page_offset=400&response_fields=id%2C_oqmd_entry_id%2Clattice_vectors%2Ccartesian_site_positions%2Cspecies_at_sites%2C_oqmd_band_gap\n",
      "Success!\n",
      "\n",
      "Querying: http://oqmd.org/optimade/structures?filter=stability%3C%3D0+AND+elements+HAS+%22O%22+AND+nelements%3D3+AND+_oqmd_band_gap%3E0&page_limit=200&page_offset=600&response_fields=id%2C_oqmd_entry_id%2Clattice_vectors%2Ccartesian_site_positions%2Cspecies_at_sites%2C_oqmd_band_gap\n",
      "Success!\n",
      "\n",
      "Querying: http://oqmd.org/optimade/structures?filter=stability%3C%3D0+AND+elements+HAS+%22O%22+AND+nelements%3D3+AND+_oqmd_band_gap%3E0&page_limit=200&page_offset=800&response_fields=id%2C_oqmd_entry_id%2Clattice_vectors%2Ccartesian_site_positions%2Cspecies_at_sites%2C_oqmd_band_gap\n",
      "Success!\n"
     ]
    }
   ],
   "source": [
    "# We need more than 200 datapoints for machine learning - if more data is available\n",
    "\n",
    "# As the first step, here's the same script from the cell above, but kept inside a function \n",
    "\n",
    "def query_oqmd_optimade(query):\n",
    "    print(\"\\nQuerying: {}\".format(query))\n",
    "    response = requests.get(query)\n",
    "    if response.status_code == 200:\n",
    "        print(\"Success!\")\n",
    "        return response.json()\n",
    "    else:\n",
    "        print(\"Query failed. Status: {}\".format(response.status_code))\n",
    "        print(\"Error Message: {}\".format(response.text))\n",
    "        return \n",
    "\n",
    "    \n",
    "    \n",
    "\n",
    "# Next, we query for 1000 materials in total using 5 sequential API queries - each paginated to\n",
    "# retrieve 200 materials\n",
    "\n",
    "\n",
    "load_data_from_saved = True  \n",
    "# This is to avoid querying OQMD repeatedly for the same data, if .\n",
    "# Because the data I downloaded is already available as a JSON file in this Git repo.\n",
    "\n",
    "# But if you'd like to try out querying OQMD, set \"load_data_from_saved\" as \"False\"\n",
    "# Querying OQMD for this particular data would take about 5-10 minutes to complete\n",
    "\n",
    "\n",
    "dataset_filename = \"./raw_dataset_saved.json\"\n",
    "\n",
    "if load_data_from_saved:\n",
    "    with open(dataset_filename, 'r') as fin:\n",
    "        dataset = json.load(fin)\n",
    "else:\n",
    "    dataset = []\n",
    "    query = oqmd_optimade_query\n",
    "    for i in range(5):\n",
    "        jsondata = query_oqmd_optimade(query)\n",
    "        if jsondata is None:\n",
    "            break\n",
    "        else:\n",
    "            # Get the link to the next page and query it in next loop iteration\n",
    "            query = deepcopy(jsondata['links']['next'])\n",
    "            dataset.append(deepcopy(jsondata))\n",
    "    with open(dataset_filename, 'w') as fout:\n",
    "        json.dump(dataset, fout)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0:page_offset=0, page_limit=200\n",
      "1:page_limit=200, page_offset=200\n",
      "2:page_limit=200, page_offset=400\n",
      "3:page_limit=200, page_offset=600\n",
      "4:page_limit=200, page_offset=800\n"
     ]
    }
   ],
   "source": [
    "# Confirm the paginated response - checking reliability of server-side and client side scripts\n",
    "\n",
    "for i in range(len(dataset)):\n",
    "    query = dataset[i]['meta']['query']['representation']\n",
    "    page_params =  [param for param in query.split(\"&\") if param.startswith(\"page\")]\n",
    "    print(\"{}:{}\".format(i,\", \".join(page_params)))\n",
    "    #print(dataset[i]['meta']['query']['_oqmd_final_query'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "dict_keys(['id', 'type', 'attributes'])\n",
      "dict_keys(['lattice_vectors', 'species_at_sites', 'cartesian_site_positions', '_oqmd_entry_id', '_oqmd_band_gap'])\n",
      "dict_keys(['query', 'api_version', 'time_stamp', '_oqmd_data_in_response', 'data_returned', 'data_available', 'more_data_available', 'provider', 'warnings', 'response_message'])\n"
     ]
    }
   ],
   "source": [
    "# Inspect the response data keys and confirm all the necessary information is available\n",
    "\n",
    "print(dataset[0]['data'][0].keys())\n",
    "print(dataset[0]['data'][0]['attributes'].keys())\n",
    "print(dataset[0]['meta'].keys())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "structures\n"
     ]
    }
   ],
   "source": [
    "# As you can see, the type of the returned data is OPTIMADE's \"structures\"\n",
    "print(dataset[0]['data'][1]['type'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Here's a function to convert OPTIMADE's structure data to POSCAR. \n",
    "\n",
    "# Make sure that 'lattice_vectors', 'species_at_sites', and 'cartesian_site_positions' are\n",
    "# included in the response_fields of query URL\n",
    "\n",
    "def get_poscar_from_optimade_structure(structure):\n",
    "    if '_oqmd_entry_id' in structure['attributes'].keys():\n",
    "        poscar =  [\"REST API StructureID {}, OQMD Entry ID {}\".format(\n",
    "            structure['id'], structure['attributes']['_oqmd_entry_id']\n",
    "        )]\n",
    "        filename = \"ID-{}_OQMD-EnID-{}.poscar\".format(structure['id'],structure['attributes']['_oqmd_entry_id'])\n",
    "    else:\n",
    "        poscar =  [\"REST API StructureID {}\".format(structure['id'])]\n",
    "        filename = \"ID-{}.poscar\".format(structure['id'])\n",
    "        \n",
    "    poscar.append(\"1.0\")\n",
    "    \n",
    "    poscar += [\" \".join([str(jtem) for jtem in item]) \n",
    "               for item in structure['attributes']['lattice_vectors']\n",
    "              ]\n",
    "    \n",
    "    elems  = []\n",
    "    counts = []\n",
    "    for item in structure['attributes']['species_at_sites']:\n",
    "        if item in elems:\n",
    "            assert elems.index(item) == len(elems)-1\n",
    "            counts[-1] += 1\n",
    "        else:\n",
    "            elems.append(deepcopy(item))\n",
    "            counts.append(1)\n",
    "    poscar.append(\" \".join(elems))\n",
    "    poscar.append(\" \".join([str(item) for item in counts]))\n",
    "    \n",
    "    poscar.append(\"Cartesian\")\n",
    "    \n",
    "    poscar += [\" \".join([str(jtem) for jtem in item]) \n",
    "               for item in structure['attributes']['cartesian_site_positions']\n",
    "              ]\n",
    "    poscar = \"\\n\".join(poscar)\n",
    "    return (poscar, filename)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Call the OPTIMADE structure -> POSCAR conversion function \n",
    "# and save all structures in directory \"./input_poscars\"\n",
    "\n",
    "# Also save the bandgap values in a file \"target_properties.csv\"\n",
    "\n",
    "poscar_dir = \"./input_poscars\"\n",
    "if os.path.exists(poscar_dir):\n",
    "    shutil.rmtree(poscar_dir)\n",
    "os.mkdir(poscar_dir)\n",
    "\n",
    "properties = []\n",
    "\n",
    "for dt in dataset:\n",
    "    for st in dt['data']:\n",
    "        poscar, filename = get_poscar_from_optimade_structure(deepcopy(st))\n",
    "        target_value = deepcopy(st['attributes']['_oqmd_band_gap'])\n",
    "        properties.append(\",\".join([filename,str(target_value)]))\n",
    "        with open(os.path.join(poscar_dir,filename),\"w\") as fout:\n",
    "            fout.write(poscar)\n",
    "with open(\"target_properties.csv\",\"w\") as fout:\n",
    "    fout.write(\"filename, _oqmd_band_gap \\n\")\n",
    "    fout.write(\"\\n\".join(properties))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Part 2: Feature Generation"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### We need to generate a set of material representation vectors (input features) from the POSCAR data for Machine Learning\n",
    "\n",
    "#### The set of features generated for this tutorial is given in the file  \"features_processed_final.csv\" in the Github repo\n",
    "\n",
    "\n",
    "### Optional\n",
    "#### An example set of steps to generate Magpie material features is shown below. The following commands are to be executed on a bash shell from the same directory where this jupyter notebook resides. \n",
    "\n",
    "$ git clone git@github.com:tachyontraveler/magpie_workflow.git\n",
    "\n",
    "> Cloning into 'magpie_workflow'...\n",
    ">\n",
    "> remote: Enumerating objects: 115, done.\n",
    ">\n",
    "> remote: Counting objects: 100% (6/6), done.\n",
    ">\n",
    "> remote: Compressing objects: 100% (6/6), done.\n",
    ">\n",
    "> remote: Total 115 (delta 2), reused 0 (delta 0), pack-reused 109\n",
    ">\n",
    "> Receiving objects: 100% (115/115), 16.03 MiB | 3.25 MiB/s, done.\n",
    ">\n",
    "> Resolving deltas: 100% (31/31), done.\n",
    "\n",
    "\n",
    "$ cp -r ./input_poscars magpie_workflow/\n",
    "\n",
    "$ cd magpie_workflow\n",
    "\n",
    "$ python3 workflow.py \n",
    "\n",
    ">Initializing the workflow class\n",
    ">\n",
    "> 2021-11-22 19:51:45.614116 ::  Generating property.txt file\n",
    ">\n",
    "> 2021-11-22 19:51:45.616918 ::  Generating Magpie input commands file \n",
    ">\n",
    "> 2021-11-22 19:51:45.617025 ::  Magpie input file created as ./OUTDIR/generate-attributes.in\n",
    ">\n",
    "> 2021-11-22 19:51:45.617060 ::  Calling Magpie with the input script.. \n",
    ">\n",
    ">2021-11-22 19:51:45.617102 ::  May check out.workflow.txt for Magpie messages\n",
    ">\n",
    ">2021-11-22 19:52:36.347243 ::  Finished Magpie generation. Now post-processing\n",
    ">\n",
    ">2021-11-22 19:52:36.399678 ::  Done\n",
    "\n",
    "$ cd ../\n",
    "\n",
    "$ cp magpie_workflow/OUTDIR/features_processed_final.csv ./\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Part 3: Machine Learning\n",
    "\n",
    "#### Now we can proceed to build a small ML model using Scikit Learn library modules and fit it on the data obtained above"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Import the required libraries\n",
    "\n",
    "import numpy as np\n",
    "from sklearn.model_selection import train_test_split\n",
    "from sklearn.pipeline import Pipeline \n",
    "from sklearn.preprocessing import StandardScaler\n",
    "from sklearn.decomposition import PCA\n",
    "from sklearn.ensemble import RandomForestRegressor as RFR\n",
    "from sklearn.svm import SVR\n",
    "from sklearn.model_selection import GridSearchCV\n",
    "import matplotlib.pyplot as plt "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Load the data generated in previous sections\n",
    "\n",
    "# Input feature data - generated using Magpie on POSCAR files \n",
    "features = open(\"./features_processed_final.csv\",\"r\").read().strip().split(\"\\n\")\n",
    "feats_title = features[0].strip().split(\",\")\n",
    "\n",
    "features = np.array([item.strip().split(\",\") for item in features[1:]])\n",
    "features = features[~np.isnan(features[:,1:].astype(float)).any(axis=1)]\n",
    "I = features[:,0]\n",
    "X = features[:,1:].astype(float)\n",
    "\n",
    "\n",
    "# Target data (bandgap values)\n",
    "targets = open(\"./target_properties.csv\").read().strip().split(\"\\n\")\n",
    "targets_title = targets[0].strip().split(\",\")\n",
    "\n",
    "targets = np.array([item.strip().split(\",\") for item in targets[1:]])\n",
    "targets = dict(zip(targets[:,0],targets[:,1].astype(float)))\n",
    "\n",
    "Y = [targets[item] for item in I]\n",
    "\n",
    "# The feats_title and targets_title are not used anywhere else in this tutorial. \n",
    "# But I'd recommend keeping them for tracking features and targets in an actual ML study\n",
    "# For example, they can be be useful in feature importance analysis or in multi-target modelling"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Split the data to train and test sets\n",
    "\n",
    "xtrain, xtest, ytrain, ytest = train_test_split(X, Y, \n",
    "                                                test_size=0.2,\n",
    "                                                random_state=0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Create a Scikit Learn pipeline with feature scaling, dimension reduction, and fianlly, a regressor\n",
    "\n",
    "pipeline = Pipeline([\n",
    "    ('scaler', StandardScaler()),\n",
    "    ('pca', PCA()),\n",
    "    ('svr', SVR())\n",
    "])\n",
    "\n",
    "# Parameters to search in finding the best model\n",
    "# The following sets of parameters are just for a shallow search for optimization\n",
    "# A finer search for parameters will be required on an actual ML study\n",
    "\n",
    "params = {\n",
    "    'pca__n_components': [200,250,xtrain.shape[-1]],\n",
    "    'svr__C':[0.1,1,10,20],\n",
    "    'svr__kernel':['rbf'],\n",
    "    'svr__epsilon':[0.1],\n",
    "    'svr__gamma':['scale']\n",
    "}\n",
    "\n",
    "gridsearch = GridSearchCV(pipeline, params, cv=3)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Traindata score: 0.9784009665385437\n",
      "Testdata score: 0.7773354091561754\n"
     ]
    }
   ],
   "source": [
    "# Do the grid search and get scores\n",
    "\n",
    "gridsearch.fit(xtrain, ytrain)\n",
    "\n",
    "print('Traindata score: {}'.format(gridsearch.score(xtrain, ytrain)))\n",
    "print('Testdata score: {}'.format(gridsearch.score(xtest, ytest)))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### (I'm aware that the test data score is much lower that of traindata - implying a possible overfit.  But this would work just fine as a representational ML modeling workflow for this tutorial session.)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Pipeline(steps=[('scaler', StandardScaler()), ('pca', PCA(n_components=200)),\n",
      "                ('svr', SVR(C=10))])\n"
     ]
    }
   ],
   "source": [
    "# See which set of paramters had the best fit as of now\n",
    "print(gridsearch.best_estimator_)\n",
    "\n",
    "# Further, a finer search for most optimum set of paramters is required to get a better model.\n",
    "# But for the sake of this tutorial, I'm gonna continue to plot the testdata predictions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Get predictions on test data\n",
    "\n",
    "ytest_pred = gridsearch.predict(xtest)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYMAAAEcCAYAAAAlVNiEAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/Il7ecAAAACXBIWXMAAAsTAAALEwEAmpwYAABWs0lEQVR4nO2deXhU1fn4P+9MZrJDNlYVQsUFUVABq/60YG2t2q+7qBX3Km7UWtHaqsV9qda6a4XWDfdYq6LW4kYUV8AFARdUAsiSBLKQPZOZ8/vj3omTyUxmbjKT9f08zzzJnHvvue9M4Lz3vKsYY1AURVEGNq6eFkBRFEXpeVQZKIqiKKoMFEVRFFUGiqIoCqoMFEVRFFQZKIqiKKgyUBRFUVBloCiKogAp8Z4oIoXAwcAIYBjQApQC64HXjDFbkyGgoiiKknwkVgayiBwFXAlMAiTKaS3AW8CVxphlCZVQURRFSTodmolE5FbgYeBt4DBgDJCFtaPwAIOBXYDpwA/A2yJySjIFVhRFURJP1J2BiJwOnAScaIzZFtdkIqOBF4EzjTGfJkxKRVEUJal0pAwmAquMMT5HE4rkATnGmO8TIJ+iKIrSDXToMxCRN40xB3ejPHFRUFBgCgsLe1oMRVGUPsWyZcu2GGOGRDoWK5poiojcADxqjFmdeNE6R2FhIUuXLu1pMRRFUfoUIrI22rFYeQabgC3AsyLynojMFJHBCZVOURRF6XFiKYO/G2PuNMbsBVwA7AosF5GnReQwEdGkNUVRlH5Ah4u5MebBkN8/N8ZcAvwEeBw4C/hORG4Tkd2TK6aiKIqSTBw/2Rtj/MaYl40x04G9gXTgcxFRI76iKEofJe5yFKGISA5WDsIZwBSszOShCZNKURRF6VY6VAYiMt4Ys9L+3QUcCpwOHAGkAQ3AU8CjwBvJFVVRlP5AcUkJ9y5ZwprKSsbk5jJryhSmaqh4jxPLTDRfRMaLyG1Y5SYWAMcDHwO/BYYZY04xxrxuYhU5UhRlwFNcUsLshQspq6tjWFYWZXV1zF64kOKSkp4WbcATy0w0AViOZQb6DrgWeMwYU5JkuRRF6Yfcu2QJmV4vg1JTAVp/3rtkie4O4uGFF2DrVvjtbxM+daydQTPwL+BAY8xOxpjrVBEoitJZ1lRWkuX1thnL8npZU1nZQxL1EUpL4YQT4Jhj4F//gkAg4beIpQy+NsbMNMa8l/A7K4oy4BiTm0ttc3ObsdrmZsbk5vaQRL0cY+Cxx2DcOHjxRbjhBiguBlfiU7xizbhPpEERGSUie9m/dyoiSVGUgcesKVOoa25mW1MTAWPY1tREXXMzs6ZM6WnReh8lJXDooXD66bDbbvD553DlleDxJOV2sZLO2lQsFZFTROQ7YA3wsj38uIjcKSLRGt8oiqIAMLWwkNsPOYShmZmU1tYyNDOT2w85RP0Fofj9cPfdsPvu8P77cO+98M47sOuuSb2tk7aXpwDzgBewGt4Em9jMBu7A6oZ2Q4LlUxSlnzG1sFAX/2isWmU5hz/8EA47DP7xDxg1qltu7cTwdClwsDHmN8aYG4B6AGPMBuBM4JgkyKcoitL/aW6G666DvfaC1ath/nx45ZVuUwTgLAM51RjzfqQDxpg69R0oiqJ0go8/tnYDK1bASSfBXXfB0O4v6OBkZ5AqIsMiHbDHsxMjkqIoygCgrg5mz4b99oPKSnjpJXjqqR5RBOBMGTwPvCsiZ4rILoBbRLYTkV8DrwJPJ0VCRVGU/sabb8Iee8Df/w4zZ8LKlXDEET0qkhPTzlXAOKwkNIOVlbzOPvYycHViRVMURelnVFbCpZfCQw/BTjvBokUwdWpPSwU4UAbGmEbg1yLyC+BgoACrC9rrxpi3kiSfoihK/+D55+HCC6G8HC6/HK6+GtLTe1qqVhw7fY0xb6AVShVFUaISWpl1LxFuev55hrz2Guy5pxUltPfePS1iO6IqAxH5P+AjY0y5kwlFZDcgyxjzcVeFUxRF6WsEK7NmejzMWLaMsx9/HG9zM9//8Y/85IYbkpZB3FU62hnUAm+LyGnGmE/imcxWIH8Hfp4I4RRFUfoa9y5ZwtiqKi6fP5+9Vqxgxa67cstpp9E0dixFvVQRQAfKwBizSERuAhaLyMfAIuBbLD9BI5YDOR0YBuwK/Mr+/ShjzA9JlltRFKX34fezX1ERF7z0EgGXi/vOOov//fzn+EUo7eWVWTv0GRhjnhSRj4CLgbOBkcFD9s9gPaKvgCLgTmNM7/7EiqIkBO1YFsaKFXD22Vzy0Ue8P3Ei884+my35+QDUNjX1+sqsMR3IxpjvgN8Bv7PzC4Zj9Tv2A6XAD8aYtUmVUlGUXkWrXdzrbdOxbEAWnWtqgptvhptugsGDWXXXXVyUmkpmaipZxlDb3GxVZp02racl7RBH0UTGmK+BrxN1cxE5Dvg9kAnkAhXAXcaY+Ym6h6IoiUc7ltl8+KFVSmLVKpgxA+68k90KCrg9fNc0bVqv/156rJ6QiPwBOBU40hjzg4h4gEexchhUGShKL2ZNZSXDsrLajA2ojmW1tXDVVVap6e23t8JFDz+89XBfrMzaI8pARAqBW4ADgs5mY4xPRC7lR7+Eoii9lDG5uZTV1bXuCGAAdSx7/XWrhERJiZVEdvPNkN33S7MlvndafJwKVBljloQOGmM2GmOW9pBMiqLEyYDsWFZRAWeeCYccAqmp8O67VuOZfqAIoOeUwf5AiYgcJyLvishXIvK+iJzVQ/IoiuKAAdWxzBh47jmr9eT8+XDFFfDZZ3DAAT0tWULpKZ/BDkAhVsOcY4Ay4DjgKREZYYy5sYfkUhQlTvqiXdwxGzdapqAXXrBKSARLSvRDHO8MRGSsiMwQkcvsnzt14r5pWBFElxljNhtjAsaYIuBF4AoRyYhw35kislRElpaXO6qQoSiK4gxj4J//tHYDr70Gt94KH33UbxUBOOuBnAH8A/gN4A455BeRJ4HzjDENcU5XY//8LGz8U+BYYDegje/AGDMXmAswefJkg6IoSjL47jvLQfzWW1Z56XnzrHLT/RwnO4O7gGnAn7HCP/cGfgFcaY/f5WCur6Lc398JuRRFUbpOSwvcfrvVdGbpUnjwQUshDABFAM58BscAe0aoO/SWiDwNLANmxjnXAuAkYAKwOGR8d6ABWOlALkVRlK6xfLmVPLZ0KRx5JNx/P2y3XU9L1a04eQJfH60AnTFmHbDewVzPAEuAG0QkC0BEDgSOB240xtQ5mEtRFKVzNDXBX/4CkybB2rXwzDOWs3iAKQJwtjNYJSITjTGfhx8QkYnAJ2Fj9xhjfhdpImOMX0QOBf4KrBSRRqAJmGWMmedAJkVRBjidLpj3/vtw9tnw5Zdw2mlWP2K7sNxAxIkyWAb8V0ReAFYB24BBwHgs38EdInJayPlHYBW4i4gxpgI4x6nAiqIoQTpVMK+21soVuPde2GEH+O9/4dBDu1Xu3ogTZfA3++d5UY7fG/ZeI34URUkqjgvmvfYanHsurF8Ps2bBjTf2mwziruJEGXwJHB7zLAsBXnEujqIoSvzEXTBv61b4wx+sDOJdd4XFi2H//btR0t6PE2Vwt5O+BSJydyfkURRFiZuYBfOMgWefhd/9DiorLWfxlVdatYWUNsQdTWSMebCj4yIyy8n5iqIoXaXDgnkbNsDRR8NJJ8Ho0bBsGVx3nSqCKHSmHMUQEfmpiPws9AVclAT5FEVRohKxYN4vfsHUhQutUhKvvw5/+xt88AFMmNDT4vZqnJSjGAw8Bvxf8sRRFEVxRpuCeatXw1lnQXExHHSQVUpixx17VL5Y9JZe0k52BrcAjcCvge+Bn9uvGcDHwHUJl05RFCUeWlqsYnITJljlpf/5T3jzzT6hCGYvXEhZXV2b0NjikpJul8WJA3k/YB9jTLOI1BhjioMHROQV4NmES6coihKLzz6zSkl88onlI7jvPhjZNxom9qZe0k52Bn5jTLP9exslYoypQdtVKorSnTQ2WpFBkydbzuKiInj++T6jCMAKjc3yetuM9VQvaSfKwIhI8FveKCKt2cUi8ltgcEIlUxRFicbixVZvgZtuglNPhVWr4PjjQaSnJXPEmNxcapub24z1VC9pJ8rgDWCZiIzBKld9p4hUiUgFVp+BuckQUFEUpZVt26zOYwceaBWZ+9//4OGHIS+vpyXrFL2pl7QTZXA1sCewzhjzX6yooleAt4CzjTE3JV48RVEUm1dfhd13hwcegIsvhi++sJrT92F6Uy/puB3IxpgmoDTk/X+B/yZDKEVRlFa2bLEW/yeesHIH3n8f9t23p6VKGL2ll3TcOwMR+ST2WYqiKAnCGHjqKRg3ziopcfXVVsRQP1IEvQknoaWjReShDo4boAxYATxjjGnpkmSKogxc1q+H88+HV16BffaBf/3LMhEpScOJMtiG1aoyDagCqoEcrCiiGqACyAOygTki8itjTEkCZVUUpb8TCFi9hy+/HPx+uOMOq8ic293TkvV7nDiQL8dKLNveGJNnjBljjMkFRmG1sTzPGDMYKMTqeva3qDMpiqKE8803MG0aXHAB/PSnsGKF5StQRdAtOFEGFwNnGmM2hg7afZHPB660368DTgP2SZCMiqL0Z3w+uOUWq5TEF1/AQw/BwoUwZkxPSzagcGImGmmMidi9zO5pPCrkvU9EtnZZOkVR+jeffGKVkvjsMzjuOKsV5fDhSblVbykI11txsjOoFZE5IuIJHRQRj4hcjeVTCI7tAmjRcEVRItPQAH/6k+Uc3rwZ/v1veO65pCqC3lIQrrfiZGfwR+AF4GIRWYnlRM4FdgMygCMBROQa4PdYfgRFUZS2vPMOnH22VW76t7+F226DJJdf6E0F4XorTpLOXhWRScCfgX2BScAm4FXgRmPMl/ap84H/AOsTLKuiKH2AaOaYxV98Qf3s2Rzy+utsHjaM0scfZ+KMGd0iU9y9kgcwTnYGGGO+AE6Occ53XZJIUZQ+S9Ack+n1tjHHzKmu5qc33URBdTXPH3YYDx51FFvr6ri9pKRbnsxj9kpWnCkDRVG6Rl9zYjqVN9Qcs7W+npr165nzn/9w5PLlrB45kr9ecgmrx47FC2Q2NXWbmWbWlCnMXrgQsHYEtc3NVkG4adOSfu++guMeyIqidI6+5sTsjLzB+vxb6+rYfeFCFt12G0etXMmcgw5iv/PO48OQXgPdaabpTQXheiu6M1CUbqKvOTE7I++Y3FxYt45L/vlPDv7qK5aMHs3vjzuOT3JzSU1JoaSqivyMDKD7zTS9pSBcb0WVgaJ0E33NielY3kCAm778kpE33YQxhiuOPJJ5++2HD/hJdjYbamqosev2q5mm96HKQFG6ib7mxHQk71dfwTnnsNPixVQceCCH/+xnrMzMJNvjYWxODgUZGXhcLrY2NFBaW2v5H6ZN0yf1XkTCfAZ2foGiKFHoTV2t4iGWvMUlJZz01FPcd/TR+PbYA98XX8Cjj5JXXMxfzz6bXfLz2Sk/n7z0dLY1NeFxu3ni2GNZOnMmRdOnqyLoZUiUChOIyM8czvVPY8zOXRcpNpMnTzZLly7tjlspSkLpL9FExSUlzJs3j6vnz2en9et5a8oUbjnhBK48/vjWzxN+7QE77MDi9ev7zGfvj4jIMmPM5IjHOlAGAaweBa1DYe/bYYzplvKCqgwUpQepr+fF3/yG/1uwgKqcHB4480w+mjyZbU1NDM3MpGj69HaXhOYfhIZ2akRP99KRMujIZ/AdcHbI+6HAFUARVgObYD+D3bGqlN6YCGEVRel5ou5gFi2Cc87hqG+/5bWDDuKRk0+mLjMTgOaWFt74/nsmz53b7sm/r0VSDUQ6UgYPG2OKg29E5GngaGPM2rDzXhSRJ4G7gceSIKOiKN1IpCziOS++yFMffcTIp56CHXfkmjlzeHvMmNZFfWt9PSvKykjzeNrkJASf/PtaJNVAJKoD2RhzU9jQuAiKIHjuGmB0IgVTFKVnCH2Kd4nwixUreO6aaxj+9NNw2WWwfDkHnXlmG+fy6q1WxfqxeXm4RBiUmkqm18u9S5YAVmRSbXNzm/v05kiqgYiTaKLhIrJHpAMiMgHLjKQoSh8nmEWcU13NH+++m7/cfjs1gwZx+p/+BLfeChkZ7TJ6WwIBxg8dSoGdUAZtn/zjiaQqLilhelERk+fOZXpRUa/NzO6vOMkzeAx4R0SeAD7lxxLWewO/AR5MuHSKonQ5Asnp9WNyctjr9deZ9fTTpDc2Mv+EE3j4l78kf/DgNueFZvROLyqirK6uzfHQJ/+g8mgjR0ieQbQCd+pg7j6iRhO1O1HEjeUknoXVv8BgRRjVAXcBVxtj/EmSsw0aTaQMFLoaheP4+rVrqTjtNPLeeYflO+3EfTNn8lVBQcx7dlXOoDIJTXDrKDpJ6RwdRRPFbSYyxviNMX8CtgemYe0GpgLbG2Ou6i5FoCgDiXD7fbgtPmHX+/1wzz0wfjx5n3zC6muv5Ybrr2fZoEFxFXXraiG4oGkqFHUwdy+Oy1EYY6qAd8LHReRYY8zziRBKURSLrkbhxHX9qlVW57EPPoDDDoN//IOdRo3iWYeydqUQXF8r1dEfcVyOQkRSRGQ7ERkV+gKuS4J8itKr6G4nZ1ejcDq8vrkZrr8e9toLvvkG5s+HV16BUaMSJn+89LVSHf2RuHcGIjIOmAfslzxxFKX30hNOzq42ZYl2/Z9zc2HyZPjiCzjpJLjrLhj6Y0BgJKczkLRSGrEczEryceJA/gBoAhYAFUAg9DBwrTGmW3IN1IGs9AQ95eTsajTRXR9+yG3vv09lQwMj3G6eWb6cSc88AyNGwAMPwBFHtLtfuDN4c00NBhiRna3lJPownS1HEc5wYEdjTCDSQREZ1hnhFKWvkMgs2ngW+PBzOrPwFpeUMH/5cnbMy2P/1auZNW8e25WXs3HGDEbedx+EhYtC5NIRq8rKQIRdCgpax4LnqjLoHzjxGXwdTRHYvNxVYRSlN5OoLNp42knGOide38W9S5YwrKWFKx99lJtvugncbi66/HJ+f9RRERUBRI7safb7afa3DRjUaJ/+hZOdwVUicgfwELDaGNMYdnw+VgJapxCRd4EDgDHGmJLOzqMoiST06TzL62VTTQ3Qtabq8RRt6+gcIG7fxei33uLPTz9NzrZtPHfEETx13HE0ejyUdrCIR4rs8brdINLmPI326V84UQYfYyWaXQQgYf8wuoKIHIelCBSl1xBqO09xuVhZXk5NUxOVjY3kpqUxcfjwTjk54zE3dXROXBVAN2+G3/2Ovz33HKtHjeL6yy7juzFj2Fpfz+rSUloCAaYXFUU0TwWdzlUNDZTV1bXuhganprKtqalLilDpvThRBqXAP6IcE2BmZwQQES9wC/AqcHhn5lCUZBBcdJv9fr7asgW3y0WGxwPGMCg1tdPRNLFi6otLSthcW8vXW7aQnZpKYU4O+RkZred0pCiK16xhxd/+ximPPEK6z8eS887jkgkTSEtPp7mujhVlZQCMHzo06o5iamEhp06YwDWLFtFiDFmpqQzNyKCxpQVA21b2U5wog0+MMddGOygiYzspw4XAEuAbVBkovYjgovvppk24XS5SXJaLramlpTWLtzOLYUfhosHdSH56OjXNzdT5fKwsK6MwJweP282sadO4d8mSiMpkbHU17sMO48Kvv+b9wkIuOf54KkaN4sK99mLx+vW88f33pHk8jM3La1NQLtLnWLx+PROGD48cOXXGGY4/s9L7cVKO4tcxTjnD6c1FJA+4DPiz02sVJdkEHcb1Ph9u2yzaEgiQ4fF0yXnaUemGe5cswef3U15fT4vfT1NLC75AgK0NDa3nhCdo1TY0cOSrr/LwVVcxcc0aLj36aI457zxWDxnCppoanv/qK4qmT2fH3Fz22W67qJVFQ9HyEAMPx+UoOuBjnDuQ5wCPG2PWJtIHoSiJIPgE73G78RsDxuAPBCjMz++y8zS0dEPQST174UK+3rIFfyBAutdLVmoqLYEALX4/mV5v6/mhCVqyciV3PPkkO69ezf922YVLjz2Wsrw8wP7PnZLCJ5s2Mb2oiO8qK1lbXd1mZxDtc2h5iIGHI2UgInsDvwV+AqSGHXZkJhKRnYATgHFOrlOU7iK46M5ZtIilGzaQ4fGwc0EBHrc7Yc7T8KzmL0pLaTEGbyBAim2aagkEqAsLaZ06ciRTV62Cm26yQkSfeIKT1q1DRPCEnNfc0kJDSwtldXXsmp/PirIyVpSWMn7oUFJTUqhrbuaAXXdlelFRm5yHrmY+K32PuM1EInIY8C7W0/8BWE5jAUZgVTH9yuG9/wrcYoypjvP+M0VkqYgsLS8vd3grRek8QzMz2X7QIFJTUtja0OC4ImdHXL1oEWurqvh882Y+3bQJlwgC1Pt8GGNoCVipPRmekCX+o49g773hmmtg+nSr0NzJJzNx+HB8fj8tgQDNfj+VDQ3U+nwI4PP7KcjMZI9hw0hLSeHrLVsYmpnJqRMmMH/58nb5DECXqpAqfQ8nO4M5wMHGmA9F5FNjzEHBAyIyHdg33olE5EBgd+DEeK8xxswF5oJVjiJuqRWlk4Q+tY+1TUPB4mmJWBSLS0pYsmEDqSkppKak0OT3W7sClwufMTT7/WR4PGyfnc3OBQVQVwdXXWXVEdpuO3j5Zfj1j668a6dNY+aCBWyuq6OmqQm3rVjSUlJYVV7ObkOGkJ+RQW56OqW1tRRNn870oqKoYapF06fr4j+AcJKBnG6M+dD+vY2B3xhTBOzlYK5fAm5giYh8JiKfAefZx161xzSySOlRutpLIJ75U1wuapqaKK+ro7KhgUAgQJPfT05aGv9v1Ch2ys/H43ZzVVMT7L473HknnH8+rFzZRhGAZdaae8QRZHu9ZHq9DMnMJDctDW9KCm6Xi5KqKqCt7V8dxUoQJzuD0Fz0FhEZYYzZBCAiOcCu8U5kjJmDtdNoRUSuAa4GDtcM5PjpahEzJTqJrEUUic83b6bZ3g0ECf42MiuL0tpadvd4uPmNNxhRVAQ77wzvvAMHHhh1zqmFhQzPymLi8OG4RNhSX8+X5eW4Xa7WCKRQ2786ipUgTnYGG0TkOhFJBYqBhSLyexH5PbAIWJUMAZXoxFPjRuk8iapFFI16nw9fIHK5r9z0dJbm5zP3j39k6L//zUOHHsrJc+ZQvMMOjuQuyMhg3JAhuERwuVztbP/aR0AJ4kQZ/B0YYr9uBGqAO+xXFvC7zgggIodHMhN1Zq6BRrLNGAOdZC+UmR4PfmNaIzGCtteRNTXM/tvf4Pjj+S4jg5lXX82/Tz2VDT5fXMo+XG6v283owYN54cQT2/kButquUuk/xN3PoN2FVmLAWKwQ06+MMS2JFKwjtJ+BxeS5cxmWlYUrJEcjYAyltbUsndmp6iBKGMk0w00vKuL5VataTUNiDGd9+im3LVxIaksLz514Ig8fdBCZIUli8fZPUPOhEolE9TNog7G0yOpOS6V0GbX3Jp+u9PWNxawpU3ht9WrqfT7GVlbywEsv8fM1a3insJB/nHsu3+TlMSw9vc018foskim30j+JaSYSkRwROUxEfiUiKfbYgSKyRETqROR7Ebkk+aIq4ai9t28ztbCQG6dO5fKPPuLT++9n0saNXHzUUZxz8cWce9JJSfdZKEooHe4M7Izj/wIF9tBSO6dgAdAALAe2B24TkQpjzCNJlFUJQ/vG9i6imWbCxw/YYQcWr1+Pa8UKbnziCcZ+9x1LJ01izvTpZP7kJ8wNMemEZgH/UF3Num3bGJKR0a78tJqFlK7Soc9ARF7DyjB+GPAA5wLfA18AlwU7n4nIH4EjjTHd0pNAfQZKdxHvIltcUsLMBQuoaGig2e/H63aTl57OBVOmMH/58tZ+wl9v2cLGLVu46t13ueydd6jJyOC2k07i0CuuYOqYMVHvv7y0lLLaWkYNHsz2gwe36UEMtOtZrP2JlUh05DOIpQx+APYyxpTb73cFVgDZxpiGkPME+N4Y0/5fcxJQZaB0B5Eaw0dbZKc98gjLNm5ERKzcATtkNMvrZaJdCnprfT0pH37IP154gd22bGH+xIlcethhZI8YwV4jRnToFJ5eVNTOPxR0JgNRj3U0p+4mBh5dcSDXBhUBgDHmKxFZE6oI7HEjInHVGFKUvsKcRYtYW12Nzy4LMTonJ2ofg89LSwFoaGlBALcIfmOobmqiuaWF+ro6Tnn8cc778EPWDx7MoTNm8PYuuxAAtlZV4XZ17L6LlQDnNDkuvEBeR60zlYFBLGVQH2GsJsq5kbNnFKUPUlxSwtING0jzePC63TT5/XxZXs6uBQWRF1ljaPL7rXwBO9RX7PER77/Pnf/5D9tXV3PfPvtwxcEHU5uaSgrgEolYlTScWJFjTqPK4mqdqQwoYimDbLuoXGgtoqwIY2AlnilKnyTcZFJeV9daN8hgPel73G6+q6xkWoTFcsLw4Sxet641PM8AuXV13Pm//3HyZ5/xZUEBB5x5Jh+MGtV6TdCUFFQaHRGrpHS85aaDn/PVb75hUFoahTk5rb0NtCbRwCaWMtgRq9REKBJlTCuJKn2SoMnE5/dTVl/P11u3UtfcjFuktW6QH6sMdHNLS7vQ3eKSEgQr4c8ALmM4YdUq7nzlFXIbGrjhZz/jr1OnUu92R7y/x+WirqWF4pKSqE/lwcixqxct4v3168EYJgwf3uZYrKiyUNPQoNRU6n0+viwvZ9yQIRSE9FhWBiaxlMF6wgrKRUGAa7osjaL0AMFWk2urq3G7XKSnpFDb3EyLMVZ2tzEEsJ52AlhP4UGHK/wYyTMuP5/akhLueflljvr6a1aMGsWxRx/NFyNG4HW5yHC5qPf5WnMHXECm14tLhFGDBsVloqltbmbCsGGtO4BQO3+sa0NNQ2Nyc1lVXo4BSqqq8CawYY/SN4mlDD4zxjwaz0QickwC5FGUbmdNZSVl9fVtmt4HMcbgcbnwG4PfGAKBQBuHa5bXi8/v59vyco7/4AOuf/VVUgIBrjjsMJ446CBKm5rIEKGppQVSUsjweKhrbsYA6XYv5cKcHHLT02OaaLpq5w91QudnZLDbkCGsqaxkW2MjQzMzNUdlgNOhMjDGHBXvRE7OVZTexJjcXL7eupX0lMj/HQJYSsGF5Rz+dNMm6n0+PG43NU1NjK2o4In//IdpJSW8PWYMFxx9NN/n5rLvkCE0b91Kvc/XOldTSwsel4t0j4cDR49uHd/W1BTTRNPVktrhTuj8jAw8bndctY6U/o+TqqWK0i+ZNWUKKSI0+a2WHS2BgLXwAykuF7lpaVaEkAjGjhryut1ISwu/e+cdltxzD3tu3sy5Rx3FIWecwVeDB+MPBBiUmkphTo4Vaupy4XW5mDh8OIU5OQzNzHRcRiTe8hTFJSVMLypi8ty5TC8qaq1yquVLlI5QZaAMeKYWFnLNtGm4RKhvbsbrclE4eDAel4sMj4cmvx+3rQgyvF5SXC7Gb9rEwvvu47bXX+e1sWOZOGsWD02ahD/E4bylvr61n0CGx9OaCDb3iCOYd8QRjstGx7OYx+pxkeX1sry0lPfXrQPQvAKllU6XsO5JNANZiRcnWbbRaggFW0N+sH49eSJc/vbbXLhoEZUZGcw6/HCeGzeOFLcbXyDQJqxucGpqa6TOtqYmBBiSmdmljN+OPk9xSQkznn+e6sZGsu1dSb59b4C65mYtWTHA6XQ5it6KKgMlHiKVk9hcU8Pw7GzqbPOKkwX5d1dcwUX//Cc7lZfzzOTJPHTKKXzU0EBtc3ObnhL+QAARwet2EzAGj9sNxjA4NZUd8/OTshgHP2vQ9+E3Bn8gwG5DhpCbns7769YxwS6LESTe3ghK/6EjZaBmIqXfEt4Jzuf3s7GmhlXl5XG3CS0uKeHUxx7j2YMO4p6bbybd72f27Nk89oc/0Dh4MEMzMxmRnU3AjjQCK0po++xsKy/B7yfL48HlclmlKfz+pHSlC37WbK8XvzGkuFy4XS5KqqosP4OINr5XOiRhysBuaK8ovYageSdISVUVHrcbX5wLcnFJCc/ffju3XnwxxxcXc9/++zNp1izeHDu21dY/74gjePyYYxielUW610tBRgZ7DBtGcyBAhtfL8KwsJm+3HQJ4U1JYW1XVOn/4YhzN8evks47OycEfCNASCOAWocb2K0wcNkx7IygdEjW0VER+5nCuk9HEM6UDurtKZngoZb3Ph9vlIiMkhDTq0/GWLcipp3LX4sV8OXQoZ5x/Pp+NGUNjSwtrq6p44aST2sj+xLHHtjFJrSgrwwUU5uQAkOHx0NjS0ibMNHQx7mrhuOBnDTqs11ZVUdPczOC0tDZlroOfuaOSFcrApKM8g0W0LTGhJSeUTtNdVTJDFU6m18vmGquuYpbXi8florGlhZ3z81vPb/d0bAw8/TRcdBH7VVbyt1/+ktunTaPB5aKxsbG1ntCcRYsoPuOM1svCS0IMTk0lPz2dfLvuT2FODl+UlpLm8RAwpt1i3NWEstDaRXnp6a0ZxaHfrzZCUjqiI2XwHXB2yPuhwBVAEVZPg2ogB9gdOA24MTkiKv2B7qiSGa5wau1MXwFKa2sZP3QoG2tqWh274QvyBx98gGfWLCZ/8gmrx47lT+edx4seD6alBWMMglVlVERYumFDu1pCoSUhgrJsa2qyFJHbzcjsbIZnZ1NaW9tuMXaSUBZthxVrsU9EX2TtgdB/6UgZPGyMKQ6+EZGngaONMWvDzntRRJ4E7gYeS4KMSj+gq9mzkQhfmMrq6iIqnCGZmVw7bRpzFi1iU00NJVVVpKaksPeIEdaT86hRfHPjjexxww24AwHmnXIK9++7L19WVmL8/tba7AbAGDI9HlLd7g4VWaTF+bqQXcCayspWX8XUwsKYJapDP3NHO6xkLszaA6F/E1UZGGNuChsaF0ERBM9dIyKjIx1TFIhdj98pkRampRs2sMewYRByj2CS1TkLFrCppgZvSkqrE3lzTQ3p338Pp5/Ozu+8w9LdduMf55xD6bBhlG7cSGpKCqluN9tsx2twZ+ASYWxeXkxFFr44d7SYxipRHbx+xvPPU93URLbtLA6Wn+6OPgTaA6F/4ySaaLiI7BHpgIhMwDIjKUpEEl0KITxsNNhqcunGjSzbuJEt9VZfpuCiWtnYiDclhRSXC4/bTYYIv3nlFSb86ldsW7qUC487jlkXX0zpsGGA5WxOdbsROyQz6DDzG8N22dl4U1IcK7JwmUOjmYI7iWhZyUFFUt3YSHpKSmuznS319d0WIhoenQUantqfiFW1NJTHgHdE5AngU6AKyAX2Bn4DPJhw6ZR+Q7w190PpyD4danbaUl/Pl+XlpNqlp+t8PlaVl1slJdxuMjweapqbSbX7Cey+YQN/f/ZZJm7cyAvjx/Pc+efzv7o6GsvL2d3loiAjw6ou6vORYnciCyaVuURYV12NLxBoNfvESyxTWUdmntY8gtRUmvz+1uqqa+3y090RIpro3Z3Su3CiDP6EVXJlFnABtPrm6oC7gKsTLp3Sr4jHph1UAMs3b6asro5Rgwez/eDB7ezToQvTWruHsBvITU8n1a4murWhgSeOPZZ7lyxhUUkJnuZmLn/zTS4oLqY8I4PpJ57Ie1OmMCkvj53S0viitJRvKyrIS09nSEYG1RUViNtNakoKLhEaWlrwut2kpaQwMjvbsWmkK4tpUJEU5uSwqtxqS+4WoaYbQ0TjMWUpfZe4zUTGGL8x5k/A9sBBWLuBqcD2xpirjDH+JMk4oOhK4lFfJ7TI2ramJgLA2upqKhoa2iWIhZqd6nw+jDE02jH89T4fWampZHg8TC0sZNaUKfxywwYW3XknF739Nk/vvTfjL7yQVydMYLSdB5CfkcHuQ4fSEghQWlvLzgUF3HTwwYgILX4/WV4vU0aO5KAxY9hn++3bJXDFQ1dMZcGKpcE+BKluNw0+H4NTU7vNgRvLlKX0bbQ2US8iUi2dgVRMbHpRUeuT87tr15KakkJLIECq282kkSMJGENpbS1LZ84EftxFvPHddwTs5jOpKSm47XLULhEWHH44BzzwADzwABuHDOGCI4+keKedSHG5GJmVxShbGUDkWj2hMgVZV1XF1oYGhmdlOQ6v7Gxo5kD/t6Ekho5qEzkxEyEiO2OZi6YBGGN+IiLXYXVEe76rgg50elO0Rk/Ek4fa1ENLRwezdsNNKkGzU3FJCUc//TSIkOJy0WJXD52xdi27HnQQVFbCJZcw8rrreCEzs/XzheYBRDN5hJtGfqiu5tuKCsbm5XUqvLKz4Z+d8bkoihPiVgYiMgV4G6gEvgJ2tA+9B9wpImKM+XfiRRw4JCMWvzP0VDx5qE09aBtvgdZeANHs01MLC1ubxTS0tLB9UxN/f/VVDvn4Y74dOZKCDz6AffZpd008i2v4eVsbGhibl9e6o+hOhZ3sPAJlYONkZ3ALlpP4DmNMQEQ+ATDG/E9EDgGeBlQZdIHeEq3RnTuUaOUjctPTGT14MOu2bSPb643Zo3fC8OGU1dZyxNKlzHzsMdIbGnjo6KN5/YQTeCpMEQSJd3ENPW/y3Lm9QmErSqJxogxGGWNuj3TAGLNeRNISJNOApbdEa3TXDiVW+YidCwq4+7DD4mpGU1hTw7UPPcS0Vav4cuxYbjnjDFYOGcLt++2XUJl7i8JWlETjRBl4RMRljAmEHxARD1CQOLEGJr3FLuxkweuMb6HV8fv996S4XIzNy2tNwgKrfMSikCJw0eaYvXAhWSkpnPXBB5zx9NMQCHD3ySfz+IEHMjo/n9uT4OfoLQpbURKNE2XwEfCciMw2xqwJDopIDnAnsDixog1MeoNdON4FrzO+hdBrAoEAARG+LC9vbQ8Z7w7k3iVL2HXrVv786KOM//prPtljD/562mkERo/m4yR27uotCltREo0TZXAplrP4WxEpAwaJyLdYeQcbgQOSIJ/SA8S74HXGtxB6TabXa0UMuVysraqiICMjPpOLz8e0J57gnFdfpcnr5Y7zzuOtAw8kAJR2g+2+NyhsRUk0cSsD2y+wJ3AJcDCWWWgL8CSWU1k9aP2IeBa8zvgWQq8JRgy5RKjz+TqMGGpl2TL47W+58PPPeXvKFB466yyqBg8GoLapqU/b7rU8tNKTOGp7aYypsLON9zPG7GSM2Q+4A8iKda3S/whmxYYS68k+9JpgNq1bBBd0nNHa0ACXXw4//SmUlbHiwQe57OyzWZeWFjGbt69lcodmX8fbn1lREkncykBEno1yaArwlYhclRiRlL5CZ8orhF/jcbsZnZPDCyedRNH06ZEVwaJFMGEC3HornHUWrFrF7jNnRi2N0NHC2luVREcVTRWlO4i7HIWIfGqM2SvKsWzgA2PM7okULhr9tRxFX6Qr0UQxr6muZuP55zPyqadYP2QI/zr/fA4688yY80cqIbGtqYnqxkbWVlfT4vfjcbkIAD6/nynbbce1PewEDuYvBKujAu3KbyhKV+l0OQoRGYTV2hKs0NIdsMLA25yG5UTO6KKcSi8m2gLeGWdqXNe89BJNM2cyrKyMpw49lKITTmCrCC/HkQkdyZfR3NLCirIyMr1e3C5Xa8Oa9JQUVpaV9XjHLs1fUHqaWGaiPwAlwBpgXMjvoa/vgXeAN5IlpNKzdKs9u6wMTjoJjjqKTWlpnP+Xv/DkaafhS0vr0HQSav7ZXFvLD9XVbY5/W1EBQKrbTWNLC2J3LPMFAvgCgR43ySS6+Y+iOCVWNNELWApAgGuBORHO8QFrjDEfJFSyPkpXqlImO5Kks/folvIUxsDjj8PFF0NtLdxwAyfm5VEweHCbJ5ZI0Up3ffgh1yxaRIsxZHm9ZLjdrYv/9oMHU9vcTL3Px6DUVPx2dVMBqzx1IECGxxN3fkOy/k6av6D0NE58BjcYY3qFk7i3+gw6W2Y42eWJi0tKuHrRIpZs2ECG18uOubmkpqTEfY+k27PXroVzz4X//Q/23x/++U8YNy6q7T+0zHRxSQlHP/MMAWNIdbutxT4QID8tjaZAoLXMdFldHVUNDaytrqbe7n+A/XmmjByJx+1uV746HC0jrfR1OvIZOGlu0ysUQW+msxEhyYwkCS5gK8vKSE1JIWAMX23ZQrPfH/EekaJtOhNCGheBANxzD4wfD++9Z/3+7rswbhzwo+lkXVUVSzdsYNGaNSwvLeWAHXZoneLeJUto8ftJS0lB7BLWbpeLer+f4VlZLJ05k6Lp07lu2jQrcmnwYLK8XgKAMYad8/LwuN1xmWQ04kfpzzgJLT1cRD4RkQ/Cxt8QketEJNyxPODobMPwzl4XT5hkcAHzBQJ43O7WxXJtVVW7e0TzDRywww6Jt2d/+SUceCBcdJH1c8UKmDULXD/+k5xaWMipEyawrrqaWrt72ahBg5i/fHnrZw1+dy2BH0tm+QMByuvq+K6ysvV7CZphdi4o4Ce5uRwwahQHjBpFuscTd8cubQiv9GeclKOYCXxN+17HvwduAq4Cro9nIjuT+UKsEhYtgBvLAX29MabcgUy9is5GhHTmunjrAgUja4LNYlJcrtaGMeH3iOYbWLx+fRt7dpbXS6bXy+yFCyPazYtLSpizaBHLN28GESYOG/Zj6GZzs5UvcP31kJUF8+fDjBmtJptwFq9fz4Thw9uZioL+ijG5ufj8ftbaDuOAMdQ0NSEi7Jqf3+576Yo5RyN+lP6MkwzkMcBpxphvQgeNMSuBk4DjHMz1NJAHTDbG7AH8EjgEeE9E0h3M06vobERIZ64LLtzNfj+fbtrE55s3s7a6mjmLFrU5L2jiKczJwR8I0BII4LeTvcLv0dGT79TCQoqmT+f2Qw5pNRlFiiwqLinhnAUL+GTjRlwuFwIs27iRmQsW8MmLL8LkyfCXv8Cxx1q7g1NOiaoIYskU/O6C5h+vy9WqCHbJz6cgMzOhphyN+FH6M06UQYoxxhfpgDGmAfBGOtYBlxtj6uzrNwC3ATsBhzucp9fQ2YbhnbluTWUlZbW1LN24kbK6Ohp8PppaWli6YUMbc1FwAfO43YwrKMAF1DU3Y+zF7N4lS1rPj8c3ELp7qGhoYPXWrXy9dSvHFxUx9ZFHOPqZZ1hbVdVqv/e43QwOBLioqIg9jzkGKirgpZfgqadg6NCY32ksmULNP9sNGkS218vkkSP5SV5e6/mJMuVoQ3ilP+PETOQXkZ8aYz4KPyAi+wB+B3NNMMY0h41ttH/26T13V3rcOrku0+vls82bEcAtggHqfb7Wp+DgXOEhi+OHDmVjTQ0jsrPJ8nrbmFHiKV0dNDttqa/ny/Jy3LbZaWt9PZ9s3Ijf3nnU+3y4XS5+vmYNtxcVUVhRwb9/9jOOe+klsAvLxUM8MoV+d8EIpFASacrRiqVKf8WJMrgDeENEngSWAhVYC/dk4DdYvoO4iKAIAHYGDFYCmxIDwYqGERHLzGKHCLtF2j0Fhy+WBiLmDATNQB3Fugft5t9s3WqFaGI5bN0ieO1wVZcIg+rrufWllzhr2TK+y8/nqHPPxXvwwRznQBEEZY83/r64pITyujqWRgih1eYzitIxTkpYPywiecA1wDnQ2qGwDrjaGPNIZ4UQETfwW+Bf4T6J/kKik5Vqm5vJSUujprmZgDG4RayoIb+/w6fgWGWnYz35zpoyhXMWLKCqoQGXCGLvSgJYSiHF5eLXK1dyx4IFDK2r486pU7nxoIPIy8tjbidt6/E8jYc61CcMG8a3FRV8UVrK5O22U1OOosSBk50BxpjbRWQusD8/9jN43xhT00U5/oKVyXxxtBNEZCZWRBOjRo3q4u3iJxGLeGc6gsUiNIomYAxNLS1sa2zE7XK1icOPdF1XImKmFhYyMjubtVVV+G0l5HG5MEB2RQX3/+9/HLF8OctHjODoGTNYXVjYNpooBgnJkk5NpSAzszVBTRWBosTGkTIAsBf+/4WPi8hZxpiHnM4nImcCJwDTgg7lKPedC8wFKwPZ6X06Q6IW8WSUcwja0vPT0li3bRsBwOVysUN2NvOXL2fP4cMjzp2IHr61zc3sNWIEX23ZgtvlIuD3c8xHH/G3hQvJ9PmYe9xx/Ovgg7ktrJl9MBP689JSMIYJw4dzXYiS6Mr37bTRjjaSUZS2OEk6G9XRC6sDmiNE5FRgNvBzY0yZ0+uTTUcZp07q4icjWSloS28KBEjzeMhJTWVQaipb7TLN4SGm4dcBvL9uHctLS/H5/Vy9aFHcNf7H2Lb4cUOGsGN1Nf9+6CEeevFFvh42jJP/8hdeP/FEbo2gCGYuWMCyjRsRLMX1ycaNnLNgQev9upLh6yRLWhvJKEp7OlQGIrJERIrstyW0r1ga+hrn5MYicgpwOfALY8xme+z/bHNQryDaIv755s2OFpNklXOYWljI8Kwsxg8Zgt92IHvdbgLGtAsxDaeuuZkJw4ezS34+31dWsmzjRjwuV1wL46wpU2hobOQ3b77JO7ffzr4//MDNJ59M85tv8u85cyI2qbl3yRIqGhrwpqS0ZkJ7U1KobGxsXey7ojSd5ABoWQlFaU+sncFirMghgPXAWVFev7WPx4WIzADmAY8AvxCRU2zlcAQw0oH8SSXaIh4M4Yx3MUlmstKY3Fy+q6zE7XKR4nIRrAqS4fFElefeJUvw+f2s3rqVZZs20eT3IyKsra6Oa2HM+OYbHr7tNi566imKx4zh9BtuYP8bb2TqT34S9Zo1lZU02xnQQdwiNPv9rYt9V5SmkxwALSuhKO3p0GdgjPlDyNsFxphHo50rIhEr4UXhHiANK9EsnGsdzJNUotnXM+2Sx6F0tJgkszzxrClTePWbb0hNsf6ULYEA/kCAnQsKosrz4fr1bKqtBazyDQI0+HwE7N1FpJpF9y5Zwvrycs5duJAZL79MfXo6t15wAa9MnkydL2IuYhvG5OaytqqKFjviCMBvDF63u3Wx76o/I94cAC0roSjtcRJaOivGKXE7j40xebHP6nmiLeL3LlnieDFJVrLS1MJCpmy3HSvLymhqaSHD46EwPx+P2832mZntzi8uKaHUTsoKOn+DJd78drG30M8StK9PXr+e2x55hMIffuCZiRN55IwzSBk2jEEAIjGd4bOmTGH55s1srKnBuN2ICM0tLYzIzm7dIXVXTf9EONEVpb/hOJqoA/4J7J3A+XoF0RbxcxYsYFV5Oc1+P163m9y0NK7rocXk2mnTmL1wIT6/n7K6Or4oLSXF7eaaqVPbnXvvkiW4RQgYQ8AYXFip4wFjaPb7WbRmTZtr5737Ln/89785/vXX2Zqbywmnn86bu+1Gqt/PJHvOeEwsUwsLmXvEEa3RRCYQYO+RI9tEEwXPS3ZUjzaSUZT2RFUGIvK9w7l6ja2/OxCwsn7tV0/W7w6Wem7t9pWaytCMjIghpmsqK8lJT6e2uRmf34/f5cIVCBDA6vwVeu20b7/llksuYfstW3jll7/kjiOP5L2qKpobGnC5XGytryc/I8ORXX/RGWck7XtwgpaVUJS2dLQzGAy8FDZ2OFAFrASqgRxgNyxF8FTixeud3LtkCcOzs9m5oKB1LLSsck8Qq9RzkGCyWoPPh9frxS1CZWMjYgz7bLcd+RkZZNbWctpjjzFx8WI2jhjB7/78Z5btuCOrysvxuN34bHPSyrIyCnNy8LjdnTaxaLy/ovQOOlIGq40xZwbfiMhlwAd28lcbRORcIHraaz+jswlOyzdvps7nI8PjYeLw4Qld+MJl2lpfz5rKSpZu2MD0oqJWu3x5XR3fVlTgcbtpbmmhoaUFvzFkezwA7P/xx5z7yCMM3raNhw49lLF33skHxcWs3boVlwhul4tUY0hzu2kOBNja0MATxx7bqc+RjMxsRVE6R1RlYIzZN2zoWGPMflHOfVBElmA1uOn3OIlGCS54Pr+fzbW1IEJNczPfbNnieOHr6Ck6VKat9fWsKi+3CtKlpVFWV8fMBQswwIjsbCYMG8aqLVuo8/nI9npJcbnIr67mT08+yZErV/JtYSGX/eEPbBs3jqJdduH21FSOfvppAkBmSgo75+dTkJHR2ge5swt3MjKzFUXpHE76GewoIhGVh4h4gdGJEan305kEp/L6elLcbtJSUkhxuSivr3eU6BQrazZUpjWVla1VBAtzclp7D1Q2NjLIrtuT6fGQnZpKpsfDJV99xcd33MEvv/6a6w8/nJlXXcXn9s4FLPv6L3bckT2HD2fSyJEUZGQAXQ/HTEa8v5PMcEVRfsSJMlgOLBCRyXaVUUQkxe5l8CLwWRLk65V0JsGp3udrja8Ptp10svDFypoNlWlbUxMZHg/jhgxpXbib/X6a/T+2nKhubGREaSmPPfAAVz7yCN9svz0Hz57NLfvvT/6gQe0+zwE77MDy0lIWrVnD0g0bWFdV1eXEuURnZmuZCUXpPE5CS88HXgc+AhCReiDDPlaC1bpywOA0wSm0B7HfGDI8HkcLXzx+iqBMwQYvoWYsr9vd2l6yoraWC959l2vfeosWl4tLjjmGhydPZlRuLocXFFA0fXqb+xSXlDB/+XJGDRpEWX09tc3NNFZXc00XwzETHe+vZidF6TxOks5Wi8jOwBnAfsAwYBPwAfBotJaYA53ggjckI4MSOwMXYPvsbEcLnxM/RaRFNi89HQMUfPst186bx97r1/Pyzjtz6THHUFlQQKClhXXV1dx92GFAW//E5tpa8tPTGZWTw6icHMCKVFq8fj2/3zfctRQ/8cT7O4k2curYVxTlR5z2M2jGKiPdLqJIiUzogtfo91PX3EyGx8POBQWOoomcPEVHWmSv328/Rt1/P9vfdx+VaWmcN2MGi/fbj22NjTT5fGSkpDAoLY2phYXtony+3rKFmuZmMrzeVrNTIvsKR/sOnEYbaZkJRek8jjOQReRnwDQg3RjzZxGZCnySgAY3/ZZEJDg5zZptc88PPoBjjoEvv4RTTuHyadP4zu1mTGoqY+zzg41goL25JTs1lTqfj7VVVQlzHseDU7OPlplQlM7jpJ9Btoi8DizCan15un3oMGC5iIyJcumAJ1ERLlMLC1v7FAPMXriw4/lqa+H3v4f/9/+s3199FebP5/SDD+4wGio8yqcwJweMaW2xmciqqx3hNNrIiWNfUZS2iDHxNQ0TkfuAPYHrsDKQXzbG7GkfOwk4xhhzYnLEbMvkyZPN0qVLY5/YSRKZFRtq6gh9Wu3sIhX3fAsXwsyZsG4dXHgh3HQTZGfH9RkjOaDXVVWxtaGB4VlZ3ZYpHEmO4A4m3MmtKEpsRGSZMSZihWknZqJDgYnGmFp70mCxS4wxT4vIpV0Ts3eQ6KzYREe4xJyvogIuuQQefRR23RXefdfaGYTRkekqkrnF43Z3OtO4s6jZR1G6Dyd5Br6gIohCThdl6RUkugtWohOros5XUQFFRTBuHDzxBFx5JXz6aURFECSa+aq3mFt6ixyKMhBwsjOoE5HjjDH/Dj8gIocDFYkTq+dIdHhitAiXTK+X6UVFjk1RkeZLLS3l/iefhKVLYdIky0Q0cWKH88TaAcXj9O6OInNaXVRRugcnO4MbgGdFpFhEbgaGiMhVdo/k5+lFHcq6QnhW7Jb6ej7esIHvKis75fyNVLpiU00Nm2tqOpUp22a+QICfLVzIo1deyaQvvoDbboMPP6R48OCYDuuu7oA021dR+hdxO5ABRORE4FbaVihdB8yOtGNIFsl0IIc+MTe1tLCyrAyA3YcOxZuS0uqsBeJ+Kg5/gi6vq7OKyMXhGA2/9oAdduD5r76iZsUK7vzPf/jZ999Tte++5MyfD2PHxu1gnjx3LsOysnDJj50YgoXnls6cGfN7UueuovQ9OnIgO1IGIRPuAhQAW4wxX3dRPsd0VzTRG999R4rLxU75+eTb8fXbmpoQfjT1dCZCKN6FOHxh/6G6mjXl5Vz92WfMevllfCkp3HnccRx43XVMHWNF9sa7SDtdzMOV0uebN7NTfn6nlUmkObWXgaIkl4REE4lIJdAE7G8rgG5XAt1F0E49ee5cPC4XJVVVrCovJ8PjYfTgwXy1dSsThg3rdIRQvJmy4ZFDQ779lkeeeYa9Nmzgw0mTeODMMynJzOTzpUtblUG8Po8DdtiBa4qLafH7yfJ6GZqZGbVJTST/Qnl9PekpKa3lKaJ9hmhoLwNF6V048Rn4sUJLnbbD7LNkeb2sKCujye8nNSWFJr+fFWVlrQto+LnxOpnjLYEdjBzyNDdzyrPP8t8772T76mpOP/lkbrzkEiry8trdN55KoKGF57JSU6n1+VhXXc2pEyZEXIgj+RdGDRrEuurquMp4RyLRUVuKonQNJ8pglTGmNNpBETkqAfL0KkINaKHmtNSUlC6VXo43ZHJMbi6FK1dy1xVXcOILL/D83nsz5eKLeX3vvVsrkIbfNx5Fc++SJfj8fsrr62mwG9yMGjyYxevXR5Q3Ujjr9oMHMzQrq9Nhn8noZaAoSudxElr6oIj8CfibMaYlwvGrsfoa9BvqmpsZP3Qo66qrqbfbVe6Ul0dlYyN1tjLobDJUzJDJmhrufPllRsyfT2leHn+5/HJeGj2a8ooKxtpdxiLdN54aRp9v3kxpXR0pLhdet5smv5+S6moaWiL9WaObtSYMG9ZpZ7EWlVOU3oUTZXA2MA64TES+A8IT0MYmTKpeQnDBmjRyZOvYtqYmJgwbxqwpU+IuGueY//4Xzj2X7X74gU+nT+fEPfdknd+PVFTgdbn4pqKCr7ZsYUhmJpftv3+7+wYVTdBBO3vhwjYO2nqfVW082GwnRYSWQKB1PJxkZAJrdrGi9C6cKIPJQGgIj0Q7sb/Q0YLVmWSomNEzW7bAH/4Ajz8O48bx6XPPcc6WLaT7/aTavRBqfT5S3W7SPB7y09OZv3w5ew4f3k6Wjhy0mR4PNU1NtAQCpLhcVo8FY8gMM9sEcVoxNR6SMaeiKJ3HSaG6T40xe3X2eCJJdmhpKIkKf7zrww8jRu/cfsghTB09Gp59Fn73O6ishCuugCuuYPpLL1FWV8fqrVtpsnsh+AMB3C4XmV4vqW43O+XnRwwH7Sh0FOCbLVsor69vNX8Nychg5whdzhRF6T8kqlDd4TGO7+dgrj5DIsohFJeUcM2iRQSADK+X5kCAtdXVjB48mCdfe42pr74KCxbAlCnw5puwxx7Aj2Gi9T4fqSkp+I1BRPAbE7OPckchprcfcgizFy5kp/z8tjueJJekVhSl9xJTGYjIdsCZwHAR+QqrxWW7RjbGmMYkyNcvuHfJElqMIT3F+rpTXC4kEOBXb7zBnFdftSKDbr/d6j3gdrdeF/RZuF0uKhsa8AcCGMATRx/lWA7aTK+X5Zs3gwgThw3T+H5FGeB0GFoqImOBL7B6GFwA3I3VyGZYN8jWq+hKg5pgGKXfNsmNKS/nxXnzuPX55/l6zBguvPVWJmdlMf3559vMO2vKFDbX1FDb1ITPVgQAvkCAxpYWhmRkRH2ijxZiesAOO7T6QfYfNYoJw4a1C5NVFGXgESvP4BasRja/wookOgqoAuYkV6zeRVeLso3JzWVoRgb4fJz/1lu8/fe/M37jRs4/5hhOPu88VmRlRZx3amEhabZ5yIXlsQ++/H4/OxcURH2ij5bLsHj9ek32UhSlHbHMRHsCk4wx1fb7r0XkC+A/SZWql9HVBjWzpkzhwX/9i6cefZTd1q/nxfHjufzII0kbNYrhaWkdzruuuprBaWmtYaAALYEAgUAgprM3kr9j9sKFCS3RrShK/yDWzqAhRBEAYIwpIYISEZF9EyhXr6JL2bKNjUydO5cnbr6Zkdu28cdzz+Xxq6/mwZkzSXG5Ys8rQnjElzGmNQPZKfGUq1AUZeARa2cQOQsp8vj9wN5dE6d3EB5OGoy4cZwt++67cPbZ8M03yFlnkfO3v3FryDXxZOFOHDaMZRs3IiKtOQE+v79NIpwTNNlLUZRIxFIGI0TkL7RPMBseYXx4QiXrISIla22qqWn9oHEtoNu2wZ//DPffD2PGwOuvwy9+0e60eBbma6dNY+aCBVQ0NNDo8+F1uxmZnc21nVy8NdlLUZRIdJh0Ftr0Pg6MMcYd+7Suk8yks2jJWgBDMzNjJ5+98gqcdx5s3GiFil5/PdiJXpGIJ6lN6/4ripIIupJ09iWxk83A2iG84lSw3ki0ZK3S2lqKzjgj+oXl5XDxxfDkkzB+PDz3HPz0pzHvF09Sm/YBVhQl2cRSBk8bY9bGM5GIPJ0AeXocx9U0jYGnnrJ2AdXVcO21vHPSSdzz+ees+fzzTj/J625AUZTupMNoImPM9fFO5OTc3ky8jWcAWL8ejjgCZsyAsWPh008pPu00Llm0qEuN4rXZvKIo3Y2T5jZtEJG5iRSktxCr8UxxSQknPPMMN8+YQf0uu+B/6y24805YvBjGj09IB6++1AWsK5nZiqL0HpwUqgsnohOiPxDNRl9cUsLdjz3GnCeeYOI337Bk/HiuP/lkZh91FFPtmkLx9iDuiETM0R1oH2NF6T90RRkMLHw+frjiCp4sKqI5NZU7zjuPtw48kG3NzW0yhhPRwasnu4A58VV0NTNbUZTeQ6fNRAyA5jatfPIJ7LMPM556io8nTeKC227jrZ/9DETaPbE78jlEIRFzdAanvgrtY6wo/YeuKIOzEyZFb6WhAS6/HPbZB0pLufXSS7nqvPOoyslpPSX8iT3eZvcdkYg5OoNTX4WWtlCU/kOnzUTGmGUikgLMA6u6sjHmrEQJ1uMUF8M558Dq1VZJidtu46dVVTwbRymHROQF9ERugVNfhZa2UJT+Q1d2BgABYBFQDPxfl6XpDVRXWxnE06aB3291Hps3D3JyeuyJvbtw+qTf378PRRlIxN0DOeZEImuMMWMSMlkMklaOYsECOP982LQJLrkErr0WMjISf59eSmh0UOiTvi7witI/6KgcRVd3BqE40ioiMlREnhCRr+3XcyKyfQLliZ+yMjjpJDjySMjLgw8/hNtuG1CKAPRJX1EGMj0SWioiXuB14BtgPJYieQh4W0T2MsbUdosgxsDjj1s1hWprraJyf/wjhEXIDCS0DpKiDExi9UD+d5LuezowAbjcGNNijPEDlwM/Ac5P0j3bsnYtHH44nHYa7LorfPYZXHXVgFYEiqIMXGLtDCaLyKnEl1OQFfuUVo4D1hljvg8OGGM2i8gq+9htDuZyRiBg9Rn405+s9/fcAxdcAK5EWsw6hxanUxSlp4ilDHYAHiE+ZeDEZzABy0QUzhrgYAfzOOOHHyzfwHvvwaGHwj/+AaNHJ+12TtDSDoqi9CSxlMF3xJdcJsCTDu5bACyLML4NyBCRdGNMQ5sbiMwEZgKMGjXKwa1CyM2FpiZ47DE45ZRO9xFOBlraQVGUniSWMnjDGFMcz0Qi8mYC5ImKMWYuMBes0NJOTZKZCR9/3KuUQJC+UpxOUZT+Sax+BnE7c40xpzq47xYgO8L4IKA+fFeQUHqhIgAt7aAoSs8SK5podrwT2WaceFkOFEYYHwN84WCefkNPFadTFEWB2ElnMxzMdZ6Dc58HRotIYXBARIYB44BkhbP2ajThS1GUnqTDchQiEsBBlJAxxh3XTa2ks6XAl1gKJwD8CzgAiJl0lrRyFIqiKP2YjspRxHIg/wBsB7wGlHV0DxwUqjPGNIvIL4E7gFVYCmcF8PNuyz5WFEVRWomlDH4CnAr8DigBbjPGlEQ6UUQ+dXJjY0wpcLKTaxRFUZTkECuaqMUY8zAwCXgTeMouLrd7hNOTlyymKIqiJJW4ajAYi+eNMfsBDwN3isjLInJgyDkVyRJSURRFSS6Oq5YaY94A3hCRfYBHRWQrcLIxZl3CpVMURVG6hU41txGRXwN/BvYH6oADjDGfJ1i2ju5fDqztwhQFWIlvSnLR77l70O85+fSX73i0MWZIpANxKwMREeAk4E/A7kAFcC9wT18zEYnI0mjhVUri0O+5e9DvOfkMhO84ppnIzgk4E7gMK7poPXAJMM8YU59c8RRFUZTuIFY5ij9ihZTeDzRhKYUdjTF3hSuCJDbCURRFUZJMrJ3BLVjZwS/aLwOcLO2LvQlW+GlfYW5PCzBA0O+5e9DvOfn0++84VjmKrcBLcc51hDGmICFSKYqiKN1KrJ3BOmPMmfFM5DQDWVEURek9xEo6O9HBXId0RRBFUTqHiLwrIia0CrCiOCVWOYpIfYqjnVvedXGSh4gMtUtpfG2/nhOR7Xtarv6EiOwpIvNE5EsR+UJEVonI3SISMa5Z6ToichxWtV8lCYjIcSLyjogsE5HvRWSpiDhp5NVniKscRV/HDo99HfAC44HdsJLl3haRrI6uVRzxNJAHTDbG7AH8EmvH+J6IpPeoZP0Q+9/1LcCrPS1Lf0RE/gBciVVhYRKwC/AN/bQO24BQBsDpwATgcrv4nh+4HCtvIu7WnkpcXG6MqQMwxmwAbgN2Ag7vUan6JxcCS+yXkkBsk9stwLnGmB8AjDE+4FKsZNt+x0BRBsdhOcO/Dw4YYzZj9VI4rsek6n9MMMZ8Gza20f6pzZwTiIjkYSWC/rmnZemnnApUGWPaKFpjzEZjTL/srDVQlMEEYE2E8TXAHt0sS7/FGNMcYXhnrPyUd7pZnP7OHOBxY0xXanQp0dkfKLF9Bu+KyFci8r6InNXTgiULx1VL+ygFwLII49uADBFJN8Y0dLNM/R4RcQO/Bf7lJBhB6RgR2Qk4AatnuJIcdgAKscxCx2B1ejwOq6fLCGPMjT0oW1IYKDsDpWf4C+ADLu5hOfobfwVuMcZU97Qg/Zg0IBO4zBiz2RgTMMYUYVViuEJEMnpWvMQzUJTBFiA7wvggoF53BYlHRM7Eeno9LOhQVrqO3VBqd+CBnpaln1Nj//wsbPxTIAMrIrFfMVDMRMuBXSOMjwG+6GZZ+j12HPZs4OfGmLKelqef8UvADSwJqRE23P75qog0A1cYYzTctGt8BexJ+wdmv/2z3z1I97sPFIXngdGhGZoiMgzL5qrVVhOIiJyCFbb7CztiCxH5PxGZ2bOS9Q+MMXOMMTsaY/YMvoB/2IcPt8dUEXSdBfbPCWHjuwMNwMruFSf5DBRl8AjWDuCvIpIiIi6sGOI16HY7YYjIDGAe1vf9CxE5xVYORwAje1I2RXHIM1j5GzcEE1NtE93xwI390fTZqbaXfRF7J3AHMBkr1HEFcLExZn2PCtaPEJEKoucTXGuMuaYbxen3iMjhwE1YZqJhwJdAs71bULqIncvxV6ws+kasni73GGPm9ahgSWLAKANFURQlOgPFTKQoiqJ0gCoDRVEURZWBoiiKospAURRFQZWBoiiKgioDRVEUBVUGiqIoCqoMFCVhiMgeIvKZiDSLyCM9LU9nEJEhInJFku9xjYjkJPMeinNUGShREZGh9uJWISLG/v1zu9HHuyIyO7yHtIgcHrIg1tq/h79qRWRayJyfRbhP6NgZUeQ72z7H2Od9JiLLReQHEVkhIrO65YuyMcZ8YWf/box1bm9ERIYD7wGrHVxzi4hssP8Ga0TkyrDjI+2/S6OIrLRLOqzF6oudn9hPoHQJY4y+9NXhC6vWkAkb2xt4F/gO2DXCNSXAoijzLQKmRZgz0n0eAc6IIZ8BHgl578Yq02CA03rg+yoJlaevvICFwH2duO7/2d/17VGOD8WqA+YKGXsYeLGnP7O+fnzpzkDpFMaYT4CDgSqs0slOmn1chVUb6n9xnPsFsMGhbH7gNvvtEU6uHaiIyCSs8ti3O73WGPMe8C0wQ0QilcWfATxtjAmEjN0GHCkiu3dGXiXxqDJQOo2xeh7PweoL8dtY54vIGSLyiDFmsTFmizHm0Djucbsx5vVOiBdclFoXIBEpEJF7bLPFp7ZJ6SYRSQs55zwRWWWbPS4Ukbm2aawkktnJNoutFJG1IrJYRA6LV0ARGS4i/7ZNXMtE5H7b7BI0lR1pn3eciLwpIkttmReJyL4h8wRNMbX2sYtF5GMR2SQiH4nIlDjEOR5Yb4z5PoKcx4jIJyLyjW0Kul9EBoWd9hhWsbxIf9PTsXZ4rRhjVgGbsRogKb2Bnt6a6Kv3v4hgvgk5lg60AK+GjZcQZiYCzqAD80lH94khX7iZKBOYC9QDB4aM74tVhz7Xfj8Yy9R1d9h8hfacy4Ex9thMLMWyS8h5e2C19bzWfu/GKuFd09HnDLn+A6xOWkF5DrOvNWHnvQbMDHl/HFAL7BB23iKgDrjBfu/CMsdUAQUxZPkQeCvC+In25z7Bfj8IeB94E7vQZch3FgCeC7t+IvBhlHsuAt7r6X/f+rL/Hj0tgL56/yvWIo31hLcybKzEXrA+C3mtS6IyqLDvsdxWAhuAE8POy4ywgJ5nL6DhC5sB7gwZK7DHQhflJ+2FNjVkbDvClFMUmQ+xzwuX8fUIymAnwB3hO788bGwRUA2khYyNwFLWN8SQpwR4NmxM7L/Z4rDxQ23Zp0W4fxOQFzL2d+D8KPcsAr7v6X/f+rJeaiZSEoFEGV9q2nbkmpNEGV6y7zMByMLahcwTkbuCJxirIclU25SyQkQ+A67A6mk7PMKc34T8XmH/HBYyth+WEmwKuccGLAURi/3sn5+Eja+IcG4j8IBt2vrcljsP+EmEc78zxjSGyLMJK7ppvwjnhjLUvk8ouwA7YEUYRZJxWtj4o4AX+A2A7T84Dng6yj0baft9Kj2IKgOlS9iO43yscMEOMcY8Yow5I9kyGWMCxvIzPAhcJCJ7gOUPwLJt322M2T1MQaVGmKo+dE77V3fI8eFEXvir4xAzqHzCr29zrYhkAm8DO2P1lJ5ofgxfjSTztghjlcTuNOenvVIvsH+eEhruC7wClGIp0VCew/rOTrffHwp8bIypjHJPwdq1KL2ASJ5/RXFCsEF7b+y7+539czesqKTTgBXGmOcTNP8mInd2y4nzWuzryzu49v8BOwJ/7GBRDSXcsQvWLuKbCOOhlGL5f0LZYv980BhzXawbG2NqROQ/WFFFu2EphUc7uCTDvq/SC9CdgdJp7Cic67BiyP/Vw+JEYgf7Z3DBScWydYcSyTwULx8Au4lI6xO6iGyP5ZiO51qw8jVCGR/2Pjh3q9xi9fAeEmXeHcOio0Zg7Qo+iHJ+kB9ob7L5GstnMDH8ZBH5q4gcFGGe4OJ/CfBTLOd3NIYD2na2l6DKQOkUIjIZeAvLPn+4Maahh0Vqg20aOgfLobzYHn4F2ENEfmmfM9w+p7PciPV0e4U9n9sea+roIgDbjPUe8CcRybWvPwwrQimU97FMSReELPKzaW+iCeIDrrTnc9ny1AB3xhDpv1g+glAZjX2vI0Tk/4LjInIScBLt/R1gRRltwAo1fs4YE9EMJCIC7Ir1N1F6Az3twdZX731hORU/w3KeGn6MCvoKa4G9FMgOu+Zw+5xmfowmuizGffIj3OfMOOQ72z43NJroMyyTyNfAXYSEVAJpWIviBmAZsMB+b4BVWLH2J9q/G6yn4iuxzEzB+2wGHguZ81Ash+o64GOsRbIkKE8M+Ydh2dkrgKVYiVg30D6a6AB77g1Y/oM5WE/yFcAbIectsl9nYoWKbrKvmxLHd7kjVmjopAjHjgCWAN9jKYAiYMcO5rrZ/q4mdnDOAVh+ih1iyaav7nmJ/YdRFKUXICJ3Y4WbOo6yEZFFAMaYaZ289zwg3xhzbGeud3iv14BvjDEXJfteSnyomUhReggR+Veofd9md6zEuJ7gImBQeLG5RCMif8WKIro0mfdRnKHKQFF6jknAxcE3InIUVuz+jT0hjLH8PodimcmSyRLgSGOVM1F6CWomUpQeQkTOx0qOy8QK66wBbjbGPONwnpFYob1j7aFvgYONMVsTJ63S31FloCiKoqiZSFEURVFloCiKoqDKQFEURUGVgaIoioIqA0VRFAVVBoqiKArw/wGzhf4ZZO0YbgAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Plot the predictions\n",
    "\n",
    "plt.rcParams['font.size'] = 16\n",
    "plt.rcParams[\"font.family\"] = \"serif\"\n",
    "\n",
    "plt.scatter(ytest, ytest_pred, \n",
    "            color='teal', \n",
    "            alpha=0.7)\n",
    "plt.xlabel(\"DFT Band gap (eV)\")\n",
    "plt.ylabel(\"ML-Predicted Band gap (eV)\")\n",
    "xmax = max([max(ytest),max(ytest_pred)])\n",
    "plt.plot([0,xmax],[0,xmax],color='red')\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Close the plots\n",
    "plt.close()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### That's it for now, folks!\n",
    "I hope this notebook helped to get started on retrieving OQMD data via OPTIMADE API and using it to build a quick ML model.\n",
    "\n",
    "I will try to add more descriptive information and enhancements to this tutorial in the future.\n",
    "\n",
    "Let me know if you have any questions!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python [conda env:optimade]",
   "language": "python",
   "name": "conda-env-optimade-py"
  },
  "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.7.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}