Prev Next

@(@\newcommand{\W}[1]{ \; #1 \; } \newcommand{\R}[1]{ {\rm #1} } \newcommand{\B}[1]{ {\bf #1} } \newcommand{\D}[2]{ \frac{\partial #1}{\partial #2} } \newcommand{\DD}[3]{ \frac{\partial^2 #1}{\partial #2 \partial #3} } \newcommand{\Dpow}[2]{ \frac{\partial^{#1}}{\partial {#2}^{#1}} } \newcommand{\dpow}[2]{ \frac{ {\rm d}^{#1}}{{\rm d}\, {#2}^{#1}} }@)@This is cppad-20221105 documentation. Here is a link to its current documentation .
Enable Use of Eigen Linear Algebra Package with CppAD

Syntax
# include <cppad/example/cppad_eigen.hpp>
Purpose
Enables the use of the eigen linear algebra package with the type AD<Base> ; see custom scalar types .

Example
The files eigen_array.cpp and eigen_det.cpp contain an example and test of this include file. They return true if they succeed and false otherwise.

CppAD Declarations
First declare some items that are defined by cppad.hpp: 
namespace CppAD {
    // AD<Base>
    template <class Base> class AD;
    // numeric_limits<Float>
    template <class Float>  class numeric_limits;
}

std Declarations
Next declare some template specializations in std namespace: 
namespace std {
    template <class Base> bool isinf(const CppAD::AD<Base> &x);
    template <class Base> bool isfinite(const CppAD::AD<Base> &x);
    template <class Base> bool isnan(const CppAD::AD<Base> &x);
}

Include Eigen/Core
Next define the eigen plugin and then include Eigen/Core: 


# define EIGEN_MATRIXBASE_PLUGIN <cppad/example/eigen_plugin.hpp>
# include <Eigen/Core>

Eigen NumTraits
Eigen needs the following definitions, in the Eigen namespace, to work properly with AD<Base> scalars: 
namespace Eigen {
    template <class Base> struct NumTraits< CppAD::AD<Base> >
    {   // type that corresponds to the real part of an AD<Base> value
        typedef CppAD::AD<Base>   Real;
        // type for AD<Base> operations that result in non-integer values
        typedef CppAD::AD<Base>   NonInteger;
        //  type to use for numeric literals such as "2" or "0.5".
        typedef CppAD::AD<Base>   Literal;
        // type for nested value inside an AD<Base> expression tree
        typedef CppAD::AD<Base>   Nested;

        enum {
            // does not support complex Base types
            IsComplex             = 0 ,
            // does not support integer Base types
            IsInteger             = 0 ,
            // only support signed Base types
            IsSigned              = 1 ,
            // must initialize an AD<Base> object
            RequireInitialization = 1 ,
            // computational cost of the corresponding operations
            ReadCost              = 1 ,
            AddCost               = 2 ,
            MulCost               = 2
        };

        // machine epsilon with type of real part of x
        // (use assumption that Base is not complex)
        static CppAD::AD<Base> epsilon(void)
        {   return CppAD::numeric_limits< CppAD::AD<Base> >::epsilon(); }

        // relaxed version of machine epsilon for comparison of different
        // operations that should result in the same value
        static CppAD::AD<Base> dummy_precision(void)
        {   return 100. *
                CppAD::numeric_limits< CppAD::AD<Base> >::epsilon();
        }

        // minimum normalized positive value
        static CppAD::AD<Base> lowest(void)
        {   return CppAD::numeric_limits< CppAD::AD<Base> >::min(); }

        // maximum finite value
        static CppAD::AD<Base> highest(void)
        {   return CppAD::numeric_limits< CppAD::AD<Base> >::max(); }

        // number of decimal digits that can be represented without change.
        static int digits10(void)
        {   return CppAD::numeric_limits< CppAD::AD<Base> >::digits10; }

        // not a number
        static CppAD::AD<Base> quiet_NaN(void)
        {   return CppAD::numeric_limits< CppAD::AD<Base> >::quiet_NaN(); }

        // positive infinite value
        static CppAD::AD<Base> infinity(void)
        {   return CppAD::numeric_limits< CppAD::AD<Base> >::infinity(); }
    };
}

Eigen ScalarBinaryOpTraits 
namespace Eigen {
    // Inform Eigen that a binary operations between Base and AD<Base>
    // are allowed and thate the return type is AD<Base>
    template<typename Base, typename BinOp>
    struct ScalarBinaryOpTraits<CppAD::AD<Base>, Base, BinOp>{
        typedef CppAD::AD<Base> ReturnType;
    };
    template<typename Base, typename BinOp>
    struct ScalarBinaryOpTraits<Base, CppAD::AD<Base>, BinOp>
    {
        typedef CppAD::AD<Base> ReturnType;
    };
}

CppAD Namespace
Eigen needs the following definitions, in the CppAD namespace, to work properly with AD<Base> scalars: 
namespace CppAD {
        // functions that return references
        template <class Base> const AD<Base>& conj(const AD<Base>& x)
        {   return x; }
        template <class Base> const AD<Base>& real(const AD<Base>& x)
        {   return x; }

        // functions that return values (note abs is defined by cppad.hpp)
        template <class Base> AD<Base> imag(const AD<Base>& x)
        {   return CppAD::AD<Base>(0.); }
        template <class Base> AD<Base> abs2(const AD<Base>& x)
        {   return x * x; }
}

eigen_vector
The class CppAD::eigen_vector is a wrapper for Eigen column vectors so that they are simple vectors . To be specific, it converts Eigen::Index arguments and return values to size_t. 
namespace CppAD {
    template <class Scalar>
    class eigen_vector : public Eigen::Matrix<Scalar, Eigen::Dynamic, 1> {
    private:
        // base_class
        typedef Eigen::Matrix<Scalar, Eigen::Dynamic, 1> base_class;
    public:
        // constructor
        eigen_vector(size_t n) : base_class( Eigen::Index(n) )
        { }
        eigen_vector(void) : base_class()
        { }
        // operator[]
        Scalar& operator[](size_t i)
        {   return base_class::operator[]( Eigen::Index(i) ); }
        const Scalar& operator[](size_t i) const
        {   return base_class::operator[]( Eigen::Index(i) ); }
        // size
        size_t size(void) const
        {   return size_t( base_class::size() ); }
        // resize
        void resize(size_t n)
        {   base_class::resize( Eigen::Index(n) ); }
    };
}

Include cppad.hpp 

# include <cppad/cppad.hpp>

std Definitions
The definitions below use cppad.hpp. Note that Value function can only be used with a constant parameter argument. 
namespace std {
    template <class Base> bool isinf(const CppAD::AD<Base> &x)
    {   return isinf(CppAD::Value( CppAD::Var2Par(x) ) ); }

    template <class Base> bool isfinite(const CppAD::AD<Base> &x)
    {   return isfinite(CppAD::Value( CppAD::Var2Par(x) ) ); }

    template <class Base> bool isnan(const CppAD::AD<Base> &x)
    {   return isnan(CppAD::Value( CppAD::Var2Par(x) ) ); }
}
Input File: include/cppad/example/cppad_eigen.hpp