code-generation

Question

Hi Folks,

Here is a simple example of using Cadabra to generate skeleton C-code for a given typical tensor expression. Note that the C99CodePrinter requires sympy 1.1.1.

There is a small bug somewhere -- using integers causes a syntax error in sympy. But writing the same integer as a float avoids the problem.

Cheers, Leo

{a,b,c}::Indices(position=fixed,values={x,y,z}).
{x,y,z}::Coordinate.

ex := 2.0 M A_{a} A^{b};   # okay
# ex := 2 A_{a} A^{b};     # sympy parse error: 2Ax**2

display(ex)

rl1 := {A_{x}=Ax,A_{y}=Ay,A_{z}=Az};
rl2 := {A^{x}=Ax,A^{y}=Ay,A^{z}=Az};

evaluate(ex,rl1+rl2,rhsonly=True);

def writecode (the_ex,name,filename,num):

   import sympy as sym

   # following requires sympy 1.1.1
   from sympy.printing.ccode import C99CodePrinter as printer
   from sympy.printing.codeprinter import Assignment

   idx=[]  # will contain all indices in the form [{x, x}, {x, y} ...]
   lst=[]  # will contain the corresponding expressions [termxx, termxy, ...]

   for i in range( len(the_ex[num]) ):       # num = the number of free indices
       idx.append( str(the_ex[num][i][0]) )  # indices for this term
       lst.append( str(the_ex[num][i][1]) )  # the term matching these indices

   mat = sym.Matrix([lst])

   sub_exprs, simplified_rhs = sym.cse(mat)

   with open(filename, 'w') as f:

      for var, sub_expr in sub_exprs:
         f.write(printer().doprint(Assignment(var, sub_expr))+'\n')

      for index, term in enumerate (simplified_rhs[0]):
         var = sym.Symbol(name+' ('+idx[index][1:-1]+')')
         f.write(printer().doprint(Assignment(var, term))+'\n')

writecode (ex,'gab','foo.c',2)

Answer 1

The problem with the '2' has to do with the fact that you use str to convert the Ex objects to strings which you then feed to Python. You can call _sympy_() on the Ex instead to get a Sympy object, but the cleanest would be to have a method which gives a sympy parseable string.

Comments

Hi Kasper,

Thanks for the tip. I've made the change and now there's no problem with using 2 or 2.0. But another small problem has cropped up. When I use

ex := 2 M A_{a} A^{b};

everything works fine. But when I use

ex := 2 N A_{a} A^{b};

the c-dode contains "sympyN" rather than "N". I can easily post-process the c-ode to fix it but I wonder if I can avoid that? I've attached the modified code below.

Cheers, Leo

{a,b,c}::Indices(position=fixed,values={x,y,z}).
{x,y,z}::Coordinate.

ex := 2 M A_{a} A^{b};       # correct c-code
# ex := 2 N A_{a} A^{b};     # c-code contains "sympyN" not "N"

display(ex)

rl1 := {A_{x}=Ax,A_{y}=Ay,A_{z}=Az};
rl2 := {A^{x}=Ax,A^{y}=Ay,A^{z}=Az};

evaluate(ex,rl1+rl2,rhsonly=True);

def writecode (the_ex,name,filename,num):

   import sympy as sym

   # following requires sympy 1.1.1
   from sympy.printing.ccode import C99CodePrinter as printer
   from sympy.printing.codeprinter import Assignment

   idx=[]  # will contain all indices in the form [{x, x}, {x, y} ...]
   lst=[]  # will contain the corresponding expressions [termxx, termxy, ...]

   for i in range( len(the_ex[num]) ):       # num = the number of free indices
       idx.append( str(the_ex[num][i][0]) )  # indices for this term
       lst.append( (the_ex[num][i][1])._sympy_() )  # the term matching these indices

   mat = sym.Matrix([lst])

   sub_exprs, simplified_rhs = sym.cse(mat)

   with open(filename, 'w') as f:

      for var, sub_expr in sub_exprs:
         f.write(printer().doprint(Assignment(var, sub_expr))+'\n')

      for index, term in enumerate (simplified_rhs[0]):
         var = sym.Symbol(name+' ('+idx[index][1:-1]+')')
         f.write(printer().doprint(Assignment(var, term))+'\n')

writecode (ex,'gab','foo.c',2)

---

That's an inevitable consequence of ending up with sympy expressions. If N is a reserved symbol in sympy, you cannot use it in your own expressions.

---

That makes sense. I can live with that (it's a trivial sed-script to clean the c-code). Thanks again.