Merge pull request #29 from unfoldtoolbox/multichannel

multichannel

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"

Project.toml

@@ -1,31 +1,41 @@
 name = "UnfoldSim"
 uuid = "ed8ae6d2-84d3-44c6-ab46-0baf21700804"
-authors = ["Benedikt Ehinger","Luis Lips", "Judith Schepers","Maanik Marathe"]
+authors = ["Benedikt Ehinger", "Luis Lips", "Judith Schepers", "Maanik Marathe"]
 version = "0.1.6"
 
 [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"
 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"

docs/Project.toml

@@ -5,6 +5,7 @@ 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"
 Parameters = "d96e819e-fc66-5662-9728-84c9c7592b0a"
 StableRNGs = "860ef19b-820b-49d6-a774-d7a799459cd3"

docs/make.jl

@@ -34,13 +34,15 @@ authors="Luis Lips, Benedikt Ehinger, Judith Schepers",
                 "Poweranalysis" => "literate/tutorials/poweranalysis.md",
         ],
         "Reference"=>[
-                "NoiseTypes" => "./literate/reference/noisetypes.md",
+                "Toolbox Overview" =>"./literate/reference/overview.md",
+                "NoiseTypes" =>      "./literate/reference/noisetypes.md",
                 "ComponentBasisTypes" => "./literate/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"
+                "New Duration/Shift-dependent Component" => "./literate/HowTo/newComponent.md",
+                "Multi Channel Data" =>"./literate/HowTo/multichannel.md",
         ],
         "DocStrings" => "api.md",
     ],

docs/src/literate/HowTo/multichannel.jl

@@ -0,0 +1,49 @@
+
+using UnfoldSim
+using UnfoldMakie
+using CairoMakie
+using DataFrames
+using Random
+
+## Specifcy 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]);
+
+# 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 + Plotting
+# 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, NoNoise())
+size(data)
+
+# Let's plot using Butterfly & Topoplot
+# first we convert the electrodes to positions usable in TopoPlots.jl
+pos3d = hart.electrodes["pos"]
+pos3d = pos3d ./ (4*maximum(pos3d,dims=1))#.+0.5  # standardize
+pos2d = to_positions(pos3d')
+pos2d = [Point2f(p[1]+0.5,p[2]+0.5) for p in pos2d]
+
+# let's 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
+

docs/src/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)
+
+

src/UnfoldSim.jl

@@ -13,18 +13,22 @@ 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")
    include("noise.jl")
    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

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

src/headmodel.jl

@@ -0,0 +1,107 @@
+
+
+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)]
+                x["pos"] = x["pos"][Not(remove_indices),:]
+            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