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.md at
// the top level directory of deal.II.
//
// ---------------------------------------------------------------------

#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
      {
        using type = Number;
      };


      template <typename Number>
      struct RemoveComplexWrapper<std::complex<Number>>
      {
        using type = typename RemoveComplexWrapper<Number>::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


      using scalar_type = ScalarType;


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


      using complex_type = std::complex<real_type>;


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

      using derivative_type = 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;


      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>
    {
      using scalar_type = ScalarType;

      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