はじめに

本稿では、PythonとGoogleのortoolsを用いて最適化を学んでいきます。
私のようなプログラマーが理解しやすいよう、用語も意図的にそちらへ寄せています。
プログラムコードもなるべくそちらへ寄せています。
第二回の今回はイラストロジック（白黒）をCP-SATで解きます。

準備

本題ではないためGoogleのColaboratoryを使用することで、
ローカル環境への構築などはすっ飛ばすことにします。

記事投稿段階では
Python3.7
ortools9.4
を使用しています。

ortoolsインストール

Colab上で以下を実行します。

!pip install --upgrade --user ortools

※tensorflow2.8の依存関係でエラーが出ますが、本稿では問題とならないため無視します。

画面の指示通りRuntimeをRestartしておきましょう。

Colabでは、Runtimeの再起動が起こるとインストールパッケージは
失われます。このpipもColab上に記載しておいて必要に応じて
動かすと後々楽になります。

確認

今回もCP-SATソルバーを利用します。
importで確認。

from ortools.sat.python import cp_model

実行してエラーが出なければひとまずOK.

イラストロジック（イラロジ）

イラロジ（白黒）の基本ルールは以下になります。

ルール① 二次元のマス上に黒か白が入ります。
ルール② 行方向、列方向それぞれに複数の数字が記載されています
ルール③ 数字の数分、連続して黒マスが続き、1つ以上の白マスを挟んで次の数字分の黒マスとつながります。

問題

今回は以下の問題を使用しています。

※自作画像から問題を起こしており、一意に回答が決まらない可能性があります。

元画像はこちらです。

問題配列はこちら。

row_problem = [[3], [4], [5], [7], [4, 3], [3, 3], [3, 3], [3, 3], [3, 4], [4, 3], [3, 3, 2], [3, 3, 4], [4, 3, 2, 3, 5], [3, 5, 5, 3, 6], [20, 7], [21, 8], [23, 9], [24, 9], [25, 9], [26, 9], [27, 10], [29, 9], [30, 9], [31, 9], [32, 9], [33, 9], [35, 9], [36, 9], [37, 9], [37, 9]]
col_problem = [[3], [4], [6], [8], [9], [11], [13], [14], [16], [18], [19], [21], [5, 16], [5, 17], [6, 18], [5, 18], [5, 18], [5, 17], [4, 16], [5, 17], [5, 18], [6, 18], [5, 17], [5, 17], [5, 16], [20], [18], [17], [15], [14], [12], [2, 10], [6, 9], [10, 7], [12, 5], [14, 4], [14, 2], [15], [16], [14], [13], [11], [10], [8], [6], [5], [3]]

シンプルな制約を与えてみる

まずは答えを格納する領域として、2次元のBool変数を用意します。

# モデルインスタンス生成と答え変数の作成

# 必要ライブラリインポート
from ortools.sat.python import cp_model
import collections

# 確認のために使用する構造体の宣言
item_type_r = collections.namedtuple('item_type', 'row value start end size interval')
item_type_c = collections.namedtuple('item_type', 'col value start end size interval')

# modelインスタンス生成
model = cp_model.CpModel();

# row_problemの要素数が領域の行数
# col_problemの要素数が領域の列数
r_cnt = len(row_problem)
c_cnt = len(col_problem)

# 答えを格納する変動変数領域を作成する
answer = {}
for i in range(r_cnt):
  for j in range(c_cnt):
    answer[(i, j)] = model.NewBoolVar("ans_%i_%i" %(i, j))

制約（ルールの適用）

まずは一番単純な、以下の３制約を適用してみます。

①全マスの黒（＝True）の総数は、行または列の合算値に等しい
　※行の総和と列の総和は等しい

black=\sum_{i,k} (r_{i,k}) \\
\sum_{i,k} (r_{i,k}) = \sum_{j,k} (c_{j,k}) \\

②各行の黒の総数は、各行の問題の総和に等しい

black_i = \sum_{k}(r_{i,k})

③各列の黒の総和は、各列の問題の総和に等しい

black_j = \sum_{k}(c_{j,k})

# 基本制約の設定

# 答え変数に対して以下の制約を適用する
# 本来、行、列のループは一度で済むが、説明の便宜上分けたまま実装する。
# ①全黒マス（True）数は、行側または列側の総和に等しい
sum_black = sum(sum(row) for row in row_problem)
# 答え変数内でTrueになる総和は、行の総和（＝列の総和）に等しい
model.Add( sum_black == sum(answer[(i, j)] for i in range(r_cnt) for j in range(c_cnt)) )

# ②各行の総和は、問題に記載された各行の総和に等しい
for i in range(r_cnt):
  sum_value = sum(row_problem[i])
  # 答え変数の該当行に対し、行総和を制約しておく
  model.Add( sum_value == sum(answer[(i, j)] for j in range(c_cnt)) )

# ③各列の総和は、問題に記載された各列の総和に等しい
for j in range(c_cnt):
  sum_value = sum(col_problem[j])
  # 答え変数の該当列に対し、列総和を制約しておく
  model.Add( sum_value == sum(answer[(i, j)] for i in range(r_cnt)) )

これだけで解けるならば話が早いので試しに実行してみます。

実行

# モデルからsolverを作成し実行する

# solver生成
solver = cp_model.CpSolver()
status = solver.Solve(model)

確認

どんな結果が出たか確認します。

# 実行結果の確認を行う
ans = [[False] * len(col_problem) for j in range(len(row_problem))]

if ((status == cp_model.OPTIMAL) or (status == cp_model.FEASIBLE)):
  print("OPTIMAL or FEASIBLE")
  print(solver.ResponseStats())

  # 答え変数の内容を写す
  for i in range(r_cnt):
    for j in range(c_cnt):
      ans[i][j] = solver.Value(answer[(i, j)])
else:
    print("could not optimize")

ResponseStats

OPTIMAL or FEASIBLE
CpSolverResponse summary:
status: OPTIMAL
objective: 0
best_bound: 0
booleans: 0
conflicts: 0
branches: 0
propagations: 0
integer_propagations: 0
restarts: 0
lp_iterations: 0
walltime: 0.122243
usertime: 0.122244
deterministic_time: 0.00445035
gap_integral: 0

得られた結果を視覚的に確認します。

# 答えを表示する
# 表示用のoutput_answer関数は本稿末尾に掲載しています
HTML(output_answer(ans, row_problem, col_problem))

面影はありますが、やはりこの制約だけでは無理でした。

連続する黒を担保する

行方向、列方向で、問題に即して連続する黒を担保する必要があります。
本稿では、CP-SATのoverlap制約を使用して実現します。

overlap制約は、範囲を意味するinterval変数のリストを受け取り、
各interval変数の開始、終了値が重ならないよう制約してくれる制約になります。
今回は各行、各列毎にoverlap制約を用います。

参考）https://developers.google.com/optimization/reference/python/sat/python/cp_model#addnooverlap

行方向のoverlap制約

# 各行/列毎にoverlapを用いて黒マスと白マスが重ならないよう配置するため、
# 各要素に該当するIntervalVarを生成する。
interval_var_row_lists = []

# 答え変数へ制約を転写するため、行単位でリスト化するための変数
cond_for_rows = {}

# TODO:本来不要、実装中の値確認用に変数を格納しておく
all_val_row = {}
all_i = 0

for i in range(r_cnt):
  ivr_buf = []
  sum_value = sum(row_problem[i])

  # 現在行の要素数
  len_value = len(row_problem[i])

  # 一行分の条件を格納する変数
  cond_for_row = {}

  # 和が0より大きければ条件がある。
  if sum_value > 0:
    for k in range(len_value):

      # 各要素間の白マスは最低1マス
      # Intervalとoverlapの都合上、長さを-1するか、start <= x < end として
      # 扱う必要がある。本実装では後者を採用するため、start, endの取りえる範囲は
      # 問題サイズ+1となる。⇒ 0 ～ 問題サイズを取り合える範囲とする
      s_var = model.NewIntVar(0, c_cnt, "rsv_b_%i_%i"%(i,k))

      # 二つ目以降の要素の場合、一つ前の終了位置から1マス以上開けた位置が開始位置になる。
      if k > 0:
        model.Add( s_var >= e_var + 1)

      e_var = model.NewIntVar(0, c_cnt, "rev_b_%i_%i"%(i,k))

      # Interval変数にしリストに格納する
      # 黒要素のサイズは各要素の値となる
      buf = model.NewIntervalVar(s_var, row_problem[i][k], e_var, "rrv_b_%i_%i"%(i,k))
      ivr_buf.append(buf)

      # 答え変数への転写用に黒マスの要素を行毎にリスト化しておく
      item = item_type_r(row=i, value=row_problem[i][k], start=s_var, end=e_var, size=row_problem[i][k], interval=buf)
      cond_for_row[(k)] = item

      # TODO: 本来不要（確認用）
      all_val_row[all_i] = item
      all_i += 1

    # 1行分の情報が出来たのでこれをoverlap制約に加える
    model.AddNoOverlap(ivr_buf)
    interval_var_row_lists.append(ivr_buf)

    # 1行分の条件を格納しておく
    cond_for_rows[i] = cond_for_row

答え変数への転写

上で作成したoverlap制約を、答え変数として用意した二次元変数へ転写します。

# 行側のInterval変数を答え変数へ転写するための制約を設定する。
# 同様に列側も同じ答え変数へ転写することで、行側、列側が同時に満たされる制約となる。

# 行毎の各要素に応じた中間バッファ領域用の変数
buf_middle_r_cond = {}

for i in range(r_cnt):
  for j in range(c_cnt):
    # 該当マスがTrueになりえるのは、要素内の何れか一つと、
    # start <= answer[i,j] < end
    # が成り立つ場合のみ。
    # And,Orの組み合わせ制約を直接入れるのは難しいが、
    # それぞれなら適用できるので、要素の数と等しいbool面を用意し、
    # 各要素毎に転写する。
    # 例）3 4 □□□□□□□□□□
    #  に対し、3⇒start=1, end=3
    #          4⇒start=5, end=8
    # 転写後は 3に対する面⇒■■■□□□□□□□
    #          4に対する面⇒□□□□■■■■□□
    # これを合算し答え変数に転写する二段階の処理とする。
    for k in range(len(cond_for_rows[(i)])):
      # 中間バッファ領域を用意し各要素の該当バッファに対して制約する
      buf_middle_r_cond[(i,j,k)] = model.NewBoolVar("buf_middle_r_cond%i_%i_%i" %(i, j, k))
      model.Add( buf_middle_r_cond[(i,j,k)] * (c_cnt - j - 1) <= c_cnt - cond_for_rows[(i)][(k)].start - 1)
      model.Add( buf_middle_r_cond[(i,j,k)] * j < cond_for_rows[(i)][(k)].end)

  # 各行で総和を制約しておくことで余計なTrueや、Falseを発生させないようにする
  for k in range(len(cond_for_rows[(i)])):
    model.Add( cond_for_rows[(i)][(k)].size == sum(buf_middle_r_cond[(i,j,k)] for j in range(c_cnt)) )

# 各要素でバッファリングしたもののAndに答え変数を制約する
for i in range(r_cnt):
  for j in range(c_cnt):
    model.Add( answer[(i,j)] == sum( buf_middle_r_cond[(i,j,k)] for k in range(len(cond_for_rows[(i)]))) )

列方向のoverlap制約

同様に列方向からも制約します。
本稿では敢えて別コードにしていますが、実際は共通関数化が可能です。

# 列側も同様に作成する。
# 本来は関数化し共通処理とするか、または1ループで完結するようにした方が
# 効率的だが、説明の便宜上敢えて別で作成することとする。

# 各行/列毎にoverlapを用いて黒マスと白マスが重ならないよう配置するため、
# 各要素をIntervalVarを用いて生成する。
interval_var_col_lists = []

# 答え変数へ制約を転写するため、列単位でリスト化するための変数
cond_for_cols = {}

for j in range(c_cnt):
  ivc_buf = []
  sum_value = sum(col_problem[j])

  # 現在列の要素数
  len_value = len(col_problem[j])

  # 一列分の条件を格納する変数
  cond_for_col = {}

  if sum_value > 0:
    for k in range(len_value):

      # 要素間の白マスは最低1マス
      s_var = model.NewIntVar(0, c_cnt, "csv_b_%i_%i"%(j,k))

      # 二つ目以降の要素の場合、一つ前の白マスのend以降でしか始められないよう制約する
      if k > 0:
        model.Add( s_var >= e_var + 1)

      e_var = model.NewIntVar(0, c_cnt, "cev_b_%i_%i"%(j,k))

      # 範囲変数にしリストに格納する
      # 要素のサイズは問題定義値そのもので固定
      buf = model.NewIntervalVar(s_var, col_problem[j][k], e_var, "crv_b_%i_%i"%(j,k))
      ivc_buf.append(buf)

      # 答え変数への転写用に黒マスの要素を列毎にリスト化しておく
      item = item_type_c(col=j, value=col_problem[j][k], start=s_var, end=e_var, size=col_problem[j][k], interval=buf)
      cond_for_col[(k)] = item

    # 1列分の情報が出来たのでこれをoverlap制約に加える
    model.AddNoOverlap(ivc_buf)
    interval_var_col_lists.append(ivc_buf)

    # 1列分の条件を格納しておく
    cond_for_cols[j] = cond_for_col

# 列毎の要素に応じた中間バッファ領域用の変数
buf_middle_c_cond = {}

for j in range(c_cnt):
  for i in range(r_cnt):
    for k in range(len(cond_for_cols[(j)])):
      # 中間バッファ領域を用意し各要素の該当バッファ対して制約する
      buf_middle_c_cond[(i,j,k)] = model.NewBoolVar("buf_middle_c_cond%i_%i_%i" %(i, j, k))
      model.Add( buf_middle_c_cond[(i,j,k)] * (r_cnt - i - 1) <= r_cnt - cond_for_cols[(j)][(k)].start - 1)
      model.Add( buf_middle_c_cond[(i,j,k)] * i < cond_for_cols[(j)][(k)].end)

  # 各列にて、総和を制約しておくことで余計なTrueや、Falseを発生させないようにする
  for k in range(len(cond_for_cols[(j)])):
    model.Add( cond_for_cols[(j)][(k)].size == sum(buf_middle_c_cond[(i,j,k)] for i in range(r_cnt)) )

# 各要素でバッファリングしたもののAndに答え変数を制約する
for i in range(r_cnt):
  for j in range(c_cnt):
    model.Add( answer[(i,j)] == sum( buf_middle_c_cond[(i,j,k)] for k in range(len(cond_for_cols[(j)]))) )

確認

実行して結果を確認します。

無事に解くことができました。

処理時間（walltime）をはじめとするRessponseStatsは以下でした。

OPTIMAL or FEASIBLE
CpSolverResponse summary:
status: OPTIMAL
objective: 0
best_bound: 0
booleans: 0
conflicts: 0
branches: 0
propagations: 0
integer_propagations: 0
restarts: 0
lp_iterations: 0
walltime: 0.162361
usertime: 0.162361
deterministic_time: 0.00504747
gap_integral: 0

解説

overlap制約を用いる場合、まずinterval変数を用意します。
今回のやり方では、問題に記載されている要素の数分、
intervalVarを用意し、当該要素の開始位置と終了位置を算出させています。

例えば10マスの盤目に3 5の行問題が定義されているとします。

この場合以下の３パターンが考えられます。

3、5それぞれのInterval変数を用意します。
Interval変数は開始、終了、サイズを保有しています。
このうち、開始と終了を更に整数変数で定義し、サイズは固定値としています。

概念的には以下のようになります。

こうして定義したInterval変数を、各行、各列ごとにリストへ入れ、
overlap制約へ与えることで、各行、各列のInterval変数が重なり合わないよう制約されます。

では実際のコードで解説します。

s_var = model.NewIntVar(0, c_cnt, "rsv_b_%i_%i"%(i,k))

開始位置を格納する整数変数を用意しています。
kによるループの最初の一回目は、各行の一つ目の問題要素になるため、
開始位置は0以降となります。

二つ目以降の要素では、一つ前の終了位置+1（白マス）以上から開始されます。
kのループが2回目以降は、e_var変数に一回前の終了位置を意味する整数変数が格納されているため、

# 二つ目以降の要素の場合、一つ前の終了位置から1マス以上開けた位置が開始位置になる。
if k > 0:
  model.Add( s_var >= e_var + 1)

として、一つ前の終了位置+1以上　という制約を付与しています。

e_var = model.NewIntVar(0, c_cnt, "rev_b_%i_%i"%(i,k))

終了位置変数e_varを上書きにて作成しなおします。

# Interval変数にしリストに格納する
# 黒要素のサイズは各要素の値となる
buf = model.NewIntervalVar(s_var, row_problem[i][k], e_var, "rrv_b_%i_%i"%(i,k))
ivr_buf.append(buf)

IntervalVarを生成しリストへ格納しています。
IntervalVarでは、開始位置、終了位置と、サイズの定義ができ、
複数のIntervalVarを同一リストに格納した後に、overlap制約に与えることで、
各IntervalVarが重なり合わないよう制約されます。

# 1行分の情報が出来たのでこれをoverlap制約に加える
model.AddNoOverlap(ivr_buf)

overlap制約として追加しました。

行、列それぞれでこの制約を起こすことで、各要素が取りえる位置が決まります。

ただし、このままでは行列各々が収まる位置を求めるだけで相互に作用しません。
このため、最初に用意した二次元の答え変数を使用し、行側、列側の制約を
それぞれ転写することで、同時に満たす答えを求めるようにしてみました。

ここがシンプルにはできず、ずいぶんと考えた結果、変数を増やしてしまいました。
後述するパフォーマンスが出ない要因の一端がここにあります。

for i in range(r_cnt):
  for j in range(c_cnt):
    # 該当マスがTrueになりえるのは、要素内の何れか一つと、
    # start <= answer[i,j] < end
    # が成り立つ場合のみ。
    # And,Orの組み合わせ制約を直接入れるのは難しいが、
    # 要素それぞれなら適用できるので、要素の数と等しいbool面を用意し、
    # 各要素毎に転写する。
    # 例）3 4 □□□□□□□□□□
    #  に対し、3⇒start=1, end=3
    #          4⇒start=5, end=8
    # 転写後は 3に対する面⇒■■■□□□□□□□
    #          4に対する面⇒□□□□■■■■□□
    # これを合算し答え変数に転写する二段階の処理とする。
    for k in range(len(cond_for_rows[(i)])):
      # 中間バッファ領域を用意し各要素の該当バッファに対して制約する
      buf_middle_r_cond[(i,j,k)] = model.NewBoolVar("buf_middle_r_cond%i_%i_%i" %(i, j, k))
      model.Add( buf_middle_r_cond[(i,j,k)] * (c_cnt - j - 1) <= c_cnt - cond_for_rows[(i)][(k)].start - 1)
      model.Add( buf_middle_r_cond[(i,j,k)] * j < cond_for_rows[(i)][(k)].end)

  # 当該の行面に対し、総和を制約しておくことで余計なTrueや、Falseを発生させないようにする
  for k in range(len(cond_for_rows[(i)])):
    model.Add( cond_for_rows[(i)][(k)].size == sum(buf_middle_r_cond[(i,j,k)] for j in range(c_cnt)) )

一度中間のバッファ領域に要素数と等しい次元数の変数を用意しています。
Interval変数の各値を、一次元のバッファー領域(buf_middle_r_cond*)へ受け、
バッファー領域の合計値を答え変数へ制約する二段階の方法で実現しています。

図中の緑枠に該当する処理が以下になります。

# 各要素でバッファリングしたもののAndに答え変数を制約する
for i in range(r_cnt):
  for j in range(c_cnt):
    model.Add( answer[(i,j)] == sum( buf_middle_r_cond[(i,j,k)] for k in range(len(cond_for_rows[(i)]))) )

列方向も同様に制約することで、行側、列側が同時に満たされるケースを探査させることが出来ます。

行側、列側が制約を同時に満たしているということは、必然的に、
それぞれの行、それぞれの列、および盤面すべての総和が満たされます。
始めに作った各行列の総和制約はいらないこととなりますが、
速度に対して貢献する場合もあるため、ひとまず残しておき、
他の問題を解いてみます。

元画像

問題配列

row_problem = [[3], [5], [5], [4], [3, 1, 3], [5, 5], [5, 5], [4, 7, 5], [2, 11, 3], [13], [14], [15], [16], [4, 17, 4], [4, 17, 4], [5, 17, 5], [4, 17, 4], [2, 16, 2], [15], [15], [13], [3, 11, 3], [5, 9, 5], [5, 3, 5], [4, 5], [2, 1, 2], [4], [5], [5], [3]]
col_problem = [[4], [5], [5], [4], [1], [3, 3], [5, 5], [5, 6, 5], [4, 10, 4], [2, 12, 2], [14], [15], [16], [3, 16, 3], [4, 17, 4], [5, 17, 5], [4, 17, 4], [2, 16, 2], [16], [15], [14], [12], [3, 9, 3], [5, 4, 5], [5, 5], [5, 4], [3, 1, 3], [4], [5], [5], [4]]

結果

OPTIMAL or FEASIBLE
CpSolverResponse summary:
status: OPTIMAL
objective: 0
best_bound: 0
booleans: 0
conflicts: 0
branches: 0
propagations: 0
integer_propagations: 0
restarts: 0
lp_iterations: 0
walltime: 0.238011
usertime: 0.238011
deterministic_time: 0.0197176
gap_integral: 0

無事解けているようです。
ここで、問題規模を大きくしてみます。

元画像

少しわかりづらいですが我が家で昔飼っていた猫の写真から作成しています。
※猫だと思ってみていただければ、かろうじて寝ている猫に見えるでしょうか・・・。

問題配列

row_problem = [[5, 16, 7, 8], [5, 24, 6, 8], [1, 2, 8, 11, 10, 2, 7, 1], [1, 1, 7, 4, 10, 2, 8, 1], [4, 5, 11, 3, 7, 1], [2, 1, 9, 4, 3, 9, 3, 6, 2], [1, 1, 3, 8, 1, 2, 11, 10, 2], [1, 1, 7, 17, 12], [1, 7, 20, 16], [3, 8, 27, 7, 5], [2, 9, 27, 7, 8], [41, 12, 3], [1, 40, 3, 4], [38, 2, 1, 5], [32, 6, 7], [2, 32, 1, 6, 8], [35, 8, 1, 10], [61], [62], [36, 26], [35, 1, 27], [37, 28], [1, 38, 28], [35, 1, 1, 28], [36, 1, 28], [4, 28, 1, 1, 28], [19, 13, 2, 27], [1, 1, 8, 5, 14, 7, 19, 3], [1, 2, 9, 16, 4, 23], [31, 5, 21], [15, 12, 1, 4, 16], [15, 12, 1, 11], [29, 2, 9], [32, 11], [35, 1, 9], [30, 5, 1, 1, 2], [33, 5, 1, 6], [36, 2, 1, 4], [19, 18, 1, 5], [17, 1, 17, 3], [1, 41, 2], [1, 31, 1, 1, 5], [1, 30, 9], [2, 20, 6, 2, 12], [2, 21, 5, 12], [3, 18, 2, 3, 13], [20, 1, 1, 2, 6], [20, 1, 10], [20, 1, 1, 2, 6, 3], [31, 1, 2, 2], [32, 4], [2, 29, 7], [2, 17, 6, 10], [7, 10, 1, 10], [3, 1, 8, 4, 9], [1, 1, 1, 11, 9, 1], [4, 7, 1, 9, 1], [2, 1, 3, 11], [1, 2, 1, 1, 1, 1, 9], [1, 1, 6, 1, 8], [1, 1, 1, 6, 9], [1, 1, 4, 6, 3], [2, 4], [1, 3, 5], [2, 4, 4], [2, 5, 3], [6, 2], [5, 2], [4, 2], [3, 2], [1, 5], [1, 3], [4], [1, 5], [1, 6], [1, 7], [1, 6], [7], [8, 1], [7, 1, 12]]
col_problem = [[8, 1, 13, 1], [1, 2, 2, 11, 1], [1, 1, 1, 6, 2, 1], [7, 6, 3, 1], [3, 10, 1, 2, 2, 1], [2, 13, 2, 2, 1], [4, 15, 3, 1, 1], [1, 4, 20], [6, 17, 1, 1], [1, 5, 1, 18, 2, 1, 1], [7, 1, 1, 18, 1, 1, 1, 1], [6, 1, 20, 2, 3], [4, 1, 3, 21, 1, 5], [3, 2, 4, 5, 22, 1, 6], [2, 9, 7, 22, 1, 3, 6], [2, 9, 11, 24, 1, 4, 7], [3, 8, 10, 1, 24, 1, 3, 7], [3, 9, 9, 2, 26, 3, 3, 3], [3, 9, 12, 30, 1, 3, 3, 3], [3, 48, 3, 2, 4, 2, 2], [3, 23, 17, 4, 3, 3, 7, 1], [3, 18, 17, 4, 4, 1, 4, 7, 2], [3, 1, 17, 16, 3, 2, 3, 4, 5, 1], [3, 28, 5, 4, 3, 2, 4, 1, 1], [3, 28, 7, 4, 3, 8], [2, 28, 3, 2, 4, 4, 3], [3, 27, 5, 4, 3], [1, 1, 16, 11, 4, 4, 2], [3, 34, 4, 2, 1], [3, 35, 4, 1], [2, 37, 3], [45, 2], [27, 14, 3], [5, 19, 2, 14], [1, 3, 38], [1, 40, 1], [1, 39, 4], [41, 4], [2, 39, 3, 3], [2, 38, 4, 2, 2], [2, 46, 6], [1, 31, 4, 7, 6], [3, 31, 4, 8, 8], [1, 2, 31, 12, 10], [2, 2, 31, 10, 10], [3, 1, 31, 9, 11], [6, 16, 1, 4, 4, 6, 14], [4, 21, 3, 2, 5, 14], [7, 15, 3, 4, 3, 12], [5, 14, 3, 2, 5], [7, 6, 5, 4, 7], [13, 4, 1, 1, 1, 1], [10, 9, 4, 3], [10, 15, 2], [7, 15, 1], [5, 15], [3, 1, 12, 1], [1, 2, 2, 16, 1], [2, 24, 1, 1, 1], [26, 2, 2], [7, 18, 2, 3], [4, 21, 2], [4, 2, 18, 2], [2, 1, 21], [4, 18, 1], [1, 2, 20], [4, 20], [3, 19], [3, 20], [2, 20], [23], [21], [20], [19], [17], [16], [9, 2], [11], [9], [8]]

結果

一応解けましたが、実行時間が長すぎました。

OPTIMAL or FEASIBLE
CpSolverResponse summary:
status: OPTIMAL
objective: 0
best_bound: 0
booleans: 63537
conflicts: 8904
branches: 87236
propagations: 13329436
integer_propagations: 5088507
restarts: 56474
lp_iterations: 7313050
walltime: 2353.19
usertime: 2353.19
deterministic_time: 930.615
gap_integral: 0

40分弱・・・。
さきほど、始めにつけた総和の制約はいらないと記載しましたが、
まずはこれを取り払って確認してみます。

実行・・・

## 総和制約除去後
OPTIMAL or FEASIBLE
CpSolverResponse summary:
status: OPTIMAL
objective: 0
best_bound: 0
booleans: 63540
conflicts: 5036
branches: 78228
propagations: 11508787
integer_propagations: 4427437
restarts: 58098
lp_iterations: 251699
walltime: 241.675
usertime: 241.675
deterministic_time: 23.491
gap_integral: 0

実行時間が1/10になりました。
総和制約が足を引っ張っていたようです。

探査範囲を絞り込める制約は、無駄が無ければいい方向に働きますが、
無駄な制約の場合lp_iterationsが増えすぎて悪化方向に働くようです。

※lp_iterationsは動かす度に多少の差が現れます。
※この規模になると同じ処理でも実行時間に差が出る場合があります。

また、今回はなるべく読みやすいプログラムを掲載するため、
最適化処理としては回りくどいことをしています。
答え変数への転写のためにバッファー変数を増やしているなど、
改善の余地が多くあります。

今後さらなる改善が出来た折にはまたご紹介しようと思います。

おまけ

確認関数
視覚的確認のための描画関数。

from IPython.display import HTML

# 描画用関数
def output_answer(a_ans, p_row, p_col):
  table_style = "<style>\
                 <!--\
                 table.rp{\
                   border:#1175ac solid 3px;\
                   font-size:10px;\
                   border-collapse: collapse;\
                   table-layout: fixed;\
                   width: " + str(13 * len(p_col) + 121) + "px;\
                 }\
                 table.rp tr{\
                   height: 12px;\
                 }\
                 table.rp td{\
                   border:#1175ac solid 1px;\
                   width: 12px;\
                   padding: 0;\
                   text-align:center !important;\
                   background-color:white;\
                 }\
                 table.rp td.bl{\
                   background-color:#ddd;\
                 }\
                 table.rp th{\
                   border: #1175ac solid 1px;\
                   padding: 0;\
                   background-color:white;\
                   vertical-align: bottom;\
                   width: 12px;\
                 }\
                 table.rp th:nth-of-type(1) {\
                   text-align:right !important;\
                   vertical-align: middle;\
                   white-space: nowrap;\
                   width: 120px;\
                 }\
                 -->\
                 </style>"
  # 1行目は列問題、1,1は空白
  ret = "<tr><th></th>"
  for j in range(len(p_col)):
    buf = " ".join([str(_) for _ in p_col[j]])
    ret += "<th>" + buf + "</th>"
  ret += "</tr>"
  for i in range(len(p_row)):
    # 先頭に行問題付与
    buf = " ".join([str(_) for _ in p_row[i]])
    ret += "<tr><th>" + buf + "</th>"
    for j in range(len(a_ans[i])):
      bgcolor = ""
      if (a_ans[i][j]):
         bgcolor = ' class="bl"'
      ret += "<td"+bgcolor+">" + "</td>"
    ret += "</tr>"

  return table_style + "<table class=\"rp\">" + ret + "</table>"