Skip to content

GitLab

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

Proof savepoint

parent 9eb9696d
Expand all Hide whitespace changes
Inline Side-by-side
Scripts/python/GP/main.py 0 → 100644
View file @ 91340630
import cvxpy as cp
# 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 = [2,2,2,2,2,3,3,3,3,3]
X = 1
Y = 4
R = lambda a, b: 1
dist = [2,2]
def get_ind(p, c):
return numb_crops * p + c
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
numb_crops = len(dist)
pixels = [(i, j) for i in range(X) for j in range(Y)]
def main():
x = cp.Variable(shape=(X * Y * numb_crops,), pos=True, name="x")
constraints = []
# constrain crop amounts in a pixel
for p in range(X * Y):
monomial = 1
for c in range(numb_crops):
monomial *= x[p * numb_crops + c]
constraints.append(x[p * numb_crops + c] <= 1)
constraints.append(monomial == 1)
# constrain global crop amounts
for c in range(numb_crops):
constr = 0
for p in range(X * Y):
constr += x[p * numb_crops + c]
constraints.append(constr <= dist[c])
# objective function
obj = 0
for i, pA in enumerate(pixels):
neighbors = get_neighbors(pA)
for pB in neighbors:
j = neighbors.index(pB)
for cA in range(numb_crops):
for cB in range(numb_crops):
obj += x[get_ind(i, cA)] * x[get_ind(j, cB)] * R(cA, cB + 0.01)
constraints.append(obj <= 0.000001)
problem = cp.Problem(cp.Maximize(x[0]), constraints)
problem.solve(gp=True)
print("Optimal value:", problem.value)
print(x, ":", x.value)
main()
\ No newline at end of file
Scripts/python/LP_Iter/debug.py
View file @ 91340630
This diff is collapsed.
Scripts/python/LP_opt/qp.py
View file @ 91340630
import numpy as np
import gurobipy as gp
from gurobipy import GRB
# test if the quadratic constraint is spd
A = np.array([[9,7],[6,14]])
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]
# auto calculate some vars
numb_crops = len(dist)
pixels = [(i, j) for i in range(X) for j in range(Y)]
Rnew = lambda a, b: R(i, j) + 99999999999999999
# print(np.linalg.eigvals(A) >= -0.000001)
#
# Helper functions
#
def get_ind(p, c):
return numb_crops * p + c
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 stuff
#
def main():
# define the score function
obj = []
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)))
obj = np.array(obj)
n = len(obj)
# Build the QP
m = gp.Model("qp")
x = m.addMVar(n, lb=[0.0 for _ in obj], ub=[1.0 for _ in obj])
m.setObjective(obj @ x)
np.all(np.linalg.eigvals(A) >= 0)
\ No newline at end of file
main()
\ No newline at end of file
Scripts/python/QP/gurobi.py
View file @ 91340630
......@@ -20,19 +20,21 @@ np.set_printoptions(threshold=np.inf)
# dist = [4,5]
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 = [2,2,2,2,2,3,3,3,3,3]
# 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)
......@@ -41,9 +43,8 @@ pixels = [(i, j) for i in range(X) for j in range(Y)]
#
# Helper functions
#
dp = {}
def get_ind(pA, pB, c1, c2):
return dp.get(str(pA) + "|" + str(pB) + "|" + str(c1) + "|" + str(c2))
def get_ind(p, c):
return numb_crops * pixels.index(p) + c
def get_neighbors(p):
neighbors = []
......@@ -66,74 +67,33 @@ def main():
# define the score function
obj = []
ind = 0
obj = [[0 for _ in range(X * Y * numb_crops)] for _ in range(X * Y * numb_crops)]
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)))
if ind < 150:
print(ind, (pA, pB, cA, cB))
# memoization
dp[str(pA) + "|" + str(pB) + "|" + str(cA) + "|" + str(cB)] = ind
dp[str(pB) + "|" + str(pA) + "|" + str(cB) + "|" + str(cA)] = ind
ind += 1
v = 1 if R(cA, cB) + R(cB, cA) > 0.5 else 0 # (-R(cA, cB) - R(cB, cA)) / 2 + 1
obj[get_ind(pA, cA)][get_ind(pA, cA)] = 1
obj[get_ind(pA, cA)][get_ind(pB, cB)] = v
obj[get_ind(pB, cB)][get_ind(pA, cA)] = v
obj[get_ind(pB, cB)][get_ind(pB, cB)] = 1
# print(obj)
obj = np.array(obj)
n = len(obj)
# Build the QP
m = gp.Model("qpc")
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 pB in neighbors:
for cA in range(numb_crops):
constr = [0 for _ in obj]
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, "neighbors agree")
# All interpixel relations sum to 1
totals_one = []
totals_one_b = []
for pA in pixels:
neighbors = get_neighbors(pA)
for pB in neighbors:
constr = [0 for _ in obj]
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, "totals one")
m.setObjective(x @ obj @ x)
# all pixel amounts sum to 1
pixels_totals_one = []
pixels_totals_one_b = []
for pA in pixels:
pB = get_neighbors(pA).pop()
for p in pixels:
constr = [0 for _ in obj]
for cA in range(numb_crops):
for cB in range(numb_crops):
constr[get_ind(pA, pB, cA, cB)] = 1
for c in range(numb_crops):
constr[get_ind(p, c)] = 1
pixels_totals_one.append(constr)
pixels_totals_one_b.append(1)
pixels_totals_one = np.array(pixels_totals_one)
......@@ -143,111 +103,27 @@ def main():
# Enforce the global crop distribution
global_crops = []
global_crops_b = []
for cTest in range(numb_crops):
for c in range(numb_crops):
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
for p in pixels:
constr[get_ind(p, c)] = 1
global_crops.append(constr)
global_crops_b.append(dist[cTest])
global_crops_b.append(dist[c])
global_crops = np.array(global_crops)
global_crops_b = np.array(global_crops_b)
m.addConstr(global_crops @ x == global_crops_b, "global crops")
# Quadratic constraint, the linear version (C = 2)
quad_lin = []
quad_lin_b = []
for pA in pixels:
neighbors = get_neighbors(pA)
for pB in neighbors:
for cA in range(numb_crops):
for cB in range(numb_crops):
constr = [0 for _ in obj]
if cA == cB:
continue
constr[get_ind(pA, pB, cA, cA)] = -1
constr[get_ind(pA, pB, cB, cB)] = -1
constr[get_ind(pA, pB, cA, cB)] = 1
constr[get_ind(pA, pB, cB, cA)] = 1
quad_lin.append(constr)
quad_lin_b.append(0)
quad_lin = np.array(quad_lin)
quad_lin_b = np.array(quad_lin_b)
m.addConstr(quad_lin @ x == quad_lin_b, "quad lin")
# quadratic constraint (relaxed)
# 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))
# extra = 0
# 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
# O[get_ind(pB, pA, cB, cBn)][get_ind(pA, pB, cA, cAn)] += 1
# # for i in range(n):
# # if O[i][i] == 0:
# # for j in range(n):
# # if O[i][j] > 0:
# # O[i][i] = 1
# # extra += x[i]
# # break
# print(O)
# O[4][4] = 1
# O[4][8] = 1
# O[8][4] = 1
# O[4][12] = 1
# O[12][4] = 1
# O[8][8] = 1
# O[12][12] = 1
# O[1][1] = 1
# O[1][2] = 1
# O[1][3] = 1
# O[2][1] = 1
# O[2][2] = 1
# O[2][3] = 1
# O[3][1] = 1
# O[3][2] = 1
# O[3][3] = 1
# print(O)
# w, v = np.linalg.eigh(O)
# print(np.all(w >= -0.00000001))
# # print(v)
# # exit()
# m.addConstr(x @ O @ x <= 2 * x[get_ind(pA, pB, cA, cB)] + extra, "quad constraint")
# m.addConstr(extra <= 0.5, "quad constraint")
# quadratic constraints (strict)
# m.params.NonConvex = 2
# enforce no same crops next to each other
# 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
# # print(O)
# m.addConstr(x @ O @ x == x[get_ind(pA, pB, cA, cB)], "quad constraint")
# for c in range(numb_crops):
# constr = np.zeros((n, n))
# constr[get_ind(pA, c)][get_ind(pA, c)] = 1
# constr[get_ind(pB, c)][get_ind(pB, c)] = 1
# m.addConstr(x @ constr @ x <= 1, "quad constr")
# opt_sol = np.array([1,0,0,0,0,1,0,0,0,0,0,1])
# print("Optimal Score:", obj @ opt_sol)
# print(quad_lin @ opt_sol)
# print(opt_sol)
# print(quad_lin)
# ask the oracle
print("Ask Oracle")
......@@ -257,23 +133,10 @@ def main():
sol = []
for v in m.getVars():
sol.append(v.X)
for pA in pixels:
neighbors = get_neighbors(pA)
for cA in range(numb_crops):
cA_amount = 0
pB = get_neighbors(pA).pop()
for cB in range(numb_crops):
cA_amount += sol[get_ind(pA, pB, cA, cB)]
print(cA_amount, end=", ")
for p in pixels:
for c in range(numb_crops):
print(round(sol[get_ind(p, c)], 5), end=", ")
print("")
# print(sol)
print(obj @ sol)
# for i in range(len(obj)):
# if obj[i] * sol[i] < 0:
# print(i, obj[i], sol[i])
# print(quad_lin @ sol)
......
Scripts/python/QP/like-dislike.py
View file @ 91340630
......@@ -34,5 +34,7 @@ for x in range(len(dist)):
matrix[x][y] += (R(x, y) + R(y, x)) / 2 + 1
matrix[y][x] += (R(x, y) + R(y, x)) / 2 + 1
matrix[y][y] += (R(x, y) + R(y, x)) / 2 + 1
for row in matrix:
print(row)
w, v = LA.eig(matrix)
print(w)
Scripts/python/graph/main.py 0 → 100644
View file @ 91340630
import itertools
import copy
# sol: (1, 3, 7, 8, 9, 4, 6, 2, 0, 5) 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]
def score(path):
s = 0
for i in range(len(dist)-1):
s += R(path[i], path[i+1]) + R(path[i+1], path[i])
return s
best_path = [i for i in range(len(dist))]
best_score = -100000
base = [i for i in range(len(dist))]
for i, perm in enumerate(itertools.permutations(base)):
if i % 100000 == 0:
print(i)
s = score(perm)
if best_score < s:
best_path = copy.deepcopy(perm)
best_score = s
print(best_path, s)
\ No newline at end of file
Thesis/helpers.tex
View file @ 91340630
......@@ -129,5 +129,201 @@ statement with a condition on $R$ might be helpful:
met to ensure that swapping in either direction
even among neighboring pixels is never detrimental.
l
Let us consider the change in the score function
when swapping $c_\alpha, c_\beta \in \Cps$ between
two neighboring pixels $u$ and $v$. Let $F$ denote the field
before the swap and $F'$ after. Further let
$p_c$ denote the amount of crop $c$ within pixel $p$
in the field before the swap. $p'_c$ is defined analogously
for the field after the swap. Lastly let $N(p)$ denote
the set containing all of pixel $p$'s neighbors.
\begin{align*}
\Delta S
&= S(F', R) - S(F, R)\\
&= \sum_{x \in F}
\sum_{y \in N(x)}
\sum_{c_i \in \Cps}
\sum_{c_j \in \Cps}
R(c_i, c_j) \cdot x'_{c_i} \cdot y'_{c_j}\\
&- \sum_{x \in F}
\sum_{y \in N(x)}
\sum_{c_i \in \Cps}
\sum_{c_j \in \Cps}
R(c_i, c_j) \cdot x_{c_i} \cdot y_{c_j}\\
&= \sum_{x \in F}
\sum_{y \in N(x)}
\sum_{c_i \in \Cps}
\sum_{j \in \Cps}
R(c_i, c_j) \cdot (x'_{c_i} \cdot y'_{c_j} - x_{c_i} \cdot y_{c_j})
\end{align*}
Because we are only swapping between two pixels,
we can drop the terms where $p'_c = p_c$ as these are $0$.
We therefore only need to sum up the relationships
involving at least one of either $u$ or $v$:
\begin{align*}
\Delta S
&= \sum_{y \in N(u)}
\sum_{c_i \in \Cps}
\sum_{c_j \in \Cps}
(R(c_i, c_j) + R(c_j, c_i)) \cdot (u'_{c_i} \cdot y'_{c_j} - u_{c_i} \cdot y_{c_j})\\
&+ \sum_{y \in N(v) \setminus x}
\sum_{c_i \in \Cps}
\sum_{c_j \in \Cps}
(R(c_i, c_j) + R(c_j, c_i)) \cdot (v'_{c_i} \cdot y'_{c_j} - v_{c_i} \cdot y_{c_j})
\end{align*}
In a next step we separate out the crops we are swapping
from the sum: $c_\alpha$ and $c_\beta$. Of particular
interest in these long formulas are mainly the sets
over which we are summing.
\begin{align*}
\Delta S
&= \sum_{y \in N(u)} \
\sum_{c_i \in \{c_\alpha, c_\beta\}} \
\sum_{c_j \in \{c_\alpha, c_\beta\}} \
(R(c_i, c_j) + R(c_j, c_i)) \cdot (u'_{c_i} \cdot y'_{c_j} - u_{c_i} \cdot y_{c_j})\\
&+ \sum_{y \in N(u)} \
\sum_{c_i \in \{c_\alpha, c_\beta\}} \
\sum_{c_j \in \Cps \setminus \{c_\alpha, c_\beta\}} \
(R(c_i, c_j) + R(c_j, c_i)) \cdot (u'_{c_i} \cdot y'_{c_j} - u_{c_i} \cdot y_{c_j})\\
&+ \sum_{y \in N(v) \setminus x} \
\sum_{c_i \in \{c_\alpha, c_\beta\}} \
\sum_{c_j \in \{c_\alpha, c_\beta\}} \
(R(c_i, c_j) + R(c_j, c_i)) \cdot (v'_{c_i} \cdot y'_{c_j} - v_{c_i} \cdot y_{c_j}) \\
&+ \sum_{y \in N(v) \setminus x} \
\sum_{c_i \in \{c_\alpha, c_\beta\}} \
\sum_{c_j \in \Cps \setminus \{c_\alpha, c_\beta\}} \
(R(c_i, c_j) + R(c_j, c_i)) \cdot (v'_{c_i} \cdot y'_{c_j} - v_{c_i} \cdot y_{c_j})
\end{align*}
Which lets us simplify even further some variables that
stay the same:
\begin{align*}
\Delta S
&= \sum_{y \in N(u)} \
\sum_{c_i \in \{c_\alpha, c_\beta\}} \
\sum_{c_j \in \{c_\alpha, c_\beta\}} \
(R(c_i, c_j) + R(c_j, c_i)) \cdot (u'_{c_i} \cdot y'_{c_j} - u_{c_i} \cdot y_{c_j})\\
&+ \sum_{y \in N(u)} \
\sum_{c_i \in \{c_\alpha, c_\beta\}} \
\sum_{c_j \in \Cps \setminus \{c_\alpha, c_\beta\}} \
(R(c_i, c_j) + R(c_j, c_i)) \cdot ((u'_{c_i} - u_{c_i}) \cdot y_{c_j})\\
&+ \sum_{y \in N(v) \setminus x} \
\sum_{c_i \in \{c_\alpha, c_\beta\}} \
\sum_{c_j \in \{c_\alpha, c_\beta\}} \
(R(c_i, c_j) + R(c_j, c_i)) \cdot (v'_{c_i} \cdot y'_{c_j} - v_{c_i} \cdot y_{c_j}) \\
&+ \sum_{y \in N(v) \setminus x} \
\sum_{c_i \in \{c_\alpha, c_\beta\}} \
\sum_{c_j \in \Cps \setminus \{c_\alpha, c_\beta\}} \
(R(c_i, c_j) + R(c_j, c_i)) \cdot ((v'_{c_i} - v_{c_i}) \cdot y_{c_j})
\end{align*}
And from the way we shift crops we know:
\begin{align*}
u'_{c_\alpha} &= u_{c_\alpha} - \lambda & u'_{c_\beta} &= u_{c_\beta} + \lambda\\
v'_{c_\alpha} &= v_{c_\alpha} + \lambda & v'_{c_\beta} &= v_{c_\beta} - \lambda
\end{align*}
Therefore we can conclude that the second and fourth
term are linear functions in $\lambda$. Next we further
split up the $c_\alpha, c_\beta$ sums to be able to apply our
$p'$ definitions from above:
\begin{align*}
\Delta S
&= \sum_{y \in N(u)} \
\sum_{c_j \in \{c_\alpha, c_\beta\}} \
(R(c_\alpha, c_j) + R(c_j, c_\alpha)) \cdot (u'_{c_\alpha} \cdot y'_{c_j} - u_{c_\alpha} \cdot y_{c_j})\\