code generation

[Loïc Gouarin]

Hi,

I took a look at the GSoC 2017 ideas and I saw that one of the topics was the code generation. You said that there is work on updating the system ongoing. I'm really interested to know what this work is about.

I will try to explain why. I develop a Python package call pyLBM (https://www.math.u-psud.fr/pyLBM/) to use lattice Boltzmann methods. The idea is that a user specifies his symbolic scheme for LBM using sympy and then we generate the numerical code associated. In the master of pyLBM, we implement our own template to write cython and numpy versions. But, few weeks ago, I have decided to look at the code generation in sympy because we want to add other implementation like pythran or loopy easily and it is a lot of pain with our template. I think also that with sympy it could be easy to add some strategies during the code generation like cache blocking, SIMD, loop unrolling, ...

The codeprinters in sympy are very useful and save a lot of time. On the other side, I had some difficulties with the original codegen to generate a numerical code so I have to rewrite it. For the moment, I have a cython version but it is very extensible. I will implement the numpy, pythran and loopy versions in few days.

The current implementation of codegen support only MatrixSymbol, so I had to add Matrix expression. I also supposed that all expression must be an equality (for the moment). It was also difficult to write loop with only the Idx so I add a sympy expression called Loop to set a loop block in my codegen.

I show you just two examples of the input and output.

First example
==========

Input:

from pyLBM.generator import codegen
from sympy import *

M = Matrix([[1/4,  1/2,     0,  1/4],
            [1/4,     0,  1/2, -1/4],
            [1/4, -1/2,     0,  1/4],
            [1/4,     0, -1/2, -1/4]])

f = MatrixSymbol('f', 4, 1)
m = MatrixSymbol('m', 4, 1)

res = codegen(('m2f', Eq(f, M*m)), "cython", header=False)
print(res[0][1]

Output:

def m2f(double[::1] f, double[::1] m):

    f[0] = 0.25*m[0] + 0.5*m[1] + 0.25*m[3]
    f[1] = 0.25*m[0] + 0.5*m[2] - 0.25*m[3]
    f[2] = 0.25*m[0] - 0.5*m[1] + 0.25*m[3]
    f[3] = 0.25*m[0] - 0.5*m[2] - 0.25*m[3]

Second example
============

Input:

from pyLBM.generator import codegen, Loop
from sympy import *
from sympy.abc import z

n, m = symbols('n m', integer=True)
A = IndexedBase('A', shape=(m, n, 4))
ff = IndexedBase('f', shape=(m, n, 4))
i = Idx('i', m)
j = Idx('j', n)
k = Idx('k', 4)

m = MatrixSymbol('m', 4, 1)

mnew = Matrix([A[i, j, kk] for kk in range(4)])
fnew = Matrix([ff[i, j, kk] for kk in range(4)])

eq = Matrix([[m[0]],
             [m[1]],
             [m[2]],
             [ 0.0]])

for ii in range(mnew.shape[0]):
    eq = eq.subs(m[ii], mnew[ii])


l2 = Loop('l2', k, Eq(z,k, evaluate=False))
l1 = Loop('l1', (i, j), [Eq(fnew, eq, evaluate=False), l2])

res = codegen(('test', l1), "CYTHON", header=False)
print(res[0][1])

Output:

def test(double[:, :, ::1] f, double[:, :, ::1] A, double z):
    cdef int j
    cdef int k
    cdef int i


    for i in range(0, m):
        for j in range(0, n):
            f[i, j, 0] = A[i, j, 0]
            f[i, j, 1] = A[i, j, 1]
            f[i, j, 2] = A[i, j, 2]
            f[i, j, 3] = 0.0
            for k in range(0, 4):
                z = k

[Loïc Gouarin]

I published accidentally so I finish my message here. Sorry...

In the last example, I have to add n and m as local variables which depends on the shape of A or f.

The most difficult for me in this implementation is the Idx: are they local, in a loop, in the function parameters, ... And this defines the block structure of the numerical code.

Anyway, I am really interested to know what you have done in the new version of code generation in sympy.

Could you tell me where I can find this version ?

Thanks,
Loic

[Aaron Meurer]

Most of the work so far has been in improving the code printers. You
can see there are some pull requests to add support for some new
languages https://github.com/sympy/sympy/pulls?q=is%3Aopen+is%3Apr+label%3Acodegen.
There has also been work in sympy.codegen for supporting higher level
abstractions like loops.

Aaron Meurer

> --
> You received this message because you are subscribed to the Google Groups
> "sympy" group.
> To unsubscribe from this group and stop receiving emails from it, send an
> email to sympy+un...@googlegroups.com.
> To post to this group, send email to sy...@googlegroups.com.
> Visit this group at https://groups.google.com/group/sympy.
> To view this discussion on the web visit
> https://groups.google.com/d/msgid/sympy/89da10b4-a954-40ab-8f11-a17594eb7d8f%40googlegroups.com.
>
> For more options, visit https://groups.google.com/d/optout.

[Loïc Gouarin]

Thanks Aaron for the links. But it is just new languages support. So, the core of code printers in SymPy is the same.

I can add support for Cython, loopy and Pythran but I need more informations about the higher level
abstractions like loops. Is it your CodeBlocks and ast file ? What is the best way to start ?

I use SymPy from pip and then I didn't see theses changes. So it seems that it is better to work with the master branch of SymPy.

[Anthony Scopatz]

Hi Loïc, 

Yes, I think it is better to work from master for this stuff, which is still in development.

Be Well

Anthony

--
You received this message because you are subscribed to the Google Groups "sympy" group.
To unsubscribe from this group and stop receiving emails from it, send an email to sympy+un...@googlegroups.com.
To post to this group, send email to sy...@googlegroups.com.
Visit this group at https://groups.google.com/group/sympy.

To view this discussion on the web visit https://groups.google.com/d/msgid/sympy/be387a74-d6bb-48e9-921e-8b0c37365639%40googlegroups.com.

For more options, visit https://groups.google.com/d/optout.

--

Asst. Prof. Anthony Scopatz
Nuclear Engineering Program
Mechanical Engineering Dept.
University of South Carolina
sco...@cec.sc.edu
Office: (803) 777-9334
Cell: (512) 827-8239
Office: Horizon 1, #011 (ground floor, not first)
Book a meeting with me at https://scopatz.youcanbook.me/
Open up an issue: https://github.com/scopatz/me/issues
Check my calendar