OnsetFormulas, Sequences, and Component Functions - v0.4.0 -> actually v0.5.0

.github/workflows/CI.yml

@@ -5,8 +5,6 @@ on:
       - main
     tags: '*'
   pull_request:
-    branches:
-      - main
 concurrency:
   # Skip intermediate builds: always.
   # Cancel intermediate builds: only if it is a pull request build.

Project.toml

@@ -5,6 +5,7 @@ version = "0.4.0"
 
 [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"
@@ -23,7 +24,7 @@ ToeplitzMatrices = "c751599d-da0a-543b-9d20-d0a503d91d24"
 
 [compat]
 Artifacts = "1"
-DSP = "0.7"
+DSP = "0.7,0.8"
 DataFrames = "1"
 Distributions = "0.25"
 FileIO = "1"
@@ -34,7 +35,7 @@ MixedModels = "4"
 MixedModelsSim = "0.2"
 Parameters = "0.12"
 Random = "1"
-SignalAnalysis = "0.4, 0.5,0.6"
+SignalAnalysis = "0.4, 0.5,0.6,0.7,0.8,0.9,0.10"
 Statistics = "1"
 StatsModels = "0.6,0.7"
 ToeplitzMatrices = "0.7, 0.8"

docs/Project.toml

@@ -15,7 +15,5 @@ StatsBase = "2913bbd2-ae8a-5f71-8c99-4fb6c76f3a91"
 TimerOutputs = "a759f4b9-e2f1-59dc-863e-4aeb61b1ea8f"
 Unfold = "181c99d8-e21b-4ff3-b70b-c233eddec679"
 UnfoldMakie = "69a5ce3b-64fb-4f22-ae69-36dd4416af2a"
+UnfoldMixedModels = "019ae9e0-8363-565c-86e5-97a5a2fe84f4"
 UnfoldSim = "ed8ae6d2-84d3-44c6-ab46-0baf21700804"
-
-[compat]
-Unfold = "0.7"

docs/literate/HowTo/componentfunction.jl

@@ -0,0 +1,69 @@
+# # Component Basisfunctions
+# HowTo use functions that depend on the `design` and return  per-event basis-vectors, instead of the same basis vector for all events.
+
+
+# ### Setup
+# ```@raw html
+# <details>
+# <summary>Click to expand</summary>
+# ```
+## Load required packages
+using UnfoldSim
+using Unfold
+using Random
+using DSP
+using CairoMakie, UnfoldMakie
+# ```@raw html
+# </details >
+# ```
+
+
+sfreq = 100;
+
+# ## Design
+# Let's generate a design with a categorical effect and a continuous duration effect
+design = UnfoldSim.SingleSubjectDesign(;
+    conditions = Dict(
+        :category => ["dog", "cat"],
+        :duration => Int.(round.(20 .+ rand(100) .* sfreq)),
+    ),
+);
+
+
+# Instead of defining a "boring" vector basis function e.g. `[0,0,1,2,3,3,2,1,0,0,0]`, let's use function - in our case a hanning window with the size depending on the experimental design's duration.
+# !!! important
+#     Two things have to be taken care of:
+#     1. in case a rng is required to e.g. generate the design, or your basisfunction depends on it, you have to specify a two-argument basis-function: `(rng,design)->...`
+#     2. a `maxlength` has to be specified via a tuple `(function,maxlength)``
+
+mybasisfun = design -> hanning.(generate_events(design).duration)
+signal = LinearModelComponent(;
+    basis = (mybasisfun, 100),
+    formula = @formula(0 ~ 1 + category),
+    β = [1, 0.5],
+);
+
+erp = UnfoldSim.simulate_component(MersenneTwister(1), signal, design);
+
+# After simulation, we are ready to plot it. We expect that the simulated responses are scaled by the design's duration. To show it more effectively, we sort by duration.
+##---
+f = Figure()
+df = UnfoldMakie.eeg_array_to_dataframe(erp')
+df.duration = repeat(generate_events(design).duration, inner = size(erp, 1))
+df.category = repeat(generate_events(design).category, inner = size(erp, 1))
+plot_erp!(
+    f[1, 1],
+    df,
+    mapping = (; group = :group => nonnumeric, col = :category), #  color = :duration, fails right nowUnfoldMakie#353
+    layout = (; legend_position = :left),
+    colorbar = (; label = "Duration"),
+)
+plot_erpimage!(
+    f[2, 1],
+    erp,
+    sortvalues = generate_events(design).duration,
+    layout = (; legend_position = :bottom),
+)
+f
+
+# The scaling by the two `condition` effect levels and the modified event duration by the `duration` are clearly visible

docs/literate/HowTo/getGroundTruth.jl

@@ -0,0 +1,116 @@
+# # Simulate ground truth marginalized Effects
+
+# Often when testing some algorithm, we want to compare our results to a known ground truth. In the case of marginalized effects via the `Unfold.effects`/ `Effects.jl` interface, we can do this using an `EffectsDesign`.
+# You can find more on what marginalized effects are here in the [Unfold.jl documentation](https://unfoldtoolbox.github.io/Unfold.jl/dev/generated/HowTo/effects/)
+
+# ### Setup
+# ```@raw html
+# <details>
+# <summary>Click to expand</summary>
+# ```
+## Load required packages
+using UnfoldSim
+using Unfold
+using CairoMakie
+using UnfoldMakie
+using Random
+# ```@raw html
+# </details >
+# ```
+# ## Simulation
+# First let's make up a SingleSubject simulation
+
+# !!! note
+#     Getting a ground truth for a MultiSubjectDesign is not implemented yet
+
+design =
+    SingleSubjectDesign(;
+        conditions = Dict(
+            :condition => ["bike", "face"],
+            :continuous => range(0, 5, length = 10),
+        ),
+    ) |> x -> RepeatDesign(x, 5);
+
+# **n170** has a condition effect, faces are more negative than bikes
+n1 = LinearModelComponent(;
+    basis = n170(),
+    formula = @formula(0 ~ 1 + condition),
+    β = [5, 3],
+);
+# **p300** has a continuous effect, higher continuous values will result in larger P300's.
+# We include both a linear and a quadratic effect of the continuous variable.
+p3 = LinearModelComponent(;
+    basis = p300(),
+    formula = @formula(0 ~ 1 + continuous + continuous^2),
+    β = [5, 1, 0.2],
+);
+
+components = [n1, p3]
+data, evts = simulate(
+    MersenneTwister(1),
+    design,
+    components,
+    UniformOnset(; width = 0, offset = 1000),
+    PinkNoise(),
+);
+
+# ## Simulate marginalized effects directly
+# To marginalize effects we first have to specify an effects dictionary and subsequently hand this dict plus the original design to `EffectsDesign()`
+
+effects_dict = Dict(:condition => ["bike", "face"])
+
+effects_design = EffectsDesign(design, effects_dict)
+
+# !!! note
+#     We only specified the condition levels here, by default every unspecified variable will be set to a "typical" (i.e. the mean) value.
+
+# And finally we can simulate our ground truth marginal effects
+
+gt_data, gt_events = simulate(
+    MersenneTwister(1),
+    effects_design,
+    components,
+    NoOnset(),
+    NoNoise(),
+    return_epoched = true,
+);
+@show gt_events
+
+# Additionally, we can get the simulated effects into a tidy dataframe using Unfold's `result_to_table`.
+# Note that the data has to be reshaped into a channel X times X predictor form. (In our one channel example `size(gt_data) = (45,2)`, missing the channel dimension)
+
+g = reshape(gt_data, 1, size(gt_data)...)
+times = range(1, 45);
+gt_effects = Unfold.result_to_table([g], [gt_events], [times], ["effects"])
+first(gt_effects, 5)
+
+
+# ## Compare with Unfold.jl results
+
+m = fit(
+    UnfoldModel,
+    [
+        Any => (
+            @formula(0 ~ 1 + condition + spl(continuous, 4)),
+            firbasis(τ = [-0.1, 1], sfreq = 100, name = "basis"),
+        ),
+    ],
+    evts,
+    data,
+);
+
+ef = effects(effects_dict, m);
+
+# !!! note
+#      The ground truth is shorter because the ground truth typically returns values between `[0 maxlength(components)]`, whereas in our unfold-model we included a baseline period of 0.1s.
+#      If you want to actually compare results with the ground truth, you could either us `UnfoldSim.pad_array()` or set the Unfold modelling window to `τ=[0,1]`
+
+gt_effects.type .= "UnfoldSim effects"
+ef.type .= "Unfold effects"
+
+gt_effects.time = gt_effects.time ./ 100 .- 1 / 100
+ef.continuous .= 2.5 # needed to be able to easily merge the two dataframes
+comb = vcat(gt_effects, ef)
+plot_erp(comb; mapping = (; color = :type, col = :condition))
+
+# The simulated ground truth marginal effects, and the fitted marginal effects look similar as expected, but the fitted has some additional noise because of finite data (also as expected).

docs/literate/HowTo/newComponent.jl

@@ -1,13 +1,18 @@
 # # Define a new component (with variable duration and shift)
 
 # We want a new component that changes its duration and shift depending on a column in the event design. This is somewhat already implemented in the HRF + Pupil bases.
+# !!! hint
+#     if you are just interested to use duration-dependency in your simulation, check out the component-function tutorial
+
+
 
 # ### Setup
 # ```@raw html
 # <details>
 # <summary>Click to expand</summary>
 # ```
 using UnfoldSim
+import UnfoldSim.simulate_component
 using Unfold
 using Random
 using DSP
@@ -39,7 +44,7 @@ end
 Base.length(c::TimeVaryingComponent) = length(c.maxlength)
 
 # While we could have put the TimeVaryingComponent.basisfunction directly into the simulate function, I thought this is a bit more modular
-function UnfoldSim.simulate(rng, c::TimeVaryingComponent, design::AbstractDesign)
+function UnfoldSim.simulate_component(rng, c::TimeVaryingComponent, design::AbstractDesign)
     evts = generate_events(design)
     return c.basisfunction(evts, c.maxlength)
 end
@@ -62,7 +67,7 @@ function basis_shiftduration(evts, maxlength)
 end
 
 # ## Simulate data with the new component type
-erp = UnfoldSim.simulate(
+erp = UnfoldSim.simulate_component(
     MersenneTwister(1),
     TimeVaryingComponent(basis_shiftduration, 50),
     design,

docs/literate/HowTo/newDesign.jl

@@ -3,8 +3,10 @@
 # A design specifies how much data is generated, and how the event-table(s)
 # should be generated. Already implemented examples are `MultiSubjectDesign` and `SingleSubjectDesign`.
 
+
 # We need 3 things for a new design: a `struct<:AbstractDesign`, a `size` and a `generate_events` function.
 
+
 # ### Setup
 # ```@raw html
 # <details>
@@ -15,6 +17,7 @@ using StableRNGs
 using DataFrames
 using Parameters
 using Random
+
 # ```@raw html
 # </details>
 # <br />

docs/literate/HowTo/sequence.jl

@@ -0,0 +1,91 @@
+# # Sequence of events (e.g. SCR)
+
+# In this HoWTo we learn to simulate a "SR"-Sequence, a stimulus response, followed by a button press response. 
+# ### Setup
+# ```@raw html
+# <details>
+# <summary>Click to expand</summary>
+# ```
+## Load required packages
+using UnfoldSim
+using CairoMakie
+using StableRNGs
+# ```@raw html
+# </details >
+# ```
+
+
+# 
+# First we generate the minimal design of the experiment by specifying our conditions (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`" (stimulus) event and an "`R`" (response) event - only single letter events are supported! The "`_`" is a signal for the onset generator to generate a bigger pause - no overlap between adjacent "`SR`" pairs.
+design = SequenceDesign(design, "SR_")
+generate_events(StableRNG(1), 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 advanced sequences are possible as well, like "SR{1,3}", or "A[BC]". Infinite sequences are **not** possible like "AB*". 
+
+# Finally, let's repeat the current design 4 times
+design = RepeatDesign(design, 4)
+generate_events(StableRNG(1), design)
+
+# This results in 16 trials that nicely follow our sequence
+
+# !!! hint
+#     There is a difference between `SequenceDesign(RepeatDesign)` and `RepeatDesign(SequenceDesign)` for variable sequences e.g. "A[BC]", where in the former case,  one sequence is drawn e.g. "AC" and applied to all repeated rows, in the latter, one sequence for each repeat is drawn.
+
+
+# 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, 0],
+)
+
+resp = LinearModelComponent(;
+    basis = UnfoldSim.hanning(Int(0.5 * 100)), # sfreq = 100 for the other bases
+    formula = @formula(0 ~ 1 + condition),
+    β = [1, 2],
+    offset = -10,
+)
+nothing ## hide
+
+
+# We combine them into a dictionary with a sequence-`Char` as key and simulate
+components = Dict('S' => [p1, n1, p3], 'R' => [resp])
+
+data, evts = simulate(
+    StableRNG(1),
+    design,
+    components,
+    UniformOnset(offset = 40, width = 10),
+    NoNoise(),
+)
+nothing ## hide
+
+# Finally we can plot the results
+f, ax, h = lines(data)
+vlines!(ax, evts.latency[evts.event.=='S'], color = (:darkblue, 0.5))
+vlines!(ax, evts.latency[evts.event.=='R'], color = (:darkred, 0.5))
+ax.xlabel = "Time [samples]"
+ax.ylabel = "EEG [a.u]"
+xlims!(ax, 0, 500)
+f
+
+# As visible, the `R` response always follows the `S` response. Due to the "`_`" we have large breaks between the individual sequences.
+
+
+

docs/literate/reference/designtypes.jl

@@ -40,7 +40,9 @@ design_single_shuffled = SingleSubjectDesign(;
         :stimulus_type => ["natural", "artificial"],
         :contrast_level => range(0, 1, length = 3),
     ),
+
     event_order_function = shuffle,
+
 );
 # ```@raw html
 # <details>