ml_code/svm_decision_boundary.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import pandas as pd\n",
    "import requests\n",
    "import json \n",
    "import numpy as np\n",
    "import matplotlib.pyplot as pl"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import os\n",
    "import pickle\n",
    "import geopandas as gpd\n",
    "from shapely.geometry import Polygon, Point"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn.svm import SVC\n",
    "from sklearn.impute import SimpleImputer\n",
    "from sklearn.preprocessing import StandardScaler\n",
    "from sklearn.model_selection import train_test_split, GridSearchCV\n",
    "from sklearn.metrics import accuracy_score, confusion_matrix, classification_report\n",
    "from sklearn.decomposition import PCA\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "threshold = 95\n",
    "data_location = '../../new_data/threshold_{}percentile/ca_data_version_6.csv'.format(threshold)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "solar = pd.read_csv(data_location)\n",
    "# solar.drop(columns='Unnamed: 0', index=1, inplace = True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "solar"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "solar[solar['Solar_farm_present']==1]['land_cover'].values"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "solar.columns\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "solar.describe()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "solar.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "solar = solar[solar['min_distance'].notnull()]\n",
    "solar.describe()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "solar = solar[solar['Solar_farm_present'].notnull()]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "solar.describe()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "solar_existing = solar[solar['Solar_farm_present']==1]\n",
    "solar_non_existing = solar[solar['Solar_farm_present']==0]\n",
    "solar_non_existing_sample = solar_non_existing.sample(500, random_state = 42)\n",
    "solar = pd.concat([solar_existing, solar_non_existing_sample])\n",
    "solar.describe()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "req_col = ['land_cover','slope', 'elevation',  'solrad_annual','ac_annual', 'min_distance']\n",
    "X = solar[req_col]\n",
    "Y = solar['Solar_farm_present']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "X_train, X_test, y_train, y_test = train_test_split(X, Y, test_size=0.2, random_state=42)\n",
    "\n",
    "sd = StandardScaler()\n",
    "columns_ML_excluding_LC = ['slope', 'elevation',  'solrad_annual','ac_annual', 'min_distance']\n",
    "X_train_excluding_LC = X_train[columns_ML_excluding_LC ]\n",
    "temp = np.hstack((sd.fit_transform(X_train_excluding_LC),X_train['land_cover'].values.reshape((-1,1))))\n",
    "X_train_2 = pd.DataFrame(temp, columns = columns_ML_excluding_LC + ['land_cover'])\n",
    "\n",
    "X_test_excluding_LC = X_test[columns_ML_excluding_LC ]\n",
    "temp = np.hstack((sd.transform(X_test_excluding_LC),X_test['land_cover'].values.reshape((-1,1))))\n",
    "X_test_2 = pd.DataFrame(temp, columns = columns_ML_excluding_LC+['land_cover'])\n",
    "\n",
    "\n",
    "assert np.sum(np.abs(X_test_2['land_cover'].values - X_test['land_cover'].values)) == 0\n",
    "assert np.sum(np.abs(X_train_2['land_cover'].values - X_train['land_cover'].values)) == 0\n",
    "\n",
    "X_test = X_test_2\n",
    "X_train = X_train_2\n",
    "\n",
    "col_index = {'slope':0, 'elevation':1,  'solrad_annual':2,'ac_annual':3, 'min_distance':4,'land_cover':5}\n",
    "index_col = {0:'slope', 1:'elevation',  2:'solrad_annual',3:'ac_annual', 4:'min_distance',5:'land_cover'}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "len(X_train),len(X_test)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "svm = SVC(kernel='rbf', gamma=0.1, C=1.0, random_state= 113)\n",
    "\n",
    "params = {\"C\":(0.1, 1, 2,4),  \n",
    "          \"kernel\":('linear', 'rbf','poly')}\n",
    "\n",
    "svm_cv = GridSearchCV(svm, params, scoring=\"accuracy\", n_jobs=-1, verbose=2, cv=5)\n",
    "svm_cv.fit(X_train, y_train)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "svm_cv.best_params_"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# svm_tuned = SVC(**svm_cv.best_params_, random_state=113)\n",
    "svm_tuned = SVC(C=4,kernel='rbf',random_state=113)\n",
    "pca = PCA(n_components=6)\n",
    "\n",
    "Xreduced = pca.fit_transform(X_train)\n",
    "# Xreduced = np.vstack([X_train['ac_annual'].values, X_train['solrad_annual'].values]).T\n",
    "print(Xreduced.shape,y_train.shape)\n",
    "\n",
    "svm_tuned.fit(Xreduced, y_train)\n",
    "y_pred_train_svm = svm_tuned.predict(Xreduced) \n",
    "Xreducedtest = pca.transform(X_test)\n",
    "# Xreducedtest = np.vstack([X_test['ac_annual'].values, X_test['solrad_annual'].values]).T\n",
    "\n",
    "y_pred_test_svm = svm_tuned.predict(Xreducedtest)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 975,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[0.1894152  0.46534452 0.56436825 0.58545286 0.25867676 0.18111   ] [0.1894152  0.46534452 0.56436825 0.58545286 0.25867676 0.13926993] 2.3673203324495815\n",
      "[ 0.42256854  0.42453365 -0.39753816 -0.34966757  0.60014297 -0.04040139] [ 0.42256854  0.42453365 -0.39753816 -0.34966757  0.60014297 -0.02704499] 1.4400893612934966\n",
      "[ 0.86410396 -0.14498632  0.0459842  -0.00749583 -0.47788059  0.07652576] [ 0.86410396 -0.14498632  0.0459842  -0.00749583 -0.47788059  0.04198465] 0.864462777312106\n",
      "[ 0.19168326 -0.75436258  0.18226372  0.11844854  0.58551338 -0.18772732] [ 0.19168326 -0.75436258  0.18226372  0.11844854  0.58551338 -0.06417626] 0.3694196297095739\n",
      "[ 0.04254719  0.08695488 -0.0214415   0.1878255  -0.04023582 -9.00768597] [ 0.04254719  0.08695488 -0.0214415   0.1878255  -0.04023582 -0.97635624] 0.0864595404076954\n",
      "[-0.01814709  0.0747032   0.69832481 -0.69685929  0.01000275 -9.14735696] [-0.01814709  0.0747032   0.69832481 -0.69685929  0.01000275 -0.1439432 ] 0.01137040738091612\n",
      "\n",
      " Varaiance covered as we increasing the number of PCA components:  [0.46064684 0.74086773 0.90907988 0.98096368 0.99778748 1.        ]\n",
      "0.1948368443326801 [ 0.1894152   0.42256854  0.86410396  0.19168326  0.04254719 -0.01814709] [0.46064684 0.28022089 0.16821215 0.0718838  0.0168238  0.00221252]\n",
      "[[0.19483684]\n",
      " [0.19483684]\n",
      " [0.19483684]\n",
      " [0.19483684]\n",
      " [0.19483684]\n",
      " [0.02581578]] 1.0000000000000002\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "[]"
      ]
     },
     "execution_count": 975,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAASkAAAGKCAYAAABQEbfCAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8qNh9FAAAACXBIWXMAAA9hAAAPYQGoP6dpAABHY0lEQVR4nO3deVxV1Rrw8d8BZVAGURTREFAcUBAUyhDHrkJSmNlg6sXMoet0FQ1TQxwoh0oNtdCcsrqa+Fr62o1UMjEUR5QyxTERrh4uogWKCQLn/cPLed0CyuGA7I7P937258PZZ+291z7dnp619tpraXQ6nQ4hhFAps9qugBBCPIgEKSGEqkmQEkKomgQpIYSqSZASQqiaBCkhhKpJkBJCqJoEKSGEqtWp7QqIv5aSkhKuXLmCra0tGo2mtqvzyOh0Om7cuEGzZs0wM6u5/7bfvn2bwsJCo89jYWGBlZVVNdSo9kmQEga5cuUKLi4utV2NWpOZmckTTzxRI+e+ffs2jazrcQvjXwJp2rQpFy9eNIlAJUFKGMTW1haAodTDgscnkypExwZu6e+/Rq5RWMgtdIRR36jfthAdX2ZlUVhYKEFKPH5Km3gWaB6rIFXqUTRxrTDDwojrmJnY67jScS6EUDXJpIRQGTOMyx5MLfOQICWEymg0YGZEq1IDVEPfu2qYWtAVQpgYyaSEUBlp7ilJkBJCZcw0GsyMeboH0twTQohHRTIpIVRGmntKEqSEUBkzI5/umVqQMrX7EUKYGMmkhFAZae4pSZASQmU0Go1R7wia2huVEqSEUBnJpJRM7X6EECZGMikhVEae7ilJkBJCZTQYF2hMrU/K1IKuEMLESCYlhMpUy7t7JkSClBAqI0/3lEztfoQQJkYyKSFURp7uKUmQEkJlpLmnZGr3Y1KGDx/OgAEDarsaQtQqyaSEUBkzNJgZMdrJ1DIPCVJCqIz0SSmZ2v38JW3ZsgVvb2+sra1p1KgRffr0IT8/v0y5goICJk6cSJMmTbCysqJbt24cOXJE/31iYiIajYbvvvsOHx8frKys6NKlCydOnFCcJzk5mR49emBtbY2LiwsTJ04s93qidphVw1YVsbGxuLu7Y2VlhZ+fH0lJSRWW1Wq1DBkyhLZt22JmZkZ4eHi55WJiYmjbtq3+/2uTJ0/m9u3bBtVLglQt02q1DB48mBEjRpCWlkZiYiIDBw5EV85S2W+//TZff/01n3/+OceOHcPDw4Pg4GCuX7+uKDd16lQWLVrEkSNHaNKkCf379+fOnTsAnDhxguDgYAYOHMgvv/xCXFwc+/btY8KECeXWr6CggLy8PMUmTE9cXBzh4eFERkZy/PhxunfvTr9+/cjIyCi3fEFBAY0bNyYyMhIfH59yy2zYsIHp06cze/Zs0tLSWLt2LXFxccyYMcOguml05f3bIB6ZY8eO4efnR3p6Oq6urorvhg8fzh9//MG2bdvIz8/HwcGB9evXM2TIEADu3LmDm5sb4eHhTJ06lcTERHr37s2mTZsYNGgQANevX+eJJ55g/fr1vPrqqwwbNgxra2s+/fRT/XX27dtHz549yc/Px8rKSlGHOXPmMHfu3DL1foP6WJjcW2IVK0THZ+STm5uLnZ1djVwjLy8Pe3t7Flo7YGXEiPPbOh3T//zdoLp26dKFzp07s2LFCv0+T09PBgwYwIIFCx54bK9evfD19SUmJkaxf8KECaSlpbF79279vrfeeovDhw8/MEu7n2RStczHx4e//e1veHt788orr7B69Wp+//33MuUuXLjAnTt3CAwM1O+rW7cuTz31FGlpaYqyAQEB+r8bNmxI27Zt9WVSUlJYv349NjY2+i04OJiSkhIuXrxY5rozZswgNzdXv2VmZlbXrYsK3H3BWFPlrTS83Z8BFxQUlHu9wsJCUlJSCAoKUuwPCgoiOTm5yvfRrVs3UlJSOHz4MAC//fYb8fHxPPfccwadRzrOa5m5uTkJCQkkJyeza9culi9fTmRkJIcOHVKUK01475+xUafTVWoWx9IyJSUl/OMf/2DixIllyrRo0aLMPktLSywtLSt9P0I9XFxcFJ9nz57NnDlzypTLycmhuLgYJycnxX4nJyeysrKqfP3XXnuNq1ev0q1bN3Q6HUVFRYwdO5bp06cbdB4JUiqg0WgIDAwkMDCQWbNm4erqytatWxVlPDw8sLCwYN++fYrm3tGjR8t0Wh48eFAfcH7//XfOnj1Lu3btAOjcuTMnT57Ew8Oj5m9MVEl1Pd3LzMxUNPce9h+bqv4HsCKJiYnMmzeP2NhYunTpwvnz55k0aRLOzs5ERUVV+jwSpGrZoUOH2L17N0FBQTRp0oRDhw5x9epVPD09+eWXX/Tl6tevz9ixY5k6dSoNGzakRYsWfPDBB9y6dYuRI0cqzhkdHU2jRo1wcnIiMjISR0dH/aDQadOm8fTTTzN+/HhGjx5N/fr1SUtLIyEhgeXLlz/KWxcVqK4R53Z2dpXqk3J0dMTc3LxM1pSdnV0muzJEVFQUYWFhjBo1CgBvb2/y8/N58803iYyMxMyscncpfVK1zM7Ojp9++omQkBDatGnDzJkzWbx4Mf369StTduHChbz00kuEhYXRuXNnzp8/z86dO3FwcChTbtKkSfj5+aHVatm+fTsWFhYAdOzYkb1793Lu3Dm6d+9Op06diIqKwtnZ+ZHcr1AfCwsL/Pz8SEhIUOxPSEiga9euVT7vrVu3ygQic3NzdDpduU+vKyKZVC3z9PRkx44d5X63fv16xWcrKyuWLVvGsmXLHnjObt268euvv1b4/ZNPPsmuXbsMrqt4NGpjMOeUKVMICwvD39+fgIAAVq1aRUZGBmPGjAHuPkC5fPkyX3zxhf6Y1NRUAG7evMnVq1dJTU3FwsKC9u3bAxAaGsqSJUvo1KmTvrkXFRVF//79MTc3r3TdJEgJoTLGvxZj+LGDBg3i2rVrREdHo9Vq8fLyIj4+Xj8sRqvVlhkz1alTJ/3fKSkpbNy4EVdXV9LT0wGYOXMmGo2GmTNncvnyZRo3bkxoaCjz5s0zqG4yTsqElI6T+v3332nQoEGNXKN0LI+Mk6p+pb/t0voNsdZUvSfmT10Jk/Kv12hdHyXJpExIr169DGrrC3WSd/eUJEgJoTIajFvxxdTyWwlSQqiMZFJKpnY/QggTI5mUECpTG0/31EyClBAqI809JVO7HyGEiZFMSgiVuTtVi3HHmxIJUkKojAxBUJLmnhBC1SSTEkJlzDQazIyYx0me7gkhapQ095SkuSeEUDXJpIRQGcmklCRICaEyEqSUJEgJoTIajcaoBRA0JhampE9KCKFqkkkJoTLS3FOSICWEylTXklamwtTuRwhhYiSTEkJlNJq7W5WPr76qqIIEKSFURvO//xlzvCmR5p4QQtUkkxJCZeTpnpIEKSFURoKUkjT3hBCqJpmUECpjhpELMZjYItYSpIRQGXm6pyRBSggVMq0wYxzpkxJCqJpkUkKojNEjzk0sDZMgJYTKyBAEJWnuCSFUTYKUECpjhsborSpiY2Nxd3fHysoKPz8/kpKSKiyr1WoZMmQIbdu2xczMjPDw8HLL/fHHH4wfPx5nZ2esrKzw9PQkPj7eoHpJkBJCZTTVsBkqLi6O8PBwIiMjOX78ON27d6dfv35kZGSUW76goIDGjRsTGRmJj49PuWUKCwvp27cv6enpbNmyhTNnzrB69WqaN29uUN2kT0oIwZIlSxg5ciSjRo0CICYmhp07d7JixQoWLFhQprybmxtLly4FYN26deWec926dVy/fp3k5GTq1q0LgKurq8F1k0xKCJUpfbpnzAaQl5en2AoKCsq9XmFhISkpKQQFBSn2BwUFkZycXOX72L59OwEBAYwfPx4nJye8vLyYP38+xcXFBp1HgpQQKlNdzT0XFxfs7e31W3kZEUBOTg7FxcU4OTkp9js5OZGVlVXl+/jtt9/YsmULxcXFxMfHM3PmTBYvXsy8efMMOo8094QwUZmZmdjZ2ek/W1paPrD8/cto6XQ6o5bWKikpoUmTJqxatQpzc3P8/Py4cuUKH374IbNmzar0eSRICaEy1fXunp2dnSJIVcTR0RFzc/MyWVN2dnaZ7MoQzs7O1K1bF3Nzc/0+T09PsrKyKCwsxMLColLnkeaeECpjpjF+M4SFhQV+fn4kJCQo9ickJNC1a9cq30dgYCDnz5+npKREv+/s2bM4OztXOkCBBCkhVKc2hiBMmTKFNWvWsG7dOtLS0pg8eTIZGRmMGTMGgBkzZjBs2DDFMampqaSmpnLz5k2uXr1Kamoqp06d0n8/duxYrl27xqRJkzh79izfffcd8+fPZ/z48QbVTZp7QggGDRrEtWvXiI6ORqvV4uXlRXx8vH7IgFarLTNmqlOnTvq/U1JS2LhxI66urqSnpwN3O+537drF5MmT6dixI82bN2fSpElMmzbNoLppdDqdiU2RJWpSXl4e9vb2vEF9LEzuLbGKFaLjM/LJzc2tVD9PVZT+tjscnalvVvVGTn5JCc/maGu0ro+SZFJCqIxMeqckfVJCCFWTTEoIlZH5pJQkSAmhMmYY18QxteaRqd1PjUtMTESj0fDHH3/UdlWYM2cOvr6+tV0NIWqUBKm/CI1Gw7Zt2xT7IiIi2L17d+1USNSY2hgnpWbS3PsLs7GxwcbGprarIaqbRmPUO3Om1iklmVQ5dDodH3zwAS1btsTa2hofHx+2bNlSYfnk5GR69OiBtbU1Li4uTJw4kfz8fODuSN2nn366zDEdO3Zk9uzZABw5coS+ffvi6OiIvb09PXv25NixY/qybm5uALz44otoNBr95/ubeyUlJURHR/PEE09gaWmJr68vO3bs0H+fnp6ORqPhm2++oXfv3tSrVw8fHx8OHDhQ1Z9K1ADJpJQkSJVj5syZfPbZZ6xYsYKTJ08yefJk/v73v7N3794yZU+cOEFwcDADBw7kl19+IS4ujn379jFhwgQAhg4dyqFDh7hw4YL+mJMnT3LixAmGDh0KwI0bN3j99ddJSkri4MGDtG7dmpCQEG7cuAHcDWIAn332GVqtVv/5fkuXLmXx4sUsWrSIX375heDgYPr378+5c+cU5SIjI4mIiCA1NZU2bdowePBgioqKyj1nQUFBmXmJhHiUZMT5ffLz83F0dOTHH38kICBAv3/UqFHcunWLN998k969e/P777/ToEEDhg0bhrW1NZ9++qm+7L59++jZsyf5+flYWVnh4+PDyy+/TFRUFADvvPMOP/zwA4cPHy63DsXFxTg4OLBx40aef/554G6f1NatWxkwYIC+3Jw5c9i2bRupqakANG/enPHjx/POO+/oyzz11FM8+eSTfPLJJ6Snp+Pu7s6aNWsYOXIkAKdOnaJDhw6kpaXRrl27MnWZM2cOc+fOLbNfRpxXv9IR53uaNMfGiBHnN0tK6J192WRGnEsmdZ9Tp05x+/Zt+vbtq+/zsbGx4YsvvlBkQ6VSUlJYv369omxwcDAlJSVcvHgRuJtNbdiwAbjblPzqq6/0WRTcnRJjzJgxtGnTRj9B2c2bNyucX7o8eXl5XLlyhcDAQMX+wMBA0tLSFPs6duyo/9vZ2Vlfh/LMmDGD3Nxc/ZaZmVnpOomq0fyvT8qYzZRIx/l9SqeV+O6778pMGG9paVkmUJWUlPCPf/yDiRMnljlXixYtABgyZAjTp0/n2LFj/Pnnn2RmZvLaa6/pyw0fPpyrV68SExODq6srlpaWBAQEUFhYaHD9KzNxWel80/eWv3c6jXtZWlo+dLI0IWqSBKn7tG/fHktLSzIyMujZs2eZ7+8PUp07d+bkyZN4eHhUeM4nnniCHj16sGHDBv7880/69OmjmEwsKSmJ2NhYQkJCgLszKubk5CjOUbdu3QfODW1nZ0ezZs3Yt28fPXr00O9PTk7mqaeeevBNC1WpypxQ9x9vSiRI3cfW1paIiAgmT55MSUkJ3bp1Iy8vj+TkZGxsbMqsdjFt2jSefvppxo8fz+jRo6lfvz5paWkkJCSwfPlyfbmhQ4cyZ84cCgsL+eijjxTn8PDw4Msvv8Tf35+8vDymTp2KtbW1ooybmxu7d+8mMDAQS0tLHBwcytR96tSpzJ49m1atWuHr68tnn31Gamqqvqkp/ho0Zho0RkQaecH4MfDuu+8ya9YsFixYgKenJ8HBwXz77be4u7uXKduxY0f27t3LuXPn6N69O506dSIqKkrf11PqlVde4dq1a9y6dUvR+Q13l/75/fff6dSpE2FhYUycOJEmTZooyixevJiEhARcXFwU8/jca+LEibz11lu89dZbeHt7s2PHDrZv307r1q2N+0GEqEXydE8YROaTqvmne/uauxj9dK/b5UyTebonzT0hVEZmQVCS5p4QQtUkkxJCZYwd6yTjpIQQNUqae0oSpIRQGcmklKRPSgihapJJCaEy0txTkiAlhMqYaTSYGRFpjDlWjaS5J4RQNcmkhFAZae4pSZASQmU0GPl0z8ReV5LmnhBC1SSTEkJlNGZ3tyofb2JTBkiQEkJtZEkrBWnuCSFUTTIpIVRGnu4pSSYlhMrcDVLGrBZTtevGxsbi7u6OlZUVfn5+JCUlVVhWq9UyZMgQ2rZti5mZGeHh4Q8896ZNm9BoNGVmpa0MCVJCqExpJmXMZqi4uDjCw8OJjIzk+PHjdO/enX79+lW4rFpBQQGNGzcmMjISHx+fB5770qVLRERE0L17d8MrhgQpIQSwZMkSRo4cyahRo/D09CQmJgYXFxdWrFhRbnk3NzeWLl3KsGHDsLe3r/C8xcXFDB06lLlz59KyZcsq1U2ClBAqU/runjEb3J0z/d6toKCg3OsVFhaSkpJCUFCQYn9QUBDJyclG3Ut0dDSNGzfWr5hdFRKkhFCZ6mruubi46FfEtre3Z8GCBeVeLycnh+LiYsVakABOTk5kZWVV+T7279/P2rVrWb16dZXPAfJ0TwiTlZmZqVgt5mErUVdm9evKunHjBn//+99ZvXo1jo6OVTpHKQlSQqhMdc3MaWdnV6klrRwdHTE3Ny+TNWVnZ5fJrirrwoULpKenExoaqt9XUlICQJ06dThz5gytWrWq1LmkuSeEyjzqp3sWFhb4+fmRkJCg2J+QkEDXrl2rdA/t2rXjxIkTpKam6rf+/fvTu3dvUlNTcXFxqfS5JJMSQjBlyhTCwsLw9/cnICCAVatWkZGRwZgxYwCYMWMGly9f5osvvtAfk5qaCsDNmze5evUqqampWFhY0L59e6ysrPDy8lJco0GDBgBl9j+MBCkhVKY2RpwPGjSIa9euER0djVarxcvLi/j4eFxdXYG7gzfvHzPVqVMn/d8pKSls3LgRV1dX0tPTq175csgy68Igssx6zS+zftqnNbbm5lU+z43iYtr9fM5kllmXPikhhKpJc08IlZEXjJUkSAmhMrJajJIEKSFURjIpJemTEkKommRSQqhMdY04NxUSpIRQGQ1GNveqrSbqIM09IYSqSSYlhMpIc09JgpQQamPk0z1Ta+9Jc08IoWqSSQmhMtLcU5IgJYTKGL3Muom1j0zsdoQQpkYyKSFURpp7ShKkhFAbM83dzZjjTYgEKSHURt4wVpA+KSGEqkkmJYTKSJ+UkgQpIdRG+qQUpLknhFA1o4LU5cuX+emnn7h165Z+X0lJCe+//z6BgYH07duXHTt2GF1JIR4rj3p1UJUzqrkXFRXFtm3b+O9//6vfN2/ePGbPnq3/vHfvXpKTk/H39zfmUkI8NjRmGjRGNNmMOVaNjMqkDhw4QJ8+fahbty5wN4tavnw57dq1IyMjg8OHD1OvXj0WLVpULZV9kPT0dDQajX5V1erm5uZGTExMjZy7Ns2ZMwdfX9/aroYQFTIqSGm1Wtzc3PSfjx07Rk5ODv/85z954okn8Pf3Z8CAARw6dMjYegrx+JDmnoJRzb3i4mJKSkr0n5OSktBoNDzzzDP6fc2bNycrK8uYy1QLnU5HcXExderIA02hbhqNkc09EwtSRmVSLVq04PDhw/rP27Ztw9nZmbZt2+r3ZWVl0aBBg0qfc8uWLXh7e2NtbU2jRo3o06cP+fn5lJSUEB0dzRNPPIGlpSW+vr4P7JRPTExEo9Gwc+dO/P39sbS0JCkpiQsXLvDCCy/g5OSEjY0NTz75JD/88IPi2OzsbEJDQ7G2tsbd3Z0NGzZU/kcBlixZgre3N/Xr18fFxYVx48Zx8+ZN/ffr16+nQYMG7Ny5E09PT2xsbHj22WfRarX6MsOHD2fAgAEsWrQIZ2dnGjVqxPjx47lz546+jEajYdu2bYprN2jQgPXr1+s/T5s2jTZt2lCvXj1atmxJVFSU4hxCqJ1RQeqll15i//79vPLKK4SFhbFv3z4GDhyoKPPrr7/SsmXLSp1Pq9UyePBgRowYQVpaGomJiQwcOBCdTsfSpUtZvHgxixYt4pdffiE4OJj+/ftz7ty5B57z7bffZsGCBaSlpdGxY0du3rxJSEgIP/zwA8ePHyc4OJjQ0FAyMjL0xwwfPpz09HR+/PFHtmzZQmxsLNnZ2ZX+XczMzFi2bBm//vorn3/+OT/++CNvv/22osytW7dYtGgRX375JT/99BMZGRlEREQoyuzZs4cLFy6wZ88ePv/8c9avX68IQJVha2vL+vXrOXXqFEuXLmX16tV89NFHlT6+oKCAvLw8xSZqmDT3FIxq+0RERLBr1y6+/vprALy9vZkzZ47++7S0NI4cOcKMGTMqdT6tVktRUREDBw7E1dVVf06ARYsWMW3aNF577TUA3n//ffbs2UNMTAyffPJJheeMjo6mb9+++s+NGjXCx8dH//m9995j69atbN++nQkTJnD27Fm+//57Dh48SJcuXQBYu3Ytnp6elboHgPDwcP3f7u7uvPvuu4wdO5bY2Fj9/jt37rBy5UpatWoFwIQJE4iOjlacx8HBgY8//hhzc3PatWvHc889x+7duxk9enSl6zJz5kz9325ubrz11lvExcWVCZoVWbBgAXPnzq309UQ1MMPIwZzVVhNVMCpI2dnZcfDgQX799VcAPD09MTc3139vbW3N1q1bKz38wMfHh7/97W94e3sTHBxMUFAQL7/8Mubm5ly5coXAwEBF+cDAQH7++ecHnvP+a+fn5zN37lz+/e9/c+XKFYqKivjzzz/1mVRaWhp16tRRHNeuXTuDmqx79uxh/vz5nDp1iry8PIqKirh9+zb5+fnUr18fgHr16ukDFICzs3OZbK1Dhw6K39PZ2ZkTJ05Uuh5wt/kcExPD+fPnuXnzJkVFRdjZ2VX6+BkzZjBlyhT957y8PFxcXAyqgzCMvBajVC0x18vLCy8vL8W/UHD3v9wvvPACzZs3r9R5zM3NSUhI4Pvvv6d9+/YsX76ctm3bcvHiRaDsj6/T6R76D6Q0KJSaOnUqX3/9NfPmzSMpKYnU1FS8vb0pLCzUn7O8a1XWpUuXCAkJwcvLi6+//pqUlBR9pndvX1DpsI1SGo1Gf+0Hlbn3QUV5x9x7jYMHD/Laa6/Rr18//v3vf3P8+HEiIyP191oZlpaW2NnZKTYhHqVqedSVlZXFN998w+nTp7l16xZr1qwB4OrVq1y8eFHfEV4ZGo2GwMBAAgMDmTVrFq6uruzevZtmzZqxb98+evTooS+bnJzMU089ZVBdk5KSGD58OC+++CIAN2/eJD09Xf+9p6cnRUVFHD16VH/uM2fO8Mcff1Tq/EePHqWoqIjFixdjZnb3vwGbN282qI6V1bhxY0Vn+7lz5xSj//fv34+rqyuRkZH6fZcuXaqRuohqJO/uKRgdpGJjY3nrrbcoKCgA7gaZ0iCVnZ1NQEAAK1eurFQ/yqFDh9i9ezdBQUE0adKEQ4cOcfXqVTw9PZk6dSqzZ8+mVatW+Pr68tlnn5GammrwkzcPDw+++eYbQkND0Wg0REVFKbKTtm3b8uyzzzJ69GhWrVpFnTp1CA8Pr3SQbdWqFUVFRSxfvpzQ0FD279/PypUrDapjZT3zzDN8/PHHPP3005SUlDBt2jRF9uXh4UFGRgabNm3iySef5LvvvmPr1q01UhdRjWQ+KQWjmnvffvstEyZMwNvbm+3btzN27FjF9x06dKBjx45lHpNXxM7Ojp9++omQkBDatGnDzJkzWbx4Mf369WPixIm89dZbvPXWW3h7e7Njxw62b99O69atDarzRx99hIODA127diU0NJTg4GA6d+6sKPPZZ5/h4uJCz549GThwIG+++SZNmjSp1Pl9fX1ZsmQJ77//Pl5eXmzYsIEFCxYYVMfKWrx4MS4uLvTo0YMhQ4YQERFBvXr19N+/8MILTJ48mQkTJuDr60tycjJRUVE1Uhfx1xcbG4u7uztWVlb4+fmRlJRUYVmtVsuQIUNo27YtZmZmiodFpVavXk337t1xcHDAwcGBPn36KIYsVZZGd3+nhgF69OhBRkYGJ0+epH79+sydO5fo6GiKi4v1ZYYNG0ZSUpK+X0n8teXl5WFvb88b1MfC1FahfIBCdHxGPrm5uTXWL1f622pDn8KubtUbOXl3inD+9rBBdY2LiyMsLIzY2FgCAwP59NNPWbNmDadOnaJFixZlyqenp/PRRx/h5+fHRx99RM+ePcu8NjZ06FACAwPp2rUrVlZWfPDBB3zzzTecPHmy0v3UYGQmlZqaynPPPVemc/pezZs3V7yALIR4iFoYJ7VkyRJGjhzJqFGj8PT0JCYmBhcXF1asWFFueTc3N5YuXcqwYcOwt7cvt8yGDRsYN24cvr6+tGvXjtWrV1NSUsLu3bsNqptRQaqkpKTME6j7Xb16FUtLS2MuoyobNmzAxsam3K1Dhw61XT0h9O4fhFvab3y/wsJCUlJSCAoKUuwPCgoiOTm52upz69Yt7ty5Q8OGDQ06zqiO87Zt27Jv374Kvy8qKmLv3r36AZmmoH///vpBnvd7WMAWojKqa6qW+8ezzZ49WzHYulROTg7FxcU4OTkp9js5OVXre7fTp0+nefPm9OnTx6DjjApSQ4cOJSIigvfee08xshnuvnwcERHBb7/9xrRp04y5jKrY2tpia2tb29UQpqyanu5lZmYq+qQe1qKpyjjEyvrggw/46quvSExMxMrKyqBjjQpS//znP/n222+ZPXs2X375pf5HePXVVzl69Cjp6ekEBQUxcuRIYy4jhKiCyg6+dXR0xNzcvEzWlJ2dXSa7qopFixYxf/58fvjhBzp27Gjw8Ub1SdWtW5edO3cyffp0cnJy+PXXX9HpdGzZsoXr168zbdo0tm/fbnLD9IWoUaWDOY3ZDGBhYYGfnx8JCQmK/QkJCXTt2tWoW/nwww9599132bFjR5Vn5zV6MKeFhQXz5s3jvffe48yZM1y/fh07O7sy7/EJISqnNt7dmzJlCmFhYfj7+xMQEMCqVavIyMhgzJgxwN13OC9fvswXX3yhP6Z0FtybN29y9epVUlNTsbCwoH379sDdJl5UVBQbN27Ezc1Nn6mVPmiqLKOCVMuWLQkJCeHjjz9Go9HQrl07Y04nhIBaeS1m0KBBXLt2jejoaLRaLV5eXsTHx+tnI9FqtYrpjAA6deqk/zslJYWNGzfi6uqqf80sNjaWwsJCXn75ZcVxFXXgV8SoIJWTkyOdyEKYiHHjxjFu3LhyvytvHrOHjQO/951YYxgVpHx9fTl79my1VEQIUcrYietMqw/YqI7zadOm8e2337Jnz57qqo8Qj73SPiljNlNiVCZ17do1goKC6Nu3Ly+++CJPPvkkTk5O5f5Iw4YNM+ZSQojHlFEvGJuZmZU78dq9Qap0QNi9Lx2Lvy55wbjmXzDOHtwTOwsjXjAuLKLJV3trtK6PklGZ1GeffVZd9RBC/I9MH6xkVJB6/fXXq6seQghRLlkpUwi1kemDFYwKUvcP7nqQ8ibOEkKUQ6YPVjAqSLm5uVWq/avRaCgqKjLmUkKIx5RRQWrYsGHlBqnc3Fx+/vlnLl68SM+ePXFzczPmMkI8VqprPilTYVSQetCS3zqdjsWLF/PBBx+wdu1aYy4jxONFmnsKNbYgs0ajISIigg4dOjB16tSauowQpqd0mfUqb7V9A9Wrxm/H39+fH3/8saYvI4QwUTU+BOHChQvSaS6EAWQwp1KNBKmSkhIuX77M+vXr+b//9//yt7/9rSYuI4RpknFSCkYFqdJ39yqi0+lo0KABH374oTGXEUI8xowKUj169Cg3SJmZmeHg4IC/vz9vvPFGtUzmLsRjQ57uKRgVpBITE6upGkIIPQlSCkY93cvIyCAvL++BZW7cuGHQ6zNCCHEvo4KUu7s7MTExDywTGxuLu7u7MZcR4jGj+f/ZVFU2E5vny6jmXmXmyzNiTj0hHk9mZnc3Y443ITV+N//5z39kRRkhRJUZnElFR0crPlfUeV5cXMx//vMfNm3aRJcuXapUOSEeS9JxrmBwkLp3UT+NRkNiYuIDn/I1a9aM999/vyp1E+LxJEFKweAgVbp8lU6n45lnnmH48OHlTiNsbm5Ow4YNadeuHWYm1kYWokZJkFIwOEj17NlT//fs2bPp3bs3PXr0qNZKCSFEKaOe7s2ePbu66iGEKCVP9xSq7QXjzMxMrly5QkFBQbnfS7YlRCVJc0/B6CD17bffMnXqVM6dO/fAcrI4qBCiKozKCxMTE3nxxRe5efMmEyZMQKfT0aNHD958803at2+PTqfjueeeY9asWdVVXyFMnzGjzY3NwlTIqCC1cOFCbGxsSElJYenSpQD07t2bFStW8MsvvzBv3jx2797NCy+8UC2VFeKxIEFKwaggdeTIEQYMGKCYiqWkpAS4O4ZqxowZdOrUSTIpIUSVGRWkbt26RfPmzfWfLS0ty8yK8PTTT7N//35jLiPE46X06Z4xmwkx6m6aNm3K1atX9Z+bN2/OyZMnFWWuXbsmneZCGKKWmnulM5ZYWVnh5+dHUlJShWW1Wi1Dhgyhbdu2mJmZER4eXm65r7/+mvbt22NpaUn79u3ZunWrwfUyKkj5+Pjw66+/6j/37t2bPXv2sGnTJvLz89m5cydxcXF07NjRmMsIIWpYXFwc4eHhREZGcvz4cbp3706/fv0qnAuuoKCAxo0bExkZiY+PT7llDhw4wKBBgwgLC+Pnn38mLCyMV199lUOHDhlUN43OiLlU1q1bx4QJE0hLS8PV1ZWLFy/i5+dHbm6uvkydOnVISEiQcVImIi8vD3t7e96gPhYmNm/RgxSi4zPyyc3Nxc7OrkauUfrbXpv6KnaWFlU/T0EhjT7cbFBdu3TpQufOnVmxYoV+n6enJwMGDGDBggUPPLZXr174+vqWmVtu0KBB5OXl8f333+v3Pfvsszg4OPDVV19V+n6MyqRGjBjBrVu3cHV1Be5OgnfkyBHGjBlDUFAQo0eP5tChQxKghDBENTX38vLyFFtFA60LCwtJSUkhKChIsT8oKIjk5OQq38aBAwfKnDM4ONjgc1b7klatWrXik08+qe7TCvHY0JiZoTGi87v0WBcXF8X+2bNnK2YxKZWTk0NxcXGZBVOcnJzIysqqcj2ysrKq5ZzVGqSuX79Ofn5+mR9HCPHoZWZmKpp7lpaWDyx//8pPOp3O6IVGq+OcRj+rzM3NZdKkSTg5OdG4cWPFfOaHDh0iJCSElJQUYy8jxGOkeuY4t7OzU2wVBSlHR0fMzc3LZDjZ2dlGLUfXtGnTajmnUUHq+vXrdOnSheXLl+Pi4oKnp6diTvOOHTuyf/9+NmzYYMxlRA3TaDRs27attqshSj3iIQgWFhb4+fmRkJCg2J+QkEDXrl2rfBsBAQFlzrlr1y6Dz2lUkJozZw5nz57lq6++4ujRo7zyyiuK762trenZsyc//vijMZcRQtSwKVOmsGbNGtatW0daWhqTJ08mIyODMWPGADBjxgyGDRumOCY1NZXU1FRu3rzJ1atXSU1N5dSpU/rvJ02axK5du3j//fc5ffo077//Pj/88EOFY6oqYlSf1Pbt23n++ecZNGhQhWVcXV2NekIgxGOnFqZqGTRoENeuXSM6OhqtVouXlxfx8fH6J/darbbMmKlOnTrp/05JSWHjxo24urqSnp4OQNeuXdm0aRMzZ84kKiqKVq1aERcXZ/CaB0ZlUlqtlvbt2z+wjJWVFfn5+cZcpsbs2LGDbt260aBBAxo1asTzzz/PhQsX9N//5z//4bXXXqNhw4bUr18ff3//Sg1Eu3DhAi+88AJOTk7Y2Njw5JNP8sMPPyjKuLm5MX/+fEaMGIGtrS0tWrRg1apV+u/T09PRaDR888039O7dm3r16uHj48OBAwf0ZebMmYOvr6/ivDExMbi5uek/HzlyhL59++Lo6Ii9vT09e/bk2LFjBv5S4pGqpddixo0bR3p6OgUFBaSkpCiGDq1fv77MWgY6na7MVhqgSr388sucPn2awsJC0tLSGDhwoMH1MipINWrUiMzMzAeWOX36NM7OzsZcpsbk5+czZcoUjhw5wu7duzEzM+PFF1+kpKSEmzdv0rNnT65cucL27dv5+eefefvtt/UvUD/IzZs3CQkJ4YcffuD48eMEBwcTGhpa5r9Eixcvxt/fn+PHjzNu3DjGjh3L6dOnFWUiIyOJiIggNTWVNm3aMHjwYIqKiip9jzdu3OD1118nKSmJgwcP0rp1a0JCQrhx40alji8oKCgz3kaIR8mo5l6PHj3Yvn07ly9fVrxoXOrUqVPs2LGDN954w5jL1JiXXnpJ8Xnt2rU0adKEU6dOkZyczNWrVzly5AgNGzYEwMPDo1Ln9fHxUbwq8N5777F161a2b9/OhAkT9PtDQkIYN24cANOmTeOjjz4iMTGRdu3a6ctERETw3HPPATB37lw6dOjA+fPnFWUe5JlnnlF8/vTTT3FwcGDv3r08//zzDz1+wYIFzJ07t1LXEtVEZuZUMCqTioyMpKioiMDAQDZu3EhOTg4AaWlprF27lmeeeQZLS0umTp1aLZWtbhcuXGDIkCG0bNkSOzs7/fCJjIwMUlNT6dSpkz5AGSI/P5+3336b9u3b06BBA2xsbDh9+nSZTOredxo1Gg1NmzYlOzu7wjKlGen9ZR4kOzubMWPG0KZNG+zt7bG3t+fmzZsVvpN1vxkzZpCbm6vfHpY5i2og80kpGJVJeXt7ExcXR1hYGGFhYcDddqqXlxc6nQ5bW1s2b95M69atq6Wy1S00NBQXFxdWr15Ns2bNKCkpwcvLi8LCQqytrat83qlTp7Jz504WLVqEh4cH1tbWvPzyyxQWFirK1a1bV/FZo9GUaU7eW6Z0EFxpGTMzszLL2N+5c0fxefjw4Vy9epWYmBhcXV2xtLQkICCgTF0qYmlp+dBBgELUJKNHnPfv35+LFy/yxRdfcPDgQa5fv46dnR1dunThjTfewNHRsTrqWe2uXbtGWloan376Kd27dwdg3759+u87duzImjVruH79usHZVFJSEsOHD+fFF18E7vZR3d+hWB0aN25MVlaWYhRvampqmbrExsYSEhIC3B2FXJrxCpWS5p6CwUHK3NycOXPmEBUVpd93/vx5zMzM2LRpU7VWriY5ODjQqFEjVq1ahbOzMxkZGUyfPl3//eDBg5k/f77+LXBnZ2eOHz9Os2bNCAgIeOC5PTw8+OabbwgNDUWj0RAVFVWpDndD9erVi6tXr/LBBx/w8ssvs2PHDr7//nvFqxAeHh58+eWX+Pv7k5eXx9SpU43KEsUjIEtaKRh8N6WPGu/1/fffM3ny5Gqr1KNQGlRTUlLw8vJi8uTJfPjhh/rvLSws2LVrF02aNCEkJARvb28WLlyIubn5Q8/90Ucf4eDgQNeuXQkNDSU4OJjOnTtX+z14enoSGxvLJ598go+PD4cPHyYiIkJRZt26dfz+++906tSJsLAwJk6cSJMmTaq9LqIaSZ+UgsHzSZmZmTFnzhzFvOVz584lOjpaZuB8DMh8UjU/n9T16FHYWRkxn9TtQhrOWlOjdX2Uqn2qFiGEkaRPSsG0Gq+PSIcOHbCxsSl3k5ephdFkIQYFyaSqID4+vsyj/lLGTG0hhCirSkHqX//6FwcPHtR/Pn/+PID+Mff9NBoN3333XVUupUqlL10KUSM0GNncq7aaqEKVgtT58+f1geleO3bsKLe8sbP7CfFYkT4pBYOD1MWLF2uiHkIIUS6Dg5Q0dYSoYZJJKUjHuRBqozHyCZ3GtJ7umdbdCCFMjmRSQqiNNPcUJEgJoTYSpBQkSAmhNhoz4/qVpE9KCCEeHcmkhFAbM83dzZjjTYgEKSHURpp7CqZ1N0IIkyOZlBBqI0/3FCRICaE2Mse5gmndjRDC5EgmJYTaSHNPQYKUEGojT/cUTOtuhBAmRzIpIdRGpg9WkCAlhNrI0z0FCVJCqI10nCuYVsgVQlRZbGws7u7uWFlZ4efnR1JS0gPL7927Fz8/P6ysrGjZsiUrV64sUyYmJoa2bdtibW2Ni4sLkydP5vbt2wbVS4KUEGpT+nTPmM1AcXFxhIeHExkZyfHjx+nevTv9+vUjIyOj3PIXL14kJCSE7t27c/z4cd555x0mTpzI119/rS+zYcMGpk+fzuzZs0lLS2Pt2rXExcUxY8YMg+omzT0h1EZj5CwIVWjuLVmyhJEjRzJq1Cjgbga0c+dOVqxYwYIFC8qUX7lyJS1atCAmJgYAT09Pjh49yqJFi3jppZcAOHDgAIGBgQwZMgQANzc3Bg8ezOHDhw2qm2RSQpiovLw8xVZQUFBuucLCQlJSUggKClLsDwoKIjk5udxjDhw4UKZ8cHAwR48e1a/u3a1bN1JSUvRB6bfffiM+Pp7nnnvOoPuQICWE2lRTc8/FxQV7e3v9Vl5GBJCTk0NxcTFOTk6K/U5OTmRlZZV7TFZWVrnli4qKyMnJAeC1117j3XffpVu3btStW5dWrVrRu3dvpk+fbtDPIc09IdSmmp7uZWZmYmdnp99taWn5kMOU19TpdA9cfby88vfuT0xMZN68ecTGxtKlSxfOnz/PpEmTcHZ2JioqqtK3I0FKCBNlZ2enCFIVcXR0xNzcvEzWlJ2dXSZbKtW0adNyy9epU4dGjRoBEBUVRVhYmL6fy9vbm/z8fN58800iIyMxq+R4LmnuCaE2j/jpnoWFBX5+fiQkJCj2JyQk0LVr13KPCQgIKFN+165d+Pv7U7duXQBu3bpVJhCZm5uj0+n0WVdlSCYlhNrUwhznU6ZMISwsDH9/fwICAli1ahUZGRmMGTMGgBkzZnD58mW++OILAMaMGcPHH3/MlClTGD16NAcOHGDt2rV89dVX+nOGhoayZMkSOnXqpG/uRUVF0b9/f8zNzStdNwlSQggGDRrEtWvXiI6ORqvV4uXlRXx8PK6urgBotVrFmCl3d3fi4+OZPHkyn3zyCc2aNWPZsmX64QcAM2fORKPRMHPmTC5fvkzjxo0JDQ1l3rx5BtVNozMk7xKPvby8POzt7XmD+liY2pusD1CIjs/IJzc3t1L9PFVR+tteXz8Pu3pWVT/Prds0HB5Zo3V9lCSTEkJtZD4pBQlSQqiNrLunYFohVwhhciSTEkJtNBojm3umlUlJkBJCbWQ+KQVp7gkhVE0yKSHURp7uKUiQEkJt5OmegmmFXCGEyZFMSgi1keaeggQpIdRGnu4pmFbIFUKYHMmkhFAbWRxUQYKUEKpjZHPPxGanMK2Q+z+JiYloNBr++OOPGrvG8OHDGTBggP5zr169CA8Pr7HricdILay7p2YmmUl17doVrVaLvb39I7vmN998o5829WF69eqFr6+vfs0yIUTFTDJIWVhY0LRp00d6zYYNGz7S6wkTJk/3FP4SeWGvXr345z//SXh4OA4ODjg5ObFq1Sry8/N54403sLW1pVWrVnz//fdA2ebe+vXradCgATt37sTT0xMbGxueffZZtFptpa5fXFzMlClTaNCgAY0aNeLtt98uM5H8/c292NhYWrdujZWVFU5OTrz88svA3Wbi3r17Wbp0KRqNBo1GQ3p6OsXFxYwcORJ3d3esra1p27YtS5cuVVyjtIm5aNEinJ2dadSoEePHj9cvxghQUFDA22+/jYuLC5aWlrRu3Zq1a9fqvz916hQhISHY2Njg5OREWFiYfp00oRKlHefGbCbkL3M3n3/+OY6Ojhw+fJh//vOfjB07lldeeYWuXbty7NgxgoODCQsL49atW+Uef+vWLRYtWsSXX37JTz/9REZGBhEREZW69uLFi1m3bh1r165l3759XL9+na1bt1ZY/ujRo0ycOJHo6GjOnDnDjh076NGjBwBLly4lICCA0aNHo9Vq0Wq1uLi4UFJSwhNPPMHmzZs5deoUs2bN4p133mHz5s2Kc+/Zs4cLFy6wZ88ePv/8c9avX8/69ev13w8bNoxNmzaxbNky0tLSWLlyJTY2NsDdeap79uyJr68vR48eZceOHfz3v//l1VdfrfBeCgoKyqyEK8Sj9Jdp7vn4+DBz5kzg7soVCxcuxNHRkdGjRwMwa9YsVqxYwS+//FLu8Xfu3GHlypW0atUKgAkTJhAdHV2pa8fExDBjxgz9JPMrV65k586dFZbPyMigfv36PP/889ja2uLq6kqnTp0AsLe3x8LCgnr16imapObm5sydO1f/2d3dneTkZDZv3qwIIg4ODnz88ceYm5vTrl07nnvuOXbv3s3o0aM5e/YsmzdvJiEhgT59+gDQsmVL/bErVqygc+fOzJ8/X79v3bp1uLi4cPbsWdq0aVPmXhYsWKCol3gEpLmn8JfJpDp27Kj/29zcnEaNGuHt7a3fV7qIYXZ2drnH16tXTx+gAJydnSsse6/c3Fy0Wi0BAQH6fXXq1MHf37/CY/r27YurqystW7YkLCyMDRs2VJjh3WvlypX4+/vTuHFjbGxsWL16tWKFDoAOHToolgO69z5SU1MxNzenZ8+e5Z4/JSWFPXv2YGNjo9/atWsHwIULF8o9ZsaMGeTm5uq3zMzMh96HMFLppHdV3kwrSP1lMqn7n5xpNBrFvtKlnUtKSip9fE0tlGNra8uxY8dITExk165dzJo1izlz5nDkyBEaNGhQ7jGbN29m8uTJLF68mICAAGxtbfnwww85dOjQQ++j9J6tra0fWK+SkhJCQ0N5//33y3zn7Oxc7jGWlpYPXZ5biJr0l8mkaou9vT3Ozs4cPHhQv6+oqIiUlJQHHlenTh369OnDBx98wC+//EJ6ejo//vgjcPfpY3FxsaJ8UlISXbt2Zdy4cXTq1AkPD48Ks5uKeHt7U1JSwt69e8v9vnPnzpw8eRI3Nzc8PDwUW/369Q26lqhBpc09YzYTIkGqEiZNmsTChQvZunUrp0+fZty4cQ8cKPrvf/+bZcuWkZqayqVLl/jiiy8oKSmhbdu2ALi5uXHo0CHS09PJycmhpKQEDw8Pjh49ys6dOzl79ixRUVEcOXLEoHq6ubnx+uuvM2LECLZt28bFixdJTEzUd76PHz+e69evM3jwYA4fPsxvv/3Grl27GDFiRJmgKWqRDOZUMK27qSFvvfUWw4YNY/jw4fqm2Isvvlhh+QYNGvDNN9/wzDPP4OnpycqVK/nqq6/o0KEDABEREZibm9O+fXsaN26sX8564MCBDBo0iC5dunDt2jXGjRtncF1XrFjByy+/zLhx42jXrh2jR48mPz8fgGbNmrF//36Ki4sJDg7Gy8uLSZMmYW9vj5mJPbYWpkNWMBYGkRWMH8EKxvHrsKtfr+rnyb9Fw5ARsoKxEKKGyKR3ChKkQD/YsTzff/893bt3f4S1EY89GSelIEGKu+OLKtK8efNHVxEhRBkSpAAPD4/aroIQ/5809xQkSAmhMqUvnhtzvCkxrZArhDA5EqSEUJtaGswZGxuLu7s7VlZW+Pn5kZSU9MDye/fuxc/PDysrK1q2bMnKlSvLlPnjjz8YP348zs7OWFlZ4enpSXx8vEH1kuaeEGpTC31ScXFxhIeHExsbS2BgIJ9++in9+vXj1KlTtGjRokz5ixcvEhISwujRo/nXv/7F/v37GTduHI0bN9bPFlJYWEjfvn1p0qQJW7Zs4YknniAzMxNbW1uD6iZBSgjBkiVLGDlyJKNGjQLuTk+0c+dOVqxYwYIFC8qUX7lyJS1atNBPge3p6cnRo0dZtGiRPkitW7eO69evk5ycrH8x3tXV1eC6SXNPCLXRaMDMiO1/Hef3T1ZYUFBQ7uUKCwtJSUkhKChIsT8oKIjk5ORyjzlw4ECZ8sHBwRw9elQ/U+z27dsJCAhg/PjxODk54eXlxfz58w1+T1SClBBqU019Ui4uLtjb2+u38jIigJycHIqLi/VzspVycnIiKyur3GOysrLKLV9UVKSfjvq3335jy5YtFBcXEx8fz8yZM1m8eDHz5s0z6OeQ5p4QJiozM1Px7t7D5gW7f+iCTqd74HCG8srfu7+kpIQmTZqwatUqzM3N8fPz48qVK3z44YfMmjWr0vchQUoItamm12Ls7Owq9YKxo6Mj5ubmZbKm7OzsMtlSqaZNm5Zbvk6dOjRq1Ai4O5Fi3bp1FTPJenp6kpWVRWFhIRYWFpW6HWnuCaE2j3j6YAsLC/z8/EhISFDsT0hIoGvXruUeExAQUKb8rl278Pf313eSBwYGcv78ecVsuWfPnsXZ2bnSAQokSAmhPrUwM+eUKVNYs2YN69atIy0tjcmTJ+vnOYO7c90PGzZMX37MmDFcunSJKVOmkJaWpl9N6d4VmMaOHcu1a9eYNGkSZ8+e5bvvvmP+/PmMHz/eoLpJc08IwaBBg7h27RrR0dFotVq8vLyIj4/XDxnQarWKRUHc3d2Jj49n8uTJfPLJJzRr1oxly5bphx/A3Y77Xbt2MXnyZDp27Ejz5s2ZNGkS06ZNM6huMumdMIhMelfzk979/tM27GyqPud83s18HHoMkEnvhBA1pHS8kzHHmxDpkxJCqJpkUkKojcwnpSBBSgi1kemDFUwr5AohTI5kUkKojTT3FCRIiSqJ0aaZxOPtysrLy+Mz57LzKtUIae4pmFbIFUKYHMmkhFAbae4pSJASQm3MzO5uxhxvQiRICaEysqSVkmmFXCGEyZFMSgi1KZ1PypjjTYgEKSHURoYgKEhzTwihapJJCaE6Rg5BMLHcQ4KUEGojzT0F0wq5QgiTI5mUEGojgzkVJEgJoTbS3FOQICWE2si7ewqmdTdCCJMjmZQQaiPNPQUJUkKojuZ/mzHHmw5p7gkhVE0yKSHURpp7ChKkhFAbCVIK0twTQqiaZFJCqI50nN9LgpQQaiPNPQVp7gkhVE0yKSHURlp7Co9FJtWrVy/Cw8Nr9Bpz5szB19e3Rq8hHheaathMx2MRpIT4SyntkzJmq4LY2Fjc3d2xsrLCz8+PpKSkB5bfu3cvfn5+WFlZ0bJlS1auXFlh2U2bNqHRaBgwYIDB9ZIgJQAoLCys7SqIWhQXF0d4eDiRkZEcP36c7t27069fPzIyMsotf/HiRUJCQujevTvHjx/nnXfeYeLEiXz99ddlyl66dImIiAi6d+9epbo9dkHqX//6F/7+/tja2tK0aVOGDBlCdna2/vvExEQ0Gg27d+/G39+fevXq0bVrV86cOaM4z8KFC3FycsLW1paRI0dy+/Ztg+qxbt06OnTogKWlJc7OzkyYMEH/XUZGBi+88AI2NjbY2dnx6quv8t///heAM2fOoNFoOH36tOJ8S5Yswc3NDZ1OB8CpU6cICQnBxsYGJycnwsLCyMnJ0Zfv1asXEyZMYMqUKTg6OtK3b1+D6i9qkAYjMynDL7lkyRJGjhzJqFGj8PT0JCYmBhcXF1asWFFu+ZUrV9KiRQtiYmLw9PRk1KhRjBgxgkWLFinKFRcXM3ToUObOnUvLli2r8GM8hkGqsLCQd999l59//plt27Zx8eJFhg8fXqZcZGQkixcv5ujRo9SpU4cRI0bov9u8eTOzZ89m3rx5HD16FGdnZ2JjYytdhxUrVjB+/HjefPNNTpw4wfbt2/Hw8ABAp9MxYMAArl+/zt69e0lISODChQsMGjQIgLZt2+Ln58eGDRsU59y4cSNDhgxBo9Gg1Wrp2bMnvr6+HD16lB07dvDf//6XV199VXHM559/Tp06ddi/fz+ffvppuXUtKCggLy9PsYmaVj19Uvf/cysoKCj3aoWFhaSkpBAUFKTYHxQURHJycrnHHDhwoEz54OBgjh49yp07d/T7oqOjady4MSNHjjTkB1DSPQZ69uypmzRpUrnfHT58WAfobty4odPpdLo9e/boAN0PP/ygL/Pdd9/pAN2ff/6p0+l0uoCAAN2YMWMU5+nSpYvOx8enUvVp1qyZLjIystzvdu3apTM3N9dlZGTo9508eVIH6A4fPqzT6XS6JUuW6Fq2bKn//syZMzpAd/LkSZ1Op9NFRUXpgoKCFOfNzMzUAbozZ87odLq7v4mvr+9D6zp79mwdUGbL1WbodPl/PDZbrjbj7n3n5lbqn3FV5Obm6gDdH2eP60q056u8/XH2eLn/zGbPnl3udS9fvqwDdPv371fsnzdvnq5NmzblHtO6dWvdvHnzFPv279+vA3RXrlzR6XQ63b59+3TNmzfXXb16VafT6XSvv/667oUXXjD4d3nsMqnjx4/zwgsv4Orqiq2tLb169QIo0/bu2LGj/m9nZ2cAfbMwLS2NgIAARfn7P1ckOzubK1eu8Le//a3c79PS0nBxccHFxUW/r3379jRo0IC0tDQAXnvtNS5dusTBgwcB2LBhA76+vrRv3x6AlJQU9uzZg42NjX5r164dABcuXNCf19/f/6H1nTFjBrm5ufotMzOzUvcpjFBNHeeZmZmKf3YzZsx4yGWV7USdTldm38PKl+6/ceMGf//731m9ejWOjo5V+RX0HqtxUvn5+QQFBREUFMS//vUvGjduTEZGBsHBwWU6juvWrav/u/QfRklJidF1sLa2fuD3Ff0f4979zs7O9O7dm40bN/L000/z1Vdf8Y9//ENftqSkhNDQUN5///0y5ykNuAD169d/aH0tLS2xtLR8aDlRnapnoJSdnR12dnYPLe3o6Ii5uTlZWVmK/dnZ2Tg5OZV7TNOmTcstX6dOHRo1asTJkydJT08nNDRU/33pvz916tThzJkztGrVqlJ381hlUqdPnyYnJ4eFCxfSvXt32rVrp+g0ryxPT099FlPq/s8VsbW1xc3Njd27d5f7ffv27cnIyFBkLKdOnSI3NxdPT0/9vqFDhxIXF8eBAwe4cOECr732mv67zp07c/LkSdzc3PDw8FBslQlM4vFiYWGBn58fCQkJiv0JCQl07dq13GMCAgLKlN+1axf+/v7UrVuXdu3aceLECVJTU/Vb//796d27N6mpqYqWwsM8VkGqRYsWWFhYsHz5cn777Te2b9/Ou+++a/B5Jk2axLp161i3bh1nz55l9uzZnDx5stLHz5kzh8WLF7Ns2TLOnTvHsWPHWL58OQB9+vShY8eODB06lGPHjnH48GGGDRtGz549Fc2zgQMHkpeXx9ixY+nduzfNmzfXfzd+/HiuX7/O4MGDOXz4ML/99hu7du1ixIgRFBcXG3y/4hGrhXFSU6ZMYc2aNaxbt460tDQmT55MRkYGY8aMAe42+4cNG6YvP2bMGC5dusSUKVNIS0tj3bp1rF27loiICACsrKzw8vJSbA0aNMDW1hYvLy8sLCwqXbfHKkg1btyY9evX83/+z/+hffv2LFy4sMwj08oYNGgQs2bNYtq0afj5+XHp0iXGjh1b6eNff/11YmJiiI2NpUOHDjz//POcO3cOuNu03LZtGw4ODvTo0YM+ffrQsmVL4uLiFOews7MjNDSUn3/+maFDhyq+a9asGfv376e4uJjg4GC8vLyYNGkS9vb2mJnYmmwmqRaC1KBBg4iJiSE6OhpfX19++ukn4uPjcXV1BUCr1Sr6bd3d3YmPjycxMRFfX1/effddli1bxksvvVRtP0Mpja60t0uISsjLy8Pe3p5cbUal+jtMRV5eHvbOLcjNza2x+y79bf+4cAI7W9uqn+fGDRq08q7Ruj5Kj1XHuRB/DdXTcW4qJEjVABsbmwq/+/7776v8eoB4PGg0mgc++q/M8aZEglQNSE1NrfC7ezu4hSiXTHqnIEGqBpS+4iKEMJ4EKSFUR/qk7iVBSgjVMbK5Z2JBSgbNCCFUTTIpIdRGOs4VJEgJoTrSJ3Uvae4JIVRNMikh1EaaewoSpIRQG2ntKUhzTwihapJJCaE6kkrdS4KUEGojfVIKEqSEUBsJUgrSJyWEUDXJpIRQHemTupcEKSHURoORzb1qq4kqSHNPCKFqkkkJoTbSca4gQUoI1ZE+qXtJc08IoWqSSQmDlC7TmHfjRi3X5NEqvd9HsUxl3s2bRjXZ8m7erMba1D4JUsIgN/73L6tLmw61XJPacePGDezt7Wvk3BYWFjRt2rRaftumTZsatJS5mskKxsIgJSUlXLlyBVtbW5Nb3+1BdDodN27coFmzZjW6VP3t27cpLCw0+jwWFhZYWVlVQ41qnwQpIYSqSce5EELVJEgJIVRNgpQQQtUkSAkhVE2ClBBC1SRICSFUTYKUEELVJEgJIVRNgpQQQtUkSAkhVE2ClBBC1SRICSFU7f8BBn0yw7dxBMoAAAAASUVORK5CYII=",
      "text/plain": [
       "<Figure size 640x480 with 2 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "components = pca.components_ #eigenvectors of pca\n",
    "num_com = components.shape[0]\n",
    "eigvals = pca.explained_variance_ #eigenvalues of pca\n",
    "covariance = np.dot(X_train.T,X_train)/(X_train.shape[0]-1)\n",
    "\n",
    "# checking the SVD decompostion\n",
    "for i in range(components.shape[0]):\n",
    "    temp = np.dot(covariance,components[i])/eigvals[i]\n",
    "    print(temp,components[i],eigvals[i])\n",
    "    # assert np.mean(np.abs(temp-components[i])) < 1e-8\n",
    "    # eigvals.append(np.dot(temp,components[i]))\n",
    "\n",
    "var_ratio = pca.explained_variance_ratio_\n",
    "var = pca.explained_variance_\n",
    "print(\"\\n Varaiance covered as we increasing the number of PCA components: \", np.cumsum(var_ratio))\n",
    "\n",
    "print(np.sum(components[:,0]**2*var_ratio),components[:,0],var_ratio)\n",
    "fig,ax = pl.subplots()\n",
    "weighted_contri = np.abs(components[:,:])**2*var_ratio.reshape((num_com,1))\n",
    "weighted_sum = np.sum(weighted_contri,axis=0).reshape((weighted_contri.shape[1],1))\n",
    "print(weighted_sum,np.sum(weighted_sum))\n",
    "xdata = weighted_sum\n",
    "img = ax.imshow(xdata,cmap='Reds')\n",
    "pl.colorbar(img)\n",
    "ax.set_yticks([0,1,2,3,4,5])\n",
    "ax.set_yticklabels([index_col[i] for i in range(components.shape[1])],fontsize=10)\n",
    "ax.set_ylabel('Features',fontsize=14)\n",
    "\n",
    "# ax.set_xticks([0,1,2,3,4,5])\n",
    "ax.set_xticks([])\n",
    "# ax.set_xticklabels([i+1 for i in [0,1,2,3,4,5]],fontsize=10)\n",
    "# ax.set_xlabel('PCA compenents',fontsize=14)\n",
    "# _ = ax.set_title('Weights of components on features \\n (multiplied by the variance ratio of the components)',fontsize=10)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "how much variance in PCA? Weights in different features using PCA.\n",
    "Conclusion:\n",
    "Decision boundaries, non-trivial points, key take away points, future improvements."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df = pd.DataFrame(classification_report(y_train, y_pred_train_svm, output_dict=True))\n",
    "fig, ax = pl.subplots(figsize=(15, 8))\n",
    "ax.axis('off')\n",
    "table = ax.table(cellText=df.values, colLabels=df.columns, rowLabels=df.index, loc='center', bbox=[0, 0, 1, 1], fontsize=14)\n",
    "pl.subplots_adjust(left=0.1, right=0.8, top=0.9, bottom=0.6)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df = pd.DataFrame(classification_report(y_test, y_pred_test_svm, output_dict=True))\n",
    "fig, ax = pl.subplots(figsize=(15, 8))\n",
    "ax.axis('off')\n",
    "table = ax.table(cellText=df.values, colLabels=df.columns, rowLabels=df.index, loc='center', bbox=[0, 0, 1, 1], fontsize=14)\n",
    "pl.subplots_adjust(left=0.1, right=0.8, top=0.9, bottom=0.6)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Predicting on new solar farms"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "file_loc = 'new_farms.csv'\n",
    "new_data = pd.read_csv(file_loc)\n",
    "new_data = new_data.loc[new_data['Solar_farm_present']==1]\n",
    "print(new_data.columns)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "new_farms = new_data[columns_ML_excluding_LC+['land_cover']]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sd.transform(new_farms[columns_ML_excluding_LC])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "temp = np.hstack((sd.transform(new_farms[columns_ML_excluding_LC]),new_farms['land_cover'].values.reshape((-1,1))))\n",
    "new_farms = pd.DataFrame(temp, columns = columns_ML_excluding_LC+['land_cover'])\n",
    "new_farms_reduced = pca.transform(new_farms)\n",
    "# Xreducedtest = np.vstack([X_test['ac_annual'].values, X_test['solrad_annual'].values]).T\n",
    "\n",
    "new_y_pred_test_svm = svm_tuned.predict(new_farms_reduced)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "new_y_pred_test_svm"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df = pd.DataFrame(classification_report([1]*len(new_y_pred_test_svm), new_y_pred_test_svm, output_dict=True))\n",
    "fig, ax = pl.subplots(figsize=(15, 8))\n",
    "ax.axis('off')\n",
    "table = ax.table(cellText=df.values, colLabels=df.columns, rowLabels=df.index, loc='center', bbox=[0, 0, 1, 1], fontsize=14)\n",
    "pl.subplots_adjust(left=0.1, right=0.8, top=0.9, bottom=0.6)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Decision boundary"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": 976,
   "metadata": {},
   "outputs": [],
   "source": [
    "def make_meshgrid(x, y, h=.02):\n",
    "    x_min, x_max = x.min() - 1, x.max() + 1\n",
    "    y_min, y_max = y.min() - 1, y.max() + 1\n",
    "    xx, yy = np.meshgrid(np.arange(x_min, x_max, h), np.arange(y_min, y_max, h))\n",
    "    return xx, yy\n",
    "\n",
    "def plot_contours(ax, clf, xx, yy, **params):\n",
    "    Z = clf.predict(np.c_[xx.ravel(), yy.ravel()])\n",
    "    Z = Z.reshape(xx.shape)\n",
    "    out = ax.contourf(xx, yy, Z, **params)\n",
    "    return out"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 988,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([4.33747858e+01, 3.73230175e+05, 2.82138881e-01, 2.01669464e+10,\n",
       "       1.48580537e+04])"
      ]
     },
     "execution_count": 988,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "sd.var_"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 989,
   "metadata": {},
   "outputs": [
    {
     "ename": "",
     "evalue": "",
     "output_type": "error",
     "traceback": [
      "\u001b[1;31mThe Kernel crashed while executing code in the the current cell or a previous cell. Please review the code in the cell(s) to identify a possible cause of the failure. Click <a href='https://aka.ms/vscodeJupyterKernelCrash'>here</a> for more info. View Jupyter <a href='command:jupyter.viewOutput'>log</a> for further details."
     ]
    },
    {
     "ename": "",
     "evalue": "",
     "output_type": "error",
     "traceback": [
      "\u001b[1;31mCanceled future for execute_request message before replies were done"
     ]
    }
   ],
   "source": [
    "pca = PCA(n_components=2)\n",
    "Xreduced = pca.fit_transform(X_train)\n",
    "Xreduced = np.vstack([X_train['ac_annual'].values, X_train['min_distance'].values]).T\n",
    "\n",
    "means = sd.mean_[col_index['ac_annual']],sd.mean_[col_index['min_distance']]\n",
    "stds = sd.var_[col_index['ac_annual']]**0.5,sd.var_[col_index['min_distance']]**0.5\n",
    "X0, X1 = Xreduced[:,0]*stds[0]+means[0], Xreduced[:,1]*stds[1]+means[1]\n",
    "xx, yy = make_meshgrid(X0, X1)\n",
    "\n",
    "model = SVC(kernel='rbf',C=4)\n",
    "clf = model.fit(Xreduced, y_train)\n",
    "\n",
    "farms_reduced = np.vstack([new_farms['ac_annual'].values,new_farms['min_distance']]).T\n",
    "\n",
    "# farms_reduced = pca.transform(new_farms)\n",
    "print(farms_reduced.shape,Xreduced.shape)\n",
    "new_y_pred_test_svm = model.predict(farms_reduced)\n",
    "\n",
    "fig, ax = pl.subplots()\n",
    "plot_contours(ax, clf, xx, yy, cmap=pl.cm.coolwarm, alpha=0.8)\n",
    "ax.scatter(X0, X1, c=y_train, cmap=pl.cm.coolwarm, s=20, edgecolors='k')\n",
    "\n",
    "print(new_farms_reduced[:,col_index['ac_annual']].shape)\n",
    "\n",
    "# ax.scatter(new_farms_reduced[:,col_index['ac_annual']],new_farms_reduced[:,col_index['elevation']], c=new_y_pred_test_svm, cmap=pl.cm.coolwarm, s=20, edgecolors='k',marker='s')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "Xreduced"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(np.sum(np.abs(y_pred_train_svm - y_train.values))/len(y_train))\n",
    "# thresholds\n",
    "solar_slopes = solar.loc[solar['Solar_farm_present']==1]['slope'].values.copy()\n",
    "solor_elevations = solar.loc[solar['Solar_farm_present']==1]['elevation'].values.copy()\n",
    "solar_distance = solar.loc[solar['Solar_farm_present']==1]['min_distance'].values.copy()\n",
    "solar_ac_annual = solar.loc[solar['Solar_farm_present']==1]['ac_annual'].values.copy()\n",
    "\n",
    "slope_threshold = np.percentile(solar_slopes,threshold)\n",
    "elevations_threshold = np.percentile(solor_elevations,threshold)\n",
    "distance_threshold = np.percentile(solar_distance,threshold)\n",
    "ac_threshold = np.percentile(solar_ac_annual,threshold)\n",
    "print('slope_th:{}, elevation_th:{}, distance_th:{}, ac_th:{}'.format(slope_threshold,elevations_threshold,distance_threshold,ac_threshold))\n",
    "\n",
    "print(np.sum(np.abs(y_pred_test_svm - y_test.values))/len(y_test))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "X_train"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fake_data = "
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3.9.13 ('base')",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.9.13"
  },
  "orig_nbformat": 4,
  "vscode": {
   "interpreter": {
    "hash": "46d4d6e35a705a453addec98352089b1c875d6193677971858623155cb10eb88"
   }
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}