From 59f53643a4d91f1b15f00ca05f8bc4ab39145222 Mon Sep 17 00:00:00 2001 From: "behinger (s-ccs 001)" Date: Fri, 20 Sep 2024 13:35:36 +0000 Subject: [PATCH 1/5] initial version of TRF simulation --- docs/literate/tutorials/TRF.jl | 29 +++++++++++++++++++++++++++++ docs/make.jl | 1 + src/UnfoldSim.jl | 4 ++-- src/component.jl | 28 ++++++++++++++++++++++++++++ src/onset.jl | 15 +++++++++++++++ 5 files changed, 75 insertions(+), 2 deletions(-) create mode 100644 docs/literate/tutorials/TRF.jl diff --git a/docs/literate/tutorials/TRF.jl b/docs/literate/tutorials/TRF.jl new file mode 100644 index 00000000..8ee07f41 --- /dev/null +++ b/docs/literate/tutorials/TRF.jl @@ -0,0 +1,29 @@ +using UnfoldSim +using CairoMakie +using Random + +# # Temporal response functions +# So far, we simulated event-related potentials. In this tutorial, we will switch gears a bit, and simulate a continuous response to a continuous feature vector. +# +# One good example is the visual response to a grayscale circle, which continuously changes its luminance (this goes back to VESPA, described in the ground breaking [Lalor 2006 paper](10.1016/j.neuroimage.2006.05.054). The brain will react to this luminance change continuously. TRFs typically describe the process to recover the brain response (in terms of a filter response). Here we set out to simulate such a dataset first and foremost. +# +# ## Simulation +# we start with the simplest possible design, one condition +design = SingleSubjectDesign(conditions=Dict(:dummy=>["A"])); + +# next we define the convolution kernel that the feature signal should be convolved with (= the brain response == the TRF) +brain_response = [LinearModelComponent(basis=p100(),formula=@formula(0~1),β=[1]), + LinearModelComponent(basis=n170(),formula=@formula(0~1),β=[1]), + LinearModelComponent(basis=p300(),formula=@formula(0~1),β=[1])]; + +# !!! hint +# For a fancy way to write the above code, you could use `LinearModelComponent.([p100,n170,p300],@formula(0~1),Ref([1]))`, notice the `.(...)` indicating broadcasting + +# Now we can simulate our feature signal, 10_000 samples of random gray scale values +feature = rand(10_000) + +# Next we have to nest the response in a `TRFComponent` and add ou +trf_component = TRFComponent(brain_response,feature) + +# Finally, when simulating, we have only a single "event" (that is, TRF-) onset, the first sample. Therefore, we use `TRFOnset` to indicate this. +dat,evts = simulate(design,trf_component,UnfoldSim.TRFOnset()) \ No newline at end of file diff --git a/docs/make.jl b/docs/make.jl index 2b24ff62..bf25d27c 100644 --- a/docs/make.jl +++ b/docs/make.jl @@ -37,6 +37,7 @@ makedocs(; "Simulate event-related potentials (ERPs)" => "generated/tutorials/simulateERP.md", "Power analysis" => "generated/tutorials/poweranalysis.md", "Multi-subject simulation" => "generated/tutorials/multisubject.md", + "Temporal Response Function" => "generated/tutorials/TRF.md", ], "Reference" => [ "Overview of functionality" => "./generated/reference/overview.md", diff --git a/src/UnfoldSim.jl b/src/UnfoldSim.jl index cd1b0287..a833ad80 100644 --- a/src/UnfoldSim.jl +++ b/src/UnfoldSim.jl @@ -43,7 +43,7 @@ export create_re export Simulation # component types -export MixedModelComponent, LinearModelComponent +export MixedModelComponent, LinearModelComponent, TRFComponent # export designs export MultiSubjectDesign, SingleSubjectDesign, RepeatDesign @@ -64,7 +64,7 @@ export simulate, export pad_array, convert # export Offsets -export UniformOnset, LogNormalOnset, NoOnset +export UniformOnset, LogNormalOnset, NoOnset, TRFOnset # re-export StatsModels export DummyCoding, EffectsCoding diff --git a/src/component.jl b/src/component.jl index 6c94e8b1..8a4b6de2 100644 --- a/src/component.jl +++ b/src/component.jl @@ -285,6 +285,34 @@ function weight_σs(σs::Dict, b_σs::Float64, σ_lmm::Float64) return named_random_effects end + +#--- TRF Component + +""" + TRFComponent(components::Vector{<:AbstractComponents},features::AbstractArray) + + This component can be used to convolve a `response` of a component-vector with a feature-array. + +Each column of the feature-array will be convolved with one generated response from the AbstractComponent-vector (that is, each row of the respective `AbstractDesign`) + +If only a single TRF is needed, a vector can be provided. +""" +struct TRFComponent <: AbstractComponent + components + features::AbstractArray +end +UnfoldSim.length(t::TRFComponent) = size(t.features,1) + + +function UnfoldSim.simulate_component(rng,c::TRFComponent,design::AbstractDesign) + kernel = UnfoldSim.simulate_responses(rng,c.components,UnfoldSim.Simulation(design,c,NoOnset(),NoNoise())) + + @assert size(kernel,2) == size(c.features,2) "if the $(typeof(design)) generates multiple columns (sz=$(size(kernel))), the $(typeof(c)) needs to have multiple columns as well (sz=$(size(c.signal)))" + x = reduce(hcat,UnfoldSim.conv.(eachcol(c.features),eachcol(kernel)))[1:size(c.features,1),:] + return x + +end + #---- """ diff --git a/src/onset.jl b/src/onset.jl index 4ad97aa7..fb693389 100644 --- a/src/onset.jl +++ b/src/onset.jl @@ -31,6 +31,21 @@ In the case that the user directly wants no overlap to be simulated (=> epoched """ struct NoOnset <: AbstractOnset end + +""" + struct TRFOnset <: AbstractOnset end +In the case that a TRF is simulated (via `TRFComponent`). This onset returns a vector of zero-latencies, indicating that the TRF starts at the beginning of the signal. +""" +struct TRFOnset <: AbstractOnset end + +function UnfoldSim.simulate_interonset_distances(rng,onset::TRFOnset,design) + sz = size(design) + return Int.(zeros(sz)) +end + + + + """ simulate_interonset_distances(rng, onset::UniformOnset, design::AbstractDesign) simulate_interonset_distances(rng, onset::LogNormalOnset, design::AbstractDesign) From 4f93a5471b6644ba7bd119f54c12ed97648aed35 Mon Sep 17 00:00:00 2001 From: "behinger (s-ccs 001)" Date: Fri, 20 Sep 2024 14:10:26 +0000 Subject: [PATCH 2/5] improved tutorial --- docs/literate/tutorials/TRF.jl | 58 ++++++++++++++++++++++++++++++++-- src/onset.jl | 12 +++---- 2 files changed, 61 insertions(+), 9 deletions(-) diff --git a/docs/literate/tutorials/TRF.jl b/docs/literate/tutorials/TRF.jl index 8ee07f41..5ab44d60 100644 --- a/docs/literate/tutorials/TRF.jl +++ b/docs/literate/tutorials/TRF.jl @@ -20,10 +20,62 @@ brain_response = [LinearModelComponent(basis=p100(),formula=@formula(0~1),β=[1] # For a fancy way to write the above code, you could use `LinearModelComponent.([p100,n170,p300],@formula(0~1),Ref([1]))`, notice the `.(...)` indicating broadcasting # Now we can simulate our feature signal, 10_000 samples of random gray scale values -feature = rand(10_000) +feature = rand(1_000) # Next we have to nest the response in a `TRFComponent` and add ou -trf_component = TRFComponent(brain_response,feature) +trf_component = TRFComponent(brain_response,feature); # Finally, when simulating, we have only a single "event" (that is, TRF-) onset, the first sample. Therefore, we use `TRFOnset` to indicate this. -dat,evts = simulate(design,trf_component,UnfoldSim.TRFOnset()) \ No newline at end of file +dat,evts = simulate(design,trf_component,UnfoldSim.TRFOnset()); + +# Let's plot the feature signal and the TRF response +f,ax,h =lines(dat) +lines!(feature) +f + +# ## Multivariate TRFs +# Now TRFs can depend on multiple variables e.g. the luminance and the size of the circle. +feature_luminance = rand(1_000) +feature_size = rand(1_000) + +# We could call the `simulate` function twice and simply add the results: +dat_l,_ = simulate(design,TRFComponent(brain_response,feature_luminance),UnfoldSim.TRFOnset()); +dat_r,_ = simulate(design,TRFComponent(brain_response,feature_size),UnfoldSim.TRFOnset()); + +# let's plot and compare to the previous plot +f,ax,h = lines(dat_l .+ dat_r) +lines!(dat) +f +# as visible, the blue line (luminance+size) has ~double the amplitude. This is because we simulated two brain responses and simply added them. + +# A bit more convenient way is to do the following +dat_combined,_ = simulate(design,[TRFComponent(brain_response,feature_size),TRFComponent(brain_response,feature_luminance)],UnfoldSim.TRFOnset()); +f,ax,h = lines(dat_combined) +lines!(dat_l .+ dat_r) +f +# where you can see that the results are equivalent. + +# ## Another cool feature is to modulate the feature vector based on the design +design_mod = SingleSubjectDesign(conditions=Dict(:condition=>["A","B"])); + +# Let's take only a single component for simplicity. Note how the `formula` has been changed. The β allows to control the amplitude. In this linear model component, the default contrast-function is `Dummy` (also known as `Reference` coding), which means, the second beta indicated a "difference" +brain_response_mod = LinearModelComponent(basis=p100(),formula=@formula(0~1+condition),β=[1,1]); + +# let's simulate another feature signal, but this time, we simulate a Matrix. +feature_mod = rand(1000,2) +feature_mod[:,2] .= 0 +feature_mod[500:600,2] .= 1; + +# to better understand how our (experimental) feature now looks like, let's visualize it +series(feature_mod') +# As visible, the first column has again a random signal, indicating e.g. luminance changes. The second temporal feature indicates some offset (a colorchange?) between 500 and 600 samples. + + +dat_mod,_ = simulate(design_mod,TRFComponent([brain_response_mod],feature_mod),UnfoldSim.TRFOnset()); +lines(dat_mod) + +# !!! hint +# Excourse: now you rightfully might ask why the jump is to >10 a.u.? The reason is that you are effectively convolving a feature-signal above 1 (feature_mod * (β_0 + β_1) = 1 * (1+1)), and with a kernel with maximum = 1, these 2's add up. The kernel has a rough width of ~5 which results in additional 5*2 => 10 per affected sample + +# ## Combination with a event-related design +# Right now there is no easy interface to do this. You have to simulate a TRF signal, and an rERP signal and then add the two signals. diff --git a/src/onset.jl b/src/onset.jl index e674ff98..4c03fd5d 100644 --- a/src/onset.jl +++ b/src/onset.jl @@ -42,7 +42,7 @@ In the case that a TRF is simulated (via `TRFComponent`). This onset returns a v """ struct TRFOnset <: AbstractOnset end -function UnfoldSim.simulate_interonset_distances(rng,onset::TRFOnset,design) +function simulate_interonset_distances(rng,onset::TRFOnset,design) sz = size(design) return Int.(zeros(sz)) end @@ -150,10 +150,10 @@ end function simulate_interonset_distances(rng, o::UniformOnsetFormula, design::AbstractDesign) events = generate_events(design) widths = - UnfoldSim.generate_designmatrix(o.width_formula, events, o.width_contrasts) * + generate_designmatrix(o.width_formula, events, o.width_contrasts) * o.width_β offsets = - UnfoldSim.generate_designmatrix(o.offset_formula, events, o.offset_contrasts) * + generate_designmatrix(o.offset_formula, events, o.offset_contrasts) * o.offset_β return Int.( @@ -208,10 +208,10 @@ function simulate_interonset_distances( events = generate_events(design) - μs = UnfoldSim.generate_designmatrix(o.μ_formula, events, o.μ_contrasts) * o.μ_β - σs = UnfoldSim.generate_designmatrix(o.σ_formula, events, o.σ_contrasts) * o.σ_β + μs = generate_designmatrix(o.μ_formula, events, o.μ_contrasts) * o.μ_β + σs = generate_designmatrix(o.σ_formula, events, o.σ_contrasts) * o.σ_β offsets = - UnfoldSim.generate_designmatrix(o.offset_formula, events, o.offset_contrasts) * + generate_designmatrix(o.offset_formula, events, o.offset_contrasts) * o.offset_β From 26ad7fb7405f46fe3bc6f70aee5929c427e3ce33 Mon Sep 17 00:00:00 2001 From: "behinger (s-ccs 001)" Date: Fri, 20 Sep 2024 14:39:08 +0000 Subject: [PATCH 3/5] fix doilink --- docs/literate/tutorials/TRF.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/docs/literate/tutorials/TRF.jl b/docs/literate/tutorials/TRF.jl index 5ab44d60..0d2fbc49 100644 --- a/docs/literate/tutorials/TRF.jl +++ b/docs/literate/tutorials/TRF.jl @@ -5,7 +5,7 @@ using Random # # Temporal response functions # So far, we simulated event-related potentials. In this tutorial, we will switch gears a bit, and simulate a continuous response to a continuous feature vector. # -# One good example is the visual response to a grayscale circle, which continuously changes its luminance (this goes back to VESPA, described in the ground breaking [Lalor 2006 paper](10.1016/j.neuroimage.2006.05.054). The brain will react to this luminance change continuously. TRFs typically describe the process to recover the brain response (in terms of a filter response). Here we set out to simulate such a dataset first and foremost. +# One good example is the visual response to a grayscale circle, which continuously changes its luminance (this goes back to VESPA, described in the ground breaking [Lalor 2006 paper](https://doi.org/10.1016/j.neuroimage.2006.05.054). The brain will react to this luminance change continuously. TRFs typically describe the process to recover the brain response (in terms of a filter response). Here we set out to simulate such a dataset first and foremost. # # ## Simulation # we start with the simplest possible design, one condition From 8ea515b968bacc94461715107aaa9eccffd043ef Mon Sep 17 00:00:00 2001 From: "behinger (s-ccs 001)" Date: Fri, 20 Sep 2024 14:41:22 +0000 Subject: [PATCH 4/5] maybe fix the interonset problem? --- src/onset.jl | 16 ++++++++++------ 1 file changed, 10 insertions(+), 6 deletions(-) diff --git a/src/onset.jl b/src/onset.jl index 4c03fd5d..edfa64a0 100644 --- a/src/onset.jl +++ b/src/onset.jl @@ -42,12 +42,6 @@ In the case that a TRF is simulated (via `TRFComponent`). This onset returns a v """ struct TRFOnset <: AbstractOnset end -function simulate_interonset_distances(rng,onset::TRFOnset,design) - sz = size(design) - return Int.(zeros(sz)) -end - - """ @@ -74,6 +68,16 @@ function simulate_interonset_distances(rng, onset::LogNormalOnset, design::Abstr end +""" + simulate_interonset_distances(rng, onset::UniformOnset, design::AbstractDesign) +In the case that a TRF is simulated (via `TRFComponent`). This onset returns a vector of zero-latencies, indicating that the TRF starts at the beginning of the signal. +""" +function simulate_interonset_distances(rng,onset::TRFOnset,design) + sz = size(design) + return Int.(zeros(sz)) +end + + #function simulate_interonset_distances(rng, onset::AbstractOnset,design::) From 73425faae5aeaa9b9f3adc93d8671fdfd9d84868 Mon Sep 17 00:00:00 2001 From: Benedikt Ehinger Date: Mon, 4 Nov 2024 11:49:28 +0100 Subject: [PATCH 5/5] Apply suggestions from code review Co-authored-by: github-actions[bot] <41898282+github-actions[bot]@users.noreply.github.com> --- src/UnfoldSim.jl | 3 ++- src/component.jl | 21 ++++++++++++++------- src/onset.jl | 6 +++--- 3 files changed, 19 insertions(+), 11 deletions(-) diff --git a/src/UnfoldSim.jl b/src/UnfoldSim.jl index 2f763f69..c545678d 100644 --- a/src/UnfoldSim.jl +++ b/src/UnfoldSim.jl @@ -67,7 +67,8 @@ export pad_array, convert # export Offsets -export UniformOnset, LogNormalOnset, NoOnset, UniformOnsetFormula, LogNormalOnsetFormula, TRFOnset +export UniformOnset, + LogNormalOnset, NoOnset, UniformOnsetFormula, LogNormalOnsetFormula, TRFOnset # re-export StatsModels export DummyCoding, EffectsCoding diff --git a/src/component.jl b/src/component.jl index a2e5ead7..c65d7924 100644 --- a/src/component.jl +++ b/src/component.jl @@ -411,17 +411,24 @@ Each column of the feature-array will be convolved with one generated response f If only a single TRF is needed, a vector can be provided. """ struct TRFComponent <: AbstractComponent - components + components::Any features::AbstractArray end -UnfoldSim.length(t::TRFComponent) = size(t.features,1) +UnfoldSim.length(t::TRFComponent) = size(t.features, 1) -function UnfoldSim.simulate_component(rng,c::TRFComponent,design::AbstractDesign) - kernel = UnfoldSim.simulate_responses(rng,c.components,UnfoldSim.Simulation(design,c,NoOnset(),NoNoise())) - - @assert size(kernel,2) == size(c.features,2) "if the $(typeof(design)) generates multiple columns (sz=$(size(kernel))), the $(typeof(c)) needs to have multiple columns as well (sz=$(size(c.signal)))" - x = reduce(hcat,UnfoldSim.conv.(eachcol(c.features),eachcol(kernel)))[1:size(c.features,1),:] +function UnfoldSim.simulate_component(rng, c::TRFComponent, design::AbstractDesign) + kernel = UnfoldSim.simulate_responses( + rng, + c.components, + UnfoldSim.Simulation(design, c, NoOnset(), NoNoise()), + ) + + @assert size(kernel, 2) == size(c.features, 2) "if the $(typeof(design)) generates multiple columns (sz=$(size(kernel))), the $(typeof(c)) needs to have multiple columns as well (sz=$(size(c.signal)))" + x = reduce(hcat, UnfoldSim.conv.(eachcol(c.features), eachcol(kernel)))[ + 1:size(c.features, 1), + :, + ] return x end diff --git a/src/onset.jl b/src/onset.jl index edfa64a0..31f6c735 100644 --- a/src/onset.jl +++ b/src/onset.jl @@ -72,9 +72,9 @@ end simulate_interonset_distances(rng, onset::UniformOnset, design::AbstractDesign) In the case that a TRF is simulated (via `TRFComponent`). This onset returns a vector of zero-latencies, indicating that the TRF starts at the beginning of the signal. """ -function simulate_interonset_distances(rng,onset::TRFOnset,design) - sz = size(design) - return Int.(zeros(sz)) +function simulate_interonset_distances(rng, onset::TRFOnset, design) + sz = size(design) + return Int.(zeros(sz)) end