We need ODEFunction(::DiffEqSystem) that generates the ODE from the system. It should check to see that there's a single derivative in each equation and then build the function.

It should be able to catch lines like

z == x*y + x
D*x == x*y + z
D*y == z/y

and automatically handle them.

Ref: https://discourse.julialang.org/t/package-compatibility-caps/15301

I am wondering about for modeling toolkit. Is there a mathematical formalism for “what is a model?” I’ve worked with different groups of scientists and it seems like every field has an informal definition for their field.

I’m working on making scientific models more amenable to analysis and machine reasoning about the model semantics. I’d like rigorously analyze the problems and algorithms we work with, and I think formalizing the objects under investigation.

This is related to why I asked about “what does context mean” (#83)

Is anyone working on anchoring modeling toolkit to some formalization of general modeling?

Many models are characterized by dynamic inputs (a.k.a external forcings). Using functors is a nice way of incorporating such inputs into a model such that, inside the, one uses forc(t) and the logic for generating the forcing is described elsewhere (this is often not a property of the model but of the scenario being simulated). Alternatively, one may interpolate from a timetable of values (e.g. with Interpolations.jl) and then do forc[t].

At the moment it is not possible to incorporate such syntax when build the Operation as Julia would try to call the Variable or apply the getindex method. I am thinking this trick may work:

t = IndependentVariable(:t)
forc = Variable(:forc, :Forcing, t)

And add the (t) or [t] during code generation in a later stage. And, of course, have a Forcing constructor wrapper.

But maybe I am missing something or there is a better design for this?

We have ~ for operation equality. What syntax for operation inequalities?

We can define these automatically like is done with the @ode_def macro. This would require the complete derivatives system.

We should be able to automatically expand out derivatives. It should have an easy way to users to plug in and define the derivative of a function, and it should be able to simplify values to 1.

order_lowering fails after some recent changes. I think the culprit is this hack:

https://github.com/JuliaDiffEq/SciCompDSL.jl/blob/master/src/systems/diffeqs/first_order_transform.jl#L10

This method loses all of the context of the variable which then makes it run into trouble when anything like derivatives are applied to it. It should instead use the constructor from 0543563 to create the new variable with matching context. This would make the new variable have the right dependence and thus not zero when the differential is applied. However, it looks like there will be some structural change necessary to get the right Variable passed in there.

Additionally, I think this may need to be updated:

https://github.com/JuliaDiffEq/SciCompDSL.jl/blob/master/src/systems/diffeqs/first_order_transform.jl#L41

to use ~

Arithmetic involving literals (e.g. 2*v1 where v1 is a Variable) does not work. I cannot imagine a convenient DSL without them. But we cannot extend the registered functions to deal with Real, so I am not sure at the moment how to handle this (using Constant is possible, but all the nice syntax goes away...).

When try to run

julia> @IVar t

I get the following error message

ERROR: MethodError: Cannot `convert` an object of type Void to an object of type ModelingToolkit.AbstractDomain
This may have arisen from a call to the constructor ModelingToolkit.AbstractDomain(...),
since type constructors fall back to convert methods.
Stacktrace:
 [1] ModelingToolkit.Variable(::Symbol, ::Void, ::Type{T} where T, ::Symbol, ::Void, ::Array{ModelingToolkit.Variable,1}, ::String, ::Bool, ::Void, ::Void, ::Void) at /Users/mason/.julia/v0.6/ModelingToolkit/src/variables.jl:3
 [2] (::Core.#kw#Type)(::Array{Any,1}, ::Type{ModelingToolkit.Variable}, ::Symbol, ::Void, ::Type{T} where T) at ./<missing>:0
 [3] #IndependentVariable#5(::Array{Any,1}, ::Function, ::Symbol) at /Users/mason/.julia/v0.6/ModelingToolkit/src/variables.jl:37
 [4] (::ModelingToolkit.#kw##IndependentVariable)(::Array{Any,1}, ::ModelingToolkit.#IndependentVariable, ::Symbol) at ./<missing>:0

By default the DSL generates noise term likes:

g[n,m] = stoch_nm * sqrt(rate_m)

It is possible the declare a special noise scaling parameter eta so that the term rather becomes

g[n,m] = eta * stoch_nm * sqrt(rate_m)

An improvement as mentioned in the previous PR would be to have affine noise a la

g[n,m] = eta1 * stoch_nm * sqrt(rate_m) + eta2

(I hope that I understood this right)

the new call would be

rn = @reaction_network rType (eta1,eta2) begin
    ...
end

alternatively

rn = @reaction_network rType begin
    ...
end

if no noise customization is desired.

On the same topic. Now noise is automatically stored as a parameter (default the last parameter in the parameter array). This is advantageous since you can the change the noise rate without redoing the DSL call. However in some cases you want a fixed value. I think it would be advantageous to allow for number as well, like:

rn = @reaction_network rType (2.0,eta2) begin
    ...
end

now the multiplicative scaling would be 2.0, with no changing that. This can be especially usefull if you are only interested in using parts of the affine noise, like:

rn = @reaction_network rType (0.0,eta2) begin
    ...
end

SciCompDSL.jl isn't really a DSL which is causing confusion. It's pretty much a compiler IR for models which is a nice target for DSLs. So we should probably rename, but to what?

• SciCompIR.jl
• ScientificIR.jl
• ModelingIR.jl

etc? No name is really vibing too much for me.

I want the names that a system assumes to be decoupled from the actual usage of the system itself. The reason is because the importance of this package is really to be an intermediate of model transformations that DSLs can plug into. We are already starting to build many different DSLs which utilize this intermediate, and one thing I've noticed is that we don't want to be constrained in subtype naming due to the systems that we want to target. For example, we want to have some biological model (like @AleMorales 's) and name things Reactants and stuff like that, but when we build the ODE we want to just say "these variables are the dependent variables, these are the independent variables, etc." and have it run to generate ODEs from there, as opposed to requiring DependentVariable, IndependentVariable, etc. (these won't go away, but will just be the defaults).

Most of the change was pretty simple. 2132199 did most. However, there are two problems to fix up:

1. We need variables to know their dependence themselves in order to simplify derivatives to zero ( #47 )
2. We need to have extract_derivatives handle aliasing. For example, right here is a hack that allows both :Variable and :DependentVariable to be used in the nonlinear system ( https://github.com/JuliaDiffEq/SciCompDSL.jl/blob/master/src/systems/nonlinear/nonlinear_system.jl#L8-L13 ), but we should make this more general and allow extract_elements to take in a vector of symbols and extract both subtypes into that form. extract_elements should also have a place of putting "all other names", i.e. not giving them special treatment. For example, in the DiffEqSystem we want :Variable to be used for just some intermediate calculations, but in reality it's just anything that's not an independent or dependent variable that should be treated like that.

If we have this, then I think that we can build DSLs with its own naming systems and all of that, but still pull from the same transformation, Jacobian building, and compilation system, which is the ultimate goal.

@YingboMa can you look at (1)? @AleMorales can you take a look at (2)?

Index reduction algorithms for DAEs are described here:

https://www.modelica.org/events/modelica2017/proceedings/html/submissions/ecp17132565_OtterElmqvist.pdf

This package looks so cool. I was trying to use it based on the readme and came across this error.

It is probably user error, but why does this segfault?

# This causes a seg fault.
using ModelingToolkit

# Define some variables
@IVar t
@DVar x(t) y(t) z(t)
@Param σ ρ β

f(x,y,z) = z*x+y
Derivative(f, x, 1)

   _       _ _(_)_     |  Documentation: https://docs.julialang.org
  (_)     | (_) (_)    |
   _ _   _| |_  __ _   |  Type "?" for help, "]?" for Pkg help.
  | | | | | | |/ _` |  |
  | | |_| | | | (_| |  |  Version 1.0.2 (2018-11-08)
 _/ |\__'_|_|_|\__'_|  |  Official https://julialang.org/ release
|__/                   |

julia> include("test/derivatives_fail.jl")
Segmentation fault: 11

Is the distinction between dependent and independent variables really necessary? Doesn't that arise at the moment of specifying the DE? I will want to use this stuff in other contexts (not DEs) where that distinction doesn't make sense.

State and time-triggered events are common in dynamical models, so modelling frameworks that generate SciComp IR may want to pass this information along, rather than hook the callbacks after code generation. After all, events are part of a model definition even if they are stored separately. I don't how abstract the event support should be, maybe for now we could allow storing a CallbackSet?

This right here:

https://github.com/JuliaDiffEq/SciCompDSL.jl/blob/master/src/operators.jl#L22-L27

I was really close to finishing the multi-arg function derivatives. Someone just needs to make the match the full if statement form below (i.e. it needs to make sure it's an operation before it indexes it like one and stuff like that).

#62 (comment)

Problem right now is that Operation needs a Vector{Expression} so without JuliaLang/julia#269 we cannot have Expression have Operation in it.

Spectral methods:

https://uu.diva-portal.org/smash/get/diva2:117236/FULLTEXT01.pdf
http://citeseerx.ist.psu.edu/viewdoc/summary?doi=10.1.1.142.2891
http://www.it.uu.se/research/publications/lic/2006-007/paperB.pdf
https://vtechworks.lib.vt.edu/handle/10919/30018
http://epubs.siam.org/doi/abs/10.1137/070689759
http://ieeexplore.ieee.org/document/6425804/

Currently, variables are equal exactly when they have the same name/subtype. However, this can be problematic in various cases, with variables created in two different contexts being treated as equal (e.g. first-order transformation).

The most obvious solution would be to treat each variable instance as strictly unique, attaching an ID at creation. At first glance, this doesn't seem conducive to testing when functions create new variables, but it's possible that some workarounds (either internally or in the test suite) could trivially solve this.

I suppose this is planned all along. ;-)

A useful extension of SciCompDSL would be to specify functions of parameters to allow coding of within and across equation restrictions. For instance, the first order ODE with speed of adjustment a

Dy = abX + ac - ay written as Dy = F1X + F2 -ay with F1 = ab and F2 = ac ; a, b, c are parameters

https://link.springer.com/article/10.1007/s10928-018-9570-4

Right now you have a single subtype::Symbol to differentiate types of variables. Since it is a Symbol it seems very generic, but it also means you can only categorise by one type.

I would like to be able to put in multiple levels of categorization given information known at compile-time. Could either the subtype be extended to support a named tuple of symbols, or could there be an additional flags which allows a set of Symbol flags (also preferably as a named tuple).

When there are multiple array variables and parameters in the model you have, you need to match the sizes of them. Rather than manually verifying the consistency of them, a better approach would be to encode the semantics of each "axis" in the computation graph. It enables SciCompDSL.jl to catch some errors at graph creation time.

Example (from my previous comment: SciML/DifferentialEquations.jl#261 (comment)):

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.

Possible benefits:

• 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.
• When the numeric arrays are passed to the problem constructor as the parameter values, constructor can check the consistency of them before running the graph. In the above example, e.g., size(Maa, 1) == size(Mab, 1) has to be satisfied.

Related project:

• Intel® Nervana™ graph (a deep learning library) https://ngraph.nervanasys.com/docs/latest/axes.html

The precedence level of == is only 6 (according to Base.operator_precedence(:(==)), whereas the arrows (like → or ⟺) have a precedence level of 5. This is an issue with the ReactionDSL (AleMorales/ReactionDSL.jl#4) but I think it could appear with other DSL that use SciComp as target for code lowering.

Since the == is actually being used as a form of assignment (to variables, derivatives, etc), I think it would make more sense to use one of the operators from the assignment-like operators prec-assignment.

A potential candidate could be ~ which is not used as binary operator in Base. There is precedent for this symbol, as this was the choice in probabilistic programming languages like JAGS, Stan or Turing to assign distributions to stochastic variables.

It would be nice to build the transformation equations for these kinds of things https://docs.sciml.ai/latest/extras/sensitivity_math/ directly from the IR since that would make it much faster by avoiding differencing.

Trying to reproduce all of the @ode_def functionality. Looks like invW is failing:

using ModelingToolkit
using Base.Test

# Define some variables
@IVar t
@DVar x(t) y(t) z(t)
@Deriv D'~t # Default of first derivative, Derivative(t,1)
@Param σ ρ β
@Const c=0
@Var a

# Define a differential equation
eqs = [D*x ~ σ*(y-x),
       D*y ~ x*(ρ-z)-y,
       D*z ~ x*y - β*z]
de = DiffEqSystem(eqs,[t],[x,y,z],Variable[],[σ,ρ,β])

u = [4.2,2.0,3.0]
p = [10.0,28.0,8/3]
gam = 5.3
t = 0.0
invW = zeros(3,3)
jac = zeros(3,3)

test_jac = eval(ModelingToolkit.generate_ode_jacobian(de))
test_jac(jac,big.(u),big.(p),t)
jac

test_w = eval(ModelingToolkit.generate_ode_iW(de)[1])
test_w(invW,big.(u),big.(p),big(gam),t)
invW

using ParameterizedFunctions

g = @ode_def LorenzExample begin
  dx = σ*(y-x)
  dy = x*(ρ-z) - y
  dz = x*y - β*z
end σ ρ β
jac2 = zeros(3,3)
g(Val{:jac},jac2,u,p,t)
jac2
jac2-jac

invW2 = zeros(3,3)
g(Val{:invW},invW2,u,p,gam,t)
invW2
invW2 - invW

Right now order lowering doesn't fit into any of the design spots. It should be a transformation of a DiffEqSystem to a DiffEqSystem, but instead right now it's a function on eqs. This fact has made it fragile and had it break when other variable contexts changed. @YingboMa

A function registration macro should be added that adds dispatches on Expression to create an Operation.

This should be a high priority to fix:

julia> WARNING: Method definition max(ModelingToolkit.Expression, ModelingToolkit.Expression) in mod
ule ModelingToolkit at C:\Users\Chris\.julia\v0.6\MacroTools\src\utils.jl:260 overwritten at C:\User
s\Chris\.julia\v0.6\MacroTools\src\utils.jl:260.
WARNING: Method definition max(Number, ModelingToolkit.Expression) in module ModelingToolkit at C:\U
sers\Chris\.julia\v0.6\MacroTools\src\utils.jl:260 overwritten at C:\Users\Chris\.julia\v0.6\MacroTo
ols\src\utils.jl:260.
WARNING: Method definition max(ModelingToolkit.Expression, Number) in module ModelingToolkit at C:\U
sers\Chris\.julia\v0.6\MacroTools\src\utils.jl:260 overwritten at C:\Users\Chris\.julia\v0.6\MacroTo
ols\src\utils.jl:260.
WARNING: Method definition min(ModelingToolkit.Expression, ModelingToolkit.Expression) in module Mod
elingToolkit at C:\Users\Chris\.julia\v0.6\MacroTools\src\utils.jl:260 overwritten at C:\Users\Chris
\.julia\v0.6\MacroTools\src\utils.jl:260.
WARNING: Method definition min(Number, ModelingToolkit.Expression) in module ModelingToolkit at C:\U
sers\Chris\.julia\v0.6\MacroTools\src\utils.jl:260 overwritten at C:\Users\Chris\.julia\v0.6\MacroTo
ols\src\utils.jl:260.
WARNING: Method definition min(ModelingToolkit.Expression, Number) in module ModelingToolkit at C:\U
sers\Chris\.julia\v0.6\MacroTools\src\utils.jl:260 overwritten at C:\Users\Chris\.julia\v0.6\MacroTo
ols\src\utils.jl:260.

This DSL has enough information that we could add detection of higher order ODEs and then make it convert down to a system of first order ODEs for the ODE solvers.

We need a conversion function that sends operations to a Julia AST for building the functions.

Each variable should have a sizing that determines if it's an array or not, but also what kind of array it is.

For variables, etc. we need to be able to assign domains. I would love to see this as the same one used for the differential operators DSL (which I think was your plan)

While the interface could be left very general, the key cases to test for are:

1. Defined on the reals, maybe as a default Reals
2. Defined on a compact set of the reals, RealsInterval(-1.0, 1.0). or even not compact RealsInterval(0.0, inf)
3. Defined on integers or Integers and IntegersInterval(-3, 4) and IntegersInterval(0, inf)
4. Defined on a discrete subset of a set, X = [1 4 5 6] then FiniteDomain(X)
5. Finally, it would be nice to have domains on finite numbers of variables with named states. A few things come to mind:

• In an enumeration @enum Fruit apple=1 orange=2 kiwi=3 and EnumDomain(Fruit). The last domain would be especially useful when working with markov chains and it would be great to be able to write code in terms of the names of states rather than position in a state vector.
• I don't really understand enums well in julia, so it could be that it is better to use a finite domain of symbols, Fruits = (:apple1, :apple2, :apple3) and FiniteDomain(Fruits)

@YingboMa I'm not sure that new display mechanism is all that helpful since it grows in complexity very fast. IMO operations should use a more succinct display than the variables just because they should be able to quickly show what the computational graph comes out to be.

There must be a way to do simple vs complex expressions to display?

The readme says:

The core idea behind ModelingToolkit.jl is that mathematical equations require context, and thus any symbolic manipulations and full model specifications requires the ability to handle such context.

Does context have a specific definition for this package? Is it the variable names and the information about what those variables are for? Like Dvar IVar, Parameter.

Hi, @ChrisRackauckas , @YingboMa
I am working with a DB where ODE models (parameters, diffeqs, init values, etc) are stored in a language neutral form. So, I can get all the components of ODEProblem from DB and convert (parse) them in Julia to Expressions like:

ex = quote
    v1 = k1*A*comp1
    dA = -v1/comp1
    dB = v1/comp1
end
params = (:k1, :comp1)
depend_vars = (:A, :B)

The interface I am looking for is a way to turn such input to DiffEqSystem and ODEProblem. Is it already supported in macro-free version?

As per this issue, I am going to list the things needed for an optimal control problem. Really there are three things, states, controls and parameters.

State Variables

This is equal to the number of differential equations, typically x1, ....xn were n is the number of states

Control Variables

Typically u1, ....um were m is the number of controls

Parameters

These might be p1.....pq where q is the number of parameters.

There are also a few more things that I have listed here, but I am not sure if I need different API for things like the initial state etc.

@ChrisRackauckas is this what you were asking for?

We need more routines to simplify equations down to common forms. It doesn't need to cover all of the fancy cases, but we should be able to have -x be the same as -1 * x.

@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

would be very nice. Basically it just assumes you want the name to match the variable name, which is very reasonable and probably recommended. It can also allow for having the value placed after the names, so like:

@DV x 2.0 y 3.0 z 4.0

It would make the upconversion of arrays automatic, but then people would have to know to duck-type for derivative tracing.

As mentioned in the OP of SciML/DifferentialEquations.jl#261 , it would be nice to be able to support Components. An AbstractComponent is a miniature system, and then you can compose components to build complicated systems. The Modelica example is that you can build components for things like Resistor and Capacitor, and then piece these together to build a circuit really easily. We would like to have something like that.

I am thinking that the base form could be:

struct Lorenz <: AbstractComponent
  x::Variable
  y::Variable
  z::Variable
  eqs::Vector{Expression}
  flags::Vector{:Symbol}
end

lorenz_eqs = [D*x == σ*(y-x)
                       D*y == x*(ρ-z)-y
                       D*z == x*y - β*z]
Lorenz() = Lorenz(@Var(x),@Var(y),@Var(z),lorenz_eqs,[:none,:none,:flow])

which could then have a macro form:

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

I'm not so familiar with Modelica so I am hoping that someone like @tshort or @smldis could help bikeshed this to make sure this could be feature-complete with it. Here's the Modia paper for reference.

Then of course we'd define some functions on AbstractComponent, but I'd need to find out what that connect stuff is supposed to be doing...

In line 18 of the master-branch there is

Base.convert(::Type{Operation},x::Variable) = Operation(identity,x)

and at line 63 there is

Base.convert(::Type{Operation},x::Variable) = Operation(identity,Expression[x])

which causes a warning when the package is loaded.

This library hit an interesting point. We have all of the features to recreate things like @ode_def, including computation of things like the inverse of W = I - gamma*J where J is the Jacobian.

But at the same time I have to admit defeat. The reason is because the computations are getting slow. Now, inverting matrices in computational graphs is quite a stress test so that's okay, but it's taking 10's of seconds instead of the about a second it took with ParameterizedFunctions.jl and SymEngine.jl, so I am afraid the current architecture will not scale well at all.

The issue goes back to our discussion and hope for context. We want to know a ton of information about our variables, but then actually keeping that information inside of the computational graph is too much. We wanted it so that way properties could be well-defined, but then what has ended up happening is that the context is used to gather the variables into specific arrays in the system's construction, and then only used from those arrays. So the computational graphs for the equations holds information it doesn't use only to automate the finding of variables, which in many cases is overridden so that way the variables end up in a different order than what the automated finder gives (since it's not all that great...).

Now let me be clear. Graph construction is cheap, but graph operations are bad.

I think we have a clear path around this though. The true format for the equations should be Julia AST earlier. Definitely operations on the graph should be done on Julia ASTs (something @YingboMa was considering before), but the question is how far back do we go.

1. Do we have the user give us vectors of expressions and then vectors of the variables for each time, more manually building the system?
2. During system construction, do we turn the eqs into Julia AST?

If we do 1, the interface changes quite a bit. If we do 2, we can keep the same user interface and switch up the internals to make it more optimized just for this goal. It would also let us keep the automation. I am inclined to go with route 2.

Either way, we would need to re-write our graph computations. But since they are all in Julia ASTs we can better make use of some existing tooling. I am considering using Reduce.jl here, though the only issue is Windows packaging which I've discussed with @chakravala. That would give us simplification, differentiation, etc. routines for free and it would be better than what we had. Thus, from route 2, we can keep the same user-facing tooling but redo the internal computation part.

With a Julia AST, we will need to enforce that the symbols used for the names are unique. That's one thing to note. We can make the AST use componentized names, i.e. comp1.a would be the name of the variable in component 1 and that would solve #36 . It is actually a valid Julia symbol: x = Symbol("comp.a"). Building the variable vectors would then require doing this kind of renaming.

@AleMorales might have some useful input here. I've also mentioned this to @jrevels but I don't know if Cassette really does all that much here (?).

I was really excited to find this package, I have been working on something very similar in Python for a long time but ended up switching to Julia recently for performance reasons. Getting started has been somewhat slow due to the lack of documentation, though. I can try to contribute by creating a skeleton for the docs if that would be helpful.