Reference documentation for deal.II version 9.5.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\}}\)
Loading...
Searching...
No Matches
symengine_number_visitor_internal.h
Go to the documentation of this file.
1// ---------------------------------------------------------------------
2//
3// Copyright (C) 2020 - 2023 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# 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
41namespace 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 */
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 (const auto &expression : intermediate_symbols_exprs)
668 {
669 const SymEngine::RCP<const SymEngine::Basic> &cse_symbol =
670 expression.first;
671 const SymEngine::RCP<const SymEngine::Basic> &cse_expr =
672 expression.second;
673 Assert(substitution_value_map.find(cse_symbol) ==
674 substitution_value_map.end(),
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,
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(),
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:472
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:473
#define Assert(cond, exc)
static ::ExceptionBase & ExcInternalError()
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