modeling: shaped Variable has no reductions; sum(xs) builds an O(n)-deep chain

[jkitchin]

Summary

A shaped Variable (m.continuous(name, shape=(n,))) supports indexing,
iteration and .shape/.size, but has no reduction methods. There is also
no module-level equivalent, so the only way to sum a vector of variables is the
Python builtin sum() — which builds a left-deep chain whose DAG depth is
O(n).

Measured on branch fix/lp-std-form-consolidation (worktree wt-stdform):

Variable has : __iter__ __len__ __getitem__ shape size

expression	result
xs.sum()	❌ AttributeError: 'Variable' object has no attribute 'sum'
X.sum() (2-D)	❌ AttributeError
np.sum(xs)	❌ TypeError: operand 'Variable' does not support ufuncs
sum(xs) (builtin)	✅ BinaryOp — but see depth below
xs[0] + xs[1]	✅ BinaryOp
xs @ xs	✅ MatMulExpression
len(xs)	✅ int

[a for a in dir(discopt) if 'sum' in a.lower()] is empty — no module-level
helper either.

The depth cost

sum() folds left, so each term adds a nesting level:

m  = Model()
xs = m.continuous("xs", shape=(200,), lb=0, ub=1)
e  = sum(xs)
# root = BinaryOp, DAG depth = 200

Measured DAG depth is exactly 200 for n = 200, i.e. depth grows 1:1 with
the vector length. A native .sum() emitting a single n-ary node would be depth
1 for the same model.

This is not only an ergonomics gap. Deep expression chains are the shape that
has already caused trouble elsewhere in the tree — cf. the deep-sum export
timeout addressed in aa8a2b9c — and construction cost in this layer is the
subject of #1215, which measured discopt at ~2.1x Pyomo's build time and ~1.9x
its memory. A left-deep chain of n-1 BinaryOp objects is n-1 Python
objects where one n-ary node would do, and every downstream traversal
(DAG compile, AD tape build, export) pays the depth in recursion.

Why np.sum cannot serve as the workaround

Expression.__array_ufunc__ = None (python/discopt/modeling/core.py:619) is a
deliberate NEP-13 opt-out and should stay (see #1234). Its consequence here is
that np.sum(xs) — which reduces via np.add.reduce — is refused, so NumPy
habits do not transfer and the builtin is the only route.

Suggested fix

Add reduction methods to the shaped Variable / Expression surface, each
emitting a single n-ary node rather than a fold:

• .sum(axis=None) — the primary ask
• .prod(axis=None), .mean(axis=None) — natural companions

and consider the shape manipulators NumPy users will also reach for, which are
likewise absent: .T, .reshape(), .flatten(). xs @ xs already works, so
matmul is supported while transpose is not, which is its own small
inconsistency.

If an n-ary sum node does not exist in the expression IR yet, a balanced
(pairwise) fold is a strictly better interim than the current left-deep one:
depth O(log n) instead of O(n) with no IR change.

Test

A regression test asserting xs.sum() exists, agrees numerically with
sum(xs), and produces DAG depth independent of n (or O(log n) for the
pairwise interim). Pinning the depth is what keeps a future refactor from
quietly reverting to a fold.

Verification

Probe counts: PROBES_EXECUTED 7 for the capability table, CHECKS_EXECUTED 1
for the depth measurement; no skipped comparisons.

Related: #1233 (intrinsics missing from the top-level namespace), #1234
(object-dtype ndarray arithmetic), #1215 (model construction cost).