tsfeatureExtend: Time Series Features Calculations

#' @title mean
#'
#' @export
tsf_mean <- function(x, par.vals = list()) {
  f <- purrr::safely(.f = function(x) {
    mean(x, na.rm = T)
  })
  c(mean = f(x)$result)
}


#' @title Variance
#'
#' @export
tsf_var <- function(x, par.vals = list()) {
  f <- purrr::safely(.f = function(x) {
    var(x, na.rm = T)
  })
  c(var = f(x)$result)
}


#' @title Minimum
#'
#' @export
tsf_min <- function(x, par.vals = list()) {
  f <- purrr::safely(.f = function(x) {
    min(x, na.rm = T)
  })
  c(min = f(x)$result)
}


#' @title Maximum
#'
#' @export
tsf_max <- function(x, par.vals = list()) {
  f <- purrr::safely(.f = function(x) {
    max(x, na.rm = T)
  })
  c(max = f(x)$result)
}


#' @title Burst
#'
#' @export
tsf_rmeaniqmean <- function(x, par.vals = list()) {
  f <- purrr::safely(.f = function(x) {
    mean(x, trim = 0.5, na.rm = T) / mean(x, na.rm = T)
  })
  c(rmeaniqmean = f(x)$result)
}


#' @title Moments
#'
#' A function to calculate the moments Returns the moment of the distribution (the measure moments, m of the distribution,
#' for m = 3,4,5,...,11.) of the input time series, normalizes by the standard deviation. Output This operation Uses the
#' moments package in R measure29 - moments m of the distribution
#'
#' @export
tsf_moments <- function(x, par.vals = list()) {
  f <- purrr::safely(.f = function(x) {
    momentssd <- 0
    momentssd <- moments::moment(x, order = 3, na.rm = TRUE)
    momentssd / sd(x, na.rm = TRUE)
  })
  c(moment = f(x)$result)
}


#' @title Ratio of val
#'
#' A function to compare the means of data that is below and upper the global mean Calculates a statistic related to the
#' mean of the time series data that is above the (global) time-series mean compared to the mean of the data that is below
#' the global time-series mean. measure32 - compare the means of data that is below and upper the global mean
#'
#' @export
tsf_highlowmu <- function(x, par.vals = list()) {
  f <- purrr::safely(.f = function(x) {
    mu <- mean(x, na.rm = T)
    mhi <- mean(x[x > mu], na.rm = T)
    mlo <- mean(x[x < mu], na.rm = T)
    mhi - mu / (mu - mlo)
  })
  c(highlowmu = f(x)$result)
}


#' @title Burst
#'
#' @export
tsf_burst <- function(x, par.vals = list()) {
  f <- purrr::safely(.f = function(x) {
    (sd(x, na.rm = T)^2) / mean(x, na.rm = T)
  })
  c(burst = f(x)$result)
}


#' @title relative first index of max
#'
#' @export
tsf_first_index_max <- function(x, par.vals = list()) {
  f <- purrr::safely(.f = function(x) {
    which.max(x) / length(x)
  })
  c(first_index_max = f(x)$result)
}


#' @title relative first index of min
#'
#' @export
tsf_first_index_min <- function(x, par.vals = list()) {
  f <- purrr::safely(.f = function(x) {
    which.min(x) / length(x)
  })
  c(first_index_min = f(x)$result)
}


#' @title relative last index of max
#'
#' @export
tsf_last_index_max <- function(x, par.vals = list()) {
  f <- purrr::safely(.f = function(x) {
    n <- length(x)
    (n - which.max(rev(x)) + 1) / n
  })
  c(last_index_max = f(x)$result)
}


#' @title relative last index of min
#'
#' @export
tsf_last_index_min <- function(x, par.vals = list()) {
  f <- purrr::safely(.f = function(x) {
    n <- length(x)
    (n - which.min(rev(x)) + 1) / n
  })
  c(last_index_min = f(x)$result)
}


#' @title longest strike above mean
#'
#' @export
tsf_longest_strike_above_mean <- function(x, par.vals = list()) {
  f <- purrr::safely(.f = function(x) {
    seqx <- as.numeric(x > mean(x, na.rm=T))
    rlex <- rle(seqx)
    max(rlex$lengths[rlex$values == 1])
  })
  c(longest_strike_above_mean = f(x)$result)
}


#' @title longest strike below mean
#'
#' @export
tsf_longest_strike_below_mean <- function(x, par.vals = list()) {
  f <- purrr::safely(.f = function(x) {
    seqx <- as.numeric(x < mean(x, na.rm=T))
    rlex <- rle(seqx)
    max(rlex$lengths[rlex$values == 1])
  })
  c(longest_strike_below_mean = f(x)$result)
}


#' @title longest strike above median
#'
#' @export
tsf_longest_strike_above_median <- function(x, par.vals = list()) {
  f <- purrr::safely(.f = function(x) {
    seqx <- as.numeric(x > median(x, na.rm=T))
    rlex <- rle(seqx)
    max(rlex$lengths[rlex$values == 1])
  })
  c(longest_strike_above_median = f(x)$result)
}


#' @title longest strike below median
#'
#' @export
tsf_longest_strike_below_median <- function(x, par.vals = list()) {
  f <- purrr::safely(.f = function(x) {
    seqx <- as.numeric(x < median(x, na.rm=T))
    rlex <- rle(seqx)
    max(rlex$lengths[rlex$values == 1])
  })
  c(longest_strike_below_median = f(x)$result)
}


#' @title longest strike negative
#'
#' @export
tsf_longest_strike_negative <- function(x, par.vals = list()) {
  f <- purrr::safely(.f = function(x) {
    seqx <- as.numeric(x < 0)
    rlex <- rle(seqx)
    rlex <- rlex$lengths[rlex$values == 1]
    if (length(rlex)) {
      return(max(rlex))
    } else {
      return(0)
    }
  })
  c(longest_strike_negative = f(x)$result)
}


#' @title longest strike positive
#'
#' @export
tsf_longest_strike_positive <- function(x, par.vals = list()) {
  f <- purrr::safely(.f = function(x) {
    seqx <- as.numeric(x > 0)
    rlex <- rle(seqx)
    rlex <- rlex$lengths[rlex$values == 1]
    if (length(rlex)) {
      return(max(rlex))
    } else {
      return(0)
    }
  })
  c(longest_strike_positive = f(x)$result)
}


#' @title longest strike zero
#'
#' @export
tsf_longest_strike_zero <- function(x, par.vals = list()) {
  f <- purrr::safely(.f = function(x) {
    seqx <- as.numeric(x == 0)
    rlex <- rle(seqx)
    rlex <- rlex$lengths[rlex$values == 1]
    if (length(rlex)) {
      return(max(rlex))
    } else {
      return(0)
    }
  })
  c(longest_strike_zero = f(x)$result)
}


#' @title mean absolute change
#'
#' @export
tsf_mean_absolute_change <- function(x, par.vals = list()) {
  f <- purrr::safely(.f = function(x) {
    mean(abs(diff(x)), na.rm = T)
  })
  c(mean_absolute_change = f(x)$result)
}


#' @title mean change
#'
#' @export
tsf_mean_change <- function(x, par.vals = list()) {
  f <- purrr::safely(.f = function(x) {
    mean(diff(x), na.rm = T)
  })
  c(mean_change = f(x)$result)
}



#' @title mean second derivative central
#'
#' @export
tsf_mean_second_derivative_central <- function(x, par.vals = list()) {
  f <- purrr::safely(.f = function(x) {
    tibble::tibble(x = x) %>%
      dplyr::mutate(xlag = dplyr::lag(x, 1),
                    xlead = dplyr::lead(x, 1)) %>%
      dplyr::mutate(v = 0.5 * (xlag - 2 * x + xlead)) %>%
      na.omit() %>%
      .[["v"]] %>%
      mean()
  })
  c(mean_second_derivative_central = f(x)$result)
}


#' @title time reversal asymmetric statistics
#'
#' @export
tsf_time_reversal_asymmetric_statistics <- function(x, par.vals = list(l = 1)) {
  f <- purrr::safely(.f = function(x, l = 1) {
    tibble::tibble(x = x) %>%
      dplyr::mutate(xlead2 = dplyr::lead(x, 2 * l),
                    xlead = dplyr::lead(x, l)) %>%
      dplyr::mutate(v = 0.5 * (xlead^2 * xlead - xlead * x^2)) %>%
      na.omit() %>%
      .[["v"]] %>%
      mean()
  })
  c(time_reversal_asymmetric_statistics = f(x, l = par.vals[["l"]])$result)
}


#' @title mean autocorrelation
#'
#' @export
tsf_mean_autocorrelation <- function(x, par.vals = list()) {
  f <- purrr::safely(.f = function(x) {
    acfAll <- acf(x, plot = F, lag.max = length(x) - 1, demean = T)
    mean(acfAll$acf[-1]) / (sd(x, na.rm=T))^2
  })
  c(mean_autocorrelation = f(x)$result)
}


#' @title Crossing Points
#'
#' @export
tsf_crossing_points <- function(x, par.vals = list()) {
  f <- purrr::safely(.f = function(x) {
    midline <- median(x, na.rm = TRUE)
    ab <- x <= midline
    lenx <- length(x)
    p1 <- ab[1:(lenx - 1)]
    p2 <- ab[2:lenx]
    cross <- (p1 & !p2) | (p2 & !p1)
    sum(cross, na.rm = T)
  })
  c(crossing_points = f(x)$result)
}


#' @title Lumpiness
#'
#' @export
tsf_lumpiness <- function(x, par.vals = list(width = 10)) {
  f <- purrr::safely(.f = function(x, width) {
    x <- scalets(x)
    nr <- length(x)
    lo <- seq(1, nr, by = width)
    up <- seq(width, nr + width, by = width)
    nsegs <- nr/width
    varx <- purrr::map_dbl(seq_len(nsegs), function(idx) var(x[lo[idx]:up[idx]],
                                                             na.rm = TRUE))
    if (length(x) < 2 * width)
      lumpiness <- 0
    else lumpiness <- var(varx, na.rm = TRUE)
    lumpiness
  }, otherwise = NA_real_)
  c(lumpiness = f(x, par.vals[["width"]])$result)
}


#' @title Stability
#'
#' @export
tsf_stability <- function(x, par.vals = list(width = 10)) {
  f <- purrr::safely(.f = function(x, width) {
    x <- scalets(x)
    nr <- length(x)
    lo <- seq(1, nr, by = width)
    up <- seq(width, nr + width, by = width)
    nsegs <- nr / width
    meanx <- purrr::map_dbl(seq_len(nsegs), function(idx) mean(x[lo[idx]:up[idx]], na.rm = TRUE))
    if (length(x) < 2 * width)
      stability <- 0
    else stability <- var(meanx, na.rm = TRUE)
    stability
  }, otherwise = NA_real_)
  c(stability = f(x, par.vals[["width"]])$result)
}


#' @title Maximal Variance Shift
#'
#' @export
tsf_max_var_shift <- function(x, par.vals = list(width = 10)) {
  f <- purrr::safely(.f = function(x, width) {
    rollvar <- RcppRoll::roll_var(x, width, na.rm = TRUE)
    vars <- abs(diff(rollvar, width))
    if (length(vars) == 0L) {
      maxVar <- 0
    } else if (all(is.na(vars))) {
      maxVar <- NA_real_
    } else {
      maxVar <- max(vars, na.rm = TRUE)
    }
    maxVar
  }, otherwise = NA_real_)
  c(max_var_shift = f(x, par.vals[["width"]])$result)
}


#' @title Maximal Level shift
#'
#' @export
tsf_max_level_shift <- function(x, par.vals = list(width = 10)) {
  f <- purrr::safely(.f = function(x, width) {
    rollmean <- RcppRoll::roll_mean(x, width, na.rm = TRUE)
    means <- abs(diff(rollmean, width))
    if (length(means) == 0L) {
      maxmeans <- 0
    }
    else if (all(is.na(means))) {
      maxmeans <- NA_real_
    }
    else {
      maxmeans <- max(means, na.rm = TRUE)
    }
    maxmeans
  }, otherwise = NA_real_)
  c(max_level_shift = f(x, par.vals[["width"]])$result)
}


#' @title Maximal KL Shift
#'
#' @export
tsf_max_kl_shift <- function(x, par.vals = list(width = 10)) {
  f <- purrr::safely(.f = function(x, width) {
    gw <- 100
    xgrid <- seq(min(x, na.rm = TRUE), max(x, na.rm = TRUE),
                 length = gw)
    grid <- xgrid[2L] - xgrid[1L]
    tmpx <- x[!is.na(x)]
    bw <- bw.nrd0(tmpx)
    lenx <- length(x)
    if (lenx <= (2 * width)) {
      return(c(max_kl_shift = NA_real_, time_kl_shift = NA_real_))
    }
    dens.mat <- matrix(NA, nrow = lenx, ncol = gw)
    for (i in 1L:lenx) {
      dens.mat[i, ] <- dnorm(xgrid, mean = x[i], sd = bw)
    }
    dens.mat <- pmax(dens.mat, dnorm(38))
    rmean <- RcppRoll::roll_mean(dens.mat, n = width, na.rm = TRUE,
                                 fill = NA, align = "right")
    lo <- seq(1, lenx - width + 1)
    hi <- seq(width + 1, lenx)
    seqidx <- min(length(lo), length(hi))
    kl <- sapply(1:seqidx, function(i) sum(rmean[lo[i], ] * (log(rmean[lo[i],]) - log(rmean[hi[i], ])) * grid, na.rm = TRUE))
    diffkl <- diff(kl, na.rm = TRUE)
    if (length(diffkl) == 0L) {
      diffkl <- 0
    }
    max(diffkl)
  }, otherwise = NA_real_)
  c(max_kl_shift = f(x, par.vals[["width"]])$result)
}


#' @title Arch Stat
#'
#' @export
tsf_arch_stat <- function(x, par.vals = list(lags = 12, demean = T)) {
  f <- purrr::safely(.f = function(x, lags = 12, demean = T) {
    if (length(x) <= 13) {
      S <- NA_real_
    }
    if (demean) {
      x <- x - mean(x, na.rm = TRUE)
    }

    mat <- embed(x^2, lags + 1)
    fit <- lm(mat[, 1] ~ mat[, -1])
    arch.lm <- summary(fit)
    S <- arch.lm$r.squared
    S
  }, otherwise = NA_real_)
  c(arch_stat = f(x, par.vals[["lags"]], par.vals[["demean"]])$result)
}


#' flat Spots
#'
#' @export
tsf_flat_spots <- function(x, par.vals = list()) {
  f <- purrr::safely(.f = function(x) {
    cutx <- cut(x, breaks = 10, include.lowest = TRUE, labels = FALSE)
    rlex <- rle(cutx)
    max(rlex$lengths)
  }, otherwise = NA_real_)
  c(flat_spots = f(x)$result)
}


#' Heterogenity
#'
#'
#' @export
tsf_herterogenity <- function(x, par.vals = list()) {
  f <- purrr::safely(.f = function(x, lags = 12, demean = T) {
    x.whitened <- na.contiguous(ar(x)$resid)
    x.archtest <- arch_stat(x.whitened)
    LBstat <- sum(acf(x.whitened^2, lag.max = 12L, plot = FALSE)$acf[-1L]^2)
    garch.fit <- suppressWarnings(tseries::garch(x.whitened, trace = FALSE))
    garch.fit.std <- residuals(garch.fit)
    x.garch.archtest <- arch_stat(garch.fit.std)
    LBstat2 <- sum(acf(na.contiguous(garch.fit.std^2), lag.max = 12L, plot = FALSE)$acf[-1L]^2)
    output <- c(arch_acf = LBstat, garch_acf = LBstat2, arch_r2 = unname(x.archtest),
                garch_r2 = unname(x.garch.archtest))
    output[is.na(output)] <- 1
    output[1]
  }, otherwise = NA_real_)
  c(herterogenity = f(x)$result)
}


#' @title Ratio beyond r time sigma
#'
#' @export
tsf_ratio_beyond_r_sigma <- function(x, par.vals = list(r = 1.5)) {
  f <- purrr::safely(.f = function(x, r) {
    sum(abs(x - mean(x, na.rm = T)) > r * sd(x, na.rm = T), na.rm = T) / length(x)
  }, otherwise = NA_real_)
  c(ratio_beyond_r_sigma = f(x, par.vals[["r"]])$result)
}


#' @title Large Standard Deviation
#'
#' @export
tsf_large_sd <- function(x, par.vals = list(r = 1.5)) {
  f <- purrr::safely(.f = function(x, r) {
    as.numeric(sd(x, na.rm = T) > r * (max(x, na.rm = T) - min(x, na.rm = T)))
  }, otherwise = NA_real_)
  c(large_sd = f(x, par.vals[["r"]])$result)
}


#' @title Has duplicated Maximum
#'
#' @export
tsf_has_duplicated_max <- function(x, par.vals = list()) {
  f <- purrr::safely(.f = function(x) {
    m <- max(x, na.rm = T)
    as.numeric(sum(x == m, na.rm = T))
  }, otherwise = NA_real_)
  c(has_duplicated_max = f(x)$result)
}


#' @title Has Duplicated Minimum
#'
#' @export
tsf_has_duplicated_min <- function(x, par.vals = list()) {
  f <- purrr::safely(.f = function(x, r) {
    m <- min(x, na.rm = T)
    as.numeric(sum(x == m, na.rm = T))
  }, otherwise = NA_real_)
  c(has_duplicated_min = f(x)$result)
}


#' @title Has Duplicated
#'
#' @export
tsf_has_duplicate <- function(x, par.vals = list()) {
  f <- purrr::safely(.f = function(x) {
    as.numeric(any(duplicated(x)))
  }, otherwise = NA_real_)
  c(has_duplicate = f(x)$result)
}


#' @title Sum of Series
#'
#' @export
tsf_sum_value <- function(x, par.vals = list()) {
  f <- purrr::safely(.f = function(x) {
    sum(x, na.rm = T)
  }, otherwise = NA_real_)
  c(sum_value = f(x)$result)
}


#' @title Time Series Complexity
#'
#' @export
tsf_cid_ce <- function(x, par.vals = list(normalize = T)) {
  f <- purrr::safely(.f = function(x, normalize) {
    if (normalize) {
      x <- (x - mean(x, na.rm = T)) / sd(x, na.rm = T)
    }
    tibble::tibble(x = x) %>%
      dplyr::mutate(xl = dplyr::lead(x, 1)) %>%
      dplyr::mutate(r = (x - xl)^2) %>%
      na.omit() %>%
      .[["r"]] %>%
      mean()
  }, otherwise = NA_real_)
  c(cid_ce = f(x, par.vals[["normalize"]])$result)
}


#' @title Returns the sum of all values, that are present in the time series
#' more than once.
#'
#' @export
tsf_sum_of_reoccurring_values <- function(x, par.vals = list()) {
  f <- purrr::safely(.f = function(x, normalize) {
    dx <- data.table::data.table(x = x)
    sum(dx[, .(n = .N), by = x][n == 1, n := 0][n > 1, n := 1][["n"]], na.rm = T)
  }, otherwise = NA_real_)
  c(sum_of_reoccurring_values = f(x, normalize)$result)
}


#' @title Returns the sum of all data points, that are present in the time
#' series more than once.
#'
#' @export
tsf_sum_of_reoccurring_data_points <- function(x, par.vals = list()) {
  f <- purrr::safely(.f = function(x, normalize) {
    dx <- data.table::data.table(x = x)
    sum(dx[, .(n = .N), by = x][n == 1, n := 0][["n"]], na.rm = T)
  }, otherwise = NA_real_)
  c(sum_of_reoccurring_data_points = f(x, normalize)$result)
}


#' @title Returns a factor which is 1 if all values in the time series occur
#' only once, and below one if this is not the case.
#'
#' @export
tsf_ratio_value_number_to_time_series_length <- function(x, par.vals = list()) {
  f <- purrr::safely(.f = function(x, normalize) {
    x %>%
      dplyr::n_distinct() -> nx
    nx / length(x)
  }, otherwise = NA_real_)
  c(ratio_value_number_to_time_series_length = f(x, normalize)$result)
}


#' @title Count observed values within the interval [min, max).
#'
#' @export
tsf_range_count <- function(x, par.vals = list()) {
  f <- purrr::safely(.f = function(x, normalize) {
    sum(x >= min(x, na.rm = T) & x < max(x, na.rm = T), na.rm = T)
  }, otherwise = NA_real_)
  c(range_count = f(x, normalize)$result)
}


#' @title absolute energy
#'
#' @export
tsf_absolute_energy <- function(x, par.vals = list()) {
  f <- purrr::safely(.f = function(x) {
    sum(x^2, na.rm = T)
  }, otherwise = NA_real_)
  c(absolute_energy = f(x)$result)
}


#' @title spectral entropy
#'
#' @export
tsf_entropy <- function(x, par.vals = list()) {
  f <- purrr::safely(.f = function(x) {
    ForeCA::spectral_entropy(na.contiguous(x))[1L]
  }, otherwise = NA_real_)
  c(entropy = f(x)$result)
}


#' @title Binned Entropy
#'
#' @description This feature calculator bins the values of the time series  into m
#' equidistant bins.
#'
#' @export
tsf_binned_entropy <- function(x, par.vals = list(m = 10)) {
  f <- purrr::safely(.f = function(x, m = 10) {
    cutx <- data.table::data.table(bins = cut(x, breaks = m, include.lowest = T, labels = F))
    entropy <- cutx[, .(n = .N), by = bins][, n := n / sum(n, na.rm=T)][n > 0, p := n * log(n)]
    sum(entropy$p, na.rm = T)
  })
  c(binned_entropy = f(x, par.vals[["m"]])$result)
}


#' @title Has Large Std Dev
#'
#' @export
tsf_has_large_standard_deviation <- function(x, par.vals = list()) {
  f <- purrr::safely(.f = function(x) {
    std <- sd(x, na.rm = T)
    as.numeric(std > 0.5 * (max(x, na.rm = T) - min(x, na.rm = T)))
  })
  c(has_large_standard_deviation = f(x)$result)
}


#' @title has larger variance than std dev
#'
#' @export
tsf_has_larger_var_than_std <- function(x, par.vals = list()) {
  f <- purrr::safely(.f = function(x) {
    std <- sd(x, na.rm = T)
    as.numeric(std^2 > std)
  })
  c(has_larger_var_than_std = f(x)$result)
}


#' @title is symmetric looking
#'
#' @export
tsf_is_symmetric_looking <- function(x, par.vals = list()) {
  f <- purrr::safely(.f = function(x) {
    as.numeric(abs(median(x, na.rm=T) - mean(x, na.rm = T)) < 0.5 * (max(x, na.rm = T) - min(x, na.rm = T)))
  })
  c(is_symmetric_looking = f(x)$result)
}


#' @title Elements Above Mean
#'
#' @export
tsf_n_above_mean <- function(x, par.vals = list()) {
  f <- purrr::safely(.f = function(x) {
    m <- mean(x, na.rm = T)
    sum(x > m, na.rm = T)
  })
  c(n_above_mean = f(x)$result)
}


#' @title Elements Below Mean
#'
#' @export
tsf_n_below_mean <- function(x, par.vals = list()) {
  f <- purrr::safely(.f = function(x) {
    m <- mean(x, na.rm = T)
    sum(x < m, na.rm = T)
  })
  c(n_below_mean = f(x)$result)
}


#' @title Elements Above Mean
#'
#' @export
tsf_n_above_median <- function(x, par.vals = list()) {
  f <- purrr::safely(.f = function(x) {
    m <- median(x, na.rm = T)
    sum(x > m, na.rm = T)
  })
  c(n_above_median = f(x)$result)
}


#' @title Elements Below Median
#'
#' @export
tsf_n_below_median <- function(x, par.vals = list()) {
  f <- purrr::safely(.f = function(x) {
    m <- median(x, na.rm = T)
    sum(x < m, na.rm = T)
  })
  c(n_below_median = n_below_medianf(x)$result)
}


#' @title Hurst Coefficient
#'
#' @description \url{https://pdfs.semanticscholar.org/a4e5/f5f6ef9db520cf2c5d355eceb6b439a0a852.pdf}
#'
#' @export
tsf_hurst_coefficient <- function(x, par.vals = list()) {
  f <- purrr::safely(.f = function(x) {
    fracdiff::fracdiff(na.contiguous(x), 0, 0)[["d"]] + 0.5
  })
  c(hurst_coefficient = f(x)$result)
}


#' @title Trend
#'
#' @export
tsf_trend <- function(x, par.vals = list()) {
  f <- purrr::safely(.f = function(x) {
    ts_info(x)$trend
  })
  c(trend = f(x)$result)
}


#' @title Spike
#'
#' @export
tsf_spike <- function(x, par.vals = list()) {
  f <- purrr::safely(.f = function(x) {
    ts_info(x)$spike
  })
  c(spike = f(x)$result)
}


#' @title Linearity
#'
#' @export
tsf_linearity <- function(x, par.vals = list()) {
  f <- purrr::safely(.f = function(x) {
    ts_info(x)$linearity
  })
  c(linearity = f(x)$result)
}


#' @title Curvature
#'
#' @export
tsf_curvature <- function(x, par.vals = list()) {
  f <- purrr::safely(.f = function(x) {
    ts_info(x)$curvature
  })
  c(curvature = f(x)$result)
}


#' Timeseries calculations
#' borrowed from tsfeatures
#'
#' @export
ts_info <- function(x) {
  if (!is(x, 'ts')) {
    x <- as.ts(x)
  }
  tspx <- tsp(x)
  freq <- tspx[3]
  contx <- try(na.contiguous(x), silent = TRUE)
  len.contx <- length(contx)
  if (length(contx) < 2 * freq || class(contx) == "try-error") {
    trend <- linearity <- curvature <- season <- spike <- peak <- trough <- NA
  } else {
    if (freq > 1L) {
      all.stl <- stl(contx, s.window = "periodic", robust = TRUE)
      starty <- start(contx)[2L]
      pk <- (starty + which.max(all.stl$time.series[, "seasonal"]) - 1L)%%freq
      th <- (starty + which.min(all.stl$time.series[, "seasonal"]) - 1L)%%freq
      pk <- ifelse(pk == 0, freq, pk)
      th <- ifelse(th == 0, freq, th)
      trend0 <- all.stl$time.series[, "trend"]
      fits <- trend0 + all.stl$time.series[, "seasonal"]
      adj.x <- contx - fits
      v.adj <- var(adj.x, na.rm = TRUE)
      detrend <- contx - trend0
      deseason <- contx - all.stl$time.series[, "seasonal"]
      peak <- pk * max(all.stl$time.series[, "seasonal"], na.rm = TRUE)
      trough <- th * min(all.stl$time.series[, "seasonal"], na.rm = TRUE)
      remainder <- all.stl$time.series[, "remainder"]
      season <- ifelse(var(detrend, na.rm = TRUE) < 1e-10, 0, max(0, min(1, 1 - v.adj/var(detrend, na.rm = TRUE))))
    } else {
      # No seasonal component
      tt <- 1:len.contx
      trend0 <- fitted(mgcv::gam(contx ~ s(tt)))
      remainder <- contx - trend0
      deseason <- contx - trend0
      v.adj <- var(trend0, na.rm = TRUE)
    }
    trend <- ifelse(var(deseason, na.rm = TRUE) < 1e-10, 0, max(0, min(1, 1 - v.adj/var(deseason, na.rm = TRUE))))
    n <- length(remainder)
    v <- var(remainder, na.rm = TRUE)
    d <- (remainder - mean(remainder, na.rm = TRUE))^2
    varloo <- (v * (n - 1) - d)/(n - 2)
    spike <- var(varloo, na.rm = TRUE)
    pl <- poly(1:len.contx, degree = 2)
    tren.coef <- coef(lm(trend0 ~ pl))[2:3]
    linearity <- tren.coef[1]
    curvature <- tren.coef[2]
  }
  if (freq > 1) {
    return(list(trend = trend, season = season, spike = spike, peak = peak, trough = trough, linearity = linearity, curvature = curvature))
  } else {
    # No seasonal component
    return(list(trend = trend, spike = spike, linearity = linearity, curvature = curvature))
  }
}