[Numpy-discussion] Numpy and PEP 343

[David M. Cooke]

Tim Hochberg <tim.hochberg at cox.net> writes:

> David M. Cooke wrote:
>
>>Wow. On par with weave.blitz, and no C++ compiler!!! :-)
>>
>>
> That's awesome! I was also tempted by this, but I never got beyond
> prototyping some stuff in Python.
>
>>You use it like this:
>>
>>from numexpr import numexpr
>>
>>func = numexpr("2*a+3*b")
>>
>>a = arange(1e6)
>>b = arange(1e6)
>>c = func(a, b)
>>
>>
> Does this just uses the order that variable are initially used in the
> expression to determine the input order? I'm not sure I like that. It
> is very convenient for simple expressions though.

By default, it does lexical order (whatever .sort() on the list of
variable names gives). You can specify the order with

func = numexpr("2*a+3*b", input_order=('a', 'b'))

>>Alternatively, you can use it like weave.blitz, like this:
>>
>>from numexpr import evaluate
>>
>>a = arange(1e6)
>>b = arange(1e6)
>>c = evaluate("2*a+3*b")
>>
>>
> That's pretty sweet. And I see that you cache the expressions, so it
> should be pretty fast if you need to loop.
>
> [snip details]
>
>>If people think it's useful enough, I can check it into scipy's sandbox.
>>
>>
> Definitely check it in. It won't compile here with VC7, but I'll see
> if I can figure out why.

It's in. Check the setup.py -- I left some gcc'isms in there
(-funroll-loops). That should be conditional on the compiler.

> This is probably thinking two far ahead, but an interesting thing to
> try would be adding conditional expressions:
>
> c = evaluate("(2*a + b) if (a > b) else (2*b + a)")
>
> If that could be made to work, and work fast, it would save both
> memory and time in those cases where you have to vary the computation
> based on the value. At present I end up computing the full arrays for
> each case and then choosing which result to use based on a mask, so it
> takes three times as much space as doing it element by element.

It's on my list of things to think about, along with implementing
reductions (sum, in particular). It'd probably look more like

c = evaluate("where(a > b, 2*a+b, 2*b+a)")

because of the vector nature of the code. Doing the "if" elementwise
would mean the bytecode program would have to be run for each element,
instead of on a vector of some fixed length.

I wrote up a bit on my blog at http://isobaric.blogspot.com/ (finally,
something to blog about :-)

-- 
|>|\/|<
/--------------------------------------------------------------------------\
|David M. Cooke                      http://arbutus.physics.mcmaster.ca/dmc/
|cookedm at physics.mcmaster.ca