eigen_in_interval.h

// #include <benchmark/benchmark.h>
#include <Eigen/Core>
#include <iostream>

// using namespace Eigen;

template<typename Middle, typename LhsScalar, typename RhsScalar>
struct scalar_bracket;
template<typename Middle, typename LhsScalar, typename RhsScalar>
struct Eigen::internal::functor_traits<scalar_bracket<Middle, LhsScalar, RhsScalar>> {
  enum {
    Cost = (NumTraits<LhsScalar>::AddCost+NumTraits<RhsScalar>::AddCost)/2
            +NumTraits<Middle>::AddCost,
    PacketAccess = false
  };
};
template<typename Lhs, typename Rhs, typename Middle>
struct Eigen::internal::result_of<scalar_bracket<Lhs, Middle, Rhs>(Lhs, Middle, Rhs)> {
  typedef bool type;
};

template<typename Middle, typename LhsScalar, typename RhsScalar>
struct scalar_bracket {
  typedef bool result_type;
  inline scalar_bracket() {}
  inline scalar_bracket(const scalar_bracket&) {}
  inline bool operator()(const LhsScalar& a, const Middle& b, const RhsScalar& c) const {
    return a < b && b < c;
  }
};

template <typename Input>
using InputScalar = typename Input::Scalar;
template <typename Input>
using OutputScalar = typename Eigen::CwiseNullaryOp<Eigen::internal::scalar_constant_op<typename Input::Scalar>, typename Eigen::DenseBase<Input>::PlainObject>;
template <typename Input>
using PlainerObject = typename Eigen::DenseBase<Input>::PlainObject;

using Principal = typename Eigen::Matrix<double, -1, -1>;
using Scalar = typename Principal::Scalar;
using UsableScalar = typename Eigen::CwiseNullaryOp<Eigen::internal::scalar_constant_op<Scalar>, Eigen::DenseBase<Principal>::PlainObject>;
using PlainObject = typename Eigen::DenseBase<Principal>::PlainObject;
// CwiseBinaryOp<scalar_bracket, const Principal, const UsableScalar> flup(const Principal& p, const double val) {
//   return CwiseBinaryOp<scalar_bracket, const Principal, const UsableScalar> (p, PlainObject::Constant(p.rows(), p.cols(), val));
// }

Eigen::CwiseTernaryOp<scalar_bracket<Scalar, Scalar, Scalar>, const UsableScalar, const Principal, const UsableScalar>
bracket(const Principal& p, Scalar const low, Scalar const up) {
  return  Eigen::CwiseTernaryOp<scalar_bracket<Scalar, Scalar, Scalar>, const UsableScalar, const Principal, const UsableScalar>(
    PlainObject::Constant(p.rows(), p.cols(), low),
    p,
    PlainObject::Constant(p.rows(), p.cols(), up)
    );
}

template <typename Input>
Eigen::CwiseTernaryOp<scalar_bracket<InputScalar<Input>, InputScalar<Input>, InputScalar<Input>>, const OutputScalar<Input>, const Input, const OutputScalar<Input>>
templated_bracket(const Input& p, InputScalar<Input> const low, InputScalar<Input> const up) {
  return  Eigen::CwiseTernaryOp<scalar_bracket<InputScalar<Input>, InputScalar<Input>, InputScalar<Input>>, const OutputScalar<Input>, const Input, const OutputScalar<Input>>(
    PlainerObject<Input>::Constant(p.rows(), p.cols(), low),
    p,
    PlainerObject<Input>::Constant(p.rows(), p.cols(), up)
    );
}

int main(){
  Eigen::Matrix<double, 3, 3> vals = Eigen::MatrixXd::Random(3, 3);
  // std::cout << vals << '\n';
  std::cout << vals - vals.transpose() + vals*vals.transpose() << '\n';
  [[maybe_unused]] auto flags = templated_bracket(vals - vals.transpose() + vals*vals.transpose(), -.5, .5).eval();
  std::cout << flags;
}





// inline Eigen::CwiseTernaryOp<
//       scalar_bracket,
//       const Principal,
//       const UsableScalar,
//       const UsableScalar
//       >
//   bracket(const Principal& p, Scalar const low, Scalar const up) {
//   return Eigen::CwiseTernaryOp<scalar_bracket, const Principal, const UsableScalar, const UsableScalar>(p, low, up);
// }

// Eigen::Matrix<double, -1, 1>::Scalar k(Eigen::Matrix<double, -1, -1> const& m) {
//   return m(1,1);
// }
// CwiseBinaryOp<scalar_bracket, const Principal, const CwiseNullaryOp<internal::scalar_constant_op<Scalar>, PlainObject>>