]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Add dictionary (normal + CSE) visitors for SymEngine
authorJean-Paul Pelteret <jppelteret@gmail.com>
Thu, 2 Apr 2020 09:23:57 +0000 (11:23 +0200)
committerJean-Paul Pelteret <jppelteret@gmail.com>
Mon, 4 May 2020 13:46:03 +0000 (15:46 +0200)
include/deal.II/differentiation/sd.h
include/deal.II/differentiation/sd/symengine_number_types.h
include/deal.II/differentiation/sd/symengine_number_visitor_internal.h [new file with mode: 0644]
source/differentiation/sd/CMakeLists.txt
source/differentiation/sd/symengine_number_types.cc
source/differentiation/sd/symengine_number_visitor_internal.cc [new file with mode: 0644]
source/differentiation/sd/symengine_number_visitor_internal.inst.in [new file with mode: 0644]

index 44dca06e9e2fefd2e77961176b43f45d5cd1c065..b35105fc7f38a352bad429d51cc70a7a54dd2e1c 100644 (file)
@@ -1,6 +1,6 @@
 // ---------------------------------------------------------------------
 //
-// Copyright (C) 2019 by the deal.II authors
+// Copyright (C) 2019-2020 by the deal.II authors
 //
 // This file is part of the deal.II library.
 //
@@ -41,7 +41,14 @@ namespace Differentiation
    * @ingroup auto_symb_diff
    */
   namespace SD
-  {}
+  {
+    /**
+     * This namespace defines the classes and functions that help provide a
+     * structured interface to symbolic numbers and operations.
+     */
+    namespace internal
+    {}
+  } // namespace SD
 } // namespace Differentiation
 
 DEAL_II_NAMESPACE_CLOSE
index 246f8884d9b5469c6a9338087861d7b6a14b94aa..a20bed3afac7153807caa470297fac0da913cad7 100644 (file)
@@ -709,6 +709,25 @@ namespace Differentiation
       Expression
       substitute(const types::substitution_map &substitution_values) const;
 
+      /**
+       * Perform substitution of all symbols found in this object's @p expression
+       * that match a key in the @p substitution_values map.
+       *
+       * This function is like the one above, but takes in a SymEngine map
+       * (one that maps a `SymEngine::RCP<const SymEngine::Basic>` to another
+       * `SymEngine::RCP<const SymEngine::Basic>`) as an argument.
+       *
+       * @note The replacement value (the entry in the @p substitution_values
+       * that is paired with a key) need not necessarily be numerical, but may
+       * also be another symbolic type.
+       *
+       * @note With dictionary substitution, partial substitution is allowed
+       * (i.e. an incomplete substitution map can be used and the return type
+       * can be symbolic).
+       */
+      Expression
+      substitute(const SymEngine::map_basic_basic &substitution_values) const;
+
       /**
        * Perform substitution of all symbols found in this object's @p expression
        * that match the @p symbol. Each @p symbol will be substituted with
@@ -746,6 +765,22 @@ namespace Differentiation
       substitute_and_evaluate(
         const types::substitution_map &substitution_values) const;
 
+      /**
+       * Full substitution and evaluation. This creates a Expression by
+       * symbol substitution and then immediately computes its numerical value.
+       *
+       * This function is like the one above, but takes in a SymEngine map
+       * (one that maps a `SymEngine::RCP<const SymEngine::Basic>` to another
+       * `SymEngine::RCP<const SymEngine::Basic>`) as an argument.
+       *
+       * @note All symbols must be resolved by the substitution map in order
+       * for this function to return successfully.
+       */
+      template <typename ReturnType>
+      ReturnType
+      substitute_and_evaluate(
+        const SymEngine::map_basic_basic &substitution_values) const;
+
       //@}
 
       /**
@@ -1256,6 +1291,15 @@ namespace Differentiation
     }
 
 
+    template <typename ReturnType>
+    ReturnType
+    Expression::substitute_and_evaluate(
+      const SymEngine::map_basic_basic &substitution_values) const
+    {
+      return static_cast<ReturnType>(substitute(substitution_values));
+    }
+
+
     template <typename NumberType>
     Expression &
     Expression::operator=(const NumberType &rhs)
diff --git a/include/deal.II/differentiation/sd/symengine_number_visitor_internal.h b/include/deal.II/differentiation/sd/symengine_number_visitor_internal.h
new file mode 100644 (file)
index 0000000..498477c
--- /dev/null
@@ -0,0 +1,1072 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2020 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_number_visitor_internal_h
+#define dealii_differentiation_sd_symengine_number_visitor_internal_h
+
+#include <deal.II/base/config.h>
+
+#ifdef DEAL_II_WITH_SYMENGINE
+
+DEAL_II_DISABLE_EXTRA_DIAGNOSTICS
+// Low level
+#  include <symengine/basic.h>
+#  include <symengine/dict.h>
+#  include <symengine/symengine_exception.h>
+#  include <symengine/symengine_rcp.h>
+
+// Visitor
+#  include <symengine/visitor.h>
+
+DEAL_II_ENABLE_EXTRA_DIAGNOSTICS
+
+#  include <deal.II/base/exceptions.h>
+#  include <deal.II/base/numbers.h>
+
+#  include <deal.II/differentiation/sd/symengine_number_types.h>
+#  include <deal.II/differentiation/sd/symengine_utilities.h>
+
+#  include <boost/serialization/split_member.hpp>
+
+
+DEAL_II_NAMESPACE_OPEN
+
+
+namespace Differentiation
+{
+  namespace SD
+  {
+    namespace internal
+    {
+      /**
+       * A class that implements common subexpression elimination
+       * for dictionary visitor classes.
+       *
+       * It is intended that this class only be used in conjunction
+       * with the DictionarySubstitutionVisitor.
+       *
+       * @author Jean-Paul Pelteret, Isuru Fernando, 2017, 2020
+       */
+      template <typename ReturnType, typename ExpressionType>
+      class CSEDictionaryVisitor
+      {
+        using symbol_vector_pair =
+          std::vector<std::pair<SD::Expression, SD::Expression>>;
+
+      public:
+        /*
+         * Constructor.
+         */
+        CSEDictionaryVisitor() = default;
+
+        /*
+         * Destructor.
+         */
+        virtual ~CSEDictionaryVisitor() = default;
+
+        /**
+         * Initialize and perform common subexpression elimination.
+         *
+         * Here we build the reduced expressions for the @p dependent_functions
+         * as well as a list of intermediate, repeating symbolic expressions
+         * that are extracted @p dependent_functions. This operation leads to
+         * the elimination of repeated expressions, so they only have to be
+         * evaluated once.
+         *
+         * @param dependent_functions A vector of expressions that represent
+         * the dependent functions that CSE will be performed on.
+         */
+        void
+        init(const types::symbol_vector &dependent_functions);
+
+        /**
+         * Evaluate the (reduced) dependent functions using the common
+         * subexpressions generated during the call to init().
+         *
+         * The @p output_values are the numerical result of substituting
+         * each of the @p substitution_values for their corresponding entry
+         * of the @p independent_symbols. Specifically, we first calculate
+         * the numerical values of each reduced subexpression, and feed those
+         * into a reduced equivalent of the dependent variables. The result of
+         * this is then written into the @p output_values array.
+         *
+         * @param[out] output_values A pointer to the first element in an array
+         * of type @p ReturnType. After this call, the underlying array will
+         * hold the numerical result of substituting the @p substitution_values
+         * into the pre-registered dependent functions.
+         * @param[in]  independent_symbols A vector of symbols that represent
+         * the independent variables that are arguments to the previously
+         * defined dependent functions.
+         * @param[in]  substitution_values A pointer to the first element in an
+         * array of type @p ReturnType. Each entry in this array stores the
+         * numerical value that an independent variable is to take for the
+         * purpose of value substitution.
+         *
+         * @note It is expected that both the @p output_values and
+         * @p substitution_values arrays be correctly dimensioned, as there is
+         * no range checking performed on these data structures.
+         * The @p output_values array should have the same number of elements as
+         * the dependent functions first passed to this class in the call to
+         * init(). Similarly, the @p substitution_values array should have the
+         * same number of elements as the @p independent_symbols vector has.
+         *
+         * @note It is expected that there be a 1-1 correspondence between the
+         * entries in @p independent_symbols and @p substitution_values. This
+         * is not checked within this function, so it is up to the user to
+         * ensure that the relationships between the independent variables and
+         * their numerical value be correctly set up and maintained.
+         */
+        void
+        call(ReturnType *                output_values,
+             const types::symbol_vector &independent_symbols,
+             const ReturnType *          substitution_values);
+
+        /**
+         * Write the data of this object to a stream for the purpose
+         * of serialization.
+         */
+        template <class Archive>
+        void
+        save(Archive &archive, const unsigned int version) const;
+
+        /**
+         * Read the data for this object from a stream for the purpose
+         * of serialization.
+         */
+        template <class Archive>
+        void
+        load(Archive &archive, const unsigned int version);
+
+#  ifdef DOXYGEN
+        /**
+         * Write and read the data of this object from a stream for the purpose
+         * of serialization.
+         */
+        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
+
+        /**
+         * Print all of the intermediate reduced expressions, as well as
+         * the final reduced expressions for dependent variables to the
+         * @p stream.
+         */
+        template <typename StreamType>
+        void
+        print(StreamType &stream) const;
+
+        /**
+         * Return a flag stating whether we've performed CSE or not.
+         */
+        bool
+        executed() const;
+
+        /**
+         * The number of intermediate expressions that must be evaluated as
+         * part of the collection of common subexpressions.
+         */
+        unsigned int
+        n_intermediate_expressions() const;
+
+        /**
+         * The size of the final set of reduced expressions.
+         */
+        unsigned int
+        n_reduced_expressions() const;
+
+      protected:
+        /**
+         * Initialize and perform common subexpression elimination.
+         *
+         * This function performs the same action as the other init() function,
+         * except that it works with native SymEngine data types;
+         * the single argument is  a vector of
+         * `SymEngine::RCP<const SymEngine::Basic>`.
+         *
+         * @param dependent_functions A vector of expressions that represent
+         * the dependent functions that CSE will be performed on.
+         */
+        void
+        init(const SymEngine::vec_basic &dependent_functions);
+
+        /**
+         * Evaluate the (reduced) dependent functions using the common
+         * subexpressions generated during the call to init().
+         *
+         * This function performs the same action as the other init() function,
+         * except that it works with native SymEngine data types.
+         * The @p independent_symbols is represented by a vector of
+         * `SymEngine::RCP<const SymEngine::Basic>`.
+         *
+         * @param[out] output_values A pointer to the first element in an array
+         * of type @p ReturnType. After this call, the underlying array will
+         * hold the numerical result of substituting the @p substitution_values
+         * into the pre-registered dependent functions.
+         * @param[in]  independent_symbols A vector of symbols that represent
+         * the independent variables that are arguments to the previously
+         * defined dependent functions.
+         * @param[in]  substitution_values A pointer to the first element in an
+         * array of type @p ReturnType. Each entry in this array stores the
+         * numerical value that an independent variable is to take for the
+         * purpose of value substitution.
+         *
+         * @note The caveats described in the documentation of the other call()
+         * function apply here as well.
+         */
+        void
+        call(ReturnType *                output_values,
+             const SymEngine::vec_basic &independent_symbols,
+             const ReturnType *          substitution_values);
+
+      private:
+        // Note: It would be more efficient to store this data in native
+        // SymEngine types, as it would prevent some copying of the data
+        // structures. However, this makes serialization more difficult,
+        // so we use our own serializable types instead, and lose a bit
+        // of efficiency.
+
+        /**
+         * Intermediate symbols and their definition.
+         */
+        symbol_vector_pair intermediate_symbols_exprs;
+
+        /**
+         * Final reduced expressions.
+         */
+        types::symbol_vector reduced_exprs;
+      };
+
+
+
+      /**
+       * A class to perform dictionary-based substitution as if it were an
+       * optimizer of the "lambda" or "LLVM" variety.
+       *
+       * This class is only really useful to assist in the easy switching
+       * between different optimizers and, more importantly, for integrating
+       * CSE into a dictionary substitution scheme. It is therefore only
+       * intended to be created and used by a BatchOptimizer.
+       *
+       * @author Jean-Paul Pelteret, Isuru Fernando, 2017, 2020
+       */
+      template <typename ReturnType, typename ExpressionType>
+      class DictionarySubstitutionVisitor
+        : public SymEngine::BaseVisitor<
+            DictionarySubstitutionVisitor<ReturnType, ExpressionType>>
+      {
+      public:
+        /*
+         * Constructor.
+         */
+        DictionarySubstitutionVisitor() = default;
+
+        /*
+         * Destructor.
+         */
+        virtual ~DictionarySubstitutionVisitor() override = default;
+
+        /**
+         * Initialization, and registration of the independent and dependent
+         * variables.
+         *
+         * This variation of the initialization function registers a single
+         * dependent expression. If the @p use_cse is set to `true`, then common
+         * subexpression elimination is also performed at the time of
+         * initialization.
+         *
+         * @param independent_symbols A vector of symbols that represent
+         * the independent variables that are arguments to the
+         * @p dependent_function.
+         * @param dependent_function A single symbolic expression that
+         * represents a dependent variable.
+         * @param use_cse A flag to indicate whether or not to use common
+         * subexpression elimination.
+         *
+         * @note After this function is called, no further registration of
+         * dependent functions can be performed. If it is desired that multiple
+         * dependent expressions be registered, then the other variant of this
+         * function that takes in a vector of dependent expressions should be
+         * used.
+         */
+        void
+        init(const types::symbol_vector &independent_symbols,
+             const Expression &          dependent_function,
+             const bool                  use_cse = false);
+
+        /**
+         * Initialization, and registration of the independent and dependent
+         * variables.
+         *
+         * This function performs the same action as the other init() function,
+         * described above, except that it works with native SymEngine data
+         * types. The @p independent_symbols are represented by a vector of
+         * `SymEngine::RCP<const SymEngine::Basic>`, and the
+         * @p dependent_function is of type
+         * `SymEngine::RCP<const SymEngine::Basic>`.
+         *
+         * @param independent_symbols A vector of symbols that represent
+         * the independent variables that are arguments to the
+         * @p dependent_function.
+         * @param dependent_function A single symbolic expression that
+         * represents a dependent variable.
+         * @param use_cse A flag to indicate whether or not to use common
+         * subexpression elimination.
+         *
+         * @note The caveats described in the documentation of the other init()
+         * function apply here as well.
+         */
+        // The following definition is required due to base class CRTP.
+        void
+        init(const SymEngine::vec_basic &independent_symbols,
+             const SymEngine::Basic &    dependent_function,
+             const bool                  use_cse = false);
+
+        /**
+         * Initialization, and registration of the independent and dependent
+         * variables.
+         *
+         * This variation of the initialization function registers a vector of
+         * dependent expressions. If the @p use_cse is set to `true`, then
+         * common subexpression elimination is also performed at the time of
+         * initialization.
+         *
+         * @param independent_symbols A vector of symbols that represent
+         * the independent variables that are arguments to the
+         * @p dependent_functions.
+         * @param dependent_functions A vector of symbolic expressions that
+         * represent the dependent variables.
+         * @param use_cse A flag to indicate whether or not to use common
+         * subexpression elimination.
+         *
+         * @note After this function is called, no further registration of
+         * dependent functions can be performed.
+         */
+        void
+        init(const types::symbol_vector &independent_symbols,
+             const types::symbol_vector &dependent_functions,
+             const bool                  use_cse = false);
+
+
+        /**
+         * Initialization, and registration of the independent and dependent
+         * variables.
+         *
+         * This function performs the same action as the other init() function,
+         * described above, except that it works with native SymEngine data
+         * types. Both the @p independent_symbols and @p dependent_functions
+         * are represented by a vector of
+         * `SymEngine::RCP<const SymEngine::Basic>`.
+         *
+         * @param independent_symbols A vector of symbols that represent
+         * the independent variables that are arguments to the
+         * @p dependent_functions.
+         * @param dependent_functions A vector of symbolic expressions that
+         * represent the dependent variables.
+         * @param use_cse A flag to indicate whether or not to use common
+         * subexpression elimination.
+         *
+         * @note The caveats described in the documentation of the other init()
+         * function apply here as well.
+         */
+        // The following definition is required due to base class CRTP.
+        void
+        init(const SymEngine::vec_basic &independent_symbols,
+             const SymEngine::vec_basic &dependent_functions,
+             const bool                  use_cse = false);
+
+        /**
+         * Evaluate the dependent functions that were registered at
+         * initializtion time.
+         *
+         * The @p output_values are the numerical result of substituting
+         * each of the @p substitution_values for their corresponding entry
+         * of the pre-registered independent variables.
+         *
+         * @param[out] output_values A pointer to the first element in an array
+         * of type @p ReturnType. After this call, the underlying array will
+         * hold the numerical result of substituting the @p substitution_values
+         * into the pre-registered dependent functions.
+         * @param[in]  substitution_values A pointer to the first element in an
+         * array of type @p ReturnType. Each entry in this array stores the
+         * numerical value that an independent variable is to take for the
+         * purpose of value substitution.
+         *
+         * @note It is expected that both the @p output_values and
+         * @p substitution_values arrays be correctly dimensioned, as there is
+         * no range checking performed on these data structures.
+         * The @p output_values array should have the same number of elements as
+         * the dependent functions first passed to this class in the call to
+         * init(). Similarly, the @p substitution_values array should have the
+         * same number of elements as the there were independent variables
+         * that were registered at the time of initialization.
+         *
+         * @note It is expected that there be a 1-1 correspondence between the
+         * entries in independent variables and @p substitution_values. This
+         * is not checked within this function, so it is up to the user to
+         * ensure that the relationships between the independent variables and
+         * their numerical value be correctly set up and maintained.
+         */
+        void
+        call(ReturnType *output_values, const ReturnType *substitution_values);
+
+        /**
+         * Evaluate the dependent function that were registered at
+         * initializtion time.
+         *
+         * The purpose of this function is the same as the other call()
+         * functions, but
+         *
+         * @param substitution_values A vector that stores the numerical values
+         * that each independent variable that was registered with the class
+         * instance is to take for the purpose of value substitution.
+         * @return ReturnType The numerical value associated with the single
+         * dependent function that is registered with this class.
+         *
+         * @note It is expected that there be a 1-1 correspondence between the
+         * entries in independent variables and @p substitution_values.
+         */
+        // The following definition is required due to base class CRTP.
+        ReturnType
+        call(const std::vector<ReturnType> &substitution_values);
+
+        /**
+         * Write the data of this object to a stream for the purpose
+         * of serialization.
+         */
+        template <class Archive>
+        void
+        save(Archive &archive, const unsigned int version) const;
+
+        /**
+         * Read the data for this object from a stream for the purpose
+         * of serialization.
+         */
+        template <class Archive>
+        void
+        load(Archive &archive, const unsigned int version);
+
+#  ifdef DOXYGEN
+        /**
+         * Write and read the data of this object from a stream for the purpose
+         * of serialization.
+         */
+        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
+
+        /**
+         * Print some information on state of the internal data
+         * structures stored in the class.
+         *
+         * @tparam Stream The type for the output stream.
+         * @param stream The output stream to print to.
+         * @param print_independent_symbols A flag to indicate if the independent
+         * variables should be outputted to the @p stream.
+         * @param print_dependent_functions A flag to indicate if the dependent
+         * expressions should be outputted to the @p stream.
+         * @param print_cse_reductions A flag to indicate whether or not all
+         * common subexpressions should be printed to the @p stream.
+         */
+        template <typename StreamType>
+        void
+        print(StreamType &stream,
+              const bool  print_independent_symbols = false,
+              const bool  print_dependent_functions = false,
+              const bool  print_cse_reductions      = false) const;
+
+#  ifndef DOXYGEN
+        // The following definitions are required due to base class CRTP.
+        // Since these are not used, and therefore not important to
+        // understand, we'll define them in the most concise manner possible.
+        // We also won't bother to document their existence, since they cannot
+        // be used.
+#    define IMPLEMENT_DSV_BVISIT(Argument)       \
+      void bvisit(const Argument &)              \
+      {                                          \
+        AssertThrow(false, ExcNotImplemented()); \
+      }
+
+        IMPLEMENT_DSV_BVISIT(SymEngine::Basic)
+        IMPLEMENT_DSV_BVISIT(SymEngine::Symbol)
+        IMPLEMENT_DSV_BVISIT(SymEngine::Constant)
+        IMPLEMENT_DSV_BVISIT(SymEngine::Integer)
+        IMPLEMENT_DSV_BVISIT(SymEngine::Rational)
+        IMPLEMENT_DSV_BVISIT(SymEngine::RealDouble)
+        IMPLEMENT_DSV_BVISIT(SymEngine::ComplexDouble)
+        IMPLEMENT_DSV_BVISIT(SymEngine::Add)
+        IMPLEMENT_DSV_BVISIT(SymEngine::Mul)
+        IMPLEMENT_DSV_BVISIT(SymEngine::Pow)
+        IMPLEMENT_DSV_BVISIT(SymEngine::Log)
+        IMPLEMENT_DSV_BVISIT(SymEngine::Sin)
+        IMPLEMENT_DSV_BVISIT(SymEngine::Cos)
+        IMPLEMENT_DSV_BVISIT(SymEngine::Tan)
+        IMPLEMENT_DSV_BVISIT(SymEngine::Csc)
+        IMPLEMENT_DSV_BVISIT(SymEngine::Sec)
+        IMPLEMENT_DSV_BVISIT(SymEngine::Cot)
+        IMPLEMENT_DSV_BVISIT(SymEngine::ASin)
+        IMPLEMENT_DSV_BVISIT(SymEngine::ACos)
+        IMPLEMENT_DSV_BVISIT(SymEngine::ATan)
+        IMPLEMENT_DSV_BVISIT(SymEngine::ATan2)
+        IMPLEMENT_DSV_BVISIT(SymEngine::ACsc)
+        IMPLEMENT_DSV_BVISIT(SymEngine::ASec)
+        IMPLEMENT_DSV_BVISIT(SymEngine::ACot)
+        IMPLEMENT_DSV_BVISIT(SymEngine::Sinh)
+        IMPLEMENT_DSV_BVISIT(SymEngine::Cosh)
+        IMPLEMENT_DSV_BVISIT(SymEngine::Tanh)
+        IMPLEMENT_DSV_BVISIT(SymEngine::Csch)
+        IMPLEMENT_DSV_BVISIT(SymEngine::Sech)
+        IMPLEMENT_DSV_BVISIT(SymEngine::Coth)
+        IMPLEMENT_DSV_BVISIT(SymEngine::ASinh)
+        IMPLEMENT_DSV_BVISIT(SymEngine::ACosh)
+        IMPLEMENT_DSV_BVISIT(SymEngine::ATanh)
+        IMPLEMENT_DSV_BVISIT(SymEngine::ACsch)
+        IMPLEMENT_DSV_BVISIT(SymEngine::ACoth)
+        IMPLEMENT_DSV_BVISIT(SymEngine::ASech)
+        IMPLEMENT_DSV_BVISIT(SymEngine::Abs)
+        IMPLEMENT_DSV_BVISIT(SymEngine::Gamma)
+        IMPLEMENT_DSV_BVISIT(SymEngine::LogGamma)
+        IMPLEMENT_DSV_BVISIT(SymEngine::Erf)
+        IMPLEMENT_DSV_BVISIT(SymEngine::Erfc)
+        IMPLEMENT_DSV_BVISIT(SymEngine::Max)
+        IMPLEMENT_DSV_BVISIT(SymEngine::Min)
+
+#    undef IMPLEMENT_DSV_BVISIT
+#  endif // DOXYGEN
+
+      private:
+        // Note: It would be more efficient to store this data in native
+        // SymEngine types, as it would prevent some copying of the data
+        // structures. However, this makes serialization more difficult,
+        // so we use our own serializable types instead, and lose a bit
+        // of efficiency.
+
+        /**
+         * A vector of symbols that represent the independent variables.
+         */
+        SD::types::symbol_vector independent_symbols;
+
+        /**
+         * A vector of expressions that represent dependent functions.
+         */
+        SD::types::symbol_vector dependent_functions;
+
+        /**
+         * A data structure that may be used to invoke common subexpression
+         * elimination on the dependent functions, with the aim to decrease
+         * the time taken to evaluate them.
+         */
+        CSEDictionaryVisitor<ReturnType, ExpressionType> cse;
+      };
+
+
+
+      /* ------------------ inline and template functions ------------------ */
+
+
+#  ifndef DOXYGEN
+
+      /* -------------- CommonSubexpressionEliminationVisitor -------------- */
+
+
+      template <typename ReturnType, typename ExpressionType>
+      void
+      CSEDictionaryVisitor<ReturnType, ExpressionType>::init(
+        const SD::types::symbol_vector &dependent_functions)
+      {
+        init(Utilities::convert_expression_vector_to_basic_vector(
+          dependent_functions));
+      }
+
+
+
+      template <typename ReturnType, typename ExpressionType>
+      void
+      CSEDictionaryVisitor<ReturnType, ExpressionType>::init(
+        const SymEngine::vec_basic &dependent_functions)
+      {
+        // After the next call, the data stored in replacements is structured
+        // as follows:
+        //
+        // replacements[i] := [f, f(x)]
+        // replacements[i].first  = intermediate function label "f"
+        // replacements[i].second = intermediate function definition "f(x)"
+        //
+        // It is to be evaluated top down (i.e. index 0 to
+        // replacements.size()), with the results going back into the
+        // substitution map for the next levels. So for each "i", "x" are the
+        // superset of the input values and the previously evaluated [f_0(x),
+        // f_1(x), ..., f_{i-1}(x)].
+        //
+        // The final result is a set of reduced expressions
+        // that must be computed after the replacement
+        // values have been computed.
+        SymEngine::vec_pair  se_replacements;
+        SymEngine::vec_basic se_reduced_exprs;
+        SymEngine::cse(se_replacements, se_reduced_exprs, dependent_functions);
+
+        intermediate_symbols_exprs =
+          Utilities::convert_basic_pair_vector_to_expression_pair_vector(
+            se_replacements);
+        reduced_exprs = Utilities::convert_basic_vector_to_expression_vector(
+          se_reduced_exprs);
+      }
+
+
+
+      template <typename ReturnType, typename ExpressionType>
+      void
+      CSEDictionaryVisitor<ReturnType, ExpressionType>::call(
+        ReturnType *                    output_values,
+        const SD::types::symbol_vector &independent_symbols,
+        const ReturnType *              substitution_values)
+      {
+        call(output_values,
+             Utilities::convert_expression_vector_to_basic_vector(
+               independent_symbols),
+             substitution_values);
+      }
+
+
+
+      template <typename ReturnType, typename ExpressionType>
+      void
+      CSEDictionaryVisitor<ReturnType, ExpressionType>::call(
+        ReturnType *                output_values,
+        const SymEngine::vec_basic &independent_symbols,
+        const ReturnType *          substitution_values)
+      {
+        Assert(n_reduced_expressions() > 0, ExcInternalError());
+
+        // First we add the input values into the substitution map...
+        SymEngine::map_basic_basic substitution_value_map;
+        for (unsigned i = 0; i < independent_symbols.size(); ++i)
+          substitution_value_map[independent_symbols[i]] =
+            static_cast<const SymEngine::RCP<const SymEngine::Basic> &>(
+              ExpressionType(substitution_values[i]));
+
+        // ... followed by any intermediate evaluations due to the application
+        // of CSE. These are fed directly back into the substitution map...
+        for (unsigned i = 0; i < intermediate_symbols_exprs.size(); ++i)
+          {
+            const SymEngine::RCP<const SymEngine::Basic> &cse_symbol =
+              intermediate_symbols_exprs[i].first;
+            const SymEngine::RCP<const SymEngine::Basic> &cse_expr =
+              intermediate_symbols_exprs[i].second;
+            Assert(substitution_value_map.find(cse_symbol) ==
+                     substitution_value_map.end(),
+                   ExcMessage(
+                     "Reduced symbol already appears in substitution map. "
+                     "Is there a clash between the reduced symbol name and "
+                     "the symbol used for an independent variable?"));
+            substitution_value_map[cse_symbol] =
+              static_cast<const SymEngine::RCP<const SymEngine::Basic> &>(
+                ExpressionType(ExpressionType(cse_expr)
+                                 .template substitute_and_evaluate<ReturnType>(
+                                   substitution_value_map)));
+          }
+
+        // ... followed by the final reduced expressions
+        for (unsigned i = 0; i < reduced_exprs.size(); ++i)
+          output_values[i] = ExpressionType(reduced_exprs[i])
+                               .template substitute_and_evaluate<ReturnType>(
+                                 substitution_value_map);
+      }
+
+
+
+      template <typename ReturnType, typename ExpressionType>
+      template <class Archive>
+      void
+      CSEDictionaryVisitor<ReturnType, ExpressionType>::save(
+        Archive &ar,
+        const unsigned int /*version*/) const
+      {
+        // The reduced expressions depend on the intermediate expressions,
+        // so we serialize the latter before the former.
+        ar &intermediate_symbols_exprs;
+        ar &reduced_exprs;
+      }
+
+
+
+      template <typename ReturnType, typename ExpressionType>
+      template <class Archive>
+      void
+      CSEDictionaryVisitor<ReturnType, ExpressionType>::load(
+        Archive &ar,
+        const unsigned int /*version*/)
+      {
+        Assert(intermediate_symbols_exprs.empty(), ExcInternalError());
+        Assert(reduced_exprs.empty(), ExcInternalError());
+
+        // The reduced expressions depend on the intermediate expressions,
+        // so we deserialize the latter before the former.
+        ar &intermediate_symbols_exprs;
+        ar &reduced_exprs;
+      }
+
+
+
+      template <typename ReturnType, typename ExpressionType>
+      template <typename StreamType>
+      void
+      CSEDictionaryVisitor<ReturnType, ExpressionType>::print(
+        StreamType &stream) const
+      {
+        stream << "Common subexpression elimination: \n";
+        stream << "  Intermediate reduced expressions: \n";
+        for (unsigned i = 0; i < intermediate_symbols_exprs.size(); ++i)
+          {
+            const SymEngine::RCP<const SymEngine::Basic> &cse_symbol =
+              intermediate_symbols_exprs[i].first;
+            const SymEngine::RCP<const SymEngine::Basic> &cse_expr =
+              intermediate_symbols_exprs[i].second;
+            stream << "  " << i << ": " << cse_symbol << " = " << cse_expr
+                   << "\n";
+          }
+
+        stream << "  Final reduced expressions for dependent variables: \n";
+        for (unsigned i = 0; i < reduced_exprs.size(); ++i)
+          stream << "  " << i << ": " << reduced_exprs[i] << "\n";
+
+        stream << std::flush;
+      }
+
+
+
+      template <typename ReturnType, typename ExpressionType>
+      bool
+      CSEDictionaryVisitor<ReturnType, ExpressionType>::executed() const
+      {
+        // For dictionary substitution, the CSE algorithm moves
+        // ownership of the dependent function expression definition
+        // to the entries in reduced_exprs. So its size thus determines
+        // whether CSE has been executed or not.
+        return (n_reduced_expressions() > 0) ||
+               (n_intermediate_expressions() > 0);
+      }
+
+
+
+      template <typename ReturnType, typename ExpressionType>
+      unsigned int
+      CSEDictionaryVisitor<ReturnType,
+                           ExpressionType>::n_intermediate_expressions() const
+      {
+        return intermediate_symbols_exprs.size();
+      }
+
+
+
+      template <typename ReturnType, typename ExpressionType>
+      unsigned int
+      CSEDictionaryVisitor<ReturnType, ExpressionType>::n_reduced_expressions()
+        const
+      {
+        return reduced_exprs.size();
+      }
+
+
+
+      /* ------------------ DictionarySubstitutionVisitor ------------------ */
+
+
+      // Perform some extra checks in debug mode. This will give some verbose
+      // output if something goes wrong.
+      // #define DEBUG_EXTRA_CSE
+
+
+      template <typename ReturnType, typename ExpressionType>
+      void
+      DictionarySubstitutionVisitor<ReturnType, ExpressionType>::init(
+        const types::symbol_vector &inputs,
+        const SD::Expression &      output,
+        const bool                  use_cse)
+      {
+        init(inputs, types::symbol_vector{output}, use_cse);
+      }
+
+
+
+      template <typename ReturnType, typename ExpressionType>
+      void
+      DictionarySubstitutionVisitor<ReturnType, ExpressionType>::init(
+        const SymEngine::vec_basic &inputs,
+        const SymEngine::Basic &    output,
+        const bool                  use_cse)
+      {
+        init(Utilities::convert_basic_vector_to_expression_vector(inputs),
+             SD::Expression(output.rcp_from_this()),
+             use_cse);
+      }
+
+
+
+      template <typename ReturnType, typename ExpressionType>
+      void
+      DictionarySubstitutionVisitor<ReturnType, ExpressionType>::init(
+        const SymEngine::vec_basic &inputs,
+        const SymEngine::vec_basic &outputs,
+        const bool                  use_cse)
+      {
+        init(Utilities::convert_basic_vector_to_expression_vector(inputs),
+             Utilities::convert_basic_vector_to_expression_vector(outputs),
+             use_cse);
+      }
+
+
+
+      template <typename ReturnType, typename ExpressionType>
+      void
+      DictionarySubstitutionVisitor<ReturnType, ExpressionType>::init(
+        const types::symbol_vector &inputs,
+        const types::symbol_vector &outputs,
+        const bool                  use_cse)
+      {
+        independent_symbols.clear();
+        dependent_functions.clear();
+
+        independent_symbols = inputs;
+
+        // Perform common subexpression elimination if requested
+        // Note: After this is done, the results produced by
+        // dependent_functions and cse.reduced_exprs should be
+        // the same. We could keep the former so that we can print
+        // out the original expressions if we wish to do so.
+        if (use_cse == false)
+          dependent_functions = outputs;
+        else
+          {
+            cse.init(outputs);
+
+#    if defined(DEBUG) && defined(DEBUG_EXTRA_CSE)
+            // Check the result of CSE by directly evaluating
+            // the dependent variables on an individual basis
+            dependent_variables = outputs;
+#    endif
+          }
+      }
+
+
+
+      template <typename ReturnType, typename ExpressionType>
+      ReturnType
+      DictionarySubstitutionVisitor<ReturnType, ExpressionType>::call(
+        const std::vector<ReturnType> &substitution_values)
+      {
+        Assert(
+          dependent_functions.size() == 1,
+          ExcMessage(
+            "Cannot use this call function when more than one symbolic expression is to be evaluated."));
+        Assert(
+          substitution_values.size() == independent_symbols.size(),
+          ExcMessage(
+            "Input substitution vector does not match size of symbol vector."));
+
+        ReturnType out = dealii::internal::NumberType<ReturnType>::value(0.0);
+        call(&out, substitution_values.data());
+        return out;
+      }
+
+
+
+      template <typename ReturnType, typename ExpressionType>
+      void
+      DictionarySubstitutionVisitor<ReturnType, ExpressionType>::call(
+        ReturnType *      output_values,
+        const ReturnType *substitution_values)
+      {
+        // Check to see if CSE has been performed
+        if (cse.executed())
+          {
+            cse.call(output_values, independent_symbols, substitution_values);
+
+#    if defined(DEBUG) && defined(DEBUG_EXTRA_CSE)
+            // Directly check the result of the CSE calculations
+            Assert(dependent_variables.size() > 0, ExcInternalError());
+            SymEngine::map_basic_basic substitution_value_map;
+            for (unsigned i = 0; i < independent_variables.size(); ++i)
+              substitution_value_map[independent_variables[i]] =
+                static_cast<const SymEngine::RCP<const SymEngine::Basic> &>(
+                  ExpressionType(substitution_values[i]));
+            const ReturnType tol =
+              dealii::internal::NumberType<ReturnType>::value(1e-9);
+
+            for (unsigned i = 0; i < dependent_variables.size(); ++i)
+              {
+                const ReturnType result =
+                  ExpressionType(dependent_variables[i])
+                    .template substitute_and_evaluate<ReturnType>(
+                      substitution_value_map);
+
+                // Note: Need to do an abs() on the tolerance to cater for
+                // the complex value case.
+                const bool value_correct =
+                  (std::abs(output_values[i] - result) / std::abs(result) <
+                   std::abs(tol));
+                if (value_correct == false)
+                  {
+                    std::ostringstream stream;
+                    stream
+                      << "A CSE based evaluation is incorrect. The result should have been "
+                      << result << " but was rather calculated to be "
+                      << output_values[i] << ".\n Substitution map: \n";
+                    for (typename SymEngine::map_basic_basic::const_iterator
+                           it = substitution_value_map.begin();
+                         it != substitution_value_map.end();
+                         ++it)
+                      {
+                        stream << *(it->first) << " = " << *(it->second)
+                               << "\n";
+                      }
+                    stream << std::endl;
+                    const std::string msg = stream.str();
+                    Assert(value_correct, ExcMessage(msg));
+                  }
+              }
+#    endif
+          }
+        else
+          {
+            // Build a substitution map.
+            SymEngine::map_basic_basic substitution_value_map;
+            for (unsigned i = 0; i < independent_symbols.size(); ++i)
+              substitution_value_map[independent_symbols[i]] =
+                static_cast<const SymEngine::RCP<const SymEngine::Basic> &>(
+                  ExpressionType(substitution_values[i]));
+
+            // Since we don't know how to definitively evaluate the
+            // input number type, we create a generic Expression
+            // with the given symbolic expression and ask it to perform
+            // substitution and evaluation for us.
+            Assert(dependent_functions.size() > 0, ExcInternalError());
+            for (unsigned i = 0; i < dependent_functions.size(); ++i)
+              output_values[i] =
+                ExpressionType(dependent_functions[i])
+                  .template substitute_and_evaluate<ReturnType>(
+                    substitution_value_map);
+          }
+      }
+
+
+
+      template <typename ReturnType, typename ExpressionType>
+      template <class Archive>
+      void
+      DictionarySubstitutionVisitor<ReturnType, ExpressionType>::save(
+        Archive &          ar,
+        const unsigned int version) const
+      {
+        // Add some dynamic information to determine if CSE has been used,
+        // without relying on the CSE class when deserializing.
+        // const bool used_cse = cse.executed();
+        // ar &used_cse;
+
+        // CSE and dependent variables both require the independent
+        // symbols, so we serialize them first. The dependent variables
+        // might depend on the outcome of CSE, so we have to serialize
+        // them last.
+        ar &independent_symbols;
+        cse.save(ar, version);
+        ar &dependent_functions;
+      }
+
+
+
+      template <typename ReturnType, typename ExpressionType>
+      template <class Archive>
+      void
+      DictionarySubstitutionVisitor<ReturnType, ExpressionType>::load(
+        Archive &          ar,
+        const unsigned int version)
+      {
+#    if !(defined(DEBUG) && defined(DEBUG_EXTRA_CSE))
+        Assert(dependent_variables.empty(), ExcInternalError());
+#    endif
+        Assert(cse.executed() == false, ExcInternalError());
+        Assert(cse.n_intermediate_expressions() == 0, ExcInternalError());
+        Assert(cse.n_reduced_expressions() == 0, ExcInternalError());
+
+        // CSE and dependent variables both require the independent
+        // symbols, so we deserialize them first. The dependent variables
+        // might depend on the outcome of CSE, so we have to deserialize
+        // them last.
+        ar &independent_symbols;
+        cse.load(ar, version);
+        ar &dependent_functions;
+      }
+
+
+
+      template <typename ReturnType, typename ExpressionType>
+      template <typename StreamType>
+      void
+      DictionarySubstitutionVisitor<ReturnType, ExpressionType>::print(
+        StreamType &stream,
+        const bool  print_independent_symbols,
+        const bool  print_dependent_functions,
+        const bool  print_cse_reductions) const
+      {
+        if (print_independent_symbols)
+          {
+            stream << "Independent variables: \n";
+            for (unsigned i = 0; i < independent_symbols.size(); ++i)
+              stream << "  " << i << ": " << independent_symbols[i] << "\n";
+
+            stream << std::flush;
+          }
+
+        // Check to see if CSE has been performed
+        if (print_cse_reductions && cse.executed())
+          {
+            cse.print(stream);
+          }
+        else
+          {
+            Assert(dependent_functions.size() > 0, ExcInternalError());
+
+            if (print_dependent_functions)
+              {
+                stream << "Dependent variables: \n";
+                for (unsigned i = 0; i < dependent_functions.size(); ++i)
+                  stream << "  " << i << dependent_functions[i] << "\n";
+
+                stream << std::flush;
+              }
+          }
+      }
+
+#  endif // DOXYGEN
+
+    } // namespace internal
+  }   // namespace SD
+} // namespace Differentiation
+
+
+DEAL_II_NAMESPACE_CLOSE
+
+#endif // DEAL_II_WITH_SYMENGINE
+
+#endif // dealii_differentiation_sd_symengine_number_visitor_internal_h
index 8f8bf1d714daa25081e76abfffb3fcff177ea78a..129f3da777208080644961ba041daaadf602c0ae 100644 (file)
@@ -18,6 +18,7 @@ INCLUDE_DIRECTORIES(BEFORE ${CMAKE_CURRENT_BINARY_DIR})
 SET(_src
   symengine_math.cc
   symengine_number_types.cc
+  symengine_number_visitor_internal.cc
   symengine_scalar_operations.cc
   symengine_tensor_operations.cc
   symengine_types.cc
@@ -25,6 +26,7 @@ SET(_src
   )
 
 SET(_inst
+  symengine_number_visitor_internal.inst.in
   symengine_tensor_operations.inst.in
   )
 
index 051771555e568cab7580d4ceb9c5f502f77ee48f..20806592a169462c5e21c725ba5d8aec36856bbb 100644 (file)
@@ -282,12 +282,20 @@ namespace Differentiation
     /* ------------- Dictionary-based substitution ------------------------- */
 
 
+    Expression
+    Expression::substitute(
+      const SymEngine::map_basic_basic &substitution_values) const
+    {
+      return Expression(get_expression().subs(substitution_values));
+    }
+
+
     Expression
     Expression::substitute(
       const types::substitution_map &substitution_values) const
     {
-      return Expression(get_expression().subs(
-        Utilities::convert_expression_map_to_basic_map(substitution_values)));
+      return substitute(
+        Utilities::convert_expression_map_to_basic_map(substitution_values));
     }
 
 
diff --git a/source/differentiation/sd/symengine_number_visitor_internal.cc b/source/differentiation/sd/symengine_number_visitor_internal.cc
new file mode 100644 (file)
index 0000000..9667ad2
--- /dev/null
@@ -0,0 +1,30 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2020 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.
+//
+// ---------------------------------------------------------------------
+
+#include <deal.II/base/config.h>
+
+#ifdef DEAL_II_WITH_SYMENGINE
+
+#  include <deal.II/differentiation/sd/symengine_number_visitor_internal.h>
+
+DEAL_II_NAMESPACE_OPEN
+
+/* --- Explicit instantiations --- */
+#  include "symengine_number_visitor_internal.inst"
+
+
+DEAL_II_NAMESPACE_CLOSE
+
+#endif // DEAL_II_WITH_SYMENGINE
diff --git a/source/differentiation/sd/symengine_number_visitor_internal.inst.in b/source/differentiation/sd/symengine_number_visitor_internal.inst.in
new file mode 100644 (file)
index 0000000..b951e5e
--- /dev/null
@@ -0,0 +1,49 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2020 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.
+//
+// ---------------------------------------------------------------------
+
+
+for (number : REAL_SCALARS)
+  {
+    namespace Differentiation
+    \{
+      namespace SD
+      \{
+        namespace internal
+        \{
+
+          template class CSEDictionaryVisitor<number, SD::Expression>;
+          template class DictionarySubstitutionVisitor<number, SD::Expression>;
+
+        \}
+      \}
+    \}
+  }
+
+for (number : COMPLEX_SCALARS)
+  {
+    namespace Differentiation
+    \{
+      namespace SD
+      \{
+        namespace internal
+        \{
+
+          template class CSEDictionaryVisitor<number, SD::Expression>;
+          template class DictionarySubstitutionVisitor<number, SD::Expression>;
+
+        \}
+      \}
+    \}
+  }

In the beginning the Universe was created. This has made a lot of people very angry and has been widely regarded as a bad move.

Douglas Adams


Typeset in Trocchi and Trocchi Bold Sans Serif.