Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

ENH: Add shape to *_like() array creation #13046

Merged
merged 24 commits into from Apr 25, 2019
Merged
Show file tree
Hide file tree
Changes from 17 commits
Commits
Show all changes
24 commits
Select commit Hold shift + click to select a range
f885c08
ENH: Added shape argument to *_like() array creation functions
pentschev Feb 26, 2019
5ad61b6
ENH: C backend adjustments for shape argument on *_like()
pentschev Feb 26, 2019
fd3e270
TST: Added test for shape argument in *_like() functions
pentschev Feb 26, 2019
b7202a7
ENH: Added PyArray_NewLikeArrayWithShape()
pentschev Feb 26, 2019
013cbce
BUG: Fix for PyArray_NewLikeArrayWithShape strides and ndim == 0
pentschev Feb 27, 2019
807f512
REL: Updates for C-API, version 1.17.x
pentschev Feb 27, 2019
95bbfd0
Revert "REL: Updates for C-API, version 1.17.x"
pentschev Mar 13, 2019
43b9828
Revert exposing PyArray_NewLikeArrayWithShape on C-API
pentschev Mar 13, 2019
fa7fd75
DOC: fix versionadded for *_like() shape argument
pentschev Mar 13, 2019
d57e6d3
STY: add missing spaces in array initializers
pentschev Mar 13, 2019
f1b3e91
ENH: empty_like raises ValueError
pentschev Mar 16, 2019
d24ac10
TST: test for exception of *_like() functions
pentschev Mar 16, 2019
ec26417
DOC: release note for shape argument in *_like() functions
pentschev Mar 16, 2019
40e7e9e
DOC: fix *_like() documentation on raises
pentschev Mar 17, 2019
928952d
BUG: *_like() raises for non-C/F-layout arrays
pentschev Mar 17, 2019
e394356
TST: change *_like() shapes to prevent NPY_RELAXED_STRIDE_DEBUG=1 fai…
pentschev Mar 17, 2019
d938fb9
Move empty_like() exception to C implementation
pentschev Mar 19, 2019
99a55e0
Merge branch 'master' into add-shape-to-like-array-creation
pentschev Mar 28, 2019
c087801
Update *_like() ValueError documentation
pentschev Apr 4, 2019
db7614b
Merge remote-tracking branch 'upstream/master' into add-shape-to-like…
pentschev Apr 5, 2019
669ea71
Rearrange stride computation for *_like() if new shape and order='K'
pentschev Apr 10, 2019
f53afbe
Change handling of order= for *_like()
pentschev Apr 19, 2019
c3ac08e
Merge branch 'master' into add-shape-to-like-array-creation
pentschev Apr 19, 2019
695b836
Fix *_like() tests
pentschev Apr 19, 2019
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Jump to
Jump to file
Failed to load files.
Diff view
Diff view
9 changes: 9 additions & 0 deletions doc/release/1.17.0-notes.rst
Expand Up @@ -75,6 +75,15 @@ divmod operation is now supported for two ``timedelta64`` operands
The divmod operator now handles two ``np.timedelta64`` operands, with
type signature mm->qm.

``np.empty_like`` and related functions now accept a ``shape`` argument
-----------------------------------------------------------------------
``np.empty_like``, ``np.full_like``, ``np.ones_like`` and ``np.zeros_like`` now
accept a ``shape`` keyword argument, which can be used to create a new array
as the prototype, overriding its shape as well. This is particularly useful
when combined with the ``__array_function__`` protocol, allowing the creation
of new arbitrary-shape arrays from NumPy-like libraries when such an array
is used as the prototype.


Improvements
============
Expand Down
14 changes: 12 additions & 2 deletions numpy/core/multiarray.py
Expand Up @@ -71,9 +71,9 @@


@array_function_from_c_func_and_dispatcher(_multiarray_umath.empty_like)
def empty_like(prototype, dtype=None, order=None, subok=None):
def empty_like(prototype, dtype=None, order=None, subok=None, shape=None):
"""
empty_like(prototype, dtype=None, order='K', subok=True)
empty_like(prototype, dtype=None, order='K', subok=True, shape=None)
Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I just noticed order and subok both default to None rather than the ('K' and True) values documentation proposes and other _like() functions have. Which one is the correct? It seems to me that the function definition should be order='K' and subok=True, am I correct?

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

It's a little confusing, but the arguments for this function (the dispatcher) need to default to None. We actually have a test that verifies this.

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Should we then update the documentation to reflect the real defaults?


Return a new array with the same shape and type as a given array.

Expand All @@ -97,13 +97,23 @@ def empty_like(prototype, dtype=None, order=None, subok=None):
If True, then the newly created array will use the sub-class
type of 'a', otherwise it will be a base-class array. Defaults
to True.
shape : int or sequence of ints, optional.
Overrides the shape of the result.

pentschev marked this conversation as resolved.
Show resolved Hide resolved
.. versionadded:: 1.17.0

Returns
-------
out : ndarray
Array of uninitialized (arbitrary) data with the same
shape and type as `prototype`.

Raises
------
ValueError
If len(shape) different from prototype.ndim and order is 'K', or
prototype is not a C/F-layout array

See Also
--------
ones_like : Return an array of ones with shape and type of input.
Expand Down
48 changes: 39 additions & 9 deletions numpy/core/numeric.py
Expand Up @@ -106,12 +106,12 @@ class ComplexWarning(RuntimeWarning):
pass


def _zeros_like_dispatcher(a, dtype=None, order=None, subok=None):
def _zeros_like_dispatcher(a, dtype=None, order=None, subok=None, shape=None):
return (a,)


@array_function_dispatch(_zeros_like_dispatcher)
def zeros_like(a, dtype=None, order='K', subok=True):
def zeros_like(a, dtype=None, order='K', subok=True, shape=None):
"""
Return an array of zeros with the same shape and type as a given array.

Expand All @@ -135,12 +135,22 @@ def zeros_like(a, dtype=None, order='K', subok=True):
If True, then the newly created array will use the sub-class
type of 'a', otherwise it will be a base-class array. Defaults
to True.
shape : int or sequence of ints, optional.
Overrides the shape of the result.

pentschev marked this conversation as resolved.
Show resolved Hide resolved
.. versionadded:: 1.17.0

Returns
-------
out : ndarray
Array of zeros with the same shape and type as `a`.

Raises
------
ValueError
If len(shape) different from a.ndim and order is 'K', or
a is not a C/F-layout array
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

What does C/F-layout mean here?

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I assume you mean the order argument has a invalid value? If so, maybe

        If `order` has an invalid value or `order` is 'K' and ``len(shape)``
        is different from ``a.ndim``.

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I think what you mean with your suggestion is that order could be any random value different than CFAK, right? If that's the case, that's not what I my intent with the documentation, what I mean is that there's no particular memory layout, maybe I got the terminology wrong?

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

And by no particular memory layout, I mean something like the result of advanced indexing, for example, stealing that terminology from indexing documentation:

The memory layout of an advanced indexing result is optimized for each indexing operation and no particular memory order can be assumed.

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

OK, I see what you mean. That is difficult to express exactly, maybe just brute force it:

        If `order` is 'K'  and any of the following is true

            - ``len(shape) != a.ndim``,
            - `a` is neither C nor F contiguous.

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Note that strictly speaking, it should be possible to add 1's to the shape without affecting the contiguity.

In [1]: a = ones((2,2), order='C')

In [2]: a.flags.c_contiguous
Out[2]: True

In [3]: a[None,:,:].flags.c_contiguous
Out[3]: True

In [4]: a = ones((2,2), order='F')

In [5]: a.flags.f_contiguous
Out[5]: True

In [6]: a[None,:,:].flags.f_contiguous
Out[6]: True

So possibly the correct check would be that the number of number of dimensions > 1 must match. @seberg Thoughts?

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

You're right, I can do that modification if necessary, but I'm not sure this is of utmost importance. Particularly, for our use cases I don't think that affects us.


See Also
--------
empty_like : Return an empty array with shape and type of input.
Expand All @@ -166,7 +176,7 @@ def zeros_like(a, dtype=None, order='K', subok=True):
array([0., 0., 0.])

"""
res = empty_like(a, dtype=dtype, order=order, subok=subok)
res = empty_like(a, dtype=dtype, order=order, subok=subok, shape=shape)
# needed instead of a 0 to get same result as zeros for for string dtypes
z = zeros(1, dtype=res.dtype)
multiarray.copyto(res, z, casting='unsafe')
rgommers marked this conversation as resolved.
Show resolved Hide resolved
Expand Down Expand Up @@ -226,12 +236,12 @@ def ones(shape, dtype=None, order='C'):
return a


def _ones_like_dispatcher(a, dtype=None, order=None, subok=None):
def _ones_like_dispatcher(a, dtype=None, order=None, subok=None, shape=None):
return (a,)


@array_function_dispatch(_ones_like_dispatcher)
def ones_like(a, dtype=None, order='K', subok=True):
def ones_like(a, dtype=None, order='K', subok=True, shape=None):
"""
Return an array of ones with the same shape and type as a given array.

Expand All @@ -255,12 +265,22 @@ def ones_like(a, dtype=None, order='K', subok=True):
If True, then the newly created array will use the sub-class
type of 'a', otherwise it will be a base-class array. Defaults
to True.
shape : int or sequence of ints, optional.
Overrides the shape of the result.

.. versionadded:: 1.17.0

Returns
-------
out : ndarray
Array of ones with the same shape and type as `a`.

Raises
------
ValueError
If len(shape) different from a.ndim and order is 'K', or
a is not a C/F-layout array

See Also
--------
empty_like : Return an empty array with shape and type of input.
Expand All @@ -286,7 +306,7 @@ def ones_like(a, dtype=None, order='K', subok=True):
array([1., 1., 1.])

"""
res = empty_like(a, dtype=dtype, order=order, subok=subok)
res = empty_like(a, dtype=dtype, order=order, subok=subok, shape=shape)
multiarray.copyto(res, 1, casting='unsafe')
rgommers marked this conversation as resolved.
Show resolved Hide resolved
return res

Expand Down Expand Up @@ -338,12 +358,12 @@ def full(shape, fill_value, dtype=None, order='C'):
return a


def _full_like_dispatcher(a, fill_value, dtype=None, order=None, subok=None):
def _full_like_dispatcher(a, fill_value, dtype=None, order=None, subok=None, shape=None):
return (a,)


@array_function_dispatch(_full_like_dispatcher)
def full_like(a, fill_value, dtype=None, order='K', subok=True):
def full_like(a, fill_value, dtype=None, order='K', subok=True, shape=None):
"""
Return a full array with the same shape and type as a given array.

Expand All @@ -365,12 +385,22 @@ def full_like(a, fill_value, dtype=None, order='K', subok=True):
If True, then the newly created array will use the sub-class
type of 'a', otherwise it will be a base-class array. Defaults
to True.
shape : int or sequence of ints, optional.
Overrides the shape of the result.

.. versionadded:: 1.17.0

Returns
-------
out : ndarray
Array of `fill_value` with the same shape and type as `a`.

Raises
------
ValueError
If len(shape) different from a.ndim and order is 'K', or
a is not a C/F-layout array

See Also
--------
empty_like : Return an empty array with shape and type of input.
Expand All @@ -395,7 +425,7 @@ def full_like(a, fill_value, dtype=None, order='K', subok=True):
array([0.1, 0.1, 0.1, 0.1, 0.1, 0.1])

"""
res = empty_like(a, dtype=dtype, order=order, subok=subok)
res = empty_like(a, dtype=dtype, order=order, subok=subok, shape=shape)
multiarray.copyto(res, fill_value, casting='unsafe')
rgommers marked this conversation as resolved.
Show resolved Hide resolved
return res

Expand Down
61 changes: 49 additions & 12 deletions numpy/core/src/multiarray/ctors.c
Expand Up @@ -1183,28 +1183,43 @@ PyArray_NewFromDescrAndBase(
flags, obj, base, 0, 0);
}

/*NUMPY_API
/*
* Creates a new array with the same shape as the provided one,
* with possible memory layout order and data type changes.
* with possible memory layout order, data type and shape changes.
*
* prototype - The array the new one should be like.
* order - NPY_CORDER - C-contiguous result.
* NPY_FORTRANORDER - Fortran-contiguous result.
* NPY_ANYORDER - Fortran if prototype is Fortran, C otherwise.
* NPY_KEEPORDER - Keeps the axis ordering of prototype.
* dtype - If not NULL, overrides the data type of the result.
* ndim - If not 0 and dims not NULL, overrides the shape of the result.
* dims - If not NULL and ndim not 0, overrides the shape of the result.
* subok - If 1, use the prototype's array subtype, otherwise
* always create a base-class array.
*
* NOTE: If dtype is not NULL, steals the dtype reference. On failure or when
* dtype->subarray is true, dtype will be decrefed.
*/
NPY_NO_EXPORT PyObject *
PyArray_NewLikeArray(PyArrayObject *prototype, NPY_ORDER order,
PyArray_Descr *dtype, int subok)
PyArray_NewLikeArrayWithShape(PyArrayObject *prototype, NPY_ORDER order,
PyArray_Descr *dtype, int ndim, npy_intp *dims, int subok)
{
PyObject *ret = NULL;
int ndim = PyArray_NDIM(prototype);

int new_shape = 1;

if (dims == NULL) {
ndim = PyArray_NDIM(prototype);
dims = PyArray_DIMS(prototype);
new_shape = 0;
}
else if (order == NPY_KEEPORDER && (ndim != PyArray_NDIM(prototype) ||
Copy link
Member

@charris charris Apr 4, 2019

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This breaks previous behavior which did not require the prototype to be contiguous.

In [1]: a = ones((3,4))[:,::2]

In [2]: a.flags
Out[2]: 
  C_CONTIGUOUS : False
  F_CONTIGUOUS : False
  OWNDATA : False
  WRITEABLE : True
  ALIGNED : True
  WRITEBACKIFCOPY : False
  UPDATEIFCOPY : False

In [3]: ones_like(a, order='K')
Out[3]: 
array([[1., 1.],
       [1., 1.],
       [1., 1.]])

I'm not sure why an error should be raised, why not just override the shape and go with C order when a new shape is specified and order == 'K'?

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@charris what you're saying is like my original implementation. I modified this based on @eric-wieser's request here #13046 (comment). For me it's fine either way.

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@charris: This shouldn't break any existing code once my comment below is addressed. The second check here makes no sense to me either.

The thinking is that requesting the order (of the strides) to be kept is nonsensical if the the number of strides is not the same.

!(PyArray_IS_C_CONTIGUOUS(prototype) || PyArray_IS_F_CONTIGUOUS(prototype)))) {
PyErr_SetString(PyExc_ValueError,
"mismatching ndim or non-C/F-layout arrays can not keep order");
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Why are you checking PyArray_IS_C_CONTIGUOUS and PyArray_IS_F_CONTIGUOUS here?

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Because we can't keep the order when there's no particular memory layout, just like in the discussion here:
#13046 (comment)

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

My current understanding of 'K' is that we try to keep the memory order if possible, not in all circumstances. Indeed, the former code fell back on recomputing the strides when the prototype was neither c/f continguous.

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

First of all, we should make sure that this does not switch to an error suddenly. I agree that in most cases K probably means if possible.

I think the default may actually be K and not something like None (to mean K if possible), so that K is the most common option? If K was passed explicitly, there might be an argument.

There are some strides where it would be unambiguous, but I am not sure it is worth to bother about those.

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Unless I'm miserably misunderstanding some of the comments, I think that there are two clashing opinions:

  1. Keep order if possible (as my original implementation did), otherwise, fallback to something like C-contiguous (currently backed by @charris and @seberg);
  2. Raise an exception if the order is K but not possible to maintain it (backed by @eric-wieser, previously discussed in ENH: Add shape to *_like() array creation #13046 (comment)).

Feel free to correct me if I'm really misunderstanding it, otherwise, it would be nice to have here @charris, @eric-wieser and @seberg commenting here to reach an agreement before I do any further changes regarding this matter.

Copy link
Contributor Author

@pentschev pentschev Apr 5, 2019

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@eric-wieser could you let me know about your findings then? I'm sorry, but I don't understand what you mean with your previous statement, it may be getting a little too complex for my understanding of NumPy handling of strides and orders.

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Sure, here's the case I'm thinking of:

>>> a = np.zeros((20, 30, 40)).transpose((1, 0, 2))[::10,::10,::10]

>>> a.strides
(3200, 96000, 80)

>>> np.copy(a, order='K').strides
(32, 96, 8)

note how the order of the strides is kept.

This is the behavior I would expect for empty_like too

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@eric-wieser I believe the new commit I pushed covers now the case you're looking for. Extending the case from your previous comment, this is what we get with this new commit:

>>> np.empty_like(a, shape=(3, 2, 4), order='K').strides
(32, 96, 8)

Can you confirm this is what you had in mind?

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

That behavior looks correct to me now. Well take one more look at this tonight.

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@eric-wieser just checking, have you had the chance to have a look at this? Do you think we're good to merge this now?

return ret;
}

/* If no override data type, use the one from the prototype */
if (dtype == NULL) {
Expand Down Expand Up @@ -1232,12 +1247,12 @@ PyArray_NewLikeArray(PyArrayObject *prototype, NPY_ORDER order,
break;
}

/* If it's not KEEPORDER, this is simple */
if (order != NPY_KEEPORDER) {
/* If it's not KEEPORDER, or there is a shape change, this is simple */
if (order != NPY_KEEPORDER || new_shape) {
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I think if the number of dims match, matching the stride order is still a good idea

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

But passing shape alone doesn't provide any information of strides, it will not be more than an assumption that keeping strides unchanged is important. Plus, we could only keep the strides for the dimensions that don't change.

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

The stride information comes from the source. The purpose of KEEPORDER is to preserve np.argsort(arr.strides).

I think that the behavior should be:

  • Same ndim - respect the stride ordering through keeporder
  • different ndim - reject keeporder since it is meaningless, throw an exception asking the use to use corder or Forder.

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

But even if ndim remains the same, suppose source shape=(5, 2) and new shape=(10000, 10), in such a case we can't respect strides. Another possibility is source shape=(10000, 10) and new shape=(5, 2), we could respect the strides at the cost of huge memory waste, which is most likely going to be an unexpected behavior from the user.

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

keeporder means keep the relative order of the strides (np.argsort(out.strides) == np.argsort(arr.strides)), not the strides themselves (out.strides == arr.strides).

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Oh, I'm sorry, I got what you mean now. Yes, I agree with that, let me change that.

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I think the new commits should cover this case now.

ret = PyArray_NewFromDescr(subok ? Py_TYPE(prototype) : &PyArray_Type,
dtype,
ndim,
PyArray_DIMS(prototype),
dims,
NULL,
NULL,
order,
Expand All @@ -1246,11 +1261,10 @@ PyArray_NewLikeArray(PyArrayObject *prototype, NPY_ORDER order,
/* KEEPORDER needs some analysis of the strides */
else {
npy_intp strides[NPY_MAXDIMS], stride;
npy_intp *shape = PyArray_DIMS(prototype);
npy_stride_sort_item strideperm[NPY_MAXDIMS];
int idim;

PyArray_CreateSortedStridePerm(PyArray_NDIM(prototype),
PyArray_CreateSortedStridePerm(ndim,
PyArray_STRIDES(prototype),
pentschev marked this conversation as resolved.
Show resolved Hide resolved
strideperm);

Expand All @@ -1259,14 +1273,14 @@ PyArray_NewLikeArray(PyArrayObject *prototype, NPY_ORDER order,
for (idim = ndim-1; idim >= 0; --idim) {
npy_intp i_perm = strideperm[idim].perm;
strides[i_perm] = stride;
stride *= shape[i_perm];
stride *= dims[i_perm];
}
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

My 1.ii proposal would look like changing everything between the new line 1261 and here to:

if (PyArray_NDIM(prototype) >= ndim) {
    int leading_dims = PyArray_NDIM(prototype) - ndim;

    /* Use only the trailing strides */
    PyArray_CreateSortedStridePerm(ndim,
                                   PyArray_STRIDES(prototype) + leading_dims,
                                   strideperm);
    /* Build the new strides */
    stride = dtype->elsize;
    for (idim = ndim-1; idim >= 0; --idim) {
        npy_intp i_perm = strideperm[idim].perm;
        strides[i_perm] = stride;
        stride *= shape[i_perm];
    }
}
else {
    int leading_dims = ndim - PyArray_NDIM(prototype);
    /* Use all the strides */
    PyArray_CreateSortedStridePerm(PyArray_NDIM(prototype),
                                   PyArray_STRIDES(prototype),
                                   strideperm);

    /* Build the new trailing strides */
    stride = dtype->elsize;
    for (idim = PyArray_NDIM(prototype)-1; idim >= 0; --idim) {
        npy_intp i_perm = strideperm[idim].perm + leading_dims;
        strides[i_perm] = stride;
        stride *= shape[i_perm];
    }

    /* Create remaining leading strides as C order */
    for (idim = leading_dims; idim >= 0; --idim) {
        strides[idim] = stride;
        stride *= shape[idim];
    }
}

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Or trading branching for state:

/* not sure if this is clearer via min/max */
int leading_src_dims = 0;  // max(src.ndim - dst.ndim, 0)
int leading_dst_dims = 0;  // max(dst.ndim - src.ndim, 0)
int shared_dims;  // min(src.ndim, dst.ndim
if (PyArray_NDIM(prototype) >= ndim) {
    shared_dims = ndim;
    leading_src_dims = PyArray_NDIM(prototype) - ndim;
}
else {
    shared_dims = PyArray_NDIM(prototype);
    leading_dst_dims = ndim - PyArray_NDIM(prototype);
}

/* Use only the trailing strides from the source */
PyArray_CreateSortedStridePerm(shared_dims,
                               PyArray_STRIDES(prototype) + leading_src_dims,
                               strideperm);

/* Build the destrination trailing strides */
stride = dtype->elsize;
for (idim = ndim-1; idim >= 0; --idim) {
    npy_intp i_perm = strideperm[idim].perm + leading_dst_dims;
    strides[i_perm] = stride;
    stride *= shape[i_perm];
}

/* Create remaining leading strides as C order */
for (idim = leading_dst_dims; idim >= 0; --idim) {
    strides[idim] = stride;
    stride *= shape[idim];
}


/* Finally, allocate the array */
ret = PyArray_NewFromDescr(subok ? Py_TYPE(prototype) : &PyArray_Type,
dtype,
ndim,
shape,
dims,
strides,
NULL,
0,
Expand All @@ -1276,6 +1290,29 @@ PyArray_NewLikeArray(PyArrayObject *prototype, NPY_ORDER order,
return ret;
}

/*NUMPY_API
* Creates a new array with the same shape as the provided one,
* with possible memory layout order and data type changes.
*
* prototype - The array the new one should be like.
* order - NPY_CORDER - C-contiguous result.
* NPY_FORTRANORDER - Fortran-contiguous result.
* NPY_ANYORDER - Fortran if prototype is Fortran, C otherwise.
* NPY_KEEPORDER - Keeps the axis ordering of prototype.
* dtype - If not NULL, overrides the data type of the result.
* subok - If 1, use the prototype's array subtype, otherwise
* always create a base-class array.
*
* NOTE: If dtype is not NULL, steals the dtype reference. On failure or when
* dtype->subarray is true, dtype will be decrefed.
*/
NPY_NO_EXPORT PyObject *
PyArray_NewLikeArray(PyArrayObject *prototype, NPY_ORDER order,
PyArray_Descr *dtype, int subok)
{
return PyArray_NewLikeArrayWithShape(prototype, order, dtype, 0, NULL, subok);
}

/*NUMPY_API
* Generic new array creation routine.
*/
Expand Down
4 changes: 4 additions & 0 deletions numpy/core/src/multiarray/ctors.h
Expand Up @@ -18,6 +18,10 @@ PyArray_NewFromDescr_int(PyTypeObject *subtype, PyArray_Descr *descr, int nd,
int flags, PyObject *obj, PyObject *base, int zeroed,
int allow_emptystring);

NPY_NO_EXPORT PyObject *
PyArray_NewLikeArrayWithShape(PyArrayObject *prototype, NPY_ORDER order,
PyArray_Descr *dtype, int ndim, npy_intp *dims, int subok);

NPY_NO_EXPORT PyObject *PyArray_New(PyTypeObject *, int nd, npy_intp *,
int, npy_intp *, void *, int, int, PyObject *);

Expand Down
23 changes: 14 additions & 9 deletions numpy/core/src/multiarray/multiarraymodule.c
Expand Up @@ -1758,7 +1758,7 @@ static PyObject *
array_copyto(PyObject *NPY_UNUSED(ignored), PyObject *args, PyObject *kwds)
{

static char *kwlist[] = {"dst","src","casting","where",NULL};
static char *kwlist[] = {"dst", "src", "casting", "where", NULL};
PyObject *wheremask_in = NULL;
PyArrayObject *dst = NULL, *src = NULL, *wheremask = NULL;
NPY_CASTING casting = NPY_SAME_KIND_CASTING;
Expand Down Expand Up @@ -1803,7 +1803,7 @@ static PyObject *
array_empty(PyObject *NPY_UNUSED(ignored), PyObject *args, PyObject *kwds)
{

static char *kwlist[] = {"shape","dtype","order",NULL};
static char *kwlist[] = {"shape", "dtype", "order", NULL};
PyArray_Descr *typecode = NULL;
PyArray_Dims shape = {NULL, 0};
NPY_ORDER order = NPY_CORDER;
Expand Down Expand Up @@ -1846,23 +1846,28 @@ static PyObject *
array_empty_like(PyObject *NPY_UNUSED(ignored), PyObject *args, PyObject *kwds)
{

static char *kwlist[] = {"prototype","dtype","order","subok",NULL};
static char *kwlist[] = {"prototype", "dtype", "order", "subok", "shape", NULL};
PyArrayObject *prototype = NULL;
PyArray_Descr *dtype = NULL;
NPY_ORDER order = NPY_KEEPORDER;
PyArrayObject *ret = NULL;
int subok = 1;
PyArray_Dims shape = {NULL, 0};

if (!PyArg_ParseTupleAndKeywords(args, kwds, "O&|O&O&i:empty_like", kwlist,
if (!PyArg_ParseTupleAndKeywords(args, kwds, "O&|O&O&iO&:empty_like", kwlist,
&PyArray_Converter, &prototype,
&PyArray_DescrConverter2, &dtype,
&PyArray_OrderConverter, &order,
&subok)) {
&subok,
&PyArray_IntpConverter, &shape)) {
goto fail;
}
/* steals the reference to dtype if it's not NULL */
ret = (PyArrayObject *)PyArray_NewLikeArray(prototype,
order, dtype, subok);
ret = (PyArrayObject *)PyArray_NewLikeArrayWithShape(prototype, order, dtype,
shape.len, shape.ptr, subok);
if (!ret) {
goto fail;
}
Py_DECREF(prototype);

return (PyObject *)ret;
Expand All @@ -1881,7 +1886,7 @@ static PyObject *
array_scalar(PyObject *NPY_UNUSED(ignored), PyObject *args, PyObject *kwds)
{

static char *kwlist[] = {"dtype","obj", NULL};
static char *kwlist[] = {"dtype", "obj", NULL};
PyArray_Descr *typecode;
PyObject *obj = NULL, *tmpobj = NULL;
int alloc = 0;
Expand Down Expand Up @@ -1957,7 +1962,7 @@ array_scalar(PyObject *NPY_UNUSED(ignored), PyObject *args, PyObject *kwds)
static PyObject *
array_zeros(PyObject *NPY_UNUSED(ignored), PyObject *args, PyObject *kwds)
{
static char *kwlist[] = {"shape","dtype","order",NULL};
static char *kwlist[] = {"shape", "dtype", "order", NULL};
PyArray_Descr *typecode = NULL;
PyArray_Dims shape = {NULL, 0};
NPY_ORDER order = NPY_CORDER;
Expand Down