{"metadata":{"kernelspec":{"language":"python","display_name":"Python 3","name":"python3"},"language_info":{"name":"python","version":"3.10.13","mimetype":"text/x-python","codemirror_mode":{"name":"ipython","version":3},"pygments_lexer":"ipython3","nbconvert_exporter":"python","file_extension":".py"},"kaggle":{"accelerator":"none","dataSources":[{"sourceId":60095,"databundleVersionId":6542333,"sourceType":"competition"}],"dockerImageVersionId":30673,"isInternetEnabled":true,"language":"python","sourceType":"notebook","isGpuEnabled":false}},"nbformat_minor":4,"nbformat":4,"cells":[{"cell_type":"code","source":"# This Python 3 environment comes with many helpful analytics libraries installed\n# It is defined by the kaggle/python Docker image: https://github.com/kaggle/docker-python\n# For example, here's several helpful packages to load\n\nimport numpy as np # linear algebra\nimport pandas as pd # data processing, CSV file I/O (e.g. pd.read_csv)\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\nimport os\nfor dirname, _, filenames in os.walk('/kaggle/input/smartphone-decimeter-2023/train'):\n    for filename in filenames:\n        \n#         print(os.path.join(dirname, filename))\n            print(dirname)\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":{"execution":{"iopub.status.busy":"2024-03-26T16:37:04.170045Z","iopub.execute_input":"2024-03-26T16:37:04.170471Z","iopub.status.idle":"2024-03-26T16:37:05.444334Z","shell.execute_reply.started":"2024-03-26T16:37:04.170435Z","shell.execute_reply":"2024-03-26T16:37:05.44303Z"},"trusted":true},"execution_count":1,"outputs":[]},{"cell_type":"code","source":"!pip install nvector\n\nimport nvector as nv\nfrom nvector import rad, deg\n\nlat_EA, lon_EA = rad(37.416314), rad(-122.0804662)\nn_EA_E = nv.lat_lon2n_E(lat_EA, lon_EA)\nprint(n_EA_E)","metadata":{"execution":{"iopub.status.busy":"2024-03-26T16:40:44.100438Z","iopub.execute_input":"2024-03-26T16:40:44.10112Z","iopub.status.idle":"2024-03-26T16:41:02.943177Z","shell.execute_reply.started":"2024-03-26T16:40:44.101085Z","shell.execute_reply":"2024-03-26T16:41:02.94177Z"},"trusted":true},"execution_count":2,"outputs":[{"name":"stdout","text":"Collecting nvector\n  Downloading nvector-0.7.7-py2.py3-none-any.whl.metadata (43 kB)\n\u001b[2K     \u001b[90m━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━\u001b[0m \u001b[32m43.8/43.8 kB\u001b[0m \u001b[31m631.0 kB/s\u001b[0m eta \u001b[36m0:00:00\u001b[0m \u001b[36m0:00:01\u001b[0m\n\u001b[?25hRequirement already satisfied: numpy>=1.9 in /opt/conda/lib/python3.10/site-packages (from nvector) (1.26.4)\nRequirement already satisfied: scipy>=0.9 in /opt/conda/lib/python3.10/site-packages (from nvector) (1.11.4)\nRequirement already satisfied: geographiclib in /opt/conda/lib/python3.10/site-packages (from nvector) (2.0)\nRequirement already satisfied: matplotlib in /opt/conda/lib/python3.10/site-packages (from nvector) (3.7.5)\nRequirement already satisfied: cartopy in /opt/conda/lib/python3.10/site-packages (from nvector) (0.22.0)\nRequirement already satisfied: shapely>=1.7 in /opt/conda/lib/python3.10/site-packages (from cartopy->nvector) (1.8.5.post1)\nRequirement already satisfied: packaging>=20 in /opt/conda/lib/python3.10/site-packages (from cartopy->nvector) (21.3)\nRequirement already satisfied: pyshp>=2.1 in /opt/conda/lib/python3.10/site-packages (from cartopy->nvector) (2.3.1)\nRequirement already satisfied: pyproj>=3.1.0 in /opt/conda/lib/python3.10/site-packages (from cartopy->nvector) (3.6.1)\nRequirement already satisfied: contourpy>=1.0.1 in /opt/conda/lib/python3.10/site-packages (from matplotlib->nvector) (1.2.0)\nRequirement already satisfied: cycler>=0.10 in /opt/conda/lib/python3.10/site-packages (from matplotlib->nvector) (0.12.1)\nRequirement already satisfied: fonttools>=4.22.0 in /opt/conda/lib/python3.10/site-packages (from matplotlib->nvector) (4.47.0)\nRequirement already satisfied: kiwisolver>=1.0.1 in /opt/conda/lib/python3.10/site-packages (from matplotlib->nvector) (1.4.5)\nRequirement already satisfied: pillow>=6.2.0 in /opt/conda/lib/python3.10/site-packages (from matplotlib->nvector) (9.5.0)\nRequirement already satisfied: pyparsing>=2.3.1 in /opt/conda/lib/python3.10/site-packages (from matplotlib->nvector) (3.1.1)\nRequirement already satisfied: python-dateutil>=2.7 in /opt/conda/lib/python3.10/site-packages (from matplotlib->nvector) (2.9.0.post0)\nRequirement already satisfied: certifi in /opt/conda/lib/python3.10/site-packages (from pyproj>=3.1.0->cartopy->nvector) (2024.2.2)\nRequirement already satisfied: six>=1.5 in /opt/conda/lib/python3.10/site-packages (from python-dateutil>=2.7->matplotlib->nvector) (1.16.0)\nDownloading nvector-0.7.7-py2.py3-none-any.whl (73 kB)\n\u001b[2K   \u001b[90m━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━\u001b[0m \u001b[32m73.7/73.7 kB\u001b[0m \u001b[31m1.1 MB/s\u001b[0m eta \u001b[36m0:00:00\u001b[0ma \u001b[36m0:00:01\u001b[0m\n\u001b[?25hInstalling collected packages: nvector\nSuccessfully installed nvector-0.7.7\n[[-0.42182948]\n [-0.67296336]\n [ 0.60760201]]\n","output_type":"stream"}]},{"cell_type":"code","source":"!pip install pymap3d\n!pip install pykalman\n\nimport pykalman\nimport numpy as np\nimport pandas as pd\nimport pymap3d as pm\nimport pymap3d.vincenty as pmv\nimport matplotlib.pyplot as plt\nimport glob as gl\nimport scipy.optimize\nimport scipy\nfrom tqdm.auto import tqdm\nfrom scipy.interpolate import InterpolatedUnivariateSpline\nfrom scipy.spatial import distance\n\n# Constants\nCLIGHT = 299_792_458   # speed of light (m/s)\nRE_WGS84 = 6_378_137   # earth semimajor axis (WGS84) (m)\nOMGE = 7.2921151467E-5  # earth angular velocity (IS-GPS) (rad/s)","metadata":{"execution":{"iopub.status.busy":"2024-03-26T16:41:44.216673Z","iopub.execute_input":"2024-03-26T16:41:44.217147Z","iopub.status.idle":"2024-03-26T16:42:16.64396Z","shell.execute_reply.started":"2024-03-26T16:41:44.217106Z","shell.execute_reply":"2024-03-26T16:42:16.642614Z"},"trusted":true},"execution_count":3,"outputs":[{"name":"stdout","text":"Collecting pymap3d\n  Downloading pymap3d-3.1.0-py3-none-any.whl.metadata (8.7 kB)\nDownloading pymap3d-3.1.0-py3-none-any.whl (60 kB)\n\u001b[2K   \u001b[90m━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━\u001b[0m \u001b[32m60.6/60.6 kB\u001b[0m \u001b[31m790.9 kB/s\u001b[0m eta \u001b[36m0:00:00\u001b[0m--:--\u001b[0m\n\u001b[?25hInstalling collected packages: pymap3d\nSuccessfully installed pymap3d-3.1.0\nRequirement already satisfied: pykalman in /opt/conda/lib/python3.10/site-packages (0.9.5)\n","output_type":"stream"}]},{"cell_type":"code","source":"def satellite_selection(df, column):\n    \"\"\"\n    Args:\n        df : DataFrame from device_gnss.csv\n        column : Column name\n    Returns:\n        df: DataFrame with eliminated satellite signals\n    \"\"\"\n    idx = df[column].notnull()\n    idx &= df['CarrierErrorHz'] < 2.0e6  # carrier frequency error (Hz)\n    idx &= df['SvElevationDegrees'] > 10.0  # elevation angle (deg)\n    idx &= df['Cn0DbHz'] > 15.0  # C/N0 (dB-Hz)\n    idx &= df['MultipathIndicator'] == 0 # Multipath flag\n    \n    idx &= df['PseudorangeUncertainty'] < 150\n\n    return df[idx]","metadata":{"execution":{"iopub.status.busy":"2024-03-26T16:42:30.537007Z","iopub.execute_input":"2024-03-26T16:42:30.537502Z","iopub.status.idle":"2024-03-26T16:42:30.547356Z","shell.execute_reply.started":"2024-03-26T16:42:30.537464Z","shell.execute_reply":"2024-03-26T16:42:30.5458Z"},"trusted":true},"execution_count":4,"outputs":[]},{"cell_type":"code","source":"# Compute line-of-sight vector from user to satellite\ndef los_vector(xusr, xsat):\n    \"\"\"\n    Args:\n        xusr : user position in ECEF (m)\n        xsat : satellite position in ECEF (m)\n    Returns:\n        u: unit line-of-sight vector in ECEF (m)\n        rng: distance between user and satellite (m)\n    \"\"\"\n    u = xsat - xusr\n    rng = np.linalg.norm(u, axis=1).reshape(-1, 1)\n    u /= rng\n    \n    return u, rng.reshape(-1)\n\n\n# Compute Jacobian matrix\ndef jac_pr_residuals(x, xsat, pr, W):\n    \"\"\"\n    Args:\n        x : current position in ECEF (m)\n        xsat : satellite position in ECEF (m)\n        pr : pseudorange (m)\n        W : weight matrix\n    Returns:\n        W*J : Jacobian matrix\n    \"\"\"\n    u, _ = los_vector(x[:3], xsat)\n    J = np.hstack([-u, np.ones([len(pr), 1])])  # J = [-ux -uy -uz 1]\n\n    return W @ J\n\n\n# Compute pseudorange residuals\ndef pr_residuals(x, xsat, pr, W):\n    \"\"\"\n    Args:\n        x : current position in ECEF (m)\n        xsat : satellite position in ECEF (m)\n        pr : pseudorange (m)\n        W : weight matrix\n    Returns:\n        residuals*W : pseudorange residuals\n    \"\"\"\n    u, rng = los_vector(x[:3], xsat)\n\n    # Approximate correction of the earth rotation (Sagnac effect) often used in GNSS positioning\n    rng += OMGE * (xsat[:, 0] * x[1] - xsat[:, 1] * x[0]) / CLIGHT\n\n    # Add GPS L1 clock offset\n    residuals = rng - (pr - x[3])\n\n    return residuals @ W\n\n\n# Compute Jacobian matrix\ndef jac_prr_residuals(v, vsat, prr, x, xsat, W):\n    \"\"\"\n    Args:\n        v : current velocity in ECEF (m/s)\n        vsat : satellite velocity in ECEF (m/s)\n        prr : pseudorange rate (m/s)\n        x : current position in ECEF (m)\n        xsat : satellite position in ECEF (m)\n        W : weight matrix\n    Returns:\n        W*J : Jacobian matrix\n    \"\"\"\n    u, _ = los_vector(x[:3], xsat)\n    J = np.hstack([-u, np.ones([len(prr), 1])])\n\n    return W @ J\n\n\n# Compute pseudorange rate residuals\ndef prr_residuals(v, vsat, prr, x, xsat, W):\n    \"\"\"\n    Args:\n        v : current velocity in ECEF (m/s)\n        vsat : satellite velocity in ECEF (m/s)\n        prr : pseudorange rate (m/s)\n        x : current position in ECEF (m)\n        xsat : satellite position in ECEF (m)\n        W : weight matrix\n    Returns:\n        residuals*W : pseudorange rate residuals\n    \"\"\"\n    u, rng = los_vector(x[:3], xsat)\n    rate = np.sum((vsat-v[:3])*u, axis=1) \\\n          + OMGE / CLIGHT * (vsat[:, 1] * x[0] + xsat[:, 1] * v[0]\n                           - vsat[:, 0] * x[1] - xsat[:, 0] * v[1])\n\n    residuals = rate - (prr - v[3])\n\n    return residuals @ W\n\n# Carrier smoothing of pseudarange\ndef carrier_smoothing(gnss_df):\n    \"\"\"\n    Args:\n        df : DataFrame from device_gnss.csv\n    Returns:\n        df: DataFrame with carrier-smoothing pseudorange 'pr_smooth'\n    \"\"\"\n    carr_th = 1.5 # carrier phase jump threshold [m] ** 2.0 -> 1.5 **\n    pr_th =  20.0 # pseudorange jump threshold [m]\n\n    prsmooth = np.full_like(gnss_df['RawPseudorangeMeters'], np.nan)\n    # Loop for each signal\n    for (i, (svid_sigtype, df)) in enumerate((gnss_df.groupby(['Svid', 'SignalType']))):\n        df = df.replace(\n            {'AccumulatedDeltaRangeMeters': {0: np.nan}})  # 0 to NaN\n\n        # Compare time difference between pseudorange/carrier with Doppler\n        drng1 = df['AccumulatedDeltaRangeMeters'].diff() - df['PseudorangeRateMetersPerSecond']\n        drng2 = df['RawPseudorangeMeters'].diff() - df['PseudorangeRateMetersPerSecond']\n\n        # Check cycle-slip\n        slip1 = (df['AccumulatedDeltaRangeState'].to_numpy() & 2**1) != 0  # reset flag\n        slip2 = (df['AccumulatedDeltaRangeState'].to_numpy() & 2**2) != 0  # cycle-slip flag\n        slip3 = np.fabs(drng1.to_numpy()) > carr_th # Carrier phase jump\n        slip4 = np.fabs(drng2.to_numpy()) > pr_th # Pseudorange jump\n\n        idx_slip = slip1 | slip2 | slip3 | slip4\n        idx_slip[0] = True\n\n        # groups with continuous carrier phase tracking\n        df['group_slip'] = np.cumsum(idx_slip)\n\n        # Psudorange - carrier phase\n        df['dpc'] = df['RawPseudorangeMeters'] - df['AccumulatedDeltaRangeMeters']\n\n        # Absolute distance bias of carrier phase\n        meandpc = df.groupby('group_slip')['dpc'].mean()\n        df = df.merge(meandpc, on='group_slip', suffixes=('', '_Mean'))\n\n        # Index of original gnss_df\n        idx = (gnss_df['Svid'] == svid_sigtype[0]) & (\n            gnss_df['SignalType'] == svid_sigtype[1])\n\n        # Carrier phase + bias\n        prsmooth[idx] = df['AccumulatedDeltaRangeMeters'] + df['dpc_Mean']\n\n    # If carrier smoothing is not possible, use original pseudorange\n    idx_nan = np.isnan(prsmooth)\n    prsmooth[idx_nan] = gnss_df['RawPseudorangeMeters'][idx_nan]\n    gnss_df['pr_smooth'] = prsmooth\n\n    return gnss_df","metadata":{"execution":{"iopub.status.busy":"2024-03-26T16:43:07.149107Z","iopub.execute_input":"2024-03-26T16:43:07.149608Z","iopub.status.idle":"2024-03-26T16:43:07.176957Z","shell.execute_reply.started":"2024-03-26T16:43:07.149572Z","shell.execute_reply":"2024-03-26T16:43:07.175284Z"},"trusted":true},"execution_count":5,"outputs":[]},{"cell_type":"code","source":"# Compute distance by Vincenty's formulae\ndef vincenty_distance(llh1, llh2):\n    \"\"\"\n    Args:\n        llh1 : [latitude,longitude] (deg)\n        llh2 : [latitude,longitude] (deg)\n    Returns:\n        d : distance between llh1 and llh2 (m)\n    \"\"\"\n    d, az = np.array(pmv.vdist(llh1[:, 0], llh1[:, 1], llh2[:, 0], llh2[:, 1]))\n\n    return d\n\n\n# Compute score\ndef calc_score(llh, llh_gt):\n    \"\"\"\n    Args:\n        llh : [latitude,longitude] (deg)\n        llh_gt : [latitude,longitude] (deg)\n    Returns:\n        score : (m)\n    \"\"\"\n    d = vincenty_distance(llh, llh_gt)\n    score = np.mean([np.quantile(d, 0.50), np.quantile(d, 0.95)])\n\n    return score","metadata":{"execution":{"iopub.status.busy":"2024-03-26T16:43:26.398354Z","iopub.execute_input":"2024-03-26T16:43:26.398845Z","iopub.status.idle":"2024-03-26T16:43:26.406982Z","shell.execute_reply.started":"2024-03-26T16:43:26.398808Z","shell.execute_reply":"2024-03-26T16:43:26.405684Z"},"trusted":true},"execution_count":6,"outputs":[]},{"cell_type":"code","source":"# GNSS single point positioning using pseudorange\ndef point_positioning(gnss_df):\n    # Add nominal frequency to each signal\n    # Note: GLONASS is an FDMA signal, so each satellite has a different frequency\n    CarrierFrequencyHzRef = gnss_df.groupby(['Svid', 'SignalType'])[\n        'CarrierFrequencyHz'].median()\n    gnss_df = gnss_df.merge(CarrierFrequencyHzRef, how='left', on=[\n                            'Svid', 'SignalType'], suffixes=('', 'Ref'))\n    gnss_df['CarrierErrorHz'] = np.abs(\n        (gnss_df['CarrierFrequencyHz'] - gnss_df['CarrierFrequencyHzRef']))\n    \n    gnss_df['PseudorangeUncertainty'] = scipy.constants.c * (1e-9 * gnss_df['ReceivedSvTimeUncertaintyNanos'])\n\n    # Carrier smoothing\n    gnss_df = carrier_smoothing(gnss_df)\n\n    # GNSS single point positioning\n    utcTimeMillis = gnss_df['utcTimeMillis'].unique()\n    nepoch = len(utcTimeMillis)\n    x0 = np.zeros(4)  # [x,y,z,tGPSL1]\n    v0 = np.zeros(4)  # [vx,vy,vz,dtGPSL1]\n    x_wls = np.full([nepoch, 3], np.nan)  # For saving position\n    v_wls = np.full([nepoch, 3], np.nan)  # For saving velocity\n    cov_x = np.full([nepoch, 3, 3], np.nan) # For saving position covariance\n    cov_v = np.full([nepoch, 3, 3], np.nan) # For saving velocity covariance\n\n    # Loop for epochs\n    for i, (t_utc, df) in enumerate(tqdm(gnss_df.groupby('utcTimeMillis'), total=nepoch)):\n        # Valid satellite selection\n        df_pr = satellite_selection(df, 'pr_smooth')\n        df_prr = satellite_selection(df, 'PseudorangeRateMetersPerSecond')\n\n        # Corrected pseudorange/pseudorange rate\n        pr = (df_pr['pr_smooth'] + df_pr['SvClockBiasMeters'] - df_pr['IsrbMeters'] -\n              df_pr['IonosphericDelayMeters'] - df_pr['TroposphericDelayMeters']).to_numpy()\n        prr = (df_prr['PseudorangeRateMetersPerSecond'] +\n               df_prr['SvClockDriftMetersPerSecond']).to_numpy()\n\n        # Satellite position/velocity\n        xsat_pr = df_pr[['SvPositionXEcefMeters', 'SvPositionYEcefMeters',\n                         'SvPositionZEcefMeters']].to_numpy()\n        xsat_prr = df_prr[['SvPositionXEcefMeters', 'SvPositionYEcefMeters',\n                           'SvPositionZEcefMeters']].to_numpy()\n        vsat = df_prr[['SvVelocityXEcefMetersPerSecond', 'SvVelocityYEcefMetersPerSecond',\n                       'SvVelocityZEcefMetersPerSecond']].to_numpy()\n\n        # Weight matrix for peseudorange/pseudorange rate\n        Wx = np.diag(1 / df_pr['RawPseudorangeUncertaintyMeters'].to_numpy())\n        Wv = np.diag(1 / df_prr['PseudorangeRateUncertaintyMetersPerSecond'].to_numpy())\n\n        # Robust WLS requires accurate initial values for convergence,\n        # so perform normal WLS for the first time\n        if len(df_pr) >= 4:\n            # Normal WLS\n            if np.all(x0 == 0):\n                opt = scipy.optimize.least_squares(\n                    pr_residuals, x0, jac_pr_residuals, args=(xsat_pr, pr, Wx))\n                x0 = opt.x \n            # Robust WLS for position estimation\n            opt = scipy.optimize.least_squares(\n                 pr_residuals, x0, jac_pr_residuals, args=(xsat_pr, pr, Wx), loss='soft_l1')\n            if opt.status < 1 or opt.status == 2:\n                print(f'i = {i} position lsq status = {opt.status}')\n            else:\n                # Covariance estimation\n                cov = np.linalg.inv(opt.jac.T @ Wx @ opt.jac)\n                cov_x[i, :, :] = cov[:3, :3]\n                x_wls[i, :] = opt.x[:3]\n                x0 = opt.x\n                 \n        # Velocity estimation\n        if len(df_prr) >= 4:\n            if np.all(v0 == 0): # Normal WLS\n                opt = scipy.optimize.least_squares(\n                    prr_residuals, v0, jac_prr_residuals, args=(vsat, prr, x0, xsat_prr, Wv))\n                v0 = opt.x\n            # Robust WLS for velocity estimation\n            opt = scipy.optimize.least_squares(\n                prr_residuals, v0, jac_prr_residuals, args=(vsat, prr, x0, xsat_prr, Wv), loss='soft_l1')\n            if opt.status < 1:\n                print(f'i = {i} velocity lsq status = {opt.status}')\n            else:\n                # Covariance estimation\n                cov = np.linalg.inv(opt.jac.T @ Wv @ opt.jac)\n                cov_v[i, :, :] = cov[:3, :3]\n                v_wls[i, :] = opt.x[:3]\n                v0 = opt.x\n\n    return utcTimeMillis, x_wls, v_wls, cov_x, cov_v\n\n# Simple outlier detection and interpolation\ndef exclude_interpolate_outlier(x_wls, v_wls, cov_x, cov_v):\n    # Up velocity / height threshold\n    v_up_th = 2.6  # m/s  2.0 -> 2.6\n    height_th = 200.0 # m\n    v_out_sigma = 3.0 # m/s\n    x_out_sigma = 30.0 # m\n    \n    # Coordinate conversion\n    x_llh = np.array(pm.ecef2geodetic(x_wls[:, 0], x_wls[:, 1], x_wls[:, 2])).T\n    x_llh_mean = np.nanmean(x_llh, axis=0)\n    v_enu = np.array(pm.ecef2enuv(\n        v_wls[:, 0], v_wls[:, 1], v_wls[:, 2], x_llh_mean[0], x_llh_mean[1])).T\n\n    # Up velocity jump detection\n    # Cars don't jump suddenly!\n    idx_v_out = np.abs(v_enu[:, 2]) > v_up_th\n    idx_v_out |= np.isnan(v_enu[:, 2])\n    v_wls[idx_v_out, :] = np.nan\n    cov_v[idx_v_out] = v_out_sigma**2 * np.eye(3)\n    print(f'Number of velocity outliers {np.count_nonzero(idx_v_out)}')\n\n    # Height check\n    hmedian = np.nanmedian(x_llh[:, 2])\n    idx_x_out = np.abs(x_llh[:, 2] - hmedian) > height_th\n    idx_x_out |= np.isnan(x_llh[:, 2])\n    x_wls[idx_x_out, :] = np.nan\n    cov_x[idx_x_out] = x_out_sigma**2 * np.eye(3)\n    print(f'Number of position outliers {np.count_nonzero(idx_x_out)}')\n\n    # Interpolate NaNs at beginning and end of array\n    x_df = pd.DataFrame({'x': x_wls[:, 0], 'y': x_wls[:, 1], 'z': x_wls[:, 2]})\n    x_df = x_df.interpolate(limit_area='outside', limit_direction='both')\n\n    # Interpolate all NaN data\n    v_df = pd.DataFrame({'x': v_wls[:, 0], 'y': v_wls[:, 1], 'z': v_wls[:, 2]})\n    v_df = v_df.interpolate(limit_area='outside', limit_direction='both')\n    v_df = v_df.interpolate('spline', order=3)\n\n    return x_df.to_numpy(), v_df.to_numpy(), cov_x, cov_v","metadata":{"execution":{"iopub.status.busy":"2024-03-26T16:43:46.888004Z","iopub.execute_input":"2024-03-26T16:43:46.888792Z","iopub.status.idle":"2024-03-26T16:43:46.926295Z","shell.execute_reply.started":"2024-03-26T16:43:46.88874Z","shell.execute_reply":"2024-03-26T16:43:46.92474Z"},"trusted":true},"execution_count":7,"outputs":[]},{"cell_type":"code","source":"# Kalman filter\ndef Kalman_filter(zs, us, cov_zs, cov_us, phone):\n    # Parameters\n    sigma_mahalanobis = 30.0  # Mahalanobis distance for rejecting innovation\n\n    n, dim_x = zs.shape\n    F = np.eye(3)  # Transition matrix\n    H = np.eye(3)  # Measurement function\n\n    # Initial state and covariance\n    x = zs[0, :3].T  # State\n    P = 5.0**2 * np.eye(3)  # State covariance\n    I = np.eye(dim_x)\n\n    x_kf = np.zeros([n, dim_x])\n    P_kf = np.zeros([n, dim_x, dim_x])\n\n    # Kalman filtering\n    for i, (u, z) in enumerate(zip(us, zs)):\n        # First step\n        if i == 0:\n            x_kf[i] = x.T\n            P_kf[i] = P\n            continue\n\n        # Prediction step\n        Q = cov_us[i] # Estimated WLS velocity covariance\n        x = F @ x + u.T\n        P = (F @ P) @ F.T + Q\n\n        # Check outliers for observation\n        d = distance.mahalanobis(z, H @ x, np.linalg.inv(P))\n\n        # Update step\n        if d < sigma_mahalanobis:\n            R = cov_zs[i] # Estimated WLS position covariance\n            y = z.T - H @ x\n            S = (H @ P) @ H.T + R\n            K = (P @ H.T) @ np.linalg.inv(S)\n            x = x + K @ y\n            P = (I - (K @ H)) @ P\n        else:\n            # If observation update is not available, increase covariance\n            P += 10**2*Q\n\n        x_kf[i] = x.T\n        P_kf[i] = P\n\n    return x_kf, P_kf\n\n\n# Forward + backward Kalman filter and smoothing\ndef Kalman_smoothing(x_wls, v_wls, cov_x, cov_v, phone):\n    n, dim_x = x_wls.shape\n\n    # For some unknown reason, the speed estimation is wrong only for XiaomiMi8\n    # so the variance is increased\n    if phone == 'XiaomiMi8':\n        v_wls = np.vstack([(v_wls[:-1, :] + v_wls[1:, :])/2, np.zeros([1, 3])])\n        cov_v = 1000.0**2 * cov_v\n         \n    # Forward\n    v = np.vstack([np.zeros([1, 3]), (v_wls[:-1, :] + v_wls[1:, :])/2])\n    x_f, P_f = Kalman_filter(x_wls, v, cov_x, cov_v, phone)\n\n    # Backward\n    v = -np.flipud(v_wls)\n    v = np.vstack([np.zeros([1, 3]), (v[:-1, :] + v[1:, :])/2])\n    cov_xf = np.flip(cov_x, axis=0)\n    cov_vf = np.flip(cov_v, axis=0)\n    x_b, P_b = Kalman_filter(np.flipud(x_wls), v, cov_xf, cov_vf, phone)\n\n    # Smoothing\n    x_fb = np.zeros_like(x_f)\n    P_fb = np.zeros_like(P_f)\n    for (f, b) in zip(range(n), range(n-1, -1, -1)):\n        P_fi = np.linalg.inv(P_f[f])\n        P_bi = np.linalg.inv(P_b[b])\n\n        P_fb[f] = np.linalg.inv(P_fi + P_bi)\n        x_fb[f] = P_fb[f] @ (P_fi @ x_f[f] + P_bi @ x_b[b])\n\n    return x_fb, x_f, np.flipud(x_b)","metadata":{"execution":{"iopub.status.busy":"2024-03-26T16:44:03.731193Z","iopub.execute_input":"2024-03-26T16:44:03.732359Z","iopub.status.idle":"2024-03-26T16:44:03.752482Z","shell.execute_reply.started":"2024-03-26T16:44:03.732311Z","shell.execute_reply":"2024-03-26T16:44:03.751507Z"},"trusted":true},"execution_count":8,"outputs":[]},{"cell_type":"code","source":"!pip install nvector\n\nimport nvector as nv\n\ndef interpolate_timestamp(utcMillis, utc, data_points):\n    output = np.empty((0,2), float)\n    wgs84 = nv.FrameE(name='WGS84')\n    for x in utcMillis:\n        upper_bound = np.where(utc < x)\n        if len(upper_bound[0]) == 0:\n            print(\"Upper bound not there\")\n            output = np.append(output, np.array([data_points[0,:]]), axis=0)\n            continue\n        else:\n            upper_datapoint = data_points[upper_bound[0][-1]]\n            \n        lower_bound = np.where(utc > x)\n        if len(lower_bound[0]) == 0:\n            print(\"Lower bound not there\")\n            output = np.append(output, np.array([data_points[len(data_points)-1,:]]), axis=0)\n            continue\n        else:\n            lower_datapoint = data_points[lower_bound[0][0]]\n            \n        equal_timestamp = np.where(utc == x)\n        if len(equal_timestamp[0]) != 0:\n            output = np.append(output, np.array([data_points[equal_timestamp[0][0]]]), axis=0)\n            continue\n            \n        n_EB_E_t0 = wgs84.GeoPoint(upper_datapoint[0], upper_datapoint[0], degrees=True).to_nvector()\n        n_EB_E_t1 = wgs84.GeoPoint(lower_datapoint[0], lower_datapoint[1], degrees=True).to_nvector()\n        \n        path = nv.GeoPath(n_EB_E_t0, n_EB_E_t1)\n        g_EB_E_ti = path.interpolate(x).to_geo_point()\n        \n        lat_ti, lon_ti = g_EB_E_ti.latitude_deg, g_EB_E_ti.longitude_deg\n        output = np.append(output, np.array([lat_ti, lon_ti]), axis=0)\n        \n    return output","metadata":{"execution":{"iopub.status.busy":"2024-03-26T16:44:22.671688Z","iopub.execute_input":"2024-03-26T16:44:22.67211Z","iopub.status.idle":"2024-03-26T16:44:37.730663Z","shell.execute_reply.started":"2024-03-26T16:44:22.67208Z","shell.execute_reply":"2024-03-26T16:44:37.72898Z"},"trusted":true},"execution_count":9,"outputs":[{"name":"stdout","text":"Requirement already satisfied: nvector in /opt/conda/lib/python3.10/site-packages (0.7.7)\nRequirement already satisfied: numpy>=1.9 in /opt/conda/lib/python3.10/site-packages (from nvector) (1.26.4)\nRequirement already satisfied: scipy>=0.9 in /opt/conda/lib/python3.10/site-packages (from nvector) (1.11.4)\nRequirement already satisfied: geographiclib in /opt/conda/lib/python3.10/site-packages (from nvector) (2.0)\nRequirement already satisfied: matplotlib in /opt/conda/lib/python3.10/site-packages (from nvector) (3.7.5)\nRequirement already satisfied: cartopy in /opt/conda/lib/python3.10/site-packages (from nvector) (0.22.0)\nRequirement already satisfied: shapely>=1.7 in /opt/conda/lib/python3.10/site-packages (from cartopy->nvector) (1.8.5.post1)\nRequirement already satisfied: packaging>=20 in /opt/conda/lib/python3.10/site-packages (from cartopy->nvector) (21.3)\nRequirement already satisfied: pyshp>=2.1 in /opt/conda/lib/python3.10/site-packages (from cartopy->nvector) (2.3.1)\nRequirement already satisfied: pyproj>=3.1.0 in /opt/conda/lib/python3.10/site-packages (from cartopy->nvector) (3.6.1)\nRequirement already satisfied: contourpy>=1.0.1 in /opt/conda/lib/python3.10/site-packages (from matplotlib->nvector) (1.2.0)\nRequirement already satisfied: cycler>=0.10 in /opt/conda/lib/python3.10/site-packages (from matplotlib->nvector) (0.12.1)\nRequirement already satisfied: fonttools>=4.22.0 in /opt/conda/lib/python3.10/site-packages (from matplotlib->nvector) (4.47.0)\nRequirement already satisfied: kiwisolver>=1.0.1 in /opt/conda/lib/python3.10/site-packages (from matplotlib->nvector) (1.4.5)\nRequirement already satisfied: pillow>=6.2.0 in /opt/conda/lib/python3.10/site-packages (from matplotlib->nvector) (9.5.0)\nRequirement already satisfied: pyparsing>=2.3.1 in /opt/conda/lib/python3.10/site-packages (from matplotlib->nvector) (3.1.1)\nRequirement already satisfied: python-dateutil>=2.7 in /opt/conda/lib/python3.10/site-packages (from matplotlib->nvector) (2.9.0.post0)\nRequirement already satisfied: certifi in /opt/conda/lib/python3.10/site-packages (from pyproj>=3.1.0->cartopy->nvector) (2024.2.2)\nRequirement already satisfied: six>=1.5 in /opt/conda/lib/python3.10/site-packages (from python-dateutil>=2.7->matplotlib->nvector) (1.16.0)\n","output_type":"stream"}]},{"cell_type":"code","source":"path = '/kaggle/input/smartphone-decimeter-2023/sdc2023'\n\nsample_df = pd.read_csv(f'{path}/sample_submission.csv')\ntest_dfs = []\n\n# Loop for each trip\nfor i, dirname in enumerate(tqdm(sorted(gl.glob(f'{path}/test/*/*/')))):\n    drive, phone = dirname.split('/')[-3:-1]\n    tripID = f'{drive}/{phone}'\n    print(tripID)\n\n    # Read data\n    gnss_df = pd.read_csv(f'{dirname}/device_gnss.csv')\n\n    # Point positioning\n    utc, x_wls, v_wls, cov_x, cov_v = point_positioning(gnss_df)\n\n    # Exclude velocity outliers\n    x_wls, v_wls, cov_x, cov_v = exclude_interpolate_outlier(x_wls, v_wls, cov_x, cov_v)\n\n    # Kalman smoothing\n    x_kf, _, _ = Kalman_smoothing(x_wls, v_wls, cov_x, cov_v, phone)\n    assert np.all(~np.isnan(x_kf))\n\n    # Convert to latitude and longitude\n    llh_kf = np.array(pm.ecef2geodetic(x_kf[:, 0], x_kf[:, 1], x_kf[:, 2])).T\n\n    # Interpolation for submission\n    UnixTimeMillis = sample_df[sample_df['tripId'] == tripID]['UnixTimeMillis'].to_numpy()\n    lat_lng = interpolate_timestamp(UnixTimeMillis, utc, llh_kf[:,:2])\n#     lat = InterpolatedUnivariateSpline(utc, llh_kf[:,0], ext=3)(UnixTimeMillis)\n#     lng = InterpolatedUnivariateSpline(utc, llh_kf[:,1], ext=3)(UnixTimeMillis)\n    trip_df = pd.DataFrame({\n        'tripId' : tripID,\n        'UnixTimeMillis': UnixTimeMillis,\n        'LatitudeDegrees': lat_lng[:,0],\n        'LongitudeDegrees':  lat_lng[:,1]\n        })\n\n    test_dfs.append(trip_df)\n\n# Write submission.csv\ntest_df = pd.concat(test_dfs)\ntest_df.to_csv('submission.csv', index=False)","metadata":{"execution":{"iopub.status.busy":"2024-03-26T16:45:32.391796Z","iopub.execute_input":"2024-03-26T16:45:32.392637Z","iopub.status.idle":"2024-03-26T17:35:40.219737Z","shell.execute_reply.started":"2024-03-26T16:45:32.392588Z","shell.execute_reply":"2024-03-26T17:35:40.218314Z"},"trusted":true},"execution_count":11,"outputs":[{"output_type":"display_data","data":{"text/plain":"  0%|          | 0/40 [00:00<?, ?it/s]","application/vnd.jupyter.widget-view+json":{"version_major":2,"version_minor":0,"model_id":"46822ca207ea4dc789bc93fbedc6b7ef"}},"metadata":{}},{"name":"stdout","text":"2020-12-11-19-30-us-ca-mtv-e/pixel4xl\n","output_type":"stream"},{"output_type":"display_data","data":{"text/plain":"  0%|          | 0/1190 [00:00<?, ?it/s]","application/vnd.jupyter.widget-view+json":{"version_major":2,"version_minor":0,"model_id":"8096fcad309a4744ba710dd70bfee198"}},"metadata":{}},{"name":"stdout","text":"Number of velocity outliers 1\nNumber of position outliers 0\nUpper bound not there\nLower bound not there\n2021-08-17-20-37-us-ca-mtv-g/pixel5\n","output_type":"stream"},{"output_type":"display_data","data":{"text/plain":"  0%|          | 0/1676 [00:00<?, ?it/s]","application/vnd.jupyter.widget-view+json":{"version_major":2,"version_minor":0,"model_id":"c1b21eca63fa44d6a6e1dfc31264a168"}},"metadata":{}},{"name":"stdout","text":"Number of velocity outliers 0\nNumber of position outliers 0\nLower bound not there\n2021-08-31-20-37-us-ca-mtv-e/sm-g988b\n","output_type":"stream"},{"name":"stderr","text":"/tmp/ipykernel_33/3205060769.py:13: DtypeWarning: Columns (21) have mixed types. Specify dtype option on import or set low_memory=False.\n  gnss_df = pd.read_csv(f'{dirname}/device_gnss.csv')\n","output_type":"stream"},{"output_type":"display_data","data":{"text/plain":"  0%|          | 0/1141 [00:00<?, ?it/s]","application/vnd.jupyter.widget-view+json":{"version_major":2,"version_minor":0,"model_id":"2a3a4c72920142c3bdc6681ae07e01d1"}},"metadata":{}},{"name":"stdout","text":"Number of velocity outliers 9\nNumber of position outliers 0\nUpper bound not there\nLower bound not there\n2021-09-14-20-32-us-ca-mtv-k/pixel4\n","output_type":"stream"},{"output_type":"display_data","data":{"text/plain":"  0%|          | 0/1274 [00:00<?, ?it/s]","application/vnd.jupyter.widget-view+json":{"version_major":2,"version_minor":0,"model_id":"321074785afa4fb4b8ce4e8e64fe45dd"}},"metadata":{}},{"name":"stdout","text":"Number of velocity outliers 1\nNumber of position outliers 0\nUpper bound not there\nLower bound not there\n2021-09-20-19-03-us-ca-mtv-l/pixel4\n","output_type":"stream"},{"output_type":"display_data","data":{"text/plain":"  0%|          | 0/1797 [00:00<?, ?it/s]","application/vnd.jupyter.widget-view+json":{"version_major":2,"version_minor":0,"model_id":"eede9a376b5b404e800b1d5cffbd5759"}},"metadata":{}},{"name":"stdout","text":"Number of velocity outliers 0\nNumber of position outliers 0\nLower bound not there\n2021-09-28-21-56-us-ca-mtv-a/mi8\n","output_type":"stream"},{"name":"stderr","text":"/tmp/ipykernel_33/3205060769.py:13: DtypeWarning: Columns (21) have mixed types. Specify dtype option on import or set low_memory=False.\n  gnss_df = pd.read_csv(f'{dirname}/device_gnss.csv')\n","output_type":"stream"},{"output_type":"display_data","data":{"text/plain":"  0%|          | 0/2479 [00:00<?, ?it/s]","application/vnd.jupyter.widget-view+json":{"version_major":2,"version_minor":0,"model_id":"628bb109366d489581e6db4a1f5ff7d3"}},"metadata":{}},{"name":"stdout","text":"i = 2358 position lsq status = 2\ni = 2453 position lsq status = 2\nNumber of velocity outliers 10\nNumber of position outliers 2\nUpper bound not there\nLower bound not there\n2021-11-05-18-28-us-ca-mtv-m/pixel6pro\n","output_type":"stream"},{"output_type":"display_data","data":{"text/plain":"  0%|          | 0/1446 [00:00<?, ?it/s]","application/vnd.jupyter.widget-view+json":{"version_major":2,"version_minor":0,"model_id":"8958277768df48abb235b48a412285b5"}},"metadata":{}},{"name":"stdout","text":"Number of velocity outliers 0\nNumber of position outliers 0\nUpper bound not there\n2021-11-30-20-59-us-ca-mtv-m/mi8\n","output_type":"stream"},{"name":"stderr","text":"/tmp/ipykernel_33/3205060769.py:13: DtypeWarning: Columns (21) have mixed types. Specify dtype option on import or set low_memory=False.\n  gnss_df = pd.read_csv(f'{dirname}/device_gnss.csv')\n","output_type":"stream"},{"output_type":"display_data","data":{"text/plain":"  0%|          | 0/1395 [00:00<?, ?it/s]","application/vnd.jupyter.widget-view+json":{"version_major":2,"version_minor":0,"model_id":"a00e8ad97f53484fa570eb7197acbce6"}},"metadata":{}},{"name":"stdout","text":"i = 15 position lsq status = 0\ni = 101 position lsq status = 2\ni = 118 position lsq status = 2\ni = 120 position lsq status = 2\ni = 121 position lsq status = 2\ni = 123 position lsq status = 2\ni = 124 position lsq status = 2\ni = 125 position lsq status = 2\ni = 126 position lsq status = 2\ni = 127 position lsq status = 2\ni = 128 position lsq status = 2\ni = 129 position lsq status = 2\ni = 130 position lsq status = 2\ni = 131 position lsq status = 2\ni = 133 position lsq status = 2\ni = 135 position lsq status = 2\ni = 137 position lsq status = 2\ni = 139 position lsq status = 2\ni = 140 position lsq status = 2\ni = 141 position lsq status = 2\ni = 142 position lsq status = 2\ni = 143 position lsq status = 2\nNumber of velocity outliers 2\nNumber of position outliers 22\nUpper bound not there\nLower bound not there\n2022-02-08-22-04-us-ca-sjc-r/pixel5\n","output_type":"stream"},{"output_type":"display_data","data":{"text/plain":"  0%|          | 0/1665 [00:00<?, ?it/s]","application/vnd.jupyter.widget-view+json":{"version_major":2,"version_minor":0,"model_id":"74bfbe9791a84bb7a61e90d034cd0fb5"}},"metadata":{}},{"name":"stdout","text":"Number of velocity outliers 2\nNumber of position outliers 0\nLower bound not there\n2022-02-23-17-46-us-ca-lax-n/pixel5\n","output_type":"stream"},{"output_type":"display_data","data":{"text/plain":"  0%|          | 0/2407 [00:00<?, ?it/s]","application/vnd.jupyter.widget-view+json":{"version_major":2,"version_minor":0,"model_id":"cb806c36e6ce4b1a81c77ab80a58f4a9"}},"metadata":{}},{"name":"stdout","text":"Number of velocity outliers 0\nNumber of position outliers 0\nLower bound not there\n2022-02-23-22-35-us-ca-lax-m/pixel5\n","output_type":"stream"},{"output_type":"display_data","data":{"text/plain":"  0%|          | 0/2805 [00:00<?, ?it/s]","application/vnd.jupyter.widget-view+json":{"version_major":2,"version_minor":0,"model_id":"5e378dccd34e4ff9bf914ae57dddc554"}},"metadata":{}},{"name":"stdout","text":"Number of velocity outliers 7\nNumber of position outliers 5\nUpper bound not there\nLower bound not there\n2022-02-24-15-10-us-ca-lax-p/pixel5\n","output_type":"stream"},{"output_type":"display_data","data":{"text/plain":"  0%|          | 0/4515 [00:00<?, ?it/s]","application/vnd.jupyter.widget-view+json":{"version_major":2,"version_minor":0,"model_id":"e5ff5f5ba16048f3aa13c302c2aeffc5"}},"metadata":{}},{"name":"stdout","text":"i = 2752 position lsq status = 0\ni = 2752 velocity lsq status = 0\nNumber of velocity outliers 10\nNumber of position outliers 5\nLower bound not there\n2022-02-24-22-14-us-ca-lax-i/pixel5\n","output_type":"stream"},{"output_type":"display_data","data":{"text/plain":"  0%|          | 0/3582 [00:00<?, ?it/s]","application/vnd.jupyter.widget-view+json":{"version_major":2,"version_minor":0,"model_id":"d23f302def3a47d5815cc0ebdf3a3eba"}},"metadata":{}},{"name":"stdout","text":"Number of velocity outliers 4\nNumber of position outliers 0\nLower bound not there\n2022-03-14-19-57-us-ca-mtv-pe1/xiaomimi8\n","output_type":"stream"},{"name":"stderr","text":"/tmp/ipykernel_33/3205060769.py:13: DtypeWarning: Columns (21) have mixed types. Specify dtype option on import or set low_memory=False.\n  gnss_df = pd.read_csv(f'{dirname}/device_gnss.csv')\n","output_type":"stream"},{"output_type":"display_data","data":{"text/plain":"  0%|          | 0/1679 [00:00<?, ?it/s]","application/vnd.jupyter.widget-view+json":{"version_major":2,"version_minor":0,"model_id":"7de83244b87e490c982012941a3409db"}},"metadata":{}},{"name":"stdout","text":"Number of velocity outliers 5\nNumber of position outliers 0\nUpper bound not there\n2022-03-17-20-16-us-ca-sjc-q/sm-g988b\n","output_type":"stream"},{"output_type":"display_data","data":{"text/plain":"  0%|          | 0/1172 [00:00<?, ?it/s]","application/vnd.jupyter.widget-view+json":{"version_major":2,"version_minor":0,"model_id":"55019a6132f64fe18e587a0c271df651"}},"metadata":{}},{"name":"stdout","text":"Number of velocity outliers 1\nNumber of position outliers 0\nUpper bound not there\n2022-03-22-18-44-us-ca-mtv-pe1/pixel5\n","output_type":"stream"},{"output_type":"display_data","data":{"text/plain":"  0%|          | 0/2112 [00:00<?, ?it/s]","application/vnd.jupyter.widget-view+json":{"version_major":2,"version_minor":0,"model_id":"5d0bb4e477f849d3bd2fcec81550888f"}},"metadata":{}},{"name":"stdout","text":"Number of velocity outliers 4\nNumber of position outliers 3\nLower bound not there\n2022-04-04-16-31-us-ca-lax-x/pixel5\n","output_type":"stream"},{"output_type":"display_data","data":{"text/plain":"  0%|          | 0/2171 [00:00<?, ?it/s]","application/vnd.jupyter.widget-view+json":{"version_major":2,"version_minor":0,"model_id":"f9392a1b388f4c4f9ce49cf4f87f0aa7"}},"metadata":{}},{"name":"stdout","text":"Number of velocity outliers 11\nNumber of position outliers 0\nLower bound not there\n2022-04-22-20-11-us-ca-ebf-y/pixel5\n","output_type":"stream"},{"output_type":"display_data","data":{"text/plain":"  0%|          | 0/1400 [00:00<?, ?it/s]","application/vnd.jupyter.widget-view+json":{"version_major":2,"version_minor":0,"model_id":"057121fb76e04e90b42c4e885e44a993"}},"metadata":{}},{"name":"stdout","text":"Number of velocity outliers 8\nNumber of position outliers 2\nLower bound not there\n2022-04-25-21-04-us-ca-ebf-x/mi8\n","output_type":"stream"},{"name":"stderr","text":"/tmp/ipykernel_33/3205060769.py:13: DtypeWarning: Columns (21) have mixed types. Specify dtype option on import or set low_memory=False.\n  gnss_df = pd.read_csv(f'{dirname}/device_gnss.csv')\n","output_type":"stream"},{"output_type":"display_data","data":{"text/plain":"  0%|          | 0/1924 [00:00<?, ?it/s]","application/vnd.jupyter.widget-view+json":{"version_major":2,"version_minor":0,"model_id":"73eb35ab8a6e498bafd2b5abaa9cf9ee"}},"metadata":{}},{"name":"stdout","text":"i = 10 position lsq status = 2\ni = 700 position lsq status = 2\nNumber of velocity outliers 3\nNumber of position outliers 2\nUpper bound not there\n2022-04-25-22-36-us-ca-ebf-z/pixel5\n","output_type":"stream"},{"output_type":"display_data","data":{"text/plain":"  0%|          | 0/1587 [00:00<?, ?it/s]","application/vnd.jupyter.widget-view+json":{"version_major":2,"version_minor":0,"model_id":"f9274d0969084542bf450cc9b5f89e68"}},"metadata":{}},{"name":"stdout","text":"Number of velocity outliers 3\nNumber of position outliers 2\nLower bound not there\n2022-04-27-18-16-us-ca-ebf-zz/pixel5\n","output_type":"stream"},{"output_type":"display_data","data":{"text/plain":"  0%|          | 0/1315 [00:00<?, ?it/s]","application/vnd.jupyter.widget-view+json":{"version_major":2,"version_minor":0,"model_id":"851783e6deef4b1f86070e019ea6e1a1"}},"metadata":{}},{"name":"stdout","text":"Number of velocity outliers 2\nNumber of position outliers 2\nUpper bound not there\nLower bound not there\n2022-04-27-19-23-us-ca-ebf-xx/pixel5\n","output_type":"stream"},{"output_type":"display_data","data":{"text/plain":"  0%|          | 0/1382 [00:00<?, ?it/s]","application/vnd.jupyter.widget-view+json":{"version_major":2,"version_minor":0,"model_id":"0cd74f46a228442ea7d3a5c8135466e2"}},"metadata":{}},{"name":"stdout","text":"Number of velocity outliers 4\nNumber of position outliers 1\nUpper bound not there\nLower bound not there\n2022-04-27-21-55-us-ca-ebf-ww/mi8\n","output_type":"stream"},{"name":"stderr","text":"/tmp/ipykernel_33/3205060769.py:13: DtypeWarning: Columns (21) have mixed types. Specify dtype option on import or set low_memory=False.\n  gnss_df = pd.read_csv(f'{dirname}/device_gnss.csv')\n","output_type":"stream"},{"output_type":"display_data","data":{"text/plain":"  0%|          | 0/1344 [00:00<?, ?it/s]","application/vnd.jupyter.widget-view+json":{"version_major":2,"version_minor":0,"model_id":"3285279654434385994c1a10d1078327"}},"metadata":{}},{"name":"stdout","text":"i = 324 position lsq status = 2\ni = 331 position lsq status = 2\ni = 393 position lsq status = 2\ni = 508 position lsq status = 2\nNumber of velocity outliers 0\nNumber of position outliers 11\nUpper bound not there\nLower bound not there\n2022-05-12-20-19-us-ca-mtv-pe1/samsunga325g\n","output_type":"stream"},{"output_type":"display_data","data":{"text/plain":"  0%|          | 0/1435 [00:00<?, ?it/s]","application/vnd.jupyter.widget-view+json":{"version_major":2,"version_minor":0,"model_id":"cb1ce69cd37b46b79c226ab9f3199362"}},"metadata":{}},{"name":"stdout","text":"Number of velocity outliers 8\nNumber of position outliers 2\nUpper bound not there\n2022-06-22-20-12-us-ca-lax-hh/samsunga325g\n","output_type":"stream"},{"output_type":"display_data","data":{"text/plain":"  0%|          | 0/1862 [00:00<?, ?it/s]","application/vnd.jupyter.widget-view+json":{"version_major":2,"version_minor":0,"model_id":"5f44043cf7a743138ca5f6ccb84650e7"}},"metadata":{}},{"name":"stdout","text":"Number of velocity outliers 12\nNumber of position outliers 2\nUpper bound not there\n2022-06-28-20-56-us-ca-sjc-r/samsunga32\n","output_type":"stream"},{"output_type":"display_data","data":{"text/plain":"  0%|          | 0/1838 [00:00<?, ?it/s]","application/vnd.jupyter.widget-view+json":{"version_major":2,"version_minor":0,"model_id":"f6428c4772b64094b3e4896a91ae181b"}},"metadata":{}},{"name":"stdout","text":"Number of velocity outliers 9\nNumber of position outliers 2\nUpper bound not there\n2022-07-12-18-37-us-ca-mtv-b/sm-a325f\n","output_type":"stream"},{"output_type":"display_data","data":{"text/plain":"  0%|          | 0/1782 [00:00<?, ?it/s]","application/vnd.jupyter.widget-view+json":{"version_major":2,"version_minor":0,"model_id":"c0965d68a9274f66a3f57a36948db1d4"}},"metadata":{}},{"name":"stdout","text":"Number of velocity outliers 3\nNumber of position outliers 0\nUpper bound not there\nLower bound not there\n2022-10-06-20-46-us-ca-sjc-r/sm-a205u\n","output_type":"stream"},{"output_type":"display_data","data":{"text/plain":"  0%|          | 0/1699 [00:00<?, ?it/s]","application/vnd.jupyter.widget-view+json":{"version_major":2,"version_minor":0,"model_id":"addf54884ea1440aa0c37cac1b3c77b5"}},"metadata":{}},{"name":"stdout","text":"Number of velocity outliers 9\nNumber of position outliers 2\nUpper bound not there\nLower bound not there\n2023-04-27-19-25-us-ca-mtv-pe1/pixel5\n","output_type":"stream"},{"output_type":"display_data","data":{"text/plain":"  0%|          | 0/1357 [00:00<?, ?it/s]","application/vnd.jupyter.widget-view+json":{"version_major":2,"version_minor":0,"model_id":"dfdd9a14337346e5910fb592e607a463"}},"metadata":{}},{"name":"stdout","text":"Number of velocity outliers 2\nNumber of position outliers 1\nLower bound not there\n2023-04-27-20-55-us-ca-sjc-q/pixel5\n","output_type":"stream"},{"output_type":"display_data","data":{"text/plain":"  0%|          | 0/1380 [00:00<?, ?it/s]","application/vnd.jupyter.widget-view+json":{"version_major":2,"version_minor":0,"model_id":"6e65591e0cff4a728d81d5ab13d0fb0d"}},"metadata":{}},{"name":"stdout","text":"Number of velocity outliers 0\nNumber of position outliers 0\nLower bound not there\n2023-05-02-19-24-us-ca-sjc-we1/pixel7pro\n","output_type":"stream"},{"output_type":"display_data","data":{"text/plain":"  0%|          | 0/2110 [00:00<?, ?it/s]","application/vnd.jupyter.widget-view+json":{"version_major":2,"version_minor":0,"model_id":"b1ff45857f964e64a83c8831e1af3785"}},"metadata":{}},{"name":"stdout","text":"Number of velocity outliers 2\nNumber of position outliers 0\nUpper bound not there\nLower bound not there\n2023-05-02-20-33-us-ca-mtv-xe1/pixel7pro\n","output_type":"stream"},{"name":"stderr","text":"/tmp/ipykernel_33/3205060769.py:13: DtypeWarning: Columns (21) have mixed types. Specify dtype option on import or set low_memory=False.\n  gnss_df = pd.read_csv(f'{dirname}/device_gnss.csv')\n","output_type":"stream"},{"output_type":"display_data","data":{"text/plain":"  0%|          | 0/2145 [00:00<?, ?it/s]","application/vnd.jupyter.widget-view+json":{"version_major":2,"version_minor":0,"model_id":"df782ce0705e4550b27a7154a095968f"}},"metadata":{}},{"name":"stdout","text":"Number of velocity outliers 1\nNumber of position outliers 0\nLower bound not there\n2023-05-09-23-10-us-ca-sjc-r/sm-a505u\n","output_type":"stream"},{"output_type":"display_data","data":{"text/plain":"  0%|          | 0/2384 [00:00<?, ?it/s]","application/vnd.jupyter.widget-view+json":{"version_major":2,"version_minor":0,"model_id":"996839af7c3e42ad88d97b9204f900f5"}},"metadata":{}},{"name":"stdout","text":"i = 1571 position lsq status = 2\ni = 1573 position lsq status = 2\nNumber of velocity outliers 5\nNumber of position outliers 2\nUpper bound not there\nLower bound not there\n2023-05-23-21-06-us-ca-mtv-de1/pixel5\n","output_type":"stream"},{"output_type":"display_data","data":{"text/plain":"  0%|          | 0/1975 [00:00<?, ?it/s]","application/vnd.jupyter.widget-view+json":{"version_major":2,"version_minor":0,"model_id":"a9b58db568e34bbc9963108bc2c310fe"}},"metadata":{}},{"name":"stdout","text":"Number of velocity outliers 2\nNumber of position outliers 0\nLower bound not there\n2023-05-23-22-16-us-ca-mtv-ie2/pixel6pro\n","output_type":"stream"},{"output_type":"display_data","data":{"text/plain":"  0%|          | 0/1020 [00:00<?, ?it/s]","application/vnd.jupyter.widget-view+json":{"version_major":2,"version_minor":0,"model_id":"44482717f8c94d47a4afa87903f7b2f4"}},"metadata":{}},{"name":"stdout","text":"Number of velocity outliers 0\nNumber of position outliers 0\nLower bound not there\n2023-05-25-17-32-us-ca-pao-j/pixel6pro\n","output_type":"stream"},{"name":"stderr","text":"/tmp/ipykernel_33/3205060769.py:13: DtypeWarning: Columns (21) have mixed types. Specify dtype option on import or set low_memory=False.\n  gnss_df = pd.read_csv(f'{dirname}/device_gnss.csv')\n","output_type":"stream"},{"output_type":"display_data","data":{"text/plain":"  0%|          | 0/1292 [00:00<?, ?it/s]","application/vnd.jupyter.widget-view+json":{"version_major":2,"version_minor":0,"model_id":"c64e832c25d348dfa8ba20654ac33aaf"}},"metadata":{}},{"name":"stdout","text":"Number of velocity outliers 0\nNumber of position outliers 0\nLower bound not there\n2023-05-25-21-50-us-ca-sjc-ke2/sm-s908b\n","output_type":"stream"},{"output_type":"display_data","data":{"text/plain":"  0%|          | 0/1728 [00:00<?, ?it/s]","application/vnd.jupyter.widget-view+json":{"version_major":2,"version_minor":0,"model_id":"4f85b972a9aa46298a0f7b9fcd0180e2"}},"metadata":{}},{"name":"stdout","text":"Number of velocity outliers 3\nNumber of position outliers 0\nUpper bound not there\nLower bound not there\n2023-05-26-21-23-us-ca-sjc-be2/pixel5\n","output_type":"stream"},{"output_type":"display_data","data":{"text/plain":"  0%|          | 0/1482 [00:00<?, ?it/s]","application/vnd.jupyter.widget-view+json":{"version_major":2,"version_minor":0,"model_id":"a5eeca3aae98403bb3ad5e20735b37c6"}},"metadata":{}},{"name":"stdout","text":"Number of velocity outliers 0\nNumber of position outliers 0\nUpper bound not there\nLower bound not there\n2023-06-06-22-43-us-ca-sjc-he2/pixel5\n","output_type":"stream"},{"output_type":"display_data","data":{"text/plain":"  0%|          | 0/1608 [00:00<?, ?it/s]","application/vnd.jupyter.widget-view+json":{"version_major":2,"version_minor":0,"model_id":"eb99900d635948edb26678e0aed19ccd"}},"metadata":{}},{"name":"stdout","text":"Number of velocity outliers 0\nNumber of position outliers 0\nUpper bound not there\nLower bound not there\n2023-06-15-18-49-us-ca-sjc-ce1/pixel7pro\n","output_type":"stream"},{"name":"stderr","text":"/tmp/ipykernel_33/3205060769.py:13: DtypeWarning: Columns (21) have mixed types. Specify dtype option on import or set low_memory=False.\n  gnss_df = pd.read_csv(f'{dirname}/device_gnss.csv')\n","output_type":"stream"},{"output_type":"display_data","data":{"text/plain":"  0%|          | 0/1495 [00:00<?, ?it/s]","application/vnd.jupyter.widget-view+json":{"version_major":2,"version_minor":0,"model_id":"c07c9248b1df4c3db2aba6549c7a0c75"}},"metadata":{}},{"name":"stdout","text":"i = 628 position lsq status = 2\ni = 647 position lsq status = 2\nNumber of velocity outliers 1\nNumber of position outliers 3\nUpper bound not there\nLower bound not there\n","output_type":"stream"}]},{"cell_type":"code","source":"import numpy as np # linear algebra\nimport pandas as pd # data processing, CSV file I/O (e.g. pd.read_csv)\nimport os\nimport random\nimport tensorflow as tf\nfrom sklearn.model_selection import train_test_split\nfrom sklearn.preprocessing import StandardScaler\n\n# data_date = []\n# base_dir = \"/kaggle/input/smartphone-decimeter-2022/train/\"\nscaler = StandardScaler()\nfor dirname, _, filenames in os.walk('/kaggle/input/smartphone-decimeter-2023/sdc2023/train'):\n    for filename in filenames:\n        if filename == 'device_gnss.csv':\n            data = pd.read_csv(os.path.join(dirname, filename))\n            preprocess_data = data[[\"SvClockBiasMeters\", \"IsrbMeters\", \"IonosphericDelayMeters\", \"TroposphericDelayMeters\", \"SvPositionXEcefMeters\", \"SvPositionYEcefMeters\", \"SvPositionZEcefMeters\", \"RawPseudorangeUncertaintyMeters\"]].to_numpy()\n            preprocess_data = preprocess_data[~np.isnan(preprocess_data).any(axis=1)]\n            train_data_len = int(len(preprocess_data) * 0.8)\n\n            train_data_select = scaler.fit_transform(preprocess_data[:train_data_len,:]) \n            test_data_select = scaler.fit_transform(preprocess_data[train_data_len: ,:])\n\n            min_val = tf.reduce_min(preprocess_data)\n            max_val = tf.reduce_max(preprocess_data)\n\n            train_data = (train_data_select - min_val) / (max_val - min_val)\n            test_data = (test_data_select - min_val) / (max_val - min_val)\n\n            train_data = tf.cast(train_data, tf.float32)\n            test_data = tf.cast(test_data, tf.float32)\n            \n            history = autoencoder.fit(train_data, train_data, \n                  epochs=20, \n                  batch_size=512,\n                  validation_data=(test_data, test_data),\n                  shuffle=False)\n#         dirname_list = dirname.split(\"/\")\n#         data_date.append((dirname_list)[len(dirname_list)-2])\n        \n# random.shuffle(data_date)\n# train_percent = int(len(data_date) * 0.8)\n# test_percent = int(len(data_date) * 0.2)\n\n# train_dir = data_date[:train_percent]\n# test_dir = data_date[train_percent:(train_percent + test_percent)]","metadata":{"execution":{"iopub.status.busy":"2024-03-26T17:35:51.273526Z","iopub.execute_input":"2024-03-26T17:35:51.274678Z","iopub.status.idle":"2024-03-26T17:36:08.70422Z","shell.execute_reply.started":"2024-03-26T17:35:51.27464Z","shell.execute_reply":"2024-03-26T17:36:08.702794Z"},"trusted":true},"execution_count":12,"outputs":[{"name":"stderr","text":"2024-03-26 17:35:53.863449: E external/local_xla/xla/stream_executor/cuda/cuda_dnn.cc:9261] Unable to register cuDNN factory: Attempting to register factory for plugin cuDNN when one has already been registered\n2024-03-26 17:35:53.863617: E external/local_xla/xla/stream_executor/cuda/cuda_fft.cc:607] Unable to register cuFFT factory: Attempting to register factory for plugin cuFFT when one has already been registered\n2024-03-26 17:35:54.063547: E external/local_xla/xla/stream_executor/cuda/cuda_blas.cc:1515] Unable to register cuBLAS factory: Attempting to register factory for plugin cuBLAS when one has already been registered\n","output_type":"stream"},{"traceback":["\u001b[0;31m---------------------------------------------------------------------------\u001b[0m","\u001b[0;31mNameError\u001b[0m                                 Traceback (most recent call last)","Cell \u001b[0;32mIn[12], line 32\u001b[0m\n\u001b[1;32m     29\u001b[0m             train_data \u001b[38;5;241m=\u001b[39m tf\u001b[38;5;241m.\u001b[39mcast(train_data, tf\u001b[38;5;241m.\u001b[39mfloat32)\n\u001b[1;32m     30\u001b[0m             test_data \u001b[38;5;241m=\u001b[39m tf\u001b[38;5;241m.\u001b[39mcast(test_data, tf\u001b[38;5;241m.\u001b[39mfloat32)\n\u001b[0;32m---> 32\u001b[0m             history \u001b[38;5;241m=\u001b[39m \u001b[43mautoencoder\u001b[49m\u001b[38;5;241m.\u001b[39mfit(train_data, train_data, \n\u001b[1;32m     33\u001b[0m                   epochs\u001b[38;5;241m=\u001b[39m\u001b[38;5;241m20\u001b[39m, \n\u001b[1;32m     34\u001b[0m                   batch_size\u001b[38;5;241m=\u001b[39m\u001b[38;5;241m512\u001b[39m,\n\u001b[1;32m     35\u001b[0m                   validation_data\u001b[38;5;241m=\u001b[39m(test_data, test_data),\n\u001b[1;32m     36\u001b[0m                   shuffle\u001b[38;5;241m=\u001b[39m\u001b[38;5;28;01mFalse\u001b[39;00m)\n\u001b[1;32m     37\u001b[0m \u001b[38;5;66;03m#         dirname_list = dirname.split(\"/\")\u001b[39;00m\n\u001b[1;32m     38\u001b[0m \u001b[38;5;66;03m#         data_date.append((dirname_list)[len(dirname_list)-2])\u001b[39;00m\n\u001b[1;32m     39\u001b[0m         \n\u001b[0;32m   (...)\u001b[0m\n\u001b[1;32m     44\u001b[0m \u001b[38;5;66;03m# train_dir = data_date[:train_percent]\u001b[39;00m\n\u001b[1;32m     45\u001b[0m \u001b[38;5;66;03m# test_dir = data_date[train_percent:(train_percent + test_percent)]\u001b[39;00m\n","\u001b[0;31mNameError\u001b[0m: name 'autoencoder' is not defined"],"ename":"NameError","evalue":"name 'autoencoder' is not defined","output_type":"error"}]},{"cell_type":"code","source":"import keras\nfrom keras import layers\nfrom tensorflow.keras.models import Model\nimport tensorflow as tf\nfrom tensorflow.keras.models import Model\n\n#Train data attributes for position estimation are Satellite position, velocity, RawPseudorangeUncertaintyMeters\n#pr_smooth, SvClockBiasMeters, IsrbMeters, IonosphericDelayMeters, TroposphericDelayMeters\nencoding_dim = 3\n\nclass AnomalyDetector(Model):\n  def __init__(self):\n    super(AnomalyDetector, self).__init__()\n    self.encoder = tf.keras.Sequential([\n      layers.Dense(32, activation=\"relu\"),\n      layers.Dense(16, activation=\"relu\"),\n      layers.Dense(3, activation=\"relu\")])\n\n    self.decoder = tf.keras.Sequential([\n      layers.Dense(16, activation=\"relu\"),\n      layers.Dense(32, activation=\"relu\"),\n      layers.Dense(8, activation=\"sigmoid\")])\n\n  def call(self, x):\n    encoded = self.encoder(x)\n    decoded = self.decoder(encoded)\n    return decoded\n\nautoencoder = AnomalyDetector()\nautoencoder.compile(optimizer='adam', loss=tf.keras.losses.MeanSquaredError())","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"","metadata":{},"execution_count":null,"outputs":[]}]}