Scipy: documentation of scipy.spatial.HalfspaceIntersection: wrong method to compute feasible point

[Ulfgard]

I am reading the docs of:
https://docs.scipy.org/doc/scipy/reference/generated/scipy.spatial.HalfspaceIntersection.html

the example states an algorithm for computing a feasible point for the halfspace intersection algorithm by setting up a linear problem. However, the algorithm uses

res = scipy.optimize.linprog(c, A_ub=A, b_ub=b)

to solve the problem, which additionally adds an implicit constraint x>0. This works in the example as the solution indeed produces positive values only. But in the code example below, flipping the one constraint in the example around to produce negative x-values will result in the wrong solution.

The correct way is to explicitly disable all box-constraints:
res = scipy.optimize.linprog(c, A_ub=A, b_ub=b, bounds = (None,None))

Reproducing code example:

from scipy.spatial import HalfspaceIntersection
from scipy.optimize import linprog
import numpy as np

halfspaces = np.array([[1, 0., 0.], # example uses -1
                       [0., -1., 0.],
                       [2., 1., -4.],
                       [-0.5, 1., -2.]])

norm_vector = np.reshape(np.linalg.norm(halfspaces[:, :-1], axis=1),
    (halfspaces.shape[0], 1))
c = np.zeros((halfspaces.shape[1],))
c[-1] = -1
A = np.hstack((halfspaces[:, :-1], norm_vector))
b = - halfspaces[:, -1:]
#wrong:
res = linprog(c, A_ub=A, b_ub=b)
print(res.x)
#correct
res = linprog(c, A_ub=A, b_ub=b, bounds = (None,None))
print(res.x)

Scipy/Numpy/Python version information:

1.4.1 1.17.4 sys.version_info(major=3, minor=7, micro=0, releaselevel='final', serial=0)

[rgommers]

Thanks for the report @Ulfgard, this is fixed now.