Merge branch 'componentfunction' into sequFormulaOnset

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)

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

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",

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,56 @@ 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)
+        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
+
+
+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 +230,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
 
 
@@ -217,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)
@@ -256,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

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

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...,
 )
 

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(