codac2_matrices.h

Go to the documentation of this file.

 
#pragma once
 
/* On MSVC arm64, Eigen does not enable NEON and its matrix product kernel
 * fails for mixed scalar types (e.g., Matrix * IntervalMatrix).
 * Use the pmadd path to support these products.
 */
#if defined(_MSC_VER) && defined(_M_ARM64)
#ifndef EIGEN_HAS_SINGLE_INSTRUCTION_MADD
#define EIGEN_HAS_SINGLE_INSTRUCTION_MADD
#endif
#endif
 
/* Eigen 5.0.0 overstates byte alignment in std::assume_aligned, triggering
 * UBSan on macOS/Clang. Disable this hint without changing SIMD loads/stores.
 */
#if defined(__APPLE__) && defined(__clang__)
#ifndef EIGEN_ASSUME_ALIGNED
#define EIGEN_ASSUME_ALIGNED(PTR, ALIGN_BYTES)
#endif
#endif
 
#include <type_traits>
#include "codac2_Interval.h"
#include "codac2_Interval_operations.h"
#include "codac2_assert.h"
 
namespace Eigen
{
  template<int R,int C>
  concept IsVectorOrRow = (C == 1 || R == 1);
 
  template<typename Scalar>
  concept IsIntervalDomain = std::is_same_v<Scalar,codac2::Interval>;
}
 
#define EIGEN_MATRIXBASE_PLUGIN "codac2_MatrixBase_addons_include.h"
#define EIGEN_MATRIX_PLUGIN "codac2_Matrix_addons_include.h"
// The option #define EIGEN_INITIALIZE_MATRICES_BY_ZERO is not used because
// of the interval cases: an IntervalMatrix would have its components initialized
// to 0 instead of [-oo,oo].
 
#ifndef EIGEN_NO_DEBUG
/* Disables Eigen's assertions if defined.
 * Not defined by default, unless the NDEBUG macro is defined
 * (this is a standard C++ macro which disables all asserts).
 * https://eigen.tuxfamily.org/dox/TopicPreprocessorDirectives.html
 */
//#define EIGEN_NO_DEBUG // uncomment to disable Eigen's assertions
#endif
 
#include <Eigen/Core>
#include <Eigen/Dense>
 
namespace Eigen
{
  template<> struct NumTraits<codac2::Interval>
   : NumTraits<double> // permits to get the epsilon, dummy_precision, lowest, highest functions
  {
    typedef codac2::Interval Real;
    typedef codac2::Interval NonInteger;
    typedef codac2::Interval Nested;
    typedef codac2::Interval Scalar;
    typedef double RealScalar;
   
    enum {
      IsComplex = 0,
      IsInteger = 0,
      IsSigned = 1,
      RequireInitialization = 1,
      ReadCost = 1,
      AddCost = 3,
      MulCost = 3
    };
  };
 
  template<typename BinOp>
  struct ScalarBinaryOpTraits<codac2::Interval,double,BinOp>
  { typedef codac2::Interval ReturnType; };
 
  template<typename BinOp>
  struct ScalarBinaryOpTraits<double,codac2::Interval,BinOp>
  { typedef codac2::Interval ReturnType; };
}
 
namespace codac2
{
  using Eigen::Dynamic;
  
  inline const Interval& conj(const Interval& x) { return x; }
  inline const Interval& real(const Interval& x) { return x; }
  inline Interval imag(const Interval&)          { return 0.; }
  //inline Interval abs(const Interval& x)         { return codac2::abs(x); }
  inline Interval abs2(const Interval& x)        { return codac2::sqr(x); }
 
  template<typename Scalar,int R,int C>
  using Mat = Eigen::Matrix<Scalar,R,C>;
 
  template<int R,int C>
  struct is_interval_based<Eigen::Matrix<double,R,C>> : std::false_type {};
 
  template<int R,int C>
  struct is_interval_based<Eigen::Matrix<Interval,R,C>> : std::true_type {};
 
  template<typename T,int R,int C>
  struct is_ctc<Eigen::Matrix<T,R,C>> : std::false_type {};
 
  template<typename T,int R,int C>
  struct is_sep<Eigen::Matrix<T,R,C>> : std::false_type {};
 
  template<typename T>
  concept is_matrix_base_double =
    requires {
      typename Eigen::internal::traits<std::decay_t<T>>::Scalar;
      requires std::same_as<
        typename Eigen::internal::traits<std::decay_t<T>>::Scalar,
        double
      >;
    };
}
 
namespace codac2
{
  template<typename OtherDerived>
  inline auto abs(const Eigen::MatrixBase<OtherDerived>& x)
  {
    using M = Eigen::MatrixBase<OtherDerived>;
    Eigen::Matrix<typename M::Scalar,M::RowsAtCompileTime,M::ColsAtCompileTime> a(x.rows(),x.cols());
 
    for(Index i = 0 ; i < x.rows() ; i++)
      for(Index j = 0 ; j < x.cols() ; j++)
      {
        if constexpr(std::is_same_v<typename M::Scalar,double>)
          a(i,j) = fabs(x(i,j));
        else
          a(i,j) = abs(x(i,j));
      }
 
    return a;
  }

  template<typename OtherDerived>
    requires (!Eigen::IsIntervalDomain<typename OtherDerived::Scalar>)
  inline auto floor(const Eigen::MatrixBase<OtherDerived>& x)
  {
    return x.array().floor().matrix();
  }

  template<typename OtherDerived>
    requires (!Eigen::IsIntervalDomain<typename OtherDerived::Scalar>)
  inline auto ceil(const Eigen::MatrixBase<OtherDerived>& x)
  {
    return x.array().ceil().matrix();
  }

  template<typename OtherDerived>
    requires (!Eigen::IsIntervalDomain<typename OtherDerived::Scalar>)
  inline auto round(const Eigen::MatrixBase<OtherDerived>& x)
  {
    return x.array().round().matrix();
  }

  inline Eigen::IOFormat codac_row_fmt()
  {
    return Eigen::IOFormat(Eigen::StreamPrecision, Eigen::DontAlignCols, " ", "", "", "", "[ ", " ]");
  }

  inline Eigen::IOFormat codac_vector_fmt()
  {
    return Eigen::IOFormat(Eigen::StreamPrecision, Eigen::DontAlignCols, "", " ; ", "", "", "[ ", " ]");
  }
    
  inline Eigen::IOFormat codac_matrix_fmt()
  {
    return Eigen::IOFormat(Eigen::StreamPrecision, 0, " , ", "\n", "[ ", " ]", "[", "]");
  }
}