Checking for zero-volume cells: Fixes issue #79

[unoebauer]

• don't split parent cell if v_inner_boundary equal to an existing cell boundary
• check if same error type may also affect outer boundary
• insert final check for zero volume cells
• add testing scenarios

[unoebauer]

Very strange things are happening: stumbled across a small bug in test_ascii_reader.py while running the test scripts. Commit acd72ea fixes this. This erroneous line has been in there since commit 13b0d2d from 07.02.2105 and it has never triggered a failure during the Travis-Ci process. @wkerzendorf, any idea why this is the case?

Offline, the testing failed (rightfully so) with the following error:

_______________________________________________________________ test_ascii_reader_density_boundaries[v_inner_boundary3-v_outer_boundary3-actual_v_inner3-actual_v_outer3-1-16] ________________________________________________________________

v_inner_boundary = <Quantity 1600.0 km / s>, v_outer_boundary = <Quantity 6000.0 km / s>, actual_v_inner = <Quantity 1600.0 km / s>, actual_v_outer = <Quantity 6000.0 km / s>, inner_index = 1, outer_index = 16

    @pytest.mark.parametrize("v_inner_boundary, v_outer_boundary, actual_v_inner, actual_v_outer, inner_index, outer_index", [
        (0.0 * u.km/u.s, np.inf * u.km/u.s, np.nan, np.nan, None, None),
        (500 * u.km/u.s, 6000 * u.km/u.s, np.nan, 6000 * u.km/u.s, None, 16),
        (1300 * u.km/u.s, 6000 * u.km/u.s, 1300 * u.km/u.s, 6000 * u.km/u.s, 0, 16),
        (1600 * u.km/u.s, 6000 * u.km/u.s, 1600 * u.km/u.s, 6000 * u.km/u.s, 1, 16)
    ])
    def test_ascii_reader_density_boundaries(v_inner_boundary, v_outer_boundary, actual_v_inner, actual_v_outer,
                                             inner_index, outer_index):
        v_inner, v_outer, mean_densities, inner_boundary_index, outer_boundary_index = \
            io.read_density_file(data_path('artis_model.dat'), 'artis', 19 * u.day, v_inner_boundary, v_outer_boundary)

        assert inner_boundary_index == inner_index
        assert outer_boundary_index == outer_index

        if not np.isnan(actual_v_inner):
            npt.assert_allclose(v_inner[0].value,
                                actual_v_inner.to(v_inner[0].unit).value)

        if not np.isnan(actual_v_outer):
            npt.assert_allclose(v_outer[-1],
>                               actual_v_outer.to(v_outer[-1].unit).value)

tardis/io/tests/test_ascii_readers.py:59: 
_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _
/home/ulrich/python-virtualenv/tardis/local/lib/python2.7/site-packages/numpy/testing/utils.py:1347: in assert_allclose
    verbose=verbose, header=header)
/home/ulrich/python-virtualenv/tardis/local/lib/python2.7/site-packages/numpy/testing/utils.py:689: in assert_array_compare
    val = safe_comparison(x, y)
/home/ulrich/python-virtualenv/tardis/local/lib/python2.7/site-packages/numpy/testing/utils.py:636: in safe_comparison
    return comparison(*args, **kwargs)
/home/ulrich/python-virtualenv/tardis/local/lib/python2.7/site-packages/numpy/testing/utils.py:1342: in compare
    equal_nan=equal_nan)
/home/ulrich/python-virtualenv/tardis/local/lib/python2.7/site-packages/numpy/core/numeric.py:2354: in isclose
    return within_tol(x, y, atol, rtol)
/home/ulrich/python-virtualenv/tardis/local/lib/python2.7/site-packages/numpy/core/numeric.py:2337: in within_tol
    result = less_equal(abs(x-y), atol + rtol * abs(y))
_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _

self = <Quantity [  6.00000000e+08] cm / s>, obj = array([ 1.]), context = (<ufunc 'subtract'>, (<Quantity [  6.00000000e+08] cm / s>, array([  6.00000000e+08])), 0)

    def __array_prepare__(self, obj, context=None):
        # This method gets called by Numpy whenever a ufunc is called on the
        # array. The object passed in ``obj`` is an empty version of the
        # output array which we can e.g. change to an array sub-class, add
        # attributes to, etc. After this is called, then the ufunc is called
        # and the values in this empty array are set.

        # If no context is set, just return the input
        if context is None:
            return obj

        # Find out which ufunc is being used
        function = context[0]

        from .quantity_helper import UNSUPPORTED_UFUNCS, UFUNC_HELPERS

        # Check whether we even support this ufunc
        if function in UNSUPPORTED_UFUNCS:
            raise TypeError("Cannot use function '{0}' with quantities"
                            .format(function.__name__))

        # Now find out what arguments were passed to the ufunc, usually, this
        # will include at least the present object, and another, which could
        # be a Quantity, or a Numpy array, etc. when using two-argument ufuncs.
        args = context[1][:function.nin]
        units = [getattr(arg, 'unit', None) for arg in args]

        # If the ufunc is supported, then we call a helper function (defined
        # in quantity_helper.py) which returns the scale by which the inputs
        # should be multiplied before being passed to the ufunc, as well as
        # the unit the output from the ufunc will have.
        if function in UFUNC_HELPERS:
            converters, result_unit = UFUNC_HELPERS[function](function, *units)
        else:
            raise TypeError("Unknown ufunc {0}.  Please raise issue on "
                            "https://github.com/astropy/astropy"
                            .format(function.__name__))

        if any(converter is False for converter in converters):
            # for two-argument ufuncs with a quantity and a non-quantity,
            # the quantity normally needs to be dimensionless, *except*
            # if the non-quantity can have arbitrary unit, i.e., when it
            # is all zero, infinity or NaN.  In that case, the non-quantity
            # can just have the unit of the quantity
            # (this allows, e.g., `q > 0.` independent of unit)
            maybe_arbitrary_arg = args[converters.index(False)]
            try:
                if _can_have_arbitrary_unit(maybe_arbitrary_arg):
                    converters = [None, None]
                else:
                    raise UnitsError("Can only apply '{0}' function to "
                                     "dimensionless quantities when other "
                                     "argument is not a quantity (unless the "
                                     "latter is all zero/infinity/nan)"
>                                    .format(function.__name__))
E                                    UnitsError: Can only apply 'subtract' function to dimensionless quantities when other argument is not a quantity (unless the latter is all zero/infinity/nan)

/home/ulrich/python-virtualenv/tardis/local/lib/python2.7/site-packages/astropy/units/quantity.py:317: UnitsError

[unoebauer]

Travis seems to be happy. Could you, @ssim or @wkerzendorf, have a look at this and merge it?

[ssim]

Looks like a good fix to me - merging!