caseyk - 1 year ago 61

C++ Question

I have noticed that

`RcppArmadillo`

`ifft2`

`RcppArmadillo`

`mvfft(..., inverse = TRUE)`

I have debugged the issue specifically to the

`ifft(arma::cx_mat input)`

Example:

`ifft2`

`[1] 0.513297156-0.423498014i -0.129250939+0.300225299i`

0.039722228-0.093052563i -0.007956237+0.018643534i 0.001181177-0.002768473i

`mvfft`

`[1] 0.278131988-0.633838170i -0.195699114+0.445980950i`

0.060070320-0.136894940i -0.011924932+0.027175865i 0.001754788-0.003999007i

Questions

- Is the FFT still in development?
`RcppArmadillo`

- Is this a common issue across FFT variants(numerical deviations outside of FP or DP noise)?
- Is there a 'low-level' function call from Rcpp or RcppArmadillo to call R's native FFT?

Reproducibility - Below I condensed the problem as much as I could and reproduced the issue.

Rcpp code:

`#include <RcppArmadillo.h>`

// [[Rcpp::depends(RcppArmadillo)]]

using namespace Rcpp;

// [[Rcpp::export]]

//profile is the dependent variable of a given variable x,

//q is a vector containing complex valued information for a single column after a tcrossprod

//Size is a scalar value which the FFT depends upon.

arma::cx_mat DebugLmnCPP( arma::cx_vec Profile, arma::cx_vec q) {

std::complex<double> oneeye (0,1);//Cmplx number (0 + 1i)

arma::cx_mat qFFT = ifft2( exp( oneeye * (Profile * q.st() ) ) );

return(qFFT );

}

// [[Rcpp::export]]

//For pedagogical purposes

arma::cx_mat DebugIFFTRCPP( arma::cx_mat input) {

arma::cx_mat qFFT = ifft2( input );

return( qFFT );

}

RCode (sorry this is sloppy)

`library(Rcpp)`

library(RcppArmadillo)

sourceCpp("/home/FILE.cpp")

#Use C++ function

qt <- c(6.0+0i, 5.95+0i, 0.10+0i)

prof <- 0.25* sin( (1:512)*(2*3.1415)/512 ) + 0.25#Offset Sine wave

Debug1 <- DebugLmnCPP( Profile = prof, q = qt )

#Use R function

FFTSize <- 2^9

DebugLmnR <- function(Profile, q) {

g <- (0+1i)*(as.matrix(Profile ) %*% t(q))

qFFT <- mvfft( exp(g) , inverse = TRUE) / FFTSize

return( qFFT )

}

#Call function

Debug2 <- DebugLmnR( Profile = prof, q = qt )

#Use R and C++

DebugLmnRC <- function(Profile, q) {

g <- (0+1i)*(as.matrix(Profile ) %*% t(q))

qFFT <- DebugIFFTRCPP(exp(g))

return( qFFT )

}

#Call function

Debug3 <- DebugLmnRC( Profile = prof, q = qt )

#Compare Results

Debug1[1:5,1] #CPP

Debug2[1:5,1] #R

Debug3[1:5,1] #R and CPP

yields :

`> Debug1[1:5,1]`

[1] 0.359632774+0.35083419i -0.037254305-0.36995074i 0.015576046+0.15288379i -0.004552119-0.03992962i

[5] 0.000967252+0.00765564i

> Debug2[1:5,1]

[1] 0.03620451+0.51053116i -0.04624384-0.55604273i 0.02204910+0.23101589i -0.00653108-0.06061692i

[5] 0.00140213+0.01167389i

> Debug3[1:5,1]

[1] 0.359632774+0.35083419i -0.037254305-0.36995074i 0.015576046+0.15288379i -0.004552119-0.03992962i

[5] 0.000967252+0.00765564i

Answer Source

I don't particularly like your example:

- as it is still way too complex
- you are transforming data in the functions you are comparing -- generally a bad idea
- so I would suggest you fix your inputs

Here is a simpler example. `help(fft)`

in R leads of with this example

```
fftR> x <- 1:4
fftR> fft(x)
[1] 10+0i -2+2i -2+0i -2-2i
fftR> fft(fft(x), inverse = TRUE)/length(x)
[1] 1+0i 2+0i 3+0i 4+0i
```

which we can easily reproduce using RcppArmadillo:

```
R> cppFunction("arma::cx_mat armafft(arma::vec x) { return fft(x); }",
+ depends="RcppArmadillo")
R> armafft(1:4)
[,1]
[1,] 10+0i
[2,] -2+2i
[3,] -2+0i
[4,] -2-2i
R>
```

and adding the inverse

```
R> cppFunction("arma::cx_mat armaifft(arma::cx_mat x) { return ifft(x); }",
+ depends="RcppArmadillo")
R> armaifft(armafft(1:4))
[,1]
[1,] 1+0i
[2,] 2+0i
[3,] 3+0i
[4,] 4+0i
R>
```

recovering our input as in the R example.

No bug as far as I can tell, and I have no reason to believe this is any different for the 2d case...

**Edit/Followup:** The error is with the OP, and not with Armadillo. The primary issues here are

- not reading the documentation carefully
- not working with a minimal examples

The main issue here is that *Armadillo's fft() can work on vectors or matrices* and does hence (in the matrix case) correspond to R's

`mvfft()`

. Armadillo's `fft2()`

is simply something else and not relevant here.Let us continue / extend our previous example. We redefine our accessor to use complex matrix values:

```
R> cppFunction("arma::cx_mat armafft(arma::cx_mat x) { return fft(x); }",
+ depends="RcppArmadillo")
R>
```

and then define a complex array of dimension 5 x 2 which we feed to it:

```
R> z <- array(1:10 + 1i, dim=c(5,2))
R> z
[,1] [,2]
[1,] 1+1i 6+1i
[2,] 2+1i 7+1i
[3,] 3+1i 8+1i
[4,] 4+1i 9+1i
[5,] 5+1i 10+1i
R>
R> armafft(z)
[,1] [,2]
[1,] 15.0+5.00000i 40.0+5.00000i
[2,] -2.5+3.44095i -2.5+3.44095i
[3,] -2.5+0.81230i -2.5+0.81230i
[4,] -2.5-0.81230i -2.5-0.81230i
[5,] -2.5-3.44095i -2.5-3.44095i
R>
```

This is the same output we would get from running the function separately on each column. And that is also what R does for `mvfft()`

(cf `help(fft)`

)

```
R> mvfft(z)
[,1] [,2]
[1,] 15.0+5.00000i 40.0+5.00000i
[2,] -2.5+3.44095i -2.5+3.44095i
[3,] -2.5+0.81230i -2.5+0.81230i
[4,] -2.5-0.81230i -2.5-0.81230i
[5,] -2.5-3.44095i -2.5-3.44095i
R>
```

Same result, different libraries / packages, no bug as far as I can see.