Implementing the Social Force Model from Helbing 1995

Hey everyone, I’m trying to implement the social force model for pedestrian dynamics, which models pedestrians using physics. I’m not sure what the best course of action is, specifically in terms of which packages I could make use of. I know that it would be possible to stay strictly within the DifferentialEquations.jl package by doing everything with Arrays, but I wanted to ask whether there were some convenient functions or features from other packages in the JuliaDynamics ecosystem ( DynamicalSystems.jl , Agents.jl ) that I could use to make the code more readable and modular.

I did some testing* earlier and it seemed that using structs in DifferentialEquations.jl as opposed to arrays had a time penalty of about 2x, which I’m not sure justifies the nicer interface, but if someone can look at the functions and see something obvious that I’m doing wrong I’d appreciate it.

Edit: removed the link to my repo and made it private by request of my supervisor. I will make the repo public after submitting in March 2022

Interesting. Changing your MWE to

using BenchmarkTools

f(x) = 2x^2

struct Agent{T}
    x::T
    y::T
    v::T
    w::T
end

import Base.similar, Base.copy;
copy(u::Agent{T}) where T = Agent{T}(u.x,u.y,u.v,u.w);
similar(u::Agent{T}) where T = copy(u);

function test!(du::T, u::T, N::Int) where {T<:Array{<:Real}}
    for i=1:N
        du[1, i] = f(u[1, i])
        du[2, i] = f(u[2, i])
        du[3, i] = f(u[3, i])
        du[4, i] = f(u[4, i])
    end
end

function test!(du::T, u::T, N::Int) where {T<:Array{<:Agent}}
    for i=1:N
        du[i] = Agent(f(u[i].x), f(u[i].y), f(u[i].v), f(u[i].w))
    end
end


function wrappedAllocations(f, args...)
    f(args...)
    res = @timed f(args...)
    return res
end

function main()
    N = 1000000;

    u0 = rand(4, N); # Array{Float64, 2}
    U0 = [Agent(u0[:,i]...) for i=1:N]; # Array{Agent{Float64}, 1}

    du = similar(u0); # Array{Float64, 2}
    dU = similar.(U0); # Array{Agent{Float64}, 1}

    print("Arrays: ")
    @btime test!($du, $u0, $N)
    print("Agents: ")
    @btime test!($dU, $U0, $N)
end

main();

yields

Arrays:   3.408 ms (0 allocations: 0 bytes)
Agents:   3.068 ms (0 allocations: 0 bytes)

for me. I believe there a way better optimizations available for immutable structures.

1 Like

Hmm, when I do the same it starts allocating like crazy:

Arrays:   5.357 ms (0 allocations: 0 bytes)
Agents:   15.557 ms (1000000 allocations: 45.78 MiB)

Why does yours not allocate a new agent every time it calls Agent?

Immutable or versioninfo()?

Julia Version 1.7.0-rc3
Commit 3348de4ea6 (2021-11-15 08:22 UTC)        
Platform Info:
  OS: Windows (x86_64-w64-mingw32)
  CPU: Intel(R) Core(TM) i7-10710U CPU @ 1.10GHz
  WORD_SIZE: 64
  LIBM: libopenlibm
  LLVM: libLLVM-12.0.1 (ORCJIT, skylake)
Environment:
  JULIA_NUM_THREADS = 6

It was indeed the immutable part that was causing the issue. I didn’t realise that you changed that as well.

1 Like

Hello,

Regarding the approach of @goerch : This is something I have considered a lot in Agents.jl: whether having agents as immutable types was worth it or not. Measuring performance benefits, if any, of immutable Agent type · Issue #421 · JuliaDynamics/Agents.jl · GitHub . The conclusion was that even though there were some measurable performance benefits, the huge downsides in convenience and intuitive rationalizing about code were simply too much to consider it useful in practice.

@jacobusmmsmit for your original question. First of all, I wasn’t even aware that DifferentialEquations.jl can work with vectors of arbitrary types. I think you would still need to define addition between vectors of agents, etc. Are you sure this “works”, or is it as easy as you think?

More importantly now. For your social force model, what you need to consider is the importance of accuracy. Is an Euler stepping approach accurate enough? Then you can go ahead with Agents.jl. If you need better accuracy for time evolution, you can consider DifferentialEquations.jl. No matter the level of accuracy, DiffEq will of course be always faster than Agents.jl. Things you should consider having a look at:

1 Like

No, but I’m going to go away and try before I reply any further…

For using DynamicalSystems.jl: This package is only useful if you want to do an analysis of differential equations that goes beyond integrating these equations. I am assuming that the entire point of the Social Force Model is to set up some differential equations and then solve them. In this scenario DynamicalSystems.jl offers nothing more besides wrapping DiffEq (and hence, just use DiffEq directly).

I’m not sure the functional programming community would like to agree? They usually prefer immutable to be able to proof properties of their programs.

Actually, looking the paper in more detail it seems like the continuous space of Agents.jl may be the best fit, or the more intuitive fit to be precise. You can also be inspired by the work done in this PR: add Force Based Motion Planning (FMP) by rmcsqrd · Pull Request #456 · JuliaDynamics/Agents.jl · GitHub on Force Based Motion planning.

(I need to get around to merging this PR…)

I dont think this really connects to to what I said about rationalizing about the code representation of an agent based model though. The in-place version with mutable types is more sensible and represents how one thinks about agent based models. When you move an agent you don’t “create a new one at a new location”, but you rather “change the location of an existing one”. If you drive your mom at the supermarket you don’t create a new mom at the supermarket while killing the one at home. I guess the functional programming community would generally agree with this line of thought.

Buuuut regardless, I think this is becoming too off-topic for this post. To not derail further, feel free to write more in the dynamics-bridged channel on slack or the issue link I posted in the original post!

2 Likes

Yes, this is correct. You have to define the relevant arithmetic functions for it to work which I should have seen coming.

1 Like

Molly.jl has functionality for this sort of thing with custom forces, though I don’t know the details of the model in question.