Identity crisis: performance with context

[ChrisRackauckas]

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 (?).

[chakravala]

Thanks for this explanation, I am in agreement with what is being proposed. You mentioned the ordering of variables, by the way Reduce comes with facilities for that using the order command, which is one of the functions included with the current development branch.

The possibilities with Reduce.jl are still in their beginning phases for Julia, I am still working on improvements to speed up the interface and to handle special cases. If I knew more details about what your requirements are, then I have a better idea of what to focus on. I'm also currently working on a new package that I have not released yet to help with code generation for basis functions. It comes with some innovative new analysis tools to get information about floating point error too, for which I am writing a paper that I will submit to the Julia Special Edition from the AES journal.

Also, for unique names in Julia you can use gensym() if I'm not mistaken? I've used this in some of my polynomial error analysis code that I'm currently working on. I'm still creatively working out the details, but I am experimenting with efficient methods for selecting optimal syntax trees for use as the generated code in the basis functions of approximations. Basically, I am looking at various charactersitic values calculated from the AST provided in order to select for the optimal form.

Please let me know in what ways I could get involved with this.

[ChrisRackauckas]

If I knew more details about what your requirements are, then I have a better idea of what to focus on.

The discussion that led to the repo is SciML/DifferentialEquations.jl#261 (this starts late into the process, but good enough). The README examples are pretty informative of what we want (https://github.com/JuliaDiffEq/ModelingToolkit.jl#introduction-by-examples), and then the extras are things like computing Jacobians and then inverting matrices of these symbolic expressions. That's in the tests.

https://github.com/JuliaDiffEq/ModelingToolkit.jl/blob/master/test/system_construction.jl this test file pretty much has what we want for now, and from there we want to build tooling like #18, #39, #55

[chakravala]

It's not necessarily the type system that is causing the instability, the symplify_constant method is recursively calling itself, and rewrites itself until it stops changing. Thus 99% of the time consumption is likely due to the fact that the simplify algorithm takes a huge amount of steps before settling down, especially when the input tree already has an enormous number of vertices.

This issue will be resolved in my next pull request.

[ChrisRackauckas]

Yeah, that is a pretty awful way to do simplification. It should likely be replaced by a call to a CAS which more appropriately implements a simplify method.

[ChrisRackauckas]

Need to do some benchmarks to see the final diff against ParameterizedFunctions.jl.

[ChrisRackauckas]

using ModelingToolkit
using Base.Test

function f3()
  # 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[],[σ,ρ,β])
  ModelingToolkit.generate_ode_function(de)
  jac_expr = ModelingToolkit.generate_ode_jacobian(de)
  ModelingToolkit.generate_ode_iW(de)
end

@time f3()
@time f3()
@time f3()

using Reduce

@time f3()
@time f3()
@time f3()

using ParameterizedFunctions
@time @ode_def_nohes Lorenz begin
  dx = σ*(y-x)
  dy = x*(ρ-z)-y
  dz = x*y - β*z
end σ ρ β

  5.412143 seconds (3.79 M allocations: 192.504 MiB, 1.00% gc time)
  0.464278 seconds (1.02 M allocations: 46.824 MiB, 2.46% gc time)
  0.474265 seconds (1.02 M allocations: 46.824 MiB, 3.14% gc time)
INFO: Recompiling stale cache file C:\Users\Chris\.julia\lib\v0.6\Reduce.ji formodule Reduce.
Reduce (Free CSL version, revision 4534), 05-Apr-18 ...
INFO: Recompiling stale cache file C:\Users\Chris\.julia\lib\v0.6\ReduceLinAlg.ji for module ReduceL
inAlg.
 10.116556 seconds (5.66 M allocations: 253.838 MiB, 1.60% gc time)
  1.058139 seconds (314.03 k allocations: 16.334 MiB, 0.64% gc time)
  1.047559 seconds (314.07 k allocations: 16.207 MiB, 0.58% gc time)
INFO: Precompiling module ParameterizedFunctions.
  0.034925 seconds (27.91 k allocations: 1.730 MiB)
  0.028652 seconds (20.30 k allocations: 1.268 MiB)
  0.028483 seconds (20.30 k allocations: 1.270 MiB)

Current results:

• ModelingToolkit is still about 1 order of magnitude slower than ParameterizedFunctions
• Reduce is slowing down the computation? Type-instability?

[ChrisRackauckas]

@eeshan9815

[ChrisRackauckas]

A lot of the v0.7 compiler updates were helpful. A newer test script is:

using DiffEqBase, LinearAlgebra
using ModelingToolkit
using Test

function f3()
  # 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[],[σ,ρ,β])
  ModelingToolkit.generate_ode_function(de)
  jac_expr = ModelingToolkit.generate_ode_jacobian(de)
  #ModelingToolkit.generate_ode_iW(de)
end

@time f3()
@time f3()
@time f3()

using Profile
@profile for i in 1:100 f3() end
Juno.profiler()

using ParameterizedFunctions
@time f = @ode_def_nohes begin
  dx = σ*(y-x)
  dy = x*(ρ-z)-y
  dz = x*y - β*z
end σ ρ β

macro ode_def_nohes2(name,ex,params...)
  opts = Dict{Symbol,Bool}(
  :build_tgrad => true,
  :build_jac => true,
  :build_expjac => false,
  :build_invjac => false,
  :build_invW => false,
  :build_hes => false,
  :build_invhes => false,
  :build_dpfuncs => false)
  name isa Expr ? ode_def_opts(gensym(),opts,name,params...) :
    ode_def_opts(name,opts,ex,params...)
end

@time f = @ode_def_nohes2 begin
  dx = σ*(y-x)
  dy = x*(ρ-z)-y
  dz = x*y - β*z
end σ ρ β

which shows

  0.001384 seconds (6.50 k allocations: 531.906 KiB)
  0.001642 seconds (6.50 k allocations: 531.906 KiB)
  0.001334 seconds (6.50 k allocations: 531.906 KiB)
  0.040115 seconds (66.50 k allocations: 3.344 MiB, 23.43% gc time)
  0.029821 seconds (66.50 k allocations: 3.344 MiB)
  0.029640 seconds (66.50 k allocations: 3.344 MiB)

and ModelingToolkit.jl is faster for calculations with Jacobians but no iWs. The profile shows that almost all of the time is spent in my admittedly bad simplification routine, so we may want a specialist to come in and help that part out, but I'm not worried about speed here anymore. The iW speed is really bad, about an order of magnitude worse than ParameterizedFunctions, but that's because the Operation type is bad inside of the generic linear algebra routines and that was probably an unnecessary abuse anyways. We can always fix that (and specialize on small matrices) in the future.

One other thing to mention is that right now we use Vector{Expression} with Expression being an abstract type. In the future when JuliaLang/julia#269 goes through or just by using

struct M
  x::Vector{Union{Nothing,M}}
end
M([M([nothing])])

we can make use of small union optimizations and handle this even better. Because of the inference we do get on this, it will be better and more flexible than using an AST so I think this is a good reason to continue on this route.