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

[Guy Nason]

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