{"metadata":{"kernelspec":{"language":"python","display_name":"Python 3","name":"python3"},"language_info":{"name":"python","version":"3.10.13","mimetype":"text/x-python","codemirror_mode":{"name":"ipython","version":3},"pygments_lexer":"ipython3","nbconvert_exporter":"python","file_extension":".py"},"kaggle":{"accelerator":"none","dataSources":[{"sourceId":56537,"databundleVersionId":8015876,"sourceType":"competition"}],"dockerImageVersionId":30732,"isInternetEnabled":false,"language":"python","sourceType":"notebook","isGpuEnabled":false}},"nbformat_minor":4,"nbformat":4,"cells":[{"cell_type":"markdown","source":"## 1. コンペ概要\n- **タスク**<br>\nマルチスケール気候モデルであるE3SM-MMFの中で、嵐、雲、乱流、降雨、放射などのサブグリッド大気プロセスをエミュレート。\n\n- **モチベーション**<br>\nMLエミュレータはMMFよりもはるかに安価に推論できるため、この面での進歩は、高解像度で物理的に信頼できる長期気候予測が広く利用できるようになり、気候変動に伴う危険性がより明確になり、政策立案者がそれを緩和するために必要な知識を得られるようになる未来の実現に役立つ。\n\n- **評価関数**<br>\n重みづけを用いたカスタム$R^2$スコアを使用する。<br>\n予測を提出する前に、「サンプル提出」と「重み付けファイル」の両方の役割を果たす「sample_submission.csv」にあるデータで、予測データを要素ごとに掛ける必要がある。<br>\n$R^2$スコアは以下のように定義される：<br>\n$$ R^2 = 1 - \\frac{SS_{res}}{SS_{tot}} $$\n$SS_{res}$ は、残差平方和。<br>\n$SS_{tot}$ は、全平方和。つまり、観測値と全観測データの平均値との差の平方和。<br>\n誤差のない予測だとスコア値1。全ての観測値の平均値を予測すると0となり、それより悪い予測だと負の値を取る。<br>\nソースコード：https://www.kaggle.com/code/jerrylin96/r2-score-default","metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19"}},{"cell_type":"markdown","source":"## 2. データ\n学習データは、10,091,520件。<br>\nテストデータ、Sample_submissionは、625,000件。","metadata":{}},{"cell_type":"code","source":"import numpy as np\nimport polars as pl\nimport matplotlib.pyplot as plt\nimport seaborn as sns","metadata":{"execution":{"iopub.status.busy":"2024-06-08T00:17:49.036974Z","iopub.execute_input":"2024-06-08T00:17:49.037348Z","iopub.status.idle":"2024-06-08T00:17:52.69236Z","shell.execute_reply.started":"2024-06-08T00:17:49.037319Z","shell.execute_reply":"2024-06-08T00:17:52.69119Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"INPUT_DIR = '/kaggle/input/leap-atmospheric-physics-ai-climsim'\ntrain = pl.scan_csv(f'{INPUT_DIR}/train.csv')\nsample_submission = pl.scan_csv(f'{INPUT_DIR}/sample_submission.csv')","metadata":{"execution":{"iopub.status.busy":"2024-06-08T00:17:52.694861Z","iopub.execute_input":"2024-06-08T00:17:52.695504Z","iopub.status.idle":"2024-06-08T00:17:52.712756Z","shell.execute_reply.started":"2024-06-08T00:17:52.695459Z","shell.execute_reply":"2024-06-08T00:17:52.710915Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# 60次元確認用の5件\ntrain_5 = train.head(5).collect()\ntrain_5","metadata":{"execution":{"iopub.status.busy":"2024-06-08T00:17:52.719568Z","iopub.execute_input":"2024-06-08T00:17:52.720031Z","iopub.status.idle":"2024-06-08T00:17:53.103363Z","shell.execute_reply.started":"2024-06-08T00:17:52.719997Z","shell.execute_reply":"2024-06-08T00:17:53.101908Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# 先頭から1,000,000件サンプリング（約10%）\ntrain = train.head(1000000).collect()","metadata":{"execution":{"iopub.status.busy":"2024-06-08T00:17:53.105162Z","iopub.execute_input":"2024-06-08T00:17:53.105632Z","iopub.status.idle":"2024-06-08T00:19:17.968569Z","shell.execute_reply.started":"2024-06-08T00:17:53.105587Z","shell.execute_reply":"2024-06-08T00:19:17.967331Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## 2-1. 説明変数","metadata":{}},{"cell_type":"markdown","source":"## 2-1-1. state\n大気の状態。`state_ps`を除き60次元。次元は垂直レベルを表す。垂直レベルの単位は不明。\n- `state_t` - 大気の温度(air temperature)。単位:$K$\n- `state_q0001` - 比湿（specific humidity）。大気中に含まれる水蒸気量の表現方法の1つで、湿潤空気（水蒸気を含む空気）の質量に対する水蒸気の質量の割合。単位:$kg/kg$\n- `state_q0002` - 雲水混合比(cloud liquid mixing ratio)。単位:$kg/kg$\n- `state_q0003` - 雲氷混合比(cloud ice mixing ratio)。単位:$kg/kg$\n- `state_u` - 帯状風速(zonal wind speed)。単位:$m/s$\n- `state_v` - 子午線風速(meridional wind speed)。単位:$m/s$\n- `state_ps` 表面圧(surface pressure)。単位:$Pa$","metadata":{}},{"cell_type":"code","source":"col_heads = ['state_t_','state_q0001_','state_q0002_','state_q0003_','state_u_','state_v_']\nlabels = ['air temperature [$K$]','specific humidity [$kg/kg$]','cloud liquid mixing ratio [$kg/kg$]','cloud ice mixing ratio [$kg/kg$]','zonal wind speed [$m/s$]','zonal wind speed [$m/s$]']\nplt.figure(figsize=(20,16))\nfor i, col_head in enumerate(col_heads):\n    cols = [c for c in train_5.columns if col_head in c]\n    for train_id, row in enumerate(train_5[cols].iter_rows()):\n        plt.subplot(2,3,i+1)\n        plt.plot(row, range(60), label=f'train_{train_id}')\n        plt.xlabel(labels[i])\n        plt.ylabel('vertical level')\n        plt.ylim([59,0])\n        plt.legend()\nplt.show()","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2024-06-08T00:19:17.970344Z","iopub.execute_input":"2024-06-08T00:19:17.970709Z","iopub.status.idle":"2024-06-08T00:19:19.99234Z","shell.execute_reply.started":"2024-06-08T00:19:17.970678Z","shell.execute_reply":"2024-06-08T00:19:19.991056Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"plt.hist(train['state_ps'],bins=100)\nplt.xlabel('surface pressure [$Pa$]')\nplt.show()","metadata":{"_kg_hide-input":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## 2-1-2. pbuf\n- `pbuf_SOLIN` - 日射量(solar insolation)。単位:$W/m^2$\n- `pbuf_LHFLX` - 地表潜熱フラックス(surface latent heat flux)。単位:$W/m^2$\n- `pbuf_SHFLX` - 地表顕熱フラックス(surface sensible heat flux)。単位:$W/m^2$\n- `pbuf_TAUX` - 帯状面応力(zonal surface stress)。単位:$N/m^2$\n- `pbuf_TAUY` - 子午面応力(meridional surface stress)。単位:$N/m^2$\n- `pbuf_COSZRS` - 太陽天頂角の余弦(cosine of solar zenith angle)。単位:$N/m^2$","metadata":{}},{"cell_type":"code","source":"cols = ['pbuf_SOLIN','pbuf_LHFLX','pbuf_SHFLX','pbuf_TAUX','pbuf_TAUY','pbuf_COSZRS']\nlabels = ['solar insolation [$W/m^2$]','surface latent heat flux [$W/m^2$]','surface sensible heat flux [$W/m^2$]','zonal surface stress [$N/m^2$]','meridional surface stress [$N/m^2$]','cosine of solar zenith angle [$N/m^2$]']\nplt.figure(figsize=(20,8))\nfor i, col in enumerate(cols):\n    plt.subplot(2,3,i+1)\n    plt.hist(train[col],bins=100)\n    plt.xlabel(labels[i])\nplt.show()","metadata":{"_kg_hide-input":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## 2-1-3. cam_in\nアルベドとは、物体表面で反射される光の割合。\n- `cam_in_ALDIF` - 拡散長波放射のアルベド(albedo for diffuse longwave radiation)。\n- `cam_in_ALDIR` - 直接長波放射のアルベド(albedo for direct longwave radiation)。\n- `cam_in_ASDIF` - 拡散短波放射のアルベド(albedo for diffuse shortwave radiation)。\n- `cam_in_ASDIR` - 短波直接放射のアルベド(albedo for direct shortwave radiation)。\n- `cam_in_LWUP` 上向き長波フラックス(upward longwave flux)。単位:$W/m^2$\n- `cam_in_ICEFRAC` - 海氷面積率(sea-ice areal fraction)。\n- `cam_in_LANDFRAC` - 陸域面積率(land areal fraction)。\n- `cam_in_OCNFRAC` - 海洋面積率(ocean areal fraction)。\n- `cam_in_SNOWHLAND` - 陸地の積雪深(snow depth over land)。単位:$m$","metadata":{}},{"cell_type":"code","source":"cols = ['cam_in_ALDIF','cam_in_ALDIR','cam_in_ASDIF','cam_in_ASDIR','cam_in_LWUP','cam_in_ICEFRAC','cam_in_LANDFRAC','cam_in_OCNFRAC','cam_in_SNOWHLAND']\nlabels = ['albedo for diffuse longwave radiation','albedo for direct longwave radiation','albedo for diffuse shortwave radiation','albedo for direct shortwave radiation','upward longwave flux [$W/m^2$]','sea-ice areal fraction','land areal fraction','ocean areal fraction','snow depth over land [$m$]']\nplt.figure(figsize=(20,12))\nfor i, col in enumerate(cols):\n    plt.subplot(3,3,i+1)\n    plt.hist(train[col],bins=100)\n    plt.xlabel(labels[i])\nplt.show()","metadata":{"_kg_hide-input":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## 2-1-4. pbuf\n- `pbuf_ozone` - オゾン体積混合比(ozone volume mixing ratio)。単位:$mol/mol$\n- `pbuf_CH4` - メタン体積混合比(methane volume mixing ratio)。単位:$mol/mol$\n- `pbuf_N2O` - 亜酸化窒素体積混合比(nitrous oxide volume mixing ratio)。単位:$mol/mol$","metadata":{}},{"cell_type":"code","source":"col_heads = ['pbuf_ozone','pbuf_CH4','pbuf_N2O']\nlabels = ['ozone volume mixing ratio [$mol/mol$]','methane volume mixing ratio [$mol/mol$]','nitrous oxide volume mixing ratio [$mol/mol$]']\nplt.figure(figsize=(20,8))\nfor i, col_head in enumerate(col_heads):\n    cols = [c for c in train_5.columns if col_head in c]\n    for train_id, row in enumerate(train_5[cols].iter_rows()):\n        plt.subplot(1,3,i+1)\n        plt.plot(row, range(60), label=f'train_{train_id}')\n        plt.xlabel(labels[i])\n        plt.ylabel('vertical level')\n        plt.ylim([59,0])\n        plt.legend()\nplt.show()","metadata":{"_kg_hide-input":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## 2-2. 目的変数","metadata":{}},{"cell_type":"markdown","source":"## 2-2-1. ptend\n全て60次元。\n- `ptend_t` - 加熱傾向(heating tendency)。単位:$K/s$\n- `ptend_q0001` - 湿潤傾向(moistening tendency)。単位:$kg/kg/s$\n- `ptend_q0002` - 雲水混合比の経時変化(cloud liquid mixing ratio change over time)。単位:$kg/kg/s$\n- `ptend_q0003` - 雲氷混合比の経時変化(cloud ice mixing ratio change over time)。単位:$kg/kg/s$\n- `ptend_u` - 西風加速度(zonal wind acceleration)。単位:$m/s^2$\n- `ptend_v` - 南風加速度(meridional wind acceleration)。単位:$m/s^2$","metadata":{}},{"cell_type":"code","source":"col_heads = ['ptend_t_','ptend_q0001_','ptend_q0002_','ptend_q0003_','ptend_u_','ptend_v_']\nlabels = ['heating tendency [$K/s$]','moistening tendency [$kg/kg/s$]','cloud liquid mixing ratio change over time [$kg/kg/s$]','cloud ice mixing ratio change over time [$kg/kg/s$]','zonal wind acceleration [$m/s/s$]','zonal wind acceleration [$m/s/s$]']\nplt.figure(figsize=(20,16))\nfor i, col_head in enumerate(col_heads):\n    cols = [c for c in train_5.columns if col_head in c]\n    for train_id, row in enumerate(train_5[cols].iter_rows()):\n        plt.subplot(2,3,i+1)\n        plt.plot(row, range(60), label=f'train_{train_id}')\n        plt.xlabel(labels[i])\n        plt.ylabel('vertical level')\n        plt.ylim([59,0])\n        plt.legend()\nplt.show()","metadata":{"_kg_hide-input":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## 2-2-2. cam_out\n- `cam_out_NETSW` - 地表の正味短波長フラックス(net shortwave flux at surface)。単位:$W/m^2$\n- `cam_out_FLWDS` - 地表面下向き長波フラックス(downward longwave flux at surface)。単位:$W/m^2$\n- `cam_out_PRECSC` - 液体水換算の積雪率(snow rate (liquid water equivalent))。単位:$m/s$\n- `cam_out_PRECC` - 降雨率(rain rate)。単位:$m/s$\n- `cam_out_SOLS` - 地表面の下向き可視直達日射量(downward visible direct solar flux to surface)。単位:$W/m^2$\n- `cam_out_SOLL` - 地表面への下向き近赤外直接日射量(downward near-infrared direct solar flux to surface)。単位:$W/m^2$\n- `cam_out_SOLSD` - 地表への下向き拡散日射量(downward diffuse solar flux to surface)。単位:$W/m^2$\n- `cam_out_SOLLD` - 地表への下向き拡散近赤外太陽フラックス(downward diffuse near-infrared solar flux to surface)。単位:$W/m^2$","metadata":{}},{"cell_type":"code","source":"cols = ['cam_out_NETSW','cam_out_FLWDS','cam_out_PRECSC','cam_out_PRECC','cam_out_SOLS','cam_out_SOLL','cam_out_SOLSD','cam_out_SOLLD']\nlabels = ['net shortwave flux at surface [$W/m^2$]','downward longwave flux at surface [$W/m^2$]','snow rate (liquid water equivalent) [$m/s$]','rain rate [$m/s$]','downward visible direct solar flux to surface [$W/m^2$]','downward near-infrared direct solar flux to surface [$W/m^2$]','downward diffuse solar flux to surface [$W/m^2$]','downward diffuse near-infrared solar flux to surface [$W/m^2$]']\nplt.figure(figsize=(20,12))\nfor i, col in enumerate(cols):\n    plt.subplot(3,3,i+1)\n    plt.hist(train[col],bins=100)\n    plt.xlabel(labels[i])\nplt.show()","metadata":{"_kg_hide-input":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## 3. サンプル提出方法\n予測した値に sample_submission.csv の値を列ごとに掛けて提出する必要がある。<br>\n試しに、学習データの目的変数の平均値を出力するコードを書いてみる。<br>\nこの場合、$R^2$スコアの定義から、スコアは0付近になることが予想される。","metadata":{}},{"cell_type":"code","source":"# 本来はsample_idごとに予測する必要があるが、今回は一律で平均値を予測値とする\ntarget_cols = sample_submission.columns[1:]\npred = train.select(pl.col(target_cols)).mean()\npred","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"sample_submission = sample_submission.collect()\nsample_submission","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"for col in pred.columns:\n    sample_submission = sample_submission.with_columns(pl.col(col) * pred[col].max())","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"sample_submission.write_csv(\"submission.csv\")\nsample_submission","metadata":{},"execution_count":null,"outputs":[]}]}