Summary:  Be more flexible regarding mixing types  

Product:  Eigen  Reporter:  Gael Guennebaud <gael.guennebaud>  
Component:  Core  general  Assignee:  Nobody <eigen.nobody>  
Status:  RESOLVED FIXED  
Severity:  Unknown  CC:  carp, chtz, eigen2, gael.guennebaud, jacob.benoit.1  
Priority:    
Version:  3.0  
Hardware:  All  
OS:  All  
Whiteboard:  
Bug Depends on:  
Bug Blocks:  558  
Attachments: 

Description
Gael Guennebaud
20110526 19:37:21 UTC
I'm a bit wary of automatic type promotion in general because it's often vague what the promoted type should be. For example, PromoteType<int,float> really shouldn't be defined. Now if you tell me that you're going to be conservative and only define PromoteType<A,B> when A is a subset of B or conversely so that PromoteType<A,B>==A or B in all cases where it's defined, then I think I'm OK. yes, PromoteType will be very conservative. As I said by default PromoteType<A,B> with A!=B will be defined only for std::complex. Then it is the responsibility of the user of a custom scalar type to properly specialize PromoteType for his/her scalar type. (I'm the one who requested this feature  thanks much for taking it up and summarizing.) The context is that we have a class VAR for automatic differentiation. We've overridden all of the scalar ops in C++, cmath lib, and Boost's C99 extensions where appropriate to deal with mixing VAR and double (i.e., all double values in original functions). When used on its own VAR works fine with Eigen in all the cases we've tested from multiplications to eigenvalues. To answer Gael and Benoit's questions: 1. In all cases we're dealing with, if VAR is the autodif variable class, (VAR _op_ double) and (double _op_ VAR) should return VAR type for all _op_ types. 2. Yes, we've overridden all of the binary ops, at least those that are built into C++, cmath and most of Boost's C99 extensions. That includes atan2(VAR,double), atan2(double,VAR), and atan2(VAR,VAR), so you shouldn't need to do anything to deal with this. Although not critical (unless it's being used internally), it'd be nice to have (VAR += double) supported, but (double += VAR) doesn't make any sense. 3. In the templated probability density library I've written for the other part of our project, I use boost::math::promote_args, which has the behavior that I need in terms of promotion. It deals with integer types, too, but it's intended for contexts where int gets cast to a float or double for arithmetic purposes, which is what we need, as we're only dealing with floatingpoint and VAR matrices and operations. I should've clarified: we have defined (VAR += double) but it'd be nice if it'd work with Eigen::Matrix<VAR,...> += Eigen::Matrix<double,...>. I'd like to add an additional use case here: Boost.Units compatibility. When working scientifically, Boost.Units is a great boon. However, since Eigen is fairly typepromotionagnostic, units don't work on multiplication. Given a quantity in meters of type double (MD), the result of MD * MD is a quantity in square meters of type double. However, Eigen's multiplication routine expects it to be an MD. Blitz++ solves this problem fairly elegantly by using a specializable class for each basic operation. Double type promotion, Bob's differentiation and Boost.Units would be almost the same issue here, with the exception that the return type of unit operations do depend on the used units. I'd be willing to help and lead in the implementation when I get a few pointers. However, before deciding about any goahead, it is probably important to know about Eigen's plans for C++11  decltype would make any implementation easier by an order of magnitude. Steve, unfortunately, in order to support as many as possible compilers and platforms we don't plan to allow any C++11 features, especially in core. That said, we already have a template struct to control the returned scalar type of mixed products, e.g.: namespace Eigen{ namespace internal{ template<typename T> struct scalar_product_traits<std::complex<T>, T> { typedef std::complex<T> ReturnType; }; } } For 3.3 we should really support real+complex and the likes. Finally, some progress for scalar multiple: https://bitbucket.org/eigen/eigen/commits/d896d06d500d/ Summary: Implement generic scalar*expr and expr*scalar operator based on scalar_product_traits. This is especially useful for custom scalar types, e.g., to enable float*expr<multi_prec> without conversion. In order to go further, I propose to generalize scalar_product_traits<Lhs,Rhs> as follows: namespace Eigen { template<typename LeftScalar,typename RightScalar,typename BinaryOp> ScalarBinaryOpTraits { static const bool Defined = true/false; // or enum { Defined = 1/0 }; typedef ... ReturnType; }; } In order to be backward compatible with code playing with internal::scalar_product_traits, I propose to only make it deprecated, and make the generic ScalarBinaryOpTraits inherits scalar_product_traits (regardless of BinaryOp, i.e., not only for BinaryOp="mul"). For the record, here is the current list of binary operators: *,/,,+,min,max,pow,zeta,polygamma,atan2 some comparisons: <,<=,>,>=,==,!= and for bool only: &&,,^ Generally that sounds good to me. I'm not sure if it makes sense to let all traits inherit from internal::scalar_product_traits, however. Furthermore, could you clarify what a user would need to define to allow (e.g.) directly adding double to a custom Float type? template<> struct ScalarBinaryOpTraits<double, MyFloat, scalar_add_op> { static const bool Defined=true; typedef MyFloat ReturnType; }; Especially, shall the third parameter be scalar_add_op<MyFloat>, or will there be a new scalar_add_op<double, MyFloat> and the third parameter is actually: template<typename LeftScalar, typename RightScalar, typename template<typename,typename> BinaryOp> struct BinaryOpTraits {....}; Currently, the scalar_product_traits are referenced by the scalar_product_op, I'm not sure if we can recursively reference from the ScalarBinaryOpTraits here (could be possible, since the scalar_product_op is forward declared). Also, will we have different traits for + and  (or for == and !=, etc)? My current impl. looks like this: template<typename T, typename U> struct scalar_product_traits { enum { Defined = 0 }; }; // all other specializations of scalar_product_traits have been removed. template<typename ScalarA, typename ScalarB, typename BinaryOp> struct ScalarBinaryOpTraits #ifndef EIGEN_PARSED_BY_DOXYGEN // for backward compatibility, use the hints given by the (deprecated) internal::scalar_product_traits class. : internal::scalar_product_traits<ScalarA,ScalarB> #endif // EIGEN_PARSED_BY_DOXYGEN {}; template<typename T, typename BinaryOp> struct ScalarBinaryOpTraits<T,T,BinaryOp> { enum { Defined = 1 }; typedef T ReturnType; }; template<typename T, typename BinaryOp> struct ScalarBinaryOpTraits<T,std::complex<T>,BinaryOp> { enum { Defined = 1 }; typedef std::complex<T> ReturnType; }; template<typename T, typename BinaryOp> struct ScalarBinaryOpTraits<std::complex<T>, T,BinaryOp> { enum { Defined = 1 }; typedef std::complex<T> ReturnType; }; Then, scalar_product_op does: typedef typename ScalarBinaryOpTraits<LhsScalar,RhsScalar>::ReturnType result_type; and I also generalized scalar_sum_op to take two scalars and does the same as scalar_product_op. Finally, EIGEN_CHECK_BINARY_COMPATIBILIY(BINOP,LHS,RHS) boils down to ScalarBinaryOpTraits<LHS, RHS,BINOP>::Defined As you probably noticed, so far I did not bother about the binaryop as it is still unclear to me what should be passed. So the first question to answer is whether we really need to allows specializations based on the underlying op? If no, then let's just remove the third template parameter of ScalarBinaryOpTraits. If yes, then we have to find an elegant solution because it would be too much trouble to specialize based on the exact binary functors, there are too many variants. Only for products: scalar_product_op scalar_conj_product_op mul_assign_op scalar_multiple2_op Perhaps, we could add generic overloads of ScalarBinaryOpTraits for all of our functors that would inherits some generic tag, e.g.: ScalarBinaryOpTraits<A,B,scalar_conj_product_op<A,B> > : ScalarBinaryOpTraits<A,B,product_tag> {}; Same for the comparisons. This way we can introduce some king of hierarchy which can be fine tuned at any level... Wouldn't it be possible to make default implementations of (at least) scalar_multiple_op and mul_assign_op based on scalar_product_op? Generally, they should have the same Defined and PacketAccess, and the Cost should just differ by perhaps a ReadCost. For the OP_assign_op<T1,T2> we must also check if is_same<T1,scalar_OP_op<T1,T2>::ReturnType> is true, i.e., adding a CustomFloat to an ordinary float usually should be impossible without prior casting. ScalarBinaryOpTraits could then simply just check if the corresponding operator is defined (with some SFINAE magic, perhaps. Or we make default implementations for every operator which just contain an enum {Defined=is_same<T1,T2>::value};). And EIGEN_CHECK_BINARY_COMPATIBILIY(BINOP,LHS,RHS) is simply: ScalarBinaryOpTraits<BINOP<LHS,RHS> >::Defined Comparison operators may need some special handling, if they shall remain in the same scalar_cmp_op (which itself makes sense, and which will make implementing packetcomparison easier later). If we want to simplify defining custom type mixing, we can provide a set of macros which implement specializations of a certain set of operators each. (We can declare the macros as experimental for 3.3.0, in case we are not sure about the API). I think that we should really be able to say that: "A op B" is fine regardless of op, which is not possible with only: ScalarBinaryOpTraits<BINOP<LHS,RHS> >. Created attachment 713 [details]
Work in progress patch
For the record, and to help discussion, here is the current status of my refactoring.
Recursively calling ScalarBinaryOpTraits<A,B,OP> from OP<A,B> seems to be OK.
It should also be noted that if someone enable "A op B" for any op but that some are not supported, then the user will still get a quite nice error message. For instance, we the attached patch, it is OK for Eigen to perform real_array = complex_array because we enabled "real op complex" for any op, but the compiler still generate:
AssignmentFunctors.h:24:102: error: assigning to 'double' from incompatible type 'const std::__1::complex<double>'
Which is probably better than what we could do on our own...
Ok, if that means that as long as there are specializations for A OP B (or func(A,B)) everything works as expected and otherwise a meaningful error is generated we can't get much better. I hope I can have a closer look at your patch over the weekend ... For completeness, here is the current error message for: real_array.pow(complex_array) (yes, this is a bad exemple as this configuration is supported by the standard, and we thus should support it as well) ../eigen/Eigen/src/Core/functors/BinaryFunctors.h:254:83: error: no matching function for call to 'pow' inline Scalar operator() (const Scalar& a, const OtherScalar& b) const { return numext::pow(a, b); } ^~~~~~~~~~~ ../eigen/Eigen/src/Core/CoreEvaluators.h:503:12: note: in instantiation of member function 'Eigen::internal::scalar_binary_pow_op<double, std::__1::complex<double> >::operator()' requested here return m_functor(m_lhsImpl.coeff(index), m_rhsImpl.coeff(index)); [... very long list of template instantiations...] If we add: template<typename T> struct ScalarBinaryOpTraits<T,std::complex<T>, internal::scalar_binary_pow_op<T, std::complex<T> > > { enum { Defined=0 }; }; then we get the same error message preceded by another one, much shorter: ../eigen/Eigen/src/Core/CwiseBinaryOp.h:106:7: error: no member named 'YOU_MIXED_DIFFERENT_NUMERIC_TYPES__YOU_NEED_TO_USE_THE_CAST_METHOD_OF_MATRIXBASE_TO_CAST_NUMERIC_TYPES_EXPLICITLY' in 'Eigen::internal::static_assertion<false>' EIGEN_CHECK_BINARY_COMPATIBILIY(BinaryOp,typename Lhs::Scalar,typename Rhs::Scalar); ^ ../eigen/Eigen/src/Core/util/XprHelper.h:659:5: note: expanded from macro 'EIGEN_CHECK_BINARY_COMPATIBILIY' YOU_MIXED_DIFFERENT_NUMERIC_TYPES__YOU_NEED_TO_USE_THE_CAST_METHOD_OF_MATRIXBASE_TO_CAST_NUMERIC_TYPES_EXPLICITLY) ^ ../eigen/Eigen/src/Core/../plugins/ArrayCwiseBinaryOps.h:89:10: note: in instantiation of member function 'Eigen::CwiseBinaryOp<Eigen::internal::scalar_binary_pow_op<double, std::__1::complex<double> >, const Eigen::Array<double, 1, 1, 0, 1, 1>, const Eigen::Array<std::__1::complex<double>, 1, 1, 0, 1, 1> >::CwiseBinaryOp' requested here return CwiseBinaryOp<internal::scalar_binary_pow_op<Scalar, typename ExponentDerived::Scalar>, const Derived, const ExponentDerived>( ^ mixing.cpp:21:12: note: in instantiation of function template specialization 'Eigen::ArrayBase<Eigen::Array<double, 1, 1, 0, 1, 1> >::pow<Eigen::Array<std::__1::complex<double>, 1, 1, 0, 1, 1> >' requested here acd = ad.pow(acd); So it might still be a good idea to explicitly disable the very few remaining incompatible configurations. I think that the set of required changes it ready to be merged, see this PR for details: https://bitbucket.org/eigen/eigen/pullrequests/194 The PR has been merged, thus closing this entry.  GitLab Migration Automatic Message  This bug has been migrated to gitlab.com's GitLab instance and has been closed from further activity. You can subscribe and participate further through the new bug through this link to our GitLab instance: https://gitlab.com/libeigen/eigen/issues/279. 