{"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":"## Summary\n\nThis is my first submission (LB: 4.906) program.\n\nInstead of using the Kalman filter, this program formulates location estimation as a quadratic programming problem.","metadata":{}},{"cell_type":"markdown","source":"## Libraries","metadata":{}},{"cell_type":"code","source":"%%writefile constants.py\n\nimport json\nimport datetime\nfrom collections import defaultdict\nimport numpy as np\n\nGPS_ORIGIN_DAY       = datetime.date(1980, 1, 6)\nGPS_ORIGIN_DATETIME  = datetime.datetime(1980, 1, 6)\nGLONASS_LEAP_SECONDS = 18\nBEIDOU_LEAP_SECONDS  = 14\nTZ_MSK = datetime.timezone(datetime.timedelta(hours=+3), 'MSK')\n\nWGS84_SEMI_MAJOR_AXIS = 6378137.0\nWGS84_SEMI_MINOR_AXIS = 6356752.314245\nWGS84_SQUARED_FIRST_ECCENTRICITY  = 6.69437999013e-3\nWGS84_SQUARED_SECOND_ECCENTRICITY = 6.73949674226e-3\nWGS84_FIRST_ECCENTRICITY  = np.sqrt(WGS84_SQUARED_FIRST_ECCENTRICITY)\nWGS84_SECOND_ECCENTRICITY = np.sqrt(WGS84_SQUARED_SECOND_ECCENTRICITY)\n\nLIGHT_SPEED = 299792458.0\n\nOMEGA_EARTH = 7.2921151467e-5\nMU_EARTH    = 3.986005e+14\n\nFREQ_GPS_L1  = 1.575420e+09\nFREQ_GPS_L5  = 1.176450e+09\nFREQ_GAL_E1  = FREQ_GPS_L1\nFREQ_GAL_E5A = FREQ_GPS_L5\nFREQ_QZS_J1  = FREQ_GPS_L1\nFREQ_QZS_J5  = FREQ_GPS_L5\nFREQ_BDS_B1I = 1.561098e+09\nFREQ_GLO_G1_NOMINAL = 1602.00 * 1e+6\nFREQ_GLO_G1_DELTA   = 562.5 * 1e+3\n\nCONSTELLATION_TYPE_MAP = {\n    'GPS'     : 1,\n    'GLONASS' : 3,\n    'QZSS'    : 4,\n    'BEIDOU'  : 5,\n    'GALILEO' : 6,\n}\n\nRAW_STATE_BIT_MAP = {\n     0: \"Code Lock\",\n     1: \"Bit Sync\",\n     2: \"Subframe Sync\",\n     3: \"Time Of Week Decoded State\",\n     4: \"Millisecond Ambiguity\",\n     5: \"Symbol Sync\",\n     6: \"GLONASS String Sync\",\n     7: \"GLONASS Time Of Day Decoded\",\n     8: \"BEIDOU D2 Bit Sync\",\n     9: \"BEIDOU D2 Subframe Sync\",\n    10: \"Galileo E1BC Code Lock\",\n    11: \"Galileo E1C 2^nd^ Code Lock\",\n    12: \"Galileo E1B Page Sync\",\n    13: \"SBAS Sync\",\n    14: \"Time Of Week Known\",\n    15: \"GLONASS Time Of Day Known\",\n}\nRAW_STATE_BIT_INV_MAP = { value : key for key, value in RAW_STATE_BIT_MAP.items() }\n\nSYSTEM_NAME_MAP = {\n    'GPS'     : 'G',\n    'GLONASS' : 'R',\n    'GALILEO' : 'E',\n    'BEIDOU'  : 'C',\n    'QZSS'    : 'J',\n}\n\nGLONASS_FREQ_CHANNEL_MAP = {\n    1 : 1,\n    2 : -4,\n    3 : 5,\n    4 : 6,\n    5 : 1,\n    6 : -4,\n    7 : 5,\n    8 : 6,\n    9 : -2,\n    10 : -7,\n    11 : 0,\n    12 : -1,\n    13 : -2,\n    14 : -7,\n    15 : 0,\n    16 : -1,\n    17 : 4,\n    18 : -3,\n    19 : 3,\n    20 : 2,\n    21 : 4,\n    22 : -3,\n    23 : 3,\n    24 : 2,\n}\n\nQZSS_PRN_SVID_MAP = {\n    193 : 1,\n    194 : 2,\n    199 : 3,\n    195 : 4,\n}\n\nINIT_B = np.deg2rad(  37.5)\nINIT_L = np.deg2rad(-122.2)\nINIT_H = 0.0\n\nFREQ_TOL = 100.0\nCn0DbHz_THRESHOLD = 20.0\nReceivedSvTimeUncertaintyNanos_THRESHOLD = 100\nRAW_PSEUDO_RANGE_THRESHOLD = 50_000 * 1e+3\n\nCLOCK_TIME_MARGIN = datetime.timedelta(seconds=90)\nORBIT_TIME_MARGIN = datetime.timedelta(hours=3)\nIONO_TIME_MARGIN  = datetime.timedelta(hours=2)\n\nEPSILON_M = 0.01\nELEVATION_CUTOFF = np.deg2rad(7.0)\nDEFAULT_TROPO_DELAY_M = 2.48\n\nHAVERSINE_RADIUS = 6_371_000","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2021-08-03T14:17:35.672731Z","iopub.execute_input":"2021-08-03T14:17:35.673268Z","iopub.status.idle":"2021-08-03T14:17:35.686521Z","shell.execute_reply.started":"2021-08-03T14:17:35.673183Z","shell.execute_reply":"2021-08-03T14:17:35.685641Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"%%writefile transform.py\n\nimport numpy as np\nfrom dataclasses import dataclass\n\nimport constants as C\n\n@dataclass\nclass ECEF:\n    x: np.array\n    y: np.array\n    z: np.array\n\n    def to_numpy(self):\n        return np.stack([self.x, self.y, self.z], axis=0)\n\n    @staticmethod\n    def from_numpy(pos):\n        x, y, z = [np.squeeze(w) for w in np.split(pos, 3, axis=-1)]\n        return ECEF(x=x, y=y, z=z)\n\n@dataclass\nclass BLH:\n    lat : np.array\n    lng : np.array\n    hgt : np.array\n\n@dataclass\nclass ENU:\n    east  : np.array\n    north : np.array\n    up    : np.array\n\n@dataclass\nclass AZEL:\n    elevation : np.array\n    azimuth   : np.array\n    zenith    : np.array\n\ndef BLH_to_ECEF(blh):\n    a  = C.WGS84_SEMI_MAJOR_AXIS\n    e2 = C.WGS84_SQUARED_FIRST_ECCENTRICITY\n    sin_B = np.sin(blh.lat)\n    cos_B = np.cos(blh.lat)\n    sin_L = np.sin(blh.lng)\n    cos_L = np.cos(blh.lng)\n    n = a / np.sqrt(1 - e2*sin_B**2)\n    x = (n + blh.hgt) * cos_B * cos_L\n    y = (n + blh.hgt) * cos_B * sin_L\n    z = ((1 - e2) * n + blh.hgt) * sin_B\n    return ECEF(x=x, y=y, z=z)\n\ndef ECEF_to_BLH_approximate(ecef):\n    a = C.WGS84_SEMI_MAJOR_AXIS\n    b = C.WGS84_SEMI_MINOR_AXIS\n    e2  = C.WGS84_SQUARED_FIRST_ECCENTRICITY\n    e2_ = C.WGS84_SQUARED_SECOND_ECCENTRICITY\n    x = ecef.x\n    y = ecef.y\n    z = ecef.z\n    r = np.sqrt(x**2 + y**2)\n    t = np.arctan2(z * (a/b), r)\n    B = np.arctan2(z + (e2_*b)*np.sin(t)**3, r - (e2*a)*np.cos(t)**3)\n    L = np.arctan2(y, x)\n    n = a / np.sqrt(1 - e2*np.sin(B)**2)\n    H = (r / np.cos(B)) - n\n    return BLH(lat=B, lng=L, hgt=H)\n\nECEF_to_BLH = ECEF_to_BLH_approximate\n\ndef ECEF_to_ENU(pos, base):\n    dx = pos.x - base.x\n    dy = pos.y - base.y\n    dz = pos.z - base.z\n    base_blh = ECEF_to_BLH(base)\n    sin_B = np.sin(base_blh.lat)\n    cos_B = np.cos(base_blh.lat)\n    sin_L = np.sin(base_blh.lng)\n    cos_L = np.cos(base_blh.lng)\n    e = -sin_L*dx + cos_L*dy\n    n = -sin_B*cos_L*dx - sin_B*sin_L*dy + cos_B*dz\n    u =  cos_B*cos_L*dx + cos_B*sin_L*dy + sin_B*dz\n    return ENU(east=e, north=n, up=u)\n\ndef ENU_to_AZEL(enu):\n    e = enu.east\n    n = enu.north\n    u = enu.up\n    elevation = np.arctan2(u, np.sqrt(e**2 + n**2))\n    azimuth   = np.arctan2(e, n)\n    zenith    = (0.5 * np.pi) - elevation\n    return AZEL(elevation=elevation,\n                azimuth=azimuth,\n                zenith=zenith)\n\ndef ECEF_to_AZEL(pos, base):\n    return ENU_to_AZEL(ECEF_to_ENU(pos, base))\n\ndef haversine_distance(blh_1, blh_2):\n    dlat = blh_2.lat - blh_1.lat\n    dlng = blh_2.lng - blh_1.lng\n    a = np.sin(dlat/2)**2 + np.cos(blh_1.lat) * np.cos(blh_2.lat) * np.sin(dlng/2)**2\n    dist = 2 * C.HAVERSINE_RADIUS * np.arcsin(np.sqrt(a))\n    return dist\n\ndef hubenys_distance(blh_1, blh_2):\n    Rx = C.WGS84_SEMI_MAJOR_AXIS\n    Ry = C.WGS84_SEMI_MINOR_AXIS\n    E2 = C.WGS84_SQUARED_FIRST_ECCENTRICITY\n    num_M = Rx * (1 - E2)\n    Dy = blh_1.lat - blh_2.lat\n    Dx = blh_1.lng - blh_2.lng\n    P  = 0.5 * (blh_1.lat + blh_2.lat)\n    W  = np.sqrt(1 - E2 * np.sin(P)**2)\n    M  = num_M / W**3\n    N  = Rx / W\n    d2 = (Dy * M)**2 + (Dx * N * np.cos(P))**2\n    d  = np.sqrt(d2)\n    return d\n\ndef jacobian_BLH_to_ECEF(blh):\n    a  = C.WGS84_SEMI_MAJOR_AXIS\n    e2 = C.WGS84_SQUARED_FIRST_ECCENTRICITY\n    B = blh.lat\n    L = blh.lng\n    H = blh.hgt\n    cos_B = np.cos(B)\n    sin_B = np.sin(B)\n    cos_L = np.cos(L)\n    sin_L = np.sin(L)\n    N = a / np.sqrt(1 - e2*sin_B**2)\n    dNdB = a * e2 * sin_B * cos_B * (1 - e2*sin_B**2)**(-3/2)\n    N_plus_H = N + H\n    cos_B_cos_L = cos_B * cos_L\n    cos_B_sin_L = cos_B * sin_L\n    sin_B_cos_L = sin_B * cos_L\n    sin_B_sin_L = sin_B * sin_L\n\n    dXdB = dNdB*cos_B_cos_L - N_plus_H*sin_B_cos_L\n    dYdB = dNdB*cos_B_sin_L - N_plus_H*sin_B_sin_L\n    dZdB = (1-e2)*dNdB*sin_B + (1-e2)*N_plus_H*cos_B\n\n    dXdL = - N_plus_H * cos_B_sin_L\n    dYdL =   N_plus_H * cos_B_cos_L\n    dZdL = np.zeros_like(dXdL)\n\n    dXdH = cos_B_cos_L\n    dYdH = cos_B_sin_L\n    dZdH = sin_B\n\n    J = np.stack([[dXdB, dXdL, dXdH],\n                  [dYdB, dYdL, dYdH],\n                  [dZdB, dZdL, dZdH]], axis=0)\n    axes = list(range(2, J.ndim)) + [0, 1]\n    J = np.transpose(J, axes)\n    return J\n\ndef jacobian_ECEF_to_ENU(blh):\n    B = blh.lat\n    L = blh.lng\n    cos_B = np.cos(B)\n    sin_B = np.sin(B)\n    cos_L = np.cos(L)\n    sin_L = np.sin(L)\n    \n    dEdX = -sin_L\n    dEdY =  cos_L\n    dEdZ = np.zeros_like(dEdX)\n    \n    dNdX = -sin_B*cos_L\n    dNdY = -sin_B*sin_L\n    dNdZ =  cos_B\n\n    dUdX = cos_B*cos_L\n    dUdY = cos_B*sin_L\n    dUdZ = sin_B\n\n    J = np.stack([[dEdX, dEdY, dEdZ],\n                  [dNdX, dNdY, dNdZ],\n                  [dUdX, dUdY, dUdZ]], axis=0)\n    axes = list(range(2, J.ndim)) + [0, 1]\n    J = np.transpose(J, axes)\n    return J\n\ndef pd_haversine_distance(df1, df2):\n    blh1 = BLH(\n        lat=np.deg2rad(df1['latDeg'].values),\n        lng=np.deg2rad(df1['lngDeg'].values),\n        hgt=0,\n    )\n    blh2 = BLH(\n        lat=np.deg2rad(df2['latDeg'].values),\n        lng=np.deg2rad(df2['lngDeg'].values),\n        hgt=0,\n    )\n    return haversine_distance(blh1, blh2)","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2021-08-03T14:17:35.687836Z","iopub.execute_input":"2021-08-03T14:17:35.688099Z","iopub.status.idle":"2021-08-03T14:17:35.714796Z","shell.execute_reply.started":"2021-08-03T14:17:35.688073Z","shell.execute_reply":"2021-08-03T14:17:35.713619Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"%%writefile area_prediction.py\n\nimport numpy as np\nimport pandas as pd\nfrom pathlib import Path\nfrom glob import glob\nfrom sklearn.neighbors import KNeighborsClassifier\n\nBASE_DIR = Path('../input/google-smartphone-decimeter-challenge')\n\ntrain_base = pd.read_csv(BASE_DIR / 'baseline_locations_train.csv')\ntrain_base = train_base.sort_values([\n    \"collectionName\", \"phoneName\", \"millisSinceGpsEpoch\"\n]).reset_index(drop=True)\n\ntrain_base['area'] = train_base['collectionName'].map(lambda x: x.split('-')[4])\n\ntrain_name = np.array(sorted(path.split('/')[-1] for path in glob(f'{BASE_DIR}/train/*')))\ntrain_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]\ntrain_tree     = train_name[np.array([22,23,25,26,28]) - 1]\ntrain_downtown = train_name[np.array([24,27,29]) - 1]\n\ntrain_base['area_target'] = -1\ntrain_base.loc[train_base['collectionName'].isin(train_highway),  'area_target'] = 0\ntrain_base.loc[train_base['collectionName'].isin(train_tree),     'area_target'] = 1\ntrain_base.loc[train_base['collectionName'].isin(train_downtown), 'area_target'] = 2\n\ndef 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\ntrain = processing_downtown(train_base, is_train=True)\ntrain['downtown_target'] = (train['area_target']==2).astype(int)\n\ndowntown_model_knn = KNeighborsClassifier(n_neighbors=1)\ndowntown_model_knn.fit(\n    train[['latDeg', 'lngDeg']],\n    train['downtown_target'],\n)\n\ndef 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\ntrain = processing_highway_tree(train_base, is_train=True)\n\nhighway_tree_model_knn = KNeighborsClassifier(n_neighbors=1)\nhighway_tree_model_knn.fit(\n    train.loc[train['area_target']!=2, ['latDeg', 'lngDeg']],\n    train.loc[train['area_target']!=2, 'area_target'],\n)\n\ndef 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)","metadata":{"execution":{"iopub.status.busy":"2021-08-03T14:17:35.826881Z","iopub.execute_input":"2021-08-03T14:17:35.827439Z","iopub.status.idle":"2021-08-03T14:17:35.834063Z","shell.execute_reply.started":"2021-08-03T14:17:35.827388Z","shell.execute_reply":"2021-08-03T14:17:35.833129Z"},"_kg_hide-input":true,"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## main.py","metadata":{}},{"cell_type":"code","source":"import multiprocessing\nimport glob\nimport numpy as np\nimport pandas as pd\nimport scipy.sparse\nimport scipy.sparse.linalg\nfrom tqdm.notebook import tqdm\n\nimport transform\nimport area_prediction\n\nINPUT_PATH = '../input/google-smartphone-decimeter-challenge'\n\nBASELINE_DF = pd.concat([pd.read_csv(f'{INPUT_PATH}/baseline_locations_train.csv'),\n                         pd.read_csv(f'{INPUT_PATH}/baseline_locations_test.csv'),\n                         ], axis=0)\ndef get_baseline(collection_name):\n    df = BASELINE_DF[BASELINE_DF['collectionName'] == collection_name].copy()\n    df.reset_index(drop=True, inplace=True)\n    return df\n\ndef get_optimization_constants(base_df, sigma_y):\n    const = dict()\n    dt = 1.0\n    t0 = base_df['millisSinceGpsEpoch'].min()\n    TIME_y = 1e-3 * (base_df['millisSinceGpsEpoch'] - t0).values\n    N_y = TIME_y.shape[0]\n    N_x = int(np.ceil(np.max(TIME_y) / dt) + 1)\n    const['N_y'] = N_y\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    sigma_u = 1.0\n    diag_R  = np.full(2*N_x - 2, 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 = alpha * (alpha - 1)**2\n    coeff_v1 = 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['C_orig']  = scipy.sparse.bmat(C, format='csr')\n\n    diag_L = np.full(2*N_y, sigma_y**(-2))\n    const['L_orig'] = scipy.sparse.spdiags(diag_L, [0], 2*N_y, 2*N_y, format='csr')\n\n    const['Y_orig'] = base_df[['latDeg', 'lngDeg']].values.flatten()\n    \n    return const\n\ndef solve_QP(const, valid):\n    A = const['A']\n    B = const['B']\n    R = const['R']\n    C_orig = const['C_orig']\n    L_orig = const['L_orig']\n    Y_orig = const['Y_orig']\n    valid2 = np.stack([valid, valid], axis=1).flatten()\n    C = C_orig[valid2, :]\n    L = L_orig[np.ix_(valid2, valid2)]\n    Y = Y_orig[valid2]\n    \n    BRB = B @ scipy.sparse.linalg.spsolve(R, B.T)\n    CLC = C.T @ (L @ C)\n    CLY = C.T @ (L @ Y)\n    A_sys  = scipy.sparse.bmat([[CLC, A.T], [A, -BRB]], format='csc')\n    b_sys  = np.concatenate([CLY, np.zeros(A.shape[0])], axis=0)\n    x_sys  = scipy.sparse.linalg.spsolve(A_sys, b_sys)\n    X_star = x_sys[0:A.shape[1]]\n    Y_star = C_orig @ X_star\n    return Y_star\n\ndef do_postprocess(args):\n    collection_name, params = args\n    base_df = get_baseline(collection_name)\n    const   = get_optimization_constants(base_df, params['sigma_y'])\n    valid   = np.full(const['N_y'], True)\n    for loop in range(3):\n        Y_star = solve_QP(const, valid)\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        d = transform.pd_haversine_distance(pp_df, base_df)\n        valid = (d < params['reject_m'])\n    return pp_df\n\ndef make_postprocessing_df(config):\n    args_list = []\n    for collection_list, params in config:\n        for collection_name in collection_list:\n            args_list.append((collection_name, params))\n    processes = multiprocessing.cpu_count()\n    with multiprocessing.Pool(processes=processes) as pool:\n        df_list = pool.imap_unordered(do_postprocess, args_list)\n        df_list = tqdm(df_list, total=len(args_list))\n        df_list = list(df_list)\n    output_df = pd.concat(df_list, axis=0).sort_values(['phone', 'millisSinceGpsEpoch'])\n    return output_df\n\ndef print_score(output_df):\n    score_list = []\n    for gid, phone_df in output_df.groupby('phone'):\n        drive, phone = gid.split('_')\n        gt_df = pd.read_csv(f'{INPUT_PATH}/train/{drive}/{phone}/ground_truth.csv')\n        d = transform.pd_haversine_distance(phone_df, gt_df)\n        score = np.mean([np.quantile(d, 0.50), np.quantile(d, 0.95)])\n        score_list.append(score)\n    score = np.mean(score_list)\n    print(f'train score: {score:.3f}')\n    return\n\ndef main():\n    params_highway = { 'sigma_y'  : 3.0,\n                       'reject_m' : 7.0,\n                      }\n    params_treeway = { 'sigma_y'  : 6.0,\n                       'reject_m' : 11.0,\n                      }\n    params_downtown = { 'sigma_y'  : 30.0,\n                        'reject_m' : 19.0,\n                       }\n\n    collection_list_all = np.array(sorted(path.split('/')[-1] for path in glob.glob(f'{INPUT_PATH}/train/*')))\n    collection_list_highway  = collection_list_all[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    collection_list_treeway  = collection_list_all[np.array([22,23,25,26,28]) - 1]\n    collection_list_downtown = collection_list_all[np.array([24,27,29]) - 1]\n    config = [\n        (collection_list_highway,  params_highway),\n        (collection_list_treeway,  params_treeway),\n        (collection_list_downtown, params_downtown),\n    ]\n    train_pp_df = make_postprocessing_df(config)\n    print_score(train_pp_df)\n    \n    test_base = pd.read_csv(f'{INPUT_PATH}/baseline_locations_test.csv')\n    collection_list_highway, collection_list_treeway, collection_list_downtown = area_prediction.predict_area(test_base)\n    config = [\n        (collection_list_highway,  params_highway),\n        (collection_list_treeway,  params_treeway),\n        (collection_list_downtown, params_downtown),\n    ]\n    test_pp_df = make_postprocessing_df(config)\n\n    train_pp_df.to_csv('smoothing_1st_train.csv', index=False)\n    test_pp_df.to_csv('smoothing_1st_test.csv', index=False)\n\n    columns = ['phone', 'millisSinceGpsEpoch', 'latDeg', 'lngDeg']\n    sub_df = test_pp_df[columns]\n    sub_df.to_csv('submission.csv', index=False)\n    return","metadata":{"execution":{"iopub.status.busy":"2021-08-03T14:17:35.835586Z","iopub.execute_input":"2021-08-03T14:17:35.836124Z","iopub.status.idle":"2021-08-03T14:17:38.394474Z","shell.execute_reply.started":"2021-08-03T14:17:35.836094Z","shell.execute_reply":"2021-08-03T14:17:38.393546Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"main()","metadata":{"execution":{"iopub.status.busy":"2021-08-03T14:17:38.39592Z","iopub.execute_input":"2021-08-03T14:17:38.396168Z","iopub.status.idle":"2021-08-03T14:19:04.923945Z","shell.execute_reply.started":"2021-08-03T14:17:38.396143Z","shell.execute_reply":"2021-08-03T14:19:04.922086Z"},"trusted":true},"execution_count":null,"outputs":[]}]}