#ifndef   math_utility_H__
#define   math_utility_H__

#include <cstdlib>
#include <vector>
#include <algorithm>
#include <cmath>

namespace meow{

//! 圓周率...
static const double PI = 3.14159265358979323846264338327950288;

/*!
 * @brief 如果abs(輸入的數值) < eps, 則回傳0, 否則回傳輸入的數值
 */
template<class T>
inline T noEPS(T value, T eps = 1e-9){
  T epsp((eps < T(0)) ? -eps : eps);
  return ((value < -epsp || value > epsp) ? value : T(0));
}

/*!
 * @brief \c (value-lower)/(upper-lower)
 */
template<class T>
inline T normalize(T lower, T upper, T value){
  return (value - lower) / (upper - lower);
}

/*!
 * @brief \c (lower+_ratio*(upper-lower))
 */
template<class T>
inline T denormalize(T lower, T upper, T _ratio){
  return lower + _ratio * (upper - lower);
}

/*!
 * @brief \c denormalize(l2,u2,normalize(l1,u1,m1))
 */
template<class T>
inline T ratioMapping(T l1, T u1, T m1, T l2, T u2){
  return denormalize(l2, u2, normalize(l1, u1, m1));
}

/*!
 * @brief \c std::min(mx,std::max(mn,v))
 */
template<class T>
inline T inRange(T const& mn, T const& mx, T const& v){
  return std::min(mx, std::max(mn, v));
}

/*!
 * @brief \c x*x
 */
template<class T>
inline T squ(T const& x){
  return x * x;
}

/*!
 * @brief \c x*x*x
 */
template<class T>
inline T cub(T const& x){
  return x * x * x;
}

/*!
 * @brief 只將 \c sigs 個標準差以內的數據拿來取平均
 */
template<class T>
inline double average(T const& beg, T const& end, double sigs){
  int N = 0;
  double av = 0;
  for(T it = beg; it != end; it++, N++){
    av += *it;
  }
  av /= N;
  double sig = 0;
  for(T it = beg; it != end; it++){
    sig += (*it - av) * (*it - av);
  }
  sig = sqrt(sig / N);
  double lower = av - sig * sigs, upper = av + sig * sigs;
  double ret = 0, retn = 0;
  for(T it = beg; it != end; it++){
    if(lower <= *it && *it <= upper){
      ret += *it;
      retn++;
    }
  }
  return ret / retn;
}

/*!
 * @brief 只將 \c sigs 個標準差以內的數據拿來取平均, 不過這次用 \c p 來加權平均
 */
template<class T>
inline double average(T const& beg, T const& end, T const& p, double sigs){
  int N = 0;
  double ps = 0;
  for(T it = beg, ip = p; it != end; it++, N++, ip++){
    ps += *ip;
  }
  double av = 0;
  for(T it = beg, ip = p; it != end; it++, ip++){
    av += *it * *ip / ps;
  }
  double sig = 0;
  for(T it = beg, ip = p; it != end; it++, ip++){
    sig += *ip / ps * (*it - av) * (*it - av);
  }
  sig = sqrt(sig);
  double lower = av - sig * sigs, upper = av + sig * sigs;
  double ret = 0, retn = 0;
  for(T it = beg, ip = p; it != end; it++, ip++){
    if(lower <= *it && *it <= upper){
      ret  += *it * *ip;
      retn += *ip;
    }
  }
  if(retn <= 1e-10) return av;
  return ret / retn;
}

/*!
 * @brief 就只是個取絕對值
 */
template<class T>
inline T tAbs(T const& t){
  return (t < 0 ? -t : t);
}

}

#endif // math_utility_H__