Skip to content

Commit

Permalink
Merge pull request #166 from unfoldtoolbox/formatEverything
Browse files Browse the repository at this point in the history
JuliaFormatter
  • Loading branch information
behinger authored Feb 2, 2024
2 parents dc0990d + bbe51fa commit 8737bb7
Show file tree
Hide file tree
Showing 54 changed files with 1,889 additions and 1,517 deletions.
29 changes: 17 additions & 12 deletions docs/literate/HowTo/effects.jl
Original file line number Diff line number Diff line change
Expand Up @@ -14,14 +14,14 @@ using UnfoldMakie

# # Setup things
# Let's generate some data and fit a model of a 2-level categorical and a continuous predictor with interaction.
data,evts = UnfoldSim.predef_eeg(;noiselevel=8)
data, evts = UnfoldSim.predef_eeg(; noiselevel = 8)

basisfunction = firbasis= (-0.1, 0.5), sfreq = 100, name = "basisA")


f = @formula 0 ~ 1+condition+continuous # 1
f = @formula 0 ~ 1 + condition + continuous # 1

m = fit(UnfoldModel, Dict(Any=>(f,basisfunction)), evts, data,eventcolumn="type")
m = fit(UnfoldModel, Dict(Any => (f, basisfunction)), evts, data, eventcolumn = "type")

# Plot the results
plot_erp(coeftable(m))
Expand All @@ -30,23 +30,28 @@ plot_erp(coeftable(m))
# ### Effects
# A convenience function is `effects`. It allows to specify effects on specific levels, while setting non-specified ones to a typical value (usually the mean)

eff = effects(Dict(:condition => ["car","face"]),m)
plot_erp(eff;mapping=(;color=:condition,))
eff = effects(Dict(:condition => ["car", "face"]), m)
plot_erp(eff; mapping = (; color = :condition,))

# We can also generate continuous predictions
eff = effects(Dict(:continuous => -5:0.5:5),m)
plot_erp(eff;mapping=(;color=:continuous,group=:continuous=>nonnumeric),categorical_color=false,categorical_group=false)
eff = effects(Dict(:continuous => -5:0.5:5), m)
plot_erp(
eff;
mapping = (; color = :continuous, group = :continuous => nonnumeric),
categorical_color = false,
categorical_group = false,
)

# or split it up by condition
eff = effects(Dict(:condition=>["car","face"],:continuous => -5:2:5),m)
plot_erp(eff;mapping=(;color=:condition,col=:continuous=>nonnumeric))
eff = effects(Dict(:condition => ["car", "face"], :continuous => -5:2:5), m)
plot_erp(eff; mapping = (; color = :condition, col = :continuous => nonnumeric))

# ## What is typical anyway?
# The user can specify the typical function applied to the covariates/factors that are marginalized over. This offers even greater flexibility.
# Note that this is rarely necessary, in most applications the mean will be assumed.
eff_max = effects(Dict(:condition=>["car","face"]),m;typical=maximum) # typical continuous value fixed to 1
eff_max = effects(Dict(:condition => ["car", "face"]), m; typical = maximum) # typical continuous value fixed to 1
eff_max.typical .= :maximum
eff = effects(Dict(:condition=>["car","face"]),m)
eff = effects(Dict(:condition => ["car", "face"]), m)
eff.typical .= :mean # mean is the default

plot_erp(vcat(eff,eff_max);mapping=(;color=:condition,col=:typical))
plot_erp(vcat(eff, eff_max); mapping = (; color = :condition, col = :typical))
2 changes: 1 addition & 1 deletion docs/literate/HowTo/juliacall_unfold.jl
Original file line number Diff line number Diff line change
Expand Up @@ -37,4 +37,4 @@
# ```

# ## Example: Unfold model fitting from Python
# In this [notebook](https://github.com/unfoldtoolbox/Unfold.jl/blob/main/docs/src/HowTo/juliacall_unfold.ipynb), you can find a more detailed example of how to use Unfold from Python to load data, fit an Unfold model and visualise the results in Python.
# In this [notebook](https://github.com/unfoldtoolbox/Unfold.jl/blob/main/docs/src/HowTo/juliacall_unfold.ipynb), you can find a more detailed example of how to use Unfold from Python to load data, fit an Unfold model and visualise the results in Python.
1 change: 1 addition & 0 deletions docs/literate/HowTo/timesplines.jl
Original file line number Diff line number Diff line change
@@ -0,0 +1 @@

16 changes: 8 additions & 8 deletions docs/literate/HowTo/unfold_io.jl
Original file line number Diff line number Diff line change
Expand Up @@ -10,14 +10,14 @@

# ### Simulate some example data using UnfoldSim.jl
using UnfoldSim
data, events = UnfoldSim.predef_eeg(; n_repeats=10)
first(events,5)
data, events = UnfoldSim.predef_eeg(; n_repeats = 10)
first(events, 5)

# ### Fit an Unfold model
using Unfold
basisfunction = firbasis= (-0.5,1.0), sfreq = 100, name = "stimulus")
f = @formula 0 ~ 1 + condition + continuous
bfDict = Dict(Any => (f,basisfunction))
basisfunction = firbasis= (-0.5, 1.0), sfreq = 100, name = "stimulus")
f = @formula 0 ~ 1 + condition + continuous
bfDict = Dict(Any => (f, basisfunction))
m = fit(UnfoldModel, bfDict, events, data);

# ```@raw html
Expand All @@ -28,9 +28,9 @@ m = fit(UnfoldModel, bfDict, events, data);

# The following code saves the model in a compressed .jld2 file. The default option of the `save` function is `compress=false`.
# For memory efficiency the designmatrix is set to missing. If needed, it can be reconstructed when loading the model.
save_path = mktempdir(; cleanup=false) # create a temporary directory for the example
save(joinpath(save_path, "m_compressed.jld2"), m; compress=true);
save_path = mktempdir(; cleanup = false) # create a temporary directory for the example
save(joinpath(save_path, "m_compressed.jld2"), m; compress = true);

# The `load` function allows to retrieve the model again.
# By default, the designmatrix is reconstructed. If it is not needed set `generate_Xs=false`` which improves time-efficiency.
m_loaded = load(joinpath(save_path, "m_compressed.jld2"), UnfoldModel, generate_Xs=true);
m_loaded = load(joinpath(save_path, "m_compressed.jld2"), UnfoldModel, generate_Xs = true);
70 changes: 35 additions & 35 deletions docs/literate/explanations/nonlinear_effects.jl
Original file line number Diff line number Diff line change
Expand Up @@ -2,7 +2,7 @@



using BSplineKit,Unfold
using BSplineKit, Unfold
using CairoMakie
using DataFrames
using Random
Expand All @@ -14,42 +14,42 @@ using Colors
# We start with generating
rng = MersenneTwister(2) # make repeatable
n = 20 # datapoints
evts = DataFrame(:x=>rand(rng,n))
signal = -(3*(evts.x .-0.5)).^2 .+ 0.5 .* rand(rng,n)
evts = DataFrame(:x => rand(rng, n))
signal = -(3 * (evts.x .- 0.5)) .^ 2 .+ 0.5 .* rand(rng, n)

plot(evts.x,signal)
plot(evts.x, signal)
#
# Looks perfectly non-linear - great!
#
# # Compare linear & non-linear fit
# First we have to reshape it to fit to Unfold format, a 3d array, 1 channel x 1 timepoint x 20 datapoints
signal = reshape(signal,length(signal),1,1)
signal = permutedims(signal,[3,2,1])
signal = reshape(signal, length(signal), 1, 1)
signal = permutedims(signal, [3, 2, 1])
size(signal)

# Next we define three different models. **linear** **4 splines** and **10 splines**
# note the different formulas, one `x` the other `spl(x,4)`
design_linear = Dict(Any=> (@formula(0~1+x),[0]))
design_spl3 = Dict(Any=> (@formula(0~1+spl(x,4)),[0]))
design_spl10 = Dict(Any=> (@formula(0~1+spl(x,10)),[0])) #hide
design_linear = Dict(Any => (@formula(0 ~ 1 + x), [0]))
design_spl3 = Dict(Any => (@formula(0 ~ 1 + spl(x, 4)), [0]))
design_spl10 = Dict(Any => (@formula(0 ~ 1 + spl(x, 10)), [0])) #hide

# and we fit the parameters
uf_linear = fit(UnfoldModel,design_linear,evts,signal);
uf_spl3 = fit(UnfoldModel,design_spl3,evts,signal);
uf_spl10 = fit(UnfoldModel,design_spl10,evts,signal); #hide
uf_linear = fit(UnfoldModel, design_linear, evts, signal);
uf_spl3 = fit(UnfoldModel, design_spl3, evts, signal);
uf_spl10 = fit(UnfoldModel, design_spl10, evts, signal); #hide

# extract the fitted values using Unfold.effects

p_linear = Unfold.effects(Dict(:x => range(0,stop=1,length=100)),uf_linear);
p_spl3 = Unfold.effects(Dict(:x => range(0,stop=1,length=100)),uf_spl3);
p_spl10 = Unfold.effects(Dict(:x => range(0,stop=1,length=100)),uf_spl10);
first(p_linear,5)
p_linear = Unfold.effects(Dict(:x => range(0, stop = 1, length = 100)), uf_linear);
p_spl3 = Unfold.effects(Dict(:x => range(0, stop = 1, length = 100)), uf_spl3);
p_spl10 = Unfold.effects(Dict(:x => range(0, stop = 1, length = 100)), uf_spl10);
first(p_linear, 5)

# And plot them
pl = plot(evts.x,signal[1,1,:])
lines!(p_linear.x,p_linear.yhat)
lines!(p_spl3.x,p_spl3.yhat)
lines!(p_spl10.x,p_spl10.yhat)
pl = plot(evts.x, signal[1, 1, :])
lines!(p_linear.x, p_linear.yhat)
lines!(p_spl3.x, p_spl3.yhat)
lines!(p_spl10.x, p_spl10.yhat)
pl

# What we can clearly see here, that the linear effect (blue) underfits the data, the `spl(x,10)` overfits it, but the `spl(x,4)` fits it perfectly.
Expand All @@ -72,46 +72,46 @@ typeof(term_spl)
# a special type generated in Unfold.jl
#
# it has a field .fun which is the spline function. We can evaluate it at a point
const splFunction = Base.get_extension(Unfold,:UnfoldBSplineKitExt).splFunction
splFunction([0.2],term_spl)
const splFunction = Base.get_extension(Unfold, :UnfoldBSplineKitExt).splFunction
splFunction([0.2], term_spl)

# each column of this 1 row matrix is a coefficient for our regression model.
lines(disallowmissing(splFunction([0.2],term_spl))[1,:])
lines(disallowmissing(splFunction([0.2], term_spl))[1, :])

# Note: We have to use disallowmissing, because our splines return a `missing` whenever we ask it to return a value outside it's defined range, e.g.:
splFunction([-0.2],term_spl)
splFunction([-0.2], term_spl)

# because it never has seen any data outside and can't extrapolate!

# Okay back to our main issue: Let's plot the whole basis set
basisSet = splFunction(0.:0.01:1,term_spl)
basisSet = disallowmissing(basisSet[.!any(ismissing.(basisSet),dims=2)[:,1],:]) # remove missings
ax = Axis(Figure()[1,1])
[lines!(ax,basisSet[:,k]) for k in 1:size(basisSet,2)]
basisSet = splFunction(0.0:0.01:1, term_spl)
basisSet = disallowmissing(basisSet[.!any(ismissing.(basisSet), dims = 2)[:, 1], :]) # remove missings
ax = Axis(Figure()[1, 1])
[lines!(ax, basisSet[:, k]) for k = 1:size(basisSet, 2)]
current_figure()

# notice how we flipped the plot around, i.e. now on the x-axis we do not plot the coefficients, but the `x`-values
# Now each line is one basis-function of the spline.
#
# Unfold returns us one coefficient per basis-function
β=coef(uf_spl10)[1,1,:]
β = coef(uf_spl10)[1, 1, :]
β = Float64.(disallowmissing(β))

# but because we used an intercept, we have to do some remodelling in the basisSet
X = hcat(ones(size(basisSet,1)), basisSet[:,1:5], basisSet[:,7:end])
X = hcat(ones(size(basisSet, 1)), basisSet[:, 1:5], basisSet[:, 7:end])

# now we can weight the spline by the basisfunction
weighted =.*X')
weighted =.* X')

# Plotting them makes for a nice looking plot!
ax = Axis(Figure()[1,1])
[lines!(weighted[k,:]) for k = 1:10]
ax = Axis(Figure()[1, 1])
[lines!(weighted[k, :]) for k = 1:10]
current_figure()


# and sum them up
lines(sum(weighted,dims=1)[1,:])
plot!(X*β,color="gray") #(same as matrixproduct X*β directly!)
lines(sum(weighted, dims = 1)[1, :])
plot!(X * β, color = "gray") #(same as matrixproduct X*β directly!)
current_figure()

# And this is how you can think about splines
Loading

0 comments on commit 8737bb7

Please sign in to comment.