Geek Answers Handbook

Optimize sympy code

Using SymPy to find the derivative (see this question: https://math.stackexchange.com/questions/726104/apply-chain-rule-to-vector-function-with-chained-dot-and-cross-product ) I came up with this code:

from sympy import * from sympy.physics.mechanics import * from sympy.printing import print_ccode from sympy.utilities.codegen import codegen x1, x2, x3 = symbols('x1 x2 x3') y1, y2, y3 = symbols('y1 y2 y3') z1, z2, z3 = symbols('z1 z2 z3') u = ReferenceFrame('u') u1=(ux*x1 + uy*y1 + uz*z1) u2=(ux*x2 + uy*y2 + uz*z2) u3=(ux*x3 + uy*y3 + uz*z3) s1=(u1-u2).normalize() s2=(u2-u3).normalize() v=cross(s1, s2) f=dot(v,v) df_dy2=diff(f, y2) print_ccode(df_dy2, assign_to='df_dy2') [(c_name, c_code), (h_name, c_header)] = codegen( ("df_dy2", df_dy2), "C", "test", header=False, empty=False) print c_code

What gives this beauty:

 #include "test.h" #include <math.h> double df_dy2(double x1, double x2, double x3, double y1, double y2, double y3, double z1, double z2, double z3) { return ((x1 - x2)*(y2 - y3)/(sqrt(pow(x1 - x2, 2) + pow(y1 - y2, 2) + pow(z1 - z2, 2))*sqrt(pow(x2 - x3, 2) + pow(y2 - y3, 2) + pow(z2 - z3, 2))) - (x2 - x3)*(y1 - y2)/(sqrt(pow(x1 - x2, 2) + pow(y1 - y2, 2) + pow(z1 - z2, 2))*sqrt(pow(x2 - x3, 2) + pow(y2 - y3, 2) + pow(z2 - z3, 2))))*(2*(x1 - x2)*(y1 - y2)*(y2 - y3)/(pow(pow(x1 - x2, 2) + pow(y1 - y2, 2) + pow(z1 - z2, 2), 3.0L/2.0L)*sqrt(pow(x2 - x3, 2) + pow(y2 - y3, 2) + pow(z2 - z3, 2))) + 2*(x1 - x2)*(-y2 + y3)*(y2 - y3)/(sqrt(pow(x1 - x2, 2) + pow(y1 - y2, 2) + pow(z1 - z2, 2))*pow(pow(x2 - x3, 2) + pow(y2 - y3, 2) + pow(z2 - z3, 2), 3.0L/2.0L)) + 2*(x1 - x2)/(sqrt(pow(x1 - x2, 2) + pow(y1 - y2, 2) + pow(z1 - z2, 2))*sqrt(pow(x2 - x3, 2) + pow(y2 - y3, 2) + pow(z2 - z3, 2))) - 2*(x2 - x3)*pow(y1 - y2, 2)/(pow(pow(x1 - x2, 2) + pow(y1 - y2, 2) + pow(z1 - z2, 2), 3.0L/2.0L)*sqrt(pow(x2 - x3, 2) + pow(y2 - y3, 2) + pow(z2 - z3, 2))) - 2*(x2 - x3)*(y1 - y2)*(-y2 + y3)/(sqrt(pow(x1 - x2, 2) + pow(y1 - y2, 2) + pow(z1 - z2, 2))*pow(pow(x2 - x3, 2) + pow(y2 - y3, 2) + pow(z2 - z3, 2), 3.0L/2.0L)) + 2*(x2 - x3)/(sqrt(pow(x1 - x2, 2) + pow(y1 - y2, 2) + pow(z1 - z2, 2))*sqrt(pow(x2 - x3, 2) + pow(y2 - y3, 2) + pow(z2 - z3, 2)))) + (-(x1 - x2)*(z2 - z3)/(sqrt(pow(x1 - x2, 2) + pow(y1 - y2, 2) + pow(z1 - z2, 2))*sqrt(pow(x2 - x3, 2) + pow(y2 - y3, 2) + pow(z2 - z3, 2))) + (x2 - x3)*(z1 - z2)/(sqrt(pow(x1 - x2, 2) + pow(y1 - y2, 2) + pow(z1 - z2, 2))*sqrt(pow(x2 - x3, 2) + pow(y2 - y3, 2) + pow(z2 - z3, 2))))*(-2*(x1 - x2)*(y1 - y2)*(z2 - z3)/(pow(pow(x1 - x2, 2) + pow(y1 - y2, 2) + pow(z1 - z2, 2), 3.0L/2.0L)*sqrt(pow(x2 - x3, 2) + pow(y2 - y3, 2) + pow(z2 - z3, 2))) - 2*(x1 - x2)*(-y2 + y3)*(z2 - z3)/(sqrt(pow(x1 - x2, 2) + pow(y1 - y2, 2) + pow(z1 - z2, 2))*pow(pow(x2 - x3, 2) + pow(y2 - y3, 2) + pow(z2 - z3, 2), 3.0L/2.0L)) + 2*(x2 - x3)*(y1 - y2)*(z1 - z2)/(pow(pow(x1 - x2, 2) + pow(y1 - y2, 2) + pow(z1 - z2, 2), 3.0L/2.0L)*sqrt(pow(x2 - x3, 2) + pow(y2 - y3, 2) + pow(z2 - z3, 2))) + 2*(x2 - x3)*(-y2 + y3)*(z1 - z2)/(sqrt(pow(x1 - x2, 2) + pow(y1 - y2, 2) + pow(z1 - z2, 2))*pow(pow(x2 - x3, 2) + pow(y2 - y3, 2) + pow(z2 - z3, 2), 3.0L/2.0L))) + ((y1 - y2)*(z2 - z3)/(sqrt(pow(x1 - x2, 2) + pow(y1 - y2, 2) + pow(z1 - z2, 2))*sqrt(pow(x2 - x3, 2) + pow(y2 - y3, 2) + pow(z2 - z3, 2))) - (y2 - y3)*(z1 - z2)/(sqrt(pow(x1 - x2, 2) + pow(y1 - y2, 2) + pow(z1 - z2, 2))*sqrt(pow(x2 - x3, 2) + pow(y2 - y3, 2) + pow(z2 - z3, 2))))*(2*pow(y1 - y2, 2)*(z2 - z3)/(pow(pow(x1 - x2, 2) + pow(y1 - y2, 2) + pow(z1 - z2, 2), 3.0L/2.0L)*sqrt(pow(x2 - x3, 2) + pow(y2 - y3, 2) + pow(z2 - z3, 2))) + 2*(y1 - y2)*(-y2 + y3)*(z2 - z3)/(sqrt(pow(x1 - x2, 2) + pow(y1 - y2, 2) + pow(z1 - z2, 2))*pow(pow(x2 - x3, 2) + pow(y2 - y3, 2) + pow(z2 - z3, 2), 3.0L/2.0L)) - 2*(y1 - y2)*(y2 - y3)*(z1 - z2)/(pow(pow(x1 - x2, 2) + pow(y1 - y2, 2) + pow(z1 - z2, 2), 3.0L/2.0L)*sqrt(pow(x2 - x3, 2) + pow(y2 - y3, 2) + pow(z2 - z3, 2))) - 2*(-y2 + y3)*(y2 - y3)*(z1 - z2)/(sqrt(pow(x1 - x2, 2) + pow(y1 - y2, 2) + pow(z1 - z2, 2))*pow(pow(x2 - x3, 2) + pow(y2 - y3, 2) + pow(z2 - z3, 2), 3.0L/2.0L)) - 2*(z1 - z2)/(sqrt(pow(x1 - x2, 2) + pow(y1 - y2, 2) + pow(z1 - z2, 2))*sqrt(pow(x2 - x3, 2) + pow(y2 - y3, 2) + pow(z2 - z3, 2))) - 2*(z2 - z3)/(sqrt(pow(x1 - x2, 2) + pow(y1 - y2, 2) + pow(z1 - z2, 2))*sqrt(pow(x2 - x3, 2) + pow(y2 - y3, 2) + pow(z2 - z3, 2)))); }

There are several multiple occurrences of sqrts and pows ​​of the same numbers that could be calculated once to improve readability and runtime. But I dont know how...

Q1: Do you know how to make sympy do this automatically?

Q2: Do you know how to process this code with another tool?

Q3: Can gcc optimize this at compile time? Why?

+6
source share
2 answers

gcc will probably optimize this, but if you want to do it yourself, take a look at cse . http://docs.sympy.org/latest/modules/simplify/simplify.html#module-sympy.simplify.cse_main

+3
source

Here is my own little script based on the asmeurers hint:

 def sympyToC( symname, symfunc ): tmpsyms = numbered_symbols("tmp") symbols, simple = cse(symfunc, symbols=tmpsyms) symbolslist = map(lambda x:str(x), list(symfunc.atoms(Symbol)) ) symbolslist.sort() varstring=",".join( " double "+x for x in symbolslist ) c_code = "double "+str(symname)+"("+varstring+" )\n" c_code += "{\n" for s in symbols: #print s c_code += " double " +ccode(s[0]) + " = " + ccode(s[1]) + ";\n" c_code += " r = " + ccode(simple[0])+";\n" c_code += " return r;\n" c_code += "}\n" return c_code

And for python3. 5+:

 def sympyToC( symname, symfunc ): tmpsyms = numbered_symbols("tmp") symbols, simple = cse(symfunc, symbols=tmpsyms) symbolslist = sorted(map(lambda x:str(x), list(symfunc.atoms(Symbol)))) varstring=",".join( " double "+x for x in symbolslist ) c_code = "double "+str(symname)+"("+varstring+" )\n" c_code += "{\n" for s in symbols: #print s c_code += " double " +ccode(s[0]) + " = " + ccode(s[1]) + ";\n" c_code += " r = " + ccode(simple[0])+";\n" c_code += " return r;\n" c_code += "}\n" return c_code
+3
source

More articles:


All Articles