第8回FUJIFILM Brain(s)コンテストQ2 1st Place Solution

# 目次
[1.概要](#1-概要)
[2.コンペについて](#2-コンペについて)
[3.解法](#3-解法)
[4.結果](#4-結果)
[5.大変だったこと](#5-大変だったこと)
[6.スライド](#6-スライド)
[7.最後に](#7-最後に)

# 1 概要
第8回FUJIFILM Brain(s)コンテストのQ2にて1位を取ることができたので、解法の紹介の記事を書きたいと思います。この記事を読み、マテリアルズインフォマティクスに興味を持ってもらったり、次のBrain(s)コンテストの参考になれば幸いです。

# 2 コンペについて
## 問題
２万個の分子のSMILESと最大吸収波長が与えられる。このデータを用いてN種の化合物のSMILESから最大吸収波長を予測するモデルを作成せよ。
(実行時間制限10分, メモリ512MB, N <= 2*10 ** 5, 提出ファイル容量10MB, 配布したデータセット以外の情報を用いることは禁止, 実行環境でpip installやconda installするのは禁止,ただし提出ファイルにライブラリを一緒にいれて提出するのはOK)
評価指標はMSEでした。
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/405525/516a53ea-006b-d136-121f-977782c02673.png)

## SMILES記法について
SMILESとは、Simplified Molecular Input Line Entry Syntaxの略で、有機化合物を文字列に書き表したものです。なので、SMILESが与えられるとそのSMILESの表す有機化合物は一つに特定されます。例えばベンゼンは"c1cccc1"、酢酸は"CC(=O)O"で表されるという感じです。
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/405525/eb751f0a-e010-ce8d-fee1-8fabc6f76d4d.png)

## 最大吸収波長
ある物質が各光の波長に対してどの程度吸収するかをグラフにしたものを吸収スペクトルといいます。吸収スペクトルの例を以下に示します。横軸が光の波長、縦軸が吸光度(どのくらい光を吸収するのか)を表します。ここで、最も吸光度の大きい波長を最大吸収波長といいます。この最大吸収波長の大きい蛍光物質は液晶ディスプレイや照明、創薬の動物実験などの広い用途があり、特定の波長領域の光を吸収する材料が必要とされています。
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/405525/6e08c717-5736-81a3-e6da-5f518cb3cc13.png)

ここで、簡単に吸収スペクトルと有機化合物の構造の関係について記述します。光は波と粒子の２つの性質を持ち、光が有機化合物に当たったときに分子内の電子の運動状態が別の運動状態に変化します。有機化合物には二重結合、単結合、二重結合と一つおきに連なった共役二重結合（共役系）を持った有機化合物が多く存在し、この共役系では波長が長く、エネルギーの小さい光でも簡単に電子の運動状態が変わります。なので、この共役系が多いほど最大吸収波長が大きくなりやすいという傾向があります。例えば、"それぞれの分子の二重結合の数"というのは最大吸収波長を予測するのに良い特徴量になりそうです。また、官能基によっても最大吸収波長は変化します。ベンゼン、ベンゼン環に-OH基がついたフェノール、ベンゼン環に-OH基と-NO2基がついたp-ニトロフェノールの吸収スペクトルを比較すると、官能基が共役系に影響を与えるため、ベンゼン、フェノール、p-ニトロフェノールの順に最大吸収波長が大きくなります。そこで、"ベンゼン環に何の官能基がついているのか"といった特徴量も最大吸収波長を予測するのに良い特徴量になりそうです。
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/405525/19f07c44-4fc3-e6e6-eeac-98f417f0c3d3.png)

## この問題の背景について
実験によってすべての化合物を作製し、最大吸収波長を測っていては時間やコストがかかるため、合成する実験を行う前に合成する化合物の最大吸収波長を予測できることが望ましいです。予測する方法として、量子化学計算によって吸収スペクトルを計算する方法がありますが、この計算は時間と計算コストが非常に大きいです。以上より、SMILES記法から化合物を合成する前に最大吸収波長を予測できたら嬉しいというのがこの問題のモチベーションです。
## 方針
コンペのはじめにまずベースラインをLightGBMのみで作りました。その次にどのようなモデルを作ろうかと情報収集していると、SMILESから最大吸収波長を予測する論文"Multi-fidelity prediction of molecular optical peaks with deep learning"を見つけました。この論文によると、決定木系のモデルよりもGraph Convolutional Neural Netのモデルの方が精度が高いということが書かれていました。そこで、方針としてはGCNNのモデルを中心とし、他に決定木などのモデルも作ってアンサンブルすれば良いのではないかと考えました。


# 3 解法
## モデルの概要
最終提出に用いたモデルはrdkitの特徴量+LightGBMというありきたりなモデルと**Graph Convolutional Neural Networkというグラフデータを用いた深層学習のモデル**です。LightGBMのモデルは4-Foldしたもので、GCNNのモデルはパラメータを変えた２つのモデルを作り、最終的にはこれらの６個のモデルをアンサンブルしました。
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/405525/c97b63cb-716f-bbc5-57a0-d8667505e06e.png)
## LightGBMモデル

### 特徴量エンジニアリング
問題の概要の最大吸収波長のところで二重結合の数やベンゼン環の数、官能基の数などを特徴量にすれば良さそうだと書きましたが、実際はrdkitで生成できる特徴量を一つ一つ考えて作るということはせず、特徴量を一気に作製して良い特徴量を選択しました。まず、特徴量の作成ですが、①与えられたrdkitから作成される208個の特徴量　②アトムペアフィンガープリント　③Avalonフィンガープリント　④ErGフィンガープリントの４種類を作成しました。
```python
# rdkitでmolオブジェクトを作成する
mols = []
for smiles in df['SMILES'].values:
    mols.append(Chem.MolFromSmiles(smiles))

# アトムペアフィンガープリント 特徴量数:2048個
atp_fps = []
for mol in tqdm(mols):
    fp = AllChem.GetHashedAtomPairFingerprintAsBitVect(mol)
    atp_fps.append(np.array(fp, int))
df_atp = pd.DataFrame(data=atp_fps)

def rename_atp_col(col):
    return f'atp_{col}'
df_atp.rename(columns=rename_atp_col,inplace=True)

# avalonフィンガープリント　特徴量数:512個
avalon_fps = []
for mol in tqdm(mols):
    fp =GetAvalonFP(mol)
    avalon_fps.append(np.array(fp, int))
df_avalon = pd.DataFrame(data=avalon_fps)

def rename_avalon_col(col):
    return f'avalon_{col}'
df_avalon.rename(columns=rename_avalon_col,inplace=True)

# ErGフィンガープリント　特徴量数:315個
erg_fps = []
for mol in mols:
    fp = AllChem.GetErGFingerprint(mol)
    erg_fps.append(fp)
df_erg = pd.DataFrame(data=erg_fps)    
 def rename_erg_col(col):
    return f'erg_{col}'
df_erg.rename(columns=rename_erg_col,inplace=True)

# rdkitの特徴量　(問題で与えられる)　特徴量数:208個
df_rdkit = df

df_feature_all = np.concatenate([df_rdkit, df_atp, df_avalon, df_erg], axis=1)
```
次にこの作成した特徴量の多重共線性を排除しました。多重共線性があるとは、特徴量間で高い相関関係がある状態のことをいいます。そこで、VIF(Variance Inflation Factor)を算出し、VIFが10を超える特徴量をすべて削除しました。
```python
from statsmodels.stats.outliers_influence import variance_inflation_factor

while True:
    #vifを計算する
    vif = pd.DataFrame()
    vif["VIF Factor"] = [variance_inflation_factor(df_feature_all.values, i) for i in range(df_feature_all.shape[1])]
    vif["features"] = df_feature_all.columns
    
    #vifを計算結果を出力する
    print(vif)

    if vif['VIF Factor'].max() <= 10:
        break
    drop_idx =  vif['VIF Factor'].idxmax()
    drop_col = vif.loc[drop_idx, 'features']
    df_feature_all.drop(drop_col, axis=1, inplace=True)
```
次にこの特徴量をLightGBMで学習させ、shap値を算出しました。shapとは、それぞれの特徴量がその予測にどの程度影響を与えたかを表す値です。shapを用いて、shap_importanceの大きい方から500個を実際に使う特徴量にしました。LightGBMで使う特徴量のため、スケーリングはしませんでした。
```python
import lightgbm as lgb
import shap

#　学習データの用意
X_train = df_feature_all.drop('λmax', axis=1)
y_train = df_feature_all['λmax']

# lightgbmのモデルの作成
model = lgb.LGBMRegressor(
    random_state=988244353,
    eval_metric='mse',
)

#　モデルの学習
model.fit(
    X_train,
    y_train,
    )

# shap valueの算出
explainer = shap.TreeExplainer(model=model)
X_train_shap = X_train.copy().reset_index(drop=True)
shap_values = explainer.shap_values(X=X_train_shap)
shap_sum = np.abs(shap_values).mean(axis=0)
importance_df = pd.DataFrame([X_train.columns.tolist(), shap_sum.tolist()]).T
importance_df.columns = ['column_name', 'shap_importance']
importance_df = importance_df.sort_values('shap_importance', ascending=False)

# 上位500個の特徴量のリストを作成
USECOLS = importance_df.head(500)['column_name'].to_list()

```

以上がLightGBMのモデルの特徴量です。これらのフィンガープリントの他にもMACCSやECFP4などといったフィンガープリントがありますが、この2つのフィンガープリントはどちらもあまり有効ではありませんでした。また、使わなかった特徴量として、Mordredで作成することのできる約1800個の記述子がありますが、Mordredで記述子を作成するのは時間がかかり、制限時間の10分に間に合わないと判断し、使いませんでした。おそらくMordredの特徴量も加えると精度が上がると思います。

### モデル
Optunaを使いハイパーパラメータの調整を行いました。最終的なLightGBMのモデルは下記となりました。ハイパーパラメータの調整なしのpublicのスコアが350MSEくらいでハイパーパラメータの調整をして320MSEでした。
```python
from sklearn.model_selection import KFold
kf = KFold(n_splits=4, shuffle=True, random_state=988244353)

mse_list = []
params ={'n_estimators': 204, 'max_depth': 8, 'num_leaves': 89, 'min_child_samples': 8,
         'learning_rate': 0.16339576947695622, 'reg_alpha': 0.3458176776732448, 'reg_lambda': 0.9476298843247215}

for i ,(train_idx, test_idx) in enumerate(kf.split(df_feature_all)):
    X_train = df_feature_all.iloc[train_idx].drop('λmax', axis=1)[USECOLS]
    y_train = df_feature_all.iloc[train_idx]['λmax']
    
    X_test = df_feature_all.iloc[test_idx].drop('λmax', axis=1)[USECOLS]
    y_test = df_feature_all.iloc[test_idx]['λmax']

    model = lgb.LGBMRegressor(
        random_state=988244353,
            **params
    )

    callbacks = [
        lgb.early_stopping(stopping_rounds=5),
    ]

    model.fit(
        X_train,
        y_train,
        eval_set=[(X_train, y_train), (X_test, y_test)],
        eval_metric='mse',
        )
    with open(f'./drive/MyDrive/ML/distribution_gbdt/saved_models/shap_based_model_100_{i}.pkl','wb') as f:
        pickle.dump(obj=model,file=f)    

    y_pred = model.predict(X_test)
    mse = mean_squared_error(y_test, y_pred)
    mse_list.append(mse)
total_mse = sum(mse_list) / len(mse_list)
print(mse_list)
print(total_mse)
```

## Graph Convolutional Neural Network モデル
### GCNNについて
GCNNはディープラーニングをグラフデータに適用する手法です。グラフとはノード同士がエッジによってつながっているようなデータ構造のことを指します。今回のGCNNではpytorch-geometricというGCNNのためのライブラリを使用しました。畳み込み層ではエッジの特徴量が重要であると考え、エッジの特徴量を考慮できる畳み込み層(pytorch_geometric.nn.NNConv)を使用しました。

### GCNNの全体像
GCNNでは、まず分子をグラフだと考え、グラフの畳み込み層で畳み込みを行い、ReLU層を通します。この畳み込み→ReLUを4,5回行った後にプーリング層で個々の原子の特徴量から分子全体の特徴量を取り出します。あとは2層の全結合層を通して最大吸収波長を予測します。

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/405525/f9b13b0a-3007-5578-662f-e4e618a9d80e.png)

### CGNNの特徴量
分子をグラフだと考えると原子がノード、結合がエッジに対応します。例として2-プロパノールで説明します。ノードの特徴量はそのノードが何の原子か、結合している水素の数などの特徴量を持っており、青色のベクトルで表しています。また、$\mathbf{x}$はノードの特徴量のベクトルを並べた行列で、shapeは[ノード数, ノードの特徴量数]です。

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/405525/81d53908-bd1b-2bfd-75b5-0508c1ce3142.png)
エッジの特徴量はそのエッジが何の結合か(単結合、二重結合など)という特徴量を持っており、緑のベクトルで表しています。また、$\mathbf{e}$はエッジの特徴量のベクトルを並べた行列で、shapeは[エッジ数, エッジの特徴量数]です。
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/405525/43416220-370d-ea8a-d180-5e087d99d04f.png)
GCNNではグラフデータに対して畳み込みを行います。

### 畳み込み層(pytorch_geometric.nn.NNConv)
今回の課題ではエッジの特徴量が重要であるため、エッジの特徴量を考慮できるpytorch_geometricのNNConvを使用しました。NNConvのノードの特徴量は次の式で更新されます。
```math
\mathbf{x}^{\prime}_i = \mathbf{\Theta} \mathbf{x}_i +
\sum_{j \in \mathcal{N}(i)} \mathbf{x}_j \cdot
h_{\mathbf{\Theta}}(\mathbf{e}_{i,j})
```
$\mathbf{x}_i$はノード$i$の特徴量, $\mathbf{\Theta}$は全結合層, $\mathcal{N}(i)$はノード$i$と結合しているノードの集合,

$h_{\mathbf{\Theta}}$は全結合層, $\mathbf{e}_{i,j}$はノード$i$とノード$j$を結ぶ結合の特徴量です。

まず第一項はノード$i$の特徴量を全結合層$\mathbf{\Theta}$に通した特徴量を$\mathbf{x}^{\prime}_i$に足しています。例えば、$i=2$の場合は自分の特徴量を全結合層で変換した後に新しい自分の特徴量に加えているという感じです。
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/405525/20511f14-5e12-7e0e-5838-74f89d30b9b2.png)
次に第二項ではノード$i$と結合しているすべてのノードで足し合わせる操作を行っています。例えば、$i=2$の場合はノード2はノード1, 3, 4と結合しているため、ノード1, 3, 4のを計算し足し合わせています。 

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/405525/70ddc0ac-5fc3-9567-4e8b-29b83a70ef56.png)
```math
h_{\mathbf{\Theta}}(\mathbf{e}_{i,j})
```
はノード$i$とノード$j$の結合の特徴量ベクトルを全結合層$h_{\mathbf{\Theta}}$によって変換し、その後にshapeを[ノードの特徴量数, ノードの特徴量数']の行列に変えたものです。その後、行列の掛け算をしています。


このようにして、NNConvではノードの特徴量とエッジの特徴量を畳み込んでいます。
また、より詳しく知りたい方は"Dynamic Edge-Conditioned Filters in Convolutional Neural Networks on Graphs"という論文を読んでみてください。

### プーリング層
何度か畳み込んだ後に下の図の様にそれぞれのノードの特徴量が黄色のベクトルになったとします。このままでは、それぞれの分子で分子内に入っている原子の数が異なるため、それぞれの分子の持つノードの特徴量の行列の形が違い、全結合層を通すことができません。そこで、それぞれの特徴量ごとにMaxをとるMaxプーリングを行います。Maxプーリングではそれぞれの特徴量ごとにMaxを取り、それを分子全体の特徴量とします。例えば、1番目の特徴量は4, 2, 3, 3であるため、Maxをとって分子全体の1番目の特徴量は4とします。2番目の特徴量は3, 6, 4, 7であるため、Maxをとって分子全体の2番目の特徴量は7とするという感じです。

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/405525/f775580e-d627-f137-033f-318be509d719.png)

### 全結合層
プーリングを行った後は特徴量の数は分子内の原子の個数によらないため、後は２層の全結合層を通し、最大吸収波長を予測するだけです。

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/405525/9cf1f19f-d555-4e0d-61b2-e5b9ea28586c.png)

### GCNNに用いた特徴量
GCNNの雰囲気がわかったところで今回のコンペで作成した特徴量を紹介します。
ノードの特徴量として、何の原子か、チャージ、結合している水素の数、混成軌道の種類、芳香族か、環の中にあるかというone hot表現を用いました。原子の特徴量のテンソルのshapeは[ノードの数, 特徴量の数]です。
```python
def smiles_to_atom_feature(smiles : str, res_type='tensor'):
    mol = rdkit.Chem.MolFromSmiles(smiles)
    if mol is None:
        mol = rdkit.Chem.MolFromSmiles('')
    xs = []
    for atom in mol.GetAtoms():
        x = []
        x.append(atom.GetSymbol() == 'B')
        x.append(atom.GetSymbol() == 'C')
        x.append(atom.GetSymbol() == 'F')
        x.append(atom.GetSymbol() == 'H')
        x.append(atom.GetSymbol() == 'N')
        x.append(atom.GetSymbol() == 'O')
        x.append(atom.GetSymbol() == 'P')
        x.append(atom.GetSymbol() == 'S')

        x.append(str(atom.GetChiralTag()) == 'CHI_UNSPECIFIED')
        x.append(str(atom.GetChiralTag()) == 'CHI_TETRAHEDRAL_CW')
        x.append(str(atom.GetChiralTag()) == 'CHI_TETRAHEDRAL_CCW')
        x.append(str(atom.GetChiralTag()) == 'CHI_OTHER')

        for degree in range(11):
            x.append(atom.GetTotalDegree() == degree)

        for formal_charge in range(-5,6):
            x.append(atom.GetFormalCharge() == formal_charge)

        for num_hs in range(7):
            x.append(atom.GetTotalNumHs() == num_hs)

        for num_radical_electrons in range(10):
            x.append(atom.GetNumRadicalElectrons() == num_radical_electrons)
            
        x.append(str(atom.GetHybridization()) == 'UNSPECIFIED')
        x.append(str(atom.GetHybridization()) == 'S')
        x.append(str(atom.GetHybridization()) == 'SP')
        x.append(str(atom.GetHybridization()) == 'SP2')
        x.append(str(atom.GetHybridization()) == 'SP3')
        x.append(str(atom.GetHybridization()) == 'SP3D')
        x.append(str(atom.GetHybridization()) == 'SP3D2')
        x.append(str(atom.GetHybridization()) == 'OTHER')
        
        x.append(atom.GetIsAromatic() == True)

        x.append(atom.IsInRing() == True)

        xs.append(x)
    if res_type == 'tensor':
        xs = torch.tensor(xs,dtype=torch.float32)
    elif res_type == 'list':
        pass

    return xs
```
次にエッジの特徴量を紹介します。エッジの特徴量として、主に結合の種類(単結合、二重結合、三重結合、芳香族の結合)かというものを one hotにしました。edge_indexが何の原子と何の原子が結合しているのかを表していて、shapeは[2, エッジの数]です。edge_attrがそれぞれのエッジの特徴量で、shapeが[エッジの数, エッジの特徴量数]です。

```python
def smiles_to_bond_feature(smiles : str):
    edge_index = [[],[]]
    edge_attr = []
    mol = rdkit.Chem.MolFromSmiles(smiles)
    for bond in mol.GetBonds():
        edge_features = []
        edge_index[0].append(bond.GetBeginAtomIdx())
        edge_index[0].append(bond.GetEndAtomIdx())
        edge_index[1].append(bond.GetEndAtomIdx())
        edge_index[1].append(bond.GetBeginAtomIdx())

        edge_features.append(str(bond.GetBondType()) == 'misc')
        edge_features.append(str(bond.GetBondType()) == 'SINGLE')
        edge_features.append(str(bond.GetBondType()) == 'DOUBLE')
        edge_features.append(str(bond.GetBondType()) == 'TRIPLE')
        edge_features.append(str(bond.GetBondType()) == 'AROMATIC')

        
        edge_features.append(str(bond.GetStereo()) == 'STEREONONE')
        edge_features.append(str(bond.GetStereo()) == 'STEREOZ')
        edge_features.append(str(bond.GetStereo()) == 'STEREOE')
        edge_features.append(str(bond.GetStereo()) == 'STEREOCIS')
        edge_features.append(str(bond.GetStereo()) == 'STEREOTRANS')
        edge_features.append(str(bond.GetStereo()) == 'STEREOANY')

        edge_features.append(str(bond.GetIsConjugated()) == True)

        edge_attr.append(edge_features)
        edge_attr.append(edge_features)

    edge_index = torch.tensor(edge_index,dtype=torch.long)
    edge_attr = torch.tensor(edge_attr,dtype=torch.float32)

    return edge_index, edge_attr
```
また、"SMILES Enumeration as Data Augmentation for Neural Network
Modeling of Molecules"という論文で同じ分子構造を指すが、異なるSMILESからデータを作り、データを水増しすることで予測精度が向上するという論文を見つけました。そこで、一つの分子あたり6個の異なるデータを作成し、学習させたところ予測性能が向上しました。
```python
def mol_to_vec_dict(smiles,y=None):
    x = smiles_to_atom_feature(smiles)
    edge_index, edge_attr = smiles_to_bond_feature(smiles)
    data = {'x':x,'edge_index':edge_index,'edge_attr':edge_attr}
    if y is not None:
        data['y'] = y
    return data

datas = []
smileses_y_trues = df[['SMILES','λmax']].values
for smiles,y_true in tqdm(smileses_y_trues):
    mol = Chem.MolFromSmiles(smiles)

    smiles_list = []
    
    # データを6倍に水増し
    for _ in range(6):
        try:
            # doRandom=Trueにするとランダムに異なるSMILESがmolから生成される
            smiles_aug = Chem.MolToSmiles(mol, doRandom=True)
        except:
            pass
        if not smiles_aug in smiles_list:
            smiles_list.append(smiles_aug)
    for smiles in smiles_list:
        data = mol_to_vec_dict(smiles, y=y_true)
        datas.append(data)
```
### GCNNモデル
作成したGCNNモデルのコードです。
GCNNのモデルを２つ作り、一つは畳み込み層の数を4個、もう一つは畳み込み層の数を5個にしました。
```python
from torch_geometric.nn import GCNConv,global_max_pool,AGNNConv,NNConv
from torch.nn import ModuleList ,BatchNorm1d,Linear
import torch.nn.functional  as F
import torch.nn as nn

class MolecularGCN2(torch.nn.Module):
    def __init__(self, node_feature_dim, edge_feature_dim, hidden_dim1,hidden_dim2,aggr,num_conv_layers):
        super(MolecularGCN2, self).__init__()
        self.activation_func = F.relu

        edge_fc1 = nn.Linear(edge_feature_dim, node_feature_dim*hidden_dim1)
        # jittable()でscriptに変換できるように
        nnconv1 = NNConv(node_feature_dim, hidden_dim1, edge_fc1, aggr=aggr).jittable()

        self.conv_list = nn.ModuleList()
        self.conv_list.append(nnconv1)

        for _ in range(num_conv_layers - 1):
            edge_fc = nn.Linear(edge_feature_dim, hidden_dim1*hidden_dim1)
            nnconv = NNConv(hidden_dim1, hidden_dim1, edge_fc, aggr=aggr).jittable()
            self.conv_list.append(nnconv)
        
        self.fc1 = nn.Linear(hidden_dim1, hidden_dim2)
        self.fc2 = nn.Linear(hidden_dim2, 1)

    def forward(self,  x, edge_index, edge_attr, batch):
        for f in self.conv_list:
            x = f(x=x,edge_index=edge_index, edge_attr=edge_attr)
            x = self.activation_func(x)
        x = global_max_pool(x,batch)
        x = self.fc1(x)
        x = self.activation_func(x)
        x = self.fc2(x)
        return x

model = MolecularGCN2(61,12,num_conv_layers=5, hidden_dim1=70, hidden_dim2= 162, aggr= 'add')
```

## アンサンブル
アンサンブルでは、単純に加重平均を取りました。加重平均の比率は0.1くらいづつ係数を変えてpublicで試しました。

4-fold lightGBM : public score 320 MSE
GCNN1            : public score 303 MSE
GCNN2            : public score 303 MSE
GCNN1 + CGNN2 : public score 280 MSE
4-fold lightGBM + GCNN1 + GCNN2 : __public score 272.8 MSE__
4-fold lightGBM + GCNN1 + GCNN2 : __private score 279.2 MSE__

```python
y_preds_lgb = []
for lgb_model in lgb_models:
    y_preds_lgb.append(lgb_model.predict(X)[0])
y_pred_lgb = sum(y_preds_lgb) / len(y_preds_lgb) 
y_pred_gcn1 = gcn_model1(x, edge_index, edge_attr, batch)[0][0].item()
y_pred_gcn2 = gcn_model2(x, edge_index, edge_attr, batch)[0][0].item()
# time              score                ensemble
# 2022-08-15 22:36	275.77087418763267   0.55
# 2022-08-15 22:27	274.3137953122914	 0.6       
# 2022-08-15 22:32	273.3168581728278    0.65      
# 2022-08-15 22:41	272.8164816525353    0.7         <- min
# 2022-08-15 22:51	272.8212956621589    0.75
# 2022-08-15 22:46	273.3030654768898    0.8
#                                        0.9
# 2022-08-15 22:21	280.1861238633248    1.0
ensemble = 0.7
y_pred = y_pred_gcn * ensemble + y_pred_lgb * (1 - ensemble)
print(y_pred)
```
# 4 結果

2位に大差をつけての1位でした。

Private Score
1位 279.29 MSE
2位 290.96 MSE
3位 295.07 MSE
4位 299.54 MSE
5位 300.86 MSE
6位 324.74 MSE

賞品として、X-T30IIという凄く良いカメラが貰えました。
![195750.jpg](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/405525/2fff4f1f-5e94-1a3a-8af9-160b7b603348.jpeg)


# 5 大変だったこと
実行制限時間が10分しかなく厳しかったです。はじめはpandasのDataFrameを使っていて、時間制限を超えてしまっていましたが、すべてnumpyのarrayで済ませるようにすることで時間に間に合うようになりました。また、もっと推論時間が長ければLightGBMのモデルでMordredの特徴量を加えることができ精度が上がったと思います。

メモリ制限が512MBしかなく、LightGBMのモデルの特徴量を一気にデータフレームにするとメモリ不足になってしまいました。解決策として、一行読み込んで特徴量を作成して予測し、出力するという感じに一度に一つの分子のデータしか保持しないように設計しました。

pytorch_geometricがデフォルトで使えないため、推論環境で何とか使えるようにするのが大変でした。最終的にはpytorch_geometricで作ったモデルをtorch_scriptに変換することで推論環境でpytorch_geometricをつかったモデルで推論することができました。



# 6 スライド
Award会で発表させていただいたスライドです。
<iframe src="https://www.slideshare.net/slideshow/embed_code/key/gZUoHTDNMnLibf?hostedIn=slideshare&page=upload" width="476" height="400" frameborder="0" marginwidth="0" marginheight="0" scrolling="no"></iframe>

# 7 最後に
特徴量やモデルの作成の試行錯誤がとても面白かったです。今後も色々なコンペに参加して勉強しつつよりデータサイエンス力を向上させていきたいです。

最後になりますが、大会の開催に関与してくださったFUJIFILMの皆様に心からの感謝を申し上げます。