Re: Advice on converting Numarray C extension?

_______________________________________________
Numpy-discussion mailing list
Numpy-discussion <at> scipy.org
http://mail.scipy.org/mailman/listinfo/numpy-discussion

Re: Tutorial topics for SciPy'09 Conference

Hi,

On Mon, Jun 1, 2009 at 10:20 PM, Fernando Perez<fperez.net <at> gmail.com> wrote:
> The time for the Scipy'09 conference is rapidly approaching, and we
> would like to both announce the plan for tutorials and solicit
> feedback from everyone on topics of interest.

rather than rehash much here, where it's not easy to paste a table,
I've posted a note with the poll results here:

http://fdoperez.blogspot.com/2009/06/scipy-advanced-tutorials-results.html

The short and plain-text-friendly version is the final topic ranking:

1	Advanced topics in matplotlib use
2	Advanced numpy
3	Designing scientific interfaces with Traits
4	Mayavi/TVTK
5	Cython
6	Symbolic computing with sympy
7	Statistics with Scipy
8	Using GPUs with PyCUDA
9	Testing strategies for scientific codes
10	Parallel computing in Python and mpi4py
11	Sparse Linear Algebra with Scipy
12	Structured and record arrays in numpy
13	Design patterns for efficient iterator-based scientific codes
14	Sage
15	The TimeSeries scikit
16	Hermes: high order Finite Element Methods
17	Graph theory with NetworkX

We're currently contacting speakers, and we'll let you know once a
final list is made with confirmed speakers.

Cheers,

f

np.memmap and memory usage

	Hello,

	I'm using numpy.memmap to open big 3-D arrays of Xray tomography
data. After I have created a new array using memmap, I modify the
contrast of every Z-slice (along the first dimension) inside a for loop,
for a better visualization of the data. Although I call memmap.flush
after each modification of a Z-slice, the memory used by Ipython keeps
increasing at every new iteration. At the end of the loop, the memory
used by Ipython is of the order of magnitude of the size of the data file
(1.8Go !). I would have expected that the maximum amount of memory used
would corresponde to only one Z-slice of the 3D array. See the code
snapshots below for more details.

	Is this an expected behaviour? How can I reduce the amount of
memory used by Ipython and still process my data?

import numpy as np
imin, imax = -2, 2
data = np.memmap(filename = 'data.vol', mode='r+', dtype=float32,
shape=(512,981, 981))
for sl in data:
    sl[sl<imin] = imin
    sl[sl>imax] = imax
    data.flush() # tried with and without this line

	Thanks in advance,

	Emmanuelle

Numpy complex types, packing and C99

Hi,

    I would like to add an explicit configuration test to check that our
complex type is compatible with the C99 complex type (when available).
Is this ok ?

    As currently defined, complex c types (npy_cfloat, etc...) are not
defined in a way such as they are binary compatible with the C99 complex
type. Strictly speaking, packing the complex type is an ABI break, but
we already make the assumption in ufunc, so we would have wrong
results/crashes currently if it were not packed, so I believe the check
is harmless by itself. The rationale is that I would like to add support
for complex math in npy_math (cabs, cexp, etc...). As I would like to
use the C99 complex type (when available), this requires that numpy
complex type is ABI compatible with C99 type.

cheers,

David

Re: permutation symbol

On Tue, 30 Jun 2009 13:51:15 -0600
  Charles R Harris <charlesr.harris <at> gmail.com> wrote:
> On Tue, Jun 30, 2009 at 12:26 PM, Nils Wagner
> <nwagner <at> iam.uni-stuttgart.de>wrote:
> 
>> On Tue, 30 Jun 2009 11:10:39 -0600
>>   Charles R Harris <charlesr.harris <at> gmail.com> wrote:
>> > On Tue, Jun 30, 2009 at 10:56 AM, Charles R Harris <
>> > charlesr.harris <at> gmail.com> wrote:
>> >
>> >>
>> >>
>> >> On Tue, Jun 30, 2009 at 10:40 AM, Nils Wagner <
>> >> nwagner <at> iam.uni-stuttgart.de> wrote:
>> >>
>> >>> On Tue, 30 Jun 2009 10:27:05 -0600
>> >>>  Charles R Harris <charlesr.harris <at> gmail.com> wrote:
>> >>> > On Tue, Jun 30, 2009 at 5:11 AM, Nils Wagner
>> >>> > <nwagner <at> iam.uni-stuttgart.de>wrote:
>> >>> >
>> >>> >> On Tue, 30 Jun 2009 11:22:34 +0200
>> >>> >>  "Nils Wagner" <nwagner <at> iam.uni-stuttgart.de> 
>>wrote:
>> >>> >>
>> >>> >>>  Hi all,
>> >>> >>>
>> >>> >>> How can I build the following product with numpy
>> >>> >>>
>> >>> >>> q_i = \varepsilon_{ijk} q_{kj}
>> >>> >>>
>> >>> >>> where  \varepsilon_{ijk} denotes the permutation
>> >>>symbol.
>> >>> >>>
>> >>> >>> Nils
>> >>> >>>
>> >>> >>  Sorry for replying to myself.
>> >>> >> The permutation symbol is also known as the
>> >>>Levi-Civita
>> >>> >>symbol.
>> >>> >> I found an explicit expression at
>> >>> >> http://en.wikipedia.org/wiki/Levi-Civita_symbol
>> >>> >>
>> >>> >> How do I build the product of the Levi-Civita 
>>symbol
>> >>> >>\varepsilon_{ijk} and
>> >>> >> the two dimensional array
>> >>> >> q_{kj}, i,j,k = 1,2,3 ?
>> >>> >>
>> >>> >
>> >>> > Write it out explicitly. It essentially
>> >>>antisymmetrizes
>> >>> >q and the three off
>> >>> > diagonal elements can then be treated as a vector.
>> >>> >Depending on how q is
>> >>> > formed and the resulting vector is used there may 
>>be
>> >>> >other things you can do
>> >>> > when you use it in a more general expression. If 
>>this
>> >>>is
>> >>> >part of a general
>> >>> > calculation there might be other ways of 
>>expressing
>> >>>it.
>> >>> >
>> >>> > Chuck
>> >>>
>> >>> Hi Chuck,
>> >>>
>> >>> Thank you for your response.
>> >>> The problem at hand is described in a paper by 
>>Angeles
>> >>> namely equation (17c) in
>> >>> "Automatic computation of the screw parameters of
>> >>> rigid-body motions.
>> >>> Part I: Finitely-separated positions"
>> >>> Journal of Dynamic systems, Measurement and Control,
>> >>>Vol.
>> >>> 108 (1986) pp. 32-38
>> >>>
>> >>
>> >> You can solve this problem using quaternions also, in
>> >>which case it reduces
>> >> to an eigenvalue problem. You will note that such 
>>things
>> >>as PCA are used in
>> >> the papers that reference the cited work so you can't
>> >>really get around that
>> >> bit of inefficiency.
>> >>
>> >
>> > Here's a reference to the quaternion approach:
>> > 
>>http://people.csail.mit.edu/bkph/papers/Absolute_Orientation.pdf.
>> >You can
>> > get the translation part from the motion of the
>> >centroid.
>> >
>> > If you are into abstractions you will note that the
>> >problem reduces to
>> > minimising a quadratic form in the quaternion
>> >components. The rest is just
>> > algebra ;)
>> >
>> > Chuck
>>
>> It turns out that the product is simply an invariant of 
>>a
>> 3 \times 3 matrix.
>>
>> from numpy import array, zeros, identity
>> from numpy.linalg import norm
>>
>>
>> def vect(A):
>>     """ linear invariant of a 3 x 3 matrix """
>>     tmp = zeros(3,float)
>>     tmp[0] = 0.5*(A[2,1]-A[1,2])
>>     tmp[1] = 0.5*(A[0,2]-A[2,0])
>>     tmp[2] = 0.5*(A[1,0]-A[0,1])
>>
>>     return tmp
> 
> 
> Out of curiosity, where did the .5 come from? It is not 
>normally part of the
> Levi-Civita symbol.
> 
> Chuck

Hi Chuck,

It's my fault. The components of the invariant q are given 
by

q_i = 0.5 \varepsilon_{ijk} q_{kj}

Nils

Re: np.memmap and memory usage

Wed, 01 Jul 2009 10:17:51 +0200, Emmanuelle Gouillart kirjoitti:
> 	I'm using numpy.memmap to open big 3-D arrays of Xray tomography
> data. After I have created a new array using memmap, I modify the
> contrast of every Z-slice (along the first dimension) inside a for loop,
> for a better visualization of the data. Although I call memmap.flush
> after each modification of a Z-slice, the memory used by Ipython keeps
> increasing at every new iteration. At the end of the loop, the memory
> used by Ipython is of the order of magnitude of the size of the data
> file (1.8Go !). I would have expected that the maximum amount of memory
> used would corresponde to only one Z-slice of the 3D array. See the code
> snapshots below for more details.
> 
> 	Is this an expected behaviour? How can I reduce the amount of
> memory used by Ipython and still process my data?

How do you measure the memory used? Note that on Linux, "top" includes 
the size of OS caches for the memmap in the RSS size of a process.
You can try to monitor "free" instead:

$ free
             total       used       free     shared    buffers     cached
Mem:      12300488   11485664     814824          0     642928    7960736
-/+ buffers/cache:    2882000    9418488
Swap:      7847712       2428    7845284

If the memory is used by OS caches, the "used" number on the "-/+ buffers/
cache" line should stay constant while the program runs.

In this case, what is most likely actually taking up memory is the OS 
buffering the data in memory, before writing it to disk. Linux has at 
least some system-wide parameters available that tune the aggressiveness 
of data cachine. I suppose there may also be some file-specific settings, 
but I have no idea what they are.

--

-- 
Pauli Virtanen

Multi thread loading data

Is it possible to use loadtxt in a mult thread way? Basically, I want
to process a very large CSV file (100+ million records) and instead of
loading thousand elements into a buffer process and then load another
1 thousand elements and process and so on...

I was wondering if there is a technique where I can use multiple
processors to do this faster.

TIA

Re: np.memmap and memory usage

A Wednesday 01 July 2009 10:17:51 Emmanuelle Gouillart escrigué:
> 	Hello,
>
> 	I'm using numpy.memmap to open big 3-D arrays of Xray tomography
> data. After I have created a new array using memmap, I modify the
> contrast of every Z-slice (along the first dimension) inside a for loop,
> for a better visualization of the data. Although I call memmap.flush
> after each modification of a Z-slice, the memory used by Ipython keeps
> increasing at every new iteration. At the end of the loop, the memory
> used by Ipython is of the order of magnitude of the size of the data file
> (1.8Go !). I would have expected that the maximum amount of memory used
> would corresponde to only one Z-slice of the 3D array. See the code
> snapshots below for more details.
>
> Is this an expected behaviour? 

I think so, yes.  As all the file (1.8 GB) has to be processed, the OS will 
need to load everything in RAM at least once.  The memory containing already 
processed data is not 'freed' by the OS even if you try to flush() the slices 
(this just force memory data to be saved on-disk), but rather when the OS 
needs it again (but still, they will be swapped out to disk, effectively 
consuming resources).  I'm afraid that the only way to actually 'free' the 
memory is to close() the *complete* dataset.

> How can I reduce the amount of
> memory used by Ipython and still process my data?

There are a number of interfaces that can deal with binary data on-disk so 
that they can read&write slices of data easily.  By using this approach, you 
only have to load the appropriate slice, operate and write it again.  That's 
similar to how numpy.memmap works, but it uses far less memory (just the 
loaded slice).

For example, by using the PyTables interface [1], you could do:

        f = tb.openFile(filename+".h5", "r+")
        data = f.root.data
        for nrow, sl in enumerate(data):
            sl[sl<imin] = imin
            sl[sl>imax] = imax
            data[nrow] = sl
        f.close()

which is similar to your memmap-based approach:

        data = np.memmap(filename+".bin", mode='r+', dtype='f4', shape=shape)
        for sl in data:
            sl[sl<imin] = imin
            sl[sl>imax] = imax

just that it takes far less memory (47 MB vs 1.9 GB) and besides, it has no 
other limitation than your available space on disk (compare this with the 
limit of your virtual memory when using memmap).  The speeds are similar too:

Using numpy.memmap                                                            
Time creating data file: 8.534                                                
Time processing data file: 32.742                                             

Using tables                                                                  
Time creating data file: 2.88                                                 
Time processing data file: 32.615                                             

However, you can still speed-up out-of-core computations by using the recently 
introduced tables.Expr class (PyTables 2.2b1, see [2]), which uses a 
combination of the Numexpr [3] and PyTables advanced computing capabilities:

        f = tb.openFile(filename+".h5", "r+")
        data = f.root.data
        expr = tb.Expr("where(data<imin, imin, data)")
        expr.setOutput(data)
        expr.eval()
        expr = tb.Expr("where(data>imax, imax, data)")
        expr.setOutput(data)
        expr.eval()
        f.close()

and the timings for this venue are:

Using tables.Expr                                                             
Time creating data file: 2.393                                                
Time processing data file: 18.25                                              

which is around a 75% faster than a pure memmap/PyTables approach.

Further, if your data is compressible, you can probably achieve additional 
speed-ups by using a fast compressor (like LZO, which is supported by PyTables 
right out-of-the-box).

I'm attaching the script I've used for producing the above timings.  You may 
find it useful for trying this out against your own data.

[1] http://www.pytables.org
[2] http://www.pytables.org/download/preliminary/
[3] http://code.google.com/p/numexpr/

HTH,

--

-- 
Francesc Alted

_______________________________________________
Numpy-discussion mailing list
Numpy-discussion <at> scipy.org
http://mail.scipy.org/mailman/listinfo/numpy-discussion

Re: np.memmap and memory usage

	Hi Pauli,

	thank you for your answer! I was indeed measuring the memory used
with top, which is not the best tool for understanding what really
happens. I monitored "free" during the execution of my program and
indeed, the used numbers on the "+/-buffers/cache" line stays roughly
constant (fluctuations are smaller than say 5Mo). When I close my program
this number gets smaller by a number of the order of magnitude of one
Z-slice of my array (35 Mo), so it is possible that the buffer is always
full. However, as I call memmap.flush, the buffer should be emptied at
each iteration, right? 

	I observed other strange things, for example if I create a first
array with np.memmap, I get a MemoryError if I try to create a second
array (without doing anything with the first array but creating it). The
sum of the sizes of the two files is 3.8 Go, much more than my RAM size,
but as I didn't load anything explicitely into the memory, I don't
understand why it is not possible to have both arrays defined...

	Anyway, I need to read a bit of documentation about data caching
in Linux, as you suggested.

	Thanks again,

	Emmanuelle

> > 	Is this an expected behaviour? How can I reduce the amount of
> > memory used by Ipython and still process my data?

> How do you measure the memory used? Note that on Linux, "top" includes 
> the size of OS caches for the memmap in the RSS size of a process.
> You can try to monitor "free" instead:

> $ free
>              total       used       free     shared    buffers     cached
> Mem:      12300488   11485664     814824          0     642928    7960736
> -/+ buffers/cache:    2882000    9418488
> Swap:      7847712       2428    7845284

> If the memory is used by OS caches, the "used" number on the "-/+ buffers/
> cache" line should stay constant while the program runs.

> In this case, what is most likely actually taking up memory is the OS 
> buffering the data in memory, before writing it to disk. Linux has at 
> least some system-wide parameters available that tune the aggressiveness 
> of data cachine. I suppose there may also be some file-specific settings, 
> but I have no idea what they are.

Re: np.memmap and memory usage

	Hi Francesc,

	many thanks for this very detailed and informative answer! This
list is really great :D.

	I'm going to install pytables at once and I will try your scripts
with my data. As your computations were made with files of the same size
as mine, hopefully it should run also quite well on my computer. I'll let
you know! 

	Thanks again,

	Emmanuelle

On Wed, Jul 01, 2009 at 03:04:08PM +0200, Francesc Alted wrote:
> A Wednesday 01 July 2009 10:17:51 Emmanuelle Gouillart escrigué:
> > 	Hello,

> > 	I'm using numpy.memmap to open big 3-D arrays of Xray tomography
> > data. After I have created a new array using memmap, I modify the
> > contrast of every Z-slice (along the first dimension) inside a for loop,
> > for a better visualization of the data. Although I call memmap.flush
> > after each modification of a Z-slice, the memory used by Ipython keeps
> > increasing at every new iteration. At the end of the loop, the memory
> > used by Ipython is of the order of magnitude of the size of the data file
> > (1.8Go !). I would have expected that the maximum amount of memory used
> > would corresponde to only one Z-slice of the 3D array. See the code
> > snapshots below for more details.

> > Is this an expected behaviour? 

> I think so, yes.  As all the file (1.8 GB) has to be processed, the OS will 
> need to load everything in RAM at least once.  The memory containing already 
> processed data is not 'freed' by the OS even if you try to flush() the slices 
> (this just force memory data to be saved on-disk), but rather when the OS 
> needs it again (but still, they will be swapped out to disk, effectively 
> consuming resources).  I'm afraid that the only way to actually 'free' the 
> memory is to close() the *complete* dataset.

> > How can I reduce the amount of
> > memory used by Ipython and still process my data?

> There are a number of interfaces that can deal with binary data on-disk so 
> that they can read&write slices of data easily.  By using this approach, you 
> only have to load the appropriate slice, operate and write it again.  That's 
> similar to how numpy.memmap works, but it uses far less memory (just the 
> loaded slice).

> For example, by using the PyTables interface [1], you could do:

>         f = tb.openFile(filename+".h5", "r+")
>         data = f.root.data
>         for nrow, sl in enumerate(data):
>             sl[sl<imin] = imin
>             sl[sl>imax] = imax
>             data[nrow] = sl
>         f.close()

> which is similar to your memmap-based approach:

>         data = np.memmap(filename+".bin", mode='r+', dtype='f4', shape=shape)
>         for sl in data:
>             sl[sl<imin] = imin
>             sl[sl>imax] = imax

> just that it takes far less memory (47 MB vs 1.9 GB) and besides, it has no 
> other limitation than your available space on disk (compare this with the 
> limit of your virtual memory when using memmap).  The speeds are similar too:

> Using numpy.memmap                                                            
> Time creating data file: 8.534                                                
> Time processing data file: 32.742                                             

> Using tables                                                                  
> Time creating data file: 2.88                                                 
> Time processing data file: 32.615                                             

> However, you can still speed-up out-of-core computations by using the recently 
> introduced tables.Expr class (PyTables 2.2b1, see [2]), which uses a 
> combination of the Numexpr [3] and PyTables advanced computing capabilities:

>         f = tb.openFile(filename+".h5", "r+")
>         data = f.root.data
>         expr = tb.Expr("where(data<imin, imin, data)")
>         expr.setOutput(data)
>         expr.eval()
>         expr = tb.Expr("where(data>imax, imax, data)")
>         expr.setOutput(data)
>         expr.eval()
>         f.close()

> and the timings for this venue are:

> Using tables.Expr                                                             
> Time creating data file: 2.393                                                
> Time processing data file: 18.25                                              

> which is around a 75% faster than a pure memmap/PyTables approach.

> Further, if your data is compressible, you can probably achieve additional 
> speed-ups by using a fast compressor (like LZO, which is supported by PyTables 
> right out-of-the-box).

> I'm attaching the script I've used for producing the above timings.  You may 
> find it useful for trying this out against your own data.

> [1] http://www.pytables.org
> [2] http://www.pytables.org/download/preliminary/
> [3] http://code.google.com/p/numexpr/

> HTH,