## -----------------------------------------------------------------------------
library(MassSpecWavelet)

## -----------------------------------------------------------------------------
onetwothree <- c(1,2,3)
x <- c(0, onetwothree, onetwothree, 0, 0, 0, onetwothree, 1, onetwothree, onetwothree, 0)
plot(x, type = "o", main = "Five peaks are expected")

## -----------------------------------------------------------------------------
options("MassSpecWavelet.localMaximum.algorithm" = "new")
local_max <- which(localMaximum(x, winSize = 4) > 0)
plot(x, type = "o", main = "With the new algorithm, 5/5 peaks are found")
points(local_max, x[local_max], col = "red", pch = 20)


## -----------------------------------------------------------------------------
options("MassSpecWavelet.localMaximum.algorithm" = "classic")
local_max <- which(localMaximum(x, winSize = 4) > 0)
plot(x, type = "o", main = "With the classic algorithm, 3/5 peaks are found")
points(local_max, x[local_max], col = "red", pch = 20)

## -----------------------------------------------------------------------------
x <- c(0, 1, 2, 3, 3, 3, 3, 2, 1, 0, 3, 0, 3, 3, 3, 3, 3, 0)
x <- c(x, 0, 1, 2, 3, 3, 3, 2, 1, 0, 3, 0, 0, 0, 3, 3, 3, 0, 0)
options("MassSpecWavelet.localMaximum.algorithm" = "classic")
local_max_classic <- which(localMaximum(x, winSize = 5) > 0)
options("MassSpecWavelet.localMaximum.algorithm" = "new")
local_max_new <- which(localMaximum(x, winSize = 5) > 0)
par(mfrow = c(2, 1))
plot(x, type = "o", main = "With the classic algorithm, 2/6 peaks are found")
points(local_max_classic, x[local_max_classic], col = "red", pch = 20)
plot(x, type = "o", main = "With the new algorithm, 6/6 peaks are found")
points(local_max_new, x[local_max_new], col = "blue", pch = 20)


## ----benchmark, fig.cap="CPU time vs length, for several algorithms and windows"----
# Run this interactively
set.seed(5413L)
winSizes <- c(5, 31, 301)
xlengths <- c(20, 200, 2000, 20000, 200000)
out <- vector("list", length(winSizes) * length(xlengths))
i <- 0L
for (winSize in winSizes) {
    for (xlength in xlengths) {
        i <- i + 1L
        x <- round(10*runif(xlength), 1)*10
        bm1 <- as.data.frame(
                bench::mark(
                classic = {
                    options(MassSpecWavelet.localMaximum.algorithm = "classic")
                    localMaximum(x, winSize = winSize)
                    
                },
                faster = {
                    options(MassSpecWavelet.localMaximum.algorithm = "faster")
                    localMaximum(x, winSize = winSize)
                    
                },
                check = TRUE,
                filter_gc = FALSE,
                time_unit = "ms"
            )
        )
        bm2 <- as.data.frame(
                bench::mark(
                new = {
                    options(MassSpecWavelet.localMaximum.algorithm = "new")
                    localMaximum(x, winSize = winSize)
                    
                },
                check = FALSE,
                filter_gc = FALSE,
                time_unit = "ms"
            )
        )
        out[[i]] <- data.frame(
            algorithm = c(as.character(bm1$expression), as.character(bm2$expression)),
            min_cpu_time_ms = c(as.numeric(bm1$min), as.numeric(bm2$min)),
            xlength = xlength,
            winSize = winSize
        )
    }
}
out2 <- do.call(rbind, out)
plot(
    x = out2$xlength,
    y = out2$min_cpu_time_ms,
    col = ifelse(out2$algorithm == "new", "red", ifelse(out2$algorithm == "faster", "darkgreen", "blue")),
    pch = ifelse(out2$winSize == 5, 1, ifelse(out2$winSize == 31, 2, 3)),
    log = "xy",
    xlab = "Signal length",
    ylab = "Min CPU time (ms)"
)
legend(
    "topleft",
    legend = c("New algorithm", "Classic algorithm", "Faster algorithm",
               "winSize = 5", "winSize = 31", "winSize = 301"),
    lty = c(1,1,1, 0, 0, 0),
    pch = c(NA, NA, NA,1, 2, 3),
    col = c("red", "blue", "darkgreen", "black", "black", "black")
)

## -----------------------------------------------------------------------------
data(exampleMS)
scales <- seq(1, 64, 1)
wCoefs <- cwt(exampleMS, scales = scales, wavelet = "mexh")

## -----------------------------------------------------------------------------
plotRange <- c(8000, 9000)

## -----------------------------------------------------------------------------
plot(exampleMS, xlim = plotRange, type = "l")

## -----------------------------------------------------------------------------
matplot(
    wCoefs[,ncol(wCoefs):1], 
    type = "l",
    col = rev(rainbow(64, start = 0.7, end = 0.1, alpha = 0.5)[scales]),
    lty = 1,
    xlim = c(8000, 9000),
    xlab = "m/z index",
    ylab = "CWT coefficients"
)
legend(
    x = "topright",
    legend = sprintf("scales = %d", scales[seq(1, length(scales), length.out = 4)]),
    lty = 1,
    col = rainbow(64, start = 0.7, end = 0.1)[scales[seq(1, length(scales), length.out = 4)]]
)


## -----------------------------------------------------------------------------
options(MassSpecWavelet.localMaximum.algorithm = "new")
localMax_new <- getLocalMaximumCWT(wCoefs)
options(MassSpecWavelet.localMaximum.algorithm = "classic")
localMax_classic <- getLocalMaximumCWT(wCoefs)

## ----localMax, width=10, height=8,fig.align='center',fig.cap="Identified local maxima of CWT coefficients at each scale"----
par(mfrow = c(2,1))
plotLocalMax(localMax_classic, NULL, range=c(plotRange[1], plotRange[2]), main = "classic")
plotLocalMax(localMax_new, NULL, range=c(plotRange[1], plotRange[2]), main = "new")

## -----------------------------------------------------------------------------
sessionInfo()