[Cython] aritmetic with arrays in Cython

[Stefan Behnel]

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.

I'm not sure yet how the code generation would look like. In simple cases,
you may be able to implement this into the existing AST arithmetic operator
nodes. But my guess is that you want to have a new set of nodes (probably
in a new module in Cython/Compiler/, e.g. EigenNodes.py) that the
analyse_types() methods of the generic arithmetic AST nodes return instead
of themselves when they determine that they are operating on memory views.
Those would then know how to generate matrix specific arithmetic C++/Eigen
code, potentially even inherit from the general arithmetic nodes to
override only the code generation method. There is a
NewNodeClass.from_node() class method to create a new node from an existing
one.

These specialised AST nodes will be optional in the end, Eigen should only
be needed if in C++ mode and this feature is enabled and used, otherwise,
array arithmetic will continue to be a compile time error.

Try to get something simple like binary matrix multiplication working for
now, or even just matrix x scalar. Map them to Ceygen's API to get a quick
start. Once that's done, try building up the code generation for more
complex nested expressions.

Stefan