# To elaborate an example
set.seed(120)
<- rnorm(10, mean = 10, sd = 2)
x <- rexp(10, rate = 1.5)
y
# For Benchmarking
set.seed(120)
<- rnorm(5e3, mean = 10, sd = 2)
xx <- rexp(5e3, rate = 1.5) yy
5.1 Convolution
I have 2 types of data here, one being used as an example if they are actually equal or the function is working. While the other one is used for benchmarking.
#include <stdio.h>
#include <stdlib.h>
void convolve_c(double *x, double *y, int *n, int *m, double *result) {
int result_size = *n + *m - 1;
for (int i = 0; i < result_size; i++) {
int j_min = (i >= *m - 1) ? i - (*m - 1) : 0;
int j_max = (i < *n - 1) ? i : *n - 1;
[i] = 0;
resultfor (int j = j_min; j <= j_max; j++) {
[i] += x[j] * y[i - j];
result}
}
}
<- function(x, y) {
convolve_c <- length(x) + length(y) - 1
totalSize <- double(totalSize) # Initialize result vector
result <- .C("convolve_c",
result as.double(x), as.double(y),
as.integer(length(x)), as.integer(length(y)),
result=double(totalSize))$result
return(result)
}
convolve_c(x, y) |> head(10)
[1] 1.866342 10.531917 12.357334 14.054101 34.520841 39.345087 41.662725
[8] 55.926328 68.759844 65.767219
#include <Rcpp.h>
using namespace Rcpp;
// [[Rcpp::export]]
(NumericVector x, NumericVector y) {
NumericVector convolve_cppint n = x.size();
int m = y.size();
int result_size = n + m - 1;
(result_size);
NumericVector result
for (int i = 0; i < result_size; i++) {
int j_min = (i >= m - 1) ? i - (m - 1) : 0;
int j_max = (i < n - 1) ? i : n - 1;
double sum = 0.0;
// Unroll the inner loop for better performance
for (int j = j_min; j <= j_max - 4; j += 4) {
+= x[j] * y[i - j] +
sum [j + 1] * y[i - (j + 1)] +
x[j + 2] * y[i - (j + 2)] +
x[j + 3] * y[i - (j + 3)];
x}
// Handle remaining elements if any
for (int j = j_max - (j_max - j_min) % 4; j <= j_max; j++) {
+= x[j] * y[i - j];
sum }
[i] = sum;
result}
return result;
}
convolve_cpp(x, y) |> head(10)
[1] 1.866342 10.531917 12.357334 14.054101 34.520841 39.345087 41.662725
[8] 55.926328 68.759844 65.767219
function convolve_vector(x::Vector{Float64}, y::Vector{Float64})
= length(x)
n = length(y)
m = zeros(Float64, n + m - 1)
result = zeros(Float64, n + m - 1)
paddedX = zeros(Float64, n + m - 1)
paddedY
1:n] = x
paddedX[1:m] = y
paddedY[
for i in 1:(n + m - 1)
for j in 1:i
+= paddedX[j] * paddedY[i - j + 1]
result[i] end
end
return result
end
convolve_vector (generic function with 1 method)
<- JuliaCall::julia_eval("convolve_vector")
convolve_jl convolve_jl(x, y) |> head(10)
[1] 1.866342 10.531917 12.357334 14.054101 34.520841 39.345087 41.662725
[8] 55.926328 68.759844 65.767219
use extendr_api::prelude::*;
#[extendr]
fn convolve_rs(x: Vec<f64>, y: Vec<f64>) -> Vec<f64> {
let n = x.len();
let m = y.len();
let result_size = n + m - 1;
let mut result = vec![0.0; result_size];
for i in 0..result_size {
let j_min = if i >= m - 1 { i - (m - 1) } else { 0 };
let j_max = if i < n - 1 { i } else { n - 1 };
let mut sum = 0.0;
let mut j = j_min;
while j + 3 <= j_max {
+= x[j] * y[i - j]
sum + x[j + 1] * y[i - (j + 1)]
+ x[j + 2] * y[i - (j + 2)]
+ x[j + 3] * y[i - (j + 3)];
+= 4;
j }
for j in j..=j_max {
+= x[j] * y[i - j];
sum }
= sum;
result[i] }
result}
convolve_rs(x, y) |> head(10)
[1] 1.866342 10.531917 12.357334 14.054101 34.520841 39.345087 41.662725
[8] 55.926328 68.759844 65.767219
subroutine convolve_fortran(x, n, y, m, result)
implicit none
integer, intent(in) :: n, m
real(8), intent(in) :: x(n), y(m)
real(8), intent(out) :: result(n + m - 1)
real(8) :: paddedX(n + m - 1), paddedY(n + m - 1)
integer :: i, j
! Initialize paddedX and paddedY with zeros
paddedX = 0.0d0
paddedY = 0.0d0
! Copy elements of x and y to paddedX and paddedY, respectively
paddedX(1:n) = x
paddedY(1:m) = y
! Compute convolution
do i = 1, n + m - 1
do j = 1, i
result(i) = result(i) + paddedX(j) * paddedY(i - j + 1)
end do
end do
end subroutine convolve_fortran
<- function(x, y) {
convolve_fortran <- length(x) + length(y) - 1
totalSize <- .Fortran("convolve_fortran",
result as.double(x), as.integer(length(x)),
as.double(y), as.integer(length(y)),
result=double(totalSize))$result
return(result)
}
convolve_fortran(x, y) |> head(10)
[1] 1.866342 10.531917 12.357334 14.054101 34.520841 39.345087 41.662725
[8] 55.926328 68.759844 65.767219
<- function(x, y) {
convolve_r <- length(x)
n <- length(y)
m <- rep(0, n + m - 1)
result <- rep(0, n + m - 1)
paddedX <- rep(0, n + m - 1)
paddedY
1:n] <- x
paddedX[1:m] <- y
paddedY[
for (i in 1:(n + m - 1)) {
for (j in 1:i) {
<- result[i] + paddedX[j] * paddedY[i - j + 1]
result[i]
}
}
return(result)
}
convolve_r(x, y) |> head(10)
[1] 1.866342 10.531917 12.357334 14.054101 34.520841 39.345087 41.662725
[8] 55.926328 68.759844 65.767219
import numpy as np
def convolve_py(x, y):
= len(x)
n = len(y)
m = np.zeros(n + m - 1)
result = np.zeros(n + m - 1)
paddedX = np.zeros(n + m - 1)
paddedY
= x
paddedX[:n] = y
paddedY[:m]
for i in range(1, n + m):
for j in range(1, i+1):
-1] += paddedX[j-1] * paddedY[i-j]
result[i
return result
<- reticulate::py$convolve_py
convolve_py convolve_py(x, y) |> head(10)
[1] 1.866342 10.531917 12.357334 14.054101 34.520841 39.345087 41.662725
[8] 55.926328 68.759844 65.767219
Benchmark
::mark(
benchC = convolve_c(xx, yy),
`C++` = convolve_cpp(xx, yy),
Julia = convolve_jl(xx, yy),
Rust = convolve_rs(xx, yy),
FORTRAN = convolve_fortran(xx, yy),
R = convolve_r(xx, yy),
Python = convolve_py(xx, yy),
check = F
)
# A tibble: 7 × 6
expression min median `itr/sec` mem_alloc `gc/sec`
<bch:expr> <bch:tm> <bch:tm> <dbl> <bch:byt> <dbl>
1 C 40.1ms 41.48ms 23.8 356.2KB 0
2 C++ 12.65ms 15.71ms 62.5 85.8KB 0
3 Julia 80.2ms 87.65ms 11.2 83.9KB 0
4 Rust 591.53ms 591.53ms 1.69 83KB 0
5 FORTRAN 75.31ms 75.8ms 13.2 274KB 0
6 R 9.21s 9.21s 0.109 273.6KB 0
7 Python 47.62s 47.62s 0.0210 206.2KB 0