Fast fourier transform

[Michael Spence]

Hi

I have been trying to implement a fast Fourier transform in TMB. I tired to create an example below but I'm struggling. ChatGPT suggest that there's an issue that Eigen::FFT does not directly support the AD:

#include <TMB.hpp>
#include <complex>
#include <Eigen/Dense>
#include "C:/Users/MS23/Documents/cpp_packages/FFTW/fftw3.h" // location of FFTW package
#include <unsupported/Eigen/FFT>

template <class Type>
struct FFTTransformer {
  template <class Mat>
  Mat apply_fft(const Mat& X) {
    Mat X_transformed(X.rows(), X.cols());
    Eigen::FFT<typename Eigen::Array<std::complex<Type>, Eigen::Dynamic, 1>> fft;
    for (int i = 0; i < X.rows(); i++) {
      Eigen::Array<std::complex<Type>, Eigen::Dynamic, 1> fft_result(X.cols());
      fft.fwd(fft_result, X.row(i).template cast<std::complex<Type>>());
      X_transformed.row(i) = fft_result.matrix();
    }
    return X_transformed;
  }
};


template <class Type>
Type objective_function<Type>::operator()() {
  // Load data
  DATA_MATRIX(X);
 
  // Apply FFT transformation
  FFTTransformer<Type> transformer;
  matrix<Type> X_transformed = transformer.apply_fft(X);
 
  // Report transformed data
  REPORT(X_transformed);
 
  // Add your model code here
  PARAMETER(a);
  Type nll = pow(a, Type(2.0));
  return nll;
}

Is it possible to do this?

Cheers

Mike

[Ben Bolker]

I'm skimming the surface here, but is this a case where you can
write a new atomic function that defines the FFT *and* its derivative
(since everything is linear the derivative should be easy?
https://kaskr.github.io/adcomp/AtomicFunctions.html

> --
> To post to this group, send email to us...@tmb-project.org. Before
> posting, please check the wiki and issuetracker at
> https://github.com/kaskr/adcomp/ <https://github.com/kaskr/adcomp/>.
> Please try to create a simple repeatable example to go with your
> question (e.g issues 154, 134, 51). Use the issuetracker to report bugs.
> ---
> You received this message because you are subscribed to the Google
> Groups "TMB Users" group.
> To unsubscribe from this group and stop receiving emails from it, send
> an email to tmb-users+...@googlegroups.com
> <mailto:tmb-users+...@googlegroups.com>.
> To view this discussion on the web visit
> https://groups.google.com/d/msgid/tmb-users/d3e3f2df-2d50-4700-85a2-50982fd9cb25n%40googlegroups.com <https://groups.google.com/d/msgid/tmb-users/d3e3f2df-2d50-4700-85a2-50982fd9cb25n%40googlegroups.com?utm_medium=email&utm_source=footer>.

[Kasper Kristensen]

I agree with Ben that fft should be made as an atomic function. It is essentially a special case of atomic matrix multiply where the left operand is constant, so 'atomic::matmul' can be used as inspiration.
I would expect the adjoint code requires the inverse fft, so one would actually need two new atomics 'fft' and 'ifft'.
Also note that these atomics do not accept complex input/output. Manual splitting/unsplitting real- and imag parts would be needed.
Given that this is not completely straight-forward, I'm thinking it might be a useful addition to TMB. For portability reasons it's probably worth considering R's fft C-side entry point (if available).

[Ben Bolker]

For what it's worth, _Numerical Recipes_ (Press et al) give an
implementation of FFT (and inverse-FFT) for real, symmetric series (so
that all imaginary parts are zero going in either direction) that
packs/unpacks the input arrays so that the algorithm doesn't
unnecessarily carry the redundant half/zero imaginary values around.

The code should be available from the old editions of the books on
their web site; licensing doesn't allow you to directly cut and paste
the code, but if this is your actual use case it shouldn't be hard to
modify it enough to make it legal.

I don't know if that would be better/easier than using R's
implementation and dealing with the imaginary parts.

> <https://github.com/kaskr/adcomp/> <https://github.com/kaskr/adcomp/

> <https://github.com/kaskr/adcomp/>>.
> > Please try to create a simple repeatable example to go with your
> > question (e.g issues 154, 134, 51). Use the issuetracker to
> report bugs.
> > ---
> > You received this message because you are subscribed to the Google
> > Groups "TMB Users" group.
> > To unsubscribe from this group and stop receiving emails from it,
> send
> > an email to tmb-users+...@googlegroups.com
> > <mailto:tmb-users+...@googlegroups.com>.
> > To view this discussion on the web visit
> >

> https://groups.google.com/d/msgid/tmb-users/d3e3f2df-2d50-4700-85a2-50982fd9cb25n%40googlegroups.com <https://groups.google.com/d/msgid/tmb-users/d3e3f2df-2d50-4700-85a2-50982fd9cb25n%40googlegroups.com> <https://groups.google.com/d/msgid/tmb-users/d3e3f2df-2d50-4700-85a2-50982fd9cb25n%40googlegroups.com?utm_medium=email&utm_source=footer <https://groups.google.com/d/msgid/tmb-users/d3e3f2df-2d50-4700-85a2-50982fd9cb25n%40googlegroups.com?utm_medium=email&utm_source=footer>>.

>
> --
> To post to this group, send email to us...@tmb-project.org. Before
> posting, please check the wiki and issuetracker at
> https://github.com/kaskr/adcomp/ <https://github.com/kaskr/adcomp/>.
> Please try to create a simple repeatable example to go with your
> question (e.g issues 154, 134, 51). Use the issuetracker to report bugs.
> ---
> You received this message because you are subscribed to the Google
> Groups "TMB Users" group.
> To unsubscribe from this group and stop receiving emails from it, send
> an email to tmb-users+...@googlegroups.com
> <mailto:tmb-users+...@googlegroups.com>.
> To view this discussion on the web visit

> https://groups.google.com/d/msgid/tmb-users/add9260e-adac-4465-a2e1-31dfc20fa928n%40googlegroups.com <https://groups.google.com/d/msgid/tmb-users/add9260e-adac-4465-a2e1-31dfc20fa928n%40googlegroups.com?utm_medium=email&utm_source=footer>.

[Kasper Kristensen]

Mike: I think the 'Eigen::FFT'  works fine with AD. Although not complete, if I change your code to something like:

   Eigen::FFT<Type> fft;

    for (int i = 0; i < X.rows(); i++) {
      Eigen::Array<std::complex<Type>, Eigen::Dynamic, 1> fft_result(X.cols());

      Eigen::Array<Type, Eigen::Dynamic, 1> Xrow = X.row(i);
      fft.fwd(fft_result.data(), Xrow.data(), X.cols());
      //X_transformed.row(i) = fft_result.matrix();
    }

it compiles. Now you just need to fix the line that assigns complex fft_result to the real output X_transformed.

[Michael Spence]

Thanks a lot Kasper. I'll give that a try.

[Michael Spence]

Thanks a lot Kasper. I am interested in both the real and the imaginary parts. I added:

struct FFTTransformer {
  template <class Mat>
  Mat apply_fft(const Mat& X) {

    Mat X_transformed(X.rows(), 2 * X.cols()); // first half are real, next half imaginary

Eigen::FFT<Type> fft;
    for (int i = 0; i < X.rows(); i++) {
      Eigen::Array<std::complex<Type>, Eigen::Dynamic, 1> fft_result(X.cols());
      Eigen::Array<Type, Eigen::Dynamic, 1> Xrow = X.row(i);
      fft.fwd(fft_result.data(), Xrow.data(), X.cols());
 

  for (int j = 0; j < X.cols(); j++){
X_transformed(i,j) = fft_result[j].real();
X_transformed(i,X.cols() + j) = fft_result[j].imag();
  }
    }
    return X_transformed;
  }

};

and it worked fine.

Thanks a lot for your help!!

Mike