Re: [eigen] Using custom scalar types.

[Bob Carpenter <carp@xxxxxxxxxxx>]

On 4/28/13 5:44 AM, Brad Bell wrote:
> I was able to fix this problem by defining a partial specialization for
>      template<typename Scalar, bool IsInteger>
>      struct significant_decimals_default_impl
> which is defined in Eigen/src/Core/IO.h. The partial specialization can be in the change set
>      https://projects.coin-or.org/CppAD/changeset/2792#file10
>
> Perhaps there is some other way (without defining a cast from AD<Base> to int) that the same thing can be accomplished
> and that is documented and hence part of the Eigen API ?

Thanks for bringing this up!

Here's what we did in the past for Stan's (http://mc-stan.org/)
auto-diff variables, which indeed involves the nasty ceil(stan::agrad::var)
operation (which we define with derivative 0):

namespace Eigen {
  namespace internal {
    ...
    template<>
    struct significant_decimals_default_impl<stan::agrad::var,false>
    {
      static inline int run()
      {
        using std::ceil;
        return cast<double,int>(ceil(-log(NumTraits<stan::agrad::var>::epsilon().val())
                                     /log(10.0)));
      }
    };
    ...

But having read your question, I realize we should've just
plugged in std::numeric_limits<double>::epsilon() instead
of the NumTraits epsilon().  That's because our NumTraits is
defined as:

namespace Eigen {
 ...
 template <>
 struct NumTraits<stan::agrad::var>  {
    ...
    typedef stan::agrad::var Real;
    ...
    inline static Real epsilon() {
      return std::numeric_limits<double>::epsilon();
    }
  ...

There's an implicit constructor stan::agrad::var(int);

I just changed our code to cut out the middelman:

    template<>
    struct significant_decimals_default_impl<stan::agrad::var,false>
    {
      static inline int run()
      {
        using std::ceil;
        using std::log;
        return cast<double,int>(ceil(-log(std::numeric_limits<double>::epsilon())
                                     / log(10.0)));
      }
    };

Our execution targets don't actually use printing of matrices,
but it could've caused a memory leak on the C++ side for an
API user if the body of this function was executed repeatedly.

- Bob


> On 04/22/2013 07:15 PM, Brad Bell wrote:
> > If I compile the program:
> > ====================
> > # include <cppad/cppad.hpp>
> > # include <cppad/example/cppad_eigen.hpp>
> >
> > int main()
> > {
> > typedef Eigen::Matrix< CppAD::AD<double>, 1, 1> MatrixAD;
> > MatrixAD X;
> > X(0, 0) = 1;
> >
> > std::cout << X << std::endl;
> > return 0;
> > }
> > ====================
> > using eigen-3.1.3 with the trunk for the package
> > http://www.coin-or.org/CppAD/
> > I get the error message
> > .../include/Eigen/src/Core/IO.h:132:97: error: no matching function for call to ‘ceil(CppAD::AD<double>)’
> > Looking at the corresponding code I see the code
> > return cast<RealScalar,int>(ceil(-log(NumTraits<RealScalar>::epsilon())/log(RealScalar(10))));
> >
> > On the page
> > http://eigen.tuxfamily.org/dox-devel/TopicCustomizingEigen.html
> > under the heading
> > Using custom scalar types
> > I see mention of
> > "define the math functions that makes sense for your type."
> > and perhaps ceil could be considered one such function. But I do not see mention of conversion from the scalar type to
> > int, which also seems to be required ?
> >
> > The problem is that, for AD types, I think that implicit conversion to int is dangerous because the AD information is
> > lost and hence CppAD provides an explicit function for this purpose
> > http://www.coin-or.org/CppAD/Doc/integer.xml
> >
> > Perhaps it would be OK if one restricts the casting operator to be explicit, but as I understand this is a c++11
> > extension and not available in c++98.