// @(#)root/base:$Id: TVirtualFFT.cxx 22967 2008-04-03 14:32:25Z rdm $
// Author: Anna Kreshuk  10/04/2006

/*************************************************************************
 * Copyright (C) 1995-2000, Rene Brun and Fons Rademakers.               *
 * All rights reserved.                                                  *
 *                                                                       *
 * For the licensing terms see $ROOTSYS/LICENSE.                         *
 * For the list of contributors see $ROOTSYS/README/CREDITS.             *
 *************************************************************************/

//////////////////////////////////////////////////////////////////////////
//
// TVirtualFFT
//
// TVirtualFFT is an interface class for Fast Fourier Transforms.
//
//
//
// The default FFT library is FFTW. To use it, FFTW3 library should already
// be installed, and ROOT should be have fftw3 module enabled, with the directories
// of fftw3 include file and library specified (see installation instructions).
// Function SetDefaultFFT() allows to change the default library.
//
// Available transform types:
// FFT:
// - "C2CFORWARD" - a complex input/output discrete Fourier transform (DFT)
//                  in one or more dimensions, -1 in the exponent
// - "C2CBACKWARD"- a complex input/output discrete Fourier transform (DFT)
//                  in one or more dimensions, +1 in the exponent
// - "R2C"        - a real-input/complex-output discrete Fourier transform (DFT)
//                  in one or more dimensions,
// - "C2R"        - inverse transforms to "R2C", taking complex input
//                  (storing the non-redundant half of a logically Hermitian array)
//                  to real output
// - "R2HC"       - a real-input DFT with output in ¡Èhalfcomplex¡É format,
//                  i.e. real and imaginary parts for a transform of size n stored as
//                  r0, r1, r2, ..., rn/2, i(n+1)/2-1, ..., i2, i1
// - "HC2R"       - computes the reverse of FFTW_R2HC, above
// - "DHT"        - computes a discrete Hartley transform
//
// Sine/cosine transforms:
// Different types of transforms are specified by parameter kind of the SineCosine() static
// function. 4 different kinds of sine and cosine transforms are available
//  DCT-I  (REDFT00 in FFTW3 notation)- kind=0
//  DCT-II (REDFT10 in FFTW3 notation)- kind=1
//  DCT-III(REDFT01 in FFTW3 notation)- kind=2
//  DCT-IV (REDFT11 in FFTW3 notation)- kind=3
//  DST-I  (RODFT00 in FFTW3 notation)- kind=4
//  DST-II (RODFT10 in FFTW3 notation)- kind=5
//  DST-III(RODFT01 in FFTW3 notation)- kind=6
//  DST-IV (RODFT11 in FFTW3 notation)- kind=7
// Formulas and detailed descriptions can be found in the chapter
// "What FFTW really computes" of the FFTW manual
//
// NOTE: FFTW computes unnormalized transforms, so doing a transform, followed by its
//       inverse will give the original array, multiplied by normalization constant
//       (transform size(N) for FFT, 2*(N-1) for DCT-I, 2*(N+1) for DST-I, 2*N for
//       other sine/cosine transforms)
//
// How to use it:
// Call to the static function FFT returns a pointer to a fast fourier transform
// with requested parameters. Call to the static function SineCosine returns a
// pointer to a sine or cosine transform with requested parameters. Example:
// {
//    Int_t N = 10; Double_t *in = new Double_t[N];
//    TVirtualFFT *fftr2c = TVirtualFFT::FFT(1, &N, "R2C");
//    fftr2c->SetPoints(in);
//    fftr2c->Transform();
//    Double_t re, im;
//    for (Int_t i=0; i<N; i++)
//       fftr2c->GetPointComplex(i, re, im);
//    ...
//    fftr2c->SetPoints(in2);
//    ...
//    fftr2c->SetPoints(in3);
//    ...
// }
// Different options are explained in the function comments
//
//
//
//
//
//////////////////////////////////////////////////////////////////////////

#include "TROOT.h"
#include "TVirtualFFT.h"
#include "TPluginManager.h"
#include "TEnv.h"
#include "TError.h"

TVirtualFFT *TVirtualFFT::fgFFT    = 0;
TString      TVirtualFFT::fgDefault   = "";

ClassImp(TVirtualFFT)

//_____________________________________________________________________________
TVirtualFFT::~TVirtualFFT()
{
   //destructor
   if (this==fgFFT)
      fgFFT = 0;
}

//_____________________________________________________________________________
TVirtualFFT* TVirtualFFT::FFT(Int_t ndim, Int_t *n, Option_t *option)
{
//Returns a pointer to the FFT of requested size and type.
//Parameters:
// -ndim : number of transform dimensions
// -n    : sizes of each dimension (an array at least ndim long)
// -option : consists of 3 parts - flag option and an option to create a new TVirtualFFT
//         1) transform type option:
//           Available transform types are:
//           C2CForward, C2CBackward, C2R, R2C, R2HC, HC2R, DHT
//           see class description for details
//         2) flag option: choosing how much time should be spent in planning the transform:
//           Possible options:
//           "ES" (from "estimate") - no time in preparing the transform,
//                                  but probably sub-optimal  performance
//           "M"  (from "measure")  - some time spend in finding the optimal way
//                                  to do the transform
//           "P" (from "patient")   - more time spend in finding the optimal way
//                                  to do the transform
//           "EX" (from "exhaustive") - the most optimal way is found
//           This option should be chosen depending on how many transforms of the
//           same size and type are going to be done.
//           Planning is only done once, for the first transform of this size and type.
//         3) option allowing to choose between the global fgFFT and a new TVirtualFFT object
//           ""  - default, changes and returns the global fgFFT variable
//           "K" (from "keep")- without touching the global fgFFT,
//           creates and returns a new TVirtualFFT*. User is then responsible for deleting it.
// Examples of valid options: "R2C ES K", "C2CF M", "DHT P K", etc.


   Int_t inputtype=0, currenttype=0;
   TString opt = option;
   opt.ToUpper();
   //find the tranform flag
   Option_t *flag;
   flag = "ES";
   if (opt.Contains("ES")) flag = "ES";
   if (opt.Contains("M"))  flag = "M";
   if (opt.Contains("P"))  flag = "P";
   if (opt.Contains("EX")) flag = "EX";

   Int_t ndiff = 0;

   if (!opt.Contains("K")) {
      if (fgFFT){
         //if the global transform exists, check if it should be changed
         if (fgFFT->GetNdim()!=ndim)
            ndiff++;
         else {
            Int_t *ncurrent = fgFFT->GetN();
            for (Int_t i=0; i<ndim; i++){
               if (n[i]!=ncurrent[i])
                  ndiff++;
            }
         }
         Option_t *t = fgFFT->GetType();
         if (!opt.Contains(t)) {
            if (opt.Contains("HC") || opt.Contains("DHT"))
               inputtype = 1;
            if (strcmp(t,"R2HC")==0 || strcmp(t,"HC2R")==0 || strcmp(t,"DHT")==0)
               currenttype=1;

            if (!(inputtype==1 && currenttype==1))
               ndiff++;
         }
         if (ndiff>0){
            delete fgFFT;
            fgFFT = 0;
         }
      }
   }

   Int_t sign = 0;
   if (opt.Contains("C2CB") || opt.Contains("C2R"))
      sign = 1;
   if (opt.Contains("C2CF") || opt.Contains("R2C"))
      sign = -1;

   TVirtualFFT *fft = 0;
   if (opt.Contains("K") || !fgFFT) {
      TPluginHandler *h;
      TString pluginname;
      if (fgDefault.Length()==0) fgDefault="fftw";
      if (strcmp(fgDefault.Data(),"fftw")==0) {
         if (opt.Contains("C2C")) pluginname = "fftwc2c";
         if (opt.Contains("C2R")) pluginname = "fftwc2r";
         if (opt.Contains("R2C")) pluginname = "fftwr2c";
         if (opt.Contains("HC") || opt.Contains("DHT")) pluginname = "fftwr2r";
         if ((h=gROOT->GetPluginManager()->FindHandler("TVirtualFFT", pluginname))) {
            if (h->LoadPlugin()==-1) {
               ::Error("TVirtualFFT::FFT", "handler not found");
               return 0;
            }
            fft = (TVirtualFFT*)h->ExecPlugin(3, ndim, n, kFALSE);
            if (!fft) {
               ::Error("TVirtualFFT::FFT", "plugin failed to create TVirtualFFT object");
               return 0;
            }
            Int_t *kind = new Int_t[1];
            if (pluginname=="fftwr2r") {
               if (opt.Contains("R2HC")) kind[0] = 10;
               if (opt.Contains("HC2R")) kind[0] = 11;
               if (opt.Contains("DHT")) kind[0] = 12;
            }
            fft->Init(flag, sign, kind);
            if (!opt.Contains("K")) {
               fgFFT = fft;
            }
            delete [] kind;
            return fft;
         }
         else {
            ::Error("TVirtualFFT::FFT", "plugin not found");
            return 0;
         }
      }
   } else {
      //if the global transform already exists and just needs to be reinitialised
      //with different parameters
      if (fgFFT->GetSign()!=sign || !opt.Contains(fgFFT->GetTransformFlag()) || !opt.Contains(fgFFT->GetType())) {
         Int_t *kind = new Int_t[1];
         if (inputtype==1) {
            if (opt.Contains("R2HC")) kind[0] = 10;
            if (opt.Contains("HC2R")) kind[0] = 11;
            if (opt.Contains("DHT")) kind[0] = 12;
         }
         fgFFT->Init(flag, sign, kind);
         delete [] kind;
      }
   }
   return fgFFT;
}

//_____________________________________________________________________________
TVirtualFFT* TVirtualFFT::SineCosine(Int_t ndim, Int_t *n, Int_t *r2rkind, Option_t *option)
{
//Returns a pointer to a sine or cosine transform of requested size and kind
//
//Parameters:
// -ndim    : number of transform dimensions
// -n       : sizes of each dimension (an array at least ndim long)
// -r2rkind : transform kind for each dimension
//     4 different kinds of sine and cosine transforms are available
//     DCT-I    - kind=0
//     DCT-II   - kind=1
//     DCT-III  - kind=2
//     DCT-IV   - kind=3
//     DST-I    - kind=4
//     DST-II   - kind=5
//     DST-III  - kind=6
//     DST-IV   - kind=7
// -option : consists of 2 parts - flag option and an option to create a new TVirtualFFT
//         - flag option: choosing how much time should be spent in planning the transform:
//           Possible options:
//           "ES" (from "estimate") - no time in preparing the transform,
//                                  but probably sub-optimal  performance
//           "M"  (from "measure")  - some time spend in finding the optimal way
//                                  to do the transform
//           "P" (from "patient")   - more time spend in finding the optimal way
//                                  to do the transform
//           "EX" (from "exhaustive") - the most optimal way is found
//           This option should be chosen depending on how many transforms of the
//           same size and type are going to be done.
//           Planning is only done once, for the first transform of this size and type.
//         - option allowing to choose between the global fgFFT and a new TVirtualFFT object
//           ""  - default, changes and returns the global fgFFT variable
//           "K" (from "keep")- without touching the global fgFFT,
//           creates and returns a new TVirtualFFT*. User is then responsible for deleting it.
// Examples of valid options: "ES K", "EX", etc

   TString opt = option;
   //find the tranform flag
   Option_t *flag;
   flag = "ES";
   if (opt.Contains("ES")) flag = "ES";
   if (opt.Contains("M"))  flag = "M";
   if (opt.Contains("P"))  flag = "P";
   if (opt.Contains("EX")) flag = "EX";

   if (!opt.Contains("K")) {
      if (fgFFT){
         Int_t ndiff = 0;
         if (fgFFT->GetNdim()!=ndim || strcmp(fgFFT->GetType(),"R2R")!=0)
            ndiff++;
         else {
            Int_t *ncurrent = fgFFT->GetN();
            for (Int_t i=0; i<ndim; i++) {
               if (n[i] != ncurrent[i])
                  ndiff++;
            }

         }
         if (ndiff>0) {
            delete fgFFT;
            fgFFT = 0;
         }
      }
   }
   TVirtualFFT *fft = 0;
   if (!fgFFT || opt.Contains("K")) {
      TPluginHandler *h;
      TString pluginname;
      if (fgDefault.Length()==0) fgDefault="fftw";
      if (strcmp(fgDefault.Data(),"fftw")==0) {
         pluginname = "fftwr2r";
         if ((h=gROOT->GetPluginManager()->FindHandler("TVirtualFFT", pluginname))) {
            if (h->LoadPlugin()==-1){
               ::Error("TVirtualFFT::SineCosine", "handler not found");
               return 0;
            }
            fft = (TVirtualFFT*)h->ExecPlugin(3, ndim, n, kFALSE);
            if (!fft) {
               ::Error("TVirtualFFT::SineCosine", "plugin failed to create TVirtualFFT object");
               return 0;
            }
            fft->Init(flag, 0, r2rkind);
            if (!opt.Contains("K"))
               fgFFT = fft;
            return fft;
         } else {
            ::Error("TVirtualFFT::SineCosine", "handler not found");
            return 0;
         }
      }
   }

   //if (fgFFT->GetTransformFlag()!=flag)
   fgFFT->Init(flag,0, r2rkind);
   return fgFFT;
}

//_____________________________________________________________________________
TVirtualFFT* TVirtualFFT::GetCurrentTransform()
{
// static: return current fgFFT

   if (fgFFT)
      return fgFFT;
   else{
      ::Warning("TVirtualFFT::GetCurrentTransform", "fgFFT is not defined yet");
      return 0;
   }
}

//_____________________________________________________________________________
void TVirtualFFT::SetTransform(TVirtualFFT* fft)
{
// static: set the current transfrom to parameter

   fgFFT = fft;
}

//_____________________________________________________________________________
const char *TVirtualFFT::GetDefaultFFT()
{
// static: return the name of the default fft

   return fgDefault.Data();
}

//______________________________________________________________________________
void TVirtualFFT::SetDefaultFFT(const char *name)
{
   // static: set name of default fft

   if (fgDefault == name) return;
   delete fgFFT;
   fgFFT = 0;
   fgDefault = name;
}