[ADMB Users] Implementing a convolution in the GLOBAL_SECTION
Øigård Tor Arne
tor.arne.oeigaard at imr.no
Fri May 25 01:41:32 PDT 2012
Hi all,
I tried to implement a convolution in the GLOBAL_SECTION so that I can call a function conv(x,h) in the PROCEDURE_SECTION. The program compiles, but does not run, get the error message "size of file cmpdiff.tmp = 0" in the log.
I read in two vectors x and h in the DATA_SECTION and want to create a vector y = conv(x,h) in the PROCEDURE_SECTION. The variable y has been declared as a vector in the PARAMETER_SECTION. In the GLOBAL_SECTION i added:
#include <admodel.h>
dvar_vector conv1(const dvar_vector& x, const dvar_vector& h)
{
int k,j,m,n;
dvector lowerlimit(1,2);
dvector upperlimit(1,2);
dvar_vector y(1,n+m-1);
m=size_count(x);
n=size_count(h);
for(k=1;k<=(int) (m+n-1);k++)
{
y(k) = 0;
lowerlimit(1) = 1;
lowerlimit(2) = k+1-n;
upperlimit(1) = k;
upperlimit(2) = m;
for(j=max(lowerlimit);j<=(int) min(upperlimit);j++)
{
y(k) = y(k) + x(j)*h(k-j+1);
}
}
return y;
}
So the local variable y in the function is declared as a dvar_vector, and the function is declared as a dvar_vector. I wonder if the problem has something to do with that.
I tried to make another conv-function (below) which only returns one element of the output vector y, i.e., y_i, instead of the whole vector y, and then run a loop in the PROCEDURE_SECTION (the first loop in the conv1-function was moved to the PROCEDURE_SECTION). In this scenario the function conv2 is declared as a dvariable, and the local variable y (scalar) is declared as a dvariable. In this case it worked like a charm. Hence, it looks like I do something wrong when declaring the variables for the conv1-function. I also tried to replace the dvar_vector declaration for the conv1 and the local y vector with a dvector, but that did not help much.
dvariable conv2(const dvar_vector& x, const dvar_vector& h, const double& k)
{
int j,m,n;
dvector lowerlimit(1,2);
dvector upperlimit(1,2);
dvariable y;
m=size_count(x);
n=size_count(h);
y = 0;
lowerlimit(1) = 1;
lowerlimit(2) = k+1-n;
upperlimit(1) = k;
upperlimit(2) = m;
for(j=max(lowerlimit);j<=(int) min(upperlimit);j++)
{
y = y + x(j)*h(k-j+1);
}
return y;
}
Any good advice is very much appreciated.
Kind regards,
Tor Arne
More information about the Users
mailing list