Sets & Indexing for Sparse Models

discopt’s modeling API has always supported dense, NumPy-style array variables (m.continuous("x", shape=(3, 4))) indexed by integer position. This notebook introduces the named-set / indexed layer for sparse, label-driven models — the style familiar from AMPL Fourer et al. [2003], GAMS, Pyomo [Hart et al., 2017], and JuMP [Dunning et al., 2017, Lubin et al., 2023].

The whole layer is pure-Python sugar: a named set is the authoritative index for variables/parameters/constraints, and everything desugars to the same flat model the solver already consumes. No solver or backend changes are involved.

Named sets

A Set holds an order-stable, de-duplicated collection of hashable members. Members may be scalars (dimen == 1) or fixed-arity tuples (dimen == N).

import discopt.modeling as dm

m = dm.Model("transport")
plants = m.set("plants", ["pitt", "sf", "nyc"])
markets = m.set("markets", ["a", "b", "c", "d"])
print(plants, "| dimen =", plants.dimen, "| len =", len(plants))

Set('plants', dimen=1, len=3) | dimen = 1 | len = 3

Set algebra

Sets compose with | (union), & (intersection), - (difference), * (cross product), and .where(pred) (filter) — the operators Pyomo [Hart et al., 2017] and JuMP [Dunning et al., 2017] use to express sparse index structure. Cross products are lazy, and filtering a product materializes only the members that pass the predicate (JuMP’s SparseAxisArray-via-condition).

# Sparse subset of the dense plants x markets grid.
links = (plants * markets).where(lambda p, k: (p, k) != ("nyc", "a"))
print("dense grid :", len(plants * markets))
print("sparse links:", len(links))
print(("nyc", "a") in links, ("pitt", "a") in links)

dense grid : 12
sparse links: 11
False True

Indexed variables

over=SET makes a variable indexed by set member. It is backed by a single flat variable, so ship["pitt", "a"] desugars to the ordinary indexing the solver already understands. Bounds may be a scalar, a dict keyed by member, or a callable fn(member).

supply = {"pitt": 20.0, "sf": 30.0, "nyc": 25.0}
ship = m.continuous("ship", over=links, lb=0, ub=1000)
cap = m.continuous("cap", over=plants, lb=0, ub=supply)  # dict bounds
print(ship)
print("ship[pitt, a] ->", ship["pitt", "a"])  # an indexing expression

IndexedVar('ship', over=(plants*markets)|where, len=11)
ship[pitt, a] -> ship[11][0]

Indexed constraints

m.constraint(SET, rule, name=) generates one constraint per member, named name[key]. Tuple members are unpacked into the rule’s arguments, and a rule may return dm.Skip to omit a member (Pyomo’s Constraint.Skip [Hart et al., 2017]).

demand = {"a": 15.0, "b": 20.0, "c": 20.0, "d": 10.0}
cost = {(p, k): 1.0 + (hash((p, k)) % 5) for (p, k) in links}

m.minimize(dm.sum(cost[(p, k)] * ship[p, k] for (p, k) in links))
m.constraint(plants,
             lambda p: dm.sum(ship[p, k] for k in markets if (p, k) in links) <= supply[p],
             name="supply")
m.constraint(markets,
             lambda k: dm.sum(ship[p, k] for p in plants if (p, k) in links) >= demand[k],
             name="demand")
print("variables:", len(m._variables), " sets:", len(m._sets))

variables: 2  sets: 2

Aggregation over sets

dm.sum and dm.prod accept generators over sets, or a callable plus over=. For dimen > 1 sets, tuple members are unpacked: dm.sum(lambda p, k: ..., over=links).

total_out = dm.sum(lambda p, k: ship[p, k], over=links)
print(total_out)

Σ[11 terms]

End-to-end: balanced transportation

A complete transportation model [Fourer et al., 2003] built from named sets, then solved. Values come back keyed by set member via IndexedVar.value.

t = dm.Model("transportation")
P = t.set("P", ["pitt", "sf"])
K = t.set("K", ["a", "b", "c"])
sup = {"pitt": 20.0, "sf": 30.0}
dem = {"a": 10.0, "b": 25.0, "c": 15.0}
c = {("pitt", "a"): 4, ("pitt", "b"): 6, ("pitt", "c"): 8,
     ("sf", "a"): 5, ("sf", "b"): 3, ("sf", "c"): 7}
x = t.continuous("x", over=P * K, lb=0, ub=1000)
t.minimize(dm.sum(c[(p, k)] * x[p, k] for p in P for k in K))
t.constraint(P, lambda p: dm.sum(x[p, k] for k in K) <= sup[p], name="supply")
t.constraint(K, lambda k: dm.sum(x[p, k] for p in P) >= dem[k], name="demand")
res = t.solve()
print(res.status, round(res.objective, 3))
print(x.value(res))

optimal 230.0
{('pitt', 'a'): 10.0, ('pitt', 'b'): 0.0, ('pitt', 'c'): 10.0, ('sf', 'a'): 0.0, ('sf', 'b'): 25.0, ('sf', 'c'): 5.0}

Dense vs sparse, and the linear fast path

Following JuMP’s tightest-container philosophy [Dunning et al., 2017], pick the representation that fits:

• Dense shape=(m, n) — rectangular data, integer position indexing.

• Indexed over a set — arbitrary labels, and sparse subsets via .where so you never allocate entries you won’t use.

When an indexed constraint family is affine in a single variable with a uniform sense, m.constraint(..., fast=True) (the default) emits it as one sparse-matrix call into the Rust builder instead of thousands of expression objects. This is a pure performance path — the resulting model and solution are identical, and any family that isn’t single-variable-affine transparently falls back.

Migrating from positional shape=

Positional (dense)	Named-set (sparse-friendly)
x = m.continuous("x", shape=(3,))	x = m.continuous("x", over=S)
x[0]	x["pitt"]
for i in range(3): m.subject_to(x[i] <= 1)	m.constraint(S, lambda i: x[i] <= 1)
dm.sum(x[i] for i in range(3))	dm.sum(x[i] for i in S)

Both produce the same flat model — named sets add labels, sparsity, and set algebra without changing the math.