savepoint

Scripts/python/QP/gurobi.py

+import numpy as np
+import gurobipy as gp
+from gurobipy import GRB
+
+
+
+# Toy Example
+X = 1
+Y = 9
+R = lambda a, b: 1 if abs(a - b) <= 1 else 0
+dist = [1,1,1,1,1,1,1,1,1]
+
+# auto calculate some vars
+numb_crops = len(dist)
+pixels = [(i, j) for i in range(X) for j in range(Y)]
+
+
+#
+# HELPERS
+#
+dp = {}
+def get_ind(pA, pB, c1, c2):
+    return dp.get(str(pA) + "|" + str(pB) + "|" + str(c1) + "|" + str(c2))
+
+def get_neighbors(p):
+    neighbors = []
+    for n1 in range(-1, 2):
+        for n2 in range(-1, 2):
+            if n1 == 0 and n2 == 0:
+                continue
+            n = (p[0] + n1, p[1] + n2)
+
+            if 0 <= n[0] and 0 <= n[1] and n[0] < X and n[1] < Y:
+                neighbors.append(n)
+    return neighbors
+            
+#
+# The Interesting Part
+#
+
+def main():
+
+    numb_connections = 0
+
+    # build the (actual) objective function
+    ind = 0
+    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 += 1
+    
+    # setup the QP
+    m = gp.Model("qpc")
+    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.setObjective(
+        sum([ 
+            (sum( 
+                O[i][j] * x[i] * x[j] for i in range(len(obj))
+            ) - x[j])
+            for j in range(len(obj)) ])
+        # 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)
+        n0 = neighbors.pop()
+        for n in neighbors:
+            for c1 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")
+    
+    # constrain the global crop distribution
+    global_crops = np.zeros((numb_crops, len(obj)))
+    global_crops_b = np.zeros(numb_crops)
+    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
+    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("")
+    
+    print("\n\nScore: ", score)
+
+    # for v in m.getVars():
+    #     print('%s %g' % (v.VarName, v.X))
+
+    # print('Obj: %g' % obj.getValue())
+
+main()

Scripts/python/QP/main.py

@@ -2,29 +2,151 @@
 import cvxpy as cp
 import numpy as np
 
-# Generate a random non-trivial quadratic program.
-m = 15
-n = 10
-p = 5
-np.random.seed(1)
-P = np.random.randn(n, n)
-P = P.T @ P
-q = np.random.randn(n)
-G = np.random.randn(m, n)
-h = G @ np.random.randn(n)
-A = np.random.randn(p, n)
-b = np.random.randn(p)
-
-# Define and solve the CVXPY problem.
-x = cp.Variable(n)
-prob = cp.Problem(cp.Minimize((1/2)*cp.quad_form(x, P) + q.T @ x),
-                 [G @ x <= h,
-                  A @ x == b])
-prob.solve()
-
-# Print result.
-print("\nThe optimal value is", prob.value)
-print("A solution x is")
-print(x.value)
-print("A dual solution corresponding to the inequality constraints is")
-print(prob.constraints[0].dual_value)
\ No newline at end of file
+# Toy Example
+# x = 1
+# y = 2
+# R = lambda a, b: 1 if abs(a - b) <= 1 else 0
+# dist = [1,1]
+
+# larger example
+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]
+
+# Cotrini Data (Sum = 26)
+# x = 15
+# y = 15
+# R = lambda a, b: [
+#     [-0.5, 1, 1, 1, 1, 1, -1, 0.5, 0.5, 1],
+#     [1, -0.5, -1, 1, 1, 0.5, -1, 1, 1, 1],
+#     [1, -1, -0.5, -1, 1, 0, 0.5, 0, 0, 0],
+#     [1, 1, -1, -0.5, 0, 0, -1, 1, -1, 0],
+#     [1, 1, 1, 0, -0.5, 0, 1, 1, 1, 1],
+#     [1, 0.5, 0, 0, 0, -0.5, 0, 0, 0, 0],
+#     [-1, -1, 0.5, -1, 1, 0, -0.5, 0, 0.5, 0],
+#     [0.5, 1, 0, 1, 1, 0, 0, -0.5, 1, 0],
+#     [0.5, 1, 0, -1, 1, 0, 0.5, 1, -0.5, 1],
+#     [1, 1, 0, 0, 1, 0, 0, 0, 1, -0.5]][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)]
+
+#
+# HELPERS
+#
+dp = {}
+def get_ind(pA, pB, c1, c2):
+    return dp.get(str(pA) + "|" + str(pB) + "|" + str(c1) + "|" + str(c2))
+
+def get_neighbors(p):
+    neighbors = []
+    for n1 in range(-1, 2):
+        for n2 in range(-1, 2):
+            if n1 == 0 and n2 == 0:
+                continue
+            n = (p[0] + n1, p[1] + n2)
+
+            if 0 <= n[0] and 0 <= n[1] and n[0] < x and n[1] < y:
+                neighbors.append(n)
+    return neighbors
+            
+#
+# END HELPERS
+#
+
+def main():
+
+    # build the objective function
+    ind = 0
+    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):
+                    if get_ind(p, n, c1, c2) is None:
+                        obj.append(-(R(c1, c2) + R(c2, c1))) # the solver minimizes
+
+                        # 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 += 1
+
+    print(len(obj))
+    
+    # build the constraints
+    A = []
+    b = []
+
+    # constrain that all neighborhood relationships sum to 1
+    for p in pixels:
+        neighbors = get_neighbors(p)
+        for n in neighbors:
+            constr = [0 for _ in obj]
+            for c1 in range(numb_crops):
+                for c2 in range(numb_crops):
+                    constr[get_ind(p, n, c1, c2)] = 1
+            A.append(constr)
+            b.append(1)
+    
+    # constrain that all neighbors imply the same pixel
+    for p in pixels:
+        neighbors = get_neighbors(p)
+        n0 = neighbors.pop()
+        for n in neighbors:
+            for c1 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)
+                b.append(0)
+    
+
+    # constrain the global crop distribution
+    for cTest in range(numb_crops):
+        constr = [0 for _ in obj]
+        for p in pixels:
+            n0 = get_neighbors(p).pop()
+            for c2 in range(numb_crops):
+                constr[get_ind(p, n0, cTest, c2)] = 1
+        A.append(constr)
+        b.append(dist[cTest])
+
+    # Define and solve the CVXPY problem.
+    P = np.ones((len(obj), len(obj)))
+    q = np.ones(len(obj))
+    x = cp.Variable(len(obj))
+    G = np.matrix(obj)
+    h = np.matrix([-1.5])
+    A = np.matrix(A)
+    b = np.array(b).T
+    print(np.all(np.linalg.eigvals(P) >= 0))
+    print(np.linalg.eigvals(P))
+    prob = cp.Problem(cp.Minimize(cp.quad_form(x, P) - q.T @ x),
+                    [G @ x <= h,
+                    0 <= x, x <= 1,
+                    A @ x == b])
+    prob.solve()
+
+    # Print result.
+    print("\nThe optimal value is", prob.value)
+    print("A solution x is")
+    print(list(x.value))
+    print("A dual solution corresponding to the inequality constraints is")
+    print(prob.constraints[0].dual_value)
+
+main()
\ No newline at end of file

Scripts/python/QP/test.py

+import numpy as np
+
+# Toy Example
+# x = 3
+# y = 3
+# R = lambda a, b: 1 if abs(a - b) <= 1 else 0
+# dist = [1 for _ in range(9)]
+
+# larger example
+x = 5
+y = 5
+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 = [3, 3, 3, 3, 3, 2, 2, 2, 2, 2]
+
+# Cotrini Data (Sum = 26)
+# x = 15
+# y = 15
+# R = lambda a, b: [
+#     [-0.5, 1, 1, 1, 1, 1, -1, 0.5, 0.5, 1],
+#     [1, -0.5, -1, 1, 1, 0.5, -1, 1, 1, 1],
+#     [1, -1, -0.5, -1, 1, 0, 0.5, 0, 0, 0],
+#     [1, 1, -1, -0.5, 0, 0, -1, 1, -1, 0],
+#     [1, 1, 1, 0, -0.5, 0, 1, 1, 1, 1],
+#     [1, 0.5, 0, 0, 0, -0.5, 0, 0, 0, 0],
+#     [-1, -1, 0.5, -1, 1, 0, -0.5, 0, 0.5, 0],
+#     [0.5, 1, 0, 1, 1, 0, 0, -0.5, 1, 0],
+#     [0.5, 1, 0, -1, 1, 0, 0.5, 1, -0.5, 1],
+#     [1, 1, 0, 0, 1, 0, 0, 0, 1, -0.5]][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)]
+
+#
+# HELPERS
+#
+dp = {}
+def get_ind(pA, pB, c1, c2):
+    return dp.get(str(pA) + "|" + str(pB) + "|" + str(c1) + "|" + str(c2))
+
+def get_neighbors(p):
+    neighbors = []
+    for n1 in range(-1, 2):
+        for n2 in range(-1, 2):
+            if n1 == 0 and n2 == 0:
+                continue
+            n = (p[0] + n1, p[1] + n2)
+
+            if 0 <= n[0] and 0 <= n[1] and n[0] < x and n[1] < y:
+                neighbors.append(n)
+    return neighbors
+
+def is_sem_pos_def(x):
+    return np.all(np.linalg.eigvals(x) >= -0.00000000001)
+            
+#
+# END HELPERS
+#
+
+
+
+def main():
+    #
+    # build the LP
+    #
+    obj = []
+
+    # build the objective function
+    ind = 0
+    for p in pixels:
+        neighbors = get_neighbors(p)
+        for n in neighbors:
+            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 minimizes
+
+                        # 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 += 1
+    
+    # test if positive semi definite
+    A = 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):
+                                A[get_ind(p, n, c1, c2)][get_ind(p, n, cA, cB)] = 1
+    
+    # B = np.zeros((numb_crops, numb_crops))
+    # for i in range(numb_crops):
+    #     for j in range(numb_crops):
+    #         B[i][j] = (R(i, j) + R(j, i)) / 2 + 5
+
+
+    # print(A)
+    # for l in A:
+    #     print(list(l))
+    print(is_sem_pos_def(A))
+
+    # seen = {}
+    # for i in range(1, 2000):
+    #     if is_sem_pos_def(np.ones((i,i))):
+    #         print(i)
+
+main()
+
+
+