BixiLR/BixiLinearAIDA-RWS.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In this Python worflow we explore the Montreal Bixi biking data set for the year 2017 https://www.kaggle.com/aubertsigouin/biximtl/data\n",
    "\n",
    "We have additionally enriched this data set with the biking distance/duration available via Google map API as gmdata2017\n",
    "\n",
    "Our objective is to predict the \"trip duration\", given the distance between two stations."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Import AIDA components"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "from aida.aida import *;"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Connection information to AIDA's server"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "host='cerberus'; dbname='bixi'; user='bixi'; passwd='bixi'; jobName='bixiLinear'; port=55660;"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Establish a connection and get a handle to the database workspace."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "dw = AIDA.connect(host, dbname, user, passwd, jobName, port);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let us see what tables we have in the database"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "      tablename\n",
      "0  stations2017\n",
      "1  tripdata2017\n",
      "2    gmdata2017\n"
     ]
    }
   ],
   "source": [
    "print(dw._tables());"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let us take a peek into tripdata2017.\n",
    "\n",
    "We can see the attributes and explore some sample data.This can be accomplished via the head() or tail() functions (similar to Pandas API) provided by TabularData that sends a sample of data from the server to the client side.\n",
    "\n",
    "Further we can also use the describe() function to look at the data distribution characteristics. This semantic is very similar to the functionality provided by pandas DataFrame to get a summary of the overall distribution of each attribute."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "   id                  starttm  stscode                    endtm  endscode  \\\n",
      "0   0  2017-04-15 00:00:00.000     7060  2017-04-15 00:31:00.000      7060   \n",
      "1   1  2017-04-15 00:01:00.000     6173  2017-04-15 00:10:00.000      6173   \n",
      "2   2  2017-04-15 00:01:00.000     6203  2017-04-15 00:04:00.000      6204   \n",
      "3   3  2017-04-15 00:01:00.000     6104  2017-04-15 00:06:00.000      6114   \n",
      "4   4  2017-04-15 00:01:00.000     6174  2017-04-15 00:11:00.000      6174   \n",
      "\n",
      "   duration  ismember  \n",
      "0      1841         1  \n",
      "1       553         1  \n",
      "2       195         1  \n",
      "3       285         1  \n",
      "4       569         1  \n",
      "                            id                 starttm  \\\n",
      "count   [          4018721.00]  [          4018721.00]   \n",
      "unique  [          4018721.00]  [           227660.00]   \n",
      "nulls   [                0.00]  [                0.00]   \n",
      "max     [          4018721.00]  [ 2017-09-30 23:59:00]   \n",
      "min     [                0.00]  [ 2017-04-15 00:00:00]   \n",
      "avg     [          2009360.21]  [                    ]   \n",
      "median  [          2009360.00]  [                    ]   \n",
      "25%     [          1004680.00]  [                    ]   \n",
      "50%     [          2009360.00]  [                    ]   \n",
      "75%     [          3014040.00]  [                    ]   \n",
      "stddev  [          1160105.11]  [                    ]   \n",
      "\n",
      "                       stscode                   endtm  \\\n",
      "count   [          4018721.00]  [          4018721.00]   \n",
      "unique  [              546.00]  [           227784.00]   \n",
      "nulls   [                0.00]  [                0.00]   \n",
      "max     [            10002.00]  [ 2017-10-01 00:55:00]   \n",
      "min     [             5002.00]  [ 2017-04-15 00:04:00]   \n",
      "avg     [             6324.82]  [                    ]   \n",
      "median  [             6203.00]  [                    ]   \n",
      "25%     [             6105.00]  [                    ]   \n",
      "50%     [             6203.00]  [                    ]   \n",
      "75%     [             6389.00]  [                    ]   \n",
      "stddev  [              375.86]  [                    ]   \n",
      "\n",
      "                      endscode                duration                ismember  \n",
      "count   [          4018721.00]  [          4018721.00]  [          4018721.00]  \n",
      "unique  [              546.00]  [             7115.00]  [                2.00]  \n",
      "nulls   [                0.00]  [                0.00]  [                0.00]  \n",
      "max     [            10002.00]  [             7199.00]  [                1.00]  \n",
      "min     [             5002.00]  [               61.00]  [                0.00]  \n",
      "avg     [             6319.87]  [              837.45]  [                0.80]  \n",
      "median  [             6195.00]  [              670.00]  [                1.00]  \n",
      "25%     [             6092.00]  [              382.00]  [                1.00]  \n",
      "50%     [             6195.00]  [              670.00]  [                1.00]  \n",
      "75%     [             6394.00]  [             1121.00]  [                1.00]  \n",
      "stddev  [              383.28]  [              657.71]  [                0.40]  \n"
     ]
    }
   ],
   "source": [
    "print(dw.tripdata2017.head());\n",
    "print(dw.tripdata2017.describe());"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "So we have 4 million + records in tripdata2017. Also, the station codes are labels. We may have to enrich this information.\n",
    "Let us take a look at the contents of stations2017."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "   scode                      sname  slatitude  slongitude  sispublic\n",
      "0   7060  \"de l'Église / de Verdun\"  45.463001  -73.571569          1\n",
      "1   6173         \"Berri / Cherrier\"  45.519088  -73.569509          1\n",
      "2   6203   \"Hutchison / Sherbrooke\"  45.507810  -73.572080          1\n",
      "3   6204        \"Milton / Durocher\"  45.508144  -73.574772          1\n",
      "4   6104    \"Wolfe / René-Lévesque\"  45.516818  -73.554188          1\n",
      "                         scode                   sname  \\\n",
      "count   [              546.00]  [              546.00]   \n",
      "unique  [              546.00]  [              546.00]   \n",
      "nulls   [                0.00]  [                0.00]   \n",
      "max     [            10002.00]   [Émile-Journault / d]   \n",
      "min     [             5002.00]  [\"10e Avenue / Rosemo]   \n",
      "avg     [             6412.74]  [                    ]   \n",
      "median  [             6305.00]  [                    ]   \n",
      "25%     [             6143.00]  [                    ]   \n",
      "50%     [             6305.00]  [                    ]   \n",
      "75%     [             6723.00]  [                    ]   \n",
      "stddev  [              405.02]  [                    ]   \n",
      "\n",
      "                     slatitude              slongitude               sispublic  \n",
      "count   [              546.00]  [              546.00]  [              546.00]  \n",
      "unique  [              545.00]  [              546.00]  [                2.00]  \n",
      "nulls   [                0.00]  [                0.00]  [                0.00]  \n",
      "max     [               45.58]  [              -73.50]  [                1.00]  \n",
      "min     [               45.43]  [              -73.67]  [                0.00]  \n",
      "avg     [               45.52]  [              -73.58]  [                0.98]  \n",
      "median  [               45.52]  [              -73.58]  [                1.00]  \n",
      "25%     [               45.50]  [              -73.60]  [                1.00]  \n",
      "50%     [               45.52]  [              -73.58]  [                1.00]  \n",
      "75%     [               45.54]  [              -73.56]  [                1.00]  \n",
      "stddev  [                0.03]  [                0.03]  [                0.13]  \n"
     ]
    }
   ],
   "source": [
    "print(dw.stations2017.head());\n",
    "print(dw.stations2017.describe());"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This is good, we have the longitude and latitude associated with each station, which can be used to enrich the tripdata."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Since we have 546 stations, this gives the possibility of 546 x 546 = 298116 possible scenarios for trips. However, we need not be concerned with trips that started and ended at the same station as those are noise. Also, to weed out any further fluctuations in the input data set, we will limit ourselves to only those station combinations which has at the least 50 trips.\n",
    "\n",
    "We can use AIDA's powerful relational API to accomplish this."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "   stscode  endscode  numtrips\n",
      "0     6203      6204       101\n",
      "1     6104      6114       308\n",
      "2     6719      6354        91\n",
      "3     6175      6118        81\n",
      "4     6280      6160        50\n",
      "                       stscode                endscode                numtrips\n",
      "count   [            19300.00]  [            19300.00]  [            19300.00]\n",
      "unique  [              544.00]  [              544.00]  [              684.00]\n",
      "nulls   [                0.00]  [                0.00]  [                0.00]\n",
      "max     [            10002.00]  [            10002.00]  [             2200.00]\n",
      "min     [             5002.00]  [             5002.00]  [               50.00]\n",
      "avg     [             6302.13]  [             6291.78]  [              116.91]\n",
      "median  [             6194.00]  [             6180.00]  [               81.00]\n",
      "25%     [             6100.00]  [             6078.00]  [               62.00]\n",
      "50%     [             6194.00]  [             6180.00]  [               81.00]\n",
      "75%     [             6350.00]  [             6362.00]  [              125.00]\n",
      "stddev  [              349.87]  [              351.21]  [              117.25]\n"
     ]
    }
   ],
   "source": [
    "freqStations = dw.tripdata2017.filter(Q('stscode', 'endscode', CMP.NE)) \\\n",
    "    .aggregate(('stscode','endscode',{COUNT('*'):'numtrips'}), ('stscode','endscode')) \\\n",
    "    .filter(Q('numtrips',C(50), CMP.GTE));\n",
    "print(freqStations.head());\n",
    "print(freqStations.describe());"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can see that there are 19,300 station combinations that is of interest to us.\n",
    "Next stop, we need to include the longitude and latitude information of the start and end stations.\n",
    "\n",
    "This can be done by joining with the station information using AIDA's relational join operator."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "   stscode  endscode  numtrips      stlat     stlong      enlat     enlong\n",
      "0     6203      6204       101  45.507810 -73.572080  45.508144 -73.574772\n",
      "1     6104      6114       308  45.516818 -73.554188  45.523530 -73.551990\n",
      "2     6719      6354        91  45.460729 -73.634073  45.471743 -73.613924\n",
      "3     6175      6118        81  45.520541 -73.567751  45.525048 -73.560036\n",
      "4     6280      6160        50  45.524505 -73.594142  45.532977 -73.581222\n"
     ]
    }
   ],
   "source": [
    "freqStationsCord = freqStations \\\n",
    "    .join(dw.stations2017, ('stscode',), ('scode',), COL.ALL, ({'slatitude':'stlat'}, {'slongitude':'stlong'})) \\\n",
    "    .join(dw.stations2017, ('endscode',), ('scode',), COL.ALL, ({'slatitude':'enlat'}, {'slongitude':'enlong'}));\n",
    "print(freqStationsCord.head());"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "It would be easier if we can translate the coordinates to a distance metric. Python's geopy module supports this computation using Vincenty's formula. This provides us with a distance as crow flies between two coordiantes. This might be a reasonable approximation of actual distance travelled in a trip. \n",
    "\n",
    "Using TabularData's user transform operator, we can generate a dataset which also includes this distance metric. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "   stscode  endscode  numtrips      stlat     stlong      enlat     enlong  \\\n",
      "0     6203      6204       101  45.507810 -73.572080  45.508144 -73.574772   \n",
      "1     6104      6114       308  45.516818 -73.554188  45.523530 -73.551990   \n",
      "2     6719      6354        91  45.460729 -73.634073  45.471743 -73.613924   \n",
      "3     6175      6118        81  45.520541 -73.567751  45.525048 -73.560036   \n",
      "4     6280      6160        50  45.524505 -73.594142  45.532977 -73.581222   \n",
      "\n",
      "   vdistm  \n",
      "0     213  \n",
      "1     765  \n",
      "2    1995  \n",
      "3     783  \n",
      "4    1380  \n"
     ]
    }
   ],
   "source": [
    "def computeDist(tblrData):\n",
    "    import geopy.distance;     #We will use this module to compute distance.\n",
    "    import copy, numpy as np;\n",
    "    #We are going to keep all the columns of the source tabularData object.\n",
    "    data = copy.copy(tblrData.rows); #This only makes a copy of the metadata, but retains original column data\n",
    "    vdistm = data['vdistm'] = np.empty(tblrData.numRows, dtype=int); #add a new empty column to hold distance.\n",
    "    #These are the inputs to Vincenty's formula.\n",
    "    stlat = data['stlat']; stlong = data['stlong']; enlat = data['enlat']; enlong = data['enlong'];\n",
    "    for i in range(0, tblrData.numRows): #populate the distance metric using longitude/latitude of coordinates.\n",
    "        vdistm[i] = int(geopy.distance.distance((stlat[i],stlong[i]), (enlat[i],enlong[i])).meters);\n",
    "    return data;\n",
    "\n",
    "freqStationsDist = freqStationsCord._U(computeDist); #Execute the user transform\n",
    "print(freqStationsDist.head());                      #Take a peek at a sample data."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can next enrich our trip data set with the distance information by joining these computed distances with each trip."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "   id  duration  vdistm\n",
      "0   2       195     213\n",
      "1   3       285     765\n",
      "2   5       620    1995\n",
      "3  12       395     783\n",
      "4  13      1085    1380\n",
      "                            id                duration                  vdistm\n",
      "count   [          2256283.00]  [          2256283.00]  [          2256283.00]\n",
      "unique  [          2256283.00]  [             6663.00]  [             3652.00]\n",
      "nulls   [                0.00]  [                0.00]  [                0.00]\n",
      "max     [          4018721.00]  [             7199.00]  [             9074.00]\n",
      "min     [                2.00]  [               61.00]  [               47.00]\n",
      "avg     [          2003455.03]  [              630.45]  [             1369.89]\n",
      "median  [          2005024.00]  [              482.00]  [             1117.00]\n",
      "25%     [           994510.00]  [              298.00]  [              711.00]\n",
      "50%     [          2005024.00]  [              482.00]  [             1117.00]\n",
      "75%     [          3010571.00]  [              793.00]  [             1755.00]\n",
      "stddev  [          1162662.25]  [              533.44]  [              921.29]\n"
     ]
    }
   ],
   "source": [
    "tripData = dw.tripdata2017.join(freqStationsDist, ('stscode','endscode'), ('stscode', 'endscode')\n",
    "                                             , ('id', 'duration'), ('vdistm',));\n",
    "print(tripData.head());\n",
    "print(tripData.describe());"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "So we have trip duration for each trip and the distance as crow flies, between the two stations involved in the trip.\n",
    "\n",
    "Also, we have about 2 million trips for which we have distance between stations metric.\n",
    "Given that there are only a few thousand unique values for distance, we might want to keep some values of distance apart for testing.\n",
    "For this purpose, we will first get distinct values for distance and then sort it."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "   vdistm\n",
      "0      47\n",
      "1      71\n",
      "2      81\n",
      "3      85\n",
      "4      88\n",
      "   vdistm\n",
      "0    8529\n",
      "1    8752\n",
      "2    8860\n",
      "3    9031\n",
      "4    9074\n",
      "                        vdistm\n",
      "count   [             3652.00]\n",
      "unique  [             3652.00]\n",
      "nulls   [                0.00]\n",
      "max     [             9074.00]\n",
      "min     [               47.00]\n",
      "avg     [             2203.77]\n",
      "median  [             2013.00]\n",
      "25%     [             1096.00]\n",
      "50%     [             2013.00]\n",
      "75%     [             3068.00]\n",
      "stddev  [             1386.05]\n"
     ]
    }
   ],
   "source": [
    "uniqueTripDist = tripData[:,['vdistm']].distinct().order('vdistm');\n",
    "print(uniqueTripDist.head());\n",
    "print(uniqueTripDist.tail());\n",
    "print(uniqueTripDist.describe());"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We will keep some data apart for testing. A rule of thumb is 30%. The neat trick below sets apart 33%, across the entire range of distance values. close enough."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "   vdistm\n",
      "0      47\n",
      "1      85\n",
      "2     110\n",
      "3     126\n",
      "4     148\n",
      "   vdistm\n",
      "0    6796\n",
      "1    7057\n",
      "2    7530\n",
      "3    8752\n",
      "4    9074\n",
      "                        vdistm\n",
      "count   [             1218.00]\n",
      "unique  [             1218.00]\n",
      "nulls   [                0.00]\n",
      "max     [             9074.00]\n",
      "min     [               47.00]\n",
      "avg     [             2205.06]\n",
      "median  [             2012.00]\n",
      "25%     [             1095.00]\n",
      "50%     [             2012.00]\n",
      "75%     [             3069.00]\n",
      "stddev  [             1390.92]\n"
     ]
    }
   ],
   "source": [
    "testTripDist = uniqueTripDist[::3];\n",
    "print(testTripDist.head());\n",
    "print(testTripDist.tail());\n",
    "print(testTripDist.describe());"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now let us get the remaining values for distances to be used for training."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "   vdistm\n",
      "0      71\n",
      "1      81\n",
      "2      88\n",
      "3      94\n",
      "4     120\n",
      "   vdistm\n",
      "0    7307\n",
      "1    7650\n",
      "2    8529\n",
      "3    8860\n",
      "4    9031\n",
      "                        vdistm\n",
      "count   [             2434.00]\n",
      "unique  [             2434.00]\n",
      "nulls   [                0.00]\n",
      "max     [             9031.00]\n",
      "min     [               71.00]\n",
      "avg     [             2203.13]\n",
      "median  [             2013.00]\n",
      "25%     [             1096.00]\n",
      "50%     [             2013.00]\n",
      "75%     [             3068.00]\n",
      "stddev  [             1383.60]\n"
     ]
    }
   ],
   "source": [
    "trainTripDist = uniqueTripDist.filter(Q('vdistm', testTripDist, CMP.NOTIN));\n",
    "print(trainTripDist.head());\n",
    "print(trainTripDist.tail());\n",
    "print(trainTripDist.describe());"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let us now extract the fields of interest to us for the training data, which is just the distance of each trip and it's duration"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "   vdistm  duration\n",
      "0     213       195\n",
      "1    1995       620\n",
      "2     495       565\n",
      "3     499       280\n",
      "4     802       252\n",
      "   vdistm  duration\n",
      "0    1010       364\n",
      "1    1318       288\n",
      "2    1300       692\n",
      "3     668       179\n",
      "4     956       306\n",
      "                        vdistm                duration\n",
      "count   [          1503408.00]  [          1503408.00]\n",
      "unique  [             2434.00]  [             6267.00]\n",
      "nulls   [                0.00]  [                0.00]\n",
      "max     [             9031.00]  [             7199.00]\n",
      "min     [               71.00]  [               61.00]\n",
      "avg     [             1370.76]  [              629.42]\n",
      "median  [             1109.00]  [              480.00]\n",
      "25%     [              706.00]  [              295.00]\n",
      "50%     [             1109.00]  [              480.00]\n",
      "75%     [             1774.00]  [              793.00]\n",
      "stddev  [              926.28]  [              536.08]\n"
     ]
    }
   ],
   "source": [
    "trainData = tripData.project(('vdistm', 'duration')).filter(Q('vdistm', trainTripDist, CMP.IN));\n",
    "print(trainData.head());\n",
    "print(trainData.tail());\n",
    "print(trainData.describe());"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As the values are huge, we should normalize the data attributes. First get the max values for these attributes."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "9074\n",
      "7199\n"
     ]
    }
   ],
   "source": [
    "maxdist = uniqueTripDist.max('vdistm');\n",
    "print(maxdist);\n",
    "maxduration = tripData.max('duration');\n",
    "print(maxduration);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now let us normalize the training data. As we are working with integer data, we will also have to convert it to float. That can be accomplished by multiplying with 1.0."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "     vdistm  duration\n",
      "0  0.023474  0.027087\n",
      "1  0.219859  0.086123\n",
      "2  0.054551  0.078483\n",
      "3  0.054992  0.038894\n",
      "4  0.088384  0.035005\n",
      "     vdistm  duration\n",
      "0  0.111307  0.050563\n",
      "1  0.145250  0.040006\n",
      "2  0.143266  0.096124\n",
      "3  0.073617  0.024865\n",
      "4  0.105356  0.042506\n"
     ]
    }
   ],
   "source": [
    "trainData = trainData.project((1.0*F('vdistm')/maxdist, 1.0*F('duration')/maxduration));\n",
    "\n",
    "print(trainData.head());\n",
    "print(trainData.tail());"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Our linear regression equation is of the form.\n",
    "\n",
    "dur = a + b*dist\n",
    "\n",
    "we will re-organize the training data set to fit this format and also setup our initial parameters for a and b. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "    x0    vdistm\n",
      "0  1.0  0.023474\n",
      "1  1.0  0.219859\n",
      "2  1.0  0.054551\n",
      "3  1.0  0.054992\n",
      "4  1.0  0.088384\n",
      "   duration\n",
      "0  0.027087\n",
      "1  0.086123\n",
      "2  0.078483\n",
      "3  0.038894\n",
      "4  0.035005\n",
      "OrderedDict([('a', array([1.])), ('b', array([1.]))])\n"
     ]
    }
   ],
   "source": [
    "trainDataSet = dw._ones((trainData.numRows, 1), (\"x0\",)).hstack(trainData[:,['vdistm']]);\n",
    "print(trainDataSet.head());\n",
    "trainDataSetDuration = trainData[:,['duration']];\n",
    "print(trainDataSetDuration.head());\n",
    "params = dw._ones((1,2), (\"a\",\"b\"));\n",
    "print(params.rows);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let us try to run a prediction using these parameters."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "VirtualOrderedColumnsDict([(0, <aidas.dborm.DBTable.Column object at 0x7f3741179940>)])\n",
      "   r_0000000000\n",
      "0      1.023474\n",
      "1      1.219859\n",
      "2      1.054551\n",
      "3      1.054992\n",
      "4      1.088384\n"
     ]
    }
   ],
   "source": [
    "pred = trainDataSet @ params.T;\n",
    "print(pred.columns);\n",
    "print(pred.head());"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We need to compute the squared error for the predictions. Since we will be reusing them, we might as well store it as a function."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def squaredErr(actual, predicted):\n",
    "    return ((predicted-actual)**2).sum()/(2*(actual.shape[0]));"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let us see what is the error for the first iteration."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.5694865536695626\n"
     ]
    }
   ],
   "source": [
    "sqerr = squaredErr(trainDataSetDuration, pred);\n",
    "print(sqerr);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We need to perform a gradient descent based on the squared errors. We will write another function to perform this."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def gradDesc(actual, predicted, indata):\n",
    "    return (predicted-actual).T @ indata / actual.shape[0];"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let us update our params using gradient descent using the error we got. We also need to use a learning rate, alpha (arbitrarily chosen)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "OrderedDict([('a', array([0.89363664])), ('b', array([0.98330567]))])\n"
     ]
    }
   ],
   "source": [
    "alpha = 0.1;\n",
    "\n",
    "params = params - alpha * gradDesc(trainDataSetDuration, pred, trainDataSet);\n",
    "print(params.rows);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now let us try to use the updated params to train the model again and see if the error is decreasing."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "   r_0000000000\n",
      "0      0.916718\n",
      "1      1.109825\n",
      "2      0.947277\n",
      "3      0.947711\n",
      "4      0.980546\n",
      "0.4594973598186901\n"
     ]
    }
   ],
   "source": [
    "pred = trainDataSet @ params.T;\n",
    "print(pred.head());\n",
    "sqerr = squaredErr(trainDataSetDuration, pred);\n",
    "print(sqerr);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Before we proceed, may be we should check if google maps API's distance metric gives a better learning rate. Let us see what fields we can use from Google."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "   stscode  endscode  gdistm  gduration\n",
      "0     6406      6052    3568        596\n",
      "1     6050      6406    3821        704\n",
      "2     6148      6173    1078        293\n",
      "3     6110      6114    1319        337\n",
      "4     6123      6114     725        177\n",
      "                       stscode                endscode  \\\n",
      "count   [            19516.00]  [            19516.00]   \n",
      "unique  [              544.00]  [              544.00]   \n",
      "nulls   [                0.00]  [                0.00]   \n",
      "max     [            10002.00]  [            10002.00]   \n",
      "min     [             5002.00]  [             5002.00]   \n",
      "avg     [             6302.62]  [             6291.74]   \n",
      "median  [             6194.00]  [             6180.00]   \n",
      "25%     [             6100.00]  [             6078.00]   \n",
      "50%     [             6194.00]  [             6180.00]   \n",
      "75%     [             6350.00]  [             6362.00]   \n",
      "stddev  [              350.30]  [              350.46]   \n",
      "\n",
      "                        gdistm               gduration  \n",
      "count   [            19516.00]  [            19516.00]  \n",
      "unique  [             4903.00]  [             1540.00]  \n",
      "nulls   [                0.00]  [                0.00]  \n",
      "max     [            14530.00]  [             3083.00]  \n",
      "min     [               18.00]  [                4.00]  \n",
      "avg     [             2079.29]  [              516.19]  \n",
      "median  [             1766.00]  [              459.00]  \n",
      "25%     [             1118.00]  [              300.00]  \n",
      "50%     [             1766.00]  [              459.00]  \n",
      "75%     [             2711.00]  [              671.00]  \n",
      "stddev  [             1345.59]  [              299.18]  \n"
     ]
    }
   ],
   "source": [
    "print(dw.gmdata2017.head());\n",
    "print(dw.gmdata2017.describe());"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can build a new data set for the trips between frequently used station combination that includes google's distance."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "   id  duration  gdistm  gduration\n",
      "0   2       195     288        218\n",
      "1   3       285    1007        296\n",
      "2   5       620    2587        538\n",
      "3  12       395    1615        322\n",
      "4  13      1085    1710        352\n",
      "                            id                duration  \\\n",
      "count   [          2256283.00]  [          2256283.00]   \n",
      "unique  [          2256283.00]  [             6663.00]   \n",
      "nulls   [                0.00]  [                0.00]   \n",
      "max     [          4018721.00]  [             7199.00]   \n",
      "min     [                2.00]  [               61.00]   \n",
      "avg     [          2003455.03]  [              630.45]   \n",
      "median  [          2005024.00]  [              482.00]   \n",
      "25%     [           994510.00]  [              298.00]   \n",
      "50%     [          2005024.00]  [              482.00]   \n",
      "75%     [          3010571.00]  [              793.00]   \n",
      "stddev  [          1162662.25]  [              533.44]   \n",
      "\n",
      "                        gdistm               gduration  \n",
      "count   [          2256283.00]  [          2256283.00]  \n",
      "unique  [             4880.00]  [             1532.00]  \n",
      "nulls   [                0.00]  [                0.00]  \n",
      "max     [            14530.00]  [             3083.00]  \n",
      "min     [               48.00]  [               11.00]  \n",
      "avg     [             1834.73]  [              454.68]  \n",
      "median  [             1497.00]  [              394.00]  \n",
      "25%     [              960.00]  [              257.00]  \n",
      "50%     [             1497.00]  [              394.00]  \n",
      "75%     [             2359.00]  [              591.00]  \n",
      "stddev  [             1236.26]  [              273.36]  \n"
     ]
    }
   ],
   "source": [
    "gtripData = dw.gmdata2017 \\\n",
    "    .join(dw.tripdata2017, ('stscode','endscode'), ('stscode', 'endscode'), COL.ALL, COL.ALL) \\\n",
    "    .join(freqStations, ('stscode','endscode'), ('stscode', 'endscode') \\\n",
    "                          , ('id', 'duration', 'gdistm', 'gduration') );\n",
    "print(gtripData.head());\n",
    "print(gtripData.describe());"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Google also provides its estimated duration for the trip. We will have to see in the end if our trained model is able to predict the trip duration better than google's estimate. So we will also save Google's estimate for the trip duration for that comparison."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Next up, we need to format this dataset the same way we did the first one."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "14530\n",
      "7199\n"
     ]
    }
   ],
   "source": [
    "guniqueTripDist = gtripData[:,['gdistm']].distinct().order('gdistm');\n",
    "gtestTripDist = guniqueTripDist[::3];\n",
    "gtrainTripDist = guniqueTripDist.filter(Q('gdistm', gtestTripDist, CMP.NOTIN));\n",
    "gtrainData = gtripData.project(('gdistm', 'duration')).filter(Q('gdistm', gtrainTripDist, CMP.IN));\n",
    "\n",
    "gmaxdist = guniqueTripDist.max('gdistm');\n",
    "print(gmaxdist);\n",
    "gmaxduration = gtripData.max('duration');\n",
    "print(gmaxduration);\n",
    "gtrainData = gtrainData.project((1.0*F('gdistm')/gmaxdist, 1.0*F('duration')/gmaxduration));\n",
    "\n",
    "gtrainDataSet = dw._ones((gtrainData.numRows, 1), (\"x0\",)).hstack(gtrainData[:,['gdistm']]);\n",
    "gtrainDataSetDuration = gtrainData[:,['duration']];\n",
    "gparams = dw._ones((1,2), (\"a\",\"b\"));"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let us see how the error rate is progressing for the new dataset."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.5419675497480458\n",
      "0.4379084758530774\n"
     ]
    }
   ],
   "source": [
    "gpred = gtrainDataSet @ gparams.T;\n",
    "gsqerr = squaredErr(gtrainDataSetDuration, gpred);\n",
    "print(gsqerr);\n",
    "gparams = gparams - alpha * gradDesc(gtrainDataSetDuration, gpred, gtrainDataSet);\n",
    "gpred = gtrainDataSet @ gparams.T;\n",
    "gsqerr = squaredErr(gtrainDataSetDuration, gpred);\n",
    "print(gsqerr);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "It looks like using Google maps' distance is giving us a slight advantage. That makes sense, since Vincenty's formula computes distances as a crow flies, where as Google maps' distance metric is based on the actual road network distances. Better data gives better prediction results !\n",
    "\n",
    "We are done with the feature selection and feature engineering phase for now.\n",
    "\n",
    "Next we will proceed to train our linear regression model using the training data set.\n",
    "\n",
    "Meanwhile, we will also let it printout the error rate at frequent intervals so that we know it is decreasing."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Error rate after 100 iterations is 0.002414885788154281\n",
      "Error rate after 200 iterations is 0.0023518940367528206\n",
      "Error rate after 300 iterations is 0.0022972381907587756\n",
      "Error rate after 400 iterations is 0.0022498149319043433\n",
      "Error rate after 500 iterations is 0.0022086671747070705\n",
      "Error rate after 600 iterations is 0.0021729644844393505\n",
      "Error rate after 700 iterations is 0.002141986317482605\n",
      "Error rate after 800 iterations is 0.0021151074794789697\n",
      "Error rate after 900 iterations is 0.002091785507800193\n",
      "Error rate after 1000 iterations is 0.00207154972368958\n",
      "OrderedDict([('a', array([0.00295944])), ('b', array([0.67162878]))])\n",
      "0.00207154972368958\n"
     ]
    }
   ],
   "source": [
    "def trainModel(w, numiters, alpha):\n",
    "    gtrainDataSet = w.gtrainDataSet; gtrainDataSetDuration=w.gtrainDataSetDuration;\n",
    "    gparams = w.gparams; gradDesc = w.gradDesc;\n",
    "    for i in range(0, numiters):\n",
    "        gpred = gtrainDataSet @ gparams.T;\n",
    "        gparams = gparams - alpha*gradDesc(gtrainDataSetDuration, gpred, gtrainDataSet);\n",
    "    w.gpred = gpred; w.gparams = gparams;\n",
    "\n",
    "#Export any objects and functions required for execution in the remote workspace.\n",
    "dw.gtrainDataSet = gtrainDataSet; dw.gtrainDataSetDuration = gtrainDataSetDuration; \n",
    "dw.gparams = gparams; dw.gradDesc = gradDesc;\n",
    "\n",
    "for i in range(0, 10):\n",
    "    dw._X(trainModel, 100, alpha);\n",
    "    print(\"Error rate after {} iterations is {}\".format((i+1)*100, squaredErr(gtrainDataSetDuration, dw.gpred)))\n",
    "\n",
    "gparams = dw.gparams; gpred = dw.gpred;\n",
    "print(gparams.rows);\n",
    "gsqerr = squaredErr(gtrainDataSetDuration, gpred);\n",
    "print(gsqerr);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let us see how our model performs in predictions against the test data set we had kept apart."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "99215.98424571125\n"
     ]
    }
   ],
   "source": [
    "gtestData = gtripData.project(('gdistm', 'duration', 'gduration')).filter(Q('gdistm', gtestTripDist, CMP.IN));\n",
    "gtestData = gtestData.project((1.0*F('gdistm')/gmaxdist, 1.0*F('duration')/gmaxduration, 'gduration'));\n",
    "gtestDataSet = dw._ones((gtestData.numRows, 1), (\"x0\",)).hstack(gtestData[:,['gdistm']]);\n",
    "gtestDataSetDuration = gtestData[:,['duration']];\n",
    "\n",
    "gtestpred = gtestDataSet @ gparams.T;\n",
    "\n",
    "gtestsqerr1 = squaredErr(gtestDataSetDuration*gmaxduration, gtestpred*gmaxduration);\n",
    "print(gtestsqerr1);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We would also like to check how the duration provided by Google maps' API hold up to the test data set."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "111763.37983591038\n"
     ]
    }
   ],
   "source": [
    "gtestsqerr2 = squaredErr(gtestDataSetDuration*gmaxduration, gtestData[:,['gduration']]);\n",
    "print(gtestsqerr2);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "So yes, our model is able to do a better job."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "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.5.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}