Skip to content

Commit

Permalink
Merge pull request #67 from unfoldtoolbox/format-everything
Browse files Browse the repository at this point in the history
Format everything
  • Loading branch information
behinger authored Jan 18, 2024
2 parents 50c7b37 + 53d3d6f commit 21e46a4
Show file tree
Hide file tree
Showing 30 changed files with 655 additions and 533 deletions.
1 change: 1 addition & 0 deletions .JuliaFormatter.toml
Original file line number Diff line number Diff line change
@@ -0,0 +1 @@
style = "default"
46 changes: 31 additions & 15 deletions docs/literate/HowTo/multichannel.jl
Original file line number Diff line number Diff line change
Expand Up @@ -9,32 +9,33 @@ using Random
# ## Specifying a design

# We are using a one-level design for testing here.
design = SingleSubjectDesign(conditions=Dict(:condA=>["levelA"]))
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]);
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])
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")
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);
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))
data, events =
simulate(MersenneTwister(1), design, [mc, mc2], onset, PinkNoise(noiselevel = 0.05))
size(data)


Expand All @@ -46,13 +47,28 @@ size(data)
# 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];
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))))
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

41 changes: 24 additions & 17 deletions docs/literate/HowTo/newComponent.jl
Original file line number Diff line number Diff line change
Expand Up @@ -11,43 +11,50 @@ sfreq = 100;

# ## 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))
design = UnfoldSim.SingleSubjectDesign(;
conditions = Dict(
:shift => rand(100) .* sfreq / 5,
:duration => 20 .+ rand(100) .* sfreq / 5,
),
)


# We also need a new AbstractComponent
struct TimeVaryingComponent <: AbstractComponent
basisfunction
maxlength
basisfunction::Any
maxlength::Any
end

# 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
function UnfoldSim.simulate(rng,c::TimeVaryingComponent,design::AbstractDesign)
function UnfoldSim.simulate(rng, c::TimeVaryingComponent, design::AbstractDesign)
evts = generate(design)
return c.basisfunction(evts,c.maxlength)
return c.basisfunction(evts, c.maxlength)
end

# finally, the actual function that does the shifting + duration
function basis_shiftduration(evts,maxlength)
function basis_shiftduration(evts, maxlength)
basis = hanning.(Int.(round.(evts.duration))) ## hanning as long as duration
if "shift" names(evts)
basis = padarray.(basis,Int.(round.(.-evts.shift)),0) ## shift by adding 0 in front
basis = padarray.(basis, Int.(round.(.-evts.shift)), 0) ## shift by adding 0 in front
end
## we should make sure that all bases have maxlength by appending / truncating
difftomax = maxlength .- length.(basis)
if any(difftomax.<0)
if any(difftomax .< 0)
@warn "basis longer than max length in at least one case. either increase maxlength or redefine function. Trying to truncate the basis"
basis[difftomax .>0] = padarray.(basis[difftomax .> 0],difftomax[difftomax .> 0],0)
basis[difftomax.>0] = padarray.(basis[difftomax.>0], difftomax[difftomax.>0], 0)
return [b[1:maxlength] for b in basis]
else
return padarray.(basis,difftomax,0)
return padarray.(basis, difftomax, 0)
end
end


erp = UnfoldSim.simulate(MersenneTwister(1),TimeVaryingComponent(basis_shiftduration,50),design)
plot_erpimage(hcat(erp...),sortvalues=generate(design).shift)


erp = UnfoldSim.simulate(
MersenneTwister(1),
TimeVaryingComponent(basis_shiftduration, 50),
design,
)
plot_erpimage(hcat(erp...), sortvalues = generate(design).shift)
12 changes: 6 additions & 6 deletions docs/literate/HowTo/newDesign.jl
Original file line number Diff line number Diff line change
Expand Up @@ -27,16 +27,16 @@ size(design::ImbalanceSubjectDesign) = (design.nTrials,);
# We need a type `generate(design::ImbalanceSubjectDesign)` function. This function should return the actual table as a `DataFrame`
function generate(design::ImbalanceSubjectDesign)
nA = Int(round.(design.nTrials .* design.balance))
nB = Int(round.(design.nTrials .* (1-design.balance)))
@assert nA + nB design.nTrials
levels = vcat(repeat(["levelA"],nA),repeat(["levelB"],nB))
return DataFrame(Dict(:condition=>levels))
nB = Int(round.(design.nTrials .* (1 - design.balance)))
@assert nA + nB design.nTrials
levels = vcat(repeat(["levelA"], nA), repeat(["levelB"], nB))
return DataFrame(Dict(:condition => levels))
end;

# Finally, we can test the function and see whether it returns a Design-DataFrame as we requested
design = ImbalanceSubjectDesign(;nTrials=6,balance=0.2)
design = ImbalanceSubjectDesign(; nTrials = 6, balance = 0.2)
generate(design)

# !!! 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.
# safe a RNG object in your struct and use it in your generate function.
61 changes: 30 additions & 31 deletions docs/literate/reference/basistypes.jl
Original file line number Diff line number Diff line change
Expand Up @@ -12,54 +12,53 @@ using StableRNGs
# ## EEG
# 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]
lines!(ax,b(),label=string(b))
scatter!(ax,b(),label=string(b))
ax = f[1, 1] = Axis(f)
for b in [p100, n170, p300, n400]
lines!(ax, b(), label = string(b))
scatter!(ax, b(), label = string(b))
end
axislegend(ax,merge=true)
axislegend(ax, merge = true)
f

# ## fMRI
# default hrf TR is 1. Get to know all your favourite shapes!
##--
f = Figure()
plotConfig = (:peak=>1:3:10,
:psUnder=>10:5:30,
:amplitude=>2:5,
:shift=>0:3:10,
:peak_width => 0.1:0.5:1.5,
:psUnder_width => 0.1:0.5:1.5,
)
plotConfig = (
:peak => 1:3:10,
:psUnder => 10:5:30,
:amplitude => 2:5,
:shift => 0:3:10,
:peak_width => 0.1:0.5:1.5,
:psUnder_width => 0.1:0.5:1.5,
)

for (ix,pl) = enumerate(plotConfig)
col = (ix-1)%3 +1
row = Int(ceil(ix/3))
ax = f[row,col] = Axis(f)
for (ix, pl) in enumerate(plotConfig)
col = (ix - 1) % 3 + 1
row = Int(ceil(ix / 3))

ax = f[row, col] = Axis(f)
cfg = collect(pl)
for k = cfg[2]
lines!(ax,UnfoldSim.hrf(;TR=0.1,(cfg[1]=>k,)...),label=string(k))
for k in cfg[2]
lines!(ax, UnfoldSim.hrf(; TR = 0.1, (cfg[1] => k,)...), label = string(k))
end
axislegend(string(cfg[1]);merge=true,)

axislegend(string(cfg[1]); merge = true)
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.
f = Figure()
plotConfig = (:n=>5:3:15,
:tmax=>0.5:0.2:1.1,
)
plotConfig = (:n => 5:3:15, :tmax => 0.5:0.2:1.1)

for (ix,pl) = enumerate(plotConfig)
ax = f[1,ix] = Axis(f)
for (ix, pl) in enumerate(plotConfig)
ax = f[1, ix] = Axis(f)
cfg = collect(pl)
for k = cfg[2]
lines!(ax,UnfoldSim.PuRF(;(cfg[1]=>k,)...),label=string(k))
for k in cfg[2]
lines!(ax, UnfoldSim.PuRF(; (cfg[1] => k,)...), label = string(k))
end
axislegend(string(cfg[1]);merge=true,)

axislegend(string(cfg[1]); merge = true)
end
f
f
2 changes: 1 addition & 1 deletion docs/literate/reference/noisetypes.jl
Original file line number Diff line number Diff line change
Expand Up @@ -34,4 +34,4 @@ 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.
# 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.
2 changes: 1 addition & 1 deletion docs/literate/reference/onsettypes.jl
Original file line number Diff line number Diff line change
Expand Up @@ -276,4 +276,4 @@ end
# - 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

# [^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#
# [^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#
2 changes: 0 additions & 2 deletions docs/literate/reference/overview.jl
Original file line number Diff line number Diff line change
Expand Up @@ -18,5 +18,3 @@ subtypes(AbstractOnset)
# ## Noise
# Choose the noise you need!
subtypes(AbstractNoise)


28 changes: 14 additions & 14 deletions docs/literate/tutorials/quickstart.jl
Original file line number Diff line number Diff line change
Expand Up @@ -8,39 +8,39 @@ using CairoMakie

# ## "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"])
) |> x->RepeatDesign(x,10);
design =
SingleSubjectDesign(; conditions = Dict(:condA => ["levelA", "levelB"])) |>
x -> RepeatDesign(x, 10);

# #### Component / Signal
# Define a simple component and ground truth simulation formula. Akin to ERP components, we call one simulation signal a component.
#
# !!! 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],
formula = @formula(0~1+condA),
β = [1,0.5]
);
basis = [0, 0, 0, 0.5, 1, 1, 0.5, 0, 0],
formula = @formula(0 ~ 1 + condA),
β = [1, 0.5],
);

# #### Onsets and Noise
# We will start with a uniform (but overlapping, `offset` < `length(signal.basis)`) onset-distribution
onset = UniformOnset(;width=20,offset=4);
onset = UniformOnset(; width = 20, offset = 4);

# And we will use some noise
noise = PinkNoise(;noiselevel=0.2);
noise = PinkNoise(; noiselevel = 0.2);

# ## Combine & Generate
# We will put it all together in one `Simulation` type
simulation = Simulation(design, signal, onset, noise);
simulation = Simulation(design, signal, onset, noise);

# finally, we will simulate some data
data,events = simulate(MersenneTwister(1),simulation);
data, events = simulate(MersenneTwister(1), simulation);
# Data is a `n-sample` Vector (but could be a Matrix for e.g. `MultiSubjectDesign`).

# events is a DataFrame that contains a column `latency` with the onsets of events.

# ## Plot them!
lines(data;color="black")
vlines!(events.latency;color=["orange","teal"][1 .+ (events.condA.=="levelB")])
current_figure()
lines(data; color = "black")
vlines!(events.latency; color = ["orange", "teal"][1 .+ (events.condA.=="levelB")])
current_figure()
2 changes: 1 addition & 1 deletion docs/literate/tutorials/simulateERP.jl
Original file line number Diff line number Diff line change
Expand Up @@ -73,4 +73,4 @@ f = plot_erp(
## Workaround to separate legend and colorbar (will be fixed in a future UnfoldMakie version)
legend = f.content[2]
f[:, 1] = legend
current_figure()
current_figure()
Loading

0 comments on commit 21e46a4

Please sign in to comment.