import numpy as np
import scipy

# ind: Spaltenindex zu Variablenname
ind = {}
ind["sa"] = 0
ind["sb"] = 1
ind["ab"] = 2
ind["ac"] = 3
ind["ad"] = 4
ind["bd"] = 5
ind["bf"] = 6
ind["cd"] = 7
ind["ce"] = 8
ind["dc"] = 9
ind["de"] = 10
ind["df"] = 11
ind["ef"] = 12
ind["et"] = 13
ind["fe"] = 14
ind["ft"] = 15

n = len(ind)		# Anzahl der Variablen

# name: Variablenname zu Spaltenindex
name = [ "sa", "sb", "ab", "ac", "ad", "bd", "bf", "cd",
         "ce", "dc", "de", "df", "ef", "et", "fe", "ft" ]

# Schranken fuer die Variablen
bounds = [
    (0, 12),	# sa
    (0,  8),	# sb
    (0,  7),	# ab
    (0,  3),	# ac
    (0,  9),	# ad
    (0,  5),	# bd
    (0,  2),	# bf
    (0,  5),	# cd
    (0,  6),	# ce
    (0,  2),	# dc
    (0,  2),	# de
    (0, 11),	# df
    (0,  8),	# ef
    (0, 16),	# et
    (0,  6),	# fe
    (0, 11)	# ft
]

# Definition der Zielfunktion
c = np.zeros(n)
c[ind["sa"]] = -1	# min -x_sa - x_sb
c[ind["sb"]] = -1

# Definition der Gleichheits-Nebenbedingungen
Aeq = np.zeros((6,n))
beq = np.zeros(6)

# 1. Nebenbedingung: x_s_a - x_a_b - x_a_c - x_a_d = 0
Aeq[0][ind["sa"]] = 1
Aeq[0][ind["ab"]] = -1
Aeq[0][ind["ac"]] = -1
Aeq[0][ind["ad"]] = -1

# 2. Nebenbedingung: x_s_b + x_a_b - x_b_d - x_b_f = 0
Aeq[1][ind["sb"]] = 1
Aeq[1][ind["ab"]] = 1
Aeq[1][ind["bd"]] = -1
Aeq[1][ind["bf"]] = -1

# 3. Nebenbedingung: x_a_c + x_d_c - x_c_d - x_c_e = 0
Aeq[2][ind["ac"]] = 1
Aeq[2][ind["dc"]] = 1
Aeq[2][ind["cd"]] = -1
Aeq[2][ind["ce"]] = -1

# 4. Nebenbedingung: x_a_d + x_b_d + x_c_d - x_d_c - x_d_e - x_d_f = 0
Aeq[3][ind["ad"]] = 1
Aeq[3][ind["bd"]] = 1
Aeq[3][ind["cd"]] = 1
Aeq[3][ind["dc"]] = -1
Aeq[3][ind["de"]] = -1
Aeq[3][ind["df"]] = -1

# 5. Nebenbedingung: x_c_e + x_d_e + x_f_e - x_e_f - x_e_t = 0
Aeq[4][ind["ce"]] = 1
Aeq[4][ind["de"]] = 1
Aeq[4][ind["fe"]] = 1
Aeq[4][ind["ef"]] = -1
Aeq[4][ind["et"]] = -1

# 6. Nebenbedingung: x_b_f + x_d_f + x_e_f - x_f_e - x_f_t = 0
Aeq[5][ind["bf"]] = 1
Aeq[5][ind["df"]] = 1
Aeq[5][ind["ef"]] = 1
Aeq[5][ind["fe"]] = -1
Aeq[5][ind["ft"]] = -1

res = scipy.optimize.linprog(c, None, None, Aeq, beq, bounds, method='highs-ds')

print(f'maximaler Flusswert = {res.fun * (-1)}')
print('Fluss:')
for i in range(n):
    print(f'    x_{name[i]} = {res.x[i]}')