import numpy as np
from itertools import combinations


def get_jb(A):
    m,n = A.shape
    for jb in combinations(range(n), m):
        if abs(np.linalg.det(A[:, jb])) > 0.000001:
            return list(jb)
    raise ValueError('no basic plan')


def dual_simplex(A, b, c, d_lo, d_hi, Jb=None):
    if Jb is None:
        Jb = get_jb(A)
    m, n = A.shape
    Jn = list(set(range(n)) - set(Jb))
    # print('Jb',Jb)
    B = np.linalg.inv(A[:,list(Jb)])
    y = c[list(Jb)] @ B
    delta = y @ A - c

    J_plus = [x for x in Jn if delta[x] >= 0]
    J_minus = [x for x in Jn if delta[x] < 0]

    iter = 0
    while iter < 100:
        vector = np.zeros(n)
        vector[J_plus] = d_lo[J_plus]
        vector[J_minus] = d_hi[J_minus]
        vector[Jb] = B @ (b - np.sum(A[:,list(Jn)]*vector[Jn], axis=1))

        jk = None
        for i in Jb:
            if not d_lo[i] <= vector[i] <= d_hi[i]:
                jk = i
                break
        if jk is None:
            return vector, Jb

        myu = 1 if vector[jk] < d_lo[jk] else -1
        k = Jb.index(jk)
        m_vector = B[k] * myu
        myu_n = m_vector @ A[:, Jn]
        steps = np.array([-delta[j] / myu_n[idx]
                 if (j in J_plus and myu_n[idx] < 0) or (j in J_minus and myu_n[idx] > 0)
                 else np.inf
                 for idx, j in enumerate(Jn)])

        j_star_idx = np.argmin(steps)
        j_star = Jn[j_star_idx]
        min_step = steps[j_star_idx]
        if min_step == np.inf:
            return None, None
        delta[Jn] += min_step * myu_n
        delta[jk] += min_step * myu
        Jb[Jb.index(jk)] = j_star

        M = np.identity(B.shape[0])
        z = B @ A[:, j_star]
        zs = z[k]
        z[k] = -1
        z = -1 / zs * z
        M[:, k] = z
        B = np.dot(M, B)

        Jn[j_star_idx] = jk
        if j_star in J_plus:
            if myu == 1:
                J_plus.remove(j_star)
                J_plus.append(jk)
            else:
                J_plus.remove(j_star)
        else:
            if myu == 1:
                J_plus.append(jk)
            else:
                pass
        J_minus = list(set(Jn) - set(J_plus))
        iter += 1
    raise RuntimeError('too many iters')


if __name__ == '__main__':

    # A = np.array([
    #     np.array([2, 1, -1, 0, 0, 1], dtype=np.float),
    #     np.array([1, 0,  1, 1, 0, 1], dtype=np.float),
    #     np.array([0, 1,  0, 0, 1, 0], dtype=np.float)
    # ])
    # b = np.array([2, 5, 0], dtype=np.float)
    # c = np.array([3, 2, 0, 3, -2, -4], dtype=np.float)
    # d_low = np.array([0, -1, 2, 1, -1, 0], dtype=np.float)
    # d_high = np.array([2, 4, 4, 3, 3, 5], dtype=np.float)
    #
    # j_b= [3, 4, 5]
    # x0 = dual_simplex(A, b, c, d_low, d_high, j_b)
    # print(x0)
    # print(c @ x0[0])
    #
    # A = np.array([
    #     np.array([1, -5, 3, 1, 0, 0], dtype=np.float),
    #     np.array([4, -1, 1, 0, 1, 0], dtype=np.float),
    #     np.array([2,  4, 2, 0, 0, 1], dtype=np.float)
    # ])
    # b = np.array([-7, 22, 30], dtype=np.float)
    # c = np.array([7, -2, 6, 0, 5, 2], dtype=np.float)
    # d_low = np.array([2, 1, 0, 0, 1, 1], dtype=np.float)
    # d_high = np.array([6, 6, 5, 2, 4, 6], dtype=np.float)
    #
    # x, jb = dual_simplex(A, b, c, d_low, d_high)
    # print(x, jb)
    # print(c @ x)

    A = np.array([
        np.array([1, 7, 2, 0, 1, -1, 4], dtype=np.float),
        np.array([0, 5, 6, 1, 0, -3,-2], dtype=np.float),
        np.array([3, 2, 2, 1, 1,  1, 5], dtype=np.float)
    ])
    b = np.array([1, 4, 7], dtype=np.float)
    c = np.array([1, 2, 1, -3, 3, 1, 0], dtype=np.float)
    d_low = np.array([-1, 1, -2, 0, 1, 2, 4], dtype=np.float)
    d_high = np.array([3, 2, 2, 5, 3, 4, 5], dtype=np.float)

    x, jb = dual_simplex(A, b, c, d_low, d_high)
    print(x, jb)


    # A = np.array([
    #     [1, -3, 2, 0, 1, -1, 4, -1, 0],
    #     [1, -1, 6, 1, 0, -2, 2, 2, 0],
    #     [2, 2, -1, 1, 0, -3, 8, -1, 1],
    #     [4, 1, 0, 0, 1, -1, 0, -1, 1],
    #     [1, 1, 1, 1, 1, 1, 1, 1, 1]
    # ], dtype=np.float32)
    # b = np.array([3, 9, 9, 5, 9], dtype=np.float32)
    # c = np.array([-1, 5, -2, 4, 3, 1, 2, 8, 3], dtype=np.float32)
    #
    # d_low = np.array([0, 0, 0, 0, 0, 0, 0, 0, 0],dtype=float)
    # d_high = np.array([5, 5, 5, 5, 5, 5, 5, 5, 5], dtype=float)
    #
    # # j_base = np.array([3, 4, 5, 6, 7])
    # x, jb = dual_simplex(A, b, c, d_low, d_high)
    # print(x, x@c)