From 8405a1a92ec00c776814ad10b93ed3a53eef877b Mon Sep 17 00:00:00 2001 From: "behinger (s-ccs 001)" Date: Fri, 23 Feb 2024 22:57:59 +0000 Subject: [PATCH 01/37] FormulaOnsets --- docs/literate/reference/onsettypes.jl | 22 +++++++ src/component.jl | 22 +++++-- src/onset.jl | 84 +++++++++++++++++++++++++++ 3 files changed, 122 insertions(+), 6 deletions(-) diff --git a/docs/literate/reference/onsettypes.jl b/docs/literate/reference/onsettypes.jl index 279a2226..e9967316 100644 --- a/docs/literate/reference/onsettypes.jl +++ b/docs/literate/reference/onsettypes.jl @@ -282,3 +282,25 @@ end # - if `offset` < `length(signal.basis)` -> there might be overlap, depending on the other parameters of the onset distribution # [^1]: Wikipedia contributors. (2023, December 5). Log-normal distribution. In Wikipedia, The Free Encyclopedia. Retrieved 12:27, December 7, 2023, from https://en.wikipedia.org/w/index.php?title=Log-normal_distribution&oldid=1188400077# + + + + +# ## Design-dependent `FormulaXOnset` + +# For additional control we provide `FormulaUniformOnset` and `FormulaLogNormalOnset` types, that allow to control all parameters by specifying formulas +o = UnfoldSim.FormulaUniformOnset( + width_formula = @formula(0 ~ 1 + cond), + width_β = [50, 20], +) +events = generate_events(design) +onsets = UnfoldSim.simulate_interonset_distances(MersenneTwister(42), o, design) + +f = Figure() +ax = f[1, 1] = Axis(f) +hist!(ax, onsets[events.cond.=="A"], bins = range(0, 100, step = 1), label = "cond: A") +hist!(ax, onsets[events.cond.=="B"], bins = range(0, 100, step = 1), label = "cond: B") +axislegend(ax) +f + +# Voila - the inter-onset intervals are `20` samples longer for condition `B`, exactly as specified.` \ No newline at end of file diff --git a/src/component.jl b/src/component.jl index a2d87fc1..ca31d99f 100644 --- a/src/component.jl +++ b/src/component.jl @@ -145,22 +145,32 @@ julia> simulate_component(StableRNG(1),c,design) """ function simulate_component(rng, c::LinearModelComponent, design::AbstractDesign) events = generate_events(design) + X = generate_designmatrix(c.formula, events, c.contrasts) + y = X * β + return y' .* c.basis +end + + +""" +Helper function to generate a designmatrix from formula, events and contrasts. +""" +function generate_designmatrix(formula, events, contrasts) # special case, intercept only # https://github.com/JuliaStats/StatsModels.jl/issues/269 - if c.formula.rhs == ConstantTerm(1) + if formula.rhs == ConstantTerm(1) X = ones(nrow(events), 1) else - if isempty(c.contrasts) - m = StatsModels.ModelFrame(c.formula, events) + if isempty(contrasts) + m = StatsModels.ModelFrame(formula, events) else - m = StatsModels.ModelFrame(c.formula, events; contrasts = c.contrasts) + m = StatsModels.ModelFrame(formula, events; contrasts = contrasts) end X = StatsModels.modelmatrix(m) end - y = X * c.β - return y' .* c.basis + return X end + """ simulate_component(rng, c::MixedModelComponent, design::AbstractDesign) Generates a MixedModel and simulates data according to c.β and c.σs. diff --git a/src/onset.jl b/src/onset.jl index 5e7d364a..22aa37b1 100644 --- a/src/onset.jl +++ b/src/onset.jl @@ -7,6 +7,8 @@ Provide a Uniform Distribution of the inter-event-distances. `width` is the width of the uniform distribution (=> the jitter). Since the lower bound is 0, `width` is also the upper bound. `offset` is the minimal distance. The maximal distance is `offset + width`. + +For a more advanced parameter specification, see `FormulaUniformOnset``, which allows to specify the onset-parameters depending on the `Design` employed via a linear regression model """ @with_kw struct UniformOnset <: AbstractOnset width = 50 # how many samples jitter? @@ -17,6 +19,8 @@ end Log-normal inter-event distances using the `Distributions.jl` truncated LogNormal distribution. Be careful with large `μ` and `σ` values, as they are on logscale. σ>8 can quickly give you out-of-memory sized signals! + +For a more advanced parameter specification, see `FormulaLogNormalOnset, which allows to specify the onset-parameters depending on the `Design` employed via linear regression model """ @with_kw struct LogNormalOnset <: AbstractOnset μ::Any # mean @@ -34,6 +38,8 @@ struct NoOnset <: AbstractOnset end """ simulate_interonset_distances(rng, onset::UniformOnset, design::AbstractDesign) simulate_interonset_distances(rng, onset::LogNormalOnset, design::AbstractDesign) + simulate_interonset_distances(rng, onset::FormulaUniformOnset, design::AbstractDesign) + simulate_interonset_distances(rng, onset::FormulaLogNormalOnset, design::AbstractDesign) Generate the inter-event-onset vector in samples (returns Int). """ @@ -69,3 +75,81 @@ function simulate_onsets(rng, onset::AbstractOnset, simulation::Simulation) return onsets_accum end + +""" + FormulaUniformOnset <: AbstractOnset +provide a Uniform Distribution of the inter-event-distances, but with regression formulas. +This is helpful if your overlap/event-distribution should be dependend on some condition, e.g. more overlap in cond='A' than cond='B'. + + - `width`: is the width of the uniform distribution (=> the jitter). Since the lower bound is 0, `width` is also the upper bound. + -`width_formula`: choose a formula depending on your `Design` + -`width_β`: Choose a vector of betas, number needs to fit the formula chosen + -`width_contrasts` (optional): Choose a contrasts-dictionary according to the StatsModels specifications + `offset` is the minimal distance. The maximal distance is `offset + width`. + -`offset_formula`: choose a formula depending on your `design` + -`offset_β`: Choose a vector of betas, number needs to fit the formula chosen + -`offset_contrasts` (optional): Choose a contrasts-dictionary according to the StatsModels specifications + +See `UniformOnset` for a simplified version without linear regression specifications +""" +@with_kw struct FormulaUniformOnset <: AbstractOnset + width_formula = @formula(0 ~ 1) + width_β::Vector = [50] + width_contrasts::Dict = Dict() + offset_formula = @formula(0 ~ 1) + offset_β::Vector = [0] + offset_contrasts::Dict = Dict() +end + + +function simulate_interonset_distances(rng, o::FormulaUniformOnset, design::AbstractDesign) + events = generate_events(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 + + +@with_kw struct FormulaLogNormalOnset <: AbstractOnset + μ_formula = @formula(0 ~ 1) + μ_β::Vector = [0] + μ_contrasts::Dict = Dict() + σ_formula = @formula(0 ~ 1) + σ_β::Vector = [0] + σ_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::FormulaLogNormalOnset, + design::AbstractDesign, +) + events = generate_events(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) + fun = truncated.(fun; upper = o.truncate_upper) + end + return Int.(round.(offsets .+ rand.(deepcopy(rng), funs, 1))) +end + + From afc69d7437547943d9c44626d6e9844ef1327c62 Mon Sep 17 00:00:00 2001 From: "behinger (s-ccs 001)" Date: Sun, 25 Feb 2024 00:00:04 +0000 Subject: [PATCH 02/37] initial sequence tryout --- Project.toml | 1 + docs/literate/HowTo/sequence.jl | 52 +++++++++++++++++++++++++++ docs/make.jl | 1 + src/UnfoldSim.jl | 4 ++- src/component.jl | 63 ++++++++++++++++++++++++++++----- src/design.jl | 44 +++++++++++++++++++++++ src/onset.jl | 23 ++++++++++++ src/sequence.jl | 48 +++++++++++++++++++++++++ src/types.jl | 5 ++- 9 files changed, 230 insertions(+), 11 deletions(-) create mode 100644 docs/literate/HowTo/sequence.jl create mode 100644 src/sequence.jl diff --git a/Project.toml b/Project.toml index 8a8913b0..4f2b1613 100644 --- a/Project.toml +++ b/Project.toml @@ -5,6 +5,7 @@ version = "0.3.1" [deps] Artifacts = "56f22d72-fd6d-98f1-02f0-08ddc0907c33" +Automa = "67c07d97-cdcb-5c2c-af73-a7f9c32a568b" DSP = "717857b8-e6f2-59f4-9121-6e50c889abd2" DataFrames = "a93c6f00-e57d-5684-b7b6-d8193f3e46c0" Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f" diff --git a/docs/literate/HowTo/sequence.jl b/docs/literate/HowTo/sequence.jl new file mode 100644 index 00000000..e66e5e95 --- /dev/null +++ b/docs/literate/HowTo/sequence.jl @@ -0,0 +1,52 @@ +using UnfoldSim +using CairoMakie +using StableRNGs + +# ## Stimulus - Response design + +# let's say we want to simulate a stimulus response, followed by a button press response. +# First we generate the minimal design of the experiment by specifying our conditins (a one-condition-two-levels design in our case) +design = SingleSubjectDesign(conditions = Dict(:condition => ["one", "two"])) +generate_events(design) +# next we use the `SequenceDesign` and nest our initial design in it. "SR_" is code for an "S" event and an "R" event - only single letter events are supported! The `_` is a signal for the Onset generator to generate a bigger pause - no overlap between adjacend `SR` pairs +design = SequenceDesign(design, "SRR_") +generate_events(design) +# The main thing that happened is that the design was repeated for every event (each 'letter') of the sequence, and an `eventtype` column was added. +# !!! hint +# more advaned sequences exist, like "SR{1,3}", or "A[BC]" etc. + +# Finally, let's repeat the design 2 times +design = RepeatDesign(design, 2) +generate_events(design) + +# This results in 20 trials that nicely follow our sequence + +# Next we have to specify for both events `S` and `R` what the responses should look like. + +p1 = LinearModelComponent(; + basis = p100(), + formula = @formula(0 ~ 1 + condition), + β = [1, 0.5], +); + +n1 = LinearModelComponent(; + basis = n170(), + formula = @formula(0 ~ 1 + condition), + β = [1, 0.5], +); +p3 = LinearModelComponent(; + basis = UnfoldSim.hanning(Int(0.5 * 100)), # sfreq = 100 for the other bases + formula = @formula(0 ~ 1 + condition), + β = [1, 2], +); + +components = Dict('S' => [p1, n1], 'R' => [p3]) + +data, evts = simulate( + StableRNG(1), + design, + components, + UniformOnset(offset = 10, width = 100), + NoNoise(), +) +lines(data) \ No newline at end of file diff --git a/docs/make.jl b/docs/make.jl index e97757f1..8cf95422 100644 --- a/docs/make.jl +++ b/docs/make.jl @@ -48,6 +48,7 @@ makedocs(; "Define a new duration & jitter component" => "./generated/HowTo/newComponent.md", "Generate multi channel data" => "./generated/HowTo/multichannel.md", "Use predefined design / onsets data" => "./generated/HowTo/predefinedData.md", + "Produce specific sequences of events" => "./generated/HowTo/sequence.md", ], "DocStrings" => "api.md", ], diff --git a/src/UnfoldSim.jl b/src/UnfoldSim.jl index cd1b0287..126264b8 100644 --- a/src/UnfoldSim.jl +++ b/src/UnfoldSim.jl @@ -14,6 +14,7 @@ using LinearAlgebra using ToeplitzMatrices # for AR Expo. Noise "Circulant" using StatsModels using HDF5, Artifacts, FileIO +using Automa # for sequence using LinearAlgebra # headmodel @@ -30,6 +31,7 @@ include("onset.jl") include("predefinedSimulations.jl") include("headmodel.jl") include("helper.jl") +include("sequence.jl") include("bases.jl") export size, length @@ -46,7 +48,7 @@ export Simulation export MixedModelComponent, LinearModelComponent # export designs -export MultiSubjectDesign, SingleSubjectDesign, RepeatDesign +export MultiSubjectDesign, SingleSubjectDesign, RepeatDesign, SequenceDesign # noise functions export PinkNoise, RedNoise, WhiteNoise, NoNoise, ExponentialNoise #,RealNoise (not implemented yet) diff --git a/src/component.jl b/src/component.jl index ca31d99f..310b4dd8 100644 --- a/src/component.jl +++ b/src/component.jl @@ -99,10 +99,16 @@ For a vector of `MultichannelComponent`s, return the first but asserts all are o """ function n_channels(c::Vector{<:AbstractComponent}) all_channels = n_channels.(c) - @assert length(unique(all_channels)) == 1 "Error - projections of different channels cannot be different from eachother" + @assert length(unique(all_channels)) == 1 "Error - projections of different components have to be of the same output (=> channel) dimension" return all_channels[1] end +function n_channels(components::Dict) + all_channels = [n_channels(c) for c in values(components)] + @assert length(unique(all_channels)) == 1 "Error - projections of different components have to be of the same output (=> channel) dimension" + return all_channels[1] + +end """ simulate_component(rng,c::MultichannelComponent,design::AbstractDesign) Return the projection of a component from source to "sensor" space. @@ -124,10 +130,13 @@ Base.length(c::AbstractComponent) = length(c.basis) """ maxlength(c::Vector{AbstractComponent}) = maximum(length.(c)) + maxlength(components::Dict) maximum of individual component lengths """ -maxlength(c::Vector{AbstractComponent}) = maximum(length.(c)) +maxlength(c::Vector{<:AbstractComponent}) = maximum(length.(c)) + +maxlength(components::Dict) = maximum([maximum(length.(c)) for c in values(components)]) """ simulate_component(rng, c::AbstractComponent, simulation::Simulation) By default call `simulate_component` with `(::Abstractcomponent,::AbstractDesign)` instead of the whole simulation. This allows users to provide a hook to do something completely different :) @@ -146,7 +155,7 @@ julia> simulate_component(StableRNG(1),c,design) function simulate_component(rng, c::LinearModelComponent, design::AbstractDesign) events = generate_events(design) X = generate_designmatrix(c.formula, events, c.contrasts) - y = X * β + y = X * c.β return y' .* c.basis end @@ -296,19 +305,55 @@ 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(components, simulation.design) + simulate_responses!(rng, epoch_data, components, simulation) + return epoch_data +end +function simulate_responses!( + rng, + epoch_data::AbstractArray, + components::Vector, + simulation::Simulation, +) for c in components simulate_and_add!(epoch_data, c, simulation, rng) end return epoch_data end +function init_epoch_data(components, design) + if n_channels(components) > 1 + epoch_data = zeros(n_channels(components), maxlength(components), length(design)) + else + epoch_data = zeros(maxlength(components), length(design)) + end + return epoch_data +end +function simulate_responses(rng, components::Dict, s::Simulation) + epoch_data = init_epoch_data(components, s.design) + evts = generate_events(s.design) + multichannel = n_channels(components) > 1 + for key in keys(components) + if key == '_' + continue + end + s_key = Simulation( + s.design |> x -> SubselectDesign(x, key), + components[key], + s.onset, + s.noisetype, + ) + ix = evts.event .== key + if multichannel + simulate_responses!(rng, @view(epoch_data[:, :, ix]), components[key], s_key) + else + #@debug sum(ix), size(simulate_responses(rng, components[key], s_key)), key + simulate_responses!(rng, @view(epoch_data[:, ix]), components[key], s_key) + end + end + return epoch_data +end """ simulate_and_add!(epoch_data::AbstractMatrix, c, simulation, rng) diff --git a/src/design.jl b/src/design.jl index fffff39a..1ea0b177 100644 --- a/src/design.jl +++ b/src/design.jl @@ -147,6 +147,28 @@ design = RepeatDesign(designOnce,4); repeat::Int = 1 end + +@with_kw struct SequenceDesign{T} <: AbstractDesign + design::T + sequence::String = "" +end + +UnfoldSim.generate_events(design::SequenceDesign{MultiSubjectDesign}) = + error("not yet implemented") + +function UnfoldSim.generate_events(design::SequenceDesign) + df = generate_events(design.design) + nrows_df = size(df, 1) + currentsequence = sequencestring(design.sequence) + currentsequence = replace(currentsequence, "_" => "") + df = repeat(df, inner = length(currentsequence)) + df.event .= repeat(collect(currentsequence), nrows_df) + + return df + +end + + """ UnfoldSim.generate_events(design::RepeatDesign{T}) @@ -160,6 +182,28 @@ function UnfoldSim.generate_events(design::RepeatDesign) return df end + + +""" +Internal helper design to subset a sequence design in its individual components +""" +struct SubselectDesign{T} <: AbstractDesign + design::T + key::Char +end + +function generate_events(design::SubselectDesign) + return subset(generate_events(design.design), :event => x -> x .== design.key) +end + + Base.size(design::RepeatDesign{MultiSubjectDesign}) = size(design.design) .* (design.repeat, 1) Base.size(design::RepeatDesign{SingleSubjectDesign}) = size(design.design) .* design.repeat +Base.size(design::SequenceDesign) = + size(design.design) .* length(replace(design.sequence, "_" => "")) + +Base.size(design::RepeatDesign{SequenceDesign{SingleSubjectDesign}}) = + size(design.design) .* design.repeat + +Base.size(design::SubselectDesign) = size(generate_events(design), 1) \ No newline at end of file diff --git a/src/onset.jl b/src/onset.jl index 22aa37b1..4ef74702 100644 --- a/src/onset.jl +++ b/src/onset.jl @@ -59,6 +59,14 @@ function simulate_interonset_distances(rng, onset::LogNormalOnset, design::Abstr end +#function simulate_interonset_distances(rng, onset::AbstractOnset,design::) + + +contains_design(d::AbstractDesign, target::Type) = false +contains_design(d::Union{RepeatDesign,SequenceDesign,SubselectDesign}, target::Type) = + d.design isa target ? true : contains_design(d.design, target) + +sequencestring(d::RepeatDesign) = sequencestring(d.design) """ simulate_onsets(rng, onset::AbstractOnset, simulation::Simulation) @@ -70,6 +78,21 @@ 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 = sequencestring(simulation.design) + if !isnothing(findfirst("_", currentsequence)) + + @assert currentsequence[end] == '_' "the blank-indicator '_' has to be the last sequence element" + df = generate_events(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 # accumulate them onsets_accum = accumulate(+, onsets, dims = 1) diff --git a/src/sequence.jl b/src/sequence.jl new file mode 100644 index 00000000..bf56cd93 --- /dev/null +++ b/src/sequence.jl @@ -0,0 +1,48 @@ +function rand_re(machine::Automa.Machine) + out = IOBuffer() + node = machine.start + + while true + if node.state ∈ machine.final_states + (rand() ≤ 1 / (length(node.edges) + 1)) && break + end + + edge, node = rand(node.edges) + label = rand(collect(edge.labels)) + print(out, Char(label)) + end + + return String(take!(out)) +end + +sequencestring(d::SequenceDesign) = sequencestring(d.sequence) +function sequencestring(str::String) + #match curly brackets and replace them + @assert isnothing(findfirst("*", str)) && isnothing(findfirst("+", str)) "'infinite' sequences 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(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] + @info str + end + return rand_re(Automa.compile(RE(str))) +end \ No newline at end of file diff --git a/src/types.jl b/src/types.jl index fb29ef53..7154340c 100644 --- a/src/types.jl +++ b/src/types.jl @@ -13,7 +13,10 @@ abstract type AbstractHeadmodel end struct Simulation design::AbstractDesign - components::Vector{AbstractComponent} + components::Union{ + <:Dict{<:Char,<:Vector{<:AbstractComponent}}, + <:Vector{<:AbstractComponent}, + } onset::AbstractOnset noisetype::AbstractNoise end From 49db0ba0ebf53442c90aee76ff614b481c3b849a Mon Sep 17 00:00:00 2001 From: "behinger (s-ccs 001)" Date: Sun, 25 Feb 2024 21:50:24 +0000 Subject: [PATCH 03/37] fix bug in predef_eeg --- src/predefinedSimulations.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/predefinedSimulations.jl b/src/predefinedSimulations.jl index 874bb870..5a1c2a65 100644 --- a/src/predefinedSimulations.jl +++ b/src/predefinedSimulations.jl @@ -79,7 +79,7 @@ function predef_eeg( kwargs..., ) - components = [] + components = AbstractComponent[] for c in comps append!(components, [T(c...)]) end From cd489c4fef1479b4c613ba426c5e5bb13058999e Mon Sep 17 00:00:00 2001 From: "behinger (s-ccs 001)" Date: Wed, 28 Feb 2024 09:56:51 +0000 Subject: [PATCH 04/37] forgot the end --- src/onset.jl | 1 + 1 file changed, 1 insertion(+) diff --git a/src/onset.jl b/src/onset.jl index f2e2469c..a30ce539 100644 --- a/src/onset.jl +++ b/src/onset.jl @@ -93,6 +93,7 @@ function simulate_onsets(rng, onset::AbstractOnset, simulation::Simulation) 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?" From 6c35f3c903b9150df89547fad6410196f5706fc4 Mon Sep 17 00:00:00 2001 From: "behinger (s-ccs 001)" Date: Wed, 28 Feb 2024 13:41:49 +0000 Subject: [PATCH 05/37] half way through to success for sequence designs or something --- src/component.jl | 87 ++++++++++++++++++++++++++++++++++++++--------- src/design.jl | 15 ++++---- src/onset.jl | 4 +-- src/sequence.jl | 16 ++++----- src/simulation.jl | 27 +++++++++++++-- 5 files changed, 113 insertions(+), 36 deletions(-) diff --git a/src/component.jl b/src/component.jl index 75f1e54b..802a38aa 100644 --- a/src/component.jl +++ b/src/component.jl @@ -27,6 +27,7 @@ MixedModelComponent(; β::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 """ @@ -55,9 +56,26 @@ LinearModelComponent(; formula::Any # e.g. 0~1+cond - left side must be "0" β::Vector contrasts::Dict = Dict() + offset::Int = 0 end +""" + 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 Integer +""" +offset(c::AbstractComponent)::Int = 0 +offset(c::LinearModelComponent)::Int = c.offset +offset(c::MixedModelComponent)::Int = c.offset + +maxoffset(c::Vector{<:AbstractComponent}) = maximum(offset.(c)) +maxoffset(d::Dict{<:Char,<:Vector{<:AbstractComponent}}) = maximum(maxoffset.(values(d))) +minoffset(c::Vector{<:AbstractComponent}) = minimum(offset.(c)) +minoffset(d::Dict{<:Char,<:Vector{<:AbstractComponent}}) = minimum(minoffset.(values(d))) + + + """ Wrapper for an `AbstractComponent` to project it to multiple target-channels via `projection`. optional adds `noise` to the source prior to projection. """ @@ -94,6 +112,7 @@ For `MultichannelComponent` return the length of the projection vector. """ n_channels(c::MultichannelComponent) = length(c.projection) + """ For a vector of `MultichannelComponent`s, return the first but asserts all are of equal length. """ @@ -145,7 +164,7 @@ simulate_component(rng, c::AbstractComponent, simulation::Simulation) = simulate_component(rng, c, simulation.design) """ - simulate_component(rng, c::AbstractComponent, simulation::Simulation) + simulate_component(rng, c::LinearModelComponent, design::AbstractDesign) Generate a linear model design matrix, weight it by c.β and multiply the result with the given basis vector. julia> c = UnfoldSim.LinearModelComponent([0,1,1,0],@formula(0~1+cond),[1,2],Dict()) @@ -192,13 +211,13 @@ Currently, it is not possible to use a different basis for fixed and random effe julia> design = MultiSubjectDesign(;n_subjects=2,n_items=50,items_between=(;:cond=>["A","B"])) julia> c = UnfoldSim.MixedModelComponent([0.,1,1,0],@formula(0~1+cond+(1|subject)),[1,2],Dict(:subject=>[2],),Dict()) -julia> simulate(StableRNG(1),c,design) +julia> simulate_component(StableRNG(1),c,design) """ function simulate_component( rng, c::MixedModelComponent, - design::AbstractDesign; + design::AbstractDesign, return_parameters = false, ) events = generate_events(design) @@ -335,34 +354,51 @@ function simulate_responses!( return epoch_data end function init_epoch_data(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), length(design)) + epoch_data = zeros( + n_channels(components), + maxlength(components) + range_offset, + length(design), + ) else - epoch_data = zeros(maxlength(components), length(design)) + epoch_data = zeros(maxlength(components) + range_offset, length(design)) end return epoch_data end -function simulate_responses(rng, components::Dict, s::Simulation) - epoch_data = init_epoch_data(components, s.design) +function simulate_responses(rng, event_component_dict::Dict, s::Simulation) + epoch_data = init_epoch_data(event_component_dict, s.design) evts = generate_events(s.design) - multichannel = n_channels(components) > 1 - for key in keys(components) + 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), - components[key], + event_component_dict, s.onset, s.noisetype, ) ix = evts.event .== key if multichannel - simulate_responses!(rng, @view(epoch_data[:, :, ix]), components[key], s_key) + simulate_responses!( + rng, + @view(epoch_data[:, :, ix]), + event_component_dict[key], + s_key, + ) else - #@debug sum(ix), size(simulate_responses(rng, components[key], s_key)), key - simulate_responses!(rng, @view(epoch_data[:, ix]), components[key], s_key) + #@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 @@ -373,11 +409,28 @@ end simulate_and_add!(epoch_data::AbstractArray, c, simulation, rng) Helper function to call `simulate_component` and add it to a provided Array. """ -function simulate_and_add!(epoch_data::AbstractMatrix, c, simulation, rng) +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 = offset(component) - minoffset(simulation.components) + @debug off + @debug 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 = 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 5e826895..f9d86932 100644 --- a/src/design.jl +++ b/src/design.jl @@ -156,15 +156,18 @@ end @with_kw struct SequenceDesign{T} <: AbstractDesign design::T sequence::String = "" + sequencelength::Int = 0 + rng = nothing end +SequenceDesign(design, sequence) = SequenceDesign(design = design, sequence = sequence) -UnfoldSim.generate_events(design::SequenceDesign{MultiSubjectDesign}) = - error("not yet implemented") +generate_events(design::SequenceDesign{MultiSubjectDesign}) = error("not yet implemented") -function UnfoldSim.generate_events(design::SequenceDesign) + +function generate_events(design::SequenceDesign) df = generate_events(design.design) nrows_df = size(df, 1) - currentsequence = sequencestring(design.sequence) + currentsequence = sequencestring(deepcopy(design.rng), design.sequence) currentsequence = replace(currentsequence, "_" => "") df = repeat(df, inner = length(currentsequence)) df.event .= repeat(collect(currentsequence), nrows_df) @@ -175,11 +178,11 @@ end """ - UnfoldSim.generate_events(design::RepeatDesign{T}) + generate_events(rng,design::RepeatDesign{T}) In a repeated design, iteratively calls the underlying {T} Design and concatenates. In case of MultiSubjectDesign, sorts by subject. """ -function UnfoldSim.generate_events(design::RepeatDesign) +function generate_events(design::RepeatDesign) df = map(x -> generate_events(design.design), 1:design.repeat) |> x -> vcat(x...) if isa(design.design, MultiSubjectDesign) sort!(df, [:subject]) diff --git a/src/onset.jl b/src/onset.jl index a30ce539..ea607506 100644 --- a/src/onset.jl +++ b/src/onset.jl @@ -66,7 +66,7 @@ contains_design(d::AbstractDesign, target::Type) = false contains_design(d::Union{RepeatDesign,SequenceDesign,SubselectDesign}, target::Type) = d.design isa target ? true : contains_design(d.design, target) -sequencestring(d::RepeatDesign) = sequencestring(d.design) +sequencestring(rng, d::RepeatDesign) = sequencestring(rng, d.design) """ simulate_onsets(rng, onset::AbstractOnset, simulation::Simulation) @@ -80,7 +80,7 @@ function simulate_onsets(rng, onset::AbstractOnset, simulation::Simulation) if contains_design(simulation.design, SequenceDesign) - currentsequence = sequencestring(simulation.design) + currentsequence = sequencestring(deepcopy(rng), simulation.design) if !isnothing(findfirst("_", currentsequence)) @assert currentsequence[end] == '_' "the blank-indicator '_' has to be the last sequence element" diff --git a/src/sequence.jl b/src/sequence.jl index bf56cd93..79ff3021 100644 --- a/src/sequence.jl +++ b/src/sequence.jl @@ -1,4 +1,4 @@ -function rand_re(machine::Automa.Machine) +function rand_re(rng::AbstractRNG, machine::Automa.Machine) out = IOBuffer() node = machine.start @@ -7,16 +7,16 @@ function rand_re(machine::Automa.Machine) (rand() ≤ 1 / (length(node.edges) + 1)) && break end - edge, node = rand(node.edges) - label = rand(collect(edge.labels)) + edge, node = rand(rng, node.edges) + label = rand(rng, collect(edge.labels)) print(out, Char(label)) end return String(take!(out)) end -sequencestring(d::SequenceDesign) = sequencestring(d.sequence) -function sequencestring(str::String) +sequencestring(rng, d::SequenceDesign) = sequencestring(rng, d.sequence) +function sequencestring(rng, str::String) #match curly brackets and replace them @assert isnothing(findfirst("*", str)) && isnothing(findfirst("+", str)) "'infinite' sequences currently not supported" crly = collect(eachmatch(r"(\{[\d],[\d]\})", str)) @@ -36,13 +36,13 @@ function sequencestring(str::String) repeat_string = string(str[c.offset-1]) end - replacement_string = repeat(repeat_string, rand(rep_minimum:rep_maximum)) + 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] - @info str + #@info str end - return rand_re(Automa.compile(RE(str))) + return rand_re(rng, Automa.compile(RE(str))) end \ No newline at end of file diff --git a/src/simulation.jl b/src/simulation.jl index 64652a49..7a292330 100755 --- a/src/simulation.jl +++ b/src/simulation.jl @@ -47,6 +47,21 @@ simulate( kwargs..., ) = simulate(rng, Simulation(design, signal, onset, noise); kwargs...) +function simulate( + rng::AbstractRNG, + design::SequenceDesign, + signal, + onset::AbstractOnset, + noise::AbstractNoise = NoNoise(); + kwargs..., +) + design = + isnothing(design.rng) ? + SequenceDesign(design.design, design.sequence, design.sequencelength, rng) : design + simulate(rng, Simulation(design, signal, onset, noise); kwargs...) +end + + function simulate(rng::AbstractRNG, simulation::Simulation; return_epoched::Bool = false) (; design, components, onset, noisetype) = simulation @@ -127,11 +142,13 @@ function create_continuous_signal(rng, responses, simulation) # combine responses with onsets max_length_component = maxlength(components) - max_length_continuoustime = Int(ceil(maximum(onsets))) .+ max_length_component + offset_range = maxoffset(simulation.components) - minoffset(simulation.components) + max_length_continuoustime = + Int(ceil(maximum(onsets))) .+ max_length_component .+ offset_range signal = zeros(n_chan, max_length_continuoustime, n_subjects) - + @debug size(signal), offset_range for e = 1:n_chan for s = 1:n_subjects for i = 1:n_trials @@ -141,7 +158,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 @@ -172,6 +191,8 @@ function add_responses!(signal, responses::Vector, e, s, tvec, erpvec) @views signal[e, tvec, s] .+= responses[:, erpvec] end function add_responses!(signal, responses::Matrix, e, s, tvec, erpvec)# + @debug size(signal), size(responses), e, s, size(tvec), size(erpvec) + @debug tvec, erpvec @views signal[e, tvec, s] .+= responses[:, erpvec] end function add_responses!(signal, responses::AbstractArray, e, s, tvec, erpvec) From 809dc2e9bf6c6abf1fd0305f2689d7cc55b3ceb2 Mon Sep 17 00:00:00 2001 From: "behinger (s-ccs 001)" Date: Wed, 28 Feb 2024 15:13:12 +0000 Subject: [PATCH 06/37] everythinig except sequencelength seems to be working now --- docs/literate/HowTo/sequence.jl | 27 +++++++++++++------ src/component.jl | 8 ++++-- src/design.jl | 46 +++++++++++++++++++++++++++------ src/onset.jl | 1 - src/sequence.jl | 9 ++++--- src/simulation.jl | 36 +++++++++++++++----------- 6 files changed, 90 insertions(+), 37 deletions(-) diff --git a/docs/literate/HowTo/sequence.jl b/docs/literate/HowTo/sequence.jl index e66e5e95..d9da5c4f 100644 --- a/docs/literate/HowTo/sequence.jl +++ b/docs/literate/HowTo/sequence.jl @@ -6,20 +6,20 @@ using StableRNGs # 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"])) +design = SingleSubjectDesign(conditions = Dict(:condition => ["one", "two"]))#|>x->RepeatDesign(x,4) generate_events(design) # next we use the `SequenceDesign` and nest our initial design in it. "SR_" is code for an "S" event and an "R" event - only single letter events are supported! The `_` is a signal for the Onset generator to generate a bigger pause - no overlap between adjacend `SR` pairs -design = SequenceDesign(design, "SRR_") +design = SequenceDesign(design, "SR{1,2}_", 0, 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 exist, like "SR{1,3}", or "A[BC]" etc. # Finally, let's repeat the design 2 times -design = RepeatDesign(design, 2) +design = RepeatDesign(design, 4) generate_events(design) -# This results in 20 trials that nicely follow our sequence +# This results in 12 trials that nicely follow our sequence # Next we have to specify for both events `S` and `R` what the responses should look like. @@ -37,16 +37,27 @@ n1 = LinearModelComponent(; p3 = LinearModelComponent(; basis = UnfoldSim.hanning(Int(0.5 * 100)), # sfreq = 100 for the other bases formula = @formula(0 ~ 1 + condition), - β = [1, 2], + β = [1, 0], ); -components = Dict('S' => [p1, n1], 'R' => [p3]) +resp = LinearModelComponent(; + basis = UnfoldSim.hanning(Int(0.5 * 100)), # sfreq = 100 for the other bases + formula = @formula(0 ~ 1 + condition), + β = [1, 2], + offset = -10, +); +components = Dict('S' => [p1, n1, p3], 'R' => [resp]) +#components = [p1, n1, resp] data, evts = simulate( StableRNG(1), design, components, - UniformOnset(offset = 10, width = 100), + UniformOnset(offset = 40, width = 10), NoNoise(), ) -lines(data) \ No newline at end of file + +lines(data) +vlines!(evts.latency, color = (:gray, 0.5)) +xlims!(0, 500) +current_figure() \ No newline at end of file diff --git a/src/component.jl b/src/component.jl index 802a38aa..d3a0c5e7 100644 --- a/src/component.jl +++ b/src/component.jl @@ -370,8 +370,12 @@ function init_epoch_data(components, design) end function simulate_responses(rng, event_component_dict::Dict, s::Simulation) + #@debug rng.state epoch_data = init_epoch_data(event_component_dict, s.design) + #@debug rng.state evts = generate_events(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 == '_' @@ -418,8 +422,8 @@ function simulate_and_add!( @debug "matrix" off = offset(component) - minoffset(simulation.components) - @debug off - @debug 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 f9d86932..5e315fe9 100644 --- a/src/design.jl +++ b/src/design.jl @@ -153,23 +153,45 @@ design = RepeatDesign(designOnce,4); 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" + return s +end @with_kw struct SequenceDesign{T} <: AbstractDesign design::T sequence::String = "" sequencelength::Int = 0 rng = nothing + + SequenceDesign{T}(d, s, sl, r) where {T<:AbstractDesign} = + new(d, check_sequence(s), sl, r) end + SequenceDesign(design, sequence) = SequenceDesign(design = design, sequence = sequence) generate_events(design::SequenceDesign{MultiSubjectDesign}) = error("not yet implemented") -function generate_events(design::SequenceDesign) +generate_events(rng, design::AbstractDesign) = generate_events(design) +generate_events(design::SequenceDesign) = generate_events(deepcopy(design.rng), design) + +function generate_events(rng, design::SequenceDesign) df = generate_events(design.design) nrows_df = size(df, 1) - currentsequence = sequencestring(deepcopy(design.rng), design.sequence) + + rng = if isnothing(rng) + @warn "Could not (yet) find an rng for `SequenceDesign` - ignore this message if you called `generate_events` yourself, be worried if you called `simulate` and still see this message. Surpress this message by defining the `rng` when creating the `SequenceDesign`" + MersenneTwister(1) + else + rng + end + # @debug design.sequence + currentsequence = sequencestring(rng, design.sequence) + # @debug currentsequence currentsequence = replace(currentsequence, "_" => "") df = repeat(df, inner = length(currentsequence)) + df.event .= repeat(collect(currentsequence), nrows_df) return df @@ -177,13 +199,19 @@ function generate_events(design::SequenceDesign) end +get_rng(design::AbstractDesign) = nothing +get_rng(design::SequenceDesign) = design.rng + """ generate_events(rng,design::RepeatDesign{T}) In a repeated design, iteratively calls the underlying {T} Design and concatenates. In case of MultiSubjectDesign, sorts by subject. """ function generate_events(design::RepeatDesign) - df = map(x -> generate_events(design.design), 1:design.repeat) |> x -> vcat(x...) + design = deepcopy(design) + df = + map(x -> generate_events(get_rng(design.design), design.design), 1:design.repeat) |> + x -> vcat(x...) if isa(design.design, MultiSubjectDesign) sort!(df, [:subject]) end @@ -208,10 +236,12 @@ end Base.size(design::RepeatDesign{MultiSubjectDesign}) = size(design.design) .* (design.repeat, 1) Base.size(design::RepeatDesign{SingleSubjectDesign}) = size(design.design) .* design.repeat -Base.size(design::SequenceDesign) = - size(design.design) .* length(replace(design.sequence, "_" => "")) +#Base.size(design::SequenceDesign) = +#size(design.design) .* length(replace(design.sequence, "_" => "",r"\{.*\}"=>"")) -Base.size(design::RepeatDesign{SequenceDesign{SingleSubjectDesign}}) = - size(design.design) .* design.repeat +#Base.size(design::) = size(design.design) .* design.repeat -Base.size(design::SubselectDesign) = size(generate_events(design), 1) \ No newline at end of file +# No way to find out what size it is without actually generating first... +Base.size( + design::Union{<:SequenceDesign,<:SubselectDesign,<:RepeatDesign{<:SequenceDesign}}, +) = size(generate_events(design), 1) \ No newline at end of file diff --git a/src/onset.jl b/src/onset.jl index ea607506..fc921d75 100644 --- a/src/onset.jl +++ b/src/onset.jl @@ -66,7 +66,6 @@ contains_design(d::AbstractDesign, target::Type) = false contains_design(d::Union{RepeatDesign,SequenceDesign,SubselectDesign}, target::Type) = d.design isa target ? true : contains_design(d.design, target) -sequencestring(rng, d::RepeatDesign) = sequencestring(rng, d.design) """ simulate_onsets(rng, onset::AbstractOnset, simulation::Simulation) diff --git a/src/sequence.jl b/src/sequence.jl index 79ff3021..8322455f 100644 --- a/src/sequence.jl +++ b/src/sequence.jl @@ -16,10 +16,11 @@ function rand_re(rng::AbstractRNG, machine::Automa.Machine) end sequencestring(rng, d::SequenceDesign) = sequencestring(rng, d.sequence) + function sequencestring(rng, str::String) #match curly brackets and replace them @assert isnothing(findfirst("*", str)) && isnothing(findfirst("+", str)) "'infinite' sequences currently not supported" - crly = collect(eachmatch(r"(\{[\d],[\d]\})", str)) + crly = collect(eachmatch(r"(\{[\d]+,[\d]+\})", str)) for c in reverse(crly) m = replace(c.match, "{" => "", "}" => "") rep_minimum, rep_maximum = parse.(Int, split(m, ",")) @@ -42,7 +43,9 @@ function sequencestring(rng, str::String) str[1:bracket_start_idx-1] * replacement_string * str[bracket_end_idx+length(c.match)+1:end] - #@info str + # @debug str end return rand_re(rng, Automa.compile(RE(str))) -end \ No newline at end of file +end + +sequencestring(rng, d::RepeatDesign) = sequencestring(rng, d.design) diff --git a/src/simulation.jl b/src/simulation.jl index 7a292330..4392424c 100755 --- a/src/simulation.jl +++ b/src/simulation.jl @@ -29,6 +29,8 @@ Some remarks to how the noise is added: - If `return_epoched = true` and `onset =NoOnset()` the noise is added to the epoched data matrix - If `onset` is not `NoOnset`, a continuous signal is created and the noise is added to this i.e. this means that the noise won't be the same as in the `onset = NoOnset()` case even if `return_epoched = true`. - The case `return_epoched = false` and `onset = NoOnset()` is not possible and therefore covered by an assert statement + - `simulate(rng,design::SequenceDesign,...)` + If no `design.rng` was defined for `SequenceDesign`, we replace it with the `simulation`-function call `rng` object """ @@ -38,30 +40,34 @@ function simulate(args...; kwargs...) simulate(MersenneTwister(1), args...; kwargs...) end -simulate( - rng::AbstractRNG, - design::AbstractDesign, - signal, - onset::AbstractOnset, - noise::AbstractNoise = NoNoise(); - kwargs..., -) = simulate(rng, Simulation(design, signal, onset, noise); kwargs...) - function simulate( rng::AbstractRNG, - design::SequenceDesign, + design::AbstractDesign, signal, onset::AbstractOnset, noise::AbstractNoise = NoNoise(); kwargs..., ) - design = - isnothing(design.rng) ? - SequenceDesign(design.design, design.sequence, design.sequencelength, rng) : design + + if is_SequenceDesign(design) + design = sequencedesign_add_rng(rng, design) + end simulate(rng, Simulation(design, signal, onset, noise); kwargs...) end +sequencedesign_add_rng(rng, design::AbstractDesign) = design +sequencedesign_add_rng(rng, design::RepeatDesign) = + RepeatDesign(sequencedesign_add_rng(rng, design.design), design.repeat) +sequencedesign_add_rng(rng, design::SequenceDesign) = + isnothing(design.rng) ? + SequenceDesign(design.design, design.sequence, design.sequencelength, rng) : design + + +is_SequenceDesign(d::AbstractDesign) = false +is_SequenceDesign(d::RepeatDesign) = is_SequenceDesign(d.design) +is_SequenceDesign(d::SequenceDesign) = true + function simulate(rng::AbstractRNG, simulation::Simulation; return_epoched::Bool = false) (; design, components, onset, noisetype) = simulation @@ -191,8 +197,8 @@ function add_responses!(signal, responses::Vector, e, s, tvec, erpvec) @views signal[e, tvec, s] .+= responses[:, erpvec] end function add_responses!(signal, responses::Matrix, e, s, tvec, erpvec)# - @debug size(signal), size(responses), e, s, size(tvec), size(erpvec) - @debug tvec, erpvec + # @debug size(signal), size(responses), e, s, size(tvec), size(erpvec) + #@debug tvec, erpvec @views signal[e, tvec, s] .+= responses[:, erpvec] end function add_responses!(signal, responses::AbstractArray, e, s, tvec, erpvec) From 31c69e6611ec30eeec020fad724396d97890c1b0 Mon Sep 17 00:00:00 2001 From: "behinger (s-ccs 001)" Date: Wed, 28 Feb 2024 15:13:52 +0000 Subject: [PATCH 07/37] added string sequence tests --- test/runtests.jl | 1 + test/sequence.jl | 11 +++++++++++ 2 files changed, 12 insertions(+) create mode 100644 test/sequence.jl diff --git a/test/runtests.jl b/test/runtests.jl index 723fc891..2cbc1de5 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -8,4 +8,5 @@ include("setup.jl") include("onset.jl") include("simulation.jl") include("helper.jl") + include("sequence.jl") end diff --git a/test/sequence.jl b/test/sequence.jl new file mode 100644 index 00000000..e4a9ace0 --- /dev/null +++ b/test/sequence.jl @@ -0,0 +1,11 @@ +@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") + +end +@test length(UnfoldSim.sequencestring(StableRNG(1), "A{10,10}")) == 10 +@test length(UnfoldSim.sequencestring(StableRNG(1), "A{10,10}B")) == 11 +@test length(UnfoldSim.sequencestring(StableRNG(1), "A{10,20}")) >= 10 \ No newline at end of file From 54e933490464f99a915f8ed9149026445d5c5fbb Mon Sep 17 00:00:00 2001 From: "behinger (s-ccs 001)" Date: Fri, 1 Mar 2024 18:59:45 +0000 Subject: [PATCH 08/37] small doc update --- docs/literate/HowTo/sequence.jl | 9 +++++++-- src/design.jl | 1 + 2 files changed, 8 insertions(+), 2 deletions(-) diff --git a/docs/literate/HowTo/sequence.jl b/docs/literate/HowTo/sequence.jl index d9da5c4f..6a07d15f 100644 --- a/docs/literate/HowTo/sequence.jl +++ b/docs/literate/HowTo/sequence.jl @@ -1,3 +1,4 @@ +using Base: add_sum using UnfoldSim using CairoMakie using StableRNGs @@ -13,12 +14,16 @@ design = SequenceDesign(design, "SR{1,2}_", 0, 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 exist, like "SR{1,3}", or "A[BC]" etc. +# 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 design 2 times +# Finally, let's repeat the design 2 times - because we can design = RepeatDesign(design, 4) generate_events(design) +#design = UnfoldSim.AddSaccadeAmplitudeDesign4(design,:rt,Normal(0,1),MersenneTwister(1)) +#generate_events(design) + + # This results in 12 trials that nicely follow our sequence # Next we have to specify for both events `S` and `R` what the responses should look like. diff --git a/src/design.jl b/src/design.jl index 5e315fe9..7779c423 100644 --- a/src/design.jl +++ b/src/design.jl @@ -199,6 +199,7 @@ function generate_events(rng, design::SequenceDesign) end + get_rng(design::AbstractDesign) = nothing get_rng(design::SequenceDesign) = design.rng From 7114cd65e5017db7c0bfe0ba55cbd96d3b54a71a Mon Sep 17 00:00:00 2001 From: "behinger (s-ccs 001)" Date: Sat, 9 Mar 2024 14:05:06 +0000 Subject: [PATCH 09/37] added jitter to the '_' trial divisor --- src/onset.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/onset.jl b/src/onset.jl index fc921d75..a6989c26 100644 --- a/src/onset.jl +++ b/src/onset.jl @@ -89,7 +89,7 @@ function simulate_onsets(rng, onset::AbstractOnset, simulation::Simulation) # 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) + onsets[stepsize+1:stepsize:end] .+= 2 .* maxlength(simulation.components) #@debug onsets[stepsize:stepsize:end] end end From 85f037b3f04745a77ebcfceb826defc749d794f0 Mon Sep 17 00:00:00 2001 From: jschepers Date: Wed, 10 Jul 2024 16:15:44 +0000 Subject: [PATCH 10/37] Improve documentation especially quickstart, home and power analysis --- docs/literate/tutorials/poweranalysis.jl | 2 +- docs/literate/tutorials/quickstart.jl | 23 ++++++++++++++--------- docs/make.jl | 2 +- docs/src/index.md | 8 ++++---- src/noise.jl | 2 +- 5 files changed, 21 insertions(+), 16 deletions(-) diff --git a/docs/literate/tutorials/poweranalysis.jl b/docs/literate/tutorials/poweranalysis.jl index bbe7f1a5..5202ad66 100644 --- a/docs/literate/tutorials/poweranalysis.jl +++ b/docs/literate/tutorials/poweranalysis.jl @@ -5,7 +5,7 @@ using DataFrames using Random -# ## Simple Poweranalysis Script +# ## Simple Power analysis Script # For a power analysis, we will repeatedly simulate data, and check whether we can find a significant effect. # # We perform the power analysis on epoched data. diff --git a/docs/literate/tutorials/quickstart.jl b/docs/literate/tutorials/quickstart.jl index 86896185..5811931b 100644 --- a/docs/literate/tutorials/quickstart.jl +++ b/docs/literate/tutorials/quickstart.jl @@ -1,18 +1,23 @@ +# # Quickstart + using UnfoldSim -using Random -using CairoMakie +using Random # to get an RNG +using CairoMakie # for plotting +# To get started with data simulation, the user needs to provide four ingredients: an experimental design, an event basis function (component), an inter-onset distribution and a noise specification. # !!! tip # Use `subtypes(AbstractNoise)` (or `subtypes(AbstractComponent)` etc.) to find already implemented building blocks. -# ## "Experimental" Design +# ## Specify the simulation ingredients + +# ### Experimental Design # Define a 1 x 2 design with 20 trials. That is, one condition (`condaA`) with two levels. design = SingleSubjectDesign(; conditions = Dict(:condA => ["levelA", "levelB"])) |> x -> RepeatDesign(x, 10); -# #### Component / Signal +# ### Event basis function (Component) # Define a simple component and ground truth simulation formula. Akin to ERP components, we call one simulation signal a component. # # !!! note @@ -23,19 +28,19 @@ signal = LinearModelComponent(; β = [1, 0.5], ); -# #### Onsets and Noise -# We will start with a uniform (but overlapping, `offset` < `length(signal.basis)`) onset-distribution +# ### Onsets and Noise +# We will start with a uniform (but overlapping, `offset` < `length(signal.basis)`) inter-onset distribution. onset = UniformOnset(; width = 20, offset = 4); # And we will use some noise noise = PinkNoise(; noiselevel = 0.2); # ## Combine & Generate -# Finally, we will simulate some data +# Finally, we will combine all ingredients and simulate some data. data, events = simulate(MersenneTwister(1), design, signal, onset, noise); -# `Data` is a `n-sample` Vector (but could be a Matrix for e.g. `MultiSubjectDesign`). +# `data` is a `n-sample` Vector (but could be a Matrix for e.g. `MultiSubjectDesign` or epoched data). -# `Events` is a DataFrame that contains a column `latency` with the onsets of events. +# `events` is a DataFrame that contains a column `latency` with the onsets of events (in samples). # ## Plot them! lines(data; color = "black") diff --git a/docs/make.jl b/docs/make.jl index ffba630f..3403cbd7 100644 --- a/docs/make.jl +++ b/docs/make.jl @@ -34,7 +34,7 @@ makedocs(; "Tutorials" => [ "Quickstart" => "generated/tutorials/quickstart.md", "Simulate ERPs" => "generated/tutorials/simulateERP.md", - "Poweranalysis" => "generated/tutorials/poweranalysis.md", + "Power analysis" => "generated/tutorials/poweranalysis.md", "Multi-subject simulation" => "generated/tutorials/multisubject.md", ], "Reference" => [ diff --git a/docs/src/index.md b/docs/src/index.md index 415199d4..40ff2f14 100644 --- a/docs/src/index.md +++ b/docs/src/index.md @@ -2,14 +2,14 @@ CurrentModule = UnfoldSim ``` -# UnfoldSim +# UnfoldSim.jl -Documentation for [UnfoldSim](https://github.com/behinger/UnfoldSim.jl). +Documentation for [UnfoldSim.jl](https://github.com/unfoldtoolbox/UnfoldSim.jl). UnfoldSim.jl is a Julia package for simulating multivariate timeseries data with a special focus on EEG data. ## Start simulating timeseries -We offer some predefined signals, check them out! +We offer some predefined (EEG) signals, check them out! -For instance an P1/N170/P300 complex. +For instance a P1/N170/P300 complex (containing three typical ERP components). ```@example main using UnfoldSim using CairoMakie # plotting diff --git a/src/noise.jl b/src/noise.jl index 6b7b8592..28421ccd 100644 --- a/src/noise.jl +++ b/src/noise.jl @@ -43,7 +43,7 @@ end """ RealisticNoise <: AbstractNoise -Not implemented - planned to use Artefacts.jl to provide real EEG data to add. +Not implemented - planned to use Artifacts.jl to provide real EEG data to add. """ @with_kw struct RealisticNoise <: AbstractNoise noiselevel = 1 From 2c65a7fa69605d018bca7d02980616ba1f3a123f Mon Sep 17 00:00:00 2001 From: jschepers Date: Thu, 18 Jul 2024 06:09:35 +0000 Subject: [PATCH 11/37] adapted the order of reference overviews and adapted titles --- docs/literate/reference/basistypes.jl | 4 ++-- docs/literate/reference/noisetypes.jl | 5 +++-- docs/literate/reference/onsettypes.jl | 12 ++++++------ docs/make.jl | 8 ++++---- 4 files changed, 15 insertions(+), 14 deletions(-) diff --git a/docs/literate/reference/basistypes.jl b/docs/literate/reference/basistypes.jl index c04d73cc..de4db4a6 100644 --- a/docs/literate/reference/basistypes.jl +++ b/docs/literate/reference/basistypes.jl @@ -3,14 +3,14 @@ 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. # ## 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] diff --git a/docs/literate/reference/noisetypes.jl b/docs/literate/reference/noisetypes.jl index fe3efc4a..26af0e1e 100644 --- a/docs/literate/reference/noisetypes.jl +++ b/docs/literate/reference/noisetypes.jl @@ -1,10 +1,11 @@ +# # Overview: Noise types +# There are several noise types directly implemented. Here is a comparison: + using UnfoldSim using CairoMakie using DSP using StableRNGs import StatsBase.autocor -# ## What's the noise? -# There are several noise-types directly implemented. Here is a comparison: f = Figure() ax_sig = diff --git a/docs/literate/reference/onsettypes.jl b/docs/literate/reference/onsettypes.jl index 279a2226..0c0980b1 100644 --- a/docs/literate/reference/onsettypes.jl +++ b/docs/literate/reference/onsettypes.jl @@ -1,4 +1,4 @@ -# # Onset types +# # 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`. @@ -58,11 +58,11 @@ let # hide design, # hide ) # hide - hist!( - ax, - distances, - bins = range(0, 100, step = 1), - label = "($width, $offset)", + hist!( # hide + ax, # hide + distances, # hide + bins = range(0, 100, step = 1), # hide + label = "($width, $offset)", # hide ) # hide if label == "offset" && offset != 0 # hide diff --git a/docs/make.jl b/docs/make.jl index 3403cbd7..efe43840 100644 --- a/docs/make.jl +++ b/docs/make.jl @@ -38,10 +38,10 @@ makedocs(; "Multi-subject simulation" => "generated/tutorials/multisubject.md", ], "Reference" => [ - "Overview: Toolbox Functions" => "./generated/reference/overview.md", - "Overview: NoiseTypes" => "./generated/reference/noisetypes.md", - "Overview: OnsetTypes" => "./generated/reference/onsettypes.md", - "Overview: Components (EEG, fMRI, Pupil)" => "./generated/reference/basistypes.md", + "Overview of functionality" => "./generated/reference/overview.md", + "Overview: Basis function (component) types" => "./generated/reference/basistypes.md", + "Overview: Onset types" => "./generated/reference/onsettypes.md", + "Overview: Noise types" => "./generated/reference/noisetypes.md", ], "HowTo" => [ "Define a new, (imbalanced) design" => "./generated/HowTo/newDesign.md", From ee9e00ad8b4565df05fdb4c9b67e3f5a68e87e92 Mon Sep 17 00:00:00 2001 From: jschepers Date: Mon, 22 Jul 2024 13:13:16 +0000 Subject: [PATCH 12/37] Updated quickstart page --- docs/literate/tutorials/quickstart.jl | 31 ++++++++++++++++++++------- 1 file changed, 23 insertions(+), 8 deletions(-) diff --git a/docs/literate/tutorials/quickstart.jl b/docs/literate/tutorials/quickstart.jl index 5811931b..4d93bb68 100644 --- a/docs/literate/tutorials/quickstart.jl +++ b/docs/literate/tutorials/quickstart.jl @@ -1,9 +1,5 @@ # # Quickstart -using UnfoldSim -using Random # to get an RNG -using CairoMakie # for plotting - # To get started with data simulation, the user needs to provide four ingredients: an experimental design, an event basis function (component), an inter-onset distribution and a noise specification. # !!! tip @@ -11,10 +7,24 @@ using CairoMakie # for plotting # ## Specify the simulation ingredients +# ### Setup +# ```@raw html +#
+# Click to expand +# ``` +## Load required packages +using UnfoldSim +using Random # to get an RNG +using CairoMakie # for plotting + +# ```@raw html +#
+# ``` + # ### Experimental Design -# Define a 1 x 2 design with 20 trials. That is, one condition (`condaA`) with two levels. +# Define a 1 x 2 design with 20 trials. That is, one condition (`cond_A`) with two levels. design = - SingleSubjectDesign(; conditions = Dict(:condA => ["levelA", "levelB"])) |> + SingleSubjectDesign(; conditions = Dict(:cond_A => ["level_A", "level_B"])) |> x -> RepeatDesign(x, 10); # ### Event basis function (Component) @@ -24,7 +34,7 @@ design = # You could easily specify multiple components by providing a vector of components, which are automatically added at the same onsets. This procedure simplifies to generate some response that is independent of simulated condition, whereas other depends on it. signal = LinearModelComponent(; basis = [0, 0, 0, 0.5, 1, 1, 0.5, 0, 0], - formula = @formula(0 ~ 1 + condA), + formula = @formula(0 ~ 1 + cond_A), β = [1, 0.5], ); @@ -44,5 +54,10 @@ data, events = simulate(MersenneTwister(1), design, signal, onset, noise); # ## Plot them! lines(data; color = "black") -vlines!(events.latency; color = ["orange", "teal"][1 .+ (events.condA.=="levelB")]) +vlines!(events.latency; color = ["orange", "teal"][1 .+ (events.cond_A.=="level_B")]) + +current_axis().title = "Simulated data" +current_axis().xlabel = "Time [samples]" +current_axis().ylabel = "Amplitude [μV]" + current_figure() From ac0d3bd736e0393c0855e1855a992fe225e8f94d Mon Sep 17 00:00:00 2001 From: jschepers Date: Mon, 22 Jul 2024 13:13:43 +0000 Subject: [PATCH 13/37] minor changes --- README.md | 2 +- docs/literate/reference/onsettypes.jl | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/README.md b/README.md index cdfb3b45..cac86db8 100644 --- a/README.md +++ b/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. diff --git a/docs/literate/reference/onsettypes.jl b/docs/literate/reference/onsettypes.jl index 0c0980b1..bbadef15 100644 --- a/docs/literate/reference/onsettypes.jl +++ b/docs/literate/reference/onsettypes.jl @@ -7,7 +7,7 @@ #
# Click to expand # ``` - +## Load required packages using UnfoldSim using CairoMakie using Random From 373e535b3d09e93cd5d50b7867bfc5490897ef10 Mon Sep 17 00:00:00 2001 From: jschepers Date: Tue, 23 Jul 2024 16:18:16 +0000 Subject: [PATCH 14/37] fixed docstrings for predef_eeg and predef_2x2 --- src/predefinedSimulations.jl | 83 ++++++++++++++++++------------------ 1 file changed, 41 insertions(+), 42 deletions(-) diff --git a/src/predefinedSimulations.jl b/src/predefinedSimulations.jl index 874bb870..8adf21f0 100644 --- a/src/predefinedSimulations.jl +++ b/src/predefinedSimulations.jl @@ -7,36 +7,36 @@ predef_eeg(; kwargs...) = predef_eeg(MersenneTwister(1); kwargs...) # without rn predef_eeg(nsubjects::Int; kwargs...) = predef_eeg(MersenneTwister(1), nsubjects; kwargs...) # without rng always call same one """ -predef_eeg(;kwargs...) -predef_eeg(rng;kwargs...) -predef_eeg(rng,n_subjects;kwargs...) + predef_eeg(; kwargs...) + predef_eeg(rng; kwargs...) + predef_eeg(rng, n_subjects; kwargs...) Generate a P1/N1/P3 complex. In case `n_subjects` is defined - `MixedModelComponents` are generated, else `LinearModelComponents`. -The most used `kwargs` is: `return_epoched=true` which returns already epoched data. If you want epoched data without overlap, specify `onset=NoOnset()` and `return_epoched=true` +The most used `kwargs` is: `return_epoched=true` which returns already epoched data. If you want epoched data without overlap, specify `onset = NoOnset()` and `return_epoched = true` +## Default parameters: -## Default params: +#### Design +- `n_repeats = 100`, +- `event_order_function = x -> shuffle(deepcopy(rng), x)`, # random trial order +- `conditions = Dict(...)`, -- n_repeats = 100 -- event_order_function = x->shuffle(deepcopy(rng),x # random trial order -- conditions = Dict(...), +#### Component / Signal +- `sfreq = 100`, +- `p1 = (p100(; sfreq = sfreq), @formula(0 ~ 1), [5], Dict())`, # P1 amp 5, no effects +- `n1 = (n170(; sfreq = sfreq), @formula(0 ~ 1 + condition), [5,-3], Dict())`, # N1 amp 5, dummy-coded condition effect (levels "car", "face") of -3 +- `p3 = (p300(; sfreq = sfreq), @formula(0 ~ 1 + continuous), [5,1], Dict())`, # P3 amp 5, continuous effect range [-5,5] with slope 1 -#### component / signal -- sfreq = 100, -- p1 = (p100(;sfreq=sfreq), @formula(0~1),[5],Dict()), # P1 amp 5, no effects -- n1 = (n170(;sfreq=sfreq), @formula(0~1+condition),[5,-3],Dict()), # N1 amp 5, dummycoded condition effect (levels "car", "face") of -3 -- p3 = (p300(;sfreq=sfreq), @formula(0~1+continuous),[5,1],Dict()), # P3 amp 5, continuous effect range [-5,5] with slope 1 +#### Onset +- `overlap = (0.5,0.2)`, # offset + width/length of Uniform noise. put offset to 1 for no overlap. put width to 0 for no jitter +- `onset = UniformOnset(; offset = sfreq * 0.5 * overlap[1], width = sfreq * 0.5 * overlap[2])`, -#### noise -- noiselevel = 0.2, -- noise = PinkNoise(;noiselevel=noiselevel), - -#### onset -- overlap = (0.5,0.2), # offset + width/length of Uniform noise. put offset to 1 for no overlap. put width to 0 for no jitter -- onset=UniformOnset(;offset=sfreq*0.5*overlap[1],width=sfreq*0.5*overlap[2]), +#### Noise +- `noiselevel = 0.2`, +- `noise = PinkNoise(; noiselevel = noiselevel)`, """ function predef_eeg( rng; @@ -139,34 +139,33 @@ function predef_eeg( end """ - predef_2x2(rng::AbstractRNG;kwargs...) + predef_2x2(rng::AbstractRNG; kwargs...) -The most used `kwargs` is: `return_epoched=true` which returns already epoched data. If you want epoched data without overlap, specify `onset=NoOnset()` and `return_epoched=true` +The most used `kwargs` is: `return_epoched = true` which returns already epoched data. If you want epoched data without overlap, specify `onset = NoOnset()` and `return_epoched = true` -#### design -- `n_items`=100, -- `n_subjects`=1, -- `conditions` = Dict(:A=>["a_small","a_big"],:B=>["b_tiny","b_large"]), -- `event_order_function` = x->shuffle(deepcopy(rng),x), +#### Design +- `n_items = 100`, +- `n_subjects = 1`, +- `conditions = Dict(:A => ["a_small","a_big"], :B => ["b_tiny","b_large"])`, +- `event_order_function = x -> shuffle(deepcopy(rng), x)`, -#### component / signal -- `signalsize` = 100, length of simulated hanning window -- `basis`` = hanning(signalsize), the actual "function", `signalsize` is only used here -- `β` = [1,-0.5,.5,+1], the parameters -- `σs` = Dict(:subject=>[1,0.5,0.5,0.5],:item=>[1]), - only in n_subjects>=2 case, specifies the random effects -- `contrasts` = Dict(:A=>EffectsCoding(),:B=>EffectsCoding()) - effect coding by default -- `formula` = n_subjects==1 ? @formula(0~1+A*B) : @formula(dv~1+A*B+(A*B|subject)+(1|item)), +#### Component / Signal +- `signalsize = 100`, # Length of simulated hanning window +- `basis = hanning(signalsize)`, # The actual "function", `signalsize` is only used here +- `β = [1, -0.5, .5, +1]`, # The parameters +- `σs = Dict(:subject => [1, 0.5, 0.5, 0.5],:item => [1])`, # Only in n_subjects >= 2 case, specifies the random effects +- `contrasts = Dict(:A => EffectsCoding(), :B => EffectsCoding())` # Effect coding by default +- `formula = n_subjects == 1 ? @formula(0 ~ 1 + A*B) : @formula(dv ~ 1 + A*B + (A*B|subject) + (1|item))`, -#### noise -- `noiselevel` = 0.2, -- `noise` = PinkNoise(;noiselevel=noiselevel), +#### Onset +- `overlap = (0.5,0.2)`, +- `onset = UniformOnset(; offset = signalsize * overlap[1], width = signalsize * overlap[2])`, # Put offset to 1 for no overlap. put width to 0 for no jitter -#### onset -- `overlap` = (0.5,0.2), -- `onset`=UniformOnset(;offset=signalsize*overlap[1],width=signalsize*overlap[2]), #put offset to 1 for no overlap. put width to 0 for no jitter +#### Noise +- `noiselevel = 0.2`, +- `noise = PinkNoise(; noiselevel = noiselevel)`, - -Careful if you modify n_items with n_subjects = 1, n_items has to be a multiple of 4 (or your equivalent conditions factorial, e.g. all combinations length) +Be careful if you modify n_items with n_subjects = 1, n_items has to be a multiple of 4 (or your equivalent conditions factorial, e.g. all combinations length) """ function predef_2x2( rng::AbstractRNG; From b073526492785ef1c0dda358a1e36b498b16a2ec Mon Sep 17 00:00:00 2001 From: jschepers Date: Tue, 23 Jul 2024 18:28:58 +0000 Subject: [PATCH 15/37] added draft of design types reference page --- docs/literate/HowTo/repeatTrials.jl | 2 +- docs/literate/reference/designtypes.jl | 100 +++++++++++++++++++++++++ docs/make.jl | 1 + 3 files changed, 102 insertions(+), 1 deletion(-) create mode 100644 docs/literate/reference/designtypes.jl diff --git a/docs/literate/HowTo/repeatTrials.jl b/docs/literate/HowTo/repeatTrials.jl index 046bfb75..1add60e6 100644 --- a/docs/literate/HowTo/repeatTrials.jl +++ b/docs/literate/HowTo/repeatTrials.jl @@ -1,7 +1,7 @@ using UnfoldSim -# ## Repeating Design entries +# # [Repeating Design entries](@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. diff --git a/docs/literate/reference/designtypes.jl b/docs/literate/reference/designtypes.jl new file mode 100644 index 00000000..d21693bf --- /dev/null +++ b/docs/literate/reference/designtypes.jl @@ -0,0 +1,100 @@ +# # Overview: Experimental design types + +# The experimental design specifies the experimental conditions and other variables that are supposed to have an influence on the simulated data. +# Currently, there are three types of designs implemented: `SingleSubjectDesign`, `MultiSubjectDesign` and `RepeatDesign`. + +# ## Setup +# ```@raw html +#
+# Click to expand +# ``` +## Load required packages +using UnfoldSim +using Random +# ```@raw html +#
+# ``` + +# ## Single-subject designs +# As the name suggests, the `SingleSubjectDesign` type can be used to specify the experimental design for a single subject. Using the `conditions` arguments, +# the user can specify all relevant conditions or predictors and their levels or value range. + +# The current implementation assumes a full factorial design (also called fully crossed design) +# in which each level of a factor occurs with each level of the other factors. Moreover, in the current implementation, there is exactly one instance of each of these factor combinations. + +# Example: +design_single = SingleSubjectDesign(; + conditions = Dict( + :stimulus_type => ["natural", "artificial"], + :contrast_level => range(0, 1, length = 3), + ), +); + +# In order to inspect the design, we can use the `generate_events` function to create an event table based on the design we specified. +generate_events(design_single) + +# To change the order of the trials e.g. to sort or shuffle them, one can use the `event_order_function` argument. +# Example: Randomize the order of trials +design_single_shuffled = SingleSubjectDesign(; + conditions = Dict( + :stimulus_type => ["natural", "artificial"], + :contrast_level => range(0, 1, length = 3), + ), + event_order_function = x -> shuffle(MersenneTwister(42), x), +); +# ```@raw html +#
+# Click to expand event table +# ``` +generate_events(design_single_shuffled) +# ```@raw html +#
+# ``` + +# ## Multi-subject designs +# The `MultiSubjectDesign` type can be used to simulate data for an experiment with multiple subjects. +# One needs to specify the number of subjects `n_subjects` and the number of items `n_items` i.e. stimuli. +# In addition, one needs to decide for every experimental factor whether it should be between- or within-subject (and item). + +# !!! note +# For factors that are not listed in `items_between` it is assumed that they vary within-item () + +design_multi = MultiSubjectDesign( + n_subjects = 6, + n_items = 4, + items_between = Dict(:colour => ["red", "blue"]), + subjects_between = Dict(:age_group => ["young", "old"]), + both_within = Dict(:luminance => range(0, 1, length = 3)), +); +# Caution, designs have to be balanced + +# ```@raw html +#
+# Click to expand event table +# ``` +generate_events(design_multi) +# ```@raw html +#
+# ``` + +# ## Repeat designs +# The `RepeatDesign` type is a functionality to encapsulate single- or multi-subject designs. It allows to repeat a generated event table multiple times. +# In other words, the `RepeatDesign` type allows to have multiple instances of the same item/subject/factor level combination. + +# Example: +# Assume, we have the following single-subject design from above: +# ```@raw html +#
+# Click to expand event table +# ``` +generate_events(design_single) +# ```@raw html +#
+# ``` + + +# But instead of having only one instance of the factor combinations e.g. `stimulus_type`: `natural` and `contrast_level`: `0`, we will repeat the design three times such that there are three occurrences of each combination. +design_repeated = RepeatDesign(design_single, 3); +generate_events(design_repeated) + +# [Here](@ref howto_repeat_design) one can find another example of how to repeat design entries for multi-subject designs. \ No newline at end of file diff --git a/docs/make.jl b/docs/make.jl index efe43840..72dc6219 100644 --- a/docs/make.jl +++ b/docs/make.jl @@ -39,6 +39,7 @@ makedocs(; ], "Reference" => [ "Overview of functionality" => "./generated/reference/overview.md", + "Overview: Experimental design types" => "./generated/reference/designtypes.md", "Overview: Basis function (component) types" => "./generated/reference/basistypes.md", "Overview: Onset types" => "./generated/reference/onsettypes.md", "Overview: Noise types" => "./generated/reference/noisetypes.md", From 89c09a1e71f3c9fc22e608a25933e52e3ad43b45 Mon Sep 17 00:00:00 2001 From: Benedikt Ehinger Date: Wed, 24 Jul 2024 13:05:46 +0200 Subject: [PATCH 16/37] Update quickstart.jl --- docs/literate/tutorials/quickstart.jl | 6 +++++- 1 file changed, 5 insertions(+), 1 deletion(-) diff --git a/docs/literate/tutorials/quickstart.jl b/docs/literate/tutorials/quickstart.jl index 4d93bb68..66d735a9 100644 --- a/docs/literate/tutorials/quickstart.jl +++ b/docs/literate/tutorials/quickstart.jl @@ -1,6 +1,10 @@ # # Quickstart -# To get started with data simulation, the user needs to provide four ingredients: an experimental design, an event basis function (component), an inter-onset distribution and a noise specification. +# To get started with data simulation, the user needs to provide four ingredients: +# 1) an experimental design, defining which conditions exist and how many events/"trials" +# 2) an event basis function, defining the simulated event-related response for every event (e.g. the ERP shape in EEG) +# 3) an inter-onset event distribution, defining the distances in time of the event sequence +# 4) a noise specification, defining, well, the noise :) # !!! tip # Use `subtypes(AbstractNoise)` (or `subtypes(AbstractComponent)` etc.) to find already implemented building blocks. From 1e2694f4b7b10bdf6b3700aadb322b9e471427c0 Mon Sep 17 00:00:00 2001 From: Judith Schepers Date: Wed, 24 Jul 2024 18:37:56 +0200 Subject: [PATCH 17/37] Add UnfoldSim logo to the documentation --- docs/assets/logo.svg | 28 ++++++++++++++++++++++++++++ 1 file changed, 28 insertions(+) create mode 100644 docs/assets/logo.svg diff --git a/docs/assets/logo.svg b/docs/assets/logo.svg new file mode 100644 index 00000000..340ebbd2 --- /dev/null +++ b/docs/assets/logo.svg @@ -0,0 +1,28 @@ + + + + + + + From 9ab04aada02234dc611823279c0cea04dce33e22 Mon Sep 17 00:00:00 2001 From: jschepers Date: Wed, 24 Jul 2024 17:12:26 +0000 Subject: [PATCH 18/37] Finished experimental design reference page --- docs/literate/reference/designtypes.jl | 13 ++++++++++--- docs/{ => src}/assets/logo.svg | 0 2 files changed, 10 insertions(+), 3 deletions(-) rename docs/{ => src}/assets/logo.svg (100%) diff --git a/docs/literate/reference/designtypes.jl b/docs/literate/reference/designtypes.jl index d21693bf..3d4cb98c 100644 --- a/docs/literate/reference/designtypes.jl +++ b/docs/literate/reference/designtypes.jl @@ -52,12 +52,12 @@ generate_events(design_single_shuffled) # ``` # ## Multi-subject designs -# The `MultiSubjectDesign` type can be used to simulate data for an experiment with multiple subjects. +# The `MultiSubjectDesign` type can be used to simulate data for an experiment with multiple subjects. Internally, it uses the [MixedModelsSim.jl package](https://github.com/RePsychLing/MixedModelsSim.jl). # One needs to specify the number of subjects `n_subjects` and the number of items `n_items` i.e. stimuli. # In addition, one needs to decide for every experimental factor whether it should be between- or within-subject (and item). # !!! note -# For factors that are not listed in `items_between` it is assumed that they vary within-item () +# For factors that are not listed in `items_between` it is assumed that they vary within-item (accordingly for `subjects_between`). design_multi = MultiSubjectDesign( n_subjects = 6, @@ -66,7 +66,6 @@ design_multi = MultiSubjectDesign( subjects_between = Dict(:age_group => ["young", "old"]), both_within = Dict(:luminance => range(0, 1, length = 3)), ); -# Caution, designs have to be balanced # ```@raw html #
@@ -75,8 +74,15 @@ design_multi = MultiSubjectDesign( generate_events(design_multi) # ```@raw html #
+#
# ``` +# As with the `SingleSubjectDesign` one can use the `event_order_function` argument to determine the order of events/trials. + +# !!! important +# The number of subjects/items has to be a divisor of the number of factor level combinations, i.e. it is assumed that the design is balanced +# which means that there is an equal number of observations for all possible factor level combinations. + # ## Repeat designs # The `RepeatDesign` type is a functionality to encapsulate single- or multi-subject designs. It allows to repeat a generated event table multiple times. # In other words, the `RepeatDesign` type allows to have multiple instances of the same item/subject/factor level combination. @@ -90,6 +96,7 @@ generate_events(design_multi) generate_events(design_single) # ```@raw html #
+#
# ``` diff --git a/docs/assets/logo.svg b/docs/src/assets/logo.svg similarity index 100% rename from docs/assets/logo.svg rename to docs/src/assets/logo.svg From 5671e2e71d602f016a90d89d6650c6eca3a414b3 Mon Sep 17 00:00:00 2001 From: Judith Schepers Date: Thu, 25 Jul 2024 15:09:27 +0200 Subject: [PATCH 19/37] Replace logo file --- docs/src/assets/logo.svg | 29 +++++++++++++++-------------- 1 file changed, 15 insertions(+), 14 deletions(-) diff --git a/docs/src/assets/logo.svg b/docs/src/assets/logo.svg index 340ebbd2..676efd96 100644 --- a/docs/src/assets/logo.svg +++ b/docs/src/assets/logo.svg @@ -1,26 +1,27 @@ - - - >1pe z!mT2nXn60-LTv1Zt}1rsQ*lS-*PR#O+U}`uQGQ<3adGr2u4ZVo)pf*EH+g(96DkT4 z2Hy)9G^FMTGrpp-iGq``N+v2?vsAdn*Z4%C;9>Dp(D!ZhnOEzKhZOqS9ZT~41W|9f z)^+aIZe|UO3q(JVuc(}3N`*4R54EbN4*g~fH2`rYc1kNRzPJ5V+OXU3pSS|?#TLLeP`!eglJn`3M)8)$PZkBXzc)DWnUDr@F zhjcMRKzD~gk`2=cr&R0u0o!42ntz=;EGkXjUhX%c6S4}@XqndavhA6fU;32X9){n& z!lw>Xt%)LGa@<=uy&}K9pFi>D51Jg!DR7m7n{+?dx+r0^FUmCqwU~U@Y8xqzrYlP4d&`jZ*)TP0x=;RdHxIpRz0^| z(Th2iEdQ}gN6G%Sg`&ryY>%mL7rdzjM)iP1PVM4LX?_PcFqRFC3oQUcOa?<)X0n+j zO8I>2DR=<>v0`eP--`Z628z{AGwumK5j+ul^O(uKKQY9`42BJhN5}|zSEX|z1Rc@$ z_FEn3n{mpdz}V{Gv%bG~NHOr!4|4x8?u*2>Or6AiPf^9NB4XF^Jv6ij1r>#loc!74 zHn6EHU{kovHnXI*C~G<%vXK)x>Dvayr@8O}1Wg^27EC{f!pv6e$}tRLt5{Z^O~&@+ zHHO?6h(`D+1Kl|U;<$l}A{dO#!C~Z_qRc^6p(!hFodKs&nQYQ4Z3)k(FRC_dxNL4% zgr9qDg7eIBCgI?8Zmx@RrIM(71v{Je`SGSUEbr5e17@WAAwnj8N|=pz0adU;j;o-C z5aP8EtfixaDKXR;Q5c@eX&Wju9SCiP`r1ND>vC(vu?{Mcw;6rDi-Wb^rLL`%1{KniG24D!2obV zIB=;5sj-Y~O=4B>%nXL#h(j^a5~~8Q+g`g4-cwE`5R&(PX01Li+$_3Fz9YG@G4PfA zYEzIJshmZ|+C0{fAaH?G^`2k&-4TV;01+;eqBAPH9B?>=@SVWC-O$WYwkyI&qPQA@ zPlBpyGF8f2lT~DauO3^b*zP$r(LP`H7jlss2t1Dl%FxK5lHu_@&%>AYfWa~~i&ud2 z$nxu*95`FYZBbTweiwL7e`WMTLy_3uXmbxpY4ts?NQ)V)TY%4d2gL{H@3C$_RRx(9`XOENy7(1Ln7+o|foL<| zbyPYJeuNf7i48<^v$mDmPJnI_DJH-8|GbTMzAk#)aq` z0Qu)2<}kMYYQQ$aQRJ8~MKhWK)`c#G=T2+pC}eg^=bMCACXyPjURaMP)c3^fH}VC- zqS_nNb~E(B!21fO>Pj!I^MPNY`+h5pQMa)!8Vx=%7r{6JQSrvUeZZ2^a=n4JSSx4o zY)3))s_&fJj+rxq?#LE`Zbp={^pwtSmvdY&ZXVt+kBL4giQ@V{xukU2_UX*`$f?tM z2n8}4h4Xrml;9YypWIj=*_NYKd@7FRFb*_`6ZTxHs!($QfgI(wSYz9{Pos>x2zYb2 zsqNj&zf5KvX6^0&HyMSFUE-u-fghns>z62fR|gwz*3Odnsl* zsLVJ1X~{FX<0%CTor~Pt)vxW6AY)^y&11jB74}LES!J8Y#y#oICkO6dHIFH#Y#NG& z$okZ->`++k6qkrw$Y=Ru3P3ZnM{lMvh@7`SzcnP>ROsa6HMc3^t_(jrQ%=X!i^J`uRIu^=rgyY~XE`UJToDaBISnz*;3tcRrm+5A==w$hwAm zYzdj2-cr}(r2wazE?*b0S245&9M)vd+{1$!hCXtCY0;6`-IU;k-qF^*UJCDNPn>%q zYBRTPuo@2fT|@v^K^UTqW~#u4N@O#xhEu6^9X#`;V_wB{If*xPF+|yBi~!I<@2iim zRqok;0u4|%!!_pKYsg{ZaA6ejHm1q%~uZ&`tUOE)9#!r9@_yk(T~C ziTYBgV4YRuQBY0_VA2IWoq%U4JMUc^wy&-4u9cY@>t_n+Bwddd(-gvkQcO|K387JBoTv zN&bwn{q{p$+-a7VFD#MPX$bLTV`d~6t?w7=L~hr^!Qp!J5G&BN4a5eN4w-%(`+~lZ zsEPTQ#q+@D*sC1!XrBX9j_B?ZB`BGW--S7xuT0+|c0B+eYdV^DR}l}fVt>`&Z$NZ8 z6Orp)=VD{!gWSa5zKFUIhLB7{Lo@NJZFUMuLg#fwUZuuTb0rkRpP5cl`4sEHKZm^4 z_rpacesmk;$dYA>djBKTN(o44ppB}Vyn}6GLNdnrvgh~^U;Gf{!>6Bh-fV5hLPHM1UTb&qgx{r>F@9acvkT5KRhrW`4mi4DxISgSYZ;c-$`t^aJrM2QU66}PhD4*GW)v(riGl7Bf?W(xMT_h zMMLXbYk`)mHt43mx_|;KKzky3SO#qxtO-khi|mm>OQ*PNMe*0rHcv+w*Eb-E z+c45UQL#m=CUB4vC&M71M9j)xpay&?!XRUG`7}n zfdYxqBk01J_P31Lm6m;Wg{C|e(pMqf zk74xT$jX(3G?(!T%|+L*Rrv18rxKaD;6w>z94~%!=MRL*6LpD;8`YO@!EJi{kcQc6 zHx`yQf~cTm&kY|RKY@8^PgZI@5%5#9ggkV1#=VGR@(7JAK%hP$KPI=FVlX= zdFq4f!?1Sl<>!y%V?x|1BTYID93JFHi=tICaCN5T;8G+)2IUMTG_>L6)x0c}TN_{$ zvqW*Pq}?9FBOEi=AFR`Yatj;JlS+foeO}%P-X5l&8B4YplftqRNjeR!@3Lazm2cP& zgDWX_&*_ieB7Ay$0x`!Bs^>~o?99b18$*%q!A#ZS+t0`AbNsZ0@~_^X($mm5eg;Ax zhGNip`zsFI$)kg~>D!})`olgEVXyF+U4Zp7)9`m^#FH|O?cjH6K^oGsucto-2)0K? zL`1yy*N=63mjc(AvXJW$uVN4-B`xjuivg|nVCXB!J2bQaQy@b z`|&82o1I?rSXA!)s!{!T5j{L{^2~B3^c~3a=pIuJe~@w*!ackw6V?qQ3pvEW!CKH( zyLptvET z5hp=#=O(ITyv?3VC_G6aT_$VfHGdkEBGz^k~cmU_5DX;g>Fu@P-Sk#<; z62KZhJzfKSocC$jV@o z(*%t0`+r^Mw&Xx>cJ`@NJ6*3h>VIFic(?{NqPr-D7W|e8ebk2k90>=y=hrEp!Fl=r z&v_1nxX0NF{O?({r157;T!jvS@Z!Vc7_pRb79wWoig=StN>HV?hH=)YrAhPI@AdIW!o{~^qlA@_w=0n$Yph5sGQW=q-{1l zlr9uG^nD0&NMri3cOZ}J!TS_!4e6FxenYrR#zQK9A|wa`C6(*0YLv-`d#_PsM|CZP+=m0QQ&WHjK-Mo=baZW+O&q zj>Aro?v34EKDA>AebQbOSpwebONhr1Q#hR1y=&Lig-PCqV@$JlBUGK%+!v$9lOJQ- zXR)EIB8J# zqutHl+xk)(Po-u`!rahJQN;ASx;}&iBLu^<)7@4ylPk5)YQYo}jP2=0+l?*A!Vz-G z*4Pi>aye>x+xTz9R(O-FZn~Rq>};$}e%P1HAxFB0^9gl;O>hz~9jc-TH*w0UWAZ~M zOZYmo1sZi~*H|=2>&C9JZY5V2mNRT_=%w%~l2HgYOL9|5d}=jARka+0*j=68wV@_b z`s=(P_n$)5d2{=T?VW_Ui}JZ09qx3(bbgL$7>e<}=x_vn1Qzj9&o(Y>umH=MtiI_F*I}b}D4Gavq93IaSH3{t4uiUzNg#T@8$9uUF_?8@eq}OLZ>`jkQN_2$|)h+6vvfB@68KsN| zzbibUjNms!6?k%~m5rN?b_(y0=F=9V6m8kylPwVeCoVDH;ze zq6Ij@iR_K2z^M$qB}gyy;*C_cKb7DE<< z2xdYS>|dm*G;jHDu#WJ|H%7kwrmBn94~=K?t_Pb6KXE2;ZV2HQj>3y5=C>g^v4zPx zd$sTzGKQZI9-Se~n?zQ>OPPR7n4REQELagbH@N7pZR;#K_&R&hxzL@VvXv@7M|#7K zN_hk+liIpWWK$d+PwDTiPpC#%i|R z#>$V1sC>$^XAHxK`N14u(@^{{-t~KQ5pETJDtPuW_o+b14f>?05$%lj8sDXTg_V+8 zHjCYxAb`rYEZFzSGeq1z5lB3)ojrtCK?=gzlWzJMgtwx1irdRBD8hXdytKDfcj~Pk zdl)}JB07>;uHZ8g;E>fSrs@+&p`_R6ZfGVBdjccv@2vu~d|Lw^8^8(A_qid$&JDKtj$OfH`9Zkd2xX}BqsG6K z{s~Vf3@4Zc{p&vO+PfP~$oi?_6F1<%#X7@WTL~dYI`arNme?7>WIjp^a}jtpOcOqm z#Wjlxj8xin5JgT$nWoP=Ur|YYs=k_XS<=|+c4p72-*?%sdSdQ8z z4kJXD>}MUrn`&({YvY}H^yg+9Io%v~i3c}f{PO~?2ldx-q7_$YbQb)MCu@c|>~IH; zFjss08st(61f%ILH^S`yZ~-LzwbLC9T=_2Qqe0tdIFoO?oj$BnRJF8nlWaz z;OGk@a&Azkv={ZS9m9|1*0+*EMeqBN7;4zdoSq*05EXW8tvsb4%XuF|xLM)s1KZy< z?lm2Q(Eq^Ndyi({!Y;D9OhtE}?P@M{!c2jv_babyYYk%TXi`py%FrvPlCZ{F_L_xS;3&OhU-*{< zjLY?4P!?rl*{*Q~z5$LF9L zVHoXeEYIXipoDkA81d0S^^S#J^=D|KIITRHdNVs{gt?0jf$-$r{b@{^O$;zF& zBPC44xxr}FXgqOrkgq*`I+Ne$iV^Fe7#*7E*_Fc)L_ujb58StcynW_3fHA$uJddLaMXll^r?xViK?D-$CC zIRY6$d}L-;PG&aol2o4ouCs{OMo@zky6_>sowU}p-H~X^L`rA-b#qFT@{ndwyPbP0 z%Heuobv?_2!mZTlwQp+Jo{nw&CO=6l6{?P6mT(<-W6*1m5z&oiqX(qyv9@Km%@XfA zw&L^^TRQm5PqT@o`f&OStzhnct4Cogr5jz8)XP7Mj&5a@BDB9P$fe%IGhKGUo57rE zn)+;2&m&P!nRgmEEXAY4`wjeB6vIJT!1UuR!w9=`_qgTQy%@e;e$GF)RE3)-WDe~P zfV&LL-wk9K6-9MOt-d)4o9N+Hj4Ry4dkDH8InC01l+RZPZOK~llIQAz1>2TxO9aR* z98fPG_Sn!=g3Z`r^V>swuDsw+uBjvK`e#U4IjQ%H7jTrZdnP1$4NH9O zvniIrIRX5xPSX$;sTZTzbfXk10@IPoYKA%VDY>t+Zdhi&Ujq#AZ+%^21yvw@-mO#p zt6xEn@{!)DXDN?-xHZMRPGSeC`)g;x*m1gYIA*BIHgfyI*Y88MG_4!d9YLzp6VGi; zEvKz~lNJ@SGG)(Y7WCxXbp$zT9M@!Pv1v*aH#(3!6y(QD>H{t3P34eK-xNrA4? zP)K;m9uF)z-22U)LZS$`m@Uq2>J~hdgzn#=uP0fkXz?A}*hW5BaC=QZZBi0{g)=e} zD7CuadRE5|ik>Iv@VjYL{Nub7AC46#P~eHJdtT1OI|0yW z)^0^49E8lm^s;7!xaC%>y`|pF@&!xcI5<@6J>k*Jo|>?)M$4#>j}GF3t`$YI72?Qm z)>(+Z-m^A{lC-n)L2WH{KiW)bt=gxKqVAch-GD$Z*YNOitgO*4(mj;`VFX0~ z4}Rom3xu;`9Lyy18nz^t&N5!@lKm(w(yI&h&;ycIe$^q<*rCi3nSlM>#!xCgLrOly zYb(F$=*7%}L}5!={oUUIMN=pf7CmAlt)0sbzB8=byEeqD@V&fkvu~8^0&!>UZ4b6b z-pnTX)DDCZR9CNUdTk43{{u6Dz%rNC-_3no#a_L0!U|vs+we;c1jqGQreB-V_i3}m zR+pXOQVH0Pl030WEFDZ=ZU)8D4R48mP@A*Ytx?RYS8x~^5E<;zu|5$)G@7=|)F`8% znRimi@^7{?+~ZC3sco6&+4b5+zAAG|^-Xf17|ZFd&9cX#v66p>K9A2q()w~8GCmaT zg#d|j5EpelD4-Em$!d)tsjUkKK!>pf_cFp8!t_xkMQUPtHrt2I{)(%J)Pp@weJfR} z;`|uTkT@VwZdZ2u_9deiqoUGYk+6|WhkL!%L*x+6wpnso{;#>jg#@`q-X5z9J<73i zk7n+tW9k`F0nS&phvpP(TxUjyB_#Qx`P-nI}IXPpJ@2Zn6$Fkt$NA93&x_IIi@(> z^@%U3V)=Z%zd}#;;JrkdR(1`$oMDzv${Aaj#>jgG?gfhiJO)JjKR@sHEmNaKb?7V% z>Ca6*?C~<5oQFs6*H;8*@ZNL1y^9wOLSBX%A7PymX1~r`L^g|WdETo>SZuf{z-T9cBFgf z>#P%*WH_?I2YA?Ze??EovTtn}*4v829Em^1>|TLpqa}{ag!4-cDD>Kn{&D>$Anrt< zR&1GTu>EPU)ClQ4$5asU?PRp@H<3(33~i;I-$V{7{ZzfKu>6ii#%nQA;NtMp2Un8g z0m;nR)REi+Q{I~0W!tL}E+fgdfGzqZhpA7!KbhW43g>+5h$N|A`{S{iUfHa;O#JY1 z?>dS9#wO8G{8d{9Jd|?D!4uB8fyYJqk60zal+@HU>``|$_?pvUUO^8He>Zwi z17NPj$WQ&qKEc=-p<`jZ?so#oJUk$N>+Bmr7$UUKkDG87|lIe-n z>>hMC^JkKnX4laBILtY80qcijko8iAQ&Cp*&X#OzJBi)S*#4JE%Cih+&;6h_$`+Vj z{&D{=PkeWc2LHiM_Yq(^H#69BTh*7s8)T%!1Q~zqImC45uTqWhGvg^ob{OUv-=7gN zBnP$nEAr$IHukyHCZnfyFDcyx@mN^O-Kx1iS`lCqfI2wM$YmsQ;07(yr@vckI9gl{ z#+FuBDWbkFSPY%gT`8xx!OKXN90aA-Kkhs~A%x0%^3;=3hv=5i=Nz_2T6+d1sR;!n z_jeEi;37@);c@7Xg|+a|T^r*_y=ONWpI$GP40WHifOTsYVa~WNgippU`RtlyU&>GT zbr#tvBomGw5z}(PQQd6p^t^(rNiBXUtQ-X zhBp~5j(dkjvYT5rCuYRumFOf0X=v6vCxwbDQpos>^ssV`uGRDgD-v4RO6RW~3fV%= zQA6LO)15UkBt?|pFI|!zt@Erk=!hN%FpRcugTx)lvXuLaIFWCEHq^xTl|del(ciu`s<qwV~oU8>d|FaoKm^~OvJL*g_J ziyOj%SekyJj^MPcrCboD`ql2E!Po4Z3{=7yB3N&Xr$*b}NAFv89`yQ+#1#DMTSkYR zK{`@y;4uptkEiZ~kY~z5JzvkSUj-#sRoHq)-Y(=Eb~g*R&Bldx^)>JAhIz927h`;U zTuT2yGXqI4S!U5mtptT!56X$+BwL~pQnM@jxdL?&_n$QxVa>Ed&ks@mq8T}*z9Q2< zfP2#9AhEhIvgCj)h5l}G$7hQs3^v>tdI9h#z^kwddmdb9MYO2MXu;>U@!--tlTpLo zKHC;;@}Ow-W`2GeEuCjlh3qgi(%sv?c0N>9cYkgob3YqkWH@WPmiokhOxZo3NHR&L zo&LD}{Hool{Ea(vFALIG!JZk9%-)>aZ1>~cFj*@z9E~v|AE2c z=1-e2sXEs=k!-H66Ose+E%jm0dZ#^puZed-=f;HR*Ib4q5(?Y4X(NZ*0&Av15}yaK z^_|JVVLDDDrwh{&z#a(%*jy`e=vZoE6U+OR#)MSZu|z&ui0{5^d9h#h0N;U_m1Gh8 zIPZz~ZL%Qr>b{39$iFzc1zhzQlM=ZhD(Xyb29JGy_i9UQ_Nx=TiByjhB~fx{x~ zj#b0)iMWQ(d37P%A}RVVhQC#o&I0o#d>T_t`_oFn2w<%KBRxy}**y0q^{g%!RSEp_ zZLwz#emer&BPocn6q`-}lH4WSr0g#>eGCa~qxmy2^{AS#eVUMy9c1NGZLi_07&xr~ zDRX~1(%p(flSWi3bckdhq<{)bF@Cip5vi8LU8~FGbU<16M>)CX^?hZ?tbamhc;*%*0p4t{H?uE!o5njJ8Ln0Jz^7|9xZANEX$m zb5tNy(M%J<%Acd{r!xmCf{{M3zZ{Y7`uYs|Di*+J{m2+|!R_AFcbVK<>)jz+$Ce?L4~V@I0D4Zrl@!;ZsOrh$_7O*G5m9 zxKCA+lnAkQXV1U-Gr7!O#csuaO7issDM!^KD%8g?YGvHd4yxqXASjMkBRW9?Rq2&5 zb7NRkbZ(3wl}j?eA&7gArBLUxPC^h%8zwc_8(Nh1uz^%6>s2PocaO^z1KQ`oYZ-Pi ze*B5t`4pe&S9{2dhqh0oYUz;M`CA>S0F}LuiKLafiq}O}{qCZQ;Tc6!)^esYqoA4y z@$uTGPi6Bzr}{FLOKOp9Yaxh`Sck(8&~@N;Wn#nxuy-mD$Rz|GiW`VlR!32IdHAPX zYxX)l_XMsyVu`;5y(*Gbi+YNdD1JU#?EH@^4HwGY-n^3dwDbWRmD`yPi97eXQh-0D zPxF%dz^TF7LfDRTSLmR&ZLUU@#cVUG1%U?5V|-9Y6^V#s44exqsF&~j5AxJQs zqz5C_K1sG!9R%)2-nHc5(DQPq=1-j|Z9$WWn97jf9YX?ehV?pT19$le*gsG*n_R&E zp2km1sp8GleI~7BV|UM4S_K{{R7j)#xru76{v7}6{^}u_e>M6}H3!sphd{-nC>~tI znZl94Rh96b;5ZFPG@{}alEJ$bGIf7GlB5Z0isL1<1Xf-DclQ%X$p)-i{qZqhlyh5q z**1^Qo|FXrLpHJwWb3#3t zZ*;!4=RtHq9!!MO;Ft9|zT}7zUf{+%-0;dqaZgQwQ}>~?lCH!VilI$x}(8b3|+K<2UAlSe%-K@Hep`F_dX&nE zC`vMDHlW$>9`*z$wbX*%KaMOYS0}iN@kw3VlbZk4_;HqF@V;5INr|_G5ts`vW>I^Y z43j}Ng0;2xTwV{?Ey~|D8*s~zZ_pe)aXD7>B9{&-^F@840c|%8V>&`%Ra@!d(w1Wcs zZA&fd*>Ri6ku+pvZdw%L=0cpC6p=IA1EV_ar(~?<4#!Dn3fI>c?G01DCq?tC%u+AX~g9sbR)>hp)N~QC-esPg}k> zNBglA+53rXndF!lY|*Uq)kdUBs%NNC!G+Vt_DFbR16?`3!#U_mXnqcG$s(F0WAu*i z()Hom19(3IkJQ0lm|n%K2`SWc46*w!>nVP2xhr_#@{!F`E}y2nQ6hu=vx;NB0$9%TCXIFZPiHChF5mRsp zn(ZZ|bhoU`R4;qPNv~GfEi#urgF@(j74(mq99dUCoTmNYv`a&Ej$6I$;P$n-Y)v!g zcOaAs9^&VqrLNU93Ga_6*BM{BD_z>n^tuX)6E1B7x6Ls4%mlZsh3C`(}A{f?H-+#Kn(tem;1>n`J}cdEGyzT51t!1K2D-XfK0ngIn2f zcwH%i1sDOnQT*rvMQQ_64lRdri}qSijx@yQvIlU>n1NzT1|>Y*&neZHQ9rl81QDAy z%7eE6E;Lsve6~8N%L z$F#Q!Bw4Qe&lc<$_+weQhEQzvp_QAZGzDt1Z3JTUB(FeOIyFo+4~fp~osZ|-G1%Si zytj$*b;1z1L60G1WM!8#(-}C>jgN}i3gkWqvDWBaR3JZ)t9uOvE?I+J)h>((ExsB- z@(`CyJg9CpiZa)F1wxf7JWIJ{8ZP}O9E5y|Sd*#=tJj!_n-V5V%vKL^b=J3-uTl7v zZ#w6PKjl(IT6VokzAag^gPL-+A9GCdMW0o#)(=KA13{bb1i^Jqh@u(6vOgV8)T()< z$%XO_b>l4*Kkl!`?JTA{EL|VbeShPCJ7;h}(?`XH%v;S3T{+eSZ;oyIKbRa>^nRHj zF228zncdK!v@(s*$k1kpz!}+ZCn{}4^rOErZRV;JSLTkeire0O3@3ew9rsu;ln%Jx zKNT1JCe|<3S~fM8U%T$_HHm}Jo^Xq)&PXJ`{?Fc!?gPC^v8f=~Cf)emt{YkYL%mYSToPP4 z$54*X=oAlHR8k^KbV>O5jQMIc|Hch0n&#Chlo(nbY80gdoQK`jQoFhOqvUfQ7s0N| zPuK-BImWfQt2tRdf$CdLDF)TVH@qZFLWmxfgsl5GB@5I*S)s{V8z6TCnHD% zh?l{qRo25f>D*6)(kuqYmxa@u%)`FP|x-2@ZGTe`e#x z&Ry2EUD}4%iFcY?oleTLXjdo?h5DNgx_t{Of0>YvQn6~cA1_T0_>Yf)s*G;M$Ol*AKkACU=kj!W9XIHfp!aVnd zpVdT+^DSl3%})GCe8>9$w5vhk>^sjbo@Vz`NbWzXE-XxisIT07*ubJ%iA~5pxJLJIy{a33-n^A_^wowFnDc#zwMnQ22V1db8_JVm|JK3Yr$Rqq(PVP| zvntQd+z*mF_jDis`8{Qy*uUM;Ws8n&|d_3 zt#T&r_#v=37HjA*IWPSujWu*N*VV(7yW>1cE6QQtsddp|b*?0%=weZ!Y#jzlMIq_~ zvmAfys}!6Ht!gtPy+#*2d0xrL)ct#LBWRLi%_)MUnfT)?-Ve5rUb)?(i`R?U@|p4+ zTUZOJzLTLLoq24?px~)t-NlIdjtC0_M>j+pWIw@oJU*;;{&hI1xFg!p$`tx`qqijtH|N8I;N$sfD%xA46fom6S z)>d-=}J)L*^%u@>GOR%N4O00=pA+hePuimBS#S=65^{?T?Zg1a5#aU z7w`S;Cr7F_-e)&niH6fPIkXW=okC*^au){0mTw8O%GyayswAluP1$kTSHJW?bGZiY z4Eyba6)lbas)p4Pbm4UJ4eT3Fu+oB;Qv@>{C{vGs0chicp4It{On8J-<;xp}`t%*d z%lFSaTT-eYOpag=(uBi-3Yiw{Q1@9qbcu5$&9S@dPbB^aC6wzvH?FgEG@=Oc<55;+ zq!G3nE42%N4<3%ZIbPHya(ejsjpsJCL*?qa(%Z6{j(>M|WWgDavc7TsTsnGggMD^- z?^QTm)h%*S{?fM)0pj41Pg%rVf8WlTD9U!*KL%^TP8gh7JnC4hfv{5YpE{?d|yz48v>l5E0 z>4pnRey<}^P(Qx|{#lr+IQLxpfx$`#H$?oo^R2^}`$v#+8|ZBEE_o$hIN{y87_%`cs%aI&(#%Mnq=?EL`rv9? zEC$@dulE=#{v3smRh+oAN{88PAgFf2g&|j=l^d-y(EDERGo4Rt(FL{jzpEO}zYn4I zZ7#BbYzldF_Fu*=L^U?QOma8~s8oleUsguSG3`JlTKF@P^b3|DFVGG?gh}qd!stOq z)DNYE>0Ivs#-Uz(Q-CrIk*lWT)pZ5&p|8VK4{}BS-Y=<(LJbrMgg3$L^2OUO7EMTT zuSYFDr`U9&wEs&aqq=u3BY1}))S235jcnLG>Iej9&%9Pb!gO-~PY2wuv zryKJ6HzK!lk@(`+P!7``{`<1-byv3mdPsuRvx1bgpPvRJK|b|i$AZif^BzR{3`YsmaT-vg;aVT2E-_~`YlrL zX+Da1O}t8S>rsoXfna5Ncz0cM+r$7FXFarfuNGjm$*zZGO6n+ zV53ajZB|fyox%vb9n;$1t>RxzY-(n+BdlI{GZMiUr0 zQ>WJ19w=Su4(fLUL9HiFpVy+AwF}qv;@>&8d5#{!vUQq^N^W($N?Hd>KjCkbg)nPT zLGOD&sBy*K$uZxbx@}FOE|mFi9Rl^c>H%eP=5=CVkDE*eb z7wiTjl_puGB;mm#tah5x!4edPlLR>%K6X<~dpLbgQ?nu(tWI?mftvV*76t!h@KOgu z0PvStU#(OIS1{*r{WQ)ETu3eZDF#N;z088rMT z{|NtjY7OODJ2Qt*+0v#&8TUyjNUuPHw!4v32Tuq>Mv|~X>QvVZEsVVV!YFrGAKu}x zzWEk2`iWX%^xiuAhO;Eo%`I5xkZ+c@m{KHdkt*M4I9V}dm*cSoo#6EJ^+g;V{KcT) zM=!!m4YS4^6%MW9Bm%Z!%LV6uzYti!X=hu94lQs=eh{_FD}5DT;no@zbWTK18Q6_j zHn%bRT9v~+=A=K;bYIm8VtMKnm)_3D(J%O2{+gT$g+)ZXTMwilZ`wFMseqrSZMd?t z`)-Vd2oQFCs~OX1mSkL;XAuV4<;v7q^enu49 zXbLXRRlkz`E%sy|SvMB`-H3`eFtIUG6Lq&wNYnfkrn{#Udm6Tsn+?mbebUx+Uw5za z?O8rRlj*HB-k71RQq{)J1}v`XoS4c3W$5l-TYAwA8%>S!Xc1P-@7 z0*o!sG1OZXtDNK zzOL4MuB0Frs8T=<7=VdgAga@!M5#Au`shylvk^N>SiLDvby%Q;a#)^Vl_heHdlBmZ zdg1WZUHskZY2&};mFgx&6gYqLqw>VfQ9}EGwI`*OP}Va&ig*lOcRd#f@mU>;JCjP$fmk=iS~V`dZGx#FkONK!EOL{xdBhx^shj2$Gy+7R2YvLt1Wx;b`oyN6 zyG5a0JZO;*s3z=DfO|gyfCeCviJienH*GiH!cYM_0Tik=^fvsbWhRbbpL<=l3PywE zck12KQC1P(hCl0CZhEB_?R0&6bkW?I0wC1Yd&&X#PgH1qxx*a>vnfvh+rt~{+>jK- zDEZ!F@`C3hRHjg4W?f_ZKcM9Wtu(f4JLZlEi+$_s&Ea%Vje@6gf_cC;xmK0Oe@hie z(5{NGszKX;Y^xU6n%vx#PxBPo2xht^1eHt=VM_DxmAglqbv;#pVR?6QBqcrs{W84xP2-n10tlw=A!&<#Db z5B=y+9#IGB7fQHzCO_ADaATauqpoSl_w{I4A6vwjZR-ffAVti*9~t&I-%33$n>jmr z+~!<>D$UFpv|FdCQ5n~%TfWA%3m5}am~L67_5&+9{sX?98+2@7fAn49&X#uw(l1@# zZYD2zJFxB?NR+LaRq=N}oRGf#8GSyWteCcWx%D4m2U*UY%}J)q*n)ra?;8wG^vy@% zl{6lDgf#>dECw&+4TQF>=RW zZs>4V94pvyutFhZ1vfEQkH(or3K5g$0+9ad{+LO7$4}cg0T&hC=R1#zLDkdSc;sI; znLXU?Qa*PyZmJ~HpjM8*-RZFX3*+U>1ndVRm1}ahR+irAG7>;a{HsENl@It# zvRLdF7#;TcLkF2b&hTooD&CJyTaKXU3DwDMI9$iJ&dcrX9$#tN@bTcnf1^?QTA=Oa zSF4j*K}m+zimO{hc;Vx+KQo5_HV9EYq`HnlLAr#9=l{Ycig2P*p;_QkWB!p_%%G3p zJC@N!#+|KTi)8IBR<5@CIjvS7&i=OZ%c-odlN;dsQ9LZ0M8x;``!3C%Bj4-yKSoN_ z)Jk>&y~O~sFh^m_8SZt&ucSWgxeaUk{o%-o!H#%oTlg!u=l+L#i))^+Evrd;xwiR{ z;){Lx;iShgR%+i@0(hEPyQm~X%dzKf^aYOQT9b!QNemGDZhx?BN zKQif*G1km@1M)_7gqm#C3pI^?(FDXn+WEiZ&o#`r{3%3mp>-@G_PYpc(PzW?^4x$W zaJ&A$VS=t4fSw_5koP&oKwU|n`}9=9{C`db&B|K~m?j1&%~cRMCeo-}2Xl&8bw}US7l8 xzDP`VO!l|@a~?ZQ5sdfXWQO^#>C402`_Q=XPx%%u99Re#W+s-G%Z#1l{ts%4>tg@_ literal 0 HcmV?d00001 From 34ce887c74d2a7aae0604368a37a27badd75e8f3 Mon Sep 17 00:00:00 2001 From: jschepers Date: Tue, 30 Jul 2024 13:58:26 +0000 Subject: [PATCH 23/37] Added intro paragraph for Simulate ERP tutorial --- docs/literate/tutorials/simulateERP.jl | 27 +++++++++++++++++++++++--- 1 file changed, 24 insertions(+), 3 deletions(-) diff --git a/docs/literate/tutorials/simulateERP.jl b/docs/literate/tutorials/simulateERP.jl index b3798096..c71f613d 100644 --- a/docs/literate/tutorials/simulateERP.jl +++ b/docs/literate/tutorials/simulateERP.jl @@ -1,12 +1,33 @@ +# # Simulate event-related potentials (ERPs) + +# One subfield of EEG research focuses on so-called event-related potentials (ERPs) which are defined as brain responses time-locked to a certain event e.g. stimulus onset. +# The waveform of an ERP usually consists of multiple ERP components which denote the peaks and troughs of the waveform. + +# ERP components are characterized (and named) by their timing relative to the event, their polarity (positive or negative) and their scalp topography. +# For example, the N170 describes a negative deflection which occurrs roughly 170 ms after the onset of (certain) visual stimuli. +# Often, researchers are interested how a component (e.g. its amplitude or timing) changes depending on certain experimental factors. +# For example, N170 has been shown to be related to face processing and its amplitude is modulated by whether the stimulus is a face or an object e.g. a car. +# ([Source](https://neuraldatascience.io/7-eeg/components.html)) + +# Here we will learn how to simulate a typical ERP complex with P100, N170, P300. + +# ### Setup +# ```@raw html +#
+# Click to expand +# ``` +## Load required packages using UnfoldSim using CairoMakie using Random using Unfold using UnfoldMakie -# ## ERP Complex -# Here we will learn how to simulate a typical ERP complex with P100, N170, P300. +# ```@raw html +#
+# ``` -# Let's grab a SingleSubjectDesign and add a continuous predictor +# ## Simulation +# Let's grab a `SingleSubjectDesign` and add a continuous predictor design = SingleSubjectDesign(; conditions = Dict( From b5a60bfe9aa67f5a7d46e2865bc0b59f7506626a Mon Sep 17 00:00:00 2001 From: jschepers Date: Tue, 30 Jul 2024 15:27:36 +0000 Subject: [PATCH 24/37] Improved docstrings for single- and multi-subject design --- src/design.jl | 27 ++++++++++++++++++++++----- 1 file changed, 22 insertions(+), 5 deletions(-) diff --git a/src/design.jl b/src/design.jl index 42e02328..66955bd9 100644 --- a/src/design.jl +++ b/src/design.jl @@ -1,6 +1,9 @@ """ - MultiSubjectDesign + MultiSubjectDesign <: AbstractDesign +A type for specifying the experimental design for multiple subjects (based on the given random-effects structure). + +### Fields - `n_subjects`::Int -> number of subjects - `n_items`::Int -> number of items (sometimes ≈trials) - `subjects_between` = Dict{Symbol,Vector} -> effects between subjects, e.g. young vs old @@ -8,10 +11,10 @@ - `both_within` = Dict{Symbol,Vector} -> effects completly crossed - `event_order_function` = `x->x`; # can be used to sort, or e.g. `x->shuffle(MersenneTwister(42),x)` - be sure to fix/update the rng accordingly!! -tipp: check the resulting dataframe using `generate_events(design)` - +Tip: Check the resulting dataframe using `generate_events(design)` +### Example ```julia # declaring same condition both sub-between and item-between results in a full between subject/item design design = MultiSubjectDesign(; @@ -33,6 +36,11 @@ end """ + SingleSubjectDesign <: AbstractDesign + +A type for specifying the experimental for a single subject (based on the given conditions). + +### Fields - conditions = Dict{Symbol,Vector} of conditions, e.g. `Dict(:A=>["a_small","a_big"],:B=>["b_tiny","b_large"])` - `event_order_function` = x->x; # can be used to sort, or x->shuffle(MersenneTwister(42),x) - be sure to fix/update the rng accordingly!! @@ -42,7 +50,16 @@ To increase the number of repetitions simply use `RepeatDesign(SingleSubjectDesi If conditions are omitted (or set to `nothing`), a single trial is simulated with a column `:dummy` and content `:dummy` - this is for convenience. -tipp: check the resulting dataframe using `generate_events(design)` +Tip: Check the resulting dataframe using `generate_events(design)` + +### Example +```julia +design = SingleSubjectDesign(; + conditions = Dict( + :stimulus_type => ["natural", "artificial"], + :contrast_level => range(0, 1, length = 5), +); +``` """ @with_kw struct SingleSubjectDesign <: AbstractDesign conditions::Dict{Symbol,Vector} = Dict() @@ -133,7 +150,7 @@ length(design::AbstractDesign) = *(size(design)...) # ---- """ - RepeatDesign{T} + RepeatDesign{T} <: AbstractDesign Repeat a design DataFrame multiple times to mimick repeatedly recorded trials. ```julia From b81b60556d2f154ff833c3c53c7dcb4ba3bf7e6d Mon Sep 17 00:00:00 2001 From: jschepers Date: Tue, 30 Jul 2024 16:07:24 +0000 Subject: [PATCH 25/37] Fixed simulate docstring --- src/simulation.jl | 15 +++++++-------- 1 file changed, 7 insertions(+), 8 deletions(-) diff --git a/src/simulation.jl b/src/simulation.jl index 2376e8ff..945c7a67 100755 --- a/src/simulation.jl +++ b/src/simulation.jl @@ -6,9 +6,15 @@ Simulation( noisetype::AbstractNoise, ) = Simulation(design, [component], onset, noisetype) + +function simulate(design::AbstractDesign, signal, onset::AbstractOnset, args...; kwargs...) + @warn "No random generator defined, used the default (`Random.MersenneTwister(1)`) with a fixed seed. This will always return the same results and the user is strongly encouraged to provide their own random generator!" + simulate(MersenneTwister(1), design, signal, onset, args...; kwargs...) +end + """ simulate( - rng, + [rng::AbstractRNG,] design::AbstractDesign, signal, onset::AbstractOnset, @@ -31,13 +37,6 @@ Some remarks to how the noise is added: - The case `return_epoched = false` and `onset = NoOnset()` is not possible and therefore covered by an assert statement """ - - -function simulate(design::AbstractDesign, signal, onset::AbstractOnset, args...; kwargs...) - @warn "No random generator defined, used the default (`Random.MersenneTwister(1)`) with a fixed seed. This will always return the same results and the user is strongly encouraged to provide their own random generator!" - simulate(MersenneTwister(1), design, signal, onset, args...; kwargs...) -end - simulate( rng::AbstractRNG, design::AbstractDesign, From 369faff6b9cd8b0ca1d9e2baa73b7044a45d5ecc Mon Sep 17 00:00:00 2001 From: jschepers Date: Wed, 31 Jul 2024 11:37:32 +0000 Subject: [PATCH 26/37] Added cross references in docstrings --- src/design.jl | 3 +++ 1 file changed, 3 insertions(+) diff --git a/src/design.jl b/src/design.jl index 66955bd9..e7bcdbcc 100644 --- a/src/design.jl +++ b/src/design.jl @@ -24,6 +24,7 @@ design = MultiSubjectDesign(; items_between = Dict(:cond => ["levelA", "levelB"]), ); ``` +See also [`SingleSubjectDesign`](@ref), [`RepeatDesign`](@ref) """ @with_kw struct MultiSubjectDesign <: AbstractDesign n_subjects::Int @@ -60,6 +61,7 @@ design = SingleSubjectDesign(; :contrast_level => range(0, 1, length = 5), ); ``` +See also [`MultiSubjectDesign`](@ref), [`RepeatDesign`](@ref) """ @with_kw struct SingleSubjectDesign <: AbstractDesign conditions::Dict{Symbol,Vector} = Dict() @@ -163,6 +165,7 @@ designOnce = MultiSubjectDesign(; design = RepeatDesign(designOnce,4); ``` +See also [`SingleSubjectDesign`](@ref), [`MultiSubjectDesign`](@ref) """ @with_kw struct RepeatDesign{T} <: AbstractDesign design::T From 541cf6712e01278bf6c4e4767d017f942732b0a8 Mon Sep 17 00:00:00 2001 From: jschepers Date: Wed, 31 Jul 2024 12:54:54 +0000 Subject: [PATCH 27/37] Added intro sentences, matched titles and sidebar, reordered pages and added collapsible setup blocks --- docs/literate/HowTo/multichannel.jl | 18 +++++++++++++++--- docs/literate/HowTo/newComponent.jl | 18 ++++++++++++++---- docs/literate/HowTo/newDesign.jl | 23 +++++++++++++++++------ docs/literate/HowTo/predefinedData.jl | 10 ++++++++++ docs/literate/HowTo/repeatTrials.jl | 6 ++---- docs/literate/reference/basistypes.jl | 19 ++++++++++++++----- docs/literate/reference/designtypes.jl | 2 +- docs/literate/reference/noisetypes.jl | 7 ++++++- docs/literate/reference/onsettypes.jl | 2 +- docs/literate/reference/overview.jl | 10 ++++++++++ docs/literate/tutorials/multisubject.jl | 15 ++++++++++++++- docs/literate/tutorials/poweranalysis.jl | 22 ++++++++++++++++------ docs/literate/tutorials/quickstart.jl | 4 ++-- docs/make.jl | 15 ++++++++------- docs/src/index.md | 4 ++-- 15 files changed, 132 insertions(+), 43 deletions(-) diff --git a/docs/literate/HowTo/multichannel.jl b/docs/literate/HowTo/multichannel.jl index 8aa5111f..50d4b5fa 100644 --- a/docs/literate/HowTo/multichannel.jl +++ b/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 +#
+# Click to expand +# ``` +## Load required packages using UnfoldSim using UnfoldMakie using CairoMakie using DataFrames using Random - +# ```@raw html +#
+# ``` # ## 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); diff --git a/docs/literate/HowTo/newComponent.jl b/docs/literate/HowTo/newComponent.jl index bb8bee72..91da1245 100644 --- a/docs/literate/HowTo/newComponent.jl +++ b/docs/literate/HowTo/newComponent.jl @@ -1,6 +1,12 @@ -# # 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. + +# ### Setup +# ```@raw html +#
+# Click to expand +# ``` using UnfoldSim using Unfold using Random @@ -8,6 +14,9 @@ using DSP using CairoMakie, UnfoldMakie sfreq = 100; +# ```@raw html +#
+# ``` # ## Design # Let's generate a design with two columns, shift + duration @@ -19,7 +28,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 @@ -51,7 +61,7 @@ function basis_shiftduration(evts, maxlength) end end - +# ## Simulate data with the new component type erp = UnfoldSim.simulate( MersenneTwister(1), TimeVaryingComponent(basis_shiftduration, 50), diff --git a/docs/literate/HowTo/newDesign.jl b/docs/literate/HowTo/newDesign.jl index 955044d2..ec0fc256 100644 --- a/docs/literate/HowTo/newDesign.jl +++ b/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 +#
+# Click to expand +# ``` 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 +#
+#
+# ``` + # #### 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 # diff --git a/docs/literate/HowTo/predefinedData.jl b/docs/literate/HowTo/predefinedData.jl index 9a394e9b..5ed27f74 100644 --- a/docs/literate/HowTo/predefinedData.jl +++ b/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 +#
+# Click to expand +# ``` +## Load required packages using UnfoldSim using DataFrames using Random using CairoMakie # for plotting +# ```@raw html +#
+#
+# ``` # 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]) diff --git a/docs/literate/HowTo/repeatTrials.jl b/docs/literate/HowTo/repeatTrials.jl index 1add60e6..b85db6b0 100644 --- a/docs/literate/HowTo/repeatTrials.jl +++ b/docs/literate/HowTo/repeatTrials.jl @@ -1,12 +1,10 @@ -using UnfoldSim - - -# # [Repeating Design entries](@id howto_repeat_design) +# # [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, diff --git a/docs/literate/reference/basistypes.jl b/docs/literate/reference/basistypes.jl index de4db4a6..188951bd 100644 --- a/docs/literate/reference/basistypes.jl +++ b/docs/literate/reference/basistypes.jl @@ -1,14 +1,23 @@ -using UnfoldSim -using CairoMakie -using DSP -using StableRNGs - # # 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 +#
+# Click to expand +# ``` +## Load required packages +using UnfoldSim +using CairoMakie +using DSP +using StableRNGs +# ```@raw html +#
+# ``` + # ## EEG # By default, the EEG bases assume a sampling rate of 100, which can easily be changed by e.g. p100(; sfreq=300) f = Figure() diff --git a/docs/literate/reference/designtypes.jl b/docs/literate/reference/designtypes.jl index 3d4cb98c..7d571f47 100644 --- a/docs/literate/reference/designtypes.jl +++ b/docs/literate/reference/designtypes.jl @@ -3,7 +3,7 @@ # The experimental design specifies the experimental conditions and other variables that are supposed to have an influence on the simulated data. # Currently, there are three types of designs implemented: `SingleSubjectDesign`, `MultiSubjectDesign` and `RepeatDesign`. -# ## Setup +# ### Setup # ```@raw html #
# Click to expand diff --git a/docs/literate/reference/noisetypes.jl b/docs/literate/reference/noisetypes.jl index 26af0e1e..264e3439 100644 --- a/docs/literate/reference/noisetypes.jl +++ b/docs/literate/reference/noisetypes.jl @@ -1,5 +1,10 @@ # # Overview: Noise types -# There are several noise types directly implemented. Here is a comparison: + +# There are different types of noise signals which differ in their power spectra. +# If you are not familiar with different types/colors of noise yet, have a look at this [source](https://en.wikipedia.org/wiki/Colors_of_noise). + +# There are several noise types directly implemented in UnfoldSim.jl. Here is a comparison: + using UnfoldSim using CairoMakie diff --git a/docs/literate/reference/onsettypes.jl b/docs/literate/reference/onsettypes.jl index bbadef15..6618706b 100644 --- a/docs/literate/reference/onsettypes.jl +++ b/docs/literate/reference/onsettypes.jl @@ -2,7 +2,7 @@ # 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`. -# ## Setup +# ### Setup # ```@raw html #
# Click to expand diff --git a/docs/literate/reference/overview.jl b/docs/literate/reference/overview.jl index 80145bcf..a4f4b611 100644 --- a/docs/literate/reference/overview.jl +++ b/docs/literate/reference/overview.jl @@ -1,7 +1,17 @@ # # Overview of functionality # UnfoldSim has many modules, here we try to collect them to provide you with an overview. + +# ### Setup +# ```@raw html +#
+# Click to expand +# ``` +## Load required packages using UnfoldSim using InteractiveUtils +# ```@raw html +#
+# ``` # ## Design # Designs define the experimental design. They can be nested, e.g. `RepeatDesign(SingleSubjectDesign,10)` would repeat the generated design-dataframe 10x. diff --git a/docs/literate/tutorials/multisubject.jl b/docs/literate/tutorials/multisubject.jl index b63e58d4..61570240 100644 --- a/docs/literate/tutorials/multisubject.jl +++ b/docs/literate/tutorials/multisubject.jl @@ -1,10 +1,23 @@ +# # Multi-subject simulation + +# In this tutorial, you will learn how to simulate data for multiple subjects. In particular, you will learn how to specify fixed and random effects and what their influence on the simulated data looks like. + +# ### Setup +# ```@raw html +#
+# Click to expand +# ``` +## Load required packages using UnfoldSim using Unfold using CairoMakie using UnfoldMakie using DataFrames +# ```@raw html +#
+#
+# ``` -# # Multi-subject simulation # Similar to the single subject case, multi-subject simulation depends on: # - `Design` (typically a `MultiSubjectDesign`) # - `Components` (typically a `MixedModelComponent`) diff --git a/docs/literate/tutorials/poweranalysis.jl b/docs/literate/tutorials/poweranalysis.jl index 5202ad66..fbc79b85 100644 --- a/docs/literate/tutorials/poweranalysis.jl +++ b/docs/literate/tutorials/poweranalysis.jl @@ -1,14 +1,24 @@ +# # Power analysis +# For a power analysis, we will repeatedly simulate data, and check whether we can find a significant effect. +# We perform the power analysis on epoched data. + +# ### Setup +# ```@raw html +#
+# Click to expand +# ``` +## Load required packages using UnfoldSim using Statistics using HypothesisTests using DataFrames using Random +# ```@raw html +#
+#
+# ``` - -# ## Simple Power analysis Script -# For a power analysis, we will repeatedly simulate data, and check whether we can find a significant effect. -# -# We perform the power analysis on epoched data. +# ## Simulation loop pvals = fill(NaN, 100) @time for seed in eachindex(pvals) ## Simulate data of 30 subjects @@ -33,5 +43,5 @@ pvals = fill(NaN, 100) pvals[seed] = pvalue(OneSampleTTest(y_big, y_small)) end -# let's calculate the power +# Let's calculate the power power = mean(pvals .< 0.05) * 100 diff --git a/docs/literate/tutorials/quickstart.jl b/docs/literate/tutorials/quickstart.jl index 66d735a9..0270a9a8 100644 --- a/docs/literate/tutorials/quickstart.jl +++ b/docs/literate/tutorials/quickstart.jl @@ -1,10 +1,10 @@ # # Quickstart # To get started with data simulation, the user needs to provide four ingredients: -# 1) an experimental design, defining which conditions exist and how many events/"trials" +# 1) an experimental design, defining which conditions and how many events/"trials" exist # 2) an event basis function, defining the simulated event-related response for every event (e.g. the ERP shape in EEG) # 3) an inter-onset event distribution, defining the distances in time of the event sequence -# 4) a noise specification, defining, well, the noise :) +# 4) a noise specification, defining the type of noise signal that is added to the simulated signal (e.g. pink noise) # !!! tip # Use `subtypes(AbstractNoise)` (or `subtypes(AbstractComponent)` etc.) to find already implemented building blocks. diff --git a/docs/make.jl b/docs/make.jl index 72dc6219..2b24ff62 100644 --- a/docs/make.jl +++ b/docs/make.jl @@ -19,7 +19,7 @@ DocMeta.setdocmeta!(UnfoldSim, :DocTestSetup, :(using UnfoldSim); recursive = tr makedocs(; modules = [UnfoldSim], - authors = "Luis Lips, Benedikt Ehinger, Judith Schepers", + authors = "Judith Schepers, Luis Lips, Maanik Marathe, Benedikt Ehinger", #repo="https://github.com/unfoldtoolbox/UnfoldSim.jl/blob/{commit}{path}#{line}", repo = Documenter.Remotes.GitHub("unfoldtoolbox", "UnfoldSim.jl"), sitename = "UnfoldSim.jl", @@ -27,13 +27,14 @@ makedocs(; prettyurls = get(ENV, "CI", "false") == "true", canonical = "https://unfoldtoolbox.github.io/UnfoldSim.jl", edit_link = "main", + sidebar_sitename = false, assets = String[], ), pages = [ "Home" => "index.md", "Tutorials" => [ "Quickstart" => "generated/tutorials/quickstart.md", - "Simulate ERPs" => "generated/tutorials/simulateERP.md", + "Simulate event-related potentials (ERPs)" => "generated/tutorials/simulateERP.md", "Power analysis" => "generated/tutorials/poweranalysis.md", "Multi-subject simulation" => "generated/tutorials/multisubject.md", ], @@ -45,13 +46,13 @@ makedocs(; "Overview: Noise types" => "./generated/reference/noisetypes.md", ], "HowTo" => [ - "Define a new, (imbalanced) design" => "./generated/HowTo/newDesign.md", - "Repeating a design" => "./generated/HowTo/repeatTrials.md", - "Define a new duration & jitter component" => "./generated/HowTo/newComponent.md", + "Define a new (imbalanced) design" => "./generated/HowTo/newDesign.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 predefined design / onsets data" => "./generated/HowTo/predefinedData.md", + "Use existing experimental designs & onsets in the simulation" => "./generated/HowTo/predefinedData.md", ], - "API / DocStrings" => "api.md", + "API / Docstrings" => "api.md", ], ) diff --git a/docs/src/index.md b/docs/src/index.md index 40ff2f14..e18f03e6 100644 --- a/docs/src/index.md +++ b/docs/src/index.md @@ -4,9 +4,9 @@ CurrentModule = UnfoldSim # UnfoldSim.jl -Documentation for [UnfoldSim.jl](https://github.com/unfoldtoolbox/UnfoldSim.jl). UnfoldSim.jl is a Julia package for simulating multivariate timeseries data with a special focus on EEG data. +Documentation for [UnfoldSim.jl](https://github.com/unfoldtoolbox/UnfoldSim.jl): a Julia package for simulating multivariate timeseries data with a special focus on EEG data. -## Start simulating timeseries +## Start simulating time series data We offer some predefined (EEG) signals, check them out! For instance a P1/N170/P300 complex (containing three typical ERP components). From f449083fdb72eec8b94e6eb83a5cad0c7ce2e1ce Mon Sep 17 00:00:00 2001 From: Benedikt Ehinger Date: Fri, 26 Jul 2024 13:26:21 +0200 Subject: [PATCH 28/37] Update noise.jl --- src/noise.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/noise.jl b/src/noise.jl index 28421ccd..ba709e17 100644 --- a/src/noise.jl +++ b/src/noise.jl @@ -72,7 +72,7 @@ Noise with exponential decay in AR spectrum. `noiselevel` is used to scale the noise !!! warning - Current implementation: Cholesky of NxN matrix needs to be calculated, which might need lots of RAM. + With the current implementation we try to get exponential decay over the whole AR spectrum, which is N-Samples long. This involves the inversion of a cholesky matrix of size NxN matrix, which will need lots of RAM for non-trivial problems. """ @with_kw struct ExponentialNoise <: AbstractNoise noiselevel = 1 From fe00ec2eb057f9b520cea8d056a3b7fb0df4c3b9 Mon Sep 17 00:00:00 2001 From: Benedikt Ehinger Date: Fri, 26 Jul 2024 14:59:57 +0200 Subject: [PATCH 29/37] Update src/noise.jl --- src/noise.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/noise.jl b/src/noise.jl index ba709e17..983f8270 100644 --- a/src/noise.jl +++ b/src/noise.jl @@ -72,7 +72,7 @@ Noise with exponential decay in AR spectrum. `noiselevel` is used to scale the noise !!! warning - With the current implementation we try to get exponential decay over the whole AR spectrum, which is N-Samples long. This involves the inversion of a cholesky matrix of size NxN matrix, which will need lots of RAM for non-trivial problems. + With the current implementation we try to get exponential decay over the whole autocorrelation (AR) spectrum, which is N-Samples (the total number of samples) long. This involves the inversion of a cholesky matrix of size NxN matrix, which will need lots of RAM for non-trivial problems. """ @with_kw struct ExponentialNoise <: AbstractNoise noiselevel = 1 From f93a3c19333a26bf7b95f77f9b39d22848c65cbc Mon Sep 17 00:00:00 2001 From: Benedikt Ehinger Date: Fri, 26 Jul 2024 15:00:22 +0200 Subject: [PATCH 30/37] Update src/noise.jl --- src/noise.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/noise.jl b/src/noise.jl index 983f8270..7ba1ddea 100644 --- a/src/noise.jl +++ b/src/noise.jl @@ -72,7 +72,7 @@ Noise with exponential decay in AR spectrum. `noiselevel` is used to scale the noise !!! warning - With the current implementation we try to get exponential decay over the whole autocorrelation (AR) spectrum, which is N-Samples (the total number of samples) long. This involves the inversion of a cholesky matrix of size NxN matrix, which will need lots of RAM for non-trivial problems. + With the current implementation we try to get exponential decay over the whole autoregressive (AR) spectrum, which is N-Samples (the total number of samples) long. This involves the inversion of a cholesky matrix of size NxN matrix, which will need lots of RAM for non-trivial problems. """ @with_kw struct ExponentialNoise <: AbstractNoise noiselevel = 1 From aca3fd2415221be8ca99183d2e4adf649f4de7b7 Mon Sep 17 00:00:00 2001 From: Judith Schepers Date: Wed, 31 Jul 2024 15:03:15 +0200 Subject: [PATCH 31/37] Update src/noise.jl --- src/noise.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/noise.jl b/src/noise.jl index 7ba1ddea..fcd376ec 100644 --- a/src/noise.jl +++ b/src/noise.jl @@ -72,7 +72,7 @@ Noise with exponential decay in AR spectrum. `noiselevel` is used to scale the noise !!! warning - With the current implementation we try to get exponential decay over the whole autoregressive (AR) spectrum, which is N-Samples (the total number of samples) long. This involves the inversion of a cholesky matrix of size NxN matrix, which will need lots of RAM for non-trivial problems. + With the current implementation we try to get exponential decay over the whole autoregressive (AR) spectrum, which is N samples (the total number of samples in the signal) long. This involves the inversion of a Cholesky matrix of size NxN matrix, which will need lots of RAM for non-trivial problems. """ @with_kw struct ExponentialNoise <: AbstractNoise noiselevel = 1 From e868a14f80a3ca37eb7303911d7fd8c9fa2cdb23 Mon Sep 17 00:00:00 2001 From: jschepers Date: Wed, 31 Jul 2024 13:36:38 +0000 Subject: [PATCH 32/37] add empty line for formatting reasons --- docs/literate/reference/designtypes.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/docs/literate/reference/designtypes.jl b/docs/literate/reference/designtypes.jl index 7d571f47..eca80e26 100644 --- a/docs/literate/reference/designtypes.jl +++ b/docs/literate/reference/designtypes.jl @@ -104,4 +104,4 @@ generate_events(design_single) design_repeated = RepeatDesign(design_single, 3); generate_events(design_repeated) -# [Here](@ref howto_repeat_design) one can find another example of how to repeat design entries for multi-subject designs. \ No newline at end of file +# [Here](@ref howto_repeat_design) one can find another example of how to repeat design entries for multi-subject designs. From 43dd57c9f63f59b570c1a5c3d7eb9b0054b8a911 Mon Sep 17 00:00:00 2001 From: Judith Schepers Date: Thu, 1 Aug 2024 10:22:17 +0200 Subject: [PATCH 33/37] Update docs/literate/reference/noisetypes.jl Co-authored-by: Benedikt Ehinger --- docs/literate/reference/noisetypes.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/docs/literate/reference/noisetypes.jl b/docs/literate/reference/noisetypes.jl index 264e3439..a14551eb 100644 --- a/docs/literate/reference/noisetypes.jl +++ b/docs/literate/reference/noisetypes.jl @@ -1,7 +1,7 @@ # # Overview: Noise types # There are different types of noise signals which differ in their power spectra. -# If you are not familiar with different types/colors of noise yet, have a look at this [source](https://en.wikipedia.org/wiki/Colors_of_noise). +# If you are not familiar with different types/colors of noise yet, have a look at the[colored noise wikipedia page](https://en.wikipedia.org/wiki/Colors_of_noise). # There are several noise types directly implemented in UnfoldSim.jl. Here is a comparison: From 68ce72283268e3f8c87dd616ffa4a03783e7ca4a Mon Sep 17 00:00:00 2001 From: Judith Schepers Date: Thu, 1 Aug 2024 10:22:27 +0200 Subject: [PATCH 34/37] Update docs/literate/reference/overview.jl Co-authored-by: Benedikt Ehinger --- docs/literate/reference/overview.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/docs/literate/reference/overview.jl b/docs/literate/reference/overview.jl index a4f4b611..f34c7f4e 100644 --- a/docs/literate/reference/overview.jl +++ b/docs/literate/reference/overview.jl @@ -1,5 +1,5 @@ # # Overview of functionality -# UnfoldSim has many modules, here we try to collect them to provide you with an overview. +# A UnfoldSim simulation has four ingredients: Design, Component, Onset and Noise. Here we provide a short overview of the implemented functions. # ### Setup # ```@raw html From 2dfd44becaa12a35d5ad854bbc83b37bed450b89 Mon Sep 17 00:00:00 2001 From: Judith Schepers Date: Thu, 1 Aug 2024 10:36:41 +0200 Subject: [PATCH 35/37] Update docs/literate/reference/overview.jl --- docs/literate/reference/overview.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/docs/literate/reference/overview.jl b/docs/literate/reference/overview.jl index f34c7f4e..74cf3bc1 100644 --- a/docs/literate/reference/overview.jl +++ b/docs/literate/reference/overview.jl @@ -1,5 +1,5 @@ # # Overview of functionality -# A UnfoldSim simulation has four ingredients: Design, Component, Onset and Noise. Here we provide a short overview of the implemented functions. +# A UnfoldSim simulation has four ingredients: Design, Component, Onset and Noise. Here we provide a short overview of the implemented types. # ### Setup # ```@raw html From 2c34e5e7bf9674073cf900cf0b4fff04e9289843 Mon Sep 17 00:00:00 2001 From: Judith Schepers Date: Thu, 1 Aug 2024 10:36:50 +0200 Subject: [PATCH 36/37] Update docs/literate/reference/noisetypes.jl --- docs/literate/reference/noisetypes.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/docs/literate/reference/noisetypes.jl b/docs/literate/reference/noisetypes.jl index a14551eb..03d35b4b 100644 --- a/docs/literate/reference/noisetypes.jl +++ b/docs/literate/reference/noisetypes.jl @@ -1,7 +1,7 @@ # # Overview: Noise types # There are different types of noise signals which differ in their power spectra. -# If you are not familiar with different types/colors of noise yet, have a look at the[colored noise wikipedia page](https://en.wikipedia.org/wiki/Colors_of_noise). +# If you are not familiar with different types/colors of noise yet, have a look at the[colors of noise Wikipedia page](https://en.wikipedia.org/wiki/Colors_of_noise). # There are several noise types directly implemented in UnfoldSim.jl. Here is a comparison: From 3a3bd3d74e7c3ef5f9b75c321a318ac8763f5137 Mon Sep 17 00:00:00 2001 From: "behinger (s-ccs 001)" Date: Fri, 9 Aug 2024 07:28:13 +0000 Subject: [PATCH 37/37] added docstring --- docs/literate/HowTo/sequence.jl | 39 ++++++++++++++++++++------------- src/component.jl | 8 +++++-- src/design.jl | 31 +++++++++++++++++++++++++- 3 files changed, 60 insertions(+), 18 deletions(-) diff --git a/docs/literate/HowTo/sequence.jl b/docs/literate/HowTo/sequence.jl index 6a07d15f..986446a6 100644 --- a/docs/literate/HowTo/sequence.jl +++ b/docs/literate/HowTo/sequence.jl @@ -5,18 +5,19 @@ using StableRNGs # ## Stimulus - Response design -# let's say we want to simulate a stimulus response, followed by a button press response. +# 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"]))#|>x->RepeatDesign(x,4) +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{1,2}_", 0, StableRNG(1)) +# 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*" +# 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 design 2 times - because we can +# Finally, let's repeat the current design 4 times design = RepeatDesign(design, 4) generate_events(design) @@ -24,36 +25,39 @@ generate_events(design) #generate_events(design) -# This results in 12 trials that nicely follow our sequence +# This results in 16 trials that nicely follow our sequence # Next we have to specify for both events `S` and `R` what the responses should look like. - p1 = LinearModelComponent(; basis = p100(), formula = @formula(0 ~ 1 + condition), β = [1, 0.5], -); +) n1 = LinearModelComponent(; basis = n170(), formula = @formula(0 ~ 1 + condition), β = [1, 0.5], -); +) + p3 = LinearModelComponent(; basis = UnfoldSim.hanning(Int(0.5 * 100)), # sfreq = 100 for the other bases formula = @formula(0 ~ 1 + condition), β = [1, 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]) -#components = [p1, n1, resp] + data, evts = simulate( StableRNG(1), design, @@ -61,8 +65,13 @@ data, evts = simulate( UniformOnset(offset = 40, width = 10), NoNoise(), ) +nothing ## hide +# Finally we can plot the results lines(data) -vlines!(evts.latency, color = (:gray, 0.5)) +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() \ No newline at end of file +current_figure() + +# As visible, the `R` response always follows the `S` response. Due to the "`_`" we have large breaks between the individual sequences. \ No newline at end of file diff --git a/src/component.jl b/src/component.jl index d3a0c5e7..65ef5e38 100644 --- a/src/component.jl +++ b/src/component.jl @@ -21,6 +21,7 @@ MixedModelComponent(; ``` """ +# backwards compatability after introducing the `offset` field` @with_kw struct MixedModelComponent <: AbstractComponent basis::Any formula::Any # e.g. 0~1+cond @@ -29,7 +30,8 @@ MixedModelComponent(; contrasts::Dict = Dict() offset::Int = 0 end - +MixedModelComponent(basis, formula, β, σs, contrasts) = + MixedModelComponent(basis, formula, β, σs, contrasts, 0) """ A multiple regression component for one subject @@ -51,6 +53,7 @@ LinearModelComponent(; ``` """ +# backwards compatability after introducing the `offset` field @with_kw struct LinearModelComponent <: AbstractComponent basis::Any formula::Any # e.g. 0~1+cond - left side must be "0" @@ -59,7 +62,8 @@ LinearModelComponent(; offset::Int = 0 end - +LinearModelComponent(basis, formula, β, contrasts) = + LinearModelComponent(basis, formula, β, contrasts, 0) """ offset(AbstractComponent) diff --git a/src/design.jl b/src/design.jl index d35ea440..740217bb 100644 --- a/src/design.jl +++ b/src/design.jl @@ -178,6 +178,34 @@ function check_sequence(s::String) @assert length(blankfind) <= 1 && (length(blankfind) == 0 || length(s) == blankfind[1]) "the blank-indicator '_' has to be the last sequence element" return s end + + +""" + SequenceDesign{T} <: AbstractDesign +Enforce 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. + + +```julia +design = SingleSubjectDesign(conditions = Dict(:condition => ["one", "two"])) +design = SequenceDesign(design, "SCR_", StableRNG(1)) +``` +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 +``` + +See also [`SingleSubjectDesign`](@ref), [`MultiSubjectDesign`](@ref), [`RepeatDesign`](@ref) +""" @with_kw struct SequenceDesign{T} <: AbstractDesign design::T sequence::String = "" @@ -187,7 +215,8 @@ end SequenceDesign{T}(d, s, sl, r) where {T<:AbstractDesign} = new(d, check_sequence(s), sl, r) end - +SequenceDesign(design, sequence, rng::AbstractRNG) = + SequenceDesign(design = design, sequence = sequence, rng = rng) SequenceDesign(design, sequence) = SequenceDesign(design = design, sequence = sequence) generate_events(design::SequenceDesign{MultiSubjectDesign}) = error("not yet implemented")