<jupyter_start><jupyter_text>The Knapsack Problem by <NAME> is licensed under a Creative Commons Attribution 4.0 International License.Based on a work at https://github.com/mathcoding/opt4ds. **NOTE:** Run the following script whenever running this script on a Google Colab.<jupyter_code>import shutil
import sys
import os.path

if not shutil.which("pyomo"):
    !pip install -q pyomo
    assert(shutil.which("pyomo"))

if not (shutil.which("glpk") or os.path.isfile("glpk")):
    if "google.colab" in sys.modules:
        !apt-get install -y -qq glpk-utils
    else:
        try:
            !conda install -c conda-forge glpk 
        except:
            pass<jupyter_output><empty_output><jupyter_text>$n$-Queens ProblemThe $n$-Queens puzzle is the problem of placing eight chess queens on an $n \times n$ chessboard so that no two queens threaten each other; thus, a solution requires that no two queens share the same row, column, or diagonal (source: [wikipedia](https://en.wikipedia.org/wiki/Eight_queens_puzzle)). A solution exists for all natural numbers n with the exception of $n = 2$ and $n = 3$.  **Example:** For $n=8$, we have the following solution:```1 . . . . . Q . .2 . . . Q . . . . 3 . . . . . . Q .4 Q . . . . . . .5 . . . . . . . Q6 . Q . . . . . .7 . . . . Q . . .8 . . Q . . . . .  a b c d e f g h ```  Integer Linear Programming ModelThe $n$-Queens problem can be formalized with the following **ILP** model.**Data:** Size of the board $n\times n$. Let $I=:\{1,\dots,n\}$ a set of indices.**Decision Variables:** The variable $x_{ij} \in \{0,1\}$ is equal to 1 if we place a queen in position $(i,j)$  on the chessboard.**Objective function:** Since the problem is a feasibility problem, we can set the objective function equal to any constant value.**Constraints:** We need the following linear constraints, which encode the puzzle rules:1. Each queen appears once per row:$$    \sum_{j \in I} x_{ij} = 1, \forall i \in I$$2. Each queen appears once per column:$$    \sum_{i \in I} x_{ij} = 1, \forall j \in I$$3. Each queen appears once per main diagonals:$$    \sum_{(i,j) \in D_k} x_{ij} \leq 1, D_k \mbox{ main diagonals}$$4. Each queen appears once per off-diagonals:$$    \sum_{(i,j) \in O_k} x_{ij} \leq 1, O_k \mbox{ off diagonals}$$  Main Diagonals $D_k$Since we need to specify the pairs of indices that define as a function of $n$, we first defined the following nested loop:<jupyter_code>n = 5
for j in range(-n+2,n-1):
    for i in range(1, n+1):
        if 0 < j+i <= n:
            print(i, j+i, end='\t')
        else:
            print('   ', end='\t')
    print()<jupyter_output><empty_output><jupyter_text>Off Diagonals $_k$Similarly, we can define the off diagonals as follows:<jupyter_code>for i in reversed(range(-n+3, n)):
    for j in range(1, n):
        if 0 < n - j+i <= n:
            print(j, n-j+i, end='\t')
        else:
            print('   ', end='\t')
    print()<jupyter_output><empty_output><jupyter_text>Full Model defined in PyomoIf we put all the definitions together, we can solve the $n$-Queens problem with the script below.Please, note the following Pyomo syntax used to define variable $x_{ij}$ over the [RangeSet](https://pyomo.readthedocs.io/en/stable/library_reference/aml/index.htmlpyomo.environ.RangeSet) $I$ and $J$:```model.I = RangeSet(1, n)model.J = RangeSet(1, n)model.x = Var(model.I, model.J, within=Binary)```Notice also the syntax used to define the row and column constraints, which uses `lambda` function to define constraint rules:```model.row = Constraint(model.I,                        rule = lambda mod, i: sum(mod.x[i,j] for j in mod.J) == 1)```Finally, to define the main and of diagonals, we use the [ConstraintList](https://pyomo.readthedocs.io/en/stable/working_models.html) class:```model.mainD = ConstraintList()...   model.mainD.add( expr <= 1 )```The complete Pyomo script is as follows.<jupyter_code># Import the libraries
from pyomo.environ import ConcreteModel, Var, Objective, Constraint, SolverFactory
from pyomo.environ import maximize, Binary, RangeSet, ConstraintList


n = 8

# Create concrete model
model = ConcreteModel()

model.I = RangeSet(1, n)
model.J = RangeSet(1, n)

# Variables
model.x = Var(model.I, model.J, within=Binary)

# Objective Function: Maximize Profit
model.obj = Objective(expr = n, sense = maximize)

# 1. Row constraints
def VincoloRighe(mod, i):
    return sum(mod.x[i,j] for j in mod.J) == 1

model.row = Constraint(model.I, 
                       rule = VincoloRighe)

# 2. Column constraints
model.column = Constraint(model.J, 
                          rule = lambda mod, j: sum(mod.x[i,j] for i in mod.I) == 1)

# 3. Main Diagonal constraints
model.mainD = ConstraintList()
# Build the list of possible pairs
for j in range(-n+2,n-1):
    expr = 0
    for i in model.I:
        if 0 < j+i <= n:
            expr += model.x[i, j+i] 
    model.mainD.add( expr <= 1 )

# 4. Off Diagonal constraints
model.offD = ConstraintList()
# Build the list of possible pairs
for i in range(-n+3,n+1):
    expr = 0
    for j in model.J:
        if 0 < n-j+i <= n:
            expr += model.x[j, n-j+i] 
    model.offD.add( expr <= 1 )<jupyter_output><empty_output><jupyter_text>To solve the script, we use a solver factory, specifying the GLPK solver, and we inspect the Solver **status** (infeasible, unbounded, or optimal).<jupyter_code># Solve the model
sol = SolverFactory('glpk').solve(model)

# Basic info about the solution process
for info in sol['Solver']:
    print(info)<jupyter_output><empty_output><jupyter_text>We aspect the optimal decision variables (only the positive variables).<jupyter_code># Report solution value
print("Optimal solution value: z =", model.obj())
print("Decision variables:")
for i in model.I:
    for j in model.J:
        if model.x[i,j]() > 0:
            print("x({},{}) = {}".format(i, j, model.x[i,j]()))<jupyter_output><empty_output><jupyter_text>And finally, we print a solution on a simplified chessboard $n\times n$.<jupyter_code>print('\nChessboard Solution:')
for i in model.I:
    for j in model.J:
        if model.x[i,j]() > 0:
            print('Q', end=' ')
        else:
            print('.', end=' ')
    print()<jupyter_output><empty_output><jupyter_text>Plotting a solution with a Chessboard<jupyter_code># CREDIT: Solution original appeared on Stackoverflow at:
# https://stackoverflow.com/questions/60608055/insert-queen-on-a-chessboard-with-pyplot

def PlotSolution(n, x, size=6):
    import matplotlib.pyplot as plt
    import numpy as np

    chessboard = np.zeros((n, n))

    chessboard[1::2,0::2] = 1
    chessboard[0::2,1::2] = 1

    plt.figure(figsize=(size, size))
    plt.imshow(chessboard, cmap='binary')

    for i, j in x:
        if x[i,j]() > 0:
            plt.text(i-1, j-1, '♕', color='darkorange',
                     fontsize=56*size/n, fontweight='bold', ha='center', va='center')
    plt.xticks([])
    plt.yticks([])
    plt.show()
PlotSolution(n, model.x)<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.14083702328238823, 'fraction_numerical': 0.011998285959148694, 'mean_word_length': 3.251366120218579, 'pattern_counts': {'":': 0, '<': 45, '<?xml version=': 0, '>': 37, 'https://': 5, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 3, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '8', 'dir': 'jupyter-structured-clean-dedup', 'id': '3148', 'license': 'MIT', 'n_tokens_mistral': 2450, 'n_tokens_neox': 2347, 'n_words': 771, 'path': 'aa2021/notebooks/N-Queens.ipynb', 'repo_name': 'gianlucacovini/opt4ds', 'size': '7011'}