I’m writing a simulation that solves a system of N differential equations. I’m using fixed step Euler integration, so it’s always a simply for loop that iterates over many times. I have two implementations, and one is 4X faster than the other, and has 30% the allocations. I’d like to understand why:
# fast code
f1(V::Float64,dt) = some function of V
f2(V::Float64,dt) = some function of V
...
fN(V::Float64,dt) = some function of V
integrate(t::CustomType,dt::Float64)
var1::Float64
var2::Float64
...
varN::Float64
for i = 2:nsteps
var1 = f1(var1,dt)
...
varN = fN(varN,dt)
end
end
This is rigid code, in the sense that if I want to modify my system by omitting a few variables and corresponding functions, it is laborious to rewrite. Instead, I wrote these structs and associated functions. This is nicer conceptually, but runs much slower:
# slow code
abstract Type X
end
mutable struct X1 <: X
v1
end
function integrate(x::X1,dt::Float64)
x.v1 = f1(x.v1,dt)
end
mutable struct Wrapper
all_X::Array{X,1}
end
# now I can define my core integration loop nicely:
function integrate(w::Wrapper,dt)
for x in w.all_X
integrate(x,dt)
end
end
I really want to write code like the second case – it’s a lot more flexible, and prevents me from repeating code. But the second implementation is 4 times slower than the first. Any tips on why would be appreciated, as would help on how to get it to run faster.
This is a mutable struct which doesn’t concretely type its fields. This makes w.all_X not type stable, which then makes integrate(x,dt) dynamically dispatch. Instead, use a type parameter
mutable struct Wrapper{T<:X}
all_X::Array{T,1}
end
and it’ll be much better.
Another thing to mention is that the way you are writing the functions you are assuming the differential equations are all independent. This is not true in general nor is it true in most cases (though it may be for your case, I don’t know), so you should be cautious about this approach unless it’s for a very specific purpose.
I don’t see the abstract type ‘X’ references anywhere here. Should it be
mutable struct Wrapper{T<:X}
How do I create a new object of this type? I tried
w = Wrapper()
ERROR: MethodError: no method matching Wrapper()
w = Wrapper({X})
ERROR: syntax: { } vector syntax is discontinued
w = Wrapper{X}()
ERROR: MethodError: no method matching Wrapper{X}()
thanks a ton!
update
I managed to get this working, but without an inner constructor (I used an external constructor)
Sadly, there’s no change in the speed of the code.
Can you show what you have now? Are you sure everything is strictly typed and inferred? And for mutable structs, did you make sure you’re only mutating them and not creating them a bunch? (If you want to create them a bunch, just use a struct)
Here’s a zip file that has two implementations of my simulations, one with mutable structs and one without. The one without runs 5X faster than the one with.