// X-Midas APIs (midaslib) implemented in C++
//
// Author:  Neon Ngo

#include <string>
#include <cctype>
#include <stdexcept>
#include <sstream>
#include <iomanip>
#include <limits>

#include "midaslib.h"
#include "midas.h"
#include "xmtime.h"

using namespace std;

/* Function Prototypes */
int_4 radix_check_slide (int_4 nfft, int_4 searchdir);
int_4 m_fftsize_valid (int_4 nfft, int_4 rcflag, int_4 searchdir);

// XM410 fjackets.cc implementation until auto generated (BUG 2369)
// CXXfunc I4 m$fftsize_valid <I4 <I4 <I4
extern "C" int_4 m$fftsize_valid_ (const int_4& a0, const int_4& a1, const int_4& a2)
{
  try {
    return m_fftsize_valid(a0, a1, a2);
  } catch (const std::exception& e) {
    m$error_(e.what(), strlen(e.what()));
  } catch (...) {
    m$error_("Unknown exception caught...", 27);
  }
  return 0;
}

// based on XM445 impl in utilities.cc
int_4 m_bytes_reqd (int_4 bpe, int_4 lper)
{
  // Check to see if this would overflow an int_4. Do the arithmetic
  // in int_8, which will safely contain the result, and check before
  // casting down.
  int_8 size = int_8(bpe) * int_8(lper);

  // If bpe is negative, that will cause 'size' to be negative,
  // and thus we must compute an 8-byte aligned quantity of bytes
  // that will hold that many bits.
  if (size > 0) {
    size = (size + 7) / 8 * 8;
  } else {
    size = ((7 - size)/8 + 7) / 8 * 8;
  }

  if (size > std::numeric_limits<int_4>::max()  ||
      size < std::numeric_limits<int_4>::min()) {
    ostringstream errmsg;
    errmsg << "m_bytes_reqd(" << bpe << ", " << lper << "): "
	   << "return value would exceed range of int_4";
    throw invalid_argument(errmsg.str());
  }
  return static_cast<int_4>(size);
}

//***********************************************************************
// FUNCTION:
//    M$FFTSIZE_VALID()
//    m_fftsize_valid()
//
//  SYNOPSIS:
//    integer*4 function M$FFTSIZE_VALID(nfft,type,search)
//
//  ARGUMENTS:
//    integer*4 nfft   --- desired FFT size
//    integer*4 type   --- FFT type (0=Complex, 1=Real)
//    integer*4 search --- direction to search for valid FFT size
//                           Negative -> return size <= nfft
//                           Positive -> return size >= nfft
//                           Zero     -> return size closest to nfft
//
//  PURPOSE:
//    Returns nearest valid FFT size (relative to the specified FFT
//    size) which can be computed by the current MIDAS FFT library for calls
//    to RFFT and CFFT.  This function is based on the math libraries and how
//    they are filtered by the X-Midas RFFT/CFFT layers.  If changes are made
//    to either the math libraries or X-Midas at the RFFT/CFFT level it may
//    effect the performance of this function.
//
//  REVISIONS:
//    04/08/08    cberster     new code
//
int_4 m_fftsize_valid (int_4 nfft, int_4 rcflag, int_4 searchdir)
// check for low else accept input
{
  int_4 valid_argument;

#if defined(_FFTW) || defined(_IPP)
  //all values are acceptable except those < 2
  if (nfft < 2) {
    valid_argument = 2;
  } else {
    valid_argument = nfft;
  }
  return (valid_argument);

#elif defined(_NMFFT)
  // check range and rail accordingly else force to nearest even number based on direction
  if (nfft < 2) {
    valid_argument = 2;
  } else {
    valid_argument = nmfft_radix_check_slide(nfft, rcflag, searchdir);
  }
  return (valid_argument);

#elif defined(_IMKL)
  if (nfft < 2) {
    valid_argument = 2;
  } else {
    if (imkl_check_all_ffts() == USE_IMKL_NONPWR2) {
      // MKL can handle all FFT Sizes up to 2^31 - 1
      valid_argument = nfft;
    } else {
      int_4 status = m_power2ge(nfft); // check for power of 2
      status = pow(2.0,status);
      if (status == nfft) {
        valid_argument = nfft;
      } else {
        // not a power of 2 so use nmfft based checks
        valid_argument = nmfft_radix_check_slide(nfft, rcflag, searchdir);
      }
    }
  }
  return (valid_argument);

#elif defined(_MLIB)
  valid_argument = nfft;
  m_warning( "MLIB library not supported by FFTSIZE_VALID");
  return (valid_argument);

#elif defined(_CXML)
  valid_argument = nfft;
  m_warning("CXML library not supported by FFTSIZE_VALID");
  return (valid_argument);

#elif defined(_ACML)
  valid_argument = nfft;
  m_warning("ACML library not supported by FFTSIZE_VALID");
  return (valid_argument);

#else // Legacy
  if (nfft < 2) {
    valid_argument = 2;
  } else {
    if (searchdir > 0) {
      valid_argument = radix_check_slide(nfft, 1);
    } else if (searchdir < 0) {
      valid_argument = radix_check_slide(nfft, -1);
    } else if (searchdir == 0) {
      int_4 status = radix_check_slide(nfft, 1);
      int_4 status1 = radix_check_slide(nfft, -1);

      if ((status - nfft) > (nfft - status1))
        valid_argument = status1;
      else
        valid_argument = status;
    }
  }
  return (valid_argument);
#endif
}

#if defined(_IMKL)
using namespace midas;

template<class RETYPE, class CXTYPE>
IMKL_FFTPlan<RETYPE,CXTYPE>::IMKL_FFTPlan (int_4 points, int_4 mflag, int_4 datatype) :
  MidasFFTPlan(points, mflag, datatype),
  handle_(NULL),
  workvec_(NULL)
{
  float scale;
  long status;

  DFTI_CONFIG_VALUE mkl_precision;
  DFTI_CONFIG_VALUE mkl_domain;

  switch (datatype_) {
  case (MidasFFTPlan::REAL):
    mkl_domain = DFTI_REAL;
    mkl_precision = DFTI_SINGLE;
    break;
  case (MidasFFTPlan::COMPLEX):
    mkl_domain = DFTI_COMPLEX;
    mkl_precision = DFTI_SINGLE;
    break;
  case (MidasFFTPlan::DOUBLE):
    mkl_domain = DFTI_REAL;
    mkl_precision = DFTI_DOUBLE;
    break;
  case (MidasFFTPlan::CMPLX_DOUBLE):
    mkl_domain = DFTI_COMPLEX;
    mkl_precision = DFTI_DOUBLE;
    break;
  default:
    throw invalid_argument("Invalid datatype specified for IMKL FFTs.");
    break;
  }

  status = DftiCreateDescriptor(&(handle_), mkl_precision, mkl_domain, 1, points_);
  if (!DftiErrorClass(status, DFTI_NO_ERROR)) {
    throw runtime_error("Error creating IMKL FFT descriptor:" +
                        string(DftiErrorMessage(status)));
  }

  // set scale for tranform
  if (mflag_ < 0) {
    status=DftiSetValue(handle_, DFTI_BACKWARD_SCALE, 1.0);
  } else {
    if (mkl_domain == DFTI_REAL)
      scale = 2.0;
    else
      scale = 1.0;
    status = DftiSetValue(handle_, DFTI_FORWARD_SCALE, scale);
  }
  if (!DftiErrorClass(status,DFTI_NO_ERROR)) {
    throw runtime_error("Error setting scale value for IMKL FFT descriptor:" +
                        string(DftiErrorMessage(status)));
  }

  status = DftiSetValue(handle_, DFTI_PACKED_FORMAT, DFTI_CCS_FORMAT);
  if (!DftiErrorClass(status,DFTI_NO_ERROR)) {
    throw runtime_error("Error setting storage format for IMKL FFT descriptor:" +
                        string(DftiErrorMessage(status)));
  }

  status = DftiCommitDescriptor(handle_);
  if (!DftiErrorClass(status,DFTI_NO_ERROR)) {
    throw runtime_error("Error committing IMKL FFT descriptor:" +
                        string(DftiErrorMessage(status)));
  } else {
    // allocate work buffer and initially zero it...
    if (mkl_domain == DFTI_REAL) {
      workvec_ = new RETYPE[points_+2];
      if (workvec_ == NULL) {
        throw runtime_error("Unable to allocate sufficient memory for FFT working buffer");
      }
      std::fill_n(workvec_, (points_+2), 0.0);
    }
  }
}

template<class RETYPE, class CXTYPE>
IMKL_FFTPlan<RETYPE,CXTYPE>::~IMKL_FFTPlan (void)
{
  if (handle_ != NULL) {
    DftiFreeDescriptor(&(handle_));
  }
  if (workvec_ != NULL) {
    delete[] workvec_;
  }
}

template<class RETYPE, class CXTYPE>
void IMKL_FFTPlan<RETYPE,CXTYPE>::execute_real (void* afft)
{
  long status;
  RETYPE* rbuf = reinterpret_cast<RETYPE*>(afft);

  // if we're unpacked, copy to work buffer, remembering to copy *in* any possible nyquist sample;
  if (mflag_ > 0) {
    if (mflag_ > FFTPLAN_FWD_REAL_UNPACKED) {
      // DR 798551 - Removed the use of IMKL cblas_scpoy() which was using creating
      // multiple threads and using excessive CPU, replaced with memcpy()
      memcpy(workvec_, rbuf, points_*sizeof(RETYPE));
      status = DftiComputeForward(handle_, workvec_);
      memcpy(rbuf, workvec_, points_*sizeof(RETYPE));
      if (mflag_ == FFTPLAN_FWD_REAL_PACKED) {
        rbuf[1] = workvec_[points_]; // want packed data, stuff the nyquist bin back down...
      }
    } else {
      status = DftiComputeForward(handle_, rbuf);
    }
  } else {
    if (mflag_ < FFTPLAN_INV_REAL_UNPACKED) {
      memcpy(workvec_, rbuf, points_*sizeof(RETYPE));
      if (mflag_ == FFTPLAN_INV_REAL_PACKED) {
        workvec_[points_] = rbuf[1];
        workvec_[points_ + 1] = 0.0;
        rbuf[1] = 0.0;
      }
      status = DftiComputeBackward(handle_, workvec_);
      memcpy(rbuf, workvec_, points_*sizeof(RETYPE));
    } else {
      status = DftiComputeBackward(handle_, rbuf);
    }
  }
  if (!DftiErrorClass(status, DFTI_NO_ERROR)) {
    m_debug("Error computing m_rfft:" + string(DftiErrorMessage(status)));
  }
}

template<class RETYPE, class CXTYPE>
void IMKL_FFTPlan<RETYPE,CXTYPE>::execute_complex (void* afft)
{
  long status;
  if (mflag_ < 0) {
    status = DftiComputeBackward(handle_, afft);
  } else {
    status = DftiComputeForward(handle_, afft);
  }
  if (!DftiErrorClass(status, DFTI_NO_ERROR)) {
    m_debug("Error computing m_cfft:" + string(DftiErrorMessage(status)));
  }
}


MidasFFTPlan*
MidasFFTPlan::Create (void*, int_4 points,
                      int_4 mflag, int_4 datatype)
{
  MidasFFTPlan* plan;
  if (datatype == MidasFFTPlan::REAL || datatype == MidasFFTPlan::COMPLEX) {
    plan = new IMKL_FFTPlan<real_4,complex_8>(points, mflag, datatype);
  } else {
    plan = new IMKL_FFTPlan<real_8,complex_16>(points, mflag, datatype);
  }
  return plan;
}

void
MidasFFTPlan::Delete (const MidasFFTPlan* plan)
{
  delete plan;
}

void m_imkl_rfft (real_4* afft, int_4 nfft, int_4 mflag)
{
  MidasFFTPlan* plan = MidasFFTPlan::Find( nfft, mflag,
                                           MidasFFTPlan::REAL);
  plan->execute_real(afft);
}

void m_imkl_cfft (complex_8* afft, int_4 nfft, int_4 mflag)
{
  MidasFFTPlan* plan = MidasFFTPlan::Find( nfft, mflag,
                                           MidasFFTPlan::COMPLEX);
  plan->execute_complex(afft);
}

void m_imkl_cleanup()
{
  MidasFFTPlan::Clear();
}

#endif  // _IMKL

// based on XM461 impl in pipeio.cc
bool_4 m_filename_is_pipe (const string& filename)
{
  return (filename.size() > 0 && Mc->mode != 0 && filename[0] == '_');
}

//***********************************************************************
// FUNCTION:
//
//    radix_check_slide()
//
//  SYNOPSIS:
//    integer*4 function radix_check_slide(nfft,searchdir)
//
//  ARGUMENTS:
//    integer*4 nfft   --- desired FFT size
//    integer*4 search --- direction to search for valid FFT size
//                           Negative -> return size <= nfft
//                           Positive -> return size >= nfft
//
//
//  PURPOSE:
//    Returns nearest valid FFT size (relative to the specified FFT
//    size) in the searchdir direction using the radix 2,3,4,5 to MOD 0 algorithm
//
//  REVISIONS:
//    04/08/08    cberster     new code
//

int_4 radix_check_slide (int_4 nfft, int_4 searchdir)
{

  int_4 status, nfft_temp, nfft_temp_save, valid_argument;

  // must be divisible by 2 and then (1 or 3 or 4 or 5 or 2) to MOD 0

  status = 0;
  nfft_temp = nfft;

  while (status == 0) {

    if ((nfft_temp % 2) != 0) { //make sure we are first, always dealing with an even number
      if (searchdir > 0)
        nfft_temp = nfft_temp + 1;
      else
        nfft_temp = nfft_temp - 1;
    }

    nfft_temp_save = nfft_temp;

    while ( (nfft_temp % 2) == 0 ) nfft_temp /= 2;
    while ( (nfft_temp % 3) == 0 ) nfft_temp /= 3;
    while ( (nfft_temp % 4) == 0 ) nfft_temp /= 4;
    while ( (nfft_temp % 5) == 0 ) nfft_temp /= 5;

    status = nfft_temp == 1; // if nfft ==1  ->> 1, or if nfft=anything else ->> 0

    if (status == 0) {
      if (searchdir > 0)
        nfft_temp = nfft_temp_save + 1;
      else
        nfft_temp = nfft_temp_save - 1;
    }

  } //endwhile

  valid_argument = nfft_temp_save;

  return (valid_argument);
}

// m_utc ()
//
// Returns current time in UTC, adjusted from Unix Epoch of 
// January 1, 1970 to X-Midas Epoch of J1950.

void m_utc (real_8& seconds, real_8& fseconds)
{

  struct timespec ts;
  if (clock_gettime(CLOCK_REALTIME, &ts))
    throw std::runtime_error("unable to determine current time in "
                             "clock_gettime()");

  seconds = ts.tv_sec + Seconds_1950to1970;
  // Adjust microseconds to fractional seconds
  fseconds = ts.tv_nsec / 1e9;
   
}
