Addressing the reviews for JOSS submission

.gitignore

@@ -4,4 +4,4 @@ Manifest.toml
 *.css
 *style.jl
 docs/src/examples/*.md
-paper/wip
+paper/data/*.csv

ext/TransitionVisualizations.jl

@@ -51,8 +51,6 @@ function TransitionsInTimeseries.plot_indicator_changes(res::SegmentedWindowResu
     fig, n = init_rowwise_visualisation(res, colors, indicator_names,
         chametric_names, accent_linewidth)
     for k in eachindex(res.t_indicator)
-        Makie.lines!(contents(fig[1, 1])[1], res.t[k], res.x[k], color = colors[1],
-            linewidth = accent_linewidth)
         lineplot_metrics!(fig, n, res.config, res.t_indicator[k], res.x_indicator[k],
             res.t_change[k], res.x_change[k, :], colors, accent_linewidth)
     end
@@ -65,17 +63,17 @@ function init_rowwise_visualisation(res, colors, indicator_names,
     chametric_names, accent_linewidth)
 
     fig = Makie.Figure(size = (700, 450), fontsize = 12)
-    rlabels = vcat([""], indicator_names)
-    llabels = vcat(["Input"], chametric_names)
+    llabels = vcat(["Input"], indicator_names)
+    rlabels = vcat([""], chametric_names)
     n = length(llabels)
 
     # rowaspect = 5
-    raxs = [Makie.Axis(fig[i, 1], ylabel = rlabels[i],
-        xticklabelsvisible = false, yaxisposition = :right, ygridvisible = false,
-        ylabelcolor = colors[2], yticklabelcolor = colors[2]) for i in 1:n]
     laxs = [Makie.Axis(fig[i, 1], ylabel = llabels[i],
         xticklabelsvisible = false, ygridvisible = false,
         ylabelcolor = colors[1], yticklabelcolor = colors[1]) for i in 1:n]
+    raxs = [Makie.Axis(fig[i, 1], ylabel = rlabels[i],
+        xticklabelsvisible = false, yaxisposition = :right, ygridvisible = false,
+        ylabelcolor = colors[2], yticklabelcolor = colors[2]) for i in 1:n]
     Makie.linkxaxes!(laxs..., raxs...)
 
     hidedecorations!(raxs[1])
@@ -90,7 +88,7 @@ function init_rowwise_visualisation(res, colors, indicator_names,
     labels = ["original signal", "surro signals"]
     width = 0.5
     if length(res.t_indicator[1]) > 1
-        elements = vcat(elements, [MarkerElement(marker = :circle, color = colors[1],
+        elements = vcat(elements, [MarkerElement(marker = :circle, color = colors[2],
             strokecolor = :transparent, markersize = ms) for ms in [10, 5]])
         labels = vcat(labels, ["original change metric", "surro change metrics"])
         width = 1
@@ -106,14 +104,14 @@ function lineplot_metrics!(fig, n, config, t_ind, x_ind, t_cha, x_cha,
     for i in 2:n
         j = i-1
         if !isnothing(config.indicators)
-            Makie.lines!(contents(fig[i, 1])[2], t_ind, x_ind[:, j], color = colors[2],
+            Makie.lines!(contents(fig[i, 1])[1], t_ind, x_ind[:, j], color = colors[1],
                 linewidth = accent_linewidth)
         end
         if length(t_cha) > 1
-            lines!(contents(fig[i, 1])[1], t_cha, x_cha[:, j], color = colors[1],
+            lines!(contents(fig[i, 1])[2], t_cha, x_cha[:, j], color = colors[2],
                 linewidth = accent_linewidth)
         else
-            Makie.scatter!(contents(fig[i, 1])[1], t_cha, x_cha[j], color = colors[1],
+            Makie.scatter!(contents(fig[i, 1])[2], t_cha, x_cha[j], color = colors[2],
                 markersize = 10)
         end
     end
@@ -166,26 +164,28 @@ function lines_over_surro!(fig, nsurro, t, t_ind, t_cha, x, signif, config,
     for ns in 1:nsurro
         s = TimeseriesSurrogates.surrogate(x, signif.surrogate)
         Makie.lines!(contents(fig[1, 1])[1], t, s; color = (colors[1], 2/nsurro), linewidth = 1)
-        for (i, cha) in enumerate(config.change_metrics)
 
+        for (i, cha) in enumerate(config.change_metrics)
             if isnothing(config.indicators)
                 p = s
             else
                 p = windowmap(config.indicators[i], s; width = config.width_ind,
                     stride = config.stride_ind)
-                Makie.lines!(contents(fig[i+1, 1])[2], t_ind, p; color = (colors[2], 2/nsurro),
+                Makie.lines!(contents(fig[i+1, 1])[1], t_ind, p; color = (colors[1], 2/nsurro),
                     linewidth = 1)
             end
+
             if length(t_cha) > 1
                 q = windowmap(cha, p; width = config.width_cha, stride = config.stride_cha)
-                Makie.lines!(contents(fig[i+1, 1])[1], t_cha, q; color = (colors[1], 2/nsurro),
+                Makie.lines!(contents(fig[i+1, 1])[2], t_cha, q; color = (colors[2], 2/nsurro),
                     linewidth = 1)
             else
                 cha = precompute(cha, eachindex(p))
                 q = windowmap(cha, p; width = length(p), stride = 1)
-                Makie.scatter!(contents(fig[i+1, 1])[1], t_cha, q[1], color = (colors[1], 2/nsurro),
+                Makie.scatter!(contents(fig[i+1, 1])[2], t_cha, q[1], color = (colors[2], 2/nsurro),
                     markersize = 5)
             end
+
             if !isnothing(flags) && ns == 1 && length(t_cha) > 1
                 Makie.vlines!(contents(fig[i+1, 1])[1], t_cha[flags[:, i]];
                     color = (:black, 0.1))

paper/code/Project.toml

@@ -1,4 +1,4 @@
 [deps]
+BenchmarkTools = "6e4b80f9-dd63-53aa-95a3-0cdb28fa8baf"
 CairoMakie = "13f3f980-e62b-5c42-98c6-ff1f3baf88f0"
 DelimitedFiles = "8bb1440f-4735-579b-a4ab-409b98df4dab"
-TransitionsInTimeseries = "5f5b98ec-1183-43e0-887a-12fdc55c52f7"

paper/code/ewstools-tuto-1.py

@@ -3,58 +3,74 @@
 import matplotlib.pyplot as plt
 import ewstools
 from ewstools.models import simulate_ricker
+from importlib.metadata import version
 
-# Set seed for reproducibility
-np.random.seed(0)
+if version('ewstools') != '2.1.0':
+    raise ValueError('Please install version 2.1.0 of ewstools')
 
-# Initialize time series and spectrum computation
-series = simulate_ricker(tmax=1000, F=[0,2.7])
-ts = ewstools.TimeSeries(data=series, transition=860)
-ts.detrend(method='Lowess', span=0.2)
-ts.state[['state','smoothing']].plot()
-ts.compute_spectrum(rolling_window=0.5, ham_length=40)
+def ewstools_setup():
+    # Initialize time series and spectrum computation
+    series = simulate_ricker(tmax=1000, F=[0,2.7])
+    ts = ewstools.TimeSeries(data=series, transition=860)
+    ts.detrend(method='Lowess', span=0.2)
+    ts.state[['state','smoothing']].plot()
+    ts.compute_spectrum(rolling_window=0.5, ham_length=40)
+    return ts
 
-# Initialize parameters for timing functions
-rw = 0.5
-n = 100
-t_elapsed = np.zeros(10)
+np.random.seed(0)           # Set random seed for reproducibility
+rw = 0.5                    # Rolling window for variance computation
+ts = ewstools_setup()       # Setup time series
+n = 100                     # Number of iterations for timing
+t_elapsed = np.zeros(n)     # Array to store elapsed times
+t_minruntime = np.zeros(8) # Array to store minimum runtime
 
-# Time functions (in a not very elegant way)
-t0 = time.time()
 for i in range(n):
+    t0 = time.time()
     ts.compute_var(rolling_window=rw)
-t_elapsed[0] = time.time() - t0
+    t_elapsed[i] = time.time() - t0
+t_minruntime[0] = min(t_elapsed)
 
-t0 = time.time()
 for i in range(n):
+    t0 = time.time()
     ts.compute_cv()
-t_elapsed[1] = time.time() - t0
+    t_elapsed[i] = time.time() - t0
+t_minruntime[1] = min(t_elapsed)
 
-t0 = time.time()
 for i in range(n):
+    t0 = time.time()
     ts.compute_skew(rolling_window=rw)
-t_elapsed[2] = time.time() - t0
+    t_elapsed[i] = time.time() - t0
+t_minruntime[2] = min(t_elapsed)
 
-t0 = time.time()
 for i in range(n):
+    t0 = time.time()
     ts.compute_kurt()
-t_elapsed[3] = time.time() - t0
+    t_elapsed[i] = time.time() - t0
+t_minruntime[3] = min(t_elapsed)
 
-t0 = time.time()
 for i in range(n):
+    t0 = time.time()
     ts.compute_auto(rolling_window=rw, lag=1)
-t_elapsed[4] = time.time() - t0
+    t_elapsed[i] = time.time() - t0
+t_minruntime[4] = min(t_elapsed)
 
-t0 = time.time()
 for i in range(n):
+    t0 = time.time()
     ts.compute_smax()
-t_elapsed[5] = time.time() - t0
+    t_elapsed[i] = time.time() - t0
+t_minruntime[5] = min(t_elapsed)
 
-t0 = time.time()
 for i in range(n):
+    t0 = time.time()
     ts.compute_ktau()
-t_elapsed[6] = time.time() - t0
+    t_elapsed[i] = time.time() - t0
+t_minruntime[6] = min(t_elapsed)
 
-t0 = time.time()
-surro = ewstools.core.block_bootstrap(ts.state.residuals, n, bs_type='Stationary', block_size=10)
-t_elapsed[7] = time.time() - t0
+for i in range(n):
+    t0 = time.time()
+    ewstools.core.block_bootstrap(ts.state.residuals, 1, bs_type='Stationary', block_size=10)
+    t_elapsed[i] = time.time() - t0
+t_minruntime[7] = min(t_elapsed)
+
+np.savetxt('../data/ewstools_ricker.csv', ts.state[['residuals']].values)
+np.savetxt('../data/ewstools_perfo.csv', t_minruntime, delimiter=',')

paper/code/figure1.jl

@@ -1,6 +1,8 @@
 using TransitionsInTimeseries, DelimitedFiles, CairoMakie, Random
 
-x = readdlm("ewstools-tuto-1.csv", ',')[:, end]
+# Run ewstools-tuto-1.csv first to generate ewstools_ricker.csv
+x = readdlm("../data/ewstools_ricker.csv", ',')[:, end]
+x = x[isnan.(x) .== 0]
 t = eachindex(x)
 
 # Choose the indicators and how to measure their change over time
@@ -12,5 +14,4 @@ results = estimate_changes(config, x, t)
 signif = SurrogatesSignificance(n = 1000, tail = [:right, :right], rng = Xoshiro(1995))
 flags = significant_transitions(results, signif)
 fig = plot_changes_significance(results, signif)
-ylims!(contents(fig[2, 1])[2], (0.037, 0.045))
 save("../figures/figure1.png", fig)

paper/code/figure2.jl

@@ -1,69 +1,60 @@
-using TransitionsInTimeseries, DelimitedFiles, CairoMakie, Random
+using CairoMakie
+using DelimitedFiles
+using Random
+using TransitionsInTimeseries
 
 coefficient_of_variation(x) = std(x) / mean(x)
 
 function main()
-    x = readdlm("ewstools-tuto-1.csv", ',')[:, end]
+    # Run ewstools-tuto-1.csv first to generate ewstools_ricker.csv and ewstools_perfo.csv
+    x = readdlm("../data/ewstools_ricker.csv", ',')[:, end]
+    x = x[isnan.(x) .== 0]
     t = eachindex(x)
 
     # Choose the indicators and how to measure their change over time
     indicators = (var, coefficient_of_variation, skewness, kurtosis,
         ar1_whitenoise, LowfreqPowerSpectrum())
     stride = [1, 1, 1, 1, 1, 40]
-    n, m = 100, length(indicators)
+    m = length(indicators)
     t_elapsed = zeros(m+2)
 
-
     for (i, ind) in enumerate(indicators)
         # Build configuration with adequate parameters of the sliding window
         config = SegmentedWindowConfig((ind, ind), (nothing, nothing), [t[1]], [t[end]];
             width_ind = length(x) ÷ 2, stride_ind = stride[i], whichtime = last,
             min_width_cha = 1)
 
-        t0 = time()
-        for i in 1:n
-            # Compute the metrics over sliding windows and their significance
-            results = estimate_changes(config, x, t)
-        end
-        t_elapsed[i] = (time() - t0) / 2
+        t_elapsed[i] = @belapsed estimate_changes($config, $x, $t)
     end
 
     config = SegmentedWindowConfig((nothing, nothing), (kendalltau, kendalltau), [t[1]], [t[end]];
         width_ind = length(x) ÷ 2, stride_ind = 1, whichtime = last, min_width_cha = 1)
-    t0 = time()
-    for i in 1:n
-        results = estimate_changes(config, x, t)
-    end
-    t_elapsed[m+1] = (time() - t0) / 2
+    t_elapsed[m+1] = @belapsed estimate_changes($config, $x, $t)
 
     sgen = surrogenerator(x, BlockShuffle(), Xoshiro(1995))
-    t0 = time()
-    for i in 1:n
-        s = sgen()
-    end
-    t_elapsed[m+2] = time() - t0
+    t_elapsed[m+2] = @belapsed $sgen()
 
     return t_elapsed
 end
 
-t_tt = main()
-t_et = [0.03840542, 0.05554581, 0.03895116, 0.04029274, 7.96556187,
-    2.73067856, 0.39529872, 0.02751493]
-
-# [0.04681492, 8.13679838, 0.04035759, 0.09219241]
-inds = eachindex(t_et)
+t_transitionsintimeries = main()
+t_ewstools = readdlm("../data/ewstools_perfo.csv", ',')[:, end]
+inds = eachindex(t_ewstools)
 w = 0.4
 
-fig, ax = barplot(inds .- 0.5*w, t_et, label = L"ewstools $\,$", width = w,
-    fillto = 1e-5)
-barplot!(ax, inds .+ 0.5*w, t_tt, label = L"TransitionsInTimeseries.jl $\,$",
-    width = w, fillto = 1e-5)
+fig = Figure(resolution = (600, 400))
+ax = Axis(fig[1, 1])
+ylims!(ax, 1e-7, 0.1)
+barplot!(ax, inds .- 0.5*w, t_ewstools, label = L"ewstools $\,$", width = w,
+    fillto = 1e-8)
+barplot!(ax, inds .+ 0.5*w, t_transitionsintimeries, label = L"TransitionsInTimeseries.jl $\,$",
+    width = w, fillto = 1e-8)
 ax.yscale = log10
 ax.xticks = (1:8, [L"Variance $\,$", L"Coeff. of variation $\,$", L"Skewness $\,$",
     L"Kurtosis $\,$", L"Lag-1 autocorr. $\,$", L"Spectral $\,$",
     L"Kendall $\tau$ corr. coeff.", L"Block bootstrap $\,$"])
-ax.ylabel = L"Run time (s) of 100 computations on Ricker model data $\,$"
-ax.yticks = (10.0 .^ (-5:1), [L"10^{%$e}" for e in -5:1])
+ax.ylabel = L"Run time (s) on Ricker model data $\,$"
+ax.yticks = (10.0 .^ (-8:1), [L"10^{%$e}" for e in -8:1])
 ax.xgridvisible = false
 ax.ygridvisible = false
 ax.xticklabelrotation = π / 4