data = pd.read_csv('data/bike_oneweekfrom20140505.csv', index_col=0, parse_dates=True)

data.info()

data.columns

"""
Explanation: Use the bash =)
End of explanation
"""


new_column_names = ['trip_time', 'pickup_datetime', 'dropoff_datetime', 'start_station_id',
       'start_station_name', 'pickup_latitude',
       'pickup_longitude', 'end_station_id', 'end_station_name',
       'dropoff_latitude', 'dropoff_longitude', 'bikeid', 'usertype',
       'birth year', 'gender']

data.columns = new_column_names

data.describe()

"""
Explanation: Rename the columns. They're now streamlined with the taxi input.
End of explanation
"""


data['pickup_datetime'] =pd.to_datetime(data['pickup_datetime'], format = '%Y-%m-%d %H:%M:%S')
data['dropoff_datetime'] =pd.to_datetime(data['dropoff_datetime'], format = '%Y-%m-%d %H:%M:%S')
data['trip_time'] = pd.to_timedelta(data['trip_time'], 's')

data.describe().transpose()

data.head()

data.info()

"""
Explanation: Parse the dates.
End of explanation
"""


group_by_start = data.groupby(data.start_station_id, )
usage_freq = group_by_start.trip_time.count()
print('Amount of stations: ' + str(len(usage_freq)))
plt.hist(usage_freq, bins = 50)
plt.title('Frequency of starts from start stations')

"""
Explanation: Have a look at every station. How many pickups are there per station?
End of explanation
"""


group_by_start_name = data.groupby(data.start_station_name, )
usage_freq_name = group_by_start_name.trip_time.count()
usage_freq_name.sort_values().keys()[-10:]

usage_freq_name.sort_values().keys()

"""
Explanation: When you take a look at the map, you won't be surprised that the Grand Central Station is the one with the most frequent usage.
End of explanation
"""


group_by_end_name = data.groupby(data.end_station_name, )
usage_freq_name = group_by_end_name.trip_time.count()
usage_freq_name.sort_values().keys()

"""
Explanation: So we have 326 different start stations.
End of explanation
"""


group_by_start_id = data.groupby(data.start_station_id, )
start_ids = group_by_start_id.trip_time.count().sort_values().keys()
group_by_end_id = data.groupby(data.end_station_id, )
end_ids = group_by_end_id.trip_time.count().sort_values().keys()
len(set(start_ids).intersection(end_ids))

"""
Explanation: We have 326 different end stations, too. But are all end_stations contained in start_stations and vive versa?
End of explanation
"""


data.isnull().sum()

"""
Explanation: So this is the proof, that all stations are used as start and end stations.
Check for missing and false data
End of explanation
"""


(data==0).sum()

"""
Explanation: So there is not that much data missing. That's quite surprising, maybe it's wrong.
End of explanation
"""


from collections import Counter
print(Counter(data.gender))
Counter(data['birth year'])

"""
Explanation: So we only have many zeros in the gender-feature. Gender can have 3 values (0=unknown, 1=male, 2=female)
End of explanation
"""


plt.hist((data['trip_time'] / np.timedelta64(1, 'm')), bins=30, range=[0, 100])
plt.title('Distribution of trip_time in minutes')
plt.xlabel('trip_time')
plt.ylabel('frequency')
plt.savefig('figures/bike_trip_time.png', format='png', dpi=300)

"""
Explanation: So we see, we have also missing values in "gender" and "birth year". But we are not interested in this features, therefore we can ignore this missing values at this point.
Quick pverview about the trip_times
End of explanation
"""


sns.boxplot((data['trip_time'] / np.timedelta64(1, 'm')))

"""
Explanation: We are interested in the trip time in minutes.
End of explanation
"""


data.trip_time.quantile([0.01, 0.05, 0.1, 0.25, 0.5, 0.75, 0.9, 0.95, 0.99])

len(data.trip_time.value_counts().values)

"""
Explanation: We have lots of outliers.
So as we can see, we have many outliers.
End of explanation
"""


anomaly = data.loc[(data['dropoff_longitude'].isnull()) | (data['dropoff_latitude'].isnull()) | 
                       (data['pickup_longitude'].isnull()) | (data['pickup_latitude'].isnull())]
data = data.drop(anomaly.index)

anomaly['flag'] = 'geo_NA'

data.isnull().sum()

"""
Explanation: Identify the the cases without geo data and remove them from our data to be processed.
End of explanation
"""


len(data)/(len(data)+len(anomaly))

anomaly.tail()

"""
Explanation: So how many percent of data are left to be processed?
End of explanation
"""


plt.hist(data.trip_time.values / np.timedelta64(1, 'm'), bins=50, range=[0,100])

"""
Explanation: <font color = 'blue'> So we dropped nothing because of missing geo tags. </font color>
End of explanation
"""


data['trip_dist'] = -1 # Init trip_dist

"""
Explanation: We sometimes have some unreasonably small trip_times.
Calculate the distance of a trip. Use the difference of lon / lat data and parse it to metrics.
End of explanation
"""


# inpout for vincenty:(location.latitude, location.longitude)
from geopy.distance import vincenty
for i in range(0, (len(data)-1)):
    pickup = (data.iloc[i]['pickup_latitude'], data.iloc[i]['pickup_longitude'])
    dropoff = (data.iloc[i]['dropoff_latitude'], data.iloc[i]['dropoff_longitude'])
    data.set_value(i, 'trip_dist', vincenty(pickup,dropoff).meters)

data.trip_dist

"""
Explanation: Be aware, this operation is quite slow.
End of explanation
"""


data.columns

"""
Explanation: We cannot use miles because vincenty ceils / floors the results
End of explanation
"""


# data.to_csv('data/bike_20140505_with_dist.csv')

data['avg_velocity'] = data.trip_dist.values/(data.trip_time / (np.timedelta64(1, 'h')))

data.avg_velocity

"""
Explanation: Check if new column is present
End of explanation
"""


data.avg_velocity = data.avg_velocity/1000

data.avg_velocity

"""
Explanation: Parse to km/h
End of explanation
"""


data['avg_velocity'] = data.avg_velocity*0.627371

# data.to_csv('data/bike_20140505_with_dist_and_avg_velo.csv')

plt.hist(data.avg_velocity, bins=60, range=[0,15])
plt.title('Distribution of avg_velocity in mph')
plt.xlabel('avg_velocity')
plt.ylabel('frequency')
plt.savefig('figures/bike_avg_vel.png', format='png', dpi=300)

data.avg_velocity.describe() # in mph

data.head()

data.avg_velocity.quantile([.0001,.01, .5, .75, .95, .975, .99, .995])

"""
Explanation: Convert avg_velocity from km/h to miles/h. (1km = 0,621371 miles)
End of explanation
"""


lb = 5
ub = 15
anomaly = data.loc[(data['avg_velocity'] < lb) | (data['avg_velocity'] > ub)] 

# be careful! Maybe adjust to append.
#anomaly.loc[
anomaly.loc[data.loc[(data['avg_velocity'] < lb)].index,'flag'] = 'too_slow'
anomaly.loc[data.loc[(data['avg_velocity'] > ub)].index,'flag'] = 'too_fast'
data = data.drop(anomaly.index, errors='ignore') # ignore uncontained labels / indices
print(1-len(data)/(len(data)+len(anomaly)))

anomaly.head()

"""
Explanation: Be aware: We calculated the avg_velocity based on the vincenty distance (equals direct connection). A real driver would never drive that fast. We take all avg_velos >= 5mph (approx. >5%-quantile) into account. We use an upper bound, too. It might not be possible to drive "too fast". If someone drives too fast, the technology has had a malfunction.
End of explanation
"""


data.avg_velocity.describe()

anomaly.tail()

"""
Explanation: We dropped about 42% of the data!
End of explanation
"""


x = data.pickup_longitude
y = data.pickup_latitude

bins = 100;
H, xedges, yedges = np.histogram2d(x, y, bins=bins)

plt.jet()
fig = plt.figure(figsize=(20, 10))
#ax.set_title('pcolormesh: exact bin edges')
#mesh = plt.pcolormesh(X, Y, H)
plt.hist2d(x,y,bins=bins)
plt.colorbar()
plt.scatter(xedges[79], yedges[61], marker='x')
ax = fig.gca()
ax.grid(False)
plt.xlabel('Longitude')
plt.ylabel('Latitude')
plt.title('Pickup Frequency')
#ax.set_aspect('equal')
#plt.savefig('figure.pdf', format='pdf')
#plt.savefig('figure.png', format='png')
plt.savefig('figures/bike_stations.png', format='png', dpi=300)

"""
Explanation: A little drawing
End of explanation
"""


np.amax(H)

np.where(H==2001)

"""
Explanation: Where is a pickup maximum?
End of explanation
"""


yedges[79],xedges[61], 

type(H)

print(H.shape)
print(H.size)
print(H.max())

print('Current bin width:')
print(vincenty((xedges[0], yedges[0]), (xedges[1], yedges[1])).meters)

distances = (-1)*np.ones(len(xedges)-1)
for x in range(0,len(xedges)-1, 1):
    distances[x] = vincenty((xedges[x], yedges[x]), (xedges[x+1], yedges[x+1])).meters

print('Sizes of all bins in meters:')
distances

(H==0).sum()/H.size

"""
Explanation: The area with the maximum pickups is around the Grand Central Terminal. Not very surprising.
End of explanation
"""


jfk_geodata = (40.641547, -73.778118)
ridgefield_geodata = (40.856406, -74.020642)
data_in_box = data.loc[(data['dropoff_latitude'] > jfk_geodata[0]) & 
                       (data['dropoff_longitude'] < jfk_geodata[1]) &
                       (data['dropoff_latitude'] < ridgefield_geodata[0]) & 
                       (data['dropoff_longitude'] > ridgefield_geodata[1]) & 
                       (data['pickup_latitude'] > jfk_geodata[0]) & 
                       (data['pickup_longitude'] < jfk_geodata[1]) &
                       (data['pickup_latitude'] < ridgefield_geodata[0]) & 
                       (data['pickup_longitude'] > ridgefield_geodata[1])         
                       ]
# taxidata = taxidata.drop(anomaly.index)

len(data_in_box)/len(data)

"""
Explanation: So we know that about 97% of the bins have 0 pickups in it. This was expected, because we have static stations only.
Only look at trips in a given bounding box
End of explanation
"""


data_in_box.head()

"""
Explanation: We lost no trips through the bounding box.
End of explanation
"""


h = data_in_box.avg_velocity.values
plt.figure(figsize=(20,10))
plt.hist(h, normed=False,  bins=150)
    #, histtype='stepfilled')
#plt.yscale('log')
#plt.ylabel('log(freq x)', fontsize=40)
#plt.xlabel('x = avg_amount_per_minute', fontsize=40)
#print('Min:' +  str(min(h)) + '\nMax:' +  str(max(h)))
#plt.locator_params(axis = 'x', nbins = 20)
plt.show()

data_in_box.head()

"""
Explanation: Let's take a first look at the distribution of the cleaned target variable which we want to estimate:
End of explanation
"""


time_regression_df = pd.DataFrame([ #data_in_box['pickup_datetime'].dt.day, # leave this out
                          data_in_box['pickup_datetime'].dt.dayofweek,
                          data_in_box['pickup_datetime'].dt.hour,
                          data_in_box['pickup_latitude'],
                          data_in_box['pickup_longitude'],
                          data_in_box['dropoff_latitude'],
                          data_in_box['dropoff_longitude'],
                          np.ceil(data_in_box['avg_velocity'])
                         ]).T

time_regression_df.columns = [#'pickup_datetime_day', 
                              'pickup_datetime_dayofweek', 'pickup_datetime_hour',
                                 'pickup_latitude', 'pickup_longitude', 'dropoff_latitude', 'dropoff_longitude',
                                 'avg_velocity_mph']

"""
Explanation: Make a new dataframe with features and targets
End of explanation
"""


time_regression_df.head()

time_regression_df.tail()

time_regression_df.ix[:,0:7].describe()

print(time_regression_df.avg_velocity_mph.value_counts())
print(len(time_regression_df.avg_velocity_mph.value_counts()))

"""
Explanation: Use minutes for prediction instead of seconds (ceil the time). Definitley more robust than seconds!
End of explanation
"""


from sklearn import cross_validation as cv
y = time_regression_df['avg_velocity_mph']
X = time_regression_df.ix[:,0:6]
X_train, X_test, y_train, y_test = cv.train_test_split(X, y,test_size=0.05,random_state=0) 

from sklearn import cross_validation as cv
time_regression_df_train, time_regression_df_test = cv.train_test_split(time_regression_df, test_size=0.05, random_state=99)
y_train = time_regression_df_train['avg_velocity_mph']
x_train = time_regression_df_train.ix[:, 0:6]
y_test = time_regression_df_test['avg_velocity_mph']
x_test = time_regression_df_test.ix[:, 0:6]

"""
Explanation: So we hace 10 different velocities to predict.
Split the data into a training dataset and a test dataset. Evaluate the performance of the decision tree on the test data
End of explanation
"""


x_test.head()

y_test.head()

xy_test = pd.concat([x_test, y_test], axis=1)

xy_test.head()

"""
Explanation: Changed test size to 5%
End of explanation
"""


#Xy_test.to_csv('taxi_tree_test_Xy_20130506-12.csv')
#X_test.to_csv('taxi_tree_test_X_20130506-12.csv')
#y_test.to_csv('taxi_tree_test_y_20130506-12.csv')
# Xy_test.to_csv('bike_tree_test_Xy_20140505-11.csv')
# X_test.to_csv('bike_tree_test_X_20140505-11.csv')
# y_test.to_csv('bike_tree_test_y_20140506-11.csv')

# Xy_test_sample = Xy_test.sample(10000, random_state=99)

# Xy_test_sample.to_csv('taxi_tree_test_Xy_sample.csv')

# Xy_test_sample.head()

"""
Explanation: If you want to export something...
End of explanation
"""


print(x_train.shape)
print(x_train.size)
print(x_test.shape)
print(X.shape)
print(x_train.shape[0]+x_test.shape[0])

import time
# Import the necessary modules and libraries
from sklearn.tree import DecisionTreeRegressor
import numpy as np
import matplotlib.pyplot as plt

regtree = DecisionTreeRegressor(min_samples_split=100, random_state=99, max_depth=20)# formerly 15. 15 is reasonable
                                                                                     # random states: 99
regtree.fit(x_train, y_train)

regtree.score(x_test, y_test)


y_pred = regtree.predict(x_test)
diff_regtree = (y_pred-y_test)
# plt.figure(figsize=(12,10)) # not needed. set values globally
plt.hist(diff_regtree.values, bins=100, range=[-6,6])
print('Perzentile(%): ', [1,5,10,15,25,50,75,90,95,99], '\n', np.percentile(diff_regtree.values, [1,5,10,15,25,50,75,90,95,99]))
print('Absolute time deviation (in 1k): ', sum(abs(diff_regtree))/1000)
plt.title('Simple Decision Tree Regressor')
plt.xlabel('deviation in mph')
plt.ylabel('frequency')
plt.savefig('figures/bike_tree.png', format='png', dpi=300)

diff_regtree.describe()

"""
Explanation: Just to be sure
End of explanation
"""


from sklearn.ensemble import RandomForestRegressor
n_est_list = list(range(2,32,8))
min_sam_leaf_list = [2,50,100,150]
max_depth_list = list(range(2,32,8))
results = np.empty([0,4])
for n_est in n_est_list:
    for min_sam_leaf in min_sam_leaf_list:
        for max_depth in max_depth_list:
            rd_regtree = RandomForestRegressor(n_estimators=n_est,n_jobs=6,min_samples_leaf=min_sam_leaf, random_state=99, max_depth=max_depth)
            rd_regtree.fit(x_train, y_train)
            score = rd_regtree.score(x_test, y_test)
            results = np.vstack((results, [n_est, min_sam_leaf, max_depth,score]))

best = np.where(results == max(results[:,3]))[0]
results[best,:]

"""
Explanation: So this is a bad result. Let's try random forest now.
Train several random forests and look at their scores
End of explanation
"""


results[np.argsort(results[:, 3])][-10:,:]

results = np.vstack((results, [n_est, min_sam_leaf, max_depth,score]))
#results = np.vstack((results, [n_est, min_sam_leaf, max_depth,score]))
#results = np.vstack((results, [n_est, min_sam_leaf, max_depth,score]))
results[:,3]

"""
Explanation: Watch at the best results.
End of explanation
"""


from sklearn.ensemble import RandomForestRegressor

#rd_regtree = RandomForestRegressor(n_estimators=10,n_jobs=6,min_samples_leaf=4, random_state=99, max_depth=20)
rd_regtree = RandomForestRegressor(n_estimators=26,n_jobs=6,min_samples_leaf=2, random_state=99, max_depth=26)
#total sum of diff: 1132
#rd_regtree = RandomForestRegressor(n_estimators=40,n_jobs=-1,min_samples_split=3, random_state=99, max_depth=11)
#total sum of diff: 1129
rd_regtree.fit(x_train, y_train)
print('R²: ', rd_regtree.score(x_test, y_test))

rd_regtree.feature_importances_

tree_list = rd_regtree.estimators_
for i in range(0,len(tree_list)):
    print((tree_list[i].feature_importances_))


from mpl_toolkits.mplot3d import Axes3D
fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
ax.scatter(results[:,0],results[:,2],results[:,3], marker='o')
plt.show()

"""
Explanation: Train the best random forest again seperately
End of explanation
"""


y_pred = rd_regtree.predict(x_test)
np.linalg.norm(np.ceil(y_pred)-y_test)
diff_rd = (y_pred-y_test)
# plt.figure(figsize=(12,10)) # not needed. set values globally
plt.hist(diff_rd.values, bins=100, range=[-6,6])
print('Perzentile(%): ', [1,5,10,15,25,50,75,90,95,99], '\n', np.percentile(diff_rd.values, [1,5,10,15,25,50,75,90,95,99]))
print('Absolute time deviation (in 1k): ', sum(abs(diff_rd))/1000)
plt.title('Random Forest Regressor')
plt.xlabel('deviation in mph')
plt.ylabel('frequency')
plt.savefig('figures/bike_randomforest.png', format='png', dpi=150)

diff_rd.describe()

"""
Explanation: So this plot is not particularly useful in this case. But it's a hinto to how to visualize the results.
End of explanation
"""


from sklearn.externals import joblib
joblib.dump(rd_regtree, 'randforlib/bike_regtree_26x_depth_26_mss_2_PY27.pkl', protocol=2)

! cd randforlib && zip bike_regtree_26x_depth_26_mss_2_PY27.zip bike_regtree_26x_depth_26_mss_2_PY27.pkl*

"""
Explanation: Dump the best random forest...
End of explanation
"""


