Advanced constraint programming modeling

advancedModelingWeek5.py

from ortools.sat.python import cp_model

''' course: https://www.coursera.org/learn/advanced-modeling#about 
These models have different classes of symmetries. 
Week6 is about symmetries.
Currently ortools don't have methods for symmetry breaking '''



def crossbowTrapProblem(): # nqueens problem
    model = cp_model.CpModel()

    sizeOfField = 7

    trapVarDict = {}
    for i in range(sizeOfField):
        for j in range(sizeOfField):
            trapVarDict[(i, j)] = model.NewBoolVar("trapIsOn"+str(i)+str(j))

    model.Add(sum([trapVarDict[(i,j)] for i in range(sizeOfField) for j in range(sizeOfField)]) == sizeOfField)

    for i in range(sizeOfField):
        varList = []
        for j in range(sizeOfField):
            varList.append(trapVarDict[(i,j)])
        model.Add(sum(varList) <= 1)

    for j in range(sizeOfField):
        varList = []
        for i in range(sizeOfField):
            varList.append(trapVarDict[(i,j)])
        model.Add(sum(varList) <= 1)


    for k in range(1-sizeOfField, sizeOfField):  # forward diagonals
        varList = []
        for i in range(sizeOfField):
            for j in range(sizeOfField):
                if k == i-j:
                    varList.append(trapVarDict[(i,j)])
        model.Add(sum(varList) <= 1)


    for k in range((2*sizeOfField)-1):  # backward diagonals
        varList = []
        for i in range(sizeOfField):
            for j in range(sizeOfField):
                if k == i+j:
                    varList.append(trapVarDict[(i,j)])
        model.Add(sum(varList) <= 1)


    solver = cp_model.CpSolver()
    status = solver.Solve(model)

    if status == cp_model.INFEASIBLE:
        print("Infeasible model")

    if status == cp_model.OPTIMAL:
        for i,j in trapVarDict:
            if solver.Value(trapVarDict[i,j]) > 0.0:
                print("trap is placed at "+str(i)+", "+str(j))
    print("----------------------------")



#----------------------------------------------------------------------------------------------------------------------#




def theLampModel():
    model = cp_model.CpModel()

    numOfRows = 7
    numOfColumns = 56
    numOfLampsOnEachRow = 24
    numOfLampsOnEachColumn = 3
    numOfCommonLitLampsInAnyTwoRows = 8

    lampsLitVarDict = {}
    for row in range(numOfRows):
        for col in range(numOfColumns):
            lampsLitVarDict[(row, col)] = model.NewBoolVar("lampLitAt"+str(row)+str(col))

    for row in range(numOfRows):
        varList = []
        for col in range(numOfColumns):
            varList.append(lampsLitVarDict[(row, col)])
        model.Add(sum(varList) == numOfLampsOnEachRow)

    for col in range(numOfColumns):
        varList = []
        for row in range(numOfRows):
            varList.append(lampsLitVarDict[(row, col)])
        model.Add(sum(varList) == numOfLampsOnEachColumn)

    for row1 in range(numOfRows-1):
        for row2 in range(row1+1, numOfRows):
            varList = []
            for col in range(numOfColumns):
                isCommonLitLamp = model.NewBoolVar("isCommonLitLampAt"+str(col)+"for"+str(row1)+str(row2))
                model.AddBoolAnd([lampsLitVarDict[(row1, col)], lampsLitVarDict[(row2, col)]]).OnlyEnforceIf(isCommonLitLamp)
                varList.append(isCommonLitLamp)
            model.Add(sum(varList) == numOfCommonLitLampsInAnyTwoRows)

    solver = cp_model.CpSolver()
    status = solver.Solve(model)

    if status == cp_model.INFEASIBLE:
        print("Infeasible model")

    if status == cp_model.OPTIMAL:
        for row in range(numOfRows):
            for col in range(numOfColumns):
                if solver.Value(lampsLitVarDict[(row, col)]) > 0.0:
                    print("Lamp lit at row_"+str(row)+", col_"+str(col))
    print("----------------------------")



#----------------------------------------------------------------------------------------------------------------------#



def windPrayingModel():
    model = cp_model.CpModel()

    numOfVertices = 6
    numOfPolygon = 3
    tributeList = ["LETTUCE", "TURNIP", "CARROT", "BANANA", "PEACH", "GRAPE", "SNAKE", "CHICKEN", "FISHHEAD", "RICEWINE",
                   "GAOLIANG", "GRAPEWINE", "GOLD", "SILVER", "BRONZE", "GREEN", "RED", "BLUE"]
    rankList = [1, 1, 1, 2, 2, 2, 3, 3, 3, 4, 4, 4, 5, 5, 5, 6, 6, 6]
    rankTributesDict = {1: [0,1,2], 2: [3,4,5], 3:[6,7,8], 4:[9,10,11], 5:[12,13,14], 6:[15,16,17]}
    numOfTributes = len(tributeList)
    numOfRanks = max(rankList)


    tributeVarList = []
    polygonVertexIndexDict, counter = {}, 0
    for polygon in range(numOfPolygon):
        for vertex in range(numOfVertices):
            tributeVarList.append(model.NewIntVar(0, numOfTributes-1, "tributeFor"+str(polygon)+str(vertex)))
            polygonVertexIndexDict[(polygon, vertex)] = counter
            counter += 1

    numCoinsVarList = []
    for polygon in range(numOfPolygon-1):
        for vertex in range(numOfVertices):
            numCoinsVarList.append(model.NewIntVar(0, numOfRanks, "numCoinsfor"+str(polygon)+str(vertex)))

    model.AddAllDifferent(tributeVarList)

    rankVarList = []
    for polygon in range(numOfPolygon):
        for vertex in range(numOfVertices):
            rankVarList.append(model.NewIntVar(1, numOfRanks, "rankFor"+str(polygon)+str(vertex)))
            index = polygonVertexIndexDict[(polygon, vertex)]
            model.AddElement(tributeVarList[index], rankList, rankVarList[index])

    positionVarList = []
    for tribute in range(numOfTributes):
        positionVarList.append(model.NewIntVar(0, (numOfPolygon*numOfVertices)-1, "positionForTribute"+str(tribute)))

    model.AddInverse(tributeVarList, positionVarList)

    for polygon in range(numOfPolygon-1):
        for vertex in range(numOfVertices):
            diff = model.NewIntVar(-numOfRanks, numOfRanks, "diffInRanksFor"+str(polygon)+str(polygon+1)+str(vertex))
            index = polygonVertexIndexDict[(polygon, vertex)]
            index2 = polygonVertexIndexDict[(polygon+1, vertex)]
            model.Add(diff == rankVarList[index] - rankVarList[index2])
            model.AddAbsEquality(numCoinsVarList[index], diff)

    for vertex in range(numOfVertices-1):  # removing variable symmetry
        index = polygonVertexIndexDict[(0, vertex)]
        index2 = polygonVertexIndexDict[(0, vertex+1)]
        model.Add(tributeVarList[index] < tributeVarList[index2])

    index = polygonVertexIndexDict[(0, 0)]
    index2 = polygonVertexIndexDict[(numOfPolygon-1, 0)]
    model.Add(tributeVarList[index] < tributeVarList[index2]) # removing variable symmetry

    for rank, tributeList in rankTributesDict.items(): # removing variable symmetry/value symmetry
        for tribute in range(len(tributeList)-1):
            model.Add(positionVarList[tribute] < positionVarList[tribute+1])

    model.Maximize(sum(numCoinsVarList))

    solver = cp_model.CpSolver()
    status = solver.Solve(model)

    if status == cp_model.INFEASIBLE:
        print("Infeasible model")

    if status == cp_model.OPTIMAL:
        print('Max number of coins: %i' % solver.ObjectiveValue())
    print("----------------------------")




#----------------------------------------------------------------------------------------------------------------------#




def generalAllocationProblem ():
    model = cp_model.CpModel()

    numOfGenerals = 14
    numOfRoads = 14
    coopList = [[0, 15, 3, 13, 4, 1, 5, 4, 4, 1, 1, 11, 10, 1],
                [15, 0, 3, 13, 4, 11, 10, 9, 11, 7, 1, 5, 4, 1],
                [3, 3, 0, 13, 4, 11, 10, 9, 10, 7, 1, 5, 4, 1],
                [13, 13, 13, 0, 4, 11, 10, 9, 9, 7, 1, 5, 4, 1],
                [4, 4, 4, 4, 0, 1, 5, 4, 10, 1, 1, 11, 10, 1],
                [1, 11, 11, 11, 1, 0, 10, 9, 8, 7, 1, 5, 4, 1],
                [5, 10, 10, 10, 5, 10, 0, 9, 8, 7, 1, 5, 4, 1],
                [4, 9, 9, 9, 4, 9, 9, 0, 8, 7, 1, 5, 4, 1],
                [4, 11, 10, 9, 10, 7, 8, 8, 0, 7, 1, 5, 4, 1],
                [1, 7, 7, 7, 1, 7, 7, 7, 7, 0, 1, 5, 4, 1],
                [1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 5, 4, 1],
                [11, 5, 5, 5, 11, 5, 5, 5, 5, 5, 5, 0, 4, 4],
                [10, 4, 4, 4, 10, 4, 4, 4, 4, 4, 4, 4, 0, 3],
                [1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 4, 3, 0]]
    coopListFlattened = sum(coopList, [])
    maxCoop = max(coopListFlattened)

    roadVarList = []
    for road in range(numOfRoads):
        roadVarList.append(model.NewIntVar(0, numOfGenerals-1, "generalForRoad"+str(road)))

    model.AddAllDifferent(roadVarList)

    coopVarList = []
    for road in range(numOfRoads-1):
        indexVar = model.NewIntVar(0, (numOfGenerals*numOfGenerals)-1, "indexFor"+str(road)+str(road+1))
        coopVarList.append(model.NewIntVar(0, maxCoop, "coopFor"+str(road)+str(road+1)))
        model.Add(indexVar == (roadVarList[road] * numOfGenerals) + roadVarList[road+1])
        model.AddElement(indexVar, coopListFlattened, coopVarList[road])


    model.Add(roadVarList[0] < roadVarList[numOfRoads-1]) # variable symmetry breaking

    for road1 in range(1, numOfRoads-2): # dominance breaking constraint
        for road2 in range(road1+1, numOfRoads-1):

            coop1Var = model.NewIntVar(0, maxCoop, "coop1For"+str(road1-1)+str(road2))
            coop2Var = model.NewIntVar(0, maxCoop, "coop2For"+str(road1)+str(road2+1))
            index1Var = model.NewIntVar(0, (numOfGenerals*numOfGenerals)-1, "index1For"+str(road1-1)+str(road2))
            index2Var = model.NewIntVar(0, (numOfGenerals*numOfGenerals)-1, "index2For" +str(road1) + str(road2+1))
            model.Add(index1Var == (roadVarList[road1-1] * numOfGenerals) + roadVarList[road2])
            model.Add(index2Var == (roadVarList[road1] * numOfGenerals) + roadVarList[road2+1])
            model.AddElement(index1Var, coopListFlattened, coop1Var)
            model.AddElement(index2Var, coopListFlattened, coop2Var)

            model.Add(coopVarList[road1-1] + coopVarList[road2] >= coop1Var + coop2Var)


    # solList = [10, 12, 4, 8, 9, 7, 6, 5, 2, 3, 1, 0, 11, 13]
    # for i,var in enumerate(roadVarList):
    #     model.Add(var == solList[i])

    model.Maximize(sum(coopVarList))

    solver = cp_model.CpSolver()
    status = solver.Solve(model)

    if status == cp_model.INFEASIBLE:
        print("Infeasible model")

    if status == cp_model.OPTIMAL:
        print('Max coop: %i' % solver.ObjectiveValue())
    print("----------------------------")




#----------------------------------------------------------------------------------------------------------------------#



if __name__ == '__main__':


    crossbowTrapProblem()
    theLampModel()
    windPrayingModel()
    generalAllocationProblem()