modeling: object-dtype ndarray times a scalar Expression raises an opaque ValueError

[jkitchin]

Summary

Mixing an object-dtype ndarray of expressions with a scalar Expression
raises an opaque NumPy-internals error instead of either working elementwise or
refusing clearly:

import numpy as np
from discopt import Model

m = Model()
x = m.continuous("x", lb=0.1, ub=2.0)
y = m.continuous("y", lb=0.1, ub=2.0)
arr = np.array([x, y], dtype=object)

arr * x     # ValueError: setting an array element with a sequence.
x * arr     # ValueError: setting an array element with a sequence.
arr + x     # ValueError: setting an array element with a sequence.

ValueError: setting an array element with a sequence gives no indication of
what the user did wrong or what to do instead.

Why it is inconsistent rather than simply unsupported

The neighbouring cases all work, which is what makes the failure surprising:

expression	result
x * np.array([1.0, 2.0]) (float ndarray ⊗ scalar expr)	✅ BinaryOp
np.array([1.0, 2.0]) * x	✅ BinaryOp
arr * arr (object ndarray ⊗ object ndarray)	✅ ndarray of exprs
A @ arr / np.dot(A, arr)	✅ ndarray of exprs
np.sum(arr)	✅ BinaryOp
arr + 1.0	✅ ndarray of exprs
arr * x (object ndarray ⊗ scalar expr)	❌ ValueError

So object arrays of expressions are supported for array⊗array and for
matrix products, and scalar-expression⊗ndarray is supported for numeric
ndarrays. Only the object-ndarray ⊗ scalar-expression corner falls through.

Mechanism

Expression.__array_ufunc__ = None (python/discopt/modeling/core.py:619) is
the documented NEP-13 opt-out: it makes NumPy return NotImplemented for
ndarray * Expression so the operation reflects to Expression.__rmul__, which
then consumes the whole array into a single expression node rather than letting
NumPy loop elementwise. That mechanism works and should be kept — it is why
x * np.array([1.0, 2.0]) yields one BinaryOp instead of an object array of
two.

The bug is that __mul__ / __rmul__ / __add__ assume any ndarray operand is
numeric data. With dtype=object the wrap fails deep inside NumPy while trying
to build the constant operand, surfacing as ValueError: setting an array element with a sequence. Confirmed by calling the dunder directly:

x.__rmul__(arr)   # ValueError: setting an array element with a sequence.

so the error originates inside discopt's own operator, not in NumPy dispatch.

Scope (deliberately narrow)

This does not affect the normal modeling path. m.continuous(name, shape=(n,)) returns a native Variable carrying .shape, not a NumPy object
array, and the shaped-variable operations are fine:

xs = m.continuous("xs", shape=(3,), lb=0, ub=1)
t  = m.continuous("t", lb=0, ub=1)
xs * t, t * xs, xs + t, xs * 2.0     # all OK -> BinaryOp

The failure requires a user to assemble np.array([...], dtype=object) by hand,
which is a plausible thing to do when collecting expressions from a
comprehension:

terms = np.array([some_expr(i) for i in range(n)], dtype=object)
scaled = terms * t          # ValueError

Suggested fix

In the arithmetic dunders, branch on operand.dtype == object for ndarray
operands. Either:

1. Elementwise — return an object ndarray whose element i is
   arr[i] <op> scalar. This is what NumPy would have produced had
   __array_ufunc__ not been set to None, and matches the already-working
   arr * arr behaviour; or
2. Refuse clearly — raise TypeError naming the object-dtype array and
   pointing at shaped variables (m.continuous(..., shape=...)) as the
   supported vectorized form.

(1) is more consistent with the table above. (2) is acceptable if object arrays
are considered out of contract, but the current opaque ValueError is not.

A regression test should pin whichever behaviour is chosen for all of
* / + - in both operand orders.

Verification

All results above were executed on branch fix/lp-std-form-consolidation
(worktree wt-stdform); the two probes reported PROBES_EXECUTED 7 and
PROBES_EXECUTED 5 respectively, with no skipped comparisons.

Related: #1233 (intrinsic functions missing from the top-level namespace).