*Very fast* Hadamard transform for single vector at a time.

    Parameters
    ----------
    Y: ndarray
        the n*2^k data to be 1d hadamard transformed
    ordering: bool, optional
        reorder from Walsh to sequence

    Returns
    -------
    H: ndarray
        hadamard transformed data.

    Examples
    --------
    from https://en.wikipedia.org/wiki/Hadamard_transform with normalisation

    >>> y = np.array([[1, 0, 1, 0, 0, 1, 1, 0]])
    >>> hadamard(y, ordering=False)
    array([[ 0.5 ,  0.25,  0.  , -0.25,  0.  ,  0.25,  0.  ,  0.25]])
    >>> hadamard(y, ordering=True)
    array([[ 0.5 ,  0.  ,  0.  ,  0.  , -0.25,  0.25,  0.25,  0.25]])

def hadamard(Y, ordering=True):
    """
    *Very fast* Hadamard transform for single vector at a time.

    Parameters
    ----------
    Y: ndarray
        the n*2^k data to be 1d hadamard transformed
    ordering: bool, optional
        reorder from Walsh to sequence

    Returns
    -------
    H: ndarray
        hadamard transformed data.

    Examples
    --------
    from https://en.wikipedia.org/wiki/Hadamard_transform with normalisation

    >>> y = np.array([[1, 0, 1, 0, 0, 1, 1, 0]])
    >>> hadamard(y, ordering=False)
    array([[ 0.5 ,  0.25,  0.  , -0.25,  0.  ,  0.25,  0.  ,  0.25]])
    >>> hadamard(y, ordering=True)
    array([[ 0.5 ,  0.  ,  0.  ,  0.  , -0.25,  0.25,  0.25,  0.25]])
    """
    # dot is a sum product over the last axis of a and the second-to-last of b.
    # Transpose - can specify axes, default transposes a[0] and a[1] hmm
    n_vectors, n_Y = Y.shape
    matching = (n_vectors, 2, int(n_Y / 2))
    H = np.array([[1, 1], [1, -1]]) / 2.  # Julia uses 2 and not sqrt(2)?
    steps = int(np.log(n_Y) / np.log(2))
    assert(2**steps == n_Y)  # required
    for _ in range(steps):
        Y = np.transpose(Y.reshape(matching), (0, 2, 1)).dot(H)
    Y = Y.reshape((n_vectors, n_Y))
    if ordering:
        Y = Y[:, _sequency(n_Y)]
    return Y