Unexpected issue building model in parallel #1549
Description
Reported via the gpkit-users mailing list:
I'm trying to solve a sequence of bigger problems and thus want to generate the constraint set in parallel, using a multiprocessing pool in Python. For some reason, the .localsolve() routine to solve the resulting model (the problem is signomial indeed) throws an exception telling me to run localsolve instead (I am though, the constraints are explicitly signomial ones). If I do exactly the same but generate the constraints not in parallel but by a simple loop, everything works well. So I seemingly input exactly the same model, but once it works and once it does not. As soon as I use any of the constraints generated in parallel instead of its analogue generated by the loop, the localsolve breaks down.
Example code:
import numpy as np
import itertools
import multiprocessing as mp
from gpkit import VectorVariable, Model, SignomialsEnabled
T = 2
K = 2
r = np.array([2, 2])
A = np.array([[1, 0, 1], [0, 1, 1]])
x = [[np.ones(3).tolist()] + [np.zeros(3).tolist()] for t in range(T)]
v = VectorVariable((T + 1) * K, "v")
def generate_constraint(t, m):
b1 = A * x[t - 1][m]
b2 = np.power(np.transpose(np.broadcast_to(v[np.arange(t * K, (t + 1) * K)],
(b1.shape[1], b1.shape[0]))),
np.transpose(np.broadcast_to(r, (b1.shape[1], b1.shape[0]))) - b1)
summation_comps = np.prod(b2, axis=0)
with SignomialsEnabled():
constraint = [np.prod(np.power(v[np.arange((t - 1) * K, t * K)], r)) <=
sum(summation_comps)]
return constraint
if __name__ == '__main__':
index_list = [(1, 0), (1, 1), (2, 0), (2, 1)]
objective = 1.0 / (np.prod(np.power(v[np.arange(0, K)], r)))
# generate constraint set with simple loops
constraints = [v[i] <= v[i + K] for i in range(0, T * K)]
constraints += [v[T * K + k] == 1.0 for k in range(0, K)]
for index_pair in index_list:
result = generate_constraint(index_pair[0], index_pair[1])
constraints += result
# now generate constraint using a multiprocessing pool
constraints2 = [v[i] <= v[i + K] for i in range(0, T * K)]
constraints2 += [v[T * K + k] == 1.0 for k in range(0, K)]
pool = mp.Pool()
result2 = pool.starmap_async(generate_constraint, index_list)
constraints2 += list(itertools.chain.from_iterable(result2.get()))
pool.close()
pool.join()
pool.terminate()
# define corresponding optimization model
m = Model(objective, constraints)
m2 = Model(objective, constraints2)
# localsolve() works well for model m
sol = m.localsolve()
# ... however does not work for seemingly equivalent model m2
sol2 = m2.localsolve()
# now try the following: replace only one of the signomial inequalities in model "m" by the seemingly equivalent
# inequality from model m2:
constraints[8] = constraints2[8]
m3 = Model(objective, constraints)
# then it breaks down
sol3 = m3.localsolve()