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