{"metadata":{"kernelspec":{"language":"python","display_name":"Python 3","name":"python3"},"language_info":{"pygments_lexer":"ipython3","nbconvert_exporter":"python","version":"3.6.4","file_extension":".py","codemirror_mode":{"name":"ipython","version":3},"name":"python","mimetype":"text/x-python"}},"nbformat_minor":4,"nbformat":4,"cells":[{"cell_type":"code","source":"import numpy as np\nimport pandas as pd\nimport datetime\nfrom sklearn import preprocessing\nimport os\nfrom datetime import datetime\nimport matplotlib.pyplot as plt\nimport math\n\npd.options.display.max_columns = 3000\npd.options.display.max_rows = 3000","metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"acc_df = pd.read_csv('../input/features-imu-preprocessing-clean/cleaned/imu/imu_acc_train.csv')\nmag_df = pd.read_csv('../input/features-imu-preprocessing-clean/cleaned/imu/imu_magnet_train.csv')\n# wypt_df = pd.read_csv('../input/data-pre-processing-to-csv/cleaned/waypoint_train.csv')\n# wifi_df = pd.read_csv('../input/data-pre-processing-to-csv/cleaned/wifi_train.csv')","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Accelerometer - calculate velocity \nreferences: https://www.kaggle.com/iamleonie/intro-to-indoor-location-navigation\n\nhttps://www.kaggle.com/museas/estimate-the-walk-with-acce-and-magn/comments\n\naccelerometer measures change in velocity \n\nunit: m/s^2 \n\nmean of z = gravity \n\nEquation: https://www.calculatorsoup.com/calculators/physics/displacement_v_a_t.php\n\n1. convert timestamp to seconds: ms/1000\n2. calculte time difference delta\n3. convert time delta to seconds\n4. velocity = acceleration * timedelta in seconds\n5. position = p_0 + velocity * timedelta in seconds\n\n* calculate position and velocity from acceleration = calculate acceleration and velocity from positon\n\n* But since wypt and accelerometer data won't match (has lots of NaN and only few intersections), integrating over the acceleration will also integrate the error for each sample which can quickly accumulate the cause large deviations. \n\n* Thus, can use velocity as a feature in pre-processing (as to calculate velocity we don't need position)","metadata":{}},{"cell_type":"code","source":"def calc_from_acce(timestamp, acce, p_0):\n    df = pd.DataFrame({'timestamp' : timestamp, 'acceleration' : acce})\n    df['timestamp_ms'] = df['timestamp'].apply(lambda x: datetime.fromtimestamp(x/1000.0))\n    df['timedelta_ms'] = df['timestamp_ms'].diff()\n    df['timedelta_s'] = df['timedelta_ms'].apply(lambda x: x.total_seconds()).fillna(0)\n    df['velocity'] = (df['acceleration']*df['timedelta_s']).cumsum()\n    df['position'] = p_0 + (df['velocity']*df['timedelta_s']).cumsum()\n    return df[['timestamp', 'timestamp_ms', 'timedelta_s', 'position', 'velocity', 'acceleration']]","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# referenced from https://www.kaggle.com/iamleonie/intro-to-indoor-location-navigation\n# and made some edits \ndef calc_velocity(timestamp, acce):\n    df = pd.DataFrame({'timestamp' : timestamp, 'acceleration' : acce})\n    df['timestamp_ms'] = df['timestamp'].apply(lambda x: datetime.fromtimestamp(x/1000.0))\n    df['timedelta_ms'] = df['timestamp_ms'].diff()\n    df['timedelta_s'] = df['timedelta_ms'].apply(lambda x: x.total_seconds()).fillna(0)\n    df['velocity'] = (df['acceleration']*df['timedelta_s'])\n    #df['position'] = p_0 + (df['velocity']*df['timedelta_s']).cumsum()\n    return df[['timestamp', 'timestamp_ms', 'timedelta_s', 'velocity', 'acceleration']]","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"acc_df = acc_df.sort_values(['building_id', 'floor', 'path_id', 'Time'])\nacc_df = acc_df[['building_id','floor','path_id', 'Time','acce_x','acce_y','acce_z','acce_accuracy']]\nacc_df.head()","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Try sample accelerometer\nacc_samp = acc_df[(acc_df['building_id'] == '5c3c44b80379370013e0fd2b') & \n                  (acc_df['floor'] == -1) & \n                  (acc_df['path_id']=='5d0761b54cae4f000a2db567')]\n\nacc_sampe_v = calc_velocity(acc_samp['Time'],(-1)*acc_samp['acce_x'])\nacc_sampe_v.drop_duplicates()","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Accelerometer - infer Number of steps\n\nThe idea of infer # of steps from accelerometer was from this notebook\nhttps://www.kaggle.com/museas/estimate-the-walk-with-acce-and-magn/comments\n, which the notebook refers from a github focuses on research on step detection algorithm https://github.com/danielmurray/adaptiv\n\n\nGeneral idea:\n1. using accelerometer define crest and trough\n2. calculate the distance between crest to trough \n3. define a treshold = how large the distence would you define a single step \n\n\n\nLimitation: as mentioend in the github, the crossings varies person by person, thus hard to define a treshold \n\n\ndefault treshold: crossing mean * 1.1\n\n\n\nDefinitions used: \n1. mixed all three accelerometers: x, y, z \n2. used butterworth highpass filter to eliminate low frequency components and generate a smaller (less noised) mixed accelerometer \n3. calculate crossing\n4. determine number of steps \n","metadata":{}},{"cell_type":"code","source":"def mix_acce(timestamp,x,y,z):\n    \"\"\"This function mixes the three accelerometers\"\"\"\n    mix_acce = np.sqrt(x**2 + y**2 + z**2)\n#     mix_acce = np.concatenate(timestamp, mix_acce, 1)\n    mix_acce = pd.DataFrame(pd.np.column_stack([timestamp,mix_acce]))\n    mix_df = pd.DataFrame(mix_acce)\n    mix_df.columns = [\"timestamp\",\"acce\"]\n    return mix_df","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"mix_acce = mix_acce(acc_samp['Time'],acc_samp['acce_x'],acc_samp['acce_y'],acc_samp['acce_z'])\nmix_acce","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"from scipy.signal import butter, lfilter\n\ndef butter_lowpass(cutoff, fs, order=5):\n    nyq = 0.5 * fs\n    normal_cutoff = cutoff / nyq\n    b, a = butter(order, normal_cutoff, btype='low', analog=False)\n    return b, a\n\ndef butter_lowpass_filter(data, cutoff, fs, order=5):\n    b, a = butter_lowpass(cutoff, fs, order=order)\n    y = lfilter(b, a, data)\n    return y\n\norder = 3\nfs = 50.0  # sample rate, Hz\n# fs = 100\n# cutoff = 3.667  # desired cutoff frequency of the filter, Hz\ncutoff = 3\n\n# filtered = butter_lowpass_filter(mix_df[\"acce\"], cutoff, fs, order)\n# filtered","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"acc_filter = butter_lowpass_filter(mix_acce['acce'],cutoff,fs,order)\nacc_filter","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# referenced from https://www.kaggle.com/museas/estimate-the-walk-with-acce-and-magn/comments\n# and made some edits \ndef peak_accel_threshold(data, timestamps, threshold):\n    \"\"\" This function used to calculate the crossings (difference between crest and trough) based on certain threshold\"\"\"\n    d_acc = []\n    last_state = 'below'\n    # below - less than threshold\n    # above - above the threshold\n    crest_troughs = 0\n    crossings = []\n\n    for i, datum in enumerate(data):\n        \n        current_state = last_state\n        if datum < threshold:\n            current_state = 'below'\n        elif datum > threshold:\n            current_state = 'above'\n\n        if current_state is not last_state:\n            if current_state is 'above':\n                crossing = [timestamps[i], threshold]\n                crossings.append(crossing)\n            else:\n                crossing = [timestamps[i], threshold]\n                crossings.append(crossing)\n\n            crest_troughs += 1\n        last_state = current_state\n\n    return np.array(crossings)","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"\nthreshold = acc_filter.mean() * 1.1\nprint(threshold)\ncrossings = peak_accel_threshold(acc_filter, mix_acce[\"timestamp\"], threshold)\n\nfig, ax = plt.subplots(nrows=1, ncols=1, figsize=(20, 5))\nplt.plot(mix_acce[\"timestamp\"], acc_filter, 'b-', color='cornflowerblue')\nplt.plot(crossings.T[0], crossings.T[1], 'ro', linewidth=2)\nplt.xlabel('Time [sec]')\nplt.grid()\nplt.legend()\nplt.show()\n\nprint(\"sum of crossings: \", len(crossings))\nprint(\"infer number of steps: \", len(crossings)/5/10)\n ","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# try another path ","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Using Magneic field to infer direction (calculating the euler angles) use as a compass\n\nnotebook referenced: https://www.kaggle.com/museas/estimate-the-walk-with-acce-and-magn/comments\n\nEuler Angle Equation: https://www.researchgate.net/figure/The-Euler-angles-Roll-Pitch-and-Yaw-Prior-to-magnetic-compensation-the-recorded_fig1_275771788\n\nEuler Angle Definition: Euler Angle describes the orientation of a rigid body with respect to a fixed coordinate system. It can also represents the orientation of a mobile frame of reference in physics or the orientation of a general basis in 3-dimensional linear algebra\n\nDefinition of roll, pitch: https://howthingsfly.si.edu/flight-dynamics/roll-pitch-and-yaw\n\nReferences: https://math.stackexchange.com/questions/2466949/get-magnetic-field-values-from-euler-angle\n\nGeneral idea: \n\n1. Use acceleromer to calculate roll: the theta angel between ay, az\n2. Use accelerometer to calculate pitch\n3. Define m_trans (a threshold we have to define )\n4. calcule euler angel based on equation ","metadata":{}},{"cell_type":"code","source":"mag_df = mag_df[['building_id','floor','path_id','Time','magnet_x','magnet_y','magnet_z']]\nacc_df = acc_df[['building_id','floor','path_id','Time','acce_x','acce_y','acce_z']]","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"mag_samp = mag_df[mag_df['path_id']=='5dc7f1951cda370006031a78']","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"acc_samp = acc_df[acc_df['path_id'] == '5dc7f1951cda370006031a78']","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"del mag_df, acc_df","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"print(acc_samp.shape)\nprint(mag_samp.shape)","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# even just one path subset fails.... \nacc_mag = pd.merge(mag_samp,acc_samp,how='inner',\n                   on=['building_id', 'floor', 'path_id', 'Time']).sort_values(['building_id', 'floor', 'path_id', 'Time'])","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"acc_mag.shape","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Tried using outer but the notebook crasshed. \n# used inner crashed again.... still crashed \n# tried using just one path ... still crashed...\n# Find out path that are both existed in magneticfield and accelerometer --- to pick samples \n# import itertools\n# path_acc = set(acc_df['path_id'])\n# path_mag = set(mag_df['path_id'])\n# path = path_acc.intersection(path_mag)\n#path\n# '5dce74bb94e49000061250ca' --- failed\n#  '5dc7f1951cda370006031a78',\n#  '5da6e6b0ab074400064325d4',\n#  '5dada49baa1d300006faaa08',\n#  '5d074cbab53a8d0008dd4928',\n#  '5dc7887317ffdd0006f11ed5',\n#  '5da9aae7df065a00069be7f5',\n#  '5da833e44091ff000604890d',\n#  '5dd23c4494e4900006126ae4',\n#  '5da54fbb6b51bc0006438eb8',\n#  '5dd10a5594e4900006126256',\n#  '5dc520fd21dceb00061147ed',","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"acc_mag","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"test = acc_mag.drop_duplicates()\nprint(test.shape)","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"acc_mag.iloc[1][4]","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Edited from https://www.kaggle.com/museas/estimate-the-walk-with-acce-and-magn/comments\ndef calc_direction(mag,acc):\n    \"\"\"This Function calculates the direction the surveyor is facing at each timestamp using euler angle \"\"\"\n    m_trans = -5 #0\n    time_di_list = []\n    df = pd.merge(mag,acc,how='inner',\n                   on=['building_id', 'floor', 'path_id', 'Time']).sort_values(['building_id', 'floor', 'path_id', 'Time'])\n    print(df.shape)\n    for i in df.iterrows():\n        gx,gy,gz = i[1][4],i[1][5],i[1][6]\n        ax,ay,az = i[1][7],i[1][8],i[1][9]\n\n        roll = math.atan2(ay,az)\n        pitch = math.atan2(-1*ax , (ay * math.sin(roll) + az * math.cos(roll)))\n\n        q =  m_trans - math.degrees(math.atan2(\n                                (gz*math.sin(roll)-gy*math.cos(roll)),(gx*math.cos(pitch) + gy*math.sin(roll)*math.sin(pitch) \n                                                                       + gz*math.sin(pitch)*math.cos(roll))\n                                )) -90\n        if q <= 0:\n            q += 360\n        \n        time_di_list.append((i[1][3],q))\n    d_list = [x[1] for x in time_di_list]\n    print(len(d_list))\n        \n    df['Direction'] = d_list\n    return df","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"direction_samp = calc_direction(mag_samp,acc_samp)\ndirection_samp","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"dir_nodup = direction_samp.drop_duplicates(subset=['Time'])\ndir_nodup","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Visualize \n# Copied from Github Helper function: https://github.com/location-competition/indoor-location-competition-20/blob/master/visualize_f.py\ndef visualize_trajectory(trajectory, floor_plan_filename, width_meter, height_meter, title=None, mode='lines + markers + text', show=False):\n    \"\"\"\n    Copied from from https://github.com/location-competition/indoor-location-competition-20/blob/master/visualize_f.py\n\n    \"\"\"\n    fig = go.Figure()\n\n    # add trajectory\n    size_list = [6] * trajectory.shape[0]\n    size_list[0] = 10\n    size_list[-1] = 10\n\n    color_list = ['rgba(4, 174, 4, 0.5)'] * trajectory.shape[0]\n    color_list[0] = 'rgba(12, 5, 235, 1)'\n    color_list[-1] = 'rgba(235, 5, 5, 1)'\n\n    position_count = {}\n    text_list = []\n    for i in range(trajectory.shape[0]):\n        if str(trajectory[i]) in position_count:\n            position_count[str(trajectory[i])] += 1\n        else:\n            position_count[str(trajectory[i])] = 0\n        text_list.append('        ' * position_count[str(trajectory[i])] + f'{i}')\n    text_list[0] = 'Start 0'\n    text_list[-1] = f'End {trajectory.shape[0] - 1}'\n\n    fig.add_trace(go.Scatter(x=pred_px, y=pred_py, name=\"pred\"))\n    \n    fig.add_trace(\n        go.Scattergl(\n            x=trajectory[:, 0],\n            y=trajectory[:, 1],\n            mode=mode,\n            marker=dict(size=size_list, color=color_list),\n            line=dict(shape='linear', color='lightgrey', width=3, dash='dash'),\n            text=text_list,\n            textposition=\"top center\",\n            name='trajectory',\n        ))\n\n    # add floor plan\n    floor_plan = Image.open(floor_plan_filename)\n    fig.update_layout(images=[\n        go.layout.Image(\n            source=floor_plan,\n            xref=\"x\",\n            yref=\"y\",\n            x=0,\n            y=height_meter,\n            sizex=width_meter,\n            sizey=height_meter,\n            sizing=\"contain\",\n            opacity=1,\n            layer=\"below\",\n        )\n    ])\n\n    # configure\n    fig.update_xaxes(autorange=False, range=[0, width_meter])\n    fig.update_yaxes(autorange=False, range=[0, height_meter], scaleanchor=\"x\", scaleratio=1)\n    fig.update_layout(\n        title=go.layout.Title(\n            text=title or \"No title.\",\n            xref=\"paper\",\n            x=0,\n        ),\n        autosize=True,\n        width=800,\n        height=  800 * height_meter / width_meter,\n        template=\"plotly_white\",\n    )\n\n    if show:\n        fig.show()\n\n    return fig\n\ndef visualize_train_trajectory(path):\n    \"\"\"\n    Edited from \n    https://www.kaggle.com/ihelon/indoor-location-exploratory-data-analysis\n    \"\"\"\n    _id, floor = path.split(\"/\")[:2]\n    \n    train_floor_data = read_data_file(f\"../input/indoor-location-navigation/train/{path}\")\n    with open(f\"../input/indoor-location-navigation/metadata/{_id}/{floor}/floor_info.json\") as f:\n        train_floor_info = json.load(f)\n\n    return visualize_trajectory(\n        train_floor_data.waypoint[:, 1:3], \n        f\"../input/indoor-location-navigation/metadata/{_id}/{floor}/floor_image.png\",\n        train_floor_info[\"map_info\"][\"width\"], \n        train_floor_info[\"map_info\"][\"height\"],\n        f\"Visualization of {path}\"\n    )\n#visualize_train_trajectory(a)","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"build_id = '5d2709e003f801723c32d896'\nfloor = 'F4'\npath_id = '5dc7f1951cda370006031a78'\np = build_id + '/' + floor + '/' + path_id + '.txt'\np","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"!git clone https://github.com/location-competition/indoor-location-competition-20.git","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"%cd /kaggle/working/indoor-location-competition-20/","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"from io_f import read_data_file","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# visualize_train_trajectory('5cd56b77e2acfd2d33b5b22b/F5/5d09fa831a9498000952a1a4.txt')","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"","metadata":{}}]}