Skip to content

Commit

Permalink
Update docs and use Literate.jl (#64)
Browse files Browse the repository at this point in the history
* update docs with Literate.jl

* reorganize sections

* update examples

* update docs

* temporary fix
  • Loading branch information
mastrof authored Dec 28, 2023
1 parent 0daee90 commit 0e3745c
Show file tree
Hide file tree
Showing 19 changed files with 1,205 additions and 218 deletions.
4 changes: 3 additions & 1 deletion docs/Project.toml
Original file line number Diff line number Diff line change
@@ -1,6 +1,8 @@
[deps]
Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4"
DocumenterTools = "35a29f4d-8980-5a13-9543-d66fff28ecb8"
Literate = "98b081ad-f1c9-55d3-8b20-4c87d4299306"
Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"

[compat]
Documenter = "0.27"
Documenter = "1"
78 changes: 67 additions & 11 deletions docs/make.jl
Original file line number Diff line number Diff line change
@@ -1,20 +1,76 @@
#=
cd(@__DIR__)
using MicrobeAgents
using Documenter
=#
# temporary
push!(LOAD_PATH, "../src/")
using Documenter, MicrobeAgents
ENV["JULIA_DEBUG"] = "Documenter"
CI = get(ENV, "CI", nothing) == "true" || get(ENV, "GITHUB_TOKEN", nothing) !== nothing
import Literate
using Plots

indir_base = joinpath(@__DIR__, "..", "examples")
sections = ("RandomWalks", "Chemotaxis")
outdir_base = joinpath(@__DIR__, "src", "examples")
indir = Dict(
section => joinpath(indir_base, section)
for section in sections
)
outdir = Dict(
section => joinpath(outdir_base, section)
for section in sections
)
rm(outdir_base; force=true, recursive=true) # clean up previous examples
mkpath(outdir_base)
for section in sections
mkpath(outdir[section])
for file in readdir(indir[section])
Literate.markdown(joinpath(indir[section], file), outdir[section]; credit=false)
end
end


# convert camelcase directory names to space-separated section names
function namify(s)
indices = findall(isuppercase, s)
if length(indices) <= 1
return s
else
s1 = s[indices[1] : indices[2]-1]
s2 = lowercase.(s[indices[2] : end])
return join((s1, namify(s2)), " ")
end
end

pages = [
"Home" => "index.md",
"Introduction" => ["firststeps.md", "randomwalks.md", "chemotaxis.md"],
"Examples" => [
[namify(section) => [joinpath.("examples", section, readdir(outdir[section]))...]
for section in sections]...
],
"Validation" => "validation.md",
"API" => "api.md"
]

makedocs(
sitename = "MicrobeAgents.jl",
authors = "Riccardo Foffi",
modules = [MicrobeAgents],
pages = [
"Home" => "index.md",
"Tutorial" => ["firststeps.md", "randomwalks.md", "chemotaxis.md"],
"Validation" => "validation.md",
"API" => "api.md"
],
pages = pages,
expandfirst = ["index.md"],
format = Documenter.HTML(
prettyurls = get(ENV, "CI", nothing) == "true"
)
prettyurls = CI,
),
warnonly = true,
)

deploydocs(;
repo = "github.com/mastrof/MicrobeAgents.jl"
)
if CI
deploydocs(;
repo = "github.com/mastrof/MicrobeAgents.jl.git",
target = "build",
push_preview = true
)
end
29 changes: 28 additions & 1 deletion docs/src/api.md
Original file line number Diff line number Diff line change
@@ -1,5 +1,32 @@
# API

## [Microbes](@id Microbes)
```@docs
AbstractMicrobe
Microbe
BrownBerg
Brumley
Celani
Xie
```

## [Motility](@id Motility)
```@docs
AbstractMotility
MotilityOneStep
MotilityTwoStep
RunTumble
RunReverse
RunReverseFlick
```

## [Utils](@id Utils)
```@docs
vectorize_adf_measurement
random_speed
random_velocity
```

## [Data analysis](@id Analysis)
```@docs
msd
Expand All @@ -10,4 +37,4 @@ detect_turns
rundurations
mean_runduration
mean_turnrate
```
```
56 changes: 56 additions & 0 deletions docs/src/examples/Chemotaxis/celani_gauss2D.md
Original file line number Diff line number Diff line change
@@ -0,0 +1,56 @@
```@meta
EditURL = "../../../../examples/Chemotaxis/celani_gauss2D.jl"
```

# Noisy chemotaxis towards Gaussian source

In this example we set up a static Gaussian source and observe the chemotactic behavior
of the `Celani` model, in the presence of sensing noise (via the `chemotactic_precision`).
Playing with the `chemotactic_precision`, it can be seen that the clustering of bacteria
at the source becomes stronger with decreasing noise (decreasing chemotactic precision).

````@example celani_gauss2D
using MicrobeAgents
using Plots
function concentration_field(pos, model)
C = model.C
σ = model.σ
p₀ = model.p₀
concentration_field(pos, p₀, C, σ)
end
concentration_field(pos, p₀, C, σ) = C * exp(-sum(abs2.(pos.-p₀))/(2*σ^2))
timestep = 0.1 # s
extent = ntuple(_ -> 1000.0, 2) # μm
space = ContinuousSpace(extent; periodic=false)
p₀ = extent./2 # μm
C = 1.0 # μM
σ = 100.0 # μm
properties = Dict(
:concentration_field => concentration_field,
:C => C,
:σ => σ,
:p₀ => p₀,
)
model = StandardABM(Celani{2}, space, timestep; properties)
foreach(_ -> add_agent!(model; chemotactic_precision=6.0), 1:300)
nsteps = 5000
adata = [position]
adf, = run!(model, nsteps; adata)
traj = vectorize_adf_measurement(adf, :position)
xmesh = range(0, first(spacesize(model)); length=100)
ymesh = range(0, last(spacesize(model)); length=100)
c = [concentration_field(p, p₀, C, σ) for p in Iterators.product(xmesh, ymesh)]
heatmap(xmesh, ymesh, c', cbar=false, ratio=1, c=:bone, axis=false)
x = getindex.(traj,1)[end-100:4:end, :]
y = getindex.(traj,2)[end-100:4:end, :]
a = axes(x,1) ./ size(x,1)
plot!(x, y,
lab=false, lims=(0,1000), lw=1, alpha=a
)
````

123 changes: 123 additions & 0 deletions docs/src/examples/Chemotaxis/linear_ramp.md
Original file line number Diff line number Diff line change
@@ -0,0 +1,123 @@
```@meta
EditURL = "../../../../examples/Chemotaxis/linear_ramp.jl"
```

# Linear concentration ramp

In this example we will setup an in-silico version of a typical laboratory assay,
with chemotactic bacteria moving in a linear concentration ramp, i.e. a concentration
field of the form
```math
C(x) = C_0 + (C_1 - C_0)\dfrac{x}{L_x}, \qquad x \in [0,L_x].
```

We will create a closed two-dimensional domain, mimicking a thin microfluidic chamber,
with a length of 3 mm along the `x` direction, and 1.5 mm along the `y` direction.
We will then setup the concentration field along the `x` direction and observe
chemotactic microbes as they drift towards the high-concentration region of the chamber.

The first thing we have to do is define two functions for the `concentration_field` and the
`concentration_gradient`. They must take as arguments the position of a single microbe,
and the model (from which we can access other properties of the system).
Of course, for our convenience we can dispatch these functions on whatever arguments we want,
as long as they have a method whose signature matches the MicrobeAgents API.

Importantly, the `concentration_field` must return a scalar, non-negative value.
Since the gradient is a vector quantity, the `concentration_gradient` should instead return
an iterable with length equal to the system dimensionality; a `SVector` is the recommended
choice, but `NTuple`s, `Vector`s, etc.. work just fine.

All the parameters that we need to evaluate the concentration field and gradient, in our case
the two concentration values `C₀` and `C₁` and the chamber length `Lx`, should be extracted
from the `model`.

````@example linear_ramp
using MicrobeAgents
using Plots
@inline function concentration_field(pos, model)
C₀ = model.C₀
C₁ = model.C₁
Lx = first(spacesize(model))
concentration_field(pos,Lx,C₀,C₁)
end
@inline concentration_field(pos,Lx,C₀,C₁) = C₀ + (C₁-C₀)*pos[1]/Lx
@inline function concentration_gradient(pos, model)
C₀ = model.C₀
C₁ = model.C₁
Lx = first(spacesize(model))
concentration_gradient(pos,Lx,C₀,C₁)
end
@inline concentration_gradient(pos,Lx,C₀,C₁) = SVector{length(pos)}(i==1 ? (C₁-C₀)/Lx : 0.0 for i in eachindex(pos))
````

Now as usual we define the simulation domain and the integration timestep, but we also define
a `properties` dictionary, which we pass as a keyword argument to `StandardABM`.
This dictionary will contain all the information regarding our concentration field.

Note that the `:C₀` and `:C₁` keys have been defined by us; we could have chosen different
names for them.
The `concentration_field` and `concentration_gradient` functions instead **must** be assigned
to the `:concentration_field` and `:concentration_gradient` keys respectively; this is required
by MicrobeAgents and assigning these functions to any other key will not produce the desired results.

To observe chemotaxis, we must use a microbe type for which chemotactic behavior is implemented.
If we used the base `Microbe`, no matter what field we define, we would only observe a random walk
since no chemotactic behavior is implemented.
The most classic model of chemotaxis is implemented in the `BrownBerg` type;
we will not modify its parameters here and just stick to the default values.

````@example linear_ramp
Lx, Ly = 3000, 1500 # domain size (μm)
periodic = false
space = ContinuousSpace((Lx,Ly); periodic)
Δt = 0.1 # timestep (s)
# model setup
C₀, C₁ = 0.0, 20.0 # μM
properties = Dict(
:C₀ => C₀,
:C₁ => C₁,
:concentration_field => concentration_field,
:concentration_gradient => concentration_gradient
)
model = StandardABM(BrownBerg{2}, space, Δt; properties)
n = 100 # number of microbes
for i in 1:n
add_agent!(model) # add at random positions
end
model
````

Now that the model is created, we just run it as usual, collecting the position
of the microbes at each timestep.
The visualization is slightly more involved since we want to plot
the microbe trajectories on top of the concentration field shown as a heatmap,
but there is really no difference from what we have seen in the random walk examples.

In the figure, we will see that all the microbes drift towards the right,
where the concentration of the attractant is higher.

````@example linear_ramp
T = 120 # simulation time (s)
nsteps = round(Int, T/Δt)
adata = [position]
adf, mdf = run!(model, nsteps; adata)
traj = vectorize_adf_measurement(adf, :position)
x = first.(traj)
y = last.(traj)
ts = unique(adf.step) .* Δt
lw = eachindex(ts) ./ length(ts) .* 3
xmesh = range(0,Lx,length=100)
ymesh = range(0,Ly,length=100)
xn = @view x[:,1:10]
yn = @view y[:,1:10]
c = concentration_field.(Iterators.product(xmesh,ymesh),Lx,C₀,C₁)
heatmap(xmesh, ymesh, c', cbar=false, ratio=1, axis=false, c=:bone)
plot!(xn, yn, lab=false, lw=lw, lc=(1:n)')
scatter!(xn[end,:], yn[end,:], lab=false, m=:c, mc=1:n, msw=0.5, ms=8)
````

Loading

0 comments on commit 0e3745c

Please sign in to comment.