symengine_number_types.h

// ---------------------------------------------------------------------
//
// Copyright (C) 2019 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_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 <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 = typename std::enable_if<
                  std::is_arithmetic<NumberType>::value>::type>
      explicit Expression(const NumberType &value);

      template <typename NumberType,
                typename = typename std::enable_if<
                  std::is_arithmetic<NumberType>::value>::type>
      explicit Expression(const std::complex<NumberType> &value);

      Expression(const SymEngine::integer_class &value);

      template <typename NumberType,
                typename = typename std::enable_if<
                  std::is_integral<NumberType>::value>::type>
      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
      serialize(Archive &archive, const unsigned int version);



      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);

      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 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 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 = typename std::enable_if<
                std::is_constructible<Expression, NumberType>::value>::type>
    inline Expression
    operator+(const NumberType &lhs, const Expression &rhs)
    {
      return Expression(lhs) + rhs;
    }

    template <typename NumberType,
              typename = typename std::enable_if<
                std::is_constructible<Expression, NumberType>::value>::type>
    inline Expression
    operator+(const Expression &lhs, const NumberType &rhs)
    {
      return lhs + Expression(rhs);
    }

    template <typename NumberType,
              typename = typename std::enable_if<
                std::is_constructible<Expression, NumberType>::value>::type>
    inline Expression
    operator-(const NumberType &lhs, const Expression &rhs)
    {
      return Expression(lhs) - rhs;
    }

    template <typename NumberType,
              typename = typename std::enable_if<
                std::is_constructible<Expression, NumberType>::value>::type>
    inline Expression
    operator-(const Expression &lhs, const NumberType &rhs)
    {
      return lhs - Expression(rhs);
    }

    template <typename NumberType,
              typename = typename std::enable_if<
                std::is_constructible<Expression, NumberType>::value>::type>
    inline Expression operator*(const NumberType &lhs, const Expression &rhs)
    {
      return Expression(lhs) * rhs;
    }

    template <typename NumberType,
              typename = typename std::enable_if<
                std::is_constructible<Expression, NumberType>::value>::type>
    inline Expression operator*(const Expression &lhs, const NumberType &rhs)
    {
      return lhs * Expression(rhs);
    }

    template <typename NumberType,
              typename = typename std::enable_if<
                std::is_constructible<Expression, NumberType>::value>::type>
    inline Expression
    operator/(const NumberType &lhs, const Expression &rhs)
    {
      return Expression(lhs) / rhs;
    }

    template <typename NumberType,
              typename = typename std::enable_if<
                std::is_constructible<Expression, NumberType>::value>::type>
    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::serialize(Archive &ar, const unsigned int /*version*/)
    {
      // This is a bit tricky... The SymEngine expression class is not
      // serializable, and we're not sure if we're writing out or
      // reading in here. So what we'll do is try to synchronise the
      // current data with the stream, predicting if we're likely trying
      // to write out data. We do this by checking if we're currently working
      // with an object that has a trivial state. If not, then we're likely
      // to output this object's contents.
      const Expression default_constructed;
      const bool       likely_writing_out =
        (get_RCP()->__eq__(*default_constructed.get_RCP()) == false);

      std::stringstream sstream;
      std::string       expr = default_constructed.get_RCP()->__str__();

      // Output
      if (likely_writing_out)
        {
          sstream << *this;
          expr = sstream.str();
        }

      // Serialise
      ar &expr;

      // Input
      if (!likely_writing_out)
        {
          sstream.clear();
          sstream << expr;
          sstream >> *this;

          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 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 = typename std::enable_if<std::is_arithmetic<T>::value>::type>
    static Differentiation::SD::Expression
    value(const T &t)
    {
      return Differentiation::SD::Expression(t);
    }

    template <
      typename T,
      typename = typename std::enable_if<std::is_arithmetic<T>::value>::type>
    static Differentiation::SD::Expression
    value(T &&t)
    {
      return Differentiation::SD::Expression(t);
    }

    template <
      typename T,
      typename = typename std::enable_if<std::is_arithmetic<T>::value>::type>
    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