巡回セールスマン問題を遺伝的アルゴリズムを使ってHaskellで解く！

大学の授業で巡回セールスマン問題を解く課題があったので、Haskellでこの問題を解いてみました！

# そもそも巡回セールスマン問題とは？

複数の都市があって、セールスマンはその場所を一回ずつ通ります。最短経路で全ての場所を通って戻ってくるにはどのルートを通ればいいですか、という問題です。

[巡回セールスマン問題](https://www.wikiwand.com/ja/%E5%B7%A1%E5%9B%9E%E3%82%BB%E3%83%BC%E3%83%AB%E3%82%B9%E3%83%9E%E3%83%B3%E5%95%8F%E9%A1%8C)

# まずデータ構造を定義する

入力データ（場所の情報）は、x座標とy座標です。例えば、１つ目の場所は座標（90, 100）にあります。

```input.txt
90 100 
40 15 
30 30 
```

## 座標の型を定義する (Coord)

座標が出てきたので、これを型にしてしまいましょう。

```Haskell
newtype Coord = Coord (Double, Double) deriving (Eq)
```

２点間の距離、例えば（90, 100）、（40, 15）の距離を図りたいときは、ユークリッド距離を使います。

Haskellでかくとこのような感じ。

```Haskell
distance :: Coord -> Coord -> Double
distance (Coord (x1, y1)) (Coord (x2, y2)) = sqrt $ (x1 - x2)^2 + (y1 - y2)^2
```

## 巡回路を定義する (Path)

また、座標を複数格納できるPath型も定義しましょう。

```Haskell
import Data.Vector
newtype Path = Vector Coord
```

ここではリストではなく、Vectorを使うことにします。

また、PathとPathを比較できるように、Ord型クラスのインスタンスにしましょう。

```Haskell
instance Ord Path where
  compare p1 p2
    | totalCost p1 > totalCost p2 = GT
    | totalCost p1 == totalCost p2 = EQ
    | totalCost p1 < totalCost p2 = LT
```

## 合計の距離の計算をどうするか？

例えば、3地点を(90,100), (40,15), (30,30)という順番で通るとします。そのときの移動距離は、

```
(90,100)と(40,15)の距離 + (40,15)と(30,30)の距離  + (30,30)と(90,100)の距離
```

となります。そこで、totalCostという関数を作ります。

```Haskell
eachTuple :: Vector Coord -> Vector (Coord, Coord)
eachTuple path =
  let tuples = do
        index <- enumFromTo 0 (length path - 2)
        return (path ! index, path ! (index + 1))
  in tuples `snoc` (last path, head path)

totalCost :: Path -> Double
totalCost (Path path) = sum $ uncurry distance <$> eachTuple path
```

eachTupleで、２点間の座標をタプルにしています。つまり、`[ ((90,100), (40,15)), ((40,15), (30,30)), ((30,30), (90,100)) ]`のようになります。それをdistanceにそれぞれ適用させて、あとはそれらを合計しています。

あとは、これらを`TSP.hs`にまとめてしまいましょう。(TSPは、traveling salesman problemの略です。)

```Haskell:TSP.hs
module TSP (
  totalCost
, Path(..)
, Coord
) where

import Data.Vector
import Prelude hiding (
  enumFromTo,
  length,
  last,
  head,
  sum
  )

newtype Coord = Coord (Double, Double) deriving (Eq, Show, Ord)
newtype Path = Path (Vector Coord) deriving (Eq, Show)

instance Ord Path where
  compare p1 p2
    | totalCost p1 > totalCost p2 = GT
    | totalCost p1 == totalCost p2 = EQ
    | totalCost p1 < totalCost p2 = LT

eachTuple :: Vector Coord -> Vector (Coord, Coord)
eachTuple path =
  let tuples = do
        index <- enumFromTo 0 (length path - 2)
        return (path ! index, path ! (index + 1))
  in tuples `snoc` (last path, head path)

totalCost :: Path -> Double
totalCost (Path path) = sum $ uncurry distance <$> eachTuple path

distance :: Coord -> Coord -> Double
distance (Coord (x1, y1)) (Coord (x2, y2)) = sqrt $ (x1 - x2)^2 + (y1 - y2)^2
```

# どのアルゴリズムで解くの？

今回は**遺伝的アルゴリズム**で解きます！


# 遺伝的アルゴリズム
簡単に言えば、生物がどのようにして進化したかをアルゴリズムにしたものです。詳しくは↓
[遺伝的アルゴリズム](https://www.wikiwand.com/ja/%E9%81%BA%E4%BC%9D%E7%9A%84%E3%82%A2%E3%83%AB%E3%82%B4%E3%83%AA%E3%82%BA%E3%83%A0)

遺伝的アルゴリズムは、以下のアルゴリズムで構成されています。

1. N個の個体(巡回路)をランダムに初期化
2. 親を選んで交叉する (crossover)
3. 低確率で突然変異 (mutation)
4. 次の世代の個体（巡回路）を選択する (selection)

今回は、ランダムな要素も入っているので、それぞれの関数を`IO`にしました。遺伝的アルゴリズムの全体像はこんな感じです。
また、最良解を保存するために「エリート保存」をしています。

```Haskell:GA.hs
mating :: Vector Path -> IO (Vector Path)
mating = undefined

mutation :: Vector Path -> IO (Vector Path)
mutation = undefined

selectChildren :: Vector Path -> IO (Vector Path)
selectChildren = undefined

geneticLoop :: Int -> Vector Path -> Path -> IO Path
geneticLoop 0 _ elite = return elite
geneticLoop generation parents elite = do

  newParents <- mating parents >>=
                mutation >>=
                selectChildren

  let minPath = minimum newParents :: Path
      newElite = min minPath elite :: Path

  printf "%d generation , %f , %f\n" generation (totalCost newElite) (totalCost minPath)

  geneticLoop (generation - 1) newParents newElite
```

## 1. N個の個体(巡回路)をランダムに初期化

ここでは、巡回路をシャッフルすることで初期化します。Vectorにはshuffle関数はないようなので、自分で実装しました。

型はこのような感じにしました。

```Haskell
shuffleM :: (Eq a) => Vector a -> IO (Vector a)
shuffleM x = undefined
```

ここでは、`a`が`Coord`型になります。

shuffleの実装は、ある要素を１つずつランダムにポップして、それらを連結するかんじです。

例えば以下のようにシャッフルします。

1. `[0,1,2,3]`をシャッフル
2. ランダムにインデックスを選ぶ（`1`のインデックスが選ばれたとする）
3. `[0,1,2,3]`から`１`をポップして、その`1`を結果に入れる。
4. 新たに`[0,2,3]`をシャッフルする

これらを繰り返します。

`Shuffle`モジュールを作ったので、それを載せておきます。

```Haskell:Shuffle.hs
module Shuffle (
  shufflePath,
  shuffle,
  shuffleM
) where

import TSP
import System.Random
import Data.Vector
import Prelude hiding (head, tail, length)

shufflePath :: Path -> StdGen -> Path
shufflePath (Path coords) gen = Path $ shuffle coords gen

shuffle :: (Eq a) => Vector a -> StdGen -> Vector a
shuffle x gen
  | x == empty = empty
  | otherwise  =
      let (index,newGen) = randomR (0, length x - 1) gen :: (Int, StdGen)
          targetElm      = x ! index
      in targetElm `cons` shuffle (deleteFirst targetElm x) newGen

shuffleM :: (Eq a) => Vector a -> IO (Vector a)
shuffleM x = do
  gen <- newStdGen
  return $ shuffle x gen

-- | 削除したい要素の中で最初に出てきたものを取り除く
-- >>> deleteFirst 2 $ fromList [1..5]
-- [1,3,4,5]
-- >>> deleteFirst 1 $ fromList [1,2,3,4,1]
-- [2,3,4,1]
deleteFirst :: (Eq a) => a -> Vector a -> Vector a
deleteFirst x vector
  | x == head vector = tail vector
  | otherwise        = head vector `cons` deleteFirst x (tail vector)
```

あとは、以下のようにして使いましょう。

```Haskell:Main.hs
import Shuffle (shuffleM)
import TSP

inputPath :: IO (Vector Coord)
inputPath = do
  contents <- getContents
  return $ fromList $ toCoord . words <$> lines contents
    where
      toCoord :: [String] -> Coord
      toCoord [s1, s2] = Coord (read s1, read s2)
      toCoord _ = error "it is illegal file format!"

main :: IO ()
main = do
  coords <- inputPath
  paths <- replicateM poolSize $ Path <$> shuffleM coords
```

## 2. 親を選んで交叉する

交叉とは、親の遺伝子が交わることを指します。人間も同じように親の遺伝子を受け継いでいますよね。
具体的には、以下のようにしました。

1. ２つの巡回路（親）をランダムに選択
2. 順序交叉（order crossover）

### ２つの巡回路（親）をランダムに選択

まず、複数の巡回路から2つの巡回路をランダムに選ぶ関数を作りましょう。

```Haskell
randomRSt :: (RandomGen g, Random a) => (a, a) -> State g a
randomRSt (a,b) = state $ randomR (a,b)

twoRandomNumState :: (Random a) => (a,a) -> State StdGen (a, a)
twoRandomNumState fromTo = do
  n1 <- randomRSt fromTo
  n2 <- randomRSt fromTo
  return (n1, n2)

twoRandomNum :: (Random a) => (a,a) -> StdGen -> (a,a)
twoRandomNum fromTo = evalState (twoRandomNumState fromTo)

randomSelect :: Vector Path -> IO (Path, Path)
randomSelect paths = do
  gen <- newStdGen
  let n = length paths
      (r1, r2) = twoRandomNum (0, n - 1) gen
  return (paths ! r1, paths ! r2)
```

2つのランダムな数字を安全に生成するために（誤って同じ乱数ジェネレータを使わないように）、`State`モナドを使いました。`twoRandomNum`関数は、ある範囲の中で２つの乱数を作ります。

### 順序交叉（order crossover）

[順序交叉](http://www.sist.ac.jp/~kanakubo/research/evolutionary_computing/genetic_algorithms.html)

さきほど選んだ２つを交叉させます。順序交叉は以下のように交叉します。

![Screenshot from 2017-12-04 10-23-05.png](https://qiita-image-store.s3.amazonaws.com/0/136124/f3155e49-b638-35cf-a913-9e7ee9a79343.png)


1. ランダムに交叉点を決める
2. 子１は、親１から左側をそのまま受け継ぎ、親２からは「親１から引き継いだものを除いたもの」を受け継ぐ。
3. 同様に、子2は、親２から左側をそのまま受け継ぎ、親２からは「親１から引き継いだものを除いたもの」を受け継ぐ。

といった感じです。

```Haskell
orderCrossOver :: Path -> Path -> IO (Vector Path)
orderCrossOver papa@(Path coords) mama = do
  gen <- newStdGen
  let n = length coords
      (dividePoint,_) = randomR (1, n - 1) gen :: (Int, StdGen)
  return $ divCpChromo dividePoint mama papa

divCpChromo :: Int -> Path -> Path -> Vector Path
divCpChromo divPoint (Path papa) (Path mama) =
  let (p1, p2) = splitAt divPoint papa :: (Vector Coord, Vector Coord)
      m1 = filter (`notElem` p1) mama :: Vector Coord
      m2 = filter (`notElem` p2) mama :: Vector Coord
  in Path (p1 ++ m1) `cons` singleton (Path (m2 ++ p2))
```

交叉点をインデックスの1からn-1としているのは、０からだと何も入れ替わらないからです。

あとは、今説明した処理を繰り返せばいいだけです。

```Haskell
mating :: Vector Path -> IO (Vector Path)
mating parents = do
  childList2 <- replicateM (length parents) $ do
    (papa, mama) <- matingSelection parents
    crossOver papa mama :: IO (Vector Path)
  return $ concatVector childList2

concatVector :: Vector (Vector a) -> Vector a
concatVector = concatMap id
```

これでmating関数ができあがりました。

## 3. 低確率で突然変異

今回は１０％の確率で、巡回路がシャッフルされます。Vectorの要素を１つずつ確認するような処理なので、foldlを使ってみました。

```Haskell
mutationRate = 0.1

mutation :: Vector Path -> IO (Vector Path)
mutation children = do
  gen <- newStdGen
  return $ fst $ foldl mutate (empty, gen) children
  where
    mutate :: (Vector Path, StdGen) -> Path -> (Vector Path, StdGen)
    mutate (ps, gen) p = do
      let (r, newG) = randomR (0, 1) gen :: (Double, StdGen)
      if r < mutationRate
        then (shufflePath p newG `cons` ps, newG)
        else (p `cons` ps, newG)
```

## 4. 次の世代の個体（巡回路）を選択する (selection)

今回は、トーナメント選択を使います。

![Screenshot from 2017-12-04 10-56-37.png](https://qiita-image-store.s3.amazonaws.com/0/136124/af4b8b2b-7bb7-4fc6-4604-85ef2fb72f87.png)

1. ランダムにｎ個選ぶ
2. その中からcostが最小のものを選ぶ
3. 指定した分だけ繰り返す

まず、複数の巡回路から１つの巡回路をランダムに選択する関数を作りましょう。

```Haskell
randomSelectOne :: Vector Path -> IO Path
randomSelectOne paths = do
  gen <- newStdGen
  let n = length paths
      (index, _) = randomR (0, n - 1) gen :: (Int, StdGen)
  return $ paths ! index
```

そして、トーナメント選択をする関数を作ります。

```Haskell
tournamentSelect :: Int -> Int -> Vector Path -> IO (Vector Path)
tournamentSelect poolSize n paths = replicateM poolSize $ do
  randomSelectPath <- replicateM n (randomSelectOne paths)
  return $ minimum randomSelectPath
```

ここでの`n`は、トーナメントの大きさです。上の図では`n=3`になっています。
`poolSize`は、「指定した分だけ繰り返す」ものです。

## 全体のコード

あとは、main関数で`geneticLoop`関数を呼べばいいですね。

```Haskell
generationNum = 100

main = do
  coords <- inputPath

  paths <- replicateM poolSize $ Path <$> shuffleM coords
  elite <- geneticLoop generationNum paths (head paths)

  print $ totalCost elite
```

Githubに全体のコードを載せました。
[Genetic-Algorithm-Haskell](https://github.com/v97ug/Genetic-Algorithm-Haskell)

# 実験

[TSP Data for the Traveling Salesperson Problem](https://people.sc.fsu.edu/~jburkardt/datasets/tsp/tsp.html)からデータを取りました。[att48.tsp](https://people.sc.fsu.edu/~jburkardt/datasets/tsp/att48_xy.txt)というデータセットを使います。

レポジトリをダウンロードして、以下のコマンドで実行できます。

```bash
$ stack build
$ stack exec genetic-algorithm-vector < input/48node.txt 
```

結果があまりよくない場合は、`poolSize`の値を増やしましょう。

なお、最適解は`33523.70850743559`だと思われます。

# 結果

`generationNum = 100`、`poolSize=6000`で試しました。

```
100 generation	117502.3637	117502.3637
99 generation	111547.6981	111547.6981
98 generation	98044.69584	98044.69584
97 generation	95280.17597	95280.17597
96 generation	95280.17597	100666.9047
95 generation	95280.17597	101406.8389
94 generation	95280.17597	100290.5332
.
.
.
7 generation	34714.01486	34714.01486
6 generation	34714.01486	34714.01486
5 generation	34714.01486	34714.01486
4 generation	34714.01486	34714.01486
3 generation	34714.01486	34714.01486
2 generation	34714.01486	34714.01486
1 generation	34714.01486	34714.01486
34714.014857028
100.689914082s
```

研究室のマシンで試したので、実行速度はかなり速いです。
グラフにしてみると、こんな感じです。

![best-result-of-GA.png](https://qiita-image-store.s3.amazonaws.com/0/136124/6052782c-825e-2bcb-3f4a-b02c5e7d854b.png)

最適解が`33523.70850743559`なので、100秒で`34714.014857028`出せるのは優秀ですね。

# 実は遺伝的アルゴリズムのライブラリがHaskellにある

[moo](https://github.com/astanin/moo)というライブラリがあるようですね。

# あとがき

がんばりました。
Haskellで遺伝的アルゴリズムを実装するのは大変かなと思いましたが、案外かけました。Haskellは型がしっかりしているし、とてもプログラムが設計しやすく、かきやすい言語だなと感じました。

Haskellでいろいろ考えながらコーディングするのが楽しかったです。おわり。

# 参考

[巡回セールスマン問題](https://www.wikiwand.com/ja/%E5%B7%A1%E5%9B%9E%E3%82%BB%E3%83%BC%E3%83%AB%E3%82%B9%E3%83%9E%E3%83%B3%E5%95%8F%E9%A1%8C)

[遺伝的アルゴリズム(wiki)](https://www.wikiwand.com/ja/%E9%81%BA%E4%BC%9D%E7%9A%84%E3%82%A2%E3%83%AB%E3%82%B4%E3%83%AA%E3%82%BA%E3%83%A0)

[遺伝的アルゴリズム](http://www.sist.ac.jp/~kanakubo/research/evolutionary_computing/genetic_algorithms.html)

[様々な選択・交叉・突然変異](http://www.sist.ac.jp/~kanakubo/research/evolutionary_computing/ga_operators.html)

[TSP Data for the Traveling Salesperson Problem](https://people.sc.fsu.edu/~jburkardt/datasets/tsp/tsp.html)