Dear all,
Somebody kindly pointed out a problem in my WaveThresh3 code.
I can't figure out what is wrong.
I have whittled down the code quite a bit into an example case which
repeats the problem.
There is one R function called "ScalingFunction" and one C function
called "CScalFn.c".
The idea is that the R function calls the C routine repeatedly to
compute a wavelet scaling function.
It *seems* that after some iterations of the C routine later on in the
code some NAs mysteriously enter into some of the vectors (e.g. v) which
then gets rejected by the next call to the C routine.
I have no idea where these ghostly NAs have arrived from.
I feel as though I've done something stupid but I cannot see what it is.
Code below. I would appreciate some R brainbox letting me know where
I've gone wrong. (I already looked at the FAQ, the mailing lists and the
bug lists [a bit]).
My setup: R1.7.0 on RedHat Linux (2.4.18-14smp)
gcc version 3.2 20020903 (Red Hat Linux 8.0 3.2-7)
Although the problem was reported to me on R1.6.2 in Linux 7.1
Here is the C code (save as CScalFn.c and compiled using R CMD SHLIB,
then dyn.load("CScalFn.so") )
Many thanks,
Guy
----------------------
#include <R.h>
#define MAX(a,b) ( (a) < (b) ? (b) : (a))
#define MIN(a,b) ( (a) < (b) ? (a) : (b))
void CScalFn(Sfloat *v, Sfloat *ans, Sint *res, Sfloat *H, Sint *lengthH)
{
int k,n;
Sfloat sum;
Sint b,e;
for(n=0; n< *res; ++n) {
sum = 0.0;
b = MAX(0, (n+1- *lengthH)/2);
e = MIN(*res, n/2);
/* if (n < 100) Rprintf("%d %d\n", b,e); */
for(k=b; k<= e; ++k) {
sum += *(H+n-2*k) * *(v+k);
}
*(ans+n) = sum;
/*
if (sum > 0.0)
Rprintf("sum %lf\n", sum);*/
}
}
---------------------------------------------------------
Here is the R function (I've inserted various "cat" statements to
count
the number of NA, nan and Infs at various stages)
---------------------------------------------------------
"ScalingFunction" <-
function(resolution = 4096, itlevels = 50)
{
res <- 4 * resolution #
#
# Select filter and work out some fixed constants
#
H <- c(-1/sqrt(2), 1/sqrt(2))
lengthH <- length(H)
ll <- lengthH
v <- rep(0, res) #
#
# Set initial coefficient to 1 in 2nd position on 1st level
#
v[2] <- 1 #
#
# Now iterate the successive filtering operations to build up the scaling
# function. The actual filtering is carried out by the C routine CScalFn.
#
for(it in 1:itlevels) {
ans <- rep(0, res)
cat("Before: ", sum(is.na(v)), " ", sum(is.nan(v)),
" ",
sum(is.infinite(v)),"\n")
z <- .C("CScalFn",
v = as.double(v),
ans = as.double(ans),
res = as.integer(res),
H = as.double(H),
lengthH = as.integer(lengthH)) #
cat("After: ", sum(is.na(v)), " ", sum(is.nan(v)),
" ",
sum(is.infinite(v)),"\n")
#
# We only ever take the first half of the result
#
v <- z$ans[1:(res/2)] #
cat("After taking first half: ", sum(is.na(v)), " ",
sum(is.nan(v)), " ", sum(is.infinite(v)),"\n")
#
# Set all other coefficients equal to zero. (This is because
# rounding errors sometimes cause small values to appear).
#
v[ - ((2^it + 1):(2^it + ll))] <- 0
cat("After Setting all other coefficients equal to zero: ",
sum(is.na(v)), " ", sum(is.nan(v)), " ",
sum(is.infinite(v)),"\n")
v <- sqrt(2) * v
cat("After mult by root 2: ", sum(is.na(v)), " ",
sum(is.nan(v)), " ", sum(is.infinite(v)),"\n")
llbef <- ll
vbef <- v #
#
# Check to see if the next iteration would send the number
# of coefficients over the resolution that we can have.
# Exit the loop if it does.
#
if(2^(it + 1) + lengthH + ll * 2 - 2 > res/2) {
cit <- it
cat("Breaking : ", sum(is.na(v)), " ",
sum(is.nan(v)), " ",
sum(is.infinite(v)),"\n")
break
}
#
#
# ll is the number of coefficients that are nonzero in
# any particular run. This formula updates ll for next time
# round.
#
ll <- lengthH + ll * 2 - 2 #
#
# Add some zeroes to v to make it the right length.
#
v <- c(v, rep(0, res - length(v)))
cat("After adding some zeros : ", sum(is.na(v)), "
",
sum(is.nan(v)), " ", sum(is.infinite(v)),"\n")
}
list(x = seq(from = 0, to = 1, length = llbef), y
= vbef[(2^cit + 1):(2^cit + llbef)])
}