2

I'm curious about the right way to use the sympy function compilation routines (like autowrap) to create a function that takes array inputs.

With lambdify, I can do either one of these:

import numpy as np
import sympy as sp
x, y, a,b = sp.symbols('x,y,a,b')
f=a*x**2 + b*y**3

A)

L1 = sp.lambdify([x,y,a,b],f)
v1=np.array([1,2])
v2=np.array([3,4])
L1(*np.concatenate([v1,v2]))

35

B)

M1=sp.Matrix([x,y])
M2=sp.Matrix([a,b])
L2 = sp.lambdify([M1,M2],f)
L2(v1,v2)

35

Presumably under the hood (B) does the same thing as (A).

But with autowrap I only have option A as far as I can tell.

from sympy.utilities.autowrap import autowrap
A1 = autowrap(f,args=[x,y,a,b],backend='cython')
A1(*np.concatenate([v1,v2]))

35.0

A2 = autowrap(f,args=[M1,M2],backend='cython')

CodeGenArgumentListError: ("Argument list didn't specify: a, b, x, y ", [InputArgument(a), InputArgument(b), InputArgument(x), InputArgument(y)])

Is there a way to do this? My use case is compiling the dy/dx function for ODE simulation. The function (much much messier than the example above) needs to take an array of x values and an array of parameters. Right now I'm using the code above, with np.concatenate and python list expansion, and if I profile the code, the concatenation actually takes more time than the computation. So I would love to get around it, and have sympy take the last statement and generate some C code that looks like

double f(double* M1,double* M2) {
   return M1[0]*pow(M2[0],2)+M1[1]*pow(M2[1],3)
}

Presumably then the numpy arrays could get passed directly.

Is that possible? Is there a better way to do this? Or am I overthinking/overoptimizing here?

Victor Chubukov
  • 1,345
  • 1
  • 10
  • 18
  • I'm not 100% sure, but this might be a duplicate of https://stackoverflow.com/questions/35430479/convert-sympy-expressions-to-function-of-numpy-arrays. – asmeurer Mar 17 '16 at 18:51
  • You may also be interested in https://github.com/sympy/sympy/pull/10640 – asmeurer Mar 17 '16 at 18:52
  • You're absolutely right, that question is asking the exact same thing. It looks like there is no perfect solution, but I will investigate the examples linked in the github thread. – Victor Chubukov Mar 18 '16 at 03:01

1 Answers1

0

The new scipy codegen tutorial covers a way to do this.

An example I was able to follow and modify is here: http://www.sympy.org/scipy-2017-codegen-tutorial/notebooks/40-chemical-kinetics-cython.html

The overview of the course is here: http://www.sympy.org/scipy-2017-codegen-tutorial/

Victor Chubukov
  • 1,345
  • 1
  • 10
  • 18