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 1/5] 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 fc921d7..a6989c2 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 b4a8ac8e8415aa433e34b44eb2c97dcffb2acada Mon Sep 17 00:00:00 2001 From: "behinger (s-ccs 001)" Date: Sun, 10 Mar 2024 19:13:52 +0000 Subject: [PATCH 2/5] generalized LinearModelComponent to arbitrary functions instead of vectors --- docs/literate/HowTo/componentfunction.jl | 39 ++++++++++++++++ docs/literate/HowTo/newComponent.jl | 3 ++ docs/make.jl | 3 +- src/component.jl | 57 ++++++++++++++++++++++-- 4 files changed, 97 insertions(+), 5 deletions(-) create mode 100644 docs/literate/HowTo/componentfunction.jl diff --git a/docs/literate/HowTo/componentfunction.jl b/docs/literate/HowTo/componentfunction.jl new file mode 100644 index 0000000..7d36c62 --- /dev/null +++ b/docs/literate/HowTo/componentfunction.jl @@ -0,0 +1,39 @@ +# # Component Functions +# HowTo put arbitrary functions into components + +using UnfoldSim +using Unfold +using Random +using DSP +using CairoMakie, UnfoldMakie + +sfreq = 100; + +# ## Design +# Let's generate a design with a categorical effect and a continuous duration effect +design = UnfoldSim.SingleSubjectDesign(; + conditions = Dict( + :category => ["dog", "cat"], + :duration => Int.(round.(20 .+ rand(100) .* sfreq)), + ), +); + + +# Instead of defining a boring vector basis function e.g. `[0,0,1,2,3,3,2,1,0,0,0]`, let's use function, generating random values for now. +# !!! important +# because any function depending on `design` can be used, two things have to be taken care of: +# +# 1. in case a random component exist, specify a `RNG`, the basis might be evaluated multiple times inside `simulate` +# 2. a `maxlength` has to be specified via a tuple `(function.maxlength)`` +mybasisfun = design -> hanning.(generate_events(design).duration) +signal = LinearModelComponent(; + basis = (mybasisfun, 100), + formula = @formula(0 ~ 1 + category), + β = [1, 0.5], +); + +erp = UnfoldSim.simulate_component(MersenneTwister(1), signal, design); + + +# Finally, let's plot it, sorted by duration +plot_erpimage(erp, sortvalues = generate_events(design).duration) diff --git a/docs/literate/HowTo/newComponent.jl b/docs/literate/HowTo/newComponent.jl index bb8bee7..80ad9a2 100644 --- a/docs/literate/HowTo/newComponent.jl +++ b/docs/literate/HowTo/newComponent.jl @@ -1,6 +1,9 @@ # # New component: Duration + Shift # We want a new component that changes its duration and shift depending on a column in the event-design. This is somewhat already implemented in the HRF + Pupil bases +# !!! hint +# if you are just interested to use duration-dependency in your simulation, check out the component-function tutorial + using UnfoldSim using Unfold using Random diff --git a/docs/make.jl b/docs/make.jl index 35c5ff1..e539d98 100644 --- a/docs/make.jl +++ b/docs/make.jl @@ -45,8 +45,9 @@ makedocs(; ], "HowTo" => [ "Define a new, (imbalanced) design" => "./generated/HowTo/newDesign.md", + "Use a component-basis-function (duration)" => "./generated/HowTo/componentfunction.md", "Repeating a design" => "./generated/HowTo/repeatTrials.md", - "Define a new duration & jitter component" => "./generated/HowTo/newComponent.md", + "Define a new 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", diff --git a/src/component.jl b/src/component.jl index d3a0c5e..c9ba43c 100644 --- a/src/component.jl +++ b/src/component.jl @@ -52,11 +52,18 @@ LinearModelComponent(; ``` """ @with_kw struct LinearModelComponent <: AbstractComponent - basis::Any - formula::Any # e.g. 0~1+cond - left side must be "0" + basis::Union{Tuple{Function,Int},Array} + formula::FormulaTerm # e.g. 0~1+cond - left side must be "0" β::Vector contrasts::Dict = Dict() offset::Int = 0 + function LinearModelComponent(basis, formula, β, contrasts, offset) + @assert isa(basis, Tuple{Function,Int}) ".basis needs to be an `::Array` or a `Tuple(function::Function,maxlength::Int)`" + @assert basis[2] > 0 "`maxlength` needs to be longer than 0" + new(basis, formula, β, contrasts, offset) + end + LinearModelComponent(basis::Array, formula, β, contrasts, offset) = + new(basis, formula, β, contrasts, offset) end @@ -143,9 +150,51 @@ function simulate_component(rng, c::MultichannelComponent, design::AbstractDesig return reshape(y_proj, length(c.projection), size(y)...) end +""" + basis(c::AbstractComponent) + +returns the basis of the component (typically `c.basis`) +""" +basis(c::AbstractComponent) = c.basis + +""" + basis(c::AbstractComponent,design) +evaluates the basis, if basis is a vector, directly returns it. if basis is a tuple `(f::Function,maxlength::Int)`, evaluates the function with input `design`. Cuts the resulting vector or Matrix at `maxlength` +""" +basis(c::AbstractComponent, design) = basis(basis(c), design) + + +basis(b::AbstractVector, design) = b +function basis(basis::Tuple{Function,Int}, design) + f = basis[1] + maxlength = basis[2] + basis_out = f(design) + if isa(basis_out, AbstractVector{<:AbstractVector}) || isa(basis_out, AbstractMatrix) + @assert length(basis_out) == size(design)[1] "Component basis function needs to either return a Vector, a Matrix with dim(2) == size(design,1), or a Vector of Vectors with length(b) == size(design,1)" + end + limit_basis(basis_out, maxlength) +end + + +function limit_basis(b::AbstractVector{<:AbstractVector}, maxlength) + + # first cut off maxlength + b = limit_basis.(b, maxlength) + # now fill up with 0's + Δlengths = maxlength .- length.(b) + + b = pad_array.(b, Δlengths, 0) + basis_out = reduce(hcat, b) + + + return basis_out +end +limit_basis(b::AbstractVector{<:Number}, maxlength) = b[1:min(length(b), maxlength)] +limit_basis(b::AbstractMatrix, maxlength) = b[1:min(length(b), maxlength), :] + +Base.length(c::AbstractComponent) = isa(basis(c), Tuple) ? basis(c)[2] : length(basis(c)) -Base.length(c::AbstractComponent) = length(c.basis) """ maxlength(c::Vector{<:AbstractComponent}) = maximum(length.(c)) @@ -176,7 +225,7 @@ function simulate_component(rng, c::LinearModelComponent, design::AbstractDesign X = generate_designmatrix(c.formula, events, c.contrasts) y = X * c.β - return y' .* c.basis + return y' .* basis(c, design) end From 0c3a11dfdb9fb3b5d51999069e154ae9232b9d22 Mon Sep 17 00:00:00 2001 From: "behinger (s-ccs 001)" Date: Sun, 10 Mar 2024 19:14:17 +0000 Subject: [PATCH 3/5] bugfix with endless loop due to multiple dispatch --- src/simulation.jl | 10 ++++++++-- 1 file changed, 8 insertions(+), 2 deletions(-) diff --git a/src/simulation.jl b/src/simulation.jl index 4392424..4d8ea43 100755 --- a/src/simulation.jl +++ b/src/simulation.jl @@ -35,9 +35,15 @@ Some remarks to how the noise is added: """ -function simulate(args...; kwargs...) +function simulate( + design::AbstractDesign, + signal, + onset::AbstractOnset, + noise::AbstractNoise; + 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), args...; kwargs...) + simulate(MersenneTwister(1), design, signal, onset, noise; kwargs...) end function simulate( From de8c2a8990bc2da394da4677839f69970d66889c Mon Sep 17 00:00:00 2001 From: "behinger (s-ccs 001)" Date: Sun, 10 Mar 2024 19:30:36 +0000 Subject: [PATCH 4/5] function component for multi-subject --- src/component.jl | 18 +++++++++++++++--- 1 file changed, 15 insertions(+), 3 deletions(-) diff --git a/src/component.jl b/src/component.jl index c9ba43c..9ecef25 100644 --- a/src/component.jl +++ b/src/component.jl @@ -170,7 +170,12 @@ function basis(basis::Tuple{Function,Int}, design) maxlength = basis[2] basis_out = f(design) if isa(basis_out, AbstractVector{<:AbstractVector}) || isa(basis_out, AbstractMatrix) - @assert length(basis_out) == size(design)[1] "Component basis function needs to either return a Vector, a Matrix with dim(2) == size(design,1), or a Vector of Vectors with length(b) == size(design,1)" + if isa(basis_out, AbstractMatrix) + l = size(basis_out, 2) + else + l = length(basis_out) # vector of vector case + end + @assert l == size(generate_events(design))[1] "Component basis function needs to either return a Vector of vectors or a Matrix with dim(2) == size(design,1) [l / $(size(design,1))], or a Vector of Vectors with length(b) == size(design,1) [$l / $(size(design,1))]. " end limit_basis(basis_out, maxlength) end @@ -266,7 +271,7 @@ julia> simulate_component(StableRNG(1),c,design) function simulate_component( rng, c::MixedModelComponent, - design::AbstractDesign, + design::AbstractDesign; return_parameters = false, ) events = generate_events(design) @@ -305,8 +310,15 @@ function simulate_component( rethrow(e) end + @debug size(basis(c, design)) # in case the parameters are of interest, we will return those, not them weighted by basis - epoch_data_component = kron(return_parameters ? [1.0] : c.basis, m.y') + b = return_parameters ? [1.0] : basis(c, design) + @debug :b, typeof(b), size(b), :m, size(m.y') + if isa(b, AbstractMatrix) + epoch_data_component = ((m.y' .* b)) + else + epoch_data_component = kron(b, m.y') + end return epoch_data_component #= else From 57a7f68eb10a992d8fc8fb5d21dc7c62a8db225e Mon Sep 17 00:00:00 2001 From: "behinger (s-ccs 001)" Date: Mon, 11 Mar 2024 09:29:15 +0000 Subject: [PATCH 5/5] forgot to define offset in LinearModelFunction --- src/predefinedSimulations.jl | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/src/predefinedSimulations.jl b/src/predefinedSimulations.jl index 5a1c2a6..85292ff 100644 --- a/src/predefinedSimulations.jl +++ b/src/predefinedSimulations.jl @@ -46,9 +46,9 @@ function predef_eeg( # component / signal sfreq = 100, - p1 = (p100(; sfreq = sfreq), @formula(0 ~ 1), [5], Dict()), - n1 = (n170(; sfreq = sfreq), @formula(0 ~ 1 + condition), [5, -3], Dict()), - p3 = (p300(; sfreq = sfreq), @formula(0 ~ 1 + continuous), [5, 1], Dict()), + p1 = (p100(; sfreq = sfreq), @formula(0 ~ 1), [5], Dict(), 0), + n1 = (n170(; sfreq = sfreq), @formula(0 ~ 1 + condition), [5, -3], Dict(), 0), + p3 = (p300(; sfreq = sfreq), @formula(0 ~ 1 + continuous), [5, 1], Dict(), 0), kwargs..., )