Reference documentation for deal.II version GIT relicensing-462-g03c82a14bd 2024-04-21 16:20:02+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# ifdef DEBUG
173 // Ensure that all of the keys in the map are actually symbolic
174 // in nature
175 for (const auto &entry : substitution_map)
176 {
177 const SD::Expression &symbol = entry.first;
178 Assert(SymEngine::is_a<SymEngine::Symbol>(*(symbol.get_RCP())),
179 ExcMessage("Key entry in map is not a symbol."));
180 }
181# endif
182 // Merge the two maps, in the process ensuring that there is no
183 // duplication of symbols
184 independent_variables_symbols.insert(substitution_map.begin(),
185 substitution_map.end());
186 }
187
188
189
190 template <typename ReturnType>
191 void
193 const SymEngine::map_basic_basic &substitution_map)
194 {
195 register_symbols(
197 }
198
199
200
201 template <typename ReturnType>
202 void
204 const SD::types::symbol_vector &symbols)
205 {
206 Assert(optimized() == false,
208 "Cannot register symbols once the optimizer is finalized."));
209
210 for (const auto &symbol : symbols)
211 {
212 Assert(independent_variables_symbols.find(symbol) ==
213 independent_variables_symbols.end(),
214 ExcMessage("Symbol is already in the map."));
215 independent_variables_symbols.insert(
216 std::make_pair(symbol, SD::Expression(0.0)));
217 }
218 }
219
220
221
222 template <typename ReturnType>
223 void
225 const SymEngine::vec_basic &symbols)
226 {
227 register_symbols(
229 }
230
231
232
233 template <typename ReturnType>
236 {
237 return Utilities::extract_symbols(independent_variables_symbols);
238 }
239
240
241
242 template <typename ReturnType>
243 std::size_t
245 {
246 return independent_variables_symbols.size();
247 }
248
249
250
251 template <typename ReturnType>
252 void
254 {
255 Assert(optimized() == false,
257 "Cannot register functions once the optimizer is finalized."));
258
259 register_scalar_function(function);
260 }
261
262
263
264 template <typename ReturnType>
265 void
267 const SD::types::symbol_vector &functions)
268 {
269 Assert(optimized() == false,
271 "Cannot register functions once the optimizer is finalized."));
272
273 register_vector_functions(functions);
274 }
275
276
277
278 template <typename ReturnType>
279 void
281 const SymEngine::vec_basic &functions)
282 {
283 register_functions(
285 }
286
287
288
289 template <typename ReturnType>
292 {
293 return dependent_variables_functions;
294 }
295
296
297
298 template <typename ReturnType>
299 std::size_t
301 {
302 if (has_been_serialized == false)
303 {
304 // If we've had to augment our map after serialization, then
305 // this check, unfortunately, cannot be performed.
306 Assert(map_dep_expr_vec_entry.size() ==
307 dependent_variables_functions.size(),
309 }
310 return dependent_variables_functions.size();
311 }
312
313
314
315 template <typename ReturnType>
316 void
318 {
319 Assert(optimized() == false,
320 ExcMessage("Cannot call optimize() more than once."));
321
322 // Create and configure the optimizer
323 create_optimizer(optimizer);
324 Assert(optimizer, ExcNotInitialized());
325
326 const SD::types::symbol_vector symbol_vec =
327 Utilities::extract_symbols(independent_variables_symbols);
329 *opt = dynamic_cast<typename internal::DictionaryOptimizer<
330 ReturnType>::OptimizerType *>(optimizer.get()))
331 {
332 Assert(optimization_method() == OptimizerType::dictionary,
334 internal::OptimizerHelper<ReturnType,
336 initialize(opt,
338 symbol_vec),
340 dependent_variables_functions),
341 optimization_flags());
342 }
344 *opt = dynamic_cast<typename internal::LambdaOptimizer<
345 ReturnType>::OptimizerType *>(optimizer.get()))
346 {
347 Assert(optimization_method() == OptimizerType::lambda,
349 internal::OptimizerHelper<ReturnType,
351 initialize(opt,
353 symbol_vec),
355 dependent_variables_functions),
356 optimization_flags());
357 }
358# ifdef HAVE_SYMENGINE_LLVM
359 else if (typename internal::LLVMOptimizer<ReturnType>::OptimizerType
360 *opt = dynamic_cast<typename internal::LLVMOptimizer<
361 ReturnType>::OptimizerType *>(optimizer.get()))
362 {
363 Assert(optimization_method() == OptimizerType::llvm,
365 internal::OptimizerHelper<ReturnType,
366 internal::LLVMOptimizer<ReturnType>>::
367 initialize(opt,
369 symbol_vec),
371 dependent_variables_functions),
372 optimization_flags());
373 }
374# endif
375 else
376 {
377 AssertThrow(false, ExcMessage("Unknown optimizer type."));
378 }
379
380 // The size of the outputs is now fixed, as is the number and
381 // order of the symbols to be substituted.
382 // Note: When no optimisation is actually used (i.e. optimization_method()
383 // == off and use_symbolic_CSE() == false), we could conceptually go
384 // without this data structure. However, since the user expects to perform
385 // substitution of all dependent variables in one go, we still require it
386 // for intermediate storage of results.
387 dependent_variables_output.resize(n_dependent_variables());
388 }
389
390
391
392 template <typename ReturnType>
393 void
395 const SD::types::substitution_map &substitution_map) const
396 {
397 Assert(
398 optimized() == true,
400 "The optimizer is not configured to perform substitution. "
401 "This action can only performed after optimize() has been called."));
402 Assert(optimizer, ExcNotInitialized());
403
404 // Check that the registered symbol map and the input map are compatible
405 // with one another
406# ifdef DEBUG
407 const SD::types::symbol_vector symbol_sub_vec =
408 Utilities::extract_symbols(substitution_map);
409 const SD::types::symbol_vector symbol_vec =
410 Utilities::extract_symbols(independent_variables_symbols);
411 Assert(symbol_sub_vec.size() == symbol_vec.size(),
412 ExcDimensionMismatch(symbol_sub_vec.size(), symbol_vec.size()));
413 for (unsigned int i = 0; i < symbol_sub_vec.size(); ++i)
414 {
415 Assert(numbers::values_are_equal(symbol_sub_vec[i], symbol_vec[i]),
417 "The input substitution map is either incomplete, or does "
418 "not match that used in the register_symbols() call."));
419 }
420# endif
421
422 // Extract the values from the substitution map, and use the other
423 // function
424 const std::vector<ReturnType> values =
425 Utilities::extract_values<ReturnType>(substitution_map);
426 substitute(values);
427 }
428
429
430
431 template <typename ReturnType>
432 void
434 const SymEngine::map_basic_basic &substitution_map) const
435 {
438 }
439
440
441
442 template <typename ReturnType>
443 void
445 const SD::types::symbol_vector &symbols,
446 const std::vector<ReturnType> &values) const
447 {
448 // Zip the two vectors and use the other function call
449 // This ensures the ordering of the input vectors matches that of the
450 // stored map.
451 substitute(make_substitution_map(symbols, values));
452 }
453
454
455
456 template <typename ReturnType>
457 void
459 const SymEngine::vec_basic &symbols,
460 const std::vector<ReturnType> &values) const
461 {
463 symbols),
464 values);
465 }
466
467
468
469 template <typename ReturnType>
470 void
472 const std::vector<ReturnType> &substitution_values) const
473 {
474 Assert(
475 optimized() == true,
477 "The optimizer is not configured to perform substitution. "
478 "This action can only performed after optimize() has been called."));
479 Assert(optimizer, ExcNotInitialized());
480 Assert(substitution_values.size() == independent_variables_symbols.size(),
481 ExcDimensionMismatch(substitution_values.size(),
482 independent_variables_symbols.size()));
483
485 *opt = dynamic_cast<typename internal::DictionaryOptimizer<
486 ReturnType>::OptimizerType *>(optimizer.get()))
487 {
488 Assert(optimization_method() == OptimizerType::dictionary,
490 internal::OptimizerHelper<ReturnType,
492 substitute(opt, dependent_variables_output, substitution_values);
493 }
495 *opt = dynamic_cast<typename internal::LambdaOptimizer<
496 ReturnType>::OptimizerType *>(optimizer.get()))
497 {
498 Assert(optimization_method() == OptimizerType::lambda,
500 internal::OptimizerHelper<ReturnType,
502 substitute(opt, dependent_variables_output, substitution_values);
503 }
504# ifdef HAVE_SYMENGINE_LLVM
505 else if (typename internal::LLVMOptimizer<ReturnType>::OptimizerType
506 *opt = dynamic_cast<typename internal::LLVMOptimizer<
507 ReturnType>::OptimizerType *>(optimizer.get()))
508 {
509 Assert(optimization_method() == OptimizerType::llvm,
511 internal::OptimizerHelper<ReturnType,
512 internal::LLVMOptimizer<ReturnType>>::
513 substitute(opt, dependent_variables_output, substitution_values);
514 }
515# endif
516 else
517 {
519 }
520
521 ready_for_value_extraction = true;
522 }
523
524
525
526 template <typename ReturnType>
527 const std::vector<ReturnType> &
529 {
530 Assert(
531 values_substituted() == true,
533 "The optimizer is not configured to perform evaluation. "
534 "This action can only performed after substitute() has been called."));
535
536 return dependent_variables_output;
537 }
538
539
540
541 template <typename ReturnType>
542 ReturnType
544 const Expression &func,
545 const std::vector<ReturnType> &cached_evaluation) const
546 {
547 // TODO[JPP]: Find a way to fix this bug that crops up in serialization
548 // cases, e.g. symengine/batch_optimizer_05. Even though the entry is
549 // in the map, it can only be found by an exhaustive search and string
550 // comparison. Why? Because the leading zero coefficient may seemingly
551 // be dropped (or added) at any time.
552 //
553 // Just this should theoretically work:
554 const typename map_dependent_expression_to_vector_entry_t::const_iterator
555 it = map_dep_expr_vec_entry.find(func);
556
557 // But instead we are forced to live with this abomination, and its
558 // knock-on effects:
559 if (has_been_serialized && it == map_dep_expr_vec_entry.end())
560 {
561 // Some SymEngine operations might return results with a zero leading
562 // coefficient. Upon serialization, this might be dropped, meaning
563 // that when we reload the expressions they now look somewhat
564 // different to as before. If all data that the user uses is
565 // guaranteed to either have been serialized or never serialized, then
566 // there would be no problem. However, users might rebuild their
567 // dependent expression and just reload the optimizer. This is
568 // completely legitimate. But in this scenario we might be out of sync
569 // with the expressions. This is not great. So we take the nuclear
570 // approach, and run everything through a serialization operation to
571 // see if we can homogenize all of the expressions such that they look
572 // the same in string form.
573 auto serialize_and_deserialize_expression =
574 [](const Expression &old_expr) {
575 std::ostringstream oss;
576 {
577 boost::archive::text_oarchive oa(oss,
578 boost::archive::no_header);
579 oa << old_expr;
580 }
581
582 Expression new_expr;
583 {
584 std::istringstream iss(oss.str());
585 boost::archive::text_iarchive ia(iss,
586 boost::archive::no_header);
587
588 ia >> new_expr;
589 }
590
591 return new_expr;
592 };
593
594 const Expression new_func =
595 serialize_and_deserialize_expression(func);
596
597 // Find this in the map, while also making sure to compactify all map
598 // entries. If we find the entry that we're looking for, then we
599 // (re-)add the input expression into the map, and do the proper
600 // search again. We should only need to do this once per invalid
601 // entry, as the corrected entry is then cached in the map.
602 for (const auto &e : map_dep_expr_vec_entry)
603 {
604 const Expression new_map_expr =
605 serialize_and_deserialize_expression(e.first);
606
607 // Add a new map entry and re-search. This is guaranteed to
608 // return a valid entry. Note that we must do a string comparison,
609 // because the data structures that form the expressions might
610 // still be different.
611 if (new_func.get_value().__str__() ==
612 new_map_expr.get_value().__str__())
613 {
614 map_dep_expr_vec_entry[func] = e.second;
615 return extract(func, cached_evaluation);
616 }
617 }
618
620 false,
622 "Still cannot find map entry, and there's no hope to recover from this situation."));
623 }
624
625 Assert(it != map_dep_expr_vec_entry.end(),
626 ExcMessage("Function has not been registered."));
627 Assert(it->second < n_dependent_variables(), ExcInternalError());
628
629 return cached_evaluation[it->second];
630 }
631
632
633
634 template <typename ReturnType>
635 ReturnType
637 {
638 Assert(
639 values_substituted() == true,
641 "The optimizer is not configured to perform evaluation. "
642 "This action can only performed after substitute() has been called."));
643
644 return extract(func, dependent_variables_output);
645 }
646
647
648
649 template <typename ReturnType>
650 std::vector<ReturnType>
652 const std::vector<Expression> &funcs,
653 const std::vector<ReturnType> &cached_evaluation) const
654 {
655 std::vector<ReturnType> out;
656 out.reserve(funcs.size());
657
658 for (const auto &func : funcs)
659 out.emplace_back(extract(func, cached_evaluation));
660
661 return out;
662 }
663
664
665
666 template <typename ReturnType>
667 std::vector<ReturnType>
669 const std::vector<Expression> &funcs) const
670 {
671 Assert(
672 values_substituted() == true,
674 "The optimizer is not configured to perform evaluation. "
675 "This action can only performed after substitute() has been called."));
676 return extract(funcs, dependent_variables_output);
677 }
678
679
680
681 template <typename ReturnType>
682 bool
684 const SD::Expression &func) const
685 {
686 return is_valid_nonunique_dependent_variable(func.get_RCP());
687 }
688
689
690
691 template <typename ReturnType>
692 bool
694 const SymEngine::RCP<const SymEngine::Basic> &func) const
695 {
696 // SymEngine's internal constants are the valid
697 // reusable return types for various derivative operations
698 // See
699 // https://github.com/symengine/symengine/blob/master/symengine/constants.h
700 if (SymEngine::is_a<SymEngine::Constant>(*func))
701 return true;
702 if (&*func == &*SymEngine::zero)
703 return true;
704 if (&*func == &*SymEngine::one)
705 return true;
706 if (&*func == &*SymEngine::minus_one)
707 return true;
708 if (&*func == &*SymEngine::I)
709 return true;
710 if (&*func == &*SymEngine::Inf)
711 return true;
712 if (&*func == &*SymEngine::NegInf)
713 return true;
714 if (&*func == &*SymEngine::ComplexInf)
715 return true;
716 if (&*func == &*SymEngine::Nan)
717 return true;
718
719 return false;
720 }
721
722
723
724 template <typename ReturnType>
725 void
727 const SD::Expression &func)
728 {
729 Assert(
730 dependent_variables_output.empty(),
732 "Cannot register function as the optimizer has already been finalized."));
733 dependent_variables_output.reserve(n_dependent_variables() + 1);
734 const bool entry_registered =
735 (map_dep_expr_vec_entry.find(func) != map_dep_expr_vec_entry.end());
736# ifdef DEBUG
737 if (entry_registered == true &&
738 is_valid_nonunique_dependent_variable(func) == false)
739 Assert(entry_registered,
740 ExcMessage("Function has already been registered."));
741# endif
742 if (entry_registered == false)
743 {
744 dependent_variables_functions.push_back(func);
745 map_dep_expr_vec_entry[func] =
746 dependent_variables_functions.size() - 1;
747 }
748 }
749
750
751
752 template <typename ReturnType>
753 void
755 const SD::types::symbol_vector &funcs)
756 {
757 Assert(
758 dependent_variables_output.empty(),
760 "Cannot register function as the optimizer has already been finalized."));
761 const std::size_t n_dependents_old = n_dependent_variables();
762 dependent_variables_output.reserve(n_dependents_old + funcs.size());
763 dependent_variables_functions.reserve(n_dependents_old + funcs.size());
764
765 for (const auto &func : funcs)
766 {
767 const bool entry_registered =
768 (map_dep_expr_vec_entry.find(func) != map_dep_expr_vec_entry.end());
769# ifdef DEBUG
770 if (entry_registered == true &&
771 is_valid_nonunique_dependent_variable(func) == false)
772 Assert(entry_registered,
773 ExcMessage("Function has already been registered."));
774# endif
775 if (entry_registered == false)
776 {
777 dependent_variables_functions.push_back(func);
778 map_dep_expr_vec_entry[func] =
779 dependent_variables_functions.size() - 1;
780 }
781 }
782 }
783
784
785
786 template <typename ReturnType>
787 void
789 std::unique_ptr<SymEngine::Visitor> &optimizer)
790 {
791 Assert(!optimizer, ExcMessage("Optimizer has already been created."));
792
793 if (optimization_method() == OptimizerType::dictionary ||
794 optimization_method() == OptimizerType::dictionary)
795 {
796 using Optimizer_t =
798 optimizer.reset(new Optimizer_t());
799 }
800 else if (optimization_method() == OptimizerType::lambda)
801 {
802 using Optimizer_t =
804 optimizer.reset(new Optimizer_t());
805 }
806 else if (optimization_method() == OptimizerType::llvm)
807 {
808# ifdef HAVE_SYMENGINE_LLVM
809 if (internal::LLVMOptimizer<ReturnType>::supported_by_LLVM)
810 {
811 using Optimizer_t =
812 typename internal::LLVMOptimizer<ReturnType>::OptimizerType;
813 optimizer.reset(new Optimizer_t());
814 }
815 else
816 {
818 }
819# else
821# endif
822 }
823 else
824 {
825 AssertThrow(false, ExcMessage("Unknown optimizer selected."));
826 }
827 }
828
829 } // namespace SD
830} // namespace Differentiation
831
832
833/* --- Explicit instantiations --- */
834# include "symengine_optimizer.inst"
835
836
838
839#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:502
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:503
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:920