取得データ量を最大化するセンサー設置箇所を求める

# <i class='fa fa-building' /> やりたいこと
下図のようなエリアがあります。

![image](https://qiita-image-store.s3.amazonaws.com/0/13955/f5e1be52-1f8a-7c5f-12ae-ec5f408dfe74.png)

- この三角形の中央にセンサーを設置して、データ収集をします。
- **総取得データ量を最大化**することが目的です。
- 1つの三角形に1つまでしかセンサーは置けません。
- 1つのセンサーを置くと、データ量 **3** を取得できます。
- 予め、センサーを設置しなければいけない箇所と、センサーを設置してはいけない箇所が決まっています。
- 下表のように、隣りに別のセンサーがあると、干渉を起こして、取得データ量が各々 **1** 減ります。(合計 **2** 減少)
    - 従って、周りに2以上のセンサーがある場合は、設置しない方がよいことになります。

|周りのセンサー数|取得データ量|干渉データ量|全体として増えるデータ量|
|:--|:--|:--|:--|
|0|3|0|3|
|1|3|2|1|
|2|3|4|-1|
|3|3|6|-3|

# <i class='fa fa-building' /> 混合整数最適化問題による最適化(数理問題によるアプローチ)
混合整数最適化問題(Mixed Integer optimization Problem: MIP)としてモデル化してみましょう。
定式化や考え方については、[組合せ最適化を使おう](http://qiita.com/SaitoTsutomu/items/bfbf4c185ed7004b5721)を参考にしてください。

## 定式化
<table>
<tr><td>最大化</td><td>$取得データ量 \\ = 3 \times 設置するかどうか - 2 \times 干渉するかどうか$</td></tr>
<tr><td rowspan="2">変数</td><td>$(センサーを)設置するかどうか \in\{0，1\} \forall 設置箇所$</td></tr>
<tr><td>$干渉するかどうか \ge 0 \forall 隣接箇所$</td></tr>
<tr><td>制約条件</td><td>$設置するかどうか_1 + 設置するかどうか_2 - 干渉するかどうか \le 1$</td></tr>
</table>

# <i class='fa fa-building' /> 最小カット問題による最適化(典型問題によるアプローチ)
グラフ理論で[最小カット問題](https://ja.wikipedia.org/wiki/%E3%82%AB%E3%83%83%E3%83%88_(%E3%82%B0%E3%83%A9%E3%83%95%E7%90%86%E8%AB%96)#.E6.9C.80.E5.B0.8F.E3.82.AB.E3.83.83.E3.83.88)というものがあります。

以下のようなs-tグラフを考えます。この有向グラフのs-t最小カットを求めると、元のセンサー設置問題の最適解を得ることができます。
元の問題は最大化問題ですが、"3 － 取得データ量"をコストと考えることにより、最小化問題としてとらえます。
詳しくは、後述の参考サイトをご覧ください。

## s-tグラフ
![image](https://qiita-image-store.s3.amazonaws.com/0/13955/a083a4e6-4429-1763-6592-3dfb5932c1c5.png)

- (0,0),(0,1),(1,0)の3つの例を上図に示します。
- 始端ノード"s"から、各箇所(三角形)に辺を引きます。また、各箇所から終端ノード"t"にも辺を引きます。
- 上三角の箇所(上図では(0,0)と(1,0))では、sからの辺は"設置"を表し、tへの辺は"非設置"を表します。
- 下三角の箇所(上図では(0,1))では、sからの辺は"非設置"を表し、tへの辺は"設置"を表します。
- 設置の辺の重みは0とし、非設置の辺の重みは3とします。(重みは、"3 － 取得データ量"です。)
- 全ての隣接箇所に対して、下三角から上三角に、重み2の辺を引きます。
- "設置"で固定する場合、"設置"の辺は削除(カット済み)とし、"非設置"の辺の重みを十分大きくします。
- "非設置"で固定する場合、"非設置"の辺は削除(カット済み)とし、"設置"の辺の重みを十分大きくします。

## s-tグラフの意味

- s-tカットするためには、各々の箇所において、設置または非設置の辺のどちらかをカットすることになります。カットされた方を選んだものとみなします。
- s-tカットが成り立つためには、隣接箇所共に設置した場合、干渉の辺もカットする必要があり、干渉量を表現できます。(上三角と下三角で分けたことで上手くいきます。)

## 最小カット問題について
最小カット問題は、多項式時間で解ける解きやすい問題です。後述のプログラムでは、Pythonのnetworkxを用いて解きます。
ちなみに、最小カット問題の[双対問題](http://qiita.com/SaitoTsutomu/items/d1812ff9b5ccf0ecc716)は、[最大流問題](http://qiita.com/SaitoTsutomu/items/80e70da6717acacefa00)です。

# <i class='fa fa-building' /> pythonによる実行例
準備をします。

```py3:python3
import numpy as np, networkx as nx
from pulp import *
def addvar(lowBound=0, var_count=[0], *args, **kwargs):
    """変数作成用ユーティリティ"""
    var_count[0] += 1
    return LpVariable('v%d' % var_count[0], lowBound=lowBound, *args, **kwargs)
def calc(a, r):
    """rを設置箇所として取得データ量を計算"""
    b = a.copy()
    b[b > 1] = 0
    for x, y in r:
        b[y, x] = 1
    s = b.sum() * 3
    for y in range(0, b.shape[0], 2):
        for x in range(b.shape[1]):
            s -= 2 * b[y, x] * b[y+1,x]
            if x:
                s -= 2 * b[y, x] * b[y+1,x-1]
            if y:
                s -= 2 * b[y, x] * b[y-1,x]
    return s
```

solve_by_mipは、MIPによる解法です。設置箇所を返します。

```py3:python3
def solve_by_mip(a):
    """MIPで問題を解き、設置箇所を返す"""
    nm, nn = a.shape
    b = a.astype(object)
    vv1 = [addvar(cat=LpBinary) for _ in range((b > 1).sum())]
    b[b > 1] = vv1
    vv2 = []
    m = LpProblem(sense=LpMaximize)
    for y in range(0, nm, 2):
        for x in range(nn):
            chk(m, vv2, b[y,x] + b[y+1,x])
            if x: chk(m, vv2, b[y,x] + b[y+1,x-1])
            if y: chk(m, vv2, b[y,x] + b[y-1,x])
    m += 3 * lpSum(vv1) - 2 * lpSum(vv2)
    m.solve()
    return [(x, y) for x in range(nn) for y in range(nm)
            if isinstance(b[y,x], LpVariable) and value(b[y, x]) > 0.5]
def chk(m, vv2, e):
    """eが変数を含むならば、共に1であれば目的関数を2減らす制約を追加"""
    if isinstance(e, LpAffineExpression):
        v = addvar()
        vv2.append(v)
        m += e - v <= 1
```

solve_by_graphは、最小カットによる解法です。同じく設置箇所を返します。


```py3:python3
def solve_by_graph(a):
    """最小カット問題で問題を解き、設置箇所を返す"""
    nm, nn = a.shape
    g = nx.DiGraph()
    for y in range(0, nm, 2):
        for x in range(nn):
            if a[y, x] == 0: # off
                g.add_edge('s', (x,y), capacity=7)
            elif a[y, x] == 1: # on
                g.add_edge((x,y), 't', capacity=7)
            else:
                g.add_edge('s', (x,y), capacity=0)
                g.add_edge((x,y), 't', capacity=3)
            if a[y+1, x] == 0: # off
                g.add_edge((x,y+1), 't', capacity=7)
            elif a[y+1, x] == 1: # on
                g.add_edge('s', (x,y+1), capacity=7)
            else:
                g.add_edge('s', (x,y+1), capacity=3)
                g.add_edge((x,y+1), 't', capacity=0)
            g.add_edge((x,y+1), (x,y), capacity=2)
            if x:
                g.add_edge((x-1,y+1), (x,y), capacity=2)
            if y:
                g.add_edge((x,y-1), (x,y), capacity=2)
    r = []
    for s in nx.minimum_cut(g, 's', 't')[1]:
        b = 's' in s
        for t in s:
            if isinstance(t, str): continue
            x, y = t
            if a[y, x] > 1 and b == (y%2 != 0):
                r.append((x, y))
    return sorted(r)
```

40×80のサイズで、固定箇所をランダムに設置して、結果を比較してみましょう。

```py3:python3
nn, nm = 40, 80 # 横、縦
np.random.seed(1)
a = np.random.randint(0, 6, (nm, nn)) # 0; fix off, 1: fix on, ow:select

%time rmip = calc(a, solve_by_mip(a))
%time rgrp = calc(a, solve_by_graph(a))
print(rmip == rgrp)
>>>
Wall time: 455 ms
Wall time: 185 ms

True
```

- どちらの手法も同じ取得データ量になっている(rmip == rgrp)のが確認できます。
- MIPの方が2倍強、遅いです。
    - 一般に、汎用ソルバーより専用ソルバーの方が計算が速いです。
    - 別途確認したところ、MIPソルバー単体の計算時間は、最小カットによる計算時間を少し上回るくらいでした。
- また、理屈はわかっていないですが、MIPを線形緩和しても取得データ量は変わらず、計算時間は若干早くなります。

参考サイト

- [Problem E. The Year of Code Jam](https://code.google.com/codejam/contest/32011/dashboard#s=p4): 参考にした問題
- [最小カットを使って「燃やす埋める問題」を解く](http://www.slideshare.net/shindannin/project-selection-problem): 参考にした解法

以上