[Cython] aritmetic with arrays in Cython

[Ian Henriksen]

On Fri, Aug 8, 2014 at 3:09 PM, Ian Henriksen <
insertinterestingnamehere at gmail.com> wrote:

> On Fri, Aug 8, 2014 at 12:02 PM, Julian Taylor <
> jtaylor.debian at googlemail.com> wrote:
>
>> 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.
>>
>>
> Sorry about the top post.
>
> Stefan, thanks for the direction. I think I know how to proceed now. It'll
> probably take me some time to learn how to do all this, but I'll start
> reading/working on it now.
>
> Julian, from a dependency standpoint I'd really prefer to use numpy. As
> near as I can tell, the people most likely to use this feature already use
> numpy. That would also remove the C++ dependency for this. I think the main
> reason for using eigen is that it will be easy to interface with and will
> eliminate temporaries automatically. In the long term, it could be good to
> offer multiple backends for these sorts of operations. Numpy would be a
> primary candidate for something like that. I'll spend some time working
> with numpy's ufunc api to see if I can figure out a good way to use that
> instead.
> Thanks!
> -Ian
>

Maybe I should clarify a little about why eigen is a good place to start.
According to http://eigen.tuxfamily.org/dox/TopicLazyEvaluation.html it
already takes care of things like the elimination of temporary variables
and common subexpression reduction at compile time. This should make it
possible to compile array expressions in Cython without having to
re-implement those sorts of optimizations. Ideally we would just have to
map memory view operations to corresponding equivalents from eigen. It's
not yet clear to me how to do things with arbitrary-dimensional arrays or
broadcasting, but, given some more time, a solution may present itself.
-Ian
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://mail.python.org/pipermail/cython-devel/attachments/20140808/3066a46d/attachment-0001.html>