{"metadata":{"kernelspec":{"language":"python","display_name":"Python 3","name":"python3"},"language_info":{"pygments_lexer":"ipython3","nbconvert_exporter":"python","version":"3.6.4","file_extension":".py","codemirror_mode":{"name":"ipython","version":3},"name":"python","mimetype":"text/x-python"}},"nbformat_minor":4,"nbformat":4,"cells":[{"cell_type":"code","source":"import numpy as np\nimport pandas as pd\nfrom pathlib import Path\nfrom tqdm import tqdm\nimport matplotlib.pyplot as plt\nfrom sklearn.model_selection import train_test_split\n","metadata":{"execution":{"iopub.status.busy":"2021-05-29T16:19:33.882574Z","iopub.execute_input":"2021-05-29T16:19:33.882954Z","iopub.status.idle":"2021-05-29T16:19:34.834732Z","shell.execute_reply.started":"2021-05-29T16:19:33.882870Z","shell.execute_reply":"2021-05-29T16:19:34.833527Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"**Loading train, test & submission files**","metadata":{}},{"cell_type":"code","source":"trainfile = pd.read_csv(\"../input/google-smartphone-decimeter-challenge/baseline_locations_train.csv\")\ntestfile = pd.read_csv(\"../input/google-smartphone-decimeter-challenge/baseline_locations_test.csv\")\nsubmission = pd.read_csv(\"../input/google-smartphone-decimeter-challenge/sample_submission.csv\")","metadata":{"execution":{"iopub.status.busy":"2021-05-29T16:19:34.836627Z","iopub.execute_input":"2021-05-29T16:19:34.837076Z","iopub.status.idle":"2021-05-29T16:19:35.480672Z","shell.execute_reply.started":"2021-05-29T16:19:34.837030Z","shell.execute_reply":"2021-05-29T16:19:35.479462Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"trainfile","metadata":{"execution":{"iopub.status.busy":"2021-05-29T16:19:58.729704Z","iopub.execute_input":"2021-05-29T16:19:58.730121Z","iopub.status.idle":"2021-05-29T16:19:58.769773Z","shell.execute_reply.started":"2021-05-29T16:19:58.730071Z","shell.execute_reply":"2021-05-29T16:19:58.768406Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"**Extracting ground truths and aligning with train inputs**","metadata":{}},{"cell_type":"code","source":"datapath = Path(\"../input/google-smartphone-decimeter-challenge\")\ntruths = (datapath / 'train').rglob('ground_truth.csv')\n\n\ncols = ['collectionName', 'phoneName', 'millisSinceGpsEpoch','phone','heightAboveWgs84EllipsoidM','latDeg',\n       'lngDeg',]\ntruth_list =[]\nfor filepath in tqdm(truths, total=73):\n    file = pd.read_csv(filepath, usecols=cols)\n    truth_list.append(file)\n    \ntruth_data = pd.concat(truth_list, ignore_index=True)\n\ntrain = trainfile[cols]\n\ntrain = train.merge(truth_data.iloc[:,3:], suffixes=(\"_current\",\"_truth\",\"_truth\"))","metadata":{"execution":{"iopub.status.busy":"2021-05-29T16:38:51.558101Z","iopub.execute_input":"2021-05-29T16:38:51.558482Z","iopub.status.idle":"2021-05-29T16:38:51.609436Z","shell.execute_reply.started":"2021-05-29T16:38:51.558436Z","shell.execute_reply":"2021-05-29T16:38:51.606433Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"train","metadata":{"execution":{"iopub.status.busy":"2021-05-29T16:37:18.824815Z","iopub.execute_input":"2021-05-29T16:37:18.825184Z","iopub.status.idle":"2021-05-29T16:37:18.838242Z","shell.execute_reply.started":"2021-05-29T16:37:18.825144Z","shell.execute_reply":"2021-05-29T16:37:18.836555Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"train.info()","metadata":{"execution":{"iopub.status.busy":"2021-05-29T16:37:58.264089Z","iopub.execute_input":"2021-05-29T16:37:58.264666Z","iopub.status.idle":"2021-05-29T16:37:58.282376Z","shell.execute_reply.started":"2021-05-29T16:37:58.264609Z","shell.execute_reply":"2021-05-29T16:37:58.281026Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"**Checking the correlation between current & truth coordinates**","metadata":{}},{"cell_type":"code","source":"train.iloc[:,3:].corr()","metadata":{"execution":{"iopub.status.busy":"2021-05-28T08:53:52.77168Z","iopub.execute_input":"2021-05-28T08:53:52.772247Z","iopub.status.idle":"2021-05-28T08:53:52.79727Z","shell.execute_reply.started":"2021-05-28T08:53:52.772208Z","shell.execute_reply":"2021-05-28T08:53:52.796303Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"* In above correlation table, we can see that ***latDeg_current is highly correlated with latDeg_truth*** and ***lngDeg_current is highly correlated with lngDeg_truth***\n* That would help us to build a baseline model","metadata":{}},{"cell_type":"code","source":"test = testfile.copy()\n\nprint(\"############### collectionName unique values ##############################\")\nprint(\"train: {}\".format(train.collectionName.nunique()))\nprint(train.collectionName.unique())\nprint(\"----------------------------------------------\")\nprint(\"test: {}\".format(test.collectionName.nunique()))\nprint(test.collectionName.unique())\nprint(\"----------------------------------------------\")\n\nprint(\"\\n\")\n\nprint(\"############### phoneName unique values ##############################\")\nprint(\"train: {}\".format(train.phoneName.nunique()))\nprint(train.phoneName.unique())\nprint(\"----------------------------------------------\")\nprint(\"test: {}\".format(test.phoneName.nunique()))\nprint(test.phoneName.unique())\nprint(\"----------------------------------------------\")","metadata":{"execution":{"iopub.status.busy":"2021-05-28T08:53:53.732483Z","iopub.execute_input":"2021-05-28T08:53:53.732843Z","iopub.status.idle":"2021-05-28T08:53:53.810316Z","shell.execute_reply.started":"2021-05-28T08:53:53.732811Z","shell.execute_reply":"2021-05-28T08:53:53.809214Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"**Note:** Here, collectionName would be useless to consider as it would be different in train and test data but **phoneName** might be helpful","metadata":{}},{"cell_type":"markdown","source":"**One-hot encoding for 'phoneName' of train and test data**","metadata":{}},{"cell_type":"code","source":"train_phone = pd.get_dummies(train.loc[:,\"phoneName\"])\ntest_phone = pd.get_dummies(test.loc[:,\"phoneName\"])\n\nprint(\"train_phone shape:{}\".format(train_phone.shape))\nprint(\"test_phone shape:{}\".format(test_phone.shape))","metadata":{"execution":{"iopub.status.busy":"2021-05-28T08:53:53.813493Z","iopub.execute_input":"2021-05-28T08:53:53.813813Z","iopub.status.idle":"2021-05-28T08:53:53.838966Z","shell.execute_reply.started":"2021-05-28T08:53:53.813782Z","shell.execute_reply":"2021-05-28T08:53:53.83795Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"**Aligning train & test columns**","metadata":{}},{"cell_type":"code","source":"train_phone1, test_phone1 = train_phone.align(test_phone, join=\"outer\",axis=1, fill_value=0)\nprint(\"Updated train shape {}\".format(train_phone.shape))\n\nprint(\"Updated test shape {}\".format(test_phone.shape))","metadata":{"execution":{"iopub.status.busy":"2021-05-28T08:53:53.84105Z","iopub.execute_input":"2021-05-28T08:53:53.841418Z","iopub.status.idle":"2021-05-28T08:53:53.849046Z","shell.execute_reply.started":"2021-05-28T08:53:53.84138Z","shell.execute_reply":"2021-05-28T08:53:53.847837Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"* Basically, there is no change is one-hot encoded phoneName columns for train or test data","metadata":{}},{"cell_type":"markdown","source":"**Train-test data after addition of one-hot encoded columns**","metadata":{}},{"cell_type":"code","source":"train1 = pd.concat([train.iloc[:,3:], train_phone1], axis=1, ignore_index=False)\ntest1 = pd.concat([test.iloc[:,3:5], test_phone1], axis=1, ignore_index=False)","metadata":{"execution":{"iopub.status.busy":"2021-05-28T08:53:53.850869Z","iopub.execute_input":"2021-05-28T08:53:53.851268Z","iopub.status.idle":"2021-05-28T08:53:53.868253Z","shell.execute_reply.started":"2021-05-28T08:53:53.851228Z","shell.execute_reply":"2021-05-28T08:53:53.867451Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"print(\"train1_shape:\",train1.shape)\ntrain1.columns","metadata":{"execution":{"iopub.status.busy":"2021-05-28T08:53:53.869494Z","iopub.execute_input":"2021-05-28T08:53:53.869869Z","iopub.status.idle":"2021-05-28T08:53:53.877961Z","shell.execute_reply.started":"2021-05-28T08:53:53.869815Z","shell.execute_reply":"2021-05-28T08:53:53.876775Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"print(\"test1-shape:\",test1.shape)\ntest1.columns","metadata":{"execution":{"iopub.status.busy":"2021-05-28T08:53:53.879709Z","iopub.execute_input":"2021-05-28T08:53:53.880332Z","iopub.status.idle":"2021-05-28T08:53:53.889112Z","shell.execute_reply.started":"2021-05-28T08:53:53.880292Z","shell.execute_reply":"2021-05-28T08:53:53.887862Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"**Plotting current vs truth coordinates**","metadata":{}},{"cell_type":"code","source":"plt.figure(figsize=[10,5])\nplt.plot(train1[\"latDeg_current\"][:200],train1[\"lngDeg_current\"][:200],\"bo\",label=\"current\")\nplt.plot(train1[\"latDeg_truth\"][:200],train1[\"lngDeg_truth\"][:200],\"r*\",label=\"truth\")\nplt.title(\"current vs truth\", fontweight=\"bold\")\nplt.xlabel(\"latDeg\")\nplt.ylabel(\"lngDeg\")\nplt.legend()\n\nplt.figure(figsize=[15,5])\nplt.subplot(1,2,1)\nplt.plot(train1[\"latDeg_current\"][:200],train1[\"latDeg_truth\"][:200],\"bo\")\nplt.title(\"lat (current vs truth)\", fontweight=\"bold\")\n\nplt.subplot(1,2,2)\nplt.plot(train1[\"lngDeg_current\"][:200],train1[\"lngDeg_truth\"][:200],\"bo\")\nplt.title(\"lng (current vs truth)\", fontweight=\"bold\")","metadata":{"execution":{"iopub.status.busy":"2021-05-28T08:53:53.890763Z","iopub.execute_input":"2021-05-28T08:53:53.891142Z","iopub.status.idle":"2021-05-28T08:53:54.345724Z","shell.execute_reply.started":"2021-05-28T08:53:53.891106Z","shell.execute_reply":"2021-05-28T08:53:54.344807Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"**Input-output**","metadata":{}},{"cell_type":"code","source":"# for lat\ncolumns1 = list(train1.columns[4:])\nX1,y1 = train1.loc[:,[train1.columns[0]]+columns1],train1[\"latDeg_truth\"].values\n  \n#for lng\ncolumns2 = list(train1.columns[4:])\nX2,y2 = train1.loc[:,[train1.columns[1]]+columns2],train1[\"lngDeg_truth\"].values\n\n\nXt1 = test1.loc[:,[test1.columns[0]]+columns1] #for lat\nXt2 = test1.loc[:,[test1.columns[1]]+columns2] #for lng\n\nprint(\"X1 columns:\", X1.columns.tolist())\nprint(\"Xt1 columns:\", Xt1.columns.tolist())\nprint(\"\\n\")\nprint(\"X2 columns:\", X2.columns.tolist())\nprint(\"Xt2 columns:\", Xt2.columns.tolist())","metadata":{"execution":{"iopub.status.busy":"2021-05-28T08:53:54.347132Z","iopub.execute_input":"2021-05-28T08:53:54.347482Z","iopub.status.idle":"2021-05-28T08:53:54.366946Z","shell.execute_reply.started":"2021-05-28T08:53:54.347444Z","shell.execute_reply":"2021-05-28T08:53:54.365953Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"xtr1,xval1,ytr1,yval1 = train_test_split(X1, y1, test_size=0.3, random_state=10)\nxtr2,xval2,ytr2,yval2 = train_test_split(X2, y2, test_size=0.3, random_state=10)\n\nprint(\"xtr1 shape:{}; xval1 shape:{}\".format(xtr1.shape,xval1.shape))\nprint(\"xtr2 shape:{}; xval2 shape:{}\".format(xtr2.shape,xval2.shape))","metadata":{"execution":{"iopub.status.busy":"2021-05-28T08:53:54.368271Z","iopub.execute_input":"2021-05-28T08:53:54.368633Z","iopub.status.idle":"2021-05-28T08:53:54.402706Z","shell.execute_reply.started":"2021-05-28T08:53:54.368596Z","shell.execute_reply":"2021-05-28T08:53:54.401732Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"**Defining a function to get prediction and ground truth distance estimation (in meters)**\n\n[Haversine formula for distance estimation using GPS co-ordinates](https://stackoverflow.com/questions/15736995/how-can-i-quickly-estimate-the-distance-between-two-latitude-longitude-points)","metadata":{}},{"cell_type":"code","source":"from math import radians, cos, sin, asin, sqrt\ndef lat_lon_dist(df):\n    \"\"\"\n    Calculate the great circle distance between two points \n    on the earth (specified in decimal degrees)\n    \"\"\"\n    dist_list = []\n    for i in tqdm(range(df.shape[0]),total=100):\n        lat1 = df[\"latDeg_truth\"][i]\n        lon1 = df[\"lngDeg_truth\"][i]\n        lat2 = df[\"latDeg_pred\"][i]\n        lon2 = df[\"lngDeg_pred\"][i]\n        # convert decimal degrees to radians \n        lon1, lat1, lon2, lat2 = map(radians, [lon1, lat1, lon2, lat2])\n        # haversine formula \n        dlon = lon2 - lon1 \n        dlat = lat2 - lat1 \n        a = sin(dlat/2)**2 + cos(lat1) * cos(lat2) * sin(dlon/2)**2\n        c = 2 * asin(sqrt(a)) \n        # Radius of earth in kilometers is 6371\n        mdist = 6371* c*1000\n        dist_list.append(mdist)\n    \n    return dist_list","metadata":{"execution":{"iopub.status.busy":"2021-05-28T08:53:54.404092Z","iopub.execute_input":"2021-05-28T08:53:54.404458Z","iopub.status.idle":"2021-05-28T08:53:54.41284Z","shell.execute_reply.started":"2021-05-28T08:53:54.404403Z","shell.execute_reply":"2021-05-28T08:53:54.411713Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"idx_Mi8=np.where(xtr1[\"Mi8\"]==1)[0]\nidx_Pixel4=np.where(xtr1[\"Pixel4\"]==1)[0]\nidx_Pixel4Modded=np.where(xtr1[\"Pixel4Modded\"]==1)[0]\nidx_Pixel4XL=np.where(xtr1[\"Pixel4XL\"]==1)[0]\nidx_Pixel4XLModded=np.where(xtr1[\"Pixel4XLModded\"]==1)[0]\nidx_Pixel5=np.where(xtr1[\"Pixel5\"]==1)[0]\nidx_SamsungS20Ultra=np.where(xtr1[\"SamsungS20Ultra\"]==1)[0]","metadata":{"execution":{"iopub.status.busy":"2021-05-28T08:53:54.414408Z","iopub.execute_input":"2021-05-28T08:53:54.414773Z","iopub.status.idle":"2021-05-28T08:53:54.429563Z","shell.execute_reply.started":"2021-05-28T08:53:54.414733Z","shell.execute_reply":"2021-05-28T08:53:54.428609Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"xtr1","metadata":{"execution":{"iopub.status.busy":"2021-05-28T08:57:17.187568Z","iopub.execute_input":"2021-05-28T08:57:17.18807Z","iopub.status.idle":"2021-05-28T08:57:17.251892Z","shell.execute_reply.started":"2021-05-28T08:57:17.187963Z","shell.execute_reply":"2021-05-28T08:57:17.250401Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"**Model building**","metadata":{}},{"cell_type":"code","source":"#lr1 = LinearRegression() #selected as a starting point\n#model_lat = lr1.fit(xtr1.to_numpy()[idx_Mi8,0].reshape(-1,1),ytr1)\n#pred_yval1 = model_lat.predict(xval1) # prediction for val data (lat)\n#lr2 = LinearRegression()\n#model_lng = lr2.fit(xtr2,ytr2)\n#pred_yval2 = model_lng.predict(xval2) # prediction for val data (long)\n\nfrom sklearn.gaussian_process import GaussianProcessRegressor\nfrom sklearn.gaussian_process.kernels import RBF, ConstantKernel as C\n\n# Instantiate a Gaussian Process model\ngp1_Mi8 = GaussianProcessRegressor()\ngp1_Mi8.fit(xtr1.to_numpy()[idx_Mi8,0].reshape(-1,1),np.array(ytr1[idx_Mi8].ravel()).reshape(-1,1))\n\n\n\n","metadata":{"execution":{"iopub.status.busy":"2021-05-28T08:53:54.431649Z","iopub.execute_input":"2021-05-28T08:53:54.432179Z","iopub.status.idle":"2021-05-28T08:54:16.692681Z","shell.execute_reply.started":"2021-05-28T08:53:54.432143Z","shell.execute_reply":"2021-05-28T08:54:16.691761Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"**Model evaluation on validation data**","metadata":{}},{"cell_type":"code","source":"val_df = pd.concat([xval1[[\"latDeg_current\"]], xval2], ignore_index=False, axis=1).reset_index(drop=[\"index\"])","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"**Adding truth & predicted lat-long values to val_df**","metadata":{}},{"cell_type":"code","source":"#truth\nval_df[\"latDeg_truth\"] = yval1\nval_df[\"lngDeg_truth\"] = yval2\n\n#pred\nval_df[\"latDeg_pred\"] = pred_yval1\nval_df[\"lngDeg_pred\"] = pred_yval2","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"val_df[\"dist\"] = lat_lon_dist(val_df)","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"phone = val_df.iloc[:,2:-5].idxmax(axis=1) #Reversing one-hot decoding for phoneName\n\nval_df1 = pd.concat([val_df.iloc[:,:2],val_df.iloc[:,-3:]], axis=1, ignore_index=False)\nval_df1[\"phoneName\"] = phone\n\nval_df1 = val_df1[val_df1.columns[-1:].tolist()+val_df1.columns[:-1].tolist()]","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"**Box-plot analysis for dist  analysis for each phone**","metadata":{}},{"cell_type":"code","source":"import seaborn as sns\nplt.figure(figsize=[15,7])\n\n# ax, fig = plt.subplots(figsize=[15,7])\nsns.boxplot(x=\"phoneName\", y=\"dist\",data=val_df1)\nplt.ylabel(\"Dist (m)\") # distance in meters\n#plt.ylim([0,30]) # for better visualization","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"**Preparing evaluation score for each phone (50th & 95th percentile)**","metadata":{}},{"cell_type":"code","source":"val_df2 = pd.DataFrame()\nval_df2[\"phoneName\"] =  val_df1.phoneName.unique().tolist()\nval_df2[\"dist_50\"] = [np.percentile(val_df1[val_df1.phoneName==ph][\"dist\"],50) for ph in val_df2[\"phoneName\"].tolist()]\nval_df2[\"dist_95\"] = [np.percentile(val_df1[val_df1.phoneName==ph][\"dist\"],95) for ph in val_df2[\"phoneName\"].tolist()]\nval_df2[\"avg_dist_50_95\"] = np.mean(np.array(val_df2.iloc[:,1:]),axis=1)\nprint(\"Val evaluation details:\\n\",val_df2)\n\nprint(\"\\n\")\nprint(\"------------------------------------------------------\")\nprint(\"Final val evaluation score: {}\".format(val_df2.iloc[:,-1].mean()))\nprint(\"------------------------------------------------------\")","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"**Training the model on complete data and predict for test data**","metadata":{}},{"cell_type":"code","source":"lr1 = LinearRegression()\nmodel_lat = lr1.fit(X1,y1)\npred_yt1 = model_lat.predict(Xt1) # prediction for test data (lat)\n\nlr2 = LinearRegression()\nmodel_lng = lr2.fit(X2,y2)\npred_yt2 = model_lng.predict(Xt2) # prediction for test data (long)","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"submission = test[['phone','millisSinceGpsEpoch']]\npd.options.mode.chained_assignment = None  # default='warn'\nsubmission['latDeg'] = pred_yt1.tolist()\nsubmission['lngDeg'] = pred_yt2.tolist()","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"submission.to_csv(\"./submission.csv\",index=False)","metadata":{"trusted":true},"execution_count":null,"outputs":[]}]}