def assert_array_compare(comparison, x, y, err_msg='', verbose=True,
                         header=''):
    from numpy.core import asarray, isnan, any
    from numpy import isreal, iscomplex
    x = asarray(x)
    y = asarray(y)

    def isnumber(x):
        return x.dtype.char in '?bhilqpBHILQPfdgFDG'

    try:
        cond = (x.shape==() or y.shape==()) or x.shape == y.shape
        if not cond:
            msg = build_err_msg([x, y],
                                err_msg
                                + '\n(shapes %s, %s mismatch)' % (x.shape,
                                                                  y.shape),
                                verbose=verbose, header=header,
                                names=('x', 'y'))
            if not cond :
                raise AssertionError(msg)

        if (isnumber(x) and isnumber(y)) and (any(isnan(x)) or any(isnan(y))):
            # Handling nan: we first check that x and y have the nan at the
            # same locations, and then we mask the nan and do the comparison as
            # usual.
            xnanid = isnan(x)
            ynanid = isnan(y)
            try:
                assert_array_equal(xnanid, ynanid)
            except AssertionError:
                msg = build_err_msg([x, y],
                                    err_msg
                                    + '\n(x and y nan location mismatch %s, ' \
                                    '%s mismatch)' % (xnanid, ynanid),
                                    verbose=verbose, header=header,
                                    names=('x', 'y'))
            val = comparison(x[~xnanid], y[~ynanid])
        else:
            val = comparison(x,y)
        if isinstance(val, bool):
            cond = val
            reduced = [0]
        else:
            reduced = val.ravel()
            cond = reduced.all()
            reduced = reduced.tolist()
        if not cond:
            match = 100-100.0*reduced.count(1)/len(reduced)
            msg = build_err_msg([x, y],
                                err_msg
                                + '\n(mismatch %s%%)' % (match,),
                                verbose=verbose, header=header,
                                names=('x', 'y'))
            if not cond :
                raise AssertionError(msg)
    except ValueError:
        msg = build_err_msg([x, y], err_msg, verbose=verbose, header=header,
                            names=('x', 'y'))
        raise ValueError(msg)

def assert_array_compare(comparison, x, y, err_msg='', verbose=True,
                         header=''):
    from numpy.core import asarray
    x = asarray(x)
    y = asarray(y)
    try:
        cond = (x.shape==() or y.shape==()) or x.shape == y.shape
        if not cond:
            msg = build_err_msg([x, y],
                                err_msg
                                + '\n(shapes %s, %s mismatch)' % (x.shape,
                                                                  y.shape),
                                verbose=verbose, header=header,
                                names=('x', 'y'))
            assert cond, msg
        val = comparison(x,y)
        if isinstance(val, bool):
            cond = val
            reduced = [0]
        else:
            reduced = val.ravel()
            cond = reduced.all()
            reduced = reduced.tolist()
        if not cond:
            match = 100-100.0*reduced.count(1)/len(reduced)
            msg = build_err_msg([x, y],
                                err_msg
                                + '\n(mismatch %s%%)' % (match,),
                                verbose=verbose, header=header,
                                names=('x', 'y'))
            assert cond, msg
    except ValueError:
        msg = build_err_msg([x, y], err_msg, verbose=verbose, header=header,
                            names=('x', 'y'))
        raise ValueError(msg)

def check_conservation(initial_condition, flux):
 """!
 @brief Check whether the Riemman intergral of the solution is equal to the integral of the initial data, i.e. check if the numerical
 method preserved the conservation property of the linear advection, up to the floating point error.
 
 @param initial_condition One of the valid initial conditions. @see initial_conditions

 @param flux One of the valid fluxes. @see fluxes
 """
 
# Path to the hdf5 file containing results of computations
 computations_database = h5py.File(database_path, "r")
 
# Path to the data group containing the computational data
 group_path = initial_condition + "/" + flux
 
 for k in range (6, 17):
  dataset_initial_path = initial_condition + "/k = " + str(k) + " initial_data"
  dataset_path = group_path + "/k = " + str(k)
  
  h = 2**(-k)
  initial_data = asarray(computations_database[dataset_initial_path])
  computed_solution = asarray(computations_database[dataset_path])
# Conservation error is the absolute value of the difference between the integrals of the initial data and the computed solution
  conservation_error = abs(numpy.sum(initial_data) - numpy.sum(computed_solution))*h

# Criterion of preserving conservation must take into account the floating point error accumulation during the computation
# Experience shows that 0.1 is the upper bound for the floating point error for a conservative numerical method.
  if (conservation_error > 1e-1):
   print("ERROR: \n\t" + dataset_path + ": conservation is not preserved")
   print("\tConservation floating-point error: {0:.1e}".format(conservation_error))

def tensorsolve(a, b, axes=None):
    """Solves the tensor equation a x = b for x

    where it is assumed that all the indices of x are summed over in
    the product.

    a can be N-dimensional.  x will have the dimensions of A subtracted from
    the dimensions of b.
    """
    a = asarray(a)
    b = asarray(b)
    an = a.ndim

    if axes is not None:
        allaxes = range(0, an)
        for k in axes:
            allaxes.remove(k)
            allaxes.insert(an, k)
        a = a.transpose(allaxes)

    oldshape = a.shape[-(an-b.ndim):]
    prod = 1
    for k in oldshape:
        prod *= k

    a = a.reshape(-1, prod)
    b = b.ravel()
    res = solve(a, b)
    res.shape = oldshape
    return res

def tensorsolve(a, b, axes=None):
    """Solve the tensor equation a x = b for x

    It is assumed that all indices of x are summed over in the product,
    together with the rightmost indices of a, similarly as in
    tensordot(a, x, axes=len(b.shape)).

    Parameters
    ----------
    a : array-like, shape b.shape+Q
        Coefficient tensor. Shape Q of the rightmost indices of a must
        be such that a is 'square', ie., prod(Q) == prod(b.shape).
    b : array-like, any shape
        Right-hand tensor.
    axes : tuple of integers
        Axes in a to reorder to the right, before inversion.
        If None (default), no reordering is done.

    Returns
    -------
    x : array, shape Q

    Examples
    --------
    >>> from numpy import *
    >>> a = eye(2*3*4)
    >>> a.shape = (2*3,4,  2,3,4)
    >>> b = random.randn(2*3,4)
    >>> x = linalg.tensorsolve(a, b)
    >>> x.shape
    (2, 3, 4)
    >>> allclose(tensordot(a, x, axes=3), b)
    True

    """
    a = asarray(a)
    b = asarray(b)
    an = a.ndim

    if axes is not None:
        allaxes = range(0, an)
        for k in axes:
            allaxes.remove(k)
            allaxes.insert(an, k)
        a = a.transpose(allaxes)

    oldshape = a.shape[-(an-b.ndim):]
    prod = 1
    for k in oldshape:
        prod *= k

    a = a.reshape(-1, prod)
    b = b.ravel()
    res = wrap(solve(a, b))
    res.shape = oldshape
    return res

    def array_equal(cls, mx, other_mx):
        try:
            mx, other_mx = asarray(mx), asarray(other_mx)
        except:
            return False
        if mx.shape != other_mx.shape:
            return False

        both_equal = np.equal(mx, other_mx)
        both_nan   = np.logical_and(np.isnan(mx), np.isnan(other_mx))
        return bool(np.logical_or(both_equal, both_nan).all())

def ifftshift(x,axes=None):
    """
    Inverse of fftshift.

    Parameters
    ----------
    x : array_like
        Input array.
    axes : int or shape tuple, optional
        Axes over which to calculate.  Defaults to None which is over all axes.

    See Also
    --------
    fftshift

    """
    tmp = asarray(x)
    ndim = len(tmp.shape)
    if axes is None:
        axes = range(ndim)
    y = tmp
    for k in axes:
        n = tmp.shape[k]
        p2 = n-(n+1)/2
        mylist = concatenate((arange(p2,n),arange(p2)))
        y = take(y,mylist,k)
    return y

def fftshift(x,axes=None):
    """
    Shift zero-frequency component to center of spectrum.

    This function swaps half-spaces for all axes listed (defaults to all).
    If len(x) is even then the Nyquist component is y[0].

    Parameters
    ----------
    x : array_like
        Input array.
    axes : int or shape tuple, optional
        Axes over which to shift.  Default is None which shifts all axes.

    See Also
    --------
    ifftshift

    """
    tmp = asarray(x)
    ndim = len(tmp.shape)
    if axes is None:
        axes = range(ndim)
    y = tmp
    for k in axes:
        n = tmp.shape[k]
        p2 = (n+1)/2
        mylist = concatenate((arange(p2,n),arange(p2)))
        y = take(y,mylist,k)
    return y

def det(a):
    """Compute the determinant of a matrix

    Parameters
    ----------
    a : array-like, shape (M, M)

    Returns
    -------
    det : float or complex
        Determinant of a

    Notes
    -----
    The determinant is computed via LU factorization, LAPACK routine z/dgetrf.
    """
    a = asarray(a)
    _assertRank2(a)
    _assertSquareness(a)
    t, result_t = _commonType(a)
    a = _fastCopyAndTranspose(t, a)
    n = a.shape[0]
    if isComplexType(t):
        lapack_routine = lapack_lite.zgetrf
    else:
        lapack_routine = lapack_lite.dgetrf
    pivots = zeros((n,), fortran_int)
    results = lapack_routine(n, n, a, n, pivots, 0)
    info = results['info']
    if (info < 0):
        raise TypeError, "Illegal input to Fortran routine"
    elif (info > 0):
        return 0.0
    sign = add.reduce(pivots != arange(1, n+1)) % 2
    return (1.-2.*sign)*multiply.reduce(diagonal(a), axis=-1)

def ihfft(a, n=None, axis=-1):
    """
    Compute the inverse fft of a signal whose spectrum has Hermitian symmetry.

    Parameters
    ----------
    a : array_like
        Input array.
    n : int, optional
        Length of the ihfft.
    axis : int, optional
        Axis over which to compute the ihfft.

    See also
    --------
    rfft, hfft

    Notes
    -----
    These are a pair analogous to rfft/irfft, but for the
    opposite case: here the signal is real in the frequency domain and has
    Hermite symmetry in the time domain. So here it's hermite_fft for which
    you must supply the length of the result if it is to be odd.

    ihfft(hfft(a), len(a)) == a
    within numerical accuracy.

    """

    a = asarray(a).astype(float)
    if n is None:
        n = shape(a)[axis]
    return conjugate(rfft(a, n, axis))/n

def hfft(a, n=None, axis=-1):
    """
    Compute the fft of a signal which spectrum has Hermitian symmetry.

    Parameters
    ----------
    a : array
        input array
    n : int
        length of the hfft
    axis : int
        axis over which to compute the hfft

    See also
    --------
    rfft
    ihfft

    Notes
    -----
    These are a pair analogous to rfft/irfft, but for the
    opposite case: here the signal is real in the frequency domain and has
    Hermite symmetry in the time domain. So here it's hermite_fft for which
    you must supply the length of the result if it is to be odd.

    ihfft(hfft(a), len(a)) == a
    within numerical accuracy.

    """

    a = asarray(a).astype(complex)
    if n is None:
        n = (shape(a)[axis] - 1) * 2
    return irfft(conjugate(a), n, axis) * n

def irfftn(a, s=None, axes=None):
    """
    Compute the n-dimensional inverse fft of a real array.

    Parameters
    ----------
    a : array (real)
        input array
    s : sequence (int)
        shape of the inverse fft
    axis : int
        axis over which to compute the inverse fft

    Notes
    -----
    The transform implemented in ifftn is applied along
    all axes but the last, then the transform implemented in irfft is performed
    along the last axis. As with irfft, the length of the result along that
    axis must be specified if it is to be odd.

    """

    a = asarray(a).astype(complex)
    s, axes = _cook_nd_args(a, s, axes, invreal=1)
    for ii in range(len(axes) - 1):
        a = ifft(a, s[ii], axes[ii])
    a = irfft(a, s[-1], axes[-1])
    return a

def rfftn(a, s=None, axes=None):
    """
    Compute the n-dimensional fft of a real array.

    Parameters
    ----------
    a : array (real)
        input array
    s : sequence (int)
        shape of the fft
    axis : int
        axis over which to compute the fft

    Notes
    -----
    A real transform as rfft is performed along the axis specified by the last
    element of axes, then complex transforms as fft are performed along the
    other axes.

    """

    a = asarray(a).astype(float)
    s, axes = _cook_nd_args(a, s, axes)
    a = rfft(a, s[-1], axes[-1])
    for ii in range(len(axes) - 1):
        a = fft(a, s[ii], axes[ii])
    return a

def rfft(a, n=None, axis=-1):
    """
    Compute the one-dimensional fft for real input.

    Return the n point discrete Fourier transform of the real valued
    array a. n defaults to the length of a. n is the length of the
    input, not the output.

    Parameters
    ----------

    a : array
        input array with real data type
    n : int
        length of the fft
    axis : int
        axis over which to compute the fft

    Notes
    -----

    The returned array will be the nonnegative frequency terms of the
    Hermite-symmetric, complex transform of the real array. So for an 8-point
    transform, the frequencies in the result are [ 0, 1, 2, 3, 4]. The first
    term will be real, as will the last if n is even. The negative frequency
    terms are not needed because they are the complex conjugates of the
    positive frequency terms. (This is what I mean when I say
    Hermite-symmetric.)

    This is most efficient for n a power of two.

    """

    a = asarray(a).astype(float)
    return _raw_fft(a, n, axis, fftpack.rffti, fftpack.rfftf, _real_fft_cache)

def _raw_fft(a, n=None, axis=-1, init_function=fftpack.cffti,
             work_function=fftpack.cfftf, fft_cache = _fft_cache ):
    a = asarray(a)

    if n == None: n = a.shape[axis]

    if n < 1: raise ValueError("Invalid number of FFT data points (%d) specified." % n)

    try:
        wsave = fft_cache[n]
    except(KeyError):
        wsave = init_function(n)
        fft_cache[n] = wsave

    if a.shape[axis] != n:
        s = list(a.shape)
        if s[axis] > n:
            index = [slice(None)]*len(s)
            index[axis] = slice(0,n)
            a = a[index]
        else:
            index = [slice(None)]*len(s)
            index[axis] = slice(0,s[axis])
            s[axis] = n
            z = zeros(s, a.dtype.char)
            z[index] = a
            a = z

    if axis != -1:
        a = swapaxes(a, axis, -1)
    r = work_function(a, wsave)
    if axis != -1:
        r = swapaxes(r, axis, -1)
    return r

def tensorinv(a, ind=2):
    """
    Find the 'inverse' of a N-d array

    The result is an inverse corresponding to the operation
    tensordot(a, b, ind), ie.,

        x == tensordot(tensordot(tensorinv(a), a, ind), x, ind)
          == tensordot(tensordot(a, tensorinv(a), ind), x, ind)

    for all x (up to floating-point accuracy).

    Parameters
    ----------
    a : array_like
        Tensor to 'invert'. Its shape must 'square', ie.,
        prod(a.shape[:ind]) == prod(a.shape[ind:])
    ind : integer > 0
        How many of the first indices are involved in the inverse sum.

    Returns
    -------
    b : array, shape a.shape[:ind]+a.shape[ind:]

    Raises LinAlgError if a is singular or not square

    Examples
    --------
    >>> a = np.eye(4*6)
    >>> a.shape = (4,6,8,3)
    >>> ainv = np.linalg.tensorinv(a, ind=2)
    >>> ainv.shape
    (8, 3, 4, 6)
    >>> b = np.random.randn(4,6)
    >>> np.allclose(np.tensordot(ainv, b), np.linalg.tensorsolve(a, b))
    True

    >>> a = np.eye(4*6)
    >>> a.shape = (24,8,3)
    >>> ainv = np.linalg.tensorinv(a, ind=1)
    >>> ainv.shape
    (8, 3, 24)
    >>> b = np.random.randn(24)
    >>> np.allclose(np.tensordot(ainv, b, 1), np.linalg.tensorsolve(a, b))
    True

    """
    a = asarray(a)
    oldshape = a.shape
    prod = 1
    if ind > 0:
        invshape = oldshape[ind:] + oldshape[:ind]
        for k in oldshape[ind:]:
            prod *= k
    else:
        raise ValueError, "Invalid ind argument."
    a = a.reshape(prod, -1)
    ia = inv(a)
    return ia.reshape(*invshape)

def _raw_fftnd(a, s=None, axes=None, function=fft):
    a = asarray(a)
    s, axes = _cook_nd_args(a, s, axes)
    itl = range(len(axes))
    itl.reverse()
    for ii in itl:
        a = function(a, n=s[ii], axis=axes[ii])
    return a

def fftshift(x, axes=None):
    """
    Shift the zero-frequency component to the center of the spectrum.

    This function swaps half-spaces for all axes listed (defaults to all).
    Note that ``y[0]`` is the Nyquist component only if ``len(x)`` is even.

    Parameters
    ----------
    x : array_like
        Input array.
    axes : int or shape tuple, optional
        Axes over which to shift.  Default is None, which shifts all axes.

    Returns
    -------
    y : ndarray
        The shifted array.

    See Also
    --------
    ifftshift : The inverse of `fftshift`.

    Examples
    --------
    >>> freqs = np.fft.fftfreq(10, 0.1)
    >>> freqs
    array([ 0.,  1.,  2.,  3.,  4., -5., -4., -3., -2., -1.])
    >>> np.fft.fftshift(freqs)
    array([-5., -4., -3., -2., -1.,  0.,  1.,  2.,  3.,  4.])

    Shift the zero-frequency component only along the second axis:

    >>> freqs = np.fft.fftfreq(9, d=1./9).reshape(3, 3)
    >>> freqs
    array([[ 0.,  1.,  2.],
           [ 3.,  4., -4.],
           [-3., -2., -1.]])
    >>> np.fft.fftshift(freqs, axes=(1,))
    array([[ 2.,  0.,  1.],
           [-4.,  3.,  4.],
           [-1., -3., -2.]])

    """
    tmp = asarray(x)
    ndim = len(tmp.shape)
    if axes is None:
        axes = list(range(ndim))
    elif isinstance(axes, (int, nt.integer)):
        axes = (axes,)
    y = tmp
    for k in axes:
        n = tmp.shape[k]
        p2 = (n+1)//2
        mylist = concatenate((arange(p2,n),arange(p2)))
        y = take(y,mylist,k)
    return y

def fftshift(x, axes=None):
    """
    Shift the zero-frequency component to the center of the spectrum.

    This function swaps half-spaces for all axes listed (defaults to all).
    Note that ``y[0]`` is the Nyquist component only if ``len(x)`` is even.

    Parameters
    ----------
    x : array_like
        Input array.
    axes : int or shape tuple, optional
        Axes over which to shift.  Default is None, which shifts all axes.

    Returns
    -------
    y : ndarray
        The shifted array.

    See Also
    --------
    ifftshift : The inverse of `fftshift`.

    Examples
    --------
    >>> freqs = np.fft.fftfreq(10, 0.1)
    >>> freqs
    array([ 0.,  1.,  2., ..., -3., -2., -1.])
    >>> np.fft.fftshift(freqs)
    array([-5., -4., -3., -2., -1.,  0.,  1.,  2.,  3.,  4.])

    Shift the zero-frequency component only along the second axis:

    >>> freqs = np.fft.fftfreq(9, d=1./9).reshape(3, 3)
    >>> freqs
    array([[ 0.,  1.,  2.],
           [ 3.,  4., -4.],
           [-3., -2., -1.]])
    >>> np.fft.fftshift(freqs, axes=(1,))
    array([[ 2.,  0.,  1.],
           [-4.,  3.,  4.],
           [-1., -3., -2.]])

    """
    x = asarray(x)
    if axes is None:
        axes = tuple(range(x.ndim))
        shift = [dim // 2 for dim in x.shape]
    elif isinstance(axes, integer_types):
        shift = x.shape[axes] // 2
    else:
        shift = [x.shape[ax] // 2 for ax in axes]

    return roll(x, shift, axes)

def ihfft(a, n=None, axis=-1):
    """
    Compute the inverse FFT of a signal which has Hermitian symmetry.

    Parameters
    ----------
    a : array_like
        Input array.
    n : int, optional
        Length of the inverse FFT.
        Number of points along transformation axis in the input to use.
        If `n` is smaller than the length of the input, the input is cropped.
        If it is larger, the input is padded with zeros. If `n` is not given,
        the length of the input along the axis specified by `axis` is used.
    axis : int, optional
        Axis over which to compute the inverse FFT. If not given, the last
        axis is used.

    Returns
    -------
    out : complex ndarray
        The truncated or zero-padded input, transformed along the axis
        indicated by `axis`, or the last one if `axis` is not specified.
        If `n` is even, the length of the transformed axis is ``(n/2)+1``.
        If `n` is odd, the length is ``(n+1)/2``.

    See also
    --------
    hfft, irfft

    Notes
    -----
    `hfft`/`ihfft` are a pair analogous to `rfft`/`irfft`, but for the
    opposite case: here the signal has Hermitian symmetry in the time domain
    and is real in the frequency domain. So here it's `hfft` for which
    you must supply the length of the result if it is to be odd:
    ``ihfft(hfft(a), len(a)) == a``, within numerical accuracy.

    Examples
    --------
    >>> spectrum = np.array([ 15, -4, 0, -1, 0, -4])
    >>> np.fft.ifft(spectrum)
    array([ 1.+0.j,  2.-0.j,  3.+0.j,  4.+0.j,  3.+0.j,  2.-0.j])
    >>> np.fft.ihfft(spectrum)
    array([ 1.-0.j,  2.-0.j,  3.-0.j,  4.-0.j])

    """

    a = asarray(a).astype(float)
    if n is None:
        n = shape(a)[axis]
    return conjugate(rfft(a, n, axis))/n