mtazzari / galario

Gpu Accelerated Library for Analysing Radio Interferometer Observations
https://mtazzari.github.io/galario/
GNU Lesser General Public License v3.0
31 stars 15 forks source link

Optimize rotix core #7

Closed mtazzari closed 7 years ago

mtazzari commented 7 years ago

The rotix_core function can be highly simplified by noting that:

    indu[i] = index + (u[i] - umin - index*du)/du =
            = index + (u[i] - umin)/du - index = (u[i] - umin)/du

And the same applies for indv. I think we can even delete the rotix kernels and put it directly in interpolate_core. Check if the indu and indv arrays are needed only in the interpolate. If yes, then it's a go for deleting rotix.

#ifdef __CUDACC__
__host__ __device__ inline void rotix_core
#else
inline void rotix_core
#endif
        (int const i, int const nx, dreal const umin, dreal const du, dreal const* const u, dreal const* const v, dreal* const __restrict__ indu, dreal*  const __restrict__ indv) {
    // u
    int index = floor((u[i] - umin) / du);
    dreal center_ind = umin + index * du;
    indu[i] = index + (u[i] - center_ind) / du;

    // v
    index = floor((v[i] - umin) / du);
    center_ind = umin + index * du;
    indv[i] = index + (v[i] - center_ind) / du;
}
fredRos commented 7 years ago

Turns out to be irrelevant when profiling

fredRos commented 7 years ago

Through git bisect I saw that commit a3e4003 breaks several unit tests related to rotix, for example test_loss. @mtazzari Could you look into fixing that? Now I'm not sure other changes I make can be tested properly.

mtazzari commented 7 years ago

Sorry, I don't get what you are trying to do. commit a3e4003 just changed the way du was computed: see line 551 in the commit page. I don't understand what rotix tests were broken and how this is an issue now. My latest commit have all the tests passing. To make things faster, we shall talk.

fredRos commented 7 years ago

on my machine with the latest master. The errors are the same as with a3e4003

python/py.test.sh -s python_package/tests/test_galario.py --gpu=0

============================= test session starts ==============================
platform linux -- Python 3.6.1, pytest-3.0.7, py-1.4.33, pluggy-0.4.0
rootdir: /home/beaujean/workspace/protoplanetary/galario/build3, inifile:
collected 33 items

python_package/tests/test_galario.py FF................FFsize:256, minuv:113.91199493408203, maxuv:19333.672265625
.size:2048, minuv:17.164953054780607, maxuv:19563.900024738796
.size:256, minuv:113.91199493408203, maxuv:19333.672265625
.size:2048, minuv:17.164953054780607, maxuv:19563.900024738796
.size:256, minuv:113.91199493408203, maxuv:19333.672265625
.size:2048, minuv:17.164953054780607, maxuv:19563.900024738796
.size:2048, minuv:17.164953231811523, maxuv:19563.899999999998
.size:2048, minuv:17.164953054780607, maxuv:19563.900024738796
.size:2048, minuv:17.164953231811523, maxuv:19563.899999999998
.size:2048, minuv:17.164953054780607, maxuv:19563.900024738796
.size:2048, minuv:17.164953231811523, maxuv:19563.899999999998
.size:2048, minuv:17.164953054780607, maxuv:19563.900024738796
.size:4096, minuv:17.164953054780607, maxuv:19563.900024738796
.

=================================== FAILURES ===================================
________________________________ test_rotix[SP] ________________________________

x_id = array([ True,  True,  True,  True,  True,  True,  True,  True,  True,  True], dtype=bool)
y_id = array([False, False, False, False, False, False, False, False, False, False], dtype=bool)
hasval = '+inf'

    def chk_same_position(x_id, y_id, hasval='nan'):
        """Handling nan/inf: check that x and y have the nan/inf at the same
            locations."""
        try:
>           assert_array_equal(x_id, y_id)

../../../../.local/miniconda3/envs/galario3/lib/python3.6/site-packages/numpy/testing/utils.py:729: 
_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ 

x = array([ True,  True,  True,  True,  True,  True,  True,  True,  True,  True], dtype=bool)
y = array([False, False, False, False, False, False, False, False, False, False], dtype=bool)
err_msg = '', verbose = True

    def assert_array_equal(x, y, err_msg='', verbose=True):
        """
        Raises an AssertionError if two array_like objects are not equal.

        Given two array_like objects, check that the shape is equal and all
        elements of these objects are equal. An exception is raised at
        shape mismatch or conflicting values. In contrast to the standard usage
        in numpy, NaNs are compared like numbers, no assertion is raised if
        both objects have NaNs in the same positions.

        The usual caution for verifying equality with floating point numbers is
        advised.

        Parameters
        ----------
        x : array_like
            The actual object to check.
        y : array_like
            The desired, expected object.
        err_msg : str, optional
            The error message to be printed in case of failure.
        verbose : bool, optional
            If True, the conflicting values are appended to the error message.

        Raises
        ------
        AssertionError
            If actual and desired objects are not equal.

        See Also
        --------
        assert_allclose: Compare two array_like objects for equality with desired
                         relative and/or absolute precision.
        assert_array_almost_equal_nulp, assert_array_max_ulp, assert_equal

        Examples
        --------
        The first assert does not raise an exception:

        >>> np.testing.assert_array_equal([1.0,2.33333,np.nan],
        ...                               [np.exp(0),2.33333, np.nan])

        Assert fails with numerical inprecision with floats:

        >>> np.testing.assert_array_equal([1.0,np.pi,np.nan],
        ...                               [1, np.sqrt(np.pi)**2, np.nan])
        ...
        <type 'exceptions.ValueError'>:
        AssertionError:
        Arrays are not equal
        <BLANKLINE>
        (mismatch 50.0%)
         x: array([ 1.        ,  3.14159265,         NaN])
         y: array([ 1.        ,  3.14159265,         NaN])

        Use `assert_allclose` or one of the nulp (number of floating point values)
        functions for these cases instead:

        >>> np.testing.assert_allclose([1.0,np.pi,np.nan],
        ...                            [1, np.sqrt(np.pi)**2, np.nan],
        ...                            rtol=1e-10, atol=0)

        """
        assert_array_compare(operator.__eq__, x, y, err_msg=err_msg,
>                            verbose=verbose, header='Arrays are not equal')

../../../../.local/miniconda3/envs/galario3/lib/python3.6/site-packages/numpy/testing/utils.py:871: 
_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ 

comparison = <built-in function eq>
x = array([ True,  True,  True,  True,  True,  True,  True,  True,  True,  True], dtype=bool)
y = array([False, False, False, False, False, False, False, False, False, False], dtype=bool)
err_msg = '', verbose = True, header = 'Arrays are not equal', precision = 6
equal_nan = True

    def assert_array_compare(comparison, x, y, err_msg='', verbose=True,
                             header='', precision=6, equal_nan=True):
        __tracebackhide__ = True  # Hide traceback for py.test
        from numpy.core import array, isnan, isinf, any, all, inf
        x = array(x, copy=False, subok=True)
        y = array(y, copy=False, subok=True)

        def safe_comparison(*args, **kwargs):
            # There are a number of cases where comparing two arrays hits special
            # cases in array_richcompare, specifically around strings and void
            # dtypes. Basically, we just can't do comparisons involving these
            # types, unless both arrays have exactly the *same* type. So
            # e.g. you can apply == to two string arrays, or two arrays with
            # identical structured dtypes. But if you compare a non-string array
            # to a string array, or two arrays with non-identical structured
            # dtypes, or anything like that, then internally stuff blows up.
            # Currently, when things blow up, we just return a scalar False or
            # True. But we also emit a DeprecationWarning, b/c eventually we
            # should raise an error here. (Ideally we might even make this work
            # properly, but since that will require rewriting a bunch of how
            # ufuncs work then we are not counting on that.)
            #
            # The point of this little function is to let the DeprecationWarning
            # pass (or maybe eventually catch the errors and return False, I
            # dunno, that's a little trickier and we can figure that out when the
            # time comes).
            # Note: Put a lock around warning filter (comment at lock definition)
            with _assert_comparison_lock, suppress_warnings() as sup:
                sup.filter(DeprecationWarning, ".*==")
                sup.filter(FutureWarning, ".*==")
                return comparison(*args, **kwargs)

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

        def chk_same_position(x_id, y_id, hasval='nan'):
            """Handling nan/inf: check that x and y have the nan/inf at the same
            locations."""
            try:
                assert_array_equal(x_id, y_id)
            except AssertionError:
                msg = build_err_msg([x, y],
                                    err_msg + '\nx and y %s location mismatch:'
                                    % (hasval), verbose=verbose, header=header,
                                    names=('x', 'y'), precision=precision)
                raise AssertionError(msg)

        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'), precision=precision)
                if not cond:
                    raise AssertionError(msg)

            if isnumber(x) and isnumber(y):
                if equal_nan:
                    x_isnan, y_isnan = isnan(x), isnan(y)
                    # Validate that NaNs are in the same place
                    if any(x_isnan) or any(y_isnan):
                        chk_same_position(x_isnan, y_isnan, hasval='nan')

                x_isinf, y_isinf = isinf(x), isinf(y)

                # Validate that infinite values are in the same place
                if any(x_isinf) or any(y_isinf):
                    # Check +inf and -inf separately, since they are different
                    chk_same_position(x == +inf, y == +inf, hasval='+inf')
                    chk_same_position(x == -inf, y == -inf, hasval='-inf')

                # Combine all the special values
                x_id, y_id = x_isinf, y_isinf
                if equal_nan:
                    x_id |= x_isnan
                    y_id |= y_isnan

                # Only do the comparison if actual values are left
                if all(x_id):
                    return

                if any(x_id):
                    val = safe_comparison(x[~x_id], y[~y_id])
                else:
                    val = safe_comparison(x, y)
            else:
                val = safe_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'), precision=precision)
                if not cond:
>                   raise AssertionError(msg)
E                   AssertionError: 
E                   Arrays are not equal
E                   
E                   (mismatch 100.0%)
E                    x: array([ True,  True,  True,  True,  True,  True,  True,  True,  True,  True], dtype=bool)
E                    y: array([False, False, False, False, False, False, False, False, False, False], dtype=bool)

../../../../.local/miniconda3/envs/galario3/lib/python3.6/site-packages/numpy/testing/utils.py:796: AssertionError

During handling of the above exception, another exception occurred:

size = 1024, real_type = 'float32', tol = 0.0001
acc_lib = <module 'galario.single' from '/home/beaujean/workspace/protoplanetary/galario/build3/python/galario/single/__init__.py'>

    @pytest.mark.parametrize("size, real_type, tol, acc_lib",
                             [(1024, 'float32', 1.e-4, g_single),
                              (1024, 'float64', 1.e-13, g_double)],
                             ids=["SP", "DP"])
    def test_rotix(size, real_type, tol, acc_lib):
        nsamples = 10
        maxuv = 1000.

        uv = pixel_coordinates(maxuv, size).astype(real_type)
        udat, vdat = create_sampling_points(nsamples, maxuv/4.8)
        assert len(udat) == nsamples
        assert len(vdat) == nsamples
        udat = udat.astype(real_type)
        vdat = vdat.astype(real_type)

        ui, vi = get_rotix_n(uv, uv, udat, vdat, len(uv))
        ui = ui.astype(real_type)
        vi = vi.astype(real_type)

        ui1, vi1 = acc_lib.acc_rotix(uv, udat, vdat)

>       np.testing.assert_allclose(ui1, ui, rtol=tol)

python_package/tests/test_galario.py:286: 
_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ 

x_id = array([ True,  True,  True,  True,  True,  True,  True,  True,  True,  True], dtype=bool)
y_id = array([False, False, False, False, False, False, False, False, False, False], dtype=bool)
hasval = '+inf'

    def chk_same_position(x_id, y_id, hasval='nan'):
        """Handling nan/inf: check that x and y have the nan/inf at the same
            locations."""
        try:
            assert_array_equal(x_id, y_id)
        except AssertionError:
            msg = build_err_msg([x, y],
                                err_msg + '\nx and y %s location mismatch:'
                                % (hasval), verbose=verbose, header=header,
                                names=('x', 'y'), precision=precision)
>           raise AssertionError(msg)
E           AssertionError: 
E           Not equal to tolerance rtol=0.0001, atol=0
E           
E           x and y +inf location mismatch:
E            x: array([ inf,  inf,  inf,  inf,  inf,  inf,  inf,  inf,  inf,  inf], dtype=float32)
E            y: array([ 458.470459,  610.98407 ,  365.234619,  323.449005,  555.142395,
E                   307.449402,  653.842163,  376.2453  ,  428.476685,  482.963196], dtype=float32)

../../../../.local/miniconda3/envs/galario3/lib/python3.6/site-packages/numpy/testing/utils.py:735: AssertionError
________________________________ test_rotix[DP] ________________________________

x_id = array([ True,  True,  True,  True,  True,  True,  True,  True,  True,  True], dtype=bool)
y_id = array([False, False, False, False, False, False, False, False, False, False], dtype=bool)
hasval = '+inf'

    def chk_same_position(x_id, y_id, hasval='nan'):
        """Handling nan/inf: check that x and y have the nan/inf at the same
            locations."""
        try:
>           assert_array_equal(x_id, y_id)

../../../../.local/miniconda3/envs/galario3/lib/python3.6/site-packages/numpy/testing/utils.py:729: 
_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ 

x = array([ True,  True,  True,  True,  True,  True,  True,  True,  True,  True], dtype=bool)
y = array([False, False, False, False, False, False, False, False, False, False], dtype=bool)
err_msg = '', verbose = True

    def assert_array_equal(x, y, err_msg='', verbose=True):
        """
        Raises an AssertionError if two array_like objects are not equal.

        Given two array_like objects, check that the shape is equal and all
        elements of these objects are equal. An exception is raised at
        shape mismatch or conflicting values. In contrast to the standard usage
        in numpy, NaNs are compared like numbers, no assertion is raised if
        both objects have NaNs in the same positions.

        The usual caution for verifying equality with floating point numbers is
        advised.

        Parameters
        ----------
        x : array_like
            The actual object to check.
        y : array_like
            The desired, expected object.
        err_msg : str, optional
            The error message to be printed in case of failure.
        verbose : bool, optional
            If True, the conflicting values are appended to the error message.

        Raises
        ------
        AssertionError
            If actual and desired objects are not equal.

        See Also
        --------
        assert_allclose: Compare two array_like objects for equality with desired
                         relative and/or absolute precision.
        assert_array_almost_equal_nulp, assert_array_max_ulp, assert_equal

        Examples
        --------
        The first assert does not raise an exception:

        >>> np.testing.assert_array_equal([1.0,2.33333,np.nan],
        ...                               [np.exp(0),2.33333, np.nan])

        Assert fails with numerical inprecision with floats:

        >>> np.testing.assert_array_equal([1.0,np.pi,np.nan],
        ...                               [1, np.sqrt(np.pi)**2, np.nan])
        ...
        <type 'exceptions.ValueError'>:
        AssertionError:
        Arrays are not equal
        <BLANKLINE>
        (mismatch 50.0%)
         x: array([ 1.        ,  3.14159265,         NaN])
         y: array([ 1.        ,  3.14159265,         NaN])

        Use `assert_allclose` or one of the nulp (number of floating point values)
        functions for these cases instead:

        >>> np.testing.assert_allclose([1.0,np.pi,np.nan],
        ...                            [1, np.sqrt(np.pi)**2, np.nan],
        ...                            rtol=1e-10, atol=0)

        """
        assert_array_compare(operator.__eq__, x, y, err_msg=err_msg,
>                            verbose=verbose, header='Arrays are not equal')

../../../../.local/miniconda3/envs/galario3/lib/python3.6/site-packages/numpy/testing/utils.py:871: 
_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ 

comparison = <built-in function eq>
x = array([ True,  True,  True,  True,  True,  True,  True,  True,  True,  True], dtype=bool)
y = array([False, False, False, False, False, False, False, False, False, False], dtype=bool)
err_msg = '', verbose = True, header = 'Arrays are not equal', precision = 6
equal_nan = True

    def assert_array_compare(comparison, x, y, err_msg='', verbose=True,
                             header='', precision=6, equal_nan=True):
        __tracebackhide__ = True  # Hide traceback for py.test
        from numpy.core import array, isnan, isinf, any, all, inf
        x = array(x, copy=False, subok=True)
        y = array(y, copy=False, subok=True)

        def safe_comparison(*args, **kwargs):
            # There are a number of cases where comparing two arrays hits special
            # cases in array_richcompare, specifically around strings and void
            # dtypes. Basically, we just can't do comparisons involving these
            # types, unless both arrays have exactly the *same* type. So
            # e.g. you can apply == to two string arrays, or two arrays with
            # identical structured dtypes. But if you compare a non-string array
            # to a string array, or two arrays with non-identical structured
            # dtypes, or anything like that, then internally stuff blows up.
            # Currently, when things blow up, we just return a scalar False or
            # True. But we also emit a DeprecationWarning, b/c eventually we
            # should raise an error here. (Ideally we might even make this work
            # properly, but since that will require rewriting a bunch of how
            # ufuncs work then we are not counting on that.)
            #
            # The point of this little function is to let the DeprecationWarning
            # pass (or maybe eventually catch the errors and return False, I
            # dunno, that's a little trickier and we can figure that out when the
            # time comes).
            # Note: Put a lock around warning filter (comment at lock definition)
            with _assert_comparison_lock, suppress_warnings() as sup:
                sup.filter(DeprecationWarning, ".*==")
                sup.filter(FutureWarning, ".*==")
                return comparison(*args, **kwargs)

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

        def chk_same_position(x_id, y_id, hasval='nan'):
            """Handling nan/inf: check that x and y have the nan/inf at the same
            locations."""
            try:
                assert_array_equal(x_id, y_id)
            except AssertionError:
                msg = build_err_msg([x, y],
                                    err_msg + '\nx and y %s location mismatch:'
                                    % (hasval), verbose=verbose, header=header,
                                    names=('x', 'y'), precision=precision)
                raise AssertionError(msg)

        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'), precision=precision)
                if not cond:
                    raise AssertionError(msg)

            if isnumber(x) and isnumber(y):
                if equal_nan:
                    x_isnan, y_isnan = isnan(x), isnan(y)
                    # Validate that NaNs are in the same place
                    if any(x_isnan) or any(y_isnan):
                        chk_same_position(x_isnan, y_isnan, hasval='nan')

                x_isinf, y_isinf = isinf(x), isinf(y)

                # Validate that infinite values are in the same place
                if any(x_isinf) or any(y_isinf):
                    # Check +inf and -inf separately, since they are different
                    chk_same_position(x == +inf, y == +inf, hasval='+inf')
                    chk_same_position(x == -inf, y == -inf, hasval='-inf')

                # Combine all the special values
                x_id, y_id = x_isinf, y_isinf
                if equal_nan:
                    x_id |= x_isnan
                    y_id |= y_isnan

                # Only do the comparison if actual values are left
                if all(x_id):
                    return

                if any(x_id):
                    val = safe_comparison(x[~x_id], y[~y_id])
                else:
                    val = safe_comparison(x, y)
            else:
                val = safe_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'), precision=precision)
                if not cond:
>                   raise AssertionError(msg)
E                   AssertionError: 
E                   Arrays are not equal
E                   
E                   (mismatch 100.0%)
E                    x: array([ True,  True,  True,  True,  True,  True,  True,  True,  True,  True], dtype=bool)
E                    y: array([False, False, False, False, False, False, False, False, False, False], dtype=bool)

../../../../.local/miniconda3/envs/galario3/lib/python3.6/site-packages/numpy/testing/utils.py:796: AssertionError

During handling of the above exception, another exception occurred:

size = 1024, real_type = 'float64', tol = 1e-13
acc_lib = <module 'galario.double' from '/home/beaujean/workspace/protoplanetary/galario/build3/python/galario/double/__init__.py'>

    @pytest.mark.parametrize("size, real_type, tol, acc_lib",
                             [(1024, 'float32', 1.e-4, g_single),
                              (1024, 'float64', 1.e-13, g_double)],
                             ids=["SP", "DP"])
    def test_rotix(size, real_type, tol, acc_lib):
        nsamples = 10
        maxuv = 1000.

        uv = pixel_coordinates(maxuv, size).astype(real_type)
        udat, vdat = create_sampling_points(nsamples, maxuv/4.8)
        assert len(udat) == nsamples
        assert len(vdat) == nsamples
        udat = udat.astype(real_type)
        vdat = vdat.astype(real_type)

        ui, vi = get_rotix_n(uv, uv, udat, vdat, len(uv))
        ui = ui.astype(real_type)
        vi = vi.astype(real_type)

        ui1, vi1 = acc_lib.acc_rotix(uv, udat, vdat)

>       np.testing.assert_allclose(ui1, ui, rtol=tol)

python_package/tests/test_galario.py:286: 
_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ 

x_id = array([ True,  True,  True,  True,  True,  True,  True,  True,  True,  True], dtype=bool)
y_id = array([False, False, False, False, False, False, False, False, False, False], dtype=bool)
hasval = '+inf'

    def chk_same_position(x_id, y_id, hasval='nan'):
        """Handling nan/inf: check that x and y have the nan/inf at the same
            locations."""
        try:
            assert_array_equal(x_id, y_id)
        except AssertionError:
            msg = build_err_msg([x, y],
                                err_msg + '\nx and y %s location mismatch:'
                                % (hasval), verbose=verbose, header=header,
                                names=('x', 'y'), precision=precision)
>           raise AssertionError(msg)
E           AssertionError: 
E           Not equal to tolerance rtol=1e-13, atol=0
E           
E           x and y +inf location mismatch:
E            x: array([ inf,  inf,  inf,  inf,  inf,  inf,  inf,  inf,  inf,  inf])
E            y: array([ 458.470451,  610.984082,  365.23462 ,  323.449008,  555.142405,
E                   307.449384,  653.842193,  376.245319,  428.47669 ,  482.963208])

../../../../.local/miniconda3/envs/galario3/lib/python3.6/site-packages/numpy/testing/utils.py:735: AssertionError
______________________________ test_loss[SP_par1] ______________________________

nsamples = 100, real_type = 'float32', complex_type = 'complex64', rtol = 0.0001
atol = 0.001
acc_lib = <module 'galario.single' from '/home/beaujean/workspace/protoplanetary/galario/build3/python/galario/single/__init__.py'>
pars = {'wle_m': 0.003, 'x0_arcsec': 0.4, 'y0_arcsec': 4.0}

    @pytest.mark.parametrize("nsamples, real_type, complex_type, rtol, atol, acc_lib, pars",
                             [(100, 'float32', 'complex64',  1e-4,  1e-3, g_single, par1),
                              (1000, 'float64', 'complex128', 1e-14, 1e-10, g_double, par1)],
                             ids=["SP_par1", "DP_par1"])
    def test_loss(nsamples, real_type, complex_type, rtol, atol, acc_lib, pars):
        # try to find out where precision is lost

        wle_m = pars.get('wle_m', 0.003)
        x0_arcsec = pars.get('x0_arcsec', 0.4)
        y0_arcsec = pars.get('y0_arcsec', 10.)

        # generate the samples
        maxuv_generator = 3.e3
        udat, vdat = create_sampling_points(nsamples, maxuv_generator, dtype=real_type)

        # compute the matrix size and maxuv
        size, minuv, maxuv = matrix_size(udat, vdat)
        uv = pixel_coordinates(maxuv, size)

        # create model complex image (it happens to have 0 imaginary part)
        reference_image = create_reference_image(size=size, kernel='gaussian', dtype=complex_type)
        ref_complex = reference_image.copy()

        ###
        # shift - FFT - shift
        ###
        cpu_shift_fft_shift = np.fft.fftshift(np.fft.fft2(np.fft.fftshift(reference_image)))

        shift_acc = reference_image.copy()
        acc_lib.fftshift_fft2d_fftshift(shift_acc)

        np.testing.assert_allclose(cpu_shift_fft_shift.real, shift_acc.real, rtol, atol)
        np.testing.assert_allclose(cpu_shift_fft_shift.imag, shift_acc.imag, rtol, atol)

        ###
        # phase
        ###
        factor = sec2rad/wle_m*maxuv
        fourier_shifted = Fourier_shift_static(cpu_shift_fft_shift, x0_arcsec, y0_arcsec, wle_m, maxuv)
        acc_lib.apply_phase_2d(shift_acc, x0_arcsec*factor, y0_arcsec*factor)

        # lose some absolute precision here  --> not anymore
        # atol *= 2
        np.testing.assert_allclose(fourier_shifted.real, shift_acc.real, rtol, atol)
        np.testing.assert_allclose(fourier_shifted.imag, shift_acc.imag, rtol, atol)
        # but continue with previous tolerance
        # atol /= 2

        ###
        # rotation indices
        ###
        uroti, vroti = get_rotix_n(uv, uv, udat, vdat, size)
        ui1, vi1 = acc_lib.acc_rotix(uv.astype(real_type), udat.astype(real_type), vdat.astype(real_type))

>       np.testing.assert_allclose(ui1, uroti, rtol, atol)
E       AssertionError: 
E       Not equal to tolerance rtol=0.0001, atol=0.001
E       
E       (mismatch 100.0%)
E        x: array([ 118.854355,  147.450668,  101.372635,   93.537834,  136.980347,
E                90.53791 ,  155.486557,  103.437141,  113.230522,  123.446747,
E               137.839371,  112.262718,  125.376747,  104.865051,  136.284317,...
E        y: array([ 118.0326  ,  146.431198,  100.671753,   92.891121,  136.033279,
E                89.911934,  154.41153 ,  102.721985,  112.447655,  122.593246,
E               136.886368,  111.486542,  124.509895,  104.140022,  135.342056,...

python_package/tests/test_galario.py:499: AssertionError
______________________________ test_loss[DP_par1] ______________________________

nsamples = 1000, real_type = 'float64', complex_type = 'complex128'
rtol = 1e-14, atol = 1e-10
acc_lib = <module 'galario.double' from '/home/beaujean/workspace/protoplanetary/galario/build3/python/galario/double/__init__.py'>
pars = {'wle_m': 0.003, 'x0_arcsec': 0.4, 'y0_arcsec': 4.0}

    @pytest.mark.parametrize("nsamples, real_type, complex_type, rtol, atol, acc_lib, pars",
                             [(100, 'float32', 'complex64',  1e-4,  1e-3, g_single, par1),
                              (1000, 'float64', 'complex128', 1e-14, 1e-10, g_double, par1)],
                             ids=["SP_par1", "DP_par1"])
    def test_loss(nsamples, real_type, complex_type, rtol, atol, acc_lib, pars):
        # try to find out where precision is lost

        wle_m = pars.get('wle_m', 0.003)
        x0_arcsec = pars.get('x0_arcsec', 0.4)
        y0_arcsec = pars.get('y0_arcsec', 10.)

        # generate the samples
        maxuv_generator = 3.e3
        udat, vdat = create_sampling_points(nsamples, maxuv_generator, dtype=real_type)

        # compute the matrix size and maxuv
        size, minuv, maxuv = matrix_size(udat, vdat)
        uv = pixel_coordinates(maxuv, size)

        # create model complex image (it happens to have 0 imaginary part)
        reference_image = create_reference_image(size=size, kernel='gaussian', dtype=complex_type)
        ref_complex = reference_image.copy()

        ###
        # shift - FFT - shift
        ###
        cpu_shift_fft_shift = np.fft.fftshift(np.fft.fft2(np.fft.fftshift(reference_image)))

        shift_acc = reference_image.copy()
        acc_lib.fftshift_fft2d_fftshift(shift_acc)

        np.testing.assert_allclose(cpu_shift_fft_shift.real, shift_acc.real, rtol, atol)
        np.testing.assert_allclose(cpu_shift_fft_shift.imag, shift_acc.imag, rtol, atol)

        ###
        # phase
        ###
        factor = sec2rad/wle_m*maxuv
        fourier_shifted = Fourier_shift_static(cpu_shift_fft_shift, x0_arcsec, y0_arcsec, wle_m, maxuv)
        acc_lib.apply_phase_2d(shift_acc, x0_arcsec*factor, y0_arcsec*factor)

        # lose some absolute precision here  --> not anymore
        # atol *= 2
        np.testing.assert_allclose(fourier_shifted.real, shift_acc.real, rtol, atol)
        np.testing.assert_allclose(fourier_shifted.imag, shift_acc.imag, rtol, atol)
        # but continue with previous tolerance
        # atol /= 2

        ###
        # rotation indices
        ###
        uroti, vroti = get_rotix_n(uv, uv, udat, vdat, size)
        ui1, vi1 = acc_lib.acc_rotix(uv.astype(real_type), udat.astype(real_type), vdat.astype(real_type))

>       np.testing.assert_allclose(ui1, uroti, rtol, atol)
E       AssertionError: 
E       Not equal to tolerance rtol=1e-14, atol=1e-10
E       
E       (mismatch 100.0%)
E        x: array([ 1003.243414,  1241.545963,   857.562428,   792.272409,
E               1154.293343,   767.272998,  1308.511762,   874.766646,
E                956.378163,  1041.513347,  1161.451931,   948.3131  ,...
E        y: array([  945.199198,  1169.71438 ,   807.946813,   746.43425 ,
E               1087.509896,   722.881219,  1232.805768,   824.155654,
E                901.045409,   981.254964,  1094.254314,   893.446963,...

python_package/tests/test_galario.py:499: AssertionError
===================== 4 failed, 29 passed in 19.32 seconds =====================
mtazzari commented 7 years ago

I have tested on Python 2 and 3 on my linux and all 33 tests pass.

The code was on the latest master commit bb8f3a2d36e41877c77786b32a43fe113e18b1f9 Author: Frederik Beaujean Frederik.Beaujean@lmu.de Date: Fri Jul 7 13:43:02 2017 +0200 [docs] Extract numpy docstring into html

Perhaps you need to cleanup your build directory or so..?

On 7 Jul 2017, at 12:57, Frederik Beaujean notifications@github.com wrote:

============================= test session starts ============================== platform linux -- Python 3.6.1, pytest-3.0.7, py-1.4.33, pluggy-0.4.0 rootdir: /home/beaujean/workspace/protoplanetary/galario/build3, inifile: collected 33 items

python_package/tests/test_galario.py FF................FFsize:256, minuv:113.91199493408203, maxuv:19333.672265625 .size:2048, minuv:17.164953054780607, maxuv:19563.900024738796 .size:256, minuv:113.91199493408203, maxuv:19333.672265625 .size:2048, minuv:17.164953054780607, maxuv:19563.900024738796 .size:256, minuv:113.91199493408203, maxuv:19333.672265625 .size:2048, minuv:17.164953054780607, maxuv:19563.900024738796 .size:2048, minuv:17.164953231811523, maxuv:19563.899999999998 .size:2048, minuv:17.164953054780607, maxuv:19563.900024738796 .size:2048, minuv:17.164953231811523, maxuv:19563.899999999998 .size:2048, minuv:17.164953054780607, maxuv:19563.900024738796 .size:2048, minuv:17.164953231811523, maxuv:19563.899999999998 .size:2048, minuv:17.164953054780607, maxuv:19563.900024738796 .size:4096, minuv:17.164953054780607, maxuv:19563.900024738796 .

=================================== FAILURES =================================== ____ test_rotix[SP] ____

x_id = array([ True, True, True, True, True, True, True, True, True, True], dtype=bool) y_id = array([False, False, False, False, False, False, False, False, False, False], dtype=bool) hasval = '+inf'

def chk_same_position(x_id, y_id, hasval='nan'):
    """Handling nan/inf: check that x and y have the nan/inf at the same
        locations."""
    try:
      assert_array_equal(x_id, y_id)

../../../../.local/miniconda3/envs/galario3/lib/python3.6/site-packages/numpy/testing/utils.py:729:


x = array([ True, True, True, True, True, True, True, True, True, True], dtype=bool) y = array([False, False, False, False, False, False, False, False, False, False], dtype=bool) err_msg = '', verbose = True

def assert_array_equal(x, y, err_msg='', verbose=True):
    """
    Raises an AssertionError if two array_like objects are not equal.

    Given two array_like objects, check that the shape is equal and all
    elements of these objects are equal. An exception is raised at
    shape mismatch or conflicting values. In contrast to the standard usage
    in numpy, NaNs are compared like numbers, no assertion is raised if
    both objects have NaNs in the same positions.

    The usual caution for verifying equality with floating point numbers is
    advised.

    Parameters
    ----------
    x : array_like
        The actual object to check.
    y : array_like
        The desired, expected object.
    err_msg : str, optional
        The error message to be printed in case of failure.
    verbose : bool, optional
        If True, the conflicting values are appended to the error message.

    Raises
    ------
    AssertionError
        If actual and desired objects are not equal.

    See Also
    --------
    assert_allclose: Compare two array_like objects for equality with desired
                     relative and/or absolute precision.
    assert_array_almost_equal_nulp, assert_array_max_ulp, assert_equal

    Examples
    --------
    The first assert does not raise an exception:

    >>> np.testing.assert_array_equal([1.0,2.33333,np.nan],
    ...                               [np.exp(0),2.33333, np.nan])

    Assert fails with numerical inprecision with floats:

    >>> np.testing.assert_array_equal([1.0,np.pi,np.nan],
    ...                               [1, np.sqrt(np.pi)**2, np.nan])
    ...
    <type 'exceptions.ValueError'>:
    AssertionError:
    Arrays are not equal
    <BLANKLINE>
    (mismatch 50.0%)
     x: array([ 1.        ,  3.14159265,         NaN])
     y: array([ 1.        ,  3.14159265,         NaN])

    Use `assert_allclose` or one of the nulp (number of floating point values)
    functions for these cases instead:

    >>> np.testing.assert_allclose([1.0,np.pi,np.nan],
    ...                            [1, np.sqrt(np.pi)**2, np.nan],
    ...                            rtol=1e-10, atol=0)

    """
    assert_array_compare(operator.__eq__, x, y, err_msg=err_msg,
                       verbose=verbose, header='Arrays are not equal')

../../../../.local/miniconda3/envs/galario3/lib/python3.6/site-packages/numpy/testing/utils.py:871:


comparison = x = array([ True, True, True, True, True, True, True, True, True, True], dtype=bool) y = array([False, False, False, False, False, False, False, False, False, False], dtype=bool) err_msg = '', verbose = True, header = 'Arrays are not equal', precision = 6 equal_nan = True

def assert_array_compare(comparison, x, y, err_msg='', verbose=True,
                         header='', precision=6, equal_nan=True):
    __tracebackhide__ = True  # Hide traceback for py.test
    from numpy.core import array, isnan, isinf, any, all, inf
    x = array(x, copy=False, subok=True)
    y = array(y, copy=False, subok=True)

    def safe_comparison(*args, **kwargs):
        # There are a number of cases where comparing two arrays hits special
        # cases in array_richcompare, specifically around strings and void
        # dtypes. Basically, we just can't do comparisons involving these
        # types, unless both arrays have exactly the *same* type. So
        # e.g. you can apply == to two string arrays, or two arrays with
        # identical structured dtypes. But if you compare a non-string array
        # to a string array, or two arrays with non-identical structured
        # dtypes, or anything like that, then internally stuff blows up.
        # Currently, when things blow up, we just return a scalar False or
        # True. But we also emit a DeprecationWarning, b/c eventually we
        # should raise an error here. (Ideally we might even make this work
        # properly, but since that will require rewriting a bunch of how
        # ufuncs work then we are not counting on that.)
        #
        # The point of this little function is to let the DeprecationWarning
        # pass (or maybe eventually catch the errors and return False, I
        # dunno, that's a little trickier and we can figure that out when the
        # time comes).
        # Note: Put a lock around warning filter (comment at lock definition)
        with _assert_comparison_lock, suppress_warnings() as sup:
            sup.filter(DeprecationWarning, ".*==")
            sup.filter(FutureWarning, ".*==")
            return comparison(*args, **kwargs)

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

    def chk_same_position(x_id, y_id, hasval='nan'):
        """Handling nan/inf: check that x and y have the nan/inf at the same
        locations."""
        try:
            assert_array_equal(x_id, y_id)
        except AssertionError:
            msg = build_err_msg([x, y],
                                err_msg + '\nx and y %s location mismatch:'
                                % (hasval), verbose=verbose, header=header,
                                names=('x', 'y'), precision=precision)
            raise AssertionError(msg)

    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'), precision=precision)
            if not cond:
                raise AssertionError(msg)

        if isnumber(x) and isnumber(y):
            if equal_nan:
                x_isnan, y_isnan = isnan(x), isnan(y)
                # Validate that NaNs are in the same place
                if any(x_isnan) or any(y_isnan):
                    chk_same_position(x_isnan, y_isnan, hasval='nan')

            x_isinf, y_isinf = isinf(x), isinf(y)

            # Validate that infinite values are in the same place
            if any(x_isinf) or any(y_isinf):
                # Check +inf and -inf separately, since they are different
                chk_same_position(x == +inf, y == +inf, hasval='+inf')
                chk_same_position(x == -inf, y == -inf, hasval='-inf')

            # Combine all the special values
            x_id, y_id = x_isinf, y_isinf
            if equal_nan:
                x_id |= x_isnan
                y_id |= y_isnan

            # Only do the comparison if actual values are left
            if all(x_id):
                return

            if any(x_id):
                val = safe_comparison(x[~x_id], y[~y_id])
            else:
                val = safe_comparison(x, y)
        else:
            val = safe_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'), precision=precision)
            if not cond:
              raise AssertionError(msg)

E AssertionError: E Arrays are not equal E
E (mismatch 100.0%) E x: array([ True, True, True, True, True, True, True, True, True, True], dtype=bool) E y: array([False, False, False, False, False, False, False, False, False, False], dtype=bool)

../../../../.local/miniconda3/envs/galario3/lib/python3.6/site-packages/numpy/testing/utils.py:796: AssertionError

During handling of the above exception, another exception occurred:

size = 1024, real_type = 'float32', tol = 0.0001 acc_lib = <module 'galario.single' from '/home/beaujean/workspace/protoplanetary/galario/build3/python/galario/single/init.py'>

@pytest.mark.parametrize("size, real_type, tol, acc_lib",
                         [(1024, 'float32', 1.e-4, g_single),
                          (1024, 'float64', 1.e-13, g_double)],
                         ids=["SP", "DP"])
def test_rotix(size, real_type, tol, acc_lib):
    nsamples = 10
    maxuv = 1000.

    uv = pixel_coordinates(maxuv, size).astype(real_type)
    udat, vdat = create_sampling_points(nsamples, maxuv/4.8)
    assert len(udat) == nsamples
    assert len(vdat) == nsamples
    udat = udat.astype(real_type)
    vdat = vdat.astype(real_type)

    ui, vi = get_rotix_n(uv, uv, udat, vdat, len(uv))
    ui = ui.astype(real_type)
    vi = vi.astype(real_type)

    ui1, vi1 = acc_lib.acc_rotix(uv, udat, vdat)
  np.testing.assert_allclose(ui1, ui, rtol=tol)

python_package/tests/test_galario.py:286:


x_id = array([ True, True, True, True, True, True, True, True, True, True], dtype=bool) y_id = array([False, False, False, False, False, False, False, False, False, False], dtype=bool) hasval = '+inf'

def chk_same_position(x_id, y_id, hasval='nan'):
    """Handling nan/inf: check that x and y have the nan/inf at the same
        locations."""
    try:
        assert_array_equal(x_id, y_id)
    except AssertionError:
        msg = build_err_msg([x, y],
                            err_msg + '\nx and y %s location mismatch:'
                            % (hasval), verbose=verbose, header=header,
                            names=('x', 'y'), precision=precision)
      raise AssertionError(msg)

E AssertionError: E Not equal to tolerance rtol=0.0001, atol=0 E
E x and y +inf location mismatch: E x: array([ inf, inf, inf, inf, inf, inf, inf, inf, inf, inf], dtype=float32) E y: array([ 458.470459, 610.98407 , 365.234619, 323.449005, 555.142395, E 307.449402, 653.842163, 376.2453 , 428.476685, 482.963196], dtype=float32)

../../../../.local/miniconda3/envs/galario3/lib/python3.6/site-packages/numpy/testing/utils.py:735: AssertionError ____ test_rotix[DP] ____

x_id = array([ True, True, True, True, True, True, True, True, True, True], dtype=bool) y_id = array([False, False, False, False, False, False, False, False, False, False], dtype=bool) hasval = '+inf'

def chk_same_position(x_id, y_id, hasval='nan'):
    """Handling nan/inf: check that x and y have the nan/inf at the same
        locations."""
    try:
      assert_array_equal(x_id, y_id)

../../../../.local/miniconda3/envs/galario3/lib/python3.6/site-packages/numpy/testing/utils.py:729:


x = array([ True, True, True, True, True, True, True, True, True, True], dtype=bool) y = array([False, False, False, False, False, False, False, False, False, False], dtype=bool) err_msg = '', verbose = True

def assert_array_equal(x, y, err_msg='', verbose=True):
    """
    Raises an AssertionError if two array_like objects are not equal.

    Given two array_like objects, check that the shape is equal and all
    elements of these objects are equal. An exception is raised at
    shape mismatch or conflicting values. In contrast to the standard usage
    in numpy, NaNs are compared like numbers, no assertion is raised if
    both objects have NaNs in the same positions.

    The usual caution for verifying equality with floating point numbers is
    advised.

    Parameters
    ----------
    x : array_like
        The actual object to check.
    y : array_like
        The desired, expected object.
    err_msg : str, optional
        The error message to be printed in case of failure.
    verbose : bool, optional
        If True, the conflicting values are appended to the error message.

    Raises
    ------
    AssertionError
        If actual and desired objects are not equal.

    See Also
    --------
    assert_allclose: Compare two array_like objects for equality with desired
                     relative and/or absolute precision.
    assert_array_almost_equal_nulp, assert_array_max_ulp, assert_equal

    Examples
    --------
    The first assert does not raise an exception:

    >>> np.testing.assert_array_equal([1.0,2.33333,np.nan],
    ...                               [np.exp(0),2.33333, np.nan])

    Assert fails with numerical inprecision with floats:

    >>> np.testing.assert_array_equal([1.0,np.pi,np.nan],
    ...                               [1, np.sqrt(np.pi)**2, np.nan])
    ...
    <type 'exceptions.ValueError'>:
    AssertionError:
    Arrays are not equal
    <BLANKLINE>
    (mismatch 50.0%)
     x: array([ 1.        ,  3.14159265,         NaN])
     y: array([ 1.        ,  3.14159265,         NaN])

    Use `assert_allclose` or one of the nulp (number of floating point values)
    functions for these cases instead:

    >>> np.testing.assert_allclose([1.0,np.pi,np.nan],
    ...                            [1, np.sqrt(np.pi)**2, np.nan],
    ...                            rtol=1e-10, atol=0)

    """
    assert_array_compare(operator.__eq__, x, y, err_msg=err_msg,
                       verbose=verbose, header='Arrays are not equal')

../../../../.local/miniconda3/envs/galario3/lib/python3.6/site-packages/numpy/testing/utils.py:871:


comparison = x = array([ True, True, True, True, True, True, True, True, True, True], dtype=bool) y = array([False, False, False, False, False, False, False, False, False, False], dtype=bool) err_msg = '', verbose = True, header = 'Arrays are not equal', precision = 6 equal_nan = True

def assert_array_compare(comparison, x, y, err_msg='', verbose=True,
                         header='', precision=6, equal_nan=True):
    __tracebackhide__ = True  # Hide traceback for py.test
    from numpy.core import array, isnan, isinf, any, all, inf
    x = array(x, copy=False, subok=True)
    y = array(y, copy=False, subok=True)

    def safe_comparison(*args, **kwargs):
        # There are a number of cases where comparing two arrays hits special
        # cases in array_richcompare, specifically around strings and void
        # dtypes. Basically, we just can't do comparisons involving these
        # types, unless both arrays have exactly the *same* type. So
        # e.g. you can apply == to two string arrays, or two arrays with
        # identical structured dtypes. But if you compare a non-string array
        # to a string array, or two arrays with non-identical structured
        # dtypes, or anything like that, then internally stuff blows up.
        # Currently, when things blow up, we just return a scalar False or
        # True. But we also emit a DeprecationWarning, b/c eventually we
        # should raise an error here. (Ideally we might even make this work
        # properly, but since that will require rewriting a bunch of how
        # ufuncs work then we are not counting on that.)
        #
        # The point of this little function is to let the DeprecationWarning
        # pass (or maybe eventually catch the errors and return False, I
        # dunno, that's a little trickier and we can figure that out when the
        # time comes).
        # Note: Put a lock around warning filter (comment at lock definition)
        with _assert_comparison_lock, suppress_warnings() as sup:
            sup.filter(DeprecationWarning, ".*==")
            sup.filter(FutureWarning, ".*==")
            return comparison(*args, **kwargs)

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

    def chk_same_position(x_id, y_id, hasval='nan'):
        """Handling nan/inf: check that x and y have the nan/inf at the same
        locations."""
        try:
            assert_array_equal(x_id, y_id)
        except AssertionError:
            msg = build_err_msg([x, y],
                                err_msg + '\nx and y %s location mismatch:'
                                % (hasval), verbose=verbose, header=header,
                                names=('x', 'y'), precision=precision)
            raise AssertionError(msg)

    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'), precision=precision)
            if not cond:
                raise AssertionError(msg)

        if isnumber(x) and isnumber(y):
            if equal_nan:
                x_isnan, y_isnan = isnan(x), isnan(y)
                # Validate that NaNs are in the same place
                if any(x_isnan) or any(y_isnan):
                    chk_same_position(x_isnan, y_isnan, hasval='nan')

            x_isinf, y_isinf = isinf(x), isinf(y)

            # Validate that infinite values are in the same place
            if any(x_isinf) or any(y_isinf):
                # Check +inf and -inf separately, since they are different
                chk_same_position(x == +inf, y == +inf, hasval='+inf')
                chk_same_position(x == -inf, y == -inf, hasval='-inf')

            # Combine all the special values
            x_id, y_id = x_isinf, y_isinf
            if equal_nan:
                x_id |= x_isnan
                y_id |= y_isnan

            # Only do the comparison if actual values are left
            if all(x_id):
                return

            if any(x_id):
                val = safe_comparison(x[~x_id], y[~y_id])
            else:
                val = safe_comparison(x, y)
        else:
            val = safe_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'), precision=precision)
            if not cond:
              raise AssertionError(msg)

E AssertionError: E Arrays are not equal E
E (mismatch 100.0%) E x: array([ True, True, True, True, True, True, True, True, True, True], dtype=bool) E y: array([False, False, False, False, False, False, False, False, False, False], dtype=bool)

../../../../.local/miniconda3/envs/galario3/lib/python3.6/site-packages/numpy/testing/utils.py:796: AssertionError

During handling of the above exception, another exception occurred:

size = 1024, real_type = 'float64', tol = 1e-13 acc_lib = <module 'galario.double' from '/home/beaujean/workspace/protoplanetary/galario/build3/python/galario/double/init.py'>

@pytest.mark.parametrize("size, real_type, tol, acc_lib",
                         [(1024, 'float32', 1.e-4, g_single),
                          (1024, 'float64', 1.e-13, g_double)],
                         ids=["SP", "DP"])
def test_rotix(size, real_type, tol, acc_lib):
    nsamples = 10
    maxuv = 1000.

    uv = pixel_coordinates(maxuv, size).astype(real_type)
    udat, vdat = create_sampling_points(nsamples, maxuv/4.8)
    assert len(udat) == nsamples
    assert len(vdat) == nsamples
    udat = udat.astype(real_type)
    vdat = vdat.astype(real_type)

    ui, vi = get_rotix_n(uv, uv, udat, vdat, len(uv))
    ui = ui.astype(real_type)
    vi = vi.astype(real_type)

    ui1, vi1 = acc_lib.acc_rotix(uv, udat, vdat)
  np.testing.assert_allclose(ui1, ui, rtol=tol)

python_package/tests/test_galario.py:286:


x_id = array([ True, True, True, True, True, True, True, True, True, True], dtype=bool) y_id = array([False, False, False, False, False, False, False, False, False, False], dtype=bool) hasval = '+inf'

def chk_same_position(x_id, y_id, hasval='nan'):
    """Handling nan/inf: check that x and y have the nan/inf at the same
        locations."""
    try:
        assert_array_equal(x_id, y_id)
    except AssertionError:
        msg = build_err_msg([x, y],
                            err_msg + '\nx and y %s location mismatch:'
                            % (hasval), verbose=verbose, header=header,
                            names=('x', 'y'), precision=precision)
      raise AssertionError(msg)

E AssertionError: E Not equal to tolerance rtol=1e-13, atol=0 E
E x and y +inf location mismatch: E x: array([ inf, inf, inf, inf, inf, inf, inf, inf, inf, inf]) E y: array([ 458.470451, 610.984082, 365.23462 , 323.449008, 555.142405, E 307.449384, 653.842193, 376.245319, 428.47669 , 482.963208])

../../../../.local/miniconda3/envs/galario3/lib/python3.6/site-packages/numpy/testing/utils.py:735: AssertionError __ test_loss[SP_par1] __

nsamples = 100, real_type = 'float32', complex_type = 'complex64', rtol = 0.0001 atol = 0.001 acc_lib = <module 'galario.single' from '/home/beaujean/workspace/protoplanetary/galario/build3/python/galario/single/init.py'> pars = {'wle_m': 0.003, 'x0_arcsec': 0.4, 'y0_arcsec': 4.0}

@pytest.mark.parametrize("nsamples, real_type, complex_type, rtol, atol, acc_lib, pars",
                         [(100, 'float32', 'complex64',  1e-4,  1e-3, g_single, par1),
                          (1000, 'float64', 'complex128', 1e-14, 1e-10, g_double, par1)],
                         ids=["SP_par1", "DP_par1"])
def test_loss(nsamples, real_type, complex_type, rtol, atol, acc_lib, pars):
    # try to find out where precision is lost

    wle_m = pars.get('wle_m', 0.003)
    x0_arcsec = pars.get('x0_arcsec', 0.4)
    y0_arcsec = pars.get('y0_arcsec', 10.)

    # generate the samples
    maxuv_generator = 3.e3
    udat, vdat = create_sampling_points(nsamples, maxuv_generator, dtype=real_type)

    # compute the matrix size and maxuv
    size, minuv, maxuv = matrix_size(udat, vdat)
    uv = pixel_coordinates(maxuv, size)

    # create model complex image (it happens to have 0 imaginary part)
    reference_image = create_reference_image(size=size, kernel='gaussian', dtype=complex_type)
    ref_complex = reference_image.copy()

    ###
    # shift - FFT - shift
    ###
    cpu_shift_fft_shift = np.fft.fftshift(np.fft.fft2(np.fft.fftshift(reference_image)))

    shift_acc = reference_image.copy()
    acc_lib.fftshift_fft2d_fftshift(shift_acc)

    np.testing.assert_allclose(cpu_shift_fft_shift.real, shift_acc.real, rtol, atol)
    np.testing.assert_allclose(cpu_shift_fft_shift.imag, shift_acc.imag, rtol, atol)

    ###
    # phase
    ###
    factor = sec2rad/wle_m*maxuv
    fourier_shifted = Fourier_shift_static(cpu_shift_fft_shift, x0_arcsec, y0_arcsec, wle_m, maxuv)
    acc_lib.apply_phase_2d(shift_acc, x0_arcsec*factor, y0_arcsec*factor)

    # lose some absolute precision here  --> not anymore
    # atol *= 2
    np.testing.assert_allclose(fourier_shifted.real, shift_acc.real, rtol, atol)
    np.testing.assert_allclose(fourier_shifted.imag, shift_acc.imag, rtol, atol)
    # but continue with previous tolerance
    # atol /= 2

    ###
    # rotation indices
    ###
    uroti, vroti = get_rotix_n(uv, uv, udat, vdat, size)
    ui1, vi1 = acc_lib.acc_rotix(uv.astype(real_type), udat.astype(real_type), vdat.astype(real_type))
  np.testing.assert_allclose(ui1, uroti, rtol, atol)

E AssertionError: E Not equal to tolerance rtol=0.0001, atol=0.001 E
E (mismatch 100.0%) E x: array([ 118.854355, 147.450668, 101.372635, 93.537834, 136.980347, E 90.53791 , 155.486557, 103.437141, 113.230522, 123.446747, E 137.839371, 112.262718, 125.376747, 104.865051, 136.284317,... E y: array([ 118.0326 , 146.431198, 100.671753, 92.891121, 136.033279, E 89.911934, 154.41153 , 102.721985, 112.447655, 122.593246, E 136.886368, 111.486542, 124.509895, 104.140022, 135.342056,...

python_package/tests/test_galario.py:499: AssertionError __ test_loss[DP_par1] __

nsamples = 1000, real_type = 'float64', complex_type = 'complex128' rtol = 1e-14, atol = 1e-10 acc_lib = <module 'galario.double' from '/home/beaujean/workspace/protoplanetary/galario/build3/python/galario/double/init.py'> pars = {'wle_m': 0.003, 'x0_arcsec': 0.4, 'y0_arcsec': 4.0}

@pytest.mark.parametrize("nsamples, real_type, complex_type, rtol, atol, acc_lib, pars",
                         [(100, 'float32', 'complex64',  1e-4,  1e-3, g_single, par1),
                          (1000, 'float64', 'complex128', 1e-14, 1e-10, g_double, par1)],
                         ids=["SP_par1", "DP_par1"])
def test_loss(nsamples, real_type, complex_type, rtol, atol, acc_lib, pars):
    # try to find out where precision is lost

    wle_m = pars.get('wle_m', 0.003)
    x0_arcsec = pars.get('x0_arcsec', 0.4)
    y0_arcsec = pars.get('y0_arcsec', 10.)

    # generate the samples
    maxuv_generator = 3.e3
    udat, vdat = create_sampling_points(nsamples, maxuv_generator, dtype=real_type)

    # compute the matrix size and maxuv
    size, minuv, maxuv = matrix_size(udat, vdat)
    uv = pixel_coordinates(maxuv, size)

    # create model complex image (it happens to have 0 imaginary part)
    reference_image = create_reference_image(size=size, kernel='gaussian', dtype=complex_type)
    ref_complex = reference_image.copy()

    ###
    # shift - FFT - shift
    ###
    cpu_shift_fft_shift = np.fft.fftshift(np.fft.fft2(np.fft.fftshift(reference_image)))

    shift_acc = reference_image.copy()
    acc_lib.fftshift_fft2d_fftshift(shift_acc)

    np.testing.assert_allclose(cpu_shift_fft_shift.real, shift_acc.real, rtol, atol)
    np.testing.assert_allclose(cpu_shift_fft_shift.imag, shift_acc.imag, rtol, atol)

    ###
    # phase
    ###
    factor = sec2rad/wle_m*maxuv
    fourier_shifted = Fourier_shift_static(cpu_shift_fft_shift, x0_arcsec, y0_arcsec, wle_m, maxuv)
    acc_lib.apply_phase_2d(shift_acc, x0_arcsec*factor, y0_arcsec*factor)

    # lose some absolute precision here  --> not anymore
    # atol *= 2
    np.testing.assert_allclose(fourier_shifted.real, shift_acc.real, rtol, atol)
    np.testing.assert_allclose(fourier_shifted.imag, shift_acc.imag, rtol, atol)
    # but continue with previous tolerance
    # atol /= 2

    ###
    # rotation indices
    ###
    uroti, vroti = get_rotix_n(uv, uv, udat, vdat, size)
    ui1, vi1 = acc_lib.acc_rotix(uv.astype(real_type), udat.astype(real_type), vdat.astype(real_type))
  np.testing.assert_allclose(ui1, uroti, rtol, atol)

E AssertionError: E Not equal to tolerance rtol=1e-14, atol=1e-10 E
E (mismatch 100.0%) E x: array([ 1003.243414, 1241.545963, 857.562428, 792.272409, E 1154.293343, 767.272998, 1308.511762, 874.766646, E 956.378163, 1041.513347, 1161.451931, 948.3131 ,... E y: array([ 945.199198, 1169.71438 , 807.946813, 746.43425 , E 1087.509896, 722.881219, 1232.805768, 824.155654, E 901.045409, 981.254964, 1094.254314, 893.446963,...

python_package/tests/test_galario.py:499: AssertionError ===================== 4 failed, 29 passed in 19.32 seconds =====================

— You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub https://github.com/mtazzari/galario/issues/7#issuecomment-313662103, or mute the thread https://github.com/notifications/unsubscribe-auth/AFvcgkpKo3twRjsWEx15PO1P_M94xw2nks5sLh0wgaJpZM4NjwHC.

mtazzari commented 7 years ago

It does not show up in the profiling for the moment, so we lower the priority for future releases.

fredRos commented 7 years ago

no more indices needed as of 1485b07