{"cells":[{"metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","trusted":true},"cell_type":"code","source":"from datetime import datetime\nimport os\n\nimport numpy as np # linear algebra\nimport pandas as pd # data processing, CSV file I/O (e.g. pd.read_csv)\nfrom scipy.stats import skew\n\nfrom IPython.core.display import display\nfrom tqdm import tqdm\ntqdm.pandas()\n\nprint(os.listdir(\"../input\"))","execution_count":null,"outputs":[]},{"metadata":{"_cell_guid":"79c7e3d0-c299-4dcb-8224-4455121ee9b0","collapsed":true,"_uuid":"d629ff2d2480ee46fbb7e2d37f6b5fab8052498a","trusted":false},"cell_type":"markdown","source":"pivot_tableを用いた高速データ処理\n==========================\n\n# TL;DR\nカテゴリごとの集計をしたいときにpivot_tableを用いると直感的で高速な処理ができる"},{"metadata":{"trusted":true,"_uuid":"f39cf48a80a00768a46754f2a894e5f90148e841"},"cell_type":"code","source":"# 冒頭のpivot_tableを用いたバージョン\n# 欠損値が含まれることに注意する\n\ntick = datetime.now()\ntrain_df = pd.read_csv(\"../input/training_set.csv\", dtype={\"object_id\": np.uint32,\n                                                           \"mjd\": np.float64,\n                                                           \"passband\": np.uint8,\n                                                           \"flux\": np.float32,\n                                                           \"flux_err\": np.float32,\n                                                           \"detected\": np.uint8})\ntrain_meta_df = pd.read_csv('../input/training_set_metadata.csv')\ntock = datetime.now()\nprint(\"load_data: {} ms\".format((tock - tick).seconds * 1000 + ((tock - tick).microseconds / 1000)))\n\ntick = datetime.now()\n\n# pivot_tableのindexをrankを用いて作成する\ntrain_df[\"rank\"] = train_df.groupby([\"object_id\", \"passband\"])[\"mjd\"].rank()\n\nflux = train_df.pivot_table(columns=[\"object_id\", \"passband\"],\n                            index=\"rank\",\n                            values=\"flux\",\n                            aggfunc=\"mean\")\ndflux = train_df.pivot_table(columns=[\"object_id\", \"passband\"],\n                             index=\"rank\",\n                             values=\"flux_err\",\n                             aggfunc=\"mean\")\n\n# 列にNaNが含まれるので扱いに注意する\nflux_mean = np.sum(flux*np.square(flux/dflux), axis=0)/np.sum(np.square(flux/dflux), axis=0)\nflux_std = np.std(flux/flux_mean, ddof = 1, axis=0)\nflux_amp = (np.max(flux, axis=0) - np.min(flux, axis=0))/flux_mean\nflux_mad = np.nanmedian(np.abs((flux - np.nanmedian(flux, axis=0))/flux_mean), axis=0) # array\nflux_beyond = np.sum(np.abs(flux - flux_mean) > np.std(flux, ddof = 1, axis=0), axis=0)/flux.count()\nflux_skew = skew(flux, nan_policy=\"omit\", axis=0)  # masked_array\n\nresult_df = pd.concat([flux_mean.reset_index(name=\"flux_mean\"),\n                      flux_std.reset_index(name=\"flux_std\").iloc[:, 2:],\n                      flux_amp.reset_index(name=\"flux_amp\").iloc[:, 2:],\n                      flux_beyond.reset_index(name=\"flux_beyond\").iloc[:, 2:]], axis=1)\nresult_df[\"flux_mad\"] = flux_mad\nresult_df[\"flux_skew\"] = flux_skew\ncolnames = [\"flux_mean\", \"flux_std\", \"flux_amp\", \"flux_beyond\", \"flux_mad\", \"flux_skew\"]\n\nfor j in range(6):\n    train_meta_df = train_meta_df.merge(result_df.loc[result_df.passband == j, :]\n                                                 .rename(columns={colname: \"{}_{}\".format(colname, j) for colname in colnames})\n                                                 .drop(\"passband\", axis=1),\n                                        how=\"left\",\n                                        on=[\"object_id\"])\ntock = datetime.now()\nprint(\"processing_time: {} sec\".format((tock - tick).seconds))\n\ntrain_meta_df.head()","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"5c30ffebd6eeded97eef71efb563ad09826b3c9b"},"cell_type":"markdown","source":"## ダミーデータでの解説\n### 基本的な演算"},{"metadata":{"trusted":true,"_uuid":"a98725a45f59ee8125fa6c3b2571c1d254dd3317"},"cell_type":"code","source":"# 以下のようなデータを用意する\ndammy_dics = []\nfor i in range(5):\n    for j in range(10):\n        dammy_dics.append({\"time\": i, \"category\": j, \"price\": 10*i + j})\n\ndammy_df = pd.DataFrame(dammy_dics)\ndammy_df.head(10)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"1c0ff52e5487945eb9c52ad51c80c682638bf132"},"cell_type":"code","source":"# DataFrame.pivot_table()でクロス集計表を作れる\ndammy_piv = dammy_df.pivot_table(index=\"time\",\n                                 columns=\"category\",\n                                 values=\"price\",\n                                 aggfunc=\"sum\")\ndisplay(dammy_piv)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"29b466a72f59f006c5f359727cfe215df30a471d"},"cell_type":"code","source":"# pivot_tableは行列として計算することができる\n# 各数値を二乗する\nprint(\"piv^2\")\ndisplay(np.square(dammy_piv))\n\n# スカラーで割る\"\nprint(\"piv / 10\")\ndisplay(dammy_piv / 10)\n\n# pivot_table同士を足す\nprint(\"piv + piv^2\")\ndisplay(dammy_piv + np.square(dammy_piv))","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"5e41ed83e5e3730c107de077ffcf5d40a3fa34cc"},"cell_type":"code","source":"# 列方向への集計\n# axisを指定しないと自動的に列方向の集計になり、Seriesが返ってくる\ndisplay(dammy_piv.mean())\n\n# pivot_tableに対してSeriesで計算するとと列方向にbroadcastされる\ndisplay(dammy_piv - dammy_piv.mean())","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"a2575b10efc92320391252dd08c4e70b186ffb55"},"cell_type":"code","source":"# \"行方向への集計も可能だが\"\ndisplay(dammy_piv.mean(axis=1))\n\n# いい感じにbroadcastしてくれない\nprint(\"piv - seires\")\ndisplay(dammy_piv - dammy_piv.mean(axis=1))\n\n# 転値を使うくらいしか良い方法が思い浮かばないので良い方法があれば教えてください\nprint(\"(piv.T - series).T\")\ndisplay((dammy_piv.T - dammy_piv.mean(axis=1)).T)","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"e3aa7310675c356266f125ef31ad46aa506e75fb"},"cell_type":"markdown","source":"### 差分"},{"metadata":{"trusted":true,"_uuid":"1243fdccff3f61cf86c98dd5f28c3746477213f5"},"cell_type":"code","source":"# piv.shift()でひとつ前の値をとれる\ndammy_piv.shift(1)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"c41fc69d6e4a76d8ca181c5cd7c3dfd19a88f345"},"cell_type":"code","source":"# これを活用すると、ひとつ前との差分をとることができる\ndammy_piv - dammy_piv.shift(1)","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"ecf275cd24487aab093034ceee2e1ed662d8b25c"},"cell_type":"markdown","source":"### 移動平均"},{"metadata":{"trusted":true,"_uuid":"f9fdfec66930fa3770df4a650698ad165455c987"},"cell_type":"code","source":"# rolling関数で、移動平均等をとることができる\n# 以下のコードは自信を含めた三つの期間分の平均\ndammy_piv.rolling(window=3, center=False).mean()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"4f993eef4f1c0f7304a6838a1de58fad5fd3c8a5"},"cell_type":"code","source":"# shiftと組み合わせることで、一つ前からn個前までの平均といった特徴量を作ることができる\ndammy_piv.rolling(window=3, center=False).mean().shift(1)","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"6138efbd0a40799f87a8ece3b5cf34b0c035dde1"},"cell_type":"markdown","source":"### ある時点までの合計を計算する"},{"metadata":{"trusted":true,"_uuid":"6ca7cee958a3fb8e904d8baf8a0a2d53a573414d"},"cell_type":"code","source":"# cum〇〇系の関数はそれまでの合計を計算できる\n# 合計\ndisplay(dammy_piv.cumsum())","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"90bfa38f17dc9ede0d56b1d303d8a93e7f2f91f0"},"cell_type":"code","source":"# 上記までのテクニックを駆使すると、leak無しに時系列のmean_encodingができる\ncum_sum = dammy_df.pivot_table(index=\"time\",\n                               columns=\"category\",\n                               values=\"price\",\n                               aggfunc=\"sum\").cumsum()\ncum_count = dammy_df.pivot_table(index=\"time\",\n                                 columns=\"category\",\n                                 values=\"price\",\n                                 aggfunc=\"count\").cumsum()\ncum_mean = cum_sum / cum_count\ncum_mean_without_leakage = cum_mean.shift(1)\ncum_mean_without_leakage","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"663d4e45988e827fbcb00dc0755785d473a36fca"},"cell_type":"markdown","source":"# PLAsTiCCのデータを用いた実例"},{"metadata":{"_uuid":"1db6c34ba154f8d94b659388139f2224c4e9e87c"},"cell_type":"markdown","source":"[Starter Kit](http://www.kaggle.com/michaelapers/the-plasticc-astronomy-starter-kit)  の3章のLightCurve内にある特徴量を計算する。"},{"metadata":{"trusted":true,"_uuid":"c0264cdf385f7c93bb1d0af72c04b48cc89d516d"},"cell_type":"code","source":"# データのロード\ntrain_df = pd.read_csv(\"../input/training_set.csv\", dtype={\"object_id\": np.uint32,\n                                                           \"mjd\": np.float64,\n                                                           \"passband\": np.uint8,\n                                                           \"flux\": np.float32,\n                                                           \"flux_err\": np.float32,\n                                                           \"detected\": np.uint8})\ntrain_meta_df = pd.read_csv('../input/training_set_metadata.csv')\ntest_meta_df = pd.read_csv('../input/test_set_metadata.csv')","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"3e0ce8c96c21356fa3ebae8af6fd5a87422b19c9"},"cell_type":"code","source":"# train_dfを集計してtrain_metaに結合したい\ndisplay(train_df.head())\ndisplay(train_meta_df.head())","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"18f3255e59c5258419bcb2fbad6a59f26a2bbdc0"},"cell_type":"code","source":"print(\"train_meta: \", train_meta_df.shape)\nprint(\"test_meta: \", test_meta_df.shape)\nprint(\"テストデータは訓練データの{:.4}倍\".format(test_meta_df.shape[0] / train_meta_df.shape[0]))","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"196d463206e1ca35b1295bef7a94efed098d04b2"},"cell_type":"markdown","source":"## 愚直に計算する"},{"metadata":{"trusted":true,"_uuid":"5f50284566f46f535111c7995ef1df7a336b76ba"},"cell_type":"code","source":"# groupby無しに毎回取り出そうとするととてつもない時間がかかるので1/100だけ計算\nbands = [train_df.passband == b for b in train_df.passband.unique()]\nfor id_ in tqdm(train_df.object_id.unique()[:78]):\n    for band in bands:\n        idx = train_df[(train_df.object_id == id_) & band].index\n        flux, dflux = train_df.loc[idx, \"flux\"], train_df.loc[idx, \"flux_err\"]\n        train_df.loc[idx, \"flux_mean\"] = np.sum(flux*np.square(flux/dflux))/np.sum(np.square(flux/dflux))\n        fluxm = train_df.loc[idx, \"flux_mean\"]\n\n        train_df.loc[idx, \"flux_std\"] = np.std(flux/fluxm, ddof = 1)\n        train_df.loc[idx, \"flux_amp\"] = (np.max(flux) - np.min(flux))/fluxm\n        train_df.loc[idx, \"flux_mad\"] = np.median(np.abs((flux - np.median(flux))/fluxm))\n        train_df.loc[idx, \"flux_beyond\"] = sum(np.abs(flux - fluxm) > np.std(flux, ddof = 1))/len(flux)\n        train_df.loc[idx, \"flux_skew\"] = skew(flux)","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"4c0252dfa0d35e5ff69113e27c6908942b249b40"},"cell_type":"markdown","source":"これだとtrain_dataの処理でも一時間以上かかるので、その450倍もあるtest_dataを処理することはできない\n\n## groupbyを使って計算する"},{"metadata":{"trusted":true,"_uuid":"ee84a1830c1d0e9c6f5fea026b16388337144fd5"},"cell_type":"code","source":"# 2. groupbyを使って計算する\ntick = datetime.now()\ntrain_df = pd.read_csv(\"../input/training_set.csv\", dtype={\"object_id\": np.uint32,\n                                                           \"mjd\": np.float64,\n                                                           \"passband\": np.uint8,\n                                                           \"flux\": np.float32,\n                                                           \"flux_err\": np.float32,\n                                                           \"detected\": np.uint8})\ntrain_meta_df = pd.read_csv('../input/training_set_metadata.csv')\ntock = datetime.now()\nprint(\"load_data: {} ms\".format((tock - tick).seconds * 1000 + ((tock - tick).microseconds / 1000)))\n\ntick = datetime.now()\n\ndef agg_func(x):\n    d = {}\n    flux, dflux = x[\"flux\"], x[\"flux_err\"]\n    flux_mean = np.sum(flux*np.square(flux/dflux))/np.sum(np.square(flux/dflux))\n    d[\"flux_mean\"] = flux_mean\n    d[\"flux_std\"] = np.std(flux/flux_mean, ddof = 1)\n    d[\"flux_amp\"] = (np.max(flux) - np.min(flux))/flux_mean\n    d[\"flux_beyond\"] = np.sum(np.abs(flux - flux_mean) > np.std(flux, ddof = 1))/flux.shape[0]\n    d[\"flux_mad\"] = np.median(np.abs((flux - np.median(flux))/flux_mean))\n    d[\"flux_skew\"] = skew(flux)\n    return pd.Series(d, index = [\"flux_mean\", \"flux_std\", \"flux_amp\", \"flux_mad\", \"flux_beyond\", \"flux_skew\"])\n\nresult_df = train_df.groupby([\"object_id\", \"passband\"]).progress_apply(agg_func).reset_index()\n\ncolnames = [\"flux_mean\", \"flux_std\", \"flux_amp\", \"flux_mad\", \"flux_beyond\", \"flux_skew\"]\nfor j in range(6):\n    train_meta_df = train_meta_df.merge(result_df.loc[result_df.passband == j, :]\n                                                 .rename(columns={colname: \"{}_{}\".format(colname, j) for colname in colnames})\n                                                 .drop(\"passband\", axis=1),\n                                        how=\"left\",\n                                        on=[\"object_id\"])\n\ntock = datetime.now()\ntmp = print(\"total_processing: {} sec\".format((tock - tick).seconds))\ntrain_meta_df.head()","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"d207265beb7246ca7298cb2a16bd0442c9de3203"},"cell_type":"markdown","source":"一時間以上かかった処理を二分半で終えることができたが、testデータだと900分 = 15時間かかるのでまだまだ高速化したい\n\n## pivot_tableを使って計算する"},{"metadata":{"trusted":true,"_uuid":"0e1078f14e4e268c08b27e81109c629e86ddeca6"},"cell_type":"code","source":"# 欠損値が含まれることに注意する\n\ntick = datetime.now()\ntrain_df = pd.read_csv(\"../input/training_set.csv\", dtype={\"object_id\": np.uint32,\n                                                           \"mjd\": np.float64,\n                                                           \"passband\": np.uint8,\n                                                           \"flux\": np.float32,\n                                                           \"flux_err\": np.float32,\n                                                           \"detected\": np.uint8})\ntrain_meta_df = pd.read_csv('../input/training_set_metadata.csv')\ntock = datetime.now()\nprint(\"load_data: {} ms\".format((tock - tick).seconds * 1000 + ((tock - tick).microseconds / 1000)))\n\ntick = datetime.now()\n\n# pivot_tableのindexをrankを用いて作成する\ntrain_df[\"rank\"] = train_df.groupby([\"object_id\", \"passband\"])[\"mjd\"].rank()\n\nflux = train_df.pivot_table(columns=[\"object_id\", \"passband\"],\n                            index=\"rank\",\n                            values=\"flux\",\n                            aggfunc=\"mean\")\ndflux = train_df.pivot_table(columns=[\"object_id\", \"passband\"],\n                             index=\"rank\",\n                             values=\"flux_err\",\n                             aggfunc=\"mean\")\n\n# 列にNaNが含まれるので扱いに注意する\nflux_mean = np.sum(flux*np.square(flux/dflux), axis=0)/np.sum(np.square(flux/dflux), axis=0)\nflux_std = np.std(flux/flux_mean, ddof = 1, axis=0)\nflux_amp = (np.max(flux, axis=0) - np.min(flux, axis=0))/flux_mean\nflux_mad = np.nanmedian(np.abs((flux - np.nanmedian(flux, axis=0))/flux_mean), axis=0) # array\nflux_beyond = np.sum(np.abs(flux - flux_mean) > np.std(flux, ddof = 1, axis=0), axis=0)/flux.count()\nflux_skew = skew(flux, nan_policy=\"omit\", axis=0)  # masked_array\n\nresult_df = pd.concat([flux_mean.reset_index(name=\"flux_mean\"),\n                      flux_std.reset_index(name=\"flux_std\").iloc[:, 2:],\n                      flux_amp.reset_index(name=\"flux_amp\").iloc[:, 2:],\n                      flux_beyond.reset_index(name=\"flux_beyond\").iloc[:, 2:]], axis=1)\nresult_df[\"flux_mad\"] = flux_mad\nresult_df[\"flux_skew\"] = flux_skew\ncolnames = [\"flux_mean\", \"flux_std\", \"flux_amp\", \"flux_beyond\", \"flux_mad\", \"flux_skew\"]\n\nfor j in range(6):\n    train_meta_df = train_meta_df.merge(result_df.loc[result_df.passband == j, :]\n                                                 .rename(columns={colname: \"{}_{}\".format(colname, j) for colname in colnames})\n                                                 .drop(\"passband\", axis=1),\n                                        how=\"left\",\n                                        on=[\"object_id\"])\ntock = datetime.now()\nprint(\"processing_time: {} sec\".format((tock - tick).seconds))\n\ntrain_meta_df.head()","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"d4f6144c6cdcfd53953cac2fffd6439dc7982d5e"},"cell_type":"markdown","source":"groupbyで二分半ほどかかっていた処理を、わずか4秒で処理することができた!!\n\ntestデータは大きすぎるので一度に計算しようとするとメモリに乗り切らないが、私の環境(RAM 32GB)だと10分割して計算しおおよそ30分くらいで処理が終わった。\n"}],"metadata":{"kernelspec":{"display_name":"Python 3","language":"python","name":"python3"},"language_info":{"name":"python","version":"3.6.6","mimetype":"text/x-python","codemirror_mode":{"name":"ipython","version":3},"pygments_lexer":"ipython3","nbconvert_exporter":"python","file_extension":".py"}},"nbformat":4,"nbformat_minor":1}