symengine_optimizer.h

Go to the documentation of this file.

// ---------------------------------------------------------------------
//
// Copyright (C) 2020 - 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_sd_symengine_optimizer_h
#define dealii_differentiation_sd_symengine_optimizer_h
 
#include <deal.II/base/config.h>
 
#ifdef DEAL_II_WITH_SYMENGINE
 
// Low level
#  include <symengine/basic.h>
#  include <symengine/dict.h>
#  include <symengine/symengine_exception.h>
#  include <symengine/symengine_rcp.h>
 
// Optimization
#  include <symengine/lambda_double.h>
#  include <symengine/visitor.h>
#  ifdef HAVE_SYMENGINE_LLVM
#    include <symengine/llvm_double.h>
#  endif
 
#  include <deal.II/base/exceptions.h>
#  include <deal.II/base/logstream.h>
#  include <deal.II/base/utilities.h>
 
#  include <deal.II/differentiation/sd/symengine_number_types.h>
#  include <deal.II/differentiation/sd/symengine_number_visitor_internal.h>
#  include <deal.II/differentiation/sd/symengine_scalar_operations.h>
#  include <deal.II/differentiation/sd/symengine_tensor_operations.h>
#  include <deal.II/differentiation/sd/symengine_utilities.h>
 
#  include <boost/serialization/split_member.hpp>
#  include <boost/type_traits.hpp>
 
#  include <algorithm>
#  include <map>
#  include <memory>
#  include <type_traits>
#  include <utility>
#  include <vector>
 
 
DEAL_II_NAMESPACE_OPEN
 
 
namespace Differentiation
{
  namespace SD
  {
    DeclExceptionMsg(ExcSymEngineLLVMNotAvailable,
                     "SymEngine has not been built with LLVM support.");
 
    DeclExceptionMsg(ExcSymEngineLLVMReturnTypeNotSupported,
                     "The SymEngine LLVM optimizer does not (yet) support the "
                     "selected return type.");
 
    // Forward declarations
    template <typename ReturnType>
    class BatchOptimizer;
 
 
    enum class OptimizerType
    {
      dictionary,
      lambda,
      llvm
    };
 
 
    template <class StreamType>
    inline StreamType &
    operator<<(StreamType &s, OptimizerType o)
    {
      if (o == OptimizerType::dictionary)
        s << "dictionary";
      else if (o == OptimizerType::lambda)
        s << "lambda";
      else if (o == OptimizerType::llvm)
        s << "llvm";
      else
        {
          Assert(false, ExcMessage("Unknown optimization method."));
        }
 
      return s;
    }
 
 
    enum class OptimizationFlags : unsigned char
    {
      optimize_default = 0,
      optimize_cse = 0x0001,
      optimize_aggressive = 0x0002,
      optimize_all = optimize_cse | optimize_aggressive
    };
 
 
    // This operator exists since if it did not then the result of the bit-or
    // <tt>operator |</tt> would be an integer which would in turn trigger a
    // compiler warning when we tried to assign it to an object of type
    // OptimizationFlags.
    inline OptimizationFlags
    operator|(const OptimizationFlags f1, const OptimizationFlags f2)
    {
      return static_cast<OptimizationFlags>(static_cast<unsigned int>(f1) |
                                            static_cast<unsigned int>(f2));
    }
 
 
    inline OptimizationFlags &
    operator|=(OptimizationFlags &f1, const OptimizationFlags f2)
    {
      f1 = f1 | f2;
      return f1;
    }
 
 
    // This operator exists since if it did not then the result of the bit-or
    // <tt>operator |</tt> would be an integer which would in turn trigger a
    // compiler warning when we tried to assign it to an object of type
    // OptimizationFlags.
    inline OptimizationFlags
    operator&(const OptimizationFlags f1, const OptimizationFlags f2)
    {
      return static_cast<OptimizationFlags>(static_cast<unsigned int>(f1) &
                                            static_cast<unsigned int>(f2));
    }
 
 
    inline OptimizationFlags &
    operator&=(OptimizationFlags &f1, const OptimizationFlags f2)
    {
      f1 = f1 & f2;
      return f1;
    }
 
 
    namespace internal
    {
      inline bool
      use_symbolic_CSE(const enum OptimizationFlags &flags)
      {
        return static_cast<int>(flags & OptimizationFlags::optimize_cse);
      }
 
      inline int
      get_LLVM_optimization_level(const enum OptimizationFlags &flags)
      {
        // With the LLVM compiler there exists the opportunity to tune
        // the level of optimizations performed during compilation.
        // By default SymEngine sets this at "opt_level=2", which one
        // presumes targets -O2. Here we are a bit more specific about
        // want we want it to do:
        // - Normal compilation: -02 (default settings)
        // - Aggressive mode: -03 (the whole lot!)
        // In theory we could also target
        // - Debug mode: -O0 (no optimizations)
        // but this doesn't make much sense since SymEngine is a
        // tested external library.
        const bool use_agg_opt =
          static_cast<int>(flags & OptimizationFlags::optimize_aggressive);
        const int opt_level = (use_agg_opt ? 3 : 2);
        return opt_level;
      }
    } // namespace internal
 
 
    template <class StreamType>
    inline StreamType &
    operator<<(StreamType &s, OptimizationFlags o)
    {
      s << " OptimizationFlags|";
      if (static_cast<unsigned int>(o & OptimizationFlags::optimize_cse))
        s << "cse|";
 
      // LLVM optimization level
      s << "-O" + std::to_string(internal::get_LLVM_optimization_level(o)) +
             "|";
 
      return s;
    }
 
 
    namespace internal
    {
      template <typename ReturnType, typename T = void>
      struct DictionaryOptimizer;
 
 
      template <typename ReturnType, typename T = void>
      struct LambdaOptimizer;
 
 
#  ifdef HAVE_SYMENGINE_LLVM
      template <typename ReturnType, typename T = void>
      struct LLVMOptimizer;
#  endif // HAVE_SYMENGINE_LLVM
 
 
      template <typename ReturnType, typename Optimizer, typename T = void>
      struct OptimizerHelper;
 
 
#  ifndef DOXYGEN
 
 
      /* ----------- Specializations for the Optimizers ----------- */
 
 
      // A helper struct to type trait detection for the optimizers that
      // will be defined next.
      template <typename ReturnType_, typename T = void>
      struct SupportedOptimizerTypeTraits
      {
        static const bool is_supported = false;
 
        using ReturnType = void;
      };
 
 
 
      // Specialization for arithmetic types
      template <typename ReturnType_>
      struct SupportedOptimizerTypeTraits<
        ReturnType_,
        std::enable_if_t<std::is_arithmetic<ReturnType_>::value>>
      {
        static const bool is_supported = true;
 
        using ReturnType =
          typename std::conditional<std::is_same<ReturnType_, float>::value,
                                    float,
                                    double>::type;
      };
 
 
 
      // Specialization for complex arithmetic types
      template <typename ReturnType_>
      struct SupportedOptimizerTypeTraits<
        ReturnType_,
        std::enable_if_t<
          boost::is_complex<ReturnType_>::value &&
          std::is_arithmetic<typename ReturnType_::value_type>::value>>
      {
        static const bool is_supported = true;
 
        using ReturnType = typename std::conditional<
          std::is_same<ReturnType_, std::complex<float>>::value,
          std::complex<float>,
          std::complex<double>>::type;
      };
 
 
 
      template <typename ReturnType_>
      struct DictionaryOptimizer<ReturnType_,
                                 std::enable_if_t<SupportedOptimizerTypeTraits<
                                   ReturnType_>::is_supported>>
      {
        using ReturnType =
          typename SupportedOptimizerTypeTraits<ReturnType_>::ReturnType;
        using OptimizerType =
          internal::DictionarySubstitutionVisitor<ReturnType, SD::Expression>;
 
 
        static void
        initialize(OptimizerType &               optimizer,
                   const SymEngine::vec_basic &  independent_symbols,
                   const SymEngine::vec_basic &  dependent_functions,
                   const enum OptimizationFlags &optimization_flags)
        {
          const bool use_symbolic_cse = use_symbolic_CSE(optimization_flags);
          optimizer.init(independent_symbols,
                         dependent_functions,
                         use_symbolic_cse);
        }
 
 
 
        template <class Archive>
        static void
        save(Archive &          archive,
             const unsigned int version,
             OptimizerType &    optimizer)
        {
          optimizer.save(archive, version);
        }
 
 
 
        template <class Archive>
        static void
        load(Archive &          archive,
             const unsigned int version,
             OptimizerType &    optimizer,
             const SymEngine::vec_basic & /*independent_symbols*/,
             const SymEngine::vec_basic & /*dependent_functions*/,
             const enum OptimizationFlags & /*optimization_flags*/)
        {
          optimizer.load(archive, version);
        }
 
 
 
        template <typename Stream>
        static void
        print(Stream &             stream,
              const OptimizerType &optimizer,
              const bool           print_independent_symbols = false,
              const bool           print_dependent_functions = false,
              const bool           print_cse_reductions      = true)
        {
          optimizer.print(stream,
                          print_independent_symbols,
                          print_dependent_functions,
                          print_cse_reductions);
        }
      };
 
 
 
      template <typename ReturnType_>
      struct LambdaOptimizer<ReturnType_,
                             std::enable_if_t<SupportedOptimizerTypeTraits<
                               ReturnType_>::is_supported>>
      {
        using ReturnType =
          typename std::conditional<!boost::is_complex<ReturnType_>::value,
                                    double,
                                    std::complex<double>>::type;
        using OptimizerType = typename std::conditional<
          !boost::is_complex<ReturnType_>::value,
          SymEngine::LambdaRealDoubleVisitor,
          SymEngine::LambdaComplexDoubleVisitor>::type;
 
 
        static void
        initialize(OptimizerType &               optimizer,
                   const SymEngine::vec_basic &  independent_symbols,
                   const SymEngine::vec_basic &  dependent_functions,
                   const enum OptimizationFlags &optimization_flags)
        {
          const bool use_symbolic_cse = use_symbolic_CSE(optimization_flags);
          optimizer.init(independent_symbols,
                         dependent_functions,
                         use_symbolic_cse);
        }
 
 
 
        template <class Archive>
        static void
        save(Archive & /*archive*/,
             const unsigned int /*version*/,
             OptimizerType & /*optimizer*/)
        {}
 
 
        template <class Archive>
        static void
        load(Archive & /*archive*/,
             const unsigned int /*version*/,
             OptimizerType &               optimizer,
             const SymEngine::vec_basic &  independent_symbols,
             const SymEngine::vec_basic &  dependent_functions,
             const enum OptimizationFlags &optimization_flags)
        {
          initialize(optimizer,
                     independent_symbols,
                     dependent_functions,
                     optimization_flags);
        }
 
 
 
        template <typename StreamType>
        static void
        print(StreamType & /*stream*/,
              const OptimizerType & /*optimizer*/,
              const bool /*print_independent_symbols*/ = false,
              const bool /*print_dependent_functions*/ = false,
              const bool /*print_cse_reductions*/      = true)
        {
          // No built-in print function
        }
      };
 
 
 
#    ifdef HAVE_SYMENGINE_LLVM
      template <typename ReturnType_>
      struct LLVMOptimizer<
        ReturnType_,
        std::enable_if_t<std::is_arithmetic<ReturnType_>::value>>
      {
        using ReturnType =
          typename std::conditional<std::is_same<ReturnType_, float>::value,
                                    float,
                                    double>::type;
        using OptimizerType =
          typename std::conditional<std::is_same<ReturnType_, float>::value,
                                    SymEngine::LLVMFloatVisitor,
                                    SymEngine::LLVMDoubleVisitor>::type;
 
        static const bool supported_by_LLVM = true;
 
 
        static void
        initialize(OptimizerType &               optimizer,
                   const SymEngine::vec_basic &  independent_symbols,
                   const SymEngine::vec_basic &  dependent_functions,
                   const enum OptimizationFlags &optimization_flags)
        {
          const int opt_level = get_LLVM_optimization_level(optimization_flags);
          const bool use_symbolic_cse = use_symbolic_CSE(optimization_flags);
          optimizer.init(independent_symbols,
                         dependent_functions,
                         use_symbolic_cse,
                         opt_level);
        }
 
 
 
        template <class Archive>
        static void
        save(Archive &archive,
             const unsigned int /*version*/,
             OptimizerType &optimizer)
        {
          const std::string llvm_compiled_function = optimizer.dumps();
          archive &         llvm_compiled_function;
        }
 
 
 
        template <class Archive>
        static void
        load(Archive &archive,
             const unsigned int /*version*/,
             OptimizerType &optimizer,
             const SymEngine::vec_basic & /*independent_symbols*/,
             const SymEngine::vec_basic & /*dependent_functions*/,
             const enum OptimizationFlags & /*optimization_flags*/)
        {
          std::string llvm_compiled_function;
          archive &   llvm_compiled_function;
          optimizer.loads(llvm_compiled_function);
        }
 
 
 
        template <typename StreamType>
        static void
        print(StreamType & /*stream*/,
              const OptimizerType & /*optimizer*/,
              const bool /*print_independent_symbols*/ = false,
              const bool /*print_dependent_functions*/ = false,
              const bool /*print_cse_reductions*/      = true)
        {
          // No built-in print function
        }
      };
 
 
      // There is no LLVM optimizer built with complex number support.
      // So we fall back to the LambdaDouble case as a type (required
      // at compile time), but offer no implementation. We expect that
      // the calling class does not create this type: This can be done by
      // checking the `supported_by_LLVM` flag.
      template <typename ReturnType_>
      struct LLVMOptimizer<
        ReturnType_,
        std::enable_if_t<
          boost::is_complex<ReturnType_>::value &&
          std::is_arithmetic<typename ReturnType_::value_type>::value>>
      {
        // Since there is no working implementation, these are dummy types
        // that help with templating in the calling function.
        using ReturnType = typename LambdaOptimizer<ReturnType_>::ReturnType;
        using OptimizerType =
          typename LambdaOptimizer<ReturnType_>::OptimizerType;
 
        static const bool supported_by_LLVM = false;
 
 
        static void
        initialize(OptimizerType & /*optimizer*/,
                   const SymEngine::vec_basic & /*independent_symbols*/,
                   const SymEngine::vec_basic & /*dependent_functions*/,
                   const enum OptimizationFlags & /*optimization_flags*/)
        {
          AssertThrow(false, ExcNotImplemented());
        }
 
 
 
        template <class Archive>
        static void
        save(Archive & /*archive*/,
             const unsigned int /*version*/,
             OptimizerType & /*optimizer*/)
        {
          AssertThrow(false, ExcNotImplemented());
        }
 
 
 
        template <class Archive>
        static void
        load(Archive & /*archive*/,
             const unsigned int /*version*/,
             OptimizerType & /*optimizer*/,
             const SymEngine::vec_basic & /*independent_symbols*/,
             const SymEngine::vec_basic & /*dependent_functions*/,
             const enum OptimizationFlags & /*optimization_flags*/)
        {
          AssertThrow(false, ExcNotImplemented());
        }
 
 
 
        template <typename StreamType>
        static void
        print(StreamType & /*stream*/,
              const OptimizerType & /*optimizer*/,
              const bool /*print_independent_symbols*/ = false,
              const bool /*print_dependent_functions*/ = false,
              const bool /*print_cse_reductions*/      = true)
        {
          AssertThrow(false, ExcNotImplemented());
        }
      };
#    endif // HAVE_SYMENGINE_LLVM
 
 
      /* ----------- Specializations for OptimizerHelper ----------- */
 
 
      template <typename ReturnType, typename Optimizer>
      struct OptimizerHelper<
        ReturnType,
        Optimizer,
        std::enable_if_t<
          std::is_same<ReturnType, typename Optimizer::ReturnType>::value>>
      {
        static void
        initialize(typename Optimizer::OptimizerType *optimizer,
                   const SymEngine::vec_basic &       independent_symbols,
                   const SymEngine::vec_basic &       dependent_functions,
                   const enum OptimizationFlags &     optimization_flags)
        {
          Assert(optimizer, ExcNotInitialized());
 
          // Some optimizers don't have the same interface for
          // initialization, we filter them out through the specializations
          // of the Optimizer class
          Optimizer::initialize(*optimizer,
                                independent_symbols,
                                dependent_functions,
                                optimization_flags);
        }
 
 
 
        static void
        substitute(typename Optimizer::OptimizerType *optimizer,
                   std::vector<ReturnType> &          output_values,
                   const std::vector<ReturnType> &    substitution_values)
        {
          Assert(optimizer, ExcNotInitialized());
          optimizer->call(output_values.data(), substitution_values.data());
        }
 
 
 
        template <class Archive>
        static void
        save(Archive &                          archive,
             const unsigned int                 version,
             typename Optimizer::OptimizerType *optimizer)
        {
          Assert(optimizer, ExcNotInitialized());
 
          // Some optimizers don't have the same interface for
          // serialization, we filter them out through the specializations
          // of the Optimizer class
          Optimizer::save(archive, version, *optimizer);
        }
 
 
 
        template <class Archive>
        static void
        load(Archive &                          archive,
             const unsigned int                 version,
             typename Optimizer::OptimizerType *optimizer,
             const SymEngine::vec_basic &       independent_symbols,
             const SymEngine::vec_basic &       dependent_functions,
             const enum OptimizationFlags &     optimization_flags)
        {
          Assert(optimizer, ExcNotInitialized());
 
          // Some optimizers don't have the same interface for
          // serialization, we filter them out through the specializations
          // of the Optimizer class
          Optimizer::load(archive,
                          version,
                          *optimizer,
                          independent_symbols,
                          dependent_functions,
                          optimization_flags);
        }
 
 
 
        template <typename Stream>
        static void
        print(Stream &                           stream,
              typename Optimizer::OptimizerType *optimizer,
              const bool print_independent_symbols = false,
              const bool print_dependent_functions = false,
              const bool print_cse_reductions      = true)
        {
          Assert(optimizer, ExcNotInitialized());
 
          // Some optimizers don't have a print function, so
          // we filter them out through the specializations of
          // the Optimizer class
          Optimizer::print(stream,
                           *optimizer,
                           print_independent_symbols,
                           print_dependent_functions,
                           print_cse_reductions);
        }
      };
 
      template <typename ReturnType, typename Optimizer>
      struct OptimizerHelper<
        ReturnType,
        Optimizer,
        std::enable_if_t<
          !std::is_same<ReturnType, typename Optimizer::ReturnType>::value>>
      {
        static void
        initialize(typename Optimizer::OptimizerType *optimizer,
                   const SymEngine::vec_basic &       independent_symbols,
                   const SymEngine::vec_basic &       dependent_functions,
                   const enum OptimizationFlags &     optimization_flags)
        {
          Assert(optimizer, ExcNotInitialized());
 
          const bool use_symbolic_cse = use_symbolic_CSE(optimization_flags);
          optimizer->init(independent_symbols,
                          dependent_functions,
                          use_symbolic_cse);
        }
 
 
 
        static void
        substitute(typename Optimizer::OptimizerType *optimizer,
                   std::vector<ReturnType> &          output_values,
                   const std::vector<ReturnType> &    substitution_values)
        {
          Assert(optimizer, ExcNotInitialized());
 
          // Intermediate values to accommodate the difference in
          // value types.
          std::vector<typename Optimizer::ReturnType> int_outputs(
            output_values.size());
          std::vector<typename Optimizer::ReturnType> int_inputs(
            substitution_values.size());
 
          std::copy(substitution_values.begin(),
                    substitution_values.end(),
                    int_inputs.begin());
          optimizer->call(int_outputs.data(), int_inputs.data());
          std::copy(int_outputs.begin(),
                    int_outputs.end(),
                    output_values.begin());
        }
 
 
 
        template <class Archive>
        static void
        save(Archive &                          archive,
             const unsigned int                 version,
             typename Optimizer::OptimizerType *optimizer)
        {
          Assert(optimizer, ExcNotInitialized());
          Optimizer::save(archive, version, *optimizer);
        }
 
 
 
        template <class Archive>
        static void
        load(Archive &                          archive,
             const unsigned int                 version,
             typename Optimizer::OptimizerType *optimizer,
             const SymEngine::vec_basic &       independent_symbols,
             const SymEngine::vec_basic &       dependent_functions,
             const enum OptimizationFlags &     optimization_flags)
        {
          Assert(optimizer, ExcNotInitialized());
 
          // Some optimizers don't have the same interface for
          // serialization, we filter them out through the specializations
          // of the Optimizer class
          Optimizer::load(archive,
                          version,
                          *optimizer,
                          independent_symbols,
                          dependent_functions,
                          optimization_flags);
        }
 
 
 
        template <typename Stream>
        static void
        print(Stream &                           stream,
              typename Optimizer::OptimizerType *optimizer,
              const bool                         print_cse_reductions = true,
              const bool print_independent_symbols                    = false,
              const bool print_dependent_functions                    = false)
        {
          Assert(optimizer, ExcNotInitialized());
 
          optimizer->print(stream,
                           print_independent_symbols,
                           print_dependent_functions,
                           print_cse_reductions);
        }
      };
 
#  endif // DOXYGEN
 
 
      /* -------------------- Utility functions ---------------------- */
 
 
      template <typename NumberType,
                int rank,
                int dim,
                template <int, int, typename>
                class TensorType>
      TensorType<rank, dim, NumberType>
      tensor_evaluate_optimized(
        const TensorType<rank, dim, Expression> &symbol_tensor,
        const std::vector<NumberType> &          cached_evaluation,
        const BatchOptimizer<NumberType> &       optimizer)
      {
        TensorType<rank, dim, NumberType> out;
        for (unsigned int i = 0; i < out.n_independent_components; ++i)
          {
            const TableIndices<rank> indices(
              out.unrolled_to_component_indices(i));
            out[indices] =
              optimizer.extract(symbol_tensor[indices], cached_evaluation);
          }
        return out;
      }
 
 
      template <typename NumberType, int dim>
      SymmetricTensor<4, dim, NumberType>
      tensor_evaluate_optimized(
        const SymmetricTensor<4, dim, Expression> &symbol_tensor,
        const std::vector<NumberType> &            cached_evaluation,
        const BatchOptimizer<NumberType> &         optimizer)
      {
        SymmetricTensor<4, dim, NumberType> out;
        for (unsigned int i = 0;
             i < SymmetricTensor<2, dim>::n_independent_components;
             ++i)
          for (unsigned int j = 0;
               j < SymmetricTensor<2, dim>::n_independent_components;
               ++j)
            {
              const TableIndices<4> indices =
                make_rank_4_tensor_indices<dim>(i, j);
              out[indices] =
                optimizer.extract(symbol_tensor[indices], cached_evaluation);
            }
        return out;
      }
 
 
      template <typename NumberType, typename T>
      void
      register_functions(BatchOptimizer<NumberType> &optimizer,
                         const T &                   function)
      {
        optimizer.register_function(function);
      }
 
 
      template <typename NumberType, typename T>
      void
      register_functions(BatchOptimizer<NumberType> &optimizer,
                         const std::vector<T> &      functions)
      {
        for (const auto &function : functions)
          register_functions(optimizer, function);
      }
 
 
      template <typename NumberType, typename T, typename... Args>
      void
      register_functions(BatchOptimizer<NumberType> &optimizer,
                         const T &                   function,
                         const Args &...other_functions)
      {
        register_functions(optimizer, function);
        register_functions(optimizer, other_functions...);
      }
 
 
      template <int rank,
                int dim,
                template <int, int, typename>
                class TensorType>
      types::symbol_vector
      unroll_to_expression_vector(
        const TensorType<rank, dim, Expression> &symbol_tensor)
      {
        SD::types::symbol_vector out;
        out.reserve(symbol_tensor.n_independent_components);
        for (unsigned int i = 0; i < symbol_tensor.n_independent_components;
             ++i)
          {
            const TableIndices<rank> indices(
              symbol_tensor.unrolled_to_component_indices(i));
            out.push_back(symbol_tensor[indices].get_RCP());
          }
        return out;
      }
 
 
      template <int dim>
      types::symbol_vector
      unroll_to_expression_vector(
        const SymmetricTensor<4, dim, Expression> &symbol_tensor)
      {
        SD::types::symbol_vector out;
        out.reserve(symbol_tensor.n_independent_components);
        for (unsigned int i = 0;
             i < SymmetricTensor<2, dim>::n_independent_components;
             ++i)
          for (unsigned int j = 0;
               j < SymmetricTensor<2, dim>::n_independent_components;
               ++j)
            {
              const TableIndices<4> indices =
                make_rank_4_tensor_indices<dim>(i, j);
              out.push_back(symbol_tensor[indices].get_RCP());
            }
        return out;
      }
 
    } // namespace internal
 
 
 
    template <typename ReturnType>
    class BatchOptimizer
    {
    public:
      BatchOptimizer();
 
      BatchOptimizer(const enum OptimizerType &    optimization_method,
                     const enum OptimizationFlags &optimization_flags =
                       OptimizationFlags::optimize_all);
 
      BatchOptimizer(const BatchOptimizer &other);
 
      BatchOptimizer(BatchOptimizer &&) noexcept = default;
 
      ~BatchOptimizer() = default;
 
      void
      copy_from(const BatchOptimizer &other);
 
      template <typename Stream>
      void
      print(Stream &stream, const bool print_cse = false) const;
 
      template <class Archive>
      void
      save(Archive &archive, const unsigned int version) const;
 
      template <class Archive>
      void
      load(Archive &archive, const unsigned int version);
 
#  ifdef DOXYGEN
      template <class Archive>
      void
      serialize(Archive &archive, const unsigned int version);
#  else
      // This macro defines the serialize() method that is compatible with
      // the templated save() and load() method that have been implemented.
      BOOST_SERIALIZATION_SPLIT_MEMBER()
#  endif
 
      void
      register_symbols(const types::substitution_map &substitution_map);
 
      void
      register_symbols(const SymEngine::map_basic_basic &substitution_map);
 
      void
      register_symbols(const types::symbol_vector &symbols);
 
      void
      register_symbols(const SymEngine::vec_basic &symbols);
 
      types::symbol_vector
      get_independent_symbols() const;
 
      std::size_t
      n_independent_variables() const;
 
      void
      register_function(const Expression &function);
 
      template <int rank, int dim>
      void
      register_function(const Tensor<rank, dim, Expression> &function_tensor);
 
      template <int rank, int dim>
      void
      register_function(
        const SymmetricTensor<rank, dim, Expression> &function_tensor);
 
      void
      register_functions(const types::symbol_vector &functions);
 
      void
      register_functions(const SymEngine::vec_basic &functions);
 
      template <typename T>
      void
      register_functions(const std::vector<T> &functions);
 
      template <typename T, typename... Args>
      void
      register_functions(const T &functions, const Args &...other_functions);
 
      const types::symbol_vector &
      get_dependent_functions() const;
 
      std::size_t
      n_dependent_variables() const;
 
      void
      set_optimization_method(const enum OptimizerType &    optimization_method,
                              const enum OptimizationFlags &optimization_flags =
                                OptimizationFlags::optimize_all);
 
      enum OptimizerType
      optimization_method() const;
 
      enum OptimizationFlags
      optimization_flags() const;
 
      bool
      use_symbolic_CSE() const;
 
      void
      optimize();
 
      bool
      optimized() const;
 
      void
      substitute(const types::substitution_map &substitution_map) const;
 
      void
      substitute(const SymEngine::map_basic_basic &substitution_map) const;
 
      void
      substitute(const types::symbol_vector &   symbols,
                 const std::vector<ReturnType> &values) const;
 
      void
      substitute(const SymEngine::vec_basic &   symbols,
                 const std::vector<ReturnType> &values) const;
 
      bool
      values_substituted() const;
 
      const std::vector<ReturnType> &
      evaluate() const;
 
      ReturnType
      evaluate(const Expression &func) const;
 
      std::vector<ReturnType>
      evaluate(const std::vector<Expression> &funcs) const;
 
      template <int rank, int dim>
      Tensor<rank, dim, ReturnType>
      evaluate(const Tensor<rank, dim, Expression> &funcs) const;
 
 
      template <int rank, int dim>
      SymmetricTensor<rank, dim, ReturnType>
      evaluate(const SymmetricTensor<rank, dim, Expression> &funcs) const;
 
 
      ReturnType
      extract(const Expression &             func,
              const std::vector<ReturnType> &cached_evaluation) const;
 
 
      std::vector<ReturnType>
      extract(const std::vector<Expression> &funcs,
              const std::vector<ReturnType> &cached_evaluation) const;
 
 
      template <int rank, int dim>
      Tensor<rank, dim, ReturnType>
      extract(const Tensor<rank, dim, Expression> &funcs,
              const std::vector<ReturnType> &      cached_evaluation) const;
 
 
      template <int rank, int dim>
      SymmetricTensor<rank, dim, ReturnType>
      extract(const SymmetricTensor<rank, dim, Expression> &funcs,
              const std::vector<ReturnType> &cached_evaluation) const;
 
    private:
      enum OptimizerType method;
 
      enum OptimizationFlags flags;
 
      types::substitution_map independent_variables_symbols;
 
      types::symbol_vector dependent_variables_functions;
 
      bool
      is_valid_nonunique_dependent_variable(
        const SD::Expression &function) const;
 
      bool
      is_valid_nonunique_dependent_variable(
        const SymEngine::RCP<const SymEngine::Basic> &function) const;
 
      mutable std::vector<ReturnType> dependent_variables_output;
 
      using map_dependent_expression_to_vector_entry_t =
        std::map<SD::Expression,
                 std::size_t,
                 SD::types::internal::ExpressionKeyLess>;
 
      mutable map_dependent_expression_to_vector_entry_t map_dep_expr_vec_entry;
 
      mutable std::unique_ptr<SymEngine::Visitor> optimizer;
 
      mutable bool ready_for_value_extraction;
 
      mutable bool has_been_serialized;
 
      void
      register_scalar_function(const SD::Expression &function);
 
      void
      register_vector_functions(const types::symbol_vector &functions);
 
      void
      create_optimizer(std::unique_ptr<SymEngine::Visitor> &optimizer);
 
      void
      substitute(const std::vector<ReturnType> &substitution_values) const;
    };
 
 
 
    /* -------------------- inline and template functions ------------------ */
 
 
#  ifndef DOXYGEN
 
 
    template <typename ReturnType>
    template <typename Stream>
    void
    BatchOptimizer<ReturnType>::print(Stream &stream,
                                      const bool /*print_cse*/) const
    {
      // Settings
      stream << "Method? " << optimization_method() << '\n';
      stream << "Flags: " << optimization_flags() << '\n';
      stream << "Optimized? " << (optimized() ? "Yes" : "No") << '\n';
      stream << "Values substituted? " << values_substituted() << "\n\n";
 
      // Independent variables
      stream << "Symbols (" << n_independent_variables()
             << " independent variables):" << '\n';
      int cntr = 0;
      for (SD::types::substitution_map::const_iterator it =
             independent_variables_symbols.begin();
           it != independent_variables_symbols.end();
           ++it, ++cntr)
        {
          stream << cntr << ": " << it->first << '\n';
        }
      stream << '\n' << std::flush;
 
      // Dependent functions
      stream << "Functions (" << n_dependent_variables()
             << " dependent variables):" << '\n';
      cntr = 0;
      for (typename SD::types::symbol_vector::const_iterator it =
             dependent_variables_functions.begin();
           it != dependent_variables_functions.end();
           ++it, ++cntr)
        {
          stream << cntr << ": " << (*it) << '\n';
        }
      stream << '\n' << std::flush;
 
      // Common subexpression
      if (optimized() == true && use_symbolic_CSE() == true)
        {
          Assert(optimizer, ExcNotInitialized());
          const bool print_cse_reductions      = true;
          const bool print_independent_symbols = false;
          const bool print_dependent_functions = false;
 
          if (optimization_method() == OptimizerType::dictionary)
            {
              Assert(dynamic_cast<typename internal::DictionaryOptimizer<
                       ReturnType>::OptimizerType *>(optimizer.get()),
                     ExcMessage("Cannot cast optimizer to Dictionary type."));
 
              internal::OptimizerHelper<
                ReturnType,
                internal::DictionaryOptimizer<ReturnType>>::
                print(stream,
                      dynamic_cast<typename internal::DictionaryOptimizer<
                        ReturnType>::OptimizerType *>(optimizer.get()),
                      print_independent_symbols,
                      print_dependent_functions,
                      print_cse_reductions);
 
              stream << '\n' << std::flush;
            }
          else if (optimization_method() == OptimizerType::lambda)
            {
              Assert(dynamic_cast<typename internal::LambdaOptimizer<
                       ReturnType>::OptimizerType *>(optimizer.get()),
                     ExcMessage("Cannot cast optimizer to Lambda type."));
 
              internal::OptimizerHelper<ReturnType,
                                        internal::LambdaOptimizer<ReturnType>>::
                print(stream,
                      dynamic_cast<typename internal::LambdaOptimizer<
                        ReturnType>::OptimizerType *>(optimizer.get()),
                      print_independent_symbols,
                      print_dependent_functions,
                      print_cse_reductions);
            }
#    ifdef HAVE_SYMENGINE_LLVM
          else if (optimization_method() == OptimizerType::llvm)
            {
              Assert(dynamic_cast<typename internal::LLVMOptimizer<
                       ReturnType>::OptimizerType *>(optimizer.get()),
                     ExcMessage("Cannot cast optimizer to LLVM type."));
 
              internal::OptimizerHelper<ReturnType,
                                        internal::LLVMOptimizer<ReturnType>>::
                print(stream,
                      dynamic_cast<typename internal::LLVMOptimizer<
                        ReturnType>::OptimizerType *>(optimizer.get()),
                      print_independent_symbols,
                      print_dependent_functions,
                      print_cse_reductions);
            }
#    endif // HAVE_SYMENGINE_LLVM
          else
            {
              AssertThrow(false, ExcMessage("Unknown optimizer type."));
            }
        }
 
      if (values_substituted())
        {
          stream << "Evaluated functions:" << '\n';
          stream << std::flush;
          cntr = 0;
          for (typename std::vector<ReturnType>::const_iterator it =
                 dependent_variables_output.begin();
               it != dependent_variables_output.end();
               ++it, ++cntr)
            {
              stream << cntr << ": " << (*it) << '\n';
            }
          stream << '\n' << std::flush;
        }
    }
 
 
 
    template <typename ReturnType>
    template <class Archive>
    void
    BatchOptimizer<ReturnType>::save(Archive &          ar,
                                     const unsigned int version) const
    {
      // Serialize enum classes...
      {
        const auto m =
          static_cast<typename std::underlying_type<OptimizerType>::type>(
            method);
        ar &m;
      }
      {
        const auto f =
          static_cast<typename std::underlying_type<OptimizationFlags>::type>(
            flags);
        ar &f;
      }
 
      // Important: Independent variables must always be
      // serialized before the dependent variables.
      ar &independent_variables_symbols;
      ar &dependent_variables_functions;
 
      ar &dependent_variables_output;
      ar &map_dep_expr_vec_entry;
      ar &ready_for_value_extraction;
 
      // Mark that we've saved this class at some point.
      has_been_serialized = true;
      ar &has_been_serialized;
 
      // When we serialize the optimizer itself, we have to (unfortunately)
      // provide it with sufficient information to rebuild itself from scratch.
      // This is because only two of the three optimization classes support
      // real serialization (i.e. have save/load capability).
      const SD::types::symbol_vector symbol_vec =
        Utilities::extract_symbols(independent_variables_symbols);
      if (typename internal::DictionaryOptimizer<ReturnType>::OptimizerType
            *opt = dynamic_cast<typename internal::DictionaryOptimizer<
              ReturnType>::OptimizerType *>(optimizer.get()))
        {
          Assert(optimization_method() == OptimizerType::dictionary,
                 ExcInternalError());
          internal::OptimizerHelper<
            ReturnType,
            internal::DictionaryOptimizer<ReturnType>>::save(ar, version, opt);
        }
      else if (typename internal::LambdaOptimizer<ReturnType>::OptimizerType
                 *opt = dynamic_cast<typename internal::LambdaOptimizer<
                   ReturnType>::OptimizerType *>(optimizer.get()))
        {
          Assert(optimization_method() == OptimizerType::lambda,
                 ExcInternalError());
          internal::OptimizerHelper<
            ReturnType,
            internal::LambdaOptimizer<ReturnType>>::save(ar, version, opt);
        }
#    ifdef HAVE_SYMENGINE_LLVM
      else if (typename internal::LLVMOptimizer<ReturnType>::OptimizerType
                 *opt = dynamic_cast<typename internal::LLVMOptimizer<
                   ReturnType>::OptimizerType *>(optimizer.get()))
        {
          Assert(optimization_method() == OptimizerType::llvm,
                 ExcInternalError());
          internal::OptimizerHelper<
            ReturnType,
            internal::LLVMOptimizer<ReturnType>>::save(ar, version, opt);
        }
#    endif
      else
        {
          AssertThrow(false, ExcMessage("Unknown optimizer type."));
        }
    }
 
 
 
    template <typename ReturnType>
    template <class Archive>
    void
    BatchOptimizer<ReturnType>::load(Archive &ar, const unsigned int version)
    {
      Assert(independent_variables_symbols.empty(), ExcInternalError());
      Assert(dependent_variables_functions.empty(), ExcInternalError());
      Assert(dependent_variables_output.empty(), ExcInternalError());
      Assert(map_dep_expr_vec_entry.empty(), ExcInternalError());
      Assert(ready_for_value_extraction == false, ExcInternalError());
 
      // Deserialize enum classes...
      {
        typename std::underlying_type<OptimizerType>::type m;
        ar &                                               m;
        method = static_cast<OptimizerType>(m);
      }
      {
        typename std::underlying_type<OptimizationFlags>::type f;
        ar &                                                   f;
        flags = static_cast<OptimizationFlags>(f);
      }
 
      // Important: Independent variables must always be
      // deserialized before the dependent variables.
      ar &independent_variables_symbols;
      ar &dependent_variables_functions;
 
      ar &dependent_variables_output;
      ar &map_dep_expr_vec_entry;
      ar &ready_for_value_extraction;
 
      ar &has_been_serialized;
 
      // If we're reading in data, then create the optimizer
      // and then deserialize it.
      Assert(!optimizer, ExcInternalError());
 
      // Create and configure the optimizer
      create_optimizer(optimizer);
      Assert(optimizer, ExcNotInitialized());
 
      // When we deserialize the optimizer itself, we have to (unfortunately)
      // provide it with sufficient information to rebuild itself from scratch.
      // This is because only two of the three optimization classes support
      // real serialization (i.e. have save/load capability).
      const SD::types::symbol_vector symbol_vec =
        Utilities::extract_symbols(independent_variables_symbols);
      if (typename internal::DictionaryOptimizer<ReturnType>::OptimizerType
            *opt = dynamic_cast<typename internal::DictionaryOptimizer<
              ReturnType>::OptimizerType *>(optimizer.get()))
        {
          Assert(optimization_method() == OptimizerType::dictionary,
                 ExcInternalError());
          internal::OptimizerHelper<ReturnType,
                                    internal::DictionaryOptimizer<ReturnType>>::
            load(ar,
                 version,
                 opt,
                 Utilities::convert_expression_vector_to_basic_vector(
                   symbol_vec),
                 Utilities::convert_expression_vector_to_basic_vector(
                   dependent_variables_functions),
                 optimization_flags());
        }
      else if (typename internal::LambdaOptimizer<ReturnType>::OptimizerType
                 *opt = dynamic_cast<typename internal::LambdaOptimizer<
                   ReturnType>::OptimizerType *>(optimizer.get()))
        {
          Assert(optimization_method() == OptimizerType::lambda,
                 ExcInternalError());
          internal::OptimizerHelper<ReturnType,
                                    internal::LambdaOptimizer<ReturnType>>::
            load(ar,
                 version,
                 opt,
                 Utilities::convert_expression_vector_to_basic_vector(
                   symbol_vec),
                 Utilities::convert_expression_vector_to_basic_vector(
                   dependent_variables_functions),
                 optimization_flags());
        }
#    ifdef HAVE_SYMENGINE_LLVM
      else if (typename internal::LLVMOptimizer<ReturnType>::OptimizerType
                 *opt = dynamic_cast<typename internal::LLVMOptimizer<
                   ReturnType>::OptimizerType *>(optimizer.get()))
        {
          Assert(optimization_method() == OptimizerType::llvm,
                 ExcInternalError());
          internal::OptimizerHelper<ReturnType,
                                    internal::LLVMOptimizer<ReturnType>>::
            load(ar,
                 version,
                 opt,
                 Utilities::convert_expression_vector_to_basic_vector(
                   symbol_vec),
                 Utilities::convert_expression_vector_to_basic_vector(
                   dependent_variables_functions),
                 optimization_flags());
        }
#    endif
      else
        {
          AssertThrow(false, ExcMessage("Unknown optimizer type."));
        }
    }
 
 
 
    template <typename ReturnType>
    template <int rank, int dim>
    void
    BatchOptimizer<ReturnType>::register_function(
      const Tensor<rank, dim, Expression> &function_tensor)
    {
      Assert(optimized() == false,
             ExcMessage(
               "Cannot register functions once the optimizer is finalised."));
 
      register_vector_functions(
        internal::unroll_to_expression_vector(function_tensor));
    }
 
 
 
    template <typename ReturnType>
    template <int rank, int dim>
    void
    BatchOptimizer<ReturnType>::register_function(
      const SymmetricTensor<rank, dim, Expression> &function_tensor)
    {
      Assert(optimized() == false,
             ExcMessage(
               "Cannot register functions once the optimizer is finalised."));
 
      register_vector_functions(
        internal::unroll_to_expression_vector(function_tensor));
    }
 
 
 
    template <typename ReturnType>
    template <typename T, typename... Args>
    void
    BatchOptimizer<ReturnType>::register_functions(
      const T &functions,
      const Args &...other_functions)
    {
      internal::register_functions(*this, functions);
      internal::register_functions(*this, other_functions...);
    }
 
 
 
    template <typename ReturnType>
    template <typename T>
    void
    BatchOptimizer<ReturnType>::register_functions(
      const std::vector<T> &functions)
    {
      internal::register_functions(*this, functions);
    }
 
 
 
    template <typename ReturnType>
    template <int rank, int dim>
    Tensor<rank, dim, ReturnType>
    BatchOptimizer<ReturnType>::extract(
      const Tensor<rank, dim, Expression> &funcs,
      const std::vector<ReturnType> &      cached_evaluation) const
    {
      return internal::tensor_evaluate_optimized(funcs,
                                                 cached_evaluation,
                                                 *this);
    }
 
 
 
    template <typename ReturnType>
    template <int rank, int dim>
    Tensor<rank, dim, ReturnType>
    BatchOptimizer<ReturnType>::evaluate(
      const Tensor<rank, dim, Expression> &funcs) const
    {
      Assert(
        values_substituted() == true,
        ExcMessage(
          "The optimizer is not configured to perform evaluation. "
          "This action can only performed after substitute() has been called."));
 
      return extract(funcs, dependent_variables_output);
    }
 
 
 
    template <typename ReturnType>
    template <int rank, int dim>
    SymmetricTensor<rank, dim, ReturnType>
    BatchOptimizer<ReturnType>::extract(
      const SymmetricTensor<rank, dim, Expression> &funcs,
      const std::vector<ReturnType> &               cached_evaluation) const
    {
      return internal::tensor_evaluate_optimized(funcs,
                                                 cached_evaluation,
                                                 *this);
    }
 
 
 
    template <typename ReturnType>
    template <int rank, int dim>
    SymmetricTensor<rank, dim, ReturnType>
    BatchOptimizer<ReturnType>::evaluate(
      const SymmetricTensor<rank, dim, Expression> &funcs) const
    {
      Assert(
        values_substituted() == true,
        ExcMessage(
          "The optimizer is not configured to perform evaluation. "
          "This action can only performed after substitute() has been called."));
 
      return extract(funcs, dependent_variables_output);
    }
 
#  endif // DOXYGEN
 
  } // namespace SD
} // namespace Differentiation
 
 
DEAL_II_NAMESPACE_CLOSE
 
#endif // DEAL_II_WITH_SYMENGINE
 
#endif