<SYSTEM_TASK:>
Solve the following problem using Python, implementing the functions described below, one line at a time
<END_TASK>
<USER_TASK:>
Description:
 def levinson_1d(r, order):
 """Levinson-Durbin recursion, to efficiently solve symmetric linear systems with toeplitz structure. Parameters --------- r : array-like input array to invert (since the matrix is symmetric Toeplitz, the corresponding pxp matrix is defined by p items only). Generally the autocorrelation of the signal for linear prediction coefficients estimation. The first item must be a non zero real. Notes ---- This implementation is in python, hence unsuitable for any serious computation. Use it as educational and reference purpose only. Levinson is a well-known algorithm to solve the Hermitian toeplitz equation: _ _ : : : : * : : : : _ * : _ with respect to a ( is the complex conjugate). Using the special symmetry in the matrix, the inversion can be done in O(p^2) instead of O(p^3). """

    r = np.atleast_1d(r)
    if r.ndim > 1:
        raise ValueError("Only rank 1 are supported for now.")

    n = r.size
    if n < 1:
        raise ValueError("Cannot operate on empty array !")
    elif order > n - 1:
        raise ValueError("Order should be <= size-1")

    if not np.isreal(r[0]):
        raise ValueError("First item of input must be real.")
    elif not np.isfinite(1 / r[0]):
        raise ValueError("First item should be != 0")

    # Estimated coefficients
    a = np.empty(order + 1, 'float32')
    # temporary array
    t = np.empty(order + 1, 'float32')
    # Reflection coefficients
    k = np.empty(order, 'float32')

    a[0] = 1.
    e = r[0]

    for i in range(1, order + 1):
        acc = r[i]
        for j in range(1, i):
            acc += a[j] * r[i - j]
        k[i - 1] = -acc / e
        a[i] = k[i - 1]

        for j in range(order):
            t[j] = a[j]

        for j in range(1, i):
            a[j] += k[i - 1] * np.conj(t[i - j])

        e *= 1 - k[i - 1] * np.conj(k[i - 1])

    return a, e, k