Re: [Numpy-discussion] Any interest in a 'heaviside' ufunc?
On Tue, Feb 3, 2015 at 12:58 PM, Warren Weckesser warren.weckes...@gmail.com wrote: I have an implementation of the Heaviside function as numpy ufunc. Is there any interest in adding this to numpy? The function is simply: 0if x 0 heaviside(x) = 0.5 if x == 0 1if x 0 I don't think there's anything like it in numpy. Wouldn't scipy.special be a better home for it? Jaime Warren ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion -- (\__/) ( O.o) ( ) Este es Conejo. Copia a Conejo en tu firma y ayúdale en sus planes de dominación mundial. ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion
Re: [Numpy-discussion] Any interest in a 'heaviside' ufunc?
Warren Weckesser warren.weckes...@gmail.com wrote: 0if x 0 heaviside(x) = 0.5 if x == 0 1if x 0 This is not correct. The discrete form of the Heaviside step function has the value 1 for x == 0. heaviside = lambda x : 1 - (x 0).astype(int) Sturla ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion
Re: [Numpy-discussion] Any interest in a 'heaviside' ufunc?
On Tue, Feb 3, 2015 at 11:14 PM, Sturla Molden sturla.mol...@gmail.com wrote: Warren Weckesser warren.weckes...@gmail.com wrote: 0if x 0 heaviside(x) = 0.5 if x == 0 1if x 0 This is not correct. The discrete form of the Heaviside step function has the value 1 for x == 0. heaviside = lambda x : 1 - (x 0).astype(int) By discrete form, do you mean discrete time (i.e. a function defined on the integers)? Then I agree, the discrete time unit step function is defined as u(k) = 0 k 0 1 k = 0 for integer k. The domain of the proposed Heaviside function is not discrete; it is defined for arbitrary floating point (real) arguments. In this case, the choice heaviside(0) = 0.5 is a common convention. See for example, * http://mathworld.wolfram.com/HeavisideStepFunction.html * http://www.mathworks.com/help/symbolic/heaviside.html * http://en.wikipedia.org/wiki/Heaviside_step_function, in particular http://en.wikipedia.org/wiki/Heaviside_step_function#Zero_argument Other common conventions are the right-continuous version that you prefer (heavisde(0) = 1), or the left-continuous version (heaviside(0) = 0). We can accommodate the alternatives with an additional argument that sets the value at 0: heaviside(x, zero_value=0.5) Warren Sturla ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion
Re: [Numpy-discussion] Any interest in a 'heaviside' ufunc?
On Wed, Feb 4, 2015 at 12:18 AM, Warren Weckesser warren.weckes...@gmail.com wrote: On Tue, Feb 3, 2015 at 11:14 PM, Sturla Molden sturla.mol...@gmail.com wrote: Warren Weckesser warren.weckes...@gmail.com wrote: 0if x 0 heaviside(x) = 0.5 if x == 0 1if x 0 This is not correct. The discrete form of the Heaviside step function has the value 1 for x == 0. heaviside = lambda x : 1 - (x 0).astype(int) By discrete form, do you mean discrete time (i.e. a function defined on the integers)? Then I agree, the discrete time unit step function is defined as u(k) = 0 k 0 1 k = 0 for integer k. The domain of the proposed Heaviside function is not discrete; it is defined for arbitrary floating point (real) arguments. In this case, the choice heaviside(0) = 0.5 is a common convention. See for example, * http://mathworld.wolfram.com/HeavisideStepFunction.html * http://www.mathworks.com/help/symbolic/heaviside.html * http://en.wikipedia.org/wiki/Heaviside_step_function, in particular http://en.wikipedia.org/wiki/Heaviside_step_function#Zero_argument Other common conventions are the right-continuous version that you prefer (heavisde(0) = 1), or the left-continuous version (heaviside(0) = 0). We can accommodate the alternatives with an additional argument that sets the value at 0: heaviside(x, zero_value=0.5) What's the usecase for a heaviside function? I don't think I have needed one since I was using mathematica or maple. (x 0).astype(...) (x = 0).astype(...) np.sign(x, dtype) look useful enough for most cases, or not? (What I wish numpy had is conditional place that doesn't calculate all the values. (I think there is a helper function in scipy.stats for that)) Josef Warren Sturla ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion
[Numpy-discussion] advanced indexing question
The numpy reference manual, array objects/indexing/advance indexing, says: Advanced indexing always returns a copy of the data (contrast with basic slicing that returns a view). If I run the following code: import numpy as np d=range[2] x=np.arange(36).reshape(3,2,3,2) y=x[:,d,:,d] y+=1 print x x[:,d,:,d]+=1 print x then the first print x shows that x is unchanged as it should be since y was a copy, not a view, but the second print x shows that all the elements of x with 1st index = 3rd index are now 1 bigger. Why did the left side of x[:,d,:,d]+=1 act like a view and not a copy? Thanks, David ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion
Re: [Numpy-discussion] Views of a different dtype
On Tue, Feb 3, 2015 at 1:28 AM, Sebastian Berg sebast...@sipsolutions.net wrote: On Mo, 2015-02-02 at 06:25 -0800, Jaime Fernández del Río wrote: On Sat, Jan 31, 2015 at 1:17 AM, Sebastian Berg sebast...@sipsolutions.net wrote: On Fr, 2015-01-30 at 19:52 -0800, Jaime Fernández del Río wrote: On Thu, Jan 29, 2015 at 8:57 AM, Nathaniel Smith n...@pobox.com wrote: On Thu, Jan 29, 2015 at 12:56 AM, Jaime Fernández del Río jaime.f...@gmail.com wrote: [...] snip Could we make it more like: check to see if the last dimension works. If not, raise an error (and let the user transpose some other dimension there if that's what they wanted)? Or require the user to specify which dimension will absorb the shape change? (If we were doing this from scratch, then it would be tempting to just say that we always add a new dimension at the end with newtype.itemsize / oldtype.itemsize entries, or absorb such a dimension if shrinking. As a bonus, this would always work, regardless of contiguity! Except that when shrinking the last dimension would have to be contiguous, of course.) When we roll @ in and people start working with stacks of matrices, we will probably find ourselves having to create an alias, similar to .T, for .swapaxes(-1, -2). Searching for the smallest stride allows to take views of such arrays, which does not work right now because the array is no longer contiguous globally. That is true, but I agree with Nathaniel at least as far as that I would prefer a user to be able to safely use `view` even he has not even an inkling about what his memory layout is. One option would be an `axis=-1` default (maybe FutureWarn this from `axis=None` which would look at order, see below -- or maybe have axis='A', 'C' and 'F' and default to 'A' for starters). This even now could start creating bugs when enabling relaxed strides :(, because your good old fortran order complex array being viewed as a float one could expand along the wrong axis, and even without such arrays swap order pretty fast when operating on them, which can create impossibly to find bugs, because even a poweruser is likely to forget about such things. Of course you could argue that view is a poweruser feature and a user using it should keep these things in mind Though if you argue that, you can almost just use `np.ndarray` directly ;) -- ok, not really considering how cumbersome it is, but still. I have been giving this some thought, and am willing to concede that my first proposal may have been too ambitious. So even though the knob goes to 11, we can always do things incrementally. I am also wary of adding new keywords when it seems obvious that we do not have the functionality completely figured out, so here's my new proposal: * The objective is that a view of an array that is the result of slicing a contiguous array should be possible, if it remains contiguous (meaning stride == itemsize) along its original contiguous (first or last) dimension. This eliminates axis transposition from the previous proposal, although reversing the axes completely would also work. * To verify this, unless the C contiguous or Fortran contiguous flags are set, we would still need to look at the strides. An array would be C contiguous if, starting from the last stride it is equal to the itemsize, and working backwards every next stride is larger or equal than the product of the previous stride by the previous dimension. dimensions of size 1 would be ignored for these, except for the last one, which would be taken to have stride = itemsize. The Fortran case is of course the same in reverse. * I think the above combined with the current preference of C contiguousness over Fortran, would actually allow the views to always be reversible, which is also a nice thing to have. This eliminates most of the weirdness, but extends
Re: [Numpy-discussion] Views of a different dtype
On Di, 2015-02-03 at 07:18 -0800, Jaime Fernández del Río wrote: snip Do you have a concrete example of what a non (1, 1) array that fails with relaxed strides would look like? If we used, as right now, the array flags as a first choice point, and only if none is set try to determine it from the strides/dimensions information, I fail to imagine any situation where the end result would be worse than now. I don't think that a little bit of predictable surprising in an advanced functionality is too bad. We could start raising on the face of ambiguity, we refuse to guess errors, even for the current behavior you show above, but that is more likely to trip people by not giving them any simple workaround, that it seems to me would be add a .T if all dimensions are 1 in some particular situations. Or are you thinking of something more serious than a shape mismatch when you write about breaking current code? Yes, I am talking only about wrongly shaped results for some fortran order arrays. A (20, 1) fortran order complex array being viewed as float, will with relaxed strides become a (20, 2) array instead of a (40, 1) one. If there are any real loopholes in expanding this functionality, then lets not do it, but we know we have at least one user unsatisfied with the current performance, so I really think it is worth trying. Plus, I'll admit to that, messing around with some of these stuff deep inside the guts of the beast is lots of fun! ;) I do not think there are loopholes with expanding this functionality. I think there have regressions when we put relaxed strides to on, because suddenly the fortran order array might be expanded along a 1-sized axis, because it is also C order. So I wonder if we can fix these regressions and at the same time maybe provide a more intuitive approach then using the memory order blindly - Sebastian Jaime - Sebastian Jaime - Sebastian I guess the main consideration for this is that we may be stuck with stuff b/c of backwards compatibility. Can you maybe say a little bit about what is allowed now, and what constraints that puts on things? E.g. are we already grovelling around in strides and picking random dimensions in some cases? Just to restate it: right now we only allow new views if the array is globally contiguous, so either along the first or last dimension. Jaime -n -- Nathaniel J. Smith Postdoctoral researcher - Informatics - University of Edinburgh http://vorpus.org ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion -- (\__/) ( O.o) ( ) Este es Conejo. Copia a Conejo en tu firma y ayúdale en sus planes de dominación mundial. ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion -- (\__/) ( O.o) ( ) Este es Conejo. Copia a Conejo en tu firma y ayúdale en sus planes de dominación mundial. ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion
[Numpy-discussion] F2PY and multi-dimension arrays - unexpected array size error
Hi I am helping out with a Python and Fortran project. Let me give you some background: * Fortran source: C Bergstrom FCC C User subroutine VUMAT subroutine VUMAT( C Read only - * nblock, ndir, nshr, nstatev, nprops, * stepTime, dt, * props, * density, strainInc, * tempOld, * stressOld, stateOld, enerInternOld, enerInelasOld, C Write only - * stressNew, stateNew, enerInternNew, enerInelasNew ) C include 'vaba_param.inc' C dimension props(nprops), 1 density(nblock), strainInc(nblock,ndir+nshr), 2 tempOld(nblock), 5 stressOld(nblock,ndir+nshr), stateOld(nblock,nstatev), 6 enerInternOld(nblock), enerInelasOld(nblock), 2 stressNew(nblock,ndir+nshr), stateNew(nblock,nstatev), 3 enerInternNew(nblock), enerInelasNew(nblock) * Corresponding .pyf integer :: nblock integer :: ndir integer :: nshr integer :: nstatev integer :: nprops real :: steptime real :: dt real dimension(nprops) :: props real dimension(nblock) :: density real dimension(nblock,ndir+nshr) :: straininc real dimension(nblock) :: tempold real dimension(nblock,ndir+nshr) :: stressold real dimension(nblock,nstatev) :: stateold real dimension(nblock) :: enerinternold real dimension(nblock) :: enerinelasold real dimension(nblock,ndir+nshr),intent(out) :: stressnew real dimension(nblock,nstatev),intent(out) :: statenew real dimension(nblock),intent(out) :: enerinternnew real dimension(nblock),intent(out) :: enerinelasnew * Python source with call of Fortran routine: nblock = 1 ndir = 3 nshr = 3 nstatev = 3 nprops = 11 stepTime = 1 dt = 1 props = np.array([10, 0.5, 1e10, 5, 1e12, 3e-6, 8e-6, 27, 2], float) density = np.array([[7.8e3]], float) strainInc = np.array([[1,-0.5,-0.5,0,0,0]], float) tempOld = np.array([[1]], float) stressOld = np.array([[1,1,1,1,1,1]], float) stateOld = np.array([[1,1,1]], float) enerInternOld = np.array([1], float) enerInelasOld = np.array([1], float) stressNew = np.array([[]], float) stateNew = np.array([[]], float) enerInternNew = np.array([[]], float) enerInelasNew = np.array([[]], float) stressNew, stateNew, enerInternNew, enerInelasNew = vumat(nblock, ndir, nshr, nstatev, nprops, stepTime, dt, props, density, strainInc, tempOld, stressOld, stateOld, enerInternOld, enerInelasOld) When trying to run with Python 2.7 I get: olof@ubuntu:~$ ./demo.py unexpected array size: new_size=4, got array with arr_size=1 Traceback (most recent call last): File ./demo.py, line 33, in module main() File ./demo.py, line 30, in main stressNew, stateNew, enerInternNew, enerInelasNew = vumat(nblock, ndir, nshr, nstatev, nprops, stepTime, dt, props, density, strainInc, tempOld, stressOld, stateOld, enerInternOld, enerInelasOld) VUMAT_Bergstrom_FCC.error: failed in converting 9th argument `stressold' of VUMAT_Bergstrom_FCC.vumat to C/Fortran array Other stuff: * python 2.7.6 * numpy/f2py 1.8.2 * gcc/gfortran 4.8.2 * ubuntu 14.04 LTS 32-bit I have tried to google, read the f2py manual, fortran tutorials etc, but to no avail. I must also admit that my knowledge in python is so-so and fortran even less(!). What is the missing statement/syntax that I can’t get correct? Your humble programmer, Olof ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion
Re: [Numpy-discussion] Views of a different dtype
On Tue, Feb 3, 2015 at 8:59 AM, Sebastian Berg sebast...@sipsolutions.net wrote: On Di, 2015-02-03 at 07:18 -0800, Jaime Fernández del Río wrote: snip Do you have a concrete example of what a non (1, 1) array that fails with relaxed strides would look like? If we used, as right now, the array flags as a first choice point, and only if none is set try to determine it from the strides/dimensions information, I fail to imagine any situation where the end result would be worse than now. I don't think that a little bit of predictable surprising in an advanced functionality is too bad. We could start raising on the face of ambiguity, we refuse to guess errors, even for the current behavior you show above, but that is more likely to trip people by not giving them any simple workaround, that it seems to me would be add a .T if all dimensions are 1 in some particular situations. Or are you thinking of something more serious than a shape mismatch when you write about breaking current code? Yes, I am talking only about wrongly shaped results for some fortran order arrays. A (20, 1) fortran order complex array being viewed as float, will with relaxed strides become a (20, 2) array instead of a (40, 1) one. That is a limitation of the current implementation too, and happens already whenever relaxed strides are in place. Which is the default for 1.10, right? Perhaps giving 'view' an 'order' or 'axis' kwarg could make sense after all? It should probably be more of a hint of what to do (fortran vs c) when in doubt. C would prioritize last axis, F the first, and we could even add a raise option to have it fail if the axis cannot be inferred from the strides and shape. Current behavior would is equivalent to what C would do. Jaime If there are any real loopholes in expanding this functionality, then lets not do it, but we know we have at least one user unsatisfied with the current performance, so I really think it is worth trying. Plus, I'll admit to that, messing around with some of these stuff deep inside the guts of the beast is lots of fun! ;) I do not think there are loopholes with expanding this functionality. I think there have regressions when we put relaxed strides to on, because suddenly the fortran order array might be expanded along a 1-sized axis, because it is also C order. So I wonder if we can fix these regressions and at the same time maybe provide a more intuitive approach then using the memory order blindly - Sebastian Jaime - Sebastian Jaime - Sebastian I guess the main consideration for this is that we may be stuck with stuff b/c of backwards compatibility. Can you maybe say a little bit about what is allowed now, and what constraints that puts on things? E.g. are we already grovelling around in strides and picking random dimensions in some cases? Just to restate it: right now we only allow new views if the array is globally contiguous, so either along the first or last dimension. Jaime -n -- Nathaniel J. Smith Postdoctoral researcher - Informatics - University of Edinburgh http://vorpus.org ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion -- (\__/) ( O.o) ( ) Este es Conejo. Copia a Conejo en tu firma y ayúdale en sus planes de dominación mundial. ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion
[Numpy-discussion] Any interest in a 'heaviside' ufunc?
I have an implementation of the Heaviside function as numpy ufunc. Is there any interest in adding this to numpy? The function is simply: 0if x 0 heaviside(x) = 0.5 if x == 0 1if x 0 Warren ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion
Re: [Numpy-discussion] Any interest in a 'heaviside' ufunc?
That seems useful to me. On Tue, Feb 3, 2015 at 3:58 PM, Warren Weckesser warren.weckes...@gmail.com wrote: I have an implementation of the Heaviside function as numpy ufunc. Is there any interest in adding this to numpy? The function is simply: 0if x 0 heaviside(x) = 0.5 if x == 0 1if x 0 Warren ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion
Re: [Numpy-discussion] Views of a different dtype
On Mo, 2015-02-02 at 06:25 -0800, Jaime Fernández del Río wrote: On Sat, Jan 31, 2015 at 1:17 AM, Sebastian Berg sebast...@sipsolutions.net wrote: On Fr, 2015-01-30 at 19:52 -0800, Jaime Fernández del Río wrote: On Thu, Jan 29, 2015 at 8:57 AM, Nathaniel Smith n...@pobox.com wrote: On Thu, Jan 29, 2015 at 12:56 AM, Jaime Fernández del Río jaime.f...@gmail.com wrote: [...] snip Could we make it more like: check to see if the last dimension works. If not, raise an error (and let the user transpose some other dimension there if that's what they wanted)? Or require the user to specify which dimension will absorb the shape change? (If we were doing this from scratch, then it would be tempting to just say that we always add a new dimension at the end with newtype.itemsize / oldtype.itemsize entries, or absorb such a dimension if shrinking. As a bonus, this would always work, regardless of contiguity! Except that when shrinking the last dimension would have to be contiguous, of course.) When we roll @ in and people start working with stacks of matrices, we will probably find ourselves having to create an alias, similar to .T, for .swapaxes(-1, -2). Searching for the smallest stride allows to take views of such arrays, which does not work right now because the array is no longer contiguous globally. That is true, but I agree with Nathaniel at least as far as that I would prefer a user to be able to safely use `view` even he has not even an inkling about what his memory layout is. One option would be an `axis=-1` default (maybe FutureWarn this from `axis=None` which would look at order, see below -- or maybe have axis='A', 'C' and 'F' and default to 'A' for starters). This even now could start creating bugs when enabling relaxed strides :(, because your good old fortran order complex array being viewed as a float one could expand along the wrong axis, and even without such arrays swap order pretty fast when operating on them, which can create impossibly to find bugs, because even a poweruser is likely to forget about such things. Of course you could argue that view is a poweruser feature and a user using it should keep these things in mind Though if you argue that, you can almost just use `np.ndarray` directly ;) -- ok, not really considering how cumbersome it is, but still. I have been giving this some thought, and am willing to concede that my first proposal may have been too ambitious. So even though the knob goes to 11, we can always do things incrementally. I am also wary of adding new keywords when it seems obvious that we do not have the functionality completely figured out, so here's my new proposal: * The objective is that a view of an array that is the result of slicing a contiguous array should be possible, if it remains contiguous (meaning stride == itemsize) along its original contiguous (first or last) dimension. This eliminates axis transposition from the previous proposal, although reversing the axes completely would also work. * To verify this, unless the C contiguous or Fortran contiguous flags are set, we would still need to look at the strides. An array would be C contiguous if, starting from the last stride it is equal to the itemsize, and working backwards every next stride is larger or equal than the product of the previous stride by the previous dimension. dimensions of size 1 would be ignored for these, except for the last one, which would be taken to have stride = itemsize. The Fortran case is of course the same in reverse. * I think the above combined with the current preference of C contiguousness over Fortran, would actually allow the views to always be reversible, which is also a nice thing to have. This eliminates most of the weirdness, but extends current functionality to cover cases like Jens reported a few days back. Does this sound better? It seems fine as such, but I still worry
Re: [Numpy-discussion] Views of a different dtype
On Tue Feb 03 2015 at 1:47:34 PM Jaime Fernández del Río jaime.f...@gmail.com wrote: On Tue, Feb 3, 2015 at 8:59 AM, Sebastian Berg sebast...@sipsolutions.net wrote: On Di, 2015-02-03 at 07:18 -0800, Jaime Fernández del Río wrote: snip Do you have a concrete example of what a non (1, 1) array that fails with relaxed strides would look like? If we used, as right now, the array flags as a first choice point, and only if none is set try to determine it from the strides/dimensions information, I fail to imagine any situation where the end result would be worse than now. I don't think that a little bit of predictable surprising in an advanced functionality is too bad. We could start raising on the face of ambiguity, we refuse to guess errors, even for the current behavior you show above, but that is more likely to trip people by not giving them any simple workaround, that it seems to me would be add a .T if all dimensions are 1 in some particular situations. Or are you thinking of something more serious than a shape mismatch when you write about breaking current code? Yes, I am talking only about wrongly shaped results for some fortran order arrays. A (20, 1) fortran order complex array being viewed as float, will with relaxed strides become a (20, 2) array instead of a (40, 1) one. That is a limitation of the current implementation too, and happens already whenever relaxed strides are in place. Which is the default for 1.10, right? Perhaps giving 'view' an 'order' or 'axis' kwarg could make sense after all? It should probably be more of a hint of what to do (fortran vs c) when in doubt. C would prioritize last axis, F the first, and we could even add a raise option to have it fail if the axis cannot be inferred from the strides and shape. Current behavior would is equivalent to what C would do. Jaime IMHO, the best option would be something like this: - When changing to a type with smaller itemsize, add a new axis after the others so the resulting array is C contiguous (unless a different axis is specified by a keyword argument). The advantage here is that if you index the new view using the old indices for an entry, you get an array showing its representation in the new type. - No shape change for views with the same itemsize - When changing to a type with a larger itemsize, collapse along the last axis unless a different axis is specified, throwing an error if the axis specified does not match the axis specified. The last point essentially is just adding an axis argument. I like that idea because it gives users the ability to do all that the array's memory layout allows in a clear and concise way. Throwing an error if the default axis doesn't work would be a good way to prevent strange bugs from happening when the default behavior is expected. The first point would be a break in backwards compatibility, so I'm not sure if it's feasible at this point. The advantage would be that all all arrays returned when using this functionality would be contiguous along the last axis. The shape of the new array would be independent of the memory layout of the original one. This would also be a much cleaner way to ensure that views of a different type are always reversible while still allowing for relaxed strides. Either way, thanks for looking into this. It's a great feature to have available. -Ian Henriksen If there are any real loopholes in expanding this functionality, then lets not do it, but we know we have at least one user unsatisfied with the current performance, so I really think it is worth trying. Plus, I'll admit to that, messing around with some of these stuff deep inside the guts of the beast is lots of fun! ;) I do not think there are loopholes with expanding this functionality. I think there have regressions when we put relaxed strides to on, because suddenly the fortran order array might be expanded along a 1-sized axis, because it is also C order. So I wonder if we can fix these regressions and at the same time maybe provide a more intuitive approach then using the memory order blindly - Sebastian Jaime - Sebastian Jaime - Sebastian I guess the main consideration for this is that we may be stuck with stuff b/c of backwards compatibility. Can you maybe say a little bit about what is allowed now, and what constraints that puts on things? E.g. are we already grovelling around in strides and picking random
