{"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":"## 5th Place Post-Processing\nThis is our team's post-processing pipeline for [Saito](https://www.kaggle.com/saitodevel01)'s baseline\n\n(Some of the input from Saito is private, but see the following notebook for that part.)\n* https://www.kaggle.com/saitodevel01/gsdc-optimization-based-smoothing-1st-version\n* https://www.kaggle.com/saitodevel01/gsdc-bias-eda\n* https://www.kaggle.com/saitodevel01/dsdc-unified-post-processing\n* https://www.kaggle.com/saitodevel01/gsdc-bias-correction\n* https://www.kaggle.com/saitodevel01/gsdc-vehicle-speed-estimation-by-doppler-shift","metadata":{}},{"cell_type":"code","source":"!pip install pandarallel","metadata":{"papermill":{"duration":10.868465,"end_time":"2021-07-19T08:08:06.246487","exception":false,"start_time":"2021-07-19T08:07:55.378022","status":"completed"},"tags":[],"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\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\nfrom pandarallel import pandarallel\npandarallel.initialize()\n\nimport scipy.sparse\nimport scipy.sparse.linalg\nimport multiprocessing\n\nfrom sklearn.neighbors import KNeighborsClassifier\nfrom utils import * \nimport simdkalman\n\nimport saito_io_f_v2 as io_f\nimport constants \nimport transform \nimport qpsolver\nimport signal_f_v2 as signal_f\nimport map_matching\nimport design_filter\nfrom scipy.interpolate import InterpolatedUnivariateSpline","metadata":{"_cell_guid":"be16aef1-c9d1-4b10-b5e2-682a68aac86d","_uuid":"e019fbb5-3f3b-495c-8bcb-06f2488a3375","collapsed":false,"jupyter":{"outputs_hidden":false},"papermill":{"duration":57.220809,"end_time":"2021-07-19T08:09:03.500753","exception":false,"start_time":"2021-07-19T08:08:06.279944","status":"completed"},"tags":[],"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def area_prediction():    \n    BASE_DIR = Path('../input/google-smartphone-decimeter-challenge')\n\n    train_base = pd.read_csv(BASE_DIR / 'baseline_locations_train.csv')\n    train_base = train_base.sort_values([\n        \"collectionName\", \"phoneName\", \"millisSinceGpsEpoch\"\n    ]).reset_index(drop=True)\n    train_base['area'] = train_base['collectionName'].map(lambda x: x.split('-')[4])\n    \n    test_base = pd.read_csv(BASE_DIR / 'baseline_locations_test.csv')\n    test_base = test_base.sort_values([\n        \"collectionName\", \"phoneName\", \"millisSinceGpsEpoch\"\n    ]).reset_index(drop=True)\n    test_base['area'] = test_base['collectionName'].map(lambda x: x.split('-')[4])\n\n    train_name = np.array(sorted(path.split('/')[-1] for path in glob(f'{BASE_DIR}/train/*')))\n    train_highway  = train_name[np.array([1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21]) - 1]\n    train_tree     = train_name[np.array([22,23,25,26,28]) - 1]\n    train_downtown = train_name[np.array([24,27,29]) - 1]\n\n    train_base['area_target'] = -1\n    train_base.loc[train_base['collectionName'].isin(train_highway),  'area_target'] = 0\n    train_base.loc[train_base['collectionName'].isin(train_tree),     'area_target'] = 1\n    train_base.loc[train_base['collectionName'].isin(train_downtown), 'area_target'] = 2\n\n    def processing_downtown(input_df: pd.DataFrame, is_train=False):\n        output_df = input_df.groupby('collectionName')[['latDeg', 'lngDeg']].std()\n        if is_train:\n            output_df = output_df.merge(\n                input_df.groupby('collectionName')[['area_target']].first(),\n                on='collectionName')\n        output_df = output_df.merge(\n            input_df.groupby('collectionName')['area'].first(),\n            on='collectionName')\n        output_df = output_df.merge(\n            input_df.groupby('collectionName')['phoneName'].unique().apply(list),\n            on='collectionName')\n        return output_df\n\n    train = processing_downtown(train_base, is_train=True)\n    train['downtown_target'] = (train['area_target']==2).astype(int)\n\n    downtown_model_knn = KNeighborsClassifier(n_neighbors=1)\n    downtown_model_knn.fit(\n        train[['latDeg', 'lngDeg']],\n        train['downtown_target'],\n    )\n\n    def processing_highway_tree(input_df: pd.DataFrame, is_train=False):\n        output_df = input_df.groupby('collectionName')[['latDeg', 'lngDeg']].min()\n        if is_train:\n            output_df = output_df.merge(\n                input_df.groupby('collectionName')[['area_target']].first(),\n                on='collectionName')\n        output_df = output_df.merge(\n            input_df.groupby('collectionName')['area'].first(),\n            on='collectionName')\n        output_df = output_df.merge(\n            input_df.groupby('collectionName')['phoneName'].unique().apply(list),\n            on='collectionName')\n        return output_df\n\n    train = processing_highway_tree(train_base, is_train=True)\n\n    highway_tree_model_knn = KNeighborsClassifier(n_neighbors=1)\n    highway_tree_model_knn.fit(\n        train.loc[train['area_target']!=2, ['latDeg', 'lngDeg']],\n        train.loc[train['area_target']!=2, 'area_target'],\n    )\n\n    def predict_area(test_base):\n        test_base = test_base.copy()\n        test_base = test_base.sort_values([\n            \"collectionName\", \"phoneName\", \"millisSinceGpsEpoch\"\n        ]).reset_index(drop=True)\n        test_base['area'] = test_base['collectionName'].map(lambda x: x.split('-')[4])\n\n        test = processing_downtown(test_base)\n        downtown_pred = downtown_model_knn.predict(test[['latDeg', 'lngDeg']])\n\n        test = processing_highway_tree(test_base)\n        test.loc[downtown_pred==1, 'area_pred'] = 2\n        pred = highway_tree_model_knn.predict(test.loc[test['area_pred'].isnull(), ['latDeg', 'lngDeg']])\n        test.loc[test['area_pred'].isnull(), 'area_pred'] = pred\n        test['area_pred'] = test['area_pred'].astype(int)\n        test['collectionName'] = test.index\n\n        test_highway  = []\n        test_tree     = []\n        test_downtown = []\n        for collection, area_pred in test[['collectionName', 'area_pred']].itertuples(index=False):\n            if area_pred == 0:\n                test_highway.append(collection)\n            elif area_pred == 1:\n                test_tree.append(collection)\n            else:\n                test_downtown.append(collection)\n        return (test_highway, test_tree, test_downtown)\n    return train_highway, train_tree, train_downtown, predict_area(test_base)\nTRAIN_HIGHWAY, TRAIN_TREEWAY, TRAIN_DOWNTOWN, \\\n    (TEST_HIGHWAY, TEST_TREEWAY, TEST_DOWNTOWN) = area_prediction()","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Generalized functions of LightGBM\ndef fit_lgbm(X, y, train_df, params: dict=None, verbose=100, seed: int=42, N_SPLITS: int=5):\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        model.fit(x_train, y_train, \n            eval_set=[(x_valid, y_valid)],  \n            early_stopping_rounds=verbose, \n            eval_metric='logloss',\n            verbose=0)\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":{"papermill":{"duration":0.051905,"end_time":"2021-07-19T08:09:16.604119","exception":false,"start_time":"2021-07-19T08:09:16.552214","status":"completed"},"tags":[],"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# KnnHeight\nMap the latitude and longitude from the ground truth to the correct altitude. <br>\nIf there is a discrepancy between the predicted altitude and the output of this model, the accuracy of the satellite positioning can be expected to be poor.<br>","metadata":{}},{"cell_type":"code","source":"class KnnHeight:\n    name = \"Knn\"\n    isPrep = True\n    \n    def __init__(self):\n        pass\n    \n    def main(self, train, test):\n        knn_model = self.fit_knn_height(\n            train[['t_latDeg', 't_lngDeg']],\n            train['t_heightAboveWgs84EllipsoidM'])\n\n        train['heightKNN'] = self.pred_knn_height(train[['latDeg', 'lngDeg']], knn_model)\n        test['heightKNN']  = self.pred_knn_height(test[['latDeg', 'lngDeg']] , knn_model)\n\n        train['heightDiff'] = np.abs(train['heightAboveWgs84EllipsoidM'] - train['heightKNN'])\n        test['heightDiff']  = np.abs(test['heightAboveWgs84EllipsoidM']  - test['heightKNN'])\n        \n        return train, test\n    \n    def fit_knn_height(self, 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    def pred_knn_height(self, X: pd.DataFrame, model):\n        return model.predict(X.values)","metadata":{"papermill":{"duration":1.100312,"end_time":"2021-07-19T08:09:16.444892","exception":false,"start_time":"2021-07-19T08:09:15.34458","status":"completed"},"tags":[],"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Outlier Detection\nPredict the probability of being an outlier using the relative coordinates (calculated from the absolute coordinates) of the surrounding 50 seconds as a feature.<br>\nIn the case of Downtown, the distance to the ground truth is added to the feature value.<br>\nIf the absolute position prediction is outputting the altitude, add the difference from the kNN output.<br>","metadata":{}},{"cell_type":"code","source":"class OutlierDetection:\n    name = \"OutlierDetection\"\n    isPrep = True\n    \n    def __init__(self, sjc_loss_threshold, other_loss_threshold):\n        self.sjc_loss_threshold   = sjc_loss_threshold\n        self.other_loss_threshold = other_loss_threshold\n        \n    def main(self, train, test):\n        # distance from ground truth\n        train_snapped = snap_to_grid(train, 1e10)\n        test_snapped = snap_to_grid(test, 1e10)\n        train['dist'] = calc_dist(train[['latDeg','lngDeg']], train_snapped[['latDeg','lngDeg']])\n        test['dist']  = calc_dist(test[['latDeg','lngDeg']] , test_snapped[['latDeg','lngDeg']])\n        del train_snapped,test_snapped\n\n        # add area\n        train['area'] = train['phone'].apply(lambda s : s.split('-')[4])\n        test['area'] = test['phone'].apply(lambda s : s.split('-')[4])\n        \n        params = {\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        \n        # Downtown\n        loss_threshold = self.sjc_loss_threshold\n        visibility     = 22\n        train['t_isOutlier'] = (train.meter > loss_threshold).astype(int)\n        N_SPLITS = 3\n        use_train_index = train.collectionName.isin(TRAIN_DOWNTOWN)\n        use_test_index = test.collectionName.isin(TEST_DOWNTOWN)\n        oof, models = fit_lgbm(self.processing_sjc(train[use_train_index].reset_index(drop=True), visibility), \n                               train.loc[use_train_index, 't_isOutlier'], \n                               train_df=train[use_train_index],\n                               params=params, N_SPLITS=N_SPLITS)\n        pred = predict_lgbm(models, self.processing_sjc(test[use_test_index].reset_index(drop=True), visibility))\n        train.loc[use_train_index, 'outlier_rate'] = oof\n        test.loc[use_test_index, 'outlier_rate'] = pred\n        print('score', accuracy_score((oof>0.5).astype(int), train.loc[use_train_index, 't_isOutlier']))\n        print(confusion_matrix((oof>0.5).astype(int), train.loc[use_train_index, 't_isOutlier']))\n\n        # others\n        loss_threshold = self.other_loss_threshold\n        visibility     = 26 \n        train['t_isOutlier'] = (train.meter > loss_threshold).astype(int)\n        N_SPLITS = 5\n        use_train_index = ~train.collectionName.isin(TRAIN_DOWNTOWN)\n        use_test_index = ~test.collectionName.isin(TEST_DOWNTOWN)\n        oof, models = fit_lgbm(self.processing_notsjc(train[use_train_index].reset_index(drop=True), visibility), \n                               train.loc[use_train_index, 't_isOutlier'], \n                               train_df=train[use_train_index],\n                               params=params, N_SPLITS=N_SPLITS)\n        pred = predict_lgbm(models, self.processing_notsjc(test[use_test_index].reset_index(drop=True), visibility))\n        train.loc[use_train_index, 'outlier_rate'] = oof\n        test.loc[use_test_index, 'outlier_rate'] = pred\n        print('score', accuracy_score((oof>0.5).astype(int), train.loc[use_train_index, 't_isOutlier']))\n        print(confusion_matrix((oof>0.5).astype(int), train.loc[use_train_index, 't_isOutlier']))\n        \n        return train, test        \n        \n\n    def processing_sjc(self, input_df: pd.DataFrame, L=25):\n        output_df = pd.DataFrame()\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                input_df.groupby('phone')[['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    def processing_notsjc(self, input_df: pd.DataFrame, L=25):\n        output_df = pd.DataFrame()\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":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Interpolate Outliers by Relpositions\nWe change the threshold and the direction of interpolation as follows, and after 20 predictions, we take a weighted average using the threshold.<br>\n\n1. take 10 threshold values at equal intervals within 0.1-0.5.\n    1. where the threshold is exceeded, it is considered as an outlier and replaced by nan.\n    2. accumulate the relative coordinates from the nearest non-nan point and prepare one new absolute coordinate. (There are two ways to do this: forward and reverse.)\n2. replace the predicted value with a weighted average of these 20 with \"1/threshold\".","metadata":{}},{"cell_type":"code","source":"class InterpolateByRelpos:\n    name = \"InterpolateByRelpos\"\n    isPrep = False\n    \n    def __init__(self, sjc_start=0.1, sjc_end=0.5, other_start=0.1, other_end=0.5):\n        # downtown\n        self.sjc_start = sjc_start\n        self.sjc_end = sjc_end\n        \n        # others\n        self.other_start = other_start\n        self.other_end = other_end\n    \n    def main(self, train, test):\n        train['area'] = train['phone'].apply(lambda s : s.split('-')[4])\n        test['area'] = test['phone'].apply(lambda s : s.split('-')[4])\n        \n        train = pd.concat([\n            self.interpolate_by_relpos(train[train.collectionName.isin(TRAIN_DOWNTOWN)], self.sjc_start, self.sjc_end),\n            self.interpolate_by_relpos(train[~train.collectionName.isin(TRAIN_DOWNTOWN)], self.other_start, self.other_end),\n        ]).sort_values([\"phone\", \"millisSinceGpsEpoch\"]).reset_index(drop=True)\n        \n        \n        test_sjc    =  test.collectionName.isin(TEST_DOWNTOWN)\n        test_nonsjc = ~test.collectionName.isin(TEST_DOWNTOWN)\n        test = pd.concat([\n            self.interpolate_by_relpos(test[test_sjc], self.sjc_start, self.sjc_end),\n            self.interpolate_by_relpos(test[test_nonsjc], self.other_start, self.other_end),\n        ]).sort_values([\"phone\", \"millisSinceGpsEpoch\"]).reset_index(drop=True)\n        \n        train[\"isOutlier\"] = 0\n        test[\"isOutlier\"] = 0\n        \n        return train, test\n        \n        \n    def interpolate_by_relpos(self, df, start, end):                                   \n        dfs = []\n        for phone, phone_df in tqdm(df.groupby(\"phone\")):\n            phone_df = phone_df.reset_index(drop=True)\n\n            phone_lats = []\n            phone_lngs = []\n\n            THRESH_LIST = np.linspace(start, end, 10)\n            for THRESH in THRESH_LIST:\n                \"\"\"\n                 Interpolate in the forward direction \n                \"\"\"\n                acsend_df = phone_df.copy()\n                acsend_df.loc[(acsend_df.index > 0) & (acsend_df.outlier_rate > THRESH), [\"latDeg\", \"lngDeg\"]] = np.nan\n                lats = acsend_df.latDeg.tolist()\n                lngs = acsend_df.lngDeg.tolist()\n                d_lats = acsend_df.delta_latDeg.tolist()\n                d_lngs = acsend_df.delta_lngDeg.tolist()\n                for i, (lat, lng, d_lat, d_lng) in enumerate(zip(lats, lngs, d_lats, d_lngs)):\n                    if lat != lat and i > 0:\n                        lats[i] = lats[i-1] + d_lat\n                        lngs[i] = lngs[i-1] + d_lng\n                acsend_df[\"latDeg\"] = lats\n                acsend_df[\"lngDeg\"] = lngs        \n\n\n                \"\"\"\n                 Interpolate in the reverse direction\n                \"\"\"\n                reverse_df = phone_df.copy()\n                reverse_df.loc[(reverse_df.index < len(reverse_df)-1) & (reverse_df.outlier_rate > THRESH), [\"latDeg\", \"lngDeg\"]] = np.nan\n                lats = reverse_df.latDeg.tolist()\n                lngs = reverse_df.lngDeg.tolist()\n                d_lats = reverse_df.delta_latDeg.tolist()\n                d_lngs = reverse_df.delta_lngDeg.tolist()\n                for i in range(len(lats)-2, -1, -1):    \n                    if lats[i] != lats[i] and i + 1 < len(lats):\n                        lats[i] = lats[i+1] - d_lats[i+1]\n                        lngs[i] = lngs[i+1] - d_lngs[i+1]\n\n                reverse_df[\"latDeg\"] = lats\n                reverse_df[\"lngDeg\"] = lngs        \n\n                phone_lats.append((acsend_df.latDeg + reverse_df.latDeg) / 2)\n                phone_lngs.append((acsend_df.lngDeg + reverse_df.lngDeg) / 2)\n\n            \"\"\"\n             average 20 cordinates\n            \"\"\"\n            lats = np.average(phone_lats, axis=0, weights=1/THRESH_LIST)\n            lngs = np.average(phone_lngs, axis=0, weights=1/THRESH_LIST)\n\n            phone_df[\"latDeg\"] = lats\n            phone_df[\"lngDeg\"] = lngs\n\n            dfs.append(phone_df)\n        df = pd.concat(dfs)\n        \n        return df.reset_index(drop=True)","metadata":{"papermill":{"duration":173.029017,"end_time":"2021-07-19T08:13:32.727751","exception":false,"start_time":"2021-07-19T08:10:39.698734","status":"completed"},"tags":[],"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# StopMean\nAs in outlier detection, the stop point is predicted from relative coordinates.","metadata":{}},{"cell_type":"code","source":"class StopMean:\n    name = \"StopMean\"\n    isPrep = False\n    \n    def __init__(self):\n        pass\n    \n    def main(self, train, test):\n        train['target_stop'] = (train['speedMps']==0).astype(int)\n\n        params = {\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\n        N_SPLITS = 5\n        oof, models = fit_lgbm(self.processing_stopmean(train), train['target_stop'], train_df=train, params=params, N_SPLITS=N_SPLITS)\n        pred = predict_lgbm(models, self.processing_stopmean(test))\n\n        train['stop'] = (oof > 0.5).astype(int)\n        test['stop'] = (pred > 0.5).astype(int)\n\n        print('score', accuracy_score((oof>0.5).astype(int), train['target_stop']))\n        print(confusion_matrix((oof>0.5).astype(int), train['target_stop']))\n        \n        train = stopmean(train)\n        test = stopmean(test)\n        \n        return train, test\n        \n    def processing_stopmean(self, input_df: pd.DataFrame):\n        output_df = pd.DataFrame()\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        ","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Cost Minimization","metadata":{}},{"cell_type":"code","source":"class SaitoOptimization:\n    name = \"SaitoOptimization\"\n    isPrep = False\n    \n    def __init__(self, N_SNAP_TO_GRID_ITERATION):\n        self.params_highway = { 'sigma_u'  : 1.0,\n                       'sigma_p'  : 3.0,\n                       'sigma_a'  : 2.0 * 1e+5,\n                       'sigma_v'  : 4.0 * 1e+5,\n                       'sigma_d'  : 0.16 * 1e+5,\n                       'reject_p' : 7.0,   # [m]\n                       'reject_d' : 1.0,   # [m/s]\n                       'vmin'     : -0.05, # [m/s]\n                       'vmax'     : 50.0,  # [m/s]\n                       'Mi8_velocity_timeshift' : 0.46,\n                       'use_not_go_back_constraint' : False,\n                       'use_map'  : False,\n                      }\n        self.params_treeway = { 'sigma_u'  : 1.0,\n                       'sigma_p'  : 6.0,\n                       'sigma_a'  : 0.8 * 1e+5,\n                       'sigma_v'  : 3.0 * 1e+5,\n                       'sigma_d'  : 0.12 * 1e+5,\n                       'reject_p' : 12.0,  # [m]\n                       'reject_d' : 1.0,   # [m/s]\n                       'vmin'     : -0.05, # [m/s]\n                       'vmax'     : 50.0,  # [m/s]\n                       'Mi8_velocity_timeshift' : 0.30,\n                       'use_not_go_back_constraint' : False,\n                       'use_map'  : False,\n                      }\n        self.params_downtown = { 'sigma_u'  : 1.0,\n                        'sigma_p'  : 20.0,\n                        'sigma_a'  : 0.4 * 1e+5,\n                        'sigma_v'  : 1.0 * 1e+5,\n                        'sigma_d'  : 1.3 * 1e+5,\n                        'reject_p' : 20.0,  # [m]\n                        'reject_d' : 3.0,   # [m/s]\n                        'vmin'     : -0.05, # [m/s]\n                        'vmax'     : 50.0,  # [m/s]\n                        'Mi8_velocity_timeshift' : 0.0,\n                        'use_not_go_back_constraint' : False,\n                        'use_map'  : True,\n                        'threshold_distance_to_nearest_neighbor' : 8.0,\n                        'sigma_p_stage2' : 3.0,\n                        'num_stage2_iterations' : N_SNAP_TO_GRID_ITERATION,\n                       }\n        self.DT_X = 1.0\n    \n    def get_optimization_constants(self, base_df, velocity_df, sensor_df, params, use_sensor):\n        const = dict()\n        dt    = self.DT_X\n        TIME_y = base_df['Time'].values\n        TIME_d = velocity_df['Time'].values\n        N_y = TIME_y.shape[0]\n        N_d = TIME_d.shape[0]\n        N_x = int(np.ceil(np.max(TIME_y) / dt) + 1)\n        const['N_y'] = N_y\n        const['N_d'] = N_d\n        const['N_x'] = N_x\n\n        a = np.array([[1, dt, (1/2)*dt**2],\n                      [0,  1,  dt],\n                      [0,  0,  1]])\n        e3 = scipy.sparse.eye(3)\n        A = np.empty(shape=(2*(N_x-1), 2*N_x), dtype=np.object)\n        for i_x in range(N_x-1):\n            A[2*i_x  , 2*i_x  ] = a\n            A[2*i_x+1, 2*i_x+1] = a\n            A[2*i_x  , 2*i_x+2] = -e3\n            A[2*i_x+1, 2*i_x+3] = -e3\n        const['A'] = scipy.sparse.bmat(A, format='csr')\n\n        b = np.array([[(1/6)*dt**3,\n                       (1/2)*dt**2,\n                       dt]]).T\n        const['B'] = scipy.sparse.block_diag([b for _ in range(2*(N_x-1))], format='csr')\n\n        diag_R  = np.full(2*N_x - 2, params['sigma_u']**(-2) * dt)\n        const['R'] = scipy.sparse.spdiags(diag_R, [0], 2*N_x - 2, 2*N_x - 2, format='csc')\n\n        x_index  = np.floor(TIME_y / dt).astype(int)\n        alpha    = (TIME_y / dt) - x_index\n        coeff_y0 = 1 - 3*alpha**2 + 2*alpha**3\n        coeff_y1 =     3*alpha**2 - 2*alpha**3\n        coeff_v0 = dt * alpha * (alpha - 1)**2\n        coeff_v1 = dt * alpha**2 * (alpha - 1)\n        C = np.empty(shape=(2*N_y, 2*N_x), dtype=np.object)\n        for i_x in range(N_x):\n            C[0, 2*i_x  ] = scipy.sparse.coo_matrix((1, 3))\n            C[0, 2*i_x+1] = scipy.sparse.coo_matrix((1, 3))\n        for i_y in range(N_y):\n            i_x = x_index[i_y]\n            c_i = np.array([[coeff_y0[i_y], coeff_v0[i_y], 0]])\n            C[2*i_y,   2*i_x]   = c_i\n            C[2*i_y+1, 2*i_x+1] = c_i\n            if i_x < N_x - 1:\n                c_iplus = np.array([[coeff_y1[i_y], coeff_v1[i_y], 0]])\n                C[2*i_y,   2*i_x+2] = c_iplus\n                C[2*i_y+1, 2*i_x+3] = c_iplus\n        const['Cp_orig']  = scipy.sparse.bmat(C, format='csr')\n\n        diag_Lp = np.full(2*N_y, params['sigma_p']**(-2))\n        const['Lp_orig'] = scipy.sparse.spdiags(diag_Lp, [0], 2*N_y, 2*N_y, format='csr')\n        const['Yp_orig'] = base_df[['latDeg', 'lngDeg']].values.flatten()\n\n        BLH = transform.BLH(\n            lat=np.deg2rad(base_df['latDeg'].values),\n            lng=np.deg2rad(base_df['lngDeg'].values),\n            hgt=np.zeros(N_y),\n        )\n        DEG2RAD = np.pi / 180.0\n        J = transform.jacobian_BL_to_EN(BLH) * DEG2RAD\n        J = np.mean(J, axis=0)\n        J[0, 0] = 0\n        J[1, 1] = 0\n        JJ = scipy.sparse.block_diag([J, J], format='csr')\n        const['J'] = J\n\n        # ドップラ速度に関するパラメータ\n        x_index  = np.floor(TIME_d / dt).astype(int)\n        alpha    = (TIME_d / dt) - x_index\n        coeff_y0 = 1 - 3*alpha**2 + 2*alpha**3\n        coeff_y1 =     3*alpha**2 - 2*alpha**3\n        coeff_v0 = dt * alpha * (alpha - 1)**2\n        coeff_v1 = dt * alpha**2 * (alpha - 1)\n        C = np.empty(shape=(N_d, N_x), dtype=np.object)\n        for i_x in range(N_x):\n            C[0, i_x] = scipy.sparse.coo_matrix((2, 6))\n        for i_d in range(N_d):\n            i_x = x_index[i_d]\n            c = np.array([[0, coeff_y0[i_d], coeff_v0[i_d], 0, 0, 0],\n                          [0, 0, 0, 0, coeff_y0[i_d], coeff_v0[i_d]]])\n            C[i_d, i_x] = J @ c\n            if i_x < N_x - 1:\n                c = np.array([[0, coeff_y1[i_d], coeff_v1[i_d], 0, 0, 0],\n                              [0, 0, 0, 0, coeff_y1[i_d], coeff_v1[i_d]]])\n                C[i_d, i_x+1] = J @ c\n        const['Cd_orig']  = scipy.sparse.bmat(C, format='csr')\n\n        diag_Ld = np.full(2*N_d, params['sigma_d']**(-2))\n        const['Ld_orig'] = scipy.sparse.spdiags(diag_Ld, [0], 2*N_d, 2*N_d, format='csr')\n        const['Yd_orig'] = velocity_df[['v_east', 'v_north']].values.flatten()\n\n        if sensor_df is None:\n            const['use_sensor'] = False\n            const['use_inquality'] = False\n            return const\n\n        TIME_s = sensor_df['Time'].values\n        N_s = TIME_s.shape[0]\n        const['N_s'] = N_s\n        const['use_sensor'] = use_sensor\n        const['use_inquality'] = (use_sensor and params['use_not_go_back_constraint'])\n        x_index = np.round(TIME_s / dt).astype(int)\n        const['x_index_sensor'] = x_index\n        if not use_sensor:\n            return const\n\n        # 速度制約・速度コストに関するパラメータ\n        COS_TH = sensor_df['cos_th'].values\n        SIN_TH = sensor_df['sin_th'].values\n        CV = np.empty(shape=(N_s, N_x), dtype=np.object)\n        GV = np.empty(shape=(N_s, N_x), dtype=np.object)\n        cv = np.array([[0, 1, 0, 0, 0, 0],\n                       [0, 0, 0, 0, 1, 0]], dtype=np.float64)\n        for i_x in range(N_x):\n            CV[0, i_x] = scipy.sparse.coo_matrix((1, 6))\n            GV[0, i_x] = scipy.sparse.coo_matrix((1, 6))\n        for i_s in range(N_s):\n            i_x = x_index[i_s]\n            k = np.array([[SIN_TH[i_s], -COS_TH[i_s]]])\n            CV[i_s, i_x] = k @ J @ cv\n            k = np.array([[-COS_TH[i_s], -SIN_TH[i_s]]])\n            GV[i_s, i_x] = k @ J @ cv\n        const['Cv'] = scipy.sparse.bmat(CV, format='csr')\n        const['Gv'] = scipy.sparse.bmat(GV, format='csr')\n        const['hv'] = np.full((N_s, ), -params['vmin'])\n\n        diag_Lv = np.full(N_s, params['sigma_v']**(-2))\n        const['Lv'] = scipy.sparse.spdiags(diag_Lv, [0], N_s, N_s, format='csr')\n\n        # 加速度コストに関するパラメータ\n        DOT_V_COS_TH = sensor_df['dotV'] * sensor_df['cos_th'].values\n        DOT_V_SIN_TH = sensor_df['dotV'] * sensor_df['sin_th'].values\n        OMEGA = sensor_df['omega'].values\n        CA = np.empty(shape=(N_s, N_x), dtype=np.object)\n        ca = np.array([[0, 1, 0, 0, 0, 0],\n                       [0, 0, 0, 0, 1, 0],\n                       [0, 0, 1, 0, 0, 0],\n                       [0, 0, 0, 0, 0, 1]], dtype=np.float64)\n        for i_x in range(N_x):\n            CA[0, i_x] = scipy.sparse.coo_matrix((2, 6))\n        for i_s in range(N_s):\n            i_x = x_index[i_s]\n            k = np.array([[0,  OMEGA[i_s], 1, 0],\n                          [-OMEGA[i_s], 0, 0, 1]])\n            CA[i_s, i_x] = k @ JJ @ ca\n        const['Ca'] = scipy.sparse.bmat(CA, format='csr')\n        const['Ya'] = np.stack([DOT_V_COS_TH, DOT_V_SIN_TH], axis=1).flatten()\n\n        diag_La = np.full(2*N_s, params['sigma_a']**(-2))\n        const['La'] = scipy.sparse.spdiags(diag_La, [0], 2*N_s, 2*N_s, format='csr')\n\n        return const\n\n\n    def solve_QP(self, const, p_valid, d_valid):\n        A = const['A']\n        B = const['B']\n        R = const['R']\n        Cp_orig = const['Cp_orig']\n        Lp_orig = const['Lp_orig']\n        Yp_orig = const['Yp_orig']\n        Cd_orig = const['Cd_orig']\n        Ld_orig = const['Ld_orig']\n        Yd_orig = const['Yd_orig']\n\n        p_valid2 = np.stack([p_valid, p_valid], axis=1).flatten()\n        Cp = Cp_orig[p_valid2, :]\n        Lp = Lp_orig[np.ix_(p_valid2, p_valid2)]\n        Yp = Yp_orig[p_valid2]\n\n        d_valid2 = np.stack([d_valid, d_valid], axis=1).flatten()\n        Cd = Cd_orig[d_valid2, :]\n        Ld = Ld_orig[np.ix_(d_valid2, d_valid2)]\n        Yd = Yd_orig[d_valid2]\n\n        CLC_p = Cp.T @ (Lp @ Cp)\n        CLC_d = Cd.T @ (Ld @ Cd)\n\n        CLY_p = Cp.T @ (Lp @ Yp)\n        CLY_d = Cd.T @ (Ld @ Yd)\n\n        if const['use_sensor']:\n            Cv = const['Cv']\n            Lv = const['Lv']\n            Ca = const['Ca']\n            La = const['La']\n            Ya = const['Ya']\n            CLC_v = Cv.T @ (Lv @ Cv)\n            CLC_a = Ca.T @ (La @ Ca)\n            Q     = CLC_p + CLC_d + CLC_v + CLC_a\n\n            CLY_a = Ca.T @ (La @ Ya)\n            q     = CLY_p + CLY_d + CLY_a\n        else:\n            Q = CLC_p + CLC_d\n            q = CLY_p + CLY_d\n\n        if const['use_inquality']:\n            G = const['Gv']\n            h = const['hv']\n            X_star = qpsolver.solve_qp_with_inequality(R=R, Q=Q, q=q, A=A, B=B, G=G, h=h)\n        else:\n            X_star = qpsolver.solve_qp(R=R, Q=Q, q=q, A=A, B=B)\n        return X_star\n\n\n    def get_baseline(self, collection_name):\n        df = self.BASELINE_DF[self.BASELINE_DF['collectionName'] == collection_name].copy()\n        df.reset_index(drop=True, inplace=True)\n        return df\n\n    def get_velocity(self, collection_name):\n        df = self.VELOCITY_DF[self.VELOCITY_DF['collectionName'] == collection_name].copy()\n        df.reset_index(drop=True, inplace=True)\n        return df\n\n    def apply_costmin(self, base_df, velocity_df, sensor_df, params, N_LOOP):\n        const = self.get_optimization_constants(base_df, velocity_df, sensor_df, params, use_sensor=True)\n\n        if params['use_map']:\n            distance = map_matching.distance_to_nearest_neighbor(base_df)\n            default_p_valid = (distance < params['threshold_distance_to_nearest_neighbor'])\n            p_valid = default_p_valid\n        else:\n            default_p_valid = np.full(const['N_y'], True)\n            p_valid = default_p_valid\n\n        V = np.sqrt(np.sum(velocity_df[['v_east', 'v_north']].values**2, axis=1))\n        default_d_valid = (V < params['vmax'])\n        d_valid = default_d_valid\n\n        for loop in range(N_LOOP):\n            X_star = self.solve_QP(const, p_valid, d_valid)\n            Y_star = const['Cp_orig'] @ X_star\n            Y_star = np.reshape(Y_star, (-1, 2))\n            pp_df  = base_df.copy()\n            pp_df['latDeg'] = Y_star[:, 0]\n            pp_df['lngDeg'] = Y_star[:, 1]\n            distance = transform.pd_haversine_distance(pp_df, base_df)\n            p_valid = default_p_valid & (distance < params['reject_p'])\n\n            dXYdt = const['Cd_orig'] @ X_star\n            dXYdt = np.reshape(dXYdt, (-1, 2))\n            v_err = dXYdt - velocity_df[['v_east', 'v_north']].values\n            v_err = np.sqrt(np.sum(v_err**2, axis=1))\n            d_valid = default_d_valid & (v_err < params['reject_d'])\n\n        return pp_df\n    \n    \n    def recalibrate_sensor_by_vehicle_motion(self, base_df, velocity_df, sensor_df, params):\n        const = self.get_optimization_constants(base_df, velocity_df, sensor_df, params, use_sensor=False)\n\n        if params['use_map']:\n            distance = map_matching.distance_to_nearest_neighbor(base_df)\n            p_valid  = (distance < params['threshold_distance_to_nearest_neighbor'])\n        else:\n            p_valid = np.full(const['N_y'], True)\n\n        V = np.sqrt(np.sum(velocity_df[['v_east', 'v_north']].values**2, axis=1))\n        d_valid = (V < params['vmax'])\n\n        X_star = self.solve_QP(const, p_valid, d_valid)\n        X_mat  = np.reshape(X_star, (-1, 6))\n        dotB   = X_mat[const['x_index_sensor'], 1]\n        dotL   = X_mat[const['x_index_sensor'], 4]\n        dotXY  = const['J'] @ np.stack([dotB, dotL], axis=0) # shape = (2, N)\n        dotX   = dotXY[0, :]\n        dotY   = dotXY[1, :]\n        V = np.sqrt(np.sum(dotXY**2, axis=0))\n        cond = (V > (20 / 3.6))\n        trig_moving_direction = signal_f.Trig.from_data(dotX[cond], dotY[cond])\n        trig_theta   = signal_f.Trig.from_rad(sensor_df['theta'].values[cond])\n        trig_offset  = trig_theta - trig_moving_direction\n        angle_offset = np.arctan2(np.mean(trig_offset.sin), np.mean(trig_offset.cos))\n        sensor_df['theta']  = sensor_df['theta'] - angle_offset\n        sensor_df['cos_th'] = np.cos(sensor_df['theta'].values)\n        sensor_df['sin_th'] = np.sin(sensor_df['theta'].values)\n\n        return sensor_df\n\n    \n    def do_postprocess(self, args):\n        train_or_test, collection, params = args\n\n        base_df = self.get_baseline(collection)\n        t_ref   = base_df['millisSinceGpsEpoch'].min()\n        base_df['Time'] = 1e-3 * (base_df['millisSinceGpsEpoch'] - t_ref).values\n\n        velocity_df = self.get_velocity(collection)\n        velocity_df['Time'] = (1e-3 * (velocity_df['millisSinceGpsEpoch'] - t_ref).values\n                               -  params['Mi8_velocity_timeshift'] * (velocity_df['phoneName'] == 'Mi8').astype(float)\n                               )\n        velocity_df = velocity_df[(  velocity_df['Time'] >= base_df['Time'].min())\n                                  & (velocity_df['Time'] <= base_df['Time'].max())]\n        velocity_df.reset_index(drop=True, inplace=True)\n\n        phone_list = [path.split('/')[-1] for path in sorted(glob(f'{BASE_DIR}/{train_or_test}/{collection}/*'))]\n        sensor_df_list   = []\n        dt_up   = 2.5 * 1e-3\n        dt_down = self.DT_X\n        FLT = design_filter.make_sinc_filter(F_cutoff=2.0, dt=dt_up)\n        for phone in phone_list:\n            gnss_log_filename = f'{BASE_DIR}/{train_or_test}/{collection}/{phone}/{phone}_GnssLog.txt'\n            sensor_df_orig = io_f.read_GnssLog_sensors(gnss_log_filename)\n            if signal_f.check_sensor_availability(sensor_df_orig):\n                sensor_df = signal_f.preprocess_sensor_data(sensor_df_orig, t_ref, dt_up, dt_down, FLT)\n                sensor_df = signal_f.remove_different_posture(sensor_df)\n                sensor_df = sensor_df[(  sensor_df['Time'] >= base_df['Time'].min())\n                                      & (sensor_df['Time'] <= base_df['Time'].max())].copy()\n                sensor_df.reset_index(drop=True, inplace=True)\n                sensor_df_list.append(sensor_df)\n        if len(sensor_df_list) > 0:\n            time_list = [df['Time'].max() - df['Time'].min() for df in sensor_df_list]\n            idx = np.argmax(time_list)\n            sensor_df = sensor_df_list[idx]\n            sensor_df = signal_f.add_calibrated_signals(sensor_df, dt_down)\n            sensor_df = self.recalibrate_sensor_by_vehicle_motion(base_df, velocity_df, sensor_df, params)\n        else:\n            sensor_df = None\n\n        pp_df = base_df\n        pp_df = self.apply_costmin(pp_df, velocity_df, sensor_df, params, N_LOOP=3)\n        if params['use_map']:\n            params_stage2 = dict(params)\n            params_stage2['sigma_p'] = params['sigma_p_stage2']\n            for _ in range(params['num_stage2_iterations']):\n                pp_df = map_matching.snap_to_nearest_neighbor(pp_df)\n                pp_df = self.apply_costmin(pp_df, velocity_df, sensor_df, params_stage2, N_LOOP=1)\n        return pp_df\n    \n\n    def main(self, train, test):\n        VELOCITY_PATH = '../input/vehicle-speed-estimation/_doppler_velocity'\n        \n        # train\n        self.train_or_test = \"train\"\n        self.BASELINE_DF = train\n        self.VELOCITY_DF = pd.read_csv(f'{VELOCITY_PATH}/doppler_velocity_train.csv')\n        \n        collection_list_highway  = TRAIN_HIGHWAY\n        collection_list_treeway  = TRAIN_TREEWAY\n        collection_list_downtown = TRAIN_DOWNTOWN\n        config = [\n            (collection_list_highway,  self.params_highway),\n            (collection_list_treeway,  self.params_treeway),\n            (collection_list_downtown, self.params_downtown),\n        ]\n        args_list = []\n        for collection_list, params in config:\n            for collection in collection_list:\n                args_list.append( ('train', collection, params) )\n        train = pd.merge(\n            train.drop(columns=[\"latDeg\", \"lngDeg\"]),\n            self.single_main(args_list)[[\"phone\", \"millisSinceGpsEpoch\", \"latDeg\", \"lngDeg\"]],\n            on=[\"phone\", \"millisSinceGpsEpoch\"], how=\"inner\")\n        \n        # test      \n        self.train_or_test = \"test\"      \n        self.BASELINE_DF = test\n        self.VELOCITY_DF = pd.read_csv(f'{VELOCITY_PATH}/doppler_velocity_test.csv')\n        \n        collection_list_all = np.array(sorted(path.split('/')[-1] for path in glob(f'{BASE_DIR}/test/*')))\n        collection_list_highway  = TEST_HIGHWAY\n        collection_list_treeway  = TEST_TREEWAY\n        collection_list_downtown = TEST_DOWNTOWN\n        config = [\n            (collection_list_highway,  self.params_highway),\n            (collection_list_treeway,  self.params_treeway),\n            (collection_list_downtown, self.params_downtown),\n        ]\n        args_list = []\n        for collection_list, params in config:\n            for collection in collection_list:\n                args_list.append( ('test', collection, params) )\n                \n        test = pd.merge(\n            test.drop(columns=[\"latDeg\", \"lngDeg\"]),\n            self.single_main(args_list)[[\"phone\", \"millisSinceGpsEpoch\", \"latDeg\", \"lngDeg\"]],\n            on=[\"phone\", \"millisSinceGpsEpoch\"], how=\"inner\")\n        \n        return train, test\n        \n    def single_main(self, args_list):\n        processes = 8\n        with multiprocessing.Pool(processes=processes) as pool:\n            df_list = pool.imap_unordered(self.do_postprocess, args_list)\n            df_list = tqdm(df_list, total=len(args_list))\n            df_list = list(df_list)\n\n        columns = ['phone', 'millisSinceGpsEpoch', 'latDeg', 'lngDeg', 'collectionName', 'phoneName']\n        sub_df = pd.concat(df_list)\n        sub_df = sub_df[columns]\n        sub_df = sub_df.sort_values(['phone', 'millisSinceGpsEpoch'])\n\n        return sub_df","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Bias Correcetion","metadata":{}},{"cell_type":"code","source":"class BiasCorrection:\n    name = 'bias_correction'\n    isPrep = False\n    \n    def __init__(self, bias_x=0.0, bias_y=0.5):\n        self.bias_x = bias_x\n        self.bias_y = bias_y\n    \n    def main(self, train, test,bias_x=0.0,bias_y=0.5):\n        train = self._bias_correction(train)\n        test = self._bias_correction(test)\n        return train, test\n\n    def _bias_correction_phone(self, args):\n        phone, phone_df = args\n\n        B = np.deg2rad(phone_df['latDeg'].values)\n        L = np.deg2rad(phone_df['lngDeg'].values)\n        H = np.zeros_like(B)\n        BLH = transform.BLH(lat=B, lng=L, hgt=H)\n        J = transform.jacobian_BL_to_EN(BLH)\n\n        t_ref  = phone_df['millisSinceGpsEpoch'].min()\n        TIME   = 1e-3 * (phone_df['millisSinceGpsEpoch'] - t_ref).values\n        dotB   = InterpolatedUnivariateSpline(TIME, B, k=3).derivative()(TIME)\n        dotL   = InterpolatedUnivariateSpline(TIME, L, k=3).derivative()(TIME)\n        dotBL  = np.stack([dotB, dotL], axis=1)\n        dotEN  = np.einsum('nij,nj->ni', J, dotBL)\n        absV   = np.sqrt(np.sum(dotEN**2, axis=1))\n        th_az  = np.arctan2(dotEN[:, 0], dotEN[:, 1])\n\n        cos_az = np.cos(th_az)\n        sin_az = np.sin(th_az)\n        valid  = (absV > (5 / 3.6))\n        cos_az = InterpolatedUnivariateSpline(TIME[valid], cos_az[valid], k=1, ext=3)(TIME)\n        sin_az = InterpolatedUnivariateSpline(TIME[valid], sin_az[valid], k=1, ext=3)(TIME)\n        th_az  = np.arctan2(sin_az, cos_az)\n        cos_az = np.cos(th_az)\n        sin_az = np.sin(th_az)\n\n        delta_X  = - self.bias_x\n        delta_Y  = - self.bias_y\n        delta_E  = (  cos_az * delta_X) + (sin_az * delta_Y)\n        delta_N  = (- sin_az * delta_X) + (cos_az * delta_Y)\n        delta_EN = np.stack([delta_E, delta_N], axis=0) # shape = (2, N)\n        Jinv = np.linalg.inv(np.mean(J, axis=0))\n        delta_BL_rad = Jinv @ delta_EN\n        delta_BL_deg = np.rad2deg(delta_BL_rad)\n\n        output_df = pd.DataFrame({\n            'phone'               : phone_df['phone'],\n            'millisSinceGpsEpoch' : phone_df['millisSinceGpsEpoch'],\n            'latDeg'              : phone_df['latDeg'] + delta_BL_deg[0, :],\n            'lngDeg'              : phone_df['lngDeg'] + delta_BL_deg[1, :],\n        })\n        return output_df\n    \n    def _bias_correction(self, base_df):\n        output_df = base_df.sort_values(['phone', 'millisSinceGpsEpoch']).reset_index(drop=True).copy()\n        output_df_list = map(self._bias_correction_phone, base_df.groupby('phone'))\n        _df = pd.concat(output_df_list, axis=0)\n        _df = _df.sort_values(['phone', 'millisSinceGpsEpoch']).reset_index(drop=True)\n        output_df[['latDeg','lngDeg']] = _df[['latDeg','lngDeg']]\n        return output_df","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Pipeline","metadata":{}},{"cell_type":"code","source":"class Pipeline:\n        \n    def __init__(self, train, test, processes):\n        self.train = train\n        self.test = test\n        self.processes = processes\n        \n        self.pipe_name = \"\"     \n        self.process_name = \"\" \n        self.process_names = []\n        self.cvs = []    \n        \n    def main(self):\n        global PROCESSED_DFS\n        \n        print(f\"baseline...\")\n        self.train, cv = self.check_meter(self.train)\n        self.process_names.append(\"baseline\")\n        self.cvs.append(cv)        \n        print(\"\\n\")\n        \n        for process in self.processes:\n            if process.name == \"skip\":\n                continue\n            \n            print(f\"{process.name}...\")\n            \n            # reuse\n            if process.isPrep == False:                \n                self.process_names.append(process.name)\n                self.process_name += \"_\" + process.name\n            self.pipe_name += \"_\" + process.name\n            \n            if self.pipe_name in PROCESSED_DFS:\n                print(f\"loaded past result: {self.pipe_name}\")\n                self.train = PROCESSED_DFS[self.pipe_name][0].copy()\n                self.test  = PROCESSED_DFS[self.pipe_name][1].copy()\n            else:\n                self.train, self.test = process.main(self.train, self.test)\n                PROCESSED_DFS[self.pipe_name] = [self.train.copy(), self.test.copy()]\n            \n            if process.isPrep == False:\n                self.train, cv = self.check_meter(self.train)\n                self.cvs.append(cv)\n                \n            print(\"\\n\")\n            self.train.to_csv(f\"train_after{process.name}.csv\", index=False)\n            self.test.to_csv(f\"test_after{process.name}.csv\", index=False)\n            \n        self.train.to_csv(f\"train{self.process_name}_FIN.csv\", index=False)\n        self.test.to_csv(f\"test{self.process_name}_FIN.csv\", index=False)\n            \n        self.plot()\n        return self.train, self.test\n                        \n            \n    def check_meter(self, 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}')\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}') \n\n        return output_df, score\n    \n    def plot(self):\n        plt.subplots(figsize=(8, 3))\n        plt.plot(self.process_names, self.cvs, marker=\"o\")\n        \n        for l, c in zip(self.process_names, self.cvs):\n            plt.text(l, c+0.05, f\"{c:.2f}\")\n        plt.grid()\n        plt.ylabel(\"CV\")\n        plt.title(self.process_name[1:])\n        plt.show()\nPROCESSED_DFS = {}","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Input","metadata":{}},{"cell_type":"code","source":"BASE_DIR = Path('../input/google-smartphone-decimeter-challenge')\nSAITO_DIR = Path('../input/gsdc-improved-raw-gnss-baseline-result')\nOUTPUT_DIR = Path('./output/')\nos.makedirs(OUTPUT_DIR, exist_ok=True)\n\n# aggregate ground_truth(common)\ntrain = pd.read_csv(SAITO_DIR / 'raw_gnss_train.csv')\ntrain[\"collectionName\"] = train[\"phone\"].apply(lambda x:x.split(\"_\")[0])\ntrain[\"phoneName\"] = train[\"phone\"].apply(lambda x:x.split(\"_\")[1])\n\ntest = pd.read_csv(SAITO_DIR / 'raw_gnss_test.csv')\ntest[\"collectionName\"] = test[\"phone\"].apply(lambda x:x.split(\"_\")[0])\ntest[\"phoneName\"] = test[\"phone\"].apply(lambda x:x.split(\"_\")[1])\n\ntrain_gt = pd.concat([pd.read_csv(path) for path in tqdm(BASE_DIR.glob('*/*/*/ground_truth.csv'),total=73)])\ntrain_gt.rename(columns={'latDeg':'t_latDeg',\n                         'lngDeg':'t_lngDeg',\n                         'heightAboveWgs84EllipsoidM':'t_heightAboveWgs84EllipsoidM'}\n               ,inplace=True)\ntrain = train.merge(\n    train_gt, on=['collectionName', 'phoneName', 'millisSinceGpsEpoch']\n)\n\nrel_pos = pd.read_csv(\"../input/chris-baseline/gps-delta-lat-lng.csv\").rename(columns={\"latDeg\": \"delta_latDeg\", \"lngDeg\":\"delta_lngDeg\"})\ntrain = train.merge(rel_pos, on=[\"phone\", \"millisSinceGpsEpoch\"])\ntest = test.merge(rel_pos, on=[\"phone\", \"millisSinceGpsEpoch\"])\n\ntrain = check_meter(train)\ndel train_gt","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Apply","metadata":{}},{"cell_type":"code","source":"pipe = Pipeline(\n    train=train, test=test,\n\n    processes=[\n#         KnnHeight(),        \n        \n        OutlierDetection(\n            sjc_loss_threshold=9.5, \n            other_loss_threshold=5.6,\n        ),     \n                \n        InterpolateByRelpos(\n            sjc_start=0.1,\n            sjc_end=0.5,\n            other_start=0.1,\n            other_end=0.5,\n        ),\n        \n        StopMean(),\n\n        SaitoOptimization(\n            N_SNAP_TO_GRID_ITERATION=1\n        ),\n        \n        BiasCorrection(\n            bias_x=0.30,\n            bias_y=0.5\n        ),\n    ])\n\ntrain, test = pipe.main()\ntest[[\"phone\", \"millisSinceGpsEpoch\", \"latDeg\", \"lngDeg\"]].to_csv(\"submission.csv\", index=False)","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# VISUALIZE","metadata":{}},{"cell_type":"code","source":"visualize_trafic(test, color='phone', savepath=str(OUTPUT_DIR /'submission.html'))","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def evaluate_train(train):\n    dfs = []\n    for c, c_df in train.groupby(\"collectionName\"):\n        data = {}\n        for phone, phone_df in c_df.groupby(\"phoneName\"):\n            p_50 = np.percentile(phone_df.meter, 50)\n            p_95 = np.percentile(phone_df.meter, 95)\n            data[phone] = [(p_50+p_95)/2]\n\n        df = pd.DataFrame(data, index=[c])\n        dfs.append(df)        \n\n    score_df = pd.concat(dfs).fillna(\"-\")\n    score_df.index=train.collectionName.unique()\n    return score_df.copy()\nscore_df = evaluate_train(train)\nscore_df.to_csv(OUTPUT_DIR / \"train_CV.csv\")\nscore_df","metadata":{"trusted":true},"execution_count":null,"outputs":[]}]}