Commit 9994d964 authored by Michael Keller's avatar Michael Keller
Browse files

working version

parent a6063461
......@@ -5,11 +5,25 @@ from gurobipy import GRB
# Toy Example
X = 1
Y = 9
X = 3
Y = 3
R = lambda a, b: 1 if abs(a - b) <= 1 else 0
dist = [1,1,1,1,1,1,1,1,1]
# X = 15
# Y = 15
# R = lambda a, b: [[-1, 0, 0.5, -0.5, 1, 0, 0.5, 1, 0.5, 0],
# [0, -1, 0, 1, 0, 0, 0.5, 0.5, 0.5, 0],
# [0.5, 0, -1, 0, 0, 0, 0.5, 0.5, 0, -0.5],
# [-0.5, 1, 0, -1, 0, 0, -1, 0.5, 0.5, 0.5],
# [0, 0, 0, 0, -1, 0, 0, 0, 0, 0],
# [0, 0, 0, 0, 0, -1, 1, 0.5, 0, 0],
# [0.5, 0.5, 0.5, -0.5, 0, 0.5, -1, 0.5, 0.5, 0.5],
# [0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
# [0.5, 0.5, 0, 1, 0, 0, 0.5, 0.5, -1, 0],
# [0, 0, -0.5, 0.5, 0, 0, 1, 0.5, 0.5, -1]][a][b]
# dist = [23, 23, 23, 23, 23, 22, 22, 22, 22, 22]
# auto calculate some vars
numb_crops = len(dist)
pixels = [(i, j) for i in range(X) for j in range(Y)]
......@@ -60,7 +74,7 @@ def main():
n = len(obj)
# Build the QP
m = gp.Model("lp")
m = gp.Model("qpc")
x = m.addMVar(n, lb=[0.0 for _ in obj], ub=[1.0 for _ in obj])
m.setObjective(obj @ x)
......@@ -80,7 +94,7 @@ def main():
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)
m.addConstr(neighbors_agree @ x == neighbors_agree_b, "neighbors agree")
# All pixels contain a total of 1 crop
totals_one = []
......@@ -95,7 +109,7 @@ def main():
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)
m.addConstr(totals_one @ x == totals_one_b, "totals one")
# Enforce the global crop distribution
global_crops = []
......@@ -110,9 +124,41 @@ def main():
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)
m.addConstr(global_crops @ x == global_crops_b, "global crops")
# quadratic constraint (relaxed)
# O = np.zeros((n, n))
# for pA in pixels:
# neighbors = get_neighbors(pA)
# for pB in neighbors:
# for cA in range(numb_crops):
# for cB in range(numb_crops):
# for cAn in range(numb_crops):
# for cBn in range(numb_crops):
# O[get_ind(pA, pB, cA, cAn)][get_ind(pB, pA, cB, cBn)] += 1
# m.addConstr(x @ O @ x <= sum(x), "quad constraint")
# quadratic constraints (strict)
m.params.NonConvex = 2
for pA in pixels:
neighbors = get_neighbors(pA)
for pB in neighbors:
for cA in range(numb_crops):
for cB in range(numb_crops):
O = np.zeros((n, n))
for cAn in range(numb_crops):
for cBn in range(numb_crops):
O[get_ind(pA, pB, cA, cAn)][get_ind(pB, pA, cB, cBn)] = 1
m.addConstr(x @ O @ x == x[get_ind(pA, pB, cA, cB)], "quad constraint")
# ask the oracle
print("Ask Oracle")
m.optimize()
# the solution is
......@@ -121,12 +167,13 @@ def main():
sol.append(v.X)
for pA in pixels:
pB = get_neighbors(pA).pop()
neighbors = get_neighbors(pA)
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=", ")
for pB in neighbors:
for cB in range(numb_crops):
cA_amount += sol[get_ind(pA, pB, cA, cB)]
print(cA_amount / len(neighbors), end=", ")
print("")
......
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