{"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":"**[NOTE] This notebook's method is not good, so I'd like to know what information would help me reproduce the baseline.**\n\n\nWe are currently working on reproducing the baseline using the derived file, and Parrot (@hyperc) has a great notebook that can be used as a reference, but it doesn't score as well as the actual baseline.\n\nAccording to [this discussion](https://www.kaggle.com/c/google-smartphone-decimeter-challenge/discussion/238583),\n> It has 4+N states, where 4 refers to the user's position in ECEF and clock offset (x, y, z, t), and N states are inter-signal biases (ISB) for the number of non-GPS-L1 signal types. For instance, if the device measures signals of GPS L1 frequency, GLO G1 frequency, GPS L5 frequency, GAL E1 frequency at the same epoch, the number of non-GPS-L1 signal types equals 3 (i.e. N=3).\n\nAnd according to [data overview](https://www.kaggle.com/c/google-smartphone-decimeter-challenge/data),\n> [train/test]/[drive_id]/[phone_name]/[phone_name]_derived.csv - GNSS intermediate values derived from raw GNSS measurements, provided for convenience.\n> The baseline locations are computed using correctedPrM and the satellite positions, using a standard Weighted Least Squares (WLS) solver, with the phone's position (x, y, z), clock bias (t), and isrbM for each unique signal type as states for each epoch.\n\n\nSo, based on the notebook published by Parrot (@hyperc), I fixed it to also estimate isrbM for each signal type, **but the score got worse**. \n  \n**I would be grateful if you could point out what I am doing wrong in reproducing the baseline.**\n\n### Reference\n-  Parrot(@hyperc)'s notebook \nhttps://www.kaggle.com/hyperc/gsdc-reproducing-baseline-wls-on-one-measurement \n\n- Information provided by host\nhttps://www.kaggle.com/c/google-smartphone-decimeter-challenge/discussion/238583\n\n- Data Overview\nhttps://www.kaggle.com/c/google-smartphone-decimeter-challenge/data  ","metadata":{}},{"cell_type":"code","source":"import os\nimport pandas as pd\nimport numpy as np\npd.set_option('display.max_columns', 100)\nimport scipy.optimize as opt\nfrom tqdm import tqdm\nimport warnings\nwarnings.filterwarnings('ignore')\nimport multiprocessing\nfrom pathlib import Path\nimport pathlib\nfrom tqdm.notebook import tqdm\nINPUT = '../input/google-smartphone-decimeter-challenge/'\nroot = Path(INPUT)","metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","execution":{"iopub.status.busy":"2021-07-27T23:52:10.0198Z","iopub.execute_input":"2021-07-27T23:52:10.020506Z","iopub.status.idle":"2021-07-27T23:52:10.681423Z","shell.execute_reply.started":"2021-07-27T23:52:10.0204Z","shell.execute_reply":"2021-07-27T23:52:10.680351Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def ecef2lla(x, y, z):\n    # x, y and z are scalars or vectors in meters\n    x = np.array([x]).reshape(np.array([x]).shape[-1], 1)\n    y = np.array([y]).reshape(np.array([y]).shape[-1], 1)\n    z = np.array([z]).reshape(np.array([z]).shape[-1], 1)\n\n    a=6378137\n    a_sq=a**2\n    e = 8.181919084261345e-2\n    e_sq = 6.69437999014e-3\n\n    f = 1/298.257223563\n    b = a*(1-f)\n\n    # calculations:\n    r = np.sqrt(x**2 + y**2)\n    ep_sq  = (a**2-b**2)/b**2\n    ee = (a**2-b**2)\n    f = (54*b**2)*(z**2)\n    g = r**2 + (1 - e_sq)*(z**2) - e_sq*ee*2\n    c = (e_sq**2)*f*r**2/(g**3)\n    s = (1 + c + np.sqrt(c**2 + 2*c))**(1/3.)\n    p = f/(3.*(g**2)*(s + (1./s) + 1)**2)\n    q = np.sqrt(1 + 2*p*e_sq**2)\n    r_0 = -(p*e_sq*r)/(1+q) + np.sqrt(0.5*(a**2)*(1+(1./q)) - p*(z**2)*(1-e_sq)/(q*(1+q)) - 0.5*p*(r**2))\n    u = np.sqrt((r - e_sq*r_0)**2 + z**2)\n    v = np.sqrt((r - e_sq*r_0)**2 + (1 - e_sq)*z**2)\n    z_0 = (b**2)*z/(a*v)\n    h = u*(1 - b**2/(a*v))\n    phi = np.arctan((z + ep_sq*z_0)/r)\n    lambd = np.arctan2(y, x)\n\n    return phi*180/np.pi, lambd*180/np.pi, h\n\ndef calc_haversine(lat1, lon1, lat2, lon2):\n    \"\"\"Calculates the great circle distance between two points\n    on the earth. Inputs are array-like and specified in decimal degrees.\n    \"\"\"\n    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","metadata":{"execution":{"iopub.status.busy":"2021-07-27T23:52:10.682854Z","iopub.execute_input":"2021-07-27T23:52:10.683169Z","iopub.status.idle":"2021-07-27T23:52:10.700788Z","shell.execute_reply.started":"2021-07-27T23:52:10.683138Z","shell.execute_reply":"2021-07-27T23:52:10.699708Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def make_gt(recal):\n    if recal:\n        p = pathlib.Path(INPUT)\n        gt_files = list(p.glob('train/*/*/ground_truth.csv'))\n        print('ground_truth.csv count : ', len(gt_files))\n\n        gts = []\n        for gt_file in tqdm(gt_files):\n            gts.append(pd.read_csv(gt_file))\n        ground_truth = pd.concat(gts)\n        ground_truth.to_csv('gt.csv',index=False)\n    else:\n        ground_truth = pd.read_csv('gt.csv')\n    return ground_truth\n    \ngt = make_gt(recal=True)","metadata":{"execution":{"iopub.status.busy":"2021-07-27T23:52:12.481671Z","iopub.execute_input":"2021-07-27T23:52:12.482248Z","iopub.status.idle":"2021-07-27T23:52:15.288103Z","shell.execute_reply.started":"2021-07-27T23:52:12.482208Z","shell.execute_reply":"2021-07-27T23:52:15.287086Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def percentile50(x):\n    return np.percentile(x, 50)\ndef percentile95(x):\n    return np.percentile(x, 95)\n\ndef get_train_score(df, gt):\n    gt = gt.rename(columns={'latDeg':'latDeg_gt', 'lngDeg':'lngDeg_gt'})\n    df = df.merge(gt, on=['collectionName', 'phoneName', 'millisSinceGpsEpoch'], how='inner')\n    # calc_distance_error\n    df['err'] = calc_haversine(df['latDeg_gt'], df['lngDeg_gt'], df['latDeg'], df['lngDeg'])\n    # calc_evaluate_score\n    df['phone'] = df['collectionName'] + '_' + df['phoneName']\n    res = df.groupby('phone')['err'].agg([percentile50, percentile95])\n    res['p50_p90_mean'] = (res['percentile50'] + res['percentile95']) / 2 \n    score = res['p50_p90_mean'].mean()\n    return score\n\ndef get_train_score_df(df, gt):\n    gt = gt.rename(columns={'latDeg':'latDeg_gt', 'lngDeg':'lngDeg_gt'})\n    df = df.merge(gt, on=['collectionName', 'phoneName', 'millisSinceGpsEpoch'], how='left')\n    # calc_distance_error\n    df['err'] = calc_haversine(df['latDeg_gt'], df['lngDeg_gt'], df['latDeg'], df['lngDeg'])\n    return df","metadata":{"execution":{"iopub.status.busy":"2021-07-27T23:52:18.032361Z","iopub.execute_input":"2021-07-27T23:52:18.032702Z","iopub.status.idle":"2021-07-27T23:52:18.042701Z","shell.execute_reply.started":"2021-07-27T23:52:18.032673Z","shell.execute_reply":"2021-07-27T23:52:18.041514Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def gnss_log_to_dataframes(path):\n    '''Load GNSS Log'''\n    print('Loading ' + path, flush = True)\n    gnss_section_names = {'Raw', 'UncalAccel', 'UncalGyro', 'UncalMag', 'Fix', 'Status', 'OrientationDeg'}\n    with open(path) as f_open:\n        datalines = f_open.readlines()\n\n    datas = {k: [] for k in gnss_section_names}\n    gnss_map = {k: [] for k in gnss_section_names}\n    for dataline in datalines:\n        is_header = dataline.startswith('#')\n        dataline = dataline.strip('#').strip().split(',')\n        # skip over notes, version numbers, etc\n        if is_header and dataline[0] in gnss_section_names:\n            gnss_map[dataline[0]] = dataline[1:]\n        elif not is_header:\n            datas[dataline[0]].append(dataline[1:])\n\n    results = dict()\n    for k, v in datas.items():\n        results[k] = pd.DataFrame(v, columns=gnss_map[k])\n    # pandas doesn't properly infer types from these lists by default\n    for k, df in results.items():\n        for col in df.columns:\n            if col == 'CodeType':\n                continue\n            results[k][col] = pd.to_numeric(results[k][col])\n\n    return results","metadata":{"execution":{"iopub.status.busy":"2021-07-27T23:52:22.048697Z","iopub.execute_input":"2021-07-27T23:52:22.049057Z","iopub.status.idle":"2021-07-27T23:52:22.059435Z","shell.execute_reply.started":"2021-07-27T23:52:22.049013Z","shell.execute_reply":"2021-07-27T23:52:22.058433Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def apply_tips1(raw_df, derived_df):\n    # Create a new column in df_raw that corresponds to derivedの['millisSinceGpsEpoch']\n    raw_df['millisSinceGpsEpoch'] = np.floor((raw_df['TimeNanos'] - raw_df['FullBiasNanos']) / 1000000.0).astype(int)\n        \n    # Change each value in df_derived['MillisSinceGpsEpoch'] to be the prior epoch.\n    raw_timestamps = raw_df['millisSinceGpsEpoch'].unique()\n    derived_timestamps = derived_df['millisSinceGpsEpoch'].unique()\n\n    # The timestamps in derived are one epoch ahead. We need to map each epoch\n    # in derived to the prior one (in Raw).\n    indexes = np.searchsorted(raw_timestamps, derived_timestamps)\n    from_t_to_fix_derived = dict(zip(derived_timestamps, raw_timestamps[indexes-1]))\n    derived_df['millisSinceGpsEpoch'] = np.array(list(map(lambda v: from_t_to_fix_derived[v], derived_df['millisSinceGpsEpoch'])))\n    return derived_df\n\ndef apply_tips5(derived_df):\n    delta_millis = derived_df['millisSinceGpsEpoch'] - derived_df['receivedSvTimeInGpsNanos'] / 1e6\n    where_good_signals = (delta_millis > 0) & (delta_millis < 300)\n    return derived_df[where_good_signals]","metadata":{"execution":{"iopub.status.busy":"2021-07-27T23:52:39.074589Z","iopub.execute_input":"2021-07-27T23:52:39.074977Z","iopub.status.idle":"2021-07-27T23:52:39.082749Z","shell.execute_reply.started":"2021-07-27T23:52:39.074943Z","shell.execute_reply":"2021-07-27T23:52:39.081677Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"output_dir = './'\nos.makedirs(output_dir, exist_ok=True)","metadata":{"execution":{"iopub.status.busy":"2021-07-27T23:52:41.618456Z","iopub.execute_input":"2021-07-27T23:52:41.618805Z","iopub.status.idle":"2021-07-27T23:52:41.623068Z","shell.execute_reply.started":"2021-07-27T23:52:41.618773Z","shell.execute_reply":"2021-07-27T23:52:41.622256Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"base_train = pd.read_csv(root/ 'baseline_locations_train.csv')\nbase_train.loc[:,['px','py','pz']] = 0\nbase_train.head()","metadata":{"execution":{"iopub.status.busy":"2021-07-27T23:52:44.628202Z","iopub.execute_input":"2021-07-27T23:52:44.628586Z","iopub.status.idle":"2021-07-27T23:52:44.982833Z","shell.execute_reply.started":"2021-07-27T23:52:44.628541Z","shell.execute_reply":"2021-07-27T23:52:44.981815Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## WLS\ndo only SJC area to reduce calculation time","metadata":{}},{"cell_type":"code","source":"# original distance function(# Parrot's method)\ndef distance(x, **kwargs):\n    satx = kwargs[\"xSatPosMRotated\"] - x[0]\n    saty = kwargs[\"ySatPosMRotated\"] - x[1]\n    satz = kwargs[\"zSatPosMRotated\"] - x[2]\n    weight = kwargs[\"uncertaintyWeight\"]\n    prm = kwargs[\"correctedPrM\"]\n\n    d = weight * (np.sqrt(satx**2 + saty**2 +satz**2) + x[3] - prm)\n    return d\n\n# Set up least squares methods\ndef distance_v2(x, **kwargs):\n    satx = kwargs[\"xSatPosMRotated\"] - x[0]\n    saty = kwargs[\"ySatPosMRotated\"] - x[1]\n    satz = kwargs[\"zSatPosMRotated\"] - x[2]\n    weight = kwargs[\"uncertaintyWeight\"]\n    prm = kwargs[\"correctedPrM\"]\n    \n    isrbms = [k for k in kwargs.keys() if \"_isrbM\" in k]\n    N = len(isrbms)\n    isrbms_loss = 0\n    for i in range(N):\n        isrbms_loss += x[4+i] - kwargs[isrbms[i]]\n        # isrbms_loss += x[4+i]\n\n    d = weight * (np.sqrt(satx**2 + saty**2 +satz**2) + x[3] - prm + isrbms_loss)\n    return d","metadata":{"execution":{"iopub.status.busy":"2021-07-28T00:17:57.132744Z","iopub.execute_input":"2021-07-28T00:17:57.133181Z","iopub.status.idle":"2021-07-28T00:17:57.143945Z","shell.execute_reply.started":"2021-07-28T00:17:57.133141Z","shell.execute_reply":"2021-07-28T00:17:57.142972Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Parrot's method\ndef estimate_train_position_by_derived(args):\n    (collection_name, phone_name), base_df = args\n    if \"SJC\" not in collection_name:\n        return base_df\n    # Train df here only contains one collection and one measurement\n    derived_df = pd.read_csv(root / f\"train/{collection_name}/{phone_name}/{phone_name}_derived.csv\")\n    gnss_df = gnss_log_to_dataframes(str(root / f\"train/{collection_name}/{phone_name}/{phone_name}_GnssLog.txt\"))\n    raw_df = gnss_df['Raw']\n    \n    # fixed epoch\n    derived_df = apply_tips1(raw_df, derived_df)\n    derived_df = apply_tips5(derived_df)\n    derived_df = derived_df.sort_values('millisSinceGpsEpoch')\n    \n    # fixed pseudorange\n    derived_df['correctedPrM'] = derived_df.apply(lambda r: r.rawPrM + r.satClkBiasM - r.isrbM - r.ionoDelayM - r.tropoDelayM,axis=1)\n    \n    # transmission time=pseudorange/light speed\n    # diff between received time and send time\n    light_speed = 299_792_458\n    derived_df['transmissionTimeSeconds'] = derived_df['correctedPrM'] / light_speed\n\n    # Compute true sat positions at arrival time\n    omega_e = 7.2921151467e-5\n    derived_df['xSatPosMRotated'] = \\\n        np.cos(omega_e * derived_df['transmissionTimeSeconds']) * derived_df['xSatPosM'] \\\n        + np.sin(omega_e * derived_df['transmissionTimeSeconds']) * derived_df['ySatPosM']\n\n    derived_df['ySatPosMRotated'] = \\\n        - np.sin(omega_e * derived_df['transmissionTimeSeconds']) * derived_df['xSatPosM'] \\\n        + np.cos(omega_e * derived_df['transmissionTimeSeconds']) * derived_df['ySatPosM']\n    derived_df['zSatPosMRotated'] = derived_df['zSatPosM']\n    \n    # weight for WLS\n    derived_df['uncertaintyWeight'] = 1 / derived_df['rawPrUncM']\n\n    output_df = pd.DataFrame()\n    d_list = []\n    x_list = []\n    y_list = []\n    z_list = []\n    epoch_list = []\n    \n    # calc position each epoch\n    for epoch, df in derived_df.groupby('millisSinceGpsEpoch'): \n        # estimate position by WLS\n        x0 = [0]*4     \n        opt_res = opt.least_squares(distance, x0, kwargs=df.to_dict(orient=\"list\"))\n\n        # Optimiser yields a position in the ECEF coordinates\n        opt_res_pos = opt_res.x\n        d = distance(opt_res_pos, **df.to_dict(orient=\"list\"))\n\n        # ECEF position to lat/long\n        wls_estimated_pos = ecef2lla(*opt_res_pos[:3])\n        wls_estimated_pos = np.squeeze(wls_estimated_pos)\n        d_list.append(d)\n        x_list.append(wls_estimated_pos[0])\n        y_list.append(wls_estimated_pos[1])\n        z_list.append(wls_estimated_pos[2])\n        epoch_list.append(epoch)\n\n    output_df[\"latDeg\"] = x_list\n    output_df[\"lngDeg\"] = y_list\n    output_df['heightAboveWgs84EllipsoidM'] = z_list\n    output_df[\"dist\"] = d_list\n    output_df['millisSinceGpsEpoch'] = epoch_list\n    output_df['collectionName'] = collection_name\n    output_df['phoneName'] = phone_name\n\n    output_df.to_csv(output_dir + f'{collection_name}_{phone_name}_derived.csv', index=False)\n    return output_df","metadata":{"execution":{"iopub.status.busy":"2021-07-28T00:17:57.446353Z","iopub.execute_input":"2021-07-28T00:17:57.446737Z","iopub.status.idle":"2021-07-28T00:17:57.462197Z","shell.execute_reply.started":"2021-07-28T00:17:57.446702Z","shell.execute_reply":"2021-07-28T00:17:57.461143Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# fixed method(but not good)\ndef estimate_train_position_by_derived_v2(args):\n    (collection_name, phone_name), base_df = args\n    if \"SJC\" not in collection_name:\n        return base_df\n    # Train df here only contains one collection and one measurement\n    derived_df = pd.read_csv(root / f\"train/{collection_name}/{phone_name}/{phone_name}_derived.csv\")\n    gnss_df = gnss_log_to_dataframes(str(root / f\"train/{collection_name}/{phone_name}/{phone_name}_GnssLog.txt\"))\n    raw_df = gnss_df['Raw']\n    \n    # fixed epoch\n    derived_df = apply_tips1(raw_df, derived_df)\n    derived_df = apply_tips5(derived_df)\n    derived_df = derived_df.sort_values('millisSinceGpsEpoch')\n    \n    # fixed pseudorange\n    derived_df['correctedPrM'] = derived_df.apply(lambda r: r.rawPrM + r.satClkBiasM - r.isrbM - r.ionoDelayM - r.tropoDelayM,axis=1)\n    \n    # transmission time=pseudorange/light speed\n    # diff between received time and send time\n    light_speed = 299_792_458\n    derived_df['transmissionTimeSeconds'] = derived_df['correctedPrM'] / light_speed\n\n    # Compute true sat positions at arrival time\n    omega_e = 7.2921151467e-5\n    derived_df['xSatPosMRotated'] = \\\n        np.cos(omega_e * derived_df['transmissionTimeSeconds']) * derived_df['xSatPosM'] \\\n        + np.sin(omega_e * derived_df['transmissionTimeSeconds']) * derived_df['ySatPosM']\n\n    derived_df['ySatPosMRotated'] = \\\n        - np.sin(omega_e * derived_df['transmissionTimeSeconds']) * derived_df['xSatPosM'] \\\n        + np.cos(omega_e * derived_df['transmissionTimeSeconds']) * derived_df['ySatPosM']\n    derived_df['zSatPosMRotated'] = derived_df['zSatPosM']\n    \n    # weight for WLS\n    derived_df['uncertaintyWeight'] = 1 / derived_df['rawPrUncM']\n\n    output_df = pd.DataFrame()\n    d_list = []\n    x_list = []\n    y_list = []\n    z_list = []\n    epoch_list = []\n    \n    # calc position each epoch\n    for epoch, df in derived_df.groupby('millisSinceGpsEpoch'): \n        \n        ################  \n        #Fixed point ↓\n        ################\n        # the number of signal type (not GPS_L1)\n        N = len([i for i in df[\"signalType\"].unique() if i != \"GPS_L1\"])\n\n        # estimate position by WLS\n        # X -> 4 + N (4 = x,y,z,t)\n        x0 = [0]*(4 + N)  \n    \n        for signal_type in df[\"signalType\"].unique():\n            if signal_type != \"GPS_L1\":\n                df[f\"{signal_type}_isrbM\"] = 0\n                df.loc[df[\"signalType\"]==signal_type, f\"{signal_type}_isrbM\"] = df.loc[df[\"signalType\"]==signal_type, \"isrbM\"].values   \n        \n        opt_res = opt.least_squares(distance_v2, x0, kwargs=df.to_dict(orient=\"list\"))\n\n        # Optimiser yields a position in the ECEF coordinates\n        opt_res_pos = opt_res.x\n        d = distance_v2(opt_res_pos, **df.to_dict(orient=\"list\"))\n        #################\n        # Fixed point ↑\n        #################\n        \n        # ECEF position to lat/long\n        wls_estimated_pos = ecef2lla(*opt_res_pos[:3])\n        wls_estimated_pos = np.squeeze(wls_estimated_pos)\n        d_list.append(d)\n        x_list.append(wls_estimated_pos[0])\n        y_list.append(wls_estimated_pos[1])\n        z_list.append(wls_estimated_pos[2])\n        epoch_list.append(epoch)\n\n    output_df[\"latDeg\"] = x_list\n    output_df[\"lngDeg\"] = y_list\n    output_df['heightAboveWgs84EllipsoidM'] = z_list\n    output_df[\"dist\"] = d_list\n    output_df['millisSinceGpsEpoch'] = epoch_list\n    output_df['collectionName'] = collection_name\n    output_df['phoneName'] = phone_name\n\n    output_df.to_csv(output_dir + f'{collection_name}_{phone_name}_derived.csv', index=False)\n    return output_df","metadata":{"execution":{"iopub.status.busy":"2021-07-28T00:17:57.643941Z","iopub.execute_input":"2021-07-28T00:17:57.644355Z","iopub.status.idle":"2021-07-28T00:17:57.663344Z","shell.execute_reply.started":"2021-07-28T00:17:57.644317Z","shell.execute_reply":"2021-07-28T00:17:57.662248Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## original method by parrot","metadata":{}},{"cell_type":"code","source":"import multiprocessing\n\ngr = base_train.groupby(['collectionName','phoneName'])\nprocesses = multiprocessing.cpu_count()\nwith multiprocessing.Pool(processes=processes) as pool:\n    dfs = pool.imap_unordered(estimate_train_position_by_derived, gr)\n    dfs = tqdm(dfs, total=len(gr))\n    dfs = list(dfs)\nall_derived_df = pd.concat(dfs).sort_values(['collectionName', 'phoneName', 'millisSinceGpsEpoch']).reset_index(drop=True)     ","metadata":{"execution":{"iopub.status.busy":"2021-07-28T00:17:58.003282Z","iopub.execute_input":"2021-07-28T00:17:58.003791Z","iopub.status.idle":"2021-07-28T00:21:17.374123Z","shell.execute_reply.started":"2021-07-28T00:17:58.003757Z","shell.execute_reply":"2021-07-28T00:21:17.373256Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"from tqdm.notebook import tqdm\ndf_list = []\ncount = 0\nfor (collection_name, phone_name), base_df in tqdm(base_train.groupby(['collectionName','phoneName'])):\n    # break\n    if \"SJC\" in collection_name:\n        print(f\"\\n{collection_name} {phone_name}\")\n        base_df = base_df.sort_values('millisSinceGpsEpoch')\n        target_gt = gt[(gt['collectionName']==collection_name)&(gt['phoneName']==phone_name)].sort_values('millisSinceGpsEpoch').reset_index(drop=True)\n        \n        # \n        derived_df = pd.read_csv(output_dir + f'{collection_name}_{phone_name}_derived.csv')\n        derived_df = derived_df[~derived_df[\"millisSinceGpsEpoch\"].duplicated()].sort_values(\"millisSinceGpsEpoch\").reset_index(drop=True)\n        derived_df = derived_df.rename(columns={\"latDeg\":\"_latDeg\", \"lngDeg\":\"_lngDeg\"})\n        print(base_df.shape, target_gt.shape, derived_df.shape)\n        base_score = get_train_score(base_df, target_gt)\n        print(\"baseline:\", base_score)\n\n        derived_df = pd.merge_asof(base_df, derived_df[[\"millisSinceGpsEpoch\", \"_latDeg\", \"_lngDeg\"]], on=[\"millisSinceGpsEpoch\"], tolerance=10, direction='nearest')\n\n        # replace if data is nan \n        derived_df.loc[derived_df[\"_latDeg\"].isna(), \"_latDeg\"] = derived_df.loc[derived_df[\"_latDeg\"].isna(), \"latDeg\"].values\n        derived_df.loc[derived_df[\"_lngDeg\"].isna(), \"_lngDeg\"] = derived_df.loc[derived_df[\"_lngDeg\"].isna(), \"lngDeg\"].values\n\n        derived_df = derived_df.drop([\"latDeg\", \"lngDeg\"], axis=1).rename(columns={\"_latDeg\":\"latDeg\",\"_lngDeg\":\"lngDeg\"})\n        df_list.append(derived_df)\n        derived_score = get_train_score(derived_df, target_gt)\n        print(\"derived:\", derived_score)\ncorrected_base_df = pd.concat(df_list).reset_index(drop=True)","metadata":{"execution":{"iopub.status.busy":"2021-07-28T00:21:17.375424Z","iopub.execute_input":"2021-07-28T00:21:17.375838Z","iopub.status.idle":"2021-07-28T00:21:18.424461Z","shell.execute_reply.started":"2021-07-28T00:21:17.375807Z","shell.execute_reply":"2021-07-28T00:21:18.423513Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"There is a difference of ~10 from the actual score.  \nNext I'll try a fixed version that estimates isrbM for each signal type.","metadata":{}},{"cell_type":"markdown","source":"## Fixed version(but not good)","metadata":{}},{"cell_type":"code","source":"gr = base_train.groupby(['collectionName','phoneName'])\nprocesses = multiprocessing.cpu_count()\nwith multiprocessing.Pool(processes=processes) as pool:\n    dfs = pool.imap_unordered(estimate_train_position_by_derived_v2, gr)\n    dfs = tqdm(dfs, total=len(gr))\n    dfs = list(dfs)\nall_derived_df = pd.concat(dfs).sort_values(['collectionName', 'phoneName', 'millisSinceGpsEpoch']).reset_index(drop=True)     ","metadata":{"execution":{"iopub.status.busy":"2021-07-28T00:22:52.367758Z","iopub.execute_input":"2021-07-28T00:22:52.368163Z","iopub.status.idle":"2021-07-28T00:31:25.135881Z","shell.execute_reply.started":"2021-07-28T00:22:52.368127Z","shell.execute_reply":"2021-07-28T00:31:25.135092Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"from tqdm.notebook import tqdm\ndf_list = []\ncount = 0\nfor (collection_name, phone_name), base_df in tqdm(base_train.groupby(['collectionName','phoneName'])):\n    # break\n    if \"SJC\" in collection_name:\n        print(f\"\\n{collection_name} {phone_name}\")\n        base_df = base_df.sort_values('millisSinceGpsEpoch')\n        target_gt = gt[(gt['collectionName']==collection_name)&(gt['phoneName']==phone_name)].sort_values('millisSinceGpsEpoch').reset_index(drop=True)\n        \n        # \n        derived_df = pd.read_csv(output_dir + f'{collection_name}_{phone_name}_derived.csv')\n        derived_df = derived_df[~derived_df[\"millisSinceGpsEpoch\"].duplicated()].sort_values(\"millisSinceGpsEpoch\").reset_index(drop=True)\n        derived_df = derived_df.rename(columns={\"latDeg\":\"_latDeg\", \"lngDeg\":\"_lngDeg\"})\n        print(base_df.shape, target_gt.shape, derived_df.shape)\n        base_score = get_train_score(base_df, target_gt)\n        print(\"baseline:\", base_score)\n\n        derived_df = pd.merge_asof(base_df, derived_df[[\"millisSinceGpsEpoch\", \"_latDeg\", \"_lngDeg\"]], on=[\"millisSinceGpsEpoch\"], tolerance=10, direction='nearest')\n\n        # replace if data is nan \n        derived_df.loc[derived_df[\"_latDeg\"].isna(), \"_latDeg\"] = derived_df.loc[derived_df[\"_latDeg\"].isna(), \"latDeg\"].values\n        derived_df.loc[derived_df[\"_lngDeg\"].isna(), \"_lngDeg\"] = derived_df.loc[derived_df[\"_lngDeg\"].isna(), \"lngDeg\"].values\n\n        derived_df = derived_df.drop([\"latDeg\", \"lngDeg\"], axis=1).rename(columns={\"_latDeg\":\"latDeg\",\"_lngDeg\":\"lngDeg\"})\n        df_list.append(derived_df)\n        derived_score = get_train_score(derived_df, target_gt)\n        print(\"derived:\", derived_score)\ncorrected_base_df = pd.concat(df_list).reset_index(drop=True)","metadata":{"execution":{"iopub.status.busy":"2021-07-28T00:32:18.150254Z","iopub.execute_input":"2021-07-28T00:32:18.150783Z","iopub.status.idle":"2021-07-28T00:32:19.146895Z","shell.execute_reply.started":"2021-07-28T00:32:18.150749Z","shell.execute_reply":"2021-07-28T00:32:19.145832Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"The error was larger than the original one.  \nAccording to the discussion, we need to estimate isrbM to reproduce baseline, but I don't know how to do that.  \n**I would be very happy if someone could tell me**.  ","metadata":{}},{"cell_type":"code","source":"","metadata":{},"execution_count":null,"outputs":[]}]}