diff --git a/.all-contributorsrc b/.all-contributorsrc new file mode 100644 index 0000000..853e1aa --- /dev/null +++ b/.all-contributorsrc @@ -0,0 +1,91 @@ +{ + "projectName": "UnfoldSim.jl", + "projectOwner": "unfoldtoolbox", + "files": [ + "README.md" + ], + "commitType": "docs", + "commitConvention": "angular", + "contributorsPerLine": 7, + "contributors": [ + { + "login": "maanikmarathe", + "name": "Maanik Marathe", + "avatar_url": "https://avatars.githubusercontent.com/u/66105649?v=4", + "profile": "https://github.com/maanikmarathe", + "contributions": [ + "doc", + "code" + ] + }, + { + "login": "behinger", + "name": "Benedikt Ehinger", + "avatar_url": "https://avatars.githubusercontent.com/u/10183650?v=4", + "profile": "http://www.benediktehinger.de", + "contributions": [ + "bug", + "code", + "doc", + "ideas", + "infra", + "maintenance", + "review", + "test", + "tutorial" + ] + }, + { + "login": "llips", + "name": "Luis", + "avatar_url": "https://avatars.githubusercontent.com/u/38983684?v=4", + "profile": "https://github.com/llips", + "contributions": [ + "bug", + "code", + "doc", + "ideas" + ] + }, + { + "login": "jschepers", + "name": "Judith Schepers", + "avatar_url": "https://avatars.githubusercontent.com/u/22366977?v=4", + "profile": "https://github.com/jschepers", + "contributions": [ + "ideas", + "bug", + "doc", + "tutorial", + "code" + ] + }, + { + "login": "vladdez", + "name": "Vladimir Mikheev", + "avatar_url": "https://avatars.githubusercontent.com/u/33777074?v=4", + "profile": "https://github.com/vladdez", + "contributions": [ + "bug" + ] + }, + { + "login": "ReboreExplore", + "name": "Manpa Barman", + "avatar_url": "https://avatars.githubusercontent.com/u/43548330?v=4", + "profile": "https://reboreexplore.github.io/", + "contributions": [ + "infra" + ] + }, + { + "login": "ReneSkukies", + "name": "René Skukies", + "avatar_url": "https://avatars.githubusercontent.com/u/57703446?v=4", + "profile": "https://reneskukies.de/", + "contributions": [ + "doc" + ] + } + ] +} diff --git a/.github/workflows/darus-dataverse.yml b/.github/workflows/darus-dataverse.yml new file mode 100644 index 0000000..0ce6ec8 --- /dev/null +++ b/.github/workflows/darus-dataverse.yml @@ -0,0 +1,19 @@ +name: DARUS dataverse sync + +on: + pull_request: + types: [labeled, opened, synchronize, reopened] + + +jobs: + dv-sync: + runs-on: ubuntu-latest + if: contains(github.event.pull_request.labels.*.name, 'sync darus') + steps: + - uses: actions/checkout@v2 + - name: Synchronize to DV + uses: SimTech-Research-Data-Management/dataverse-sync@main + with: + dataverse_url: "https://darus.uni-stuttgart.de/" + api_token: ${{ secrets.DV_API_TOKEN }} + persistent_id: "doi:10.18419/darus-3804" diff --git a/.github/workflows/format.yml b/.github/workflows/format.yml new file mode 100644 index 0000000..33acf68 --- /dev/null +++ b/.github/workflows/format.yml @@ -0,0 +1,8 @@ +name: Format suggestions +on: + pull_request: +jobs: + code-style: + runs-on: ubuntu-latest + steps: + - uses: julia-actions/julia-format@v2.0.0 diff --git a/.github/workflows/joss.yml b/.github/workflows/joss.yml new file mode 100644 index 0000000..c08542a --- /dev/null +++ b/.github/workflows/joss.yml @@ -0,0 +1,30 @@ +name: JOSS generation +on: + push: + paths: + - joss_paper/paper.md + - joss_paper/paper.bib + - .github/workflows/joss.yml + +jobs: + paper: + runs-on: ubuntu-latest + name: Paper Draft + steps: + - name: Checkout + uses: actions/checkout@v4 + - name: Build draft PDF + uses: openjournals/openjournals-draft-action@master + with: + journal: joss + # This should be the path to the paper within your repo. + paper-path: joss_paper/paper.md + - name: Upload + uses: actions/upload-artifact@v3 + with: + name: paper + # This is the output path where Pandoc will write the compiled + # PDF. Note, this should be the same directory as the input + # paper.md + path: joss_paper/ paper.pdf + diff --git a/.gitignore b/.gitignore index 35ee841..2222544 100644 --- a/.gitignore +++ b/.gitignore @@ -6,3 +6,4 @@ /docs/Manifest.toml /docs/build/ /docs/src/literate/*/*.md +/docs/src/generated diff --git a/Artifacts.toml b/Artifacts.toml new file mode 100644 index 0000000..b4390bb --- /dev/null +++ b/Artifacts.toml @@ -0,0 +1,6 @@ +[hartmut] +git-tree-sha1 = "07d7f285b4fed914b2ab5dc25dae3e108b8ce35e" + + [[hartmut.download]] + sha256 = "8fa88ed124e2307ad757956fcdf4f670193a6c58974f3adb8473b56afeb995a7" + url = "https://gist.github.com/behinger/96973aae77e56a73ddeea83792bc52ad/raw/07d7f285b4fed914b2ab5dc25dae3e108b8ce35e.tar.gz" diff --git a/CITATION.bib b/CITATION.bib deleted file mode 100644 index f107738..0000000 --- a/CITATION.bib +++ /dev/null @@ -1,8 +0,0 @@ -@misc{UnfoldSim.jl, - author = {Benedikt Ehinger, Luis Lips, Judith Schepers}, - title = {UnfoldSim.jl}, - url = {https://github.com/unfoldtoolbox/UnfoldSim.jl}, - version = {v0.1.0}, - year = {2023}, - month = {1} -} diff --git a/CITATION.cff b/CITATION.cff index a4b84ba..9a0df5c 100644 --- a/CITATION.cff +++ b/CITATION.cff @@ -2,21 +2,20 @@ # Visit https://bit.ly/cffinit to generate yours today! cff-version: 1.2.0 -title: UnfoldSim +title: UnfoldSim.jl message: 'If you use this page, please cite it as below.' type: software authors: + - given-names: Judith + family-names: Schepers + - 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' - - given-names: Maanik - family-names: Marathe - - given-names: Luis - family-names: Lips - identifiers: - type: doi - value: 10.5281/zenodo.8375512 + value: 10.5281/zenodo.7738651 repository-code: 'https://github.com/unfoldtoolbox/UnfoldSim.jl' -version: 0.1.6 -date-released: '2023-09-25' diff --git a/Project.toml b/Project.toml index f68b3cc..922df7e 100644 --- a/Project.toml +++ b/Project.toml @@ -1,12 +1,15 @@ name = "UnfoldSim" uuid = "ed8ae6d2-84d3-44c6-ab46-0baf21700804" -authors = ["Benedikt Ehinger","Luis Lips", "Judith Schepers","Maanik Marathe"] -version = "0.1.6" +authors = ["Benedikt Ehinger", "Luis Lips", "Judith Schepers", "Maanik Marathe"] +version = "0.1.7" [deps] +Artifacts = "56f22d72-fd6d-98f1-02f0-08ddc0907c33" DSP = "717857b8-e6f2-59f4-9121-6e50c889abd2" DataFrames = "a93c6f00-e57d-5684-b7b6-d8193f3e46c0" Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f" +FileIO = "5789e2e9-d7fb-5bc7-8068-2c6fae9b9549" +HDF5 = "f67ccb44-e63f-5c2f-98bd-6dc0ccc4ba2f" ImageFiltering = "6a3955dd-da59-5b1f-98d4-e7296123deb5" LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" MixedModels = "ff71e718-51f3-5ec2-a782-8ffcbfa3c316" @@ -14,18 +17,25 @@ MixedModelsSim = "d5ae56c5-23ca-4a1f-b505-9fc4796fc1fe" Parameters = "d96e819e-fc66-5662-9728-84c9c7592b0a" Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c" SignalAnalysis = "df1fea92-c066-49dd-8b36-eace3378ea47" +Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2" StatsModels = "3eaba693-59b7-5ba5-a881-562e759f1c8d" ToeplitzMatrices = "c751599d-da0a-543b-9d20-d0a503d91d24" [compat] +Artifacts = "1" DSP = "0.7" DataFrames = "1" Distributions = "0.25" +FileIO = "1" +HDF5 = "0.17" ImageFiltering = "0.7" +LinearAlgebra = "1" MixedModels = "4" MixedModelsSim = "0.2" Parameters = "0.12" +Random = "1" SignalAnalysis = "0.4, 0.5,0.6" +Statistics = "1" StatsModels = "0.6,0.7" ToeplitzMatrices = "0.7, 0.8" julia = "1" diff --git a/README.md b/README.md index 283cb11..d80d000 100644 --- a/README.md +++ b/README.md @@ -10,10 +10,10 @@ |---|---|---|---|---|---| | | ||||| -A package to simulate single timeseries model-based ERPs, fMRI activity, pupil dilation etc. -If you have one channel, it is a timeseries of (overlapping) event-related activity and some noise - you might have fun here! +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 available in the documentation! +Many Tutorials, Guides, How-Tos and References available in the [documentation](https://unfoldtoolbox.github.io/UnfoldSim.jl/dev/)! ![grafik](https://user-images.githubusercontent.com/10183650/213565922-90feec23-3b51-4602-b50c-31561dbfc261.png) @@ -49,9 +49,9 @@ Once the toolbox is registered this will translate to ```Pkg.add("UnfoldSim")``` # Quickstart ```julia -data,evts = UnfoldSim.predef_eeg(;n_trials=20,noiselevel=0.8) +data,evts = UnfoldSim.predef_eeg(;n_repeats=1,noiselevel=0.8) ``` -Produces continuous "EEG" with PinkNoise and some Overlap between 20 events +Produces continuous "EEG" with PinkNoise and some Overlap between 20 events (2 conditions * 10 levels of continuous variable). ## Slightly longer ```julia @@ -82,12 +82,34 @@ You are very welcome to raise issues and start pull requests! 2. Literate.jl converts the `.jl` file to a `.md` automatically and places it in `docs/src/generated/FOLDER/FILENAME.md`. 3. Edit [make.jl](https://github.com/unfoldtoolbox/Unfold.jl/blob/main/docs/make.jl) with a reference to `docs/src/generated/FOLDER/FILENAME.md`. -## Contributors (alphabetically) +## Contributors + + + + + + + + + + + + + + + +
Maanik Marathe
Maanik Marathe
📖 💻		Benedikt Ehinger
Benedikt Ehinger
🐛 💻 📖 🤔 🚇 🚧 👀 ⚠️ ✅		Luis
Luis
🐛 💻 📖 🤔		Judith Schepers
Judith Schepers
🤔 🐛 📖 ✅ 💻		Vladimir Mikheev
Vladimir Mikheev
🐛		Manpa Barman
Manpa Barman
🚇		René Skukies
René Skukies
📖
-- **Benedikt Ehinger** -- **Luis Lips** -- **Maanik Marathe** -- **Judith Schepers** + + + + + + + +This project follows the [all-contributors](https://allcontributors.org/docs/en/specification) specification. + +Contributions of any kind welcome! ## Citation diff --git a/docs/Project.toml b/docs/Project.toml index c859cb9..9ee1484 100644 --- a/docs/Project.toml +++ b/docs/Project.toml @@ -5,7 +5,9 @@ DataFrames = "a93c6f00-e57d-5684-b7b6-d8193f3e46c0" Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4" Glob = "c27321d9-0574-5035-807b-f59d2c89b15c" HypothesisTests = "09f84164-cd44-5f33-b23f-e6b0d136a0d5" +InteractiveUtils = "b77e0a4c-d291-57a0-90e8-8db25a27a240" Literate = "98b081ad-f1c9-55d3-8b20-4c87d4299306" +LiveServer = "16fef848-5104-11e9-1b77-fb7a48bbb589" Parameters = "d96e819e-fc66-5662-9728-84c9c7592b0a" StableRNGs = "860ef19b-820b-49d6-a774-d7a799459cd3" StatsBase = "2913bbd2-ae8a-5f71-8c99-4fb6c76f3a91" diff --git a/docs/literate/HowTo/multichannel.jl b/docs/literate/HowTo/multichannel.jl new file mode 100644 index 0000000..15a0361 --- /dev/null +++ b/docs/literate/HowTo/multichannel.jl @@ -0,0 +1,58 @@ + +using UnfoldSim +using UnfoldMakie +using CairoMakie +using DataFrames +using Random + + +# ## Specifying a design + +# We are using a one-level design for testing here. +design = SingleSubjectDesign(conditions=Dict(:condA=>["levelA"])) + +# Next we generate two simple components at two different times without any formula attached (we have a single condition anyway) +c = LinearModelComponent(;basis=p100(),formula = @formula(0~1),β = [1]); +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.: +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) +hart = headmodel(type="hartmut") +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 +# +# finally we need to define the onsets of the signal +onset = UniformOnset(;width=20,offset=4); + +# ## Simulation + +# Now as usual we simulate data. Inspecting data shows our result is now indeed ~230 Electrodes large! Nice! +data,events = simulate(MersenneTwister(1),design, [mc,mc2], onset, PinkNoise(noiselevel=0.05)) +size(data) + + +# !!! hint +# The noise declared in the `simulate` function is added after mixing to channels, each channel receives independent noise. It is also possible to add noise to each individual component+source prior to projection. This would introduce correlated noise. +# +# ## Plotting +# Let's plot using Butterfly & Topoplot +# first we convert the electrodes to positions usable in TopoPlots.jl +pos3d = hart.electrodes["pos"]; +pos2d = to_positions(pos3d') +pos2d = [Point2f(p[1]+0.5,p[2]+0.5) for p in pos2d]; + +# now plot! +f = Figure() +df = DataFrame(:estimate => data[:],:channel => repeat(1:size(data,1),outer=size(data,2)),:time => repeat(1:size(data,2),inner=size(data,1))) +plot_butterfly!(f[1,1:2],df;positions=pos2d) +plot_topoplot!(f[2,1],df[df.time .== 28,:];positions=pos2d,visual=(;enlarge=0.5,label_scatter=false),axis=(;limits=((0,1),(0,0.9)))) +plot_topoplot!(f[2,2],df[df.time .== 48,:];positions=pos2d,visual=(;enlarge=0.5,label_scatter=false),axis=(;limits=((0,1),(0,0.9)))) +f + diff --git a/docs/src/literate/HowTo/newComponent.jl b/docs/literate/HowTo/newComponent.jl similarity index 78% rename from docs/src/literate/HowTo/newComponent.jl rename to docs/literate/HowTo/newComponent.jl index 65d6955..23332fe 100644 --- a/docs/src/literate/HowTo/newComponent.jl +++ b/docs/literate/HowTo/newComponent.jl @@ -1,18 +1,16 @@ # # New component: Duration + Shift -# We want a new component that changes it's duration and shift depending on a column in the event-design. This is somewhat already implemented in the HRF + Pupil bases -begin - using UnfoldSim - using Unfold - using Random - using DSP - using CairoMakie -end +# 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 +using UnfoldSim +using Unfold +using Random +using DSP +using CairoMakie, UnfoldMakie -sfreq = 100 +sfreq = 100; -## Design -# Let's genrate a design with two columns, shift + duration +# ## Design +# Let's generate a design with two columns, shift + duration design = UnfoldSim.SingleSubjectDesign(;conditions= Dict( :shift => rand(100).*sfreq/5, :duration=>20 .+rand(100).*sfreq/5)) @@ -24,7 +22,7 @@ struct TimeVaryingComponent <: AbstractComponent maxlength end -# we have to define the length of a component +# We have to define the length of a component Base.length(c::TimeVaryingComponent) = length(c.maxlength) # While we could have put the TimeVaryingComponent.basisfunction directly into the simulate function, I thought this is a bit more modular @@ -52,4 +50,4 @@ end erp = UnfoldSim.simulate(MersenneTwister(1),TimeVaryingComponent(basis_shiftduration,50),design) -heatmap(hcat(erp...)) \ No newline at end of file +plot_erpimage(hcat(erp...),sortvalues=generate(design).shift) \ No newline at end of file diff --git a/docs/src/literate/HowTo/newDesign.jl b/docs/literate/HowTo/newDesign.jl similarity index 90% rename from docs/src/literate/HowTo/newDesign.jl rename to docs/literate/HowTo/newDesign.jl index 0b8ea37..9ea5761 100644 --- a/docs/src/literate/HowTo/newDesign.jl +++ b/docs/literate/HowTo/newDesign.jl @@ -4,7 +4,7 @@ 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` +# 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 # @@ -20,7 +20,7 @@ end; # we need a `size(design::ImbalanceSubjectDesign)` function to tell how many events we will have. # This is used at different places, e.g. in the Default onset implementation -## note the trailling , to make it a Tuple +## note the trailing , to make it a Tuple size(design::ImbalanceSubjectDesign) = (design.nTrials,); # #### 3) `generate` @@ -37,6 +37,6 @@ end; design = ImbalanceSubjectDesign(;nTrials=6,balance=0.2) generate(design) -# !!! important -# it is the users task to ensure that each run is reproducible. So if you have a random process (e.g. shuffling), be sure to +# !!! warning "Important" +# It is the users task to ensure that each run is reproducible. So if you have a random process (e.g. shuffling), be sure to # safe a RNG object in your struct and use it in your generate function. \ No newline at end of file diff --git a/docs/src/literate/HowTo/repeatTrials.jl b/docs/literate/HowTo/repeatTrials.jl similarity index 59% rename from docs/src/literate/HowTo/repeatTrials.jl rename to docs/literate/HowTo/repeatTrials.jl index 8e984cd..e544e16 100644 --- a/docs/src/literate/HowTo/repeatTrials.jl +++ b/docs/literate/HowTo/repeatTrials.jl @@ -2,24 +2,24 @@ using UnfoldSim # ## Repeating Design entries -# Sometimes we want to repeat a design, that is, have multiple trials with identical values, but it is not always straight forward to implement +# 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 +# If you struggle with this problem `RepeatDesign` is an easy tool for you: designOnce = MultiSubjectDesign(; - n_items=2, - n_subjects = 2, - subjects_between =Dict(:cond=>["levelA","levelB"]), - items_between =Dict(:cond=>["levelA","levelB"]), - ); + n_items = 2, + n_subjects = 2, + subjects_between = Dict(:cond => ["levelA", "levelB"]), + items_between = Dict(:cond => ["levelA", "levelB"]), +); -design = RepeatDesign(designOnce,4); +design = RepeatDesign(designOnce, 4); generate(design) -# As you can see, the design was simply repeated. As always, you can ignore the `dv` column, it is for internal consistency with MixedModelsSim.jl +# As you can see, the design was simply repeated. # !!! note -# if you implemented your own `AbstractDesign`, you need to define the size function accordingly. E.g.: +# If you implemented your own `AbstractDesign`, you need to define the size function accordingly. E.g.: # `Base.size(design::RepeatDesign{SingleSubjectDesign}) = size(design.design).*design.repeat` diff --git a/docs/src/literate/reference/basistypes.jl b/docs/literate/reference/basistypes.jl similarity index 85% rename from docs/src/literate/reference/basistypes.jl rename to docs/literate/reference/basistypes.jl index c14b83e..290abd2 100644 --- a/docs/src/literate/reference/basistypes.jl +++ b/docs/literate/reference/basistypes.jl @@ -2,15 +2,15 @@ using UnfoldSim using CairoMakie using DSP using StableRNGs -## + # ## Basistypes -# There are several bases types directly implemented. They can be easily used for the `components`. +# 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 +# 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] @@ -46,8 +46,8 @@ for (ix,pl) = enumerate(plotConfig) end f -## pupil -# we use the simplified PuRF from Hoeks & Levelt, 1993. Note that https://www.science.org/doi/10.1126/sciadv.abi9979 show some evidence in their supplementary material, that the convolution model is not fully applicable. +# ## Pupil +# We use the simplified PuRF from Hoeks & Levelt, 1993. Note that https://www.science.org/doi/10.1126/sciadv.abi9979 show some evidence in their supplementary material, that the convolution model is not fully applicable. f = Figure() plotConfig = (:n=>5:3:15, :tmax=>0.5:0.2:1.1, diff --git a/docs/literate/reference/noisetypes.jl b/docs/literate/reference/noisetypes.jl new file mode 100644 index 0000000..a79a1e5 --- /dev/null +++ b/docs/literate/reference/noisetypes.jl @@ -0,0 +1,37 @@ +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 = f[1, 1:2] = Axis(f; title = "1.000 samples of noise") +ax_spec = f[2, 1] = Axis(f; title = "Welch Periodigram") +ax_auto = f[2, 2] = Axis(f; title = "Autocorrelogram (every 10th lag)") +for n in [PinkNoise RedNoise WhiteNoise NoNoise ExponentialNoise] + + ## generate + noisevec = gen_noise(StableRNG(1), n(), 10000) + + ## plot 1000 samples + lines!(ax_sig, noisevec[1:1000]; label = string(n)) + + ## calc spectrum + perio = welch_pgram(noisevec) + + ## plot spectrum + lines!(ax_spec, freq(perio), log10.(power(perio))) + + lags = 0:10:500 + autocor_vec = autocor(noisevec, lags) + lines!(ax_auto, lags, autocor_vec) + +end +f[1:2, 3] = Legend(f, ax_sig, "NoiseType") +f + + +# !!! hint +# We recommed for smaller signals the `ExponentialNoise`, maybe with a removed DC offset or a HighPass filter. For long signals, this Noise requires lots of memory though. maybe Pinknoise is a better choice then. \ No newline at end of file diff --git a/docs/literate/reference/onsettypes.jl b/docs/literate/reference/onsettypes.jl new file mode 100644 index 0000000..0634e05 --- /dev/null +++ b/docs/literate/reference/onsettypes.jl @@ -0,0 +1,281 @@ +# # 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`. + +# ## Setup +# ```@raw html +#
+# Click to expand +# ``` + +using UnfoldSim +using CairoMakie +using Random + +## Define a simple design and repeat it 10000. +## This will result in 20000 events i.e. event onsets. +design = + SingleSubjectDesign(conditions = Dict(:cond => ["A", "B"])) |> + x -> RepeatDesign(x, 10000); + +# ```@raw html +#
+# ``` + +# ## UniformOnset +# The `UniformOnset` is based on a uniform distribution and has two parameters: `width` and `offset`. + +# Example: +onset_uniform = UniformOnset(; width = 50, offset = 0); + +# The `width` parameter defines the upper bound of the interval of the uniform distribution (its lower bound is 0) i.e. all values between 0 and `width` are equally probable. + +# The `offset` parameter determines the minimal distance between two events and its value is added to the value sampled from the uniform distribution i.e. it shifts the distribution. +# Its default value is `0`, i.e. no offset. + +# In the figure below, it is illustrated how the onset distribution changes when changing one of its parameters. +let # hide + f = Figure() # hide + + ## Define parameter combinations # hide + parameters = [ # hide + (((50, 0), (80, 0)), "width"), # hide + (((50, 0), (50, 20)), "offset"), # hide + ] # hide + + axes_list = Array{Any}(undef, length(parameters)) # hide + + ## Create a subplot for each parameter i.e. one for width and one for offset # hide + for (index, (combinations, label)) in enumerate(parameters) # hide + ax = Axis(f[index, 1], title = "Parameter: $label") # hide + axes_list[index] = ax # hide + + ## Go through all parameter combinations and plot a histogram of the sampled onsets # hide + for (width, offset) in combinations # hide + onsets = UnfoldSim.rand_onsets( # hide + MersenneTwister(42), # hide + UniformOnset(; width = width, offset = offset), # hide + design, # hide + ) # hide + + hist!(ax, onsets, bins = range(0, 100, step = 1), label = "($width, $offset)") # hide + + if label == "offset" && offset != 0 # hide + vlines!(offset, color = "black") # hide + end # hide + end # hide + hideydecorations!(ax) # hide + hidespines!(ax, :t, :r) # hide + axislegend( # hide + ax, # hide + framevisible = false, # hide + labelsize = 12, # hide + markersize = 5, # hide + patchsize = (10, 10), # hide + ) # hide + end # hide + axes_list[end].xlabel = "Time between events [samples]" # hide + linkyaxes!(axes_list...) # hide + current_figure() # hide +end # hide + +## Note: The code is repeated because I did not manage to show the figure but make the code collapsible # hide +# ```@raw html +#
+# Click to show the code for the figure above +# ``` +let + f = Figure() + + ## Define parameter combinations + parameters = [(((50, 0), (80, 0)), "width"), (((50, 0), (50, 20)), "offset")] + + axes_list = Array{Any}(undef, length(parameters)) + + ## Create a subplot for each parameter i.e. one for width and one for offset + for (index, (combinations, label)) in enumerate(parameters) + ax = Axis(f[index, 1], title = "Parameter: $label") + axes_list[index] = ax + + ## Go through all parameter combinations and plot a histogram of the sampled onsets + for (width, offset) in combinations + onsets = UnfoldSim.rand_onsets( + MersenneTwister(42), + UniformOnset(; width = width, offset = offset), + design, + ) + + hist!(ax, onsets, bins = range(0, 100, step = 1), label = "($width, $offset)") + + if label == "offset" && offset != 0 + vlines!(offset, color = "black") + end + end + hideydecorations!(ax) + hidespines!(ax, :t, :r) + axislegend( + ax, + framevisible = false, + labelsize = 12, + markersize = 5, + patchsize = (10, 10), + ) + end + axes_list[end].xlabel = "Time between events [samples]" + linkyaxes!(axes_list...) +end +# ```@raw html +#
+# ``` + +# ## LogNormalOnset +# The `LogNormalOnset` is based on a log-normal distribution and has four parameters: `μ`, `σ`, `offset` and `truncate_upper`. + +# Example: +onset_lognormal = LogNormalOnset(; μ = 3, σ = 0.25, offset = 0, truncate_upper = nothing); + +# The parameters `μ` and `σ` are the location and scale parameter of the log-normal distribution. However, they are not identical to its mean and standard deviation. +# If a variable $X$ is log-normally distributed then $Y = ln(X)$ is normally distributed with mean `μ` and standard deviation `σ`[^1]. + +# The `offset` parameter determines the minimal distance between two events and its value is added to the value sampled from the log-normal distribution i.e. it shifts the distribution. +# Its default value is `0`, i.e. no offset. + +# The `truncate_upper` parameter allows to truncate the distribution at a certain sample value. Its default value is `nothing`, i.e. no truncation. + +# In the figure below, it is illustrated how the onset distribution changes when changing one of its parameters. +let # hide + f = Figure(size = (600, 800)) # hide + + ## Define parameter combinations # hide + parameters = [ # hide + (((3, 0.25, 0, nothing), (2.5, 0.25, 0, nothing)), "μ"), # hide + (((3, 0.25, 0, nothing), (3, 0.35, 0, nothing)), "σ"), # hide + (((3, 0.25, 0, nothing), (3, 0.25, 30, nothing)), "offset"), # hide + (((3, 0.25, 0, nothing), (3, 0.25, 0, 25)), "truncate_upper"), # hide + ] # hide + + axes_list = Array{Any}(undef, length(parameters)) # hide + + ## Create a subplot for each parameter i.e. one for μ, one for σ etc # hide + for (index, (combinations, label)) in enumerate(parameters) # hide + ax = Axis(f[index, 1], title = "Parameter: $label") # hide + axes_list[index] = ax # hide + + ## Go through all parameter combinations and plot a histogram of the sampled onsets # hide + for (μ, σ, offset, truncate_upper) in combinations # hide + onsets = UnfoldSim.rand_onsets( # hide + MersenneTwister(42), # hide + LogNormalOnset(; # hide + μ = μ, # hide + σ = σ, # hide + offset = offset, # hide + truncate_upper = truncate_upper, # hide + ), # hide + design, # hide + ) # hide + + hist!( # hide + ax, # hide + onsets, # hide + bins = range(0, 100, step = 1), # hide + label = "($μ,$σ,$offset,$truncate_upper)", # hide + ) # hide + + if label == "offset" && offset !== 0 # hide + vlines!(offset, color = "black") # hide + elseif label == "truncate_upper" && truncate_upper !== nothing # hide + vlines!(truncate_upper, color = "black") # hide + end # hide + end # hide + hideydecorations!(ax) # hide + hidespines!(ax, :t, :r) # hide + axislegend( # hide + ax, # hide + framevisible = false, # hide + labelsize = 12, # hide + markersize = 5, # hide + patchsize = (10, 10), # hide + ) # hide + end # hide + axes_list[end].xlabel = "Time between events [samples]" # hide + linkyaxes!(axes_list...) # hide + current_figure() # hide +end # hide + + +## Note: The code is repeated because I did not manage to show the figure but make the code collapsible # hide +# ```@raw html +#
+# Click to show the code for the figure above +# ``` +let + f = Figure(size = (600, 800)) + + ## Define parameter combinations + parameters = [ + (((3, 0.25, 0, nothing), (2.5, 0.25, 0, nothing)), "μ"), + (((3, 0.25, 0, nothing), (3, 0.35, 0, nothing)), "σ"), + (((3, 0.25, 0, nothing), (3, 0.25, 30, nothing)), "offset"), + (((3, 0.25, 0, nothing), (3, 0.25, 0, 25)), "truncate_upper"), + ] + + axes_list = Array{Any}(undef, length(parameters)) + + ## Create a subplot for each parameter i.e. one for μ, one for σ etc + for (index, (combinations, label)) in enumerate(parameters) + ax = Axis(f[index, 1], title = "Parameter: $label") + axes_list[index] = ax + + ## Go through all parameter combinations and plot a histogram of the sampled onsets + for (μ, σ, offset, truncate_upper) in combinations + onsets = UnfoldSim.rand_onsets( + MersenneTwister(42), + LogNormalOnset(; + μ = μ, + σ = σ, + offset = offset, + truncate_upper = truncate_upper, + ), + design, + ) + + hist!( + ax, + onsets, + bins = range(0, 100, step = 1), + label = "($μ,$σ,$offset,$truncate_upper)", + ) + + if label == "offset" && offset !== 0 + vlines!(offset, color = "black") + elseif label == "truncate_upper" && truncate_upper !== nothing + vlines!(truncate_upper, color = "black") + end + end + hideydecorations!(ax) + hidespines!(ax, :t, :r) + axislegend( + ax, + framevisible = false, + labelsize = 12, + markersize = 5, + patchsize = (10, 10), + ) + end + axes_list[end].xlabel = "Time between events [samples]" + linkyaxes!(axes_list...) +end +# ```@raw html +#
+# ``` + +# # Overlap of subsequent events +# !!! note +# The overlap of subsequent events can be indirectly controlled by setting the `offset` parameter relative to the length of the component basis. +# Assuming that `signal` is a component e.g. `LinearModelComponent`, +# - if `offset` > `length(signal.basis)` -> no overlap +# - if `offset` < `length(signal.basis)` -> there might be overlap, depending on the other parameters of the onset distribution + + +## Footnotes # hide +# [^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# \ No newline at end of file diff --git a/docs/literate/reference/overview.jl b/docs/literate/reference/overview.jl new file mode 100644 index 0000000..6114079 --- /dev/null +++ b/docs/literate/reference/overview.jl @@ -0,0 +1,22 @@ +# # Overview of functionality +# UnfoldSim has many modules, here we try to collect them to provide you with an overview. +using UnfoldSim +using InteractiveUtils + +# ## Design +# Designs define the experimental design. They can be nested, e.g. `RepeatDesign(SingleSubjectDesign,10)` would repeat the generated design-dataframe 10x. +subtypes(AbstractDesign) + +# ## Component +# Components define a signal. Some components can be nested, e.g. `LinearModelComponent|>MultichannelComponent`, see the multi-channel tutorial for more information. +subtypes(AbstractComponent) + +# ## Onsets +# Onsets define the distance between events in the continuous signal. +subtypes(AbstractOnset) + +# ## Noise +# Choose the noise you need! +subtypes(AbstractNoise) + + diff --git a/docs/src/literate/tutorials/poweranalysis.jl b/docs/literate/tutorials/poweranalysis.jl similarity index 100% rename from docs/src/literate/tutorials/poweranalysis.jl rename to docs/literate/tutorials/poweranalysis.jl diff --git a/docs/src/literate/tutorials/quickstart.jl b/docs/literate/tutorials/quickstart.jl similarity index 92% rename from docs/src/literate/tutorials/quickstart.jl rename to docs/literate/tutorials/quickstart.jl index 76a9f2d..04a9d62 100644 --- a/docs/src/literate/tutorials/quickstart.jl +++ b/docs/literate/tutorials/quickstart.jl @@ -3,10 +3,10 @@ using Random using CairoMakie -# !!! tipp -# Use `subtypes(AbstractNoise)` (or `subtypes(AbstractComponent)` etc.) to find already implemented building blocks +# !!! tip +# Use `subtypes(AbstractNoise)` (or `subtypes(AbstractComponent)` etc.) to find already implemented building blocks. -# #### "Experimental" Design +# ## "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"]) @@ -15,7 +15,7 @@ design = SingleSubjectDesign(; # #### Component / Signal # Define a simple component and ground truth simulation formula. Akin to ERP components, we call one simulation signal a component. # -# !!! highlight +# !!! note # 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], @@ -30,7 +30,7 @@ onset = UniformOnset(;width=20,offset=4); # And we will use some noise noise = PinkNoise(;noiselevel=0.2); -# #### Combine & Generate +# ## Combine & Generate # We will put it all together in one `Simulation` type simulation = Simulation(design, signal, onset, noise); diff --git a/docs/literate/tutorials/simulateERP.jl b/docs/literate/tutorials/simulateERP.jl new file mode 100644 index 0000000..f71a80b --- /dev/null +++ b/docs/literate/tutorials/simulateERP.jl @@ -0,0 +1,76 @@ +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. + +# Let's grab a SingleSubjectDesign and add a continuous predictor +design = + SingleSubjectDesign(; + conditions = Dict( + :condition => ["car", "face"], + :continuous => range(-5, 5, length = 10), + ), + ) |> x -> RepeatDesign(x, 100); + +# Let's make use of the prespecified basis functions, but use different formulas + parameters for each! + +# **p100** is unaffected by our design and has amplitude of 5 +p1 = LinearModelComponent(; basis = p100(), formula = @formula(0 ~ 1), β = [5]); + +# **n170** has a condition effect, faces are more negative than cars +n1 = LinearModelComponent(; + basis = n170(), + formula = @formula(0 ~ 1 + condition), + β = [5, -3], +); +# **p300** has a continuous effect, higher continuous values will result in larger P300's. +# We include both a linear and a quadratic effect of the continuous variable. +p3 = LinearModelComponent(; + basis = p300(), + formula = @formula(0 ~ 1 + continuous + continuous^2), + β = [5, 1, 0.2], +); + +# Now we can simply combine the components and simulate +components = [p1, n1, p3] +data, evts = simulate( + MersenneTwister(1), + design, + components, + UniformOnset(; width = 0, offset = 1000), + PinkNoise(), +); + + +# ## Analysis +# Let's check that everything worked out well, by using Unfold +m = fit( + UnfoldModel, + Dict( + Any => ( + @formula(0 ~ 1 + condition + spl(continuous, 4)), + firbasis(τ = [-0.1, 1], sfreq = 100, name = "basis"), + ), + ), + evts, + data, +); + +# first the "pure" beta/linear regression parameters +plot_erp(coeftable(m)) + +# and now beautifully visualized as marginal betas / predicted ERPs +f = plot_erp( + effects(Dict(:condition => ["car", "face"], :continuous => -5:5), m); + mapping = (:color => :continuous, linestyle = :condition, group = :continuous), + legend = (; valign = :top, halign = :right, tellwidth = false), + categorical_color = false, +); + +## Workaround to separate legend and colorbar (will be fixed in a future UnfoldMakie version) +legend = f.content[2] +f[:, 1] = legend +current_figure() \ No newline at end of file diff --git a/docs/make.jl b/docs/make.jl index 2f347d4..74d8e12 100644 --- a/docs/make.jl +++ b/docs/make.jl @@ -4,21 +4,24 @@ using Glob using Literate -GENERATED = joinpath(@__DIR__, "src", "literate") +GENERATED = joinpath(@__DIR__, "src", "generated") +SOURCE = joinpath(@__DIR__, "literate") + for subfolder ∈ ["explanations","HowTo","tutorials","reference"] - local SOURCE_FILES = Glob.glob(subfolder*"/*.jl", GENERATED) + local SOURCE_FILES = Glob.glob(subfolder*"/*.jl", SOURCE) #config=Dict(:repo_root_path=>"https://github.com/unfoldtoolbox/UnfoldSim") foreach(fn -> Literate.markdown(fn, GENERATED*"/"*subfolder), SOURCE_FILES) end -DocMeta.setdocmeta!(UnfoldSim, :DocTestSetup, :(using UnfoldSim); recursive=true) +DocMeta.setdocmeta!(UnfoldSim, :DocTestSetup, :(using UnfoldSim); recursive = true) makedocs(; - modules=[UnfoldSim], -authors="Luis Lips, Benedikt Ehinger, Judith Schepers", - repo="https://github.com/unfoldtoolbox/UnfoldSim.jl/blob/{commit}{path}#{line}", + modules = [UnfoldSim], + authors = "Luis Lips, Benedikt Ehinger, Judith Schepers", + #repo="https://github.com/unfoldtoolbox/UnfoldSim.jl/blob/{commit}{path}#{line}", + repo=Documenter.Remotes.GitHub("unfoldtoolbox", "UnfoldSim.jl"), sitename="UnfoldSim.jl", format=Documenter.HTML(; prettyurls=get(ENV, "CI", "false") == "true", @@ -29,18 +32,21 @@ authors="Luis Lips, Benedikt Ehinger, Judith Schepers", pages=[ "Home" => "index.md", "Tutorials"=>[ - "Quickstart" => "literate/tutorials/quickstart.md", - "Simulate ERPs" => "literate/tutorials/simulateERP.md", - "Poweranalysis" => "literate/tutorials/poweranalysis.md", + "Quickstart" => "generated/tutorials/quickstart.md", + "Simulate ERPs" => "generated/tutorials/simulateERP.md", + "Poweranalysis" => "generated/tutorials/poweranalysis.md", ], - "Reference"=>[ - "NoiseTypes" => "./literate/reference/noisetypes.md", - "ComponentBasisTypes" => "./literate/reference/basistypes.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", ], "HowTo" => [ - "New Experimental Design" => "./literate/HowTo/newDesign.md", - "Repeating Trials within a Design" => "./literate/HowTo/repeatTrials.md", - "New Duration/Shift-dependent Component" => "./literate/HowTo/newComponent.md" + "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", + "Generate multi channel data" =>"./generated/HowTo/multichannel.md", ], "DocStrings" => "api.md", ], diff --git a/docs/run_liveserver.jl b/docs/run_liveserver.jl new file mode 100644 index 0000000..efb6d45 --- /dev/null +++ b/docs/run_liveserver.jl @@ -0,0 +1,7 @@ +using LiveServer + +dr = filter(isdir,readdir(joinpath("src","generated"),join=true)) +push!(dr,"./build") +servedocs(skip_dirs=dr,literate_dir=joinpath("literate"),foldername=".",host="0.0.0.0") + + diff --git a/docs/src/literate/reference/noisetypes.jl b/docs/src/literate/reference/noisetypes.jl deleted file mode 100644 index cbbb909..0000000 --- a/docs/src/literate/reference/noisetypes.jl +++ /dev/null @@ -1,37 +0,0 @@ -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 = f[1,1:2] = Axis(f;title="1.000 samples of noise") -ax_spec = f[2,1] = Axis(f;title="Welch Periodigram") -ax_auto = f[2,2] = Axis(f;title="Autocorrelogram (every 10th lag)") -for n = [PinkNoise RedNoise WhiteNoise NoNoise ExponentialNoise] - - ## generate - noisevec = gen_noise(StableRNG(1),n(),10000) - - ## plot 1000 samples - lines!(ax_sig,noisevec[1:1000];label=string(n)) - - ## calc spectrum - perio = welch_pgram(noisevec) - - ## plot spectrum - lines!(ax_spec,freq(perio),log10.(power(perio))) - - lags = 0:10:500 - autocor_vec = autocor(noisevec,lags) - lines!(ax_auto,lags,autocor_vec) - -end -f[1:2,3] = Legend(f,ax_sig,"NoiseType") -f - - -# !!! Recommendation -# We recommed for smaller signals the `ExponentialNoise`, maybe with a removed DC offset or a HighPass filter. For long signals, this Noise requires lot's of memory though. maybe Pinknoise is a better choice \ No newline at end of file diff --git a/docs/src/literate/tutorials/simulateERP.jl b/docs/src/literate/tutorials/simulateERP.jl deleted file mode 100644 index 809b95f..0000000 --- a/docs/src/literate/tutorials/simulateERP.jl +++ /dev/null @@ -1,51 +0,0 @@ -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 - -# let's grab a SingleSubjectDesign and add a continuous predictor -design = SingleSubjectDesign(; - conditions=Dict(:condition=>["car","face"],:continuous=>range(-5,5,length=10)) - ) |> x->RepeatDesign(x,100); - -# let's make use of the prespecified basis functions, but use different formulas + parameters for each! - -# **p100** is unaffected by our design and has amplitude of 5 -p1 = LinearModelComponent(; - basis = p100(), - formula = @formula(0~1), - β = [5] - ); - -# **n170** has a condition effect, faces are more negative than cars -n1 = LinearModelComponent(; - basis = n170(), - formula = @formula(0~1+condition), - β = [5,-3] - ); -# **p300** has a continuous effect, higher continuous values will result in larger P300's -p3 = LinearModelComponent(; - basis = p300(), - formula = @formula(0~1+continuous), - β = [5,1] - ); - -# now we can simply combine the components and simulate -components = [p1,n1,p3] -data,evts = simulate(MersenneTwister(1),design,[p1,n1,p3],UniformOnset(;width=0,offset=1000),PinkNoise()); - -## Analysis -# Let's check that everything worked out well, by using Unfold - -m = fit(UnfoldModel,Dict(Any=>(@formula(0~1+condition+continuous),firbasis(τ=[-0.1,1],sfreq=100,name="basis"))),evts,data); - -# first the "pure" beta/linear regression parameters -plot_erp(coeftable(m)) - -# and now beautifully visualized as marginal betas / predicted ERPs -plot_erp(effects(Dict(:condition=>["car","face"],:continuous=>-5:5),m); - mapping=(:color=>:continuous,linestyle=:condition,group=:continuous), - extra=(;categoricalColor=false)) diff --git a/joss_paper/paper.bib b/joss_paper/paper.bib new file mode 100644 index 0000000..404e716 --- /dev/null +++ b/joss_paper/paper.bib @@ -0,0 +1,181 @@ +@article{Julia-2017, + title={Julia: A fresh approach to numerical computing}, + author={Bezanson, Jeff and Edelman, Alan and Karpinski, Stefan and Shah, Viral B}, + journal={SIAM {R}eview}, + volume={59}, + number={1}, + pages={65--98}, + year={2017}, + publisher={SIAM}, + doi={10.1137/141000671}, + url={https://epubs.siam.org/doi/10.1137/141000671} +} + +@article{JSSv107i04, + title={DataFrames.jl: Flexible and Fast Tabular Data in Julia}, + volume={107}, + url={https://www.jstatsoft.org/index.php/jss/article/view/v107i04}, + doi={10.18637/jss.v107.i04}, + abstract={DataFrames.jl is a package written for and in the Julia language offering flexible and efficient handling of tabular data sets in memory. Thanks to Julia’s unique strengths, it provides an appealing set of features: Rich support for standard data processing tasks and excellent flexibility and efficiency for more advanced and non-standard operations. We present the fundamental design of the package and how it compares with implementations of data frames in other languages, its main features, performance, and possible extensions. We conclude with a practical illustration of typical data processing operations.}, + number={4}, + journal={Journal of Statistical Software}, + author={Bouchet-Valat, Milan and Kamiński, Bogumił}, + year={2023}, + pages={1--32} +} + +@misc{Distributions.jl-2019, + author = {Dahua Lin and + John Myles White and + Simon Byrne and + Douglas Bates and + Andreas Noack and + John Pearson and + Alex Arslan and + Kevin Squire and + David Anthoff and + Theodore Papamarkou and + Mathieu Besançon and + Jan Drugowitsch and + Moritz Schauer and + other contributors}, + title = {{JuliaStats/Distributions.jl: a Julia package for probability distributions and associated functions}}, + month = jul, + year = 2019, + doi = {10.5281/zenodo.2647458}, + url = {https://doi.org/10.5281/zenodo.2647458} +} + +@article{JSSv098i16, +author = {Mathieu Besançon and Theodore Papamarkou and David Anthoff and Alex Arslan and Simon Byrne and Dahua Lin and John Pearson}, +title = {Distributions.jl: Definition and Modeling of Probability Distributions in the JuliaStats Ecosystem}, +journal = {Journal of Statistical Software}, +volume = {98}, +number = {16}, +year = {2021}, +keywords = {Julia; distributions; modeling; interface; mixture; KDE; sampling; probabilistic programming; inference}, +issn = {1548-7660}, +pages = {1--30}, +doi = {10.18637/jss.v098.i16}, +url = {https://www.jstatsoft.org/v098/i16} +} + +@software{simon_kornblith_2023_8344531, +author = {Simon Kornblith and +Galen Lynch and +Martin Holters and +João Felipe Santos and +Spencer Russell and +Jay Kickliter and +Jeff Bezanson and +Gudmundur Adalsteinsson and +Alex Arslan and +Ryuichi Yamamoto and +jordancluts and +Matti Pastell and +Tony Kelman and +Ben Arthur and +Tom Krauss and +HDictus and +Hamza El-Saawy and +Jared Kofron and +Eric Hanson and +Rob Luke and +Viral B. Shah and +Aleco Kastanos and +Bill and +Clemens and +Elliot Saba and +ibadr and +Jake Bolewski and +Jordan Smith}, +title = {JuliaDSP/DSP.jl: v0.7.9}, +month = sep, +year = 2023, +publisher = {Zenodo}, +version = {v0.7.9}, +doi = {10.5281/zenodo.8344531}, +url = {https://doi.org/10.5281/zenodo.8344531} +} + +@article{DanischKrumbiegel2021, + doi = {10.21105/joss.03349}, + url = {https://doi.org/10.21105/joss.03349}, + year = {2021}, + publisher = {The Open Journal}, + volume = {6}, + number = {65}, + pages = {3349}, + author = {Simon Danisch and Julius Krumbiegel}, + title = {{Makie.jl}: Flexible high-performance data visualization for {Julia}}, + journal = {Journal of Open Source Software} +} + +@software{douglas_bates_2023_10268806, + author = {Douglas Bates and + Phillip Alday and + Dave Kleinschmidt and + José Bayoán Santiago Calderón and + Likan Zhan and + Andreas Noack and + Milan Bouchet-Valat and + Alex Arslan and + Tony Kelman and + Antoine Baldassari and + Benedikt Ehinger and + Daniel Karrasch and + Elliot Saba and + Jacob Quinn and + Michael Hatherly and + Morten Piibeleht and + Patrick Kofod Mogensen and + Simon Babayan and + Tim Holy and + Yakir Luc Gagnon and + Yoni Nazarathy}, + title = {JuliaStats/MixedModels.jl: v4.22.3}, + month = dec, + year = 2023, + publisher = {Zenodo}, + version = {v4.22.3}, + doi = {10.5281/zenodo.10268806}, + url = {https://doi.org/10.5281/zenodo.10268806} +} + +@software{phillip_alday_2022_7407741, + author = {Phillip Alday and + Douglas Bates and + Lisa DeBruine and + José Bayoán Santiago Calderón, PhD and + June Choe and + Lisa Schwetlick}, + title = {RePsychLing/MixedModelsSim.jl: v0.2.7}, + month = dec, + year = 2022, + publisher = {Zenodo}, + version = {v0.2.7}, + doi = {10.5281/zenodo.7407741}, + url = {https://doi.org/10.5281/zenodo.7407741} +} + +@article{Ehinger_Unfold_an_integrated, +author = {Ehinger, Benedikt}, +doi = {https://doi.org/10.7717/peerj.7838}, +journal = {PeerJ}, +title = {{Unfold: an integrated toolbox for overlap correction, non-linear modeling, and regression-based EEG analysis}} +} + +@software{mikheev_2023_10235220, + author = {Mikheev, Vladimir and + Döring, Sören and + Gärtner, Niklas and + Baumgartner, Daniel and + Ehinger, Benedikt}, + title = {UnfoldMakie}, + month = nov, + year = 2023, + publisher = {Zenodo}, + version = {v0.4.0}, + doi = {10.5281/zenodo.10235220}, + url = {https://doi.org/10.5281/zenodo.10235220} +} \ No newline at end of file diff --git a/joss_paper/paper.md b/joss_paper/paper.md new file mode 100644 index 0000000..05c5ad6 --- /dev/null +++ b/joss_paper/paper.md @@ -0,0 +1,110 @@ + +--- +title: 'Unfold.jl: Regression models for Cognitive Neuroscience' +tags: + - Julia + - EEG + - neuroimaging +authors: + - name: Benedikt V. Ehinger + orcid: 0000-0002-6276-3332 + equal-contrib: false + affiliation: "1, 2" # (Multiple affiliations must be quoted) + +affiliations: + - name: Stuttgart Center for Simulation Science, University of Stuttgart, Germany + index: 1 + - name: Institute for Visualization and Interactive Systems, University of Stuttgart, Germany + index: 2 +date: 06 November 2023 +bibliography: paper.bib + + +--- +# Summary + +UnfoldSim.jl is used to simulate multivariate timeseries, with a focus on EEG, especially event-related potentials (ERPs). The user provides four ingredients: 1) an experimental design, with both categorical and continuous variables, 2) event basis functions specified via linear or hierarchical models, 3) an inter-event onset distribution, and 4) a noise specification. Unfold.jl simulates continuous EEG signals with potentially overlapping events. Multi-channel support via EEG-forward models is available as well. UnfoldSim.jl is modular, providing intuitive entrance points for custom requirements. For instance, support for other modalities, e.g. single-voxel fMRI or pupil dilation signals is easily provided. + +# Statement of Need +Simulated EEG data is necessary to test preprocessing and analysis tools, to illustrate issues and to test functions. While other simulation tools exist, they are dominantly based in Matlab (below) and they do not adress our unique challenges and requirements. In our work (e.g. Ehinger & Dimigen 2019), we focus in regression-based deconvolution of ERPs (Smith & Kutas 2015 ab). In short, multiple-regression is used for linear overlap correction, non-linear (hierarchical) effects fitting (similar to GAMMs, Wood et al), often applied to use-cases where events overlap in time (e.g. stimulus and button press or consecutive eye-fixations). The tools used up to now (e.g. mTRF-toolbox, Unfold, Unfold.jl, fitgrid, - add citations!) + +# Functionality +The toolbox provides four abstract components: AbstractDesign, AbstractComponent, AbstractOnset and AbstractNoise. + +## Concrete Designs +Currently we support a single, and a multi-subject design. They are used to generate an experimental design containing the conditions and levels of all predictors. Randomisation is possible via a user-defined function which is applied after design-generation. Designs could be nested, like our RepeatDesign which simply repeats the generated design multiple times. + +## Concrete Components +We provide a LinearModelComponent and a MixedModelComponent for multi-subject simulation respectively. For the components model-formulae, fixed-effects ($\betas$) and random effects need to be specified. Further the coding-schema can be provided following StatsModels.jl. The component MultichannelComponent can be applied to any component and allows for projecting the simulated source-component to multi-channel electrodes via a headmodel. Using Artifacts.jl we provide on-demand access to the Hartmut (cite) model. + +## Concrete Noise +We provide different noise-types "White","Red" and "Pink", but also an exponentially declining Autoregressive noise type. + +## Concerete Onsets +The onsets define the distance between events for the continuous EEG. Currently UniformOnset and LogNormalOnset are implemented. + + + +# Related tools +Not many toolboxes for simulating EEG data exist. Nearly all toolboxes have been developed in MATLAB (e.g. EEGg - Vaziri, SimMEEG - Herdman 2021, SEED-G-Toolbox Anzolin 2021) but are largely abandoned. MNE-Python provides basic tutorials to simulate EEG data as well, but no dedicated functionality. Thus we highlight in the following two other excellent matlab-based tools: Brainstorm (Tadel 2021) and SEREEGA (Krol 2017). Both toolboxes are based in matlab and provide forward-simulation of EEG signals. Brainstorm especially excells at the visualization of the forward-model, and provides interesting capabilities to generate ERPs based on phase-aligned oscillations. SEREEGA provides the most complete simulation capabilities with a greater focus on ERP-component simulation, tools for benchmarking like signal-to-noise specification, and more realistic noise simulation (e.g. via random sources). + +In relation to these tools, Unfold uniquely focuses on the regression-ERP aspect, providing functions to simulate multi-condition experiments, uniquely allows to model hierarchical, that is, multi-subject EEG datasets, and offers unique support to model continuous EEG data with overlap events. + +Due to its different focus, UnfoldSim.jl currently lacks advanced visualizations of leadfields and does not provide any tools for simulating oscillations or phase-based simulation of ERPs. + + + +# Other notes +SEEREGA - Matlab, best in class, no continuous data +BRAINSTORM - Matlab, good for forward-model simulation +EEGg - Vaziri - Matlab, very limited functionality +SimMEEG - Herdman 2021 - Matlab, Discontinued +SEED-G-toolbox - Anzolin 2021 - Matlab, Discontinued, + +MNE-python - matlab, only basic tutorials + +## JOSS Stuff +Single dollars ($) are required for inline mathematics e.g. $f(x) = e^{\pi/x}$ + +Double dollars make self-standing equations: + +$$\Theta(x) = \left\{\begin{array}{l} +0\textrm{ if } x < 0\cr +1\textrm{ else} +\end{array}\right.$$ + +You can also use plain \LaTeX for equations +\begin{equation}\label{eq:fourier} +\hat f(\omega) = \int_{-\infty}^{\infty} f(x) e^{i\omega x} dx +\end{equation} +and refer to \autoref{eq:fourier} from text. + +# Citations + +Citations to entries in paper.bib should be in +[rMarkdown](http://rmarkdown.rstudio.com/authoring_bibliographies_and_citations.html) +format. + +If you want to cite a software repository URL (e.g. something on GitHub without a preferred +citation) then you can do it with the example BibTeX entry below for @fidgit. + +For a quick reference, the following citation commands can be used: +- `@author:2001` -> "Author et al. (2001)" +- `[@author:2001]` -> "(Author et al., 2001)" +- `[@author1:2001; @author2:2001]` -> "(Author1 et al., 2001; Author2 et al., 2002)" + +# Figures + +Figures can be included like this: +![Caption for example figure.\label{fig:example}](figure.png) +and referenced from text using \autoref{fig:example}. + +Figure sizes can be customized by adding an optional second parameter: +![Caption for example figure.](figure.png){ width=20% } + +# Acknowledgements + +We acknowledge contributions from Brigitta Sipocz, Syrtis Major, and Semyeong +Oh, and support from Kathryn Johnston during the genesis of this project. + +# References diff --git a/paper.md b/paper.md deleted file mode 100644 index 8b13789..0000000 --- a/paper.md +++ /dev/null @@ -1 +0,0 @@ - diff --git a/src/UnfoldSim.jl b/src/UnfoldSim.jl index 6541a73..e8756fe 100644 --- a/src/UnfoldSim.jl +++ b/src/UnfoldSim.jl @@ -13,11 +13,14 @@ module UnfoldSim using LinearAlgebra using ToeplitzMatrices # for AR Expo. Noise "Circulant" using StatsModels - + using HDF5,Artifacts,FileIO + + using LinearAlgebra # headmodel import DSP.hanning import Base.length import Base.size + import Base.show include("types.jl") include("design.jl") include("component.jl") @@ -25,6 +28,7 @@ module UnfoldSim include("simulation.jl") include("onset.jl") include("predefinedSimulations.jl") + include("headmodel.jl") include("helper.jl") include("bases.jl") @@ -61,4 +65,10 @@ module UnfoldSim # export bases export p100,n170,p300,n400,hrf + + # headmodel + export AbstractHeadmodel,Hartmut,headmodel,leadfield,orientation,magnitude + + # multichannel + export MultichannelComponent end diff --git a/src/component.jl b/src/component.jl index 11b76a2..f147d69 100644 --- a/src/component.jl +++ b/src/component.jl @@ -58,6 +58,54 @@ LinearModelComponent(; end +""" +Wrapper for an `AbstractComponent` to project it to multiple target-channels via `projection`. optional adds `noise` to the source prior to projection. +""" +@with_kw struct MultichannelComponent <:AbstractComponent + component::AbstractComponent + projection::AbstractVector + noise::AbstractNoise # optional +end + +MultichannelComponent(c::AbstractComponent,p) = MultichannelComponent(c::AbstractComponent,p,NoNoise()) + +function MultichannelComponent(c::AbstractComponent,p::Pair{<:AbstractHeadmodel,String},n::AbstractNoise) + ix = closest_src(p[1],p[2]) + mg = magnitude(p[1]) + return MultichannelComponent(c,mg[:,ix],n) +end +Base.length(c::MultichannelComponent) = length(c.component) + +""" +Returns the number of channels. By default = 1 +""" +n_channels(c::AbstractComponent) = 1 + +""" +for `MultichannelComponent` returns the length of the projection vector +""" +n_channels(c::MultichannelComponent) = length(c.projection) + + +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" + return all_channels[1] +end + +function simulate(rng,c::MultichannelComponent,design::AbstractDesign) + y = simulate(rng,c.component,design) + + for tr = 1:size(y,2) + y[:,tr] .= y[:,tr] .+ gen_noise(rng,c.noise,size(y,1)) + end + + y_proj = kron(y,c.projection) + return reshape(y_proj,length(c.projection),size(y)...,) +end + + + Base.length(c::AbstractComponent) = length(c.basis) maxlength(c::Vector{AbstractComponent}) = maximum(length.(c)) diff --git a/src/design.jl b/src/design.jl index a22c05f..1b8f0bf 100644 --- a/src/design.jl +++ b/src/design.jl @@ -10,11 +10,11 @@ tipp: check the resulting dataframe using `generate(design)` ```julia # declaring same condition both sub-between and item-between results in a full between subject/item design -design = MultiSubjectDesignjectDesign(; - n_items=10, +design = MultiSubjectDesign(; + n_items = 10, n_subjects = 30, - subjects_between=Dict(:cond=>["levelA","levelB"]), - items_between =Dict(:cond=>["levelA","levelB"]), + subjects_between = Dict(:cond => ["levelA", "levelB"]), + items_between = Dict(:cond => ["levelA", "levelB"]), ); ``` """ @@ -24,7 +24,7 @@ design = MultiSubjectDesignjectDesign(; subjects_between = nothing items_between = nothing both_within = nothing - tableModifyFun = x->x; # can be used to sort, or x->shuffle(rng,x) + tableModifyFun = x -> x # can be used to sort, or x->shuffle(rng,x) end @@ -39,8 +39,8 @@ To increase the number of repetitions simply use `RepeatDesign(SingleSubjectDesi tipp: check the resulting dataframe using `generate(design)` """ @with_kw struct SingleSubjectDesign <: AbstractDesign - conditions = nothing - tableModifyFun = x->x; + conditions = nothing + tableModifyFun = x -> x end @@ -57,13 +57,14 @@ julia> d = SingleSubjectDesign(;conditions= Dict(:A=>nlevels(5),:B=>nlevels(2))) julia> generate(d) """ function generate(expdesign::SingleSubjectDesign) - # we get a Dict(:A=>["1","2"],:B=>["3","4"]), but needed a list - # of named tuples for MixedModelsSim.factorproduct function. - evts = factorproduct( - ((;k=>v) for (k,v) in pairs(expdesign.conditions))...) |> DataFrame - - # by default does nothing - return expdesign.tableModifyFun(evts) + # we get a Dict(:A=>["1","2"],:B=>["3","4"]), but needed a list + # of named tuples for MixedModelsSim.factorproduct function. + evts = + factorproduct(((; k => v) for (k, v) in pairs(expdesign.conditions))...) |> + DataFrame + + # by default does nothing + return expdesign.tableModifyFun(evts) end """ @@ -80,32 +81,33 @@ julia> d = MultiSubjectDesign(;n_subjects = 10,n_items=20,both_within= Dict(:A=> julia> generate(d) """ function generate(expdesign::MultiSubjectDesign) - #generate(expdesign::AbstractDesign) = generate(MersenneTwister(1),expdesign) - - # check that :dv is not in any condition - allconditions = [expdesign.subjects_between,expdesign.items_between,expdesign.both_within] - @assert :dv ∉ keys(merge(allconditions[.!isnothing.(allconditions)]...)) "due to technical limitations in MixedModelsSim.jl, `:dv` cannot be used as a factorname" - - - data = DataFrame( - MixedModelsSim.simdat_crossed( - expdesign.n_subjects, - expdesign.n_items, - subj_btwn=expdesign.subjects_between, - item_btwn=expdesign.items_between, - both_win=expdesign.both_within - ) - ) - rename!(data,:subj => :subject) - select!(data,Not(:dv)) # remove the default column from MixedModelsSim.jl - we don't need it in UnfoldSim.jl - # by default does nothing - data = expdesign.tableModifyFun(data) - - # sort by subject - data = sort!(data,(order(:subject))) - - return data - + #generate(expdesign::AbstractDesign) = generate(MersenneTwister(1),expdesign) + + # check that :dv is not in any condition + allconditions = + [expdesign.subjects_between, expdesign.items_between, expdesign.both_within] + @assert :dv ∉ keys(merge(allconditions[.!isnothing.(allconditions)]...)) "due to technical limitations in MixedModelsSim.jl, `:dv` cannot be used as a factorname" + + + data = DataFrame( + MixedModelsSim.simdat_crossed( + expdesign.n_subjects, + expdesign.n_items, + subj_btwn = expdesign.subjects_between, + item_btwn = expdesign.items_between, + both_win = expdesign.both_within, + ), + ) + rename!(data, :subj => :subject) + select!(data, Not(:dv)) # remove the default column from MixedModelsSim.jl - we don't need it in UnfoldSim.jl + # by default does nothing + data = expdesign.tableModifyFun(data) + + # sort by subject + data = sort!(data, (order(:subject))) + + return data + end @@ -131,17 +133,18 @@ design = RepeatDesign(designOnce,4); ``` """ @with_kw struct RepeatDesign{T} <: AbstractDesign - design::T - repeat::Int = 1 + design::T + repeat::Int = 1 end -function UnfoldSim.generate(design::RepeatDesign) - df = map(x->generate(design.design),1:design.repeat) |>x->vcat(x...) - if isa(design.design,MultiSubjectDesign) - sort!(df,[:subject]) - end - return df - +function UnfoldSim.generate(design::RepeatDesign) + df = map(x -> generate(design.design), 1:design.repeat) |> x -> vcat(x...) + if isa(design.design, MultiSubjectDesign) + sort!(df, [:subject]) + end + return df + 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::RepeatDesign{MultiSubjectDesign}) = + size(design.design) .* (design.repeat, 1) +Base.size(design::RepeatDesign{SingleSubjectDesign}) = size(design.design) .* design.repeat diff --git a/src/headmodel.jl b/src/headmodel.jl new file mode 100644 index 0000000..e046de1 --- /dev/null +++ b/src/headmodel.jl @@ -0,0 +1,109 @@ + + +struct Hartmut <: AbstractHeadmodel + artefactual + cortical + electrodes +end + +function Base.show(io::IO,h::Hartmut) + src_label = h.cortical["label"] + ele_label = h.electrodes["label"] + art_label = h.artefactual["label"] + + println(io,"""HArtMuT-Headmodel + $(length(ele_label)) electrodes: ($(ele_label[1]),$(ele_label[2])...) - hartmut.electrodes + $(length(src_label)) source points: ($(src_label[1]),...) - hartmut.cortical + $(length(art_label)) source points: ($(art_label[1]),...) - hartmut.artefactual + + In addition to UnfoldSim.jl please cite: + $(hartmut_citation()) + """) +end + +"Returns citation-string for HArtMuT" +hartmut_citation() = "HArtMuT: Harmening Nils, Klug Marius, Gramann Klaus and Miklody Daniel - 10.1088/1741-2552/aca8ce" + +"Returns the leadfield" +leadfield(hart::Hartmut;type="cortical") = type == "cortical" ? hart.cortical["leadfield"] : hart.artefactual["leadfield"] +orientation(hart::Hartmut;type="cortical") = type == "cortical" ? hart.cortical["orientation"] : hart.artefactual["orientation"] + + +""" +Load a headmodel, using Artifacts.jl automatically downloads the required files + +Currently only `type="hartmut"` is implemented +""" +function headmodel(;type="hartmut") + if type == "hartmut" + println("""Please cite: $(hartmut_citation())""") + path = joinpath(artifact"hartmut", "hartmut.h5") + h = h5open(path) + + + weirdchan = ["Nk1","Nk2","Nk3","Nk4"] + ## getting index of these channels from imported hartmut model data, exclude them in the topoplot + remove_indices = findall(l -> l ∈ weirdchan, h["electrodes"]|>read|>x->x["label"]); + + function sel_chan(x) + + if "leadfield" ∈ keys(x) + x["leadfield"] = x["leadfield"][Not(remove_indices),:,:].*10e3 # this scaling factor seems to generate potentials with +-1 as max + else + x["label"] = x["label"][Not(remove_indices)] + pos3d= x["pos"][Not(remove_indices),:] + pos3d = pos3d ./ (4*maximum(pos3d,dims=1)) + x["pos"] = pos3d + end + return x + end + headmodel = Hartmut( + h["artefacts"]|>read|>sel_chan, + h["cortical"]|>read|>sel_chan, + h["electrodes"]|>read|>sel_chan, + ) + else + error("unknown headmodel. currently only 'hartmut' allowed") + end + + return headmodel +end + +""" +Extracts magnitude of the orientation-including leadfield. + +By default uses the orientation specified in the headmodel + +Fallback: along the third dimension using `norm` - the maximal projection +""" +magnitude(headmodel::AbstractHeadmodel) = magnitude(leadfield(headmodel)) + +""" +Extract magnitude of 3-orientation-leadfield, +`type` (default: "perpendicular") => uses the provided source-point orientations - otherwise falls back to `norm` +""" +magnitude(headmodel::Hartmut;type="perpendicular") = type == "perpendicular" ? magnitude(leadfield(headmodel),orientation(headmodel)) : magnitude(leadfield(headmodel)) + + +function magnitude(lf::AbstractArray{T,3},orientation::AbstractArray{T,2}) where {T<:Real} + si = size(lf); + magnitude = fill(NaN,si[1:2]); + for e=1:si[1] + for s=1:si[2] + magnitude[e,s] = lf[e,s,:]' * orientation[s,:] + end + end + return magnitude +end + + +function magnitude(lf::AbstractArray{T,3}) where {T<:Real} + si = size(lf); + magnitude = fill(NaN,si[1:2]); + for e=1:si[1] + for s=1:si[2] + magnitude[e,s] = norm(lf[e,s,:]); + end + end + return magnitude +end diff --git a/src/helper.jl b/src/helper.jl old mode 100644 new mode 100755 index 66e97af..bab5704 --- a/src/helper.jl +++ b/src/helper.jl @@ -10,42 +10,78 @@ function padarray(arr::Vector, len::Int, val) end +# TODO: Transfer function to Unfold.jl """ Function to convert output similar to unfold (data, evts) """ -function convert(eeg, onsets, design;reshape=true) +function convert(eeg, onsets, design,n_ch,;reshape=true) evt = UnfoldSim.generate(design) + @debug size(eeg) if reshape + n_subj = length(size(design))==1 ? 1 : size(design)[2] + + if n_ch == 1 data = eeg[:,] + evt.latency = (onsets' .+ range(0,size(eeg,2)-1).*size(eeg,1) )'[:,] + elseif n_subj == 1 + data = eeg + @debug size(onsets) + evt.latency = onsets + else # multi subject + multi channel + data = eeg[:,:,] + evt.latency = (onsets' .+ range(0,size(eeg,3)-1).*size(eeg,2) )'[:,] + end else data = eeg end - if :d ∈ names(evt) - select!(evt, Not([:dv])) - end - return data,evt end + """ -Takes an array of 'm' target coordinate arrays (1-by-3) and a matrix (n-by-3) of all available positions, and returns an array of size 'm' containing the indices of the respective items in 'pos' that are nearest to each of the target coordinates. + closest_src(coords_list::AbstractVector{<:AbstractVector}, pos) + closest_src(coords::Vector{<:Real}, pos) + +Takes an array of 'm' target coordinate vector (size 3) (or vector of vectors) and a matrix (n-by-3) of all available positions, and returns an array of size 'm' containing the indices of the respective items in 'pos' that are nearest to each of the target coordinates. """ -function closest_srcs(coords_list, pos) - out = []; +closest_src(coords_list::AbstractVector{<:AbstractVector}, pos) = closest_src.(coords_list, Ref(pos)) + +function closest_src(coords::Vector{<:Real}, pos) s = size(pos); dist = zeros(s[1]); diff = zeros(s[2]); - for coords in coords_list for i=1:s[1] for j=1:s[2] diff[j] = pos[i,j] - coords[j]; end dist[i] = norm(diff); end - push!(out,findmin(dist)[2]) #retain only the index of the minimum difference. - end - return out + return findmin(dist)[2] + + end + +""" + closest_src(head::Hartmut,label::String) +Returns src-`ix` of the Headmodel `Hartmut` which is closest to the average of the `label`. + +!!! important + We use the average in eucledean space, but the cortex is a curved surface. In most cases they will not overlap. Ideally we would calculate the average on the surface, but this is a bit more complex to do (you'd need to calculate the vertices etc.) + +```julia +hartmut = headmodel() +pos = closest_src(hartmut=>"Left Middle Temporal Gyrus, posterior division") +``` +""" +function closest_src(head::Hartmut,label::String) + + pos = head.cortical["pos"] + ix = findall(head.cortical["label"] .== label) + @assert sum(ix)>0 """could not find label $label in hartmut.cortical["label"] - try unique(hartmut.cortical["label"]) for a list""" + + ix = UnfoldSim.closest_src(mean(pos[ix,:],dims=1)[1,:],pos) + return ix +end \ No newline at end of file diff --git a/src/predefinedSimulations.jl b/src/predefinedSimulations.jl index 4c20756..af0f0ad 100644 --- a/src/predefinedSimulations.jl +++ b/src/predefinedSimulations.jl @@ -7,30 +7,35 @@ predef_eeg(;kwargs...) = predef_eeg(MersenneTwister(1);kwargs...) # without rng predef_eeg(nsubjects::Int;kwargs...) = predef_eeg(MersenneTwister(1),nsubjects;kwargs...) # without rng always call same one """ -Generates a P1/N1/P3 complex. predef_eeg(;kwargs...) predef_eeg(rng;kwargs...) predef_eeg(rng,n_subjects;kwargs...) -In case of n_subjects, MixedModelComponents are generated - -Default params: - n_repeats=100 - tableModifyFun = 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, 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 - - # 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]), +Gene +rates a P1/N1/P3 complex. +In case `n_subjects` is defined - `MixedModelComponents`` are generated, else `LinearModelComponents` + +Most used kwargs is: `return_epoched=true` to ignore the overlap/onset bits and return already epoched data + +## Default params: + +. n_repeats=100 +- tableModifyFun = 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, 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 + +#### 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]), """ function predef_eeg(rng; # design @@ -74,12 +79,6 @@ function predef_eeg(rng::AbstractRNG,design::AbstractDesign,T::Type{<:AbstractCo return simulate(rng,design, components, onset, noise;kwargs...); end - -""" -predef_eeg(rng,n_subjects;kwargs...) -Runs predef_eeg(rng;kwargs...) n_subject times and concatenates the results. - -""" function predef_eeg(rng::AbstractRNG,n_subjects; # design n_items=100, @@ -105,7 +104,33 @@ function predef_eeg(rng::AbstractRNG,n_subjects; end """ -todo + + predef_2x2(rng::AbstractRNG;kwargs...) + +Most used kwargs is: `return_epoched=true` to ignore the overlap/onset bits and return already epoched data + +#### design +- `n_items`=100, +- `n_subjects`=1, +- `conditions` = Dict(:A=>["a_small","a_big"],:B=>["b_tiny","b_large"]), +- `tableModifyFun` = 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)), + +#### 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 + 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) """ diff --git a/src/simulation.jl b/src/simulation.jl old mode 100644 new mode 100755 index d81e77c..f0d0ff4 --- a/src/simulation.jl +++ b/src/simulation.jl @@ -2,12 +2,17 @@ # helper to move input ::Component to ::Vector{Component} Simulation(design::AbstractDesign,component::AbstractComponent,onset::AbstractOnset,noisetype::AbstractNoise) = Simulation(design,[component],onset,noisetype) +# by default no noise +Simulation(design::AbstractDesign,component,onset::AbstractOnset) = Simulation(design,component,onset,NoNoise()) + """ Simulate eeg data given a simulation design, effect sizes and variances make use of `return_epoched=true` to skip the Onset-calculation + conversion to continuous data and get the epoched data directly """ -simulate(rng,design, signal, onset, noise;kwargs...) = simulate(rng,Simulation(design, signal, onset, noise);kwargs...) + +simulate(rng, design::AbstractDesign, signal, onset::AbstractOnset, noise::AbstractNoise=NoNoise(); kwargs...) = simulate(rng, Simulation(design, signal, onset, noise); kwargs...) + function simulate(rng, simulation::Simulation;return_epoched::Bool=false) # unpacking fields @@ -15,57 +20,99 @@ function simulate(rng, simulation::Simulation;return_epoched::Bool=false) # create epoch data / erps erps = simulate(deepcopy(rng), components,simulation) +# @debug size(erps) + + n_subj = length(size(design))==1 ? 1 : size(design)[2] + n_trial = size(design)[1] + n_ch = n_channels(components) + + # create events data frame + events = UnfoldSim.generate(design) if !return_epoched # we only need to simulate onsets & pull everything together, if we # want a continuous EEG onsets = generate(deepcopy(rng),onset,simulation) - - # XXX todo: Separate Subjects in Time by adding offset to onsets!! + + # save the onsets in the events df + events.latency = onsets[:,] # combine erps with onsets maxlen = maxlength(components) max_length = Int(ceil(maximum(onsets))) .+ maxlen - n_subj = length(size(design))==1 ? 1 : size(design)[2] - n_trial = size(design)[1] - eeg = zeros(max_length,n_subj) + + eeg = zeros(n_ch,max_length,n_subj) + # not all designs have multiple subjects + if n_subj == 1; eeg = dropdims(eeg,dims=3); end + + # not all designs have multiple channels + if n_ch == 1; eeg = dropdims(eeg,dims=1); end + - for s in 1:n_subj - for i in 1:n_trial - one_onset = onsets[CartesianIndex(i, s)] - @views eeg[one_onset:one_onset+maxlen-1,s] .+= erps[:, (s-1)*n_trial+i] - end - end + for e in 1:n_ch + for s in 1:n_subj + for i in 1:n_trial + one_onset = onsets[CartesianIndex(i, s)] + adderp!(eeg,erps,e,s,one_onset:one_onset+maxlen-1,(s-1)*n_trial+i) + end + end + end else eeg = erps - onsets = [] # this is still a bit ugly + end + + add_noise!(deepcopy(rng),noisetype,eeg) - return convert(eeg,onsets,design;reshape=!return_epoched) + return eeg, events end """ -Simulates erp data given the specified parameters +Helper function to add inplace the erps to the EEG, but for both 2D (1 channel) and 3D (X channel case) """ -function simulate(rng, components::Vector{<:AbstractComponent},simulation::Simulation) +function adderp!(eeg,erps::Vector,e,s,tvec,erpvec) + @views eeg[tvec] .+= erps[:, erpvec] +end +function adderp!(eeg,erps::Matrix,e,s,tvec,erpvec) + @views eeg[tvec,s] .+= erps[:, erpvec] +end +function adderp!(eeg,erps::AbstractArray,e,s,tvec,erpvec) + @views eeg[e,tvec,s] .+= erps[e,:, erpvec] +end - epoch_data = zeros(maxlength(components), length(simulation.design)) - # Simulate each component - for c in components - # add them up +""" +Simulates erp data given the specified parameters +""" +function simulate(rng, 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 - @views epoch_data[1:length(c),:] .+= simulate(rng,c,simulation) + for c in components + simulateandadd!(epoch_data,c,simulation,rng) end return epoch_data end +function simulateandadd!(epoch_data::AbstractMatrix,c,simulation,rng) + @debug "matrix" + @views epoch_data[1:length(c),:] .+= simulate(rng,c,simulation) +end +function simulateandadd!(epoch_data::AbstractArray,c,simulation,rng) + @debug "3D Array" + @views epoch_data[:,1:length(c),:] .+= simulate(rng,c,simulation) +end + + function add_noise!(rng,noisetype::AbstractNoise,eeg) diff --git a/src/types.jl b/src/types.jl index e257259..fadfd6e 100644 --- a/src/types.jl +++ b/src/types.jl @@ -8,7 +8,7 @@ abstract type AbstractComponent end # find other types in onset.jl and noise.jl # and in design.jl and component.jl - +abstract type AbstractHeadmodel end struct Simulation diff --git a/test/Project.toml b/test/Project.toml index e6b8ec6..a12cea7 100644 --- a/test/Project.toml +++ b/test/Project.toml @@ -6,3 +6,4 @@ Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c" StableRNGs = "860ef19b-820b-49d6-a774-d7a799459cd3" Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2" Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" +UnfoldSim = "ed8ae6d2-84d3-44c6-ab46-0baf21700804" diff --git a/test/artifactheadmodel.jl b/test/artifactheadmodel.jl new file mode 100644 index 0000000..6e9b6b8 --- /dev/null +++ b/test/artifactheadmodel.jl @@ -0,0 +1,50 @@ +using Base: make_atomicreplace +using MAT +using HDF5 + + +dlpath = download("https://github.com/harmening/HArtMuT/raw/main/HArtMuTmodels/HArtMuT_NYhead_small.mat"); +file = matopen(dlpath); +hartmut = read(file, "HArtMuT"); +close(file) + +fname = tempname(); # temporary file + +fid = h5open(fname, "w") + + +label = string.(hartmut["electrodes"]["label"][:,1]) +chanpos = Float64.(hartmut["electrodes"]["chanpos"]) + +cort_label = string.(hartmut["cortexmodel"]["labels"][:,1]) +cort_orient = hartmut["cortexmodel"]["orientation"] +cort_leadfield = hartmut["cortexmodel"]["leadfield"] +cort_pos = hartmut["cortexmodel"]["pos"] + + +art_label = string.(hartmut["artefactmodel"]["labels"][:,1]) +art_orient = hartmut["artefactmodel"]["orientation"] +art_leadfield = hartmut["artefactmodel"]["leadfield"] +art_pos = hartmut["artefactmodel"]["pos"] + +e = create_group(fid, "electrodes") +e["label"] = label +e["pos"] = chanpos +c = create_group(fid, "cortical") +c["label"] = cort_label +c["orientation"] = cort_orient +c["leadfield"] = cort_leadfield +c["pos"] = cort_pos +a = create_group(fid, "artefacts") +a["label"] = art_label +a["orientation"] = art_orient +a["leadfield"] = art_leadfield +a["pos"] = art_pos + +close(fid) +mkdir(joinpath(tempdir(),"artifact")) +mv(fname,joinpath(tempdir(),"artifact","hartmut.h5")) +using ArtifactUtils +artifact_id = artifact_from_directory("/tmp/artifact") +gist = upload_to_gist(artifact_id) +add_artifact!("Artifacts.toml", "hartmut", gist) diff --git a/test/headmodel.jl b/test/headmodel.jl new file mode 100644 index 0000000..7dcbab0 --- /dev/null +++ b/test/headmodel.jl @@ -0,0 +1,23 @@ +hart = UnfoldSim.headmodel() +@testset "hartmut" begin + @test length(hart.electrodes["label"]) == 231-4 + @test length(hart.cortical["label"]) == 2004 + @test length(hart.artefactual["label"]) == 4260 +end +@testset "leadfield/magnitude" begin + lf = leadfield(hart) + @test isa(lf,AbstractArray) + mg = magnitude(hart,type="norm") + mg_lf = magnitude(lf) + @test mg == mg_lf + + or = orientation(hart) + + @test isa(lf,AbstractArray) + mg_man = magnitude(lf,or) + mg = magnitude(hart) + @test mg_man == mg + @test size(mg) == (231-4,2004) + + +end \ No newline at end of file diff --git a/test/helper.jl b/test/helper.jl index 71a3c4b..846d8cd 100644 --- a/test/helper.jl +++ b/test/helper.jl @@ -8,4 +8,11 @@ @test padarray([2,2],2,-1) == [2,2,-1,-1] end + + @testset "closest_src" begin + @test UnfoldSim.closest_src([0,0,1],[0 0 0.5; -1 -1 -1; -3 -3 -3]) == 1 + @test UnfoldSim.closest_src([[0,0,1],[-1,-1,0]],[0 0 0.5; -1 -1 -1; -3 -3 -3]) == [1,2] + + + end end \ No newline at end of file diff --git a/test/multichannel.jl b/test/multichannel.jl new file mode 100644 index 0000000..ab58b1d --- /dev/null +++ b/test/multichannel.jl @@ -0,0 +1,41 @@ +using UnfoldSim +using StableRNGs +using Test +#--- +hart = UnfoldSim.headmodel() +mg = magnitude(hart) + + +design = SingleSubjectDesign(conditions=Dict(:condA=>["levelA","levelB"])) |> x->RepeatDesign(x,5); +signal = LinearModelComponent(;basis=[0,1,2,3,0],formula = @formula(0~1+condA),β = [2,5]); +onset = UniformOnset(;width=20,offset=4); + + +@testset "multichannel" begin + mc = UnfoldSim.MultichannelComponent(signal,[-2,-1,0,1,2,3,4]) + s = simulate(StableRNG(1),mc,design) + @test size(s) == (7,5,10) # 7 chanels, 5 timepoints, 10 trials + @ŧest all(s[3,:,:] .== 0) # 3 is 0 in the projection + @test s[4,:,3] == [0,1,2,3,0].*2 # basis * β[1] + + # test different projection size, should error + mcA = UnfoldSim.MultichannelComponent(signal,[-2,-1,0,1,2,3,4]) + mcB = UnfoldSim.MultichannelComponent(signal,[-2,-1,0,1,2,3,4,5]) + @test_throws AssertionError UnfoldSim.n_channels([mcA,mcB]) + @test UnfoldSim.n_channels([mcA,mcA]) == 7 + + simulation = Simulation(design, mc, onset, NoNoise()); + data,events = simulate(StableRNG(1),simulation) + + @test all(data[:,events.latency[1]+1] .≈ [-4. -2. 0 2. 4. 6. 8.]') + + @test data[1,:] == data[2,:]*2 + @test all(data[3,:] .== 0) +end + +@testset "multichannel with helper" begin + + mc = UnfoldSim.MultichannelComponent(signal, hart=>"Left Middle Temporal Gyrus, posterior division") + data,events = simulate(StableRNG(1),design, mc, onset, NoNoise()) + @test size(data,1) == (231 - 4) +end \ No newline at end of file diff --git a/test/simulation.jl b/test/simulation.jl index 2163c3e..70248e8 100644 --- a/test/simulation.jl +++ b/test/simulation.jl @@ -1,3 +1,102 @@ @testset "simulation" begin + ## Define elements for the simulation + + # Define experimental factors + conditions = Dict(:cond => ["A", "B"]) + + # Create design for one subject with 20 trials (10 in each of the two factor levels) + repetitions = 10 + design_single_subject = SingleSubjectDesign(; + conditions = conditions) |> x -> RepeatDesign(x, repetitions) + + # Create design for multiple subjects with conditions as a between-items factor + n_subjects = 5 + n_items = 4 + design_multiple_subjects = MultiSubjectDesign(;n_subjects=n_subjects, n_items=n_items, items_between = conditions) + + # Linear component for the single-subject simulation + signal_linear = LinearModelComponent(; + basis = p100(), + formula = @formula(0 ~ 1 + cond), + β = [1, 0.5]) + + # Mixed-model component for the multi-subject simulation + signal_mixed = MixedModelComponent(; + basis = p100(), + formula = @formula(0 ~ 1 + cond + (1 + cond|subject)), + β = [1, 0.5], + σs = Dict(:subject => [0.2,0.1])) + + # Define headmodel and MultichannelComponent for multi-channel cases + hartmut_model = headmodel(type="hartmut") + signal_linear_multichannel = MultichannelComponent(signal_linear, hartmut_model => "Left Central Opercular Cortex") + signal_mixed_multichannel = MultichannelComponent(signal_mixed, hartmut_model => "Left Central Opercular Cortex") + + # Overlap since offset