Reference documentation for deal.II version GIT 5983d193e2 2023-05-27 16:50:02+00:00
\(\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 - 2022 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 
23 # include <deal.II/base/exceptions.h>
24 # include <deal.II/base/numbers.h>
25 
28 
29 # include <boost/serialization/split_member.hpp>
30 
31 # include <symengine/basic.h>
32 # include <symengine/dict.h>
33 # include <symengine/symengine_exception.h>
34 # include <symengine/symengine_rcp.h>
35 # include <symengine/visitor.h>
36 
37 
39 
40 
41 namespace Differentiation
42 {
43  namespace SD
44  {
45  namespace internal
46  {
54  template <typename ReturnType, typename ExpressionType>
56  {
58  std::vector<std::pair<SD::Expression, SD::Expression>>;
59 
60  public:
61  /*
62  * Constructor.
63  */
64  CSEDictionaryVisitor() = default;
65 
66  /*
67  * Destructor.
68  */
69  virtual ~CSEDictionaryVisitor() = default;
70 
83  void
84  init(const types::symbol_vector &dependent_functions);
85 
123  void
124  call(ReturnType * output_values,
125  const types::symbol_vector &independent_symbols,
126  const ReturnType * substitution_values);
127 
133  template <class Archive>
134  void
135  save(Archive &archive, const unsigned int version) const;
136 
142  template <class Archive>
143  void
144  load(Archive &archive, const unsigned int version);
145 
146 # ifdef DOXYGEN
152  template <class Archive>
153  void
154  serialize(Archive &archive, const unsigned int version);
155 # else
156  // This macro defines the serialize() method that is compatible with
157  // the templated save() and load() method that have been implemented.
158  BOOST_SERIALIZATION_SPLIT_MEMBER()
159 # endif
160 
166  template <typename StreamType>
167  void
168  print(StreamType &stream) const;
169 
173  bool
174  executed() const;
175 
180  unsigned int
182 
186  unsigned int
188 
189  protected:
201  void
202  init(const SymEngine::vec_basic &dependent_functions);
203 
228  void
229  call(ReturnType * output_values,
230  const SymEngine::vec_basic &independent_symbols,
231  const ReturnType * substitution_values);
232 
233  private:
234  // Note: It would be more efficient to store this data in native
235  // SymEngine types, as it would prevent some copying of the data
236  // structures. However, this makes serialization more difficult,
237  // so we use our own serializable types instead, and lose a bit
238  // of efficiency.
239 
244 
249  };
250 
251 
252 
262  template <typename ReturnType, typename ExpressionType>
264  : public SymEngine::BaseVisitor<
265  DictionarySubstitutionVisitor<ReturnType, ExpressionType>>
266  {
267  public:
268  /*
269  * Constructor.
270  */
272 
273  /*
274  * Destructor.
275  */
276  virtual ~DictionarySubstitutionVisitor() override = default;
277 
301  void
303  const Expression & dependent_function,
304  const bool use_cse = false);
305 
328  // The following definition is required due to base class CRTP.
329  void
330  init(const SymEngine::vec_basic &independent_symbols,
331  const SymEngine::Basic & dependent_function,
332  const bool use_cse = false);
333 
354  void
357  const bool use_cse = false);
358 
359 
381  // The following definition is required due to base class CRTP.
382  void
383  init(const SymEngine::vec_basic &independent_symbols,
384  const SymEngine::vec_basic &dependent_functions,
385  const bool use_cse = false);
386 
419  void
420  call(ReturnType *output_values, const ReturnType *substitution_values);
421 
438  // The following definition is required due to base class CRTP.
439  ReturnType
440  call(const std::vector<ReturnType> &substitution_values);
441 
447  template <class Archive>
448  void
449  save(Archive &archive, const unsigned int version) const;
450 
456  template <class Archive>
457  void
458  load(Archive &archive, const unsigned int version);
459 
460 # ifdef DOXYGEN
466  template <class Archive>
467  void
468  serialize(Archive &archive, const unsigned int version);
469 # else
470  // This macro defines the serialize() method that is compatible with
471  // the templated save() and load() method that have been implemented.
472  BOOST_SERIALIZATION_SPLIT_MEMBER()
473 # endif
474 
488  template <typename StreamType>
489  void
490  print(StreamType &stream,
491  const bool print_independent_symbols = false,
492  const bool print_dependent_functions = false,
493  const bool print_cse_reductions = false) const;
494 
495 # ifndef DOXYGEN
496  // The following definitions are required due to base class CRTP.
497  // Since these are not used, and therefore not important to
498  // understand, we'll define them in the most concise manner possible.
499  // We also won't bother to document their existence, since they cannot
500  // be used.
501 # define IMPLEMENT_DSV_BVISIT(Argument) \
502  void bvisit(const Argument &) \
503  { \
504  AssertThrow(false, ExcNotImplemented()); \
505  }
506 
507  IMPLEMENT_DSV_BVISIT(SymEngine::Basic)
508  IMPLEMENT_DSV_BVISIT(SymEngine::Symbol)
509  IMPLEMENT_DSV_BVISIT(SymEngine::Constant)
510  IMPLEMENT_DSV_BVISIT(SymEngine::Integer)
511  IMPLEMENT_DSV_BVISIT(SymEngine::Rational)
512  IMPLEMENT_DSV_BVISIT(SymEngine::RealDouble)
513  IMPLEMENT_DSV_BVISIT(SymEngine::ComplexDouble)
514  IMPLEMENT_DSV_BVISIT(SymEngine::Add)
515  IMPLEMENT_DSV_BVISIT(SymEngine::Mul)
516  IMPLEMENT_DSV_BVISIT(SymEngine::Pow)
517  IMPLEMENT_DSV_BVISIT(SymEngine::Log)
518  IMPLEMENT_DSV_BVISIT(SymEngine::Sin)
519  IMPLEMENT_DSV_BVISIT(SymEngine::Cos)
520  IMPLEMENT_DSV_BVISIT(SymEngine::Tan)
521  IMPLEMENT_DSV_BVISIT(SymEngine::Csc)
522  IMPLEMENT_DSV_BVISIT(SymEngine::Sec)
523  IMPLEMENT_DSV_BVISIT(SymEngine::Cot)
524  IMPLEMENT_DSV_BVISIT(SymEngine::ASin)
525  IMPLEMENT_DSV_BVISIT(SymEngine::ACos)
526  IMPLEMENT_DSV_BVISIT(SymEngine::ATan)
527  IMPLEMENT_DSV_BVISIT(SymEngine::ATan2)
528  IMPLEMENT_DSV_BVISIT(SymEngine::ACsc)
529  IMPLEMENT_DSV_BVISIT(SymEngine::ASec)
530  IMPLEMENT_DSV_BVISIT(SymEngine::ACot)
531  IMPLEMENT_DSV_BVISIT(SymEngine::Sinh)
532  IMPLEMENT_DSV_BVISIT(SymEngine::Cosh)
533  IMPLEMENT_DSV_BVISIT(SymEngine::Tanh)
534  IMPLEMENT_DSV_BVISIT(SymEngine::Csch)
535  IMPLEMENT_DSV_BVISIT(SymEngine::Sech)
536  IMPLEMENT_DSV_BVISIT(SymEngine::Coth)
537  IMPLEMENT_DSV_BVISIT(SymEngine::ASinh)
538  IMPLEMENT_DSV_BVISIT(SymEngine::ACosh)
539  IMPLEMENT_DSV_BVISIT(SymEngine::ATanh)
540  IMPLEMENT_DSV_BVISIT(SymEngine::ACsch)
541  IMPLEMENT_DSV_BVISIT(SymEngine::ACoth)
542  IMPLEMENT_DSV_BVISIT(SymEngine::ASech)
543  IMPLEMENT_DSV_BVISIT(SymEngine::Abs)
544  IMPLEMENT_DSV_BVISIT(SymEngine::Gamma)
545  IMPLEMENT_DSV_BVISIT(SymEngine::LogGamma)
546  IMPLEMENT_DSV_BVISIT(SymEngine::Erf)
547  IMPLEMENT_DSV_BVISIT(SymEngine::Erfc)
548  IMPLEMENT_DSV_BVISIT(SymEngine::Max)
549  IMPLEMENT_DSV_BVISIT(SymEngine::Min)
550 
551 # undef IMPLEMENT_DSV_BVISIT
552 # endif // DOXYGEN
553 
554  private:
555  // Note: It would be more efficient to store this data in native
556  // SymEngine types, as it would prevent some copying of the data
557  // structures. However, this makes serialization more difficult,
558  // so we use our own serializable types instead, and lose a bit
559  // of efficiency.
560 
565 
570 
577  };
578 
579 
580 
581  /* ------------------ inline and template functions ------------------ */
582 
583 
584 # ifndef DOXYGEN
585 
586  /* -------------- CommonSubexpressionEliminationVisitor -------------- */
587 
588 
589  template <typename ReturnType, typename ExpressionType>
590  void
592  const SD::types::symbol_vector &dependent_functions)
593  {
595  dependent_functions));
596  }
597 
598 
599 
600  template <typename ReturnType, typename ExpressionType>
601  void
603  const SymEngine::vec_basic &dependent_functions)
604  {
605  // After the next call, the data stored in replacements is structured
606  // as follows:
607  //
608  // replacements[i] := [f, f(x)]
609  // replacements[i].first = intermediate function label "f"
610  // replacements[i].second = intermediate function definition "f(x)"
611  //
612  // It is to be evaluated top down (i.e. index 0 to
613  // replacements.size()), with the results going back into the
614  // substitution map for the next levels. So for each "i", "x" are the
615  // superset of the input values and the previously evaluated [f_0(x),
616  // f_1(x), ..., f_{i-1}(x)].
617  //
618  // The final result is a set of reduced expressions
619  // that must be computed after the replacement
620  // values have been computed.
621  SymEngine::vec_pair se_replacements;
622  SymEngine::vec_basic se_reduced_exprs;
623  SymEngine::cse(se_replacements, se_reduced_exprs, dependent_functions);
624 
625  intermediate_symbols_exprs =
627  se_replacements);
629  se_reduced_exprs);
630  }
631 
632 
633 
634  template <typename ReturnType, typename ExpressionType>
635  void
637  ReturnType * output_values,
638  const SD::types::symbol_vector &independent_symbols,
639  const ReturnType * substitution_values)
640  {
641  call(output_values,
643  independent_symbols),
644  substitution_values);
645  }
646 
647 
648 
649  template <typename ReturnType, typename ExpressionType>
650  void
652  ReturnType * output_values,
653  const SymEngine::vec_basic &independent_symbols,
654  const ReturnType * substitution_values)
655  {
656  Assert(n_reduced_expressions() > 0, ExcInternalError());
657 
658  // First we add the input values into the substitution map...
659  SymEngine::map_basic_basic substitution_value_map;
660  for (unsigned i = 0; i < independent_symbols.size(); ++i)
661  substitution_value_map[independent_symbols[i]] =
662  static_cast<const SymEngine::RCP<const SymEngine::Basic> &>(
663  ExpressionType(substitution_values[i]));
664 
665  // ... followed by any intermediate evaluations due to the application
666  // of CSE. These are fed directly back into the substitution map...
667  for (unsigned i = 0; i < intermediate_symbols_exprs.size(); ++i)
668  {
669  const SymEngine::RCP<const SymEngine::Basic> &cse_symbol =
670  intermediate_symbols_exprs[i].first;
671  const SymEngine::RCP<const SymEngine::Basic> &cse_expr =
672  intermediate_symbols_exprs[i].second;
673  Assert(substitution_value_map.find(cse_symbol) ==
674  substitution_value_map.end(),
675  ExcMessage(
676  "Reduced symbol already appears in substitution map. "
677  "Is there a clash between the reduced symbol name and "
678  "the symbol used for an independent variable?"));
679  substitution_value_map[cse_symbol] =
680  static_cast<const SymEngine::RCP<const SymEngine::Basic> &>(
681  ExpressionType(ExpressionType(cse_expr)
682  .template substitute_and_evaluate<ReturnType>(
683  substitution_value_map)));
684  }
685 
686  // ... followed by the final reduced expressions
687  for (unsigned i = 0; i < reduced_exprs.size(); ++i)
688  output_values[i] = ExpressionType(reduced_exprs[i])
689  .template substitute_and_evaluate<ReturnType>(
690  substitution_value_map);
691  }
692 
693 
694 
695  template <typename ReturnType, typename ExpressionType>
696  template <class Archive>
697  void
699  Archive &ar,
700  const unsigned int /*version*/) const
701  {
702  // The reduced expressions depend on the intermediate expressions,
703  // so we serialize the latter before the former.
704  ar &intermediate_symbols_exprs;
705  ar &reduced_exprs;
706  }
707 
708 
709 
710  template <typename ReturnType, typename ExpressionType>
711  template <class Archive>
712  void
714  Archive &ar,
715  const unsigned int /*version*/)
716  {
717  Assert(intermediate_symbols_exprs.empty(), ExcInternalError());
718  Assert(reduced_exprs.empty(), ExcInternalError());
719 
720  // The reduced expressions depend on the intermediate expressions,
721  // so we deserialize the latter before the former.
722  ar &intermediate_symbols_exprs;
723  ar &reduced_exprs;
724  }
725 
726 
727 
728  template <typename ReturnType, typename ExpressionType>
729  template <typename StreamType>
730  void
732  StreamType &stream) const
733  {
734  stream << "Common subexpression elimination: \n";
735  stream << " Intermediate reduced expressions: \n";
736  for (unsigned i = 0; i < intermediate_symbols_exprs.size(); ++i)
737  {
738  const SymEngine::RCP<const SymEngine::Basic> &cse_symbol =
739  intermediate_symbols_exprs[i].first;
740  const SymEngine::RCP<const SymEngine::Basic> &cse_expr =
741  intermediate_symbols_exprs[i].second;
742  stream << " " << i << ": " << cse_symbol << " = " << cse_expr
743  << '\n';
744  }
745 
746  stream << " Final reduced expressions for dependent variables: \n";
747  for (unsigned i = 0; i < reduced_exprs.size(); ++i)
748  stream << " " << i << ": " << reduced_exprs[i] << '\n';
749 
750  stream << std::flush;
751  }
752 
753 
754 
755  template <typename ReturnType, typename ExpressionType>
756  bool
758  {
759  // For dictionary substitution, the CSE algorithm moves
760  // ownership of the dependent function expression definition
761  // to the entries in reduced_exprs. So its size thus determines
762  // whether CSE has been executed or not.
763  return (n_reduced_expressions() > 0) ||
764  (n_intermediate_expressions() > 0);
765  }
766 
767 
768 
769  template <typename ReturnType, typename ExpressionType>
770  unsigned int
771  CSEDictionaryVisitor<ReturnType,
772  ExpressionType>::n_intermediate_expressions() const
773  {
774  return intermediate_symbols_exprs.size();
775  }
776 
777 
778 
779  template <typename ReturnType, typename ExpressionType>
780  unsigned int
782  const
783  {
784  return reduced_exprs.size();
785  }
786 
787 
788 
789  /* ------------------ DictionarySubstitutionVisitor ------------------ */
790 
791 
792  template <typename ReturnType, typename ExpressionType>
793  void
795  const types::symbol_vector &inputs,
796  const SD::Expression & output,
797  const bool use_cse)
798  {
799  init(inputs, types::symbol_vector{output}, use_cse);
800  }
801 
802 
803 
804  template <typename ReturnType, typename ExpressionType>
805  void
807  const SymEngine::vec_basic &inputs,
808  const SymEngine::Basic & output,
809  const bool use_cse)
810  {
812  SD::Expression(output.rcp_from_this()),
813  use_cse);
814  }
815 
816 
817 
818  template <typename ReturnType, typename ExpressionType>
819  void
821  const SymEngine::vec_basic &inputs,
822  const SymEngine::vec_basic &outputs,
823  const bool use_cse)
824  {
827  use_cse);
828  }
829 
830 
831 
832  template <typename ReturnType, typename ExpressionType>
833  void
835  const types::symbol_vector &inputs,
836  const types::symbol_vector &outputs,
837  const bool use_cse)
838  {
839  independent_symbols.clear();
840  dependent_functions.clear();
841 
842  independent_symbols = inputs;
843 
844  // Perform common subexpression elimination if requested
845  // Note: After this is done, the results produced by
846  // dependent_functions and cse.reduced_exprs should be
847  // the same. We could keep the former so that we can print
848  // out the original expressions if we wish to do so.
849  if (use_cse == false)
850  dependent_functions = outputs;
851  else
852  {
853  cse.init(outputs);
854  }
855  }
856 
857 
858 
859  template <typename ReturnType, typename ExpressionType>
860  ReturnType
862  const std::vector<ReturnType> &substitution_values)
863  {
864  Assert(
865  dependent_functions.size() == 1,
866  ExcMessage(
867  "Cannot use this call function when more than one symbolic expression is to be evaluated."));
868  Assert(
869  substitution_values.size() == independent_symbols.size(),
870  ExcMessage(
871  "Input substitution vector does not match size of symbol vector."));
872 
873  ReturnType out = ::internal::NumberType<ReturnType>::value(0.0);
874  call(&out, substitution_values.data());
875  return out;
876  }
877 
878 
879 
880  template <typename ReturnType, typename ExpressionType>
881  void
883  ReturnType * output_values,
884  const ReturnType *substitution_values)
885  {
886  // Check to see if CSE has been performed
887  if (cse.executed())
888  {
889  cse.call(output_values, independent_symbols, substitution_values);
890  }
891  else
892  {
893  // Build a substitution map.
894  SymEngine::map_basic_basic substitution_value_map;
895  for (unsigned i = 0; i < independent_symbols.size(); ++i)
896  substitution_value_map[independent_symbols[i]] =
897  static_cast<const SymEngine::RCP<const SymEngine::Basic> &>(
898  ExpressionType(substitution_values[i]));
899 
900  // Since we don't know how to definitively evaluate the
901  // input number type, we create a generic Expression
902  // with the given symbolic expression and ask it to perform
903  // substitution and evaluation for us.
904  Assert(dependent_functions.size() > 0, ExcInternalError());
905  for (unsigned i = 0; i < dependent_functions.size(); ++i)
906  output_values[i] =
907  ExpressionType(dependent_functions[i])
908  .template substitute_and_evaluate<ReturnType>(
909  substitution_value_map);
910  }
911  }
912 
913 
914 
915  template <typename ReturnType, typename ExpressionType>
916  template <class Archive>
917  void
919  Archive & ar,
920  const unsigned int version) const
921  {
922  // Add some dynamic information to determine if CSE has been used,
923  // without relying on the CSE class when deserializing.
924  // const bool used_cse = cse.executed();
925  // ar &used_cse;
926 
927  // CSE and dependent variables both require the independent
928  // symbols, so we serialize them first. The dependent variables
929  // might depend on the outcome of CSE, so we have to serialize
930  // them last.
931  ar &independent_symbols;
932  cse.save(ar, version);
933  ar &dependent_functions;
934  }
935 
936 
937 
938  template <typename ReturnType, typename ExpressionType>
939  template <class Archive>
940  void
942  Archive & ar,
943  const unsigned int version)
944  {
945  Assert(cse.executed() == false, ExcInternalError());
946  Assert(cse.n_intermediate_expressions() == 0, ExcInternalError());
947  Assert(cse.n_reduced_expressions() == 0, ExcInternalError());
948 
949  // CSE and dependent variables both require the independent
950  // symbols, so we deserialize them first. The dependent variables
951  // might depend on the outcome of CSE, so we have to deserialize
952  // them last.
953  ar &independent_symbols;
954  cse.load(ar, version);
955  ar &dependent_functions;
956  }
957 
958 
959 
960  template <typename ReturnType, typename ExpressionType>
961  template <typename StreamType>
962  void
964  StreamType &stream,
965  const bool print_independent_symbols,
966  const bool print_dependent_functions,
967  const bool print_cse_reductions) const
968  {
969  if (print_independent_symbols)
970  {
971  stream << "Independent variables: \n";
972  for (unsigned i = 0; i < independent_symbols.size(); ++i)
973  stream << " " << i << ": " << independent_symbols[i] << '\n';
974 
975  stream << std::flush;
976  }
977 
978  // Check to see if CSE has been performed
979  if (print_cse_reductions && cse.executed())
980  {
981  cse.print(stream);
982  }
983  else
984  {
985  Assert(dependent_functions.size() > 0, ExcInternalError());
986 
987  if (print_dependent_functions)
988  {
989  stream << "Dependent variables: \n";
990  for (unsigned i = 0; i < dependent_functions.size(); ++i)
991  stream << " " << i << dependent_functions[i] << '\n';
992 
993  stream << std::flush;
994  }
995  }
996  }
997 
998 # endif // DOXYGEN
999 
1000  } // namespace internal
1001  } // namespace SD
1002 } // namespace Differentiation
1003 
1004 
1006 
1007 #endif // DEAL_II_WITH_SYMENGINE
1008 
1009 #endif // dealii_differentiation_sd_symengine_number_visitor_internal_h
void serialize(Archive &archive, const unsigned int version)
void call(ReturnType *output_values, const SymEngine::vec_basic &independent_symbols, const ReturnType *substitution_values)
void call(ReturnType *output_values, const types::symbol_vector &independent_symbols, const ReturnType *substitution_values)
void save(Archive &archive, const unsigned int version) const
void init(const SymEngine::vec_basic &dependent_functions)
void init(const types::symbol_vector &dependent_functions)
void load(Archive &archive, const unsigned int version)
std::vector< std::pair< SD::Expression, SD::Expression > > symbol_vector_pair
ReturnType call(const std::vector< ReturnType > &substitution_values)
void init(const types::symbol_vector &independent_symbols, const Expression &dependent_function, const bool use_cse=false)
void init(const SymEngine::vec_basic &independent_symbols, const SymEngine::Basic &dependent_function, const bool use_cse=false)
void init(const types::symbol_vector &independent_symbols, const types::symbol_vector &dependent_functions, const bool use_cse=false)
void call(ReturnType *output_values, const ReturnType *substitution_values)
void print(StreamType &stream, const bool print_independent_symbols=false, const bool print_dependent_functions=false, const bool print_cse_reductions=false) const
void save(Archive &archive, const unsigned int version) const
void serialize(Archive &archive, const unsigned int version)
void init(const SymEngine::vec_basic &independent_symbols, const SymEngine::vec_basic &dependent_functions, const bool use_cse=false)
void load(Archive &archive, const unsigned int version)
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:474
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:475
static ::ExceptionBase & ExcInternalError()
#define Assert(cond, exc)
Definition: exceptions.h:1614
static ::ExceptionBase & ExcMessage(std::string arg1)
std::vector< std::pair< Expression, Expression > > convert_basic_pair_vector_to_expression_pair_vector(const SymEngine::vec_pair &symbol_value_vector)
SD::types::symbol_vector convert_basic_vector_to_expression_vector(const SymEngine::vec_basic &symbol_vector)
SymEngine::vec_basic convert_expression_vector_to_basic_vector(const SD::types::symbol_vector &symbol_vector)
std::vector< SD::Expression > symbol_vector
void call(const std::function< RT()> &function, internal::return_value< RT > &ret_val)
static constexpr DEAL_II_HOST_DEVICE_ALWAYS_INLINE const T & value(const T &t)
Definition: numbers.h:702