Dot products of vectors to perform 1D convolution

General
1

(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!

Thanks so much in advance!

1 Like
2

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
3

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]);

instead of

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!

4

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

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

Sorry for resuming an old thread, but I am currently facing the same issue and I am not sure I understand correctly the solution. Doesn’t tail just report the last elements of the vector?
Is there anything wrong in doing instead the following?

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

I haven’t tried it but yeah that looks like it should work.

And yeah I don’t think tail(b, i) is equivalent to b[i:1], so I’m not sure if that really solved the original question (but maybe I’m missing something). I think your b[sort_desc(1:i)] makes more sense.

2 Likes
7

I had the same issue, I agree tail(b, i) is not same as b[i:1] in R syntax. dot_product(a[1:i], b[sort_desc(1:i)]); looks right, but I tried and got error for syntax saying PARSER EXPECTED.
image

Cmdstan indeed has a function reverse() 5.16 Reverse functions | Stan Functions Reference but not available in rstan

8

This seems working but I am not sure if it will slow down the calculation

  for(i in 2:lxy){
    int asc[i] = 1:i;
    cvxy[i] = dot_product(x[asc], y[sort_desc(asc)]);
  }

What confuses me is the following also works even if asc is not used,

  for(i in 2:lxy){
    int asc[i] = 1:i;
    cvxy[i] = dot_product(x[1:i], y[sort_desc(1:i)]);
  }

but as long as I comment out line int asc[i] = 1:i;, I start having syntax error. I guess I am doing a silly mistake here, but could not find why…

1 Like
9

Sorry, I just noticed the following did not report syntax error in Rstudio IDE but will when compiling it. It seems still not legitimate. How should I code a reversed vector with the latest version of rstan? The only way is to do a for loop? Thanks.

  for(i in 2:lxy){
    int asc[i] = 1:i;
    cvxy[i] = dot_product(x[asc], y[sort_desc(asc)]);
  }
1 Like
10

Sorry for being slow to the party! I solved the reversal issue with the following code as a function

vector stan_rev(int nPoints, vector a) {
  
  vector [nPoints] a_rev;
  
  for(i in 1:nPoints)  
    a_rev[i] = a[nPoints - i + 1];
 
  return(a_rev);
}

… and then used that in a convolution function

vector stan_convolve(int nPoints, vector a,  vector b, real stepsize) {
  
  vector [nPoints] out;
  vector [nPoints] b_rev;
  
  b_rev = stan_rev(nPoints, b);
  
  for(i in 1:nPoints)  
    out[i] = dot_product(head(a, i), tail(b_rev, i)) * stepsize;
 
  return(out);
}

I know there are different ways of performing the convolution in different scenarios, but this is the one that I needed. Hope this helps!!

My main issue with using the convolution in general was that it ended up being really slow, as the whole analytical convolution would need to be performed over 1024 points (in my case) for each iteration of each sample . So everything got out of hand quite quickly.

4 Likes
11

Since Stan 2.24 there is also a built-in reverse function for vector, row_vector and arrays : 5.16 Reverse functions | Stan Functions Reference

If not, you can simplify the function to not require nPoints:

vector stan_rev(vector a) {
  int a_size = num_elements(a);
  vector [a_size] a_rev;
  
  for(i in 1:a_size)  
    a_rev[i] = a[a_size - i + 1];
 
  return a_rev;
}
4 Likes
12

Thanks for replying @mathesong and @rok_cesnovar. Yes, the latest rstan has a reverse() function, the syntax error I posted here were also disappeared after I updated to 2.26. @martinmodrak helped to resolve this problem yesterday in this post Syntax errors for a user defined function - #9 by martinmodrak. Sorry, I forgot to update the solution here.

1 Like