Commit a6063461 authored by Michael Keller's avatar Michael Keller
Browse files

clean LP in gurobi

parent 6aeb910a
......@@ -14,9 +14,8 @@ dist = [1,1,1,1,1,1,1,1,1]
numb_crops = len(dist)
pixels = [(i, j) for i in range(X) for j in range(Y)]
#
# HELPERS
# Helper functions
#
dp = {}
def get_ind(pA, pB, c1, c2):
......@@ -33,131 +32,104 @@ def get_neighbors(p):
if 0 <= n[0] and 0 <= n[1] and n[0] < X and n[1] < Y:
neighbors.append(n)
return neighbors
#
# The Interesting Part
# The interesting stuff
#
def main():
numb_connections = 0
# build the (actual) objective function
ind = 0
# define the score function
obj = []
for p in pixels:
neighbors = get_neighbors(p)
for n in neighbors:
numb_connections += 1
for c1 in range(numb_crops):
for c2 in range(numb_crops):
if get_ind(p, n, c1, c2) is None:
obj.append(R(c1, c2) + R(c2, c1)) # the solver MAXimizes
# speeds up index queries
dp[str(p) + "|" + str(n) + "|" + str(c1) + "|" + str(c2)] = ind
dp[str(n) + "|" + str(p) + "|" + str(c2) + "|" + str(c1)] = ind
ind = 0
for pA in pixels:
neighbors = get_neighbors(pA)
for pB in neighbors:
for cA in range(numb_crops):
for cB in range(numb_crops):
if get_ind(pA, pB, cA, cB) is None:
obj.append(-(R(cA, cB) + R(cB, cA)))
# memoization
dp[str(pA) + "|" + str(pB) + "|" + str(cA) + "|" + str(cB)] = ind
dp[str(pB) + "|" + str(pA) + "|" + str(cB) + "|" + str(cA)] = ind
ind += 1
# setup the QP
m = gp.Model("qpc")
m.params.NonConvex = 2
x = m.addMVar(len(obj), lb=[0.0 for _ in obj], ub=[1.0 for _ in obj], name="x")
# the quadratic constraint / objective function (ensures a valid field)
O = np.zeros((len(obj), len(obj)))
for p in pixels:
neighbors = get_neighbors(p)
for n in neighbors:
for c1 in range(numb_crops):
for c2 in range(numb_crops):
for cA in range(numb_crops):
for cB in range(numb_crops):
O[get_ind(p, n, c1, c2)][get_ind(p, n, cA, cB)] = 1
# O = np.ones((len(obj), len(obj)))
m.addConstr(
x @ O @ x - sum(x) == 0,
"Quad_constr",
# sum([
# (sum(
# O[i][j] * x[i] * x[j] for i in range(len(obj))
# ) - x[j])
# for j in range(len(obj)) ])
# x @ O @ x - np.ones(len(obj)) @ x
)
# constrain that all pixels contain 1 total crop
pixel_sums = np.zeros((len(pixels), len(obj)))
pixel_sums_b = np.ones(len(pixels))
for i, p in enumerate(pixels):
n0 = get_neighbors(p).pop()
for c1 in range(numb_crops):
for c2 in range(numb_crops):
pixel_sums[i][get_ind(p, n0, c1, c2)] = 1
m.addConstr(pixel_sums @ x == pixel_sums_b, "pixel_sums_eq")
# all neighborhoods imply the same pixel
A = []
for p in pixels:
neighbors = get_neighbors(p)
obj = np.array(obj)
n = len(obj)
# Build the QP
m = gp.Model("lp")
x = m.addMVar(n, lb=[0.0 for _ in obj], ub=[1.0 for _ in obj])
m.setObjective(obj @ x)
# All neighbors have the same view of a pixel
neighbors_agree = []
neighbors_agree_b = []
for pA in pixels:
neighbors = get_neighbors(pA)
n0 = neighbors.pop()
for n in neighbors:
for c1 in range(numb_crops):
for pB in neighbors:
for cA in range(numb_crops):
constr = [0 for _ in obj]
for c2 in range(numb_crops):
constr[get_ind(p, n, c1, c2)] = -1
constr[get_ind(p, n0, c1, c2)] = 1
A.append(constr)
same_pixel_constr = np.zeros((len(A), len(obj)))
same_pixel_constr_b = np.zeros(len(A))
for i, l in enumerate(A):
for j, v in enumerate(l):
same_pixel_constr[i][j] = v
m.addConstr(same_pixel_constr @ x == same_pixel_constr_b, "same_pixel_eq")
for cB in range(numb_crops):
constr[get_ind(pA, pB, cA, cB)] = 1
constr[get_ind(pA, n0, cA, cB)] = -1
neighbors_agree.append(constr)
neighbors_agree_b.append(0)
neighbors_agree = np.array(neighbors_agree)
neighbors_agree_b = np.array(neighbors_agree_b)
m.addConstr(neighbors_agree @ x == neighbors_agree_b)
# All pixels contain a total of 1 crop
totals_one = []
totals_one_b = []
for pA in pixels:
constr = [0 for _ in obj]
pB = get_neighbors(pA).pop()
for cA in range(numb_crops):
for cB in range(numb_crops):
constr[get_ind(pA, pB, cA, cB)] = 1
totals_one.append(constr)
totals_one_b.append(1)
totals_one = np.array(totals_one)
totals_one_b = np.array(totals_one_b)
m.addConstr(totals_one @ x == totals_one_b)
# constrain the global crop distribution
global_crops = np.zeros((numb_crops, len(obj)))
global_crops_b = np.zeros(numb_crops)
# Enforce the global crop distribution
global_crops = []
global_crops_b = []
for cTest in range(numb_crops):
for p in pixels:
n0 = get_neighbors(p).pop()
for c2 in range(numb_crops):
global_crops[cTest][get_ind(p, n0, cTest, c2)] = 1
global_crops_b[cTest] = dist[cTest]
m.addConstr(global_crops @ x == global_crops_b, "global_crops_eq")
# constrain the score
objective = np.zeros(len(obj))
for i, v in enumerate(obj):
objective[i] = -v
m.setObjective(objective @ x)
# solve the QP
constr = [0 for _ in obj]
for pA in pixels:
pB = get_neighbors(pA).pop()
for cB in range(numb_crops):
constr[get_ind(pA, pB, cTest, cB)] = 1
global_crops.append(constr)
global_crops_b.append(dist[cTest])
global_crops = np.array(global_crops)
global_crops_b = np.array(global_crops_b)
m.addConstr(global_crops @ x == global_crops_b)
# ask the oracle
m.optimize()
score = 0
sol = list(m.getVars())
for p in pixels:
neighbors = get_neighbors(p)
n0 = neighbors.pop()
for c1 in range(numb_crops):
s = 0
for c2 in range(numb_crops):
s += sol[get_ind(p, n0, c1, c2)].X
score += sol[get_ind(p, n0, c1, c2)].X * obj[get_ind(p, n0, c1, c2)]
print(s, end=", ")
print("")
# the solution is
sol = []
for v in m.getVars():
sol.append(v.X)
print("\n\nScore: ", score)
for pA in pixels:
pB = get_neighbors(pA).pop()
for cA in range(numb_crops):
cA_amount = 0
for cB in range(numb_crops):
cA_amount += sol[get_ind(pA, pB, cA, cB)]
print(cA_amount, end=", ")
print("")
# for v in m.getVars():
# print('%s %g' % (v.VarName, v.X))
# print('Obj: %g' % obj.getValue())
main()
main()
\ No newline at end of file
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment