Using le or ge with scalar left hand side creates unsized formula array

Question

When creating a formula with le or ge with a scalar left hand side, the returned formula has .shape == (), i.e. is an unsized object. This causes problems when adding a constraint with the formula.

I'm not sure if this is a bug.

prog = MathematicalProgram()
q = prog.NewContinuousVariables(1, 'q')
z = prog.NewContinuousVariables(2, 'z')
r = prog.NewContinuousVariables(rows=2, cols=3, name='r')
constraint = prog.AddConstraint(le(q, r[0,2] + 2*r[1,0])) # works
constraint2 = prog.AddConstraint(z[1] <= r[0,2] + 2*r[1,0]) # works
# constraint2 = prog.AddConstraint(le(z[1], r[0,2] + 2*r[1,0])) # fails
constraint2 = prog.AddConstraint(le([z[1]], r[0,2] + 2*r[1,0])) # works

formula = le(z[1], r[0,2] + 2*r[1,0])
print(formula.shape)
# > ()

Answer 1

TL;DR

For now, your available workarounds:

• Reformulate to be strictly vector-valued (as you had done)
• Pass the direct object value (using ndarray.item)
• Reshaping / unsqueeze the scalar array into a 1-d array

Will also post Drake issue.

Longform

Soonho's post is generally correct, with some minor corrections:

1. le(x, y) is another way of writing np.asarray(x) <= np.asarray(y). (See: drake#11171, numpy.vectorize docs)
2. While it is true that it isn't iterable, that's not the important part. The more important part is that <Formula> object output gets wrapped into a "scalar array" array(<Formula>, dtype=object).
3. As mentioned, this is a 0-dimensional array (empty shape); however, "empty array" may be a bit of a misnomer, since it has data. You can retrieve contents by using np.item()
4. Also correct. However, the nuance here is that np.vectorize will return "scalar arrays", and those in turn cannot be interpreted by our Python bindings using pybind11.

For point (4), I re-ran your example on the nightly drake-20201104-bionic.tar.gz, and got this error message:

TypeError: AddConstraint(): incompatible function arguments. The following argument types are supported:
    1. (self: pydrake.solvers.mathematicalprogram.MathematicalProgram, func: function, lb: numpy.ndarray[numpy.float64[m, 1]], ub: numpy.ndarray[numpy.float64[m, 1]], vars: numpy.ndarray[object[m, 1]], description: str = '') -> drake::solvers::Binding<drake::solvers::Constraint>
    2. (self: pydrake.solvers.mathematicalprogram.MathematicalProgram, arg0: pydrake.symbolic.Expression, arg1: float, arg2: float) -> drake::solvers::Binding<drake::solvers::Constraint>
    3. (self: pydrake.solvers.mathematicalprogram.MathematicalProgram, arg0: pydrake.symbolic.Formula) -> drake::solvers::Binding<drake::solvers::Constraint>
    4. (self: pydrake.solvers.mathematicalprogram.MathematicalProgram, constraint: drake::solvers::Constraint, vars: numpy.ndarray[object[m, 1]]) -> drake::solvers::Binding<drake::solvers::Constraint>
    5. (self: pydrake.solvers.mathematicalprogram.MathematicalProgram, formulas: numpy.ndarray[object[m, n], flags.f_contiguous]) -> drake::solvers::Binding<drake::solvers::Constraint>
Invoked with: <pydrake.solvers.mathematicalprogram.MathematicalProgram object at 0x7f6eb081cdf0>, array(<Formula "(z(1) <= (2 * r(1,0) + r(0,2)))">, dtype=object)

The most important thing to note is that overload (3) could have caught a raw scalar <Formula>, and (5) could have catch a 1- or 2-dimensional array, but neither are set up to catch array(<Formula>).

This is more-or-less a pydrake and/or pybind11 issue. I will file a Drake issue for now.

---

FWIW I've posted a NumPy issue asking about the right terminology for a 0-dim array:
https://github.com/numpy/numpy/issues/17744

Answer 2

1. le(x,y) is another way of writing np.array(x <= y).
2. x <= y forms a symbolic formula which is not iterable.
3. As a result, np.array(x <= y) creates an empty array. You can tell this by checking its shape (()) or trying np.array(x <= y)[0] (which gives you an error).
4. This is not Drake-symbolic specific. np.array(42) gives you an empty array.