ad_number_traits.h

// ---------------------------------------------------------------------
//
// Copyright (C) 2016 - 2017 by the deal.II authors
//
// This file is part of the deal.II library.
//
// The deal.II library is free software; you can use it, redistribute
// it, and/or modify it under the terms of the GNU Lesser General
// Public License as published by the Free Software Foundation; either
// version 2.1 of the License, or (at your option) any later version.
// The full text of the license can be found in the file LICENSE at
// the top level of the deal.II distribution.
//
// ---------------------------------------------------------------------

#ifndef dealii_differentiation_ad_ad_number_traits_h
#define dealii_differentiation_ad_ad_number_traits_h

#include <deal.II/base/exceptions.h>
#include <deal.II/differentiation/ad/ad_number_types.h>

#include <boost/type_traits.hpp>

#include <complex>
#include <type_traits>

DEAL_II_NAMESPACE_OPEN

namespace Differentiation
{
  namespace AD
  {

    template<typename ScalarType, enum NumberTypes ADNumberTypeCode, typename T = void>
    struct NumberTraits;



    template<typename ADNumberType, typename T = void>
    struct ADNumberTraits;


    namespace internal
    {

      // The following three classes, namely ADNumberInfoFromEnum, Marking, and
      // ExtractData, are those that need to be implemented for each new
      // auto-differentiable number type. This information is then used by
      // NumberTraits and ADNumberTraits to provide a uniform interface, as used
      // by our drivers, to the underlying number types.


      template<typename ScalarType, enum NumberTypes ADNumberTypeCode, typename = void>
      struct ADNumberInfoFromEnum;


      template<typename ADNumberType, typename T = void>
      struct Marking;


      template<typename ADNumberType, typename T = void>
      struct ExtractData;


      template <typename ADNumberTrait, typename T = void>
      struct HasRequiredADInfo;


      template <typename T>
      struct NumberType;


      template<typename Number>
      struct RemoveComplexWrapper;

    } // namespace internal


    template <typename NumberType>
    struct is_ad_number;


    template <typename NumberType, typename = void>
    struct is_taped_ad_number;


    template <typename NumberType, typename = void>
    struct is_tapeless_ad_number;


    template <typename NumberType, typename = void>
    struct is_real_valued_ad_number;


    template <typename NumberType, typename = void>
    struct is_complex_valued_ad_number;

  } // namespace AD
} // namespace Differentiation


/* --------------------------- inline and template functions and specializations ------------------------- */


#ifndef DOXYGEN


namespace Differentiation
{
  namespace AD
  {
    namespace internal
    {
      template <typename ADNumberTrait, typename>
      struct HasRequiredADInfo
        : std::false_type
      { };


      template <typename ADNumberTrait>
      struct HasRequiredADInfo <ADNumberTrait, decltype(
        (void) ADNumberTrait::type_code,
        (void) ADNumberTrait::is_taped,
        (void) std::declval<typename ADNumberTrait::real_type>(),
        (void) std::declval<typename ADNumberTrait::derivative_type>(),
  void() )> : std::true_type
      { };


      template<typename ScalarType>
      struct Marking<ScalarType, typename std::enable_if<
        std::is_arithmetic<ScalarType>::value
        >::type>
      {
        template<typename ADNumberType>
        static void
        independent_variable(const ScalarType   &in,
                             const unsigned int ,
                             const unsigned int ,
                             ADNumberType       &out)
        {
          out = in;
        }

        /*
         * Initialize the state of a dependent variable.
         */
        template<typename ADNumberType>
        static void
        dependent_variable(ADNumberType &,
                           const ScalarType &)
        {
          AssertThrow(false, ExcMessage("Floating point numbers cannot be marked as dependent variables."));
        }
      };


      template<typename ADNumberType>
      struct Marking<ADNumberType, typename std::enable_if<
        boost::is_complex<ADNumberType>::value
        >::type>
      {

        /*
         * Initialize the state of an independent variable.
         */
        template<typename ScalarType>
        static void
        independent_variable(const ScalarType   &in,
                             const unsigned int ,
                             const unsigned int ,
                             ADNumberType       &out)
        {
          AssertThrow(false, ExcMessage("Marking for complex numbers has not yet been implemented."));
          out = in;
        }

        /*
         * Initialize the state of a dependent variable.
         */
        template<typename ScalarType>
        static void
        dependent_variable(ADNumberType &,
                           const ScalarType &)
        {
          AssertThrow(false, ExcMessage("Marking for complex numbers has not yet been implemented."));
        }
      };

    } // namespace internal


    template <typename NumberType, typename>
    struct is_taped_ad_number
      : std::false_type
    {};


    template <typename NumberType, typename>
    struct is_tapeless_ad_number
      : std::false_type
    {};


    template <typename NumberType, typename>
    struct is_real_valued_ad_number
      : std::false_type
    {};


    template <typename NumberType, typename>
    struct is_complex_valued_ad_number
      : std::false_type
    {};


    template <typename NumberType>
    struct is_ad_number
      : internal::HasRequiredADInfo< ADNumberTraits<typename std::decay<NumberType>::type> >
    {};


    template <typename NumberType>
    struct is_taped_ad_number<NumberType, typename std::enable_if<
      ADNumberTraits<typename std::decay<NumberType>::type>::is_taped
      >::type>
      : std::true_type
    {};


    template <typename NumberType>
    struct is_tapeless_ad_number<NumberType, typename std::enable_if<
      ADNumberTraits<typename std::decay<NumberType>::type>::is_tapeless
      >::type>
      : std::true_type
    {};


    template <typename NumberType>
    struct is_real_valued_ad_number<NumberType, typename std::enable_if<
      ADNumberTraits<typename std::decay<NumberType>::type>::is_real_valued
      >::type>
      : std::true_type
    {};


    template <typename NumberType>
    struct is_complex_valued_ad_number<NumberType, typename std::enable_if<
      ADNumberTraits<typename std::decay<NumberType>::type>::is_complex_valued
      >::type>
      : std::true_type
    {};


    namespace internal
    {

      template<typename Number>
      struct RemoveComplexWrapper
      {
        typedef Number type;
      };


      template<typename Number>
      struct RemoveComplexWrapper<std::complex<Number> >
      {
        typedef typename RemoveComplexWrapper<Number>::type type;
      };


      template<typename NumberType>
      struct ExtractData<NumberType, typename std::enable_if<
        std::is_arithmetic<NumberType>::value
        >::type>
      {
        static const NumberType &
        value (const NumberType &x)
        {
          return x;
        }


        static unsigned int
        n_directional_derivatives (const NumberType &)
        {
          return 0;
        }


        static NumberType
        directional_derivative  (const NumberType &,
                                 const unsigned int )
        {
          return 0.0;
        }
      };



      template<typename ADNumberType>
      struct ExtractData< std::complex<ADNumberType> >
      {
        static_assert(
          Differentiation::AD::is_ad_number<ADNumberType>::value,
          "Expected an auto-differentiable number."
        );


        static std::complex<typename ADNumberTraits<ADNumberType>::scalar_type>
        value (const std::complex<ADNumberType> &x)
        {
          return std::complex<typename ADNumberTraits<ADNumberType>::scalar_type>(
                   ExtractData<ADNumberType>::value(x.real()),
                   ExtractData<ADNumberType>::value(x.imag()) );
        }


        static unsigned int
        n_directional_derivatives (const std::complex<ADNumberType> &x)
        {
          return ExtractData<ADNumberType>::n_directional_derivatives(x.real());
        }


        static std::complex<typename ADNumberTraits<ADNumberType>::derivative_type>
        directional_derivative  (const std::complex<ADNumberType> &x,
                                 const unsigned int                direction)
        {
          return std::complex<typename ADNumberTraits<ADNumberType>::derivative_type>(
                   ExtractData<ADNumberType>::directional_derivative(x.real(), direction),
                   ExtractData<ADNumberType>::directional_derivative(x.imag(), direction) );
        }
      };


      template <typename T>
      struct NumberType
      {
        template<typename F>
        static auto
        value (const F &f,
               typename std::enable_if<
               !is_ad_number<F>::value
               >::type * = nullptr)
        -> decltype (::internal::NumberType<T>::value(f))
        {
          // We call the other function defined in the numbers
          // header to take care of all of the usual cases.
          return ::internal::NumberType<T>::value(f);
        }

        template<typename F>
        static T
        value (const F &f,
               typename std::enable_if<
               is_ad_number<F>::value &&
               std::is_arithmetic<T>::value
               >::type * = nullptr)
        {
          // We recursively call this function in case the AD number is a
          // nested one. The recursion ends when the extracted value is
          // a floating point number.
          return NumberType<T>::value(ExtractData<F>::value(f));
        }

        template<typename F>
        static T
        value (const F &f,
               typename std::enable_if<
               is_ad_number<F>::value &&
               is_ad_number<T>::value
               >::type * = nullptr)
        {
          return T(f);
        }
      };

      template <typename T>
      struct NumberType< std::complex<T> >
      {
        template<typename F>
        static auto
        value (const F &f,
               typename std::enable_if<
               !is_ad_number<F>::value
               >::type * = nullptr)
        -> decltype (::internal::NumberType< std::complex<T> >::value(f))
        {
          // We call the other function defined in the numbers
          // header to take care of all of the usual cases.
          return ::internal::NumberType< std::complex<T> >::value(f);
        }


        template<typename F>
        static std::complex<T>
        value (const F &f,
               typename std::enable_if<
               is_ad_number<F>::value &&
               std::is_arithmetic<T>::value
               >::type * = nullptr)
        {
          // We recursively call this function in case the AD number is a
          // nested one. The recursion ends when the extracted value is
          // a floating point number.
          return std::complex<T>(NumberType<T>::value(ExtractData<F>::value(f)));
        }

        template <typename F>
        static std::complex<T>
        value (const std::complex<F> &f)
        {
          // Deal with the two parts of the input complex
          // number individually.
          return std::complex<T>(
                   NumberType<T>::value(f.real()),
                   NumberType<T>::value(f.imag()));
        }
      };

    } // namespace internal



    template<typename ScalarType, enum NumberTypes ADNumberTypeCode>
    struct NumberTraits<ScalarType, ADNumberTypeCode, typename std::enable_if<
      std::is_floating_point<ScalarType>::value ||
  (boost::is_complex<ScalarType>::value &&
     std::is_floating_point<typename internal::RemoveComplexWrapper<ScalarType>::type>::value)
    >::type>
    {
      static constexpr enum NumberTypes  type_code
      = ADNumberTypeCode;

      // The clang compiler does not seem to like these
      // variables being defined as constant expressions
      // (the tests <adolc|sacado>/ad_number_traits_02 will
      // fail with linking errors). However, GCC complains
      // about the use of non-constant expressions in
      // std::conditional.
#ifdef __clang__

      static const bool                 is_taped;


      static const bool                 is_tapeless;


      static const bool                 is_real_valued;


      static const bool                 is_complex_valued;


      static const unsigned int         n_supported_derivative_levels;

#else

      static constexpr bool             is_taped
      = internal::ADNumberInfoFromEnum<
        typename internal::RemoveComplexWrapper<ScalarType>::type, ADNumberTypeCode
      >::is_taped;


      static constexpr bool             is_tapeless
      = !(NumberTraits<ScalarType,ADNumberTypeCode>::is_taped);


      static constexpr bool             is_real_valued
      = (!boost::is_complex<ScalarType>::value);


      static constexpr bool             is_complex_valued
      = !(NumberTraits<ScalarType,ADNumberTypeCode>::is_real_valued);


      static constexpr unsigned int     n_supported_derivative_levels
      = internal::ADNumberInfoFromEnum<
        typename internal::RemoveComplexWrapper<ScalarType>::type, ADNumberTypeCode
      >::n_supported_derivative_levels;

#endif


      typedef ScalarType                scalar_type;


      typedef typename internal::ADNumberInfoFromEnum<
      typename internal::RemoveComplexWrapper<ScalarType>::type,ADNumberTypeCode
      >::real_type                      real_type;


      typedef std::complex<real_type>   complex_type;


      typedef typename std::conditional<
      is_real_valued,real_type,complex_type
      >::type                            ad_type;

      typedef typename std::conditional<
      is_real_valued,
      typename internal::ADNumberInfoFromEnum<
      typename internal::RemoveComplexWrapper<ScalarType>::type,ADNumberTypeCode
      >::derivative_type,
      std::complex<typename internal::ADNumberInfoFromEnum<
      typename internal::RemoveComplexWrapper<ScalarType>::type,ADNumberTypeCode
      >::derivative_type>
      >::type          derivative_type;


      static scalar_type
      get_scalar_value(const ad_type &x)
    {
      // Some tricky conversion cases to consider here:
      // - Nested AD numbers
      // - std::complex<double> --> std::complex<float>
      //   e.g. when ScalarType = float and ADNumberTypeCode = adolc_taped
      // Therefore, we use the internal casting mechanism
      // provided by the internal::NumberType struct.
      return internal::NumberType<scalar_type>::value(internal::ExtractData<ad_type>::value(x));
    }


    static derivative_type
    get_directional_derivative(const ad_type      &x,
                               const unsigned int  direction)
    {
      return internal::ExtractData<ad_type>::directional_derivative(x, direction);
    }


    static unsigned int
    n_directional_derivatives(const ad_type &x)
    {
      return internal::ExtractData<ad_type>::n_directional_derivatives(x);
    }


    static_assert(
      (is_real_valued==true ?
       std::is_same<ad_type,real_type>::value :
       std::is_same<ad_type,complex_type>::value),
      "Incorrect template type selected for ad_type");

    static_assert(
      (is_complex_valued==true ?
       boost::is_complex<scalar_type>::value :
       true),
      "Expected a complex float_type");

    static_assert(
      (is_complex_valued==true ?
       boost::is_complex<ad_type>::value :
       true),
      "Expected a complex ad_type");

    };

#ifdef __clang__

    template<typename ScalarType, enum NumberTypes ADNumberTypeCode>
    const bool
    NumberTraits<ScalarType, ADNumberTypeCode, typename std::enable_if<
    std::is_floating_point<ScalarType>::value ||
    (boost::is_complex<ScalarType>::value &&
     std::is_floating_point<typename internal::RemoveComplexWrapper<ScalarType>::type>::value)
    >::type>::is_taped
      = internal::ADNumberInfoFromEnum<
        typename internal::RemoveComplexWrapper<ScalarType>::type, ADNumberTypeCode
        >::is_taped;


    template<typename ScalarType, enum NumberTypes ADNumberTypeCode>
    const bool
    NumberTraits<ScalarType, ADNumberTypeCode, typename std::enable_if<
    std::is_floating_point<ScalarType>::value ||
    (boost::is_complex<ScalarType>::value &&
     std::is_floating_point<typename internal::RemoveComplexWrapper<ScalarType>::type>::value)
    >::type>::is_tapeless
      = !(NumberTraits<ScalarType,ADNumberTypeCode>::is_taped);


    template<typename ScalarType, enum NumberTypes ADNumberTypeCode>
    const bool
    NumberTraits<ScalarType, ADNumberTypeCode, typename std::enable_if<
    std::is_floating_point<ScalarType>::value ||
    (boost::is_complex<ScalarType>::value &&
     std::is_floating_point<typename internal::RemoveComplexWrapper<ScalarType>::type>::value)
    >::type>::is_real_valued
      = (!boost::is_complex<ScalarType>::value);


    template<typename ScalarType, enum NumberTypes ADNumberTypeCode>
    const bool
    NumberTraits<ScalarType, ADNumberTypeCode, typename std::enable_if<
    std::is_floating_point<ScalarType>::value ||
    (boost::is_complex<ScalarType>::value &&
     std::is_floating_point<typename internal::RemoveComplexWrapper<ScalarType>::type>::value)
    >::type>::is_complex_valued
      = !(NumberTraits<ScalarType,ADNumberTypeCode>::is_real_valued);


    template<typename ScalarType, enum NumberTypes ADNumberTypeCode>
    const unsigned int
    NumberTraits<ScalarType, ADNumberTypeCode, typename std::enable_if<
    std::is_floating_point<ScalarType>::value ||
    (boost::is_complex<ScalarType>::value &&
     std::is_floating_point<typename internal::RemoveComplexWrapper<ScalarType>::type>::value)
    >::type>::n_supported_derivative_levels
      = internal::ADNumberInfoFromEnum<
        typename internal::RemoveComplexWrapper<ScalarType>::type, ADNumberTypeCode
        >::n_supported_derivative_levels;

#endif


    template<typename ScalarType>
    struct ADNumberTraits<ScalarType, typename std::enable_if<
      std::is_arithmetic<ScalarType>::value
      >::type>
    {
      typedef ScalarType                scalar_type;

      static ScalarType
      get_directional_derivative(const ScalarType   &/*x*/,
                                 const unsigned int  /*direction*/)
      {
        // If the AD drivers are correctly implemented then we should not get here.
        // This is essentially a dummy for when the ADNumberTypeCode for the original
        // AD number (from which one is getting a derivative >= 2) is one that specified
        // Adol-C taped and tapeless numbers, or a non-nested Sacado number.
        AssertThrow(false, ExcMessage("Floating point numbers have no directional derivatives."));
        return 0.0;
      }
    };

  } // namespace AD
} // namespace Differentiation

#endif // DOXYGEN


namespace numbers
{

  template<typename ADNumberType>
  bool
  is_nan (const typename std::enable_if<
          Differentiation::AD::is_ad_number<ADNumberType>::value,
          ADNumberType>::type &x)
  {
    return is_nan(Differentiation::AD::ADNumberTraits<ADNumberType>::get_value(x));
  }

} // namespace numbers


DEAL_II_NAMESPACE_CLOSE

#endif