Tumble-detection (#19)

* routines for turn detection and run statistics

* turn detection tests

* example exponential run time distribution

examples/Analysis/run_distribution.jl

@@ -0,0 +1,23 @@
+using MicrobeAgents
+using Distributions
+using Plots
+
+L = 1000
+extent = (L,L)
+dt = 0.01
+
+model = UnremovableABM(Microbe{2}, extent, dt)
+n = 500
+Drot = 0.1
+τ = 1.0
+for i in 1:n
+    add_agent!(model; turn_rate=1/τ, rotational_diffusivity=Drot)
+end
+
+nsteps = 5000
+adf, _ = run!(model, nsteps; adata=[:vel])
+
+ᾱ = 4 * sqrt(2*Drot*dt)
+runs = rundurations(adf, dt; threshold_angle=ᾱ)
+histogram(vcat(runs...), bins=dt/2:10dt:5τ, normalize=:pdf, lab="Simulation", lw=0)
+plot!(5dt:dt:5τ, t -> pdf(Exponential(τ), t), lw=4, lab="exp(-t/τ)/τ")

src/MicrobeAgents.jl

@@ -57,5 +57,6 @@ include("analysis/postprocess.jl")
 include("analysis/msd.jl")
 include("analysis/correlation_functions.jl")
 include("analysis/drift.jl")
+include("analysis/turn_detection.jl")
 
 end

src/analysis/turn_detection.jl

@@ -0,0 +1,89 @@
+export rundurations, mean_runduration, mean_turnrate, detect_turns, has_turned
+
+
+"""
+    rundurations(adf, Δt; threshold_angle=0.0)
+Evaluate the duration of all runs observed during a simulation.
+The dataframe `adf` should contain the field `:vel`.
+`Δt` is the integration timestep of the simulation.
+`threshold_angle` defines the threshold (in radians) to define a reorientation.
+The function returns a vector of run durations for each microbe.
+
+Also works if `adf` is an `AbstractMatrix`, with unit velocity vectors
+for each microbe sorted by column, or if `adf` is an `AbstractVector`
+(unit velocity vectors for a single microbe).
+"""
+function rundurations(adf, Δt; threshold_angle=0.0)
+    turns = detect_turns(adf; threshold_angle)
+    map(col -> diff([0; findall(col)]) .* Δt, eachcol(turns))
+end
+
+
+"""
+    mean_runduration(adf, Δt; threshold_angle=0.0)
+Evaluate the mean run duration sampled by the simulation.
+The dataframe `adf` should contain the field `:vel`.
+`Δt` is the integration timestep of the simulation.
+`threshold_angle` defines the threshold (in radians) to define a reorientation.
+
+Also works if `adf` is an `AbstractMatrix`, with unit velocity vectors
+for each microbe sorted by column, or if `adf` is an `AbstractVector`
+(unit velocity vectors for a single microbe).
+"""
+mean_runduration(adf, Δt; threshold_angle=0.0) = 1.0 / mean_turnrate(adf, Δt; threshold_angle)
+
+"""
+    mean_turnrate(adf, Δt; threshold_angle=0.0)
+Evaluate the mean turn rate sampled by the simulation.
+The dataframe `adf` should contain the field `:vel`.
+`Δt` is the integration timestep of the simulation.
+`threshold_angle` defines the threshold (in radians) to define a reorientation.
+
+Also works if `adf` is an `AbstractMatrix`, with unit velocity vectors
+for each microbe sorted by column, or if `adf` is an `AbstractVector`
+(unit velocity vectors for a single microbe).
+"""
+function mean_turnrate(adf, Δt; threshold_angle=0.0)
+    turns = detect_turns(adf; threshold_angle)
+    count(turns) / (length(turns) * Δt)
+end
+
+"""
+    detect_turns(adf; threshold_angle=0.0)
+Detect reorientations in the microbe trajectories in `adf`.
+Requires `adf` to have a `:vel` field containing the unit-norm velocity vector.
+In order to be detected, a reorientation must be larger than the `threshold_angle`
+(in radians); the threshold defaults to 0.
+
+Also works if `adf` is an `AbstractMatrix`, with unit velocity vectors
+for each microbe sorted by column, or if `adf` is an `AbstractVector`
+(unit velocity vectors for a single microbe).
+"""
+function detect_turns(adf; threshold_angle=0.0)
+    vel = vectorize_adf_measurement(adf, :vel)
+    detect_turns(vel; threshold_angle)
+end
+
+function detect_turns(vel::AbstractMatrix; threshold_angle=0.0)
+    mapslices(u -> detect_turns(u; threshold_angle), vel, dims=1)
+end
+
+function detect_turns(vel::AbstractVector; threshold_angle=0.0)
+    a₀ = UnitRange(1, length(vel)-1)
+    a₁ = UnitRange(2, length(vel))
+    has_turned.(view(vel,a₀), view(vel,a₁); threshold_angle)
+end
+
+"""
+    has_turned(v₁, v₂; threshold_angle=0.0)
+Determine if `v₂` is rotated w.r.t `v₁` by comparing the angle between
+the two vectors with a `threshold_angle` (defaults to 0).
+"""
+has_turned(v₁, v₂; threshold_angle=0.0) = safe_acos(dot(v₁,v₂)) > abs(threshold_angle)%π
+
+"""
+    safe_acos(x)
+Numerically safe version of `acos`.
+Useful when `x` exceeds [-1,1] due to floating point errors.
+"""
+safe_acos(x::T) where T = x ≈ 1 ? zero(x) : x ≈ -1 ? T(π) : acos(x)

test/analysis.jl

@@ -27,6 +27,38 @@ using LinearAlgebra: norm
         end
     end
 
+    @testset "Turn detection and run statistics" begin
+        dt = 0.1
+        L = 100
+        for D in 1:3
+            extent = ntuple(_ -> L, D)
+            model = StandardABM(Microbe{D}, extent, dt)
+            add_agent!(model; turn_rate=Inf) # turns every step
+            nsteps = 10
+            adf, _ = run!(model, nsteps; adata=[:vel])
+            @test mean_turnrate(adf,dt) ≈ 1/dt
+            @test mean_runduration(adf,dt) ≈ dt
+            @test detect_turns(adf) == BitMatrix(ones(nsteps,1))
+            @test rundurations(adf,dt) == [repeat([dt],nsteps)]
+
+            D == 1 && continue # only test rotational diffusion for D>1
+            model = StandardABM(Microbe{D}, extent, dt)
+            # never turns but makes small deviations due to brownian noise
+            Drot = 0.1
+            add_agent!(model; turn_rate=0.0, rotational_diffusivity=Drot)
+            nsteps = 10
+            adf, _ = run!(model, nsteps; adata=[:vel])
+            # detects turns at each step if threshold is not set
+            @test detect_turns(adf) == BitMatrix(ones(nsteps,1))
+            # set threshold at 4σ to ignore rotational diffusion
+            threshold_angle = 4 * sqrt(2*Drot*dt)
+            @test detect_turns(adf;threshold_angle) == BitMatrix(zeros(nsteps,1))
+            @test mean_turnrate(adf,dt;threshold_angle) == 0.0
+            @test mean_runduration(adf,dt;threshold_angle) == Inf
+            @test rundurations(adf,dt;threshold_angle) == [Float64[]]
+        end
+    end
+
     @testset "Autocorrelation function" begin
         for D in 1:3
             dt = 0.1