predef eeg multichannel

Project.toml

@@ -1,10 +1,11 @@
 name = "UnfoldSim"
 uuid = "ed8ae6d2-84d3-44c6-ab46-0baf21700804"
 authors = ["Benedikt Ehinger", "Luis Lips", "Judith Schepers", "Maanik Marathe"]
-version = "0.3.2"
+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"
@@ -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"
 Statistics = "1"
 StatsModels = "0.6,0.7"
 ToeplitzMatrices = "0.7, 0.8"

README.md

@@ -76,7 +76,7 @@ data, events = simulate(
     PinkNoise(),
 );    
 ```
-All components (design, components, onsets, noise) can be easily modified and you can simply plugin your own!
+All simulation ingredients (design, components, onsets, noise) can be easily modified and you can simply plugin your own!
 
 ## Contributions
 Contributions of any kind are very welcome. Please have a look at [CONTRIBUTING.md](https://github.com/unfoldtoolbox/UnfoldSim.jl/blob/main/CONTRIBUTING.md) for guidance on contributing to UnfoldSim.jl.

docs/literate/HowTo/componentfunction.jl

@@ -0,0 +1,59 @@
+# # Component Functions
+# HowTo put arbitrary functions into components
+
+using UnfoldSim
+using Unfold
+using Random
+using DSP
+using CairoMakie, UnfoldMakie
+
+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, modulating a hanning windows by the experimental design's duration.
+# !!! important
+#     because any function depending on `design` can be used, two things have to be taken care of:
+#     
+#     1. in case a random component exist in the function, specify a `<:AbstractRNG` within the function call , the basis might be evaluated multiple times inside `simulate`
+#     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);
+
+
+# Finally, let's plot it, sorted by duration
+
+f = Figure()
+df = UnfoldMakie.eeg_array_to_dataframe(erp')
+df.duration = repeat(generate_events(design).duration, inner = size(erp, 1))
+plot_erp!(
+    f[1, 1],
+    df,
+    mapping = (; group = :duration, color = :duration),
+    categorical_color = false,
+    categorical_group = true,
+    layout = (; legend_position = :left),
+)
+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/multichannel.jl

@@ -1,10 +1,22 @@
+# # Generate multi channel data
 
+# Here you will learn how to simulate EEG data for multiple channels/electrodes.
+# The idea is to specify a signal on source level and then use a head model or a manual projection matrix to project the source signal to a number of electrodes.
+
+# ### Setup
+# ```@raw html
+# <details>
+# <summary>Click to expand</summary>
+# ```
+## Load required packages
 using UnfoldSim
 using UnfoldMakie
 using CairoMakie
 using DataFrames
 using Random
-
+# ```@raw html
+# </details >
+# ```
 
 # ## Specifying a design
 
@@ -17,7 +29,7 @@ c2 = LinearModelComponent(; basis = p300(), formula = @formula(0 ~ 1), β = [1])
 
 
 # ## The multichannel component
-# next similar to the nested design above, we can nest the component in a `MultichannelComponent`. We could either provide the projection marix manually, e.g.:
+# Next, similar to the nested design above, we can nest the component in a `MultichannelComponent`. We could either provide the projection matrix manually, e.g.:
 mc = UnfoldSim.MultichannelComponent(c, [1, 2, -1, 3, 5, 2.3, 1])
 
 # or maybe more convenient: use the pair-syntax: Headmodel=>Label which makes use of a headmodel (HaRTmuT is currently easily available in UnfoldSim)
@@ -26,7 +38,7 @@ mc = UnfoldSim.MultichannelComponent(c, hart => "Left Postcentral Gyrus")
 mc2 = UnfoldSim.MultichannelComponent(c2, hart => "Right Occipital Pole")
 
 # !!! hint
-#     You could also specify a noise-specific component which is applied prior to projection & summing with other components
+#     You could also specify a noise-specific component which is applied prior to projection & summing with other components.
 # 
 # finally we need to define the onsets of the signal
 onset = UniformOnset(; width = 20, offset = 4);

docs/literate/HowTo/newComponent.jl

@@ -1,13 +1,26 @@
-# # New component: Duration + Shift
+# # 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
+# 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
 using CairoMakie, UnfoldMakie
 
 sfreq = 100;
+# ```@raw html
+# </details >
+# ```
 
 # ## Design
 # Let's generate a design with two columns, shift + duration
@@ -19,7 +32,8 @@ design = UnfoldSim.SingleSubjectDesign(;
 )
 
 
-# We also need a new AbstractComponent
+# ## Implement a new AbstractComponent
+# We also need a new `AbstractComponent`
 struct TimeVaryingComponent <: AbstractComponent
     basisfunction::Any
     maxlength::Any
@@ -29,7 +43,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
@@ -51,8 +65,8 @@ function basis_shiftduration(evts, maxlength)
     end
 end
 
-
-erp = UnfoldSim.simulate(
+# ## Simulate data with the new component type
+erp = UnfoldSim.simulate_component(
     MersenneTwister(1),
     TimeVaryingComponent(basis_shiftduration, 50),
     design,

docs/literate/HowTo/newDesign.jl

@@ -1,13 +1,24 @@
+# # Define a new (imbalanced) design
+
+# 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` function.
+
+# ### Setup
+# ```@raw html
+# <details>
+# <summary>Click to expand</summary>
+# ```
 using UnfoldSim
 using StableRNGs
 using DataFrames
 using Parameters
-# ## Define a new Design
-# 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` function
-#
+# ```@raw html
+# </details>
+# <br />
+# ```
+
 # #### 1) `type`
 # We need a `ImbalanceSubjectDesign` struct. You are free to implement it as you wish, as long as the other two functions are implemented
 #

docs/literate/HowTo/predefinedData.jl

@@ -2,10 +2,20 @@
 
 # Let's say you want to use the events data frame (containing the levels of the experimental variables and the event onsets (latencies)) from a previous study in your simulation.
 
+# ### Setup
+# ```@raw html
+# <details>
+# <summary>Click to expand</summary>
+# ```
+## Load required packages
 using UnfoldSim
 using DataFrames
 using Random
 using CairoMakie # for plotting
+# ```@raw html
+# </details >
+# <br />
+# ```
 
 # From a previous study, we (somehow, e.g. by using [pyMNE.jl](https://unfoldtoolbox.github.io/Unfold.jl/dev/HowTo/pymne/)) imported an event data frame like this:
 my_events = DataFrame(:condition => [:A, :B, :B, :A, :A], :latency => [7, 13, 22, 35, 41])

docs/literate/HowTo/repeatTrials.jl

@@ -1,12 +1,10 @@
-using UnfoldSim
-
-
-# ## Repeating Design entries
+# # [Get multiple trials with identical subject/item combinations](@id howto_repeat_design)
 # Sometimes we want to repeat a design, that is, have multiple trials with identical values, but it is not always straight forward to implement. 
 # For instance, there is no way to easily modify `MultiSubjectDesign` to have multiple identical subject/item combinations,
 # without doing awkward repetitions of condition-levels or something.
 
 # If you struggle with this problem `RepeatDesign` is an easy tool for you:
+using UnfoldSim
 
 designOnce = MultiSubjectDesign(;
     n_items = 2,

docs/literate/HowTo/sequence.jl

@@ -0,0 +1,77 @@
+using Base: add_sum
+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, "SR_", StableRNG(1))
+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 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(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
+lines(data)
+vlines!(evts.latency[evts.event.=='S'], color = (:darkblue, 0.5))
+vlines!(evts.latency[evts.event.=='R'], color = (:darkred, 0.5))
+xlims!(0, 500)
+current_figure()
+
+# As visible, the `R` response always follows the `S` response. Due to the "`_`" we have large breaks between the individual sequences.

docs/literate/reference/basistypes.jl

@@ -1,16 +1,25 @@
-using UnfoldSim
-using CairoMakie
-using DSP
-using StableRNGs
-
-# ## Basistypes
+# # Overview: Basis function (component) types
 # There are several basis types directly implemented. They can be easily used for the `components`.
 #
 # !!! note
 #       You can use any arbitrary shape defined by yourself! We often make use of `hanning(50)` from the DSP.jl package.
 
+# ### Setup
+# ```@raw html
+# <details>
+# <summary>Click to expand</summary>
+# ```
+## Load required packages
+using UnfoldSim
+using CairoMakie
+using DSP
+using StableRNGs
+# ```@raw html
+# </details >
+# ```
+
 # ## EEG
-# By default, the EEG bases assume a sampling rate of 100, which can easily be changed by e.g. p100(;sfreq=300)
+# By default, the EEG bases assume a sampling rate of 100, which can easily be changed by e.g. p100(; sfreq=300)
 f = Figure()
 ax = f[1, 1] = Axis(f)
 for b in [p100, n170, p300, n400]