diff --git a/.github/workflows/CI.yml b/.github/workflows/CI.yml
index 6788ea41..546b8452 100644
--- a/.github/workflows/CI.yml
+++ b/.github/workflows/CI.yml
@@ -5,8 +5,6 @@ on:
- main
tags: '*'
pull_request:
- branches:
- - main
concurrency:
# Skip intermediate builds: always.
# Cancel intermediate builds: only if it is a pull request build.
@@ -59,3 +57,10 @@ jobs:
using UnfoldSim
DocMeta.setdocmeta!(UnfoldSim, :DocTestSetup, :(using UnfoldSim); recursive=true)
doctest(UnfoldSim)'
+ - name: Update UnfoldDocs
+ if: github.ref == 'refs/heads/main'
+ uses: peter-evans/repository-dispatch@v3
+ with:
+ token: ${{ secrets.UNFOLDDOCSPAT }}
+ repository: unfoldtoolbox/UnfoldDocs
+ event-type: individual-docs-updated
diff --git a/CITATION.cff b/CITATION.cff
index 9a0df5cc..22070ecc 100644
--- a/CITATION.cff
+++ b/CITATION.cff
@@ -2,20 +2,52 @@
# Visit https://bit.ly/cffinit to generate yours today!
cff-version: 1.2.0
-title: UnfoldSim.jl
-message: 'If you use this page, please cite it as below.'
+title: >-
+ UnfoldSim.jl: Simulating continuous event-based time
+ series data for EEG and beyond
+message: "If you use this package, please cite it as below."
type: software
authors:
- given-names: Judith
family-names: Schepers
+ orcid: "https://orcid.org/0009-0000-9270-730X"
- given-names: Luis
family-names: Lips
- given-names: Maanik
family-names: Marathe
- given-names: Benedikt
family-names: Ehinger
- orcid: 'https://orcid.org/0000-0002-6276-3332'
+ orcid: "https://orcid.org/0000-0002-6276-3332"
identifiers:
- type: doi
value: 10.5281/zenodo.7738651
-repository-code: 'https://github.com/unfoldtoolbox/UnfoldSim.jl'
+ description: Zenodo
+ - type: doi
+ value: 10.21105/joss.06641
+ description: JOSS paper
+repository-code: "https://github.com/unfoldtoolbox/UnfoldSim.jl"
+preferred-citation:
+ authors:
+ - family-names: Schepers
+ given-names: Judith
+ orcid: "https://orcid.org/0009-0000-9270-730X"
+ - family-names: Lips
+ given-names: Luis
+ - family-names: Marathe
+ given-names: Maanik
+ - family-names: Ehinger
+ given-names: Benedikt V.
+ orcid: "https://orcid.org/0000-0002-6276-3332"
+ date-published: "2025-03-14"
+ doi: 10.21105/joss.06641
+ issn: 2475-9066
+ issue: 107
+ journal: Journal of Open Source Software
+ publisher:
+ name: Open Journals
+ start: 6641
+ title: "UnfoldSim.jl: Simulating continuous event-based time series
+ data for EEG and beyond"
+ type: article
+ url: "https://joss.theoj.org/papers/10.21105/joss.06641"
+ volume: 10
diff --git a/Project.toml b/Project.toml
index 96d86d6b..05e4ba08 100644
--- a/Project.toml
+++ b/Project.toml
@@ -1,7 +1,7 @@
name = "UnfoldSim"
uuid = "ed8ae6d2-84d3-44c6-ab46-0baf21700804"
authors = ["Benedikt Ehinger", "Luis Lips", "Judith Schepers", "Maanik Marathe"]
-version = "0.4.0"
+version = "0.4.1"
[deps]
Artifacts = "56f22d72-fd6d-98f1-02f0-08ddc0907c33"
@@ -17,6 +17,7 @@ MixedModels = "ff71e718-51f3-5ec2-a782-8ffcbfa3c316"
MixedModelsSim = "d5ae56c5-23ca-4a1f-b505-9fc4796fc1fe"
Parameters = "d96e819e-fc66-5662-9728-84c9c7592b0a"
Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
+SequentialSamplingModels = "0e71a2a6-2b30-4447-8742-d083a85e82d1"
SignalAnalysis = "df1fea92-c066-49dd-8b36-eace3378ea47"
Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
StatsModels = "3eaba693-59b7-5ba5-a881-562e759f1c8d"
@@ -35,6 +36,7 @@ MixedModels = "4"
MixedModelsSim = "0.2"
Parameters = "0.12"
Random = "1"
+SequentialSamplingModels = "0.12"
SignalAnalysis = "0.4, 0.5,0.6,0.7,0.8,0.9,0.10"
Statistics = "1"
StatsModels = "0.6,0.7"
diff --git a/README.md b/README.md
index 82132f27..733dfae7 100644
--- a/README.md
+++ b/README.md
@@ -1,19 +1,20 @@
# [](https://github.com/unfoldtoolbox/UnfoldSim.jl/tree/main)
-[](https://unfoldtoolbox.github.io/UnfoldSim.jl/stable/)
-[](https://unfoldtoolbox.github.io/UnfoldSim.jl/dev/)
+[](https://unfoldtoolbox.github.io/UnfoldDocs/UnfoldSim.jl/stable/)
+[](https://unfoldtoolbox.github.io/UnfoldDocs/UnfoldSim.jl/dev/)
[](https://github.com/unfoldtoolbox/UnfoldSim.jl/actions/workflows/CI.yml?query=branch%3Amain)
[](https://codecov.io/gh/unfoldtoolbox/UnfoldSim.jl)
[](https://zenodo.org/badge/latestdoi/413455526)
+[](https://doi.org/10.21105/joss.06641)
-|rERP|EEG visualisation|EEG Simulations|BIDS pipeline|Decode EEG data|Statistical testing|
-|---|---|---|---|---|---|
-| 
| 
|
|
|
|
|
+|Estimation|Visualisation|Simulation|BIDS pipeline|Decoding|Statistics|MixedModelling|
+|---|---|---|---|---|---|---|
+| | |||||
|
A Julia package to simulate multivariate time series, e.g. model-based ERPs, fMRI activity, pupil dilation etc.
UnfoldSim.jl provides multi-channel support via EEG-forward models. Moreover, it is possible to simulate overlapping event-related activity and to add noise of a certain type e.g. Pink noise.
-Many tutorials, guides, how-tos and references are available in the [documentation](https://unfoldtoolbox.github.io/UnfoldSim.jl/)!
+Many tutorials, guides, how-tos and references are available in the [documentation](https://unfoldtoolbox.github.io/UnfoldDocs/UnfoldSim.jl/)!
@@ -79,7 +80,7 @@ data, events = simulate(
EEG researchers often analyze data containing (temporally) overlapping events (e.g. stimulus onset and button press, or consecutive eye-fixations), non-linear effects, and complex experimental designs. For a multitude of reasons, we often need to simulate such kinds of data: Simulated EEG data is useful to test preprocessing and analysis tools, validate statistical methods, illustrate conceptual issues, test toolbox functionalities, and find limitations of traditional analysis workflows. For instance, such simulation tools allow for testing the assumptions of new analysis algorithms and testing their robustness against any violation of these assumptions.
## Contributions
@@ -115,8 +116,38 @@ Please reach out, if you have contributed to UnfoldSim.jl but we have not listed
## Citation
-TBA
+If you use UnfoldSim.jl, please acknowledge and support our work by citing
+
+**1\. Our [JOSS paper](https://doi.org/10.21105/joss.06641):**
+
+```
+Schepers et al., (2025). UnfoldSim.jl: Simulating continuous event-based time series data for EEG and beyond.
+Journal of Open Source Software, 10(107), 6641, https://doi.org/10.21105/joss.06641
+```
+
+
+ BibTeX entry:
+
+```bib
+@article{Schepers2025,
+ doi = {10.21105/joss.06641},
+ url = {https://doi.org/10.21105/joss.06641},
+ year = {2025},
+ publisher = {The Open Journal},
+ volume = {10},
+ number = {107},
+ pages = {6641},
+ author = {Judith Schepers and Luis Lips and Maanik Marathe and Benedikt V. Ehinger},
+ title = {UnfoldSim.jl: Simulating continuous event-based time series data for EEG and beyond},
+ journal = {Journal of Open Source Software}
+ }
+```
+
+
+and
+
+**2\. The corresponding [Zenodo DOI](https://doi.org/10.5281/zenodo.7738651)** for the specific UnfoldSim.jl version that you are using in your work.
## Acknowledgements
-Funded by Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) under Germany´s Excellence Strategy – EXC 2075 – 390740016. Furthermore, the authors thank the International Max Planck Research School for Intelligent Systems (IMPRS-IS) for supporting Judith Schepers.
+Funded by Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) under Germany's Excellence Strategy – EXC 2075 – 390740016. Furthermore, the authors thank the International Max Planck Research School for Intelligent Systems (IMPRS-IS) for supporting Judith Schepers.
diff --git a/docs/Project.toml b/docs/Project.toml
index 7257e3f5..40ca016c 100644
--- a/docs/Project.toml
+++ b/docs/Project.toml
@@ -10,6 +10,7 @@ Literate = "98b081ad-f1c9-55d3-8b20-4c87d4299306"
LiveServer = "16fef848-5104-11e9-1b77-fb7a48bbb589"
Parameters = "d96e819e-fc66-5662-9728-84c9c7592b0a"
Revise = "295af30f-e4ad-537b-8983-00126c2a3abe"
+SequentialSamplingModels = "0e71a2a6-2b30-4447-8742-d083a85e82d1"
StableRNGs = "860ef19b-820b-49d6-a774-d7a799459cd3"
StatsBase = "2913bbd2-ae8a-5f71-8c99-4fb6c76f3a91"
TimerOutputs = "a759f4b9-e2f1-59dc-863e-4aeb61b1ea8f"
diff --git a/docs/literate/tutorials/simulateDriftOverlap.jl b/docs/literate/tutorials/simulateDriftOverlap.jl
new file mode 100644
index 00000000..adb640ed
--- /dev/null
+++ b/docs/literate/tutorials/simulateDriftOverlap.jl
@@ -0,0 +1,152 @@
+# # Simulate an Evidence Accumulation Overlap and Deconvolution
+
+# We want to use the SequenceDesign with three components to simulate an full evidence accumulation process with a Overlap between the evidence accumulation and the response signal. Afterwards we use the Deconvolution Method from Unfold to unfold the signals again.
+
+
+
+# ### Setup
+# ```@raw html
+#
+# Click to expand
+# ```
+using UnfoldSim
+using Unfold
+using StableRNGs
+using CairoMakie, UnfoldMakie
+using SequentialSamplingModels
+
+fs = 500
+Δt = 1 / fs; # time step
+tEnd = 1.0 # trial Duration
+time_vec = 0:Δt:tEnd; # time base - let's make it a typical stimulus duration
+max_length = tEnd / Δt
+
+c_stimulus = LinearModelComponent(;
+ basis = ones(Int(0.1 * fs)),
+ formula = @formula(0 ~ 1 + condition),
+ β = [1, 0],
+)
+c_response = LinearModelComponent(;
+ basis = ones(Int(0.2 * fs)),
+ formula = @formula(0 ~ 1 + condition),
+ β = [-0.5, 0],
+)
+
+# ```@raw html
+#
+# ```
+# ## Design
+# Let's generate a single design with two drift_rates
+design_single =
+ SingleSubjectDesign(conditions = Dict(:drift_rate => [1, 2], :condition => [1]))
+design_seq = SequenceDesign(design_single, "SCR_")
+design_rep = RepeatDesign(design_seq, 100) # 100 SCR-repeats
+
+# # Simplified S-C-R drift-model with LBA
+# ## LBA
+# The following parameters are described in the [SequentialSamplingModels.LBA model documentation](https://itsdfish.github.io/SequentialSamplingModels.jl/dev/lba/)
+v = [0.8] ## Drift Rate, needs to be in [ ] for compatability with `SequentialSamplingModels.jl``
+A = 0.2 ## The starting point is sampled uniformly between `[0, A]`
+k = 0.4 ## Difference of A with the bound/threshold ($A + k = b$), where b is the decision threshold.
+
+lba_parameter = Dict(:ν => v, :A => A, :k => k, :τ => 0.0) # we didnt specify σ thus it is 1 by default,
+
+# !!! important
+# The non-decision time has to be specified here as "τ=>0." because it is not defined whether non-decision time applies to stimulus encoding or response execution. If you have a non-0 value here, the events following the DriftComponent will be shifted by τ.
+drift = DriftComponent(max_length, fs, LBA, lba_parameter)
+
+# Let's have a quick look at the simulation
+series(
+ hcat(simulate_component.(StableRNG.(1:10), Ref(drift), Ref(design_single))...)',
+ solid_color = (:black, 0.5),
+)
+# Super! Our start values are sampled between [0, 0.2] as indicated, and the boundary is at A+k = 0.2 + 0.4 = 0.6
+
+# ### Simulate EEG data
+components = Dict('S' => [c_stimulus], 'C' => [drift], 'R' => [c_response])
+
+# Next we define a SequenceOnset ingredient, this allows us fine-grained control on when events should happen.
+
+seq_onset = SequenceOnset(
+ Dict(
+ 'S' => UniformOnset(width = 0, offset = 0.2 * fs), # distance S<->C - we could choose 0.1*fs to immediately start after the stimulus component
+ 'C' => DriftOnset(), ## automatically determines which is the respective DriftComponent
+ 'R' => UniformOnset(width = 0, offset = 2 * fs), # distance R<->S
+ ),
+)
+
+# ## Simulate data with the created design, components and onsets
+data, evts = UnfoldSim.simulate(
+ StableRNG(12),
+ design_rep,
+ components,
+ seq_onset,
+ NoNoise(), # PinkNoise(noiselevel=1)
+)
+# ## Plot EEG to see the overlap as it appears
+f_data, _ax, _h = lines(data[1:5000])
+vlines!(evts.latency, color = Int.(evts.event), alpha = 0.3)
+xlims!(700, 3500)
+f_data
+
+# # Complex S-C-R drift-model ala `Kelly et al.`
+
+# In addition to changing the model to a more realistic version, we will also make use of the `drift_rate` specification in our `design``
+
+v = "drift_rate" ## take the drift rate directly from the design
+
+# The KellyModel has _many_ default parameters, let's grab them and immediately change a few
+
+kelly_model =
+ UnfoldSim.KellyModel(drift_rate = v, motor_delay = 0.4, sensor_encoding_delay = 0.2)
+
+kelly_model_parameter = UnfoldSim.create_kelly_parameters_dict(kelly_model)
+drift = DriftComponent(max_length, fs, UnfoldSim.KellyModel, kelly_model_parameter)
+components = Dict('S' => [c_stimulus], 'C' => [drift], 'R' => [c_response])
+
+# We are ready to simulate!
+data, evts = UnfoldSim.simulate(
+ StableRNG(12),
+ design_rep,
+ components,
+ seq_onset,
+ NoNoise(), # PinkNoise(noiselevel=1)
+)
+# ## Plot EEG to see the overlap as it appears
+f_data, _ax, _h = lines(data[1:5000])
+vlines!(evts.latency, color = Int.(evts.event), alpha = 0.3)
+xlims!(700, 3500)
+f_data
+
+# !!! note
+# There are varying gaps between the stimulus offset <-> driftstart and driftend <-> response start. This is because the `KellyModel` by default has a `sensor_encoding_delay` and a `motor_delay` (including a `motor_delay_variability` etc.)
+
+# # ERP & Deconvolution
+# Let's calculate traditional ERPs and deconvolved ones!
+evts.event = string.(evts.event) # due to a Unfold-bug we have to convert :S
+data_epochs, times_epoch =
+ Unfold.epoch(data = data, tbl = evts, τ = (-1.3, 1.3), sfreq = fs);
+f = @formula(0 ~ 1 + drift_rate)
+m = fit(UnfoldModel, ["S" => (f, times_epoch), "R" => (f, times_epoch)], evts, data_epochs);
+plot_erp(
+ effects(Dict(:drift_rate => [2, 3]), m);
+ mapping = (;
+ color = :drift_rate => nonnumeric,
+ col = :eventname => UnfoldMakie.sorter(["S", "R"]),
+ ),
+)
+
+# ## Deconvolution of the overlap
+fir_stim = firbasis(τ = (-0.1, 1.3), sfreq = fs)
+fir_resp = firbasis(τ = (-1, 0.3), sfreq = fs)
+
+m = fit(UnfoldModel, ["S" => (f, fir_stim), "R" => (f, fir_resp)], evts, data);
+plot_erp(
+ effects(Dict(:drift_rate => [2, 3]), m);
+ mapping = (;
+ color = :drift_rate => nonnumeric,
+ col = :eventname => UnfoldMakie.sorter(["S", "R"]),
+ ),
+)
+# Deconvolution sucessfully removed the stimulus and response overlaps from our estimates, but it kept the drift-diffusion aspect nearly untouched. It is important to stress the "nearly", because depending on the `delay` parameters, this can potentially lead to problematic situation where it is unclear whether the drift activity should be assigned to S or R. In most realistic situations we found to be safe.
+# Further note that the drift-diffusion amplitude can be reduced because it is split up between `S` and `R` events.`
\ No newline at end of file
diff --git a/docs/make.jl b/docs/make.jl
index f5a9c249..ee120eb3 100644
--- a/docs/make.jl
+++ b/docs/make.jl
@@ -40,6 +40,7 @@ makedocs(;
"Simulate event-related potentials (ERPs)" => "generated/tutorials/simulateERP.md",
"Power analysis" => "generated/tutorials/poweranalysis.md",
"Multi-subject simulation" => "generated/tutorials/multisubject.md",
+ "S-C-R & Evidence accumulation models" => "./generated/tutorials/simulateDriftOverlap.md",
],
"Reference" => [
"Overview of functionality" => "./generated/reference/overview.md",
@@ -57,6 +58,7 @@ makedocs(;
"Use existing experimental designs & onsets in the simulation" => "./generated/HowTo/predefinedData.md",
"Simulated marginal effects" => "./generated/HowTo/getGroundTruth.md",
"Sequence of events (e.g. SCR)" => "./generated/HowTo/sequence.md",
+ #"Simulate LBA evidence accumulation models" => "./generated/HowTo/driftComponent.md",
],
"Developer documentation" => "developer_docs.md",
"API / Docstrings" => "api.md",
diff --git a/src/UnfoldSim.jl b/src/UnfoldSim.jl
index 655c385d..3f35a362 100644
--- a/src/UnfoldSim.jl
+++ b/src/UnfoldSim.jl
@@ -17,6 +17,8 @@ using HDF5, Artifacts, FileIO
using Automa # for sequence
using LinearAlgebra # headmodel
+using SequentialSamplingModels # for SequentialSamplingModels
+
import DSP.hanning
import Base.length
import Base.size
@@ -32,6 +34,7 @@ include("headmodel.jl")
include("helper.jl")
include("sequence.jl")
include("bases.jl")
+include("sequentialSamplingModels.jl")
export size, length
export AbstractComponent, AbstractNoise, AbstractOnset, AbstractDesign
@@ -78,4 +81,7 @@ export AbstractHeadmodel, Hartmut, headmodel, leadfield, orientation, magnitude
# multichannel
export MultichannelComponent
+
+# evidence accumulation
+export DriftComponent, DriftOnset, SequenceOnset, KellyModel, LBA, DDM
end
diff --git a/src/component.jl b/src/component.jl
index f0b56d57..437f6a55 100644
--- a/src/component.jl
+++ b/src/component.jl
@@ -706,3 +706,128 @@ function simulate_and_add!(
@views epoch_data[:, 1+off:length(component)+off, :] .+=
simulate_component(rng, component, simulation)
end
+
+
+
+"""
+ DriftComponent <: AbstractComponent
+A component that adds an evidence accumulation trace, following a specified sequential sampling model.
+
+Works best with [`SequenceDesign`](@ref).
+
+# Fields
+- `max_length::Int`: maximum length of the simulated trace.
+- `sfreq::Real`: sample frequency used to simulate the trace, used to evaluate the parameters in the model_type.
+- `model_type`: Model struct which defines the model to use to generate the trace, e.g. `KellyModel`
+- `model_parameters`: Dict. Containing the parameters for the sequential sampling simulation model specified in `model_type`.
+
+# Examples
+```julia-repl
+# use the KellyModel and its default parameters to simulate traces from 0:1/500:1.0
+sfreq=500
+tEnd = 1.0
+max_length = tEnd*sfreq;
+model_parameter = Dict()
+DriftComponent(
+ max_length,
+ sfreq,
+ KellyModel,
+ model_parameter)
+```
+
+"""
+struct DriftComponent <: AbstractComponent
+ max_length::Int
+ sfreq::Real
+ model_type::Any
+ model_parameters::Dict
+end
+"""
+ simulate_component(rng, c::DriftComponent, design::AbstractDesign)
+
+Generate evidence accumulation traces based on the model specified in `c::DriftComponent` and (optionally) parameters from `design::AbtractDesign`
+
+# Returns
+- `Matrix{Float64}`: Simulated component for each event in the events data frame. The output dimensions are `length(c.time_vec) x size(events, 1)`.
+
+# Examples
+```julia-repl
+# use the KellyModel and its default parameters to simulate traces from 0:1/500:1.0
+julia> model_parameter = Dict();
+
+julia> c = DriftComponent(500, 500, KellyModel, model_parameter);
+
+julia> design_single = SingleSubjectDesign(conditions = Dict(:drift_rate => [0.5, 0.8], :condition => [1]));
+
+julia> design_seq = SequenceDesign(design_single,"SCR_");
+
+julia> simulate_component(StableRNG(1),c,design_seq)
+501x6 Matrix{Float64}:
+0.0 0.0 0.0 0.0 0.0 0.0
+0.0 0.0 0.0 0.0 0.0 0.0
+⋮ ⋮
+```
+"""
+function UnfoldSim.simulate_component(rng, c::DriftComponent, design::AbstractDesign)
+ _, traces = simulate_drift_component(deepcopy(rng), c, design)
+ return traces
+end
+"""
+ calculate_response_times_for_ssm(rng, component::DriftComponent, design::AbstractDesign)
+
+Generate response times of the evidence accumulation by using the model and its parameters specified in the component.
+
+Important: specify the same rng as in `simulate_component(rng, component::DriftComponent, design::AbstractDesign)` to ensure that the response times match the generated traces.
+
+# Arguments
+- `rng::AbstractRNG`: Random seed to ensure the same traces are created as in the use of the [`simulate_component`](@ref) function.
+- `component::DriftComponent`: Component to specify the model and its parameters to simulate the evidence accumulation.
+- `design::UnfoldSim.SubselectDesign`: Subselection of the Sequence Design to ensure we only generate rt for the drift component events.
+
+# Returns
+- `Vector{Float64}`: Simulated response times for each event in the events data frame. The output dimension is `size(events, 1)`.
+
+```julia-repl
+# use the KellyModel and its default parameters to simulate traces from 0:1/500:1.0
+julia> model_parameter = Dict();
+
+julia> c = DriftComponent(500, 500, KellyModel, model_parameter);
+
+julia> design_single = SingleSubjectDesign(conditions = Dict(:drift_rate => [0.5, 0.8], :condition => [1]));
+
+julia> design_seq = SequenceDesign(design_single,"SCR_");
+
+julia> calculate_response_times_for_ssm(StableRNG(1),c,design_seq)
+2-element Vector{Float64}:
+ 260.70134768436486
+ 360.1329203034039
+```
+"""
+function calculate_response_times_for_ssm(rng, component::DriftComponent, design::UnfoldSim.SubselectDesign)
+ rts, _ = simulate_drift_component(deepcopy(rng), component, design)
+ return rts
+end
+Base.length(c::DriftComponent) = c.max_length
+"""
+ get_model_parameter(rng, evt, d::Dict)
+
+Construct a parameter dictionary of parameter names and values.
+"""
+function get_model_parameter(rng, evt, d::Dict)
+ result_parameter = Dict()
+ for key in keys(d)
+ result_parameter[key] = get_model_parameter(rng, evt, d[key])
+ end
+ return result_parameter
+end
+"""
+ get_model_parameter(rng, evt, val)
+
+Recursive default call return the model parameter as it is.
+"""
+get_model_parameter(rng, evt, val) = val
+"""
+ get_model_parameter(rng, evt, val::String)
+Recursive call if the model parameter specified as string from the conditions in the events data frame.
+"""
+get_model_parameter(rng, evt, val::String) = evt[val]
diff --git a/src/onset.jl b/src/onset.jl
index 4d632a5b..1a105673 100644
--- a/src/onset.jl
+++ b/src/onset.jl
@@ -332,3 +332,155 @@ function simulate_interonset_distances(
#@debug reduce(hcat, rand.(deepcopy(rng), funs, 1))
return Int.(round.(offsets .+ reduce(vcat, rand.(deepcopy(rng), funs, 1))))
end
+"""
+ SequenceOnset <: AbstractOnset
+Allows to specify inter-onset-distance functions per event.
+
+All fields are mandatory. Works best with [`SequenceDesign`](@ref).
+
+# Fields
+- `onset::Dict`: for each Sequence event a Onset is defined.
+
+# Examples
+```julia-repl
+sequence_onset = SequenceOnset(
+ Dict('S'=>UniformOnset(width=0,offset=85*fs/100),
+ 'C'=>(DriftOnset(), UniformOnset(width=0, offset=150)),
+ 'R'=>UniformOnset(width=0,offset=120*fs/100)))
+```
+"""
+struct SequenceOnset <: AbstractOnset
+ onset::Dict
+end
+SequenceOnset(args::Pair...) = SequenceOnset(Dict(args...))
+"""
+ DriftOnset<: AbstractOnset
+A type that defines the onsets for a [`DriftComponent`](@ref).
+
+Works best with [`DriftComponent`](@ref).
+
+# Fields
+- `onset::Dict`: onset for the DriftComponent.
+
+# Examples
+```julia-repl
+drift_onset = DriftOnset()
+```
+"""
+struct DriftOnset{T} <: AbstractOnset
+ onset::T
+end
+DriftOnset() = DriftOnset(UniformOnset(width = 0, offset = 0))
+
+"""
+ UnfoldSim.simulate_interonset_distances(rng, onset::AbstractOnset, design::AbstractDesign, components::AbstractComponent)
+
+Base function to simulate Abstract Onset between components in an [`SequenceDesign`](@ref).
+
+# Arguments
+- `rng::StableRNG`: Random seed to ensure reproducibility.
+- `onset::AbstractOnset`: Onset of type AbstractOnset which defines how the onset is created.
+- `design::AbstractDesign`: Design for which the onsets are simulated.
+- `components::AbstractComponent`: The Component for which the onset is simulated.
+
+# Returns
+- `simulate_interonset_distances`: function call.
+"""
+UnfoldSim.simulate_interonset_distances(rng, onset::AbstractOnset, design::AbstractDesign, components::AbstractComponent) = UnfoldSim.simulate_interonset_distances(rng,onset,design)
+
+"""
+ UnfoldSim.simulate_interonset_distances(rng, onset::DriftOnset, design::AbstractDesign, components::AbstractComponent)
+
+Generates list of onsets for multiple [`DriftComponent`](@ref) in an [`SequenceDesign`](@ref), one for each sequence. Onsets are rounded, returned as `Int`
+
+# Arguments
+- `rng::StableRNG`: Random seed to ensure reproducibility.
+- `onset::DriftOnset`: DriftOnset defines to create onsets for a [`DriftComponent`](@ref).
+- `design::AbstractDesign`: Design for which the onsets are simulated.
+- `components::AbstractComponent`: The Component for which the onset is simulated.
+
+# Returns
+- `Vector{Int}`: the generated onsets for the drift components in the SequenceDesign.
+"""
+function UnfoldSim.simulate_interonset_distances(rng, onset::DriftOnset, design::AbstractDesign, components::AbstractComponent)
+ rts = calculate_response_times_for_ssm(deepcopy(rng), components, design)
+ return Int.(round.(rts))
+end
+
+"""
+ UnfoldSim.simulate_interonset_distances(rng, onset::Tuple{DriftOnset, UniformOnset}, design::AbstractDesign, components::AbstractComponent)
+
+Generates list of onsets for multiple [`DriftComponent`](@ref) in an [`SequenceDesign`](@ref) and possibility to ad an [`UniformOnset`](@ref).
+
+# Arguments
+- `rng::StableRNG`: Random seed to ensure reproducibility.
+- `onset::Tuple{DriftOnset, UniformOnset}`: DriftOnset defines to create onsets for a [`DriftComponent`](@ref) on top with an [`UniformOnset`](@ref).
+- `design::AbstractDesign`: Design for which the onsets are simulated.
+- `components::AbstractComponent`: The Component for which the onset is simulated.
+
+# Returns
+- `Vector{Float64}`: the generated onsets for the drift components in the SequenceDesign.
+"""
+function UnfoldSim.simulate_interonset_distances(rng, onset::Tuple, design::AbstractDesign, components::AbstractComponent)
+ return Int.(round.(reduce(.+, simulate_interonset_distances.(deepcopy(rng), onset, Ref(design), Ref(components)))))
+end
+
+"""
+ UnfoldSim.simulate_interonset_distances(rng, onset::Char, design::AbstractDesign, components::AbstractComponent)
+
+Generates list of onsets for the end of a sequence in an [`SequenceDesign`](@ref).
+
+# Arguments
+- `rng::StableRNG`: Random seed to ensure reproducibility.
+- `onset::Char`: Defines to simulates onsets at the end of a sequence.
+- `design::AbstractDesign`: Design for which the onsets are simulated.
+- `components::AbstractComponent`: The Component for which the onset is simulated.
+
+# Returns
+- `Vector{Float64}`: the generated onsets for the end of a sequence in the SequenceDesign.
+"""
+function UnfoldSim.simulate_interonset_distances(rng, onset::Char, design::AbstractDesign, components::AbstractComponent)
+ @assert onset == '_'
+ df = generate_events(rng, design)
+ nrows_df = Int(size(df, 1))
+ onsets = repeat([UnfoldSim.maxlength([components])], nrows_df)
+ return onsets
+end
+
+"""
+ UnfoldSim.simulate_onsets(rng, onset::SequenceOnset, simulation::Simulation)
+
+Generates list of onsets for all events of an [`SequenceDesign`](@ref), how to simulate the onsets is defined in the [`SequenceOnset`](@ref).
+
+# Arguments
+- `rng::StableRNG`: Random seed to ensure reproducibility.
+- `onset::SequenceOnset`: onset definition for each event in the sequence design.
+- `simulation::Simulation`: Simulation which contains the design and other elements for the experiment.
+
+# Returns
+- `Vector{Float64}`: the generated onsets for all events in the SequenceDesign.
+"""
+function UnfoldSim.simulate_onsets(rng, onset::SequenceOnset, simulation::Simulation)
+ @assert isa(simulation.design.design, SequenceDesign) "`SequenceOnset`is currently only compatible with a `SequenceDesign`"
+ events = generate_events(deepcopy(rng), simulation.design)
+ onset_map = Dict()
+ onset_counter = Dict()
+ for k in keys(onset.onset)
+ sub_design = UnfoldSim.SubselectDesign(simulation.design, k)
+ onsets_for_k = simulate_interonset_distances(deepcopy(rng), onset.onset[k], sub_design, simulation.components[k][1])
+ onset_map[k] = onsets_for_k
+ onset_counter[k] = 1
+ end
+ final_onsets = []
+ for (i, evt_k) in enumerate(events.event)
+ push!(final_onsets, onset_map[evt_k][onset_counter[evt_k]])
+ onset_counter[evt_k] += 1
+ end
+ final_onsets = vcat(final_onsets[end], final_onsets[1:end-1])
+ if maximum(final_onsets) > 10000
+ @warn "Number of simulated inter-event-distances was $(maximum(final_onsets)) - are you sure this is what you want?"
+ end
+ onsets_accum = accumulate(+, final_onsets, dims = 1, init = 1)
+ return onsets_accum
+
+end
diff --git a/src/sequentialSamplingModels.jl b/src/sequentialSamplingModels.jl
new file mode 100644
index 00000000..cc45d271
--- /dev/null
+++ b/src/sequentialSamplingModels.jl
@@ -0,0 +1,294 @@
+"""
+ KellyModel
+
+A advanced drift diffusion Model which can be used to simulate evidence accumulation.
+
+All fields can be named. Is used with [`DriftComponent`](@ref).
+The fields can be specified as a string, as this allows a reference to be made for selecting values from the design as parameters for the model used in the simulation.
+For example, different drift_rate values can be used depending on the design specification, which enables a kind of subject difference in the hole process.
+# Fields T::Union{Real,String}
+- `drift_rate::T`: defines the amount of evidence accumulated per time step. (roughly the steepness of the trace)
+- `sensor_encoding_delay::T`: constant event onset delay in seconds. (mimics sensory evidence)
+- `sensor_encoding_delay_variability::T`: variability (σ) in the delay of the event onset in seconds, added ontop of `event_onset`. (mimics sensory encoding delay)
+- `event_onset_distribution::Distribution{Univariate, Continuous}`: By default: Normal distribution using `sensor_encoding_delay` and `sensor_encoding_delay_variability` as μ and σ for sensor encoding delay. Any `Univariate` distribution from e.g. `Distributions.jl` can be used here.
+- `accumulative_level_noise::T`: sigma of the normal-distributed noise added to the accumulation process.
+- `boundary::T`: the threshold of evidence needed to make a decision. See also `urgency` for collapsing bounds
+- `start_point::T`: fixed delay between boundary reached and response time in seconds. (mimics motor time)
+- `start_point_variability::T`: variability in delay between boundary reached and response time in seconds. (mimics different reaction times of participants)
+- `start_point_distribution::Distribution{Univariate, Continuous}`: By default: Normal distribution using `start_point` and `start_point_variability` as μ and σ to add noise to the start-point. Any `Univariate` distribution from e.g. `Distributions.jl` can be used here. The default Normal distribution is truncated at the `boundary`
+- `motor_delay::T`: fixed delay between boundary reached and response time in seconds. (mimics motor time)
+- `motor_delay_variability::T`: variability in delay between boundary reached and response time in seconds. (mimics different reaction times of participants)
+- `motor_onset_distribution::Distribution{Univariate, Continuous}`: By default: Normal distribution using `motor_delay` and `motor_delay_variability` as μ and σ to add a motor delay. Any `Univariate` distribution from e.g. `Distributions.jl` can be used here.
+- `urgency::T`: fixed delay between boundary reached and response time in seconds. (mimics motor time)
+- `urgency_variability::T`: variability in delay between boundary reached and response time in seconds. (mimics different reaction times of participants)
+- `urgency_distribution::Distribution{Univariate, Continuous}`: By default: Normal distribution using `urgency` and `urgency_variability` as μ and σ for an urgency signal, truncated to be positive at 0 (see `post_accumulation_distribution`). Effectively one value is sampled per trial and used as a slope value times the time-vector. Any `Univariate` distribution from e.g. `Distributions.jl` can be used here.
+- `post_accumulation_duration::T`: fixed time the accumulation process resumes after boundary reached in seconds. (mimics evidence overshoot)
+- `post_accumulation_duration_variability::T`: variability in time the accumulation process resumes after boundary reached in seconds.
+- `post_accumulation_distribution::Distribution{Univariate, Continuous}`: By default: Normal distribution, trunacted to be positive at 0, using `post_accumulation_duration` and `post_accumulation_duration_variability` as μ (untruncated, effectively the mean is therefore higher and a function of both μ and σ, use e.g. `mean(truncated(Normal(100,10),lower=0))` to calculate the mean) and σ to sample a time per trial for how long the post-decision accumulation should take time. Effectively one value is sampled per trial and used as a slope value times the time-vector. Any `Univariate` distribution from e.g. `Distributions.jl` can be used here.
+- `ramp_down_duration::T`: duration (in s) of the post accumulation ramp down process.
+
+Notes: If a `..._distribution` is specified, then the other parameters for that distribution are no longer used.
+
+# Examples
+```julia-repl
+julia> KellyModel();
+KellyModel{Float64}(6.0, 0.2, 0.4, 0.5, 1.0, 0.1, 0.4, 0.1, 0.2, 0.1)
+```
+
+See also [`LinearModelComponent`](@ref), [`MultichannelComponent`](@ref).
+"""
+Base.@kwdef mutable struct KellyModel <: SequentialSamplingModels.SSM2D
+ drift_rate::Union{Real,String} = 6.0 # drift rate
+ sensor_encoding_delay::Union{Real,String} = 0.2 # mean onset (sensory evidence)
+ sensor_encoding_delay_variability::Union{Real,String} = 0.1 # sensory encoding delay
+ event_onset_distribution::Distribution{Univariate,Continuous} =
+ Normal(sensor_encoding_delay, sensor_encoding_delay_variability) # Normal distribution for sensor encoding delay
+ accumulative_level_noise::Union{Real,String} = 0.5 # accumulation level noise
+ boundary::Union{Real,String} = 1.0 # boundary height
+ start_point::Union{Real,String} = 0
+ start_point_variability::Union{Real,String} = 0.2
+ start_point_distribution::Distribution{Univariate,Continuous} =
+ truncated(Normal(start_point, start_point_variability), upper = boundary)
+ motor_delay::Union{Real,String} = 0.4 # mean motor onset
+ motor_delay_variability::Union{Real,String} = 0.1 # motor delay
+ motor_onset_distribution::Distribution{Univariate,Continuous} =
+ Normal(motor_delay, motor_delay_variability) # Normal distribution for motor delay
+ urgency::Union{Real,String} = 1.0 # mean slope of urgency signal
+ urgency_variability::Union{Real,String} = 0.17 # σ of urgency rate between trials
+ urgency_distribution::Distribution{Univariate,Continuous} =
+ truncated(Normal(urgency, urgency_variability), lower = 0) # Normal distribution for motor delay
+ post_accumulation_duration::Union{Real,String} = 0.1 # mean post-decision duration
+ post_accumulation_duration_variability::Union{Real,String} = 0.001 # variability post-decision
+ post_accumulation_distribution::Distribution{Univariate,Continuous} = truncated(
+ Normal(post_accumulation_duration, post_accumulation_duration_variability),
+ lower = 0,
+ ) # Normal distribution for post accumulation duration
+ ramp_down_duration::Union{Real,String} = 0.1 # CPPrampdown duration
+end
+
+"""
+ create_kelly_parameters_dict(model::KellyModel)
+
+Convert a `KellyModel` instance into a dictionary containing its parameters.
+
+# Arguments
+- `model::KellyModel`: The Kelly model instance whose parameters will be extracted.
+
+# Returns
+- `Dict{Symbol, Any}`: A dictionary where the keys are the parameter names as symbols, and the values are the corresponding parameter values.
+
+# Examples
+```julia-repl
+julia> model = KellyModel(drift_rate=5.5, boundary=1.2)
+KellyModel(5.5, 0.2, 0.4, 0.5, 1.2, 0.1, 0.4, 0.1, 0.2, 0.1)
+
+julia> create_kelly_parameters_dict(model)
+Dict{Symbol, Any}(:drift_rate => 5.5, :event_onset => 0.2, :sensor_encoding_delay => 0.4,
+ :accumulative_level_noise => 0.5, :boundary => 1.2, :motor_onset => 0.1,
+ :motor_delay => 0.4, :post_accumulation_duration => 0.1,
+ :post_accumulation_duration_variability => 0.2, :ramp_down_duration => 0.1)
+"""
+function create_kelly_parameters_dict(model::KellyModel)
+ return Dict(name => getfield(model, name) for name in fieldnames(typeof(model)))
+end
+
+
+"""
+ KellyModel_simulate_cpp(rng, model::KellyModel, time_vec, Δt)
+
+Generate a single response time, and evidence trace of an evidence accumulation process using the Kelly model.
+
+# Arguments
+- `rng::StableRNG`: Random seed to ensure the same traces are created for reconstruction.
+- `model::KellyModel`: specifies the model and its parameters to simulate the evidence accumulation.
+- `time_vec::StepRangeLen`: range of time steps for which the evidence is accumulated.
+- `Δt::Float64`: size of the time steps
+
+# Returns
+- `Float64`: Simulated response time for this trial.
+- `Vector{Float64}`: evidence values over time. The output dimension is `length(time_vec)`.
+
+# Examples
+```julia-repl
+# use the KellyModel and its default parameters to simulate traces from 0:1/500:1.0
+julia> KellyModel_simulate_cpp(StableRNG(1), KellyModel(), 0:1/500:1.0, 1/500)
+(260.70134768436486, [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 … 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0])
+```
+"""
+function KellyModel_simulate_cpp(rng, model::KellyModel, time_vec, Δt)
+ evidence = zeros(length(time_vec))
+ evidence[time_vec.>=rand(rng, model.event_onset_distribution)] .= 1
+ startAccT = time_vec[findfirst(evidence .== 1)]
+
+ noise = vcat(
+ zeros(sum(time_vec .< startAccT)),
+ randn(rng, sum(time_vec .>= startAccT)) .* model.accumulative_level_noise .*
+ sqrt(Δt),
+ )
+
+ cum_evidence = cumsum(evidence .* model.drift_rate .* Δt .+ noise) # This is the cumulative differential evidence, just as in a 1d DDM.
+
+
+ # terminate the decision process on boundary crossing, record threshold-crossing samplepoint:
+ clamp!(cum_evidence, 0, Inf)
+
+
+ urgency_slope = rand(rng, model.urgency_distribution)
+ urgency_slope > 0 ? "" :
+ @warn "urgency slope random sample was $urgency_slope, that is, smaller than 0, this is theoretically not possible. Maybe you forgot a `truncate(...,lower=0)`?"
+ dti = findfirst(cum_evidence .> model.boundary .- urgency_slope .* time_vec) # finding the sample point of threshold crossing of each, then will pick the earlier as the winner
+ if isnothing(dti) # Check if no crossing was found
+ dti = length(time_vec) # Set to the last time step
+ end
+ # now record RT in sec after adding motor time, with variability
+ rt = time_vec[dti] + rand(rng, model.motor_onset_distribution)
+
+ # now make the CPP peak and go down linearly after a certain amount of post-dec accum time for this trial:
+ post_acc_duration = rand(rng, model.post_accumulation_distribution)
+ # so post_acc_duration is the post accumulation duration time, where the accumulation spikes over the threshold
+
+ # acc_stop_index is the accumulation Stop index which is the index from the time Vector where the accumulation really stops
+ acc_stop_index = dti + (post_acc_duration ÷ Δt) |> Int
+ # Take the absolute value of the accumulations
+ cum_evidence = abs.(cum_evidence)
+ if acc_stop_index < length(time_vec)
+ nT = length(time_vec)
+ tmp =
+ cum_evidence[acc_stop_index] .-
+ (1:(nT-acc_stop_index)) .* cum_evidence[acc_stop_index] .*
+ (Δt ./ model.ramp_down_duration)
+ cum_evidence[(acc_stop_index+1):end] .= max.(Ref(0), tmp)
+ end
+ return rt / Δt, cum_evidence[1:end]
+end
+
+
+"""
+ simulate_drift_component(rng, component::DriftComponent, design::AbstractDesign)
+
+Generate response times and evidence Vectors of an given [`AbstractDesign`](@ref) with a [`DriftComponent`](@ref) which contains the model used for the simulation.
+
+# Arguments
+- `rng::StableRNG`: Random seed to ensure the same traces are created for reconstruction.
+- `component::DriftComponent`: Component to specify the model and its parameters to simulate the evidence accumulation.
+- `design::AbstractDesign`: design of the experiment preferable SequenceDesign.
+
+# Returns
+- `Vector{Float64}`: Simulated response times for the trials.
+- `Matrix{Float64}`: evidence values over time for each trial. The output dimensions are `c.max_length x size(events, 1)`.
+
+# Examples
+```julia-repl
+julia> model_parameter = Dict(:motor_onset => 0.4, :event_onset => 0.2);
+
+julia> c = DriftComponent(500, 500, KellyModel, model_parameter);
+
+julia> design_single = SingleSubjectDesign(conditions = Dict(:drift_rate => [0.5, 0.8], :condition => [1]));
+
+julia> design_seq = SequenceDesign(design_single,"SCR_");
+
+julia> simulate_component(StableRNG(1),c,design_seq)
+Vector{Float64}, 501x6 Matrix{Float64}:
+([96.65745162948949, 273.7368235451535, 271.86040880709123, 128.41057786118193, 342.35208862144276, 237.14773586760617], [0.0 0.0 … 0.0 0.0; 0.0 0.0 … 0.0 0.0; … ; 0.0 0.0 … 0.0 0.0; 0.0 0.0 … 0.0 0.0])
+```
+"""
+function simulate_drift_component(rng, component::DriftComponent, design::AbstractDesign)
+ events = generate_events(deepcopy(rng), design)
+ traces = Matrix{Float64}(undef, component.max_length, size(events, 1))
+ rts = Vector{Float64}(undef, size(events, 1))
+ for (i, evt) in enumerate(eachrow(events))
+ parameters = get_model_parameter(rng, evt, component.model_parameters)
+ model =
+ component.model_type(; (key => parameters[key] for key in keys(parameters))...)
+ rt, evidence = SSM_Simulate(rng, model, component.sfreq, component.max_length)
+
+ rts[i] = rt
+ traces[:, i] = evidence[1:component.max_length]
+ end
+ return rts, traces
+end
+
+"""
+ SSM_Simulate(rng, model::SequentialSamplingModels.SSM2D, sfreq, max_length)
+
+Generate response time and evidence Vector of component.max_length by using a SequentialSamplingModels.SSM2D model (tested with DDM and LBA) for simulation.
+
+# Arguments
+- `rng::StableRNG`: Random seed to ensure the same traces are created for reconstruction.
+- `model::SequentialSamplingModels.SSM2D`: SequentialSamplingModel to simulate the evidence and response time.
+- `sfreq::Real`: sample frequency used to simulate the signal.
+- `max_length::Int`: maximum length of the simulated signal.
+
+# Returns
+- `Float64`: Simulated response time for the trial.
+- `Vector{Float64}`: evidence values over time. The output dimension is `component.max_length`.
+
+# Examples
+```julia-repl
+julia> model = DDM()
+julia> SSM_Simulate(StableRNG(1), model, 500, 500)
+Float64, Vector{Float64}:
+(96.65745162948949, [0.0 0.0 … 0.0 0.0])
+```
+"""
+function SSM_Simulate(rng, model::SequentialSamplingModels.SSM2D, sfreq, max_length)
+ if isa(model, LBA) && !(model.τ ≈ 0.0)
+ @warn(
+ "LBA Model with non-0 non-decision. Given we do not know if non-decision time is encoding or response generation, we put everyhing to response generation. We recommend to use τ=0"
+ )
+ end
+ Δt = 1 / sfreq
+
+
+ time_steps, evidence = SequentialSamplingModels.simulate(rng, model; Δt)
+ if !(evidence isa Vector{Float64})
+ evidence = hcat(evidence...)
+ evidence = collect(vec(evidence))
+ end
+ # Store results for this trial
+ rt = (time_steps[end] + model.τ) / Δt
+ if length(evidence) < max_length
+ final_value = 0
+ append!(evidence, fill(final_value, max_length - length(evidence)))
+ else
+ rt = max_length + (model.τ / Δt)
+ evidence = evidence[1:max_length]
+ end
+ return rt, evidence
+end
+
+"""
+ SSM_Simulate(rng, model::KellyModel, sfreq, max_length)
+
+Generate response time and evidence Vector of max_length by using the Kelly Model for the simulation.
+
+# Arguments
+- `rng::StableRNG`: Random seed to ensure the same traces are created for reconstruction.
+- `model::KellyModel`: SequentialSamplingModel to simulate the evidence and response time.
+- `sfreq::Real`: sample frequency used to simulate the trace.
+- `max_length::Int`: maximum length of the simulated trace.
+
+# Returns
+- `Float64`: Simulated response time for the trial.
+- `Vector{Float64}`: evidence values over time. The output dimension is `c.max_length`.
+
+# Examples
+```julia-repl
+julia> model = KellyModel()
+
+julia> SSM_Simulate(StableRNG(1), model, 500, 500)
+Float64, Vector{Float64}:
+(96.65745162948949, [0.0 0.0 … 0.0 0.0])
+```
+"""
+function SSM_Simulate(rng, model::KellyModel, sfreq, max_length)
+ Δt = 1 / sfreq
+ time_vec = 0:Δt:max_length*Δt
+ rt, evidence = KellyModel_simulate_cpp(rng, model, time_vec, Δt)
+ if length(evidence) < max_length
+ final_value = 0
+ append!(evidence, fill(final_value, max_length - length(evidence)))
+ end
+ evidence = evidence[1:max_length]
+ return rt, evidence
+end
\ No newline at end of file
diff --git a/src/types.jl b/src/types.jl
index f0b80467..b94595b2 100644
--- a/src/types.jl
+++ b/src/types.jl
@@ -1,13 +1,20 @@
+"""
+"""
abstract type AbstractOnset end
+"""
+"""
abstract type AbstractNoise end
-
+"""
+"""
abstract type AbstractDesign end
-
+"""
+"""
abstract type AbstractComponent end
# find other types in onset.jl and noise.jl
# and in design.jl and component.jl
-
+"""
+"""
abstract type AbstractHeadmodel end
"""
diff --git a/test/component.jl b/test/component.jl
index 19944924..a3352a8d 100644
--- a/test/component.jl
+++ b/test/component.jl
@@ -44,8 +44,63 @@
@test UnfoldSim.weight_σs(Dict(:subj => [1, 2, [1 0.5; 0.5 1]]), 1.0, 2.0).subj ==
create_re(1, 2; corrmat = [1 0.5; 0.5 1]) ./ 2
end
+
+ @testset "DriftComponent" begin
+ # Test UnfoldSim.simulate_component(rng, c::DriftComponent, design::AbstractDesign)
+ boundary = 1.0
+ model_parameter = Dict(:boundary => boundary)
+ c = DriftComponent(500, 500, KellyModel, model_parameter)
+ design_single = UnfoldSim.SingleSubjectDesign(conditions = Dict(:condition => [1]))
+ design_seq = UnfoldSim.SequenceDesign(design_single, "SCR_")
+ result_traces = UnfoldSim.simulate_component(StableRNG(1), c, design_seq)
+
+ @test size(result_traces) == (500, 3)
+ @test any(result_traces .== 0)
+ @test any(result_traces .>= boundary)
+
+ # Test UnfoldSim.simulate_component(rng, c::DriftComponent, design::AbstractDesign)
+ boundary = 1.0
+ model_parameter = Dict(:boundary => boundary)
+ c = DriftComponent(500, 500, KellyModel, model_parameter)
+ design_single = UnfoldSim.SingleSubjectDesign(
+ conditions = Dict(:drift_rate => [0.5, 0.8], :condition => [1]),
+ )
+ design_seq = UnfoldSim.SequenceDesign(design_single, "SCR_")
+ result_traces = UnfoldSim.simulate_component(StableRNG(1), c, design_seq)
+
+ @test size(result_traces) == (500, 6)
+ @test any(result_traces .== 0)
+ @test any(result_traces .>= boundary)
+
+ # Test calculate_response_times_for_ssm(rng, component::DriftComponent, design::AbstractDesign)
+ model_parameter = Dict()
+ c = DriftComponent(500, 500, KellyModel, model_parameter)
+ design_single = UnfoldSim.SingleSubjectDesign(
+ conditions = Dict(:drift_rate => [0.5, 0.8], :condition => [1]),
+ )
+ design_seq = UnfoldSim.SequenceDesign(design_single, "SCR_")
+ sub_design = UnfoldSim.SubselectDesign(design_seq, 'C')
+ result_rts = UnfoldSim.calculate_response_times_for_ssm(StableRNG(1), c, sub_design)
+ @test size(result_rts) == (2,)
+ @test isapprox(result_rts, [393.3925219909412, 361.18737039119077], atol = 1e-8)
+
+ # Test get_model_parameter(rng, evt, d::Dict)
+ rng = StableRNG(1)
+ model_parameter = Dict(:drift_rate => "drift_rate")
+ c = DriftComponent(500, 500, KellyModel, model_parameter)
+ drift_rates = [0.5, 0.8]
+ design_single = UnfoldSim.SingleSubjectDesign(
+ conditions = Dict(:drift_rate => drift_rates, :condition => [1]),
+ )
+ events = UnfoldSim.generate_events(rng, design_single)
+ for (i, evt) in enumerate(eachrow(events))
+ parameters = UnfoldSim.get_model_parameter(rng, evt, c.model_parameters)
+ @test parameters[:drift_rate] == drift_rates[i]
+ end
+ end
@testset "get_basis" begin
+
rng = StableRNG(1)
design = UnfoldSim.SingleSubjectDesign(; conditions = Dict(:duration => 10:-1:5))
mybasisfun =
@@ -104,5 +159,7 @@
@test length(d) > 10_000
@test e.latency[1] > 10_000
+
+
end
end
diff --git a/test/onset.jl b/test/onset.jl
index d2620c5a..68050fa8 100644
--- a/test/onset.jl
+++ b/test/onset.jl
@@ -103,4 +103,71 @@
end
+
+ @testset "Sequence/Drift Onset" begin
+ rng = StableRNG(1)
+ fs = 500
+ p3 = LinearModelComponent(;
+ basis = UnfoldSim.hanning(Int(0.5 * fs)),
+ formula = @formula(0 ~ 1 + condition),
+ β = [1.0, 0],
+ )
+
+ resp = LinearModelComponent(;
+ basis = UnfoldSim.hanning(Int(0.5 * fs)),
+ formula = @formula(0 ~ 1 + condition),
+ β = [0.5, 0],
+ offset = -10,
+ )
+ sequence_onset = SequenceOnset(
+ Dict(
+ 'S' => UniformOnset(width = 0, offset = 80),
+ 'C' => DriftOnset(),
+ 'R' => UniformOnset(width = 0, offset = 120),
+ ),
+ )
+ model_parameter = Dict(:drift_rate => "drift_rate")
+ drift = UnfoldSim.DriftComponent(500, 500, KellyModel, model_parameter)
+ components = Dict('S' => [p3], 'C' => [drift], 'R' => [resp])
+ design_single = UnfoldSim.SingleSubjectDesign(
+ conditions = Dict(:drift_rate => [0.5, 0.8], :condition => [1]),
+ )
+ design_seq = UnfoldSim.SequenceDesign(design_single, "SCR_")
+ design_rep = UnfoldSim.RepeatDesign(design_seq, 10)
+ simulation = UnfoldSim.Simulation(
+ design_rep,
+ components,
+ sequence_onset,
+ UnfoldSim.NoNoise(),
+ )
+
+ result_onsets = simulate_onsets(rng, sequence_onset, simulation)
+
+ size(result_onsets) == (60,)
+ result_onsets[1] == 121
+ result_onsets[2] == 201
+ result_onsets[3] == 906
+
+ # Test DriftOnset combined with UniformOnset
+ sequence_onset = SequenceOnset(
+ Dict(
+ 'S' => UniformOnset(width = 0, offset = 80),
+ 'C' => (DriftOnset(), UniformOnset(width = 0, offset = 140)),
+ 'R' => UniformOnset(width = 0, offset = 120),
+ ),
+ )
+ simulation = UnfoldSim.Simulation(
+ design_rep,
+ components,
+ sequence_onset,
+ UnfoldSim.NoNoise(),
+ )
+
+ result_onsets = simulate_onsets(rng, sequence_onset, simulation)
+ size(result_onsets) == (60,)
+ result_onsets[1] == 121
+ result_onsets[2] == 201
+ result_onsets[3] == 1046
+ end
+
end
diff --git a/test/sequentialSamplingModelSimulation.jl b/test/sequentialSamplingModelSimulation.jl
new file mode 100644
index 00000000..4337e251
--- /dev/null
+++ b/test/sequentialSamplingModelSimulation.jl
@@ -0,0 +1,76 @@
+@testset "sequentialSamplingModelSimulation" begin
+ fs = 500
+ Δt = 1 / fs # time step
+ tEnd = 1.0 # trial Duration
+ time_vec = 0:Δt:tEnd # time base
+ max_length = tEnd / Δt
+ rng = StableRNG(1)
+ @testset "KellyModel" begin
+ assert_event_onset = 0.663
+ assert_drift_rate = "drift_rate"
+ km = KellyModel(event_onset = assert_event_onset, drift_rate = assert_drift_rate)
+
+ @test km.event_onset == assert_event_onset
+ @test km.drift_rate == assert_drift_rate
+ end
+
+ @testset "KellyModel_simulate_cpp" begin
+ boundary = 1.0
+ result_rt, result_trace = UnfoldSim.KellyModel_simulate_cpp(
+ rng,
+ KellyModel(boundary = boundary),
+ time_vec,
+ Δt,
+ )
+ @test size(result_rt) == ()
+ @test size(result_trace) == (501,)
+ @test isapprox(result_rt, 399.6903067274333, atol = 1e-8)
+ @test any(result_trace .== 0)
+ @test any(result_trace .>= boundary)
+
+ result_sim_rt, result_sim_trace = UnfoldSim.SSM_Simulate(rng, KellyModel(), fs, max_length)
+ @test result_rt == result_sim_rt
+ @test result_trace == result_sim_trace
+ end
+
+ @testset "simulate_drift_component" begin
+ boundary = 1.0
+ model_parameter = Dict(:boundary => boundary);
+ c = UnfoldSim.DriftComponent(
+ max_length,
+ fs,
+ KellyModel,
+ model_parameter,
+ )
+ design_single = UnfoldSim.SingleSubjectDesign(
+ conditions = Dict(:drift_rate => [0.5, 0.8], :condition => [1]),
+ )
+ design_seq = UnfoldSim.SequenceDesign(design_single, "SCR_")
+
+ result_rts, result_traces = UnfoldSim.simulate_drift_component(rng, c, design_seq)
+ @test size(result_rts) == (6,)
+ @test size(result_traces) == (501, 6)
+ @test any(result_traces .>= 1.0)
+ end
+
+ @testset "SSM_Simulate" begin
+ result_rt, result_trace = UnfoldSim.SSM_Simulate(deepcopy(rng), KellyModel(), fs, max_length)
+ @test size(result_rt) == ()
+ @test size(result_trace) == (501,)
+ @test isapprox(result_rt, 399.6903067274333, atol = 1e-8)
+ @test any(result_trace .== 0)
+ @test any(result_trace .>= boundary)
+
+ result_rt, result_trace = UnfoldSim.SSM_Simulate(deepcopy(rng), DDM(), fs, max_length)
+ @test size(result_rt) == ()
+ @test size(result_trace) == (501,)
+ @test isapprox(result_rt, 223.00000000000003, atol = 1e-8)
+ @test any(result_trace .== 0)
+
+ result_rt, result_trace = UnfoldSim.SSM_Simulate(deepcopy(rng), LBA(), fs, max_length)
+ @test size(result_rt) == ()
+ @test size(result_trace) == (501,)
+ @test isapprox(result_rt, 397.0, atol = 1e-8)
+ @test any(result_trace .== 0)
+ end
+end