>>>>> GILLIBERT, Andre
>>>>> on Wed, 20 Oct 2021 08:10:00 +0000 writes:
> Hello,
> That sounds like a good diagnosis!
> Indeed, R vectors are passed "by reference" to C code, but
the semantic must be "by value", i.e. the C function must NOT change
the contents of the vector, except in very specific cases.
> A good program that has to work on a vector, must first duplicate the
vector, unless the only reference to the vector is the reference inside the C
function.
> This can be tested by the MAYBE_REFERENCED() macro in Rinternals.h.
> A good example can be found in the fft() function in
src/library/stats/src/fourier.c in R source code:
> switch (TYPEOF(z)) {
> case INTSXP:
> case LGLSXP:
> case REALSXP:
> z = coerceVector(z, CPLXSXP);
> break;
> case CPLXSXP:
> if (MAYBE_REFERENCED(z)) z = duplicate(z);
> break;
> default:
> error(_("non-numeric argument"));
> }
> PROTECT(z);
> This code coerces non-complex vectors to complex. Since this makes a
copy, there is no need to duplicate.
> Complex vectors are duplicated, unless they are not referenced by
anything but the fft() function.
> Now, the z vector can be modified "in place" without
inconsistency.
> Properly using R vectors in C code is tricky. You have to understand.
> 1) When you are allowed or not to modify vectors
> 2) When to PROTECT()vectors
> 3) How the garbage collector works and when it can trigger (answer :
basically, when you call any internal R function)
> Chapter 5 of "Writing R Extensions" documentation is quite
extensive:
>
https://cran.r-project.org/doc/manuals/r-release/R-exts.html#Handling-R-objects-in-C
> --
> Sincerely
> Andr? GILLIBERT
Thank you, Andr? , that's very good.
Just to state the obvious conclusion:
If Ben's suggestion is correct (and Andr? has explained *how*
that could happen) this would mean a
SEVERE BUG in package ravetools's mvfftw() function.
and it would have been (yet another) case of gaining speed by
killing correctness...
... but then ravetools is not even a CRAN package, so why
should you dare to use it for anything serious ?
... yes, being grouchy ..
> -----Message d'origine-----
> De?: R-devel <r-devel-bounces at r-project.org> De la part de Ben
Bolker
> Envoy??: mercredi 20 octobre 2021 03:27
> ??: r-devel at r-project.org
> Objet?: Re: [Rd] stats::fft produces inconsistent results
> This is a long shot, but here's a plausible scenario:
> as part of its pipeline, ravetools::mvfftw computes the mean of the
> input vector **and then centers it to a mean of zero** (intentionally
or
> accidentally?)
> because variables are passed to compiled code by reference (someone
> can feel free to correct my terminology), this means that the original
> vector in R now has a mean of zero
> the first element of fft() is mean(x)*length(x), so if mean(x) has
> been forced to zero, that would explain your issue.
> I don't know about the non-reproducibility part.
> On 10/19/21 7:06 PM, Dipterix Wang wrote:
>> Dear R-devel Team,
>>
>> I'm developing a neuroscience signal pipeline package in R
(https://github.com/dipterix/ravetools) and I noticed a weird issue that failed
my unit test.
>>
>> Basically I was trying to use `fftw3` library to implement fast
multivariate fft function in C++. When I tried to compare my results with
stats::fft, the test result showed the first element of **expected** (which was
produced by stats::fft) was zero, which, I am pretty sure, is wrong, and I can
confirm that my function produces correct results.
>>
>> However, somehow I couldn?t reproduce this issue on my personal
computer (osx, M1, R4.1.1), the error simply went away.
>>
>> The catch is my function produced consistent and correct results
but stats::fft was not. This does not mean `stats::fft` has bugs. Instead, I
suspect there could be some weird interactions between my code and stats::fft at
C/C++ level, but I couldn?t figure it out why.
>>
>> +++ Details:
>>
>> Here?s the code I used for the test:
>>
>>
https://github.com/dipterix/ravetools/blob/4dc35d64763304aff869d92dddad38a7f2b30637/tests/testthat/test-fftw.R#L33-L41
>>
>> ????????Test code????????
>> set.seed(1)
>> x <- rnorm(1000)
>> dim(x) <- c(100,10)
>> a <- ravetools:::mvfftw_r2c(x, 0)
>> c <- apply(x, 2, stats::fft)[1:51,]
>> expect_equal(a, c)
>> ????????????????????????
>>
>> Here are the tests that gave me the errors:
>>
>> The test logs on win-builder
>> https://win-builder.r-project.org/07586ios8AbL/00check.log
>>
>> Test logs on GitHub
>>
https://github.com/dipterix/ravetools/runs/3944874310?check_suite_focus=true
>>
>>
>> ?????????????? Failed tests ??????????????
>> -- Failure (test-fftw.R:41:3): mvfftw_r2c
--------------------------------------
>> `a` (`actual`) not equal to `c` (`expected`).
>>
>> actual vs expected
>> [,1] [,2] [,3]
[,4] ...
>> - actual[1, ] 10.8887367+ 0.0000000i -3.7808077+ 0.0000000i
2.967354+ 0.000000i 5.160186+ 0.000000i ...
>> + expected[1, ] 0.0000000+ 0.0000000i -3.7808077+ 0.0000000i
2.967354+ 0.000000i 5.160186+ 0.000000i...
>>
>> ????????????????????????
>>
>> The first columns are different, `actual` is the results I produced
via `ravetools:::mvfftw_r2c`, and `expected` was produced by `stats::fft`
>>
>>
>> Any help or attention is very much appreciated.
>> Thanks,
>> - Zhengjia