 Dot products of vectors to perform 1D convolution

(I feel slightly embarrassed with this question, as I feel I’m probably just missing something stupid as I’m quite new to defining everything in STAN syntax directly, so my apologies in advance. If there’s someplace that there are a bunch of examples involving doing a bunch of manipulation of data in STAN syntax someplace where I might be able to get some pointers, I’d be very appreciative of the tip.)

I’m trying to write a function for performing 1D convolution of two curves for fitting of specific pharmacokinetic models where I need the interpolation over the sampled time points. In R, I usually use an FFT, but in STAN I’ll just have to resort to doing it in a for-loop. I’ve written a function for it, and it seems to compile fine, but I run into an error as soon as I try to execute it.

(I’m going to get the linear interpolation done outside of this function using Torsten’s linear_interpolation(), so by the time the curves arrive here, they should be interpolated so that they are equally spaced and measured at the same time on their x (time) axes. (Speaking of which, is there an equivalent of the seq() function in R? Otherwise, I’ll get back to that later…)).

R

So here could be some example input data (in R)

``````a <- c(0:10, 9:0)
b <- c( rep(0,5), rep(1,11), rep(0, 5))
nPoints = length(a)
``````

Here’s how I can do it in the for loop in R

``````r_manual_convolve <- function(a, b, stepsize) {

out <- rep(0, length(a))

for(i in 1:length(a)) {
out[i] <- sum( a[1:i] * b[i:1] )
}
return(out)
}

r_manual_convolve(a, b, 1)
``````

STAN

So here’s how I’ve tried to get the same function working in STAN syntax

``````functions {
vector stan_convolve(int nPoints, vector a,  vector b,  real stepsize) {

vector [nPoints] out;

for(i in 1:nPoints) {

vector [i] a_vec;
vector [i] b_vec;
vector [i] abprod;

a_vec = a[1:i];
b_vec = b[i:1];

abprod = a_vec .* b_vec;

out[i] = sum ( abprod );
}

return(out);

}
}
data{}
model{}
``````

This compiles fine. However, when I expose it to test it…

``````expose_stan_functions(stan.conv.code)
stan_convolve(nPoints, a, b, 1)
``````

… it seems to fail. I get the following error

``````Error in stan_convolve(nPoints, a, b, 1) :
Exception: assign: Rows of left-hand-side (2) and rows of right-hand-side (0) must match in size  (in 'unknown file name' at line 47)
``````

I just really don’t know where to begin in debugging this, whether it’s in the input data from R, or the information in the function itself. Any help would be thoroughly appreciated!

You could write this as

``````functions {
vector stan_convolve(int nPoints, vector a,  vector b) {

vector [nPoints] out;

for(i in 1:nPoints)
out[i] = dot_product(a[1:i], b[1:i]);

return(out);

}
``````
1 Like

Thanks so much for the response! Much smarter to use the dot_product() command - that’s a great improvement! And I was being silly to have forgotten to include the stepsize in both the R and STAN syntax (should just be a scalar to multiply the out vector), but I’ll fix that.

I tried your code and got slightly different answers. Then realised that you’d suggested

``````dot_product(a[1:i], b[1:i]);
``````

``````dot_product(a[1:i], b[i:1]); // note the reversed order of the b vector i:1
``````

So I merrily just switched this around, and the function no longer works, giving the following error:

``````Exception: dot_product: size of v1 (2) and size of v2 (0) must match in size
``````

So this has been massively helpful in making progress towards figuring out that the i:1 syntax is likely the issue here. I see that there was a GitHub issue and merged pull request (https://github.com/stan-dev/math/pull/479) about adding a `reverse()` function, but it doesn’t seem to be included in the latest stable release yet. Is there some kind of workaround for getting, e.g. indices 5:1 other than simply writing a whole for-loop to build it up item-by-item at the start? Or perhaps I just deal with a for-loop for now until the next release.

This has been extremely helpful! Thanks so much!

In Stan (but not R) the colon only yields increasing sequences. You could do

``````dot_product(head(a, i), tail(b, i));
``````
1 Like