report.Rmd

---
title: "Report"
author: "Francisco Bischoff"
date: "on Jun 06, 2023"
output:
  bookdown::html_document2:
    base_format: workflowr::wflow_html
    toc: true
    fig_caption: yes
    number_sections: yes
bibliography: [../papers/references.bib]
link-citations: true
csl: ../thesis/csl/ama.csl
css: style.css
editor_options:
  chunk_output_type: console
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(
  echo = FALSE, fig.align = "center", autodep = TRUE,
  fig.height = 5, fig.width = 10,
  tidy = "styler",
  tidy.opts = list(strict = TRUE)
)

if (knitr::is_latex_output()) {
  knitr::opts_chunk$set(dev = "pdf")
} else {
  knitr::opts_chunk$set(dev = "svg")
}

my_graphics <- function(image_name, base_path = here::here("docs", "figure")) {
  file_path <- glue::glue("{base_path}/{image_name}")

  if (knitr::is_latex_output()) {
    if (file.exists(glue::glue("{file_path}.pdf"))) {
      file_path <- glue::glue("{file_path}.pdf")
    } else if (file.exists(glue::glue("{file_path}.png"))) {
      file_path <- glue::glue("{file_path}.png")
    } else {
      file_path <- glue::glue("{file_path}.jpg")
    }
  } else {
    if (file.exists(glue::glue("{file_path}.svg"))) {
      file_path <- glue::glue("{file_path}.svg")
    } else if (file.exists(glue::glue("{file_path}.png"))) {
      file_path <- glue::glue("{file_path}.png")
    } else {
      file_path <- glue::glue("{file_path}.jpg")
    }
  }

  knitr::include_graphics(file_path)
}

my_kable <- function(title, label, content) {
  res <- glue(r"(<br><table class="tg"><caption>)", "(\\#tab:{label}) {title}", r"(</caption>{content}</table>)")
  out <- structure(res, format = "html", class = "knitr_kable")
  attr(out, "format") <- "html"
  out
}

tkplot <- function(object, interactive = FALSE, res = 50) {
  ecg <- read_ecg_with_atr(here::here("inst/extdata/afib_regimes", object$record), resample_from = 200, resample_to = res)
  value <- ecg[[1]]$II
  prop <- 250 / res
  mask <- seq.int(50, 100)
  value[1:5] <- median(value[mask])
  value[(length(value) - 5):length(value)] <- median(value[mask])
  time <- seq(1, floor(length(value) * prop), length.out = length(value))
  data <- tibble::tibble(time = time, value = value)
  min_data <- min(data$value)
  max_data <- max(data$value)
  truth <- clean_truth(floor(attr(ecg[[1]], "regimes") * prop), floor(length(value) * prop)) # object$truth[[1]]
  preds <- object$pred[[1]]

  title <- glue::glue(
    "Recording: {object$record} ",
    "#truth: {length(truth)}, ",
    "#preds: {length(preds)}, ",
    "length: {floor(length(value)*prop)} ",
    "FLOSS Score: {round(object$score, 3)}"
  )

  subtitle <- glue::glue(
    "Parameters: ",
    "MP window: {object$window_size}, ",
    "MP threshold: {object$mp_threshold}, ",
    "Time constraint: {object$time_constraint}, ",
    "Regime threshold: {object$regime_threshold}, ",
    "Regime landmark: {object$regime_landmark}"
  )


  plot <- data %>%
    timetk::plot_time_series(
      time, value,
      .title = glue::glue(title, "<br><sup>{subtitle}</sup>"),
      .interactive = interactive,
      .smooth = FALSE,
      .line_alpha = 0.3,
      .line_size = 0.2,
      .plotly_slider = interactive
    )

  if (interactive) {
    plot <- plot %>%
      plotly::add_segments(
        x = preds, xend = preds, y = min_data,
        yend = max_data * 1.1,
        line = list(width = 2.5, color = "#0108c77f"),
        name = "Predicted"
      ) %>%
      plotly::add_segments(
        x = truth, xend = truth, y = min_data,
        yend = max_data,
        line = list(width = 2.5, color = "#ff00007f"),
        name = "Truth"
      )
  } else {
    plot <- plot +
      ggplot2::geom_segment(
        data = tibble::tibble(tru = truth),
        aes(
          x = tru, xend = tru,
          y = min_data, yend = max_data - (max_data - min_data) * 0.1
        ), linewidth = 2, color = "#ff00007f"
      ) +
      ggplot2::geom_segment(
        data = tibble::tibble(pre = preds),
        aes(
          x = pre, xend = pre,
          y = min_data, yend = max_data
        ), linewidth = 1, color = "#0108c77f"
      ) +
      ggplot2::theme_bw() +
      ggplot2::theme(
        legend.position = "none",
        plot.margin = margin(0, 0, 0, 10)
      ) +
      ggplot2::labs(title = title, subtitle = subtitle, y = ggplot2::element_blank())
  }
  plot
}

options(dplyr.summarise.inform = FALSE)

# tinytex::tlmgr_install("tabu")
library(here)
library(glue)
library(visNetwork)
library(tibble)
library(kableExtra)
library(patchwork)
library(targets)
library(ggplot2)

source(here::here("scripts", "common", "read_ecg.R"))

# knitr::opts_knit$set(
#   root.dir = here("docs"),
#   base.dir = here("protocol"),
#   verbose = TRUE
# )
```

Last Updated: 2023-06-12 12:51:54 UTC

# Objectives and the research question

While this research was inspired by the CinC/Physionet Challenge 2015, its purpose is not to beat
the state-of-the-art on that challenge, but to identify, on streaming data, abnormal hearth electric
patterns, specifically those which are life-threatening, using low CPU and low memory requirements
in order to be able to generalize the use of such information on lower-end devices, outside the ICU, as ward devices, home devices, and wearable devices.

The main question is: can we accomplish this objective using a minimalist approach (low CPU, low
memory) while maintaining robustness?

# Principles

This research is being conducted using the Research Compendium principles [@compendium2019]:

1. Stick with the convention of your peers;
2. Keep data, methods, and output separated;
3. Specify your computational environment as clearly as you can.

Data management follows the FAIR principle (findable, accessible, interoperable, reusable)
[@wilkinson2016]. Concerning these principles, the dataset was converted from Matlab's format to
CSV format, allowing more interoperability. Additionally, all the project, including the dataset, is
in conformity with the Codemeta Project [@CodeMeta2017].

# Materials and methods

## Softwares

### Pipeline management

All process steps are managed using the R package `targets` [@landau2021], from data extraction to
the final report. An example of a pipeline visualization created with `targets` is shown in Fig.
\@ref(fig:targets). This package helps to record the random seeds (allowing reproducibility),
changes in some part of the code (or dependencies), and then run only the branches that need to be
updated and several other features to keep a reproducible workflow avoiding unnecessary repetitions.

```{r targets, echo=FALSE, out.width="100%"}
#| fig.cap="Example of pipeline visualization using `targets`.
#|  From left to right we see 'Stems' (steps that do not create branches) and 'Patterns'
#|  (that contains two or more branches) and the flow of the information.
#|  The green color means that the step is up to date to the current code and dependencies."

my_graphics("targets", "figure")
```

### Reports management

The report is available on the main webpage [@franz_website], allowing inspection of previous
versions managed by the R package `workflowr`[@workflowr2021]. This package complements the
`targets` package by taking care of the versioning of every report. It is like a Log Book that keeps
track of every important milestone of the project while summarizing the computational environment
where it was run. Fig. \@ref(fig:workflowr) shows only a fraction of the generated website,
where we can see that this version passed the required checks (the system is up-to-date, no caches,
session information was recorded, and others), and we see a table of previous versions.

```{r workflowr, echo=FALSE, out.width="100%"}
#| fig.cap="Fraction of the website generated by `workflowr`.
#|  On top we see that this version passed all checks, and in the middle we see a table
#|  referring to the previous versions of the report."
my_graphics("workflowr_print", "figure")
```

### Modeling and parameter tuning

The well-known package used for data science in R is the `caret` (short for **C**lassification
**A**nd **RE**gression **T**raining) [@JSSv028i05]. Nevertheless, the author of `caret` recognizes
several limitations of his (great) package and is now in charge of developing the `tidymodels`
[@tidymodels2020] collection. For sure, there are other available frameworks and opinions
[@Thompson2020]. Notwithstanding, this project will follow the `tidymodels` road. Three significant
arguments 1) constantly improving and constantly being re-checked for bugs; large community
contribution; 2) allows to plug in a custom modeling algorithm that, in this case, will be the one
needed for developing this work; 3) `caret` is not in active development.

### Continuous integration

Meanwhile, the project pipeline has been set up on GitHub, Inc. [@bischoffrepo2021], leveraging
Github Actions [@gitactions2021] for the Continuous Integration lifecycle. The repository is
available at [@bischoffrepo2021], and the resulting report is available at [@franz_website]. This
thesis's roadmap and tasks status are also publicly available on Zenhub [@zenhub2021].

## Developed software

### Matrix Profile {#matrixprofile}

Matrix Profile (MP) [@Yeh2017a] is a state-of-the-art [@DePaepe2020; @Feremans2020] time series
analysis technique that, once computed, allows us to derive frameworks to all sorts of tasks, as
motif discovery, anomaly detection, regime change detection, and others [@Yeh2017a].

Before MP, time series analysis relied on the *distance matrix* (DM), a matrix that stores all the
distances between two time series (or itself, in the case of a Self-Join). This was very
power-consuming, and several pruning and dimensionality reduction methods were researched
[@Lin2007].

For brevity, let's just understand that the MP and the companion Profile Index (PI) are two vectors
that hold one floating-point value and one integer value, respectively, regarding the original time
series: (1) the similarity distance between that point on time (let's call these points "indexes")
and its first nearest-neighbor (1-NN), (2) The index where this 1-NN is located. The original paper
has more detailed information [@Yeh2017a]. It is computed using a rolling window, but instead of
creating a whole DM, only the minimum values and the index of these minima are stored (in the MP and
PI, respectively). We can have an idea of the relationship of both on Fig. \@ref(fig:thematrix).

```{r thematrix, echo=FALSE}
#| fig.cap="A distance matrix (top), and a matrix profile (bottom). The matrix profile stores only
#|  the minimum values of the distance matrix."
my_graphics("mp_1", "figure")
```

This research has already yielded two R packages concerning the MP algorithms from UCR [@mpucr]. The
first package is called `tsmp`, and a paper has also been published in the R Journal [@RJ-2020-021]
(Journal Impact Factor™, 2020 of 3.984). The second package is called `matrixprofiler` and enhances
the first one, using low-level language to improve computational speed. The author has also joined
the Matrix Profile Foundation as a co-founder with contributors from Python and Go languages
[@mpf2020; @VanBenschoten2020].

This implementation in R is being used for computing the MP and MP-based algorithms of this thesis.

## The data {#the-data}

The current dataset used is the CinC/Physionet Challenge 2015 public dataset, modified to include
only the actual data and the header files in order to be read by the pipeline and is hosted by
Zenodo [@bischoff2021] under the same license as Physionet.

The dataset is composed of 750 patients with at least five minutes records. All signals have been
resampled (using anti-alias filters) to 12 bit, 250 Hz, FIR band-pass (0.05 to 40Hz), and mains
notch filters applied to remove noise. Pacemaker and other artifacts are still present on the ECG
[@Clifford2015]. Furthermore, this dataset contains at least two ECG derivations and one or more
variables like arterial blood pressure, photoplethysmograph readings, and respiration movements.

The _events_ we seek to identify are the life-threatening arrhythmias as defined by Physionet in
Table \@ref(tab:alarms).

```{r alarms, echo=FALSE}
alarms <- tribble(
  ~Alarm, ~Definition,
  "Asystole", "No QRS for at least 4 seconds",
  "Extreme Bradycardia", "Heart rate lower than 40 bpm for 5 consecutive beats",
  "Extreme Tachycardia", "Heart rate higher than 140 bpm for 17 consecutive beats",
  "Ventricular Tachycardia", "5 or more ventricular beats with heart rate higher than 100 bpm",
  "Ventricular Flutter/Fibrillation", "Fibrillatory, flutter, or oscillatory waveform for at least 4 seconds"
)

kbl(alarms,
  booktabs = TRUE,
  caption = "Definition of the five alarm types used in CinC/Physionet Challenge 2015.",
  align = "ll",
  position = "ht",
  linesep = "\\addlinespace"
) |>
  row_spec(0, bold = TRUE) |>
  kable_styling(full_width = TRUE)
```

The fifth minute is precisely where the alarm has been triggered on the original recording set. To
meet the ANSI/AAMI EC13 Cardiac Monitor Standards [@AAMI2002], the onset of the event is within 10
seconds of the alarm (i.e., between 4:50 and 5:00 of the record). That doesn't mean that there have
been no other arrhythmias before.

For comparison, on Table \@ref(tab:challenge) we collected the score of the five best participants
of the challenge [@plesinger2015; @kalidas2015; @couto2015; @fallet2015; @hoogantink2015].

```{r challenge, echo=FALSE}
challenge <- tribble(
  ~Score, ~Authors,
  "81.39", "Filip Plesinger, Petr Klimes, Josef Halamek, Pavel Jurak",
  "79.44", "Vignesh Kalidas",
  "79.02", "Paula Couto, Ruben Ramalho, Rui Rodrigues",
  "76.11", "Sibylle Fallet, Sasan Yazdani, Jean-Marc Vesin",
  "75.55", "Christoph Hoog Antink, Steffen Leonhardt"
)

kbl(challenge,
  booktabs = TRUE,
  caption = "Challenge Results on real-time data. The scores were multiplied by 100.",
  align = "cl",
  position = "ht"
) |>
  row_spec(0, bold = TRUE) |>
  # column_spec(1, width = "5em") |>
  # column_spec(2, width = "30em") |>
  kable_styling(full_width = TRUE)
```

The Equation used on this challenge to compute the score of the algorithms is in the Equation
$\eqref{score}$. This Equation is the accuracy formula, with the penalization of the false
negatives. The reasoning pointed out by the authors [@Clifford2015] is the clinical impact of
existing an actual life-threatening event that was considered unimportant. Accuracy is known to be
misleading when there is a high class imbalance [@Akosa2017].

\
$$
Score = \frac{TP+TN}{TP+TN+FP+5*FN} \tag{1} \label{score}
$$
\

Assuming that this is a finite dataset, the pathologic cases (1) $\lim_{TP \to \infty}$ (whenever
there is an event, it is positive) or (2) $\lim_{TN \to \infty}$ (whenever there is an event, it is
false), cannot happen. This dataset has 292 True alarms and 458 False alarms. Experimentally, this
equation yields:

  - 0.24 if all guesses are on False class
  - 0.28 if random guesses
  - 0.39 if all guesses are on True class
  - 0.45 if no false positives plus random on True class
  - 0.69 if no false negatives plus random on False class

This small experiment (knowing the data in advance) shows that "a single line of code and a few
minutes of effort" [@Wu2020] algorithm could achieve at most a score of 0.39 in this challenge (the
last two lines, the algorithm must be very good on one class).

Nevertheless, this Equation will only be helpful to allow us to compare the results of this thesis
with other algorithms.

## Work structure

### Project start

The project started with a literature survey on the databases Scopus, PubMed, Web of Science, and
Google Scholar with the following query (the syntax was adapted for each database):

\
    TITLE-ABS-KEY ( algorithm  OR  'point of care'  OR  'signal processing'  OR  'computer
    assisted'  OR  'support vector machine'  OR  'decision support system*'  OR  'neural
    network*'  OR  'automatic interpretation' OR 'machine learning')  AND  TITLE-ABS-KEY
    ( electrocardiography  OR  cardiography  OR  'electrocardiographic tracing'  OR  ecg
    OR  electrocardiogram  OR  cardiogram )  AND  TITLE-ABS-KEY ( 'Intensive care unit' OR
    'cardiologic care unit'  OR  'intensive care center'  OR  'cardiologic care center' )
\

The inclusion and exclusion criteria were defined as in Table \@ref(tab:criteria).

```{r criteria, echo=FALSE}
criteria <- tribble(
  ~"Inclusion criteria", ~"Exclusion criteria",
  "ECG automatic interpretation", "Manual interpretation",
  "ECG anomaly detection", "Publication older than ten years",
  "ECG context change detection", "Do not attempt to identify life-threatening arrhythmias, namely asystole, extreme bradycardia, extreme tachycardia, ventricular tachycardia, and ventricular flutter/fibrillation",
  "Online Stream ECG analysis", "No performance measurements reported",
  "Specific diagnosis (like a flutter, hyperkalemia, etc.)", ""
)

kbl(criteria,
  booktabs = TRUE,
  caption = "Literature review criteria.",
  align = "ll",
  position = "ht",
  linesep = "\\addlinespace"
) |>
  row_spec(0, bold = TRUE) |>
  kable_styling(full_width = TRUE)
```


The survey is being conducted with peer review; all articles on full-text phase were obtained and
assessed for the extraction phase, except 5 articles that were not available. Due to external
factors, the survey is currently stalled in the Data Extraction phase.

Fig. \@ref(fig:prisma) shows the flow diagram of the resulting screening using PRISMA format.

```{r prisma, echo=FALSE, out.width="70%", fig.cap="Flowchart of the literature survey."}
my_graphics("PRISMA", "figure")
```

The peer review is being conducted by the author of this thesis and another colleague, Dr. Andrew
Van Benschoten, from the Matrix Profile Foundation [@mpf2020].

Table. \@ref(tab:kappa) shows the Inter-rater Reliability (IRR) of the screening phases, using
Cohen's $\kappa$ statistic. The bottom line shows the estimated accuracy after correcting
possible confounders [@Bakeman2011].

```{r kappa, echo=FALSE, eval=!knitr::is_latex_output()}
my_kable(
  title = "Inter-rater Reliability on the literature survey process.",
  label = "kappa",
  content = r"(<thead> <tr> <th class="tg-top" colspan="2"> </th> <th class="tg-top"
colspan="2"> Title-Abstract<br>(2388 articles) </th> <th class="tg-top"> </th> <th class="tg-top"
colspan="2"> Full-Review<br>(303 articles) </th> </tr></thead> <tbody> <tr> <td
colspan="2"> </td><td class="tg-second-top" colspan="2"> Reviewer #2 </td><td></td><td
class="tg-second-top" colspan="2"> Reviewer #2 </td></tr><tr> <td colspan="2"> </td><td
class="tg-second-top"> Include </td><td class="tg-second-top"> Exclude </td><td></td><td
class="tg-second-top"> Include </td><td class="tg-second-top"> Exclude </td></tr><tr> <td class="tg-cross-top-low"
rowspan="2"> Reviewer #1 </td><td class="tg-cross-top"> Include </td><td class="tg-cross-top"> 185 </td><td
class="tg-cross-top"> 381 </td><td class="tg-cross-top"></td><td class="tg-cross-top"> 63 </td><td
class="tg-cross-top"> 58 </td></tr><tr> <td class="tg-cross-low"> Exclude </td><td class="tg-cross-low"> 129
</td><td class="tg-cross-low"> 1693 </td><td class="tg-cross-low"> </td><td class="tg-cross-low"> 13
</td><td class="tg-cross-low"> 169 </td></tr><tr> <td class="tg-body-left" colspan="2"> Cohen’s omnibus
<span class="math inline">\(\kappa\)</span> </td><td class="tg-body-center" colspan="2"> 0.30 </td><td
class="tg-body-center"> </td><td class="tg-body-center" colspan="2"> 0.48 </td></tr><tr> <td class="tg-body-left"
colspan="2"> Maximum possible <span class="math inline">\(\kappa\)</span> </td><td class="tg-body-center"
colspan="2"> 0.66 </td><td class="tg-body-center"> </td><td class="tg-body-center" colspan="2"> 0.67
</td></tr><tr> <td class="tg-body-left" colspan="2"> Std Err for <span class="math
inline">\(\kappa\)</span> </td><td class="tg-body-center" colspan="2"> 0.02 </td><td class="tg-body-center">
</td><td class="tg-body-center" colspan="2"> 0.05 </td></tr><tr> <td class="tg-body-left" colspan="2"> Observed
Agreement </td><td class="tg-body-center" colspan="2"> 79% </td><td class="tg-body-center"> </td><td
class="tg-body-center" colspan="2"> 77% </td></tr><tr> <td class="tg-body-left" colspan="2"> Random Agreement
</td><td class="tg-body-center" colspan="2"> 69% </td><td class="tg-body-center"> </td><td class="tg-body-center"
colspan="2"> 55% </td></tr><tr> <td class="tg-bottom-left" colspan="2"> Agreement corrected with KappaAcc
</td><td class="tg-bottom" colspan="2"> 82% </td><td class="tg-bottom"> </td><td class="tg-bottom"
colspan="2"> 85% </td></tr></tbody>)"
)
```
\
```{r kappaa, echo = FALSE, eval=knitr::is_latex_output(), results="asis"}
cat(r"(
\begin{table}[ht]
\caption{\label{tab:kappa}Inter-rater Reliability on the literature survey process.}
\centering
\begin{tabular}{llcclcc}
\toprule
                                                  &                             & \multicolumn{2}{c}{\textbf{\begin{tabular}[c]{@{}c@{}}Title-Abstract\\ (2388 articles)\end{tabular}}} & \textbf{} & \multicolumn{2}{c}{\textbf{\begin{tabular}[c]{@{}c@{}}Full-Review\\ (303 articles)\end{tabular}}} \\ \cline{3-4} \cline{6-7}
                                                  &                             & \multicolumn{2}{c}{Reviewer \#2}                                                                      &           & \multicolumn{2}{c}{Reviewer \#2}                                                                  \\ \cline{3-4} \cline{6-7}
                                                  &                             & \multicolumn{1}{l}{Include}                       & \multicolumn{1}{l}{Exclude}                       &           & \multicolumn{1}{l}{Include}                     & \multicolumn{1}{l}{Exclude}                     \\ \hline
\multicolumn{1}{r}{\multirow{2}{*}{Reviewer \#1}} & \multicolumn{1}{r}{Include} & 185                                               & 381                                               &           & 63                                              & 58                                              \\
\multicolumn{1}{r}{}                              & \multicolumn{1}{r}{Exclude} & 129                                               & 1693                                              &           & 13                                              & 169                                             \\ \hline
Cohen's omnibus  $\kappa$                         &                             & \multicolumn{2}{c}{0.30}                                                                              &           & \multicolumn{2}{c}{0.48}                                                                          \\
Maximum possible $\kappa$                         &                             & \multicolumn{2}{c}{0.66}                                                                              &           & \multicolumn{2}{c}{0.67}                                                                          \\
Std Err for $\kappa$                              &                             & \multicolumn{2}{c}{0.02}                                                                              &           & \multicolumn{2}{c}{0.05}                                                                          \\
Observed Agreement                                &                             & \multicolumn{2}{c}{79\%}                                                                              &           & \multicolumn{2}{c}{77\%}                                                                          \\
Random Agreement                                  &                             & \multicolumn{2}{c}{69\%}                                                                              &           & \multicolumn{2}{c}{55\%}                                                                          \\ \hline\addlinespace
\multicolumn{2}{l}{\textbf{Agreement corrected with KappaAcc}}          & \multicolumn{2}{c}{\textbf{82\%}}                                                                     & \textbf{} & \multicolumn{2}{c}{\textbf{85\%}}                                                                 \\ \bottomrule
\end{tabular}
\end{table}
)")
```

The purpose of using Cohen's $\kappa$ in such a review is to allow us to gauge the agreement of both
reviewers on the task of selecting the articles according to the goal of the survey. The most naive
way to verify this would be simply to measure the overall agreement (the number of articles included
and excluded by both, divided by the total number of articles). Nevertheless, this would not take
into account the agreement we could expect purely by chance.

However, the $\kappa$ statistic must be assessed carefully. This topic is beyond the scope of this work
therefore it will be explained briefly.

While it is widely used, the $\kappa$ statistic is well criticized. The direct interpretation of
its value depends on several assumptions that are often violated. (1) It is assumed that both
reviewers have the same level of experience; (2) The "codes" (include, exclude) are identified with
same accuracy; (3) The "codes" prevalences are the same; (4) There is no reviewer bias towards one of
the choices [@Sim2005; @Bakeman1997].

In addition, the number of "codes" affects the relation between the value of $\kappa$ and the actual
agreement between the reviewers. For example, given equiprobable "codes" and reviewers who are 85%
accurate, the value of $\kappa$ are 0.49, 0.60, 0.66, and 0.69 when number of codes is 2, 3, 5, and
10, respectively [@Bakeman1997; @Morgan2019].

To take into account these limitations, the agreement between reviewers was calculated using the
KappaAcc [@Bakeman2011] from Professor Emeritus Roger Bakeman, Georgia State University, which
computes the estimated accuracy of simulated reviewers.

### RAW data

To better understand the data acquisition, it has been acquired a Single Lead Heart Rate Monitor
breakout from Sparkfun™ [@sparkfun2021] using the AD8232 [@AnalogDevices2020] microchip from Analog
Devices Inc., compatible with Arduino^®^ [@arduino2021], for an in-house experiment (Fig.
\@ref(fig:ad8232)).

```{r ad8232, echo=FALSE, out.width="40%", fig.show="hold", fig.cap="Single Lead Heart Rate Monitor."}
knitr::include_graphics(c("figure/sparkfun.jpg", "figure/FullSetup.jpg"))
```

The output gives us a RAW signal, as shown in Fig. \@ref(fig:rawsignal).

```{r rawsignal, echo=FALSE, out.width="50%", fig.cap="RAW output from Arduino at ~300hz."}
my_graphics("arduino_plot", "figure")
```

After applying the same settings as the Physionet database (collecting the data at 500hz, resample
to 250hz, pass-filter, and notch filter), the signal is much better, as shown in Fig.
\@ref(fig:filtersignal).

```{r filtersignal, echo=FALSE, out.width="90%", fig.cap="Gray is RAW, Red is filtered."}
my_graphics("filtered_ecg", "figure")
```

### Preparing the data

Usually, data obtained by sensors needs to be "cleaned" for proper evaluation. That is different
from the initial filtering process where the purpose is to enhance the signal. Here we are dealing
with artifacts, disconnected cables, wandering baselines and others.

Several SQIs (Signal Quality Indexes) are used in the literature [@eerikainen2015], some trivial
measures as _kurtosis_, _skewness_, median local noise level, other more complex as pcaSQI (the
ratio of the sum of the five largest eigenvalues associated with the principal components over the
sum of all eigenvalues obtained by principal component analysis applied to the time aligned ECG
segments in the window). An assessment of several different methods to estimate electrocardiogram
signal quality can was performed by Del Rio, *et al* [@DelRio2011].

 By experimentation (yet to be validated), a simple formula gives us the "complexity" of the signal
 and correlates well with the noisy data is shown in Equation $\eqref{complex}$ [@Batista2014].

\
$$
\sqrt{\sum_{i=1}^w((x_{i+1}-x_i)^2)}, \quad \text{where}\; w \; \text{is the window size} \tag{2} \label{complex}
$$
\

Fig. \@ref(fig:sqi) shows some SQIs and their relation with the data.

```{r sqi, echo=FALSE, out.width="100%", fig.cap="Green line is the \"complexity\" of the signal."}
my_graphics("noise", "figure")
```

```{r createfilter, include=FALSE}
source(here("scripts", "common", "read_ecg.R"))
source(here("scripts", "common", "sqi.R"))
filter_w <- 200
limit <- 8
file <- "a104s"
size_w <- 16
size_h <- 5

if (file.exists(here("output/createfilter.rds"))) {
  data <- readRDS(here("output/createfilter.rds"))
} else {
  cli::cli_abort("File not found.")
  data <- read_ecg_csv(here(glue("inst/extdata/physionet/{file}.hea")))
  data <- data[[file]]$II
  saveRDS(data, here("output/createfilter.rds"), compress = "xz")
}
norm_data <- tsmp:::znorm(data)
filter <- win_complex(norm_data, filter_w)
filter <- filter > limit

if (knitr::is_latex_output()) {
  grDevices::pdf(here("docs/figure/regime_filter.pdf"),
    width = size_w, height = size_h
  )
} else {
  svglite::svglite(here("docs/figure/regime_filter.svg"),
    width = size_w, height = size_h
  )
}
plot(norm_data, main = "", type = "l", ylab = "", xlab = "index", lwd = 0.2)
points(cbind(which(filter), 0), col = "blue", pch = 19)
dev.off()
```

Fig. \@ref(fig:datafilter) shows that noisy data (probably patient muscle movements) are marked
with a blue point and thus are ignored by the algorithm.

```{r datafilter, echo=FALSE, out.width="100%", fig.cap="Noisy data marked by the \"complexity\" filter."}
my_graphics("regime_filter", "figure")
```

Although this step of "cleaning" the data is often used, this step will also be tested if it is
really necessary, and the performance with and without "cleaning" will be reported.

### Detecting regime changes

The regime change approach will be using the _Arc Counts_ concept, used on the FLUSS (Fast Low-cost
Unipotent Semantic Segmentation) algorithm, as explained by Gharghabi, _et al._,[@gharghabi2018].

The FLUSS (and FLOSS, the online version) algorithm is built on top of the Matrix Profile
(MP)[@Yeh2017a], described on section \@ref(matrixprofile). Recalling that the MP and the
companion Profile Index (PI) are two vectors holding information about the 1-NN. One can imagine
several "arcs" starting from one "index" to another. This algorithm is based on the assumption that
between two regimes, the most similar shape (its nearest neighbor) is located on "the same side", so
the number of "arcs" decreases when there is a change in the regime and increases again. As show on
Fig. \@ref(fig:arcsoriginal). This drop on the _Arc Counts_ is a signal that a change in the
shape of the signal has happened.

```{r arcsoriginal, echo=FALSE, out.width="100%", fig.cap="FLUSS algorithm, using arc counts."}
my_graphics("fluss_arcs", "figure")
```

The choice of the FLOSS algorithm (the online version of FLUSS) is founded on the following arguments:

-  **Domain Agnosticism:** the algorithm makes no assumptions about the data as opposed to most
   available algorithms to date.
-  **Streaming:** the algorithm can provide real-time information.
-  **Real-World Data Suitability:** the objective is not to _explain_ all the data. Therefore, areas
   marked as "don't know" areas are acceptable.
-  **FLOSS is not:** a change point detection algorithm [@aminikhanghahi2016]. The interest here is
   changes in the shapes of a sequence of measurements.

Other algorithms we can cite are based on Hidden Markov Models (HMM) that require at least two
parameters to be set by domain experts: cardinality and dimensionality reduction. The most
attractive alternative could be the Autoplait [@Matsubara2014], which is also domain agnostic and
parameter-free. It segments the time series using Minimum Description Length (MDL) and recursively
tests if the region is best modeled by one or two HMM. However, Autoplait is designed for batch
operation, not streaming, and also requires discrete data. FLOSS was demonstrated to be superior in
several datasets in its original paper. In addition, FLOSS is robust to several changes in data like
downsampling, bit depth reduction, baseline wandering, noise, smoothing, and even deleting 3% of the
data and filling with simple interpolation. Finally, the most important, the algorithm is light and
suitable for low-power devices.

In the MP domain, it is worth also mentioning another possible algorithm: the Time Series Snippets
[@Imani2018], based on MPdist [@gharghabi2018b]. The MPdist measures the distance between two
sequences, considering how many similar sub-sequences they share, no matter the matching order. It
proved to be a useful measure (not a metric) for meaningfully clustering similar sequences. Time
Series Snippets exploits MPdist properties to summarize a dataset extracting the $k$ sequences
representing most of the data. The final result seems to be an alternative for detecting regime
changes, but it is not. The purpose of this algorithm is to find which pattern(s) explains most of
the dataset. Also, it is not suitable for streaming data. Lastly, MPdist is quite expensive compared
to the trivial Euclidean distance.

The regime change detection will be evaluated following the criteria explained in section
\@ref(evaluation).

### Classification of the new regime {#classregime}

The next step towards the objective of this work is to verify if the new regime detected by the
previous step is indeed a life-threatening pattern that we should trigger the alarm.

First, let's dismiss some apparent solutions: (1) Clustering. It is well understood that we cannot
cluster time series subsequences meaningfully with any distance measure or with any algorithm
[@Keogh2005]. The main argument is that in a meaningful algorithm, the output depends on the input,
and this has been proven to not happen in time series subsequence clustering [@Keogh2005]. (2)
Anomaly detection. In this work, we are not looking for surprises but for patterns that are known to
be life-threatening. (3) Forecasting. We may be tempted to make predictions, but this clearly is not
the idea here.

The method of choice is classification. The simplest algorithm could be a `TRUE`/`FALSE` binary
classification. Nevertheless, the five life-threatening patterns have well-defined characteristics
that may seem more plausible to classify the new regime using some kind of ensemble of binary
classifiers or a "six-class" classifier (the sixth class being the `FALSE` class).

Since the model doesn't know which life-threatening pattern will be present in the regime (or if it
will be a `FALSE` case), the model will need to check for all five `TRUE` cases, and if none of
these cases are identified, it will classify the regime as `FALSE`.

To avoid exceeding processor capacity, an initial set of shapelets [@Rakthanmanon2013] can be
sufficient to build the `TRUE`/`FALSE` classifier. And to build such a set of shapelets, leveraging
on the MP, we will use the Contrast Profile [@Mercer2021].

The Contrast Profile (CP) looks for patterns that are at the same time very *similar* to its
neighbors in class *A* while is very *different* from the nearest neighbor from class *B*. In other
words, this means that such a pattern represents well class *A* and may be taken as a "signature" of
that class.

In this case, we need to compute two MP, one self-join MP using the *positive* class $MP^{(++)}$
(the class that has the signature we want to find) and one AB-join MP using the *positive* and
*negative* classes $MP^{(+-)}$. Then we subtract the first $MP^{(++)}$ from the last $MP^{(+-)}$,
resulting in the $CP$. The high values on $CP$ are the locations for the signature candidates we
look for (the author of CP calls these segments *Plato's*).

Due to the nature of this approach, the MP's (containing values in Euclidean Distance) are truncated
for values above $\sqrt{2w}$, where $w$ is the window size. This is because values above this
threshold are negatively correlated in the Pearson Correlation space. Finally, we normalize the
values by $\sqrt{2w}$. The formula $\eqref{contrast}$ synthesizes this computation.

\
$$
CP_w = \frac{MP_{w}^{(+-)} - MP_{w}^{(++)}}{\sqrt{2w}} \quad \text{where}\; w \; \text{is the window size} \tag{3} \label{contrast}
$$
\

For a more complete understanding of the process, Fig. \@ref(fig:contrast) shows a practical example
from the original article [@Mercer2021].

\
```{r contrast, echo=FALSE, out.width="100%"}
#| fig.cap =  "Top to bottom: two weakly-labeled snippets of a larger time series. T(-) contains
#|  only normal beats. T(+) also contains PVC (premature ventricular contractions).
#|  Next, two Matrix Profiles with window size 91; AB-join is in red and self-join in blue.
#|  Bottom, the Contrast Profile showing the highest location."

my_graphics("contrast", "figure")
```


After extracting candidates for each class signature, a classification algorithm will be fitted and
evaluated using the criteria explained on section \@ref(evaluation).

### Summary of the methodology {#methodology}

To summarize the steps taken on this thesis to accomplish the main objective, Figs.
\@ref(fig:regimedetection), \@ref(fig:shapelets) and \@ref(fig:fullmodel) show the
overview of the processes involved.

First, let us introduce the concept of Nested Resampling [@Bischl2012]. It is known that when
increasing model complexity, overfitting on the training set becomes more likely to happen
[@Hastie2009]. This is an issue that this work has to countermeasure as many steps require parameter
tuning, even for almost parameter-free algorithms, like the MP.

The rule that must be followed is simple: *do not* evaluate a model on the same resampling split
used to perform its own parameter tuning. Using simple cross-validation, the information about the
test set "leaks" into the evaluation, leading to overfitting/overtuning, and gives us an optimistic
biased estimate of the performance. Bernd Bischl, 2012 [@Bischl2012] describes more deeply these
factors and also provides us with a countermeasure for that: (1) from preprocessing the data to
model selection use the training set; (2) the test set should be touched once, on the evaluation
step; (3) repeat. This guarantees that a "new" separated data is only used *after* the model is
trained/tuned.

Fig. \@ref(fig:nestedresampling) shows us this principle. The steps (1) and (2) described above
are part of the **Outer resampling**, which in each loop splits the data into two sets: the training
set and the test set. The training set is then used in the **Inner resampling** where, for example,
the usual cross-validation may be used (creating an *Analysis set* and an *Assessment set* to avoid
terminology conflict), and the best model/parameters are selected. Then, this best model is
evaluated against the unseen test set created for this resampling.

The resulting (aggregated) performance of all outer samples gives us a more honest estimative
of the expected performance on new data.

```{r nestedresampling, echo=FALSE, out.width="70%"}
#| fig.cap =  "Nested resampling.
#|  The full dataset is resampled several times (outer resampling), so each branch has its
#|  own Test set (yellow). On each branch, the Training set is used as if it were a full dataset,
#|  being resampled again (inner resampling); here, the Assessment set (blue) is used to test the
#|  learning model and tune parameters. The best model is finally evaluated on its own
#|  Test set."

my_graphics("draw-nested-resampling", "figure")
```
\

After the understanding of the Nested Resampling [@Bischl2012], the following flowcharts can be
better interpreted. Fig. \@ref(fig:regimedetection) starts with the "Full Dataset" that contains
all time series from the dataset described in section \@ref(the-data). Each time series
represents one file from the database and represents one patient.

The regime change detection will use subsampling (bootstrapping can lead to substantial bias toward
more complex models) in the Outer resampling and cross-validation in the Inner resampling. How the
evaluation will be performed and why the use of cross-validation will be explained in section
\@ref(evaluation).

```{r regimedetection, echo=FALSE, out.width="90%"}
#| fig.cap = "Pipeline for regime change detection.
#|  The full dataset (containing several patients) is divided into a Training set and a Test set.
#|  The Training set is then resampled in an Analysis set and an Assessment set. The former is
#|  used for training/parameter tuning and the latter for assessing the result. The best parameters
#|  are then used for evaluation on the Test set. This may be repeated several times."

my_graphics("draw-regime-model", "figure")
```

Fig. \@ref(fig:shapelets) shows the processes for training the classification model. First, the
last ten seconds of each time series will be identified (the event occurs in this segment). Then the
dataset will be grouped by class (type of event) and `TRUE`/`FALSE` (alarm), so the Outer/Inner
resampling will produce a Training/Analysis set and Test/Assessment set with similar frequency to
the full dataset.

The next step will be to extract shapelet candidates using the Contrast Profile and train the
classifier.

This pipeline will use subsampling (for the same reason above) in the Outer resampling and
cross-validation in the Inner resampling. How the evaluation will be performed and why the use of
cross-validation will be explained in section \@ref(evaluation).

```{r shapelets, echo=FALSE, out.width="60%"}
#| fig.cap = "Pipeline for alarm classification.
#|  The full dataset (containing several patients) is grouped by class and by TRUE/FALSE alarm.
#|  This grouping allows resampling to keep a similar frequency of classes and TRUE/FALSE of the full dataset.
#|  Then the full dataset is divided on a Training set and a Test set.
#|  The Training set is then resampled in an Analysis set and an Assessment set. The former is
#|  used for extracting shapelets, training the model and parameter tuning; the latter for assessing
#|  the performance of the model. Finally, the best model is evaluated on the Test set.
#|  This may be repeated several times."
my_graphics("draw-classif-model", "figure")
```

Finally, Fig. \@ref(fig:fullmodel) shows how the final model will be used on the field. In a
streaming scenario, the data will be collected and processed in real-time to maintain an up to date
Matrix Profile. The FLOSS algorithm will be looking for a regime change. When a regime change is
detected, a sample of this new regime will be presented to the trained classifier that will evaluate
if this new regime is a life-threatening condition or not.

```{r fullmodel, echo=FALSE, out.width="60%"}
#| fig.cap = "Pipeline of the final process.
#|  The streaming data, coming from one patient, is processed to create its Matrix Profile.
#|  Then, the FLOSS algorithm is computed for detecting a regime change. When a new regime is
#|  detected, a sample of this new regime is analysed by the model and a decision is made. If
#|  the new regime is life-threatening, the alarm will be fired."
my_graphics("draw-global-model", "figure")
```


## Evaluation of the algorithms {#evaluation}

The subsampling method used on both algorithms, regime change, and classification, will be
the Cross-Validation, as the learning task will be in batches.

Other options dismissed [@Bischl2012]:

* Leave-One-Out Cross-Validation: has better properties for regression than for classification. It
  has a high variance as an estimator of the mean loss. It also is asymptotically inconsistent and
  tends to select too complex models. It is demonstrated empirically that 10-fold CV is often
  superior.

* Bootstrapping: while it has low variance, it may be optimistic-biased on more complex models.
  Also, its resampling method with replacement can leak information into the assessment set.

* Subsampling: is like bootstrapping, but without replacement. The only argument for not choosing it
  is that we make sure all the data is used for analysis and assessment with Cross-Validation.


### Regime change

A detailed discussion about the evaluation process of segmentation algorithms is made by the
FLUSS/FLOSS author [@gharghabi2018]. Previous researches have used precision/recall or derived
measures for performance. The main issue is how to assume that the algorithm was correct? Is this a
miss if the ground truth says the change occurred at location 10,000, and the algorithm detects a
change at location 10,001?

As pointed out by the author, several independent researchers have suggested a temporal tolerance,
that solves one issue but has a hard time penalizing any tiny miss beyond this tolerance.

The second issue is an over-penalization of an algorithm in which most of the detections are good,
but just one (or a few) is poor.

The author proposes the solution depicted in Fig. \@ref(fig:flosseval). It gives 0 as the best
score and 1 as the worst. The function sums the distances between the ground truth locations and the
locations suggested by the algorithm. The sum is then divided by
<!--the product of the number of segments, and then -->
the length of the time series to normalize the range to [0, 1].

The goal is to minimize this score.

<!--
TODO: review the problem when there are too many detections
-->

```{r flosseval, echo=FALSE, out.width="100%"}
#| fig.cap="Regime change evaluation. The top line illustrates the ground truth, and the
#|  bottom line the locations reported by the algorithm. Note that multiple proposed locations
#|  can be mapped to a single ground truth point."
my_graphics("floss_eval", "figure")
```

### Classification {#classification-criteria}

As described in section \@ref(classregime), the model for classification will use a set of shapelets
to identify if we have a `TRUE` (life-threatening) regime or a `FALSE` (non-life-threatening) regime.

Although the implementation of the final process will be using streaming data, the classification
algorithm will work in batches because it will not be applied on every single data point but on
samples that are extracted when a regime change is detected. During the training phase, the data is
also analyzed in batches.

One important factor we must consider is that, in real-world, the majority of regime changes will be
`FALSE` (i.e., not life-threatening). Thus, a performance measure that is robust to class imbalance
is needed if we want to assess the model after it was trained on the field.

It is well known that the *Accuracy* measure is not reliable for unbalanced data
[@Akosa2017; @Bekkar2013] as it returns optimistic results for a classifier on the majority class. A
description of common measures used on classification is available [@Akosa2017; @Chicco2020]. Here
we will focus on three candidate measures that can be used: F-score (well discussed on
[@Chicco2020]), Matthew's Correlation Coefficient (MCC) [@Matthews1975] and $\kappa_m$ statistic
[@Bifet2015].

The F-score (let's abbreviate to F~1~ as this is the more common setting) is widely used on
*information retrieval*, where the classes are usually classified as "relevant" and "irrelevant",
and combines the *recall* (also known as sensitivity) and the *precision* (the positive predicted
value). *Recall* assess how well the algorithm retrieves relevant examples among the (usually few)
relevant items in the dataset. In contrast, *precision* assesses the proportion of indeed relevant
items which are contained in the retrieved examples. It ranges from [0, 1]. It completely ignores
the irrelevant items that were not retrieved (usually, this set contains many items). In
classification tasks, its main weakness is not evaluating the True Negatives. If the proportion of a
random classifier gets towards the `TRUE` class (increasing the False Positives significantly), this
score actually gets better, thus not suitable to our case. The F~1~ score is defined on Equation
$\eqref{score}$.

$$
F_1 score = \frac{2 \cdot TP}{2 \cdot TP + FP + FN} = 2 \cdot \frac{precision \cdot recall}{precision + recall} \tag{4} \label{fscore}
$$

The MCC is a good alternative to the F~1~ when we do care about the True negatives (both were
considered to "provide more realistic estimates of real-world model performance" [@Dubey2018]). It
is a method to compute the *Pearson product-moment correlation coefficient* [@Delgado2019] between
the actual and predicted values. It ranges from [-1, 1]. The MCC is the only binary classification
rate that only gives a high score if the binary classifier correctly classified the majority of the
positive and negative instances [@Chicco2020]. One may argue that Cohen's $\kappa$ has the same
behavior. Still, there are two main differences (1) MCC is *undefined* in the case of a *majority
voter*. At the same time, Cohen's $\kappa$ doesn't discriminate this case from the random classifier
($\kappa$ is zero for both cases). (2) It is proven that in a special case when the classifier is
increasing the False Negatives, Cohen's $\kappa$ doesn't get worse as expected, MCC doesn't have
this issue [@Delgado2019]. MCC is defined on equation $\eqref{mccval}$.

$$
MCC = \frac{TP \cdot TN - FP \cdot FN}{\sqrt{(TP + FP) \cdot (TP + FN) \cdot (TN + FP) \cdot (TN + FN)}} \tag{5} \label{mccval}
$$

The $\kappa_m$ statistic [@Bifet2015] is a measure that considers not the *random classifier* but
the *majority voter* (a classifier that only votes on the larger class). It was introduced by Bifet
*et al.* [@Bifet2015] for being used in online settings, where the class balance may change over
time. It is defined on Equation $\eqref{kappam}$, where $p_0$ is the observed accuracy, and $p_m$ is
the accuracy of the majority voter. Theoretically, the score ranges from ($-\infty$, 1]. Still, in
practice, you see negative numbers if the classifier performs worse than the majority voter and
positive numbers if performing better than the majority number until the maximum of 1 when the
classifier is optimal.

$$
\kappa_m = \frac{p_0 - p_m}{1 - p_m} \tag{6} \label{kappam}
$$

In the inner resampling (model training/tuning), the classification will be binary, and in our case,
we know that the data is slightly unbalanced (60% false alarms). For this step, the metric for model
selection will be the MCC. Nevertheless, during the optimization process, the algorithm will seek to
minimize the False Negative Rate ($FNR = \frac{FN}{TP+FN}$), and between ties, the smaller FNR wins.

In the outer resampling, the MCC and $\kappa_m$ of all winning models will aggregate and report
using the median and interquartile range.

For different classifiers, we will use Wilcoxon's signed-rank test for comparing their performances.
This method is known to have low Type I and Type II errors in this kind of comparison [@Bifet2015].

### Full model (streaming setting)

For the final assessment, the best and the average model of the previous pipelines will be assembled
and tested using the whole original dataset.

The algorithm will be tested in each of the five life-threatening events split individually in order
to evaluate its strengths and weakness.

For more transparency, the whole confusion matrix will be reported, as well as the MCC, $\kappa_m$, and
the FLOSS evaluation.

# Current results

## Regime change detection

As we have seen previously, the FLOSS algorithm is built on top of the Matrix Profile (MP). Thus,
we have proposed several parameters that may or not impact the FLOSS prediction performance.

The variables for building the MP are:

-   **`mp_threshold`**: the minimum similarity value to be considered for 1-NN.
-   **`time_constraint`**: the maximum distance to look for the nearest neighbor.
-   **`window_size`**: the default parameter always used to build an MP.

Later, the FLOSS algorithm also has a parameter that needs tuning to optimize the prediction:

-   **`regime_threshold`**: the threshold below which a regime change is considered.
-   **`regime_landmark`**: the point in time where the regime threshold is applied.

Using the `tidymodels` framework, we performed a basic grid search on all these parameters.

Fig. \@ref(fig:thepipeline) shows the workflow using Nested resamplig as described on section \@ref(methodology).
Fig. \@ref(fig:flossregime) shows an example of the regime change detection pipeline. The graph on top shows
 the ECG streaming; the blue line marks the ten seconds before the original alarm was fired; the red line
 marks the time constraint used on the example; the dark red line marks the limit for taking a decision
 in this case of Asystole; The blue horizontal line represents the size of the sliding window. The graph
  on the middle shows the Arc counts as seen by the algorithm (with the corrected distribution); the red
  line marks the current minimum value and its index; the blue horizontal line shows the minimum value
  seen until then. The graph on the bottom shows the computed Arc counts (raw) and the red line is the
   theoretical distribution used for correction.

```{r floss_cached, echo=FALSE, cache=FALSE}
network <- readRDS(here::here("output", "regime_network.rds"))
net <- network |>
  visNetwork::visPhysics(hierarchicalRepulsion = list(
    springLength = 1,
    avoidOverlap = 0.5,
    nodeDistance = 120
  ))
```

```{r thepipeline, out.width="90%", fig.cap="FLOSS pipeline."}
if (knitr::is_latex_output()) {
  my_graphics("the-network")
} else {
  visNetwork::visInteraction(net, hover = TRUE, multiselect = TRUE, tooltipDelay = 100)
}
```

```{r flossregime, echo=FALSE, out.height="80%", out.width="80%"}
#| fig.cap="Regime change detection example.
#|  The graph on top shows the ECG streaming; the blue line marks the ten seconds
#|  before the original alarm was fired; the red line marks the time constraint of 1250;
#|  the dark red line marks the limit for taking a decision in this case of Asystole
#|  the blue horizontal line represents the size of the sliding window.
#|  The graph on the middle shows the Arc counts as seen by the algorithm (with the corrected
#|  distribution); the red line marks the current minimum value and its index; the blue
#|  horizontal line shows the minimum value seen until then.
#|  The graph on the bottom shows the computed Arc counts (raw), and the red line is the
#|  theoretical distribution used for correction."
my_graphics("floss_regime", "figure")
```

The dataset used for working with the Regime Change algorithm was the "Paroxysmal Atrial Fibrillation Events Detection from Dynamic ECG Recordings: The 4th China Physiological Signal Challenge 2021" hosted by Zenodo [@bischoff2021afib] under the same license as Physionet.

The selected records were those that contain paroxysmal atrial fibrillation events, a total of 229 records. The records were split in a proportion of 3/4 for the training set (inner resampling) and 1/4 for the test set (outer resampling). The inner resampling was performed using a 5-fold cross-validation, which accounts for 137 records for fitting the models and 92 records for assessing them in the inner resampling.

The following parameters were used:

-   The MP parameters were explored using the following values:
    -   **`mp_threshold`**: 0.0 to 0.9, by 0.1 steps;
    -   **`time_constraint`**: 0, 800 and 1500;
    -   **`window_size`**: 25 to 350, by 25 steps;
-   The FLOSS parameters were explored using the following values:
    -   **`regime_threshold`**: 0.05 to 0.90, by 0.05 steps;
    -   **`regime_landmark`**: 1 to 10, by 0.5 steps.

### Parameters analysis

The above process was an example of parameter tuning seeking the best model for a given set of parameters. It used a nested cross-validation procedure that aims to find the best combination of parameters and avoid overfitting.

While this process is powerful and robust, it does not show us the importance of each parameter. At least one parameter has been introduced by reasoning about the problem (`mp_threshold`), but how important it (and other parameters) is for predicting regime changes?

For example, the process above took 4 days, 20 hours, and 15 minutes to complete the grid search using an Intel(R) Xeon(R) Silver 4210R \@ 2.40 GHz server. Notice that about 133 different combinations of parameters were tested on computing the MP (not FLOSS, the `regime_threshold`), 5 folds, 2 times each. That sums up about 35.2 x 10^9^ all-pairs Euclidean distances computed on less than 5 days (on CPU, not GPU). Not bad.

Another side note on the above process, it is not a "release" environment, so we must consider lots of overhead in computation and memory usage that must be taken into account during these five days of grid search. Thus, much time can be saved if we know what parameters are essential for the problem.

In order to check the effect of the parameters on the model, we need to compute the *importance* of each parameter.

Wei *et al.* published a comprehensive review on variable importance analysis [@Wei2015].

Our case is not a typical case of variable importance analysis, where a set of *features* are tested against an *outcome*. Instead, we have to proxy our analysis by using as *outcome* the FLOSS performance score and as *features* (or *predictors*) the tuning parameters that lead to that score.

That is accomplished by fitting a model using the tuning parameters to predict the FLOSS score and then applying the techniques to compute the importance of each parameter.

For this matter, a Bayesian Additive Regression Trees (BART) model was chosen after an experimental trial with a set of regression models (including glmnet, gbm, mlp) and for its inherent characteristics, which allows being used for model-free variable selection [@Chipman2010]. The best BART model was selected using 10-fold cross-validation repeated 3 times, having great predictive power with an RMSE around 0.2 and an R^2^ around 0.99. With this fitted model, we could evaluate each parameter's importance.

### Interactions

Before starting the parameter importance analysis, we need to consider the parameter interactions since this is usually the weak spot of the analysis techniques.

The first BART model was fitted using the following parameters:

$$
\begin{aligned}
E( score ) &= \alpha + time\_constraint\\
 &\quad + mp\_threshold + window\_size\\
 &\quad + regime\_threshold + regime\_landmark
\end{aligned} \tag{7} \label{eq-first}
$$

After checking the interactions, this is the refitted model:

$$
\begin{aligned}
E( score ) &= \alpha + time\_constraint\\
&\quad + mp\_threshold + window\_size\\
&\quad + regime\_threshold + regime\_landmark\\
&\quad + \left(regime\_threshold \times regime\_landmark\right)\\
&\quad + \left(mp\_threshold \times regime\_landmark\right)\\
&\quad + \left(mp\_threshold \times window\_size\right)
\end{aligned} \tag{8} \label{eq-fitted}
$$

Fig. \@ref(fig:interaction) shows the variable interaction strength between pairs of variables. That allows us to verify if there are any significant interactions between the variables. Using the information from the first model fit, equation $\eqref{eq-first}$, we see that `regime_threshold` interacts strongly with `regime_landmark`. This interaction was already expected, and we see that even after refitting the model, equation $\eqref{eq-fitted}$, this interaction is still strong.

This is not a problem *per se* but a signal we must be aware of when exploring the parameters.

```{r interaction, fig.height = 7, fig.width= 8, out.width="100%"}
#| fig.cap="Variable interactions strength using feature importance ranking measure (FIRM) approach [@Greenwell2018].
#|  A) Shows strong interaction between `regime_threshold` and `regime_landmark`, `mp_threshold` and `window_size`,
#|     `mp_threshold` and `regime_landmark`.
#|  B) Refitting the model with these interactions taken into account, the strength is substantially reduced, except
#|     for the first, showing that indeed there is a strong correlation between those variables."

interactions <- readRDS(here("output", "importances_lmk.rds"))
importance_firm <- interactions$importance_firm
importance_perm <- interactions$importance_perm
importance_shap <- interactions$importance_shap
shap_html_test <- interactions$shap_html_test
shap_fastshap_all_test <- interactions$shap_fastshap_all_test
interactions <- interactions$interactions

interactions2 <- readRDS(here("output", "importances2_lmk.rds"))
importance_firm2 <- interactions2$importance_firm2
importance_perm2 <- interactions2$importance_perm2
importance_shap2 <- interactions2$importance_shap2
shap_fastshap_all_test2 <- interactions2$shap_fastshap_all_test2
shap_html_test2 <- interactions2$shap_html_test2
interactions2 <- interactions2$interactions2

interactions_plot <- ggplot2::ggplot(interactions, ggplot2::aes(
  x = reorder(Variables, Interaction),
  y = Interaction, fill = Variables
)) +
  ggplot2::geom_col(color = "grey35", linewidth = 0.2) +
  ggplot2::coord_flip() +
  ggplot2::labs(
    title = "Normal fit",
    y = ggplot2::element_blank(),
    x = ggplot2::element_blank()
  ) +
  ggplot2::ylim(0, 1.65) +
  ggplot2::theme_bw() +
  ggplot2::theme(legend.position = "none")

interactions2_plot <- ggplot2::ggplot(interactions2, ggplot2::aes(
  x = reorder(Variables, Interaction),
  y = Interaction, fill = Variables
)) +
  ggplot2::geom_col(color = "grey35", linewidth = 0.2) +
  ggplot2::coord_flip() +
  ggplot2::labs(
    title = "Taking into account the interactions",
    y = "Interaction strength",
    x = ggplot2::element_blank()
  ) +
  ggplot2::ylim(0, 1.65) +
  ggplot2::theme_bw() +
  ggplot2::theme(legend.position = "none")

inter <- interactions_plot / interactions2_plot
inter + plot_annotation(
  title = "Variable Interaction Strength",
  tag_levels = c("A", "1"),
  theme = ggplot2::theme_bw()
)
```

### Importance

After evaluating the interactions, we then can perform the analysis of the variable importance. The goal is to understand how the FLOSS score behaves when we change the parameters.

Here is a brief overview of the different techniques:

#### Feature Importance Ranking Measure (FIRM)

The FIRM is a variance-based method. This implementation uses the ICE curves to quantify each feature effect which is more robust than partial dependance plots (PDP) [@Greenwell2020].

It is also helpful to inspect the ICE curves to uncover some heterogeneous relationships with the outcome [@Molnar2022].

**Advantages:**

-   Has a causal interpretation (for the model, not for the real world)
-   ICE curves can uncover heterogeneous relationships

**Disadvantages:**

-   The method does not take into account interactions.

#### Permutation

The Permutation method was introduced by Breiman in 2001 [@Breiman2001] for Random Forest, and the implementation used here is a model-agnostic version introduced by Fisher *et al.* in 2019 [@Fisher2018]. A feature is "unimportant" if shuffling its values leaves the model error unchanged, assuming that the model has ignored the feature for the prediction.

**Advantages:**

-   Easy interpretation: the importance is the increase in model error when the feature's information is destroyed.
-   No interactions: the interaction effects are also destroyed by permuting the feature values.

**Disadvantages:**

-   It is linked to the model error: not a disadvantage *per se*, but may lead to misinterpretation if the goal is to understand how the output varies, regardless of the model's performance. For example, if we want to measure the robustness of the model when someone tampers the features, we want to know the *model variance* explained by the features. Model variance (explained by the features) and feature importance correlate strongly when the model generalizes well (it is not overfitting).
-   Correlations: If features are correlated, the permutation feature importance can be biased by unrealistic data instances. Thus we need to be careful if there are strong correlations between features.

#### SHAP

The SHAP feature importance [@Lundberg2017] is an alternative to permutation feature importance. The difference between both is that Permutation feature importance is based on the decrease in model performance, while SHAP is based on the magnitude of feature attributions.

**Advantages:**

-   It is not linked to the model error: as the underlying concept of SHAP is the Shapley value, the value attributed to each feature is related to its contribution to the output value. If a feature is important, its addition will significantly affect the output.

**Disadvantages:**

-   Computer time: Shapley value is a computationally expensive method and usually is computed using Montecarlo simulations.
-   The Shapley value can be misinterpreted: The Shapley value of a feature value **is not** the difference of the predicted value after removing the feature from the model training. The interpretation of the Shapley value is: "Given the current set of feature values, the contribution of a feature value to the difference between the actual prediction and the mean prediction is the estimated Shapley value" [@Molnar2022].
-   Correlations: As with other permutation methods, the SHAP feature importance can be biased by unrealistic data instances when features are correlated.

### Importance analysis

Using the three techniques simultaneously allows a broad comparison of the model behavior [@Greenwell2020]. All three methods are model-agnostic (separates interpretation from the model), but as we have seen, each method has its advantages and disadvantages [@Molnar2022].

Fig. \@ref(fig:importance) then shows the variable importance using three methods: Feature Importance Ranking Measure (FIRM) using Individual Conditional Expectation (ICE), Permutation-based, and Shapley Additive explanations (SHAP). The first line of this figure shows an interesting result that probably comes from the main disadvantage of the FIRM method: *the method does not take into account interactions*. We see that FIRM is the only one that disagrees with the other two methods, giving much importance to `window_size`.

In the second line, taking into account the interactions, we see that all methods somewhat agree with each other, accentuating the importance of `regime_threshold`, which makes sense as it is the most evident parameter we need to set to determine if the *Arc Counts* are low enough to indicate a regime change.

```{r importance, fig.height = 7, fig.width= 15, out.width="100%"}
#| fig.cap="Variables importances using three different methods. A) Feature Importance Ranking Measure
#|  using ICE curves. B) Permutation method. C) SHAP (400 iterations). Line 1 refers to the original
#|  fit, and line 2 to the re-fit, taking into account the interactions between variables
#|  (@fig-interaction)."


importance_firm_plot <- ggplot2::ggplot(importance_firm, aes(
  x = reorder(Variable, Importance),
  y = Importance, fill = Variable
)) +
  ggplot2::geom_col(colour = "grey35", linewidth = 0.8, show.legend = FALSE) +
  ggplot2::coord_flip() +
  ggplot2::labs(
    title = "Feature Importance Ranking Measure",
    subtitle = "Individual Conditional Expectation",
    x = "",
    y = ggplot2::element_blank()
  ) +
  ggplot2::ylim(0, 3.5) +
  ggplot2::theme_bw() +
  ggplot2::theme(
    legend.position = "none",
    plot.margin = margin(0, 0, 0, 10)
  )

importance_perm_plot <- ggplot2::ggplot(importance_perm, aes(
  x = reorder(Variable, Importance),
  y = Importance, fill = Variable
)) +
  ggplot2::geom_boxplot(colour = "grey35", linewidth = 0.5, show.legend = FALSE) +
  ggplot2::coord_flip() +
  ggplot2::labs(
    title = "Permutation-based (100x)",
    x = "",
    y = ggplot2::element_blank()
  ) +
  ggplot2::ylim(2, 15) +
  ggplot2::theme_bw() +
  ggplot2::theme(
    legend.position = "none",
    plot.margin = margin(0, 0, 0, 10)
  )

importance_shap_plot <- ggplot2::ggplot(importance_shap, aes(
  x = reorder(Variable, Importance),
  y = Importance, fill = Variable
)) +
  ggplot2::geom_col(colour = "grey35", linewidth = 0.8, show.legend = FALSE) +
  ggplot2::coord_flip() +
  ggplot2::labs(
    title = "SHAP (400 iterations)",
    x = "",
    y = ggplot2::element_blank()
  ) +
  ggplot2::ylim(0, 1.6) +
  ggplot2::theme_bw() +
  ggplot2::theme(
    legend.position = "none",
    plot.margin = margin(0, 0, 0, 10)
  )

importance_firm2_plot <- ggplot2::ggplot(importance_firm2, aes(
  x = reorder(Variable, Importance),
  y = Importance, fill = Variable
)) +
  ggplot2::geom_col(colour = "grey35", linewidth = 0.8, show.legend = FALSE) +
  ggplot2::coord_flip() +
  ggplot2::labs(
    x = "",
    y = "Importance"
  ) +
  ggplot2::ylim(0, 3.5) +
  ggplot2::theme_bw() +
  ggplot2::theme(
    legend.position = "none",
    plot.margin = margin(0, 0, 0, 10)
  )

importance_perm2_plot <- ggplot2::ggplot(importance_perm2, aes(
  x = reorder(Variable, Importance),
  y = Importance, fill = Variable
)) +
  ggplot2::geom_boxplot(colour = "grey35", linewidth = 0.5, show.legend = FALSE) +
  ggplot2::coord_flip() +
  ggplot2::labs(
    x = "",
    y = "Importance"
  ) +
  ggplot2::ylim(2, 15) +
  ggplot2::theme_bw() +
  ggplot2::theme(
    legend.position = "none",
    plot.margin = margin(0, 0, 0, 10)
  )

importance_shap2_plot <- ggplot2::ggplot(importance_shap2, aes(
  x = reorder(Variable, Importance),
  y = Importance, fill = Variable
)) +
  ggplot2::geom_col(colour = "grey35", linewidth = 0.8, show.legend = FALSE) +
  ggplot2::coord_flip() +
  ggplot2::labs(
    x = "",
    y = "Importance"
  ) +
  ggplot2::ylim(0, 1.6) +
  ggplot2::theme_bw() +
  ggplot2::theme(
    legend.position = "none",
    plot.margin = margin(0, 0, 0, 10)
  )


all <- (importance_firm_plot / importance_firm2_plot + plot_layout(tag_level = "new")) |
  (importance_perm_plot / importance_perm2_plot + plot_layout(tag_level = "new")) |
  (importance_shap_plot / importance_shap2_plot + plot_layout(tag_level = "new")) +
    plot_layout(guides = "collect")
all + plot_annotation(
  title = "Variable importances",
  tag_levels = c("A", "1"),
  theme = ggplot2::theme_bw() + ggplot2::theme(
    plot.title = ggplot2::element_text(size = 20)
  )
)
```

Fig. \@ref(fig:importanceshap) and Fig. \@ref(fig:importanceshap2) show the effect of each feature on the FLOSS score. The more evident difference is the shape of the effect of `time_constraint` that initially suggested better results with larger values. However, removing the interactions seems to be a flat line.

Based on Fig. \@ref(fig:importance) and Fig. \@ref(fig:importanceshap2) we can infer that:

-   **`regime_threshold`**: is the most important feature, has an optimal value to be set, and since the high interaction with the `regime_landmark`, both must be tuned simultaneously. In this setting, high thresholds significantly impact the score, probably due to an increase in false positives starting on \>0.65 the overall impact is mostly negative.

-   **`regime_landmark`**: is not as important as the `regime_threshold,` but since there is a high interaction, it must not be underestimated. It is known that the *Arc Counts* have more uncertainty as we approach the margin of the streaming, and this becomes evident looking at how the score is negatively affected for values below 3.5s.

-   **`window_size`**: has a near zero impact on the score when correctly set. Nevertheless, for higher window values, the score is negatively affected. This high value probably depends on the data domain. In this setting, the model is being tuned towards the changes from atrial fibrillation/non-fibrillation; thus, the "shape of interest" is small compared to the whole heartbeat waveform. Window sizes smaller than 150 are more suitable in this case. As Beyer *et al.* noted, "as dimensionality increases, the distance to the nearest data point approaches the distance to the farthest data point" [@Beyer1999], which means that the bigger the window size, the smaller will be the contrast between different regimes.

-   **`mp_threshold`**: has a fair impact on the score, but primarily by *not using it*. We start to see a negative impact on the score with values above 0.60, while a constant positive impact with lower values.

-   **`time_constraint`**: is a parameter that must be interpreted cautiously. The 0 (zero) value means **no constraint**, which is equivalent to the size of the FLOSS history buffer (in our setting, 5000). We can see that this parameter's impact throughout the possible values is constantly near zero.

In short, for the MP computation, the parameter that is worth tuning is the `window_size`, while for the FLOSS computation, both `regime_threshold` (mainly) and `regime_landmark` shall be tuned.

```{r importanceshap, message=FALSE, fig.height = 6, fig.width= 10, out.width="100%"}
#| fig.cap="This shows the effect each variable has on the FLOSS score. This plot doesn't take into account the
#|  variable interactions."

trained_model <- readRDS(here("output", "dbarts_fitted_lmk.rds"))
train_data <- trained_model$training_data
testing_data <- trained_model$testing_data
predictors_names <- c("time_constraint", "regime_threshold", "mp_threshold", "window_size", "regime_landmark")
outcome_name <- "mean"

testing_data2 <- testing_data %>%
  dplyr::mutate(
    int_rt_rl = regime_threshold * regime_landmark,
    int_mp_w = mp_threshold * window_size,
    int_mp_rl = mp_threshold * regime_landmark,
    # int_mp_rt = mp_threshold * regime_threshold
    # int_mp_tc = mp_threshold * time_constraint
    .before = mean
  )

layout <- "
AABB
CCDD
#EE#
"
# all <- t1 + t2 + t3 + t4 + t5 +
#   plot_layout(design = layout, guides = "collect")

# all + plot_annotation(
#   title = "Shapley value vs. variable values",
#   subtitle = "Original fit",
#   theme = ggplot2::theme_bw()
# )

d1 <- shapviz::shapviz(shap_fastshap_all_test, X = testing_data[, predictors_names], baseline = mean(testing_data[, outcome_name]$mean))
a1 <- shapviz::sv_dependence(d1, "window_size", color_var = "auto") +
  ggplot2::geom_smooth(method = loess, colour = "#0000ff44", alpha = 0.2) +
  ggplot2::labs(y = ggplot2::element_blank()) +
  ggplot2::theme_bw()
a2 <- shapviz::sv_dependence(d1, "regime_threshold", color_var = "auto") +
  ggplot2::geom_smooth(method = loess, colour = "#0000ff44", alpha = 0.2) +
  ggplot2::labs(y = ggplot2::element_blank()) +
  ggplot2::theme_bw()
a3 <- shapviz::sv_dependence(d1, "regime_landmark", color_var = "auto") +
  ggplot2::geom_smooth(method = loess, colour = "#0000ff44", alpha = 0.2) +
  ggplot2::labs(y = ggplot2::element_blank()) +
  ggplot2::theme_bw()
a4 <- shapviz::sv_dependence(d1, "mp_threshold", color_var = "auto") +
  ggplot2::geom_smooth(method = loess, colour = "#0000ff44", alpha = 0.2) +
  ggplot2::labs(y = ggplot2::element_blank()) +
  ggplot2::theme_bw()
a5 <- shapviz::sv_dependence(d1, "time_constraint", color_var = "auto") +
  ggplot2::geom_smooth(method = loess, colour = "#0000ff44", alpha = 0.2) +
  ggplot2::labs(y = ggplot2::element_blank()) +
  ggplot2::theme_bw()

all2 <- ((a1 + a5 + a3 + a4 +
  plot_layout(guides = "collect")) +
  (a2 + plot_layout(guides = "collect")))
all2 + plot_layout(design = layout) + plot_annotation(
  title = "Shapley value vs. variable values",
  theme = ggplot2::theme_bw()
)
```

```{r importanceshap2, message=FALSE, fig.height = 6, fig.width= 10, out.width="100%"}
#| fig.cap="This shows the effect each variable has on the FLOSS score, taking into account the interactions."

# t1 <- autoplot(shap_fastshap_all_test2,
#   type = "dependence",
#   X = testing_data2, feature = predictors_names[2], alpha = 0.2
# ) + ggplot2::geom_smooth(method = loess) + ggplot2::labs(y = ggplot2::element_blank()) + ggplot2::theme_bw()
# t2 <- autoplot(shap_fastshap_all_test2,
#   type = "dependence",
#   X = testing_data2, feature = predictors_names[1], alpha = 0.2
# ) + ggplot2::geom_smooth(method = loess) + ggplot2::labs(y = ggplot2::element_blank()) + ggplot2::theme_bw()
# t3 <- autoplot(shap_fastshap_all_test2,
#   type = "dependence",
#   X = testing_data2, feature = predictors_names[4], alpha = 0.2
# ) + ggplot2::geom_smooth(method = loess) + ggplot2::labs(y = ggplot2::element_blank()) + ggplot2::theme_bw()
# t4 <- autoplot(shap_fastshap_all_test2,
#   type = "dependence",
#   X = testing_data2, feature = predictors_names[3], alpha = 0.2
# ) + ggplot2::geom_smooth(method = loess) + ggplot2::labs(y = ggplot2::element_blank()) + ggplot2::theme_bw()
# t5 <- autoplot(shap_fastshap_all_test,
#   type = "dependence",
#   X = testing_data2, feature = predictors_names[5], alpha = 0.2
# ) + ggplot2::geom_smooth(method = loess) + ggplot2::labs(y = ggplot2::element_blank()) + ggplot2::theme_bw()


layout <- "
AABB
CCDD
#EE#
"


d1 <- shapviz::shapviz(shap_fastshap_all_test2, X = testing_data2[, predictors_names], baseline = mean(testing_data2[, outcome_name]$mean))
a1 <- shapviz::sv_dependence(d1, "window_size", color_var = "auto") +
  ggplot2::geom_smooth(method = loess, colour = "#0000ff44", alpha = 0.2) +
  ggplot2::labs(y = ggplot2::element_blank()) +
  ggplot2::theme_bw()
a2 <- shapviz::sv_dependence(d1, "regime_threshold", color_var = "auto") +
  ggplot2::geom_smooth(method = loess, colour = "#0000ff44", alpha = 0.2) +
  ggplot2::labs(y = ggplot2::element_blank()) +
  ggplot2::theme_bw()
a3 <- shapviz::sv_dependence(d1, "regime_landmark", color_var = "auto") +
  ggplot2::geom_smooth(method = loess, colour = "#0000ff44", alpha = 0.2) +
  ggplot2::labs(y = ggplot2::element_blank()) +
  ggplot2::theme_bw()
a4 <- shapviz::sv_dependence(d1, "mp_threshold", color_var = "auto") +
  ggplot2::geom_smooth(method = loess, colour = "#0000ff44", alpha = 0.2) +
  ggplot2::labs(y = ggplot2::element_blank()) +
  ggplot2::theme_bw()
a5 <- shapviz::sv_dependence(d1, "time_constraint", color_var = "auto") +
  ggplot2::geom_smooth(method = loess, colour = "#0000ff44", alpha = 0.2) +
  ggplot2::labs(y = ggplot2::element_blank()) +
  ggplot2::theme_bw()

all2 <- ((a1 + a5 + a3 + a4 +
  plot_layout(guides = "collect")) +
  (a2 + plot_layout(guides = "collect")))
all2 + plot_layout(design = layout) + plot_annotation(
  title = "Shapley value vs. variable values",
  theme = ggplot2::theme_bw()
)
```

According to the FLOSS paper [@gharghabi2018], the `window_size` is indeed a feature that can be tuned; nevertheless, the results appear to be similar in a reasonably wide range of window sizes, up to a limit, consistent with our findings.

### Visualizing the predictions

At this point, the grid search tested a total of 23,389 models with resulting (individual) scores from 0.0002 to 1669.83 (Q25: 0.9838, Q50: 1.8093, Q75: 3.3890).

#### By recording

First, we will visualize how the models (in general) performed throughout the individual recordings.

Fig. \@ref(fig:global) shows a violin plot of equal areas clipped to the minimum value. The blue color indicates the recordings with a small IQR (interquartile range) of model scores. We see on the left half 10% of the recordings with the worst minimum score, and on the right half, 10% of the recordings with the best minimum score.

Next, we will visualize some of these predictions to understand why some recordings were difficult to segment. For us to have a simple baseline: a recording with just one regime change, and the model predicts exactly one regime change, but far from the truth, the score will be roughly 1.

```{r global, eval=TRUE, fig.height=5, fig.width=10, out.width="100%"}
#| fig.cap="Violin plot showing the distribution of the FLOSS score achieved by all tested models by
#|  recording.  The left half shows the recordings that were difficult to predict (10% overall), whereas
#|  the right half shows the recordings that at least one model could achieve a good prediction (10%
#|  overall).  The recordings are sorted (left-right) by the minimum (best) score achieved in descending
#|  order, and ties are sorted by the median of all recording scores.  The blue color highlights
#|  recordings where models had an IQR variability of less than one.  As a simple example, a recording
#|  with just one regime change, and the model predicts exactly one change, far from the truth, the
#|  score will be roughly 1."

all_scores <- readRDS(here::here("output/regime_outputs_lmk.rds"))

scores_stats <- all_scores %>%
  dplyr::select(record, score) %>%
  dplyr::group_by(record) %>%
  dplyr::reframe(
    score = score, min = min(score), q25 = quantile(score, 0.25),
    median = quantile(score, 0.5), q75 = quantile(score, 0.75),
    mean = mean(score), max = max(score)
  )

records_factors <- forcats::as_factor(scores_stats$record)

scores_stats$id <- sprintf("%03d", (as.numeric(records_factors)))

scores_stats %>%
  dplyr::mutate(low_iqr = q75 - q25 < 1) %>%
  dplyr::filter(min > quantile(min, 0.9) | min < quantile(min, 0.1)) %>%
  ggplot2::ggplot(aes(x = reorder(reorder(id, -median), -min), y = score, colour = low_iqr)) +
  ggplot2::geom_violin() +
  ggplot2::coord_cartesian(ylim = c(0, 5)) +
  ggplot2::theme_bw() +
  ggplot2::labs(title = "Scores by recording", colour = "IQR < 1", x = "Recording ID", y = "Score distribution")
```

Fig. \@ref(fig:worst) shows the best effort in predicting the most complex recordings. One information not declared before is that if the model does not predict any change, it will put a mark on the zero position. On the other side, the truth markers positioned at the beginning and the end of the recording were removed, as these locations lack information and do not represent a streaming setting.

```{r worst, eval=TRUE, fig.height=15, fig.width=12, out.width="100%", dev="png"}
#| fig.cap="Prediction of the worst 10% of recordings (red is the truth, blue are the predictions)."

worst <- scores_stats %>%
  dplyr::filter(min > quantile(min, 0.9)) %>%
  dplyr::group_by(record) %>%
  dplyr::slice_head() %>%
  dplyr::ungroup() %>%
  dplyr::arrange(desc(min), desc(median)) %>%
  dplyr::slice_head(n = 10) %>%
  dplyr::pull(record)
# "data_25_1.par"  "data_32_12.par" "data_85_1.par"  "data_90_1.par"  "data_68_2.par"

worst_data <- all_scores %>%
  dplyr::filter(record %in% worst) %>%
  dplyr::group_by(record) %>%
  dplyr::slice_min(n = 1, order_by = score, with_ties = FALSE) %>%
  dplyr::arrange(desc(score)) %>%
  dplyr::ungroup()

plots <- list()
for (i in seq_len(nrow(worst_data))) {
  plots[[i]] <- tkplot(worst_data[i, ], FALSE, 50)
}

wrap_plots(plots, ncol = 1)
```

Fig. \@ref(fig:best) shows the best performances of the best recordings. Notice that there are recordings with a significant duration and few regime changes, making it hard for a "trivial model" to predict randomly.

```{r best, eval=TRUE, fig.height=15, fig.width=12, out.width="100%", dev="png"}
#| fig.cap="Prediction of the best 10% of recordings (red is the truth, blue are the predictions)."

bests <- scores_stats %>%
  dplyr::filter(min < quantile(min, 0.1)) %>%
  dplyr::group_by(record) %>%
  dplyr::slice_head() %>%
  dplyr::ungroup() %>%
  dplyr::arrange(desc(min), desc(median)) %>%
  dplyr::slice_head(n = 10) %>%
  dplyr::pull(record)

bests_data <- all_scores %>%
  dplyr::filter(record %in% bests) %>%
  dplyr::group_by(record) %>%
  dplyr::slice_min(n = 1, order_by = score, with_ties = FALSE) %>%
  dplyr::arrange(desc(score)) %>%
  dplyr::ungroup()

plots <- list()
for (i in seq_len(nrow(bests_data))) {
  plots[[i]] <- tkplot(bests_data[i, ], FALSE, 50)
}

wrap_plots(plots, ncol = 1)
```

#### By model

Fig. \@ref(fig:globalmodel) shows the distribution of the FLOSS score of the 10% worst (left side) and 10% best models across the recordings (right side). The bluish color highlights the models with SD below 3 and IQR below 1.

```{r globalmodel, eval=TRUE, fig.height=5, fig.width=10, out.width="100%"}
#| fig.cap="Violin plot showing the distribution of the FLOSS score achieved by all tested models during the
#|  inner ressample.  The left half shows the models with the worst performances (10% overall), whereas
#|  the right half shows the models with the best performances (10% overall).
#|  The models are sorted (left-right) by the mean score (top) and by the median (below). Ties are
#|  sorted by the SD and IQR, respectively.  The bluish colors highlights models with an SD below 3
#|  and IQR below 1."

# cores_stats_model <- all_scores %>% dplyr::mutate(across(all_of(predictors_names), as.factor, .unpack = FALSE))

if (file.exists(here("output", "scores_stats_model_rep.rds"))) {
  scores_stats_model <- readRDS(here("output", "scores_stats_model_rep.rds"))
} else {
  cli::cli_abort("File not found.")
  scores_stats_model <- all_scores %>%
    dplyr::group_by(across(all_of(predictors_names))) %>%
    dplyr::mutate(model = glue("{window_size}_{time_constraint}_{mp_threshold}_{regime_threshold}_{regime_landmark}")) %>%
    dplyr::reframe(
      record = record,
      model = model,
      score = score, min = min(score), q25 = quantile(score, 0.25),
      median = quantile(score, 0.5), q75 = quantile(score, 0.75),
      iqr = q75 - q25,
      mean = mean(score), max = max(score),
      sd = sd(score)
    )
  saveRDS(scores_stats_model, file = here("output", "scores_stats_model_rep.rds"))
}


scores_stats_model$id <- (sprintf("%05d", (as.numeric(as.factor(scores_stats_model$model)))))
scores_stats_model$id_text <- (sprintf("Model_%05d", (as.numeric(as.factor(scores_stats_model$model)))))
scores_stats_model$record <- (sprintf("%03d", (as.numeric(factor(scores_stats_model$record, labels = levels(records_factors))))))
scores_stats_model <- scores_stats_model %>% dplyr::select(-model)

low <- head(sort(unique(scores_stats_model$mean)), 20)
high <- tail(sort(unique(scores_stats_model$mean)), 20)

model_mean <- scores_stats_model %>%
  dplyr::mutate(low_sd = sd < 3) %>%
  dplyr::filter(mean > high | mean < low) %>%
  ggplot2::ggplot(aes(x = reorder(reorder(id, -sd), -mean), y = score, colour = low_sd)) +
  ggplot2::scale_colour_manual(values = c("FALSE" = "#ff0000c2", "TRUE" = "#0000ffb5")) +
  ggplot2::geom_violin() +
  ggplot2::coord_cartesian(ylim = c(0, 3)) +
  ggplot2::theme_bw() +
  ggplot2::theme(axis.text.x = element_text(size = 8, angle = 90, vjust = 0.5, hjust = 1)) +
  ggplot2::labs(subtitle = "Ordered by Mean and SD", colour = "SD < 3", x = ggplot2::element_blank(), y = "Score distribution")

low <- head(sort(unique(scores_stats_model$median)), 20)
high <- tail(sort(unique(scores_stats_model$median)), 20)

model_median <- scores_stats_model %>%
  dplyr::mutate(low_iqr = q75 - q25 < 1) %>%
  dplyr::filter(median > high | median < low) %>%
  ggplot2::ggplot(aes(x = reorder(reorder(id, -iqr), -median), y = score, colour = low_iqr)) +
  ggplot2::geom_violin() +
  ggplot2::coord_cartesian(ylim = c(0, 3)) +
  ggplot2::theme_bw() +
  ggplot2::theme(axis.text.x = element_text(size = 8, angle = 90, vjust = 0.5, hjust = 1)) +
  ggplot2::labs(subtitle = "Ordered by Median and IQR", colour = "IQR < 1", x = "Model ID", y = "Score distribution")

(model_mean / model_median) + plot_layout(guides = "auto") +
  plot_annotation(
    title = "Scores grouped by model",
    theme = ggplot2::theme_bw()
  )
```

Fig. \@ref(fig:bestmodels) the performance of the six best models. They are ordered from left to right, from the worst record to the best record. The top model is the one with the lowest mean across the scores. The blue line indicates the mean score, and the red line the median score. The scores above 3 are squished in the plot and colored according to the scale in the legend.

```{r bestmodels, eval=TRUE, fig.height=22, fig.width=18, out.width="100%"}
#| fig.cap="Performances of the best 6 models across all inner resample of recordings.
#|  The recordings are ordered by score, from the worst to the best.
#|  Each plot shows one model, starting from the best one.
#|  The red line indicates the median score of the model. The blue line indicates the mean score of the model.
#|  The gray line limits the zero-score region. The plot is limited on the \"y\" axis, and the scores above this
#|  limit are shown in color."

best_models <- scores_stats_model %>%
  dplyr::filter(mean < quantile(mean, 0.1)) %>%
  dplyr::arrange(mean, sd) %>%
  dplyr::pull(id_text) %>%
  unique() %>%
  .[1:6]

plots <- list()
for (i in seq_len(length(best_models))) {
  dd <- scores_stats_model %>% dplyr::filter(id_text == best_models[i])
  plots[[i]] <- ggplot2::ggplot(dd, aes(x = reorder(record, -score), y = score, colour = score)) +
    ggplot2::geom_point(size = 2) +
    ggplot2::geom_hline(aes(yintercept = median), colour = "red") +
    ggplot2::geom_hline(aes(yintercept = mean), colour = "blue") +
    ggplot2::geom_hline(aes(yintercept = 0), colour = "gray50") +
    ggplot2::scale_y_continuous(
      limits = c(0, 3),
      oob = scales::oob_squish,
      expand = c(0.1, 0.05, 0.1, -0.1)
    ) +
    ggplot2::scale_color_gradientn(colors = c("#ffd500db", "#ff8800", "#ff5e00", "#ff0000"), limits = c(3.1, 100)) +
    ggplot2::theme_bw(base_size = 15) +
    ggplot2::theme(axis.text.x = element_text(size = 9, angle = 90, vjust = 0.5, hjust = 1)) +
    ggplot2::labs(
      title = best_models[i],
      colour = "Score out\nof bounds",
      x = ifelse(i == length(best_models), "Record ID", ""),
      y = ggplot2::element_blank()
    )
}

wrap_plots(plots, ncol = 1, guides = "collect") + plot_annotation(
  title = "Performances of the 6 best models",
  theme = ggplot2::theme_bw()
)
```

#### The Holdout

Finally, Table \@ref(tab:tblbestparam) shows a summary of the best five models across all the inner resample (cross-validation). The column `mean`
shows the average score, and column `std_err` shows the standard error of the mean. The column `holdout` shows the
final score of this model on the holdout set (outer resample).


```{r tblbestparam, echo=FALSE}
if (file.exists(here("output", "regime_outputs_holdout.rds"))) {
  regime_outputs_holdout <- readRDS(here("output", "regime_outputs_holdout.rds"))

  best_parameters <- regime_outputs_holdout$combined %>%
    dplyr::group_by(
      window_size, regime_threshold, regime_landmark, .metric, .estimator
    ) %>%
    dplyr::summarise(mean = mean(.estimate), std_err = sd(.estimate), n = dplyr::n()) %>%
    dplyr::ungroup() %>%
    dplyr::arrange(mean, std_err) %>%
    dplyr::slice_min(n = 5, order_by = mean) %>%
    dplyr::select(-.metric, -.estimator, -n)

  holdout <- regime_outputs_holdout$results %>%
    dplyr::select(.estimate) %>%
    dplyr::rename(holdout = .estimate)

  best_parameters <- dplyr::bind_cols(best_parameters, holdout)

  kbl(best_parameters,
    booktabs = TRUE,
    longtable = TRUE,
    caption = "Summary of the five best models. The `mean` shows the inner resample average score. The `holdout` shows the final score of the model on the holdout set (outer resample).",
    align = "cccccc",
    digits = 2,
    position = "ht",
    linesep = ""
  ) %>%
    row_spec(0, bold = TRUE) %>%
    column_spec(6, bold = TRUE) |>
    kable_styling(full_width = TRUE)
}
```


## Classification

As described in section \@ref(classregime), the classification algorithm is based on the Contrast Profile (CP) [@Mercer2021]. The CP allows us
to detect sequences that at the same time very *similar* to its neighbors in class *A* while is very *different* from the nearest neighbor from class *B*.

The variables parameters we can tune are the following:

-   **`shapelet_size`**: that we will use interchangeably with the term *window_size* as is defines the size of the rolling window used to compute the CP. It was used a series of shapelet sizes with exponential distribution, resulting on `20` values from `21` to `401`.
-   **`top_k`**: how many shapelets we will select on each CP, being the first one the shapelet with the highest constrast value. It was chosen the default value of `10`.
-   **`max_shapelets`**: the maximum number of shapelets we will allow on selecting the shapelet set. It was set as `20` to allow more freedom on the selection of shapelets.
-   **`max_redundance`**: the maximum number of "redundance" we will allow on selecting the shapelet set. The redundance means more than one shapelet can classify correctly the same observation. It was set as `10`, to also allow more freedom on the selection of shapelets.

Fig. \@ref(fig:classthepipeline) shows the workflow using Nested resamplig as described on \@ref(methodology).

```{r class_cached, echo=FALSE, cache=FALSE}
network <- readRDS(here::here("output", "classification_network.rds"))
net <- visNetwork::visPhysics(network, hierarchicalRepulsion = list(
  springLength = 1,
  avoidOverlap = 0.5,
  nodeDistance = 120
))
```

```{r classthepipeline, out.width="100%", fig.cap="Classification pipeline."}
if (knitr::is_latex_output()) {
  my_graphics("theclass-network")
} else {
  visNetwork::visInteraction(net, hover = TRUE, multiselect = TRUE, tooltipDelay = 100)
}
```

The dataset used for working Classification algorithm was the CinC/Physionet Challenge 2015 was about "Reducing False Arrhythmia Alarms in the ICU [@Clifford2015].

The selected records were those that contain ventricular tachycardia. The last 10 seconds (at 250hz) of all records were selected and grouped as `TRUE` alarm and `FALSE` alarm. A total of 331 records were used, being 245 `FALSE` alarms, and 86 `TRUE` alarms.

The records were split in a proportion of 3/4 (248) for the training set (inner resampling) and 1/4 (83) for the test set (outer resampling). The proportions of `TRUE` and `FALSE` alarms were similar to the original dataset: 184 `FALSE` alarms and 64 `TRUE` alarms in the training set, and 61 `FALSE` alarms and 22 `TRUE` alarms in the test set. The inner resampling was performed using a 5-fold cross-validation.

In order to compute the Contrast Profile (CP), on each fold, the `TRUE` alarms were concatenated in one single time-series with a small gap of 300 observation of random noise in order to isolate each alarm. The same was done for the `FALSE` alarms.

The following steps were performed for each fold:

1.  The CP was computed with several shapelet sizes (from 21 to 401), and top 10 best shapelet candidates were stored based on the CP values.
2.  Every shapelet candidate was evaluated for the hability to classify as `TRUE` or `FALSE` alarm along all the concatenated time-series.
    a.  First, the distance profile of the shapelet was computed against the `FALSE` time-series, and threshold was set in order to not detect any `FALSE` alarm as `TRUE` alarm.
    b.  Second, the distance profile of the shapelet was computed against the `TRUE` time-series, and using the threshold computed in the previous step, the number of `TRUE` alarms detected was recorded and called the "coverage" of the shapelet.
3.  With this information, a Beam Search was performed in order to select the best set of shapelets that maximize the coverage and minimize the number of shapelets. The set of shapelets may contain shapelets of different sizes.
4.  The confusion matrix of these sets of shapelets was computed. Other metrics were also computed, such as the F1-score, the accuracy, the precision, the specificity, the MCC, the $\kappa_m$ and the Kappa.

An example of candidates for ventricular tachycardia is presented on \@ref(fig:holdout).

```{r holdout, echo=FALSE, out.width="70%", fig.height=8, fig.width=8}
#| fig-cap: "Set of shapelets for the classification task. Above shows the model number, the coverage (the proportion of true alarms that were detected) and redundancy during the analysis_split (inner resample),
#|  The Precision, the Specificity, and the Km during the testing_split (outer resample)."

data <- readRDS(here("output", "classification_holdout.rds"))

data$plots[[4]] # + ggplot2::facet_wrap(~name, ncol = 1, scales = "free_y")
```


After the Inner Resampling is done, the best sets of shapelets are selected and evaluated on the Test Set without retraining a new Contrast Profile. Thus assessing the generalization of the shapelet set on new data.

The criteria to select the best sets of shapelets was described on section \@ref(classification-criteria) being the Precision the ranking criteria. It was also required that the set being present on more than one fold and in both repetitions. Also, the sets of shapelets that had a negative $\kappa_m$ were discarded.

The following results were obtained:

```{r tabmodel, echo=FALSE, out.width="70%"}
colnames <- c("tp", "fp", "tn", "fn", "precision", "recall", "specificity", "accuracy", "f1", "mcc", "km", "kappa")

kbl(data$model %>% dplyr::select(dplyr::all_of(colnames)),
  booktabs = TRUE,
  longtable = TRUE,
  caption = "Performance of the best five sets of shapelet on the inner resample.",
  # align = "cccccccccc",
  digits = 2,
  position = "ht",
  linesep = ""
) %>%
  row_spec(0, bold = TRUE) |>
  kable_styling(full_width = TRUE)
```

```{r tabmetric, echo=FALSE, out.width="70%"}
kbl(data$metric %>% dplyr::select(dplyr::all_of(colnames)),
  booktabs = TRUE,
  longtable = TRUE,
  caption = "Performance of the best five sets of shapelet on outer resample.",
  # align = "cccccccccc",
  digits = 2,
  position = "ht",
  linesep = ""
) %>%
  row_spec(0, bold = TRUE) |>
  kable_styling(full_width = TRUE)
```

```{r taboverall, echo=FALSE, out.width="70%"}
colnames <- c("precision", "recall", "specificity", "accuracy", "f1_micro", "f1_macro", "mcc", "km", "kappa")

kbl(data$overall %>% dplyr::select(dplyr::all_of(colnames)),
  booktabs = TRUE,
  longtable = TRUE,
  caption = "The aggregated performance of the best five sets of shapelet on the outer resample. Median, Q25 and Q75",
  # align = "cccccccccc",
  digits = 2,
  position = "ht",
  linesep = ""
) %>%
  row_spec(0, bold = TRUE) |>
  kable_styling(full_width = TRUE)
```

## Feasibility trial

A side-project called "false.alarm.io" has been derived from this work (an unfortunate mix of
"false.alarm" and "PlatformIO" [@PlatformIO], the IDE chosen to interface the panoply of embedded
systems we can experiment with). The current results of this side-project are very enlightening and
show that the final algorithm can indeed be used in small hardware. Further data will be available
in the future.

A brief mentioning, linking back to the objectives of this work, an initial trial was done using an
ESP32 MCU (Fig. \@ref(fig:esp32)) in order to be sure if such a small device can handle the task.

```{r esp32, echo=FALSE, out.width="50%", fig.cap="ESP32 MCU."}
my_graphics("esp32", "figure")
```

Current results show that such device has enough computation power to handle the task in real-time
using just one of its two microprocessors. The main limitation seen in advance is the on-chip SRAM
that must be well managed.

# Scientific contributions

## Matrix Profile

Since the first paper that presented this new concept [@Yeh2017a], many investigations have been
made to speed its computation. It is notable how all computations are not dependent on the _rolling
window size_ as previous works not using Matrix Profile. Aside from this, we can see that the first
STAMP [@Yeh2017a] algorithm has the time complexity of $O(n^2log{n})$ while STOMP [@zhu2016]
$O(n^2)$ (a significant improvement), but STOMP lacks the "any-time" property. Later SCRIMP
[@zhu2018] solves this problem keeping the same time complexity of $O(n^2)$. Here we are in the
"exact" algorithms domain, and we will not extend the scope for conciseness.

The main issue with the algorithms above is the dependency on a fast Fourier transform (FFT)
library. FFT has been extensively optimized and architecture/CPU bounded to exploit the most of
speed. Also, padding data to some power of 2 happens to increase the algorithm's efficiency. We can
argue that time complexity doesn't mean "faster" when we can exploit low-level instructions. In our
case, using FFT in a low-power device is overkilling. For example, a quick search over the internet
gives us a hint that computing FFT on 4096 data in an ESP32 takes about 21ms (~47 computations in 1
second). This means ~79 seconds for computing all FFT's (~3797) required for STAMP using a window of 300.
Currently, we can compute a full matrix of 5k data in about 9 seconds in an ESP32 MCU (Fig.
\@ref(fig:esp32)), and keep updating it as fast as 1 min of data (at 250hz) in just 6 seconds.

Recent works using _exact_ algorithms are using an unpublished algorithm called **MPX**, which
computes the Matrix Profile using cross-correlation methods ending up faster and is easily
portable.

**On computing the Matrix Profile:** the contribution of this work on this area is adding the
*Online* capability to MPX, which means we can update the Matrix Profile as new data comes in.

**On extending the Matrix Profile:** the contribution of this work on this area is the use of an
unexplored constraint that we could apply on building the Matrix Profile we are calling _Similarity
Threshold_ (ST). The original work outputs the similarity values in Euclidean Distance (ED) values,
while MPX naturally outputs the Pearson's correlation coefficients (CC) values. Both ED and CC
are interchangeable using the Equation $\eqref{edcc}$. However, we may argue that it is easier to
compare values that do not depend on the window size during an exploratory phase. MPX happens to
naturally return values in CC, saving a few more computation time. The ST is an interesting factor
that we can use, especially when detecting pattern changes during time. The FLOSS algorithm relies
on counting references between indexes in the time series. ST can help remove "noise" from these
references since only similar patterns above a certain threshold are referenced, and changes have
more impact on these counts. The best ST threshold is still to be determined.

\
$$
CC = 1 - \frac{ED}{(2 \times WindowSize)} \tag{9} \label{edcc}
$$
\

## Regime change detection

In the original paper, in chapter 3.5, the authors of FLOSS wisely introduce the **temporal
constraint**, which improves the sensitivity of regime change detection on situations where
a regime may alternate in short periods.

Nevertheless, the authors declare the correction curve typically used on the algorithm as "simply a
uniform distribution", but this is not an accurate statement. The _Arc Counts_ of newly incoming
data is truncated by the same amount of temporal constraint. This prevents detecting a regime change
in the last 10 seconds completely, as this thesis requires.

The main contribution of this work in this area is overcoming this issue by computing the
theoretical distribution using the temporal constraint parameters beforehand. as shown in Fig.
\@ref(fig:distributions). That gives us enough data to evaluate a regime change accurately utilizing
a minimum of $2 \times WindowSize$ datapoints.

```{r dist-data, message=FALSE, warning=FALSE, include=FALSE}
source(here("scripts", "common", "compute_floss.R"))

get_dist <- function(mp_const = 1250, floss_const = 0) {
  set.seed(2021)
  iac <- list()
  pro_size <- 5000
  mp_time_constraint <- mp_const
  floss_time_constraint <- floss_const
  for (i in 1:100) {
    iac[[i]] <- get_asym(pro_size, mp_time_constraint, floss_time_constraint)
  }

  aic_avg <- rowMeans(as.data.frame(iac))

  data.frame(index = 1:5000, counts = aic_avg)
}

data_5000 <- get_dist(5000)
data_4250 <- get_dist(4250)
data_2500 <- get_dist(2500)
data_1250 <- get_dist(1250)

floss_data_5000 <- get_dist(0, 5000)
floss_data_4250 <- get_dist(0, 4250)
floss_data_2500 <- get_dist(0, 2500)
floss_data_1250 <- get_dist(0, 1250)
```

```{r distributions, echo=FALSE, fig.cap="1D-IAC distributions for earlier temporal constraint (on Matrix Profile).", message=FALSE, warning=FALSE}
floss_dist <- ggplot(data_5000, aes(index, counts)) +
  geom_line(size = 0.1) +
  ggtitle("a) No constraint") +
  theme_grey(base_size = 7)

floss_4250 <- ggplot(data_4250, aes(index, counts)) +
  geom_line(size = 0.1) +
  ggtitle("b) Constraint of 4250") +
  theme_grey(base_size = 7)

floss_2500 <- ggplot(data_2500, aes(index, counts)) +
  geom_line(size = 0.1) +
  annotate("segment", y = 0, yend = max(data_2500$counts), x = 2500, xend = 2500, linetype = 2, size = 0.1) +
  annotate("text", x = 2500 - 80, y = 40, label = "start", color = "black", size = 2, angle = 90, hjust = 0) +
  annotate("segment", y = 0, yend = max(data_2500$counts), x = 5000 - 2500 * 0.9, xend = 5000 - 2500 * 0.9, linetype = 2, size = 0.1) +
  annotate("text", x = 5000 - 2500 * 0.9 - 80, y = 40, label = "end", color = "black", size = 2, angle = 90, hjust = 0) +
  ggtitle("c) Constraint of 2500") +
  theme_grey(base_size = 7)

floss_1250 <- ggplot(data_1250, aes(index, counts)) +
  geom_line(size = 0.1) +
  annotate("segment", y = 0, yend = max(data_1250$counts), x = 1250, xend = 1250, linetype = 2, size = 0.1) +
  annotate("text", x = 1250 - 80, y = 40, label = "start", color = "black", size = 2, angle = 90, hjust = 0) +
  annotate("segment", y = 0, yend = max(data_1250$counts), x = 5000 - 1250 * 0.9, xend = 5000 - 1250 * 0.9, linetype = 2, size = 0.1) +
  annotate("text", x = 5000 - 1250 * 0.9 - 80, y = 40, label = "end", color = "black", size = 2, angle = 90, hjust = 0) +
  ggtitle("d) Constraint of 1250") +
  theme_grey(base_size = 7)

gg <- gridExtra::arrangeGrob(floss_dist, floss_4250, floss_2500, floss_1250,
  nrow = 2,
  bottom = grid::textGrob(paste("The plot a) shows the distribution used for the arc count correction when there is no time constraint.", "\n", "b) Shows a constraint of 3/4 of the total. c) 1/2 of the total. d) 1/4 of the total; here we see clearly the flat line.", "\n", "The dashed line marks the start and the end of the uniform zone."), just = "center", gp = grid::gpar(fontsize = 7))
)

grid::grid.draw(gg)
```

# Scientific outcomes

This research has already yielded two R packages concerning the MP algorithms from UCR [@mpucr]. The
first package is called `tsmp`, and a paper was also published in the R Journal [@RJ-2020-021]
(Journal Impact Factor™, 2020 of 3.984). The second package is called `matrixprofiler` and enhances
the first one, using low-level language to improve computational speed. The author has also joined
the Matrix Profile Foundation is a co-founder with contributors from Python and Go languages
[@mpf2020; @VanBenschoten2020]. The benchmarks of the R implementation are available online
[@Bischoff2021a].

Additionally to the above publication and the publication of the ongoing literature survey, two
articles about this thesis subject will be published. The first regards the application of the FLOSS
algorithm on real-time ECG showing its potential on using on low-power devices. The second is
regarding the use of combined shapelets for relevant ECG patterns identification.

# Expected results and outcomes

In the end, this thesis will provide a framework for identifying life-threatening conditions using
biological streaming data on devices with low CPU and low memory specifications. We expect to
achieve a high-quality model on identifying these pathological conditions maintaining their
robustness in the presence of noise and artifacts seen on real-world applications.

# References