diff --git a/Project.toml b/Project.toml
index 889f5c68..45178896 100644
--- a/Project.toml
+++ b/Project.toml
@@ -5,6 +5,7 @@ version = "0.4.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"
diff --git a/docs/Project.toml b/docs/Project.toml
index c3cb1bdf..7257e3f5 100644
--- a/docs/Project.toml
+++ b/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"
\ No newline at end of file
diff --git a/docs/literate/HowTo/componentfunction.jl b/docs/literate/HowTo/componentfunction.jl
new file mode 100644
index 00000000..be6a1d06
--- /dev/null
+++ b/docs/literate/HowTo/componentfunction.jl
@@ -0,0 +1,69 @@
+# # [Define design-dependent component basis functions](@id componentfunction)
+# Here you will learn how to specify a component basis that uses a function that depends on the `design` and returns per-event basis vectors, instead of the same basis vector for all events.
+
+
+# ### Setup
+# ```@raw html
+#
+# Click to expand
+# ```
+## Load required packages
+using UnfoldSim
+using Unfold
+using Random
+using DSP
+using CairoMakie, UnfoldMakie
+# ```@raw html
+#
+# ```
+
+
+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 a 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 basis function 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 length of the simulated responses is 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.
diff --git a/docs/literate/HowTo/getGroundTruth.jl b/docs/literate/HowTo/getGroundTruth.jl
new file mode 100644
index 00000000..e9618dbb
--- /dev/null
+++ b/docs/literate/HowTo/getGroundTruth.jl
@@ -0,0 +1,116 @@
+# # Simulate ground truth marginal effects
+
+# Often when testing some algorithm, we want to compare our results to a known ground truth. In the case of marginal 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
+#
+# Click to expand
+# ```
+## Load required packages
+using UnfoldSim
+using Unfold
+using CairoMakie
+using UnfoldMakie
+using Random
+# ```@raw html
+#
+# ```
+# ## Simulation
+# First let's set up a SingleSubject simulation
+
+# !!! note
+# An `EffectsDesign` 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 marginal effects directly
+# To simulate marginal 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, size(gt_data, 1));
+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 signal 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).
diff --git a/docs/literate/HowTo/newComponent.jl b/docs/literate/HowTo/newComponent.jl
index 91da1245..a1b39cb7 100644
--- a/docs/literate/HowTo/newComponent.jl
+++ b/docs/literate/HowTo/newComponent.jl
@@ -1,6 +1,10 @@
# # 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](@ref componentfunction).
+
+
# ### Setup
# ```@raw html
@@ -8,6 +12,7 @@
# Click to expand
# ```
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,
diff --git a/docs/literate/HowTo/newDesign.jl b/docs/literate/HowTo/newDesign.jl
index 8d8bd6de..dc896511 100644
--- a/docs/literate/HowTo/newDesign.jl
+++ b/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
#
@@ -15,6 +17,7 @@ using StableRNGs
using DataFrames
using Parameters
using Random
+
# ```@raw html
#
#
diff --git a/docs/literate/HowTo/sequence.jl b/docs/literate/HowTo/sequence.jl
new file mode 100644
index 00000000..4e86f6ac
--- /dev/null
+++ b/docs/literate/HowTo/sequence.jl
@@ -0,0 +1,93 @@
+# # Sequence of events (e.g. SCR)
+
+# In this HowTo you will learn to simulate a "SR"-Sequence, a stimulus response, followed by a button press response.
+
+# ### Setup
+# ```@raw html
+#
+# Click to expand
+# ```
+## Load required packages
+using UnfoldSim
+using CairoMakie
+using StableRNGs
+# ```@raw html
+#
+#
+# ```
+
+# ## Create sequence design
+# 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 `event` 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.
+
+# ## Specify components for sequence events
+# 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
+components = Dict('S' => [p1, n1, p3], 'R' => [resp])
+
+# ## Simulate data
+# Given the design and the components. we specify onset and noise and simulate data
+
+
+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.
diff --git a/docs/literate/reference/designtypes.jl b/docs/literate/reference/designtypes.jl
index a0be4cc5..8f5330c0 100644
--- a/docs/literate/reference/designtypes.jl
+++ b/docs/literate/reference/designtypes.jl
@@ -40,6 +40,7 @@ design_single_shuffled = SingleSubjectDesign(;
:stimulus_type => ["natural", "artificial"],
:contrast_level => range(0, 1, length = 3),
),
+
event_order_function = shuffle,
);
# ```@raw html
diff --git a/docs/literate/reference/onsettypes.jl b/docs/literate/reference/onsettypes.jl
index 6618706b..18d10f81 100644
--- a/docs/literate/reference/onsettypes.jl
+++ b/docs/literate/reference/onsettypes.jl
@@ -1,6 +1,6 @@
# # Overview: Onset types
# The onset types determine the distances between event onsets in the continuous EEG signal. The distances are sampled from a certain probability distribution.
-# Currently, there are two types of onset distributions implemented: `UniformOnset` and `LogNormalOnset`.
+# Currently, there are two types of onset distributions implemented: `UniformOnset` and `LogNormalOnset`. Both are accompanied by their `UniformOnsetFormula` and `LogNormalOnsetFormula` conterparts, which allows to modify the overlap based on the design.
# ### Setup
# ```@raw html
@@ -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 `X-OnsetFormula`
+
+# For additional control, we provide `UniformOnsetFormula` and `LogNormalOnsetFormula` types, which allow to control all distribution parameters by specifying formulas based on the design
+o = UnfoldSim.UniformOnsetFormula(
+ 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.
diff --git a/docs/literate/tutorials/multisubject.jl b/docs/literate/tutorials/multisubject.jl
index 61570240..b91d95ab 100644
--- a/docs/literate/tutorials/multisubject.jl
+++ b/docs/literate/tutorials/multisubject.jl
@@ -10,9 +10,11 @@
## Load required packages
using UnfoldSim
using Unfold
+import UnfoldMixedModels
using CairoMakie
using UnfoldMakie
using DataFrames
+using StableRNGs
# ```@raw html
#
#
@@ -63,7 +65,8 @@ signal = MixedModelComponent(;
)
# and simulate!
-data, evts = simulate(design, signal, NoOnset(), NoNoise(), return_epoched = true);
+data, evts =
+ simulate(StableRNG(1), design, signal, NoOnset(), NoNoise(), return_epoched = true);
# We get data with 50 samples (our `basis` from above), with `4` items and 20 subjects. We get items and subjects separately because we chose no-overlap (via `NoOnset`) and `return_epoched = true``.
size(data)
@@ -94,7 +97,8 @@ f
# Let's continue our tutorial and simulate overlapping signals instead.
#
# We replace the `NoOnset` with an `UniformOnset` with 20 to 70 samples between subsequent events. We further remove the `return_epoched`, because we want to have continuous data for now.
-data, evts = simulate(design, signal, UniformOnset(offset = 20, width = 50), NoNoise());
+data, evts =
+ simulate(StableRNG(1), design, signal, UniformOnset(offset = 20, width = 50), NoNoise());
size(data)
# with the first dimension being continuous data, and the latter still the subjects.
@@ -108,6 +112,7 @@ series(data', solid_color = :black)
# # Analyzing these data with Unfold.jl
# We will analyze these data using the `Unfold.jl` toolbox. While preliminary support for deconvolution (overlap correction) for mixed models is available, here we will not make use of it, but rather apply a MixedModel to each timepoint, following the Mass-univariate approach.
data, evts = simulate(
+ StableRNG(1),
design,
signal,
UniformOnset(offset = 20, width = 50),
diff --git a/docs/literate/tutorials/simulateERP.jl b/docs/literate/tutorials/simulateERP.jl
index b6b94d04..3b8b9998 100644
--- a/docs/literate/tutorials/simulateERP.jl
+++ b/docs/literate/tutorials/simulateERP.jl
@@ -11,6 +11,7 @@
# Here we will learn how to simulate a typical ERP complex with P100, N170, P300.
+
# ### Setup
# ```@raw html
#
@@ -32,6 +33,7 @@ using UnfoldMakie # For plotting
# ### 1. Experimental design
# We specify an experimental design with one subject in two experimental conditions including a continuous variable with 10 values.
# To mimic randomization in an experiment, we shuffle the trials using the `event_order_function` argument. To generate more trials we repeat the design 100 times which results in 2000 trials in total.
+
design =
SingleSubjectDesign(;
conditions = Dict(
diff --git a/docs/make.jl b/docs/make.jl
index 12cd7c56..7feb7d86 100644
--- a/docs/make.jl
+++ b/docs/make.jl
@@ -50,10 +50,13 @@ makedocs(;
],
"HowTo" => [
"Define a new (imbalanced) design" => "./generated/HowTo/newDesign.md",
+ "Define design-dependent component basis functions" => "./generated/HowTo/componentfunction.md",
"Get multiple trials with identical subject/item combinations" => "./generated/HowTo/repeatTrials.md",
"Define a new component (with variable duration and shift)" => "./generated/HowTo/newComponent.md",
"Generate multi channel data" => "./generated/HowTo/multichannel.md",
"Use existing experimental designs & onsets in the simulation" => "./generated/HowTo/predefinedData.md",
+ "Simulate ground truth marginal effects" => "./generated/HowTo/getGroundTruth.md",
+ "Sequence of events (e.g. SCR)" => "./generated/HowTo/sequence.md",
],
"Developer documentation" => "developer_docs.md",
"API / Docstrings" => "api.md",
diff --git a/src/UnfoldSim.jl b/src/UnfoldSim.jl
index 80944294..13a95257 100644
--- a/src/UnfoldSim.jl
+++ b/src/UnfoldSim.jl
@@ -14,7 +14,7 @@ using LinearAlgebra
using ToeplitzMatrices # for AR Expo. Noise "Circulant"
using StatsModels
using HDF5, Artifacts, FileIO
-
+using Automa # for sequence
using LinearAlgebra # headmodel
import DSP.hanning
@@ -30,6 +30,7 @@ include("onset.jl")
include("predefinedSimulations.jl")
include("headmodel.jl")
include("helper.jl")
+include("sequence.jl")
include("bases.jl")
export size, length
@@ -46,7 +47,7 @@ export Simulation
export MixedModelComponent, LinearModelComponent
# export designs
-export MultiSubjectDesign, SingleSubjectDesign, RepeatDesign
+export MultiSubjectDesign, SingleSubjectDesign, RepeatDesign, SequenceDesign, EffectsDesign
# noise functions
export PinkNoise, RedNoise, WhiteNoise, NoNoise, ExponentialNoise #,RealNoise (not implemented yet)
@@ -64,7 +65,8 @@ export simulate,
export pad_array
# export Offsets
-export UniformOnset, LogNormalOnset, NoOnset
+export UniformOnset,
+ LogNormalOnset, NoOnset, UniformOnsetFormula, LogNormalOnsetFormula, ShiftOnsetByOne
# re-export StatsModels
export DummyCoding, EffectsCoding
diff --git a/src/bases.jl b/src/bases.jl
index 3f3fe0a2..ff4a2054 100644
--- a/src/bases.jl
+++ b/src/bases.jl
@@ -135,8 +135,8 @@ Generate a (potentially shifted) hanning window with a certain duration.
Note: This function extends the `DSP.hanning` function using multiple dispatch.
# Arguments
-- `duration`: in s.
-- `offset`: in s, defines hanning peak i.e. shift of the hanning window.
+- `width`: in s.
+- `offset`: in s, defines the location of the hanning peak i.e. shift of the hanning window. Must be > `round(width/2)` (otherwise the left part of the curve would be cut off).
- `sfreq`: Sampling rate in Hz.
# Returns
@@ -145,7 +145,7 @@ Note: This function extends the `DSP.hanning` function using multiple dispatch.
# Examples
```julia-repl
julia> UnfoldSim.hanning(0.1, 0.3, 100)
-25-element Vector{Float64}:
+35-element Vector{Float64}:
0.0
0.0
0.0
@@ -154,14 +154,22 @@ julia> UnfoldSim.hanning(0.1, 0.3, 100)
⋮
0.9698463103929542
0.75
- 0.41317591116653485
- 0.11697777844051105
+ 0.4131759111665348
+ 0.116977778440511
0.0
```
"""
-function DSP.hanning(duration, offset, sfreq)
- signal = hanning(Int(round(duration * sfreq)))
- return pad_array(signal, -Int(round(offset * sfreq / 2)), 0)
+function DSP.hanning(width, offset, sfreq)
+ width = width * sfreq
+ offset = offset * sfreq
+ signal = hanning(Int(round(width)))
+ pad_by = Int(round(offset - (length(signal)+1) / 2))
+
+ pad_by < 0 ?
+ error(
+ "The offset has to be > round((width+1)/2). Otherwise, the left part of the curve would be cut off. To create a component which starts before the event onset, one can use the `offset` parameter of a component.",
+ ) : ""
+ return pad_array(signal, -pad_by, 0)
end
diff --git a/src/component.jl b/src/component.jl
index 9d06e660..62934ef8 100644
--- a/src/component.jl
+++ b/src/component.jl
@@ -6,7 +6,7 @@ A component that adds a hierarchical relation between parameters according to a
All fields can be named. Works best with [`MultiSubjectDesign`](@ref).
# Fields
-- `basis::Any`: an object, if accessed, provides a 'basis function', e.g. `hanning(40)::Vector`, this defines the response at a single event. It will be weighted by the model prediction. Future versions will allow for functions, as of v0.3 this is restricted to array-like objects
+- `basis::Any`: an object, if accessed, provides a 'basis function', e.g. `hanning(40)::Vector`, this defines the response at a single event. It will be weighted by the model prediction. It is also possible to provide a function that evaluates to an `Vector` of `Vectors`, with the `design` as input to the function, the outer vector has to have `nrows(design)`, one for each event. The inner vector represents the basis functions which can be of different size (a ragged array). Alternatively, one can also return a Matrix with the second dimension representing `nrows(design)`. In the case of providing a function, one has to specify the `maxlength` as well in a tuple. E.g. `basis=(myfun,40)`, which would automatically cut the output of `myfun` to 40 samples. If your design depends on `rng`, e.g. because of `event_order_function=shuffle` or some special `SequenceDesign`, then you can provide a two-arguments function `(rng,design)->...`.
- `formula::Any`: Formula-object in the style of MixedModels.jl e.g. `@formula 0 ~ 1 + cond + (1|subject)`. The left-hand side is ignored.
- `β::Vector` Vector of betas (fixed effects), must fit the formula.
- `σs::Dict` Dict of random effect variances, e.g. `Dict(:subject => [0.5, 0.4])` or to specify correlation matrix `Dict(:subject=>[0.5,0.4,I(2,2)],...)`. Technically, this will be passed to the MixedModels.jl `create_re` function, which creates the θ matrices.
@@ -37,8 +37,10 @@ See also [`LinearModelComponent`](@ref), [`MultichannelComponent`](@ref).
β::Vector
σs::Dict # Dict(:subject=>[0.5,0.4]) or to specify correlationmatrix Dict(:subject=>[0.5,0.4,I(2,2)],...)
contrasts::Dict = Dict()
+ offset::Int = 0
end
-
+MixedModelComponent(basis, formula, β, σs, contrasts) =
+ MixedModelComponent(basis, formula, β, σs, contrasts, 0)
"""
LinearModelComponent <: AbstractComponent
@@ -47,11 +49,11 @@ A multiple regression component for one subject.
All fields can be named. Works best with [`SingleSubjectDesign`](@ref).
# Fields
-- `basis::Any`: an object, if accessed, provides a 'basis function', e.g. `hanning(40)::Vector`, this defines the response at a single event. It will be weighted by the model prediction. Future versions will allow for functions, as of v0.3 this is restricted to array-like objects
+- `basis::Any`: an object, if accessed, provides a 'basis function', e.g. `hanning(40)::Vector`, this defines the response at a single event. It will be weighted by the model prediction. It is also possible to provide a function that evaluates to an `Vector` of `Vectors`, with the `design` as input to the function, the outer vector has to have `nrows(design)`, one for each event. The inner vector represents the basis functions which can be of different size (a ragged array). Alternatively, one can also return a Matrix with the second dimension representing `nrows(design)`. In the case of providing a function, one has to specify the `maxlength` as well in a tuple. E.g. `basis=(myfun,40)`, which would automatically cut the output of `myfun` to 40 samples. If your design depends on `rng`, e.g. because of `event_order_function=shuffle` or some special `SequenceDesign`, then you can provide a two-arguments function `(rng,design)->...`.
- `formula::Any`: StatsModels `formula` object, e.g. `@formula 0 ~ 1 + cond` (left-hand side must be 0).
- `β::Vector` Vector of betas/coefficients, must fit the formula.
-- `contrasts::Dict` (optional): Determines which coding scheme to use for which categorical variables. Default is empty which corresponds to dummy coding.
- For more information see .
+- `contrasts::Dict` (optional): Determines which coding scheme to use for which categorical variables. Default is empty which corresponds to dummy coding. For more information see .
+- `offset::Int = 0`: Can be used to shift the basis function in time.
# Examples
```julia-repl
@@ -75,8 +77,38 @@ See also [`MixedModelComponent`](@ref), [`MultichannelComponent`](@ref).
formula::Any # e.g. 0 ~ 1 + cond - left side must be "0"
β::Vector
contrasts::Dict = Dict()
+ offset::Int = 0
+ function LinearModelComponent(basis, formula, β, contrasts, offset)
+ @assert isa(basis, Tuple{Function,Int}) "`basis` needs to be an `::Array` or a `Tuple(function::Function,maxlength::Int)`"
+ @assert basis[2] > 0 "`maxlength` needs to be longer than 0"
+ new(basis, formula, β, contrasts, offset)
+ end
+ LinearModelComponent(basis::Array, formula, β, contrasts, offset) =
+ new(basis, formula, β, contrasts, offset)
end
+# backwards compatability after introducing the `offset` field
+LinearModelComponent(basis, formula, β, contrasts) =
+ LinearModelComponent(basis, formula, β, contrasts, 0)
+
+"""
+ get_offset(AbstractComponent)
+
+Should the `basis` be shifted? Returns c.offset for most components, if not implemented for a type, returns 0. Can be positive or negative, but has to be an Integer.
+"""
+get_offset(c::AbstractComponent)::Int = 0
+get_offset(c::LinearModelComponent)::Int = c.offset
+get_offset(c::MixedModelComponent)::Int = c.offset
+get_offset(c::Vector{<:AbstractComponent}) = get_offset.(c)
+get_offset(d::Dict{<:Char,<:Vector{<:AbstractComponent}}) =
+ Dict(k => get_offset(v) for (k, v) in d)
+
+maxoffset(c::Vector{<:AbstractComponent}) = maximum(get_offset.(c))
+maxoffset(d::Dict{<:Char,<:Vector{<:AbstractComponent}}) = maximum(maxoffset.(values(d)))
+minoffset(c::Vector{<:AbstractComponent}) = minimum(get_offset.(c))
+minoffset(d::Dict{<:Char,<:Vector{<:AbstractComponent}}) = minimum(minoffset.(values(d)))
+
+
"""
MultichannelComponent <: AbstractComponent
@@ -139,7 +171,7 @@ function MultichannelComponent(
return MultichannelComponent(component, mg[:, ix], noise)
end
-Base.length(c::AbstractComponent) = length(c.basis)
+
Base.length(c::MultichannelComponent) = length(c.component)
"""
@@ -157,10 +189,12 @@ For `MultichannelComponent` return the length of the projection vector.
"""
n_channels(c::MultichannelComponent) = length(c.projection)
+
"""
n_channels(c::Vector{<:AbstractComponent})
+ n_channels(c::Dict)
-For a vector of `MultichannelComponent`s, return the number of channels for the first component but assert all are of equal length.
+For a vector of `MultichannelComponent`s or a Dict of components, return the number of channels for the first component but assert all are of equal length.
"""
function n_channels(c::Vector{<:AbstractComponent})
all_channels = n_channels.(c)
@@ -168,12 +202,83 @@ function n_channels(c::Vector{<:AbstractComponent})
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
+
+
+"""
+ get_basis([rng],c::AbstractComponent)
+
+Return the basis of the component (typically `c.basis`). rng is optional and ignored, but exists to have the same interface as `get_basis(c,design)`.
+"""
+get_basis(c::AbstractComponent) = c.basis
+get_basis(rng::AbstractRNG, c::AbstractComponent) = get_basis(c)
+
+
+"""
+ get_basis(c::AbstractComponent, design)
+ get_basis(rng, c::AbstractComponent, design)
+
+Evaluate the basis, if basis is a vector, directly returns it. If basis is a tuple `(f::Function, maxlength::Int)`, evaluate the function with input `design`. Cut the resulting vector or Matrix at `maxlength`.
+
+To ensure the same design can be generated, `rng` should be the global simulation `rng`. If not specified, `MersenneTwister(1)` is used.
+"""
+get_basis(basis, design) = get_basis(MersenneTwister(1), basis, design)
+get_basis(rng, c::AbstractComponent, design) = get_basis(rng, get_basis(c), design)
+get_basis(rng, b::AbstractVector, design) = b
+
+
+_get_basis_length(basis_out::AbstractMatrix) = size(basis_out, 2)
+_get_basis_length(basis_out::AbstractVector{<:AbstractVector}) = length(basis_out)
+_get_basis_length(basis_out) = error(
+ "Component basis function needs to either return a Vector of vectors or a Matrix ",
+)
+
+function get_basis(rng::AbstractRNG, basis::Tuple{Function,Int}, design)
+ f = basis[1]
+ maxlength = basis[2]
+ basis_out = applicable(f, rng, design) ? f(rng, design) : f(design)
+ l = _get_basis_length(basis_out)
+
+ @assert l == length(design) "Component basis function needs to either return a Vector of vectors or a Matrix with dim(2) == length(design) [$l / $(length(design))], or a Vector of Vectors with length(b) == length(design) [$l / $(length(design))]. "
+ limit_basis(basis_out, maxlength)
+end
+
+"""
+ limit_basis(b::AbstractVector{<:AbstractVector}, maxlength)
+
+ Cut all basis vectors to `maxlength` and pad them with 0s if they are shorter than `maxlength`.
+"""
+function limit_basis(b::AbstractVector{<:AbstractVector}, maxlength)
+
+ # first cut off after maxlength
+ b = limit_basis.(b, maxlength)
+ # now fill up with 0's
+ Δlengths = maxlength .- length.(b)
+
+ b = pad_array.(b, Δlengths, 0)
+ basis_out = reduce(hcat, b)
+
+ return basis_out
+end
+limit_basis(b::AbstractVector{<:Number}, maxlength) = b[1:min(length(b), maxlength)]
+limit_basis(b::AbstractMatrix, maxlength) = b[1:min(length(b), maxlength), :]
+
+Base.length(c::AbstractComponent) =
+ isa(get_basis(c), Tuple) ? get_basis(c)[2] : length(get_basis(c))
+
"""
maxlength(c::Vector{<:AbstractComponent}) = maximum(length.(c))
+ maxlength(components::Dict)
-Return the maximum of the individual components' lengths.
+Maximum of individual component lengths
"""
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)
@@ -189,7 +294,7 @@ simulate_component(rng, c::AbstractComponent, simulation::Simulation) =
Generate a linear model design matrix, weight it by the coefficients `c.β` and multiply the result with the given basis vector.
# Returns
-- `Matrix{Float64}`: Simulated component for each event in the events data frame. The output dimensions are `length(c.basis) x length(design)`.
+- `Matrix{Float64}`: Simulated component for each event in the events data frame. The output dimensions are `length(get_basis(c.basis)) x length(design)`.
# Examples
```julia-repl
@@ -209,21 +314,30 @@ julia> simulate_component(StableRNG(1), c, design)
"""
function simulate_component(rng, c::LinearModelComponent, design::AbstractDesign)
events = generate_events(deepcopy(rng), design)
+ X = generate_designmatrix(c.formula, events, c.contrasts)
+ y = X * c.β
+ return y' .* get_basis(deepcopy(rng), c, design)
+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
"""
@@ -236,7 +350,7 @@ Generate a MixedModel and simulate data according to the given parameters `c.β`
- `return_parameters::Bool = false`: Can be used to return the per-event parameters used to weight the basis function. Sometimes useful to inspect what is simulated.
# Returns
-- `Matrix{Float64}`: Simulated component for each event in the events data frame. The output dimensions are `length(c.basis) x length(design)`.
+- `Matrix{Float64}`: Simulated component for each event in the events data frame. The output dimensions are `length(get_basis(basis)) x length(design)`.
# Notes
1) MixedModels/Sim does not allow simulation of data without white noise of the residuals. Because we want our own noise, we use the following trick to remove the MixedModels-Noise:
@@ -311,8 +425,17 @@ function simulate_component(
rethrow(e)
end
+ @debug size(get_basis(deepcopy(rng), c, design))
# in case the parameters are of interest, we will return those, not them weighted by basis
- epoch_data_component = kron(return_parameters ? [1.0] : c.basis, m.y')
+ b = return_parameters ? [1.0] : get_basis(deepcopy(rng), c, design)
+ @debug :b, typeof(b), size(b), :m, size(m.y')
+ # in case get_basis returns a Matrix, it will be trial x time (or the transpose of it, we didnt check when writing this comment), thus we only need to scale each row by the scaling factor from the LMM
+ # in case get_basis returns a Vector of length time, it needs to be "repeated" to a trial x time matrix and then scaled again. The kronecker product efficiently does that.
+ if isa(b, AbstractMatrix)
+ epoch_data_component = ((m.y' .* b))
+ else
+ epoch_data_component = kron(b, m.y')
+ end
return epoch_data_component
#=
else
@@ -348,7 +471,7 @@ end
Return the projection of a `MultichannelComponent c` from "source" to "sensor" space.
# Returns
-- `Array{Float64,3}`: Projected simulated component for each event in the events data frame. The output dimensions are `length(c.projection) x length(c.basis) x length(design)`.
+- `Array{Float64,3}`: Projected simulated component for each event in the events data frame. The output dimensions are `length(c.projection) x length(get_basis(c)) x length(design)`.
# Examples
```julia-repl
@@ -482,13 +605,17 @@ function simulate_responses(
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(deepcopy(rng), 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, deepcopy(rng)) # TODO: `deepcopy` necessary?
end
@@ -496,17 +623,93 @@ function simulate_responses(
end
+"""
+Initializes an Array with zeros. Returns either a 2-dimensional for component-length x length(design), or a 3-D for channels x component-length x length(design)
+
+"""
+function init_epoch_data(rng, components, design)
+ max_offset = maxoffset(components)
+ min_offset = minoffset(components)
+ range_offset = (max_offset - min_offset)
+ if n_channels(components) > 1
+ epoch_data = zeros(
+ n_channels(components),
+ maxlength(components) + range_offset,
+ length(deepcopy(rng), design),
+ )
+ else
+ epoch_data = zeros(maxlength(components) + range_offset, length(rng, design))
+ end
+ return epoch_data
+end
+
+function simulate_responses(rng, event_component_dict::Dict, s::Simulation)
+ #@debug rng.state
+ epoch_data = init_epoch_data(deepcopy(rng), event_component_dict, s.design)
+ #@debug rng.state
+ evts = generate_events(deepcopy(rng), s.design)
+ #@debug rng.state
+ @debug size(epoch_data), size(evts)
+ multichannel = n_channels(event_component_dict) > 1
+ for key in keys(event_component_dict)
+ if key == '_'
+ continue
+ end
+ s_key = Simulation(
+ s.design |> x -> SubselectDesign(x, key),
+ event_component_dict,
+ s.onset,
+ s.noisetype,
+ )
+ ix = evts.event .== key
+ if multichannel
+ simulate_responses!(
+ rng,
+ @view(epoch_data[:, :, ix]),
+ event_component_dict[key],
+ s_key,
+ )
+ else
+ #@debug sum(ix), size(simulate_responses(rng, event_component_dict[key], s_key)), key
+ simulate_responses!(
+ rng,
+ @view(epoch_data[:, ix]),
+ event_component_dict[key],
+ s_key,
+ )
+ end
+ end
+ return epoch_data
+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 `epoch_data`.
"""
-function simulate_and_add!(epoch_data::AbstractMatrix, c, simulation, rng)
+function simulate_and_add!(
+ epoch_data::AbstractMatrix,
+ component::AbstractComponent,
+ simulation,
+ rng,
+)
@debug "matrix"
- @views epoch_data[1:length(c), :] .+= simulate_component(rng, c, simulation)
+
+ off = get_offset(component) - minoffset(simulation.components)
+
+
+ @views epoch_data[(1+off):(length(component)+off), :] .+=
+ simulate_component(rng, component, simulation)
end
-function simulate_and_add!(epoch_data::AbstractArray, c, simulation, rng)
+function simulate_and_add!(
+ epoch_data::AbstractArray,
+ component::AbstractComponent,
+ simulation,
+ rng,
+)
@debug "3D Array"
- @views epoch_data[:, 1:length(c), :] .+= simulate_component(rng, c, simulation)
+ off = get_offset(component) - minoffset(simulation.components)
+ @views epoch_data[:, (1+off):(length(component)+off), :] .+=
+ simulate_component(rng, component, simulation)
end
diff --git a/src/design.jl b/src/design.jl
index f50470f8..5780c0ac 100644
--- a/src/design.jl
+++ b/src/design.jl
@@ -178,16 +178,6 @@ end
#----
# Design helper functions
-"Return the dimensions of the experiment design."
-size(design::MultiSubjectDesign) = (design.n_items, design.n_subjects)
-size(design::SingleSubjectDesign) = (*(length.(values(design.conditions))...),)
-
-Base.size(design::RepeatDesign{MultiSubjectDesign}) =
- size(design.design) .* (design.repeat, 1)
-Base.size(design::RepeatDesign{SingleSubjectDesign}) = size(design.design) .* design.repeat
-
-"Length is the product of all dimensions and equals the number of events in the corresponding events dataframe."
-length(design::AbstractDesign) = *(size(design)...)
"""
apply_event_order_function(fun, rng, events)
@@ -352,7 +342,128 @@ function generate_events(rng::AbstractRNG, design::MultiSubjectDesign)
end
+# ----
+
+
+function check_sequence(s::String)
+ blankfind = findall('_', s)
+ @assert length(blankfind) <= 1 && (length(blankfind) == 0 || length(s) == blankfind[1]) "the blank-indicator '_' has to be the last sequence element, and only one can exist"
+ return s
+end
+
+
+"""
+ SequenceDesign{T} <: AbstractDesign
+
+Create a sequence of events for each entry of a provided `AbstractDesign`.
+
+The sequence string can contain any number of `char`, but the `_` character is used to indicate a break between events without any overlap and has to be at the end of the sequence string. There can only be one `_` character in a sequence string.
+Important: The exact same variable sequence is used for current rows of a design. Only, if you later nest in a `RepeatDesign` then each `RepeatDesign` repetition will gain a new variable sequence. If you need imbalanced designs, please refer to the [`ImbalancedDesign`](https://unfoldtoolbox.github.io/UnfoldDocs/UnfoldSim.jl/stable/generated/HowTo/newDesign/) tutorial.
+
+# Fields
+- `design::AbstractDesign`: The design that is generated for every sequence event.
+- `sequence::String = ""` (optional): A string of characters depicting sequences.
+ A variable sequence is defined using `[]`. For example, `S[ABC]` could result in any one sequence `SA`, `SB`, `SC`.
+ Experimental: It is also possible to define variable length sequences using `{}`. For example, `A{10,20}` would result in a sequence of 10 to 20 `A`s.
+
+# Examples
+```julia
+design = SingleSubjectDesign(conditions = Dict(:condition => ["one", "two"]))
+design = SequenceDesign(design, "SCR_")
+```
+Would result in a `generate_events(design)`
+```repl
+6×2 DataFrame
+ Row │ condition event
+ │ String Char
+─────┼──────────────────
+ 1 │ one S
+ 2 │ one C
+ 3 │ one R
+ 4 │ two S
+ 5 │ two C
+ 6 │ two R
+```
+
+## Combination of SequenceDesign and RepeatDesign
+
+### Sequence -> Repeat
+```julia
+design = SingleSubjectDesign(conditions = Dict(:condition => ["one", "two"]))
+design = SequenceDesign(design, "[AB]")
+design = RepeatDesign(design,2)
+```
+
+```julia-repl
+julia> generate_events(design)
+4×2 DataFrame
+ Row │ condition event
+ │ String Char
+─────┼──────────────────
+ 1 │ one B
+ 2 │ two B
+ 3 │ one A
+ 4 │ two A
+```
+Sequence -> Repeat: If a sequence design is repeated, then for each repetition a sequence is generated and applied. Events have different values.
+
+### Repeat -> Sequence
+```julia
+design = SingleSubjectDesign(conditions = Dict(:condition => ["one", "two"]))
+design = RepeatDesign(design,2)
+design = SequenceDesign(design, "[AB]")
+```
+
+```julia-repl
+julia> generate_events(design)
+4×2 DataFrame
+ Row │ condition event
+ │ String Char
+─────┼──────────────────
+ 1 │ one B
+ 2 │ two B
+ 3 │ one B
+ 4 │ two B
+```
+Repeat -> Sequence: The design is first repeated, then for that design one sequence is generated and applied. All events are the same.
+
+See also [`SingleSubjectDesign`](@ref), [`MultiSubjectDesign`](@ref), [`RepeatDesign`](@ref)
+"""
+@with_kw struct SequenceDesign{T} <: AbstractDesign
+ design::T
+ sequence::String = ""
+ SequenceDesign{T}(d, s) where {T<:AbstractDesign} = new(d, check_sequence(s))
+end
+
+generate_events(rng, design::SequenceDesign{MultiSubjectDesign}) =
+ error("not yet implemented")
+
+
"""
+ generate_events(rng, design::SequenceDesign)
+
+First generates a sequence-string using `evaluate_sequencestring(rng,design.sequence)`. Then generates events from the nested `design.design` and repeats them according to the length of the sequence string.
+Finally, assigns the sequence-string-characters to the `:event` column in `events`.
+"""
+function generate_events(rng, design::SequenceDesign)
+ df = generate_events(deepcopy(rng), design.design)
+ nrows_df = size(df, 1)
+
+ # @debug design.sequence
+ currentsequence = evaluate_sequencestring(rng, design.sequence)
+ # @debug currentsequence
+ currentsequence = replace(currentsequence, "_" => "")
+ df = repeat(df, inner = length(currentsequence))
+
+ df.event .= repeat(collect(currentsequence), nrows_df)
+
+ return df
+
+end
+
+
+"""
+
UnfoldSim.generate_events([rng::AbstractRNG, ]design::RepeatDesign{T})
For a `RepeatDesign`, iteratively call `generate_events` for the underlying {T} design and concatenate the results.
@@ -363,9 +474,131 @@ As a result, the order of events will differ for each repetition.
"""
function UnfoldSim.generate_events(rng::AbstractRNG, design::RepeatDesign)
df = map(x -> generate_events(rng, design.design), 1:design.repeat) |> x -> vcat(x...)
+
if isa(design.design, MultiSubjectDesign)
sort!(df, [:subject])
end
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(rng, design::SubselectDesign)
+ return subset(generate_events(rng, design.design), :event => x -> x .== design.key)
+end
+
+
+
+# ---
+# Effects
+
+"""
+ EffectsDesign <: AbstractDesign
+
+This design evaluates the nested design at the marginalized effects specified in the `effects_dict`. That is, it calculates all combination of the variables in the `effects_dict` while setting all other variables to a "typical" value (i.e. the mean for numerical variables).
+
+# Fields
+- `design::AbstractDesign`: The design of your (main) simulation.
+- `effects_dict::Dict`: Effects.jl style dictionary specifying variable effects. See also [Unfold.jl marginal effects](https://unfoldtoolbox.github.io/UnfoldDocs/Unfold.jl/stable/generated/HowTo/effects/)
+
+# Examples
+```julia-repl
+design =
+ SingleSubjectDesign(;
+ conditions = Dict(
+ :condition => ["bike", "face"],
+ :continuous => range(0, 5, length = 10),
+ ),
+ ) |> x -> RepeatDesign(x, 5);
+
+effects_dict = Dict(:condition => ["bike", "face"])
+effects_design = EffectsDesign(design, effects_dict)
+```
+"""
+struct EffectsDesign <: AbstractDesign
+ design::AbstractDesign
+ effects_dict::Dict
+end
+EffectsDesign(design::MultiSubjectDesign, effects_dict::Dict) = error("not yet implemented")
+UnfoldSim.size(rng::AbstractRNG, t::EffectsDesign) = size(generate_events(rng, t), 1)
+
+"""
+ expand_grid(design)
+
+calculate all possible combinations of the key/value pairs of the design-dict. Copied from Effects.jl
+"""
+function expand_grid(design)
+ colnames = tuple(Symbol.(keys(design))...)
+ rowtab = NamedTuple{colnames}.(Base.Iterators.product(values(design)...))
+
+ return DataFrame(vec(rowtab))
+end
+
+typical_value(v::Vector{<:Number}) = [mean(v)]
+typical_value(v) = unique(v)
+
+"""
+ generate_events(rng::AbstractRNG, design::EffectsDesign)
+
+Generate events to simulate marginal effects using an Effects.jl reference-grid dictionary. Every covariate that is in the `EffectsDesign` but not in the `effects_dict` will be set to a `typical_value` (i.e. the mean).
+
+# Examples
+```julia-repl
+julia> effects_dict = Dict(:conditionA => [0, 1]);
+
+julia> design = SingleSubjectDesign(; conditions = Dict(:conditionA => [0, 1, 2]));
+
+julia> eff_design = EffectsDesign(design,effects_dict);
+
+julia> generate_events(MersenneTwister(1),eff_design)
+2×1 DataFrame
+ Row │ conditionA
+ │ Int64
+─────┼────────────
+ 1 │ 0
+ 2 │ 1
+```
+"""
+function generate_events(rng::AbstractRNG, t::EffectsDesign)
+ effects_dict = Dict{Any,Any}(t.effects_dict)
+ #effects_dict = t.effects_dict
+ current_design = generate_events(deepcopy(rng), t.design)
+ to_be_added = setdiff(names(current_design), string.(keys(effects_dict)))
+ for tba in to_be_added
+ effects_dict[tba] = typical_value(current_design[:, tba])
+ end
+ return expand_grid(effects_dict)
+end
+
+
+"""
+ Base.size([rng],design::AbstractDesign)
+
+Return the dimensions of the experiment design. For some designs (SequenceDesign), rng is required, as the size is only determined by the design, which in turn can be probabilistic.
+"""
+Base.size(design::MultiSubjectDesign) = (design.n_items, design.n_subjects)
+Base.size(design::SingleSubjectDesign) = (*(length.(values(design.conditions))...),)
+
+Base.size(design::RepeatDesign{MultiSubjectDesign}) =
+ size(design.design) .* (design.repeat, 1)
+Base.size(design::RepeatDesign{SingleSubjectDesign}) = size(design.design) .* design.repeat
+
+Base.size(rng::AbstractRNG, design::AbstractDesign) = size(design) # by default we drop the RNG
+
+Base.size(
+ rng::AbstractRNG,
+ design::Union{<:SequenceDesign,<:SubselectDesign,<:RepeatDesign{<:SequenceDesign}},
+) = size(generate_events(rng, design), 1)
+
+
+"Length is the product of all dimensions and equals the number of events in the corresponding events dataframe."
+length(design::AbstractDesign) = *(size(design)...)
+length(rng::AbstractRNG, design::AbstractDesign) = *(size(rng, design)...)
diff --git a/src/onset.jl b/src/onset.jl
index d313bf8b..0b013f48 100644
--- a/src/onset.jl
+++ b/src/onset.jl
@@ -21,7 +21,7 @@ UniformOnset
offset: Int64 5
```
-See also [`LogNormalOnset`](@ref), [`NoOnset`](@ref).
+See also [`LogNormalOnset`](@ref UnfoldSim.LogNormalOnset), [`NoOnset`](@ref).
"""
@with_kw struct UniformOnset <: AbstractOnset
width = 50 # how many samples jitter?
@@ -51,7 +51,7 @@ LogNormalOnset
truncate_upper: Int64 25
```
-See also [`UniformOnset`](@ref), [`NoOnset`](@ref).
+See also [`UniformOnset`](@ref UnfoldSim.UniformOnset), [`NoOnset`](@ref).
"""
@with_kw struct LogNormalOnset <: AbstractOnset
μ::Any # mean
@@ -71,11 +71,43 @@ julia> onset_distribution = NoOnset()
NoOnset()
```
-See also [`UniformOnset`](@ref), [`LogNormalOnset`](@ref).
+See also [`UniformOnset`](@ref UnfoldSim.UniformOnset), [`LogNormalOnset`](@ref UnfoldSim.LogNormalOnset).
"""
struct NoOnset <: AbstractOnset end
+"""
+ ShiftOnsetByOne <:AbstractOnset
+
+This container AbstractOnset shifts the ShiftOnsetByOne.onset::AbstractOnset inter-onset-distance vector by one, adding a `0` to the front and removing the last `inter-onset distance`.
+
+This is helpful in combination with [`LogNormalOnsetFormula`](@ref) or [`UniformOnsetFormula`](@ref), to generate biased distances not of the previous, but of the next event.
+
+Visualized:
+
+|__1__| A |__2__| B |__3__| C
+Right now, the inter-onset distances are assigned in the order 1,2,3 inbetween the events A,B,C. After ShiftOnsetByOne we would have
+
+|__0__| A |__1__| B |__2__| C
+
+with 0 being a new distance of `0`, and the 3 removed (it would describe the distance after C, because there is nothing coming, the signal is not further prolonged).
+
+
+# Examples
+```julia-repl
+julia> o = UniformOnset(10,20)
+julia> d = SingleSubjectDesign(conditions=Dict(:trial=>[1,2,3]))
+julia> simulate_interonset_distances(MersenneTwister(1),o,d)'
+> 26 30 27
+julia> simulate_interonset_distances(MersenneTwister(1),ShiftOnsetByOne(o),d)'
+> 0 26 30
+```
+
+"""
+struct ShiftOnsetByOne <: AbstractOnset
+ onset::AbstractOnset
+end
+
#-----------------------------
# Onset simulation functions
#-----------------------------
@@ -125,12 +157,18 @@ function simulate_interonset_distances end
function simulate_interonset_distances(rng, onset::UniformOnset, design::AbstractDesign)
return Int.(
- round.(rand(deepcopy(rng), onset.offset:(onset.offset+onset.width), size(design)))
+ round.(
+ rand(
+ deepcopy(rng),
+ onset.offset:(onset.offset+onset.width),
+ size(deepcopy(rng), design),
+ ),
+ ),
)
end
function simulate_interonset_distances(rng, onset::LogNormalOnset, design::AbstractDesign)
- s = size(design)
+ s = size(deepcopy(rng), design)
fun = LogNormal(onset.μ, onset.σ)
if !isnothing(onset.truncate_upper)
fun = truncated(fun; upper = onset.truncate_upper)
@@ -139,6 +177,13 @@ function simulate_interonset_distances(rng, onset::LogNormalOnset, design::Abstr
end
+"""
+Returns true if the design, or any nested design contains the target design type
+"""
+contains_design(d::AbstractDesign, target::Type) = false
+contains_design(d::Union{RepeatDesign,SequenceDesign,SubselectDesign}, target::Type) =
+ (d isa target || d.design isa target) ? true : contains_design(d.design, target)
+
"""
simulate_onsets(rng, onset::AbstractOnset, simulation::Simulation)
@@ -148,6 +193,8 @@ Call `simulate_interonset_distances` to generate distances between events and th
Please note that this function is mainly for internal use in the context of `simulate` function calls. \n
Also note that the accumulation of onsets starts at 1 to avoid indexing problems in the case that the first sampled onset is 0.
+In case of a SequenceDesign with a '_' no-overlap indicator, we use twice the `maxlength(components)` as the distance following that sequence character.
+
# Arguments
- `rng`: Random number generator (RNG) to make the process reproducible.
- `onset::AbstractOnset`: Inter-onset distance distribution which is passed to `simulate_interonset_distances`.
@@ -189,11 +236,181 @@ function simulate_onsets(rng, onset::AbstractOnset, simulation::Simulation)
# sample different onsets
onsets = simulate_interonset_distances(rng, onset, simulation.design)
+
+ if contains_design(simulation.design, SequenceDesign)
+ currentsequence = evaluate_sequencestring(deepcopy(rng), simulation.design)
+ if !isnothing(findfirst("_", currentsequence))
+
+ @assert currentsequence[end] == '_' "the blank-indicator '_' has to be the last sequence element"
+ df = generate_events(deepcopy(rng), simulation.design)
+ nrows_df = size(df, 1)
+ stepsize = length(currentsequence) - 1
+ # add to every stepsize onset the maxlength of the response
+ #@debug onsets[stepsize:stepsize:end]
+ @debug stepsize
+ onsets[(stepsize+1):stepsize:end] .+= 2 .* maxlength(simulation.components)
+ #@debug onsets[stepsize:stepsize:end]
+ end
+ end
+
if maximum(onsets) > 10000
@warn "Maximum of inter-event distances was $(maximum(onsets)) - are you sure this is what you want?"
end
# accumulate them
onsets_accum = accumulate(+, onsets, dims = 1, init = 1)
+ # If the minimum component offset is negative, the onsets are shifted towards later in time to avoid that a component starts before the continuous signal starts.
+ onsets_accum = onsets_accum .- min(minoffset(simulation.components), 0)
return onsets_accum
end
+
+"""
+ simulate_interonset_distances(rng, onsets::ShiftOnsetByOne, design)
+
+Same functionality as `simulate_interonset_distances(rng,onsets::AbstractOnset)` except that it shifts the resulting vector by one, adding a `0` to the front and removing the last simuluated distance.
+"""
+UnfoldSim.simulate_interonset_distances(rng, onsets::ShiftOnsetByOne, design) =
+ vcat(0, UnfoldSim.simulate_interonset_distances(rng, onsets.onset, design)[1:(end-1)])
+
+
+"""
+ UniformOnsetFormula <: AbstractOnset
+
+Provide a Uniform Distribution for the inter-event distances, but with regression formulas for the distribution's parameters `offset` and `width`.
+
+This is helpful if your overlap/event-distribution should be dependent on some condition, e.g. more overlap in cond = 'A' than cond = 'B'.
+`Offset` affects the minimal distance. The maximal distance is `offset + width`.
+
+# Fields
+
+- `offset_formula = @formula(0~1)`: Choose a formula depending on your `design`.
+- `offset_β::Vector = [0] `(optional): Choose a `Vector` of betas. The number of betas needs to fit the formula chosen.
+- `offset_contrasts::Dict = Dict()` (optional): Choose a contrasts-`Dict`ionary according to the StatsModels specifications.
+- `width_formula = `@formula(0~1)`: Choose a formula depending on your `Design`.
+- `width_β::Vector`: Choose a `Vector` of betas, number needs to fit the formula chosen.
+- `width_contrasts::Dict = Dict()` (optional) : Choose a contrasts-`Dict`ionary according to the StatsModels specifications.
+
+# Combined with [ShiftOnsetByOne](@ref)
+Sometimes one wants to bias not the inter-onset distance prior to the current event, but after the current event.
+This is possible by using `ShiftOnsetByOne(UniformOnset(...))`, effectively shifting the inter-onset-distance vector by one. See `?ShiftOnsetByOne` for a visualization.
+
+
+# Examples
+```julia-repl
+julia> o = UnfoldSim.UniformOnsetFormula(
+ width_formula = @formula(0 ~ 1 + cond),
+ width_β = [50, 20],
+ )
+UniformOnsetFormula
+ width_formula: StatsModels.FormulaTerm{StatsModels.ConstantTerm{Int64}, Tuple{StatsModels.ConstantTerm{Int64}, StatsModels.Term}}
+ width_β: Array{Int64}((2,)) [50, 20]
+ width_contrasts: Dict{Any, Any}
+ offset_formula: StatsModels.FormulaTerm{StatsModels.ConstantTerm{Int64}, StatsModels.ConstantTerm{Int64}}
+ offset_β: Array{Int64}((1,)) [0]
+ offset_contrasts: Dict{Any, Any}
+```
+
+See also [`UniformOnset`](@ref UnfoldSim.UniformOnset) for a simplified version without linear regression specifications.
+"""
+@with_kw struct UniformOnsetFormula <: AbstractOnset
+ width_formula = @formula(0 ~ 1)
+ width_β::Vector
+ width_contrasts::Dict = Dict()
+ offset_formula = @formula(0 ~ 1)
+ offset_β::Vector = [0]
+ offset_contrasts::Dict = Dict()
+end
+
+
+function simulate_interonset_distances(rng, o::UniformOnsetFormula, design::AbstractDesign)
+ events = generate_events(deepcopy(rng), design)
+ widths =
+ UnfoldSim.generate_designmatrix(o.width_formula, events, o.width_contrasts) *
+ o.width_β
+ offsets =
+ UnfoldSim.generate_designmatrix(o.offset_formula, events, o.offset_contrasts) *
+ o.offset_β
+
+ return Int.(
+ round.(reduce(vcat, rand.(deepcopy(rng), range.(offsets, offsets .+ widths), 1)))
+ )
+end
+
+
+"""
+
+ LogNormalOnsetFormula <: AbstractOnset
+
+Provide a Log-normal Distribution of the inter-event distances, but with regression formulas for the distribution's parameters `offset`, `μ` and `σ`.
+
+This is helpful if your overlap/event-distribution should be dependent on some condition, e.g. more overlap in cond = 'A' than cond = 'B'.
+
+μ: The mean of the log-transformed variable (the log-normal random variable's logarithm follows a normal distribution).
+σ: The standard deviation of the log-transformed variable.
+offset: The minimal distance between events - aka a shift of the LogNormal distribution.
+
+# Fields
+
+- `μ_formula = @formula(0~1)` (optional): Choose a formula depending on your `design`
+- `μ_β::Vector`: Choose a `Vector` of betas, number needs to fit the formula chosen.
+- `μ_contrasts::Dict = Dict()` (optional): Choose a contrasts-`Dict`ionary according to the StatsModels specifications.
+- `σ_formula = @formula(0~1)` (optional): Choose a formula depending on your `Design`.
+- `σ_β::Vector`: Choose a `Vector` of betas, number needs to fit the formula chosen.
+- `σ_contrasts::Dict = Dict()` (optional) : Choose a contrasts-`Dict`ionary according to the StatsModels specifications.
+- `offset_formula = @formula(0~1)` (optional): Choose a formula depending on your `design` for the offset.
+- `offset_β::Vector = [0] ` (optional): Choose a `Vector` of betas. The number of betas needs to fit the formula chosen.
+- `offset_contrasts::Dict = Dict()` (optional): Choose a contrasts-`Dict`ionary according to the StatsModels specifications.
+- `truncate_upper::nothing` (optional): Upper limit (in samples) at which the distribution is truncated (formula for truncation currently not implemented)
+
+# Combined with [ShiftOnsetByOne](@ref)
+
+Sometimes one wants to bias not the inter-onset distance prior to the current event, but after the current event.
+This is possible by using `ShiftOnsetByOne(LogNormalOnset(...))`, effectively shifting the inter-onset-distance vector by one. See `?ShiftOnsetByOne` for a visualization.
+
+
+# Examples
+```julia-repl
+julia> o = LogNormalOnsetFormula(
+ σ_formula = @formula(0 ~ 1 + cond),
+ σ_β = [0.25, 0.5],
+ μ_β = [2],
+)
+```
+
+See also [`LogNormalOnset`](@ref UnfoldSim.LogNormalOnset) for a simplified version without linear regression specifications.
+"""
+@with_kw struct LogNormalOnsetFormula <: AbstractOnset
+ μ_formula = @formula(0 ~ 1)
+ μ_β::Vector
+ μ_contrasts::Dict = Dict()
+ σ_formula = @formula(0 ~ 1)
+ σ_β::Vector
+ σ_contrasts::Dict = Dict()
+ offset_formula = @formula(0 ~ 1)
+ offset_β::Vector = [0]
+ offset_contrasts::Dict = Dict()
+ truncate_upper = nothing # truncate at some sample?
+end
+
+function simulate_interonset_distances(
+ rng,
+ o::LogNormalOnsetFormula,
+ design::AbstractDesign,
+)
+ events = generate_events(deepcopy(rng), design)
+
+
+ μs = UnfoldSim.generate_designmatrix(o.μ_formula, events, o.μ_contrasts) * o.μ_β
+ σs = UnfoldSim.generate_designmatrix(o.σ_formula, events, o.σ_contrasts) * o.σ_β
+ offsets =
+ UnfoldSim.generate_designmatrix(o.offset_formula, events, o.offset_contrasts) *
+ o.offset_β
+
+
+ funs = LogNormal.(μs, σs)
+ if !isnothing(o.truncate_upper)
+ funs = truncated.(funs; upper = o.truncate_upper)
+ end
+ #@debug reduce(hcat, rand.(deepcopy(rng), funs, 1))
+ return Int.(round.(offsets .+ reduce(vcat, rand.(deepcopy(rng), funs, 1))))
+end
diff --git a/src/predefinedSimulations.jl b/src/predefinedSimulations.jl
index e269c5dc..0038a264 100644
--- a/src/predefinedSimulations.jl
+++ b/src/predefinedSimulations.jl
@@ -90,9 +90,9 @@ function predef_eeg(
# component / signal
sfreq = 100,
- p1 = (p100(; sfreq = sfreq), @formula(0 ~ 1), [5], Dict()),
- n1 = (n170(; sfreq = sfreq), @formula(0 ~ 1 + condition), [5, 3], Dict()),
- p3 = (p300(; sfreq = sfreq), @formula(0 ~ 1 + continuous), [5, 1], Dict()),
+ p1 = (p100(; sfreq = sfreq), @formula(0 ~ 1), [5], Dict(), 0),
+ n1 = (n170(; sfreq = sfreq), @formula(0 ~ 1 + condition), [5, 3], Dict(), 0),
+ p3 = (p300(; sfreq = sfreq), @formula(0 ~ 1 + continuous), [5, 1], Dict(), 0),
kwargs...,
)
@@ -124,7 +124,7 @@ function predef_eeg(
kwargs...,
)
- components = []
+ components = AbstractComponent[]
for c in comps
append!(components, [T(c...)])
end
diff --git a/src/sequence.jl b/src/sequence.jl
new file mode 100644
index 00000000..602940a2
--- /dev/null
+++ b/src/sequence.jl
@@ -0,0 +1,94 @@
+"""
+ rand_re(rng::AbstractRNG, machine::Automa.Machine)
+
+Mimicks a reverse-regex, generating strings from regex instead of matching. Based on Automata.jl
+
+# Arguments
+- `machine::Automa.Machine`: A Automa.Machine, typically output of `Automa.Compile(RE("mystring"))`
+
+# Returns
+- `result::String` : A string following the rules in `Automa.Machine`. `{}` are not supported, but e.g. `+`, `*`
+
+# Examples
+```julia-repl
+julia> using Automa
+julia> machine = Automa.compile(Automa.RegExp.RE("b+l+a+"))
+julia> rand_re(MersenneTwister(2),machine)
+"bbbbblllaaa"
+```
+"""
+function rand_re(rng::AbstractRNG, machine::Automa.Machine)
+ out = IOBuffer()
+ node = machine.start
+
+ while true
+ if node.state ∈ machine.final_states
+ (rand(rng) ≤ 1 / (length(node.edges) + 1)) && break
+ end
+
+ edge, node = rand(rng, node.edges)
+ label = rand(rng, collect(edge.labels))
+ print(out, Char(label))
+ end
+
+ return String(take!(out))
+end
+
+evaluate_sequencestring(rng, d::SequenceDesign) = evaluate_sequencestring(rng, d.sequence)
+
+
+"""
+ evaluate_sequencestring(rng, str::String)
+ evaluate_sequencestring(rng, dS::SequenceDesign)
+ evaluate_sequencestring(rng, dR::RepeatDesign)
+
+Generate a sequence based on the reverse regex style string in `str`, `dS.sequence` or `dR.design.sequence`.
+
+Directly converting to Automa.Compileis not possible, as we first need to match & evaluate the curly brackets. We simply detect and expand them.
+
+# Arguments
+- `str::String`: a string mimicking a regex, e.g. "b[lL]{3,4}a" should evaluate to e.g. "bLlLLa". E.g. "b+l*a{3,4}" should in principle evaluate to "bbbbaaa" or "bllllllllllllaaaa" - but right now we disallow `+` and `*` - we should revisit why exactly though.
+
+# Returns
+- `result::String` : a simulated string
+
+# Examples
+```julia-repl
+julia> evaluate_sequencestring(MersenneTwister(1),"bla{3,4}")
+"blaaaa"
+```
+
+See also [`rand_re`](@ref)
+"""
+function evaluate_sequencestring(rng, str::String)
+ #match curly brackets and replace them
+ @assert isnothing(findfirst("*", str)) && isnothing(findfirst("+", str)) "'Infinite' sequences are currently not supported."
+ crly = collect(eachmatch(r"(\{[\d]+,[\d]+\})", str))
+ for c in reverse(crly)
+ m = replace(c.match, "{" => "", "}" => "")
+ rep_minimum, rep_maximum = parse.(Int, split(m, ","))
+ #@info str[c.offset-1]
+ if str[c.offset-1] == ']'
+ #@info "brackets found"
+ bracket_end_idx = c.offset - 1
+ bracket_start_idx = findlast("[", str[1:bracket_end_idx])[1]
+ #@info bracket_end_idx,bracket_start_idx
+ repeat_string = "[" * str[bracket_start_idx+1:bracket_end_idx-1] * "]"
+ else
+ bracket_start_idx = c.offset - 1
+ bracket_end_idx = c.offset - 1
+ repeat_string = string(str[c.offset-1])
+ end
+
+ replacement_string = repeat(repeat_string, rand(rng, rep_minimum:rep_maximum))
+ #@info "rep" replacement_string
+ str =
+ str[1:bracket_start_idx-1] *
+ replacement_string *
+ str[bracket_end_idx+length(c.match)+1:end]
+ # @debug str
+ end
+ return rand_re(rng, Automa.compile(RE(str)))
+end
+
+evaluate_sequencestring(rng, d::RepeatDesign) = evaluate_sequencestring(rng, d.design)
diff --git a/src/simulation.jl b/src/simulation.jl
index 07539dc0..098d87ae 100755
--- a/src/simulation.jl
+++ b/src/simulation.jl
@@ -18,6 +18,7 @@ function simulate(
simulate(MersenneTwister(1), design, components, onset, args...; kwargs...)
end
+
"""
simulate(
rng::AbstractRNG,
@@ -127,7 +128,6 @@ simulate(
kwargs...,
) = simulate(rng, Simulation(design, components, onset, noise); kwargs...)
-
function simulate(rng::AbstractRNG, simulation::Simulation; return_epoched::Bool = false)
(; design, components, onset, noisetype) = simulation
@@ -145,7 +145,8 @@ function simulate(rng::AbstractRNG, simulation::Simulation; return_epoched::Bool
# such that the resulting dimensions are dimensions: channels x times x trials x subjects
# TODO: This assumes a balanced design, but create_continuous_signal also assumes this, so we should be fine ;)
size_responses = size(responses)
- signal = reshape(responses, size_responses[1:end-1]..., size(design)...)
+ signal =
+ reshape(responses, size_responses[1:(end-1)]..., size(deepcopy(rng), design)...)
else # if there is an onset distribution given the next step is to create a continuous signal
signal, latencies = create_continuous_signal(deepcopy(rng), responses, simulation)
events.latency = latencies
@@ -241,8 +242,9 @@ function create_continuous_signal(rng, responses, simulation)
(; design, components, onset, noisetype) = simulation
- n_subjects = length(size(design)) == 1 ? 1 : size(design)[2]
- n_trials = size(design)[1]
+ n_subjects =
+ length(size(deepcopy(rng), design)) == 1 ? 1 : size(deepcopy(rng), design)[2]
+ n_trials = size(deepcopy(rng), design)[1]
n_chan = n_channels(components)
# we only need to simulate onsets & pull everything together, if we
@@ -252,13 +254,26 @@ function create_continuous_signal(rng, responses, simulation)
# flatten onsets (since subjects are concatenated in the events df)
latencies = onsets[:,]
- # combine responses with onsets
+ # calculate the required length of the continuous signal
+
+ # Reasoning:
+ # (1) We want the last event onset time to be within the signal -> lowerbound of max_length_continuous is `maximum(onsets)`
+ # (2) We want to extend the signal for those cases, where the response is longer than the last event onset time, without offset, the upper bound is: maximum(onsets)+max_length_component
+ # (3) In cases where we have a positive offset, the largest offset needs to be added => maximum(onsets) + max_length_component + max(maxoffset(components), 0)
+ # (4) In cases where we have a negative offset, `max_length_component` might be reduced, by maximally the largest minoffset => maximum(onsets) + max_length_component + max(maxoffset(components), 0) + maximum(min.(get_offset.(components),0))
+
+ last_onset = maximum(onsets)
max_length_component = maxlength(components)
- max_length_continuoustime = Int(ceil(maximum(onsets))) .+ max_length_component
+ calculated_onset = maximum(onsets) + max_length_component # add the signal (ideally, we'd add the longest signal of the last event - but it's not so easy). (2)
+ calculated_onset += max(maxoffset(components), 0) # if the largest offset is positive, add it (3)
+ calculated_onset += maximum(min.(vcat(values(get_offset(components))...), 0)) # add maximum of offsets that is smaller than 0 (4)
+
+ max_length_continuoustime = max(last_onset, calculated_onset) # ensure that maximum(onsets) is lowerbound (1)
signal = zeros(n_chan, max_length_continuoustime, n_subjects)
+ # combine responses with onsets
for e = 1:n_chan
for s = 1:n_subjects
for i = 1:n_trials
@@ -268,7 +283,9 @@ function create_continuous_signal(rng, responses, simulation)
responses,
e,
s,
- one_onset:one_onset+max_length_component-1,
+ (one_onset+minoffset(simulation.components)):(one_onset+max_length_component-1+maxoffset(
+ simulation.components,
+ )),
(s - 1) * n_trials + i,
)
end
diff --git a/src/types.jl b/src/types.jl
index 314f32ef..f0b80467 100644
--- a/src/types.jl
+++ b/src/types.jl
@@ -65,7 +65,18 @@ julia> data
"""
struct Simulation
design::AbstractDesign
- components::Vector{AbstractComponent}
+ components::Union{
+ <:Dict{<:Char,<:Vector{<:AbstractComponent}},
+ <:Vector{<:AbstractComponent},
+ }
onset::AbstractOnset
noisetype::AbstractNoise
end
+
+
+Simulation(
+ design::AbstractDesign,
+ components::Dict{<:Char,<:Vector},
+ onset::AbstractOnset,
+ noisetype::AbstractNoise,
+) = Simulation(design, Dict{Char,Vector{<:AbstractComponent}}(components), onset, noisetype)
diff --git a/test/bases.jl b/test/bases.jl
new file mode 100644
index 00000000..bed6067b
--- /dev/null
+++ b/test/bases.jl
@@ -0,0 +1,16 @@
+using UnfoldSim
+@testset "hanning" begin
+ @test UnfoldSim.hanning(0.021, 0.04, 1000)[40] == 1.0
+ @test UnfoldSim.hanning(0.011, 0.04, 1000)[40] == 1.0
+ @test UnfoldSim.hanning(0.011, 0.02, 1000)[20] == 1.0
+ @test isapprox(argmax(UnfoldSim.hanning(0.021, 0.04, 256)) / 256, 0.0390625)
+ @test_throws Exception UnfoldSim.hanning(0.011, 0.0, 1000)
+end
+
+@testset "p100,N170,p300,n400" begin
+ sfreq = 1000
+ @test argmax(p100(; sfreq)) == 0.1 * sfreq
+ @test argmin(n170(; sfreq)) == 0.169 * sfreq # Why not 0.17? Because the peak of the function is in between samples (169 and 170)
+ @test argmax(p300(; sfreq)) == 0.3 * sfreq
+ @test argmin(n400(; sfreq)) == 0.4 * sfreq
+end
diff --git a/test/component.jl b/test/component.jl
index 3893a391..9c476225 100644
--- a/test/component.jl
+++ b/test/component.jl
@@ -1,5 +1,35 @@
@testset "component" begin
+ @testset "componentfunction" begin
+ design = UnfoldSim.SingleSubjectDesign(; conditions = Dict(:duration => 10:-1:5))
+
+ mybasisfun = design -> (collect.(range.(1, generate_events(design).duration)))
+ signal = LinearModelComponent(;
+ basis = (mybasisfun, 15),
+ formula = @formula(0 ~ 1),
+ β = [1],
+ )
+
+ erp = UnfoldSim.simulate_component(StableRNG(1), signal, design)
+
+ @test size(erp) == (15, 6)
+ @test all(erp[11:15, :] .== 0)
+ @test erp[1:9, 2] == collect(1.0:9)
+
+ # test shorter cut
+ signal = LinearModelComponent(;
+ basis = (mybasisfun, 5),
+ formula = @formula(0 ~ 1),
+ β = [1],
+ )
+
+ erp = UnfoldSim.simulate_component(StableRNG(1), signal, design)
+ @test size(erp) == (5, 6)
+ @test !any(erp .== 0)
+
+
+
+ end
@testset "LMM" begin
@test UnfoldSim.weight_σs(Dict(:subj => [1, 2]), 0.5, 1.0).subj ==
LowerTriangular([0.5 0; 0 1.0])
@@ -14,5 +44,155 @@
@test UnfoldSim.weight_σs(Dict(:subj => [1, 2, [1 0.5; 0.5 1]]), 1.0, 2.0).subj ==
create_re(1, 2; corrmat = [1 0.5; 0.5 1]) ./ 2
end
+ @testset "get_basis" begin
+
+ rng = StableRNG(1)
+ design = UnfoldSim.SingleSubjectDesign(; conditions = Dict(:duration => 10:-1:5))
+ mybasisfun =
+ (rng, design) -> (collect.(range.(1, generate_events(rng, design).duration)))
+ signal = LinearModelComponent(;
+ basis = (mybasisfun, 15),
+ formula = @formula(0 ~ 1),
+ β = [1],
+ )
+ @test UnfoldSim.get_basis(deepcopy(rng), signal, design) ==
+ UnfoldSim.get_basis(signal, design)
+
+ shuffle_design = UnfoldSim.SingleSubjectDesign(;
+ conditions = Dict(:duration => 10:-1:5),
+ event_order_function = shuffle,
+ )
+ # with same seed => equal result
+ @test UnfoldSim.get_basis(StableRNG(1), signal, shuffle_design) ==
+ UnfoldSim.get_basis(StableRNG(1), signal, shuffle_design)
+ # with different seed => unequal result
+ @test UnfoldSim.get_basis(StableRNG(1), signal, shuffle_design) !=
+ UnfoldSim.get_basis(StableRNG(2), signal, shuffle_design)
+
+ end
+
+ @testset "max/min offset" begin
+ # test max/min offset
+ smin10 = LinearModelComponent(;
+ basis = [1, 2, 3],
+ formula = @formula(0 ~ 1),
+ β = [1],
+ offset = -10,
+ )
+ splus5 = LinearModelComponent(;
+ basis = [1, 2, 3],
+ formula = @formula(0 ~ 1),
+ β = [1],
+ offset = 5,
+ )
+ @test UnfoldSim.get_offset(smin10) == -10
+ @test UnfoldSim.get_offset([smin10, splus5]) == [-10, 5]
+ @test UnfoldSim.maxoffset([smin10, splus5]) == 5
+ @test UnfoldSim.minoffset([smin10, splus5]) == -10
+ @test UnfoldSim.minoffset(Dict('A' => [smin10, splus5])) == -10
+ @test UnfoldSim.maxoffset(Dict('A' => [smin10, smin10], 'B' => [splus5, splus5])) ==
+ 5
+ # test that you can have a super large negative offset and don't run into errors (e.g. an event cannot even run in the issue to start before simulation time = 0)
+
+ smin10000 = LinearModelComponent(;
+ basis = [1, 2, 3],
+ formula = @formula(0 ~ 1),
+ β = [1],
+ offset = -10_000,
+ )
+ design = UnfoldSim.SingleSubjectDesign(; conditions = Dict(:duration => 10:-1:5))
+ d, e = simulate(design, smin10000, UniformOnset(50, 0))
+ @test length(d) > 10_000
+ @test e.latency[1] > 10_000
+ @test d[e.latency[1]-10_000] == 1
+ smax10000 = LinearModelComponent(;
+ basis = [1, 2, 3],
+ formula = @formula(0 ~ 1),
+ β = [1],
+ offset = +10_000,
+ )
+ d, e = simulate(design, smax10000, UniformOnset(50, 0))
+ @test length(d) > 10_000
+ @test e.latency[1] < 100
+ @test d[e.latency[1]+10_000] == 1
+
+
+ # if we go back -10_000 and front +10_000, we should get a signal measuring 20_000
+ d, e = simulate(design, [smax10000, smin10000], UniformOnset(50, 0))
+ @test length(d) > 20_000
+ @test length(d) < 25_000 # earlier tests had the signal at 30_000, a bit too long
+ @test d[e.latency[1]+10_000] == 1
+ @test d[e.latency[1]-10_000] == 1
+
+
+
+ smax1000 = LinearModelComponent(;
+ basis = [1, 2, 3],
+ formula = @formula(0 ~ 1),
+ β = [1],
+ offset = +1000,
+ )
+ smax2000 = LinearModelComponent(;
+ basis = [1, 2, 3],
+ formula = @formula(0 ~ 1),
+ β = [1],
+ offset = +2000,
+ )
+
+ d, e = simulate(design, [smax1000, smax2000], UniformOnset(50, 0))
+ @test d[e.latency[1]+1000] == 1
+ @test d[e.latency[1]+2000] == 1
+
+ smin1000 = LinearModelComponent(;
+ basis = [1, 2, 3],
+ formula = @formula(0 ~ 1),
+ β = [1],
+ offset = -1000,
+ )
+ smin2000 = LinearModelComponent(;
+ basis = [1, 2, 3],
+ formula = @formula(0 ~ 1),
+ β = [1],
+ offset = -2000,
+ )
+
+ d, e = simulate(design, [smin1000, smin2000], UniformOnset(50, 0))
+ @test d[e.latency[1]-1000] == 1
+ @test d[e.latency[1]-2000] == 1
+
+ # Sequences with component offsets
+ design =
+ SingleSubjectDesign(conditions = Dict(:condition => ["one", "two"])) |>
+ d -> RepeatDesign(SequenceDesign(d, "SR_"), 4)
+
+ components = Dict('S' => [smin10, splus5], 'R' => [smin10000])
+
+ @test UnfoldSim.get_offset(components) == Dict('R' => [-10_000], 'S' => [-10, 5])
+
+ o_width = 20
+ o_offset = 0
+ minoffset_shift = -1 * min(UnfoldSim.minoffset(components), 0) # latencies should be shifted to the right if minoffset is negative
+
+ for seed in range(1, 10)
+ d, e = simulate(
+ StableRNG(seed),
+ design,
+ components,
+ UniformOnset(offset = o_offset, width = o_width),
+ NoNoise(),
+ )
+ sequence_length =
+ length(UnfoldSim.evaluate_sequencestring(StableRNG(seed), design)) - 1 # without _
+
+ # Test onset shifts with component offsets and sequences (in particular inter-event-block distances) combined
+ @test minoffset_shift + 1 <= e.latency[1] <= minoffset_shift + 1 + o_width
+ @test minoffset_shift + 1 + (sequence_length + 1) * o_offset <=
+ e.latency[sequence_length+1] <=
+ minoffset_shift +
+ 1 +
+ (sequence_length + 1) * o_width +
+ 2 * UnfoldSim.maxlength(components) # TODO: This part will fail once we implement a different way to specify the inter-event-block distances. Should be adapted then.
+ end
+ end
end
diff --git a/test/design.jl b/test/design.jl
index 8a251887..b10879ba 100644
--- a/test/design.jl
+++ b/test/design.jl
@@ -188,4 +188,39 @@
@test length(unique(i1.cond)) == 1
end
+ @testset "Effects Design" begin
+ # Begin with simulation design
+ design = SingleSubjectDesign(;
+ conditions = Dict(
+ :condition => ["car", "face"],
+ :continuous => range(0, 5, length = 10),
+ ),
+ )
+
+ # Effects dictionary
+ effects_dict_1 = Dict(:condition => ["car", "face"])
+ effects_dict_2 = Dict(:condition => ["car", "face"], :continuous => [2, 3, 4])
+
+ # Generate effects design
+ ef_design_1 = EffectsDesign(design, effects_dict_1)
+ ef_design_2 = EffectsDesign(design, effects_dict_2)
+
+ # Generate events
+ ef_events_1 = generate_events(ef_design_1)
+ ef_events_2 = generate_events(ef_design_2)
+
+ # SingleSubject tests
+ @test size(ef_events_1, 1) == 2 # Test correct length of events df
+ @test unique(ef_events_1[!, :continuous])[1] ≈ mean(range(0, 5, length = 10)) # Test that average is calculated correctly and only one value is present in df
+ @test size(ef_events_2, 1) == 6 # Test correct length of events df when one inputs values for continuous variable
+
+ # MultiSubjectDesign -> not implemented yet, so should error
+ design = MultiSubjectDesign(
+ n_subjects = 20,
+ n_items = 8,
+ items_between = Dict(:condition => ["car", "face"], :continuous => [1, 2]),
+ )
+ @test_throws ErrorException EffectsDesign(design, effects_dict_1)
+
+ end
end
diff --git a/test/onset.jl b/test/onset.jl
index 0e6ec715..a21f08a3 100644
--- a/test/onset.jl
+++ b/test/onset.jl
@@ -60,4 +60,91 @@
@test accumulated_onset[1] >= 1
end
+ @testset "OnsetFormula" begin
+
+ design =
+ SingleSubjectDesign(conditions = Dict(:cond => ["A", "B"])) |>
+ x -> RepeatDesign(x, 10000)
+
+
+ o = UniformOnsetFormula(width_formula = @formula(0 ~ 1 + cond), width_β = [50, 20])
+ events = generate_events(design)
+ onsets = UnfoldSim.simulate_interonset_distances(StableRNG(1), o, design)
+ @test minimum(onsets[1:2:end]) == 0
+ @test maximum(onsets[1:2:end]) == 50
+ @test minimum(onsets[2:2:end]) == 0
+ @test maximum(onsets[2:2:end]) == 70
+
+ o = UniformOnsetFormula(
+ offset_formula = @formula(0 ~ 1 + cond),
+ offset_β = [50, 20],
+ width_β = [50],
+ )
+ events = generate_events(design)
+ onsets = UnfoldSim.simulate_interonset_distances(StableRNG(1), o, design)
+ @test minimum(onsets[1:2:end]) == 50
+ @test maximum(onsets[1:2:end]) == 100
+ @test minimum(onsets[2:2:end]) == 70
+ @test maximum(onsets[2:2:end]) == 120
+
+
+ o = LogNormalOnsetFormula(
+ μ_formula = @formula(0 ~ 1 + cond),
+ μ_β = [1, 1],
+ σ_β = [1],
+ )
+ events = generate_events(design)
+ onsets = UnfoldSim.simulate_interonset_distances(StableRNG(1), o, design)
+ @test minimum(onsets[1:2:end]) == 0
+ @test maximum(onsets[1:2:end]) < 150
+ @test minimum(onsets[2:2:end]) == 0
+ @test maximum(onsets[2:2:end]) > 300
+
+
+
+ end
+
+ @testset "ShiftOnset" begin
+ design =
+ SingleSubjectDesign(conditions = Dict(:cond => ["A", "B"])) |>
+ x -> RepeatDesign(x, 100)
+
+ o = UniformOnset(width = 50, offset = 10)
+
+ without = UnfoldSim.simulate_interonset_distances(StableRNG(1), o, design)
+ with = UnfoldSim.simulate_interonset_distances(
+ StableRNG(1),
+ ShiftOnsetByOne(o),
+ design,
+ )
+ # ShiftOnsetByOne adds a 0 to the front, thereby the first "non-0" "real" simulated inter onset distance is used for the second event
+ @test with[1] == 0
+
+ @test without[1:(end-1)] == with[2:end]
+
+
+ end
end
+
+@testset "contains design" begin
+ @test UnfoldSim.contains_design(
+ RepeatDesign(SequenceDesign(SingleSubjectDesign(), "ABC"), 1),
+ SequenceDesign,
+ )
+ @test UnfoldSim.contains_design(
+ RepeatDesign(SequenceDesign(SingleSubjectDesign(), "ABC"), 1),
+ SingleSubjectDesign,
+ )
+ @test UnfoldSim.contains_design(
+ RepeatDesign(SequenceDesign(SingleSubjectDesign(), "ABC"), 1),
+ RepeatDesign,
+ )
+ @test UnfoldSim.contains_design(
+ SequenceDesign(SingleSubjectDesign(), "ABC"),
+ SequenceDesign,
+ )
+ @test !UnfoldSim.contains_design(
+ SequenceDesign(SingleSubjectDesign(), "ABC"),
+ RepeatDesign,
+ )
+end
\ No newline at end of file
diff --git a/test/runtests.jl b/test/runtests.jl
index ab9284e4..8927505b 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -2,12 +2,14 @@ using UnfoldSim
include("setup.jl")
@testset "UnfoldSim.jl" begin
+ include("bases.jl")
include("component.jl")
include("design.jl")
+ include("headmodel.jl")
+ include("helper.jl")
+ include("multichannel.jl")
include("noise.jl")
include("onset.jl")
+ include("sequence.jl")
include("simulation.jl")
- include("helper.jl")
- include("headmodel.jl")
- include("multichannel.jl")
end
diff --git a/test/sequence.jl b/test/sequence.jl
new file mode 100644
index 00000000..b12cadd6
--- /dev/null
+++ b/test/sequence.jl
@@ -0,0 +1,83 @@
+using Automa, Random, Test
+@testset "Check Sequences" begin
+ @test isa(UnfoldSim.check_sequence("bla_"), String)
+ @test isa(UnfoldSim.check_sequence("bla"), String)
+ @test_throws AssertionError UnfoldSim.check_sequence("b_la_")
+ @test_throws AssertionError UnfoldSim.check_sequence("b_la")
+ @test_throws AssertionError UnfoldSim.check_sequence("_bla")
+
+ @test length(UnfoldSim.evaluate_sequencestring(StableRNG(1), "A{10,10}")) == 10
+ @test length(UnfoldSim.evaluate_sequencestring(StableRNG(1), "A{10,10}B")) == 11
+ @test length(UnfoldSim.evaluate_sequencestring(StableRNG(1), "A{10,20}")) >= 10
+end
+
+@testset "Simulate Sequences" begin
+
+
+
+ design = SingleSubjectDesign(conditions = Dict(:condition => ["one", "two"]))
+ design = SequenceDesign(design, "SCR_")
+ evt = generate_events(StableRNG(1), design)
+ @test size(evt, 1) == 6
+ @test evt.event == ['S', 'C', 'R', 'S', 'C', 'R']
+
+ design = RepeatDesign(design, 2)
+ evt = generate_events(StableRNG(1), design)
+ @test size(evt, 1) == 12
+ @test evt.event == ['S', 'C', 'R', 'S', 'C', 'R', 'S', 'C', 'R', 'S', 'C', 'R']
+
+
+ # repeat first, then sequence => same sequence
+ design = SingleSubjectDesign(conditions = Dict(:condition => ["A", "B"]))
+ design = RepeatDesign(design, 2)
+ design = SequenceDesign(design, "S[ABCD]")
+ evt = generate_events(StableRNG(2), design)
+
+ @test all(evt.event[2:2:end] .== 'B')
+
+
+ # sequence first, then repeat => different sequence for each repetition
+ design = SingleSubjectDesign(conditions = Dict(:condition => ["A", "B"]))
+ design = SequenceDesign(design, "S[ABCD]")
+ design = RepeatDesign(design, 2)
+ evt = generate_events(StableRNG(2), design)
+ @test !all(evt.event[2:2:end] .== 'B')
+end
+
+@testset "simulate_sequence" begin
+ design = SingleSubjectDesign(conditions = Dict(:condition => ["one", "two"]))
+ design = SequenceDesign(design, "SCR_")
+ c = LinearModelComponent(;
+ basis = UnfoldSim.hanning(40),
+ formula = @formula(0 ~ 1 + condition),
+ β = [1.0, 2.0],
+ contrasts = Dict(:cond => EffectsCoding()),
+ )
+ s, e = simulate(design, c, NoOnset(); return_epoched = true)
+ @test size(s) == (40, 6)
+ @test s[:, 1] == s[:, 2] # If no component dict is specified, all events have the same component
+
+end
+
+@testset "rand_re" begin
+
+ machine = Automa.compile(Automa.RegExp.RE("b+l+a+"))
+ @test UnfoldSim.rand_re(MersenneTwister(2), machine) == "bbbbblllaa"
+ # trivial single-character regex
+ @test UnfoldSim.rand_re(MersenneTwister(1), Automa.compile(Automa.RegExp.RE("a"))) ==
+ "a"
+
+ # different seeds should (very likely) produce different strings
+ r1 = UnfoldSim.rand_re(MersenneTwister(1), machine)
+ r2 = UnfoldSim.rand_re(MersenneTwister(2), machine)
+ @test r1 != r2
+
+ # produced string should match the regex pattern
+ for k = 1:10
+ s = UnfoldSim.rand_re(MersenneTwister(k), machine)
+ @test occursin(r"^b+l+a+$", s)
+ end
+
+ # integration: evaluate_sequencestring expands curly braces and delegates to rand_re
+ @test UnfoldSim.evaluate_sequencestring(MersenneTwister(1), "bla{3,4}") == "blaaaa"
+end
\ No newline at end of file
diff --git a/test/simulation.jl b/test/simulation.jl
index f2817a38..a05fd36e 100644
--- a/test/simulation.jl
+++ b/test/simulation.jl
@@ -1,3 +1,4 @@
+using Base: AbstractCartesianIndex
@testset "simulation" begin
@testset "general_test_simulate" begin
@@ -215,4 +216,37 @@
end
end
+ @testset "multi-component sequence #124" begin
+ struct MyLinearModelComponent1 <: AbstractComponent
+ comp::Any
+ end
+ MyLinearModelComponent1(b, f, β) =
+ MyLinearModelComponent1(LinearModelComponent(; basis = b, formula = f, β))
+ UnfoldSim.simulate_component(
+ rng,
+ c::MyLinearModelComponent1,
+ design::UnfoldSim.SubselectDesign,
+ ) = simulate_component(rng, c.comp, design)
+ UnfoldSim.length(c::MyLinearModelComponent1) = length(c.comp)
+ UnfoldSim.size(c::MyLinearModelComponent1) = size(c.comp)
+ sim = Simulation(
+ SingleSubjectDesign(conditions = Dict(:event => ['A', 'B'])),
+ Dict(
+ 'A' => [
+ LinearModelComponent(
+ basis = p100(),
+ formula = @formula(0 ~ 1),
+ β = [1],
+ ),
+ ],
+ 'B' => [MyLinearModelComponent1(p100(), @formula(0 ~ 1), [2])],
+ ),
+ NoOnset(),
+ NoNoise(),
+ )
+ d, e = simulate(UnfoldSim.MersenneTwister(1), sim; return_epoched = true)
+ @test d[10, 1] > 0.9 # 1 if the hanning would hit perfectly (currently the peak is between samples)
+ @test d[10, 2] > 1.9 # 2 if the hanning would hit perfectly (currently the peak is between samples)
+ end
+
end