sequence

Project.toml

@@ -5,6 +5,7 @@ version = "0.3.1"
 
 [deps]
 Artifacts = "56f22d72-fd6d-98f1-02f0-08ddc0907c33"
+Automa = "67c07d97-cdcb-5c2c-af73-a7f9c32a568b"
 DSP = "717857b8-e6f2-59f4-9121-6e50c889abd2"
 DataFrames = "a93c6f00-e57d-5684-b7b6-d8193f3e46c0"
 Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f"

docs/literate/HowTo/sequence.jl

@@ -0,0 +1,52 @@
+using UnfoldSim
+using CairoMakie
+using StableRNGs
+
+# ## Stimulus - Response design
+
+# let's say we want to simulate a stimulus response, followed by a button press response. 
+# First we generate the minimal design of the experiment by specifying our conditins (a one-condition-two-levels design in our case)
+design = SingleSubjectDesign(conditions = Dict(:condition => ["one", "two"]))
+generate_events(design)
+# next we use the `SequenceDesign` and nest our initial design in it. "SR_" is code for an "S" event and an "R" event - only single letter events are supported! The `_` is a signal for the Onset generator to generate a bigger pause - no overlap between adjacend `SR` pairs
+design = SequenceDesign(design, "SRR_")
+generate_events(design)
+# The main thing that happened is that the design was repeated for every event (each 'letter') of the sequence, and an `eventtype` column was added.
+# !!! hint
+#     more advaned sequences exist, like "SR{1,3}", or "A[BC]" etc.
+
+# Finally, let's repeat the design 2 times
+design = RepeatDesign(design, 2)
+generate_events(design)
+
+# This results in 20 trials that nicely follow our sequence
+
+# Next we have to specify for both events `S` and `R` what the responses should look like.
+
+p1 = LinearModelComponent(;
+    basis = p100(),
+    formula = @formula(0 ~ 1 + condition),
+    β = [1, 0.5],
+);
+
+n1 = LinearModelComponent(;
+    basis = n170(),
+    formula = @formula(0 ~ 1 + condition),
+    β = [1, 0.5],
+);
+p3 = LinearModelComponent(;
+    basis = UnfoldSim.hanning(Int(0.5 * 100)), # sfreq = 100 for the other bases
+    formula = @formula(0 ~ 1 + condition),
+    β = [1, 2],
+);
+
+components = Dict('S' => [p1, n1], 'R' => [p3])
+
+data, evts = simulate(
+    StableRNG(1),
+    design,
+    components,
+    UniformOnset(offset = 10, width = 100),
+    NoNoise(),
+)
+lines(data)

docs/literate/reference/onsettypes.jl

@@ -282,3 +282,25 @@ end
 #        - if `offset` < `length(signal.basis)` -> there might be overlap, depending on the other parameters of the onset distribution
 
 # [^1]: Wikipedia contributors. (2023, December 5). Log-normal distribution. In Wikipedia, The Free Encyclopedia. Retrieved 12:27, December 7, 2023, from https://en.wikipedia.org/w/index.php?title=Log-normal_distribution&oldid=1188400077# 
+
+
+
+
+# ## Design-dependent `FormulaXOnset`
+
+# For additional control we provide `FormulaUniformOnset` and `FormulaLogNormalOnset` types, that allow to control all parameters by specifying formulas
+o = UnfoldSim.FormulaUniformOnset(
+    width_formula = @formula(0 ~ 1 + cond),
+    width_β = [50, 20],
+)
+events = generate_events(design)
+onsets = UnfoldSim.simulate_interonset_distances(MersenneTwister(42), o, design)
+
+f = Figure()
+ax = f[1, 1] = Axis(f)
+hist!(ax, onsets[events.cond.=="A"], bins = range(0, 100, step = 1), label = "cond: A")
+hist!(ax, onsets[events.cond.=="B"], bins = range(0, 100, step = 1), label = "cond: B")
+axislegend(ax)
+f
+
+# Voila - the inter-onset intervals are `20` samples longer for condition `B`, exactly as specified.`

docs/make.jl

@@ -48,6 +48,7 @@ makedocs(;
             "Define a new duration & jitter component" => "./generated/HowTo/newComponent.md",
             "Generate multi channel data" => "./generated/HowTo/multichannel.md",
             "Use predefined design / onsets data" => "./generated/HowTo/predefinedData.md",
+            "Produce specific sequences of events" => "./generated/HowTo/sequence.md",
         ],
         "DocStrings" => "api.md",
     ],

src/UnfoldSim.jl

@@ -14,6 +14,7 @@ using LinearAlgebra
 using ToeplitzMatrices # for AR Expo. Noise "Circulant"
 using StatsModels
 using HDF5, Artifacts, FileIO
+using Automa # for sequence
 
 using LinearAlgebra # headmodel
 
@@ -30,6 +31,7 @@ include("onset.jl")
 include("predefinedSimulations.jl")
 include("headmodel.jl")
 include("helper.jl")
+include("sequence.jl")
 include("bases.jl")
 
 export size, length
@@ -46,7 +48,7 @@ export Simulation
 export MixedModelComponent, LinearModelComponent
 
 # export designs
-export MultiSubjectDesign, SingleSubjectDesign, RepeatDesign
+export MultiSubjectDesign, SingleSubjectDesign, RepeatDesign, SequenceDesign
 
 # noise functions
 export PinkNoise, RedNoise, WhiteNoise, NoNoise, ExponentialNoise #,RealNoise (not implemented yet)

src/component.jl

@@ -99,10 +99,16 @@
 """
 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"
+    @assert length(unique(all_channels)) == 1 "Error - projections of different components have to be of the same output (=> channel) dimension"
     return all_channels[1]
 end
 
+function n_channels(components::Dict)
+    all_channels = [n_channels(c) for c in values(components)]
+    @assert length(unique(all_channels)) == 1 "Error - projections of different components have to be of the same output (=> channel) dimension"
+    return all_channels[1]
+
+end
 """
     simulate_component(rng,c::MultichannelComponent,design::AbstractDesign)
 Return the projection of a component from source to "sensor" space.
@@ -124,10 +130,13 @@
 
 """
     maxlength(c::Vector{AbstractComponent}) = maximum(length.(c))
+    maxlength(components::Dict) 
 maximum of individual component lengths
 """
-maxlength(c::Vector{AbstractComponent}) = maximum(length.(c))
+maxlength(c::Vector{<:AbstractComponent}) = maximum(length.(c))
 
+
+maxlength(components::Dict) = maximum([maximum(length.(c)) for c in values(components)])
 """
     simulate_component(rng, c::AbstractComponent, simulation::Simulation)
 By default call `simulate_component` with `(::Abstractcomponent,::AbstractDesign)` instead of the whole simulation. This allows users to provide a hook to do something completely different :)
@@ -145,22 +154,32 @@
 """
 function simulate_component(rng, c::LinearModelComponent, design::AbstractDesign)
     events = generate_events(design)
+    X = generate_designmatrix(c.formula, events, c.contrasts)
+    y = X * c.β
+
+    return y' .* c.basis
+end
+
 
+"""
+Helper function to generate a designmatrix from formula, events and contrasts.
+"""
+function generate_designmatrix(formula, events, contrasts)
     # special case, intercept only 
     # https://github.com/JuliaStats/StatsModels.jl/issues/269
-    if c.formula.rhs == ConstantTerm(1)
+    if formula.rhs == ConstantTerm(1)
         X = ones(nrow(events), 1)
     else
-        if isempty(c.contrasts)
-            m = StatsModels.ModelFrame(c.formula, events)
+        if isempty(contrasts)
+            m = StatsModels.ModelFrame(formula, events)
         else
-            m = StatsModels.ModelFrame(c.formula, events; contrasts = c.contrasts)
+            m = StatsModels.ModelFrame(formula, events; contrasts = contrasts)
         end
         X = StatsModels.modelmatrix(m)
     end
-    y = X * c.β
-    return y' .* c.basis
+    return X
 end
+
 """
     simulate_component(rng, c::MixedModelComponent, design::AbstractDesign)
 Generates a MixedModel and simulates data according to c.β and c.σs.
@@ -286,19 +305,55 @@
     components::Vector{<:AbstractComponent},
     simulation::Simulation,
 )
-    if n_channels(components) > 1
-        epoch_data =
-            zeros(n_channels(components), maxlength(components), length(simulation.design))
-    else
-        epoch_data = zeros(maxlength(components), length(simulation.design))
-    end
+    epoch_data = init_epoch_data(components, simulation.design)
+    simulate_responses!(rng, epoch_data, components, simulation)
+    return epoch_data
+end
 
+function simulate_responses!(
+    rng,
+    epoch_data::AbstractArray,
+    components::Vector,
+    simulation::Simulation,
+)
     for c in components
         simulate_and_add!(epoch_data, c, simulation, rng)
     end
     return epoch_data
 end
+function init_epoch_data(components, design)
+    if n_channels(components) > 1
+        epoch_data = zeros(n_channels(components), maxlength(components), length(design))
+    else
+        epoch_data = zeros(maxlength(components), length(design))
+    end
+    return epoch_data
+end
 
+function simulate_responses(rng, components::Dict, s::Simulation)
+    epoch_data = init_epoch_data(components, s.design)
+    evts = generate_events(s.design)
+    multichannel = n_channels(components) > 1
+    for key in keys(components)
+        if key == '_'
+            continue
+        end
+        s_key = Simulation(
+            s.design |> x -> SubselectDesign(x, key),
+            components[key],
+            s.onset,
+            s.noisetype,
+        )
+        ix = evts.event .== key
+        if multichannel
+            simulate_responses!(rng, @view(epoch_data[:, :, ix]), components[key], s_key)
+        else
+            #@debug sum(ix), size(simulate_responses(rng, components[key], s_key)), key
+            simulate_responses!(rng, @view(epoch_data[:, ix]), components[key], s_key)
+        end
+    end
+    return epoch_data
+end
 
 """
     simulate_and_add!(epoch_data::AbstractMatrix, c, simulation, rng)

src/design.jl

@@ -147,6 +147,28 @@
     repeat::Int = 1
 end
 
+
+@with_kw struct SequenceDesign{T} <: AbstractDesign
+    design::T
+    sequence::String = ""
+end
+
+UnfoldSim.generate_events(design::SequenceDesign{MultiSubjectDesign}) =
+    error("not yet implemented")
+
+function UnfoldSim.generate_events(design::SequenceDesign)
+    df = generate_events(design.design)
+    nrows_df = size(df, 1)
+    currentsequence = sequencestring(design.sequence)
+    currentsequence = replace(currentsequence, "_" => "")
+    df = repeat(df, inner = length(currentsequence))
+    df.event .= repeat(collect(currentsequence), nrows_df)
+
+    return df
+
+end
+
+
 """
     UnfoldSim.generate_events(design::RepeatDesign{T})
 
@@ -160,6 +182,28 @@
     return df
 
 end
+
+
+"""
+Internal helper design to subset a sequence design in its individual components
+"""
+struct SubselectDesign{T} <: AbstractDesign
+    design::T
+    key::Char
+end
+
+function generate_events(design::SubselectDesign)
+    return subset(generate_events(design.design), :event => x -> x .== design.key)
+end
+
+
 Base.size(design::RepeatDesign{MultiSubjectDesign}) =
     size(design.design) .* (design.repeat, 1)
 Base.size(design::RepeatDesign{SingleSubjectDesign}) = size(design.design) .* design.repeat
+Base.size(design::SequenceDesign) =
+    size(design.design) .* length(replace(design.sequence, "_" => ""))
+
+Base.size(design::RepeatDesign{SequenceDesign{SingleSubjectDesign}}) =
+    size(design.design) .* design.repeat
+
+Base.size(design::SubselectDesign) = size(generate_events(design), 1)