I'm struggling with the fact that elements of sympy.MatrixSymbol don't seem to interact well with sympy's differentiation routines.
The fact that I'm trying to work with elements of sympy.MatrixSymbol rather than "normal" sympy symbols is beacause I want to autowrap a large function, and this seems that this is the only way to overcome argument limitations and enable input of a single array.
To give the reader a picture of the restrictions on possible solutions, I'll start with an overview of my intentions; however, the hasty reader might as well jump to the codeblocks below, which illustrate my problem.
Declare a vector or array of variables of some sort.
Build some expressions out of the elements of the former; these expressions are to make up the components of a vector valued function of said vector. In addition to this function, I'd like to obtain the Jacobian w.r.t. the vector.
Use autowrap (with the cython backend) to get numerical implementations of the vector function and its Jacobian. This puts some limitations on the former steps: (a) it is desired that the input of the function is given as a vector, rather than a list of symbols. (Both because there seems to be a limit to the number of inputs for an autwrapped function, and to ease interaction with scipy later on, i.e. avoid having to unpack numpy vectors to lists often).
On my journey, I ran into 2 issues:
- Cython does not seem to like some sympy function, among them
sympy.Max, upon which I heavily rely. The "helpers" kwarg of autowrap seems unable to handle multiple helpers at once. - This is by itself not a big deal, as I learned to circumvent it using abs() or sign(), which cython readily understands.
(see also this question on the above)
- As stated before, autowrap/cython do not accept more than 509 arguments in form of symbols, at least not in my compiler setup. (See also here)
As I would prefer to give a vector rather than a list as input to the function anyways, I looked for a way to get the wrapped function to take a numpy array as input (comparable to DeferredVector + lambdify). It seems the natural way to do this is sympy.MatrixSymbol. (See thread linked above. I'm not sure there'd be an alternative, if so, suggestions are welcome.)
My latest problem then starts here: I realized that the elements of sympy.MatrixSymbol in many ways do not behave like "other" sympy symbols. One has to assign the properties real and commutative individually, which then seems to work fine though. However, my real trouble starts when trying to get the Jacobian; sympy seems not to get derivatives of the elements right out of the box:
import sympy
X= sympy.MatrixSymbol("X",10,1)
for element in X:
element._assumptions.update({"real":True, "commutative":True})
X[0].diff(X[0])
Out[2]: Derivative(X[0, 0], X[0, 0])
X[1].diff(X[0])
Out[15]: Derivative(X[1, 0], X[0, 0])
The following block is a minimal example of what I'd like to do, but here using normal symbols: (I think it captures all I need, if I forgot something I'll add that later.)
import sympy
from sympy.utilities.autowrap import autowrap
X = sympy.symbols("X:2", real = True)
expr0 = X[1]*( (X[0] - abs(X[0]) ) /2)**2
expr1 = X[0]*( (X[1] - abs(X[1]) ) /2)**2
F = sympy.Matrix([expr0, expr1])
J = F.jacobian([X[0],X[1]])
J_num = autowrap(J, args = [X[0],X[1]], backend="cython")
And here is my (currently) best guess using sympy.MatrixSymbol, which then of course fails because the Derivative-expressions within J:
X= sympy.MatrixSymbol("X",2,1)
for element in X:
element._assumptions.update({"real":True, "commutative":True, "complex":False})
expr0 = X[1]*( (X[0] - abs(X[0]) ) /2)**2
expr1 = X[0]*( (X[1] - abs(X[1]) ) /2)**2
F = sympy.Matrix([expr0, expr1])
J = F.jacobian([X[0],X[1]])
J_num = autowrap(J, args = [X], backend="cython")
Here is what Jlooks like after running the above:
J
Out[50]:
Matrix([
[(1 - Derivative(X[0, 0], X[0, 0])*X[0, 0]/Abs(X[0, 0]))*(-Abs(X[0, 0])/2 + X[0, 0]/2)*X[1, 0], (-Abs(X[0, 0])/2 + X[0, 0]/2)**2],
[(-Abs(X[1, 0])/2 + X[1, 0]/2)**2, (1 - Derivative(X[1, 0], X[1, 0])*X[1, 0]/Abs(X[1, 0]))*(-Abs(X[1, 0])/2 + X[1, 0]/2)*X[0, 0]]])
Which, unsurprisingly, autowrap does not like:
[...]
wrapped_code_2.c(4): warning C4013: 'Derivative' undefined; assuming extern returning int
[...]
wrapped_code_2.obj : error LNK2001: unresolved external symbol Derivative
How can I tell sympy that X[0].diff(X[0])=1 and X[0].diff(X[1])=0? And perhaps even that abs(X[0]).diff(X[0]) = sign(X[0]).
Or is there any way around using sympy.MatrixSymbol and still get a cythonized function, where the input is a single vector rather than a list of symbols?
Would be greatful for any input, might well be a workaround at any step of the process described above. Thanks for reading!
Edit:
One short remark: One solution I could come up with myself is this:
Construct F and J using normal symbols; then replace the symbols in both expressions by the elements of some sympy.MatrixSymbol. This seems to get the job done, but the replacement takes considerable time, as J can reach dimensions of ~1000x1000 and above. I therefore would prefer to avoid such an approach.
