From ripley at stats.ox.ac.uk  Sat Jul 28 07:54:10 2012
From: ripley at stats.ox.ac.uk (Prof Brian Ripley)
Date: Sat, 28 Jul 2012 06:54:10 +0100
Subject: [Rd] C code validation
In-Reply-To: <CAJ=0CtAQw+J5xMBHvnK4ttyRNfzWM2Ak=Q3u8m9HGZt9G0fnmQ@mail.gmail.com>
References: <CAJ=0CtAQw+J5xMBHvnK4ttyRNfzWM2Ak=Q3u8m9HGZt9G0fnmQ@mail.gmail.com>
Message-ID: <50137E82.8040706@stats.ox.ac.uk>

On 27/07/2012 23:52, Adrian Du?a wrote:
> Dear R-devel,
>
> I'm trying to validate the results from a C function, against a (trial
> and tested) older R function. For reasons unknown to me, the C
> function seems to give different result sometimes at each trial, even
> with the very same data.

See the advice in 'Writing R Extensions'.  In particular, the most 
common cause is the use of initialized values, so run under valgrind.

Also check that your C code does not change any of its arguments: if it 
does you may end up with self-modifying R code.  (That's probably rare, 
but I spent a couple of hours on such an example last night.)


>
> These are the relevant outputs from R:
>
>
>> library(QCA)
> Loading required package: lpSolve
>> benchmark <- function(x, y) {
> +     index <- 0
> +     while((index <- index + 1) < length(x)) {
> +         x <- setdiff(x, findSubsets(y + 1, x[index], max(x)))
> +     }
> +     return(x)
> + }
>>
>> set.seed(12345)
>> exprnec <- sort(sample(13:19683, 500))
>> noflevels <- rep(2, 9)
>> mbase <- rev(c(1, cumprod(rev(noflevels + 1))))[-1]
>> foo1 <- benchmark(exprnec, noflevels)
>> foo2 <- .Call("removeRedundants", exprnec, noflevels, mbase)
>>
>> all.equal(foo1, foo2)
> [1] "Numeric: lengths (184, 181) differ"
>>
>>
>> set.seed(12345)
>> exprnec <- sort(sample(13:19683, 500))
>> foo2 <- .Call("removeRedundants", exprnec, noflevels, mbase)
>> all.equal(foo1, foo2)
> [1] TRUE
>>
>>
>> set.seed(12345)
>> exprnec <- sort(sample(13:19683, 500))
>> foo2 <- .Call("removeRedundants", exprnec, noflevels, mbase)
>> all.equal(foo1, foo2)
> [1] TRUE
>>
>>
>> set.seed(12345)
>> exprnec <- sort(sample(13:19683, 500))
>> foo2 <- .Call("removeRedundants", exprnec, noflevels, mbase)
>> all.equal(foo1, foo2)
> [1] "Numeric: lengths (184, 173) differ"
>
>
> The relevant C code can be found here:
> http://codeviewer.org/view/code:287e
>
> As I am too baffled by these (to me strange) results, any hint is much
> appreciated.
>
> Thanks very much in advance,
> Adrian
>


-- 
Brian D. Ripley,                  ripley at stats.ox.ac.uk
Professor of Applied Statistics,  http://www.stats.ox.ac.uk/~ripley/
University of Oxford,             Tel:  +44 1865 272861 (self)
1 South Parks Road,                     +44 1865 272866 (PA)
Oxford OX1 3TG, UK                Fax:  +44 1865 272595



[Rd] C code validation
Prof Brian Ripley