
ModelingToolkit.jl, an IR and Compiler for Scientific Models - ChrisRackauckas
http://juliadiffeq.org/DiffEqTutorials.jl/html/ode_extras/ModelingToolkit.html
======
eigenspace
Awesome work, it’s really cool to see the progress being made in
ModellingToolkit.jl.

As someone who spends more time on symbolic mathematics than numerical
mathematics, I’ve often wished that Julia’s symbolic libraries had the same
sort of energy behind them that the numerical libraries had. But it’s great to
see that a lot of the work put into the numerical libraries is relevant to the
symbolic ones.

~~~
lewis500
Me too. I think that because of the meta programming capabilities we’re going
to see a lot more development in that direction. Recently I started using
Turing.jl, which uses a sort of declarative syntax that was easy to write and
understand.

------
marmaduke
as a research engineer in data science/modeling, posts like this just look too
good to be true. there's no free lunch. what is the catch with all this stuff?
who is going to maintain all this when Chris Rackauckas gets plucked by
MathWorks or Wolfram?

(this is a compliment! it's great stuff)

~~~
ChrisRackauckas
If you're looking for a catch, I guess the issue is that the modeling domain
is relatively simple, restricted to "mathematical models", i.e. no
conditional-based control flow like while loops. Like systems like Jax, for
loops are supported by simply unrolling the expression into something without
control flow. This lets the resulting IR be relatively simple and then all of
the math is easy to define, but of course means the domain for this is "math"
and not "programs". That said, it's still large enough for things like
Modelica models, which is a large enough domain for a lot of the algorithms I
want to explore like automatic index reduction of DAEs.

------
marmaduke
I have the obligatory does it work with GPUs question. I have lamented writing
out jacobians for GPUs by hand for ages, and thought, SymPy should be able to
do this.

~~~
new4thaccount
As someone whose entire industry runs off of a Jacobian matrix, I wonder how
easy it is to run on the GPU? I assumed this was difficult as splitting up the
matrix would impact the solution.

~~~
ChrisRackauckas
DifferentialEquations.jl supports GPUs. The question is what is the right way
to utilize DiffEq's GPU support here. That's actually a much deeper question
than you might at first think. The reason is because sparse factorizations are
not an inherently parallel operation, so CUSOLVER is actually quite bad at
doing them. It even seems that CUSOLVER does some CPU operations in there...
ouch.

But what can we do with ModelingToolkit to start pushing more into this
domain? Well, here's a few things:

1) We plan to have new compilation targets, kind of like the SArray one that
was shown, which target parallelism. For example, for large Jacobians, we want
a target so you can say "generate the code that fills the Jacobian optimally
with 16 threads", and then what it will do is split up the Jacobian into 16
parts with approximately equally many operations and thread between those.
Since it's all scalar operations then this should be something that gets as
efficient as you'd want to write by hand! Then a similar compilation target
can be made for the GPU, where it would need to know the number of GPU threads
you'd be using and it would split apart that Jacobian filling and write a
gigantic `if threadid() == ...` setup, and use CUDAnative.jl to compile this
to a .ptx kernel so that all of those threads are parallelized. Both of these
trivially extend to sparse matrices doing what I showed in the writing. So
filling sparse matrices is something we can and plan to have nice compilation
targets for.

2) We can also have compilation targets for `J*v` operations, so that way it
can accelerate Newton-Krylov solvers, which is probably the better way of
using sparse matrices on the GPU since it would avoid sparse factorizations.

3) We can also implement a generic sparse factorization in Julia, to then
build analytical solutions for sparse factorizations, run Rewrite.jl
simplification passes on these, and generate the parallelized kernels for
them. This would then be something that does sparse factorizations on the GPU
with parallelized code but does not need to use CUSOLVER, and so that would be
quite ideal. The issue is that it would require writing generic sparse
factorization tools, since Julia currently relies on SuiteSparse for this with
no generic fallback. This is not an easy task, but I think Yingbo will be up
for it in a few years since we have a plan to not use standard sparse linear
algebra routines in Julia's differential equation suites because there is a
way to speed up operations specifically in this domain by using a new
factorization, and so if we go this direction our first task will be to learn
how the current sparse LU is implemented.

~~~
new4thaccount
Wow! A lot going on there and my sincere thanks for your part in what be a
massive amount of free work that I get to benefit from.

I'm really excited for what y'all have coming in the future and will stay
tuned. There is some literature on Newton-Krylov solvers for my domain (power
systems) so I'll have to brush up on that math.

Help me understand one thing. I suppose if we could one day compile Julia to a
binary and sell software (just curious...I don't sell software today), would
you have to secure a license for suitesparse? Is that another benefit from
rewriting?

~~~
ChrisRackauckas
Yes, you would need a commercial license for suite sparse I think, so that
would be a good reason to rewrite its core methods to a new license.

~~~
new4thaccount
Thanks for clearing that up. I think someone else on here pointed that out as
well as the fact that Matlab and some other tools provide that for you as part
of your license/subscription fee. I started digging around and found out that
it seems like the suitesparse author is one of the top authorities in this
area and therefore, it might not be super realistic to get something as fast
or faster on your own. Maybe a big enough open-source project could do it
though. Any thoughts on this?

