// GSL fft test with either std::vector or ublas::vector
// 
// Compile for std::vector :
// g++ gsl_fft_boost_vector_main.cpp -lgsl -L/usr/lib64/atlas -latlas -lcblas
// Compile for boost::numeric::ublas::vector :
// g++ -DBOOST gsl_fft_boost_vector_main.cpp -lgsl -L/usr/lib64/atlas -latlas -lcblas
//
#include <iostream>
#include <cmath>
#include <complex>
#include <algorithm>
#include <fstream>
#include "gsl_fft.hpp"
#include <iomanip>

using namespace std;
typedef complex<double> dtype ;

#ifdef BOOST
#include <boost/numeric/ublas/vector.hpp>
#include <boost/numeric/ublas/io.hpp>
#include <boost/numeric/ublas/vector_proxy.hpp>
typedef boost::numeric::ublas::vector<dtype> vec; 
#else
#include <vector>
typedef vector<dtype> vec; 
#endif


void vector_out(ostream& stream, vec v){
  vec::iterator vpos;
  stream<<setprecision(8);
  stream<<scientific;
  int i=0;
  for(vpos=v.begin();vpos!=v.end();++vpos){
    stream<<i++<<" "<<(*vpos).real()<<" "<<(*vpos).imag()<<endl;
  }
  stream<<endl;
}

void create_pulse(vec& v){
  const int nones= 10; //n/4; less than n!
  // pulse: 
  //     begin       nones+1                      end-nones        end
  //1    -------------------                      -----------------
  //0                       ----------------------
  // alternatively:
  //for(unsigned i=0;i<v.size();++i) v[i]=0.0;
  //for(unsigned i=0;i<nones+1;++i) v[i]=1.0;
  //for(unsigned i=v.size()-nones;i<v.size();++i) v[i]=1.0;
  //
  transform(v.begin(),v.begin()+nones+1,v.begin(),bind2nd(plus<dtype>(),1.0));
  transform(v.end()-nones,v.end(),v.end()-nones,bind2nd(plus<dtype>(),1.0));  
  //
}

void normalize_fft(vec& v){
   // alternatively:
  //double norm = sqrt(v.size());
  //for(unsigned i=0;i<v.size();++i) v[i] /= norm;
  transform(v.begin(),v.end(),v.begin(),bind2nd(divides<dtype>(),sqrt(v.size())));  
}

int main (void){
  const int n=128;
 
  vec data(n,0.0);
  ofstream myfile("result");

  create_pulse(data);
   
  cout<<"before fft :\n";
  vector_out(cout,data);
  vector_out(myfile,data);
  // ----------------------
  // forward FFT
  int status = fft(data,1);
  // ----------------------
  if(status!=0){
    cout << "fft fails\n";
    return 1;
  }
  cout<<"after forward fft (divided by sqrt(n)) :\n";
  normalize_fft(data);
  vector_out(cout,data);
  vector_out(myfile,data);

  // ----------------------
  // backward FFT
  status = fft(data,-1);
  // ----------------------
  if(status!=0){
    cout << "fft fails\n";
    return 1;
  }
  normalize_fft(data);
  cout<<"after backward fft (divided by sqrt(n)) :\n";
  vector_out(cout,data);

  myfile.close();
  return 0;
}