[Cython] aritmetic with arrays in Cython

[Julian Taylor]

On 08.08.2014 19:36, Stefan Behnel wrote:
> Hi,
> 
> please don't top-post.
> 
> Ian Henriksen schrieb am 08.08.2014 um 18:47:
>> I'd like to work on it. That was why I asked. It's been bothering me for
>> some time that we don't have this feature yet. This will take me some time
>> though since I'm not at all familiar with Cython's internals. I've used it
>> in relatively small settings for speeding up array operations and wrapping
>> external functions, but I've never had occasion to dig through much of what
>> it does internally. Some help navigating Cython's internals would be
>> greatly appreciated.
> 
> Here's a little intro.
> 
> https://github.com/cython/cython/wiki/HackerGuide#getting-started
> 
> I recommend to start by writing a simple test that does some integer
> arithmetic and enable code generation traces for AST nodes in the generated
> C code. "cython -a" should make this quite navigable.
> 
> 
>> I'll also have to take some time to become more familiar with eigen itself.
>> Most of my work has used numpy arrays.
>>
>> Eigen seems to be ideal for this since it allows us to offload the
>> expression analysis to another package where it is already
>> well-implemented. It looks like they have an array class for n-dimensional
>> array operations as well. My guess is that that would allow a fairly
>> natural translation from memory views to eigen arrays.
> 
> You should dig into the code that Ceygen uses for the mapping. That's most
> likely the best start you can get.
> 
> 
>> Would it be more helpful for me to start working on Cython bindings for
>> eigen, or to work on adding this directly as a part of Cython? Where should
>> I look in Cython's source code to do that? What kinds of modifications
>> would be necessary?
> 
> My guess is that this is best implemented with a set of matrix specific AST
> nodes for the arithmetic operators. There is a recursive tree processing
> phase in Cython's pipeline called "AnalyseExpressionsTransform" or type
> analysis, where it figures out what variables reference memory views, what
> the result type of an arithmetic expression is, etc. This phase can easily
> replace nodes by returning new nodes from the analyse_types() methods of
> AST nodes. Currently, arithmetic with memory views is not allowed, so this
> needs to be enabled in the corresponding NumBinopNode AST subclasses
> (AddNode, MatMulNode, etc. in ExprNodes.py). Write a test and debug it to
> see where it fails.

If someone wants to work on an ast transformer for array expressions I
would not restrict its output to something ceygen can parse but also numpy.
E.g. transforming:

cpdef whateverfunction(double[:] a, float[:] b):
    return a + a * b + .5

to

cpdef whateverfunction(double[:] a, float[:] b):
    bs = 10000
    r = empty_like(a)
    for i in range(0, a.size, bs):
        u = min(i + bs, a.size)
        r[i:u] = a[i:u] * b[i:u]
        r[i:u] += a[i:u]
        r[i:u] += 0.5
    return r

this is already very efficient in numpy.