I've read that integer programming is either very tricky or not possible with SciPy and that I probably need to use something like zibopt to do it in Python . But I really thought I could do it by creating one "is binary" constraint for each element in a vector being optimized by SciPy.
To do that, I utilized the closure trick from http://docs.python-guide.org/en/latest/writing/gotchas/#late-binding-closures and created one constraint function for each element, like so:
def get_binary_constraints(vector, indices_to_make_binary=None):
indices_to_make_binary = indices_to_make_binary or range(len(vector))
for i in indices_to_make_binary:
def ith_element_is_binary(vector, index=i):
return vector[index] == 0 or vector[index] == 1
yield ith_element_is_binary
test_vector = scipy.array([0.5, 1, 3])
constraints = list(get_binary_constraints(test_vector))
for constraint in constraints:
print constraint(test_vector)
which prints:
False
True
False
Then I modified get_binary_constraints for fmin_cobyla, whose constraints are a "sequence of functions that all must be >=0".
def get_binary_constraints(vector, indices_to_make_binary=None):
indices_to_make_binary = indices_to_make_binary or range(len(vector))
for i in indices_to_make_binary:
def ith_element_is_binary(vector, index=i):
return int(vector[index] == 0 or vector[index] == 1) - 1
yield ith_element_is_binary
which prints the following for the same test vector [0.5, 1, 3]:
-1
0
-1
So, only the 2nd value in the array will meet the condition >= 0.
Then, I set up a very simple optimization problem as follows:
from scipy import optimize
import scipy
def get_binary_constraints(vector, indices_to_make_binary=None):
indices_to_make_binary = indices_to_make_binary or range(len(vector))
for i in indices_to_make_binary:
def ith_element_is_binary(vector, index=i):
return int(vector[index] == 0 or vector[index] == 1) - 1
yield ith_element_is_binary
def objective_function(vector):
return scipy.sum(vector)
def main():
guess_vector = scipy.zeros(3)
constraints = list(get_binary_constraints(guess_vector))
result = optimize.fmin_cobyla(objective_function, guess_vector, constraints)
print result
if __name__ == '__main__':
main()
And this is what I get:
Return from subroutine COBYLA because the MAXFUN limit has been reached.
NFVALS = 1000 F =-8.614066E+02 MAXCV = 1.000000E+00
X =-2.863657E+02 -2.875204E+02 -2.875204E+02
[-286.36573349 -287.52043407 -287.52043407]
Before I go use R's LPSolve package or install zipobt for this, I'd really like to see if I can just use SciPy.
Am I doing something wrong, or is this just not possible to do in SciPy?
