{"metadata":{"kernelspec":{"language":"python","display_name":"Python 3","name":"python3"},"language_info":{"pygments_lexer":"ipython3","nbconvert_exporter":"python","version":"3.6.4","file_extension":".py","codemirror_mode":{"name":"ipython","version":3},"name":"python","mimetype":"text/x-python"}},"nbformat_minor":4,"nbformat":4,"cells":[{"cell_type":"markdown","source":"This method is a post-processing applied to some of our team's models.\n\nThis notebook incorporated geojson into [Saito's post-processing](https://www.kaggle.com/saitodevel01/indoor-post-processing-by-cost-minimization) and solved it as a nonlinear cost minimization problem.\n\nTo combine Saito's post-processing with geojson, I defined cost function as follows,\n$$\nL(X_{1:N}) = \\sum_{i=1}^{N} \\alpha_i \\| X_i - \\hat{X}_i \\|^2 + \\sum_{i=1}^{N-1} \\beta_i \\| (X_{i+1} - X_{i}) - \\Delta \\hat{X}_i \\|^2 + \\sum_{i=1}^{N-1} \\ \\sum_{j=1}^{M} \\gamma_{ij} \\| Dist(X_{ij}) \\| ^2\n$$\nwhere $X_{ij}$ is an evenly spaced point on the line segment $X_{i}$$X_{i+1}$, and $Dist(X_{ij})$ represents the distance to the wall, which is represented as the following distance image created from geojson.","metadata":{}},{"cell_type":"code","source":"import os\nfrom glob import glob\nimport numpy as np\nimport math\nimport pandas as pd\nfrom tqdm.notebook import tqdm\nimport copy\nimport pickle as pkl\nimport json\nimport matplotlib.pyplot as plt\nimport multiprocessing\n\nfrom sklearn.model_selection import GroupKFold\nfrom scipy.ndimage import gaussian_filter1d\nfrom scipy.spatial.distance import cdist\nimport scipy\n\nimport torch\nfrom torch import nn\nimport torch.nn.functional as F\nfrom torch.utils.data import Dataset, DataLoader\nfrom torch.optim import AdamW, lr_scheduler\n\nimport gc\nimport warnings\nwarnings.filterwarnings('ignore')\n\nimport glob\nfrom scipy import optimize\nfrom scipy.optimize import minimize\nimport math","metadata":{"_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","papermill":{"duration":2.319228,"end_time":"2021-03-19T16:05:54.752689","exception":false,"start_time":"2021-03-19T16:05:52.433461","status":"completed"},"tags":[],"_kg_hide-input":true,"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"site = '5d27096c03f801723c31e5e0'\nfloorNo = 'F3'\n\nfloor = plt.imread(f\"../input/indoor-location-navigation/metadata/{site}/{floorNo}/floor_image.png\")\nbw = plt.imread(f\"../input/indoorlocationnavigationbwdist/metadata/{site}/{floorNo}/geojson_map_bw_filled_x16.png\")\ndist = plt.imread(f\"../input/indoorlocationnavigationbwdist/metadata/{site}/{floorNo}/geojson_map_cv_dist_normalized_x16.png\")\n\nplt.figure(figsize=(18,12))\nplt.subplot(1,3,1)\nplt.title('Map')\nplt.imshow(floor)\nplt.subplot(1,3,2)\nplt.title('Wall')\nplt.imshow(bw[::-1,:])\nplt.subplot(1,3,3)\nplt.title('Dist')\nplt.imshow(dist[::-1,:])\nplt.show()","metadata":{"_kg_hide-input":true,"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"To make it easier to converge to the correct position, the scale of the distance image is multiplied by 16.\n\nSince this equation is nonlinear, it can be solved using the COBYLA method of scipy's minimize function.\n\nThe raw prediction from machine learning is based on [Tri's transformer model](https://www.kaggle.com/shinomoriaoshi/iln-transformer-inference).\n\n\n\n* public score : 5.17309 -> 4.38225\n* private score : 5.50094 -> 4.95899","metadata":{}},{"cell_type":"code","source":"train_waypoints = pd.read_csv('../input/indoor-location-train-waypoints/train_waypoints.csv')\nss = pd.read_csv('../input/indoorlocationnavigationbwdist/metadata/Tri_raw_prediction.csv')\n\ntmp = ss['site_path_timestamp'].apply(lambda s : pd.Series(s.split('_')))\nss['site'] = tmp[0]\nss['path'] = tmp[1]\nss['timestamp'] = tmp[2].astype(float)\n\nss = ss.merge(train_waypoints[['site','floor','floorNo']].drop_duplicates())","metadata":{"papermill":{"duration":0.268891,"end_time":"2021-03-19T16:08:09.525941","exception":false,"start_time":"2021-03-19T16:08:09.25705","status":"completed"},"tags":[],"_kg_hide-input":false,"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Saito's post-processing with geojson","metadata":{"papermill":{"duration":0.017685,"end_time":"2021-03-19T16:08:09.561605","exception":false,"start_time":"2021-03-19T16:08:09.54392","status":"completed"},"tags":[]}},{"cell_type":"code","source":"!git clone --depth 1 https://github.com/location-competition/indoor-location-competition-20 indoor_location_competition_20\n!rm -rf indoor_location_competition_20/data\n\nfrom indoor_location_competition_20.io_f import read_data_file\nimport indoor_location_competition_20.compute_f as compute_f","metadata":{"papermill":{"duration":59.174274,"end_time":"2021-03-19T16:09:10.982742","exception":false,"start_time":"2021-03-19T16:08:11.808468","status":"completed"},"tags":[],"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def compute_rel_positions(acce_datas, ahrs_datas):\n    step_timestamps, step_indexs, step_acce_max_mins = compute_f.compute_steps(acce_datas)\n    headings = compute_f.compute_headings(ahrs_datas)\n    stride_lengths = compute_f.compute_stride_length(step_acce_max_mins)\n    step_headings = compute_f.compute_step_heading(step_timestamps, headings)\n    rel_positions = compute_f.compute_rel_positions(stride_lengths, step_headings)\n    return rel_positions","metadata":{"papermill":{"duration":0.092476,"end_time":"2021-03-19T16:09:11.145491","exception":false,"start_time":"2021-03-19T16:09:11.053015","status":"completed"},"tags":[],"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def calc_points(x1_realnum, y1_realnum, x2_realnum, y2_realnum):\n    \n    x1 = math.floor(x1_realnum)\n    y1 = math.floor(y1_realnum)\n    x2 = math.floor(x2_realnum)\n    y2 = math.floor(y2_realnum)\n    \n    points = []\n    \n    if (x1 == x2) and (y1 == y2):\n        points.append([x1_realnum,y1_realnum])\n        points.append([x2_realnum,y2_realnum])\n        return np.array(points)\n\n    x1, y1 = x1_realnum, y1_realnum\n    x2, y2 = x2_realnum, y2_realnum\n    \n    step_x = np.sign(x2 - x1)\n    step_y = np.sign(y2 - y1)\n    a = (y2-y1)/(x2-x1)\n    dx = 1/((a**2+1)**0.5)\n    dy = a/((a**2+1)**0.5)\n    n = math.floor(((y2-y1)**2 + (x2-x1)**2)**0.5)\n    \n    x = x1_realnum\n    y = y1_realnum\n    points.append([x1_realnum, y1_realnum])\n    for i in range(n-1):\n        x = x + step_x * dx\n        y = y + step_y * dy\n        points.append([x, y])\n    points.append([x2_realnum, y2_realnum])\n    \n    return np.array(points)","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def f(v,path_df, alpha_, beta_, xy_hat_, delta_xy_hat_, dist, delta_t): \n\n    length = int(len(v)/2)-1\n    dx,dy = v[-2], v[-1]\n    x = v[0:-2:2] + dx\n    y = v[1:-1:2] + dy\n    v = v[0:-2]\n\n    delta_x = delta_xy_hat_[0::2]\n    delta_y = delta_xy_hat_[1::2]\n    \n    cost_alpha = sum(alpha_ * (v - xy_hat_)**2)\n    cost_beta = sum(beta_ * (v[2:2*length] - v[0:2*length-2] - delta_xy_hat_)**2)\n    \n    cost_gamma = 0\n    for i in range(length-1):\n        points = calc_points(x[i],y[i],x[i+1],y[i+1])\n        for j, point in enumerate(points):\n            gamma = 100\n            if j==0 or j==len(points)-1:\n                start_and_end_coeff = 20\n            else:\n                start_and_end_coeff = 1\n            \n            if point[0] < 0 or point[1] < 0 or point[0] > dist.shape[1] or point[1] > dist.shape[0]:\n                cost_gamma += gamma * start_and_end_coeff\n            else:\n                try:\n                    px = point[0]\n                    py = point[1]\n                    pxf = math.floor(px)\n                    pyf = math.floor(py)\n                    dist_bilinear = (pxf+1-px)*((pyf+1-py)*dist[pyf][pxf]+(py-pyf)*dist[pyf+1][pxf])+(px-pxf)*((pyf+1-py)*dist[pyf][pxf+1]+(py-pyf)*dist[pyf+1][pxf+1])\n                    cost_gamma += start_and_end_coeff * gamma*dist_bilinear\n                except:\n                    cost_gamma += gamma * start_and_end_coeff\n        \n    return cost_alpha + cost_beta + cost_gamma","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"show_path_list = ['07db4eec1fc8df5040d86602','0412d582bb8a2c89400a1ffb', '08a29ac71fe460a44ca73b38', '0969b6db91bb1697f99a8a21', '0fb75bf676e1cddace3590d5', '1747b3e80ba80910158c0cef', '1e85728385bc83e6bb650798', '270a66500675e14fb3154c63','947e17f82dbddfbdb4cb2447']","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def correct_path(args):\n    path, path_df = args\n    site = path_df['site'].iloc[0]\n    floorNo = path_df['floorNo'].iloc[0]\n    \n    T_ref  = path_df['timestamp'].values\n    xy_hat = path_df[['x', 'y']].values\n    \n    example_path = glob.glob(f'../input/indoor-location-navigation/**/{path}.txt', recursive=True)\n    example = read_data_file(example_path[0])\n    rel_positions = compute_rel_positions(example.acce, example.ahrs)\n\n    if T_ref[-1] > rel_positions[-1, 0]:\n        rel_positions = [np.array([[0, 0, 0]]), rel_positions, np.array([[T_ref[-1], 0, 0]])]\n    else:\n        rel_positions = [np.array([[0, 0, 0]]), rel_positions]\n    rel_positions = np.concatenate(rel_positions)\n    \n    T_rel = rel_positions[:, 0]\n    \n    try:\n        delta_xy_hat = np.diff(scipy.interpolate.interp1d(T_rel, np.cumsum(rel_positions[:, 1:3],\n        axis=0), axis=0)(T_ref), axis=0)\n    except:\n        return path_df.drop(['site','path','timestamp'],axis=1)\n    \n    if sum(delta_xy_hat.flatten()) == 0:\n        return path_df.drop(['site','path','timestamp'],axis=1)\n    \n    N = xy_hat.shape[0]\n    delta_t = np.diff(T_ref)\n    alpha = (8.1)**(-2) * np.ones(N)\n    beta  = (0.4 + 0.4 * 1e-3 * delta_t)**(-2)\n    \n    alpha_ = np.ravel(np.stack([alpha, alpha],1))\n    beta_ = np.ravel(np.stack([beta, beta],1))\n    \n    xy_hat_ = np.ravel(xy_hat)\n    delta_xy_hat_ = np.ravel(delta_xy_hat)\n    \n    dist = plt.imread(f\"../input/indoorlocationnavigationbwdist/metadata/{site}/{floorNo}/geojson_map_cv_dist_normalized_x16.png\")\n    \n    x = xy_hat_[0::2]\n    y = xy_hat_[1::2]\n    \n    initial_v = np.ravel(path_df[['x', 'y']].values)\n    initial_v = np.append(initial_v, [0,0])\n    \n    v = minimize( f, x0=initial_v*16, args=(path_df, alpha_, beta_, xy_hat_*16, delta_xy_hat_*16, dist, delta_t), method=\"cobyla\",options={'rhobeg': 10.0, 'maxiter': 1000, 'disp': False, 'catol': 0.0002})\n    \n    cost = f(v['x']*16,path_df, alpha_, beta_, xy_hat_*16, delta_xy_hat_*16, dist, delta_t)\n    \n    dx,dy = v['x'][-2], v['x'][-1]\n    x_ = v['x'][0:-2:2] + dx\n    y_ = v['x'][1:-1:2] + dy\n    \n    # visualize \n    if path in show_path_list:\n        print(site,floorNo,path,cost, flush=True)\n        plt.figure(figsize=(12,12))\n        plt.title(f'{site}, {floorNo}, {path}, {cost}')\n        plt.imshow(dist)\n        plt.plot(x*16, y*16, '.-',label=f\"before_pp\")\n        plt.plot(x_, y_, '.-',label=f\"after_pp\")\n        plt.plot(dist.shape[1]/2 + np.cumsum(np.append(0,delta_xy_hat[:,0])*16), dist.shape[0]/2+ np.cumsum(np.append(0,delta_xy_hat[:,1])*16), '.-',label=\"delta\")\n        plt.legend(bbox_to_anchor=(1.05, 1), loc='upper left', borderaxespad=0)\n        plt.show()\n    \n    return pd.DataFrame({\n        'site_path_timestamp' : path_df['site_path_timestamp'],\n        'floor' : path_df['floor'],\n        'x' : x_/16,\n        'y' : y_/16\n    })","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"To compare delta and prediction, the shape of delta is shown in the center of the image.","metadata":{}},{"cell_type":"code","source":"processes = multiprocessing.cpu_count()\nwith multiprocessing.Pool(processes = processes) as pool:\n    dfs = pool.imap_unordered(correct_path, ss.groupby('path'))\n    dfs = tqdm(dfs)\n    dfs = list(dfs)\nss = pd.concat(dfs).sort_values('site_path_timestamp')","metadata":{"papermill":{"duration":351.258931,"end_time":"2021-03-19T16:15:02.477024","exception":false,"start_time":"2021-03-19T16:09:11.218093","status":"completed"},"tags":[],"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Submission","metadata":{"papermill":{"duration":0.071204,"end_time":"2021-03-19T16:15:10.517213","exception":false,"start_time":"2021-03-19T16:15:10.446009","status":"completed"},"tags":[]}},{"cell_type":"code","source":"ss.to_csv('submission.csv', index = None)","metadata":{"papermill":{"duration":0.155056,"end_time":"2021-03-19T16:15:10.742762","exception":false,"start_time":"2021-03-19T16:15:10.587706","status":"completed"},"tags":[],"trusted":true},"execution_count":null,"outputs":[]}]}