---
title: "Migrating from fftwtools"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{Migrating from fftwtools}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r, include = FALSE}
knitr::opts_chunk$set(collapse = TRUE, comment = "#>")
has_fftwtools <- requireNamespace("fftwtools", quietly = TRUE)
```

zufft computes the same transforms as fftwtools without the FFTW system
library, and under a permissive licence. This vignette maps every fftwtools
call to its zufft equivalent. Where fftwtools is installed, each example also
checks that both give the same values.

```{r}
library(zufft)
```

## Two conventions to know

**Scaling.** fftwtools, like `stats::fft()`, never scales: an inverse after
a forward transform multiplies by the length. zufft follows `scipy.fft` and
divides the inverse by the length by default (`norm = "backward"`). Use
`norm = "none"` to get fftwtools' numbers exactly.

```{r}
x <- c(1, 4, 2, 8, 5, 7)
fft_c2c(fft_c2c(x), inverse = TRUE)                 # x again
fft_c2c(fft_c2c(x), inverse = TRUE, norm = "none")  # 6 * x, as fftwtools
```

**Half spectra.** fftwtools' `HermConj = 1` (its default) returns the full
spectrum of real input, and `HermConj = 0` the non-redundant half. In zufft
that is the choice between `fft_c2c()` and `fft_r2c()`.

```{r}
fft_c2c(x)  # full spectrum, as fftw(x) and fftw_r2c(x)
fft_r2c(x)  # first n %/% 2 + 1 values, as fftw_r2c(x, HermConj = 0)
```

To go back from a half spectrum you need the original length, because an
even and an odd length give the same half size. Both packages take it as
`n`:

```{r}
X <- fft_r2c(1:7)
fft_c2r(X, n = 7, norm = "none")   # as fftw_c2r(X, HermConj = 0, n = 7)
```

## One-dimensional transforms

| fftwtools | zufft |
| --- | --- |
| `fftw(x)` | `fft_c2c(x)` |
| `fftw(x, inverse = 1)` | `fft_c2c(x, inverse = TRUE, norm = "none")` |
| `fftw(x, HermConj = 0)` (real `x`) | `fft_r2c(x)` |
| `fftw(X, inverse = 1, HermConj = 0, n = n)` | `fft_c2r(X, n = n, norm = "none")` |
| `fftw_r2c(x)` | `fft_c2c(x)` |
| `fftw_r2c(x, HermConj = 0)` | `fft_r2c(x)` |
| `fftw_c2r(X, HermConj = 0, n = n)` | `fft_c2r(X, n = n, norm = "none")` |
| `fftw_c2r(X)` (full spectrum) | `Re(fft_c2c(X, inverse = TRUE, norm = "none"))` |
| `fftw_c2c(x, inverse = i)` | `fft_c2c(x, inverse = i, norm = "none")` |

`inverse = 0` and `inverse = 1` work in zufft too.

```{r, eval = has_fftwtools}
z <- complex(real = rnorm(8), imaginary = rnorm(8))
all.equal(fft_c2c(z, inverse = 1, norm = "none"), fftwtools::fftw_c2c(z, inverse = 1))
all.equal(fft_r2c(x), fftwtools::fftw_r2c(x, HermConj = 0))
all.equal(fft_c2r(X, n = 7, norm = "none"), fftwtools::fftw_c2r(X, HermConj = 0, n = 7))
```

## Columns of a matrix

fftwtools' `mvfftw*` functions transform each column. In zufft that is
`axes = 1`.

| fftwtools | zufft |
| --- | --- |
| `mvfftw(m)`, `mvfftw_c2c(m)` | `fft_c2c(m, axes = 1)` |
| `mvfftw_r2c(m, HermConj = 0)` | `fft_r2c(m, axes = 1)` |
| `mvfftw_c2r(M, HermConj = 0, n = n)` | `fft_c2r(M, n = n, axes = 1, norm = "none")` |
| `mvfftw(..., fftplanopt = 1)` | nothing to do: zufft caches its plans |

```{r, eval = has_fftwtools}
m <- matrix(rnorm(24), 6)
all.equal(fft_c2c(m, axes = 1), fftwtools::mvfftw(m))
M <- fft_r2c(m, axes = 1)
all.equal(fft_c2r(M, n = 6, axes = 1, norm = "none"), fftwtools::mvfftw_c2r(M, HermConj = 0, n = 6))
```

## Two, three and more dimensions

| fftwtools | zufft |
| --- | --- |
| `fftw2d(x)`, `fftw_c2c_2d(x)`, `fftw3d(x)`, `fftw_c2c_3d(x)`, `fftw_c2c_xd(x)` | `fft_c2c(x)` |
| `fftw_r2c_2d(x, HermConj = 0)` | `fft_r2c(x)` |
| `fftw_r2c_3d(x, HermConj = 0)` | `fft_r2c(x)` |

FFTW stores arrays row by row, so fftwtools' real 2-D and 3-D transforms
halve R's *first* dimension. zufft's real transforms do the same by default:
in R the first dimension is the contiguous one, and halving it is the fast
choice.

```{r, eval = has_fftwtools}
x2 <- matrix(rnorm(35), 5, 7)
all.equal(fft_r2c(x2), fftwtools::fftw_r2c_2d(x2, HermConj = 0))
all.equal(fft_c2c(x2), fftwtools::fftw2d(x2))
```

## What you gain

- **No system library.** zufft compiles its bundled ducc0 code, so
  `install.packages("zufft")` is enough, everywhere.
- **No GPL.** zufft is MIT-licensed and its bundled code BSD-3, so a package
  that imports it keeps its own licence.
- **Speed.** On the benchmark grid in the package's `bench/RESULTS.md`,
  zufft is faster than fftwtools on most cases and within 1.5 times
  everywhere, because it caches its plans while fftwtools plans on every
  call. Multi-dimensional complex transforms gain the most.
- **Threads.** `workers = 2` (or `fft_workers(2)`) splits a transform over
  threads; fftwtools has none.
- **More transforms.** `fft_dct()` and `fft_dst()` (types 1 to 4), any
  subset of axes, cropping and zero-padding with `n`, and `fft_next_len()`
  for fast lengths.
- **Untouched inputs.** `fft_c2r()` never modifies its input, while FFTW's
  complex-to-real transform destroys it.