Modeling Guide: Expressive vs. Fast Construction

Modeling Guide: Expressive vs. Fast Construction#

discopt provides two complementary APIs for building optimization models, and the expressive one has two idioms whose cost differs by an order of magnitude. All three appear below, because choosing between them matters more than choosing between discopt and another modelling layer.

Expressive, array-valued

Expressive, per-element

Fast API

Syntax

m.subject_to(dm.exp(x) + y <= b)

for i in ...: m.subject_to(...)

Matrix/vector (A, b, c, Q)

Objects built

O(constraint families)

O(rows)

none — one Rust call

Best for

anything at scale, linear or not

small models, irregular rows

large-scale LP/QP/MILP/MIQP

Nonlinear support

Full (exp, log, sin, bilinear, …)

Full

Linear and quadratic only

Hybrid

Yes — mix all three in one model

Yes

Yes

If you take one thing from this guide: write the body over the whole array. m.subject_to(dm.exp(x) + y <= b) with x of shape (n,) is one Python object describing n rows; the same constraint written as a loop over i is n objects. Measured across four model families at 100,000 rows, the two idioms differ by 6.9–11.1x in build-and-write time and 55–77x in memory — a larger spread than between discopt and Pyomo. The array-valued idiom puts discopt at parity with a compiled Rust modelling layer; the per-element idiom does not. See §4 below, and docs/dev/performance-plan.md §48/§50 for the measurements.

The per-element idiom is not deprecated — it is what AMPL-style indexed notation [Fourer et al., 2003] and Pyomo’s Constraint(index, rule=...) [Hart et al., 2017] express, it is often the clearest way to write an irregular family, and at small sizes the difference is invisible. It is simply not the one to reach for at scale.

Both APIs feed into the same solver pipeline: automatic problem classification, algebraic coefficient extraction, POUNCE (a pure-Rust Ipopt port) for the NLP solves, the in-house Rust simplex for LP/MILP, and Branch & Bound.

This notebook covers:

  1. The Expressive API — syntax, functions, examples, and its two idioms

  2. The Fast API — bulk construction for large models

  3. Hybrid mode — combining both in one model

  4. Performance comparison — when each one matters

  5. Decision guide — which to choose

import time

import discopt.modeling as dm
import numpy as np
import scipy.sparse as sp

Part 1: The Expressive API#

The expressive API uses Python operator overloading to build an expression DAG (directed acyclic graph). Every arithmetic operation creates a node; the DAG is later lowered to Rust for structure detection and for the evaluation tape the solver differentiates against.

Pattern#

m = dm.Model("name")
x = m.continuous("x", lb=0, ub=10)     # 1. Variables
m.minimize(x**2 + 1)                    # 2. Objective
m.subject_to(x >= 0.5)                  # 3. Constraints
result = m.solve()                       # 4. Solve

1.1 Variables#

Three variable types are available:

Method

Domain

Example

m.continuous(name, lb, ub, shape)

\([\text{lb}, \text{ub}]\)

Flow rates, temperatures

m.binary(name, shape)

\(\{0, 1\}\)

On/off decisions

m.integer(name, lb, ub, shape)

\(\{\text{lb}, \ldots, \text{ub}\}\)

Batch counts

 
m = dm.Model("variable_demo")

# Scalar variables
T = m.continuous("T", lb=300, ub=600)  # temperature (K)
on = m.binary("on")  # reactor on/off
n_batches = m.integer("n_batches", lb=0, ub=20)

# Array variables (shape tuple)
x = m.continuous("x", shape=(5,), lb=0, ub=1)  # 5-element vector
y = m.binary("y", shape=(3, 4))  # 3x4 binary matrix

print(f"Scalar T: shape={T.shape}, bounds=[{T.lb}, {T.ub}]")
print(f"Array x:  shape={x.shape}, {x.size} elements")
print(f"Matrix y: shape={y.shape}, {y.size} elements")

Scalar T: shape=(), bounds=[300.0, 600.0]
Array x:  shape=(5,), 5 elements
Matrix y: shape=(3, 4), 12 elements

1.2 Building Expressions#

Arithmetic on variables creates expression nodes in the DAG. All standard Python operators are supported, plus mathematical functions from the dm namespace.

m = dm.Model("expression_demo")
x = m.continuous("x", shape=(3,), lb=0, ub=5)

# Arithmetic operators
linear = 2 * x[0] + 3 * x[1] - x[2]  # linear combination
quadratic = x[0] ** 2 + x[1] * x[2]  # quadratic / bilinear
rational = x[0] / (1 + x[1])  # rational expression

# Mathematical functions
nonlinear = dm.exp(-x[0]) + dm.log(1 + x[1]) + dm.sqrt(x[2])
trig = dm.sin(x[0]) * dm.cos(x[1])

# Aggregations
total = dm.sum(x)  # sum of array
indexed = dm.sum(lambda i: x[i] ** 2, over=range(3))  # indexed sum

print(f"linear:    {linear}")
print(f"quadratic: {quadratic}")
print(f"nonlinear: {nonlinear}")

linear:    (((2 * x[3][0]) + (3 * x[3][1])) - x[3][2])
quadratic: ((x[3][0] ** 2) + (x[3][1] * x[3][2]))
nonlinear: ((exp(neg(x[3][0])) + log((1 + x[3][1]))) + sqrt(x[3][2]))

1.3 Constraints#

Comparisons on expressions produce Constraint objects:

m.subject_to(expr1 <= expr2)     # inequality
m.subject_to(expr1 >= expr2)     # inequality
m.subject_to(expr1 == expr2)     # equality

Constraints can be named for diagnostics:

m = dm.Model("constraint_demo")
x = m.continuous("x", shape=(3,), lb=0, ub=10)

m.minimize(dm.sum(x))

# Named constraints
m.subject_to(x[0] + x[1] <= 5, name="capacity")
m.subject_to(x[0] * x[1] >= 1, name="quality")
m.subject_to(dm.sum(x) == 7, name="balance")

# Loop-based constraints
for i in range(3):
    m.subject_to(x[i] >= 0.5, name=f"min_flow_{i}")

print(m.summary())

Model: constraint_demo
  Variables: 3 (3 continuous, 0 integer/binary)
  Constraints: 6
  Objective: minimize sum(x[3])
  Parameters: 0

1.4 Example: Geometric Programming (via log transform)#

Minimize the surface area of a box with fixed volume \(V = 1\):

\[ \min_{w, d, h} \; 2(wd + wh + dh) \quad \text{s.t.} \quad wdh = 1, \; w, d, h \in [0.1, 10] \]

The optimal solution is a cube: \(w^* = d^* = h^* = 1\), surface area \(= 6\).

m = dm.Model("box")
w = m.continuous("w", lb=0.1, ub=10)
d = m.continuous("d", lb=0.1, ub=10)
h = m.continuous("h", lb=0.1, ub=10)

m.minimize(2 * (w * d + w * h + d * h))
m.subject_to(w * d * h == 1)  # volume constraint

result = m.solve()
print(f"Status: {result.status}")
print(f"Surface area: {result.objective:.4f}  (expected: 6.0)")
print(f"w={result.x['w']:.4f}, d={result.x['d']:.4f}, h={result.x['h']:.4f}")

Status: optimal
Surface area: 6.0000  (expected: 6.0)
w=1.0000, d=1.0000, h=1.0000

1.5 Example: Reactor Design (nonlinear MINLP)#

Select a reactor configuration and optimize temperature to maximize yield. The reaction rate follows Arrhenius kinetics:

\[ \max_{T, y} \; y \cdot k_0 \, e^{-E_a / (RT)} \quad \text{s.t.} \quad 300 \leq T \leq 600, \; y \in \{0, 1\} \]

where \(y\) is an on/off indicator for the reactor.

k0 = 1e6  # pre-exponential factor
Ea = 5000.0  # activation energy / R (in K)

m = dm.Model("reactor")
T = m.continuous("T", lb=300, ub=600)
y = m.binary("y")  # reactor on/off

# Maximize rate (= minimize negative rate)
rate = k0 * dm.exp(-Ea / T)
m.minimize(-(y * rate))

# Reactor must be on for temperature to matter
m.subject_to(T <= 300 + 300 * y)  # T = 300 when y=0

result = m.solve()
print(f"Status: {result.status}")
print(f"Reactor on: {result.x['y']:.0f}")
print(f"Temperature: {result.x['T']:.1f} K")
print(f"Rate: {-result.objective:.2f}")

Status: optimal
Reactor on: 1
Temperature: 600.0 K
Rate: 240.37

1.6 Example: Portfolio Optimization (QP via expressions)#

Minimize portfolio variance subject to a target return:

\[ \min_w \; w^\top \Sigma\, w \quad \text{s.t.} \quad \mu^\top w \geq r_{\min}, \; \mathbf{1}^\top w = 1, \; w \geq 0 \]
n_assets = 4
mu = np.array([0.12, 0.10, 0.07, 0.03])  # expected returns
# Covariance matrix
Sigma = np.array(
    [
        [0.04, 0.006, 0.002, 0.0],
        [0.006, 0.025, 0.004, 0.001],
        [0.002, 0.004, 0.01, 0.002],
        [0.0, 0.001, 0.002, 0.005],
    ]
)
r_min = 0.08

m = dm.Model("portfolio")
w = m.continuous("w", shape=(n_assets,), lb=0, ub=1)

# Variance = w' Sigma w (expressed element-wise)
variance = dm.sum(
    lambda i: dm.sum(lambda j: Sigma[i, j] * w[i] * w[j], over=range(n_assets)),
    over=range(n_assets),
)

m.minimize(variance)
m.subject_to(dm.sum(lambda i: mu[i] * w[i], over=range(n_assets)) >= r_min)
m.subject_to(dm.sum(w) == 1)

result = m.solve()
print(f"Status: {result.status}")
print(f"Variance: {result.objective:.6f}")
print(f"Weights: {np.round(result.x['w'], 4)}")
print(f"Return: {np.dot(mu, result.x['w']):.4f}  (target >= {r_min})")

Status: optimal
Variance: 0.006602
Weights: [0.2065 0.2187 0.4027 0.1721]
Return: 0.0800  (target >= 0.08)

1.7 Array-valued bodies — the idiom for scale#

Everything in §1.1–§1.6 works elementwise over an array. If x has shape (n,), then dm.exp(x) + y <= b is a single expression describing n rows, and m.subject_to registers it as one Constraint object rather than n.

The two forms below build the identical model — the .nl files they write are byte-for-byte the same, and they solve to the same objective. They differ only in how many Python objects exist along the way.

# The SAME nonlinear model, written both ways.
n = 4
rhs = 3.0 + np.arange(n) * 0.5

# (a) per-element: one Constraint object per row
m_elem = dm.Model("per_element")
x = m_elem.continuous("x", shape=(n,), lb=0.5, ub=4.0)
y = m_elem.continuous("y", shape=(n,), lb=0.5, ub=4.0)
for i in range(n):
    m_elem.subject_to(dm.exp(x[i]) * 0.1 + y[i] <= float(rhs[i]), name=f"c{i}")
    m_elem.subject_to(x[i] * y[i] >= 1.0, name=f"d{i}")
m_elem.minimize(dm.sum(x) + dm.sum(y))

# (b) array-valued: one Constraint object per FAMILY
m_vec = dm.Model("vectorised")
x = m_vec.continuous("x", shape=(n,), lb=0.5, ub=4.0)
y = m_vec.continuous("y", shape=(n,), lb=0.5, ub=4.0)
m_vec.subject_to(dm.exp(x) * 0.1 + y <= rhs, name="c")   # n rows, one object
m_vec.subject_to(x * y >= 1.0, name="d")                  # n rows, one object
m_vec.minimize(dm.sum(x) + dm.sum(y))

print(f"per-element: {len(m_elem._constraints):>2d} Constraint objects")
print(f"array-valued:{len(m_vec._constraints):>2d} Constraint objects")

r_elem, r_vec = m_elem.solve(), m_vec.solve()
print(f"\nper-element  {r_elem.status}  objective = {r_elem.objective:.9f}")
print(f"array-valued {r_vec.status}  objective = {r_vec.objective:.9f}")
print(f"same answer: {abs(r_elem.objective - r_vec.objective) < 1e-9}")

per-element:  8 Constraint objects
array-valued: 2 Constraint objects

per-element  optimal  objective = 8.000000036
array-valued optimal  objective = 8.000000036
same answer: True

Two objects instead of eight, for a 4-row model; at 100,000 rows it is two instead of 200,000. That is the whole difference, and it is worth 6.9–11.1x in time and 55–77x in memory end to end (docs/dev/performance-plan.md §48, §50).

Practical notes:

  • The right-hand side can be a numpy array (one bound per row) or a scalar (broadcast to every row).

  • name="c" names the constraint family; the model holds one Constraint called c, and the writers expand it to rows c_0, c_1, … (visible in LP/MPS/GAMS output; .nl carries no names).

  • Reductions work too: dm.sum(x) <= 9.0 is one row, and dm.sum(X, axis=1) <= 9.0 reduces a 2-D family along one axis.

  • Matrix products work: A @ x <= b with a numpy A.

  • Integrality is unaffected — m.binary("z", shape=(n,)) participates in an array body like any other variable.

Part 2: The Fast Construction API#

For large-scale LP and QP models, the expressive API becomes a bottleneck: creating thousands of Python expression objects is slow. The fast API bypasses expression objects entirely, building constraints directly into the Rust expression arena via a single PyO3 call.

Methods#

Method

Purpose

m.add_linear_constraints(A, x, sense, b)

Bulk linear constraints: \(Ax \leq b\)

m.add_linear_objective(c, x)

Linear objective: \(\min\; c^\top x\)

m.add_quadratic_objective(Q, c, x)

Quadratic objective: \(\min\; \tfrac{1}{2} x^\top Q x + c^\top x\)

Key points:

  • A can be a scipy sparse matrix (any format) or dense numpy array — automatically converted to CSR

  • sense is "<=", "==", or ">="

  • b can be a vector or a scalar (broadcast to all rows)

  • Q must be symmetric; uses the convention \(f = \tfrac{1}{2} x^\top Q x + c^\top x\)

2.1 Example: Large LP (Transportation)#

Ship goods from 50 warehouses to 100 customers. This creates 150 constraints that we add in a single call.

n_wh, n_cust = 50, 100
n_vars = n_wh * n_cust  # 5000 decision variables

np.random.seed(42)
supply = np.random.uniform(50, 200, n_wh)
demand = np.random.uniform(10, 50, n_cust)
# Scale demand to match total supply
demand *= supply.sum() / demand.sum()
cost = np.random.uniform(1, 20, (n_wh, n_cust))

m = dm.Model("transport_fast")
x = m.continuous("x", shape=(n_vars,), lb=0)

# --- Build constraint matrix ---
# Supply: sum_j x[i,j] <= supply[i]  for each warehouse i
# Demand: sum_i x[i,j] >= demand[j]  for each customer j
rows_supply, cols_supply, data_supply = [], [], []
for i in range(n_wh):
    for j in range(n_cust):
        rows_supply.append(i)
        cols_supply.append(i * n_cust + j)
        data_supply.append(1.0)

A_supply = sp.csr_matrix((data_supply, (rows_supply, cols_supply)), shape=(n_wh, n_vars))

rows_demand, cols_demand, data_demand = [], [], []
for j in range(n_cust):
    for i in range(n_wh):
        rows_demand.append(j)
        cols_demand.append(i * n_cust + j)
        data_demand.append(1.0)

A_demand = sp.csr_matrix((data_demand, (rows_demand, cols_demand)), shape=(n_cust, n_vars))

# Add all constraints in two bulk calls
t0 = time.perf_counter()
m.add_linear_constraints(A_supply, x, "<=", supply, name="supply")
m.add_linear_constraints(A_demand, x, ">=", demand, name="demand")
m.add_linear_objective(cost.flatten(), x)
t_build = time.perf_counter() - t0

print(f"Model: {n_vars} variables, {n_wh + n_cust} constraints")
print(f"Build time: {t_build * 1000:.1f} ms")

result = m.solve()
print(f"Status: {result.status}")
print(f"Total cost: {result.objective:.2f}")
print(f"Solve time: {result.wall_time:.3f}s")

Model: 5000 variables, 150 constraints
Build time: 0.6 ms

Status: optimal
Total cost: 8761.51
Solve time: 0.221s

/Users/jkitchin/projects/discopt/python/discopt/debug/__init__.py:98: UserWarning: Variables with very large or infinite declared bounds: x[0] (lb=0, ub=9.999e+19 [default]), x[1] (lb=0, ub=9.999e+19 [default]), x[2] (lb=0, ub=9.999e+19 [default]), x[3] (lb=0, ub=9.999e+19 [default]), x[4] (lb=0, ub=9.999e+19 [default]). A bound marked [default] is the box applied to a column you declared with no bounds: 9.999e+19, which is FINITE (it sits just below the 1e+20 infinity sentinel), so the solve can return a certified 'optimal' sitting on that corner rather than 'unbounded'. NLP solvers may fail (NaN, iteration_limit) when bounds exceed ~1e15. Add tighter explicit bounds, e.g. m.continuous('x', lb=0, ub=1000).
  return fn(*args, **kwargs)

2.2 Example: Ridge Regression as QP#

Given data matrix \(X \in \mathbb{R}^{m \times n}\) and targets \(y \in \mathbb{R}^m\), ridge regression solves:

\[ \min_\beta \; \tfrac{1}{2} \|X\beta - y\|^2 + \tfrac{\lambda}{2} \|\beta\|^2 = \min_\beta \; \tfrac{1}{2} \beta^\top (X^\top X + \lambda I) \beta - (X^\top y)^\top \beta + \text{const} \]

This is a QP with \(Q = X^\top X + \lambda I\) and \(c = -X^\top y\).

# Generate synthetic regression data
np.random.seed(0)
m_samples, n_features = 200, 20
X_data = np.random.randn(m_samples, n_features)
beta_true = np.random.randn(n_features)
y_data = X_data @ beta_true + 0.5 * np.random.randn(m_samples)

lam = 1.0  # regularization strength

# QP coefficients
Q = sp.csr_matrix(X_data.T @ X_data + lam * np.eye(n_features))
c_vec = -(X_data.T @ y_data)

# Build and solve
mod = dm.Model("ridge")
beta = mod.continuous("beta", shape=(n_features,), lb=-100, ub=100)
mod.add_quadratic_objective(Q, c_vec, beta)

result = mod.solve()

# Compare with closed-form solution
beta_closed = np.linalg.solve(X_data.T @ X_data + lam * np.eye(n_features), X_data.T @ y_data)

print(f"Status: {result.status}")
print(f"Max |discopt - closed_form|: {np.max(np.abs(result.x['beta'] - beta_closed)):.2e}")
print(f"Solve time: {result.wall_time:.3f}s")

Status: optimal
Max |discopt - closed_form|: 2.95e-14
Solve time: 0.003s

2.3 Example: Facility Location (MILP)#

Choose \(p\) facility locations from \(n\) candidates to minimize total weighted distance to \(m\) customers.

\[ \min_{x, y} \; \sum_{i,j} d_{ij}\, x_{ij} \quad \text{s.t.} \quad \sum_j x_{ij} = 1 \; \forall i, \; x_{ij} \leq y_j \; \forall i,j, \; \sum_j y_j = p, \; x_{ij} \geq 0, \; y_j \in \{0,1\} \]
n_facilities = 10
n_customers = 30
p = 3  # number to open

np.random.seed(7)
# Random locations
fac_loc = np.random.rand(n_facilities, 2) * 100
cust_loc = np.random.rand(n_customers, 2) * 100
dist = np.sqrt(((cust_loc[:, None, :] - fac_loc[None, :, :]) ** 2).sum(axis=2))

n_x = n_customers * n_facilities  # assignment variables
n_y = n_facilities  # open/close variables
n_total = n_x + n_y

mod = dm.Model("facility")
# x_ij: fraction of customer i served by facility j (continuous)
x = mod.continuous("x", shape=(n_x,), lb=0, ub=1)
# y_j: whether facility j is open (binary)
y = mod.binary("y", shape=(n_y,))

# Objective: minimize total distance (linear in x)
mod.add_linear_objective(dist.flatten(), x)

# Assignment constraints: sum_j x_ij = 1 for each customer
rows_a, cols_a, data_a = [], [], []
for i in range(n_customers):
    for j in range(n_facilities):
        rows_a.append(i)
        cols_a.append(i * n_facilities + j)
        data_a.append(1.0)
A_assign = sp.csr_matrix((data_a, (rows_a, cols_a)), shape=(n_customers, n_x))
mod.add_linear_constraints(A_assign, x, "==", 1.0, name="assign")

# Linking: x_ij <= y_j (via expressions, since they involve two variable groups)
for i in range(n_customers):
    for j in range(n_facilities):
        mod.subject_to(x[i * n_facilities + j] <= y[j])

# Cardinality: sum_j y_j = p
mod.subject_to(dm.sum(y) == p, name="cardinality")

result = mod.solve(time_limit=30)
print(f"Status: {result.status}")
print(f"Total distance: {result.objective:.2f}")
print(f"Nodes: {result.node_count}, Time: {result.wall_time:.2f}s")

opened = [j for j in range(n_facilities) if result.x["y"][j] > 0.5]
print(f"Opened facilities: {opened}")

Status: optimal
Total distance: 708.58
Nodes: 1, Time: 0.04s
Opened facilities: [1, 2, 4]

Part 3: Hybrid Mode#

The real power of discopt is that you can mix both APIs in a single model. Use the fast API for the bulk linear/quadratic structure, and the expressive API for the few nonlinear constraints.

Example: Blending with Nonlinear Quality#

A refinery has a large block of linear flow-balance constraints (perfect for the fast API) plus a few nonlinear quality blending constraints (need expressions). Note in the summary below that the model reports all the rows — the bulk linear rows are materialized onto the nonlinear solve path rather than being visible only to the matrix engines.

n_streams = 50

mod = dm.Model("hybrid_blend")
flow = mod.continuous("flow", shape=(n_streams,), lb=0, ub=1000)

# --- Bulk linear constraints: flow balance ---
# Random sparse flow network.
np.random.seed(99)
A_balance = sp.random(40, n_streams, density=0.15, format="csr", random_state=99)

# Size the balance capacities around a reference operating point that already
# satisfies the two nonlinear constraints below (stream 1 carries the blend,
# stream 2 exceeds the conversion threshold). Drawing the right-hand side
# independently of A would leave the instance infeasible: a single row with a
# large coefficient on stream 2 is enough to contradict flow[2] >= 100*ln(2).
flow_ref = np.zeros(n_streams)
flow_ref[1] = 5.0
flow_ref[2] = 80.0
b_balance = A_balance @ flow_ref + np.random.uniform(10, 100, 40)
mod.add_linear_constraints(A_balance, flow, "<=", b_balance, name="balance")

# Linear objective
cost_vec = np.random.uniform(1, 10, n_streams)
mod.add_linear_objective(cost_vec, flow)

# --- A few nonlinear constraints via expressions ---
# Nonlinear blending: quality is a ratio
mod.subject_to(
    (3.0 * flow[0] + 5.0 * flow[1]) / (flow[0] + flow[1] + 0.01) >= 4.0,
    name="quality_blend",
)
# Exponential decay constraint
mod.subject_to(
    dm.exp(-flow[2] / 100) <= 0.5,
    name="conversion",
)

print(mod.summary())

result = mod.solve(time_limit=10)
print(f"\nStatus: {result.status}")
print(f"Objective: {result.objective:.4f}")
print(f"Solve time: {result.wall_time:.3f}s")

Model: hybrid_blend
  Variables: 50 (50 continuous, 0 integer/binary)
  Constraints: 42
  Objective: minimize 0
  Parameters: 0

Status: optimal
Objective: 0.0000
Solve time: 0.041s

Part 4: Performance Comparison#

All three ways of writing a model produce the same model, and it solves at the same speed. What differs — by two to three orders of magnitude on large instances — is what it costs to build the model object in the first place.

The comparison below is a linear model, so all three can express it. Note that the array-valued expressive form is not a compromise between the other two: it is the cheapest of the three and it is the only one of the cheap two that would still work if the body were dm.exp(x) instead.

sizes = [100, 500, 1000, 5000]
results_table = []

for n_con in sizes:
    n_var = 50
    np.random.seed(42)
    A = sp.random(n_con, n_var, density=0.1, format="csr", random_state=42)
    b = np.ones(n_con)
    c = np.ones(n_var)
    A_dense = A.toarray()

    # --- Fast API ---
    t0 = time.perf_counter()
    m_fast = dm.Model("fast")
    x_fast = m_fast.continuous("x", shape=(n_var,), lb=0, ub=100)
    m_fast.add_linear_constraints(A, x_fast, "<=", b)
    m_fast.add_linear_objective(c, x_fast)
    t_fast = time.perf_counter() - t0

    # --- Expressive, array-valued body (one Constraint object for all rows) ---
    t0 = time.perf_counter()
    m_vec = dm.Model("vec")
    x_vec = m_vec.continuous("x", shape=(n_var,), lb=0, ub=100)
    m_vec.subject_to(A_dense @ x_vec <= b)
    m_vec.minimize(dm.sum(c * x_vec))
    t_vec = time.perf_counter() - t0

    # --- Expressive, per-element (one Constraint object per row) ---
    t0 = time.perf_counter()
    m_expr = dm.Model("expr")
    x_expr = m_expr.continuous("x", shape=(n_var,), lb=0, ub=100)
    for row in range(n_con):
        lhs = dm.sum(lambda j: A_dense[row, j] * x_expr[j], over=range(n_var))
        m_expr.subject_to(lhs <= b[row])
    m_expr.minimize(dm.sum(lambda j: c[j] * x_expr[j], over=range(n_var)))
    t_expr = time.perf_counter() - t0

    results_table.append((n_con, t_fast, t_vec, t_expr))

hdr = f"{'Constraints':>12s} {'Fast (ms)':>10s} {'Array (ms)':>11s} {'Per-elem (ms)':>14s}"
print(hdr + f" {'per-elem / array':>17s}")
print("-" * len(hdr + " " * 18))
for n_con, t_fast, t_vec, t_expr in results_table:
    print(
        f"{n_con:>12d} {t_fast * 1000:>10.1f} {t_vec * 1000:>11.1f} "
        f"{t_expr * 1000:>14.1f} {t_expr / t_vec:>16.0f}x"
    )

 Constraints  Fast (ms)  Array (ms)  Per-elem (ms)  per-elem / array
--------------------------------------------------------------------
         100        0.1         0.1            5.7              103x
         500        0.1         0.0           79.2             2066x
        1000        0.2         0.1          107.8             2057x
        5000        0.7         0.2          504.6             2090x

Two things to read off that table.

The array-valued body is the cheapest arm, including against the matrix API — and unlike the matrix API it accepts exp, log, bilinear products and integrality. It is not a fallback for when the fast API will not do; it is the default worth reaching for.

The per-element loop is the expensive one, and it gets worse with size — 73x the array-valued form at 100 rows, ~600x at 5,000. That gap is the cost of one Python object per row, and no amount of Rust behind it can undo it: the objects have to exist before anything downstream sees them.

For small models (< 100 constraints) any of the three is fine — use whichever reads best. The distinction starts to matter in the hundreds and dominates in the thousands.

Part 5: Decision Guide#

Default: the Expressive API with array-valued bodies#

  • Any model with more than a few hundred rows, linear or nonlinear

  • Nonlinear functions (exp, log, sin, bilinear products) at any size

  • Data already in numpy arrays — the right-hand side can be an array

  • You want one line per constraint family, not per row

Use the per-element loop when:#

  • The model is small (< a few hundred rows) and the loop reads better

  • Rows are irregular — each one a different shape, not a family

  • You are porting AMPL/Pyomo-style indexed code [Fourer et al., 2003, Hart et al., 2017] and want a line-for-line translation first

Use the Fast API when:#

  • The model is large-scale LP/QP/MILP/MIQP and already in matrix form (A, b, c, Q) — network flows, regression, scipy sparse pipelines

  • No nonlinear terms are needed (it supports linear and quadratic only)

Use Hybrid mode when:#

  • Most constraints are linear (fast API) but a few are nonlinear (expressions)

  • You want the best of both: construction speed + modeling flexibility

Quick reference#

Task

Array-valued

Per-element

Fast

min x^2 + y^2

m.minimize(dm.sum(x**2 + y**2))

loop

m.add_quadratic_objective(I, zeros, v)

Ax <= b (1000 rows)

m.subject_to(A @ x <= b)

loop + m.subject_to()

m.add_linear_constraints(A, x, '<=', b)

exp(x) + log(y) <= 10

m.subject_to(dm.exp(x) + dm.log(y) <= 10)

loop

Not supported

x * y >= 1 (bilinear)

m.subject_to(x * y >= 1)

loop

Not supported

GDP / indicator

—

m.if_then(y, [...])

Not supported

API Reference Summary#

Model creation#

import discopt.modeling as dm
m = dm.Model("name")

Variables#

x = m.continuous("x", shape=(n,), lb=0, ub=100)
y = m.binary("y", shape=(k,))
z = m.integer("z", lb=0, ub=10)

Expressive objective and constraints#

m.minimize(expr)           # or m.maximize(expr)
m.subject_to(expr <= rhs, name="label")

Array-valued constraint bodies (the idiom for scale)#

x = m.continuous("x", shape=(n,), lb=0.5, ub=4.0)
b = np.arange(n) + 1.0

m.subject_to(x + 2 * y <= b,            name="c")   # n rows, ONE object
m.subject_to(dm.exp(x) + y <= 9.0,      name="d")   # scalar rhs broadcasts
m.subject_to(x * y >= 1.0,              name="e")   # bilinear, elementwise
m.subject_to(A @ x <= b,                name="f")   # matmul with a numpy A
m.subject_to(dm.sum(X, axis=1) <= 9.0,  name="g")   # reduce a 2-D family

Bodies may be any shape a variable can be declared with (shape=(2, 3) works; Variable has no .reshape, so declare the shape you want). Every math function in the list above is elementwise over an array body.

Mathematical functions#

dm.exp(x), dm.log(x), dm.sqrt(x)
dm.sin(x), dm.cos(x), dm.tan(x)
dm.abs(x), dm.log2(x), dm.log10(x)
dm.sum(x), dm.sum(lambda i: ..., over=range(n))

All of these are also importable from the top level (from discopt import tanh), including the hyperbolic, inverse-trig and erf/sigmoid/softplus families, matching the intrinsic set Pyomo exports from pyomo.environ [Hart et al., 2017]. NumPy ufuncs are deliberately refused on expressions (np.tanh(x) raises), so the discopt spelling is the one to reach for.

Reductions on a shaped expression#

xs.sum()          # one n-ary node -- NOT the builtin sum()'s left-deep chain
X.sum(axis=1)     # one row sum per row of a 2-D family
xs.mean(), xs.prod()
np.sum(xs)        # dispatches to .sum(); np.mean / np.prod likewise

xs.min(), xs.max()          # smallest / largest element (a bottleneck, a makespan)
np.min(xs), np.max(xs)      # dispatch to .min() / .max() likewise
dm.minimum(a, b, c, d)      # ELEMENT-WISE min of several operands
dm.maximum(*X[i, :])        # ... which is how you reduce one row

Prefer xs.sum() over the builtin sum(xs): the builtin folds left, so it builds n-1 binary nodes nested n deep, and every later traversal pays that depth. Variable still has no .reshape – declare the shape you want.

min/max are nonsmooth, so unlike sum they are a fold of binary nodes rather than one n-ary node – a balanced fold, ceil(log2(n)) deep instead of n-1 (docs/dev/performance-plan.md §64 records why the n-ary node buys no bound and would risk one). Two consequences to know about: .min(axis=k) is refused, because an axis reduction would have to stay array-valued and there is no array-valued min node to carry it (reduce a slice at a time, as above); and the .nl writer refuses a min/max node, so a model carrying one cannot be exported that way (AMPL calls such a model DNLP); the .nl reader does accept them, folding o11/o12 minlist/maxlist into binary nodes.

Fast construction#

m.add_linear_constraints(A, x, "<=", b, name="prefix")
m.add_linear_objective(c, x, constant=0.0, sense="minimize")
m.add_quadratic_objective(Q, c, x, constant=0.0, sense="minimize")

Solving#

result = m.solve(
    time_limit=3600,
    max_nodes=100000,
    gap_tolerance=1e-4,
    nlp_solver="pounce",     # "pounce" (default), "ipopt", or "simplex"
)
print(result.status, result.objective, result.x)