{"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":"# はじめてのG2N\n\nG2Nをトライするにあたり、知っていないといけないことが多そうなのでここにつらつらとまとめる。","metadata":{}},{"cell_type":"markdown","source":"# ライブラリ読み込み","metadata":{}},{"cell_type":"code","source":"!pip -qq install git+https://github.com/PyFstat/PyFstat@python37\n!pip install dictdiffer\n##== h5webはうまく機能せず。\n#!pip install jupyterlab_h5web\n \nimport os\nimport h5py # import to read hdf5\nimport numpy as np\nimport pandas as pd\nfrom glob import glob\nfrom pathlib import Path\nimport matplotlib.pyplot as plt\nfrom datetime import datetime\n%matplotlib inline\n\nimport pyfstat\nfrom pyfstat.utils import get_sft_as_arrays\n\n##== hdf5 h5webはうまく機能せず。\n#from jupyterlab_h5web import H5Web\n#!jupyter serverextension enable jupyterlab_h5web\n\n##== dict diff\nfrom dictdiffer import diff ","metadata":{"_kg_hide-input":true,"_kg_hide-output":true,"execution":{"iopub.status.busy":"2022-11-27T07:26:09.109201Z","iopub.execute_input":"2022-11-27T07:26:09.110519Z","iopub.status.idle":"2022-11-27T07:26:36.573011Z","shell.execute_reply.started":"2022-11-27T07:26:09.110464Z","shell.execute_reply":"2022-11-27T07:26:36.571175Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Utils","metadata":{}},{"cell_type":"code","source":"def h5_struct(h5,l=0,struct=dict()):\n l += 1\n if type(h5) == h5py._hl.dataset.Dataset :\n print(\" \"*l + f\"Dataset shape : {h5.shape}\")\n else:\n for k in h5.keys():\n print(\" \"*l + k)\n struct[k] = dict()\n h5_struct(h5[k],l,struct[k])\n return struct\n\ndef h5_struct2dict(h5,l=0,struct=dict()):\n l += 1\n if type(h5) == h5py._hl.dataset.Dataset :\n struct = h5.shape\n return struct\n else:\n for k in h5.keys():\n struct[k] = dict()\n struct[k] = h5_struct2dict(h5[k],l,struct[k])\n return struct\n\ndef h5_show_freq_time(h5,l=0,key=str()):\n l += 1\n if type(h5) == h5py._hl.dataset.Dataset :\n if key == \"timestamps_GPS\":\n print(\" \"*l + \"head data : \")\n print([datetime.utcfromtimestamp(i).strftime('%Y-%m-%d %H:%M:%S') for i in h5[()][:5]])\n print(\" \"*l + \"tail data : \")\n print([datetime.utcfromtimestamp(i).strftime('%Y-%m-%d %H:%M:%S') for i in h5[()][-5:]])\n elif key == \"frequency_Hz\":\n print(\" \"*l + f\"head data : {h5[()][:5]}\")\n print(\" \"*l + f\"tail data : {h5[()][-5:]}\")\n else:\n pass\n else:\n for k in h5.keys():\n print(\" \"*l + k)\n h5_show_freq_time(h5[k],l,key=k)\n\ndef h5_show_dataset(h5,l=0,key=str()):\n l += 1\n if type(h5) == h5py._hl.dataset.Dataset :\n if key == \"timestamps_GPS\":\n print(\" \"*l + \"head data : \")\n print([datetime.utcfromtimestamp(i).strftime('%Y-%m-%d %H:%M:%S') for i in h5[()][:5]])\n print(\" \"*l + \"tail data : \")\n print([datetime.utcfromtimestamp(i).strftime('%Y-%m-%d %H:%M:%S') for i in h5[()][-5:]])\n else:\n print(\" \"*l + f\"head data : {h5[()][:5]}\")\n print(\" \"*l + f\"tail data : {h5[()][-5:]}\")\n else:\n for k in h5.keys():\n print(\" \"*l + k)\n h5_show_dataset(h5[k],l,key=k)","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2022-11-27T07:36:03.813276Z","iopub.execute_input":"2022-11-27T07:36:03.814565Z","iopub.status.idle":"2022-11-27T07:36:03.833129Z","shell.execute_reply.started":"2022-11-27T07:36:03.814509Z","shell.execute_reply":"2022-11-27T07:36:03.831468Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### 2021年にもG2Netのコンペがあった\n\n2021年開催のコンペ \n参考になるはず \nhttps://www.kaggle.com/c/g2net-gravitational-wave-detection ","metadata":{}},{"cell_type":"markdown","source":"### これは読んでおいたほうがよさそう\n\n* SFTsについて \n [前回コンペのディスカッションから](https://www.kaggle.com/competitions/g2net-detecting-continuous-gravitational-waves/discussion/358444)","metadata":{}},{"cell_type":"markdown","source":"### 重力波ってそもそもなによ\n \n質量を持った物体が加速度運動することで放射される波。 \n観測できるほどの大きな振幅の重力波を発生させるには、高密度で非常に大きな質量の物体が加速度運動する必要。 \n\n重力波の発生源としては以下のような天体運動、天体現象が挙げられる。 \n\n\n\n[国立天文台　重力波プロジェクト推進室](https://gwpo.nao.ac.jp/about_gw)　から\n \n \n","metadata":{}},{"cell_type":"markdown","source":"### どうやって重力波を観測？\n\nレーザー干渉計型検出器を使って検出\n\n\n\n[国立天文台　重力波プロジェクト推進室](https://gwpo.nao.ac.jp/about_gw) から","metadata":{}},{"cell_type":"markdown","source":"### 観測しているのはどこ？\n\n* LIGO「レーザー干渉計重力波観測所」 \n [天文学辞典](https://astro-dic.jp/laser-interferometer-gravitational-wave-observatory/) \n [wikipedia](https://ja.wikipedia.org/wiki/LIGO#:~:text=LIGO%EF%BC%88%E3%83%A9%E3%82%A4%E3%82%B4%E3%80%81%E8%8B%B1%E8%AA%9E%3A%20Laser,%E6%B3%A2%E8%A6%B3%E6%B8%AC%E6%89%80%E3%80%8D%E3%81%A8%E3%81%AA%E3%82%8B%E3%80%82) \n [【Python】自力で重力波解析](http://yasketballclub.work/2020/02/24/%E3%80%90python%E3%80%91%E8%87%AA%E5%8A%9B%E3%81%A7%E9%87%8D%E5%8A%9B%E6%B3%A2%E8%A7%A3%E6%9E%90/) \n \n* Virgo \n [天文学辞典](https://astro-dic.jp/virgo-interferometer/) \n [Wikipedia](https://ja.wikipedia.org/wiki/Virgo#:~:text=Virgo%E3%81%AF10%20Hz%E3%81%8B%E3%82%89,%E6%B3%A2%E3%81%AB%E6%84%9F%E5%BA%A6%E3%82%92%E6%8C%81%E3%81%A4%E3%80%82) \nなどがある","metadata":{}},{"cell_type":"markdown","source":"### コンペデータと観測所\n\n今回のコンペで使用するデータは、LIGOの二箇所の観測所LIGO Hanford & LIGO Livingstonのデータ。 \n２つの施設は3002 km 離れており、光速度で伝播する重力波の到達時間として約10ミリ秒の差がある。 \n波源からの2つの施設への重力波の到達時間の違いから、三角測量を応用して波源の位置を知ることができる。\n\n[wikipedia](https://ja.wikipedia.org/wiki/LIGO#:~:text=LIGO%EF%BC%88%E3%83%A9%E3%82%A4%E3%82%B4%E3%80%81%E8%8B%B1%E8%AA%9E%3A%20Laser,%E6%B3%A2%E8%A6%B3%E6%B8%AC%E6%89%80%E3%80%8D%E3%81%A8%E3%81%AA%E3%82%8B%E3%80%82) から \n\n\n","metadata":{}},{"cell_type":"markdown","source":"### どれくらい検出されている？\n\n年に数回程度らしい。 \n[wikipedia](https://ja.wikipedia.org/wiki/%E9%87%8D%E5%8A%9B%E6%B3%A2%E3%81%AE%E5%88%9D%E6%A4%9C%E5%87%BA) から \n\n観測されたデータは以下のネーミングルールで名前が振られている。 \nGW150914 (重力波 Gravitational Waveの頭文字と、観測された日付2015年09月14日) \n[wikipedia](https://ja.wikipedia.org/wiki/%E9%87%8D%E5%8A%9B%E6%B3%A2%E3%81%AE%E5%88%9D%E6%A4%9C%E5%87%BA) から \nコンペデータのIDに使われているかと思ったが、そうではなかった。残念。","metadata":{}},{"cell_type":"markdown","source":"### コンペデータのフォーマットHDF5\n\nよく聞くけど使うのは初めて。 \nHierarchical Data Formatsの略。 \n \n[Hierarchical Data Formats - What is HDF5?](https://www.neonscience.org/resources/learning-hub/tutorials/about-hdf5) \n\npythonでHDF5を扱うときの説明はこちら。 \n[意外と奥が深い、HDFの世界（Python・h5py入門）](https://qiita.com/simonritchie/items/23db8b4cb5c590924d95) \n\n[け日記](https://ohke.hateblo.jp/entry/2020/05/16/230000)\n\n\n","metadata":{}},{"cell_type":"markdown","source":"### HDF5フォーマットに触れてみる\n\nこちらのコードを参考にしながら \nhttps://www.kaggle.com/code/ayuraj/g2net-understand-the-data\n\n","metadata":{}},{"cell_type":"code","source":"train_files = glob(\"/kaggle/input/g2net-detecting-continuous-gravitational-waves/train/*.hdf5\")\nfile = Path(train_files[0])\nh5 = h5py.File(file, 'r')\nh5_struct(h5)","metadata":{"execution":{"iopub.status.busy":"2022-11-27T07:47:43.504408Z","iopub.execute_input":"2022-11-27T07:47:43.505178Z","iopub.status.idle":"2022-11-27T07:47:43.537387Z","shell.execute_reply.started":"2022-11-27T07:47:43.505134Z","shell.execute_reply":"2022-11-27T07:47:43.536054Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### コメント\n\nH1とL1でHDF5データの構造が違う。。。 \nもしかしてデータの中に違う構造があったりする？？ ","metadata":{}},{"cell_type":"code","source":"train_data_struct = dict()\nfor f in train_files:\n h5f = Path(f)\n h5 = h5py.File(h5f, 'r')\n h5_struct2dict(h5,0,train_data_struct)","metadata":{"execution":{"iopub.status.busy":"2022-11-27T07:08:07.024764Z","iopub.execute_input":"2022-11-27T07:08:07.025244Z","iopub.status.idle":"2022-11-27T07:08:14.634082Z","shell.execute_reply.started":"2022-11-27T07:08:07.025203Z","shell.execute_reply":"2022-11-27T07:08:14.63256Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### コメント\n\nすべて構造は同じ。以下のCellを実行することで確認可能。 \nデータのサイズはどれも違う。 \nfrequency_Hzはどれも同じ様子。 ","metadata":{}},{"cell_type":"code","source":"kamisama_dict_struct = dict()\nfor k,v in train_data_struct.items():\n if len(kamisama_dict_struct) == 0:\n kamisama_dict_struct = v\n else:\n print(list(diff(kamisama_dict_struct,v)))","metadata":{"_kg_hide-output":true,"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### 時刻と周波数は全部同じ？\n\nどうもバラバラ。全ファイルの時刻と周波数を、前後5つずつ確認","metadata":{}},{"cell_type":"code","source":"for f in train_files[:10]:\n h5 = h5py.File(f, 'r')\n h5_show_freq_time(h5)","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### 同じ開始時刻のデータがあるか知りたくなる。\n\n全データから、開始時刻が同じデータが複数あるか確認。 \n結果、ばんらばら。 \n興味があれば次を実行。 ","metadata":{}},{"cell_type":"code","source":"h1_timestanp_dict = dict()\nl1_timestanp_dict = dict()\nfor f in train_files:\n h5 = h5py.File(f, 'r')\n for k in h5.keys():\n if str(h5[k][\"H1\"][\"timestamps_GPS\"][0]) not in h1_timestanp_dict:\n h1_timestanp_dict[f't{str(h5[k][\"H1\"][\"timestamps_GPS\"][0])}'] = list()\n \n if str(h5[k][\"L1\"][\"timestamps_GPS\"][0]) not in h1_timestanp_dict:\n l1_timestanp_dict[f't{str(h5[k][\"L1\"][\"timestamps_GPS\"][0])}'] = list()\n \n h1_timestanp_dict[f't{str(h5[k][\"H1\"][\"timestamps_GPS\"][0])}'].append(k)\n l1_timestanp_dict[f't{str(h5[k][\"L1\"][\"timestamps_GPS\"][0])}'].append(k)\n \nprint(l1_timestanp_dict)","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### SFTsって何よ？\n\nデータセットの説明にあるように、 \nShort-time Fourier Transforms (SFTs)のこと。\n \n各データセットの前5つ、後ろ5をチラ見。 \nSFTsには実数成分、虚数成分がある。 \n","metadata":{}},{"cell_type":"code","source":"file = Path(train_files[0])\nh5 = h5py.File(file, 'r')\nh5_show_dataset(h5)","metadata":{"execution":{"iopub.status.busy":"2022-11-27T07:26:36.576851Z","iopub.execute_input":"2022-11-27T07:26:36.57741Z","iopub.status.idle":"2022-11-27T07:26:36.656249Z","shell.execute_reply.started":"2022-11-27T07:26:36.577355Z","shell.execute_reply":"2022-11-27T07:26:36.654209Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"","metadata":{},"execution_count":null,"outputs":[]}]}