HI all,
There has been some recent discussion going on on the limitations that numpy imposes to taking views of an array with a different dtype.
As of right now, you can basically only take a view of an array if it has no Python objects and neither the old nor the new dtype are structured. Furthermore, the array has to be either C or Fortran contiguous.
This seem to be way too strict, but the potential for disaster getting a loosening of the restrictions wrong is big, so it should be handled with care.
Allan Haldane and myself have been looking into this separately and discussing some of the details over at github, and we both think that the only true limitation that has to be imposed is that the offsets of Python objects within the new and old dtypes remain compatible. I have expanded Allan's work from here:
to make it as flexible as I have been able. An implementation of the algorithm in Python, with a few tests, can be found here:
I would appreciate getting some eyes on it for correctness, and to make sure that it won't break with some weird dtype.
I am also trying to figure out what the ground rules for stride and shape conversions when taking a view with a different dtype should be. I submitted a PR (gh-5508) a couple for days ago working on that, but I am not so sure that the logic is completely sound. Again, to get more eyes on it, I am going to reproduce my thoughts here on the hope of getting some feedback.
The objective would be to always allow a view of a different dtype (given that the dtypes be compatible as described above) to be taken if:
- The itemsize of the dtype doesn't change.
- The itemsize changes, but the array being viewed is the result of slicing and transposing axes of a contiguous array, and it is still contiguous, defined as stride == dtype.itemsize, along its smallest-strided dimension, and the itemsize of the newtype exactly divides the size of that dimension.
- Ideally taking a view should be a reversible process, i.e. if oldtype = arr.dtype, then doing arr.view(newtype).view(oldtype) should give you back a view of arr with the same original shape, strides and dtype.
This last point can get tricky if the minimal stride dimension has size 1, as there could be several of those, e.g.:
>>> a = np.ones((3, 4, 1), dtype=float)[:, :2, :].transpose(0, 2, 1)
>>> a.strides # the stride of the size 1 dimension could be anything, ignore it!
b = a.view(complex) # this fails right now, but should work
>>> b.strides # the stride of the size 1 dimensions could be anything, ignore them!
c = b.view(float) # which of the two size 1 dimensions should we expand?
"In the face of ambiguity refuse the temptation to guess" dictates that last view should raise an error, unless we agree and document some default. Any thoughts?
Then there is the endless complication one could get into with arrays created with as_strided. I'm not smart enough to figure when and when not those could work, but am willing to retake the discussion if someone wiser si interested.
With all these in mind, my proposal for the new behavior is that taking a view of an array with a different dtype would require:
- That the newtype and oldtype be compatible, as defined by the algorithm checking object offsets linked above.
- If newtype.itemsize == oldtype.itemsize no more checks are needed, make it happen!
- If the array is C/Fortran contiguous, check that the size in bytes of the last/first dimension is evenly divided by newtype.itemsize. If it does, go for it.
- For non-contiguous arrays:
- Ignoring dimensions of size 1, check that no stride is smaller than either oldtype.itemsize or newtype.itemsize. If any is found this is an as_strided product, sorry, can't do it!
- Ignoring dimensions of size 1, find a contiguous dimension, i.e. stride == oldtype.itemsize
- If found, check that it is the only one with that stride, that it is the minimal stride, and that the size in bytes of that dimension is evenly divided by newitem,itemsize.
- If none is found, check if there is a size 1 dimension that is also unique (unless we agree on a default, as mentioned above) and that newtype.itemsize evenly divides oldtype.itemsize.
Apologies for the long, dense content, but any thought or comments are very welcome.
Jaime
--
(\__/)
( O.o)
( > <) Este es Conejo. Copia a Conejo en tu firma y ayúdale en sus planes de dominación mundial.