deal.II version GIT relicensing-6839-g338455934c 2026-10-02 12:10:01+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_optimizer.cc
Go to the documentation of this file.
1// -----------------------------------------------------------------------------
2//
3// SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception OR LGPL-2.1-or-later
4// Copyright (C) 2020 - 2026 by the deal.II authors
5//
6// This file is part of the deal.II library.
7//
8// Detailed license information governing the source code and contributions
9// can be found in LICENSE.md and CONTRIBUTING.md at the top level directory.
10//
11// -----------------------------------------------------------------------------
12
13#include <deal.II/base/config.h>
14
15#ifdef DEAL_II_WITH_SYMENGINE
16
19
20# include <boost/archive/text_iarchive.hpp>
21# include <boost/archive/text_oarchive.hpp>
22
23# include <utility>
24
25
26#endif // DEAL_II_WITH_SYMENGINE
27
29
30#ifdef DEAL_II_WITH_SYMENGINE
31
32
33namespace Differentiation
34{
35 namespace SD
36 {
37 template <typename ReturnType>
39 : method(OptimizerType::dictionary)
41 , ready_for_value_extraction(false)
42 , has_been_serialized(false)
43 {}
44
45
46
47 template <typename ReturnType>
49 const enum OptimizerType &optimization_method,
50 const enum OptimizationFlags &optimization_flags)
52 {
54 }
55
56
57
58 template <typename ReturnType>
60 const BatchOptimizer<ReturnType> &other)
61 : method(other.method)
62 , flags(other.flags)
63 , independent_variables_symbols(other.independent_variables_symbols)
64 , dependent_variables_functions(other.dependent_variables_functions)
65 , dependent_variables_output(0)
66 , map_dep_expr_vec_entry(other.map_dep_expr_vec_entry)
67 , ready_for_value_extraction(false)
68 , has_been_serialized(false)
69 {}
70
71
72
73 template <typename ReturnType>
74 void
76 const BatchOptimizer<ReturnType> &other)
77 {
78 method = other.method;
79 flags = other.flags;
80 independent_variables_symbols = other.independent_variables_symbols;
81 dependent_variables_functions = other.dependent_variables_functions;
82 dependent_variables_output.clear();
83 map_dep_expr_vec_entry = other.map_dep_expr_vec_entry;
84 ready_for_value_extraction = false;
85 has_been_serialized = false;
86 }
87
88
89
90 template <typename ReturnType>
91 void
93 const enum OptimizerType &optimization_method,
94 const enum OptimizationFlags &optimization_flags)
95 {
96 Assert(
97 optimized() == false,
99 "Cannot call set_optimization_method() once the optimizer is finalized."));
100
101# ifndef HAVE_SYMENGINE_LLVM
102 if (optimization_method == OptimizerType::llvm)
103 {
105 }
106# endif
107 method = optimization_method;
108 flags = optimization_flags;
109 }
110
111
112
113 template <typename ReturnType>
114 enum OptimizerType
116 {
117 return method;
118 }
119
120
121
122 template <typename ReturnType>
125 {
126 return flags;
127 }
128
129
130
131 template <typename ReturnType>
132 bool
137
138
139
140 template <typename ReturnType>
141 bool
143 {
144 if (dependent_variables_output.size() > 0)
145 {
146 Assert(dependent_variables_output.size() ==
147 dependent_variables_functions.size(),
149 return true;
150 }
151
152 return false;
153 }
154
155
156
157 template <typename ReturnType>
158 bool
160 {
161 return ready_for_value_extraction;
162 }
163
164
165
166 template <typename ReturnType>
167 void
169 const SD::types::substitution_map &substitution_map)
170 {
171 Assert(optimized() == false,
173 "Cannot register symbols once the optimizer is finalized."));
174
175 if constexpr (running_in_debug_mode())
176 {
177 // Ensure that all of the keys in the map are actually symbolic
178 // in nature
179 for (const auto &entry : substitution_map)
180 {
181 const SD::Expression &symbol = entry.first;
182 Assert(SymEngine::is_a<SymEngine::Symbol>(*(symbol.get_RCP())),
183 ExcMessage("Key entry in map is not a symbol."));
184 }
185 }
186 // Merge the two maps, in the process ensuring that there is no
187 // duplication of symbols
188 independent_variables_symbols.insert(substitution_map.begin(),
189 substitution_map.end());
190 }
191
192
193
194 template <typename ReturnType>
195 void
197 const SymEngine::map_basic_basic &substitution_map)
198 {
199 register_symbols(
201 }
202
203
204
205 template <typename ReturnType>
206 void
208 const SD::types::symbol_vector &symbols)
209 {
210 Assert(optimized() == false,
212 "Cannot register symbols once the optimizer is finalized."));
213
214 for (const auto &symbol : symbols)
215 {
216 Assert(independent_variables_symbols.find(symbol) ==
217 independent_variables_symbols.end(),
218 ExcMessage("Symbol is already in the map."));
219 independent_variables_symbols.insert(
220 std::make_pair(symbol, SD::Expression(0.0)));
221 }
222 }
223
224
225
226 template <typename ReturnType>
227 void
229 const SymEngine::vec_basic &symbols)
230 {
231 register_symbols(
233 }
234
235
236
237 template <typename ReturnType>
240 {
241 return Utilities::extract_symbols(independent_variables_symbols);
242 }
243
244
245
246 template <typename ReturnType>
247 std::size_t
249 {
250 return independent_variables_symbols.size();
251 }
252
253
254
255 template <typename ReturnType>
256 void
258 {
259 Assert(optimized() == false,
261 "Cannot register functions once the optimizer is finalized."));
262
263 register_scalar_function(function);
264 }
265
266
267
268 template <typename ReturnType>
269 void
271 const SD::types::symbol_vector &functions)
272 {
273 Assert(optimized() == false,
275 "Cannot register functions once the optimizer is finalized."));
276
277 register_vector_functions(functions);
278 }
279
280
281
282 template <typename ReturnType>
283 void
285 const SymEngine::vec_basic &functions)
286 {
287 register_functions(
289 }
290
291
292
293 template <typename ReturnType>
296 {
297 return dependent_variables_functions;
298 }
299
300
301
302 template <typename ReturnType>
303 std::size_t
305 {
306 if (has_been_serialized == false)
307 {
308 // If we've had to augment our map after serialization, then
309 // this check, unfortunately, cannot be performed.
310 Assert(map_dep_expr_vec_entry.size() ==
311 dependent_variables_functions.size(),
313 }
314 return dependent_variables_functions.size();
315 }
316
317
318
319 template <typename ReturnType>
320 void
322 {
323 Assert(optimized() == false,
324 ExcMessage("Cannot call optimize() more than once."));
325
326 // Create and configure the optimizer
327 create_optimizer(optimizer);
328 Assert(optimizer, ExcNotInitialized());
329
330 const SD::types::symbol_vector symbol_vec =
331 Utilities::extract_symbols(independent_variables_symbols);
333 *opt = dynamic_cast<typename internal::DictionaryOptimizer<
334 ReturnType>::OptimizerType *>(optimizer.get()))
335 {
336 Assert(optimization_method() == OptimizerType::dictionary,
338 internal::OptimizerHelper<ReturnType,
340 initialize(opt,
342 symbol_vec),
344 dependent_variables_functions),
345 optimization_flags());
346 }
348 *opt = dynamic_cast<typename internal::LambdaOptimizer<
349 ReturnType>::OptimizerType *>(optimizer.get()))
350 {
351 Assert(optimization_method() == OptimizerType::lambda,
353 internal::OptimizerHelper<ReturnType,
355 initialize(opt,
357 symbol_vec),
359 dependent_variables_functions),
360 optimization_flags());
361 }
362# ifdef HAVE_SYMENGINE_LLVM
363 else if (typename internal::LLVMOptimizer<ReturnType>::OptimizerType
364 *opt = dynamic_cast<typename internal::LLVMOptimizer<
365 ReturnType>::OptimizerType *>(optimizer.get()))
366 {
367 Assert(optimization_method() == OptimizerType::llvm,
369 internal::OptimizerHelper<ReturnType,
370 internal::LLVMOptimizer<ReturnType>>::
371 initialize(opt,
373 symbol_vec),
375 dependent_variables_functions),
376 optimization_flags());
377 }
378# endif
379 else
380 {
381 AssertThrow(false, ExcMessage("Unknown optimizer type."));
382 }
383
384 // The size of the outputs is now fixed, as is the number and
385 // order of the symbols to be substituted.
386 // Note: When no optimisation is actually used (i.e. optimization_method()
387 // == off and use_symbolic_CSE() == false), we could conceptually go
388 // without this data structure. However, since the user expects to perform
389 // substitution of all dependent variables in one go, we still require it
390 // for intermediate storage of results.
391 dependent_variables_output.resize(n_dependent_variables());
392 }
393
394
395
396 template <typename ReturnType>
397 void
399 const SD::types::substitution_map &substitution_map) const
400 {
401 Assert(
402 optimized() == true,
404 "The optimizer is not configured to perform substitution. "
405 "This action can only performed after optimize() has been called."));
406 Assert(optimizer, ExcNotInitialized());
407
408 // Check that the registered symbol map and the input map are compatible
409 // with one another
410 if constexpr (running_in_debug_mode())
411 {
412 const SD::types::symbol_vector symbol_sub_vec =
413 Utilities::extract_symbols(substitution_map);
414 const SD::types::symbol_vector symbol_vec =
415 Utilities::extract_symbols(independent_variables_symbols);
416 Assert(symbol_sub_vec.size() == symbol_vec.size(),
417 ExcDimensionMismatch(symbol_sub_vec.size(),
418 symbol_vec.size()));
419 for (unsigned int i = 0; i < symbol_sub_vec.size(); ++i)
420 {
421 Assert(
422 numbers::values_are_equal(symbol_sub_vec[i], symbol_vec[i]),
424 "The input substitution map is either incomplete, or does "
425 "not match that used in the register_symbols() call."));
426 }
427 }
428
429 // Extract the values from the substitution map, and use the other
430 // function
431 const std::vector<ReturnType> values =
432 Utilities::extract_values<ReturnType>(substitution_map);
433 substitute(values);
434 }
435
436
437
438 template <typename ReturnType>
439 void
441 const SymEngine::map_basic_basic &substitution_map) const
442 {
445 }
446
447
448
449 template <typename ReturnType>
450 void
452 const SD::types::symbol_vector &symbols,
453 const std::vector<ReturnType> &values) const
454 {
455 // Zip the two vectors and use the other function call
456 // This ensures the ordering of the input vectors matches that of the
457 // stored map.
458 substitute(make_substitution_map(symbols, values));
459 }
460
461
462
463 template <typename ReturnType>
464 void
466 const SymEngine::vec_basic &symbols,
467 const std::vector<ReturnType> &values) const
468 {
470 symbols),
471 values);
472 }
473
474
475
476 template <typename ReturnType>
477 void
479 const std::vector<ReturnType> &substitution_values) const
480 {
481 Assert(
482 optimized() == true,
484 "The optimizer is not configured to perform substitution. "
485 "This action can only performed after optimize() has been called."));
486 Assert(optimizer, ExcNotInitialized());
487 Assert(substitution_values.size() == independent_variables_symbols.size(),
488 ExcDimensionMismatch(substitution_values.size(),
489 independent_variables_symbols.size()));
490
492 *opt = dynamic_cast<typename internal::DictionaryOptimizer<
493 ReturnType>::OptimizerType *>(optimizer.get()))
494 {
495 Assert(optimization_method() == OptimizerType::dictionary,
497 internal::OptimizerHelper<ReturnType,
499 substitute(opt, dependent_variables_output, substitution_values);
500 }
502 *opt = dynamic_cast<typename internal::LambdaOptimizer<
503 ReturnType>::OptimizerType *>(optimizer.get()))
504 {
505 Assert(optimization_method() == OptimizerType::lambda,
507 internal::OptimizerHelper<ReturnType,
509 substitute(opt, dependent_variables_output, substitution_values);
510 }
511# ifdef HAVE_SYMENGINE_LLVM
512 else if (typename internal::LLVMOptimizer<ReturnType>::OptimizerType
513 *opt = dynamic_cast<typename internal::LLVMOptimizer<
514 ReturnType>::OptimizerType *>(optimizer.get()))
515 {
516 Assert(optimization_method() == OptimizerType::llvm,
518 internal::OptimizerHelper<ReturnType,
519 internal::LLVMOptimizer<ReturnType>>::
520 substitute(opt, dependent_variables_output, substitution_values);
521 }
522# endif
523 else
524 {
526 }
527
528 ready_for_value_extraction = true;
529 }
530
531
532
533 template <typename ReturnType>
534 const std::vector<ReturnType> &
536 {
537 Assert(
538 values_substituted() == true,
540 "The optimizer is not configured to perform evaluation. "
541 "This action can only performed after substitute() has been called."));
542
543 return dependent_variables_output;
544 }
545
546
547
548 template <typename ReturnType>
549 ReturnType
551 const Expression &func,
552 const std::vector<ReturnType> &cached_evaluation) const
553 {
554 // TODO[JPP]: Find a way to fix this bug that crops up in serialization
555 // cases, e.g. symengine/batch_optimizer_05. Even though the entry is
556 // in the map, it can only be found by an exhaustive search and string
557 // comparison. Why? Because the leading zero coefficient may seemingly
558 // be dropped (or added) at any time.
559 //
560 // Just this should theoretically work:
561 const typename map_dependent_expression_to_vector_entry_t::const_iterator
562 it = map_dep_expr_vec_entry.find(func);
563
564 // But instead we are forced to live with this abomination, and its
565 // knock-on effects:
566 if (has_been_serialized && it == map_dep_expr_vec_entry.end())
567 {
568 // Some SymEngine operations might return results with a zero leading
569 // coefficient. Upon serialization, this might be dropped, meaning
570 // that when we reload the expressions they now look somewhat
571 // different to as before. If all data that the user uses is
572 // guaranteed to either have been serialized or never serialized, then
573 // there would be no problem. However, users might rebuild their
574 // dependent expression and just reload the optimizer. This is
575 // completely legitimate. But in this scenario we might be out of sync
576 // with the expressions. This is not great. So we take the nuclear
577 // approach, and run everything through a serialization operation to
578 // see if we can homogenize all of the expressions such that they look
579 // the same in string form.
580 auto serialize_and_deserialize_expression =
581 [](const Expression &old_expr) {
582 std::ostringstream oss;
583 {
584 boost::archive::text_oarchive oa(oss,
585 boost::archive::no_header);
586 oa << old_expr;
587 }
588
589 Expression new_expr;
590 {
591 std::istringstream iss(oss.str());
592 boost::archive::text_iarchive ia(iss,
593 boost::archive::no_header);
594
595 ia >> new_expr;
596 }
597
598 return new_expr;
599 };
600
601 const Expression new_func =
602 serialize_and_deserialize_expression(func);
603
604 // Find this in the map, while also making sure to compactify all map
605 // entries. If we find the entry that we're looking for, then we
606 // (re-)add the input expression into the map, and do the proper
607 // search again. We should only need to do this once per invalid
608 // entry, as the corrected entry is then cached in the map.
609 for (const auto &e : map_dep_expr_vec_entry)
610 {
611 const Expression new_map_expr =
612 serialize_and_deserialize_expression(e.first);
613
614 // Add a new map entry and re-search. This is guaranteed to
615 // return a valid entry. Note that we must do a string comparison,
616 // because the data structures that form the expressions might
617 // still be different.
618 if (new_func.get_value().__str__() ==
619 new_map_expr.get_value().__str__())
620 {
621 map_dep_expr_vec_entry[func] = e.second;
622 return extract(func, cached_evaluation);
623 }
624 }
625
627 false,
629 "Still cannot find map entry, and there's no hope to recover from this situation."));
630 }
631
632 Assert(it != map_dep_expr_vec_entry.end(),
633 ExcMessage("Function has not been registered."));
634 Assert(it->second < n_dependent_variables(), ExcInternalError());
635
636 return cached_evaluation[it->second];
637 }
638
639
640
641 template <typename ReturnType>
642 ReturnType
644 {
645 Assert(
646 values_substituted() == true,
648 "The optimizer is not configured to perform evaluation. "
649 "This action can only performed after substitute() has been called."));
650
651 return extract(func, dependent_variables_output);
652 }
653
654
655
656 template <typename ReturnType>
657 std::vector<ReturnType>
659 const std::vector<Expression> &funcs,
660 const std::vector<ReturnType> &cached_evaluation) const
661 {
662 std::vector<ReturnType> out;
663 out.reserve(funcs.size());
664
665 for (const auto &func : funcs)
666 out.emplace_back(extract(func, cached_evaluation));
667
668 return out;
669 }
670
671
672
673 template <typename ReturnType>
674 std::vector<ReturnType>
676 const std::vector<Expression> &funcs) const
677 {
678 Assert(
679 values_substituted() == true,
681 "The optimizer is not configured to perform evaluation. "
682 "This action can only performed after substitute() has been called."));
683 return extract(funcs, dependent_variables_output);
684 }
685
686
687
688 template <typename ReturnType>
689 bool
691 const SD::Expression &func) const
692 {
693 return is_valid_nonunique_dependent_variable(func.get_RCP());
694 }
695
696
697
698 template <typename ReturnType>
699 bool
701 const SymEngine::RCP<const SymEngine::Basic> &func) const
702 {
703 // SymEngine's internal constants are the valid
704 // reusable return types for various derivative operations
705 // See
706 // https://github.com/symengine/symengine/blob/master/symengine/constants.h
707 if (SymEngine::is_a<SymEngine::Constant>(*func))
708 return true;
709 if (&*func == &*SymEngine::zero)
710 return true;
711 if (&*func == &*SymEngine::one)
712 return true;
713 if (&*func == &*SymEngine::minus_one)
714 return true;
715 if (&*func == &*SymEngine::I)
716 return true;
717 if (&*func == &*SymEngine::Inf)
718 return true;
719 if (&*func == &*SymEngine::NegInf)
720 return true;
721 if (&*func == &*SymEngine::ComplexInf)
722 return true;
723 if (&*func == &*SymEngine::Nan)
724 return true;
725
726 return false;
727 }
728
729
730
731 template <typename ReturnType>
732 void
734 const SD::Expression &func)
735 {
736 Assert(
737 dependent_variables_output.empty(),
739 "Cannot register function as the optimizer has already been finalized."));
740 dependent_variables_output.reserve(n_dependent_variables() + 1);
741 const bool entry_registered =
742 (map_dep_expr_vec_entry.find(func) != map_dep_expr_vec_entry.end());
743 if constexpr (running_in_debug_mode())
744 {
745 if (entry_registered == true &&
746 is_valid_nonunique_dependent_variable(func) == false)
747 Assert(entry_registered,
748 ExcMessage("Function has already been registered."));
749 }
750 if (entry_registered == false)
751 {
752 dependent_variables_functions.push_back(func);
753 map_dep_expr_vec_entry[func] =
754 dependent_variables_functions.size() - 1;
755 }
756 }
757
758
759
760 template <typename ReturnType>
761 void
763 const SD::types::symbol_vector &funcs)
764 {
765 Assert(
766 dependent_variables_output.empty(),
768 "Cannot register function as the optimizer has already been finalized."));
769 const std::size_t n_dependents_old = n_dependent_variables();
770 dependent_variables_output.reserve(n_dependents_old + funcs.size());
771 dependent_variables_functions.reserve(n_dependents_old + funcs.size());
772
773 for (const auto &func : funcs)
774 {
775 const bool entry_registered =
776 (map_dep_expr_vec_entry.find(func) != map_dep_expr_vec_entry.end());
777 if constexpr (running_in_debug_mode())
778 {
779 if (entry_registered == true &&
780 is_valid_nonunique_dependent_variable(func) == false)
781 Assert(entry_registered,
782 ExcMessage("Function has already been registered."));
783 }
784 if (entry_registered == false)
785 {
786 dependent_variables_functions.push_back(func);
787 map_dep_expr_vec_entry[func] =
788 dependent_variables_functions.size() - 1;
789 }
790 }
791 }
792
793
794
795 template <typename ReturnType>
796 void
798 std::unique_ptr<SymEngine::Visitor> &optimizer)
799 {
800 Assert(!optimizer, ExcMessage("Optimizer has already been created."));
801
802 if (optimization_method() == OptimizerType::dictionary ||
803 optimization_method() == OptimizerType::dictionary)
804 {
805 using Optimizer_t =
807 optimizer.reset(new Optimizer_t());
808 }
809 else if (optimization_method() == OptimizerType::lambda)
810 {
811 using Optimizer_t =
813 optimizer.reset(new Optimizer_t());
814 }
815 else if (optimization_method() == OptimizerType::llvm)
816 {
817# ifdef HAVE_SYMENGINE_LLVM
818 if (internal::LLVMOptimizer<ReturnType>::supported_by_LLVM)
819 {
820 using Optimizer_t =
821 typename internal::LLVMOptimizer<ReturnType>::OptimizerType;
822 optimizer.reset(new Optimizer_t());
823 }
824 else
825 {
827 }
828# else
830# endif
831 }
832 else
833 {
834 AssertThrow(false, ExcMessage("Unknown optimizer selected."));
835 }
836 }
837
838 } // namespace SD
839} // namespace Differentiation
840
841
842/* --- Explicit instantiations --- */
843# include "differentiation/sd/symengine_optimizer.inst"
844
845
846
847#endif // DEAL_II_WITH_SYMENGINE
types::substitution_map independent_variables_symbols
void substitute(const types::substitution_map &substitution_map) const
void register_scalar_function(const SD::Expression &function)
const types::symbol_vector & get_dependent_functions() const
void create_optimizer(std::unique_ptr< SymEngine::Visitor > &optimizer)
enum OptimizerType optimization_method() const
void copy_from(const BatchOptimizer &other)
void set_optimization_method(const enum OptimizerType &optimization_method, const enum OptimizationFlags &optimization_flags=OptimizationFlags::optimize_all)
enum OptimizationFlags optimization_flags() const
void register_functions(const types::symbol_vector &functions)
void register_symbols(const types::substitution_map &substitution_map)
const std::vector< ReturnType > & evaluate() const
void register_function(const Expression &function)
ReturnType extract(const Expression &func, const std::vector< ReturnType > &cached_evaluation) const
bool is_valid_nonunique_dependent_variable(const SD::Expression &function) const
void register_vector_functions(const types::symbol_vector &functions)
types::symbol_vector get_independent_symbols() const
map_dependent_expression_to_vector_entry_t map_dep_expr_vec_entry
const SymEngine::RCP< const SymEngine::Basic > & get_RCP() const
const SymEngine::Basic & get_value() const
#define DEAL_II_NAMESPACE_OPEN
Definition config.h:38
constexpr bool running_in_debug_mode()
Definition config.h:76
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:39
static ::ExceptionBase & ExcNotImplemented()
#define Assert(cond, exc)
static ::ExceptionBase & ExcSymEngineLLVMReturnTypeNotSupported()
static ::ExceptionBase & ExcSymEngineLLVMNotAvailable()
static ::ExceptionBase & ExcInternalError()
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
static ::ExceptionBase & ExcNotInitialized()
static ::ExceptionBase & ExcMessage(std::string arg1)
#define AssertThrow(cond, exc)
SD::types::symbol_vector extract_symbols(const SD::types::substitution_map &substitution_values)
SD::types::substitution_map convert_basic_map_to_expression_map(const SymEngine::map_basic_basic &substitution_map)
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)
bool use_symbolic_CSE(const enum OptimizationFlags &flags)
std::vector< SD::Expression > symbol_vector
std::map< SD::Expression, SD::Expression, internal::ExpressionKeyLess > substitution_map
Expression substitute(const Expression &expression, const types::substitution_map &substitution_map)
types::substitution_map make_substitution_map(const Expression &symbol, const Expression &value)
constexpr bool values_are_equal(const Number1 &value_1, const Number2 &value_2)
Definition numbers.h:858