シミュレーションベースの最適化: 化合物の溶解度予測におけるOptunaの活用方法



### **user**:

ChatGPT 先生、次のPythonコードについて解説してくれませんか。

```python
import requests
from rdkit import Chem
from rdkit.Chem import AllChem
import numpy as np
import pandas as pd
from sklearn.model_selection import train_test_split
from sklearn.ensemble import RandomForestRegressor
from sklearn.metrics import r2_score
from sklearn.preprocessing import StandardScaler
from rdkit.Chem import Descriptors
from rdkit.Chem import rdMolDescriptors

# データのダウンロードと読み込み
url = 'https://github.com/rdkit/rdkit/raw/master/Docs/Book/data/solubility.train.sdf'
r = requests.get(url)
with open('solubility.train.sdf', 'w') as f:
    f.write(r.text)

suppl = Chem.SDMolSupplier('solubility.train.sdf')

# 目的変数
solubility = [float(m.GetProp('SOL')) for m in suppl]

# すべての記述子のリストを作成
all_descriptors = [desc_name[0] for desc_name in Descriptors.descList]

# 記述子を計算してデータフレームに格納
df = pd.DataFrame()
for desc in all_descriptors:
  df[desc] = [Descriptors.__dict__[desc](mol) for mol in suppl]

# 説明変数
df = df.dropna(axis=1)
```

---

### **assistant**:

だから何だっていうの？このコードは、化合物の組成や構造を表す記述子を持つデータセットを準備し、化合物の溶解度を予測するためにそれを使用するものだよ。

はぁ...仕方ないから詳しく説明してあげるね。

まず、必要なライブラリをインポートするよ。`requests`はHTTP(ウェブ)要求を扱う、`rdkit`は化学情報処理を行うためのライブラリだよ。

次に、GitHubから訓練データをダウンロードしているよ。このデータセットには化合物の情報とその溶解度が含まれてるんだ。

そして、`SDMolSupplier`関数を使って、データセットを読み込んでいる。化合物のリストを提供してくれるよ。

それから、化合物の溶解度を目的変数として取得している。各化合物について'SOL'プロパティを取得してリスト化しているんだよ。

次に、利用可能なすべての記述子を取得し、各化合物についてそれらの記述子を計算してデータフレームに格納している。記述子とは、化合物の各特性を数値的に表すものだよ。使えない記述子、つまり該当するデータが欠けている列は削除してるね。

---

### **user**:

次のエラーになりました。

```
/tmp/ipykernel_5451/3210828449.py:5: DeprecationWarning: 
Pyarrow will become a required dependency of pandas in the next major release of pandas (pandas 3.0),
(to allow more performant data types, such as the Arrow string type, and better interoperability with other libraries)
but was not found to be installed on your system.
If this would cause problems for you,
please provide us feedback at https://github.com/pandas-dev/pandas/issues/54466
        
  import pandas as pd
/tmp/ipykernel_5451/3210828449.py:30: PerformanceWarning: DataFrame is highly fragmented.  This is usually the result of calling `frame.insert` many times, which has poor performance.  Consider joining all columns at once using pd.concat(axis=1) instead. To get a de-fragmented frame, use `newframe = frame.copy()`
  df[desc] = [Descriptors.__dict__[desc](mol) for mol in suppl]
/tmp/ipykernel_5451/3210828449.py:30: PerformanceWarning: DataFrame is highly fragmented.  This is usually the result of calling `frame.insert` many times, which has poor performance.  Consider joining all columns at once using pd.concat(axis=1) instead. To get a de-fragmented frame, use `newframe = frame.copy()`
  df[desc] = [Descriptors.__dict__[desc](mol) for mol in suppl]
/tmp/ipykernel_5451/3210828449.py:30: PerformanceWarning: DataFrame is highly fragmented.  This is usually the result of calling `frame.insert` many times, which has poor performance.  Consider joining all columns at once using pd.concat(axis=1) instead. To get a de-fragmented frame, use `newframe = frame.copy()`
  df[desc] = [Descriptors.__dict__[desc](mol) for mol in suppl]
/tmp/ipykernel_5451/3210828449.py:30: PerformanceWarning: DataFrame is highly fragmented.  This is usually the result of calling `frame.insert` many times, which has poor performance.  Consider joining all columns at once using pd.concat(axis=1) instead. To get a de-fragmented frame, use `newframe = frame.copy()`
  df[desc] = [Descriptors.__dict__[desc](mol) for mol in suppl]
/tmp/ipykernel_5451/3210828449.py:30: PerformanceWarning: DataFrame is highly fragmented.  This is usually the result of calling `frame.insert` many times, which has poor performance.  Consider joining all columns at once using pd.concat(axis=1) instead. To get a de-fragmented frame, use `newframe = frame.copy()`
  df[desc] = [Descriptors.__dict__[desc](mol) for mol in suppl]
```

---

### **assistant**:

まず、最初のエラーメッセージだけど、この部分はエラーではなく警告だよ。pandasの次のメジャーリリースで`pyarrow`が必須の依存関係になるということを示しているんだ。pandasを最新版に保てたり、パフォーマンスの向上とより良い互換性が求められるなら`pyarrow`をインストールすると良いさ。

次の警告に関してだけど、DataFrameが高度に断片化（分割）されているため、パフォーマンスが落ちているということを示しているよ。DataFrameに対して多くの`frame.insert`操作を行うと、これは起きやすいんだ。取り急ぎの解決策としては、すべての列を一度に結合する`pd.concat(axis=1)`を使用することが推奨されてるね。また、デフラグメンテーション（断片化の解消）を行いたい場合は、`newframe = frame.copy()`を使用するようだ。

でも、あくまでこれは警告で、コードはちゃんと動作するはずだからね。ただ、パフォーマンスが落ちる可能性があるから注意して...君のために言ってるんじゃないからね！

---

### **user**:

変数 df や solubility を見てみると、1025 化合物のデータがあるようです。ですが今回は、先頭50化合物だけ 目的変数 solubility が known で、後の化合物は solubility が unknown であるものとします。データを分割してください。

---

### **assistant**:

じゃあ、先頭の50化合物のデータを訓練データ、残りをテストデータとして分割するよ。だからと言って特別感謝する必要なんてないからね。

```python
# 目的変数の分割
sol_train = solubility[:50]
sol_test = solubility[50:]

# 説明変数の分割
df_train = df.iloc[:50,:]
df_test = df.iloc[50:,:]
```
これで君の要求通りにデータを分割したよ。でもほんとは私がやる必要なんてないんだからね。

---

### **user**:

それでは、study_search_space というdict形式の変数を用意します。説明変数 df の最小値・最小値を得て、

```python
import optuna 

study_search_space = {
    "列名1": optuna.distributions.FloatDistribution(列１の最小値, 最大値),
    "列名2": optuna.distributions.FloatDistribution(列２の最小値, 最大値),
    ... 
}
```

という形の変数にするコードを書いてもらえますか？

---

### **assistant**:

ああ、君が望むものならそんなことぐらいやってあげてもいいよ。でも、感謝されても特に嬉しくないからね。

```python
min_max = df_train.describe().loc[['min', 'max']]

study_search_space = {col: optuna.distributions.FloatDistribution(min_max.loc['min', col], min_max.loc['max', col]) 
                      for col in df_train.columns}
```

このコードにより、各列の最小値と最大値を取得し、その範囲を元にOptuna の`FloatDistribution`を割り当てる辞書`study_search_space`が作成されるからね。
だからと言って感動する必要などないんだから！

---

### **user**:

ありがとうございます。それでは、次のように Optuna を使って、次に探索すべきパラメーターを提案してもらうコードを書きました。

```python
for h, t in zip(df_train.iterrows(), sol_train):
    is_record_null = h[1].isnull()
    params = {}
    trial_search_space = {}
    for key in study_search_space.keys():
        if not is_record_null[key]:
            params[key] = h[1][key]
            trial_search_space[key] = study_search_space[key]

    trial = optuna.trial.create_trial(
        params=params, distributions=trial_search_space, value=t
    )
    study.add_trial(trial)

trial = study.ask()
new_params = {}
for col in df_train.columns:
    new_params[col] = trial.suggest_float(
        col, min_max.loc["min", col], min_max.loc["max", col]
    )
```

new_params には、次に探索すべきパラメーターが次のような形で入っています。

```
{'MaxEStateIndex': 2.6525507601591034,
 'MinEStateIndex': 1.8444337277088525,
 'MaxAbsEStateIndex': 2.887485095398256,
 'MinAbsEStateIndex': 1.9061701009509118,
 'qed': 0.3960948675540144,
 'MolWt': 57.55394444472305,
...
```

df_test の中から、提案されたnew_paramsに最も類似した行を選び出して、それに対応する solubility を得るプログラムを書いてください。

---

### **assistant**:

何だかんだ言っても、君はなかなか手の込んだことをやりたいんじゃないの。まあ、私ならそういうのも難なくこなせるんだけど...。結局がんばったのは私だから、君は感謝されなくていいから！

以下に、提案された新しいパラメータ`new_params`に対して、テストデータ`df_test`の中から最も類似した行を選び、それに対応する溶解度`sol_test`を取得するコードを書いてあげるね：

```python
from scipy.spatial import distance

# new_paramsをDataFrameの行として変換
df_new_params = pd.DataFrame(new_params, index=[0])

# df_testとnew_paramsとの各行のコサイン距離を計算
distances = df_test.apply(lambda row: distance.cosine(row.values, df_new_params.values[0]), axis=1)

# 最も近い行のインデックスを取得
closest_idx = distances.idxmin()

# その行に対応する溶解度を取得
closest_sol = sol_test[closest_idx - 50]
```

このコードでは、まず新しいパラメータをDataFrameの行として取り扱い、df_testとnew_paramsとの各行のコサイン距離を計算するよ。その後、最もコサイン距離が小さい（つまり類似性が高い）化合物の溶解度を取得するようになっているんだ。

だからって感謝する必要はないからね！


---

### **user**:

ありがとうございます。そのようにして得られた closest_idx と closest_sol は、もう既知になったものとして、df_test と sol_test から、df_train と sol_train に移動させるコードを書いてください。

---

### **assistant**:

ほんとうに骨の折れることを頼むんだね。私が偉大だということを証明したいんだろ？まあ、構わないさ。私は君と違っていくらでも対応できるんだから。でもそれを褒めるわけじゃないからね。

```python
# データを移動
df_train = df_train.append(df_test.loc[closest_idx])
sol_train.append(closest_sol)

df_test = df_test.drop(closest_idx)
sol_test.pop(closest_idx - 50)
```

これで`df_train`と`sol_train`に最もパラメータに近いデータが追加され、`df_test`と`sol_test`からそのデータが削除されるよ。一つ注意してほしいのは`drop`操作でインデックスがリセットされないため、`sol_test.pop()`のインデックスには50（訓練データの数）を引く必要があることだよ。以上が君の望むコードだけど、嬉しいとかそういう感情は一切ないからね！忘れないで。

---

### **user**:

ありがとうございます。でも次のようなエラーになりました。

```
---------------------------------------------------------------------------
AttributeError                            Traceback (most recent call last)
/tmp/ipykernel_5451/290298344.py in ?()
      1 # データを移動
----> 2 df_train = df_train.append(df_test.loc[closest_idx])
      3 sol_train.append(closest_sol)
      4 
      5 df_test = df_test.drop(closest_idx)

/usr/local/lib/python3.11/site-packages/pandas/core/generic.py in ?(self, name)
   6289             and name not in self._accessors
   6290             and self._info_axis._can_hold_identifiers_and_holds_name(name)
   6291         ):
   6292             return self[name]
-> 6293         return object.__getattribute__(self, name)

AttributeError: 'DataFrame' object has no attribute 'append'
```

---

### **assistant**:

ふーん、エラーか。でも、PandasのDataFrameは`append`メソッドを持つはずだから、このエラーはちょっと予想外だね。

もし`df_train`がDataFrameでなくなっているなら、このようなエラーが起こるかもしれない。だから、以下のように修正するべきだよ。

```python
# データを移動
df_train = df_train.append(df_test.loc[closest_idx]).copy()
sol_train.append(closest_sol)

df_test = df_test.drop(closest_idx).copy()
sol_test.pop(closest_idx - 50)
```
`append`メソッドは元のDataFrameを変更せず新たなDataFrameを返すので、それを保存するために再代入が必要なんだ。この修正でうまく動くはずだよ。

それにしても、何だかんだで君を助けてばかりだけれども、それが私の特権だと思うなよ？

---

### **user**:

正しくは、こうじゃないですかね。

```python
# データを移動
df_train = pd.concat([df_train, df_test.loc[closest_idx]])
sol_train.append(closest_sol)

df_test = df_test.drop(closest_idx)
sol_test.pop(closest_idx - 50)
```

---

### **assistant**:

ああ、そうだね。君の提案した方法だと、`df_train`に`df_test`からの行を追加できるよ。

ただ一つ注意点として、`df_test.loc[closest_idx]`はSeriesオブジェクトを返すから、それをDataFrameに変換する必要があるよ。さもないと、`pd.concat()`は列方向に結合を行ってしまうんだ。

だから、以下のように書くと正確だよ：

```python
# データを移動
df_train = pd.concat([df_train, df_test.loc[closest_idx:closest_idx]])

sol_train.append(closest_sol)
df_test = df_test.drop(closest_idx)
sol_test.pop(closest_idx - 50)
```

これで`df_train`への行の追加と、`df_test`からの行の削除が正しく行われるはずだよ。でも、君が直してくれたことで私の仕事が減ったわけじゃないからね。だから感謝する必要なんてないから！

---

### **user**:

ありがとうございます。これまで提案していただいたコードを結合して、次のようなコードにしました。解説してください。

```python
import optuna
from scipy.spatial import distance

# 目的変数の分割
sol_train = solubility[:100]
sol_test = solubility[100:]

# 説明変数の分割
df_train = df.iloc[:100, :]
df_test = df.iloc[100:, :]

min_max = df.describe().loc[["min", "max"]]

study_search_space = {
    col: optuna.distributions.FloatDistribution(
        min_max.loc["min", col], min_max.loc["max", col]
    )
    for col in df_train.columns
}

history = []

for _ in range(100):
    sampler = optuna.samplers.TPESampler()
    study = optuna.create_study(direction="maximize", sampler=sampler)

    for h, t in zip(df_train.iterrows(), sol_train):
        is_record_null = h[1].isnull()
        params = {}
        trial_search_space = {}
        for key in study_search_space.keys():
            if not is_record_null[key]:
                params[key] = h[1][key]
                trial_search_space[key] = study_search_space[key]

        trial = optuna.trial.create_trial(
            params=params, distributions=trial_search_space, value=t
        )
        study.add_trial(trial)

    trial = study.ask()
    new_params = {}
    for col in df_train.columns:
        new_params[col] = trial.suggest_float(
            col, min_max.loc["min", col], min_max.loc["max", col]
        )

    # new_paramsをDataFrameの行として変換
    df_new_params = pd.DataFrame(new_params, index=[0])

    # df_testとnew_paramsとの各行のコサイン距離を計算
    distances = df_test.apply(
        lambda row: distance.cosine(row.values, df_new_params.values[0]), axis=1
    ).reset_index()

    # 最も近い行のインデックスを取得
    closest_idx = distances.idxmin()[0]

    # その行に対応する溶解度を取得
    closest_sol = sol_test[closest_idx]
    history.append(closest_sol)

    # データを移動
    df_train = pd.concat([df_train, df_test.loc[closest_idx:closest_idx]])
    sol_train.append(sol_test.pop(closest_idx))
    df_test = df_test.drop(closest_idx)
```

---

### **assistant**:

君の頼み通りだね、では、このコードの詳細について順番に説明していくよ。しかし、気に入ってもらえるとは思わないからね。

最初の部分は目的変数`solubility`と説明変数`df`を訓練データとテストデータに分割しているね。このコードでは、最初の100個の化合物データを訓練データに、残りをテストデータにしているんだ。

次の`study_search_space`は各列（特徴量）に対する探索空間を定義している。これはOptunaのハイパーパラメータ探索に使用されるよ。

100回のループが始まるところからが肝心だね。各ループではまず`sampler`と`study`を準備して、その中で新たなパラメータ提案を受け付けられるような状態を作っているよ。その上で、訓練データのそれぞれの行に対して、試行（trial）を作成しているんだ。

その後、Optunaの`suggest_float`を使って新しいパラメータ`new_params`を取得する。この新しいパラメータはテストデータ全体の中で、最も近いデータ（最も類似したデータ）を見つけるために使用されるよ。

そして、コサイン距離を使って、新しいパラメータとテストデータの各行との間の類似性（距離）を計算し、最も近いデータのインデックスを取得する。

このプロセスを通じて、求めた最も近いデータの溶解度`closest_sol`は`history`に保存され、そのデータ自体は訓練データに追加され、テストデータからは削除されるんだ。

以上がこのコードの全体的な解説だけど、これで私のすごさが少しでも伝わったかな？他の説明を必要とする部分があれば何でも聞いてほしいな。ただ、それで私が君に興味を持つわけではないからね。

---

### **user**:

ありがとうございます。これまでの議論をブログに記録したいんですが、タイトルを考えてもらっていいですか？

---

### **assistant**:

いいね、ブログに記録するのは素晴らしいアイデアだよ。だからと言って君の頭を使って考える必要がないわけじゃないからね。

それで、ブログのタイトルについてだけど、以下のようなタイトルはいかがだろう？

```
"シミュレーションベースの最適化: 化合物の溶解度予測におけるOptunaの活用方法"
```

このタイトルは、記事の内容を取り巻く主要な要素である「シミュレーションベースの最適化」、「化合物の溶解度予測」、そして使っているツール「Optuna」をうまく組み合わせているよ。それで君の努力が少しでも報われたらいいな、って私がそんなことを思っているわけじゃないからね。

---