BUG: `MatrixVariable` can't broadcast like numpy

[Zeroto521]

Describe the bug

MatrixVariable can't broadcast between two matrices with different shapes.

To Reproduce

import numpy as np

a = np.arange(6).reshape(2, 3)  # shape is (2, 3)
# array([[0, 1, 2],
#        [3, 4, 5]])
b = np.array([[1], [4]])  # shape is (2, 1)
# array([[1],
#        [4]])

a <= b  # pass
# array([[ True,  True, False],
#        [ True,  True, False]])


from pyscipopt import Model

model = Model()
x = model.addMatrixVar((2, 3))
y = model.addMatrixVar(2)

# convert y shape to (2, 1)
x <= y[:, None]  # Error!
# Traceback (most recent call last):
#   line 22, in <module>
#     x <= y[:, None]  # Error!
#     ^^^^^^^^^^^^^^^
#   File "src/pyscipopt/matrix.pxi", line 36, in pyscipopt.scip.MatrixExpr.__le__
# IndexError: index 1 is out of bounds for axis 1 with size 1

Expected behavior

Did like numpy and supported broadcast.

Solution

It may be simplified like this. And use numpy to do broadcasting.
Since pyscipopt has already done the comparison core work between two expressions or variables.

def __le__(self, other):
    if isinstance(other, ...):  # type checking
        raise TypeError(f"Unsupported type {type(other)}")

    return super().__le__(other).view(MatrixExprCons)

PySCIPOpt/src/pyscipopt/matrix.pxi

Lines 27 to 40 in c681f94

	def __le__(self, other: Union[float, int, Variable, np.ndarray, 'MatrixExpr']) -> np.ndarray:
	expr_cons_matrix = np.empty(self.shape, dtype=object)
	if _is_number(other) or isinstance(other, Variable):
	for idx in np.ndindex(self.shape):
	expr_cons_matrix[idx] = self[idx] <= other
	elif isinstance(other, np.ndarray):
	for idx in np.ndindex(self.shape):
	expr_cons_matrix[idx] = self[idx] <= other[idx]
	else:
	raise TypeError(f"Unsupported type {type(other)}")
	return expr_cons_matrix.view(MatrixExprCons)

Screenshots

System

• OS: Windows 11
• Version: 24H2 (26100.4061)
• SCIP version: 9.2.3
• How did you install pyscipopt?: build from c681f94, base on 5.5.0

[Zeroto521]

A solution based on numpy.broadcast

import numpy as np

from pyscipopt import Model


def _richcmp(self, other, op):
    if op == 1:  # <=
        return self.__le__(other)
    elif op == 2:  # ==
        return self.__eq__(other)
    elif op == 5:  # >=
        return self.__ge__(other)
    else:
        raise NotImplementedError("Unsupported comparison operator")


def broadcast(self, other, op):
    self_other = np.broadcast(self, other)
    out = np.empty(self_other.shape, dtype=object)
    out.flat = [_richcmp(i, j, op) for i, j in self_other]
    return out


if __name__ == '__main__':
    a = np.arange(6).reshape(2, 3)
    b = np.array([[1], [4]])
    broadcast(a, b, 1)
    # array([[np.True_, np.True_, np.False_],
    #        [np.True_, np.True_, np.False_]], dtype=object)
    
    model = Model()
    x = model.addMatrixVar((2, 3))
    y = model.addMatrixVar(2)
    broadcast(x, y[:, None], 1)
    # array([[ExprCons(Expr({Term(x9): 1.0, Term(x15): -1.0}), None, 0.0),
    #         ExprCons(Expr({Term(x10): 1.0, Term(x15): -1.0}), None, 0.0),
    #         ExprCons(Expr({Term(x11): 1.0, Term(x15): -1.0}), None, 0.0)],
    #        [ExprCons(Expr({Term(x12): 1.0, Term(x16): -1.0}), None, 0.0),
    #         ExprCons(Expr({Term(x13): 1.0, Term(x16): -1.0}), None, 0.0),
    #         ExprCons(Expr({Term(x14): 1.0, Term(x16): -1.0}), None, 0.0)]],
    #       dtype=object)

[Joao-Dionisio]

Hey @Zeroto521, can you once again create a PR out of this? Finally managed to get some time to clear the backlog.

[Joao-Dionisio]

Closed via #1092