Reference documentation for deal.II version 9.3.0
\(\newcommand{\dealvcentcolon}{\mathrel{\mathop{:}}}\) \(\newcommand{\dealcoloneq}{\dealvcentcolon\mathrel{\mkern-1.2mu}=}\) \(\newcommand{\jump}[1]{\left[\!\left[ #1 \right]\!\right]}\) \(\newcommand{\average}[1]{\left\{\!\left\{ #1 \right\}\!\right\}}\)
symengine_number_visitor_internal.h
Go to the documentation of this file.
1 // ---------------------------------------------------------------------
2 //
3 // Copyright (C) 2020 - 2021 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 #ifndef dealii_differentiation_sd_symengine_number_visitor_internal_h
17 #define dealii_differentiation_sd_symengine_number_visitor_internal_h
18 
19 #include <deal.II/base/config.h>
20 
21 #ifdef DEAL_II_WITH_SYMENGINE
22 
24 // Low level
25 # include <symengine/basic.h>
26 # include <symengine/dict.h>
27 # include <symengine/symengine_exception.h>
28 # include <symengine/symengine_rcp.h>
29 
30 // Visitor
31 # include <symengine/visitor.h>
32 
34 
35 # include <deal.II/base/exceptions.h>
36 # include <deal.II/base/numbers.h>
37 
40 
41 # include <boost/serialization/split_member.hpp>
42 
43 
45 
46 
47 namespace Differentiation
48 {
49  namespace SD
50  {
51  namespace internal
52  {
60  template <typename ReturnType, typename ExpressionType>
62  {
63  using symbol_vector_pair =
64  std::vector<std::pair<SD::Expression, SD::Expression>>;
65 
66  public:
67  /*
68  * Constructor.
69  */
70  CSEDictionaryVisitor() = default;
71 
72  /*
73  * Destructor.
74  */
75  virtual ~CSEDictionaryVisitor() = default;
76 
89  void
90  init(const types::symbol_vector &dependent_functions);
91 
129  void
130  call(ReturnType * output_values,
131  const types::symbol_vector &independent_symbols,
132  const ReturnType * substitution_values);
133 
139  template <class Archive>
140  void
141  save(Archive &archive, const unsigned int version) const;
142 
148  template <class Archive>
149  void
150  load(Archive &archive, const unsigned int version);
151 
152 # ifdef DOXYGEN
153 
158  template <class Archive>
159  void
160  serialize(Archive &archive, const unsigned int version);
161 # else
162  // This macro defines the serialize() method that is compatible with
163  // the templated save() and load() method that have been implemented.
164  BOOST_SERIALIZATION_SPLIT_MEMBER()
165 # endif
166 
172  template <typename StreamType>
173  void
174  print(StreamType &stream) const;
175 
179  bool
180  executed() const;
181 
186  unsigned int
188 
192  unsigned int
193  n_reduced_expressions() const;
194 
195  protected:
207  void
208  init(const SymEngine::vec_basic &dependent_functions);
209 
234  void
235  call(ReturnType * output_values,
236  const SymEngine::vec_basic &independent_symbols,
237  const ReturnType * substitution_values);
238 
239  private:
240  // Note: It would be more efficient to store this data in native
241  // SymEngine types, as it would prevent some copying of the data
242  // structures. However, this makes serialization more difficult,
243  // so we use our own serializable types instead, and lose a bit
244  // of efficiency.
245 
250 
255  };
256 
257 
258 
268  template <typename ReturnType, typename ExpressionType>
270  : public SymEngine::BaseVisitor<
271  DictionarySubstitutionVisitor<ReturnType, ExpressionType>>
272  {
273  public:
274  /*
275  * Constructor.
276  */
277  DictionarySubstitutionVisitor() = default;
278 
279  /*
280  * Destructor.
281  */
282  virtual ~DictionarySubstitutionVisitor() override = default;
283 
307  void
308  init(const types::symbol_vector &independent_symbols,
309  const Expression & dependent_function,
310  const bool use_cse = false);
311 
334  // The following definition is required due to base class CRTP.
335  void
336  init(const SymEngine::vec_basic &independent_symbols,
337  const SymEngine::Basic & dependent_function,
338  const bool use_cse = false);
339 
360  void
361  init(const types::symbol_vector &independent_symbols,
362  const types::symbol_vector &dependent_functions,
363  const bool use_cse = false);
364 
365 
387  // The following definition is required due to base class CRTP.
388  void
389  init(const SymEngine::vec_basic &independent_symbols,
390  const SymEngine::vec_basic &dependent_functions,
391  const bool use_cse = false);
392 
425  void
426  call(ReturnType *output_values, const ReturnType *substitution_values);
427 
444  // The following definition is required due to base class CRTP.
445  ReturnType
446  call(const std::vector<ReturnType> &substitution_values);
447 
453  template <class Archive>
454  void
455  save(Archive &archive, const unsigned int version) const;
456 
462  template <class Archive>
463  void
464  load(Archive &archive, const unsigned int version);
465 
466 # ifdef DOXYGEN
467 
472  template <class Archive>
473  void
474  serialize(Archive &archive, const unsigned int version);
475 # else
476  // This macro defines the serialize() method that is compatible with
477  // the templated save() and load() method that have been implemented.
478  BOOST_SERIALIZATION_SPLIT_MEMBER()
479 # endif
480 
494  template <typename StreamType>
495  void
496  print(StreamType &stream,
497  const bool print_independent_symbols = false,
498  const bool print_dependent_functions = false,
499  const bool print_cse_reductions = false) const;
500 
501 # ifndef DOXYGEN
502  // The following definitions are required due to base class CRTP.
503  // Since these are not used, and therefore not important to
504  // understand, we'll define them in the most concise manner possible.
505  // We also won't bother to document their existence, since they cannot
506  // be used.
507 # define IMPLEMENT_DSV_BVISIT(Argument) \
508  void bvisit(const Argument &) \
509  { \
510  AssertThrow(false, ExcNotImplemented()); \
511  }
512 
513  IMPLEMENT_DSV_BVISIT(SymEngine::Basic)
514  IMPLEMENT_DSV_BVISIT(SymEngine::Symbol)
515  IMPLEMENT_DSV_BVISIT(SymEngine::Constant)
516  IMPLEMENT_DSV_BVISIT(SymEngine::Integer)
517  IMPLEMENT_DSV_BVISIT(SymEngine::Rational)
518  IMPLEMENT_DSV_BVISIT(SymEngine::RealDouble)
519  IMPLEMENT_DSV_BVISIT(SymEngine::ComplexDouble)
520  IMPLEMENT_DSV_BVISIT(SymEngine::Add)
521  IMPLEMENT_DSV_BVISIT(SymEngine::Mul)
522  IMPLEMENT_DSV_BVISIT(SymEngine::Pow)
523  IMPLEMENT_DSV_BVISIT(SymEngine::Log)
524  IMPLEMENT_DSV_BVISIT(SymEngine::Sin)
525  IMPLEMENT_DSV_BVISIT(SymEngine::Cos)
526  IMPLEMENT_DSV_BVISIT(SymEngine::Tan)
527  IMPLEMENT_DSV_BVISIT(SymEngine::Csc)
528  IMPLEMENT_DSV_BVISIT(SymEngine::Sec)
529  IMPLEMENT_DSV_BVISIT(SymEngine::Cot)
530  IMPLEMENT_DSV_BVISIT(SymEngine::ASin)
531  IMPLEMENT_DSV_BVISIT(SymEngine::ACos)
532  IMPLEMENT_DSV_BVISIT(SymEngine::ATan)
533  IMPLEMENT_DSV_BVISIT(SymEngine::ATan2)
534  IMPLEMENT_DSV_BVISIT(SymEngine::ACsc)
535  IMPLEMENT_DSV_BVISIT(SymEngine::ASec)
536  IMPLEMENT_DSV_BVISIT(SymEngine::ACot)
537  IMPLEMENT_DSV_BVISIT(SymEngine::Sinh)
538  IMPLEMENT_DSV_BVISIT(SymEngine::Cosh)
539  IMPLEMENT_DSV_BVISIT(SymEngine::Tanh)
540  IMPLEMENT_DSV_BVISIT(SymEngine::Csch)
541  IMPLEMENT_DSV_BVISIT(SymEngine::Sech)
542  IMPLEMENT_DSV_BVISIT(SymEngine::Coth)
543  IMPLEMENT_DSV_BVISIT(SymEngine::ASinh)
544  IMPLEMENT_DSV_BVISIT(SymEngine::ACosh)
545  IMPLEMENT_DSV_BVISIT(SymEngine::ATanh)
546  IMPLEMENT_DSV_BVISIT(SymEngine::ACsch)
547  IMPLEMENT_DSV_BVISIT(SymEngine::ACoth)
548  IMPLEMENT_DSV_BVISIT(SymEngine::ASech)
549  IMPLEMENT_DSV_BVISIT(SymEngine::Abs)
550  IMPLEMENT_DSV_BVISIT(SymEngine::Gamma)
551  IMPLEMENT_DSV_BVISIT(SymEngine::LogGamma)
552  IMPLEMENT_DSV_BVISIT(SymEngine::Erf)
553  IMPLEMENT_DSV_BVISIT(SymEngine::Erfc)
554  IMPLEMENT_DSV_BVISIT(SymEngine::Max)
555  IMPLEMENT_DSV_BVISIT(SymEngine::Min)
556 
557 # undef IMPLEMENT_DSV_BVISIT
558 # endif // DOXYGEN
559 
560  private:
561  // Note: It would be more efficient to store this data in native
562  // SymEngine types, as it would prevent some copying of the data
563  // structures. However, this makes serialization more difficult,
564  // so we use our own serializable types instead, and lose a bit
565  // of efficiency.
566 
571 
576 
583  };
584 
585 
586 
587  /* ------------------ inline and template functions ------------------ */
588 
589 
590 # ifndef DOXYGEN
591 
592  /* -------------- CommonSubexpressionEliminationVisitor -------------- */
593 
594 
595  template <typename ReturnType, typename ExpressionType>
596  void
598  const SD::types::symbol_vector &dependent_functions)
599  {
601  dependent_functions));
602  }
603 
604 
605 
606  template <typename ReturnType, typename ExpressionType>
607  void
609  const SymEngine::vec_basic &dependent_functions)
610  {
611  // After the next call, the data stored in replacements is structured
612  // as follows:
613  //
614  // replacements[i] := [f, f(x)]
615  // replacements[i].first = intermediate function label "f"
616  // replacements[i].second = intermediate function definition "f(x)"
617  //
618  // It is to be evaluated top down (i.e. index 0 to
619  // replacements.size()), with the results going back into the
620  // substitution map for the next levels. So for each "i", "x" are the
621  // superset of the input values and the previously evaluated [f_0(x),
622  // f_1(x), ..., f_{i-1}(x)].
623  //
624  // The final result is a set of reduced expressions
625  // that must be computed after the replacement
626  // values have been computed.
627  SymEngine::vec_pair se_replacements;
628  SymEngine::vec_basic se_reduced_exprs;
629  SymEngine::cse(se_replacements, se_reduced_exprs, dependent_functions);
630 
633  se_replacements);
635  se_reduced_exprs);
636  }
637 
638 
639 
640  template <typename ReturnType, typename ExpressionType>
641  void
643  ReturnType * output_values,
644  const SD::types::symbol_vector &independent_symbols,
645  const ReturnType * substitution_values)
646  {
647  call(output_values,
649  independent_symbols),
650  substitution_values);
651  }
652 
653 
654 
655  template <typename ReturnType, typename ExpressionType>
656  void
658  ReturnType * output_values,
659  const SymEngine::vec_basic &independent_symbols,
660  const ReturnType * substitution_values)
661  {
663 
664  // First we add the input values into the substitution map...
665  SymEngine::map_basic_basic substitution_value_map;
666  for (unsigned i = 0; i < independent_symbols.size(); ++i)
667  substitution_value_map[independent_symbols[i]] =
668  static_cast<const SymEngine::RCP<const SymEngine::Basic> &>(
669  ExpressionType(substitution_values[i]));
670 
671  // ... followed by any intermediate evaluations due to the application
672  // of CSE. These are fed directly back into the substitution map...
673  for (unsigned i = 0; i < intermediate_symbols_exprs.size(); ++i)
674  {
675  const SymEngine::RCP<const SymEngine::Basic> &cse_symbol =
677  const SymEngine::RCP<const SymEngine::Basic> &cse_expr =
678  intermediate_symbols_exprs[i].second;
679  Assert(substitution_value_map.find(cse_symbol) ==
680  substitution_value_map.end(),
681  ExcMessage(
682  "Reduced symbol already appears in substitution map. "
683  "Is there a clash between the reduced symbol name and "
684  "the symbol used for an independent variable?"));
685  substitution_value_map[cse_symbol] =
686  static_cast<const SymEngine::RCP<const SymEngine::Basic> &>(
687  ExpressionType(ExpressionType(cse_expr)
688  .template substitute_and_evaluate<ReturnType>(
689  substitution_value_map)));
690  }
691 
692  // ... followed by the final reduced expressions
693  for (unsigned i = 0; i < reduced_exprs.size(); ++i)
694  output_values[i] = ExpressionType(reduced_exprs[i])
695  .template substitute_and_evaluate<ReturnType>(
696  substitution_value_map);
697  }
698 
699 
700 
701  template <typename ReturnType, typename ExpressionType>
702  template <class Archive>
703  void
705  Archive &ar,
706  const unsigned int /*version*/) const
707  {
708  // The reduced expressions depend on the intermediate expressions,
709  // so we serialize the latter before the former.
711  ar &reduced_exprs;
712  }
713 
714 
715 
716  template <typename ReturnType, typename ExpressionType>
717  template <class Archive>
718  void
720  Archive &ar,
721  const unsigned int /*version*/)
722  {
725 
726  // The reduced expressions depend on the intermediate expressions,
727  // so we deserialize the latter before the former.
729  ar &reduced_exprs;
730  }
731 
732 
733 
734  template <typename ReturnType, typename ExpressionType>
735  template <typename StreamType>
736  void
738  StreamType &stream) const
739  {
740  stream << "Common subexpression elimination: \n";
741  stream << " Intermediate reduced expressions: \n";
742  for (unsigned i = 0; i < intermediate_symbols_exprs.size(); ++i)
743  {
744  const SymEngine::RCP<const SymEngine::Basic> &cse_symbol =
746  const SymEngine::RCP<const SymEngine::Basic> &cse_expr =
747  intermediate_symbols_exprs[i].second;
748  stream << " " << i << ": " << cse_symbol << " = " << cse_expr
749  << "\n";
750  }
751 
752  stream << " Final reduced expressions for dependent variables: \n";
753  for (unsigned i = 0; i < reduced_exprs.size(); ++i)
754  stream << " " << i << ": " << reduced_exprs[i] << "\n";
755 
756  stream << std::flush;
757  }
758 
759 
760 
761  template <typename ReturnType, typename ExpressionType>
762  bool
764  {
765  // For dictionary substitution, the CSE algorithm moves
766  // ownership of the dependent function expression definition
767  // to the entries in reduced_exprs. So its size thus determines
768  // whether CSE has been executed or not.
769  return (n_reduced_expressions() > 0) ||
771  }
772 
773 
774 
775  template <typename ReturnType, typename ExpressionType>
776  unsigned int
777  CSEDictionaryVisitor<ReturnType,
778  ExpressionType>::n_intermediate_expressions() const
779  {
780  return intermediate_symbols_exprs.size();
781  }
782 
783 
784 
785  template <typename ReturnType, typename ExpressionType>
786  unsigned int
788  const
789  {
790  return reduced_exprs.size();
791  }
792 
793 
794 
795  /* ------------------ DictionarySubstitutionVisitor ------------------ */
796 
797 
798  template <typename ReturnType, typename ExpressionType>
799  void
801  const types::symbol_vector &inputs,
802  const SD::Expression & output,
803  const bool use_cse)
804  {
805  init(inputs, types::symbol_vector{output}, use_cse);
806  }
807 
808 
809 
810  template <typename ReturnType, typename ExpressionType>
811  void
813  const SymEngine::vec_basic &inputs,
814  const SymEngine::Basic & output,
815  const bool use_cse)
816  {
818  SD::Expression(output.rcp_from_this()),
819  use_cse);
820  }
821 
822 
823 
824  template <typename ReturnType, typename ExpressionType>
825  void
827  const SymEngine::vec_basic &inputs,
828  const SymEngine::vec_basic &outputs,
829  const bool use_cse)
830  {
833  use_cse);
834  }
835 
836 
837 
838  template <typename ReturnType, typename ExpressionType>
839  void
841  const types::symbol_vector &inputs,
842  const types::symbol_vector &outputs,
843  const bool use_cse)
844  {
845  independent_symbols.clear();
846  dependent_functions.clear();
847 
848  independent_symbols = inputs;
849 
850  // Perform common subexpression elimination if requested
851  // Note: After this is done, the results produced by
852  // dependent_functions and cse.reduced_exprs should be
853  // the same. We could keep the former so that we can print
854  // out the original expressions if we wish to do so.
855  if (use_cse == false)
856  dependent_functions = outputs;
857  else
858  {
859  cse.init(outputs);
860  }
861  }
862 
863 
864 
865  template <typename ReturnType, typename ExpressionType>
866  ReturnType
868  const std::vector<ReturnType> &substitution_values)
869  {
870  Assert(
871  dependent_functions.size() == 1,
872  ExcMessage(
873  "Cannot use this call function when more than one symbolic expression is to be evaluated."));
874  Assert(
875  substitution_values.size() == independent_symbols.size(),
876  ExcMessage(
877  "Input substitution vector does not match size of symbol vector."));
878 
879  ReturnType out = ::internal::NumberType<ReturnType>::value(0.0);
880  call(&out, substitution_values.data());
881  return out;
882  }
883 
884 
885 
886  template <typename ReturnType, typename ExpressionType>
887  void
889  ReturnType * output_values,
890  const ReturnType *substitution_values)
891  {
892  // Check to see if CSE has been performed
893  if (cse.executed())
894  {
895  cse.call(output_values, independent_symbols, substitution_values);
896  }
897  else
898  {
899  // Build a substitution map.
900  SymEngine::map_basic_basic substitution_value_map;
901  for (unsigned i = 0; i < independent_symbols.size(); ++i)
902  substitution_value_map[independent_symbols[i]] =
903  static_cast<const SymEngine::RCP<const SymEngine::Basic> &>(
904  ExpressionType(substitution_values[i]));
905 
906  // Since we don't know how to definitively evaluate the
907  // input number type, we create a generic Expression
908  // with the given symbolic expression and ask it to perform
909  // substitution and evaluation for us.
910  Assert(dependent_functions.size() > 0, ExcInternalError());
911  for (unsigned i = 0; i < dependent_functions.size(); ++i)
912  output_values[i] =
913  ExpressionType(dependent_functions[i])
914  .template substitute_and_evaluate<ReturnType>(
915  substitution_value_map);
916  }
917  }
918 
919 
920 
921  template <typename ReturnType, typename ExpressionType>
922  template <class Archive>
923  void
925  Archive & ar,
926  const unsigned int version) const
927  {
928  // Add some dynamic information to determine if CSE has been used,
929  // without relying on the CSE class when deserializing.
930  // const bool used_cse = cse.executed();
931  // ar &used_cse;
932 
933  // CSE and dependent variables both require the independent
934  // symbols, so we serialize them first. The dependent variables
935  // might depend on the outcome of CSE, so we have to serialize
936  // them last.
937  ar &independent_symbols;
938  cse.save(ar, version);
939  ar &dependent_functions;
940  }
941 
942 
943 
944  template <typename ReturnType, typename ExpressionType>
945  template <class Archive>
946  void
948  Archive & ar,
949  const unsigned int version)
950  {
951  Assert(cse.executed() == false, ExcInternalError());
952  Assert(cse.n_intermediate_expressions() == 0, ExcInternalError());
953  Assert(cse.n_reduced_expressions() == 0, ExcInternalError());
954 
955  // CSE and dependent variables both require the independent
956  // symbols, so we deserialize them first. The dependent variables
957  // might depend on the outcome of CSE, so we have to deserialize
958  // them last.
959  ar &independent_symbols;
960  cse.load(ar, version);
961  ar &dependent_functions;
962  }
963 
964 
965 
966  template <typename ReturnType, typename ExpressionType>
967  template <typename StreamType>
968  void
970  StreamType &stream,
971  const bool print_independent_symbols,
972  const bool print_dependent_functions,
973  const bool print_cse_reductions) const
974  {
975  if (print_independent_symbols)
976  {
977  stream << "Independent variables: \n";
978  for (unsigned i = 0; i < independent_symbols.size(); ++i)
979  stream << " " << i << ": " << independent_symbols[i] << "\n";
980 
981  stream << std::flush;
982  }
983 
984  // Check to see if CSE has been performed
985  if (print_cse_reductions && cse.executed())
986  {
987  cse.print(stream);
988  }
989  else
990  {
991  Assert(dependent_functions.size() > 0, ExcInternalError());
992 
993  if (print_dependent_functions)
994  {
995  stream << "Dependent variables: \n";
996  for (unsigned i = 0; i < dependent_functions.size(); ++i)
997  stream << " " << i << dependent_functions[i] << "\n";
998 
999  stream << std::flush;
1000  }
1001  }
1002  }
1003 
1004 # endif // DOXYGEN
1005 
1006  } // namespace internal
1007  } // namespace SD
1008 } // namespace Differentiation
1009 
1010 
1012 
1013 #endif // DEAL_II_WITH_SYMENGINE
1014 
1015 #endif // dealii_differentiation_sd_symengine_number_visitor_internal_h
std::vector< std::pair< SD::Expression, SD::Expression > > symbol_vector_pair
void init(const types::symbol_vector &dependent_functions)
static constexpr const T & value(const T &t)
Definition: numbers.h:693
void init(const types::symbol_vector &independent_symbols, const Expression &dependent_function, const bool use_cse=false)
SD::types::symbol_vector convert_basic_vector_to_expression_vector(const SymEngine::vec_basic &symbol_vector)
void load(Archive &archive, const unsigned int version)
void save(Archive &archive, const unsigned int version) const
std::vector< std::pair< Expression, Expression > > convert_basic_pair_vector_to_expression_pair_vector(const SymEngine::vec_pair &symbol_value_vector)
#define DEAL_II_DISABLE_EXTRA_DIAGNOSTICS
Definition: config.h:408
static ::ExceptionBase & ExcMessage(std::string arg1)
#define Assert(cond, exc)
Definition: exceptions.h:1465
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:395
void call(ReturnType *output_values, const types::symbol_vector &independent_symbols, const ReturnType *substitution_values)
SymEngine::vec_basic convert_expression_vector_to_basic_vector(const SD::types::symbol_vector &symbol_vector)
void serialize(Archive &archive, const unsigned int version)
std::vector< SD::Expression > symbol_vector
void load(Archive &archive, const unsigned int version)
void call(ReturnType *output_values, const ReturnType *substitution_values)
#define DEAL_II_ENABLE_EXTRA_DIAGNOSTICS
Definition: config.h:445
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:394
void print(StreamType &stream, const bool print_independent_symbols=false, const bool print_dependent_functions=false, const bool print_cse_reductions=false) const
static ::ExceptionBase & ExcInternalError()
void save(Archive &archive, const unsigned int version) const