openBF: 1D blood flow simulations with Julia

[alemelis]

Hello there, I’d like to introduce to the Julia community my openBF package.

openBF is a finite-volume solver for one-dimensional blood flow equation in networks of elastic (we are working on visco-elastic behaviour as well) arteries. It has a number of useful features and we hope it will be used for research as well as employed in clinical contex.

I wrote a small description of it on Medium [link]
here it is the GitHub repository [link]
and we have also an official website to host the documentation [link]

I’ve been developing it throughout my phd and it is currently used in a number of research project (mainly in my university), but I think that others may find it useful/interesting.

I am now re-organising the code to make it compliant with official guidelines as I’d like to submit it to METADATA. Also, I’d like to improve the solver performance, so any contribution is well accepted; enjoy!

[rdeits]

Thanks for sharing this!

One thing I noticed looking through your code is that there seem to be quite a few places where you construct very small arrays or matrices, such as https://github.com/INSIGNEO/openBF/blob/master/src/conjunctions.jl#L208 and https://github.com/INSIGNEO/openBF/blob/master/src/conjunctions.jl#L333 If you’re constructing those small arrays in a loop (and not modifying them), then you might benefit from using SVector and SMatrix from https://github.com/JuliaArrays/StaticArrays.jl

For example:

julia> using StaticArrays

julia> using BenchmarkTools

julia> @benchmark [1, 2, 3]
BenchmarkTools.Trial:
  memory estimate:  112 bytes
  allocs estimate:  1
  --------------
  minimum time:     37.542 ns (0.00% GC)
  median time:      41.233 ns (0.00% GC)
  mean time:        48.009 ns (5.93% GC)
  maximum time:     1.055 μs (81.68% GC)
  --------------
  samples:          10000
  evals/sample:     991

julia> @benchmark SVector(1, 2, 3)
BenchmarkTools.Trial:
  memory estimate:  0 bytes
  allocs estimate:  0
  --------------
  minimum time:     1.738 ns (0.00% GC)
  median time:      1.856 ns (0.00% GC)
  mean time:        1.969 ns (0.00% GC)
  maximum time:     12.985 ns (0.00% GC)
  --------------
  samples:          10000
  evals/sample:     1000

SVectors can also be efficiently stored inline in an array or a struct, so they’re great if you’re storing hundreds or thousands of small vectors. The only downside is that they’re necessarily immutable. There’s a mutable MVector also in StaticArrays.jl, although it’s more expensive than the immutable variant.

[alemelis]

Hey thanks for pointing me to StaticArrays!
I have a lot of arrays/matrices that are simply initialised at the beginning and I see how I can make them all immutable.

[kristoffer.carlsson]

In JuAFEM.jl (finite element library), we use https://github.com/KristofferC/Tensors.jl which is quite similar to StaticArrays (it was basically a wrapper over it once upon a time). The point I want to make is that using something like StaticArrays can have tremendous performance impact over constantly allocating fat Julia Arrays. It also allows for cleaner code since there is no need to preallocate output and use in-place operations. The following is an exhibit from some very hot code in JuAFM:

fecv_J = zero(Tensor{2,dim})
for j in 1:n_geom_basefuncs
    fecv_J += x[j] ⊗ cv.dMdξ[j, i]
end
detJ = det(fecv_J)
cv.detJdV[i] = detJ * w
Jinv = inv(fecv_J)
for j in 1:n_func_basefuncs
    cv.dNdx[j, i] = cv.dNdξ[j, i] ⋅ Jinv
end

Since nothing is mutated, we can use nice infix operators which is just great imo.

[Paulo_Jabardo]

Very cool. I will definetely take a close look! I am interested in wave propagation along pressure tap tubing. Working with water in flexible tubes is most probably directly applicable. Air (or other gases) where the waves are due to compressibility (entirely or partly) is a little more complex but maybe rewriting the differential equations carefully, the same basic solver could also be used.

Paulo

[alemelis]

@rdeits and @kristoffer.carlsson thanks for pointing me to static arrays.
I’m not 100% I am using them correctly.

In my code I have several Array{Float64, 1} that are initialised inside a struct. The values in those arrays are not going to be changed, so it seems there is room for making them immutable (am I right?).

In the current implementation, I have something like this

struct MyType
    a :: Array{Float64,1}
end

function useArray(M::Int)
    b = [i*i for i = 1:M]
    return MyType(b)
end

julia> a = useArray(10)
MyType([1.0, 4.0, 9.0, 16.0, 25.0, 36.0, 49.0, 64.0, 81.0, 100.0])

The StaticArrays test doesn’t work in the same way though.

struct MyTypeStatic
    a :: AbstractArray # I guess it should be SVector{M,Float64}
end

function useStaticArrays(M::Int)
    static_b = @SArray [i*i for i = 1:M]
    return MyTypeStatic(static_b)
end

By simply running this, I get ERROR: UndefVarError: M not defined. Does it mean that I should know the final size of the static array I want to allocate at compilation time?
By defining M before the function, I manage to make it work…

M = 10
function useStaticArrays(M::Int)
    static_b = @SArray [i*i for i = 1:M]
    return MyTypeStatic(static_b)
end

julia> useStaticArrays(10)
MyTypeStatic([1.0, 4.0, 9.0, 16.0, 25.0, 36.0, 49.0, 64.0, 81.0, 100.0])

This wouldn’t be a problem in parts of the code where I have to allocate small fixed size arrays. However, the size of a large part of the static-wannabe-arrays is computed at the beginning of the simulation and is case-specific, i.e. cannot be fixed to a certain value beforehand.

I haven’t looked at Tensors.jl yet, but I’m curious to check my understanding on StaticArrays anyway.

[alemelis]

Hi Paulo, I’m glad you found some similarities. The hyperbolic system at the base of openBF can be used for various applications as shallow-water equations, traffic models, and any other physical system presenting shock-like behaviour. The first version of oBF solver was based on Godunov’s scheme (first-order in time and space) which was originally developed to solve the shock-tube problem in gas-dynamics. Please keep me update if you decide to give oBF a try: it would be very interesting (from my side) to expand the solver capabilities and make it useful in other research fields.

[rdeits]

Does it mean that I should know the final size of the static array I want to allocate at compilation time?

You can construct an SVector whose size depends on a value like M in your code, but constructing it will be more expensive because the type of the returned SVector will depend on the value of M. For example, you can do:

SVector([i*i for i in 1:M])

but the compiler won’t be able to infer the type of the result. If the value of M is really something that can’t be known ahead of time, then it’s possible that it’s not worth using a static array for this particular case.

If you want to make it easier to play around with different array implementations, you could parameterize your type on the kind of array it contains, e.g.:

struct MyType{A <: AbstractArray}
  a::A 
end

which will let you use any kind of array you want without changing your type definition. Then you can try static vs. non-static arrays and see how your program’s performance changes.

[PetrKryslUCSD]

This looks interesting. Perhaps you have a paper in mind? If so, consider [ANN] Special issue on the Julia programming language

[alemelis]

I see, therefore I may see some benefit in using SArray when allocating a lot of arrays of fixed sized.

(I’m learning more in two days over discourse than in the previous four years, thanks!)

[sdanisch]

@alemelis

However, the size of a large part of the static-wannabe-arrays is computed at the beginning of the simulation and is case-specific, i.e. cannot be fixed to a certain value beforehand.

As long as your code has this structure, you’re good:

function expensive_long_running(fixed_size_over_run)
     ...
end

function setup_sim()
     N = calculate_dims()
     expensive_long_running(N)
end

You can write the above, so that you can force Julia to specialize code in expensive_long_running do your N.
This means, that for each N you have additional compilation costs, which can be a trade off well worth it.

You can write the above as:

struct MyTypeStatic{T <: AbstractArray} 
    a::T
end
# Or maybe just:
struct MyTypeStatic{N, T}
    a::SVector{N, T}
end
# move N into type domain, to make it statically inferable by julia
function expensive_long_running(n::Val{N}) where N
    x = MyTypeStatic(SVector(ntuple(i-> i * i, n))) # there is a fast ntuple(f, Val{N})
    while not_done
        x = simulate_stuff(x)
    end
end
function setup_sim()
    N = calculate_dims()
    # inside setup_sim Val{N} won't be inferrable by the compiler, so 
    # so setup_sim itsel won't be as performant - but it will call 
    # a `expensive_long_running` specialized to the value of N!
    expensive_long_running(Val{N}())
end

# Alternatively:

function expensive_long_running(v::SVector{N}) where N
    v_other = SVector(ntuple(i-> i * i, Val{N}())) 
    x = MyTypeStatic(v_other) # or just use v directly?
    while not_done
        x = simulate_stuff(x)
    end
end

function setup_sim()
    N = calculate_dims()
    v = SVector([i*i for i in 1:N]) # the slow creation, that will set the example for following SVectors
    expensive_long_running(v)
end

Take away: as long as the dimensions of your arrays stay fixed across a long enough run of your program, you can compile a faster version specialized to your dimensions. You will need to amortize the time lost by compilation times, though

[ChrisRackauckas]

@sdanisch

You will need to amortize the time lost by compilation times, though

Well, for actual applications of differential equations I really wouldn’t worry too much about this. I mean, in this example the simulation already lasts 20 minutes, so extreme optimizations which kick compile time up to a minute still make sense. Hyper optimization with static vectors seems to make the compile time “seconds more” or “tens of seconds more”, so if that cuts runtime by a factor like 2x, that pays off big time even if there’s a single run.

[sairus7]

Do you plan to add some breath and blood oxygenation models to estimate breath disorders?

Or do you know any working breath models in other languages?

[alemelis]

I haven’t looked into this yet, but I’m planning to add several other physiological processes. For example, I’d like to include a feedback loop to represent baroreceptors action and the autoregulation effect. One idea is to couple openBF with a model like this. However, this is a very long term goal

[sairus7]

I think, one of the promising directions in this field is estimation of real biological signals (from medical sensors) with dynamic physiological models. This can be used to improve diagnostics.

[alemelis]

bumping this thread to announce the latest paper about openBF

perhaps of interest to this community is this performance benchmark

It’s to take with a pinch of salt (different langs to implement different numerical schemes and run on different hardware), but I’d say it’s a good testament of

• ease of development in Julia
• performance for free you get by simple implementation (I tried to keep the code easy to understand vs. over optimising it)

[TheCedarPrince]

This is really interesting! Would you want to present about this work sometime at one of our JuliaHealth monthly meetings @alemelis?

[alemelis]

Hey, thanks for the interest…DMing you!