is there a C-API function for numpy which implements Python's
multidimensional indexing? Say, I have a 2d-array
PyArrayObject * M;
and an index
how do I extract the i-th row or column M[i,:] respectively M[:,i]?
I am looking for a function which gives again a PyArrayObject * and
which is a view to M (no copied data; the result should be another
PyArrayObject whose data and strides points to the correct memory
portion of M).
I searched the API documentation, Google and mailing lists for quite a
long time but didn't find anything. Can you help me?
quite often I work with block matrices. Matlab offers the convenient notation
[ a b; c d ]
to stack matrices. The numpy equivalent is kinda clumsy:
I wrote the little function `stack` that does exactly that:
stack([[a, b], [c, d]])
In my case `stack` replaced `hstack` and `vstack` almost completely.
If you're interested in including it in numpy I created a pull request
. I'm looking forward to getting some feedback!
I have just sent a PR (https://github.com/numpy/numpy/pull/5015), adding
the possibility of having frozen dimensions in gufunc signatures. As a
proof of concept, I have added a `cross1d` gufunc to
In : import numpy as np
In : from numpy.core.umath_tests import cross1d
In : cross1d.signature
In : a = np.random.rand(1000, 3)
In : b = np.random.rand(1000, 3)
In : np.allclose(np.cross(a, b), cross1d(a, b))
In : %timeit np.cross(a, b)
10000 loops, best of 3: 76.1 us per loop
In : %timeit cross1d(a, b)
100000 loops, best of 3: 13.1 us per loop
In : c = np.random.rand(1000, 2)
In : d = np.random.rand(1000, 2)
In : cross1d(c, d)
ValueError Traceback (most recent call last)
<ipython-input-11-72c66212e40c> in <module>()
----> 1 cross1d(c, d)
ValueError: cross1d: Operand 0 has a mismatch in its core dimension 0, with
gufunc signature (3),(3)->(3) (size 2 is different from 3)
The speed up over `np.cross` is nice, and while `np.cross` is not the best
of examples, as it needs to handle more sizes, in many cases this will
allow producing gufuncs that work without a Python wrapper redoing checks
that are best left to the iterator, such as dimension sizes.
It still needs tests, but before embarking on fully developing those, I
wanted to make sure that there is an interest on this.
I would also like to further enhance gufuncs providing computed dimensions,
e.g. making it possible to e.g. define `pairwise_cross` with signature '(n,
3)->($m, 3)', where the $ indicates that m is a computed dimension, that
would have to be calculated by a function passed to the gufunc constructor
and stored in the gufunc object, based on the other core dimensions. In
this case it would make $m be n*(n-1), so that all pairwise cross products
between 3D vectors could be computed.
The syntax with '$' is kind of crappy, so any suggestions on how to better
express this in the signature are more than welcome, as well as any
feedback on the merits (or lack of them) of implementing this.
( > <) Este es Conejo. Copia a Conejo en tu firma y ayúdale en sus planes
de dominación mundial.
I am on Windows 8.1, 64-bit, using the usual version of Anaconda-64-bit for
I just come from the mailing list of matplotlib: my PC is not able to do 3D
plot with matplotlib.
Python crashes completely if I try to plot something in 3D and to show it:
"Python.exe has stopped working..."
The main problem in resolving my issue is that the python traceback and the
faulthandler package do not succeed in catching the errors.
Only the Windows crash report allows us to see what happened before the
The windows report tells that the file _dotblas.pyd from numpy is likely
responsible for the crash.
So I tried a simple:
import numpy; numpy.test()
numpy.test() does not work! I dont get any traceback from python. It just
crashes as before: "Python.exe has stopped working..."
And again, the Windows report blames _dotblas.pyd :)
Here the new Windows report for the numpy.test(), I highlighted the lines
which blame _dotblas.pyd:
*Sig.Name=Fault Module NameSig.Value=_dotblas.pyd*
Sig.Name=Fault Module Version
Sig.Name=Fault Module Timestamp
DynamicSig.Name=Additional Information 1
DynamicSig.Name=Additional Information 2
DynamicSig.Name=Additional Information 3
DynamicSig.Name=Additional Information 4
UI=python.exe has stopped working
UI=Windows can check online for a solution to the problem.
UI=Check online for a solution and close the program
UI=Check online for a solution later and close the program
UI=Close the program
Any idea, of what I can do?
I try to reinstall anaconda,
I also try conda remove numpy and conda install numpy.
That does not seem to change anything.
I'm a novice with respect to scientific computing using python. I've read
that numpy.array isn't arranged according to the 'right-hand-rule'
(right-hand-rule => thumb = +x; index finger = +y, bend middle finder =
+z). This is also confirmed by an old message I dug up from the mailing
list archives. (see message below)
I guess this has consequences for certain algorithms (been a while, should
actually revise some algebra textbooks). At least it does when I try to
mentally visualize a 3D array when I'm not sure which dimensions to use...
What are the consequences in using 3D arrays in for example transformation
algorithms or other algorithms which expect certain shape? Should I always
'reshape' before using them or do most algorithms take this into account in
the underlying algebra?
Does the 'Fortran contiguous' shape always respect the 'right-hand-rule'
(for 3D arrays)? And just to be sure: Is the information telling if an
array is C-contiguous or Fortran-contiguous always stored within the array?
/ Old message
From: Anne Archibald <peridot.faceted <at> gmail.com>
Subject: Re: dimension aligment
Date: 2008-05-20 18:04:46 GMT (6 years, 35 weeks, 5 days, 4 hours and 34
2008/5/20 Thomas Hrabe <thrabe <at> burnham.org>:
> given a *3d* *array*
> a =
> returns (4,2,3)
> so I assume the first digit is the 3rd dimension, second is 2nd dim and
> third is the first.
> how is the data aligned in memory now?
> according to the strides it should be
> if I had an *array* of more dimensions, the first digit returned by shape
> should always be the highest dim.
You are basically *right*, but this is a surprisingly subtle issue for numpy.
A numpy *array* is basically a block of memory and some description. One
piece of that description is the type of data it contains (i.e., how
to interpret each chunk of memory) for example int32, float64, etc.
Another is the sizes of all the various dimensions. A third piece,
which makes many of the things numpy does possible, is the "strides".
The way numpy works is that basically it translates
into a lookup of the item in the memory block at position
This means, if you have an *array* A and you want every second element
(A[::2]), all numpy needs to do is *hand* you back a new *array* pointing
to the same data block, but with strides doubled. Similarly if you
want to transpose a two-dimensional *array*, all it needs to do is
exchange strides and strides; no data need be moved.
This means, though, that if you are *hand*ed a numpy *array*, the elements
can be arranged in memory in quite a complicated fashion. Sometimes
this is no problem - you can always use the strides to find it all.
But sometimes you need the data arranged in a particular way. numpy
defines two particular ways: "C contiguous" and "FORTRAN contiguous".
"C contiguous" *array*s are what you describe, and they're what numpy
produces by default; they are arranged so that the *right*most index has
the smallest stride. "FORTRAN contiguous" *array*s are arranged the
other way around; the leftmost index has the smallest stride. (This is
how FORTRAN *array*s are arranged in memory.)
There is also a special case: the reshape() function changes the shape
of the *array*. It has an "order" argument that describes not how the
elements are arranged in memory but how you want to think of the
elements as arranged in memory for the reshape operation.
I took time to create mingw-w64 based wheels of numpy-1.9.1 and
scipy-0.15.1 source distributions and put them on
https://bitbucket.org/carlkl/mingw-w64-for-python/downloads as well as on
binstar.org. The test matrix is python-2.7 and 3.4 for both 32bit and
Feedback is welcome.
The wheels can be pip installed with:
pip install -i https://pypi.binstar.org/carlkl/simple numpy
pip install -i https://pypi.binstar.org/carlkl/simple scipy
Some technical details: the binaries are build upon OpenBLAS as accelerated
BLAS/Lapack. OpenBLAS itself is build with dynamic kernels (similar to MKL)
and automatic runtime selection depending on the CPU. The minimal requested
feature supplied by the CPU is SSE2. SSE1 and non-SSE CPUs are not
supported with this builds. This is the default for 64bit binaries anyway.
OpenBLAS is deployed as part of the numpy wheel. That said, the scipy
wheels mentioned above are dependant on the installation of the OpenBLAS
based numpy and won't work i.e. with an installed numpy-MKL.
For the numpy 32bit builds there are 3 failures for special FP value tests,
due to a bug in mingw-w64 that is still present. All scipy versions show up
7 failures with some numerical noise, that could be ignored (or corrected
with relaxed asserts in the test code).
PR's for numpy and scipy are in preparation. The mingw-w64 compiler used
for building can be found at
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
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
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
- 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)
(3, 1, 2)
>>> a.strides # the stride of the size 1 dimension could be anything,
(32, 8, 8)
b = a.view(complex) # this fails right now, but should work
(3, 1, 1)
>>> b.strides # the stride of the size 1 dimensions could be anything,
(32, 16, 16)
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.
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
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:
1. That the newtype and oldtype be compatible, as defined by the
algorithm checking object offsets linked above.
2. If newtype.itemsize == oldtype.itemsize no more checks are needed,
make it happen!
3. 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.
4. For non-contiguous arrays:
1. 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!
2. Ignoring dimensions of size 1, find a contiguous dimension, i.e.
stride == oldtype.itemsize
1. 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
evenly divided by newitem,itemsize.
2. 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
( > <) Este es Conejo. Copia a Conejo en tu firma y ayúdale en sus planes
de dominación mundial.
I tried to use numpy (version 1.9.1, installed by `pip install numpy`)
on cygwin64. I encountered the following weird bug:
>>> import numpy
>>> with numpy.errstate(all='raise'):
... print 1/float64(0.0)
I was expecting a FloatingPointError, but it didn't show up. Curiously,
with different numerical types (all intxx, or float128), I indeed get
Same thing with the most recent git version, or with 1.7.1 provided as a
precompiled package by cygwin. This behavior does not happen on cygwin32
(I always get the FloatingPointError there).
I wonder if there is something weird with my config, or if this is a
genuine reproducible bug. If so, where should I start looking if I want
to fix it? (I don't know anything about numpy's code)