Integrating code into bestie template

.gitignore

@@ -10,3 +10,7 @@ coverage
 docs/build/
 env
 node_modules
+*.txt
+/*.png
+/*.h5
+*.so

Project.toml

@@ -1,11 +1,18 @@
 name = "UnfoldRIDE"
 uuid = "348b5554-d1f6-52bd-94b1-729601d351bf"
-authors = [
-    "Till Prölß",
-    " René Skukies",
-    " Benedikt Ehinger",
-]
+authors = ["Till Prölß", " René Skukies", " Benedikt Ehinger"]
 version = "0.1.0"
 
-[compat]
-julia = "1.10"
+[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"
+Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f"
+FFTW = "7a1cc6ca-52ef-59f5-83cd-3a7055c09341"
+Parameters = "d96e819e-fc66-5662-9728-84c9c7592b0a"
+Peaks = "18e31ff7-3703-566c-8e60-38913d67486b"
+Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
+StableRNGs = "860ef19b-820b-49d6-a774-d7a799459cd3"
+Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
+Unfold = "181c99d8-e21b-4ff3-b70b-c233eddec679"

README.md

@@ -6,11 +6,62 @@
 [![Lint workflow Status](https://github.com/unfoldtoolbox/UnfoldRIDE.jl/actions/workflows/Lint.yml/badge.svg?branch=main)](https://github.com/unfoldtoolbox/UnfoldRIDE.jl/actions/workflows/Lint.yml?query=branch%3Amain)
 [![Docs workflow Status](https://github.com/unfoldtoolbox/UnfoldRIDE.jl/actions/workflows/Docs.yml/badge.svg?branch=main)](https://github.com/unfoldtoolbox/UnfoldRIDE.jl/actions/workflows/Docs.yml?query=branch%3Amain)
 [![Coverage](https://codecov.io/gh/unfoldtoolbox/UnfoldRIDE.jl/branch/main/graph/badge.svg)](https://codecov.io/gh/unfoldtoolbox/UnfoldRIDE.jl)
-[![DOI](https://zenodo.org/badge/DOI/FIXME)](https://doi.org/FIXME)
+<!--[![DOI](https://zenodo.org/badge/DOI/FIXME)](https://doi.org/FIXME)-->
 [![Contributor Covenant](https://img.shields.io/badge/Contributor%20Covenant-2.1-4baaaa.svg)](CODE_OF_CONDUCT.md)
 [![All Contributors](https://img.shields.io/github/all-contributors/unfoldtoolbox/UnfoldRIDE.jl?labelColor=5e1ec7&color=c0ffee&style=flat-square)](#contributors)
 [![BestieTemplate](https://img.shields.io/endpoint?url=https://raw.githubusercontent.com/JuliaBesties/BestieTemplate.jl/main/docs/src/assets/badge.json)](https://github.com/JuliaBesties/BestieTemplate.jl)
 
+A re-implementation of the [RIDE](https://cns.hkbu.edu.hk/RIDE.htm) algorithm in Julia with an extension to replace the RIDEs iterative decomposition with an [Unfold](https://github.com/unfoldtoolbox/Unfold.jl) deconvolution.
+
+## Install
+
+### Installing Julia
+
+<details>
+<summary>Click to expand</summary>
+
+The recommended way to install julia is [juliaup](https://github.com/JuliaLang/juliaup).
+It allows you to, e.g., easily update Julia at a later point, but also test out alpha/beta versions etc.
+
+TL:DR; If you dont want to read the explicit instructions, just copy the following command
+
+#### Windows
+
+AppStore -> JuliaUp,  or `winget install julia -s msstore` in CMD
+
+#### Mac & Linux
+
+`curl -fsSL https://install.julialang.org | sh` in any shell
+</details>
+
+### Installing UnfoldRIDE
+
+```julia
+using Pkg
+Pkg.add("UnfoldRIDE")
+```
+
+## Quickstart
+
+```Julia
+#config for ride algorithm
+cfg = RideConfig(
+    #sfreq is the sampling frequency of the data
+    sfreq = 100,
+    #ranges for the individual components have to be determined through manual inspection of the data
+    s_range = [-0.1, 0.3],
+    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
+    c_estimation_range = [0, 0.9],
+    #the range for one epoch
+    epoch_range = [-0.1, 1]
+)
+#run the ride algorithm
+resultsClassic = ride_algorithm(ClassicMode, data_noisy, evts_without_c, cfg)
+resultsUnfold = ride_algorithm(UnfoldMode, data_noisy, evts_without_c, cfg)
+```
+
 ## How to Cite
 
 If you use UnfoldRIDE.jl in your work, please cite using the reference given in [CITATION.cff](https://github.com/unfoldtoolbox/UnfoldRIDE.jl/blob/main/CITATION.cff).

dev/Project.toml

@@ -0,0 +1,25 @@
+[deps]
+BenchmarkTools = "6e4b80f9-dd63-53aa-95a3-0cdb28fa8baf"
+CairoMakie = "13f3f980-e62b-5c42-98c6-ff1f3baf88f0"
+DSP = "717857b8-e6f2-59f4-9121-6e50c889abd2"
+DataFrames = "a93c6f00-e57d-5684-b7b6-d8193f3e46c0"
+DataFramesMeta = "1313f7d8-7da2-5740-9ea0-a2ca25f37964"
+Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f"
+FFTW = "7a1cc6ca-52ef-59f5-83cd-3a7055c09341"
+HDF5 = "f67ccb44-e63f-5c2f-98bd-6dc0ccc4ba2f"
+LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
+PackageCompiler = "9b87118b-4619-50d2-8e1e-99f35a4d4d9d"
+Parameters = "d96e819e-fc66-5662-9728-84c9c7592b0a"
+Peaks = "18e31ff7-3703-566c-8e60-38913d67486b"
+Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
+Revise = "295af30f-e4ad-537b-8983-00126c2a3abe"
+SignalAnalysis = "df1fea92-c066-49dd-8b36-eace3378ea47"
+StableRNGs = "860ef19b-820b-49d6-a774-d7a799459cd3"
+Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
+Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
+Unfold = "181c99d8-e21b-4ff3-b70b-c233eddec679"
+UnfoldSim = "ed8ae6d2-84d3-44c6-ab46-0baf21700804"
+UnicodePlots = "b8865327-cd53-5732-bb35-84acbb429228"
+
+[sources]
+UnfoldSim = {rev = "v4.0", url = "https://github.com/unfoldtoolbox/UnfoldSim.jl"}

dev/create_sysimage.jl

@@ -0,0 +1,27 @@
+# This script creates a sysimage for the RIDE project
+# WARNING: This will take a LONG time (1-2 hours)
+# The sysimage massively speeds up julia startup time
+# After the sysimage is created, you need to add:
+# "-JRIDE_Sysimage.so" to your "julia.additionalArgs" in your vscode settings
+# You can use Base.loaded_modules to confirm that the sysimage is loaded correctly
+
+include("./runner.jl")
+# Get all packages from the current project
+all_packages = String[]
+for (uuid, dep) in Pkg.dependencies()
+    dep.is_direct_dep || continue
+    dep.version === nothing && continue
+    push!(all_packages, dep.name)
+end
+
+# Remove unneeded packages
+do_not_include = ["PackageCompiler"]
+package_list = filter(x -> x ∉ do_not_include, all_packages)
+
+using PackageCompiler
+
+PackageCompiler.create_sysimage(
+    package_list;
+    sysimage_path = "RIDE_Sysimage.so",
+    precompile_execution_file = "./dev/runner.jl",
+)

dev/plotting_methods.jl

@@ -0,0 +1,78 @@
+function plot_c_latency_estimation_four_epochs(data_epoched, c_latencies, c_erp)
+    #plotting the estimated c latencies on a couple of epochs
+    f = Figure()
+    c_latency_from_epoch_start = c_latencies .- round(Int, cfg.epoch_range[1] * cfg.sfreq)
+    for a = 1:2, b = 1:2
+        i = (a - 1) * 2 + b
+        Axis(f[a, b], title = "Estimated C latency epoch $i")
+        lines!(f[a, b], data_epoched[1, :, i]; color = "black")
+        y = c_latency_from_epoch_start[i]:(c_latency_from_epoch_start[i]+length(c_erp)-1)
+        lines!(f[a, b], y, c_erp; color = "red")
+    end
+    return f
+end
+
+function plot_first_three_epochs_of_raw_data(data, evts)
+    f = Figure()
+    Axis(f[1, 1], title = "First three epochs")
+    evts_s = @subset(evts, :event .== 'S')
+    evts_r = @subset(evts, :event .== 'R')
+    graph = lines!(first(vec(data), evts_s.latency[4]); color = "black")
+    graph_r = vlines!(first(evts_r.latency, 3), color = "blue")
+    graph_s = vlines!(first(evts_s.latency, 3), color = "red")
+
+    Legend(
+        f[1, 2],
+        [graph, graph_r, graph_s],
+        ["Data", "Reaction Times", "Stimulus Onsets"],
+    )
+    display(f)
+end
+
+function plot_first_three_epochs_of_raw_data(data)
+    f = Figure()
+    Axis(f[1, 1], title = "First 600 data points")
+    graph = lines!(first(vec(data), 600); color = "black")
+    Legend(f[1, 2], [graph], ["Data"])
+    display(f)
+end
+
+function plot_interim_results(data, evts, results, cfg)
+    evts_s = @subset(evts, :event .== 'S')
+    data_epoched =
+        Unfold.epoch(data = data, tbl = evts_s, τ = cfg.epoch_range, sfreq = cfg.sfreq)[1]
+    data_epoched = Unfold.drop_missing_epochs(evts_s, data_epoched)[2]
+    raw_erp = mean(data_epoched, dims = 3)[1, :, 1]
+    for (i, r) in enumerate(vcat(results.interim_results, results))
+        f = plot_c_latency_estimation_four_epochs(
+            data_epoched,
+            r.c_latencies,
+            r.c_erp_unpadded,
+        )
+        Label(f[0, :], text = "Estimated C latency, Iteration $(i-1)", halign = :center)
+        display(f)
+    end
+    for (i, r) in enumerate(vcat(results.interim_results, results))
+        f = plot_simple_results(
+            raw_erp,
+            r.s_erp,
+            r.r_erp,
+            r.c_erp;
+            legend = true,
+            title = "Calculated Components, Iteration $(i-1)",
+        )
+        display(f)
+    end
+end
+
+function plot_simple_results(raw_erp, s, r, c; legend = false, f = Figure(), title = "")
+    ax = Axis(f[1, 1], yticks = -100:100, title = title)
+    raw = lines!(raw_erp; color = "black", linewidth = 3, label = "ERP")
+    s = lines!(s; color = "blue", label = "S")
+    c = lines!(c; color = "red", label = "C")
+    r = lines!(r; color = "green", label = "R")
+    if legend
+        axislegend(ax)
+    end
+    return f
+end

dev/runner.jl

@@ -0,0 +1,130 @@
+using Pkg
+Pkg.activate("./dev")
+
+using Revise
+includet("../src/UnfoldRIDE.jl")
+includet("../test/simulate_test_data.jl")
+includet("./plotting_methods.jl")
+using .UnfoldRIDE
+using CairoMakie
+using StableRNGs
+using BenchmarkTools
+#simulate data
+begin
+    sim_inputs = simulation_inputs()
+    sim_inputs.noise = PinkNoise(; noiselevel = 1)
+    #sim_inputs.c_width = 30
+    #sim_inputs.c_offset = 10
+    #sim_inputs.r_width = 60
+    #sim_inputs.r_offset = 15
+    #sim_inputs.c_beta = -4
+    #sim_inputs.r_beta = 6
+    #sim_inputs.r_offset = 0
+    #sim_inputs.c_beta = -5
+    #sim_inputs.c_offset = 50
+    #sim_inputs.s_offset = 70
+    data, evts, data_clean, evts_clean, data_clean_s, data_clean_r, data_clean_c =
+        simulate_default_plus_clean(sim_inputs)
+    plot_first_three_epochs_of_raw_data(data_clean_s, evts)
+    plot_first_three_epochs_of_raw_data(data, evts)
+end
+
+#run the ride algorithm on the simulated data
+begin
+    #ENV["JULIA_DEBUG"] = "UnfoldRIDE"
+    #config for ride algorithm
+    cfg = RideConfig(
+        sfreq = 100,
+        s_range = [-0.2, 0.4],
+        r_range = [0, 0.8],
+        c_range = [-0.4, 0.4],
+        c_estimation_range = [-0.1, 0.9],
+        epoch_range = [-0.3, 1.6],
+        iteration_limit = 5,
+        heuristic1 = true,
+        heuristic2 = true,
+        heuristic3 = true,
+        save_interim_results = true,
+    )
+
+    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')
+
+    #run the ride algorithm
+    results = ride_algorithm(ClassicMode, data, evts_without_c, cfg)
+    results2 = ride_algorithm(UnfoldMode, data, evts_without_c, cfg)
+    s_erp = results[1].s_erp
+    r_erp = results[1].r_erp
+    c_erp = results[1].c_erp
+    c_latencies = results[1].c_latencies
+
+    plot_interim_results(data, evts, results[1], cfg)
+end
+
+# calculate and plot clean erps from the simulated data
+# these represent the optimal result from the algorithm
+begin
+    evts_clean_s = @subset(evts_clean, :event .== 'S')
+    data_epoched_clean = Unfold.epoch(
+        data = data_clean,
+        tbl = evts_clean_s,
+        τ = cfg.epoch_range,
+        sfreq = cfg.sfreq,
+    )[1]
+    n, data_epoched_clean = Unfold.drop_missing_epochs(evts_clean_s, data_epoched_clean)
+    erp_clean = mean(data_epoched_clean, dims = 3)
+
+    data_epoched_clean = Unfold.epoch(
+        data = data_clean_s,
+        tbl = evts_clean_s,
+        τ = cfg.epoch_range,
+        sfreq = cfg.sfreq,
+    )[1]
+    n, data_epoched_clean = Unfold.drop_missing_epochs(evts_clean_s, data_epoched_clean)
+    erp_clean_s = mean(data_epoched_clean, dims = 3)
+
+    data_epoched_clean = Unfold.epoch(
+        data = data_clean_c,
+        tbl = evts_clean_s,
+        τ = cfg.epoch_range,
+        sfreq = cfg.sfreq,
+    )[1]
+    n, data_epoched_clean = Unfold.drop_missing_epochs(evts_clean_s, data_epoched_clean)
+    erp_clean_c = mean(data_epoched_clean, dims = 3)
+
+    data_epoched_clean = Unfold.epoch(
+        data = data_clean_r,
+        tbl = evts_clean_s,
+        τ = cfg.epoch_range,
+        sfreq = cfg.sfreq,
+    )[1]
+    n, data_epoched_clean = Unfold.drop_missing_epochs(evts_clean_s, data_epoched_clean)
+    erp_clean_r = mean(data_epoched_clean, dims = 3)
+
+    #plot the results
+    f = Figure()
+    Axis(f[1, 1], yticks = -100:100)
+    s = lines!(f[1, 1], erp_clean_s[1, :, 1], color = "blue")
+    c = lines!(f[1, 1], erp_clean_c[1, :, 1], color = "red")
+    r = lines!(f[1, 1], erp_clean_r[1, :, 1], color = "green")
+    Legend(f[1, 2], [s, r, c], ["S", "R", "C"])
+    Label(f[0, :], text = "Expected Results")
+    display(f)
+    save("actual_erps.png", f)
+end
+
+if true == false
+    @benchmark ride_algorithm(UnfoldMode, data, evts_without_c, cfg)
+
+    @benchmark ride_algorithm(ClassicMode, data, evts_without_c, cfg)
+end