Reference documentation for deal.II version 9.1.1
\(\newcommand{\dealcoloneq}{\mathrel{\vcenter{:}}=}\)
symengine_number_types.cc
1 // ---------------------------------------------------------------------
2 //
3 // Copyright (C) 2019 by the deal.II authors
4 //
5 // This file is part of the deal.II library.
6 //
7 // The deal.II library is free software; you can use it, redistribute
8 // it, and/or modify it under the terms of the GNU Lesser General
9 // Public License as published by the Free Software Foundation; either
10 // version 2.1 of the License, or (at your option) any later version.
11 // The full text of the license can be found in the file LICENSE at
12 // the top level of the deal.II distribution.
13 //
14 // ---------------------------------------------------------------------
15 
16 #include <deal.II/base/config.h>
17 
18 #ifdef DEAL_II_WITH_SYMENGINE
19 
20 # include <deal.II/base/signaling_nan.h>
21 
22 # include <deal.II/differentiation/sd/symengine_number_types.h>
23 # include <deal.II/differentiation/sd/symengine_types.h>
24 # include <deal.II/differentiation/sd/symengine_utilities.h>
25 
26 # include <symengine/complex_double.h>
27 # include <symengine/logic.h>
28 # include <symengine/number.h>
29 # include <symengine/parser.h>
30 # include <symengine/real_double.h>
31 # include <symengine/symbol.h>
32 # include <symengine/symengine_exception.h>
33 
34 # include <string>
35 
36 DEAL_II_NAMESPACE_OPEN
37 
38 namespace Differentiation
39 {
40  namespace SD
41  {
42  namespace SE = ::SymEngine;
43 
44 
45  /* ---------------------------- Constructors -------------------------- */
46 
47 
49  : expression()
50  {}
51 
52 
53  Expression::Expression(const bool &value)
54  : expression(SE::boolean(value))
55  {}
56 
57 
58  Expression::Expression(const SE::integer_class &value)
59  : expression(value)
60  {}
61 
62 
63  Expression::Expression(const SE::rational_class &value)
64  : expression(value)
65  {}
66 
67 
69  const Expression &expression_if_true,
70  const Expression &expression_if_false)
71  {
72  Assert(SE::is_a_Boolean(condition.get_value()),
73  ExcMessage(
74  "The conditional expression must return a boolean type."));
75 
76  const SE::RCP<const SE::Boolean> condition_rcp =
77  SE::rcp_static_cast<const SE::Boolean>(condition.get_RCP());
78  expression =
79  SE::piecewise({{expression_if_true.get_RCP(), condition_rcp},
80  {expression_if_false.get_RCP(), SE::boolTrue}});
81  }
82 
83 
84  Expression::Expression(const std::vector<std::pair<Expression, Expression>>
85  & condition_expression,
86  const Expression &expression_otherwise)
87  {
88  SE::PiecewiseVec piecewise_function;
89  piecewise_function.reserve(condition_expression.size() + 1);
90 
91  // Add tested conditional entries
92  for (const auto &entry : condition_expression)
93  {
94  Assert(SE::is_a_Boolean(entry.first.get_value()),
95  ExcMessage(
96  "The conditional expression must return a boolean type."));
97  piecewise_function.push_back(
98  {entry.second.get_RCP(),
99  SE::rcp_static_cast<const SE::Boolean>(entry.first.get_RCP())});
100  }
101 
102  // Add default value
103  piecewise_function.push_back(
104  {expression_otherwise.get_RCP(), SE::boolTrue});
105 
106  // Initialize
107  expression = SE::piecewise(std::move(piecewise_function));
108  }
109 
110 
111  Expression::Expression(const std::vector<std::pair<Expression, Expression>>
112  &condition_expression)
113  {
114  // Use the other constructor with a fatal termination point
115  // ensuring that an error is thrown if none of the conditions
116  // are met.
117  *this = Expression(condition_expression,
118  Expression(numbers::signaling_nan<double>()));
119  }
120 
121 
122  Expression::Expression(const char *symbol)
123  : expression(SE::symbol(symbol))
124  {}
125 
126 
127  Expression::Expression(const std::string &str,
128  const bool parse_as_expression)
129  {
130  try
131  {
132  expression = (parse_as_expression ?
133  SE::parse(str) // The string is a symbolic "name"
134  :
135  SE::rcp_static_cast<const SE::Basic>(SE::symbol(
136  str))); // The string is a symbolic "expression"
137  }
138  catch (...)
139  {
141  }
142  }
143 
144 
145  Expression::Expression(const std::string & symbol_func,
146  const types::symbol_vector &arguments)
147  : expression(SE::function_symbol(
148  symbol_func,
149  Utilities::convert_expression_vector_to_basic_vector(arguments)))
150  {}
151 
152 
153  Expression::Expression(const SE::Expression &rhs)
154  : expression(rhs)
155  {}
156 
157 
158  Expression::Expression(const SE::RCP<const SE::Basic> &rhs)
159  : expression(rhs)
160  {}
161 
162 
163  Expression::Expression(SE::RCP<const SE::Basic> &&rhs)
164  : expression(rhs)
165  {}
166 
167 
168  /* ------------------------------ Utilities ---------------------------- */
169 
170 
171  Expression &
172  Expression::parse(const std::string &expression)
173  {
174  *this = Expression(expression, true /*parse_as_expression*/);
175  return *this;
176  }
177 
178 
179  std::ostream &
180  Expression::print(std::ostream &os) const
181  {
182  os << *this;
183  return os;
184  }
185 
186 
187  void
188  Expression::save(std::ostream &os) const
189  {
190  // We write each expression on a new line.
191  // Note: SymEngine outputs a non-terminating string
192  os << *this;
193  os << std::endl;
194  }
195 
196 
197  void
198  Expression::load(std::istream &is)
199  {
200  // Need to make sure that we read the entire line in,
201  // and then subsequently parse it.
202  std::string expr;
203  std::getline(is, expr);
204  Assert(!is.bad(), ExcIO());
205  parse(expr);
206  }
207 
208 
209  /* ------------------------------- Values ----------------------------- */
210 
211 
212  const SE::Expression &
214  {
215  return expression;
216  }
217 
218 
219  SE::Expression &
221  {
222  return expression;
223  }
224 
225 
226  const SE::Basic &
228  {
229  return *get_RCP();
230  }
231 
232 
233  const SE::RCP<const SE::Basic> &
235  {
236  return get_expression().get_basic();
237  }
238 
239 
240  /* --------------------------- Differentiation ------------------------- */
241 
242 
243  Expression
244  Expression::differentiate(const SE::RCP<const SE::Symbol> &symbol) const
245  {
246  return Expression(SE::diff(get_RCP(), symbol));
247  }
248 
249 
250  Expression
251  Expression::differentiate(const SE::RCP<const SE::Basic> &symbol) const
252  {
253  // Potential symbol
254  return Expression(SE::sdiff(get_RCP(), symbol));
255  }
256 
257 
258  Expression
260  {
261  return differentiate(symbol.get_RCP());
262  }
263 
264 
265  /* ------------- Conversion operators ------------------------- */
266 
267 
268  Expression::operator const SE::Expression &() const
269  {
270  return get_expression();
271  }
272 
273 
274  Expression::operator const SE::RCP<const SE::Basic> &() const
275  {
276  return get_expression().get_basic();
277  }
278 
279 
280  /* ------------- Dictionary-based substitution ------------------------- */
281 
282 
283  Expression
285  const types::substitution_map &substitution_values) const
286  {
287  return Expression(get_expression().subs(
288  Utilities::convert_expression_map_to_basic_map(substitution_values)));
289  }
290 
291 
292  Expression
294  const Expression &value) const
295  {
296  Assert(SE::is_a<SE::Symbol>(symbol.get_value()),
297  ExcMessage(
298  "Substitution with a number that does not represent a symbol."));
299 
300  types::substitution_map sub_vals;
301  sub_vals[symbol] = value;
302  return substitute(sub_vals);
303  }
304 
305 
306  /* -------------------- Math and relational operators ------------------ */
307 
308 
309  Expression &
311  {
312  if (this != &rhs)
313  this->expression = rhs.get_expression();
314 
315  return *this;
316  }
317 
318 
319  Expression &
321  {
322  if (this != &rhs)
323  this->expression = std::move(rhs.expression);
324 
325  return *this;
326  }
327 
328 
329  Expression
331  {
332  return Expression(-get_expression());
333  }
334 
335 
336  Expression &
338  {
339  this->expression += rhs.get_expression();
340  return *this;
341  }
342 
343 
344  Expression &
346  {
347  this->expression -= rhs.get_expression();
348  return *this;
349  }
350 
351 
352  Expression &
354  {
355  this->expression *= rhs.get_expression();
356  return *this;
357  }
358 
359 
360  Expression &
362  {
363  this->expression /= rhs.get_expression();
364  return *this;
365  }
366 
367 
368  std::ostream &
369  operator<<(std::ostream &stream, const Expression &expr)
370  {
371  stream << expr.get_expression();
372  return stream;
373  }
374 
375 
376  std::istream &
377  operator>>(std::istream &stream, Expression &expr)
378  {
379  std::string str;
380  stream >> str;
381  expr.parse(str);
382  return stream;
383  }
384 
385 
386  Expression
387  operator==(const Expression &lhs, const Expression &rhs)
388  {
389  return Expression(SE::Eq(lhs.get_RCP(), rhs.get_RCP()));
390  }
391 
392 
393  Expression
394  operator!=(const Expression &lhs, const Expression &rhs)
395  {
396  return Expression(SE::Ne(lhs.get_RCP(), rhs.get_RCP()));
397  }
398 
399 
400  Expression
401  operator<(const Expression &lhs, const Expression &rhs)
402  {
403  return Expression(SE::Lt(lhs.get_RCP(), rhs.get_RCP()));
404  }
405 
406 
407  Expression
408  operator>(const Expression &lhs, const Expression &rhs)
409  {
410  return Expression(SE::Gt(lhs.get_RCP(), rhs.get_RCP()));
411  }
412 
413 
414  Expression
415  operator<=(const Expression &lhs, const Expression &rhs)
416  {
417  return Expression(SE::Le(lhs.get_RCP(), rhs.get_RCP()));
418  }
419 
420 
421  Expression
422  operator>=(const Expression &lhs, const Expression &rhs)
423  {
424  return Expression(SE::Ge(lhs.get_RCP(), rhs.get_RCP()));
425  }
426 
427 
428  Expression operator!(const Expression &expression)
429  {
430  Assert(SE::is_a_Boolean(expression.get_value()),
431  ExcMessage("The expression must return a boolean type."));
432 
433  const SE::RCP<const SE::Boolean> expression_rcp =
434  SE::rcp_static_cast<const SE::Boolean>(expression.get_RCP());
435 
436  return Expression(SE::logical_not(expression_rcp));
437  }
438 
439 
440  Expression operator&(const Expression &lhs, const Expression &rhs)
441  {
442  Assert(SE::is_a_Boolean(lhs.get_value()),
443  ExcMessage("The lhs expression must return a boolean type."));
444  Assert(SE::is_a_Boolean(rhs.get_value()),
445  ExcMessage("The rhs expression must return a boolean type."));
446 
447  const SE::RCP<const SE::Boolean> lhs_rcp =
448  SE::rcp_static_cast<const SE::Boolean>(lhs.get_RCP());
449  const SE::RCP<const SE::Boolean> rhs_rcp =
450  SE::rcp_static_cast<const SE::Boolean>(rhs.get_RCP());
451 
452  return Expression(SE::logical_and({lhs_rcp, rhs_rcp}));
453  }
454 
455 
456  Expression
457  operator|(const Expression &lhs, const Expression &rhs)
458  {
459  Assert(SE::is_a_Boolean(lhs.get_value()),
460  ExcMessage("The lhs expression must return a boolean type."));
461  Assert(SE::is_a_Boolean(rhs.get_value()),
462  ExcMessage("The rhs expression must return a boolean type."));
463 
464  const SE::RCP<const SE::Boolean> lhs_rcp =
465  SE::rcp_static_cast<const SE::Boolean>(lhs.get_RCP());
466  const SE::RCP<const SE::Boolean> rhs_rcp =
467  SE::rcp_static_cast<const SE::Boolean>(rhs.get_RCP());
468 
469  return Expression(SE::logical_or({lhs_rcp, rhs_rcp}));
470  }
471 
472 
473  Expression
474  operator^(const Expression &lhs, const Expression &rhs)
475  {
476  Assert(SE::is_a_Boolean(lhs.get_value()),
477  ExcMessage("The lhs expression must return a boolean type."));
478  Assert(SE::is_a_Boolean(rhs.get_value()),
479  ExcMessage("The rhs expression must return a boolean type."));
480 
481  const SE::RCP<const SE::Boolean> lhs_rcp =
482  SE::rcp_static_cast<const SE::Boolean>(lhs.get_RCP());
483  const SE::RCP<const SE::Boolean> rhs_rcp =
484  SE::rcp_static_cast<const SE::Boolean>(rhs.get_RCP());
485 
486  return Expression(SE::logical_xor({lhs_rcp, rhs_rcp}));
487  }
488 
489 
490  Expression
491  operator&&(const Expression &lhs, const Expression &rhs)
492  {
493  return lhs & rhs;
494  }
495 
496 
497  Expression
498  operator||(const Expression &lhs, const Expression &rhs)
499  {
500  return lhs | rhs;
501  }
502 
503 
504  Expression
505  operator+(Expression lhs, const Expression &rhs)
506  {
507  lhs += rhs;
508  return lhs;
509  }
510 
511 
512  Expression
513  operator-(Expression lhs, const Expression &rhs)
514  {
515  lhs -= rhs;
516  return lhs;
517  }
518 
519 
521  {
522  lhs *= rhs;
523  return lhs;
524  }
525 
526 
527  Expression
528  operator/(Expression lhs, const Expression &rhs)
529  {
530  lhs /= rhs;
531  return lhs;
532  }
533 
534 
535  } // namespace SD
536 } // namespace Differentiation
537 
538 
539 DEAL_II_NAMESPACE_CLOSE
540 
541 #endif // DEAL_II_WITH_SYMENGINE
std::ostream & operator<<(std::ostream &stream, const Expression &expression)
Expression operator+(Expression lhs, const Expression &rhs)
Expression operator<(const Expression &lhs, const Expression &rhs)
Expression differentiate(const Expression &symbol) const
std::istream & operator>>(std::istream &stream, Expression &expression)
Expression & operator/=(const Expression &rhs)
Expression operator &&(const Expression &lhs, const Expression &rhs)
Expression operator>=(const Expression &lhs, const Expression &rhs)
static ::ExceptionBase & ExcIO()
Expression & parse(const std::string &expression)
#define AssertThrow(cond, exc)
Definition: exceptions.h:1519
std::ostream & print(std::ostream &stream) const
Expression substitute(const types::substitution_map &substitution_values) const
const SymEngine::RCP< const SymEngine::Basic > & get_RCP() const
Expression operator/(Expression lhs, const Expression &rhs)
static ::ExceptionBase & ExcMessage(std::string arg1)
Expression operator<=(const Expression &lhs, const Expression &rhs)
#define Assert(cond, exc)
Definition: exceptions.h:1407
Expression & operator*=(const Expression &rhs)
Expression operator|(const Expression &lhs, const Expression &rhs)
static ::ExceptionBase & ExcSymEngineParserError(std::string arg1)
const SymEngine::Basic & get_value() const
const SymEngine::Expression & get_expression() const
Expression & operator=(const Expression &rhs)
Expression operator>(const Expression &lhs, const Expression &rhs)
void save(std::ostream &stream) const
Expression operator*(Expression lhs, const Expression &rhs)
Expression operator!(const Expression &expression)
Definition: cuda.h:31
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 &rhs)
Expression operator-(Expression lhs, const Expression &rhs)
Expression & operator+=(const Expression &rhs)
Expression operator!=(const Expression &lhs, const Expression &rhs)
Expression operator &(const Expression &lhs, const Expression &rhs)