import numpy as np from IPython.display import Image import pandas as pd import seaborn as sns import geopandas as gpd from shapely.geometry import Point, Polygon from geopandas import GeoDataFrame import matplotlib.pyplot as plt import sklearn import sklearn.utils from sklearn.pipeline import Pipeline from sklearn.preprocessing import StandardScaler from sklearn.cluster import DBSCAN import matplotlib.pyplot as plt import scipy.stats as stats from datetime import datetime import datetime from datetime import timedelta from sklearn.model_selection import cross_val_score, KFold, cross_validate, RandomizedSearchCV, GridSearchCV from sklearn.model_selection import train_test_split from sklearn.neighbors import KNeighborsClassifier from sklearn import metrics import xgboost as xgb from sklearn.metrics import accuracy_score from xgboost import XGBClassifier,plot_tree from sklearn.multioutput import MultiOutputClassifier from sklearn.svm import SVC

df = pd.read_csv("/work/Data/newdataset.csv", index_col= [0,1]) df = df.dropna() df = df.reset_index() df = df.drop(columns = ["level_0"]) df = df.set_index(["vin","sequence"]) df = df.sort_index() df.sort_values('dwell time duration (seconds)', ascending = False) df['datetime'] = pd.to_datetime(df['last'], unit = 's').dt.hour df['daylight'] = ((df.datetime > 7) & (df.datetime < 19)) * 1 df

## Geofencing again to only keep points in the city, set boundary coords for mapping df = df[(df["latitude"] <= 40.5) \ & (df["latitude"] >= 39.5) \ & (df["longitude"] >= -84)\ & (df["longitude"] <= -82.5)] y_boundaries = max(df['latitude']), min(df['latitude']) x_boundaries = min(df['longitude']), max(df['longitude'])

## Drop columns that interfere with dbscan clustering dbs_data = df[['latitude', 'longitude', 'dwell time duration (seconds)','daylight' ]] dbs_data.head(1)

## Drop values with dwell times shorter than 5 minutes, drop dwell times with values greater than 200000 dbs_data = dbs_data[dbs_data["dwell time duration (seconds)"] >= 300] dbs_data = dbs_data[dbs_data["dwell time duration (seconds)"] <= 200000] dbs_data.sort_values('dwell time duration (seconds)', ascending = False)

sns.histplot(data = dbs_data, x = 'dwell time duration (seconds)', kde = True)

sns.histplot(data = dbs_data, x = 'daylight')

def gen_dbs_clusters(df1, eps, min_samples): """ Generates the clusters based off of location and dwell time using DBScan Parameters: _ _ _ _ _ _ _ _ _ _ _ _ _ df1 = DataFrame, containing lat, long, dwell time interval eps = Float, epsilon value for DBScan min_sample = Integer, mininimum number of samples in a given cluster for DBScan Returns _ _ _ _ _ _ _ _ _ _ _ _ _ A DataFrame containing lat, long, dwell time interval, associated DBS Cluster """ df = df1.copy() ## eps .14, .2, .17 ## min_samples = 2,3,4,5 scaled = StandardScaler().fit_transform(df) dbs = DBSCAN(eps = eps, min_samples = min_samples).fit(df) labels = dbs.labels_ df['dbs_cluster'] = labels n_clusters = len(set(labels)) - (1 if -1 in labels else 0) print(n_clusters, 'clusters') n_noise = list(labels).count(-1) print(n_noise, 'noisey values') # n_clusters, "clusters" n_noise, 'noisey values' return df # eps_vals = [.2, .15]

def plot_clusters(df1): """ Plots the DBScan clusters on a geographical map Parameters: _ _ _ _ _ _ _ _ _ _ _ _ _ df1 = DataFrame, containing lat, long, dwell time interval, associated dbs cluster for each point Returns _ _ _ _ _ _ _ _ _ _ _ _ _ A geographic visualisation of the data points and their associated cluster """ df = df1.copy() geometry = [Point(xy) for xy in zip(df['longitude'], df['latitude'])] gdf = GeoDataFrame(df, crs="EPSG:4326", geometry=geometry) ohio_shape = gpd.read_file('/work/Data/cb_2018_39_bg_500k.zip') ax = ohio_shape.boundary.plot(figsize=(15, 15), alpha=0.1,) gdf.reset_index().plot(ax = ax,cmap = 'jet', alpha = 0.2,column='dbs_cluster',legend=True) plt.xlim(x_boundaries[0], x_boundaries[1]) plt.ylim(y_boundaries[0], y_boundaries[1])

def gen_bins_pointwise(df1, min_interval, max_interval): """ Assigns interval bins with ranges (0, min_interval),(min_interval, max_interval),(max_interval,∞) to each cluster in the inputted DataFrame Parameters: _ _ _ _ _ _ _ _ _ _ _ _ _ df1 = DataFrame, containing lat,long, dwell time duration (seconds) min_interval = Integer, minimum interval value max_interval = Integer, maximum interval value Returns _ _ _ _ _ _ _ _ _ _ _ _ _ A DataFrame containing lat, long, dwell time interval, associated interval bins """ df = df1.copy() df['0-{}hr'.format(min_interval)] = 1* (df["dwell time duration (seconds)"]/3600 <= min_interval) df['{}-{}hr'.format(min_interval, max_interval)] = 1* ((df["dwell time duration (seconds)"]/3600 > min_interval) & (df["dwell time duration (seconds)"]/3600 <= max_interval)) df['+{}hr'.format(max_interval)] = 1 * (df["dwell time duration (seconds)"]/3600 > max_interval) return df

dbs_df = gen_bins_pointwise(dbs_data, 3, 6) dbs_interval = dbs_df.drop(columns = ['dwell time duration (seconds)']) dbs_interval = gen_dbs_clusters(dbs_interval, .0075, 10) dbs_interval_noiseless = dbs_interval.loc[dbs_interval['dbs_cluster'] != -1] plot_clusters(dbs_interval_noiseless) dbs_interval_noiseless

def cluster_metrics(df): """ Generates the clusters based off of location and dwell time using DBScan Parameters: _ _ _ _ _ _ _ _ _ _ _ _ _ df1 = DataFrame, containing lat,long, dwell time interval eps = Float, epsilon value for DBScan min_sample = Integer, mininimum number of samples in a given cluster for DBScan Returns _ _ _ _ _ _ _ _ _ _ _ _ _ A DataFrame containing lat, long, dwell time interval, associated DBS Cluster """ ## Takes in DataFrame with dbs cluster labels and return summary metrics per Cluster cluster_hour_info = pd.DataFrame() # Grouped by labels to find the average dwell time in hours for each cluster cluster_hour_info["avg duration (hr)"] = (df.groupby("dbs_cluster")["dwell time duration (seconds)"].agg('mean'))/3600 # Grouped by labels to find the median dwell time in hours for each cluster cluster_hour_info["median duration (hr)"] = (df.groupby("dbs_cluster")["dwell time duration (seconds)"].agg('median'))/3600 # Grouped by labels to find the min and max dwell times in hours for each cluster cluster_hour_info["min duration (hr)"] = (df.groupby("dbs_cluster")["dwell time duration (seconds)"].agg('min'))/3600 cluster_hour_info["max duration (hr)"] = (df.groupby("dbs_cluster")["dwell time duration (seconds)"].agg('max'))/3600 cluster_hour_info["std duration (hr)"] = (df.groupby("dbs_cluster")["dwell time duration (seconds)"].agg('std'))/3600 cluster_hour_info["cluster size"] = df.groupby("dbs_cluster").size() cluster_hour_info.sort_values(by = 'std duration (hr)') return cluster_hour_info

def train_pred_model(dbs_df1, knn = True, svm = False): """ Trains a provided predictive model on dataset for cluster-prediction. Performs hyper-parameter tuning using a randomised Cross-Validaiton Search. Parameters: _ _ _ _ _ _ _ _ _ _ _ _ _ dbs_df1 = DataFrame, containing lat, long, associated DBS Cluster knn = Boolean, specifies if model is knn. If not knn or svm, then XGBoost svm = Boolean, specifies if model is svm. If not svm or knn, then XGBoost Returns _ _ _ _ _ _ _ _ _ _ _ _ _ A trained predictedive model with tuned hyperparameters to predict dbs_cluster, prints CV score. BEWARE MAY BE COMPUTATIONALLY EXPENSIVE DUE TO HYPER-PARAMETER TUNING, USE OF GPU ADVISED """ dbs_df = dbs_df1.copy() X = dbs_df.loc[: , ['latitude','longitude', 'daylight']] y = dbs_df.loc[:, 'dbs_cluster'] X_train, X_test, y_train, y_test = train_test_split(X,y, test_size=0.3, random_state=0) if svm: kf = KFold(4, shuffle=True, random_state=0) pipe = Pipeline([('scaler', StandardScaler()), ('svc', SVC())]) param_grid = { 'svc__kernel': ['poly', 'rbf'], 'svc__coef0' : np.logspace(-2,2,9), 'svc__degree': [1,2,3,4,5,6], 'svc__C': np.logspace(-2,2,9), 'svc__cache_size': [ 12000] } search = RandomizedSearchCV(pipe, param_grid, cv =kf, n_jobs = -1) best_model = search.fit(X_train, y_train) print("Best parameter (CV score=%0.3f):" % search.best_score_) print(search.best_params_) cv_scores = cross_val_score(best_knn,X_train,y_train, cv=3) print ("Cross-Validation Accuracies:", cv_scores) print ("Overall CV score is:", np.mean(cv_scores)) return best_model if knn: print(KNeighborsClassifier().get_params().keys()) N_grid = range(3,50) kf = KFold(4, shuffle=True, random_state=0) param_grid = { 'n_neighbors': N_grid, 'algorithm':['ball_tree', 'kd_tree','brute'], 'p':[1,2, 3,4, 5, 6, 7, 8, 9, 10, 15, 20], 'leaf_size':np.arange(10,80), 'weights': ['uniform', 'distance'], 'n_jobs': [-1]} knn = KNeighborsClassifier() search = GridSearchCV(knn, param_grid, cv = kf, n_jobs = -1) best_model = search.fit(X_train, y_train) best_knn = best_model.best_estimator_ cv_scores = cross_val_score(best_knn,X_train,y_train, cv=3) print ("Cross-Validation Accuracies:", cv_scores) print ("Overall CV score is:", np.mean(cv_scores)) return best_knn else: kf = KFold(4, shuffle=True, random_state=0) params = { 'min_child_weight': range(2,8), 'eta' : [.3, .2, .1, .05, .01, .005], 'subsample': [0.65, 0.8, 1.0], 'gamma': [0.5, 1, 1.5, 2, 3.5, 5], 'colsample_bytree': [0.6, 0.8, 1.0], 'max_depth': [3, 4, 5, 7, 10], 'njobs': [-1]} model = XGBClassifier() model.get_params().keys() search = GridSearchCV(model, params, cv = kf, n_jobs = -1) best_model = search.fit(X_train, y_train) best_xg = best_model.best_estimator_ cv_scores = cross_val_score(best_xg,X_train,y_train, cv=3) print ("Cross-Validation Accuracies:", cv_scores) print ("Overall CV score is:", np.mean(cv_scores)) return best_xg def predicted_cluster_vals(df1,best_model): """ Fits best_model and predicts the cluster of each datapoint in the dataframe using either KNN or XGBoost Parameters _ _ _ _ _ _ _ _ _ _ _ _ _ df1 = DataFrame, containing lat, long, daylight, dbs_cluster pred_model = KNN or XGBoost model Returns _ _ _ _ _ _ _ _ _ _ _ _ _ A dataframe containing Lat,Long, predicted_cluster """ X = df1.loc[: , ['latitude','longitude', 'daylight']] y = df1.loc[:, 'dbs_cluster'] X_train, X_test, y_train, y_test = train_test_split(X,y, test_size=0.3) best_model = best_model.fit(X_train, y_train) pred_df = df1.loc[: , ['latitude','longitude','daylight']] pred_df['predicted_cluster'] = best_model.predict(pred_df) return pred_df

##CHANGE TO GENERAL CASE FOR INTERVAL BINS (FROM 0-3, 3-6, +6) def get_interval_pred(pred, actual): """ Connects your prediction DataFrame with the associated cluster intervals to retrieve interval bin predictions Parameters: _ _ _ _ _ _ _ _ _ _ _ _ _ pred = DataFrame, containing lat,long, predicted DBS Cluster actual = DataFrame, containing cluster, lat, long, dwell time interval bins Returns _ _ _ _ _ _ _ _ _ _ _ _ _ A DataFrame containing lat, long, dwell time interval, associated interval bins (Predictions for interval bins) """ actual = dbs_interval_noiseless.drop_duplicates(subset ="dbs_cluster") actual = actual.copy() pred = pred.copy() merged_data = pd.merge(pred.reset_index(),actual,how='inner',left_on=['predicted_cluster'],right_on=['dbs_cluster']) merged_data = merged_data.set_index(["vin","sequence"]).sort_index() sortedvals = merged_data[["latitude_x","longitude_x","predicted_cluster","0-3hr","3-6hr","+6hr"]] return sortedvals.sort_index()

Best_KNN_Model = train_pred_model(dbs_interval_noiseless, knn = True) Best_KNN_Model

KNeighborsClassifier(algorithm='kd_tree', leaf_size=69, n_jobs=-1, n_neighbors=20, p=6, weights='distance')

#Add this to get_interval_pred pred = dbs_interval_noiseless[['latitude', 'longitude', 'daylight']] actual_clusters_with_intervals = dbs_interval_noiseless.drop_duplicates(subset ="dbs_cluster") dbs_interval_noiseless

pred_dataframe = predicted_cluster_vals(dbs_interval_noiseless, Best_KNN_Model)

def accuracy(pred_df,actual_df): """ Creates an accuracy measurement for the predicted cluster dwell times Parameters: _ _ _ _ _ _ _ _ _ _ _ _ _ pred_df = DataFrame, containing lat,long, predicted DBS Cluster, predicted interval bin actual_df = DataFrame, containing cluster, lat, long, dwell time interval bins, actual interval bin Returns _ _ _ _ _ _ _ _ _ _ _ _ _ A Float value of the accuracy of the predictions """ pred = pred_df[['0-3hr', '3-6hr','+6hr']].values actual = actual_df.loc[pred_df.index.values][['0-3hr', '3-6hr','+6hr']].values return metrics.accuracy_score(pred, actual) def interval_pred_accuracy(cluster_pred_df, actual): """ Generates interval predictions for predicted clusters and returns accuracy for interval predictions Parameters _ _ _ _ _ _ _ _ _ _ _ _ _ cluster_pred_df = DataFrame, containing lat, long, predicted DBS Cluster actual = DataFrame, containing lat, long, actual DBS Cluster Returns _ _ _ _ _ _ _ _ _ _ _ _ _ An accuracy score for interval predictions """ pred = get_interval_pred(cluster_pred_df, actual) return accuracy(pred, actual)

interval_pred_accuracy(pred_dataframe, dbs_interval_noiseless)

## Binary > 10 hours classification, no pre-clustering dbs_copy = dbs_data.copy() dbs_copy['more than 10 hr'] = 1 * 1* (dbs_data["dwell time duration (seconds)"]/3600 > 10) X = dbs_copy.loc[: , ['latitude','longitude']] y = dbs_copy.loc[:, ['more than 10 hr']] X_train, X_test, y_train, y_test = train_test_split(X,y, random_state=0) kf = KFold(4, shuffle=True, random_state=0) params = { 'min_child_weight': range(2,8), 'eta' : [.3, .2, .1, .05, .01, .005], 'subsample': [0.65, 0.8, 1.0], 'gamma': [0.5, 1, 1.5, 2, 3.5, 5], 'colsample_bytree': [0.6, 0.8, 1.0], 'max_depth': [3, 4, 5, 7, 10], 'n_jobs': [-1]} model = XGBClassifier() search = RandomizedSearchCV(model, params, cv = kf, n_jobs = -1) best_model = search.fit(X_train, y_train) best_xg = best_model.best_estimator_ cv_scores = cross_val_score(best_xg,X_train,y_train, cv=3) print ("Cross-Validation Accuracies:", cv_scores) print ("Overall CV score is:", np.mean(cv_scores)) metrics.accuracy_score(best_xg.predict(X_train), y_train)