JuliaDynamics / DynamicalSystems.jl

Award winning software library for nonlinear dynamics and nonlinear timeseries analysis
https://juliadynamics.github.io/DynamicalSystemsDocs.jl/dynamicalsystems/dev/
Other
838 stars 93 forks source link

Decide on what package will be used for Symbolic computations #1

Closed Datseris closed 7 years ago

Datseris commented 7 years ago

There is SymEngine.jl, SymPy and ParameterizedFunctions (which relies on SymEngine but defines a convenient interface for an ode problem and it's jacobian)

ChrisRackauckas commented 7 years ago

Relevant information is here:

https://github.com/Datseris/DynamicalSystems.jl/issues/2#issuecomment-305846620

ParameterizedFunctions uses SymEngine under the hood because it's much much faster than SymPy, and while it's less complete of a package, it has everything that's needed for calculating system Jacobians. So I definitely think that SymEngine > SymPy for this application.

Whether to use ParameterizedFunctions or directly use SymEngine? Depends. If you're just getting the Jacobian, ParameterizedFunctions.jl computes that for you, and FYI you can get the expression that it creates via f.Jex from the returned "function". The Jacobian it returns is accessed via f(Val{:jac},t,u,J), which you can just put a closure over to match your format. If that's precisely what you need, then you also get the benefit that you're building off a well-established DSL, so letting other users give their own ParameterizedFunction is an option. Also, it does have the option of calculating a lot more, including Hessians, in case you ever need more.

That said, using SymEngine directly would give more control. I don't think it would be any faster than ParameterizedFunctions (because it would likely turn out to be very very similar code), but if you need to do some odd symbolic calculation you can. I'm not entirely sure what else in the realm of dynamical systems you would need though (I built it to include the things needed for stochastic dynamical systems).

If there is anything extra you need, you may want to put your suggestions into DiffEqDSL.jl. I am building that DSL to handle the specification and symbolic calculations for both continuous and discrete dynamics, since mixing discrete stochastic (Gillespie) simulations and continuous equations (SDEs, PDEs, ODEs) is what I am currently doing a lot of work in. If there's something you'd find useful there, feel free to put in a feature request. That's still a ways off, and ParameterizedFunctions is what we have for now.

https://github.com/JuliaDiffEq/DiffEqDSL.jl/issues/1

In any case, symbolic computing is limited because it requires the actual function definition and it requires that the function satisfies a sufficiently restrained DSL that it can do computations on. This means it really does need to be implemented as a macro, or you need the user to pass expressions (which isn't that nice of an API: that tends to confuse most people). There are ways to grab a function definition from an anonymous function (and this is what is used in Gallium.jl), but it's undocumented and changes every Julia version, so it's inherently fragile. In fact, if you look at this issue:

https://github.com/Keno/Gallium.jl/issues/153

LambdaInfo is a field inside of a field for "standard data types" which would contain the code that defined the function. Gallium uses this to analyze the function in some way, but the fact that this changed --> Julia v0.6 is a major reason why Gallium has been incompatible with v0.6 for all of these months. I think that's sufficient to show that hacking away at the undocumented function internals to grab the code is not a good idea unless it's absolutely necessary.

Autodifferentiation via ForwardDiff is less restrained and except for stiff equations with hyper-symbolic optimizations, it will do almost as well. So #2 discusses the alternative there, which is likely a better idea for this use case.

Datseris commented 7 years ago

Hey Chris. thanks for taking the effort to write down so much stuff.

As you said at the end of your message, I will actually completely avoid Symbolic computation. The reason is that if one uses SVectors, the automatic jacobians is almost as fast as the user provided ones (at maximum up to a factor of 2). This mind-blowing fact simply led me to restructure completely how discrete systems are written. (See the examles in the benchmark folder, and specifically this one: jrevels_more.jl This one really blew my mind!!!

So I am closing this one. If at some point in the future symbolic computation is really necessary, I guess I will have to learn SymEngine.