{"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":"markdown","source":"This notebook is the most accurate notebook that [colum2131](https://www.kaggle.com/columbia2131), [tubo213](https://www.kaggle.com/tubotubo) and [penguin46](https://www.kaggle.com/ryotayoshinobu) created before merging with [chris](https://www.kaggle.com/chris62) and [\nAkio Saito](https://www.kaggle.com/saitodevel01). Descriptions of each process can be found in [this discussion](https://www.kaggle.com/c/google-smartphone-decimeter-challenge/discussion/261739). Please comment if there are any unclear points.","metadata":{}},{"cell_type":"code","source":"%%capture\n!pip install pandarallel \n!pip install vincenty\n!pip install simdkalman","metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","execution":{"iopub.status.busy":"2021-08-05T09:08:18.453097Z","iopub.execute_input":"2021-08-05T09:08:18.453493Z","iopub.status.idle":"2021-08-05T09:08:49.581298Z","shell.execute_reply.started":"2021-08-05T09:08:18.453462Z","shell.execute_reply":"2021-08-05T09:08:49.579912Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"import warnings\nwarnings.simplefilter('ignore')\n\nfrom contextlib import contextmanager\nfrom glob import glob\nfrom time import time\nimport pickle\n\nimport lightgbm as lgbm\nimport matplotlib.pyplot as plt\nimport numpy as np\nimport pandas as pd\nimport plotly.express as px\nimport seaborn as sns\nfrom pathlib import Path\nimport torch\nfrom scipy import interpolate\nfrom sklearn.metrics import accuracy_score, confusion_matrix\nfrom sklearn.model_selection import GroupKFold, KFold, StratifiedKFold\nfrom sklearn.neighbors import KNeighborsRegressor\nfrom tqdm.notebook import tqdm\nfrom tqdm._tqdm_notebook import tqdm_notebook as tqdm\ntqdm.pandas()\n\nimport simdkalman\nfrom pandarallel import pandarallel\npandarallel.initialize()\nfrom vincenty import vincenty\nimport cupy as cp","metadata":{"execution":{"iopub.status.busy":"2021-08-05T09:22:53.569237Z","iopub.execute_input":"2021-08-05T09:22:53.569588Z","iopub.status.idle":"2021-08-05T09:22:53.581609Z","shell.execute_reply.started":"2021-08-05T09:22:53.569557Z","shell.execute_reply":"2021-08-05T09:22:53.580161Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"import pyproj\nfrom pyproj import Proj, transform\n\ndef calc_haversine(lat1, lon1, lat2, lon2):\n    \"\"\"Calculates the great circle distance between two points\n    on the earth. Inputs are array-like and specified in decimal degrees.\n    \"\"\"\n    lat1, lon1, lat2, lon2 = map(np.radians, [lat1, lon1, lat2, lon2])\n    dlat = lat2 - lat1\n    dlon = lon2 - lon1\n    a = np.sin(dlat/2.0)**2 + np.cos(lat1) * np.cos(lat2) * np.sin(dlon/2.0)**2\n\n    c = 2 * np.arcsin(a**0.5)\n    dist = 6_367_000 * c\n    return dist\n\ndef WGS84_to_ECEF(lat, lon, alt):\n    # convert to radians\n    rad_lat = lat * (np.pi / 180.0)\n    rad_lon = lon * (np.pi / 180.0)\n    a = 6378137.0\n    # f is the flattening factor\n    finv = 298.257223563\n    f = 1 / finv   \n    # e is the eccentricity\n    e2 = 1 - (1 - f) * (1 - f)    \n    # N is the radius of curvature in the prime vertical\n    N = a / np.sqrt(1 - e2 * np.sin(rad_lat) * np.sin(rad_lat))\n    x = (N + alt) * np.cos(rad_lat) * np.cos(rad_lon)\n    y = (N + alt) * np.cos(rad_lat) * np.sin(rad_lon)\n    z = (N * (1 - e2) + alt)        * np.sin(rad_lat)\n    return x, y, z\n\n\ntransformer = pyproj.Transformer.from_crs(\n    {\"proj\":'geocent', \"ellps\":'WGS84', \"datum\":'WGS84'},\n    {\"proj\":'latlong', \"ellps\":'WGS84', \"datum\":'WGS84'},)\n\n\ndef ECEF_to_WGS84(x,y,z):\n    lng, lat, alt = transformer.transform(x,y,z,radians=False)\n    return lat, lng, alt\n\n\nW2E = lambda r: WGS84_to_ECEF(r.latDeg, r.lngDeg, r.heightKNN)\nE2W = lambda r: ECEF_to_WGS84(r.x, r.y, r.z)\n\n\ndef visualize_trafic(df:pd.DataFrame,lat='latDeg',lng='lngDeg',color=None,savepath=None,zoom=9,center={\"lat\":37.423576, \"lon\":-122.094132}):\n    fig = px.scatter_mapbox(df,\n                            \n                            # Here, plotly gets, (x,y) coordinates\n                            lat=lat,\n                            lon=lng,                          \n                            #Here, plotly detects color of series\n                            color=color,\n                            zoom=zoom,\n                            center=center,\n                            height=600,\n                            width=800)\n    fig.update_layout(mapbox_style='stamen-terrain')\n    fig.update_layout(margin={\"r\": 0, \"t\": 0, \"l\": 0, \"b\": 0})\n    fig.update_layout(title_text=\"GPS trafic\")\n    fig.show()\n    if savepath != None:\n        fig.write_html(savepath)","metadata":{"execution":{"iopub.status.busy":"2021-08-05T09:08:55.376518Z","iopub.execute_input":"2021-08-05T09:08:55.376873Z","iopub.status.idle":"2021-08-05T09:08:55.530637Z","shell.execute_reply.started":"2021-08-05T09:08:55.376843Z","shell.execute_reply":"2021-08-05T09:08:55.529363Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Reading-data","metadata":{}},{"cell_type":"code","source":"def vincenty_meter(r, lat='latDeg', lng='lngDeg', tlat='t_latDeg', tlng='t_lngDeg'):\n    return vincenty((r[lat], r[lng]), (r[tlat], r[tlng])) * 1000\n\n\ndef check_meter(input_df: pd.DataFrame, save=False):\n    output_df = input_df.copy()\n    \n    output_df['meter'] = input_df.parallel_apply(vincenty_meter, axis=1)\n    if save == True:\n        output_df.to_csv('train_output.csv', index=False)\n\n    meter_score = output_df['meter'].mean()\n    print(f'meter: {meter_score}') # 2.533116208067488\n\n    scores = []\n    for phone in output_df['phone'].unique():\n        p_50 = np.percentile(output_df.loc[output_df['phone']==phone, 'meter'], 50)\n        p_95 = np.percentile(output_df.loc[output_df['phone']==phone, 'meter'], 95)\n        scores.append(p_50)\n        scores.append(p_95)\n\n    score = sum(scores) / len(scores)\n    print(f'CV: {score}') # 3.53009109589041\n    \n    return output_df\n\n\ndef get_groundtruth():\n    output_df = pd.DataFrame()\n    \n    for path in glob(str(BASE_DIR / 'train/*/*/ground_truth.csv')):\n        _df = pd.read_csv(path)\n        output_df = pd.concat([output_df, _df])\n    output_df = output_df.reset_index(drop=True)\n    \n    _columns = ['latDeg', 'lngDeg', 'heightAboveWgs84EllipsoidM']\n    output_df[['t_'+col for col in _columns]] = output_df[_columns]\n    output_df = output_df.drop(columns=_columns, axis=1)\n    return output_df","metadata":{"execution":{"iopub.status.busy":"2021-08-05T09:08:55.532907Z","iopub.execute_input":"2021-08-05T09:08:55.533554Z","iopub.status.idle":"2021-08-05T09:08:55.547867Z","shell.execute_reply.started":"2021-08-05T09:08:55.5335Z","shell.execute_reply":"2021-08-05T09:08:55.546299Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"DEVICE = torch.device(\"cuda:0\" if torch.cuda.is_available() else \"cpu\")\nBASE_DIR = Path('../input/google-smartphone-decimeter-challenge')\n\ntrain_base = pd.read_csv(BASE_DIR / 'baseline_locations_train.csv')\ntest_base = pd.read_csv(BASE_DIR / 'baseline_locations_test.csv')\nsample = pd.read_csv(BASE_DIR / 'sample_submission.csv')\n\ntrain_base = train_base.merge(\n    get_groundtruth(), on=['collectionName', 'phoneName', 'millisSinceGpsEpoch']\n)\ntrain_base = check_meter(train_base)","metadata":{"execution":{"iopub.status.busy":"2021-08-05T09:08:55.550023Z","iopub.execute_input":"2021-08-05T09:08:55.550578Z","iopub.status.idle":"2021-08-05T09:09:07.050017Z","shell.execute_reply.started":"2021-08-05T09:08:55.55052Z","shell.execute_reply":"2021-08-05T09:09:07.048071Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# kNN(height)","metadata":{}},{"cell_type":"code","source":"def fit_knn_height(X: pd.DataFrame, y: pd.Series, n_neighbors=15):\n    model = KNeighborsRegressor(n_neighbors=n_neighbors, weights='distance')\n    model.fit(X.values, y.values)\n    return model\n\n\ndef pred_knn_height(X: pd.DataFrame, model):\n    return model.predict(X.values)","metadata":{"execution":{"iopub.status.busy":"2021-08-05T09:09:07.052054Z","iopub.execute_input":"2021-08-05T09:09:07.052556Z","iopub.status.idle":"2021-08-05T09:09:07.060712Z","shell.execute_reply.started":"2021-08-05T09:09:07.052506Z","shell.execute_reply":"2021-08-05T09:09:07.059003Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"knn_model = fit_knn_height(\n    train_base[['t_latDeg', 't_lngDeg']],\n    train_base['t_heightAboveWgs84EllipsoidM'])\n\ntrain_base['heightKNN'] = pred_knn_height(train_base[['latDeg', 'lngDeg']], knn_model)\ntest_base['heightKNN'] = pred_knn_height(test_base[['latDeg', 'lngDeg']], knn_model)\ntrain_base['heightDiff'] = np.abs(train_base['heightAboveWgs84EllipsoidM'] - train_base['heightKNN'])\ntest_base['heightDiff']  = np.abs(test_base['heightAboveWgs84EllipsoidM'] - test_base['heightKNN'])\n\ndel knn_model","metadata":{"execution":{"iopub.status.busy":"2021-08-05T09:09:07.062506Z","iopub.execute_input":"2021-08-05T09:09:07.063454Z","iopub.status.idle":"2021-08-05T09:09:08.098576Z","shell.execute_reply.started":"2021-08-05T09:09:07.063409Z","shell.execute_reply":"2021-08-05T09:09:08.097493Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# outlier detection","metadata":{}},{"cell_type":"code","source":"def calc_haversine(lat1, lon1, lat2, lon2):\n    lat1, lon1, lat2, lon2 = map(np.radians, [lat1, lon1, lat2, lon2])\n    dlat = lat2 - lat1\n    dlon = lon2 - lon1\n    a = np.sin(dlat/2.0)**2 + np.cos(lat1) * np.cos(lat2) * np.sin(dlon/2.0)**2\n\n    c = 2 * np.arcsin(a**0.5)\n    dist = 6_367_000 * c\n    return dist\n\ndef generate_waypoints(train=train_base.copy()):\n    train['area'] = train['phone'].apply(lambda s : s.split('-')[4])\n    train[\"phone\"] = train.collectionName.str.cat(train.phoneName, sep=\"_\")\n    train = train[train.area==\"SJC\"]\n\n    # augment by interpolate\n    dfs = []\n    for c, c_df in tqdm(train.groupby(\"collectionName\")):    \n        c_df = c_df.sort_values(\"millisSinceGpsEpoch\")  \n\n        f_lat = interpolate.interp1d(c_df.millisSinceGpsEpoch, c_df.t_latDeg, kind='linear')\n        f_lng = interpolate.interp1d(c_df.millisSinceGpsEpoch, c_df.t_lngDeg, kind='linear')\n\n        start_time = c_df.millisSinceGpsEpoch.min()\n        end_time = c_df.millisSinceGpsEpoch.max()\n\n        times = range(start_time, end_time, 19)\n        lats = f_lat(times)\n        lngs = f_lng(times)\n\n        dfs.append(pd.DataFrame({\n            \"collectionName\": c,\n            \"latDeg\": lats,\n            \"lngDeg\": lngs,          \n            \"gt_time\": times,\n        })) \n\n    return pd.concat(dfs).drop_duplicates().reset_index(drop=True).copy()\n\ndef closest_point(point, waypoints):\n    point = cp.radians(cp.array([point]))\n    waypoints  = cp.array(waypoints)\n    _waypoints = cp.radians(cp.array(waypoints))\n    diffs = point - _waypoints\n    a = cp.sin(diffs[:,0]/2.0)**2 + cp.cos(_waypoints[:,0])*cp.cos(point[:,0])*cp.sin(diffs[:,1]/2.0)**2\n    c = 2 * cp.arcsin(a**0.5)\n    dist = 6_367_000 * c\n    return waypoints[dist.argmin()].get().tolist()\n\ndef calc_dist(pred, waypoints):\n    pred = cp.radians(cp.array(pred))\n    waypoints = cp.radians(cp.array(waypoints))\n    diffs = pred - waypoints\n    a = cp.sin(diffs[:,0]/2.0)**2 + cp.cos(waypoints[:,0])*cp.cos(pred[:,0])*cp.sin(diffs[:,1]/2.0)**2\n    c = 2 * cp.arcsin(a**0.5)\n    dist = 6_367_000 * c\n    return dist.get()\n\ndef add_latlng(df):\n    df['latlng'] = [[x, y] for x,y in zip(df['latDeg'], df['lngDeg'])]\n    return df\n\ndef snap_to_grid(sub, train_waypoints, threshhold=1e10):\n    sub = add_latlng(sub)\n    sub['matched_point'] = sub['latlng'].apply(lambda x:closest_point(x,train_waypoints))\n    sub['dist'] = calc_dist(sub['latlng'].tolist(),sub['matched_point'].tolist())\n\n    sub['_latDeg_'] = sub['latlng'].apply(lambda x:x[0])\n    sub['_lngDeg_'] = sub['latlng'].apply(lambda x:x[1])\n\n    sub.loc[sub['dist'] < threshhold,'_latDeg_'] = sub['matched_point'].apply(lambda x:x[0])\n    sub.loc[sub['dist'] < threshhold,'_lngDeg_'] = sub['matched_point'].apply(lambda x:x[1])\n    return sub[[\"phone\", \"millisSinceGpsEpoch\", \"dist\"]].copy()\n\ndef processing(input_df: pd.DataFrame, L=25):\n    output_df = pd.DataFrame(dtype=np.float64)\n    shift_list = list(range(-L, L+1, 1))\n\n    for i in shift_list:\n        output_df = pd.concat([\n            output_df,\n            input_df.groupby('phone')[['heightDiff', 'dist']]\\\n                .shift(i).add_prefix(f'shift{i}_')\n        ], axis=1)\n\n    for i in shift_list:\n        if i == 0: continue\n        output_df = pd.concat([\n            output_df,\n            input_df.groupby('phone')[['latDeg', 'lngDeg']]\\\n                .diff(i).add_prefix(f'diff{i}_')\n        ], axis=1)\n    output_df = pd.concat([\n        output_df,\n        input_df.groupby('phone')[['latDeg', 'lngDeg']]\\\n            .pct_change(i).add_prefix(f'change{i}_')\n    ], axis=1)\n\n    return output_df\n\n@contextmanager\ndef timer(logger=None, format_str='{:.3f}[s]', prefix=None, suffix=None):\n    if prefix: format_str = str(prefix) + format_str\n    if suffix: format_str = format_str + str(suffix)\n    start = time()\n    yield\n    d = time() - start\n    out_str = format_str.format(d)\n    if logger:\n        logger.info(out_str)\n    else:\n        print(out_str)\n\ndef fit_lgbm(X, y, train_df=train_base, params: dict=None, verbose=100, seed: int=42):\n    models = []\n    oof_pred = np.zeros(len(y), dtype=np.float64)\n    \n    kf = GroupKFold(n_splits=N_SPLITS)\n    for i, (idx_train, idx_valid) in enumerate(kf.split(X, y, train_df['collectionName'].reset_index(drop=True))):\n        x_train, y_train = X.iloc[idx_train], y.iloc[idx_train]\n        x_valid, y_valid = X.iloc[idx_valid], y.iloc[idx_valid]\n\n        model = lgbm.LGBMClassifier(**params)\n        with timer(prefix='fit fold={} '.format(i + 1)):\n            model.fit(x_train, y_train, \n                eval_set=[(x_valid, y_valid)],  \n                early_stopping_rounds=verbose, \n                eval_metric='logloss',\n                verbose=verbose)\n            \n        pred_i = model.predict_proba(x_valid)[:, 1]\n        oof_pred[x_valid.index] = pred_i\n        models.append(model)\n\n    return oof_pred, models\n\ndef predict_lgbm(models, feat_df):\n    pred = np.array([model.predict_proba(feat_df.values)[:, 1] for model in models])\n    pred = np.mean(pred, axis=0)\n    return pred","metadata":{"execution":{"iopub.status.busy":"2021-08-05T09:09:08.102514Z","iopub.execute_input":"2021-08-05T09:09:08.102995Z","iopub.status.idle":"2021-08-05T09:09:08.153347Z","shell.execute_reply.started":"2021-08-05T09:09:08.10295Z","shell.execute_reply":"2021-08-05T09:09:08.152021Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### Only SJC","metadata":{}},{"cell_type":"code","source":"# Generate roads\ntrain_waypoints = generate_waypoints()\ntrain_base['area'] = train_base['phone'].apply(lambda s : s.split('-')[4])\ntest_base['area'] = test_base['phone'].apply(lambda s : s.split('-')[4])\n\ntrain_base = train_base.merge(\n    snap_to_grid(\n        train_base.loc[train_base.area==\"SJC\"],\n        train_waypoints[[\"latDeg\", \"lngDeg\"]]\n    ), on=[\"phone\", \"millisSinceGpsEpoch\"], how=\"outer\")\n\ntest_base = test_base.merge(\n    snap_to_grid(\n        test_base.loc[test_base.area==\"SJC\"],\n        train_waypoints[[\"latDeg\", \"lngDeg\"]]\n    ), on=[\"phone\", \"millisSinceGpsEpoch\"], how=\"outer\")","metadata":{"execution":{"iopub.status.busy":"2021-08-05T09:10:54.423297Z","iopub.execute_input":"2021-08-05T09:10:54.423757Z","iopub.status.idle":"2021-08-05T09:12:56.137042Z","shell.execute_reply.started":"2021-08-05T09:10:54.423725Z","shell.execute_reply":"2021-08-05T09:12:56.135842Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"loss_threshold = 9.5\nvisibility     = 22\n\ntrain_base['t_isOutlier'] = (train_base.meter > loss_threshold).astype(int)\n\nparams = {\n 'reg_alpha': 0.01,\n 'reg_lambda': 0.01, \n 'num_leaves': 40,\n 'n_estimators': 10000,\n 'learning_rate': 0.1,\n 'random_state': 42,\n 'max_depth': -1\n}\n\nN_SPLITS = 3\noof, models = fit_lgbm(processing(train_base[train_base.area==\"SJC\"].reset_index(drop=True), visibility), \n                       train_base.loc[train_base.area==\"SJC\", 't_isOutlier'], \n                       train_df=train_base[train_base.area==\"SJC\"],\n                       params=params)\npred = predict_lgbm(models, processing(test_base[test_base.area==\"SJC\"].reset_index(drop=True), visibility))\n\ntrain_base.loc[train_base.area==\"SJC\", 'isOutlier'] = (oof > 0.5).astype(int)\ntest_base.loc[test_base.area==\"SJC\", 'isOutlier'] = (pred > 0.5).astype(int)\ntest_base.loc[test_base.collectionName==\"2021-04-02-US-SJC-1\", 'isOutlier'] = 0\n\n\nprint('score', accuracy_score((oof>0.5).astype(int), train_base.loc[train_base.area==\"SJC\", 't_isOutlier']))\nprint(confusion_matrix((oof>0.5).astype(int), train_base.loc[train_base.area==\"SJC\", 't_isOutlier']))","metadata":{"execution":{"iopub.status.busy":"2021-08-05T09:12:56.139111Z","iopub.execute_input":"2021-08-05T09:12:56.139612Z","iopub.status.idle":"2021-08-05T09:13:37.882084Z","shell.execute_reply.started":"2021-08-05T09:12:56.13957Z","shell.execute_reply":"2021-08-05T09:13:37.880837Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### OtherSJC","metadata":{}},{"cell_type":"code","source":"def processing(input_df: pd.DataFrame, L=25):\n    output_df = pd.DataFrame(dtype=np.float64)\n    shift_list = list(range(-L, L+1, 1))\n\n    for i in shift_list:\n        if i == 0: continue\n        output_df = pd.concat([\n            output_df,\n            input_df.groupby('phone')[['latDeg', 'lngDeg']]\\\n                .diff(i).add_prefix(f'diff{i}_')\n        ], axis=1)\n    output_df = pd.concat([\n        output_df,\n        input_df.groupby('phone')[['latDeg', 'lngDeg']]\\\n            .pct_change(i).add_prefix(f'change{i}_')\n    ], axis=1)\n\n    return output_df","metadata":{"execution":{"iopub.status.busy":"2021-08-05T09:13:37.885071Z","iopub.execute_input":"2021-08-05T09:13:37.885975Z","iopub.status.idle":"2021-08-05T09:13:37.894688Z","shell.execute_reply.started":"2021-08-05T09:13:37.88593Z","shell.execute_reply":"2021-08-05T09:13:37.893489Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"loss_threshold = 5.6 \nvisibility     = 26   \n\ntrain_base['t_isOutlier'] = (train_base.meter > loss_threshold).astype(int)\n\nparams = {\n 'reg_alpha': 0.01,\n 'reg_lambda': 0.01, \n 'num_leaves': 40,\n 'n_estimators': 10000,\n 'learning_rate': 0.1,\n 'random_state': 42,\n 'max_depth': -1\n}\n\nN_SPLITS = 5\nuse_train_index = train_base.area!=\"SJC\"\nuse_test_index = (test_base.area!=\"SJC\") | (test_base.collectionName==\"2021-04-02-US-SJC-1\")\noof, models = fit_lgbm(processing(train_base[use_train_index].reset_index(drop=True), visibility), \n                       train_base.loc[use_train_index, 't_isOutlier'].reset_index(drop=True), \n                       train_df=train_base[use_train_index].reset_index(drop=True),\n                       params=params)\npred = predict_lgbm(models, processing(test_base[use_test_index].reset_index(drop=True), visibility))\n\ntrain_base.loc[use_train_index, 'isOutlier'] = (oof > 0.5).astype(int)\ntest_base.loc[use_test_index, 'isOutlier'] = (pred > 0.5).astype(int)\n\n\nprint('score', accuracy_score((oof>0.5).astype(int), train_base.loc[use_train_index, 't_isOutlier']))\nprint(confusion_matrix((oof>0.5).astype(int), train_base.loc[use_train_index, 't_isOutlier']))","metadata":{"execution":{"iopub.status.busy":"2021-08-05T09:13:37.897402Z","iopub.execute_input":"2021-08-05T09:13:37.897886Z","iopub.status.idle":"2021-08-05T09:16:08.812322Z","shell.execute_reply.started":"2021-08-05T09:13:37.897841Z","shell.execute_reply":"2021-08-05T09:16:08.810907Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### Interpolate outliers","metadata":{}},{"cell_type":"code","source":"train_base.loc[train_base.isOutlier==1, [\"latDeg\", \"lngDeg\", \"heightAboveWgs84EllipsoidM\"]] = np.nan\ntest_base.loc[test_base.isOutlier==1, [\"latDeg\", \"lngDeg\", \"heightAboveWgs84EllipsoidM\"]] = np.nan\n\ndfs = []\nfor phone, phone_df in train_base.groupby(\"phone\"):\n    phone_df.loc[:, [\"latDeg\", \"lngDeg\", 'heightAboveWgs84EllipsoidM']] \\\n        = phone_df.loc[:, [\"latDeg\", \"lngDeg\", \"heightAboveWgs84EllipsoidM\"]].interpolate(limit_area=None, limit_direction='both')\n    dfs.append(phone_df.copy())\ntrain_base = pd.concat(dfs).reset_index(drop=True)\n\ndfs = []\nfor phone, phone_df in test_base.groupby(\"phone\"):\n    phone_df.loc[:, [\"latDeg\", \"lngDeg\", \"heightAboveWgs84EllipsoidM\"]] \\\n         = phone_df.loc[:, [\"latDeg\", \"lngDeg\", \"heightAboveWgs84EllipsoidM\"]].interpolate(limit_area=None, limit_direction='both')\n    dfs.append(phone_df.copy())\ntest_base = pd.concat(dfs).reset_index(drop=True)\n\ntrain_base[\"isOutlier\"] = 0\ntest_base[\"isOutlier\"] = 0","metadata":{"execution":{"iopub.status.busy":"2021-08-05T09:16:08.8144Z","iopub.execute_input":"2021-08-05T09:16:08.814864Z","iopub.status.idle":"2021-08-05T09:16:09.565694Z","shell.execute_reply.started":"2021-08-05T09:16:08.81482Z","shell.execute_reply":"2021-08-05T09:16:09.564583Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"train_base = check_meter(train_base)","metadata":{"execution":{"iopub.status.busy":"2021-08-05T09:16:09.567413Z","iopub.execute_input":"2021-08-05T09:16:09.568172Z","iopub.status.idle":"2021-08-05T09:16:18.568781Z","shell.execute_reply.started":"2021-08-05T09:16:09.568126Z","shell.execute_reply":"2021-08-05T09:16:18.566792Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# StopMean","metadata":{}},{"cell_type":"code","source":"def processing(input_df: pd.DataFrame):\n    output_df = pd.DataFrame(dtype=np.float64)\n    shift_list = list(range(-15, 16, 1))\n    \n    for i in shift_list:\n        if i == 0: continue\n        output_df = pd.concat([\n            output_df,\n            input_df.groupby('phone')[['latDeg', 'lngDeg']]\\\n                .diff(i).add_prefix(f'diff{i}_')\n        ], axis=1)\n        \n        output_df = pd.concat([\n            output_df,\n            input_df.groupby('phone')[['latDeg', 'lngDeg']]\\\n                .pct_change(i).add_prefix(f'change{i}_')\n        ], axis=1)\n    \n    return output_df\n\n@contextmanager\ndef timer(logger=None, format_str='{:.3f}[s]', prefix=None, suffix=None):\n    if prefix: format_str = str(prefix) + format_str\n    if suffix: format_str = format_str + str(suffix)\n    start = time()\n    yield\n    d = time() - start\n    out_str = format_str.format(d)\n    if logger:\n        logger.info(out_str)\n    else:\n        print(out_str)\n        \n\ndef fit_lgbm(X, y, params: dict=None, verbose=100, seed: int=42):\n    models = []\n    oof_pred = np.zeros(len(y), dtype=np.float64)\n    \n    kf = GroupKFold(n_splits=N_SPLITS)\n    for i, (idx_train, idx_valid) in enumerate(kf.split(X, y, train_base['collectionName'])):\n        x_train, y_train = X.iloc[idx_train], y.iloc[idx_train]\n        x_valid, y_valid = X.iloc[idx_valid], y.iloc[idx_valid]\n\n        model = lgbm.LGBMClassifier(**params)\n        with timer(prefix='fit fold={} '.format(i + 1)):\n            model.fit(x_train, y_train, \n                eval_set=[(x_valid, y_valid)],  \n                early_stopping_rounds=verbose, \n                eval_metric='logloss',\n                verbose=verbose)\n            \n        pred_i = model.predict_proba(x_valid)[:, 1]\n        oof_pred[x_valid.index] = pred_i\n        models.append(model)\n\n    return oof_pred, models\n\n\ndef predict_lgbm(models, feat_df):\n    pred = np.array([model.predict_proba(feat_df.values)[:, 1] for model in models])\n    pred = np.mean(pred, axis=0)\n    return pred\n\n\ndef stopmean(input_df: pd.DataFrame):\n    output_df = input_df.copy()\n    stop_index = []\n    \n    for _, sub_df in input_df.groupby('phone'): \n        _index = sub_df.index\n        _stops = sub_df['stop'].tolist()\n        for i in range(1, len(_stops)-1):\n            if _stops[i-1] == 1 and  _stops[i] == 0 and  _stops[i+1] == 1:\n                _stops[i] = 1\n\n        before = 1 if _stops[0] == 1 else 0\n\n        tmp = []\n        if before:\n            tmp.append(_index[0])\n\n        for i, flag in enumerate(_stops):\n            if flag == 1 and before == 0:\n                tmp.append(_index[i]) \n            elif flag == 0 and before == 1:\n                tmp.append(_index[i])\n                stop_index.append(tmp)\n                tmp = []\n            before = flag\n            \n        if tmp:\n            tmp.append(_index[-1]+1)\n            stop_index.append(tmp)\n\n    output_df['stop_id'] = 0\n    for i, (indexi, indexj) in enumerate(stop_index):\n        output_df.iloc[indexi:indexj]['stop_id'] = i+1\n\n    stopid2deg = output_df.loc[output_df['stop_id']!=0].groupby('stop_id')[['latDeg', 'lngDeg']].mean()\n    stopid2deg.columns = ['stop_latDeg', 'stop_lngDeg']\n\n    output_df = output_df.merge(stopid2deg, on='stop_id', how='left')\n    output_df.loc[output_df['stop_id']!=0, ['latDeg', 'lngDeg']] = output_df.loc[output_df['stop_id']!=0, ['stop_latDeg', 'stop_lngDeg']].values\n    \n    output_df = output_df.drop(columns=['stop_latDeg', 'stop_lngDeg'], axis=1)\n    \n    return output_df","metadata":{"execution":{"iopub.status.busy":"2021-08-05T09:16:18.574574Z","iopub.execute_input":"2021-08-05T09:16:18.575092Z","iopub.status.idle":"2021-08-05T09:16:18.620729Z","shell.execute_reply.started":"2021-08-05T09:16:18.575041Z","shell.execute_reply":"2021-08-05T09:16:18.61899Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"train_base['target_stop'] = (train_base['speedMps']==0).astype(int)\n\nparams = {\n 'reg_alpha': 0.01,\n 'reg_lambda': 0.01, \n 'num_leaves': 40,\n 'n_estimators': 10000,\n 'learning_rate': 0.1,\n 'random_state': 42,\n 'max_depth': -1\n}\n\nN_SPLITS = 5\noof, models = fit_lgbm(processing(train_base), train_base['target_stop'], params=params)\npred = predict_lgbm(models, processing(test_base))\n\ntrain_base['stop'] = (oof > 0.5).astype(int)\ntest_base['stop'] = (pred > 0.5).astype(int)\n\nprint('score', accuracy_score((oof>0.5).astype(int), train_base['target_stop']))\nprint(confusion_matrix((oof>0.5).astype(int), train_base['target_stop']))\n\ndel train_base['target_stop'] \ndel models","metadata":{"execution":{"iopub.status.busy":"2021-08-05T09:16:18.626675Z","iopub.execute_input":"2021-08-05T09:16:18.62699Z","iopub.status.idle":"2021-08-05T09:18:09.561525Z","shell.execute_reply.started":"2021-08-05T09:16:18.62696Z","shell.execute_reply":"2021-08-05T09:18:09.5604Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"train_base = stopmean(train_base)\ntest_base = stopmean(test_base)","metadata":{"execution":{"iopub.status.busy":"2021-08-05T09:18:09.564242Z","iopub.execute_input":"2021-08-05T09:18:09.564791Z","iopub.status.idle":"2021-08-05T09:18:10.258736Z","shell.execute_reply.started":"2021-08-05T09:18:09.564747Z","shell.execute_reply":"2021-08-05T09:18:10.257647Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"train_base = check_meter(train_base)","metadata":{"execution":{"iopub.status.busy":"2021-08-05T09:18:10.260909Z","iopub.execute_input":"2021-08-05T09:18:10.261566Z","iopub.status.idle":"2021-08-05T09:18:18.509232Z","shell.execute_reply.started":"2021-08-05T09:18:10.261511Z","shell.execute_reply":"2021-08-05T09:18:18.507075Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Kalman Filter","metadata":{}},{"cell_type":"code","source":"def interpolate_beforekm(input_df: pd.DataFrame):\n    first_dict = dict(input_df.groupby('phone')['millisSinceGpsEpoch'].first())\n    last_dict = dict(input_df.groupby('phone')['millisSinceGpsEpoch'].last())\n    columns = ['phone', 'millisSinceGpsEpoch', 'latDeg', 'lngDeg', 'isOutlier']\n\n    time_df = pd.DataFrame()\n    for phone in input_df['phone'].unique():\n        _list = np.arange(first_dict[phone], last_dict[phone]+1000, 1000, dtype=int)\n        times = input_df.loc[input_df['phone']==phone, 'millisSinceGpsEpoch'].unique()\n        \n        _df = pd.DataFrame({\n            'phone': phone,\n            'millisSinceGpsEpoch': _list})\n        _df = pd.concat([\n            _df.loc[~_df['millisSinceGpsEpoch'].isin(times)],\n            input_df.loc[input_df['phone']==phone, columns]])\n\n        _df['target'] = 0\n        _df.loc[_df['millisSinceGpsEpoch'].isin(_list), 'target'] = 1\n\n        # interpolate\n        _df = _df.sort_values('millisSinceGpsEpoch')\n        _df.index =_df['millisSinceGpsEpoch'].values\n        _df['latDeg'] = _df['latDeg'].interpolate(method='index', limit_direction='both')\n        _df['lngDeg'] = _df['lngDeg'].interpolate(method='index', limit_direction='both')\n        _df['isOutlier'] = _df['isOutlier'].interpolate(method='index', limit_direction='both')\n\n        time_df = pd.concat([time_df, _df])\n\n    time_df = time_df.sort_values(['phone', 'millisSinceGpsEpoch']).reset_index(drop=True)\n    return time_df\n\n\ndef kalmanfilter_interpolate(input_df: pd.DataFrame, base_df: pd.DataFrame, params: list):\n    _index = (input_df['target']==1)\n    _df = apply_kf_smoothing(\n        input_df.loc[_index].reset_index(drop=True),\n        params\n    )\n    input_df.loc[_index, ['latDeg', 'lngDeg']] = _df[['latDeg', 'lngDeg']].values\n    \n    input_df.loc[input_df['target']==0, 'latDeg'] = np.nan\n    input_df.loc[input_df['target']==0, 'lngDeg'] = np.nan\n    \n    input_df.index = input_df['millisSinceGpsEpoch'].values\n    for phone in input_df['phone'].unique():\n        _index = (input_df['phone']==phone)\n        input_df.loc[_index, 'latDeg'] = input_df.loc[_index, 'latDeg'].interpolate(method='index', limit_direction='both')\n        input_df.loc[_index, 'lngDeg'] = input_df.loc[_index, 'lngDeg'].interpolate(method='index', limit_direction='both')\n    \n    base_df = base_df.drop(columns=['latDeg', 'lngDeg', 'isOutlier'], axis=1)\n    output_df = base_df.merge(input_df, on=['phone', 'millisSinceGpsEpoch'])\n    return output_df\n\n\n###KALMAN FILTER###\ndef make_shifted_matrix(vec):\n    matrix = []\n    size = len(vec)\n    for i in range(size):\n        row = [0] * i + vec[:size-i]\n        matrix.append(row)\n    return np.array(matrix)\n\n\ndef make_state_vector(T, size):\n    vector = [1, 0]\n    step = 2\n    for i in range(size - 2):\n        if i % 2 == 0:\n            vector.append(T)\n            T *= T / step\n            step += 1\n        else:\n            vector.append(0)\n    return vector\n\n\ndef make_noise_vector(noise, size):\n    noise_vector = []\n    for i in range(size):\n        if i > 0 and i % 2 == 0:\n            noise *= 0.5\n        noise_vector.append(noise)\n    return noise_vector\n\n\ndef make_kalman_filter(T, size, noise, obs_noise):\n    vec = make_state_vector(T, size)\n    state_transition = make_shifted_matrix(vec)\n    process_noise = np.diag(make_noise_vector(noise, size)) + np.ones(size) * 1e-9\n    observation_model = np.array([[1] + [0] * (size - 1), [0, 1] + [0] * (size - 2)])\n    observation_noise = np.diag([obs_noise] * 2) + np.ones(2) * 1e-9\n    kf = simdkalman.KalmanFilter(\n            state_transition = state_transition,\n            process_noise = process_noise,\n            observation_model = observation_model,\n            observation_noise = observation_noise)\n    return kf\n\n\ndef apply_kf_smoothing(df, params):     \n    for name in df['phone'].unique():  \n        df.loc[df.isOutlier>0.5, 'latDeg'] = np.nan\n        df.loc[df.isOutlier>0.5, 'lngDeg'] = np.nan\n        \n    T, half_size, noise, obs_noise = params\n    size = half_size * 2\n    kf = make_kalman_filter(T, size, noise, obs_noise)\n    \n    unique_paths = df['phone'].unique()\n    for phone in tqdm(unique_paths):\n        _index = (df['phone']==phone)\n        data = df.loc[_index,['latDeg', 'lngDeg']].to_numpy()\n        data = data.reshape(1, len(data), 2)\n        smoothed = kf.smooth(data)\n        df.loc[_index, 'latDeg'] = smoothed.states.mean[0, :, 0]\n        df.loc[_index, 'lngDeg'] = smoothed.states.mean[0, :, 1]\n    \n    return df","metadata":{"execution":{"iopub.status.busy":"2021-08-05T09:18:18.511237Z","iopub.execute_input":"2021-08-05T09:18:18.511558Z","iopub.status.idle":"2021-08-05T09:18:18.54512Z","shell.execute_reply.started":"2021-08-05T09:18:18.511525Z","shell.execute_reply":"2021-08-05T09:18:18.543737Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"train_timedf = interpolate_beforekm(train_base)\ntest_timedf = interpolate_beforekm(test_base)\n\nparams = [1.5, 2, 1.3376883684997819e-07, 9.861453983492513e-07]\ntrainkm_df = kalmanfilter_interpolate(train_timedf, train_base, params)\ntestkm_df = kalmanfilter_interpolate(test_timedf, test_base, params)\n\ndel train_timedf\ndel test_timedf\ndel params","metadata":{"execution":{"iopub.status.busy":"2021-08-05T09:18:18.546858Z","iopub.execute_input":"2021-08-05T09:18:18.547515Z","iopub.status.idle":"2021-08-05T09:21:01.368316Z","shell.execute_reply.started":"2021-08-05T09:18:18.547464Z","shell.execute_reply":"2021-08-05T09:21:01.367029Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"trainkm_df = check_meter(trainkm_df)","metadata":{"execution":{"iopub.status.busy":"2021-08-05T09:21:01.369962Z","iopub.execute_input":"2021-08-05T09:21:01.37047Z","iopub.status.idle":"2021-08-05T09:21:10.467351Z","shell.execute_reply.started":"2021-08-05T09:21:01.37041Z","shell.execute_reply":"2021-08-05T09:21:10.465838Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# collectionNameMean","metadata":{}},{"cell_type":"code","source":"def collectionNamemean(df,cols:list,weight_col=None):\n    df_copy = df.copy()\n\n    if weight_col != None:\n        tmp_base = pd.pivot_table(data=df_copy,index=['collectionName','millisSinceGpsEpoch'],columns=['phoneName'],values=cols+[weight_col])\\\n        .reset_index(level='collectionName').groupby('collectionName')\\\n        .apply(lambda x:x.interpolate(limit_area='inside',method='index'))\n        for col in cols:\n            tmp_base[col] = tmp_base[col] * tmp_base[weight_col]\n            tmp_base[col] = tmp_base[col].values / tmp_base[weight_col].sum(axis=1).values.reshape(-1,1)\n\n        meta_df = tmp_base[cols[0]].sum(axis=1,skipna=True).reset_index().copy()\n        meta_df.columns = ['millisSinceGpsEpoch',cols[0]]\n\n        if len(cols) >=2:\n            for i in range(1,len(cols)):\n                meta_df[cols[i]] = tmp_base[cols[i]].sum(axis=1,skipna=True).reset_index()[0]\n\n    else:\n        tmp_base = pd.pivot_table(data=df_copy,index=['collectionName','millisSinceGpsEpoch'],columns=['phoneName'],values=cols)\\\n        .reset_index(level='collectionName').groupby('collectionName')\\\n        .apply(lambda x:x.interpolate(limit_area='inside',method='index'))\n    \n        meta_df = tmp_base[cols[0]].mean(axis=1,skipna=True).reset_index().copy()\n        meta_df.columns = ['millisSinceGpsEpoch',cols[0]]\n    \n        if len(cols) >=2:\n            for i in range(1,len(cols)):\n                meta_df[cols[i]] = tmp_base[cols[i]].mean(axis=1,skipna=True).reset_index()[0]\n\n    output_df = pd.merge(df_copy.drop(columns=cols),meta_df,how='left',on=['millisSinceGpsEpoch'])\n    return output_df\n","metadata":{"execution":{"iopub.status.busy":"2021-08-05T09:21:10.469449Z","iopub.execute_input":"2021-08-05T09:21:10.470165Z","iopub.status.idle":"2021-08-05T09:21:10.486622Z","shell.execute_reply.started":"2021-08-05T09:21:10.470114Z","shell.execute_reply":"2021-08-05T09:21:10.485537Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"trainkm_df = trainkm_df.sort_values([\"collectionName\", \"phoneName\", \"millisSinceGpsEpoch\"])\ntrainkm_df = trainkm_df.reset_index(drop=True)\n\ntestkm_df = testkm_df.sort_values([\"phone\", \"millisSinceGpsEpoch\"])\ntestkm_df = testkm_df.reset_index(drop=True)\n\nsatelite_train = pd.read_csv(\"../input/satelite/train_pseudorange.csv\").fillna(np.inf)\nrange_colums = [f\"pseudoranges_sigma_{i}\" for i in range(100)]\ntrainkm_df[\"acc_satelite\"] = 0\nfor c in range_colums:\n    trainkm_df[\"acc_satelite\"] += 1 / satelite_train[c]**2\n\nsatelite_test = pd.read_csv(\"../input/satelite/test_pseudorange.csv\")\nsatelite_test = satelite_test.sort_values([\"phone\", \"millisSinceGpsEpoch\"])\nsatelite_test = satelite_test.reset_index(drop=True).fillna(np.inf)\nrange_colums = [f\"pseudoranges_sigma_{i}\" for i in range(100)]\ntestkm_df[\"acc_satelite\"] = 0\nfor c in range_colums:\n    testkm_df[\"acc_satelite\"] += 1 / satelite_test[c]**2\n\ntrainmean_df = collectionNamemean(trainkm_df, ['latDeg','lngDeg'], 'acc_satelite')\ntestmean_df = collectionNamemean(testkm_df, ['latDeg','lngDeg'], 'acc_satelite')","metadata":{"execution":{"iopub.status.busy":"2021-08-05T09:21:10.488466Z","iopub.execute_input":"2021-08-05T09:21:10.489036Z","iopub.status.idle":"2021-08-05T09:21:31.337041Z","shell.execute_reply.started":"2021-08-05T09:21:10.488992Z","shell.execute_reply":"2021-08-05T09:21:31.335613Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"trainmean_df = check_meter(trainmean_df)","metadata":{"execution":{"iopub.status.busy":"2021-08-05T09:21:31.342718Z","iopub.execute_input":"2021-08-05T09:21:31.343167Z","iopub.status.idle":"2021-08-05T09:21:39.827048Z","shell.execute_reply.started":"2021-08-05T09:21:31.343124Z","shell.execute_reply":"2021-08-05T09:21:39.82496Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# back shift","metadata":{}},{"cell_type":"code","source":"def get_Dm_back(lat1, lon1, lat2, lon2, D):\n    L = calc_haversine(lat1, lon1, lat2, lon2)\n    if L == 0:\n        return lat2, lon2\n    while abs(D - L) > 0.0001:\n        L = calc_haversine(lat1, lon1, lat2, lon2)\n        lat1, lon1 = (lat2 - D/L * (lat2 - lat1), lon2 - D/L * (lon2 - lon1))\n    return lat1, lon1\n\n\ndef backshift(input_df: pd.DataFrame, d=0.42):\n    dfs = []\n    for pName, phone_df in tqdm(input_df.groupby([\"phone\"])):\n        phone_df = phone_df.reset_index(drop=True).copy()\n        phone_df[\"latDeg_moved\"] = phone_df.latDeg\n        phone_df[\"lngDeg_moved\"] = phone_df.lngDeg\n        phone_df.loc[1:, \"latDeg_moved\"], phone_df.loc[1:, \"lngDeg_moved\"] = np.vectorize(get_Dm_back)(\n            phone_df.latDeg[:-1],\n            phone_df.lngDeg[:-1],\n            phone_df.latDeg[1:],\n            phone_df.lngDeg[1:],\n            d\n        )\n        dfs.append(phone_df.copy()) \n\n    output_df = pd.concat(dfs).reset_index(drop=True)\n    output_df = output_df.drop(columns=['latDeg', 'lngDeg'], axis=1)\n    output_df = output_df.rename(columns={'latDeg_moved': 'latDeg', 'lngDeg_moved': 'lngDeg'})\n    \n    return output_df","metadata":{"execution":{"iopub.status.busy":"2021-08-05T09:21:39.829045Z","iopub.execute_input":"2021-08-05T09:21:39.82987Z","iopub.status.idle":"2021-08-05T09:21:39.843252Z","shell.execute_reply.started":"2021-08-05T09:21:39.829813Z","shell.execute_reply":"2021-08-05T09:21:39.841978Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# stop mean\ntrainmean_df = stopmean(trainmean_df)\ntestmean_df = stopmean(testmean_df)\n\n# back shift\ntrainshift_df = backshift(trainmean_df)\ntestshift_df = backshift(testmean_df)","metadata":{"execution":{"iopub.status.busy":"2021-08-05T09:21:39.84494Z","iopub.execute_input":"2021-08-05T09:21:39.845651Z","iopub.status.idle":"2021-08-05T09:21:57.544722Z","shell.execute_reply.started":"2021-08-05T09:21:39.845606Z","shell.execute_reply":"2021-08-05T09:21:57.543686Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"trainshift_df = check_meter(trainshift_df)","metadata":{"execution":{"iopub.status.busy":"2021-08-05T09:21:57.546305Z","iopub.execute_input":"2021-08-05T09:21:57.546793Z","iopub.status.idle":"2021-08-05T09:22:11.692541Z","shell.execute_reply.started":"2021-08-05T09:21:57.546753Z","shell.execute_reply":"2021-08-05T09:22:11.691397Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# snap to grid & submission","metadata":{}},{"cell_type":"code","source":"def snap_to_grid(df,threshhold):\n    dfs = []\n    \n    no_area = 0\n    if 'area' not in df.columns:\n        no_area = 1\n        df['area'] = df['phone'].apply(lambda x:x.split('-')[-2])\n    \n    for area,gdf in tqdm(df.groupby('area')):\n        if area in ['MTV', 'SF', 'RWC', 'SVL', 'SJC']:\n            gdf['neighbored'] = knn_dict[area].predict(gdf[['latDeg','lngDeg']]).astype(int)\n        else:\n            gdf['neighbored'] = 0 #testにしかないareaに関してはとりあえず0で置く,thresholdで消されるからOK\n        dfs.append(gdf['neighbored'])\n    neighbored = pd.concat(dfs).sort_index().values\n    \n    df[[\"latDeg_near\", \"lngDeg_near\", \"millisSinceGpsEpoch_near\"]] = wps.loc[neighbored, [\"latDeg\", \"lngDeg\", \"millisSinceGpsEpoch\"]].values\n    df[\"d_near\"] = np.vectorize(calc_haversine)(df.latDeg, df.lngDeg, df.latDeg_near, df.lngDeg_near)\n    df.loc[df['d_near']<threshhold,'latDeg'] = df.loc[df['d_near']<threshhold,'latDeg_near']\n    df.loc[df['d_near']<threshhold,'lngDeg'] = df.loc[df['d_near']<threshhold,'lngDeg_near']\n\n    if no_area:\n        df.drop(columns=['area'],inplace=True)\n    \n    return df.drop(columns=['latDeg_near','lngDeg_near','millisSinceGpsEpoch_near','d_near']).copy()\n\n\nwps = pd.read_csv('../input/snap-data/train_waypoint.csv')\nwith open('../input/snap-data/knn_dict.pkl', 'rb') as dic:\n    knn_dict = pickle.load(dic)","metadata":{"execution":{"iopub.status.busy":"2021-08-05T09:23:01.177816Z","iopub.execute_input":"2021-08-05T09:23:01.178223Z","iopub.status.idle":"2021-08-05T09:23:25.829516Z","shell.execute_reply.started":"2021-08-05T09:23:01.178159Z","shell.execute_reply":"2021-08-05T09:23:25.828086Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"sjc_train_index = trainshift_df['collectionName'].map(lambda x: 'SJC' in x)\ntrainsnap_df = pd.concat([\n    snap_to_grid(trainshift_df[sjc_train_index], 100),\n    snap_to_grid(trainshift_df[~sjc_train_index], 2)\n]).reset_index(drop=True)\n\nsjc_test_index = testshift_df['collectionName'].map(lambda x: 'SJC' in x)\ntestsnap_df = pd.concat([\n    snap_to_grid(testshift_df[sjc_test_index], 100),\n    snap_to_grid(testshift_df[~sjc_test_index], 2)\n]).reset_index(drop=True)","metadata":{"execution":{"iopub.status.busy":"2021-08-05T09:25:28.546527Z","iopub.execute_input":"2021-08-05T09:25:28.547038Z","iopub.status.idle":"2021-08-05T09:26:06.544142Z","shell.execute_reply.started":"2021-08-05T09:25:28.546995Z","shell.execute_reply":"2021-08-05T09:26:06.54303Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"trainsnap_df = check_meter(trainsnap_df)","metadata":{"execution":{"iopub.status.busy":"2021-08-05T09:26:06.546118Z","iopub.execute_input":"2021-08-05T09:26:06.54662Z","iopub.status.idle":"2021-08-05T09:26:16.356085Z","shell.execute_reply.started":"2021-08-05T09:26:06.546575Z","shell.execute_reply":"2021-08-05T09:26:16.354966Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"sub = sample.drop(columns=['latDeg', 'lngDeg'], axis=1).merge(\n        testsnap_df[['phone', 'millisSinceGpsEpoch', 'latDeg', 'lngDeg']],\n        on=['phone', 'millisSinceGpsEpoch']\n)\nsub.to_csv('submission.csv', index=False)","metadata":{"execution":{"iopub.status.busy":"2021-08-05T09:26:19.784732Z","iopub.execute_input":"2021-08-05T09:26:19.785109Z","iopub.status.idle":"2021-08-05T09:26:19.844966Z","shell.execute_reply.started":"2021-08-05T09:26:19.785074Z","shell.execute_reply":"2021-08-05T09:26:19.843888Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"visualize_trafic(trainshift_df, color='collectionName')","metadata":{"execution":{"iopub.status.busy":"2021-08-05T09:26:20.628025Z","iopub.execute_input":"2021-08-05T09:26:20.628436Z","iopub.status.idle":"2021-08-05T09:26:22.841511Z","shell.execute_reply.started":"2021-08-05T09:26:20.628398Z","shell.execute_reply":"2021-08-05T09:26:22.840307Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"visualize_trafic(testshift_df, color='collectionName')","metadata":{"execution":{"iopub.status.busy":"2021-08-05T09:26:22.843348Z","iopub.execute_input":"2021-08-05T09:26:22.843805Z","iopub.status.idle":"2021-08-05T09:26:23.985225Z","shell.execute_reply.started":"2021-08-05T09:26:22.843761Z","shell.execute_reply":"2021-08-05T09:26:23.982428Z"},"trusted":true},"execution_count":null,"outputs":[]}]}