ad_drivers.h

Go to the documentation of this file.

// ---------------------------------------------------------------------
//
// Copyright (C) 2017 - 2023 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_drivers_h
#define dealii_differentiation_ad_ad_drivers_h
 
#include <deal.II/base/config.h>
 
#include <deal.II/base/exceptions.h>
#include <deal.II/base/types.h>
 
#include <deal.II/differentiation/ad/ad_number_traits.h>
#include <deal.II/differentiation/ad/ad_number_types.h>
 
#include <deal.II/lac/full_matrix.h>
#include <deal.II/lac/vector.h>
 
#ifdef DEAL_II_WITH_ADOLC
 
#  include <adolc/internal/usrparms.h>
 
#endif // DEAL_II_WITH_ADOLC
 
#include <limits>
#include <vector>
 
 
DEAL_II_NAMESPACE_OPEN
 
 
namespace Differentiation
{
  namespace AD
  {
    DeclExceptionMsg(
      ExcRequiresADNumberSpecialization,
      "This function is called in a class that is expected to be specialized "
      "for auto-differentiable numbers.");
 
    DeclExceptionMsg(
      ExcRequiresADOLC,
      "This function is only available if deal.II is compiled with ADOL-C.");
 
    DeclException2(
      ExcSupportedDerivativeLevels,
      std::size_t,
      std::size_t,
      << "The number of derivative levels that this auto-differentiable number type supports is "
      << arg1
      << ", but to perform the intended operation the number must support at least "
      << arg2 << " levels.");
 
    template <typename ADNumberType, typename T = void>
    struct Types
    {
      using tape_index = unsigned int;
 
      using tape_buffer_sizes = unsigned int;
    }; // struct Types
 
 
    template <typename ADNumberType, typename T = void>
    struct Numbers
    {
      // Note: We stipulate in the documentation for the helper classes that the
      // valid tape index range is (invalid_tape_index,max_tape_index).
      static const typename Types<ADNumberType>::tape_index invalid_tape_index =
        0u;
 
      static const typename Types<ADNumberType>::tape_index max_tape_index =
        std::numeric_limits<typename Types<ADNumberType>::tape_index>::max();
    }; // struct Numbers
 
 
    template <typename ADNumberType, typename ScalarType, typename T = void>
    struct TapedDrivers
    {
      // This dummy class definition safely supports compilation
      // against tapeless numbers or taped number types that have
      // not yet been implemented.
 
      bool
      is_recording() const;
 
      typename Types<ADNumberType>::tape_index
      active_tape_index() const;
 
      bool
      is_registered_tape(
        const typename Types<ADNumberType>::tape_index tape_index) const;
 
      bool
      keep_independent_values() const;
 
      void
      set_tape_buffer_sizes(
        const typename Types<ADNumberType>::tape_buffer_sizes obufsize = 64 *
                                                                         1024 *
                                                                         1024,
        const typename Types<ADNumberType>::tape_buffer_sizes lbufsize = 64 *
                                                                         1024 *
                                                                         1024,
        const typename Types<ADNumberType>::tape_buffer_sizes vbufsize = 64 *
                                                                         1024 *
                                                                         1024,
        const typename Types<ADNumberType>::tape_buffer_sizes tbufsize = 64 *
                                                                         1024 *
                                                                         1024);
 
      void
      start_taping(const typename Types<ADNumberType>::tape_index tape_index,
                   const bool keep_independent_values);
 
      void
      stop_taping(
        const typename Types<ADNumberType>::tape_index active_tape_index,
        const bool                                     write_tapes_to_file);
 
      std::vector<typename Types<ADNumberType>::tape_index>
      get_registered_tape_indices() const;
 
      void
      activate_tape(const typename Types<ADNumberType>::tape_index tape_index);
 
      bool
      requires_retaping(
        const typename Types<ADNumberType>::tape_index tape_index) const;
 
      bool
      last_action_requires_retaping() const;
 
      void
      remove_tape(const typename Types<ADNumberType>::tape_index tape_index);
 
      void
      reset(const bool clear_registered_tapes);
 
      void
      print(std::ostream &stream) const;
 
      void
      print_tape_stats(
        const typename Types<ADNumberType>::tape_index tape_index,
        std::ostream &                                 stream) const;
 
      ScalarType
      value(const typename Types<ADNumberType>::tape_index active_tape_index,
            const std::vector<ScalarType> &independent_variables) const;
 
      void
      gradient(const typename Types<ADNumberType>::tape_index active_tape_index,
               const std::vector<ScalarType> &independent_variables,
               Vector<ScalarType> &           gradient) const;
 
      void
      hessian(const typename Types<ADNumberType>::tape_index active_tape_index,
              const std::vector<ScalarType> &independent_variables,
              FullMatrix<ScalarType> &       hessian) const;
 
      void
      values(const typename Types<ADNumberType>::tape_index active_tape_index,
             const unsigned int             n_dependent_variables,
             const std::vector<ScalarType> &independent_variables,
             Vector<ScalarType> &           values) const;
 
      void
      jacobian(const typename Types<ADNumberType>::tape_index active_tape_index,
               const unsigned int             n_dependent_variables,
               const std::vector<ScalarType> &independent_variables,
               FullMatrix<ScalarType> &       jacobian) const;
 
    };
 
 
 
    template <typename ADNumberType, typename ScalarType, typename T = void>
    struct TapelessDrivers
    {
      // This dummy class definition safely supports compilation
      // against taped numbers or tapeless number types that have
      // not yet been implemented.
 
      static void
      initialize_global_environment(const unsigned int n_independent_variables);
 
      void
      allow_dependent_variable_marking();
 
      void
      prevent_dependent_variable_marking();
 
      bool
      is_dependent_variable_marking_allowed() const;
 
      ScalarType
      value(const std::vector<ADNumberType> &dependent_variables) const;
 
      void
      gradient(const std::vector<ADNumberType> &independent_variables,
               const std::vector<ADNumberType> &dependent_variables,
               Vector<ScalarType> &             gradient) const;
 
      void
      hessian(const std::vector<ADNumberType> &independent_variables,
              const std::vector<ADNumberType> &dependent_variables,
              FullMatrix<ScalarType> &         hessian) const;
 
      void
      values(const std::vector<ADNumberType> &dependent_variables,
             Vector<ScalarType> &             values) const;
 
      void
      jacobian(const std::vector<ADNumberType> &independent_variables,
               const std::vector<ADNumberType> &dependent_variables,
               FullMatrix<ScalarType> &         jacobian) const;
 
    };
 
  } // namespace AD
} // namespace Differentiation
 
 
 
/* --------------------- inline and template functions --------------------- */
 
 
#ifndef DOXYGEN
 
namespace Differentiation
{
  namespace AD
  {
#  ifdef DEAL_II_WITH_ADOLC
 
    template <typename ADNumberType>
    struct Types<ADNumberType,
                 std::enable_if_t<ADNumberTraits<ADNumberType>::type_code ==
                                  NumberTypes::adolc_taped>>
    {
      using tape_index = unsigned short;
 
      using tape_buffer_sizes = unsigned int;
    }; // struct Types
 
 
    template <typename ADNumberType>
    struct Numbers<ADNumberType,
                   std::enable_if_t<ADNumberTraits<ADNumberType>::type_code ==
                                    NumberTypes::adolc_taped>>
    {
      static const typename Types<ADNumberType>::tape_index invalid_tape_index =
        0;
 
      // Note: This value is a limitation of ADOL-C, and not something that we
      // have control over. See test adolc/helper_tape_index_01.cc for
      // verification that we cannot use or exceed this value. This value is
      // defined as TBUFNUM; see
      // https://gitlab.com/adol-c/adol-c/blob/master/ADOL-C/include/adolc/internal/usrparms.h#L34
#    ifdef __clang__
      static const typename Types<ADNumberType>::tape_index max_tape_index =
        TBUFNUM;
#    else
      // For some reason, the test adolc/helper_tape_index_01 indicates that
      // ADOL-C does not reliably perform correct computations for the full
      // range of tape indices when GCC is the compiler. So we limit this number
      // according to the results of the test.
      static const typename Types<ADNumberType>::tape_index max_tape_index =
        TBUFNUM - 2;
#    endif
    }; // struct Numbers
 
 
    template <typename ADNumberType>
    struct TapedDrivers<
      ADNumberType,
      double,
      std::enable_if_t<ADNumberTraits<ADNumberType>::type_code ==
                       NumberTypes::adolc_taped>>
    {
      using scalar_type = double;
 
      TapedDrivers();
 
 
      bool
      is_recording() const;
 
      typename Types<ADNumberType>::tape_index
      active_tape_index() const;
 
      bool
      keep_independent_values() const;
 
      bool
      is_registered_tape(
        const typename Types<ADNumberType>::tape_index tape_index) const;
 
      void
      set_tape_buffer_sizes(
        const typename Types<ADNumberType>::tape_buffer_sizes in_obufsize,
        const typename Types<ADNumberType>::tape_buffer_sizes in_lbufsize,
        const typename Types<ADNumberType>::tape_buffer_sizes in_vbufsize,
        const typename Types<ADNumberType>::tape_buffer_sizes in_tbufsize);
 
      void
      start_taping(const typename Types<ADNumberType>::tape_index tape_index,
                   const bool keep_independent_values);
 
      void
      stop_taping(
        const typename Types<ADNumberType>::tape_index active_tape_index,
        const bool                                     write_tapes_to_file);
 
      std::vector<typename Types<ADNumberType>::tape_index>
      get_registered_tape_indices() const;
 
      void
      activate_tape(const typename Types<ADNumberType>::tape_index tape_index);
 
      bool
      requires_retaping(
        const typename Types<ADNumberType>::tape_index tape_index) const;
 
      bool
      last_action_requires_retaping() const;
 
      void
      remove_tape(const typename Types<ADNumberType>::tape_index tape_index);
 
      void
      reset(const bool clear_registered_tapes);
 
      void
      print(std::ostream &stream) const;
 
      void
      print_tape_stats(
        const typename Types<ADNumberType>::tape_index tape_index,
        std::ostream &                                 stream) const;
 
      scalar_type
      value(const typename Types<ADNumberType>::tape_index active_tape_index,
            const std::vector<scalar_type> &independent_variables) const;
 
      void
      gradient(const typename Types<ADNumberType>::tape_index active_tape_index,
               const std::vector<scalar_type> &independent_variables,
               Vector<scalar_type> &           gradient) const;
 
      void
      hessian(const typename Types<ADNumberType>::tape_index active_tape_index,
              const std::vector<scalar_type> &independent_variables,
              FullMatrix<scalar_type> &       hessian) const;
 
      void
      values(const typename Types<ADNumberType>::tape_index active_tape_index,
             const unsigned int              n_dependent_variables,
             const std::vector<scalar_type> &independent_variables,
             Vector<scalar_type> &           values) const;
 
      void
      jacobian(const typename Types<ADNumberType>::tape_index active_tape_index,
               const unsigned int              n_dependent_variables,
               const std::vector<scalar_type> &independent_variables,
               FullMatrix<scalar_type> &       jacobian) const;
 
    protected:
      typename Types<ADNumberType>::tape_index active_tape;
 
      bool keep_values;
 
      bool is_recording_flag;
 
      mutable std::map<typename Types<ADNumberType>::tape_index, int> status;
 
      bool use_stored_taped_buffer_sizes;
 
      typename Types<ADNumberType>::tape_buffer_sizes obufsize;
 
      typename Types<ADNumberType>::tape_buffer_sizes lbufsize;
 
      typename Types<ADNumberType>::tape_buffer_sizes vbufsize;
 
      typename Types<ADNumberType>::tape_buffer_sizes tbufsize;
    };
 
#  else
 
    template <typename ADNumberType>
    struct TapedDrivers<
      ADNumberType,
      double,
      std::enable_if_t<ADNumberTraits<ADNumberType>::type_code ==
                       NumberTypes::adolc_taped>>
    {
      using scalar_type = double;
 
      bool
      is_recording() const;
 
      typename Types<ADNumberType>::tape_index
      active_tape_index() const;
 
      bool
      keep_independent_values() const;
 
      bool
      is_registered_tape(
        const typename Types<ADNumberType>::tape_index tape_index) const;
 
      void
      set_tape_buffer_sizes(
        const typename Types<ADNumberType>::tape_buffer_sizes,
        const typename Types<ADNumberType>::tape_buffer_sizes,
        const typename Types<ADNumberType>::tape_buffer_sizes,
        const typename Types<ADNumberType>::tape_buffer_sizes);
 
      void
      start_taping(const typename Types<ADNumberType>::tape_index, const bool);
 
      void
      stop_taping(const typename Types<ADNumberType>::tape_index, const bool);
 
      std::vector<typename Types<ADNumberType>::tape_index>
      get_registered_tape_indices() const;
 
      void
      activate_tape(const typename Types<ADNumberType>::tape_index);
 
      bool
      requires_retaping(const typename Types<ADNumberType>::tape_index) const;
 
      bool
      last_action_requires_retaping() const;
 
      void
      remove_tape(const typename Types<ADNumberType>::tape_index);
 
      void
      reset(const bool);
 
      void
      print(std::ostream &stream) const;
 
      void
      print_tape_stats(const typename Types<ADNumberType>::tape_index,
                       std::ostream &) const;
 
      scalar_type
      value(const typename Types<ADNumberType>::tape_index,
            const std::vector<scalar_type> &) const;
 
      void
      gradient(const typename Types<ADNumberType>::tape_index,
               const std::vector<scalar_type> &,
               Vector<scalar_type> &) const;
 
      void
      hessian(const typename Types<ADNumberType>::tape_index,
              const std::vector<scalar_type> &,
              FullMatrix<scalar_type> &) const;
 
      void
      values(const typename Types<ADNumberType>::tape_index,
             const unsigned int,
             const std::vector<scalar_type> &,
             Vector<scalar_type> &) const;
 
      void
      jacobian(const typename Types<ADNumberType>::tape_index,
               const unsigned int,
               const std::vector<scalar_type> &,
               FullMatrix<scalar_type> &) const;
 
    };
 
#  endif // DEAL_II_WITH_ADOLC
 
    template <typename ADNumberType>
    struct TapedDrivers<
      ADNumberType,
      float,
      std::enable_if_t<ADNumberTraits<ADNumberType>::type_code ==
                       NumberTypes::adolc_taped>>
    {
      using scalar_type = float;
 
      bool
      is_recording() const;
 
      typename Types<ADNumberType>::tape_index
      active_tape_index() const;
 
      bool
      keep_independent_values() const;
 
      bool
      is_registered_tape(
        const typename Types<ADNumberType>::tape_index tape_index) const;
 
      void
      set_tape_buffer_sizes(
        const typename Types<ADNumberType>::tape_buffer_sizes obufsize,
        const typename Types<ADNumberType>::tape_buffer_sizes lbufsize,
        const typename Types<ADNumberType>::tape_buffer_sizes vbufsize,
        const typename Types<ADNumberType>::tape_buffer_sizes tbufsize);
 
      void
      start_taping(const typename Types<ADNumberType>::tape_index tape_index,
                   const bool keep_independent_values);
 
      void
      stop_taping(
        const typename Types<ADNumberType>::tape_index active_tape_index,
        const bool                                     write_tapes_to_file);
 
      std::vector<typename Types<ADNumberType>::tape_index>
      get_registered_tape_indices() const;
 
      void
      activate_tape(const typename Types<ADNumberType>::tape_index tape_index);
 
      bool
      requires_retaping(
        const typename Types<ADNumberType>::tape_index tape_index) const;
 
      bool
      last_action_requires_retaping() const;
 
      void
      remove_tape(const typename Types<ADNumberType>::tape_index tape_index);
 
      void
      reset(const bool clear_registered_tapes);
 
      void
      print(std::ostream &stream) const;
 
      void
      print_tape_stats(
        const typename Types<ADNumberType>::tape_index tape_index,
        std::ostream &                                 stream) const;
 
      scalar_type
      value(const typename Types<ADNumberType>::tape_index active_tape_index,
            const std::vector<scalar_type> &independent_variables) const;
 
      void
      gradient(const typename Types<ADNumberType>::tape_index active_tape_index,
               const std::vector<scalar_type> &independent_variables,
               Vector<scalar_type> &           gradient) const;
 
      void
      hessian(const typename Types<ADNumberType>::tape_index active_tape_index,
              const std::vector<scalar_type> &independent_variables,
              FullMatrix<scalar_type> &       hessian) const;
 
      void
      values(const typename Types<ADNumberType>::tape_index active_tape_index,
             const unsigned int              n_dependent_variables,
             const std::vector<scalar_type> &independent_variables,
             Vector<scalar_type> &           values) const;
 
      void
      jacobian(const typename Types<ADNumberType>::tape_index active_tape_index,
               const unsigned int              n_dependent_variables,
               const std::vector<scalar_type> &independent_variables,
               FullMatrix<scalar_type> &       jacobian) const;
 
    private:
      std::vector<double>
      vector_float_to_double(const std::vector<float> &in) const;
 
      TapedDrivers<ADNumberType, double> taped_driver;
    };
 
 
    // -------------   TapelessDrivers   -------------
 
 
    template <typename ADNumberType, typename ScalarType>
    struct TapelessDrivers<
      ADNumberType,
      ScalarType,
      std::enable_if_t<ADNumberTraits<ADNumberType>::type_code ==
                         NumberTypes::sacado_rad ||
                       ADNumberTraits<ADNumberType>::type_code ==
                         NumberTypes::sacado_rad_dfad>>
    {
      TapelessDrivers();
 
      static void
      initialize_global_environment(const unsigned int n_independent_variables);
 
      void
      allow_dependent_variable_marking();
 
      void
      prevent_dependent_variable_marking();
 
      bool
      is_dependent_variable_marking_allowed() const;
 
      ScalarType
      value(const std::vector<ADNumberType> &dependent_variables) const;
 
      void
      gradient(const std::vector<ADNumberType> &independent_variables,
               const std::vector<ADNumberType> &dependent_variables,
               Vector<ScalarType> &             gradient) const;
 
      void
      hessian(const std::vector<ADNumberType> &independent_variables,
              const std::vector<ADNumberType> &dependent_variables,
              FullMatrix<ScalarType> &         hessian) const;
 
      void
      values(const std::vector<ADNumberType> &dependent_variables,
             Vector<ScalarType> &             values) const;
 
      void
      jacobian(const std::vector<ADNumberType> &independent_variables,
               const std::vector<ADNumberType> &dependent_variables,
               FullMatrix<ScalarType> &         jacobian) const;
 
    private:
      bool dependent_variable_marking_safe;
    };
 
 
    template <typename ADNumberType, typename ScalarType>
    struct TapelessDrivers<
      ADNumberType,
      ScalarType,
      std::enable_if_t<ADNumberTraits<ADNumberType>::type_code ==
                         NumberTypes::adolc_tapeless ||
                       ADNumberTraits<ADNumberType>::type_code ==
                         NumberTypes::sacado_dfad ||
                       ADNumberTraits<ADNumberType>::type_code ==
                         NumberTypes::sacado_dfad_dfad>>
    {
      TapelessDrivers();
 
      static void
      initialize_global_environment(const unsigned int n_independent_variables);
 
      void
      allow_dependent_variable_marking();
 
      void
      prevent_dependent_variable_marking();
 
      bool
      is_dependent_variable_marking_allowed() const;
 
      ScalarType
      value(const std::vector<ADNumberType> &dependent_variables) const;
 
      void
      gradient(const std::vector<ADNumberType> &independent_variables,
               const std::vector<ADNumberType> &dependent_variables,
               Vector<ScalarType> &             gradient) const;
 
      void
      hessian(const std::vector<ADNumberType> &independent_variables,
              const std::vector<ADNumberType> &dependent_variables,
              FullMatrix<ScalarType> &         hessian) const;
 
      void
      values(const std::vector<ADNumberType> &dependent_variables,
             Vector<ScalarType> &             values) const;
 
      void
      jacobian(const std::vector<ADNumberType> &independent_variables,
               const std::vector<ADNumberType> &dependent_variables,
               FullMatrix<ScalarType> &         jacobian) const;
 
    private:
      bool dependent_variable_marking_safe;
    };
 
  } // namespace AD
} // namespace Differentiation
 
 
#endif // DOXYGEN
 
DEAL_II_NAMESPACE_CLOSE
 
 
#endif