{"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":"# 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'):\n    for filename in filenames:\n        print(os.path.join(dirname, filename))\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":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","_kg_hide-output":true,"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Prepare paths:\nimport glob\nfrom pathlib import Path\ninpath = '../input/indoor-location-navigation/'\nmetapath = inpath + 'metadata/'\ntrainpath = inpath + 'train/'\ntestpath = inpath + 'test/'\n\n# Extract testing files, buildings and sites:\nos.system(f'grep SiteID {testpath}/* > test_buildings.txt' )\ntest_buildings = pd.read_csv('test_buildings.txt',sep='\\t',header=None,names=['file','building','site'])\ntest_buildings['file'] = test_buildings['file'].apply(lambda x: x[:-2])\ntest_buildings['building'] = test_buildings['building'].apply(lambda x: x[7:])\n\n# How many buildings in the testing set?\nbuildings = np.unique(test_buildings['building'])\nprint('There are',len(buildings),'buildings in the testing set.')\n\ntest_buildings.head()","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Floor dictionary as defined by Devin in his nice work @ https://www.kaggle.com/devinanzelmo/wifi-features","metadata":{}},{"cell_type":"code","source":"floors = {\"B3\":-3,\"B2\":-2,\"B1\":-1,\"F1\":0,\"1F\":0,\"F2\":1,\"2F\":1,\"F3\":2,\"3F\":2,\"F4\":3,\"4F\":3,\n          \"F5\":4,\"5F\":4,\"F6\":5,\"6F\":5,\"F7\":6,\"7F\":6,\"F8\":7,\"8F\": 7,\"F9\":8,\"9F\":8,\"F10\":9}\n\n# Other floors:\nfloorsB = {\"B\":0,\"BF\":1,\"BM\":2}\nfloorsG = {\"G\":0}\nfloorsM = {\"M\":0}\nfloorsP = {\"P1\":0,\"P2\":1}\nfloorsL = {\"LG2\":-2,\"LG1\":-1,\"LG\":0,\"LM\":0,\"L1\":1,\"L2\":2,\"L3\":3,\"L4\":4,\"L5\":5,\"L6\":6,\"L7\":7,\"L8\":8,\"L9\":9,\"L10\":10,\"L11\":11}","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"A custom C++ code I wrote to fix the data quality issues reported by the hosts and to reformat the dataset without losing any information:","metadata":{}},{"cell_type":"code","source":"# Compile C++ pre-processing code:\ner=os.system(\"g++ /kaggle/input/indoor-cpp/1_preprocess.cpp -std=c++11 -o preprocess\")\nif(er): print(\"Error\")\n\n# Fix data quality issues and reformat dataset:\nos.system('mkdir train')\nfor building in buildings:\n    for floor in os.listdir(trainpath+building):\n        path_filenames = list(Path(trainpath+building+f'/{floor}/').resolve().glob(\"*.txt\"))\n        for path_filename in path_filenames:\n            if floor in floors:\n                er=os.system(f'./preprocess {path_filename} train {building} {floors[floor]}')\n            if(er): print(\"Error:\",path_filename)","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Visualization example:","metadata":{}},{"cell_type":"code","source":"# Column names in the C++ ordering:\ncols = ['t_start','xyz_time','x_acce','y_acce','z_acce','a_acce',\n        'x_acce_uncali','y_acce_uncali','z_acce_uncali',\n        'x2_acce_uncali','y2_acce_uncali','z2_acce_uncali','a_acce_uncali',\n        'x_gyro','y_gyro','z_gyro','a_gyro',\n        'x_gyro_uncali','y_gyro_uncali','z_gyro_uncali',\n        'x2_gyro_uncali','y2_gyro_uncali','z2_gyro_uncali','a_gyro_uncali',\n        'x_magn','y_magn','z_magn','a_magn',\n        'x_magn_uncali','y_magn_uncali','z_magn_uncali',\n        'x2_magn_uncali','y2_magn_uncali','z2_magn_uncali','a_magn_uncali',\n        'x_ahrs','y_ahrs','z_ahrs','a_ahrs',\n        't1_wifi','ssid_wifi','bssid_wifi','rssid_wifi','freq_wifi','t2_wifi',\n        't1_beac','id_beac','power_beac','rssi_beac','dist_beac','mac_beac','t2_beac',\n        't_waypoint','floor','x_waypoint','y_waypoint']\n\n# Column names by data type:\nint8col = ['floor','a_acce','a_gyro','a_magn','a_magn_uncali','a_ahrs',\n           'a_acce_uncali','a_gyro_uncali','a_magn_uncali']\nint16col = ['rssid_wifi','freq_wifi','power_beac','rssi_beac']\nint32col = ['t_start','xyz_time','t1_wifi','t2_wifi','t1_beac','t2_beac','t_waypoint']\nfloatcol = ['x_acce','y_acce','z_acce','x_acce_uncali','y_acce_uncali','z_acce_uncali',\n            'x2_acce_uncali','y2_acce_uncali','z2_acce_uncali',\n            'x_gyro','y_gyro','z_gyro','x_gyro_uncali','y_gyro_uncali','z_gyro_uncali',\n            'x2_gyro_uncali','y2_gyro_uncali','z2_gyro_uncali',\n            'x_magn','y_magn','z_magn','x_magn_uncali','y_magn_uncali','z_magn_uncali',\n            'x2_magn_uncali','y2_magn_uncali','z2_magn_uncali',\n            'x_ahrs','y_ahrs','z_ahrs','dist_beac','x_waypoint','y_waypoint']\nstringcol = ['ssid_wifi','bssid_wifi','id_beac','mac_beac']\ntarcol = ['t_waypoint','floor','x_waypoint','y_waypoint']\n\n# Choose any building:\nbuilding = '5a0546857ecc773753327266'\nxtrain = pd.read_csv(f'train/{building}.txt',index_col=0,header=None,names=cols,dtype='object')\nxtrain = xtrain.fillna('')\nxtrain[int8col] = xtrain[int8col].apply(lambda x: x.apply(lambda y: np.array(y.split(),dtype=np.int8)))\nxtrain[int16col] = xtrain[int16col].apply(lambda x: x.apply(lambda y: np.array(y.split(),dtype=np.int16)))\nxtrain[int32col] = xtrain[int32col].apply(lambda x: x.apply(lambda y: np.array(y.split(),dtype=np.int32)))\nxtrain[floatcol] = xtrain[floatcol].apply(lambda x: x.apply(lambda y: np.array(y.split(),dtype=np.float32)))\nxtrain[stringcol] = xtrain[stringcol].apply(lambda x: x.apply(lambda y: np.array(y.split())))\nxtrain.head()","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"import matplotlib.pyplot as plt\n\n# Choose any path:\npathid = '5e157315a280850006f3d00d'\n# XYZ:\nplt.figure(figsize=(25,5))\nplt.subplot(141)\n# https://developer.android.com/guide/topics/sensors/sensors_overview#sensors-coords\nfor col in ['x_acce','y_acce','z_acce','a_acce']: #y_acce is the acceleration in the walking direction, z_acce is in the vertical direction (i.e. including gravity) ?\n    plt.plot(xtrain.loc[pathid,'xyz_time'],xtrain.loc[pathid,col],label=col)\n    plt.legend(loc=\"upper right\")\n    plt.xlabel(\"xyz_time\")\n    plt.title('Acce')\nplt.subplot(142)\nfor col in ['x_gyro','y_gyro','z_gyro','a_gyro']: #z_gyro is the angular velocity around the axis pointing to the sky (i.e. the change in walking direction) ?\n    plt.plot(xtrain.loc[pathid,'xyz_time'],xtrain.loc[pathid,col],label=col)\n    plt.legend(loc=\"upper right\")\n    plt.xlabel(\"xyz_time\")\n    plt.title('Gyro')\nplt.subplot(143)\nfor col in ['x_magn','y_magn','z_magn','a_magn']: #maybe need to use magnitude? (x_magn and y_magn will depend on the orientation, z_magn should be independent if the phone is held correctly) ?\n    plt.plot(xtrain.loc[pathid,'xyz_time'],xtrain.loc[pathid,col],label=col)\n    plt.legend(loc=\"upper right\")\n    plt.xlabel(\"xyz_time\")\n    plt.title('Magn')\nplt.subplot(144)\n# https://developer.android.com/guide/topics/sensors/sensors_motion#sensors-motion-rotate\nfor col in ['x_ahrs','y_ahrs','z_ahrs','a_ahrs']: #z_ahrs is the rotation around the axis pointing to the sky (i.e. the walking direction) ?\n    plt.plot(xtrain.loc[pathid,'xyz_time'],xtrain.loc[pathid,col],label=col)\n    plt.legend(loc=\"upper right\")\n    plt.xlabel(\"xyz_time\")\n    plt.title('Ahrs')\nplt.show()\n#XYZ Uncali:\nplt.figure(figsize=(25,5))\nplt.subplot(141)\nfor col in ['x_acce_uncali','y_acce_uncali','z_acce_uncali','a_acce_uncali']:\n    plt.plot(xtrain.loc[pathid,'xyz_time'],xtrain.loc[pathid,col],label=col)\n    plt.legend(loc=\"upper right\")\n    plt.xlabel(\"xyz_time\")\n    plt.title('Acce_uncali')\nplt.subplot(142)\nfor col in ['x_gyro_uncali','y_gyro_uncali','z_gyro_uncali','a_gyro_uncali']:\n    plt.plot(xtrain.loc[pathid,'xyz_time'],xtrain.loc[pathid,col],label=col)\n    plt.legend(loc=\"upper right\")\n    plt.xlabel(\"xyz_time\")\n    plt.title('Gyro_uncali')\nplt.subplot(143)\nfor col in ['x_magn_uncali','y_magn_uncali','z_magn_uncali','a_magn_uncali']:\n    plt.plot(xtrain.loc[pathid,'xyz_time'],xtrain.loc[pathid,col],label=col)\n    plt.legend(loc=\"upper right\")\n    plt.xlabel(\"xyz_time\")\n    plt.title('Magn_uncali')\nplt.subplot(144)\nfor col in ['x_waypoint','y_waypoint']:\n    plt.plot(xtrain.loc[pathid,'t_waypoint'],xtrain.loc[pathid,col],'-+',label=col)\n    plt.legend(loc=\"right\")\n    plt.xlabel(\"t_waypoint\")\n    plt.title('Waypoint')\nplt.show()\n#XYZ Uncali 2:\nplt.figure(figsize=(25,5))\nplt.subplot(131)\nfor col in ['x2_acce_uncali','y2_acce_uncali','z2_acce_uncali']:\n    plt.plot(xtrain.loc[pathid,'xyz_time'],xtrain.loc[pathid,col],label=col)\n    plt.legend(loc=\"upper right\")\n    plt.xlabel(\"xyz_time\")\n    plt.title('Acce_uncali_2')\nplt.subplot(132)\nfor col in ['x2_gyro_uncali','y2_gyro_uncali','z2_gyro_uncali']:\n    plt.plot(xtrain.loc[pathid,'xyz_time'],xtrain.loc[pathid,col],label=col)\n    plt.legend(loc=\"upper right\")\n    plt.xlabel(\"xyz_time\")\n    plt.title('Gyro_uncali_2')\nplt.subplot(133)\nfor col in ['x2_magn_uncali','y2_magn_uncali','z2_magn_uncali']:\n    plt.plot(xtrain.loc[pathid,'xyz_time'],xtrain.loc[pathid,col],label=col)\n    plt.legend(loc=\"upper right\")\n    plt.xlabel(\"xyz_time\")\n    plt.title('Magn_uncali_2')\nplt.show()\n# Wifi:\nplt.figure(figsize=(25,5))\nax = plt.subplot(131)\npd.DataFrame({'rssid_wifi':xtrain.loc[pathid,'rssid_wifi'],'bssid_wifi':xtrain.loc[pathid,'bssid_wifi']},\n    index=xtrain.loc[pathid,'t1_wifi']).groupby('bssid_wifi')['rssid_wifi'].plot(ax=ax,title='rssid_wifi',legend=False)\nax = plt.subplot(132)\npd.DataFrame({'freq_wifi':xtrain.loc[pathid,'freq_wifi'],'bssid_wifi':xtrain.loc[pathid,'bssid_wifi']},\n    index=xtrain.loc[pathid,'t1_wifi']).groupby('bssid_wifi')['freq_wifi'].plot(ax=ax,title='freq_wifi',legend=False)\nax = plt.subplot(133)\npd.DataFrame({'t2_wifi':xtrain.loc[pathid,'t2_wifi'],'bssid_wifi':xtrain.loc[pathid,'bssid_wifi']},\n    index=xtrain.loc[pathid,'t1_wifi']).groupby('bssid_wifi')['t2_wifi'].plot(ax=ax,title='t2_wifi',legend=False)\nplt.show()\n# Beacon. Note: use mac to split ids, but use ids to train model (i.e. multiple sequences for each id, one for each mac) ?\nplt.figure(figsize=(25,5))\nax = plt.subplot(131)\npd.DataFrame({'rssi_beac':xtrain.loc[pathid,'rssi_beac'],'mac_beac':xtrain.loc[pathid,'mac_beac']},\n    index=xtrain.loc[pathid,'t1_beac']).groupby('mac_beac')['rssi_beac'].plot(ax=ax,title='rssi_beac',legend=False)\nax = plt.subplot(132)\npd.DataFrame({'dist_beac':xtrain.loc[pathid,'dist_beac'],'mac_beac':xtrain.loc[pathid,'mac_beac']},\n    index=xtrain.loc[pathid,'t1_beac']).groupby('mac_beac')['dist_beac'].plot(ax=ax,title='dist_beac',legend=False)\nax = plt.subplot(133)\npd.DataFrame({'power_beac':xtrain.loc[pathid,'power_beac'],'id_beac':xtrain.loc[pathid,'id_beac']},\n    index=xtrain.loc[pathid,'t1_beac']).groupby('id_beac')['power_beac'].plot(ax=ax,title='power_beac',legend=False)\nplt.show()\nplt.figure(figsize=(25,5))\nax = plt.subplot(131)\npd.DataFrame({'rssi_beac':xtrain.loc[pathid,'rssi_beac'],'id_beac':xtrain.loc[pathid,'id_beac']},\n    index=xtrain.loc[pathid,'t1_beac']).groupby('id_beac')['rssi_beac'].plot(ax=ax,title='rssi_beac',legend=False)\nax = plt.subplot(132)\npd.DataFrame({'dist_beac':xtrain.loc[pathid,'dist_beac'],'id_beac':xtrain.loc[pathid,'id_beac']},\n    index=xtrain.loc[pathid,'t1_beac']).groupby('id_beac')['dist_beac'].plot(ax=ax,title='dist_beac',legend=False)\nax = plt.subplot(133)\npd.DataFrame({'power_beac':xtrain.loc[pathid,'power_beac'],'id_beac':xtrain.loc[pathid,'id_beac']},\n    index=xtrain.loc[pathid,'t1_beac']).groupby('id_beac')['power_beac'].plot(ax=ax,title='power_beac',legend=False)\nplt.show()","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Run the GitHub scripts shared by the hosts on the training set @ https://github.com/location-competition/indoor-location-competition-20 Copyright (c) 2017-2020 XYZ10, Inc. https://dangwu.com/","metadata":{}},{"cell_type":"code","source":"import numpy as np\nimport scipy.signal as signal\n\n\ndef split_ts_seq(ts_seq, sep_ts):\n    \"\"\"\n\n    :param ts_seq:\n    :param sep_ts:\n    :return:\n    \"\"\"\n    tss = ts_seq[:, 0].astype(float)\n    unique_sep_ts = np.unique(sep_ts)\n    ts_seqs = []\n    start_index = 0\n    for i in range(0, unique_sep_ts.shape[0]):\n        end_index = np.searchsorted(tss, unique_sep_ts[i], side='right')\n        if start_index == end_index:\n            continue\n        ts_seqs.append(ts_seq[start_index:end_index, :].copy())\n        start_index = end_index\n\n    # tail data\n    if start_index < ts_seq.shape[0]:\n        ts_seqs.append(ts_seq[start_index:, :].copy())\n\n    return ts_seqs\n\n\ndef correct_trajectory(original_xys, end_xy):\n    \"\"\"\n\n    :param original_xys: numpy ndarray, shape(N, 2)\n    :param end_xy: numpy ndarray, shape(1, 2)\n    :return:\n    \"\"\"\n    corrected_xys = np.zeros((0, 2))\n\n    A = original_xys[0, :]\n    B = end_xy\n    Bp = original_xys[-1, :]\n\n    angle_BAX = np.arctan2(B[1] - A[1], B[0] - A[0])\n    angle_BpAX = np.arctan2(Bp[1] - A[1], Bp[0] - A[0])\n    angle_BpAB = angle_BpAX - angle_BAX\n    AB = np.sqrt(np.sum((B - A) ** 2))\n    ABp = np.sqrt(np.sum((Bp - A) ** 2))\n\n    corrected_xys = np.append(corrected_xys, [A], 0)\n    for i in np.arange(1, np.size(original_xys, 0)):\n        angle_CpAX = np.arctan2(original_xys[i, 1] - A[1], original_xys[i, 0] - A[0])\n\n        angle_CAX = angle_CpAX - angle_BpAB\n\n        ACp = np.sqrt(np.sum((original_xys[i, :] - A) ** 2))\n\n        AC = ACp * AB / ABp\n\n        delta_C = np.array([AC * np.cos(angle_CAX), AC * np.sin(angle_CAX)])\n\n        C = delta_C + A\n\n        corrected_xys = np.append(corrected_xys, [C], 0)\n\n    return corrected_xys\n\n\ndef correct_positions(rel_positions, reference_positions):\n    \"\"\"\n\n    :param rel_positions:\n    :param reference_positions:\n    :return:\n    \"\"\"\n    rel_positions_list = split_ts_seq(rel_positions, reference_positions[:, 0])\n    if len(rel_positions_list) != reference_positions.shape[0] - 1:\n        # print(f'Rel positions list size: {len(rel_positions_list)}, ref positions size: {reference_positions.shape[0]}')\n        del rel_positions_list[-1]\n    assert len(rel_positions_list) == reference_positions.shape[0] - 1\n\n    corrected_positions = np.zeros((0, 3))\n    for i, rel_ps in enumerate(rel_positions_list):\n        start_position = reference_positions[i]\n        end_position = reference_positions[i + 1]\n        abs_ps = np.zeros(rel_ps.shape)\n        abs_ps[:, 0] = rel_ps[:, 0]\n        # abs_ps[:, 1:3] = rel_ps[:, 1:3] + start_position[1:3]\n        abs_ps[0, 1:3] = rel_ps[0, 1:3] + start_position[1:3]\n        for j in range(1, rel_ps.shape[0]):\n            abs_ps[j, 1:3] = abs_ps[j-1, 1:3] + rel_ps[j, 1:3]\n        abs_ps = np.insert(abs_ps, 0, start_position, axis=0)\n        corrected_xys = correct_trajectory(abs_ps[:, 1:3], end_position[1:3])\n        corrected_ps = np.column_stack((abs_ps[:, 0], corrected_xys))\n        if i == 0:\n            corrected_positions = np.append(corrected_positions, corrected_ps, axis=0)\n        else:\n            corrected_positions = np.append(corrected_positions, corrected_ps[1:], axis=0)\n\n    corrected_positions = np.array(corrected_positions)\n\n    return corrected_positions\n\n\ndef init_parameters_filter(sample_freq, warmup_data, cut_off_freq=2):\n    order = 4\n    filter_b, filter_a = signal.butter(order, cut_off_freq / (sample_freq / 2), 'low', False)\n    zf = signal.lfilter_zi(filter_b, filter_a)\n    _, zf = signal.lfilter(filter_b, filter_a, warmup_data, zi=zf)\n    _, filter_zf = signal.lfilter(filter_b, filter_a, warmup_data, zi=zf)\n\n    return filter_b, filter_a, filter_zf\n\n\ndef get_rotation_matrix_from_vector(rotation_vector):\n    q1 = rotation_vector[0]\n    q2 = rotation_vector[1]\n    q3 = rotation_vector[2]\n\n    if rotation_vector.size >= 4:\n        q0 = rotation_vector[3]\n    else:\n        q0 = 1 - q1*q1 - q2*q2 - q3*q3\n        if q0 > 0:\n            q0 = np.sqrt(q0)\n        else:\n            q0 = 0\n\n    sq_q1 = 2 * q1 * q1\n    sq_q2 = 2 * q2 * q2\n    sq_q3 = 2 * q3 * q3\n    q1_q2 = 2 * q1 * q2\n    q3_q0 = 2 * q3 * q0\n    q1_q3 = 2 * q1 * q3\n    q2_q0 = 2 * q2 * q0\n    q2_q3 = 2 * q2 * q3\n    q1_q0 = 2 * q1 * q0\n\n    R = np.zeros((9,))\n    if R.size == 9:\n        R[0] = 1 - sq_q2 - sq_q3\n        R[1] = q1_q2 - q3_q0\n        R[2] = q1_q3 + q2_q0\n\n        R[3] = q1_q2 + q3_q0\n        R[4] = 1 - sq_q1 - sq_q3\n        R[5] = q2_q3 - q1_q0\n\n        R[6] = q1_q3 - q2_q0\n        R[7] = q2_q3 + q1_q0\n        R[8] = 1 - sq_q1 - sq_q2\n\n        R = np.reshape(R, (3, 3))\n    elif R.size == 16:\n        R[0] = 1 - sq_q2 - sq_q3\n        R[1] = q1_q2 - q3_q0\n        R[2] = q1_q3 + q2_q0\n        R[3] = 0.0\n\n        R[4] = q1_q2 + q3_q0\n        R[5] = 1 - sq_q1 - sq_q3\n        R[6] = q2_q3 - q1_q0\n        R[7] = 0.0\n\n        R[8] = q1_q3 - q2_q0\n        R[9] = q2_q3 + q1_q0\n        R[10] = 1 - sq_q1 - sq_q2\n        R[11] = 0.0\n\n        R[12] = R[13] = R[14] = 0.0\n        R[15] = 1.0\n\n        R = np.reshape(R, (4, 4))\n\n    return R\n\n\ndef get_orientation(R):\n    flat_R = R.flatten()\n    values = np.zeros((3,))\n    if np.size(flat_R) == 9:\n        values[0] = np.arctan2(flat_R[1], flat_R[4])\n        values[1] = np.arcsin(-flat_R[7])\n        values[2] = np.arctan2(-flat_R[6], flat_R[8])\n    else:\n        values[0] = np.arctan2(flat_R[1], flat_R[5])\n        values[1] = np.arcsin(-flat_R[9])\n        values[2] = np.arctan2(-flat_R[8], flat_R[10])\n\n    return values\n\n\ndef compute_steps(acce_datas):\n    step_timestamps = np.array([])\n    step_indexs = np.array([], dtype=int)\n    step_acce_max_mins = np.zeros((0, 4))\n    sample_freq = 50\n    window_size = 22\n    low_acce_mag = 0.6\n    step_criterion = 1\n    interval_threshold = 250\n\n    acce_max = np.zeros((2,))\n    acce_min = np.zeros((2,))\n    acce_binarys = np.zeros((window_size,), dtype=int)\n    acce_mag_pre = 0\n    state_flag = 0\n\n    warmup_data = np.ones((window_size,)) * 9.81\n    filter_b, filter_a, filter_zf = init_parameters_filter(sample_freq, warmup_data)\n    acce_mag_window = np.zeros((window_size, 1))\n\n    # detect steps according to acceleration magnitudes\n    for i in np.arange(0, np.size(acce_datas, 0)):\n        acce_data = acce_datas[i, :]\n        acce_mag = np.sqrt(np.sum(acce_data[1:] ** 2))\n\n        acce_mag_filt, filter_zf = signal.lfilter(filter_b, filter_a, [acce_mag], zi=filter_zf)\n        acce_mag_filt = acce_mag_filt[0]\n\n        acce_mag_window = np.append(acce_mag_window, [acce_mag_filt])\n        acce_mag_window = np.delete(acce_mag_window, 0)\n        mean_gravity = np.mean(acce_mag_window)\n        acce_std = np.std(acce_mag_window)\n        mag_threshold = np.max([low_acce_mag, 0.4 * acce_std])\n\n        # detect valid peak or valley of acceleration magnitudes\n        acce_mag_filt_detrend = acce_mag_filt - mean_gravity\n        if acce_mag_filt_detrend > np.max([acce_mag_pre, mag_threshold]):\n            # peak\n            acce_binarys = np.append(acce_binarys, [1])\n            acce_binarys = np.delete(acce_binarys, 0)\n        elif acce_mag_filt_detrend < np.min([acce_mag_pre, -mag_threshold]):\n            # valley\n            acce_binarys = np.append(acce_binarys, [-1])\n            acce_binarys = np.delete(acce_binarys, 0)\n        else:\n            # between peak and valley\n            acce_binarys = np.append(acce_binarys, [0])\n            acce_binarys = np.delete(acce_binarys, 0)\n\n        if (acce_binarys[-1] == 0) and (acce_binarys[-2] == 1):\n            if state_flag == 0:\n                acce_max[:] = acce_data[0], acce_mag_filt\n                state_flag = 1\n            elif (state_flag == 1) and ((acce_data[0] - acce_max[0]) <= interval_threshold) and (\n                    acce_mag_filt > acce_max[1]):\n                acce_max[:] = acce_data[0], acce_mag_filt\n            elif (state_flag == 2) and ((acce_data[0] - acce_max[0]) > interval_threshold):\n                acce_max[:] = acce_data[0], acce_mag_filt\n                state_flag = 1\n\n        # choose reasonable step criterion and check if there is a valid step\n        # save step acceleration data: step_acce_max_mins = [timestamp, max, min, variance]\n        step_flag = False\n        if step_criterion == 2:\n            if (acce_binarys[-1] == -1) and ((acce_binarys[-2] == 1) or (acce_binarys[-2] == 0)):\n                step_flag = True\n        elif step_criterion == 3:\n            if (acce_binarys[-1] == -1) and (acce_binarys[-2] == 0) and (np.sum(acce_binarys[:-2]) > 1):\n                step_flag = True\n        else:\n            if (acce_binarys[-1] == 0) and acce_binarys[-2] == -1:\n                if (state_flag == 1) and ((acce_data[0] - acce_min[0]) > interval_threshold):\n                    acce_min[:] = acce_data[0], acce_mag_filt\n                    state_flag = 2\n                    step_flag = True\n                elif (state_flag == 2) and ((acce_data[0] - acce_min[0]) <= interval_threshold) and (\n                        acce_mag_filt < acce_min[1]):\n                    acce_min[:] = acce_data[0], acce_mag_filt\n        if step_flag:\n            step_timestamps = np.append(step_timestamps, acce_data[0])\n            step_indexs = np.append(step_indexs, [i])\n            step_acce_max_mins = np.append(step_acce_max_mins,\n                                           [[acce_data[0], acce_max[1], acce_min[1], acce_std ** 2]], axis=0)\n        acce_mag_pre = acce_mag_filt_detrend\n\n    return step_timestamps, step_indexs, step_acce_max_mins\n\n\ndef compute_stride_length(step_acce_max_mins):\n    K = 0.4\n    K_max = 0.8\n    K_min = 0.4\n    para_a0 = 0.21468084\n    para_a1 = 0.09154517\n    para_a2 = 0.02301998\n\n    stride_lengths = np.zeros((step_acce_max_mins.shape[0], 2))\n    k_real = np.zeros((step_acce_max_mins.shape[0], 2))\n    step_timeperiod = np.zeros((step_acce_max_mins.shape[0] - 1, ))\n    stride_lengths[:, 0] = step_acce_max_mins[:, 0]\n    window_size = 2\n    step_timeperiod_temp = np.zeros((0, ))\n\n    # calculate every step period - step_timeperiod unit: second\n    for i in range(0, step_timeperiod.shape[0]):\n        step_timeperiod_data = (step_acce_max_mins[i + 1, 0] - step_acce_max_mins[i, 0]) / 1000\n        step_timeperiod_temp = np.append(step_timeperiod_temp, [step_timeperiod_data])\n        if step_timeperiod_temp.shape[0] > window_size:\n            step_timeperiod_temp = np.delete(step_timeperiod_temp, [0])\n        step_timeperiod[i] = np.sum(step_timeperiod_temp) / step_timeperiod_temp.shape[0]\n\n    # calculate parameters by step period and acceleration magnitude variance\n    k_real[:, 0] = step_acce_max_mins[:, 0]\n    k_real[0, 1] = K\n    for i in range(0, step_timeperiod.shape[0]):\n        k_real[i + 1, 1] = np.max([(para_a0 + para_a1 / step_timeperiod[i] + para_a2 * step_acce_max_mins[i, 3]), K_min])\n        k_real[i + 1, 1] = np.min([k_real[i + 1, 1], K_max]) * (K / K_min)\n\n    # calculate every stride length by parameters and max and min data of acceleration magnitude\n    stride_lengths[:, 1] = np.max([(step_acce_max_mins[:, 1] - step_acce_max_mins[:, 2]),\n                                   np.ones((step_acce_max_mins.shape[0], ))], axis=0)**(1 / 4) * k_real[:, 1]\n\n    return stride_lengths\n\n\ndef compute_headings(ahrs_datas):\n    headings = np.zeros((np.size(ahrs_datas, 0), 2))\n    for i in np.arange(0, np.size(ahrs_datas, 0)):\n        ahrs_data = ahrs_datas[i, :]\n        rot_mat = get_rotation_matrix_from_vector(ahrs_data[1:])\n        azimuth, pitch, roll = get_orientation(rot_mat)\n        around_z = (-azimuth) % (2 * np.pi)\n        headings[i, :] = ahrs_data[0], around_z\n    return headings\n\n\ndef compute_step_heading(step_timestamps, headings):\n    step_headings = np.zeros((len(step_timestamps), 2))\n    step_timestamps_index = 0\n    for i in range(0, len(headings)):\n        if step_timestamps_index < len(step_timestamps):\n            if headings[i, 0] == step_timestamps[step_timestamps_index]:\n                step_headings[step_timestamps_index, :] = headings[i, :]\n                step_timestamps_index += 1\n        else:\n            break\n    assert step_timestamps_index == len(step_timestamps)\n\n    return step_headings\n\n\ndef compute_rel_positions(stride_lengths, step_headings):\n    rel_positions = np.zeros((stride_lengths.shape[0], 3))\n    for i in range(0, stride_lengths.shape[0]):\n        rel_positions[i, 0] = stride_lengths[i, 0]\n        rel_positions[i, 1] = -stride_lengths[i, 1] * np.sin(step_headings[i, 1])\n        rel_positions[i, 2] = stride_lengths[i, 1] * np.cos(step_headings[i, 1])\n\n    return rel_positions\n\n\ndef compute_step_positions(acce_datas, ahrs_datas, posi_datas):\n    step_timestamps, step_indexs, step_acce_max_mins = compute_steps(acce_datas)\n    headings = compute_headings(ahrs_datas)\n    stride_lengths = compute_stride_length(step_acce_max_mins)\n    step_headings = compute_step_heading(step_timestamps, headings)\n    rel_positions = compute_rel_positions(stride_lengths, step_headings)\n    step_positions = correct_positions(rel_positions, posi_datas)\n\n    return step_positions\n\ndef extract_magnetic_strength(mwi_datas):\n    magnetic_strength = {}\n    for position_key in mwi_datas:\n        # print(f'Position: {position_key}')\n        magnetic_data = mwi_datas[position_key]['magnetic']\n        magnetic_s = np.mean(np.sqrt(np.sum(magnetic_data[:, 1:4] ** 2, axis=1)))\n        magnetic_strength[position_key] = magnetic_s\n    return magnetic_strength\ndef extract_wifi_rssi(mwi_datas):\n    wifi_rssi = {}\n    for position_key in mwi_datas:\n        # print(f'Position: {position_key}')\n        wifi_data = mwi_datas[position_key]['wifi']\n        for wifi_d in wifi_data:\n            bssid = wifi_d[2]\n            rssi = int(wifi_d[3])\n            if bssid in wifi_rssi:\n                position_rssi = wifi_rssi[bssid]\n                if position_key in position_rssi:\n                    old_rssi = position_rssi[position_key][0]\n                    old_count = position_rssi[position_key][1]\n                    position_rssi[position_key][0] = (old_rssi * old_count + rssi) / (old_count + 1)\n                    position_rssi[position_key][1] = old_count + 1\n                else:\n                    position_rssi[position_key] = np.array([rssi, 1])\n            else:\n                position_rssi = {}\n                position_rssi[position_key] = np.array([rssi, 1])\n\n            wifi_rssi[bssid] = position_rssi\n    return wifi_rssi\ndef extract_ibeacon_rssi(mwi_datas):\n    ibeacon_rssi = {}\n    for position_key in mwi_datas:\n        # print(f'Position: {position_key}')\n        ibeacon_data = mwi_datas[position_key]['ibeacon']\n        for ibeacon_d in ibeacon_data:\n            ummid = ibeacon_d[1]\n            rssi = int(ibeacon_d[2])\n            if ummid in ibeacon_rssi:\n                position_rssi = ibeacon_rssi[ummid]\n                if position_key in position_rssi:\n                    old_rssi = position_rssi[position_key][0]\n                    old_count = position_rssi[position_key][1]\n                    position_rssi[position_key][0] = (old_rssi * old_count + rssi) / (old_count + 1)\n                    position_rssi[position_key][1] = old_count + 1\n                else:\n                    position_rssi[position_key] = np.array([rssi, 1])\n            else:\n                position_rssi = {}\n                position_rssi[position_key] = np.array([rssi, 1])\n            ibeacon_rssi[ummid] = position_rssi\n    return ibeacon_rssi\ndef extract_wifi_count(mwi_datas):\n    wifi_counts = {}\n    for position_key in mwi_datas:\n        # print(f'Position: {position_key}')\n        wifi_data = mwi_datas[position_key]['wifi']\n        count = np.unique(wifi_data[:, 2]).shape[0]\n        wifi_counts[position_key] = count\n    return wifi_counts","metadata":{"_kg_hide-input":true,"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"os.system('mkdir indoor_xy_floor')\n\n# Column names by type:\ncolacce = ['xyz_time','x_acce','y_acce','z_acce']\ncolmagn = ['xyz_time','x_magn','y_magn','z_magn']\ncolahrs = ['xyz_time','x_ahrs','y_ahrs','z_ahrs']\ncolwifi = ['t1_wifi','ssid_wifi','bssid_wifi','rssid_wifi','t2_wifi']\ncolbeac = ['t1_beac','id_beac','rssi_beac']\ncolwayp = ['t_waypoint','x_waypoint','y_waypoint']\n\nfor building in buildings:\n    train = pd.read_csv(f'train/{building}.txt',index_col=0,header=None,names=cols,dtype='object')\n    os.remove(f'train/{building}.txt') # due to the disk space limit\n\n    for floor in np.unique(train.floor):\n        xtrain = train.loc[train.floor==floor]\n        print(f'{building},{floor}',xtrain.shape)\n\n        xtrain = xtrain.fillna('')\n        xtrain[int8col] = xtrain[int8col].apply(lambda x: x.apply(lambda y: np.array(y.split(),dtype=np.int8)))\n        xtrain[int16col] = xtrain[int16col].apply(lambda x: x.apply(lambda y: np.array(y.split(),dtype=np.int16)))\n        xtrain[int32col] = xtrain[int32col].apply(lambda x: x.apply(lambda y: np.array(y.split(),dtype=np.int32)))\n        xtrain[floatcol] = xtrain[floatcol].apply(lambda x: x.apply(lambda y: np.array(y.split(),dtype=np.float32)))\n        xtrain[stringcol] = xtrain[stringcol].apply(lambda x: x.apply(lambda y: np.array(y.split())))\n\n        mwi_datas = {}\n        for path_id in xtrain.index:\n            acce_datas = np.array(xtrain.loc[path_id,colacce].apply(lambda x: pd.Series(x)).T)\n            magn_datas = np.array(xtrain.loc[path_id,colmagn].apply(lambda x: pd.Series(x)).T)\n            ahrs_datas = np.array(xtrain.loc[path_id,colahrs].apply(lambda x: pd.Series(x)).T)\n            wifi_datas = np.array(xtrain.loc[path_id,colwifi].apply(lambda x: pd.Series(x)).T)\n            ibeacon_datas = np.array(xtrain.loc[path_id,colbeac].apply(lambda x: pd.Series(x)).T)\n            posi_datas = np.array(xtrain.loc[path_id,colwayp].apply(lambda x: pd.Series(x)).T)\n            step_positions = compute_step_positions(acce_datas, ahrs_datas, posi_datas)\n            if wifi_datas.size != 0:\n                sep_tss = np.unique(wifi_datas[:, 0].astype(float))\n                wifi_datas_list = split_ts_seq(wifi_datas, sep_tss)\n                for wifi_ds in wifi_datas_list:\n                    diff = np.abs(step_positions[:, 0] - float(wifi_ds[0, 0]))\n                    index = np.argmin(diff)\n                    target_xy_key = tuple(step_positions[index, 1:3])\n                    if target_xy_key in mwi_datas:\n                        mwi_datas[target_xy_key]['wifi'] = np.append(\n                            mwi_datas[target_xy_key]['wifi'], wifi_ds, axis=0)\n                    else:\n                        mwi_datas[target_xy_key] = {\n                            'magnetic': np.zeros((0, 4)),'wifi': wifi_ds,'ibeacon': np.zeros((0, 3))}\n            if ibeacon_datas.size != 0:\n                sep_tss = np.unique(ibeacon_datas[:, 0].astype(float))\n                ibeacon_datas_list = split_ts_seq(ibeacon_datas, sep_tss)\n                for ibeacon_ds in ibeacon_datas_list:\n                    diff = np.abs(step_positions[:, 0] - float(ibeacon_ds[0, 0]))\n                    index = np.argmin(diff)\n                    target_xy_key = tuple(step_positions[index, 1:3])\n                    if target_xy_key in mwi_datas:\n                        mwi_datas[target_xy_key]['ibeacon'] = np.append(\n                            mwi_datas[target_xy_key]['ibeacon'], ibeacon_ds, axis=0)\n                    else:\n                        mwi_datas[target_xy_key] = {\n                            'magnetic': np.zeros((0, 4)),'wifi': np.zeros((0, 5)),'ibeacon': ibeacon_ds}\n            sep_tss = np.unique(magn_datas[:, 0].astype(float))\n            magn_datas_list = split_ts_seq(magn_datas, sep_tss)\n            for magn_ds in magn_datas_list:\n                diff = np.abs(step_positions[:, 0] - float(magn_ds[0, 0]))\n                index = np.argmin(diff)\n                target_xy_key = tuple(step_positions[index, 1:3])\n                if target_xy_key in mwi_datas:\n                    mwi_datas[target_xy_key]['magnetic'] = np.append(\n                        mwi_datas[target_xy_key]['magnetic'], magn_ds, axis=0)\n                else:\n                    mwi_datas[target_xy_key] = {\n                        'magnetic': magn_ds,'wifi': np.zeros((0, 5)),'ibeacon': np.zeros((0, 3))}\n\n        magnetic_strength = extract_magnetic_strength(mwi_datas)\n        wifi_counts = extract_wifi_count(mwi_datas)\n        wifi_rssi = extract_wifi_rssi(mwi_datas)\n        ibeacon_rssi = extract_ibeacon_rssi(mwi_datas)\n        print(f'This floor has {len(wifi_rssi.keys())} wifi aps')\n        print(f'This floor has {len(ibeacon_rssi.keys())} ibeacons')\n\n        x = [x for (x,y) in mwi_datas.keys()]\n        y = [y for (x,y) in mwi_datas.keys()]\n        xy = pd.DataFrame({'x':x,'y':y})\n        xy['count'] = [wifi_counts[(x,y)] if wifi_counts[(x,y)]>0 else np.nan for (x,y) in zip(xy.x,xy.y)]\n        xy['magn'] = [magnetic_strength[(x,y)] for (x,y) in zip(xy.x,xy.y)]\n        for bssid in wifi_rssi.keys():\n            xy[bssid] = [wifi_rssi[bssid][(x,y)][0] if (x,y) in wifi_rssi[bssid].keys() else np.nan for (x,y) in zip(xy.x,xy.y)]\n        for beacid in ibeacon_rssi.keys():\n            xy[beacid] = [ibeacon_rssi[beacid][(x,y)][0] if (x,y) in ibeacon_rssi[beacid].keys() else np.nan for (x,y) in zip(xy.x,xy.y)]\n        xy.to_csv(f'/kaggle/working/indoor_xy_floor/{building}_{floor}.csv')\n","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"building = '5a0546857ecc773753327266'\nfloor = \"0\"\nxy = pd.read_csv(f'indoor_xy_floor/{building}_{floor}.csv',index_col=0)\nbssid = xy.columns[4]\nbeaid = xy.columns[-1]\nxyw = xy.loc[~np.isnan(xy['count'])]\nplt.figure(figsize=(20,20))\nplt.subplot(221)\nplt.scatter(xyw.x,xyw.y,c='black')\nplt.scatter(xyw.x,xyw.y,c=xyw[bssid],cmap='viridis')\nplt.title(f'wifi id {bssid}')\nplt.subplot(222)\nplt.scatter(xyw.x,xyw.y,c='black')\nplt.scatter(xyw.x,xyw.y,c=xyw['count'],cmap='viridis')\nplt.title(f'wifi count {bssid}')\nplt.subplot(223)\nplt.scatter(xy.x,xy.y,c='black')\nplt.scatter(xy.x,xy.y,c=xy['magn'],cmap='viridis')\nplt.title(f'magnetic strength')\nplt.subplot(224)\nplt.scatter(xy.x,xy.y,c='black')\nplt.scatter(xy.x,xy.y,c=xy[beaid],cmap='viridis')\nplt.title(f'beacon id {beaid}')\nplt.show()","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"The output is stored at \"indoor_xy_floor\" @ https://www.kaggle.com/oxzplvifi/indoor-xy-floor Let's take a look at how the hosts' GitHub scripts created the data inside the \"indoor_xy_floor\" directory:","metadata":{}},{"cell_type":"markdown","source":"Step I. step detection and stride length:","metadata":{}},{"cell_type":"code","source":"# Step detection:\nstep_timestamps, step_indexs, step_acce_max_mins = compute_steps(acce_datas)\nstride_lengths = compute_stride_length(step_acce_max_mins)\n# Step detection is to detect at which timepoints the person walked:\nplt.figure(figsize=(20,10))\nplt.plot(acce_datas[:,0],acce_datas[:,3],label='az')\nplt.plot(step_acce_max_mins[:,0],step_acce_max_mins[:,1],'*-',label='max')\nplt.plot(step_acce_max_mins[:,0],step_acce_max_mins[:,2],'*-',label='min')\nplt.plot(step_acce_max_mins[:,0],step_acce_max_mins[:,3],'*-',label='variance')\nplt.legend()\nplt.title('Step detection')\nplt.show()\n# Stride length is related to the distance traveled in each step:\nplt.figure(figsize=(20,10))\nplt.plot(stride_lengths[:,0],stride_lengths[:,1],'*-')\nplt.title('Stride length')\nplt.show()","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Step II. Walking direction (angle around z):","metadata":{}},{"cell_type":"code","source":"headings = compute_headings(ahrs_datas)\nstep_headings = compute_step_heading(step_timestamps, headings)\nplt.figure(figsize=(20,10))\nplt.plot(ahrs_datas[:,0],ahrs_datas[:,1],label='x_ahrs')\nplt.plot(ahrs_datas[:,0],ahrs_datas[:,2],label='y_ahrs')\nplt.plot(ahrs_datas[:,0],ahrs_datas[:,3],label='z_ahrs')\nplt.legend()\nplt.show()\nplt.figure(figsize=(20,10))\nplt.plot(headings[:,0],headings[:,1],label='around_z')\nplt.plot(step_headings[:,0],step_headings[:,1],'*-',label='steps around_z')\nplt.legend()\nplt.show()","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Step III. Waypoint interpolation:","metadata":{}},{"cell_type":"code","source":"rel_positions = compute_rel_positions(stride_lengths, step_headings)\nstep_positions = correct_positions(rel_positions, posi_datas)\nplt.figure(figsize=(20,10))\nplt.plot(rel_positions[:,0],rel_positions[:,1],'*-',label='stride*-sin')\nplt.plot(rel_positions[:,0],rel_positions[:,2],'*-',label='stride*cos')\nplt.legend()\nplt.show()\n\nplt.figure(figsize=(20,10))\nplt.plot(posi_datas[:,0],posi_datas[:,1],'*-',label='x waypoint')\nplt.plot(step_positions[:,0],step_positions[:,1],'+-',label='x position')\nplt.legend()\nplt.show()\n\nplt.figure(figsize=(20,10))\nplt.plot(posi_datas[:,0],posi_datas[:,2],'*-',label='y waypoint')\nplt.plot(step_positions[:,0],step_positions[:,2],'+-',label='y position')\nplt.legend()\nplt.show()","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"The training set thus consists of a 4D grid map of (X,Y,Floor,RSSI) values. On the other hand, the testing set is 2D consisting of (time,RSSI) values.","metadata":{}},{"cell_type":"markdown","source":"Thank you for reading! Let me know if you have any question or suggestion.","metadata":{}},{"cell_type":"code","source":"","metadata":{"trusted":true},"execution_count":null,"outputs":[]}]}