### R code from vignette source 'jss5361.Rnw'

###################################################
### code chunk number 1: preliminaries
###################################################
options(prompt = "R> ", continue = "+  ", width = 70, useFancyQuotes = FALSE)
library("doRNG")
library("fastcpd")


###################################################
### code chunk number 2: mean-shift-model-with-gaussian-family-setup
###################################################
set.seed(1)
p <- 3
mean_data <- rbind(
  mvtnorm::rmvnorm(100, mean = rep(0, p), sigma = diag(100, p)),
  mvtnorm::rmvnorm(300, mean = rep(50, p), sigma = diag(100, p)),
  mvtnorm::rmvnorm(100, mean = rep(2, p), sigma = diag(100, p))
)


###################################################
### code chunk number 3: mean-shift-model-with-gaussian-family-custom
###################################################
result <- fastcpd.mean(mean_data)
summary(result)


###################################################
### code chunk number 4: variance-change-model-setup
###################################################
set.seed(1)
p <- 3
variance_data <- rbind(
  mvtnorm::rmvnorm(500, rep(0, p), crossprod(matrix(runif(p^2) * 2 - 1, p))),
  mvtnorm::rmvnorm(400, rep(0, p), crossprod(matrix(runif(p^2) * 2 - 1, p))),
  mvtnorm::rmvnorm(100, rep(0, p), crossprod(matrix(runif(p^2) * 2 - 1, p)))
)


###################################################
### code chunk number 5: variance-change-model-summary
###################################################
result <- fastcpd.variance(variance_data)
summary(result)


###################################################
### code chunk number 6: mean-and-variance-change-model-setup
###################################################
set.seed(1)
p <- 3
mean_variance_data <- rbind(
  mvtnorm::rmvnorm(300, mean = rep(0, p), sigma = diag(1, p)),
  mvtnorm::rmvnorm(100, mean = rep(30, p), sigma = diag(1, p)),
  mvtnorm::rmvnorm(300, mean = rep(30, p), sigma = diag(100, p)),
  mvtnorm::rmvnorm(100, mean = rep(0, p), sigma = diag(1, p))
)


###################################################
### code chunk number 7: mean and variance change model with gaussian family
###################################################
meanvariance_result <- fastcpd.meanvariance(mean_variance_data)
meanvariance_result@cp_set


###################################################
### code chunk number 8: variance-change-model-mean-and-variance-change
###################################################
meanvariance_result@thetas[seq_len(p), ]
lapply(
  meanvariance_result@thetas[seq_len(p^2) + p, ],
  function(thetas) matrix(thetas, p)
)


###################################################
### code chunk number 9: linear-regression-with-one-dimensional-covariate-example-setup
###################################################
set.seed(1)
p <- 1
x <- mvtnorm::rmvnorm(300, rep(0, p), diag(p))
theta_0 <- matrix(c(1, -1, 0.5))
y <- c(
  x[1:100, ] * theta_0[1, ] + rnorm(100, 0, 1),
  x[101:200, ] * theta_0[2, ] + rnorm(100, 0, 1),
  x[201:300, ] * theta_0[3, ] + rnorm(100, 0, 1)
)


###################################################
### code chunk number 10: linear-regression-with-one-dimensional-covariate-example-summary
###################################################
result <- fastcpd.lm(cbind(y, x))
summary(result)


###################################################
### code chunk number 11: linear-regression-with-one-dimensional-covariate-example-plot
###################################################
plot(result)


###################################################
### code chunk number 12: linear-regression-large-variance-example-setup
###################################################
set.seed(1)
p <- 3
n <- 1000
cp <- c(300, 700)
x <- mvtnorm::rmvnorm(n, rep(0, p), diag(p))
theta_0 <- rbind(c(10, 1.2, -1), c(-1, 8, 0.5), c(0.5, -3, 0.2))
y <- c(
  x[1:cp[1], ] %*% theta_0[1, ] + rnorm(cp[1], 0, sd = 10),
  x[(cp[1] + 1):cp[2], ] %*% theta_0[2, ] + rnorm(cp[2] - cp[1], 0, sd = 10),
  x[(cp[2] + 1):n, ] %*% theta_0[3, ] + rnorm(n - cp[2], 0, sd = 10)
)


###################################################
### code chunk number 13: linear-regression-large-variance-example-variance-estimation
###################################################
(variance_estimator <- variance.lm(cbind(y, x)))


###################################################
### code chunk number 14: penalized-linear-regression-example-setup
###################################################
set.seed(1)
n <- 480
p_true <- 5
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)
)


###################################################
### code chunk number 15: penalized-linear-regression-example-cp-set
###################################################
result <- fastcpd.lasso(cbind(y, x))
result@cp_set


###################################################
### code chunk number 16: penalized-regression-example-plot
###################################################
plot(result)


###################################################
### code chunk number 17: penalized-regression-example-coefficients-plot
###################################################
thetas <- result@thetas
thetas <- cbind.data.frame(thetas, t(theta_0))
names(thetas) <- c("segment 1", "segment 2", "segment 3", "segment 4", "segment 1 truth", "segment 2 truth", "segment 3 truth", "segment 4 truth")
thetas$coordinate <- c(seq_len(p_true), rep("rest", p - p_true))
molten <- reshape2::melt(thetas, id.vars = "coordinate")
molten <- dplyr::mutate(
  molten,
  segment = gsub("segment ", "", variable),
  segment = gsub(" truth", "", segment),
  height = as.numeric(gsub("segment.*", "", segment)) + 0.2 * as.numeric(grepl("truth", variable)),
  parameter = ifelse(grepl("truth", variable), "truth", "estimated")
)
ggplot2::ggplot() +
  ggplot2::geom_point(data = molten, ggplot2::aes(x = value, y = height, shape = coordinate, color = parameter), size = 4) +
  ggplot2::ylim(0.8, 4.4) +
  ggplot2::ylab("segment") +
  ggplot2::theme_bw()


###################################################
### code chunk number 18: logistic-regression-example-setup
###################################################
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)


###################################################
### code chunk number 19: logistic-regression-example-summary
###################################################
result <- fastcpd.binomial(cbind(y, x))
summary(result)


###################################################
### code chunk number 20: poisson-regression-example-setup
###################################################
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:n, ] %*% (theta_0 - delta)))
)


###################################################
### code chunk number 21: poisson-regression-example-summary
###################################################
result <- fastcpd.poisson(cbind(y, x), epsilon = 1e-5)
summary(result)


###################################################
### code chunk number 22: ar3-example-setup
###################################################
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)
}
ar3_data <- x[4:(n + 3)]


###################################################
### code chunk number 23: ar3-example-summary
###################################################
result <- fastcpd.ar(ar3_data, order = 3)
summary(result)


###################################################
### code chunk number 24: ar3-example-plot
###################################################
plot(result)


###################################################
### code chunk number 25: arma-32-model-setup
###################################################
set.seed(1)
n <- 500
w <- rnorm(n + 3, 0, 3)
x <- rep(0, n + 3)
for (i in 1:200) {
  x[i + 3] <- 0.1 * x[i + 2] - 0.3 * x[i + 1] + 0.1 * x[i] +
    0.1 * w[i + 2] + 0.5 * w[i + 1] + w[i + 3]
}
for (i in 201:n) {
  x[i + 3] <- 0.3 * x[i + 2] + 0.1 * x[i + 1] - 0.3 * x[i] -
    0.6 * w[i + 2] - 0.1 * w[i + 1] + w[i + 3]
}
arma32_data <- x[4:(n + 3)]


###################################################
### code chunk number 26: arma-32-model-summary
###################################################
result <- fastcpd.arma(
  data = arma32_data,
  order = c(3, 2),
  segment_count = 3,
  lower = c(-1, -1, -1, -1, -1, 1e-10),
  upper = c(1, 1, 1, 1, 1, Inf),
  line_search = c(1, 0.1, 1e-2)
)
summary(result)


###################################################
### code chunk number 27: arma-32-model-plot
###################################################
plot(result)


###################################################
### code chunk number 28: garch-model-setup
###################################################
set.seed(1)
n <- 1501
sigma_2 <- rep(1, n + 1)
x <- rep(0, n + 1)
for (i in seq_len(750)) {
  sigma_2[i + 1] <- 20 + 0.8 * x[i]^2 + 0.1 * sigma_2[i]
  x[i + 1] <- rnorm(1, 0, sqrt(sigma_2[i + 1]))
}
for (i in 751:n) {
  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]))
}
garch11_data <- x[2:(n + 1)]


###################################################
### code chunk number 29: garch-model-summary
###################################################
result <- fastcpd.garch(garch11_data, c(1, 1))
print(result)


###################################################
### code chunk number 30: garch-model-plot
###################################################
plot(result)


###################################################
### code chunk number 31: var-model-setup
###################################################
set.seed(1)
n <- 300
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:200) {
  x[i + 2, ] <- theta_1 %*% c(x[i + 1, ], x[i, ]) + rnorm(p, 0, 1)
}
for (i in 201:n) {
  x[i + 2, ] <- theta_2 %*% c(x[i + 1, ], x[i, ]) + rnorm(p, 0, 1)
}
var_data <- x[3:(n + 2), ]


###################################################
### code chunk number 32: var-model-summary
###################################################
result <- fastcpd.var(var_data, 2)
summary(result)


###################################################
### code chunk number 33: huber-regression-model-setup
###################################################
set.seed(1)
n <- 400 + 300 + 500
huber_p <- 5
x <- mvtnorm::rmvnorm(n, mean = rep(0, huber_p), sigma = diag(huber_p))
theta <- rbind(
  mvtnorm::rmvnorm(1, mean = rep(0, huber_p - 3), sigma = diag(huber_p - 3)),
  mvtnorm::rmvnorm(1, mean = rep(5, huber_p - 3), sigma = diag(huber_p - 3)),
  mvtnorm::rmvnorm(1, mean = rep(9, huber_p - 3), sigma = diag(huber_p - 3))
)
theta <- cbind(theta, matrix(0, 3, 3))
theta <- theta[rep(seq_len(3), c(400, 300, 500)), ]
y_true <- rowSums(x * theta)
factor <- c(
  2 * stats::rbinom(400, size = 1, prob = 0.95) - 1,
  2 * stats::rbinom(300, size = 1, prob = 0.95) - 1,
  2 * stats::rbinom(500, size = 1, prob = 0.95) - 1
)
y <- factor * y_true + stats::rnorm(n)
huber_data <- cbind.data.frame(y, x)
huber_threshold <- 1
huber_loss <- function(data, theta) {
  residual <- data[, 1] - data[, -1, drop = FALSE] %*% theta
  indicator <- abs(residual) <= huber_threshold
  sum(
    residual^2 / 2 * indicator +
      huber_threshold * (abs(residual) - huber_threshold / 2) * (1 - indicator)
  )
}
huber_loss_gradient <- function(data, theta) {
  residual <- c(data[nrow(data), 1] - data[nrow(data), -1] %*% theta)
  if (abs(residual) <= huber_threshold) {
    -residual * data[nrow(data), -1]
  } else {
    -huber_threshold * sign(residual) * data[nrow(data), -1]
  }
}
huber_loss_hessian <- function(data, theta) {
  residual <- c(data[nrow(data), 1] - data[nrow(data), -1] %*% theta)
  if (abs(residual) <= huber_threshold) {
    outer(data[nrow(data), -1], data[nrow(data), -1])
  } else {
    0.01 * diag(length(theta))
  }
}


###################################################
### code chunk number 34: huber-regression-model-summary
###################################################
fastcpd(
  formula = y ~ . - 1,
  data = huber_data,
  cost = huber_loss,
  cost_gradient = huber_loss_gradient,
  cost_hessian = huber_loss_hessian
)@cp_set


###################################################
### code chunk number 35: well-log-data-summary
###################################################
result <- fastcpd.mean(well_log, trim = 0.002)
plot(result)


###################################################
### code chunk number 36: well-log-data-summary-median
###################################################
(sigma2 <- variance.median(well_log))
median_loss <- function(data) {
  sum(abs(data - matrixStats::colMedians(data))) / sqrt(sigma2) / 2
}
result <- fastcpd(
  formula = ~ x - 1,
  data = cbind.data.frame(x = well_log),
  cost = median_loss,
  trim = 0.002
)


###################################################
### code chunk number 37: well-log-data-summary-larger-beta
###################################################
segment_starts <- c(1, result@cp_set)
segment_ends <- c(result@cp_set - 1, length(well_log))
residual <- NULL
for (segment_index in seq_along(segment_starts)) {
  segment <- well_log[segment_starts[segment_index]:segment_ends[segment_index]]
  residual <- c(residual, segment - median(segment))
}

result@residuals <- matrix(residual)
result@family <- "mean"
result@data <- data.frame(x = c(well_log))
plot(result)


###################################################
### code chunk number 38: micro-array-acgh-data-summary
###################################################
result <- fastcpd.mean(transcriptome$"10", trim = 0.005)
plot(result)


###################################################
### code chunk number 39: micro-array-acgh-data-summary-whole-data-set
###################################################
result_all <- fastcpd.mean(transcriptome, trim = 0.0005)


###################################################
### code chunk number 40: micro-array-acgh-data-summary-whole-data-set-plot
###################################################
plots <- lapply(
  seq_len(ncol(transcriptome)), function(i) {
    ggplot2::ggplot(
      data = data.frame(
        x = seq_along(transcriptome[, i]), y = transcriptome[, i]
      ),
      ggplot2::aes(x = x, y = y)
    ) +
      ggplot2::geom_line(color = "steelblue") +
      ggplot2::geom_vline(
        xintercept = result_all@cp_set,
        color = "red",
        linetype = "dotted",
        linewidth = 0.5,
        alpha = 0.7
      ) +
      ggplot2::theme_void()
  }
)
gridExtra::grid.arrange(grobs = plots, ncol = 1, nrow = ncol(transcriptome))


###################################################
### code chunk number 41: bitcoin-market-price-summary
###################################################
result <- fastcpd.garch(
  diff(log(bitcoin$price[600:900])), c(1, 1),
  beta = "BIC", cost_adjustment = "BIC"
)
plot(result)


###################################################
### code chunk number 42: bitcoin-market-price-original-data
###################################################
cp_dates <- bitcoin[600 + result@cp_set + 1, "date"]
ggplot2::ggplot(
  data = data.frame(
    x = bitcoin$date[600:900], y = bitcoin$price[600:900]
  ),
  ggplot2::aes(x = x, y = y)
) +
  ggplot2::geom_line(color = "steelblue") +
  ggplot2::geom_vline(
    xintercept = cp_dates,
    color = "red",
    linetype = "dotted",
    linewidth = 0.5,
    alpha = 0.7
  ) +
  ggplot2::labs(
    x = "Year",
    y = "Bitcoin price in USD"
  ) +
  ggplot2::annotate(
    "text",
    x = cp_dates,
    y = 2000,
    label = as.character(cp_dates),
    color = "steelblue"
  ) +
  ggplot2::theme_bw()


###################################################
### code chunk number 43: road-casualties-in-great-britain-from-1969-to-1984-ar-plot
###################################################
result_ar <- fastcpd.ar(
  data = diff(uk_seatbelts[, "drivers"], lag = 12), order = 1,
  beta = "BIC", cost_adjustment = "BIC"
)
plot(result_ar)


###################################################
### code chunk number 44: road-casualties-in-great-britain-from-1969-to-1984-lm-setup
###################################################
result_lm <- fastcpd.lm(
  diff(uk_seatbelts[, c("drivers", "kms", "PetrolPrice", "law")], 12)
)


###################################################
### code chunk number 45: road-casualties-in-great-britain-from-1969-to-1984-lm-plot
###################################################
cp_dates <- as.Date("1969-01-01", format = "%Y-%m-%d")
cp_dates <- cp_dates + lubridate::period(month = 1 + result_lm@cp_set + 12)
cp_dates <- zoo::as.yearmon(cp_dates)

dates <- zoo::as.yearmon(time(uk_seatbelts))

uk_seatbelts_df <- data.frame(
  dates = dates,
  drivers = c(uk_seatbelts[, "drivers"]),
  color = as.factor((dates < cp_dates[1]) + (dates < cp_dates[2]))
)

ggplot2::ggplot() +
  ggplot2::geom_line(
    data = uk_seatbelts_df,
    mapping = ggplot2::aes(x = dates, y = drivers, color = color)
  ) +
  ggplot2::geom_vline(
    xintercept = cp_dates,
    linetype = "dashed",
    color = "red"
  ) +
  zoo::scale_x_yearmon() +
  ggplot2::annotate(
    "text",
    x = cp_dates,
    y = 1025,
    label = as.character(cp_dates),
    color = "blue"
  ) +
  ggplot2::theme_bw() +
  ggplot2::theme(legend.position = "none")


###################################################
### code chunk number 46: vanilla-pelt-setup
###################################################
set.seed(1)
n <- 300
p_true <- 4
p <- 40
x <- mvtnorm::rmvnorm(n, rep(0, p), diag(p))
theta_0 <- rbind(
  runif(p_true, 1, 4),
  runif(p_true, -5, 5),
  runif(p_true, 1, 4),
  runif(p_true, -3, 3)
)
theta_0 <- cbind(theta_0, matrix(0, ncol = p - p_true, nrow = 4))
y <- c(
  x[1:(n * 0.25), ] %*% theta_0[1, ] + rnorm(n * 0.25),
  x[(n * 0.25 + 1):(n * 0.5), ] %*% theta_0[2, ] + rnorm(n * 0.25),
  x[(n * 0.5 + 1):(n * 0.75), ] %*% theta_0[3, ] + rnorm(n * 0.25),
  x[(n * 0.75 + 1):n, ] %*% theta_0[4, ] + rnorm(n * 0.25)
)
small_lasso <- cbind.data.frame(y, x)


###################################################
### code chunk number 47: vanilla-pelt-0-result1
###################################################
fastcpd.lasso(small_lasso, segment_count = 2)@cp_set


###################################################
### code chunk number 48: vanilla-pelt-0-result2
###################################################
fastcpd.lasso(
  small_lasso,
  segment_count = 2, vanilla_percentage = 0.5
)@cp_set


###################################################
### code chunk number 49: vanilla-pelt-0-result-plot
###################################################
microbenchmark:::autoplot.microbenchmark(microbenchmark::microbenchmark(
  "vanilla_percentage = 0" = fastcpd::fastcpd.lasso(small_lasso, segment_count = 2),
  "vanilla_percentage = 0.5" = fastcpd::fastcpd.lasso(small_lasso, segment_count = 2, vanilla_percentage = 0.5),
  "vanilla_percentage = 1" = fastcpd::fastcpd.lasso(small_lasso, segment_count = 2, vanilla_percentage = 1),
  times = 10
), log = FALSE) + ggplot2::theme_bw()


###################################################
### code chunk number 50: adaptive-number-of-epochs
###################################################
multiple_epochs <- function(segment_length) {
  if (segment_length < 50) 1 else 0
}


###################################################
### code chunk number 51: adaptive-number-of-epochs-result
###################################################
fastcpd.lasso(
  small_lasso,
  segment_count = 2,
  multiple_epochs = multiple_epochs
)@cp_set