codac2_Matrix_addons_IntervalMatrixBase.h

Go to the documentation of this file.

 
// does not work.. template<typename OtherDerived1,typename OtherDerived2,typename U=Scalar,int R=RowsAtCompileTime,int C=ColsAtCompileTime>
// does not work..   requires (
// does not work..        IsIntervalDomain<U>
// does not work..     && (!IsIntervalDomain<typename OtherDerived1::Scalar>)
// does not work..     && (!IsIntervalDomain<typename OtherDerived2::Scalar>)
// does not work..     && OtherDerived1::RowsAtCompileTime == R
// does not work..     && OtherDerived2::RowsAtCompileTime == R
// does not work..     && OtherDerived1::ColsAtCompileTime == C
// does not work..     && OtherDerived2::ColsAtCompileTime == C)
// does not work.. Matrix(MatrixBase<OtherDerived1> lb, MatrixBase<OtherDerived2> ub)
// does not work..   : Matrix<U,R,C>(lb)

template<int R=RowsAtCompileTime,int C=ColsAtCompileTime>
  requires IsIntervalDomain<Scalar>
Matrix(const Matrix<double,R,C>& lb, const Matrix<double,R,C>& ub)
  : Matrix<Scalar,R,C>(lb)
{
  assert_release(lb.size() == ub.size());
 
  for(Index i = 0 ; i < this->size() ; i++)
  {
    auto& lbi = *(this->data()+i);
    const auto& ubi = *(ub.data()+i);
 
    if(lbi.lb() > ubi)
    {
      set_empty();
      break;
    }
 
    else
      lbi |= ubi;
  }
}

template<int R=RowsAtCompileTime,int C=ColsAtCompileTime>
  requires IsIntervalDomain<Scalar>
Matrix(int r, int c, const double bounds[][2])
  : Matrix<Scalar,R,C>(r,c)
{
  assert_release(r > 0 && c > 0);
 
  Index k = 0;
  for(Index i = 0 ; i < this->rows() ; i++)
    for(Index j = 0 ; j < this->cols() ; j++)
    {
      (*this)(i,j) = codac2::Interval(bounds[k][0],bounds[k][1]);
      k++;
    }
  assert_release(k == this->size() && "incorrect array size");
}

inline auto& init()
{
  for(Index i = 0 ; i < this->size() ; i++)
    (this->data()+i)->init();
  return *this;
}

template<int R=RowsAtCompileTime,int C=ColsAtCompileTime,typename OtherDerived>
  requires IsIntervalDomain<Scalar>
inline bool operator==(const MatrixBase<OtherDerived>& x) const
{
  return operator==(x.eval().template cast<codac2::Interval>());
}

inline void set_empty()
  requires IsIntervalDomain<Scalar>
{
  this->init(codac2::Interval::empty());
}

inline auto& inflate(double r)
  requires IsIntervalDomain<Scalar>
{
  assert_release(r >= 0.);
 
  for(Index i = 0 ; i < this->size() ; i++)
    (this->data()+i)->inflate(r);
  return *this;
}

template<typename OtherDerived>
  requires IsIntervalDomain<Scalar>
inline auto& inflate(const MatrixBase<OtherDerived>& r)
{
  assert_release(this->size() == r.size());
  assert_release(r.min_coeff() >= 0.);
 
  for(Index i = 0 ; i < this->rows() ; i++)
    for(Index j = 0 ; j < this->cols() ; j++)
      (*this)(i,j).inflate(r(i,j));
  return *this;
}

template<typename OtherDerived>
inline auto& operator&=(const MatrixBase<OtherDerived>& x)
{
  assert_release(this->size() == x.size());
 
  if constexpr(std::is_same_v<typename MatrixBase<OtherDerived>::Scalar,codac2::Interval>)
  {
    if(x.is_empty())
    {
      set_empty();
      return *this;
    }
  }
  
  for(Index i = 0 ; i < this->rows() ; i++)
    for(Index j = 0 ; j < this->cols() ; j++)
      (*this)(i,j) &= x(i,j);
  return *this;
}

template<typename OtherDerived>
inline auto& operator|=(const MatrixBase<OtherDerived>& x)
{
  assert_release(this->size() == x.size());
 
  if constexpr(std::is_same_v<typename MatrixBase<OtherDerived>::Scalar,codac2::Interval>)
  {
    if(x.is_empty())
      return *this;
  }
 
  for(Index i = 0 ; i < this->rows() ; i++)
    for(Index j = 0 ; j < this->cols() ; j++)
      (*this)(i,j) |= x(i,j);
  return *this;
}

template<typename OtherDerived>
inline auto operator&(const MatrixBase<OtherDerived>& x) const
{
  auto y = *this;
  return y &= x;
}

template<typename OtherDerived>
inline auto operator|(const MatrixBase<OtherDerived>& x) const
{
  auto y = *this;
  return y |= x.template cast<codac2::Interval>();
}

template<int R=RowsAtCompileTime,int C=ColsAtCompileTime>
  requires IsIntervalDomain<Scalar>
inline static auto empty(Index r, Index c)
{
  assert_release(r >= 0 && c >= 0);
  return Matrix<Scalar,R,C>::constant(r,c,codac2::Interval::empty());
}

template<int R=RowsAtCompileTime,int C=ColsAtCompileTime>
  requires IsIntervalDomain<Scalar>
inline auto bisect(Index i, float ratio = 0.49) const
{
  assert_release(i >= 0 && i < this->size());
  assert_release((this->data()+i)->is_bisectable());
  assert_release(codac2::Interval(0,1).interior_contains(ratio));
 
  auto p = std::make_pair(*this,*this);
  auto pi = (this->data()+i)->bisect(ratio);
  *(p.first.data()+i) = pi.first;
  *(p.second.data()+i) = pi.second;
  return p;
}

template<int R=RowsAtCompileTime,int C=ColsAtCompileTime>
  requires IsIntervalDomain<Scalar>
inline auto bisect_largest(float ratio = 0.49) const
{
  return bisect(this->max_diam_index(), ratio);
}