From: Jean-Paul Pelteret Date: Sat, 30 Mar 2019 10:19:55 +0000 (+0100) Subject: Add basis for SymEngine (SD) wrappers X-Git-Tag: v9.1.0-rc1~218^2~2 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=9386bc34a4c2b7f1821d070b3f8136eb8ac861ac;p=dealii.git Add basis for SymEngine (SD) wrappers - Types - Expression class - Utilities --- diff --git a/doc/news/changes/major/20190329Jean-PaulPelteret b/doc/news/changes/major/20190329Jean-PaulPelteret new file mode 100644 index 0000000000..d5f571d8fa --- /dev/null +++ b/doc/news/changes/major/20190329Jean-PaulPelteret @@ -0,0 +1,16 @@ +New: A new class Differentiation::SD::Expression has been added. which allows +the creation and manipulation of scalar symbolic functions using the SymEngine +library. It will later be used as the basis of symbolic tensor calculus, amongst +other things. It incorporates numerous features, which currently include: + +
+(Jean-Paul Pelteret, 2019/03/29) diff --git a/include/deal.II/differentiation/sd/symengine_number_types.h b/include/deal.II/differentiation/sd/symengine_number_types.h new file mode 100644 index 0000000000..b3a487dcb2 --- /dev/null +++ b/include/deal.II/differentiation/sd/symengine_number_types.h @@ -0,0 +1,1302 @@ +// --------------------------------------------------------------------- +// +// 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 + +#ifdef DEAL_II_WITH_SYMENGINE + +// Low level +# include +# include +# include +# include + +// Number types +# include +# include +# include +# include + +// Number operations +# include +# include +# include +# include + +// Evaluation +# include +# include +# include + +// Differentiation +# include + +# include +# include + +# include + +# include +# include +# include +# include +# include +# include + +DEAL_II_NAMESPACE_OPEN + +namespace Differentiation +{ + namespace SD + { + namespace SE = ::SymEngine; + + /** + * @addtogroup Exceptions + * @{ + */ + + /** + * An exception to indicate that the string sent to the SymEngine + * parser is not valid. + */ + 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?"); + + //@} + + /** + * A class to wrap SymEngine expressions. + * + * With this number type, SymEngine numbers can be used to perform scalar + * and tensor mathematics in deal.II. It (or the SymEngine::Expression + * class, of which it stores an instance) therefore forms the basis of + * symbolic computation via SymEngine in deal.II. With it one can + * perform symbolic differentiation and subsequent substitution with both + * scalars and deal.II's native Tensor and SymmetricTensor types. + * + * The symbolic features that this class supports includes: + * - expression parsing, + * - comparison operations, + * - logical operations, + * - math operations, + * - conditional expression construction, + * - differentiation, + * - substitution (partial and complete), and + * - serialization. + * + * A simple example of how this class may be used is as follows: + * @code + * // Constructing a symbolic expression: + * // This is a symbol, which we will treat as an argument to a symbolic + * // function. + * const Expression x("x"); + * const Expression y("y"); + * // This is a symbolic expression, which is an expression constructed + * // from individual symbols. + * const Expression f = (x + y)*(x + y); + * + * // Value substitution + * types::substitution_map substitution_map; + * substitution_map[x] = Expression(1); + * substitution_map[y] = Expression(2.5); + * const double evaluated_f = + * f.substitute_and_evaluate(substitution_map); + * // We could also have performed substitution of each individual + * // argument, if we wanted to. This means that one can partially + * // substitute an expression at any time. + * @endcode + * + * A more intricate example of conditional evaluation is as follows: + * @code + * // Construct symbolic expressions + * const Expression x("x"); + * const Expression y("y"); + * const Expression f_plus = (x + y)*(x + y); + * // Parsing expressions from a string is also possible. Its arguments + * // must have been previously declared through (and in scope). + * const Expression f_minus ("(x-y)*(x-y)", true); + * + * // Constructing a conditional expression + * const SD_number_t f((x > Expression(0.0)), f_plus, f_minus); + * + * // Value substitution + * types::substitution_map substitution_map; + * substitution_map[x] = Expression(1); + * substitution_map[y] = Expression(2.5); + * const double evaluated_f = + * f.substitute_and_evaluate(substitution_map); + * // Since the substituted value for x was greater than zero, we expect + * // that the returned result now in evaluated_f was evaluated from + * // the function f_plus. + * @endcode + * + * Lastly, here is an example using symbolic differentiation: + * @code + * // Construct symbolic expressions + * const Expression x("x"); + * const Expression f("x**2", true); + * + * // Now perform differentiation. Specifically, we differentiate the + * // function "f" with respect to the symbolic variable "x". + * // The result should be the expression "2*x". + * const Expression df_dx = f.differentiate(x); + * + * // Value substitution + * types::substitution_map substitution_map; + * substitution_map[x] = Expression(10.0); + * const double evaluated_df_dx = + * evaluated_df_dx.substitute_and_evaluate(substitution_map); + * // We can expect the above to evaluate to "2*10" which is, + * // of course, the numeric value 20. + * @endcode + * + * @author Jean-Paul Pelteret, 2019 + */ + class Expression + { + public: + /** + * @name Constructors + */ + //@{ + + /** + * Default constructor. + */ + Expression(); + + /** + * Constructor for boolean types. + * + * @note This constructor is marked as explicit so that there are no + * potential ambiguities related to implicit conversions in either user + * code or math functions that are loaded into the standard namespace. + */ + explicit Expression(const bool &value); + + /** + * Constructor for arithmetic number types. + * + * @note This constructor is marked as explicit so that there are no + * potential ambiguities related to implicit conversions in either user + * code or math functions that are loaded into the standard namespace. + */ + template ::value>::type> + explicit Expression(const NumberType &value); + + /** + * Constructor for complex numbers templated on arithmetic number types. + * + * @note This constructor is marked as explicit so that there are no + * potential ambiguities related to implicit conversions in either user + * code or math functions that are loaded into the standard namespace. + */ + template ::value>::type> + explicit Expression(const std::complex &value); + + /** + * Constructor for integer types. + */ + Expression(const SE::integer_class &value); + + /** + * Constructor for rational types. + * + * It is expected that both the @p numerator and @p denominator + * be integral types. + */ + template ::value>::type> + Expression(const NumberType &numerator, const NumberType &denominator); + + /** + * Constructor for rational types. + */ + Expression(const SE::rational_class &value); + + /** + * Constructor for a piecewise defined function. + * + * The generated expression may be interpreted as the result of the + * teniary operator, i.e. + * (condition ? expression_if_true : expression_if_false). + * + * The @p condition can be any expression that renders an expression that + * is convertible to a SymEngine::Boolean operator. This includes: + * - the logical operators of the deal.II Expression class + * - SymEngine::boolean() + * - SymEngine::contains() + * - SymEngine::Eq() + * - SymEngine::Ne() + * - SymEngine::Ge() + * - SymEngine::Gt() + * - SymEngine::Le() + * - SymEngine::Lt() + * - SymEngine::logical_and() + * - SymEngine::logical_nand() + * - SymEngine::logical_or() + * - SymEngine::logical_not() + * - SymEngine::logical_nor() + * - SymEngine::logical_xor() + * - SymEngine::logical_xnor() + * - ... + * + * An example of this constructor's use is as follows: + * @code + * const Expression x("x"); + * const Expression y("y"); + * + * // Construct a conditional expression using the symbolic variables. + * const Expression f ((x < Expression(0.0)), x+y, x-y); + * @endcode + */ + Expression(const Expression &condition, + const Expression &expression_if_true, + const Expression &expression_if_false); + + /** + * Constructor for a piecewise defined function. + * + * The generated expression may be interpreted as the result of the + * set of nested if-elseif-else statements, i.e. (in pseudo-code) + * @code + * if (condition_expression[0].first == true) + * return condition_expression[0].second; + * else if (condition_expression[1].first == true) + * return condition_expression[1].second; + * else if (...) + * return ...; + * else + * return expression_otherwise; + * @endcode + * if the input vector has more than 2 elements. + * + * This variant takes the piecewise evaluated conditions and its results + * as the first argument, and the default return value as the second + * argument. + */ + Expression(const std::vector> + & condition_expression, + const Expression &expression_otherwise); + + /** + * Constructor for a piecewise defined function. + * + * The generated expression may be interpreted as the result of the + * set of nested if-elseif statements, i.e. (in pseudo-code) + * @code + * if (condition_expression[0].first == true) + * return condition_expression[0].second; + * else if (condition_expression[1].first == true) + * return condition_expression[1].second; + * else if (...) + * return ...; + * @endcode + * if the input vector has more than 2 elements. + * + * This variant takes only the piecewise evaluated conditions and its + * results. If none of the conditions are met upon evaluation then the + * returned result will be NaN. + */ + Expression(const std::vector> + &condition_expression); + + + /** + * Constructor for symbolic types. + * + * This constructor initializes a symbolic type with a character array + * representing its symbolic value. + */ + Expression(const char *symbol); + + /** + * Constructor for symbolic types. + * + * This constructor initializes a symbolic type with a string + * representing its symbolic value. + * If the @p parse_as_expression flag is false, then the + * @p symb_expr (potentially composed of mulitple characters) will be + * interpreted as a single symbol. + * If the @p parse_as_expression flag is true, then the + * @p symb_expr will be parsed as a symbolic expression (potentially + * composed of multiple symbols, constants, etc.). + */ + Expression(const std::string &symb_expr, + const bool parse_as_expression = false); + + /** + * Constructor for function symbol types. + * + * This constructor initializes a function symbol with a string + * representing its symbolic name. + */ + Expression(const std::string & symbol_func, + const types::symbol_vector &arguments); + + /** + * Copy constructor. + */ + Expression(const Expression &rhs) = default; + + /** + * Copy constructor. + * + * @note This constructor is marked as explicit to prevent any ambiguities + * from implicit conversion when both the deal.II and SymEngine namespaces + * are imported. + */ + explicit Expression(const SE::Expression &rhs); + + /** + * Copy constructor. + * + * This allows us to create our class straight from the results of + * SymEngine operations. This is especially important for operations like + * "diff", because the returned result is not primitive, but rather a set + * of compound operations. + */ + Expression(const SE::RCP &rhs); + + /** + * Move constructor. + */ + Expression(Expression &&rhs) = default; + + /** + * Move constructor. + * + * This allows us to create our class straight from the results of + * SymEngine operations. This is especially important for operations like + * "diff", because the returned result is not primitive, but rather a set + * of compound operations. + */ + Expression(SE::RCP &&rhs); + + /** + * Destructor. + */ + virtual ~Expression() = default; + + //@} + + /** + * Utilities + */ + //@{ + + /** + * Parse an expression from a string representing a symbolic @p expression. + * This overwrites any existing value or expression that this object + * represents. + */ + Expression & + parse(const std::string &expression); + + /** + * Print the value stored by this object. + * + * Since the stored value could be one of a number of types, we leave + * SymEngine to cast and output the correct representation of the data. + */ + std::ostream & + print(std::ostream &stream) const; + + /** + * Save the value stored by this object to the @p stream. + * + * Each expression will be saved on a new line of the @p stream. + */ + void + save(std::ostream &stream) const; + + /** + * Load the value stored in the @p stream into this object. + * + * It is expected that each expression appears on its own on single + * line of @p stream. + * + * @note When loading a symbolic expression, it is imperative that + * you first create or load all of the symbolic variables used in + * the saved expression. + */ + void + load(std::istream &stream); + + /** + * Write and read the data of this object from a stream for the purpose + * of serialization. + * + * This effecively saves or loads the value stored into/out of the + * @p archive with the given @p version number into this object. + * If deserializing data, then the previous contents of this object + * are thrown away. + * + * @note When deserializing a symbolic expression, it is imperative that + * you first create or deserialize all of the symbolic variables used in + * the serialized expression. + */ + template + void + serialize(Archive &archive, const unsigned int version); + + //@} + + /** + * @name Values + */ + //@{ + + /** + * Return the value or expression that this class instance represents. + */ + const SE::Expression & + get_expression() const; + + /** + * Return the primitive SymEngine data type that stores the value or + * expression represented by this object. + */ + const SE::Basic & + get_value() const; + + /** + * Return the pointer to the primitive SymEngine data type that stores + * the value or expression represented by this object. + */ + const SE::RCP & + get_RCP() const; + + //@} + + /** + * @name Math and relational operators with (potentially) symbolic types + */ + //@{ + + /** + * Assignment operator. + * + * Sets the data of this object's @p expression equal + * to that of the @p rhs object. + */ + Expression & + operator=(const Expression &rhs); + + /** + * Assignment operator. + * + * Sets the data of this object's @p expression equal + * to that of the @p rhs object. + */ + Expression & + operator=(Expression &&rhs); + + /** + * Addition assignment. + * + * The @p rhs @p expression is added in-place to that of + * this object's @p expression. + */ + Expression & + operator+=(const Expression &rhs); + + /** + * Subtraction assignment. + * + * The @p rhs @p expression is subtracted in-place from that of + * this object's @p expression. + */ + Expression & + operator-=(const Expression &rhs); + + /** + * Multiplication assignment. + * + * This object's @p expression is multiplied in-place by that of + * the @p rhs @p expression. + */ + Expression & + operator*=(const Expression &rhs); + + /** + * Division assignment. + * + * This object's @p expression is divided in-place by that of + * the @p rhs @p expression. + */ + Expression & + operator/=(const Expression &rhs); + + //@} + + /** + * @name Math and relational operators with numeric types + */ + //@{ + + /** + * Assignment operator. + * + * Set the data of this object's @p expression equal + * to the numerical value of the @p rhs. + */ + template + Expression & + operator=(const NumberType &rhs); + + /** + * Negation operator. + * + * Return a the result of pre-multipying this object's @p expression + * by -1. + * + * @note This operation is not performed in-place. + */ + Expression + operator-() const; + + /** + * Addition assignment. + * + * The numerical value of the @p rhs is added in-place to that of + * this object's @p expression. + */ + template + Expression & + operator+=(const NumberType &rhs); + + /** + * Subtraction assignment. + * + * The numerical value of the @p rhs is subtracted in-place from that of + * this object's @p expression. + */ + template + Expression & + operator-=(const NumberType &rhs); + + /** + * Multiplication assignment. + * + * This object's @p expression is multiplied in-place by that of + * the numerical value of the @p rhs. + */ + template + Expression & + operator*=(const NumberType &rhs); + + /** + * Division assignment. + * + * This object's @p expression is divided in-place by that of + * the numerical value of the @p rhs. + */ + template + Expression & + operator/=(const NumberType &rhs); + + //@} + + /** + * @name Differentiation + */ + //@{ + + /** + * Return the derivative of this object's @p expression + * with respect to the given @p symbol. + */ + Expression + differentiate(const Expression &symbol) const; + + /** + * Return the derivative of this object's @p expression + * with respect to the given @p symbol. + */ + Expression + differentiate(const SE::RCP &symbol) const; + + /** + * Return the derivative of this object's @p expression + * with respect to the potential @p symbol. + */ + Expression + differentiate(const SE::RCP &symbol) const; + + //@} + + /** + * @name Dictionary-based substitution + */ + //@{ + + /** + * Perform substitution of all symbols found in this object's @p expression + * that match a key in the @p substitution_values map. + * + * @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 types::substitution_map &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 + * the given @p value. + * + * @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 Expression &symbol, const Expression &value) 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 + * the given @p value. + * + * @note With dictionary substitution, partial substitution is allowed + * (i.e. an incomplete substitution map can be used and the return type + * can be symbolic). + */ + template + Expression + substitute(const Expression &symbol, const NumberType &value) const; + + /** + * Full substitution and evaluation. This creates a Expression by + * symbol substitution and then immediately computes its numerical value. + * + * @note All symbols must be resolved by the substitution map in order + * for this function to return successfully. + */ + template + ReturnType + substitute_and_evaluate( + const types::substitution_map &substitution_values) const; + + //@} + + /** + * @name Conversion operators + */ + //@{ + + /** + * Conversion operator for real integer or floating point values, and + * complex integer or floating point values. + * + * @note This function is marked explicit so that the conversion must be + * performed using a static_cast. In normal use, one would have expected + * (Expression)*(double) --> (Expression) + * If this function were not marked as explicit, then we could potentially + * have + * (Expression)*(double) --> (double) + * So, to get out a value on needs to do the following: + * + * + * const NumberType val = static_cast(Expression); + * + * + * or, probably less desirably, + * + * + * const NumberType val = NumberType(Expression); + * + * + * @note If the underlying number is a custom type (i.e. encapsulated by a + * NumberWrapper), then it is necessary to derive a new class from + * Expression and define a specialized conversion operator + * that calls an Evaluator that is specialized for this custom number + * type. This could be achieved with an overriding conversion function + * in the base class, for example: + * + * + * class MyNumber : public Expression + * { + * ... + * template + * explicit operator ResultType() const + * { + * if (this->get_value()->get_type_code() == SE::NUMBER_WRAPPER) + * { + * // Implement custom evaluation function + * const ResultType result = ...; + * return result; + * } + * else + * // Call base class conversion operator + * return Expression::operator ResultType(); + * } + * } + * + */ + template + explicit operator ResultType() const; + + /** + * Conversion operator that returns the value or expression that this + * class instance represents. + */ + explicit operator const SE::Expression &() const; + + /** + * Conversion operator that returns a SymEngine reference counted pointer + * to the fundamental type. + */ + operator const SE::RCP &() const; + + //@} + + protected: + /** + * Return the value or expression that this class instance represents. + */ + SE::Expression & + get_expression(); + + private: + /** + * The value or expression that this instance of this class is to + * represent. + */ + SE::Expression expression; + }; + + /** + * @name Type traits + */ + //@{ + + /** + * A struct to indicate whether a given @p NumberType is a supported + * symbolically differentiable number or not. + * This is a specialization for the deal.II Expression class. + * + * @author Jean-Paul Pelteret, 2019 + */ + template <> + struct is_symengine_number : std::true_type + {}; + + + /** + * A struct to indicate whether a given @p NumberType is a supported + * SymEngine number or not. + * This is a specialization for the deal.II Expression class. + * + * @author Jean-Paul Pelteret, 2019 + */ + template <> + struct is_sd_number : std::true_type + {}; + + //@} + + /** + * @name Bitwise operators + */ + //@{ + + /** + * Bitwise left shift operator. + * + * This is used to output the @p expression to the input @p stream. + */ + std::ostream & + operator<<(std::ostream &stream, const Expression &expression); + + /** + * Bitwise right shift operator. + * + * This is used to read in an @p expression from the input @p stream. + */ + std::istream & + operator>>(std::istream &stream, Expression &expression); + + //@} + + /** + * @name Comparison operators + */ + //@{ + + /** + * Equality operator. + * + * Return whether the @p lhs is equal to the @p rhs. + */ + Expression + operator==(const Expression &lhs, const Expression &rhs); + + /** + * Non-equality operator. + * + * Return whether the @p lhs is not equal to the @p rhs. + */ + Expression + operator!=(const Expression &lhs, const Expression &rhs); + + /** + * Less than operator. + * + * Return whether the @p lhs is less than the @p rhs. + */ + Expression + operator<(const Expression &lhs, const Expression &rhs); + + /** + * Greater than operator. + * + * Return whether the @p lhs is greater than the @p rhs. + */ + Expression + operator>(const Expression &lhs, const Expression &rhs); + + /** + * Less than or equals operator. + * + * Return whether the @p lhs is less than, or equal to, the @p rhs. + */ + Expression + operator<=(const Expression &lhs, const Expression &rhs); + + /** + * Greater than or equals operator. + * + * Return whether the @p lhs is greater than, or equal to, the @p rhs. + */ + Expression + operator>=(const Expression &lhs, const Expression &rhs); + + //@} + + /** + * @name Logical operators + */ + //@{ + + /** + * Logical not operator. + * + * @note This operator can only be applied on boolean and conditional + * expressions. + */ + Expression operator!(const Expression &expression); + + /** + * Logical and operator. + * + * @note This operator can only be applied on boolean and conditional + * expressions. + */ + Expression operator&(const Expression &lhs, const Expression &rhs); + + /** + * Logical inclusive or operator. + * + * @note This operator can only be applied on boolean and conditional + * expressions. + */ + Expression + operator|(const Expression &lhs, const Expression &rhs); + + /** + * Logical exclusive or (xor) operator. + * + * @note This operator can only be applied on boolean and conditional + * expressions. + */ + Expression + operator^(const Expression &lhs, const Expression &rhs); + + /** + * And operator. + * + * @note This operator can only be applied on boolean and conditional + * expressions. + * This operator is a convenience wrapper for the logical and operator. + */ + Expression + operator&&(const Expression &lhs, const Expression &rhs); + + /** + * Inclusive or operator. + * + * @note This operator can only be applied on boolean and conditional + * expressions. + * This operator is a convenience wrapper for the logical or operator. + */ + Expression + operator||(const Expression &lhs, const Expression &rhs); + + //@} + + /** + * @name Mathematical operators + */ + //@{ + + /** + * Addition operator. + * + * Return the result of adding the @p rhs to the @p lhs. + */ + Expression + operator+(Expression lhs, const Expression &rhs); + + /** + * Subtraction operator. + * + * Return the result of subtracting the @p rhs from the @p lhs. + */ + Expression + operator-(Expression lhs, const Expression &rhs); + + /** + * Multiplication operator. + * + * Return the result of multiplying the @p lhs by the @p rhs. + */ + Expression operator*(Expression lhs, const Expression &rhs); + + /** + * Division operator. + * + * Return the result of dividing the @p lhs by the @p rhs. + */ + Expression + operator/(Expression lhs, const Expression &rhs); + + /** + * General addition operator. + * + * Return the result of adding the @p rhs to the @p lhs. + * The @p lhs type may be any supported number type, and the result + * is promoted to a Expression. The type conversion makes writing + * scalar expressions using Expression more natural. + */ + template ::value>::type> + inline Expression + operator+(const NumberType &lhs, const Expression &rhs) + { + return Expression(lhs) + rhs; + } + + /** + * General addition operator. + * + * Return the result of adding the @p rhs to the @p lhs. + * The @p rhs type may be any supported number type, and the result + * is promoted to a Expression. The type conversion makes writing + * scalar expressions using Expression more natural. + */ + template ::value>::type> + inline Expression + operator+(const Expression &lhs, const NumberType &rhs) + { + return lhs + Expression(rhs); + } + + /** + * General subtraction operator. + * + * Return the result of subtracting the @p rhs from the @p lhs. + * The @p lhs type may be any supported number type, and the result + * is promoted to a Expression. The type conversion makes writing + * scalar expressions using Expression more natural. + */ + template ::value>::type> + inline Expression + operator-(const NumberType &lhs, const Expression &rhs) + { + return Expression(lhs) - rhs; + } + + /** + * General subtraction operator. + * + * Return the result of subtracting the @p rhs from the @p lhs. + * The @p rhs type may be any supported number type, and the result + * is promoted to a Expression. The type conversion makes writing + * scalar expressions using Expression more natural. + */ + template ::value>::type> + inline Expression + operator-(const Expression &lhs, const NumberType &rhs) + { + return lhs - Expression(rhs); + } + + /** + * General multiplication operator. + * + * Return the result of multiplying the @p lhs by the @p rhs. + * The @p lhs type may be any supported number type, and the result + * is promoted to a Expression. The type conversion makes writing + * scalar expressions using Expression more natural. + */ + template ::value>::type> + inline Expression operator*(const NumberType &lhs, const Expression &rhs) + { + return Expression(lhs) * rhs; + } + + /** + * General multiplication operator. + * + * Return the result of multiplying the @p lhs by the @p rhs. + * The @p rhs type may be any supported number type, and the result + * is promoted to a Expression. The type conversion makes writing + * scalar expressions using Expression more natural. + */ + template ::value>::type> + inline Expression operator*(const Expression &lhs, const NumberType &rhs) + { + return lhs * Expression(rhs); + } + + /** + * General division operator. + * + * Return the result of dividing the @p lhs by the @p rhs. + * The @p lhs type may be any supported number type, and the result + * is promoted to a Expression. The type conversion makes writing + * scalar expressions using Expression more natural. + */ + template ::value>::type> + inline Expression + operator/(const NumberType &lhs, const Expression &rhs) + { + return Expression(lhs) / rhs; + } + + /** + * General division operator. + * + * Return the result of dividing the @p lhs by the @p rhs. + * The @p lhs type may be any supported number type, and the result + * is promoted to a Expression. The type conversion makes writing + * scalar expressions using Expression more natural. + */ + template ::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 + Expression::Expression(const NumberType &value) + : expression(value) + {} + + + template + Expression::Expression(const std::complex &value) + : expression(value) + {} + + + template + Expression::Expression(const NumberType &numerator, + const NumberType &denominator) + : expression(SE::Rational::from_two_ints(*SE::integer(numerator), + *SE::integer(denominator))) + {} + + + template + 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 + Expression + Expression::substitute(const Expression &symbol, + const NumberType &value) const + { + Assert(SE::is_a(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 + ReturnType + Expression::substitute_and_evaluate( + const types::substitution_map &substitution_values) const + { + return static_cast(substitute(substitution_values)); + } + + + template + Expression & + Expression::operator=(const NumberType &rhs) + { + *this = Expression(rhs); + return *this; + } + + + template + Expression & + Expression::operator+=(const NumberType &rhs) + { + *this = Expression(SE::add(get_RCP(), Expression(rhs).get_RCP())); + return *this; + } + + + template + Expression & + Expression::operator-=(const NumberType &rhs) + { + *this = Expression(SE::sub(get_RCP(), Expression(rhs).get_RCP())); + return *this; + } + + + template + Expression & + Expression::operator*=(const NumberType &rhs) + { + *this = Expression(SE::mul(get_RCP(), Expression(rhs).get_RCP())); + return *this; + } + + + template + Expression & + Expression::operator/=(const NumberType &rhs) + { + *this = Expression(SE::div(get_RCP(), Expression(rhs).get_RCP())); + return *this; + } + + + template + Expression::operator ResultType() const + { + return static_cast(get_expression()); + } + + } // namespace SD +} // namespace Differentiation + + +DEAL_II_NAMESPACE_CLOSE + + +# endif // DOXYGEN + +#endif // DEAL_II_WITH_SYMENGINE + +#endif // dealii_differentiation_sd_symengine_number_types_h diff --git a/include/deal.II/differentiation/sd/symengine_types.h b/include/deal.II/differentiation/sd/symengine_types.h new file mode 100644 index 0000000000..1f493c007c --- /dev/null +++ b/include/deal.II/differentiation/sd/symengine_types.h @@ -0,0 +1,81 @@ +// --------------------------------------------------------------------- +// +// 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_types_h +#define dealii_differentiation_sd_symengine_types_h + +#include + +#ifdef DEAL_II_WITH_SYMENGINE + +# include +# include + +DEAL_II_NAMESPACE_OPEN + +namespace Differentiation +{ + namespace SD + { + // Forward declarations + class Expression; + + + namespace types + { + namespace internal + { + /** + * A comparator for Expressions used as keys in maps. + */ + struct ExpressionKeyLess + { + bool + operator()(const SD::Expression &lhs, + const SD::Expression &rhs) const; + }; + } // namespace internal + + /** + * Type definition for a value substitution map. + * + * This serves the same purpose as a + * SymEngine::map_basic_basic, which is equivalent to a + * std::map, SE::RCP>. + */ + using substitution_map = + std::map; + + /** + * Type definition for a vector of symbols. + * + * This serves the same purpose as a SymEngine::vec_basic, + * which is equivalent to a std::vector> + */ + using symbol_vector = std::vector; + + } // namespace types + + } // namespace SD +} // namespace Differentiation + + +DEAL_II_NAMESPACE_CLOSE + +#endif // DEAL_II_WITH_SYMENGINE + +#endif // dealii_differentiation_sd_symengine_types_h diff --git a/include/deal.II/differentiation/sd/symengine_utilities.h b/include/deal.II/differentiation/sd/symengine_utilities.h new file mode 100644 index 0000000000..e9f64c5716 --- /dev/null +++ b/include/deal.II/differentiation/sd/symengine_utilities.h @@ -0,0 +1,62 @@ +// --------------------------------------------------------------------- +// +// 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_utilities_h +#define dealii_differentiation_sd_symengine_utilities_h + +#include + +#ifdef DEAL_II_WITH_SYMENGINE + +# include + +# include +# include + +DEAL_II_NAMESPACE_OPEN + +namespace Differentiation +{ + namespace SD + { + namespace SE = ::SymEngine; + + namespace Utilities + { + /** + * Convert a map of Expressions to its SymEngine counterpart. + */ + SE::map_basic_basic + convert_expression_map_to_basic_map( + const SD::types::substitution_map &substitution_map); + + /** + * Convert a vector of Expressions to its SymEngine counterpart. + */ + SE::vec_basic + convert_expression_vector_to_basic_vector( + const SD::types::symbol_vector &symbol_vector); + + } // namespace Utilities + + } // namespace SD +} // namespace Differentiation + + +DEAL_II_NAMESPACE_CLOSE + +#endif // DEAL_II_WITH_SYMENGINE + +#endif // dealii_differentiation_sd_symengine_utilities_h diff --git a/source/differentiation/CMakeLists.txt b/source/differentiation/CMakeLists.txt index 28015cc9d5..3d27c72a33 100644 --- a/source/differentiation/CMakeLists.txt +++ b/source/differentiation/CMakeLists.txt @@ -1,6 +1,6 @@ ## --------------------------------------------------------------------- ## -## Copyright (C) 2017 by the deal.II authors +## Copyright (C) 2017-2019 by the deal.II authors ## ## This file is part of the deal.II library. ## @@ -14,3 +14,4 @@ ## --------------------------------------------------------------------- ADD_SUBDIRECTORY(ad) +ADD_SUBDIRECTORY(sd) diff --git a/source/differentiation/sd/CMakeLists.txt b/source/differentiation/sd/CMakeLists.txt new file mode 100644 index 0000000000..f6730585d8 --- /dev/null +++ b/source/differentiation/sd/CMakeLists.txt @@ -0,0 +1,32 @@ +## --------------------------------------------------------------------- +## +## Copyright (C) 2017 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_DIRECTORIES(BEFORE ${CMAKE_CURRENT_BINARY_DIR}) + +SET(_src + symengine_number_types.cc + symengine_types.cc + symengine_utilities.cc + ) + +SET(_inst + ) + +FILE(GLOB _header + ${CMAKE_SOURCE_DIR}/include/deal.II/differentiation/sd/*.h + ) + +DEAL_II_ADD_LIBRARY(obj_differentiation_sd OBJECT ${_src} ${_header} ${_inst}) +EXPAND_INSTANTIATIONS(obj_differentiation_sd "${_inst}") diff --git a/source/differentiation/sd/symengine_number_types.cc b/source/differentiation/sd/symengine_number_types.cc new file mode 100644 index 0000000000..fce494713b --- /dev/null +++ b/source/differentiation/sd/symengine_number_types.cc @@ -0,0 +1,541 @@ +// --------------------------------------------------------------------- +// +// 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. +// +// --------------------------------------------------------------------- + +#include + +#ifdef DEAL_II_WITH_SYMENGINE + +# include + +# include +# include +# include + +# include +# include +# include +# include +# include +# include +# include + +# include + +DEAL_II_NAMESPACE_OPEN + +namespace Differentiation +{ + namespace SD + { + namespace SE = ::SymEngine; + + + /* ---------------------------- Constructors -------------------------- */ + + + Expression::Expression() + : expression() + {} + + + Expression::Expression(const bool &value) + : expression(SE::boolean(value)) + {} + + + Expression::Expression(const SE::integer_class &value) + : expression(value) + {} + + + Expression::Expression(const SE::rational_class &value) + : expression(value) + {} + + + Expression::Expression(const Expression &condition, + const Expression &expression_if_true, + const Expression &expression_if_false) + { + Assert(SE::is_a_Boolean(condition.get_value()), + ExcMessage( + "The conditional expression must return a boolean type.")); + + const SE::RCP condition_rcp = + SE::rcp_static_cast(condition.get_RCP()); + expression = + SE::piecewise({{expression_if_true.get_RCP(), condition_rcp}, + {expression_if_false.get_RCP(), SE::boolTrue}}); + } + + + Expression::Expression(const std::vector> + & condition_expression, + const Expression &expression_otherwise) + { + SE::PiecewiseVec piecewise_function; + piecewise_function.reserve(condition_expression.size() + 1); + + // Add tested conditional entries + for (const auto &entry : condition_expression) + { + Assert(SE::is_a_Boolean(entry.first.get_value()), + ExcMessage( + "The conditional expression must return a boolean type.")); + piecewise_function.push_back( + {entry.second.get_RCP(), + SE::rcp_static_cast(entry.first.get_RCP())}); + } + + // Add default value + piecewise_function.push_back( + {expression_otherwise.get_RCP(), SE::boolTrue}); + + // Initialize + expression = SE::piecewise(std::move(piecewise_function)); + } + + + Expression::Expression(const std::vector> + &condition_expression) + { + // Use the other constructor with a fatal termination point + // ensuring that an error is thrown if none of the conditions + // are met. + *this = Expression(condition_expression, + Expression(numbers::signaling_nan())); + } + + + Expression::Expression(const char *symbol) + : expression(SE::symbol(symbol)) + {} + + + Expression::Expression(const std::string &str, + const bool parse_as_expression) + { + try + { + expression = (parse_as_expression ? + SE::parse(str) // The string is a symbolic "name" + : + SE::rcp_static_cast(SE::symbol( + str))); // The string is a symbolic "expression" + } + catch (...) + { + AssertThrow(false, ExcSymEngineParserError(str)); + } + } + + + Expression::Expression(const std::string & symbol_func, + const types::symbol_vector &arguments) + : expression(SE::function_symbol( + symbol_func, + Utilities::convert_expression_vector_to_basic_vector(arguments))) + {} + + + Expression::Expression(const SE::Expression &rhs) + : expression(rhs) + {} + + + Expression::Expression(const SE::RCP &rhs) + : expression(rhs) + {} + + + Expression::Expression(SE::RCP &&rhs) + : expression(rhs) + {} + + + /* ------------------------------ Utilities ---------------------------- */ + + + Expression & + Expression::parse(const std::string &expression) + { + *this = Expression(expression, true /*parse_as_expression*/); + return *this; + } + + + std::ostream & + Expression::print(std::ostream &os) const + { + os << *this; + return os; + } + + + void + Expression::save(std::ostream &os) const + { + // We write each expression on a new line. + // Note: SymEngine outputs a non-terminating string + os << *this; + os << std::endl; + } + + + void + Expression::load(std::istream &is) + { + // Need to make sure that we read the entire line in, + // and then subsequently parse it. + std::string expr; + std::getline(is, expr); + Assert(!is.bad(), ExcIO()); + parse(expr); + } + + + /* ------------------------------- Values ----------------------------- */ + + + const SE::Expression & + Expression::get_expression() const + { + return expression; + } + + + SE::Expression & + Expression::get_expression() + { + return expression; + } + + + const SE::Basic & + Expression::get_value() const + { + return *get_RCP(); + } + + + const SE::RCP & + Expression::get_RCP() const + { + return get_expression().get_basic(); + } + + + /* --------------------------- Differentiation ------------------------- */ + + + Expression + Expression::differentiate(const SE::RCP &symbol) const + { + return Expression(SE::diff(get_RCP(), symbol)); + } + + + Expression + Expression::differentiate(const SE::RCP &symbol) const + { + // Potential symbol + return Expression(SE::sdiff(get_RCP(), symbol)); + } + + + Expression + Expression::differentiate(const Expression &symbol) const + { + return differentiate(symbol.get_RCP()); + } + + + /* ------------- Conversion operators ------------------------- */ + + + Expression::operator const SE::Expression &() const + { + return get_expression(); + } + + + Expression::operator const SE::RCP &() const + { + return get_expression().get_basic(); + } + + + /* ------------- Dictionary-based substitution ------------------------- */ + + + Expression + Expression::substitute( + const types::substitution_map &substitution_values) const + { + return Expression(get_expression().subs( + Utilities::convert_expression_map_to_basic_map(substitution_values))); + } + + + Expression + Expression::substitute(const Expression &symbol, + const Expression &value) const + { + Assert(SE::is_a(symbol.get_value()), + ExcMessage( + "Substitution with a number that does not represent a symbol.")); + + types::substitution_map sub_vals; + sub_vals[symbol] = value; + return substitute(sub_vals); + } + + + /* -------------------- Math and relational operators ------------------ */ + + + Expression & + Expression::operator=(const Expression &rhs) + { + if (this != &rhs) + this->expression = rhs.get_expression(); + + return *this; + } + + + Expression & + Expression::operator=(Expression &&rhs) + { + if (this != &rhs) + this->expression = std::move(rhs.expression); + + return *this; + } + + + Expression + Expression::operator-() const + { + return Expression(-get_expression()); + } + + + Expression & + Expression::operator+=(const Expression &rhs) + { + this->expression += rhs.get_expression(); + return *this; + } + + + Expression & + Expression::operator-=(const Expression &rhs) + { + this->expression -= rhs.get_expression(); + return *this; + } + + + Expression & + Expression::operator*=(const Expression &rhs) + { + this->expression *= rhs.get_expression(); + return *this; + } + + + Expression & + Expression::operator/=(const Expression &rhs) + { + this->expression /= rhs.get_expression(); + return *this; + } + + + std::ostream & + operator<<(std::ostream &stream, const Expression &expr) + { + stream << expr.get_expression(); + return stream; + } + + + std::istream & + operator>>(std::istream &stream, Expression &expr) + { + std::string str; + stream >> str; + expr.parse(str); + return stream; + } + + + Expression + operator==(const Expression &lhs, const Expression &rhs) + { + return Expression(SE::Eq(lhs.get_RCP(), rhs.get_RCP())); + } + + + Expression + operator!=(const Expression &lhs, const Expression &rhs) + { + return Expression(SE::Ne(lhs.get_RCP(), rhs.get_RCP())); + } + + + Expression + operator<(const Expression &lhs, const Expression &rhs) + { + return Expression(SE::Lt(lhs.get_RCP(), rhs.get_RCP())); + } + + + Expression + operator>(const Expression &lhs, const Expression &rhs) + { + return Expression(SE::Gt(lhs.get_RCP(), rhs.get_RCP())); + } + + + Expression + operator<=(const Expression &lhs, const Expression &rhs) + { + return Expression(SE::Le(lhs.get_RCP(), rhs.get_RCP())); + } + + + Expression + operator>=(const Expression &lhs, const Expression &rhs) + { + return Expression(SE::Ge(lhs.get_RCP(), rhs.get_RCP())); + } + + + Expression operator!(const Expression &expression) + { + Assert(SE::is_a_Boolean(expression.get_value()), + ExcMessage("The expression must return a boolean type.")); + + const SE::RCP expression_rcp = + SE::rcp_static_cast(expression.get_RCP()); + + return Expression(SE::logical_not(expression_rcp)); + } + + + Expression operator&(const Expression &lhs, const Expression &rhs) + { + Assert(SE::is_a_Boolean(lhs.get_value()), + ExcMessage("The lhs expression must return a boolean type.")); + Assert(SE::is_a_Boolean(rhs.get_value()), + ExcMessage("The rhs expression must return a boolean type.")); + + const SE::RCP lhs_rcp = + SE::rcp_static_cast(lhs.get_RCP()); + const SE::RCP rhs_rcp = + SE::rcp_static_cast(rhs.get_RCP()); + + return Expression(SE::logical_and({lhs_rcp, rhs_rcp})); + } + + + Expression + operator|(const Expression &lhs, const Expression &rhs) + { + Assert(SE::is_a_Boolean(lhs.get_value()), + ExcMessage("The lhs expression must return a boolean type.")); + Assert(SE::is_a_Boolean(rhs.get_value()), + ExcMessage("The rhs expression must return a boolean type.")); + + const SE::RCP lhs_rcp = + SE::rcp_static_cast(lhs.get_RCP()); + const SE::RCP rhs_rcp = + SE::rcp_static_cast(rhs.get_RCP()); + + return Expression(SE::logical_or({lhs_rcp, rhs_rcp})); + } + + + Expression + operator^(const Expression &lhs, const Expression &rhs) + { + Assert(SE::is_a_Boolean(lhs.get_value()), + ExcMessage("The lhs expression must return a boolean type.")); + Assert(SE::is_a_Boolean(rhs.get_value()), + ExcMessage("The rhs expression must return a boolean type.")); + + const SE::RCP lhs_rcp = + SE::rcp_static_cast(lhs.get_RCP()); + const SE::RCP rhs_rcp = + SE::rcp_static_cast(rhs.get_RCP()); + + return Expression(SE::logical_xor({lhs_rcp, rhs_rcp})); + } + + + Expression + operator&&(const Expression &lhs, const Expression &rhs) + { + return lhs & rhs; + } + + + Expression + operator||(const Expression &lhs, const Expression &rhs) + { + return lhs | rhs; + } + + + Expression + operator+(Expression lhs, const Expression &rhs) + { + lhs += rhs; + return lhs; + } + + + Expression + operator-(Expression lhs, const Expression &rhs) + { + lhs -= rhs; + return lhs; + } + + + Expression operator*(Expression lhs, const Expression &rhs) + { + lhs *= rhs; + return lhs; + } + + + Expression + operator/(Expression lhs, const Expression &rhs) + { + lhs /= rhs; + return lhs; + } + + + } // namespace SD +} // namespace Differentiation + + +DEAL_II_NAMESPACE_CLOSE + +#endif // DEAL_II_WITH_SYMENGINE diff --git a/source/differentiation/sd/symengine_types.cc b/source/differentiation/sd/symengine_types.cc new file mode 100644 index 0000000000..4c628ae238 --- /dev/null +++ b/source/differentiation/sd/symengine_types.cc @@ -0,0 +1,55 @@ +// --------------------------------------------------------------------- +// +// 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. +// +// --------------------------------------------------------------------- + +#include + +#ifdef DEAL_II_WITH_SYMENGINE + +# include +# include + +# include +# include + + +DEAL_II_NAMESPACE_OPEN + +namespace Differentiation +{ + namespace SD + { + namespace SE = ::SymEngine; + + namespace types + { + namespace internal + { + bool + ExpressionKeyLess::operator()(const SD::Expression &lhs, + const SD::Expression &rhs) const + { + return SE::RCPBasicKeyLess()(lhs.get_RCP(), rhs.get_RCP()); + } + } // namespace internal + + } // namespace types + + } // namespace SD +} // namespace Differentiation + + +DEAL_II_NAMESPACE_CLOSE + +#endif // DEAL_II_WITH_SYMENGINE diff --git a/source/differentiation/sd/symengine_utilities.cc b/source/differentiation/sd/symengine_utilities.cc new file mode 100644 index 0000000000..9697316831 --- /dev/null +++ b/source/differentiation/sd/symengine_utilities.cc @@ -0,0 +1,63 @@ +// --------------------------------------------------------------------- +// +// 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. +// +// --------------------------------------------------------------------- + +#include + +#ifdef DEAL_II_WITH_SYMENGINE + +# include +# include + +DEAL_II_NAMESPACE_OPEN + +namespace Differentiation +{ + namespace SD + { + namespace SE = ::SymEngine; + + namespace Utilities + { + SE::map_basic_basic + convert_expression_map_to_basic_map( + const SD::types::substitution_map &substitution_map) + { + SE::map_basic_basic sub_map; + for (const auto &entry : substitution_map) + sub_map[entry.first.get_RCP()] = entry.second.get_RCP(); + return sub_map; + } + + + SE::vec_basic + convert_expression_vector_to_basic_vector( + const SD::types::symbol_vector &symbol_vector) + { + SE::vec_basic symb_vec; + symb_vec.reserve(symbol_vector.size()); + for (const auto &entry : symbol_vector) + symb_vec.push_back(entry.get_RCP()); + return symb_vec; + } + + } // namespace Utilities + + } // namespace SD +} // namespace Differentiation + + +DEAL_II_NAMESPACE_CLOSE + +#endif // DEAL_II_WITH_SYMENGINE