{"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":"# Correcting baseline with lightgbm ","metadata":{}},{"cell_type":"markdown","source":"## Additive functions","metadata":{}},{"cell_type":"code","source":"import numpy as np \nimport pandas as pd \nfrom sklearn.preprocessing import MinMaxScaler\nimport matplotlib.pyplot as plt\nfrom tqdm.notebook import tqdm\nimport pathlib\nfrom sklearn.model_selection import train_test_split\nfrom sklearn.ensemble import RandomForestRegressor as RF\nfrom sklearn.model_selection import train_test_split\nfrom sklearn.model_selection import KFold\nfrom xgboost import XGBRegressor\nimport itertools\nfrom catboost import CatBoostRegressor\nimport lightgbm as lgb\nfrom sklearn.ensemble import StackingRegressor\nfrom sklearn.linear_model import LinearRegression\nfrom catboost import CatBoostRegressor\nfrom lightgbm import LGBMRegressor\nfrom xgboost import XGBRegressor\nfrom sklearn.ensemble import RandomForestRegressor\nfrom sklearn.model_selection import GroupKFold\nimport optuna\n\nimport plotly.express as px\nimport torch\nfrom torch.utils.data import Dataset\nfrom torchvision import datasets\nfrom torchvision.transforms import ToTensor\nfrom torch.utils.data import DataLoader\nfrom sklearn.preprocessing import MinMaxScaler as MMScaler, StandardScaler as SSScaler\nfrom sklearn.preprocessing import OneHotEncoder\n\nimport os\nfrom torchvision.io import read_image\n\n# Input data files are available in the read-only \"../input/\" directory\n# For example, running this (by clicking run or pressing Shift+Enter) will list all files under the input directory\n\n# import os\n# for dirname, _, filenames in os.walk('/kaggle/input'):\n#     for filename in filenames:\n#         print(os.path.join(dirname, filename))\n\n#print(os.listdir())\n\n# You can write up to 20GB to the current directory (/kaggle/working/) that gets preserved as output when you create a version using \"Save & Run All\" \n# You can also write temporary files to /kaggle/temp/, but they won't be saved outside of the current session","metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","execution":{"iopub.status.busy":"2021-08-04T21:17:23.296278Z","iopub.execute_input":"2021-08-04T21:17:23.296843Z","iopub.status.idle":"2021-08-04T21:17:23.715262Z","shell.execute_reply.started":"2021-08-04T21:17:23.296794Z","shell.execute_reply":"2021-08-04T21:17:23.714507Z"},"jupyter":{"source_hidden":true},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def 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    RADIUS = 6_367_000\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)**2 + \\\n        np.cos(lat1) * np.cos(lat2) * np.sin(dlon/2)**2\n    dist = 2 * RADIUS * np.arcsin(a**0.5)\n    return dist\n\ndef compute_dist_hav(df, cols):\n    df = df.copy()\n    #cols = ['latDeg', 'lngDeg']\n    prev_cols = [col+'_prev' for col in cols]\n    \n    df.loc[:, ['dist']] = np.nan\n    \n    for phone in tqdm(df['phone'].unique()):\n        ser = df[df['phone'] == phone]\n        \n        ser.loc[:, prev_cols] = ser.loc[:, cols].shift(1).values\n        #ser['latDeg_prev'] = ser['latDeg'].shift(1)\n        #ser['lngDeg_prev'] = ser['lngDeg'].shift(1)\n        \n        #display(ser)\n        f_ind = ser.index[0]\n        ser.loc[f_ind, prev_cols] = ser.loc[f_ind, cols].values\n        #ser.loc[f_ind, 'latDeg_prev'] = ser.loc[f_ind, 'latDeg']\n        #ser.loc[f_ind, 'lngDeg_prev'] = ser.loc[f_ind, 'lngDeg']\n        \n        df.loc[(df['phone'] == phone), prev_cols] = ser[prev_cols]\n        #df.loc[(df['phone'] == phone), ['latDeg_prev', 'lngDeg_prev']] = ser[['latDeg_prev', 'lngDeg_prev']]\n        \n    df['dist'] = calc_haversine(df[cols[0]], df[cols[1]], df[prev_cols[0]], df[prev_cols[1]])\n    return df\n\ndef get_deltas(t):\n    t = t.copy()\n    \n    t.loc[:, ['d_latDeg', 'd_lngDeg']] = np.nan\n    #display(t)\n    \n    for phone in tqdm(t['phone'].unique()):\n        ser = t[t['phone'] == phone]\n        \n        ser['d_latDeg'] = ser['latDeg'] - ser['latDeg'].shift(1)\n        ser['d_lngDeg'] = ser['lngDeg'] - ser['lngDeg'].shift(1)\n        ser = ser.fillna(0)\n        \n        #display(t.loc[t['phone'] == phone])\n        t.loc[(t['phone'] == phone), ['d_latDeg', 'd_lngDeg']] = ser[['d_latDeg', 'd_lngDeg']]\n    \n    return t\n\ndef get_moving_average(t, w=5):\n    t = t.copy()\n    t2 = pd.DataFrame([])\n    for phone in t['phone'].unique():\n        t3 = t[t['phone'] == phone].rolling(w).mean().shift(-w//2)\n        t3 = t3.interpolate(method='linear', limit_direction='both', axis=0)\n\n        if t2.shape[0] == 0:\n            t2 = t3\n        else:\n            t2 = pd.concat([t2, t3], axis=0)\n\n    t = t.loc[:, [col for col in t if col not in t2.columns]]\n    t = pd.concat([t, t2], axis=1)\n    \n    return t\n\ndef weighted_average(df, data_col, weight_col, by_col):\n    df['_data_times_weight'] = df[data_col] * df[weight_col]\n    df['_weight_where_notnull'] = df[weight_col] * pd.notnull(df[data_col])\n    g = df.groupby(by_col)\n    result = g['_data_times_weight'].sum() / g['_weight_where_notnull'].sum()\n    del df['_data_times_weight'], df['_weight_where_notnull']\n    return result","metadata":{"execution":{"iopub.status.busy":"2021-08-04T21:17:23.809624Z","iopub.execute_input":"2021-08-04T21:17:23.810227Z","iopub.status.idle":"2021-08-04T21:17:23.830277Z","shell.execute_reply.started":"2021-08-04T21:17:23.810178Z","shell.execute_reply":"2021-08-04T21:17:23.829473Z"},"jupyter":{"source_hidden":true},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def get_score(df_to_score, ground_truth):\n    df_to_score = df_to_score.copy()\n    \n    df_to_score.sort_values(by=['phone', 'time'], inplace=True)\n    ground_truth.sort_values(by=['phone', 'time'], inplace=True)\n\n    df_to_score['t_latDeg'] = ground_truth['latDeg'].values\n    df_to_score['t_lngDeg'] = ground_truth['lngDeg'].values\n    \n    meter_score, score = _check_score(df_to_score)\n    return meter_score, score\n\ndef _check_score(input_df, silent=True):\n    output_df = input_df.copy()\n    \n    output_df['meter'] = calc_haversine(\n            input_df.latDeg, input_df.lngDeg, input_df.t_latDeg, input_df.t_lngDeg\n        )\n\n    meter_score = output_df['meter'].mean()\n    \n    if not silent:\n        print(f'error meter: {meter_score}')\n\n    scores_50 = []\n    scores_95 = []\n    for phone in output_df['phone'].unique():\n        _index = output_df['phone']==phone\n        p_50 = np.percentile(output_df.loc[_index, 'meter'], 50)\n        p_95 = np.percentile(output_df.loc[_index, 'meter'], 95)\n        scores_50.append(p_50)\n        scores_95.append(p_95)\n    \n    scores = scores_50 + scores_95\n    \n    score_50 = sum(scores_50) / len(scores_50)\n    score_95 = sum(scores_95) / len(scores_95)\n    score = sum(scores) / len(scores)\n    \n    if not silent:\n        print(f'score 50: {score_50}')\n        print(f'score 95: {score_95}')\n        print(f'score   : {score}')\n    \n    return meter_score, score","metadata":{"execution":{"iopub.status.busy":"2021-08-04T21:17:23.998599Z","iopub.execute_input":"2021-08-04T21:17:23.999149Z","iopub.status.idle":"2021-08-04T21:17:24.009963Z","shell.execute_reply.started":"2021-08-04T21:17:23.999113Z","shell.execute_reply":"2021-08-04T21:17:24.008825Z"},"jupyter":{"source_hidden":true},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# В некоторых местах преобразование не точно на 1 секунду. Но в рамках этой истории - это не важно\ndef utc_to_gps(time):#+ timedelta(seconds=1092121243.0 - (35 - 19))\n    return time - 315964800000 + 18000\n\ndef gps_to_utc(time):\n    return time + 315964800000 - 18000","metadata":{"execution":{"iopub.status.busy":"2021-08-04T21:17:24.307141Z","iopub.execute_input":"2021-08-04T21:17:24.307513Z","iopub.status.idle":"2021-08-04T21:17:24.311899Z","shell.execute_reply.started":"2021-08-04T21:17:24.307483Z","shell.execute_reply":"2021-08-04T21:17:24.31119Z"},"jupyter":{"source_hidden":true},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def visualize_trafic(df, center, zoom=15):\n    fig = px.scatter_mapbox(df,\n                            \n                            # Here, plotly gets, (x,y) coordinates\n                            lat=\"latDeg\",\n                            lon=\"lngDeg\",\n                            \n                            #Here, plotly detects color of series\n                            color=\"phoneName\",\n                            labels=\"phoneName\",\n                            \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    \ndef visualize_collection(df):\n    target_df = df\n    lat_center = target_df['latDeg'].mean()\n    lng_center = target_df['lngDeg'].mean()\n    center = {\"lat\":lat_center, \"lon\":lng_center}\n    \n    visualize_trafic(target_df, center)","metadata":{"execution":{"iopub.status.busy":"2021-08-04T21:17:24.475771Z","iopub.execute_input":"2021-08-04T21:17:24.476262Z","iopub.status.idle":"2021-08-04T21:17:24.483746Z","shell.execute_reply.started":"2021-08-04T21:17:24.476227Z","shell.execute_reply":"2021-08-04T21:17:24.482726Z"},"jupyter":{"source_hidden":true},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def get_lag_features(df_train, columns, count=15):\n    \"\"\" Get lag features ONLY for separated collections+phones\"\"\"\n    new_feat_df = dict()\n    \n    t_size = df_train.shape[0]\n    for col in columns:\n        t = df_train[col]\n        for i in range(1, count + 1):\n            new_feat_df[f'd_next_{col}_{i}'] = t - t.shift(i).fillna(method='bfill')\n        for i in range(1, count + 1):\n            new_feat_df[f'd_prev_{col}_{i}'] = t - t.shift(-i).fillna(method='ffill')\n            \n    new_feat_df = pd.DataFrame(new_feat_df)\n    \n    return pd.concat( [df_train.reset_index(drop=True), new_feat_df.reset_index(drop=True)], axis=1)\n\n\ndef get_collections_list(df):\n    \"\"\" Separate dataframe on collections\"\"\"\n    df_list = []\n    for phone in df['phone'].unique():\n        df_list.append( df[df['phone'] == phone] )\n        \n    return df_list","metadata":{"execution":{"iopub.status.busy":"2021-08-04T21:17:24.645889Z","iopub.execute_input":"2021-08-04T21:17:24.646226Z","iopub.status.idle":"2021-08-04T21:17:24.655808Z","shell.execute_reply.started":"2021-08-04T21:17:24.646197Z","shell.execute_reply":"2021-08-04T21:17:24.65449Z"},"jupyter":{"source_hidden":true},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Prepare input","metadata":{}},{"cell_type":"code","source":"TRAIN_INPUT_DIR = '/kaggle/input/gps-datasets/' \ntrain_fname = 'train_clean.csv'\n#TRAIN_INPUT_DIR = '/kaggle/input/data-derived/' \n#train_fname = 'df_train_derived.csv'\ntrain_base = pd.read_csv(TRAIN_INPUT_DIR + train_fname)\ntrain_base.head(2)","metadata":{"execution":{"iopub.status.busy":"2021-08-04T21:17:27.052751Z","iopub.execute_input":"2021-08-04T21:17:27.053141Z","iopub.status.idle":"2021-08-04T21:17:27.48039Z","shell.execute_reply.started":"2021-08-04T21:17:27.053108Z","shell.execute_reply":"2021-08-04T21:17:27.479271Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"train_base.drop(['d_latDeg', 'd_lngDeg', 'dist',\n                 'latDeg_prev', 'lngDeg_prev'], axis=1, inplace=True)","metadata":{"execution":{"iopub.status.busy":"2021-08-04T21:17:30.809291Z","iopub.execute_input":"2021-08-04T21:17:30.809686Z","iopub.status.idle":"2021-08-04T21:17:30.823332Z","shell.execute_reply.started":"2021-08-04T21:17:30.809651Z","shell.execute_reply":"2021-08-04T21:17:30.822095Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"p = pathlib.Path('../input/google-smartphone-decimeter-challenge')\ngt_files = list(p.glob('train/*/*/ground_truth.csv'))\n\ngts = []\nfor gt_file in gt_files:\n    gts.append(pd.read_csv(gt_file))\n    \nground_truth = pd.concat(gts)\nground_truth['phone'] = ground_truth['collectionName'].astype(str) + '_' + ground_truth['phoneName']\nground_truth['time'] = pd.to_datetime(gps_to_utc(ground_truth['millisSinceGpsEpoch'])//1000, unit='s')\nground_truth = compute_dist_hav(ground_truth, ['latDeg', 'lngDeg'])","metadata":{"execution":{"iopub.status.busy":"2021-08-04T21:17:34.12208Z","iopub.execute_input":"2021-08-04T21:17:34.122584Z","iopub.status.idle":"2021-08-04T21:17:38.584776Z","shell.execute_reply.started":"2021-08-04T21:17:34.122552Z","shell.execute_reply":"2021-08-04T21:17:38.583698Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"train_base.sort_values(by=['phone', 'time'], inplace=True)\nground_truth.sort_values(by=['phone', 'time'], inplace=True)\n\ntrain_base.reset_index(drop=True, inplace=True)\nground_truth.reset_index(drop=True, inplace=True)","metadata":{"execution":{"iopub.status.busy":"2021-08-04T21:17:38.586599Z","iopub.execute_input":"2021-08-04T21:17:38.586996Z","iopub.status.idle":"2021-08-04T21:17:38.784835Z","shell.execute_reply.started":"2021-08-04T21:17:38.586955Z","shell.execute_reply":"2021-08-04T21:17:38.783768Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"train_base_list = get_collections_list(train_base)\nlag_columns = ['latDeg', 'lngDeg']\nlag_count = 10\nlag_dfs = [get_lag_features(df, lag_columns, lag_count) for df in train_base_list]\ntrain_base = pd.concat(lag_dfs, axis=0)\ntrain_base.reset_index(drop=True, inplace=True)","metadata":{"execution":{"iopub.status.busy":"2021-08-04T21:17:39.305336Z","iopub.execute_input":"2021-08-04T21:17:39.305715Z","iopub.status.idle":"2021-08-04T21:17:42.393766Z","shell.execute_reply.started":"2021-08-04T21:17:39.305684Z","shell.execute_reply":"2021-08-04T21:17:42.39264Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"GT = pd.DataFrame({\n    'target_d_lat' : (train_base['latDeg'].to_numpy() - ground_truth['latDeg'].to_numpy()),\n    'target_d_lng' : (train_base['lngDeg'].to_numpy() - ground_truth['lngDeg'].to_numpy())\n})","metadata":{"execution":{"iopub.status.busy":"2021-08-04T21:17:44.222707Z","iopub.execute_input":"2021-08-04T21:17:44.223128Z","iopub.status.idle":"2021-08-04T21:17:44.230289Z","shell.execute_reply.started":"2021-08-04T21:17:44.223088Z","shell.execute_reply":"2021-08-04T21:17:44.22956Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"train_cols = ['latDeg', 'lngDeg'] + list(filter(lambda w : w.startswith('d_'), train_base.columns))","metadata":{"execution":{"iopub.status.busy":"2021-08-04T21:17:44.534647Z","iopub.execute_input":"2021-08-04T21:17:44.535002Z","iopub.status.idle":"2021-08-04T21:17:44.540021Z","shell.execute_reply.started":"2021-08-04T21:17:44.534971Z","shell.execute_reply":"2021-08-04T21:17:44.538912Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"enc = OneHotEncoder()\nenc.fit(train_base.phoneName.unique().reshape(-1, 1))\nphone_feat = enc.transform(train_base.phoneName.values.reshape(-1, 1))\nphone_cols = [f'oh{i}' for i in range(phone_feat.shape[1])]\nphone_feat_df = pd.DataFrame( phone_feat.toarray(), columns=phone_cols )\ntrain_base = pd.concat([train_base, phone_feat_df], axis=1)\ntrain_cols += phone_cols","metadata":{"execution":{"iopub.status.busy":"2021-08-04T21:17:46.448917Z","iopub.execute_input":"2021-08-04T21:17:46.449283Z","iopub.status.idle":"2021-08-04T21:17:46.59941Z","shell.execute_reply.started":"2021-08-04T21:17:46.449251Z","shell.execute_reply":"2021-08-04T21:17:46.598271Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Optuna parameters search","metadata":{}},{"cell_type":"code","source":"def train_evaluate(param, trial):\n    \n    kfold = GroupKFold(n_splits=min(5, max(2, len(np.unique(groups))) ))\n    scores = []\n    \n    for fold_id, (trn_idx, val_idx) in enumerate(kfold.split(X, y, groups=groups)):\n        all_y_pred = []\n        \n        for ind in range(2):\n            \n            X_train = X[trn_idx]\n            X_val = X[val_idx]\n            y_train = y[trn_idx, ind]\n            y_val = y[val_idx, ind] \n            \n            dtrain = lgb.Dataset(X_train, y_train)\n            \n            dtest = lgb.Dataset(X_val, y_val)\n            gbm = lgb.train(param, dtrain, num_boost_round=1000, valid_sets=dtest, verbose_eval=500)\n            \n            y_pred = gbm.predict(X_val)\n            \n            all_y_pred.append( y_pred )\n\n        gt_val = ground_truth.iloc[val_idx]\n        ans_base = train_base.iloc[val_idx].copy()\n        \n        all_y_pred = np.array(all_y_pred).T\n        ans_base[['latDeg', 'lngDeg']] -= all_y_pred\n        meter_score, score = get_score(ans_base, gt_val)\n        \n        scores.append(score)\n    \n    return np.mean(scores)\n\n\ndef objective(trial):\n    \n    param = {'num_leaves': trial.suggest_int('num_leaves', 24, 1024),\n            'boosting_type': 'gbdt',\n            'objective': 'regression',\n            'metric': 'mae',\n            #'eval_metric': 'mae', \n            'feature_fraction': trial.suggest_uniform('feature_fraction', 0.4, 1.0),\n            'bagging_fraction': trial.suggest_uniform('bagging_fraction', 0.4, 1.0),\n            'bagging_freq': trial.suggest_int('bagging_freq', 1, 7),\n            'min_child_samples': trial.suggest_int('min_child_samples', 5, 100),\n            'lambda_l1': trial.suggest_loguniform('lambda_l1', 1e-8, 10.0),\n            'lambda_l2': trial.suggest_loguniform('lambda_l2', 1e-8, 10.0),\n            'early_stopping_rounds': 200,\n              #'device' : 'gpu',\n             'verbosity' : -1,\n             #'verbose_eval': -1\n            }\n    \n    return train_evaluate(param, trial)","metadata":{"execution":{"iopub.status.busy":"2021-08-04T21:20:51.464708Z","iopub.execute_input":"2021-08-04T21:20:51.465067Z","iopub.status.idle":"2021-08-04T21:20:51.479287Z","shell.execute_reply.started":"2021-08-04T21:20:51.465036Z","shell.execute_reply":"2021-08-04T21:20:51.478247Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def find_params():\n    X = train_base[train_cols].values\n    y = GT.values\n    groups = train_base.collectionName.values\n    \n    study = optuna.create_study(direction='minimize')\n    study.optimize(objective, n_trials=10)\n    \n    print('Number of finished trials: {}'.format(len(study.trials)))\n    print('Best trial:')\n    trial = study.best_trial\n\n    print('  Value: {}'.format(trial.value))\n    print('  Params: ')\n    for key, value in trial.params.items():\n        print('    {}: {}'.format(key, value))","metadata":{"execution":{"iopub.status.busy":"2021-08-04T21:20:51.929808Z","iopub.execute_input":"2021-08-04T21:20:51.930171Z","iopub.status.idle":"2021-08-04T21:20:51.948374Z","shell.execute_reply.started":"2021-08-04T21:20:51.930138Z","shell.execute_reply":"2021-08-04T21:20:51.947239Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# uncomment to find best params\n# find_params()","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Cross validation","metadata":{}},{"cell_type":"code","source":"params = {\n    \n    'num_leaves': 99,\n    'feature_fraction': 0.7544601013793489,\n    'bagging_fraction': 0.9508442472495611,\n    'bagging_freq': 5,\n    'min_child_samples': 51,\n    'lambda_l1': 0.00019899635179289397,\n    'lambda_l2': 0.2466948417767759,\n    \n    'early_stopping_rounds': 200,\n    'boosting_type': 'gbdt',\n    'objective': 'regression',\n    'metric': 'mae',\n    'verbosity' : -1,\n}","metadata":{"execution":{"iopub.status.busy":"2021-08-04T22:56:59.653279Z","iopub.execute_input":"2021-08-04T22:56:59.653672Z","iopub.status.idle":"2021-08-04T22:56:59.659545Z","shell.execute_reply.started":"2021-08-04T22:56:59.653639Z","shell.execute_reply":"2021-08-04T22:56:59.658862Z"},"jupyter":{"source_hidden":true},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def coord_cv_oof(params, groups):\n    n_splits = min(5, max(2, len(np.unique(groups))))\n    kfold = GroupKFold(n_splits=n_splits)\n    \n    scores = []\n    oofs = train_base.copy()\n    \n    models_lat = []\n    models_lng = []\n    \n    for fold_id, (trn_idx, val_idx) in enumerate(kfold.split(X, y, groups=groups)):\n        print(fold_id, '/', n_splits)\n        \n        X_train = X[trn_idx]\n        X_val = X[val_idx]\n        all_y_pred = []\n        for ind in range(2):\n            \n            y_train = y[trn_idx, ind]\n            y_val = y[val_idx, ind]\n            \n            dtrain = lgb.Dataset(X_train, y_train)\n            \n            dtest = lgb.Dataset(X_val, y_val)\n            model = lgb.train(params[ind], dtrain, num_boost_round=1000, valid_sets=dtest, verbose_eval=500)\n            y_pred = model.predict(X_val)\n            \n            all_y_pred.append( y_pred )\n            \n            if ind == 0:\n                models_lat.append(model)\n            else:\n                models_lng.append(model)\n        \n        gt_val = ground_truth.iloc[val_idx]\n        ans_base = train_base.iloc[val_idx].copy()\n        all_y_pred = np.array(all_y_pred).T\n        ans_base[['latDeg', 'lngDeg']] -= all_y_pred\n        #ans_base[['latDeg_bl', 'lngDeg_bl']] -= all_y_pred\n        meter_score, score = get_score(ans_base, gt_val)\n        scores.append( score )\n        \n        # There was Danil's bug with iloc\n        oofs.loc[val_idx, ['d_latDeg', 'd_lngDeg']] = all_y_pred\n\n        print(f'Fold {fold_id}: {scores[-1]}')\n        \n    score = np.mean(scores)\n    \n    print(params)\n    print('Validation:', score)\n    print('-' * 60)\n    \n    return score, scores, oofs, models_lat, models_lng","metadata":{"execution":{"iopub.status.busy":"2021-08-04T22:57:01.645241Z","iopub.execute_input":"2021-08-04T22:57:01.64562Z","iopub.status.idle":"2021-08-04T22:57:01.658723Z","shell.execute_reply.started":"2021-08-04T22:57:01.645589Z","shell.execute_reply":"2021-08-04T22:57:01.65773Z"},"jupyter":{"source_hidden":true},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def coord_cv_oof_zone(params, train_base, ground_truth, X, y, groups):\n    n_splits = min(5, max(2, len(np.unique(groups))))\n    kfold = GroupKFold(n_splits=n_splits)\n    \n    scores = []\n    oofs = train_base.copy()\n    \n    models_lat = []\n    models_lng = []\n    \n    for fold_id, (trn_idx, val_idx) in enumerate(kfold.split(X, y, groups=groups)):\n        print(fold_id, '/', n_splits)\n        \n        X_train = X[trn_idx]\n        X_val = X[val_idx]\n        all_y_pred = []\n        for ind in range(2):\n            \n            y_train = y[trn_idx, ind]\n            y_val = y[val_idx, ind]\n            \n            model = lgb.LGBMRegressor(**params)\n            model = model.fit(X_train, \n                            y_train,\n                            eval_metric=params['metric'])\n\n            y_pred = model.predict(X_val, num_iteration = model.best_iteration_)\n            all_y_pred.append( y_pred )\n            \n            if ind == 0:\n                models_lat.append(model)\n            else:\n                models_lng.append(model)\n        \n        gt_val = ground_truth.iloc[val_idx]\n        ans_base = train_base.iloc[val_idx].copy()\n        all_y_pred = np.array(all_y_pred).T\n        ans_base[['latDeg', 'lngDeg']] -= all_y_pred\n        meter_score, score = get_score(ans_base, gt_val)\n        scores.append( score )\n        \n        oofs.iloc[val_idx].loc[:, ['d_latDeg', 'd_lngDeg']] = all_y_pred\n\n        print(f'Fold {fold_id}: {scores[-1]}')\n        \n    score = np.mean(scores)\n    \n    print(params)\n    print('Validation:', score)\n    print('-' * 60)\n    \n    return scores, oofs, models_lat, models_lng","metadata":{"execution":{"iopub.status.busy":"2021-07-29T20:49:38.623143Z","iopub.execute_input":"2021-07-29T20:49:38.623476Z","iopub.status.idle":"2021-07-29T20:49:38.634487Z","shell.execute_reply.started":"2021-07-29T20:49:38.623445Z","shell.execute_reply":"2021-07-29T20:49:38.633354Z"},"jupyter":{"source_hidden":true},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def simple_crossval_scoring():\n    results = []\n    for road_type_num in range(3):\n        road_ind = train_base['road'] == road_names[road_type_num]\n        X = train_base[road_ind][train_cols].values\n        y = GT[road_ind].values\n        groups = train_base[road_ind].collectionName.values\n        print('type:', road_names[road_type_num])\n        results.append(\n            coord_cv_oof_zone(params_zone[road_type_num], train_base[road_ind], ground_truth[road_ind], X, y, groups)\n        )\n\n    return np.mean( np.mean([res[0] for res in results], axis=0) ), results\n    # result - 3.478843607529768\n    \n# # Use simple cross validation for predicting oof\n# score, results = simple_crossval_scoring()\n# oof_full = pd.concat([res[1] for res in results], axis=0).reset_index(drop=True)\n# print(f'val score: {score}, train score: {get_score(oof_full, ground_truth)}')","metadata":{"execution":{"iopub.status.busy":"2021-07-29T20:49:39.106654Z","iopub.execute_input":"2021-07-29T20:49:39.109554Z","iopub.status.idle":"2021-07-29T20:49:39.121307Z","shell.execute_reply.started":"2021-07-29T20:49:39.109386Z","shell.execute_reply":"2021-07-29T20:49:39.120117Z"},"jupyter":{"source_hidden":true},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def validate_one_zone(params, road_type_num, train_base_input, ground_truth_input):\n    \n    road_ind = train_base_input['road'] == road_names[road_type_num]\n    train_base = train_base_input[road_ind]\n    ground_truth = ground_truth_input[road_ind]\n    \n    X = train_base[train_cols].values\n    y = GT[road_ind].values\n    groups = train_base.collectionName.values\n    param = params[road_type_num]\n    \n    kfold = GroupKFold(n_splits=min(5, max(2, len(np.unique(groups))) ))\n    scores = []\n    \n    for fold_id, (trn_idx, val_idx) in enumerate(kfold.split(X, y, groups=groups)):\n        print(f'Fold {fold_id} ')\n        all_y_pred = []\n        \n        for ind in range(2):\n            \n            X_train = X[trn_idx]\n            X_val = X[val_idx]\n            y_train = y[trn_idx, ind]\n            y_val = y[val_idx, ind] \n            \n            dtrain = lgb.Dataset(X_train, y_train)\n            \n            #dtest = lgb.Dataset(X_val, y_val)\n            #gbm = lgb.train(param, dtrain, num_boost_round=1000, valid_sets=dtest, verbose_eval=500)\n            \n            gbm = lgb.train(param, dtrain, num_boost_round=1000)\n            y_pred = gbm.predict(X_val)\n            \n            all_y_pred.append( y_pred )\n            \n\n        gt_val = ground_truth.iloc[val_idx]\n        ans_base = train_base.iloc[val_idx].copy()\n        all_y_pred = np.array(all_y_pred).T\n        ans_base[['latDeg', 'lngDeg']] -= all_y_pred\n        meter_score, score = get_score(ans_base, gt_val)\n        \n        scores.append(score)\n        \n        print(scores[-1])\n    \n    return np.mean(scores)","metadata":{"execution":{"iopub.status.busy":"2021-07-29T21:18:42.551817Z","iopub.execute_input":"2021-07-29T21:18:42.552125Z","iopub.status.idle":"2021-07-29T21:18:42.562319Z","shell.execute_reply.started":"2021-07-29T21:18:42.552096Z","shell.execute_reply":"2021-07-29T21:18:42.561255Z"},"jupyter":{"source_hidden":true},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def coord_cv_oof_zone_real(params, train_base, ground_truth):\n    \n    results = []\n    n_splits = 3 #min(3, max(2, len(np.unique(groups))))\n    n_roads = 1\n    val_idx_list = [] # check correctness\n    \n    train_base.reset_index(drop=True, inplace=True)\n    ground_truth.reset_index(drop=True, inplace=True)\n    \n    oofs = train_base.copy()\n    \n    for road_type_num in range(n_roads):\n        road_ind = train_base['road'] == road_names[road_type_num]\n        X = train_base[road_ind][train_cols].values\n        y = GT[road_ind].values\n        groups = train_base[road_ind].collectionName.values\n    \n        kfold = GroupKFold(n_splits=n_splits)\n\n        scores = []\n        results_split = [None] * n_splits\n        \n        for fold_id, (trn_idx, val_idx) in enumerate(kfold.split(X, y, groups=groups)):\n            print(fold_id, '/', n_splits)\n\n            X_train = X[trn_idx]\n            X_val = X[val_idx]\n            all_y_pred = []\n            models = []\n            \n            for ind in range(2):\n\n                y_train = y[trn_idx, ind]\n                y_val = y[val_idx, ind]\n                dtrain = lgb.Dataset(X_train, y_train)\n                dtest = lgb.Dataset(X_val, y_val)\n\n                model = lgb.train(params[road_type_num], dtrain, num_boost_round=1000, valid_sets=dtest, verbose_eval=500)\n                #model = lgb.train(params[road_type_num], dtrain, num_boost_round=1000)\n                y_pred = model.predict(X_val)\n\n                all_y_pred.append( y_pred )\n                models.append(model)\n            \n            global_ind_val = train_base.index[road_ind][val_idx].values \n            # may be global_ind_val = train_base.index[road_ind].values[val_idx], didn't check\n            \n            gt_val = ground_truth.iloc[global_ind_val].copy()\n            ans_base = train_base.iloc[global_ind_val].copy()\n            all_y_pred = np.array(all_y_pred).T\n            ans_base[['latDeg', 'lngDeg']] -= all_y_pred\n            \n            oofs.loc[global_ind_val, ['d_latDeg', 'd_lngDeg']] = all_y_pred\n            val_idx_list += global_ind_val.tolist()\n            \n            results_split[fold_id] = (ans_base, gt_val, models[0], models[1])\n        \n        results.append(results_split)\n        \n    scores = []\n    for i in range(n_splits):\n        ans_base = pd.concat([results[road_type][i][0] for road_type in range(n_roads)], axis=0)\n        gt_val = pd.concat([results[road_type][i][1] for road_type in range(n_roads)], axis=0)\n        meter_score, score = get_score(ans_base, gt_val)\n        print(ans_base.shape, gt_val.shape)\n        scores.append( score )\n        \n    score = np.mean(scores)\n    \n    print(params)\n    print('Validation:', score)\n    print('-' * 60)\n    \n    return score, scores, results, oofs, val_idx_list","metadata":{"execution":{"iopub.status.busy":"2021-07-29T20:55:17.233662Z","iopub.execute_input":"2021-07-29T20:55:17.233993Z","iopub.status.idle":"2021-07-29T20:55:17.251087Z","shell.execute_reply.started":"2021-07-29T20:55:17.233957Z","shell.execute_reply":"2021-07-29T20:55:17.250242Z"},"jupyter":{"source_hidden":true},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"X = train_base[train_cols].values\ny = GT.values\ngroups = train_base.collectionName.values\n\nscore, scores, oofs, models_lat, models_lng = coord_cv_oof([params, params], groups)\nscore","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"oofs_copy = oofs.copy()\noofs_copy.loc[:, ['latDeg']] = oofs_copy['latDeg'] - oofs_copy['d_latDeg']\noofs_copy.loc[:, ['lngDeg']] = oofs_copy['lngDeg'] - oofs_copy['d_lngDeg']\nprint(get_score(oofs_copy, ground_truth))\noofs_copy.to_csv('oofs_4p303.csv', index=False)","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"get_score(train_base, ground_truth)","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Test prediction","metadata":{}},{"cell_type":"code","source":"TEST_INPUT_DIR = '/kaggle/input/google-smartphone-decimeter-challenge/'\ntest_fname = 'baseline_locations_test.csv'\n\ntest_base = pd.read_csv(TEST_INPUT_DIR + test_fname)","metadata":{"execution":{"iopub.status.busy":"2021-08-04T23:19:10.011712Z","iopub.execute_input":"2021-08-04T23:19:10.012114Z","iopub.status.idle":"2021-08-04T23:19:10.145909Z","shell.execute_reply.started":"2021-08-04T23:19:10.012083Z","shell.execute_reply":"2021-08-04T23:19:10.144837Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"test_base_list = get_collections_list(test_base)\n\nlag_dfs = [get_lag_features(df, lag_columns, lag_count) for df in test_base_list]\ntest_base = pd.concat(lag_dfs, axis=0)\ntest_base.reset_index(drop=True, inplace=True)","metadata":{"execution":{"iopub.status.busy":"2021-08-04T23:19:13.79508Z","iopub.execute_input":"2021-08-04T23:19:13.795442Z","iopub.status.idle":"2021-08-04T23:19:15.58122Z","shell.execute_reply.started":"2021-08-04T23:19:13.795388Z","shell.execute_reply":"2021-08-04T23:19:15.580302Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"phone_feat_test = enc.transform(test_base.phoneName.values.reshape(-1, 1))\nphone_feat_test_df = pd.DataFrame( phone_feat_test.toarray(), columns=phone_cols )\ntest_base = pd.concat([test_base, phone_feat_test_df], axis=1)","metadata":{"execution":{"iopub.status.busy":"2021-08-04T23:19:17.572931Z","iopub.execute_input":"2021-08-04T23:19:17.573267Z","iopub.status.idle":"2021-08-04T23:19:17.684728Z","shell.execute_reply.started":"2021-08-04T23:19:17.573238Z","shell.execute_reply":"2021-08-04T23:19:17.683736Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"test_base_deltas = get_deltas(test_base)[train_cols]\nX_test = test_base_deltas[train_cols].values\nall_y_test = []\nfor model_lat, model_lng in zip(models_lat, models_lng):\n    y_test = np.array([model_lat.predict(X_test), model_lng.predict(X_test)]).T\n    all_y_test.append(y_test)","metadata":{"execution":{"iopub.status.busy":"2021-08-04T23:19:18.460951Z","iopub.execute_input":"2021-08-04T23:19:18.461307Z","iopub.status.idle":"2021-08-04T23:19:27.904319Z","shell.execute_reply.started":"2021-08-04T23:19:18.461272Z","shell.execute_reply":"2021-08-04T23:19:27.903438Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"all_y_test = np.array(all_y_test)\nall_y_test.shape","metadata":{"execution":{"iopub.status.busy":"2021-08-04T23:19:27.905977Z","iopub.execute_input":"2021-08-04T23:19:27.906382Z","iopub.status.idle":"2021-08-04T23:19:27.917999Z","shell.execute_reply.started":"2021-08-04T23:19:27.906341Z","shell.execute_reply":"2021-08-04T23:19:27.917118Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"y_test = np.median(all_y_test, axis=0)\ny_test","metadata":{"execution":{"iopub.status.busy":"2021-08-04T23:19:35.810954Z","iopub.execute_input":"2021-08-04T23:19:35.811476Z","iopub.status.idle":"2021-08-04T23:19:35.819535Z","shell.execute_reply.started":"2021-08-04T23:19:35.811444Z","shell.execute_reply":"2021-08-04T23:19:35.818525Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"test_base[['latDeg', 'lngDeg']] -= y_test\ntest_base.head(2)","metadata":{"execution":{"iopub.status.busy":"2021-08-04T23:19:55.539786Z","iopub.execute_input":"2021-08-04T23:19:55.540131Z","iopub.status.idle":"2021-08-04T23:19:55.600502Z","shell.execute_reply.started":"2021-08-04T23:19:55.540103Z","shell.execute_reply":"2021-08-04T23:19:55.599616Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"sub_num=40\ntest_base[['phone', 'millisSinceGpsEpoch', 'latDeg', 'lngDeg']].to_csv(f'subm{sub_num}.csv', index=False)\ntest_base.to_csv(f'baseline_locations_test_upd_{sub_num}.csv', index=False)","metadata":{"execution":{"iopub.status.busy":"2021-08-04T23:20:18.346228Z","iopub.execute_input":"2021-08-04T23:20:18.346611Z","iopub.status.idle":"2021-08-04T23:20:28.893446Z","shell.execute_reply.started":"2021-08-04T23:20:18.346577Z","shell.execute_reply":"2021-08-04T23:20:28.892366Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Visualization","metadata":{}},{"cell_type":"code","source":"def plot_change():\n    #TRAIN_INPUT_DIR = '/kaggle/input/gps-datasets/' \n    #train_fname = 'train_clean.csv'\n    INPUT_DIR = '/kaggle/input/gps-datasets/test_subm_boost.csv'\n    test_base_best = pd.read_csv(INPUT_DIR)\n    test_base_best.head(2)\n    \n    cname = 'MTV-1'\n    pname = 'Pixel4'\n\n    ind1 = test_base_best.collectionName.apply(lambda w : cname in w).values & test_base_best.phoneName.apply(lambda w : pname == w).values\n    collection1 = test_base_best[ind1]\n\n    ind2 = test_base.collectionName.apply(lambda w : cname in w).values & test_base.phoneName.apply(lambda w : pname == w).values\n    collection2 = test_base[ind2].copy()\n\n    collection2['phoneName'] = collection2['phoneName'].apply(lambda w : w + '_new')\n    visualize_collection( pd.concat([collection1, collection2]) )","metadata":{"execution":{"iopub.status.busy":"2021-08-04T20:57:38.980764Z","iopub.execute_input":"2021-08-04T20:57:38.981267Z","iopub.status.idle":"2021-08-04T20:57:39.359944Z","shell.execute_reply.started":"2021-08-04T20:57:38.981235Z","shell.execute_reply":"2021-08-04T20:57:39.359261Z"},"trusted":true},"execution_count":null,"outputs":[]}]}