Clean up toolbox

Project.toml

@@ -4,7 +4,6 @@ authors = ["Till Prölß", " René Skukies", " Benedikt Ehinger"]
 version = "0.1.0"
 
 [deps]
-CairoMakie = "13f3f980-e62b-5c42-98c6-ff1f3baf88f0"
 DSP = "717857b8-e6f2-59f4-9121-6e50c889abd2"
 DataFrames = "a93c6f00-e57d-5684-b7b6-d8193f3e46c0"
 DataFramesMeta = "1313f7d8-7da2-5740-9ea0-a2ca25f37964"
@@ -16,3 +15,16 @@ Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
 StableRNGs = "860ef19b-820b-49d6-a774-d7a799459cd3"
 Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
 Unfold = "181c99d8-e21b-4ff3-b70b-c233eddec679"
+
+[compat]
+DSP = "0.8"
+DataFrames = "1"
+DataFramesMeta = "0.15"
+Distributions = "0.25"
+FFTW = "1"
+Parameters = "0.12"
+Peaks = "0.5, 0.6"
+StableRNGs = "1"
+Statistics = "1"
+Unfold = "0.8"
+julia = "1.11, 1.12"

README.md

@@ -53,6 +53,7 @@ cfg = RideConfig(
     r_range = [0, 0.4],
     c_range = [-0.4, 0.4],
     #the range in which the initial peak estimation for the C component is performed
+    formulas = [@formula(0 ~ 1), @formula(0 ~ 1), @formula(0 ~ 1)] # formulas used for S, R, and C component (in this order!); if not specified these will default to `@formula(0~1)`
     c_estimation_range = [0, 0.9],
     #the range for one epoch
     epoch_range = [-0.1, 1]

dev/Project.toml

@@ -18,6 +18,7 @@ StableRNGs = "860ef19b-820b-49d6-a774-d7a799459cd3"
 Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 Unfold = "181c99d8-e21b-4ff3-b70b-c233eddec679"
+UnfoldBIDS = "b54e767b-1ebd-4480-ac2a-f5f4d8853074"
 UnfoldSim = "ed8ae6d2-84d3-44c6-ab46-0baf21700804"
 UnicodePlots = "b8865327-cd53-5732-bb35-84acbb429228"
 

dev/runner.jl

@@ -12,7 +12,7 @@ using BenchmarkTools
 #simulate data
 begin
     sim_inputs = simulation_inputs()
-    sim_inputs.noise = PinkNoise(; noiselevel = 1)
+    sim_inputs.noise = PinkNoise(; noiselevel = 0)
     #sim_inputs.c_width = 30
     #sim_inputs.c_offset = 10
     #sim_inputs.r_width = 60
@@ -47,15 +47,15 @@ begin
         save_interim_results = true,
     )
 
-    save_to_hdf5_ride_format(
+    #=save_to_hdf5_ride_format(
         "./dev/data/simulated_data.h5",
         data,
         evts,
         cfg.epoch_range,
         'S',
         'R',
         cfg.sfreq,
-    )
+    )=#
 
     #remove the C events from the evts table, these will be estimated by the ride algorithm
     evts_without_c = @subset(evts, :event .!= 'C')
@@ -120,7 +120,7 @@ begin
     Legend(f[1, 2], [s, r, c], ["S", "R", "C"])
     Label(f[0, :], text = "Expected Results")
     display(f)
-    save("actual_erps.png", f)
+    #save("actual_erps.png", f)
 end
 
 if true == false

docs/src/11-running_ride.md

@@ -27,6 +27,7 @@ begin
         r_range = [0, 0.4],
         c_range = [-0.4, 0.4],
         #the range in which the initial peak estimation for the C component is performed
+        formulas = [@formula(0 ~ 1), @formula(0 ~ 1), @formula(0 ~ 1)] # formulas used for S, R, and C component; if not specified these will default to `@formula(0~1)`
         c_estimation_range = [0, 0.9],
         #the range for one epoch
         epoch_range = [-0.1, 1]

src/ride/ride_classic_algorithm.jl

@@ -28,31 +28,10 @@ function ride_algorithm(
 
     ## data_preparation
     data_reshaped = reshape(data, (1, :))
-    evts_s = @subset(evts, :event .== 'S')
-    evts_r = @subset(evts, :event .== 'R')
-    interim_results = Vector{RideResults}()
-
-    #epoch data with the cfg.epoch_range to see how many epochs we have
-    #cut evts to match the determined number of epochs
-    #the resulting data_epoched is also used for the c latency estimation
-    data_epoched, data_epoched_times = Unfold.epoch(
-        data = data_reshaped,
-        tbl = evts_s,
-        τ = cfg.epoch_range,
-        sfreq = cfg.sfreq,
-    )
-    n, data_epoched = Unfold.drop_missing_epochs(evts_s, data_epoched)
-    number_epochs = size(data_epoched, 3)
+    data_epoched, evts_s, evts_r, evts, number_epochs =
+        prepare_epoch_info(data_reshaped, evts, cfg)
     raw_erp = mean(data_epoched, dims = 3)[1, :, 1]
-    #@assert size(evts) == (number_epochs * 2) "Size of evts is $(size(evts)) but should be $(number_epochs * 2)"
-    evts_s = evts_s[1:number_epochs, :]
-    evts_r = evts_r[1:number_epochs, :]
-
-    #reduce evts to the number of epochs
-    while size(evts, 1) > number_epochs * 2
-        deleteat!(evts, size(evts, 1))
-    end
-    @assert size(evts, 1) == number_epochs * 2 "Size of evts is $(size(evts,1)) but should be $(number_epochs*2)"
+    interim_results = Vector{RideResults}()
     ##
 
     ## initial C latency estimation
@@ -134,8 +113,8 @@ function ride_algorithm(
             [(evts_s, s_erp, cfg.s_range), (evts_r, r_erp, cfg.r_range)],
             cfg.sfreq,
         )
-        if cfg.filtering
-            data_subtracted_s_and_r = dspfilter(data_subtracted_s_and_r[1, :], 5, 20)
+        if cfg.filtering # TODO: Check if this is correct; also check for filter artefacts
+            data_subtracted_s_and_r = dspfilter(data_subtracted_s_and_r[1, :], 5, cfg.sfreq)
         end
         data_epoched_subtracted_s_and_r, n = Unfold.epoch(
             data = data_subtracted_s_and_r,

src/ride/ride_classic_methods.jl

@@ -2,11 +2,12 @@ function subtract_to_data(data, others_evts_erp_tuples, sfreq)
     data_subtracted = copy(data)
     for (evts, erp, range) in others_evts_erp_tuples
         for i in evts.latency
-            sub_range =
-                i+round(Int, range[1] * sfreq):i+round(Int, (range[1] * sfreq))+size(
-                    erp[1, :, 1],
-                )[1]-1
-            if sub_range[end] > length(data_subtracted)
+            sub_start = round(Int, i + range[1] * sfreq)
+            sub_end = round(Int, i + (range[1] * sfreq) + size(erp[1, :, 1], 1) - 1)
+            sub_range = sub_start:sub_end
+
+            # Check both start and end bounds
+            if sub_start < 1 || sub_end > length(data_subtracted)
                 continue
             end
             data_subtracted[1, sub_range] -= erp[1, :, 1]
@@ -23,14 +24,13 @@ function subtract_to_data_epoched(
     sfreq,
 )
     data_subtracted = subtract_to_data(data, others_evts_erp_tuples, sfreq)
-    data_epoched_subtracted, n = Unfold.epoch(
+    data_epoched_subtracted, _ = Unfold.epoch(
         data = data_subtracted,
         tbl = target_evts,
         τ = target_range,
         sfreq = sfreq,
     )
-    n, data_epoched_subtracted =
+    _, data_epoched_subtracted =
         Unfold.drop_missing_epochs(target_evts, data_epoched_subtracted)
-    #new_erp = median(data_epoched_subtracted, dims = 3)
     return data_epoched_subtracted
-end
+end