ad_drivers.h

// ---------------------------------------------------------------------
//
// Copyright (C) 2017 - 2018 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

DEAL_II_DISABLE_EXTRA_DIAGNOSTICS
#  include <adolc/internal/usrparms.h>
DEAL_II_ENABLE_EXTRA_DIAGNOSTICS

#endif // DEAL_II_WITH_ADOLC

#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,
      typename std::enable_if<ADNumberTraits<ADNumberType>::type_code ==
                              NumberTypes::adolc_taped>::type>
    {
      using tape_index = unsigned short;

      using tape_buffer_sizes = unsigned int;
    }; // struct Types


    template <typename ADNumberType>
    struct Numbers<
      ADNumberType,
      typename std::enable_if<ADNumberTraits<ADNumberType>::type_code ==
                              NumberTypes::adolc_taped>::type>
    {
      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,
      typename std::enable_if<ADNumberTraits<ADNumberType>::type_code ==
                              NumberTypes::adolc_taped>::type>
    {
      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,
      typename std::enable_if<ADNumberTraits<ADNumberType>::type_code ==
                              NumberTypes::adolc_taped>::type>
    {
      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,
      typename std::enable_if<ADNumberTraits<ADNumberType>::type_code ==
                              NumberTypes::adolc_taped>::type>
    {
      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,
      typename std::enable_if<ADNumberTraits<ADNumberType>::type_code ==
                                NumberTypes::sacado_rad ||
                              ADNumberTraits<ADNumberType>::type_code ==
                                NumberTypes::sacado_rad_dfad>::type>
    {
      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,
      typename std::enable_if<ADNumberTraits<ADNumberType>::type_code ==
                                NumberTypes::adolc_tapeless ||
                              ADNumberTraits<ADNumberType>::type_code ==
                                NumberTypes::sacado_dfad ||
                              ADNumberTraits<ADNumberType>::type_code ==
                                NumberTypes::sacado_dfad_dfad>::type>
    {
      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