A = np.array([[2,-1,0,0,0,0,0,0,0,0],
              [-1,2,-1,0,0,0,0,0,0,0],
              [0,-1,2,-1,0,0,0,0,0,0],
              [0,0,-1,2,-1,0,0,0,0,0],
              [0,0,0,-1,2,-1,0,0,0,0],
              [0,0,0,0,-1,2,-1,0,0,0],
              [0,0,0,0,0,-1,2,-1,0,0],
              [0,0,0,0,0,0,-1,2,-1,0],
              [0,0,0,0,0,0,0,-1,2,-1],
              [0,0,0,0,0,0,0,0,-1,2]])


B = np.array([[1 if j <= i else 0 for j in range(10)] for i in range(10)])
q = np.array([1 for i in range(10)]).reshape((10,1))
d = np.array([i + 1 for i in range(10)]).reshape((10,1))



def get_x(l):
    return np.linalg.inv(A) @ (q - B @ l) / 2


def get_l(x, l, r):
    l = l + r * (B @ x - d)
#     print(np.maximum(np.zeros((10, 1)), l).shape)
#     print(l.shape)
    return np.maximum(np.zeros((10, 1)), l)


def f(x):
    print(B @ x - d)
    return (x.T @ A) @ x - q.T @ x


def ADMM(r=1, eps=0.01):
    x = np.random.random(size=(10,1))
    l = np.random.random(size=(10,1))
    stopper = 0
    while True:
        if stopper > 10:
            print('Stopped in %d iteration' % stopper)
            break
        l = get_l(x, l, r)
        x_next = get_x(x)
        if np.linalg.norm(x - x_next) < eps:
            print('Stopped in %d iteration' % stopper)
            break
        x = x_next
        stopper += 1
        if stopper % 5 == 0:
            print('%d iter f(x) = %f' % (stopper, f(x)))

    return x

ADMM(r=1, eps=0.1)