ベイズ最適化をステップ実行してみた。(実施例付き)

# 概要

ベイズ最適化とは、未知の関数`f(x)`に対する、最大値もしくは最小値を求めるための手法の一つです。
その性質上、次のような問題ケースに向いているとされます。

- `f(x)`の形状が未知である。凸包かどうか、双峰性を持つかなども分からないものとする
- ある入力`x`に対する関数値`f(x)`を1回計算するだけでもコストが嵩む(時間が掛かる、など)
- 大域的な最適解を求めたいが、逐次的に「推定される`f(x)`の最大値(最小値)」を求めたい

具体的なアルゴリズムや数学的解説については、次のスライドが分かりやすかったです。

　[機械学習のためのベイズ最適化入門 - Slideshare](https://www.slideshare.net/hoxo_m/ss-77421091)
　[01.ベイズ最適化概説](https://colab.research.google.com/github/TatsuyaKatayama/OpenFOAMandBaysianOpt_Notebooks/blob/master/01.%E3%83%99%E3%82%A4%E3%82%BA%E6%9C%80%E9%81%A9%E5%8C%96%E6%A6%82%E8%AA%AC.ipynb)

……で、ちょっと最適化したい事象を見つけたので、それに向けたプログラムを構築してみました。

# ベイズ最適化に使用したライブラリとサンプルコード

前述したPython Notebook資料でも使われていた、[GPyOpt](https://github.com/SheffieldML/GPyOpt)ライブラリを使用します。
これはPythonで簡単にベイズ最適化計算が行えるライブラリで、サンプルコードを書くとこんな感じ。

```python
import GPyOpt
from numpy.core.multiarray import ndarray


def func(x: float, y: float) -> float:
    # Rosenbrock function
    # 今回の最適化対象。厳密解は(x, y) = (1, 1) の時にf(x, y) = 0
    val = 100.0 * (y - x * x) ** 2 + (1 - x) ** 2
    print(f'f({x}, {y}) = {val}')
    return val


def func_wrapper(arr: ndarray) -> float:
    # デフォルトでは、全引数がndarray型で纏まって渡されて分かりにくいので、
    # 実際の関数との間に立ってラップするための関数
    return func(arr[:, 0][0], arr[:, 1][0])


if __name__ == '__main__':
    # 各変数における、名称(変数名)・連続値か離散値か・値の範囲。
    # type=continuousではなくtype=discreteとすると、離散値も扱える
    bounds = [
        {'name': 'x', 'type': 'continuous', 'domain': (-2, 2)},
        {'name': 'y', 'type': 'continuous', 'domain': (-1, 3)},
        # {'name': 'y', 'type': 'discrete', 'domain': tuple(range(-1, 3))},
        ]

    # 問題定義。maximize=Trueなら最大化を目指す
    problem = GPyOpt.methods.BayesianOptimization(func_wrapper, domain=bounds, maximize=False, verbosity=True)

    # 実行。verbosity=Trueにすると、経過した計算時間も表示される
    max_iter = 25
    problem.run_optimization(max_iter, verbosity=True)

    # 最適化処理での、各施行における評価値と、最良値の推移をグラフ表示する
    problem.plot_convergence()

    # 最良解(今まで見た中で最も成績が良かった評価値)を表示
    print(problem.x_opt, problem.fx_opt)

    # 推定された関数形状のグラフ(平均・分散)と、次にどこを探索する予定かを示すグラフ(獲得関数)を表示する
    problem.plot_acquisition()
```

**各施行における評価値と、最良値の推移**
![Figure_1.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/29802/cd14057e-1ff0-b447-22e7-d163ad364d17.png)

**推定された関数形状のグラフ(平均・分散)と、次にどこを探索する予定かを示すグラフ(獲得関数)**
![Figure_2.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/29802/12ca3e5e-84a3-c346-923f-9b5d18e500eb.png)

# 最適化したい事象って？

**レンズの色味を合わせるためのカメラのホワイトバランス設定です。**

……と言っても大多数の人には伝わらないと思いますので、段階を踏んで説明していきます。

## ホワイトバランスとは

おおよそ感覚で察していただけると思いますが、たとえカメラで撮る物体が同じだったとしても、
その物体が「どのような光(照明)に晒されているか」によって、その色味は変化します。
例えば**晴れている時**。例えば**曇り空な時**。例えば**電球照明の室内にいる時**……などなど。

その、照明の色味のことを「色温度」と呼びます。
また、色温度から逆算して写真を補正し、正しい色味に戻すことを「ホワイトバランス調整」と呼びます。

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/29802/2bae11ab-4e96-4e18-3158-18b30a10bcaf.png)
(https://av.jpn.support.panasonic.com/support/dsc/knowhow/knowhow31.html より引用)

## レンズの色味とは

レンズというのは基本的にガラスで出来ています。
また、板ガラスの側面を見れば分かるように、ガラスには僅かながら「色」が付いています。
すると、レンズを含むガラスを通った光は、ほんの僅かに「色味」が付いてしまうことになります。

さらに、昨今のカメラレンズには、薄いコーティングが施されています。
これにより、写真が更に綺麗に映ることになるのですが、このコーティングにも僅かながら「色」が付いています。

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/29802/bef6d5e2-b003-21fb-8acd-ae70fe644a8d.png)
(https://av.jpn.support.panasonic.com/support/dsc/knowhow/knowhow17.html より引用)

加えて、全てのメーカーが同一のガラス・コーティングを使用しているわけではありませんし、
使用しているレンズの種類・枚数はカメラレンズ毎に異なります。
よって、**カメラレンズはそれぞれ僅かながら色味の差が存在する**のです。

……ただ、メーカー固有の技術となるコーティングの方が影響が強く、「レンズメーカーの差＞ガラスの差」といった状態です。
(まあ同一メーカーでもブランドによって発色が違ったりしますけどね)

## 今回の案件

普段は[Panasonic LUMIX G9 PRO](https://panasonic.jp/dc/products/g_series/g9pro.html)にPanasonicの[LEICA DG LENS](https://panasonic.jp/dc/lens/leica_dg_lens.html)を付けているのですが、[M.ZUIKO PROレンズ](https://cameras.olympus.com/zuiko/ja-jp/mzuikopro/)を使うこともあります。
ただ両者は**微妙に色味が違う**ので、どれぐらい補正すればいいのかが気にかかっていました。

そこで、前者で白っぽい壁を撮影し(※ピントを意図的にボカして一様な絵にする)、後者でも同じ被写体を**ホワイトバランスを都度変更しながら撮影**することで、どのホワイトバランス調整値が補正として適切かを判断します。

上記の画像では「色温度」という1次元しか調整軸がありませんが、G9 PROのホワイトバランス調整は「A(オレンジ)～B(青)軸」と「G(緑)～M(赤)軸」の2軸があり、これらを±9の範囲内で調整することで最適解を探索します。

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/29802/367cfb8a-e1d9-64db-5793-0067d36ea760.png)

(https://panasonic.jp/p-db/contents/manualdl/1428353073217.pdf より引用)

# ベイズ最適化をステップ実行だ！

まず最適化したい関数を考えましょう。……とその前に、「画像を読み込んで色味を(R,G,B)で算出する」関数を作成しておきます。

```python
from typing import Tuple
from PIL import Image
from PIL.Image import BOX
from PIL.MpoImagePlugin import MpoImageFile

def get_sample_color(path: str) -> Tuple[int, int, int]:
    """指定したパスの画像ファイルを読み込み、縦横1/10の空間を中央から切り取り、その平均色を取得して返す

    :param path: ファイル名
    :return: (R, G, B)
    """
    im: MpoImageFile = Image.open(path)
    image_size = im.size
    image_size2 = (image_size[0] / 10, image_size[1] / 10)
    left_upper_pos = ((image_size[0] - image_size2[0]) / 2, (image_size[1] - image_size2[1]) / 2)
    cropped_im = im.crop((left_upper_pos[0], left_upper_pos[1],
                          left_upper_pos[0] + image_size2[0], left_upper_pos[1] + image_size2[1]))
    resized_im = cropped_im.resize((1, 1), resample=BOX)
    data_r = resized_im.getcolors()[0][1][0]
    data_g = resized_im.getcolors()[0][1][1]
    data_b = resized_im.getcolors()[0][1][2]
    return data_r, data_g, data_b
```

また、色味の差を算出する関数も必要ですね。RGB式空間ではなく、YCbCr色空間に変換し、CbとCrについてL2ノルムを計算します。

```python
def calc_color_diff(color1: Tuple[int, int, int], color2: Tuple[int, int, int]) -> float:
    """2つの色の差を計算する

    :param color1: 色1
    :param color2: 色2
    :return: YCrCr空間における2つの色の差
    """
    # JPEGのYCbCr式のdiff
    cb_1 = -0.1687 * color1[0] - 0.3313 * color1[1] + 0.5 * color1[2] + 128
    cr_1 = 0.5 * color1[0] - 0.4187 * color1[1] - 0.0813 * color1[2] + 128
    cb_2 = -0.1687 * color2[0] - 0.3313 * color2[1] + 0.5 * color2[2] + 128
    cr_2 = 0.5 * color2[0] - 0.4187 * color2[1] - 0.0813 * color2[2] + 128
    return (cb_1 - cb_2) ** 2 + (cr_1 - cr_2) ** 2
```

そして最適化したい関数ですが、一般的なPyhon関数と異なり、
「写真を撮影して読み込める状態になるまで**処理を一旦止める**」処理が必要です。
……一番手っ取り早いのは`input()`関数だと思ったので、そちらを使って実装しました。

```python
import os

def get_score(a: float, g: float) -> float:
    """指定した設定におけるスコア(＝基準となる色味との間の距離)を計算する

    :param a: パラメーターA(アンバー)の値。負数だとB(ブルー)側の色味になる
    :param g: パラメーターG(グリーン)の値。負数だとM(マゼンタ)側の色味になる
    :return: 色味の間の距離。小さいほど良い
    """

    # 読み込みたいファイルの名称を生成する
    if a >= 0.0:
        temp1 = f'A{int(a)}'
    else:
        temp1 = f'B{int(-a)}'
    if g >= 0.0:
        temp2 = f'G{int(g)}'
    else:
        temp2 = f'M{int(-g)}'
    file_name = f'MZD_{temp1}{temp2}.jpg'

    # ファイルが存在しないかを確認し、存在しないなら作成させる
    while not os.path.exists(file_name):
        _ = input(f'{file_name}を作成したら、Enterキーを押してください.')

    # ファイルを読み込み、中央付近の平均色(＝測定値)を取得する
    mzd_rgb = get_sample_color(file_name)

    # サンプルデータの色と比較する
    score = calc_color_diff(mzd_rgb, base_rgb)
    print(f'a={a} g={g} score={score}')
    return score
```

後の流れは大体同じですね。注意したいのは、今回の入力(ホワイトバランス設定値)は離散値なので、離散値用に`'type': 'discrete'`としておく必要があることです。まあ前述のように、`'type': 'discrete'`でも整数値(`int型`)が最適化したい関数に渡されるわけではなく、常に`float`型なのですが。

```python
from numpy.core.multiarray import ndarray

def get_score2(x: ndarray) -> float:
    """指定した設定におけるスコア(＝基準となる色味との間の距離)を計算する"""
    return get_score(x[:, 0][0], x[:, 1][0])


if __name__ == '__main__':
    bounds = [{'name': 'a', 'type': 'discrete', 'domain': tuple(range(-9, 10))},
              {'name': 'g', 'type': 'discrete', 'domain': tuple(range(-9, 10))}]

    myProblem = GPyOpt.methods.BayesianOptimization(get_score2,
                                                    domain=bounds,
                                                    maximize=False,
                                                    verbosity=True)
    max_iter = 25
    myProblem.run_optimization(max_iter, verbosity=True)
    myProblem.plot_convergence()
    print(myProblem.x_opt, myProblem.fx_opt)
    myProblem.plot_acquisition()

```

# 最適化結果

おおよそこんな感じ。評価値が若干ベタっとしていたので収束性は悪いですが、概ね満足できる結果が得られました。

![Figure_1.jpg](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/29802/554b999a-d8e5-de5c-6d53-3dbc5534cdea.jpeg)
![Figure_2.jpg](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/29802/b98dabd3-61f9-a5f6-fa55-72cc826c31a5.jpeg)

# さらに改造するなら

「色味」というのはホワイトバランスだけではなく、絵のコントラストも重要です。
例えば、「LUMIXブランドかLEICAブランドかでコントラストが違う……」と感じる方もいるそうですし。

今回は「白っぽいボカした壁(≒18％グレー)」で調整したのでホワイトバランスしか調整しようがありませんが、
「段階的な白～黒の色が印刷された物体(グレーチャート)」で調整することにより、
中間色の度合い……つまりコントラストについてもチューニングが可能になります。

ベイズ最適化はもちろん3次元以上の入力にも対応していますので、
より厳密な調整を行いたい場合は、上記プログラムを改造してみるといいでしょう。