Commit 875ae1e9 authored by Michael Keller's avatar Michael Keller
Browse files

Worked on possibly crazy idea

parent 3b22b98e
......@@ -46,7 +46,7 @@ def score(field):
score += scorePixel(field, (i, j))
return score
x0 = [0.125, 0.109375, 0.111328125, 0.111083984375, 0.111114501953125, 0.11111068725585938, 0.11111116409301758, 0.0970017569405692, 0.11287478038242885, 0.125, 0.109375, 0.111328125, 0.111083984375, 0.111114501953125, 0.11111068725585938, 0.11111116409301758, 0.0970017569405692, 0.11287478038242885, 0.125, 0.109375, 0.111328125, 0.111083984375, 0.111114501953125, 0.11111068725585938, 0.11111116409301758, 0.0970017569405692, 0.11287478038242885, 0.125, 0.109375, 0.111328125, 0.111083984375, 0.111114501953125, 0.11111068725585938, 0.11111116409301758, 0.0970017569405692, 0.11287478038242885, 0.125, 0.109375, 0.111328125, 0.111083984375, 0.111114501953125, 0.11111068725585938, 0.11111116409301758, 0.0970017569405692, 0.11287478038242885, 0.125, 0.109375, 0.111328125, 0.111083984375, 0.111114501953125, 0.11111068725585938, 0.11111116409301758, 0.0970017569405692, 0.11287478038242885, 0.125, 0.109375, 0.111328125, 0.111083984375, 0.111114501953125, 0.11111068725585938, 0.11111116409301758, 0.0970017569405692, 0.11287478038242885, 0.125, 0.109375, 0.111328125, 0.111083984375, 0.111114501953125, 0.11111068725585938, 0.11111116409301758, 0.20987653732299805, 0.0, 0.0, 0.125, 0.109375, 0.111328125, 0.111083984375, 0.111114501953125, 0.11111068725585938, 0.11111116409301758, 0.20987653732299805]
x0 = [1.0, 1.0, 1.0, 4.590202509439396e-12, 4.590202509439396e-12, 4.590202509439396e-12, 4.590202509439396e-12, 4.590202509439396e-12, 4.590202509439396e-12, -1.0476920291210941e-13, 8.70963648832408e-12, -1.0476920291210941e-13, -1.0476920291210941e-13, -1.0476920291210941e-13, -1.0476920291210941e-13, -1.0476920291210941e-13, -1.0476920291210941e-13, -1.0476920291210941e-13, 1.0000000281552288, 1.0000000096375863, 0.9999999622071887, 4.590202509557559e-12, 4.590202509415376e-12, 4.590202509733667e-12, 4.590202509529271e-12, 4.590202509636449e-12, 4.590202509723023e-12, -1.0476920291210941e-13, -1.0476920291210941e-13, -1.0476920291210941e-13, -1.0476920291210941e-13, -1.0476920291210941e-13, -1.0476920291210941e-13, -1.0476920291210941e-13, -1.0476920291210941e-13, -1.0476920291210941e-13, -1.0476920291210941e-13, -1.0476920291210941e-13, -1.0476920291210941e-13, -1.0476920291210941e-13, -1.0476920291210941e-13, -1.0476920291210941e-13, -1.0476920291210941e-13, -1.0476920291210941e-13, -1.0476920291210941e-13, -1.0476920291210941e-13, -1.0476920291210941e-13, -1.0476920291210941e-13, -1.0476920291210941e-13, -1.0476920291210941e-13, -1.0476920291210941e-13, -1.0476920291210941e-13, -1.0476920291210941e-13, -1.0476920291210941e-13, -1.0476920291210941e-13, -1.0476920291210941e-13, -1.0476920291210941e-13, -1.0476920291210941e-13, -1.0476920291210941e-13, -1.0476920291210941e-13, -1.0476920291210941e-13, -1.0476920291210941e-13, -1.0476920291210941e-13, -1.0476920291210941e-13, -1.0476920291210941e-13, -1.0476920291210941e-13, -1.0476920291210941e-13, -1.0476920291210941e-13, -1.0476920291210941e-13, -1.0476920291210941e-13, -1.0476920291210941e-13, -1.0476920291210941e-13, -1.0476920291210941e-13, -1.0476920291210941e-13, -1.0476920291210941e-13, -1.0476920291210941e-13, -1.0476920291210941e-13, -1.0476920291210941e-13, -1.0476920291210941e-13, -1.0476920291210941e-13, -1.0476920291210941e-13]
print("Sum: ", sum(x0))
......
from scipy.optimize import linprog
import numpy as np
# Toy Example
......@@ -74,6 +75,8 @@ def main():
# build the LP
#
obj = []
rhs_leq = []
lhs_leq = []
rhs_eq = []
lhs_eq = []
......@@ -92,39 +95,44 @@ def main():
dp[str(n) + "|" + str(p) + "|" + str(c2) + "|" + str(c1)] = ind
ind += 1
# constraint that, for a given pixel, all variables imply the same pixel
# constrain that all neighborhood relationships sum to 1
for p in pixels:
for c1 in range(numb_crops):
neighbors = get_neighbors(p)
n0 = neighbors.pop()
for n in neighbors:
constr = [0 for _ in obj]
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, n0, c1, c2)] = 1
constr[get_ind(p, n, c1, c2)] = -1
lhs_eq.append(constr)
rhs_eq.append(0)
constr[get_ind(p, n, c1, c2)] = 1
lhs_eq.append(constr)
rhs_eq.append(1)
# constrain the total number of planted crops in a pixel
# constrain that all neighbors imply the same pixel
for p in pixels:
constr = [0 for _ in obj]
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
lhs_eq.append(constr)
rhs_eq.append(len(neighbors))
constr[get_ind(p, n, c1, c2)] = -1
constr[get_ind(p, n0, c1, c2)] = 1
lhs_eq.append(constr)
rhs_eq.append(0)
# follow the crop distribution
for cCheck in range(numb_crops):
# 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, cCheck, c2)] = 1
constr[get_ind(p, n0, cTest, c2)] = 1
lhs_eq.append(constr)
rhs_eq.append(dist[cCheck])
rhs_eq.append(dist[cTest])
......@@ -141,14 +149,13 @@ def main():
sol = list(opt.x)
for p in pixels:
neighbors = get_neighbors(p)
if len(neighbors) > 1:
neighbors.pop()
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)]
print(s, end=", ")
print("")
......
import numpy as np
# test if the quadratic constraint is spd
A = np.array([[9,7],[6,14]])
np.all(np.linalg.eigvals(A) >= 0)
\ No newline at end of file
# Import packages.
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
import cvxpy as cp
import numpy as np
# problem definition
x = 1
y = 3
R = lambda a, b: 1 if abs(a - b) <= 1 else 0
dist = [1,1,1]
# autocalculated variables
numb_crops = len(dist)
n = x * y * numb_crops
pixels = [(i, j) for i in range(x) for j in range(y)]
# Helper functions
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 get_ind(p, c):
return (p[0] * y + p[1]) * numb_crops + c
# the important stuff
def main():
#
# generate the SDP
#
A = []
b = []
C = np.zeros((n, n))
# build the objective function
for pA in pixels:
for pB in pixels:
for cA in range(numb_crops):
for cB in range(numb_crops):
if pB in get_neighbors(pA):
C[get_ind(pA, cA)][get_ind(pB, cB)] = -R(cA, cB) # solver minimizes
# force the matrix to be symmetric
# all neighborhoods imply the same pixel
for pA in pixels:
neighbors = get_neighbors(pA)
n0 = neighbors.pop()
for pB in neighbors:
for cA in range(numb_crops):
constr = np.zeros((n, n))
for cB in range(numb_crops):
constr[get_ind(pA, cA)][get_ind(pB, cB)] = 1
constr[get_ind(pA, cA)][get_ind(n0, cB)] = -1
A.append(constr)
b.append(0)
# pixels contain a sum of 1 crop
# for pA in pixels:
# constrA = np.zeros((n, n))
# constrB = np.zeros((n, n))
# for i in range(n):
# for j in range(n):
# constrA[i][j] = 0
# constrB[i][j] = 0
# neighbors = get_neighbors(pA)
# for pB in pixels:
# for cA in range(numb_crops):
# for cB in range(numb_crops):
# if pB in neighbors:
# constrA[get_ind(pA, cA)][get_ind(pB, cB)] = 1
# else:
# constrB[get_ind(pA, cA)][get_ind(pB, cB)] = 1
# A.append(constrA)
# A.append(constrB)
# b.append(1)
# b.append(0)
# constrain total connection count
A.append(np.ones((n, n)))
b.append(4)
#
# Build the solver and run it
#
xTx = cp.Variable((n, n), symmetric=True)
constraints = [0 <= xTx, xTx <= 1]
constraints += [xTx >> 0]
constraints += [
cp.trace(A[i] @ xTx) == b[i] for i in range(len(A))
]
prob = cp.Problem(cp.Minimize(cp.trace(C @ xTx)),
constraints)
prob.solve()
print("status:", prob.status)
print("The optimal value is", abs(prob.value))
print("A solution X is")
# print(xTx.value)
for v in xTx.value:
print(list(map(lambda x: round(x, 5), list(v))))
# find the solution
largest_entry = xTx.value[0]
for v in xTx.value:
if v[0] > largest_entry[0]:
largest_entry = v
sol = []
for v in largest_entry:
sol.append(v / largest_entry[0])
# print the solution
print("\n", largest_entry, "\n")
for p in pixels:
for c in range(numb_crops):
print(round(sol[get_ind(p, c)], 5), end=", ")
print("")
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