Merge pull request #6166 from pv/sld

Better memory overlap detection

doc/release/1.11.0-notes.rst

@@ -29,6 +29,9 @@ as the argument to 'bins' results in the corresponding estimator being used.
 has been added, converting the previous vbench-based one. You can run the suite locally
 via ``python runtests.py --bench``. For more details, see ``benchmarks/README.rst``.
 
+* A new function ``np.shares_memory`` that can check exactly whether two
+arrays have memory overlap is added. ``np.may_share_memory`` also now
+has an option to spend more effort to reduce false positives.
 
 Improvements
 ============

doc/source/reference/routines.other.rst

@@ -22,3 +22,12 @@ Performance tuning
    restoredot
    setbufsize
    getbufsize
+
+Memory ranges
+-------------
+
+.. autosummary::
+   :toctree: generated/
+
+   shares_memory
+   may_share_memory

numpy/add_newdocs.py

@@ -3784,24 +3784,40 @@ def luf(lamdaexpr, *args, **kwargs):
     """))
 
 
-add_newdoc('numpy.core.multiarray', 'may_share_memory',
+add_newdoc('numpy.core.multiarray', 'shares_memory',
     """
-    Determine if two arrays can share memory
+    shares_memory(a, b, max_work=None)
 
-    The memory-bounds of a and b are computed.  If they overlap then
-    this function returns True.  Otherwise, it returns False.
-
-    A return of True does not necessarily mean that the two arrays
-    share any element.  It just means that they *might*.
+    Determine if two arrays share memory
 
     Parameters
     ----------
     a, b : ndarray
+        Input arrays
+    max_work : int, optional
+        Effort to spend on solving the overlap problem (maximum number
+        of candidate solutions to consider). The following special
+        values are recognized:
+
+        max_work=MAY_SHARE_EXACT  (default)
+            The problem is solved exactly. In this case, the function returns
+            True only if there is an element shared between the arrays.
+        max_work=MAY_SHARE_BOUNDS
+            Only the memory bounds of a and b are checked.
+
+    Raises
+    ------
+    numpy.TooHardError
+        Exceeded max_work.
 
     Returns
     -------
     out : bool
 
+    See Also
+    --------
+    may_share_memory
+
     Examples
     --------
     >>> np.may_share_memory(np.array([1,2]), np.array([5,8,9]))

numpy/core/_internal.py

@@ -759,3 +759,7 @@ def _gcd(a, b):
     while b:
         a, b = b, a % b
     return a
+
+# Exception used in shares_memory()
+class TooHardError(RuntimeError):
+    pass

numpy/core/bento.info

@@ -22,7 +22,8 @@ Library:
             src/multiarray/multiarraymodule_onefile.c
     Extension: multiarray_tests
         Sources:
-            src/multiarray/multiarray_tests.c.src
+            src/multiarray/multiarray_tests.c.src,
+            src/private/mem_overlap.c
     Extension: umath
         Sources:
             src/umath/umathmodule_onefile.c

numpy/core/bscript

@@ -484,6 +484,7 @@ def pre_build(context):
                 pjoin('src', 'multiarray', 'usertypes.c'),
                 pjoin('src', 'multiarray', 'vdot.c'),
                 pjoin('src', 'private', 'templ_common.h.src'),
+                pjoin('src', 'private', 'mem_overlap.c'),
                 ]
 
             if bld.env.HAS_CBLAS:

numpy/core/function_base.py

@@ -1,9 +1,9 @@
 from __future__ import division, absolute_import, print_function
 
-__all__ = ['logspace', 'linspace']
+__all__ = ['logspace', 'linspace', 'may_share_memory']
 
 from . import numeric as _nx
-from .numeric import result_type, NaN
+from .numeric import result_type, NaN, shares_memory, MAY_SHARE_BOUNDS, TooHardError
 
 
 def linspace(start, stop, num=50, endpoint=True, retstep=False, dtype=None):
@@ -201,3 +201,46 @@ def logspace(start, stop, num=50, endpoint=True, base=10.0, dtype=None):
     if dtype is None:
         return _nx.power(base, y)
     return _nx.power(base, y).astype(dtype)
+
+
+def may_share_memory(a, b, max_work=None):
+    """Determine if two arrays can share memory
+
+    A return of True does not necessarily mean that the two arrays
+    share any element.  It just means that they *might*.
+
+    Only the memory bounds of a and b are checked by default.
+
+    Parameters
+    ----------
+    a, b : ndarray
+        Input arrays
+    max_work : int, optional
+        Effort to spend on solving the overlap problem.  See
+        `shares_memory` for details.  Default for ``may_share_memory``
+        is to do a bounds check.
+
+    Returns
+    -------
+    out : bool
+
+    See Also
+    --------
+    shares_memory
+
+    Examples
+    --------
+    >>> np.may_share_memory(np.array([1,2]), np.array([5,8,9]))
+    False
+    >>> x = np.zeros([3, 4])
+    >>> np.may_share_memory(x[:,0], x[:,1])
+    True
+
+    """
+    if max_work is None:
+        max_work = MAY_SHARE_BOUNDS
+    try:
+        return shares_memory(a, b, max_work=max_work)
+    except (TooHardError, OverflowError):
+        # Unable to determine, assume yes
+        return True

numpy/core/numeric.py

@@ -10,6 +10,7 @@
                     ERR_DEFAULT, PINF, NAN)
 from . import numerictypes
 from .numerictypes import longlong, intc, int_, float_, complex_, bool_
+from ._internal import TooHardError
 
 if sys.version_info[0] >= 3:
     import pickle
@@ -39,8 +40,8 @@
     'getbufsize', 'seterrcall', 'geterrcall', 'errstate', 'flatnonzero',
     'Inf', 'inf', 'infty', 'Infinity', 'nan', 'NaN', 'False_', 'True_',
     'bitwise_not', 'CLIP', 'RAISE', 'WRAP', 'MAXDIMS', 'BUFSIZE',
-    'ALLOW_THREADS', 'ComplexWarning', 'may_share_memory', 'full',
-    'full_like', 'matmul',
+    'ALLOW_THREADS', 'ComplexWarning', 'full', 'full_like', 'matmul',
+    'shares_memory', 'MAY_SHARE_BOUNDS', 'MAY_SHARE_EXACT', 'TooHardError',
     ]
 
 if sys.version_info[0] < 3:
@@ -65,6 +66,8 @@ class ComplexWarning(RuntimeWarning):
 MAXDIMS = multiarray.MAXDIMS
 ALLOW_THREADS = multiarray.ALLOW_THREADS
 BUFSIZE = multiarray.BUFSIZE
+MAY_SHARE_BOUNDS = multiarray.MAY_SHARE_BOUNDS
+MAY_SHARE_EXACT = multiarray.MAY_SHARE_EXACT
 
 ndarray = multiarray.ndarray
 flatiter = multiarray.flatiter
@@ -375,7 +378,7 @@ def extend_all(module):
 fromiter = multiarray.fromiter
 fromfile = multiarray.fromfile
 frombuffer = multiarray.frombuffer
-may_share_memory = multiarray.may_share_memory
+shares_memory = multiarray.shares_memory
 if sys.version_info[0] < 3:
     newbuffer = multiarray.newbuffer
     getbuffer = multiarray.getbuffer

numpy/core/setup.py

@@ -762,6 +762,8 @@ def generate_multiarray_templated_sources(ext, build_dir):
             join('src', 'private', 'npy_config.h'),
             join('src', 'private', 'templ_common.h.src'),
             join('src', 'private', 'lowlevel_strided_loops.h'),
+            join('src', 'private', 'mem_overlap.h'),
+            join('src', 'private', 'npy_extint128.h'),
             join('include', 'numpy', 'arrayobject.h'),
             join('include', 'numpy', '_neighborhood_iterator_imp.h'),
             join('include', 'numpy', 'npy_endian.h'),
@@ -831,6 +833,7 @@ def generate_multiarray_templated_sources(ext, build_dir):
             join('src', 'multiarray', 'ucsnarrow.c'),
             join('src', 'multiarray', 'vdot.c'),
             join('src', 'private', 'templ_common.h.src'),
+            join('src', 'private', 'mem_overlap.c'),
             ]
 
     blas_info = get_info('blas_opt', 0)
@@ -959,7 +962,10 @@ def generate_umath_c(ext, build_dir):
     #######################################################################
 
     config.add_extension('multiarray_tests',
-                    sources=[join('src', 'multiarray', 'multiarray_tests.c.src')])
+                    sources=[join('src', 'multiarray', 'multiarray_tests.c.src'),
+                             join('src', 'private', 'mem_overlap.c')],
+                    depends=[join('src', 'private', 'mem_overlap.h'),
+                             join('src', 'private', 'npy_extint128.h')])
 
     #######################################################################
     #                        operand_flag_tests module                    #

numpy/core/src/multiarray/array_assign.c

@@ -23,6 +23,7 @@
 #include "array_assign.h"
 #include "common.h"
 #include "lowlevel_strided_loops.h"
+#include "mem_overlap.h"
 
 /* See array_assign.h for parameter documentation */
 NPY_NO_EXPORT int
@@ -101,27 +102,17 @@ raw_array_is_aligned(int ndim, char *data, npy_intp *strides, int alignment)
 }
 
 
-/* Gets a half-open range [start, end) which contains the array data */
-NPY_NO_EXPORT void
-get_array_memory_extents(PyArrayObject *arr,
-                         npy_uintp *out_start, npy_uintp *out_end)
-{
-    npy_intp low, upper;
-    offset_bounds_from_strides(PyArray_ITEMSIZE(arr), PyArray_NDIM(arr),
-                               PyArray_DIMS(arr), PyArray_STRIDES(arr),
-                               &low, &upper);
-    *out_start = (npy_uintp)PyArray_DATA(arr) + (npy_uintp)low;
-    *out_end = (npy_uintp)PyArray_DATA(arr) + (npy_uintp)upper;
-}
-
 /* Returns 1 if the arrays have overlapping data, 0 otherwise */
 NPY_NO_EXPORT int
 arrays_overlap(PyArrayObject *arr1, PyArrayObject *arr2)
 {
-    npy_uintp start1 = 0, start2 = 0, end1 = 0, end2 = 0;
-
-    get_array_memory_extents(arr1, &start1, &end1);
-    get_array_memory_extents(arr2, &start2, &end2);
+    mem_overlap_t result;
 
-    return (start1 < end2) && (start2 < end1);
+    result = solve_may_share_memory(arr1, arr2, NPY_MAY_SHARE_BOUNDS);
+    if (result == MEM_OVERLAP_NO) {
+        return 0;
+    }
+    else {
+        return 1;
+    }
 }

numpy/core/src/multiarray/arrayobject.c

@@ -51,6 +51,7 @@ maintainer email:  oliphant.travis@ieee.org
 #include "buffer.h"
 #include "array_assign.h"
 #include "alloc.h"
+#include "mem_overlap.h"
 
 /*NUMPY_API
   Compute the size of an array (in number of items)

numpy/core/src/multiarray/common.c

@@ -765,39 +765,6 @@ _IsWriteable(PyArrayObject *ap)
 }
 
 
-/* Gets a half-open range [start, end) of offsets from the data pointer */
-NPY_NO_EXPORT void
-offset_bounds_from_strides(const int itemsize, const int nd,
-                           const npy_intp *dims, const npy_intp *strides,
-                           npy_intp *lower_offset, npy_intp *upper_offset) {
-    npy_intp max_axis_offset;
-    npy_intp lower = 0;
-    npy_intp upper = 0;
-    int i;
-
-    for (i = 0; i < nd; i++) {
-        if (dims[i] == 0) {
-            /* If the array size is zero, return an empty range */
-            *lower_offset = 0;
-            *upper_offset = 0;
-            return;
-        }
-        /* Expand either upwards or downwards depending on stride */
-        max_axis_offset = strides[i] * (dims[i] - 1);
-        if (max_axis_offset > 0) {
-            upper += max_axis_offset;
-        }
-        else {
-            lower += max_axis_offset;
-        }
-    }
-    /* Return a half-open range */
-    upper += itemsize;
-    *lower_offset = lower;
-    *upper_offset = upper;
-}
-
-
 /**
  * Convert an array shape to a string such as "(1, 2)".
  *

numpy/core/src/multiarray/common.h

@@ -65,11 +65,6 @@ _IsAligned(PyArrayObject *ap);
 NPY_NO_EXPORT npy_bool
 _IsWriteable(PyArrayObject *ap);
 
-NPY_NO_EXPORT void
-offset_bounds_from_strides(const int itemsize, const int nd,
-                           const npy_intp *dims, const npy_intp *strides,
-                           npy_intp *lower_offset, npy_intp *upper_offset);
-
 NPY_NO_EXPORT PyObject *
 convert_shape_to_string(npy_intp n, npy_intp *vals, char *ending);
 

numpy/core/src/multiarray/getset.c

@@ -17,6 +17,7 @@
 #include "descriptor.h"
 #include "getset.h"
 #include "arrayobject.h"
+#include "mem_overlap.h"
 
 /*******************  array attribute get and set routines ******************/