Numpy Container and N-queens Problem

In the previous container section, we have introduced the tupledict container to store and manipulate multi-dimensional data.

However, due to the Bring Your Own Container (BYOC) principle, variables and constraints in PyOptInterface can just simple Python objects that can be stored in Numpy ndarrays directly as a multi-dimensional array, and you can enjoy the features of Numpy such like fancy-indexing automatically.

We will use N-queens problem as example to show how to use Numpy ndarrays as container to store 2-dimensional variables and construct optimization model.

Firstly, we import the necessary modules:

import numpy as np
import pyoptinterface as poi
from pyoptinterface import highs

model = highs.Model()

Then we create a 2-dimensional variable x with shape \(N \times N\) to represent the placement of queens. Each element of x is a binary variable that indicates whether a queen is placed at the corresponding position. We use object as the data type of x to store the binary variables. The following code snippet creates the variables:

N = 8

x = np.empty((N, N), dtype=object)
for i in range(N):
    for j in range(N):
        x[i, j] = model.add_variable(domain=poi.VariableDomain.Binary)

Running HiGHS 1.8.1 (git hash: 4a7f24a): Copyright (c) 2024 HiGHS under MIT licence terms

Next, we add the constraints to ensure that each row, each column has exact one queen, and each diagonal has at most one queen.

for i in range(N):
    # Row and column
    model.add_linear_constraint(poi.quicksum(x[i, :]), poi.Eq, 1.0)
    model.add_linear_constraint(poi.quicksum(x[:, i]), poi.Eq, 1.0)
for i in range(-N+1, N):
    # Diagonal
    model.add_linear_constraint(poi.quicksum(x.diagonal(i)), poi.Leq, 1.0)
    # Anti-diagonal
    model.add_linear_constraint(poi.quicksum(np.fliplr(x).diagonal(i)), poi.Leq, 1.0)

Finally, we solve the model.

model.optimize()

print("Termination status:", model.get_model_attribute(poi.ModelAttribute.TerminationStatus))

Coefficient ranges:
  Matrix [1e+00, 1e+00]
  Cost   [0e+00, 0e+00]
  Bound  [1e+00, 1e+00]
  RHS    [1e+00, 1e+00]
Presolving model
42 rows, 64 cols, 252 nonzeros  0s
42 rows, 64 cols, 270 nonzeros  0s
Objective function is integral with scale 1

Solving MIP model with:
   42 rows
   64 cols (64 binary, 0 integer, 0 implied int., 0 continuous)
   270 nonzeros
MIP-Timing:      0.0012 - starting analytic centre calculation

Src: B => Branching; C => Central rounding; F => Feasibility pump; H => Heuristic; L => Sub-MIP;
     P => Empty MIP; R => Randomized rounding; S => Solve LP; T => Evaluate node; U => Unbounded;
     z => Trivial zero; l => Trivial lower; u => Trivial upper; p => Trivial point

        Nodes      |    B&B Tree     |            Objective Bounds              |  Dynamic Constraints |       Work      
Src  Proc. InQueue |  Leaves   Expl. | BestBound       BestSol              Gap |   Cuts   InLp Confl. | LpIters     Time

         0       0         0   0.00%   0               inf                  inf        0      0      0         0     0.0s
         0       0         0   0.00%   0               inf                  inf        0      0      6        29     0.0s
 L       0       0         0 100.00%   0               0                  0.00%     1559     28     84       213     0.1s
         1       0         1 100.00%   0               0                  0.00%     1559     28     84       267     0.1s

Solving report
  Status            Optimal
  Primal bound      0
  Dual bound        0
  Gap               0% (tolerance: 0.01%)
  P-D integral      0
  Solution status   feasible
                    0 (objective)
                    0 (bound viol.)
                    0 (int. viol.)
                    0 (row viol.)
  Timing            0.08 (total)
                    0.00 (presolve)
                    0.00 (solve)
                    0.00 (postsolve)
  Max sub-MIP depth 1
  Nodes             1
  Repair LPs        0 (0 feasible; 0 iterations)
  LP iterations     267 (total)
                    0 (strong br.)
                    184 (separation)
                    54 (heuristics)
Termination status: TerminationStatusCode.OPTIMAL

The solution can be obtained and visualized by the following code:

get_v = np.vectorize(lambda x: model.get_value(x))
x_value = get_v(x)

print(x_value.astype(int))

[[0 0 1 0 0 0 0 0]
 [1 0 0 0 0 0 0 0]
 [0 0 0 0 0 0 1 0]
 [0 0 0 0 1 0 0 0]
 [0 0 0 0 0 0 0 1]
 [0 1 0 0 0 0 0 0]
 [0 0 0 1 0 0 0 0]
 [0 0 0 0 0 1 0 0]]