[SciPy-user] Calculating symbolic Jacobians, was [OT] Re: Pros and Cons of Python verses other array environments
rclewley at cam.cornell.edu
Thu Sep 28 18:34:37 CDT 2006
AFAIK sympy does not support multivariable expressions or user-defined
functions. These are useful in constructing Jacobians in the most
syntactically neat way.
But if you need something simple and in pure python right away, our
PyDSTool project contains a modestly-featured symbolic toolbox in pure
python, including functions to create Jacobians very easily. You don't
have to keep the whole baggage of our project to use the module
Symbolic.py. Although it currently only works with "old" SciPy/Numeric, we
are working on migrating to NumPy, etc. right now.
Also, it's not as professionally developed as swiginac/ginac and
undoubtedly not as fast, but it has worked well for me on moderately-sized
e.g. for a vector field defined in a test example of SciPy's VODE
>>> ydot0=Fun(-0.04*y0 + 1e4*y1*y2, [y0, y1, y2], 'ydot0')
>>> ydot2=Fun(3e7*y1*y1, [y0, y1, y2], 'ydot2')
>>> ydot1=Fun(-ydot0(y0,y1,y2)-ydot2(y0,y1,y2), [y0, y1, y2], 'ydot1')
>>> F = Fun([ydot0(y0,y1,y2),ydot1(y0,y1,y2),ydot2(y0,y1,y2)], [y0,y1,y2], 'F')
# Diff returns a symbolic object so would like to turn back into a
# function of y0, y1, y2
>>> jac=Fun(Diff(F,[y0,y1,y2]), [y0, y1, y2], 'Jacobian')
# Calls to this function return a symbolic object by default
>>> jac(0.1, 0.3, 0.5)
QuantSpec Jacobian (ExpFuncSpec)
# so use
>>> print jac(0.1, 0.3, 0.5)
# or use .tonumeric() because there are no free names left in the
# resulting expression
>>> jac(0.1, 0.3, 0.5).tonumeric()
array([[ -4.00000000e-02, 5.00000000e+03, 3.00000000e+03],
[ 4.00000000e-02, -1.80050000e+07, -3.00000000e+03],
[ 0.00000000e+00, 1.80000000e+07, 0.00000000e+00]])
# The beauty of this is being able to do use other symbols in the call,
# substitutions, etc.
>>> j = jac(x, x, 0.5)
>>> print j
>>> jsubs = j.eval(x=10)
array([[ -4.00000000e-02, 5.00000000e+03, 1.00000000e+05],
[ 4.00000000e-02, -6.00005000e+08, -1.00000000e+05],
[ 0.00000000e+00, 6.00000000e+08, 0.00000000e+00]])
Automatic simplification of the resulting float-only sub-expressions (like
"10000*0.5") isn't fully working yet for these "symbol arrays", but it
does work in the scalar case. Anyway, check out the wiki page Symbolic
at pydstool.sourceforge.net/ if you want more examples and documentation.
Hope this is of interest!
On Thu, 28 Sep 2006, Gael Varoquaux wrote:
> On Thu, Sep 28, 2006 at 04:06:23PM -0400, A. M. Archibald wrote:
>> As for missing features, I'd like to see even very basic symbolic
>> calculation tools, so that (for example) I could just call a symbolic
>> derivative function to supply a Jacobian to the minimizers. From a bit
>> of web research, swiginac seems like the most promising alternative.
> You can have a look at sympy too. It is very young, though.
More information about the SciPy-user