Fortran order arrays to and from numpy arrays
Hi, I'm currently puzzling over how to best convert (column major order) matlab arrays to numpy arrays and vice versa -- I'm looking for a solution that's simple, general and reasonably fast -- being also applicable to Numeric arrays would be a plus (I'd like to retain Numeric compatibility for the current version of the code). The solution that mlabwrap (http://mlabwrap.sf.net) currently employs to transfer arrays between matlab and python is to have a low level C++ extension (mlabraw) copying the data into a new numpy array created with Pyarray_SimpleNew and using a simple or nested for-loop to copy the data around (only rank up to 2 is currently handled), e.g. //matlab matrix -> numpy for (mwIndex lCol = 0; lCol < lCols; lCol++) { double *lDst = (double *)(lRetval->data) + lCol; for (mwIndex lRow = 0; lRow < lRows; lRow++, lDst += lCols) { *lDst = *lPR++;}} //numpy 2d array -> matlab npy_intp lRowDelta = pStrides[0]/sizeof(T); npy_intp lColDelta = pStrides[1]/sizeof(T); for (npy_intp lCol=0; lCol < pCols; lCol++) { T *lSrc = pSrc + lCol*lColDelta; for (npy_intp lRow=0; lRow < pRows; lRow++, lSrc += lRowDelta) { *pDst++ = *lSrc;}} Some simple timing results suggest 2 things: 1. There is an enormous overhead associated with calling into matlab (e.g. just using mlabraw to evaluate ``sin(1);`` in matlab, without accessing the result the result is about 1000 times slower than ``math.sin(1)``. I have no idea where this rather huge overhead comes from, but I suspect there isn't much one can do about it. 2. Despite this overhead, copying around large arrays (e.g. >=1e5 elements) in above way causes notable additional overhead. Whilst I don't think there's a sane way to avoid copying by sharing data between numpy and matlab the copying could likely be done better. I guess the moral from 1 might be that ctypes will be the way to go for future version, I would assume that the overhead associated with ctypes doesn't matter that much and that it will be much less painful to maintain than having a C or C++ extension module, especially given that the Mathworks seem fond of changing the C API every minor release. I haven't used ctypes yet so I'd be interested to know what people think. The moral from 2. might be just to memcpy the data wholesale where possible and transpose and reshape afterwards or, in numpy, to set the fortran-order flag (and query fortran-contiguousness on copying to matlab) -- I suppose that specifying fortran order and tranposing/reshaping are pretty much equivalent? So matlab to numpy would look roughly so: a_fortran_order = memcpy_raw_matlab_data_into_numpy_array(ptr) a = a_fortran_order.T.reshape(a_fortran_order.shape) I would assume that this ought be on average more efficient than always transforming to C-contiguous (as some of the time the returned array will anyway be transposed afterwards or accessed in a fashion were column major is faster; the case that several operations that profit from row-major storage are carried out on the same array seems unlikely). Unfortunately I don't see an easy way to use the same approach the other way (matlab doesn't seem to offer much on the C level to manipulate arrays), so I'd presumably need something like: stuff_into_matlab_array(a.T.reshape(a.shape).copy()) the question is how to avoid doing two copies. Any comments appreciated, 'as
On 23 Feb 2007 21:19:05 +0000, Alexander Schmolck <a.schmolck@gmx.net> wrote:
Hi,
I'm currently puzzling over how to best convert (column major order) matlab arrays to numpy arrays and vice versa -- I'm looking for a solution that's simple, general and reasonably fast -- being also applicable to Numeric arrays would be a plus (I'd like to retain Numeric compatibility for the current version of the code).
The solution that mlabwrap (http://mlabwrap.sf.net) currently employs to transfer arrays between matlab and python is to have a low level C++ extension (mlabraw) copying the data into a new numpy array created with Pyarray_SimpleNew and using a simple or nested for-loop to copy the data around (only rank up to 2 is currently handled), e.g.
//matlab matrix -> numpy for (mwIndex lCol = 0; lCol < lCols; lCol++) { double *lDst = (double *)(lRetval->data) + lCol; for (mwIndex lRow = 0; lRow < lRows; lRow++, lDst += lCols) { *lDst = *lPR++;}}
//numpy 2d array -> matlab npy_intp lRowDelta = pStrides[0]/sizeof(T); npy_intp lColDelta = pStrides[1]/sizeof(T); for (npy_intp lCol=0; lCol < pCols; lCol++) { T *lSrc = pSrc + lCol*lColDelta; for (npy_intp lRow=0; lRow < pRows; lRow++, lSrc += lRowDelta) { *pDst++ = *lSrc;}}
Some simple timing results suggest 2 things:
1. There is an enormous overhead associated with calling into matlab (e.g. just using mlabraw to evaluate ``sin(1);`` in matlab, without accessing the result the result is about 1000 times slower than ``math.sin(1)``. I have no idea where this rather huge overhead comes from, but I suspect there isn't much one can do about it.
2. Despite this overhead, copying around large arrays (e.g. >=1e5 elements) in above way causes notable additional overhead. Whilst I don't think there's a sane way to avoid copying by sharing data between numpy and matlab the copying could likely be done better.
I guess the moral from 1 might be that ctypes will be the way to go for future version, I would assume that the overhead associated with ctypes doesn't matter that much and that it will be much less painful to maintain than having a C or C++ extension module, especially given that the Mathworks seem fond of changing the C API every minor release. I haven't used ctypes yet so I'd be interested to know what people think.
The moral from 2. might be just to memcpy the data wholesale where possible and transpose and reshape afterwards or, in numpy, to set the fortran-order flag (and query fortran-contiguousness on copying to matlab) -- I suppose that specifying fortran order and tranposing/reshaping are pretty much equivalent? So matlab to numpy would look roughly so:
a_fortran_order = memcpy_raw_matlab_data_into_numpy_array(ptr) a = a_fortran_order.T.reshape(a_fortran_order.shape)
I would assume that this ought be on average more efficient than always transforming to C-contiguous (as some of the time the returned array will anyway be transposed afterwards or accessed in a fashion were column major is faster; the case that several operations that profit from row-major storage are carried out on the same array seems unlikely).
Unfortunately I don't see an easy way to use the same approach the other way (matlab doesn't seem to offer much on the C level to manipulate arrays), so I'd presumably need something like:
stuff_into_matlab_array(a.T.reshape(a.shape).copy())
the question is how to avoid doing two copies.
Any comments appreciated,
The easiest way to deal with the ordering is to use the order keyword in numpy: In [4]: a = array([0,1,2,3]).reshape((2,2), order='F') In [5]: a Out[5]: array([[0, 2], [1, 3]]) You would still need to get access to something to reshape, shared memory or something, but the key is that you don't have to reorder the elements, you just need the correct strides and offsets to address the elements in Fortran order. I have no idea if this works in numeric. Chuck
"Charles R Harris" <charlesr.harris@gmail.com> writes:
Unfortunately I don't see an easy way to use the same approach the other way (matlab doesn't seem to offer much on the C level to manipulate arrays), so I'd presumably need something like:
stuff_into_matlab_array(a.T.reshape(a.shape).copy())
the question is how to avoid doing two copies.
Any comments appreciated,
The easiest way to deal with the ordering is to use the order keyword in numpy:
In [4]: a = array([0,1,2,3]).reshape((2,2), order='F')
In [5]: a Out[5]: array([[0, 2], [1, 3]])
You would still need to get access to something to reshape, shared memory or something, but the key is that you don't have to reorder the elements, you just need the correct strides and offsets to address the elements in Fortran order. I have no idea if this works in numeric.
It doesn't work in Numeric, but that isn't much of any issue because I think it ought to be pretty much equivalent by transposing and reshaping. However the problem is that I *do* need to reorder the elements for numpy->matlab and I'm not sure how to best do this (without unnecessary copying and temporary numpy array creation but using numpy functionality if possible). 'as
On 25 Feb 2007 01:44:01 +0000, Alexander Schmolck <a.schmolck@gmx.net> wrote:
"Charles R Harris" <charlesr.harris@gmail.com> writes:
Unfortunately I don't see an easy way to use the same approach the other way (matlab doesn't seem to offer much on the C level to manipulate arrays), so I'd presumably need something like:
stuff_into_matlab_array(a.T.reshape(a.shape).copy())
the question is how to avoid doing two copies.
Any comments appreciated,
The easiest way to deal with the ordering is to use the order keyword in numpy:
In [4]: a = array([0,1,2,3]).reshape((2,2), order='F')
In [5]: a Out[5]: array([[0, 2], [1, 3]])
You would still need to get access to something to reshape, shared memory or something, but the key is that you don't have to reorder the elements, you just need the correct strides and offsets to address the elements in Fortran order. I have no idea if this works in numeric.
It doesn't work in Numeric, but that isn't much of any issue because I think it ought to be pretty much equivalent by transposing and reshaping. However the problem is that I *do* need to reorder the elements for numpy->matlab and I'm not sure how to best do this (without unnecessary copying and temporary numpy array creation but using numpy functionality if possible).
I don't see any way to get around a copy, but you can make numpy do the work. For example: In [12]: a = array([[0,1],[2,3]]) In [13]: b = array(a, order='f') In [14]: a.flags Out[14]: C_CONTIGUOUS : True F_CONTIGUOUS : False OWNDATA : True WRITEABLE : True ALIGNED : True UPDATEIFCOPY : False In [15]: b.flags Out[15]: C_CONTIGUOUS : False F_CONTIGUOUS : True OWNDATA : True WRITEABLE : True ALIGNED : True UPDATEIFCOPY : False F_CONTIGUOUS is what you want. The trick is to somehow use memory in the construction of the reordered array that is already designated for matlab. I don't know how to do this, but I think it might be doable. Travis is your best bet to answer that question. Chuck
Alexander Schmolck wrote:
2. Despite this overhead, copying around large arrays (e.g. >=1e5 elements) in above way causes notable additional overhead. Whilst I don't think there's a sane way to avoid copying by sharing data between numpy and matlab the copying could likely be done better.
Alex, what do you think about "hybrid arrays"? http://www.mail-archive.com/numpy-discussion@lists.sourceforge.net/msg03748....
Andrew Straw <strawman@astraw.com> writes:
Alexander Schmolck wrote:
2. Despite this overhead, copying around large arrays (e.g. >=1e5 elements) in above way causes notable additional overhead. Whilst I don't think there's a sane way to avoid copying by sharing data between numpy and matlab the copying could likely be done better.
Alex, what do you think about "hybrid arrays"?
http://www.mail-archive.com/numpy-discussion@lists.sourceforge.net/msg03748....
Oh, that's why I said no *sane* way :) I read about hybrid arrays, but as far as I can tell the only way to use them to avoid copying stuff around is to create your own hybrid array memory pool (as you suggested downthread), which seems a highly unattractive pain-for-gain trade-off, especially given that you *do* have to reorder the data as you go from numpy to matlab -- unless of course I'm missing something. I have the impression that just memcpy'ing a single chunk of data isn't that much of a performance sink, but the reordering copying that mlabwrap currently does seems rather expensive. In other words, I think going from matlab to numpy is fine (just memcpy into a newly created fortran-order numpy array, or more or less equivalently, memcpy into a C-order array, transpose and reshape), the question appears to be how to best go from numpy to matlab when the numpy array isn't fortran-contiguous. I assume that it makes more sense to rely on some numpy functionality than to use a custom reordering-copy routine, especially if I want to move to ctypes later. Is there anything better than 1. allocating a matlab array 2. transposing and reshaping the numpy array 3. allocating (or keeping around) a temporary numpy array with data pointing to the matlab array data 4. using some function (PyArray_CopyInto?) to copy from the transposed, reshaped numpy array into the temporary numpy array thereby filling the matlab array with an appropriately reordered copy of the original array ? cheers, alex
participants (3)
-
Alexander Schmolck -
Andrew Straw -
Charles R Harris