improved and added docstrings

src/bases.jl

@@ -1,10 +1,25 @@
 # here we define some commonly used basis for simulation
 
 ## EEG
-# P100, N170, P300, N400
+"""
+    p100(;sfreq=100)
+Generator for Hanning window, peak at 100ms, width 100ms, at kwargs `sfreq` (default 100)
+"""
 p100(; sfreq = 100) = hanning(0.1, 0.1, sfreq)
+"""
+    p300(;sfreq=100)
+Generator for Hanning window, peak at 300ms, width 300ms, at kwargs `sfreq` (default 100)
+"""
 p300(; sfreq = 100) = hanning(0.3, 0.3, sfreq)
+"""
+    n170(;sfreq=100)
+Generator for Hanning window, negative (!) peak at 170ms, width 150ms, at kwargs `sfreq` (default 100)
+"""
 n170(; sfreq = 100) = -hanning(0.15, 0.17, sfreq)
+"""
+    n400(;sfreq=100)
+Generator for Hanning window, negative (!) peak at 400ms, width 400ms, at kwargs `sfreq` (default 100)
+"""
 n400(; sfreq = 100) = -hanning(0.4, 0.4, sfreq)
 
 """
@@ -20,7 +35,10 @@ function DSP.hanning(duration, offset, sfreq)
 end
 
 ## pupil
-
+"""
+    PuRF()
+Default generator for PuRF Pupil Response Function
+"""
 function PuRF(; n = 10.1, tmax = 0.93, sfreq = 100)
     t = (0:1/sfreq:3*tmax)
     return PuRF(t, n, tmax) ./ PuRF(tmax, n, tmax)

src/component.jl

@@ -82,22 +82,31 @@ end
 Base.length(c::MultichannelComponent) = length(c.component)
 
 """
+    n_channels(c::AbstractComponent)
 Returns the number of channels. By default = 1
 """
 n_channels(c::AbstractComponent) = 1
 
 """
-for `MultichannelComponent` returns the length of the projection vector
+    n_channels(c::MultichannelComponent)
+for `MultichannelComponent` returns the length of the projection vector.
+
 """
 n_channels(c::MultichannelComponent) = length(c.projection)
 
-
+"""
+for a vector of multichanelcomponents, returns the first but asserts all are of equal length
+"""
 function n_channels(c::Vector{<:AbstractComponent})
     all_channels = n_channels.(c)
     @assert length(unique(all_channels)) == 1 "Error - projections of different channels cannot be different from eachother"
     return all_channels[1]
 end
 
+"""
+    simulate_component(rng,c::MultichannelComponent,design::AbstractDesign)
+Returns the projection of a component from source to "sensor" space
+"""
 function simulate_component(rng, c::MultichannelComponent, design::AbstractDesign)
     y = simulate_component(rng, c.component, design)
 
@@ -112,17 +121,22 @@ end
 
 
 Base.length(c::AbstractComponent) = length(c.basis)
+
+"""
+    maxlength(c::Vector{AbstractComponent}) = maximum(length.(c))
+"""
 maxlength(c::Vector{AbstractComponent}) = maximum(length.(c))
 
 """
-# by default call simulate with `::Abstractcomponent,::AbstractDesign``, but allow for custom types
-# making use of other information in simulation
+    simulate_component(rng, c::AbstractComponent, simulation::Simulation)
+by default call simulate_component with `(::Abstractcomponent,::AbstractDesign)` instead of whole simulation. This allows users to provide a hook to do something completly different :)
 """
 simulate_component(rng, c::AbstractComponent, simulation::Simulation) =
     simulate_component(rng, c, simulation.design)
 
 """
-simulate a linearModel
+    simulate_component(rng, c::AbstractComponent, simulation::Simulation)
+Generates a linear model design matrix and weights it by c.β
 
 julia> c = UnfoldSim.LinearModelComponent([0,1,1,0],@formula(0~1+cond),[1,2],Dict())
 julia> design = MultiSubjectDesign(;n_subjects=2,n_items=50,item_between=(;:cond=>["A","B"]))
@@ -147,7 +161,12 @@ function simulate_component(rng, c::LinearModelComponent, design::AbstractDesign
     return y' .* c.basis
 end
 """
-simulate MixedModelComponent
+    simulate_component(rng, c::MixedModelComponent, design::AbstractDesign)
+Generates a MixedModel and simulates data according to c.β and c.σs
+
+A trick is used to remove the Normal-Noise from the MixedModel which might lead to rare numerical instabilities. Practically, we upscale the σs by factor 10000, and provide a σ=0.0001. Internally this results in a normalization where the response scale is 10000 times larger than the noise.
+
+Currently it is not possible to use a different basis for fixed and random effects, but a code-stub exists (it is slow though)
 
 julia> design = MultiSubjectDesign(;n_subjects=2,n_items=50,item_between=(;:cond=>["A","B"]))
 julia> c = UnfoldSim.MixedModelComponent([0.,1,1,0],@formula(0~1+cond+(1|subject)),[1,2],Dict(:subject=>[2],),Dict())
@@ -255,6 +274,10 @@ end
 #----
 
 """
+    simulate_responses(
+        rng,
+        components::Vector{<:AbstractComponent},
+        simulation::Simulation)
 Simulates multiple component responses and accumulates them on a per-event basis
 """
 function simulate_responses(
@@ -276,6 +299,11 @@ function simulate_responses(
 end
 
 
+"""
+    simulate_and_add!(epoch_data::AbstractMatrix, c, simulation, rng)
+    simulate_and_add!(epoch_data::AbstractArray, c, simulation, rng)
+Helper function to call `simulate_component` and add it to a provided Array`
+"""
 function simulate_and_add!(epoch_data::AbstractMatrix, c, simulation, rng)
     @debug "matrix"
     @views epoch_data[1:length(c), :] .+= simulate_component(rng, c, simulation)

src/design.jl

@@ -77,6 +77,7 @@ function generate_events(design::SingleSubjectDesign)
 end
 
 """
+    generate_events(design::MultiSubjectDesign)
 Generates full factorial Dataframe according to MixedModelsSim.jl 's simdat_crossed function
 Note: n_items = you can think of it as `trials` or better, as stimuli
 
@@ -127,6 +128,7 @@ length(design::AbstractDesign) = *(size(design)...)
 # ----
 
 """
+    RepeatDesign{T}
 repeat a design DataFrame multiple times to mimick repeatedly recorded trials
 
 ```julia
@@ -145,6 +147,11 @@ design = RepeatDesign(designOnce,4);
     repeat::Int = 1
 end
 
+"""
+    UnfoldSim.generate_events(design::RepeatDesign{T})
+
+In a repeated design, iteratively calls the underlying {T} Design and concatenates. In case of MultiSubjectDesign, sorts by subject
+"""
 function UnfoldSim.generate_events(design::RepeatDesign)
     df = map(x -> generate_events(design.design), 1:design.repeat) |> x -> vcat(x...)
     if isa(design.design, MultiSubjectDesign)

src/headmodel.jl

@@ -86,14 +86,18 @@ Fallback: along the third dimension using `norm` - the maximal projection
 magnitude(headmodel::AbstractHeadmodel) = magnitude(leadfield(headmodel))
 
 """
+magnitude(headmodel::Hartmut; type = "perpendicular") =
 Extract magnitude of 3-orientation-leadfield, 
 `type` (default: "perpendicular") => uses the provided source-point orientations - otherwise falls back to `norm`
 """
 magnitude(headmodel::Hartmut; type = "perpendicular") =
     type == "perpendicular" ? magnitude(leadfield(headmodel), orientation(headmodel)) :
     magnitude(leadfield(headmodel))
 
-
+"""
+    magnitude(lf::AbstractArray{T,3}, orientation::AbstractArray{T,2}) where {T<:Real}
+Returns the magnitude along an orientation of the leadfield
+"""
 function magnitude(lf::AbstractArray{T,3}, orientation::AbstractArray{T,2}) where {T<:Real}
     si = size(lf)
     magnitude = fill(NaN, si[1:2])
@@ -105,7 +109,10 @@ function magnitude(lf::AbstractArray{T,3}, orientation::AbstractArray{T,2}) wher
     return magnitude
 end
 
-
+"""
+    magnitude(lf::AbstractArray{T,3}) where {T<:Real}
+If orientation is not specified, returns the maximal magnitude (norm of leadfield)
+"""
 function magnitude(lf::AbstractArray{T,3}) where {T<:Real}
     si = size(lf)
     magnitude = fill(NaN, si[1:2])

src/helper.jl

@@ -91,9 +91,24 @@ function closest_src(head::Hartmut, label::String)
 end
 
 
-# Adapted from Unfold.jl: https://github.com/unfoldtoolbox/Unfold.jl/blob/b3a21c2bb7e93d2f45ec64b0197f4663a6d7939a/src/utilities.jl#L40
+
 
 # One channel case
+
+"""
+    epoch(data::AbstractVector, args...; kwargs...)
+    epoch(
+        data::AbstractArray{T,2},
+        events,
+        τ::Tuple{Number,Number},
+        sfreq;
+        eventtime::Symbol = :latency,
+    ) where {T<:Union{Missing,Number}}
+Helper function to epoch data
+
+Adapted from Unfold.jl: https://github.com/unfoldtoolbox/Unfold.jl/blob/b3a21c2bb7e93d2f45ec64b0197f4663a6d7939a/src/utilities.jl#L40
+
+"""
 function epoch(data::AbstractVector, args...; kwargs...)
     data_r = reshape(data, (1, :))
     ep = epoch(data_r, args...; kwargs...)

src/noise.jl

@@ -2,38 +2,68 @@
 # Types
 #-------------
 
+
+"""
+    PinkNoise <: AbstractNoise
+Generates Pink Noise using the SignalAnalysis.jl implementation
+"""
 @with_kw struct PinkNoise <: AbstractNoise
     noiselevel = 1
     func = SignalAnalysis.PinkGaussian
 end
 
 
+"""
+    RedNoise <: AbstractNoise
+Generates Red Noise using the SignalAnalysis.jl implementation
+"""
 @with_kw struct RedNoise <: AbstractNoise
     noiselevel = 1
     func = SignalAnalysis.RedGaussian
 end
 
 
+"""
+    WhiteNoise <: WhiteNoise
+    noiselevel = 1
+    imfilter = 0
+Generates White Noise using `randn` - thus Gaussian noise.
+Using imfilter > 0 it is possible to smooth the noise using Image.imfilter
+"""
 @with_kw struct WhiteNoise <: AbstractNoise
     noiselevel = 1
     imfilter = 0
 end
 
-
+"""
+    RealisticNoise <: AbstractNoise
+Not implemented- planned to use Artefacts.jl to provide real EEG data to add
+"""
 @with_kw struct RealisticNoise <: AbstractNoise
     noiselevel = 1
 end
 
+"""
+    NoNoise <: AbstractNoise
+
+Returns zeros instead of noise
+"""
 struct NoNoise <: AbstractNoise end
 
+"""
+    AutoRegressiveNoise <: AbstractNoise
+Not implemented
+"""
 struct AutoRegressiveNoise <: AbstractNoise end
 
 
 
 """ 
-    Noise with exponential decay in AR spectrum
-    !!! warning
-        Current implementation: cholesky of NxN matrix needs to be calculated, might need lot's of RAM
+    ExponentialNoise <: AbstractNoise
+
+Noise with exponential decay in AR spectrum
+!!! warning
+    Current implementation: cholesky of NxN matrix needs to be calculated, might need lot's of RAM
 
 """
 
@@ -44,23 +74,23 @@ end
 
 
 """
-    simulate_noise(t::Union{PinkNoise, RedNoise}, n::Int)
+    simulate_noise(rng, t::Union{PinkNoise,RedNoise}, n::Int)
 
-Generate noise of a given type t and length n
+Generate Pink or Red Noise using the SignalAnalysis.jl implementation
 """
 function simulate_noise(rng, t::Union{PinkNoise,RedNoise}, n::Int)
     return t.noiselevel .* rand(rng, t.func(n, 1.0))
 end
 
+"""
+    simulate_noise(rng, t::NoNoise, n::Int)
+Returns zeros instead of noise
+"""
 function simulate_noise(rng, t::NoNoise, n::Int)
     return zeros(n)
 end
 
-"""
-    simulate_noise(t::WhiteNoise, n::Int)
 
-Generate noise of a given type t and length n
-"""
 function simulate_noise(rng, t::WhiteNoise, n::Int)
     noisevector = randn(rng, n)
     if !isnothing(t.imfilter)
@@ -70,11 +100,7 @@ function simulate_noise(rng, t::WhiteNoise, n::Int)
 end
 
 
-"""
-    simulate_noise(t::RealisticNoise, n::Int)
 
-Generate noise of a given type t and length n
-"""
 function simulate_noise(rng, t::RealisticNoise, n::Int)
     error("not implemented")
     return 0
@@ -100,8 +126,9 @@ end
 
 
 """
-Generate and add noise to the data-matrix
+    add_noise!(rng, noisetype::AbstractNoise, signal)
 
+Generate and add noise to a data matrix.
 Assumes that the signal can be linearized, that is, that the noise is stationary
 """
 function add_noise!(rng, noisetype::AbstractNoise, signal)