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