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

Gradient Like LP approach

parent 24563cc0
import random
from scipy.optimize import linprog
import copy
import multiprocessing as mp
# Toy Example
# 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]
# larger example
# maxScoreGaertner is: 986.333333
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)
#
# HELPERS
#
# i is x-coord, j is y coord, c is the crop
def get_ind(i, j, c):
if i < 0 or j < 0:
return -1
ind = int((i * y + j) * len(dist) + c)
if ind >= x * y * numb_crops:
return -1
return ind
def get_pos(var):
cPos = var % numb_crops
xPos = int(((var - cPos) / numb_crops) / y)
yPos = ((var - cPos) / numb_crops) % y
return [xPos, yPos, cPos]
def scorePixel(field, pixel):
score = 0
p_x = pixel[0]
p_y = pixel[1]
for n1 in range(-1, 2):
for n2 in range(-1, 2):
# not a neighbor
if n1 == 0 and n2 == 0:
continue
# test valid neighbor
if not (0 <= (p_x + n1) and (p_x + n1) < x):
continue
if not (0 <= (p_y + n2) and (p_y + n2) < y):
continue
# update score
numb_crops = len(dist)
for c1 in range(numb_crops):
for c2 in range(numb_crops):
indA = get_ind(p_x, p_y, c1)
indB = get_ind(p_x + n1, p_y + n2, c2)
if field[indA] > 0 and field[indB] > 0:
score += field[indA] * field[indB] * (R(c1, c2) + R(c2, c1)) /2
return score
def score(field):
score = 0
for i in range(x):
for j in range(y):
score += scorePixel(field, (i, j))
return score
def print_field(x0):
for i in range(x):
for j in range(y):
for c in range(numb_crops):
if x0[get_ind(i, j, c)] == 1:
print(c, end=" ")
print("|", end="")
print("")
#
# END HELPERS
#
def iteration(x0):
# generate pixel subset
pixels = [(i, j) for i in range(x) for j in range(y)]
# generate objective function
obj = []
for p in pixels:
for c in range(numb_crops):
s = 0
for cn in range(numb_crops):
for n1 in range(-1, 2):
for n2 in range(-1, 2):
if n1 == 0 and n2 == 0:
continue
# test valid neighbor
if not (0 <= (p[0] + n1) and (p[0] + n1) < x):
continue
if not (0 <= (p[1] + n2) and (p[1] + n2) < y):
continue
ind = get_ind(p[0] + n1, p[1] + n2, cn)
if ind > -1:
s -= (R(c, cn) + R(cn, c)) * x0[ind] # the solver minimizes
obj.append(s)
# generate constraints for pixel sum
lhs_eq = []
for p in range(len(pixels)):
const = []
for q in range(len(obj)):
if int(q / numb_crops) == p:
const.append(1)
else:
const.append(0)
lhs_eq.append(const)
rhs_eq = [1 for _ in range(len(lhs_eq))]
# generate constraints for global sum
for c in range(numb_crops):
s = 0
const = [ 0 for _ in range(len(obj)) ]
for i, p in enumerate(pixels):
const[numb_crops * i + c] = 1
s += x0[get_ind(p[0], p[1], c)]
lhs_eq.append(const)
rhs_eq.append(s)
# generate variable bounds
bnd = [(0, 1) for _ in range(len(obj))]
# ask the oracle
opt = linprog(c=obj, A_eq=lhs_eq, b_eq=rhs_eq, bounds=bnd, method="revised simplex", options={ "tol": 0.001 }).x
# insert the solution
x1 = copy.deepcopy(x0)
for i, p in enumerate(pixels):
for c in range(numb_crops):
ind = get_ind(p[0], p[1], c)
x1[ind] += opt[i * numb_crops + c] * 0.1
x1[ind] -= x0[i * numb_crops + c] * 0.1
return x1
def main():
best_x0 = []
best_x0_score = 658.0
while True:
# generate random 1d field
field_stretched = []
for i, d in enumerate(dist):
field_stretched += [i] * d
random.shuffle(field_stretched)
# generate x0 vector
x0 = []
for i in range(x):
for j in range(y):
pixel = [0 for _ in range(len(dist))]
pixel[field_stretched[i * y + j]] = 1.0
x0 += pixel
iterations = 0
prev_score = -100 * x * y
print("Initial Score: ", score(x0), "Sum: ", sum(x0))
while iterations < 100 or (iterations < 200 and score(x0) > best_x0_score):
x0 = iteration(x0)
iterations += 1
scr = score(x0)
if scr > prev_score:
prev_score = scr
iterations = 0
print("Score: ", scr)
if score(x0) > best_x0_score:
best_x0 = copy.deepcopy(x0)
best_x0_score = score(x0)
print("Write to file")
f = open("best-score-grad-like.txt", "a")
f.write("Score: " + str(best_x0_score) + "\n x0 = " + str(best_x0) + "\n\n")
f.close()
print("Terminated")
print("Score: ", score(x0))
print("Best Score: ", best_x0_score)
# print(best_x0)
main()
\ No newline at end of file
{
"cSpell.words": [
"Birkhoff",
"Polytope"
]
}
\ No newline at end of file
......@@ -66,4 +66,8 @@
\usetikzlibrary{arrows}% To get more arrow heads
%% [OPT] Added by kelmich. Allows for text immediately after figures.
\usepackage{placeins}
\ No newline at end of file
\usepackage{placeins}
%% [OPT] Added by kelmich. Displays Pseudocode Algorithms
\usepackage{algorithm}
\usepackage{algpseudocode}
\ No newline at end of file
......@@ -40,7 +40,7 @@ their preferred crop:
\end{tikzpicture}
\end{figure}
\FloatBarrier
However, the question is how we find such a multi-way swap.
The question is how we find such a multi-way swap.
If we allowed swapping between neighboring pixels this would
indeed be tricky. If instead we restrict ourselves to only
swapping between non-neighboring pixels however, then the score
......@@ -48,7 +48,7 @@ function changes only linearly and we can use linear programming
to find the optimal swap.
Concretely, we setup the following linear program. Let $P$
be a set on non-neighboring pixels. Every pixel
be a set of $n$ non-neighboring pixels. Every pixel
gets $C$ variables that denote how much of each crop
is planted within the pixel:
\begin{align*}
......@@ -95,6 +95,77 @@ or formally, a solution $x^*$ must satisfy:
\forall c \in \Cps:\quad \sum_{p \in P} x^*_{p, c} &= \sum_{p \in P} x_{p, c}
\end{align*}
In theory an optimal solution to this linear program
could result in fractional amounts of crop being
planted in a pixel. However, the following theorem
demonstrates that there is always an optimal integer solution.
\begin{theorem}
The linear program for finding the optimal swap
between non neighboring pixels (described above)
always has an optimal integer solution.
\end{theorem}
\begin{proof}
The essence of our proof is that we first show
that when we remap $k$ crops between $k$ pixels
there must exist an integer solution, because
then we are in a Birkhoff polytope. Then we will
generalize the $k$ crop to $k$ pixels remapping
to a $m$ crop to $k$ pixels remapping.
Consider the case where we select $k$ non-neighboring
pixels that all contain a different crop. Then we ask
how to best remap these $k$ crops among the $k$ pixels.
In this case our solution space is the Birkhoff polytope.
Consider the following representation of our $x$ as a
matrix:
\begin{align*}
\begin{bmatrix}
x_{1, 1} & x_{1, 2} & \dots & x_{1, C}\\
x_{2, 1} & x_{2, 2} & \dots & x_{2, C}\\
\vdots & \vdots & & \vdots\\
x_{n, 1} & x_{n, 2} & \dots & x_{n, C}
\end{bmatrix}
\end{align*}
In order for a solution to be valid every row
must sum up to $1$ (a pixel contains a total
crop amount of $1$) and each column must sum up
to $1$ (we only have $1$ of each crop to distribute
among the pixels). This is the definition of the
Birkhoff polytope. Because the vertices of the
Birkhoff polytope are the permutation matrices,
these correspond exactly to our integer solutions.
Lastly, because the function we are maximizing
is linear, we know that at least one optimal
solution will be at a vertex. Because there is an
optimal solution at a vertex we know to be integer,
we know there must be an optimal integer solution.
What when the number of crop types is smaller
then the number of involved pixels? Even then
the linear program must have an optimal integer
solution.
TODO: finish proof
\end{proof}
\section{Method}
Now that we have an understanding of the
\ No newline at end of file
Now that we have a way to find the optimal multi-way swap
between pixels, how do we use it? We tried the following
procedure:
\begin{algorithm}
\caption{Random LP Iteration Algorithm}
\begin{algorithmic}
\State $i \gets 0$
\State $x_0 \gets \text{Random initial valid field}$
\While{$i \leq \text{Max Iterations}$}
\State $P \gets \{\text{random set of non-neighboring pixels}\}$
\State $opt \gets solveLP(P, x_0)$
\State $x_0 \gets updateField(opt, x0)$
\State $i \gets i + 1$
\EndWhile\\
\Return $x_0$
\end{algorithmic}
\end{algorithm}
No preview for this file type
......@@ -4,3 +4,5 @@
\contentsline {proof}{{Proof}{2}{}}{17}{proof.2}%
\contentsline {theorem}{{Theorem}{7.{2}}{}}{18}{theorem.7.2.2}%
\contentsline {proof}{{Proof}{3}{}}{19}{proof.3}%
\contentsline {theorem}{{Theorem}{8.{1}}{}}{22}{theorem.8.1.1}%
\contentsline {proof}{{Proof}{4}{}}{23}{proof.4}%
......@@ -20,8 +20,8 @@
\contentsline {subsection}{\numberline {7.2.2}The advanced method}{18}{subsection.7.2.2}%
\contentsline {chapter}{\chapternumberline {8}The Linear Programming Method}{21}{chapter.8}%
\contentsline {section}{\numberline {8.1}Problem Setup}{21}{section.8.1}%
\contentsline {section}{\numberline {8.2}Method}{22}{section.8.2}%
\contentsline {chapter}{\chapternumberline {9}Gradient Descent}{23}{chapter.9}%
\contentsline {chapter}{\chapternumberline {10}Conclusion}{25}{chapter.10}%
\contentsline {appendix}{\chapternumberline {A}Calculations Appendix}{27}{appendix.A}%
\contentsline {chapter}{Bibliography}{29}{appendix*.4}%
\contentsline {section}{\numberline {8.2}Method}{23}{section.8.2}%
\contentsline {chapter}{\chapternumberline {9}Gradient Descent}{25}{chapter.9}%
\contentsline {chapter}{\chapternumberline {10}Conclusion}{27}{chapter.10}%
\contentsline {appendix}{\chapternumberline {A}Calculations Appendix}{29}{appendix.A}%
\contentsline {chapter}{Bibliography}{31}{appendix*.4}%
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