Reference documentation for deal.II version 9.4.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 - 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
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
36# include <deal.II/base/numbers.h>
37
40
41# include <boost/serialization/split_member.hpp>
42
43
45
46
47namespace Differentiation
48{
49 namespace SD
50 {
51 namespace internal
52 {
60 template <typename ReturnType, typename ExpressionType>
62 {
64 std::vector<std::pair<SD::Expression, SD::Expression>>;
65
66 public:
67 /*
68 * Constructor.
69 */
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
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
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 */
278
279 /*
280 * Destructor.
281 */
282 virtual ~DictionarySubstitutionVisitor() override = default;
283
307 void
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
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
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
631 intermediate_symbols_exprs =
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 {
662 Assert(n_reduced_expressions() > 0, ExcInternalError());
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 =
676 intermediate_symbols_exprs[i].first;
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(),
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.
710 ar &intermediate_symbols_exprs;
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 {
723 Assert(intermediate_symbols_exprs.empty(), ExcInternalError());
724 Assert(reduced_exprs.empty(), ExcInternalError());
725
726 // The reduced expressions depend on the intermediate expressions,
727 // so we deserialize the latter before the former.
728 ar &intermediate_symbols_exprs;
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 =
745 intermediate_symbols_exprs[i].first;
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) ||
770 (n_intermediate_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,
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(),
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
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:442
#define DEAL_II_DISABLE_EXTRA_DIAGNOSTICS
Definition: config.h:456
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:443
#define DEAL_II_ENABLE_EXTRA_DIAGNOSTICS
Definition: config.h:495
#define Assert(cond, exc)
Definition: exceptions.h:1473
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 const T & value(const T &t)
Definition: numbers.h:705