codac2_MatrixBase_addons_IntervalMatrixBase.h

Go to the documentation of this file.

template<typename U=Scalar>
  requires IsIntervalDomain<U>
inline double volume() const
{
  if(this->is_empty())
    return 0.;
 
  double v = 0.;
  for(Index i = 0 ; i < this->rows() ; i++)
    for(Index j = 0 ; j < this->cols() ; j++)
    {
      if((*this)(i,j).is_unbounded()) return codac2::oo;
      if((*this)(i,j).is_degenerated()) return 0.;
      v += std::log((*this)(i,j).diam());
    }
  return std::exp(v);
}

inline bool is_empty() const
{
  for(Index i = 0 ; i < rows() ; i++)
    for(Index j = 0 ; j < cols() ; j++)
      if((*this)(i,j).is_empty())
        return true;
  return false;
}

#define degenerate_mat(op) \
  Matrix<double,RowsAtCompileTime,ColsAtCompileTime> m(this->rows(),this->cols()); \
  \
  if(this->is_empty()) \
    m.init(std::numeric_limits<double>::quiet_NaN()); \
  \
  else \
  { \
    for(Index i = 0 ; i < this->rows() ; i++) \
      for(Index j = 0 ; j < this->cols() ; j++) \
        m(i,j) = (*this)(i,j).op(); \
  } \
  \
  return m; \


template<typename U=Scalar>
  requires IsIntervalDomain<U>
inline auto lb() const
{
  degenerate_mat(lb);
}

template<typename U=Scalar>
  requires IsIntervalDomain<U>
inline auto ub() const
{
  degenerate_mat(ub);
}

template<typename U=Scalar>
  requires IsIntervalDomain<U>
inline auto mid() const
{
  degenerate_mat(mid);
}

template<typename U=Scalar>
  requires IsIntervalDomain<U>
inline auto mag() const
{
  degenerate_mat(mag);
}

template<typename U=Scalar>
  requires IsIntervalDomain<U>
inline auto mig() const
{
  degenerate_mat(mig);
}

template<typename U=Scalar>
  requires IsIntervalDomain<U>
inline auto smag() const
{
  degenerate_mat(smag);
}

template<typename U=Scalar>
  requires IsIntervalDomain<U>
inline auto smig() const
{
  degenerate_mat(smig);
}

template<typename U=Scalar>
  requires IsIntervalDomain<U>
inline auto rand() const
{
  degenerate_mat(rand);
}

template<typename U=Scalar>
  requires IsIntervalDomain<U>
inline auto rad() const
{
  degenerate_mat(rad);
}

template<typename U=Scalar>
  requires IsIntervalDomain<U>
inline auto diam() const
{
  degenerate_mat(diam);
}

inline double min_rad(const std::vector<Index>& among_indices = {}) const
  requires IsIntervalDomain<Scalar>
{
  return coeff(this->extr_diam_index(true,among_indices)).rad();
}

inline double max_rad(const std::vector<Index>& among_indices = {}) const
  requires IsIntervalDomain<Scalar>
{
  return coeff(this->extr_diam_index(false,among_indices)).rad();
}

inline double min_diam(const std::vector<Index>& among_indices = {}) const
  requires IsIntervalDomain<Scalar>
{
  return coeff(this->extr_diam_index(true,among_indices)).diam();
}

inline double max_diam(const std::vector<Index>& among_indices = {}) const
  requires IsIntervalDomain<Scalar>
{
  return coeff(this->extr_diam_index(false,among_indices)).diam();
}

inline Index min_diam_index(const std::vector<Index>& among_indices = {}) const
  requires IsIntervalDomain<Scalar>
{
  return this->extr_diam_index(true,among_indices);
}

inline Index max_diam_index(const std::vector<Index>& among_indices = {}) const
  requires IsIntervalDomain<Scalar>
{
  return this->extr_diam_index(false,among_indices);
}

inline Index extr_diam_index(bool min, const std::vector<Index>& among_indices = {}) const
  requires IsIntervalDomain<Scalar>
{
  assert_release(
    among_indices.empty() ||
    std::ranges::all_of(among_indices, [&](Index k){ return 0 <= k && k < this->size(); })
  );
 
  // This code originates from the ibex-lib
  // See: ibex_TemplateVector.h
  // Author: Gilles Chabert
 
  double d = min ? codac2::oo : -1; // -1 to be sure that even a 0-diameter interval can be selected
  int selected_index = -1;
  bool unbounded = false;
  assert_release(!this->is_empty() && "Diameter of an empty IntervalVector is undefined");
 
  Index i;
 
  for(i = 0 ; i < this->size() ; i++) 
  {
    if(!among_indices.empty()
      && std::find(among_indices.begin(),among_indices.end(),i) == among_indices.end())
      continue; // avoiding this index
 
    if(coeff(i).is_unbounded()) 
    {
      unbounded = true;
      if(!min) break;
    }
    else
    {
      double w = coeff(i).diam();
      if(min ? w<d : w>d)
      {
        selected_index = i;
        d = w;
      }
    }
  }
 
  if(min && selected_index == -1)
  {
    assert(unbounded);
 
    if(among_indices.empty())
    {
      // the selected interval is the first one.
      i = 0;
    }
    else
      i = among_indices[0];
  }
 
  // The unbounded intervals are not considered if we look for the minimal diameter
  // and some bounded intervals have been found (selected_index!=-1)
  if(unbounded && (!min || selected_index == -1))
  {
    double pt = min ? -codac2::oo : codac2::oo; // keep the point less/most distant from +oo (we normalize if the lower bound is -oo)
    selected_index = i;
 
    for(; i < this->size() ; i++)
    {
      if(!among_indices.empty()
        && std::find(among_indices.begin(),among_indices.end(),i) == among_indices.end())
        continue; // avoiding this index
 
      if(coeff(i).lb() == -codac2::oo)
      {
        if(coeff(i).ub() == codac2::oo)
          if(!min)
          {
            selected_index = i;
            break;
          }
 
        if((min && (-coeff(i).ub() > pt)) || (!min && (-coeff(i).ub() < pt)))
        {
          selected_index = i;
          pt = -coeff(i).ub();
        }
      }
 
      else if(coeff(i).ub() == codac2::oo)
        if((min && (coeff(i).lb() > pt)) || (!min && (coeff(i).lb() < pt)))
        {
          selected_index = i;
          pt = coeff(i).lb();
        }
    }
  }
 
  return selected_index;
}

inline bool contains(const Matrix<double,RowsAtCompileTime,ColsAtCompileTime>& x) const
{
  return _contains(x);
}

template<typename OtherDerived>
inline bool contains(const MatrixBase<OtherDerived>& x) const
{
  return _contains(x);
}

template<typename T>
inline bool _contains(const T& x) const
{
  assert_release(x.size() == this->size());
 
  if(this->is_empty())
    return false;
 
  for(Index i = 0 ; i < this->rows() ; i++)
    for(Index j = 0 ; j < this->cols() ; j++)
      if(!(*this)(i,j).contains(x(i,j)))
        return false;
 
  return true;
}

inline bool interior_contains(const Matrix<double,RowsAtCompileTime,ColsAtCompileTime>& x) const
{
  return _interior_contains(x);
}

template<typename OtherDerived>
inline bool interior_contains(const MatrixBase<OtherDerived>& x) const
{
  return _interior_contains(x);
}

template<typename T>
inline bool _interior_contains(const T& x) const
{
  assert_release(x.size() == this->size());
 
  if(this->is_empty())
    return false;
 
  for(Index i = 0 ; i < this->rows() ; i++)
    for(Index j = 0 ; j < this->cols() ; j++)
      if(!(*this)(i,j).interior_contains(x(i,j)))
        return false;
 
  return true;
}

inline bool is_unbounded() const
{
  if(this->is_empty()) return false;
  for(Index i = 0 ; i < this->rows() ; i++)
    for(Index j = 0 ; j < this->cols() ; j++)
      if((*this)(i,j).is_unbounded())
        return true;
  return false;
}

inline bool is_degenerated() const
{
  for(Index i = 0 ; i < this->rows() ; i++)
    for(Index j = 0 ; j < this->cols() ; j++)
      if(!(*this)(i,j).is_degenerated())
        return false;
  return true;
}

inline bool is_flat() const
{
  if(this->is_empty()) return true;
  for(Index i = 0 ; i < this->rows() ; i++)
    for(Index j = 0 ; j < this->cols() ; j++)
      if((*this)(i,j).is_degenerated()) // don't use diam() because of roundoff
        return true;
  return false;
}

inline bool intersects(const Matrix<codac2::Interval,RowsAtCompileTime,ColsAtCompileTime>& x) const
{
  return _intersects(x);
}

template<typename OtherDerived>
inline bool intersects(const MatrixBase<OtherDerived>& x) const
{
  return _intersects(x);
}

template<typename OtherDerived>
inline bool _intersects(const MatrixBase<OtherDerived>& x) const
{
  assert_release(this->size() == x.size());
 
  if(this->is_empty())
    return false;
 
  for(Index i = 0 ; i < this->rows() ; i++)
    for(Index j = 0 ; j < this->cols() ; j++)
      if(!(*this)(i,j).intersects(x(i,j)))
        return false;
 
  return true;
}

inline bool is_disjoint(const Matrix<codac2::Interval,RowsAtCompileTime,ColsAtCompileTime>& x) const
{
  return _is_disjoint(x);
}

template<typename OtherDerived>
inline bool is_disjoint(const MatrixBase<OtherDerived>& x) const
{
  return _is_disjoint(x);
}

template<typename OtherDerived>
inline bool _is_disjoint(const MatrixBase<OtherDerived>& x) const
{
  assert_release(this->size() == x.size());
 
  if(this->is_empty())
    return true;
 
  for(Index i = 0 ; i < this->rows() ; i++)
    for(Index j = 0 ; j < this->cols() ; j++)
      if((*this)(i,j).is_disjoint(x(i,j)))
        return true;
 
  return false;
}

inline bool overlaps(const Matrix<codac2::Interval,RowsAtCompileTime,ColsAtCompileTime>& x) const
{
  return _overlaps(x);
}

template<typename OtherDerived>
inline bool overlaps(const MatrixBase<OtherDerived>& x) const
{
  return _overlaps(x);
}

template<typename OtherDerived>
inline bool _overlaps(const MatrixBase<OtherDerived>& x) const
{
  assert_release(this->size() == x.size());
 
  if(this->is_empty())
    return false;
 
  for(Index i = 0 ; i < this->rows() ; i++)
    for(Index j = 0 ; j < this->cols() ; j++)
      if(!(*this)(i,j).overlaps(x(i,j)))
        return false;
 
  return true;
}

inline bool is_subset(const Matrix<codac2::Interval,RowsAtCompileTime,ColsAtCompileTime>& x) const
{
  return _is_subset(x);
}

template<typename OtherDerived>
inline bool is_subset(const MatrixBase<OtherDerived>& x) const
{
  return _is_subset(x);
}

template<typename T>
inline bool _is_subset(const T& x) const
{
  assert_release(this->size() == x.size());
 
  if(this->is_empty())
    return true;
 
  for(Index i = 0 ; i < this->rows() ; i++)
    for(Index j = 0 ; j < this->cols() ; j++)
      if(!(*this)(i,j).is_subset(x(i,j)))
        return false;
 
  return true;
}

inline bool is_strict_subset(const Matrix<codac2::Interval,RowsAtCompileTime,ColsAtCompileTime>& x) const
{
  return _is_strict_subset(x);
}

template<typename OtherDerived>
inline bool is_strict_subset(const MatrixBase<OtherDerived>& x) const
{
  return _is_strict_subset(x);
}

template<typename T>
inline bool _is_strict_subset(const T& x) const
{
  assert_release(this->size() == x.size());
 
  if(this->is_empty())
    return true;
 
  if(!is_subset(x))
    return false;
 
  for(Index i = 0 ; i < this->rows() ; i++)
    for(Index j = 0 ; j < this->cols() ; j++)
      if((*this)(i,j).is_strict_subset(x(i,j)))
        return true;
 
  return false;
}

inline bool is_interior_subset(const Matrix<codac2::Interval,RowsAtCompileTime,ColsAtCompileTime>& x) const
{
  return _is_interior_subset(x);
}

template<typename OtherDerived>
inline bool is_interior_subset(const MatrixBase<OtherDerived>& x) const
{
  return _is_interior_subset(x);
}

template<typename OtherDerived>
inline bool _is_interior_subset(const MatrixBase<OtherDerived>& x) const
{
  assert_release(this->size() == x.size());
 
  if(this->is_empty())
    return true;
 
  for(Index i = 0 ; i < this->rows() ; i++)
    for(Index j = 0 ; j < this->cols() ; j++)
      if(!(*this)(i,j).is_interior_subset(x(i,j)))
        return false;
 
  return true;
}

inline bool is_strict_interior_subset(const Matrix<codac2::Interval,RowsAtCompileTime,ColsAtCompileTime>& x) const
{
  return _is_strict_interior_subset(x);
}

template<typename OtherDerived>
inline bool is_strict_interior_subset(const MatrixBase<OtherDerived>& x) const
{
  return _is_strict_interior_subset(x);
}

template<typename OtherDerived>
inline bool _is_strict_interior_subset(const MatrixBase<OtherDerived>& x) const
{
  assert_release(this->size() == x.size());
 
  if(this->is_empty())
    return true;
 
  for(Index i = 0 ; i < this->rows() ; i++)
    for(Index j = 0 ; j < this->cols() ; j++)
      if(!(*this)(i,j).is_strict_interior_subset(x(i,j)))
        return false;
 
  return true;
}

inline bool is_superset(const Matrix<codac2::Interval,RowsAtCompileTime,ColsAtCompileTime>& x) const
{
  return _is_superset(x);
}

template<typename OtherDerived>
inline bool is_superset(const MatrixBase<OtherDerived>& x) const
{
  return _is_superset(x);
}

template<typename OtherDerived>
inline bool _is_superset(const MatrixBase<OtherDerived>& x) const
{
  assert_release(this->size() == x.size());
 
  if(this->is_empty())
    return false;
 
  for(Index i = 0 ; i < this->rows() ; i++)
    for(Index j = 0 ; j < this->cols() ; j++)
      if(!x(i,j).is_subset((*this)(i,j)))
        return false;
 
  return true;
}

inline bool is_strict_superset(const Matrix<codac2::Interval,RowsAtCompileTime,ColsAtCompileTime>& x) const
{
  return _is_strict_superset(x);
}

template<typename OtherDerived>
inline bool is_strict_superset(const MatrixBase<OtherDerived>& x) const
{
  return _is_strict_superset(x);
}

template<typename OtherDerived>
inline bool _is_strict_superset(const MatrixBase<OtherDerived>& x) const
{
  assert_release(this->size() == x.size());
 
  if(this->is_empty())
    return false;
 
  if(!is_superset(x))
    return false;
 
  for(Index i = 0 ; i < this->rows() ; i++)
    for(Index j = 0 ; j < this->cols() ; j++)
      if(x(i,j).is_strict_subset((*this)(i,j)))
        return true;
 
  return false;
}

inline bool is_bisectable() const
{
  for(Index i = 0 ; i < this->rows() ; i++)
    for(Index j = 0 ; j < this->cols() ; j++)
      if((*this)(i,j).is_bisectable())
        return true;
  return false;
}

inline bool has_integer_bounds() const
{
  for(Index i = 0 ; i < this->rows() ; i++)
    for(Index j = 0 ; j < this->cols() ; j++)
      if(!(*this)(i,j).has_integer_bounds())
        return false;
  return true;
}