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 https://github.com/JuliaDiffEq/DiffEqDSL.jl/issues/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:

https://github.com/JuliaDiffEq/DiffEqBase.jl/issues/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[])

[ChrisRackauckas]

@enum VariableSubType None IndependentVariable DependentVariable Parameter

abstract type AbstractVariable <: Real end
struct Variable <: AbstractVariable
    name::Symbol
    subtype::VariableSubType
end
Variable(name) = Variable(Variable,None)
(x::VariableSubType)(name) = Variable(name,x)

struct Differential <: AbstractVariable
    iv::Variable
end
struct Derivative <: AbstractVariable
    v::Variable
    iv::Variable
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{Variable}
    dvs::Vector{Variable}
    vs::Vector{Variable}
    ps::Vector{Variable}
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

# 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,Variable[],Variable[],Variable[],Variable[])

Using an enum for the subtype of a variable seems to make a lot of sense and now all things are type stable and the type information isn't growing as fast. Just cutting out some type information leads to:

abstract type AbstractVariable <: Real end
abstract type AbstractStatement end

@enum VariableSubType None IndependentVariable DependentVariable Parameter NoiseVariable JumpVariable DelayedDependentVariable

struct Variable <: AbstractVariable
    name::Symbol
    subtype::VariableSubType
end
Variable(name) = Variable(Variable,None)
(x::VariableSubType)(name) = Variable(name,x)

const Expression = Union{AbstractVariable,AbstractStatement}

struct Differential
    x::Variable
end
function Derivative end
Base.:*(D::Differential,x::Expression) = Operation(Derivative,Expression[x,D.x])

struct Operation <: AbstractStatement
    op::Function
    args::Vector{Expression}
end
struct Statement <: AbstractStatement
    lhs::Expression
    rhs::Expression
end
struct Equations
    eqs::Vector{Statement}
    ivs::Vector{Variable}
    dvs::Vector{Variable}
    vs::Vector{Variable}
    ps::Vector{Variable}
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)

# 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

# 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,Variable[],Variable[],Variable[],Variable[])

for the enum version, or

abstract type AbstractVariable <: Real end
abstract type AbstractStatement 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

const Expression = Union{AbstractVariable,AbstractStatement}

struct Differential
    x::Variable
end
function Derivative end
Base.:*(D::Differential,x::Expression) = Operation(Derivative,Expression[x,D.x])

struct Operation <: AbstractStatement
    op::Function
    args::Vector{Expression}
end
struct Statement <: AbstractStatement
    lhs::Expression
    rhs::Expression
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)

# 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

# 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[])

Basically, the big difference is that if the variable types are enum'd then they cannot be extended but all are of type Variable, while if they are AbstractVariable they can be extended. The real question is whether it ever makes sense to have a vector of these Variables where it wouldn't make sense to have an Expression. I can't find one, in which case typing it all to Variable doesn't improve concreteness anywhere in which case more extensible seems preferred.

Edit: I realized that constants need type information for units, so variables will always be type-unstable anyways, so we might as well make it fully type-based.

[ChrisRackauckas]

Started the repo at https://github.com/JuliaDiffEq/SciCompDSL.jl

[tkf]

This looks like awesome package and I can't help its release! I have some comments on a few points:

Delayed Variable

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

Does this syntax mean the DSL would make y(t) == y in some sense? Isn't it weird? Since the derivatives are defined via the operator D, how about defining it via shift operator?

T = Shift(t)
y_τ  = T(y, -5)  # or T(-5) * y
y_τ2 = T(y, -f)

I guess this approach also helps if you want to shifts in several axes of a "SDE".

Possible syntax sugars (they are not so nice):

@shift t y(t-f)  # need to tell what the independent variable is
@shift y([t]-f)  # mark it by []?
@shift y(t_-f)   # by _?

Bike-shedding @nl

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

This is a bike-shedding but... It can also be used for DiscreteProblem so naming it to be something that can cover it too would be nicer. A proper name would be @map or @mapping but that name in programming context would be confused with the map-over-collection map. How about @recurrence_relation? That's too long? Maybe too mouthful but @endomorphism also fits here. Or @endo as in Haskell.

Typed axis?

# Array Variables
x = DependentVariable(:x,1:10) # Array of variables

Are you planning to make variable shape fixed at the creation time? Wouldn't it be nice if you let user specify ndims? For extensible systems, it is necessary if you don't want to re-create the computation graphs.

Or, even better approach would be to let user declare axes types and make operation along the axes "type safe". That's what Nervana graph is doing (and maybe other deep net frameworks).

a = ArrayAxis(:a)
b = ArrayAxis(:b)
x = DependentVariable(:x, b)
y = DependentVariable(:y, a)
M = Parameter(:M, a, b)
E = y ⋅ M * x
wrong = x ⋅ M * y  # should throw at some point (maybe now)

[ChrisRackauckas]

Does this syntax mean the DSL would make y(t) == y in some sense? Isn't it weird?

I don't think that's weird. That's usually how it's interpreted. Also, we don't have another use for the call, so that works out well.

I don't get what you mean by the typed axis part.

Are you planning to make variable shape fixed at the creation time? Wouldn't it be nice if you let user specify ndims?

Yes, it would be nice to allow this too. I think the more important case though to get right first is making it easy to get n variables.

[AleMorales]

I wonder if using == will cause problems with logical statements on rhs (e.g. piecewise functions) being incorrectly transformed into a Statement.

[ChrisRackauckas]

That's a very good point. We might not want to make use of that syntax. Any other suggestions, or just a simple Statement(lhs,rhs) with a nice macro?

[AleMorales]

I would rely on dispatch as much as possible and push macros to a later stage of processing (it would also add visual noise to the code).

I think arrows might be nice. For example → or ← would not conflict with anything (at least in Base), but (i) maybe this has some meaning in mathematics that confuses people and (ii) I don't know how people feel about Unicode (though Julia makes it so easy!)

[tkf]

Okay, maybe y(t) == y is not weird. It's just me wanting to distinguish expressions/variables, functions and values.

The typed axis idea is that valid operation along the axis of tensor parameters and dependent variables have to have the same "type" (or label or whatever the name is). It helps, say, when you have discretized spatial system and don't want to mix the operation with respect to the x axis and y axis. In principle the error can be detected at graph creation time. Also, when there are several vector/matrix/tensor parameters, it is easy to verify the consistency of those parameters since it's clear what dimensions on what arrays have to match. Example:

a = ArrayAxis(:a)
b = ArrayAxis(:b)
xa = DependentVariable(:xa, a)
xb = DependentVariable(:xb, b)

# Parameters are block matrices
Maa = Parameter(:Maa, a, a)  # Na by Na matrix
Mab = Parameter(:Mab, a, b)  # Na by Nb matrix
Mbb = Parameter(:Mbb, b, b)  # Nb by Nb matrix
τa = Parameter(:τa, a)     # Na-dimensional vector
τb = Parameter(:τb, b)     # Nb-dimensional vector

de = @diffeq begin
  D*xa ./ τa = -xa .+ tanh.(Maa * xa .+ Mab * xb)
  D*xb ./ τb = -xb .+ tanh.(Mbb * xb)
end

where Na = size(Mab, 1) and Nb = size(Mab, 2) when Mab is a normal array. It is, as a whole, Na + Nb-dimensional ODE. Imagine that you happen to have test code with Na == Nb (that's not a great test but I suppose it can happen when there are many axes). Then your code will run with Mab replaced by Maa and it takes some extra time to debug it.

[ChrisRackauckas]

I see, yeah that makes sense. Open an issue.

[tkf]

@ChrisRackauckas Was it reply to my comment? Anyway, I opened an issue. Github issue is free. https://github.com/JuliaDiffEq/SciCompDSL.jl/issues/9

[jlperla]

@sglyon @albop

As discussed, here are a few thoughts on how this EDSL might be useful for economics applications, but you guys are much more equiped to think about this than I would @ChrisRackauckas : to give you some perspective, a fairly general sort of equation to think about solving is the stochastic difference equation in http://www.econforge.org/Dolo.jl/latest/model_specification.html#Euler-equation-1 You can also see the Dolo DSL there.

• You need to be careful when solving stochastic difference equations to keep track of information sets in the expectations, so part of this is the definition of "shocks" which are unknown at time t and generally normalized to be independent of the current state. I think the convention could be that these are the same as your Noise variables... though I don't really love that term. Although I don't think the economics term Shock is very general.
• For many of the algorithms for solving these sorts of models, it is important to be able to tell which variables are controls. In Dolo.jl you can see that these are listed as controls and in others like dynare these are flagged as predetermined... maybe just based on the timing convention
  • What this leads to, though, is that it would be very nice to be able to add in compile-time metadata to most of the definitions. For example I could go DependentVariable(:x, Predetermined()) or whatever.
  • To implement that sort of stuff, I can imagine that this could be handled very abstractly with traits, which could then be used for dispatch. i.e. something along the lines of

    struct DependentVariable{T} where T <: AbstractDependentVariableType
    name
    end

    where there could be a default NoMetadata type in AbstractDependentVariableType and handy constructors....

  • For economics applications, it is very important to be able to define variables defined on a discrete space (rather than the continuous one you have emphasized above).
  • I think this applies to everything, so my inclanation would be that there should be a domain as well, where when you define any type of variable you can provide the domain.
  • The default in a convenience constructor could be that it is in R or R^N if given an array, etc.
  • It would be very nice to leave this as something abstract to allow the state to be one of an enumeration
  • For discrete time applications, it would be nice to be able to specify an open-ended domain for t = 0, 1, ... as opposed to having to pass in a finite number of states. This would also be useful for simulation of markov chains with infinite transition matrices (e.g. M/M/1 queues and counting processes)
• Speaking of discrete time, for solving stochastic difference equations it would be important to specify the recursive structure difference equation directly. I think this is a generalization of your delay, but where the delay is on the discrete domain. For example

  
  t = IndependentVariable(:t, :Integers)
  c = DependentVarible(:c) #A continuous variable, defaults to R
  gamma = Parameter(:gamma)
  beta = Parameter(:beta)
  univariatenormal = Noise(:Normal)
  shock = univariatenormal * sigma #etc.

eulerequation = c(t)^(-gamma) == beta c(t+1)^(-gamma) shock(t+1)

[ChrisRackauckas]

For many of the algorithms for solving these sorts of models, it is important to be able to tell which variables are controls. In Dolo.jl you can see that these are listed as controls and in others like dynare these are flagged as predetermined... maybe just based on the timing convention

I'm thinking the variables need some kind of flags part that can hold a few flags. For example, the flow flag like in Modelica (@tshort would have to explain to me what that actually does though 👍 ).

Speaking of discrete time, for solving stochastic difference equations it would be important to specify the recursive structure difference equation directly. I think this is a generalization of your delay, but where the delay is on the discrete domain. For example

Yes, the current DelayedVariable should be able to hold that information just fine.

I think this applies to everything, so my inclanation would be that there should be a domain as well, where when you define any type of variable you can provide the domain.

Yup, it probably needs domains.

[jlperla]

Yes. I would just make sure that the flags can be dispatched on at compile-time (as necessary).

Given integer domains, does the Shift operator you are thinking about generalize a Lag operator, and would this be directly connected to the difference operator?

For example, you could use

t = IndependentVariable(:t, :Integers)
L = Shift(t)
D = Difference(t)
D == 1 - L #Thinking along these lines?

Lag operators are nice because an implementation of the DSL could do all sorts of lag operator tricks, such as polynomials of the lag operator (see https://ocw.mit.edu/courses/economics/14-384-time-series-analysis-fall-2013/lecture-notes/MIT14_384F13_lec1.pdf and https://lectures.quantecon.org/jl/arma.html ) @tpapp Thoughts on this sort of stuff? You know much more econometrics than I do.

[ChrisRackauckas]

Yes. I would just make sure that the flags can be dispatched on at compile-time (as necessary).

No, it's not a good idea to overuse compile time here.

Lag operators are nice because an implementation of the DSL could do all sorts of lag operator tricks, such as polynomials of the lag operator (see https://ocw.mit.edu/courses/economics/14-384-time-series-analysis-fall-2013/lecture-notes/MIT14_384F13_lec1.pdf and https://lectures.quantecon.org/jl/arma.html ) @tpapp Thoughts on this sort of stuff? You know much more econometrics than I do.

Yes, this would be really good because then series approaches like those from Hamilton's Time Series Analysis book can be employed to do model reductions.

[tpapp]

In principle it could be a nifty feature for notation, but you have to handle the algebra on lag (and backshift, etc) operators to make it useful in practice. Also, while the notation can be shared, AFAIK solution method methods for either time series econometrics or model solution (eg Euler equation) problems are very different from DE methods.

I am still digesting the implementation, but I am also concerned about the explosion of type space.

[ChrisRackauckas]

We can go back to having a single Variable type, and then have symbols for the subtypes. If we don't make those enumerates then it would be extendable. Then DependentVariable etc. would just be functions that create a Variable with the right subtype.

[ChrisRackauckas]

Also, I've been thinking that a Statement is really just an Operation with ==. This fixes @AleMorales 's comparison problem as well. That would make Expression = Union{Operation,Variable} which Julia can actually fully optimize... I'm going to go give that a try.

[ChrisRackauckas]

Yes, everything works just fine with just two types.

https://github.com/JuliaDiffEq/SciCompDSL.jl/pull/10

[ChrisRackauckas]

Alright it looks like the design here is pretty much settled. @mikeinnes pointed out where are limitations are. Basically if there's branching code then the tracing method will miss some branches unless we add a macro or use something like Cassette. However, for our use in describing scientific computing models, branching would have other concerns anyways since that would involve discontinuities which would require special handling. Thus I think we have arrived at something that is somewhat the simplest form that can handle the full types of models we want to describe. I'm closing this issue since the development will continue at https://github.com/JuliaDiffEq/SciCompDSL.jl . Thanks everyone for your input. For concrete things you'd like to see, like domains, please open an issue in that repo.

[ChrisRackauckas]

Also, while the notation can be shared, AFAIK solution method methods for either time series econometrics or model solution (eg Euler equation) problems are very different from DE methods.

BTW, that's not necessarily true.

[ChrisRackauckas]

For those who are curious, the proof of concept from the OP is complete. See the README.

https://github.com/JuliaDiffEq/SciCompDSL.jl#introduction

Now we just need the sugar and the extra features. If you know of any GSoC students looking for a project, direct them here. It's quite an easy project to hack away at. Thanks again.