text:
/*!
 *  Copyright (c) 2014 by Contributors
 * \file base.h
 * \brief definitions of base types, operators, macros functions
 *
 * \author Bing Xu, Tianqi Chen
 */
#ifndef MSHADOW_BASE_H_
#define MSHADOW_BASE_H_
#ifdef _MSC_VER
#ifndef _CRT_SECURE_NO_WARNINGS
#define _CRT_SECURE_NO_WARNINGS
#endif
#ifndef _CRT_SECURE_NO_DEPRECATE
#define _CRT_SECURE_NO_DEPRECATE
#endif
#define NOMINMAX
#endif
#include <cmath>
#include <cstdio>
#include <cfloat>
#include <climits>
#include <algorithm>
#include <functional>
#include <sstream>
#include <string>

#ifdef _MSC_VER
//! \cond Doxygen_Suppress
typedef signed char int8_t;
typedef __int16 int16_t;
typedef __int32 int32_t;
typedef __int64 int64_t;
typedef unsigned char uint8_t;
typedef unsigned __int16 uint16_t;
typedef unsigned __int32 uint32_t;
typedef unsigned __int64 uint64_t;
//! \endcond
#else
#include <inttypes.h>
#endif
// macro defintiions
/*!
 * \brief if this macro is define to be 1,
 * mshadow should compile without any of other libs
 */
#ifndef MSHADOW_STAND_ALONE
#define MSHADOW_STAND_ALONE 0
#endif
/*! \brief whether do padding during allocation */
#ifndef MSHADOW_ALLOC_PAD
#define MSHADOW_ALLOC_PAD true
#endif
/*!
 * \brief
 *  x dimension of data must be bigger pad_size * ratio to be alloced padded memory,
 *  otherwise use tide allocation
 *  for example, if pad_ratio=2, GPU memory alignement size is 32,
 *  then we will only allocate padded memory if x dimension > 64
 *  set it to 0 then we will always allocate padded memory
 */
#ifndef MSHADOW_MIN_PAD_RATIO
  #define MSHADOW_MIN_PAD_RATIO 2
#endif

#if MSHADOW_STAND_ALONE
  #define MSHADOW_USE_CBLAS 0
  #define MSHADOW_USE_MKL   0
  #define MSHADOW_USE_CUDA  0
#endif

/*!
 * \brief force user to use GPU stream during computation
 *  error will be shot when default stream NULL is used
 */
#ifndef MSHADOW_FORCE_STREAM
#define MSHADOW_FORCE_STREAM 1
#endif

/*! \brief use CBLAS for CBLAS */
#ifndef MSHADOW_USE_CBLAS
  #define MSHADOW_USE_CBLAS 0
#endif
/*! \brief use MKL for BLAS */
#ifndef MSHADOW_USE_MKL
  #define MSHADOW_USE_MKL   1
#endif

/*!
 * \brief use CUDA support, must ensure that the cuda include path is correct,
 * or directly compile using nvcc
 */
#ifndef MSHADOW_USE_CUDA
  #define MSHADOW_USE_CUDA   1
#endif

/*!
 * \brief use CUDNN support, must ensure that the cudnn include path is correct
 */
#ifndef MSHADOW_USE_CUDNN
  #define MSHADOW_USE_CUDNN 0
#endif

/*!
 * \brief seems CUDAARCH is deprecated in future NVCC
 * set this to 1 if you want to use CUDA version smaller than 2.0
 */
#ifndef MSHADOW_OLD_CUDA
#define MSHADOW_OLD_CUDA 0
#endif

/*!
 * \brief macro to decide existence of c++11 compiler
 */
#ifndef MSHADOW_IN_CXX11
#define MSHADOW_IN_CXX11 (defined(__GXX_EXPERIMENTAL_CXX0X__) ||\
                          __cplusplus >= 201103L || defined(_MSC_VER))
#endif

/*! \brief whether use SSE */
#ifndef MSHADOW_USE_SSE
  #define MSHADOW_USE_SSE 1
#endif
/*! \brief whether use NVML to get dynamic info */
#ifndef MSHADOW_USE_NVML
  #define MSHADOW_USE_NVML 0
#endif
// SSE is conflict with cudacc
#ifdef __CUDACC__
  #undef MSHADOW_USE_SSE
  #define MSHADOW_USE_SSE 0
#endif

#if MSHADOW_USE_CBLAS
extern "C" {
    #include <cblas.h>
}
#elif MSHADOW_USE_MKL
  #include <mkl_blas.h>
  #include <mkl_cblas.h>
  #include <mkl_vsl.h>
  #include <mkl_vsl_functions.h>
#endif

#if MSHADOW_USE_CUDA
  #include <cuda.h>
  #include <cublas_v2.h>
  #include <curand.h>
#endif

#if MSHADOW_USE_CUDNN == 1
  #include <cudnn.h>
#endif

#if MSHADOW_USE_NVML
  #include <nvml.h>
#endif

// --------------------------------
// MSHADOW_XINLINE is used for inlining template code for both CUDA and CPU code
#ifdef MSHADOW_XINLINE
  #error "MSHADOW_XINLINE must not be defined"
#endif
#ifdef _MSC_VER
#define MSHADOW_FORCE_INLINE __forceinline
#pragma warning(disable : 4068)
#else
#define MSHADOW_FORCE_INLINE inline __attribute__((always_inline))
#endif
#ifdef __CUDACC__
  #define MSHADOW_XINLINE MSHADOW_FORCE_INLINE __device__ __host__
#else
  #define MSHADOW_XINLINE MSHADOW_FORCE_INLINE
#endif
/*! \brief cpu force inline */
#define MSHADOW_CINLINE MSHADOW_FORCE_INLINE

#if defined(__GXX_EXPERIMENTAL_CXX0X) ||\
    defined(__GXX_EXPERIMENTAL_CXX0X__) || __cplusplus >= 201103L
  #define MSHADOW_CONSTEXPR constexpr
#else
  #define MSHADOW_CONSTEXPR const
#endif

/*!
 * \brief default data type for tensor string
 *  in code release, change it to default_real_t
 *  during development, change it to empty string so that missing
 *  template arguments can be detected
 */
#ifndef MSHADOW_DEFAULT_DTYPE
#define MSHADOW_DEFAULT_DTYPE = default_real_t
#endif

/*!
 * \brief DMLC marco for logging
 */
#ifndef MSHADOW_USE_GLOG
#define MSHADOW_USE_GLOG DMLC_USE_GLOG
#endif  // MSHADOW_USE_GLOG

#if DMLC_USE_CXX11
#define MSHADOW_THROW_EXCEPTION noexcept(false)
#define MSHADOW_NO_EXCEPTION  noexcept(true)
#else
#define MSHADOW_THROW_EXCEPTION
#define MSHADOW_NO_EXCEPTION
#endif

/*!
 * \brief Protected cuda call in mshadow
 * \param func Expression to call.
 * It checks for CUDA errors after invocation of the expression.
 */
#define MSHADOW_CUDA_CALL(func)                                    \
  {                                                                \
    cudaError_t e = (func);                                        \
    if (e == cudaErrorCudartUnloading) {                           \
      throw dmlc::Error(cudaGetErrorString(e));                    \
    }                                                              \
    CHECK(e == cudaSuccess)                                        \
        << "CUDA: " << cudaGetErrorString(e);                      \
  }

/*!
 * \brief Run function and catch error, log unknown error.
 * \param func Expression to call.
 */
#define MSHADOW_CATCH_ERROR(func)                                     \
  {                                                                   \
    try {                                                             \
      (func);                                                         \
    } catch (const dmlc::Error &e) {                                    \
      std::string what = e.what();                                      \
      if (what.find("driver shutting down") == std::string::npos) {     \
        LOG(ERROR) << "Ignore CUDA Error " << what;                     \
      }                                                                 \
    }                                                                   \
  }

#include "./half.h"
#include "./half2.h"
#include "./logging.h"
/*! \brief namespace for mshadow */
namespace mshadow {
/*! \brief buffer size for each random number generator */
const unsigned kRandBufferSize = 1000000;
/*! \brief pi  */
const float kPi = 3.1415926f;
/*! \brief type that will be used for index */
typedef unsigned index_t;

#ifdef _WIN32
  /*! \brief openmp index for windows */
  typedef int64_t openmp_index_t;
#else
  /*! \brief openmp index for linux */
  typedef index_t openmp_index_t;
#endif

/*! \brief float point type that will be used in default by mshadow */
typedef float default_real_t;

/*! \brief data type flag */
enum TypeFlag {
  kFloat32 = 0,
  kFloat64 = 1,
  kFloat16 = 2,
  kUint8 = 3,
  kInt32 = 4,
  kInt8  = 5
};

template<typename DType>
struct DataType;

template<>
struct DataType<float> {
  static const int kFlag = kFloat32;
  static const int kLanes = 1;
#if MSHADOW_USE_CUDA
#if (CUDA_VERSION >= 8000)
  static const cudaDataType_t kCudaFlag = CUDA_R_32F;
#endif
#if MSHADOW_USE_CUDNN
  static const cudnnDataType_t kCudnnFlag = CUDNN_DATA_FLOAT;
  typedef float ScaleType;
#endif
#endif
};
template<>
struct DataType<double> {
  static const int kFlag = kFloat64;
  static const int kLanes = 1;
#if MSHADOW_USE_CUDA
#if (CUDA_VERSION >= 8000)
  static const cudaDataType_t kCudaFlag = CUDA_R_64F;
#endif
#if MSHADOW_USE_CUDNN
  static const cudnnDataType_t kCudnnFlag = CUDNN_DATA_DOUBLE;
  typedef double ScaleType;
#endif
#endif
};
template<>
struct DataType<half::half_t> {
  static const int kFlag = kFloat16;
  static const int kLanes = 1;
#if MSHADOW_USE_CUDA
#if (CUDA_VERSION >= 8000)
  static const cudaDataType_t kCudaFlag = CUDA_R_16F;
#endif
#if MSHADOW_USE_CUDNN
  static const cudnnDataType_t kCudnnFlag = CUDNN_DATA_HALF;
  typedef float ScaleType;
#endif
#endif
};
template<>
struct DataType<half::half2_t> {
  static const int kFlag = kFloat16;
  static const int kLanes = 2;
};
template<>
struct DataType<uint8_t> {
  static const int kFlag = kUint8;
  static const int kLanes = 1;
#if MSHADOW_USE_CUDA
#if (CUDA_VERSION >= 8000)
  static const cudaDataType_t kCudaFlag = CUDA_R_8U;
#endif
#if (MSHADOW_USE_CUDNN == 1 && CUDNN_MAJOR >= 6)
  // no uint8 in cudnn for now
  static const cudnnDataType_t kCudnnFlag = CUDNN_DATA_INT8;
  typedef uint8_t ScaleType;
#endif
#endif
};
template<>
struct DataType<int8_t> {
  static const int kFlag = kInt8;
  static const int kLanes = 1;
#if MSHADOW_USE_CUDA
#if (CUDA_VERSION >= 8000)
  static const cudaDataType_t kCudaFlag = CUDA_R_8I;
#endif
#if (MSHADOW_USE_CUDNN == 1 && CUDNN_MAJOR >= 6)
  static const cudnnDataType_t kCudnnFlag = CUDNN_DATA_INT8;
  typedef int8_t ScaleType;
#endif
#endif
};
template<>
struct DataType<int32_t> {
  static const int kFlag = kInt32;
  static const int kLanes = 1;
#if MSHADOW_USE_CUDA
#if (CUDA_VERSION >= 8000)
  static const cudaDataType_t kCudaFlag = CUDA_R_32I;
#endif
#if (MSHADOW_USE_CUDNN == 1 && CUDNN_MAJOR >= 6)
  static const cudnnDataType_t kCudnnFlag = CUDNN_DATA_INT32;
  typedef int32_t ScaleType;
#endif
#endif
};

/*! \brief type enum value for default real type */
const int default_type_flag = DataType<default_real_t>::kFlag;

/*! layout flag */
enum LayoutFlag {
  kNCHW = 0,
  kNHWC,
  kCHWN,

  kNCW = 1 << 3,
  kNWC,
  kCWN,

  kNCDHW = 1 << 5,
  kNDHWC,
  kCDHWN
};

template<int layout>
struct LayoutType;

template<>
struct LayoutType<kNCHW> {
  static const index_t kNdim = 4;
#if (MSHADOW_USE_CUDA && MSHADOW_USE_CUDNN == 1 && CUDNN_MAJOR >= 4)
  static const cudnnTensorFormat_t kCudnnFlag = CUDNN_TENSOR_NCHW;
#else
  static const int kCudnnFlag = -1;
#endif
};

template<>
struct LayoutType<kNHWC> {
  static const index_t kNdim = 4;
#if (MSHADOW_USE_CUDA && MSHADOW_USE_CUDNN == 1 && CUDNN_MAJOR >= 4)
  static const cudnnTensorFormat_t kCudnnFlag = CUDNN_TENSOR_NHWC;
#else
  static const int kCudnnFlag = -1;
#endif
};

/*! \brief default layout for 4d tensor */
const int default_layout = kNCHW;

template<>
struct LayoutType<kNCDHW> {
  static const index_t kNdim = 5;
#if (MSHADOW_USE_CUDA && MSHADOW_USE_CUDNN == 1 && CUDNN_MAJOR >= 4)
  static const cudnnTensorFormat_t kCudnnFlag = CUDNN_TENSOR_NCHW;
#else
  static const int kCudnnFlag = -1;
#endif
};

template<>
struct LayoutType<kNDHWC> {
  static const index_t kNdim = 5;
#if (MSHADOW_USE_CUDA && MSHADOW_USE_CUDNN == 1 && CUDNN_MAJOR >= 4)
  static const cudnnTensorFormat_t kCudnnFlag = CUDNN_TENSOR_NHWC;
#else
  static const int kCudnnFlag = -1;
#endif
};

/*! \brief default layout for 5d tensor */
const int default_layout_5d = kNCDHW;

/*! \brief namespace for operators */
namespace op {
// binary operator
/*! \brief mul operator */
struct mul{
  /*! \brief map a, b to result using defined operation */
  template<typename DType>
  MSHADOW_XINLINE static DType Map(DType a, DType b) {
    return a * b;
  }
};
/*! \brief plus operator */
struct plus {
  /*! \brief map a, b to result using defined operation */
  template<typename DType>
  MSHADOW_XINLINE static DType Map(DType a, DType b) {
    return a + b;
  }
};
/*! \brief minus operator */
struct minus {
  /*! \brief map a, b to result using defined operation */
  template<typename DType>
  MSHADOW_XINLINE static DType Map(DType a, DType b) {
    return a - b;
  }
};
/*! \brief divide operator */
struct div {
  /*! \brief map a, b to result using defined operation */
  template<typename DType>
  MSHADOW_XINLINE static DType Map(DType a, DType b) {
    return a / b;
  }
};
/*! \brief get rhs */
struct right {
  /*! \brief map a, b to result using defined operation */
  template<typename DType>
  MSHADOW_XINLINE static DType Map(DType a, DType b) {
    return b;
  }
};
// unary operator/ function: example
// these operators can be defined by user,
// in the same style as binary and unary operator
// to use, simply write F<op::identity>( src )
/*! \brief identity function that maps a real number to it self */
struct identity{
  /*! \brief map a to result using defined operation */
  template<typename DType>
  MSHADOW_XINLINE static DType Map(DType a) {
    return a;
  }
};
}  // namespace op
/*! \brief namespace for savers */
namespace sv {
/*! \brief save to saver: = */
struct saveto {
  /*! \brief save b to a using save method */
  template<typename DType>
  MSHADOW_XINLINE static void Save(DType &a, DType b) { // NOLINT(*)
    a = b;
  }
  /*! \brief helper constant to use BLAS, alpha */
  inline static default_real_t AlphaBLAS(void) { return 1.0f; }
  /*! \brief helper constant to use BLAS, beta */
  inline static default_real_t BetaBLAS(void) { return 0.0f; }
  /*! \brief corresponding binary operator type */
  typedef op::right OPType;
};
/*! \brief save to saver: += */
struct plusto {
  /*! \brief save b to a using save method */
  template<typename DType>
  MSHADOW_XINLINE static void Save(DType &a, DType b) { // NOLINT(*)
    a += b;
  }
  /*! \brief helper constant to use BLAS, alpha */
  inline static default_real_t AlphaBLAS(void) { return 1.0f; }
  /*! \brief helper constant to use BLAS, beta */
  inline static default_real_t BetaBLAS(void) { return 1.0f; }
  /*! \brief corresponding binary operator type */
  typedef op::plus OPType;
};
/*! \brief minus to saver: -= */
struct minusto {
  /*! \brief save b to a using save method */
  template<typename DType>
  MSHADOW_XINLINE static void Save(DType &a, DType b) { // NOLINT(*)
    a -= b;
  }
  /*! \brief helper constant to use BLAS, alpha */
  inline static default_real_t AlphaBLAS(void) { return -1.0f; }
  /*! \brief helper constant to use BLAS, beta */
  inline static default_real_t BetaBLAS(void) { return 1.0f; }
  /*! \brief corresponding binary operator type */
  typedef op::minus OPType;
};
/*! \brief multiply to saver: *= */
struct multo {
  /*! \brief save b to a using save method */
  template<typename DType>
  MSHADOW_XINLINE static void Save(DType &a, DType b) { // NOLINT(*)
    a *= b;
  }
  /*! \brief corresponding binary operator type */
  typedef op::mul OPType;
};
/*! \brief divide to saver: /= */
struct divto {
  /*! \brief save b to a using save method */
  template<typename DType>
  MSHADOW_XINLINE static void Save(DType& a, DType b) { // NOLINT(*)
    a /= b;
  }
  /*! \brief corresponding binary operator type */
  typedef op::div OPType;
};
}  // namespace sv
/*! \brief namespace for potential reducer operations */
namespace red {
namespace limits {
/*!
 * \brief minimum value of certain types
 * \tparam DType data type
 */
template<typename DType>
MSHADOW_XINLINE DType MinValue(void);
/*! \brief minimum value of float */
template<>
MSHADOW_XINLINE float MinValue<float>(void) {
  return -FLT_MAX;
}
/*! \brief minimum value of double */
template<>
MSHADOW_XINLINE double MinValue<double>(void) {
  return -DBL_MAX;
}
/*! \brief minimum value of half */
template<>
MSHADOW_XINLINE half::half_t MinValue<half::half_t>(void) {
  return MSHADOW_HALF_MIN;
}
/*! \brief minimum value of uint8_t */
template<>
MSHADOW_XINLINE uint8_t MinValue<uint8_t>(void) {
  return 0;
}
/*! \brief minimum value of int8_t */
template<>
MSHADOW_XINLINE int8_t MinValue<int8_t>(void) {
  return SCHAR_MIN;
}
/*! \brief minimum value of int32_t */
template<>
MSHADOW_XINLINE int MinValue<int32_t>(void) {
  return INT_MIN;
}

/*!
 * \brief maximum value of certain types
 * \tparam DType data type
 */
template<typename DType>
MSHADOW_XINLINE DType MaxValue(void);
/*! \brief maximum value of float */
template<>
MSHADOW_XINLINE float MaxValue<float>(void) {
  return FLT_MAX;
}
/*! \brief maximum value of double */
template<>
MSHADOW_XINLINE double MaxValue<double>(void) {
  return DBL_MAX;
}
/*! \brief maximum value of uint8_t */
template<>
MSHADOW_XINLINE uint8_t MaxValue<uint8_t>(void) {
  return UCHAR_MAX;
}
/*! \brief maximum value of int8_t */
template<>
MSHADOW_XINLINE int8_t MaxValue<int8_t>(void) {
  return SCHAR_MAX;
}
/*! \brief maximum value of int32_t */
template<>
MSHADOW_XINLINE int MaxValue<int32_t>(void) {
  return INT_MAX;
}
}  // namespace limits

/*! \brief sum reducer */
struct sum {
  /*! \brief do reduction into dst */
  template<typename DType>
  MSHADOW_XINLINE static void Reduce(volatile DType& dst,  volatile DType src) { // NOLINT(*)
    dst += src;
  }
  /*!
   *\brief calculate gradient of redres with respect to redsrc,
   * redres: reduced result, redsrc: one of reduction element
   */
  template<typename DType>
  MSHADOW_XINLINE static DType PartialGrad(DType redres, DType redsrc) {
    return 1;
  }
  /*!
   *\brief set the initial value during reduction
   */
  template<typename DType>
  MSHADOW_XINLINE static void SetInitValue(DType &initv) { // NOLINT(*)
    initv = 0;
  }
};
/*! \brief maximum reducer */
struct maximum {
  /*! \brief do reduction into dst */
  template<typename DType>
  MSHADOW_XINLINE static void Reduce(volatile DType& dst,  volatile DType src) { // NOLINT(*)
    using namespace std;
#ifdef __CUDACC__
    dst = ::max(dst, src);
#else
    dst = max(dst, src);
#endif  // __CUDACC__
  }
  /*!
   * \brief calculate gradient of redres with respect to redsrc,
   * redres: reduced result, redsrc: one of reduction element
   */
  template<typename DType>
  MSHADOW_XINLINE static DType PartialGrad(DType redres, DType redsrc) {
    return redres == redsrc ? 1: 0;
  }
  /*!
   *\brief set the initial value during reduction
   */
  template<typename DType>
  MSHADOW_XINLINE static void SetInitValue(DType &initv) { // NOLINT(*)
    initv = limits::MinValue<DType>();
  }
};
/*! \brief minimum reducer */
struct minimum {
  /*! \brief do reduction into dst */
  template<typename DType>
  MSHADOW_XINLINE static void Reduce(volatile DType& dst,  volatile DType src) { // NOLINT(*)
    using namespace std;
#ifdef __CUDACC__
    dst = ::min(dst, src);
#else
    dst = min(dst, src);
#endif  // __CUDACC__
  }
  /*!
   * \brief calculate gradient of redres with respect to redsrc,
   * redres: reduced result, redsrc: one of reduction element
   */
  template<typename DType>
  MSHADOW_XINLINE static DType PartialGrad(DType redres, DType redsrc) {
    return redres == redsrc ? 1: 0;
  }
  /*!
   *\brief set the initial value during reduction
   */
  template<typename DType>
  MSHADOW_XINLINE static void SetInitValue(DType &initv) { // NOLINT(*)
    initv = -limits::MinValue<DType>();
  }
};
}  // namespace red

#define MSHADOW_TYPE_SWITCH(type, DType, ...)       \
  switch (type) {                                   \
  case mshadow::kFloat32:                           \
    {                                               \
      typedef float DType;                          \
      {__VA_ARGS__}                                 \
    }                                               \
    break;                                          \
  case mshadow::kFloat64:                           \
    {                                               \
      typedef double DType;                         \
      {__VA_ARGS__}                                 \
    }                                               \
    break;                                          \
  case mshadow::kFloat16:                           \
    {                                               \
      typedef mshadow::half::half_t DType;          \
      {__VA_ARGS__}                                 \
    }                                               \
    break;                                          \
  case mshadow::kUint8:                             \
    {                                               \
      typedef uint8_t DType;                        \
      {__VA_ARGS__}                                 \
    }                                               \
    break;                                          \
  case mshadow::kInt8:                              \
    {                                               \
      typedef int8_t DType;                         \
      {__VA_ARGS__}                                 \
    }                                               \
    break;                                          \
  case mshadow::kInt32:                             \
    {                                               \
      typedef int32_t DType;                        \
      {__VA_ARGS__}                                 \
    }                                               \
    break;                                          \
  default:                                          \
    LOG(FATAL) << "Unknown type enum " << type;     \
  }

#define MSHADOW_TYPE_SWITCH_WITH_HALF2(type, DType, ...)  \
  switch (type) {                                         \
  case mshadow::kFloat32:                                 \
    {                                                     \
      typedef float DType;                                \
      {__VA_ARGS__}                                       \
    }                                                     \
    break;                                                \
  case mshadow::kFloat64:                                 \
    {                                                     \
      typedef double DType;                               \
      {__VA_ARGS__}                                       \
    }                                                     \
    break;                                                \
  case mshadow::kFloat16:                                 \
    {                                                     \
      typedef mshadow::half::half2_t DType;               \
      {__VA_ARGS__}                                       \
    }                                                     \
    break;                                                \
  case mshadow::kUint8:                                   \
    {                                                     \
      typedef uint8_t DType;                              \
      {__VA_ARGS__}                                       \
    }                                                     \
    break;                                                \
  case mshadow::kInt32:                                   \
    {                                                     \
      typedef int32_t DType;                              \
      {__VA_ARGS__}                                       \
    }                                                     \
    break;                                                \
  default:                                                \
    LOG(FATAL) << "Unknown type enum " << type;           \
  }

#define MSHADOW_REAL_TYPE_SWITCH(type, DType, ...)  \
  switch (type) {                                   \
  case mshadow::kFloat32:                           \
    {                                               \
      typedef float DType;                          \
      {__VA_ARGS__}                                 \
    }                                               \
    break;                                          \
  case mshadow::kFloat64:                           \
    {                                               \
      typedef double DType;                         \
      {__VA_ARGS__}                                 \
    }                                               \
    break;                                          \
  case mshadow::kFloat16:                           \
    {                                               \
      typedef mshadow::half::half_t DType;          \
      {__VA_ARGS__}                                 \
    }                                               \
    break;                                          \
  case mshadow::kUint8:                             \
    LOG(FATAL) << "This operation only support "    \
                  "floating point types not uint8"; \
    break;                                          \
  case mshadow::kInt8:                              \
    LOG(FATAL) << "This operation only support "    \
                  "floating point types not int8";  \
    break;                                          \
  case mshadow::kInt32:                             \
    LOG(FATAL) << "This operation only support "    \
                  "floating point types, not int32";\
    break;                                          \
  default:                                          \
    LOG(FATAL) << "Unknown type enum " << type;     \
  }

#define MSHADOW_REAL_TYPE_SWITCH_EX(type$, DType$, DLargeType$, ...)  \
  switch (type$) {                                  \
  case mshadow::kFloat32:                           \
    {                                               \
      typedef float DType$;                         \
      typedef float DLargeType$;                    \
      {__VA_ARGS__}                                 \
    }                                               \
    break;                                          \
  case mshadow::kFloat64:                           \
    {                                               \
      typedef double DType$;                        \
      typedef double DLargeType$;                   \
      {__VA_ARGS__}                                 \
    }                                               \
    break;                                          \
  case mshadow::kFloat16:                           \
    {                                               \
      typedef mshadow::half::half_t DType$;         \
      typedef float DLargeType$;                    \
      {__VA_ARGS__}                                 \
    }                                               \
    break;                                          \
  case mshadow::kUint8:                             \
    LOG(FATAL) << "This operation only support "    \
                  "floating point types not uint8"; \
    break;                                          \
  case mshadow::kInt8:                              \
    LOG(FATAL) << "This operation only support "    \
                  "floating point types not int8";  \
    break;                                          \
  case mshadow::kInt32:                             \
    LOG(FATAL) << "This operation only support "    \
                  "floating point types, not int32";\
    break;                                          \
  default:                                          \
    LOG(FATAL) << "Unknown type enum " << type$;    \
  }

#define MSHADOW_LAYOUT_SWITCH(layout, Layout, ...)  \
  switch (layout) {                                 \
  case mshadow::kNCHW:                              \
    {                                               \
      const int Layout = kNCHW;                     \
      {__VA_ARGS__}                                 \
    }                                               \
    break;                                          \
  case mshadow::kNHWC:                              \
    {                                               \
      const int Layout = kNHWC;                     \
      {__VA_ARGS__}                                 \
    }                                               \
    break;                                          \
  case mshadow::kNCDHW:                             \
    {                                               \
      const int Layout = kNCDHW;                    \
      {__VA_ARGS__}                                 \
    }                                               \
    break;                                          \
  case mshadow::kNDHWC:                             \
    {                                               \
      const int Layout = kNDHWC;                    \
      {__VA_ARGS__}                                 \
    }                                               \
    break;                                          \
  default:                                          \
    LOG(FATAL) << "Unknown layout enum " << layout; \
  }

/*! \brief get data type size from type enum */
inline size_t mshadow_sizeof(int type) {
  int size = 0;
  MSHADOW_TYPE_SWITCH(type, DType, size = sizeof(DType););
  return size;
}

}  // namespace mshadow
#endif  // MSHADOW_BASE_H_

max_stars_repo_name: withinnoitatpmet/mxnetsource

max_issues_repo_name: withinnoitatpmet/mxnetsource

max_forks_repo_name: withinnoitatpmet/mxnetsource


text:
#include <mpi.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_integration.h>

#include "allvars.h"
#include "proto.h"

static double logTimeBegin;
static double logTimeMax;

double drift_integ(double a, void *param)
{
  double h;

  h = hubble_function(a);

  return 1 / (h * a * a * a);
}

double gravkick_integ(double a, void *param)
{
  double h;

  h = hubble_function(a);

  return 1 / (h * a * a);
}

double hydrokick_integ(double a, void *param)
{
  double h;

  h = hubble_function(a);

  return 1 / (h * pow(a, 3 * GAMMA_MINUS1) * a);
}

double growthfactor_integ(double a, void *param)
{
  double s;

  s = hubble_function(a) / All.Hubble * sqrt(a * a * a);

  return pow(sqrt(a) / s, 3);
}

void init_drift_table(void)
{
#define WORKSIZE 100000
  int i;
  double result, abserr;

  /*---------- DEBUG ! ------------*/
#ifdef DARKENERGY_DEBUG
  FILE *FdDKfac;
  char mode[2], buf[200];

  strcpy(mode, "w");
  sprintf(buf, "%s%s", All.OutputDir, "driftkickfac.txt");
  FdDKfac = fopen(buf, mode);
  fprintf(FdDKfac, "i a drift GravKick HydroKick\n");
  /*---------- DEBUG ! ------------*/
#endif

  gsl_function F;
  gsl_integration_workspace *workspace;

  logTimeBegin = log(All.TimeBegin);
  logTimeMax = log(All.TimeMax);

  workspace = gsl_integration_workspace_alloc(WORKSIZE);

  for(i = 0; i < DRIFT_TABLE_LENGTH; i++)
    {
      F.function = &drift_integ;
      gsl_integration_qag(&F, exp(logTimeBegin), exp(logTimeBegin + ((logTimeMax - logTimeBegin) / DRIFT_TABLE_LENGTH) * (i + 1)), All.Hubble,	/* note: absolute error just a dummy */
			  1.0e-8, WORKSIZE, GSL_INTEG_GAUSS41, workspace, &result, &abserr);
      DriftTable[i] = result;

      F.function = &gravkick_integ;
      gsl_integration_qag(&F, exp(logTimeBegin), exp(logTimeBegin + ((logTimeMax - logTimeBegin) / DRIFT_TABLE_LENGTH) * (i + 1)), All.Hubble,	/* note: absolute error just a dummy */
			  1.0e-8, WORKSIZE, GSL_INTEG_GAUSS41, workspace, &result, &abserr);
      GravKickTable[i] = result;

      F.function = &hydrokick_integ;
      gsl_integration_qag(&F, exp(logTimeBegin), exp(logTimeBegin + ((logTimeMax - logTimeBegin) / DRIFT_TABLE_LENGTH) * (i + 1)), All.Hubble,	/* note: absolute error just a dummy */
			  1.0e-8, WORKSIZE, GSL_INTEG_GAUSS41, workspace, &result, &abserr);
      HydroKickTable[i] = result;

#ifdef DARKENERGY_DEBUG
  /*---------- DEBUG ! ------------*/
      fprintf(FdDKfac, "%d %e %e %e %e \n", i,
	      exp(logTimeBegin + ((logTimeMax - logTimeBegin) / DRIFT_TABLE_LENGTH) * (i + 1)), DriftTable[i],
	      GravKickTable[i], HydroKickTable[i]);
  /*---------- DEBUG ! ------------*/
#endif
    }

  gsl_integration_workspace_free(workspace);

#ifdef DARKENERGY_DEBUG
  /*---------- DEBUG ! ------------*/
  fclose(FdDKfac);
  /*---------- DEBUG ! ------------*/
#endif
}

/*! This function integrates the cosmological prefactor for a drift
 *   step between time0 and time1. The value returned is
 *  \f[ \int_{a_0}^{a_1} \frac{{\rm d}a}{H(a)}
 *  \f]
 *
 *  A lookup-table is used for reasons of speed.
 */
double get_drift_factor(int time0, int time1)
{
  double a1, a2, df1, df2, u1, u2;
  int i1, i2;
  static int last_time0 = -1, last_time1 = -1;
  static double last_value;

  if(time0 == last_time0 && time1 == last_time1)
    return last_value;

  /* note: will only be called for cosmological integration */

  a1 = logTimeBegin + time0 * All.Timebase_interval;
  a2 = logTimeBegin + time1 * All.Timebase_interval;

  if(logTimeMax > logTimeBegin)
    u1 = (a1 - logTimeBegin) / (logTimeMax - logTimeBegin) * DRIFT_TABLE_LENGTH;
  else
    u1 = 0;
  i1 = (int) u1;
  if(i1 >= DRIFT_TABLE_LENGTH)
    i1 = DRIFT_TABLE_LENGTH - 1;

  if(i1 <= 1)
    df1 = u1 * DriftTable[0];
  else
    df1 = DriftTable[i1 - 1] + (DriftTable[i1] - DriftTable[i1 - 1]) * (u1 - i1);

  if(logTimeMax > logTimeBegin)
    u2 = (a2 - logTimeBegin) / (logTimeMax - logTimeBegin) * DRIFT_TABLE_LENGTH;
  else
    u2 = 0;
  i2 = (int) u2;
  if(i2 >= DRIFT_TABLE_LENGTH)
    i2 = DRIFT_TABLE_LENGTH - 1;

  if(i2 <= 1)
    df2 = u2 * DriftTable[0];
  else
    df2 = DriftTable[i2 - 1] + (DriftTable[i2] - DriftTable[i2 - 1]) * (u2 - i2);

  last_time0 = time0;
  last_time1 = time1;

  return last_value = (df2 - df1);
}

double get_gravkick_factor(int time0, int time1)
{
  double a1, a2, df1, df2, u1, u2;
  int i1, i2;
  static int last_time0 = -1, last_time1 = -1;
  static double last_value;

  if(time0 == last_time0 && time1 == last_time1)
    return last_value;

  /* note: will only be called for cosmological integration */

  a1 = logTimeBegin + time0 * All.Timebase_interval;
  a2 = logTimeBegin + time1 * All.Timebase_interval;

  if(logTimeMax > logTimeBegin)
    u1 = (a1 - logTimeBegin) / (logTimeMax - logTimeBegin) * DRIFT_TABLE_LENGTH;
  else
    u1 = 0;
  i1 = (int) u1;
  if(i1 >= DRIFT_TABLE_LENGTH)
    i1 = DRIFT_TABLE_LENGTH - 1;

  if(i1 <= 1)
    df1 = u1 * GravKickTable[0];
  else
    df1 = GravKickTable[i1 - 1] + (GravKickTable[i1] - GravKickTable[i1 - 1]) * (u1 - i1);

  if(logTimeMax > logTimeBegin)
    u2 = (a2 - logTimeBegin) / (logTimeMax - logTimeBegin) * DRIFT_TABLE_LENGTH;
  else
    u2 = 0;
  i2 = (int) u2;
  if(i2 >= DRIFT_TABLE_LENGTH)
    i2 = DRIFT_TABLE_LENGTH - 1;

  if(i2 <= 1)
    df2 = u2 * GravKickTable[0];
  else
    df2 = GravKickTable[i2 - 1] + (GravKickTable[i2] - GravKickTable[i2 - 1]) * (u2 - i2);

  last_time0 = time0;
  last_time1 = time1;

  return last_value = (df2 - df1);
}

double get_hydrokick_factor(int time0, int time1)
{
  double a1, a2, df1, df2, u1, u2;
  int i1, i2;
  static int last_time0 = -1, last_time1 = -1;
  static double last_value;

  if(time0 == last_time0 && time1 == last_time1)
    return last_value;

  /* note: will only be called for cosmological integration */

  a1 = logTimeBegin + time0 * All.Timebase_interval;
  a2 = logTimeBegin + time1 * All.Timebase_interval;

  if(logTimeMax > logTimeBegin)
    u1 = (a1 - logTimeBegin) / (logTimeMax - logTimeBegin) * DRIFT_TABLE_LENGTH;
  else
    u1 = 0;
  i1 = (int) u1;
  if(i1 >= DRIFT_TABLE_LENGTH)
    i1 = DRIFT_TABLE_LENGTH - 1;

  if(i1 <= 1)
    df1 = u1 * HydroKickTable[0];
  else
    df1 = HydroKickTable[i1 - 1] + (HydroKickTable[i1] - HydroKickTable[i1 - 1]) * (u1 - i1);

  if(logTimeMax > logTimeBegin)
    u2 = (a2 - logTimeBegin) / (logTimeMax - logTimeBegin) * DRIFT_TABLE_LENGTH;
  else
    u2 = 0;
  i2 = (int) u2;
  if(i2 >= DRIFT_TABLE_LENGTH)
    i2 = DRIFT_TABLE_LENGTH - 1;

  if(i2 <= 1)
    df2 = u2 * HydroKickTable[0];
  else
    df2 = HydroKickTable[i2 - 1] + (HydroKickTable[i2] - HydroKickTable[i2 - 1]) * (u2 - i2);

  last_time0 = time0;
  last_time1 = time1;

  return last_value = (df2 - df1);
}

ext: c

max_stars_repo_name: egpbos/egp

max_issues_repo_name: egpbos/egp

max_forks_repo_name: egpbos/egp


text:
#pragma once

#include <list>
#include <map>
#include <memory>
#include <string>
#include <gsl/gsl_math.h>

#include "Expressions/Expression.h"
#include "Scanner/scanner.h"
#include "Parser/parser.h"
#include "Utils/Exception.h"

class MathEngine {

    std::unique_ptr<Parser> parser;
    Variables variables;

    public:
        MathEngine();

        std::list<Scanner::Token> scan(const std::string& input);

        expression parse(const std::string& input);
        expression operator()(const std::string& input);

        expression eval(const std::string& input);
        expression eval(expression inputExpr);

        std::string formatInput(const std::string& input, int& cursorPosition);
        std::string evaluateOutput(const std::string& input, const std::string& output);

};

ext: h

max_stars_repo_name: antoniojkim/CalcPlusPlus

max_issues_repo_name: antoniojkim/CalcPlusPlus

max_forks_repo_name: antoniojkim/CalcPlusPlus


text:
/**
 * \author Sylvain Marsat, University of Maryland - NASA GSFC
 *
 * \brief C header for EOBNRv2HM reduced order model (non-spinning version).
 * See CQG 31 195010, 2014, arXiv:1402.4146 for details on the reduced order method.
 * See arXiv:1106.1021 for the EOBNRv2HM model.
 *
 * Borrows from the SEOBNR ROM LAL code written by Michael Puerrer and John Veitch.
 *
 * Put the untared data in the directory designated by the environment variable ROM_DATA_PATH.
 *
 * Parameter range:
 *   q = 1-12 (almost)
 *   No spin
 *   Mtot >= 10Msun for fstart=8Hz
 *
 */

#ifndef _EOBNRV2HMROM_H
#define _EOBNRV2HMROM_H

#define _XOPEN_SOURCE 500

#ifdef __GNUC__
#define UNUSED __attribute__ ((unused))
#else
#define UNUSED
#endif

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <complex.h>
#include <time.h>
#include <unistd.h>
#include <getopt.h>
#include <stdbool.h>
#include <string.h>

#include <gsl/gsl_errno.h>
#include <gsl/gsl_bspline.h>
#include <gsl/gsl_blas.h>
#include <gsl/gsl_min.h>
#include <gsl/gsl_spline.h>
#include <gsl/gsl_complex.h>

#include "constants.h"
#include "struct.h"
#include "EOBNRv2HMROMstruct.h"

#if defined(__cplusplus)
extern "C" {
#elif 0
} /* so that editors will match preceding brace */
#endif

   /********General definitions********/

#define nk_amp 10  /* number of SVD-modes == number of basis functions for amplitude */
#define nk_phi 20  /* number of SVD-modes == number of basis functions for phase */

/********External array for the list of modes********/
#define nbmodemax 5
extern const int listmode[nbmodemax][2];

/**************************************************/
/**************** Prototypes **********************/

/* Functions to load, initalize and cleanup data */
int EOBNRv2HMROM_Init_DATA(void);
int EOBNRv2HMROM_Init(const char dir[]);

void EOBNRHMROMdata_Init(EOBNRHMROMdata **data);
void EOBNRHMROMdata_interp_Init(EOBNRHMROMdata_interp **data_interp);
void EOBNRHMROMdata_coeff_Init(EOBNRHMROMdata_coeff **data_coeff);

void EOBNRHMROMdata_Cleanup(EOBNRHMROMdata *data);
void EOBNRHMROMdata_interp_Cleanup(EOBNRHMROMdata_interp *data_interp);
void EOBNRHMROMdata_coeff_Cleanup(EOBNRHMROMdata_coeff *data_coeff);

/* Function to read data */
int Read_Data_Mode(const char dir[], const int mode[2], EOBNRHMROMdata *data);

/* Functions to interpolate the data and to evaluate the interpolated data for a given q */

int Evaluate_Spline_Data(
  const double q,                           /* Input: q-value for which projection coefficients should be evaluated */
  const EOBNRHMROMdata_interp* data_interp,  /* Input: data in interpolated form */
  EOBNRHMROMdata_coeff* data_coeff           /* Output: vectors of projection coefficients and shifts in time and phase */
);

int Interpolate_Spline_Data(
  const EOBNRHMROMdata *data,           /* Input: data in vector/matrix form to interpolate */
  EOBNRHMROMdata_interp *data_interp    /* Output: interpolated data */
);

/* Functions for waveform reconstruction */

int EOBNRv2HMROMCore(
  struct tagListmodesCAmpPhaseFrequencySeries **listhlm,
  int nbmode,
  double tRef,
  double phiRef,
  double fRef,
  double Mtot_sec,
  double q,
  double distance,
  int setphiRefattRef);

int SimEOBNRv2HMROM(
  ListmodesCAmpPhaseFrequencySeries **listhlm,  /* Output: list of modes in Frequency-domain amplitude and phase form */
  int nbmode,                                    /* Number of modes to generate (starting with the 22) */
  double tRef,                                   /* Time shift with respect to the 22-fit removed waveform (s) */
  double phiRef,                                 /* Phase at reference frequency */
  double fRef,                                   /* Reference frequency (Hz); 0 defaults to fLow */
  double m1SI,                                   /* Mass of companion 1 (kg) */
  double m2SI,                                   /* Mass of companion 2 (kg) */
  double distance,                               /* Distance of source (m) */
  int setphiRefattRef);                          /* Flag for adjusting the FD phase at phiRef at the given fRef, which depends also on tRef - if false, treat phiRef simply as an orbital phase shift (minus an observer phase shift) */

int SimEOBNRv2HMROMExtTF2(
  ListmodesCAmpPhaseFrequencySeries **listhlm,   /* Output: list of modes in Frequency-domain amplitude and phase form */
  int nbmode,                                    /* Number of modes to generate (starting with the 22) */
  double Mf_match,                               /* Minimum frequency using EOBNRv2HMROM in inverse total mass units*/
  double minf,                                   /* Minimum frequency required */
  int tagexthm,                                  /* Tag to decide whether or not to extend the higher modes as well */
  double deltatRef,                              /* Time shift so that the peak of the 22 mode occurs at deltatRef */
  double phiRef,                                 /* Phase at reference frequency */
  double fRef,                                   /* Reference frequency (Hz); 0 defaults to fLow */
  double m1SI,                                   /* Mass of companion 1 (kg) */
  double m2SI,                                   /* Mass of companion 2 (kg) */
  double distance,                               /* Distance of source (m) */
  int setphiRefattRef);                          /* Flag for adjusting the FD phase at phiRef at the given fRef, which depends also on tRef - if false, treat phiRef simply as an orbital phase shift (minus an observer phase shift) */

#if 0
{ /* so that editors will match succeeding brace */
#elif defined(__cplusplus)
}
#endif

#endif /* _EOBNRV2HMROM_H */

max_stars_repo_name: JohnGBaker/flare

max_issues_repo_name: JohnGBaker/flare

max_forks_repo_name: JohnGBaker/flare


text:
/*
 * Copyright (c) 2016-2021 lymastee, All rights reserved.
 * Contact: lymastee@hotmail.com
 *
 * This file is part of the gslib project.
 *
 * Permission is hereby granted, free of charge, to any person obtaining a copy
 * of this software and associated documentation files (the "Software"), to deal
 * in the Software without restriction, including without limitation the rights
 * to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
 * copies of the Software, and to permit persons to whom the Software is
 * furnished to do so, subject to the following conditions:
 *
 * The above copyright notice and this permission notice shall be included in all
 * copies or substantial portions of the Software.
 *
 * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
 * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
 * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
 * AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
 * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
 * OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
 * SOFTWARE.
 */

#pragma once

#ifndef graph_46d3170b_a19a_4663_bbf9_c7372932671b_h
#define graph_46d3170b_a19a_4663_bbf9_c7372932671b_h

#include <assert.h>
#include <gslib/std.h>

__gslib_begin__

template<class _meta>
struct _mono_rel
{
    _meta*          _next;

    _mono_rel() { _next = nullptr; }
    _meta* next() const { return _next; }
    void set_next(_meta* ptr) { _next = ptr; }
    _meta* operator+(int d) const
    {
        assert(d >= 0);
        _meta* p = static_cast<_meta*>(this);
        for( ; d > 0; d --)
            p = p->_next;
        return p;
    }
    int operator-(const _meta& that) const
    {
        int c = 0;
        _meta* p = static_cast<_meta*>(this);
        for( ; p->_next; p = p->_next, c ++) {
            if(p == &that)
                return c;
        }
        assert(!"unexpected.");
        return -1;
    }
    int count_to_tail() const
    {
        int c = 0;
        for(_meta* p = static_cast<_meta*>(this); p->_next; p = p->_next, c ++);
        return c;
    }
    static void shake(_meta* p1, _meta* p2) { p1->_next = p2; }
};

template<class _meta>
struct _dual_rel
{
    _meta           *_prev, *_next;

    _dual_rel() { _prev = _next = nullptr; }
    _meta* prev() const { return _prev; }
    _meta* next() const { return _next; }
    void set_prev(_meta* ptr) { _prev = ptr; }
    void set_next(_meta* ptr) { _next = ptr; }
    _meta* operator+(int d) const
    {
        if(d < 0)
            return operator-(-d);
        _meta* p = static_cast<_meta*>(this);
        for( ; d > 0; d --)
            p = p->_next;
        return p;
    }
    _meta* operator-(int d) const
    {
        if(d < 0)
            return operator+(-d);
        _meta* p = static_cast<_meta*>(this);
        for( ; d > 0; d --)
            p = p->_prev;
        return p;
    }
    int operator-(const _meta& that) const
    {
        int c = 0;
        _meta* p = static_cast<_meta*>(this);
        for( ; p->_next; p = p->_next, c ++) {
            if(p == &that)
                return c;
        }
        assert(!"unexpected.");
        return -1;
    }
    int count_to_head() const
    {
        int c = 0;
        for(_meta* p = static_cast<_meta*>(this); p->_next; p = p->_prev, c ++);
        return c;
    }
    int count_to_tail() const
    {
        int c = 0;
        for(_meta* p = static_cast<_meta*>(this); p->_next; p = p->_next, c ++);
        return c;
    }
    static void shake(_meta* p1, _meta* p2) { p1->_next = p2, p2->_prev = p1; }
};

template<class _meta>
struct _pack_rel
{
    _meta       *_from, *_to;
    int         _size;

    _pack_rel() { _from = _to = nullptr, _size = 0; }
    _pack_rel(_meta* from, _meta* to) { _from = from, _to = to, _size = 0; }
    _meta* from() const { return _from; }
    _meta* to() const { return _to; }
    bool empty() const { return !_from; }
    int size() const { return _size; }
    void set_from(_meta* ptr) { _from = ptr; }
    void set_to(_meta* ptr) { _to = ptr; }
    void increase_size() { _size ++; }
    void decrease_size() { _size --; }
    void set_size(int s) { _size = s; }
    bool check_mono() const
    {
        if(!_from)
            return false;
        _meta* p = static_cast<_meta*>(_from);
        for( ; p->_next; p = p->_next);
        return p == _to;
    }
    bool check_dual() const
    {
        if(!_from || _from->_prev)
            return false;
        _meta* p = static_cast<_meta*>(_from);
        for( ; p->_next; p = p->_next) {
            if(p->_prev && p != p->_prev->_next)
                return false;
        }
        return p == _to;
    }
    bool is_belong_to(const _meta* ptr) const
    {
        _meta* p = static_cast<_meta*>(_from);
        for( ; p; p = p->_next) {
            if(ptr == p)
                return true;
        }
        return false;
    }
};

class _graph_ortho_rel:
    public _dual_rel<_graph_ortho_rel>
{
public:
    typedef _pack_rel<_graph_ortho_rel> hierarchic_rel;
    typedef _dual_rel<_graph_ortho_rel> peer_rel;
    typedef hierarchic_rel upper_rel;
    typedef hierarchic_rel lower_rel;

protected:
    hierarchic_rel  _parents, _children;

public:
    hierarchic_rel& parents() { return _parents; }
    hierarchic_rel& children() { return _children; }
    peer_rel& siblings() { return static_cast<peer_rel&>(*this); }
    const hierarchic_rel& const_parents() const { return _parents; }
    const hierarchic_rel& const_children() const { return _children; }
    const peer_rel& const_siblings() const { return static_cast<const peer_rel&>(*this); }
    int parents_count() const { return _parents.size(); }
    int children_count() const { return _children.size(); }
};

struct _graph_trait_copy {};
struct _graph_trait_detach {};

template<class _ty,
    class _rel = _graph_ortho_rel
    >
struct _graph_ortho_cpy_wrapper:
    public _ty,
    public _rel
{
    typedef _ty value;
    typedef _rel relation;
    typedef _graph_ortho_cpy_wrapper<_ty, _rel> wrapper;
    typedef _graph_trait_copy tsf_behavior;

    value* get_ptr() { return static_cast<value*>(this); }
    const value* const_ptr() const { return static_cast<const value*>(this); }
    value& get_ref() { return static_cast<value&>(*this); }
    const value& const_ref() const { return static_cast<const value&>(*this); }
    template<class _cst>
    void born() {}  /* maybe a problem. */
    void born() {}
    void kill() {}
    void copy(const wrapper* a) { assign<value>(a); }
    void attach(wrapper* a) { !!error!! }
    int get_parents_count() const { return parents_count(); }
    int get_children_count() const { return children_count(); }
    template<class c>
    void assign(const wrapper* a) { static_cast<c&>(*this) = static_cast<c&>(*a); }
};

template<class _ty,
    class _rel = _graph_ortho_rel
    >
struct _graph_ortho_wrapper:
    public _rel
{
    typedef _ty value;
    typedef _rel relation;
    typedef _graph_ortho_wrapper<_ty, _rel> wrapper;
    typedef _graph_trait_detach tsf_behavior;

    value*          _value;

    wrapper() { _value = nullptr; }
    value* get_ptr() { return _value; }
    const value* const_ptr() const { return _value; }
    value& get_ref() { return *_value; }
    const value& const_ref() const { return *_value; }
    relation& get_relation() { return static_cast<relation&>(*this); }
    const relation& const_relation() const { return static_cast<const relation&>(*this); }
    template<class _cst>
    void born() { _value = new _cst; }
    void born() { _value = new value; }
    void kill() { delete _value; }
    void copy(const wrapper* a) { assert(!"prevent."); }
    void attach(wrapper* a)
    {
        assert(a && a->_value);
        kill();
        _value = a->_value;
        a->_value = nullptr;
    }
    int get_parents_count() const { return parents_count(); }
    int get_children_count() const { return children_count(); }
};

template<class _wrapper>
struct _graph_allocator
{
    typedef _wrapper wrapper;
    static wrapper* born() { return new wrapper; }
    static void kill(wrapper* w) { delete w; }
};

template<class _val>
struct _graph_ortho_val
{
    typedef _val value;
    typedef const _val const_value;
    union
    {
        value*          _vptr;
        const_value*    _cvptr;
    };
    value* get_wrapper() const { return _vptr; }
    value* trace_front(int d) const { return _vptr->siblings() + d; }
    value* trace_back(int d) const { return _vptr->siblings() - d; }
    value* trace_up(int d) const
    {
        if(d < 0)
            return trace_down(-d);
        value* p = _vptr;
        for( ; d > 0; d --) {
            assert(!p->parents().empty());
            p = p->parents().from();
        }
        return p;
    }
    value* trace_down(int d) const
    {
        if(d < 0)
            return trace_up(-d);
        value* p = _vptr;
        for( ; d > 0; d --) {
            assert(!p->children().empty());
            p = p->children().from();
        }
        return p;
    }
    value* trace_front() const
    {
        assert(_vptr);
        return _vptr->next();
    }
    value* trace_back() const
    {
        assert(_vptr);
        return _vptr->prev();
    }
};

template<class _ty,
    class _wrapper = _graph_ortho_cpy_wrapper<_ty>
    >
class _graph_ortho_const_iterator:
    public _graph_ortho_val<_wrapper>
{
public:
    typedef _ty value;
    typedef _wrapper wrapper;
    typedef _graph_ortho_const_iterator<_ty, _wrapper> iterator;

public:
    iterator(const wrapper* w = nullptr) { _cvptr = w; }
    bool is_valid() const { return _cvptr != nullptr; }
    const value* get_ptr() const { return _cvptr->const_ptr(); }
    const value* operator->() const { return _cvptr->const_ptr(); }
    const value& operator*() const { return _cvptr->const_ref(); }
    iterator operator+(int d) const { return iterator(trace_front(d)); }
    iterator operator-(int d) const { return iterator(trace_back(d)); }
    iterator operator<<(int d) const { return iterator(trace_up(d)); }
    iterator operator>>(int d) const { return iterator(trace_down(d)); }
    int operator-(iterator i) const { return _cvptr->peer_rel::operator-(*(i._cvptr)); }
};

template<class _ty,
    class _wrapper = _graph_ortho_cpy_wrapper<_ty>
    >
class _graph_ortho_iterator:
    public _graph_ortho_const_iterator<_ty, _wrapper>
{
public:
    typedef _ty value;
    typedef _wrapper wrapper;
    typedef _graph_ortho_const_iterator<_ty, _wrapper> superref;
    typedef _graph_ortho_const_iterator<_ty, _wrapper> const_iterator;
    typedef _graph_ortho_iterator<_ty, _wrapper> iterator;

public:
    iterator(wrapper* w = nullptr) { _vptr = w; }
    value* get_ptr() const { return _vptr->get_ptr(); }
    value* operator->() const { return _vptr->get_ptr(); }
    value& operator*() const { return _vptr->get_ref(); }
    bool operator==(const iterator& that) const { return _vptr == that._vptr; }
    bool operator!=(const iterator& that) const { return _vptr != that._vptr; }
    bool operator==(const const_iterator& that) const { return _vptr == that._vptr; }
    bool operator!=(const const_iterator& that) const { return _vptr != that._vptr; }
    void operator++() { _vptr = trace_front(); }
    void operator--() { _vptr = trace_back(); }
    void operator+=(int d) { _vptr = trace_front(d); }
    void operator-=(int d) { _vptr = trace_back(d); }
    void operator<<=(int d) { _vptr = trace_up(d); }
    void operator>>=(int d) { _vptr = trace_down(d); }
    operator const_iterator() { return const_iterator(_cvptr); }
};

template<class _ty,
    class _wrapper = _graph_ortho_wrapper<_ty>,
    class _alloc = _graph_allocator<_wrapper>
    >
class ortho_graph:
    public _graph_ortho_val<_wrapper>
{
public:
    typedef _ty value;
    typedef _wrapper wrapper;
    typedef _alloc alloc;
    typedef ortho_graph<_ty, _wrapper, _alloc> myref;
    typedef _graph_ortho_const_iterator const_iterator;
    typedef _graph_ortho_iterator iterator;

public:
    ortho_graph() { _vptr = nullptr; }
    ~ortho_graph() { destroy(); }
    void destroy()
    {
        for(wrapper* i = _vptr; i; ) {
            wrapper* n = i->next();
            _erase(i);
            i = n;
        }
    }
    void clear() { destroy(); }
    iterator get_root() { return iterator(_vptr); }
    const_iterator const_root() const { return const_iterator(_cvptr); }
    bool is_valid() const { return !_cvptr; }
    iterator insert_before(iterator i) { return insert_before<value>(i); }
    iterator insert_after(iterator i) { return insert_after<value>(i); }
    iterator create_child(iterator from, iterator to = from) { return create_child<value>(from, to); }
    iterator erase(iterator i)
    {
        iterator r(i->next());
        _erase(i.get_wrapper());
        return r;
    }

public:
    template<class _cst>
    iterator insert_before(iterator i)
    {
        if(!i.is_valid())
            return _cvptr ? get_root() : _init<_cst>();
        wrapper* n = alloc::born();
        n->born<_cst>();
        _fix_peer_rel(i->prev(), n, *i);
        _fix_upper_rel(n, *i);
        _fix_lower_rel(n, *i);
        /* fix root, keep the root always be the top left corner. */
        if(i.get_wrapper() == _vptr)
            _vptr = n;
        return iterator(n);
    }
    template<class _cst>
    iterator insert_after(iterator i)
    {
        if(!i.is_valid())
            return _cvptr ? get_root() : _init<_cst>();
        wrapper* n = alloc::born();
        n->born<_cst>();
        _fix_peer_rel(*i, n, i->next());
        _fix_upper_rel(n, *i);
        _fix_lower_rel(n, *i);
        return iterator(n);
    }
    template<class _cst>
    iterator create_child(iterator from, iterator to = from)
    {
        assert(from.is_valid() && from->children().empty());
        wrapper* n = alloc::born();
        n->born<_cst>();
        _fix_upper_rel(n, *from, *to);
        return iterator(n);
    }
    template<class _lambda>
    void _for_range(iterator from, iterator to, _lambda lam)
    {
        for(iterator i = from; ; i ++) {
            lam(i->get_wrapper());
            if(i == to)
                break;
        }
    }
    template<class _lambda>
    void for_range(iterator from, iterator to, _lambda lam)
    {
        for(iterator i = from; ; i ++) {
            lam(*i);
            if(i == to)
                break;
        }
    }

protected:
    template<class _cst>
    iterator _init()
    {
        _vptr = alloc::born();
        _vptr->born<_cst>();
        return iterator(_vptr);
    }
    void _erase(wrapper* w)
    {
        assert(w);
        _fix_upper_rel(w);
        _fix_lower_rel(w);
        peer_rel::shake(w->prev(), w->next());
        lower_rel& rel = w->children();
        if(!rel.empty())
            _for_range(rel.from(), rel.to(), [](wrapper* w) { _erase(w); });
        if(w == _vptr)
            _vptr = w->next();
        w->kill();
        alloc::kill(w);
    }
    void _set_range_upper(wrapper* from, wrapper* to, wrapper* r1, wrapper* r2, int size)
    {
        _for_range(iterator(from), iterator(to), [&r1, &r2, &size](wrapper* w) {
            upper_rel& rel = w->parents();
            rel.set_from(r1);
            rel.set_to(r2);
            rel.set_size(size);
        });
    }
    void _set_range_lower(wrapper* from, wrapper* to, wrapper* r1, wrapper* r2, int size)
    {
        _for_range(iterator(from), iterator(to), [&r1, &r2, &size](wrapper* w) {
            lower_rel& rel = w->children();
            rel.set_from(r1);
            rel.set_to(r2);
            rel.set_size(size);
        });
    }
    void _fix_peer_rel(wrapper* w1, wrapper* w2, wrapper* w3)
    {
        assert(w2 && (w1 || w3));
        if(w1 != nullptr)
            peer_rel::shake(w1, w2);
        if(w3 != nullptr)
            peer_rel::shake(w2, w3);
    }
    void _fix_upper_rel(wrapper* n, wrapper* w)
    {
        assert(n && w);
        upper_rel& urel = w->parents();
        wrapper* u1 = urel.from();
        if(u1 == nullptr)
            return;
        wrapper* u2 = urel.to();
        lower_rel& lrel = u1->children();
        wrapper* l1 = lrel.from();
        wrapper* l2 = lrel.to();
        assert(l1 && l2);
        wrapper *from = l1, *to = l2;
        if(l1 == w && n == w->prev())
            from = n;
        else if(l2 == w && n == w->next())
            to = n;
        _set_range_lower(u1, u2, from, to, lrel.size() + 1);
        upper_rel& u = n->parents();
        u.set_from(u1);
        u.set_to(u2);
    }
    void _fix_lower_rel(wrapper* n, wrapper* w)
    {
        assert(n && w);
        lower_rel& lrel = w->children();
        wrapper* l1 = lrel.from();
        if(l1 == nullptr)
            return;
        wrapper* l2 = lrel.to();
        upper_rel& urel = l1->parents();
        wrapper* u1 = urel.from();
        wrapper* u2 = urel.to();
        assert(u1 && u2);
        wrapper *from = u1, *to = u2;
        if(u1 == w && n == w->prev())
            from = n;
        else if(u2 == w && n == w->next())
            to = n;
        _set_range_upper(l1, l2, from, to, urel.size() + 1);
        lower_rel& l = n->children();
        l.set_from(l1);
        l.set_to(l2);
    }
    void _fix_upper_rel(wrapper* w)
    {
        assert(w);
        upper_rel& urel = w->parents();
        wrapper* u1 = urel.from();
        if(u1 == nullptr)
            return;
        wrapper* u2 = urel.to();
        assert(u2);
        lower_rel& lrel = u1->children();
        wrapper* l1 = lrel.from();
        wrapper* l2 = lrel.to();
        assert(l1 && l2);
        wrapper *from = l1, *to = l2;
        if(l1 == w && l2 == w)
            from = to = nullptr;
        else if(l1 == w)
            from = w->next();
        else if(l2 == w)
            to = w->prev();
        _set_range_lower(u1, u2, from, to, lrel.size() - 1);
        urel.set_from(nullptr);
        urel.set_to(nullptr);
        urel.set_size(0);
    }
    void _fix_lower_rel(wrapper* w)
    {
        assert(w);
        lower_rel& lrel = w->children();
        wrapper* l1 = lrel.from();
        if(l1 == nullptr)
            return;
        wrapper* l2 = lrel.to();
        assert(l2);
        upper_rel& urel = l1->parents();
        wrapper* u1 = urel.from();
        wrapper* u2 = urel.to();
        assert(u1 && u2);
        wrapper *from = u1, *to = u2;
        if(u1 == w && u2 == w);
        else if(u1 == w || u2 == w) {
            u1 == w ? from = w->next() : to = w->prev();
            lrel.set_from(nullptr);
            lrel.set_to(nullptr);
            lrel.set_size(0);
        }
        _set_range_upper(l1, l2, from, to, urel.size() - 1);
    }
    void _fix_upper_rel(wrapper* w, wrapper* u1, wrapper* u2)
    {
        assert(w && u1 && u2);
        upper_rel& urel = w->parents();
        urel.set_from(u1);
        urel.set_to(u2);
        urel.set_size(u1->count_to_tail(u2) + 1);
        _set_range_lower(u1, u2, w, w, 1);
    }
};

__gslib_end__

#endif

max_stars_repo_name: lymastee/gslib

max_issues_repo_name: lymastee/gslib

max_forks_repo_name: lymastee/gslib


text:
#include <math.h>
#include <stdlib.h>

#include "fd3sep.h"
/* #include "mxfuns.h" */

#include <gsl/gsl_linalg.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_fft_real.h>
#include <gsl/gsl_fft_halfcomplex.h>

#define SVCUT 1.0e-9

/*************************************************************************/

double fd3sep ( long K, long M, long N, double **dftobs, double *sig,
   double **rvm, double **lfm, double **dftmod, double **dftres ) {

   long i, j, k, n;
   double s2;
   gsl_matrix *A, *U, *X, *V;
   gsl_vector *S, *w, *b, *x;

   /* prepare */

   A = gsl_matrix_alloc ( 2*M, 2*K );
   U = gsl_matrix_alloc ( 2*M, 2*K );
   X = gsl_matrix_alloc ( 2*K, 2*K );
   V = gsl_matrix_alloc ( 2*K, 2*K );
   S = gsl_vector_alloc ( 2*K );
   w = gsl_vector_alloc ( 2*K );
   b = gsl_vector_alloc ( 2*M );
   x = gsl_vector_alloc ( 2*K );

   /* for each DFT component */

   for ( s2 = n = 0 ; n <= N/2 ; n++ ) {

      /* assemble data vector and model matrix */

      for ( j = 0 ; j < M ; j++ ) {
         double s = *(sig+j);
         gsl_vector_set ( b, 2*j,   *(*(dftobs+j)+2*n)  /s );
         gsl_vector_set ( b, 2*j+1, *(*(dftobs+j)+2*n+1)/s );
         for ( k = 0 ; k < K ; k++ ) {
            double q, v, fv, rez, imz;
            q = 2.0 * M_PI * ((double)n)/((double)N);
            v = *(*(rvm+k)+j);
            fv = floor(v);
            rez = *(*(lfm+k)+j) *
                ( ((fv+1)-v)*cos(fv*q) + (v-fv)*cos((fv+1)*q) );
            imz = - *(*(lfm+k)+j) *
                ( ((fv+1)-v)*sin(fv*q) + (v-fv)*sin((fv+1)*q) );
            gsl_matrix_set ( A, 2*j,   2*k,     rez/s );
            gsl_matrix_set ( A, 2*j,   2*k+1, - imz/s );
            gsl_matrix_set ( A, 2*j+1, 2*k,     imz/s );
            gsl_matrix_set ( A, 2*j+1, 2*k+1,   rez/s );
         }
         gsl_matrix_memcpy ( U, A );
      }

      /* solve for model parameters */

      /* gsl_linalg_SV_decomp ( U, V, S, w ); */
      /* gsl_linalg_SV_decomp_mod ( U, X, V, S, w ); */
      /* gsl_linalg_SV_decomp_jacobi ( U, V, S ); */

      gsl_linalg_SV_decomp_mod ( U, X, V, S, w );
      for ( k = 0 ; k < 2*K-1 ; k++ )
         if ( gsl_vector_get(S,2*K-1-k)/gsl_vector_get(S,0) < SVCUT )
            gsl_vector_set ( S, 2*K-1-k, 0 );
      gsl_linalg_SV_solve ( U, V, S, b, x );

      /* copy to output */

      for ( k = 0 ; k < K ; k++ ) {
         *(*(dftmod+k)+2*n  ) = gsl_vector_get ( x, 2*k   );
         *(*(dftmod+k)+2*n+1) = gsl_vector_get ( x, 2*k+1 );
      }

      /* compute s2 */

      for ( j = 0 ; j < M ; j++ ) for ( i = 0 ; i < 2 ; i++ ) {

        double bc, db, s = *(sig+j);

        for ( bc = k = 0 ; k < K ; k++ ) {
          bc += gsl_matrix_get( A, 2*j+i, 2*k   ) * gsl_vector_get( x, 2*k   );
          bc += gsl_matrix_get( A, 2*j+i, 2*k+1 ) * gsl_vector_get( x, 2*k+1 );
        }
        db = gsl_vector_get ( b, 2*j+i ) - bc;
        s2 += db * db * ( n % ((N+1)/2) ? 2 : 1 );
        *(*(dftres+j)+2*n+i) = db * s;

      }

   }

   /* close */

   gsl_matrix_free ( A );
   gsl_matrix_free ( U );
   gsl_matrix_free ( X );
   gsl_matrix_free ( V );
   gsl_vector_free ( S );
   gsl_vector_free ( w );
   gsl_vector_free ( b );
   gsl_vector_free ( x );

   return s2;

}

/*************************************************************************/

void dft_fwd ( long m, long n, double **mxin, double **mxout ) {

   long i, j;
   double a = 1.0 / sqrt(n);
   gsl_fft_real_wavetable * dft_rewt;
   gsl_fft_real_workspace * dft_rews;

   dft_rewt = gsl_fft_real_wavetable_alloc (n);
   dft_rews = gsl_fft_real_workspace_alloc (n);

   for ( j = 0 ; j < m ; j++ ) {
      for ( i = 0; i < n; i++ ) *(*(mxout+j)+i+1) = a * *(*(mxin+j)+i);
      gsl_fft_real_transform ( *(mxout+j)+1, 1, n, dft_rewt, dft_rews );
      *(*(mxout+j)+0) = *(*(mxout+j)+1);
      *(*(mxout+j)+1) = 0;
      if ( ! n % 2 ) *(*(mxout+j)+n+1) = 0; /* if n even */
   }

   gsl_fft_real_wavetable_free ( dft_rewt );
   gsl_fft_real_workspace_free ( dft_rews );

}

/*************************************************************************/

void dft_bck ( long m, long n, double **mxin, double **mxout ) {

   long i, j;
   double a = 1.0 / sqrt(n);
   gsl_fft_halfcomplex_wavetable * dft_hcwt;
   gsl_fft_real_workspace * dft_rews;

   dft_hcwt = gsl_fft_halfcomplex_wavetable_alloc (n);
   dft_rews = gsl_fft_real_workspace_alloc (n);

   for ( j = 0 ; j < m ; j++ ) {
      *(*(mxout+j)+0) = *(*(mxin+j)+0) / a;
      for ( i = 1; i < n; i++ ) *(*(mxout+j)+i) = *(*(mxin+j)+i+1) / a;
      gsl_fft_halfcomplex_inverse ( *(mxout+j), 1, n, dft_hcwt, dft_rews );
   }

   gsl_fft_halfcomplex_wavetable_free ( dft_hcwt );
   gsl_fft_real_workspace_free ( dft_rews );

}

/*************************************************************************/

max_stars_repo_name: mrawls/FDBinary-tools

max_issues_repo_name: mrawls/FDBinary-tools

max_forks_repo_name: mrawls/FDBinary-tools


text:
/**
 * @file sad.h
 * @brief Voice activity detection.
 * @author Kenichi Kumatani and John McDonough
 */

#ifndef SAD_H
#define SAD_H

#include <stdio.h>
#include <assert.h>
#define _LOG_SAD_
#ifdef _LOG_SAD_
#include <stdarg.h>
#endif /* _LOG_SAD_ */

#include <gsl/gsl_block.h>
#include <gsl/gsl_vector.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_complex.h>
#include <gsl/gsl_complex_math.h>
#include <gsl/gsl_fit.h>

#include "common/jexception.h"
#include "stream/stream.h"
#include "feature/feature.h"
#include "sad/neural_spnsp_incl.h"

// ----- definition for abstract base class `NeuralNetVAD' -----
//
class NeuralNetVAD {
public:
  NeuralNetVAD(VectorFloatFeatureStreamPtr& cep,
				  unsigned context = 4, unsigned hiddenUnitsN = 1000, unsigned outputUnitsN = 2, float threshold = 0.1, const String& neuralNetFile = "");
  ~NeuralNetVAD();

  bool next(int frame_no = -5);
  void reset();
  void read(const String& neuralNetFile);

private:
  void shift_down_();
  void increment_() { frame_no_++; }
  void update_buffer_(int frame_no);

  const VectorFloatFeatureStreamPtr			cep_;
  const int						frame_reset_no_;
  const unsigned					cepLen_;
  bool							is_speech_;
  int							frame_no_;
  unsigned						framesPadded_;
  unsigned						context_;
  unsigned						hiddenUnitsN_;
  unsigned						outputUnitsN_;
  float							threshold_;
  MLP*							mlp_;
  float**						frame_;
};

typedef refcount_ptr<NeuralNetVAD> NeuralNetVADPtr;

// ----- definition for abstract base class `VAD' -----
//
class VAD {
 public:
  VAD(VectorComplexFeatureStreamPtr& samp);
  ~VAD();

  virtual bool next(int frame_no = -5) = 0;
  virtual void reset() { frame_no_ = frame_reset_no_; }
  const gsl_vector_complex* frame() const { return frame_; }
  virtual void next_speaker() = 0;

 protected:
  void increment_() { frame_no_++; }

  const VectorComplexFeatureStreamPtr			samp_;
  const int						frame_reset_no_;
  const unsigned					fftLen_;
  bool							is_speech_;
  int							frame_no_;
  gsl_vector_complex*					frame_;
};

typedef refcount_ptr<VAD> VADPtr;

// ----- definition for class `EnergyVAD' -----
//
class SimpleEnergyVAD : public VAD {
 public:
  SimpleEnergyVAD(VectorComplexFeatureStreamPtr& samp,
			       double threshold, double gamma = 0.995);
  ~SimpleEnergyVAD();

  virtual bool next(int frame_no = -5);
  virtual void reset();
  virtual void next_speaker();

#ifdef ENABLE_LEGACY_BTK_API__
  void nextSpeaker(){ next_speaker(); }
#endif

 private:
  const double						threshold_;
  const double						gamma_;
  double						spectral_energy_;
};

typedef Inherit<SimpleEnergyVAD, VADPtr> SimpleEnergyVADPtr;

// ----- definition for class `SimpleLikelihoodRatioVAD' -----
//
class SimpleLikelihoodRatioVAD : public VAD {
 public:
  SimpleLikelihoodRatioVAD(VectorComplexFeatureStreamPtr& samp,
					double threshold, double alpha);
  ~SimpleLikelihoodRatioVAD();

  virtual bool next(int frame_no = -5);
  virtual void reset();
  virtual void next_speaker();
  void set_variance(const gsl_vector* variance);

#ifdef ENABLE_LEGACY_BTK_API__
  void setVariance(const gsl_vector* variance);
  void nextSpeaker(){ next_speaker(); }
#endif

 private:
  double calc_Ak_(double vk, double gammak, double Rk);

  bool							variance_set_;
  gsl_vector*						noise_variance_;
  gsl_vector*						prev_Ak_;
  gsl_vector_complex*					prev_frame_;
  double						threshold_;
  double						alpha_;
};

typedef Inherit<SimpleLikelihoodRatioVAD, VADPtr> SimpleLikelihoodRatioVADPtr;

// ----- definition for class `EnergyVADFeature' -----
//
class EnergyVADFeature : public VectorFloatFeatureStream {
 public:
  EnergyVADFeature(const VectorFloatFeatureStreamPtr& source, double threshold = 0.5, unsigned bufferLength = 30, unsigned energiesN = 200, const String& nm = "Energy VAD");
  virtual ~EnergyVADFeature();

  virtual const gsl_vector_float* next(int frame_no = -5);
  virtual void reset();
  void next_speaker();

#ifdef ENABLE_LEGACY_BTK_API__
  void nextSpeaker(){ next_speaker(); }
#endif

 private:
  static int comparator_(const void* elem1, const void* elem2);
  virtual bool above_threshold_(const gsl_vector_float* vector);

  VectorFloatFeatureStreamPtr			source_;
  bool						recognizing_;

  gsl_vector_float**				buffer_;
  const unsigned				bufferLen_;
  int						bufferX_;
  unsigned					bufferedN_;

  unsigned					abovethresholdN_;
  unsigned					belowThresholdN_;

  unsigned					energiesN_;
  double*					energies_;
  double*					sorted_energies_;
  const unsigned				medianX_;
};

typedef Inherit<EnergyVADFeature, VectorFloatFeatureStreamPtr> EnergyVADFeaturePtr;

// ----- definition for abstract base class `VADMetric' -----
//
class VADMetric :  public Countable {
public:
  VADMetric();
  ~VADMetric();

  virtual double next(int frame_no = -5) = 0;
  virtual void reset() = 0;
  double score(){ return cur_score_;}
  virtual void next_speaker() = 0;

#ifdef  _LOG_SAD_
  bool openLogFile( const String & logfilename );
  int  writeLog( const char *format, ... );
  void closeLogFile();
  void initScore(){ _score=0.0; _scoreX=0; }
  void setScore(double score){ _score=score; _scoreX++; }
  double getAverageScore(){ if(_scoreX==0){return 0;} return(_score/_scoreX); }

  int   frame_no_;
#endif /* _LOG_SAD_ */

protected:
  double                                        cur_score_;

#ifdef  _LOG_SAD_
private:
  FILE *logfp_;
  double _score;
  unsigned _scoreX;
#endif /* _LOG_SAD_ */
};

typedef refcountable_ptr<VADMetric> VADMetricPtr;

// ----- definition for class `EnergyVADMetric' -----
//
class EnergyVADMetric : public VADMetric {
public:
  EnergyVADMetric(const VectorFloatFeatureStreamPtr& source, double initialEnergy = 5.0e+07, double threshold = 0.5, unsigned headN = 4,
		  unsigned tailN = 10, unsigned energiesN = 200, const String& nm = "Energy VAD Metric");

  ~EnergyVADMetric();

  virtual double next(int frame_no = -5);
  virtual void reset();
  virtual void next_speaker();
  double energy_percentile(double percentile = 50.0) const;

#ifdef ENABLE_LEGACY_BTK_API__
  void nextSpeaker(){ next_speaker(); }
  double energyPercentile(double percentile = 50.0){ return energy_percentile(percentile); }
#endif

private:
  static int comparator_(const void* elem1, const void* elem2);
  virtual bool above_threshold_(const gsl_vector_float* vector);

  VectorFloatFeatureStreamPtr			source_;
  const double					initial_energy_;

  const unsigned				headN_;
  const unsigned				tailN_;
  bool						recognizing_;

  unsigned					abovethresholdN_;
  unsigned					belowThresholdN_;

  unsigned					energiesN_;
  double*					energies_;
  double*					sorted_energies_;
  const unsigned				medianX_;
};

typedef Inherit<EnergyVADMetric, VADMetricPtr> EnergyVADMetricPtr;

// ----- definition for class `MultiChannelVADMetric' -----
//
template <typename ChannelType>
class MultiChannelVADMetric : public VADMetric {
 public:
  MultiChannelVADMetric(unsigned fftLen,
			double sampleRate, double lowCutoff, double highCutoff, const String& nm);
  ~MultiChannelVADMetric();

  void set_channel(ChannelType& chan);

#ifdef ENABLE_LEGACY_BTK_API__
  void setChannel(ChannelType& chan){ set_channel(chan); }
#endif

protected:
  unsigned set_lowX_(double lowCutoff) const;
  unsigned set_highX_(double highCutoff) const;
  unsigned set_binN_() const;

  typedef list<ChannelType>      	ChannelList_;
  ChannelList_				_channelList;

  const unsigned			fftLen_;
  const unsigned			fftLen2_;
  const double				samplerate_;
  const unsigned			lowX_;
  const unsigned			highX_;
  const unsigned			binN_;
  FILE                                 *logfp_;
};

template<> MultiChannelVADMetric<VectorFloatFeatureStreamPtr>::MultiChannelVADMetric(unsigned fftLen, double sampleRate, double lowCutoff, double highCutoff, const String& nm);
template<> MultiChannelVADMetric<VectorComplexFeatureStreamPtr>::MultiChannelVADMetric(unsigned fftLen, double sampleRate, double lowCutoff, double highCutoff, const String& nm);
template<> MultiChannelVADMetric<VectorFloatFeatureStreamPtr>::~MultiChannelVADMetric();
template<> MultiChannelVADMetric<VectorComplexFeatureStreamPtr>::~MultiChannelVADMetric();
template<> void MultiChannelVADMetric<VectorFloatFeatureStreamPtr>::set_channel(VectorFloatFeatureStreamPtr& chan);
template<> void MultiChannelVADMetric<VectorComplexFeatureStreamPtr>::set_channel(VectorComplexFeatureStreamPtr& chan);

typedef MultiChannelVADMetric<VectorFloatFeatureStreamPtr>   FloatMultiChannelVADMetric;
typedef MultiChannelVADMetric<VectorComplexFeatureStreamPtr> ComplexMultiChannelVADMetric;
typedef refcountable_ptr<FloatMultiChannelVADMetric>         FloatMultiChannelVADMetricPtr;
typedef refcountable_ptr<ComplexMultiChannelVADMetric>       ComplexMultiChannelVADMetricPtr;

// ----- definition for class `PowerSpectrumVADMetric' -----
//
/**
   @brief detect voice activity based on the energy comparison
   @usage
   1. construct the object with PowerSpectrumVADMetric().
   2. set the channel data with set_channel().
   3. call next().
   @note the first channel is associated with the target speaker.
 */
class PowerSpectrumVADMetric : public FloatMultiChannelVADMetric {
public:
  PowerSpectrumVADMetric(unsigned fftLen,
			 double sampleRate = 16000.0, double lowCutoff = -1.0, double highCutoff = -1.0,
			 const String& nm = "Power Spectrum VAD Metric");
  PowerSpectrumVADMetric(VectorFloatFeatureStreamPtr& source1, VectorFloatFeatureStreamPtr& source2,
			 double sampleRate = 16000.0, double lowCutoff = -1.0, double highCutoff = -1.0,
			 const String& nm = "Power Spectrum VAD Metric");

  ~PowerSpectrumVADMetric();

  virtual double next(int frame_no = -5);
  virtual void reset();
  virtual void next_speaker();
  gsl_vector *get_metrics() const { return powerList_;}
  void set_E0( double E0 ){ E0_ = E0;}
  void clear_channel();

#ifdef ENABLE_LEGACY_BTK_API__
  void nextSpeaker(){ next_speaker(); }
  gsl_vector *getMetrics(){ return get_metrics(); }
  void setE0( double E0 ){ set_E0(E0); }
  void clearChannel(){ clear_channel(); }
#endif

protected:
  typedef ChannelList_::iterator		ChannelIterator_;
  gsl_vector                                   *powerList_;
  double                                        E0_;
};

typedef Inherit<PowerSpectrumVADMetric, FloatMultiChannelVADMetricPtr> PowerSpectrumVADMetricPtr;

// ----- definition for class `NormalizedEnergyMetric' -----
//

/**
   @class

   @usage
   @note
*/

class NormalizedEnergyMetric : public  PowerSpectrumVADMetric {
  static double initial_energy_;
public:
  NormalizedEnergyMetric(unsigned fftLen,
		double sampleRate = 16000.0, double lowCutoff = -1.0, double highCutoff = -1.0,
		const String& nm = "TSPS VAD Metric");

  ~NormalizedEnergyMetric();

  virtual double next(int frame_no = -5);
  virtual void reset();
};

typedef Inherit<NormalizedEnergyMetric, PowerSpectrumVADMetricPtr> NormalizedEnergyMetricPtr;

// ----- definition for class `CCCVADMetric' -----
//

/**
   @class compute cross-correlation coefficients (CCC) as a function of time delays,
          and average up the n-best values.
   @usage
   @note
*/

class CCCVADMetric : public ComplexMultiChannelVADMetric {
public:
  CCCVADMetric(unsigned fftLen, unsigned nCand,
	       double sampleRate = 16000.0, double lowCutoff = -1.0, double highCutoff = -1.0,
	       const String& nm = "CCC VAD Metric");
  ~CCCVADMetric();

  virtual double next(int frame_no = -5);
  virtual void reset();
  virtual void next_speaker();
  void set_NCand(unsigned nCand);
  void set_threshold(double threshold){ threshold_ = threshold;}
  gsl_vector *get_metrics() const { return ccList_;}
  void clear_channel();

#ifdef ENABLE_LEGACY_BTK_API__
  void nextSpeaker(){ next_speaker(); }
  void setNCand(unsigned nCand){ set_NCand(nCand); }
  void setThreshold(double threshold){ set_threshold(threshold); }
  gsl_vector *getMetrics(){ return get_metrics(); }
  void clearChannel(){ clear_channel(); }
#endif

protected:
  typedef ChannelList_::iterator		ChannelIterator_;
  unsigned                                      nCand_;
  gsl_vector                                   *ccList_;
  gsl_vector_int                               *sample_delays_;
  double                                        threshold_;
  double                                       *pack_cross_spectrum_;
};

typedef Inherit<CCCVADMetric, ComplexMultiChannelVADMetricPtr> CCCVADMetricPtr;

// ----- definition for class `TSPSVADMetric' -----
//

/**
   @class

   @usage
   @note
*/

class TSPSVADMetric : public  PowerSpectrumVADMetric {
  static double initial_energy_;
public:
  TSPSVADMetric(unsigned fftLen,
		double sampleRate = 16000.0, double lowCutoff = -1.0, double highCutoff = -1.0,
		const String& nm = "TSPS VAD Metric");
  ~TSPSVADMetric();

  virtual double next(int frame_no = -5);
  virtual void reset();
};

typedef Inherit<TSPSVADMetric,PowerSpectrumVADMetricPtr> TSPSVADMetricPtr;

// ----- definition for class `NegentropyVADMetric' -----
//
class NegentropyVADMetric : public VADMetric {

protected:
  class ComplexGeneralizedGaussian_ : public Countable {
  public:
    ComplexGeneralizedGaussian_(double shapeFactor = 2.0);
    double logLhood(gsl_complex X, double scaleFactor) const;
    double shapeFactor() const { return shape_factor_; }
    double Bc() const { return Bc_; }
    double normalization() const { return normalization_; }

  protected:
    virtual double calc_Bc_() const;
    virtual double calc_normalization_() const;

    const double				shape_factor_;
    /* const */ double				Bc_;
    /* const */ double				normalization_;
  };

  typedef refcountable_ptr<ComplexGeneralizedGaussian_> ComplexGeneralizedGaussianPtr_;

  typedef list<ComplexGeneralizedGaussianPtr_>	GaussianList_;
  typedef GaussianList_::iterator		GaussianListIterator_;
  typedef GaussianList_::const_iterator		GaussianListConstIterator_;

public:
  NegentropyVADMetric(const VectorComplexFeatureStreamPtr& source, const VectorFloatFeatureStreamPtr& spectralEstimator,
		      const String& shapeFactorFileName = "", double threshold = 0.5,
		      double sampleRate = 16000.0, double lowCutoff = -1.0, double highCutoff = -1.0,
		      const String& nm = "Negentropy VAD Metric");

  ~NegentropyVADMetric();

  virtual double next(int frame_no = -5);
  virtual void reset();
  virtual void next_speaker();
  double calc_negentropy(int frame_no);

#ifdef ENABLE_LEGACY_BTK_API__
  void nextSpeaker(){ next_speaker(); }
  double calcNegentropy(int frame_no){ return calc_negentropy(frame_no); }
#endif

protected:
  virtual bool above_threshold_(int frame_no);
  unsigned set_lowX_(double lowCutoff) const;
  unsigned set_highX_(double highCutoff) const;
  unsigned set_binN_() const;

  VectorComplexFeatureStreamPtr			source_;
  VectorFloatFeatureStreamPtr			spectral_estimator_;

  GaussianList_					generalized_gaussians_;
  ComplexGeneralizedGaussianPtr_		gaussian_;

  const double					threshold_;

  const unsigned				fftLen_;
  const unsigned				fftLen2_;

  const double					samplerate_;
  const unsigned				lowX_;
  const unsigned				highX_;
  const unsigned				binN_;
};

typedef Inherit<NegentropyVADMetric, VADMetricPtr> NegentropyVADMetricPtr;

// ----- definition for class `MutualInformationVADMetric' -----
//
class MutualInformationVADMetric : public NegentropyVADMetric {
 public:
  MutualInformationVADMetric(const VectorComplexFeatureStreamPtr& source1,
                             const VectorComplexFeatureStreamPtr& source2,
                             const VectorFloatFeatureStreamPtr& spectralEstimator1,
                             const VectorFloatFeatureStreamPtr& spectralEstimator2,
                             const String& shapeFactorFileName = "",
                             double twiddle = -1.0,
                             double threshold = 1.3,
                             double beta = 0.95,
                             double sampleRate = 16000.0,
                             double lowCutoff = 187.0,
                             double highCutoff = 1000.0,
                             const String& nm = "Mutual Information VAD Metric");
  ~MutualInformationVADMetric();

  virtual double next(int frame_no = -5);
  virtual void reset();
  virtual void next_speaker();
  double calc_mutual_information(int frame_no);

#ifdef ENABLE_LEGACY_BTK_API__
  double calcMutualInformation(int frame_no){ return calc_mutual_information(frame_no); }
#endif

 protected:
  class JointComplexGeneralizedGaussian_ : public NegentropyVADMetric::ComplexGeneralizedGaussian_ {
 public:
   JointComplexGeneralizedGaussian_(const NegentropyVADMetric::ComplexGeneralizedGaussianPtr_& ggaussian);
   ~JointComplexGeneralizedGaussian_();

   double logLhood(gsl_complex X1, gsl_complex X2, double scaleFactor1, double scaleFactor2, gsl_complex rho12) const;

 private:
   static const double      sqrt_two_;
   static const gsl_complex complex_one_;
   static const gsl_complex complex_zero_;

   virtual double calc_Bc_() const;
   virtual double calc_normalization_() const;

   double lngamma_ratio_(double f) const;
   double lngamma_ratio_joint_(double f) const;
   double match_(double f) const;

   double match_score_marginal_(double f) const;
   double match_score_joint_(double fJ) const;

   static const double				tolerance_;

   gsl_vector_complex*				X_;
   gsl_vector_complex*				scratch_;
   gsl_matrix_complex*				SigmaX_inverse_;
 };

 protected:
 typedef refcountable_ptr<JointComplexGeneralizedGaussian_> JointComplexGeneralizedGaussianPtr_;

 typedef list<JointComplexGeneralizedGaussianPtr_>	JointGaussianList_;
 typedef JointGaussianList_::iterator			JointGaussianListIterator_;
 typedef JointGaussianList_::const_iterator		JointGaussianListConstIterator_;

 typedef vector<gsl_complex>				CrossCorrelationVector_;

 static const double				epsilon_;

 virtual bool above_threshold_(int frame_no);
 double calc_fixed_threshold_();
 double calc_total_threshold_() const;
 void initialize_pdfs_();

 VectorComplexFeatureStreamPtr			source2_;
 VectorFloatFeatureStreamPtr			spectral_estimator2_;

 JointGaussianList_				joint_generalized_gaussians_;
 CrossCorrelationVector_			ccs_; // cross correlations
 const double					fixed_threshold_;
 const double					twiddle_;
 const double					threshold_;
 const double					beta_;
};

typedef Inherit<MutualInformationVADMetric, NegentropyVADMetricPtr> MutualInformationVADMetricPtr;

// ----- definition for class `LikelihoodRatioVADMetric' -----
//
class LikelihoodRatioVADMetric : public NegentropyVADMetric {
 public:
  LikelihoodRatioVADMetric(const VectorComplexFeatureStreamPtr& source1, const VectorComplexFeatureStreamPtr& source2,
			   const VectorFloatFeatureStreamPtr& spectralEstimator1, const VectorFloatFeatureStreamPtr& spectralEstimator2,
			   const String& shapeFactorFileName = "", double threshold = 0.0,
			   double sampleRate = 16000.0, double lowCutoff = 187.0, double highCutoff = 1000.0,
			   const String& nm = "Likelihood VAD Metric");
  ~LikelihoodRatioVADMetric();

  virtual double next(int frame_no = -5);
  virtual void reset();
  virtual void next_speaker();
  double calc_likelihood_ratio(int frame_no);

#ifdef ENABLE_LEGACY_BTK_API__
  void nextSpeaker(){ next_speaker(); }
  double calcLikelihoodRatio(int frame_no){ return calc_likelihood_ratio(frame_no); }
#endif

 private:
  const VectorComplexFeatureStreamPtr			source2_;
  const VectorFloatFeatureStreamPtr			spectral_estimator2_;
};

typedef Inherit<LikelihoodRatioVADMetric, NegentropyVADMetricPtr> LikelihoodRatioVADMetricPtr;

// ----- definition for class `LowFullBandEnergyRatioVADMetric' -----
//
class LowFullBandEnergyRatioVADMetric : public VADMetric {
public:
  LowFullBandEnergyRatioVADMetric(const VectorFloatFeatureStreamPtr& source, const gsl_vector* lowpass, double threshold = 0.5, const String& nm = "Low- Full-Band Energy Ratio VAD Metric");
  ~LowFullBandEnergyRatioVADMetric();

  virtual double next(int frame_no = -5);
  virtual void reset();
  virtual void next_speaker();

#ifdef ENABLE_LEGACY_BTK_API__
  void nextSpeaker(){ next_speaker(); }
#endif

private:
  virtual bool above_threshold_(int frame_no);
  void calc_auto_correlation_vector_(int frame_no);
  void calc_covariance_matrix_();
  double calc_lower_band_energy_();

  VectorFloatFeatureStreamPtr			source_;
  const unsigned				_lagsN;
  gsl_vector*					_lowpass;
  gsl_vector*					scratch_;
  double*					_autocorrelation;
  gsl_matrix*					_covariance;
};

typedef Inherit<LowFullBandEnergyRatioVADMetric, VADMetricPtr> LowFullBandEnergyRatioVADMetricPtr;

// ----- definition for class `HangoverVADFeature' -----
//
class HangoverVADFeature : public VectorFloatFeatureStream {
 public:
  HangoverVADFeature(const VectorFloatFeatureStreamPtr& source, const VADMetricPtr& metric, double threshold = 0.5,
		     unsigned headN = 4, unsigned tailN = 10, const String& nm = "Hangover VAD Feature");
  virtual ~HangoverVADFeature();

  virtual const gsl_vector_float* next(int frame_no = -5);
  virtual void reset();
  virtual void next_speaker();
  int prefixN() const { return prefixN_ - headN_; }

#ifdef ENABLE_LEGACY_BTK_API__
  void nextSpeaker(){ next_speaker(); }
#endif

 protected:
  typedef pair<VADMetricPtr, double>		MetricPair_;
  typedef vector<MetricPair_>			MetricList_;
  typedef MetricList_::iterator			MetricListIterator_;
  typedef MetricList_::const_iterator		MetricListConstIterator_;

  static const unsigned EnergyVADMetricX		= 0;
  static const unsigned MutualInformationVADMetricX	= 1;
  static const unsigned LikelihoodRatioVADMetricX	= 2;

  static int comparator_(const void* elem1, const void* elem2);
  virtual bool above_threshold_(int frame_no);

  VectorFloatFeatureStreamPtr			source_;
  bool						recognizing_;

  gsl_vector_float**				buffer_;
  const unsigned				headN_;
  const unsigned				tailN_;
  int						bufferX_;
  unsigned					bufferedN_;

  unsigned					abovethresholdN_;
  unsigned					belowThresholdN_;
  unsigned					prefixN_;

  MetricList_					metricList_;
};

typedef Inherit<HangoverVADFeature, VectorFloatFeatureStreamPtr> HangoverVADFeaturePtr;

// ----- definition for class `HangoverMIVADFeature' -----
//
class HangoverMIVADFeature : public HangoverVADFeature {
  static const unsigned EnergyVADMetricX		= 0;
  static const unsigned MutualInformationVADMetricX	= 1;
  static const unsigned PowerVADMetricX			= 2;

 public:
  HangoverMIVADFeature(const VectorFloatFeatureStreamPtr& source,
		       const VADMetricPtr& energyMetric, const VADMetricPtr& mutualInformationMetric, const VADMetricPtr& powerMetric,
		       double energythreshold = 0.5, double mutualInformationThreshold = 0.5, double powerThreshold = 0.5,
		       unsigned headN = 4, unsigned tailN = 10, const String& nm = "Hangover VAD Feature");

  int decision_metric() const { return decision_metric_; }

#ifdef ENABLE_LEGACY_BTK_API__
  int decisionMetric() { return decision_metric(); }
#endif

protected:
  virtual bool above_threshold_(int frame_no);

  int						decision_metric_;
};

typedef Inherit<HangoverMIVADFeature, HangoverVADFeaturePtr> HangoverMIVADFeaturePtr;

// ----- definition for class `HangoverMultiStageVADFeature' -----
//
class HangoverMultiStageVADFeature : public HangoverVADFeature {
 public:
  HangoverMultiStageVADFeature(const VectorFloatFeatureStreamPtr& source,
			       const VADMetricPtr& energyMetric, double energythreshold = 0.5,
			       unsigned headN = 4, unsigned tailN = 10, const String& nm = "HangoverMultiStageVADFeature");
  ~HangoverMultiStageVADFeature();
  int decision_metric() const { return decision_metric_; }
  void set_metric( const VADMetricPtr& metricPtr, double threshold ){
    metricList_.push_back( MetricPair_( metricPtr, threshold ) );
  }

#ifdef ENABLE_LEGACY_BTK_API__
  int decisionMetric() { return decision_metric(); }
  void setMetric( const VADMetricPtr& metricPtr, double threshold ){ set_metric(metricPtr, threshold); }
#endif

#ifdef _LOG_SAD_
  void initScores();
  gsl_vector *getScores();
#endif /* _LOG_SAD_ */

protected:
  virtual bool above_threshold_(int frame_no);
  int          decision_metric_;
#ifdef _LOG_SAD_
  gsl_vector *_scores;
#endif /* _LOG_SAD_ */
};

typedef Inherit<HangoverMultiStageVADFeature, HangoverVADFeaturePtr> HangoverMultiStageVADFeaturePtr;

#endif