I worry this patch, by first creating an object array, and then converting to scalars, will make the common case where the output does have the correct type quite a bit slower. Also, do you have an example where this is an issue? It looks as if the function itself has to give different types of outputs depending on its inputs, which seems strange. If this is really an issue, I would suggest that for these rare cases, one should advise to wrap the function into another one that ensures the output is of the right type.

@mhvk : Did you read my PR description? All of your questions are answered there. Not to mention, your suggestion would not be feasible because the damage would already have been done in terms of element coercion when calling apply_along_axis. Wrapping this function would be moot.

Nevertheless, your point about performance is well taken. There is a noticeable performance hit:

In [1]: import numpy as np
In [2]: arr = np.arange(1000000)
In [3]: %timeit np.apply_along_axis(lambda x: x + 1, 0, arr) 

master : 5.93 ms / loop
PR : 82.8 ms / loop

In [4]: %timeit np.apply_along_axis(lambda x: [x + 1, 5], 0, arr)

master : 1.88 ms / loop
PR : 2.28 ms / loop

@gfyoung - yes, I did read your description, but as I said I do not understand how one could easily trigger the sitation you are trying to solve, i.e., how one can have a situation where the returned dtype will depend on which part of the array one is applying the function to. Could you give an example of a function for which the issue arises?

@mhvk : But then your question about the issue arising should have been answered. Just read the issue I pointed to in the description. Not to mention, my first sentence in the description clearly explains the dtype casting behavior in np.apply_along_axis IMO.

Duh... somehow I missed that link completely. But it does confirm what I felt already: that one is trying to solve corner cases that would be better dealt with in an update to the documentation rather than in a change that causes a significant performance hit to what I would think is the main use case: numerical computations.

This is not a corner case though, nor do I understand why we should preclude non-numerical computations for this functions. And if you want a numerical example, here is one:

>>> import numpy as np
>>> arr = [1, np.iinfo(np.int64).max]
>>> np.apply_along_axis(lambda x : x + 1, 0, arr)
array([2, -9223372036854775808])  # Forced casting to int64, which is wrong

To me, that is just as much a corner case. Anyway, I guess we'll have to agree to disagree, and it is up to the numpy maintainers to decide whether the performance hit is worth it. (Unsubscribing).

Is there any situation this would reasonably arise other than for strings/unicode? Perhaps only those should be special-cased?

That is a plausible idea. I did some testing with different numerical types, and everything seems to be okay. It just seems to the case with strictly strings/unicode.

I agree with @mhvk that such performance degradation for the typical case is not acceptable. A sane way to handle this would be to explicitly cast to the desired dtype in func1d

@shoyer : That's not the issue though. The problem is that we create the output array based on the dtype of the result of the first element. So in the case of strings, numpy creates an array of dtype "<U19" for example. Then, if a string of length-21 is added, numpy simply truncates this element before adding it.

This is why I am looking into creating a special case for strings/unicode as @eric-wieser suggested so that the common cases won't take a significant hit.

The problem is that we create the output array based on the dtype of the result of the first element.

If the first element of the result has the correct dtype, this isn't a problem, e.g.,

>>> arr = np.array([1, np.iinfo(np.int64).max], dtype=object)
>>> np.apply_along_axis(lambda x: x + 1, 0, arr)
array([2, 9223372036854775808L], dtype=object)

If the first element of the result has the correct dtype, this isn't a problem

The issue is that if each element of the result is a string, the length of the string is part of the dtype, and it's pretty unlikely that each string element will be the same length.

Another option would be to addng a dtype argument:

b = np.array([[111,111,0,0,0], [111,111,111,111,111]])
np.apply_along_axis(lambda x: " ".join(map(str, x)), 1, b, dtype=object).astype(str)

I'm surprised this doesn't work either:

b = np.array([[111,111,0,0,0], [111,111,111,111,111]])
def make_object_scalar(x):
    val = np.empty((), dtype=object)
    val[()] = x
    return val
np.apply_along_axis(lambda x: make_object_scalar(" ".join(map(str, x))), 1, b).astype(str)
# TypeError: len() of unsized object

@shoyer : Agreed. What do you think about @eric-wieser proposal?

I would prefer to fix handling of functions that return 0d strings arrays, rather than adding a separate dtype argument.

That TypeError has already been fixed here, so the make_object_scalar version should work on master

Rewrote patch to handle special case for str type as first element, which dramatically improved performance as seen below:

In [1]: import numpy as np
In [2]: arr = np.arange(1000000)
In [3]: %timeit np.apply_along_axis(lambda x: x + 1, 0, arr) 
master : 5.93 ms / loop
PR : 5.97 ms / loop

In [4]: %timeit np.apply_along_axis(lambda x: [x + 1, 5], 0, arr)
master : 1.88 ms / loop
PR : 1.90 ms / loop

Could someone take a look at this? The changes seem to be both robust and performant AFAICT.

☔ The latest upstream changes (presumably b1b67d6) made this pull request unmergeable. Please resolve the merge conflicts.

So, right now, this fails:

>>> np.apply_along_axis(lambda r: np.array(np.sum(r)), 0, np.eye(5))
TypeError: len() of unsized object

Where I would expect an output of np.ones(5).

But this also fails:

np.apply_along_axis(np.diag, 0, np.eye(5))

Both of these are the case when the result of f is an array with ndim != 1. I propose we should do at least one of:

• extend this to work for all ndims
• add a keepdims argument for the 0D case to preserve a unit dimension, so that np.apply_along_axis(np.sum, 0, arr, keepdims=True) is the same as np.sum(arr, axis=0, keepdims=True)

(fixed in #8441, applying the first bullet point )

Fun fact: np.ma.apply_along_axis() already does what this PR is proposing:

>>> a = [[111, 111, 0, 0, 0], [111, 111, 111, 111, 111]]
>>> np.apply_along_axis(lambda x: " ".join(map(str, x)), 1, a)
array(['111 111 0 0 0', '111 111 111 1'], 
      dtype='<U13')
>>> np.ma.apply_along_axis(lambda x: " ".join(map(str, x)), 1, a).data
array(['111 111 0 0 0', '111 111 111 111 111'], 
      dtype='<U19')

Less fun fact - the implementation of np.ma.apply_along_axis is totally independant of apply_along_axis, and probably misses most bugfixes

@eric-wieser : I'm going to close this. You seem to be moving in a much better direction with this patch than I am given all of these other issues you have uncovered with your patches.

Thanks for all of your feedback!

@gfyoung: That PR re-implementing apply_along_axis (#8441) is now merged, so you might want to take another look at this