coin-or / CppAD

A C++ Algorithmic Differentiation Package: Home Page
https://cppad.readthedocs.io
Other
469 stars 95 forks source link

Add support for Eigen 3.4.0 #140

Closed benmwebb closed 2 years ago

benmwebb commented 2 years ago

We provide a specialization of Eigen::NumTraits for CppAD::AD classes. However, in Eigen 3.4.0 as a result of commit https://gitlab.com/libeigen/eigen/-/commit/bde674164 Eigen (in Core/MathFunctionsImpl.h) now tries to use NumTraits::infinity() and NumTraits::quiet_NaN(), which our specialization does not provide. Add these two methods, and also allow Eigen to call std::isinf() and std::isnan() on AD classes. Closes #139.

CLAassistant commented 2 years ago

CLA assistant check
All committers have signed the CLA.

bradbell commented 2 years ago

This looks great. Thanks.

I have created a branch corresponding to the pull request and have fixed some documentation problems; see https://github.com/coin-or/CppAD/tree/pull_140

The documentation is currently in omhelp, but I plan to automatically convert it to sphinx to make it easier for others to use. I have a convertor, but it needs more work and I am currently working on other things.

One thing, the test_more version of num_limits.cpp is not the one in the documentation. You can see which source file is used to create any documentation section by checking the bottom of the documentation section. If you check https://coin-or.github.io/CppAD/doc/num_limits.cpp.htm you will see the following text at the bottom: Input File: example/general/num_limits.cpp

Another thing, I have avoided adding template specializations for AD types in the std namespace. See https://en.cppreference.com/w/cpp/language/extending_std where the following text appears:

Specializations of std::numeric_limits must define all members declared static const (until C++11)static constexpr (since C++11) in the primary template, in such a way thatey are usable as integral constant expressions.

I do not think that AD types correspond to integral constant expressions, but perhaps they do ?

I think it will work to have the definition of isinf and isnan in the CppAD namespace ?

benmwebb commented 2 years ago

I don't think it would work to put the isinf and isnan overloads in the CppAD namespace, since they are called by Eigen deep inside its Core/MathFunctions.h in the Eigen::internal namespace:

template<typename T>
EIGEN_DEVICE_FUNC
typename internal::enable_if<(!internal::is_integral<T>::value)&&(!NumTraits<T>::IsComplex),bool>::type
isinf_impl(const T& x)
{
  #if defined(EIGEN_GPU_COMPILE_PHASE)
    return (::isinf)(x);
  #elif EIGEN_USE_STD_FPCLASSIFY
    using std::isinf;
    return isinf EIGEN_NOT_A_MACRO (x);
  #else
    return x>NumTraits<T>::highest() || x<NumTraits<T>::lowest();
  #endif
}

(the EIGEN_USE_STD_FPCLASSIFY case is the one that's called, in my testing)

Eigen isn't going to look in the CppAD namespace to find isinf, surely. It might work though if we define isinf in the Eigen::internal namespace, if that seems cleaner to you?

bradbell commented 2 years ago

I see the problem, Eigen is using ::isinf or std::isinf instead of just isinf, so it is requiring a change to to standard namespace. I had also missed the definitions in the std namespace that you added at the end. I will change pull_140 to specialize in the standard namespace as you suggested.

bradbell commented 2 years ago

@benmwebb I have changed cppad_eigen.hpp back to the way you had it. Please test the branch https://github.com/coin-or/CppAD/tree/pull_140

and see if it works for you. If so, I will merge it into the master.

bradbell commented 2 years ago

I have not heard back on your testing and I have a bug fix to merge into the master, so I first merged this pull request using the pull_140 branch.

benmwebb commented 2 years ago

I see the problem, Eigen is using ::isinf or std::isinf instead of just isinf, so it is requiring a change to to standard namespace

Right, I wasn't super happy poking around in std either but could find no other solution. This is also what the libMesh folks did to make Eigen 3.4 work with their dual number class: libMesh/MetaPhysicL#21

Please test the branch https://github.com/coin-or/CppAD/tree/pull_140

@bradbell I just tested master and... it nearly works, in that my simple test case passes now but unfortunately not the full code that originally failed. Looks like Eigen also tries to call std::isfinite on AD classes. But that's an easy fix; I'll add a second PR for it.