10 ms·

Diffeq.jl v6.4: Full GPU ODEs, Neural ODEs with Batching on GPUs, and More

173 points by ViralBShah 7y ago

eigenspace 7y ago
Wow, that's really cool. I feel like DiffEq is one of those libraries that seems to leverage every one of Julia's cool features that I can think of. Making differential equations 'just work' on a GPU is super cool. I feel like if this were almost any other language, it'd be impossible to get this level of code-reuse where a user can just bring their own differential equation that may not even know about GPUs and then have the differential equation library solve it on the GPU by just giving a GPUarray as the initial condition. ________________________________ Also, since the blog post doesn't seem to link to it, here's a link to the actual DifferentialEquations.jl repo: https://github.com/JuliaDiffEq/DifferentialEquations.jl https://github.com/JuliaDiffEq/DifferentialEquations.jl
- ViralBShah 7y ago
  Yeah, it is amazing to see how one can really rethink a core mathematical library with new tools and ideas. My mind continues to be blown every time I hear Chris talk about Measurements.jl and DifferentialEquations.jl just work, and you get sensitivity analysis for free (and neither package had planned for the other).
- The_rationalist 7y ago
  well openmp or openacc are still simpler.
  - gugagore 7y ago
    Do you mean they are simpler to use? Because those are just models of parallel computation. A little like saying Turing machines are simpler than python.
adamnemecek 7y ago
I've been writing julia a bunch lately. It's a joy to use.
- amrrs 7y ago
  Were you a Python-convert? Could you share you transition journey like how much time or difficulty you faced.
  - speps 7y ago
    In my own recent experience, Julia is awesome but documentation could be sometimes clearer (eg. missing examples with each method, missing cross references). Also, when searching for details on some features, it's hard to figure which version someone is talking about, it seems some versions (0.6 for example) were there for ages and most forums refer to it but some stuff is now deprecated or replaced. Otherwise, I would probably use it for everything now given the performance and the very good FFI story.
    - ViralBShah 7y ago
      Please file issues when you find something is less than perfect. That's the first step towards making it better!
  - adamnemecek 7y ago
    There were some problems during the transition between 0.6 and 1.0 but that was kinda expected. After 1.0 it's been pretty breezy. Good IDE support (juno). Julia is fast. Numpy is actually c, julia is just julia. You have an amazing interop that just works. With python, matlab, c,c++,etc etc. https://github.com/JuliaPy/PyCall.jl https://github.com/JuliaPy/PyCall.jl Package manager just works. The language is a lot more powerful. Diffeq is around 1000 lines of code which is nuts.
    - ChrisRackauckas 7y ago
      DiffEq is a lot more than a 1000 lines of code. By now it's about 75 git repo and some of those repos are themselves >10,000 lines of code. And there is a single file in the test suites that is 10,000 lines defining mathematical constants for 200 Runge-Kutta tableaus. So it's way way more, and actually quite large by now. However, we still get a lot of feature bang for our lines of code buck.
sara7262 7y ago
Who is Satan and why there is evil Best video about this with proofs https://youtu.be/AWNM1-xYI5E https://youtu.be/AWNM1-xYI5E Funny monkey using phone in video Don't think just you can use phone See how he is searching and watching own monkey videos on phone http://bit.ly/2DS4p2R http://bit.ly/2DS4p2R See what happened when a lion was trying to kill man http://bit.ly/2H2zsKc http://bit.ly/2H2zsKc Best video with proofs Which showing there is creator of this universe https://youtu.be/j5HHl6-r3VU https://youtu.be/j5HHl6-r3VU Dogs are meeting their owners after a long time see in this video the reaction of dogs and try to control your cry http://bit.ly/2JmjuNw http://bit.ly/2JmjuNw Funny cats are teasing babies in video just try to control your laugh I am sure you can't http://bit.ly/2LJu58n http://bit.ly/2LJu58n
microcolonel 7y ago
Which GPUs? I'm guessing by the "CuArrays" identifier that this is CUDA-only, it'd be nice if they said NVIDIA GPUs instead of getting my hopes up. :'- (
- ChrisRackauckas 7y ago
  Well the implementation is generic to the array type, so if someone creates an array type which builds broadcast kernels on with OpenCL then it would work with OpenCL. Sadly, CLArrays.jl never got updated to Julia v1.0, so until someone does that it's just a dream. However, DifferentialEquations.jl will be waiting ready for that to happen :).
  - dnautics 7y ago
    I'm honestly surprised that amd doesn't get more aggressive about developing towards Julia as a deep learning target, it would be way easier for their devs to optimize for performance parity than say tensorflows spaghetti code. Iirc, the raw flops performance of AMD gpus is even or better than Nvidia but their unoptimized libraries kick it in the shins (I could be misremembering thus)
    - ViralBShah 7y ago
      There are folks working on the AMDgpu target for Julia, and I see that some early work has already happened in terms of building LLVM and such. I suspect the AMD toolchains are not as strong - but with the right contributors, it should be possible to leverage the work done on the Nvidia GPUs for the AMD GPUs and even other accelerators.
      - celrod 7y ago
        As someone who owns a couple of Vega (AMD) GPUs, I'm watching a few of jpsamaroo's GitHub repos. In particular: https://github.com/jpsamaroo/HSAArrays.jl https://github.com/jpsamaroo/HSAArrays.jl I don't know how far along it is (no documentation yet), but those are definite steps.
pepijndevos 7y ago
I'm contemplating writing some fun electromagnetic simulator, and so far had my eyes on Rust and/or Futhark with some prototypes in Python even Matlab, but now I may need to consider Julia too...
skdotdan 7y ago
Noob question, having played around with CUDA and NVCC a bit: in languages such as Julia or Python, why does every algorithm have to be specifically adapted to be able to run on GPUs? Couldn't algorithms be written from higher level building blocks? Eg. implement map, reduce and scanl in CUDA/OpenCL, and have a default parameter like backend='cpu' than can be set to 'gpu'.
- krastanov 7y ago
  Julia is doing exactly that (even better than setting a flag, it is based on subtypes). Algorithms are usually written to work with `AbstractArray`. The usual RAM/CPU arrays (i.e. the type `Array`) is a subtype. The CUDA array which is stored on the GPU is also a subtype. Most general purpose code does not care which type you use.
  - ViralBShah 7y ago
    Definitely look at the whole JuliaGPU ecosystem, and the CUDAnative.jl package. There are two big questions. First, how do you get code generation for targets such as GPUs from a dynamic high level language such as Julia, and that is what CUDAnative.jl achieves. The second question is what abstractions should Julia present to programmers so that increasingly larger codebases can automatically leverage GPUs. Packages such as CuArrays.jl are early answers, but there is much work to be done here.
- ChrisRackauckas 7y ago
  This is precisely correct. You just have to make sure in the thousands of lines of code that you only use higher level building blocks and make no assumptions greater than what a GPU-based array allows. GPU-based arrays don't even want to allow you to do scalar indexing for example, since every scalar indexing operation is a data transfer back to the CPU. So, with those rules in mind, getting the stiff ODE solvers to automatically work with GPU-based arrays amounted to: 1. Making sure every operation in the solvers used a higher level building block (@..) that automatically built fused broadcast kernels, and CuArrays.jl then supplies the broadcast overloads. That makes all of the non-stiff methods work. 2. Change the default choice of Jacobian to "dense array of size NxN" to "zero'd outer product of input array defined via broadcast". This translates to `J = false .* z .* z'` as a trick to make a Jacobian match whatever array type the array type library thinks is best, so J will be a GPUArray while DiffEq knows nothing about GPUs (also works for things like MPI, but GPUs is easier for people to reason about these days). So a one line change to the stiff ODE solvers made them all GPU compatible! But... 3. There's a small bug, or at least missing feature, in the GPU libraries that their QR doesn't implement backsolve. Technically the user can work around that because we allow someone to pass a linear solver routine with the `linsolve` argument, but it's easier to just handle it for everyone in the defaults. So using Requires.jl, if you have CuArrays installed, then we add a linsolve routine to the defaults for you, defined by https://github.com/JuliaDiffEq/DiffEqBase.jl/blob/master/src/init.jl#L104-L110 https://github.com/JuliaDiffEq/DiffEqBase.jl/blob/master/src... . Technically this can get upstreamed to CuArrays.jl and get deleted from DiffEqBase, but for now it'll live there so that everything is easier on users. So yes, make it high level, avoid assuming things like the Jacobian live on the CPU, and add one missing method to the GPU array library, and now when the user passes `u0 = CuArray(...)` for the initial condition, the algorithm compiles a version which does all state operations on the GPU (and all logic on the CPU, so it's a pretty good mix!).
Certhas 7y ago
Incredible work. Thank you everyone for the great release. This massively improves our usecase (sparse matrix equations) and I am looking forward to playing with it.
- ChrisRackauckas 7y ago
  Oh boy, just wait until next release! Sparse matrices are one of our big focuses for the summer: https://github.com/pkj-m/SparseDiffTools.jl/issues/2 https://github.com/pkj-m/SparseDiffTools.jl/issues/2
vsskanth 7y ago
Anyone tried using Julia's solvers to run FMI ModelExchange binaries ?
- ChrisRackauckas 7y ago
  Not yet, but we are looking into it. I'm hoping to get some funding for this, since it's a no-research software dev project though, which would be great for a just-out-of-undergrad student looking to transition to software dev, hence need money to pay such a student.
  - vsskanth 7y ago
    Thank you.
mlevental 7y ago
this happens literally every single time i try to hop on the julia bang wagon: run example -> example breaks. 1. fresh install of juno (1.1.0) 2. add latest packages 3. something goes wrong. this time it's p = destructure(model) from example page gives UndefVarError: destructure not defined top-level scope at none:0 https://i.imgur.com/8fCpVrR.png https://i.imgur.com/8fCpVrR.png how the hell does anyone get anything when the codebase is shifting this fast!?
- ChrisRackauckas 7y ago
  That's my bad. Somehow an earlier draft of the post was used. The post was updated with working code. I guess that's what happens when you lose internet and push to a Github repo from your phone :P.
  - mlevental 7y ago
    :) well thanks! edit: UndefVarError: JacVecOperator not defined :( also this https://gist.github.com/makslevental/fe9144d196cade685d333cd862c012c0 https://gist.github.com/makslevental/fe9144d196cade685d333cd... i realize it's just a warning but it shows another undef error
    - ChrisRackauckas 7y ago
      Are you on the latest versions? See what running a package update does. If that doesn't work, file an issue with your package Manifest so we can track this down. Thanks.

173 points by ViralBShah 7y ago

Wow, that's really cool. I feel like DiffEq is one of those libraries that seems to leverage every one of Julia's cool features that I can think of. Making differential equations 'just work' on a GPU is super cool. I feel like if this were almost any other language, it'd be impossible to get this level of code-reuse where a user can just bring their own differential equation that may not even know about GPUs and then have the differential equation library solve it on the GPU by just giving a GPUarray as the initial condition. ________________________________ Also, since the blog post doesn't seem to link to it, here's a link to the actual DifferentialEquations.jl repo: https://github.com/JuliaDiffEq/DifferentialEquations.jl https://github.com/JuliaDiffEq/DifferentialEquations.jl

ViralBShah 7y ago
Yeah, it is amazing to see how one can really rethink a core mathematical library with new tools and ideas. My mind continues to be blown every time I hear Chris talk about Measurements.jl and DifferentialEquations.jl just work, and you get sensitivity analysis for free (and neither package had planned for the other).
The_rationalist 7y ago
well openmp or openacc are still simpler.
- gugagore 7y ago
  Do you mean they are simpler to use? Because those are just models of parallel computation. A little like saying Turing machines are simpler than python.

adamnemecek 7y ago

I've been writing julia a bunch lately. It's a joy to use.

amrrs 7y ago
Were you a Python-convert? Could you share you transition journey like how much time or difficulty you faced.
- speps 7y ago
  In my own recent experience, Julia is awesome but documentation could be sometimes clearer (eg. missing examples with each method, missing cross references). Also, when searching for details on some features, it's hard to figure which version someone is talking about, it seems some versions (0.6 for example) were there for ages and most forums refer to it but some stuff is now deprecated or replaced. Otherwise, I would probably use it for everything now given the performance and the very good FFI story.
  - ViralBShah 7y ago
    Please file issues when you find something is less than perfect. That's the first step towards making it better!
- adamnemecek 7y ago
  There were some problems during the transition between 0.6 and 1.0 but that was kinda expected. After 1.0 it's been pretty breezy. Good IDE support (juno). Julia is fast. Numpy is actually c, julia is just julia. You have an amazing interop that just works. With python, matlab, c,c++,etc etc. https://github.com/JuliaPy/PyCall.jl https://github.com/JuliaPy/PyCall.jl Package manager just works. The language is a lot more powerful. Diffeq is around 1000 lines of code which is nuts.
  - ChrisRackauckas 7y ago
    DiffEq is a lot more than a 1000 lines of code. By now it's about 75 git repo and some of those repos are themselves >10,000 lines of code. And there is a single file in the test suites that is 10,000 lines defining mathematical constants for 200 Runge-Kutta tableaus. So it's way way more, and actually quite large by now. However, we still get a lot of feature bang for our lines of code buck.

sara7262 7y ago

Who is Satan and why there is evil Best video about this with proofs https://youtu.be/AWNM1-xYI5E https://youtu.be/AWNM1-xYI5E Funny monkey using phone in video Don't think just you can use phone See how he is searching and watching own monkey videos on phone http://bit.ly/2DS4p2R http://bit.ly/2DS4p2R See what happened when a lion was trying to kill man http://bit.ly/2H2zsKc http://bit.ly/2H2zsKc Best video with proofs Which showing there is creator of this universe https://youtu.be/j5HHl6-r3VU https://youtu.be/j5HHl6-r3VU Dogs are meeting their owners after a long time see in this video the reaction of dogs and try to control your cry http://bit.ly/2JmjuNw http://bit.ly/2JmjuNw Funny cats are teasing babies in video just try to control your laugh I am sure you can't http://bit.ly/2LJu58n http://bit.ly/2LJu58n

microcolonel 7y ago

Which GPUs? I'm guessing by the "CuArrays" identifier that this is CUDA-only, it'd be nice if they said NVIDIA GPUs instead of getting my hopes up. :'- (

ChrisRackauckas 7y ago
Well the implementation is generic to the array type, so if someone creates an array type which builds broadcast kernels on with OpenCL then it would work with OpenCL. Sadly, CLArrays.jl never got updated to Julia v1.0, so until someone does that it's just a dream. However, DifferentialEquations.jl will be waiting ready for that to happen :).
- dnautics 7y ago
  I'm honestly surprised that amd doesn't get more aggressive about developing towards Julia as a deep learning target, it would be way easier for their devs to optimize for performance parity than say tensorflows spaghetti code. Iirc, the raw flops performance of AMD gpus is even or better than Nvidia but their unoptimized libraries kick it in the shins (I could be misremembering thus)
  - ViralBShah 7y ago
    There are folks working on the AMDgpu target for Julia, and I see that some early work has already happened in terms of building LLVM and such. I suspect the AMD toolchains are not as strong - but with the right contributors, it should be possible to leverage the work done on the Nvidia GPUs for the AMD GPUs and even other accelerators.
    - celrod 7y ago
      As someone who owns a couple of Vega (AMD) GPUs, I'm watching a few of jpsamaroo's GitHub repos. In particular: https://github.com/jpsamaroo/HSAArrays.jl https://github.com/jpsamaroo/HSAArrays.jl I don't know how far along it is (no documentation yet), but those are definite steps.

pepijndevos 7y ago

I'm contemplating writing some fun electromagnetic simulator, and so far had my eyes on Rust and/or Futhark with some prototypes in Python even Matlab, but now I may need to consider Julia too...

skdotdan 7y ago

Noob question, having played around with CUDA and NVCC a bit: in languages such as Julia or Python, why does every algorithm have to be specifically adapted to be able to run on GPUs? Couldn't algorithms be written from higher level building blocks? Eg. implement map, reduce and scanl in CUDA/OpenCL, and have a default parameter like backend='cpu' than can be set to 'gpu'.

krastanov 7y ago
Julia is doing exactly that (even better than setting a flag, it is based on subtypes). Algorithms are usually written to work with `AbstractArray`. The usual RAM/CPU arrays (i.e. the type `Array`) is a subtype. The CUDA array which is stored on the GPU is also a subtype. Most general purpose code does not care which type you use.
- ViralBShah 7y ago
  Definitely look at the whole JuliaGPU ecosystem, and the CUDAnative.jl package. There are two big questions. First, how do you get code generation for targets such as GPUs from a dynamic high level language such as Julia, and that is what CUDAnative.jl achieves. The second question is what abstractions should Julia present to programmers so that increasingly larger codebases can automatically leverage GPUs. Packages such as CuArrays.jl are early answers, but there is much work to be done here.
ChrisRackauckas 7y ago
This is precisely correct. You just have to make sure in the thousands of lines of code that you only use higher level building blocks and make no assumptions greater than what a GPU-based array allows. GPU-based arrays don't even want to allow you to do scalar indexing for example, since every scalar indexing operation is a data transfer back to the CPU. So, with those rules in mind, getting the stiff ODE solvers to automatically work with GPU-based arrays amounted to: 1. Making sure every operation in the solvers used a higher level building block (@..) that automatically built fused broadcast kernels, and CuArrays.jl then supplies the broadcast overloads. That makes all of the non-stiff methods work. 2. Change the default choice of Jacobian to "dense array of size NxN" to "zero'd outer product of input array defined via broadcast". This translates to `J = false .* z .* z'` as a trick to make a Jacobian match whatever array type the array type library thinks is best, so J will be a GPUArray while DiffEq knows nothing about GPUs (also works for things like MPI, but GPUs is easier for people to reason about these days). So a one line change to the stiff ODE solvers made them all GPU compatible! But... 3. There's a small bug, or at least missing feature, in the GPU libraries that their QR doesn't implement backsolve. Technically the user can work around that because we allow someone to pass a linear solver routine with the `linsolve` argument, but it's easier to just handle it for everyone in the defaults. So using Requires.jl, if you have CuArrays installed, then we add a linsolve routine to the defaults for you, defined by https://github.com/JuliaDiffEq/DiffEqBase.jl/blob/master/src/init.jl#L104-L110 https://github.com/JuliaDiffEq/DiffEqBase.jl/blob/master/src... . Technically this can get upstreamed to CuArrays.jl and get deleted from DiffEqBase, but for now it'll live there so that everything is easier on users. So yes, make it high level, avoid assuming things like the Jacobian live on the CPU, and add one missing method to the GPU array library, and now when the user passes `u0 = CuArray(...)` for the initial condition, the algorithm compiles a version which does all state operations on the GPU (and all logic on the CPU, so it's a pretty good mix!).

Certhas 7y ago

Incredible work. Thank you everyone for the great release. This massively improves our usecase (sparse matrix equations) and I am looking forward to playing with it.

ChrisRackauckas 7y ago
Oh boy, just wait until next release! Sparse matrices are one of our big focuses for the summer: https://github.com/pkj-m/SparseDiffTools.jl/issues/2 https://github.com/pkj-m/SparseDiffTools.jl/issues/2

vsskanth 7y ago

Anyone tried using Julia's solvers to run FMI ModelExchange binaries ?

ChrisRackauckas 7y ago
Not yet, but we are looking into it. I'm hoping to get some funding for this, since it's a no-research software dev project though, which would be great for a just-out-of-undergrad student looking to transition to software dev, hence need money to pay such a student.
- vsskanth 7y ago
  Thank you.

mlevental 7y ago

this happens literally every single time i try to hop on the julia bang wagon: run example -> example breaks. 1. fresh install of juno (1.1.0) 2. add latest packages 3. something goes wrong. this time it's p = destructure(model) from example page gives UndefVarError: destructure not defined top-level scope at none:0 https://i.imgur.com/8fCpVrR.png https://i.imgur.com/8fCpVrR.png how the hell does anyone get anything when the codebase is shifting this fast!?

ChrisRackauckas 7y ago
That's my bad. Somehow an earlier draft of the post was used. The post was updated with working code. I guess that's what happens when you lose internet and push to a Github repo from your phone :P.
- mlevental 7y ago
  :) well thanks! edit: UndefVarError: JacVecOperator not defined :( also this https://gist.github.com/makslevental/fe9144d196cade685d333cd862c012c0 https://gist.github.com/makslevental/fe9144d196cade685d333cd... i realize it's just a warning but it shows another undef error
  - ChrisRackauckas 7y ago
    Are you on the latest versions? See what running a package update does. If that doesn't work, file an issue with your package Manifest so we can track this down. Thanks.