Repository: cupy/cupy
Pull request: Add `order` argument to` empty_like`,` zeros_like`, etc.
closes #1654

Problem statement:
Feature: preserve memory layout order in empty_like, etc?
One thing that was surprising to me on using `cupy.zeros_like`, `cupy.empty_like` or `cupy.ones_like` were that when the input array was Fortran-ordered, the output is C-ordered. This differs from the NumPy behavior which preserves order. I see that it the docstring does state that "this function does not support order and subok options", so I suppose this is to be expected.

Would the CuPy devs be open to a PR that would preserve order 'C' or 'F' for contiguous arrays in these functions?. I think supporting the `order` parameter for values of 'C', 'F' and 'A' should be easy to do and would be willing to make a PR for that. (The underlying `cupy.zeros`, etc. already have the necessary `order` argument).

p.s. 
Congrats on the CuPy project! It is very nice and with the recently added support for custom kernels, I think I would prefer it going forward to PyCUDA in many cases due to the more complete NumPy API coverage.

Hints from issue discussion:

Commit message: pep8 fix

Patch to review:
diff --git a/cupy/creation/basic.py b/cupy/creation/basic.py
index ba429f56078..13cff939c6b 100644
--- a/cupy/creation/basic.py
+++ b/cupy/creation/basic.py
@@ -20,7 +20,56 @@ def empty(shape, dtype=float, order='C'):
     return cupy.ndarray(shape, dtype, order=order)
 
 
-def empty_like(a, dtype=None):
+def _new_like_order_and_strides(a, dtype, order):
+    """
+    Determine order and strides as in NumPy's PyArray_NewLikeArray.
+
+    (see: numpy/core/src/multiarray/ctors.c)
+    """
+    if order == 'A':
+        if a.flags.f_contiguous:
+            order = 'F'
+        else:
+            order = 'C'
+    elif order == 'K':
+        if a.flags.c_contiguous or a.ndim <= 1:
+            order = 'C'
+        elif a.flags.f_contiguous:
+            order = 'F'
+
+    if order in ['C', 'F']:
+        return order, None, None
+
+    elif order == 'K':
+        """
+        stable sort of strides in descending order
+        axes are included so the sort will always be stable
+        (mimics numpy's PyArray_CreateSortedStridePerm)
+        """
+        tmp = numpy.empty(len(a.strides),
+                          dtype=[('strides', numpy.intp),
+                                 ('axes', numpy.intp)])
+        tmp['strides'] = a.strides
+        tmp['axes'] = numpy.arange(a.ndim, dtype=numpy.intp)
+        tmp_sorted = numpy.sort(tmp, order=['strides', 'axes'])[::-1]
+        perm = tmp_sorted['axes']
+
+        # fill in strides based on the sorted order
+        order = 'C'
+        stride = numpy.dtype(dtype).itemsize
+        strides = numpy.zeros(a.ndim, dtype=numpy.intp)
+        for idim in range(a.ndim-1, -1, -1):
+            i_perm = perm[idim]
+            strides[i_perm] = stride
+            stride *= a.shape[i_perm]
+
+        memptr = cupy.empty(a.size, dtype=dtype).data
+        return order, strides, memptr
+    else:
+        raise ValueError("unknown order, {}".format(order))
+
+
+def empty_like(a, dtype=None, order='K'):
     """Returns a new array with same shape and dtype of a given array.
 
     This function currently does not support ``order`` and ``subok`` options.
@@ -28,6 +77,10 @@ def empty_like(a, dtype=None):
     Args:
         a (cupy.ndarray): Base array.
         dtype: Data type specifier. The data type of ``a`` is used by default.
+        order ({'C', 'F', 'A', or 'K'}): Overrides the memory layout of the
+            result. 'C' means C-order, 'F' means F-order, 'A' means 'F' if
+            ``a`` is Fortran contiguous, 'C' otherwise. 'K' means match the
+            layout of ``a`` as closely as possible.
 
     Returns:
         cupy.ndarray: A new array with same shape and dtype of ``a`` with
@@ -36,10 +89,11 @@ def empty_like(a, dtype=None):
     .. seealso:: :func:`numpy.empty_like`
 
     """
-    # TODO(beam2d): Support ordering option
     if dtype is None:
         dtype = a.dtype
-    return cupy.ndarray(a.shape, dtype)
+
+    order, strides, memptr = _new_like_order_and_strides(a, dtype, order)
+    return cupy.ndarray(a.shape, dtype, memptr, strides, order)
 
 
 def eye(N, M=None, k=0, dtype=float):
@@ -106,7 +160,7 @@ def ones(shape, dtype=float):
     return a
 
 
-def ones_like(a, dtype=None):
+def ones_like(a, dtype=None, order='K'):
     """Returns an array of ones with same shape and dtype as a given array.
 
     This function currently does not support ``order`` and ``subok`` options.
@@ -114,6 +168,10 @@ def ones_like(a, dtype=None):
     Args:
         a (cupy.ndarray): Base array.
         dtype: Data type specifier. The dtype of ``a`` is used by default.
+        order ({'C', 'F', 'A', or 'K'}): Overrides the memory layout of the
+            result. 'C' means C-order, 'F' means F-order, 'A' means 'F' if
+            ``a`` is Fortran contiguous, 'C' otherwise. 'K' means match the
+            layout of ``a`` as closely as possible.
 
     Returns:
         cupy.ndarray: An array filled with ones.
@@ -121,10 +179,10 @@ def ones_like(a, dtype=None):
     .. seealso:: :func:`numpy.ones_like`
 
     """
-    # TODO(beam2d): Support ordering option
     if dtype is None:
         dtype = a.dtype
-    a = cupy.ndarray(a.shape, dtype)
+    order, strides, memptr = _new_like_order_and_strides(a, dtype, order)
+    a = cupy.ndarray(a.shape, dtype, memptr, strides, order)
     a.fill(1)
     return a
 
@@ -149,7 +207,7 @@ def zeros(shape, dtype=float, order='C'):
     return a
 
 
-def zeros_like(a, dtype=None):
+def zeros_like(a, dtype=None, order='K'):
     """Returns an array of zeros with same shape and dtype as a given array.
 
     This function currently does not support ``order`` and ``subok`` options.
@@ -157,6 +215,10 @@ def zeros_like(a, dtype=None):
     Args:
         a (cupy.ndarray): Base array.
         dtype: Data type specifier. The dtype of ``a`` is used by default.
+        order ({'C', 'F', 'A', or 'K'}): Overrides the memory layout of the
+            result. 'C' means C-order, 'F' means F-order, 'A' means 'F' if
+            ``a`` is Fortran contiguous, 'C' otherwise. 'K' means match the
+            layout of ``a`` as closely as possible.\
 
     Returns:
         cupy.ndarray: An array filled with zeros.
@@ -164,10 +226,10 @@ def zeros_like(a, dtype=None):
     .. seealso:: :func:`numpy.zeros_like`
 
     """
-    # TODO(beam2d): Support ordering option
     if dtype is None:
         dtype = a.dtype
-    a = cupy.ndarray(a.shape, dtype)
+    order, strides, memptr = _new_like_order_and_strides(a, dtype, order)
+    a = cupy.ndarray(a.shape, dtype, memptr, strides, order)
     a.data.memset_async(0, a.nbytes)
     return a
 
@@ -199,7 +261,7 @@ def full(shape, fill_value, dtype=None):
     return a
 
 
-def full_like(a, fill_value, dtype=None):
+def full_like(a, fill_value, dtype=None, order='K'):
     """Returns a full array with same shape and dtype as a given array.
 
     This function currently does not support ``order`` and ``subok`` options.
@@ -208,6 +270,10 @@ def full_like(a, fill_value, dtype=None):
         a (cupy.ndarray): Base array.
         fill_value: A scalar value to fill a new array.
         dtype: Data type specifier. The dtype of ``a`` is used by default.
+        order ({'C', 'F', 'A', or 'K'}): Overrides the memory layout of the
+            result. 'C' means C-order, 'F' means F-order, 'A' means 'F' if
+            ``a`` is Fortran contiguous, 'C' otherwise. 'K' means match the
+            layout of ``a`` as closely as possible.
 
     Returns:
         cupy.ndarray: An array filled with ``fill_value``.
@@ -215,9 +281,9 @@ def full_like(a, fill_value, dtype=None):
     .. seealso:: :func:`numpy.full_like`
 
     """
-    # TODO(beam2d): Support ordering option
     if dtype is None:
         dtype = a.dtype
-    a = cupy.ndarray(a.shape, dtype)
+    order, strides, memptr = _new_like_order_and_strides(a, dtype, order)
+    a = cupy.ndarray(a.shape, dtype, memptr, strides, order)
     a.fill(fill_value)
     return a
diff --git a/tests/cupy_tests/creation_tests/test_basic.py b/tests/cupy_tests/creation_tests/test_basic.py
index 7ee92eb1538..5ebd7383eb5 100644
--- a/tests/cupy_tests/creation_tests/test_basic.py
+++ b/tests/cupy_tests/creation_tests/test_basic.py
@@ -77,6 +77,102 @@ def test_empty_like(self, xp, dtype):
         b.fill(0)
         return b
 
+    @testing.for_all_dtypes()
+    @testing.numpy_cupy_array_equal()
+    def test_empty_like_C_order(self, xp, dtype):
+        a = testing.shaped_arange((2, 3, 4), xp, dtype)
+        a = xp.asfortranarray(a)
+        b = xp.empty_like(a, order='C')
+        b.fill(0)
+        self.assertTrue(b.flags.c_contiguous)
+        return b
+
+    @testing.for_all_dtypes()
+    @testing.numpy_cupy_array_equal()
+    def test_empty_like_F_order(self, xp, dtype):
+        a = testing.shaped_arange((2, 3, 4), xp, dtype)
+        b = xp.empty_like(a, order='F')
+        b.fill(0)
+        self.assertTrue(b.flags.f_contiguous)
+        return b
+
+    @testing.for_all_dtypes()
+    @testing.numpy_cupy_array_equal()
+    def test_empty_like_A_order(self, xp, dtype):
+        a = testing.shaped_arange((2, 3, 4), xp, dtype)
+        b = xp.empty_like(a, order='A')
+        b.fill(0)
+        self.assertTrue(b.flags.c_contiguous)
+        return b
+
+    @testing.for_all_dtypes()
+    @testing.numpy_cupy_array_equal()
+    def test_empty_like_A_order2(self, xp, dtype):
+        a = testing.shaped_arange((2, 3, 4), xp, dtype)
+        a = xp.asfortranarray(a)
+        b = xp.empty_like(a, order='A')
+        b.fill(0)
+        self.assertTrue(b.flags.f_contiguous)
+        return b
+
+    @testing.for_all_dtypes()
+    @testing.numpy_cupy_array_equal()
+    def test_empty_like_A_order3(self, xp, dtype):
+        a = testing.shaped_arange((2, 3, 4), xp, dtype)
+        b = xp.empty_like(a[:, ::2, :], order='A')
+        b.fill(0)
+        self.assertTrue(b.flags.c_contiguous)
+        return b
+
+    @testing.for_all_dtypes()
+    @testing.numpy_cupy_array_equal()
+    def test_empty_like_K_order(self, xp, dtype):
+        a = testing.shaped_arange((2, 3, 4), xp, dtype)
+        b = xp.empty_like(a, order='K')
+        b.fill(0)
+        self.assertTrue(b.flags.c_contiguous)
+        return b
+
+    @testing.for_all_dtypes()
+    @testing.numpy_cupy_array_equal()
+    def test_empty_like_K_order2(self, xp, dtype):
+        a = testing.shaped_arange((2, 3, 4), xp, dtype)
+        a = xp.asfortranarray(a)
+        b = xp.empty_like(a, order='K')
+        b.fill(0)
+        self.assertTrue(b.flags.f_contiguous)
+        return b
+
+    @testing.for_all_dtypes()
+    @testing.numpy_cupy_array_equal()
+    def test_empty_like_K_order3(self, xp, dtype):
+        a = testing.shaped_arange((2, 3, 4), xp, dtype)
+        # test strides that are both non-contiguous and non-descending
+        a = a[:, ::2, :].swapaxes(0, 1)
+        b = xp.empty_like(a, order='K')
+        b.fill(0)
+        self.assertFalse(b.flags.c_contiguous)
+        self.assertFalse(b.flags.f_contiguous)
+        return b
+
+    @testing.for_all_dtypes()
+    def test_empty_like_K_strides(self, dtype):
+        # test strides that are both non-contiguous and non-descending
+        a = testing.shaped_arange((2, 3, 4), numpy, dtype)
+        a = a[:, ::2, :].swapaxes(0, 1)
+        b = numpy.empty_like(a, order='K')
+        b.fill(0)
+
+        # GPU case
+        ag = testing.shaped_arange((2, 3, 4), cupy, dtype)
+        ag = ag[:, ::2, :].swapaxes(0, 1)
+        bg = cupy.empty_like(ag, order='K')
+        bg.fill(0)
+
+        # make sure NumPy and CuPy strides agree
+        self.assertEqual(b.strides, bg.strides)
+        return
+
     @testing.for_CF_orders()
     def test_empty_zero_sized_array_strides(self, order):
         a = numpy.empty((1, 0, 2), dtype='d', order=order)


Code review comments:
diff_hunk: @@ -20,14 +20,67 @@ def empty(shape, dtype=float, order='C'):
     return cupy.ndarray(shape, dtype, order=order)
 
 
-def empty_like(a, dtype=None):
+def _new_like_order_and_strides(a, dtype, order):
+    """
+    Determine order and strides as in NumPy's PyArray_NewLikeArray.
+
+    (see: numpy/core/src/multiarray/ctors.c)
+    """
+    if order == 'A':
+        if a.flags.f_contiguous:
+            order = 'F'
+        else:
+            order = 'C'
+    elif order == 'K':
+        if a.flags.c_contiguous or a.ndim <= 1:
+            order = 'C'
+        elif a.flags.f_contiguous:
+            order = 'F'
+
+    if order in ['C', 'F']:
+        return order, None, None
+
+    elif order == 'K':
+        """
+        stable sort of strides in descending order
+        axes are included so the sort will always be stable
+        (mimics numpy's PyArray_CreateSortedStridePerm)
+        """
+        tmp = numpy.empty(len(a.strides),
+                          dtype=[('strides', numpy.intp),
+                                 ('axes', numpy.intp)])
+        tmp['strides'] = a.strides
+        tmp['axes'] = numpy.arange(a.ndim, dtype=numpy.intp)
+        tmp_sorted = numpy.sort(tmp, order=['strides', 'axes'])[::-1]
+        perm = tmp_sorted['axes']
+
+        # fill in strides based on the sorted order
+        order = 'C'
+        stride = numpy.dtype(dtype).itemsize
+        strides = numpy.zeros(a.ndim, dtype=numpy.intp)
+        for idim in range(a.ndim-1, -1, -1):
+            i_perm = perm[idim]
+            strides[i_perm] = stride
+            stride *= a.shape[i_perm]
+
+        memptr = cupy.empty(a.size, dtype=dtype).data
+        return order, strides, memptr
+    else:
+        raise ValueError("unknown order, {}".format(order))
path: cupy/creation/basic.py
text: @user1:
Let's align to NumPy's behavior.
```suggestion
        raise TypeError('order not understood: {}'.format(order))
```

@author:
done