LogisticRegression.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Baseline algorithm: Logistic Regression"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import pandas as pd\n",
    "from sklearn.model_selection import train_test_split\n",
    "from sklearn.linear_model import LogisticRegression\n",
    "from sklearn import svm\n",
    "from sklearn.metrics import brier_score_loss"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Data pre-processing"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/Users/davidlafferty/anaconda3/lib/python3.7/site-packages/IPython/core/interactiveshell.py:3063: DtypeWarning: Columns (10,17) have mixed types.Specify dtype option on import or set low_memory=False.\n",
      "  interactivity=interactivity, compiler=compiler, result=result)\n"
     ]
    }
   ],
   "source": [
    "# Read data\n",
    "data = pd.read_csv('./data_clean.csv', index_col='date_time', parse_dates=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(381075, 28)"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "data.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Index(['station_id', 'latitude', 'longitude', 'elevation', 'date_time.1',\n",
       "       'wind_direction', 'wind_speed', 'current_wx1', 'current_wx2',\n",
       "       'current_wx3', 'low_cloud_type', 'low_cloud_level', 'medium_cloud_type',\n",
       "       'medium_cloud_level', 'high_cloud_type', 'high_cloud_level',\n",
       "       'highest_cloud_type', 'highest_cloud_level', 'cloud_coverage',\n",
       "       'air_temperature', 'dew_point_temperature', 'altimeter',\n",
       "       'present_weather', 'past_weather', 'past_weather2',\n",
       "       'air_pressure_at_sea_level', 'eastward_wind', 'northward_wind'],\n",
       "      dtype='object')"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "data.columns"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Drop NaNs or constants\n",
    "data.drop(columns = ['date_time.1', 'station_id', 'latitude', 'longitude', 'elevation', 'current_wx2', 'current_wx3', 'air_pressure_at_sea_level'], inplace = True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array(['BR', 'HZ', '-RA', '-FZRA', '-SN', 'SN', 'FG', 'RA', 'FZFG', '+SN',\n",
       "       '+RA', '+TSRA', 'TSRA', '-TSRA', 'FZRA', 'BL', 'TS', '-SNRA',\n",
       "       'VCTS', '+TS', '-TS', 'MIFG', '-SH', 'FZDZ'], dtype=object)"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Unique weather symbols\n",
    "data['current_wx1'].unique()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([10,  5, 61, 66, 71, 73, 45, 63, 49, 75, 65, 97, 95, 67,  0, 17, 68,\n",
       "       13, 12, 80, 57])"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Unique weather codes\n",
    "data['present_weather'].unique()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Choose classifications: 0 for liquid precipitation (rain, feezing rain, thunderstorm + rain, showers, freezing drizzle); 1 for solid (snow)\n",
    "# METAR info here: https://en.wikipedia.org/wiki/METAR\n",
    "# wx_to_prcp = {'RA' : 0, '-RA' : 0, '+RA' : 0, 'FZRA' : 0, '-FZRA' : 0, '+FZRA' : 0, 'SN' : 1, '-SN' : 1, '+SN' : 1, 'TSRA' : 0, '-TSRA' : 0, '+TSRA' : 0, 'SH' : 0, '+SH' : 0, '-SH' : 0, 'FZDZ' : 0, '-FZDZ' : 0, '+FZDZ' : 0}\n",
    "wx_to_prcp = {'RA' : 0, '-RA' : 0, '+RA' : 0, 'FZRA' : 0, '-FZRA' : 0, '+FZRA' : 0, 'SN' : 1, '-SN' : 1, '+SN' : 1, 'TSRA' : 0, '-TSRA' : 0, '+TSRA' : 0}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Simple dictionary lookup to classify weather symbols\n",
    "def classify_prcp(wx_code):\n",
    "    try:\n",
    "        return wx_to_prcp[wx_code]\n",
    "    except:\n",
    "        return np.NaN    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Apply to data\n",
    "data['prcp_type'] = data['current_wx1'].apply(classify_prcp)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Filter for precipitating observations\n",
    "data = data.dropna(subset = ['prcp_type'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(115301, 21)"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "data.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "# One-hot encoding for dates\n",
    "data['month'] = data.index.month\n",
    "data['day'] = data.index.day\n",
    "data['hour'] = data.index.hour"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>wind_direction</th>\n",
       "      <th>wind_speed</th>\n",
       "      <th>current_wx1</th>\n",
       "      <th>low_cloud_type</th>\n",
       "      <th>low_cloud_level</th>\n",
       "      <th>medium_cloud_type</th>\n",
       "      <th>medium_cloud_level</th>\n",
       "      <th>high_cloud_type</th>\n",
       "      <th>high_cloud_level</th>\n",
       "      <th>highest_cloud_type</th>\n",
       "      <th>...</th>\n",
       "      <th>altimeter</th>\n",
       "      <th>present_weather</th>\n",
       "      <th>past_weather</th>\n",
       "      <th>past_weather2</th>\n",
       "      <th>eastward_wind</th>\n",
       "      <th>northward_wind</th>\n",
       "      <th>prcp_type</th>\n",
       "      <th>month</th>\n",
       "      <th>day</th>\n",
       "      <th>hour</th>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>date_time</th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>2000-01-02 01:25:00</th>\n",
       "      <td>90.0</td>\n",
       "      <td>12.0</td>\n",
       "      <td>-RA</td>\n",
       "      <td>SCT</td>\n",
       "      <td>7000.0</td>\n",
       "      <td>OVC</td>\n",
       "      <td>8500.0</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>...</td>\n",
       "      <td>29.79</td>\n",
       "      <td>61</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>-12.000000</td>\n",
       "      <td>-7.347881e-16</td>\n",
       "      <td>0.0</td>\n",
       "      <td>1</td>\n",
       "      <td>2</td>\n",
       "      <td>1</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2000-01-02 01:30:00</th>\n",
       "      <td>90.0</td>\n",
       "      <td>13.0</td>\n",
       "      <td>-RA</td>\n",
       "      <td>OVC</td>\n",
       "      <td>7500.0</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>...</td>\n",
       "      <td>29.79</td>\n",
       "      <td>61</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>-13.000000</td>\n",
       "      <td>-7.960204e-16</td>\n",
       "      <td>0.0</td>\n",
       "      <td>1</td>\n",
       "      <td>2</td>\n",
       "      <td>1</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2000-01-02 01:35:00</th>\n",
       "      <td>100.0</td>\n",
       "      <td>10.0</td>\n",
       "      <td>-RA</td>\n",
       "      <td>OVC</td>\n",
       "      <td>8000.0</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>...</td>\n",
       "      <td>29.80</td>\n",
       "      <td>61</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>-9.848078</td>\n",
       "      <td>1.736482e+00</td>\n",
       "      <td>0.0</td>\n",
       "      <td>1</td>\n",
       "      <td>2</td>\n",
       "      <td>1</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2000-01-02 04:10:00</th>\n",
       "      <td>80.0</td>\n",
       "      <td>11.0</td>\n",
       "      <td>-RA</td>\n",
       "      <td>FEW</td>\n",
       "      <td>6000.0</td>\n",
       "      <td>OVC</td>\n",
       "      <td>8500.0</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>...</td>\n",
       "      <td>29.73</td>\n",
       "      <td>61</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>-10.832885</td>\n",
       "      <td>-1.910130e+00</td>\n",
       "      <td>0.0</td>\n",
       "      <td>1</td>\n",
       "      <td>2</td>\n",
       "      <td>4</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2000-01-02 04:15:00</th>\n",
       "      <td>90.0</td>\n",
       "      <td>14.0</td>\n",
       "      <td>-RA</td>\n",
       "      <td>FEW</td>\n",
       "      <td>6000.0</td>\n",
       "      <td>OVC</td>\n",
       "      <td>8000.0</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>...</td>\n",
       "      <td>29.73</td>\n",
       "      <td>61</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>-14.000000</td>\n",
       "      <td>-8.572528e-16</td>\n",
       "      <td>0.0</td>\n",
       "      <td>1</td>\n",
       "      <td>2</td>\n",
       "      <td>4</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "<p>5 rows × 24 columns</p>\n",
       "</div>"
      ],
      "text/plain": [
       "                     wind_direction  wind_speed current_wx1 low_cloud_type  \\\n",
       "date_time                                                                    \n",
       "2000-01-02 01:25:00            90.0        12.0         -RA            SCT   \n",
       "2000-01-02 01:30:00            90.0        13.0         -RA            OVC   \n",
       "2000-01-02 01:35:00           100.0        10.0         -RA            OVC   \n",
       "2000-01-02 04:10:00            80.0        11.0         -RA            FEW   \n",
       "2000-01-02 04:15:00            90.0        14.0         -RA            FEW   \n",
       "\n",
       "                     low_cloud_level medium_cloud_type  medium_cloud_level  \\\n",
       "date_time                                                                    \n",
       "2000-01-02 01:25:00           7000.0               OVC              8500.0   \n",
       "2000-01-02 01:30:00           7500.0               NaN                 NaN   \n",
       "2000-01-02 01:35:00           8000.0               NaN                 NaN   \n",
       "2000-01-02 04:10:00           6000.0               OVC              8500.0   \n",
       "2000-01-02 04:15:00           6000.0               OVC              8000.0   \n",
       "\n",
       "                    high_cloud_type  high_cloud_level highest_cloud_type  ...  \\\n",
       "date_time                                                                 ...   \n",
       "2000-01-02 01:25:00             NaN               NaN                NaN  ...   \n",
       "2000-01-02 01:30:00             NaN               NaN                NaN  ...   \n",
       "2000-01-02 01:35:00             NaN               NaN                NaN  ...   \n",
       "2000-01-02 04:10:00             NaN               NaN                NaN  ...   \n",
       "2000-01-02 04:15:00             NaN               NaN                NaN  ...   \n",
       "\n",
       "                     altimeter  present_weather  past_weather  past_weather2  \\\n",
       "date_time                                                                      \n",
       "2000-01-02 01:25:00      29.79               61             0              0   \n",
       "2000-01-02 01:30:00      29.79               61             0              0   \n",
       "2000-01-02 01:35:00      29.80               61             0              0   \n",
       "2000-01-02 04:10:00      29.73               61             0              0   \n",
       "2000-01-02 04:15:00      29.73               61             0              0   \n",
       "\n",
       "                     eastward_wind  northward_wind  prcp_type  month  day  \\\n",
       "date_time                                                                   \n",
       "2000-01-02 01:25:00     -12.000000   -7.347881e-16        0.0      1    2   \n",
       "2000-01-02 01:30:00     -13.000000   -7.960204e-16        0.0      1    2   \n",
       "2000-01-02 01:35:00      -9.848078    1.736482e+00        0.0      1    2   \n",
       "2000-01-02 04:10:00     -10.832885   -1.910130e+00        0.0      1    2   \n",
       "2000-01-02 04:15:00     -14.000000   -8.572528e-16        0.0      1    2   \n",
       "\n",
       "                     hour  \n",
       "date_time                  \n",
       "2000-01-02 01:25:00     1  \n",
       "2000-01-02 01:30:00     1  \n",
       "2000-01-02 01:35:00     1  \n",
       "2000-01-02 04:10:00     4  \n",
       "2000-01-02 04:15:00     4  \n",
       "\n",
       "[5 rows x 24 columns]"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "data.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Logistic Regression"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Choose features - this could be updated later\n",
    "data = data.filter(['wind_direction', 'wind_speed',\n",
    "                    'low_cloud_level', 'low_cloud_type', 'cloud_coverage', # Medium or high cloud data contains many many NaNs\n",
    "                    'air_temperature', 'dew_point_temperature', \n",
    "                    'altimeter', 'eastward_wind', 'northward_wind',\n",
    "                    'month', 'day', 'hour',\n",
    "                    'prcp_type']).dropna()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [],
   "source": [
    "# One-hot encoding for cloud type\n",
    "data = pd.get_dummies(data)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Separate labels\n",
    "labels = data['prcp_type'].to_numpy()\n",
    "features = data.drop(columns = ['prcp_type']).to_numpy()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Split the data into training and testing sets\n",
    "train_features, test_features, train_labels, test_labels = train_test_split(features, labels, test_size = 0.3, random_state = 42)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Training Features Shape: (78136, 17)\n",
      "Training Labels Shape: (78136,)\n",
      "Testing Features Shape: (33488, 17)\n",
      "Testing Labels Shape: (33488,)\n"
     ]
    }
   ],
   "source": [
    "print('Training Features Shape:', train_features.shape)\n",
    "print('Training Labels Shape:', train_labels.shape)\n",
    "print('Testing Features Shape:', test_features.shape)\n",
    "print('Testing Labels Shape:', test_labels.shape)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "LogisticRegression(C=1.0, class_weight=None, dual=False, fit_intercept=True,\n",
       "                   intercept_scaling=1, l1_ratio=None, max_iter=100000000000.0,\n",
       "                   multi_class='auto', n_jobs=None, penalty='l2',\n",
       "                   random_state=42, solver='lbfgs', tol=0.0001, verbose=0,\n",
       "                   warm_start=False)"
      ]
     },
     "execution_count": 20,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Train the model!\n",
    "clf = LogisticRegression(max_iter=10e10, random_state=42)\n",
    "clf.fit(train_features, train_labels)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Training data score: 0.9546687826354049\n",
      "Testing data score: 0.9555064500716675\n"
     ]
    }
   ],
   "source": [
    "# How well did we do?\n",
    "print('Training data score: ' + str(clf.score(train_features, train_labels)))\n",
    "print('Testing data score: ' + str(clf.score(test_features, test_labels)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Probabilities for logistic regression\n",
    "train_prob = clf.predict_proba(train_features)\n",
    "test_prob = clf.predict_proba(test_features)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Climatological probability\n",
    "def get_climo(labels, shape):\n",
    "    climo = (np.size((labels == 1).nonzero()))/float(np.size(labels))\n",
    "    prob  = np.zeros((shape))\n",
    "    prob[:,0] = 1-climo\n",
    "    prob[:,1] = climo\n",
    "\n",
    "    return prob"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 35,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Calculate climatological probability for training and testing data\n",
    "train_climo_prob = get_climo(train_labels, train_prob.shape)\n",
    "test_climo_prob = get_climo(test_labels, test_prob.shape)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 36,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Calculate the Brier Skill Score\n",
    "train_bss = 1 - (brier_score_loss(train_labels, train_prob[:,1])/brier_score_loss(train_labels, train_climo_prob[:,1]))\n",
    "test_bss = 1 - (brier_score_loss(test_labels, test_prob[:,1])/brier_score_loss(test_labels, test_climo_prob[:,1]))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 38,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Training Brier Skill Score: 0.8599724342674815\n",
      "Testing Brier Skill Score: 0.8607393010682034\n"
     ]
    }
   ],
   "source": [
    "# Bier skill score (USE THIS ONE!!)\n",
    "print('Training Brier Skill Score: ' + str(train_bss))\n",
    "print('Testing Brier Skill Score: ' + str(test_bss))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Training Brier score: 0.04533121736459506\n",
      "Testing Brier score: 0.044493549928332536\n"
     ]
    }
   ],
   "source": [
    "# Bier score\n",
    "print('Training Brier score: ' + str(brier_score_loss(train_labels, clf.predict(train_features))))\n",
    "print('Testing Brier score: ' + str(brier_score_loss(test_labels, clf.predict(test_features))))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Confusion Matrix\n",
    "from sklearn.metrics import confusion_matrix\n",
    "mat = confusion_matrix(test_labels, clf.predict(test_features))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 40,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAPsAAAEQCAYAAAB2hgJTAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nO3deVxU1fvA8c8Ayi6IyCKCoJi577mQZUpk7ua+pblUiuYualqmFC5fNXMtTXNLC/d9zcQNNcldUxFc0VgEZJNtfn/4Y3KcAQcYBJnn/Xr5ijn33HOfO/Fwz7n33HsVsbGxSoQQxZ5RYQcghHg1JNmFMBCS7EIYCEl2IQyEJLsQBkKSXQgDIclexNSsWZOaNWvmu52AgABsbW05evSoHqLKO33tj8i/Yp/stra2av/s7OyoUKECrVq1YuXKlWRkZBR2iK/M7du3sbW1pU2bNnprs2bNmtja2uqtvYI0ZMgQjd+HnP7p83sqCkwKO4BXxc/PD4CMjAzCwsLYuXMnwcHB/Pnnn6xataqQo/vP9u3b9dLOp59+SufOnSlfvrxe2ssrfe2PPrRp0wY3Nze1sosXL7J7925q1Kihkdwv1n3dKYr7DLqso05sbKxa+eXLl/H29iY5OZndu3fTtGnTwgjvlbp9+za1a9fGy8uLXbt26aXNmjVrcvfuXY3v93Wxbt06fH196dmzJ0uWLCnscApUse/GZ6d69ep4eXkBcPbsWVV5Vrc0JSUFf39/6tatS9myZZkwYYKqTmZmJqtXr+aDDz7Azc0NR0dHmjRpwty5c0lNTdW6vdDQUEaMGEHt2rVxdHTEw8ODFi1aMHv2bLV62sa469atw9bWloCAAE6fPk2HDh1wdXXF1dWVLl26cO7cOY3tvThmX7duHbVr1wbg+PHjat3VgIAA1Xpr166lT58+1K5dGycnJ1xdXfnggw9Yv369WvtZQ4K7d+8CZNv9zW7M/vTpU+bPn4+XlxfOzs6UL18eb29vVq9ejVKpefyxtbWlZs2aJCUlMWXKFGrUqIGDgwN169Zl3rx5WtfRh++//x5bW1sWLFigdXlkZCQODg7Ur19fFcPkyZOxtbVl27ZtbNu2jRYtWuDs7EzFihUZPHiw6jt7UVJSEnPnzqVZs2a4uLhQrlw5WrZsydq1a/WyLwbTjc+tjz/+mAsXLtCyZUtKly6Nu7s7AOnp6fTp04e9e/fi6elJ586dMTU15fjx40ybNo0jR46wadMmTEz++2oPHjzIxx9/THJyMs2bN6dTp04kJiZy9epVAgICGDdunE4xnT17lnnz5vHee+8xePBgQkND2bFjB8ePH2fr1q00atQo23Vr1qzJ559/ztKlS3F1daVXr16qZW+//bbq57Fjx1KlShWaNm2Kk5MT0dHR7N+/nyFDhnDjxg2++uorAGxsbPDz82PJkiXEx8erhknw8u5vWloanTt35tixY3h6ejJgwABSU1PZuXMnX3zxBSdOnGDp0qUa66Wnp/PRRx/x8OFDvL29MTExYdeuXXzzzTckJyczadIknb7H3Ojbty8zZszgl19+YdiwYSgUCrXla9euJTU1lf79+2ss27BhAwcPHqRjx440b96cs2fPEhgYyNGjRzlw4ACurq6quo8fP6Zdu3ZcunSJ+vXr07t3bzIyMjhw4ADDhg3j/PnzGgeG3DLYZL969SrHjx8HoF69ehrL7927x/HjxylTpoxa+bx589i7dy+DBw9mxowZGBsbA8+O9qNGjWLVqlUsX76czz//HIDo6GgGDBhASkoKgYGBeHt7a2xHVwcPHmT27NkMHjxYVbZt2zb69evHsGHDOH36tMYvXJZatWphY2PD0qVLcXNzY+LEiVrrnTx5Eg8PD7Wyp0+f0rlzZ3744QcGDhyIi4sLtra2TJw4kV9//ZX4+Phs29NmwYIFHDt2jBYtWrBhwwZKliwJPDsitmrVig0bNtCqVSs6duyotl5ERAS1atVi69atmJmZAc/OxdSvX5+lS5cybtw4SpQooXMcuihTpgwdO3Zkw4YNBAUF8e6776qWKZVKVq1ahZmZGb1799ZYd+/evWzfvp1mzZqpymbMmMGMGTOYNGkSa9asUZWPGTOGS5cuMWvWLD799FNVeXJyMl27dmXZsmW0bdtWbfu5ZTDd+ICAAAICAvD392fw4MG89957JCcn07ZtW1V3/nmTJk3SSPTMzEyWLl1K2bJlCQgIUCU6gJGREdOmTUOhUPDbb7+pyrOSoV+/fhqJDuTqBFrFihUZOHCgWlmHDh1o1KgRN27c4NSpUzq3lZ0XEx3A1NSUwYMHk56eTlBQUL63kdUt/fbbb1WJDs96C1k9h+xOms6cOVOV6ABly5alTZs2xMfHc+PGjXzHpk3WH9cVK1aolf/xxx+Eh4fToUMH7OzsNNZr1aqVWqIDjBw5krJly7Jr1y7VeY4HDx6wZcsWmjRpopboAObm5kyZMgVA7fcqLwzmyD5z5kwAFAoF1tbW1K5dm65du9K/f3+t9Rs0aKBRdvPmTaKjo/Hw8Mi2S2Vubq72S/fXX38B4OPjk889gCZNmmBkpPn3uWnTppw6dYoLFy7QuHHjfG3j7t27zJ8/nz///JP79++TnJystjwiIiJf7T958oRbt27h4OBA1apVNZZnHbnOnz+vsczGxkY1nHqei4sLoHkSVl/q169P3bp12b17N48ePcLR0RGAlStXAmj8Ac6i7SBiZmZGvXr12LdvHxcvXqRZs2acOXMGpVJJRkaG2vmTLElJSQD5/mNmMMme21+ErP+hz4uJiQEgLCxM9cfjZeLi4gAoV65crravjYODg9bysmXLAhAfH5+v9sPDw2nRogWxsbE0adKEFi1aUKpUKYyNjblz5w7r16/n6dOn+dpGVozZ7YuFhQWlSpXSui+lSpXSuk5WD6sg50wMGjQIX19f1qxZw9ixY4mIiGDv3r1Ur16dt956S+s62e1jVnnWPmb9Xp0+fZrTp09nG0NiYmJ+dsFwkj23tI19s37ZssaVurCxsQGeHRGzzobn1b///qu1PDIyUi2+vFq0aBExMTEsWrRIYwy6ceNGjTPyeZEVY3b7kpSURHx8vNZucWHq3LkzU6ZMYdWqVYwePZo1a9aQnp7OgAEDsl0nu33MKs/6LrL+O3r0aNUwpiAYzJhdH9544w1sbGw4e/ZstpfYXtSwYUMA9u/fn+/tBwcHk5mZqVF+4sQJ4NlJuJw8fzJRm1u3bgHQvn17jWVZJzOza1PXo6q1tTUVK1bk33//5dq1axrLs84J1KlTR6f2XhUzMzP69OnD3bt32bdvH6tXr8bKyopu3bplu4627ywlJYWQkBCMjIxUlySzhownT54smOD/nyR7LpiYmPD5558TGRnJ2LFjVWOp50VHR3PhwgXV5169elGqVClWrVrF4cOHNerfv39f5+2Hhoby888/q5Vt27aNU6dOUbly5RwvvQGULl0ahUKR7XXerEtmL86nP3ToEKtXr9a6TtZJzOza1KZv377As7PvaWlpqvL4+HimTZsGPLv0WdQMGDAAIyMjRo0axb179+jatSvW1tbZ1t+7d6/Gd/n9998TGRlJ69atVRO+KlSoQPv27Tl58iTz58/X+ofz9u3bqj/GeSXd+FwaN24cV65cYfXq1ezfv5933nkHFxcXoqKiCAsLIzg4mEGDBqmOsnZ2dqxYsYKPP/6Yjz76iPfee4/atWuTmJjI9evXOXr0KNHR0Tpt29vbm8mTJ3Pw4EGqV6+uus5ubm7OggULsr3slsXS0pLGjRtz8uRJunfvTp06dTAxMaFp06Z4eXkxcOBA1q1bxyeffEL79u1xdnbm6tWrHDx4kE6dOrF582aNNt977z3Onj1L37598fHxwczMDFdXV3r06JFtHL6+vhw8eJCDBw/StGlTPvjgA9LS0tixYwcPHjygR48eGpfdigJ3d3e8vb1VvbRPPvkkx/offvghH330ER07dsTV1ZWzZ89y5MgRHB0dNU7E/fDDD9y5c4evv/6aNWvW0LhxY8qWLUtERAQ3btwgJCSE+fPnU7FixTzHL8meSyYmJqxevZpNmzaxbt06Dhw4QEJCAnZ2dri6ujJq1CiNX3Rvb2/+/PNPvv/+e44cOcLRo0dV3dncTASpX78+48ePx9/fn59++gl4lmxTpkzRudu7dOlSvvzyS06cOMGBAwfIzMzEz88PLy8vatSowY4dO/D392f//v1kZGRQo0YN1qxZg42NjdZkHzNmDPHx8ezevZv58+eTnp6Ol5dXjslesmRJNm/ezJIlS/j9999Zvnw5RkZGVK1alQkTJqiO/EVR79692b9/Pw0bNnzpsKl79+706NGDefPmsXPnTszMzOjSpQtfffWV2oQaeDZDcN++faxcuZLNmzezfft2UlJScHBwwN3dnWnTptGqVat8xV7s58YXB1nzt/38/HI1eUXo37Rp05g7dy5LliyhZ8+eWutMnjyZhQsXsmrVKjp06PCKI8yejNmF0FFsbCwrV67E3t6eTp06FXY4uSbdeCFeYtu2bVy6dIndu3fz+PFjZs2apTaL73UhyS7ES2zatInt27fj5OSEn5+f2r0JrxMZswthIGTMLoSBeK278bYTPijsEASwbunjwg5B/L82sdnPrZcjuxAGQpJdCAMhyS6EgZBkF8JASLILYSAk2YUwEJLsQhgISXYhDIQkuxAGQpJdCAMhyS6EgZBkF8JASLILYSAk2YUwEJLsQhgISXYhDIQkuxAGQpJdCAMhyS6EgZBkF8JASLILYSAk2YUwEJLsQhgISXYhDIQkuxAGQpJdCAMhyS6EgZBkF8JASLILYSAk2YUwEJLsQhgISXYhDIQkuxAGQpJdCAMhyS6EgTDJboGvr2+uG1MoFCxcuDBfAQkhCka2yR4UFIRCochVY7mtL4R4dbJN9osXL77KOIQQBUzG7EIYiGyP7NkJDg4mKCiIyMhIPvvsMzw9PUlMTOTatWtUrlyZUqVKFUScQoh80jnZU1NTGTBgALt370apVKJQKGjbti2enp4YGxvTpUsXfH19GTt2bEHGK4TII5278QEBAezbt4/Zs2dz5swZlEqlapmZmRkdO3Zkz549BRKkECL/dE72wMBA+vfvz8CBA7Gzs9NYXrlyZcLDw/UZmxBCj3TuxkdGRlKzZs1sl5uampKYmKiXoIqqZp51GPt+b+q7VcHF1oH+q6azKniXarmDtR0zO/niU/UtbC2sCbrxN8N/m8vNyLsAlLYoxTdtB/N+1YZUsHMiKiGOnZeOM3n7UmIS41Xt2FpY80O30bSv1QyA7ReOMvy3OcQlJwBgalKSpb38qOdaharO7hwPvcB784a+wm+i6Kv59TBqTh2uVpb8MJItzm8D0Ev5j9b1ri9ax1/DpgFQa9oIXLu2wtLViczUNGJCrnBhynyiTv5dsMEXEJ2T3dHRMccj99mzZ6lQoYI+YiqyrEzNufQglNXBu1nd/2uN5Vs/n0mmUknHpX7EpSQwumVPDo74gWrTepKUmkI5W3tcbMsyfssirkSE4WJblsU9xrF+wHQ+WDBC1c6vA6bhVtqRDxeOQqlUsrzPJNb0n0r7Jc/OhxgbGZGS9pSFRwJpXb0pthbWr+w7eJ3EXbvFoeZ9VZ+VGRmqnzc7eanVtWtQg+Y7f+TO7/8NReP/CeMv329ICLuHibkZVUb15729y9lR2YeUf6MLfgf0TOdufPv27Vm5ciU3b95UlWVNotmzZw+BgYF89NFH+o+wCNlz+SRfblvKpr8Pk5mZqbassoMrTSrWZOj6WZy5fYXrj+4wZP0szEua0rOhDwCXH9yi808T2HHhKKGR9wi68TfjNi/A+82GWJtZAPCmkzsfVm/Cp+tmcPLWRYLDLvHZrzNpV+tt3nB0AyApNYUh62ex7Ng27sX++2q/hNeIMj2dlEdRqn9Pox6rlj1fnvIoivIdWhL/Txj/Bp1R1Qlft51HfwSTGHaPuCs3CRkdQIlSVpSuU7UwdiffdE52Pz8/XF1deffddxk0aBAKhYK5c+fi7e1N7969qVOnDiNGjHh5Q8WUqUlJAFLSU1VlSqWSp2lpvF2pdrbrlTK35Gl6KkmpTwFo4lGDJymJnLh1QVXneOh5ElKSaFox+2GU0GRV0ZWO94Jof+sQXuvnYulRXms9E0sLKvRow81lv2fbllGJEnh+2p3UuCc8Pne1oEIuUDp3462trdm/fz+LFi1i69atmJmZERwcjIeHBxMnTmT48OGYmZnpvOGZM2fqXNfPz0/nuoXl2sNwwqMj+K7DEAavDSDhaRKjWvbE1c4RZ5syWtexMbdiervPWHZ8OxmZz7qYTjZliEyI1aj7b8JjnEppb0doijp1geD+E4m/dgtTBztqTB6Cz4kN7KreltQY9e+3Qq+2GJmWIGzVFo12yrVpjteGuZhYmJMcEcnh9z95LbvwkMtJNWZmZowZM4YxY8bke8Pbtm1T+3z37l2SkpJwdnYGICIiAgsLC1xdXV+LZE/PzKDzTxP4uc+XxMzZT3pGOgevnWH3pRNa61uUNGPH0P9xPzaS8ZvVbx56/rJmFgUKreVCu4i9QWqfo4PP0/7WQSr268i1eb+oLfMc3I17Ww+pdfOzPDp8ij11OmJqXxrPwd3w+v179jfpQcrDyIIMv0DkegYdQFxcHPfu3QOgfPny2NjY5LqNEyf+S4K1a9eyYcMGlixZgqurK/As+X19fenatWteQiwUIXf+oe53H1PKzJKSJiWISoglePzP/HVHvdtnaWrObt+5ALRdPIanz3X9H8ZF42BdWqPtsla2PHoSU7A7UIylJyYRd/km1pXd1cpta79JmYY1OT9prtb1MpKSSQi9Q0LoHaJPnafd9X14DurKJf/FryBq/crV3PiTJ0/SqlUrPDw8aNasGc2aNcPDw4NWrVqpJW9uzZo1i4CAAFWiA7i6uuLv75+r7n5REZ+SSFRCLJ5lXWlQ4U22nf/vKGNlasHeYd9jbGRM64WjSXyarLbuybBLWJtZ0uS58XmTijWxMrPgxC25OSmvjExLUupND5Ij1I/Inp92JyHsHg8P6vj7a2SEkWnJAoiw4Ol8ZN+/fz+9e/fGysqKgQMH4unpiVKpJDQ0lI0bN9KhQwfWrVuHj49ProOIjIwkJSVFo/zp06fExBSdo5mlqTmeZZ+d5DEyMsLNzpHa5SsTkxjP3ceP6FKvBVEJsdyOeUjNcpWY3200W88HceDqaeBZou//Yj6lzCzpuNQPS1MzLE2fneeISYwnLSOdaw/D2XP5JD/28mPwugAUKPixlx87Lhzj+qM7qliqOrlT0qQE9la2WJmaU7t8ZQDO37vxir+Voqnu7PHc33GYxDsRmDnYUWPKUEwsLbj13Ljc2NwM997tuDprucb6JtaWVBs/mPs7/iA5IhLTsna84dsbi/JOapfnXieK2NhYnQaCXl5epKWlsW/fPkqXVu9mxsTE4OPjg6mpKcePH891ED179uT27dvMnz+fevXqARASEsLIkSOpUKECv/76q9b1bCd8kOtt5ce7levx52jN7tsvJ3fxyerpDH+vG+O8e+NYyo6IuChWn9rD9N0rSMtIz3F9gOZzh3LkRgjwbPLNi5Nqhv32P9WkGoAw/y24l3HWaEcxpHG+9zO31i3VHOsWNq/1cyn7TkNM7W15GvmYqOBzXJgyn/iroao6Fft/xFvLprPN7T2SI9QvYRqbm9F03f8o06g2pmVseRodS8yZi1z+binRpy+8uLkio03s6WyX6ZzsTk5OTJkyJdsn2CxYsIBvv/2Whw8f5jrAqKgohgwZwsGDBzE2NgYgMzOTli1bsmTJEuzt7bWu96qTXWhXFJPdUOWU7Dp3493d3XOcDpuYmJjnGXT29vYEBgZy8+ZNrl+/jlKppEqVKnh6euapPSGEJp2T3c/PjzFjxtCiRQsaNGigtuzMmTMsW7aMOXPm5CsYT09PSXAhCki2ya7tWrqTkxM+Pj7UrVuXSpUqARAaGsrff/9N1apVOXbsGB07dtRpw+PHj+frr7/G0tKS8ePH51h31qxZOrUphMhetsm+YsWKbFcKCQkhJCRErezKlStcvXqV//3vfzpt+MqVK6Slpal+zo48xFII/cg22R8/LtiTLjt37tT6sxCiYMgDJ4UwEHmaLqtvMmYXouDlKtkPHz7MggULOHfuHHFxcVpvzMjLjLcXx+zp6elcv36d9PR0atfO/vZQIYTudE72PXv20KdPHypXrkynTp1YsWIFXbt2RalUsmvXLjw9PWndunWegtA2Zk9JSWH48OE0adIkT20KIdTpPGafM2cONWvW5NixY0yaNAmA3r17s2zZMk6cOMGDBw9Ul+P0Iet22vxeuxdCPKNzsl++fJmuXbtiYmKimtKa8f/P9HJ3d2fAgAHMmzdPr8FFRUWRkJDw8opCiJfSuRtvamqKubk5AJaWligUCiIj/7td0MXFhbCwsDwF8eKbX5VKJY8ePSIwMDBPd9EJITTpnOweHh7888+zx++WKFGCKlWqsGPHDrp37w7A7t27cXJyylMQP/30k9pnIyMj7O3t6d27N6NGjcpTm0IIdTonu7e3N6tWrcLf358SJUowZMgQRowYobolNSwsjGnTpuUpiAsXiu4tg0IUFzrf4pqWlsaTJ08oXbq0agrrxo0b2bJlC8bGxnz44Yf07NmzQIN9kdziWjTILa5Fh17uZ3+ZqKgoHj58SI0aNXK9bnb3yGuzaNEi1c+S7EWDJHvRoZf72V/ml19+4bvvvsvTpJqoqChOnjyJQqGgWrVqAFy9epXMzEy5zi6EnhSJ6bKNGjXC3NycRYsWYWlpCTx7GMbw4cOpVq2avAZaCD0oEjfC/Pjjj0yYMEGV6PDs8t64ceM0ztQLIfKmSCR7YmKi1mfXPXr0iOTkZC1rCCFyq0gke7t27fD19WXTpk3cvn2b27dvs2nTJoYNG0bbtm0LOzwhioUcx+xnz57VuaEHDx7kOYi5c+cyefJkhg4dqnp6jYmJCX379mX69Ol5blcI8Z8cL709f039ZZRKJQqFIl8vdUhMTCQsLAylUknFihXVxvDayKW3okEuvRUdeb709vw17VfB0tIyT9fphRAvl2Oy9+rVq8A23KNHD3766SdKlSpFjx49cqy7YcOGAotDCENRaNfZ7ezsVEOE3AwXhBB5o7fpsoVBxuxFg4zZi45XMl02t17Wdc+iUChYv359AUcjRPFXqN14IcSrU2jJvnix9lcXCyEKRpGYQSeEKHjZHtnv3r2bpwZdXV3zHIwQouBkm+y1atXK0+Ww/MygE0IUnGyTfeHChWrJrlQqWbp0KXfu3KFbt254enqiVCq5efMmGzduxM3Njc8+++yVBC2EyL1sk713795qn7///nuSkpL4+++/Nc6kT5gwAR8fH6KjowsmSiFEvul8gm758uX0799f6yUze3t7+vXrx7Jly/QanBBCf3RO9qioKNXtp9qkp6cTFRWll6CEEPqnc7LXqlWL5cuXc/v2bY1l4eHhLF++nFq1auk1OCGE/ug8qebbb7+lU6dOvPXWW3z44Yd4enoCcOPGDfbu3YuJiQn+/v4FFqgQIn90TvaGDRty6NAh/P39OXDgANu2bQPAwsICHx8fJk2aRNWqVQssUCFE/uRqumyVKlVYs2YNmZmZREVFoVQqKVu2LEZGMhFPiKIuT3PjlUolJiYmlCpVShJdiNdErjI1JCSEjh07Uq5cOTw9PTl+/DgA0dHRdOvWjSNHjhRIkEKI/NP5yH769Gnat2+Po6MjPXr0YPXq1aplZcqUISEhgTVr1vDuu+8WSKDabPpFHppQFNje7VzYIQgd6Hxknz59OpUqVeLUqVN89dVXGsubNWvGX3/9pdfghBD6o3Oyh4SE0KdPH8zMzLTeIOPi4sKjR4/0GpwQQn90TnYjI6McT8Y9evQIc3NzvQQlhNA/nZO9Tp067N27V+uy1NRUAgMDeeutt/QWmBBCv3RO9tGjRxMUFMSwYcO4ePEiAA8fPuTgwYO0b9+esLAwxowZU2CBCiHyJ1ePkt64cSPjxo0jLi5O9bonpVKJjY0N8+fPp0OHDgUZq4ZDTtKTKApsb8jZ+KKivrVftstyNammS5cutG7dmsOHDxMaGkpmZiYeHh60bNkSKyurfAcqhCg4Oif78ePHqVKlCvb29rRp00ZjeXR0NNeuXcPLy0uvAQoh9EPnMXu7du04fPhwtsuPHDlCu3bt9BKUEEL/dE52pTLnoX1qaqrMkxeiCMuxGx8fH09cXJzqc0xMjNZHTMfGxrJp0yacnZ31H6EQQi9yTPbFixcza9Ys4Nk71yZOnMjEiRO11lUqlUyZMkX/EQoh9CLHZG/evDlmZmYolUqmTZvGRx99RM2aNdXqKBQKLCwsqFu3Lg0aNCjQYIUQeZdjsjdu3JjGjRsD8PTpU9q1a0f16tVfSWBCCP3S+dLbhAkTCjIOIUQB0/n0+cSJE6lXr162y+vXry9jdiGKMJ2Tff/+/Xz00UfZLu/UqVO2N8oIIQqfzsl+//593Nzcsl3u5ubG/fv39RKUEEL/dE52a2trwsPDs10eFhaGmZmZPmISQhQAnZP9nXfeYcWKFVoTPjw8nJUrV/LOO+/oMzYhhB7pfDZ+0qRJHDhwAC8vL3r16kW1atVQKBRcvnyZ9evXY2xszJdfflmQsQoh8kHnZK9UqRL79u1j7NixLF++XG2Zl5cXs2bNonLlynoPUAihH7m6n71q1ars2rWL6OhowsPDUSqVVKxYUetrnIUQRUue3ghTpkwZypQpo+9YhBAFKNtkz3rbS9bDKLI+v4w8vEKIoinbZG/bti0KhYKHDx9SsmRJ1efsZD2TLiYmpkACFULkT7bJvmPHDgBKliyp9lkI8XrKNtnffvvtHD8LIV4v8hwpIQxEtkd2X1/fXDemUChYuHBhvgISQhSMbJM9KChI44RccnIyUVFRANja2qJUKlXPqLO3t8fCwqIAQxVC5Ee2yZ71iqcsN27coFOnTowePZqhQ4eqrrNHR0ezaNEifv/9dzZt2lSw0Qoh8kzn1z916tQJFxeXbLvpvr6+REREsHnzZr0GmBN5/dqC1DoAABWgSURBVFPRIK9/Kjpyev2TzifoTp8+neOTaurVq8fp06dzF5kQ4pXROdmtrKw4duxYtsuPHj0qY3YhijCdk71Hjx5s2bKFESNGcPXqVdLS0khLS+Pq1at88cUXbN++nZ49exZkrEKIfND5RpjJkycTGRnJ6tWrWbNmjepMvVKpRKlU0r17d3ngpBBFWK7ezw5w5coV9u7dy71791Aqlbi5ueHj41Moz5OXE3RFg5ygKzr09n52gGrVqlGtWrV8BSSEePVynezBwcEEBQURGRnJZ599hqenJ4mJiVy7do3KlStTqlSpgohTCJFPOid7amoqAwYMYPfu3arbWdu2bYunpyfGxsZ06dIFX19fxo4dW5DxCiHySOez8QEBAezbt4/Zs2dz5swZtfe1m5mZ0bFjR/bs2VMgQQoh8k/nZA8MDKR///4MHDhQ6zPnKleunONz5YUQhUvnbnxkZKTG65qfZ2pqSmJiol6Cet2ZWFlS7esRlGvvjWnZMsSev8KFsd/x+Oyz+w3KdXgfj4Hdsa1THdOydgT59CXq6H+zDy3cXGj1zx9a2744aRY35v38SvajKLtxIZIDv9/gzvXHxEWn8PG4+jRp5Q5ARnom21dc5vLph0RGJGJmUYI36pSl06Aa2Dk+m/gV/TCRyb21v66s06c18OleBYDIBwlsWnqR0EtRpKdlUq2hI92H1aGU3bMXolw/F8m8MUFa2xn0VSPqv1tez3uedzonu6OjY45H7rNnz1KhQgV9xPTaq7fEn1I1qvDXoAkk33+IW8/2vL1rJQfqtSblwb8YW1gQHfw3d9bvoOGKWRrrJ92LYJe7+rP8yrV/nzrff8X9zfte1W4UaU+T0ynnXorG77vxy8y/1JalpmRw50YsrXq/SXlPW1IS09i45AILJhxj8nJvjI2NKF3WghmBbdTWO3fsPr/9cI5675RXbeOH8cdw8SjFiNnvoFDAjpWXWTz5BOMXvoeRkYKK1ctotHN4803+3BpK9becCvZLyCWdu/Ht27dn5cqV3Lx5U1WWNbFmz549BAYG5vjiR0NhZGZKuY4+XJ4yh6ijp0m8dYer3y4kIfQ2FQf3AuDu+m1c+24Rj/ZrPyKQmcnTR1Fq/8p1eJ9//zhB0u17r3Bviq4ajZzpOKgG9d4tj9ELj0Y0tyrBiNnNaPCeK06u1ri/aUfvUfV4eOcJD28/AcDIWIGNnZnav3NHH/BmPQfsnS0BCL0cTfTDRPqOb0D5Sja4VLShn19D7lx/zD9//wuASQkjjXb+Pnqfhi3KY2aep4c3Fxido/Hz8yMoKIh3332XRo0aoVAomDt3LtOmTSMkJIT69eszYsQIndqaOXOmzgH6+WU/SaAoMjIxwcjEhIyUp2rlGSlPKdM0+xuJcmJRoTwO7zXhVJ+R+gjRICUnpQFgYV1C6/KoiET++ftfBk1ppCpLT80ABZQoaawqMylphEKhIPRSNFXrO2q0c/1cJP/eS2DApKI34UvnZLe2tmb//v0sWrSIrVu3YmZmRnBwMB4eHkycOJHhw4fr/GLHbdu2qX2+e/cuSUlJODs7AxAREYGFhQWurq6vXbKnJyQSHRzCmxOGEH/lOikPo3Dt1pYyjeqQEHonT216DOjK06jHROw4pOdoDUN6Wiabll6kZhNnSpfVfrPWsV1hWNqYUturnKrMo1oZTM1N2PzjBToNfna+auvyS2RmKomLTtHaztFdYZSvZEOFKqX1vyP5lKt+hpmZGWPGjGHMmDH52uiJEydUP69du5YNGzawZMkSXF1dgWfJ7+vrS9euXfO1ncLy14Dx1PvxO1qHHiUzPZ3Yc1e4+/subOvkfuahwtgYt74fcWftFpTp6QUQbfGWkZHJyoDTJCekMmR6k2zrBO+/TROfChib/DeytbY1ZfBXjVn//d8Ebb+FQqGgQQtXXCvbYmSs+Vj1xPhUzh29T5chtQpsf/JDp2RPTk7GxcWFSZMm6X3SzKxZs1i3bp0q0QFcXV3x9/enV69e9O3bV6/bexUSw+5y1KcvxhbmlChlRcrDSN5aM4+k8NyPt53bvIe5swPhKwMLINLiLSMjkxX+p7kfFseoue9iZWOqtd7FkxHERafg1dpdY1m1Bo5MX9uKhLinGBkrsLAqiV+XnZRx0uwhBO+/jcJIwVst3fS9K3qh0wk6c3NzypYtWyBTYSMjI0lJ0ewSPX369LV/4URGUjIpDyMpYVsKB++3ebAz991w90+6ERl0ioSb4foPsBjLSM9k+fRT3LsVx6g572Bjl/0Q89iuMCrXtsfR1TrbOlY2plhYleTa3//yJPYptZqW06hzfHcY9d8tj7mV9vMChU3ns/GdOnViy5YtZGZm6jWA5s2bM2LECM6cOUNGRgYZGRmcOXOGkSNH0rx5c71u61Vx8H4bR593np1Ya9GUZntXk3AjjNurnz2yq0RpG2xqvUmpas/eemtVyQ2bWm9i6miv1o65qzOO778tR3UtUpLTuXszlrs3Y8lUQsy/ydy9GUvMoyQyMjJZ9k0wYVdiGPjlW6BQEBeTQlxMCqlPM9TaiXmUxJW/HvF2aw+t2zmxN5xbl6OJfJDAqQN3WP7NKVp0rozTC38Ybl6MIuL2E95u415Qu5xvOo/Z27RpQ1BQEK1ateLjjz/G3d0dc3NzjXr169fPVQALFixgyJAh+Pj4YGz87KxnZmYmLVu25IcffshVW0VFCRtrqk8bjbmLE2kxsdzftp/LX89Tjbmd27SgwbIZqvr1lnwLwFX/BVz99r9n/Ln360Ja3BPub5Fr6y+6889jtcksO1ddYeeqKzT2qUDbflU5fyICgIAh6pOTnp98A3B8TzjmliWo+46L1u08uvuEbcsvkfgklTKOlrTqXYWWXTRfTX5sVxhObtZUqmGvpZWiQef72UuXVj+7+OJjpvP7rrebN29y/fp1lEolVapUwdPT86XryP3sRYPcz1506OV+9oULF+b4Ysf88vT01CnBhRB5o3Oy9+7dW28bHT9+PF9//TWWlpaMHz8+x7qzZmlOJxVC5N5Lk/3p06fs3r2b8PBw7Ozs+OCDD3Byyt+c3ytXrpCWlqb6OTsF2ZMQwtDkmOyPHj2idevWhIWFqe5ft7Cw4Pfff8fLyyunVXO0c+dOrT8LIQpOjpfe/P39CQ8PZ+jQofz2228EBARgZmb20q63EKLoyfHI/scff9CzZ0/8/f1VZQ4ODgwaNIj79+/j4qL9csXL5OaPhYzZhdCPl3bjGzVqpFbWuHFjlEol9+7dy3Oy5zROf56M2YXQnxyTPSMjQ+NOtqzP2qa46krG6UK8ei89Gx8eHs7Zs2dVn+Pj44Fnr3C2srLSqJ/bGXTPS0hIQKFQYGlpmec2hBDa5TiDrnTp0lq70lmz5bSV5WUG3bJly5g/fz4PHjwAoFy5cowcOZJBgwbluJ7MoCsaZAZd0ZHnGXSLFi3SezAvmjNnDvPmzWPYsGE0afLsfuMTJ07wzTff8OTJE0aNGlXgMQhhCHL9rjd9q1GjBlOnTqVLly5q5b///jvTpk3j0qVL2a4rR/aiQY7sRUdOR3adb3EtKJGRkdSrp/lstvr16xMZGVkIEQlRPBV6sleqVInAQM37tQMDA+XGGCH0qNCfdTthwgQ++eQTTpw4oXpqbXBwMMeOHWPVqlWFHZ4QxUahJ3v79u05dOgQixcvZt++fSiVSt58803++OMPateuXdjhCVFsFHo3/tq1a1haWvLTTz9x5MgRpk6dCsChQ4fIyMjIeWUhhM4KPdmHDx/OhQsXALh//z59+vTh8ePHLF++XG1OvhAifwo92f/55x9Vd33r1q3Uq1ePwMBAli5dysaNGws5OiGKj0JP9szMTEqUePbo3aCgIHx8fADw8PCQS29C6FGhJ3vVqlVZsWIFJ06c4MiRI7Rs2RJ49gqoMmXKFHJ0QhQfhZ7sU6dOZdWqVbRt25bOnTtTvXp14NmbYbVNthFC5E2hX3rz8vIiNDSUJ0+eYGtrqyrv378/FhbaX8InhMi9Qk92AGNjY7VEB6hQoUIhRSNE8VTo3XghxKshyS6EgZBkF8JASLILYSAk2YUwEJLsQhgISXYhDIQkuxAGQpJdCAMhyS6EgZBkF8JASLILYSAk2YUwEJLsQhgISXYhDIQkuxAGQpJdCAMhyS6EgZBkF8JASLILYSAUsbGxysIOQghR8OTILoSBkGQXwkBIsgthICTZ82jIkCF0795d42d9atOmDePGjcuxTpMmTQgICND7tg3Ri993cfv+i8QbYV53M2bMQKnU/3nOtWvXYmIi/4sKS3H7/ovPnhQiGxubAmm3dOnSBdKu0E1x+/6lG68HL3bjk5KSGDJkCC4uLlSuXJk5c+bQvXt3hgwZoqpTs2ZNFixYoNbOy7qRkZGR9OzZEycnJ2rUqMGaNWsKcK9eT8ePH8fb2xsXFxfc3Nxo2bIlV65cAWD79u00bdoUBwcHqlevzv/+978ce2TF7fuXI3sBmDJlCn/++SerV6/G2dmZmTNncuLECdq2bZuvdocOHcrdu3fZunUr5ubmTJo0iTt37ugp6tdfeno6vXr1om/fvixbtoy0tDTOnz+PsbEx586do3///owdO5Zu3boREhLCqFGjsLa25rPPPtOp/df9+5dk17OEhATWrFnDwoULadmyJQCLFi2iWrVq+Wr35s2bHDhwgL1799K4cWMAlixZQp06dfIdc3Hx5MkT4uLiaNWqFR4eHgC88cYbAAwePBgvLy8mTZoEgKenJ6GhocyfP1+nZC8O37904/UsLCyM1NRU3nrrLVWZlZUV1atXz1e7//zzD0ZGRtSvX19V5ubmhrOzc77aLU5Kly5Nr1696Ny5M926dWPhwoXcu3cPePb9NWrUSK1+kyZNePDgAfHx8S9tuzh8/5LseqbrWXkjIyONuunp6flu19AtXryYgwcP0rRpU/bs2UODBg04dOgQSqUShUKhdZ3syp9XHL5/SXY9q1ixIiVKlODMmTOqssTERNVJoiz29vY8fPhQ9TklJYXr169n226VKlXIzMwkJCREVXb37l0iIiL0GH3xULNmTUaOHMmuXbt4++23Wb9+PW+++SbBwcFq9U6ePImLiwvW1tYvbbM4fP+S7HpmZWVF3759mTp1KocPH+bq1asMGzaMzMxMtXrvvPMOgYGBHD16VFUnpyN75cqV8fb2ZtSoUZw+fZoLFy4wdOhQzM3NC3qXXhvh4eFMnTqVU6dOcefOHYKCgrh8+TJVqlTB19eX48ePExAQwM2bN/n9999ZtGgRX3zxhU5tF4fvX07QFYDp06eTmJhInz59MDc359NPPyUpKUmtzqhRo7hz5w69e/fG0tKSMWPGvPQosXjxYr744gvat29PmTJl8PPzIyoqqiB35bViYWHBzZs36d+/P9HR0Tg4ONC1a1dGjhxJiRIl+OWXX5gxYwZz587FwcGBkSNH8umnn+rc/uv+/cstrq9I9+7dsbOzY8mSJYUdijBQ0o0XwkBIsgthIKQbL4SBkCO7EAZCkl0IAyHJLoSBkGQXRc7BgwextbVl586demvz+vXr2NrasmnTJr21+bqRZC+CbG1tdfq3bt26wg5VZcCAAbi5uRV2GCIHMoOuCPrxxx/VPv/yyy/89ddfLFy4UK38xbu4hMiJJHsR9OLDK//8809CQkJ0fqhleno6mZmZlCxZsiDCE68p6ca/5rLGoosXL2bJkiXUrVsXR0dHzp8/rxr7Pn8HHjy7w87W1pZ58+aplT948ABfX1/eeOMNHBwcaNCgAYsWLdLb7Z1hYWGMGjWKBg0a4OzsTIUKFejVqxc3btzQWj8jI4Np06ZRpUoVnJ2d6dKlC2FhYRr1/vnnH/r164eHhweOjo54eXnx22+/6SXm4kSO7MXEmjVrSE5Opl+/fpiZmWFvb09cXJzO60dERNCyZUuMjIwYOHAgDg4OHD16lC+//JLIyEimTp2a7xhPnz7NqVOn6NChA25ubty7d48VK1bQunVrgoODKVOmjFr9mTNnAvDFF18QExPD0qVLadeuHcePH1c95PPSpUu0bt0aZ2dnRowYgZWVFXv27OGzzz4jISGBgQMH5jvu4kKSvZi4f/8+ISEh2Nvbq8pCQ0N1Xn/q1KkoFAqOHj2KnZ0dAJ988gn29vYsXLiQzz77LN9PZWnXrp3GUKRLly54eXnx66+/Mnz4cLVlcXFxBAcHq+43b9SoEd26dWPp0qX4+fkBMG7cONzc3Dh06BCmpqYADBo0iJ49e+Lv70+fPn1U5YZOuvHFRIcOHdQSPTfS09PZuXMnH374IUqlkujoaNU/b29v0tPTOXnyZL5jtLCwUP2cmJhITEwM9vb2uLu7c/78eY36PXv2VHuwhI+PDx4eHhw4cACAR48ecfLkSTp37kxCQoJa3O+//z6PHz/m4sWL+Y67uJAjezGR9YDFvIiIiCAxMZGff/6Zn3/+WWudyMjIPLefJSkpCX9/fzZu3Mi///6rtkxb/JUqVdIo8/T05MKFC8Czh0ACfPPNN3zzzTdat6mPuIsLSfZiwszMTKMsu2erZWRkqH3OeopOz549sz3j7+npmc8IYfTo0QQGBjJkyBAaNmyItbU1RkZGjBkzRuNJPqA9/udPFmatM3LkSJo3b651mzVq1Mh33MWFJHsxZmtrC6Bxou7FZ507Oztjbm5OZmZmtkmjD1u3bqVfv374+/urypRKJY8fP8bd3V2jftaR+3m3bt3C1dUV+K83UKJEiQKNu7iQMXsx5u7urjrp9rxly5apfS5ZsiRt27Zly5YtGg/GBIiNjc3x+Xi6UCqVGBsba1zGW7duHTExMVrXWb9+PU+ePFF93r9/P7du3cLb2xuA8uXL06hRI37++We1h3dmeZ0eGfUqyJG9GCtTpgzt27dn0aJFZGRk4O7uzuHDh7l//75G3enTpxMcHIy3tzd9+/alatWqxMXFcfnyZXbs2MHVq1dVPYXspKamMnv2bI1yS0tLhg4dSqtWrVizZg1mZmZUqVKFc+fOsWPHDtWR+kU2Nja0atWKXr16qS69lStXjs8//1xV5/vvv+fDDz+kSZMmqmvtUVFRnDt3jhMnTuTqikRxJ8lezM2ZMwelUsnKlSsxNjbmww8/ZP369VStWlWtnpOTE4cPH2bWrFns3r2bFStWYGtri6enJ5MnT8bKyuql20pJSeHbb7/VKHdwcGDo0KHMmTMHMzMzAgMDSU5Opl69emzevJkxY8Zobc/Pz49z584xf/584uPjadq0KbNnz1b7o1O1alX+/PNPZs2axYYNG4iOjsbe3p6qVatme9LOUMmTaoQwEDJmF8JASLILYSAk2YUwEJLsQhgISXYhDIQkuxAGQpJdCAMhyS6EgZBkF8JASLILYSD+DwFP+AAwFm7eAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Plot matrix\n",
    "import matplotlib as mpl\n",
    "import matplotlib.pyplot as plt \n",
    "import seaborn as sns\n",
    "mpl.style.use('fivethirtyeight')\n",
    "\n",
    "labels_str = ['liquid', 'solid']\n",
    "sns.heatmap(mat.T, square=True, annot=True, fmt='d', cbar=False,\n",
    "            xticklabels=labels_str, yticklabels=labels_str, cmap='RdYlGn')\n",
    "plt.title('Precipitation Type')\n",
    "plt.xlabel('True Label')\n",
    "plt.ylabel('Predicted Label')\n",
    "plt.tight_layout()\n",
    "plt.savefig('./logreg_conf_mat.png', dpi=600)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## SVM"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "SVC(C=1.0, break_ties=False, cache_size=200, class_weight=None, coef0=0.0,\n",
       "    decision_function_shape='ovr', degree=3, gamma='scale', kernel='rbf',\n",
       "    max_iter=-1, probability=False, random_state=None, shrinking=True,\n",
       "    tol=0.001, verbose=False)"
      ]
     },
     "execution_count": 25,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Train the model!\n",
    "clf = svm.SVC()\n",
    "clf.fit(train_features, train_labels)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Training data score: 0.9065731544998464\n",
      "Testing data score: 0.9070413282369804\n"
     ]
    }
   ],
   "source": [
    "# How well did we do?\n",
    "print('Training data score: ' + str(clf.score(train_features, train_labels)))\n",
    "print('Testing data score: ' + str(clf.score(test_features, test_labels)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Training Brier score: 0.09342684550015358\n",
      "Testing Brier score: 0.09295867176301958\n"
     ]
    }
   ],
   "source": [
    "# Bier score\n",
    "print('Training Brier score: ' + str(brier_score_loss(train_labels, clf.predict(train_features))))\n",
    "print('Testing Brier score: ' + str(brier_score_loss(test_labels, clf.predict(test_features))))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "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.7.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}