pulp_and_numpy.md

This is an explanation for the issue Erwin Kalvelagen mentioned in his blog post.

When you add an arbitrary object to a numpy array, numpy tries to add that object to every element of the array which is known as broadcasting.

Here, for example, when you add a pulp.LpVariable to a numpy array, it is added to every element of the array:

In [36]: arr = np.array([1, 2, 3])

In [37]: x = pulp.LpVariable('x')

In [38]: arr + x
Out[38]: array([1*x + 1, 1*x + 2, 1*x + 3], dtype=object)

This behaviour is a little different when you use a list, instead of a numpy array:

In [39]: [1, 2, 3] + x

Out[39]: 1*x + 6

This is due to how the pulp.LpAffineExpression.addInPlace method is implemented:

    def addInPlace(self, other):
        if isinstance(other,int) and (other == 0): 
            return self
        if other is None: return self
        if isinstance(other,LpElement):
            self.addterm(other, 1)
        elif isinstance(other,LpAffineExpression):
            self.constant += other.constant
            for v,x in other.items():
                self.addterm(v, x)
        elif isinstance(other,dict):
            for e in other.values():
                self.addInPlace(e)
        elif (isinstance(other,list)
              or isinstance(other, Iterable)):
           for e in other:
                self.addInPlace(e)
        else:
            self.constant += other
        return self

Since the check for isinstance(other, list) returns True, it iterates over every element of the list, and adds them to x: x + 1 + 2 + 3 = x + 6.

This addInPlace method is called by the x variable's __radd__ method. __radd__ is a special method in Python that gets called when the left operand (in this case, the [1, 2, 3] list) does not now how to do addition with the right operand (the LpVariable x).

This does not happen with the numpy array because when you type arr + x, numpy assumes it knows how to add an object to an array, and starts the computation. Only when the broadcasting step is completed (i.e. when numpy has to add the first element of the array to x) the __radd__ method of x is called. If you switch the operands, you'll see something different:

In [55]: x + arr

Out[55]: 1*x + 6

Now, the result is inline with the list addition result because numpy's broadcasting doesn't happen. x sees that arr is an Iterable, and adds every element of the array to itself.

---

So why does something like prob += x + (arr - y) work but prob += (arr - y) doesn't? It is because of the addInPlace method again. Since arr - y returns an Iterable, when you add that to an LpVariable (x), every element of the Iterable is added to x, one by one, resulting in a single pulp expression. However, if the left operand is a numpy array, broadcasting happens and it propagates to every operation. arr - y yields an array, (arr - y) + x also yields an array, and so on. Since this does not result in a single expression, pulp cannot do optimization on it.

This indeed looks like unintended behaviour. I would follow Erwin's suggestion on not using numpy arrays with pulp or using it carefully.