How to use the celerite.solver.get_kernel_value function in celerite

solver = solver()
    np.random.seed(seed)
    t = np.sort(np.random.rand(500))
    diag = np.random.uniform(0.1, 0.5, len(t))
    b = np.random.randn(len(t))

    with pytest.raises(RuntimeError):
        solver.log_determinant()
    with pytest.raises(RuntimeError):
        solver.dot_solve(b)

    solver.compute(
        alpha_real, beta_real, alpha_complex_real, alpha_complex_imag,
        beta_complex_real, beta_complex_imag, t, diag
    )
    K = get_kernel_value(
        alpha_real, beta_real, alpha_complex_real, alpha_complex_imag,
        beta_complex_real, beta_complex_imag, t[:, None] - t[None, :]
    )
    K[np.diag_indices_from(K)] += diag
    assert np.allclose(solver.solve(b).T, np.linalg.solve(K, b))

    b = np.random.randn(len(t), 5)
    assert np.allclose(solver.solve(b), np.linalg.solve(K, b))

def test_dot_L(with_general, seed=42):
    solver = celerite.CholeskySolver()
    np.random.seed(seed)
    t = np.sort(np.random.rand(5))
    b = np.random.randn(len(t), 5)
    yerr = np.random.uniform(0.1, 0.5, len(t))

    alpha_real = np.array([1.3, 0.2])
    beta_real = np.array([0.5, 0.8])
    alpha_complex_real = np.array([0.1])
    alpha_complex_imag = np.array([0.0])
    beta_complex_real = np.array([1.5])
    beta_complex_imag = np.array([0.1])

    K = get_kernel_value(
        alpha_real, beta_real, alpha_complex_real, alpha_complex_imag,
        beta_complex_real, beta_complex_imag, t[:, None] - t[None, :]
    )
    K[np.diag_indices_from(K)] += yerr**2

    if with_general:
        U = np.vander(t - np.mean(t), 4).T
        V = U * np.random.rand(4)[:, None]
        A = np.sum(U * V, axis=0) + 1e-8

        K[np.diag_indices_from(K)] += A
        K += np.tril(np.dot(U.T, V), -1) + np.triu(np.dot(V.T, U), 1)
    else:
        A = np.empty(0)
        U = np.empty((0, 0))
        V = np.empty((0, 0))

def _test_log_determinant(alpha_real, beta_real, alpha_complex_real,
                          alpha_complex_imag, beta_complex_real,
                          beta_complex_imag, solver, seed=42):
    solver = solver()
    np.random.seed(seed)
    t = np.sort(np.random.rand(5))
    diag = np.random.uniform(0.1, 0.5, len(t))

    solver.compute(
        alpha_real, beta_real, alpha_complex_real, alpha_complex_imag,
        beta_complex_real, beta_complex_imag, t, diag
    )
    K = get_kernel_value(
        alpha_real, beta_real, alpha_complex_real, alpha_complex_imag,
        beta_complex_real, beta_complex_imag, t[:, None] - t[None, :]
    )
    K[np.diag_indices_from(K)] += diag
    assert np.allclose(solver.log_determinant(), np.linalg.slogdet(K)[1])

def _test_log_determinant(alpha_real, beta_real, alpha_complex_real,
                          alpha_complex_imag, beta_complex_real,
                          beta_complex_imag, seed=42):
    solver = celerite.CholeskySolver()
    np.random.seed(seed)
    t = np.sort(np.random.rand(5))
    diag = np.random.uniform(0.1, 0.5, len(t))

    solver.compute(
        0.0, alpha_real, beta_real, alpha_complex_real, alpha_complex_imag,
        beta_complex_real, beta_complex_imag,
        np.empty(0), np.empty((0, 0)), np.empty((0, 0)),
        t, diag
    )
    K = get_kernel_value(
        alpha_real, beta_real, alpha_complex_real, alpha_complex_imag,
        beta_complex_real, beta_complex_imag, t[:, None] - t[None, :]
    )
    K[np.diag_indices_from(K)] += diag
    assert np.allclose(solver.log_determinant(), np.linalg.slogdet(K)[1])

@lapack_switch
def test_dot(solver, seed=42):
    solver = solver()
    np.random.seed(seed)
    t = np.sort(np.random.rand(300))
    b = np.random.randn(len(t), 5)

    alpha_real = np.array([1.3, 0.2])
    beta_real = np.array([0.5, 0.8])
    alpha_complex_real = np.array([0.1])
    alpha_complex_imag = np.array([0.3])
    beta_complex_real = np.array([0.5])
    beta_complex_imag = np.array([3.0])

    K = get_kernel_value(
        alpha_real, beta_real, alpha_complex_real, alpha_complex_imag,
        beta_complex_real, beta_complex_imag, t[:, None] - t[None, :]
    )
    x0 = np.dot(K, b)

    x = solver.dot(
        alpha_real, beta_real, alpha_complex_real, alpha_complex_imag,
        beta_complex_real, beta_complex_imag, t, b
    )
    print(x0)
    print(x)
    assert np.allclose(x0, x)

if with_general:
        U = np.vander(t - np.mean(t), 4).T
        V = U * np.random.rand(4)[:, None]
        A = np.sum(U * V, axis=0) + 1e-8
    else:
        A = np.empty(0)
        U = np.empty((0, 0))
        V = np.empty((0, 0))

    solver.compute(
        0.0, alpha_real, beta_real, alpha_complex_real, alpha_complex_imag,
        beta_complex_real, beta_complex_imag,
        A, U, V, t, diag
    )
    K = get_kernel_value(
        alpha_real, beta_real, alpha_complex_real, alpha_complex_imag,
        beta_complex_real, beta_complex_imag, t[:, None] - t[None, :]
    )
    K[np.diag_indices_from(K)] += diag

    if len(A):
        K[np.diag_indices_from(K)] += A
        K += np.tril(np.dot(U.T, V), -1) + np.triu(np.dot(V.T, U), 1)

    assert np.allclose(solver.solve(b).T, np.linalg.solve(K, b))

    b = np.random.randn(len(t), 5)
    assert np.allclose(solver.solve(b), np.linalg.solve(K, b))

"""
        Compute the value of the term for an array of lags

        Args:
            tau (array[...]): An array of lags where the term should be
                evaluated.

        Returns:
            The value of the term for each ``tau``. This will have the same
            shape as ``tau``.

        """
        tau = np.asarray(tau)
        (alpha_real, beta_real, alpha_complex_real, alpha_complex_imag,
         beta_complex_real, beta_complex_imag) = self.coefficients
        k = get_kernel_value(
            alpha_real, beta_real,
            alpha_complex_real, alpha_complex_imag,
            beta_complex_real, beta_complex_imag,
            tau.flatten(),
        )
        return np.asarray(k).reshape(tau.shape)