MultiScaleModels Profiling about multiscalearrays.jl

[ChrisRackauckas]

v0.6 DiffEq got a lot faster while this stayed the same:

prob = ODEProblem(f,em,(0.0,1500.0))
@benchmark sol1 = solve(prob,Tsit5(),save_everystep=false)
prob = ODEProblem(f,em[:],(0.0,1500.0))
@benchmark sol2 = solve(prob,Tsit5(),save_everystep=false)


BenchmarkTools.Trial: 
  memory estimate:  24.44 KiB
  allocs estimate:  323
  --------------
  minimum time:     393.737 μs (0.00% GC)
  median time:      432.086 μs (0.00% GC)
  mean time:        452.956 μs (0.88% GC)
  maximum time:     4.695 ms (89.11% GC)
  --------------
  samples:          10000
  evals/sample:     1


BenchmarkTools.Trial: 
  memory estimate:  103.53 KiB
  allocs estimate:  1774
  --------------
  minimum time:     27.596 ms (0.00% GC)
  median time:      28.951 ms (0.00% GC)
  mean time:        28.998 ms (0.00% GC)
  maximum time:     32.709 ms (0.00% GC)
  --------------
  samples:          173
  evals/sample:     1

prob = ODEProblem(f,em,(0.0,1500.0))
@benchmark sol1 = solve(prob,Tsit5())
prob = ODEProblem(f,em[:],(0.0,1500.0))
@benchmark sol2 = solve(prob,Tsit5())

BenchmarkTools.Trial: 
  memory estimate:  13.31 MiB
  allocs estimate:  123635
  --------------
  minimum time:     93.224 ms (0.00% GC)
  median time:      102.820 ms (6.55% GC)
  mean time:        100.983 ms (4.83% GC)
  maximum time:     110.067 ms (12.34% GC)
  --------------
  samples:          50
  evals/sample:     1


BenchmarkTools.Trial: 
  memory estimate:  1.21 MiB
  allocs estimate:  6499
  --------------
  minimum time:     939.406 μs (0.00% GC)
  median time:      1.068 ms (0.00% GC)
  mean time:        1.208 ms (7.40% GC)
  maximum time:     4.706 ms (52.20% GC)
  --------------
  samples:          4123
  evals/sample:     1

Don't know what to say about that. Probably could use some optimizations. Or the "full broadcast" integrators would probably do well, since this is using the indexing versions right now.

The major change that could have caused this is that now the "non-user facing cache variables" also match the type that the user gives. Before they were transformed into contiguous arrays since they were not shown to the user. That had some weird side-effects though, and for example slows down "add_daughter" types of events. So it's somewhat a wash... broadcast will likely be the savior here.

[zahachtah]

Do you anticipate 0.7 or 1.0 do improve the 2-3x slower performance?

[ChrisRackauckas]

It'll fix it. The detail is that linear indexing is slow because it has to do a binary search, but the built in broadcast is fast. A v0.6 bug prevented broadcast in most ODE/SDE methods: SciML/OrdinaryDiffEq.jl#106

However, this bug was fixed in Julia v0.7: JuliaLang/julia#22255 . So in the v0.7 upgrade we will be making all of those algorithms internally use broadcasting which will get rid of the slow indexing that MultiScaleArrays is hitting (and it also has other nice side effects, like all of the RK methods will be GPU-compatible!).

With ArrayPartitions from RecursiveArrayTools.jl we found that grouping small broadcasts can actually be faster than contiguous loops because of smart cache handling, so there is a chance this can do extremely well. One issue this library will still have is that copying a MultiScaleArray is more expensive than a standard Array, so if you have a lot of saving going on that will be more expensive. But together I wouldn't be surprised if MultiScaleArrays is a small cost (1.3x?) or just a wash (<1.1x performance loss). The 2x-3x shouldn't exist after these changes. I am very very very happy that compiler change happened in Base.

[zahachtah]

Oh, thats very cool! I may at times have very large arrays (10000+), do you think threading the ode will be efficient too in 0.7?

[ChrisRackauckas]

Yes, for arrays of that size it would be a good idea. I want to create a package which makes broadcast multithreaded via a wrapper. That's pretty simple on v0.7, so I was waiting on that as the solution.

[zahachtah]

Hi again, read this thread (https://discourse.julialang.org/t/agent-based-modeling-in-julia/12431/9) and realized its relevant to my case. I have different types as leaves in multiscalearrays. I haven’t benchmarked yet but i guess i would have issues with mixed types? Would https://github.com/tkoolen/TypeSortedCollections.jl solve that?

[ChrisRackauckas]

Or just have tuple types for holding the nodes: #26 . TypeSortedCollections.jl will work similarly. ArrayPartition works as well. These all require indexing with literals or broadcasting to get the type-stable operations though.