Example: fast convolution

library(insitu)

#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
# Vectorised solution
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
convv <- function(x, y) {
  nx <- length(x)
  ny <- length(y)
  xy <- rbind(outer(x,y),
              matrix(0, nx, ny))
  nxy <- nx + ny - 1
  length(xy) <- nxy * ny
  dim(xy) <- c(nxy, ny)
  rowSums(xy)
}

#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
# Nested for loops - initial version
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
conv0 <- function(x,y) {
  nx <- length(x)
  ny <- length(y)
  z <- numeric(nx + ny - 1)
  for(i in seq(length = nx)) {
    xi <- x[[i]]
    for(j in seq(length = ny)) {
      ij <- i + j - 1
      z[[ij]] <- z[[ij]] + xi * y[[j]]
    }
  }
  z
}



#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
# Base R: using 'replace()'
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
conv2 <- function(x,y) {
  nx <- length(x)
  ny <- length(y)
  sy <- seq_along(y) - 1L
  z <- numeric(nx + ny - 1)
  for(i in seq(length = nx)) {
    ij <- i + sy
    z <- replace(z, ij, z[ij] + x[[i]] * y)
  }
  z
}


#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
# by-reference 3 - using 'br_copy()'
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
conv3 <- function(x,y) {
  nx <- length(x)
  ny <- length(y)
  sy <- seq_along(y) - 1L
  z <- numeric(nx + ny - 1)
  for(i in seq(length = nx)) {
    ij <- i + sy
    br_copy(z, z[ij] + x[[i]] * y, n = ny, xi = i)
  }
  z
}


#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
# by-reference 4
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
conv4 <- function(x,y) {
  nx <- length(x)
  ny <- length(y)
  sy <- seq_along(y) - 1L
  z <- numeric(nx + ny - 1)
  for(i in seq(length = nx)) {
    ij <- i + sy
    t  <- duplicate(y)
    br_fmadd(t, x[[i]], z[ij])
    br_copy(z, t, xi = i)
  }
  z
}


#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
# by-reference 5
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
conv5 <- function(x,y) {
  nx <- length(x)
  ny <- length(y)
  sy <- seq_along(y) - 1L
  z <- numeric(nx + ny - 1)
  t <- duplicate(y)
  for(i in seq(length = nx)) {
    br_copy(t, y)
    br_fmadd(t, x[[i]], z[i + sy])
    br_copy(z, t, xi = i)
  }
  z
}


#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
# by-reference 6
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
conv6 <- function(x,y) {
  nx <- length(x)
  ny <- length(y)
  sy <- seq_along(y) - 1L
  z <- numeric(nx + ny - 1)
  ty <- duplicate(y)
  tz <- duplicate(y)
  for(i in seq(length = nx)) {
    br_copy(ty, y)
    br_copy(tz, z, n = ny, xi = 1, yi = i)
    br_fmadd(ty, x[[i]], tz)
    br_copy(z, ty, xi = i)
  }
  z
}




set.seed(1)
N <- 1000
x <- runif(N)
y <- runif(N)

if (FALSE) {
  conv0(x, y)
  conv1(x, y)
  conv2(x, y)
  conv3(x, y)
  conv4(x, y)
  convv(x, y)
  convolve(x, rev(y), conj = TRUE, type = 'open')
}

conv0c <- compiler::cmpfun(conv0)
# conv1c <- compiler::cmpfun(conv1)
conv2c <- compiler::cmpfun(conv2)
conv3c <- compiler::cmpfun(conv3)
conv4c <- compiler::cmpfun(conv4)
conv5c <- compiler::cmpfun(conv5)
conv6c <- compiler::cmpfun(conv6)
convvc <- compiler::cmpfun(convv)

bench::mark(
  conv0c(x, y),
  # conv1c(x, y),
  conv2c(x, y),
  conv3c(x, y),
  conv4c(x, y),
  conv5c(x, y),
  conv6c(x, y),
  convvc(x, y),
  convolve(x, rev(y), conj = TRUE, type = 'open')
)[, 1:5] |>
  knitr::kable()



bench::mark(
  conv0c(x, y),
  conv6c(x, y),
  convolve(x, rev(y), conj = TRUE, type = 'open')
)[, 1:5] |>
  knitr::kable()