N-QueensパズルのIsingモデルをSymPyで生成し，Gurobiで解く 

# ポイント
* SymPyでスピン変数($-1/+1$)を扱う方法
* Gurobiでスピン変数を扱う方法

# N-QueensパズルのIsingモデル化
[N-QueensパズルのQUBOをSymPyで生成し，Gurobiで解く](https://qiita.com/koji_nakano/items/860bd9fb46072af986af)のIsingモデル版です．
QUBOからIsingモデルへは自動的に変換可能ですが，ここでは敢えてその変換を使用せずに，Isingモデルを直接構築します．

Isingモデルは，スピン値（$-1$および$+1$）を取る変数の二次式です．$n \times n$のスピン変数$S = (s_{i,j})$を用いて，$n \times n$のボードを表現します．ここで，$s_{i,j} = +1$は，座標$(i,j)$にクイーンが配置されていることを意味します．各行および各列に$+1$が一つだけ存在する（one-hot制約）場合，さらに，各対角線上で$+1$が0個または1個だけ存在する（zero-one-hot制約）場合，変数$S$はN-Queensパズルの解を表しています．

# Ising modelのone-hot制約とzero-one-hot制約
$n$個のスピン変数$s_0, s_1, \ldots, s_{n-1}$に対して、one-hot制約が成立するのは以下を満たした時です．
$$
s_0+s_1+\cdots +s_{n-1}=1\times (+1)+ (n-1)\times(-1)=-n+2
$$
これは，式
$$
(s_0+s_1+\cdots +s_{n-1}+n-2)^2
$$
が最小値0を取る場合に相当します．また、zero-one-hot制約を満たすのは，さらに以下を満たす時です．
$$
s_0+s_1+\cdots +s_{n-1}=0\times (+1)+ n\times(-1)=-n
$$
つまり，式
$$
(s_0+s_1+\cdots +s_{n-1}+n)(s_0+s_1+\cdots+ s_{n-1}+n-2)
$$
が最小値0をとる場合に該当します．したがって、全ての行と列がone-hot制約を満たし、さらに全ての斜めラインがzero-one-hot制約を満たす二次式を合計することで，N-Queensパズルを解くためのIsingモデルを構築できます．


# Gurobiでスピン変数を扱う方法
Gurobi Optimizerでは，スピン変数を直接扱う機能は（たぶん）提供されていません．そのため、$-1$と$+1$の値のみを取る整数変数を定義するために，2値変数を介する必要があります．具体的には，2値変数$x_{i,j}$と整数変数$s_{i,j}$の間に，次の制約を設定します。
$$
  2x_{i,j}-s_{i,j}=1
$$
この制約により，$x_{i,j}$が0の場合$s_{i,j}$は$-1$となり，$x_{i,j}$が1の場合は$s_{i,j}$が$+1$になります．

# Pythonプログラム
SymPyを利用してスピン変数を扱うために，`sp_Spin`クラスを定義しています．このクラスでは，スピン変数が指数演算により偶数乗されると1になるという性質を利用しています．`self.S`を$S$を保持するのに用います．`H`と`V`はそれぞれ`self.S`の行と列にあるスピン変数を表し，`A`と`B`は`self.S`の斜めライン上のスピン変数を表しています．`H`と`V`の各要素がone-hot制約を，AとBの各要素がzero-one-hot制約を満たす二次式の合計が`formula`となります．この`formula`を展開して得られる項のリストが`self.ising`です．

Gurobiを使用してこのIsingモデルから解を求めるために，`varX`と`varS`をGurobiの変数として用います．ここで，`varX[i,j]`は2値変数で$x_{i,j}$に対応し，`varS[self.S[i][j]]`は最小値が$-1$，最大値が$+1$の整数変数で$s_{i,j}$に対応します．最小値の指定は重要で，指定しない場合Gurobiはデフォルトで$0$が最小値となるため，明示的に$-1$を設定します。$2x_{i,j}-s_{i,j}=1$の制約を`addConstrs`メソッドを使用して追加します．最後に，`self.ising`に含まれる各一次項と二次項を，`obj`に`varS`の項として加えていきます．

```Python
#!/usr/bin/env python3

import argparse
import sympy as sp
import gurobipy as gp


class sp_Spin(sp.Symbol):
    def _eval_power(self, exp):
        if exp % 2 == 0:
            return sp.Integer(1)
        else:
            return self


class NQUEEN:
    def __init__(self, n):
        self.n = n
        self.S = [[sp_Spin(f"s_{i}_{j}") for j in range(self.n)] for i in range(self.n)]

    def generate(self):
        H = [[self.S[i][j] for j in range(self.n)] for i in range(self.n)]
        V = [[self.S[i][j] for i in range(self.n)] for j in range(self.n)]
        A = [
            [self.S[i + k][k] for k in range(self.n) if 0 <= i + k and i + k < self.n]
            for i in range(-self.n, 2 * self.n)
        ]
        B = [
            [self.S[i - k][k] for k in range(self.n) if 0 <= i - k and i - k < self.n]
            for i in range(-self.n, 2 * self.n)
        ]
        formula = (
            sum([(sum(x) + len(x) - 2) ** 2 for x in H])
            + sum([(sum(x) + len(x) - 2) ** 2 for x in V])
            + sum([(sum(x) + len(x)) * (sum(x) + len(x) - 2) for x in A if len(x) >= 2])
            + sum([(sum(x) + len(x)) * (sum(x) + len(x) - 2) for x in B if len(x) >= 2])
        )
        self.ising = formula.expand().as_coefficients_dict()

    def gurobi(self, timelimit):
        self.model = gp.Model()
        self.model.setParam("TimeLimit", timelimit)
        varX = {}
        varS = {}
        for i in range(self.n):
            for j in range(self.n):
                varX[i, j] = self.model.addVar(vtype=gp.GRB.BINARY, name=f"x_{i}_{j}")
                varS[self.S[i][j]] = self.model.addVar(
                    lb=-1.0, ub=1.0, vtype=gp.GRB.INTEGER, name=f"s_{i}_{j}"
                )
        self.optimal = -self.ising.pop(1, 0)
        self.model.addConstrs(
            (
                2 * varX[i, j] - varS[self.S[i][j]] == 1
                for i in range(self.n)
                for j in range(self.n)
            )
        )
        obj = []
        for key, val in self.ising.items():
            if isinstance(key, sp.Mul):
                quad = key.expand().args
                obj.append(val * varS[quad[0]] * varS[quad[1]])
            else:
                obj.append(val * varS[key])
        self.model.setObjective(gp.quicksum(obj), sense=gp.GRB.MINIMIZE)
        self.model.optimize()

    def print(self):
        for i in range(self.n):
            for j in range(self.n):
                print(
                    "-" if self.model.getVarByName(f"s_{i}_{j}").X < 0 else "+", end=""
                )
            print()
        print(f"optimal={self.optimal} obtained={int(self.model.ObjVal)}")


def main():
    parser = argparse.ArgumentParser(
        description="Solve N-Queens puzzle through SymPy using Gurobi"
    )
    parser.add_argument(
        "-n", default=8, type=int, help="Value of N of N-QUEENS (default:8)"
    )
    parser.add_argument(
        "-t", default=10, type=int, help="Time limit of Gurobi optimizer"
    )
    args = parser.parse_args()

    nqueen = NQUEEN(args.n)
    nqueen.generate()
    nqueen.gurobi(args.t)
    nqueen.print()


if __name__ == "__main__":
    main()

```

# 実行例
8-QueensパズルをSymPyでIsingモデル化してGurobiで解くと，値$-1264$の最適解が得られ，8-Queensパズルの正しい解になっています．
```
Gurobi Optimizer version 9.1.2 build v9.1.2rc0 (linux64)
Thread count: 16 physical cores, 32 logical processors, using up to 32 threads
Optimize a model with 64 rows, 128 columns and 128 nonzeros
Model fingerprint: 0x984ca412
Model has 728 quadratic objective terms
Variable types: 0 continuous, 128 integer (64 binary)
Coefficient statistics:
  Matrix range     [1e+00, 2e+00]
  Objective range  [4e+01, 5e+01]
  QObjective range [4e+00, 4e+00]
  Bounds range     [1e+00, 1e+00]
  RHS range        [1e+00, 1e+00]
Found heuristic solution: objective -1200.000000
Presolve removed 64 rows and 64 columns
Presolve time: 0.00s
Presolved: 728 rows, 792 columns, 2184 nonzeros
Variable types: 0 continuous, 792 integer (792 binary)

Root relaxation: objective -1.456000e+03, 64 iterations, 0.00 seconds

    Nodes    |    Current Node    |     Objective Bounds      |     Work
 Expl Unexpl |  Obj  Depth IntInf | Incumbent    BestBd   Gap | It/Node Time

     0     0 -1456.0000    0   64 -1200.0000 -1456.0000  21.3%     -    0s
H    0     0                    -1240.000000 -1456.0000  17.4%     -    0s
H    0     0                    -1248.000000 -1456.0000  16.7%     -    0s
H    0     0                    -1256.000000 -1456.0000  15.9%     -    0s
H    0     0                    -1264.000000 -1436.0000  13.6%     -    0s
     0     0 -1370.6667    0   64 -1264.0000 -1370.6667  8.44%     -    0s
     0     0 -1304.0000    0   64 -1264.0000 -1304.0000  3.16%     -    0s
     0     2 -1304.0000    0   64 -1264.0000 -1304.0000  3.16%     -    0s

Cutting planes:
  Gomory: 1
  MIR: 354
  Zero half: 149
  RLT: 289
  BQP: 227

Explored 111 nodes (2489 simplex iterations) in 0.69 seconds
Thread count was 32 (of 32 available processors)

Solution count 5: -1264 -1256 -1248 ... -1200
No other solutions better than -1264

Optimal solution found (tolerance 1.00e-04)
Best objective -1.264000000000e+03, best bound -1.264000000000e+03, gap 0.0000%
----+---
-+------
---+----
------+-
--+-----
-------+
-----+--
+-------
optimal=-1264 obtained=-1264

```