# ---
# jupyter:
#   jupytext:
#     text_representation:
#       extension: .py
#       format_name: light
#       format_version: '1.5'
#       jupytext_version: 1.15.2
#   kernelspec:
#     display_name: Python 3
#     language: python
#     name: python3
# ---

# # Linear Regression [35pts (+5 bonus)]
#
# ## Introduction
# One of the most widespread regression tools is the simple but powerful linear regression. In this notebook, you will engineer the Pittsburgh bus data into numerical features and use them to predict the number of minutes until the bus reaches the bus stop at Forbes and Morewood. 
#
# Notebook restriction: you may not use scikit-learn for this notebook.  
#
# ## Q1: Labeling the Dataset [8pts]
#
# You may have noticed that the Pittsburgh bus data has a predictions table with the TrueTime predictions on arrival time, however it does not have the true label: the actual number of minutes until a bus reaches Forbes and Morewood. You will have to generate this yourself. 
#
# Using the `all_trips` function that you implemented in homework 2, you can split the dataframe into separate trips. You will first process each trip into a form more natural for the regression setting. For each trip, you will need to locate the point at which a bus passes the bus stop to get the time at which the bus passes the bus stop. From here, you can calculate the true label for all prior datapoints, and throw out the rest. 
#
# ### Importing functions from homework 2
#
# Using the menu in Jupyter, you can import code from your notebook as a Python script using the following steps: 
# 1. Click File -> Download as -> Python (.py)
# 2. Save file (time_series.py) in the same directory as this notebook 
# 3. (optional) Remove all test code (i.e. lines between AUTOLAB_IGNORE macros) from the script for faster loading time
# 4. Import from the notebook with `from time_series import function_name`
#
# ### Specifications
#
# 1. To determine when the bus passes Morewood, we will use the Euclidean distance as a metric to determine how close the bus is to the bus stop. 
# 2. We will assume that the row entry with the smallest Euclidean distance to the bus stop is when the bus reaches the bus stop, and that you should truncate all rows that occur **after** this entry.  In the case where there are multiple entries with the exact same minimal distance, you should just consider the first one that occurs in the trip (so truncate everything after the first occurance of minimal distance). 
# 3. Assume that the row with the smallest Euclidean distance to the bus stop is also the true time at which the bus passes the bus stop. Using this, create a new column called `eta` that contains for each row, the number of minutes until the bus passes the bus stop (so the last row of every trip will have an `eta` of 0).
# 4. Make sure your `eta` is numerical and not a python timedelta object. 

import pandas as pd
import numpy as np
import scipy.linalg as la
from collections import Counter
import datetime

# AUTOLAB_IGNORE_START
from time_series import load_data, split_trips
vdf, _ = load_data('bus_train.db')
all_trips = split_trips(vdf)
# AUTOLAB_IGNORE_STOP

#https://stackoverflow.com/questions/4983258/python-how-to-check-list-monotonicity
pd.options.mode.chained_assignment = None
def label_and_truncate(trip, bus_stop_coordinates):
    """ Given a dataframe of a trip following the specification in the previous homework assignment,
        generate the labels and throw away irrelevant rows. 
        
        Args: 
            trip (dataframe): a dataframe from the list outputted by split_trips from homework 2
            stop_coordinates ((float, float)): a pair of floats indicating the (latitude, longitude) 
                                               coordinates of the target bus stop. 
            
        Return:
            (dataframe): a labeled trip that is truncated at Forbes and Morewood and contains a new column 
                         called `eta` which contains the number of minutes until it reaches the bus stop. 
        """
    result_lst =pd.DataFrame()
    lat=(trip["lat"]-bus_stop_coordinates[0])**2
    lon =(trip["lon"]-bus_stop_coordinates[1])**2
    dist = np.array(lat+lon)
    result_lst=trip[:np.argmin(dist)+1]
    passing_time = np.array(result_lst.tail(1).index)
    ongoing_time = np.array(result_lst.index)
    etas = pd.to_numeric((passing_time-ongoing_time).astype('timedelta64[m]'))
    result_lst["eta"]=etas
    return result_lst
    pass
#trip.index[:np.argmin(dist)+1]
# AUTOLAB_IGNORE_START
morewood_coordinates = (40.444671114203, -79.94356058465502) # (lat, lon)
labeled_trips = [label_and_truncate(trip, morewood_coordinates) for trip in all_trips]
# print(len(labeled_trips))
labeled_vdf = pd.concat(labeled_trips).reset_index()
# We remove datapoints that make no sense (ETA more than 10 hours)
labeled_vdf = labeled_vdf[labeled_vdf["eta"] < 10*60].reset_index(drop=True)
print(Counter([len(t) for t in labeled_trips]))
print(labeled_vdf.head())
# AUTOLAB_IGNORE_STOP

# For our implementation, this returns the following output
# ```python
# >>> Counter([len(t) for t in labeled_trips])
# Counter({1: 506, 21: 200, 18: 190, 20: 184, 19: 163, 16: 162, 22: 159, 17: 151, 23: 139, 31: 132, 15: 128, 2: 125, 34: 112, 32: 111, 33: 101, 28: 98, 14: 97, 30: 95, 35: 95, 29: 93, 24: 90, 25: 89, 37: 86, 27: 83, 39: 83, 38: 82, 36: 77, 26: 75, 40: 70, 13: 62, 41: 53, 44: 52, 42: 47, 6: 44, 5: 39, 12: 39, 46: 39, 7: 38, 3: 36, 45: 33, 47: 33, 43: 31, 48: 27, 4: 26, 49: 26, 11: 25, 50: 25, 10: 23, 51: 23, 8: 19, 9: 18, 53: 16, 54: 15, 52: 14, 55: 14, 56: 8, 57: 3, 58: 3, 59: 3, 60: 3, 61: 1, 62: 1, 67: 1}) 
# >>> labeled_vdf.head()
#                tmstmp   vid        lat        lon  hdg   pid   rt        des  \
# 0 2016-08-11 10:56:00  5549  40.439504 -79.996981  114  4521  61A  Swissvale   
# 1 2016-08-11 10:57:00  5549  40.439504 -79.996981  114  4521  61A  Swissvale   
# 2 2016-08-11 10:58:00  5549  40.438842 -79.994733  124  4521  61A  Swissvale   
# 3 2016-08-11 10:59:00  5549  40.437938 -79.991213   94  4521  61A  Swissvale   
# 4 2016-08-11 10:59:00  5549  40.437938 -79.991213   94  4521  61A  Swissvale   
#
#    pdist  spd tablockid  tatripid  eta  
# 0   1106    0  061A-164      6691   16  
# 1   1106    0  061A-164      6691   15  
# 2   1778    8  061A-164      6691   14  
# 3   2934    7  061A-164      6691   13  
# 4   2934    7  061A-164      6691   13 
# ```

# ## Q2: Generating Basic Features [8pts]
# In order to perform linear regression, we need to have numerical features. However, not everything in the bus database is a number, and not all of the numbers even make sense as numerical features. If you use the data as is, it is highly unlikely that you'll achieve anything meaningful.
#
# Consequently, you will perform some basic feature engineering. Feature engineering is extracting "features" or statistics from your data, and hopefully improve the performance if your learning algorithm (in this case, linear regression). Good features can often make up for poor model selection and improve your overall predictive ability on unseen data. In essence, you want to turn your data into something your algorithm understands. 
#
# ### Specifications
# 1. The input to your function will be a concatenation of the trip dataframes generated in Q1 with the index dropped (so same structure as the original dataframe, but with an extra column and less rows). 
# 2. Linear models typically have a constant bias term. We will encode this as a column of 1s in the dataframe. Call this column 'bias'. 
# 2. We will keep the following columns as is, since they are already numerical:  pdist, spd, lat, lon, and eta 
# 3. Time is a cyclic variable. To encode this as a numerical feature, we can use a sine/cosine transformation. Suppose we have a feature of value f that ranges from 0 to N. Then, the sine and cosine transformation would be $\sin\left(2\pi \frac{f}{N}\right)$ and $\cos\left(2\pi \frac{f}{N}\right)$. For example, the sine transformation of 6 hours would be $\sin\left(2\pi \frac{6}{24}\right)$, since there are 24 hours in a cycle. You should create sine/cosine features for the following:
#     * day of week (cycles every week, 0=Monday)
#     * hour of day (cycles every 24 hours, 0=midnight)
#     * time of day represented by total number of minutes elapsed in the day (cycles every 60*24 minutes, 0=midnight).
# 4. Heading is also a cyclic variable, as it is the ordinal direction in degrees (so cycles every 360 degrees). 
# 4. Buses run on different schedules on the weekday as opposed to the weekend. Create a binary indicator feature `weekday` that is 1 if the day is a weekday, and 0 otherwise. 
# 5. Route and destination are both categorical variables. We can encode these as indicator vectors, where each column represents a possible category and a 1 in the column indicates that the row belongs to that category. This is also known as a one hot encoding. Make a set of indicator features for the route, and another set of indicator features for the destination. 
# 6. The names of your indicator columns for your categorical variables should be exactly the value of the categorical variable. The pandas function `pd.DataFrame.get_dummies` will be useful. 

# +
def create_features(vdf):
    """ Given a dataframe of labeled and truncated bus data, generate features for linear regression. 
    
        Args:
            df (dataframe) : dataframe of bus data with the eta column and truncated rows
        Return: 
            (dataframe) : dataframe of features for each example
        """
    df = pd.DataFrame()
    df['pdist'] =vdf["pdist"] 
    df['spd']= vdf['spd']
    df['lat']= vdf['lat']
    df['lon']= vdf['lon']
    df['eta']= vdf["eta"] 
    df['sin_hdg']=np.sin(2*np.pi*vdf["hdg"]/360)
    df['cos_hdg']=np.cos(2*np.pi*vdf["hdg"]/360)
    df['sin_day_of_week']=np.sin((2*np.pi*vdf["tmstmp"].dt.dayofweek)/7)
    df['cos_day_of_week']=np.cos((2*np.pi*vdf["tmstmp"].dt.dayofweek)/7)
    df['sin_hour_of_day']=np.sin((2*np.pi*vdf["tmstmp"].dt.hour)/24)
    df['cos_hour_of_day']=np.cos((2*np.pi*vdf["tmstmp"].dt.hour)/24)
    minutes=pd.DataFrame()
    mins=[]
    for i in vdf["tmstmp"]:
        d1 = datetime.datetime.combine(i,datetime.datetime.min.time())
        secs=(i-d1).total_seconds()/60
        mins.append(secs)
    minutes["mins"]=mins
    df['sin_time_of_day']=np.sin((2*np.pi*minutes["mins"])/(60*24))
    df['cos_time_of_day']=np.cos((2*np.pi*minutes["mins"])/(60*24))
    df["weekday"]=[1 if i<4 else 0 for i in vdf["tmstmp"].dt.dayofweek]
    df['bias']=1
    set_des=set(vdf["des"])
    set_rt =set(vdf["rt"])
    for i in set_des:
        df[i]=[1 if i==j else 0 for j in vdf["des"]]
    for i in set_rt:
        df[i]=[1 if i==j else 0 for j in vdf["rt"]]
    return df
    
        
# AUTOLAB_IGNORE_START
#print(labeled_vdf["des"])
vdf_features = create_features(labeled_vdf)
vdf_features
# AUTOLAB_IGNORE_STOP
# -





# AUTOLAB_IGNORE_START
with pd.option_context('display.max_columns', 26):
    print(vdf_features.columns)
    print(vdf_features.head())
# AUTOLAB_IGNORE_STOP

# Our implementation has the following output. Verify that your code has the following columns (order doesn't matter): 
# ```python
# >>> vdf_features.columns
# Index([             u'bias',             u'pdist',               u'spd',
#                      u'lat',               u'lon',               u'eta',
#                  u'sin_hdg',           u'cos_hdg',   u'sin_day_of_week',
#          u'cos_day_of_week',   u'sin_hour_of_day',   u'cos_hour_of_day',
#          u'sin_time_of_day',   u'cos_time_of_day',           u'weekday',
#                u'Braddock ',          u'Downtown',   u'Greenfield Only',
#              u'McKeesport ', u'Murray-Waterfront',         u'Swissvale',
#                      u'61A',               u'61B',               u'61C',
#                      u'61D'],
#       dtype='object')
#    bias  pdist  spd        lat        lon  eta   sin_hdg   cos_hdg  \
# 0   1.0   1106    0  40.439504 -79.996981   16  0.913545 -0.406737   
# 1   1.0   1106    0  40.439504 -79.996981   15  0.913545 -0.406737   
# 2   1.0   1778    8  40.438842 -79.994733   14  0.829038 -0.559193   
# 3   1.0   2934    7  40.437938 -79.991213   13  0.997564 -0.069756   
# 4   1.0   2934    7  40.437938 -79.991213   13  0.997564 -0.069756   
#
#    sin_day_of_week  cos_day_of_week ...   Braddock   Downtown  \
# 0         0.433884        -0.900969 ...         0.0       0.0   
# 1         0.433884        -0.900969 ...         0.0       0.0   
# 2         0.433884        -0.900969 ...         0.0       0.0   
# 3         0.433884        -0.900969 ...         0.0       0.0   
# 4         0.433884        -0.900969 ...         0.0       0.0   
#
#    Greenfield Only  McKeesport   Murray-Waterfront  Swissvale  61A  61B  61C  \
# 0              0.0          0.0                0.0        1.0  1.0  0.0  0.0   
# 1              0.0          0.0                0.0        1.0  1.0  0.0  0.0   
# 2              0.0          0.0                0.0        1.0  1.0  0.0  0.0   
# 3              0.0          0.0                0.0        1.0  1.0  0.0  0.0   
# 4              0.0          0.0                0.0        1.0  1.0  0.0  0.0   
#
#    61D  
# 0  0.0  
# 1  0.0  
# 2  0.0  
# 3  0.0  
# 4  0.0  
#
# [5 rows x 25 columns]
# ```

# ## Q3 Linear Regression using Ordinary Least Squares [10 + 4pts]
# Now you will finally implement a linear regression. As a reminder, linear regression models the data as
#
# $$\mathbf y = \mathbf X\mathbf \beta + \mathbf \epsilon$$
#
# where $\mathbf y$ is a vector of outputs, $\mathbf X$ is also known as the design matrix, $\mathbf \beta$ is a vector of parameters, and $\mathbf \epsilon$ is noise. We will be estimating $\mathbf \beta$ using Ordinary Least Squares, and we recommending following the matrix notation for this problem (https://en.wikipedia.org/wiki/Ordinary_least_squares). 
#
# ### Specification
# 1. We use the numpy term array-like to refer to array like types that numpy can operate on (like Pandas DataFrames). 
# 1. Regress the output (eta) on all other features
# 2. Return the predicted output for the inputs in X_test
# 3. Calculating the inverse $(X^TX)^{-1}$ is unstable and prone to numerical inaccuracies. Furthermore, the assumptions of Ordinary Least Squares require it to be positive definite and invertible, which may not be true if you have redundant features. Thus, you should instead use $(X^TX + \lambda*I)^{-1}$ for identity matrix $I$ and $\lambda = 10^{-4}$, which for now acts as a numerical "hack" to ensure this is always invertible. Furthermore, instead of computing the direct inverse, you should utilize the Cholesky decomposition which is much more stable when solving linear systems. 

class LR_model():
    """ Perform linear regression and predict the output on unseen examples. 
        Attributes: 
            beta (array_like) : vector containing parameters for the features """
    
    def __init__(self, X, y):
        """ Initialize the linear regression model by computing the estimate of the weights parameter
            Args: 
                X (array-like) : feature matrix of training data where each row corresponds to an example
                y (array like) : vector of training data outputs 
            """
        self.beta = np.zeros(X.shape[1])
        x = np.array(X)
        y = np.array(y)
        lambdaa = 10**-4
        part1 = ((x.T @ x) +(lambdaa * np.identity(X.shape[1])))
        part2 =(x.T @ y)
        self.beta=np.linalg.solve(part1,part2)
        pass
        
    def predict(self, X_p): 
        """ Predict the output of X_p using this linear model. 
            Args: 
                X_p (array_like) feature matrix of predictive data where each row corresponds to an example
            Return: 
                (array_like) vector of predicted outputs for the X_p
            """
        x_arr = np.array(X_p)
        y_pred =  x_arr @ self.beta
        return y_pred
        pass



# We have provided some validation data for you, which is another scrape of the Pittsburgh bus data (but for a different time span). You will need to do the same processing to generate labels and features to your validation dataset. Calculate the mean squared error of the output of your linear regression on both this dataset and the original training dataset. 
#
# How does it perform? One simple baseline is to make sure that it at least predicts as well as predicting the mean of what you have seen so far. Does it do better than predicting the mean? Compare the mean squared error of a predictor that predicts the mean vs your linear classifier. 
#
# ### Specifications
# 1. Build your linear model using only the training data
# 2. Compute the mean squared error of the predictions on both the training and validation data. 
# 3. Compute the mean squared error of predicting the mean of the **training outputs** for all inputs. 
# 4. You will need to process the validation dataset in the same way you processed the training dataset.
# 5. You will need to split your features from your output (eta) prior to calling compute_mse

# +
# Calculate mean squared error on both the training and validation set
def compute_mse(LR, X, y, X_v, y_v):
    """ Given a linear regression model, calculate the mean squared error for the 
        training dataset, the validation dataset, and for a mean prediction
        Args:
            LR (LR_model) : Linear model
            X (array-like) : feature matrix of training data where each row corresponds to an example
            y (array like) : vector of training data outputs 
            X_v (array-like) : feature matrix of validation data where each row corresponds to an example
            y_v (array like) : vector of validation data outputs 
        Return: 
            (train_mse, train_mean_mse, 
             valid_mse, valid_mean_mse) : a 4-tuple of mean squared errors
                                             1. MSE of linear regression on the training set
                                             2. MSE of predicting the mean on the training set
                                             3. MSE of linear regression on the validation set
                                             4. MSE of predicting the mean on the validation set
    """
    yhat = LR.predict(X)
    
    mse_lr_tr = np.mean((y-yhat)**2)
    mse_me_tr = np.mean((y-np.mean(y))**2)
    
    yhat_v = LR.predict(X_v)
    
    mse_lr_v = np.mean((y_v-yhat_v)**2)
    mse_me_v = np.mean((y_v-np.mean(y))**2)
    
    return (mse_lr_tr,mse_me_tr,mse_lr_v,mse_me_v)


    

    pass



# +
# AUTOLAB_IGNORE_START
# First you should replicate the same processing pipeline as we did to the training set
vdf_valid, pdf_valid = load_data('bus_valid.db')
all_trips_valid =split_trips(vdf_valid)
labeled_trips_valid = [label_and_truncate(trip, morewood_coordinates) for trip in all_trips_valid]
labeled_vdf_valid = pd.concat(labeled_trips_valid).reset_index()
vdf_features_valid = create_features(labeled_vdf_valid)
# Separate the features from the output and pass it into your linear regression model.
y_df =vdf_features.eta
X_df = vdf_features.drop("eta",axis=1)
y_valid_df = vdf_features_valid.eta
X_valid_df =vdf_features_valid.drop("eta",axis=1)
LR = LR_model(X_df, y_df)
print(compute_mse(LR, 
                  X_df, 
                  y_df, 
                  X_valid_df, 
                  y_valid_df))

# AUTOLAB_IGNORE_STOP
# -

# As a quick check, our training data MSE is approximately 38.99. 

# ## Q4 TrueTime Predictions [5pts]
# How do you fare against the Pittsburgh Truetime predictions? In this last problem, you will match predictions to their corresponding vehicles to build a dataset that is labeled by TrueTime. Remember that we only evaluate performance on the validation set (never the training set). How did you do?
#
# ### Specification
# 1. You should use the pd.DataFrame.merge function to combine your vehicle dataframe and predictions dataframe into a single dataframe. You should drop any rows that have no predictions (see the how parameter). (http://pandas.pydata.org/pandas-docs/stable/generated/pandas.DataFrame.merge.html)
# 2. You can compute the TrueTime ETA by taking their predicted arrival time and subtracting the timestamp, and converting that into an integer representing the number of minutes. 
# 3. Compute the mean squared error for linear regression only on the rows that have predictions (so only the rows that remain after the merge). 

# +
def compare_truetime(LR, labeled_vdf, pdf):
    """ Compute the mse of the truetime predictions and the linear regression mse on entries that have predictions.
        Args:
            LR (LR_model) : an already trained linear model
            labeled_vdf (pd.DataFrame): a dataframe of the truncated and labeled bus data (same as the input to create_features)
            pdf (pd.DataFrame): a dataframe of TrueTime predictions
        Return: 
            (tt_mse, lr_mse): a tuple of the TrueTime MSE, and the linear regression MSE
        """
    featured_vdf = create_features(labeled_vdf)
    eta =featured_vdf["eta"]
    x = featured_vdf.drop("eta",axis=1)
    eta_hat = LR.predict(x)
    labeled_vdf["eta_lr"]=eta_hat
    labeled_vdf.reset_index()
    merged_df=pd.merge(labeled_vdf, pdf,how="inner")
    mins = np.array(((merged_df["prdtm"]-merged_df["tmstmp"]).dt.seconds)/60)
    merged_df["eta_tt"] =mins
    mse_lr = np.mean((merged_df["eta_lr"]-merged_df["eta"])**2)
    mse_tt = np.mean((merged_df["eta_tt"]-merged_df["eta"])**2)
    return (mse_tt,mse_lr)
    pass
    
# AUTOLAB_IGNORE_START
compare_truetime(LR, labeled_vdf_valid, pdf_valid)
# AUTOLAB_IGNORE_STOP
#50.20239900730732, 60.40782041336532
# -

# As a sanity check, your linear regression MSE should be approximately 50.20. 

# ## Q5 Feature Engineering contest (bonus)
#
# You may be wondering "why did we pick the above features?" Some of the above features may be entirely useless, or you may have ideas on how to construct better features. Sometimes, choosing good features can be the entirety of a data science problem. 
#
# In this question, you are given complete freedom to choose what and how many features you want to generate. Upon submission to Autolab, we will run linear regression on your generated features and maintain a scoreboard of best regression accuracy (measured by mean squared error). 
#
# The top scoring students will receive a bonus of 5 points. 
#
# ### Tips:
# * Test your features locally by building your model using the training data, and predicting on the validation data. Compute the mean squared error on the **validation dataset** as a metric for how well your features generalize. This helps avoid overfitting to the training dataset, and you'll have faster turnaround time than resubmitting to autolab. 
# * The linear regression model will be trained on your chosen features of the same training examples we provide in this notebook. 
# * We test your regression on a different dataset from the training and validation set that we provide for you, so the MSE you get locally may not match how your features work on the Autolab dataset. 
# * We will solve the linear regression using Ordinary Least Squares with regularization $\lambda=10^{-4}$ and a Cholesky factorization, exactly as done earlier in this notebook. 
# * Note that the argument contains **UNlabeled** data: you cannot build features off the output labels (there is no ETA column). This is in contrast to before, where we kept everything inside the same dataframe for convenience. You can produce the sample input by removing the "eta" column, which we provide code for below. 
# * Make sure your features are all numeric. Try everything!

# +
def contest_features(vdf, vdf_train):
    """ Given a dataframe of UNlabeled and truncated bus data, generate ANY features you'd like for linear regression. 
        Args:
            vdf (dataframe) : dataframe of bus data with truncated rows but unlabeled (no eta column )
                              for which you should produce features
            vdf_train (dataframe) : dataframe of training bus data, truncated and labeled 
        Return: 
            (dataframe) : dataframe of features for each example in vdf
        """
    # create your own engineered features
    pass
    
# AUTOLAB_IGNORE_START
# contest_cols = list(labeled_vdf.columns)
# contest_cols.remove("eta")
# contest_features(labeled_vdf_valid[contest_cols], labeled_vdf).head()
# AUTOLAB_IGNORE_STOP
# -








