{"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\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\n","metadata":{"_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","_kg_hide-output":true,"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","papermill":{"duration":9.495204,"end_time":"2021-02-27T17:22:28.195457","exception":false,"start_time":"2021-02-27T17:22:18.700253","status":"completed"},"scrolled":true,"tags":[],"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Thank you for reading, I would like to share the procedure I used to obtain my current best single-model results. Basically it consists of the following steps:\n1. Fix data quality issues and reformat the training data through C++ custom code @ https://www.kaggle.com/oxzplvifi/indoor-cpp The procedure is shown @ https://www.kaggle.com/oxzplvifi/indoor-preprocessing-and-eda And the result is located @ https://www.kaggle.com/oxzplvifi/indoor-xy-floor\n2. Use the hosts' GitHub scripts to generate a grid map of (floor,X,Y,rssi,magnetism) values from the \"step detection\" and \"orientation detection\" for training the models. The hosts' scripts are available @ https://github.com/location-competition/indoor-location-competition-20 Copyright (c) 2017-2020 XYZ10, Inc. https://dangwu.com/\n3. Predict the floor for each path through KMeans and GBM. The procedure is shown @ https://www.kaggle.com/oxzplvifi/indoor-kmeans-gbm-floor-prediction\n4. Predict the XY coordinates using the GBM model shared by BIZEN (hiro5299834) in his very nice work @ https://www.kaggle.com/hiro5299834/wifi-features-with-lightgbm-kfold , but using the (floor,X,Y,rssi,magnetism) grid map as input instead of the raw rssi data, and stratifying by floor.\n5. The postprocessing procedure is explained near the end of this notebook.","metadata":{"papermill":{"duration":0.112591,"end_time":"2021-02-27T17:22:28.422593","exception":false,"start_time":"2021-02-27T17:22:28.310002","status":"completed"},"tags":[]}},{"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":{"papermill":{"duration":45.003013,"end_time":"2021-02-27T17:23:13.533064","exception":false,"start_time":"2021-02-27T17:22:28.530051","status":"completed"},"tags":[],"trusted":true},"execution_count":null,"outputs":[]},{"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# Reformat the testing set:\nos.system('mkdir test')\nfor i,(path_filename,building) in enumerate(zip(test_buildings['file'],test_buildings['building'])):\n    er=os.system(f'./preprocess {path_filename} test {building} {0}') #since we do not know the floor, I put 0.\n    if(er): print(\"Error:\",path_filename)","metadata":{"papermill":{"duration":439.856888,"end_time":"2021-02-27T17:30:33.503456","exception":false,"start_time":"2021-02-27T17:23:13.646568","status":"completed"},"tags":[],"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Acceleration, magnetic and orientation testing data:\nos.system('mkdir indoor_testing_accel')\nos.system(\"g++ /kaggle/input/indoor-cpp/2_preprocess_accel.cpp -std=c++11 -o preprocess_accel\")\nfor building in buildings:\n    os.system(f'./preprocess_accel {building}')","metadata":{"papermill":{"duration":57.340848,"end_time":"2021-02-27T17:31:30.957517","exception":false,"start_time":"2021-02-27T17:30:33.616669","status":"completed"},"tags":[],"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Wifi testing data:\nos.system('mkdir test_wifi')\nos.system(\"g++ /kaggle/input/indoor-cpp/2_preprocess_wifi.cpp -std=c++11 -o preprocess_wifi\")\nfor building in buildings:\n    os.system(f'./preprocess_wifi {building}')","metadata":{"papermill":{"duration":31.422014,"end_time":"2021-02-27T17:32:02.496614","exception":false,"start_time":"2021-02-27T17:31:31.074600","status":"completed"},"tags":[],"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Procedure to produce floor results through KMeans+GBM @ https://www.kaggle.com/oxzplvifi/indoor-kmeans-gbm-floor-prediction","metadata":{"papermill":{"duration":0.111555,"end_time":"2021-02-27T17:32:02.719015","exception":false,"start_time":"2021-02-27T17:32:02.607460","status":"completed"},"tags":[]}},{"cell_type":"code","source":"floors = pd.read_csv('/kaggle/input/indoor-xy-floor/result_floor_feb22.csv',index_col=0)\nfloors.head()","metadata":{"papermill":{"duration":0.144027,"end_time":"2021-02-27T17:32:02.976471","exception":false,"start_time":"2021-02-27T17:32:02.832444","status":"completed"},"tags":[],"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Compute running average on the magnetic field and interpolate into the wifi coordinates of the testing set:","metadata":{"papermill":{"duration":0.157584,"end_time":"2021-02-27T17:32:03.244698","exception":false,"start_time":"2021-02-27T17:32:03.087114","status":"completed"},"tags":[]}},{"cell_type":"code","source":"from scipy.interpolate import interp1d\nfrom scipy.ndimage.filters import uniform_filter1d\n\nfor building in buildings:\n    print(building)\n    \n    # Testing set:\n    tfw = pd.read_csv(f'test_wifi/{building}.txt')\n    tfw['path_id'] = [building+'_'+x for x in tfw.path_id]\n    tfw = tfw.merge(floors['floor'],left_on='path_id',right_index=True)\n    tfw = tfw.pivot_table(index=['path_id','t1_wifi','floor'],columns='bssid_wifi',values='rssid_wifi')\n    tfw.fillna(-99.0,inplace=True)\n    tfw.insert(0,'count',tfw.apply(lambda x: np.sum(x!=-99.0),axis=1))\n    tfw.insert(0,'magn',0)\n    paths = np.unique([x[0] for x in tfw.index])\n    \n    # Read magnetic signature:\n    tfm = pd.read_csv(f'indoor_testing_accel/{building}.txt',index_col=0)\n    tfm.index = [building+'_'+x for x in tfm.index]\n    tfm['magn'] = np.sqrt(tfm['x_magn']**2 + tfm['y_magn']**2 + tfm['z_magn']**2)\n    tfm = tfm.sort_values(by=['xyz_time'])\n    \n    # Interpolate the magnetic signature for each wifi timestamp:\n    for path in paths:\n        tm = tfm.loc[path,'xyz_time']\n        m = uniform_filter1d(tfm.loc[path,'magn'], size=3, mode='reflect')\n        fm = interp1d(tm,m,kind='linear',fill_value=(m[0],m[-1]),bounds_error=False) #fill_value=\"extrapolate\"\n        tw = [x[0] for x in tfw.loc[path].index]\n        tfw.loc[path,'magn'] = fm(tw)\n    \n    # Save output:\n    tfw = tfw.reset_index()\n    tfw.to_csv(f'{building}.csv')\n#     break","metadata":{"papermill":{"duration":144.650422,"end_time":"2021-02-27T17:34:28.028790","exception":false,"start_time":"2021-02-27T17:32:03.378368","status":"completed"},"tags":[],"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"import matplotlib.pyplot as plt\n\nplt.figure(figsize=(10,5))\nplt.plot(tfm.loc[path,'xyz_time'],m)\nplt.plot(tw,fm(tw),'+r')\nplt.show()","metadata":{"papermill":{"duration":0.366246,"end_time":"2021-02-27T17:34:28.521514","exception":false,"start_time":"2021-02-27T17:34:28.155268","status":"completed"},"tags":[],"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"GBM parameters from BIZEN (hiro5299834) in his very nice work @ https://www.kaggle.com/hiro5299834/wifi-features-with-lightgbm-kfold","metadata":{"papermill":{"duration":0.183827,"end_time":"2021-02-27T17:34:28.828193","exception":false,"start_time":"2021-02-27T17:34:28.644366","status":"completed"},"tags":[]}},{"cell_type":"code","source":"lgb_params = {\n    'objective': 'root_mean_squared_error',\n    'boosting_type': 'gbdt', #default=gbdt\n    'n_estimators': 50000, #default=100*num_class\n    'learning_rate': 0.1, #default=0.1\n    'num_leaves': 90, #default=31\n    'colsample_bytree': 0.4, #default=1.0 (% of cols selected)\n    'subsample': 0.6, #default=1.0 (% of rows selected)\n    'subsample_freq': 2, #default=0 (perform bagging every kth iteration)\n    'bagging_seed': 42, #default=3\n    'reg_alpha': 10, #default=0.0 (L1 regularization)\n    'reg_lambda': 2, #default=0.0 (L2 regularization)\n    'random_state': 42, #default=None\n    'n_jobs': -1, #default=0\n#     'device':'gpu'\n}","metadata":{"papermill":{"duration":0.131669,"end_time":"2021-02-27T17:34:29.082021","exception":false,"start_time":"2021-02-27T17:34:28.950352","status":"completed"},"tags":[],"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Procedure to generate training tables and EDA @ https://www.kaggle.com/oxzplvifi/indoor-preprocessing-and-eda","metadata":{"papermill":{"duration":0.120094,"end_time":"2021-02-27T17:34:29.327185","exception":false,"start_time":"2021-02-27T17:34:29.207091","status":"completed"},"tags":[]}},{"cell_type":"code","source":"import lightgbm as lgb\nfrom sklearn.model_selection import StratifiedKFold\n\nresult = pd.DataFrame(columns=['path_id','t1_wifi','x','y'])\n\nfor building in buildings:\n\n    # Training set:\n    xyw = pd.DataFrame()\n    for floor in np.arange(-3,10):\n        file = f'/kaggle/input/indoor-xy-floor/{building}_{floor}.csv'\n        if Path(file).is_file():\n            xyi = pd.read_csv(file,index_col=0)\n            bcols = [c for c in xyi.columns if len(c.split('_'))==3] #beacon cols\n            wcols = [c for c in xyi.columns if c not in ['x','y','count','magn']+bcols] #wifi cols\n            xyi = xyi.loc[~np.isnan(xyi['count']),['x','y','count','magn']+wcols]\n            xyi.insert(0,'floor',floor)\n            if(len(xyw)):\n                xyw = xyw.merge(xyi,how='outer')\n            else: xyw = xyi\n    xyw.replace(np.nan,-99.0,inplace=True)\n    wcols = list(xyw.columns[5:])\n\n    # Testing set:\n    tfw = pd.read_csv(f'{building}.csv',index_col=0)\n    tfw.set_index(['path_id','t1_wifi'],inplace=True)\n    tfw = tfw.reindex(columns=['floor','magn']+wcols,fill_value=-99.0)\n\n    # Arrays:\n    dfmat = np.array(xyw[['floor','magn']+wcols])\n    mtest = np.array(tfw)\n    labs = np.array(xyw[['x','y']])\n    yvalid = pd.DataFrame(np.zeros([len(labs),2]),index=xyw.index,columns=['x','y'])\n    ytest = pd.DataFrame(np.zeros([len(tfw),2]),index=tfw.index,columns=['x','y'])\n\n    # K-fold CV of coordinates:\n    seeds, folds = 1, 10\n    skf = StratifiedKFold(n_splits=folds,random_state=0,shuffle=True)\n    for fold, (idt,idv) in enumerate(skf.split(dfmat,xyw['floor'])):\n        print('\\r',f'{fold}',end='\\t')\n        mtrain, mvalid = dfmat[idt], dfmat[idv]\n        ltrain, lvalid = labs[idt], labs[idv]\n        # X prediction:\n        modelx = lgb.LGBMRegressor(**lgb_params)\n        modelx.fit(mtrain,ltrain[:,0],eval_set=[(mvalid,lvalid[:,0])],\n            eval_metric='rmse',early_stopping_rounds=10,verbose=False)\n        yvalid.loc[xyw.index[idv],'x'] = modelx.predict(mvalid)\n#         modelx.booster_.save_model(f'lgb/{building}_{fold}_x.txt',num_iteration=modelx.best_iteration_)\n#         modelx = lgb.Booster(model_file=f'lgb/{building}_{fold}_x.txt')\n        ytest['x'] += modelx.predict(mtest) / folds\n        # Y prediction:\n        modely = lgb.LGBMRegressor(**lgb_params)\n        modely.fit(mtrain,ltrain[:,1],eval_set=[(mvalid,lvalid[:,1])],\n            eval_metric='rmse',early_stopping_rounds=10,verbose=False)\n#         modely.booster_.save_model( f'lgb/{building}_{fold}_y.txt',num_iteration=modely.best_iteration_)\n#         modely = lgb.Booster(model_file=f'lgb/{building}_{fold}_y.txt')\n        yvalid.loc[xyw.index[idv],'y'] = modely.predict(mvalid)\n        ytest['y'] += modely.predict(mtest) / folds\n\n    # Performance:\n    yvalid['xtruth'] = xyw['x']\n    yvalid['ytruth'] = xyw['y']\n    xrmse = np.sqrt((yvalid['x']-yvalid['xtruth'])**2)\n    yrmse = np.sqrt((yvalid['y']-yvalid['ytruth'])**2)\n    print(building, f'xrmse = {np.mean(xrmse)}, yrmse = {np.mean(yrmse)}')\n\n    # Prediction:\n    result = pd.concat([result,ytest.reset_index()])\n#     break\n\nresult.set_index('path_id',inplace=True)","metadata":{"papermill":{"duration":15453.441642,"end_time":"2021-02-27T21:52:02.887878","exception":false,"start_time":"2021-02-27T17:34:29.446236","status":"completed"},"tags":[],"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"The predicted XY values correspond to the WIFI timestamps, not to the submission timepoints. Therefore, we must interpolate them. Prior to interpolation, we employ the hosts' GitHub scripts containing the \"step detection\" and \"orientation detection\" procedures in order to improve the prediction @ https://github.com/location-competition/indoor-location-competition-20 Copyright (c) 2017-2020 XYZ10, Inc. https://dangwu.com/","metadata":{"papermill":{"duration":0.234942,"end_time":"2021-02-27T21:52:03.322137","exception":false,"start_time":"2021-02-27T21:52:03.087195","status":"completed"},"tags":[]}},{"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\n\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\n\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\n\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,"papermill":{"duration":0.423484,"end_time":"2021-02-27T21:52:03.939094","exception":false,"start_time":"2021-02-27T21:52:03.515610","status":"completed"},"tags":[],"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Load submission file to see which timepoints are required:\nsample_submission = pd.read_csv(inpath+'sample_submission.csv')\nsample_submission['building'] = [x.split('_')[0] for x in sample_submission['site_path_timestamp']]\nsample_submission['path_id'] = [x.split('_')[1] for x in sample_submission['site_path_timestamp']]\nsample_submission['timestamp'] = [x.split('_')[2] for x in sample_submission['site_path_timestamp']]\nsamples = pd.DataFrame(sample_submission.groupby(['building','path_id'])['timestamp'].apply(lambda x: list(x)))\nbuildings = np.unique([x[0] for x in samples.index])\nsamples.head()","metadata":{"papermill":{"duration":0.362381,"end_time":"2021-02-27T21:52:04.493525","exception":false,"start_time":"2021-02-27T21:52:04.131144","status":"completed"},"tags":[],"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"The following is my postprocessing procedure:\n1. First we remove the outliers by comparing each value to its running average (i.e. we get the distribution of differences between the original values and the running average and use 3 standard deviations as the outlier threshold).\n2. Then we run the hosts' scripts to detect steps and orientation.\n3. Then we do another running average. But we do not keep all values because if we do that, then the time interval between consecutive timepoints will be too short to make any use of the hosts' script to take into account the steps and orientation. So instead we keep a third of the data points, so that each point is an average of three original data prediction timepoints.\n4. Then we run the hosts' scripts to correct the position using the step and orientation.\n5. Finally, we interpolate to the desired timepoints.","metadata":{"papermill":{"duration":0.196003,"end_time":"2021-02-27T21:52:04.887292","exception":false,"start_time":"2021-02-27T21:52:04.691289","status":"completed"},"tags":[]}},{"cell_type":"code","source":"from scipy.interpolate import interp1d\nfrom scipy.ndimage.filters import uniform_filter1d\n\ncolacce = ['xyz_time','x_acce','y_acce','z_acce']\ncolahrs = ['xyz_time','x_ahrs','y_ahrs','z_ahrs']\n\nfor building in buildings:\n    print(building)\n    paths = samples.loc[building].index\n    # Acceleration info:\n    tfm = pd.read_csv(f'indoor_testing_accel/{building}.txt',index_col=0)\n    for path_id in paths:\n        # Original predicted values:\n        xy = result.loc[building+'_'+path_id]\n        tfmi = tfm.loc[path_id]\n        acce_datas = np.array(tfmi[colacce],dtype=np.float)\n        ahrs_datas = np.array(tfmi[colahrs],dtype=np.float)\n        posi_datas = np.array(xy[['t1_wifi','x','y']],dtype=np.float)\n        # 1.Outlier removal:\n        xyout = uniform_filter1d(posi_datas,size=3,axis=0,mode='reflect')\n        xydiff = np.abs(posi_datas-xyout)\n        xystd = np.std(xydiff,axis=0)*3\n        posi_datas = posi_datas[(xydiff[:,1]<xystd[1])&(xydiff[:,2]<xystd[2])]\n        # 2.1 Step detection:\n        step_timestamps, step_indexs, step_acce_max_mins = compute_steps(acce_datas)\n        stride_lengths = compute_stride_length(step_acce_max_mins)\n        # 2.2 Orientation detection:\n        headings = compute_headings(ahrs_datas)\n        step_headings = compute_step_heading(step_timestamps, headings)\n        rel_positions = compute_rel_positions(stride_lengths, step_headings)\n        # 3. Running average:\n        posi_datas = uniform_filter1d(posi_datas,size=3,axis=0,mode='reflect')[0::3,:]\n        # The 1st prediction timepoint should be earlier than the 1st step timepoint.\n        rel_positions = rel_positions[rel_positions[:,0]>posi_datas[0,0],:]\n        # If two consecutive predictions are in-between two step datapoints,\n        # the last one is removed, causing error (in the \"split_ts_seq\" function).\n        posi_index = [np.searchsorted(rel_positions[:,0], x, side='right') for x in posi_datas[:,0]]\n        u, i1, i2 = np.unique(posi_index, return_index=True, return_inverse=True)\n        posi_datas = np.vstack([np.mean(posi_datas[i2==i],axis=0) for i in np.unique(i2)])\n        # 4.Position correction:\n        step_positions = correct_positions(rel_positions, posi_datas)\n        # 5.Interpolate for timestamps in the testing set:\n        t = step_positions[:,0]\n        x = step_positions[:,1]\n        y = step_positions[:,2]\n        fx = interp1d(t, x, kind='linear', fill_value=(x[0],x[-1]), bounds_error=False) #fill_value=\"extrapolate\"\n        fy = interp1d(t, y, kind='linear', fill_value=(y[0],y[-1]), bounds_error=False)\n        # 6.Output result:\n        t0 = np.array(samples.loc[(building,path_id),'timestamp'],dtype=np.float64)\n        sample_submission.loc[(sample_submission.building==building)&(sample_submission.path_id==path_id),'x'] = fx(t0)\n        sample_submission.loc[(sample_submission.building==building)&(sample_submission.path_id==path_id),'y'] = fy(t0)\n        sample_submission.loc[(sample_submission.building==building)&(sample_submission.path_id==path_id),'floor'] = floors.loc[building+'_'+path_id,'floor']\n#         break\n#     break\n\nsample_submission[['site_path_timestamp','floor','x','y']].to_csv('submission.csv',index=False)\nsample_submission.head()","metadata":{"papermill":{"duration":750.175198,"end_time":"2021-02-27T22:04:35.258139","exception":false,"start_time":"2021-02-27T21:52:05.082941","status":"completed"},"tags":[],"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Visualize the postprocessing procedure:","metadata":{"papermill":{"duration":0.201458,"end_time":"2021-02-27T22:04:35.660872","exception":false,"start_time":"2021-02-27T22:04:35.459414","status":"completed"},"tags":[]}},{"cell_type":"code","source":"# Original predicted values:\nxy = result.loc[building+'_'+path_id]\nposi_datas = np.array(xy[['t1_wifi','x','y']],dtype=np.float)\n\n# Outlier removal (|Z-score|>3 of difference with running average):\nxyout = uniform_filter1d(posi_datas,size=3,axis=0,mode='reflect')\nxydiff = np.abs(posi_datas-xyout)\nxystd = np.std(xydiff,axis=0)*3\nposi_datas = posi_datas[(xydiff[:,1]<xystd[1])&(xydiff[:,2]<xystd[2])]\n\nplt.figure(figsize=(20,10))\nplt.plot(xy.t1_wifi,xy.x,'+-b')\nplt.plot(posi_datas[:,0],posi_datas[:,1],'+-r')\nplt.show()\n\nplt.figure(figsize=(20,10))\nplt.plot(xy.t1_wifi,xy.y,'+-b')\nplt.plot(posi_datas[:,0],posi_datas[:,2],'+-r')\nplt.show()","metadata":{"papermill":{"duration":0.580423,"end_time":"2021-02-27T22:04:36.442686","exception":false,"start_time":"2021-02-27T22:04:35.862263","status":"completed"},"tags":[],"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Step detection:\ntfmi = tfm.loc[path_id]\nstep_timestamps, step_indexs, step_acce_max_mins = compute_steps(acce_datas)\nstride_lengths = compute_stride_length(step_acce_max_mins)\nplt.figure(figsize=(20,10))\n\n# Step detection is to detect at which timepoints the person walked:\nplt.plot(tfmi['xyz_time'],tfmi['z_acce'],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\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":{"papermill":{"duration":1.531455,"end_time":"2021-02-27T22:04:38.186906","exception":false,"start_time":"2021-02-27T22:04:36.655451","status":"completed"},"tags":[],"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Orientation detection:\nheadings = compute_headings(ahrs_datas)\nstep_headings = compute_step_heading(step_timestamps, headings)\nrel_positions = compute_rel_positions(stride_lengths, step_headings)\n\nplt.figure(figsize=(20,10))\nplt.plot(tfmi['xyz_time'],tfmi['x_ahrs'],label='x_ahrs')\nplt.plot(tfmi['xyz_time'],tfmi['y_ahrs'],label='y_ahrs')\nplt.plot(tfmi['xyz_time'],tfmi['z_ahrs'],label='z_ahrs')\nplt.legend()\nplt.show()\n\nplt.figure(figsize=(20,10))\nplt.plot(headings[:,0],headings[:,1],label='around_z')\nplt.plot(step_headings[:,0],step_headings[:,1],'*-',label='around_z')\nplt.legend()\nplt.show()","metadata":{"papermill":{"duration":0.950363,"end_time":"2021-02-27T22:04:39.368972","exception":false,"start_time":"2021-02-27T22:04:38.418609","status":"completed"},"tags":[],"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Position correction:\nposi_avg = uniform_filter1d(posi_datas,size=3,axis=0,mode='reflect')[0::3,:]\n#The first prediction timepoint should be earlier than the first step timepoint.\nrel_positions = rel_positions[rel_positions[:,0]>posi_avg[0,0],:]\n#If two consecutive predictions are in-between two step datapoints,\n#the last one is removed, causing error (in the \"split_ts_seq\" function).\nposi_index = [np.searchsorted(rel_positions[:,0], x, side='right') for x in posi_avg[:,0]]\nu, i1, i2 = np.unique(posi_index, return_index=True, return_inverse=True)\nposi_avg = np.vstack([np.mean(posi_avg[i2==i],axis=0) for i in np.unique(i2)])\nstep_positions = correct_positions(rel_positions, posi_avg)\n\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(posi_avg[:,0],posi_avg[:,1],'*-',label='x waypoint avg')\nplt.plot(step_positions[:,0],step_positions[:,1],'-',label='x step')\nplt.legend()\nplt.show()\n\nplt.figure(figsize=(20,10))\nplt.plot(posi_datas[:,0],posi_datas[:,2],'*-',label='y waypoint')\nplt.plot(posi_avg[:,0],posi_avg[:,2],'*-',label='y waypoint avg')\nplt.plot(step_positions[:,0],step_positions[:,2],'-',label='y step')\nplt.legend()\nplt.show()","metadata":{"papermill":{"duration":0.909495,"end_time":"2021-02-27T22:04:40.506218","exception":false,"start_time":"2021-02-27T22:04:39.596723","status":"completed"},"tags":[],"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Thank you for reading! Let me know if you have any questions, and I look forward to any suggestions!","metadata":{"papermill":{"duration":0.240975,"end_time":"2021-02-27T22:04:40.991525","exception":false,"start_time":"2021-02-27T22:04:40.750550","status":"completed"},"tags":[]}},{"cell_type":"code","source":"","metadata":{"papermill":{"duration":0.230953,"end_time":"2021-02-27T22:04:41.463343","exception":false,"start_time":"2021-02-27T22:04:41.232390","status":"completed"},"tags":[]},"execution_count":null,"outputs":[]}]}