def plot_cub_as_curve(c, colors=None, plot_kwargs=None, legend_prefix='',
                      show_axis_labels=True, show_legend=False, axes=None):
    """
    Plot a cuboid (ndims <= 2) as curve(s).
    If the input is 1D: one single curve.
    If the input is 2D:
       * multiple curves are plotted: one for each domain value on the 1st axis.
       * legends are shown to display which domain value is associated
         to which curve.

    Args:
        - colors (dict <domain value>: <matplotlib color>):
            associate domain values of the 1st axis to color curves
        - plot_kwargs (dict <arg name>:<arg value>):
            dictionary of named argument passed to the plot function
        - legend_prefix (str): prefix to prepend to legend labels.

    Return:
        None
    """
    import matplotlib.pyplot as plt

    def protect_latex_str(s):
        return s.replace('_','\_')
        
    axes = axes or plt.gca()
    colors = colors or {}
    plot_kwargs = plot_kwargs or {}
    if c.get_ndims() == 1:
        dom = c.axes_domains[c.axes_names[0]]
        if np.issubsctype(dom.dtype, str):
            dom = np.arange(len(dom))
        axes.plot(dom, c.data, **plot_kwargs)
        if np.issubsctype(c.axes_domains[c.axes_names[0]], str):
            set_int_tick_labels(axes.xaxis, c.axes_domains[c.axes_names[0]],
                                rotation=30)
    elif c.get_ndims() == 2:
        for val, sub_c in c.split(c.axes_names[0]).iteritems():
            pkwargs = plot_kwargs.copy()
            col = colors.get(val, None)
            if col is not None:
                pkwargs['color'] = col
            pkwargs['label'] = protect_latex_str(legend_prefix + \
                                                 c.axes_names[0] + \
                                                 '=' + str(val))
            plot_cub_as_curve(sub_c, plot_kwargs=pkwargs, axes=axes,
                              show_axis_labels=False)
        if show_legend:
            axes.legend()

    else:
        raise Exception('xndarray has too many dims (%d), expected at most 2' \
                        %c.get_ndims())

    if show_axis_labels:
        if c.get_ndims() == 1:
            axes.set_xlabel(protect_latex_str(c.axes_names[0]))
        else:
            axes.set_xlabel(protect_latex_str(c.axes_names[1]))
        axes.set_ylabel(protect_latex_str(c.value_label))

def test_float32_input():
    # Check that an float32 input is correctly output as float32
    Yl, Yh = dtwavexfm2(lena.astype(np.float32))
    assert np.issubsctype(Yl.dtype, np.float32)
    assert np.all(list(np.issubsctype(x.dtype, np.complex64) for x in Yh))

    lena_recon = dtwaveifm2(Yl, Yh)
    assert np.issubsctype(lena_recon.dtype, np.float32)

def test_float32_input_inv():
    # Check that an float32 input is correctly output as float32
    Yl, Yh = dtwavexfm(np.array([1, 2, 3, 4]).astype(np.float32))
    assert np.issubsctype(Yl.dtype, np.float32)
    assert np.all(list(np.issubsctype(x.dtype, np.complex64) for x in Yh))

    recon = dtwaveifm(Yl, Yh)
    assert np.issubsctype(recon.dtype, np.float32)

def test_float32_recon():
    # Check that an float32 input is correctly output as float32
    Yl, Yh = dtwavexfm3(ellipsoid.astype(np.float32))
    assert np.issubsctype(Yl.dtype, np.float32)
    assert np.all(list(np.issubsctype(x.dtype, np.complex64) for x in Yh))

    recon = dtwaveifm3(Yl, Yh)
    assert np.issubsctype(recon.dtype, np.float32)

def is_range_domain(d):
    #print 'd:', d.dtype
    if not (np.issubsctype(d.dtype, np.unicode) or \
            np.issubsctype(d.dtype, np.str) or (d.size==1)):
        delta = np.diff(d)
        return (delta == delta[0]).all()
    else:
        return False

def _encode_sql_value(val):
    if np.issubsctype(type(val), np.str_):
        return sqlite3.Binary(val)
    elif np.issubsctype(type(val), np.bool_):
        return bool(val)
    elif np.issubsctype(type(val), np.integer):
        return int(val)
    elif np.issubsctype(type(val), np.floating):
        return float(val)
    else:
        return val

def _encode_sql_value(val):
    # This is a non-trivial speedup for bulk inserts (e.g. creating thousands
    # of events while loading a file):
    if type(val) in ok_types:
        return val
    if np.issubsctype(type(val), np.str_):
        return sqlite3.Binary(val)
    elif np.issubsctype(type(val), np.bool_):
        return bool(val)
    elif np.issubsctype(type(val), np.integer):
        return int(val)
    elif np.issubsctype(type(val), np.floating):
        return float(val)
    else:
        return val

def convert_value_read(value):
    """Convert attribute value from bytes to string."""
    if isinstance(value, bytes):
        return value.decode()
    elif not np.isscalar(value) and np.issubsctype(value, np.bytes_):
        return value.astype(np.str_)
    return value

def isnan(a):
    """
    isnan is equivalent to np.isnan, except that it returns False instead of
    raising a TypeError if the argument is an array of non-numeric.
    """
    if isinstance(a, np.ndarray):
        return np.issubsctype(a, np.floating) and np.isnan(a)
    else:
        return np.isnan(a)

def appropriate_complex_type_for(X):
    """Return an appropriate complex data type depending on the type of X. If X
    is already complex, return that, if it is floating point return a complex
    type of the appropriate size and if it is integer, choose an complex
    floating point type depending on the result of :py:func:`numpy.asfarray`.

    """
    X = asfarray(X)
    
    if np.issubsctype(X.dtype, np.complex64) or np.issubsctype(X.dtype, np.complex128):
        return X.dtype
    elif np.issubsctype(X.dtype, np.float32):
        return np.complex64
    elif np.issubsctype(X.dtype, np.float64):
        return np.complex128

    # God knows, err on the side of caution
    return np.complex128

def _copy_array_if_base_present(a):
    """
    Copies the array if its base points to a parent array.
    """
    if a.base is not None:
        return a.copy()
    elif np.issubsctype(a, np.float32):
        return np.array(a, dtype=np.double)
    else:
        return a

def test_ufuncs():
    # Cannot use fixture due to bug in pytest
    fn = Rn(3)

    for name, n_args, n_out, _ in UFUNCS:
        if (np.issubsctype(fn.dtype, np.floating) and
            name in ['bitwise_and', 'bitwise_or', 'bitwise_xor', 'invert',
                     'left_shift', 'right_shift']):
            # Skip integer only methods if floating point type
            continue
        yield _impl_test_ufuncs, fn, name, n_args, n_out

    def _update_colors(self, numcolors=None):
        """ Update the colors cache using our color mapper and based
        on our number of levels.  The **mode** parameter accounts for fenceposting:
          - If **mode** is "poly", then the number of colors to generate is 1
            less than the number of levels
          - If **mode** is "line", then the number of colors to generate is
            equal to the number of levels
        """
        if numcolors is None:
            numcolors = len(self._levels)

        colors = self.colors
        # If we are given no colors, set a default for all levels
        if colors is None:
            self._color_map_trait = "black"
            self._colors = [self._color_map_trait_] * numcolors

        # If we are given a single color, apply it to all levels
        elif isinstance(colors, basestring):
            self._color_map_trait = colors
            self._colors = [self._color_map_trait_] * numcolors

        # If we are given a colormap, use it to map all the levels to colors
        elif isinstance(colors, ColorMapper):
            self._colors =  []
            mapped_colors = self.color_mapper.map_screen(array(self._levels))
            for i in range(numcolors):
                self._color_map_trait = tuple(mapped_colors[i])
                self._colors.append(self._color_map_trait_)

        # A list or tuple
        # This could be a length 3 or 4 sequence of scalars, which indicates
        # a color; otherwise, this is interpreted as a list of items to
        # be converted via self._color_map_trait.
        else:
            if len(colors) in (3,4) and \
                    (isscalar(colors[0]) and issubsctype(type(colors[0]), number)):
                self._color_map_trait = colors
                self._colors = [self._color_map_trait_] * numcolors
            else:
                # if the list of colors is shorter than the list of levels, simply
                # repeat colors from the beginning of the list as needed
                self._colors = []
                for i in range(len(self._levels)):
                    self._color_map_trait = colors[i%len(colors)]
                    self._colors.append(self._color_map_trait_)

        self._colors_cache_valid = True
        return

def test_ufuncs(fn, ufunc):
    name, n_args, n_out, _ = ufunc
    if (np.issubsctype(fn.dtype, np.floating) and
            name in ['bitwise_and',
                     'bitwise_or',
                     'bitwise_xor',
                     'invert',
                     'left_shift',
                     'right_shift']):
        # Skip integer only methods if floating point type
        return

    # Get the ufunc from numpy as reference
    ufunc = getattr(np, name)

    # Create some data
    arrays, vectors = example_vectors(fn, n_args + n_out)
    in_arrays = arrays[:n_args]
    out_arrays = arrays[n_args:]
    data_vector = vectors[0]
    in_vectors = vectors[1:n_args]
    out_vectors = vectors[n_args:]

    # Out of place:
    np_result = ufunc(*in_arrays)
    vec_fun = getattr(data_vector.ufunc, name)
    odl_result = vec_fun(*in_vectors)
    assert all_almost_equal(np_result, odl_result)

    # Test type of output
    if n_out == 1:
        assert isinstance(odl_result, fn.element_type)
    elif n_out > 1:
        for i in range(n_out):
            assert isinstance(odl_result[i], fn.element_type)

    # In place:
    np_result = ufunc(*(in_arrays + out_arrays))
    vec_fun = getattr(data_vector.ufunc, name)
    odl_result = vec_fun(*(in_vectors + out_vectors))
    assert all_almost_equal(np_result, odl_result)

    # Test inplace actually holds:
    if n_out == 1:
        assert odl_result is out_vectors[0]
    elif n_out > 1:
        for i in range(n_out):
            assert odl_result[i] is out_vectors[i]

    def slice(self, columns_specifier):
        """Locate a subset of design matrix columns, specified symbolically.

        A patsy design matrix has two levels of structure: the individual
        columns (which are named), and the :ref:`terms <formulas>` in
        the formula that generated those columns. This is a one-to-many
        relationship: a single term may span several columns. This method
        provides a user-friendly API for locating those columns.

        (While we talk about columns here, this is probably most useful for
        indexing into other arrays that are derived from the design matrix,
        such as regression coefficients or covariance matrices.)

        The `columns_specifier` argument can take a number of forms:

        * A term name
        * A column name
        * A :class:`Term` object
        * An integer giving a raw index
        * A raw slice object

        In all cases, a Python :func:`slice` object is returned, which can be
        used directly for indexing.

        Example::

          y, X = dmatrices("y ~ a", demo_data("y", "a", nlevels=3))
          betas = np.linalg.lstsq(X, y)[0]
          a_betas = betas[X.design_info.slice("a")]

        (If you want to look up a single individual column by name, use
        ``design_info.column_name_indexes[name]``.)
        """
        if isinstance(columns_specifier, slice):
            return columns_specifier
        if np.issubsctype(type(columns_specifier), np.integer):
            return slice(columns_specifier, columns_specifier + 1)
        if (self.term_slices is not None
            and columns_specifier in self.term_slices):
            return self.term_slices[columns_specifier]
        if columns_specifier in self.term_name_slices:
            return self.term_name_slices[columns_specifier]
        if columns_specifier in self.column_name_indexes:
            idx = self.column_name_indexes[columns_specifier]
            return slice(idx, idx + 1)
        raise PatsyError("unknown column specified '%s'"
                            % (columns_specifier,))

def test_ufuncs():
    # Cannot use fixture due to bug in pytest
    spaces = [odl.uniform_discr([0, 0], [1, 1], [2, 2])]

    if odl.CUDA_AVAILABLE:
        spaces += [odl.uniform_discr([0, 0], [1, 1], [2, 2], impl='cuda')]

    for fn in spaces:
        for name, n_args, n_out, _ in odl.util.ufuncs.UFUNCS:
            if (np.issubsctype(fn.dtype, np.floating) and
                    name in ['bitwise_and',
                             'bitwise_or',
                             'bitwise_xor',
                             'invert',
                             'left_shift',
                             'right_shift']):
                # Skip integer only methods if floating point type
                continue
            yield _impl_test_ufuncs, fn, name, n_args, n_out

def mask_source_lonlats(source_def, mask):
    """Mask source longitudes and latitudes to match data mask."""
    source_geo_def = source_def

    # the data may have additional masked pixels
    # let's compare them to see if we can use the same area
    # assume lons and lats mask are the same
    if mask is not None and mask is not False and isinstance(source_geo_def, SwathDefinition):
        if np.issubsctype(mask.dtype, np.bool):
            # copy the source area and use it for the rest of the calculations
            LOG.debug("Copying source area to mask invalid dataset points")
            if mask.ndim != source_geo_def.lons.ndim:
                raise ValueError("Can't mask area, mask has different number "
                                 "of dimensions.")

            return SwathDefinition(source_geo_def.lons.where(~mask),
                                   source_geo_def.lats.where(~mask))
        else:
            return SwathDefinition(source_geo_def.lons.where(~xu.isnan(mask)),
                                   source_geo_def.lats.where(~xu.isnan(mask)))

    return source_geo_def

def is_real_floating_dtype(dtype):
    """`True` if ``dtype`` is real floating-point, else `False`."""
    return np.issubsctype(dtype, np.floating)

def is_complex_floating_dtype(dtype):
    """`True` if ``dtype`` is complex floating-point, else `False`."""
    return np.issubsctype(dtype, np.complexfloating)

def is_scalar_dtype(dtype):
    """`True` if ``dtype`` is scalar, else `False`."""
    return np.issubsctype(dtype, np.number)