81、数独求解（策略求解、整数规划建模求解）

81、数独求解（策略求解、整数规划建模求解）

import re
from pyscipopt import Model, quicksum

#set存value的可能性，当value值确认后，set清空
class Cell:
    def __init__(self):
        self.value = "_"
        self.set = {1,2,3,4,5,6,7,8,9}

    def __repr__(self):
        return str(self.value)

#读数据
def dataRead():
    dataPath = r"/Users/zhangchaoyu/PycharmProjects/pythonProject/Game/data/sudoku1.csv"
    data = []
    f = open(dataPath, encoding='UTF-8-sig')
    while True:
        row = f.readline().strip('n')
        if row == ',n' or row == '':
            break
        row = re.split(',', row)
        data.append(row)
    return data

def case():
    data = [["", "", "", "", 2, "", "", 8, 5],
            [6, "", 2, "", 7, "", "", "", ""],
            [3, "", "", "", "", "", "", 1, ""],
            [7, "", "", 9, "", 4, "", "", ""],
            ["", 4, "", "", "", "", "", 3, ""],
            ["", "", 9, "", "", 5, "", 2, ""],
            ["", "", 5, "", "", "", 2, 6, ""],
            [2, "", "", "", 8, 6, 1, "", ""],
            ["", 3, "", "", "", 2, "", 7, ""]]
    return data

#把该行的num可能清理
def XRowRefresh(X, idx, num):
    for j in range(9):
        if num in X[idx][j].set:
            X[idx][j].set.remove(num)

#把该列的num可能清理
def XColRefresh(X, idx, num):
    for i in range(9):
        if num in X[i][idx].set:
            X[i][idx].set.remove(num)

#把该block的num可能清理
def XBlockRefresh(X, row, col, num):
    x, y = row//3, col//3
    for i in range(3*x, 3*(x+1)):
        for j in range(3*y, 3*(y+1)):
            if num in X[i][j].set:
                X[i][j].set.remove(num)

#X初始化
def XInit(data):
    X = [[Cell() for _ in range(9)] for _ in range(9)]
    for i in range(9):
        for j in range(9):
            if data[i][j] != "":
                num = int(data[i][j])
                X[i][j].value = num
                X[i][j].set = {}
                XRowRefresh(X, i, num)
                XColRefresh(X, j, num)
                XBlockRefresh(X, i, j, num)
    return X

#打印X
def XPrint(X):
    for i in range(9):
        print(X[i])
    print("")

#贪心策略，如果一个格子只有一种可能性，则直接使用
def Greedy(X):
    flag = False
    for i in range(9):
        for j in range(9):
            if len(X[i][j].set) == 1:
                flag = True
                num = X[i][j].set.pop()
                print("i:", i, ", j:", j, ", num:", num)
                X[i][j].value = num
                XRowRefresh(X, i, num)
                XColRefresh(X, j, num)
                XBlockRefresh(X, i, j, num)
    return flag

def mustExistRow(X, num):
    flag = False
    for i in range(9):
        Xi = [X[i][j].value for j in range(9)]
        #如果该行存在则过
        if num in Xi:
            continue
        #如果不存在，则看能放下该值的位置是否唯一，唯一则更新
        idxs = []
        for j in range(9):
            if num in X[i][j].set:
                idxs.append(j)
        if len(idxs) == 1:
            flag = True
            j = idxs[0]
            print("i:", i, ", j:", j, ", num:", num)
            X[i][j].value = num
            X[i][j].set = {}
            XRowRefresh(X, i, num)
            XColRefresh(X, j, num)
            XBlockRefresh(X, i, j, num)
    return flag

def mustExistRCol(X, num):
    flag = False
    for j in range(9):
        Xj = [X[i][j].value for i in range(9)]
        # 如果该列存在则过
        if num in Xj:
            continue
        # 如果不存在，则看能放下该值的位置是否唯一，唯一则更新
        idxs = []
        for i in range(9):
            if num in X[i][j].set:
                idxs.append(i)
        if len(idxs) == 1:
            flag = True
            i = idxs[0]
            print("i:", i, ", j:", j, ", num:", num)
            X[i][j].value = num
            X[i][j].set = {}
            XRowRefresh(X, i, num)
            XColRefresh(X, j, num)
            XBlockRefresh(X, i, j, num)
    return flag

def mustExistBlock(X, num):
    flag = False
    for x in range(3):
        for y in range(3):
            Xxy = [X[i][j].value for i in range(3*x, 3*(x+1)) for j in range(3*y, 3*(y+1))]
            # 如果该列存在则过
            if num in Xxy:
                continue
            # 如果不存在，则看能放下该值的位置是否唯一，唯一则更新
            idxs = []
            for i in range(3*x, 3*(x+1)):
                for j in range(3*y, 3*(y+1)):
                    if num in X[i][j].set:
                        idxs.append((i,j))
            if len(idxs) == 1:
                flag = True
                i,j = idxs[0]
                print("i:", i, ", j:", j, ", num:", num)
                X[i][j].value = num
                X[i][j].set = {}
                XRowRefresh(X, i, num)
                XColRefresh(X, j, num)
                XBlockRefresh(X, i, j, num)
    return flag

#反向找，一行，一列，一个九宫格内只有一个格能放下某个数，那么没错，就是它了
def mustExist(X):
    flag = False
    for num in range(1, 10):
        flagRow = mustExistRow(X, num)
        flagCol = mustExistRCol(X, num)
        flagBlock = mustExistBlock(X, num)
        flag = flag|flagRow|flagCol|flagBlock
    return flag

#各种搜索
def XRefresh(X):
    XPrint(X)
    while True:
        flagGreedy = Greedy(X)
        flagmustExist = mustExist(X)
        if flagGreedy|flagmustExist:
            XPrint(X)
        else:
            break

def subIP(X, timeLimit=100, isPrint=True):
    # 建模求解
    model = Model("IP")
    Y = {(i,j):[model.addVar(vtype="B", name="Y[%s, %s, %s]"%(i, j, k)) for k in range(9)] for i in range(9) for j in range(9)}

    #随便设置个目标
    model.setObjective(Y[0,0][0], "minimize")

    #每个数每行，每列只能出现一次
    for k in range(9):
        for i in range(9):
            model.addCons(quicksum(Y[i, j][k] for j in range(9)) == 1)
        for j in range(9):
            model.addCons(quicksum(Y[i, j][k] for i in range(9)) == 1)
        for x in range(3):
            for y in range(3):
                model.addCons(quicksum(Y[i, j][k] for i in range(3*x, 3*(x+1)) for j in range(3*y, 3*(y+1))) == 1)

    #每个地方只能出现一个数
    for i in range(9):
        for j in range(9):
            model.addCons(quicksum(Y[i, j][k] for k in range(9)) == 1)

    #现有的数要满足
    for i in range(9):
        for j in range(9):
            if X[i][j].value != "_":
                k = X[i][j].value-1
                model.addCons(Y[i, j][k] == 1)
                for l in range(9):
                    if l == k:
                        continue
                    model.addCons(Y[i, j][l] == 0)
            else:
                for k in range(9):
                    if k+1 not in X[i][j].set:
                        model.addCons(Y[i, j][k] == 0)

    # 设置求解时间
    model.setRealParam("limits/time", timeLimit)
    model.hideOutput()
    model.optimize()

    if isPrint:
        print("ngap:", model.getGap())

    Y1 = {(i,j):[round(model.getVal(Y[i,j][k])) for k in range(9)] for i in range(9) for j in range(9)}

    return Y1

def IP(X):
    #建模求解
    Y = subIP(X)
    #刷新X
    for i in range(9):
        for j in range(9):
            if X[i][j].value != "_":
                continue
            for k in range(9):
                if Y[i, j][k] == 1:
                    X[i][j].value = k+1
                    break
    XPrint(X)

if __name__ == '__main__':
    #数据录入
    # data = dataRead()
    data = case()
    #X初始化
    X = XInit(data)
    #X迭代
    XRefresh(X)
    #IP求解
    IP(X)