symengine_number_types.h

Go to the documentation of this file.

// ------------------------------------------------------------------------
//
// SPDX-License-Identifier: LGPL-2.1-or-later
// Copyright (C) 2019 - 2023 by the deal.II authors
//
// This file is part of the deal.II library.
//
// Part of the source code is dual licensed under Apache-2.0 WITH
// LLVM-exception OR LGPL-2.1-or-later. Detailed license information
// governing the source code and code contributions can be found in
// LICENSE.md and CONTRIBUTING.md at the top level directory of deal.II.
//
// ------------------------------------------------------------------------
 
#ifndef dealii_differentiation_sd_symengine_number_types_h
#define dealii_differentiation_sd_symengine_number_types_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>
 
// Number types
#  include <symengine/expression.h>
#  include <symengine/integer.h>
#  include <symengine/logic.h>
#  include <symengine/number.h>
#  include <symengine/rational.h>
 
// Number operations
#  include <symengine/add.h>
#  include <symengine/functions.h>
#  include <symengine/mul.h>
#  include <symengine/pow.h>
 
// Evaluation
#  include <symengine/eval.h>
#  include <symengine/eval_arb.h>
#  include <symengine/eval_double.h>
 
// Differentiation
#  include <deal.II/base/exceptions.h>
#  include <deal.II/base/numbers.h>
 
#  include <deal.II/differentiation/sd/symengine_number_traits.h>
#  include <deal.II/differentiation/sd/symengine_types.h>
 
#  include <boost/serialization/split_member.hpp>
 
#  include <symengine/derivative.h>
 
#  include <algorithm>
#  include <memory>
#  include <sstream>
#  include <type_traits>
#  include <utility>
#  include <vector>
 
DEAL_II_NAMESPACE_OPEN
 
namespace Differentiation
{
  namespace SD
  {
    DeclException1(
      ExcSymEngineParserError,
      std::string,
      << "The string '" << arg1
      << "' could not be parsed successfully. Are you sure that (1) it "
      << "consists of legitimate operations and syntax, and (2) you've "
      << "previously declared all symbolic variables that are present "
      << "in the expression?");
 
    class Expression
    {
    public:
      Expression();
 
      explicit Expression(const bool value);
 
      template <typename NumberType,
                typename = std::enable_if_t<std::is_arithmetic_v<NumberType>>>
      explicit Expression(const NumberType &value);
 
      template <typename NumberType,
                typename = std::enable_if_t<std::is_arithmetic_v<NumberType>>>
      explicit Expression(const std::complex<NumberType> &value);
 
      Expression(const SymEngine::integer_class &value);
 
      template <typename NumberType,
                typename = std::enable_if_t<std::is_integral_v<NumberType>>>
      Expression(const NumberType &numerator, const NumberType &denominator);
 
      Expression(const SymEngine::rational_class &value);
 
      Expression(const Expression &condition,
                 const Expression &expression_if_true,
                 const Expression &expression_if_false);
 
      Expression(const std::vector<std::pair<Expression, Expression>>
                                  &condition_expression,
                 const Expression &expression_otherwise);
 
      Expression(const std::vector<std::pair<Expression, Expression>>
                   &condition_expression);
 
 
      Expression(const char *symbol);
 
      Expression(const std::string &symb_expr,
                 const bool         parse_as_expression = false);
 
      Expression(const std::string          &symbol_func,
                 const types::symbol_vector &arguments);
 
      Expression(const Expression &rhs) = default;
 
      explicit Expression(const SymEngine::Expression &rhs);
 
      Expression(const SymEngine::RCP<const SymEngine::Basic> &rhs);
 
      Expression(Expression &&rhs) = default;
 
      Expression(SymEngine::RCP<const SymEngine::Basic> &&rhs);
 
      virtual ~Expression() = default;
 
      Expression &
      parse(const std::string &expression);
 
      std::ostream &
      print(std::ostream &stream) const;
 
      void
      save(std::ostream &stream) const;
 
      void
      load(std::istream &stream);
 
      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
 
      const SymEngine::Expression &
      get_expression() const;
 
      const SymEngine::Basic &
      get_value() const;
 
      const SymEngine::RCP<const SymEngine::Basic> &
      get_RCP() const;
 
      Expression &
      operator=(const Expression &rhs);
 
      Expression &
      operator=(Expression &&rhs) noexcept;
 
      Expression &
      operator+=(const Expression &rhs);
 
      Expression &
      operator-=(const Expression &rhs);
 
      Expression &
      operator*=(const Expression &rhs);
 
      Expression &
      operator/=(const Expression &rhs);
 
      template <typename NumberType>
      Expression &
      operator=(const NumberType &rhs);
 
      Expression
      operator-() const;
 
      template <typename NumberType>
      Expression &
      operator+=(const NumberType &rhs);
 
      template <typename NumberType>
      Expression &
      operator-=(const NumberType &rhs);
 
      template <typename NumberType>
      Expression &
      operator*=(const NumberType &rhs);
 
      template <typename NumberType>
      Expression &
      operator/=(const NumberType &rhs);
 
      Expression
      differentiate(const Expression &symbol) const;
 
      Expression
      differentiate(
        const SymEngine::RCP<const SymEngine::Symbol> &symbol) const;
 
      Expression
      differentiate(const SymEngine::RCP<const SymEngine::Basic> &symbol) const;
 
      Expression
      substitute(const types::substitution_map &substitution_values) const;
 
      Expression
      substitute(const SymEngine::map_basic_basic &substitution_values) const;
 
      Expression
      substitute(const Expression &symbol, const Expression &value) const;
 
      template <typename NumberType>
      Expression
      substitute(const Expression &symbol, const NumberType &value) const;
 
      template <typename ReturnType>
      ReturnType
      substitute_and_evaluate(
        const types::substitution_map &substitution_values) const;
 
      template <typename ReturnType>
      ReturnType
      substitute_and_evaluate(
        const SymEngine::map_basic_basic &substitution_values) const;
 
      template <typename ResultType>
      explicit operator ResultType() const;
 
      explicit operator const SymEngine::Expression &() const;
 
      operator const SymEngine::RCP<const SymEngine::Basic> &() const;
 
    protected:
      SymEngine::Expression &
      get_expression();
 
    private:
      SymEngine::Expression expression;
    };
 
    template <>
    struct is_symengine_number<Expression> : std::true_type
    {};
 
 
    template <>
    struct is_sd_number<Expression> : std::true_type
    {};
 
    std::ostream &
    operator<<(std::ostream &stream, const Expression &expression);
 
    std::istream &
    operator>>(std::istream &stream, Expression &expression);
 
    Expression
    operator==(const Expression &lhs, const Expression &rhs);
 
    Expression
    operator!=(const Expression &lhs, const Expression &rhs);
 
    Expression
    operator<(const Expression &lhs, const Expression &rhs);
 
    Expression
    operator>(const Expression &lhs, const Expression &rhs);
 
    Expression
    operator<=(const Expression &lhs, const Expression &rhs);
 
    Expression
    operator>=(const Expression &lhs, const Expression &rhs);
 
    Expression
    operator!(const Expression &expression);
 
    Expression
    operator&(const Expression &lhs, const Expression &rhs);
 
    Expression
    operator|(const Expression &lhs, const Expression &rhs);
 
    Expression
    operator^(const Expression &lhs, const Expression &rhs);
 
    Expression
    operator&&(const Expression &lhs, const Expression &rhs);
 
    Expression
    operator||(const Expression &lhs, const Expression &rhs);
 
    Expression
    operator+(Expression lhs, const Expression &rhs);
 
    Expression
    operator-(Expression lhs, const Expression &rhs);
 
    Expression
    operator*(Expression lhs, const Expression &rhs);
 
    Expression
    operator/(Expression lhs, const Expression &rhs);
 
    template <typename NumberType,
              typename = std::enable_if_t<
                std::is_constructible_v<Expression, NumberType>>>
    inline Expression
    operator+(const NumberType &lhs, const Expression &rhs)
    {
      return Expression(lhs) + rhs;
    }
 
    template <typename NumberType,
              typename = std::enable_if_t<
                std::is_constructible_v<Expression, NumberType>>>
    inline Expression
    operator+(const Expression &lhs, const NumberType &rhs)
    {
      return lhs + Expression(rhs);
    }
 
    template <typename NumberType,
              typename = std::enable_if_t<
                std::is_constructible_v<Expression, NumberType>>>
    inline Expression
    operator-(const NumberType &lhs, const Expression &rhs)
    {
      return Expression(lhs) - rhs;
    }
 
    template <typename NumberType,
              typename = std::enable_if_t<
                std::is_constructible_v<Expression, NumberType>>>
    inline Expression
    operator-(const Expression &lhs, const NumberType &rhs)
    {
      return lhs - Expression(rhs);
    }
 
    template <typename NumberType,
              typename = std::enable_if_t<
                std::is_constructible_v<Expression, NumberType>>>
    inline Expression
    operator*(const NumberType &lhs, const Expression &rhs)
    {
      return Expression(lhs) * rhs;
    }
 
    template <typename NumberType,
              typename = std::enable_if_t<
                std::is_constructible_v<Expression, NumberType>>>
    inline Expression
    operator*(const Expression &lhs, const NumberType &rhs)
    {
      return lhs * Expression(rhs);
    }
 
    template <typename NumberType,
              typename = std::enable_if_t<
                std::is_constructible_v<Expression, NumberType>>>
    inline Expression
    operator/(const NumberType &lhs, const Expression &rhs)
    {
      return Expression(lhs) / rhs;
    }
 
    template <typename NumberType,
              typename = std::enable_if_t<
                std::is_constructible_v<Expression, NumberType>>>
    inline Expression
    operator/(const Expression &lhs, const NumberType &rhs)
    {
      return lhs / Expression(rhs);
    }
 
  } // namespace SD
} // namespace Differentiation
 
 
/* ---------------------- inline and template functions -------------------- */
 
 
#  ifndef DOXYGEN
 
 
namespace Differentiation
{
  namespace SD
  {
    template <typename NumberType, typename>
    Expression::Expression(const NumberType &value)
      : expression(value)
    {}
 
 
    template <typename NumberType, typename>
    Expression::Expression(const std::complex<NumberType> &value)
      : expression(value)
    {}
 
 
    template <typename NumberType, typename>
    Expression::Expression(const NumberType &numerator,
                           const NumberType &denominator)
      : expression(
          SymEngine::Rational::from_two_ints(*SymEngine::integer(numerator),
                                             *SymEngine::integer(denominator)))
    {}
 
 
    template <class Archive>
    void
    Expression::save(Archive &ar, const unsigned int /*version*/) const
    {
      std::stringstream sstream;
      sstream << *this;
      const std::string expr = sstream.str();
      ar               &expr;
    }
 
 
    template <class Archive>
    void
    Expression::load(Archive &ar, const unsigned int /*version*/)
    {
      std::string expr;
      ar         &expr;
      parse(expr);
    }
 
 
    template <typename NumberType>
    Expression
    Expression::substitute(const Expression &symbol,
                           const NumberType &value) const
    {
      Assert(SymEngine::is_a<SymEngine::Symbol>(symbol.get_value()),
             ExcMessage(
               "Substitution with a number that does not represent a symbol."));
 
      types::substitution_map sub_vals;
      sub_vals[symbol] = Expression(value);
      return substitute(sub_vals);
    }
 
 
    template <typename ReturnType>
    ReturnType
    Expression::substitute_and_evaluate(
      const types::substitution_map &substitution_values) const
    {
      return static_cast<ReturnType>(substitute(substitution_values));
    }
 
 
    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)
    {
      *this = Expression(rhs);
      return *this;
    }
 
 
    template <typename NumberType>
    Expression &
    Expression::operator+=(const NumberType &rhs)
    {
      *this = Expression(SymEngine::add(get_RCP(), Expression(rhs).get_RCP()));
      return *this;
    }
 
 
    template <typename NumberType>
    Expression &
    Expression::operator-=(const NumberType &rhs)
    {
      *this = Expression(SymEngine::sub(get_RCP(), Expression(rhs).get_RCP()));
      return *this;
    }
 
 
    template <typename NumberType>
    Expression &
    Expression::operator*=(const NumberType &rhs)
    {
      *this = Expression(SymEngine::mul(get_RCP(), Expression(rhs).get_RCP()));
      return *this;
    }
 
 
    template <typename NumberType>
    Expression &
    Expression::operator/=(const NumberType &rhs)
    {
      *this = Expression(SymEngine::div(get_RCP(), Expression(rhs).get_RCP()));
      return *this;
    }
 
 
    template <typename ResultType>
    Expression::operator ResultType() const
    {
      return static_cast<ResultType>(get_expression());
    }
 
  } // namespace SD
} // namespace Differentiation
 
 
// Specializations from numbers.h
// These are required in order to make the SD::Expression class a compatible
// number for the Tensor and SymmetricTensor classes
 
// Forward declarations:
template <int rank_, int dim, typename Number>
class Tensor;
 
namespace internal
{
  template <>
  struct NumberType<Differentiation::SD::Expression>
  {
    static const Differentiation::SD::Expression &
    value(const Differentiation::SD::Expression &t)
    {
      return t;
    }
 
    template <typename T, typename = std::enable_if_t<std::is_arithmetic_v<T>>>
    static Differentiation::SD::Expression
    value(const T &t)
    {
      return Differentiation::SD::Expression(t);
    }
 
    template <typename T, typename = std::enable_if_t<std::is_arithmetic_v<T>>>
    static Differentiation::SD::Expression
    value(T &&t)
    {
      return Differentiation::SD::Expression(t);
    }
 
    template <typename T, typename = std::enable_if_t<std::is_arithmetic_v<T>>>
    static Differentiation::SD::Expression
    value(const std::complex<T> &t)
    {
      return Differentiation::SD::Expression(t);
    }
 
    template <typename T, int dim>
    static Differentiation::SD::Expression
    value(const Tensor<0, dim, T> &t)
    {
      return Differentiation::SD::Expression(static_cast<T>(t));
    }
 
    template <typename T, int dim>
    static Differentiation::SD::Expression
    value(const Tensor<0, dim, std::complex<T>> &t)
    {
      return Differentiation::SD::Expression(static_cast<std::complex<T>>(t));
    }
  };
} // namespace internal
 
 
namespace numbers
{
  template <>
  inline bool
  value_is_zero(const Differentiation::SD::Expression &value)
  {
    if (SymEngine::is_a_Number(value.get_value()))
      {
        const SymEngine::RCP<const SymEngine::Number> number_rcp =
          SymEngine::rcp_static_cast<const SymEngine::Number>(value.get_RCP());
        return number_rcp->is_zero();
      }
 
    return false;
  }
 
  template <>
  inline bool
  values_are_equal(const Differentiation::SD::Expression &value_1,
                   const Differentiation::SD::Expression &value_2)
  {
    return (value_1.get_value().__cmp__(value_2.get_value()) == 0);
  }
 
  template <>
  inline bool
  value_is_less_than(const Differentiation::SD::Expression &value_1,
                     const Differentiation::SD::Expression &value_2)
  {
    return (value_1.get_value().__cmp__(value_2.get_value()) == -1);
  }
} // namespace numbers
 
 
#  endif // DOXYGEN
 
DEAL_II_NAMESPACE_CLOSE
 
#endif // DEAL_II_WITH_SYMENGINE
 
#endif // dealii_differentiation_sd_symengine_number_types_h