import numpy as np
np.set_printoptions(suppress=True)

c = np.array([2, 1.5, 0, 0, 0, 0, 0])
a = np.array([
     [12., 24., -1., 0., 0., 0., 0.],
     [16., 16., 0., -1., 0., 0., 0.],
     [30., 12., 0., 0., -1., 0., 0.],
     [1., 0., 0., 0., 0., 1., 0.],
     [0., 1., 0., 0., 0., 0., 1.]
    ])
b = np.array([120, 120, 120, 15, 15])

c = np.array([0,0,0,0,-27,-13,0,-15,-1,0,0,0,0,0])
b = np.array([14,20,4,7,13])
a = np.array([
    [0,0,0,0,1,1,0,1,1,-1,0,0,0,0],
    [0,0,0,0,1,0,0,1,0,0,1,0,0,0],
    [0,0,0,0,0,1,0,0,1,0,0,1,0,0],
    [0,0,0,0,1,1,0,0,0,0,0,0,1,0],
    [0,0,0,0,0,0,0,1,1,0,0,0,0,1]
])
alpha = 0.99995
eps = 0.00001
beta = 0.1



def make_mat(x, sigma, a):
    n, m = len(a), len(a[0])
    mat = np.array(a)
    at = np.transpose(mat)
    a1 = np.diagflat(-sigma / x)
    return np.block([[a1, at], [mat, np.zeros((n, n))]])

def solve(a, b, c):
    n, m = len(a), len(a[0])
    a_t = np.transpose(a)
    x = np.ones(m, dtype=float)
    pi = np.ones(n, dtype=float)
    sigma = np.ones(m, dtype=float)
    mu = beta * sum(x * sigma) / m
    while max(abs(x[j] * sigma[j] - mu) for j in range(m)) > eps:
        mu = beta * sum(x * sigma) / m
        fst = c - np.transpose(a) @ pi - mu / x 
        snd = b - a @ x
        B = np.concatenate((fst, snd))
        A = make_mat(x, sigma, a)
        # try:
        #     cho = np.linalg.cholesky(A)
        # except:
        #     print("Oh no!!!")
     
        x0 = np.linalg.solve(A, B)
        # print("sol :", x0)
        delta_x, delta_pi = x0[:m], x0[m:]
        delta_sigma = - np.array(sigma) - (sigma * delta_x - mu) / x #[(sigma[i] * delta_x[i] - mu) / x[i] for i in range(m)]
        theta_x = min((x[j] / -delta_x[j] for j in range(m) if delta_x[j] < 0), default=1)
        phi_s = min((sigma[i] / -delta_sigma[i] for i in range(m) if delta_sigma[i] < 0), default=1)
        theta = min(1, alpha * theta_x)
        phi = min(1, alpha * phi_s)
        print("total cost :", sum(c * x))
        x += theta * delta_x
        sigma += phi * delta_sigma
        pi += phi * delta_pi
    return x
x_sol = solve(a, b, c)
print(x_sol)