From 461d54e17c1243b77164ddb65768b1dd99b903f1 Mon Sep 17 00:00:00 2001
From: Anthony Li <38049831+doccstat@users.noreply.github.com>
Date: Thu, 23 May 2024 23:47:51 -0500
Subject: [PATCH] Clean up vignettes
---
DESCRIPTION | 11 -
vignettes/comparison-packages-ar_data-mcp.svg | 1275 ----
.../comparison-packages-mean_data_1-bcp.svg | 1365 -----
.../comparison-packages-mean_data_1-mcp.svg | 1264 ----
.../comparison-packages-mean_data_3-bcp.svg | 3373 -----------
.../comparison-packages-mv_data_1-mcp.svg | 2255 --------
.../comparison-packages-well_log-cpts.svg | 5102 -----------------
.../comparison-packages-well_log-time.svg | 497 --
vignettes/comparison-packages.Rmd | 2295 ++++++--
vignettes/comparison-packages.Rmd.original | 1358 +++++
vignettes/figure/ar3-mcp-result-1.png | Bin 0 -> 71750 bytes
vignettes/figure/data-setup-ar3-1.png | Bin 0 -> 34641 bytes
vignettes/figure/data-setup-garch11-1.png | Bin 0 -> 26736 bytes
vignettes/figure/data-setup-lasso-1.png | Bin 0 -> 49080 bytes
.../figure/data-setup-linear-regression-1.png | Bin 0 -> 41719 bytes
.../data-setup-logistic-regression-1.png | Bin 0 -> 52845 bytes
...ivariate-mean-and-or-variance-change-1.png | Bin 0 -> 29927 bytes
.../data-setup-multivariate-mean-change-1.png | Bin 0 -> 36866 bytes
.../data-setup-poisson-regression-1.png | Bin 0 -> 30561 bytes
.../data-setup-poisson-regression-2.png | Bin 0 -> 55363 bytes
...ivariate-mean-and-or-variance-change-1.png | Bin 0 -> 21988 bytes
.../data-setup-univariate-mean-change-1.png | Bin 0 -> 29307 bytes
vignettes/figure/data-setup-var2-1.png | Bin 0 -> 40305 bytes
.../detection-comparison-well-log-plot-1.png | Bin 0 -> 56791 bytes
.../multivariate-mean-change-bcp-result-1.png | Bin 0 -> 89528 bytes
.../time-comparison-well-log-plot-1.png | Bin 0 -> 7313 bytes
...an-and-or-variance-change-mcp-result-1.png | Bin 0 -> 53996 bytes
...an-and-or-variance-change-not-result-1.png | Bin 0 -> 32246 bytes
...an-and-or-variance-change-not-result-2.png | Bin 0 -> 33680 bytes
.../univariate-mean-change-bcp-result-1.png | Bin 0 -> 48235 bytes
.../univariate-mean-change-mcp-result-1.png | Bin 0 -> 60330 bytes
.../univariate-mean-change-not-result-1.png | Bin 0 -> 39442 bytes
32 files changed, 3104 insertions(+), 15691 deletions(-)
delete mode 100644 vignettes/comparison-packages-ar_data-mcp.svg
delete mode 100644 vignettes/comparison-packages-mean_data_1-bcp.svg
delete mode 100644 vignettes/comparison-packages-mean_data_1-mcp.svg
delete mode 100644 vignettes/comparison-packages-mean_data_3-bcp.svg
delete mode 100644 vignettes/comparison-packages-mv_data_1-mcp.svg
delete mode 100644 vignettes/comparison-packages-well_log-cpts.svg
delete mode 100644 vignettes/comparison-packages-well_log-time.svg
create mode 100644 vignettes/comparison-packages.Rmd.original
create mode 100644 vignettes/figure/ar3-mcp-result-1.png
create mode 100644 vignettes/figure/data-setup-ar3-1.png
create mode 100644 vignettes/figure/data-setup-garch11-1.png
create mode 100644 vignettes/figure/data-setup-lasso-1.png
create mode 100644 vignettes/figure/data-setup-linear-regression-1.png
create mode 100644 vignettes/figure/data-setup-logistic-regression-1.png
create mode 100644 vignettes/figure/data-setup-multivariate-mean-and-or-variance-change-1.png
create mode 100644 vignettes/figure/data-setup-multivariate-mean-change-1.png
create mode 100644 vignettes/figure/data-setup-poisson-regression-1.png
create mode 100644 vignettes/figure/data-setup-poisson-regression-2.png
create mode 100644 vignettes/figure/data-setup-univariate-mean-and-or-variance-change-1.png
create mode 100644 vignettes/figure/data-setup-univariate-mean-change-1.png
create mode 100644 vignettes/figure/data-setup-var2-1.png
create mode 100644 vignettes/figure/detection-comparison-well-log-plot-1.png
create mode 100644 vignettes/figure/multivariate-mean-change-bcp-result-1.png
create mode 100644 vignettes/figure/time-comparison-well-log-plot-1.png
create mode 100644 vignettes/figure/univariate-mean-and-or-variance-change-mcp-result-1.png
create mode 100644 vignettes/figure/univariate-mean-and-or-variance-change-not-result-1.png
create mode 100644 vignettes/figure/univariate-mean-and-or-variance-change-not-result-2.png
create mode 100644 vignettes/figure/univariate-mean-change-bcp-result-1.png
create mode 100644 vignettes/figure/univariate-mean-change-mcp-result-1.png
create mode 100644 vignettes/figure/univariate-mean-change-not-result-1.png
diff --git a/DESCRIPTION b/DESCRIPTION
index 5685c7bd..617417d1 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -36,30 +36,19 @@ Imports:
tseries
Suggests:
abind,
- breakfast,
- changepoint,
- cpm,
- CptNonPar,
dplyr,
- fpop,
ggplot2,
gridExtra,
- jointseg,
knitr,
lubridate,
matrixStats,
mockthat,
mvtnorm,
- not,
numDeriv,
RcppClock,
reshape2,
rmarkdown,
- segmented,
- stepR,
testthat (>= 3.0.0),
- VARDetect,
- wbs,
xml2,
zoo
LinkingTo:
diff --git a/vignettes/comparison-packages-ar_data-mcp.svg b/vignettes/comparison-packages-ar_data-mcp.svg
deleted file mode 100644
index 8b4cbbb1..00000000
--- a/vignettes/comparison-packages-ar_data-mcp.svg
+++ /dev/null
@@ -1,1275 +0,0 @@
-
-
diff --git a/vignettes/comparison-packages-mean_data_1-bcp.svg b/vignettes/comparison-packages-mean_data_1-bcp.svg
deleted file mode 100644
index 82c3abb7..00000000
--- a/vignettes/comparison-packages-mean_data_1-bcp.svg
+++ /dev/null
@@ -1,1365 +0,0 @@
-
-
diff --git a/vignettes/comparison-packages-mean_data_1-mcp.svg b/vignettes/comparison-packages-mean_data_1-mcp.svg
deleted file mode 100644
index 75ec9634..00000000
--- a/vignettes/comparison-packages-mean_data_1-mcp.svg
+++ /dev/null
@@ -1,1264 +0,0 @@
-
-
diff --git a/vignettes/comparison-packages-mean_data_3-bcp.svg b/vignettes/comparison-packages-mean_data_3-bcp.svg
deleted file mode 100644
index eedf390e..00000000
--- a/vignettes/comparison-packages-mean_data_3-bcp.svg
+++ /dev/null
@@ -1,3373 +0,0 @@
-
-
diff --git a/vignettes/comparison-packages-mv_data_1-mcp.svg b/vignettes/comparison-packages-mv_data_1-mcp.svg
deleted file mode 100644
index 6971c29c..00000000
--- a/vignettes/comparison-packages-mv_data_1-mcp.svg
+++ /dev/null
@@ -1,2255 +0,0 @@
-
-
diff --git a/vignettes/comparison-packages-well_log-cpts.svg b/vignettes/comparison-packages-well_log-cpts.svg
deleted file mode 100644
index 0ae82362..00000000
--- a/vignettes/comparison-packages-well_log-cpts.svg
+++ /dev/null
@@ -1,5102 +0,0 @@
-
-
diff --git a/vignettes/comparison-packages-well_log-time.svg b/vignettes/comparison-packages-well_log-time.svg
deleted file mode 100644
index 0894f11e..00000000
--- a/vignettes/comparison-packages-well_log-time.svg
+++ /dev/null
@@ -1,497 +0,0 @@
-
-
diff --git a/vignettes/comparison-packages.Rmd b/vignettes/comparison-packages.Rmd
index cbe82e77..8d371d5e 100644
--- a/vignettes/comparison-packages.Rmd
+++ b/vignettes/comparison-packages.Rmd
@@ -9,18 +9,14 @@ vignette: >
%\VignetteEncoding{UTF-8}
---
-```{r, include = FALSE}
-knitr::opts_chunk$set(
- collapse = TRUE, comment = "#>", eval = TRUE, cache = FALSE,
- warning = FALSE, fig.width = 8, fig.height = 5
-)
-```
+
# Data setup
## Univariate mean change
-```{r data-setup-univariate-mean-change}
+
+```r
# Univariate mean change
set.seed(1)
p <- 1
@@ -33,9 +29,12 @@ mean_data_1 <- rbind(
plot.ts(mean_data_1)
```
+
+
## Univariate mean and/or variance change
-```{r data-setup-univariate-mean-and-or-variance-change}
+
+```r
# Univariate mean and/or variance change
set.seed(1)
p <- 1
@@ -51,9 +50,12 @@ mv_data_1 <- rbind(
plot.ts(mv_data_1)
```
+
+
## Multivariate mean change
-```{r data-setup-multivariate-mean-change}
+
+```r
# Multivariate mean change
set.seed(1)
p <- 3
@@ -66,9 +68,12 @@ mean_data_3 <- rbind(
plot.ts(mean_data_3)
```
+
+
## Multivariate mean and/or variance change
-```{r data-setup-multivariate-mean-and-or-variance-change}
+
+```r
# Multivariate mean and/or variance change
set.seed(1)
p <- 4
@@ -84,9 +89,12 @@ mv_data_3 <- rbind(
plot.ts(mv_data_3)
```
+
+
## Linear regression
-```{r data-setup-linear-regression}
+
+```r
# Linear regression
set.seed(1)
n <- 300
@@ -103,9 +111,12 @@ lm_data <- data.frame(y = y, x = x)
plot.ts(lm_data)
```
+
+
## Logistic regression
-```{r data-setup-logistic-regression}
+
+```r
# Logistic regression
set.seed(1)
n <- 500
@@ -121,9 +132,12 @@ binomial_data <- data.frame(y = y, x = x)
plot.ts(binomial_data)
```
+
+
## Poisson regression
-```{r data-setup-poisson-regression}
+
+```r
# Poisson regression
set.seed(1)
n <- 1100
@@ -140,12 +154,20 @@ y <- c(
poisson_data <- data.frame(y = y, x = x)
plot.ts(log(poisson_data$y))
+```
+
+
+
+```r
plot.ts(poisson_data[, -1])
```
+
+
## Lasso
-```{r data-setup-lasso}
+
+```r
# Lasso
set.seed(1)
n <- 480
@@ -170,9 +192,12 @@ lasso_data <- data.frame(y = y, x = x)
plot.ts(lasso_data[, seq_len(p_true + 1)])
```
+
+
## AR(3)
-```{r data-setup-ar3}
+
+```r
# AR(3)
set.seed(1)
n <- 1000
@@ -188,9 +213,12 @@ ar_data <- x[-seq_len(3)]
plot.ts(ar_data)
```
+
+
## GARCH(1, 1)
-```{r data-setup-garch11}
+
+```r
# GARCH(1, 1)
set.seed(1)
n <- 400
@@ -209,9 +237,12 @@ garch_data <- x[-1]
plot.ts(garch_data)
```
+
+
## VAR(2)
-```{r data-setup-var2}
+
+```r
# VAR(2)
set.seed(1)
n <- 800
@@ -230,635 +261,1846 @@ var_data <- x[-seq_len(2), ]
plot.ts(var_data)
```
+
+
# Univariate mean change
The true change points are 300 and 700.
Some methods are plotted due to the un-retrievable change points.
-```{r univariate-mean-change-fastcpd}
-(fastcpd_result <- fastcpd::fastcpd.mean(mean_data_1, r.progress = FALSE)@cp_set)
-```
-```{r univariate-mean-change-fastcpd-testthat, include = FALSE}
-testthat::expect_equal(fastcpd_result, c(300, 700), tolerance = 0.2)
+```r
+results[["mean_data_1"]][["fastcpd"]] <-
+ fastcpd::fastcpd.mean(mean_data_1, r.progress = FALSE)@cp_set
```
-```{r univariate-mean-change-CptNonPar}
-(CptNonPar_result <- CptNonPar::np.mojo(mean_data_1, G = floor(length(mean_data_1) / 6))$cpts)
+
+```r
+results[["mean_data_1"]][["fastcpd"]]
+#> [1] 300 700
```
-```{r univariate-mean-change-CptNonPar-testthat, include = FALSE}
-testthat::expect_equal(CptNonPar_result, c(300, 700), tolerance = 0.2)
+
+
+
+```r
+results[["mean_data_1"]][["CptNonPar"]] <-
+ CptNonPar::np.mojo(mean_data_1, G = floor(length(mean_data_1) / 6))$cpts
```
-```{r univariate-mean-change-strucchange, eval = FALSE}
-strucchange::breakpoints(y ~ 1, data = data.frame(y = mean_data_1))$breakpoints
+
+```r
+results[["mean_data_1"]][["CptNonPar"]]
#> [1] 300 700
```
-```{r univariate-mean-change-ecp, eval = FALSE}
-ecp::e.divisive(mean_data_1)$estimates
-#> [1] 1 301 701 1001
+
+
+
+```r
+results[["mean_data_1"]][["strucchange"]] <-
+ strucchange::breakpoints(y ~ 1, data = data.frame(y = mean_data_1))$breakpoints
```
-```{r univariate-mean-change-changepoint}
-(changepoint_result <- changepoint::cpt.mean(c(mean_data_1))@cpts)
+
+```r
+results[["mean_data_1"]][["strucchange"]]
+#> [1] 300 700
```
-```{r univariate-mean-change-changepoint-testthat, include = FALSE}
-testthat::expect_equal(changepoint_result, c(300, 1000), tolerance = 0.2)
+
+
+
+```r
+results[["mean_data_1"]][["ecp"]] <- ecp::e.divisive(mean_data_1)$estimates
```
-```{r univariate-mean-change-breakfast}
-(breakfast_result <- breakfast::breakfast(mean_data_1)$cptmodel.list[[6]]$cpts)
+
+```r
+results[["mean_data_1"]][["ecp"]]
+#> [1] 1 301 701 1001
```
-```{r univariate-mean-change-breakfast-testthat, include = FALSE}
-testthat::expect_equal(breakfast_result, c(300, 700), tolerance = 0.2)
+
+
+
+```r
+results[["mean_data_1"]][["changepoint"]] <-
+ changepoint::cpt.mean(c(mean_data_1))@cpts
```
-```{r univariate-mean-change-wbs}
-(wbs_result <- wbs::wbs(mean_data_1)$cpt$cpt.ic$mbic.penalty)
+
+```r
+results[["mean_data_1"]][["changepoint"]]
+#> [1] 300 1000
```
-```{r univariate-mean-change-wbs-testthat, include = FALSE}
-testthat::expect_equal(wbs_result, c(300, 700), tolerance = 0.2)
+
+
+
+```r
+results[["mean_data_1"]][["breakfast"]] <-
+ breakfast::breakfast(mean_data_1)$cptmodel.list[[6]]$cpts
```
-```{r univariate-mean-change-mosum, eval = FALSE}
-mosum::mosum(c(mean_data_1), G = 40)$cpts.info$cpts
+
+```r
+results[["mean_data_1"]][["breakfast"]]
#> [1] 300 700
```
-```{r univariate-mean-change-fpop}
-(fpop_result <- fpop::Fpop(mean_data_1, nrow(mean_data_1))$t.est)
-```
-```{r univariate-mean-change-fpop-testthat, include = FALSE}
-testthat::expect_equal(fpop_result, c(300, 700, 1000), tolerance = 0.2)
-```
-```{r univariate-mean-change-gfpop, eval = FALSE}
-gfpop::gfpop(
- data = mean_data_1,
- mygraph = gfpop::graph(
- penalty = 2 * log(nrow(mean_data_1)) * gfpop::sdDiff(mean_data_1) ^ 2,
- type = "updown"
- ),
- type = "mean"
-)$changepoints
-#> [1] 300 700 1000
+
+```r
+results[["mean_data_1"]][["wbs"]] <-
+ wbs::wbs(mean_data_1)$cpt$cpt.ic$mbic.penalty
```
-```{r univariate-mean-change-InspectChangepoint, eval = FALSE}
-invisible(
- suppressMessages(
- capture.output(
- result_InspectChangepoint <- InspectChangepoint::inspect(
- t(mean_data_1),
- threshold = InspectChangepoint::compute.threshold(
- nrow(mean_data_1), ncol(mean_data_1)
- )
- )
- )
- )
-)
-result_InspectChangepoint$changepoints[, "location"]
+
+```r
+results[["mean_data_1"]][["wbs"]]
#> [1] 300 700
```
-```{r univariate-mean-change-jointseg}
-(jointseg_result <- jointseg::jointSeg(mean_data_1, K = 2)$bestBkp)
-```
-```{r univariate-mean-change-jointseg-testthat, include = FALSE}
-testthat::expect_equal(jointseg_result, c(300, 700), tolerance = 0.2)
-```
-```{r univariate-mean-change-Rbeast, eval = FALSE}
-Rbeast::beast(
- mean_data_1, season = "none", print.progress = FALSE, quiet = TRUE
-)$trend$cp
-#> [1] 701 301 NaN NaN NaN NaN NaN NaN NaN NaN
-```
-```{r univariate-mean-change-stepR}
-(stepR_result <- stepR::stepFit(mean_data_1, alpha = 0.5)$rightEnd)
+```r
+results[["mean_data_1"]][["mosum"]] <-
+ mosum::mosum(c(mean_data_1), G = 40)$cpts.info$cpts
```
-```{r univariate-mean-change-stepR-testthat, include = FALSE}
-testthat::expect_equal(stepR_result, c(300, 700, 1000), tolerance = 0.2)
-```
-```{r univariate-mean-change-cpm}
-(cpm_result <- cpm::processStream(mean_data_1, cpmType = "Student")$changePoints)
+```r
+results[["mean_data_1"]][["mosum"]]
+#> [1] 300 700
```
-```{r univariate-mean-change-cpm-testthat, include = FALSE}
-testthat::expect_equal(cpm_result, c(299, 699), tolerance = 0.2)
-```
-```{r univariate-mean-change-segmented}
-(segmented_result <- segmented::stepmented(
- as.numeric(mean_data_1), npsi = 2
-)$psi[, "Est."])
-```
-```{r univariate-mean-change-segmented-testthat, include = FALSE}
-testthat::expect_equivalent(segmented_result, c(298, 699), tolerance = 0.2)
+
+```r
+results[["mean_data_1"]][["fpop"]] <-
+ fpop::Fpop(mean_data_1, nrow(mean_data_1))$t.est
```
-```{r univariate-mean-change-mcp, eval = FALSE}
-plot(
- mcp::mcp(
- list(y ~ 1, ~ 1, ~ 1),
- data = data.frame(y = mean_data_1, x = seq_len(nrow(mean_data_1))),
- par_x = "x"
- )
-)
+
+```r
+results[["mean_data_1"]][["fpop"]]
+#> [1] 300 700 1000
```
-
-```{r univariate-mean-change-not}
-plot(not::not(mean_data_1, contrast = "pcwsConstMean"))
+
+
+```r
+results[["mean_data_1"]][["gfpop"]] <-
+ gfpop::gfpop(
+ data = mean_data_1,
+ mygraph = gfpop::graph(
+ penalty = 2 * log(nrow(mean_data_1)) * gfpop::sdDiff(mean_data_1) ^ 2,
+ type = "updown"
+ ),
+ type = "mean"
+ )$changepoints
```
-```{r univariate-mean-change-bcp, eval = FALSE}
-plot(bcp::bcp(mean_data_1))
+
+```r
+results[["mean_data_1"]][["gfpop"]]
+#> [1] 300 700 1000
```
-
-# Univariate mean and/or variance change
-The true change points are 300, 700, 1000, 1300 and 1700.
-Some methods are plotted due to the un-retrievable change points.
-```{r univariate-mean-and-or-variance-change-fastcpd}
-(fastcpd_result <- fastcpd::fastcpd.mv(mv_data_1, r.progress = FALSE)@cp_set)
+```r
+results[["mean_data_1"]][["InspectChangepoint"]] <- InspectChangepoint::inspect(
+ t(mean_data_1),
+ threshold = InspectChangepoint::compute.threshold(
+ nrow(mean_data_1), ncol(mean_data_1)
+ )
+)$changepoints[, "location"]
```
-```{r univariate-mean-and-or-variance-change-fastcpd-testthat, include = FALSE}
-testthat::expect_equal(fastcpd_result, c(300, 700, 1001, 1300, 1700), tolerance = 0.2)
-```
-```{r univariate-mean-and-or-variance-change-ecp, eval = FALSE}
-ecp::e.divisive(mv_data_1)$estimates
-#> [1] 1 301 701 1001 1301 1701 2001
+```r
+results[["mean_data_1"]][["InspectChangepoint"]]
+#> [1] 300 700
```
-```{r univariate-mean-and-or-variance-change-changepoint}
-(changepoint_result <- changepoint::cpt.meanvar(c(mv_data_1))@cpts)
-```
-```{r univariate-mean-and-or-variance-change-changepoint-testthat, include = FALSE}
-testthat::expect_equal(changepoint_result, c(300, 2000), tolerance = 0.2)
-```
-```{r univariate-mean-and-or-variance-change-CptNonPar}
-(CptNonPar_result <- CptNonPar::np.mojo(mv_data_1, G = floor(length(mv_data_1) / 6))$cpts)
+
+```r
+results[["mean_data_1"]][["jointseg"]] <-
+ jointseg::jointSeg(mean_data_1, K = 2)$bestBkp
```
-```{r univariate-mean-and-or-variance-change-CptNonPar-testthat, include = FALSE}
-testthat::expect_equal(CptNonPar_result, c(333, 700, 1300), tolerance = 0.2)
+
+```r
+results[["mean_data_1"]][["jointseg"]]
+#> [1] 300 700
```
-```{r univariate-mean-and-or-variance-change-cpm}
-(cpm_result <- cpm::processStream(mv_data_1, cpmType = "GLR")$changePoints)
+
+
+
+```r
+results[["mean_data_1"]][["Rbeast"]] <-
+ Rbeast::beast(
+ mean_data_1, season = "none", print.progress = FALSE, quiet = TRUE
+ )$trend$cp
```
-```{r univariate-mean-and-or-variance-change-cpm-testthat, include = FALSE}
-testthat::expect_equal(cpm_result, c(293, 300, 403, 408, 618, 621, 696, 1000, 1021, 1024, 1293, 1300, 1417, 1693, 1700, 1981), tolerance = 0.2)
+
+```r
+results[["mean_data_1"]][["Rbeast"]]
+#> [1] 701 301 NaN NaN NaN NaN NaN NaN NaN NaN
```
-```{r univariate-mean-and-or-variance-change-InspectChangepoint, eval = FALSE}
-invisible(
- suppressMessages(
- capture.output(
- result_InspectChangepoint <- InspectChangepoint::inspect(
- t(mv_data_1),
- threshold = InspectChangepoint::compute.threshold(
- nrow(mv_data_1), ncol(mv_data_1)
- )
- )
- )
- )
-)
-result_InspectChangepoint$changepoints[, "location"]
-#> [1] 300 700 701 702 704 707 708 712 715 716 717 718 721 722
-#> [15] 723 726 727 729 731 732 734 736 740 742 744 746 748 750
-#> [29] 753 755 756 757 759 760 762 764 765 768 769 771 772 774
-#> [43] 776 777 784 785 786 789 791 792 794 797 798 799 801 802
-#> [57] 803 807 809 810 813 815 817 819 826 827 828 829 831 833
-#> [71] 835 836 837 838 840 841 842 843 845 848 849 852 854 860
-#> [85] 862 864 866 868 870 872 875 879 881 884 886 887 888 889
-#> [99] 896 897 898 899 901 903 904 905 906 909 912 913 915 917
-#> [113] 919 921 922 923 927 928 932 934 936 937 940 944 945 947
-#> [127] 948 949 951 956 958 959 961 962 963 964 966 967 968 972
-#> [141] 974 976 978 979 986 988 990 992 995 1000 1300 1700 1702 1703
-#> [155] 1704 1705 1708 1710 1712 1714 1716 1717 1718 1720 1721 1723 1725 1726
-#> [169] 1727 1729 1731 1733 1735 1736 1737 1739 1742 1745 1747 1748 1752 1754
-#> [183] 1756 1758 1759 1760 1766 1768 1770 1771 1773 1775 1778 1782 1784 1785
-#> [197] 1790 1792 1793 1795 1796 1797 1799 1800 1802 1803 1804 1805 1806 1807
-#> [211] 1808 1809 1821 1824 1825 1827 1828 1829 1833 1835 1837 1840 1841 1842
-#> [225] 1848 1849 1851 1852 1854 1855 1857 1859 1860 1862 1863 1865 1867 1868
-#> [239] 1876 1878 1879 1880 1882 1883 1884 1886 1887 1889 1894 1898 1899 1905
-#> [253] 1906 1907 1908 1909 1912 1919 1926 1927 1928 1930 1933 1934 1935 1936
-#> [267] 1938 1940 1944 1947 1950 1952 1954 1955 1956 1960 1962 1963 1965 1966
-#> [281] 1967 1969 1970 1974 1976 1977 1978 1980 1985 1987 1988 1990 1997 1998
-```
-
-```{r univariate-mean-and-or-variance-change-Rbeast, eval = FALSE}
-Rbeast::beast(
- mv_data_1, season = "none", print.progress = FALSE, quiet = TRUE
-)$trend$cp
-#> [1] 1794 1855 1986 1301 301 703 1981 1769 1860 1834
-```
-
-```{r univariate-mean-and-or-variance-change-mcp, eval = FALSE}
-plot(
- mcp::mcp(
- list(y ~ 1, ~ 1, ~ 1, ~ 1, ~ 1, ~ 1),
- data = data.frame(y = mv_data_1, x = seq_len(nrow(mv_data_1))),
- par_x = "x"
- )
-)
+
+
+
+```r
+results[["mean_data_1"]][["stepR"]] <-
+ stepR::stepFit(mean_data_1, alpha = 0.5)$rightEnd
```
-
-```{r univariate-mean-and-or-variance-change-not}
-plot(not::not(mv_data_1, contrast = "pcwsConstMeanVar"))
+```r
+results[["mean_data_1"]][["stepR"]]
+#> [1] 300 700 1000
```
-# Multivariate mean change
-The true change points are 300 and 700.
-Some methods are plotted due to the un-retrievable change points.
-```{r multivariate-mean-change-fastcpd}
-(fastcpd_result <- fastcpd::fastcpd.mean(mean_data_3, r.progress = FALSE)@cp_set)
+
+```r
+results[["mean_data_1"]][["cpm"]] <-
+ cpm::processStream(mean_data_1, cpmType = "Student")$changePoints
```
-```{r multivariate-mean-change-fastcpd-testthat, include = FALSE}
-testthat::expect_equal(fastcpd_result, c(300, 700), tolerance = 0.2)
+
+```r
+results[["mean_data_1"]][["cpm"]]
+#> [1] 299 699
```
-```{r multivariate-mean-change-CptNonPar}
-(CptNonPar_result <- CptNonPar::np.mojo(mean_data_3, G = floor(nrow(mean_data_3) / 6))$cpts)
+
+
+
+```r
+results[["mean_data_1"]][["segmented"]] <-
+ segmented::stepmented(
+ as.numeric(mean_data_1), npsi = 2
+ )$psi[, "Est."]
```
-```{r multivariate-mean-change-CptNonPar-testthat, include = FALSE}
-testthat::expect_equal(CptNonPar_result, c(300, 700), tolerance = 0.2)
+
+```r
+results[["mean_data_1"]][["segmented"]]
+#> psi1.index psi2.index
+#> 298.1981 699.1524
```
-```{r multivariate-mean-change-InspectChangepoint, eval = FALSE}
-invisible(
- suppressMessages(
- capture.output(
- result_InspectChangepoint <- InspectChangepoint::inspect(
- t(mean_data_3),
- threshold = InspectChangepoint::compute.threshold(
- nrow(mean_data_3), ncol(mean_data_3)
- )
- )
- )
- )
+
+
+
+```r
+results[["mean_data_1"]][["mcp"]] <- mcp::mcp(
+ list(y ~ 1, ~ 1, ~ 1),
+ data = data.frame(y = mean_data_1, x = seq_len(nrow(mean_data_1))),
+ par_x = "x"
)
-result_InspectChangepoint$changepoints[, "location"]
-#> [1] 300 700
```
-```{r multivariate-mean-change-jointseg}
-(jointseg_result <- jointseg::jointSeg(mean_data_3, K = 2)$bestBkp)
-```
-```{r multivariate-mean-change-jointseg-testthat, include = FALSE}
-testthat::expect_equal(jointseg_result, c(300, 700), tolerance = 0.2)
+```r
+if (requireNamespace("mcp", quietly = TRUE)) {
+ plot(results[["mean_data_1"]][["mcp"]])
+}
```
-```{r multivariate-mean-change-Rbeast, eval = FALSE}
-invisible(
- capture.output(
- result_Rbeast <- Rbeast::beast123(
- mean_data_3,
- metadata = list(whichDimIsTime = 1),
- season = "none"
- )
- )
-)
-result_Rbeast$trend$cp
-#> Warning message:
-#> In Rbeast::beast123(mean_data_3, metadata = list(whichDimIsTime = 1), :
-#> WARNING: If the input data is regular and ordered in time,the times of individual datapoints are determined fully by 'metadata$startTime' and 'metadata$deltaTime'. But startTime and deltaTime are missing and a default value 1 is used for both!
-#> [,1] [,2] [,3]
-#> [1,] 301 301 701
-#> [2,] 701 701 301
-#> [3,] NaN 225 NaN
-#> [4,] NaN 684 NaN
-#> [5,] NaN NaN NaN
-#> [6,] NaN NaN NaN
-#> [7,] NaN NaN NaN
-#> [8,] NaN NaN NaN
-#> [9,] NaN NaN NaN
-#> [10,] NaN NaN NaN
+
+
+
+```r
+results[["mean_data_1"]][["not"]] <-
+ not::not(mean_data_1, contrast = "pcwsConstMean")
```
-```{r multivariate-mean-change-strucchange, eval = FALSE}
-strucchange::breakpoints(
- cbind(y.1, y.2, y.3) ~ 1, data = data.frame(y = mean_data_3)
-)$breakpoints
-#> [1] 300 700
+
+```r
+if (requireNamespace("not", quietly = TRUE)) {
+ plot(results[["mean_data_1"]][["not"]])
+}
```
-```{r multivariate-mean-change-ecp, eval = FALSE}
-ecp::e.divisive(mean_data_3)$estimates
-#> [1] 1 301 701 1001
+
+
+
+```r
+results[["mean_data_1"]][["bcp"]] <- bcp::bcp(mean_data_1)
```
-```{r multivariate-mean-change-bcp, eval = FALSE}
-plot(bcp::bcp(mean_data_3))
+
+```r
+if (requireNamespace("bcp", quietly = TRUE)) {
+ plot(results[["mean_data_1"]][["bcp"]])
+}
```
-
+
-# Multivariate mean and/or variance change
+# Univariate mean and/or variance change
The true change points are 300, 700, 1000, 1300 and 1700.
Some methods are plotted due to the un-retrievable change points.
-```{r multivariate-mean-and-or-variance-change-fastcpd}
-(fastcpd_result <- fastcpd::fastcpd.mv(mv_data_3, r.progress = FALSE)@cp_set)
+
+```r
+results[["mv_data_1"]][["fastcpd"]] <-
+ fastcpd::fastcpd.mv(mv_data_1, r.progress = FALSE)@cp_set
```
-```{r multivariate-mean-and-or-variance-change-fastcpd-testthat, include = FALSE}
-testthat::expect_equal(fastcpd_result, c(300, 700, 1000, 1300, 1700), tolerance = 0.2)
+
+```r
+results[["mv_data_1"]][["fastcpd"]]
+#> [1] 300 700 1001 1300 1700
```
-```{r multivariate-mean-and-or-variance-change-ecp, eval = FALSE}
-ecp::e.divisive(mv_data_3)$estimates
-#> [1] 1 301 701 1001 1301 1701 2001
+
+
+
+```r
+results[["mv_data_1"]][["ecp"]] <- ecp::e.divisive(mv_data_1)$estimates
```
-```{r multivariate-mean-and-or-variance-change-InspectChangepoint, eval = FALSE}
-invisible(
- suppressMessages(
- capture.output(
- result_InspectChangepoint <- InspectChangepoint::inspect(
- t(mv_data_3),
- threshold = InspectChangepoint::compute.threshold(
- nrow(mv_data_3), ncol(mv_data_3)
- )
- )
- )
- )
-)
-result_InspectChangepoint$changepoints[, "location"]
-#> [1] 300 700 701 703 705 707 708 709 711 712 714 715 717 718
-#> [15] 720 721 723 724 726 727 729 731 733 734 736 737 739 740
-#> [29] 742 743 744 746 747 749 750 752 753 754 755 756 758 760
-#> [43] 762 763 765 766 767 769 770 772 773 774 775 777 779 780
-#> [57] 782 784 786 788 790 791 793 795 797 799 801 803 804 806
-#> [71] 808 809 810 811 813 814 816 817 818 820 821 823 825 827
-#> [85] 828 830 831 833 835 836 837 838 840 842 843 845 846 848
-#> [99] 849 850 852 853 854 855 856 858 859 860 862 863 865 866
-#> [113] 868 869 871 872 874 876 877 878 879 881 883 885 887 888
-#> [127] 889 891 893 894 895 897 898 900 901 903 904 906 908 909
-#> [141] 911 913 914 916 917 918 920 921 923 924 925 927 928 929
-#> [155] 931 932 934 936 937 938 939 941 942 943 945 946 947 949
-#> [169] 950 952 954 955 956 957 958 959 961 962 964 965 967 968
-#> [183] 970 972 973 974 975 977 979 981 982 984 985 986 987 988
-#> [197] 990 991 992 994 995 997 999 1000 1300 1700 1702 1703 1704 1705
-#> [211] 1706 1708 1709 1710 1712 1713 1714 1715 1717 1719 1721 1722 1723 1725
-#> [225] 1727 1729 1730 1732 1734 1735 1737 1738 1739 1741 1742 1744 1746 1748
-#> [239] 1750 1752 1753 1754 1755 1757 1758 1759 1761 1762 1763 1764 1766 1767
-#> [253] 1769 1770 1771 1773 1774 1775 1777 1779 1781 1782 1783 1785 1786 1788
-#> [267] 1789 1791 1793 1794 1796 1798 1800 1803 1804 1805 1806 1808 1809 1811
-#> [281] 1812 1814 1815 1817 1818 1819 1821 1822 1824 1825 1827 1828 1829 1831
-#> [295] 1833 1835 1836 1838 1839 1841 1843 1844 1846 1847 1848 1850 1851 1853
-#> [309] 1854 1856 1857 1858 1859 1860 1862 1863 1864 1865 1867 1869 1870 1872
-#> [323] 1873 1874 1876 1878 1879 1881 1882 1884 1885 1887 1889 1891 1893 1894
-#> [337] 1896 1898 1899 1900 1901 1902 1904 1906 1907 1909 1911 1913 1914 1916
-#> [351] 1917 1918 1919 1921 1923 1924 1925 1927 1928 1930 1932 1933 1935 1936
-#> [365] 1938 1939 1941 1942 1944 1946 1948 1950 1951 1952 1954 1956 1957 1959
-#> [379] 1961 1963 1965 1967 1968 1970 1972 1973 1974 1976 1977 1979 1981 1982
-#> [393] 1984 1985 1987 1989 1990 1992 1993 1995 1996 1998
-```
-
-```{r multivariate-mean-and-or-variance-change-Rbeast, eval = FALSE}
-invisible(
- capture.output(
- result_Rbeast <- Rbeast::beast123(
- mv_data_3,
- metadata = list(whichDimIsTime = 1),
- season = "none"
- )
- )
-)
-result_Rbeast$trend$cp
-#> Warning message:
-#> In Rbeast::beast123(mv_data_3, metadata = list(whichDimIsTime = 1), :
-#> WARNING: If the input data is regular and ordered in time,the times of individual datapoints are determined fully by 'metadata$startTime' and 'metadata$deltaTime'. But startTime and deltaTime are missing and a default value 1 is used for both!
-#> [,1] [,2] [,3] [,4]
-#> [1,] 701 301 301 710
-#> [2,] 1301 1301 1301 301
-#> [3,] 301 701 702 1301
-#> [4,] 814 1993 1829 1986
-#> [5,] 1968 767 1822 785
-#> [6,] 1994 781 810 774
-#> [7,] 761 884 845 1912
-#> [8,] 1873 755 798 792
-#> [9,] 1865 747 771 879
-#> [10,] 1962 924 1700 1919
+
+```r
+results[["mv_data_1"]][["ecp"]]
+#> [1] 1 301 701 1001 1301 1701 2001
```
-# Linear regression
-The true change points are 100 and 200.
-```{r linear-regression-fastcpd}
-(fastcpd_result <- fastcpd::fastcpd.lm(lm_data, r.progress = FALSE)@cp_set)
-```
-```{r linear-regression-fastcpd-testthat, include = FALSE}
-testthat::expect_equal(fastcpd_result, c(97, 201), tolerance = 0.2)
+```r
+results[["mv_data_1"]][["changepoint"]] <-
+ changepoint::cpt.meanvar(c(mv_data_1))@cpts
```
-```{r linear-regression-strucchange, eval = FALSE}
-strucchange::breakpoints(y ~ . - 1, data = lm_data)$breakpoints
-#> [1] 100 201
+
+```r
+results[["mv_data_1"]][["changepoint"]]
+#> [1] 300 2000
```
-```{r linear-regression-segmented}
-(segmented_result <- segmented::segmented(
- lm(
- y ~ . - 1,
- data.frame(y = lm_data$y, x = lm_data[, -1], index = seq_len(nrow(lm_data)))
- ),
- seg.Z = ~ index
-)$psi[, "Est."])
+
+
+
+```r
+results[["mv_data_1"]][["CptNonPar"]] <-
+ CptNonPar::np.mojo(mv_data_1, G = floor(length(mv_data_1) / 6))$cpts
```
-```{r linear-regression-segmented-testthat, include = FALSE}
-testthat::expect_equivalent(segmented_result, c(233), tolerance = 0.2)
+
+```r
+results[["mv_data_1"]][["CptNonPar"]]
+#> [1] 333 700 1300
```
-# Logistic regression
-The true change point is 300.
-```{r logistic-regression-fastcpd}
-(fastcpd_result <- fastcpd::fastcpd.binomial(binomial_data, r.progress = FALSE)@cp_set)
-```
-```{r logistic-regression-fastcpd-testthat, include = FALSE}
-testthat::expect_equal(fastcpd_result, 302, tolerance = 0.2)
+```r
+results[["mv_data_1"]][["cpm"]] <-
+ cpm::processStream(mv_data_1, cpmType = "GLR")$changePoints
```
-```{r logistic-regression-strucchange, eval = FALSE}
-strucchange::breakpoints(y ~ . - 1, data = binomial_data)$breakpoints
-#> [1] 297
+
+```r
+results[["mv_data_1"]][["cpm"]]
+#> [1] 293 300 403 408 618 621 696 1000 1021 1024 1293 1300 1417 1693 1700
+#> [16] 1981
```
-# Poisson regression
-The true change points are 500, 800 and 1000.
-```{r poisson-regression-fastcpd}
-(fastcpd_result <- fastcpd::fastcpd.poisson(poisson_data, r.progress = FALSE)@cp_set)
-```
-```{r poisson-regression-fastcpd-testthat, include = FALSE}
-testthat::expect_equal(fastcpd_result, c(498, 805, 1003), tolerance = 0.2)
+```r
+results[["mv_data_1"]][["InspectChangepoint"]] <- InspectChangepoint::inspect(
+ t(mv_data_1),
+ threshold = InspectChangepoint::compute.threshold(
+ nrow(mv_data_1), ncol(mv_data_1)
+ )
+)$changepoints[, "location"]
```
-```{r poisson-regression-strucchange, eval = FALSE}
-strucchange::breakpoints(y ~ . - 1, data = poisson_data)$breakpoints
-#> [1] 935
+
+```r
+results[["mv_data_1"]][["InspectChangepoint"]]
+#> [1] 300 700 701 702 704 707 708 712 715 716 717 718 721 722 723
+#> [16] 726 727 729 731 732 734 736 740 742 744 746 748 750 753 755
+#> [31] 756 757 759 760 762 764 765 766 768 769 771 772 774 776 777
+#> [46] 784 785 786 789 791 792 794 797 798 799 801 802 803 807 809
+#> [61] 810 813 815 817 819 826 827 828 829 831 833 835 836 837 838
+#> [76] 840 841 842 843 845 848 849 852 854 860 862 864 866 868 870
+#> [91] 872 875 879 881 884 886 887 888 889 896 897 898 899 901 903
+#> [106] 904 905 906 909 910 912 913 915 917 919 921 922 923 925 927
+#> [121] 928 932 934 936 937 940 944 945 947 948 949 951 956 958 959
+#> [136] 961 962 963 964 966 967 968 972 974 976 978 979 986 988 990
+#> [151] 992 995 998 1000 1300 1700 1702 1703 1704 1705 1708 1710 1712 1714 1716
+#> [166] 1717 1718 1720 1721 1723 1725 1726 1727 1729 1731 1733 1735 1736 1737 1739
+#> [181] 1742 1745 1747 1748 1752 1754 1756 1758 1759 1760 1766 1768 1770 1771 1773
+#> [196] 1775 1778 1782 1784 1785 1790 1792 1793 1795 1796 1797 1799 1800 1802 1803
+#> [211] 1804 1805 1806 1807 1808 1809 1813 1815 1816 1818 1821 1824 1825 1827 1828
+#> [226] 1829 1833 1835 1837 1840 1841 1842 1848 1849 1851 1852 1854 1855 1857 1859
+#> [241] 1860 1862 1863 1865 1867 1868 1876 1878 1879 1880 1882 1883 1884 1886 1887
+#> [256] 1889 1894 1898 1899 1905 1906 1907 1908 1909 1912 1919 1920 1921 1923 1924
+#> [271] 1926 1927 1928 1930 1933 1934 1935 1936 1938 1940 1941 1944 1947 1950 1952
+#> [286] 1954 1955 1956 1960 1962 1963 1965 1966 1967 1969 1970 1974 1976 1977 1978
+#> [301] 1980 1985 1987 1988 1990 1996 1997 1998
```
-# Lasso
-The true change points are 80, 200 and 320.
-```{r lasso-fastcpd}
-(fastcpd_result <- fastcpd::fastcpd.lasso(lasso_data, r.progress = FALSE)@cp_set)
-```
-```{r lasso-fastcpd-testthat, include = FALSE}
-testthat::expect_true(sum(fastcpd_result - c(79, 199, 320)) <= 1)
+```r
+results[["mv_data_1"]][["Rbeast"]] <-
+ Rbeast::beast(
+ mv_data_1, season = "none", print.progress = FALSE, quiet = TRUE
+ )$trend$cp
```
-```{r lasso-strucchange, eval = FALSE}
-strucchange::breakpoints(y ~ . - 1, data = lasso_data)$breakpoints,
-#> [1] 80 200 321
+
+```r
+results[["mv_data_1"]][["Rbeast"]]
+#> [1] 1855 1794 1301 1986 301 703 1981 1769 1860 709
```
-# AR(3)
-The true change point is 600.
-Some methods are plotted due to the un-retrievable change points.
-```{r ar3-fastcpd}
-(fastcpd_result <- fastcpd::fastcpd.ar(ar_data, 3, r.progress = FALSE)@cp_set)
+
+```r
+results[["mv_data_1"]][["mcp"]] <- mcp::mcp(
+ list(y ~ 1, ~ 1, ~ 1, ~ 1, ~ 1, ~ 1),
+ data = data.frame(y = mv_data_1, x = seq_len(nrow(mv_data_1))),
+ par_x = "x"
+)
```
-```{r ar3-fastcpd-testthat, include = FALSE}
-testthat::expect_equal(fastcpd_result, c(614), tolerance = 0.2)
+
+```r
+if (requireNamespace("mcp", quietly = TRUE)) {
+ plot(results[["mv_data_1"]][["mcp"]])
+}
```
-```{r ar3-CptNonPar}
-(CptNonPar_result <- CptNonPar::np.mojo(ar_data, G = floor(length(ar_data) / 6))$cpts)
+
+
+
+```r
+results[["mv_data_1"]][["not"]] <-
+ not::not(mv_data_1, contrast = "pcwsConstMeanVar")
```
-```{r ar3-CptNonPar-testthat, include = FALSE}
-testthat::expect_equal(CptNonPar_result, numeric(0), tolerance = 0.2)
+
+```r
+if (requireNamespace("not", quietly = TRUE)) {
+ plot(results[["mv_data_1"]][["not"]])
+}
```
-```{r ar3-segmented}
-(segmented_result <- segmented::segmented(
- lm(
- y ~ x + 1, data.frame(y = ar_data, x = seq_along(ar_data))
- ),
- seg.Z = ~ x
-)$psi[, "Est."])
+
+
```
+#> Press [enter] to continue
+```
+
+
+
+# Multivariate mean change
+
+The true change points are 300 and 700.
+Some methods are plotted due to the un-retrievable change points.
+
-```{r ar3-segmented-testthat, include = FALSE}
-testthat::expect_equivalent(segmented_result, c(690), tolerance = 0.2)
+```r
+results[["mean_data_3"]][["fastcpd"]] <-
+ fastcpd::fastcpd.mean(mean_data_3, r.progress = FALSE)@cp_set
```
-```{r ar3-mcp, eval = FALSE}
-plot(
- mcp::mcp(
- list(y ~ 1 + ar(3), ~ 0 + ar(3)),
- data = data.frame(y = ar_data, x = seq_along(ar_data)),
- par_x = "x"
- )
-)
+
+```r
+results[["mean_data_3"]][["fastcpd"]]
+#> [1] 300 700
```
-
-# GARCH(1, 1)
-The true change point is 200.
-```{r garch11-fastcpd}
-(fastcpd_result <- fastcpd::fastcpd.garch(garch_data, c(1, 1), r.progress = FALSE)@cp_set)
+```r
+results[["mean_data_3"]][["CptNonPar"]] <-
+ CptNonPar::np.mojo(mean_data_3, G = floor(nrow(mean_data_3) / 6))$cpts
```
-```{r garch11-fastcpd-testthat, include = FALSE}
-testthat::expect_equal(fastcpd_result, c(205), tolerance = 0.2)
-```
-```{r garch11-CptNonPar}
-(CptNonPar_result <- CptNonPar::np.mojo(garch_data, G = floor(length(garch_data) / 6))$cpts)
+```r
+results[["mean_data_3"]][["CptNonPar"]]
+#> [1] 300 700
```
-```{r garch11-CptNonPar-testthat, include = FALSE}
-testthat::expect_equal(CptNonPar_result, c(206), tolerance = 0.2)
+
+
+
+```r
+results[["mean_data_3"]][["InspectChangepoint"]] <- InspectChangepoint::inspect(
+ t(mean_data_3),
+ threshold = InspectChangepoint::compute.threshold(
+ nrow(mean_data_3), ncol(mean_data_3)
+ )
+)$changepoints[, "location"]
```
-```{r garch11-strucchange, eval = FALSE}
-strucchange::breakpoints(x ~ 1, data = data.frame(x = garch_data))$breakpoints
-#> [1] NA
+
+```r
+results[["mean_data_3"]][["InspectChangepoint"]]
+#> [1] 300 700
```
-# VAR(2)
-The true change points is 500.
-```{r var2-fastcpd}
-(fastcpd_result <- fastcpd::fastcpd.var(
- var_data, 2, cost_adjustment = NULL, r.progress = FALSE
-)@cp_set)
+
+```r
+results[["mean_data_3"]][["jointseg"]] <-
+ jointseg::jointSeg(mean_data_3, K = 2)$bestBkp
```
-```{r var2-fastcpd-testthat, include = FALSE}
-testthat::expect_equal(fastcpd_result, c(500), tolerance = 0.2)
+
+```r
+results[["mean_data_3"]][["jointseg"]]
+#> [1] 300 700
```
-```{r var2-VARDetect}
-(VARDetect_result <- VARDetect::tbss(var_data)$cp)
+
+
+
+```r
+results[["mean_data_3"]][["Rbeast"]] <-
+ Rbeast::beast123(
+ mean_data_3,
+ metadata = list(whichDimIsTime = 1),
+ season = "none"
+ )$trend$cp
```
-```{r var2-VARDetect-testthat, include = FALSE}
-testthat::expect_equal(VARDetect_result, c(501), tolerance = 0.2)
+
+```r
+results[["mean_data_3"]][["Rbeast"]]
+#> [,1] [,2] [,3]
+#> [1,] 301 701 301
+#> [2,] 701 301 701
+#> [3,] 142 117 926
+#> [4,] 305 NaN NaN
+#> [5,] 705 NaN NaN
+#> [6,] 694 NaN NaN
+#> [7,] NaN NaN NaN
+#> [8,] NaN NaN NaN
+#> [9,] NaN NaN NaN
+#> [10,] NaN NaN NaN
```
-# Detection comparison using `well_log`
-```{r detection-comparison-well-log, eval = FALSE}
-well_log <- well_log[well_log > 1e5]
-result <- list(
- fastcpd = fastcpd.mean(well_log, trim = 0.003)@cp_set,
- changepoint = changepoint::cpt.mean(well_log)@cpts,
- CptNonPar =
- CptNonPar::np.mojo(well_log, G = floor(length(well_log) / 6))$cpts,
- strucchange = strucchange::breakpoints(
+
+```r
+results[["mean_data_3"]][["strucchange"]] <-
+ strucchange::breakpoints(
+ cbind(y.1, y.2, y.3) ~ 1, data = data.frame(y = mean_data_3)
+ )$breakpoints
+```
+
+
+```r
+results[["mean_data_3"]][["strucchange"]]
+#> [1] 300 700
+```
+
+
+
+
+```r
+results[["mean_data_3"]][["ecp"]] <- ecp::e.divisive(mean_data_3)$estimates
+```
+
+
+```r
+results[["mean_data_3"]][["ecp"]]
+#> [1] 1 301 701 1001
+```
+
+
+
+
+```r
+results[["mean_data_3"]][["bcp"]] <- bcp::bcp(mean_data_3)
+```
+
+
+```r
+if (requireNamespace("bcp", quietly = TRUE)) {
+ plot(results[["mean_data_3"]][["bcp"]])
+}
+```
+
+
+
+# Multivariate mean and/or variance change
+
+The true change points are 300, 700, 1000, 1300 and 1700.
+Some methods are plotted due to the un-retrievable change points.
+
+
+```r
+results[["mv_data_3"]][["fastcpd"]] <-
+ fastcpd::fastcpd.mv(mv_data_3, r.progress = FALSE)@cp_set
+```
+
+
+```r
+results[["mv_data_3"]][["fastcpd"]]
+#> [1] 300 700 1000 1300 1700
+```
+
+
+
+
+```r
+results[["mv_data_3"]][["ecp"]] <- ecp::e.divisive(mv_data_3)$estimates
+```
+
+
+```r
+results[["mv_data_3"]][["ecp"]]
+#> [1] 1 301 701 1001 1301 1701 2001
+```
+
+
+
+
+```r
+results[["mv_data_3"]][["InspectChangepoint"]] <- InspectChangepoint::inspect(
+ t(mv_data_3),
+ threshold = InspectChangepoint::compute.threshold(
+ nrow(mv_data_3), ncol(mv_data_3)
+ )
+)$changepoints[, "location"]
+```
+
+
+```r
+results[["mv_data_3"]][["InspectChangepoint"]]
+#> [1] 300 700 701 703 705 707 708 709 711 712 714 715 717 718 720
+#> [16] 721 723 724 726 727 729 731 733 734 736 737 739 740 742 743
+#> [31] 744 746 747 749 750 752 753 754 755 756 758 760 762 763 765
+#> [46] 766 767 769 770 772 773 774 775 777 779 780 782 784 786 788
+#> [61] 790 791 793 795 797 799 801 803 804 806 809 810 811 813 814
+#> [76] 816 817 818 820 821 823 825 827 828 830 831 833 835 836 837
+#> [91] 838 840 842 843 845 846 848 849 850 852 853 854 855 856 858
+#> [106] 859 860 862 863 865 866 868 869 871 872 874 876 877 878 879
+#> [121] 881 883 885 887 888 889 891 893 894 895 897 898 900 901 903
+#> [136] 904 906 908 909 911 913 914 916 917 918 920 921 923 924 925
+#> [151] 927 928 929 931 932 934 936 937 938 939 941 942 943 945 946
+#> [166] 947 949 950 952 954 955 956 957 958 959 961 962 964 965 967
+#> [181] 968 970 972 973 974 975 977 979 981 982 984 985 986 987 988
+#> [196] 990 991 992 994 995 997 999 1000 1300 1700 1702 1703 1704 1705 1706
+#> [211] 1708 1709 1710 1712 1713 1714 1715 1717 1719 1721 1722 1723 1725 1727 1729
+#> [226] 1730 1732 1734 1735 1737 1738 1739 1741 1742 1744 1746 1748 1750 1752 1753
+#> [241] 1754 1755 1757 1758 1759 1761 1762 1763 1764 1766 1767 1769 1770 1771 1773
+#> [256] 1774 1775 1777 1779 1781 1782 1783 1785 1786 1788 1789 1791 1793 1794 1796
+#> [271] 1798 1800 1803 1804 1805 1806 1808 1809 1811 1812 1814 1815 1817 1818 1819
+#> [286] 1821 1822 1824 1825 1827 1828 1829 1831 1833 1835 1836 1838 1839 1841 1843
+#> [301] 1844 1846 1847 1848 1850 1851 1853 1854 1856 1857 1858 1859 1860 1862 1863
+#> [316] 1864 1865 1867 1869 1870 1872 1873 1874 1876 1878 1879 1881 1882 1884 1885
+#> [331] 1887 1889 1891 1893 1894 1896 1898 1899 1900 1901 1902 1904 1906 1907 1909
+#> [346] 1911 1913 1914 1916 1917 1918 1919 1921 1923 1924 1925 1927 1928 1930 1932
+#> [361] 1933 1935 1936 1938 1939 1941 1942 1944 1946 1948 1950 1951 1952 1954 1956
+#> [376] 1957 1959 1961 1963 1965 1967 1968 1970 1972 1973 1974 1976 1977 1979 1981
+#> [391] 1982 1984 1985 1987 1989 1990 1992 1993 1995 1996 1998
+```
+
+
+
+
+```r
+results[["mv_data_3"]][["Rbeast"]] <-
+ Rbeast::beast123(
+ mv_data_3,
+ metadata = list(whichDimIsTime = 1),
+ season = "none"
+ )$trend$cp
+```
+
+
+```r
+results[["mv_data_3"]][["Rbeast"]]
+#> [,1] [,2] [,3] [,4]
+#> [1,] 701 1301 301 1301
+#> [2,] 1301 301 1301 710
+#> [3,] 301 701 1829 301
+#> [4,] 1968 1993 702 886
+#> [5,] 1994 884 1822 1975
+#> [6,] 814 755 810 1915
+#> [7,] 1962 781 845 778
+#> [8,] 1978 767 1738 1985
+#> [9,] 1870 747 1754 792
+#> [10,] 1843 722 771 953
+```
+
+
+
+# Linear regression
+
+The true change points are 100 and 200.
+
+
+```r
+results[["lm_data"]][["fastcpd"]] <-
+ fastcpd::fastcpd.lm(lm_data, r.progress = FALSE)@cp_set
+```
+
+
+```r
+results[["lm_data"]][["fastcpd"]]
+#> [1] 97 201
+```
+
+
+
+
+```r
+results[["lm_data"]][["strucchange"]] <-
+ strucchange::breakpoints(y ~ . - 1, data = lm_data)$breakpoints
+```
+
+
+```r
+results[["lm_data"]][["strucchange"]]
+#> [1] 100 201
+```
+
+
+
+
+```r
+results[["lm_data"]][["segmented"]] <-
+ segmented::segmented(
+ lm(
+ y ~ . - 1, data.frame(y = lm_data$y, x = lm_data[, -1], index = seq_len(nrow(lm_data)))
+ ),
+ seg.Z = ~ index
+ )$psi[, "Est."]
+```
+
+
+```r
+results[["lm_data"]][["segmented"]]
+#> [1] 233
+```
+
+
+
+# Logistic regression
+
+The true change point is 300.
+
+
+```r
+results[["binomial_data"]][["fastcpd"]] <-
+ fastcpd::fastcpd.binomial(binomial_data, r.progress = FALSE)@cp_set
+```
+
+
+```r
+results[["binomial_data"]][["fastcpd"]]
+#> [1] 302
+```
+
+
+
+
+```r
+results[["binomial_data"]][["strucchange"]] <-
+ strucchange::breakpoints(y ~ . - 1, data = binomial_data)$breakpoints
+```
+
+
+```r
+results[["binomial_data"]][["strucchange"]]
+#> [1] 297
+```
+
+
+
+# Poisson regression
+
+The true change points are 500, 800 and 1000.
+
+
+```r
+results[["poisson_data"]][["fastcpd"]] <-
+ fastcpd::fastcpd.poisson(poisson_data, r.progress = FALSE)@cp_set
+```
+
+
+```r
+results[["poisson_data"]][["fastcpd"]]
+#> [1] 498 805 1003
+```
+
+
+
+
+```r
+results[["poisson_data"]][["strucchange"]] <-
+ strucchange::breakpoints(y ~ . - 1, data = poisson_data)$breakpoints
+```
+
+
+```r
+results[["poisson_data"]][["strucchange"]]
+#> [1] 935
+```
+
+
+
+# Lasso
+
+The true change points are 80, 200 and 320.
+
+
+```r
+results[["lasso_data"]][["fastcpd"]] <-
+ fastcpd::fastcpd.lasso(lasso_data, r.progress = FALSE)@cp_set
+```
+
+
+```r
+results[["lasso_data"]][["fastcpd"]]
+#> [1] 79 199 321
+```
+
+
+
+
+```r
+results[["lasso_data"]][["strucchange"]] <-
+ strucchange::breakpoints(y ~ . - 1, data = lasso_data)$breakpoints
+```
+
+
+```r
+results[["lasso_data"]][["strucchange"]]
+#> [1] 80 200 321
+```
+
+
+
+# AR(3)
+
+The true change point is 600.
+Some methods are plotted due to the un-retrievable change points.
+
+
+```r
+results[["ar_data"]][["fastcpd"]] <-
+ fastcpd::fastcpd.ar(ar_data, 3, r.progress = FALSE)@cp_set
+```
+
+
+```r
+results[["ar_data"]][["fastcpd"]]
+#> [1] 614
+```
+
+
+
+
+```r
+results[["ar_data"]][["CptNonPar"]] <-
+ CptNonPar::np.mojo(ar_data, G = floor(length(ar_data) / 6))$cpts
+```
+
+
+```r
+results[["ar_data"]][["CptNonPar"]]
+#> numeric(0)
+```
+
+
+
+
+```r
+results[["ar_data"]][["segmented"]] <-
+ segmented::segmented(
+ lm(
+ y ~ x + 1, data.frame(y = ar_data, x = seq_along(ar_data))
+ ),
+ seg.Z = ~ x
+ )$psi[, "Est."]
+```
+
+
+```r
+results[["ar_data"]][["segmented"]]
+#> [1] 690
+```
+
+
+
+
+```r
+results[["ar_data"]][["mcp"]] <-
+ mcp::mcp(
+ list(y ~ 1 + ar(3), ~ 0 + ar(3)),
+ data = data.frame(y = ar_data, x = seq_along(ar_data)),
+ par_x = "x"
+ )
+```
+
+
+```r
+if (requireNamespace("mcp", quietly = TRUE)) {
+ plot(results[["ar_data"]][["mcp"]])
+}
+```
+
+
+
+# GARCH(1, 1)
+
+The true change point is 200.
+
+
+```r
+results[["garch_data"]][["fastcpd"]] <-
+ fastcpd::fastcpd.garch(garch_data, c(1, 1), r.progress = FALSE)@cp_set
+```
+
+
+```r
+results[["garch_data"]][["fastcpd"]]
+#> [1] 205
+```
+
+
+
+
+```r
+results[["garch_data"]][["CptNonPar"]] <-
+ CptNonPar::np.mojo(garch_data, G = floor(length(garch_data) / 6))$cpts
+```
+
+
+```r
+results[["garch_data"]][["CptNonPar"]]
+#> [1] 206
+```
+
+
+
+
+```r
+results[["garch_data"]][["strucchange"]] <-
+ strucchange::breakpoints(x ~ 1, data = data.frame(x = garch_data))$breakpoints
+```
+
+
+```r
+results[["garch_data"]][["strucchange"]]
+#> [1] NA
+```
+
+
+
+# VAR(2)
+
+The true change points is 500.
+
+
+```r
+results[["var_data"]][["fastcpd"]] <-
+ fastcpd::fastcpd.var(var_data, 2, r.progress = FALSE)@cp_set
+```
+
+
+```r
+results[["var_data"]][["fastcpd"]]
+#> [1] 500
+```
+
+
+
+
+```r
+results[["var_data"]][["VARDetect"]] <- VARDetect::tbss(var_data)$cp
+```
+
+
+```r
+results[["var_data"]][["VARDetect"]]
+#> [1] 501
+```
+
+
+
+# Detection comparison using `well_log`
+
+
+```r
+well_log <- fastcpd::well_log
+```
+
+
+```r
+well_log <- well_log[well_log > 1e5]
+
+results[["well_log"]] <- list(
+ fastcpd = fastcpd.mean(well_log, trim = 0.003)@cp_set,
+ changepoint = changepoint::cpt.mean(well_log)@cpts,
+ CptNonPar =
+ CptNonPar::np.mojo(well_log, G = floor(length(well_log) / 6))$cpts,
+ strucchange = strucchange::breakpoints(
+ y ~ 1, data = data.frame(y = well_log)
+ )$breakpoints,
+ ecp = ecp::e.divisive(matrix(well_log))$estimates,
+ breakfast = breakfast::breakfast(well_log)$cptmodel.list[[6]]$cpts,
+ wbs = wbs::wbs(well_log)$cpt$cpt.ic$mbic.penalty,
+ mosum = mosum::mosum(c(well_log), G = 40)$cpts.info$cpts,
+ # fpop = fpop::Fpop(well_log, length(well_log))$t.est, # meaningless
+ gfpop = gfpop::gfpop(
+ data = well_log,
+ mygraph = gfpop::graph(
+ penalty = 2 * log(length(well_log)) * gfpop::sdDiff(well_log) ^ 2,
+ type = "updown"
+ ),
+ type = "mean"
+ )$changepoints,
+ InspectChangepoint = InspectChangepoint::inspect(
+ well_log,
+ threshold = InspectChangepoint::compute.threshold(length(well_log), 1)
+ )$changepoints[, "location"],
+ jointseg = jointseg::jointSeg(well_log, K = 12)$bestBkp,
+ Rbeast = Rbeast::beast(
+ well_log, season = "none", print.progress = FALSE, quiet = TRUE
+ )$trend$cp,
+ stepR = stepR::stepFit(well_log, alpha = 0.5)$rightEnd
+)
+```
+
+
+```r
+results[["well_log"]]
+#> $fastcpd
+#> [1] 12 572 704 776 1021 1057 1198 1347 1406 1502 1660 1842 2023 2202 2384
+#> [16] 2445 2507 2567 2749 2926 3094 3107 3509 3622 3709 3820 3976
+#>
+#> $changepoint
+#> [1] 2738 3989
+#>
+#> $CptNonPar
+#> [1] 1021 1681 2022 2738
+#>
+#> $strucchange
+#> [1] 1057 1660 2568 3283
+#>
+#> $ecp
+#> [1] 1 33 315 435 567 705 803 1026 1058 1348 1503 1662 1843 2024 2203
+#> [16] 2386 2446 2508 2569 2745 2780 2922 3073 3136 3252 3465 3500 3554 3623 3710
+#> [31] 3821 3868 3990
+#>
+#> $breakfast
+#> [1] 33 310 434 572 704 779 1021 1057 1347 1502 1659 1842 2021 2032 2202
+#> [16] 2384 2446 2507 2567 2747 2779 2926 3094 3106 3125 3283 3464 3499 3622 3709
+#> [31] 3806 3835 3848 3877 3896 3976
+#>
+#> $wbs
+#> [1] 2568 1057 1661 1842 2385 2023 2445 1502 2744 6 2507 1021 3709 3820 1402
+#> [16] 434 1408 3131 3976 3509 3622 776 1197 3094 704 3104 1347 314 2921 3251
+#> [31] 3464 3848 3906 2779 1663 3636 60 3904 2202 566 12 3639 1200 7 1671
+#> [46] 706
+#>
+#> $mosum
+#> [1] 6 434 704 1017 1057 1325 1502 1661 1842 2023 2385 2445 2507 2567 2744
+#> [16] 3060 3438 3509 3610 3697 3820 3867 3976
+#>
+#> $gfpop
+#> [1] 6 7 8 12 314 434 556 560 704 776 1021 1057 1197 1200 1347
+#> [16] 1364 1405 1407 1491 1502 1661 1842 2023 2385 2445 2507 2567 2664 2747 2752
+#> [31] 2921 3094 3104 3125 3251 3464 3499 3622 3709 3820 3976 3989
+#>
+#> $InspectChangepoint
+#> [1] 6 8 15 71 314 434 521 704 706 776 830 1021 1057 1059 1347
+#> [16] 1402 1405 1408 1412 1502 1520 1659 1661 1663 1694 1842 1848 2022 2202 2384
+#> [31] 2387 2445 2507 2567 2568 2673 2738 2746 2752 2779 2921 3131 3251 3464 3509
+#> [46] 3609 3658 3709 3806 3976
+#>
+#> $jointseg
+#> [1] 6 1021 1057 1502 1661 1842 2022 2384 2445 2507 2568 2738
+#>
+#> $Rbeast
+#> [1] 1058 1662 7 1022 2447 2386 2023 1503 2745 3710
+#>
+#> $stepR
+#> [1] 7 14 314 434 566 704 776 1021 1057 1197 1200 1347 1405 1407 1502
+#> [16] 1661 1694 1842 2023 2202 2385 2445 2507 2567 2747 2752 2921 3094 3104 3125
+#> [31] 3251 3464 3499 3609 3658 3709 3820 3867 3905 3976 3989
+```
+
+
+```r
+package_list <- sort(names(results[["well_log"]]), decreasing = TRUE)
+comparison_table <- NULL
+for (package_index in seq_along(package_list)) {
+ package <- package_list[[package_index]]
+ comparison_table <- rbind(
+ comparison_table,
+ data.frame(
+ change_point = results[["well_log"]][[package]],
+ package = package,
+ y_offset = (package_index - 1) * 1000
+ )
+ )
+}
+
+most_selected <- sort(table(comparison_table$change_point), decreasing = TRUE)
+most_selected <- sort(as.numeric(names(most_selected[most_selected >= 4])))
+for (i in seq_len(length(most_selected) - 1)) {
+ if (most_selected[i + 1] - most_selected[i] < 2) {
+ most_selected[i] <- NA
+ most_selected[i + 1] <- most_selected[i + 1] - 0.5
+ }
+}
+(most_selected <- most_selected[!is.na(most_selected)])
+#> [1] 6.5 314.0 434.0 704.0 776.0 1021.0 1057.0 1347.0 1502.0 1661.0
+#> [11] 1842.0 2023.0 2202.0 2384.5 2445.0 2507.0 2567.5 2738.0 2921.0 3094.0
+#> [21] 3251.0 3464.0 3509.0 3622.0 3709.0 3820.0 3976.0
+```
+
+
+```r
+if (requireNamespace("ggplot2", quietly = TRUE)) {
+ ggplot2::ggplot() +
+ ggplot2::geom_point(
+ data = data.frame(x = seq_along(well_log), y = c(well_log)),
+ ggplot2::aes(x = x, y = y)
+ ) +
+ ggplot2::geom_vline(
+ xintercept = most_selected,
+ color = "black",
+ linetype = "dashed",
+ alpha = 0.2
+ ) +
+ ggplot2::geom_point(
+ data = comparison_table,
+ ggplot2::aes(x = change_point, y = 50000 + y_offset, color = package),
+ shape = 17,
+ size = 1.9
+ ) +
+ ggplot2::geom_hline(
+ data = comparison_table,
+ ggplot2::aes(yintercept = 50000 + y_offset, color = package),
+ linetype = "dashed",
+ alpha = 0.1
+ ) +
+ ggplot2::coord_cartesian(
+ ylim = c(50000 - 500, max(well_log) + 1000),
+ xlim = c(-200, length(well_log) + 200),
+ expand = FALSE
+ ) +
+ ggplot2::theme(
+ panel.background = ggplot2::element_blank(),
+ panel.border = ggplot2::element_rect(colour = "black", fill = NA),
+ panel.grid.major = ggplot2::element_blank(),
+ panel.grid.minor = ggplot2::element_blank()
+ ) +
+ ggplot2::xlab(NULL) + ggplot2::ylab(NULL)
+}
+```
+
+
+
+# Time comparison using `well_log`
+
+
+```r
+results[["microbenchmark"]] <- microbenchmark::microbenchmark(
+ fastcpd = fastcpd::fastcpd.mean(well_log, trim = 0.003, r.progress = FALSE),
+ changepoint = changepoint::cpt.mean(well_log, method = "PELT"),
+ CptNonPar = CptNonPar::np.mojo(well_log, G = floor(length(well_log) / 6)),
+ strucchange =
+ strucchange::breakpoints(y ~ 1, data = data.frame(y = well_log)),
+ ecp = ecp::e.divisive(matrix(well_log)),
+ breakfast = breakfast::breakfast(well_log),
+ wbs = wbs::wbs(well_log),
+ mosum = mosum::mosum(c(well_log), G = 40),
+ fpop = fpop::Fpop(well_log, nrow(well_log)),
+ gfpop = gfpop::gfpop(
+ data = well_log,
+ mygraph = gfpop::graph(
+ penalty = 2 * log(length(well_log)) * gfpop::sdDiff(well_log) ^ 2,
+ type = "updown"
+ ),
+ type = "mean"
+ ),
+ InspectChangepoint = InspectChangepoint::inspect(
+ well_log,
+ threshold = InspectChangepoint::compute.threshold(length(well_log), 1)
+ ),
+ jointseg = jointseg::jointSeg(well_log, K = 12),
+ Rbeast = Rbeast::beast(
+ well_log, season = "none", print.progress = FALSE, quiet = TRUE
+ ),
+ stepR = stepR::stepFit(well_log, alpha = 0.5),
+ not = not::not(well_log, contrast = "pcwsConstMean"),
+ times = 10
+)
+```
+
+
+```r
+results[["microbenchmark"]]
+#> Unit: milliseconds
+#> expr min lq mean median
+#> fastcpd 6.257120e+01 6.696175e+01 7.183964e+01 7.168669e+01
+#> changepoint 3.205076e+01 3.305076e+01 4.025595e+01 3.774400e+01
+#> CptNonPar 1.875995e+04 2.014073e+04 2.244910e+04 2.224506e+04
+#> strucchange 6.359889e+04 6.409690e+04 6.651444e+04 6.565439e+04
+#> ecp 7.505232e+05 7.895889e+05 8.168647e+05 8.170952e+05
+#> breakfast 9.606171e+03 9.819381e+03 1.052611e+04 1.048197e+04
+#> wbs 1.163237e+02 1.180537e+02 1.263443e+02 1.254133e+02
+#> mosum 1.059973e+00 1.147898e+00 2.175931e+00 1.409396e+00
+#> fpop 2.604002e+00 3.676689e+00 4.859312e+00 4.308841e+00
+#> gfpop 6.057383e+01 6.595090e+01 7.277643e+01 7.212125e+01
+#> InspectChangepoint 1.584727e+02 2.171007e+02 2.414921e+02 2.424665e+02
+#> jointseg 2.104403e+01 2.284920e+01 3.577861e+01 2.708083e+01
+#> Rbeast 6.545716e+02 6.913542e+02 7.993305e+02 7.504324e+02
+#> stepR 2.740003e+01 3.020269e+01 3.149624e+01 3.188329e+01
+#> not 9.715867e+01 1.000287e+02 1.101258e+02 1.074475e+02
+#> uq max neval
+#> 7.689381e+01 8.332375e+01 10
+#> 4.669485e+01 5.904173e+01 10
+#> 2.380236e+04 2.705807e+04 10
+#> 6.597650e+04 7.757488e+04 10
+#> 8.259550e+05 8.834288e+05 10
+#> 1.103289e+04 1.197235e+04 10
+#> 1.276892e+02 1.455261e+02 10
+#> 1.957312e+00 8.705062e+00 10
+#> 5.107175e+00 1.085365e+01 10
+#> 7.742432e+01 8.794816e+01 10
+#> 3.060519e+02 3.145695e+02 10
+#> 4.362760e+01 9.516090e+01 10
+#> 8.539695e+02 1.163607e+03 10
+#> 3.294534e+01 3.530981e+01 10
+#> 1.147504e+02 1.331989e+02 10
+```
+
+
+```r
+if (requireNamespace("ggplot2", quietly = TRUE) && requireNamespace("microbenchmark", quietly = TRUE)) {
+ ggplot2::autoplot(results[["microbenchmark"]])
+}
+```
+
+
+
+# Notes
+
+This document is generated by the following code:
+
+```shell
+R -e 'knitr::knit("vignettes/comparison-packages.Rmd.original", output = "vignettes/comparison-packages.Rmd")' && rm -rf vignettes/figure && mv -f figure vignettes
+```
+
+# Related issues
+
+- `mosum`: [[#4](https://github.com/doccstat/fastcpd/issues/4)].
+- `mcp`: [link 1](https://github.com/doccstat/fastcpd/actions/runs/7272426093/job/19814531091)
+ and [link 2](https://github.com/doccstat/fastcpd/actions/runs/7272426093/job/19814531284).
+- `bcp`: [[#5](https://github.com/doccstat/fastcpd/issues/5)].
+- `gfpop`: [[#10](https://github.com/doccstat/fastcpd/issues/10)].
+
+# Acknowledgements
+
+- [Dr. Vito Muggeo](https://www.unipa.it/persone/docenti/m/vito.muggeo),
+ author of the `segmented` package for the tips about the piece-wise
+ constant function.
+
+# Appendix: all code snippets
+
+
+
+
+```r
+knitr::opts_chunk$set(
+ collapse = TRUE, comment = "#>", eval = TRUE, cache = FALSE,
+ warning = FALSE, fig.width = 8, fig.height = 5
+)
+
+if (requireNamespace("microbenchmark", quietly = TRUE)) {
+ library(microbenchmark)
+}
+
+if (file.exists("comparison-packages-results.RData")) {
+ # Available at https://pcloud.xingchi.li/comparison-packages-results.RData
+ load("comparison-packages-results.RData")
+} else {
+ results <- list()
+}
+# Univariate mean change
+set.seed(1)
+p <- 1
+mean_data_1 <- rbind(
+ mvtnorm::rmvnorm(300, mean = rep(0, p), sigma = diag(100, p)),
+ mvtnorm::rmvnorm(400, mean = rep(50, p), sigma = diag(100, p)),
+ mvtnorm::rmvnorm(300, mean = rep(2, p), sigma = diag(100, p))
+)
+
+plot.ts(mean_data_1)
+# Univariate mean and/or variance change
+set.seed(1)
+p <- 1
+mv_data_1 <- rbind(
+ mvtnorm::rmvnorm(300, mean = rep(0, p), sigma = diag(1, p)),
+ mvtnorm::rmvnorm(400, mean = rep(10, p), sigma = diag(1, p)),
+ mvtnorm::rmvnorm(300, mean = rep(0, p), sigma = diag(100, p)),
+ mvtnorm::rmvnorm(300, mean = rep(0, p), sigma = diag(1, p)),
+ mvtnorm::rmvnorm(400, mean = rep(10, p), sigma = diag(1, p)),
+ mvtnorm::rmvnorm(300, mean = rep(10, p), sigma = diag(100, p))
+)
+
+plot.ts(mv_data_1)
+# Multivariate mean change
+set.seed(1)
+p <- 3
+mean_data_3 <- rbind(
+ mvtnorm::rmvnorm(300, mean = rep(0, p), sigma = diag(100, p)),
+ mvtnorm::rmvnorm(400, mean = rep(50, p), sigma = diag(100, p)),
+ mvtnorm::rmvnorm(300, mean = rep(2, p), sigma = diag(100, p))
+)
+
+plot.ts(mean_data_3)
+# Multivariate mean and/or variance change
+set.seed(1)
+p <- 4
+mv_data_3 <- rbind(
+ mvtnorm::rmvnorm(300, mean = rep(0, p), sigma = diag(1, p)),
+ mvtnorm::rmvnorm(400, mean = rep(10, p), sigma = diag(1, p)),
+ mvtnorm::rmvnorm(300, mean = rep(0, p), sigma = diag(100, p)),
+ mvtnorm::rmvnorm(300, mean = rep(0, p), sigma = diag(1, p)),
+ mvtnorm::rmvnorm(400, mean = rep(10, p), sigma = diag(1, p)),
+ mvtnorm::rmvnorm(300, mean = rep(10, p), sigma = diag(100, p))
+)
+
+plot.ts(mv_data_3)
+# Linear regression
+set.seed(1)
+n <- 300
+p <- 4
+x <- mvtnorm::rmvnorm(n, rep(0, p), diag(p))
+theta_0 <- rbind(c(1, 3.2, -1, 0), c(-1, -0.5, 2.5, -2), c(0.8, 0, 1, 2))
+y <- c(
+ x[1:100, ] %*% theta_0[1, ] + rnorm(100, 0, 3),
+ x[101:200, ] %*% theta_0[2, ] + rnorm(100, 0, 3),
+ x[201:n, ] %*% theta_0[3, ] + rnorm(100, 0, 3)
+)
+lm_data <- data.frame(y = y, x = x)
+
+plot.ts(lm_data)
+# Logistic regression
+set.seed(1)
+n <- 500
+p <- 4
+x <- mvtnorm::rmvnorm(n, rep(0, p), diag(p))
+theta <- rbind(rnorm(p, 0, 1), rnorm(p, 2, 1))
+y <- c(
+ rbinom(300, 1, 1 / (1 + exp(-x[1:300, ] %*% theta[1, ]))),
+ rbinom(200, 1, 1 / (1 + exp(-x[301:n, ] %*% theta[2, ])))
+)
+binomial_data <- data.frame(y = y, x = x)
+
+plot.ts(binomial_data)
+# Poisson regression
+set.seed(1)
+n <- 1100
+p <- 3
+x <- mvtnorm::rmvnorm(n, rep(0, p), diag(p))
+delta <- rnorm(p)
+theta_0 <- c(1, 0.3, -1)
+y <- c(
+ rpois(500, exp(x[1:500, ] %*% theta_0)),
+ rpois(300, exp(x[501:800, ] %*% (theta_0 + delta))),
+ rpois(200, exp(x[801:1000, ] %*% theta_0)),
+ rpois(100, exp(x[1001:1100, ] %*% (theta_0 - delta)))
+)
+poisson_data <- data.frame(y = y, x = x)
+
+plot.ts(log(poisson_data$y))
+plot.ts(poisson_data[, -1])
+# Lasso
+set.seed(1)
+n <- 480
+p_true <- 6
+p <- 50
+x <- mvtnorm::rmvnorm(n, rep(0, p), diag(p))
+theta_0 <- rbind(
+ runif(p_true, -5, -2),
+ runif(p_true, -3, 3),
+ runif(p_true, 2, 5),
+ runif(p_true, -5, 5)
+)
+theta_0 <- cbind(theta_0, matrix(0, ncol = p - p_true, nrow = 4))
+y <- c(
+ x[1:80, ] %*% theta_0[1, ] + rnorm(80, 0, 1),
+ x[81:200, ] %*% theta_0[2, ] + rnorm(120, 0, 1),
+ x[201:320, ] %*% theta_0[3, ] + rnorm(120, 0, 1),
+ x[321:n, ] %*% theta_0[4, ] + rnorm(160, 0, 1)
+)
+lasso_data <- data.frame(y = y, x = x)
+
+plot.ts(lasso_data[, seq_len(p_true + 1)])
+# AR(3)
+set.seed(1)
+n <- 1000
+x <- rep(0, n + 3)
+for (i in 1:600) {
+ x[i + 3] <- 0.6 * x[i + 2] - 0.2 * x[i + 1] + 0.1 * x[i] + rnorm(1, 0, 3)
+}
+for (i in 601:1000) {
+ x[i + 3] <- 0.3 * x[i + 2] + 0.4 * x[i + 1] + 0.2 * x[i] + rnorm(1, 0, 3)
+}
+ar_data <- x[-seq_len(3)]
+
+plot.ts(ar_data)
+# GARCH(1, 1)
+set.seed(1)
+n <- 400
+sigma_2 <- rep(1, n + 1)
+x <- rep(0, n + 1)
+for (i in seq_len(200)) {
+ sigma_2[i + 1] <- 20 + 0.5 * x[i]^2 + 0.1 * sigma_2[i]
+ x[i + 1] <- rnorm(1, 0, sqrt(sigma_2[i + 1]))
+}
+for (i in 201:400) {
+ sigma_2[i + 1] <- 1 + 0.1 * x[i]^2 + 0.5 * sigma_2[i]
+ x[i + 1] <- rnorm(1, 0, sqrt(sigma_2[i + 1]))
+}
+garch_data <- x[-1]
+
+plot.ts(garch_data)
+# VAR(2)
+set.seed(1)
+n <- 800
+p <- 2
+theta_1 <- matrix(c(-0.3, 0.6, -0.5, 0.4, 0.2, 0.2, 0.2, -0.2), nrow = p)
+theta_2 <- matrix(c(0.3, -0.4, 0.1, -0.5, -0.5, -0.2, -0.5, 0.2), nrow = p)
+x <- matrix(0, n + 2, p)
+for (i in 1:500) {
+ x[i + 2, ] <- theta_1 %*% c(x[i + 1, ], x[i, ]) + rnorm(p, 0, 1)
+}
+for (i in 501:n) {
+ x[i + 2, ] <- theta_2 %*% c(x[i + 1, ], x[i, ]) + rnorm(p, 0, 1)
+}
+var_data <- x[-seq_len(2), ]
+
+plot.ts(var_data)
+results[["mean_data_1"]][["fastcpd"]] <-
+ fastcpd::fastcpd.mean(mean_data_1, r.progress = FALSE)@cp_set
+results[["mean_data_1"]][["fastcpd"]]
+testthat::expect_equal(results[["mean_data_1"]][["fastcpd"]], c(300, 700), tolerance = 0.2)
+results[["mean_data_1"]][["CptNonPar"]] <-
+ CptNonPar::np.mojo(mean_data_1, G = floor(length(mean_data_1) / 6))$cpts
+results[["mean_data_1"]][["CptNonPar"]]
+testthat::expect_equal(results[["mean_data_1"]][["CptNonPar"]], c(300, 700), tolerance = 0.2)
+results[["mean_data_1"]][["strucchange"]] <-
+ strucchange::breakpoints(y ~ 1, data = data.frame(y = mean_data_1))$breakpoints
+results[["mean_data_1"]][["strucchange"]]
+testthat::expect_equal(results[["mean_data_1"]][["strucchange"]], c(300, 700), tolerance = 0.2)
+results[["mean_data_1"]][["ecp"]] <- ecp::e.divisive(mean_data_1)$estimates
+results[["mean_data_1"]][["ecp"]]
+testthat::expect_equal(results[["mean_data_1"]][["ecp"]], c(1, 301, 701, 1001), tolerance = 0.2)
+results[["mean_data_1"]][["changepoint"]] <-
+ changepoint::cpt.mean(c(mean_data_1))@cpts
+results[["mean_data_1"]][["changepoint"]]
+testthat::expect_equal(results[["mean_data_1"]][["changepoint"]], c(300, 1000), tolerance = 0.2)
+results[["mean_data_1"]][["breakfast"]] <-
+ breakfast::breakfast(mean_data_1)$cptmodel.list[[6]]$cpts
+results[["mean_data_1"]][["breakfast"]]
+testthat::expect_equal(results[["mean_data_1"]][["breakfast"]], c(300, 700), tolerance = 0.2)
+results[["mean_data_1"]][["wbs"]] <-
+ wbs::wbs(mean_data_1)$cpt$cpt.ic$mbic.penalty
+results[["mean_data_1"]][["wbs"]]
+testthat::expect_equal(results[["mean_data_1"]][["wbs"]], c(300, 700), tolerance = 0.2)
+results[["mean_data_1"]][["mosum"]] <-
+ mosum::mosum(c(mean_data_1), G = 40)$cpts.info$cpts
+results[["mean_data_1"]][["mosum"]]
+testthat::expect_equal(results[["mean_data_1"]][["mosum"]], c(300, 700), tolerance = 0.2)
+results[["mean_data_1"]][["fpop"]] <-
+ fpop::Fpop(mean_data_1, nrow(mean_data_1))$t.est
+results[["mean_data_1"]][["fpop"]]
+testthat::expect_equal(results[["mean_data_1"]][["fpop"]], c(300, 700, 1000), tolerance = 0.2)
+results[["mean_data_1"]][["gfpop"]] <-
+ gfpop::gfpop(
+ data = mean_data_1,
+ mygraph = gfpop::graph(
+ penalty = 2 * log(nrow(mean_data_1)) * gfpop::sdDiff(mean_data_1) ^ 2,
+ type = "updown"
+ ),
+ type = "mean"
+ )$changepoints
+results[["mean_data_1"]][["gfpop"]]
+testthat::expect_equal(results[["mean_data_1"]][["gfpop"]], c(300, 700, 1000), tolerance = 0.2)
+results[["mean_data_1"]][["InspectChangepoint"]] <- InspectChangepoint::inspect(
+ t(mean_data_1),
+ threshold = InspectChangepoint::compute.threshold(
+ nrow(mean_data_1), ncol(mean_data_1)
+ )
+)$changepoints[, "location"]
+results[["mean_data_1"]][["InspectChangepoint"]]
+testthat::expect_equal(results[["mean_data_1"]][["InspectChangepoint"]], c(300, 700), tolerance = 0.2)
+results[["mean_data_1"]][["jointseg"]] <-
+ jointseg::jointSeg(mean_data_1, K = 2)$bestBkp
+results[["mean_data_1"]][["jointseg"]]
+testthat::expect_equal(results[["mean_data_1"]][["jointseg"]], c(300, 700), tolerance = 0.2)
+results[["mean_data_1"]][["Rbeast"]] <-
+ Rbeast::beast(
+ mean_data_1, season = "none", print.progress = FALSE, quiet = TRUE
+ )$trend$cp
+results[["mean_data_1"]][["Rbeast"]]
+testthat::expect_equal(results[["mean_data_1"]][["Rbeast"]], c(701, 301, NaN, NaN, NaN, NaN, NaN, NaN, NaN, NaN), tolerance = 0.2)
+results[["mean_data_1"]][["stepR"]] <-
+ stepR::stepFit(mean_data_1, alpha = 0.5)$rightEnd
+results[["mean_data_1"]][["stepR"]]
+testthat::expect_equal(results[["mean_data_1"]][["stepR"]], c(300, 700, 1000), tolerance = 0.2)
+results[["mean_data_1"]][["cpm"]] <-
+ cpm::processStream(mean_data_1, cpmType = "Student")$changePoints
+results[["mean_data_1"]][["cpm"]]
+testthat::expect_equal(results[["mean_data_1"]][["cpm"]], c(299, 699), tolerance = 0.2)
+results[["mean_data_1"]][["segmented"]] <-
+ segmented::stepmented(
+ as.numeric(mean_data_1), npsi = 2
+ )$psi[, "Est."]
+results[["mean_data_1"]][["segmented"]]
+testthat::expect_equal(results[["mean_data_1"]][["segmented"]], c(298, 699), ignore_attr = TRUE, tolerance = 0.2)
+results[["mean_data_1"]][["mcp"]] <- mcp::mcp(
+ list(y ~ 1, ~ 1, ~ 1),
+ data = data.frame(y = mean_data_1, x = seq_len(nrow(mean_data_1))),
+ par_x = "x"
+)
+if (requireNamespace("mcp", quietly = TRUE)) {
+ plot(results[["mean_data_1"]][["mcp"]])
+}
+results[["mean_data_1"]][["not"]] <-
+ not::not(mean_data_1, contrast = "pcwsConstMean")
+if (requireNamespace("not", quietly = TRUE)) {
+ plot(results[["mean_data_1"]][["not"]])
+}
+results[["mean_data_1"]][["bcp"]] <- bcp::bcp(mean_data_1)
+if (requireNamespace("bcp", quietly = TRUE)) {
+ plot(results[["mean_data_1"]][["bcp"]])
+}
+results[["mv_data_1"]][["fastcpd"]] <-
+ fastcpd::fastcpd.mv(mv_data_1, r.progress = FALSE)@cp_set
+results[["mv_data_1"]][["fastcpd"]]
+testthat::expect_equal(results[["mv_data_1"]][["fastcpd"]], c(300, 700, 1001, 1300, 1700), tolerance = 0.2)
+results[["mv_data_1"]][["ecp"]] <- ecp::e.divisive(mv_data_1)$estimates
+results[["mv_data_1"]][["ecp"]]
+testthat::expect_equal(results[["mv_data_1"]][["ecp"]], c(1, 301, 701, 1001, 1301, 1701, 2001), tolerance = 0.2)
+results[["mv_data_1"]][["changepoint"]] <-
+ changepoint::cpt.meanvar(c(mv_data_1))@cpts
+results[["mv_data_1"]][["changepoint"]]
+testthat::expect_equal(results[["mv_data_1"]][["changepoint"]], c(300, 2000), tolerance = 0.2)
+results[["mv_data_1"]][["CptNonPar"]] <-
+ CptNonPar::np.mojo(mv_data_1, G = floor(length(mv_data_1) / 6))$cpts
+results[["mv_data_1"]][["CptNonPar"]]
+testthat::expect_equal(results[["mv_data_1"]][["CptNonPar"]], c(333, 700, 1300), tolerance = 0.2)
+results[["mv_data_1"]][["cpm"]] <-
+ cpm::processStream(mv_data_1, cpmType = "GLR")$changePoints
+results[["mv_data_1"]][["cpm"]]
+testthat::expect_equal(results[["mv_data_1"]][["cpm"]], c(293, 300, 403, 408, 618, 621, 696, 1000, 1021, 1024, 1293, 1300, 1417, 1693, 1700, 1981), tolerance = 0.2)
+results[["mv_data_1"]][["InspectChangepoint"]] <- InspectChangepoint::inspect(
+ t(mv_data_1),
+ threshold = InspectChangepoint::compute.threshold(
+ nrow(mv_data_1), ncol(mv_data_1)
+ )
+)$changepoints[, "location"]
+results[["mv_data_1"]][["InspectChangepoint"]]
+testthat::expect_equal(results[["mv_data_1"]][["InspectChangepoint"]], c(
+ 300, 700, 701, 702, 704, 707, 708, 712, 715, 716, 717, 718,
+ 721, 722, 723, 726, 727, 729, 731, 732, 734, 736, 740, 742,
+ 744, 746, 748, 750, 753, 755, 756, 757, 759, 760, 762, 764,
+ 765, 766, 768, 769, 771, 772, 774, 776, 777, 784, 785, 786,
+ 789, 791, 792, 794, 797, 798, 799, 801, 802, 803, 807, 809,
+ 810, 813, 815, 817, 819, 826, 827, 828, 829, 831, 833, 835,
+ 836, 837, 838, 840, 841, 842, 843, 845, 848, 849, 852, 854,
+ 860, 862, 864, 866, 868, 870, 872, 875, 879, 881, 884, 886,
+ 887, 888, 889, 896, 897, 898, 899, 901, 903, 904, 905, 906,
+ 909, 910, 912, 913, 915, 917, 919, 921, 922, 923, 925, 927,
+ 928, 932, 934, 936, 937, 940, 944, 945, 947, 948, 949, 951,
+ 956, 958, 959, 961, 962, 963, 964, 966, 967, 968, 972, 974,
+ 976, 978, 979, 986, 988, 990, 992, 995, 998, 1000, 1300, 1700,
+ 1702, 1703, 1704, 1705, 1708, 1710, 1712, 1714, 1716, 1717, 1718, 1720,
+ 1721, 1723, 1725, 1726, 1727, 1729, 1731, 1733, 1735, 1736, 1737, 1739,
+ 1742, 1745, 1747, 1748, 1752, 1754, 1756, 1758, 1759, 1760, 1766, 1768,
+ 1770, 1771, 1773, 1775, 1778, 1782, 1784, 1785, 1790, 1792, 1793, 1795,
+ 1796, 1797, 1799, 1800, 1802, 1803, 1804, 1805, 1806, 1807, 1808, 1809,
+ 1813, 1815, 1816, 1818, 1821, 1824, 1825, 1827, 1828, 1829, 1833, 1835,
+ 1837, 1840, 1841, 1842, 1848, 1849, 1851, 1852, 1854, 1855, 1857, 1859,
+ 1860, 1862, 1863, 1865, 1867, 1868, 1876, 1878, 1879, 1880, 1882, 1883,
+ 1884, 1886, 1887, 1889, 1894, 1898, 1899, 1905, 1906, 1907, 1908, 1909,
+ 1912, 1919, 1920, 1921, 1923, 1924, 1926, 1927, 1928, 1930, 1933, 1934,
+ 1935, 1936, 1938, 1940, 1941, 1944, 1947, 1950, 1952, 1954, 1955, 1956,
+ 1960, 1962, 1963, 1965, 1966, 1967, 1969, 1970, 1974, 1976, 1977, 1978,
+ 1980, 1985, 1987, 1988, 1990, 1996, 1997, 1998
+), tolerance = 0.2)
+results[["mv_data_1"]][["Rbeast"]] <-
+ Rbeast::beast(
+ mv_data_1, season = "none", print.progress = FALSE, quiet = TRUE
+ )$trend$cp
+results[["mv_data_1"]][["Rbeast"]]
+testthat::expect_equal(results[["mv_data_1"]][["Rbeast"]], c(1855, 1794, 1301, 1986, 301, 703, 1981, 1769, 1860, 709), tolerance = 0.2)
+results[["mv_data_1"]][["mcp"]] <- mcp::mcp(
+ list(y ~ 1, ~ 1, ~ 1, ~ 1, ~ 1, ~ 1),
+ data = data.frame(y = mv_data_1, x = seq_len(nrow(mv_data_1))),
+ par_x = "x"
+)
+if (requireNamespace("mcp", quietly = TRUE)) {
+ plot(results[["mv_data_1"]][["mcp"]])
+}
+results[["mv_data_1"]][["not"]] <-
+ not::not(mv_data_1, contrast = "pcwsConstMeanVar")
+if (requireNamespace("not", quietly = TRUE)) {
+ plot(results[["mv_data_1"]][["not"]])
+}
+results[["mean_data_3"]][["fastcpd"]] <-
+ fastcpd::fastcpd.mean(mean_data_3, r.progress = FALSE)@cp_set
+results[["mean_data_3"]][["fastcpd"]]
+testthat::expect_equal(results[["mean_data_3"]][["fastcpd"]], c(300, 700), tolerance = 0.2)
+results[["mean_data_3"]][["CptNonPar"]] <-
+ CptNonPar::np.mojo(mean_data_3, G = floor(nrow(mean_data_3) / 6))$cpts
+results[["mean_data_3"]][["CptNonPar"]]
+testthat::expect_equal(results[["mean_data_3"]][["CptNonPar"]], c(300, 700), tolerance = 0.2)
+results[["mean_data_3"]][["InspectChangepoint"]] <- InspectChangepoint::inspect(
+ t(mean_data_3),
+ threshold = InspectChangepoint::compute.threshold(
+ nrow(mean_data_3), ncol(mean_data_3)
+ )
+)$changepoints[, "location"]
+results[["mean_data_3"]][["InspectChangepoint"]]
+testthat::expect_equal(results[["mean_data_3"]][["InspectChangepoint"]], c(300, 700), tolerance = 0.2)
+results[["mean_data_3"]][["jointseg"]] <-
+ jointseg::jointSeg(mean_data_3, K = 2)$bestBkp
+results[["mean_data_3"]][["jointseg"]]
+testthat::expect_equal(results[["mean_data_3"]][["jointseg"]], c(300, 700), tolerance = 0.2)
+results[["mean_data_3"]][["Rbeast"]] <-
+ Rbeast::beast123(
+ mean_data_3,
+ metadata = list(whichDimIsTime = 1),
+ season = "none"
+ )$trend$cp
+results[["mean_data_3"]][["Rbeast"]]
+testthat::expect_equal(results[["mean_data_3"]][["Rbeast"]], matrix(c(
+ 301, 701, 301,
+ 701, 301, 701,
+ 142, 117, 926,
+ 305, NaN, NaN,
+ 705, NaN, NaN,
+ 694, NaN, NaN,
+ NaN, NaN, NaN,
+ NaN, NaN, NaN,
+ NaN, NaN, NaN,
+ NaN, NaN, NaN
+), nrow = 10, ncol = 3, byrow = TRUE), tolerance = 0.2)
+results[["mean_data_3"]][["strucchange"]] <-
+ strucchange::breakpoints(
+ cbind(y.1, y.2, y.3) ~ 1, data = data.frame(y = mean_data_3)
+ )$breakpoints
+results[["mean_data_3"]][["strucchange"]]
+testthat::expect_equal(results[["mean_data_3"]][["strucchange"]], c(300, 700), tolerance = 0.2)
+results[["mean_data_3"]][["ecp"]] <- ecp::e.divisive(mean_data_3)$estimates
+results[["mean_data_3"]][["ecp"]]
+testthat::expect_equal(results[["mean_data_3"]][["ecp"]], c(1, 301, 701, 1001), tolerance = 0.2)
+results[["mean_data_3"]][["bcp"]] <- bcp::bcp(mean_data_3)
+if (requireNamespace("bcp", quietly = TRUE)) {
+ plot(results[["mean_data_3"]][["bcp"]])
+}
+results[["mv_data_3"]][["fastcpd"]] <-
+ fastcpd::fastcpd.mv(mv_data_3, r.progress = FALSE)@cp_set
+results[["mv_data_3"]][["fastcpd"]]
+testthat::expect_equal(results[["mv_data_3"]][["fastcpd"]], c(300, 700, 1000, 1300, 1700), tolerance = 0.2)
+results[["mv_data_3"]][["ecp"]] <- ecp::e.divisive(mv_data_3)$estimates
+results[["mv_data_3"]][["ecp"]]
+testthat::expect_equal(results[["mv_data_3"]][["ecp"]], c(1, 301, 701, 1001, 1301, 1701, 2001), tolerance = 0.2)
+results[["mv_data_3"]][["InspectChangepoint"]] <- InspectChangepoint::inspect(
+ t(mv_data_3),
+ threshold = InspectChangepoint::compute.threshold(
+ nrow(mv_data_3), ncol(mv_data_3)
+ )
+)$changepoints[, "location"]
+results[["mv_data_3"]][["InspectChangepoint"]]
+testthat::expect_equal(results[["mv_data_3"]][["InspectChangepoint"]], c(
+ 300, 700, 701, 703, 705, 707, 708, 709, 711, 712, 714, 715,
+ 717, 718, 720, 721, 723, 724, 726, 727, 729, 731, 733, 734,
+ 736, 737, 739, 740, 742, 743, 744, 746, 747, 749, 750, 752,
+ 753, 754, 755, 756, 758, 760, 762, 763, 765, 766, 767, 769,
+ 770, 772, 773, 774, 775, 777, 779, 780, 782, 784, 786, 788,
+ 790, 791, 793, 795, 797, 799, 801, 803, 804, 806, 809, 810,
+ 811, 813, 814, 816, 817, 818, 820, 821, 823, 825, 827, 828,
+ 830, 831, 833, 835, 836, 837, 838, 840, 842, 843, 845, 846,
+ 848, 849, 850, 852, 853, 854, 855, 856, 858, 859, 860, 862,
+ 863, 865, 866, 868, 869, 871, 872, 874, 876, 877, 878, 879,
+ 881, 883, 885, 887, 888, 889, 891, 893, 894, 895, 897, 898,
+ 900, 901, 903, 904, 906, 908, 909, 911, 913, 914, 916, 917,
+ 918, 920, 921, 923, 924, 925, 927, 928, 929, 931, 932, 934,
+ 936, 937, 938, 939, 941, 942, 943, 945, 946, 947, 949, 950,
+ 952, 954, 955, 956, 957, 958, 959, 961, 962, 964, 965, 967,
+ 968, 970, 972, 973, 974, 975, 977, 979, 981, 982, 984, 985,
+ 986, 987, 988, 990, 991, 992, 994, 995, 997, 999, 1000, 1300,
+ 1700, 1702, 1703, 1704, 1705, 1706, 1708, 1709, 1710, 1712, 1713, 1714,
+ 1715, 1717, 1719, 1721, 1722, 1723, 1725, 1727, 1729, 1730, 1732, 1734,
+ 1735, 1737, 1738, 1739, 1741, 1742, 1744, 1746, 1748, 1750, 1752, 1753,
+ 1754, 1755, 1757, 1758, 1759, 1761, 1762, 1763, 1764, 1766, 1767, 1769,
+ 1770, 1771, 1773, 1774, 1775, 1777, 1779, 1781, 1782, 1783, 1785, 1786,
+ 1788, 1789, 1791, 1793, 1794, 1796, 1798, 1800, 1803, 1804, 1805, 1806,
+ 1808, 1809, 1811, 1812, 1814, 1815, 1817, 1818, 1819, 1821, 1822, 1824,
+ 1825, 1827, 1828, 1829, 1831, 1833, 1835, 1836, 1838, 1839, 1841, 1843,
+ 1844, 1846, 1847, 1848, 1850, 1851, 1853, 1854, 1856, 1857, 1858, 1859,
+ 1860, 1862, 1863, 1864, 1865, 1867, 1869, 1870, 1872, 1873, 1874, 1876,
+ 1878, 1879, 1881, 1882, 1884, 1885, 1887, 1889, 1891, 1893, 1894, 1896,
+ 1898, 1899, 1900, 1901, 1902, 1904, 1906, 1907, 1909, 1911, 1913, 1914,
+ 1916, 1917, 1918, 1919, 1921, 1923, 1924, 1925, 1927, 1928, 1930, 1932,
+ 1933, 1935, 1936, 1938, 1939, 1941, 1942, 1944, 1946, 1948, 1950, 1951,
+ 1952, 1954, 1956, 1957, 1959, 1961, 1963, 1965, 1967, 1968, 1970, 1972,
+ 1973, 1974, 1976, 1977, 1979, 1981, 1982, 1984, 1985, 1987, 1989, 1990,
+ 1992, 1993, 1995, 1996, 1998
+), tolerance = 0.2)
+results[["mv_data_3"]][["Rbeast"]] <-
+ Rbeast::beast123(
+ mv_data_3,
+ metadata = list(whichDimIsTime = 1),
+ season = "none"
+ )$trend$cp
+results[["mv_data_3"]][["Rbeast"]]
+testthat::expect_equal(results[["mv_data_3"]][["Rbeast"]], matrix(c(
+ 701, 1301, 301, 1301,
+ 1301, 301, 1301, 710,
+ 301, 701, 1829, 301,
+ 1968, 1993, 702, 886,
+ 1994, 884, 1822, 1975,
+ 814, 755, 810, 1915,
+ 1962, 781, 845, 778,
+ 1978, 767, 1738, 1985,
+ 1870, 747, 1754, 792,
+ 1843, 722, 771, 953
+), nrow = 10, ncol = 4, byrow = TRUE), tolerance = 0.2)
+results[["lm_data"]][["fastcpd"]] <-
+ fastcpd::fastcpd.lm(lm_data, r.progress = FALSE)@cp_set
+results[["lm_data"]][["fastcpd"]]
+testthat::expect_equal(results[["lm_data"]][["fastcpd"]], c(97, 201), tolerance = 0.2)
+results[["lm_data"]][["strucchange"]] <-
+ strucchange::breakpoints(y ~ . - 1, data = lm_data)$breakpoints
+results[["lm_data"]][["strucchange"]]
+testthat::expect_equal(results[["lm_data"]][["strucchange"]], c(100, 201), tolerance = 0.2)
+results[["lm_data"]][["segmented"]] <-
+ segmented::segmented(
+ lm(
+ y ~ . - 1, data.frame(y = lm_data$y, x = lm_data[, -1], index = seq_len(nrow(lm_data)))
+ ),
+ seg.Z = ~ index
+ )$psi[, "Est."]
+results[["lm_data"]][["segmented"]]
+testthat::expect_equal(results[["lm_data"]][["segmented"]], c(233), ignore_attr = TRUE, tolerance = 0.2)
+results[["binomial_data"]][["fastcpd"]] <-
+ fastcpd::fastcpd.binomial(binomial_data, r.progress = FALSE)@cp_set
+results[["binomial_data"]][["fastcpd"]]
+testthat::expect_equal(results[["binomial_data"]][["fastcpd"]], 302, tolerance = 0.2)
+results[["binomial_data"]][["strucchange"]] <-
+ strucchange::breakpoints(y ~ . - 1, data = binomial_data)$breakpoints
+results[["binomial_data"]][["strucchange"]]
+testthat::expect_equal(results[["binomial_data"]][["strucchange"]], 297, tolerance = 0.2)
+results[["poisson_data"]][["fastcpd"]] <-
+ fastcpd::fastcpd.poisson(poisson_data, r.progress = FALSE)@cp_set
+results[["poisson_data"]][["fastcpd"]]
+testthat::expect_equal(results[["poisson_data"]][["fastcpd"]], c(498, 805, 1003), tolerance = 0.2)
+results[["poisson_data"]][["strucchange"]] <-
+ strucchange::breakpoints(y ~ . - 1, data = poisson_data)$breakpoints
+results[["poisson_data"]][["strucchange"]]
+testthat::expect_equal(results[["poisson_data"]][["strucchange"]], 935, tolerance = 0.2)
+results[["lasso_data"]][["fastcpd"]] <-
+ fastcpd::fastcpd.lasso(lasso_data, r.progress = FALSE)@cp_set
+results[["lasso_data"]][["fastcpd"]]
+testthat::expect_equal(results[["lasso_data"]][["fastcpd"]], c(79, 199, 320), tolerance = 0.2)
+results[["lasso_data"]][["strucchange"]] <-
+ strucchange::breakpoints(y ~ . - 1, data = lasso_data)$breakpoints
+results[["lasso_data"]][["strucchange"]]
+testthat::expect_equal(results[["lasso_data"]][["strucchange"]], c(80, 200, 321), tolerance = 0.2)
+results[["ar_data"]][["fastcpd"]] <-
+ fastcpd::fastcpd.ar(ar_data, 3, r.progress = FALSE)@cp_set
+results[["ar_data"]][["fastcpd"]]
+testthat::expect_equal(results[["ar_data"]][["fastcpd"]], c(614), tolerance = 0.2)
+results[["ar_data"]][["CptNonPar"]] <-
+ CptNonPar::np.mojo(ar_data, G = floor(length(ar_data) / 6))$cpts
+results[["ar_data"]][["CptNonPar"]]
+testthat::expect_equal(results[["ar_data"]][["CptNonPar"]], numeric(0), tolerance = 0.2)
+results[["ar_data"]][["segmented"]] <-
+ segmented::segmented(
+ lm(
+ y ~ x + 1, data.frame(y = ar_data, x = seq_along(ar_data))
+ ),
+ seg.Z = ~ x
+ )$psi[, "Est."]
+results[["ar_data"]][["segmented"]]
+testthat::expect_equal(results[["ar_data"]][["segmented"]], c(690), ignore_attr = TRUE, tolerance = 0.2)
+results[["ar_data"]][["mcp"]] <-
+ mcp::mcp(
+ list(y ~ 1 + ar(3), ~ 0 + ar(3)),
+ data = data.frame(y = ar_data, x = seq_along(ar_data)),
+ par_x = "x"
+ )
+if (requireNamespace("mcp", quietly = TRUE)) {
+ plot(results[["ar_data"]][["mcp"]])
+}
+results[["garch_data"]][["fastcpd"]] <-
+ fastcpd::fastcpd.garch(garch_data, c(1, 1), r.progress = FALSE)@cp_set
+results[["garch_data"]][["fastcpd"]]
+testthat::expect_equal(results[["garch_data"]][["fastcpd"]], c(205), tolerance = 0.2)
+results[["garch_data"]][["CptNonPar"]] <-
+ CptNonPar::np.mojo(garch_data, G = floor(length(garch_data) / 6))$cpts
+results[["garch_data"]][["CptNonPar"]]
+testthat::expect_equal(results[["garch_data"]][["CptNonPar"]], c(206), tolerance = 0.2)
+results[["garch_data"]][["strucchange"]] <-
+ strucchange::breakpoints(x ~ 1, data = data.frame(x = garch_data))$breakpoints
+results[["garch_data"]][["strucchange"]]
+testthat::expect_equal(results[["garch_data"]][["strucchange"]], NA, tolerance = 0.2)
+results[["var_data"]][["fastcpd"]] <-
+ fastcpd::fastcpd.var(var_data, 2, r.progress = FALSE)@cp_set
+results[["var_data"]][["fastcpd"]]
+testthat::expect_equal(results[["var_data"]][["fastcpd"]], c(500), tolerance = 0.2)
+results[["var_data"]][["VARDetect"]] <- VARDetect::tbss(var_data)$cp
+results[["var_data"]][["VARDetect"]]
+testthat::expect_equal(results[["var_data"]][["VARDetect"]], c(501), tolerance = 0.2)
+well_log <- fastcpd::well_log
+well_log <- well_log[well_log > 1e5]
+
+results[["well_log"]] <- list(
+ fastcpd = fastcpd.mean(well_log, trim = 0.003)@cp_set,
+ changepoint = changepoint::cpt.mean(well_log)@cpts,
+ CptNonPar =
+ CptNonPar::np.mojo(well_log, G = floor(length(well_log) / 6))$cpts,
+ strucchange = strucchange::breakpoints(
y ~ 1, data = data.frame(y = well_log)
)$breakpoints,
ecp = ecp::e.divisive(matrix(well_log))$estimates,
@@ -884,15 +2126,15 @@ result <- list(
)$trend$cp,
stepR = stepR::stepFit(well_log, alpha = 0.5)$rightEnd
)
-
-package_list <- sort(names(result), decreasing = TRUE)
+results[["well_log"]]
+package_list <- sort(names(results[["well_log"]]), decreasing = TRUE)
comparison_table <- NULL
for (package_index in seq_along(package_list)) {
package <- package_list[[package_index]]
comparison_table <- rbind(
comparison_table,
data.frame(
- change_point = result[[package]],
+ change_point = results[["well_log"]][[package]],
package = package,
y_offset = (package_index - 1) * 1000
)
@@ -908,54 +2150,44 @@ for (i in seq_len(length(most_selected) - 1)) {
}
}
(most_selected <- most_selected[!is.na(most_selected)])
-#> [1] 6.0 314.0 434.0 704.0 776.0 1021.0 1057.0 1347.0 1405.0 1502.0 1661.0 1842.0 2023.0 2202.0
-#> [15] 2384.5 2445.0 2507.0 2567.5 2738.0 2921.0 3094.0 3251.0 3464.0 3499.0 3622.0 3709.0 3820.0 3976.0
-```
-
-```{r detection-comparison-well-log-plot, eval = FALSE}
-ggplot2::ggplot() +
- ggplot2::geom_point(
- data = data.frame(x = seq_along(well_log), y = c(well_log)),
- ggplot2::aes(x = x, y = y)
- ) +
- ggplot2::geom_vline(
- xintercept = most_selected,
- color = "black",
- linetype = "dashed",
- alpha = 0.2
- ) +
- ggplot2::geom_point(
- data = comparison_table,
- ggplot2::aes(x = change_point, y = 50000 + y_offset, color = package),
- shape = 17,
- size = 1.9
- ) +
- ggplot2::geom_hline(
- data = comparison_table,
- ggplot2::aes(yintercept = 50000 + y_offset, color = package),
- linetype = "dashed",
- alpha = 0.1
- ) +
- ggplot2::coord_cartesian(
- ylim = c(50000 - 500, max(well_log) + 1000),
- xlim = c(-200, length(well_log) + 200),
- expand = FALSE
- ) +
- ggplot2::theme(
- panel.background = ggplot2::element_blank(),
- panel.border = ggplot2::element_rect(colour = "black", fill = NA),
- panel.grid.major = ggplot2::element_blank(),
- panel.grid.minor = ggplot2::element_blank()
- ) +
- ggplot2::xlab(NULL) + ggplot2::ylab(NULL)
-```
-
-
-
-# Time comparison using `well_log`
-
-```{r time-comparison-well-log, eval = FALSE}
-microbenchmark_result <- microbenchmark::microbenchmark(
+if (requireNamespace("ggplot2", quietly = TRUE)) {
+ ggplot2::ggplot() +
+ ggplot2::geom_point(
+ data = data.frame(x = seq_along(well_log), y = c(well_log)),
+ ggplot2::aes(x = x, y = y)
+ ) +
+ ggplot2::geom_vline(
+ xintercept = most_selected,
+ color = "black",
+ linetype = "dashed",
+ alpha = 0.2
+ ) +
+ ggplot2::geom_point(
+ data = comparison_table,
+ ggplot2::aes(x = change_point, y = 50000 + y_offset, color = package),
+ shape = 17,
+ size = 1.9
+ ) +
+ ggplot2::geom_hline(
+ data = comparison_table,
+ ggplot2::aes(yintercept = 50000 + y_offset, color = package),
+ linetype = "dashed",
+ alpha = 0.1
+ ) +
+ ggplot2::coord_cartesian(
+ ylim = c(50000 - 500, max(well_log) + 1000),
+ xlim = c(-200, length(well_log) + 200),
+ expand = FALSE
+ ) +
+ ggplot2::theme(
+ panel.background = ggplot2::element_blank(),
+ panel.border = ggplot2::element_rect(colour = "black", fill = NA),
+ panel.grid.major = ggplot2::element_blank(),
+ panel.grid.minor = ggplot2::element_blank()
+ ) +
+ ggplot2::xlab(NULL) + ggplot2::ylab(NULL)
+}
+results[["microbenchmark"]] <- microbenchmark::microbenchmark(
fastcpd = fastcpd::fastcpd.mean(well_log, trim = 0.003, r.progress = FALSE),
changepoint = changepoint::cpt.mean(well_log, method = "PELT"),
CptNonPar = CptNonPar::np.mojo(well_log, G = floor(length(well_log) / 6)),
@@ -986,46 +2218,11 @@ microbenchmark_result <- microbenchmark::microbenchmark(
not = not::not(well_log, contrast = "pcwsConstMean"),
times = 10
)
-#> Unit: milliseconds
-#> expr min lq mean median uq max neval
-#> fastcpd 8.411284e+01 8.692226e+01 1.011440e+02 1.044509e+02 1.089672e+02 118.05842 10
-#> changepoint 3.206773e+01 3.377081e+01 6.843465e+01 3.857181e+01 5.243834e+01 244.76672 10
-#> CptNonPar 1.827381e+04 1.901094e+04 2.002752e+04 1.985180e+04 2.076803e+04 22511.59316 10
-#> strucchange 5.955079e+04 6.059315e+04 6.185727e+04 6.131291e+04 6.312073e+04 64638.93090 10
-#> ecp 7.590543e+05 7.707573e+05 7.859080e+05 7.830752e+05 8.093015e+05 810339.52140 10
-#> breakfast 9.170992e+03 9.344041e+03 9.628236e+03 9.382078e+03 9.628663e+03 11073.79318 10
-#> wbs 1.139078e+02 1.145472e+02 1.178167e+02 1.166746e+02 1.201676e+02 127.27064 10
-#> mosum 1.172847e+00 1.231747e+00 1.740727e+01 1.416854e+00 1.919586e+00 160.76997 10
-#> fpop 2.588228e+00 2.630407e+00 4.587742e+00 2.832556e+00 3.312986e+00 18.52067 10
-#> gfpop 5.971245e+01 6.072684e+01 6.533492e+01 6.172578e+01 6.839653e+01 87.89407 10
-#> InspectChangepoint 1.698673e+02 1.909034e+02 2.392539e+02 2.117010e+02 3.004474e+02 329.87724 10
-#> jointseg 2.000894e+01 2.136878e+01 2.551210e+01 2.167757e+01 2.403593e+01 47.98397 10
-#> Rbeast 6.533998e+02 6.625203e+02 6.783586e+02 6.792646e+02 6.875840e+02 723.45376 10
-#> stepR 2.793709e+01 2.902084e+01 4.380857e+01 3.068416e+01 3.227125e+01 164.81082 10
-#> not 9.599763e+01 9.701856e+01 1.028601e+02 1.012292e+02 1.049974e+02 120.73529 10
-```
-
-```{r time-comparison-well-log-plot, eval = FALSE}
-ggplot2::autoplot(microbenchmark_result)
-```
-
-
-
-# Related issues
-
-- `mosum`: [[#4](https://github.com/doccstat/fastcpd/issues/4)].
-- `mcp`: [link 1](https://github.com/doccstat/fastcpd/actions/runs/7272426093/job/19814531091)
- and [link 2](https://github.com/doccstat/fastcpd/actions/runs/7272426093/job/19814531284).
-- `bcp`: [[#5](https://github.com/doccstat/fastcpd/issues/5)].
-- `gfpop`: [[#10](https://github.com/doccstat/fastcpd/issues/10)].
-
-# Acknowledgements
-
-- [Dr. Vito Muggeo](https://www.unipa.it/persone/docenti/m/vito.muggeo),
- author of the `segmented` package for the tips about the piece-wise
- constant function.
-
-# Appendix: all code snippets
-
-```{r ref.label = knitr::all_labels(), echo = TRUE, eval = FALSE}
+results[["microbenchmark"]]
+if (requireNamespace("ggplot2", quietly = TRUE) && requireNamespace("microbenchmark", quietly = TRUE)) {
+ ggplot2::autoplot(results[["microbenchmark"]])
+}
+if (!file.exists("comparison-packages-results.RData")) {
+ save(results, file = "comparison-packages-results.RData")
+}
```
diff --git a/vignettes/comparison-packages.Rmd.original b/vignettes/comparison-packages.Rmd.original
new file mode 100644
index 00000000..be272d83
--- /dev/null
+++ b/vignettes/comparison-packages.Rmd.original
@@ -0,0 +1,1358 @@
+---
+title: "Comparison with other R packages"
+output: rmarkdown::html_vignette
+description: |
+ Speed and result comparison with other R packages.
+vignette: >
+ %\VignetteIndexEntry{Comparison with other R packages}
+ %\VignetteEngine{knitr::rmarkdown}
+ %\VignetteEncoding{UTF-8}
+---
+
+```{r, include = FALSE}
+knitr::opts_chunk$set(
+ collapse = TRUE, comment = "#>", eval = TRUE, cache = FALSE,
+ warning = FALSE, fig.width = 8, fig.height = 5
+)
+
+if (requireNamespace("microbenchmark", quietly = TRUE)) {
+ library(microbenchmark)
+}
+
+if (file.exists("comparison-packages-results.RData")) {
+ # Available at https://pcloud.xingchi.li/comparison-packages-results.RData
+ load("comparison-packages-results.RData")
+} else {
+ results <- list()
+}
+```
+
+# Data setup
+
+## Univariate mean change
+
+```{r data-setup-univariate-mean-change}
+# Univariate mean change
+set.seed(1)
+p <- 1
+mean_data_1 <- rbind(
+ mvtnorm::rmvnorm(300, mean = rep(0, p), sigma = diag(100, p)),
+ mvtnorm::rmvnorm(400, mean = rep(50, p), sigma = diag(100, p)),
+ mvtnorm::rmvnorm(300, mean = rep(2, p), sigma = diag(100, p))
+)
+
+plot.ts(mean_data_1)
+```
+
+## Univariate mean and/or variance change
+
+```{r data-setup-univariate-mean-and-or-variance-change}
+# Univariate mean and/or variance change
+set.seed(1)
+p <- 1
+mv_data_1 <- rbind(
+ mvtnorm::rmvnorm(300, mean = rep(0, p), sigma = diag(1, p)),
+ mvtnorm::rmvnorm(400, mean = rep(10, p), sigma = diag(1, p)),
+ mvtnorm::rmvnorm(300, mean = rep(0, p), sigma = diag(100, p)),
+ mvtnorm::rmvnorm(300, mean = rep(0, p), sigma = diag(1, p)),
+ mvtnorm::rmvnorm(400, mean = rep(10, p), sigma = diag(1, p)),
+ mvtnorm::rmvnorm(300, mean = rep(10, p), sigma = diag(100, p))
+)
+
+plot.ts(mv_data_1)
+```
+
+## Multivariate mean change
+
+```{r data-setup-multivariate-mean-change}
+# Multivariate mean change
+set.seed(1)
+p <- 3
+mean_data_3 <- rbind(
+ mvtnorm::rmvnorm(300, mean = rep(0, p), sigma = diag(100, p)),
+ mvtnorm::rmvnorm(400, mean = rep(50, p), sigma = diag(100, p)),
+ mvtnorm::rmvnorm(300, mean = rep(2, p), sigma = diag(100, p))
+)
+
+plot.ts(mean_data_3)
+```
+
+## Multivariate mean and/or variance change
+
+```{r data-setup-multivariate-mean-and-or-variance-change}
+# Multivariate mean and/or variance change
+set.seed(1)
+p <- 4
+mv_data_3 <- rbind(
+ mvtnorm::rmvnorm(300, mean = rep(0, p), sigma = diag(1, p)),
+ mvtnorm::rmvnorm(400, mean = rep(10, p), sigma = diag(1, p)),
+ mvtnorm::rmvnorm(300, mean = rep(0, p), sigma = diag(100, p)),
+ mvtnorm::rmvnorm(300, mean = rep(0, p), sigma = diag(1, p)),
+ mvtnorm::rmvnorm(400, mean = rep(10, p), sigma = diag(1, p)),
+ mvtnorm::rmvnorm(300, mean = rep(10, p), sigma = diag(100, p))
+)
+
+plot.ts(mv_data_3)
+```
+
+## Linear regression
+
+```{r data-setup-linear-regression}
+# Linear regression
+set.seed(1)
+n <- 300
+p <- 4
+x <- mvtnorm::rmvnorm(n, rep(0, p), diag(p))
+theta_0 <- rbind(c(1, 3.2, -1, 0), c(-1, -0.5, 2.5, -2), c(0.8, 0, 1, 2))
+y <- c(
+ x[1:100, ] %*% theta_0[1, ] + rnorm(100, 0, 3),
+ x[101:200, ] %*% theta_0[2, ] + rnorm(100, 0, 3),
+ x[201:n, ] %*% theta_0[3, ] + rnorm(100, 0, 3)
+)
+lm_data <- data.frame(y = y, x = x)
+
+plot.ts(lm_data)
+```
+
+## Logistic regression
+
+```{r data-setup-logistic-regression}
+# Logistic regression
+set.seed(1)
+n <- 500
+p <- 4
+x <- mvtnorm::rmvnorm(n, rep(0, p), diag(p))
+theta <- rbind(rnorm(p, 0, 1), rnorm(p, 2, 1))
+y <- c(
+ rbinom(300, 1, 1 / (1 + exp(-x[1:300, ] %*% theta[1, ]))),
+ rbinom(200, 1, 1 / (1 + exp(-x[301:n, ] %*% theta[2, ])))
+)
+binomial_data <- data.frame(y = y, x = x)
+
+plot.ts(binomial_data)
+```
+
+## Poisson regression
+
+```{r data-setup-poisson-regression}
+# Poisson regression
+set.seed(1)
+n <- 1100
+p <- 3
+x <- mvtnorm::rmvnorm(n, rep(0, p), diag(p))
+delta <- rnorm(p)
+theta_0 <- c(1, 0.3, -1)
+y <- c(
+ rpois(500, exp(x[1:500, ] %*% theta_0)),
+ rpois(300, exp(x[501:800, ] %*% (theta_0 + delta))),
+ rpois(200, exp(x[801:1000, ] %*% theta_0)),
+ rpois(100, exp(x[1001:1100, ] %*% (theta_0 - delta)))
+)
+poisson_data <- data.frame(y = y, x = x)
+
+plot.ts(log(poisson_data$y))
+plot.ts(poisson_data[, -1])
+```
+
+## Lasso
+
+```{r data-setup-lasso}
+# Lasso
+set.seed(1)
+n <- 480
+p_true <- 6
+p <- 50
+x <- mvtnorm::rmvnorm(n, rep(0, p), diag(p))
+theta_0 <- rbind(
+ runif(p_true, -5, -2),
+ runif(p_true, -3, 3),
+ runif(p_true, 2, 5),
+ runif(p_true, -5, 5)
+)
+theta_0 <- cbind(theta_0, matrix(0, ncol = p - p_true, nrow = 4))
+y <- c(
+ x[1:80, ] %*% theta_0[1, ] + rnorm(80, 0, 1),
+ x[81:200, ] %*% theta_0[2, ] + rnorm(120, 0, 1),
+ x[201:320, ] %*% theta_0[3, ] + rnorm(120, 0, 1),
+ x[321:n, ] %*% theta_0[4, ] + rnorm(160, 0, 1)
+)
+lasso_data <- data.frame(y = y, x = x)
+
+plot.ts(lasso_data[, seq_len(p_true + 1)])
+```
+
+## AR(3)
+
+```{r data-setup-ar3}
+# AR(3)
+set.seed(1)
+n <- 1000
+x <- rep(0, n + 3)
+for (i in 1:600) {
+ x[i + 3] <- 0.6 * x[i + 2] - 0.2 * x[i + 1] + 0.1 * x[i] + rnorm(1, 0, 3)
+}
+for (i in 601:1000) {
+ x[i + 3] <- 0.3 * x[i + 2] + 0.4 * x[i + 1] + 0.2 * x[i] + rnorm(1, 0, 3)
+}
+ar_data <- x[-seq_len(3)]
+
+plot.ts(ar_data)
+```
+
+## GARCH(1, 1)
+
+```{r data-setup-garch11}
+# GARCH(1, 1)
+set.seed(1)
+n <- 400
+sigma_2 <- rep(1, n + 1)
+x <- rep(0, n + 1)
+for (i in seq_len(200)) {
+ sigma_2[i + 1] <- 20 + 0.5 * x[i]^2 + 0.1 * sigma_2[i]
+ x[i + 1] <- rnorm(1, 0, sqrt(sigma_2[i + 1]))
+}
+for (i in 201:400) {
+ sigma_2[i + 1] <- 1 + 0.1 * x[i]^2 + 0.5 * sigma_2[i]
+ x[i + 1] <- rnorm(1, 0, sqrt(sigma_2[i + 1]))
+}
+garch_data <- x[-1]
+
+plot.ts(garch_data)
+```
+
+## VAR(2)
+
+```{r data-setup-var2}
+# VAR(2)
+set.seed(1)
+n <- 800
+p <- 2
+theta_1 <- matrix(c(-0.3, 0.6, -0.5, 0.4, 0.2, 0.2, 0.2, -0.2), nrow = p)
+theta_2 <- matrix(c(0.3, -0.4, 0.1, -0.5, -0.5, -0.2, -0.5, 0.2), nrow = p)
+x <- matrix(0, n + 2, p)
+for (i in 1:500) {
+ x[i + 2, ] <- theta_1 %*% c(x[i + 1, ], x[i, ]) + rnorm(p, 0, 1)
+}
+for (i in 501:n) {
+ x[i + 2, ] <- theta_2 %*% c(x[i + 1, ], x[i, ]) + rnorm(p, 0, 1)
+}
+var_data <- x[-seq_len(2), ]
+
+plot.ts(var_data)
+```
+
+# Univariate mean change
+
+The true change points are 300 and 700.
+Some methods are plotted due to the un-retrievable change points.
+
+```{r univariate-mean-change-fastcpd, eval = FALSE}
+results[["mean_data_1"]][["fastcpd"]] <-
+ fastcpd::fastcpd.mean(mean_data_1, r.progress = FALSE)@cp_set
+```
+
+```{r univariate-mean-change-fastcpd-result}
+results[["mean_data_1"]][["fastcpd"]]
+```
+
+```{r univariate-mean-change-fastcpd-testthat, include = FALSE}
+testthat::expect_equal(results[["mean_data_1"]][["fastcpd"]], c(300, 700), tolerance = 0.2)
+```
+
+```{r univariate-mean-change-CptNonPar, eval = FALSE}
+results[["mean_data_1"]][["CptNonPar"]] <-
+ CptNonPar::np.mojo(mean_data_1, G = floor(length(mean_data_1) / 6))$cpts
+```
+
+```{r univariate-mean-change-CptNonPar-result}
+results[["mean_data_1"]][["CptNonPar"]]
+```
+
+```{r univariate-mean-change-CptNonPar-testthat, include = FALSE}
+testthat::expect_equal(results[["mean_data_1"]][["CptNonPar"]], c(300, 700), tolerance = 0.2)
+```
+
+```{r univariate-mean-change-strucchange, eval = FALSE}
+results[["mean_data_1"]][["strucchange"]] <-
+ strucchange::breakpoints(y ~ 1, data = data.frame(y = mean_data_1))$breakpoints
+```
+
+```{r univariate-mean-change-strucchange-result}
+results[["mean_data_1"]][["strucchange"]]
+```
+
+```{r univariate-mean-change-strucchange-testthat, include = FALSE}
+testthat::expect_equal(results[["mean_data_1"]][["strucchange"]], c(300, 700), tolerance = 0.2)
+```
+
+```{r univariate-mean-change-ecp, eval = FALSE}
+results[["mean_data_1"]][["ecp"]] <- ecp::e.divisive(mean_data_1)$estimates
+```
+
+```{r univariate-mean-change-ecp-result}
+results[["mean_data_1"]][["ecp"]]
+```
+
+```{r univariate-mean-change-ecp-testthat, include = FALSE}
+testthat::expect_equal(results[["mean_data_1"]][["ecp"]], c(1, 301, 701, 1001), tolerance = 0.2)
+```
+
+```{r univariate-mean-change-changepoint, eval = FALSE}
+results[["mean_data_1"]][["changepoint"]] <-
+ changepoint::cpt.mean(c(mean_data_1))@cpts
+```
+
+```{r univariate-mean-change-changepoint-result}
+results[["mean_data_1"]][["changepoint"]]
+```
+
+```{r univariate-mean-change-changepoint-testthat, include = FALSE}
+testthat::expect_equal(results[["mean_data_1"]][["changepoint"]], c(300, 1000), tolerance = 0.2)
+```
+
+```{r univariate-mean-change-breakfast, eval = FALSE}
+results[["mean_data_1"]][["breakfast"]] <-
+ breakfast::breakfast(mean_data_1)$cptmodel.list[[6]]$cpts
+```
+
+```{r univariate-mean-change-breakfast-result}
+results[["mean_data_1"]][["breakfast"]]
+```
+
+```{r univariate-mean-change-breakfast-testthat, include = FALSE}
+testthat::expect_equal(results[["mean_data_1"]][["breakfast"]], c(300, 700), tolerance = 0.2)
+```
+
+```{r univariate-mean-change-wbs, eval = FALSE}
+results[["mean_data_1"]][["wbs"]] <-
+ wbs::wbs(mean_data_1)$cpt$cpt.ic$mbic.penalty
+```
+
+```{r univariate-mean-change-wbs-result}
+results[["mean_data_1"]][["wbs"]]
+```
+
+```{r univariate-mean-change-wbs-testthat, include = FALSE}
+testthat::expect_equal(results[["mean_data_1"]][["wbs"]], c(300, 700), tolerance = 0.2)
+```
+
+```{r univariate-mean-change-mosum, eval = FALSE}
+results[["mean_data_1"]][["mosum"]] <-
+ mosum::mosum(c(mean_data_1), G = 40)$cpts.info$cpts
+```
+
+```{r univariate-mean-change-mosum-result}
+results[["mean_data_1"]][["mosum"]]
+```
+
+```{r univariate-mean-change-mosum-testthat, include = FALSE}
+testthat::expect_equal(results[["mean_data_1"]][["mosum"]], c(300, 700), tolerance = 0.2)
+```
+
+```{r univariate-mean-change-fpop, eval = FALSE}
+results[["mean_data_1"]][["fpop"]] <-
+ fpop::Fpop(mean_data_1, nrow(mean_data_1))$t.est
+```
+
+```{r univariate-mean-change-fpop-result}
+results[["mean_data_1"]][["fpop"]]
+```
+
+```{r univariate-mean-change-fpop-testthat, include = FALSE}
+testthat::expect_equal(results[["mean_data_1"]][["fpop"]], c(300, 700, 1000), tolerance = 0.2)
+```
+
+```{r univariate-mean-change-gfpop, eval = FALSE}
+results[["mean_data_1"]][["gfpop"]] <-
+ gfpop::gfpop(
+ data = mean_data_1,
+ mygraph = gfpop::graph(
+ penalty = 2 * log(nrow(mean_data_1)) * gfpop::sdDiff(mean_data_1) ^ 2,
+ type = "updown"
+ ),
+ type = "mean"
+ )$changepoints
+```
+
+```{r univariate-mean-change-gfpop-result}
+results[["mean_data_1"]][["gfpop"]]
+```
+
+```{r univariate-mean-change-gfpop-testthat, include = FALSE}
+testthat::expect_equal(results[["mean_data_1"]][["gfpop"]], c(300, 700, 1000), tolerance = 0.2)
+```
+
+```{r univariate-mean-change-InspectChangepoint, eval = FALSE}
+results[["mean_data_1"]][["InspectChangepoint"]] <- InspectChangepoint::inspect(
+ t(mean_data_1),
+ threshold = InspectChangepoint::compute.threshold(
+ nrow(mean_data_1), ncol(mean_data_1)
+ )
+)$changepoints[, "location"]
+```
+
+```{r univariate-mean-change-InspectChangepoint-result}
+results[["mean_data_1"]][["InspectChangepoint"]]
+```
+
+```{r univariate-mean-change-InspectChangepoint-testthat, include = FALSE}
+testthat::expect_equal(results[["mean_data_1"]][["InspectChangepoint"]], c(300, 700), tolerance = 0.2)
+```
+
+```{r univariate-mean-change-jointseg, eval = FALSE}
+results[["mean_data_1"]][["jointseg"]] <-
+ jointseg::jointSeg(mean_data_1, K = 2)$bestBkp
+```
+
+```{r univariate-mean-change-jointseg-result}
+results[["mean_data_1"]][["jointseg"]]
+```
+
+```{r univariate-mean-change-jointseg-testthat, include = FALSE}
+testthat::expect_equal(results[["mean_data_1"]][["jointseg"]], c(300, 700), tolerance = 0.2)
+```
+
+```{r univariate-mean-change-Rbeast, eval = FALSE}
+results[["mean_data_1"]][["Rbeast"]] <-
+ Rbeast::beast(
+ mean_data_1, season = "none", print.progress = FALSE, quiet = TRUE
+ )$trend$cp
+```
+
+```{r univariate-mean-change-Rbeast-result}
+results[["mean_data_1"]][["Rbeast"]]
+```
+
+```{r univariate-mean-change-Rbeast-testthat, include = FALSE}
+testthat::expect_equal(results[["mean_data_1"]][["Rbeast"]], c(701, 301, NaN, NaN, NaN, NaN, NaN, NaN, NaN, NaN), tolerance = 0.2)
+```
+
+```{r univariate-mean-change-stepR, eval = FALSE}
+results[["mean_data_1"]][["stepR"]] <-
+ stepR::stepFit(mean_data_1, alpha = 0.5)$rightEnd
+```
+
+```{r univariate-mean-change-stepR-result}
+results[["mean_data_1"]][["stepR"]]
+```
+
+```{r univariate-mean-change-stepR-testthat, include = FALSE}
+testthat::expect_equal(results[["mean_data_1"]][["stepR"]], c(300, 700, 1000), tolerance = 0.2)
+```
+
+```{r univariate-mean-change-cpm, eval = FALSE}
+results[["mean_data_1"]][["cpm"]] <-
+ cpm::processStream(mean_data_1, cpmType = "Student")$changePoints
+```
+
+```{r univariate-mean-change-cpm-result}
+results[["mean_data_1"]][["cpm"]]
+```
+
+```{r univariate-mean-change-cpm-testthat, include = FALSE}
+testthat::expect_equal(results[["mean_data_1"]][["cpm"]], c(299, 699), tolerance = 0.2)
+```
+
+```{r univariate-mean-change-segmented, eval = FALSE}
+results[["mean_data_1"]][["segmented"]] <-
+ segmented::stepmented(
+ as.numeric(mean_data_1), npsi = 2
+ )$psi[, "Est."]
+```
+
+```{r univariate-mean-change-segmented-result}
+results[["mean_data_1"]][["segmented"]]
+```
+
+```{r univariate-mean-change-segmented-testthat, include = FALSE}
+testthat::expect_equal(results[["mean_data_1"]][["segmented"]], c(298, 699), ignore_attr = TRUE, tolerance = 0.2)
+```
+
+```{r univariate-mean-change-mcp, eval = FALSE}
+results[["mean_data_1"]][["mcp"]] <- mcp::mcp(
+ list(y ~ 1, ~ 1, ~ 1),
+ data = data.frame(y = mean_data_1, x = seq_len(nrow(mean_data_1))),
+ par_x = "x"
+)
+```
+
+```{r univariate-mean-change-mcp-result}
+if (requireNamespace("mcp", quietly = TRUE)) {
+ plot(results[["mean_data_1"]][["mcp"]])
+}
+```
+
+```{r univariate-mean-change-not, eval = FALSE}
+results[["mean_data_1"]][["not"]] <-
+ not::not(mean_data_1, contrast = "pcwsConstMean")
+```
+
+```{r univariate-mean-change-not-result}
+if (requireNamespace("not", quietly = TRUE)) {
+ plot(results[["mean_data_1"]][["not"]])
+}
+```
+
+```{r univariate-mean-change-bcp, eval = FALSE}
+results[["mean_data_1"]][["bcp"]] <- bcp::bcp(mean_data_1)
+```
+
+```{r univariate-mean-change-bcp-result}
+if (requireNamespace("bcp", quietly = TRUE)) {
+ plot(results[["mean_data_1"]][["bcp"]])
+}
+```
+
+# Univariate mean and/or variance change
+
+The true change points are 300, 700, 1000, 1300 and 1700.
+Some methods are plotted due to the un-retrievable change points.
+
+```{r univariate-mean-and-or-variance-change-fastcpd, eval = FALSE}
+results[["mv_data_1"]][["fastcpd"]] <-
+ fastcpd::fastcpd.mv(mv_data_1, r.progress = FALSE)@cp_set
+```
+
+```{r univariate-mean-and-or-variance-change-fastcpd-result}
+results[["mv_data_1"]][["fastcpd"]]
+```
+
+```{r univariate-mean-and-or-variance-change-fastcpd-testthat, include = FALSE}
+testthat::expect_equal(results[["mv_data_1"]][["fastcpd"]], c(300, 700, 1001, 1300, 1700), tolerance = 0.2)
+```
+
+```{r univariate-mean-and-or-variance-change-ecp, eval = FALSE}
+results[["mv_data_1"]][["ecp"]] <- ecp::e.divisive(mv_data_1)$estimates
+```
+
+```{r univariate-mean-and-or-variance-change-ecp-result}
+results[["mv_data_1"]][["ecp"]]
+```
+
+```{r univariate-mean-and-or-variance-change-ecp-testthat, include = FALSE}
+testthat::expect_equal(results[["mv_data_1"]][["ecp"]], c(1, 301, 701, 1001, 1301, 1701, 2001), tolerance = 0.2)
+```
+
+```{r univariate-mean-and-or-variance-change-changepoint, eval = FALSE}
+results[["mv_data_1"]][["changepoint"]] <-
+ changepoint::cpt.meanvar(c(mv_data_1))@cpts
+```
+
+```{r univariate-mean-and-or-variance-change-changepoint-result}
+results[["mv_data_1"]][["changepoint"]]
+```
+
+```{r univariate-mean-and-or-variance-change-changepoint-testthat, include = FALSE}
+testthat::expect_equal(results[["mv_data_1"]][["changepoint"]], c(300, 2000), tolerance = 0.2)
+```
+
+```{r univariate-mean-and-or-variance-change-CptNonPar, eval = FALSE}
+results[["mv_data_1"]][["CptNonPar"]] <-
+ CptNonPar::np.mojo(mv_data_1, G = floor(length(mv_data_1) / 6))$cpts
+```
+
+```{r univariate-mean-and-or-variance-change-CptNonPar-result}
+results[["mv_data_1"]][["CptNonPar"]]
+```
+
+```{r univariate-mean-and-or-variance-change-CptNonPar-testthat, include = FALSE}
+testthat::expect_equal(results[["mv_data_1"]][["CptNonPar"]], c(333, 700, 1300), tolerance = 0.2)
+```
+
+```{r univariate-mean-and-or-variance-change-cpm, eval = FALSE}
+results[["mv_data_1"]][["cpm"]] <-
+ cpm::processStream(mv_data_1, cpmType = "GLR")$changePoints
+```
+
+```{r univariate-mean-and-or-variance-change-cpm-result}
+results[["mv_data_1"]][["cpm"]]
+```
+
+```{r univariate-mean-and-or-variance-change-cpm-testthat, include = FALSE}
+testthat::expect_equal(results[["mv_data_1"]][["cpm"]], c(293, 300, 403, 408, 618, 621, 696, 1000, 1021, 1024, 1293, 1300, 1417, 1693, 1700, 1981), tolerance = 0.2)
+```
+
+```{r univariate-mean-and-or-variance-change-InspectChangepoint, eval = FALSE}
+results[["mv_data_1"]][["InspectChangepoint"]] <- InspectChangepoint::inspect(
+ t(mv_data_1),
+ threshold = InspectChangepoint::compute.threshold(
+ nrow(mv_data_1), ncol(mv_data_1)
+ )
+)$changepoints[, "location"]
+```
+
+```{r univariate-mean-and-or-variance-change-InspectChangepoint-result}
+results[["mv_data_1"]][["InspectChangepoint"]]
+```
+
+```{r univariate-mean-and-or-variance-change-InspectChangepoint-testthat, include = FALSE}
+testthat::expect_equal(results[["mv_data_1"]][["InspectChangepoint"]], c(
+ 300, 700, 701, 702, 704, 707, 708, 712, 715, 716, 717, 718,
+ 721, 722, 723, 726, 727, 729, 731, 732, 734, 736, 740, 742,
+ 744, 746, 748, 750, 753, 755, 756, 757, 759, 760, 762, 764,
+ 765, 766, 768, 769, 771, 772, 774, 776, 777, 784, 785, 786,
+ 789, 791, 792, 794, 797, 798, 799, 801, 802, 803, 807, 809,
+ 810, 813, 815, 817, 819, 826, 827, 828, 829, 831, 833, 835,
+ 836, 837, 838, 840, 841, 842, 843, 845, 848, 849, 852, 854,
+ 860, 862, 864, 866, 868, 870, 872, 875, 879, 881, 884, 886,
+ 887, 888, 889, 896, 897, 898, 899, 901, 903, 904, 905, 906,
+ 909, 910, 912, 913, 915, 917, 919, 921, 922, 923, 925, 927,
+ 928, 932, 934, 936, 937, 940, 944, 945, 947, 948, 949, 951,
+ 956, 958, 959, 961, 962, 963, 964, 966, 967, 968, 972, 974,
+ 976, 978, 979, 986, 988, 990, 992, 995, 998, 1000, 1300, 1700,
+ 1702, 1703, 1704, 1705, 1708, 1710, 1712, 1714, 1716, 1717, 1718, 1720,
+ 1721, 1723, 1725, 1726, 1727, 1729, 1731, 1733, 1735, 1736, 1737, 1739,
+ 1742, 1745, 1747, 1748, 1752, 1754, 1756, 1758, 1759, 1760, 1766, 1768,
+ 1770, 1771, 1773, 1775, 1778, 1782, 1784, 1785, 1790, 1792, 1793, 1795,
+ 1796, 1797, 1799, 1800, 1802, 1803, 1804, 1805, 1806, 1807, 1808, 1809,
+ 1813, 1815, 1816, 1818, 1821, 1824, 1825, 1827, 1828, 1829, 1833, 1835,
+ 1837, 1840, 1841, 1842, 1848, 1849, 1851, 1852, 1854, 1855, 1857, 1859,
+ 1860, 1862, 1863, 1865, 1867, 1868, 1876, 1878, 1879, 1880, 1882, 1883,
+ 1884, 1886, 1887, 1889, 1894, 1898, 1899, 1905, 1906, 1907, 1908, 1909,
+ 1912, 1919, 1920, 1921, 1923, 1924, 1926, 1927, 1928, 1930, 1933, 1934,
+ 1935, 1936, 1938, 1940, 1941, 1944, 1947, 1950, 1952, 1954, 1955, 1956,
+ 1960, 1962, 1963, 1965, 1966, 1967, 1969, 1970, 1974, 1976, 1977, 1978,
+ 1980, 1985, 1987, 1988, 1990, 1996, 1997, 1998
+), tolerance = 0.2)
+```
+
+```{r univariate-mean-and-or-variance-change-Rbeast, eval = FALSE}
+results[["mv_data_1"]][["Rbeast"]] <-
+ Rbeast::beast(
+ mv_data_1, season = "none", print.progress = FALSE, quiet = TRUE
+ )$trend$cp
+```
+
+```{r univariate-mean-and-or-variance-change-Rbeast-result}
+results[["mv_data_1"]][["Rbeast"]]
+```
+
+```{r univariate-mean-and-or-variance-change-Rbeast-testthat, include = FALSE}
+testthat::expect_equal(results[["mv_data_1"]][["Rbeast"]], c(1855, 1794, 1301, 1986, 301, 703, 1981, 1769, 1860, 709), tolerance = 0.2)
+```
+
+```{r univariate-mean-and-or-variance-change-mcp, eval = FALSE}
+results[["mv_data_1"]][["mcp"]] <- mcp::mcp(
+ list(y ~ 1, ~ 1, ~ 1, ~ 1, ~ 1, ~ 1),
+ data = data.frame(y = mv_data_1, x = seq_len(nrow(mv_data_1))),
+ par_x = "x"
+)
+```
+
+```{r univariate-mean-and-or-variance-change-mcp-result}
+if (requireNamespace("mcp", quietly = TRUE)) {
+ plot(results[["mv_data_1"]][["mcp"]])
+}
+```
+
+```{r univariate-mean-and-or-variance-change-not, eval = FALSE}
+results[["mv_data_1"]][["not"]] <-
+ not::not(mv_data_1, contrast = "pcwsConstMeanVar")
+```
+
+```{r univariate-mean-and-or-variance-change-not-result}
+if (requireNamespace("not", quietly = TRUE)) {
+ plot(results[["mv_data_1"]][["not"]])
+}
+```
+
+# Multivariate mean change
+
+The true change points are 300 and 700.
+Some methods are plotted due to the un-retrievable change points.
+
+```{r multivariate-mean-change-fastcpd, eval = FALSE}
+results[["mean_data_3"]][["fastcpd"]] <-
+ fastcpd::fastcpd.mean(mean_data_3, r.progress = FALSE)@cp_set
+```
+
+```{r multivariate-mean-change-fastcpd-result}
+results[["mean_data_3"]][["fastcpd"]]
+```
+
+```{r multivariate-mean-change-fastcpd-testthat, include = FALSE}
+testthat::expect_equal(results[["mean_data_3"]][["fastcpd"]], c(300, 700), tolerance = 0.2)
+```
+
+```{r multivariate-mean-change-CptNonPar, eval = FALSE}
+results[["mean_data_3"]][["CptNonPar"]] <-
+ CptNonPar::np.mojo(mean_data_3, G = floor(nrow(mean_data_3) / 6))$cpts
+```
+
+```{r multivariate-mean-change-CptNonPar-result}
+results[["mean_data_3"]][["CptNonPar"]]
+```
+
+```{r multivariate-mean-change-CptNonPar-testthat, include = FALSE}
+testthat::expect_equal(results[["mean_data_3"]][["CptNonPar"]], c(300, 700), tolerance = 0.2)
+```
+
+```{r multivariate-mean-change-InspectChangepoint, eval = FALSE}
+results[["mean_data_3"]][["InspectChangepoint"]] <- InspectChangepoint::inspect(
+ t(mean_data_3),
+ threshold = InspectChangepoint::compute.threshold(
+ nrow(mean_data_3), ncol(mean_data_3)
+ )
+)$changepoints[, "location"]
+```
+
+```{r multivariate-mean-change-InspectChangepoint-result}
+results[["mean_data_3"]][["InspectChangepoint"]]
+```
+
+```{r multivariate-mean-change-InspectChangepoint-testthat, include = FALSE}
+testthat::expect_equal(results[["mean_data_3"]][["InspectChangepoint"]], c(300, 700), tolerance = 0.2)
+```
+
+```{r multivariate-mean-change-jointseg, eval = FALSE}
+results[["mean_data_3"]][["jointseg"]] <-
+ jointseg::jointSeg(mean_data_3, K = 2)$bestBkp
+```
+
+```{r multivariate-mean-change-jointseg-result}
+results[["mean_data_3"]][["jointseg"]]
+```
+
+```{r multivariate-mean-change-jointseg-testthat, include = FALSE}
+testthat::expect_equal(results[["mean_data_3"]][["jointseg"]], c(300, 700), tolerance = 0.2)
+```
+
+```{r multivariate-mean-change-Rbeast, eval = FALSE}
+results[["mean_data_3"]][["Rbeast"]] <-
+ Rbeast::beast123(
+ mean_data_3,
+ metadata = list(whichDimIsTime = 1),
+ season = "none"
+ )$trend$cp
+```
+
+```{r multivariate-mean-change-Rbeast-result}
+results[["mean_data_3"]][["Rbeast"]]
+```
+
+```{r multivariate-mean-change-Rbeast-testthat, include = FALSE}
+testthat::expect_equal(results[["mean_data_3"]][["Rbeast"]], matrix(c(
+ 301, 701, 301,
+ 701, 301, 701,
+ 142, 117, 926,
+ 305, NaN, NaN,
+ 705, NaN, NaN,
+ 694, NaN, NaN,
+ NaN, NaN, NaN,
+ NaN, NaN, NaN,
+ NaN, NaN, NaN,
+ NaN, NaN, NaN
+), nrow = 10, ncol = 3, byrow = TRUE), tolerance = 0.2)
+```
+
+```{r multivariate-mean-change-strucchange, eval = FALSE}
+results[["mean_data_3"]][["strucchange"]] <-
+ strucchange::breakpoints(
+ cbind(y.1, y.2, y.3) ~ 1, data = data.frame(y = mean_data_3)
+ )$breakpoints
+```
+
+```{r multivariate-mean-change-strucchange-result}
+results[["mean_data_3"]][["strucchange"]]
+```
+
+```{r multivariate-mean-change-strucchange-testthat, include = FALSE}
+testthat::expect_equal(results[["mean_data_3"]][["strucchange"]], c(300, 700), tolerance = 0.2)
+```
+
+```{r multivariate-mean-change-ecp, eval = FALSE}
+results[["mean_data_3"]][["ecp"]] <- ecp::e.divisive(mean_data_3)$estimates
+```
+
+```{r multivariate-mean-change-ecp-result}
+results[["mean_data_3"]][["ecp"]]
+```
+
+```{r multivariate-mean-change-ecp-testthat, include = FALSE}
+testthat::expect_equal(results[["mean_data_3"]][["ecp"]], c(1, 301, 701, 1001), tolerance = 0.2)
+```
+
+```{r multivariate-mean-change-bcp, eval = FALSE}
+results[["mean_data_3"]][["bcp"]] <- bcp::bcp(mean_data_3)
+```
+
+```{r multivariate-mean-change-bcp-result}
+if (requireNamespace("bcp", quietly = TRUE)) {
+ plot(results[["mean_data_3"]][["bcp"]])
+}
+```
+
+# Multivariate mean and/or variance change
+
+The true change points are 300, 700, 1000, 1300 and 1700.
+Some methods are plotted due to the un-retrievable change points.
+
+```{r multivariate-mean-and-or-variance-change-fastcpd, eval = FALSE}
+results[["mv_data_3"]][["fastcpd"]] <-
+ fastcpd::fastcpd.mv(mv_data_3, r.progress = FALSE)@cp_set
+```
+
+```{r multivariate-mean-and-or-variance-change-fastcpd-result}
+results[["mv_data_3"]][["fastcpd"]]
+```
+
+```{r multivariate-mean-and-or-variance-change-fastcpd-testthat, include = FALSE}
+testthat::expect_equal(results[["mv_data_3"]][["fastcpd"]], c(300, 700, 1000, 1300, 1700), tolerance = 0.2)
+```
+
+```{r multivariate-mean-and-or-variance-change-ecp, eval = FALSE}
+results[["mv_data_3"]][["ecp"]] <- ecp::e.divisive(mv_data_3)$estimates
+```
+
+```{r multivariate-mean-and-or-variance-change-ecp-result}
+results[["mv_data_3"]][["ecp"]]
+```
+
+```{r multivariate-mean-and-or-variance-change-ecp-testthat, include = FALSE}
+testthat::expect_equal(results[["mv_data_3"]][["ecp"]], c(1, 301, 701, 1001, 1301, 1701, 2001), tolerance = 0.2)
+```
+
+```{r multivariate-mean-and-or-variance-change-InspectChangepoint, eval = FALSE}
+results[["mv_data_3"]][["InspectChangepoint"]] <- InspectChangepoint::inspect(
+ t(mv_data_3),
+ threshold = InspectChangepoint::compute.threshold(
+ nrow(mv_data_3), ncol(mv_data_3)
+ )
+)$changepoints[, "location"]
+```
+
+```{r multivariate-mean-and-or-variance-change-InspectChangepoint-result}
+results[["mv_data_3"]][["InspectChangepoint"]]
+```
+
+```{r multivariate-mean-and-or-variance-change-InspectChangepoint-testthat, include = FALSE}
+testthat::expect_equal(results[["mv_data_3"]][["InspectChangepoint"]], c(
+ 300, 700, 701, 703, 705, 707, 708, 709, 711, 712, 714, 715,
+ 717, 718, 720, 721, 723, 724, 726, 727, 729, 731, 733, 734,
+ 736, 737, 739, 740, 742, 743, 744, 746, 747, 749, 750, 752,
+ 753, 754, 755, 756, 758, 760, 762, 763, 765, 766, 767, 769,
+ 770, 772, 773, 774, 775, 777, 779, 780, 782, 784, 786, 788,
+ 790, 791, 793, 795, 797, 799, 801, 803, 804, 806, 809, 810,
+ 811, 813, 814, 816, 817, 818, 820, 821, 823, 825, 827, 828,
+ 830, 831, 833, 835, 836, 837, 838, 840, 842, 843, 845, 846,
+ 848, 849, 850, 852, 853, 854, 855, 856, 858, 859, 860, 862,
+ 863, 865, 866, 868, 869, 871, 872, 874, 876, 877, 878, 879,
+ 881, 883, 885, 887, 888, 889, 891, 893, 894, 895, 897, 898,
+ 900, 901, 903, 904, 906, 908, 909, 911, 913, 914, 916, 917,
+ 918, 920, 921, 923, 924, 925, 927, 928, 929, 931, 932, 934,
+ 936, 937, 938, 939, 941, 942, 943, 945, 946, 947, 949, 950,
+ 952, 954, 955, 956, 957, 958, 959, 961, 962, 964, 965, 967,
+ 968, 970, 972, 973, 974, 975, 977, 979, 981, 982, 984, 985,
+ 986, 987, 988, 990, 991, 992, 994, 995, 997, 999, 1000, 1300,
+ 1700, 1702, 1703, 1704, 1705, 1706, 1708, 1709, 1710, 1712, 1713, 1714,
+ 1715, 1717, 1719, 1721, 1722, 1723, 1725, 1727, 1729, 1730, 1732, 1734,
+ 1735, 1737, 1738, 1739, 1741, 1742, 1744, 1746, 1748, 1750, 1752, 1753,
+ 1754, 1755, 1757, 1758, 1759, 1761, 1762, 1763, 1764, 1766, 1767, 1769,
+ 1770, 1771, 1773, 1774, 1775, 1777, 1779, 1781, 1782, 1783, 1785, 1786,
+ 1788, 1789, 1791, 1793, 1794, 1796, 1798, 1800, 1803, 1804, 1805, 1806,
+ 1808, 1809, 1811, 1812, 1814, 1815, 1817, 1818, 1819, 1821, 1822, 1824,
+ 1825, 1827, 1828, 1829, 1831, 1833, 1835, 1836, 1838, 1839, 1841, 1843,
+ 1844, 1846, 1847, 1848, 1850, 1851, 1853, 1854, 1856, 1857, 1858, 1859,
+ 1860, 1862, 1863, 1864, 1865, 1867, 1869, 1870, 1872, 1873, 1874, 1876,
+ 1878, 1879, 1881, 1882, 1884, 1885, 1887, 1889, 1891, 1893, 1894, 1896,
+ 1898, 1899, 1900, 1901, 1902, 1904, 1906, 1907, 1909, 1911, 1913, 1914,
+ 1916, 1917, 1918, 1919, 1921, 1923, 1924, 1925, 1927, 1928, 1930, 1932,
+ 1933, 1935, 1936, 1938, 1939, 1941, 1942, 1944, 1946, 1948, 1950, 1951,
+ 1952, 1954, 1956, 1957, 1959, 1961, 1963, 1965, 1967, 1968, 1970, 1972,
+ 1973, 1974, 1976, 1977, 1979, 1981, 1982, 1984, 1985, 1987, 1989, 1990,
+ 1992, 1993, 1995, 1996, 1998
+), tolerance = 0.2)
+```
+
+```{r multivariate-mean-and-or-variance-change-Rbeast, eval = FALSE}
+results[["mv_data_3"]][["Rbeast"]] <-
+ Rbeast::beast123(
+ mv_data_3,
+ metadata = list(whichDimIsTime = 1),
+ season = "none"
+ )$trend$cp
+```
+
+```{r multivariate-mean-and-or-variance-change-Rbeast-result}
+results[["mv_data_3"]][["Rbeast"]]
+```
+
+```{r multivariate-mean-and-or-variance-change-Rbeast-testthat, include = FALSE}
+testthat::expect_equal(results[["mv_data_3"]][["Rbeast"]], matrix(c(
+ 701, 1301, 301, 1301,
+ 1301, 301, 1301, 710,
+ 301, 701, 1829, 301,
+ 1968, 1993, 702, 886,
+ 1994, 884, 1822, 1975,
+ 814, 755, 810, 1915,
+ 1962, 781, 845, 778,
+ 1978, 767, 1738, 1985,
+ 1870, 747, 1754, 792,
+ 1843, 722, 771, 953
+), nrow = 10, ncol = 4, byrow = TRUE), tolerance = 0.2)
+```
+
+# Linear regression
+
+The true change points are 100 and 200.
+
+```{r linear-regression-fastcpd, eval = FALSE}
+results[["lm_data"]][["fastcpd"]] <-
+ fastcpd::fastcpd.lm(lm_data, r.progress = FALSE)@cp_set
+```
+
+```{r linear-regression-fastcpd-result}
+results[["lm_data"]][["fastcpd"]]
+```
+
+```{r linear-regression-fastcpd-testthat, include = FALSE}
+testthat::expect_equal(results[["lm_data"]][["fastcpd"]], c(97, 201), tolerance = 0.2)
+```
+
+```{r linear-regression-strucchange, eval = FALSE}
+results[["lm_data"]][["strucchange"]] <-
+ strucchange::breakpoints(y ~ . - 1, data = lm_data)$breakpoints
+```
+
+```{r linear-regression-strucchange-result}
+results[["lm_data"]][["strucchange"]]
+```
+
+```{r linear-regression-strucchange-testthat, include = FALSE}
+testthat::expect_equal(results[["lm_data"]][["strucchange"]], c(100, 201), tolerance = 0.2)
+```
+
+```{r linear-regression-segmented, eval = FALSE}
+results[["lm_data"]][["segmented"]] <-
+ segmented::segmented(
+ lm(
+ y ~ . - 1, data.frame(y = lm_data$y, x = lm_data[, -1], index = seq_len(nrow(lm_data)))
+ ),
+ seg.Z = ~ index
+ )$psi[, "Est."]
+```
+
+```{r linear-regression-segmented-result}
+results[["lm_data"]][["segmented"]]
+```
+
+```{r linear-regression-segmented-testthat, include = FALSE}
+testthat::expect_equal(results[["lm_data"]][["segmented"]], c(233), ignore_attr = TRUE, tolerance = 0.2)
+```
+
+# Logistic regression
+
+The true change point is 300.
+
+```{r logistic-regression-fastcpd, eval = FALSE}
+results[["binomial_data"]][["fastcpd"]] <-
+ fastcpd::fastcpd.binomial(binomial_data, r.progress = FALSE)@cp_set
+```
+
+```{r logistic-regression-fastcpd-result}
+results[["binomial_data"]][["fastcpd"]]
+```
+
+```{r logistic-regression-fastcpd-testthat, include = FALSE}
+testthat::expect_equal(results[["binomial_data"]][["fastcpd"]], 302, tolerance = 0.2)
+```
+
+```{r logistic-regression-strucchange, eval = FALSE}
+results[["binomial_data"]][["strucchange"]] <-
+ strucchange::breakpoints(y ~ . - 1, data = binomial_data)$breakpoints
+```
+
+```{r logistic-regression-strucchange-result}
+results[["binomial_data"]][["strucchange"]]
+```
+
+```{r logistic-regression-strucchange-testthat, include = FALSE}
+testthat::expect_equal(results[["binomial_data"]][["strucchange"]], 297, tolerance = 0.2)
+```
+
+# Poisson regression
+
+The true change points are 500, 800 and 1000.
+
+```{r poisson-regression-fastcpd, eval = FALSE}
+results[["poisson_data"]][["fastcpd"]] <-
+ fastcpd::fastcpd.poisson(poisson_data, r.progress = FALSE)@cp_set
+```
+
+```{r poisson-regression-fastcpd-result}
+results[["poisson_data"]][["fastcpd"]]
+```
+
+```{r poisson-regression-fastcpd-testthat, include = FALSE}
+testthat::expect_equal(results[["poisson_data"]][["fastcpd"]], c(498, 805, 1003), tolerance = 0.2)
+```
+
+```{r poisson-regression-strucchange, eval = FALSE}
+results[["poisson_data"]][["strucchange"]] <-
+ strucchange::breakpoints(y ~ . - 1, data = poisson_data)$breakpoints
+```
+
+```{r poisson-regression-strucchange-result}
+results[["poisson_data"]][["strucchange"]]
+```
+
+```{r poisson-regression-strucchange-testthat, include = FALSE}
+testthat::expect_equal(results[["poisson_data"]][["strucchange"]], 935, tolerance = 0.2)
+```
+
+# Lasso
+
+The true change points are 80, 200 and 320.
+
+```{r lasso-fastcpd, eval = FALSE}
+results[["lasso_data"]][["fastcpd"]] <-
+ fastcpd::fastcpd.lasso(lasso_data, r.progress = FALSE)@cp_set
+```
+
+```{r lasso-fastcpd-result}
+results[["lasso_data"]][["fastcpd"]]
+```
+
+```{r lasso-fastcpd-testthat, include = FALSE}
+testthat::expect_equal(results[["lasso_data"]][["fastcpd"]], c(79, 199, 320), tolerance = 0.2)
+```
+
+```{r lasso-strucchange, eval = FALSE}
+results[["lasso_data"]][["strucchange"]] <-
+ strucchange::breakpoints(y ~ . - 1, data = lasso_data)$breakpoints
+```
+
+```{r lasso-strucchange-result}
+results[["lasso_data"]][["strucchange"]]
+```
+
+```{r lasso-strucchange-testthat, include = FALSE}
+testthat::expect_equal(results[["lasso_data"]][["strucchange"]], c(80, 200, 321), tolerance = 0.2)
+```
+
+# AR(3)
+
+The true change point is 600.
+Some methods are plotted due to the un-retrievable change points.
+
+```{r ar3-fastcpd, eval = FALSE}
+results[["ar_data"]][["fastcpd"]] <-
+ fastcpd::fastcpd.ar(ar_data, 3, r.progress = FALSE)@cp_set
+```
+
+```{r ar3-fastcpd-result}
+results[["ar_data"]][["fastcpd"]]
+```
+
+```{r ar3-fastcpd-testthat, include = FALSE}
+testthat::expect_equal(results[["ar_data"]][["fastcpd"]], c(614), tolerance = 0.2)
+```
+
+```{r ar3-CptNonPar, eval = FALSE}
+results[["ar_data"]][["CptNonPar"]] <-
+ CptNonPar::np.mojo(ar_data, G = floor(length(ar_data) / 6))$cpts
+```
+
+```{r ar3-CptNonPar-result}
+results[["ar_data"]][["CptNonPar"]]
+```
+
+```{r ar3-CptNonPar-testthat, include = FALSE}
+testthat::expect_equal(results[["ar_data"]][["CptNonPar"]], numeric(0), tolerance = 0.2)
+```
+
+```{r ar3-segmented, eval = FALSE}
+results[["ar_data"]][["segmented"]] <-
+ segmented::segmented(
+ lm(
+ y ~ x + 1, data.frame(y = ar_data, x = seq_along(ar_data))
+ ),
+ seg.Z = ~ x
+ )$psi[, "Est."]
+```
+
+```{r ar3-segmented-result}
+results[["ar_data"]][["segmented"]]
+```
+
+```{r ar3-segmented-testthat, include = FALSE}
+testthat::expect_equal(results[["ar_data"]][["segmented"]], c(690), ignore_attr = TRUE, tolerance = 0.2)
+```
+
+```{r ar3-mcp, eval = FALSE}
+results[["ar_data"]][["mcp"]] <-
+ mcp::mcp(
+ list(y ~ 1 + ar(3), ~ 0 + ar(3)),
+ data = data.frame(y = ar_data, x = seq_along(ar_data)),
+ par_x = "x"
+ )
+```
+
+```{r ar3-mcp-result}
+if (requireNamespace("mcp", quietly = TRUE)) {
+ plot(results[["ar_data"]][["mcp"]])
+}
+```
+
+# GARCH(1, 1)
+
+The true change point is 200.
+
+```{r garch11-fastcpd, eval = FALSE}
+results[["garch_data"]][["fastcpd"]] <-
+ fastcpd::fastcpd.garch(garch_data, c(1, 1), r.progress = FALSE)@cp_set
+```
+
+```{r garch11-fastcpd-result}
+results[["garch_data"]][["fastcpd"]]
+```
+
+```{r garch11-fastcpd-testthat, include = FALSE}
+testthat::expect_equal(results[["garch_data"]][["fastcpd"]], c(205), tolerance = 0.2)
+```
+
+```{r garch11-CptNonPar, eval = FALSE}
+results[["garch_data"]][["CptNonPar"]] <-
+ CptNonPar::np.mojo(garch_data, G = floor(length(garch_data) / 6))$cpts
+```
+
+```{r garch11-CptNonPar-result}
+results[["garch_data"]][["CptNonPar"]]
+```
+
+```{r garch11-CptNonPar-testthat, include = FALSE}
+testthat::expect_equal(results[["garch_data"]][["CptNonPar"]], c(206), tolerance = 0.2)
+```
+
+```{r garch11-strucchange, eval = FALSE}
+results[["garch_data"]][["strucchange"]] <-
+ strucchange::breakpoints(x ~ 1, data = data.frame(x = garch_data))$breakpoints
+```
+
+```{r garch11-strucchange-result}
+results[["garch_data"]][["strucchange"]]
+```
+
+```{r garch11-strucchange-testthat, include = FALSE}
+testthat::expect_equal(results[["garch_data"]][["strucchange"]], NA, tolerance = 0.2)
+```
+
+# VAR(2)
+
+The true change points is 500.
+
+```{r var2-fastcpd, eval = FALSE}
+results[["var_data"]][["fastcpd"]] <-
+ fastcpd::fastcpd.var(var_data, 2, r.progress = FALSE)@cp_set
+```
+
+```{r var2-fastcpd-result}
+results[["var_data"]][["fastcpd"]]
+```
+
+```{r var2-fastcpd-testthat, include = FALSE}
+testthat::expect_equal(results[["var_data"]][["fastcpd"]], c(500), tolerance = 0.2)
+```
+
+```{r var2-VARDetect, eval = FALSE}
+results[["var_data"]][["VARDetect"]] <- VARDetect::tbss(var_data)$cp
+```
+
+```{r var2-VARDetect-result}
+results[["var_data"]][["VARDetect"]]
+```
+
+```{r var2-VARDetect-testthat, include = FALSE}
+testthat::expect_equal(results[["var_data"]][["VARDetect"]], c(501), tolerance = 0.2)
+```
+
+# Detection comparison using `well_log`
+
+```{r detection-comparison-well-log-setup}
+well_log <- fastcpd::well_log
+```
+
+```{r detection-comparison-well-log, eval = FALSE}
+well_log <- well_log[well_log > 1e5]
+
+results[["well_log"]] <- list(
+ fastcpd = fastcpd.mean(well_log, trim = 0.003)@cp_set,
+ changepoint = changepoint::cpt.mean(well_log)@cpts,
+ CptNonPar =
+ CptNonPar::np.mojo(well_log, G = floor(length(well_log) / 6))$cpts,
+ strucchange = strucchange::breakpoints(
+ y ~ 1, data = data.frame(y = well_log)
+ )$breakpoints,
+ ecp = ecp::e.divisive(matrix(well_log))$estimates,
+ breakfast = breakfast::breakfast(well_log)$cptmodel.list[[6]]$cpts,
+ wbs = wbs::wbs(well_log)$cpt$cpt.ic$mbic.penalty,
+ mosum = mosum::mosum(c(well_log), G = 40)$cpts.info$cpts,
+ # fpop = fpop::Fpop(well_log, length(well_log))$t.est, # meaningless
+ gfpop = gfpop::gfpop(
+ data = well_log,
+ mygraph = gfpop::graph(
+ penalty = 2 * log(length(well_log)) * gfpop::sdDiff(well_log) ^ 2,
+ type = "updown"
+ ),
+ type = "mean"
+ )$changepoints,
+ InspectChangepoint = InspectChangepoint::inspect(
+ well_log,
+ threshold = InspectChangepoint::compute.threshold(length(well_log), 1)
+ )$changepoints[, "location"],
+ jointseg = jointseg::jointSeg(well_log, K = 12)$bestBkp,
+ Rbeast = Rbeast::beast(
+ well_log, season = "none", print.progress = FALSE, quiet = TRUE
+ )$trend$cp,
+ stepR = stepR::stepFit(well_log, alpha = 0.5)$rightEnd
+)
+```
+
+```{r detection-comparison-well-log-result}
+results[["well_log"]]
+```
+
+```{r detection-comparison-well-log-most-selected}
+package_list <- sort(names(results[["well_log"]]), decreasing = TRUE)
+comparison_table <- NULL
+for (package_index in seq_along(package_list)) {
+ package <- package_list[[package_index]]
+ comparison_table <- rbind(
+ comparison_table,
+ data.frame(
+ change_point = results[["well_log"]][[package]],
+ package = package,
+ y_offset = (package_index - 1) * 1000
+ )
+ )
+}
+
+most_selected <- sort(table(comparison_table$change_point), decreasing = TRUE)
+most_selected <- sort(as.numeric(names(most_selected[most_selected >= 4])))
+for (i in seq_len(length(most_selected) - 1)) {
+ if (most_selected[i + 1] - most_selected[i] < 2) {
+ most_selected[i] <- NA
+ most_selected[i + 1] <- most_selected[i + 1] - 0.5
+ }
+}
+(most_selected <- most_selected[!is.na(most_selected)])
+```
+
+```{r detection-comparison-well-log-plot}
+if (requireNamespace("ggplot2", quietly = TRUE)) {
+ ggplot2::ggplot() +
+ ggplot2::geom_point(
+ data = data.frame(x = seq_along(well_log), y = c(well_log)),
+ ggplot2::aes(x = x, y = y)
+ ) +
+ ggplot2::geom_vline(
+ xintercept = most_selected,
+ color = "black",
+ linetype = "dashed",
+ alpha = 0.2
+ ) +
+ ggplot2::geom_point(
+ data = comparison_table,
+ ggplot2::aes(x = change_point, y = 50000 + y_offset, color = package),
+ shape = 17,
+ size = 1.9
+ ) +
+ ggplot2::geom_hline(
+ data = comparison_table,
+ ggplot2::aes(yintercept = 50000 + y_offset, color = package),
+ linetype = "dashed",
+ alpha = 0.1
+ ) +
+ ggplot2::coord_cartesian(
+ ylim = c(50000 - 500, max(well_log) + 1000),
+ xlim = c(-200, length(well_log) + 200),
+ expand = FALSE
+ ) +
+ ggplot2::theme(
+ panel.background = ggplot2::element_blank(),
+ panel.border = ggplot2::element_rect(colour = "black", fill = NA),
+ panel.grid.major = ggplot2::element_blank(),
+ panel.grid.minor = ggplot2::element_blank()
+ ) +
+ ggplot2::xlab(NULL) + ggplot2::ylab(NULL)
+}
+```
+
+# Time comparison using `well_log`
+
+```{r time-comparison-well-log, eval = FALSE}
+results[["microbenchmark"]] <- microbenchmark::microbenchmark(
+ fastcpd = fastcpd::fastcpd.mean(well_log, trim = 0.003, r.progress = FALSE),
+ changepoint = changepoint::cpt.mean(well_log, method = "PELT"),
+ CptNonPar = CptNonPar::np.mojo(well_log, G = floor(length(well_log) / 6)),
+ strucchange =
+ strucchange::breakpoints(y ~ 1, data = data.frame(y = well_log)),
+ ecp = ecp::e.divisive(matrix(well_log)),
+ breakfast = breakfast::breakfast(well_log),
+ wbs = wbs::wbs(well_log),
+ mosum = mosum::mosum(c(well_log), G = 40),
+ fpop = fpop::Fpop(well_log, nrow(well_log)),
+ gfpop = gfpop::gfpop(
+ data = well_log,
+ mygraph = gfpop::graph(
+ penalty = 2 * log(length(well_log)) * gfpop::sdDiff(well_log) ^ 2,
+ type = "updown"
+ ),
+ type = "mean"
+ ),
+ InspectChangepoint = InspectChangepoint::inspect(
+ well_log,
+ threshold = InspectChangepoint::compute.threshold(length(well_log), 1)
+ ),
+ jointseg = jointseg::jointSeg(well_log, K = 12),
+ Rbeast = Rbeast::beast(
+ well_log, season = "none", print.progress = FALSE, quiet = TRUE
+ ),
+ stepR = stepR::stepFit(well_log, alpha = 0.5),
+ not = not::not(well_log, contrast = "pcwsConstMean"),
+ times = 10
+)
+```
+
+```{r time-comparison-well-log-result}
+results[["microbenchmark"]]
+```
+
+```{r time-comparison-well-log-plot}
+if (requireNamespace("ggplot2", quietly = TRUE) && requireNamespace("microbenchmark", quietly = TRUE)) {
+ ggplot2::autoplot(results[["microbenchmark"]])
+}
+```
+
+# Notes
+
+This document is generated by the following code:
+
+```shell
+R -e 'knitr::knit("vignettes/comparison-packages.Rmd.original", output = "vignettes/comparison-packages.Rmd")' && rm -rf vignettes/figure && mv -f figure vignettes
+```
+
+# Related issues
+
+- `mosum`: [[#4](https://github.com/doccstat/fastcpd/issues/4)].
+- `mcp`: [link 1](https://github.com/doccstat/fastcpd/actions/runs/7272426093/job/19814531091)
+ and [link 2](https://github.com/doccstat/fastcpd/actions/runs/7272426093/job/19814531284).
+- `bcp`: [[#5](https://github.com/doccstat/fastcpd/issues/5)].
+- `gfpop`: [[#10](https://github.com/doccstat/fastcpd/issues/10)].
+
+# Acknowledgements
+
+- [Dr. Vito Muggeo](https://www.unipa.it/persone/docenti/m/vito.muggeo),
+ author of the `segmented` package for the tips about the piece-wise
+ constant function.
+
+# Appendix: all code snippets
+
+```{r save-results, include = FALSE}
+if (!file.exists("comparison-packages-results.RData")) {
+ save(results, file = "comparison-packages-results.RData")
+}
+```
+
+```{r ref.label = knitr::all_labels(), echo = TRUE, eval = FALSE}
+```
diff --git a/vignettes/figure/ar3-mcp-result-1.png b/vignettes/figure/ar3-mcp-result-1.png
new file mode 100644
index 0000000000000000000000000000000000000000..00963e8349cfd47e4c16622856091b91e0c27aa1
GIT binary patch
literal 71750
zcmb5WcUX@7A3t2t9!f*e-YV^sBs8?ACP{meij<_i2c@MVqS9VOX-isYD=BT-DhdhD
z>$<U82qZ{UE#7EV!})mkWOh
zk-aUhKFz0YKLZ2=1iWkOW6Asa`yaZV5lk-3(~%=f?;fgAzO2A5o-UfC=~lbcx{i&J~ze$@nfid*e7-+s-?%N=SCSFs(s6dV!~qAMtR@Su@S
zx*N@n8#fY`WbXHTS@+vNd=XcMu3)ORJKj>&z@WUM!f@u4uwFpL2dUiGlk0wqi-*;X
zXoz1I(OVL3`v2>mEr?&g#Pu+#f+YUGhvl;+Wu?Xycj_gTBn#n62xqntBV79_|8Eap
z$YFf?G$SLUd)>bMH2-^%TL|re!&epu3$OoJed63x<{$g*pWy$7s?Zg$zQuWrYUO~u
zhf?59&DGzJa{l-EdHD=f5!%XC1)6&4r>DGYnbb|!xeZsebc*wjo>v!+?XOX8{nP)-
zWBBSmk(e{qiTh76-lwOF3ZpM-^LXwtKKkLqudS7j6g(Ol?SX%`3;Wsh>IdD}RhiBw
zWHz#}(>-4NDx_{?$gg+kT#=OQ(;m%TbHqwK)#Gtp4`N7k8M!XRJ`|2E8uz$s2{#T)CWMpJ^l@>EpXWXX4&S|`+CnqO&DHvg-Z5YIy9~w7H3Eh80$C0pbZ2$j{
zc%h@Vwszx<-)0{hY;60K2$N`><#+9xoNIrs&(izLU)~MXG|RE1jYbfrt~yeae|zTE
z@dr2lY;R)72R(wY7^n#gcyPbEy83MU^78WNA|-X4AYsna&a}RHvpjr0sEzC1y?gGJ
zGMM?ONP)vnw!gTz6%>AsR=a=x{F$0XOu_$qL2~lxt5>(yrdlxsmoHy-uk*#L|Jx*P
zDR=JN@tS#k^7s9qzqJ9sq-!a#PmDOsB063gKT53j-L-2MH5HZHfA3FASTt@aW~FD6
zxM*o;+=_+7#4_pr_sf%EWTjTMdwFC%Cr<)<4aZ=Q^H2?Dkul;j68Kcjoq2tx%GA8yot$R
zt`?)Yxw#L=eZro|d(qj`!zF1$^z_n1BZhD(lJ>xRkMxuWk<6@m16;}hBTt-fsqvnq
zP*+ovb$`FybNt_%uCDE!Uzdk&^$>QKu}po?zw8}WPEPBy$BanjO0c!3mWIn6rP%^>
zH8tnommW|lI9GR-l9KZM`}Y(%v$L~r-spxWCM3|_YA>EDbL!m
z(N*=f2g0$fQ}atp!&fy#$du+~_5S-o-N@hHA4LeC{jb|T+DhD9TC0C}_{XIE%GK{rT?g{@vLflbIy(#f
z8!uu-Waye18ymOH3`%ZqZ^v&N_;4xTFMs^_vDd+VJd~q!#%pGUGV^llz(8d)Cl{*B
zpEm#G+Y|S2DVLmmF0P>9*OADRot^FET|Mt*x6kZx$_AIdtA8-^0lrzWsgSiA&!@w^f^n|5X{CLy4)WrpCtoHC{9OkMc#e
zvK_-4_dPQGU?IkGU%_h{cYgM{2ZwTi-=k9rMRnIxPqgyt6S~J=ax8NL%&nr&qO9-g
z*oy_<7w*GldRd46Y!08tPYu_;IbHO9eC2=J&-n7ISFd#MNZNjRcW-SSo5rWNZV%7*
zt(CK5fg7`0hqlxxOQ@0vyDKuW(_t4mokM`7{G<6dZ{K>=smHNYk1z^Uh!Pr5+J0V9
zQJr-i=06AB?x3dxS6`TZ731?PSmFOP8VW4q>({UEcws)Yw6r24BdbU9stV3G{jIcb
z-%ZjIvWtAtX%9oOiUUe*p_Fj!Xj9NX6wi%s--Pt0xn&CNC#
zGYmrc%IeO4C9*&*izbaA0s)v$||2GM7r*7Z9TYvH6;XQk9q4c^bst=c)
zGcP^UHs!fr*2Qyc`NIpd{fX5f1c2efNfOfh?j7CeEw1B@0l%t7moE3q{y~#FANb4h
z@6NBz&d!UYPmR8=7t_^|sM=Nje4TJWSy}nowQJSm`_eYFG&L#q?1_3zKR&>p_{=@dgGxH8l$4KN^h|Y;yKlY6X8va_eSNY8((s6gAv@@2-pw4><|DbV7IZFHiuX-ILV|m-&|=&8_;^px9-b@T
zKEK6oAG+#WL>l{0q~Q|&+ODof+qXEN^Gu{d*vj#KGQzH6b$nIX>N
ztbZ?vWaQ<~*IhmIR7fv#u+I1E=*lgT=(cGY=XV#zYNU91AC_FCqN4JwTMhT1h
z{fGX2$i9<7MaS+!ok$-gS
z^6c5O2M%QYBwT)#7FHkcIc}VP9>NnoHa2E;?wr70t#jSCjafv^O4iobPg0D^t+KE5~ij{5dKIC;J?U;(u{JKNe+
zQRc!(1=f}69EBUdy!`XXgmRS}`sb#;FF!qXWmGppBe-|u^5>mj--;i*i!tX{R7mpi
zJwj7q_EAgG&J#5+3*UOAaSM-j_}Yd;xM8~ZIgAwTyM)ov(Fc_sZ{Ent$eh0ZrT8Tl
zIT`DL1Maa{Od;P%P8bw)#4?LwBg!~-rOdS7A`C_T3o|ATluWMnosH}i}S(%iyh7y4fsqOKzuxX{&q!`0PQ@1A0HY
z{#@pt*XMWlK0Z{B!&0|y@b~%pso_~&tMN6$gC(F~BnuhXWkrk3&(G^(QM}2JS@U`o
z%`h7TvXK1WWn0CYq&Ya~SZy>wqN}0tc#(W+-Ooz=xgxt_gi6jCe+v=7*h5~+IXO9!
zH?DelWyoK$uGSBqQY;qT2}Ct#r~9%o*9{0|@x)1Qv-Dw&mAKx3kK~0sF*RO>F^zNQ
z&SmDif6nDQ70)VVr)FzA_4m*3;{iW=QFXiY{0BE}hWXE&IWr_1X+lzKr9bdgM9W*J
zVbI~EGnP}|b+5I_mcB~tFE?4zgu6uA-$
zGEo(;#%UR)=nC$-{^*6Kp5D1<7c4zJJtrD(_^yn$tr?%D28b{{*q}XdMMWj_=1u*N
z&n`~XU;8%S#|Ny6k1H)Z2{dkCU~rDv9WH`Ovz>s;g8J4ek?*y+{yyOs5(>lmzWAYvl-$OY0xt&`;ZpLG6cjww
zP?l3$TPv$K!)D51!c#@E{Ez$5jmzHNf_j;Ro-##?C9V={xbIQ-jNR<
z^zKN8BxOBXJRwV${c8B+JtfiRzkB*0jrRBN{oQe{!N0oSAa%bHZ=Ro$wcgfj(+MrD
z&dy{Z`AeVJIXKSs-6%ib)M(QZiowL2>-k)MVIop`=+M~bw^^Vn3@Rv)cctu4m42T(
zt`kMavAnXf-VO+8zHKP0cWL~^(AZc&f3DVImsqRe%V7o|Cd#FVpFuS>HK@R+)YW}f
z$1KHIlKoSw>iMglJdv=j*BtmBPRWRVQ(S!=z>1EJE^(pj!w0JS2GKY)_B=3zkY4pI={G>t>8sM
zVsB+-RW!Ze`Q;tQo;`cQd&2E5qDT+p(-?_u;85w#@R?l4Ehn`#uRmHzsj|dl^
ziY9N8!jU86b-sIeo?cn#H@5!)B7rS#oAV4a0NUW@^QwPf;M@8-B@NBs@UX9^=cK2?
z{o_H04Y
zVn4mAz5FF``+H}1cXwwe<;O4G5%~+?TU7U;+`MXOvBhTV^0++chT4P2N3Crqu6YjFMZMv
zoh74w1swyu`a6)J;yYszPEJnzLY)GgO_P5wwp3m$`J0{|q2rC$f0T(4n^{>Y*}taI
zOsdAK(2q>l2OYtMX*KoNKVE7w9W1@u^X2P@CmNf}p?4Hrc$^E(2J1>HpFF9}
zKYXF!)uxz$fFt^$prGK9BS-rB`WhP>^QeS$baaj#J2uZ5UU-sET)a#_wYubyclGK3
z7_3s0^%E!Bw5t~WIbTEM+%&|~!@iSbD_t|@a6TA^RW)X&Lwn5L-rm@56I7(+Q2DvK
zNJb^`^Nsb}(|P_sR%x{we|()N-x7K5XheP^#j@yDeZ7(}3HRQZjg5`L+akWr!NgH}
zxU8E38*ax(7-TDwP6c+lcOQ&2Y)CEUEGukIqS~K1
zGkRMQ!{m>9I&+Pm>MM4oH8yg!V<3OIyzGcTcEH;!9NIlT1
zv^I51+-yHoIYPw)lfXn+*Im5`vIHhcLrWXt%pGLoF?(W8vF8!DYdj90eX)gbl7LFzBLW+y)CkpkMyL#H!buv`Tyb(7mV>^3<6=l2ijjI{e@xgR+%3K|y@ggg9#B1ieg
zwpESi$!nQU+662{smVFBcsJzpSHY?DzuE7Zt&lfmX)3LuFXxD=}HxKIKx%4
zr?cM7gt8BLFl^ECCd5aY%G}{ZubojUmB~ad=cs0KL*|zR5*8MF|?%f
zz}~&LQ&LREoq^~M9KBis`D6Ubq5!HQt)@h{p^J-)gF~*!^(n<Kq
zj}0htpgIkepP%UMHH7A~`Lc_~aR!T(Mck?edI2iB=(TxF0xL+>#-5asmiAN1PD^_aami#4yTdENPYu`iJ@J<2Z0nQ{lSgdj->r%2d&zFVKL$oEZPgq1ipsk}Lx^-%2t#uD`samuQIZrs)54k1X-Fx@Id&?XofaAd{UcOX`
zZhHFksRKnl0EvR4VrEuW&;8?EJh>349z7yRSs)~Rvd3(iiB6wBtyEm*V?mrx?5V)c
zE_{v$x8n(j=8NTzN~b{qvO_1dS)eE<)Chy74jfc6{G&J#8ul{H^6-=(ob^?*lU
zn_?{Ae(iY6lZB~RcJ76c6d%320(h87o$bH@)saoqG?DSn#b8q&v)d3M06_fx{01Ye
zHe+{_k^=R`l-x7|e^0cg$XV>RHFD)V9#{`e9S}=+&ll*KjdgX`AsgT^L`)PLR*3?0
zPn26X(A!=rs8+5_(bLq_gjR$RJX`Iy2`~`e3$+UZ<;vuwr0vUaXkez_gR%YY%ek|t
zMCUJnGp+e4t(ld-_?RxzPkp0@h$PRLiTIK#`({WgxsOa1Qyob$ja?oMBb)Sg==)de2hm4S)zLG&GCS>L07Q#_=m2pA{2ZI9#GJY%95}%9%jL*bhQP>Z?4(8y2gm^
zBXS^qbB5gN0@1;Mm+W?VU-DtB_A+W+`lzupHu80Uf6k*d9v&3hZY=cp{};$|TX%(m}m+;j5eNjM_7fE#}5w1w~A8xzj{
znA^bM?f1kyG0fEL^fY9Z3{DDBWd%E`gD1jsN!N{AN)CNrnGnwsm%Q{E
zvj4eHJfYZnY@{&%=kjRij}RV)?>FtNQ}Xn?w%s@y!qrBOM17C8g2t7Weq`
zL4P*iX30ZN2i_$1so*{S=lf}IZ*NpjZ||zxyxd|$d?}?O*B*DQ5lNf#cVZ6K-?V=G
zShsgw(ZMMkN^Mb<4o0FgqWpG}CO7?LUM?7M*1*{Wd#DFkkTDkq1ZT|7IN)R9WIw89o43&#aA=AXZR
zuJ;>LbMR3s55Cnf`6DWM>{tUGhvJpxk-mX}Cm@89-O9eJ?EHtvpIywltt%u#qS&xl
zmTzlk_b_18+R!i&s0Q-HvZ}$rcPI`JR(1X)K*&irC~5LoB;x$R*lAA%w}UHvV1{z;
z?=O6OM$0q1uwa>bW(e@F^`fvu_XttoSCj{EMDv9&;dfJZHo7HdEqP(WEIs;yV8H*f
z0QUA+IHeZxD=#Ccp|u5WEn7l@fxy=N;tERC-k9BgVewV$g&hAJU%GUuf8sHEbxow3
zo!vCd4MO|*`t|JyeJMt3iP!b_pzQ{R+HVBcVQSO`tme?Uz8={4u6z@(m$nQ14XMi-^(@j#JaVplUQ?z5?B
zJO#1+Q&Z3$k`nAZNq2z3wxWO9;6S2)KigU~kmxbQ#Xep1pwsY~Tz>!hwK?CX%D)$a
z;?Tqd7td4MRD*Qjout
z)+Se?^QNPtnX#II?-a!E)H-T=ZMDo$U?;avQ;arhw{r9IMJ(sik46#|o0W=)i13zo
zLL+{5;X{AJb>HyB*q$})O45F69spt{X6A^ujI+Tfj8ck=m(cc7*q`mSr}zwnEh!@t
z7ZXEl6iyFvZ~!X^D*@(Eb1|p>Ijs4+cT4i}u$oIs&pv;tL+e|=xBO5%k15xdk!#W@
z<`EtZwvY?7N^NbbGPVJj(sxSni6#wKUOantnCXStj=8`8HdgY<>%U;oK{7F0w{olu
z`WHB+3@02_+v8I`C2sR6dLEmMlcXE##;={&zgs|Z_!~B!E9=whAt!>{^@bQ?)a)v2
z`iI%gE9(OPCJ$y*gXcM?8$u;rH}y$V2(wvNte-Jmw^
zLx=P!%vL{7^O-ZW4L#nzXaj=W@f?)n3)jPzto0JJ?Iouroi4BGnk1oQo>av%e2FR(HnoZgghi|)5l+w2-wd@1
zVtKs`NP|ST5IxURK>foPM@L3SZyPlXk-y)0@YOLQ%p0kXvhYfpBj!8rI;CBM<1FA4oerLSxqF
zd_p{3!m=tCCM*aO)a9;s?*zl-Q-W`t$~9g96Nhx(Y3gHayw^H!N`z(Y?qIdep3U-R
zXjaierE<+!Ccs||D*F2RsDDB-jAk#7vUxO>m3h>j#WtSKP)Iq&A?-j>+k$PTzP;as
zjy=3LU+;i$sG*UAw59qC@1utg_wL<`(H=q!BiB=o_0z1PqM{P#o^)zslTimp$3y2{
zq|kiFE(nf`)A1SJZ(U(SoGYN7_p3MAiOfP6Ue$F!KLXB%UUmAhzA5h}i0|Z=k&hom
zx7(MNmf&-fq+Rjw@PPUNAMItFBt)0!mx43mPphko?HDQv%AeHl)9@Zdp6Veo3X!`j`m;WgPYM*3UUU0s2LB$8fP
zfo%dQMc!i!PV>pDR}UA~Cu#C~o=CLzu4Rum5aB$+&p(T8-it=+gTf(lWMy(IL>LyT
z^xdQsU09pgzP+>;E;(K1F{9>gpjQC~fx^pUz7(H%sffsQ)P%tX$~4GVY)p*nMB_0R
z4Yto*)nXZ7$H`e)*UWiJ5iHA^};8pGX`f|T_>&K6;
zeUGBDvp+$Bm;v2NsusrY5VT
z%?Y3l1{LXy3m7vX3n?k74ES;Yb75g&9PI4P^8V=LsE(&7K+T~-1FCJTiL#ISk)~y4y4lz?Ykq$e
zzPx;aF3R4+BReU{*2X6B&Ygm;Mum;|`~}o6v2Vh$w>urI#C-{U{3gdESAWn>JaX7S
zl|Si4Z=M!~fzX*o@$kdZ(|pnJ%xy@UfA2V~SO`XQWv^`xaqZjJE=eky17)abXM=hn
z=IM>!SI4rY9b;qLjY2`q3Qs@YO>z>`_VCdoRC9l%B7mmB17{fc4;^|tJiKqiLDbm9
z#Kdeo+_vflAuk5pxN+Lpxa@orUD}yM>t>943K!&Q
zxL@WkACt)vjP6}|#aUTrf#ji28`tWxqx8NebRMU-`!8Sq24bI9tnpIR3A6^fM-}fh
zL}HU8KIcvgm@wXXUvA4m*yPz?P
zq;Q-UNxi+OSM%7rb_cJKEId>1C-tmgVSZjXKeWV0_H6LFU;RYjZ`z<&Y;TF-wc_Wu
zLqk5Wbs#2~uPP*_qx=KE;tfXY~MLd^0!(ucCD`gR+yHrbd5Vb^iK{+q+nerVWSix%vLR;HCp#Q;=+|pnOssqZ
zqzhf%TVJ0pAl#8lm{nA<@&*V~?K@QpM+XP)@OO8Q%q%Pv227xMB(TYl@eEejszMX?
zt}V19NmrDkEf)k+94I_3bo6Ll1ogfm-^pi4e=?CVsgJ{Fg~D2Nef>^G23WcMzh+b}
z5bFYIQsSh8$(@kbA{#}S560fs&I89rZMb4&#Gr08+e<9+>P4FOnJ1*IDlWgbRz`kz
zw6~K7ef|6y8PiJDsQdTt+v!^q#3CRLW+o;k78VrrqSfNA3Q3qJ7(
zl!cw+jsbXZ?0}?AGeka*v6_(LBaoqqJ*+%EC2PagDLUsikdX69;F#Q|X&qzPSY
z>fRM!UoJT}CKVS<_F@7*hpda?g$rU^SJr;~IPN(~#OrWkZZSk^aV8O*@!-UB6kHJM
zEE#9I`#t}*Rt`k9Eg2>kPwY?^*z`w4M6?sg
z(@{Y4SuUX7-tRG2HkR`Hxt8-^?k+4f6~ZEOjM`i#&(xL?6Vg^yD4`qO7!TF`7V`(#
z_LP;DzNL5dAx_w515?LhN;T3Q)i7qNtt{709+5(<*IyfXX!Yscw8#kEyeTdw#((UC
zq@|^$MKSbZK3-lXN7I0SfCv*}ik$fPc&w*1fh{%OBn19|W5E6ck;I!y&SYg*&
zAMAP1#WYQR$Gl?0SABg_5G#0jU!u&2W@@sIHW5;q1qGVFJP#vt0iiV>1p?wSV1Q1V
zSjMsIhpibAZlk=!hAiA3=6(ByADj#WH}`_H6aS{3c^bPFsgDsfOORAwP)k_p046F*
z)mGv?Z{LyUz6(LlArn
zPysKzy3iyuIr%bjOVQEMJ39g1Y}y+Duh8qD`QZ&Ha%cY46WQz9GFMmOc~Y%!2D3JQ
zBZ^TYiISF-bQ!705xMz*RXqJ;fe{}p_O<&}PoL6IQ=j&*e;xV_`4a@dpjZbnKnQlt
z{N53@{-kVfWzdiKwFSKm6*V=Xi8nN1t)8{8Q1BR2K(WJ^p7Y{_8EWI<@or{)G8=-B
zguJ|Y-@&og)<+Ov0j^MqusXi`PBy`~fg}v{>K+*Q_uH2qgb~Z9%hv!QLAHSX@%ho4
zI{a(Tz@A=pbj;6oCd2Mf$&rLm48aI_z6x4C_dq04QK&J`!m*i0F6S?Xgrgl2>=(;y
zQq+B@+IPTYkxS&|;{*3I{o##=RHg6%Lx25Rkei$P;6ch828!hI2juCR&eeR@@$*yQ
z<9&$fZmZ;j?1z8|Q0q9l2=rv9H+Q_IUXnsVWc|}-=Y>F*rb&&^akJ|Od~o!KKTB9WEA-ejd$UbKU3)*
zl?lpzqE34HN7iGEwfhOufW}Ou3?=xtfVWWaqN4AkI6ieBrd#Uj{b8kS|2lrkQvn42
z*Tx)zXt4Lbyd~sA(_*3$|dkKfG$?B<&JuEN36cA94S}{|or>}olMkY+6
z*2=@p&CRq*J_(tNl=zh3;@zChYY7%&)|GYu@g+t~D_UZ>Zi!=_&At2SQL0^*$U)aCg&t=UJDcA3QMV=LHQ(SH%y?c(p_AzEK7lDIg
zKZ5JJxVj$H`klLbIrrwC2)+-ZT?WVhps4^ygjLsJ?}_+IfXf)AL`0C*aNB9ze(>M{
zT6<5XT=&GpdGYjQ4{dCN5-OuVS60VYUz=RNBj(2|gC|m;HHb
zYpagWwZ!f@KO8za&{(#l6^_4-u8{tq8MNDUql+1Ok%_w(!Z9k;2C
zFMa~`He*QZR`5KX>UBTURx3>j^S8#Tqphh5Lo6%`i
zo6{GAC?q7*LD~7`8w`FXMq`sqh^zG?q3Ka3HA#Y^{pwly3xGkR@Dj7K;5kTfbK4*}
z=;*qLoXzT|hNab2#8=mr%A2OYs^U^n!SB-~&RdtIT|)7jjGdmJFL0p<|E*6YS@V-<
z!PZKzSN6vAD@vL><3Cu5GP|E
zuR=V$`t=jX-n}F_PMHWDzj4(*WZ?(33^17>aRFG(>Y*iPz*u|8xe8g>3cDuRk(1j7%m)ySzY=_o~B8MUAd-ZDfqG#(hm;uBeuK-r|
zd?vOrx3~NQSTb
z|L9moT7_7L<}PWTw8gu}uCu)SiD!##ymcBHE;eiK@jpPd9rV%2iXoy->pU?$#ykR%
zE}E7Y3WkU&IV*&`FID(ig`y$gmp4EQ5PQSrFE^h*{{gfAm7SF3iDH$+n=MT0Azx3;
z*@0S6v)SLXo6Y03nfko$C&|~Y(JX!@mbklur1-l5dddhH3qBjMWmq@_NA;>S^rgEN
z>QdL~@ApV_*U;Ka#Vh=5AvNGEG|_$`0~fe|;Q>0XczUZI2SRtVKA>EU=mM}wA3l6&
z^1OZ8?`$(1i=LjZz@R9TT`31K9oU@~ZxjBd9jKxj-NV6g4+CY|QKxkThJmbJ>r;ko
zpGrjc<8A~5pfW|$aV*ZwZK6?d1^x`F{|vwMySIgq{iRFib27*Arxd6B{CreR=GZW}PG`93m(e*;hQUIG+WIk*m4z2zPblPx-3d7$=k^`L
z1t5n;0mbi(7z+YOEc#zoRuV2>ba_g}z$1%L2aeF1;8Ldf9L^$YAAt#}4?rp_NFUEBl9a9Y!~!xS
z>MraZL?nJgSwNoe;#iGQCDm@TViV0pLY@{OSSYivlQt?O!Y!tD6?yAQNR;CkPT0c&
z9OszK%R)Rvjn_z75HX4Y*|&4vElaUe`XpVsX!pp#$ZQiG&lvO`B9%3b=;+wDGug6B
zGMKP1aArsaBE;?_-H}RW>HWct`nAKfZ!{{TSJ25Jw*_B?FRQ-U28HFN>!uVyBmQht
z%zck;jj6HZQAh+$V9Kc6FlfVw4!u6l-Wbgf!PPK|T~W&Z0F@#yxPBq&VF5TY&Jpa%
zWiYg`DCk-u?pB5Zit-?^m!zY+yVZ9M9|vFzCvV?hiIP&SR-q3e7_h3lk@#dlF)C-F
z(;TPiD!vU!80YL4}uh=?%6e
zEkGTcm#zIEF8WqLA{)xcMN##
zlx{1?Kh@;#>gW)5-&Qe70eoJzE@FJ3$sekDwbrW4AwKLoOD;8AVpdkMUE>Vmi|WTj
z)Qxfl*|JFNXY_{5pvNG<1DHNWE1!nyGWf$RbKtpJQYUF
zq+YlYsWO2zpxK=Puj6LC+1MPO)Nq2DhiTtFazTrYLrxcZev)NhKAdPtn^=1w7(Gz*
z)c!PW8=D6>aGWfthZ7MoQ9?O}nnq0Z4QspaD&II$n-r73fM8T2*G~lpb!N>^@r0vl
zSwPibjoo@tN$KhB3g`6oZ))=I$6>`3rC>vv`pHMa+FGmn8((Z}&9Hp^Hd&GVot(6p
zbEVY!u;4NK6K8hj-h`H_c>&5x`X7)5)ucv5~7bS!1ZZT7dui9C_E`@VrunQE|K5}7T|U#nEgXS#t2u_k@Y}TGFXIF=&|O#x8#=p?`@W<^H~)Dp
zGqd-JJgUu|^x2zmk3)DX5l#lkppQ|b(ec6A3GhE{?JkL?oqnKQD1GQ2DQCBI#nbI)
zWU@Okh|^P_!2(BScZo_uj|5WSQH`8*y9)~u_YxYWnP&{)nEfJcv>I=QZW_)S0%%)U
zXqpw`BjDU3VF8%39|uEVzG|-sS$#uK-N(d)7A5ofiKIh{>5wwior;RY_L&wN1njQE
zbB0W7+ZrjwE{~(3c<_$fw@*BAdfR-X%*w`z_RLB3!VlP(!t$zCx#=bUZx#SAC3(ct
zYE$UXk*FB`{JH#WRfy~uQZ6*#9!Es90rLRaqQvubCFi@M{mnUf;~11Z0in-!Fi^yi
zwKc+)fXrb1s#O*7^lopCmk4*IP*0V?&aAQ#S2!-g@m(j7BGR>9oHT
zWQzu$B^Gg7(x1P6wHT;U;4D@CV7(bp<$g1e&=n@C=e
zCznH0fywQKx(98R=k%z|@{ptBHE!k~jxomAm7cHVqzRYem3ev!nvu4lp;WO?WY>(l
zs$HkJ=p)JO#Tle?knoH=_|{=x`7`C>D5v45?(VKa_y1J$>_WKcZ+{sx?yV1F(8
zzqv2#$v1vAM81
zu(&MLgY&4f$4-9!3&0NaHVGK%*N!?jbsN&Mhoh-*6y4-u5k0S-DNiW9ZF7O-}7|z)~
z=!t*+2mb(SbXeYH7}0zmB~&^-eF{W!n@QcPb_@zoD3ICef%nED9=qIq)n~5Zls-GV
zGF60%wB%ieBiu|>$o1d9&vlnqKp^4@hnTpyxVfN_nxtW&IDz_v60!k3h;1N5;TMbZ
z)rgITg@(cqGC62x?MqLw8->S-HFkSzWoc$c2I?4+0_xO;NX{TCkZ0_Q24P~tG;+g8
zy`J7R6K4-Y4v+cr@$o^i;o#t?CsCb{S(e;7289x?5)eiE9i5hNFhFcwvRv+214BDY
zqOaK){ON&%9>|fqN4ity7Zl7x@*+5H$~P3vhJ;lkd&75lO{tO1o#Yx48XnBQ%r?VYVSPgmV-O3=Gip
z1a4Ux8lEVkW#h9iw10|F48}w__C?UYzmN~`p4!%;$B!R(f%$d<-;AJ%TSK5V6o!0Q
z*q66i?Q4oHMJC_>Y}(&B?p1Y?#`@lKNZ@?dL#U#J?`#-s34JR1Kzd_eGZXoxm^!Gsd+7TqazI3%xeho1WGC)iXbB~CFp=Ko
zTu#c&q`F*!<40AwemE4cQDLT_15-iiHK4bdSk4+wz6c2jz)uu|#fM-;+&&3pW)A5M
zz#?!TeY)^-eP(5O8TBq_R?H_Bq!9Cq#C-4QXlBlFE=Y(`QRIiIZvVroI-uKhJ2}~_
zQ#LanV97)*RJSO4(oG)UiUP-yKtlcL{3_T)D=>Nxn-*c!W8FDBkVXACBc6_!+<5&p
z86hEZ?*k%~Yf*(?pje-nM6QfepvIm_-SW&Ckyy|AqaHFQOZ>=xzPhOX^}O7}7Lt;J
z0*W}{*oQcahc$@>4^e;#=gn4DF1owNr90+6(`G11wY#&MUg|p+65e=5fS=gtXk?Eq
z0SG$2uznR_qVf`u+;9@TPXQ-~yKO1B47SN{zjE!F^(mbYuWl52NBpg3`c@lfHu}r^aaH
zQ`nC_X!-K(zjOAR>+48|(4=?c+}s>WAt^}tugw;}?YKvF
z7Evl?e*WE6dsMYW1QvxYj)m=`%qE@Y!l^Zgu!FFKjq~;uvYwzw%`;G;z=75&kO}CV
zF}H3x*Ju=`XKty?uKPJMXBu;?8IotoxblYnb8~gY@`+V@2uud$eGJ;gg9pKPQz5?$
zx3q-NM7-X`cKFg{9pb5zOt>p(4uW-1EJBziW_#lqm5nA;sTouVhiH&RNmbYv;>>-k
zJ?tX4^I@@D=(ru%olBpv^BkNDmgATgv%TzR!<*k%h#OUH3U1q>3$
zjtWg;1=dmHxU)*{-8+ZxRlu$TVnENO_$??TWX}`6(i3gS^RoTbD>a-;LC*piKo@?=
zJ2Wy<1BWkHvf{}TXxV2`X#fH~AfIr@C-lQDo)E-Jw5{b|NALMJT7T^c>Iu$ge<6HR
zjS_pWETWOnOU^FyHX(k7ldmP`BM8Sn-Elz6#bthZnVo@w*Z}ZLxnP)8eifc_4F@ah
zowT%5`$s%`aXgA;$Dz{VE|~6@Zwc}7IE-WSj{3#Xo!qHZsb!UcD<~PvPlEp73ms}G
zY`wju{~S}SsHyo0FbrRUBxDk1?c_961AOt4(fiZnjg#zjLJ|^t-q8{7LA~F>zpH%7
zEhMBNm<;`T1Q;y*Mj`oNU!T!+9z#>^rJpIW(2GM!sUO_S&er$PZOoY^?)#y^;9YIg
zK&nDvch@#XcR)71yCLYG^44R+caWolKm&o>(fwNrr%2oRXli09B2=o(p&`nM5B-KC
zA<{^VZ|lIMBcHt361wuNB#SOC5GP0Kd$EolFhn?Ooj=b-*FERWyyxYYZ(z(e48%zPpC=+4puaTJ$
zV!w!L=Gq+f!NF100KjZs!rb1>s6~GamJbk)5=rFgO+$LZV
zP``zLjM_uALuO{?7NXr*Cky0i?HwHWzZ2*fYDs73K`?TZd5HPa894Bdi*kZwXUS86
zn($pk^|&u6+R%$&7*d)WpFe*;DUC8TRz+168=dNU%f-}A2ejV*Y~eBl{s1v6lad%_
zM5elRtfi3O!NJh6V+u_5$hHw$UXWC|cTYG!!cWRSp3Le(s#o_WjD%3=+d&G9z1`@q{
zjzK_tuYz#=5hc)4s7MfY2#hyJ99D)G5As`enD5WwTkJ=q9lMxzybz6v{b(PT0*MMT
zJ0RA{b)0vt(G{=1jUo_hpjo))SHHHlhC<)rAf1_yGz)=*<1E+Da`0Krlfu368A{hS
zAls9pkK!YMfs+4))@~r2kFQ7Am8l#3|Il>Z@l^L=+g?e>mLx*9WQ8OvOm=kt-{{Kogb?`z%PX!pM)MIJZw
z;q+`SYd-(BJ?Pn6-uJwIe%}iU3Ly3_SIucpEAP#oEv=Wk^NE{p^?Cwhwfl@}$x&SxBg128?UO?cG(or*L<
z@52F&4**?)XU1bZ(uMtFrIn{vlj?{`u?OD<{g|yFx+=>OHBBoxhZp_+2@Y0=jD~a9ifz7Ewp<#O_nK&
zW|YW&I2_O$fk`~MU+c1g1(%_l2bc$jz!yVG|6mSxq381SHZ(Jfqw?=C73I7vQ3FVl%L*sUqS)b#MPpJgU8PWO2
z@AJd;_hipe#Bbnjfw&&YB~Xme
zpCBDO#?$9PJQ$Aj)tlY{{PCB6g2vd}U`PiS#Opic~-EZCi8hii#
zJ#cb?)Ba=ei1;E#L|
zU@n?P+s^hbH4CZBs;i#^ZUiehYm=J5Eu2&!lrjuOLn5O(HN0cR1DJ$C(&U(FInz2m
zKoCb;On+}LFf1V0h9)LoG4KAoB2;yYV=Kdn9rH72Q9-e2CrpH694jv`$N2%bIdng3
zYioE`PaHC9v_;$=;JHAEG4Z-iAI_+&t8;zjzOGnev-^LWy83!Pj*rOe*=tiB@h5nX
z`dikUu5HG;%L@^XOmC?EqQpV8aF2yuddFl#jKvNsQqscg$B#c3#&<;hN4-d={Xs;x
zTcQ8$r9+B}oL{wykHwh6#8K2f$;Zzhs4Ac!X3(!?cP?wj-xg*oj|o&GwtEo(esyD6
zEi#Na_>lwWf{xF6iSBF`c&O%{E
zQJBEgtDrjFY6_wK?KbZH&$_!K)Yx86Pll%-Z1n*ndiB}n#$P{(=g^F2{kSh#m>A}>
zC_o7l0$|>;8FMG%&DyYz&Xbhc;D?iA3))$)VynZf~FOuvg~*=1wIgF7Y^e
z1B0hSeUttZJF9DJ3*=p*u3g*tXQM6c8r*Va063ATXS#RGHCQl`AV;^bcnOIVlJxRa
zd{`=S?*;7)jWKRY616dZ`9X0O9FDlCsPkEZH^aT#)q)TS0wG`PJ#kIGv4xd*AIb7sa-F
zdgifh+0wp?dv4!^Tr>gC@oMum(R-CNhqzTG6}=ruq~*`AqyzsU4FHK`Z+~*_J^Vn2
z$6F-mkix*w(1DP}Vl5%n^FL)88Y7sZNx_jf6cxgtp}}=gsS(tbzWwML0TXa}R20fc4RfWBe(3`{+r+WC3^@{+JIkchF^Yy}X{PnNYi=;h9fjQBh1DlqQ2JKpz5`4Q#xr=ozMt
zHW<_R#>T}Ru$AZGAHi+JRNL0xjt}X+9b*;mApggBwbi-=2MXSEKpf&?tQV24P@v?Q
zAz!dpfUOHBXrJuq=TbhrAD!Iro3onrWbG;h4pz#mq
zS!?g3aRNy+lPE^+=wRcM7m|}(@Lgseljeh36oNnaWN-Jh;>YJ*e0}8#6vuV!?DF(;
zAp&D}ixf7{(=(Lx)N}A^-tqE-tdi&8Nl47buBOU(LH#Nu7`!*K#JK!LO9>Z^J_nza
z%obr(IIf|d%#ys|jP=H70(5Ernllx@QtxrjIWs6@E0YaHGfr=}Q>Qq+AfMnbs0Jc?
z5THC%55k(7g3GNSk8lzJd44IM5+ylUsHffno*Li)h+7v0S?0f%zx)}WR3g4;#~+pf
zHipXq86$xcX-CepJ!oB=0#E{toYf7J?fL!s*{2K-0Utu_A=iU&uo0v`W*;uB^zXN?
z^|*(p293_?`!@(Qo_9m|ZperENFTk0)+|3K$M6Q1=h_4@yU9@!{#ccnQ`5;dl@!I{
zG&c|{oN)m$+s;;<`K73-B`LPspZ#4TkGvWH9*KCO(lq3)d1=DYOw7!`0W{;YA?dye
zE%eloZ_SNCi3w(!Y$nHa^IesWP5e92EER-HAt=YX>OE-Gik&IXKs8LXH$irck3y+Q
zasB$=fzt4hBF|EY{$D|3fl`CumoRAg&=FKxRdDnR^GZ%cCBA`K!FaiP?)F!uv6CM6
zdb4GBV2Iqc%Za#f?=LmF-iveBf_PbVYi2*kxHEw^tw)PWSFTeG0KV7YO9@|`1#WOncemnS=WiH-V&
z&Nt~jj@9zSI|{PW`&Ijk4i%MFhfWm^w$29rOU|g-y1JW({CZYbmWS4y|!_fHBO;|061b}ZF*a2Jl+)c
z4w*Hj(Qp6tUK1gUrn4RSh7KA8&nmNO8~Z}bQt=Cm28HO0L4u`teR1n-AA09zaa5V#
z9|cvgRr~QZBgKg?<#VMu+KxG*hCVdE*K2*p&*4@;K&@4r^mL3sk)v`l`5o+&O7eQk
z0MJ&vP2Uf(Q!4|t>Z)ko(B18AT-3=`aFepW%X)RI?O1YraL!$qUf7x7E&dIAEa)GK
zPtfHnNicB`MggxC03-1WjN*~`uWlWkuDRJivl-((3mr1(ME)2CI9@M5KFx8`-F;c8
z0U*KrRD7gJD_0HO7oiV);bv@?H1&y>tw@`%9uq4lj4wXI!O!KMjppC4H4rnP8$5$l
zoa|emw6v$LIN4Mz*j$tSzGxT4v8-d0oJXJ49#?*9C`^6#*v*~{PBk;L9pRs9^br*OT|
zGoo=Gi~uis*P%f{S^0{U0Es9%ZG&N?;n29)!bEDTQ$fexJT1-5(c+C|dtcl0fm6lH
zJTpI9HoX+T1~!D8p5M`YDO`A}SV8#%H+d+~BNqo9;2Ye-bqaqSSn$4W#$F~HPKvaV
zDrn^Xd%W{F{lcEAfXNh+t915YayTPi?Amn^-2icQ$$%_K-nA?OUB8RzM>{I;XNxre?G>
zE<+J6=HYAxJP;>%ZXSsPf1k0u+$0D()5MN622AS%PT^4y5xI2JnODaIYNwq2|JVCE
zUkBx=iqZzOrd@jpcZ{E9@y8Efdt&2_NJ`k0dmelB>Q&ob_RbNKgF{2$Ec|9yNnk%Y
zC5wqsgSO%u!fS0(B1w6gzuLC!g%3xIUqyU_l6hcU4z=RG)x$n-@Wem?x{?pZTqS>N
z*E?`U;A^s?TJN0We#%Vk#-ZbHrlv{KVCKz7FF*7cNt`t!-u?kR45bDfC(9m6xh3=&
zT+OEC?ELrIgGSl)>EqAGRHyfkd!3gM->#OJF&Py{sS~M9X#bp6lVbc5O$|RaWW>vlzk)ya7Qk4jO!G-x`<2N$N$l
zagW>FZO%a$rns#CeD>_wlkvwFh+j3$y0t9R!`w+uX|e-yL2YJS
zNjp)$Bo*;XfO4xE^abSZ-PZWf@N5=GK~w~=vS36i!lf%kM%5&YpM{I%>a#1kxjjbg
zZ%qo5Xs>8S6C~PMChTob9Il`Zxd5rEzi}o>N4GL^=#%IRBgOwzz8@wHFZvEW|1wf@
zQ(-pu81bkw>9FC)rAWr|_pO-eDRF&K#q>aOLK+cM%A9g=FZ8)zkA>n~@
z6w2)#EYr9CUkh;TwlbmuG&BlUY&vL-BtcSzZb2=5J)<$YgC}t|-9RuhwfxyE4EYX-
zt31wgvY~1sumP8k;(mk(X1M4pRC>TpfduSBxX(!g-Ty*mje;?SYu766JI-WbV>1ue
z>FjK0U6)f}0a?pPNuS4B=IU~SE!*?n-l(4DEXS8C6FarbtBX(0SoW#3PBx!aI((qO
zvwm)|g;V}+0riUE0oIpS3C{TAfWjrRh;)1Z2Ny=YSzepihYl~5;X2^Zf;uzHEgbvx)Uhe*H3S@GFF0jwH%DW
z%EY=8I$#<%V%Nlo*V+e5;2+LXD~APALiuhMFt-oj>mT8U;!@SZ>pylIRh6QC1TlgZ
zl7{z3w@9jagiz^nOl#JqjL?vfoGlh!awEQxy{60$UwtolyXPuf=+MF2=FldNx}s`jA9KZzIv)T`X(UOKZK<>o;l<+~Dxy;4gt%
z7oJ2hdHI1mDZs=2lYIVUG-WbLt*e2K%dfmi=dmRda<(3Tyg|#)c@deZ{@&i+Iyy9-
zO^c$gR^W_ov0jcQ{(c!Y*x1^$gOlk6KiLY6PQQ5eBSsp1Swl(T_m`&}nUCJGsxU8!
z>UXDzqQn%Hr+ZMJ@mql=0FOQWhb861hp?3+3lB}l2OPEsw2iywLUn!y?+s6)8~%V{
zmz9i3X@~p!PPI<%IQsVtu$A#M&yEL7oLK^Plaon-A7IbLRcl)Tf;c~PbK98ezFa$M
zDb39+qj@-T(VaCRf$qDoxO)~Yr)&mU;m%Ypu;URj~N@A*EimpkO123R4U
z(sO3XrrJwdkks;E$sa7evX#P8ug#
zNip)ktL!9(Qk!adQ;{?wvnC7VegGDJ;60+Sd$*2AG_e9`3hpe1UFgR%M5)MD=H>)m
zIkx~G;bya%Ta{lyKLm_IBKk4=JeBzj;S;M`rl`#x^;1jy8GH&~>qCmm4DCl*8&@b(
z)5THDpop-GW!)a74_A@V#sQ)v8n=x{D*7xh{N6qM29;pFy@t>EjolWf2Szt#xhZ3k
zc_B;s8>Mo%$N^Eo4B3QRPbr9i=Sqy@6576AvP8KgbG<1wmi+hcAI4iyq>PmZ_|1%r
zA|1!S39cCryuE>23G`jy@(tpfjy5yCHy>veX^-aJu&$mBX{0-7^mr7~*RLIs2gpgk
zQ1$(qeq{eQ8#0+re|g;r_u>nbO+o)c6W_B%fa*fpsh9n(%JEURi>?v+kG
z_>+%0h+jRML^|GqVG{s?tV2D)a_Q$gKTH%kZ`GM$RgWX{iHrR}vMh
zel5#}=|RR7qm8uKi&s9N(E)7+$sC%8si3L3IqM40hfZPXtEFMj*{4)4{5jLpW?G@$
z_7EWf6W)FG6cu%tV$S)*L>!FW`!EyCywjqU^OOA8QF5?2Yl@>pDQQ==POcg^jqKU&
z7h>+OSq@WQ)N!auP2DzYF5Y`HKrL!12r7YX2UiJU+7IG>+gqKd1{mhMs`^Z0M~w49
zWQgA*l!E-lCA%2v9^dx6J~j%M=plJEgb|={?R~OVDzgxz&8fp@d;4LEgmoV0%^~K!%%Vd2
zTQGEwO2q%;jXfatgOkc$Nz5Wk-zwsj*I@xZ@ogkxv{H0joOG~b{BRSpZ1OM2RmVd%
zFxQ{0*MpdKrf@^j<@!$ZjyFR$8~oDcdtQ$wTD`RAy2-dIYsXeepYY=+lF+5Lvu@&K
z;?}_qk;2q;T3*7ZNU+cdEN=@A^9>crBdolzJiA$EZf@=mi~)*@X{Q7)7B;qvUn=X&
znXLG^5@;@Sr@h&}?XiPg^ws4z4}#Ssrw-&;iu|GLkskF;yVF*Y
z#bMMg@db!f!X84F)yt!|_nU<{@fRoC@j242U0Ds}f61+*dd*ZelWBQ+N#hqkVRkym
zbm-Z3e`oapSx?nxm#3_>xvkTe1f|H(6`q~#<&MuNh$kHi(EG$a(=f5!D5R>~Pt+}8
zs@#caTjzWH+u^#TsOIKFAtG|cTY^6=D{+)R(0`dwz#&Z*Y7VSO{UlmQYF2BXwM}3u%*QXrgb>cy)e5$N_PiagrF(L8#3YqVki;q80C?&_1DhOaLf$j0np50Mp6McNhq8hk@eowwd9Fj2siD0eI>Le{R9
zl1?s9$q(ewwE-*)R@t{^1g+jH8cis=t^dZB#oiSlO
zxGc*z^Kglk4X7QU*Ep(R;#D~N*{U$<1`t4P>9sujF-R2!|OP3U1
zX%MSfIvVm{(iw!?$?v%s#U6jIG>416f19waXGk7qe)9=M>fBQvnx3k5;aH+keD>0>
zxNVUP&w{rQe-0>Wg%@-RXVW`rv&eZHzb^t4c)TGnhVhXX2Vc>DKCZuIn1vx9GVZ(t
z(Hdp3h8>%xy6NSL9X5{vZt6@EMP{QD)Tv@_Dt$V;sT`l0ihd59J>p~Xu~+Ad+hSR
z0Jq(mv>3Ij%Hwh0=i5P4!sgQ-P2Q3eTKU0S6o;VA7h;hUUvjt
z#Y9IJ>xuW>BL@qPK0{v)?ZbIY%V^8;Fofm%J?SX|X862oSNIr47ZhY^+65<*R=WF4TetcQmh}Kc#6TI*9Z{zAuS~wgF4_KSXM1oK36BEG#V=jYmB|6PX=gbO6_O6nX
z$LGG;Crv#YHTc~gRifAvZC$5$eQUOz`5iCjP2}iZlE401<}UQiA&tCyUAwLi06WKn
z)Ya*5ghn^t|KsZ$RtMM+mtc5DJQ~wC)E}bmOT51q5uM#&VVRgt==}6Z|465k@WNk4
z|MQO_gUS<>GRWIA{P1
zu{lKVf~}br6r(E>h5Gf(Y1wuQa;e6L3%?-7@q#S8(QL`HCe|NJF%BL)l~BTYYHL&f
z`a$Utd1Y*bqQM`+xD2-y1rMbtTZb9v3)RLKo1!R3uTrJU_C-AsQBQ7Ky!-3w-rxrl
zKE!NMqjv)SxkG8H*G_(BWjR6Pf~Q#7XOzo#mCF0i+qX8mb`^;ihUmx#pR_0-dQK}V
zHiUK8f86e;e5-vIj}Bfl^!?XSq`N&FDoA28gsIuyEE42iBF&*erwiJZ&2GFtSeI0y
z#-uHIrKrg4@K3z}QDZ5C!n4C}!wZA8-?VsE&mnB&rQKwTCUdZ8axlfJn75uR4!NE@
zH=1ck>0HSGvtH2lbm@5Qm#tM(RjtP6MlvT2Ovk;G1&o>CC0Xr()B
zsre7^R;W?ax+&L~hgtY*pIHc(zRb0{8-QIkY=^`1C}qS7h^k_vqaB`vHDKt~iAQcJ
z0us`dn#3v&QqtcVf=&t8j3(oXP+GsRlM(vw&wZ%`wa+Xp7&?PoK&|o(s$&R>t!hq5
zPJS@rU7Zt1#VG*4`mwWAYoY*B&m(?wbqL4;szi1{0YXEMdkoS?1@bQQeBaTP^RcVXhecPLCEM2C$t~Jw@vG0e`1_QPZCyw2=dp@IJN!FS
z?&@iHt^r$~{i|_{%p?qY6%+%jt8T?F`b7SYoq@++nm)LC)f$g$!!Ekeu=Elo$cO`v
zECv7Gqy0f*^!RpP)sNK6p2JC-8V@7CQ~T(a#_=T5{|T+D-Q3W|H&AT@6^^Ix$M+AL
zi2m+06Vh?I@Pm6J=ENk(9foa!VMy&NT!^MsO0ycaAZ~}_yOnX*9h-YmyCsjy@s{(p
zeC?GWBVqS)?z!|JE@zIvXgwSx90z`*u;(p=i8Rerss%<=qh~O-~cNx~TYN$ErwV76syB*ZYv^
ztvOkjeXNO-d}zxGH2`xsKYSmg6Z2Y?nxKR217^nWy3d!FcuUUWU{4B0qLTeX?*{No
zvdL)#nkCxNL?t!To0i5=rm5`AT^IUq7|slsSkXz>cK>}Ln|SQ@l|#CFf|Ca`dwIiK
zvvQala`gwR&QLHvG<@b6x;{bxD%s*h4jrOg;UlQVDzW0$`GXh8WI<`?b+3
zzEHpc=n$rV0AyI1>uHSpPLzZe3?AZw?CghFl7LgN`%1Al#I=`S7(BsYZlH^4Oq
z?Ym^CFTVA9`8#V>7S28?U9HW$Awy#*p!=7bfLsR!60Krw3-$g&-ZYq5MseW3g^;!X
z$^^YgU4aM>LxS45Rf)Psr(eWRMWts!Dofx+f+&^=Ev>95)_8NRdhLmvugN}Wo4@Z|
zjM7G<^=Kpe(P=OB=T`@Ig6$Qw>=_p;%fe(8w<12N;4q#Szl)XF$Y3F$JygoC}
z9g<*9j18_me)ETb<=phNuG_JVo$-l@BGJ+|PeL;?pKSd30(G=Si44<8A(e9C313Pk
z>HvXbs3OY?`&>rIUjQn>r(uk;)GO?YJmy+eUurv4
zy13HkBt9z~PW4GKB!aEw+p>R3lE@KTjnLa!8Pyp|x1L|8MNRNsDYg{h&u7g#sj|iL
zL%sU5KG(m1FV@zO8qOYD%2?Z{#lIvcA<;tU^EhsqVxe{-rC05=pRdD-sj{r{N4ujJ
zS65dT7h9A=z0!&^b8_NMoNe}UjK_y7y?B8p?DFDcvELRJ7T~dl?@VQ-N0)~HTa#Sa
z>Sz=7Ql6hZM(T85r2O#ksQLzKx!xHqegK3#m})gI;uuY9zB|K<*Ds$E=P1EPcw|M&;(RGVl*>c2^qq)^i`
z)F+Yb@xTF|#IX%AI^VodaNY8
zk-6>VgVqv
z2&2`nBSE;>m;^91cJ&LM3qJH2c7eqT2M_eO#{jqDquv1~0S3>8)Wz!alOK@t$zc>0B1u_fV@37FfTs
zaV|fzwD;2uw6p*mpDW7va~KvqrJq!Gv1*?hwjzT5EFNz7+xvfZwkoXIh$X
zTkI943z0}s^896<68cn!3gtD@pvxueEkalH^|{z{{8)T|iw
ziZHosU2AaSzA}HTks@z@>9*XV4P+%<*|DBQ;!^K5=n%KeQRNo&Ayye??=X?%=HoLY
ziTHRcPR0$y+QzsN<<-SEWN|{Qf6Vv8nP;-23g{t3#5YNI;d+N@T-#UI#sA{5+YeyX
z(O5izTMRh9AMXrAhot6<7kUM?Xm2vPwcudK=N?uU?7SbjDFy_@Ti8@IZD-M?+*Xb|
z4+YL{JJ*fV26-pRxAPWrkPn4h;(R+FR*e_p=DiA_0eT9e3u)hrUGoCXjF%KT^qo9B
zD&ZSp!g?4iC+Jed1MU-hG=FTHeE;|*glR4Fpy`?*QpkF>hU1u^@(u+xyznjiT<+i2
zpxqRQdAqIZLyj9iu!9<_2XY(`^PJ$%C&0#d>*eWf#v*RYubEkVF~`=7Ec$C}N8XT2
zBax)uLU-Qi@&2d)IdF#F$-k~DZQH?dajB!=haf&#GIT6m&?{K8vZ}P>9%U`IgNK$#
zl1m6;6+3D(ArpnJhIOr{fHp~$^?ms6gjH8;fKnLSdh1vv!j0VORV@#Q3NQ%+<7YU<
zxZArenqOW`?z&B^%K66g@RLCbOP+BYtgZ&SLDG4TrKETR2<|didKBs&@u%W*+|Nu2
zVQ$V-ea53L^?ncR8g3OB9!RCy_dV1li-;}3WxbkV+_g>@A18`QOZz6ZW
z!H6{11K>+V?@9vz{TxWGePwtOGG`N*sxcT)&7kM=yJ=2BdN|NIfI?P%3HWkU@9
zq}A0kvS(j^psh*x`<@g3L%{qzt}{m+Jhx20UgNinHTBnsR@=yE#j76A+YvMLNl)J^
zLZzfv^rELh)@cj)WcM{Qo#?pWNylzqc`U2Dtr$Y|OPRD?TdQZ>^&Coc=);!qjE!aM
z9z<>Qm*$ipgF!{&&o$GVIA2#kopEi|Je@%pu!Scna&#@4QUt{AodBVJ>D`WQUD!Wx
z(kU}De3*L9wd$FowL}EqKpBX~-~kvGs4H{X{(2$EOPDXPIb_Kdcqg<;iC_OKXc_7R
zsdV|r>$U>z+O-y-b;JtZ2h3f?@HLb#?jNsv_l+p2lLx_X{APQVXncOfo?!pE=BPpc
zr#ehakS__U1?rt(>hqCmXNrp&xF-CMe68aHt9TxPeO1<;p0mrWO9D`
zwb%bB?Q>S24^OTocEk2=dUNRAyT3T~=f@fpE4rG0f6T`ekjy2^uTJyj^4mQoZP*8G
z(u|0deAf_OS-a@^OE9VK}_`)0^PN;
zxHw_qL8ldf9|F=?lpWA(;*6zI$oGqxnK@b4P_dylcAn1xiW!rS%-KeLwh9{GU0T=A-m#DR?z*gq+!;QY0@gYn7Zq4@X}olEFk?1u9cNRb249(_%)S@l91}(Jvp>SCEP)oyjR?m-tzvx7GOpD-h3{^
zXp*ylW~)D-MV8Bz+GZaW#HXYv8rd+OJGSlhY~aDA4$g9J+upu`*rb${*Iqj(_6hzh
zlifew^Zfbiz~f#u(fRSgKUY?s_4jv{Ns=j74V?lYZp{1ceSCFSy1uKypTO%FQUd28
zQ8KQ5sZa4mPJyx!b0Rkn4>a1G&n<{0zOSJ^+d4~qu{5`X^$tWQ7{=PlWGA*u
z2${MxeCN11$>h8TwbIn@;19@>pq9JKex`GWUNnQ}&h$WQzAZDa_jK_Fqr<7zdIJUz
zxJLI+JO$})?A3tmfyF(((AVBRJ!E!>`-atXeC^jMJ7ikXeWQZQ_;@X_^}Wyk^`3zT
z5m|1U@a63v|6o7yTKR3$r$$IEaj*tjl~%Z^vJJLrqZ`C0Fukcy;Qc=FiBoMZCZDT#
zdUlq6;hjv%qT{3~W_uR`C4#^}x|?u|mEDrB{3IT7!|!Jh?OzxZKd}3U!S;%u_&pdQ
zV3O&sd(i0`IhP^v{Q0yjiG708{RqvGEpI+b)9M9O1>4a(Ahj`PEGslFw6^N?zGUWZ
zXu^AXT5lSixWQ`A1vGmoqjP2O=fvEu(+N4|Th|nr!2vzK0F}Nlou(+O0ki2{x_45I
z%C|~OEiU*}N6C^K&W~FAQ@`;tl?f;cQd#F#h?>$&ohW%j&cede+9{@>4z?04`jdlg
zRx6{sesN#K(l>NT{+d)H?g7cqG2h(E6*BSidXj|;*mlzP*Yfgb*!ve&Kfl_d{g^6n
zU@tR*jFnb(CQEdf?$*sXJAQ#IqpVCcfip)VRN=Mp^gX8HZ8swa_Awe}JH2C-aG=7P
z)ac*To-_ZoHWwBaqC5W!^(w4}wG%xG`N
zCbCRK-|gL{{C<%2Z=xTR3De1%#g5dKSzT
zu4;Psa2ORcqWHf3cfd%ntmP^Ro6rhbB!Yz4l0ImF4M9$yXAR)GG&ZMuIRb2|Za~iv
zj%Etn?{`lp!)s!C07k8}PFJ-v^X90P02&So4p0^
z&olRVJxIZVYHVNFG-jLq9*Mnvn3Z!Eu8`2Tm2~q-x*?Rp%j?$Dnl0`glP+Z#x2NJ?
z;57zvtpAKeRGsl}OwyM+c84jytl#~T0D+g+pl2N0f$HBqaxNDbp;rB+_i9(j&b_t{
zG2K3XlzsqJ6LX64;Z)B*`=4ZnsCbOp+f7N-5v7Fw0Oeeg2@P}Tm>UdS6glFXhO~{W
zYkB%*w$C2$d;2yRudiZM-7@X%|G87yIqPylJ>j$9CgH~TSL@QVm!!#`tV0AX9oYJB
zzI%S+Yh#~Mk{1=8LZb?8
zzjo7rj?t_n1etA_EW0P@-T;D7GVtNy;?npv7WB5=uW^ZC+zeZGUkkeBsHx4lY{>C^
z$u!f1S;GWT(q`qct2yMg?L^5Uq_NR&;i0U2@;S_w>0fk$
zjanbUtH>Z3*?>^z_3ATHmvUYMr-f;n9rxCpfN2^^af8h6FKYNDEBr(?pp(Vf|
zcS3hat?absAi@Gpn3oUV5k97>_p;etPk#vZ4v#^-fS(`w8iF_uJKTc0`K|tYm%=>b
z9M(6t>dia5Vy@H6Mef~2CQf!=C4Vv9s+*h$w^SJ0+M$
z0EUoy!HV3G7|U9^-S3_lA3s62%^&p{%84aTau$@rJNOw7b3NVocq;bKUb4oS`2i3b
z=)5SKY#h}rt1oGX+$3I+>h+i2uXq}(d~_GWa_-mnPe7rv${kX%=yj5xp1dz(ABoB0
zrs+;Qg7
z+73$VY^yk{mH*!nUAMw8AIU2uOq=8))rx7q-+%s`3k_X5H*%}Ktn;$i1q*vUpN!v5
z2!+VwOfy=4H516?-w>{S$bI5n6}pr#wg-~~|L|r|cCQwgVi4rNeo42&wQm4SpXS=c
zr9n4h&1Y9sc*q2hd}j;=59c&85(u+VdmmjfRJEJ<6?xr*G3s{gM-gG+$b{=05fQ?n
zX@pXABofu3===!P;Td+k$*Tcit+u
zHaeHyV07+#wfCNApf=>h9P3XRLziE#ujSHJd(XXVZKO-PtJ
z(AW456P6{$??-<3Y``XtKWDX|@wE9{K|3)Lddfp(aSGAVO10ADvQ@4a?q6@zr#H$kOA-nWA?rkH(?(vfHa_&9pSI*#
zWIQM9lNi_bOa>3^D@#S%69oqfpe>-Aq8_9E7-LBN=I3L^Xva{v!`9Y>j+_VgZ}{Wp
zJ##vAzy@1MlCu;7zkJV6*3i4VE5%GE@RT6bA@((!jeptm3=$p`ZRI@ZIxde?Zk0jE
z6|5TYy3vrI1#T((tpQ~$=4Z&6O--NUWdW?XHWu`PI{#(i*n`u&>&**J3E7OP5ul`z
zaY;C~`_`p!qYJGsM9pZiyX_V$x%8D6{XUUsnj~Gi3=>;{y~dK)OvPn(l_^iDsA-_E*z%iJJr0kSXTm0Cs)8?M!@o
zYi>sV)&*`F`*WKIUH&)BWUmJ!^ICoDC|f*A_>01pnaAH0c17=kf(OEQ5aOS|nuh~>
zJK%DM-cDLZCT#YKH8z1o`i$*KL}xJr{1Vr89k5O&^HbKABAs|hY9SqjC={d`4q3|D
zwa%Wk|9~OsunbIrfO4#VwXX5
zMWf}SZqpzGR!$X@_V(zu${Omn9udzxlsba7m%*w6Il3teF)lPpbrs{0#LsAaj>*i)
zG_ol7T7y~b-Kss_ZXwK^^z``(`!vNHwVTZJBhsH`0Z3Wxi(OO$cY3FKii*Gb3=~sa
z)pvyYZ__0YLradsB}7inqJWTnpAfbz8MG&LkUaAAO}>&R2R}c{Le8HS`{^#~L%-M7
z^bQ}k*k6`~=IuPHjaXtAmFl(U2(tNe_v%&N-{GU)Ds%=!T`vnFMdG3vFh)_N*7@za
zE2(#RmqgyKz)AKSTfJ*WitW1k1T;bjl2RuF?uCF(A@$+`!DZHJT|wysi?1
zS}Ho(^9PY|EM=`3k97Hkh0yv?G>8r|4F@*qW~da6pH=Ofn+P*n(O(D?*o1};@9UzW
zCV3c~&?l$Q;wMGk?^-OxKIyz0SsVtF_i&Cu6+X#-%*>lg&iUHG^V#RPIXHHcS65ap
zovb4Ar~}|KWY%89ZVL%hB4P)y!s>H4d|I{9HT6s^EK*KS*3r!&5dm9_A3@jvu_Btz8%W-S
zv{JGnZfu-)VtU#V?LHm>%i*9;K~pLj+S~VCYg6sSWWc;-3w(#^AHOC$uuw1&>@H&W
zoHibij4S&6n6vbvF_kjX4K{~Ge>JM5YX!<@ea#yG11Bx22L5R$4kZyfx&
z5V`L9`-o?`D>(Xf%#(}Xx{qkfH|i2w9JQ1p!a&`8SwuPpVhCBkAfL|Sx^5l3=|!<)TAo(hUH$(^pc|
zZ?GI`v+P#fyF2Iol`hplImXWcrbQAf_hIo|MY4Yuds;-%4YQ(+_vgRA%z$fb5^4c#
zb!1w$9g{BBh}MZ`=G-%&ysNm=6=K7k$1m|l(iE&|$Vu~Za|5+YDWbU_favQy6UV6B
zEu8%J6%B>J_85;>+9~^7Hm0s)%--?!3JK31|9)Qwh=yx}B@d}GP?m%7XGz)4|5{vJ
z18pAgZ}fE4lE3@b__tWQ%dp;oLIw`{C5}PyqOc65L1G2?`)lY)qeH(JUOG<$Gf50G
zp*OMQ8m+-ORENIb8Joe8(yK=UCas7@dpPadsu$2?DTrx+HHaHD>0
z3y(~9PfujzXUcJ*EkVX~t~IGz<9W}fFl^cA?!|j{l>jsGgA%hXZoQzY0~}%{FF;|H
z{yI9MI_Z^HG91g7FK@?IBm5k+HBwgP&NU;zij8*)ZCw2*B_u>e>P8Rqo5kx-E}puF
z_#Kc5OtEUpZfC6ZyxS_f_3YiO`;Z@^W)-)R-d;p+VeN$|w@8^M1p*vbqyWALN0wC7
z*i{Ot;K<1Ou`3+%B?tWnuYmjG#6o6kID{t&Mw@JLw7-M1TrwF1)Xx}>9e{MN6XF^?
zJWt@lwRQ|W%(edt)pI@$`(0L&qi=vun<*
zm$L{DNN=!ceA;#X>ppdJvCd@bmMf1O
zS3pofe?&Q-Aqhr6&>(6jgoN+kHLaOx-9EP7D_!OPsUEjqf=iFU#ibXZfwqw1hSZ3v
zX0MRHC-0D&=40}|id&BEF?BXL+=>nan@&IZy#dR}3&`hiu_1nT&6h0i%LX+KgZyuU@UMS2(%<B-JsP3K|=HiHSt49P3@1ZW$@tG)o*E9!(SrwpSRa
zTCuZLc-DYbQkKxoK594EH}}!fkD(w6I55Ux@Xl#ddY7MM?>3}P>wWg-O)!E*W`3}2
z>T02;Z9KdwakLbY6`!Y;WtnR03ktqhHl4a6`4OpI9`52*>UjS>JX;G
z1K0!vN%PRq5E{;~B8AiITZJ5MB!w^uMxCe`!HN4+EdU0oKpp%_i0#0MI@iWTT1Fi4
z7wq0{#dq{V3<12@o5et^Ov2CteC_r(abn?~XjOar(2(`o4*X$sTrnvSdMb;RZm`@#l^-UjG;AKz)S63
zx=cUo+>&rJpl@vxZ9}5J|01|+$;5peqPvCfs5m}92**BRA#kcLpBkWo|u#HA0YvU-=LK8=5yLZYjtFWe`kCx`Q~=ESiK);`(xd?ZPUx$
z?08MmfoKkw0X*8F;$7Ub^Nt<70)Q_oxWxZZ@>g%g_|%Nln^ifFJww-_gV`1?iUC0V8vj(`0Hvj!`{mVFFB?zr^*}+ICz)r4$U?%Ct0KvVl@RlPE
zP7vY!D2LrGuD)9*u#TCsyCsaPru&dxZ+y#kb3q-hB|_(~l6O8@T8xXh;RCdiWch8Y
zKbn2gX&6OpM;yt+7&?GJoAnXV9QF4c9%rH2q%qc^A
zQFb5rtZi|8pt^pnaFcZGt2Al5om?YrK!c9r}IWv{fu@gi?Cp5oxs+x<;;rM#)B
z2^j*w(?b@%Vd!gBA~giM)i=0mdJ50?K+G40;}Cl;EIK(k0(V2{Y7Ib0ldjV${(pT~
zc7Uchcz9mp;r5DbJlWHom7b0jYtA8lgKd+c$hk=uF_5#XWum6F@-+Oy;TWKfj$naH
zM^sF#X|ZjKeDIRdZ492jM2-@VKMDH{5v%^ktH??#5l|^|l%%9E7G3yh>+-dZBH+Mn
zN4S#ki;cf@F;@FLXYrSnpHr!hk&XnC>4`*Q1_@rP=2g$u6%d>{EY6Z*i|a@3#gem2x|N|ROqJl-FF7s#{JDW%47jOZ%?QHbT3wd8(9pPgv2OLKTB8&cfp&&d7YNN8AQ83(FLrQ)Ze~0JpXz{|c09Cq$&k#X2(!
z@$TP}x4-aOd1}TVV(->n&%`jdbK2v%&q?O0S0ex`F$Avt`GciG7eVMFN%#wvz2Oyw
ztQ3?PgOv_H4yN17Coh-Q`#}7E{@gk6=azp9Y{@+zrVj1``+|aj)bje>Vt`+8&hhbn
zZ+L3sY>FoahMe;nWDuv
zxT(G4-gsFiJi4Y(iedT_5Fs(l+!h+}D<%9%vSRb({9435tb=|Lv)?Jiz)6#0eqnAZO
z64cWpQ1ju;dGE&P5es(dple4w)CWk2pP8GpbZx+VjMxWHS0~A>QQt8USO3-}_2LB`
zaj*F&_Pe!xF|JMJ#)JlW7;i7}>LigG3WXW_R`KfdIPwg8p27=`(-K=p^Wc00D5>Hx
zp!MuPfZt?f60Pn>UAzQH&_{5Kzv;lYm6LuF_77zQwve-+H=vE8eUCDxD_8+d)4}C#
z5gjJ-FnEAmu1{Zv}-#$w)3B0(@$4ZpDbc0h*dw7(5mVb@?my
zr^>v0=Jy2m3_ZN4rA77de!u(BeJR6ZvhOu^tUJAzF|(ziks*gc>ks7nv~OHSFH>d-
zTagxHIVI*4n$BI%v9O1}3bqPB2#&7yL|aOC{{!tgUX9*obk_qGGU9|Ic9!UhC2-0=
z7lSb7$H=Lo<4D}v#>=gf#@t`9dr)SGiL>bW|Il>Z@l?0}``E|cd+!m2D3yeeS@vki
z%$7t{3fZIVkrqOBh@z-tW;BeHiXt;qRz*X<>pb7rufLwB;hfL7@B96}-q#2x8^hrH
zpL!mZg9H-0>mF15ZB*K4p#Q^VPd1r^-UF*1tvODl7?rG|{U0J{Z=Ux@$ue@nztyD{
zzDPRrg?wR)+E~W$k;9z%BC@-%F%-9+aC6g6V)2LuP=xxbP&aPoJaYtJJ&MBLIZv{v
zB-Brhl7=#sj_goSn8yujL7jnDSnl}ajbmP1{DecKd#zrh9`X?eu(*H~IiV{A$3XX-
zIMbDEmuhvEK_${eDF_T7!~gA#`8aR{oxJc4Uwiq_QyX!i!n8L}m80ni*|qyYep8i^
zysL`~*=6A5zEKoG269p8sRGmmp<`^Ov@tKerr@Hx*w87c=%IdgWh
zpG%{8HbgRH=njM%s7KkF6q$}wRk2s;+#|(xdnR?JW)lVWzwWovV(g;(hX>=K4|ikOa{6{ptYcPzZ*
zaXr10`CzkPBaui%17HA6PR|}|q6+$BAUCO&`&MuIn8vadV*Jo=pC@W%yi8#F8j}ki
zqzhS4zM`+2k1fvQvVgi+w3lh@u=fr4JVJW~CQPGXJpZ;gC{HjrD_lQ3
zF=-M78qXRhFIP!Y*5O{B$~c7X;U3gUFEv-0z6d>z4~nDqXAr@5ZSA=|skhdyXL#h%?;g50>Tz_D@MmW(=Xl$)6ss}v=!X9#nRcgDo*!@Gs
zSF^G@IZyv%ju5%_BED~WWo-=w)=kTD%VOgbpO^m0I(3etUIcj^#-rHl(T$8hi$9{3
zv~Ay^S4SGWZ8od#{`afv-Kz>a)<24lv-Hcv)Y$}oWQhrAk{$K_!Qkgb3fbWrcd36h
zH(21^&HNN^eS^^#2}7U5Hl}mZ_GTPp@uK&nT6XfQGW#=Kpm6jClrH~ZTrxfae?-B^
zUxq4*(KYYC^otF!MI14*%|Ds7?f3^j_fryfp&DXh+9s87xAs%rK`LOB&G3w{oWJ){
zYkwCtv
z$x@C>IMar`zyJJ+XcPrejlW7(CC~8g;8+0nN?wZe@#7OI@gw_RhrL3h`vp$uIy$Y;
zAERJl%_Ne-^7^trWg&G8L<1Mum)CZ(e)Yp!H>A43Adcgg+U?7K2||t!F8yw~L>?s^
zI3Q`J^v714mTgKl8nVLyh*e;US?Ooe`a|#`!&oqyqUFK|4?^$~Zq`wugu;@F<*KfA
z|MTa?)$MfA*F*LbRbbTx3J=4?uJ+%NwcTWk;ux4Ek}?6K^{
z*`a#
z^exH5cDpJ3<`AL7zkQ9elj_wS@)>c3N8|&fRixcJ%#+JI~`S9*l>U#Mo&I3jk^5~tISy*N|H8Lj#?tBrv
zN0&D2-FY?gQ_}esMsawyV6a~Y77X}W{@nRovYKz`pXRDj>so4_ZKln%xMVJ3!Fx|3
zVTG`qAhEK)hjD!biIxS$t*$lfg^O1^UFca^=`3ZekAM>esQO_re+c8~+4Q6Y_*_7#
zq#F+;QK2~YkafAW{H|S1c#&UeEb;C;IFlE7A{TmRoW|?Ip@l92@-a3>f9`e(HD;YB
zy)~~#-yDqsy?sK=!w#0XBO@cQS;BrU-1yA@L58jkA5+87XH9Sg{Bohy2)xlB;rQyH
zPGH~)){5i#FX~ffJO(#?aBVm8@^Cns3v^=3Iwg!k3R$pct2S@AkvO<++mJb7u(3lio+Wc#9*BJn=NX8a_^<(CF+x
zgl@igkK-n1=n~N$l`j^?CeTW9&R_MsG{r4F-`;c<%}-=2l4qM@TwQ7D>AH|aXe-4M
zj;!82^QZ+!KAJzoCI@9qDGW9zed5!};G?ei=uErlDiud79Q^1&Bq!ElYqs}?7Lp4N
zDE~2D9AWXNLlLtZEeeSuL>xMy#mav2(}JLX;O!?%y}A!Rw5szWS!LKOX5(fV@9)`;
z7*=sa5{YoEGn}~yC9C1BlH{N_+Cp+?w!`D+7j@2Muzlpt>qF0g!u`d~y}vQ?Y&iLG
zv-YV7_ydDjc<#%ibYa=gGF%_OOxo=d1iLOrXOz>x=^YvaoUyoyob&$vTDZHc{HX)&
z-`QVV9*tr5|Kf4$!hQJGq9p0GN26j=Zff`{b4HKo+;A9Us5+E^=7aLCfakmh4|S=h#0&2S0V6pIpAaOFw&2}@
zGU`-;192H&0RoM1Y&`I8QS~E`KtG6SRd1V5Zgm8%c{_N;rYU$KyX%Vf2Z?Yew_*9u
zi=Ivo-lI9%l3ZyGkm1PVu(R-xMxiO6$g<~F*bUHuKY2vNIEF2ElAXeUts{Az47P<%
z4k}Ar-xI>JLwWuOuEGODhc3aBs0&VDbbCB(Ip+OW+@FL@i6RXfV|@dI&Hbz_&NWT1
zY%V^egS|fLIWP`c1a~BxW2jqnq`h_T|AvsLVF2P@OHCz(cpKfU>zGPRZ~qDY6JlVD
z)T{MQg~>8Kk5Hy}BP^~1ynu&2(kyMa$8y4|?5d?MTTa|mdGy_7`gL0`(L$}<$n;J5
z_D7~B8HP&QTyjR=5tMfXMkQ4{)FMn7GvoSQYNw{BVW&_Jq7UGMr*RluWaB`@2p8+2
zZQIrB?$1)QvVu4NoV_z-#Ik1D3+ERUz&a^+GRZR@)bb$5e8=v|E@+};?(ci;`43ps
zpt0d`6KNwGtr$8MNaozSbsP{s8XBnYl7AdBH0&mGVc-x#W*Rxv72)HIIp
z)K&T^suMS%=u~TOes%@y;Fy52Oa>;Vsp*gyHVhyQO-!t4Vc%Dch%^?)sdW(!j!&<>
zrJ3qofC~-S?GX!TDmoVhyduA}5nf8pzag8tX&+A+s>IWWEku5H_eq|2Gh`FOyN)-K
z%uoRxt8r>Y+1exmOVtP+aSbGJ{^N~lAWZeql|osLVn1)^$?8nqtsG1>=r1s=DIk1_
zk0+*XsR_AMo6xZ%(u;5V38q)F($C9jVKU>4eyG(-m!3R|@j%W*_S|Rz6AC!A|Ej7S
zJvlLHU5k&tJ-(*uo+z^_FX>9c|1i>X{G2Fw(Y`(OsvrvvHVzT`u*Q-+H)^t(5ISsF
zTvcFxsVthUT?PEt#CYWFxE*ke|JnDDW2h@qdJ@l`YV>p>=gV+DC)u>p8Z<)gVppol%Ry$iUPI%I)XYh%SqIuvQFJ~KTHnI83Lf0L7!cevO%3rmJ$Ry>)>#;F>)Q?{>R{RXpln>nLtV*p_)uVf
zz+6Bf;LMV;aOB|*J^P!`yFlzuc5XhEclvbBvp;bLJz|kC3AQb>bVnbW`G}|=o@seh
zB>*}AXzQ}GW2vPF1Hz|uhV1WRf!?bO0Htkxj9wgLpxBkwXbSl+-%*<68%vzsx^Lgd
zstb04ec0YGcfw4e29G4YUaB|@qy1Ed?B9}@wNUYLjfqFYEtoOkK&GsK^7$XWetM{N
zwK(ow7Tx#3PlBd8YGUg8UBxRlXvJ!GdF{ZsB|q^_@>>ugK(ZP0f0Iry6`!tK`Tl)c
zTE8e$>OKo3_5KeI=1x29WoPQwtuIS_
zaG&tnwr!gLmuMMVA=FTf*g*b0xv@C;g!1RZ8^7uP4bRBY8%}JEMlwXvujs8gNsBvq
zZgSUB^Z
zgBS8AQTAhQGf(B+Ect?q3nV0+@rBzt#pUFP0#fqDWB&gjK*6~I6&aGtAEB;?KS11S
zCt!YTA4rLWc;pYtKTvr%IXHkAQUxWOC4aL;#kyrn1zG{jI>E?*GyCC-mq{cNh*WH!
z(WgD|eao@$FHSCW2Qy|AA{W;uQnY>%--wBWi)xJQ;Olsd5@KawAndqSwN_LOJuZsz
zm`_~(c82~v9;SS8(dloA;%r?~(dger!{;zp=jE={$m`P(Z-IL?W&|lI;W=t5=60
z0k;dNI4PX1>e8WJIEcc-25)3BFa~4!rgN!;`+L`r4mAz$+76N8J9kn9cd{huE^U8HAMe){)1ITC)~hY3ZmE5kSSO~CSe-xe=0k92gUdH}Kc8yup?#K`dQ!)HR&>eM;o(=@1TGL*22n;Y+4%Zh;=EaP*vQp50sY
zjLvKbS$>*n=s}unqerW-O~F$;rf<03%PRDibGD+AlC#Y4n>VPse!=nuA~4M1RX`Z;
zO2Zf{>O%)HDUb#J?hO2MX?$`@!fY@&`WQMW$Z
z4oSuzjnq4-xw)QBP6VdCTefVmvKlg0zapOqM|puD&;6>7FA$nB^WsLu6cO70tC^Y9
zJFl6`Nj*6X2(P`JYL6qpOl!+9{$Hw5-^oz~H-a$xMGD1y5bF4XZpuax=t$dLtT`W={
zHXSu-GLnu$E|*)*#}ull!X!$kLK9}6`liEX1X;|mW0Zcvcq?98Oj%700dvUsh&1?4
z^6ike#qbgM{HfsJF_5ebZ?_RZ7PUjRs0>qt$Qc~J5V4Upu#ktN*v&DNf-R-J8siUY
z!Sb$
z-q97wqT4AsKjuD)5JGn>9domjDkT=ld1Gfr7l)CiAckupA1|Pb1STh2_cX-Fs0Z<*
zK=<=)>A9TU?lLLdF>l-U`}xj~7&Of7im{RW
z;D1ffB<7WnS1kSr(XC->baZj>z^|L!mE2Dv5%LJ6@BC;`L2bwsrP3B-D=IF=(^WZq
z()E!&fV`6kluCw{Oflv#`}nMDsp6%+M!aiiX7}vg-7DTZP&0Yua}RrUpkWUau$|mno?3l=rY*+OfYRSGi$hiKijYIClcL5
zLPAzpkVyohP;~LtxH!67FOT0l!OZeYGP?1oKe{PoW?^bjG$lXo&*j!Gj(m*xz3C?U
zTDNddzF~R+`?sXdPf)Xvr-+)FRaHt$d2qb{jiFt9Cp<@8n)_v}G>`!f%&YhLbEKBg
zsUCq@&sLFqFM+@|G_?3RAF$zK2^Cp>#@JOnc|ZFN8GI9)
zAur)Tgh&+#xJf3Ofp|wH58!u@HH`Gx2dG}R=A0+D_qq-a{I>&Y3i)fW>tnk=_xG2&
zz6C7&5g?3!fyfz497dR6ZXgW$tc{LUW{DT=cb1DjvFFq!8Ql?#Tej!zh-4VUMXN?n
zR9NqUxos%mP!z~$kEWXh#MV=tZRY{cZ-
zW^wjLGJ+~#GKEc$O&4-7jk7YmOV6&L+?){#`x~P>0ro1FJtbz80e5gZr`RdIAUk_09Aui3rKsf&<3K5dX5?1k8URs7!)Iym-jljh_Hlls=I8h
z&=T}H9dUGc!`&tg83l%1OI6F7F?Yo!3v+MX3P&Jvt$j0`yVfh4nl?b(866vQ
zKWS&}Sn!O^6$z1z9|aZe^7FwM68O5W?gRW(jZ1M}jAWC69iXuZT@{^P4n=$=w?jTm
zOVAulm>K~c0(WS7et!Pz*Ndn*!&tW;r2hDD{S4*!JH7AA(|x>eIpmsZk4C2!7fW}T
zh?OObxp%KCtmmtW?fEWtc7ZSdA%fIu9Ek<+@L!y>qs}YA*}IfBvAfSw*$^3TX^P31
zNxp;N|J^%sPn@{wnb$SIuHNc-E2;OJ5$^0;G4F=L#{xCR6js*|B_&g(t8@ub^KD94
zPvp$rjWV?55t+12-O|F3(BFgPilx!khbj%(1pg4u2%x^BShJUK_p1th%&V2wh6U?g
zpRW{_qV}7(I3b1Yd`a&k3M5=IF#Hrn&-N6G`|;qr;fT!6M)KDu$24V&Hpu>N8pb{}
z&?-%K6tB66vx}*n_MaswtBu4WG_zO`ZW7)@v+FpZjI-zu8JamY4s_ll&;~rfvTbN+
zK>XNpmUDM=@8NA-r$CGe#CZ?a3cj2&kEne@Vygq$jQ#V5cr^5&`a3L&};D$NKZ!9O8adN4H&9$~Kh3-qpz^65WqzY7mKGbWQi&
zxyXC_WL~K&x3bp7KlrN(75Ns{*3>7sGu8-!Z>TjzSpB7z-YjqHoLmm5)1EF!_wkgv
zYsa1`r5Z07=+MhoW7#1gAJla2%BR_Qoo(I#6^o1AEiDHTom1mvrr%&tPLV_-D;vPlX#jZns@w4hx76HAKfCx_xo>5G>+5XcseuRYygUU2O(|~ktjR&5Tg&+n#n;jy&md5g{j$|ex!he^VM%LUQA3C&RdC{XQOm&a
zFv>;Q{cMCoceP8?%;o^`ACu!kL~lS8-+=Psr&2G{!CP8m*tysx4b+e&%6pSaD}0k_
znkV37`>OrtoxR2fso2w&Gi>OS4@MIlo(PHS_0p=8sxhU;@(}px1{rgbFyPbP{tMnl
zkcFU_VXM3}$evi(l_GjPE5nL9_cW=eOFe06n|kM^is)EZ65oRV9K{I?>;Ij+VFqe{
z`BDTvTVwa`hmDQOwrNxKEH>4E;K|AY>~HWIkGy4X`EGyDuLC4Ja4zpx%6+Vp_->puabA}OnDg@T9@U>zN5L=hXfDo31SXy%1#V20+z|u
z5-l9c!N##{;6h^aI*EM2Uy->nUZm;
z&#aXB`iMuxAOECcPQ4IUC?_m@1Dv!En4w{y2Uv^g$z=*TlNgbj((8B1lclz*IU5~1
z^s(NHQ!&PQgM~J^=%`SsA(g#|>t
zgkkT%K9{IQwsP7P!bpwrF%P=LnC17-rgMBvK##L$*Df1tzHLd^*dBx4y=RZ0WeAr(
zl8YgjXieuSVeK!mRGt7=UYlkaf4%mEFjZZeQ0^A5Jec9XYdOEu^LpQ3#j+$hUgQ`J
z4svX*)02p#d6d4DB>}-B8fn;Sv0yPeF{DxXv+wQc)^FCvex2tPyhx-F6^yuQYfl9Q
zQEyV34dK$KUM2W>>u(>&t1t^<{rGrIz=Uq#bo4h3ZxeqhLq&VYM4
z78MwZ;(!nCpA`j8W&S@GV6YhO3z&Ut0eE$-DygX0vZaUi$h!()BcsWk+eJco}`==
zT2kD9q6HiUc=p^?i)U`yt$G#tF4gb@1V8S;sl_x3B&@>53227UMH%%soJ>L#j1!wK
zS|1$37_G6evfhG11ZaiIw2q39urLUxEV^3C>yWXa;(#p!>Q~%!^ka;SNhnbnsL7FY
z*wpA_+F(XBUc6JX=61e5*DL`p(qSnMxEbt*LWsdV
z&9ELS;{tOH*hpMD
z2Gy92KjOH@bi-ElX{EfH5LFl9*2^~tYz9uKqC!t3dJdOo-Me=W%>_@NITL!2rYZ`z
zG~g_Yi^;6#5e2alX-B~L{n_>sIfl>=Xd?j-J^8mujtS1n;_gsDO-FbAJ2*U(cMcyq
zl#)z?2-|fa3q|+$YN@)u$b;GsKi^OaMeemdA=UUdFQ}EZe4L!}ivq|toCO)kP}=91
zA_}RJ#6H|Uu{|O^lzG0Grn>cT0IyFh?4BWHQNB=-kZ61Q)B}XMx;)AefBCnHRD`Rq
zP@@XHjx+ckN~bDiGeF%s5df_=?Z-zbIdCsAk2;SS5;)Z$IR`mr096Y|jj0rA4X*6k
zhnS~Q