pulp-lparray

PuLP x NumPy, a match made in heaven.

An lparray is a NumPy array of PuLP variables, affine expressions or constraints. Broadcasting, indexing, sums and matrix products build whole families of expressions at once; <=, >= and == build families of constraints; and a handful of methods linearise abs, max, min, argmax, top-k, logic and piecewise-linear functions. Never see a for or an index variable again.

It needs Python 3.12+ and PuLP 4. PuLP 4 ships no solver: pip install "pulp[cbc]" adds CBC.

Super Sudoku

A Super Sudoku is a Sudoku in which the cells at the same position within their boxes also hold distinct digits.

import numpy as np
from pulp import LpBinary, LpProblem

from pulp_lparray import lparray

clues = np.array(
    [
        [0, 4, 0, 0, 1, 8, 0, 0, 0],
        [0, 5, 0, 0, 0, 2, 0, 0, 0],
        [2, 0, 8, 0, 0, 0, 7, 1, 4],
        [1, 8, 2, 0, 0, 0, 5, 6, 0],
        [0, 6, 0, 0, 0, 0, 0, 0, 0],
        [4, 7, 0, 0, 0, 6, 0, 0, 2],
        [0, 0, 5, 0, 9, 7, 0, 0, 1],
        [0, 0, 4, 0, 0, 0, 0, 0, 0],
        [9, 0, 0, 0, 0, 0, 3, 0, 0],
    ]
)  # 0 marks an empty cell

prob = LpProblem("super_sudoku")
# board[R, C, r, c, n] == 1: the cell in row 3R + r and column 3C + c holds digit n + 1
board = lparray.create(prob, "board", (3, 3, 3, 3, 9), cat=LpBinary)
(board.sum(axis=4) == 1).constrain(prob, "cell")
(board.sum(axis=(1, 3)) == 1).constrain(prob, "row")
(board.sum(axis=(0, 2)) == 1).constrain(prob, "column")
(board.sum(axis=(2, 3)) == 1).constrain(prob, "box")
(board.sum(axis=(0, 1)) == 1).constrain(prob, "position")
cells = board.transpose(0, 2, 1, 3, 4).reshape(9, 9, 9)  # row, column, digit
row, col = clues.nonzero()
(cells[row, col, clues[row, col] - 1] == 1).constrain(prob, "clue")
prob.solve()

grid = cells.values.argmax(axis=2) + 1
assert (grid[row, col] == clues[row, col]).all()
assert all(sorted(line) == list(range(1, 10)) for line in [*grid, *grid.T])

Portfolio allocation

Whole numbers of shares whose values best match target weights within a budget, measured by relative error so that every position counts alike. Minimising deviation subject to deviation >= error and deviation >= -error makes it |error| without binaries; abs is for magnitudes that the objective does not push down.

from pulp import LpInteger, LpMinimize

price = np.array([31.0, 77.0, 12.5, 220.0])
weight = np.array([0.4, 0.3, 0.2, 0.1])
funds = 10_000

prob = LpProblem("allocation", LpMinimize)
shares = lparray.create(prob, "shares", 4, lowBound=0, cat=LpInteger)
(shares @ price <= funds).constrain(prob, "budget")
target = funds * weight / price
error = (shares - target) / target
deviation = lparray.create(prob, "deviation", 4)
(deviation >= error).constrain(prob, "over")
(deviation >= -error).constrain(prob, "under")
prob += deviation.sum().item()
prob.solve()

Facility logistics

Open exactly two facilities per region, ship every product to every customer, and pay an operating cost with a volume discount on each facility's throughput. The discount makes the cost concave, so lp_interp selects its segment with binaries.

region = np.array([0, 0, 0, 1, 1, 1])
capacity = np.array([12, 6, 6, 12, 6, 6])
open_cost = np.array([10, 5, 7, 10, 5, 7])
ship_cost = np.array(
    [
        [2, 2, 3, 9, 9],
        [3, 1, 2, 9, 9],
        [2, 3, 1, 9, 9],
        [9, 9, 9, 2, 2],
        [9, 9, 9, 3, 1],
        [9, 9, 9, 2, 3],
    ]
)  # facility × customer
demand = np.array([[4, 2, 1], [0, 3, 2], [1, 1, 0], [3, 0, 1], [2, 2, 2]])  # customer × product

prob = LpProblem("logistics", LpMinimize)
is_open = lparray.create(prob, "open", 6, cat=LpBinary)
flow = lparray.create(prob, "flow", (6, 5, 3), lowBound=0)
(flow.sum(axis=0) == demand).constrain(prob, "demand")
throughput = flow.sum(axis=(1, 2))
(throughput <= capacity * is_open).constrain(prob, "capacity")
((region == np.arange(2)[:, None]) @ is_open == 2).constrain(prob, "two_per_region")
# The first 5 units cost 3 each, and the rest 2 each.
operating = throughput.lp_interp(prob, "operating", xp=[0, 5, 12], fp=[0, 15, 29])
prob += ((ship_cost[:, :, None] * flow).sum() + open_cost @ is_open + operating.sum()).item()
prob.solve()

assert (is_open.values.reshape(2, 3).sum(axis=1) == 2).all()

Formulations

Methods taking (prob, name) add auxiliary variables and constraints named {name}_{role}.

bigM must bound the spread of the values involved: |x| for abs, x for logical_clip, and the range along axis for the rest. Too small a bigM silently cuts off solutions. lp_interp constrains x to lie between the first and last breakpoints rather than clamping it.

Semantics

constrain(prob, name) adds each constraint as name(i,j,…) and returns the array of LpConstraints that the problem stores.

Sums, @ and dot take time linear in the contracted length. Other contractions, such as np.einsum and np.inner, run NumPy's object loops, which chain PuLP's + and take quadratic time. Reductions return 0-d lparrays, whose item() is the PuLP expression. Ufuncs that are not affine, such as np.maximum and np.abs, raise TypeError; the methods above linearise them.

Keep the lparray on the left of + and - against a single PuLP variable or expression: PuLP treats an array on its right as an iterable and sums it into one expression.

Upgrading from 1.x

Every bigM is required, and constrain returns the stored constraints.