Derivatives about modelingtoolkit.jl HOT 14 CLOSED

Not too clutter up the ideas here, but I wanted to make sure that @sglyon and @albop knew what you are up to, as they work on DSLs.

In the longrun, it sure would be nice to have be able to handle difference equations (i.e. a difference rather than a derivative operator) in this sort of DSL for my dream Dolo2.jl which ditches python :-)

from modelingtoolkit.jl.

Hi there! Interesting project indeed and I'm curious about it.
@jlperla : I don't know what you are talking about. Dolo.jl doesn't depend on python. It depends on a serialization file in yaml though. And has some algorithm in common with dolo.py. But you are not suggesting to remove dependence on math are you ? ;-) I'll create an issue to provide an option to write models in julia or with a dsl and assign it to you :-)

from modelingtoolkit.jl.

This isn't the right place to have the discussion, but I will do it anyways!

• What you have written in Dolo is a language agnostic domain specific language (DSL) specified in yaml. You then have an interpreter written in python, julia, or whatever. It is along the same lines as AMPL or Dynare.
• Then anytime you want to add features, you need to add the features first to your language and then to the underlying implementations in julia, python. If something is easy in Julia but not python? Tough. You want to add in some new monte-carlo integration or GPU support or high precision floating points? Add it to the yaml specification and implement it in julia and python. Want to add in a new prior distribution for simulating, or support a new commerical optimizer type? Implement it yourself.
• Hence, python is important because as long as you support it, you need to be held back by it as the lowest common denominator for available features.

What I am suggesting is a radical redesign: build an EDSL (embedded DSL) in Julia itself. Think along the lines of https://github.com/JuliaApproximation/ApproxFun.jl or https://github.com/JuliaOpt/JuMP.jl or https://github.com/yebai/Turing.jl/wiki/Get-started or https://github.com/JuliaStats/GLM.jl

from modelingtoolkit.jl.

Why not have that discussion in the DSL thread?

SciML/DifferentialEquations.jl#261

This is just about derivatives...

from modelingtoolkit.jl.

https://github.com/JuliaDiff/DiffRules.jl might be useful.

It would be nice if this was generally extensible.

from modelingtoolkit.jl.

It would be nice if this was generally extensible.

I think that's the key property IMO. I had an implementation which didn't have this so I scrapped it. I think it can be done by allowing the definition of a dispatch for a function. I almost have it.

from modelingtoolkit.jl.

I put in a basic derivative expansion API but right now it requires a weird dispatch and doesn't take into account the chain rule. We likely need something slightly different here but I don't know what.

from modelingtoolkit.jl.

If the purpose is to calculate Jacobians, Hessians, etc., why not do this at a later stage by feeding Julia expressions to Calculus.jl or similar? That is, before the construction of the ODEFunction we will have expressions for the rhs of each equation, so we can just generate the corresponding derivatives, no?

from modelingtoolkit.jl.

I think it would be harder to write functions to feed everything to Calculus.jl and then convert back to our AST, rather than to just write a simple mechanism on our AST. We want to do things like index reduction of DAEs, conversion to first-order ODEs, etc. directly on our ASTs and it would be bad in the long-run to hack it onto Calculus.jl (which really doesn't offer much here: one page of derivative rules? We can just copy those over).

from modelingtoolkit.jl.

OK. A dual number system (storing Variables and their partials) would be relatively easy if we have a list of derivative rules, since most functions and operators are already defined, but we would need to interpret of the AST (i.e. using Variable as values).

from modelingtoolkit.jl.

Yeah, that's pretty much what this is doing right now. I just need to find out how to make it do the chain rule, and how to make it easy for users to declare derivatives of their function by some variable (i.e. d sin(x) / dt we need to have them say it should be cos(x)*dx/dt) . I'm not sure about the best way to do that when you scale to multiple arguments.

from modelingtoolkit.jl.

But with a dual number you do not need to worry about chain rule. For example, + and sin would be something like

+(x::Dual, y::Dual) = Dual(x.val + y.val, x.partials .+ y.partials)
sin(x::Dual) = Dual(sin(x.val), cos(x.val).*x.partials)

So in the end the derivatives are stored in partials. This allows calculating Jacobians and sensitivity equations, but maybe d sin(x)/dt needs a different approach, I am not sure.

from modelingtoolkit.jl.

That's just different and for different applications. Dual numbers actually just do the chain rule in the Dual part. But the issue with this is it doesn't build functions. We want the ASTs to do further transformations on. For example, take a 4th order ODE. What's the inverse of the Jacobian of the transformation of the 4th order ODE to a system of 1st order ODEs? You can generate functions of functions to do this, but then it'll be doing all of the calculations at runtime, and it'll have to do the full calculation at every step.

If we instead build the AST, we can build functions that directly compute the result. This is of course much much faster which is why @ode_def does well.

Of course there's a tradeoff though. Autodiff does scale better to large systems of variables. One way that we avoid the problems here is by allowing our variables to be vectors, and thus we aren't writing out the AST for every component but instead doing the vector calculus. This helps manage the explosion for many sizable models with structure. In addition, with ParameterizedFunctions.jl we saw that this scaled very well with systems that had sparse Jacobians which is quite common in diffeq and nonlinear solving applications.

When the systems get large enough we would need to make use of autodiff. That simply means that some transformations won't be possible, or we just add a separate transform function that does this via autodiff. Additionally, users could tell it to avoid differentiating complicated functions. If you do something like:

@register complicated_function(x)
d_complicated_function(x) = ForwardDiff.derivative(complicated_function,x)
@register d_complicated_function(x)
@derivative complicated_function(x) d_complicated_function(x)

would tell it to not try differentiating past complicated_function and instead use autodiff on it (and we can make that simpler with a macro of course). This way we can more precisely apply autodiff to smaller functions to make the calculations even cheaper, even when it's required. It also lets users say when 1-D numerical diff could be useful (and our simplification routines can then make it done only once a function).

from modelingtoolkit.jl.

Basic derivatives with automatic chain rule is complete.

from modelingtoolkit.jl.