From 64bb7af0153b7084202e71bebe83b977bdb8bb43 Mon Sep 17 00:00:00 2001 From: George Datseris Date: Sun, 17 Dec 2023 19:07:39 +0000 Subject: [PATCH] Add KS-statistic example (#61) * add ks initial draft * outline what we have to do * dont exceed 92 chars * add citations to Witold's work * finish KS statistic code * nicer plot * make it wokr with `nothing` * add threshold significance * finish example * add test for nothing (doesn't work!) * enforce change metrics equal in length to indicators * correct handling of number of metrics * update tutorial with change metric number * fix incorrect number of ` in tuytorial * finish code but tests still don't pass... :S * i still don't understand why the nothing test fails. * make indicators nothing instead of tuple of nothing * fix printing not working * add visualizations code since I anyways use it to test stuff --- Project.toml | 7 + docs/src/api.md | 2 + docs/src/examples/do-events.jl | 65 ++++++- docs/src/examples/ks_paleojump.jl | 168 ++++++++++++++++++ docs/src/refs.bib | 33 ++++ docs/src/tutorial.md | 62 ++++--- ext/TITVisualizations.jl | 69 +++++++ src/TransitionsInTimeseries.jl | 8 +- src/analysis/segmented_window.jl | 4 +- src/analysis/sliding_window.jl | 86 +++++---- src/significance/api_significance.jl | 3 + ...ificance.jl => basic_stat_significance.jl} | 36 ++++ src/significance/surrogates_significance.jl | 33 ++-- src/visualizations.jl | 9 + test/full_analysis.jl | 51 +++++- 15 files changed, 544 insertions(+), 92 deletions(-) create mode 100644 docs/src/examples/ks_paleojump.jl create mode 100644 ext/TITVisualizations.jl rename src/significance/{quantile_significance.jl => basic_stat_significance.jl} (71%) create mode 100644 src/visualizations.jl diff --git a/Project.toml b/Project.toml index cb5aea55..47078df6 100644 --- a/Project.toml +++ b/Project.toml @@ -16,6 +16,12 @@ Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2" StatsBase = "2913bbd2-ae8a-5f71-8c99-4fb6c76f3a91" TimeseriesSurrogates = "c804724b-8c18-5caa-8579-6025a0767c70" +[weakdeps] +Makie = "ee78f7c6-11fb-53f2-987a-cfe4a2b5a57a" + +[extensions] +TITVisualizations = "Makie" + [compat] DelimitedFiles = "1" Downloads = "1" @@ -23,6 +29,7 @@ FFTW = "^1.6" HypothesisTests = "0.11" InteractiveUtils = "1" LinearAlgebra = "1" +Makie = "≥ 0.19" Random = "1" Reexport = "1.2" Statistics = "1" diff --git a/docs/src/api.md b/docs/src/api.md index 5ec1d3da..971758ee 100644 --- a/docs/src/api.md +++ b/docs/src/api.md @@ -18,6 +18,8 @@ SegmentedWindowResults significant_transitions TransitionsSignificance SurrogatesSignificance +ThresholdSignificance +SigmaSignificance QuantileSignificance ``` diff --git a/docs/src/examples/do-events.jl b/docs/src/examples/do-events.jl index ff3974a9..6ed2fd93 100644 --- a/docs/src/examples/do-events.jl +++ b/docs/src/examples/do-events.jl @@ -1,16 +1,36 @@ #= # Dansgaard-Oescher events and Critical Slowing Down -The $\delta^{18}O$ timeseries of the North Greenland Ice Core Project ([NGRIP](https://en.wikipedia.org/wiki/North_Greenland_Ice_Core_Project)) are, to this date, the best proxy record for the Dansgaard-Oeschger events ([DO-events](https://en.wikipedia.org/wiki/Dansgaard%E2%80%93Oeschger_event)). DO-events are sudden warming episodes of the North Atlantic, reaching 10 degrees of regional warming within 100 years. They happened quasi-periodically over the last glacial cycle due to transitions between strong and weak states of the Atlantic Meridional Overturning Circulation and might be therefore be the most prominent examples of abrupt transitions in the field of climate science. We here propose to hindcast these events by applying the theory of Critical Slowing Down (CSD) on the NGRIP data, which can be found [here](https://www.iceandclimate.nbi.ku.dk/data/) in its raw format. This analysis has already been done in [boers-early-warning-2018](@cite) and we here try to reproduce Figure 2.d-f. +The $\delta^{18}O$ timeseries of the North Greenland Ice Core Project +([NGRIP](https://en.wikipedia.org/wiki/North_Greenland_Ice_Core_Project)) are, +to this date, the best proxy record for the Dansgaard-Oeschger events +([DO-events](https://en.wikipedia.org/wiki/Dansgaard%E2%80%93Oeschger_event)). +DO-events are sudden warming episodes of the North Atlantic, reaching 10 degrees +of regional warming within 100 years. They happened quasi-periodically over the +last glacial cycle due to transitions between strong and weak states of the Atlantic +Meridional Overturning Circulation and might be therefore be the most prominent +examples of abrupt transitions in the field of climate science. We here propose +to hindcast these events by applying the theory of Critical Slowing Down (CSD) +on the NGRIP data, which can be found [here](https://www.iceandclimate.nbi.ku.dk/data/) +in its raw format. This analysis has already been done in [boers-early-warning-2018](@cite) +and we here try to reproduce Figure 2.d-f. ## Preprocessing NGRIP -Data pre-processing is not part of TransitionsInTimeseries.jl, but a step the user has to do before using the package. To present an example with a complete scientific workflow, we will showcase typical data pre-processing here, that consist of the following steps: +Data pre-processing is not part of TransitionsInTimeseries.jl, +but a step the user has to do before using the package. +To present an example with a complete scientific workflow, +we will showcase typical data pre-processing here, that consist of the following steps: 1. Load the data, reverse and offset it to have time vector = time before 2000 AD. 2. Filter non-unique points in time and sort the data. 3. Regrid the data from uneven to even sampling. -The time and $\delta^{18}O$ vectors resulting from the $i$-th preprocessing step are respectively called $t_i$ and $x_i$. The final step consists in obtaining a residual $r$, i.e. the fluctuations of the system around the attractor, which, within the CSD theory, is assumed to be tracked. Over this example, it will appear that the convenience of TransitionsInTimeseries.jl leads the bulk of the code to be written for plotting and preprocessing. +The time and $\delta^{18}O$ vectors resulting from the $i$-th preprocessing step are +respectively called $t_i$ and $x_i$. The final step consists in obtaining a residual +$r$, i.e. the fluctuations of the system around the attractor, which, within the CSD +theory, is assumed to be tracked. Over this example, it will appear that the +convenience of TransitionsInTimeseries.jl leads the bulk of the code to be written +for plotting and preprocessing. ### Step 1: =# @@ -84,7 +104,9 @@ fcutoff = 0.95 * 0.01 # cutoff ≃ 0.01 yr^-1 as in (Boers 2018) t, x, xtrend, r = chebyshev_filter(t3, x3, fcutoff) #= -Let's now visualize our data in what will become our main figure. For the segmentation of the DO-events, we rely on the tabulated data from [rasmussen-stratigraphic-2014](@cite) (which will soon be available as downloadable): +Let's now visualize our data in what will become our main figure. +For the segmentation of the DO-events, we rely on the tabulated +data from [rasmussen-stratigraphic-2014](@cite) (which will soon be available as downloadable): =# using CairoMakie, Loess @@ -140,7 +162,15 @@ fig #= ## Hindcast on NGRIP data -As one can see... there is not much to see so far. Residuals are impossible to simply eye-ball and we therefore use TransitionsInTimeseries.jl to study the evolution, measured by the ridge-regression slope of the residual's variance and lag-1 autocorrelation (AC1) over time. In many examples of the literature, including [boers-early-warning-2018](@cite), the CSD analysis is performed over segments (sometimes only one) of the timeseries, such that a significance value is obtained for each segment. By using `SegmentedWindowConfig`, dealing with segments can be easily done in TransitionsInTimeseries.jl and is demonstrated here: +As one can see... there is not much to see so far. +Residuals are impossible to simply eye-ball and we therefore use +TransitionsInTimeseries.jl to study the evolution, measured by the ridge-regression +slope of the residual's variance and lag-1 autocorrelation (AC1) over time. +In many examples of the literature, including [boers-early-warning-2018](@cite), +the CSD analysis is performed over segments (sometimes only one) of the timeseries, +such that a significance value is obtained for each segment. By using +`SegmentedWindowConfig`, dealing with segments can be easily done in +TransitionsInTimeseries.jl and is demonstrated here: =# using TransitionsInTimeseries, StatsBase @@ -184,11 +214,32 @@ plot_segment_analysis!(axs, results, signif) fig #= -In [boers-early-warning-2018](@cite), 13/16 and 7/16 true positives are respectively found for the variance and AC1, with 16 referring to the total number of transitions. The timeseries actually includes 18 transition but, in [boers-early-warning-2018](@cite), some segments are considered too small to be analysed. In contrast, we here respectively find 9/16 true positives for the variance and 3/16 for AC1. We can track down the discrepancies to be in the surrogate testing, since the indicator timeseries computed here are almost exactly similar to those of [boers-early-warning-2018](@cite). This mismatch points out that packages like TransitionsInTimeseries.jl are wishful for research to be reproducible, especially since CSD is gaining attention - not only within the scientific community but also in popular media. +In [boers-early-warning-2018](@cite), 13/16 and 7/16 true positives are respectively +found for the variance and AC1, with 16 referring to the total number of transitions. +The timeseries actually includes 18 transition but, in +[boers-early-warning-2018](@cite), some segments are considered too small to be analysed. +In contrast, we here respectively find 9/16 true positives +for the variance and 3/16 for AC1. We can track down the discrepancies to be in the +surrogate testing, since the indicator timeseries computed here are almost exactly +similar to those of [boers-early-warning-2018](@cite). This mismatch points +out that packages like TransitionsInTimeseries.jl are wishful for research to be +reproducible, especially since CSD is gaining attention - not only within the +scientific community but also in popular media. ## CSD: only a necessary condition, only in some cases -For codimension-1 systems, approaching a fold, Hopf or transcritical bifurcation implies a widening of the potential $U$, which defines the deterministic term $f = -∇U$ of the SDE's right-hand-side. In the presence of noise, this leads to CSD, which is therefore a **necessary condition** for crossing one of these bifurcations - although it is not always assessable by analysing the timeseries due to practical limitations (e.g. sparse data subject to large measurement noise). It is nonetheless not given that DO-events, as many other real-life applications, can be seen as a codimension-1 fold, Hopf or transcritical bifurcations. Besides this, we emphasise that CSD is **not a sufficient condition** for assessing a transition being ahead in near future, since a resilience loss can happen without actually crossing any bifurcation. This can be illustrated on the present example by performing the same analysis only until few hundred years before the transition: +For codimension-1 systems, approaching a fold, Hopf or transcritical bifurcation implies +a widening of the potential $U$, which defines the deterministic term $f = -∇U$ of the +SDE's right-hand-side. In the presence of noise, this leads to CSD, which is therefore +a **necessary condition** for crossing one of these bifurcations - although it is not +always assessable by analysing the timeseries due to practical limitations (e.g. sparse +data subject to large measurement noise). It is nonetheless not given that DO-events, +as many other real-life applications, can be seen as a codimension-1 fold, Hopf or +transcritical bifurcations. Besides this, we emphasise that CSD is **not a sufficient +condition** for assessing a transition being ahead in near future, since a resilience +loss can happen without actually crossing any bifurcation. This can be illustrated on +the present example by performing the same analysis only until few hundred years before +the transition: =# tseg_end = t_rasmussen[2:end] .- 700 # stop analysis 500 years earlier than before diff --git a/docs/src/examples/ks_paleojump.jl b/docs/src/examples/ks_paleojump.jl new file mode 100644 index 00000000..bb2dd14f --- /dev/null +++ b/docs/src/examples/ks_paleojump.jl @@ -0,0 +1,168 @@ +# # Kolmogorov-Smirnov test for detecting transitions in paleoclimate timeseries + +# The goal of this example is to show how simple it is to re-create an analysis _similar_ +# to what was done in the paper +# "Automatic detection of abrupt transitions in paleoclimate records", +# [Bagniewski2021](@cite). The same analysis was then used to create a database +# of transitions in paleoclimate records in [Bagniewski2023](@cite) +# Using TransitionsInTimeseries.jl and HypothesisTests.jl, +# the analysis becomes a 10-lines-of-code script (for a given timeseries). + +# ## Scientific background + +# The approach of [Bagniewski2021](@cite) is based on the [two-sample +# Kolmogorov Smirnov test](https://en.wikipedia.org/wiki/Kolmogorov%E2%80%93Smirnov_test#Two-sample_Kolmogorov%E2%80%93Smirnov_test). +# It tests whether the samples from two datasets or timeseries +# are distributed according to the same cumulative density function or not. +# This can be estimated by comparing the value of the _KS-statistic_ +# versus some threshold that depends on the required confidence. + +# The application of this test for identifying transitions in timeseries is simple: + +# 1. A sliding window analysis is performed in the timesries +# 1. In each window, the KS statistic is estimated between the first half and the second +# half of the timeseries within this window. +# 1. Transitions are defined by when the KS statistic exceeds a particular value +# based on some confidence. The transition occurs in the middle of the window. + +# We should point out that in [Bagniewski2021](@cite) the authors did a +# more detailed analysis: analyzed many different window widths, added a +# conditional clause to exclude transitions that do not exceed a predefined minimum +# "jump" in the data, and also added another conditional clause that +# filtered out transitions that are grouped in time (which is a natural consequence +# of using the Kolmogorov-Smirov test for detecting transitions). +# +# Here we won't do that post processing, mainly because it is rather simple +# to include these additional conditional clauses to filter transitions after they are found. + +# ## Steps for TransitionsInTimeseries.jl + +# Doing this kind of work with TransitionsInTimeseries.jl is so easy you won't even trip! +# This analysis follows the same sliding window approach showcased in our [Tutorial](@ref), +# and it even excludes the "indicator" aspect: the change metric is estimated directly +# from the input data! + +# As such, we really only need to define/do these things before we have finished the analysis: + +# 1. Load the input data (we will use the same example as the NGRIP data of +# [Dansgaard-Oescher events and Critical Slowing Down](@ref) example) +# and set the appropriate time window. +# 3. Define the function that estimates the change metric (i.e., the KS-statistic) +# 3. Perform the sliding window analysis as in the [Tutorial](@ref) with [`estimate_indicator_changes`](@ref) +# 4. Estimate the "confident" transitions in the data by comparing the estimated +# KS-statistic with a predefined threshold. + +# ## Load timeseries and window length +# Following the Dansgaard-Oescher events example, we load +# the data after all the processing steps done in that example: + +using DelimitedFiles, CairoMakie + +tmp = Base.download("https://raw.githubusercontent.com/JuliaDynamics/JuliaDynamics/"* + "master/timeseries/NGRIP_processed.csv") +data = readdlm(tmp) +t, xtrend, xresid, xloess = collect.(eachcol(data)) + +fig, ax = lines(t, xtrend; axis = (ylabel = "NGRIP (processed)", xlabel = "time")) +lines!(ax, t, xloess; linewidth = 2) +fig + +# For the window, since we are using a sliding window here, we will be using a +# window of length 500 (which is approximately 1/2 to 1/4 the span between typical +# transitions found by [rasmussen-stratigraphic-2014](@cite)). + +window = 500 + +# ## Defining the change metric function + +# HypothesisTest.jl implements the Kolmogorov-Smirnov test, however here we are interested +# in the value of the test iself (the so-called KS-statistic), rather than a p-value. +# To this end, we define the following function to compute the statistic, +# which also normalizes it as in [Bagniewski2021](@cite). + +using HypothesisTests + +function normalized_KS_statistic(timeseries) + N = length(timeseries) + i = N÷2 + x = view(timeseries, 1:i) + y = view(timeseries, (i+1):N) + kstest = ApproximateTwoSampleKSTest(x, y) + nx = ny = i # length of each timeseries half of total + n = nx*ny/(nx + ny) # written fully for concreteness + D_KS = kstest.δ # can be compared directly with sqrt(-log(α/2)/2) + ## Rescale according to eq. (5) of the paper + rescaled = 1 - ((1 - D_KS)/(1 - sqrt(1/n))) + return rescaled +end + +N = 1000 # the statistic is independent of `N` for large enough `N`! +x = randn(N) +y = 1.8randn(N) .+ 1.0 +z = randn(N) +w = 0.6randn(N) .- 2.0 + +fig, ax = density(x; color = ("black", 0.5), strokewidth = 4.0, label = "reference distribution") +ax.title = "showcase of normalized KS-statistic" +for q in (y, z, w) + D_KS = normalized_KS_statistic(vcat(x, q)) + density!(ax, q; label = "D_KS = $(D_KS)") +end +axislegend(ax) +fig + +# ## Perform the sliding window analysis + +# This is just a straightforward call to [`estimate_indicator_changes`](@ref). +# In fact, it is even simpler than the tutorial. Here we skip completely +# the "indicator" estimation step, and we evaluate the change metric directly +# on input data. We do this by simply passing `nothing` as the indicators. + +using TransitionsInTimeseries + +config = SlidingWindowConfig(nothing, normalized_KS_statistic; width_cha = 500) + +results = estimate_indicator_changes(config, xtrend, t) + +# Which we can visualize +function visualize_results(results) + fig, ax1 = lines(t, xtrend; axis = (ylabel = "NGRIP (processed)",)) + ax2, = lines(fig[2, 1], results.t_change, vec(results.x_change), axis = (ylabel = "D_KS (normalized)", xlabel = "time")) + linkxaxes!(ax1, ax2) + hidexdecorations!(ax1; grid = false) + xloess_normed = (xloess .- minimum(xloess))./(maximum(xloess) - minimum(xloess)) + lines!(ax2, t, xloess_normed; color = ("gray", 0.5)) + fig +end + +visualize_results(results) + +# By overplotting the (smoothened) NGRIP timeseries and the +# normalized KS-statistic, it already becomes pretty clear +# that the statistic peaks when transitions occur. + +# The same thing happens if we alter the window duration + +config = SlidingWindowConfig(nothing, normalized_KS_statistic; width_cha = 200) +results = estimate_indicator_changes(config, xtrend, t) +visualize_results(results) + +# So one can easily obtain extra confidence by varying window +# size as in [Bagniewski2021](@cite). + +# ## Identifying "confident" transitions + +# As this identification here is done via a simple threshold, +# identifying the transitions is a nearly trivial call +# to [`significant_transitions`](@ref) with [`ThresholdSignificance`](@ref) + +signif = ThresholdSignificance(0.5) # or any other threshold +flags = significant_transitions(results, signif) + +fig = visualize_results(results) +axDKS = content(fig[2,1]) +vlines!(axDKS, results.t_change[vec(flags)], color = ("red", 0.25)) +fig + +# We could proceed with a lot of preprocessing as in [Bagniewski2021](@cite) +# but we skip this here for the sake of simplicity. \ No newline at end of file diff --git a/docs/src/refs.bib b/docs/src/refs.bib index 8010a9d3..8cc9690b 100644 --- a/docs/src/refs.bib +++ b/docs/src/refs.bib @@ -52,3 +52,36 @@ @article{rasmussen-stratigraphic-2014 pages = {14--28}, file = {Rasmussen et al. - 2014 - A stratigraphic framework for abrupt climatic chan.pdf:/home/jan/Zotero/storage/VSUTVWPN/Rasmussen et al. - 2014 - A stratigraphic framework for abrupt climatic chan.pdf:application/pdf}, } + +@article{Bagniewski2023, +abstract = {Tipping points (TPs) in Earth's climate system have been the subject of increasing interest and concern in recent years, given the risk that anthropogenic forcing could cause abrupt, potentially irreversible, climate transitions. Paleoclimate records are essential for identifying past TPs and for gaining a thorough understanding of the underlying nonlinearities and bifurcation mechanisms. However, the quality, resolution, and reliability of these records can vary, making it important to carefully select the ones that provide the most accurate representation of past climates. Moreover, as paleoclimate time series vary in their origin, time spans, and periodicities, an objective, automated methodology is crucial for identifying and comparing TPs. To address these challenges, we introduce the open-source PaleoJump database, which contains a collection of carefully selected, high-resolution records originating in ice cores, marine sediments, speleothems, terrestrial records, and lake sediments. These records describe climate variability on centennial, millennial and longer time scales and cover all the continents and ocean basins. We provide an overview of their spatial distribution and discuss the gaps in coverage. Our statistical methodology includes an augmented Kolmogorov–Smirnov test and Recurrence Quantification Analysis; it is applied here, for illustration purposes, to selected records in which abrupt transitions are automatically detected and the presence of potential tipping elements is investigated. These transitions are shown in the PaleoJump database along with other essential information about the records, including location, temporal scale and resolution, as well as temporal plots. This open-source database represents, therefore, a valuable resource for researchers investigating TPs in past climates.}, +archivePrefix = {arXiv}, +arxivId = {2206.06832}, +author = {Bagniewski, Witold and Rousseau, Denis Didier and Ghil, Michael}, +doi = {10.1038/s41598-023-30592-1}, +eprint = {2206.06832}, +isbn = {0123456789}, +issn = {20452322}, +journal = {Scientific Reports}, +number = {1}, +pages = {1--18}, +pmid = {36934110}, +publisher = {Nature Publishing Group UK}, +title = {{The PaleoJump database for abrupt transitions in past climates}}, +url = {https://doi.org/10.1038/s41598-023-30592-1}, +volume = {13}, +year = {2023} +} +@article{Bagniewski2021, +abstract = {Bifurcations and tipping points (TPs) are an important part of the Earth system's behavior. These critical points represent thresholds at which small changes in the system's parameters or in the forcing abruptly switch it from one state or type of behavior to another. Current concern with TPs is largely due to the potential of slow anthropogenic forcing to bring about abrupt, and possibly irreversible, change to the physical climate system and impacted ecosystems. Paleoclimate proxy records have been shown to contain abrupt transitions, or "jumps,"which may represent former instances of such dramatic climate change events. These transitions can provide valuable information for identifying critical TPs in current and future climate evolution. Here, we present a robust methodology for detecting abrupt transitions in proxy records that is applied to ice core and speleothem records of the last climate cycle. This methodology is based on the nonparametric Kolmogorov-Smirnov (KS) test for the equality, or not, of the probability distributions associated with two samples drawn from a time series, before and after any potential jump. To improve the detection of abrupt transitions in proxy records, the KS test is augmented by several other criteria and it is compared with recurrence analysis. The augmented KS test results show substantial skill when compared with more subjective criteria for jump detection. This test can also usefully complement recurrence analysis and improve upon certain aspects of its results.}, +author = {Bagniewski, W. and Ghil, M. and Rousseau, D. D.}, +doi = {10.1063/5.0062543}, +issn = {10897682}, +journal = {Chaos}, +number = {11}, +pmid = {34881579}, +publisher = {AIP Publishing LLC}, +title = {{Automatic detection of abrupt transitions in paleoclimate records}}, +volume = {31}, +year = {2021} +} diff --git a/docs/src/tutorial.md b/docs/src/tutorial.md index 4db96c03..d051e431 100644 --- a/docs/src/tutorial.md +++ b/docs/src/tutorial.md @@ -30,7 +30,7 @@ Let us load data from a bistable nonlinear model subject to noise and to a gradu with $x_{l}$ the state of the linear model, $x_{nl}$ the state of the bistable model, $f$ the forcing and $n$ the noise. For $f=0$ they both display an equilibrium point at $x=-1$. However, the bistable model also displays a further equilibrium point at $x=1$. Loading (and visualizing with [Makie](https://docs.makie.org/stable/)) such prototypical data to test some indicators can be done by simply running: -````@example tutorial +```@example tutorial using TransitionsInTimeseries, CairoMakie t, x_linear, x_nlinear = load_linear_vs_doublewell() @@ -38,7 +38,7 @@ fig, ax = lines(t, x_linear) lines!(ax, t, x_nlinear) ax.title = "raw data" fig -```` +``` ### Preprocessing @@ -49,7 +49,7 @@ fig The nonlinear system clearly displays a transition between two stability regimes. To forecast such transition, we analyze the fluctuations of the timeseries around the attractor, assumed to be tracked. Therefore, a detrending step is needed - here simply obtained by building the difference of the timeseries with lag 1. -````@example tutorial +```@example tutorial x_l_fluct = diff(x_linear) x_nl_fluct = diff(x_nlinear) tfluct = t[2:end] @@ -58,7 +58,7 @@ fig, ax = lines(tfluct, x_l_fluct) lines!(ax, tfluct, x_nl_fluct .+ 0.05) ax.title = "input timeseries" fig -```` +``` At this point, `x_l_fluct` and `x_nl_fluct` are considered the **input timeseries**. @@ -69,7 +69,7 @@ At this point, `x_l_fluct` and `x_nl_fluct` are considered the **input timeserie We can then compute the values of some "indicator" (a Julia function that inputs a timeseries and outputs a number). An indicator should be a quantity that is likely to change if a transition occurs, or is about to occur in the timeseries. We compute indicators by applying a sliding window over the **input timeseries**, determined by the width and the stride with which it is applied. Here we demonstrate this computation with the AR1-regression coefficient (under white-noise assumption), implemented as [`ar1_whitenoise`](@ref): -````@example tutorial +```@example tutorial indicator = ar1_whitenoise indicator_window = (width = 400, stride = 1) @@ -85,7 +85,7 @@ fig, ax = lines(t_indicator, indicator_l) lines!(ax, t_indicator, indicator_nl) ax.title = "indicator timeseries" fig -```` +``` The lines plotted above are the **indicator timeseries**. @@ -93,7 +93,7 @@ The lines plotted above are the **indicator timeseries**. From here, we process the **indicator timeseries** to quantify changes in it. This step is in essence the same as before: we apply some function over a sliding window of the indicator timeseries. We call this new timeseries the **change metric timeseries**. In the example here, the change metric we will employ will be the slope (over a sliding window), calculated via means of a [`RidgeRegressionSlope`](@ref): -````@example tutorial +```@example tutorial change_window = (width = 30, stride = 1) ridgereg = RidgeRegressionSlope(lambda = 0.0) precompridgereg = precompute(ridgereg, t[1:change_window.width]) @@ -106,7 +106,7 @@ fig, ax = lines(t_change, change_l) lines!(ax, t_change, change_nl) ax.title = "change metric timeseries" fig -```` +``` ### Timeseries surrogates @@ -116,7 +116,7 @@ In TransitionsIdentifiers.jl we perform significance testing using the method of To illustrate the surrogate, we compare the change metric computed from the bistable timeseries what that computed from a surrogate of the same timeseries. -````@example tutorial +```@example tutorial # Generate Fourier random-phase surrogates using Random: Xoshiro s = surrogate(x_nl_fluct, RandomFourier(), Xoshiro(123)) @@ -144,13 +144,13 @@ axs[2].title = "change metric" [xlims!(ax, 0, 50) for ax in axs] fig -```` +``` ### Quantifying significance To quantify the significance of the values of the **change metric timeseries** we perform a standard surrogate test by computing the [p-value](https://en.wikipedia.org/wiki/P-value) w.r.t. the change metrics of thousands of surrogates of the input timeseries. A low p-value (typically `p<0.05`) is commonly considered as significant. To visualize significant trends, we plot the p-value vs. time: -````@example tutorial +```@example tutorial n_surrogates = 1_000 fig, axs = gridfig(2, 2) axs[1].title = "linear" @@ -177,7 +177,7 @@ end [xlims!(ax, 0, 50) for ax in axs] fig -```` +``` As expected, the data generated by the nonlinear model displays a significant increase of the AR1-regression coefficient before the transition, which is manifested by a low p-value. In contrast, the data generated by the linear model does not show anything similar. @@ -186,12 +186,11 @@ Performing the step-by-step analysis of transition indicators is possible and mi ## [Tutorial -- TransitionsInTimeseries.jl] (@id example_fastforward) TransitionsInTimeseries.jl wraps this typical workflow into a simple, extendable, and modular API that researchers can use with little effort. In addition, it allows performing the same analysis for several indicators / change metrics in one go. +The interface is simple, and directly parallelizes the [Workflow](@ref). -The interface is simple, and directly parallelizes the [Workflow](@ref). It is based on the creation of a [`TransitionsSurrogatesConfig`](@ref), which contains a list of indicators, and corresponding metrics, to use for doing the above analysis. It also specifies what kind of surrogates to generate. - -The following blocks illustrate how the above extensive example is re-created in TransitionsInTimeseries.jl +The following blocks illustrate how the above extensive example is re-created in TransitionsInTimeseries.jl. But first, let's load the input timeseries: -````@example tutorial +```@example tutorial using TransitionsInTimeseries, CairoMakie t, x_linear, x_nlinear = load_linear_vs_doublewell() @@ -203,19 +202,19 @@ t = t[2:end] fig, ax = lines(t, input) ax.title = "input timeseries" fig -```` +``` -To perform all of the above analysis we follow a 2-step process. +To perform all of the whole [Workflow](@ref) analysis we follow a 2-step process. -Step 1, we decide what indicators and change metrics to use in [`SlidingWindowConfig`](@ref) and apply those via -a sliding window to the input timeseries using [`transition_metrics`](@ref). +**Step 1** is to provide what indicators and change metrics to use in [`SlidingWindowConfig`](@ref) and apply those via +a sliding window to the input timeseries using [`estimate_indicator_changes`](@ref). -````@example tutorial +```@example tutorial # These indicators are suitable for Critical Slowing Down indicators = (var, ar1_whitenoise) # use the ridge regression slope for both indicators -change_metrics = RidgeRegressionSlope() +change_metrics = (RidgeRegressionSlope(), RidgeRegressionSlope()) # choices go into a configuration struct config = SlidingWindowConfig(indicators, change_metrics; @@ -223,11 +222,11 @@ config = SlidingWindowConfig(indicators, change_metrics; # choices are processed results = estimate_indicator_changes(config, input, t) -```` +``` From `result` we can plot the change metric timeseries: -````@example tutorial +```@example tutorial fig, axs = gridfig(3, 1) lines!(axs[1], t, input; label = "input", color = Cycled(2)) scatter!(axs[2], results.t_change, results.x_change[:, 1]; @@ -236,19 +235,19 @@ scatter!(axs[3], results.t_change, results.x_change[:, 2]; label = "ar1 slopes", color = Cycled(4)) [xlims!(ax, 0, 50) for ax in axs] fig -```` +``` -Step 2 is to estimate significance using [`SurrogatesConfig`](@ref) -and the function [`estimate_significance!`](@ref). +**Step 2** is to estimate significance using [`SurrogatesSignificance`](@ref) +and the function [`significant_transitions`](@ref). -````@example tutorial +```@example tutorial signif = SurrogatesSignificance(n = 1000, tail = :right) flags = significant_transitions(results, signif) -```` +``` We can now plot the p-values corresponding to each time series of the change metrics. From the `flags` we can additionally obtain the time points where _both_ indicators show significance, via a simple reduction: -````@example tutorial +```@example tutorial fig, axs = gridfig(2, 1) lines!(axs[1], vcat(0.0, t), x_nlinear; label = "raw", color = Cycled(1)) lines!(axs[1], t, input; label = "input", color = Cycled(2)) @@ -263,5 +262,4 @@ vlines!(axs[1], results.t_change[flagsboth]; label = "flags", color = ("black", [axislegend(ax) for ax in axs] [xlims!(ax, 0, 50) for ax in axs] fig -```` - +``` diff --git a/ext/TITVisualizations.jl b/ext/TITVisualizations.jl new file mode 100644 index 00000000..2a2f8435 --- /dev/null +++ b/ext/TITVisualizations.jl @@ -0,0 +1,69 @@ +module TITVisualizations + +using TransitionsInTimeseries, Makie + +function TransitionsInTimeseries.plot_indicator_changes(res::SlidingWindowResults, + colors = ["#7143E0", "#0A9A84", "#191E44", "#AF9327", "#701B80", "#2E6137",], + skip = Int[], + ) + config = res.config + fig = Figure() + axts = Axis(fig[1,1]; ylabel = "input") + if isnothing(config.indicators) + axcha = Axis(fig[2,1]; ylabel = "changes") + linkxaxes!(axts, axcha) + else + axind = Axis(fig[2,1]; ylabel = "indicators") + axcha = Axis(fig[3,1]; ylabel = "changes", xlabel = "time") + linkxaxes!(axts, axind, axcha) + hidexdecorations!(axind; grid = false) + end + hidexdecorations!(axts; grid = false) + + lines!(axts, res.t, res.x; color = "black", linewidth = 3) + for (i, cha) in enumerate(config.change_metrics) + i ∈ skip && continue + if !isnothing(config.indicators) + lines!(axind, res.t_indicator, res.x_indicator[:, i]; + color = colors[i], linewidth = 3, label = string(nameof(config.indicators[i])) + ) + end + lines!(axcha, res.t_change, res.x_change[:, i]; + color = colors[i], linewidth = 3, label = string(nameof(cha)) + ) + end + + !isnothing(config.indicators) && axislegend(axind) + axislegend(axcha) + return fig +end + +function TransitionsInTimeseries.plot_significance!(fig::Figure, res::SlidingWindowResults, signif::SurrogatesSignificance; + colors = ["#7143E0", "#0A9A84", "#191E44", "#AF9327", "#701B80", "#2E6137",], + skip = Int[], + nsurro = 20, + ) + config = res.config + for (i, cha) in enumerate(config.change_metrics) + i ∈ skip && continue + + for _ in 1:nsurro + s = TimeseriesSurrogates.surrogate(res.x, signif.surrogate) + if isnothing(config.indicators) + p = s + chaj = 2 + else + p = windowmap(config.indicators[i], s; width = config.width_ind, stride = config.stride_ind) + lines!(fig[2, 1], res.t_indicator, q; color = (colors[i], 2/nsurro), linewidth = 1) + chaj = 3 + end + q = windowmap(cha, p; width = config.width_cha, stride = config.stride_cha) + lines!(fig[chaj, 1], res.t_change, q; color = (colors[i], 2/nsurro), linewidth = 1) + end + end + + return fig +end + + +end \ No newline at end of file diff --git a/src/TransitionsInTimeseries.jl b/src/TransitionsInTimeseries.jl index d191ba5d..97715dc4 100644 --- a/src/TransitionsInTimeseries.jl +++ b/src/TransitionsInTimeseries.jl @@ -28,7 +28,7 @@ include("analysis/sliding_window.jl") include("analysis/segmented_window.jl") include("significance/api_significance.jl") include("significance/surrogates_significance.jl") -include("significance/quantile_significance.jl") +include("significance/basic_stat_significance.jl") include("indicators/critical_slowing_down.jl") include("indicators/distribution_distance.jl") @@ -39,6 +39,7 @@ include("indicators/statistics.jl") include("change_metrics/slope.jl") include("change_metrics/valuediff.jl") +include("visualizations.jl") # windowing.jl export WindowViewer, windowmap, windowmap!, midpoint, midvalue @@ -58,7 +59,7 @@ export IndicatorsChangesConfig, SlidingWindowConfig, SegmentedWindowConfig export SlidingWindowResults, SegmentedWindowResults export estimate_indicator_changes, IndicatorsChangesResults export TransitionsSignificance, significant_transitions, segmented_significance -export QuantileSignificance, SigmaSignificance, SurrogatesSignificance +export ThresholdSignificance, QuantileSignificance, SigmaSignificance, SurrogatesSignificance # timeseries export isequispaced, equispaced_step @@ -66,4 +67,7 @@ export isequispaced, equispaced_step # load_data.jl export load_linear_vs_doublewell +# visualizations +export plot_indicator_changes, plot_significance! + end # module TransitionsInTimeseries \ No newline at end of file diff --git a/src/analysis/segmented_window.jl b/src/analysis/segmented_window.jl index 5b28c8df..ddb85cb3 100644 --- a/src/analysis/segmented_window.jl +++ b/src/analysis/segmented_window.jl @@ -157,10 +157,10 @@ end # Segmented and Sliding results share their show method function Base.show(io::IO, ::MIME"text/plain", res::Union{SegmentedWindowResults, SlidingWindowResults}) - println(io, "IndicatorsChangesResults") + println(io, nameof(typeof(res))) descriptors = [ "input timeseries" => summary(res.x), - "indicators" => [nameof(i) for i in res.config.indicators], + "indicators" => isnothing(res.config.indicators) ? "nothing" : [nameof(i) for i in res.config.indicators], "indicator (window, stride)" => (res.config.width_ind, res.config.stride_ind), "change metrics" => [nameof(c) for c in res.config.change_metrics], show_changemetric(res), diff --git a/src/analysis/sliding_window.jl b/src/analysis/sliding_window.jl index 25741468..1602223b 100644 --- a/src/analysis/sliding_window.jl +++ b/src/analysis/sliding_window.jl @@ -9,18 +9,22 @@ It estimates transitions by a sliding window approach: 2. Estimate changes of an indicator by sliding a window of the change metric over the indicator timeseries. -`indicators` is a tuple of indicators (or a single indicator). -`change_metrics` is also a tuple or a single function. If a single function, -the same change metric is used for all provided indicators. This way the analysis -can be efficiently repeated for many indicators and/or change metrics. - Both indicators and change metrics are **generic Julia functions** that input an `x::AbstractVector` and output an `s::Real`. Any function may be given and see [making custom indicators/change metrics](@ref own_indicator) in the documentation for more information on possible optimizations. +`indicators` can be a single function or a tuple of indicators. +Similarly, `change_metrics` can be a tuple or a single function. +If tuples, the length of `indicators` and `change_metrics` must match. +This way the analysis +can be efficiently repeated for many indicators and/or change metrics. + The results output corresponding to `SlidingWindowConfig` is [`SlidingWindowResults`](@ref). +Step 1. is skipped if `nothing` is provided as `indicators`, in which case +the change metrics are estimated directly from input data. + ## Keyword arguments - `width_ind::Int=100, stride_ind::Int=1`: width and stride given to [`WindowViewer`](@ref) @@ -63,7 +67,9 @@ function SlidingWindowConfig( ) indicators, change_metrics = sanitycheck_metrics(indicators, change_metrics) # Last step: precomputable functions, if any - indicators = map(f -> precompute(f, 1:T(width_ind)), indicators) + if !isnothing(indicators) + indicators = map(f -> precompute(f, 1:T(width_ind)), indicators) + end change_metrics = map(f -> precompute(f, 1:T(width_cha)), change_metrics) return SlidingWindowConfig( @@ -73,46 +79,58 @@ function SlidingWindowConfig( end function sanitycheck_metrics(indicators, change_metrics) - if !(indicators isa Tuple) - indicators = (indicators,) - end if !(change_metrics isa Tuple) change_metrics = (change_metrics,) end - L = length(indicators) - if length(change_metrics) ∉ (1, L) - throw(ArgumentError("The amount of change metrics must be as many as the"* - "indicators, or only 1.")) + if indicators isa Function + indicators = (indicators, ) + end + if !isnothing(indicators) && (length(change_metrics) ≠ length(indicators)) + throw(ArgumentError("The amount of change metrics and indicators must match.")) end return indicators, change_metrics end +# TODO: This function needs to be split into two as per good practices +# for performant code. One function does the initialization, +# and the other overwrites in place. This makes it easier +# to apply for surrogates as well! + function estimate_indicator_changes(config::SlidingWindowConfig, x, t = eachindex(x)) + (; indicators, change_metrics) = config # initialize time vectors - t_indicator = windowmap(config.whichtime, t; width = config.width_ind, - stride = config.stride_ind) - t_change = windowmap(config.whichtime, t_indicator; width = config.width_cha, - stride = config.stride_cha) + if isnothing(indicators) + # Skip indicators if they are nothing + t_indicator = t + else + t_indicator = windowmap(config.whichtime, t; width = config.width_ind, + stride = config.stride_ind) + end + t_change = windowmap(config.whichtime, t_indicator; + width = config.width_cha, stride = config.stride_cha + ) len_ind = length(t_indicator) len_change = length(t_change) # initialize array containers - (; indicators, change_metrics) = config X = eltype(x) - n_ind = length(indicators) - x_indicator = zeros(X, len_ind, n_ind) - x_change = zeros(X, len_change, n_ind) - one2one = length(change_metrics) == length(indicators) - - # Loop over indicators - for i in 1:n_ind - indicator = config.indicators[i] - chametric = one2one ? change_metrics[i] : change_metrics[1] - z = view(x_indicator, :, i) - windowmap!(indicator, z, x; - width = config.width_ind, stride = config.stride_ind - ) - windowmap!(chametric, view(x_change, :, i), z; + n_metrics = length(change_metrics) + x_indicator = zeros(X, len_ind, n_metrics) + x_change = zeros(X, len_change, n_metrics) + + for i in 1:n_metrics + # estimate indicator timeseries + if !isnothing(indicators) + z = view(x_indicator, :, i) + windowmap!(indicators[i], z, x; + width = config.width_ind, stride = config.stride_ind + ) + else + # Just equate x here, we're skipping the indicator estimation + z = x + end + # and then apply the change metric + windowmap!(change_metrics[i], view(x_change, :, i), z; width = config.width_cha, stride = config.stride_cha ) end @@ -143,11 +161,11 @@ It has the following fields that the user may access - `config::SlidingWindowConfig`, what was used for the analysis. """ -struct SlidingWindowResults{TT, T<:Real, X<:Real, XX<:AbstractVector{X}, +struct SlidingWindowResults{TT, T<:Real, X<:Real, XX<:AbstractVector{X}, IT, W} <: IndicatorsChangesResults t::TT # original time vector; most often it is `Base.OneTo`. x::XX - t_indicator::Vector{T} + t_indicator::IT x_indicator::Matrix{X} t_change::Vector{T} x_change::Matrix{X} diff --git a/src/significance/api_significance.jl b/src/significance/api_significance.jl index a81306dd..7a381884 100644 --- a/src/significance/api_significance.jl +++ b/src/significance/api_significance.jl @@ -5,6 +5,9 @@ Supertype used to test for significance in [`significant_transitions`](@ref). Valid subtypes are: - [`SurrogatesSignificance`](@ref). +- [`SigmaSignificance`](@ref). +- [`QuantileSignificance`](@ref). +- [`ThresholdSignificance`](@ref). """ abstract type TransitionsSignificance end diff --git a/src/significance/quantile_significance.jl b/src/significance/basic_stat_significance.jl similarity index 71% rename from src/significance/quantile_significance.jl rename to src/significance/basic_stat_significance.jl index 332cf5ac..8ea52857 100644 --- a/src/significance/quantile_significance.jl +++ b/src/significance/basic_stat_significance.jl @@ -1,3 +1,39 @@ +""" + ThresholdSignificance(threshold::Real; tail = :right) <: TransitionsSignificance + +A configuration struct for significance testing in [`significant_transitions`](@ref). +Significance is estimated by comparing the value of each change metric with +the given `threshold`. +Values that exceed the threshold (if `tail = :right`) +or subseed the threshold (if `tail = :left`) +are deemed significant. +If `tail = :both` then either condition is checked. +""" +struct ThresholdSignificance{T<:Real} + threshold::T + tail::Symbol +end +ThresholdSignificance(threshold; tail = :right) = ThresholdSignificance(threshold, tail) + +function significant_transitions(res::IndicatorsChangesResults, signif::ThresholdSignificance) + flags = similar(res.x_change, Bool) + t = signif.threshold + for (i, x) in enumerate(eachcol(res.x_change)) + flag = view(flags, :, i) + if signif.tail == :right + @. flag = x > t + elseif signif.tail == :left + @. flag = x < t + elseif signif.tail == :both + @. flag = (x < t) | (x > t) + else + error("`tail` can be only `:left, :right, :both`. Got $(tail).") + end + end + return flags +end + + """ QuantileSignificance(; p = 0.95, tail = :both) <: TransitionsSignificance diff --git a/src/significance/surrogates_significance.jl b/src/significance/surrogates_significance.jl index dfdbece0..3e930cc2 100644 --- a/src/significance/surrogates_significance.jl +++ b/src/significance/surrogates_significance.jl @@ -1,6 +1,6 @@ """ SurrogatesSignificance <: TransitionsSignificance - SurrogatesSignificance(; surrogate = RandomFourier(), n = 10_000, tail = :both, rng) + SurrogatesSignificance(; kwargs...) A configuration struct for significance testing [`significant_transitions`](@ref) using timeseries surrogates. @@ -30,7 +30,7 @@ in the field `.pvalues` (to e.g., threshold with different `p`). The p-value is simply the proportion of surrogate change metric values that exceed (for `tail = :right`) or subseed (`tail = :left`) the original change metric at each given time point. Use `tail = :left` if the surrogate data are expected to have -higher change metric, discriminatory statistic values. This is the case for statistics +higher change metric values. This is the case for statistics that quantify entropy. For statistics that quantify autocorrelation, use `tail = :right` instead. For anything else, use the default `tail = :both`. An iterable of `tail` values can also be given, in which case a specific `tail` @@ -51,7 +51,7 @@ end function SurrogatesSignificance(; surromethod = RandomFourier(), - n::Int = 10_000, + n::Int = 1_000, tail = :both, rng = Random.default_rng(), p = 0.05) return SurrogatesSignificance(surromethod, n, tail, rng, p, zeros(1,1)) @@ -62,9 +62,9 @@ function significant_transitions(res::SlidingWindowResults, signif::SurrogatesSi (; x, x_indicator, x_change) = res (; indicators, change_metrics, width_ind, stride_ind, width_cha, stride_cha) = res.config (; surrogate, n, tail, rng, p, pvalues) = signif - n_ind = length(indicators) + n_ind = length(change_metrics) sanitycheck_tail(tail, n_ind) - + # Init pvalues X = eltype(x_change) pvalues = signif.pvalues = zeros(X, size(x_change)) @@ -75,11 +75,14 @@ function significant_transitions(res::SlidingWindowResults, signif::SurrogatesSi seeds = rand(rng, 1:typemax(Int), Threads.nthreads()) sgens = [surrogenerator(x, surrogate, Xoshiro(seed)) for seed in seeds] # Dummy vals for surrogate parallelization - indicator_dummys = [x_indicator[:, 1] for _ in 1:Threads.nthreads()] + if !isnothing(indicators) + indicator_dummys = [x_indicator[:, 1] for _ in 1:Threads.nthreads()] + else + indicator_dummys = [copy(x) for _ in 1:Threads.nthreads()] + end change_dummys = [x_change[:, 1] for _ in 1:Threads.nthreads()] Threads.@threads for _ in 1:n - id = Threads.threadid() s = sgens[id]() change_dummy = change_dummys[id] @@ -87,12 +90,16 @@ function significant_transitions(res::SlidingWindowResults, signif::SurrogatesSi for i in 1:n_ind indicator, change_metric, tai = choose_metrics(indicators, change_metrics, tail, i) - c = view(x_change, :, i) # change metric timeseries + c = view(x_change, :, i) # change metric timeseries pval_right = view(pvals_right, :, i) pval_left = view(pvals_left, :, i) - windowmap!(indicator, indicator_dummys[id], s; + if !isnothing(indicator) + windowmap!(indicator, indicator_dummys[id], s; width = width_ind, stride = stride_ind) - windowmap!(change_metric, change_dummys[id], indicator_dummys[id]; + # Skip the indicator step if not provided + # (we don't need a clause, this is already the `x` timeseries) + end + windowmap!(change_metric, change_dummy, indicator_dummys[id]; width = width_cha, stride = stride_cha) accumulate_pvals!(pval_right, pval_left, tai, c, change_dummy) end @@ -164,10 +171,10 @@ function sanitycheck_tail(tail, n_ind) end end -function choose_metrics(indicators::Tuple{Vararg}, change_metrics, tail, i::Int) - i_metric = length(change_metrics) == length(indicators) ? i : 1 +function choose_metrics(indicators, change_metrics, tail, i::Int) + ind = isnothing(indicators) ? nothing : indicators[i] tai = tail isa Symbol ? tail : tail[i] - return indicators[i], change_metrics[i_metric], tai + return ind, change_metrics[i], tai end function accumulate_pvals!(pval_right, pval_left, tail, c, change_dummy) diff --git a/src/visualizations.jl b/src/visualizations.jl new file mode 100644 index 00000000..02cdc637 --- /dev/null +++ b/src/visualizations.jl @@ -0,0 +1,9 @@ +function plot_changes_significance(res::IndicatorsChangesConfig, signif::TransitionsSignificance) + fig = plot_indicator_changes() + plot_significance!() + return fig +end + +function plot_indicator_changes end + +function plot_significance! end \ No newline at end of file diff --git a/test/full_analysis.jl b/test/full_analysis.jl index ace87c44..811dc556 100644 --- a/test/full_analysis.jl +++ b/test/full_analysis.jl @@ -1,4 +1,4 @@ -using TransitionsInTimeseries, Test +using TransitionsInTimeseries, Test, Statistics, Random @testset "sliding ridge regression" begin n = 1001 @@ -6,7 +6,7 @@ using TransitionsInTimeseries, Test x = copy(t) indicators = (mean, var) - change_metric = RidgeRegressionSlope() + change_metric = (RidgeRegressionSlope(), RidgeRegressionSlope()) config = SlidingWindowConfig(indicators, change_metric; width_ind = 100, stride_ind = 1, @@ -28,4 +28,51 @@ using TransitionsInTimeseries, Test @test all(.!flags) + @testset "missmatch in length" begin + indicators = (mean, var) + change_metric = RidgeRegressionSlope() + @test_throws ArgumentError SlidingWindowConfig(indicators, change_metric) + end + +end + +@testset "nothing indicator" begin + rng = Xoshiro(123) + + N = 1000 # the statistic is independent of `N` for large enough `N`! + x = randn(rng, N) + y = 1.8randn(rng, N) .+ 4.0 + + x = vcat(x, y) + + function difference_of_maxes(x::AbstractArray) + # assumes 1-based indexing + n = length(x) + x1 = view(x, 1:n÷2) + x2 = view(x, (n÷2 + 1):n) + return abs(maximum(x1) - maximum(x2)) + end + + config = SlidingWindowConfig(nothing, (difference_of_means, difference_of_maxes); + width_cha = 100, stride_cha = 50 + ) + + res = estimate_indicator_changes(config, x) + + c = res.x_change[:, 1] + + @test maximum(c) > 1 + @test count(>(1), c) == 1 + + # surrogates + signif = SurrogatesSignificance(surromethod = RandomShuffle(), n = 1000, tail = :left, p = 0.05) + flags = significant_transitions(res, signif) + + @test count(flags) == 1 + + # Test the surrogate results + # using CairoMakie + # fig = plot_indicator_changes(res) + # plot_significance!(fig, res, signif) + end \ No newline at end of file