adolc_number_types.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_adolc_number_types_h
#define dealii_differentiation_ad_adolc_number_types_h

#include <deal.II/base/config.h>

#include <type_traits>

DEAL_II_NAMESPACE_OPEN


namespace Differentiation
{
  namespace AD
  {
    template <typename NumberType, typename = void>
    struct is_adolc_number : std::false_type
    {};


    template <typename NumberType, typename = void>
    struct is_adolc_taped_number : std::false_type
    {};


    template <typename NumberType, typename = void>
    struct is_adolc_tapeless_number : std::false_type
    {};
  } // namespace AD
} // namespace Differentiation


DEAL_II_NAMESPACE_CLOSE


#ifdef DEAL_II_WITH_ADOLC

#  include <deal.II/base/exceptions.h>
#  include <deal.II/base/numbers.h>

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

#  include <adolc/adouble.h> // Taped double
#  include <adolc/adtl.h>    // Tapeless double
#  include <adolc/internal/adolc_settings.h>
#  include <adolc/internal/adubfunc.h> // Taped double math functions

#  include <complex>
#  include <type_traits>

DEAL_II_NAMESPACE_OPEN

DeclExceptionMsg(ExcADOLCAdvancedBranching,
                 "This function has not yet been implemented for taped Adol-C "
                 "numbers when the advanced branching feature is activated.");


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


#  ifndef DOXYGEN


namespace Differentiation
{
  namespace AD
  {
    namespace internal
    {
      template <typename ScalarType>
      struct ADNumberInfoFromEnum<
        ScalarType,
        Differentiation::AD::NumberTypes::adolc_taped,
        typename std::enable_if<
          std::is_floating_point<ScalarType>::value>::type>
      {
        static const bool is_taped = true;
        using real_type            = adouble;
        using derivative_type      = double;
        static const unsigned int n_supported_derivative_levels =
          std::numeric_limits<unsigned int>::max();
      };


      template <typename ScalarType>
      struct ADNumberInfoFromEnum<
        ScalarType,
        Differentiation::AD::NumberTypes::adolc_tapeless,
        typename std::enable_if<
          std::is_floating_point<ScalarType>::value>::type>
      {
        static const bool is_taped                              = false;
        using real_type                                         = adtl::adouble;
        using derivative_type                                   = double;
        static const unsigned int n_supported_derivative_levels = 1;
      };


      template <typename ADNumberType>
      struct Marking<
        ADNumberType,
        typename std::enable_if<
          ADNumberTraits<ADNumberType>::type_code == NumberTypes::adolc_taped &&
          ADNumberTraits<ADNumberType>::is_real_valued>::type>
      {
        using scalar_type = typename ADNumberTraits<ADNumberType>::scalar_type;

        /*
         * Initialize the state of an independent variable.
         */
        static void
        independent_variable(const scalar_type &in,
                             const unsigned int,
                             const unsigned int,
                             ADNumberType &out)
        {
          out <<= in;
        }

        /*
         * Initialize the state of a dependent variable.
         *
         * @note The second argument must be writable, so we
         * simply pass a copy instead of a non-constant reference.
         */
        static void
        dependent_variable(ADNumberType &out, ADNumberType func)
        {
          // Store the value only (strip it of all sensitivities)
          out = ADNumberTraits<ADNumberType>::get_scalar_value(func);
          // Mark as a dependent variable
          scalar_type tmp;
          func >>= tmp;
        }
      };

      template <typename ADNumberType>
      struct Marking<ADNumberType,
                     typename std::enable_if<
                       ADNumberTraits<ADNumberType>::type_code ==
                         NumberTypes::adolc_tapeless &&
                       ADNumberTraits<ADNumberType>::is_real_valued>::type>
      {
        using scalar_type = typename ADNumberTraits<ADNumberType>::scalar_type;

        /*
         * Initialize the state of an independent variable.
         */
        static void
        independent_variable(const scalar_type &in,
                             const unsigned int index,
                             const unsigned int,
                             ADNumberType &out)
        {
          // It is important that the tapeless variables have their values set
          // before defining their directional derivative index
          out = in;

          // Violating this condition when will result in an Adol-C internal
          // error. We could rather always throw here in order to provide a
          // less cryptic message.
          AssertThrow(index < adtl::getNumDir(),
                      ExcMessage(
                        "The index number of the independent variable being "
                        "marked is greater than the number of independent "
                        "variables that have been declared."));
          out.setADValue(index, 1 /*seed value for first derivative*/);
        }

        /*
         * Initialize the state of a dependent variable.
         */
        static void
        dependent_variable(ADNumberType &out, const ADNumberType &func)
        {
          // Simply transfer value with sensitivities
          out = 0.0;
          out = func;
        }
      };


      template <>
      struct ExtractData<adouble>
      {
        static double
        value(const adouble &x)
        {
          return x.getValue();
        }


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


        static double
        directional_derivative(const adouble &, const unsigned int)
        {
          AssertThrow(false,
                      ExcMessage(
                        "The derivative values for taped Adol-C numbers must be"
                        "computed through the ::gradient function."));
          return 0.0;
        }
      };


      template <>
      struct ExtractData<adtl::adouble>
      {
        static double
        value(const adtl::adouble &x)
        {
          return x.getValue();
        }


        static unsigned int
        n_directional_derivatives(const adtl::adouble &)
        {
          // This is a global function call...
          return adtl::getNumDir();
        }


        static double
        directional_derivative(const adtl::adouble &x,
                               const unsigned int   direction)
        {
          Assert(
            direction < n_directional_derivatives(x),
            ExcMessage(
              "Requested directional derivative is greater than the number "
              "registered by Adol-C."));
          return x.getADValue(direction);
        }
      };

    } // namespace internal



    template <typename ADNumberType>
    struct ADNumberTraits<
      ADNumberType,
      typename std::enable_if<std::is_same<ADNumberType, adouble>::value>::type>
      : NumberTraits<double, NumberTypes::adolc_taped>
    {
      static_assert(std::is_same<ad_type, adouble>::value,
                    "Incorrect template type selected for taped ad_type");
      static_assert(is_taped == true, "Incorrect setting for taping");
    };



    template <typename ADNumberType>
    struct ADNumberTraits<
      ADNumberType,
      typename std::enable_if<
        std::is_same<ADNumberType, std::complex<adouble>>::value>::type>
      : NumberTraits<std::complex<double>, NumberTypes::adolc_taped>
    {
      static_assert(std::is_same<ad_type, std::complex<adouble>>::value,
                    "Incorrect template type selected for taped ad_type");
      static_assert(is_taped == true, "Incorrect setting for taping");
    };



    template <typename ADNumberType>
    struct ADNumberTraits<
      ADNumberType,
      typename std::enable_if<
        std::is_same<ADNumberType, adtl::adouble>::value>::type>
      : NumberTraits<double, NumberTypes::adolc_tapeless>
    {
      static_assert(std::is_same<ad_type, adtl::adouble>::value,
                    "Incorrect template type selected for tapeless ad_type");
      static_assert(is_tapeless == true, "Incorrect setting for taping");
    };



    template <typename ADNumberType>
    struct ADNumberTraits<
      ADNumberType,
      typename std::enable_if<
        std::is_same<ADNumberType, std::complex<adtl::adouble>>::value>::type>
      : NumberTraits<std::complex<double>, NumberTypes::adolc_tapeless>
    {
      static_assert(std::is_same<ad_type, std::complex<adtl::adouble>>::value,
                    "Incorrect template type selected for tapeless ad_type");
      static_assert(is_tapeless == true, "Incorrect setting for taping");
    };



    template <>
    struct NumberTraits<adouble, NumberTypes::adolc_taped>
      : NumberTraits<typename ADNumberTraits<adouble>::scalar_type,
                     NumberTypes::adolc_taped>
    {};


    template <>
    struct NumberTraits<std::complex<adouble>, NumberTypes::adolc_taped>
      : NumberTraits<
          typename ADNumberTraits<std::complex<adouble>>::scalar_type,
          NumberTypes::adolc_taped>
    {};


    template <>
    struct NumberTraits<adtl::adouble, NumberTypes::adolc_tapeless>
      : NumberTraits<typename ADNumberTraits<adtl::adouble>::scalar_type,
                     NumberTypes::adolc_tapeless>
    {};


    template <>
    struct NumberTraits<std::complex<adtl::adouble>,
                        NumberTypes::adolc_tapeless>
      : NumberTraits<
          typename ADNumberTraits<std::complex<adtl::adouble>>::scalar_type,
          NumberTypes::adolc_tapeless>
    {};


    template <typename NumberType>
    struct is_adolc_taped_number<
      NumberType,
      typename std::enable_if<
        ADNumberTraits<typename std::decay<NumberType>::type>::type_code ==
        NumberTypes::adolc_taped>::type> : std::true_type
    {};


    template <typename NumberType>
    struct is_adolc_tapeless_number<
      NumberType,
      typename std::enable_if<
        ADNumberTraits<typename std::decay<NumberType>::type>::type_code ==
        NumberTypes::adolc_tapeless>::type> : std::true_type
    {};


    template <typename NumberType>
    struct is_adolc_number<NumberType,
                           typename std::enable_if<
                             is_adolc_taped_number<NumberType>::value ||
                             is_adolc_tapeless_number<NumberType>::value>::type>
      : std::true_type
    {};

  } // namespace AD
} // namespace Differentiation


#  endif // DOXYGEN


DEAL_II_NAMESPACE_CLOSE


#endif // DEAL_II_WITH_ADOLC

#endif