{"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":"!pip install nb-black","metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","_kg_hide-output":true,"execution":{"iopub.status.busy":"2022-06-30T06:19:23.203717Z","iopub.execute_input":"2022-06-30T06:19:23.204104Z","iopub.status.idle":"2022-06-30T06:19:37.964095Z","shell.execute_reply.started":"2022-06-30T06:19:23.204021Z","shell.execute_reply":"2022-06-30T06:19:37.963025Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"from functools import partial\n\nimport pandas as pd\nimport numpy as np\nimport matplotlib.pyplot as plt\nimport seaborn as sns\n\nimport plotly.express as px\nimport plotly.offline as py\nimport plotly.graph_objects as go\nimport pyarrow.parquet as pq\nimport pyarrow as pa\n\nfrom numpy import sin, cos, deg2rad, rad2deg\nfrom plotly.offline import init_notebook_mode, iplot\nfrom sklearn.metrics.pairwise import haversine_distances\nfrom sklearn.model_selection import train_test_split\nfrom sklearn.neighbors import KNeighborsRegressor\nfrom tqdm.auto import tqdm\n\npd.set_option(\"max_colwidth\", 900)\nplt.style.use(\"ggplot\")\ninit_notebook_mode(connected=True)\n\n%load_ext lab_black","metadata":{"execution":{"iopub.status.busy":"2022-06-30T06:19:37.965852Z","iopub.execute_input":"2022-06-30T06:19:37.966435Z","iopub.status.idle":"2022-06-30T06:19:40.08849Z","shell.execute_reply.started":"2022-06-30T06:19:37.966409Z","shell.execute_reply":"2022-06-30T06:19:40.087521Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def vectorized_haversine(X):\n    lat1, lon1, lat2, lon2 = X\n    radius_km = 6371\n    dlat = lat2 - lat1\n    dlon = lon2 - lon1\n    a = (np.sin(dlat / 2) ** 2) + np.cos(lat1) * np.cos(lat2) * (np.sin(dlon / 2) ** 2)\n    c = 2 * np.arctan2(np.sqrt(a), np.sqrt(1 - a))\n    d_km = radius_km * c\n    return d_km","metadata":{"execution":{"iopub.status.busy":"2022-06-30T06:19:40.089863Z","iopub.execute_input":"2022-06-30T06:19:40.090197Z","iopub.status.idle":"2022-06-30T06:19:40.109996Z","shell.execute_reply.started":"2022-06-30T06:19:40.090161Z","shell.execute_reply":"2022-06-30T06:19:40.109018Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def make_poi_cluster(df, key=\"point_of_interest\"):\n    clusters = (\n        df.groupby(key)\n        .agg(\n            count=(\"id\", \"count\"),\n            lat_min=(\"latitude\", \"min\"),\n            lon_min=(\"longitude\", \"min\"),\n            lat_max=(\"latitude\", \"max\"),\n            lon_max=(\"longitude\", \"max\"),\n        )\n        .reset_index()\n    )\n    clusters[[\"lat_min_rad\", \"lon_min_rad\", \"lat_max_rad\", \"lon_max_rad\"]] = clusters[\n        [\"lat_min\", \"lon_min\", \"lat_max\", \"lon_max\"]\n    ].apply(np.deg2rad)\n\n    X = clusters[\n        [\"lat_min_rad\", \"lon_min_rad\", \"lat_max_rad\", \"lon_max_rad\"]\n    ].to_numpy()\n    clusters[\"size_km\"] = vectorized_haversine(X.T)\n    clusters[\"size_m\"] = clusters[\"size_km\"] * 1000.0\n\n    return clusters","metadata":{"execution":{"iopub.status.busy":"2022-06-30T06:19:40.112774Z","iopub.execute_input":"2022-06-30T06:19:40.113127Z","iopub.status.idle":"2022-06-30T06:19:40.136639Z","shell.execute_reply.started":"2022-06-30T06:19:40.113106Z","shell.execute_reply":"2022-06-30T06:19:40.135185Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def make_h3_cluster(df, key=\"h3_res4\"):\n    clusters = (\n        df.groupby(key)\n        .agg(\n            count=(\"id\", \"count\"),\n            max_poi_size_m=(\"poi_size_m\", \"max\"),\n            p99_poi_size_m=(\"poi_size_m\", partial(np.quantile, q=0.95)),\n            p95_poi_size_m=(\"poi_size_m\", partial(np.quantile, q=0.99)),\n            mean_poi_size_m=(\"poi_size_m\", \"mean\"),\n        )\n        .reset_index()\n    )\n\n    return clusters","metadata":{"execution":{"iopub.status.busy":"2022-06-30T06:19:40.137881Z","iopub.execute_input":"2022-06-30T06:19:40.138239Z","iopub.status.idle":"2022-06-30T06:19:40.156923Z","shell.execute_reply.started":"2022-06-30T06:19:40.138211Z","shell.execute_reply":"2022-06-30T06:19:40.15559Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def add_poi_size(df, clusters):\n    df = df.copy()\n    clusters = clusters.copy()\n    clusters = clusters[[\"point_of_interest\", \"size_m\"]]\n    clusters = clusters.rename({\"size_m\": \"poi_size_m\"}, axis=1)\n    df = df.merge(clusters, on=\"point_of_interest\")\n    return df","metadata":{"execution":{"iopub.status.busy":"2022-06-30T06:19:40.158218Z","iopub.execute_input":"2022-06-30T06:19:40.158722Z","iopub.status.idle":"2022-06-30T06:19:40.173557Z","shell.execute_reply.started":"2022-06-30T06:19:40.158673Z","shell.execute_reply":"2022-06-30T06:19:40.172131Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def joint_plot(df, ycol):\n    g = sns.jointplot(\n        data=df,\n        y=ycol,\n        x=\"count\",\n        marginal_kws=dict(bins=20, fill=True),\n    )\n    plt.show()\n    x = df[ycol]\n    y = df[\"count\"]\n    print(f\"correlation: {np.corrcoef(x, y)[0, 1]}\")","metadata":{"execution":{"iopub.status.busy":"2022-06-30T06:19:40.175249Z","iopub.execute_input":"2022-06-30T06:19:40.176354Z","iopub.status.idle":"2022-06-30T06:19:40.369007Z","shell.execute_reply.started":"2022-06-30T06:19:40.176322Z","shell.execute_reply":"2022-06-30T06:19:40.367844Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"train = pd.read_csv(\"../input/foursquare-location-matching/train.csv\")\nh3_df = pq.read_table(\"../input/4sq-h3/h3.parquet\").to_pandas()\ntrain = train.merge(h3_df, on=\"id\")\ndel h3_df","metadata":{"execution":{"iopub.status.busy":"2022-06-30T06:19:40.369998Z","iopub.execute_input":"2022-06-30T06:19:40.370205Z","iopub.status.idle":"2022-06-30T06:19:50.879668Z","shell.execute_reply.started":"2022-06-30T06:19:40.370181Z","shell.execute_reply":"2022-06-30T06:19:50.878663Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"len(train.drop_duplicates(\"name\")) - len(train)","metadata":{"execution":{"iopub.status.busy":"2022-06-30T06:23:55.654457Z","iopub.execute_input":"2022-06-30T06:23:55.655125Z","iopub.status.idle":"2022-06-30T06:23:56.342946Z","shell.execute_reply.started":"2022-06-30T06:23:55.655098Z","shell.execute_reply":"2022-06-30T06:23:56.34238Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"clusters = make_poi_cluster(train)","metadata":{"execution":{"iopub.status.busy":"2022-06-30T06:19:50.88108Z","iopub.execute_input":"2022-06-30T06:19:50.881389Z","iopub.status.idle":"2022-06-30T06:19:54.670177Z","shell.execute_reply.started":"2022-06-30T06:19:50.881362Z","shell.execute_reply":"2022-06-30T06:19:54.669313Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"distinct_clusters = make_poi_cluster(train.drop_duplicates(\"name\"))","metadata":{"execution":{"iopub.status.busy":"2022-06-30T06:24:38.946764Z","iopub.execute_input":"2022-06-30T06:24:38.947064Z","iopub.status.idle":"2022-06-30T06:24:42.492007Z","shell.execute_reply.started":"2022-06-30T06:24:38.947041Z","shell.execute_reply":"2022-06-30T06:24:42.491018Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"train_ext = add_poi_size(train, clusters)","metadata":{"execution":{"iopub.status.busy":"2022-06-30T06:19:54.673021Z","iopub.execute_input":"2022-06-30T06:19:54.673609Z","iopub.status.idle":"2022-06-30T06:19:58.384183Z","shell.execute_reply.started":"2022-06-30T06:19:54.673578Z","shell.execute_reply":"2022-06-30T06:19:58.383207Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"h3_clusters = make_h3_cluster(train_ext)","metadata":{"execution":{"iopub.status.busy":"2022-06-30T06:19:58.385214Z","iopub.execute_input":"2022-06-30T06:19:58.385423Z","iopub.status.idle":"2022-06-30T06:20:02.243386Z","shell.execute_reply.started":"2022-06-30T06:19:58.3854Z","shell.execute_reply":"2022-06-30T06:20:02.242426Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"clusters","metadata":{"execution":{"iopub.status.busy":"2022-06-30T06:20:02.244444Z","iopub.execute_input":"2022-06-30T06:20:02.244672Z","iopub.status.idle":"2022-06-30T06:20:02.27544Z","shell.execute_reply.started":"2022-06-30T06:20:02.244642Z","shell.execute_reply":"2022-06-30T06:20:02.274613Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"h3_clusters","metadata":{"execution":{"iopub.status.busy":"2022-06-30T06:20:02.276615Z","iopub.execute_input":"2022-06-30T06:20:02.277361Z","iopub.status.idle":"2022-06-30T06:20:02.29343Z","shell.execute_reply.started":"2022-06-30T06:20:02.277338Z","shell.execute_reply":"2022-06-30T06:20:02.292517Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# EDA on POIs\n\n* number of clusters\n* cluster size\n* number of locations per POI\n* number of locations v.s. cluster size","metadata":{}},{"cell_type":"code","source":"print(f\"#distinct clusters: {len(clusters):,}\")\nprint(f\"#distinct clusters (n > 1): {len(clusters.query('count > 1')):,}\")\nprint(\n    f\"mean location per distinct POIs: {clusters.query('count > 1')['count'].mean():.3f}\"\n)","metadata":{"execution":{"iopub.status.busy":"2022-06-30T06:37:06.475797Z","iopub.execute_input":"2022-06-30T06:37:06.476096Z","iopub.status.idle":"2022-06-30T06:37:06.573482Z","shell.execute_reply.started":"2022-06-30T06:37:06.476073Z","shell.execute_reply":"2022-06-30T06:37:06.572483Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"fig, ax = plt.subplots(figsize=(16, 4))\nsns.boxenplot(x=\"size_m\", data=clusters, ax=ax)\nax.set(xscale=\"log\", xlim=(1e1, 1e8), title=\"Letter-value plot of cluster size\")\nplt.show()","metadata":{"execution":{"iopub.status.busy":"2022-06-30T06:20:02.295323Z","iopub.execute_input":"2022-06-30T06:20:02.29606Z","iopub.status.idle":"2022-06-30T06:20:03.307183Z","shell.execute_reply.started":"2022-06-30T06:20:02.296022Z","shell.execute_reply":"2022-06-30T06:20:03.306151Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Reference\n\n* letter-value plot[[1]]\n\n[1]: https://www.semanticscholar.org/paper/Letter-Value-Plots%3A-Boxplots-for-Large-Data-Hofmann-Wickham/9f4dfd20d874e16b6288cfc463297ab4a800a032","metadata":{}},{"cell_type":"code","source":"y = [0.5]\nfor i in range(12):\n    y.append(y[i] * 0.5)\ny = np.array(y)\ny = y.cumsum()\nx = [clusters[\"size_m\"].quantile(q=q) for q in y]\n\nfig, ax = plt.subplots()\nax.plot(x, y)\nax.set(\n    xscale=\"log\",\n    xlabel=\"Cluster size [m]\",\n    ylabel=\"Proportion\",\n    title=\"Cumulative plot of size of POIs\",\n)\nplt.show()\npd.DataFrame({\"size_m\": x, \"cum_dist\": y})","metadata":{"execution":{"iopub.status.busy":"2022-06-30T06:37:57.530807Z","iopub.execute_input":"2022-06-30T06:37:57.531803Z","iopub.status.idle":"2022-06-30T06:37:58.04768Z","shell.execute_reply.started":"2022-06-30T06:37:57.53176Z","shell.execute_reply":"2022-06-30T06:37:58.046519Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"clusters[\"size_m\"].max()","metadata":{"execution":{"iopub.status.busy":"2022-06-30T06:37:58.049361Z","iopub.execute_input":"2022-06-30T06:37:58.049623Z","iopub.status.idle":"2022-06-30T06:37:58.059115Z","shell.execute_reply.started":"2022-06-30T06:37:58.049597Z","shell.execute_reply":"2022-06-30T06:37:58.057768Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"* almost all (>99.2%) of POIs have size < 26.4km","metadata":{}},{"cell_type":"code","source":"vc = clusters[\"count\"].value_counts()\nfig, ax = plt.subplots()\nax.plot(vc.cumsum() / vc.sum())\nax.set(\n    xscale=\"log\",\n    xlabel=\"#Locations per POIs\",\n    ylabel=\"Proportion\",\n    title=\"Cumulative plot of #locations per POI\",\n)\nplt.show()\nvc[:10].cumsum() / vc.sum()","metadata":{"execution":{"iopub.status.busy":"2022-06-30T06:37:48.594424Z","iopub.execute_input":"2022-06-30T06:37:48.594813Z","iopub.status.idle":"2022-06-30T06:37:48.90243Z","shell.execute_reply.started":"2022-06-30T06:37:48.594773Z","shell.execute_reply":"2022-06-30T06:37:48.901424Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"* Almost all (>99%) of POIs have locations <= 4.","metadata":{}},{"cell_type":"code","source":"joint_plot(clusters, \"size_m\")","metadata":{"execution":{"iopub.status.busy":"2022-06-25T23:49:17.034281Z","iopub.execute_input":"2022-06-25T23:49:17.035098Z","iopub.status.idle":"2022-06-25T23:49:19.158198Z","shell.execute_reply.started":"2022-06-25T23:49:17.035054Z","shell.execute_reply":"2022-06-25T23:49:19.156905Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Distinct Clusters","metadata":{}},{"cell_type":"code","source":"print(f\"#distinct clusters: {len(distinct_clusters):,}\")\nprint(f\"#distinct clusters (n > 1): {len(distinct_clusters.query('count > 1')):,}\")\nprint(\n    f\"mean location per distinct POIs: {distinct_clusters.query('count > 1')['count'].mean():.3f}\"\n)","metadata":{"execution":{"iopub.status.busy":"2022-06-30T06:37:16.999006Z","iopub.execute_input":"2022-06-30T06:37:16.999387Z","iopub.status.idle":"2022-06-30T06:37:17.061758Z","shell.execute_reply.started":"2022-06-30T06:37:16.999354Z","shell.execute_reply":"2022-06-30T06:37:17.061105Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"vc = distinct_clusters[\"count\"].value_counts()\ndf = pd.DataFrame({\"density\": (vc / vc.sum())})\ndf.index.name = \"#loc per POI\"\n\nfig, ax = plt.subplots()\nax.scatter(df.index, df)\nax.set(\n    xscale=\"log\",\n    xlabel=\"#Locations per POI\",\n    ylabel=\"Proportion\",\n    title=\"Distribution of #locations per POI\",\n)\nplt.show()\ndf.head()","metadata":{"execution":{"iopub.status.busy":"2022-06-30T06:45:55.082467Z","iopub.execute_input":"2022-06-30T06:45:55.082892Z","iopub.status.idle":"2022-06-30T06:45:55.418368Z","shell.execute_reply.started":"2022-06-30T06:45:55.082858Z","shell.execute_reply":"2022-06-30T06:45:55.417151Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# EDA on H3 cluster","metadata":{}},{"cell_type":"code","source":"joint_plot(h3_clusters, \"max_poi_size_m\")","metadata":{"execution":{"iopub.status.busy":"2022-06-25T23:48:29.292613Z","iopub.execute_input":"2022-06-25T23:48:29.29348Z","iopub.status.idle":"2022-06-25T23:48:29.920773Z","shell.execute_reply.started":"2022-06-25T23:48:29.293437Z","shell.execute_reply":"2022-06-25T23:48:29.919424Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"joint_plot(h3_clusters, \"p99_poi_size_m\")","metadata":{"execution":{"iopub.status.busy":"2022-06-25T23:48:33.71484Z","iopub.execute_input":"2022-06-25T23:48:33.71528Z","iopub.status.idle":"2022-06-25T23:48:34.298399Z","shell.execute_reply.started":"2022-06-25T23:48:33.715241Z","shell.execute_reply":"2022-06-25T23:48:34.297167Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"H3セグメント内のポイント数とPOIの最大サイズには少し相関がある~0.33。\nつまり、測位データの外れ値のノイズの大きさはポイントの密度にやや依存してそうである。\n\n逆に、H3セグメント内のポイント数とPOIの99パーセンタイルのサイズには相関がない。","metadata":{}},{"cell_type":"markdown","source":"# EDA on train data\n\n* proportion of filled rows\n* letter-value plot of length of text colums","metadata":{}},{"cell_type":"code","source":"fig, ax = plt.subplots()\ndf = (1 - (train.isna().sum(axis=0) / len(train))).reset_index()\nsns.barplot(data=df, x=0, y=\"index\", ax=ax)\nax.set(title=\"Proportion of filled rows\", xlabel=\"proportion\", ylabel=None)\nplt.show()\ndf","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"text_columns = [\n    \"name\",\n    \"address\",\n    \"city\",\n    \"state\",\n    \"zip\",\n    \"country\",\n    \"url\",\n    \"phone\",\n    \"categories\",\n]","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"train_ext = train.copy()\nfor col in text_columns:\n    train_ext[f\"{col}_len\"] = train_ext[col].str.len()","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"len_columns = [f\"{col}_len\" for col in text_columns]","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"df = train_ext[len_columns].melt()\n\nfig, ax = plt.subplots(figsize=(8, 8))\nsns.boxenplot(x=\"value\", data=df, y=\"variable\", ax=ax)\nax.set(xscale=\"linear\", title=\"Letter-value plot of length of text columns\")\nax.set_yticklabels(len_columns)\nplt.show()","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"y = [0.5]\nfor i in range(10):\n    y.append(y[i] * 0.5)\ny = np.array(y)\ny = y.cumsum()\nx = [train_ext[\"name_len\"].quantile(q=q) for q in y]\n\nfig, ax = plt.subplots()\nax.plot(x, y)\nax.set(\n    xscale=\"linear\",\n    xlabel=\"Text length\",\n    ylabel=\"Proportion\",\n    title=\"Cumulative plot of text length of name\",\n)\nplt.show()\npd.DataFrame({\"name_len\": x, \"cum_dist\": y})","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"train_ext.agg([\"min\", \"median\", \"max\"])","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"train_ext[[\"latitude\", \"longitude\"]].agg([\"min\", \"median\", \"max\"])","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"train_ext[train_ext[\"latitude\"] > 80]","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"train_ext[train_ext[\"latitude\"] <= -80]","metadata":{"trusted":true},"execution_count":null,"outputs":[]}]}