Components

[ChrisRackauckas]

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...

[AleMorales]

Not very familiar with Modelica/Modia, but would a collection of Components generate in the end a system of equations (like a DiffEqSystem)? If so I think it might be better to have a separate project (ComponentDSL or something like this) than then generates SciComp IR. That way here we could focus on the IR specification and its manipulation rather than adding layers of abstraction. But I think it is still a good idea to start the issue here.

[ChrisRackauckas]

In the end it generates equations, but the IR should be able to hold things in its component form because the component expansion is a complicated process. In fact, if you see the way things like Modelica are structured you'd notice that all of its modeling packages are essentially defining components and connections between them. Thus I think that this is a pretty fundamental building block that we should properly define here, and then leave open the definition of actual components (like circuits) to the modeling extension packages.

[tshort]

All of this activity on SciCompDSL is awesome!

For tying components together, the connect add's equations to the system for each node connected. For "potential" variables like voltage in an electrical system or temperature in a thermal system, the equations assign the potential at each node to be the same. For "flow" variables like current in an electrical system or heat flow rate in a thermal system, the equations sum the flow variables at the node to be zero.

@ChrisRackauckas, for your example formulation, why did you separate out flags from the variables themselves? It seems like variables will have several properties that might be better stored as part of the variable. That could include a flag indicating a flow, but it could also include unit information, an initial value, a string description, and more.

I'm also a little uncertain as to how SciCompDSL will fit into the ecosystem relative to Modia. Speak up if I am interpreting this wrong:

• SciCompDSL mainly focuses on the IR. This can be consumed by many packages for many different purposes.
• At a high level, the API for generating the IR focuses on helping the user with bookkeeping related to model variables (along with the benefits of having the IR like being able to generate Jacobian functions).
• Modia mainly focuses on the high-level view of a dynamic system. Modia is more monolithic relative to SciCompDSL. SciCompDSL may eventually cover many of the areas Modia covers.

I think it would be a good exercise to see if the Modia front-end language could be the front-end language for SciCompDSL. Another good exercise would be to compare the IR that SciCompDSL has with the internal representation used in Modia. Is there a chance the efforts can merge? I'm pinging @MartinOtter and @HildingElmqvist to see if they have any input and to ask if they could give @ChrisRackauckas access to Modia to facilitate further conversations on this.

[ChrisRackauckas]

@ChrisRackauckas, for your example formulation, why did you separate out flags from the variables themselves? It seems like variables will have several properties that might be better stored as part of the variable. That could include a flag indicating a flow, but it could also include unit information, an initial value, a string description, and more.

Yeah, there are flags and things like this in the Variable we are planning to use for other purposes, but I didn't realize that :flow is a flag generally on variables themselves instead of a Component. Since you say it's part of the Variable, then we could just handle it there. String description is a good idea too. Unit information would be in the value_type field. An initial value is what the value field is for.

Related: #21

I'm also a little uncertain as to how SciCompDSL will fit into the ecosystem relative to Modia.

I think it would be a good exercise to see if the Modia front-end language could be the front-end language for SciCompDSL. Another good exercise would be to compare the IR that SciCompDSL has with the internal representation used in Modia. Is there a chance the efforts can merge?

It's a very good question. To understand the motivation for this repository you have to look at ParameterizedFunctions.jl which is a very simple DSL for DiffEq. The issue that started cropping up was that the data format being used was Julia expressions, which are great because of homeoiconicity but have the limitation that variables do not contain context. In DiffEq's models we have a growing amount of contextualization which is required in order to support stochastic differential equations, jump diffusion models, differential operators for partial differential equations, differencing operators for function maps like in economic time-series models, and we needed things like domains for the PDE and optimal control models.

Altogether, when starting to implement abstract composition of lazy PDE operators (SciML/DifferentialEquations.jl#260) I realized that I cannot avoid having contextualized IR for these models. In addition, there were requests to extend ParameterizedFunctions.jl, but I realized that this would be difficult because its internal data structure is essentially the native Julia expressions which are not easy to query for information ("does this statement have any independent variables? If so, how many are there?"). So I gathered the DiffEqDSL.jl ideas that brewed over the last year and prototyped it in (SciML/DifferentialEquations.jl#261).

It became quite apparent early on that the functionality of this IR was almost a superset of what Modia.jl was aiming for. But as an IR it's focus is just the internal data structure for programmatic extension and manipulation, so this project is more along the lines of building a nice scientific computing compiler instead of a syntactic frontends (like TensorFlow). Since it's an IR, our plan is to stick different macro-based frontends onto it and write different system conversion methods for it, so that way the IR itself is disconnected from its input styles and its compilation targets (so even if it does support Modia/Modelica, it doesn't have to internally be like Modia/Modelica and instead just have at least a superset of representations). For example our reaction network library https://github.com/JuliaDiffEq/DiffEqBiological.jl has plans to be a generator for such IR SciML/Catalyst.jl#27. Thus one of the main reasons for adding components was to make the IR a full superset of Modia's IR so we could allow conversion to and in some cases from SciCompDSL (but this functionality is interesting in its own right of course). As for the compilation targets, the README shows pretty generic specifications of systems which currently have functions for targetting DiffEq's ODEProblem and NLsolve.jl with plans for doing things like detection of DAEs and automatic index reduction.

And then the "monolith" idea comes into play. There's a wide range of systems, variable subtypes, etc. defined in the repository, but in actuality there's no bounds to defining your own which is how DiffEqBiological is going to define some things. While we plan to build the basics into the DSL, transformation functions on the system IR is part of a standardized extension interface, so things like conversion to first-order ODEs and Lamperti transforms of stochastic differential equations are perfectly fine to be defined in extension libraries. Our computational graph's primitives also are constructed via a formalized interface of function registration and derivative dispatch registration, so new functions and derivatives can be easily added by users. From that it should be pretty clear that the core goal is to build an extendable system with enough basic constructs that most models can do things on their own with little overlap.

So where does Modia fit into this? Here's what I am envisioning. To me, Modia would be great as a Modelica front-end over the IR. We should try to get feature-complete with Modia so that its macro-based form of Modelica can easily generate IR code. I am planning on getting implementations of the index reduction and singularity removal mostly because I would want to find extensions to things like SDAEs and and neutral delay differential equations (a large part of creating an extendable system is that it makes it easier to do such research and share the results in a usable way!). The IR is not difficult to work on, so we have a GSoC potential @YingboMa who is looking to take on parts of this. In addition, if Modia wanted to write to Modelica code as output from the IR it could theoretically do that.

Of course, I do not have the years of experience writing a modeling language that @MartinOtter and @HildingElmqvist have. I would love to have them chime in and discuss how they think the two could interact and what changes they see necessary. Modelica has the advantage of being well-established with a large ecosystem to plug into, but we have an advantage of being highly extendable with a pretty big open-source workforce chugging away. There's benefits to be had for both of us.

[smldis]

I am afraid I cannot help a lot since i never used modelica, but reading the papers it seems there is great attention to be paid to symbolic manipulation and related physics laws (flow conservation, Ampere law, potential energies etc..). Starting by getting this core part correct and extendable would provide a great starting point to DSLs. This is why I think this package name could be made clearer (or rather factor out the symbolic part of it). With respect to MathOptInterface and Jump ecosystem I see another layer in between (with Pantelides algorithm and such things part of the middle layer).

[ChrisRackauckas]

With respect to MathOptInterface and Jump ecosystem I see another layer in between (with Pantelides algorithm and such things part of the middle layer).

Maybe we should rename, because the point of this is to be that "MathOptInterface layer" but allowing all sorts of symbolic manipulations, and leaving alternative simplified syntax (Jump) and domain-specific tooling as add-on components.

[ChrisRackauckas]

@ChrisRackauckas, for your example formulation, why did you separate out flags from the variables themselves? It seems like variables will have several properties that might be better stored as part of the variable. That could include a flag indicating a flow, but it could also include unit information, an initial value, a string description, and more.

Following @tshort 's suggestion and Modia's example, these are all part of the variable. My component interface is now shaping up to be something like:

struct Lorenz <: AbstractComponent
  x::Variable
  y::Variable
  z::Variable
  σ::Variable
  ρ::Variable
  β::Variable
  eqs::Vector{Expression}
end
function generate_lorenz_eqs(t)
  D = Differential(t)
  [D*x ~ σ*(y-x)
   D*y ~ x*(ρ-z)-y
   D*z ~ x*y - β*z]
end
Lorenz(t) = Lorenz(@DVar(x(t)),@DVar(y(t)),@DVar(z(t)),@Param(σ),@Param(ρ),@Param(β),generate_lorenz_eqs(t))

Then we could put a macro form to it:

@component struct Lorenz(t)
  # Define variables as you normally would
  @DVar x(t) y(t) z(t)
  @Param σ ρ β
end begin # Define the equations
  D = Differential(t)
  D*x ~ σ*(y-x)
  D*y ~ x*(ρ-z)-y
  D*z ~ x*y - β*z
end

[ChrisRackauckas]

Notes from our call:

Components as structs. Build equations. Define relations by equalities Heart.left_ventrical_flux ~ Lung.outside_pressure/Lung.α. Constructor could give a name prefix: Kidney(:kid1). Change equations inside of another by using an equality to a parameter: Kidney1.α ~ α2(t) / (ω * Kideny2.y(t-τ)) as adding a regulation term. Somehow add terms to an existing ODE: D(Heart.z) ~ Heart.eqs.D(z) - a*kidney2.state.

[asprionj]

In a broader sense, this is about model hierarchy, right? I mean, one wants to be able to create and/or re-use components, but also to build assemblies / submodels out of components. This way, the modelling process (for complex systems) and the collaboration on such models gets way more convenient and efficient. Instead of passing 20 components and 50 connections to someone, one would just hand over a single assembly with, say, the 3 relevant "ports" (or "signals" or "in-/outputs"). Then, connecting these 3 signals to another assembly is much simpler than finding and connecting the relevant ones in a list of 2x20 components...

The assemblies can be purely "virtual", i.e. the model can be flattened before anything more is performed. It's just for the humans working with the models.

I don't like graphical / block-diagram tools for modelling components ("visual programming" is not efficient), but for complex ones it is nice since one can just visually define/check connections and structure the model using arbitrarily many levels of hierarchy.

[ChrisRackauckas]

Yes, that's exactly what's going on here. And yes, we won't be building graphical tools. That's too much time and isn't mathematically interesting, but we will leave open the opportunity for others to build (and commercialize) GUIs based off of this tool.

[asprionj]

Sure, graphical tools are more of commercial nature; making stuff accessible / efficiently usable by less nerdy guys 😉. Should we start a meta-discussion on discourse about all this? I mean, there are many aspects and several (open-source) tools around. I think aspects like

• equation-based / symbolic modelling
• "imperative" modelling (causal modelling, "the needs of control-system engineers")
• components and hierarchical model building
• parametrised models (functions)
• calculating derivatives efficiently, for sensitivity analysis of all kinds (implicit integrators, parameter identification, linearisation, uncertainty propagation, static/dynamic optimisation, ...)

and tools / languages like

• Modelica (OpenModelica, JModelica)
• Modia
• CasADi -- they're considering a native Julia interface
• Scicos (don't know if this really is related)
• JuMP
• I'm sure there's a lot of stuff in the Python-verse
• ...

I think we should have a clear vision we're targeting at, based on what's already around (and may be (re-)used, or avoided, or taken inspiration from). Then focus all efforts, as far as reasonable, towards this vision and not do anything twice in a too similar way...

[ChrisRackauckas]

We're targeting something that you can think of as a next generation Modelica, replacing the mix of the language an internals with something that is just the intermediate representation, allowing others to define languages on top of it. Components essentially allow for acausal modeling, and of course causal modeling falls out if you don't add in any constraint equations and instead make it all pushing forward. It's different from Modia since Modia is a Julia-based Modelica: I could see Modia down the line actually use ModelingToolkit.jl for some of its internals down the line. JuMP's DSL could theoretically use MT down the line as well, and I've already talked with Miles about this as a way for handling some of the transforms people do in optimization. CasADi isn't about problem transformations, index reductions, distributed delay equations, etc., so it's not really to related.

Basically, no one let's you build a gigantic stochastic differential-algebraic equation and then automatically perform a Lamperti transform to get an additive noise version to use with DiffEq's order 1.5 additive noise methods. No one is even close, and no one really has a hackable system where I could easily add that on. Things like Modelica, CasADi, and Modia bake their transforms into their language system instead of just being a compiler IR for writing transforms, so it's a very different goal.

[ChrisRackauckas]

Should we start a meta-discussion on discourse about all this?

If it's for the purpose of letting contributors know what can be done? Sure! If it's a discussion telling us what we should do? No, that's not helpful. I know what my projects need from this, and will be building it towards what we need (and according to the funding that we are likely to be pulling in). If someone wants creative design input in the project, they can contribute or give $$$s, but at this point we don't need unfunded mandates (open source has far too many of those). There's already 1 billion things we know we need to do. Not trying to be combative about that, but meta-discussions from people who aren't putting in work are rarely helpful since most of it ends up just explaining how what they want is already going to be solved by the design that you would be working on if you weren't discussing it...

[asprionj]

We're targeting something that you can think of as a next generation Modelica [...]

I really like this approach, especially the "separation of concerns". And this also directly implies the direction in which my second thought was pointing:

If it's for the purpose of letting contributors know what can be done?

Yes. Or, from another point of view: to let people looking for "such stuff" know what's there and where it is going. Then, they can decide whether they can just start using it (and report bugs, suggest improvements, etc.), contribute to it to get out of it what they need, or start their own thing since there's no substantial overlap.

I'm personally in a situation where I developed, at work, a "(very) simple causal modelling framework without the GUI-overhead", relying heavily on DiffEq.jl. It's nice to use (for people who like code), flexible, fast and does what I need. I already started to conceptualise a more polished / systematically thought-out version in my spare time, considering all lessons learned, which I intended to publish as a completely open-source project. (I could even convince my company to go open source with this so I wouldn't have to separate the two versions anymore.) But then I noticed that I was constantly asking myself questions that probably came up during the development of many similar projects and thought to myself "well, why do it again when there's possibly battle-proven concepts and implementations for most aspects of it? Why create another framework / tool that's too similar to probably 27 already existing ones?" Contributing to an already existing base would be much more beneficial for all sides. But to be able to do that, one has to "know what's there and where it is going" -- see above.

meta-discussions from people who aren't putting in work are rarely helpful since most of it ends up just explaining how what they want is already going to be solved by the design that you would be working on if you weren't discussing it...

Agreed; constantly explaining things over and over again get's boring, tiring and eats up your time. The thing is, most of those people aren't lazy or wanting to prey on your work. Most are probably excited/interested but just don't know what's there and where it's going (that was definitely the last time I repeated this phrase 😉). If there's docs or some previous discussions, point them to it. And yes, sometimes it's good practice not to showcase anything too early to avoid too many of those discussions -- and from that point of view, starting an open meta-discussion is a very bad idea.

[asprionj]

From the 3rd post:

In the end it [the component] generates equations, but the IR should be able to hold things in its component form because the component expansion is a complicated process.

Could you explain this in some more detail? I can't follow, sorry. What about components has to be "expanded"? As far as I understand components, their connections need to be "resolved" in the parent scope, if these are not established directly by passing around expressions. For example, in some hypothetical modeling framework the user could define components as

# component 1: purely dynamic, no inputs
function component1()
    t   = Variable(:t, known=true)()
    x   = [Variable(:c1_x1)(t), Variable(:c1_x2)(t)]
    D   = Differential(t)
    eqs = [D(x[1]) ~ 1 - x[1],
           D(x[2]) ~ x[1] - x[2]]

    return x, [], eqs
end

# component 2: purely algebraic
function component2(u)
    y = [(u[1] + u[2])/2]

    return [], y, []
end

and pass the model structure (components and connections) to the framework. The framework would then instantiate and connect them (e.g. in a macro, or a model-building function) by

c1_x, c1_y, c1_eqs = component1()
c2_x, c2_y, c2_eqs = component2(c1_x)

and use c2_y as input to other component(s). (Finally, the framework could just collect all *_x and all *_eqs to build the system.) What could be gained over such a composition approach when having components in the IR? IMO this would unneccessarily complicate the "core engine" / manipulation of the IR. Especially once #136 is solved, because then multiple usage of an "output" *_y in other components would be nicely handled by the IR.

Another approach would be to allow "intermediate equations" in the list of equations, which could assign intermediate quantities to variables. These variables could then used to define "connections" by just equating two symbols (as stated in the "Notes from our call" post).

In fact, if you see the way things like Modelica are structured you'd notice that all of its modeling packages are essentially defining components and connections between them. Thus I think that this is a pretty fundamental building block that we should properly define here, and then leave open the definition of actual components (like circuits) to the modeling extension packages.

I don't know how Modelica (and Modia) handle components internally, i.e. in the compiler. Are they kept throughout the entire or at least some steps of the equation-manipulation process? Or are they expanded/resolved before touching any equations?

[jd-lara]

An alternative to handle connections is to have a topological struct that can be struct node <: Topological and others for multiterminal devices and a simple data template for any Device <: AbstractComponent that has to include a topological field.

This is pretty much what we do in PowerSystems.jl, every device that "injects" into the system has field Bus and every 2-port device has a NamedTuple (from = ::Bus, to = ::Bus). Then the device construction functions take care of adding the variables/equations accordingly.

[jamesjscully]

I am making a package for coupled systems at https://github.com/jamesjscully/CoupledODEToolkit (will move to .jl to name soon but git was being weird today)

This could probably be converted to work with modelling toolkit with relatively little work, but i found it easier to just define my own types.

It works like this:

@Component FitzughCell begin
    Dv = v - v^3 - w + Iext + gsyn, .0
    Dw = v + a -b*w, .1
    to_alpha = vpre -> ~v
    elec = gsyn -> gelec*(~v - v)
end begin
    Iext = .5
    a = :free
    b = .7
    gelec = .02
end

@Component AlphaSynapse begin
    DS = α*(1-S)/(1+exp(♉-vpre)/k) - β*S, 0
    synout = gsyn -> g*~S*(v-E)
end begin
    ♉ = -20
    k = 20
    E = -5
    g = .1
end

c1 = FitzughCell(:c1, elec = :c2)
c2 = FitzughCell(:c2, v0 = .1, w0 = .2, a = :free, elec = :c1)
s1 = AlphaSynapse(:s1)
net = Network([c1,c2])() # returns function

State variables are : D(sym) = (expr, ic)
couplings are : name = flag -> expr
The ~ Symbol is used to indicate that the variable belongs to the source, not target component.

[ChrisRackauckas]

This is what we are aiming for:

## ODE Example

# System 1
@parameters t σ ρ β f(t)
@variables x(t) y(t) z(t)
@derivatives D'~t
eqs = [D(x) ~ σ*(y-x),
       D(y) ~ x*(ρ-z)-y + f(t),
       D(z) ~ x*y - β*z]
lorenz = ODESystem(eqs)

# System 2
@parameters t α β γ δ f(t)
@variables x(t) y(t)
@derivatives D'~t
eqs = [D(x) ~ α*x - β*x*y + f(t)
       D(y) ~ -γ*y + δ*x*y]
lotka = ODESystem(eqs)

# Composed System

eqs = [lorenz.eqs;lotka.eqs;lorenz.f ~ lotka.x;lotka.f ~ lorenz.y]

#### Should be

[
D(lorenz_x) ~ lorenz_σ*(lorenz_y-lorenz_x),
D(lorenz_y) ~ lorenz_x*(lorenz_ρ-lorenz_z)-lorenz_y + lotka_y,
D(lorenz_z) ~ lorenz_x*lorenz_y - lorenz_β*lorenz_z,
D(lotka_x)  ~ lotka_α*lotka_x - lotka_β*lotka_x*lotka_y + lorenz_x,
D(lotka_y)  ~ -lotka_γ*lotka_y + lotka_δ*lotka_x*lotka_y
]

## DAE Example

# System 1
@parameters t σ ρ β f(t)
@variables x(t) y(t) z(t)
@derivatives D'~t
eqs = [D(x) ~ σ*(y-x),
       D(y) ~ x*(ρ-z)-y,
       D(z) ~ x*y - β*z]
lorenz = ODESystem(eqs)

# System 2
@parameters t α β γ δ f(t)
@variables x(t) y(t)
@derivatives D'~t
eqs = [D(x) ~ α*x - β*x*y
       D(y) ~ -γ*y + δ*x*y]
lotka = ODESystem(eqs)

# Composed System

eqs = [lorenz.eqs;lotka.eqs;lorenz.x ~ lotka.x]

#### Should be

[
D(lorenz_x) ~ lorenz_σ*(lorenz_y-lorenz_x),
D(lorenz_y) ~ lorenz_x*(lorenz_ρ-lorenz_z)-lorenz_y,
D(lorenz_z) ~ lorenz_x*lorenz_y - lorenz_β*lorenz_z,
D(lotka_x)  ~ lotka_α*lotka_x - lotka_β*lotka_x*lotka_y,
D(lotka_y)  ~ -lotka_γ*lotka_y + lotka_δ*lotka_x*lotka_y,
lorenz_x ~ lotka_x
]

[jamesjscully]

A few questions/comments:

1. what if you want to make a lot of components with common structure, and only a few differences in parameters. different types of synapses and neurons for example.
2. what if some components always couple at a flag with a certain expression that depends on the "type" of source, not the "type" of target? And if that expression also has a flag?
3. why f(t) in the eqs array? f is not callable.
4. are you suggesting to dispatch on vcat?
5. lorenz has no field x, etc

something like this would work, but its not very pretty in my opinion

function couple(systups)
    #replace the stuff appropriately
end


couple([
    (lorenz, [
        x => (lotka, f)]),
    (lotka, [
        x => (lorenz, f)
    ])])

I should be able to get this to work since i figured out how to convert Operations to Expr.

[ChrisRackauckas]

You may want to join the #symbolic-manipulations chatroom of the Slack where a good bit of this was being discussed. That was a direct comment to @HarrisonGrodin , there's a lot more context than what was shown there. Basically, our idea is to not flatten the expressions and instead keep it in a lazy form using types and a function-based (trait-based) interface.

why f(t) in the eqs array? f is not callable.

Because it depends on time.

what if some components always couple at a flag with a certain expression that depends on the "type" of source, not the "type" of target? And if that expression also has a flag?

I don't get what you mean here.

what if you want to make a lot of components with common structure, and only a few differences in parameters. different types of synapses and neurons for example.

You'd make multiple instances of that type.

are you suggesting to dispatch on vcat?

No, it's just markup syntax for @HarrisonGrodin . I am not sure the right way to write it right now.

[ChrisRackauckas]

Let me know if you want to chat about this @jamesjscully . You seem keen to get this working for something you need, so I'd like to help "onboard you" here, since there's a lot of stuff going on and so it might help to give you some of the context.

[jamesjscully]

@ChrisRackauckas Yes, I'm definately interested in contributing more and would love to chat. I don't necessarily need it for the stuff we do right now, because it's easy enough to just create a system manually, but I want to move towards working in a framework where other groups can easily reuse our code, and this seems like the right framework for that.

[ChrisRackauckas]

Feel free to ping me on the Slack in the #symbolic-manipulation channel.

[jonniedie]

Where did this leave off? Is any effort still being directed towards this? If not, I might want to pick it up. This type of modularity is the only thing that's really keeping me from being productive enough in the JuliaDiffEq ecosystem to replace a lot of my Simulink/Simscape workflow, so I'm pretty motivated to see this succeed.

[ChrisRackauckas]

This is where we are at:

#214

[ChrisRackauckas]

Completed by #265