Reference documentation for deal.II version 9.4.0
\(\newcommand{\dealvcentcolon}{\mathrel{\mathop{:}}}\) \(\newcommand{\dealcoloneq}{\dealvcentcolon\mathrel{\mkern-1.2mu}=}\) \(\newcommand{\jump}[1]{\left[\!\left[ #1 \right]\!\right]}\) \(\newcommand{\average}[1]{\left\{\!\left\{ #1 \right\}\!\right\}}\)
symengine_optimizer.cc
Go to the documentation of this file.
1// ---------------------------------------------------------------------
2//
3// Copyright (C) 2020 - 2021 by the deal.II authors
4//
5// This file is part of the deal.II library.
6//
7// The deal.II library is free software; you can use it, redistribute
8// it, and/or modify it under the terms of the GNU Lesser General
9// Public License as published by the Free Software Foundation; either
10// version 2.1 of the License, or (at your option) any later version.
11// The full text of the license can be found in the file LICENSE at
12// the top level of the deal.II distribution.
13//
14// ---------------------------------------------------------------------
15
16#include <deal.II/base/config.h>
17
18#ifdef DEAL_II_WITH_SYMENGINE
19
22
23# include <boost/archive/text_iarchive.hpp>
24# include <boost/archive/text_oarchive.hpp>
25
26# include <utility>
27
29
30
31namespace Differentiation
32{
33 namespace SD
34 {
35 template <typename ReturnType>
37 : method(OptimizerType::dictionary)
39 , ready_for_value_extraction(false)
40 , has_been_serialized(false)
41 {}
42
43
44
45 template <typename ReturnType>
47 const enum OptimizerType & optimization_method,
48 const enum OptimizationFlags &optimization_flags)
50 {
52 }
53
54
55
56 template <typename ReturnType>
58 const BatchOptimizer<ReturnType> &other)
59 : method(other.method)
60 , flags(other.flags)
61 , independent_variables_symbols(other.independent_variables_symbols)
62 , dependent_variables_functions(other.dependent_variables_functions)
63 , dependent_variables_output(0)
64 , map_dep_expr_vec_entry(other.map_dep_expr_vec_entry)
65 , ready_for_value_extraction(false)
66 , has_been_serialized(false)
67 {}
68
69
70
71 template <typename ReturnType>
72 void
74 const BatchOptimizer<ReturnType> &other)
75 {
76 method = other.method;
77 flags = other.flags;
78 independent_variables_symbols = other.independent_variables_symbols;
79 dependent_variables_functions = other.dependent_variables_functions;
80 dependent_variables_output.clear();
81 map_dep_expr_vec_entry = other.map_dep_expr_vec_entry;
82 ready_for_value_extraction = false;
83 has_been_serialized = false;
84 }
85
86
87
88 template <typename ReturnType>
89 void
91 const enum OptimizerType & optimization_method,
92 const enum OptimizationFlags &optimization_flags)
93 {
94 Assert(
95 optimized() == false,
97 "Cannot call set_optimization_method() once the optimizer is finalized."));
98
99# ifndef HAVE_SYMENGINE_LLVM
100 if (optimization_method == OptimizerType::llvm)
101 {
103 }
104# endif
105 method = optimization_method;
106 flags = optimization_flags;
107 }
108
109
110
111 template <typename ReturnType>
112 enum OptimizerType
114 {
115 return method;
116 }
117
118
119
120 template <typename ReturnType>
123 {
124 return flags;
125 }
126
127
128
129 template <typename ReturnType>
130 bool
132 {
133 return internal::use_symbolic_CSE(flags);
134 }
135
136
137
138 template <typename ReturnType>
139 bool
141 {
142 if (dependent_variables_output.size() > 0)
143 {
144 Assert(dependent_variables_output.size() ==
145 dependent_variables_functions.size(),
147 return true;
148 }
149
150 return false;
151 }
152
153
154
155 template <typename ReturnType>
156 bool
158 {
159 return ready_for_value_extraction;
160 }
161
162
163
164 template <typename ReturnType>
165 void
168 {
169 Assert(optimized() == false,
171 "Cannot register symbols once the optimizer is finalized."));
172
173# ifdef DEBUG
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# endif
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
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 {
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
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
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# ifdef DEBUG
408 const SD::types::symbol_vector symbol_sub_vec =
410 const SD::types::symbol_vector symbol_vec =
411 Utilities::extract_symbols(independent_variables_symbols);
412 Assert(symbol_sub_vec.size() == symbol_vec.size(),
413 ExcDimensionMismatch(symbol_sub_vec.size(), symbol_vec.size()));
414 for (unsigned int i = 0; i < symbol_sub_vec.size(); ++i)
415 {
416 Assert(numbers::values_are_equal(symbol_sub_vec[i], symbol_vec[i]),
418 "The input substitution map is either incomplete, or does "
419 "not match that used in the register_symbols() call."));
420 }
421# endif
422
423 // Extract the values from the substitution map, and use the other
424 // function
425 const std::vector<ReturnType> values =
426 Utilities::extract_values<ReturnType>(substitution_map);
428 }
429
430
431
432 template <typename ReturnType>
433 void
435 const SymEngine::map_basic_basic &substitution_map) const
436 {
439 }
440
441
442
443 template <typename ReturnType>
444 void
446 const SD::types::symbol_vector &symbols,
447 const std::vector<ReturnType> & values) const
448 {
449 // Zip the two vectors and use the other function call
450 // This ensures the ordering of the input vectors matches that of the
451 // stored map.
453 }
454
455
456
457 template <typename ReturnType>
458 void
460 const SymEngine::vec_basic & symbols,
461 const std::vector<ReturnType> &values) const
462 {
464 symbols),
465 values);
466 }
467
468
469
470 template <typename ReturnType>
471 void
473 const std::vector<ReturnType> &substitution_values) const
474 {
475 Assert(
476 optimized() == true,
478 "The optimizer is not configured to perform substitution. "
479 "This action can only performed after optimize() has been called."));
480 Assert(optimizer, ExcNotInitialized());
481 Assert(substitution_values.size() == independent_variables_symbols.size(),
482 ExcDimensionMismatch(substitution_values.size(),
483 independent_variables_symbols.size()));
484
486 *opt = dynamic_cast<typename internal::DictionaryOptimizer<
487 ReturnType>::OptimizerType *>(optimizer.get()))
488 {
489 Assert(optimization_method() == OptimizerType::dictionary,
491 internal::OptimizerHelper<ReturnType,
493 substitute(opt, dependent_variables_output, substitution_values);
494 }
496 *opt = dynamic_cast<typename internal::LambdaOptimizer<
497 ReturnType>::OptimizerType *>(optimizer.get()))
498 {
499 Assert(optimization_method() == OptimizerType::lambda,
501 internal::OptimizerHelper<ReturnType,
503 substitute(opt, dependent_variables_output, substitution_values);
504 }
505# ifdef HAVE_SYMENGINE_LLVM
507 *opt = dynamic_cast<typename internal::LLVMOptimizer<
508 ReturnType>::OptimizerType *>(optimizer.get()))
509 {
510 Assert(optimization_method() == OptimizerType::llvm,
512 internal::OptimizerHelper<ReturnType,
513 internal::LLVMOptimizer<ReturnType>>::
514 substitute(opt, dependent_variables_output, substitution_values);
515 }
516# endif
517 else
518 {
520 }
521
522 ready_for_value_extraction = true;
523 }
524
525
526
527 template <typename ReturnType>
528 const std::vector<ReturnType> &
530 {
531 Assert(
532 values_substituted() == true,
534 "The optimizer is not configured to perform evaluation. "
535 "This action can only performed after substitute() has been called."));
536
537 return dependent_variables_output;
538 }
539
540
541
542 template <typename ReturnType>
543 ReturnType
545 const Expression & func,
546 const std::vector<ReturnType> &cached_evaluation) const
547 {
548 // TODO[JPP]: Find a way to fix this bug that crops up in serialization
549 // cases, e.g. symengine/batch_optimizer_05. Even though the entry is
550 // in the map, it can only be found by an exhaustive search and string
551 // comparison. Why? Because the leading zero coefficient may seemingly
552 // be dropped (or added) at any time.
553 //
554 // Just this should theoretically work:
555 const typename map_dependent_expression_to_vector_entry_t::const_iterator
556 it = map_dep_expr_vec_entry.find(func);
557
558 // But instead we are forced to live with this abomination, and its
559 // knock-on effects:
560 if (has_been_serialized && it == map_dep_expr_vec_entry.end())
561 {
562 // Some SymEngine operations might return results with a zero leading
563 // coefficient. Upon serialization, this might be dropped, meaning
564 // that when we reload the expressions they now look somewhat
565 // different to as before. If all data that the user uses is
566 // guaranteed to either have been serialized or never serialized, then
567 // there would be no problem. However, users might rebuild their
568 // dependent expression and just reload the optimizer. This is
569 // completely legitimate. But in this scenario we might be out of sync
570 // with the expressions. This is not great. So we take the nuclear
571 // approach, and run everything through a serialization operation to
572 // see if we can homogenize all of the expressions such that they look
573 // the same in string form.
574 auto serialize_and_deserialize_expression =
575 [](const Expression &old_expr) {
576 std::ostringstream oss;
577 {
578 boost::archive::text_oarchive oa(oss,
579 boost::archive::no_header);
580 oa << old_expr;
581 }
582
583 Expression new_expr;
584 {
585 std::istringstream iss(oss.str());
586 boost::archive::text_iarchive ia(iss,
587 boost::archive::no_header);
588
589 ia >> new_expr;
590 }
591
592 return new_expr;
593 };
594
595 const Expression new_func =
596 serialize_and_deserialize_expression(func);
597
598 // Find this in the map, while also making sure to compactify all map
599 // entries. If we find the entry that we're looking for, then we
600 // (re-)add the input expression into the map, and do the proper
601 // search again. We should only need to do this once per invalid
602 // entry, as the corrected entry is then cached in the map.
603 for (const auto &e : map_dep_expr_vec_entry)
604 {
605 const Expression new_map_expr =
606 serialize_and_deserialize_expression(e.first);
607
608 // Add a new map entry and re-search. This is guaranteed to
609 // return a valid entry. Note that we must do a string comparison,
610 // because the data structures that form the expressions might
611 // still be different.
612 if (new_func.get_value().__str__() ==
613 new_map_expr.get_value().__str__())
614 {
615 map_dep_expr_vec_entry[func] = e.second;
616 return extract(func, cached_evaluation);
617 }
618 }
619
621 false,
623 "Still cannot find map entry, and there's no hope to recover from this situation."));
624 }
625
626 Assert(it != map_dep_expr_vec_entry.end(),
627 ExcMessage("Function has not been registered."));
628 Assert(it->second < n_dependent_variables(), ExcInternalError());
629
630 return cached_evaluation[it->second];
631 }
632
633
634
635 template <typename ReturnType>
636 ReturnType
638 {
639 Assert(
640 values_substituted() == true,
642 "The optimizer is not configured to perform evaluation. "
643 "This action can only performed after substitute() has been called."));
644
645 return extract(func, dependent_variables_output);
646 }
647
648
649
650 template <typename ReturnType>
651 std::vector<ReturnType>
653 const std::vector<Expression> &funcs,
654 const std::vector<ReturnType> &cached_evaluation) const
655 {
656 std::vector<ReturnType> out;
657 out.reserve(funcs.size());
658
659 for (const auto &func : funcs)
660 out.emplace_back(extract(func, cached_evaluation));
661
662 return out;
663 }
664
665
666
667 template <typename ReturnType>
668 std::vector<ReturnType>
670 const std::vector<Expression> &funcs) const
671 {
672 Assert(
673 values_substituted() == true,
675 "The optimizer is not configured to perform evaluation. "
676 "This action can only performed after substitute() has been called."));
677 return extract(funcs, dependent_variables_output);
678 }
679
680
681
682 template <typename ReturnType>
683 bool
685 const SD::Expression &func) const
686 {
687 return is_valid_nonunique_dependent_variable(func.get_RCP());
688 }
689
690
691
692 template <typename ReturnType>
693 bool
695 const SymEngine::RCP<const SymEngine::Basic> &func) const
696 {
697 // SymEngine's internal constants are the valid
698 // reusable return types for various derivative operations
699 // See
700 // https://github.com/symengine/symengine/blob/master/symengine/constants.h
701 if (SymEngine::is_a<SymEngine::Constant>(*func))
702 return true;
703 if (&*func == &*SymEngine::zero)
704 return true;
705 if (&*func == &*SymEngine::one)
706 return true;
707 if (&*func == &*SymEngine::minus_one)
708 return true;
709 if (&*func == &*SymEngine::I)
710 return true;
711 if (&*func == &*SymEngine::Inf)
712 return true;
713 if (&*func == &*SymEngine::NegInf)
714 return true;
715 if (&*func == &*SymEngine::ComplexInf)
716 return true;
717 if (&*func == &*SymEngine::Nan)
718 return true;
719
720 return false;
721 }
722
723
724
725 template <typename ReturnType>
726 void
728 const SD::Expression &func)
729 {
730 Assert(
731 dependent_variables_output.size() == 0,
733 "Cannot register function as the optimizer has already been finalized."));
734 dependent_variables_output.reserve(n_dependent_variables() + 1);
735 const bool entry_registered =
736 (map_dep_expr_vec_entry.find(func) != map_dep_expr_vec_entry.end());
737# ifdef DEBUG
738 if (entry_registered == true &&
739 is_valid_nonunique_dependent_variable(func) == false)
740 Assert(entry_registered,
741 ExcMessage("Function has already been registered."));
742# endif
743 if (entry_registered == false)
744 {
745 dependent_variables_functions.push_back(func);
746 map_dep_expr_vec_entry[func] =
747 dependent_variables_functions.size() - 1;
748 }
749 }
750
751
752
753 template <typename ReturnType>
754 void
756 const SD::types::symbol_vector &funcs)
757 {
758 Assert(
759 dependent_variables_output.size() == 0,
761 "Cannot register function as the optimizer has already been finalized."));
762 const std::size_t n_dependents_old = n_dependent_variables();
763 dependent_variables_output.reserve(n_dependents_old + funcs.size());
764 dependent_variables_functions.reserve(n_dependents_old + funcs.size());
765
766 for (const auto &func : funcs)
767 {
768 const bool entry_registered =
769 (map_dep_expr_vec_entry.find(func) != map_dep_expr_vec_entry.end());
770# ifdef DEBUG
771 if (entry_registered == true &&
772 is_valid_nonunique_dependent_variable(func) == false)
773 Assert(entry_registered,
774 ExcMessage("Function has already been registered."));
775# endif
776 if (entry_registered == false)
777 {
778 dependent_variables_functions.push_back(func);
779 map_dep_expr_vec_entry[func] =
780 dependent_variables_functions.size() - 1;
781 }
782 }
783 }
784
785
786
787 template <typename ReturnType>
788 void
790 std::unique_ptr<SymEngine::Visitor> &optimizer)
791 {
792 Assert(!optimizer, ExcMessage("Optimizer has already been created."));
793
794 if (optimization_method() == OptimizerType::dictionary ||
795 optimization_method() == OptimizerType::dictionary)
796 {
797 using Optimizer_t =
799 optimizer.reset(new Optimizer_t());
800 }
801 else if (optimization_method() == OptimizerType::lambda)
802 {
803 using Optimizer_t =
805 optimizer.reset(new Optimizer_t());
806 }
807 else if (optimization_method() == OptimizerType::llvm)
808 {
809# ifdef HAVE_SYMENGINE_LLVM
810 if (internal::LLVMOptimizer<ReturnType>::supported_by_LLVM)
811 {
812 using Optimizer_t =
814 optimizer.reset(new Optimizer_t());
815 }
816 else
817 {
819 }
820# else
822# endif
823 }
824 else
825 {
826 AssertThrow(false, ExcMessage("Unknown optimizer selected."));
827 }
828 }
829
830 } // namespace SD
831} // namespace Differentiation
832
833
834/* --- Explicit instantiations --- */
835# include "symengine_optimizer.inst"
836
837
839
840#endif // DEAL_II_WITH_SYMENGINE
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:442
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:443
types::substitution_map independent_variables_symbols
types::symbol_vector dependent_variables_functions
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
const SymEngine::RCP< const SymEngine::Basic > & get_RCP() const
void copy_from(const BatchOptimizer &other)
static ::ExceptionBase & ExcNotImplemented()
void set_optimization_method(const enum OptimizerType &optimization_method, const enum OptimizationFlags &optimization_flags=OptimizationFlags::optimize_all)
bool use_symbolic_CSE(const enum OptimizationFlags &flags)
enum OptimizationFlags optimization_flags() const
#define Assert(cond, exc)
Definition: exceptions.h:1473
void register_functions(const types::symbol_vector &functions)
static ::ExceptionBase & ExcSymEngineLLVMReturnTypeNotSupported()
void register_symbols(const types::substitution_map &substitution_map)
const std::vector< ReturnType > & evaluate() const
static ::ExceptionBase & ExcSymEngineLLVMNotAvailable()
static ::ExceptionBase & ExcInternalError()
void register_function(const Expression &function)
ReturnType extract(const Expression &func, const std::vector< ReturnType > &cached_evaluation) const
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
static ::ExceptionBase & ExcNotInitialized()
bool is_valid_nonunique_dependent_variable(const SD::Expression &function) const
void register_vector_functions(const types::symbol_vector &functions)
const SymEngine::Basic & get_value() const
void register_functions(BatchOptimizer< NumberType > &optimizer, const T &function)
static ::ExceptionBase & ExcMessage(std::string arg1)
types::symbol_vector get_independent_symbols() const
map_dependent_expression_to_vector_entry_t map_dep_expr_vec_entry
#define AssertThrow(cond, exc)
Definition: exceptions.h:1583
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)
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)
static const types::blas_int zero
static const types::blas_int one
SymmetricTensor< 2, dim, Number > e(const Tensor< 2, dim, Number > &F)
constexpr ReturnType< rank, T >::value_type & extract(T &t, const ArrayType &indices)
int(&) functions(const void *v1, const void *v2)
constexpr bool values_are_equal(const Number1 &value_1, const Number2 &value_2)
Definition: numbers.h:927