最適化アルゴリズムを実装していくぞ（カッコウ探索）


# はじめに

最適化アルゴリズムの実装シリーズです。
まずは[概要](https://qiita.com/pocokhc/items/07b698cc426cadb3a64e)を見てください。

コードは[github](https://github.com/pocokhc/metaheuristics)にあります。



# カッコウ探索
## 概要

カッコウ探索(Cucko Search)は、カッコウの托卵という行動をもとに作られたアルゴリズムです。
托卵は他の種類の鳥の巣に卵を産み、その巣の親に育成を頼む行為です。

また、カッコウを含めた野生動物は獲物を探すときにレヴィフライトというレヴィ分布と呼ばれる分布に従って移動することが知られています。


参考
・[Cuckoo search](https://en.wikipedia.org/wiki/Cuckoo_search)
・[進化計算アルゴリズム入門 生物の行動科学から導く最適解](https://www.amazon.co.jp/%E9%80%B2%E5%8C%96%E8%A8%88%E7%AE%97%E3%82%A2%E3%83%AB%E3%82%B4%E3%83%AA%E3%82%BA%E3%83%A0%E5%85%A5%E9%96%80-%E7%94%9F%E7%89%A9%E3%81%AE%E8%A1%8C%E5%8B%95%E7%A7%91%E5%AD%A6%E3%81%8B%E3%82%89%E5%B0%8E%E3%81%8F%E6%9C%80%E9%81%A9%E8%A7%A3-%E5%A4%A7%E8%B0%B7%E7%B4%80%E5%AD%90-ebook/dp/B07DQFVK1H)


## アルゴリズム

ランダムな巣を元に卵を作成し、別の巣を作ります。
この時に作成される卵は、レヴィ分布に従って生成されます。

作成された巣はさらに別のランダムな巣と比較し、良ければ置き換えられます。

最後に評価値の悪い巣は破棄され、新しい巣に置き換えられます。


* アルゴリズムのフロー

![draw2-Cuckoo.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/396150/3d3590ad-e1b0-b94e-5261-ae3df48ac04e.png)


* 用語の対応

|問題|カッコウ探索|
|---|---|
|入力値の配列|巣|
|入力値|卵|
|評価値|巣の評価値|



* ハイパーパラメータに関して

|変数名|意味|所感|
|---|---|---|
|scaling_rate|レヴィ分布の尺度|大きいほど遠くに移動しやすくなる|
|levy_rate|レヴィフライトの反映率|大きいほどレヴィフライトによる移動が大きい|
|bad_nest_rate|悪い巣の割合|高いほど巣の入れ替え（探索範囲）が多くなる|



## 巣の作成

ランダムに選んだ巣からレヴィ分布に従って新しい巣を作成します。

$$ x^{new}_i = x^{r}_i + \alpha s $$

$s$ はレヴィ分布に従った乱数（レヴィフライト）で、$\alpha$ はレヴィフライトの反映率となります。
この巣の生成がカッコウ探索のすべてですが、レヴィフライトがちょっと厄介です。

レヴィ分布は以下です。

$$ f(x; \mu, c) = \sqrt{ \frac{c}{2 \pi} } \frac{ \exp^{\frac{-c}{2(x - \mu)}} }{ (x - \mu)^{\frac{2}{3}}} $$

``` python
import math
def levy(x, u=0, c=1):
    if x == 0:
        return 0
    t = math.exp((-c/(2 * (x-u))))
    t /= (x-u) ** (3/2)
    return math.sqrt(c/(2*math.pi)) * t
```

このレヴィ分布に従った乱数を生成する必要があります。
レヴィ分布に従った乱数は Mantegna アルゴリズムで求めます。

$$ s = \frac{u}{|v|^{\frac{1}{\beta}}} $$

ここで $\beta$ はスケーリング指数で0～2の実数をとります。
$u$ は平均0,分散 $\sigma^2$の正規分布に従う乱数、
$v$ は標準正規分布に従う乱数です。

$\sigma$ は以下の式で求まります。

$$ \sigma = \Biggl( \frac{ \Gamma(1+\beta) \sin(\frac{\pi \beta}{2}) }{ \Gamma(\frac{1+\beta}{2}) \beta 2^{\frac{(\beta-1)}{2}} } \Biggr)^{\frac{1}{\beta}} $$

$\Gamma$ はガンマ関数を表します。

$$ \Gamma(x) = \int_{0}^{\infty} t^{x-1}e^{-t} dt $$

* 正規分布の乱数

正規分布の乱数ですが、pythonのnumpyライブラリを使うと以下で簡単に出せます。

``` python
import numpy as np
np.random.randn()  # 標準正規分布の一様乱数を生成
np.random.normal(0, sigma)  # 平均0、分散sigma^2の正規分布に従う乱数
```

が、他の言語(主にluaとか)での実装も考えているのでライブラリを使わない方法も書いておきます。

標準正規分布の乱数はボックス＝ミュラー法で出すことができます。

$$ z = \sqrt{-2 \log{X}} \cos{2\pi Y} $$ 

ここで X と Y は互いに独立した一様乱数になります。
pythonだと以下です。

``` python
import math
def random_normal():
    r1 = random.random()
    r2 = random.random()
    return math.sqrt(-2.0 * math.log(r1)) * math.cos(2*math.pi*r2)
```


* ガンマ関数

ガンマ関数も math を使うと簡単に出せます。

``` python
import math
math.gamma(x)
```

使わない場合は[ここのサイト](https://www.sist.ac.jp/~suganuma/programming/9-sho/prob/gamma/gamma.htm)のコードを使わせてもらいました。
（コードは全体コードに記載しています）


## レヴィ分布に関して

レヴィ分布のグラフです。

![plot_levy.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/396150/682eba4b-065d-e4e8-15cd-7afa4dfbffb9.png)

レヴィ分布ですが見てわかるように範囲が 0～∞ です。
これは卵を生成する場合に少しネックになります。

例えばOneMaxの問題だと範囲が0～1しかとりません。
この問題に対してレヴィ分布で生成された1以上の値をどうするかは決まっていません。
ここで、例えば1以上の値は1にするとした場合、生成された乱数がレヴィ分布に従わなくなります。
(1の生成確率だけ高くなってしまう)

これはどうなんでしょうかね。
解決方法もよくわからないので今のところは放置しています。



## コード全体

``` python
import math
import random

import numpy as np


############################################
# Γ（ｘ）の計算（ガンマ関数，近似式）
#      ier : =0 : normal
#            =-1 : x=-n (n=0,1,2,･･･)
#      return : 結果
#      coded by Y.Suganuma
# https://www.sist.ac.jp/~suganuma/programming/9-sho/prob/gamma/gamma.htm
############################################
def gamma(x):
    if x <= 0:
        raise ValueError("math domain error")

    ier = 0

    if x > 5.0 :
        v = 1.0 / x
        s = ((((((-0.000592166437354 * v + 0.0000697281375837) * v + 0.00078403922172) * v - 0.000229472093621) * v - 0.00268132716049) * v + 0.00347222222222) * v + 0.0833333333333) * v + 1.0
        g = 2.506628274631001 * math.exp(-x) * pow(x,x-0.5) * s

    else:

        err = 1.0e-20
        w   = x
        t   = 1.0

        if x < 1.5 :

            if x < err :
                k = int(x)
                y = float(k) - x
                if abs(y) < err or abs(1.0-y) < err :
                    ier = -1

            if ier == 0 :
                while w < 1.5 :
                    t /= w
                    w += 1.0

        else :
            if w > 2.5 :
                while w > 2.5 :
                    w -= 1.0
                    t *= w

        w -= 2.0
        g  = (((((((0.0021385778 * w - 0.0034961289) * w + 0.0122995771) * w - 0.00012513767) * w + 0.0740648982) * w + 0.0815652323) * w + 0.411849671) * w + 0.422784604) * w + 0.999999926
        g *= t

    return g


def random_normal():
    """
    正規分布の乱数
    ボックス＝ミュラー法
    """
    r1 = random.random()
    r2 = random.random()
    return math.sqrt(-2.0 * math.log(r1)) * math.cos(2*math.pi*r2)


def mantegna(beta):
    """
    mantegna アルゴリズム
    """
    #beta:  0.0 - 2.0
    if beta < 0.005:
        # 低すぎると OverflowError: (34, 'Result too large')
        beta = 0.005
    # siguma
    t = gamma(1+beta) * math.sin(math.pi*beta/2)
    t = t/( gamma((1+beta)/2) * beta * 2**((beta-1)/2) )
    siguma = t**(1/beta)

    u = random_normal()*siguma  # 平均0 分散siguma^2 の正規分布に従う乱数
    v = random_normal()  # 標準正規分布に従う乱数

    s = (abs(v)**(1/beta))
    if s < 0.0001:
        # 低すぎると ValueError: supplied range of [-inf, inf] is not finite
        s = 0.0001
    s = u / s
    return s


class Cuckoo():
    def __init__(self, 
        nest_max,
        scaling_rate=1.0,
        levy_rate=1.0,
        bad_nest_rate=0.1
    ):
        self.nest_max = nest_max
        self.scaling_rate = scaling_rate
        self.levy_rate = levy_rate

        # 悪い巣の割合から悪い巣の個数を算出
        self.bad_nest_num = int(nest_max * bad_nest_rate + 0.5)
        if self.bad_nest_num > nest_max-1:
            self.bad_nest_num = nest_max-1
        if self.bad_nest_num < 0:
            self.bad_nest_num = 0

    def init(self, problem):
        self.problem = problem

        self.nests = []
        for _ in range(self.nest_max):
            self.nests.append(problem.create())
    
    def step(self):
        # ランダムに巣を選択
        r = random.randint(0, self.nest_max-1)  # a<=x<=b

        # 新しい巣を作成
        arr = self.nests[r].getArray()

        for i in range(len(arr)):

            # レヴィフライで卵を作る
            arr[i] = arr[i] + self.levy_rate * mantegna(self.scaling_rate)

        new_nest = self.problem.create(arr)

        # ランダムな巣と比べてよければ変える
        r = random.randint(0, self.nest_max-1)  # a<=x<=b
        if self.nests[r].getScore() < new_nest.getScore():
            self.nests[r] = new_nest

        # 悪い巣を消して新しく作る
        self.nests.sort(key=lambda x:x.getScore())
        for i in range(self.bad_nest_num):
            self.nests[i] = self.problem.create()

```

# カッコウ探索(ε-greedy)

これは本記事オリジナルです。（探せばあるかも？）
カッコウ探索はとてもシンプルなアルゴリズムですが、レヴィフライトの実装がすごく厄介です。

そこでレヴィフライトをε-greedyに置き換えてみました。
ε-greedyに置き換えることで実装がすごい簡単になりましたね。

ε-greedyでもそこそこの精度がでるような気がします。


## コード全体

``` python
import math
import random

class Cuckoo_greedy():
    def __init__(self, 
        nest_max,
        epsilon=0.1, 
        bad_nest_rate=0.1
    ):
        self.nest_max = nest_max
        self.epsilon = epsilon

        # 悪い巣の割合から悪い巣の個数を算出
        self.bad_nest_num = int(nest_max * bad_nest_rate + 0.5)
        if self.bad_nest_num > nest_max-1:
            self.bad_nest_num = nest_max-1
        if self.bad_nest_num < 0:
            self.bad_nest_num = 0

    def init(self, problem):
        self.problem = problem

        self.nests = []
        for _ in range(self.nest_max):
            self.nests.append(problem.create())
    
    def step(self):
        # ランダムに巣を選択
        r = random.randint(0, self.nest_max-1)  # a<=x<=b

        # 新しい巣を作成
        arr = self.nests[r].getArray()

        for i in range(len(arr)):

            # ε-greedy で卵を新しく作成する
            if random.random() < self.epsilon:
                arr[i] = self.problem.randomVal()

        new_nest = self.problem.create(arr)

        # ランダムな巣と比べてよければ変える
        r = random.randint(0, self.nest_max-1)  # a<=x<=b
        if self.nests[r].getScore() < new_nest.getScore():
            self.nests[r] = new_nest

        # 悪い巣を消して新しく作る
        self.nests.sort(key=lambda x:x.getScore())
        for i in range(self.bad_nest_num):
            self.nests[i] = self.problem.create()

```



# ハイパーパラメータ例

各問題に対して optuna でハイパーパラメータを最適化した結果です。
最適化の1回の試行は、探索時間を2秒間として結果を出しています。
これを100回実行し、最適なハイパーパラメータを optuna に探してもらいました。

* カッコウ探索

|問題|bad_nest_rate|levy_rate|nest_max|scaling_rate|
|---|---|---|---|---|
|EightQueen|0.09501642206708413|0.9797131483689493|14|1.9939515457735189|
|function_Ackley|0.0006558326608885681|0.3538825414958845|4|0.9448539685962172|
|function_Griewank|0.23551408245457767|0.30150681160121073|2|0.9029863706820189|
|function_Michalewicz|0.00438839398648697|0.0004796264527609298|2|1.5288609934193742|
|function_Rastrigin|0.13347040982335695|0.031401149135082206|7|1.6949622109706082|
|function_Schwefel|0.0003926596935418525|0.02640034426449156|4|0.5809451877075759|
|function_StyblinskiTang|0.08462936367613791|0.0633939067767827|5|1.7236388666366773|
|LifeGame|0.8819375718376719|0.015175414454036936|33|1.3899842408715666|
|OneMax|0.89872646833605|0.1261650035421213|17|0.04906594355889626|
|TSP|0.024559598255857823|0.008225444982304852|4|1.8452535160497248|


* カッコウ探索(ε-greedy)

|問題|bad_nest_rate|epsilon|nest_max|
|---|---|---|---|
|EightQueen|0.004374125594794304|0.03687227169502155|7|
|function_Ackley|0.5782260075661492|0.031195954391595435|2|
|function_Griewank|0.23314007403872794|0.05206930732996057|2|
|function_Michalewicz|0.11845570554906226|0.02242832420874199|3|
|function_Rastrigin|0.009725819291390304|0.025727770986639094|3|
|function_Schwefel|0.22978641596753258|0.048159183280607774|2|
|function_StyblinskiTang|0.14184473157004032|0.01965829867603547|2|
|LifeGame|0.7358005558643367|0.9115290938258255|39|
|OneMax|0.0016700608620328905|0.006003869128710593|2|
|TSP|0.00023997215188062415|0.030790166824531992|29|


# 実際の動きの可視化

1次元は6個体、2次元は20個体で50step実行した結果です。
赤い丸がそのstepでの最高スコアを持っている個体となります。

パラメータは以下で実行しました。

``` python
Cuckoo(N, scaling_rate=1.0, levy_rate=1.0, bad_nest_rate=0.1)
Cuckoo_greedy(N, epsilon=0.5, bad_nest_rate=0.1)
```




## function_Ackley

* Cuckoo
  * 1次元

![function_Ackley_Cuckoo_2.gif](https://github.com/pocokhc/metaheuristics/blob/main/imgs/function_Ackley_Cuckoo_2.gif?raw=true)


  * 2次元

![function_Ackley_Cuckoo_3.gif](https://github.com/pocokhc/metaheuristics/blob/main/imgs/function_Ackley_Cuckoo_3.gif?raw=true)


* Cuckoo_greedy
  * 1次元

![function_Ackley_Cuckoo_greedy_2.gif](https://github.com/pocokhc/metaheuristics/blob/main/imgs/function_Ackley_Cuckoo_greedy_2.gif?raw=true)


  * 2次元

![function_Ackley_Cuckoo_greedy_3.gif](https://github.com/pocokhc/metaheuristics/blob/main/imgs/function_Ackley_Cuckoo_greedy_3.gif?raw=true)






## function_Rastrigin


* Cuckoo
  * 1次元

![ffunction_Rastrigin_Cuckoo_2.gif](https://github.com/pocokhc/metaheuristics/blob/main/imgs/function_Rastrigin_Cuckoo_2.gif?raw=true)


  * 2次元

![function_Rastrigin_Cuckoo_3.gif](https://github.com/pocokhc/metaheuristics/blob/main/imgs/function_Rastrigin_Cuckoo_3.gif?raw=true)



* Cuckoo_greedy
  * 1次元

![function_Rastrigin_Cuckoo_greedy_2.gif](https://github.com/pocokhc/metaheuristics/blob/main/imgs/function_Rastrigin_Cuckoo_greedy_2.gif?raw=true)


  * 2次元

![function_Rastrigin_Cuckoo_greedy_3.gif](https://github.com/pocokhc/metaheuristics/blob/main/imgs/function_Rastrigin_Cuckoo_greedy_3.gif?raw=true)



## function_Schwefel

* Cuckoo
  * 1次元

![function_Schwefel_Cuckoo_2.gif](https://github.com/pocokhc/metaheuristics/blob/main/imgs/function_Schwefel_Cuckoo_2.gif?raw=true)


  * 2次元

![function_Schwefel_Cuckoo_3.gif](https://github.com/pocokhc/metaheuristics/blob/main/imgs/function_Schwefel_Cuckoo_3.gif?raw=true)



* Cuckoo_greedy
  * 1次元

![function_Schwefel_Cuckoo_greedy_2.gif](https://github.com/pocokhc/metaheuristics/blob/main/imgs/function_Schwefel_Cuckoo_greedy_2.gif?raw=true)


  * 2次元

![function_Schwefel_Cuckoo_greedy_3.gif](https://github.com/pocokhc/metaheuristics/blob/main/imgs/function_Schwefel_Cuckoo_greedy_3.gif?raw=true)



## function_StyblinskiTang


* Cuckoo
  * 1次元

![function_StyblinskiTang_Cuckoo_2.gif](https://github.com/pocokhc/metaheuristics/blob/main/imgs/function_StyblinskiTang_Cuckoo_2.gif?raw=true)


  * 2次元

![function_StyblinskiTang_Cuckoo_3.gif](https://github.com/pocokhc/metaheuristics/blob/main/imgs/function_StyblinskiTang_Cuckoo_3.gif?raw=true)



* Cuckoo_greedy
  * 1次元

![function_StyblinskiTang_Cuckoo_greedy_2.gif](https://github.com/pocokhc/metaheuristics/blob/main/imgs/function_StyblinskiTang_Cuckoo_greedy_2.gif?raw=true)


  * 2次元

![function_StyblinskiTang_Cuckoo_greedy_3.gif](https://github.com/pocokhc/metaheuristics/blob/main/imgs/function_StyblinskiTang_Cuckoo_greedy_3.gif?raw=true)



## function_XinSheYang

* Cuckoo
  * 1次元

![function_XinSheYang_Cuckoo_2.gif](https://github.com/pocokhc/metaheuristics/blob/main/imgs/function_XinSheYang_Cuckoo_2.gif?raw=true)


  * 2次元

![function_XinSheYang_Cuckoo_3.gif](https://github.com/pocokhc/metaheuristics/blob/main/imgs/function_XinSheYang_Cuckoo_3.gif?raw=true)



* Cuckoo_greedy
  * 1次元

![function_XinSheYang_Cuckoo_greedy_2.gif](https://github.com/pocokhc/metaheuristics/blob/main/imgs/function_XinSheYang_Cuckoo_greedy_2.gif?raw=true)


  * 2次元

![function_XinSheYang_Cuckoo_greedy_3.gif](https://github.com/pocokhc/metaheuristics/blob/main/imgs/function_XinSheYang_Cuckoo_greedy_3.gif?raw=true)




# あとがき

アルゴリズム自体はかなり簡単で精度もかなりいいらしいです。
レヴィ分布に従った乱数の生成が一番難しかったです…。