Type-Based DSL

[ChrisRackauckas]

I think it's time that we put down a design for a type-based DSL which can specify ODEs, higher order ODEs, SDEs, regular jumps, DAEs, DDEs, PDEs, and all combinations. I think we actually get a DSL for linear and nonlinear solve problems for free. It sounds terrifying, but the specification doesn't seem all that bad after giving it a try.

I think I can just explain by examples. First there's the simple example for the fully type-based form:

# Define some variables
x = DependentVariable(:x)
y = DependentVariable(:y)
z = DependentVariable(:z)
t = IndependentVariable(:t)
D = Derivative(t) # Default of first derivative, Derivative(t,1)
σ = Parameter(:σ)
ρ = Parameter(:ρ)
β = Parameter(:β)

# Define some expressions
eqs = [D*x == σ*(y-x)
       D*y == x*(ρ-z)-y
       D*z == x*y - β*z]

# Define a differential equation
de = DiffEq(eqs)

# Now go to the DiffEq pipeline
prob = ODEProblem(de,u0,tspan,p)

Now we start adding sugar:

@DV x y z
@IV t
@Parameter σ,ρ,β
D = Derivative(t)
de = @diffeq begin
  D*x = σ*(y-x)
  D*y = x*(ρ-z)-y
  D*z = x*y - β*z
end
prob = ODEProblem(de,u0,tspan,p)

Now we start adding features.

# Array Variables
x = DependentVariable(:x,1:10) # Array of variables
@DV x[1:10]
# Can do the same on IndependentVariable, Parameters, ...

# Noise
dW = NoiseVariable(:dW,3)

# Regular Jump
dp = JumpVariable(5x) # rate(u,p,t) = 5u[1]

# Delayed Variable
y_τ  = y(t-5)
y_τ2 = y(t-f)

# Mass matrix DAE
de = @diffeq begin
  D*x + D*y = σ*(y-x)
  D*y = x*(ρ-z)-y
  D*z = x*y - β*z
end

# Fully Implicit SDAE
de = @diffeq begin
  exp(D*x + D*y) = σ*(y-x) + dW[1]
  D*y = x*(ρ-z)-y + dW[2]
  D*z = x*y - β*z + dW_3 # Same as dW[3]
end

# Define explicit orderings for variables, parameters, etc.
de = DiffEq(eqs,dep_vars,indep_vars,parameters,noises)

Now we start specifying some PDEs:

@DV u
@IV x t
D_t = Derivative(t)
D_x = Derivative(x)
D_xx = D_x*D_x
D_xx == Derivative(x,2)
D_t*u == D_xx*u

# Assume μ(u,p,t,x) <: Function and σ(u,p,t,x) <: Function are defined
L = Operator(mu*D_x + σ*D_xx)
D_t*u = L*u

And as an added bonus

# Nonlinear rootfinding problem
nl = @nl begin
  σ*(y-x)
  x*(ρ-z)-y
  x*y - β*z
end

nlprob = NonlinearProblem(nl,u0)
# Now send this over to NLsolve with pre-defined Jacobians!

# Similar for linear

# Possible extension down the line
nl = @nl begin
  σ*(y-x) + x*(ρ-z)-y + x*y - β*z
end
scalar_function = Scalar(nl) # Defines and computes Jacobians/Hessians, send to Optim!

Extra Features

The purpose of this is not to do a full physical modeling setup (though it should be possible to build one off of it), but it's to make it easy to specify/modify/optimize differential equations in a natural form that's both computer and human readable. Also, if Modia.jl is going to come out then I don't see a reason to directly build something like that, so I don't plan to spend much time in feature areas that they focus on (but of course those extensions can be readily added). But, there's definitely a few things from Modelica that we can put in here.

One thing is the possibility of initial conditions and parameter values in the model specification.

x = DependentVariable(:x,2.0)
@DV x 2.0 y 3.0 z 4.0

σ = Parameter(:σ,3.0)

Then u0 and p wouldn't be needed in the ODEProblem construction.

Model composition can come later. It could be something like:

struct Lorenz <: DiffEqComponent
  x
  y
  z::Flow
end

eqs = [D*x == σ*(y-x)
       D*y == x*(ρ-z)-y
       D*z == x*y - β*z]
lorenz = Component(Lorenz,eqs)

or

@component struct Lorenz
  x
  y
  z::Flow
end begin
  D*x = σ*(y-x)
  D*y = x*(ρ-z)-y
  D*z = x*y - β*z
end

and then have some kind of connect. This is stuff that can be added later for people to build modeling frameworks, but it's not truly essential to the original goal.

Implementation details

Most of the implementation details are pretty obvious. It's a bunch of expressions with small types. But a few things to mention. Variables, parameters, etc. have direct conversions to SymEngine.Basics. This is then makes it easy to do math on them.

The "missing" implementation detail is that any arithmetic operation forms a DEStatement which holds the components by their types by their addition (differential statements, jump statements, noise statements, base statements). == then generates a DEquation with a statement on the lhs and a statement on the rhs.

Remaining Issues

• How do you specify history functions (for DDEs)?

Deprecation Path

There won't truly be a "deprecation" of @ode_def since it is still nice. Instead, it would just output a DiffEq when this is ready, so it would be directly incorporated into the next DSL.

Remarks

Based off of SciML/DiffEqDSL.jl#1

Pinging @tshort @gabrielgellner @zahachtah @mauro3 due to how helpful you all were in round one!
@korsbo and @Gaussia probably have great ideas as well. @vjd and @simonbyrne

[gabrielgellner]

This looks absolutely amazing to me. Don't have a feel for some of the higher level stuff, but will play around with some of my problems to see how it feels. I have a mild worry about having @ode_def around as well, as it begins to have a lot of different ways to specify a problem which can be a burden for documentation / teaching.

[ChrisRackauckas]

I have a mild worry about having @ode_def around as well, as it begins to have a lot of different ways to specify a problem which can be a burden for documentation / teaching.

It'll get a deprecation warning once the upgrade path is stabilized. What I meant by the "no dep" is that there's no rush for it. It'll slowly be pushed out of the docs though.

[mauro3]

Yes, this looks cool! This has Fenics-y feel to it, presumably you took some inspiration from there too? Maybe worth comparing.

Key is (unlike Fenics), and I think that is your philosophy/plan, that the solvers keep on working with just normal functions. That way one can hook into the diff-eq system at any level.

[korsbo]

I like the idea, and I think that it would be relatively easy (although laborious) to implement. There are a few points that I would like to raise.

Local scoping

One thing that makes me hesitant to your specific suggestion is the separation of the model and variable/parameter scope.

Imagine the following scenario:

I start my work by defining a model

@DV x y z
@IV t
@Parameter a b c
D = Derivative(t)
model = @diffeq begin
    Dx = a - b*x 
    Dy = c*x/z - y
    Dz = y - z
end

I then realize that the dynamics of x is relatively uninteresting and I can replace it with a parameter x = a/b instead. So this prompts me to write:

@DV y z
@IV t
@Parameter x c
D = Derivative(t)
model2 = @diffeq begin
    Dy = c*x/z - y
    Dz = y - z
end

However, since I did not restart my kernel, I now have x defined as both a variable and a parameter. Madness ensues.

If we instead let the parameter and variable definition reside inside the @diffeq macro, we could avoid this.

Type specification

In my work, I often compare a few different models in any given project. For me, it is a massive benefit to be able to create specific methods for each of them. With @ode_def I name my own types and can then for example decide if my equilibrate function can be done analytically for some models while it is done numerically for others. This control over my type and its name is brilliant, let's keep that.

We could also expand it! I often find that my models can be divides into groups which are all able to share some methods. It would be nice if I as a user could assign an abstract supertype to them.

myModel = @diffeq MyType <: MyAbstractSupertype begin
...
end

Here, MyAbstractSupertype should automatically inherit from the common AbstractDifferentialEquation type that everything created by this macro would inherit from. And while we are at it, we could allow this in a chained fashion: MyType <: MyAbstractType <:MySecondAbstractType.

Knowing what you are going to get.

As I see it, the prototype DSL would infer the length of the parameter vector (and their relative position). I would prefer this to be explicitly done by the user. Just like the parameter definitions of @ode_def. I even like (and sometimes use) the ability to specify more parameters of my DE than I actually use.

Allowing "multi-cellular" input.

I would love to be able to use this DSL to define a reaction-diffusion model with ODEs. Could we think about the ability to supply nested or multi-dimensional u0 to the problem and just have the DSL do the right thing? Could we also provide diffusion reactions, and possibly some topology type to specify neighbourhood relations? This has some overlaps with PDEs, but I would like to be able to supply my own discretization and my own graph. Vertex based models would also be cool (although a fairly tall order).

Anyway. Loving this new DSL idea, this could be really cool!

[ChrisRackauckas]

Yes, this looks cool! This has Fenics-y feel to it, presumably you took some inspiration from there too? Maybe worth comparing.

Yes, there was inspiration pulled from FEniCS, XPPAUT, JuMP, Maple, Mathematica, and Modelica. And yeah, they key is to have this be an abstract definition of a differential equation for many different purposes (simplification and optimization, PDE operator specification, Latexify, etc.) and as a generative tool for the DE solvers, but not have things reliant on the use of the DSL. In the end, some things will never be nice to do via a DSL so it's good to keep them separate.

However, since I did not restart my kernel, I now have x defined as both a variable and a parameter. Madness ensues.

No that's fine, because x would just be a variable of a pure type.

struct DependentVariable{name} end
Base.@pure DependentVariable(name) = DependentVariable{name}()
struct Parameter{name} end
Base.@pure Parameter(name) = Parameter{name}()

x = DependentVariable(:x)
x = Parameter(:x)

is what that would come out to. That's not defining any new types (just parameterizing) and so it's fine.

Type specification

I didn't think about this part, but it'll be a little different since we do want to stay away from type generation due to the redefinition issues (and this would tie it to macro use). However, what we could do is in the DiffEq constructor allow a type parameter:

de = DiffEq{MyType}(eqs)

so that it's an AbstractDiffEq{MyType}. This would allow you to dispatch how you want. Then the syntax:

struct MyType <: MySuperType end
myModel = @diffeq MyType begin
...
end

would boil down to creating a DiffEq{MyType,...} <: AbstractDiffEq{MyType} which gives you what you need to dispatch. And we can easily make that optional, so if it's not used it's just Void.

As I see it, the prototype DSL would infer the length of the parameter vector (and their relative position). I would prefer this to be explicitly done by the user.

This is what I meant by:

# Define explicit orderings for variables, parameters, etc.
de = DiffEq(eqs,dep_vars,indep_vars,parameters,noises)

So if you pass in parameters = [e,b,f,c,a], then that's the ordering if uses. Otherwise it tries to determine it automatically. I'm not sure how that interacts with the @diffeq macro way of specifying though: maybe we allow parameters after end just like in @ode_def for this purpose, but make it optional?

I would love to be able to use this DSL to define a reaction-diffusion model with ODEs. Could we think about the ability to supply nested or multi-dimensional u0 to the problem and just have the DSL do the right thing? Could we also provide diffusion reactions, and possibly some topology type to specify neighbourhood relations?

Yesyesyes. We need this. In fact, I have been working on stiff solvers which automatically decompose on reaction-diffusion equations to have a block-diagonal Jacobian which is just the reaction Jacobian. I've been looking to create another ODEProblem type (like SplitODEProblem) that is specifically for this type of data so that way the solvers could make use of this optimization (because this obvious is highly specialized to this kind of "nonlocal linear + local nonlinear" PDE)

I think this is a good way to describe the extensibility part of the DSL. The first idea is the DiffEqFunction:

SciML/DiffEqBase.jl#52

Basically each problem type will be defined by its function type, and its function type is specific to its internal data specification. For example, ODEFunction can take jacobians, inverse-W, etc. SplitODEFunction needs to functions. So then ODEProblem(f,u0,tspan), SplitODEProblem(f1,f2,u0,tspan) just wrap your f into these. This then allows users to pack the functions into a type instead of doing the dispatch+trait specification for the performance addons. There's a few other benefits by handling the data like this for the solvers too, but in the case where the user isn't defining Jacobians and all of that it's invisible.

But then this allows us a step between the DSL and the problem. Basically,

prob = ODEProblem(de,u0,tspan,p)

is generating the problem via a dispatch that essentially does

prob = ODEProblem(ODEFunction(de),u0,tspan,p)

This function transformation is the key piece of code that does the translation from the DSL to the actual function instantiation. A dispatch here would do something like, check to make sure every statement has a single derivative or is an intermediate calculation (but not an equality constraint) and then after verifying that it's really an ODE, building the Julia functions for the f, Jacobian, etc. and returning the ODEFunction.

[ChrisRackauckas]

That last response made me notice that I missed the "intermediate variable" part of this, like:

x = DependentVariable(:x)
y = DependentVariable(:y)
z = Variable(:z)
t = IndependentVariable(:t)
D = Derivative(t)
eqs = [z == x-y
           D*x == z*x
           D*y == z*y]

This fixes the "intermediate calculation problem" we always had with @ode_def just by being able to recognize that you can substitute this in everywhere, yet still build Jacobians etc. which are smart about saving calculations.

[TorkelE]

Just reiterating that this seems like a very useful feature and something that should be done at some point. Will keep my eyes on how things develop. Some quick things:

• Defining Derivatives is something that will be done a lot. since e.g. x' is a correct Julia expression, maybe one could use that somehow to define derivatives? x'' could be the second derivative (That is two ', and not one "). Admittedly this does not specify derivative with respect to what and might be a bit confusing for PDEs, but I see a potential for streamlining notation a little bit.
• At some stage I would love to see bifurcation analysis integrated in such a DSL (although I am uncertain what the current plans for bifurcation analysis in DiffEq is).

[ChrisRackauckas]

Defining Derivatives is something that will be done a lot. since e.g. x' is a correct Julia expression, maybe one could use that somehow to define derivatives?

We can do my favorite Julia hack:

const t = IndependentVariable(:t) # global, in DiffEqDSL, unexported
ctranpose(x::DependentVariable{name}) = Derivative(t)*x # called adjoint on v0.7
eqs = [x' == σ*(y-x)
           y' == x*(ρ-z)-y
           z' == x*y - β*z]

That would then work without a macro. You'd have to be careful to use that t though, so it would be like:

t = DiffEqDSL.t
eqs = [x' == σ*(y-x)*t
           y' == x*(ρ-z)-y
           z' == x*y - β*z]

would be required to use that.

But there's another underlying question of what the purpose of the DSL is. That rides in to this question.

At some stage I would love to see bifurcation analysis integrated in such a DSL (although I am uncertain what the current plans for bifurcation analysis in DiffEq is).

The analysis of the DSL's results, and its input syntax, are distinct from the data structure. What I've outlined is (I think) pretty simple of a syntax, but the main purpose of this is to be a data structure which is extendable and can encode the idea of "what is a differential equation" for automatic manipulation. So it turns out to be not too bad notationally, but that's not its true goal in this form.

What is nice though is that, with a well-defined data structure for the DSL, you can target it for simplified languages. For example, @ode_def can just target this as its output. Whether streamlined versions that are complimentary and target the same output should exist is another question (and it seems @gabrielgellner is against that), but with this setup they can always be added. One open question for me is how streamlined @diffeq should be. It's supposed to be the "simpler version", but there's always a tension between feature-richness and simplicity.

But because this output is just a bundle of diffeq information, bifurcation analysis can be performed on this data without being directly integrated into the DSL. The hard part of bifurcation analysis is the continuation methods (every bifurcation needs special handling for them to be robust...). But for example right now we have interfacing to PyDSTool.

http://docs.juliadiffeq.org/latest/analysis/bifurcation.html

This setup has the information required to build the input to PyDSTool (it wants a string) and use that. It also has the information to build functions for XPPAUT, and then going further to weirdness, Stan, FEniCS, etc.

Bringing @AleMorales into this because he has experienced and expressed interest.

[AleMorales]

This looks really exciting! And I like that Variable is also there, big models really need it. And default parameter values are also really useful for big models.

But there is something not clear to me (probably just my ignorance) and I think it could make a big difference: does this approach aim at storing a data structure representing the data flow of the model? Or is it just doing code transformation (and maybe some code generation for Jacobians) and just storing the final functions/parameters, etc as required by the solver? The first approach opens up a lot of possibilities like model reduction, genetic programming, merging models into a single system of equations (rather than composition), being able to write model statements in any order (even when there are intermediate variables), generating documentation for models (i.e. all those equations that end up in supplementary materials), etc. At the moment it seems to me that second approach is being used as the parameters, dependent variables, etc are not being stored into an object model explicitly so part of the semantic information would be lost.

[ChrisRackauckas]

But there is something not clear to me (probably just my ignorance) and I think it could make a big difference: does this approach aim at storing a data structure representing the data flow of the model? Or is it just doing code transformation (and maybe some code generation for Jacobians) and just storing the final functions/parameters, etc as required by the solver? The first approach opens up a lot of possibilities like model reduction, genetic programming, merging models into a single system of equations (rather than composition), being able to write model statements in any order (even when there are intermediate variables), generating documentation for models (i.e. all those equations that end up in supplementary materials), etc. At the moment it seems to me that second approach is being used as the parameters, dependent variables, etc are not being stored into an object model explicitly so part of the semantic information would be lost.

That's a good way to frame the discussion. DiffEq(eqs) builds a type which stores the full data flow of the model. Functions can be written on that structure to manipulate it to other ways. Then the ODEProblem (or more concretely the ODEFunction, and then the SDEFunction etc. for all of the different possible transformations) is the code transformation step. The reason to keep them separate is exactly for the reasons you mention. If we just have the whole model stored, then people can extend this to do whatever transformation they want on model functions, and build tooling to translate the model to whatever output they need.

The proposed type-structure is the data. This can be separated from the syntax by the macros. This is separated from the end ODE function via the translation function. By focusing on an extendable data structure for DiffEq, this then leaves it untied to the input mechanism (one can add macros for allowing syntax from other languages or parsing files), the translation mechanism (writing to Stan or FEniCS), and data transformations (model reductions which take in a DiffEq and spits out another DiffEq).

Essentially it's a DiffEq IR, as in, it's a strongly typed intermediate representation for specifying a differential equation.

[AleMorales]

Ok that sounds great. It would be interesting/useful to then add some basic sanity checks (like some very basic static analysis) when constructing the data structure, to reduce errors emitted from the generated functions (e.g. missing inputs and the like).

I wonder how this approach will deal with arbitrary functions on the rhs? You mention arithmetic operations but I can imagine every math function needs to have a method for Parameter, IndependentVariable, etc for this to work. And perhaps a special type for functors that act as parameters (for time-dependent inputs and the like).

[ChrisRackauckas]

I wonder how this approach will deal with arbitrary functions on the rhs? You mention arithmetic operations but I can imagine every math function needs to have a method for Parameter, IndependentVariable, etc for this to work.

Nope, this is the beauty of Julia. It's just like AD: since every basic arithmetic operation produces an Operation (or whatever we want to call it), many arbitrary functions will just work through this. So for example:

f(x,y) = x+y

For two variable types (any of them) x+y produces an Operation, and so when you pass a Parameter and an IndependentVariable to this you get out an Operation. At the bottom level we just have to have registered functions where if it's registered, then it produces that as the call on the operation. So x+y goes to Operation(+,x,y). Function registration is just adding the dispatch f(x::AbstractVariable,y::AbstractVariable) = Operation(f,x,y), so we can have @register f(x,y) so that way users can register functions to be kept in the call graph, and all of the other functions are automatically decomposed in this mechanism.

This would allow us to define things like exp(x) --> Operation(exp,x) and so it's encoded in the call graph instead of its numerical approximation, but I suspect users might not have to touch that much.
Then one thing which is actually quite simple is that by overloading * you define D*exp(x) = exp(x)*D*x.
Calculus.jl already has a decent list of rules here:

https://github.com/johnmyleswhite/Calculus.jl/blob/master/src/differentiate.jl#L116

Honestly, that's easy enough that we might want to just drop SymEngine and use this. I don't think SymEngine really gives us anything if we calculate derivatives ourselves, and this would make it trivial for users to extend.

We might want to just call this SciCompDSL...

[ChrisRackauckas]

Here is a prototype:

abstract type AbstractVariable <: Real end
struct DependentVariable <: AbstractVariable
    name::Symbol
end
struct IndependentVariable <: AbstractVariable
    name::Symbol
end
struct Variable <: AbstractVariable
    name::Symbol
end
struct Parameter <: AbstractVariable
    name::Symbol
end
struct Differential{V<:AbstractVariable} <: AbstractVariable
    iv::V
end
struct Derivative{V<:AbstractVariable,V2<:AbstractVariable} <: AbstractVariable
    v::V
    iv::V2
end

abstract type AbstractStatement end
const Expression = Union{AbstractVariable,AbstractStatement}

struct Operation{F} <: AbstractStatement
    op::F
    args::Vector{Expression}
end
struct Statement{L<:Expression,R<:Expression} <: AbstractStatement
    lhs::L
    rhs::R
end
struct Equations
    eqs::Vector{Statement}
    ivs::Vector{IndependentVariable}
    dvs::Vector{DependentVariable}
    vs::Vector{Variable}
    ps::Vector{Parameter}
end

Base.:-(x::Expression,y::Expression) = Operation(-,Expression[x,y])
Base.:*(x::Expression,y::Expression) = Operation(*,Expression[x,y])
Base.:(==)(x::Expression,y::Expression) = Statement(x,y)

Base.:*(D::Differential,x::Expression) = Derivative(x,D.iv)

# Define some variables
x = DependentVariable(:x)
y = DependentVariable(:y)
z = DependentVariable(:z)
t = IndependentVariable(:t)
D = Differential(t) # Default of first derivative, Derivative(t,1)
σ = Parameter(:σ)
ρ = Parameter(:ρ)
β = Parameter(:β)

σ*(y-x)
D*x
D*x == σ*(y-x)
D*y == x*(ρ-z)-y

Expression[ρ,z]

ρ-z

# Define some expressions
eqs = [D*x == σ*(y-x)
       D*y == x*(ρ-z)-y
       D*z == x*y - β*z]

# Define a differential equation
de = Equations(eqs,[],[],[],[])

I think there's far too much type information in there. Maybe variable types should be enums. I'm going to play with trying to make this more inferrable.

[ChrisRackauckas]

The fully inferred version isn't much different, but the type signatures grow exponentially...

abstract type AbstractVariable <: Real end
struct DependentVariable <: AbstractVariable
    name::Symbol
end
struct IndependentVariable <: AbstractVariable
    name::Symbol
end
struct Variable <: AbstractVariable
    name::Symbol
end
struct Parameter <: AbstractVariable
    name::Symbol
end
struct Differential{V<:AbstractVariable} <: AbstractVariable
    iv::V
end
struct Derivative{V<:AbstractVariable,V2<:AbstractVariable} <: AbstractVariable
    v::V
    iv::V2
end

abstract type AbstractStatement end
const Expression = Union{AbstractVariable,AbstractStatement}

struct Operation{F,T} <: AbstractStatement
    op::F
    args::T
end
struct Statement{L<:Expression,R<:Expression} <: AbstractStatement
    lhs::L
    rhs::R
end
struct Equations{E}
    eqs::E
    ivs::Vector{IndependentVariable}
    dvs::Vector{DependentVariable}
    vs::Vector{Variable}
    ps::Vector{Parameter}
end

Base.:-(x::Expression,y::Expression) = Operation(-,(x,y))
Base.:*(x::Expression,y::Expression) = Operation(*,(x,y))
Base.:(==)(x::Expression,y::Expression) = Statement(x,y)

Base.:*(D::Differential,x::Expression) = Derivative(x,D.iv)

# Define some variables
x = DependentVariable(:x)
y = DependentVariable(:y)
z = DependentVariable(:z)
t = IndependentVariable(:t)
D = Differential(t) # Default of first derivative, Derivative(t,1)
σ = Parameter(:σ)
ρ = Parameter(:ρ)
β = Parameter(:β)

σ*(y-x)
D*x
D*x == σ*(y-x)
D*y == x*(ρ-z)-y

Expression[ρ,z]

ρ-z

# Define some expressions
eqs = (D*x == σ*(y-x),
       D*y == x*(ρ-z)-y,
       D*z == x*y - β*z)

# Define a differential equation
de = Equations(eqs,IndependentVariable[],DependentVariable[],Variable[],Parameter[])