ElementFunctions.h

Go to the documentation of this file.

#ifndef ELEMENT_FUNCTIONS_H
#define ELEMENT_FUNCTIONS_H
 
#include <cmath>
#include <complex>
#include <limits>
 
namespace casacore {
namespace arrays_internal {
 
inline bool near(unsigned val1, unsigned val2, double tol) {
  if (tol <= 0) {
    return (val1 == val2);
  }
  if (val1 == val2) {
    return true;
  } else if (val1 > val2) {
    return (double(val1 - val2) <= tol * std::max(val1, val2));
  } else {
    return (double(val2 - val1) <= tol * std::max(val1, val2));
  }
}
 
inline bool near(int val1, int val2, double tol) {
  if (tol <= 0) {
    return (val1 == val2);
  }
  if (val1 == val2) {
    return true;
  }
  if ((0 < val1) != (0 < val2)) {
    return false;
  }
  const int aval1 = std::abs(val1);
  const int aval2 = std::abs(val2);
  return (double(aval1 - aval2) <= tol * double(std::max(aval1, aval2)));
}
 
inline bool near(float val1, float val2, double tol) {
  if (tol <= 0) {
    return (val1 == val2);
  }
  if (val1 == val2) {
    return true;
  }
  if (val1 == 0) {
    return (std::abs(val2) <= (1 + tol) * std::numeric_limits<float>::min());
  } else if (val2 == 0) {
    return (std::abs(val1) <= (1 + tol) * std::numeric_limits<float>::min());
  }
  if ((0 < val1) != (0 < val2)) {
    return false;
  }
  return (std::abs(val1 - val2) <= tol * std::max(std::abs(val1), std::abs(val2)));
}
 
inline bool near(double val1, double val2, double tol) {
  if (tol <= 0) {
    return (val1 == val2);
  }
  if (val1 == val2) {
    return true;
  }
  if (val1 == 0) {
    return (std::abs(val2) <= (1 + tol) * std::numeric_limits<double>::min());
  } else if (val2 == 0) {
    return (std::abs(val1) <= (1 + tol) * std::numeric_limits<double>::min());
  }
  if ((0 < val1) != (0 < val2)) {
    return false;
  }
  return (std::abs(val1 - val2) <= tol * std::max(std::abs(val1), std::abs(val2)));
}
 
inline bool near(float val1, double val2, double tol) { return near(double(val1), val2, tol); }
 
inline bool near(double val1, float val2, double tol) { return near(val1, double(val2), tol); }
 
inline bool near(const std::complex<float> &val1, const std::complex<float> &val2,
                 double tol = 1.0e-5) {
  if (tol <= 0) return val1 == val2;
  if (val1 == val2) return true;
  if (near(val1.real(), val2.real(), tol) && near(val1.imag(), val2.imag(), tol)) return true;
  float aval1(std::abs(val1)), aval2(std::abs(val2));
  if (aval1 == 0)
    return aval2 <= (1 + tol) * std::numeric_limits<float>::min();
  else if (aval2 == 0)
    return aval1 <= (1 + tol) * std::numeric_limits<float>::min();
  std::complex<double> dval(val1);
  dval -= std::complex<double>(val2);
  return std::abs(dval) <= tol * (aval1 < aval2 ? aval2 : aval1);
}
 
inline bool near(const std::complex<double> &val1, const std::complex<double> &val2,
                 double tol = 1.0e-13) {
  if (tol <= 0) return val1 == val2;
  if (val1 == val2) return true;
  if (std::abs(val1) == 0)
    return std::abs(val2) <= (1 + tol) * std::numeric_limits<double>::min();
  else if (std::abs(val2) == 0)
    return std::abs(val1) <= (1 + tol) * std::numeric_limits<double>::min();
  double aval1(std::abs(val1)), aval2(std::abs(val2));
  return std::abs(val1 - val2) <= tol * (aval1 < aval2 ? aval2 : aval1);
}
 
inline bool nearAbs(const std::complex<float> &val1, const std::complex<float> &val2,
                    double tol = 1.0e-5) {
  return std::abs(val2 - val1) <= tol;
}
 
inline bool nearAbs(const std::complex<double> &val1, const std::complex<double> &val2,
                    double tol = 1.0e-13) {
  return std::abs(val2 - val1) <= tol;
}
 
inline bool nearAbs(unsigned val1, unsigned val2, double tol) {
  if (val1 == val2) {
    return true;
  } else if (val1 > val2) {
    return (tol >= double(val1 - val2));
  } else {
    return (tol >= double(val2 - val1));
  }
}
 
inline bool nearAbs(int val1, int val2, double tol) {
  return (tol >= double(std::abs(val2 - val1)));
}
 
inline bool nearAbs(float val1, float val2, double tol) {
  return (tol >= double(std::abs(val2 - val1)));
}
 
inline bool nearAbs(double val1, double val2, double tol) { return (tol >= std::abs(val2 - val1)); }
 
// Define a function to compare all elements of two sequences.
// It returns true if all elements compare true.
// An example compare operator is <src>std::equal_to</src>.
// <group>
template <typename InputIterator1, typename InputIterator2, typename CompareOperator>
inline bool compareAll(InputIterator1 first1, InputIterator1 last1, InputIterator2 first2,
                       CompareOperator op) {
  for (; first1 != last1; ++first1, ++first2) {
    if (!op(*first1, *first2)) return false;
  }
  return true;
}
// For use with a constant left value.
// This avoids use of bind1st or bind2nd which can fail for gcc-4.3.
// (see ArrayMath.h).
template <typename InputIterator1, typename T, typename CompareOperator>
inline bool compareAllLeft(InputIterator1 first1, InputIterator1 last1, T left,
                           CompareOperator op) {
  for (; first1 != last1; ++first1) {
    if (!op(left, *first1)) return false;
  }
  return true;
}
// For use with a constant right value.
// This avoids use of bind1st or bind2nd which can fail for gcc-4.3.
// (see ArrayMath.h).
template <typename InputIterator1, typename T, typename CompareOperator>
inline bool compareAllRight(InputIterator1 first1, InputIterator1 last1, T right,
                            CompareOperator op) {
  for (; first1 != last1; ++first1) {
    if (!op(*first1, right)) return false;
  }
  return true;
}
// </group>
 
// Define a function to compare all elements of two sequences.
// It returns true if any element compares true.
// An example compare operator is <src>std::equal_to</src>.
// <group>
template <typename InputIterator1, typename InputIterator2, typename CompareOperator>
inline bool compareAny(InputIterator1 first1, InputIterator1 last1, InputIterator2 first2,
                       CompareOperator op) {
  for (; first1 != last1; ++first1, ++first2) {
    if (op(*first1, *first2)) return true;
  }
  return false;
}
// For use with a constant left value.
// This avoids use of bind1st or bind2nd which can fail for gcc-4.3.
// (see ArrayMath.h).
template <typename InputIterator1, typename T, typename CompareOperator>
inline bool compareAnyLeft(InputIterator1 first1, InputIterator1 last1, T left,
                           CompareOperator op) {
  for (; first1 != last1; ++first1) {
    if (op(left, *first1)) return true;
  }
  return false;
}
// For use with a constant right value.
// This avoids use of bind1st or bind2nd which can fail for gcc-4.3.
// (see ArrayMath.h).
template <typename InputIterator1, typename T, typename CompareOperator>
inline bool compareAnyRight(InputIterator1 first1, InputIterator1 last1, T right,
                            CompareOperator op) {
  for (; first1 != last1; ++first1) {
    if (op(*first1, right)) return true;
  }
  return false;
}
// </group>
 
// Functor to add squared diff of right and base value to left.
// It can be used to calculate the variance.
// Note: it is specialized for complex values to handle real and imag separately.
template <typename T, typename Accum = T>
struct SumSqrDiff {
  explicit SumSqrDiff(T base) : itsBase(base) {}
  Accum operator()(Accum left, T right) const {
    return left + (right - itsBase) * (right - itsBase);
  }
 
 private:
  Accum itsBase;  // store as Accum, so subtraction results in Accum
};
 
// Specialize for complex values.
// Variance has to be taken for the absolute value of a complex value. thus
//       sum(abs((a[i] - mean)**2
// where the sqrt used in abs and the **2 cancel each other, thus can be left out.
// See also https://en.wikipedia.org/wiki/Complex_random_variable#Variance
// Note that although the sum is real, a complex value is used to have equal template types.
template <typename T>
struct SumSqrDiff<std::complex<T>> {
  explicit SumSqrDiff(std::complex<T> base) : itsBase(base) {}
  std::complex<T> operator()(std::complex<T> left, std::complex<T> right) const {
    return left + ((right.real() - itsBase.real()) * (right.real() - itsBase.real()) +
                   (right.imag() - itsBase.imag()) * (right.imag() - itsBase.imag()));
  }
 
 private:
  std::complex<T> itsBase;
};
 
template <typename T>
bool isnan(const std::complex<T> &val) {
  return std::isnan(val.real()) || std::isnan(val.imag());
}
 
template <typename T>
bool isinf(const std::complex<T> &val) {
  return std::isinf(val.real()) || std::isinf(val.imag());
}
 
template <typename T>
bool isfinite(const std::complex<T> &val) {
  return std::isfinite(val.real()) || std::isfinite(val.imag());
}
 
inline int floormod(int x, int y) {
  int r = x % y;
  if (r != 0 && (x < 0) != (y < 0)) r += y;
  return r;
}
inline long long floormod(long long x, long long y) {
  long long r = x % y;
  if (r != 0 && (x < 0) != (y < 0)) r += y;
  return r;
}
inline float floormod(float x, float y) {
  float r = std::fmod(x, y);
  if (r != 0 && (x < 0) != (y < 0)) r += y;
  return r;
}
inline double floormod(double x, double y) {
  double r = std::fmod(x, y);
  if (r != 0 && (x < 0) != (y < 0)) r += y;
  return r;
}
 
template <class T, class F>
inline void convertScalar(T &out, F in) {
  out = static_cast<T>(in);
}
 
inline void convertScalar(std::complex<float> &out, std::complex<double> in) {
  out = std::complex<float>(in.real(), in.imag());
}
 
}  // namespace arrays_internal
}  // namespace casacore
 
#endif