text:
/*******************************************************************************
 *
 * A simple class for representing intervals and performing interval
 * arithmetic.
 *
 ******************************************************************************/

#pragma once

#include <boost/optional.hpp>
#include <crab/common/stats.hpp>
#include <crab/common/types.hpp>
#include <crab/domains/linear_interval_solver.hpp>
#include <crab/numbers/bignums.hpp>

namespace ikos {

template <typename Number> class bound {
public:
  typedef bound<Number> bound_t;

private:
  bool _is_infinite;
  Number _n;

private:
  bound();

  bound(bool is_infinite, Number n) : _is_infinite(is_infinite), _n(n) {
    if (is_infinite) {
      if (n > 0)
        _n = 1;
      else
        _n = -1;
    }
  }

public:
  static bound_t min(bound_t x, bound_t y) { return (x.operator<=(y) ? x : y); }

  static bound_t min(bound_t x, bound_t y, bound_t z) {
    return min(x, min(y, z));
  }

  static bound_t min(bound_t x, bound_t y, bound_t z, bound_t t) {
    return min(x, min(y, z, t));
  }

  static bound_t max(bound_t x, bound_t y) { return (x.operator<=(y) ? y : x); }

  static bound_t max(bound_t x, bound_t y, bound_t z) {
    return max(x, max(y, z));
  }

  static bound_t max(bound_t x, bound_t y, bound_t z, bound_t t) {
    return max(x, max(y, z, t));
  }

  static bound_t plus_infinity() { return bound_t(true, 1); }

  static bound_t minus_infinity() { return bound_t(true, -1); }

public:
  bound(int n) : _is_infinite(false), _n(n) {}

  bound(std::string s) : _n(1) {
    if (s == "+oo") {
      _is_infinite = true;
    } else if (s == "-oo") {
      _is_infinite = true;
      _n = -1;
    } else {
      _is_infinite = false;
      _n = Number(s);
    }
  }

  bound(Number n) : _is_infinite(false), _n(n) {}

  bound(const bound_t &o) : _is_infinite(o._is_infinite), _n(o._n) {}

  bound_t &operator=(const bound_t &o) {
    if (this != &o) {
      _is_infinite = o._is_infinite;
      _n = o._n;
    }
    return *this;
  }

  bool is_infinite() const { return _is_infinite; }

  bool is_finite() const { return !_is_infinite; }

  bool is_plus_infinity() const { return (is_infinite() && _n > 0); }

  bool is_minus_infinity() const { return (is_infinite() && _n < 0); }

  bound_t operator-() const { return bound_t(_is_infinite, -_n); }

  bound_t operator+(bound_t x) const {
    if (is_finite() && x.is_finite()) {
      return bound_t(_n + x._n);
    } else if (is_finite() && x.is_infinite()) {
      return x;
    } else if (is_infinite() && x.is_finite()) {
      return *this;
    } else if (_n == x._n) {
      return *this;
    } else {
      CRAB_ERROR("Bound: undefined operation -oo + +oo");
    }
  }

  bound_t &operator+=(bound_t x) { return operator=(operator+(x)); }

  bound_t operator-(bound_t x) const { return operator+(x.operator-()); }

  bound_t &operator-=(bound_t x) { return operator=(operator-(x)); }

  bound_t operator*(bound_t x) const {
    if (x._n == 0)
      return x;
    else if (_n == 0)
      return *this;
    else
      return bound_t(_is_infinite || x._is_infinite, _n * x._n);
  }

  bound_t &operator*=(bound_t x) { return operator=(operator*(x)); }

  bound_t operator/(bound_t x) const {
    if (x._n == 0) {
      CRAB_ERROR("Bound: division by zero");
    } else if (is_finite() && x.is_finite()) {
      return bound_t(false, _n / x._n);
    } else if (is_finite() && x.is_infinite()) {
      if (_n > 0) {
        return x;
      } else if (_n == 0) {
        return *this;
      } else {
        return x.operator-();
      }
    } else if (is_infinite() && x.is_finite()) {
      if (x._n > 0) {
        return *this;
      } else {
        return operator-();
      }
    } else {
      return bound_t(true, _n * x._n);
    }
  }

  bound_t &operator/=(bound_t x) { return operator=(operator/(x)); }

  bool operator<(bound_t x) const { return !operator>=(x); }

  bool operator>(bound_t x) const { return !operator<=(x); }

  bool operator==(bound_t x) const {
    return (_is_infinite == x._is_infinite && _n == x._n);
  }

  bool operator!=(bound_t x) const { return !operator==(x); }

  /*	operator<= and operator>= use a somewhat optimized implementation.
   *	results include up to 20% improvements in performance in the octagon
   *domain over a more naive implementation.
   */
  bool operator<=(bound_t x) const {
    if (_is_infinite xor x._is_infinite) {
      if (_is_infinite) {
        return _n < 0;
      }
      return x._n > 0;
    }
    return _n <= x._n;
  }

  bool operator>=(bound_t x) const {
    if (_is_infinite xor x._is_infinite) {
      if (_is_infinite) {
        return _n > 0;
      }
      return x._n < 0;
    }
    return _n >= x._n;
  }

  bound_t abs() const {
    if (operator>=(0)) {
      return *this;
    } else {
      return operator-();
    }
  }

  boost::optional<Number> number() const {
    if (is_infinite()) {
      return boost::optional<Number>();
    } else {
      return boost::optional<Number>(_n);
    }
  }

  void write(crab::crab_os &o) const {
    if (is_plus_infinity()) {
      o << "+oo";
    } else if (is_minus_infinity()) {
      o << "-oo";
    } else {
      o << _n;
    }
  }

}; // class bound

template <typename Number>
inline crab::crab_os &operator<<(crab::crab_os &o, const bound<Number> &b) {
  b.write(o);
  return o;
}

typedef bound<z_number> z_bound;
typedef bound<q_number> q_bound;

namespace bounds_impl {
// Conversion between z_bound and q_bound
// template<class B1, class B2>
// inline void convert_bounds(B1 b1, B2& b2);

inline void convert_bounds(z_bound b1, z_bound &b2) { std::swap(b1, b2); }
inline void convert_bounds(q_bound b1, q_bound &b2) { std::swap(b1, b2); }
inline void convert_bounds(z_bound b1, q_bound &b2) {
  if (b1.is_plus_infinity())
    b2 = q_bound::plus_infinity();
  else if (b1.is_minus_infinity())
    b2 = q_bound::minus_infinity();
  else
    b2 = q_bound(q_number(*b1.number()));
}
inline void convert_bounds(q_bound b1, z_bound &b2) {
  if (b1.is_plus_infinity())
    b2 = z_bound::plus_infinity();
  else if (b1.is_minus_infinity())
    b2 = z_bound::minus_infinity();
  else
    b2 = z_bound((*(b1.number())).round_to_lower());
}
} // namespace bounds_impl

template <typename Number> class interval {

public:
  typedef bound<Number> bound_t;
  typedef interval<Number> interval_t;

private:
  bound_t _lb;
  bound_t _ub;

public:
  static interval_t top() {
    return interval_t(bound_t::minus_infinity(), bound_t::plus_infinity());
  }

  static interval_t bottom() { return interval_t(); }

private:
  interval() : _lb(0), _ub(-1) {}

  static Number abs(Number x) { return x < 0 ? -x : x; }

  static Number max(Number x, Number y) { return x.operator<=(y) ? y : x; }

  static Number min(Number x, Number y) { return x.operator<(y) ? x : y; }

public:
  interval(bound_t lb, bound_t ub) : _lb(lb), _ub(ub) {
    if (lb > ub) {
      _lb = 0;
      _ub = -1;
    }
  }

  interval(bound_t b) : _lb(b), _ub(b) {
    if (b.is_infinite()) {
      _lb = 0;
      _ub = -1;
    }
  }

  interval(Number n) : _lb(n), _ub(n) {}

  interval(std::string b) : _lb(b), _ub(b) {
    if (_lb.is_infinite()) {
      _lb = 0;
      _ub = -1;
    }
  }

  interval(const interval_t &i) : _lb(i._lb), _ub(i._ub) {}

  interval_t &operator=(interval_t i) {
    _lb = i._lb;
    _ub = i._ub;
    return *this;
  }

  bound_t lb() const { return _lb; }

  bound_t ub() const { return _ub; }

  bool is_bottom() const { return (_lb > _ub); }

  bool is_top() const { return (_lb.is_infinite() && _ub.is_infinite()); }

  interval_t lower_half_line() const {
    return interval_t(bound_t::minus_infinity(), _ub);
  }

  interval_t upper_half_line() const {
    return interval_t(_lb, bound_t::plus_infinity());
  }

  bool operator==(interval_t x) const {
    if (is_bottom()) {
      return x.is_bottom();
    } else {
      return (_lb == x._lb) && (_ub == x._ub);
    }
  }

  bool operator!=(interval_t x) const { return !operator==(x); }

  bool operator<=(interval_t x) const {
    if (is_bottom()) {
      return true;
    } else if (x.is_bottom()) {
      return false;
    } else {
      return (x._lb <= _lb) && (_ub <= x._ub);
    }
  }

  interval_t operator|(interval_t x) const {
    if (is_bottom()) {
      return x;
    } else if (x.is_bottom()) {
      return *this;
    } else {
      return interval_t(bound_t::min(_lb, x._lb), bound_t::max(_ub, x._ub));
    }
  }

  interval_t operator&(interval_t x) const {
    if (is_bottom() || x.is_bottom()) {
      return bottom();
    } else {
      return interval_t(bound_t::max(_lb, x._lb), bound_t::min(_ub, x._ub));
    }
  }

  interval_t operator||(interval_t x) const {
    if (is_bottom()) {
      return x;
    } else if (x.is_bottom()) {
      return *this;
    } else {
      return interval_t(x._lb < _lb ? bound_t::minus_infinity() : _lb,
                        _ub < x._ub ? bound_t::plus_infinity() : _ub);
    }
  }

  template <typename Thresholds>
  interval_t widening_thresholds(interval_t x, const Thresholds &ts) {
    if (is_bottom()) {
      return x;
    } else if (x.is_bottom()) {
      return *this;
    } else {
      bound_t lb = (x._lb < _lb ? ts.get_prev(x._lb) : _lb);
      bound_t ub = (_ub < x._ub ? ts.get_next(x._ub) : _ub);
      return interval_t(lb, ub);
    }
  }

  interval_t operator&&(interval_t x) const {
    if (is_bottom() || x.is_bottom()) {
      return bottom();
    } else {
      return interval_t(_lb.is_infinite() && x._lb.is_finite() ? x._lb : _lb,
                        _ub.is_infinite() && x._ub.is_finite() ? x._ub : _ub);
    }
  }

  interval_t operator+(interval_t x) const {
    if (is_bottom() || x.is_bottom()) {
      return bottom();
    } else {
      return interval_t(_lb + x._lb, _ub + x._ub);
    }
  }

  interval_t &operator+=(interval_t x) { return operator=(operator+(x)); }

  interval_t operator-() const {
    if (is_bottom()) {
      return bottom();
    } else {
      return interval_t(-_ub, -_lb);
    }
  }

  interval_t operator-(interval_t x) const {
    if (is_bottom() || x.is_bottom()) {
      return bottom();
    } else {
      return interval_t(_lb - x._ub, _ub - x._lb);
    }
  }

  interval_t &operator-=(interval_t x) { return operator=(operator-(x)); }

  interval_t operator*(interval_t x) const {
    if (is_bottom() || x.is_bottom()) {
      return bottom();
    } else {
      bound_t ll = _lb * x._lb;
      bound_t lu = _lb * x._ub;
      bound_t ul = _ub * x._lb;
      bound_t uu = _ub * x._ub;
      return interval_t(bound_t::min(ll, lu, ul, uu),
                        bound_t::max(ll, lu, ul, uu));
    }
  }

  interval_t &operator*=(interval_t x) { return operator=(operator*(x)); }

  interval_t operator/(interval_t x) const;

  interval_t &operator/=(interval_t x) { return operator=(operator/(x)); }

  boost::optional<Number> singleton() const {
    if (!is_bottom() && _lb == _ub) {
      return _lb.number();
    } else {
      return boost::optional<Number>();
    }
  }

  bool operator[](Number n) const {
    if (is_bottom()) {
      return false;
    } else {
      bound_t b(n);
      return (_lb <= b) && (b <= _ub);
    }
  }

  void write(crab::crab_os &o) const {
    if (is_bottom()) {
      o << "_|_";
    } else {
      o << "[" << _lb << ", " << _ub << "]";
    }
  }

  // division and remainder operations

  interval_t UDiv(interval_t x) const {
    if (is_bottom() || x.is_bottom()) {
      return bottom();
    } else {
      return top();
    }
  }

  interval_t SRem(interval_t x) const {
    if (is_bottom() || x.is_bottom()) {
      return bottom();
    } else {
      return top();
    }
  }

  interval_t URem(interval_t x) const {
    if (is_bottom() || x.is_bottom()) {
      return bottom();
    } else {
      return top();
    }
  }

  // bitwise operations
  interval_t And(interval_t x) const {
    if (is_bottom() || x.is_bottom()) {
      return bottom();
    } else {
      return top();
    }
  }

  interval_t Or(interval_t x) const {
    if (is_bottom() || x.is_bottom()) {
      return bottom();
    } else {
      return top();
    }
  }

  interval_t Xor(interval_t x) const { return Or(x); }

  interval_t Shl(interval_t x) const {
    if (is_bottom() || x.is_bottom()) {
      return bottom();
    } else {
      return top();
    }
  }

  interval_t LShr(interval_t x) const {
    if (is_bottom() || x.is_bottom()) {
      return bottom();
    } else {
      return top();
    }
  }

  interval_t AShr(interval_t x) const {
    if (is_bottom() || x.is_bottom()) {
      return bottom();
    } else {
      return top();
    }
  }

}; //  class interval

template <>
inline interval<q_number> interval<q_number>::
operator/(interval<q_number> x) const {
  if (is_bottom() || x.is_bottom()) {
    return bottom();
  } else {
    boost::optional<q_number> d = x.singleton();
    if (d && *d == 0) {
      // [_, _] / 0 = _|_
      return bottom();
    } else if (x[0]) {
      boost::optional<q_number> n = singleton();
      if (n && *n == 0) {
        // 0 / [_, _] = 0
        return interval_t(q_number(0));
      } else {
        return top();
      }
    } else {
      bound_t ll = _lb / x._lb;
      bound_t lu = _lb / x._ub;
      bound_t ul = _ub / x._lb;
      bound_t uu = _ub / x._ub;
      return interval_t(bound_t::min(ll, lu, ul, uu),
                        bound_t::max(ll, lu, ul, uu));
    }
  }
}

template <>
inline interval<z_number> interval<z_number>::
operator/(interval<z_number> x) const {
  if (is_bottom() || x.is_bottom()) {
    return bottom();
  } else {
    // Divisor is a singleton:
    //   the linear interval solver can perform many divisions where
    //   the divisor is a singleton interval. We optimize for this case.
    if (boost::optional<z_number> n = x.singleton()) {
      z_number c = *n;
      if (c == 1) {
        return *this;
      } else if (c > 0) {
        return interval_t(_lb / c, _ub / c);
      } else if (c < 0) {
        return interval_t(_ub / c, _lb / c);
      } else {
      }
    }
    // Divisor is not a singleton
    typedef interval<z_number> z_interval;
    if (x[0]) {
      z_interval l(x._lb, z_bound(-1));
      z_interval u(z_bound(1), x._ub);
      return (operator/(l) | operator/(u));
    } else if (operator[](0)) {
      z_interval l(_lb, z_bound(-1));
      z_interval u(z_bound(1), _ub);
      return ((l / x) | (u / x) | z_interval(z_number(0)));
    } else {
      // Neither the dividend nor the divisor contains 0
      z_interval a =
          (_ub < 0) ? (*this + ((x._ub < 0) ? (x + z_interval(z_number(1)))
                                            : (z_interval(z_number(1)) - x)))
                    : *this;
      bound_t ll = a._lb / x._lb;
      bound_t lu = a._lb / x._ub;
      bound_t ul = a._ub / x._lb;
      bound_t uu = a._ub / x._ub;
      return interval_t(bound_t::min(ll, lu, ul, uu),
                        bound_t::max(ll, lu, ul, uu));
    }
  }
}

template <>
inline interval<z_number> interval<z_number>::SRem(interval<z_number> x) const {
  // note that the sign of the divisor does not matter

  if (is_bottom() || x.is_bottom()) {
    return bottom();
  } else if (singleton() && x.singleton()) {
    z_number dividend = *singleton();
    z_number divisor = *x.singleton();

    if (divisor == 0) {
      return bottom();
    }

    return interval_t(dividend % divisor);
  } else if (x.ub().is_finite() && x.lb().is_finite()) {
    z_number max_divisor = max(abs(*x.lb().number()), abs(*x.ub().number()));

    if (max_divisor == 0) {
      return bottom();
    }

    if (lb() < 0) {
      if (ub() > 0) {
        return interval_t(-(max_divisor - 1), max_divisor - 1);
      } else {
        return interval_t(-(max_divisor - 1), 0);
      }
    } else {
      return interval_t(0, max_divisor - 1);
    }
  } else {
    return top();
  }
}

template <>
inline interval<z_number> interval<z_number>::URem(interval<z_number> x) const {

  if (is_bottom() || x.is_bottom()) {
    return bottom();
  } else if (singleton() && x.singleton()) {
    z_number dividend = *singleton();
    z_number divisor = *x.singleton();

    if (divisor < 0) {
      return top();
    } else if (divisor == 0) {
      return bottom();
    } else if (dividend < 0) {
      // dividend is treated as an unsigned integer.
      // we would need the size to be more precise
      return interval_t(0, divisor - 1);
    } else {
      return interval_t(dividend % divisor);
    }
  } else if (x.ub().is_finite() && x.lb().is_finite()) {
    z_number max_divisor = *x.ub().number();

    if (x.lb() < 0 || x.ub() < 0) {
      return top();
    } else if (max_divisor == 0) {
      return bottom();
    }

    return interval_t(0, max_divisor - 1);
  } else {
    return top();
  }
}

template <>
inline interval<z_number> interval<z_number>::And(interval<z_number> x) const {
  if (is_bottom() || x.is_bottom()) {
    return bottom();
  } else {
    boost::optional<z_number> left_op = singleton();
    boost::optional<z_number> right_op = x.singleton();

    if (left_op && right_op) {
      return interval_t((*left_op) & (*right_op));
    } else if (lb() >= 0 && x.lb() >= 0) {
      return interval_t(0, bound_t::min(ub(), x.ub()));
    } else {
      return top();
    }
  }
}

template <>
inline interval<z_number> interval<z_number>::Or(interval<z_number> x) const {
  if (is_bottom() || x.is_bottom()) {
    return bottom();
  } else {
    boost::optional<z_number> left_op = singleton();
    boost::optional<z_number> right_op = x.singleton();

    if (left_op && right_op) {
      return interval_t((*left_op) | (*right_op));
    } else if (lb() >= 0 && x.lb() >= 0) {
      boost::optional<z_number> left_ub = ub().number();
      boost::optional<z_number> right_ub = x.ub().number();

      if (left_ub && right_ub) {
        z_number m = (*left_ub > *right_ub ? *left_ub : *right_ub);
        return interval_t(0, m.fill_ones());
      } else {
        return interval_t(0, bound_t::plus_infinity());
      }
    } else {
      return top();
    }
  }
}

template <>
inline interval<z_number> interval<z_number>::Xor(interval<z_number> x) const {
  if (is_bottom() || x.is_bottom()) {
    return bottom();
  } else {
    boost::optional<z_number> left_op = singleton();
    boost::optional<z_number> right_op = x.singleton();

    if (left_op && right_op) {
      return interval_t((*left_op) ^ (*right_op));
    } else {
      return Or(x);
    }
  }
}

template <>
inline interval<z_number> interval<z_number>::Shl(interval<z_number> x) const {
  if (is_bottom() || x.is_bottom()) {
    return bottom();
  } else {
    if (boost::optional<z_number> shift = x.singleton()) {
      z_number k = *shift;
      if (k < 0) {
        // CRAB_ERROR("lshr shift operand cannot be negative");
        return top();
      }
      // Some crazy linux drivers generate shl instructions with
      // huge shifts.  We limit the number of times the loop is run
      // to avoid wasting too much time on it.
      if (k <= 128) {
        z_number factor = 1;
        for (int i = 0; k > i; i++) {
          factor *= 2;
        }
        return (*this) * factor;
      }
    }
    return top();
  }
}

template <>
inline interval<z_number> interval<z_number>::AShr(interval<z_number> x) const {
  if (is_bottom() || x.is_bottom()) {
    return bottom();
  } else {
    if (boost::optional<z_number> shift = x.singleton()) {
      z_number k = *shift;
      if (k < 0) {
        return top();
      }
      // Some crazy linux drivers generate ashr instructions with
      // huge shifts.  We limit the number of times the loop is run
      // to avoid wasting too much time on it.
      if (k <= 128) {
        z_number factor = 1;
        for (int i = 0; k > i; i++) {
          factor *= 2;
        }
        return (*this) / factor;
      }
    }
    return top();
  }
}

template <>
inline interval<z_number> interval<z_number>::LShr(interval<z_number> x) const {
  if (is_bottom() || x.is_bottom()) {
    return bottom();
  } else {
    if (boost::optional<z_number> shift = x.singleton()) {
      z_number k = *shift;
      if (k < 0) {
        return top();
      }
      if (lb() >= 0 && ub().is_finite()) {
        z_number lb = *this->lb().number();
        z_number ub = *this->ub().number();
        return interval<z_number>(lb >> k, ub >> k);
      }
    }
    return this->top();
  }
}

template <typename Number>
inline interval<Number> operator+(Number c, interval<Number> x) {
  return interval<Number>(c) + x;
}

template <typename Number>
inline interval<Number> operator+(interval<Number> x, Number c) {
  return x + interval<Number>(c);
}

template <typename Number>
inline interval<Number> operator*(Number c, interval<Number> x) {
  return interval<Number>(c) * x;
}

template <typename Number>
inline interval<Number> operator*(interval<Number> x, Number c) {
  return x * interval<Number>(c);
}

template <typename Number>
inline interval<Number> operator/(Number c, interval<Number> x) {
  return interval<Number>(c) / x;
}

template <typename Number>
inline interval<Number> operator/(interval<Number> x, Number c) {
  return x / interval<Number>(c);
}

template <typename Number>
inline interval<Number> operator-(Number c, interval<Number> x) {
  return interval<Number>(c) - x;
}

template <typename Number>
inline interval<Number> operator-(interval<Number> x, Number c) {
  return x - interval<Number>(c);
}

template <typename Number>
inline crab::crab_os &operator<<(crab::crab_os &o, const interval<Number> &i) {
  i.write(o);
  return o;
}

namespace linear_interval_solver_impl {

typedef interval<z_number> z_interval;
typedef interval<q_number> q_interval;

template <> inline z_interval trim_interval(z_interval i, z_interval j) {
  if (boost::optional<z_number> c = j.singleton()) {
    if (i.lb() == *c) {
      return z_interval(*c + 1, i.ub());
    } else if (i.ub() == *c) {
      return z_interval(i.lb(), *c - 1);
    } else {
    }
  }
  return i;
}

template <> inline q_interval trim_interval(q_interval i, q_interval /* j */) {
  // No refinement possible for disequations over rational numbers
  return i;
}

template <>
inline z_interval lower_half_line(z_interval i, bool /*is_signed*/) {
  return i.lower_half_line();
}

template <>
inline q_interval lower_half_line(q_interval i, bool /*is_signed*/) {
  return i.lower_half_line();
}

template <>
inline z_interval upper_half_line(z_interval i, bool /*is_signed*/) {
  return i.upper_half_line();
}

template <>
inline q_interval upper_half_line(q_interval i, bool /*is_signed*/) {
  return i.upper_half_line();
}
} // namespace linear_interval_solver_impl

} // namespace ikos

ext: hpp

max_stars_repo_name: yugeshk/crab

max_issues_repo_name: yugeshk/crab

max_forks_repo_name: yugeshk/crab


text:
/* Copyright Institute of Sound and Vibration Research - All rights reserved */

#include "options.hpp"
#include "signal_flow.hpp"

#include <libefl/denormalised_number_handling.hpp>

#include <librrl/audio_signal_flow.hpp>

#include <libvisr/version.hpp>

#include <libaudiointerfaces/audio_interface_factory.hpp>
#ifdef VISR_JACK_SUPPORT
#include <libaudiointerfaces/jack_interface.hpp>
#endif
#include <libaudiointerfaces/portaudio_interface.hpp>
#include <libaudiointerfaces/audio_interface.hpp>

#include <boost/algorithm/string.hpp> // case-insensitive string compare
#include <boost/filesystem.hpp>

#include <algorithm>
#include <cstddef>
#include <cstdlib>
#include <cstdio> // for getc(), for testing purposes
#include <iostream>
#include <memory>
#include <sstream>

int main( int argc, char const * const * argv )
{
    using namespace visr;
    using namespace visr::apps::audio_network_streamer;

    try
    {
        efl::DenormalisedNumbers::State const oldDenormNumbersState
        = efl::DenormalisedNumbers::setDenormHandling();

        Options cmdLineOptions;
        std::stringstream errMsg;
        switch (cmdLineOptions.parse(argc, argv, errMsg))
        {
            case Options::ParseResult::Failure:
                std::cerr << "Error while parsing command line options: " << errMsg.str() << std::endl;
                return EXIT_FAILURE;
            case Options::ParseResult::Help:
                cmdLineOptions.printDescription(std::cout);
                return EXIT_SUCCESS;
            case Options::ParseResult::Version:
                // TODO: Outsource the version string generation to a central file.
                std::cout << "VISR S3A Audio over UDP streaming "
                << visr::version::versionString() << std::endl;
                return EXIT_SUCCESS;
            case Options::ParseResult::Success:
                break; // carry on
        }

        std::string const audioBackend = cmdLineOptions.getDefaultedOption<std::string>("audio-backend", "default");

        SamplingFrequencyType const samplingFrequency = cmdLineOptions.getDefaultedOption<SamplingFrequencyType>( "sampling-frequency", 48000 );

        std::string const sceneSendAddresses = cmdLineOptions.getOption<std::string>("send-addresses" );

        std::size_t const blockSize = cmdLineOptions.getDefaultedOption<std::size_t>("block-size", 1024);

        // Selection of audio interface:
        // For the moment we check for the name 'NATIVE_JACK' and select the specialized audio interface and fall
        // back to PortAudio in all other cases.
        // TODO: Provide factory and backend-specific options) to make selection of audio interfaces more general and extendable.
#ifdef VISR_JACK_SUPPORT
        bool const useNativeJack = boost::iequals(audioBackend, "NATIVE_JACK");
#endif

        SignalFlowContext const context{ blockSize, samplingFrequency };

        SignalFlow topLevelGraph( context, "NetworkAudioStreamer", nullptr, sceneSendAddresses );

        rrl::AudioSignalFlow flow( topLevelGraph );

        std::size_t const numberOfSignals = flow.numberOfCaptureChannels();

        audiointerfaces::AudioInterface::Configuration baseConfig(numberOfSignals,0,samplingFrequency,blockSize);

        std::string type;
        std::string specConf;

        std::unique_ptr<visr::audiointerfaces::AudioInterface> audioInterface;

#ifdef VISR_JACK_SUPPORT
        if (useNativeJack)
        {
            std::string pconfig = "{\"ports\":[ \n { \"captbasename\": \"input_\"}, \n ]}";
            specConf = "{\"clientname\": \"S3A network audio streamer\", \"servername\": \"\", \"portsconfig\" : "+pconfig+"}";

            type = "Jack";
        }
        else
        {
#endif
            specConf = "{\"sampleformat\": 8, \"interleaved\": \"false\", \"hostapi\" : "+audioBackend+"}";
            type = "PortAudio";

#ifdef VISR_JACK_SUPPORT
        }
#endif

        audioInterface.reset( audiointerfaces::AudioInterfaceFactory::create( type, baseConfig, specConf).get());
        audioInterface->registerCallback( &rrl::AudioSignalFlow::processFunction, &flow );

        // should there be a separate start() method for the audio interface?
        audioInterface->start( );

        // Rendering runs until q<Return> is entered on the console.
        std::cout << "S3A network audio streamer running. Press \"q<Return>\" or Ctrl-C to quit." << std::endl;
        char c;
        do
        {
            c = std::getc( stdin );
        }
        while( c != 'q' );

        audioInterface->stop( );

        audioInterface->unregisterCallback( &rrl::AudioSignalFlow::processFunction );

        efl::DenormalisedNumbers::resetDenormHandling( oldDenormNumbersState );
    }
    catch( std::exception const & ex )
    {
        std::cout << "Exception caught on top level: " << ex.what() << std::endl;
        return EXIT_FAILURE;
    }

    return EXIT_SUCCESS;
}

ext: cpp


text:
// Software License for MTL
//
// Copyright (c) 2007 The Trustees of Indiana University.
//               2008 Dresden University of Technology and the Trustees of Indiana University.
//               2010 SimuNova UG (haftungsbeschränkt), www.simunova.com.
// All rights reserved.
// Authors: Peter Gottschling and Andrew Lumsdaine
//
// This file is part of the Matrix Template Library
//
// See also license.mtl.txt in the distribution.

#ifndef MTL_TWO_NORM_INCLUDE
#define MTL_TWO_NORM_INCLUDE

#include <iostream>
#include <cmath>

#include <boost/numeric/mtl/concept/collection.hpp>
#include <boost/numeric/mtl/concept/magnitude.hpp>
#include <boost/numeric/mtl/utility/tag.hpp>
#include <boost/numeric/mtl/utility/category.hpp>
#include <boost/numeric/mtl/vector/lazy_reduction.hpp>
#include <boost/numeric/mtl/vector/reduction.hpp>
#include <boost/numeric/mtl/vector/reduction_functors.hpp>
#include <boost/numeric/mtl/interface/vpt.hpp>

namespace mtl {

    namespace vec {

	template <unsigned long Unroll, typename Value>
	typename RealMagnitude<typename Collection<Value>::value_type>::type
	inline two_norm(const Value& value)
	{
	    using std::sqrt;
	    vampir_trace<2039> tracer;
	    typedef typename RealMagnitude<typename Collection<Value>::value_type>::type result_type;
	    return sqrt(reduction<Unroll, two_norm_functor, result_type>::apply(value));
	}

	/*! Two-norm for vectors: two_norm(x) \f$\rightarrow |x|_2\f$.
	    \retval The magnitude type of the respective value type, see Magnitude.
	    The norms are defined as \f$|v|_2=\sqrt{\sum_i |v_i|^2}\f$.

	    Vector norms are unrolled 8-fold by default.
	    An n-fold unrolling can be generated with two_norm<n>(x).
	    The maximum for n is 8 (it might be increased later).
	**/
	template <typename Value>
	typename RealMagnitude<typename Collection<Value>::value_type>::type
	inline two_norm(const Value& value)
	{
	    return two_norm<4>(value);
	}

	template <typename Vector>
	lazy_reduction<Vector, two_norm_functor> inline lazy_two_norm(const Vector& v)
	{  return lazy_reduction<Vector, two_norm_functor>(v); 	}

    } // namespace vector

    // two_norm for matrices not implemented (would need enable_if like one_norm)

    using vec::two_norm;
    using vec::lazy_two_norm;

} // namespace mtl

#endif // MTL_TWO_NORM_INCLUDE


text: #include <boost/xpressive/sub_match.hpp>

max_stars_repo_name: miathedev/BoostForArduino

max_issues_repo_name: miathedev/BoostForArduino

max_forks_repo_name: miathedev/BoostForArduino


text:
//
// $Id: SpectrumList_PeakFilter.hpp 6902 2014-11-17 22:34:37Z chambm $
//
//
// Original author: Matt Chambers <matt.chambers <a.t> vanderbilt.edu>
//
// Copyright 2008 Vanderbilt University - Nashville, TN 37232
//
// Licensed under the Apache License, Version 2.0 (the "License");
// you may not use this file except in compliance with the License.
// You may obtain a copy of the License at
//
// http://www.apache.org/licenses/LICENSE-2.0
//
// Unless required by applicable law or agreed to in writing, software
// distributed under the License is distributed on an "AS IS" BASIS,
// WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
// See the License for the specific language governing permissions and
// limitations under the License.
//

#ifndef _SPECTRUMLIST_PEAKFILTER_HPP_
#define _SPECTRUMLIST_PEAKFILTER_HPP_

#include <boost/scoped_ptr.hpp>
#include "pwiz/data/msdata/SpectrumListWrapper.hpp"
#include "pwiz/analysis/common/DataFilter.hpp"

namespace pwiz {
namespace analysis {

/// SpectrumList implementation that returns spectra with the specified SpectrumDataFilter operation applied
class PWIZ_API_DECL SpectrumList_PeakFilter : public msdata::SpectrumListWrapper
{
    public:

    SpectrumList_PeakFilter(const msdata::SpectrumListPtr& inner,
                            SpectrumDataFilterPtr filterFunctor);

    /// peak filters work on any SpectrumList
    static bool accept(const msdata::SpectrumListPtr& inner) {return true;}

    virtual msdata::SpectrumPtr spectrum(size_t index, bool getBinaryData = false) const;
    virtual msdata::SpectrumPtr spectrum(size_t index, msdata::DetailLevel detailLevel = msdata::DetailLevel_FullMetadata) const;

    private:
	SpectrumDataFilterPtr filterFunctor_;
};

} // namespace analysis
} // namespace pwiz

#endif // _SPECTRUMLIST_PEAKFILTER_HPP_

max_stars_repo_name: shze/pwizard-deb

max_issues_repo_name: shze/pwizard-deb

max_forks_repo_name: shze/pwizard-deb


text:
//============================================================================
// Name        : Request
// Author      : Avi
// Revision    : $Revision: #32 $
//
// Copyright 2009-2020 ECMWF.
// This software is licensed under the terms of the Apache Licence version 2.0
// which can be obtained at http://www.apache.org/licenses/LICENSE-2.0.
// In applying this licence, ECMWF does not waive the privileges and immunities
// granted to it by virtue of its status as an intergovernmental organisation
// nor does it submit to any jurisdiction.
//
// Description :
//============================================================================
#include <string>
#include <iostream>

#include <boost/test/unit_test.hpp>
#include <boost/date_time/posix_time/time_formatters.hpp>

#include "TimeSeries.hpp"
#include "Calendar.hpp"
#include "Str.hpp"

using namespace std;
using namespace ecf;
using namespace boost::posix_time;
using namespace boost::gregorian;

using namespace boost;

BOOST_AUTO_TEST_SUITE( CoreTestSuite )

BOOST_AUTO_TEST_CASE( test_default_boost_time_duration )
{
   cout << "ACore:: ...test_default_boost_time_duration\n";

   boost::posix_time::time_duration td;
   BOOST_CHECK_MESSAGE( td.is_special() == false, "Default constructed boost::posix_time::time_duration, should not be special");
   BOOST_CHECK_MESSAGE( td.hours() == 0, "Expected 0 for hours");
   BOOST_CHECK_MESSAGE( td.minutes() == 0, "Expected 0 for minutes");
   BOOST_CHECK_MESSAGE( td.seconds() == 0, "Expected 0 for seconds ");
}

BOOST_AUTO_TEST_CASE( test_time_series_min_max_time_slots )
{
   cout << "ACore:: ...test_time_series_min_max_time_slots\n";

   TimeSlot the_min;
   TimeSlot the_max;

   TimeSeries a(9,0,true);
   a.min_max_time_slots(the_min,the_max);
   BOOST_CHECK_MESSAGE( the_min == TimeSlot(9,0),"Max min time slot failed");
   BOOST_CHECK_MESSAGE( the_max == TimeSlot(9,0),"Max min time slot failed");

   TimeSeries b(TimeSlot(10,12),false);
   b.min_max_time_slots(the_min,the_max);
   BOOST_CHECK_MESSAGE( the_min == TimeSlot(9,0),"Max min time slot failed");
   BOOST_CHECK_MESSAGE( the_max == TimeSlot(10,12),"Max min time slot failed");

   TimeSeries c(TimeSlot(0,10), TimeSlot(10,13), TimeSlot(0,1));
   c.min_max_time_slots(the_min,the_max);
   BOOST_CHECK_MESSAGE( the_min == TimeSlot(0,10),"Max min time slot failed");
   BOOST_CHECK_MESSAGE( the_max == TimeSlot(10,13),"Max min time slot failed");

   TimeSeries x;
   x.min_max_time_slots(the_min,the_max);
   BOOST_CHECK_MESSAGE( the_min == TimeSlot(),"Max min time slot failed");
   BOOST_CHECK_MESSAGE( the_max == TimeSlot(10,13),"Max min time slot failed");

}

BOOST_AUTO_TEST_CASE( test_time_series_constrcution )
{
	cout << "ACore:: ...test_time_series_constrcution\n";

	{
		TimeSeries x;
		TimeSeries y;
		BOOST_CHECK_MESSAGE( x == y,"Equality operator expected to succeed");

		TimeSeries a(10,12,true);
		TimeSeries b(10,12,true);
		BOOST_CHECK_MESSAGE( a == b,"Equality operator expected to succeed");

		TimeSeries c(10,12,false);
		TimeSeries d(10,12,false);
		BOOST_CHECK_MESSAGE( c == d,"Equality operator expected to succeed");

		TimeSeries e(TimeSlot(14,59),false);
		TimeSeries f(TimeSlot(14,59),false);
 		BOOST_CHECK_MESSAGE( e == f,"Equality operator expected to succeed");

 		TimeSeries g(TimeSlot(0,10), TimeSlot(10,4), TimeSlot(0,1));
 		TimeSeries h(TimeSlot(0,10), TimeSlot(10,4), TimeSlot(0,1));
 		BOOST_CHECK_MESSAGE( g == h,"Equality operator expected to succeed");
	}
	{
		TimeSeries a(10,12,false);
		TimeSeries b(10,12,true);
 		BOOST_CHECK_MESSAGE( a != b,"Equality operator expected to fail");

		TimeSeries c(10,13,false);
		TimeSeries d(10,12,false);
		BOOST_CHECK_MESSAGE( c != d,"Equality operator expected to fail");

		TimeSeries e(TimeSlot(14,59),false);
		TimeSeries f(TimeSlot(14,10),false);
 		BOOST_CHECK_MESSAGE( e != f,"Equality operator expected to fail");

 		TimeSeries g(TimeSlot(0,10), TimeSlot(10,4), TimeSlot(0,1), true);
 		TimeSeries h(TimeSlot(0,10), TimeSlot(10,4), TimeSlot(0,1), false);
 		BOOST_CHECK_MESSAGE( g != h,"Equality operator expected to fail");
 	}

	/// Basic test for time series of getters and setters
	TimeSeries timeSeries(TimeSlot(0,1), TimeSlot(0,4), TimeSlot(0,1));
	BOOST_CHECK_MESSAGE(timeSeries.hasIncrement(),"expected increment");
	std::string errMsg; BOOST_CHECK_MESSAGE(timeSeries.checkInvariants(errMsg),errMsg);

	TimeSeries another = timeSeries;
	BOOST_CHECK_MESSAGE(another == timeSeries,"copy constructor failed for TimeSeries ");
	BOOST_CHECK_MESSAGE(another.checkInvariants(errMsg),errMsg);
	BOOST_CHECK_MESSAGE(timeSeries.checkInvariants(errMsg),errMsg);

	TimeSeries timeSeries2(TimeSlot(0,1), TimeSlot(0,4), TimeSlot(0,1));
	BOOST_CHECK_MESSAGE(timeSeries2.checkInvariants(errMsg),errMsg);
	another = timeSeries2;
	BOOST_CHECK_MESSAGE(another ==  timeSeries2,"assignment operator failed for TimeSeries ");
	BOOST_CHECK_MESSAGE(another.checkInvariants(errMsg),errMsg);

	TimeSeries timeSeries3(TimeSlot(0,1));
	BOOST_CHECK_MESSAGE(timeSeries3.checkInvariants(errMsg),errMsg);
}

BOOST_AUTO_TEST_CASE( test_time_series )
{
	cout << "ACore:: ...test_time_series\n";

	/// Basic test for time series of getters and setters
	TimeSeries timeSeries(TimeSlot(0,1), TimeSlot(0,4), TimeSlot(0,1));
	BOOST_CHECK_MESSAGE(timeSeries.hasIncrement(),"expected increment");

	std::string expected = "00:01:00";
	std::string actual = to_simple_string(timeSeries.start().duration());
	BOOST_CHECK_MESSAGE(expected == actual," expected " << expected << " but found " << actual );

	expected = "00:04:00";
	actual = to_simple_string(timeSeries.finish().duration());
	BOOST_CHECK_MESSAGE(expected == actual," expected " << expected << " but found " << actual );

	expected = "00:01:00";
	actual = to_simple_string(timeSeries.incr().duration());
	BOOST_CHECK_MESSAGE(expected == actual," expected " << expected << " but found " << actual );

	TimeSeries timeSeries2(TimeSlot(0,1));
	BOOST_CHECK_MESSAGE(!timeSeries2.hasIncrement(),"Not expecting time series");
}

BOOST_AUTO_TEST_CASE( test_time_series_increment_real )
{
	cout << "ACore:: ...test_time_series_increment_real\n";

	// Test time series with  a calendar, we update calendar then
	// test time series isFree(), and checkForRequeue
	Calendar c;
	c.init(ptime(date(2010,2,10), minutes(0)), Calendar::REAL);

	// Create a test when we can match a time series. Need to sync hour with suite time
	// at hour 1, suite time should also be 01:00, for test to work
	//
	// Create the time series: start  10:00
	//                         finish 20:00
	//                         incr    1:00
	TimeSeries timeSeries(TimeSlot(10,0), TimeSlot(20,0), TimeSlot(1,0), true/* relative */);
	TimeSeries timeSeries2(TimeSlot(11,0), TimeSlot(15,0), TimeSlot(1,0), true/* relative */);
	TimeSeries timeSeries3(TimeSlot(15,0),  true/* relative */);

   TimeSlot t1_min, t1_max,t2_min,t2_max,t3_min,t3_max;
   timeSeries.min_max_time_slots(t1_min, t1_max);
   timeSeries2.min_max_time_slots(t2_min, t2_max);
   timeSeries3.min_max_time_slots(t3_min, t3_max);
   BOOST_CHECK_MESSAGE(t1_min == TimeSlot(10,0) && t1_max == TimeSlot(20,0),"Not as expected");
   BOOST_CHECK_MESSAGE(t2_min == TimeSlot(11,0) && t2_max == TimeSlot(15,0),"Not as expected");
   BOOST_CHECK_MESSAGE(t3_min == TimeSlot(15,0) && t3_max == TimeSlot(15,0),"Not as expected");

   // Follow normal process
   timeSeries.reset( c );
   timeSeries2.reset( c );
   timeSeries3.reset( c );

   bool cmd_context = true;
	for(int hour=1; hour < 24; hour++) {
		// Update calendar every hour, then see we can match time series, *RELATIVE* to suite start
		c.update( time_duration( hours(1) ) );
		timeSeries.calendarChanged( c );
		timeSeries2.calendarChanged( c );
		timeSeries3.calendarChanged( c );

//		cerr << "hour = " << hour << " calendar_duration " << to_simple_string(timeSeries.duration(c))
//		    << " timeSeries=" << timeSeries.toString() << " timeSeries2=" << timeSeries2.toString() << " timeSeries3=" << timeSeries3.toString() << "\n";
		if (hour < timeSeries.start().hour()) {
	         BOOST_CHECK_MESSAGE(timeSeries.checkForRequeue(c,t1_min,t1_max,cmd_context)," Time series " << timeSeries.toString() << " checkForRequeue should pass at " << c.suite_time_str() );
	         BOOST_CHECK_MESSAGE(!timeSeries.checkForRequeue(c,t1_min,t1_max)," Time series " << timeSeries.toString() << " checkForRequeue should fail at " << c.suite_time_str() );
             BOOST_CHECK_MESSAGE(!timeSeries.isFree(c),"Time series " << timeSeries.toString() << " should NOT be free at hour " << hour );
		}
		else if (hour >= timeSeries.start().hour() && hour <= timeSeries.finish().hour()) {
			BOOST_CHECK_MESSAGE(timeSeries.isFree(c),"Time series " << timeSeries.toString() << " should be free at hour " << hour );

			/// At the last hour checkForRequeue should return false; This ensures that value will
			/// not get incremented and so, should leave node in the complete state.
			if ( hour < timeSeries.finish().hour()) {
				BOOST_CHECK_MESSAGE(timeSeries.checkForRequeue(c,t1_min,t1_max),"Time series " << timeSeries.toString() << " checkForRequeue should be free at hour " << hour );
			}
			else {
				BOOST_CHECK_MESSAGE(!timeSeries.checkForRequeue(c,t1_min,t1_max),"Time series " << timeSeries.toString() << "checkForRequeue should Not free at hour " << hour );
			}
		}
		else {
			BOOST_CHECK_MESSAGE(!timeSeries.isFree(c),"Time series " << timeSeries.toString() << " should NOT be free at hour " << hour );
 			BOOST_CHECK_MESSAGE(!timeSeries.checkForRequeue(c,t1_min,t1_max)," Time series " << timeSeries.toString() << " should fail at " << hour );
		}

      if (hour < timeSeries2.start().hour()) {
         BOOST_CHECK_MESSAGE(!timeSeries2.checkForRequeue(c,t1_min,t1_max)," Time series " << timeSeries2.toString() << "checkForRequeue should pass at " << hour );
         BOOST_CHECK_MESSAGE(!timeSeries2.isFree(c),"Time series " << timeSeries2.toString() << " should NOT be free at hour " << hour );
      }
      else if (hour >= timeSeries2.start().hour() && hour <=timeSeries2.finish().hour()) {
			BOOST_CHECK_MESSAGE(timeSeries2.isFree(c),"Time series " << timeSeries2.toString() << " should be free at hour " << hour );

			/// At the last hour checkForRequeue should return false;
			if ( hour < timeSeries2.finish().hour()) {
				BOOST_CHECK_MESSAGE(timeSeries2.checkForRequeue(c,t2_min,t2_max),"Time series " << timeSeries2.toString() << " checkForRequeue should be free at hour " << hour );
			}
			else {
				BOOST_CHECK_MESSAGE(!timeSeries2.checkForRequeue(c,t2_min,t2_max),"Time series " << timeSeries2.toString() << "checkForRequeue should Not free at hour " << hour );
			}
  		}
		else {
			BOOST_CHECK_MESSAGE(!timeSeries2.isFree(c),"Time series " << timeSeries2.toString() << " not be free at hour " << hour );
			BOOST_CHECK_MESSAGE(!timeSeries2.checkForRequeue(c,t2_min,t2_max)," Time series " << timeSeries2.toString() << " should fail at " << hour );
 		}

		if (hour == timeSeries3.start().hour()  ) {
			BOOST_CHECK_MESSAGE(timeSeries3.isFree(c),"Time series " << timeSeries3.toString() << " should be free at hour " << hour );
 		}
		else if (hour < timeSeries3.start().hour()) {
			BOOST_CHECK_MESSAGE(!timeSeries3.isFree(c),"Time series " << timeSeries3.toString() << " isFree should fail at hour " << hour );
 		}
		else if (hour > timeSeries3.start().hour()) {
			BOOST_CHECK_MESSAGE(!timeSeries3.isFree(c),"Time series " << timeSeries3.toString() << " is Free should fail at hour " << hour );
	 	}
		BOOST_CHECK_MESSAGE(!timeSeries3.checkForRequeue(c,t3_min,t3_max)," Time series " << timeSeries3.toString() << " checkForRequeue should fail at " << hour );
	}
}

BOOST_AUTO_TEST_CASE( test_time_series_requeueable_and_compute_next_time_slot )
{
   cout << "ACore:: ...test_time_series_requeueable_and_compute_next_time_slot\n";

   // Test time series with  a calendar, we update calendar then
   // test time series requeueable(), and compute_next_time_slot
   // This are used with the WHY command
   Calendar c;
   c.init(ptime(date(2010,2,10), minutes(0)), Calendar::REAL);

   // Create a test when we can match a time series. Need to sync hour with suite time
   // at hour 1, suite time should also be 01:00, for test to work
   //
   // Create the time series: start  10:00
   //                         finish 20:00
   //                         incr    1:00
   TimeSeries timeSeries(TimeSlot(10,0), TimeSlot(20,0), TimeSlot(1,0), true/* relative */);
   TimeSeries timeSeries2(TimeSlot(11,0), TimeSlot(15,0), TimeSlot(1,0), true/* relative */);
   TimeSeries timeSeries3(TimeSlot(15,0),  true/* relative */);

   for(int hour=1; hour < 24; hour++) {
      // Update calendar every hour, then see we can match time series,
      c.update( time_duration( hours(1) ) );
      timeSeries.calendarChanged( c );
      timeSeries2.calendarChanged( c );
      timeSeries3.calendarChanged( c );

//    cerr << "hour = " << hour << " calendar_duration " << to_simple_string(timeSeries.duration(c))
//        << " timeSeries=" << timeSeries.toString() << " timeSeries2=" << timeSeries2.toString() << " timeSeries3=" << timeSeries3.toString() << "\n";
      if (hour < timeSeries.start().hour()) {
         TimeSlot next_time_slot = timeSeries.compute_next_time_slot(c);
         TimeSlot expected(10,0);
         BOOST_CHECK_MESSAGE(next_time_slot ==expected," Time series " << timeSeries.toString() << " at " << hour << " expected next time slot at " << expected.toString() << " but found " << next_time_slot.toString());
         BOOST_CHECK_MESSAGE(timeSeries.requeueable(c),"Time series " << timeSeries.toString() << " should be requeueable at hour " << hour );
      }
      else if (hour == timeSeries.start().hour() ) {
         TimeSlot next_time_slot = timeSeries.compute_next_time_slot(c);
         TimeSlot expected(11,0);
         BOOST_CHECK_MESSAGE(next_time_slot == expected," Time series " << timeSeries.toString() << " at " << hour << " expected next time slot at " << expected.toString() << " but found " << next_time_slot.toString());
         BOOST_CHECK_MESSAGE(timeSeries.requeueable(c),"Time series " << timeSeries.toString() << " should be requeueable at hour " << hour );
      }
      else if (hour > timeSeries.start().hour() &&  hour < timeSeries.finish().hour()) {
         TimeSlot next_time_slot = timeSeries.compute_next_time_slot(c);
         TimeSlot expected(hour+1,0);
         BOOST_CHECK_MESSAGE(next_time_slot == expected," Time series " << timeSeries.toString() << " at " << hour << " expected next time slot at " << expected.toString() << " but found " << next_time_slot.toString());
         BOOST_CHECK_MESSAGE(timeSeries.requeueable(c),"Time series " << timeSeries.toString() << " should be requeueable at hour " << hour );
      }
      else if (hour == timeSeries.finish().hour()) {
         TimeSlot next_time_slot = timeSeries.compute_next_time_slot(c);
         BOOST_CHECK_MESSAGE(next_time_slot.isNULL()," Time series " << timeSeries.toString() << " at " << hour << " expected next time slot  to be NULL");
         BOOST_CHECK_MESSAGE(!timeSeries.requeueable(c),"Time series " << timeSeries.toString() << " should NOT be requeueable at hour " << hour );
      }
      else if (hour > timeSeries.finish().hour()) {
         TimeSlot next_time_slot = timeSeries.compute_next_time_slot(c);
         BOOST_CHECK_MESSAGE(next_time_slot.isNULL()," Time series " << timeSeries.toString() << " at " << hour << " expected next time slot to be NULL");
         BOOST_CHECK_MESSAGE(!timeSeries.requeueable(c),"Time series " << timeSeries.toString() << " should NOT be requeueable at hour " << hour );
      }

      if (hour < timeSeries2.start().hour()) {
         TimeSlot next_time_slot = timeSeries2.compute_next_time_slot(c);
         TimeSlot expected(11,0);
         BOOST_CHECK_MESSAGE(next_time_slot ==expected," Time series " << timeSeries2.toString() << " at " << hour << " expected next time slot at " << expected.toString() << " but found " << next_time_slot.toString());
         BOOST_CHECK_MESSAGE(timeSeries2.requeueable(c),"Time series " << timeSeries2.toString() << " should be requeueable at hour " << hour );
      }
      else if (hour == timeSeries2.start().hour() ) {
         TimeSlot next_time_slot = timeSeries2.compute_next_time_slot(c);
         TimeSlot expected(12,0);
         BOOST_CHECK_MESSAGE(next_time_slot == expected," Time series " << timeSeries2.toString() << " at " << hour << " expected next time slot at " << expected.toString() << " but found " << next_time_slot.toString());
         BOOST_CHECK_MESSAGE(timeSeries2.requeueable(c),"Time series " << timeSeries2.toString() << " should NOT be requeueable at hour " << hour );
      }
      else if (hour > timeSeries2.start().hour() &&  hour < timeSeries2.finish().hour()) {
         TimeSlot next_time_slot = timeSeries2.compute_next_time_slot(c);
         TimeSlot expected(hour+1,0);
         BOOST_CHECK_MESSAGE(next_time_slot == expected," Time series " << timeSeries2.toString() << " at " << hour << " expected next time slot at " << expected.toString() << " but found " << next_time_slot.toString());
         BOOST_CHECK_MESSAGE(timeSeries2.requeueable(c),"Time series " << timeSeries2.toString() << " should be requeueable at hour " << hour );
      }
      else if (hour == timeSeries2.finish().hour()) {
         TimeSlot next_time_slot = timeSeries2.compute_next_time_slot(c);
         BOOST_CHECK_MESSAGE(next_time_slot.isNULL()," Time series " << timeSeries2.toString() << " at " << hour << " expected next time slot  to be NULL");
         BOOST_CHECK_MESSAGE(!timeSeries2.requeueable(c),"Time series " << timeSeries2.toString() << " should NOT be requeueable at hour " << hour );
      }
      else if (hour > timeSeries2.finish().hour()) {
         TimeSlot next_time_slot = timeSeries2.compute_next_time_slot(c);
         BOOST_CHECK_MESSAGE(next_time_slot.isNULL()," Time series " << timeSeries2.toString() << " at " << hour << " expected next time slot to be NULL");
         BOOST_CHECK_MESSAGE(!timeSeries2.requeueable(c),"Time series " << timeSeries2.toString() << " should NOT be requeueable at hour " << hour );
      }

      if (hour < timeSeries3.start().hour()) {
         TimeSlot next_time_slot = timeSeries3.compute_next_time_slot(c);
         TimeSlot expected(15,0);
         BOOST_CHECK_MESSAGE(next_time_slot ==expected," Time series " << timeSeries3.toString() << " at " << hour << " expected next time slot at " << expected.toString() << " but found " << next_time_slot.toString());
         BOOST_CHECK_MESSAGE(timeSeries3.requeueable(c),"Time series " << timeSeries3.toString() << " should be requeueable at hour " << hour );
      }
      else if (hour == timeSeries3.start().hour() ) {
         TimeSlot next_time_slot = timeSeries3.compute_next_time_slot(c);
         BOOST_CHECK_MESSAGE(next_time_slot.isNULL()," Time series " << timeSeries3.toString() << " at " << hour << " expected next time slot to be NULL");
         BOOST_CHECK_MESSAGE(!timeSeries3.requeueable(c),"Time series " << timeSeries3.toString() << " should NOT be requeueable at hour " << hour );
      }
      else if (hour > timeSeries3.start().hour()) {
         TimeSlot next_time_slot = timeSeries3.compute_next_time_slot(c);
         BOOST_CHECK_MESSAGE(next_time_slot.isNULL()," Time series " << timeSeries3.toString() << " at " << hour << " expected next time slot to be NULL");
         BOOST_CHECK_MESSAGE(!timeSeries3.requeueable(c),"Time series " << timeSeries3.toString() << " should NOT be requeueable at hour " << hour );
      }
   }
}

BOOST_AUTO_TEST_CASE( test_time_series_finish_not_divisble_by_increment )
{
	cout << "ACore:: ...test_time_series_finish_not_divisble_by_increment\n";

	// HANDLE CASE WHERE FINISH MINUTES IS NOT DIVISIBLE BY THE INCREMENT

	// Test time series with  a calendar, we update calendar then
	// test time series isFree(), and checkForRequeue
	Calendar calendar;
	calendar.init(ptime(date(2008,10,8), hours(0) ), Calendar::REAL);

	// Create a test when we can match a time series.
	// NOTE: Finish minute is not a multiple of INCREMENT hence last valid time slot is 23:50
	//
	// Create the time series: start  00:00
	//                         finish 23:59
	//                         incr   00:10
	TimeSeries timeSeries(TimeSlot(0,0), TimeSlot(23,59), TimeSlot(0,10), false/* relative */);
	TimeSeries timeSeries2(TimeSlot(0,30), TimeSlot(23,59), TimeSlot(4,0), false/* relative */);

   TimeSlot t1_min, t1_max,t2_min,t2_max;
   timeSeries.min_max_time_slots(t1_min, t1_max);
   timeSeries2.min_max_time_slots(t2_min, t2_max);
   BOOST_CHECK_MESSAGE(t1_min == TimeSlot(0,0) && t1_max == TimeSlot(23,59),"Not as expected");
   BOOST_CHECK_MESSAGE(t2_min == TimeSlot(0,30) && t2_max == TimeSlot(23,59),"Not as expected");

 	time_duration last = hours(23) + minutes(50);  // last valid time is 23:50
	time_duration last2 = hours(20) + minutes(30);  // last valid time is 20:30

	// follow normal process
   timeSeries.reset(calendar);
   timeSeries2.reset(calendar);

 	for(int hour=0; hour < 24; hour++) {
 		for( int minute=0; minute<60; minute++) {

 			calendar.update( minutes(1) );
         timeSeries.calendarChanged(calendar);
         timeSeries2.calendarChanged(calendar);

 			//cout << to_simple_string(calendar.suiteTime()) << "\n";

 			if (calendar.dayChanged()) {
            BOOST_CHECK_MESSAGE(timeSeries.checkForRequeue(calendar,t1_min,t1_max)," expected " << timeSeries.toString() << " checkForRequeue to pass at " << to_simple_string(calendar.suiteTime()));
 			}
 			else if ( calendar.suiteTime().time_of_day() < timeSeries.start().duration()) {
            BOOST_CHECK_MESSAGE(timeSeries.checkForRequeue(calendar,t1_min,t1_max)," expected " << timeSeries.toString() << " checkForRequeue to pass at " << to_simple_string(calendar.suiteTime()) );
 			}
 			else if ( calendar.suiteTime().time_of_day() >= timeSeries.start().duration() && calendar.suiteTime().time_of_day() < last	) {
     			BOOST_CHECK_MESSAGE(timeSeries.checkForRequeue(calendar,t1_min,t1_max)," expected " << timeSeries.toString() << " checkForRequeue to pass at " << to_simple_string(calendar.suiteTime())  );
  			}
  			else {
     			BOOST_CHECK_MESSAGE(!timeSeries.checkForRequeue(calendar,t1_min,t1_max)," expected " << timeSeries.toString() << " checkForRequeue to fail at " << to_simple_string(calendar.suiteTime()) );
  			}

         if (calendar.dayChanged()) {
            BOOST_CHECK_MESSAGE(!timeSeries2.checkForRequeue(calendar,t2_min,t2_max)," expected " << timeSeries2.toString() << " checkForRequeue to pass at " << to_simple_string(calendar.suiteTime()));
         }
         else if ( calendar.suiteTime().time_of_day() < timeSeries2.start().duration()) {
            BOOST_CHECK_MESSAGE(!timeSeries2.checkForRequeue(calendar,t2_min,t2_max)," expected " << timeSeries2.toString() << " checkForRequeue to pass at " << to_simple_string(calendar.suiteTime()) );
         }
         else if ( calendar.suiteTime().time_of_day() >= timeSeries2.start().duration() && calendar.suiteTime().time_of_day() < last2	) {
     			BOOST_CHECK_MESSAGE(timeSeries2.checkForRequeue(calendar,t2_min,t2_max)," expected " << timeSeries2.toString() << " checkForRequeue to pass at " << to_simple_string(calendar.suiteTime()));
  			}
  			else {
     			BOOST_CHECK_MESSAGE(!timeSeries2.checkForRequeue(calendar,t2_min,t2_max)," expected " << timeSeries2.toString() << " checkForRequeue to fail at " << to_simple_string(calendar.suiteTime()));
  			}

//         if (calendar.dayChanged()) {
//            timeSeries.reset(calendar);
//            timeSeries2.reset(calendar);
//         }
 		}
 	}
}

BOOST_AUTO_TEST_CASE( test_time_series_miss_time_slot )
{
   cout << "ACore:: ...test_time_series_miss_time_slot\n";

   Calendar calendar;
   calendar.init(ptime(date(2008,10,8), hours(10) ), Calendar::REAL);

   // Create a test when we can match a time series.
   // NOTE: Finish minute is not a multiple of INCREMENT hence last valid time slot is 23:50
   //
   // Create the time series: start  10:00
   //                         finish 23:59
   //                         incr   00:10
   TimeSeries timeSeries(TimeSlot(10,0), TimeSlot(23,59), TimeSlot(0,10), false/* relative */);
   timeSeries.miss_next_time_slot(); // will increment next_time_slot
   BOOST_CHECK_MESSAGE(timeSeries.value() == TimeSlot(10,10), "miss time slot not working ");

   timeSeries.miss_next_time_slot(); // will increment next_time_slot
   BOOST_CHECK_MESSAGE(timeSeries.value() == TimeSlot(10,20), "miss time slot not working ");

   // requeue time series,
   timeSeries.requeue(calendar); // calendar time is at midnight 10:00
   BOOST_CHECK_MESSAGE(timeSeries.value() == TimeSlot(10,10),"Expected requeue to reset time series");
}

BOOST_AUTO_TEST_CASE( test_time_series_miss_time_slot_1 )
{
   cout << "ACore:: ...test_time_series_miss_time_slot_1\n";

   // Create calendar before time series
   {
      Calendar calendar;
      calendar.init(ptime(date(2008,10,8), hours(9) ), Calendar::REAL);  // 09:00

      TimeSeries timeSeries(10,0, false/* relative */);                 // 10:00

      TimeSlot t1_min, t1_max;
      timeSeries.min_max_time_slots(t1_min, t1_max);
      BOOST_CHECK_MESSAGE(t1_min == TimeSlot(10,0) && t1_max == TimeSlot(10,0),"Not as expected");

      BOOST_CHECK_MESSAGE(!timeSeries.isFree(calendar), "Expected time holding at 10:00 since calendar is 09:00");
      BOOST_CHECK_MESSAGE(!timeSeries.checkForRequeue(calendar,t1_min,t1_max), "Expected checkForRequeue to return false always");

      timeSeries.miss_next_time_slot();
      BOOST_CHECK_MESSAGE(!timeSeries.isFree(calendar), "Expected time to hold");
      BOOST_CHECK_MESSAGE(!timeSeries.checkForRequeue(calendar,t1_min,t1_max), "Expected checkForRequeue to return false");

      timeSeries.requeue(calendar);
      BOOST_CHECK_MESSAGE(!timeSeries.isFree(calendar), "Expected time holding at 10:00 since calendar is 09:00");
      BOOST_CHECK_MESSAGE(!timeSeries.checkForRequeue(calendar,t1_min,t1_max), "Expected checkForRequeue to return false after Requeue");
   }
   {
      // Create a time after the time series
      Calendar calendar;
      calendar.init(ptime(date(2008,10,8), hours(11) ), Calendar::REAL); // 11:00

      TimeSeries timeSeries(10,0, false/* relative */);                 // 10:00

      TimeSlot t1_min, t1_max;
      timeSeries.min_max_time_slots(t1_min, t1_max);
      BOOST_CHECK_MESSAGE(t1_min == TimeSlot(10,0) && t1_max == TimeSlot(10,0),"Not as expected");

      BOOST_CHECK_MESSAGE(!timeSeries.isFree(calendar), "Expected time NOT to be free, since calendar time is after time slot");
      BOOST_CHECK_MESSAGE(!timeSeries.checkForRequeue(calendar,t1_min,t1_max), "Expected checkForRequeue to return false, since calendar time > slot time");

      timeSeries.miss_next_time_slot();
      BOOST_CHECK_MESSAGE(!timeSeries.isFree(calendar), "Expected time to hold");
      BOOST_CHECK_MESSAGE(!timeSeries.checkForRequeue(calendar,t1_min,t1_max), "Expected checkForRequeue to return false");

      timeSeries.requeue(calendar);
      BOOST_CHECK_MESSAGE(!timeSeries.isFree(calendar), "Expected time to hold, since calendar time is after time slot");
      BOOST_CHECK_MESSAGE(!timeSeries.checkForRequeue(calendar,t1_min,t1_max), "Expected checkForRequeue to return false,since calendar time is after time slot");
   }
}

BOOST_AUTO_TEST_CASE( test_time_series_reset )
{
   cout << "ACore:: ...test_time_series_reset\n";

   // Create a test when we can match a time series.
   // NOTE: Finish minute is not a multiple of INCREMENT hence last valid time slot is 23:50
   //
   // Create the time series: start  10:00
   //                         finish 23:59
   //                         incr   00:10
   TimeSeries timeSeries(TimeSlot(10,0), TimeSlot(23,59), TimeSlot(0,10), false/* relative */);
   BOOST_CHECK_MESSAGE(timeSeries.value() == TimeSlot(10,0), "defaults not correct");

   { // set calendar before time series start & then reset

      Calendar calendar;
      calendar.init(ptime(date(2008,10,8), hours(9) ), Calendar::REAL);
      timeSeries.reset(calendar);
      BOOST_CHECK_MESSAGE(timeSeries.value() == TimeSlot(10,0), "Reset should set value(next_valid_time_slot) to start." << timeSeries.dump());

      timeSeries.miss_next_time_slot();

      timeSeries.requeue(calendar);
      BOOST_CHECK_MESSAGE(timeSeries.value() == TimeSlot(10,0), "requeue should set value(next_valid_time_slot) to start."<< timeSeries.dump());
   }
   { // set calendar at time series start & then reset

      Calendar calendar;
      calendar.init(ptime(date(2008,10,8), hours(10) ), Calendar::REAL);
      timeSeries.reset(calendar);
      BOOST_CHECK_MESSAGE(timeSeries.value() == TimeSlot(10,0), "Reset should set value(next_valid_time_slot) to start."<< timeSeries.dump());

      timeSeries.miss_next_time_slot();

      timeSeries.requeue(calendar);
      BOOST_CHECK_MESSAGE(timeSeries.value() == TimeSlot(10,10), "requeue should set value(next_valid_time_slot) to start."<< timeSeries.dump());
   }
   { // set calendar after time series start & before time series end, then reset

      Calendar calendar;
      calendar.init(ptime(date(2008,10,8), hours(11) ), Calendar::REAL);
      timeSeries.reset(calendar);
      BOOST_CHECK_MESSAGE(timeSeries.value() == TimeSlot(11,0), "Reset should update free slot."<< timeSeries.dump());

      timeSeries.miss_next_time_slot();

      timeSeries.requeue(calendar);
      BOOST_CHECK_MESSAGE(timeSeries.value() == TimeSlot(11,10), "requeue should update to first time slot."<< timeSeries.dump());
   }
}

BOOST_AUTO_TEST_CASE( test_time_series_parsing )
{
   cout << "ACore:: ...test_time_series_parsing\n";

   BOOST_CHECK_MESSAGE(TimeSeries::create("00:30") == ecf::TimeSeries(0,30), "Error ");
   BOOST_CHECK_MESSAGE(TimeSeries::create("+00:30") == ecf::TimeSeries(0,30,true), "Error ");

   ecf::TimeSlot start(0,30);
   ecf::TimeSlot finish(21,03);
   ecf::TimeSlot incr(1,30);
   BOOST_CHECK_MESSAGE(TimeSeries::create("00:30 21:03 01:30") == ecf::TimeSeries(start,finish,incr), "Error");
   BOOST_CHECK_MESSAGE(TimeSeries::create("+00:30 21:03 01:30") == ecf::TimeSeries(start,finish,incr,true), "Error");
}

BOOST_AUTO_TEST_CASE( test_time_series_state_parsing )
{
   cout << "ACore:: ...test_time_series_state_parsing\n";

   /// extract string like
   ///     time +00:00 20:00 00:10 # this is a comment which will be ignored. index = 1
   ///     time +20:00        // index = 1
   ///     today 20:00        // index = 1
   ///     +00:00 20:00 00:10 // index = 0
   ///     +20:00             // index = 0
   /// will throw std:runtime_error for errors
   /// will assert if index >= lineTokens.size()
   {
      size_t index = 0;
      std::string the_time = "+00:00 20:00 00:10";
      std::vector<std::string> lineTokens;
      Str::split(the_time,lineTokens);
      BOOST_CHECK_MESSAGE(TimeSeries::create(index,lineTokens) == ecf::TimeSeries(TimeSlot(0,0),TimeSlot(20,0),TimeSlot(0,10),true), "Error");
   }
   {
      size_t index = 0;
      std::string the_time = "+10:10";
      std::vector<std::string> lineTokens;
      Str::split(the_time,lineTokens);
      BOOST_CHECK_MESSAGE(TimeSeries::create(index,lineTokens) == ecf::TimeSeries(TimeSlot(10,10),true), "Error");
   }
   {
      size_t index = 0;
      std::string the_time = "+10:10 # free isValid:false";
      ecf::TimeSeries expected(TimeSlot(10,10),true);
      expected.set_isValid(false);

      std::vector<std::string> lineTokens;
      Str::split(the_time,lineTokens);
      ecf::TimeSeries parsed_ts = TimeSeries::create(index,lineTokens,true);
      BOOST_CHECK_MESSAGE(parsed_ts == expected,
                          "Expected \n'" << expected.toString() << expected.state_to_string(false) << "'" <<
                          " But found \n'" << parsed_ts.toString() << parsed_ts.state_to_string(false) << "'");
   }
   {
      size_t index = 0;
      std::string the_time = "+10:10 # free isValid:false nextTimeSlot/10:10 relativeDuration/00:00:00";
      ecf::TimeSeries expected(TimeSlot(10,10),true);
      expected.set_isValid(false);

      std::vector<std::string> lineTokens;
      Str::split(the_time,lineTokens);
      ecf::TimeSeries parsed_ts = TimeSeries::create(index,lineTokens,true);
      BOOST_CHECK_MESSAGE(parsed_ts == expected,
                          "Expected \n'" << expected.toString() << expected.state_to_string(false) << "'" <<
                          " But found\n'" << parsed_ts.toString() << parsed_ts.state_to_string(false) << "'");
   }
   {
      // Update relative duration
      size_t index = 0;
      std::string the_time = "+10:10 # free isValid:false nextTimeSlot/10:10 relativeDuration/01:00:00";
      ecf::TimeSeries expected(TimeSlot(10,10),true);
      expected.set_isValid(false);

      Calendar calendar;
      calendar.init(ptime(date(2008,10,8), hours(0) ), Calendar::REAL);
      calendar.update( time_duration( hours(1) ) );
      expected.calendarChanged( calendar );

      std::vector<std::string> lineTokens;
      Str::split(the_time,lineTokens);
      ecf::TimeSeries parsed_ts = TimeSeries::create(index,lineTokens,true);
      BOOST_CHECK_MESSAGE(parsed_ts == expected,
                          "Expected \n'" << expected.toString() << expected.state_to_string(false) << "'" <<
                          " But found\n'" << parsed_ts.toString() << parsed_ts.state_to_string(false) << "'");
   }

   {
      // Update nextTimeSlot, create calendar at 09:00 and increment time series.
      std::string the_time = "09:00 12:00 01:00 # isValid:false nextTimeSlot/10:00 relativeDuration/00:00:00";
      ecf::TimeSeries expected(TimeSlot(9,0),TimeSlot(12,0),TimeSlot(1,0),false);
      expected.set_isValid(false);

      Calendar calendar;
      calendar.init(ptime(date(2008,10,8), hours(9) ), Calendar::REAL);
      expected.requeue( calendar ); // this will reset isValid to true
      expected.set_isValid(false);

      std::vector<std::string> lineTokens;
      Str::split(the_time,lineTokens);
      size_t index = 0;
      ecf::TimeSeries parsed_ts = TimeSeries::create(index,lineTokens,true);

      BOOST_CHECK_MESSAGE(parsed_ts == expected,
                          "Expected \n'" << expected.dump() << expected.state_to_string(false) << "'" <<
                          " But found\n'" << parsed_ts.dump() << parsed_ts.state_to_string(false) << "'");
   }
}

BOOST_AUTO_TEST_SUITE_END()

max_stars_repo_name: ecmwf/ecflow

max_issues_repo_name: CoollRock/ecflow

max_forks_repo_name: CoollRock/ecflow


text:
#include <cassert>
#include <cmath>
#include <string>
#include <sstream>
#include <iostream>
#include <fstream>
#include <algorithm>
#include <random>
#include <ctime>              // required for "struct timezone" and
                              // "gettimeofday()" used to set the randomseed
using namespace std;

#ifndef DISABLE_OPENMP
#include <omp.h>       // (OpenMP-specific)
#endif

//#include <boost/math/special_functions/gamma.hpp>

#include <mrc_simple.hpp>
#include <visfd.hpp>
using namespace visfd;

#include "err.hpp"
#include "settings.hpp"

int main(int argc, char **argv) {
  try {

    Settings settings; // parse the command-line argument list from the shell
    settings.ParseArgs(argc, argv);

    #ifndef DISABLE_OPENMP
    #pragma omp parallel
    {
      int rank, nthr;
      rank = omp_get_thread_num();
      //cerr << "rank=" << rank << endl;
      if (rank == 0) {
        nthr = omp_get_num_threads();
        cerr << "   (Using " << nthr << " threads (cpu cores).  You can change this using the \"-np n\"\n"
             << "    argument, or by setting the OMP_NUM_THREADS environment variable.)" << endl;
      }
    }
    #else
    cerr << " (Serial version)" << endl;
    #endif //#ifndef DISABLE_OPENMP

    // Read the input tomogram
    MrcSimple tomo_in;

    if (settings.in_file_name != "") {
      cerr << "Reading tomogram \""<<settings.in_file_name<<"\"" << endl;
      tomo_in.Read(settings.in_file_name, false);
      // (Note: You can also use "tomo_in.Read(cin);" or "cin >> tomo;")
      tomo_in.PrintStats(cerr);      //Optional (display the tomogram size & format)
      WarnMRCSignedBytes(tomo_in, settings.in_file_name, cerr);
    }
    else {
      assert((settings.image_size[0] > 0) &&
             (settings.image_size[1] > 0) &&
             (settings.image_size[2] > 0));
      tomo_in.Resize(settings.image_size);
    }

    // ---- mask ----

    // Optional: if there is a "mask", read that too
    MrcSimple mask;
    if (settings.mask_file_name != "") {
      cerr << "Reading mask \""<<settings.mask_file_name<<"\"" << endl;
      mask.Read(settings.mask_file_name, false);
      if ((mask.header.nvoxels[0] != tomo_in.header.nvoxels[0]) ||
          (mask.header.nvoxels[1] != tomo_in.header.nvoxels[1]) ||
          (mask.header.nvoxels[2] != tomo_in.header.nvoxels[2]))
        throw InputErr("Error: The size of the mask image does not match the size of the input image.\n");
    }

    // ---- Voxel size? ----

    float voxel_width[3] = {1.0, 1.0, 1.0};
    if (settings.voxel_width > 0.0) {
      // Did the user manually specify the width of each voxel?
      voxel_width[0] = settings.voxel_width;
      voxel_width[1] = settings.voxel_width;
      voxel_width[2] = settings.voxel_width;
    }
    else {
      // Otherwise, infer it from the header of the MRC file
      voxel_width[0] = tomo_in.header.cellA[0]/tomo_in.header.nvoxels[0];
      voxel_width[1] = tomo_in.header.cellA[1]/tomo_in.header.nvoxels[1];
      voxel_width[2] = tomo_in.header.cellA[2]/tomo_in.header.nvoxels[2];
      if (settings.voxel_width_divide_by_10) {
        voxel_width[0] *= 0.1;
        voxel_width[1] *= 0.1;
        voxel_width[2] *= 0.1;
      }
      cerr << "voxel width in physical units = ("
           << voxel_width[0] << ", "
           << voxel_width[1] << ", "
           << voxel_width[2] << ")\n";
    }

    if ((voxel_width[0] <= 0.0) ||
        (voxel_width[1] <= 0.0) ||
        (voxel_width[2] <= 0.0))
      throw InputErr("Error in tomogram header: Invalid voxel width(s).\n"
                     "Use the -w argument to specify the voxel width.");

    if ((voxel_width[0] != voxel_width[1]) ||
        (voxel_width[1] != voxel_width[2]))
      throw InputErr("Error in tomogram header: Unequal voxel widths in the x, y and directions.\n"
                     "Use the -w argument to specify the voxel width.");

    if (settings.in_coords_file_name != "") {
      for (int iz = 0; iz < tomo_in.header.nvoxels[2]; iz++)
        for (int iy = 0; iy < tomo_in.header.nvoxels[1]; iy++)
          for (int ix = 0; ix < tomo_in.header.nvoxels[0]; ix++)
            tomo_in.aaafI[iz][iy][ix] = 0.0;
      fstream coords_file;
      coords_file.open(settings.in_coords_file_name, ios::in);
      if (! coords_file)
        throw InputErr("Error: unable to open \""+
                       settings.in_coords_file_name +"\" for reading.\n");
      while (coords_file) {
        float x, y, z;
        coords_file >> x;
        coords_file >> y;
        coords_file >> z;
        int ix, iy, iz;
        ix = static_cast<int>(x / voxel_width[0]);
        iy = static_cast<int>(y / voxel_width[1]);
        iz = static_cast<int>(z / voxel_width[2]);
        if (((0 <= ix) && (ix < tomo_in.header.nvoxels[0])) &&
            ((0 <= iy) && (iy < tomo_in.header.nvoxels[1])) &&
            ((0 <= iz) && (iz < tomo_in.header.nvoxels[2])))
          tomo_in.aaafI[iz][iy][ix] = 1.0;
      }
    } //if (settings.in_coords_file_name != "") {

    double vol_total = settings.compartment_volume;
    if (settings.compartment_volume < 0) {
      if (mask.aaafI) {
        for (int iz=0; iz<tomo_in.header.nvoxels[2]; iz++)
          for (int iy=0; iy<tomo_in.header.nvoxels[1]; iy++)
            for (int ix=0; ix<tomo_in.header.nvoxels[0]; ix++)
              vol_total += mask.aaafI[iz][iy][ix];
      }
      else {
        vol_total = (tomo_in.header.nvoxels[0] *
                     tomo_in.header.nvoxels[1] *
                     tomo_in.header.nvoxels[2]);
      }
      vol_total *= (voxel_width[0]*voxel_width[1]*voxel_width[2]);
    }

    if (settings.num_particles < 0) {
      // If the user did not explicitly specify the number of particles
      // (ie "ribosomes") present in the original image, we can attempt
      // to infer it from the image itself.  To do this, I assume that
      // the image is filled with voxels containing 0s everywhere except
      // at the center of where each particle is located, where the voxel is 1.
      // Consequently, the number of non-zero voxels (not excluded by the mask)
      // should equal the sum of all of the voxel entries.
      // (This should also equal the number of non-zero voxels.  However
      //  sometimes these images are themselves slightly blurred to make it
      //  easier to see each particle.  This blurring should not effect the sum)
      double sum = 0.0;
      for (int iz=0; iz<tomo_in.header.nvoxels[2]; iz++) {
        for (int iy=0; iy<tomo_in.header.nvoxels[1]; iy++) {
          for (int ix=0; ix<tomo_in.header.nvoxels[0]; ix++) {
            if (mask.aaafI)
              sum += tomo_in.aaafI[iz][iy][ix] * mask.aaafI[iz][iy][ix];
            else
              sum += tomo_in.aaafI[iz][iy][ix];
          }
        }
      }
      settings.num_particles = sum;
    }

    // DEBUG: The next if-then statement is for debugging only.
    if (settings.randomize_input_image) {
      size_t nparticles = floor(settings.num_particles);
      size_t nvoxels;
      if (mask.aaafI) {
        for (int iz=0; iz<tomo_in.header.nvoxels[2]; iz++)
          for (int iy=0; iy<tomo_in.header.nvoxels[1]; iy++)
            for (int ix=0; ix<tomo_in.header.nvoxels[0]; ix++)
              if (mask.aaafI[iz][iy][ix] != 0.0)
                nvoxels += 1;
      }
      else {
        nvoxels = (tomo_in.header.nvoxels[0] *
                   tomo_in.header.nvoxels[1] *
                   tomo_in.header.nvoxels[2]);
      }
      vector<bool> random_bit_list(nvoxels, false);
      for (size_t i=0; i < nparticles; i++)
        random_bit_list[i] = true;

      long random_seed = settings.random_seed;
      if (random_seed <= 0) {
        random_seed = time(nullptr);
        cerr << "(random_seed = " << random_seed << ")" << endl;
      }
      shuffle(random_bit_list.begin(), random_bit_list.end(),
              default_random_engine(random_seed));
      size_t i=0;
      for (int iz = 0; iz < tomo_in.header.nvoxels[2]; iz++) {
        for (int iy = 0; iy < tomo_in.header.nvoxels[1]; iy++) {
          for (int ix = 0; ix < tomo_in.header.nvoxels[0]; ix++) {
            tomo_in.aaafI[iz][iy][ix] = 0.0;
            if ((! mask.aaafI) || mask.aaafI[iz][iy][ix] != 0.0) {
              tomo_in.aaafI[iz][iy][ix] = random_bit_list[i];
              i++;
            }
          }
        }
      }
    } //if (settings.randomize_input_image)

    // First we must allocate a new array to store the filtered image.
    MrcSimple tomo_out = tomo_in;

    for (int i_sig = 0; i_sig < settings.vfSigma.size(); ++i_sig) {

      // ---- main calculation: -----

      //A Gaussian filter is a (weighted) average over a volume of nearby voxels
      //Below we estimate the effective volume of those nearby voxels, by taking
      //advantage of the fact that Gaussian is normalized (it's integral = 1).
      //Normalization effectively reduces the height of the filter by the number
      //of voxels which effectively belong to the filter.
      // Example: a cube filter
      // If we were averaging the brightness of the 27 voxels
      // belonging to a 3*3*3 cube surrounding the voxel of interest
      // then the effective volume of the filter is is 3*3*3=27, and
      // weight of each voxel's contribution to the average is 1/27.  Similarly,
      // for a more general filter, we can interpret the peak height of the
      // filter as 1/volume. Below, we calculate the height of the 3D Gaussian's
      // peak after normalization.  Normally the 3D Gaussian is:
      //   (2*pi*sigma^2)^(-3/2) * exp(-0.5*(x^2+y^2+z^2)/sigma^2)
      // and the peak height is (2*pi*sigma^2)^(-3/2),
      // However this formula does not work for small sigma << 1.
      // For small sigma, you have to normalize by computing a discrete sum
      // (over a finite range).  We do this below (in a non-intuitive way).

      float sigma = settings.vfSigma[i_sig];

      sigma /= voxel_width[0]; //scale by voxel_width (if specified)

      // How wide is the filter window (over what size window do we integrate)?
      if (settings.filter_truncate_ratio <= 0) {
        assert(settings.filter_truncate_threshold > 0.0);
        settings.filter_truncate_ratio = sqrt(-2*log(settings.filter_truncate_threshold));
      }
      int filter_truncate_halfwidth = floor(sigma *
                                            settings.filter_truncate_ratio);

      float gauss_peak_height_3D;

      {
        // GenFilterGauss1D generates a 1D discrete normalized Gaussian (located
        // in the afH[]) member. afH[0] is the height at the central peak.
        // However, for now I just use them to estimate the height of the
        // Gaussian peak.  (See below)
        Filter1D<float, int> filter1D =
          GenFilterGauss1D(sigma, filter_truncate_halfwidth);

        // A 3D Gaussian blur can be performed by expoiting the fact that
        // multidimensional Gaussians are "seperable" filters:  You can blur
        // in the X direction, then in the Y direcion, then in the Z direction.
        // The result is the same that you would get from a 3D filter Gaussian
        // which is the product of Gaussians in the X,Y,Z directions.
        // So the peak height is the product of the 3 1D Gaussian peak heights.
        gauss_peak_height_3D =
          (filter1D.afH[0] *
           filter1D.afH[0] *
           filter1D.afH[0]);
      }

      // The Gaussian peak's height is 1/volume.

      float volume_gaussian_bin = 1.0 / gauss_peak_height_3D;
      volume_gaussian_bin *= (voxel_width[0]*voxel_width[1]*voxel_width[2]);

      float num_bins = vol_total / volume_gaussian_bin;

      // Now blur the original image (if the user did not do it already)
      if ((! settings.precomputed_gaussian_blur) ||
          (settings.vfSigma.size() > 1)) {

        ApplyGauss(tomo_in.header.nvoxels,
                   tomo_in.aaafI,
                   tomo_out.aaafI, //<-store resulting image here
                   mask.aaafI,
                   sigma,
                   filter_truncate_halfwidth,
                   true,
                   &cerr);

        // Densities everywhere are assumed to be in physical units
        // (ie. of 1/Angstroms^3),  NOT   1/voxels^3
        // The Gaussian blur computes densities in 1/voxels^3
        // To compensate for this, divide the densities by voxel_width^3
        MultiplyScalarArr(static_cast<float>(
                          1.0/(voxel_width[0]*voxel_width[1]*voxel_width[2])),
                          tomo_out.header.nvoxels,
                          tomo_out.aaafI);
      }

      // After bluring the image, find the lowest density:

      float extreme_density;
      int afXextreme[3] = {-1, -1, -1};
      float global_minima;
      float global_maxima;

      // Find the voxels with the minima and maxima intensities.
      // Discard global minima or maxima which lie on the boundary
      // The safe, careful way to do this is to only consider
      // voxels which are surrounded by other voxels and are
      // local minima or maxima.
      for (int iz=0; iz < tomo_out.header.nvoxels[2]; iz++) {
        for (int iy=0; iy < tomo_out.header.nvoxels[1]; iy++) {
          for (int ix=0; ix < tomo_out.header.nvoxels[0]; ix++) {

            bool is_local_minima = true;
            bool is_local_maxima = true;

            if (! settings.extrema_on_boundary) {
              if ((ix == 0) || (ix ==  tomo_out.header.nvoxels[0]-1) ||
                  (ix == 0) || (ix ==  tomo_out.header.nvoxels[0]-1) ||
                  (ix == 0) || (ix ==  tomo_out.header.nvoxels[0]-1)) {
                is_local_minima = false;
                is_local_maxima = false;
              }
            }

            float center_val = tomo_out.aaafI[iz][iy][ix];
            for (int jz=-1; jz<=1; jz++) {
              for (int jy=-1; jy<=1; jy++) {
                for (int jx=-1; jx<=1; jx++) {
                  if (! settings.extrema_on_boundary) {
                    if (mask.aaafI && (mask.aaafI[iz+jz][iy+jy][ix+jx] == 0.0))
                    {
                      is_local_minima = false;
                      is_local_maxima = false;
                    }
                  }
                  if (tomo_out.aaafI[iz+jz][iy+jy][ix+jx] <= center_val)
                    is_local_minima = false;
                  if (tomo_out.aaafI[iz+jz][iy+jy][ix+jx] >= center_val)
                    is_local_maxima = false;
                }
              }
            }
            if ((center_val < global_minima) || (afXextreme[0] == -1)) {
              global_minima = center_val;
              if (settings.use_min_density) {
                extreme_density = global_minima;
                afXextreme[0] = ix;
                afXextreme[1] = iy;
                afXextreme[2] = iz;
              }
            }
            if ((center_val > global_maxima) || (afXextreme[0] == -1)) {
              global_maxima = center_val;
              if (! settings.use_min_density) {
                extreme_density = global_maxima;
                afXextreme[0] = ix;
                afXextreme[1] = iy;
                afXextreme[2] = iz;
              }
            }
          } //for (int ix=1; ix < tomo_out.header.nvoxels[0]-1; ix++)
        } //for (int iy=1; iy < tomo_out.header.nvoxels[1]-1; iy++)
      } //for (int iz=1; iz < tomo_out.header.nvoxels[2]-1; iz++)

      // Did we fail to find any local minima or maxima densities?
      if (afXextreme[0] == -1) {
        string extrema_type = "minina";
        if (! settings.use_min_density)
          extrema_type = "maxima";
        stringstream msg_ss;
        msg_ss << "Error: There are no local density " << extrema_type << "\n"
               << "       in the image (at scale sigma = " << sigma << ")\n"
               << "       Aborting...\n";
        throw InputErr(msg_ss.str());
      }

      float ave_density = settings.num_particles / vol_total;

      // Suppose you divide the volume of the region (ie. cell, compartment...)
      // into equal size "bins" of equal size ("volume_gaussian_bin").
      // A fixed number of particles lie within this volume.
      // ...What is the probability that the number of particles in a bin
      //    does not exceed "extreme_density"?
      //  This is given by the (cumulative) Poisson distribution:
      //     prob = lambda^k * exp(-lambda) / k!
      //          = lambda^k * exp(-lambda) / Gamma(k+1)  <-continuous version
      // "k" is the number of particles in this Gaussian-shaped "bin"
      long double k = extreme_density * volume_gaussian_bin;
      // "lambda" is the expected number of particles in this "bin"
      long double lambda = ave_density * volume_gaussian_bin;

      long double prob_cdf_obvious_way = 0.0;
      if (settings.use_min_density) {
        // Calculate the culmulative probability of seeing this many
        // particles in the bin or less.  First try it the obvious way:

        for (long i=0; i <= floor(k); i++)
          prob_cdf_obvious_way += pow(lambda, i) * exp(-lambda) / tgamma(i+1.0);
      }
      else {
      // Calculate the culmulative probability of seeing this many
      // particles in the bin or more.  First try it the obvious way:
        for (long i=0; i < floor(k); i++)
          prob_cdf_obvious_way += pow(lambda, i) * exp(-lambda) / tgamma(i+1.0);
        prob_cdf_obvious_way = 1.0 - prob_cdf_obvious_way;
      }

      cerr << "####################################" << endl;
      cerr << "## DEBUG MESSAGES (PLEASE IGNORE) ##" << endl;
      cerr << "## location (in voxels) = ("
           << afXextreme[0] << ","
           << afXextreme[1] << ","
           << afXextreme[2] << ")" << endl;
      cerr << "## num_in_bin = " << k << endl;
      cerr << "## num_expected_in_bin = " << lambda << endl;
      cerr << "## prob_cdf_obvious_way = " << prob_cdf_obvious_way << endl;
      cerr << "####################################" << endl;

      long double prob_cdf;

      // However, in our case the number of particles in this "bin" is not (k)
      // necessarily an integer.  In that case use the continuous version of
      // (the cumulative distribution) of the Poisson distributio.
      // COMMENTING OUT: This requires the BOOST libraries
      // COMMENTING OUT: if (use_min_density)
      // COMMENTING OUT:   prob_cdf = gamma_q(k+1,lambda);
      // COMMENTING OUT: else
      // COMMENTING OUT:   prob_cdf = 1.0 - gamma_q(k,lambda);
      //    (gamma_q() is defined in boost)

      // The cumulative distribution of the Poisson distribution
      // is also given by the upper incomplete Gamma function
      // COMMENTING OUT: USE THE BOOST LIBRARY
      //double prob_cdf = gamma_q(floor(k)+1, lambda);
      // Unfortunately the  BOOST library can not be easily distributed with
      // this code.  (Even if I throw away most ofthe BOOST code using the
      // "bcp" utility, the remaining BOOST code exceeds the size of this
      // the code for this project by a factor of 50.)
      // Hence we will use the integer approximation ("obvious_way"):

      prob_cdf = prob_cdf_obvious_way;

      // ...Now what is the probability that the number of particles
      //    in ANY OF THE BINS does not exceed "extreme_density"?

      long double prob_total = 1.0 - pow((1.0 - prob_cdf), num_bins);

      double effective_bin_size = (pow(volume_gaussian_bin, 1.0/3) *
                                   voxel_width[0]);
      cout << prob_total
           << " " << extreme_density
           << " " << afXextreme[0]
           << " " << afXextreme[1]
           << " " << afXextreme[2]
           << " " << effective_bin_size
           << endl;

      // Discussion:
      //      What do we expect this number to be?
      // What is the expected value of "prob_total" for a randomly distributed
      // particles in the same number of bins.
      //
      // Answer: 0.5
      //
      // The proof is not specific to the kind of distribution we are using.
      // (a Poisson distribution)
      //
      // Proof:
      //
      // Let "c" be the cumulative probability distrubution for the entire
      // system (in this "prob_total", but it could also be "prob_cdf")
      //
      // Let p(x) be the probability density of measuring x.
      //
      // Let C(x) be the cumulative distribution of measuring something <= x.
      // (in this case the probability that none of the bins have
      //  a density exceeding x)  This means that:
      //
      // C(X) = \int_{-\infty}^X  dx  p(x)
      //
      // Let C^{-1}(c)  denote the inverse of C(x)
      // ==> C^{-1}(c) = x   if   C(x) = c
      //
      // What is the average value of c?
      //
      // <c> = \int_{-infty}^{\infty} * c(x) * p(x) * dx
      //
      //     = \int_0^1         c *  p(x(c)) * dx(c)/dc * dc
      //
      //     = \int_0^1 c * p( C^{-1}(c) ) * (d/dc) C^{-1}(c)  * dc
      //
      //     = \int_0^1 c * p( C^{-1}(c) ) * ( 1 / (d/dx) C( C^{-1}(c) ) ) * dc
      //
      //     = \int_0^1 c * p( C^{-1}(c) ) * ( 1 / p( C^{-1}(c) ) )  * dc
      //
      //     = \int_0^1 c * dc
      //
      //     = 0.5
      //
      // However this number will be lower if you repeat this multiple times
      // (using different bin-sizes ("sigma" values), for example)
      // because by repeating, you give it more opportunities to find
      // a more extreme value.  (Here we scan over a range of "sigma" values.)

    } // for (int i_sig = 0; i_sig < settings.vfSigma.size(); i_sig++) {...

    if ((settings.out_file_name != "") && (settings.vfSigma.size() == 1))
    {
      // Write the image file containing the filtered version for that sigma
      cerr << "writing a tomogram containing the most recently calculated density cloud"
           << endl;

      tomo_out.Write(settings.out_file_name);
      // (You can also use "file_stream << tomo_out;")
    }

  } //try {
  catch (const std::exception& e) {
    cerr << "\n" << e.what() << endl;
    exit(1);
  }

} // main()

max_stars_repo_name: jewettaij/visfd

max_issues_repo_name: jewettaij/visfd

max_forks_repo_name: jewettaij/visfd


text:
#include <iostream>
#include <Eigen/Dense>
using namespace Eigen;
using namespace std;

MatrixXcd bl_bicggr(const MatrixXcd& A, const MatrixXcd& B, const double& tol, const int& itermax)
{
// Block BiCGGR [Tadano etal 2009 JSIAM letters]
  double Bnorm= B.norm();
  MatrixXcd X= MatrixXcd::Zero(B.rows(),B.cols()); // Initial guess of X (zeros)
  MatrixXcd R= B-A*X;
  MatrixXcd P= R;
  MatrixXcd V= A*R;
  MatrixXcd W= V;
  MatrixXcd R0til= R; //MatrixXcd::Random(n,L);
  MatrixXcd R0til_H= R0til.adjoint();
  for(int k= 0; k < itermax; ++k){
    MatrixXcd alfa= (R0til_H*V).fullPivLu().solve(R0til_H*R);
    complex<double> qsi= (W.adjoint()*R).trace()/(W.adjoint()*W).trace();
    MatrixXcd S= P-qsi*V;
    MatrixXcd U= S*alfa;
    MatrixXcd Y= A*U;
    X= X+qsi*R+U;
    MatrixXcd Rnew= R-qsi*W-Y;
    double err= Rnew.norm()/Bnorm;
    cout << "bl_bicggr: " << "iter= " << k << " relative err= " << err << endl;
    if(err < tol) break;
    W= A*Rnew;
    MatrixXcd gamma = (R0til_H*R).fullPivLu().solve(R0til_H*Rnew/qsi);
    R=Rnew;
    P= R+U*gamma;
    V= W+Y*gamma;

  }
  if((A*X-B).norm()/Bnorm > 10*tol){
      cerr << "bl_bicggr did not converge to solution within error tolerance !" << endl;
     // exit(EXIT_FAILURE);
  }

  return X;
}

max_stars_repo_name: nmoteki/block-Krylov-linear-solvers

max_issues_repo_name: nmoteki/block-Krylov-linear-solvers

max_forks_repo_name: nmoteki/block-Krylov-linear-solvers