# [R] Weird problem with my code....

Guy Nason g.p.nason at bristol.ac.uk
Fri Sep 19 13:25:43 CEST 2003

```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)])
}

```

More information about the R-help mailing list