{"metadata":{"kernelspec":{"language":"python","display_name":"Python 3","name":"python3"},"language_info":{"pygments_lexer":"ipython3","nbconvert_exporter":"python","version":"3.6.4","file_extension":".py","codemirror_mode":{"name":"ipython","version":3},"name":"python","mimetype":"text/x-python"}},"nbformat_minor":4,"nbformat":4,"cells":[{"cell_type":"code","source":"! pip3 install simdkalman\n\nimport pandas as pd\nimport numpy as np\nimport simdkalman\n\nfrom pathlib import Path\nfrom tqdm import tqdm\n\npd.set_option('expand_frame_repr', False)\n\n\ndef get_union(dir_name, file_name):\n    print(\"Reading\", file_name, \"from\", dir_name, \"...\")\n    datas = []\n    for path in tqdm(Path(dir_name).rglob(file_name)):\n        datas.append(pd.read_csv(path, low_memory=False, index_col=False))\n    return pd.concat(datas)\n\n\ndef calc_haversine(lat1, lng1, lat2, lng2):\n    \"\"\"from\n    https://www.kaggle.com/dehokanta/baseline-post-processing-by-outlier-correction\n    \"\"\"\n    RADIUS = 6_367_000\n    lat1, lng1, lat2, lng2 = map(np.radians, [lat1, lng1, lat2, lng2])\n    dlat = lat2 - lat1\n    dlng = lng2 - lng1\n    a = np.sin(dlat / 2) ** 2 + \\\n        np.cos(lat1) * np.cos(lat2) * np.sin(dlng / 2) ** 2\n    dist = 2 * RADIUS * np.arcsin(a ** 0.5)\n    return dist\n\n\nclass KF:\n    def __init__(self):\n        T = 1.0\n        state_transition = np.array([[1, 0, T, 0, 0.5 * T ** 2, 0],\n                                     [0, 1, 0, T, 0, 0.5 * T ** 2],\n                                     [0, 0, 1, 0, T, 0],\n                                     [0, 0, 0, 1, 0, T],\n                                     [0, 0, 0, 0, 1, 0],\n                                     [0, 0, 0, 0, 0, 1]])\n        process_noise = np.diag([1e-5, 1e-5, 5e-6, 5e-6, 1e-6, 1e-6]) + np.ones((6, 6)) * 1e-9\n        observation_model = np.array([[1, 0, 0, 0, 0, 0], [0, 1, 0, 0, 0, 0]])\n        observation_noise = np.diag([5e-5, 5e-5]) + np.ones((2, 2)) * 1e-9\n        self.kf = simdkalman.KalmanFilter(\n            state_transition=state_transition,\n            process_noise=process_noise,\n            observation_model=observation_model,\n            observation_noise=observation_noise)\n\n    def apply_kf_smoothing(self, df, suffix):\n        phones = df[\"phone\"].drop_duplicates().tolist()\n        for phone in tqdm(phones):\n            cond = df['phone'] == phone\n            tmp = df[cond].copy()\n            tmp[0] = tmp[\"millisSinceGpsEpoch\"] // 1000\n            tmp = tmp.merge(pd.DataFrame(range(tmp[0].min(), tmp[0].max() + 1)), on=[0], how=\"right\")\n            tmp_np = tmp[['latDeg', 'lngDeg']].to_numpy()\n            nan_idxs = tmp[tmp[\"millisSinceGpsEpoch\"].isnull()].index.to_list()\n            tmp_np = tmp_np.reshape(1, len(tmp_np), 2)\n            smoothed = self.kf.smooth(tmp_np).states.mean\n            smoothed = np.delete(smoothed, list(nan_idxs), 1)\n            df.loc[cond, 'latDeg' + suffix] = smoothed[0, :, 0]\n            df.loc[cond, 'lngDeg' + suffix] = smoothed[0, :, 1]","metadata":{"execution":{"iopub.status.busy":"2021-07-31T06:09:16.557076Z","iopub.execute_input":"2021-07-31T06:09:16.557655Z","iopub.status.idle":"2021-07-31T06:09:25.021369Z","shell.execute_reply.started":"2021-07-31T06:09:16.55755Z","shell.execute_reply":"2021-07-31T06:09:25.020425Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def parse_gnsslog(data_dir, output_dir):\n    header_of_table_named = {\n        \"UncalMag\":\n            \"utcTimeMillis,elapsedRealtimeNanos,UncalMagXMicroT,UncalMagYMicroT,UncalMagZMicroT\\n\",\n        \"UncalAccel\":\n            \"utcTimeMillis,elapsedRealtimeNanos,UncalAccelXMps2,UncalAccelYMps2,UncalAccelZMps2\\n\",\n        \"UncalGyro\":\n            \"utcTimeMillis,elapsedRealtimeNanos,UncalGyroXRadPerSec,UncalGyroYRadPerSec,UncalGyroZRadPerSec\\n\",\n        \"Status\":\n            \"UnixTimeMillis,SignalCount,SignalIndex,ConstellationType,Svid,CarrierFrequencyHz,Cn0DbHz,\"\n            \"AzimuthDegrees,ElevationDegrees,UsedInFix,HasAlmanacData,HasEphemerisData\\n\",\n        \"Raw\":\n            \"utcTimeMillis,TimeNanos,LeapSecond,TimeUncertaintyNanos,FullBiasNanos,BiasNanos,BiasUncertaintyNanos,\"\n            \"DriftNanosPerSecond,DriftUncertaintyNanosPerSecond,HardwareClockDiscontinuityCount,\"\n            \"Svid,TimeOffsetNanos,State,ReceivedSvTimeNanos,ReceivedSvTimeUncertaintyNanos,Cn0DbHz,\"\n            \"PseudorangeRateMetersPerSecond,PseudorangeRateUncertaintyMetersPerSecond,AccumulatedDeltaRangeState,\"\n            \"AccumulatedDeltaRangeMeters,AccumulatedDeltaRangeUncertaintyMeters,CarrierFrequencyHz,CarrierCycles,\"\n            \"CarrierPhase,CarrierPhaseUncertainty,MultipathIndicator,SnrInDb,ConstellationType,AgcDb\\n\",\n        \"Fix\":\n            \"Provider,LatitudeDegrees,LongitudeDegrees,AltitudeMeters,SpeedMps,AccuracyMeters,BearingDegrees,\"\n            \"UnixTimeMillis,SpeedAccuracyMps,BearingAccuracyDegrees\\n\"\n    }\n    for part in [\"train\", \"test\"]:\n        print(part, \"gnsslog parsing ...\")\n        for file_name in Path(data_dir / part).rglob(\"*GnssLog.txt\"):\n            print(file_name)\n            with open(str(file_name)) as f_open:\n                datalines = f_open.readlines()\n            f_named = {}\n            for f_name in header_of_table_named.keys():\n                dir_pth = Path('/'.join(str(file_name.parent).split(\"/\")[-3:]))\n                Path(dir_pth.parent.parent).mkdir(exist_ok=True)\n                Path(dir_pth.parent).mkdir(exist_ok=True)\n                Path(dir_pth).mkdir(exist_ok=True)\n                output_filename = dir_pth / (f_name + \".csv\")\n                if not (output_filename).exists():\n                    f_named[f_name] = open(output_filename, \"w\")\n            if f_named:\n                collectionName = file_name.parent.parent.name\n                phoneName = file_name.parent.name\n                for f_name, f in f_named.items():\n                    f.write(\"collectionName,phoneName,\" + header_of_table_named[f_name])\n                for dataline in tqdm(datalines):\n                    for f_name in f_named.keys():\n                        if dataline.startswith(f_name):\n                            f_named[f_name].write(collectionName)\n                            f_named[f_name].write(',')\n                            f_named[f_name].write(phoneName)\n                            f_named[f_name].write(',')\n                            f_named[f_name].write(dataline[len(f_name) + 1:])\n                            break\n                for f in f_named.values():\n                    f.close()\n    print(\"finish\")","metadata":{"execution":{"iopub.status.busy":"2021-07-31T06:09:25.023042Z","iopub.execute_input":"2021-07-31T06:09:25.023389Z","iopub.status.idle":"2021-07-31T06:09:25.037188Z","shell.execute_reply.started":"2021-07-31T06:09:25.02335Z","shell.execute_reply":"2021-07-31T06:09:25.035644Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"DATA_DIR = Path(\"../input/google-smartphone-decimeter-challenge\")\nOUTPUT_DIR = Path(\"./\")\n\n# 读取\ngt_data = get_union(DATA_DIR / \"train\", \"ground_truth.csv\")\ntrain_data = pd.read_csv(DATA_DIR / \"baseline_locations_train.csv\")\ntest_data = pd.read_csv(DATA_DIR / \"baseline_locations_test.csv\")\n\nparse_gnsslog(DATA_DIR, OUTPUT_DIR)\n\n# gt合并\ngt_data.rename(columns={\n    \"latDeg\": \"latDeg_truth\",\n    \"lngDeg\": \"lngDeg_truth\",\n    \"heightAboveWgs84EllipsoidM\": \"heightAboveWgs84EllipsoidM_truth\"}, inplace=True)\ntrain_data = train_data.merge(gt_data, on=['collectionName', 'phoneName', 'millisSinceGpsEpoch'])\ndel gt_data\n\n# 计算baseline和truth的距离\ntrain_data[\"dist_between_baseline_and_truth\"] = calc_haversine(\n    train_data.latDeg, train_data.lngDeg,\n    train_data.latDeg_truth, train_data.lngDeg_truth\n)\n\n# 手机类别\nidx_of_phonename = {phoneName: idx for idx, phoneName in enumerate(test_data.phoneName.drop_duplicates().to_list())}\ntrain_data[\"phoneCat\"] = train_data[\"phoneName\"].map(idx_of_phonename)\ntest_data[\"phoneCat\"] = test_data[\"phoneName\"].map(idx_of_phonename)\n\n\ndef data_processing(data):\n    # 时间周期性特征\n    data[[\"month\", \"day\"]] = data['collectionName'].str.split('-', expand=True)[[1, 2]]\n    data[\"ratio_of_year\"] = (30 * (data.month.astype(int) - 1) + data.day.astype(int)) / 365\n\n    # 相对于相邻帧平均值的差值的绝对值\n    offsets = [1, 2, 3, 4, 5]\n    groups = data[[\"phone\", \"latDeg\", \"lngDeg\"]].groupby(\"phone\")\n    for offset in offsets:\n        data[[\"mean_latDeg+-\" + str(offset), \"mean_lngDeg+-\" + str(offset)]] = \\\n            groups.rolling(window=3, min_periods=2, center=True).mean().values\n        data[\"latDeg_mean_delta_\" + str(offset)] = np.abs(data[\"mean_latDeg+-\" + str(offset)] - data[\"latDeg\"])\n        data[\"lngDeg_mean_delta_\" + str(offset)] = np.abs(data[\"mean_lngDeg+-\" + str(offset)] - data[\"lngDeg\"])\n\n    # 添加前后差值特征\n    data[[\"latDeg-1\", \"lngDeg-1\"]] = groups.shift(1)\n    data[[\"latDeg+1\", \"lngDeg+1\"]] = groups.shift(-1)\n    data[\"latDeg_pre_increment\"] = data[\"latDeg\"] - data[\"latDeg-1\"]\n    data[\"lngDeg_pre_increment\"] = data[\"lngDeg\"] - data[\"lngDeg-1\"]\n    data[\"latDeg_post_increment\"] = data[\"latDeg+1\"] - data[\"latDeg\"]\n    data[\"lngDeg_post_increment\"] = data[\"lngDeg+1\"] - data[\"lngDeg\"]\n    data[\"dist_pre\"] = calc_haversine(data[\"latDeg\"], data[\"lngDeg\"], data[\"latDeg-1\"], data[\"lngDeg-1\"])\n    data[\"dist_post\"] = calc_haversine(data[\"latDeg\"], data[\"lngDeg\"], data[\"latDeg+1\"], data[\"lngDeg+1\"])\n\n    # 没有前一帧，用后一帧来代替\n    nan_idxs = data[data.latDeg_pre_increment.isnull()].index\n    data.loc[nan_idxs, [\"latDeg_pre_increment\", \"lngDeg_pre_increment\", \"dist_pre\"]] = \\\n        data.loc[nan_idxs + 1, [\"latDeg_pre_increment\", \"lngDeg_pre_increment\", \"dist_pre\"]].values\n    nan_idxs = data[data.latDeg_post_increment.isnull()].index\n    data.loc[nan_idxs, [\"latDeg_post_increment\", \"lngDeg_post_increment\", \"dist_post\"]] = \\\n        data.loc[nan_idxs - 1, [\"latDeg_post_increment\", \"lngDeg_post_increment\", \"dist_post\"]].values\n\n    # 添加前后差值绝对均值特征\n    data[\"latDeg_pre_post_mean_abs_delta\"] = (np.abs(data[\"latDeg_pre_increment\"]) +\n                                              np.abs(data[\"latDeg_post_increment\"])) / 2\n    data[\"lngDeg_pre_post_mean_abs_delta\"] = (np.abs(data[\"lngDeg_pre_increment\"]) +\n                                              np.abs(data[\"lngDeg_post_increment\"])) / 2\n    data[\"pre_post_mean_abs_dist\"] = (np.abs(data[\"dist_pre\"]) + np.abs(data[\"dist_post\"])) / 2\n\n    # 同一时间不同设备的统计数值\n    col_names = [\"latDeg\", \"lngDeg\", \"heightAboveWgs84EllipsoidM\"]\n    fn_names = [\"max\", \"min\", \"mean\", \"sum\", \"count\"]\n    groups = data.groupby(\"millisSinceGpsEpoch\")\n    for fn_name_ in fn_names:\n        data[[fn_name_ + '_' + col_name + \"_with_same_millisSinceGpsEpoch\" for col_name in col_names]] = \\\n            groups[col_names].transform(fn_name_).values\n\n    # 与同时间所有手机位置均值的距离\n    data[\"dist_between_baseline_and_phone_mean\"] = calc_haversine(data.latDeg, data.lngDeg,\n                                                                  data.mean_latDeg_with_same_millisSinceGpsEpoch,\n                                                                  data.mean_lngDeg_with_same_millisSinceGpsEpoch)\n\ndata_processing(train_data)\ndata_processing(test_data)\n\ntrain_derived_data = get_union(DATA_DIR / \"train\", \"*_derived.csv\")\ntest_derived_data = get_union(DATA_DIR / \"test\", \"*_derived.csv\")\ndef derived_data_processing(derived_data, data):\n    # 添加correctedPrm特征\n    derived_data[\"correctedPrm\"] = derived_data[\"rawPrM\"] + derived_data[\"satClkBiasM\"] - \\\n                                   derived_data[\"isrbM\"] - derived_data[\"ionoDelayM\"] - \\\n                                   derived_data[\"tropoDelayM\"]\n\n    # 信号类型与卫星id\n    derived_data[\"signalType_svid\"] = derived_data[\"signalType\"] + '_' + derived_data[\"svid\"].astype(\"string\")\n\n    # 时间对齐\n    data[\"millisSinceGpsEpoch/1000_round\"] = np.round(data[\"millisSinceGpsEpoch\"] / 1000).astype(np.int64)\n    derived_data[\"millisSinceGpsEpoch/1000_round\"] = np.round(derived_data[\"millisSinceGpsEpoch\"] / 1000).astype(\n        np.int64)\n\n    # data和derived_data合并\n    groups = derived_data.groupby([\"collectionName\", \"phoneName\", \"millisSinceGpsEpoch/1000_round\"])\n    data = data.merge(\n        groups[\"signalType_svid\"].agg(lambda group: group.values),\n        on=[\"collectionName\", \"phoneName\", \"millisSinceGpsEpoch/1000_round\"], how=\"left\"\n    )\n    data = data.merge(\n        groups[[\"correctedPrm\", \"rawPrUncM\", \"satClkDriftMps\"]].mean().rename(columns={\n            \"correctedPrm\": \"correctedPrm_avg\",\n            \"rawPrUncM\": \"rawPrUncM_avg\",\n            \"satClkDriftMps\": \"satClkDriftMps_avg\"\n        }),\n        on=[\"collectionName\", \"phoneName\", \"millisSinceGpsEpoch/1000_round\"], how=\"left\"\n    )\n    data = data.merge(\n        groups[[\"correctedPrm\", \"rawPrUncM\", \"satClkDriftMps\"]].std().rename(columns={\n            \"correctedPrm\": \"correctedPrm_std\",\n            \"rawPrUncM\": \"rawPrUncM_std\",\n            \"satClkDriftMps\": \"satClkDriftMps_std\"\n        }),\n        on=[\"collectionName\", \"phoneName\", \"millisSinceGpsEpoch/1000_round\"], how=\"left\"\n    )\n\n    return data\n\ntrain_data = derived_data_processing(train_derived_data, train_data)\ntest_data = derived_data_processing(test_derived_data, test_data)\ndel train_derived_data\ndel test_derived_data\n\n# 分桶去null\nnum_buckets = 20\nattr_names = [\"correctedPrm\", \"rawPrUncM\", \"satClkDriftMps\"]\nconcat_data = pd.concat([train_data, test_data])\nfor fn_name in [\"avg\", \"std\"]:\n    for attr_name in attr_names:\n        concat_data[attr_name + '_' + fn_name + \"_bucketized\"] = pd.qcut(concat_data[attr_name + '_' + fn_name],\n                                                                  num_buckets, labels=False)\n        concat_data.loc[concat_data[attr_name + '_' + fn_name + \"_bucketized\"].isnull(),\n                 attr_name + '_' + fn_name + \"_bucketized\"] = num_buckets\n        concat_data[attr_name + '_' + fn_name + \"_bucketized\"] = \\\n            concat_data[attr_name + '_' + fn_name + \"_bucketized\"].astype(np.int64)\n        train_data[attr_name + '_' + fn_name + \"_bucketized\"] = \\\n            concat_data[:len(train_data)][attr_name + '_' + fn_name + \"_bucketized\"]\n        test_data[attr_name + '_' + fn_name + \"_bucketized\"] = \\\n            concat_data[-len(test_data):][attr_name + '_' + fn_name + \"_bucketized\"]\n\ntrain_raw_data = get_union(OUTPUT_DIR / \"train\", \"Raw.csv\")\ntest_raw_data = get_union(OUTPUT_DIR / \"test\", \"Raw.csv\")\ndef raw_data_processing(raw_data, data):\n    raw_data[\"millisSinceGpsEpoch\"] = np.round((raw_data.TimeNanos - raw_data.FullBiasNanos) / 1000000).\\\n        astype(np.int64)\n    raw_data[\"millisSinceGpsEpoch/1000_round\"] = np.round(raw_data[\"millisSinceGpsEpoch\"] / 1000).astype(np.int64)\n    groups = raw_data.groupby([\"collectionName\", \"phoneName\", \"millisSinceGpsEpoch/1000_round\"])\n    data = data.merge(\n        groups[[\"Cn0DbHz\", \"BiasUncertaintyNanos\"]].mean().rename(columns={\n            \"Cn0DbHz\": \"Cn0DbHz_avg\",\n            \"BiasUncertaintyNanos\": \"BiasUncertaintyNanos_avg\",\n        }),\n        on=[\"collectionName\", \"phoneName\", \"millisSinceGpsEpoch/1000_round\"], how=\"left\"\n    )\n    data = data.merge(\n        groups[[\"Cn0DbHz\", \"BiasUncertaintyNanos\"]].std().rename(columns={\n            \"Cn0DbHz\": \"Cn0DbHz_std\",\n            \"BiasUncertaintyNanos\": \"BiasUncertaintyNanos_std\",\n        }),\n        on=[\"collectionName\", \"phoneName\", \"millisSinceGpsEpoch/1000_round\"], how=\"left\"\n    )\n    data = data.merge(\n        groups[\"utcTimeMillis\"].first(),\n        on=[\"collectionName\", \"phoneName\", \"millisSinceGpsEpoch/1000_round\"], how=\"left\"\n    )\n    return data\n\ntrain_data = raw_data_processing(train_raw_data, train_data)\ntest_data = raw_data_processing(test_raw_data, test_data)\ndel train_raw_data\ndel test_raw_data\n\ntrain_gyro_data = get_union(OUTPUT_DIR / \"train\", \"UncalGyro.csv\")\ntest_gyro_data = get_union(OUTPUT_DIR / \"test\", \"UncalGyro.csv\")\ntrain_mag_data = get_union(OUTPUT_DIR / \"train\", \"UncalMag.csv\")\ntest_mag_data = get_union(OUTPUT_DIR / \"test\", \"UncalMag.csv\")\ntrain_accel_data = get_union(OUTPUT_DIR / \"train\", \"UncalAccel.csv\")\ntest_accel_data = get_union(OUTPUT_DIR / \"test\", \"UncalAccel.csv\")\ndef sensor_data_processing(gyro_data, mag_data, accel_data, data):\n    data[\"utcTimeMillis/1000_round\"] = np.round(data.utcTimeMillis / 1000).astype(np.int64)\n    axises = ['X', 'Y', 'Z']\n    prefix_suffix_df_dicts = [\n        {\"prefix\": \"UncalAccel\", \"suffix\": \"Mps2\", \"sensor_data\": accel_data},\n        {\"prefix\": \"UncalGyro\", \"suffix\": \"RadPerSec\", \"sensor_data\": gyro_data},\n        {\"prefix\": \"UncalMag\", \"suffix\": \"MicroT\", \"sensor_data\": mag_data}\n    ]\n    for prefix_suffix_df_dict in prefix_suffix_df_dicts:\n        prefix = prefix_suffix_df_dict[\"prefix\"]\n        suffix = prefix_suffix_df_dict[\"suffix\"]\n        sensor_data = prefix_suffix_df_dict[\"sensor_data\"]\n        sensor_data[\"utcTimeMillis\"] = sensor_data.utcTimeMillis.astype(np.int64)\n        sensor_data[\"utcTimeMillis/1000_round\"] = np.round(sensor_data.utcTimeMillis / 1000).astype(np.int64)\n        groups = sensor_data.groupby([\"collectionName\", \"phoneName\", \"utcTimeMillis/1000_round\"])\n        data = data.merge(groups[[prefix + axis + suffix for axis in axises]].mean().rename(columns={\n            prefix + axis + suffix: prefix + axis + suffix + \"_avg\" for axis in axises\n        }), on=[\"collectionName\", \"phoneName\", \"utcTimeMillis/1000_round\"], how=\"left\")\n        data = data.merge(groups[[prefix + axis + suffix for axis in axises]].std().rename(columns={\n            prefix + axis + suffix: prefix + axis + suffix + \"_std\" for axis in axises\n        }), on=[\"collectionName\", \"phoneName\", \"utcTimeMillis/1000_round\"], how=\"left\")\n\n    col_names = [\n        \"UncalAccelXMps2_avg\",\n        \"UncalAccelYMps2_avg\",\n        \"UncalAccelZMps2_avg\",\n        \"UncalAccelXMps2_std\",\n        \"UncalAccelYMps2_std\",\n        \"UncalAccelZMps2_std\",\n        \"UncalGyroXRadPerSec_avg\",\n        \"UncalGyroYRadPerSec_avg\",\n        \"UncalGyroZRadPerSec_avg\",\n        \"UncalGyroXRadPerSec_std\",\n        \"UncalGyroYRadPerSec_std\",\n        \"UncalGyroZRadPerSec_std\",\n        \"UncalMagXMicroT_avg\",\n        \"UncalMagYMicroT_avg\",\n        \"UncalMagZMicroT_avg\",\n        \"UncalMagXMicroT_std\",\n        \"UncalMagYMicroT_std\",\n        \"UncalMagZMicroT_std\"\n    ]\n    groups = data.groupby([\"collectionName\", \"phoneName\"])\n    data[[col_name + \"+1\" for col_name in col_names]] = groups[col_names].shift(-1)\n\n    nan_idxs = data[data[\"UncalAccelXMps2_avg+1\"].isnull()].index\n    data.loc[nan_idxs, [col_name + \"+1\" for col_name in col_names]] = data.loc[nan_idxs, col_names].values\n\n    data[\"ratio\"] = (data[\"utcTimeMillis\"] % 1000) / 1000\n    data[col_names] = data[col_names].values * (1 - data[\"ratio\"]).values.reshape(-1, 1) + \\\n        data[[col_name + \"+1\" for col_name in col_names]].values * data[\"ratio\"].values.reshape(-1, 1)\n\n    return data\n\ntrain_data = sensor_data_processing(train_gyro_data, train_mag_data, train_accel_data, train_data)\ntest_data = sensor_data_processing(test_gyro_data, test_mag_data, test_accel_data, test_data)\ndel train_gyro_data\ndel test_gyro_data\ndel train_mag_data\ndel test_mag_data\ndel train_accel_data\ndel test_accel_data\n\n# 卡尔曼相关特征\nkf = KF()\nkf.apply_kf_smoothing(train_data, suffix=\"_kf\")\nkf.apply_kf_smoothing(test_data, suffix=\"_kf\")\n\ndef add_kf_features(data):\n    data[\"dist_between_base_and_kf\"] = calc_haversine(data.latDeg, data.lngDeg, data.latDeg_kf, data.lngDeg_kf)\n    data[\"base_kf_lat_delta\"] = data[\"latDeg\"] - data[\"latDeg_kf\"]\n    data[\"base_kf_lng_delta\"] = data[\"lngDeg\"] - data[\"lngDeg_kf\"]\n    data[\"abs_base_kf_lng_delta\"] = abs(data[\"lngDeg\"] - data[\"lngDeg_kf\"])\n    data[\"abs_base_kf_lat_delta\"] = abs(data[\"latDeg\"] - data[\"latDeg_kf\"])\n\nadd_kf_features(train_data)\nadd_kf_features(test_data)\n\n# 缺失值均值填充\ncol_names_ = [\n    \"UncalAccelXMps2_avg\",\n    \"UncalAccelYMps2_avg\",\n    \"UncalAccelZMps2_avg\",\n    \"UncalAccelXMps2_std\",\n    \"UncalAccelYMps2_std\",\n    \"UncalAccelZMps2_std\",\n    \"UncalGyroXRadPerSec_avg\",\n    \"UncalGyroYRadPerSec_avg\",\n    \"UncalGyroZRadPerSec_avg\",\n    \"UncalGyroXRadPerSec_std\",\n    \"UncalGyroYRadPerSec_std\",\n    \"UncalGyroZRadPerSec_std\",\n    \"UncalMagXMicroT_avg\",\n    \"UncalMagYMicroT_avg\",\n    \"UncalMagZMicroT_avg\",\n    \"UncalMagXMicroT_std\",\n    \"UncalMagYMicroT_std\",\n    \"UncalMagZMicroT_std\"\n]\ntrain_data[col_names_] = train_data[col_names_].fillna(train_data[col_names_].mean())\ntest_data[col_names_] = test_data[col_names_].fillna(test_data[col_names_].mean())\n\n# 切分训练验证\ntrain_pd = train_data[\n    (train_data.collectionName == \"2020-05-14-US-MTV-1\") |\n    (train_data.collectionName == \"2020-05-21-US-MTV-2\") |\n    (train_data.collectionName == \"2020-05-29-US-MTV-1\") |\n    (train_data.collectionName == \"2020-06-04-US-MTV-1\") |\n    (train_data.collectionName == \"2020-06-05-US-MTV-1\") |\n    (train_data.collectionName == \"2020-07-08-US-MTV-1\") |\n    (train_data.collectionName == \"2020-07-17-US-MTV-1\") |\n    (train_data.collectionName == \"2020-08-06-US-MTV-2\") |\n    (train_data.collectionName == \"2020-09-04-US-SF-1\") |\n    (train_data.collectionName == \"2021-01-04-US-RWC-2\") |\n    (train_data.collectionName == \"2021-01-05-US-SVL-1\") |\n    (train_data.collectionName == \"2021-03-10-US-SVL-1\") |\n    (train_data.collectionName == \"2021-04-22-US-SJC-1\") |\n    (train_data.collectionName == \"2021-04-28-US-MTV-1\") |\n    (train_data.collectionName == \"2021-04-29-US-SJC-2\")\n]\nval_pd = train_data[\n    (train_data.collectionName == \"2020-05-14-US-MTV-2\") |\n    (train_data.collectionName == \"2020-05-21-US-MTV-1\") |\n    (train_data.collectionName == \"2020-05-29-US-MTV-2\") |\n    (train_data.collectionName == \"2020-06-05-US-MTV-2\") |\n    (train_data.collectionName == \"2020-06-11-US-MTV-1\") |\n    (train_data.collectionName == \"2020-07-17-US-MTV-2\") |\n    (train_data.collectionName == \"2020-08-03-US-MTV-1\") |\n    (train_data.collectionName == \"2020-09-04-US-SF-2\") |\n    (train_data.collectionName == \"2021-01-04-US-RWC-1\") |\n    (train_data.collectionName == \"2021-01-05-US-SVL-2\") |\n    (train_data.collectionName == \"2021-04-15-US-MTV-1\") |\n    (train_data.collectionName == \"2021-04-26-US-SVL-1\") |\n    (train_data.collectionName == \"2021-04-28-US-SJC-1\") |\n    (train_data.collectionName == \"2021-04-29-US-MTV-1\")\n]\n\n# 保存\ntrain_pd.to_csv(\"trainset_pandas.csv\", index=False)\nval_pd.to_csv(\"valset_pandas.csv\", index=False)\ntrain_data.to_csv(\"total_trainset_pandas.csv\", index=False)\ntest_data.to_csv(\"testset_pandas.csv\", index=False)","metadata":{"execution":{"iopub.status.busy":"2021-07-31T06:09:25.040351Z","iopub.execute_input":"2021-07-31T06:09:25.040821Z","iopub.status.idle":"2021-07-31T06:16:17.890378Z","shell.execute_reply.started":"2021-07-31T06:09:25.040776Z","shell.execute_reply":"2021-07-31T06:16:17.889009Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"from scipy.interpolate import interp1d\n\n\ndef mean_with_other_phones(df):\n    \"\"\"\n    https://www.kaggle.com/bpetrb/adaptive-gauss-phone-mean\n    \"\"\"\n    collections_list = df[['collectionName']].drop_duplicates().to_numpy()\n\n    for collection in collections_list:\n        phone_list = df[df['collectionName'].to_list() == collection][['phoneName']].drop_duplicates().to_numpy()\n\n        phone_data = {}\n        corrections = {}\n        for phone in phone_list:\n            cond = np.logical_and(df['collectionName'] == collection[0], df['phoneName'] == phone[0]).to_list()\n            phone_data[phone[0]] = df[cond][['millisSinceGpsEpoch', 'latDeg', 'lngDeg']].to_numpy()\n\n        for current in phone_data:\n            correction = np.ones(phone_data[current].shape, dtype=np.float)\n            correction[:, 1:] = phone_data[current][:, 1:]\n\n            # Telephones data don't complitely match by time, so - interpolate.\n            for other in phone_data:\n                if other == current:\n                    continue\n\n                loc = interp1d(phone_data[other][:, 0],\n                               phone_data[other][:, 1:],\n                               axis=0,\n                               kind='linear',\n                               copy=False,\n                               bounds_error=None,\n                               fill_value='extrapolate',\n                               assume_sorted=True)\n\n                start_idx = 0\n                stop_idx = 0\n                for idx, val in enumerate(phone_data[current][:, 0]):\n                    if val < phone_data[other][0, 0]:\n                        start_idx = idx\n                    if val < phone_data[other][-1, 0]:\n                        stop_idx = idx\n\n                if stop_idx - start_idx > 0:\n                    correction[start_idx:stop_idx, 0] += 1\n                    correction[start_idx:stop_idx, 1:] += loc(phone_data[current][start_idx:stop_idx, 0])\n\n            correction[:, 1] /= correction[:, 0]\n            correction[:, 2] /= correction[:, 0]\n\n            corrections[current] = correction.copy()\n\n        for phone in phone_list:\n            cond = np.logical_and(df['collectionName'] == collection[0], df['phoneName'] == phone[0]).to_list()\n\n            df.loc[cond, ['latDeg', 'lngDeg']] = corrections[phone[0]][:, 1:]\n\n    return df","metadata":{"execution":{"iopub.status.busy":"2021-07-31T06:16:17.892417Z","iopub.execute_input":"2021-07-31T06:16:17.892715Z","iopub.status.idle":"2021-07-31T06:16:18.265256Z","shell.execute_reply.started":"2021-07-31T06:16:17.892684Z","shell.execute_reply":"2021-07-31T06:16:18.264389Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"import lightgbm as lgb\n\nfeature_cols = [\n     'latDeg',\n     'lngDeg',\n     'heightAboveWgs84EllipsoidM',\n     'phoneCat',\n     'ratio_of_year',\n     'mean_latDeg+-1',\n     'mean_lngDeg+-1',\n     'latDeg_mean_delta_1',\n     'lngDeg_mean_delta_1',\n     'mean_latDeg+-2',\n     'mean_lngDeg+-2',\n     'latDeg_mean_delta_2',\n     'lngDeg_mean_delta_2',\n     'mean_latDeg+-3',\n     'mean_lngDeg+-3',\n     'latDeg_mean_delta_3',\n     'lngDeg_mean_delta_3',\n     'mean_latDeg+-4',\n     'mean_lngDeg+-4',\n     'latDeg_mean_delta_4',\n     'lngDeg_mean_delta_4',\n     'mean_latDeg+-5',\n     'mean_lngDeg+-5',\n     'latDeg_mean_delta_5',\n     'lngDeg_mean_delta_5',\n     'latDeg_pre_increment',\n     'lngDeg_pre_increment',\n     'latDeg_post_increment',\n     'lngDeg_post_increment',\n     'dist_pre',\n     'dist_post',\n     'latDeg_pre_post_mean_abs_delta',\n     'lngDeg_pre_post_mean_abs_delta',\n     'pre_post_mean_abs_dist',\n     'max_latDeg_with_same_millisSinceGpsEpoch',\n     'max_lngDeg_with_same_millisSinceGpsEpoch',\n     'max_heightAboveWgs84EllipsoidM_with_same_millisSinceGpsEpoch',\n     'min_latDeg_with_same_millisSinceGpsEpoch',\n     'min_lngDeg_with_same_millisSinceGpsEpoch',\n     'min_heightAboveWgs84EllipsoidM_with_same_millisSinceGpsEpoch',\n     'mean_latDeg_with_same_millisSinceGpsEpoch',\n     'mean_lngDeg_with_same_millisSinceGpsEpoch',\n     'mean_heightAboveWgs84EllipsoidM_with_same_millisSinceGpsEpoch',\n     'sum_latDeg_with_same_millisSinceGpsEpoch',\n     'sum_lngDeg_with_same_millisSinceGpsEpoch',\n     'sum_heightAboveWgs84EllipsoidM_with_same_millisSinceGpsEpoch',\n     'count_latDeg_with_same_millisSinceGpsEpoch',\n     'dist_between_baseline_and_phone_mean',\n     'correctedPrm_avg_bucketized',\n     'rawPrUncM_avg_bucketized',\n     'satClkDriftMps_avg_bucketized',\n     'correctedPrm_std_bucketized',\n     'rawPrUncM_std_bucketized',\n     'satClkDriftMps_std_bucketized',\n     'Cn0DbHz_avg',\n     'BiasUncertaintyNanos_avg',\n     'Cn0DbHz_std',\n     'BiasUncertaintyNanos_std',\n     'UncalAccelXMps2_avg',\n     'UncalAccelYMps2_avg',\n     'UncalAccelZMps2_avg',\n     'UncalAccelXMps2_std',\n     'UncalAccelYMps2_std',\n     'UncalAccelZMps2_std',\n     'UncalGyroXRadPerSec_avg',\n     'UncalGyroYRadPerSec_avg',\n     'UncalGyroZRadPerSec_avg',\n     'UncalGyroXRadPerSec_std',\n     'UncalGyroYRadPerSec_std',\n     'UncalGyroZRadPerSec_std',\n     'UncalMagXMicroT_avg',\n     'UncalMagYMicroT_avg',\n     'UncalMagZMicroT_avg',\n     'UncalMagXMicroT_std',\n     'UncalMagYMicroT_std',\n     'UncalMagZMicroT_std',\n     'latDeg_kf',\n     'lngDeg_kf',\n     'dist_between_base_and_kf',\n     'base_kf_lat_delta',\n     'base_kf_lng_delta',\n     'abs_base_kf_lng_delta',\n     'abs_base_kf_lat_delta'\n]\n\nxtr = train_pd[feature_cols]\nytr = train_pd[\"dist_between_baseline_and_truth\"]\nxval = val_pd[feature_cols]\nyval = val_pd[\"dist_between_baseline_and_truth\"]\n\nparams = {\n    'objective': 'mae',\n    'max_bin': 600,\n    'learning_rate': 0.02,\n    'num_leaves': 80\n}\n\nlgb_train = lgb.Dataset(xtr, ytr)\nlgb_eval = lgb.Dataset(xval, yval, reference=lgb_train)\nmodel = lgb.train(\n    params, lgb_train,\n    valid_sets=[lgb_train, lgb_eval],\n    verbose_eval=25,\n    num_boost_round=10000,\n    early_stopping_rounds=10\n)\n\nkf = KF()\n\ntest_pd = test_data\nxtest = test_pd[feature_cols]\nytest = model.predict(xtest)\n\nthreshold = 13\ntmp = test_pd.copy()\ntmp.loc[np.where(ytest > threshold)[0].tolist(), [\"latDeg\", \"lngDeg\"]] = np.nan\nkf.apply_kf_smoothing(tmp, \"\")\nsmoothed_tmp = mean_with_other_phones(tmp)\ndata_sub = pd.read_csv(DATA_DIR / \"sample_submission.csv\")\ndf_sub = data_sub[['phone', 'millisSinceGpsEpoch']].merge(\n    smoothed_tmp[['phone', 'millisSinceGpsEpoch', 'latDeg', 'lngDeg']], on=['phone', 'millisSinceGpsEpoch'],\n    how='inner')\ndf_sub.to_csv(OUTPUT_DIR / \"submission.csv\", index=False)","metadata":{"execution":{"iopub.status.busy":"2021-07-31T06:16:18.266959Z","iopub.execute_input":"2021-07-31T06:16:18.267386Z","iopub.status.idle":"2021-07-31T06:17:10.225684Z","shell.execute_reply.started":"2021-07-31T06:16:18.267339Z","shell.execute_reply":"2021-07-31T06:17:10.224431Z"},"trusted":true},"execution_count":null,"outputs":[]}]}