Possible to use numexpr with user made ufuncs/scipy ufuncs?
![](https://secure.gravatar.com/avatar/725a50eff7b51f36402b946bf786393e.jpg?s=120&d=mm&r=g)
Hello, Does anyone know whether it is possible to use numexpr with scipy ufuncs (such as those in scipy.special) or user made ufuncs? This functionality would be extremely useful. I don't see them in the list of supported functions (http://code.google.com/p/numexpr/wiki/Overview) and I get a 'TypeError: 'VariableNode' object is not callable' error when I try this. Best Regards, John Salvatier
![](https://secure.gravatar.com/avatar/bd813aa01a8be59c71e6c797194beffe.jpg?s=120&d=mm&r=g)
A Friday 25 June 2010 21:33:43 John Salvatier escrigué:
Hello,
Does anyone know whether it is possible to use numexpr with scipy ufuncs (such as those in scipy.special) or user made ufuncs? This functionality would be extremely useful. I don't see them in the list of supported functions (http://code.google.com/p/numexpr/wiki/Overview) and I get a 'TypeError: 'VariableNode' object is not callable' error when I try this.
Yeah, you need to explicitly code the support for new functions in numexpr. But another possibility, more doable, would be to code the scipy.special functions by using numexpr as a computing back-end. However, I understand that many scipy.special functions need to evaluate basic transcendental functions (trignonometrical, exponential, logarithmic...), so you should no expect big speed-ups just by using a plain Numexpr. However, you can always link Numexpr with Intel's VML (Vector Math Library), which can use SIMD and multi-core capabilities in modern Intel/AMD processors so as to accelerate the evaluation of such transcendental functions (typically between 2x and 10x, depending on the function and number of cores in your processor). -- Francesc Alted
![](https://secure.gravatar.com/avatar/da3a0a1942fbdc5ee9a9b8115ac5dae7.jpg?s=120&d=mm&r=g)
Hi, la, 2010-06-26 kello 14:24 +0200, Francesc Alted kirjoitti: [clip]
Yeah, you need to explicitly code the support for new functions in numexpr. But another possibility, more doable, would be to code the scipy.special functions by using numexpr as a computing back-end.
Would it be possible to add generic support for user-supplied ufuncs into numexpr? This would maybe be more of a convenience feature than a speed improvement, although perhaps some speed could be gained by evaluating ufuncs per-block. -- Pauli Virtanen
![](https://secure.gravatar.com/avatar/725a50eff7b51f36402b946bf786393e.jpg?s=120&d=mm&r=g)
OK, perhaps I will just continue preevaluating the expressions involving special functions. Thank you for the help Fancesc +1 to Pauli's suggestion. On Sat, Jun 26, 2010 at 6:19 AM, Pauli Virtanen <pav@iki.fi> wrote:
Hi,
la, 2010-06-26 kello 14:24 +0200, Francesc Alted kirjoitti: [clip]
Yeah, you need to explicitly code the support for new functions in numexpr. But another possibility, more doable, would be to code the scipy.special functions by using numexpr as a computing back-end.
Would it be possible to add generic support for user-supplied ufuncs into numexpr? This would maybe be more of a convenience feature than a speed improvement, although perhaps some speed could be gained by evaluating ufuncs per-block.
-- Pauli Virtanen
_______________________________________________ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion
![](https://secure.gravatar.com/avatar/bd813aa01a8be59c71e6c797194beffe.jpg?s=120&d=mm&r=g)
2010/6/26 Pauli Virtanen <pav@iki.fi>
Hi,
la, 2010-06-26 kello 14:24 +0200, Francesc Alted kirjoitti: [clip]
Yeah, you need to explicitly code the support for new functions in numexpr. But another possibility, more doable, would be to code the scipy.special functions by using numexpr as a computing back-end.
Would it be possible to add generic support for user-supplied ufuncs into numexpr? This would maybe be more of a convenience feature than a speed improvement, although perhaps some speed could be gained by evaluating ufuncs per-block.
Well, I'd say that this support can be faked in numexpr easily. For example, if one want to compute a certain ufunc called, say, sincos(x) defined as sin(cos(x)) (okay, that's very simple, but it will suffice for demonstration purposes), he can express that as sincos = "sin(cos(%s))", and then use it in a more complex expression like: "%s+1*cos(%s)" % (sincos % 'x', 'y') that will be expanded as: "sin(cos(x))+1*cos(y)" Of course, this is a bit crude, but I'd say that's more than enough for allowing the evaluation of moderately complex expressions. Having said this, one can always think in setting up a wrapper for doing a similar thing more elegantly (although I'm not sure if this is worth the effort). -- Francesc Alted
![](https://secure.gravatar.com/avatar/da3a0a1942fbdc5ee9a9b8115ac5dae7.jpg?s=120&d=mm&r=g)
Sat, 26 Jun 2010 17:51:52 +0100, Francesc Alted wrote: [clip]
Well, I'd say that this support can be faked in numexpr easily. For example, if one want to compute a certain ufunc called, say, sincos(x) defined as sin(cos(x)) (okay, that's very simple, but it will suffice for demonstration purposes), he can express that as sincos = "sin(cos(%s))", and then use it in a more complex expression like:
"%s+1*cos(%s)" % (sincos % 'x', 'y')
that will be expanded as:
"sin(cos(x))+1*cos(y)"
Of course, this is a bit crude, but I'd say that's more than enough for allowing the evaluation of moderately complex expressions.
But what if such an expression does not exist? For example, there is no finite closed form expression for the Bessel functions in terms of elementary functions. An accurate implementation is rather complicated, and usually piecewise defined. For convenience reasons, it could be useful if one could do something like numexpr.evaluate("cos(iv(0, x))", functions=dict(iv=scipy.special.iv)) and this would be translated to numexpr bytecode that would make a Python function call to obtain "iv(0, x)" for each block of data required, assuming "iv" is a vectorized function. It's of course possible to precompute the value of "iv(0, x)", but this is extra hassle and requires additional memory. -- Pauli Virtanen
![](https://secure.gravatar.com/avatar/bd813aa01a8be59c71e6c797194beffe.jpg?s=120&d=mm&r=g)
A Saturday 26 June 2010 19:17:43 Pauli Virtanen escrigué:
But what if such an expression does not exist? For example, there is no finite closed form expression for the Bessel functions in terms of elementary functions. An accurate implementation is rather complicated, and usually piecewise defined.
For convenience reasons, it could be useful if one could do something like
numexpr.evaluate("cos(iv(0, x))", functions=dict(iv=scipy.special.iv))
and this would be translated to numexpr bytecode that would make a Python function call to obtain "iv(0, x)" for each block of data required, assuming "iv" is a vectorized function. It's of course possible to precompute the value of "iv(0, x)", but this is extra hassle and requires additional memory.
Yeah, I can see the merit of implementing such a thing, mainly for avoiding additional memory consumption. But again, the nice thing would be to implement such a special functions in terms of numexpr expressions so that the evaluation itself can be faster. Admittedly, that would take a bit more time. Anyway, if someone comes with patches for implementing this, I'd glad to commit them. -- Francesc Alted
![](https://secure.gravatar.com/avatar/da3a0a1942fbdc5ee9a9b8115ac5dae7.jpg?s=120&d=mm&r=g)
ma, 2010-06-28 kello 09:48 +0200, Francesc Alted kirjoitti: [clip]
But again, the nice thing would be to implement such a special functions in terms of numexpr expressions so that the evaluation itself can be faster. Admittedly, that would take a bit more time.
Quite often, you need to evaluate a series or continued fraction expansion to get the value of a special function at some point, limiting the number of terms by stopping when a certain convergence criterion is satisfied. Also, which series to sum typically depends on the point where things are evaluated. These things don't vectorize very nicely. If I understand correctly, numexpr bytecode interpreter does not support conditionals and loops like this at the moment? The special function implementations are typically written in bare C/Fortran. Do you think numexpr could give speedups there? As I see it, speedups (at least with MKL) could come from using faster implementations of sin/cos/exp etc. basic functions. Using SIMD to maximum effect would then require an amount of cleverness in re-writing the evaluation algorithms, which sounds like a major amount of work. -- Pauli Virtanen
![](https://secure.gravatar.com/avatar/bd813aa01a8be59c71e6c797194beffe.jpg?s=120&d=mm&r=g)
A Monday 28 June 2010 10:22:31 Pauli Virtanen escrigué:
ma, 2010-06-28 kello 09:48 +0200, Francesc Alted kirjoitti: [clip]
But again, the nice thing would be to implement such a special functions in terms of numexpr expressions so that the evaluation itself can be faster. Admittedly, that would take a bit more time.
Quite often, you need to evaluate a series or continued fraction expansion to get the value of a special function at some point, limiting the number of terms by stopping when a certain convergence criterion is satisfied. Also, which series to sum typically depends on the point where things are evaluated. These things don't vectorize very nicely. If I understand correctly, numexpr bytecode interpreter does not support conditionals and loops like this at the moment?
No, it does not support loops. And conditionals are supported only via vectorized conditionals (i.e. the `where()` opcode).
The special function implementations are typically written in bare C/Fortran. Do you think numexpr could give speedups there? As I see it, speedups (at least with MKL) could come from using faster implementations of sin/cos/exp etc. basic functions. Using SIMD to maximum effect would then require an amount of cleverness in re-writing the evaluation algorithms, which sounds like a major amount of work.
Okay. I thought that special functions were more 'vectorizable', or that it could be expressed in terms of more basic functions (sin/cos/exp). But if this is not the case, then it is quite clear that the best solution would be your first suggestion (i.e. implement user-provide ufunc evaluation). And definitely, implementing special functions in terms of SIMD would really be a *major* effort, and only doable by very specialized people ;-) -- Francesc Alted
participants (3)
-
Francesc Alted
-
John Salvatier
-
Pauli Virtanen